ang-ex4.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Neural Networks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy import optimize as opt\n",
    "from scipy import io\n",
    "from scipy.special import expit\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Feedforward"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This carries over from the last exercise. We will use a pretrained network here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Training data\n",
    "data = io.loadmat(os.path.join('data', 'ex4data1.mat'))\n",
    "X, y = data['X'], data['y'].flatten()\n",
    "y[y == 10] = 0 # from matlab indexing to python indexing\n",
    "m, n = X.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABIEAAABICAYAAABhnuWuAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJztnXmcXFWVx3+vqqu6E8JkTFgGkC2iCBIWUVlcANEBDA4MjkCYIUFFBGFgRAbIKFExIyqKyiLKaMQlATPAYBiDGFkEGWBAWQaBEAEh7EJCOr3W9uaPyve+V6+q0lu9V6/I+X0++XS6u7rer+4799z7zvndczzf92UwGAwGg8FgMBgMBoPBYHh9I9NuAgaDwWAwGAwGg8FgMBgMhvhhQSCDwWAwGAwGg8FgMBgMho0AFgQyGAwGg8FgMBgMBoPBYNgIYEEgg8FgMBgMBoPBYDAYDIaNABYEMhgMBoPBYDAYDAaDwWDYCGBBIIPBYDAYDAaDwWAwGAyGjQAWBDIYDAaDwWAwGAwGg8Fg2AgwoSCQ53mHep63wvO8P3med26rSMUJ45wMjHMyMM7JwDgnA+OcDDqNc6fxlYxzUjDOycA4JwPjnAyMczIwzimH7/vj+icpK+kJSTMk5SU9KGnX8b5fEv+Ms3E2zsbZOKfzn3E2zq8HvsbZOBvn9v8zzsbZOBvnjYHzRP556z/0mOF53n6Svuj7/iHrv58nSb7vX7CBv/EzmfadQONDwyE0EF6zv8lkMm3nXKlUlM1mJUmVSmVEzjbOY8d4x9nzmv46djB34cD3r6dxbrctS03tud/3/SmNXt+pnNtpy1JTe35djXMmk/Gx/Xah0RysVCqpHedGY7z+6yu+72/e6G/S6OfW/7wp5zSOs+01Wo+NZU9nnMcO45wMNga/0e4xljaOvXO77aIZyuVy071GGF0TuMY2klaFvn9W0j7RF3med5Kkk9b/X5MmTZrAJSeGUqmkcrms7u5uSVKxWFShUKh7XZTzpptumijPMAqFgkqlkiZPnuy+HxwcrHtdlDOvbwdKpZJKpZJ6enokVcd5eHi47nVhzplMRlOnTk2UZxiFQkHFYlGbbLKJJGl4eFgDAwN1rwtzluRsqR0ol8uqVCrK5XLu+2KxWPe6qG1MmdLQlyWCYrGoYrFYY89DQ0M1r0mTz5Ca+o0a40jT/JOazsENcua17UK5XFa5XFY+n5dU/QzFYnGDnNtpy1J1XEulkrPR9fbclHMmk9Ff//VfJ080hOHhYRUKBbeuDQ0Nqa+vb4PjnLZ1u1KpqFQqPR1+XdrXk1KppOHh4aac2+03xrNudyrntK+BUrr8c3QOlkqlUe2dO5FzJ9pz2vzz63Wc0+430mQX0vj2zp3Iud370EZYu3bt0yO/ShNSAn1U0iG+75+4/vvjJb3L9/1/bvY32WzWT6Oz2lD0t6ury09jEKhSqTTlnM1m/TQ62JHGOW2b9oGBgQ2OcyaT8VsVBCqVSjXfE+3fkDqjWRBopHFO40LWbJzb7TOkpn7jZd/3t2z0+nbPP6npHGzKOZPJ+CkNAm1wnNu9+DYJAjXlnMvl/JQGgTY4zmlbt9cHgX7v+/47Gv1NK9YT9kblclnShn0yr0Xxs9523eYwFARqyrndfmM863acnBl3xjSTySiqohvvXqOT1kCpOs5pDKiMZBudyLkT5+BY/XPUX/F1JDSbg2N9rurEce40v9HudVsa39650zi32y6aYe3atU33GmFMRMP0rKRtQ9+/UdLzE3i/2OF5nsJBL9/3N/iQnQZkMpkaBz1aZ91OdOI4e57nNp1S545z2hHl3MHjXC+5ShE6kbOkjuPseV6NDXcC5yZrSmo5d+J60miM0ygZD6MTx/n1wLlT18BOG2fjHA+MczJ4vfgNpXivIXUm54lgIsfB7pX0Zs/zdpT0nKRjJR3XElYxgY1ZpVKpe+hPK7LZrCqVisrlsjKZjAqFQuqdVXScoyqXuBDN3PKV8eKoVKNgSTabdUqEJMYZDlxj6623rvn5X/7yl4avCwNnlTZ75n4Xi0WXaYE/4wznYrHY8nFmzJAkR+0P5RRfR7p+E7/xWktJtxhN5mDqOTew5wlxxhaY+42OS44G+Xze2UsY+OewPSvl49zV1VXj69bL4GPhzNziCHMzH8XYogILo9H86+7uHve9HAnYTFdXdXu05ZbVBOCaNWskBX4lvK7wf9RHKElLpZJbTzbZZJOGR7nTgnat21FwLOM973mPJOmb3/ymJOmOO+7Q5z//eUnBPUoL57EgiTWwEZrtj6IPkz09PXV8or65XeMMZ2yEbH2jOmvt4sw4c71QPbCa73O5XN2aErc983742S222EKSdOSRR0qS9thjD0lyx4ze9ra3SZJee626PMyfP1+PPfaYpGDMk36uwjajdhwdb16Xy+WcLycQb36jtVivzJUUzEnq6rR778w+Ifo8yP5HCuZDJpNJFWfGFM6+79c9U00U4w4C+b5f8jzvNEk3qVpNe6Hv+39sCauY4Hme8vm8c3BshtMMzkj29/dLqhprXJvfVsHzPHV3d7sNby6XS72TxTbWrVsnqerIGp1rThM8z1NXV5ezBxaJNCNqz51kG2G/UalU6os4pAiN5mChUEg95/U8JVXtuVQqpZ5zT0+PBgYG5Pu+8vm8hoeHU895ypQpWrt2rSTBP7WcG82/TlDVbLrppurt7ZVUXU94EEkrOnXd7kTOnboGEnzplL1zJ3LuRHvuxOeqThznTvUbnbZ37jTOE8GEdiW+7y+TtKxFXMZ6bUlBRoDvw2eQG0XKurq6ajZjSWbmmh3XGSmiF80YxBkEimbOceZsukdbAyc6znE6K96b673xjW+UJG27bfW0Ik5zt912kyRX90cKPteTTz4pSbr55pvdz5MIApFN+clPflLD9R//8R9rvm/WUSibzdb8Ls7FN1zvIvyVcceOd9ppJ0nSzJkz9atf/UpSYE886Iftua+vr6X84EEmmfvOXH/wwQclSQ8//HDN58jn803nYtSeJzIHwwUquV44ayU1v99jQZRz3PY80nHE0WQuovY8Eb8RrtVDhhNbCGerwsAfRO/H7373Oz300EM1rwFRe25UYLKViM5DvmfcGimWosjn85o2bZr7vlER/FaA9z388MMlSTvuuKOk+nXlvvvukyT94Q9/aGj7UVuOE3CD6yWXXCJJOuOMMyTJZcIZ/5122knHHHOMJGmffaq9Mb7yla9Iqn4eSS3nHlU7gonWqEty3Y6Ca+G3r7jiCknSm970JknSb3/724Y+pp2cx4u41kCpXnHCWsX12Bdtt912klRXqP7mm2+uWa9BNputKdqa5IM+19pqq60kSQcccICkqk1I0gsvvNDQb8TBmXGNvhd2x/5y5513lhSobeDxV3/1V5KqfuT//u//JNWuKXHYMzbBvT/99NMlSXvuuackaffdd5ckrV69WpL01FNPSZLe8IY3SJLe9a53Sara0sc+9rEaXiQj436uwha57vTp0yVJf/M3fyNJ2mWXXWo48xD/0EMPaeXKlTU/y2QyNQqh8PvGAWwFfx3dq0b3G83Wizj9xniAXU2dOtXZ+z333CMp2I+0cu88Fk7czxkzZkgKbJ/f77DDDnrxxRclSXfddZek4Hkrac7RZynm41ve8hZJcsmkKVOm1D1TTRTpTqcZDAaDwWAwGAwGg8FgMBhagnTrkxsgfJ5Wks4991xJQVSY7yuVygbrqSQBontE+VEpEe3j92E1E9HTpIv8ElUkS/H3f//3kqSDDjpIklwk/aKLLpIUcG6FYmGsiGa6yAAcffTRkoIsEbUc+GzbbLONpKo9UN+Be/HlL39ZUn2GOu7PwHXIxmHXZHSJTLcT0cwL58ff9773SQrUUzfeeKMk6Qtf+IIk6dBDD9VJJ50kSfqv//ovSaNTKYwVjCW2+LnPfU6SXHY+WgMIdQLKjl/84heSpOuvv74mU9RqMI4HH3ywPv3pT9dch8zbf/7nf0qS7r//fkmqkStLjZUr7S4IHu2exPlq5hbZr+hniQvM9+23314nnniiJOmII45wP2vEgc9AJvTPf/6zpGqGWapmYv74x+pp56THO5rZIhv+d3/3d5Kkt771rZKCbNa1114rKX4ftiHA9cMf/rAk6dJLL5UUqDSjSgPG9rTTTtO9994rqfXqmZEQ9cmzZ8+WJNeNBb+B7fzzP1cboR5//PHuPW644QZJwXxu9T2IdvP5zGc+I0mu09tXv/rVmtenpWZEI8CNvcSuu+4qSfre974nKbCV73znO5Kqny3cKSwNaNTBTKod96S4Yr8oTqgzyL4IP8HXzTbbTFKgXMH3HXnkkW4+NqrR1Q7gL/bdd19J0re+9S1J0jnnnCNJWrhwYcOxbzXK5bJTFqCOIUsPN8aX1zHO7P/pJPS9731PZ599dmxcQfi5AjXjaaedJimw39///vc1v3/mmWckBT7wi1/8oqTqHpr3iltxF639wz6Cff2XvvQlSdL+++8vKfAX4Y6MUnUt/+///m9J0oUXXigpUNDEZSvRPcK73/1uSdIpp5wiKVDhwWPVqlWSpCuvvFJSdS1vx7PVWMG92W+//ZxCjOeDpLt8RddvbJn1mWcsjsH/+c9/dj5y2bLqgSbWz7jXmah94DPw1eztNt98c0mB/+vu7nafi33eRPfTqQ8CMVhR+fOhhx4qKZhUTHwc7xlnnOE28nE/dETBDePowUc/+lFJQWBlhx12kBQ8aCxatEhSdVHAKSS1cWA82YDhJJFYLl68WFIQaOHnzz33nKT2dKRi8aFo3VVXXVXDjQnMES/G9NZbb5VUPdqxfPlySYETC5//TBJRu+Yz4bBeeeUVSe0JtrGIHnzwwZKCIw7IPqmf9N73vleSdMghh0gK5PyTJ092zmzJkiWS4pmLzLcPfOADkqRPfepTkqT/+I//kCT9+7//u6TA0RIoJiDw7W9/W5K01157ad68eZI2XJB7rMAeCaKdffbZTuZ5+eWXSwrG5ayzzpIU2C6bAjbl0VaUg4ODbXvQiwZXkZh/6EMfkhQcp7njjjskSeedd56k4Dheq4Ff4OjXhRdeqL333ltSsPA/8MADkoKgNkFWHtp//etfSwo2wczLdjx0Rm2QgNbcuXMlBUeNuA/YLgHXdnQ8wZ9yJJSALEfPOPZF4PXVV1+VFARU5s6d6wKgSSVxokeg2VvMmTNHUnAs6ZFHHpEU2Dmbseuuu07f+MY3JAX+OlyAtFWoVCouIHXqqadKkiuUjG8ey5hF17yk1j64RRN67OXwjRSE5jOG21O3K/AdbTKBXXPkBz/D2uj7vguuRG2CoEArPovv+zr22GMlBQ/4+F8exiju++ijj0oKfN6KFStqvj733HOpq18VDXDxPfMhKT9XqVT09re/XVKwbwDR5Ac2wPFRgI+45pprEvNxvD8BKxIG2MRPf/pTSYF/Zm4m/TwSBhyYW+94R7Xb9VFHHSUpSDA8++yzkmrXPSkIvu23337utT/72c8kBfckLjuP3k/W6n/5l3+RFCR5Obr2iU98QlJwHw444AD3uZKci80Cb+EyK2HgB2fNmuX2pu3yzSQZ2UPPnz9fUnV9loIEDc/d9913n3tmIdFw0003SZJLRMWVxGHNI+hEwpq98cUXXywpaBBEsHvmzJmaOXOmpCBhPFGkI51iMBgMBoPBYDAYDAaDwWCIFekK94cQjSYSOSV7T9aT6GO0wNbAwEDiWVCUE7NmzZIURPmIJhJlJMuM+gZJ9+23364777xT0sQLPI4EIrxIWlEkkCk44YQTJMkVriOqzrjffvvtkqoZ3aQkw83k+hyPoEgg6iWiqmSbuT/9/f11ipSo/SQFMlnI+VFYcfSgHSoPFBXYBllusnBE2ImWU/R1r732khRkbk855RSXuY3zc8CXIofcZyLlZAg4/kXRbcacTP9xxx3nJMatHH+i/2R9ttxyS6eOOf/88yUF2R7UbHAkw08GC8kzr//2t79dV2QwbjDe73//+yUFKg7GH8UdGRiUWXyWOXPm6KWXXpLUmiwXWSsyhhTy3XXXXZ3MGjk7KhOOeUULr0dbyYJ2Kh6PO+44SdLHP/5xScH8IvOJgoJMaTuPQcMZ1RJZKzLMqBRQI7AOwfnUU091WcekmjZwbdYEsoN8FpS6KINQ+JLV5RispFj9XalUcipLFFZkjr/+9a/XXLfZ9cPtwMk+PvHEE5KC9TJu1Wl0Hceuke2zfl9wwQU1r2vH0ddwa14psMl3vvOdkgLlNMWr8WuofyqVihtf9n28BgXzRGwlPNexS/whPu7555+XFByfufvuu2veI3r0IZfLpe64HT6BPQn7Uuy/q6srEX/neZ7b96DORW0a9WWobDguz/ENFFjDw8OJqDzCbdsp54DSIHpUm+cOfo4KHHv4y1/+klhB8CgH5hpq3QULFkgK7Jn9P3bAUZpZs2Y5zqjdkla6sSdFfQmwBZ5hUKu/8sorbZmDjDnH1T74wQ9Kkq6++uqa36MA4uTNgQce6PxZ0qcWWKfZa7BOY+OsI9g49nHYYYfpk5/8pKTgc6HWjwv4WtTynA768Y9/LEm67LLLJAVrMvEN1uodd9zRnapolQ2nw9MbDAaDwWAwGAwGg8FgMBhiRWqVQNGWkGTrUQLxc86sTp06VVJQ4OzZZ591Kou4QXSP4qNknYn8E+XjjCdqGwqaUYMgqQhqsVh0KiQUQJxTprgXWSSiphTQIluHUoJsdNzwfd+NMyolzsCTeaHNOplbCh5GW5jHrbIaDfgsRNLJJH7ta1+TFGTvks5YlEoll0HBjomOc24ZRRtzlAg7RcUp7Ln33ns79UtcGdxwnQyUSMwjMnF8z1d4U78BtUhcYx3OoknB/A9zItOJao16AygDKbBNgUE+25IlS/T000/Hyh9gsxRmxFZR4Hz2s5+VVC2wLQVZL2p8MP/6+vpamuXCtpjX8Mnlci4jz5izXkSLVCddN24swPdGs+JkoFFkUdyQzFiSn4nMIJyoBca6glKFDH60SQL3bKuttnLtUf/nf/6n5jVxgXGF+z/8wz9Ikn70ox9JCgrgU9uPWjVknnt6emJVIYRriqHKpKYFSsKR7jn3B9XuUUcd5fZU7E/OPPPMmuu1OhPNGOEXUHVQHJc9BzXcqKnCvG6HGo9rYiPs8bAB7Jb7gsKT2n677rqrew37PxRBrUDY7vDHqD1QwmIr7I2pSYkPZF1KUx2gqDqPvR5ZflSmrOFJ7emy2axTa//pT3+SFNhx+DVS8KxCDUIK5qL8/93vfhc/4Qi459G5FK21RU1KlMmsnwsXLnTzN+71BU6sGXBgb8weCZUN+338Fmv9okWLalRuYSTlU5opNOFDfbnbbrtNUtWu+TxJgPWB+46KBt+Fygc/yHyjCcymm26qH/7wh5KS23dEFc88mzJu1ARF5QMvnnE///nPuzUcv0KtyLifxdlL0KL+Bz/4gaSgFijji8qePcill17qVNSmBDIYDAaDwWAwGAwGg8FgMIwa6Qn9R0A2inbfBx54YM3vOWeLIoWK8UQ0pWRqIhSLRXeOk6go9TvIyBCFJnJ52GGHSQqUFlQrf/rpp2ONQHL9TTbZxF0bUMOBLAdRU6L+KIeIUJKlC9cYiJNzqVRyNS9OPvlkSUHUmkwXUVUyXUSAySCkIdMVzSxi39gx2Tp+345uZShq+EqUGgUQWXzuPTZCFD3czpyoPGqdVtcryWQyrq7Lww8/LCkYU+wDpQxzC94oVDhz/uijj7qMeit9B/MfVcSDDz7ozvpG6/lEVUv33HNPzd9SfwR11tSpU2Ofg4DroFxkfFHm4ZO5xygqaJ1Ldmnt2rVNu02MB4wV85+s/OzZs50Nc6admkC0jEVpRZ0HlBJpqInB50KpgfqVmkBk6cJdUNoB3/edn6JOGx3jUNOQ4Yxm7Bln1Ifd3d2uHgH3JimgBEJVzNyjyxIKRxRK2HDc+wzm3a677qptt91WUqCeia7XvBY/xvesidSXmjdvnlvD6dTGehlXjb+o76erJEo3ajewnrN2tEMBFFUtUXeLPR3+gxqJqEKwZ2omht8raifMmVZ9PtShjBvXo4PZpZdeKinYf1KPhHbrqLrTpJR+85vfLClQvDG+1PdL2k+XSiXX+ZA9cbReFHs4FCyAZxgUCnPnztX//u//Skq+M18zsA/9p3/6J0mBT/zNb34jqer7khpzxoR9JSo71Hgokun6xP4jut8P+7OkfAnPodhx9LmO37M/Yc3jFEMr90cbAusE6y+KXfbErBfYNmsE3Sk/8pGPSKrWx0OhlZQSiLFlveb5EEUSa2O0Dg+1pBYuXOj8Cmqy6KmRuBBdn3new06Yf//6r/8qKVA+fve73235/BvR83iet1DS4ZJe9n1/t/U/mybp55J2kPRnSUf7vr+mpcwmAI4jDAwM1Dx8Dg0Nuba5cUu4x4rzzz9fK1eulO/7jnOlUtHAwIAqlYoymYwmT56ciocTcMMNNziO0XGGc9rGeWBgQMViUZlMxjmPSqWivr4+x3nKlCmpGudisahKpSLP82qk8cVi0dlzLpdL1TiXSiV99rOf1bRp09zDE7bi+76z5zRxXrBggXp7e51/kNLvN5588kkVCgV5nuc2O2mfg77vu9b24XEuFApunPP5fKo4Dw4OqlQqyfO8mlao2LPneamz53Xr1qlQKCiTybhjmZVKRevWrVO5XFY2m9Wmm26aKl93yimn6JprrqnZBKd9Dvb19XXUmrJ06VJdfvnlmj59ujuSUigU3HpSLpdTN8ZDQ0Mql8tunknp93PDw8Mb5JzNZlPnM1ijPc+racc+PDzs5l93d3eqOJdKJRWLxY7ifPbZZ+uuu+5SqVSq2dOl2Z4vvvhivfzyy5JUxzmtvnmkOZjGcX7hhRfU398vz/Nc85JKpaLBwcGaZ8E0cb7tttvcPjScBE6zbSxZskR9fX3yPK8mIB/2G2njHAdGE36+UtKlkn4S+tm5km72ff+rnuedu/77c1pBiMgYxk+XHH5OZwXUCdTK4PznE088oS222KKuk0SxWFQ2m1U+n3cbnlZkvCqVisvQrllTjYPRsSCqKqCzB9FWujmcc845mjJlinK5nAqFgjNIugf09PRoaGhIw8PDbpEbL1epmnmj6jt1PKhGjqMkm0QdGM4msjkvl8vabbfdXEYGFAoFZbNZTZo0SYVCQYVCoSWZpXw+76LPZGaJSpNJZFwBmS42utdff716enrU19dX03lmaGhIuVxOkyZN0uDgoIaGhmKrJxVWY0lBxJ17Q3YWhB1UWOVWKpWUyWTU1dWlUqmkUqk0oQh8+Aw+Cg4W/CuvvNL9Dk7hr9RdueWWW9zPs9msTj/9dF111VXafPPNNX36dL366qtuDhaLRQ0NDU3IngHjQwYZG6UuCR0h8BEoWFBU7LLLLpIC22XeSq31G4zftdde67r7UBOKrEU0I8jc4fdkDsnAzJgxQ2vWrKnp2BHHHKxUKs5mycTdddddkqTly5dLCjqX0Q1s7ty5kgKfwn0gyMP8lQJ7zuVyKhaLY7Jn7JAxIJu2fPlyV5dhjz32kCT97d/+raSgUw51SH7+859LCtRKZNXDtQTy+XyN38A/d3d3a3h4WMPDwy3N4HFt6kSdddZZkoI6DdSawBai3S96eno0adKkGp8yODioXC6nqVOnamBgQIODgxOuPVAul51PRvWFP+GMOxyi6o6oiu34449XqVTSr3/9a/eauNbuKKJcGE+6fBL8i167u7tbPT09TpEotXZNwa733XdfV1eGWgYAW2Geo9ZjfUEl9pa3vEUHHXSQfN/XEUccoU9+8pO6+eabVSwW1dXV5Xxzq9ZtKbBH1i/8BDUcUHBS44PPwuu7urqUy+WUy+VcxlaKb68hBXs39mjRWkzMRRRArGPhhwnP8+R5ngYGBtzPBwcHa/Z0rVoDw9dmrqGM/MAHPlDzc9RL1LKhltCqVavU29urlStX1qwnzD98c1zzrxGYU/gCVHjUKGEPlM1ma2wjLs6ZTMbt89n3wJE1ne5VzFXGku6Z3//+93XwwQfrXe96l8455xzl83mVy+VY7TmqgAknjMKfgVMNdOyD+3XXXaeuri5lMhkXwJTi9c1RNR7+AdUM+1R8HUogbOPOO+90PiNpv0HXJ5RsqEqpmwlYw6lL+dRTT7m9FqosqbrXyGaz2mSTTVrmN8rlslMTs+9BJUjXLAQV2AGfi5MkdKp6/vnn1d3dXbOnS2rdBtEuvKwf7El4PsROHn30Uff8xF4T24ibM+PKuKMqpPYPaiVUupwy6u/vb7lqcMTUlO/7t0taHfnxEZJ+vP7/P5Z0ZEtZTRDIBMMolUpu8LjxaQIy7zDCBsjmLE3AGYeR9nFuFNkNLwJRR5YGNMoghxfibDab2JGg0YJsRRjhzXA+n0+dbXSi30DqH0baOTfKrJTL5RrOSbWhHS0atSAObx5yuVzqxrmROrBQKLhAVU9PT+p83Xve8566QFra7bnZOKd1TXnf+97nNsqgUqmkeoyz2WzD+Zdmzs18RnhPlzbO0YdUKf2+udH+KO2c99prL5foBmm3Z5IPYaSdcyNFcdo5ZzKZjvMbjQITaR/nRut22jnHgfGGlLb0ff8FSfJ9/wXP87ZoBRmkyJJ0zDHHSAqieGyiFi5cKCnIiHM2nqg853G7urpqgiYcQZFUpxKaCDzPc1ljHsjIEv3qV7+SFETdzzvvPElB9p9aJWQOoosZ0r9WcQ5H1uFMVp+6EgSjiAiT9Vi6dKkk6eijj5YUdJSI8oprnCuVih5//HFJQYSXz0A2PyxDlKR3vvOdkoKzoHydN2+e7r//fpXLZccxDs4bAgGcsNxT0qg3K2Qapep9nShnrnvwwQe7zMpll10mqb6mDmDMuA/YCg9xm222mTt+0tvbq1Kp5H5HRqkV4L6jTHnggQckBedpycZQW4XK/GTJibpff/31ymazTj4stdaeGb/777/fKRqZZxyZaxbl5/4QeOX89erVq+vGMs45KAUKiZkzZ0qSvvzlL0sKxpN5x3XJeqBeyOfzdZxaYc9RBcEzzzzjFI7U6EAlRB2jQw45RFKgeKSeFHUGqJPxzDPPjOif4wpp/MmCAAAcRklEQVTERju3UEOO7j+M+2gQB2ff910GmbUNhQRrWzSwgz3TMQffPDg4qMcff1zDw8M1tfSS9s9SsF6ScR1LRrCVnOHx/PPPu9o+dNBCacUDGuoO/AS+MPpwv2DBAq1bt05Lly5VT0+PfN93dtaK9SQM3ouvqPFQYKLUxQbo+EQXFNafqH+I0y54L/wxtcR23nlnSdKtt94qKainQ80G1LPNjv4l4TOiSl2y/dQzYh2key1qJxRCc+fOVW9vr+Pa6r3GaMDDF+p5QF0jlLHh48RhxMU5m826/Q7daFH+448BfoN7jDKdr9tuu61TepfL5VjXbfZIKDzYAzHXsGfsG+UsdY1QeXd1ddUEtJOYg/glFIN8Bro8sYbwLMN6iHIIpVvcnMvlslOZM77UVJo1a5akQEmIKh11GF1fm9WEjcNvhOtbMXY8S+2zzz6Sglo51Az6t3/7N0lBh0PquA0ODtb5vLjX7ahSDOUdaiViAXCmxmBYVZPNZhOJE0Q5sxdGHX/JJZe460vS4sWLJQUKVJ63G+2dJ4rYq5F5nneSpJPW/z/uy7UExjkZhDmnpV7CSAhz7hR0mm10Gl/JOCeFTufciX7Oxjk+dPo4G+f4YJyTgXFOBsY5fnQaX8k4txvjDQK95HneVutVQFtJernZC33fv0LSFZKUzWY3GMIKRybpxkANAeoiEPkngspZVqJ5RFij8DzPRVMpsjsazl1dXRvknMvlXMaKjPicOXMkBTVJ4EgkEsUEmUUyA42UFuPh3Gyc2Rg/99xzOv300yVJZ555pqTgfD4ZXDJcv/zlLyUFkV8ixM3qdbR6nMP1Pr773e9KCqKn2AScQ+8rKag98JWvfEWSXB2kRYsW1SiBxss5k8mMOSRLBJ/6Etg3EXiUCyOd+yS7Ff46Gs7Nxpn5c9hhhzl7JYtPVq7ZuemwakYKbAMlUG9vr1577TX3+7Bdj8R3JJ/B+0lBvTDOs5MBwLY5Q87YkoEmY0RdmuhnG41tjIYzf7tmzRpdeOGFkqT58+fXcPnjH/8oSTVnlKUgw/+mN71JUpAdX7t2bR2n8XAeyZapbSEF2Sw6MpJJ5ufUCKJbA91QsKtcLleXzRitPY9lnLu6uuo68JD5/M53viMpyObiq8ngUS9q//33lySdccYZevbZZ2t4R/3zaOw5l8tNOI3DPMSf4TPoILehAMh4OI9mDrKGoVCjkyDqzaiaCX/DOoSS5cEHH9RDDz2kwcHBmu5bY7XnkdbtRmAsojWm8NlwH82mr5V+A19www03OGUa2U5AzQn2GKx5rN/UyUDFcN999zkVArVrWrnXCIN7/f73v19S0N2Tn6Mio/YS486eA+6LFy/W4ODghNft0dgGY476gPHG16EYZJx33313SdUuOeHPFsVo51+U82jGOQpsgXpQfI9PJIOOChWfuMsuu+juu+92R1DGu9cYD2fmGD4OlccjjzwiKejKFu3MN971ZDycuWdR5U+zPTHcOIKJUvaXv/ylfN+vqR3ayjnIdSdNmuRqz9B1iFML2CnrXNQ/oxqjvlR0Xxqn3whfQwqel1B8sEbj46hvSg0sur92d3fXNFGJy29ks1nXmezTn/60pGDPw3Mpvg3VHZ2kqa3arGNgHOt2Pp93Ckdq70afXfF31Kph/UCNxd6ju7u7Tp3UyjWwEXg/xpy9Bvv9yy+/XFIwTzk1xD42l8u1hPNY9hq8H3bB+PFczb4J+wDRUy6txHjTZUslzV3//7mSftEaOvEhfL4vfO4vzaBQtFR9GEyq9d5E0KnjjOxxeHg4saKHE0Emk3ELCpvitMPzvBp77gTb6ER77kTO4YLcpVKpqSw6TQgfOaawbtqRz+fdw9/Q0FBH+LpOtOd8Pt9Ra0onjnGncu60NbATfXMncu5UezbO8aMT/UanjnOncZ4oRtMi/ipJB0razPO8ZyV9QdJXJS3xPO8Tkp6R9NFWkPF930X4qTtBhpEHXrLPRCpReyxbtkxS9excoVBwrx8YGHCdJYaGhlz7y1Z1cfE8zxkNlbypPk6dHSZvNPOCgqVUKrmWhb7va+3aterp6VF3d7cGBgbU29vbsNDueJHJZFz0nOwQmRci1kR8w4V8pSAiOX36dF133XVunPv7+5XP591DRn9/v2u/2AqEo+ycZSezQnY/imjtjD333FPHH3+8li9f7lqh9vT0uM4u1KFApdNKRLvhUHeHbAVnQFGzhDcthULB/d3Q0JBTOBSLRQ0PD8vzvHEHCHlfzjLvttturs4EmdjROkJsYYcddtCTTz6pk08+Wa+++qrLLlA8d926dfI8r2X2HD07TqX9gw46SFLgS6gLQx0mzmGfeOKJkqrdA3ivOP1GPp/XH/7wB0mBagYV4TnnVJsscgaYuhnnn3++pOBsNpm5F1980fmXOOcgqgEpyLBQ/4K22MxLxhMfQkcX7k/YP9Otik3O0NCQPM9r2UNztB5JNFhKLblonQwyX6jEdtllF1eHYt26dTXdwmhPHFdHQcBcpePZ7NmzJQX1YRhvalH09va6DOjq1as1efJk1y1s9erVrkV8K0DdkenTp0sK1sBwoFoKbAAFEPW6yJCdeuqpbg3CNuJau6PA77GeoEK47777JAUZR2yIz4IN+L6vNWvWaPLkya4LZSvWFMZmYGDAqWioa4Bdoy4mwwhXFFasM0uWLNGPfvQjpyTEd8S1boeBH6abKpnbz3zmM5KC/RC29La3vU2SdOyxx+qCCy5Qf3+/fN/XwMCAuru7Y+UcrQmE6vjGG2+s+YoNsMaxjr700kuurb3v++rt7a3Z07V6DQwDTux7HnzwQUlBfSj2Cuxb8eObb765Zs+e7V4f7tg4PDysUqkkz/Na1kmpEbBH1Eso/n/zm99IChQq4f0ObcClWp8RF2f8MLbRTC3Az6Nd8WbMmKHZs2dr6dKlTikdhz1z3eOOO84pN1DTU1v12GOPlRR0/8Ue2TtTz69UKsn3ffeecT9XSc1VXtg36xwqSNQTqFd22mknzZ8/XwMDA/J9X/39/bH7DdYGOvHxGVD8M76HH364pGBfh91PmjRJAwMDbrzj9BvhpOw111wjKdiPsv9CAYkSHfUV9W4Zu/DzaxK2IQX2zXP2Ntts4z6XFHR8ZU/N+hNWQibNmfGmVhQKILhh29xfYgZxYsSnO9/3Zzf51eirUSYMJKxR59yqVpxxoKenp+GDfBwBiVbhqKOO0mOPPVbnkNI8zj/96U+1bNkytyiyeWjUFSMtaPZAnObs8owZM7Ro0SIXbORoVprtOZPJuM1i2Hek2Z633357dzwy7D/SzLmZ3cb5cDFRXHTRRU5OzEYpk8lMuL16nGjm0wjapRFz5sxx9twpc7CZT0vrmvKxj31Mjz76aN3xtzSP8bx58/Szn/3MtSTG16WZc09PT0MFSprXwKuuukof//jHXbFjHpjiCry2As3WjTRzvuqqq3TmmWe6hArjnGZ7JtkYDX6lmfP555+vVatWuaBWJ4xzswBPmv1Gs7mW5nHuRM5xIFVap1wu5zaAZOGiN4ToH0ogamuwOSC7keQDcrgGhRScOcTxwImOHUQuOeuedGGpcFY/qqKJOnaCJGSfOUvMfUpKLlepVNxiTycfnOU3v/lNSfVZfrhFaztst912bnMWPo8dJ7g2DwbYL5vZq6++uoZrUkf/yLKgCJgyZYq7t6hRRsuFsSR7N23aNKfQi9NOGNsd1nfOoRsE95QuYWxsd9xxR0lBRxf4hivvx20P+Kcf/vCHkoLxQSUID7LhjC1d+jjDP23atNi6UkXBmKCegSs2i9qAmkyMN79P2s+VSiXHCXXBbbfdJimw6ajPQEFBPawvfelLkqr1Bghkogqg+0RSc5W1j+wtvptz5O08LkwmiyNQ2CScWT+oPUDNBHw6ypTf/va37m+TOsLBdbCNK664QlK1DpQUqDtQN1FHATXemjVrEuGayWTc3iJaBy8a0MG/oKggA4n/KJfLiY0vdklHVzoO8XPqoFEXj3WHPSDqx5UrV8bCjzHDx4br+UT3P/gTOpkdd9xxkoJuLyix2lkslGszntTpCnd/CwPbQYGHf0sSUcUmXWp58EXFm4bjGeyB8Xlk66Pjii9k38+zCnP4hRdeiLXTGn50t912c/OevfJee+0lKei+Fn3OIlGAwg0175o1a9xzzWg72Y4H1GaRAvvk3kfrbWHX7P/Yy+LPn3766cTtJrqvwP+yduNrOInRroBldL1ATYya6tRTT5UUKIWoYdTOBDT3HxUx+318MqpY9tb45HZy5n7TGZNacuwhqZtJXSs639ENMU77TX8hEYPBYDAYDAaDwWAwGAwGw4TR/rB6CFTyl4KIJFE9QBQQGT7ReLJLaWjXRhSYr0S0yepTy4a6IO3kPFI2kGg/mRkyt0S2k0KlUnHZCerpYANEUVEoYENk7ej4QyR42bJlLnuTVEFlOJHZYjyJELcL0c5aTz31lOvi8uEPf1hS0H0hWtcDMCdR2MydW60Z/9prr7ksb5zjzPyJnvkl64V9kO2imx+1D6ibMW3atFF1WGolZ3zDZZddJimoP8LY33XXXZKkO+64Q1LgOzhLTD2aJBH1b3wlA8dnYpyxj6QzMaVSyWVhv/jFL0oKav9QY4JOI2Rtmadkj77//e9LqqoWyOTSfQLbjltNyLyj3g7+DyUpmf52Hqdj3uCD6e5JJxwyYKyB2AK2QtZu2bJlbcv2s6ZdfPHFkgIboGbNIYccIilQqjAnzzvvPKfKiltdE51zUWC/2Ay+nLWT/VQ7GgnACdUM9SWoj4ECFZXYRz7yEUmBovCxxx6LxTaY+wceeKAk6YQTTqhTl6Lexv+iQsZv85U1qJ2NGphb4fp2UmCbUZUeyhbs/cUXX4xVobIhoA6bOXOmpED9dcstt0hqnxKIMdt666319a9/XVKwlqCywf8ydnwG1h66ZaL6vuWWW2LxF1yf/fGee+7pfBZHmlFkorJj7vHc9fa3v11SoIRF9djf3+9qTVHHMI513fd9Zwv4AdZqbIHf0zEV1QrrJHV50gDmGnX8OAbYbE+dNNjDsBeOdnSlBh2+sl18S6WS4/iNb3xDUqDUxS7Za7SbKyiXy27dQF3MekxNKPZNp512mqTgPrCOxHlCIVVBICn4kNEPixNmwaLdMxM9qWMc4wFGyIYYOSWy2zRyjiJqjEywpLhns1lXwJCCXwsWLJAUtJZFsrj11ltLkvbdd19JweLLInLnnXfWSU3jBtdjgUKOy8/ZgHIMJamjHY2CEUjff/CDH0gKFgTas65atarmb9gg0x6X41Wf+tSn3PGaOD9PNMDGwoocPmqjbI4IxLFJpohkkuB6+Ahk1xydAMw3jh/wfbNWxEmC8Z4xY4akoO1luHZOO9DT0+OO8BAApA0ym1t+j8yZTXC0re8+++xTtwax+eUexGU7BCcIDGK/PBy3u2tk2Dfz0MbGneLsPLwxhmw6Fy1aJCnwMdlsNrE5CBc2ZPgDQGCHjTBHeRnvZsdA0gA4sT9ivNu5IcYPkFQ47LDDJAXF1/FpPHDykEkx27gC3tgmQevFixdrjz32kBQcu2MOEvjjtawd0WLB7QBcCDoQtKToL22o8YUcBWKcObq9YsWKxH0KawiJJPb3d999t6TgIaldc43rvvbaa26d4+gGx1L4DOyH8IEkI0kaXHLJJZKq/iXOcYbziy++qGOOOUZSsK5hz+yl58+fLykIYvIw+t73vldSsNZtt9127qgQdhRXEgRO1O8kgAVHjovSEp5gEEXbCcy20z/j0xhHbCQNR7ilwB4IYF500UWSgjWRwtAc127XsbVw85dzzz1XUlAIGttmb8wRaJ5T2r0+l0ol5yN4Nvrc5z4nKfAJrIUkUtmHJrGu2HEwg8FgMBgMBoPBYDAYDIaNAKlTAjUDkcAtt9xSUrX9nxRkjeIsUjZRkH2jRSSFJckwtlM6PFqE29RKgZopKe7haC4tIWnZjMTygx/8oKRAHsjf0M6Qgp/33HPPiO09W41o60jar6MqQPnRrkwtGYk777zTFd6m5SYyW455kR2NFn1FIk9m8aabbkrk80RbxKPeIPvC/Sdzf/TRR0sKWsnDd/Xq1W3PzDSTVmOv+BDmIS3V2wHGnTEjI8e8pOV2O7u0MO84/oUSiKwRR6sOPrja7JJ1JGpT4aNW+D4yknEfAyPjxTE0CnKjvml3l8BMJuP8GjJ3jn2RnWU8ySiiMkTNRGY6ybWQLBtHIs477zxJwVFLjshybJRjCPfee6+k6vEkqToH2y05bwZsI7reteO4DxzIfjIHZ82aJalebYOtRItetxrcO9QmN998s1PkNnttsyOx7QTjS7FZjtmhkOYzsT6ilH7kkUckBQVWBwcHE1c0YY/MRTLnKGOjpRXahUKh4EpPHHXUUZICRVV0bqGQefzxxyUFDSoeeOABSfH5bXigQlmwYIGbU4wjCm2UsBwX4/fYAqol1vDtt9/e7T2Yr3EdaWN95Z5zxBjFN0dmOA595ZVXSpK+9a1vSQr2IblcLnF/F90/0CyDe4/KtN2FzrFR9vwcH+a4ILablk57XV1deutb3yopUCtxxArl46233iop+dMqzdDV1eVUSShIKRCOMggfwl4DhXQS60r7Vy6DwWAwGAwGg8FgMBgMBkPs6BglECCyynlUzjGTtUtDRgaQyd17770lBRnGk08+WVIQhW1nQc/RgnEn6x93FrwRiJpTc4SCnSeccIKk+ppLZP8p7kq2o7u7O/HoMONHnZoDDjhAUpANQhnU7sxAd3e3q+uBAggFFcoasvsof7gf1157raSgNaeUzHxkbOHB/SYTetZZZ0kKiqKSwaPeAGqErq6utmcNRgK+Ilrksx3Af+28886SAv+GHaRhLKNZUVpuUpeEbDg1ut797ndLqlcEeZ7nMsDUD+LzxzVnuTaKCZSjS5YskdT+godh4MdQaVIfDJUdGf2FCxdKClQe3J92rNuMH+sFGUUyzZzLx3ai5/RZC9N0H0BUOY1SC+V0O8YbTtgKezbUYyBp9VK0nmQ2m03lPR0JjBuFt6mrwzqO4gowF1EWUty9HWpYxh41Ftl8MuPtVgCFaydSk4RalBSi5St1jfCFF1xwgaTAzySl3MSGV65cqYcffrjmd3we9hPRNYw50dfXJ6nWf/C3cdd5pBYYiksKbVPzBftln4/KAoUS/Nqheow21IEzc5Hft1N5Hm6cccQRR0gKCkCjEEvLsyk2VygUdOaZZ0oKuDLvUPxTuD0tRbez2ax77kMNOGfOHEnS/vvvLyk4CUJTClRs0WLzcSA9ERODwWAwGAwGg8FgMBgMBkNs8JKMkmazWR8Fz1hBJoB6L9RHIGJJ54OxRtkHBwdVLpebpqy7urp8KnaPFUSkaXlK1pnzoePNxq1bt06lUqkp52w266PSmCiIppLt59wlLX0HBwdHlfEfGBgYcZw54zkSom1wR7JhxnmsEeG1a9ducJwzmYw/1kh5M+6tqlE0PDysSqWywXEm0z0SmHOjrbcV7XQ12s/S19fXdJzH4jOiZ/G5fjTLxefh80XrPIwG/f39v/d9/x3NOLdq/gG4ch6azNLVV1896i5FfX19TTlnMhl/rOe+qQm1++67SwpaltP54Prrr5c0sWzX4ODgBsd5tLYMsBHGczwZ5pHspbe3tynnXC7ns4aNFmQPyRZR86OVGeVXXnllg+M81nU7OsdAtIbKRPzdhuZgK9aTaKfSVnBevXp1In6DuUnNKzqb0UZ3tOq1kdbtOHzdRDGavcZY/Ubc2NAaKFXHeaz+udkcjGK8+6OhoaERbWOsnEe7ho8Xo+E8kj03Wzua7UPHO75gNHNwvM9VcWGk56qRxtnzvFHvP5v557EqCFvpN+BMp2JUeIsXL3bXCnMeL1q1d272PBIHRto7j8R5JL82nv38SBhprzFWu4B71EYn6ivCWLt2bVPOYZgSyGAwGAwGg8FgMBgMBoNhI0DH1AQiQsa5/AsvvFBSkA1td4eURiCaSo2U5cuXt5HN+EBEcsWKFZICxRXR2nbU/Wh1dihJdBL3NHU+GQ1Ge1Ydm273WeGxAK6ce2ce9vT0uPuT9Nl3bJh6VieddJIkufoDabTxqFqtE2yAtQ67TuNaF0UnjW8n+eTRAlu58cYbJQWZx7R0eTHEj06agyCJejMTRaftizoRvu9P2C+3oxYQiNYAveyyyyQF/jdtttNJPqIT/RpII/d0WaLBYDAYDAaDwWAwGAwGgyEWJFoTyPO8v0jql/RKC95usxa9z/a+72/e7JfG2Tg3+2WHcl4naUULriMlwLnFYywZZ+PcHMa5AYyzpI2Xc5JroGScG8I4256u2S+Ns+2dG/0ipT5DMs4jIRF7BokGgSTJ87z7RlOsKKn3SfJaxjmZaxnnZK5jnJO5jnFO5jrGOZnrGOdkrmOck7mOcY7/OranS+ZaxjmZ6xjnZK5jnEeGHQczGAwGg8FgMBgMBoPBYNgIYEEgg8FgMBgMBoPBYDAYDIaNAO0IAl2RsvdJ8lrGOZlrGedkrmOck7mOcU7mOsY5mesY52SuY5yTuY5xjv86tqdL5lrGOZnrGOdkrmOcR0DiNYEMBoPBYDAYDAaDwWAwGAzJw46DGQwGg8FgMBgMBoPBYDBsBEgsCOR53qGe563wPO9PnuedO8a/Xeh53sue5z0c+tk0z/OWe563cv3XNxhn42ycjbNxNs7G2Th3EueJ8O1EzmYXxtk4G2fjbJyNc/Kca+D7fuz/JGUlPSFphqS8pAcl7TqGv3+fpLdLejj0s69LOnf9/8+V9DXjbJyNs3E2zsbZOBvnTuE8Ub6dyNnswjgbZ+NsnI2zcU6Wcx2HON889KH2k3RT6Pt5kuaN8T12iAzUCklbrf//VpJWGGfjbJyNs3E2zsbZOHcK51bw7UTOZhfG2TgbZ+NsnI1zcpyj/5I6DraNpFWh759d/7OJYEvf91+QpPVft5jg+0VhnKswzvUwzlUY53oY5yqMcz2McxXGuRZx8JU6j7PZRT2McxXGuR7GuQrjXA/jXIVxHgFJBYG8Bj/zE7r2eGGck4FxTgbGORkY52RgnJOBcY4fncZXMs5JwTgnA+OcDIxzMjDOyaATOdcgqSDQs5K2DX3/RknPT/A9X/I8bytJWv/15Qm+XxTGuQrjXA/jXIVxrodxrsI418M4V2GcaxEHX6nzOJtd1MM4V2Gc62GcqzDO9TDOVRjnEZBUEOheSW/2PG9Hz/Pyko6VtHSC77lU0tz1/58r6RcTfL8ojHMVxrkexrkK41wP41yFca6Hca7CONciDr5S53E2u6iHca7CONfDOFdhnOthnKswziMhzoJDkeJHH5L0uKqVtD83xr+9StILkoqqRt4+IWm6pJslrVz/dZpxNs7G2TgbZ+NsnI1zJ3GeCN9O5Gx2YZyNs3E2zsbZOCfPOfzPW0/EYDAYDAaDwWAwGAwGg8HwOkZSx8EMBoPBYDAYDAaDwWAwGAxthAWBDAaDwWAwGAwGg8FgMBg2AlgQyGAwGAwGg8FgMBgMBoNhI4AFgQwGg8FgMBgMBoPBYDAYNgJYEMhgMBgMBoPBYDAYDAaDYSOABYEMBoPBYDAYDAaDwWAwGDYCWBDIYDAYDAaDwWAwGAwGg2EjgAWBDAaDwWAwGAwGg8FgMBg2Avw/aD9qyOqYgW0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1440x1440 with 25 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Visualize data\n",
    "nrows = 25\n",
    "\n",
    "def row_to_square(X_row):\n",
    "    return np.reshape(X_row, (20, 20)).T\n",
    "\n",
    "# Choose nrows random rows of X.\n",
    "indices = np.random.randint(m, size=nrows)\n",
    "X_sel = X[indices, :]\n",
    "\n",
    "fig, ax = plt.subplots(1, nrows, figsize=(20, 20))\n",
    "\n",
    "for i in range(nrows):\n",
    "    ax[i].imshow(row_to_square(X_sel[i, :]), cmap='gray')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Network Model\n",
    "\n",
    "The images are of size $20 \\times 20$, so there are 400 (n) input layer units. The hidden layer has 25 units, and since there are 10 digits, therefore the output layer has 10 units. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "input_layer_size = n\n",
    "hidden_layer_size = 25\n",
    "num_labels = 10\n",
    "\n",
    "weights = io.loadmat(os.path.join('data', 'ex4weights.mat'))\n",
    "theta1, theta2 = weights['Theta1'], weights['Theta2']\n",
    "theta2 = np.roll(theta2, 1, axis=0) # Matlab legacy\n",
    "nn_params = np.concatenate([theta1.ravel(), theta2.ravel()])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Cost function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 155,
   "metadata": {},
   "outputs": [],
   "source": [
    "def feed_forward(thetas, x):\n",
    "    \"\"\"\n",
    "    Propagate a the input data through the network. \n",
    "    The thetas is a array of weights for each layer of the network. \n",
    "    \"\"\"\n",
    "    a = np.copy(x.T)\n",
    "    az_arr = []\n",
    "    for theta in thetas:\n",
    "        a = np.insert(a, 0, 1, axis=0)\n",
    "        z = np.einsum('ik,kj->ij', theta, a)\n",
    "        a = expit(z)\n",
    "        az_arr.append((z, a))\n",
    "    return az_arr\n",
    "        \n",
    "def map_y(y):\n",
    "    \"\"\"\n",
    "    To train the network we need to encode the labels as vectors containing\n",
    "    0s or 1s. i.e. each entry in the original y matrix will correspond to a \n",
    "    vector, containing 0s or 1s, of length num_labels.\n",
    "    For example, in this case if the i-th entry of y is 5, then the \n",
    "    corresponding vector will be the transpose of [0, 0, 0, 0, 1, 0, 0, 0, 0].\n",
    "    \"\"\"\n",
    "    m = y.size\n",
    "    mapped_y = np.zeros((num_labels, m))\n",
    "    for i in range(m):\n",
    "        mapped_y[y[i], i] = 1\n",
    "    return mapped_y\n",
    "\n",
    "def unroll_params(nn_params):\n",
    "    theta1 = np.reshape(nn_params[:hidden_layer_size * (input_layer_size + 1)], \n",
    "                        (hidden_layer_size, (input_layer_size + 1)))\n",
    "    theta2 = np.reshape(nn_params[(hidden_layer_size * (input_layer_size + 1)):],\n",
    "                       (num_labels, (hidden_layer_size + 1)))\n",
    "    return theta1, theta2\n",
    "    \n",
    "def cost_function(nn_params, x, y, regulator=0):\n",
    "    m, n = x.shape\n",
    "    thetas = np.array(unroll_params(nn_params))\n",
    "    mapped_y = map_y(y)\n",
    "    hypothesis = feed_forward(thetas, x)[1][1]\n",
    "    term1 = -np.einsum('ij,ij->', mapped_y, np.log(hypothesis))\n",
    "    term2 = np.einsum('ij,ij->', (1 - mapped_y), np.log(1 - hypothesis))\n",
    "    cost = (term1 - term2) / m\n",
    "    cost = np.sum(cost)\n",
    "    reg_arr = np.asarray([np.sum(theta * theta) for theta in thetas])\n",
    "    reg_term = (regulator / 2 / m) * np.sum(reg_arr)\n",
    "    return cost + reg_term"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 156,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.2876291651613189"
      ]
     },
     "execution_count": 156,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cost_function(nn_params, X, y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 147,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.384487796242894"
      ]
     },
     "execution_count": 147,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cost_function(nn_params, X, y, 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Back Propagation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we will train the network using the back propagation algorithm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sigmoid_gradient(z):\n",
    "    return expit(z) * (1 - expit(z))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.19661193, 0.23500371, 0.25      , 0.23500371, 0.19661193])"
      ]
     },
     "execution_count": 81,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "z = np.linspace(-1, 1, 5)\n",
    "sigmoid_gradient(z)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gen_random_weights(lin, lout):\n",
    "    eps = np.sqrt(6) / np.sqrt(lin + lout)\n",
    "    w = np.random.rand(lout, lin + 1) * 2 * eps - eps\n",
    "    return w"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {},
   "outputs": [],
   "source": [
    "initial_theta1 = gen_random_weights(input_layer_size, hidden_layer_size)\n",
    "initial_theta2 = gen_random_weights(hidden_layer_size, num_labels)\n",
    "initial_nn_params = np.concatenate([initial_theta1.ravel(), initial_theta2.ravel()], axis=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 228,
   "metadata": {},
   "outputs": [],
   "source": [
    "def back_propagation(nn_params, x, y, regulator=0):\n",
    "    m, n = x.shape\n",
    "    \n",
    "    # Gradient arrays\n",
    "    theta1, theta2 = unroll_params(nn_params)\n",
    "    D1 = np.zeros_like(theta1)\n",
    "    D2 = np.zeros_like(theta2)\n",
    "    \n",
    "    # Feed forward\n",
    "    az_arr = feed_forward([theta1, theta2], x)\n",
    "    a3 = az_arr[1][1] # 10 x 5000\n",
    "    z3 = az_arr[1][0] # 10 x 5000\n",
    "    a2 = az_arr[0][1] # 25 x 5000\n",
    "    z2 = az_arr[0][0] # 25 x 5000\n",
    "    a1 = np.insert(x, 0, 1, axis=1) # 5000 x 401\n",
    "    \n",
    "    # Error at each node\n",
    "    mapped_y = map_y(y)\n",
    "    delta3 = a3 - mapped_y\n",
    "    delta2 = np.dot(theta2.T, delta3)[1:] * sigmoid_gradient(z2)\n",
    "    \n",
    "    # Unregularized gradient\n",
    "    D1 = (D1 + np.dot(delta2, a1)) / m # 25 x 401\n",
    "    D2[:, 1:] = (D2[:, 1:] + np.dot(delta3, a2.T)) / m # 10 x 25\n",
    "    \n",
    "    # Regularization\n",
    "    D1[:, 1:] = D1[:, 1:] + (regulator / m) * theta1[:, 1:]\n",
    "    D2[:, 1:] = D2[:, 1:] + (regulator / m) * theta2[:, 1:]\n",
    "    \n",
    "    return np.concatenate([D1.ravel(), D2.ravel()])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 256,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Training\n",
    "\n",
    "def train_network(initial_nn_params, x, y, regulator=0, iters=100):\n",
    "    result = opt.minimize(cost_function, initial_nn_params,\n",
    "                          args=(x, y, regulator),\n",
    "                          jac=back_propagation, method='TNC',\n",
    "                          options={'maxiter': iters})\n",
    "    return result.x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 254,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "95.46"
      ]
     },
     "execution_count": 254,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Train and predict\n",
    "\n",
    "nn_params = train_network(initial_nn_params, X, y)\n",
    "theta1, theta2 = unroll_params(nn_params)\n",
    "prediction = np.argmax(feed_forward([theta1, theta2], X)[1][1], axis=0)\n",
    "np.mean(prediction == y) * 100"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This network predicts the training data with an average accuracy of 95.46%. This could possibly be overfitting. One way to go around this would be to have a nonzero regulator."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 257,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "93.97999999999999"
      ]
     },
     "execution_count": 257,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nn_params = train_network(initial_nn_params, X, y, 1)\n",
    "theta1, theta2 = unroll_params(nn_params)\n",
    "prediction = np.argmax(feed_forward([theta1, theta2], X)[1][1], axis=0)\n",
    "np.mean(prediction == y) * 100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 258,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "99.02"
      ]
     },
     "execution_count": 258,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nn_params = train_network(initial_nn_params, X, y, 0.1, 150)\n",
    "theta1, theta2 = unroll_params(nn_params)\n",
    "prediction = np.argmax(feed_forward([theta1, theta2], X)[1][1], axis=0)\n",
    "np.mean(prediction == y) * 100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 259,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "93.76"
      ]
     },
     "execution_count": 259,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "nn_params = train_network(initial_nn_params, X, y, 10, 150)\n",
    "theta1, theta2 = unroll_params(nn_params)\n",
    "prediction = np.argmax(feed_forward([theta1, theta2], X)[1][1], axis=0)\n",
    "np.mean(prediction == y) * 100"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}