{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### <center>San Jose State University<br>Department of Applied Data Science</center>\n",
    "#  <center>DATA 220<br>Mathematical Methods for Data Analysis</center>\n",
    "### <center>Spring 2021<br>Instructor: Ron Mak</center>\n",
    "#  <center>Assignment #12<br>Polynomial Regression and Markov Chain Problem Set<br>SOLUTIONS</center>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import numpy as np\n",
    "from numpy.linalg import solve, inv\n",
    "import PolynomialRegression as pr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr>\n",
    "\n",
    "#### <strong>PROBLEM 1.</strong> Find the least-squares quadratic (degree 2) regression line through the following set of points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "xs = [11, 12, 13, 14, 15, 17, 18, 19, 20, 21,\n",
    "      24, 25, 27, 28, 31, 32, 33, 36, 37, 38,\n",
    "      40, 41, 44, 45, 46, 49, 50, 51, 54, 55]\n",
    "\n",
    "ys = [11.3, 15.1,  6.6, 12.9, 12.1, 18.1, 20.9, 17.6, 11.0, 24.6,\n",
    "      11.3, 18.4, 16.2, 19.5, 35.8, 37.1, 45.7, 34.8, 25.6, 26.7,\n",
    "      22.0, 26.0, 10.5, 18.6, 21.1, 11.9, 13.7, 13.7,  6.3,  1.8]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Print the quadratic equation and plot it among the points.\n",
    "\n",
    "#### <strong>SOLUTION:</strong>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([11, 12, 13, 14, 15, 17, 18, 19, 20, 21, 24, 25, 27, 28, 31, 32, 33,\n",
       "       36, 37, 38, 40, 41, 44, 45, 46, 49, 50, 51, 54, 55])"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "xs = np.array(xs)\n",
    "xs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([11.3, 15.1,  6.6, 12.9, 12.1, 18.1, 20.9, 17.6, 11. , 24.6, 11.3,\n",
       "       18.4, 16.2, 19.5, 35.8, 37.1, 45.7, 34.8, 25.6, 26.7, 22. , 26. ,\n",
       "       10.5, 18.6, 21.1, 11.9, 13.7, 13.7,  6.3,  1.8])"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ys = np.array(ys)\n",
    "ys"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[   1,   11,  121],\n",
       "       [   1,   12,  144],\n",
       "       [   1,   13,  169],\n",
       "       [   1,   14,  196],\n",
       "       [   1,   15,  225],\n",
       "       [   1,   17,  289],\n",
       "       [   1,   18,  324],\n",
       "       [   1,   19,  361],\n",
       "       [   1,   20,  400],\n",
       "       [   1,   21,  441],\n",
       "       [   1,   24,  576],\n",
       "       [   1,   25,  625],\n",
       "       [   1,   27,  729],\n",
       "       [   1,   28,  784],\n",
       "       [   1,   31,  961],\n",
       "       [   1,   32, 1024],\n",
       "       [   1,   33, 1089],\n",
       "       [   1,   36, 1296],\n",
       "       [   1,   37, 1369],\n",
       "       [   1,   38, 1444],\n",
       "       [   1,   40, 1600],\n",
       "       [   1,   41, 1681],\n",
       "       [   1,   44, 1936],\n",
       "       [   1,   45, 2025],\n",
       "       [   1,   46, 2116],\n",
       "       [   1,   49, 2401],\n",
       "       [   1,   50, 2500],\n",
       "       [   1,   51, 2601],\n",
       "       [   1,   54, 2916],\n",
       "       [   1,   55, 3025]])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = np.array([[1, xs[i], xs[i]*xs[i]] for i in range(len(xs))]) \n",
    "A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([11.3, 15.1,  6.6, 12.9, 12.1, 18.1, 20.9, 17.6, 11. , 24.6, 11.3,\n",
       "       18.4, 16.2, 19.5, 35.8, 37.1, 45.7, 34.8, 25.6, 26.7, 22. , 26. ,\n",
       "       10.5, 18.6, 21.1, 11.9, 13.7, 13.7,  6.3,  1.8])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "b = ys\n",
    "b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[   1,    1,    1,    1,    1,    1,    1,    1,    1,    1,    1,\n",
       "           1,    1,    1,    1,    1,    1,    1,    1,    1,    1,    1,\n",
       "           1,    1,    1,    1,    1,    1,    1,    1],\n",
       "       [  11,   12,   13,   14,   15,   17,   18,   19,   20,   21,   24,\n",
       "          25,   27,   28,   31,   32,   33,   36,   37,   38,   40,   41,\n",
       "          44,   45,   46,   49,   50,   51,   54,   55],\n",
       "       [ 121,  144,  169,  196,  225,  289,  324,  361,  400,  441,  576,\n",
       "         625,  729,  784,  961, 1024, 1089, 1296, 1369, 1444, 1600, 1681,\n",
       "        1936, 2025, 2116, 2401, 2500, 2601, 2916, 3025]])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "AT = A.T\n",
    "AT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[      30,      946,    35368],\n",
       "       [     946,    35368,  1473706],\n",
       "       [   35368,  1473706, 65637328]])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ATA = AT@A\n",
    "ATA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.21793903e+00, -8.11479370e-02,  1.16568021e-03],\n",
       "       [-8.11479370e-02,  5.84526527e-03, -8.75136517e-05],\n",
       "       [ 1.16568021e-03, -8.75136517e-05,  1.35199922e-06]])"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ATAinv = inv(ATA)\n",
    "ATAinv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([5.669000e+02, 1.786420e+04, 6.333008e+05])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ATb = AT@b\n",
    "ATb"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Coefficients computed with the matrix normal equations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-20.96713356,   2.99575668,  -0.04631508])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "xhat = ATAinv@ATb\n",
    "xhat"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Coefficients computed with the system of normal equations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-20.96713356,   2.99575668,  -0.04631508])"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "_, _, a = pr.polynomial_regression(2, xs, ys)\n",
    "a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAFXCAYAAABN1VJsAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAAArzUlEQVR4nO3de5zN1f7H8dcyxNBlUhSDdHEZdxEVdSSikiSVonISXRRdfmp0Ot0cmcxRSVd1lJNTThdJFNLoJsplOJSmSMmQSFPRKGbW7481dmiGmdn7u7/fvff7+XjMY+zv7L2/ny9mf75rrc9ay1hrERERAajgdwAiIhIcSgoiIhKipCAiIiFKCiIiEqKkICIiIUoKIiISUtHvAPZkjJllre1expepplZEpOxMcQeD1lI40u8AREQSWdCSgoiI+EhJQUREQpQUREQkRElBRERClBRERCRESUFEREKUFEREJERJQUREQpQUREQkRElBRERCArX2kUi8mZadS+bsHDbk5VM7JZnh3RrRq3Wq32GJlEhJQcQj07JzGTF1Bfk7CwDIzctnxNQVAEoMEljqPhLxSObsnFBC2C1/ZwGZs3N8ikjkwJQURDyyIS+/TMdFgkBJQcQjtVOSy3RcJAiUFEQ8MrxbI5IrJe11LLlSEsO7NfIpIpED00CziEd2Dyar+khiibE2OLtZGmMWW2vblvFlwbkAEZHYERPbcYqIiI+UFEREJERJQUREQpQUREQkRElBRERClBRERCRESUFEREKUFEREJERJQUREQpQUREQkRElBRERClBRERCRESUFEREKUFEREJERJQUREQpQUREQkRElBRERClBRERCRESUFEREKUFEREJERJQUREQpQUREQkRElBRERCIpYUjDFJxphsY8yMosfVjTFvG2O+LPp+eKTOJSIi3ohkS2EYsGqPx+nAO9baBsA7RY9FRCTAIpIUjDF1gHOBZ/Y4fD4wqejPk4BekTiXiIh4J1IthYeB24DCPY4dZa3dCFD0vWZxLzTGDDbGLDbGLAaOjFA8IiJSDmEnBWNMD+B7a+2S8rzeWjvBWtvWWtsW2BJuPCIiUn4VI/AeHYCexphzgCrAocaYycAmY0wta+1GY0wt4PsInEtERDwUdkvBWjvCWlvHWlsf6AtkWWv7A9OBK4uediXwerjnEhERb3k5TyED6GqM+RLoWvRYREQCzFhr/Y4hxBizuGhsoSyCcwEiIrHDFHdQM5pFRCRESUFEREKUFEREJERJQUREQpQUREQkRElBRERClBRERCRESUFEREIisfaRSMKblp1L5uwcNuTlUzslmeHdGtGrdarfYYmUmZKCSJimZecyYuoK8ncWAJCbl8+IqSsAlBgk5qj7SCRMmbNzQglht/ydBWTOzvEpIpHyU1IQCdOGvPwyHRcJMiUFkTDVTkku03GRIFNSEAnT8G6NSK6UtNex5EpJDO/WyKeIRMpPA80iYdo9mKzqI4kH2k9BpIxUfipxotj9FNRSECkDlZ9KvNOYgkgZqPxU4p2SgkgZqPxU4p2SgkgZqPxU4p2SgkgZ+FF+Oi07lw4ZWRybPpMOGVlMy8717FwiGmgWKYNol59qYFuiTSWpIgHWISOL3GLGK1JTkpmf3tmHiCSOFFuSqu4jkQDTwLZEm5KCSIBpYFuiTUlBJMC0rpJEmwaaRQJM6ypJtGmgWUQkMWmgWURE9k9JQUREQpQUREQkRElBRERClBRERCRESUFEREKUFEREJERJQUREQpQUREQkRElBRERClBRERCRESUFEREKUFEREJERJQUREQpQUREQkRElBRERClBRERCRESUFERELCTgrGmCrGmE+MMcuNMZ8aY+4tOl7dGPO2MebLou+Hhx+uiIh4KRIthd+AztbalkAroLsx5mQgHXjHWtsAeKfosYiIBFjYScE624oeVir6ssD5wKSi45OAXuGeS0REvBWRMQVjTJIxZhnwPfC2tfZj4Chr7UaAou81I3EuERHxTkSSgrW2wFrbCqgDtDPGNCvta40xg40xi40xi4EjIxGPiIiUT0Srj6y1ecC7QHdgkzGmFkDR9+9LeM0Ea21ba21bYEsk4xERkbKJRPVRDWNMStGfk4EuwOfAdODKoqddCbwe7rlERMRbFSPwHrWAScaYJFySeclaO8MYswB4yRgzEFgHXBSBc4mIiIeMtdbvGEKMMYuLupHKIjgXICISO0xxBzWjWUREQpQUREQkRElBRERClBRERCRESUFEREKUFEREJERJQUREQiIxeU1EdisshB07oEIFqFgRkpLAFFsOLhJISgoipbV1K/zvf7B6NaxZA199Bd99B5s2webNsH07/Pbbn1938MFw+OHuq3ZtqFcP6taFBg2gWTP3/aCDon89IsVQUhApzs6dsGQJvP8+fPQRZGfDunV//LxSJahfH2rVgubNoUYNOOQQSE6GKlXcc3btcu/z88+Ql+eSyoYNsHQpfL/H+pAVK7rkcPLJcMop0LEjHHdcNK9WJETLXIjs9s03MHOm+3r3Xfj1V3e8YUNo0wZat4YWLaBRI3enn5RU/nPl58MXX8DKle5r8WL4+GP45Rf38+OOg7POgrPPhq5dXbIphWnZuWTOzmFDXj61U5IZ3q0RvVqnlj9OiWfF9msqKUhiW7cO/vtfePFF1xoAOP5492HcqROcdhrUjNL+UAUFsGqVS0hz5sC8ebBtG1SrBuecAxdfDD16/NES2ce07FxGTF1B/s6C0LHkSkmM7t1ciUGKo6QgAriB4KlT4emn3QcwQLt2cNFFcN55rmUQhMHhnTtdfK++Cq+95rqcDj8c+vaFq66Ctnv/qnTIyCI3L/9Pb5Oaksz89M5RClpiiJKCJLh162D8eJg40fXvH3cc/PWvcOmlrnUQZAUFkJUFzz3nEtqOHS6R3XijS2aVK3Ns+sxifxkMsDbj3CgHLDFASUES1NKlkJkJL7/sHvfuDddcA2ec4UpHY81PP8Hzz8Ojj0JODhx1FNx8M13zm/Dljj9fT1laChqTSChKCpJgliyBe++FN96AQw+FQYNg6FBXEhoPCgth7lwYOxbmzOH3Qw7jmZbn8GSb8/m5ysFA2cYUNCaRcJQUJEF89hmMGAHTp7s++FtvhRtugMMO8zuyiNn3jn5U3R10evUZeO01fk4+hMfbX8icMy5i6HktS/2BrjGJhKNNdiTOffed6xZq3twN0I4cCV9/DX/7W9wlhBFTV5Cbl48FcvPyue7zCkz7+3hYupRDO3Uk/d3nyHpqEL3+N9e1KEphQzEJYX/HJT4pKUjs27kTHnzQzQyeOBGGDHEzju+803UbxZnM2Tl7dfEA5O8sIHN2jptL8eabbtJdnTowYICbEPfxxwd839opxc+FKOm4xCclBYlt773nPghvvRVOP911HT3yCBx5pN+ReaZUd/SnnQYLFsCkSa7q6uSTYfBg+PHHEt93eLdGJFfae0JecqUkhndrFJG4JTYoKUhs+ukn9yHXqZOb4DVtGsyY4VoLca7Ud/QVKsAVV7iZ07feCv/6F6Slucl6xYwl9mqdyujezUlNScbgxhI0yJx4NNAssWfmTDd2sHEj3HKLqzCqWtXvqKKm3FVCS5e6RLpkCVxwATz5ZPRma0sQaaBZomdadi4dMrI4Nn0mHTKymJadG/6bbt/uPtR69ICUFNc9kpmZUAkBwrijP/FEWLgQxoxxibVZMzdTWmQPailIxHlS775oEfTr55atHj4c7rsPKleOUMQJaOVK17WUne2WzHjkEbfGkiQStRQkOvZbHVNW1rrJWaee6pZ2yMqCBx5QQghXs2au1TBiBDz7LJx0EqxY4XdUEgBKChJxEat3//FH6NUL/u//oGdPt8FNp05hxydFDjoI7r/frcj6448uMTzzjN9Ric+UFCTiIlLvnp3t+sDfegvGjYNXXnHjCBJ5XbrA8uWupHfQILj6atcqk4SkpCARF3a9+4svQocObueyDz906xUFYSnreFazpkvAd9zhSldPO23vneYkYSgpSMSVuzqmoABuuw0uu8ztFbBkiVseWqIjKQlGjXJzPr74wnUnLVjgd1QSZao+kmDYts0lgzfegOuvh4ce0mb2flq1ym04tH69azn06+d3RBJ5qj6SgNqwwfVnz5zp9gh47DElBL+lpbn1kk4+Gfr3h7vuKnYWtMQfJQXx14oV0L49fPmlayUMGeJ3RLLbEUe4yqSrrnIrzl51lVt8UOJaRb8DkAT24Yeui6JqVfjgA2jVyu+IZF8HHeTKVOvVg3vucUuLvPIKHHyw35GJR9RSEH9Mnw5du7qtJD/6SAkhyIyBu+92yWHuXLeN6ZYtfkclHlFSkOj797/dPsktWrjWwjHH+B2RlMbAga4yaeVK+Mtf3FiQxB0lBYmuJ5+EK690d5vvvBPX+x7EpR493HyGdeugY0f46iu/I5IIU0mqRM9DD7mlrnv0gJdfhipV/I4ooe27z/Pwbo1Kv2DhokXQvbv7N3z33YTYxyIOqSRVfJSR4RJCnz7w6qtKCD4rbp/nEVNXlH6J85NOcslg5063HtUXX3gYrUSTkoJ4LzPTrcZ56aVuCQvNQfBdRFaybd4c5s1zy5F06gQ55VgFVwJHSUG89eCDbumKSy5xA8wVVQUdBBFbybZpU5cYCgtdYli9OvzgxFdKCuKd8ePd3sB9+sDkyUoIARKRlWx3a9LE7XOxaxeceSZ8802Y0YmflBTEG88+61Y37dULXnhBCSFgwl7Jdl9NmsDbb8PPP7vEUMpyVU+2bZWwqPpI/iSsqhRwA8kXX+w+HN54Q7ukBVTY/87F+fhjtz9D3bpulvoRR+z3/BHftlXKotjqIyUF2UvYv6hz5riS05NOcn/Wvr+J5733oFs3N0t97twSl8TokJFFbjFjGKkpycxP7+xxkIJKUqU0wqpKWbTIzVRu0sSteKqEkJj+8heYMsX9f+jTB37/vdinRWywWyJKSUH2Uu5f1NWr4dxzoUYNN+NVW2cmtl69YMIEmD0bBgxw1Un7iOhgt0RM2EnBGFPXGDPPGLPKGPOpMWZY0fHqxpi3jTFfFn0/PPxwxWvl+kXdtMnNbi0sdB8CtWp5FJ3ElIED4f773dyUESP+9OOID3ZLRESipbALuNVamwacDAwxxjQB0oF3rLUNgHeKHkvAlfkX9ddf3RjChg2uy6hhwyhEKTEjPR2uuw7GjIEnntjrR+XetlU8FfGBZmPM68CjRV+drLUbjTG1gHettfu9BdBAczCUuiqlsND1GU+b5r569ox2qBILdu1y3UlvvQWvv+5uIiQIvK8+MsbUB94HmgHrrLUpe/zsR2vtfruQlBRizPDh8M9/wsMPw7BhfkcjQbZtm5vxvGqVWy69dWu/IxKvk4Ix5mDgPWCUtXaqMSavNEnBGDMYGFz08Ehrbf0ynlpJwQ8TJsA117jtM8ePdxuxxABPavOldL77Dtq1c3s9f/KJxp78511SMMZUAmYAs621DxYdy0HdR/Hp3Xfdrmldu7od1GJktrImSwXAsmXQoQM0a+b+HyWr0shH3sxTMMYY4F/Aqt0Joch04MqiP18JvB7uuSQA1q514wgNGrha9BhJCBChlUElPK1awX/+41oKAwe6VoMESiSqjzoAlwOdjTHLir7OATKArsaYL4GuRY8llv3yixtMLix0LYRDD/U7ojLRZKmA6NULRo92paoPPOB3NLKPsG/zrLUfUkIzBDgz3PeXgCgshCuucAOFs2fDCSf4HVGZ1U5JLnZZBU2W8sHtt8Py5XDHHa710L273xFJEc1oltIZPdqVnY4d6xa6i0GaLOWfP62GumwDPPOM26jn0ku1D0OAaEE8ObBZs+Ccc+Cyy+D5532tNAq3ekjVR9G33wH+lN+hbVtXibRwYYmL54kntEqqlMPatdCmjVsKecECqFrVt1BUPRSbDrga6ty5blXVSy5xg9AxUt4cB7RKquztgBuc5Oe7VU+thalTfU0IoOqhWHXAAf4uXWDkSDfw/PjjUYxMiqOkkKB233Xn5uVjgdy8fEZMXbF3YrjpJldX/vzzcPzxPkX6B1UPxaZSLbKYnu5W2b35ZrdRj/hGSSFBHfCue/JkN2s5PR169AjEtolaajk2lWqAv0IF+Pe/ITUVLroItmyJcpSym5JCgtrvXfdnn7klLE47DUaOLF2rIgpUPRSbSr0aavXq8PLLbin2AQM0sc0nGmhOUCUN/h1X1ZD13+GwebPrOqpdO1DbJqp6KAGMHw9Dh7ry51tu8TuaeKbqI/lDSZU8b336PPVfe8FNUDvrLACOTZ9Z7F+yAdZmnBudgCVxWOsKHGbMgPnz3SJ64gVVH8kfimvSTzrkG5cQbr89lBBAffkSZcbAxIlQuzb07Qt5eX5HlFDUUhDn66/dcgONGrn17itVCv1I8wPEFwsWuHGtiy6CF17Q/IXIU0tBSrBrl5utXFjoasX3SAigbRPFJ6ecAvfe61bjnTzZ72gShloK4n7x7rnHzSa97DK/oxH5Q0EBnHGGK3pYtgyOO87viOKJBpqlGB9/7DY96ds37u7GVKkUJ775Blq2hLQ0eP/9P7VkpdzUfST72LYN+vd3E4YefdTvaCIqKHMrJAKOOQaefNItmDdqlN/RxD0lhUR2yy2wZo2bSZqS4nc0EaV1kuJM377Qrx/84x+waJHf0cQ1JYVENWMGPP003HYb/OUvfkcTcVonKQ49+igcfTRcfrlbrFE8oaSQiH74AQYNghYt3CBzHNLcijiUkgLPPQc5OTBihN/RxC0lhUR0441uwbFJk6ByZb+j8YTWSYpTXbq4/7/jxkFWlt/RxCVVHwWM5xUzr7ziJgONHAl33hm59w0gVR/FqV9/hdat4bffYMUKOOQQvyOKVSpJDTrPZw5//z00bQr167vZohUrhv+eIn746CPo2BGuvVYb85SfSlKDzvOKmRtugJ9/dt1GSggSy0491W3I88QT6kaKMCWFAPG0YmbqVLdW/d13Q5Mm4b+fiN9GjoQGDWDgQDfnRiJCSSFAPKuY2boVrr/e9cMOHx7ee4kERdWqbjXVb75xOwRKRCgpBIhnFTM33+zKUCdO1BIBEl86dnTVSI8/7vZekLBpoDlgIl4xM2sWnH22qzQaOTJygYoExbZtroCiWjXIzo7bMmsPqPoo4Wzf7n5ZqlbVL4vEt7fegnPOgbvuitsJmR4oNimoBCWe3XWX62/94AMlBIlvZ5/t1kYaPRouuohpOw8vtsWtuSsHppZCvFqyxO1tO2iQW2FSJN5t2QJpaWytVY+OPUfy664/PhqSKyVxYZtUXl2Sqx0E/6B5Cglj1y6XDGrWhIwMv6MRiY4jj4SxY6m+Yim9Fr+514/ydxbw4sffauXcUlBSiEePPOLGEMaPj7slsUX26/LL+aheC9LffY4a237c60cFJfSKaOXcvSkpxLBp2bl0yMji2PSZdMjIchvIfPutG0s491y48EK/QxSJLmMY3+cWKu/6jTuzntnrR0mm2N4SrZy7DyWFGFXSzmIbBlwDhYWulVDCL4FIPLukfxeePvUSzl/1HqetXQq4sYNL29fVyrmloKQQo4pbJ6n95wupnfWWm5Nw7LE+RSbir16tU6k35l7WHZHKyDlPcGw1N5j8j17NGd27OakpyRggNSU5kQeZS6Tqoxh1bPrMvS68ys4dzPnXEH5PqsQJG9fAQQf5FptIILz9Npx1Ftx3H/z9735HE0SqPoon+/aDXrfwFer9tIlxvW9SQhAB6NrV7R1y//2wdq3f0cQMJYUYtec6SfV+3Mi1H7/KG007cdi5Z/158FkkUT34ICQlwU03+R1JzNCM5hi1ux80c3YOd78ygYKkinx2y517Tc7ZPfi85/NF4skBZyjXqeOWi7/tNpgxA3r08C/YGKExhVg3Ywacdx5kZtJh14nkFlNznZqSzPz0zj4EJ+KdUu9U+Pvv0KqV277z00+hSpXoBxtMGlOIOzt2uGZx48YwdKi3m/SIBEypdyo86CBXov3VV647SfZL3Uc+icjCXGPHwpo1rsrioIOonZJcbEshFifnaOEyOZAy3QSdeSb07g2jRsEVV7huJSmWWgo+KGniWZkGhXNzXVVF797QpQvg4SY9URaRvx+Je2XeqXDsWDex8/bbPYwq9ikp+KDUzd79SU+HggL45z9Dh3q1To365Jxil9oIU0T+fiTulfkmqH59tx3tCy/Ahx96H2CMUveRD8Lu+1+wACZPhjvu+NPM5V6tU6PWzbLvQF+kqp00NiKlsWcFXqm7GdPT4bnnYOhQWLTIlavKXpQUfBBW339hIQwbBrVqwYgRHkRXevu7ow8nKcTT2Ih4q8w3QVWrwpgxcOmlMGkSXHWVd8HFKHUf+SCsvv/nn3d3OA88AAcf7FGEpePVHX28jI0kKi+6FCPqkkvglFPgb3+DX37xO5rAiUhSMMZMNMZ8b4xZucex6saYt40xXxZ9PzwS54oH5e77377dtQ7atXNbD/qszAN9peTH2IhERkwUCRgDDz0E332nTaiKEZHJa8aY04FtwL+ttc2Kjo0BtlprM4wx6cDh1tr9Dvtr8toB3Hsv3HMPzJ8Pp57qdzSlnzwkCaNDRlbsTKDs3x9eeQVycuCYY/yOxg/eTV6z1r4PbN3n8PnApKI/TwJ6ReJcCWvDBtcXetFFgUgIoDt6+bOYKhIYPRoqVFCJ6j68HGg+ylq7EcBau9EYU9PDc8W/O+90ey8HrLkbzWonCb6YKhKoWxf+7/9g5Ei3MsDJJ/sdUSD4PtBsjBlsjFlsjFkMHOl3PIG0bNkfZXTHHed3NCIlirkigeHD4aijXHII0DpwfvIyKWwyxtQCKPr+fXFPstZOsNa2LRpL2OJhPLHJWvcftnp1Vy0hEmAx16V4yCFurG7+fJg2ze9oAiFiq6QaY+oDM/YYaM4EfthjoLm6tfa2A7yHBpr3NWsWnH02PPywm58gIpG1axe0aOG+f/opVKrkd0TRUuxAc6Sqj14EOuG6fzYBdwPTgJeAesA64CJr7b6D0fu+j5LCngoK4MQTYds2WLVKO6qJeOWNN6BnT3j0URgyxO9oosW7pBApSgr7mDQJBgyAKVPchBsR8Ya10LkzrFzpVh4+9FC/I4oG7acQU/Lz3Wbjbdu6MlQR8Y4xkJkJW7bstchkIlJSCKrx4+Hbb93chAr6ZxLx3O4bsAcfhE2b/I7GN+o+CqIff3Slp6eeCjNn+h2NSOL44gto0gSuvdaNL8Q3dR/FjDFj4Kef3IxLEYmehg3h6qvhqadg9Wq/o/GFWgpBs3EjHH+821Ft8mS/oxFJPBs3wgknuGqkF1/0OxovqaUQE0aOhJ073YQaEYm+WrXcshdTpkB2tt/RRJ1aCkGyZg00bgyDB8Njj/kdjUjiystzuxrG97ieWgqBd9ddboLanXf6HYlIYktJcaunvvmmWwIjgSgpBMWKFa7/cuhQ13wVEX/deKNbLO9vf0uoxfKUFILi7rvd4lzDh/sdiYgAVKvmEsJ778HcuX5HEzVKCkGwZAm89hrceqtbDVVEgmHwYKhXD+64I2FaCxpoDoJzzoGPP4a1a4tdc2Vadi6Zs3PYkJdP7ZRkhndrFNyliEXizcSJMHCgu3Hr1cvvaCJJC+IF0vz50LEjPPAA3PbnlcW1D7KIz3btgrQ01520dGk8LTujpBBIZ5zhlsVes8b9p9tHTG2ELhKvJk+Gyy/nkzFPcnNBg3hptaskNXDefdd9jRhRbEKAGNsIXSReXXopv9Q/gepj7mfj1m1YIDcvnxFTVzAtO9fv6CJKScEv1rqKo1q13GBWCUra8DyQG6GLxKukJDJP6csJW9Zx7ucfhg7n7ywgc3aOj4FFnpKCX+bNg/ffd62E5JI/4GNuI3SRODW5bjtyjqzHTfNfpELhH2N88dZqV1Lwg7Vwzz1QuzYMGrTfp8bcRugicarW4dV4qGM/jt+6np6r3g8dj7dWe0W/A0hIWVnwwQduI50qVQ749F6tU5UERHw2vFsj7ti2g1Xz63PjR/9letrpVK58UNy12pUUom13KyE11a3bLiIxYfeN2eS1VzDqhfu4fN0ntL79uri7YVNJarRlZcGZZ7pdnYYM8TsaESmrwkJo0cJ9X7ECkpIO/JpgUklqINx3n6s4GjjQ70hEpDwqVIC//93NL3rlFb+jiTi1FKLpgw/g9NPh4Ydh2DC/oxGR8ioogObNXYL43/9idZazWgq+GzkSatY8YMWRiARcUpJrLXz6KUyd6nc0EaWkEC0LF8Lbb7ulsatW9TsaEQnXxRdDw4YwalRcraCqpBAtI0fCEUfAtdf6HYmIREJSkpt8umyZ26EtTigpRMPSpe4/zS23wMEH+x2NiERKv35wzDHupi9OWgtKCtFw//1w2GEqQRWJN5UqQXq62w8lK8vvaCJCScFrq1a5gagbbnCJQUTiy4ABrsx81Ci/I4kIJQWvZWS4Be9UgioSn6pUcQUk8+a5TbNinJKCl77+Gv7zH7c0do0afkcjIl4ZPNgVkmRk+B1J2JQUvDRmjJvUcuutfkciIl6qVs31BsyY4Za+iGFKCl757ju34feVV0KdOn5HIyJeGzLEVRc+8IDfkYRFScErDz8MO3fCbbf5HYmIREP16nDNNTBlCqxd63c05aak4IWffoInnoALL4QGDfyORkSi5eabXZfxP//pdyTlpqTghSefhJ9/httv9zsSEYmm1FTXZTxxImza5Hc05aKkEGk7driuoy5doE0bv6MRkWi77Tb47TcYN87vSMpFSSHSnn/eDTKnp/sdiYj4oUED13X8+OPwyy9+R1Nm2k8hkgoKoHFjN3N50SIwxS5XLiJxZlp2Lpmzc9iQl0/tlGRG1cmn0+U9YOxYt+ZZMGk/Bc9NmwarV7uxBCUEkYQwLTuXEVNXkJuXjwVy8/K5LieJzW1PgYceclWIMURJIVKshcxMOP546N3b72hEJEoyZ+eQv7Ngr2P5OwvIaHYerF/vSlRjiJJCpMyf71ZKvOWWWN7IW0TKaENefrHHp9ZsDs2auZUNAtRNfyBKCpGSmenWPhkwwO9IRCSKaqckF3/88KquEmnlSpg1K8pRlZ+SQiTk5MD06W6au7baFEkow7s1IrnS3r0DyZWSGN6tEfTt65a5ycz0KbqyU1KIhLFj3fK52kRHJOH0ap3K6N7NSU1JxgCpKcmM7t2cXq1T3SY8w4a5ZbWzs/0OtVRUkhquTZvcdnwDBriZzCIie/rpJ6hbF3r2hMmT/Y5mTypJ9cTjj8Pvvwe5FllE/HTYYXD11fDf/7pqpIDzPCkYY7obY3KMMauNMfE1zTc/3yWF886Dhg39jkZEgmrYMFeB9MgjfkdyQJ4mBWNMEvAYcDbQBLjUGNPEy3NG1eTJsGWLWgkisn/HHAN9+sBTT7nFMgPM65ZCO2C1tfYra+3vwBTgfI/PGR2FhW624oknwumn+x2NiATdrbe6hPCvf/kdyX55nRRSgW/3eLy+6Fjsmz0bVq1yrQQtaSEiB3LSSdCxo+tCKig48PN9UtHj9y/u03KvaiFjzGBgcNHDI8t6gu7du7Nly5ZyhBamL75w5WYPPuhaDKWwefNmatSo4XFg/tC1xa54vr7AXVteHnz9NTRqBCkpYb1VuNe2ZMmSWdba7n/6gbXWsy/gFGD2Ho9HACP28/zF5ThP9C1fbi1Ym5FRppe1adPGo4D8p2uLXfF8fYG7tl27rK1f39qOHcN+qwhcW7GfqV53Hy0CGhhjjjXGHAT0BaZ7fE7vPfywm7k8aJDfkYhILElKgqFD4cMPYfFiv6MplqdJwVq7C7gBmA2sAl6y1n7q5Tk99/338MILbsu96tX9jkZEYs3AgXDIIe7mMoA8n6dgrX3TWtvQWnu8tXaU1+fz3FNPua32hg4t80sHDx584CfFKF1b7Irn6wvktR16KFx1lZvMlptb7rfx6tq0zEVZ/P67qzdu1QreeitqpxWROPPVV3DCCW7b3vvv9ysKLXMRtpdecvsvDxvmdyQiEsuOOw7OPx8mTHArIwSIkkJpWQvjxrk9mM86a79Pveqqq6hZsybNmjULHdu6dStdu3alQYMGdO3alR9//NHriD3x7bffcsYZZ5CWlkbTpk0ZN24cED/Xt2PHDtq1a0fLli1p2rQpd999NxA/1wdQUFBA69at6dGjBxA/11a/fn2aN29Oq1ataNvWdTgE+tqGDoUffoAXXyzV0/Py8ujTpw+NGzcmLS2NBQsWeHJ9Sgql9dFHrlpg6FCosP+/tgEDBjBrn001MjIyOPPMM/nyyy8588wzycjI8DJaz1SsWJGxY8eyatUqFi5cyGOPPcZnn30WN9dXuXJlsrKyWL58OcuWLWPWrFksXLgwbq4PYNy4caSlpYUex9O1zZs3j2XLlrG4qLIn0NfWqZPbme2RR0q1M9uwYcPo3r07n3/+OcuXLyctLc2b6yupVtWPL4I8T+Hii61NSbF227ZSPX3t2rW2adOmoccNGza0GzZssNZau2HDBtuwYUNPwoy2nj172jlz5sTl9W3fvt22bt3aLly4MG6u79tvv7WdO3e277zzjj333HOttfHzf/OYY46xmzdv3utY4K9twgQ35+m99/b7tJ9++snWr1/fFhYW7nU8zOvzZZ5CfMjNhVdfdRUD1aqV6y02bdpErVq1AKhVqxbff/99JCP0xddff012djbt27ePq+srKCigVatW1KxZk65du8bV9d10002MGTOGCnu0duPl2owxnHXWWbRp04YJEyYAMXBt/fq50vbx4/f7tK+++ooaNWrw17/+ldatW3P11Vezfft2T65PSaE0nnzSLYCnndVCtm3bxoUXXsjDDz/MoYce6nc4EZWUlMSyZctYv349n3zyCStXrvQ7pIiYMWMGNWvWpE2bNn6H4on58+ezdOlS3nrrLR577DHef/99v0M6sN2TYF97DdatK/Fpu3btYunSpVx33XVkZ2dTrVo1z7rClBQO5LffXIXAuee6ioFyOuqoo9i4cSMAGzdupGbNmpGKMOp27tzJhRdeSL9+/ejduzcQX9e3W0pKCp06dWLWrFlxcX3z589n+vTp1K9fn759+5KVlUX//v3j4toAateuDUDNmjW54IIL+OSTT2Lj2q6/3o0pPPFEiU+pU6cOderUoX379gD06dOHpUuXenJ9SgoH8vLLbhbzjTeG9TY9e/Zk0qRJAEyaNInzz4/NFcSttQwcOJC0tDRu2WMfiXi5vs2bN5OXlwdAfn4+c+fOpXHjxnFxfaNHj2b9+vV8/fXXTJkyhc6dOzN58uS4uLbt27fzyy+/hP48Z84cmjVrFhvXVq8e9Oq13/LUo48+mrp165KTkwPAO++8Q5MmTby5vpIGG/z4IogDze3aWduokbUFBaV+Sd++fe3RRx9tK1asaFNTU+0zzzxjt2zZYjt37mxPOOEE27lzZ/vDDz94GLR3PvjgAwvY5s2b25YtW9qWLVvamTNnxs31LV++3LZq1co2b97cNm3a1N57773WWhs317fbvHnzQgPN8XBta9assS1atLAtWrSwTZo0sf/4xz+stTF0bVlZbsD52WdLfEp2drZt06aNbd68uT3//PPt1q1bw72+Yj9TNaN5fz75BNq3d4NAN9zg2WlEJMFZ68pTk5Nh0aJo7dGiGc1l9uijcPDBcMUVfkciIvHMGFfIsmSJuxn1kZJCSTZvdgtWXXmlW8BKRMRLl1/uVk999FFfw1BSKMnEiW4BvOuu8zsSEUkEhxzibkJfeskVt/hESaE4BQWuPKxTJ2ja1O9oRCRRDBnibkafeca3EJQUivPmm/DNN5qsJiLR1bgxdOnibkp37fIlBCWF4jz2GNSu7Za2FRGJpuuvh/XrYeZMX06vpLCv1ath9mwYPBgqVfI7GhFJNOedB6mp+53h7CUlhX098QRUrOiSgohItO3+/Jk9G9asifrplRT2lJ8Pzz3nppwXrTwoIhJ1V18NSUluMc4oU1LY0yuvwNatKkMVEX/Vrg0XXADPPgs7dkT11EoKe3riCWjYEM44w+9IRCTRXXed267z5Zejelolhd2WL4cFC+Daa6O17oiISMnOOAMaNYLHH4/qaZUUdnvqKahSxc0oFBHxmzHuJnXhQnfTGiVKCgC//ALPPw+XXOK2xhMRCYIrroDKld1Na5QoKQC88AJs2+aysohIUFSvDhdfDJMnu8+oKFBSsNZl4ZYt3d4JIiJBcu21rjdjypSonE5JYfFiyM6Ga67RALOIBM8pp7gNeKLUhaSkMGECVK0Kl13mdyQiIn9mjLtpXbzYbcLjscROCj//DC++CJdeCocd5nc0IiLF69/fbdUZhdZCYieFF16A7du1zpGIBFtKCvTt6z6zfv7Z01MlblLYc4D5pJP8jkZEZP8GD3Y3sR4POCduUliyBJYtc3/RGmAWkaBr394NOD/9tKenSdyk8NRTboC5Xz+/IxEROTBjYNAgN+C8bJlnp0ncpJCUBAMGaIBZRGJH//5uOR4PWwvGWuvZm5eVMWaxtbZtGV9W/guwVl1HIhJbLr8cpk+HjRtdb0f5Ffvhl7gtBVBCEJHYM2iQq0B66SVP3j6xWwoiIrHGWkhLgyOOgPnzw3kntRRERGKeMW67zo8+glWrIv72FSP+jiIi4q0BA6BNG2jcOOJvre4jEZHEpO4jERHZPyUFEREJUVIQEZEQJQUREQlRUhARkZCwkoIx5iJjzKfGmEJjTNt9fjbCGLPaGJNjjOkWXpgiIhIN4c5TWAn0BvbaDsgY0wToCzQFagNzjTENrbUFYZ5PREQ8FFZLwVq7ylqbU8yPzgemWGt/s9auBVYD7cI5l4iIeM+rMYVU4Ns9Hq8vOiYiIgF2wO4jY8xc4OhifvQ3a+3rJb2smGPFzjw2xgwGdm+SfOSB4hEREe8cMClYa7uU433XA3X3eFwH2FDC+08AJoBb5qIc5xIRkQjxqvtoOtDXGFPZGHMs0AD4xKNziYhIhIRVfWSMuQAYD9QAZhpjlllru1lrPzXGvAR8BuwChpSy8mhLecIox2tERKQYgVolVURE/KUZzSIiEqKkICIiIUoKIiISoqQgIiIhSgoiIhKipCAiIiFKCiIiEqKkICIiIf8PFOmdzVXWojIAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "The quadratic equation is y = -20.967 + 2.996x - 0.046x^2\n"
     ]
    }
   ],
   "source": [
    "pr.plot_polynomial(xs, ys, xhat, degree=2, x_min=5, x_max=60)\n",
    "\n",
    "print()\n",
    "print(f'The quadratic equation is y = {xhat[0]:.3f}', end='')\n",
    "if xhat[1] >= 0: \n",
    "    print(f' + {xhat[1]:.3f}x', end ='')\n",
    "else:\n",
    "    print(f' - {-xhat[1]:.3f}x', end ='')\n",
    "\n",
    "if xhat[2] >= 0: \n",
    "    print(f' + {xhat[2]:.3f}x^2', end ='')\n",
    "else:\n",
    "    print(f' - {-xhat[2]:.3f}x^2')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr>\n",
    "\n",
    "#### <strong>PROBLEM 2.</strong> If you have <em>n</em> points and you generate a least-squares polynomial of degree <em>n</em>-1, you've calculated the <u>interpolation</u> line that runs through all <em>n</em> points. This is also known as <u>polynomial curve fitting</u>.\n",
    "\n",
    "#### Calculate and plot the interpolation line for these points: (-2, 3), (-1, 5), (0, 1), (1, 4), and (2, 10).\n",
    "\n",
    "#### <strong>SOLUTION:</strong>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "xs = [-2, -1, 0, 1, 2]\n",
    "ys = [3, 5, 1, 4, 10]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1, -2,  4, -8, 16],\n",
       "       [ 1, -1,  1, -1,  1],\n",
       "       [ 1,  0,  0,  0,  0],\n",
       "       [ 1,  1,  1,  1,  1],\n",
       "       [ 1,  2,  4,  8, 16]])"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = np.array([[xs[i]**j for j in range(len(xs))] for i in range(len(xs))]) \n",
    "A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 3,  5,  1,  4, 10])"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "b = np.array(ys)\n",
    "b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 1.        , -1.25      ,  4.20833333,  0.75      , -0.70833333])"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a = solve(A, b)\n",
    "a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAV0AAAFUCAYAAACHh+9/AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAAAvMElEQVR4nO3dd3RVVdoG8CeQBKkGTCGFEGoSAiHIpSiO9KIytCB9kJpxRlSWAwqKWGYQEBnKCCoCGmFMUEBgEIIUBaV8cCGUCGSoQkKAUEIJAVLu98drKAOElHvOPufc57cW60jKve81ycPOPu/e283hcICIiPRRSnUBRESuhKFLRKQjhi4RkY4YukREOmLoEhHpiKFLRKQj94e8n/1k5DSdOnVCQkKC6jKI9OD2oHdwpEu6OXfunOoSiJRj6BIR6YihS0SkI4YuEZGOGLpERDpi6BIR6YihS0SkI4YuEZGOGLpERDpi6BIR6YihS0SkI4YuEZGOGLpUoCFDhsDX1xf169e/9bYLFy6gffv2qFOnDtq3b4+LFy8qrJCUy8wEfvoJsNuBkyeBmzdVV2RoDF0q0KBBg+7ZGWzSpElo27YtDh06hLZt22LSpEmKqiM9LUtMRYtJG1BjzPd4esIPsH/wMRAdDfj4AK1bA02aAMHBQLlywJAhEsB0D7eHnAbMrR0Jx48fR+fOnZGUlAQACA0NxU8//QR/f3+kpaWhVatWSE5Ofujj2Gw22O12rcslDSxLTMXYpfuQlZ2LGhdSMW3lR4hKO4Qsbz+U7RUNPPsskJMDnDkD7N0LfP454OYGvPIKMH48UKGC6pegtwdu7fiw/XSJ7nHmzBn4+/sDAPz9/XH27FnFFZHWpqxJRtbNHPTbk4BxG+biZmkPvPzH0Uhs3gG/vNnu3k8YPVrC9qOPgK1bgdWrXTF474vTC6SpOXPmwGazwWazIT09XXU5VEynLl7D+PWf44M1s7AzIBwdh3yM/9RridTLN+7/CdWrA7GxwKJFwJYtQJcuwLVr+hZtUAxdKjI/Pz+kpaUBANLS0uDr6/vAj42JiYHdbofdboePj49eJZKTvbf9awzZuQLzG3fBwN7v40xFbwBAgFfZgj/x+eeBr76SG23duwPXr2tfrMExdKnIunTpgtjYWABAbGwsunbtqrgi0tSECRj4UxwWNXoG77cdDoebxEZZj9IY3TH04Z/fvz8wfz7www/AyJHa1moCDF0qUN++ffHEE08gOTkZQUFBmDdvHsaMGYO1a9eiTp06WLt2LcaMGaO6TNLK3LnAuHHAgAEo8/lnCKxcDm4AAr3KYmKPBujWKLBwjzNoEPDaa8BnnwEbN2pZseGxe4F0w+4Fk9m3D2jaFPjDH4BVqwD3Et53z8wEGjSQx9mzByj7kKkJc+PBlERUBJmZQO/egJcXsGBByQMXAMqXB+bMAQ4dAt5/v+SPZ1IMXSK616uvAgcPAgsXAn5+znvcdu2AwYOBKVOAxETnPa6JMHSJ6G7ffAPMmwe8+SbQtq3zH3/qVKBKFcBF7wUwdInotowM4OWXZUnvu+9q8xyVK8tNtR9+AHbu1OY5DIyhS0S3vfUWcO6cdBk4Yx73Qf7yF6BSJcAF9+1g6BKR2LED+OQTYMQIoFEjbZ/r0UeBl14CliwBCrFvh5UwdIkIyM2V0aefn36dBSNHAmXKAB9+qM/zGQRDl4hkOmHnTmDaNBmF6sHXFxg6VFrSUlL0eU4DYOgSuborV+SmWcuW0purp1GjgLw8CXsXwdAlcnXTpwPp6cDkybIHrp5CQoBu3WS0m52t73MrwtAlcmXnzslChe7dgWbN1NTwwgsS+mvWqHl+nTF0iVzZxImy5Pcf/1BXQ6dOgLe37L/rAhi6RK7q5Elg1ixg4ECgXj11dXh4AP36AStWAC5wyClDl8hVvf8+4HBot/KsKAYOlFOEv/lGdSWaY+gSuaKTJ+XX+WHD5Ggd1R5/XEbbX32luhLNMXSJXNHUqTLKff111ZUINzcZ7W7ZAhw+rLoaTTF0iVxNerockd6/vzFGufkGDJDwtfhol6FL5GpmzgSysoA33lBdyd0CA4HWrYHFi1VXoimGLpEruXwZ+Ne/pC83PFx1Nffq0gU4cAA4ckR1JZph6BK5kk8+AS5dAsaOVV3J/XXuLNfvv1dbh4YYukSu4uZNWfLbrh1gs6mu5v5q1QLCwoCVK1VXohmGLpGr+PZb4PRpObXByDp3Bn76STbisSCGLpErcDhklBsaCnTsqLqagnXuLJvfrF2ruhJNMHSJXMHWrYDdLqf8ljL4j/2TT8rR7xadYjD4/30icorp0yXIBg5UXcnDeXjIJjjffy977VoMQ5fI6k6cAJYuBYYPB8qXV11N4XTuDJw9K6Nzi2HoElndrFlyHTFCbR1F0amTTINYcIqBoUtkZVlZwNy5cjpDcLDqagrvscdkbnfVKtWVOB1Dl8jKFi8GLlyQ487Npk0bIDFRFnNYCEOXyMo+/RSoWxdo1Up1JUXXqpXcSPvlF9WVOBVDl8iq9u6VrRJffFH/AyedoXlzwNMT2LhRdSVOxdAlsqpPPwXKlJGDH82obFmgaVOGLhGZwJUrcqx5795AlSqqqym+li2BnTsttSSYoUtkRXFxwNWrMrVgZi1bArm5wObNqitxGoYukRV9+inQsKHMi5rZk08C7u6WmmJg6BJZza5d0mo1fLg5b6DdqXx5oEkThi4RGdi8eXIDrV8/1ZU4R8uWwI4dQGam6kqcgqFLxTZt2jRERESgfv366Nu3L65fv666JMrKAv79byA6GqhcWXU1ztGyJZCTIzulWQBDl4olNTUVM2fOhN1uR1JSEnJzcxEfH6+6LFq6VFZwDR2quhLnadECKF1aNja3AIYuFVtOTg6ysrKQk5ODa9euISAgQHVJNG8eUKOGOVegPUjFikDjxpaZ12XoUrEEBgZi1KhRCA4Ohr+/Px599FF06NBBdVmu7cgR4McfgSFDjL9ReVE99ZRs85idrbqSErPYV4b0cvHiRSxfvhzHjh3DqVOnkJmZiYULF97zcXPmzIHNZoPNZkN6erqCSl3IF19It4JZV6AVpEkT4Pp1IClJdSUlxtClYlm3bh1q1KgBHx8feHh4oEePHtiyZcs9HxcTEwO73Q673Q4fHx8FlbqI3Fzgyy/l/LNq1VRX43xNm8p1xw61dTgBQ5eKJTg4GNu2bcO1a9fgcDiwfv16hIeHqy7LdW3YAKSmAoMHq65EGzVqyB6727errqTEGLpULM2aNUPPnj3x+OOPo0GDBsjLy0NMTIzqslxXbKycgdali+pKtOHmJlMMFhjpujkcjoLeX+A7iYrCZrPBbsEzr5S7fBmoWlUOnfz0U9XVaGf8eGDCBHm9xj/r7YFLATnSJTK7JUtkUYQVb6DdqWlT2dQ8MVF1JSXC0CUyu9hYoE4d829u8zBNmsjV5PO6DF0iMzt2TBYNvPCC+Te3eRg/Pzlc0+TzugxdIjNbuFDC9k9/Ul2JPpo04UiXiBRxOICvvgJatzbX8eol0aQJcPQocP686kqKjaFLZFbbtgGHD0vXgquwwCIJhi6RWS1cKIc39uihuhL9NG4s0ykMXSLS1c2bwKJFQNeusguXq6hUCQgLM/W8LkOXyIzWrJF5zQEDVFeiP5OvTGPoEpnRwoWAtzfgittpNmoEnDkjf0yIoUtkNpcuAStWAH36AB4eqqvRX2SkXPfuVVtHMTF0icxm6VLZW9YVpxYAoEEDuTJ0iUgX//43ULv27fYpV+PjA/j7M3SJSAepqbJ3bv/+1l/2W5CGDRm6RKSDRYtkJVq/fqorUSsyEti/35RnpjF0iczk668Bmw2oW1d1JWpFRkqv8n//q7qSImPoEplFcjKwcydHuYCpOxgYukRmERcn87i9e6uuRL3QUGmXY+gSkSYcDplaaN0aCAhQXY16np5AeDhDl4g0snMncOgQpxbuFBnJ0CUijcTFyejOlXYUe5jISCAlBbhwQXUlRcLQJTK63FwgPh549lmgcmXV1RhH/s20ffvU1lFEDF0io9u0CTh1CujbV3UlxmLSDgaGLpHRxcUBFSoAnTurrsRYqlaVndYYukTkNDdvAosXy2bl5cqprsZY3NxMuRyYoUtkZGvXAhcvyjaOdK/ISJnTzctTXUmhMXSJjCwuTm6eueJm5YVRrx6QlQWcOKG6kkJj6BIZ1bVrwLJlQM+e0i5G9woLk+uBA2rrKAKGLpFRrVwJZGZyaqEg4eFyPXhQbR1FwNAlMqr4eNmsu2VL1ZUY12OPSQcDQ5eISuTSJWDVKqBXL6B0adXVGFtYGEOXiEpo2TLgxg1OLRQGQ5eISiw+HggJAZo1U12J8YWFAWfPmmYPBoYukdGcOyf9uX36uPY5aIWV38GQnKy2jkJyV12AoWRlAd9/D/z2G5CWBpw/L3dH//AHoHFjtu2QPhYvlk1uOLVQOHe2jT3xhNpaCoGhC8hNi9mzgenT5dcUAChbFvDyAr788vbfhw4F3n4b8PVVVCi5hPh4CZL8DV2oYCEhMiAyybwupxeWLgWCg4E33wQaNQLWrZNll5mZsrPT6dMy8ujTB/jkE6BWLeDvf5fGdSJnS02VXcU4tVB4pUvLQZ0MXROYMUNW+4SHA7t2AQkJQNu2MsLN/4b38wOio4H584GkJKBdO2D8eKBFC5mGcGEZGRno2bMnwsLCEB4ejq1bt6ouyfy+/VaO5uHUQtGYqIPBNUM3Lw947TVg5EigWzdgwwYZ5T5MWBjw3Xcy73vsGNCkCfDzz1pXa1ivvvoqOnXqhIMHD2LPnj0Iz18dRMUXHy/fi6Ghqisxl/Bw4OhRabMzONcM3bffBqZNA155RUYWRd0y79lngf/7P9mIpG1bYOFCbeo0sMuXL2PTpk0YOnQoAMDT0xNeXl5qizK7Y8fk+4qj3KILC5Obj0eOqK7koVwvdFevBj74ABg2TKYXirvaJzRUfkCeegp44QVgyRLn1mlwR48ehY+PDwYPHoxGjRph2LBhyMzMVF2WuS1aJNdevdTWYUb5HQwmmGJwrdA9eRIYMEA2Pp45s+SP5+UF/Oc/0sDety+wZk3JH9MkcnJysGvXLvzlL39BYmIiypcvj0mTJt3zcXPmzIHNZoPNZkN6erqCSk0kLk5ankJCVFdiPnXrytUEu425TujevCkjiOxsmVIoW9Y5j1u+vKyRr1cP6N4d2LzZOY9rcEFBQQgKCkKz31dM9ezZE7t27brn42JiYmC322G32+Hj46N3meaxf7+cgMCpheKpUAGoVo0jXUOZPBnYtg2YNw+oU8e5j+3lBfzwAxAUJMF78qRzH9+AqlatimrVqiH591VA69evR7169RRXZWKLFgGlSnFqoSRM0sHg5nA4Cnp/ge80jbQ0CdpOnaTnVisHD0pHQ7160mtZpox2z2UAu3fvxrBhw3Dz5k3UrFkTX3zxBSoXcES4zWaD3W7XsUKTcDgkMIKCgPXrVVdjXq+8AnzxBXD5shF6nB9YgGusSHvnHZleuM+co1OFhckKtp49gVdfBT79VNvnUywqKooh6gyJicB//wuMGqW6EnMLDQWuXpVFTYGBqqt5IOtPLyQlyZTCSy8BtWtr/3zR0cDrrwOffXZ7CTFRQeLjAXd3oEcP1ZWYW/7Pt8HbxqwfuqNHA5UqAePG6fecEyYArVsDI0YY/huAFMvLk/ncDh3kFAQqPoauAaxdK0t7x43T9xva3R2IjZXrCy9I0zbR/WzbJifZ9u2ruhLzCw6WvnuGrkITJ8rNiREj9H/uatWAf/1LWsimTtX/+ckc4uKARx4BunRRXYn5eXgA1asDhw+rrqRA1g3dffuAH3+UwFXVRTBggMzTvf221EN0p5wc4JtvgOeekykwKrnatTnSVWbmTFkAMXy4uhrc3KSDwcsLGDRIfsiI8v34o+zfzKkF56lVi6GrxPnzsgnNgAFAlSpqa/HxAT7+WLaOnDVLbS1kLPHxQMWKsoESOUetWrIftoHPS7Nm6H7+OXD9OvDyy6orET17ysKMceNkk2qiGzdkk6Tu3Z23JJ1M0cFgvdDNyZERZZs2QIMGqqsRbm5SU06OLJogSkiQY6I4teBctWrJlaGro+++A1JSZEmgkdSsKTfUliyRTdDJtcXFAd7esh8zOU/NmnJl6Opo7lzp1+vcWXUl9xo1Sna4f/llmf4g15SZKVuCPv+8tDmR85QrBwQEGLptzFqhe/asbBjSv3/xNyfXkqen9O4eOyYnD5NrWrFCDjbl1II2DN7BYK3Q/fZbWf1l5G/mtm2Brl1lqXBamupqSIWvv5ZFOy1aqK7Emhi6OoqLAyIijHMD7UGmTJG713ruB0HGcP683ETr21f2zyXnq11bdhq7dk11Jfdlna/6iROy5NbIo9x8depIF8MXX0j/LrmOxYuli6VfP9WVWFd+B8PRo2rreADrhG58vFzNctzJuHG44VUFu6MHocYbK9Fi0gYsS2QPr+V9/bXcTG3YUHUl1mXwXl3rhG5cnBwQmf+vnMEtO3oVE5v3RdTxfWh7eDtSM7Iwduk+Bq+VnTghJ4r062eEkw2sKz8DDNrBYI3QPXgQ2L3bHFMLv5uyJhkLI9rhSJVAvL4xFqXzcpGVnYspa5JVl0Zayf9tjFML2qpcWf5wpKuhuDjTHep3KiMLOaXdMeXpgah7/gR6JK2/9XayqK+/Bpo3v93AT9ox8G5j1gjd5cuBp54C/P1VV1JoAV6y3j6h7pNI9A/Faz//G2Wyb9x6O1nMr78Ce/ZwlKsXA7eNmT90U1Plm/m551RXUiSjO4airEdpwM0Nk1oNgv/V84jZvRKjO4aqLo208PXXpvttzNRq1QKOHweys1VXcg/znwa8erVcTbY9XrdGclrplDXJ2I4G2BzaDK9sXwKPEJ4yYTl5eRK67doBfn6qq3ENtWrJQqkTJwx3c90aoVutmiyKMJlujQJvhS/6VZM2oo8+ktVqZB1btsio6+9/V12J66hRQ67HjxsudM09vXDzphw++eyz5m/BiYwEevcGZsyQPSTIOhYskI1YunVTXYnrqF5drr/9praO+zB36G7eDFy5AjzzjOpKnOO994CsLGDyZNWVkLPcuCHnoHXvDlSooLoa1xEUJHPox4+rruQe5g7dVatkazyr7EkaGgoMHCgbnvOECWtYtQrIyJCjo0g/Hh5AYCBHuk63ejXQsqW1RhDjx8uNF87rWsPChYCvr9xEI31Vr87QdarffpPeR6tMLeSrUQMYNkzOeTPgr0ZUBBcvAitXykpJd/PfszadkBBD/gyZN3RN2ipWKG+9JfNRHO2a2+LFcrOXUwtqVK8uR3fl5Kiu5C7mDd2EBPmXLNSCiwkCA4GYGODLL+WUCTKnBQvk+7NxY9WVuKbq1aVX99Qp1ZXcxZyhm5sLbNwo82RmbxV7kDFj5MghjnbN6ehR4OefgRdesO73qNGFhMjVYFMM5gzdpCS5I9yypepKtHPnaNegmzFTARYskLDl1II6Bu3VNWfobtok16efVluH1saMkRswHO2ai8MBfPUV0KaNrJYkNYKD5cqRrhNs3Cj/iuX/T7WqgADgz38GYmMNu2MS3cfmzfLbycCBqitxbY88AlStypFuiTkcMtK18tTCnd54Q0a7EyeqroQKKzYWKF8e6NFDdSVkwF5d84VucjKQnm79qYV8AQHA8OHyg2ywX5PoPrKyZNlvdLS1Fu2YVfXqhvu5MV/obtwoV1cJXUBGu6VKAZMmqa6EHmbFCuDyZelaIPVCQmR7x7w81ZXcYr7Q3bRJ5mnyT/x0BUFBwNChwPz5wMmTqqu5S25uLho1aoTOnTurLsUYYmPl5lmrVqorIUBGujdvAmfOqK7kFnOFrsMhI92WLV2v93HMGLkabLQ7Y8YMhIeHqy5DuWWJqeg2dhFyE9bgy9otsWxPmuqSCLjdNmagKQZzhe6xY7L7litNLeQLDgYGDwbmzjXMDmQpKSn4/vvvMWzYMNWlKLUsMRVjl+7Dk1tWobQjD/PrtMTYpfuwLNEYXyeXlr9AwkA308wVuq7Sn/sgY8fK3NSHH6quBAAwcuRIfPjhhyhVylzfRs42ZU0ysm7m4Pl9a7E1uAFOVPZHVnYupqxJVl0acaRbQps2AY89BtSrp7oSNUJCpPdzzhzg9GmlpaxcuRK+vr5o/JB9BebMmQObzQabzYb09HSdqtPXqYwsNEn5FTUupuGbBu3vejspVqECUKUKR7rF9ssvQIsWciffVY0dKzcGPvpIaRmbN2/GihUrEBISgj59+mDDhg0YcJ8lrzExMbDb7bDb7fDx8VFQqfYCvMqi1951uOJZFqtDn7zr7WQAISEM3WK5eBE4dAho3lx1JWrVrg307w988onSs9QmTpyIlJQUHD9+HPHx8WjTpg0WLlyorB6Vxj4ViOeSf8Z/wp/GdY9HAABlPUpjdEcL7oBnRgbr1TVP6Nrtcm3SRG0dRvDWW9KE/89/qq6EAHRO/gXlsm/gpyc7ww1AoFdZTOzR4PZJz6RW/qo0h0N1JQAAN0fBhRijSkA2fRk3Tka8Xl6qq1Gvb185leD4cZnnNgGbzQZ7/j+eVtKiBXDhArB/v+u1MprBjBnAyJHym6F+U1wP/EYwz0h3xw6gbl0Gbr5x44CrV4Hp01VX4tp+/RXYskWOWGLgGpPBtng0T+hu3w40baq6CuOIiAB69gRmzpS9hUmNzz8HPD257NfI8ncjTElRW8fvzBG6qalAWhrnc//XuHGyzn/mTNWVuKasLNk3t0cPwNtbdTX0IEFBcmXoFsH27XLlSPduDRsCXbsC06ZJ+JK+liyRewwxMaoroYJ4e8tvIwbZt8Qcobtjh+wpGxWluhLjefttmV74+GPVlbiezz+XFj5ubmNspUrJ8Vcc6RbB9u1AZKTsBE93a9wYeO45aR+7elV1Na7j4EFZITl8OG+gmUFQEEO30PLypEeX87kP9vbbwPnzwKxZqitxHZ9/Dnh4AIMGqa6ECoOhWwSHDgGXLnE+tyDNmgGdOsnSYI52tZeVJac0d+0K+PqqroYKIz90DbBAwvihu2OHXDnSLdg77wDnzsnyYNLWokWyGOKvf1VdCRVWUJDsWXLunOpKTBK65coB3Ci7YM2bAx07AlOmAJmZqquxttmzZac73kAzDwO1jRk/dLdvl5tF7u6qKzG+d96RQzs52tXOjh3y569/5Q00M2HoFlJuLrBnj4QuPdwTTwAdOsgm5xztamPWLNmj9U9/Ul0JFQVDt5AOH5abFg0bqq7EPPJHu7Nnq67Ees6dA+LjJXArVVJdDRWFnx9QujRD96H27JErQ7fwnnzy9miXnQzONX8+cOMGb6CZUenSQEAAQ/eh9u6V/1m8iVY0778vo7J//Ut1JdaRmwt8+qmcz1e/vupqqDgM0qtr2NBdlpiKX5ZuQHLlQLSYvoUnqxZFs2aySm3KFO7J4CwrVshp1C+/rLoSKi6G7oPlH2kdknoYB31qIDUji0daF9V778lmLNxv1zmmTZN9Wbt1U10JFZdBFkgYMnSnrEmGx5VLCLqcjoO+IQDAI62LqnFjCYh//lPCl4pv507g55+BV15h66KZBQUB164p33/akKF7KiMLYenHAQAHfGrc9XYqgnfflSXUU6eqrsTcpk+XNrGhQ1VXQiVRrZpcFU8xGDJ0A7zKIuzsMQDAgd9HuvlvpyJo2BDo1UtCQ+HJwaZ26pS0iQ0dCjz6qOpqqCQM0qtryNAd3TEUDc79hgtlK+FMBTl0kUdaF9P77wPXrwMffKC6EnOaPVs6F3gDzfwYug/WrVEg2uacxjH/mnBzc+OR1iURGirbD37yCXDihOpqzOXaNWkT69oVqFVLdTVUUlWryobmDN37yM1FlSPJaNylFY5Neg6bx7Rh4JbE+PFyfe89tXWYzbx5sk/xqFGqKyFn8PCQ4GXo3seRI7L8NzJSdSXWEBwsq6i+/BJIZgdIoWRny/7Ef/gD0KKF6mrIWQzQq2vM0OXyX+cbO1a2yHzrLdWVmENcnEzHjBmjuhJyJobuA+zdK3Mv9eqprsQ6fH3l1+QlS4Bt21RXY2x5ecDkyUCDBsAzz6iuhpyJofsAe/fKDSAeROlcf/ub7LY0apTyVTmGtnIlsH+/jHK5Z661BAXJ0niFy+ONGbp79nBqQQsVKsjNtM2bgeXLVVdjTA4HMHEiEBIiPc5kLfltY6nqthQwXuheuQL89ht3ctLK0KFAWBjwxhtys4jutmGDTL+MHs0lv1aUH7onTyorwXihm393nds5asPdXeYr//tfYO5c1dUYi8Mhm8AHBQFDhqiuhrQQECDXU6eUlWC80D1wQK4MXe388Y+yL+w77yjf/MNQ1q6VqZc33+T9BKvy95drWpqyEowZuu7uQO3aqiuxLjc32arw3Dng739XXY0x5I9yq1XjKNfKypWTPTQYunc4cEAC18NDdSXW9vjjMr87cyYXTABAQoLM5Y4bB5Qpo7oa0pK/P6cX7nLwoNzoIe394x9A2bLSSubK8ke51avLPhVkbf7+HOnekp0tJwBzPlcffn6yL8P338tIz1UtXw7s2CGjXE9P1dWQ1hi6dzh8GMjJYejq6ZVXgDp1gJEj5aRbV5OdLe1zYWEc5bqKgAAJXUULhIwVuuxc0J+n5+153Y8+Ul2N/ubOlfa5yZPZl+sq/P1lj2lFnTvGDF3O6eqrUyfg+edljvfo0UJ9ysmTJ9G6dWuEh4cjIiICM2bM0LhIDVy5IkcaPf20tNGRa1DcNmas0D14UBrTK1RQXYnrmTZNRnojRhTq1y53d3dMnToVBw4cwLZt2zBr1izs379fh0KdaMoUOcZoyhTuseBKGLp3OHCAUwuqBAZKz+7q1cDSpQ/9cH9/fzz++OMAgIoVKyI8PBypCtezF9mpU3JgZ69eQNOmqqshPeWvSnP50M3Lk5EuQ1edESOAqCg5D6wI813Hjx9HYmIimjVrpllpTve3v8nZZxMnqq6E9JY/0lXUq2uc0E1JATIzGboqubvLjaWzZ4HXXivUp1y9ehXR0dGYPn06KlWqdM/758yZA5vNBpvNhvT0dGdXXDzr18sJv2PGADVrqq6G9FaxIlC+PEe67FwwiMaNpYXqiy9kqqEA2dnZiI6ORv/+/dGjR4/7fkxMTAzsdjvsdjt8fHy0qLhobt6UEX3NmvI6yTUp7NU1TugePChXdi6oN368nNoREwNcunTfD3E4HBg6dCjCw8PxWiFHxYYwbZp8r82cKavxyDXl9+oqYJzQPXAAqFxZjpUhtcqUkUMsT5164BLhzZs3Y8GCBdiwYQOioqIQFRWFVatW6VtnUZ04Abz/vhyp/txzqqshlRSOdI3TDZ7fucDWHWNo0gR4/XVg0iQJqO7d73r3U089BYeZjvxxOIBhw+T7a/p01dWQago3vTHWSJfzucby3nuAzSa7kSncad8pPvtM9sudMkWO4iHX5u8vN+6vXNH9qY0RuhcvAunpnM81Gk9POYo8Oxvo319arMzo6FE5jLNdO+DFF1VXQ0agsFfXGKF7+LBc69RRWwfdq3ZtYPZs4OefZZmw2eTlyabkpUoB8+Zx+oqEwlVpxgpdnhZhTH/6EzBggNyEWrNGdTVFM3UqsHGjdC0EB6uuhoxC4QIJY4UuG9WN65NPgIgIoE+f218vo9u0CRg7FoiO5hE8dDeOdA/LRjfsmzSuChVks+9SpaTlSsENiCI5fRro3Vv+IZ8/n9MKdDcvLzl81KVDl1MLxlejBvDNN7L37sCBMl9qRDk5QN++srBjyRLgPsuTycW5uSnr1TVO6PImmjm0bSvzpMuWyakTRuvVdThkQcdPP8mUSIMGqisio1LUq6t+ccTly7LBCke65vHqq0BqqvS8envLRuBG8dFHssR35EjghRdUV0NG5u8P/Pqr7k+rPnTZuWBOkycD58/LAorHHpPtIFVbuFBW0fXuLaNxooIEBADr1un+tAxdKh43N1nldeHC7WmGV15RV09CAjB4MNCqFRAbKzf8iAri7y/z/llZut7EV/+dmR+6tWqprYOKzt1d9qXt3l2mHN59V80c7+LFQJcuQP36Mtdcpoz+NZD5KGobM0bo+vvLpsJkPmXKSEfD4MEy1fDyy/ouF547V6YTmjYFfvwRePRR/Z6bzE3RAgljhC6nFszN3V2W2P7tb8CsWUDHjsCZM9o+Z26uhPzw4UCHDsAPP0jvJVFhKdp/gaFLzuHmJp0D8+YBmzfLWWsbNwIAliWmosWkDdiXegktJm3AssQSHmCZliZB++670i+8fDlQrlyJXwK5mKpV5Xr6tK5PqzZ0MzPlB4g9utYxZAjwf/8nCxLatMGRvkMw8d9bkJqRBQBIzcjC2KX7ihe8DocEbMOGwNatstLsyy9lNzSionrsMaB0ae1/K/sfakP3yBG5cqRrLZGRgN0O/PnPCFkUi9Wzh6Hf7tVw+/0mW1Z2LqasSS7aY+7YAbRuDXTrBvj5yeMPHszlvVR8pUoBPj4uFrpsF7OuihWB2bPxx0HTccg7GB+smYWw9OMYtekrBF06g1O/j3wLlJUlN+n++Ee5UbZ/v8wZ79olZ7gRlZSvryzO0pHaPl22i1nepboR6O0zEa2O2nHtP1Pwl22L8det3+KYXwhwfhnQrJncRS5TRqYJTp2ScP31Vznp4coVef/48bIRecWKql8SWYmfn+4jXfWh6+vLDUksbHTHUIxdug8/1WqCtMoBeKr7W+h14Cf0u/mbbEYzd+69n1SqlPxD3LOnnFjRqpXMvRE5m58fcOiQrk+pNnQPHeLUgsV1axQIAJiyJhlpAEoFB6PG8AnwaxQoN8aOHJHjmq5fB27ckB+COnVk2z0ireWPdB0O3e4PqB/ptmmjtATSXrdGgejWKBC2xY9i85g7vt5ubvxHl9Ty85N7B1ev6jZ1pe5GWlYWkJLCHzoiUsfPT646zuuqC93jx+XKm2hEpIpLhm6NGspKICIX5+srVx3bxtSF7rFjcg0JUVYCEbk4lxvplilz+0UTEenNx0euLhO61atzs2kiUsfDQ/ZgcJnQ5dQCEamm86o0taHLm2hEpJpLhO7Vq0B6Oke6RKSen58LdC/89ptcGbpEpJqvrwuMdPN7dBm6RKSanx9w+bLs/6EDhi4RuTade3XVhO6xY7KLFHt0iUg1lwjd/HYxHrVCRKq5VOgSEanG0CWzSEhIQGhoKGrXro1JkyapLoeoeHTe9Eb/0L1yBTh/nqFrcrm5uXjppZewevVq7N+/H3Fxcdi/f7/qsoiK7pFH5Mgwy4502blgCdu3b0ft2rVRs2ZNeHp6ok+fPli+fLnqsoiKR8dVaW4Oh+OB7+zUqZPj3Llzzn3GS5fkmJ6wMKB8+Yd+eHp6OnzydwKyKDO+xosXL+Ly5cuoXr06AOD8+fPIzMxEcHDwXR+Xnp6O/O+hGzduICoqSu9SdWXGr2VRWfI1JifLNTT01ptK8jp37ty5xuFwdLrvOx0OR0F/nG/mTIcDcDjOnCnUhzdu3FiTMozEjK/xm2++cQwdOvTW37/66ivHiBEjCvyccuXKaV2Wcmb8WhaVJV9jdLTDER5+15tK+DofmKtqphfKlr29jyWZUlBQEE6ePHnr7ykpKQgICFBYEVEJ6Di9oH/oHjvGHl0LaNKkCQ4dOoRjx47h5s2biI+PR5cuXVSXRVQ8fn7AhQtAdrbmT6X/EexFbBeLiYnRrBSjMONrdHd3x8cff4yOHTsiNzcXQ4YMQURERIGf4+3trVN16pjxa1lUlnyN+W1j6enA77+xafU6C7yRBqDAdxZLlSpAnz7A7NlOf2gyNpvNBrvdrroMont99x3QowewaxfQqJEzHvGBv8rrO71w6RJw8SLbxYjIWHRclaZv6HIfXSIyIsuG7okTcv2fXs6CjB49GmFhYYiMjET37t2RkZGhTW2Kffvtt4iIiECpUqUs9yt4/nLhpKQkyy4XHjJkCHx9fVG/fn3VpWjm5MmTaN26NcLDwxEREYEZM2aoLsl5fg/d7JQUNG3aFA0bNkRERATeeecdpz+VvqGb32JUrVqhP6V9+/ZISkrC3r17UbduXUycOFGj4tSqX78+li5diqefflp1KU5153LhiIgIyy4XHjRoEBISElSXoSl3d3dMnToVBw4cwLZt2zBr1izrfC3LlwceeQTuGRnYsGED9uzZg927dyMhIQHbtm1z6lPpH7ru7kDVqoX+lA4dOsDdXZosmjdvjpSUFK2qUyo8PByhd6yGsYo7lwu7ublZdrnw008/jSpVqqguQ1P+/v54/PHHAQAVK1ZEeHg4UlNTFVflJG5ugI8P3M6dQ4UKFQAA2dnZyM7OhpuT21v1D92AAKB06WJ9+vz58/HMM884uSjSUmpqKqrd8ZtNUFCQdX5QXdjx48eRmJiIZs2aqS7Feby9gfR05ObmIioqCr6+vmjfvr3TX6O+fbopKUBQ0D1vbteuHU6fPn3P2ydMmICuXbve+m93d3f0799f8zK1UpjXaTX3a0l09siB9HX16lVER0dj+vTpqFSpkupynMfHB0hPR+nSpbF7925kZGSge/fuSEpKcupcvb6he/IkYLPd8+Z169YV+GmxsbFYuXIl1q9fb+of2Ie9TivicmFryc7ORnR0NPr3748ePXqoLse5fHxkM67feXl5oVWrVkhISHBq6Oo3veBwyEi3CDfRALnzPXnyZKxYsQLlypXTqDjSyp3LhR0OB5cLm5jD4cDQoUMRHh6O1157TXU5zufjg7yzZ291SGVlZWHdunUICwtz6tPoF7rp6cCNG0UO3REjRuDKlSto3749oqKi8OKLL2pUoFrfffcdgoKCsHXrVjz33HPo2LGj6pKc4s7lwklJSejVq9dDlwubUd++ffHEE08gOTkZQUFBmDdvnuqSnG7z5s1YsGABNmzYgKioKERFRWHVqlWqy3Ieb2+UunoVHVq2RGRkJJo0aYL27dujc+fOTn0a/ZYB79wpUwtLlwLduzvtYck8uAyYDG3OHODPf5Zp0PvceyoiAywDLkaPLhGRbvK3m01P1/RpGLpERMDt0HX2aTn/Q9/Q9fTk5uVEZEz5W49aaqQbFASUUnPqOxFRgSw5vcCpBSIyqsqVZbWsZUK3GD26RES6KVUKeOwxi8zp5uYCqakMXSIytt/3X9CSPqF75gyQk+OM3jciIu38vv+ClvQJXbaLEZEZMHSJiHTE0CUi0pGPD3DhgtyH0oh+oVu2rBy/TkRkVN7esiPihQuaPYV+oVutmhyJQURkVDoskNA3dImIjMwyocuFEURkBpYI3ZwcIC2NoUtExpe/6Y2Gq9K0D91Tp4C8PIYuERmfDjuNaR+6+e1iXI1GREbn6Qk8+qjJQzc1Va4MXSIyA40XSOgzvQAAPHabiMzA29sCc7qenlwYQUTmYImRbkAAF0YQkTlYInQDAzV/GiIip8gPXYdDk4fXb6RLRGQG3t5AdjZw+bImD8/QJSK6k8ZHsWsbuleuyB+GLhGZhcZLgbUN3bQ0uTJ0LWX06NEICwtDZGQkunfvjoyMDNUlETmPqUOXPbqW1L59eyQlJWHv3r2oW7cuJk6cqLokIudh6JLRdOjQAe7u7gCA5s2bIyUlRXFFRE6k8f4LDF0qkfnz5+OZZ5554PvnzJkDm80Gm82GdI3PniJyivLlgUce0exGmrsmj5rv1Cl5ARUravo05Hzt2rXD6dOn73n7hAkT0LVr11v/7e7ujv79+z/wcWJiYhATEwMAsNls2hRL5ExubsD+/benGZxM+9DlajRTWrduXYHvj42NxcqVK7F+/Xq48etLVlOjhmYPrW3opqZyasGCEhISMHnyZGzcuBHlypVTXQ6RqWg/p8vQtZwRI0bgypUraN++PaKiovDiiy+qLonINLQb6Toc3HfBog4fPqy6BCLT0m6km5EBXL/OkS4R0R20C122ixER3YOhS0SkI4YuEZGOtA9df3/NnoKIyGy0DV0vL4B9nEREt2gbupxaICK6C0OXiEhHDF0iIh1pE7p5eXJqBEOXiOgu2oTuuXNymiZDl4joLtqEbn67GPddICK6i7ahy5EuEdFdtAndcuWAtm2BatU0eXgiIrPSZmvHVq3kDxER3UXbTcyJiOguDF0iIh0xdImIdMTQJSLSEUOXiEhHDF0iIh0xdImIdMTQJSLSEUOXiEhHDF0iIh0xdImIdMTQJSLSEUOXiEhHbg6HQ3UN5CLc3NwSHA5HJ9V1EKnE0CUi0hGnF4iIdMTQJSLSEUOXiEhHDF0iIh0xdImIdPT/+9OtjLY0UcEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "pr.plot_polynomial(xs, ys, a, degree=len(xs)-1, x_min=-2.3, x_max=3, y_axis_pos=0, size=(6, 6))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr>\n",
    "\n",
    "#### <strong>PROBLEM 3.a.</strong> Amalgamated Widgets, Inc. ranks its employees by productivity: 1 (lowest), 2, 3, and 4 (highest). The following table shows for each rank, at the start of a new training program, the percentage of employees, the percentage productivity loss per employee, and the total percentage productivity loss.\n",
    "\n",
    "<table>\n",
    "    <tr><th>Rank</th><th>Employees</th><th>x</th><th>Productivity loss</th><th>=</th><th>Total loss</th></tr>\n",
    "    <tr><td align=center>1</td><td align=center>30%</td><td></td><td align=center>40%</td><td></td><td align=center>12.00%</td></tr>\n",
    "    <tr><td align=center>2</td><td align=center>35%</td><td></td><td align=center>25%</td><td></td><td align=center>&nbsp;8.75%</td></tr>\n",
    "    <tr><td align=center>3</td><td align=center>25%</td><td></td><td align=center>15%</td><td></td><td align=center>&nbsp;3.75%</td></tr>\n",
    "    <tr><td align=center>4</td><td align=center>10%</td><td></td><td align=center>&nbsp;5%</td><td></td><td align=center>&nbsp;0.50%</td></tr>\n",
    "    <tr><th align=right colspan=5>Overall total</th><th align=center>25.00%</th></tr>\n",
    "</table>\n",
    "\n",
    "#### The training supervisor provided the following table that predicts what percentage of employees will move from one rank to another or stay in the same rank after each hourly training session, as determined by a test given at the end of the hour:\n",
    "\n",
    "<table>\n",
    "    <tr><th colspan=5 align=center>FROM:</th></tr>\n",
    "    <tr><th>TO:</th><th>1</th><th>2</th><th>3</th><th>4</th></tr>\n",
    "    <tr><th>1</th><td align=center>50%</td><td align=center>30%</td><td align=center>10%</td><td align=center>10%</td>\n",
    "    <tr><th>2</th><td align=center>30%</td><td align=center>40%</td><td align=center>20%</td><td align=center>10%</td>\n",
    "    <tr><th>3</th><td align=center>20%</td><td align=center>30%</td><td align=center>20%</td><td align=center>10%</td>\n",
    "    <tr><th>4</th><td align=center>&nbsp;&nbsp;0%</td><td align=center>&nbsp;&nbsp;0%</td><td align=center>50%</td><td align=center>70%</td>\n",
    "</table>\n",
    "\n",
    "#### If training ends after two consecutive sessions with no rank changes to within three decimal places, how many hours will the training program last?\n",
    "\n",
    "#### <strong>SOLUTION:</strong>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.5, 0.3, 0.1, 0.1],\n",
       "       [0.3, 0.4, 0.2, 0.1],\n",
       "       [0.2, 0.3, 0.2, 0.1],\n",
       "       [0. , 0. , 0.5, 0.7]])"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Probability matrix.\n",
    "P = np.array([[0.5, 0.3, 0.1, 0.1],\n",
    "              [0.3, 0.4, 0.2, 0.1],\n",
    "              [0.2, 0.3, 0.2, 0.1],\n",
    "              [0.0, 0.0, 0.5, 0.7]])\n",
    "\n",
    "P"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.3 , 0.35, 0.25, 0.1 ])"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Initial state (percentage of employees in each rank).\n",
    "X = np.array([0.30, 0.35, 0.25, 0.10])\n",
    "\n",
    "X"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "After hour 1:\n",
      "[0.29  0.29  0.225 0.195]\n",
      "\n",
      "After hour 2:\n",
      "[0.274  0.2675 0.2095 0.249 ]\n",
      "\n",
      "After hour 3:\n",
      "[0.2631 0.256  0.2019 0.2791]\n",
      "\n",
      "After hour 4:\n",
      "[0.2564 0.2496 0.1977 0.2963]\n",
      "\n",
      "After hour 5:\n",
      "[0.2525 0.2459 0.1953 0.3062]\n",
      "\n",
      "After hour 6:\n",
      "[0.2502 0.2438 0.194  0.312 ]\n",
      "\n",
      "After hour 7:\n",
      "[0.2488 0.2426 0.1932 0.3154]\n",
      "\n",
      "After hour 8:\n",
      "[0.2481 0.2419 0.1927 0.3174]\n",
      "\n",
      "After hour 9:\n",
      "[0.2476 0.2414 0.1924 0.3185]\n",
      "\n",
      "After hour 10:\n",
      "[0.2473 0.2412 0.1923 0.3192]\n",
      "\n",
      "After hour 11:\n",
      "[0.2472 0.2411 0.1922 0.3196]\n",
      "\n",
      "After hour 12:\n",
      "[0.2471 0.241  0.1921 0.3198]\n",
      "\n",
      "After hour 13:\n",
      "[0.247  0.2409 0.1921 0.3199]\n",
      "\n",
      "After hour 14:\n",
      "[0.247  0.2409 0.1921 0.32  ]\n",
      "\n",
      "Training lasts 14 hours.\n"
     ]
    }
   ],
   "source": [
    "%precision 4\n",
    "# for printing only\n",
    "\n",
    "prev_X = np.zeros(len(X))\n",
    "hour = 0\n",
    "\n",
    "# Loop until the steady state is achieved.\n",
    "# Limit 100 in case of nonconvergence.\n",
    "for i in range(100):\n",
    "    hour += 1\n",
    "    \n",
    "    prev_X = X    \n",
    "    X = P@X\n",
    "    \n",
    "    print()\n",
    "    print(f'After hour {hour}:')\n",
    "    print(X)\n",
    "    \n",
    "    if (abs(X - prev_X) < 1.0e-4).all(): break\n",
    "\n",
    "print()\n",
    "print(f'Training lasts {hour} hours.')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr>\n",
    "\n",
    "#### <strong>PROBLEM 3.b.</strong> Compute directly (without iteration) the percentages of employees in each rank after all the training sessions have ended. Revise the first table in Problem 3.a to show the new percentage of employees in each rank. Assume the productivity loss of each rank hasn't changed, and include the new total loss percentage per rank and the overall total loss percentage. What percentage improvement resulted from the training?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-0.5,  0.3,  0.1,  0.1],\n",
       "       [ 0.3, -0.6,  0.2,  0.1],\n",
       "       [ 0.2,  0.3, -0.8,  0.1],\n",
       "       [ 0. ,  0. ,  0.5, -0.3]])"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%precision 3\n",
    "\n",
    "A = P - np.identity(len(P))\n",
    "A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-0.5,  0.3,  0.1,  0.1],\n",
       "       [ 0.3, -0.6,  0.2,  0.1],\n",
       "       [ 0.2,  0.3, -0.8,  0.1],\n",
       "       [ 0. ,  0. ,  0.5, -0.3],\n",
       "       [ 1. ,  1. ,  1. ,  1. ]])"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = np.insert(A, len(P), 1, axis=0)\n",
    "A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0., 0., 0., 0., 1.])"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "b = np.append(np.zeros(len(P)), np.array([1]))\n",
    "b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-0.5,  0.3,  0.2,  0. ,  1. ],\n",
       "       [ 0.3, -0.6,  0.3,  0. ,  1. ],\n",
       "       [ 0.1,  0.2, -0.8,  0.5,  1. ],\n",
       "       [ 0.1,  0.1,  0.1, -0.3,  1. ]])"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "AT = A.T\n",
    "AT"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1.38, 0.73, 0.85, 1.  ],\n",
       "       [0.73, 1.54, 0.67, 1.  ],\n",
       "       [0.85, 0.67, 1.94, 0.8 ],\n",
       "       [1.  , 1.  , 0.8 , 1.12]])"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ATA = AT@A\n",
    "ATA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 2.300e+00,  5.644e-01, -2.098e-01, -2.408e+00],\n",
       "       [ 5.644e-01,  1.687e+00, -1.376e-03, -2.009e+00],\n",
       "       [-2.098e-01, -1.376e-03,  7.515e-01, -3.482e-01],\n",
       "       [-2.408e+00, -2.009e+00, -3.482e-01,  5.085e+00]])"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ATAinv = inv(ATA)\n",
    "ATAinv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1., 1., 1., 1.])"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ATb = AT@b\n",
    "ATb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The steady state:\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([0.247, 0.241, 0.192, 0.32 ])"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "print(\"The steady state:\")\n",
    "\n",
    "xhat = ATAinv@ATb\n",
    "xhat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.247, 0.241, 0.192, 0.32 ])"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Verify.\n",
    "P@xhat"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.247, 0.241, 0.192, 0.319])"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Should equal what we calculated by iteration.\n",
    "X"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.9999999999999992"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Check that the percentages sum to 1.\n",
    "np.sum(xhat)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.4 , 0.25, 0.15, 0.05])"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Vector of losses per rank.\n",
    "losses = np.array([0.40, 0.25, 0.15, 0.05])\n",
    "losses"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Rank  Employees X Productivity Loss = Total Loss\n",
      "  1     24.7%         40.0%             9.88%\n",
      "  2     24.1%         25.0%             6.02%\n",
      "  3     19.2%         15.0%             2.88%\n",
      "  4     32.0%          5.0%             1.60%\n",
      "\n",
      "                       Overall total:  20.38%\n",
      "                         Improvement:  18.48%\n"
     ]
    }
   ],
   "source": [
    "print('Rank  Employees X Productivity Loss = Total Loss')\n",
    "\n",
    "old_total  = 0.25\n",
    "new_total = 0\n",
    "\n",
    "for rank in range(4):\n",
    "    rank_loss = xhat[rank]*losses[rank]\n",
    "    new_total += rank_loss\n",
    "    print(f'{rank+1:3d} {100*xhat[rank]:8.1f}% {100*losses[rank]:12.1f}% {100*rank_loss:16.2f}%')\n",
    "\n",
    "print()\n",
    "print(f'                       Overall total: {100*new_total:6.2f}%')\n",
    "print(f'                         Improvement: {100*(old_total - new_total)/old_total: 6.2f}%')"
   ]
  },
  {
   "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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}