Chapter 1 - Linear Equation System.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib widget\n",
    "from mayavi import mlab"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "import scipy as sp\n",
    "import scipy.linalg\n",
    "import sympy as sy\n",
    "sy.init_printing() \n",
    "np.set_printoptions(precision=3)\n",
    "np.set_printoptions(suppress=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mlab.init_notebook(backend='x3d')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We start from plotting basics in Python environment, in the meanwhile refresh the system of linear equations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <font face=\"gotham\" color=\"purple\"> Visualisation of A System of Two Linear Equations </font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider a linear system of two equations:\n",
    "\\begin{align}\n",
    "x+y&=6\\\\\n",
    "x-y&=-4\n",
    "\\end{align}\n",
    "Easy to solve: $(x, y)^T = (1, 5)^T$. Let's plot the linear system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "6a6cda854af743ebbc5838943dbf4d16",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x = np.linspace(-5, 5, 100)\n",
    "y1 = -x + 6\n",
    "y2 = x + 4\n",
    "\n",
    "fig, ax = plt.subplots(figsize = (12, 7))\n",
    "ax.scatter(1, 5, s = 200, zorder=5, color = 'r', alpha = .8) \n",
    "\n",
    "ax.plot(x, y1, lw =3, label = '$x+y=6$')\n",
    "ax.plot(x, y2, lw =3, label = '$x-y=-4$')\n",
    "ax.plot([1, 1], [0, 5], ls = '--', color = 'b', alpha = .5)\n",
    "ax.plot([-5, 1], [5, 5], ls = '--', color = 'b', alpha = .5)\n",
    "ax.set_xlim([-5, 5])\n",
    "ax.set_ylim([0, 12])\n",
    "\n",
    "ax.legend()\n",
    "s = '$(1,5)$'\n",
    "ax.text(1, 5.5, s, fontsize = 20)\n",
    "ax.set_title('Solution of $x+y=6$, $x-y=-4$', size = 22)\n",
    "ax.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <font face=\"gotham\" color=\"purple\"> How to Draw a Plane </font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before drawing a plane, let's refresh the logic of Matplotlib 3D plotting. This should be familiar to you if you are a MATLAB user. \n",
    "\n",
    "First, create meshgrids."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "x, y = [-1, 0, 1], [-1, 0, 1]\n",
    "X, Y = np.meshgrid(x, y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Mathematically, meshgrids are the coordinates of <font face=\"gotham\" color=\"red\">Cartesian product</font>. To illustrate, we can plot all the coordinates of these meshgrids"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "375bc50022f44bc89ee0463add9f8606",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize = (12, 7))\n",
    "ax.scatter(X, Y, s = 200, color = 'red')\n",
    "ax.axis([-2, 3.01, -2.01, 2])\n",
    "ax.spines['left'].set_position('zero') # alternative position is 'center'\n",
    "ax.spines['right'].set_color('none')\n",
    "ax.spines['bottom'].set_position('zero')\n",
    "ax.spines['top'].set_color('none')\n",
    "ax.grid()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Try a more complicated meshgrid."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "9d8b829bf13249ce8e14e2bc73b46c82",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x, y = np.arange(-3, 4, 1), np.arange(-3, 4, 1)\n",
    "X, Y = np.meshgrid(x, y)\n",
    "\n",
    "fig, ax = plt.subplots(figsize = (12, 12))\n",
    "ax.scatter(X, Y, s = 200, color = 'red', zorder = 3)\n",
    "ax.axis([-5, 5, -5, 5])\n",
    "\n",
    "ax.spines['left'].set_position('zero') # alternative position is 'center'\n",
    "ax.spines['right'].set_color('none')\n",
    "ax.spines['bottom'].set_position('zero')\n",
    "ax.spines['top'].set_color('none')\n",
    "ax.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now consider the function $z = f(x, y)$, $z$ is in the $3rd$ dimension. Though Matplotlib is not meant for delicate plotting of 3D graphics, basic 3D plotting is still acceptable. \n",
    "\n",
    "For example, we define a simple plane as \n",
    "$$z= x + y$$\n",
    "Then plot $z$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "9a42f621a1404ea5965c21db201dc2b3",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "Z = X + Y\n",
    "fig = plt.figure(figsize = (9,9))\n",
    "ax = fig.add_subplot(111, projection = '3d')\n",
    "ax.scatter(X, Y, Z, s = 100, label = '$z=x+y$')\n",
    "ax.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Or we can plot it as a surface, Matplotlib will automatically interpolate values among the Cartesian coordinates such that the graph will look like a surface."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "dd024b438d1e4be4bc39f52632f4f071",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure(figsize = (9, 9))\n",
    "ax = fig.add_subplot(111, projection = '3d')\n",
    "ax.plot_surface(X, Y, Z, cmap ='viridis') # MATLAB default color map\n",
    "ax.set_xlabel('x-axis')\n",
    "ax.set_ylabel('y-axis')\n",
    "ax.set_zlabel('z-axis')\n",
    "ax.set_title('$z=x+y$', size = 18)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <font face=\"gotham\" color=\"purple\"> Visualisation of A System of Three Linear Equations  </font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We have reviewed on plotting planes, now we are ready to plot several planes all together."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider this system of linear equations\n",
    "\\begin{align}\n",
    "x_1- 2x_2+x_3&=0\\\\\n",
    "2x_2-8x_3&=8\\\\\n",
    "-4x_1+5x_2+9x_3&=-9\n",
    "\\end{align}\n",
    "And solution is $(x_1, x_2, x_3)^T = (29, 16, 3)^T$. Let's reproduce the system visually."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "ab9ad0b1a8f84be4a9feae9f9a22039c",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "x1 = np.linspace(25, 35, 20)\n",
    "x2 = np.linspace(10, 20, 20)\n",
    "X1, X2 = np.meshgrid(x1, x2)\n",
    "\n",
    "fig = plt.figure(figsize = (9, 9))\n",
    "ax = fig.add_subplot(111, projection = '3d')\n",
    "\n",
    "X3 = 2*X2 - X1\n",
    "ax.plot_surface(X1, X2, X3, cmap ='viridis', alpha = 1) \n",
    "\n",
    "X3 = .25*X2 - 1\n",
    "ax.plot_surface(X1, X2, X3, cmap ='summer', alpha = 1)\n",
    "\n",
    "X3 = -5/9*X2 + 4/9*X1 - 1\n",
    "ax.plot_surface(X1, X2, X3, cmap ='spring', alpha = 1)\n",
    "\n",
    "ax.scatter(29, 16, 3, s = 200, color = 'black')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are certain there is a solution, however the graph does not show the intersection of planes. The problem originates from Matplotlib's rendering algorithm, which is not designed for drawing genuine 3D graphics. It merely projects 3D objects onto 2D dimension to imitate 3D features.\n",
    "\n",
    "Mayavi is much professional in rendering 3D graphics, we give an example here. If not installed, run ```conda install -c anaconda mayavi```."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mlab.clf()\n",
    "X1, X2 = np.mgrid[-10:10:21*1j, -5:10:21*1j]\n",
    "X3 = 6 - X1 - X2\n",
    "mlab.mesh(X1, X2, X3,colormap=\"spring\")\n",
    "X3 = 3 - 2*X1 + X2\n",
    "mlab.mesh(X1, X2, X3,colormap=\"winter\")\n",
    "X3 = 3*X1 + 2*X2 -4\n",
    "mlab.mesh(X1, X2, X3,colormap=\"summer\")\n",
    "mlab.axes()\n",
    "mlab.outline()\n",
    "mlab.points3d(1, 2, 3, color = (.8, 0.2, .2), )\n",
    "mlab.title('A System of Linear Equations')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <font face=\"gotham\" color=\"purple\"> Visualisation of An Inconsistent System </font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's visualise the linear system that does not have a solution."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{align}\n",
    "x+y+z&=1\\\\\n",
    "x-y-2z&=2\\\\\n",
    "2x-z&=1\n",
    "\\end{align}\n",
    "Rearrange the system to solve for $z$:\n",
    "\n",
    "\\begin{align}\n",
    "z&=1-x-y\\\\\n",
    "z&=\\frac{x}{2}-\\frac{y}{2}+1\\\\\n",
    "z&=2x-1\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mlab.clf()\n",
    "X, Y = np.mgrid[-5:5:21*1j, -5:5:21*1j]\n",
    "Z = 1 - X - Y\n",
    "mlab.mesh(X, Y, Z,colormap=\"spring\")\n",
    "\n",
    "Z = X/2 - Y/2 + 1\n",
    "mlab.mesh(X, Y, Z,colormap=\"summer\")\n",
    "\n",
    "Z = 2*X - 1\n",
    "mlab.mesh(X, Y, Z,colormap=\"autumn\")\n",
    "mlab.axes()\n",
    "mlab.outline()\n",
    "mlab.title('A Inconsistent System of Linear Equations')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <font face=\"gotham\" color=\"purple\"> Visualisation of A System With Infinite Numbers of Solutions </font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Our system of equations is given\n",
    "\n",
    "\\begin{align}\n",
    "y-z=&4\\\\\n",
    "2x+y+2z=&4\\\\\n",
    "2x+2y+z=&8\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Rearrange to solve for $z$\n",
    "\n",
    "\\begin{align}\n",
    "z=&y-4\\\\\n",
    "z=&2-x-\\frac{y}{2}\\\\\n",
    "z=&8-2x-2y\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mlab.clf()\n",
    "X, Y = np.mgrid[-2:2:21*1j, 2:6:21*1j]\n",
    "Z = Y - 4\n",
    "mlab.mesh(X, Y, Z,colormap=\"spring\")\n",
    "\n",
    "Z = 2 - X - Y/2\n",
    "mlab.mesh(X, Y, Z,colormap=\"summer\")\n",
    "\n",
    "Z = 8 - 2*X - 2*Y\n",
    "mlab.mesh(X, Y, Z,colormap=\"autumn\")\n",
    "mlab.axes()\n",
    "mlab.outline()\n",
    "mlab.title('A System of Linear Equations With Infinite Number of Solutions')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The solution of the system is $(x,y,z)=(-3z/2,z+4,z)^T$, where $z$ is a **free variable**. \n",
    "\n",
    "The solution is an infinite line in $\\mathbb{R}^3$, to visualise the solution requires setting a range of $x$ and $y$, for instance we can set\n",
    "\n",
    "\\begin{align}\n",
    "-2 \\leq x \\leq 2\\\\\n",
    "2 \\leq y \\leq 6\n",
    "\\end{align}\n",
    "\n",
    "which means\n",
    "\n",
    "\\begin{align}\n",
    "-2\\leq -\\frac32z\\leq 2\\\\\n",
    "2\\leq z+4 \\leq 6\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can pick one inequality to set the range of $z$, e.g. second inequality: $-2 \\leq z \\leq 2$. \n",
    "\n",
    "Then plot the planes and the solutions together."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mlab.clf()\n",
    "X, Y = np.mgrid[-2:2:21*1j, 2:6:21*1j]\n",
    "Z = Y - 4\n",
    "mlab.mesh(X, Y, Z,colormap=\"spring\")\n",
    "\n",
    "Z = 2 - X - Y/2\n",
    "mlab.mesh(X, Y, Z,colormap=\"summer\")\n",
    "\n",
    "Z = 8 - 2*X - 2*Y\n",
    "mlab.mesh(X, Y, Z,colormap=\"autumn\")\n",
    "\n",
    "ZL = np.linspace(-2, 2, 20) # ZL means Z for line, we have chosen the range [-2, 2]\n",
    "X = -3*ZL/2\n",
    "Y = ZL + 4\n",
    "\n",
    "mlab.plot3d(X, Y, ZL)\n",
    "\n",
    "mlab.axes()\n",
    "mlab.outline()\n",
    "mlab.title('A System of Linear Equations With Infinite Number of Solutions')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <font face=\"gotham\" color=\"purple\"> Reduced Row Echelon Form </font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For easy demonstration, we will be using SymPy frequently in lectures. SymPy is a very power symbolic computation library, we will see its basic features as the lectures move forward."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We define a SymPy matrix:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}5 & 0 & 11 & 3\\\\7 & 23 & -3 & 7\\\\12 & 11 & 3 & -4\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡5   0   11  3 ⎤\n",
       "⎢              ⎥\n",
       "⎢7   23  -3  7 ⎥\n",
       "⎢              ⎥\n",
       "⎣12  11  3   -4⎦"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "M = sy.Matrix([[5, 0, 11, 3], [7, 23, -3, 7], [12, 11, 3, -4]]); M"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Think of it as an **augmented matrix** which combines coefficients of linear system. With row operations, we can solve the system quickly. Let's turn it into a **row reduced echelon form**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left( \\left[\\begin{matrix}1 & 0 & 0 & - \\frac{2165}{1679}\\\\0 & 1 & 0 & \\frac{1358}{1679}\\\\0 & 0 & 1 & \\frac{1442}{1679}\\end{matrix}\\right], \\  \\left( 0, \\  1, \\  2\\right)\\right)$"
      ],
      "text/plain": [
       "⎛⎡         -2165 ⎤           ⎞\n",
       "⎜⎢1  0  0  ──────⎥           ⎟\n",
       "⎜⎢          1679 ⎥           ⎟\n",
       "⎜⎢               ⎥           ⎟\n",
       "⎜⎢          1358 ⎥           ⎟\n",
       "⎜⎢0  1  0   ──── ⎥, (0, 1, 2)⎟\n",
       "⎜⎢          1679 ⎥           ⎟\n",
       "⎜⎢               ⎥           ⎟\n",
       "⎜⎢          1442 ⎥           ⎟\n",
       "⎜⎢0  0  1   ──── ⎥           ⎟\n",
       "⎝⎣          1679 ⎦           ⎠"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "M_rref = M.rref(); M_rref # .rref() is the SymPy method for row reduced echelon form"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Take out the first element in the big parentheses, i.e. the rref matrix."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1, 0, 0, -2165/1679],\n",
       "       [0, 1, 0, 1358/1679],\n",
       "       [0, 0, 1, 1442/1679]], dtype=object)"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "M_rref = np.array(M_rref[0]);M_rref"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you don't like fractions, convert it into float type."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.   ,  0.   ,  0.   , -1.289],\n",
       "       [ 0.   ,  1.   ,  0.   ,  0.809],\n",
       "       [ 0.   ,  0.   ,  1.   ,  0.859]])"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "M_rref.astype(float)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The last column of the rref matrix is the solution of the system."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <font face=\"gotham\" color=\"purple\"> Example: rref and Visualisation </font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's use ```.rref()``` method to compute a solution of a system then visualise it. Consider the system:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{align}\n",
    "3x+6y+2z&=-13\\\\\n",
    "x+2y+z&=-5\\\\\n",
    "-5x-10y-2z&=19\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Extract the augmented matrix into a SymPy matrix:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}3 & 6 & 2 & -13\\\\1 & 2 & 1 & -5\\\\-5 & -10 & -2 & 19\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡3    6   2   -13⎤\n",
       "⎢                ⎥\n",
       "⎢1    2   1   -5 ⎥\n",
       "⎢                ⎥\n",
       "⎣-5  -10  -2  19 ⎦"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = sy.Matrix([[3, 6, 2, -13], [1, 2, 1, -5], [-5, -10, -2, 19]]);A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left( \\left[\\begin{matrix}1 & 2 & 0 & -3\\\\0 & 0 & 1 & -2\\\\0 & 0 & 0 & 0\\end{matrix}\\right], \\  \\left( 0, \\  2\\right)\\right)$"
      ],
      "text/plain": [
       "⎛⎡1  2  0  -3⎤        ⎞\n",
       "⎜⎢           ⎥        ⎟\n",
       "⎜⎢0  0  1  -2⎥, (0, 2)⎟\n",
       "⎜⎢           ⎥        ⎟\n",
       "⎝⎣0  0  0  0 ⎦        ⎠"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A_rref = A.rref(); A_rref"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In case you are wondering what's $(0, 2)$: they are the column number of pivot columns, in the augmented matrix above the pivot columns resides on the $0$th and $2$nd column."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Because it's not a rank matrix, therefore solutions is in general form\n",
    "\\begin{align}\n",
    "x + 2y & = -3\\\\\n",
    "z &= -2\\\\\n",
    "y &= free\n",
    "\\end{align}\n",
    "Let's pick 3 different values of $y$, for instance $(3, 5, 7)$, to calculate $3$ special solutions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ -9,   3,  -2],\n",
       "       [-13,   5,  -2],\n",
       "       [-17,   7,  -2]])"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "point1 = (-2*3-3, 3, -2)\n",
    "point2 = (-2*5-3, 5, -2)\n",
    "point3 = (-2*7-3, 7, -2)\n",
    "special_solution = np.array([point1, point2, point3]); special_solution # each row is a special solution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can visualise the general solution, and the 3 specific solutions altogether."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "y = np.linspace(2, 8, 20) # y is the free variable\n",
    "x = -3 - 2*y\n",
    "z = np.full((len(y), ), -2) # z is a constant"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "4145126158fe4b3d8111e8b470633547",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Canvas(toolbar=Toolbar(toolitems=[('Home', 'Reset original view', 'home', 'home'), ('Back', 'Back to previous …"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure(figsize = (12,9))\n",
    "ax = fig.add_subplot(111, projection='3d')\n",
    "ax.plot(x, y, z, lw = 3, color = 'red')\n",
    "ax.scatter(special_solution[:,0], special_solution[:,1], special_solution[:,2], s = 200)\n",
    "ax.set_title('General Solution and Special Solution of the Linear Sytem', size= 16)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <font face=\"gotham\" color=\"purple\"> Example: A Symbolic Solution </font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider a system where all right-hand side values are indeterminate:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{align}\n",
    "x + 2y - 3z &= a\\\\\n",
    "4x - y + 8z &= b\\\\\n",
    "2x - 6y - 4z &= c\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We define $a, b, c$ as SymPy objects, then extract the augmented matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}1 & 2 & -3 & a\\\\4 & -1 & 8 & b\\\\2 & -6 & -4 & c\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡1  2   -3  a⎤\n",
       "⎢            ⎥\n",
       "⎢4  -1  8   b⎥\n",
       "⎢            ⎥\n",
       "⎣2  -6  -4  c⎦"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a, b, c = sy.symbols('a, b, c', real = True)\n",
    "A = sy.Matrix([[1, 2, -3, a], [4, -1, 8, b], [2, -6, -4, c]]); A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can immediately achieve the symbolic solution by using ```.rref()``` method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left( \\left[\\begin{matrix}1 & 0 & 0 & \\frac{2 a}{7} + \\frac{b}{7} + \\frac{c}{14}\\\\0 & 1 & 0 & \\frac{16 a}{91} + \\frac{b}{91} - \\frac{10 c}{91}\\\\0 & 0 & 1 & - \\frac{11 a}{91} + \\frac{5 b}{91} - \\frac{9 c}{182}\\end{matrix}\\right], \\  \\left( 0, \\  1, \\  2\\right)\\right)$"
      ],
      "text/plain": [
       "⎛⎡            2⋅a   b   c    ⎤           ⎞\n",
       "⎜⎢1  0  0     ─── + ─ + ──   ⎥           ⎟\n",
       "⎜⎢             7    7   14   ⎥           ⎟\n",
       "⎜⎢                           ⎥           ⎟\n",
       "⎜⎢          16⋅a   b    10⋅c ⎥           ⎟\n",
       "⎜⎢0  1  0   ──── + ── - ──── ⎥, (0, 1, 2)⎟\n",
       "⎜⎢           91    91    91  ⎥           ⎟\n",
       "⎜⎢                           ⎥           ⎟\n",
       "⎜⎢           11⋅a   5⋅b   9⋅c⎥           ⎟\n",
       "⎜⎢0  0  1  - ──── + ─── - ───⎥           ⎟\n",
       "⎝⎣            91     91   182⎦           ⎠"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A_rref = A.rref(); A_rref"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Of course, we can substitute values of $a$, $b$ and $c$ to get a specific solution. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}1 & 0 & 0 & \\frac{31}{14}\\\\0 & 1 & 0 & - \\frac{16}{91}\\\\0 & 0 & 1 & - \\frac{69}{182}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡          31 ⎤\n",
       "⎢1  0  0   ── ⎥\n",
       "⎢          14 ⎥\n",
       "⎢             ⎥\n",
       "⎢         -16 ⎥\n",
       "⎢0  1  0  ────⎥\n",
       "⎢          91 ⎥\n",
       "⎢             ⎥\n",
       "⎢         -69 ⎥\n",
       "⎢0  0  1  ────⎥\n",
       "⎣         182 ⎦"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "vDict = {a: 3, b: 6, c: 7}\n",
    "A_rref = A_rref[0].subs(vDict);A_rref # define a dictionary for special values to substitute in"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## <font face=\"gotham\" color=\"purple\"> Example: Polynomials </font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider this question : How to find a cubic polynomial that passes through each of these points $(1,3)$,$(2, -2)$ ,$(3, -5)$, and $(4, 0)$.\n",
    "\n",
    "The form of cubic polynomial is \n",
    "\n",
    "\\begin{align}\n",
    "y=a_0+a_1x+a_2x^2+a_3x^3\n",
    "\\end{align}\n",
    "We substitute all the points:\n",
    "\n",
    "\\begin{align}\n",
    "(x,y)&=(1,3)\\qquad\\longrightarrow\\qquad \\ 2=a_0+3a_1+9a_2 +27a_3 \\\\\n",
    "(x,y)&=(2,-2)\\qquad\\longrightarrow\\qquad 3=a_0+a_1+a_2+a_3\\\\\n",
    "(x,y)&=(3,-5)\\qquad\\longrightarrow\\qquad 2=a_0-4a_1+16a_2-64a_3\\\\\n",
    "(x,y)&=(4,0)\\qquad\\longrightarrow\\qquad -2=a_0+2a_1+4a_2+8a_3\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It turns to be a linear system, the rest should be familiar already."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The augmented matrix is "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}1 & 1 & 1 & 1 & 3\\\\1 & 2 & 4 & 8 & -2\\\\1 & 3 & 9 & 27 & -5\\\\1 & 4 & 16 & 64 & 0\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡1  1  1   1   3 ⎤\n",
       "⎢                ⎥\n",
       "⎢1  2  4   8   -2⎥\n",
       "⎢                ⎥\n",
       "⎢1  3  9   27  -5⎥\n",
       "⎢                ⎥\n",
       "⎣1  4  16  64  0 ⎦"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = sy.Matrix([[1, 1, 1, 1, 3], [1, 2, 4, 8, -2], [1, 3, 9, 27, -5], [1, 4, 16, 64, 0]]); A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left( \\left[\\begin{matrix}1 & 0 & 0 & 0 & 4\\\\0 & 1 & 0 & 0 & 3\\\\0 & 0 & 1 & 0 & -5\\\\0 & 0 & 0 & 1 & 1\\end{matrix}\\right], \\  \\left( 0, \\  1, \\  2, \\  3\\right)\\right)$"
      ],
      "text/plain": [
       "⎛⎡1  0  0  0  4 ⎤              ⎞\n",
       "⎜⎢              ⎥              ⎟\n",
       "⎜⎢0  1  0  0  3 ⎥              ⎟\n",
       "⎜⎢              ⎥, (0, 1, 2, 3)⎟\n",
       "⎜⎢0  0  1  0  -5⎥              ⎟\n",
       "⎜⎢              ⎥              ⎟\n",
       "⎝⎣0  0  0  1  1 ⎦              ⎠"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A_rref = A.rref(); A_rref"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1, 0, 0, 0, 4],\n",
       "       [0, 1, 0, 0, 3],\n",
       "       [0, 0, 1, 0, -5],\n",
       "       [0, 0, 0, 1, 1]], dtype=object)"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A_rref = np.array(A_rref[0]); A_rref"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The last column is the solution, i.e. the coefficients of the cubic polynomial."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 4.,  3., -5.,  1.])"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "poly_coef = A_rref.astype(float)[:,-1]; poly_coef"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Cubic polynomial form is:\n",
    "\\begin{align}\n",
    "y = 4 + 3x - 5x^2 + x^3\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since we have the specific form of the cubic polynomial, we can plot it"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.linspace(-5, 5, 40)\n",
    "y = poly_coef[0] + poly_coef[1]*x + poly_coef[2]*x**2 + poly_coef[3]*x**3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfUAAAHWCAYAAABucBCJAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABFz0lEQVR4nO3deXhU1f3H8fchCVsg7CiCikaJiIKiQgWU4Ar8tCqLIoobFlDcgVZr60prrQItooCCBRRFFi2IgKVI3NCiuIAbCIiCsikEiREIyfn9cWaYhCRkm5ubufN5Pc883HPnzsyXy/LJnHPuucZai4iIiMS+an4XICIiItGhUBcREQkIhbqIiEhAKNRFREQCQqEuIiISEAp1ERGRgIhKqBtjnjXGbDPGfJZvX0NjzGJjzNehXxsU89ruxpjVxpi1xpi7o1GPiIhIPIrWN/UpQPeD9t0NLLHWHg8sCbULMMYkAE8CPYATgSuNMSdGqSYREZG4EpVQt9a+Bew4aPclwNTQ9lTg0iJe2gFYa61db63dB8wIvU5ERETKyMsx9cOstZsBQr82LeKY5sDGfO1NoX0iIiJSRok+f74pYl+R69YaYwYBgwBq1qx52lFHHeVlXXEvLy+PatU0j9JrOs/e0zn2XlDPce3vviNhzx4Afm3enP3Jyb7VsmbNmh+ttU1KOs7LUN9qjGlmrd1sjGkGbCvimE3AkfnaLYAfinoza+3TwNMAaWlpdvXq1dGuV/LJyMggPT3d7zICT+fZezrH3gvkOc7NhZSUSPvDD+Hww30rxxjzbWmO8/JHq3nAtaHta4G5RRzzAXC8MeYYY0x1oF/odSIiIv5ZvRqys932EUf4GuhlEa1L2l4E3gPSjDGbjDEDgb8B5xtjvgbOD7UxxhxhjFkAYK3dD9wCvA58Ccy01n4ejZpERETKbcWKyHb79v7VUUZR6X631l5ZzFPnFnHsD0DPfO0FwIJo1CEiIhIV+UP9tNP8q6OMgjezQUREpKI++iiyrVAXERGJUXl58PHHkXYMdb8r1EVERPL7+mvIynLbhx3mJsrFCIW6iIhIfgdPkjNFLalSNSnURURE8ovRSXKgUBcRESkoRifJgUJdREQkIi+vYKjH0CQ5UKiLiIhErF8PP//sths3hiOPPPTxVYxCXUREJCyGJ8mBQl1ERCQihifJgUJdREQkIoYnyYFCXURExLE2pifJgUJdRETE2bABdu502w0aQMuWflZTLgp1ERERiPlJcqBQFxERcWJ8khwo1EVERJwYnyQHCnURERE3Se7g7vcYpFAXERHZuBF++slt16sHqan+1lNOCnUREZH839JPPTUmJ8mBQl1ERCQQk+RAoS4iIhKISXKgUBcRkXgXkElyoFAXEZF49/33sG2b265TB44/3t96KkChLiIi8S1/1/upp0K12I3G2K1cREQkGgIySQ4U6iIiEu8CMkkOFOoiIhLvAjJJDhTqIiISzzZvdg+A2rUhLc3feipIoS4iIvErf9f7KadAQoJvpUSDQl1EROJXgCbJgUJdRETiWYAmyYFCXURE4lmAJsmBQl1EROLVtm2waZPbrlkTWrf2t54oUKiLiEh8yt/13q4dJCb6V0uUKNRFRCQ+BWySHCjURUQkXgVskhwo1EVEJF4FbJIcKNRFRCQe/fQTfPut265RA9q08beeKFGoi4hI/Mnf9d62LSQl+VdLFCnURUQk/gSw6x0U6iIiEo8COEkOFOoiIhKP9E1dREQkAHbuhPXr3XZSEpx0kr/1RJFCXURE4svHH0e2Tz7ZzX4PCIW6iIjEl4B2vYNCXURE4k0Al4cN8zTUjTFpxphP8j1+NsbccdAx6caYXfmOuc/LmkREJM7ln/kesG/qnt6Sxlq7GjgFwBiTAHwPvFLEoW9bay/yshYRERF27YKvv3bbiYlu4ZkAqczu93OBddbabyvxM0VERCI++SSy3aaNu496gFTmzWP7AS8W89yZxphPgR+A4dbazw8+wBgzCBgE0KRJEzIyMryqU4CsrCyd40qg8+w9nWPvxdI5bjFzJseFtjc3a8bqGKm7tIy11vsPMaY6LrDbWGu3HvRcCpBnrc0yxvQE/mmtPf5Q75eWlmZXr17tXcFCRkYG6enpfpcReDrP3tM59l5MneOrroIXXnDb48bB0KH+1lNKxpgV1trTSzqusrrfewAfHRzoANban621WaHtBUCSMaZxJdUlIiLxJMCT5KDyQv1Kiul6N8Ycbowxoe0OoZp+qqS6REQkXuzeDeFe3mrVoF07f+vxgOdj6saY2sD5wOB8+4YAWGsnAH2Am4wx+4FfgX62MsYEREQkvnz6KYTj5cQToXZtf+vxgOehbq3NBhodtG9Cvu1xwDiv6xARkTgX4JXkwrSinIiIxIcAryQXplAXEZH4EPBJcqBQFxGReLBjB3zxhdtOSIBTTvG1HK8o1EVEJPjefjsySe6006BOHX/r8YhCXUREgi//ynGxslBOOSjURUQk+BTqIiIiAbBjh7tGHdx4eufO/tbjIYW6iIgE28Hj6Skp/tbjIYW6iIgEW5x0vYNCXUREgk6hLiIiEgBxNJ4OCnUREQmyOBpPB4W6iIgEWRx1vYNCXUREgkyhLiIiEgBxNp4OCnUREQmqOBtPB4W6iIgEVZx1vYNCXUREgkqhLiIiEgBxOJ4OCnUREQmiOBxPB4W6iIgEURx2vYNCXUREgkihLiIiEgBxOp4OCnUREQmaOB1PB4W6iIgETZx2vYNCXUREgkahLiIiEgBxPJ4OCnUREQmSOB5PB4W6iIgESRx3vYNCXUREgkShLiIiEgBxPp4OCnUREQmKOB9PB4W6iIgERZx3vYNCXUREguLNNyPbCnUREZEYtXMnfPKJ247T8XRQqIuISBBoPB1QqIuISBBoPB1QqIuISBAo1AGFuoiIxDqNpx+gUBcRkdim8fQDFOoiIhLb1PV+gEJdRERim0L9AIW6iIjELo2nF6BQFxGR2KXx9AIU6iIiErvU9V6A56FujNlgjFlljPnEGPNhEc8bY8xYY8xaY8xKY0x7r2sSEZGAUKgXkFhJn9PNWvtjMc/1AI4PPToC40O/ioiIFE/j6YVUhe73S4Bp1nkfqG+MaeZ3USIiUsVpPL2Qygh1C/zHGLPCGDOoiOebAxvztTeF9omIiBRPXe+FVEb3e2dr7Q/GmKbAYmPMV9bat/I9b4p4jT14R+gHgkEATZo0ISP/H6ZEXVZWls5xJdB59p7Osff8OsenzZ9P3dD2yoYN2aE/Z4y1hfLTuw8z5gEgy1r7eL59E4EMa+2LofZqIN1au7m490lLS7OrV6/2uty4lpGRQbp+8vWczrP3dI6958s53rkTGjVy3e8JCbBjR6C7340xK6y1p5d0nKfd78aYZGNM3fA2cAHw2UGHzQOuCc2C/w2w61CBLiIiovH0onnd/X4Y8IoxJvxZL1hrFxljhgBYaycAC4CewFogG7je45pERCTWaTy9SJ6GurV2PdCuiP0T8m1bYKiXdYiISMAo1ItUFS5pExERKT1dn14shbqIiMQWjacXS6EuIiKxRV3vxVKoi4hIbFGoF6uy1n4PDmthyxb46ivIynJjOQ0b+l2ViEh80Hj6ISnUi7N/P6xf78L7yy8L/rprV+S4xETo1g169YJLL4XDD/etZBGRwNN4+iEp1LOyYPXqSGCHw/vrryEnp+TX798Pixe7x803u58ae/d2IX/UUd7XLyIST9T1fkjxEerWwtatRX/r3rix5NcfLCUFTjgB8vLgw3y3iLcW3nnHPe68E04/3YV7797QqlX0fj8iIvFKoX5IwQ11a+Hpp2HqVBfgmZllf4/mzV14t25d8NdmzcCE7kOzcSO88grMmVOwWwhc4H/4Ifzxj9CmTeQbfNu2kdeLiEjpaDy9RMEM9e3b4brrYMGCko9NTITjjisc3mlppRurOfJIuO0299i6FebOhZdfhiVLXNd82Oefu8dDD0FqaiTgzzgDqukiBBGREuX/4tS+vcbTixC8UM/IgKuugh9+KLi/bl0X2AeHd2oqJCVF57MPOwwGDXKPnTth/nz3Df7112HPnshx69bB3//uHi1awGWXuZDv0sX99CkiIoWp671EwQn1/fvh4YfdI38X+F13uccRR1Rul3eDBjBggHtkZcGiRS7g58937bBNm+CJJ9yjSRM3g75XLzjnHKhevfLqFRGp6hTqJQpGv+/GjS4EH3ooEuiNG8Nrr8GoUW5s3M8x7Dp1oE8fePFFNzTw6qtw/fWFr2/fvh2eeQZ69HDj9s88U/AHFBGReJV/PL1aNdezKYXEfqjPmwennOLGWsK6dYNPP4WePX0rq1g1a8JFF8Gzz7pFbP77X7jppsLXt+/Y4brxe/WCH3/0p1YRkapC16eXSuyG+t69cPvtcMklLgDB/fT28MPumvEjjvC3vtJISoJzz4WnnoLvv3eXwt11lxtnD/v3v+Hkk924vIhIvFLXe6nEZKhX27cPzjwTxo6N7GzRwv2h/+lPsTnZrFo1d3nGqFFuMZyh+W4xv2ULdO/ufoj59Vf/ahQR8YtCvVRiMtSTv/0WPv44suOSS1x3+1ln+VdUNNWuDePGuTkBTZtG9o8d6y6BW7myUsux1tKuXTumTp0KwNq1axk8eDDt2rUjISGB9HL+A5s9ezadOnWiUaNG1KxZk7S0NEaOHMm+ffsOHDN06FAGDhwYjd+GiMQqjaeXWkyG+oFxlerV3azxV14J5k1VevaEVavcGHzY55+7YB892q1oVwlmzpzJzp076d+/f6iEz1mwYAGtWrWiVQVWyvvpp5/o1q0bkyZNYuHChdxwww385S9/4a677jpwzIgRI5g+fTpr166t8O9DRGKUxtNLLXYvaWvVCl56yU2SC7KmTd1kwIkT3Xj7r7/Cvn0wbJhbXGfqVDe730Njx45lwIABJIWu57/44ou55JJLAOjTpw8/lnMi3+DBgwu0u3Xrxs8//8yTTz7JE088gTGGli1b0qVLF8aPH8+oUaMq9hsRkdikrvdSi8lv6jkpKbBiRfADPcwYGDIEPvrIraIUtmSJm0Q3e7ZnH7127VqWLVtGnz59Duyr5uEKeI0aNSrQ/Q7Qu3dvpk+fTl4l9UyISBWjUC+1mAz1PYcf7q79jjcnnADvvQf33BO57n7nTujb1133vnt31D9yyZIlJCcn065du6i/d1hubi7Z2dm88847jB07lptuugmTb12BTp06sXXrVlatWuVZDSJSRWk8vUxiMtTjWvXq8Ne/up9c89/adcoU13Px3ntR/bgVK1bQunVrT7+dJycnk5yczFlnnUXXrl157LHHCjzfpk0bEhISWL58uWc1iEgVpfH0MlGox6qzz3Yz/kOT1wBYv95dAfDAAwVvJlMBW7ZsoXHjxlF5r+IsW7aMt99+m1GjRjF37lxuueWWAs8nJiZSv359tmzZ4mkdIlIFqeu9TBTqsax+fZg+3T3q1XP7cnPhwQdduK9bV+GP2LNnDzVq1Kjw+xxK+/bt6dKlC3fddRdjx45l/PjxrDuo9ho1arAn/01xRCQ+KNTLRKEeBP37u2/tZ58d2ff++647/l//qtD68Q0bNiSzPPeiL6f2oYmA33zzTYH9mZmZNAziZYsiUjyNp5eZQj0ojj4a3ngDHnnE3SMe3N3gbrjBTaT76adyvW1aWlqhgPXSu+++C8AxxxxzYN/27dvJzs6u0DXxIhKDNJ5eZgr1IElIgLvvdt/S09Ii++fMgbZt3c1jyqhz58589913bN++/cC+7OxsZs+ezezZs/n+++/Zvn37gXZ2dvaB4zIyMjDGkJG/+yyf7t278/jjj7Nw4UL+85//cP/99zNs2DCuuOIKUlNTDxz34YcfYoyhU6dOZa5fRGKYut7LLHYXn5HinXaau6Z92DCYMMHt++EHOP98uPNON3u+Zs1SvVV6ejoNGzZk0aJFDBgwAIBt27bRt2/fAseF29988w0tW7YEOBDwTfMvdZvPGWecwZQpU9iwYQOJiYkce+yxPPLIIwwZMqTAcYsWLaJr1640atSoVDWLSEAo1MtM39SDqnZtGD/e3bu9SZPI/jFjoGNHKGWXevXq1bn66quZMWPGgX0tW7bEWlvkIxzoAP/73/9IT0/nxBNPLPK9H374YT777DOysrLIzMzko48+4tZbbz2wch24a9jnzJmj9d9F4o3G08tFoR50F13k1o/Pf2/5lSuha1d3CVwpjBgxgoyMDNasWVOmj162bFmBddzLY9asWdSqVYt+/fpV6H1EJMZoPL1cFOrx4LDDYP58d9/28OVpGze67qxSXPbWokULJk+ezObNm8v0sYsXL+biiy8uR8ER1lomT55MYqJGikTiirrey0X/U8YLY+Cmm+C44+C3v4U9eyLBvnSp238Ifn1TvvLKK335XBHxmUK9XPRNPd6cf74bZw9PlNu0yf2D0a1NRaSq0Hh6uSnU49F557nu+Fq1XPv77yE9nXVLNnDzzW7o6pxzupKSAjffHJWF6URESk/j6eWmUI9X555bINgXfn8ybc9vyqRnLLt3g7WG3bth0iR3ifvChT7XKyLxQ13v5aZQj2fnnAPz57Ouxon0YTbZtjY5+02BQ3JyIDsb+vTRN3YRqSQK9XJTqMe7c85h1IX/IYekQx6Wk+MucRcR8ZTG0ytEoS48v7Q5OVQ/5DE5OfDcc5VUkIjEL42nV4hCXcjKiu5xIiLlpq73ClGoC3XqRPc4EZFyy3/jKYV6mSnUhauvhqRDD6mTlGgJ3c9FRMQbX3/tlrUGt/rlWWf5W08MUqgLw4aVItT3/8qd/1e2td9FRMpk9uzIdvfuULeuf7XEKIW6kJrq/i3Vrl043JPYR21+YTa9Sb22C3z2mT9Fikjw5Q/1Pn38qyOGKdQFgB493M3bBg1yk02NsaSkwKBeP7EyuRM9WATbt7tr2xXsIhJt69fDRx+57erVoYI3g4pXCnU5IDUVxo2DXbvgjTfeZNcuGDenGamLJ0S6wbZvh27dIuNeIiLRMGdOZPuCC6BePf9qiWGehrox5khjzFJjzJfGmM+NMbcXcUy6MWaXMeaT0OM+L2uScjjzTPjPfyLXi/74owv2lSv9rUtEgkNd71Hh9Tf1/cAwa21r4DfAUGPMiUUc97a19pTQ4yGPa5Ly+M1vCgb7Tz+5rvhPP/W3LhGJfd9+C8uXu+3ERHd7aCkXT0PdWrvZWvtRaHs38CXQ3MvPFA917Fg42M89N7Kko4hIeeTvej/vPGjQwL9aYlyljakbY1oCpwL/K+LpM40xnxpjFhpj2lRWTVIOHTvC4sWR8a5wsH/8sb91iUjsyt/13revf3UEgLHhNXa9/BBj6gBvAn+x1r580HMpQJ61NssY0xP4p7X2+CLeYxAwCKBJkyanzZw50/O641lWVhZ1DrGEXN2vvqLd8OEk/vILADl16/Lp44+T1apVZZUYCCWdZ6k4nWPvVeQc19i+nTMvvxwAW60a7778Mvs1Sa6Qbt26rbDWnl7ScZ6HujEmCZgPvG6tHV2K4zcAp1trfyzumLS0NLt69eroFSmFZGRkkF7SEo0ffgjnnw+Zma7doIFb4rF9e6/LC4xSnWepEJ1j71XoHI8dC7eH5lCff74b4pNCjDGlCnWvZ78bYDLwZXGBbow5PHQcxpgOoZp+8rIuiZLTT3chHh7/2rnTdcXrcjcRKS3Neo8qr8fUOwMDgHPyXbLW0xgzxBgzJHRMH+AzY8ynwFign62MMQGJjtNOKxjsmZlu0YitW30tS0RiwObN8M47brtaNbj0Ul/LCYJEL9/cWvsOYEo4Zhwwzss6xGPt27tg79rV3Z/122/hssvgjTegZk2/qxORqurllyP3Tu/aFZo29beeANCKchId7dvDjBnup22A996DgQMj/2BFRA6mWe9Rp1CX6Pm//4NRoyLtF16AkSP9q0dEqq6tW+Gtt9y2Ma53TypMoS7RdfvtMHhwpH3ffaDLD0XkYP/+N+Tlue2zzoLDD/e1nKBQqEt0GQNPPOFmwYdde21kCUgREYBZsyLbmvUeNQp1ib6kJPcPNrwQzZ49bi3n777zty4RqRq2b4eMjEi7Vy/fSgkahbp4o0EDmD8fGjZ07a1b3aVuu3f7W5eI+G/uXMjNddudO0Nz3RIkWhTq4p3jj3c3akgMXTm5ciVcdVXkH7OIxCctOOMZhbp4Kz0dJkyItF99Fe6+27dyRMRnO3bAkiWRtrreo0qhLt4bOBBGjIi0H38cJk3yrx4R8c/cubB/v9vu2BGOOsrfegJGoS6V45FH3GS5sJtugqVL/atHRPyhrndPKdSlciQkwPTp0K6da+/fD717w9df+1uXiFSezExYvDjSVqhHnUJdKk+dOm5MPbzIxM6dcNFFboxNRILv1VchJ8dtn346tGzpazlBpFCXynXkkTBvXuRGL2vWuDWfw//QRSS41PXuOYW6VL4zzoBp0yLtN96AW27RzV9Eguznn+H11yPt3r39qyXAFOrij7594eGHI+2nn4Z//MO3ckTEY/Pnw969bvuUU+C443wtJ6gU6uKfe+91i9GEDRvm/uGLSPCo671SKNTFP8a469U7dXJta+HKK93KcyISHFlZsHBhpK17p3tGoS7+qlkTXnklMgs2K8utEb9li69liUgULVjgbuwEcPLJkZs9SdQp1MV/TZu6S13q1nXt776DSy+FX3/1tSwRiRLdZrXSKNSlajjpJHjpJagW+iv5v//BDTdoRrxIrPvlF/dNPUyh7imFulQdPXrAmDGR9owZ8NBD/tUjIhW3aBFkZ7vt1q3hxBP9rSfgFOpStdx6q1sXPuyBB+DFF30rR0QqKP+sd02Q85xCXaoWY+Cf/4Tzz4/su/56eP99/2oSkfL59deCl6mq691zCnWpepKSYOZMOOEE1967Fy65BL791t+6RKRsXn/dXdECbsb7SSf5W08cUKhL1VS/vvsJv1Ej1962zV3qFv4PQkSqvoMXnDHGv1rihEJdqq7UVHj5ZffNHWDVKhg4UDPiRWLB3r3u5k1h6nqvFAp1qdrOPhsmToy0Z850Y+4iUrUtXgy7d7vtY491672L5xTqUvVdfz0MGRJpDx8Ob7/tXz0iUrKDZ72r671SKNQlNvzjH9Chg9vOzYXLL4fNm30tSUSKsW8fzJ0baavrvdIo1CU21KjhfvJv3Ni1t2xxwZ6T429dIlLYkiWQmem2jz4aTjvN13LiiUJdYseRR7pV5sJLyb7zDvzhD/7WJCKFada7bxTqElvOPRdGjoy0x4xxa8aLSNWQkwP//nekra73SqVQl9jzhz+4xWjCBg6EL77wrx4RicjIgB073PaRR0LHjr6WE28U6hJ7qlWDqVPhuONc+5dfoFcv+Plnf+sSkYK3We3dW13vlUyhLrGpXj23ME2tWq69erW79E0L04j4Z/9+eOWVSFtd75VOoS6x6+ST4ZlnIu2XX4ZRo/yrRyTevfUW/Pij227WDM4809964pBCXWLbVVfBLbdE2n/4gxvTE5HKl3/We+/ekStVpNLojEvsGzUq8o0gLw+uuAK+/97fmkTiTW6u6y0LU9e7LxTqEvuqV3eTc5o2de1t29yylPv2+VuXSDx55x3YutVtH3YYdOlS7KHWWtq1a8fUqVMLPff9999Tp04djDFkRfGujFOmTMEYU+gxYcKEAscNHTqUgQMHRu1zK1ui3wWIREXz5u569fPOc98Y3nvPrRE/dqzflYnEh/xd7716QUJCsYfOnDmTnTt30r9//0LPjRgxgjp16vDLL794USVvvPEGtcITbIFjjz220OefcMIJ3HPPPRwXvsImhuibugRHejo88kik/cQTMH26b+WIxI28PJgzJ9Iuoet97NixDBgwgKTwbZVD3n77bRYtWsTw4cO9qBKAM844g9/85jcHHk3DPXwhLVu2pEuXLowfP96zGrykUJdgGT7cfUsIGzTI3YddRDyT8vnnkRssNW7sbplcjLVr17Js2TL6HBT8ubm53Hrrrdx33300Dt/jwSe9e/dm+vTp5OXl+VpHeSjUJViMgX/9C9LSXDs724X8rl3+1iUSYE3eeivSuOwySCx+ZHfJkiUkJyfTrl27AvsnTJjAnj17GDp0qFdlApCamkpiYiJpaWlMnDixyGM6derE1q1bWRWDXwgU6hI8KSluFm5ysmuvXQvXXuu6CEUkuvLyaPLmm5F2376HPHzFihW0bt2aavkud/vpp5/485//zOjRowt1yUdLs2bNePjhh3nuued49dVX6dixI0OGDGHMmDGFjm3Tpg0JCQksX77ck1q8pFCXYDrxRJg8OdKeOxf+/nf/6hEJmHXr4OabIaVuHrW3byGFTG6uMYl1R6Yf8nVbtmwp1L1+77330rFjR3r27OlZvRdeeCF/+tOfuOCCC+jRowfTpk3j8ssvZ+TIkYW62RMTE6lfvz5btmzxrB6vKNQluK64Au64I9K+9153n2cRqZCFC6FtW5g0CXZnJ2Kpxm7qMWnftbQ9LYmFC4t/7Z49e6hRo8aB9ueff86zzz7L/fffT2ZmJpmZmWRnZwOwa9cufv31V89+H3369GHHjh1s2LCh0HM1atRgz549nn22VzwPdWNMd2PMamPMWmPM3UU8b4wxY0PPrzTGtPe6Jokjf/975HrZvDzo1w82bvS3JpEYtm6dm9yene3usppfjk0kO9s9v25d0a9v2LAhmZmZB9pff/01OTk5nHnmmTRo0IAGDRocGFdv0aIFt956q0e/kwhTxE1nMjMzadiwoeefHW2eXqdujEkAngTOBzYBHxhj5llr898nswdwfOjRERgf+lWk4pKSYOZMaN8etmxx61L36ePWqM73bUFESmfUqMJhfrCcHBgzBsaNK/xcWloa77333oF2ly5dWLp0aYFjFi1axKOPPsqCBQsKXUceTXPmzKFx48YcffTRBfZv376d7OxsWrVq5dlne8Xrb+odgLXW2vXW2n3ADOCSg465BJhmnfeB+saYZh7XJfGkWTMX7OHFMJYvhzvv9LcmkRj1/POlC/Xnniv6uc6dO/Pdd9+xfft2ABo3bkx6enqBxwknnADAWWedRVr4ShYgIyMDYwwZxdzfYdq0aSQmJvLtt98Weq537948+uijLFy4kPnz5zNgwABeeukl7rvvvgKT9gA+/PBDjDF06tTp0L/RKsjrFeWaA/n7OjdR+Ft4Ucc0BzbnP8gYMwgYBNCkSZNi/1AlOrKysgJ3jlsMHsxxTz3lGuPH82W9emy98EJfawriea5qdI6jKyurK1DyPdJ377ZkZLxZaL8xhpSUFMaMGcMFF1xQ5Gu/+uorwC1Gk3/1t/fffx+ADRs2FPln+sUXX5Cbm8t7773HN998U+C5mjVr8uSTT7Jt2zastRx99NHcc889nHzyyYXea9KkSbRr1y4mL2nDWuvZA+gLTMrXHgA8cdAxrwFd8rWXAKcd6n1btWplxVtLly71u4Toy8uz9vLLrXV3Xbe2Zk1rP/7Y15ICeZ6rGJ3j6KpbN/JP6FCPlJTi3+O2226zPXv2LPNn33fffTY9Pb0C1Zds//79tnnz5va5557z9HPKCvjQliJ3ve5+3wQcma/dAvihHMeIVJwx7jK31q1de88ed3vInTv9rUskhlx9tZuqcihJSTBgQPHPjxgxgoyMDNasWVOmz162bBl33XVXmV5TVrNmzaJWrVr069fP08/xiteh/gFwvDHmGGNMdaAfMO+gY+YB14Rmwf8G2GWt3XzwG4lERZ06bmGaOnVce/1697+PFqYRKZVhwyApyR7ymKSkQ09badGiBZMnT2bz5rL9V7948WIuvvjiMr2mrKy1TJ48mcRDrIpXlXka6tba/cAtwOvAl8BMa+3nxpghxpghocMWAOuBtcAzwM1e1iTCCSfAlCmR9muvwV//6ls5IrEkNRVmP7aB2vxCEgVvb5yUBLVruxu2paYe+n369etH165dPay0fK688krOPsTa9VWd59epW2sXWGtbWWtTrbV/Ce2bYK2dENq21tqhoedPttZ+6HVNIvTu7W7+EnbfffD66/7VIxJDemz5FytpyyCeJiUpG2MsKSnu/kkrV0KPHn5XGL+0opzEr0cecbdrBTe3p39/KGJlKRHJZ+9emDSJVNYzjlvZ9cIC3njjTXbtctell/QNXbylUJf4lZgIM2bAEUe49o4dbmGaGFwaUqTSTJ8euc1qs2bg8Ri3lI1CXeLbYYfBrFmRW0WuWAG33eZvTSJVVV4ePP54pH3HHVqZsYpRqIt06gSjR0fazzwDzz7rXz0iVdVrr8GXX7rtunVh8GB/65FCFOoiALfc4sbUw26+GT76yL96RKqi/LcvHjwY6tXzrxYpkkJdBNzCNE8/DSed5Np797oZ8jt2+FuXSFXx3nvwzjtuOykJbr/d33qkSAp1kbDkZLcwTUqKa2/YAFddpYVpRAAeeyyy3b8/tGjhXy1SLIW6SH7HHw/TpkXaixbBQw/5V49IVbBmDfz735F2/jUepEpRqIsc7JJL4O67I+0HH4QFC/yrR8Rvo0a5tRwA/u//IsNUUuUo1EWK8vDDcO65kfZVV7l14kXizdatMHVqpD1ihH+1SIkU6iJFSUyEF1+MjBtmZrqJc7/+6mtZIpXuiSfcxFGADh0ghtdFjwcKdZHiNGni7kwRvs/kJ5+4S93soe9QJRIYWVnw5JOR9u9/764UkSpLoS5yKB07wj//GWlPmeIWpxGJB5MmuV4qgOOOg0sv9bMaKQWFukhJhgyBa66JtG+9FT74wL96RCpDTg6MGRNpDxsGCQn+1SOlolAXKYkxMH48tG3r2vv2ufH1H3/0ty4RL82cCd9957abNIFrr/W3HikVhbpIadSu7RamCS+LuXGjW4AjN9ffukS8YG3BJWFvvRVq1fKvHik1hbpIaaWmwvPPR9qLF8P99/tXj4hXFi+GlSvddu3aboKoxASFukhZXHQR/OlPkfZf/gKvvupfPSJeyP8t/cYboVEj/2qRMlGoi5TVAw/ABRdE2gMGwNq1vpUjElUffQRLlrjthAS4805/65EyUaiLlFVCAkyfDkcd5dq7drmJc9nZ/tYlEg35b9xy+eXQsqVvpUjZKdRFyqNxY7cwTfXqrr1ypbu/tBamkVj2zTdu1nuYloSNOQp1kfI64wwYNy7Sfv55d+mbSKwaPTpyq+HzzoNTT/W3HikzhbpIRdx4I1x/faR9xx3w/vu+lSNSbj/+CJMnR9q//71/tUi5KdRFKsIYtzZ2+BtNTg706QPbtvlbl0hZPfVU5IZFp5zivqlLzFGoi1RUrVowZw40aODa338P/frB/v3+1iVSWtnZ7m5sYSNG6MYtMUqhLhINxxzjZsSH/yNcurTg9ewiVdnUqZFlj48+Gvr29bceKTeFuki09OhRcIW5Rx+FV17xrx6R0sjNhccfj7Tvuityu2GJOQp1kWj6859duIddey2sWeNfPSIlefllWL/ebTdoADfc4G89UiEKdZFoqlbNXdoWXrBj927o1Qt+/tnXskSKZG3BxWaGDoU6dfyrRypMoS4SbQ0buolzNWq49uefw5VX6o5uUvW8+SZ88IHbrlHD3Y1NYppCXcQL7dvD009H2gsWwPDh/tUjUpT8N2657jpo2tS3UiQ6FOoiXrnmGrjnnkj7H/+AiRN9K0ekgFWrYOFCt20MDBvmbz0SFQp1ES+NHOnG1MOGDoX//te/ekTC8s9479ULjj/ev1okahTqIl6qVg2mTXPd8QC5uay7bDg3X7mTlBQ455yupKTAzTfDunX+lipxZONGeOGFSFs3bgkMhbqI15KTYd48OOIIFtKdtlnvMmlGMrt3g7WG3bth0iRo2zbSGyriqX/+M7Li4dlnQ8eO/tYjUaNQF6kMzZuz7qlF9GE22SSTQ/UCT+fkuJU6+/TRN3bxWGZmwbkdunFLoCjURSrJqNdPJieh1iGPycmBMWMqqSCJTxMmQFaW2z7xxIKLJUnMU6iLVJLnn4ec3EP/k8vJgeeeq6SCJP7s3eu63sNGjHDzPiQw9KcpgWStpV27dkydOhWAWbNm8dvf/pbmzZtTp04dTjvtNF588cVKrSn85ciZBfwWaA7UAU4DXix03NChQxk4cGBllShB9/zzsGWL2z7iCOjf3996JOoU6hJIM2fOZOfOnfQP/ac1evRo6tSpw5gxY5g3bx7dunWjf//+PJH/dpMeK7j65mhcmI8B5gHdgP7AEwWOGzFiBNOnT2ft2rWVVqcEVF5ewSVh77gDqlcv9nCJTYl+FyDihbFjxzJgwACSQnebevXVV2ncuPGB58855xx++OEHRo8eza2VtDTm1Ve7We45OQCvAo3zPXsO8AMwmgEDIvW0bNmSLl26MH78eEaNGlUpdUpAzZ8Pq1e77ZQUGDTI33rEE/qmLoGzdu1ali1bRp8+fQ7syx/oYaeeeirbtm2rtLqGDct/R8vC9cCpwDbuHPJrgb29e/dm+vTp5OXleVyhBFr+JWGHDIF69fyrRTyjUJfAWbJkCcnJybRr1+6Qxy1btowTTzyxkqqC1FSYPRtq1y58u+ok9pHAWxxPIqkPXuO6SkM6derE1q1bWbVqVaXVKgGzbBm8+67bTkqC22/3tx7xjEJdAmfFihW0bt2aaoeY1btkyRLmzp3L0KFDK7Eyd/XQypWu5zMlBYyxpKRAj5Onksd8/sjPLvnvv//Aa9q0aUNCQgLLly+v1FolQPKPpV99tZskJ4GkMXUJnC1bthTZ3R62YcMG+vfvzyWXXMJ1111X7s/ZtWsXmzdvLvG4E044oUA7NRXGjXOPjIw3admyJR07/olLjj2W69avdweNHAlpaXD11SQmJlK/fn22hGcti5TF6tUwd26krbsFBppnoW6MeQy4GNgHrAOut9ZmFnHcBmA3kAvst9ae7lVNEh/27NlD7dq1i3xux44d9OjRg6OOOornn3++Qp8za9Ysfve735V4nLW22Od+/vnnSD3//S/06weLFrknBw6EY46Bzp2pUaMGe/bsqVC9EqdGjYLw38GLLnILzkhgedn9vhg4yVrbFlgD3HOIY7tZa09RoEs0NGzYkMzMzEL7s7Ozueiii9i3bx+vvfYaycnJFfqcG2+8EWttiY/iZGdn88c//jFST716MGNG5D/dffvgssvgm2/IzMykYcOGFapX4tCWLRBaqwHQkrBxwLNQt9b+x1obumMA7wMtvPoskfzS0tL45ptvCuzbv38/ffv25euvv2bhwoU0bdrUp+oK1rNp06aC9dSr5y49Cg8fbN/O9u7dyc7OplWrVv4VLLFp7Fj3wyHAb34DXbr4W494rrImyt0AFHf/KQv8xxizwhijCyelwjp37sx3333H9u3bD+y7+eabWbBgAX/+85/ZsWMH77///oHH3r17DxyXkZGBMYaMjAxPawzXM2DAgML1HHEE/PvfBxYG+XDNGgzQqUMHT2uSgNm9G8aPj7RHjABj/KtHKkWFxtSNMf8FDi/iqXuttXNDx9wL7AemF/M2na21PxhjmgKLjTFfWWvfKuKzBgGDAJo0aeL5f7rxLisrK2bPsTGGlJQUxowZwwUXXADAvHnzALi9iEt5XnzxRQ4/3P01fv/99wE3mc7L33+4nnHjxjFu3Lgi6zls+HBa//WvLAK6Ar/efDMZlbRQTpDE8t/lijjq+ec5NjQMld2iBcvr1QOPzkO8nuMqqTRjguV9ANcC7wG1S3n8A8Dwko5r1aqVFW8tXbrU7xIq5LbbbrM9e/Ys8+vuu+8+m56e7kFFRSvpPO+/5x7bHOxzbqqTtU8+WTmFBUis/10ul48+sjYpyf2dAWsnTvT04+LyHFcy4ENbihz1rPvdGNMd+APwW2ttdjHHJBtj6oa3gQuAz7yqSeLHiBEjyMjIYM2aNWV63bJly7jrrrs8qqrsZp10ErXq1KFfeMdtt8HixX6WJFVddra7UYtbjxg6dIDrr/e3Jqk0Xo6pjwPq4rrUPzHGTAAwxhxhjFkQOuYw4B1jzKfAcuA1a+0iD2uSONGiRQsmT55cquvI81u8eDEXX3yxR1WVnTWGyXPmkHh66MKQ3Fzo2xe+/NLfwqTqGj4cvvrKbScnw/TphZcwlMDy7Dp1a+1xxez/AegZ2l4PHHotT5Fy6tevX8kHVXFXXnml2zj5ZPeNa9Mm2LXLXW/8v/9FZsmLALz6asHJcWPHwnFF/lcsAaVlYkViQbNmMG+eWzgeYP166NUL8s3clzi3ZQvccEOk3bu3ut3jkEJdJFaceqrrSg1flvT223DjjQVu/iJxKi8PrrsOfvzRtZs3h6ef1iVscUihLhJLLr0UHn000n7+ebecbG6ubyVJFTBuHLz+uts2BqZNA61AGJcU6iKxZvhwd5u3sClTFOzxbNWqgsu/Dh8O55zjXz3iK4W6SKwxxk2Gyj9+OnWqayvY48uePe7ytfDcilNPhYcf9rcm8ZVCXSQWVasGzzzjxtTDpk1zE6MU7PHj7rvhs9DSHrVquTkXNWr4W5P4SqEuEquqVYOJEyH/7V+fe85NmFKwB9+iRfDPf0bao0dD69b+1SNVgkJdJJZVqwYTJhQcY3/+ebj2WgV7kG3f7n54C7v4Yhg82LdypOpQqIvEumrV3Bh7/v/Up0+Ha66B/fuLf53EJmvdxMitW137sMNg8mRdviaAQl0kGKpVg6eegptuiux74QUYMEDBHjQTJ7qV48KmTIEmTXwrR6oWhbpIUFSrBk8+WTDYZ8xQsAfJV19B/hsO3XYbdO/uXz1S5SjURYLEGBfsQ4dG9s2YAVddpWCPdfv2ucvXfv3VtU86qeBCRCIo1EWCxxh44gm45ZbIvpkzC96OU2LPn/8MH3/stmvUcMMrNWv6W5NUOQp1kSAyxt2h69ZbI/tmzVKwx6o33oDHHou0H33U3blP5CAKdZGgMsZdx3z77ZF9s2fDlVcq2GPJjh3uSgZrXfvCCwv+sCaSj0JdJMiMgTFj4I47IvvmzIF+/RTsscBatwbB99+7duPG8K9/uUmRIkXQ3wyRoDPGrTZ2552RfS+/DFdc4SZf+cxaS7t27Zg6dSoAs2fPplOnTjRq1IiaNWuSlpbGyJEj2edTrVOmTMEYU+gxYcKEqLxu6NChDBw4sLg3cT+EhU2eDM2aVfB3JEGW6HcBIlIJjIFRoyIBD/DKK3D55W4SXfXqvpU2c+ZMdu7cSf/+/QH46aef6NatGyNGjKB+/fosX76cBx54gC1btjBu3Djf6nzjjTeoVavWgfaxxx4bldeNGDGCE044gXvuuYfjjjsu8sTatQW72YcMgd/+tnzFS9xQqIvEC2Pg8ccjAQ8wd67vwT527FgGDBhAUlISAIMPWu60W7du/Pzzzzz55JM88cQTGJ9WTjvjjDOoU6dO1F/XsmVLunTpwvjx4xkV/nPJyXGXIf7yi2unpUX+zEQOQd3vIvHEGDeLesSIyL65c6FPn8jtOyvR2rVrWbZsGX369DnkcY0aNfKt+70y9O7dm+nTp5OXl+d2PPQQLF/utpOS3OVrtWv7V6DEDIW6SLwxxl0S9fvfR/a9+qovwb5kyRKSk5Np165doedyc3PJzs7mnXfeYezYsdx0002+fUsHSE1NJTExkbS0NCZOnBjV13Xq1ImtW7eyatUqePtt+OtfI0+OHAnt21e0fIkT6n4XiUfGwN/+Fgl4gPnzoXdvNzGrku7JvWLFClq3bk21ImZzJycnszf0Q8Y111zDY/mv065EzZo14+GHH6ZDhw7k5uby4osvMmTIELKzs7kz/+TDCryuTZs2JCQksDwjg3ajR0P4G3u3bjB8uJe/PQkYhbpIvDIGHnnEXR71yCNu32uvQa9eLtgrYbWyLVu20Lhx4yKfW7ZsGdnZ2SxfvpyHHnqIW265haeeeqpcn7Nr1y42b95c4nEnnHBCoX0XXnghF1544YF2jx492Lt3LyNHjuT2228v8geSsr4uMTGR+vXrs2XyZPjuO7ezQQOYOlWXr0mZKNRF4pkx8Je/uF/DXb4LFrhgf/llz4N9z5491C5mrLh9qMu5S5cuNG7cmGuvvZZhw4aRmppa5s+ZNWsWv/vd70o8zoYXeClBnz59mDlzJhs2bCj1LPiSXlcjL489q1ZFdkycCEceWer3FgGNqYuIMW7c9t57I/sWLoRzz4VvvvH0oxs2bEhmZmaJx4UD/pty1nPjjTdirS3xUVblHeMv9LoNG8jcuZOG4fb110PfvuV6b4lvCnURccH+8MPupiFhy5ZBu3bw3HORJUqjLC0trVRB/e677wJwzDHHeFJHWc2ZM4fGjRtz9NFHV/x1+/ez/YoryAZaAaSmuuV9RcpB3e8i4hgDDz4IdevCPfdAbi7s3u3WHX/tNRg/3o3zRlHnzp156KGH2L59O02aNAGge/funHfeeQcmj7377ruMGjWKK664okDXe0ZGBt26dWPp0qWkp6dHta78evfuTYcOHWjbti25ubm89NJLvPTSS4wdO7bAuPi0adO44YYbWLduHUcffXSpX8ff/saHy5djgE7VqsH06e7PQKQcFOoiEmGMu4b97LPh6qvdqmYAL73kvrlPmwZRDND09HQaNmzIokWLGDBgAOAWa5kyZQobNmwgMTGRY489lkceeYQhQ4YUeG12djYATZs2jVo9RUlLS+PZZ59l48aNWGs58cQTmTZt2oF6w/Ly8sjNzT3QjV+q102fDg88wCKgK9DowQehY0dPfz8ScKUZZ6pqj1atWlnx1tKlS/0uIS5U6fO8e7e1Awda6zrf3cMYa3//e2v37o3ax9x22222Z8+eZX7dfffdZ9PT00s8rkqe46wsa6+/3lqw+8E2B/tcq1bW7t/vd2XlUiXPccAAH9pS5KPG1EWkaHXqwKRJbhZ8w9AULmvh73+H3/wGvvwyKh8zYsQIMjIyWLNmTZlet2zZMu66666o1FCpPv0UTj/d3W0NmAXUSkqi38KFkJDgb20S8xTqInJol10Gq1bBBRdE9n38sVvl7MknKzyJrkWLFkyePLlU15Hnt3jxYi6++OIKfXalshaeesp1r3/1VWR3585MnjePxDJcGidSHIW6iJTsiCPcZW7/+Edktbk9e+CWW+Cii2Dr1gq9fb9+/ejatWvF66yqdu50q/UNHRpZird2bZgyhSvffpuzu3f3tz4JDIW6iJROtWpw++3wwQdw8smR/QsWuPb8+f7VVpUtWwannOJudRvWrh2sWAHXXusmJ4pEiUJdRMrm5JPdHcTyj2dv3w4XXww33QShWelxLzfXrdJ39tmRpV/BfVt//30oYklakYpSqItI2dWs6e7vvXix65oPmzDBjbWvWOFfbVXBli1w4YVulb7cXLevQQM36XDcuEpZV1/ik0JdRMrvvPNg5Uo3Xhy2erWbHf+3v0UCLZ68/rrrXl+yJLKvc2f45BM36VDEQwp1EamYRo1g1ix49ll3GRzA/v1uVbpzzoFvv/W3vsqSkwN/+AN07w7btrl9xrhv6xkZcNRRvpYn8UGhLiIVZ4y7Ccknn7hv6WFvveW+tb7wgm+lVYpvvoGzznLX8Icdfrgbnhg5EhK1eKdUDoW6iERPaiq8/TY88EBkIZVdu+Cqq9yjFHdkizmzZrnZ7f/7X2TfhRe6RWbOPde3siQ+KdRFJLoSE+H++124519Q5YUX3Lf2t97yr7Zoys6GwYPh8svh55/dvsREeOwxd5mfx2vSixRFoS4i3jjzTNcdf/31kX3ffeduCNO9Ozz9dGTsOdZ8/jl06OB+D2HHHAPvvgvDh7tr+kV8oL95IuKdunXdBLpZsyK3bbXWzRAfPBiaNYOuXWHsWNi40d9aS8NaeOYZOOMMF+xhl1/uls7t0MG/2kRQqItIZejTx60f37Nnwf15ea47/vbb3ezwjh3h0Ucjt3ytSnbtgn79YNAg+PVXt69WLRfyM2ZAvXr+1ieCQl1EKkvz5vDaa+4St3/8w80WP3iJ1OXL4e674fjjoW1bePBB98NABW8aUyG//OKuOT/1VJg5M7K/TRu3ZO6NN2qpV6kyFOoiUrmOOsp9M3/rLdi82a1Cd8EFhS/7WrXKzaJv2xbS0lzYL19eqoBftw5uvhlSUuCcc7qSkuLa69YV8wJr3fj+m2/CxIlwxx1u3L9lS3ft/XnnucvWwgYNcrW0aVPOkyDiDc8unjTGPAD8Dtge2vVHa+2CIo7rDvwTSAAmWWv/5lVNIlLFHHaYG1sfPBh27HA3hZkzx425h+9mBvD1165b/tFH4cgj3cpsvXu7ldoOugf5woWutz8nxz3AsHu3uzX81KmW2U9soUfTFe72p19+Gfl1586S601JcW/Ut29UT4NItHi9IsIYa+3jxT1pjEkAngTOBzYBHxhj5llrv/C4LhGpaho2hGuucY+sLHdZ2Msvuy77rKzIcRs3uol1Y8e6y8YuvRR69YJu3Vi3sTp9+hR9TxkX8oY+A1NYye2ksr50dSUkuOvvO3Z0wwHHHBOV366IF/xe5qgDsNZaux7AGDMDuARQqIvEszp13Izyyy93921fvNgF/Ny5Bb9Rb9vmLit7+mmoX59RyVPIye4BVC/2rXNIYgx3Mo5bCz6RnOzunNa6dcFfjzsOqhf/fiJVidehfosx5hrgQ2CYtfbg/q3mQP7rWDYBHT2uSURiSc2a7rauF1/svm6/+abron/lFdi6NXJcZibPZ6aTc4hAB8ihOs8lXMe4QV8UDO8WLTThTWKesRWYVWqM+S9weBFP3Qu8D/wIWOBhoJm19oaDXt8XuNBae2OoPQDoYK096EdoMMYMAgYBNGnS5LSZ+WehStRlZWVRJ3xzDvGMznMF5OZS74svaPzWWzR5+21qbt1KNXKxpZj/a4zljTferIQi44P+HnuvW7duK6y1p5d0XIVCvbSMMS2B+dbakw7afybwgLX2wlD7HgBr7SOHer+0tDS7evVqj6oVgIyMDNLT0/0uI/B0nqPEWvj4Y1K6nMzuX5NKPDwlxV12LtGhv8feM8aUKtQ9u6TNGNMsX/My4LMiDvsAON4Yc4wxpjrQD5jnVU0iElDGQPv2XH1dEkklZHpSEgwYUDlliVQ2L69T/7sxZpUxZiXQDbgTwBhzhDFmAYC1dj9wC/A68CUw01r7eXFvKCJyKMOGUapQv/POyqlHpLJ5NlHOWlvkz8LW2h+AnvnaC4BC16+LiJRVairMnn3wdepOUpJ7zJ7tjhMJIq0oJyKB0qMHrFzpFn1LSXGT4lJSXHvlSve8SFAp1EUkcFJTYdw4NxnujTfeZNcu19Y3dAk6hbqIiEhAKNRFREQCQqEuIiISEAp1ERGRgFCoi4iIBIRCXUREJCAU6iIiIgGhUBcREQkIhbqIiEhAKNRFREQCQqEuIiISEAp1ERGRgFCoi4iIBIRCXUREJCAU6iIiIgGhUBcREQkIhbqIiEhAKNRFREQCQqEuIiISEAp1ERGRgFCoi4iIBIRCXUREJCAU6iIiIgGhUBcREQkIhbqIiEhAKNRFREQCQqEuIiISEAp1ERGRgFCoi4iIBIRCXUREJCAU6iIiIgGhUBcREQkIhbqIiEhAKNRFREQCQqEuIiISEAp1ERGRgFCoi4iIBIRCXUREJCAU6iIiIgGhUBcREQkIhbqIiEhAKNRFREQCQqEuIiISEAp1ERGRgEj06o2NMS8BaaFmfSDTWntKEcdtAHYDucB+a+3pXtUkIiISZJ6FurX2ivC2MWYUsOsQh3ez1v7oVS0iIiLxwLNQDzPGGOBy4ByvP0tERCSeGWuttx9gzNnA6OK61Y0x3wA7AQtMtNY+Xcxxg4BBAE2aNDlt5syZHlUsAFlZWdSpU8fvMgJP59l7Osfe0zn2Xrdu3VaUZni6QqFujPkvcHgRT91rrZ0bOmY8sNZaO6qY9zjCWvuDMaYpsBi41Vr71qE+Ny0tza5evbrcdUvJMjIySE9P97uMwNN59p7Osfd0jr1njClVqFeo+91ae14JRSQCvYDTDvEeP4R+3WaMeQXoABwy1EVERKQwry9pOw/4ylq7qagnjTHJxpi64W3gAuAzj2sSEREJJK9DvR/wYv4dxpgjjDELQs3DgHeMMZ8Cy4HXrLWLPK5JREQkkDyd/W6tva6IfT8APUPb64F2XtYgIiISL7SinIiISEAo1EVERAJCoS4iIhIQCnUREZGAUKiLiIgEhEJdREQkIBTqIiIiAaFQFxERCQiFuoiISEAo1EVERAJCoS4iIhIQCnUREZGAUKiLiIgEhEJdREQkIBTqIiIiAaFQFxERCQiFuoiISEAo1EVERAJCoS4iIhIQCnUREZGAUKiLiIgEhEJdREQkIBTqIiIiAaFQFxERCQiFuoiISEAo1EVERAJCoS4iIhIQCnUREZGAUKiLiIgEhEJdREQkIBTqIiIiAaFQFxERCQiFuoiISEAo1EVERAJCoS4iIhIQCnUREZGAUKiLiIgEhEJdREQkIBTqIiIiAaFQFxERCQiFuoiISEAo1EVERAJCoS4iIhIQCnUREZGAUKiLiIgERIVC3RjT1xjzuTEmzxhz+kHP3WOMWWuMWW2MubCY1zc0xiw2xnwd+rVBReoRERGJZxX9pv4Z0At4K/9OY8yJQD+gDdAdeMoYk1DE6+8GllhrjweWhNoiIiJSDhUKdWvtl9ba1UU8dQkww1q711r7DbAW6FDMcVND21OBSytSj4iISDzzaky9ObAxX3tTaN/BDrPWbgYI/drUo3pEREQCL7GkA4wx/wUOL+Kpe621c4t7WRH7bFkKK6KOQcCgUHOvMeaziryflKgx8KPfRcQBnWfv6Rx7T+fYe2mlOajEULfWnleOD98EHJmv3QL4oYjjthpjmllrNxtjmgHbDlHH08DTAMaYD621pxd3rFScznHl0Hn2ns6x93SOvWeM+bA0x3nV/T4P6GeMqWGMOQY4HlhezHHXhravBYr75i8iIiIlqOglbZcZYzYBZwKvGWNeB7DWfg7MBL4AFgFDrbW5oddMynf529+A840xXwPnh9oiIiJSDiV2vx+KtfYV4JVinvsL8Jci9t+Yb/sn4NxyfPTT5XiNlI3OceXQefaezrH3dI69V6pzbKyt0Pw1ERERqSK0TKyIiEhAxFSoG2O6h5adXWuM0epzHjDGPGuM2aZLBr1jjDnSGLPUGPNlaJnl2/2uKWiMMTWNMcuNMZ+GzvGDftcUVMaYBGPMx8aY+X7XElTGmA3GmFXGmE9KmgUfM93voWVm1+Am1G0CPgCutNZ+4WthAWOMORvIAqZZa0/yu54gCl2+2cxa+5Expi6wArhUf5ejxxhjgGRrbZYxJgl4B7jdWvu+z6UFjjHmLuB0IMVae5Hf9QSRMWYDcLq1tsS1AGLpm3oHYK21dr21dh8wA7fMrESRtfYtYIffdQSZtXaztfaj0PZu4EuKXnFRysk6WaFmUugRG99gYogxpgXwf8Akv2sRJ5ZCvbRLz4rEDGNMS+BU4H8+lxI4oW7hT3CLWi221uocR98/gN8DeT7XEXQW+I8xZkVoddVixVKoR33pWRE/GWPqAHOAO6y1P/tdT9BYa3OttafgVrTsYIzRcFIUGWMuArZZa1f4XUsc6GytbQ/0AIaGhkmLFEuhXtqlZ0WqvNA47xxgurX2Zb/rCTJrbSaQgbsNtERPZ+C3ofHeGcA5xpjn/S0pmKy1P4R+3YZbG6aou54CsRXqHwDHG2OOMcZUx92vfZ7PNYmUWWgS12TgS2vtaL/rCSJjTBNjTP3Qdi3gPOArX4sKGGvtPdbaFtbalrj/j9+w1l7tc1mBY4xJDk2oxRiTDFwAFHt1UsyEurV2P3AL8DpuYtHM0HK0EkXGmBeB94A0Y8wmY8xAv2sKoM7AANw3m09Cj55+FxUwzYClxpiVuC8Ei621uuRKYtFhwDvGmE9x91B5zVq7qLiDY+aSNhERETm0mPmmLiIiIoemUBcREQkIhbqIiEhAKNRFREQCQqEuIiISEAp1ERGRgFCoi4iIBIRCXUREJCD+HybkWTH7fzhLAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize = (8, 8))\n",
    "ax.plot(x, y, lw = 3, color ='red')\n",
    "ax.scatter([1, 2, 3, 4], [3, -2, -5, 0], s = 100, color = 'blue', zorder = 3)\n",
    "ax.grid()\n",
    "ax.set_xlim([0, 5])\n",
    "ax.set_ylim([-10, 10])\n",
    "\n",
    "ax.text(1, 3.5, '$(1, 3)$', fontsize = 15)\n",
    "ax.text(1.5, -2.5, '$(2, -2)$', fontsize = 15)\n",
    "ax.text(2.7, -4, '$(3, -5.5)$', fontsize = 15)\n",
    "ax.text(4.1, 0, '$(4, .5)$', fontsize = 15)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now you know the trick, try another 5 points: $(1,2)$, $(2,5)$, $(3,8)$, $(4,6)$, $(5, 9)$. And polynomial form is \n",
    "\\begin{align}\n",
    "y=a_0+a_1x+a_2x^2+a_3x^3+a_4x^4\n",
    "\\end{align}\n",
    "\n",
    "The augmented matrix is"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}1 & 1 & 1 & 1 & 1 & 2\\\\1 & 2 & 4 & 8 & 16 & 5\\\\1 & 3 & 9 & 27 & 81 & 8\\\\1 & 4 & 16 & 64 & 256 & 6\\\\1 & 5 & 25 & 125 & 625 & 9\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "⎡1  1  1    1    1   2⎤\n",
       "⎢                     ⎥\n",
       "⎢1  2  4    8   16   5⎥\n",
       "⎢                     ⎥\n",
       "⎢1  3  9   27   81   8⎥\n",
       "⎢                     ⎥\n",
       "⎢1  4  16  64   256  6⎥\n",
       "⎢                     ⎥\n",
       "⎣1  5  25  125  625  9⎦"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = sy.Matrix([[1, 1, 1, 1, 1, 2],\n",
    "               [1, 2, 4, 8, 16, 5], \n",
    "               [1, 3, 9, 27, 81, 8], \n",
    "               [1, 4, 16, 64, 256, 6], \n",
    "               [1, 5, 25,125, 625, 9]]); A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 19.   , -37.417,  26.875,  -7.083,   0.625])"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A_rref = A.rref()\n",
    "A_rref = np.array(A_rref[0])\n",
    "coef = A_rref.astype(float)[:,-1];coef"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.linspace(0, 6, 100)\n",
    "y = coef[0] + coef[1]*x + coef[2]*x**2 + coef[3]*x**3 + coef[4]*x**4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeAAAAHSCAYAAAA0ZhgzAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABCB0lEQVR4nO3deXzV1Z3/8dfJvm8QQkJYww4CQhSXqgFX1FZr7WZ17DKlU1vHOk47ttOZbr922k6t06nTmVp1StUWHa3VulMVgYrIIvueAAGSQMi+kPWe3x/35i4SSAg393uX9/Px4OH93ntz74dTypvvWY21FhEREQmtOKcLEBERiUUKYBEREQcogEVERBygABYREXGAAlhERMQBCmAREREHJITyy0aOHGknTJgQtM9ra2sjPT09aJ8X6dQegdQePmqLQGoPH7VFoGC3x8aNG09Ya/P7ey2kATxhwgQ2bNgQtM9buXIlZWVlQfu8SKf2CKT28FFbBFJ7+KgtAgW7PYwxh073mrqgRUREHKAAFhERcYACWERExAEKYBEREQcogEVERBygABYREXGAAlhERMQBCmAREREHKIBFREQcoAAWERFxgAJYRETEAQpgERERByiARUREHKAAFhERcYACWERExAEKYBEREQckOF2AiIiI01bsPMa3/7SNzLhudlHOl8tKhv07FcAiIhLzDte3c6y5k2PA4Yb2kHynuqBFRCTmVTWe9D4ek5Maku9UAIuISMyralIAi4iIhNzRxg7v4yIFsIiISGgcbfC7A85VAIuIiAy7ju5eTrR2AmCAgszkkHyvAlhERGJaTZOv+zk3xZAQH5poVACLiEhMO+o3A3pkqgnZ9yqARUQkpvkHcF6KAlhERCQk/NcAj0gJXSwqgEVEJKb5z4AeoS5oERGR0PDfhEMBLCIiEiJVfptwqAtaREQkBKy1AZOwdAcsIiISAidau+jqcQGQlZJAaoICWEREZNj5z4AO1R7QfRTAIiISs/wDuDhEe0D3GXQAG2PijTHvG2Ne9FznGWNWGGP2ef6bO3xlioiIBN/RCLkDvgfY5Xd9P/CGtXYK8IbnWkREJGKEfQAbY4qBG4BH/J6+CVjmebwMuDmolYmIiAwz/y7oMSEOYGOtHfhNxjwD/BuQCfyjtfZGY0yjtTbH7z0N1tpTuqGNMUuBpQAFBQULli9fHqzaaW1tJSMjI2ifF+nUHoHUHj5qi0BqD59Yb4vvvHOSQ83uWdDfXpjC6MSTQW2PRYsWbbTWlvb3WsJAP2yMuRE4bq3daIwpO9svt9Y+DDwMUFpaasvKzvojTmvlypUE8/MindojkNrDR20RSO3hE+tt0bJ6BdAFwI2LL2X3+++GrD0GDGDgUuAjxpjrgRQgyxjzBHDMGFNora02xhQCx4ezUBERkWA62dVLfZs7fBPjDaMyk9kdwu8fcAzYWvtNa22xtXYC8CngTWvt7cALwJ2et90JPD9sVYqIiASZ/x7Qo7NTiIsL3SYccG7rgH8MXG2M2Qdc7bkWERGJCP6nIBVlh3YCFgyuC9rLWrsSWOl5XAdcGfySREREhl/ADOgQb8IB2glLRERilJNLkEABLCIiMeqIg5twgAJYRERilO6ARUREHFDV2OF9rDtgERGREHC5LNVN/l3QKSGvQQEsIiIxp7a1k+5e91bMeelJpCWd1aKgoFAAi4hIzAk8BSn0d7+gABYRkRjk9CYcoAAWEZEYVOXwEiRQAIuISAw66vASJFAAi4hIDDpc3+59PDZPASwiIhISR/zGgItz0xypQQEsIiIxxVobEMBjFcAiIiLDr66ti5PdvQBkpiSQnZboSB0KYBERiSkB478O3f2CAlhERGJM4PivMxOwQAEsIiIx5nCD/wxo3QGLiIiExOF63QGLiIiE3JEGjQGLiIiEXMAYsEObcIACWEREYojLZQMOYnBqEw5QAIuISAw53tJJV68LgNy0RDKSQ38OcB8FsIiIxIwjYTIDGhTAIiISQw6HyQQsUACLiEgMORImS5BAASwiIjHE/w64WF3QIiIioREum3CAAlhERGLIkUaNAYuIiIRUT6+LqsYO77XugEVEREKgprmDXpcFID8zmZTEeEfrUQCLiEhM8B//Hevw3S8ogEVEJEb4b8Lh5BaUfRTAIiISEw777QE91sFDGPoogEVEJCYcqdcdsIiISMj5H0Po9BIkUACLiEiMCNgFS5OwREREhl9Xj4uaZvcaYGOgKEcBLCIiMuyqGk9i3UuAGZ2VQlKC8/HnfAUiIiLDLNzGf0EBLCIiMSDwFCTnu59hEAFsjEkxxrxnjNlijNlhjPme5/nvGmOOGmM2e35dP/zlioiInL1w24QDIGEQ7+kEFltrW40xicAaY8wrntcetNb+bPjKExEROXeVYXQMYZ8BA9haa4FWz2Wi55cdzqJERESCqbKuzft4fF543AEPagzYGBNvjNkMHAdWWGvXeV76qjFmqzHmMWNM7nAVKSIici4O+e2CNX5EuoOV+BhrB38za4zJAZ4D7gZqgRO474Z/ABRaaz/fz88sBZYCFBQULFi+fPm5V+3R2tpKRkZG0D4v0qk9Aqk9fNQWgdQePrHQFm3dlq+84Q7gxDj49dVpxBnT73uD3R6LFi3aaK0t7e+1wYwBe1lrG40xK4Hr/Md+jTG/AV48zc88DDwMUFpaasvKys7mK89o5cqVBPPzIp3aI5Daw0dtEUjt4RMLbbHtSBO8sQaACSMzWLzoitO+N5TtMZhZ0PmeO1+MManAVcBuY0yh39s+CmwflgpFRETOwaF6v/HfEeEx/guDuwMuBJYZY+JxB/bT1toXjTGPG2Pm4e6CPgh8adiqFBERGaJDdb7x33F54TH+C4ObBb0VOL+f5+8YlopERESCqLLOfwJW+NwBaycsERGJav5d0OMUwCIiIqERcAccJmuAQQEsIiJRrLOnl2rPMYRxJny2oQQFsIiIRLHD9b5jCAuzU8PiGMI+4VOJiIhIkFWG6RIkUACLiEgUOxSmM6BBASwiIlEsXNcAgwJYRESiWGW97oBFRERC7lCdxoBFRERCyuWyHG446b0Ol2MI+yiARUQkKtU0d9DV4wJgRHoSGclndQDgsFMAi4hIVAqYgBVm3c+gABYRkSgVsAY4jLag7KMAFhGRqBR4Bxxe47+gABYRkSh1qD48D2HoowAWEZGoFM5LkEABLCIiUchaq0lYIiIiodbY3k1LRw8AaUnx5GckO1zRqRTAIiISdfzHf8flpWGMcbCa/imARUQk6viP/44LwwlYoAAWEZEoVBnGxxD2UQCLiEjUCeiCDsM1wKAAFhGRKBRwB6wuaBERkdCoOOEbA544UnfAIiIiw66lo5sTrZ0AJMXHUZST6nBF/VMAi4hIVDl4InACVnxc+C1BAgWwiIhEmYoTrd7H4dr9DApgERGJMgf8x3/zFcAiIiIh4R/Ak3QHLCIiEhr+ATwhTNcAgwJYRESiiLVWXdAiIiKhVtfW5T0FKSM5ISxPQeqjABYRkahx4AMbcITjKUh9FMAiIhI1DtSG/w5YfRTAIiISNfy3oJygABYREQmNgxGyBAkUwCIiEkU+OAYczhTAIiISFVwuy4E6dUGLiIiEVFXTSbp6XACMzEgiOzXR4YrOTAEsIiJRIVJ2wOozYAAbY1KMMe8ZY7YYY3YYY77neT7PGLPCGLPP89/c4S9XRESkfwcjaPwXBncH3AksttbOBeYB1xljLgLuB96w1k4B3vBci4iIOKIiQrag7DNgAFu3vsMVEz2/LHATsMzz/DLg5uEoUEREZDAi5RSkPoMaAzbGxBtjNgPHgRXW2nVAgbW2GsDz31HDVqWIiMgADkTQJhwAxlo7+DcbkwM8B9wNrLHW5vi91mCtPWUc2BizFFgKUFBQsGD58uXnWLJPa2srGRkZQfu8SKf2CKT28FFbBFJ7+ERLW/S4LF98vZ2+RHv46jSS4s9+H+hgt8eiRYs2WmtL+3st4Ww+yFrbaIxZCVwHHDPGFFprq40xhbjvjvv7mYeBhwFKS0ttWVnZ2XzlGa1cuZJgfl6kU3sEUnv4qC0CqT18oqUtymtbsa+/DcCYnFSuuXLRkD4nlO0xmFnQ+Z47X4wxqcBVwG7gBeBOz9vuBJ4fphpFRETOKJIOYegzmDvgQmCZMSYed2A/ba190RizFnjaGPMFoBL4+DDWKSIiclqRtAVlnwED2Fq7FTi/n+frgCuHoygREZGzEUmnIPXRTlgiIhLxDpxo9T6OhCVIoAAWEZEoUBGBY8AKYBERiWjNHd0cb+kEICk+jrF5aQ5XNDgKYBERiWj7j/u6nyeOTCc+7uzX/zpBASwiIhHNP4Anj4qcTUUUwCIiEtHK/QK4RAEsIiISGroDFhERcUB5rV8A5yuARUREhl1Hdy+V9e0AGAOTIuAc4D4KYBERiVgH69pweY5AGpubRkpivLMFnQUFsIiIRCz/8d+SCLr7BQWwiIhEsEidgAUKYBERiWAKYBEREQcogEVEREKs12UDjiGcnJ/pYDVnTwEsIiIR6WjDSbp6XACMzEgmOy3R4YrOjgJYREQi0v7aFu/jSJsBDQpgERGJUJE8/gsKYBERiVAKYBEREQcogEVERELMWqsAFhERCbXa1k6aO3oAyEhOYHRWisMVnT0FsIiIRJzy4771vyX56RhjHKxmaBTAIiIScfbX+h/CEHndz6AAFhGRCFTufwpSBI7/ggJYREQiUKRPwAIFsIiIRCAFsIiISIg1neymprkDgKT4OMblpTlc0dAogEVEJKLsPebbA3pSfjqJ8ZEZZZFZtYiIxKw9Nb4Anj46so4g9KcAFhGRiOJ/BzxVASwiIhIa/nfA0woUwCIiIsPOWssevzvgaboDFhERGX61LZ00tncDkJ4Uz5icVIcrGjoFsIiIRIw9Hxj/jcQ9oPsogEVEJGJEy/gvKIBFRCSCBARwBI//ggJYREQiiP8SJN0Bi4iIhIDLZdl7zLcHdCSvAQYFsIiIRIjDDe2c7O4FYGRGEiMzkh2u6NwMGMDGmLHGmLeMMbuMMTuMMfd4nv+uMeaoMWaz59f1w1+uiIjEKv/x36kR3v0MkDCI9/QA91lrNxljMoGNxpgVntcetNb+bPjKExERcQvYgjIWAthaWw1Uex63GGN2AWOGuzARERF/e/zGfyN9BjSc5RiwMWYCcD6wzvPUV40xW40xjxljcoNdnIiISJ89Nc3ex9EQwMZaO7g3GpMBvA380Fr7R2NMAXACsMAPgEJr7ef7+bmlwFKAgoKCBcuXLw9W7bS2tpKRkRG0z4t0ao9Aag8ftUUgtYdPpLRFj8vypRXt9Hoi67+vSiM1Ifi7YAW7PRYtWrTRWlva32uDCmBjTCLwIvCatfbn/bw+AXjRWjv7TJ9TWlpqN2zYMKiiB2PlypWUlZUF7fMindojkNrDR20RSO3hEyltsaemhWv/YxUAY3JS+ev9i4fle4LdHsaY0wbwYGZBG+BRYJd/+BpjCv3e9lFg+7kWKiIi0p9oOQHJ32BmQV8K3AFsM8Zs9jz3LeDTxph5uLugDwJfGob6REREom78FwY3C3oN0F9H+8vBL0dERORUe2r8ZkBHwRIk0E5YIiISAaJtDTAogEVEJMy1dfZQWd8OQHycoWRUusMVBYcCWEREwtpuv/Hfkvx0khPiHawmeBTAIiIS1nZW+QJ4ZmGWg5UElwJYRETC2s5qXwDPKsp2sJLgUgCLiEhYC7gDLtIdsIiIyLDr6XWx2+8YwhnqghYRERl+FSfa6OxxAVCYnUJeepLDFQWPAlhERMJWtE7AAgWwiIiEMf8JWNE0/gsKYBERCWP+d8CzFMAiIiLDz1obeAdcGD1LkEABLCIiYepYcyf1bV0AZCYnUJyb6nBFwaUAFhGRsLSjqsn7eEZhFnFx/R3MF7kUwCIiEpaidQOOPgpgEREJS4HjvwpgERGRkIjmJUigABYRkTDU0tHNoTr3GcAJcYYpBRkOVxR8CmAREQk7/vs/Tx6VETVnAPtTAIuISNiJ9glYoAAWEZEwFM17QPdRAIuISNjZUe1bAzyrKLp2wOqjABYRkbDS3etib02r91p3wCIiIiGw71grXb3uM4DH5KSSnZbocEXDQwEsIiJhZdvRRu/j88ZEZ/czKIBFRCTMbDniG/+dM1YBLCIiEhLb/AN4TI5zhQwzBbCIiISNzp5edtf4liCpC1pERCQEdle30N1rAZgwIi1qJ2CBAlhERMLI1iON3sdzinMcqyMUFMAiIhI2tvqP/xZHb/czKIBFRCSMBAZwjnOFhIACWEREwkJ7Vw/7jrtPQYozMCtKD2HoowAWEZGwsKOqGZd7/hWTR2WQnpzgbEHDTAEsIiJhwb/7+bwoXv/bRwEsIiJhwX8G9Nwo3gGrjwJYRETCwraAO2AFsIiIyLBrOtlNxYk2ABLiDDOi9AhCfwpgERFx3I6jvrvfaaMzSUmMd7Ca0FAAi4iI47YejZ31v30UwCIi4rjALSijf/wXBhHAxpixxpi3jDG7jDE7jDH3eJ7PM8asMMbs8/w3d/jLFRGRaBRLW1D2GcwdcA9wn7V2BnAR8BVjzEzgfuANa+0U4A3PtYiIyFmpa+3kSMNJAJIT4phakOlwRaExYABba6uttZs8j1uAXcAY4CZgmedty4Cbh6lGERGJYv53vzOLskiMj43RUWOtHfybjZkArAJmA5XW2hy/1xqstad0QxtjlgJLAQoKChYsX778HEv2aW1tJSMjI2ifF+nUHoHUHj5qi0BqD59waItn93bx54puAK4Zn8BtM5IdqyXY7bFo0aKN1trS/l4b9EabxpgM4Fnga9baZmPMoH7OWvsw8DBAaWmpLSsrG+xXDmjlypUE8/MindojkNrDR20RSO3hEw5t8fC+d4E6AG66dA5lcwodqyWU7TGo+3xjTCLu8H3SWvtHz9PHjDGFntcLgePDU6KIiESrnl4Xmw83eq/nj89xrJZQG8wsaAM8Cuyy1v7c76UXgDs9j+8Eng9+eSIiEs32HGuhvasXgMLsFAqzUx2uKHQG0wV9KXAHsM0Ys9nz3LeAHwNPG2O+AFQCHx+WCkVEJGptOtTgfTx/fGytZh0wgK21a4DTDfheGdxyREQklmyqbPQ+nj8utgI4NuZ6i4hIWNrodwe8IMbugBXAIiLiiBOtnVTWtwPuDThmxsAJSP4iMoBrWzp5dXs1z+7r4pVt1U6XIyIiQ+A//junOJukhIiMpCEb9DrgcLJi5zG+9dw2AGx6NUvOc27NmIiIDM3GSr8JWDE2/gsRegc8q8jXTbGzqtnBSkREZKjeP9TofXy+AjgyTBudSXyce2J2xYk2Wjt7HK5IRETORlePiy1+RxDG0gYcfSIygFMS45mc79urc1e17oJFRCLJrupmOntcAIzNS2VUZorDFYVeRAYwBHZD7zjadIZ3iohIuAlYfhSD3c8QwQE80z+ANQ4sIhJRNlXG7g5YfSI2gGePyfY+VgCLiESW92N4B6w+ERvA/nfA+4630NnT62A1IiIyWDVNHRxtPAlAamI800dnOlyRMyI2gLNSEslPdc+E7u617DvW6nBFIiIyGP7jv3PHZpMQH7FRdE4i+nc9PstX/o4qTcQSEYkE6w7UeR9fMCHPwUqcFdEBPC4ggDUOLCISCdZV1HsfL5w4wsFKnBXRATxeASwiElHq27rYc6wFgIQ4E5MbcPSJ7ADO9JW/q7qZXpd1sBoRERnIe37dz3OKs0lLisgjCYIiogM4JyWO/MxkANq7ejlY1+ZwRSIicibv+nc/T4rd7meI8ACGwB2xtmtHLBGRsLbugP/4b+xOwIIoC2CdjCQiEr6a2rvZXeP+ezo+zlAawzOgISoCWDtiiYhEgvcO1mM9U3VmF2WRkRy7478QFQHsvyd0E9ZqIpaISDhaV+GbgBXr478QBQE8NjeNTM+/ohrau6lu6nC4IhER6Y/GfwNFfADHxRlm6GQkEZGw1tzR7d2x0BhifvwXoiCAAWYHjANrJrSISLjZeLCBvq0aZhZmkZ2a6GxBYSAqAjhwKZLugEVEws27fhtwxPL2k/6iI4DHaC2wiEg4C9j/eZK6nyFKAnhyfgapifEA1DR3cKxZE7FERMJFa2cP2/xuji7U+C8QJQGcEB/HeWN848BbDjc6V4yIiATYeKjBu1f/9NGZ5KYnOVxReIiKAAb3pt59th5RN7SISLh4p/yE97GWH/lETQDPHZvjfbzlSKNjdYiISKA1+3wB/KEp+Q5WEl6iJoDn+Qfw4UZcOppQRMRxJ1o7vfszxMcZLtIELK+oCeDi3FTyPOMKzR09OppQRCQM/HW/7+53/rgcMlO0/rdP1ASwMUbjwCIiYWa1f/fzZHU/+4uaAAaYW5zjfbxZM6FFRBxlrQ0Y/71s6kgHqwk/URXA8zQRS0QkbOw/3kqNZ1+GzJQE5vgtF5UoC2D/LugdVc109bgcrEZEJLb5dz9fWjKShPioipxzFlWtMSIjmeLcVAC6elzsPdbicEUiIrFrzX7/5Ufqfv6gqApgCFwPrHFgERFndPW4eLfCdwDD5Vr/e4qoC+B5fhOxtCWliIgzNlU20N7VC8C4vDTGjUhzuKLwE3UB7H8HrKVIIiLOWL2v1vv4MnU/92vAADbGPGaMOW6M2e733HeNMUeNMZs9v64f3jIHb/aYLOKM+/He4y20dvY4W5CISAwKWH6kAO7XYO6Afwtc18/zD1pr53l+vRzcsoYuLSmBqQWZAFir84FFREKtoa2LrZ6/e+MMXFyiAO7PgAFsrV0F1A/0vnAyV+PAIiKOeae8DuvZjn/u2ByyU7X9ZH+MtQMfWmCMmQC8aK2d7bn+LvBZoBnYANxnrW04zc8uBZYCFBQULFi+fHkw6gagtbWVjIyMU55/q7KbZTu7ALhgdDxfmZcStO8MZ6drj1il9vBRWwRSe/gMR1s8uq2T1Ufdw38fKUnklimRc/5vsNtj0aJFG621pf29NtQALgBOABb4AVBorf38QJ9TWlpqN2zYcBaln9nKlSspKys75fntR5u48ZdrABiTk8pf718ctO8MZ6drj1il9vBRWwRSe/gEuy1cLsuFP/oLJ1rdN0HP3XUJ54/LDdrnD7dgt4cx5rQBPKRZ0NbaY9baXmutC/gNcOG5FBhs00ZnkpoYD8DRxpPUNHU4XJGISGzYfKTRG74jM5IDhgQl0JAC2BhT6Hf5UWD76d7rhMT4OOaO9W1LueFQRA1hi4hErDd2HfM+Xjw9n7i+ZSlyisEsQ/oDsBaYZow5Yoz5AvBTY8w2Y8xWYBFw7zDXedZKx/sOfd5wsN/haRERCbI3dh33Pr5yRoGDlYS/hIHeYK39dD9PPzoMtQTVggm+MYeNhxTAIiLD7UhDO7tr3HvwJyXE8aHJWn50JlG3E1af+eNyMZ6ej53VzbRpQw4RkWHlf/d78aQRpCcPeI8X06I2gLNTE5nm2ZCj12V1MIOIyDD7i9/471UzRjlYSWSI2gAGWDDe1w2tcWARkeHT2tnDugrfhNfFGv8dUFQHcKnfOLBmQouIDJ/Ve2vp6nUBMKMwizE5qQ5XFP6iO4D9ZkK/X9lIr2vgTUdEROTs/cVv/Ffdz4MT1QFcnJvKqMxkwN09ssczO09ERIKn12V5a4+WH52tqA5gY0xAN/RGdUOLiATd5sMN1Le5d7/Kz0xmzpjsAX5CIMoDGD6wIYfWA4uIBN2Knb6738XTRmn3q0GK/gCeoJnQIiLDxVrLK9urvddXavx30KI+gGcUZgUczFDddNLhikREoseOqmYO1bUDkJGcwOVT8x2uKHJEfQAnxscxb2yO91p3wSIiwfPSNt/d71UzRpHiueGRgUV9AAMfmIilABYRCQZrLS/7BfANc4ocrCbyxEQAB+yIpZnQIiJB8cHu58um6PCFsxETATx/vN/BDFXNtHR0O1uQiEgU8O9+vnpmgbqfz1JMBHBWSiIzRmcB4LIaBxYROVfWWl7a6tf9fF6hg9VEppgIYIBLSkZ4H79TfsLBSkREIt/2o81U1ru7nzOTE7hsqrqfz1bMBPDFfgG8tqLOwUpERCLfB7ufkxPU/Xy2YiaAL5iYR9/mLDuqmmlq1ziwiMhQWGt5aVuV9/qGOep+HoqYCeCslETO8+xPai2sO6C7YBGRodh+tJnD9e5NjTKTE/iQZj8PScwEMMDFJb4/JO+UK4BFRIbiRb+736tnqft5qGIsgH3jwO9qHFhE5Kz1uiwvbPbrftbs5yGLqQAuHZ9LgmcgeHdNC3WtnQ5XJCISWd4pP0F1UwcAI9KTtPfzOYipAE5PTmCu377Q6w5oVywRkbPxzMYj3sc3zRtDYnxMxUhQxVzLXTxJ64FFRIaiuaOb13bUeK8/tmCMg9VEvpgLYP8NOdZqIpaIyKC9vLWajm4X4D7qdVZRtsMVRbaYC+D543NJ8nSZlNe2cby5w+GKREQig3/3860Lih2sJDrEXACnJMZz/rgc77V2xRIRGdiBE21s8BznmhBnuGmejh48VzEXwPCBbSnVDS0iMqBn/e5+y6aNYmRGsoPVRIeYDOBL/Dbk0B2wiMiZuVyWP25S93OwxWQAzx2bTUqi+7d+qK6dIw3tDlckIhK+1lbUUeVZ+5ublsji6aMcrig6xGQAJyfEc+FEXzf0qr1ajiQicjofXPublBCT0RF0MduKV/jt3rJyz3EHKxERCV91rZ0BRw9+bL66n4MlZgO4bJovgN8pr6Orx+VgNSIi4empDYe9fz/OLc7mvGKt/Q2WmA3gSSPTKc5NBaC1s4eNnun1IiLi1uuyPPlupff6by6e4FwxUShmA9gYE9AN/fbeWgerEREJP2/sOsbRRve5v3npSdwwRycfBVPMBjC417L10TiwiEig36095H38yQvGkpKoc3+DKaYD+OKSESTG+44nPKZtKUVEANh/vIU1+90rROIMfGbhOIcrij4xHcAZyQlcMCHPe/32HnVDi4gAPO5393vVjAKKc9McrCY6xXQAAxoHFhH5gNbOHp7ddNR7feclE5wrJorFfAD7jwOv3ldLT6+WI4lIbHtu0xFaO3sAKMlPDzjGVYJnwAA2xjxmjDlujNnu91yeMWaFMWaf57+5w1vm8JlakMHorBQAmjt62Hy40dmCREQc5HJZfvvOQe/131w8AWOMcwVFscHcAf8WuO4Dz90PvGGtnQK84bmOSMaYgE051A0tIrHs9Z01lNe2Ae55MrfMH+NwRdFrwAC21q4C6j/w9E3AMs/jZcDNwS0rtAK3pVQAi0hsstbyX2+Ve69vv2g8mSmJDlYU3YY6Blxgra0G8Pw3oo/GuHTKSOLj3F0s2442UdvS6XBFIiKht3rfCbYdbQIgOSGOL3xoosMVRTdjrR34TcZMAF601s72XDdaa3P8Xm+w1vY7DmyMWQosBSgoKFiwfPnyIJTt1traSkZGRlA+60frTrK3wT0B6/Ozk7i8OPL+1RfM9ogGag8ftUUgtYePf1v827qT7PH8PXjluATumJnsZGmOCPafjUWLFm201pb291rCED/zmDGm0FpbbYwpBE67jZS19mHgYYDS0lJbVlY2xK881cqVKwnW5+2Lq+CHL+8C4FBPDmVlFwTlc0MpmO0RDdQePmqLQGoPn762WH+wnj2vrgUgIc7w/dsuZ0xOqsPVhV4o/2wMtQv6BeBOz+M7geeDU45zrp012vt49f4T3in4IiKx4Fdv7fc+/uj5Y2IyfENtMMuQ/gCsBaYZY44YY74A/Bi42hizD7jacx3Rxo1IY/roTAC6elzaFUtEYsaOqibe8vydZwz8XVmJwxXFhgG7oK21nz7NS1cGuRbHXTNrNLtrWgD3VHyd/CEiUam8HB54AJ54gitaW+lITuMHM67gNxd8lPMun09JvsbHQyHmd8Lyd+2sAu/jN3cf9x5CLSISNV55BebMgUcegZYWjLWkdrTxyS2v8epjX+XrrgqnK4wZCmA/MwuzKM51j3u0dPSwtqLO4YpERIKovBxuvRXa26G7O+ClJFcvaT2dTPi7O93vk2GnAPZjjAmYjPXajhoHqxERCbIHHjgleE/R3Q0PPhiaemKcAvgD/AN4xc5juFwDr5MWEYkITzwxuAB+/PHQ1BPjFMAfsGB8LiPSkwCobenkfR3OICLRorU1uO+Tc6IA/oD4OMNVM3yTsV5XN7SIRIvB7vCkXcJCQgHcj2tn+wL4tR01DGa7ThGRsHf77ZA4wDa7iYlwxx2hqSfGKYD7cUnJSNKT4gE4WNfO3mPqjhGRKHDffdjBBPC994amnhinAO5HSmI8ZdN9Bzy9uLXKwWpERIKkpIRf3/0T2hOS6YqLD3wtMRHS0uCZZ6BEO2GFggL4ND48p8j7+PnNVeqGlshTXg533QVZWVyxeDFkZbmvtcYzZr29t5YfM5HrPv8Qf5h7HT0ZmVhj3H82li6FrVthyRKny4wZCuDTWDQ9n8wU906dlfXtbKpsdLYgkbPRz25HtLS4r+fMcb8uMeVkVy/f/tM2ACpzC9n0jR+Q0NLM22++CU1N8NBDuvMNMQXwaSQnxHP9bN9e0M9vPupgNSJn4Qy7HdHd7X7+1lt1Jxxj/uONvRyuPwlAdmoi/3LjTIcrEgXwGdx0vq8b+qWt1XT3am9oiQAPPIAdYLMF291NzwMPaGglRuysauaR1Qe819+6fjojM5IdrEhgEKchxbKFE0cwOiuFmuYO6tq6WLP/BIumjRr4B0WGkbWW2tZOdle3cKi+nSMN7RxtOMmRhpOcaO3k1UeXkTFAAJvubk4+uowFuR9mZEYS+ZnJ5GcmMyorhXF5aUwcmc7EkemMy0sjJTH+jJ8l4a3XZfnmc9vo9ezqt3BiHp8oHetwVQIK4DOKjzN8eG4hv/H8y/H5948qgCXk6tu6WH+wno2HGthZ1cyu6mbq2rpO+/60rpOD+tz0rpN09bqoauqgqqmj3/cYAxNHpjOrKJvZRVnMKspmzthsslIGWMoiYeM/39jHFs+OfknxcfzolvMwxjhblAAK4AHdNG+MN4Bf33mM9q4e0pLUbDJ8Wjt7WLPvBKv21fLegXr2Hz+7dehtSSlkDiKE25JSB3yPtVBR20ZFbRt/3uJejhdnYProLC6YkMsFE/O4aNIIdWeGqbXldfzyzX3e67+/crLO+g0jSpIBzCrKYvKoDPYfb6W9q5cVO49x07wxTpclUaa66SSvbq/hzd3HWVdRT9cA8w3SkuKZNjqTkvwMxuamUZybSnFuKgVZKSTW/w32t49hztQNnZhI5hc/x87vX8uJli5qWzuobemkqrGDg3VtHDjRxsG6No42nOSD55G4LOysbmZndTPL1h4C3P8/uWxKPpdPHUnp+DySEjS9xGl1rZ187an3vf/7XTQpjy+XTXa2KAmgAB6AMYab5xXxs9f3Au41wQpgCYaGti5e3l7N85urWH+wntPNh0qMN5w3JpsLJuZx/tgcZhRmMTY3jbi403Qj/tPX4cnHz3zqjWe3o7SkBMaNSGDciLR+33ayq5ddNc3sqGpmZ1UTW480sau6+ZRQ3lHlfs//vF1OZnICZdNHcc3MAsqm5ZOp7uqQs9byj/+3hWPNnQDkpSfxi0+dT/zp/syIIxTAg/CRuWO8Abxqby31bV3keU5MEjkbLpdl9f4T/H7dId7YdZye0xx3OaMwi8XT8/nQ5Hzmjc0hNeksJkKVlLh3M7r1VncI+wdxYqL71yB3O0pNimf+uFzmj8v1PtfS0c37lY2sP1jPuop6NlU2BPw+Wjp7+POWKv68pYqk+DgunTyCD88t4uqZBQrjEHl0zQHe2lPrvX7g43MpyEpxsCLpjwJ4EMaNSGP+uBw2VTbS47K8sPkon710otNlSQQ50drJ0xsO84f3Kr1rMf3FGbh08kiumz2aRdNGUZQz8PjsGS1Z4t7V6MEH4fHHsS0tmMxM9yb79957ThsuZKYkcvnUfC6fmg+4A3lteR2r9tWyck8tRxp8v7+uXhdv7anlrT21JCXEsXjaKD4yr4jF00dpdvUwWVdRx09e3e29/uJlE1k0XZNHw5ECeJA+Or/YuxvWH947zJ2XTNBMQhlQeW0rj6w+wLObjtDVc+q47ryxOdw8r4jr5xQyKjPIdyglJe7djR56iLdXrqSsrCy4n++RmZLINbNGc82s0Vhr2V3Twus7jvH6zhp2VDV739fV4+LVHTW8uqOGrJQEPjy3iI8tKOb8sTn6/1KQ7D/eytLHN9Ld6+6RmFOczdevne5wVXI6CuBBumleET96aRcnu3vZc6yFTZWNLBifO/APSkzaeKiB/3m7nL/sOnbK2G52aiIfm1/MbQvHMnlUpjMFDhNjDDMKs5hRmMU9V03hcH07L26t5s9bqthZ7Qvj5o4enlxXyZPrKpk0Mp1PXDCWj80vJj9Ts6mHqralk8/+73s0nXQPOYzMSOa/bpuvCXFhTAE8SFkpiXxkbhFPbTgMwO/XVSqA5RQbDzXwH3/Zy+p9J055bU5xNp+9ZALXn1cYM92vY/PS+HJZCV8uK2H/8VZe2FLFc+8fCeiGrzjRxo9f2c3PXtvDVTMK+OSFY7l8Sr4mDJ2F9q4e/nbZem/3f2piPI99tpSxef1PrpPwoAA+C59eOM4bwC9ureJfb5xJdpomlQhsPtzIz1fsZdXe2lNeu2rGKL542SQunJgX012tk0dl8A9XT+VrV05h/cF6nt10hJe31dDa2QNAj8t6u6jH5KRy28JxfLy0OPhd81Gm12W5Z/lmthxpAtzzCR667XzmFOc4W5gMSAF8FuYWZzOzMIud1c109rh47v0jmowV4w6eaOOnr+3m5W01Ac/HGbj5/DHcVVYSdd3M5youzrBw0ggWThrBdz8yi5e31fDU+krWH2zwvudo40n+/bU9PLhiL9fMKuD2heO5uGRETP8Dpj9dPS7ufWozK3Ye8z73vY/M4soZBQ5WJYOlAD4LxhhuWziOb/9pO6DJWLGsvq2L/3xjH0+uO+Sd8ALu4L1p3hjuXjyZSdpxaEBpSQncuqCYWxcUs/94K0+tr+SZjUdoaHePY/a4LC9vq+HlbTWU5KfzmYXj+diCYrJT1fPU0d3Ll5/YGLDcaOnlk7jj4gnOFSVnRQF8lm6aV8SPXt5Fe1ffZKwGFozPc7osCZGeXhe/W3uIB/+yl5aOnoDXbjivkHuvnsrkUQreoZg8KoN/vmEm910zjdd21PDku5W8d7De+3p5bRvff3EnP31tNzfPG8PtF41n9phsByt2TmtnD19ctoG1FXXe5z536QS+uUQzniOJAvgsZXomYy1f7x4LfnJdpQI4Rqwtr+O7L+xgz7GWgOcvmJDLt66fwfnjNCkvGFIS47lp3hhumjeGvcdaePLdQzy76ah3rLij28Xy9YdZvv4w88bmcMdF47lhTuxMbKtv6+ILy9bzvmdZJMDdiyfzD1dPVW9chFEAD8GnLxznDeCXtlbznRtnaTJWFDve3MH3X9zJi1urA56fODKd+5dM55qZBfqLb5hMLcjkezfN5hvXTef5zVX8bu1Bdtf4/gG0+XAjmw838oOXdnLr/GJuWzguqrv+N1U28JUnN1Htd3rVP103nS+XDX1jFXGOFogNwZzibGYVZQHQ2ePimU1HHK5IhoPLZXn83UNc+cDbAeGblhTPP103nVe/dhnXzhqt8A2B9OQEbls4jlfuuYxnv3wxN88rIine99dXY3s3j6w5wOIH3ua237zLn7dU0dnTG/gh5eVw112QlcUVixdDVpb7urw8xL+bs2etZdk7B/nkr9cGhO/3PjJL4RvBdAc8BMYYPrNwPN96bhsAj605wJ0XjychXv+eiRZ7alr41nPb2HioIeD5m+YV8c0lMxidraUxTjDGsGB8HgvG5/EvN3by1IbD/H5dZcD2l++U1/FOeR25aYncMr+YT14wlqmb1gTsjW0AWlrgkUdg2TL33thLljj2+zqT5o5uvv3cdl7wHAcJ7s1cHvzkXBZP12znSKYAHqKPnj+Gn72+h/q2Lo42nuSlbdU6JSkKdPW4+NXK/Tz05v6AAwYmjkznhx+dzSUlIx2sTvyNyEjmrrLJ/N3lJazaV8uT6yp5Y9cx70lNDe3dPLrmACv+/A6v//ZuUro6Tv2QvsMqbr3VvXf2OeyRHWzWWv60+Sg/fGk3J1o7vc+fNyabX31mvjbZiAIK4CFKTYrnjovG84s33IddP7yqgo/MLVJ3ZATbfrSJrz+zlV1+WyYmxhu+fEUJdy2aHDOTfCJNXJyhbNooyqaNoqrxJE9vOMz/bTjC0Ub3XfEX1z9HXM8ZjmYEdwg/+KB77+wwsPdYC//yp+2sO1Af8PxtC8fxrzfO1J/FKKEAPgd/c/F4/uftcjp7XOyoauad8jounaw7pEjT1ePioTf38auV5QF3vfPH5fCTj81hSoE20ogURTmpfO2qqdy9eAp/3X+Cp9Yf5uYH3yLJ1XvmH+zuhscfdzyAd1U385vVFbywuSrgz+LorBT+9cMzuf68Qgerk2BTAJ+DERnJfLy0mCferQTg16sqFMARZk9NC/c+tTngoICUxDi+fu10PnvJBO1HHKHi44z3yER7Rz9dz/1wtbTw5NqDLJo+iuLc0HXvulyWd8rreHh1xSlbmSbEGT7/oYn8/ZVTyEjWX9fRRv+LnqO//dAknlxXibWwam8tu6qbmVGY5XRZMgCXy/LomgP8+2t76Or1HRN44YQ8fnLrHCaOTHewOgkmk5HhnnA1gLbEVP7l+R3w/A4m5aezcGIeF0xw/yrOTQ3q8FJ3r4t1FfW8tqOG13fWcKy585T3XFLi3qpzqnpgopYC+BxNGJnOdbNG88p2917Av1lVwc8/Oc/ZouSMjjS0c9/TWwLG15IT4vjGddP53CUTiNNdb3S5/Xb3bOfu048Dd8XF88dZi7zXFbVtVNS28Yf33Ov9R2UmM70wi+mjM5lakMmUURmMzk5hRHrSgKsfOnt6qWnqYEdVM9uPNrGjqpn3Kxto/sBOauDeyvS62aP54mWTtLFLDFAAB8HSyyd5A/iFLVX847XTKMpJdbgq6c8LW6r45+e2BWwjed6YbB785FwdmhCt7rvPvdToDAGckJxM5v1fp6wtlbXldXT2uAJeP97SyfGW2lO6iOOMeyjKHcSG+Lg44g24LDS0d1Hf2kVL56lB+0G5ae4d9j7/oYmMH6Hel1ihAA6C88flcuGEPN47WE+Py/Kb1RV858OznC5L/LR29vCd53fwrN+mKXEGvrpoMndfOYVEreGOXiUl7nW+fuuAvRITITGRuGee4ZYlV3AL7kMONh9uZP2Bet47WM+mQw20dfU/ictlobalk9qWU7uQB1KYncK1s0Zz7azRXDAhV/sIxCAFcJB86YpJ3o3jn1xXyRcvm6S7YCeUl8MDD8ATT3BFaytkZHDi5o9zd9Ei1uLr0hubl8p/fPJ8FoxXN19MWLLEvc73wQfh8cexLS2YzEy44w64996A9b8pifFcNGkEF00aAbjP2z1Y18aemhZ217Swp6aZyvqTHG/uoK6ta8Cvjo8zjEhPYmpBJrPGZDGryL2T3qSR6Vq2GOPOKYCNMQeBFqAX6LHWlgajqEi0ePoo5o3NYfPhRrp6XPziL/v4ya1znC4rtrzySr+7HWU9uYxH457krpu/ycqSUm45fwzfu2kWmSnavzumlJS4lxk99BBvr1xJWVnZoH4sPs5Qkp9BSX7GKcuAunpcnGjtpKG9C5cLelwuej3Lh3LTkxiRnkRWSqLmFUi/gnEHvMhaeyIInxPRjDF849pp3PbIOgD+b+Nhll4xiZIo3hg+rJSXu8O3vf2Ul5JcvSS5evnv5/+NNc++wdUfnhf6+iQqJSXEUZSTqt4uGRINOgTRJZNH8iHPOmCXhZ+/vtfhimLIAw+ccZINQAourn7tDyEqSETkzM41gC3wujFmozFmaTAKinRfv3aa9/FL26rZfrTJwWpiyBNPDBjApm+3IxGRMGCstQO/63Q/bEyRtbbKGDMKWAHcba1d9YH3LAWWAhQUFCxYvnz5udQboLW1lYyM8Ovi/eX7HWw85p41ed7IeO4rDc3JOeHaHqFwxeLFmEH8WbbG8Pabb4agovASy382+qP28FFbBAp2eyxatGjj6eZHnVMAB3yQMd8FWq21Pzvde0pLS+2GDRuC8n0AK89iIkUo7T/ewjUPrvKeyvLU0otY6JlROZzCtT2G28o9xymdO5GMzlPHf0+RlQVNsdcrEat/Nk5H7eGjtggU7PYwxpw2gIfcBW2MSTfGZPY9Bq4Btg/186LJ5FGZ3DK/2Hv9g5d2emdGSvD09Lr4yau7+ez/rue5mWV0xQ1wQkxionvZiYhIGDiXMeACYI0xZgvwHvCStfbV4JQV+b521RSSE9zNu/1oM0+8e8jhiqJLVeNJPvXwu/z3ynIAfnPBR+mNH2BSf2Kie82niEgYGHIAW2srrLVzPb9mWWt/GMzCIl1xbhpfXTTZe/2z1/ZwvHlwp7LImb21+zg3/OdqNhxq8D43YeEcupc/DWlp7qD1l5jofv6ZZ8LqwHURiW1ahjSMll4xiUn57n1dWzp7+MFLuxyuKLJ197r4t5d38bnfrqeh3T3jOT7O8I3rpvHbz15A1i0fce92tHQpZGVhjXGP+S5d6n5+yRKHfwciIj4K4GGUnBDP/7tptvf6z1uqWL2v9gw/IadzpKGdT/x6Lb9eVeF9bnRWCsuXXsRdZZN9Ow317XbU1OSe7dzU5L7Wna+IhBkF8DC7ZPJIbp5X5L3+lz9tp6O7/43dpX+v76jh+l+s5v3KRu9zZdPyefmey7hgQp5zhYmInAMFcAj88w0zyUxxTxA6WNfOr97a73BFkaGju5fvvrCDpY9v9J6dGh9n+OaS6Tx25wXkpSc5XKGIyNApgEMgPzOZb1w33Xv9XyvL2XCw/gw/IeW1rdzyq3f47TsHvc8VZafw9Jcu5ktXlGhzexGJeArgELntwnFc6Oku7XVZ7lm+maaTZ946MRZZa3l24xE+/Ms17Kxu9j5/9cwCXr7nMh0fKCJRQwEcIvFxhgc/NY/sVPcSmaONJ/nmH7cSrJ3IokHTyW7uWb6Z+/5vC+2eA9CT4uP4/k2zePiOBeSkqctZRKKHAjiExuSk8pOP+c4IfnlbDcvXH3awovCx/mA91/9iNS9sqfI+Nyk/nee+cgl/c/EEHVwuIlFHARxi180eze0XjfNef+/PO9h3rMXBipzV3evi56/v4ZO/XsvRxpPe5z9RWsyfv/ohZhVlO1idiMjwUQA74Ns3zGRqgfu0jY5uF3/3xEYa27scrir09h9v4ZZfvcN/vrnfe3BFdmoiv/rMfH5661zSkwfYWlJEJIIpgB2QkhjPQ7fN9+4VXV7bxtLfbYyZ9cEul+WR1RVc/59r2OZ3XvJFk/J45Z7LuP68QgerExEJDQWwQ6YWZPKzj8/1Xr93sJ77/m8Lrig/Namyrp1P/+Zd/t9Lu+jqcQHuiVbfXDKdJ//2IopyUh2uUEQkNNTH56APzy2iuukkP3p5NwAvba1mTE4q37p+hsOVBV+vy/K/fz3Az17fQ0e3y/v8rKIsfv6JeUwbnelgdSIioacAdtgXL5vE0YaTLFvrPq7w4VUVjM5K4fMfmuhwZcGzp6aFbzy7lS2HG73PxccZ7ior4e7FU0hKUEeMiMQeBbDDjDH864dnUd3Uwes7jwHw/Rd3crK7l7vKSiJ6+U17Vw+/fHM/j6yuoLvX17U+fXQmP711DnOKc5wrTkTEYQrgMBAfZ/jFp87nM4+8yybPgQP//toe6lq7+PYNMyJu20VrLa9sr+EHL+6kusl3BnJSfBx/f+VkvnRFCYnxuusVkdimAA4TqUnx/O4LC1n6uw28U14HwGN/PUBDexc/vXVOxATW3mMt/ODFnazedyLg+dLxufz4Y+cxeZTGekVEQAEcVjKSE/jfz13AvU9t5uVtNQA89/5RTrR28vNPzCM/M9nhCk+vuukkD67YyzMbj+A/kXtEehLfvH4Gt5w/JuLu5EVEhpMCOMwkJ8Tzy0/PJydtO79fVwnA6n0nWPKLVfz7rXNZNH2UwxUGamzv4terKnhszQE6e3yzm+MM3HHReP7hmmne/a9FRMRHARyG4uMMP7x5NvkZyfzijX0AnGjt4nO/Xc+dF4/nm9fPICUx3v3m8nJ44AF44gmuaG2FjAy4/Xa47z4oKRm2Go81d/DI6gp+v66Stq7ADUQun5rP/ddNZ2ZR1rB9v4hIpFMAhyljDPdePZXSCbnc9/QWjrd0ArBs7SFW7TvBP14zjSWH3yfuEx+H7m7o7sYAtLTAI4/AsmXwzDOwZElQ69p/vIVH1xzg2Y1H6ep1Bbx23phs7l8ynUsnjwzqd4qIRCMFcJi7bEo+r37tcv7p2a2s8CxTOnCijZ/810tc+du7SenqOPWHPIHMrbfC1q3nfCfc3tXDS1urWb7+MBsPNZzy+tSCDL66eAo3nleocV4RkUFSAEeAvPQkHr5jAb9/r5Ifv7ybls4evrj+OeJ6us/8g93d8OCD8NBDZ/2drZ09rNpby192HWPFjmO0dPac8p7543K4q2wyi6ePUvCKiJwlBXCEMMbwmYXjueG8Qv7n7Qo++uBbJLkGOLyhuxsef3xQAVzf1sX2o01sO9rEuxV1vFtRF7B5Rp+EOMPVMwu485IJLJyYF9EbhYiIOEkBHGFy0pK4f8l0bHc/Xc/9cLW0cPfvN5GZnEC65xdAfVsnDW3d1LV1UlnXTlXTmT9v0sh0PnnBWD62oJiRGeG7HEpEJFIogCOUychwT7gaQFtiKi9trR7Sd8wszOKqGaO4ckYBc4qzdbcrIhJECuBIdfvt7tnO3acfB+6Ki+ePsxYN6uOSEuKYUZjF7KIszhuTzWVT8xmjowFFRIaNAjhS3Xefe6nRGQI4MSWZ2T/7Dr/IK6K1s4e2zh5aO3vBWvLSk8hNT2JEejKjspKZODI9Yra7FBGJBgrgSFVS4l7ne+utvmVHfRITITER88wzLFh8AQucq1JERE5DtzyRbMkS9zrfpUshKwtrDGRlua+3bg36JhwiIhI8CuBIV1LiXmbU1MTbb74JTU3u62HchlJERM6dAlhERMQBCmAREREHKIBFREQcoAAWERFxgAJYRETEAQpgERERByiARUREHKAAFhERcYACWERExAHnFMDGmOuMMXuMMfuNMfcHqygREZFoN+QANsbEA/8FLAFmAp82xswMVmEiIiLR7FzugC8E9ltrK6y1XcBy4KbglCUiIhLdziWAxwCH/a6PeJ4TERGRARhr7dB+0JiPA9daa//Wc30HcKG19u4PvG8psBSgoKBgwfLly8+tYj+tra1kZGQE7fMindojkNrDR20RSO3ho7YIFOz2WLRo0UZrbWl/ryWcw+ceAcb6XRcDVR98k7X2YeBhgNLSUltWVnYOXxlo5cqVBPPzIp3aI5Daw0dtEUjt4aO2CBTK9jiXO+AEYC9wJXAUWA/cZq3dcYafqQUODekL+zcSOBHEz4t0ao9Aag8ftUUgtYeP2iJQsNtjvLU2v78XhnwHbK3tMcZ8FXgNiAceO1P4en6m3yKGyhiz4XS39rFI7RFI7eGjtgik9vBRWwQKZXucSxc01tqXgZeDVIuIiEjM0E5YIiIiDoj0AH7Y6QLCjNojkNrDR20RSO3ho7YIFLL2GPIkLBERERm6SL8DFhERiUgRG8A6CMLHGPOYMea4MWa707U4zRgz1hjzljFmlzFmhzHmHqdrcpIxJsUY854xZounPb7ndE1OM8bEG2PeN8a86HQtTjPGHDTGbDPGbDbGbHC6HqcZY3KMMc8YY3Z7/g65eFi/LxK7oD0HQewFrsa9Ich64NPW2p2OFuYQY8zlQCvwO2vtbKfrcZIxphAotNZuMsZkAhuBm2P4z4YB0q21rcaYRGANcI+19l2HS3OMMeYfgFIgy1p7o9P1OMkYcxAotdZqHTBgjFkGrLbWPmKMSQLSrLWNw/V9kXoHrIMg/FhrVwH1TtcRDqy11dbaTZ7HLcAuYniPcuvW6rlM9PyKvH91B4kxphi4AXjE6VokvBhjsoDLgUcBrLVdwxm+ELkBrIMgZEDGmAnA+cA6h0txlKfLdTNwHFhhrY3l9vgP4BuAy+E6woUFXjfGbPTs2x/LJgG1wP96higeMcakD+cXRmoAm36ei9l/1cupjDEZwLPA16y1zU7X4yRrba+1dh7u/dovNMbE5DCFMeZG4Li1dqPTtYSRS62183Gf6/4Vz3BWrEoA5gP/ba09H2gDhnV+UaQG8KAOgpDY5BnrfBZ40lr7R6frCRee7rSVwHXOVuKYS4GPeMY9lwOLjTFPOFuSs6y1VZ7/Hgeewz28F6uOAEf8eoiewR3IwyZSA3g9MMUYM9EzUP4p4AWHa5Iw4Jl09Ciwy1r7c6frcZoxJt8Yk+N5nApcBex2tCiHWGu/aa0tttZOwP13xpvW2tsdLssxxph0z0RFPF2t1wAxu5LCWlsDHDbGTPM8dSUwrJM3z2kvaKcM5SCIaGaM+QNQBow0xhwBvmOtfdTZqhxzKXAHsM0z7gnwLc++5bGoEFjmWTkQBzxtrY355TcCQAHwnPvfrCQAv7fWvupsSY67G3jSc2NXAXxuOL8sIpchiYiIRLpI7YIWERGJaApgERERByiARUREHKAAFhERcYACWERExAEKYBEREQcogEVERBygABYREXHA/wfg0OzEKQAVagAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(figsize= (8, 8))\n",
    "ax.plot(x, y, lw =3)\n",
    "ax.scatter([1, 2, 3, 4, 5], [2, 5, 8, 6, 9], s= 100, color = 'red', zorder = 3)\n",
    "ax.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <font face=\"gotham\" color=\"purple\"> Solving The System of Linear Equations By NumPy </font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Set up the system $A  x =  b$, generate a random $A$ and $b$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "A = np.round(10 * np.random.rand(5, 5))\n",
    "b = np.round(10 * np.random.rand(5,))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-2.283,  1.431,  0.677, -0.718,  1.908])"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = np.linalg.solve(A, b);x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's verify if $ Ax = b$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-0.,  0.,  0.,  0., -0.])"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A@x - b"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "They are technically zeros, due to some round-off errors omitted, that's why there is $-$ in front $0$."
   ]
  }
 ],
 "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.8.8"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}