{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "db5748ac",
   "metadata": {},
   "source": [
    "# Linear Algebra I\n",
    "\n",
    "Linear algebra is a core topic in modern applied mathematics. Essentially every important method in statistics, data science, and machine learning is built on linear algebra. Indeed, deep neural networks, which we will discuss shortly, are built on a foundation of matrix multiplication coupled with simple nonlinear functions. \n",
    "\n",
    "<figure class=\"image\" style=\"width:30%\">\n",
    "  <img src=\"https://imgs.xkcd.com/comics/machine_learning_2x.png\" alt=\"\">\n",
    "  <figcaption><i></i></figcaption>\n",
    "</figure>\n",
    "\n",
    "In this lecture, we'll see how to perform several important operations in linear algebra using our good friend, Numpy. These operations include: \n",
    "\n",
    "- Matrix multiplication. \n",
    "- Exact and approximate solutions to linear systems. \n",
    "- Singular value and eigenvalue-eigenvector decompositions. \n",
    "\n",
    "Along the way, we'll show several examples from statistics and applied mathematics, including simulation of partial differential equations; least-squares regression; and image compression. \n",
    "\n",
    "This is probably the lecture in which things will get \"the most mathematical.\" Comfort with concepts from Math 33A or equivalent will be helpful. If you're not familiar with these concepts, that's ok -- feel free to ask questions. We'll all get through this just fine. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "051120b9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# no fancy packages this time! just our good friends numpy and matplotlib\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "np.random.seed(1234)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4b5b3ff1",
   "metadata": {},
   "source": [
    "## Basic Matrix Operations\n",
    "\n",
    "A *matrix* is a two-dimensional array of numbers. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "0155704f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[4, 7, 6, 5, 9],\n",
       "       [2, 8, 7, 9, 1],\n",
       "       [6, 1, 7, 3, 1],\n",
       "       [6, 3, 7, 4, 8],\n",
       "       [1, 1, 4, 3, 4]])"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# random matrix data to play with\n",
    "n = 5\n",
    "A = np.random.randint(1, 10, size = (n, n))\n",
    "A"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cedfe61e",
   "metadata": {},
   "source": [
    "Matrices admit several standard operations, including: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "8fc6a898",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 8, 14, 12, 10, 18],\n",
       "       [ 4, 16, 14, 18,  2],\n",
       "       [12,  2, 14,  6,  2],\n",
       "       [12,  6, 14,  8, 16],\n",
       "       [ 2,  2,  8,  6,  8]])"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# scalar multiplication\n",
    "2*A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "cc59f510",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[4, 2, 6, 6, 1],\n",
       "       [7, 8, 1, 3, 1],\n",
       "       [6, 7, 7, 7, 4],\n",
       "       [5, 9, 3, 4, 3],\n",
       "       [9, 1, 1, 8, 4]])"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# transposition\n",
    "A.T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "74293632",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[4. , 4.5, 6. , 5.5, 5. ],\n",
       "       [4.5, 8. , 4. , 6. , 1. ],\n",
       "       [6. , 4. , 7. , 5. , 2.5],\n",
       "       [5.5, 6. , 5. , 4. , 5.5],\n",
       "       [5. , 1. , 2.5, 5.5, 4. ]])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# application of transposition: a \"symmetrized\" version of A\n",
    "# symmetric matrices satisfy A = A.T\n",
    "(A + A.T)/2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0ebb0d40",
   "metadata": {},
   "source": [
    "Inversion is an especially important matrix operation. The inverse $\\mathbf{A}^{-1}$ of a square matrix $\\mathbf{A}$ satisfies $\\mathbf{A}\\mathbf{A}^{-1} = \\mathbf{I}$, where $\\mathbf{I}$ is the identity matrix. We'll see how to multiply matrices and check this in a sec. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "80fa14b2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-0.30292479,  0.12743733, -0.1810585 ,  0.58704735, -0.47910864],\n",
       "       [ 0.38857939, -0.07381616,  0.19777159, -0.42200557, -0.06128134],\n",
       "       [ 0.48746518, -0.23955432,  0.5097493 , -0.84122563,  0.51810585],\n",
       "       [-0.65529248,  0.33774373, -0.51810585,  0.88370474, -0.24791086],\n",
       "       [-0.01740947, -0.02715877, -0.12534819,  0.13718663,  0.05292479]])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# inverse of A\n",
    "A_inv = np.linalg.inv(A) \n",
    "A_inv"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ab81e0ea",
   "metadata": {},
   "source": [
    "## Matrix multiplication"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "8be3aba9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# random vector\n",
    "b = np.random.randint(1, 5, size = n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "82692921",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([98, 86, 52, 88, 47])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# matrix-vector product\n",
    "A.dot(b)\n",
    "\n",
    "# modern, convenient syntax -- same effect\n",
    "A@b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "e160aba8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ True,  True,  True,  True,  True])"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A.dot(b) == A@b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "ec830806",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[172,  91, 202, 149, 203],\n",
       "       [166,  89, 169, 108, 148],\n",
       "       [ 85,  56, 135,  90, 120],\n",
       "       [146,  86, 193, 145, 191],\n",
       "       [ 78,  42,  86,  74,  83]])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# random matrix\n",
    "B = np.random.randint(1, 10, size = (n, n))\n",
    "\n",
    "# matrix-matrix product (same as A.dot(B))\n",
    "A@B"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "06a17aa4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.00000000e+00,  2.77555756e-17,  8.32667268e-17,\n",
       "        -1.38777878e-16, -8.32667268e-17],\n",
       "       [-1.49186219e-16,  1.00000000e+00, -2.49800181e-16,\n",
       "         8.32667268e-17,  5.55111512e-17],\n",
       "       [ 7.28583860e-17, -8.32667268e-17,  1.00000000e+00,\n",
       "         8.32667268e-17,  0.00000000e+00],\n",
       "       [ 1.38777878e-16, -2.22044605e-16, -2.22044605e-16,\n",
       "         1.00000000e+00,  2.22044605e-16],\n",
       "       [-4.16333634e-17, -1.11022302e-16,  0.00000000e+00,\n",
       "         7.77156117e-16,  1.00000000e+00]])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# checking our inverse from earlier\n",
    "# observe--finite precision arithmetic! \n",
    "A @ A_inv"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "dbc4f3f5",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x7fab3a49cc40>"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPUAAAD4CAYAAAA0L6C7AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAJEUlEQVR4nO3dz2ucBR7H8c9ns7G1uCC79mCbsvUgskXYCKEIvRWh9Qd6VdCTkMsKFQTRo/+AePFSVFxQFEEPIi6hrIoIbjVqLHajUsTFUqHuivhj2dbqZw+ZQ9dNOs9MnmeezJf3CwJJJsx8CHnnmZmEZ5xEAOr4Vd8DALSLqIFiiBoohqiBYogaKObXXVzpVb+dyd49s11cdes+PbGj7wnAyP6jH3Q+57zeZZ1EvXfPrN5Z2tPFVbfu0K75vicAIzuev254GXe/gWKIGiiGqIFiiBoohqiBYogaKIaogWKIGiiGqIFiiBoohqiBYogaKIaogWKIGiiGqIFiiBoohqiBYhpFbfuw7U9sn7L9UNejAIxvaNS2ZyQ9LulmSfsk3WV7X9fDAIynyZF6v6RTST5Lcl7S85Lu6HYWgHE1iXq3pC8u+vj04HP/w/ai7WXby1/966e29gEYUZOo1zsN6f+9ql6So0kWkizs/N3M5pcBGEuTqE9Luvh8v3OSznQzB8BmNYn6XUnX2r7G9mWS7pT0crezAIxr6Mn8k1ywfZ+kJUkzkp5KcrLzZQDG0ugVOpK8KunVjrcAaAH/UQYUQ9RAMUQNFEPUQDFEDRRD1EAxRA0UQ9RAMUQNFEPUQDFEDRRD1EAxRA0UQ9RAMUQNFEPUQDGNTpIwqk9P7NChXfNdXHXrls6s9D1hJNPyfUV/OFIDxRA1UAxRA8UQNVAMUQPFEDVQDFEDxRA1UAxRA8UQNVAMUQPFEDVQDFEDxRA1UAxRA8UQNVAMUQPFDI3a9lO2z9r+aBKDAGxOkyP105IOd7wDQEuGRp3kTUlfT2ALgBbwmBooprWzidpelLQoSdu1o62rBTCi1o7USY4mWUiyMKttbV0tgBFx9xsopsmftJ6T9Lak62yftn1v97MAjGvoY+okd01iCIB2cPcbKIaogWKIGiiGqIFiiBoohqiBYogaKIaogWKIGiiGqIFiiBoohqiBYogaKIaogWKIGiiGqIFiWjvx4LQ6tGu+7wkjWTqz0veEkUzb97cCjtRAMUQNFEPUQDFEDRRD1EAxRA0UQ9RAMUQNFEPUQDFEDRRD1EAxRA0UQ9RAMUQNFEPUQDFEDRRD1EAxRA0UMzRq23tsv2571fZJ20cmMQzAeJqco+yCpAeSvG/7N5Les30syd873gZgDEOP1Em+TPL+4P3vJK1K2t31MADjGelsorb3SrpB0vF1LluUtChJ27WjjW0AxtD4iTLbV0h6UdL9Sb795eVJjiZZSLIwq21tbgQwgkZR257VWtDPJnmp20kANqPJs9+W9KSk1SSPdj8JwGY0OVIfkHSPpIO2VwZvt3S8C8CYhj5RluQtSZ7AFgAt4D/KgGKIGiiGqIFiiBoohqiBYogaKIaogWKIGiiGqIFiiBoohqiBYogaKIaogWKIGiiGqIFiiBooZqSziaJ/h3bN9z1hJEtnVvqe0Ni0fW83wpEaKIaogWKIGiiGqIFiiBoohqiBYogaKIaogWKIGiiGqIFiiBoohqiBYogaKIaogWKIGiiGqIFiiBooZmjUtrfbfsf2h7ZP2n5kEsMAjKfJ6YzOSTqY5Hvbs5Lesv2XJH/reBuAMQyNOkkkfT/4cHbwli5HARhfo8fUtmdsr0g6K+lYkuOdrgIwtkZRJ/kpybykOUn7bV//y6+xvWh72fbyjzrX8kwATY307HeSbyS9IenwOpcdTbKQZGFW29pZB2BkTZ793mn7ysH7l0u6SdLHHe8CMKYmz35fLenPtme09kvghSSvdDsLwLiaPPt9QtINE9gCoAX8RxlQDFEDxRA1UAxRA8UQNVAMUQPFEDVQDFEDxRA1UAxRA8UQNVAMUQPFEDVQDFEDxRA1UAxRA8U0OfMJMLZDu+b7ntDY0pmVvic0tv/Qvze8jCM1UAxRA8UQNVAMUQPFEDVQDFEDxRA1UAxRA8UQNVAMUQPFEDVQDFEDxRA1UAxRA8UQNVAMUQPFEDVQDFEDxTSO2vaM7Q9sv9LlIACbM8qR+oik1a6GAGhHo6htz0m6VdIT3c4BsFlNj9SPSXpQ0s8bfYHtRdvLtpd/1Lk2tgEYw9Cobd8m6WyS9y71dUmOJllIsjCrba0NBDCaJkfqA5Jut/25pOclHbT9TKerAIxtaNRJHk4yl2SvpDslvZbk7s6XARgLf6cGihnpZXeSvCHpjU6WAGgFR2qgGKIGiiFqoBiiBoohaqAYogaKIWqgGKIGiiFqoBiiBoohaqAYogaKIWqgGKIGiiFqoBiiBopxkvav1P5K0j9avtqrJP2z5evs0jTtnaat0nTt7Wrr75PsXO+CTqLugu3lJAt972hqmvZO01Zpuvb2sZW730AxRA0UM01RH+17wIimae80bZWma+/Et07NY2oAzUzTkRpAA0QNFDMVUds+bPsT26dsP9T3nkux/ZTts7Y/6nvLMLb32H7d9qrtk7aP9L1pI7a3237H9oeDrY/0vakJ2zO2P7D9yqRuc8tHbXtG0uOSbpa0T9Jdtvf1u+qSnpZ0uO8RDV2Q9ECSP0i6UdKftvD39pykg0n+KGle0mHbN/Y7qZEjklYneYNbPmpJ+yWdSvJZkvNae+XNO3retKEkb0r6uu8dTST5Msn7g/e/09oP3+5+V60va74ffDg7eNvSz/LanpN0q6QnJnm70xD1bklfXPTxaW3RH7xpZnuvpBskHe95yoYGd2VXJJ2VdCzJlt068JikByX9PMkbnYaovc7ntvRv6Glj+wpJL0q6P8m3fe/ZSJKfksxLmpO03/b1PU/akO3bJJ1N8t6kb3saoj4tac9FH89JOtPTlnJsz2ot6GeTvNT3niaSfKO1V1/dys9dHJB0u+3PtfaQ8aDtZyZxw9MQ9buSrrV9je3LtPbC9y/3vKkE25b0pKTVJI/2vedSbO+0feXg/csl3STp415HXUKSh5PMJdmrtZ/Z15LcPYnb3vJRJ7kg6T5JS1p7IueFJCf7XbUx289JelvSdbZP2763702XcEDSPVo7iqwM3m7pe9QGrpb0uu0TWvtFfyzJxP5MNE34N1GgmC1/pAYwGqIGiiFqoBiiBoohaqAYogaKIWqgmP8CUqLe49gT7FkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.imshow(A@A_inv) # looks like the identity matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "7093bd65",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ True,  True,  True,  True,  True],\n",
       "       [ True,  True,  True,  True,  True],\n",
       "       [ True,  True,  True,  True,  True],\n",
       "       [ True,  True,  True,  True,  True],\n",
       "       [ True,  True,  True,  True,  True]])"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# identity matrix\n",
    "I_n = np.eye(n)\n",
    "\n",
    "# check the result to within machine precision, entrywise\n",
    "np.isclose(I_n, A@A_inv)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "8afaf3b4",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# aggregates the above\n",
    "np.allclose(I_n, A@A_inv)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7b4d46fe",
   "metadata": {},
   "source": [
    "## Application: Simulating Heat Diffusion\n",
    "\n",
    "Matrix multiplication provides a simple way to simulate 1-dimensional partial differential equations in discrete time. For example, the 1-d heat equation reads\n",
    "\n",
    "$$\\frac{\\partial f(x, t)}{\\partial t} = \\frac{\\partial^2 f}{\\partial x^2 }\\;.$$\n",
    "\n",
    "In a discrete approximation, we can write this as \n",
    "\n",
    "$$f(x_i, t + 1) - f(x_i, t) \\approx \\epsilon\\left[f(x_{i+1}, t) - 2f(x_i, t) + f(x_{i-1}, t)\\right]\\;,$$\n",
    "\n",
    "where $\\epsilon$ is a small positive number that controls the timescale of the approximation. \n",
    "\n",
    "We can write the righthand side of this equation as a matrix-vector product:\n",
    "\n",
    "- Let $\\mathbf{v}(t)$ be the values of $f(\\mathbf{x}, t)$ -- that is, $v_i = f(x_i)$. \n",
    "- Let $\\mathbf{A}$ be the *transition operator*: \n",
    "\n",
    "$$\n",
    "\\mathbf{A} = \\left[\\begin{matrix}\n",
    "    -2 &  1 & 0 & \\cdots& 0& 0 & 0\\\\\n",
    "     1 & -2 & 1 & \\cdots & 0& 0 &  0\\\\\n",
    "     0 & 1 & -2 & \\cdots & 0& 0 &  0\\\\\n",
    "     \\vdots & \\vdots &\\vdots & \\ddots &  \\vdots & \\vdots & \\vdots \\\\ \n",
    "     0  & 0  & 0  & \\cdots   & -2     & 1 & 0\\\\\n",
    "     0  & 0  & 0  & \\cdots   & 1     & -2 & 1\\\\\n",
    "     0 & 0  & 0  & \\cdots & 0     & 1     & -2 \\\\\n",
    "\\end{matrix}\\right]\n",
    "$$\n",
    "\n",
    "This transition operator has the property that $[\\mathbf{A}\\mathbf{v}]_i$ is equal to the righthand side of the discrete approximation, for $i = 2,\\ldots,n-1$. That is, we can write \n",
    "\n",
    "$$\n",
    "\\mathbf{v}(t+1) = \\mathbf{v}(t) + \\epsilon \\mathbf{A}\\mathbf{v}(t) = (\\mathbf{I} + \\epsilon\\mathbf{A})\\mathbf{v}(t)\n",
    "$$\n",
    "\n",
    "Note that there are issues at the boundary (i.e. where $i = 1$ or $i = n$), and further modeling decisions must be made in order to handle these. In the transition operator above, we are effectively allowing heat to escape at the boundaries. \n",
    "\n",
    "To simulate heat diffusion in Python, we can just build this transition operator as a matrix (`numpy` array) and then iterate this update. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "be489814",
   "metadata": {},
   "outputs": [],
   "source": [
    "# size of simulation\n",
    "n = 201\n",
    "\n",
    "# Construct A using the handy np.diag() function\n",
    "A = -2*np.eye(n)\n",
    "A += np.diag(np.ones(n-1), 1)\n",
    "A += np.diag(np.ones(n-1), -1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "f2a563e7",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7fab3a6afc10>]"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAATcUlEQVR4nO3dbYycV3mH8ev2mqSiAQJ4iVzbwU4xtP7AS7IEJBpKRQE7bXEpVZuACk1BVqSkAqFKSYVKkegXikAVImC51ApUlNCKUFxkGtqqDR9Q2mzSvDkhwQRIjNNkQxBUBBqcufthzniPZ3e9s5t52bNcP2nlmWfO7t46M/nn7HnOc57ITCRJ7dsw6QIkScNhoEvSOmGgS9I6YaBL0jphoEvSOrFxUr9406ZNuX379kn9eklq0i233PJoZk4v9trEAn379u3Mzs5O6tdLUpMi4jtLveaUiyStEwa6JK0TBrokrRMGuiStEwa6JK0TywZ6RByMiEci4q4lXo+I+GhEHI2IOyLi/OGXKUlaziAj9GuB3ad5fQ+ws3ztAz7x1MuSJK3UsoGemV8FHjtNk73Ap7PrJuDsiNg8rAKlccpMrr/1GD9+4slJlyKt2DDm0LcAD1bPj5VjC0TEvoiYjYjZubm5IfxqabgeeOxx3vP3t/Ov9zw86VKkFRtGoMcixxa9a0ZmHsjMmcycmZ5e9MpVaaJ++mQHgBOdzoQrkVZuGIF+DNhWPd8KHB/Cz5XGrlOGIua5WjSMQD8EvK2sdnkl8IPMfGgIP1cau065JWPHWzOqQctuzhURnwVeA2yKiGPAnwNPA8jM/cBh4GLgKPA4cNmoipVGrTcyN8/VomUDPTMvXeb1BK4YWkXSBDlCV8u8UlSq9HK8Y56rQQa6VHGErpYZ6FKlF+RpoKtBBrpU6TjlooYZ6FIlnXJRwwx0qeIIXS0z0KWKc+hqmYEuVVzlopYZ6FLFdehqmYEuVRyhq2UGulTpjczNc7XIQJcqJ0fozrmoQQa6VOmtbjHO1SIDXar0ts91Dl0tMtClyvxJ0QkXIq2CgS5V5k+Kmuhqj4EuncJli2qXgS5V3MtFLTPQpYoXFqllBrpU8cIitcxAlyrphUVqmIEuVVy2qJYZ6FLFC4vUMgNdqniDC7XMQJcq7oeulhnoUsVli2qZgS5VvLBILTPQpYpz6GqZgS5V0ikXNcxAlypOuahlBrpU8aSoWjZQoEfE7oi4NyKORsTVi7z+rIj4p4i4PSKORMRlwy9VGj33clHLlg30iJgCrgH2ALuASyNiV1+zK4C7M/MlwGuAD0fEGUOuVRo559DVskFG6BcCRzPz/sx8ArgO2NvXJoFnREQAZwGPASeGWqk0Bu7lopYNEuhbgAer58fKsdrHgF8GjgN3Au/KzE7/D4qIfRExGxGzc3NzqyxZGp35k6ImutozSKDHIsf6P+1vAG4DfgF4KfCxiHjmgm/KPJCZM5k5Mz09vcJSpdFzHbpaNkigHwO2Vc+30h2J1y4Drs+uo8C3gF8aTonS+Jzcy2XB35fS2jdIoN8M7IyIHeVE5yXAob42DwCvBYiIc4AXAfcPs1BpHHo3tnDKRS3auFyDzDwREVcCNwBTwMHMPBIRl5fX9wMfAK6NiDvpTtFclZmPjrBuaSS8sEgtWzbQATLzMHC479j+6vFx4PXDLU0aP+fQ1TKvFJUqrkNXywx0qeKUi1pmoEsV93JRywx0qeJeLmqZgS5VnENXywx0qeKUi1pmoEsVp1zUMgNdqsyvQ59wIdIqGOhSJd1tUQ0z0KWKc+hqmYEuVbzBhVpmoEuV+ZOiJrraY6BLlXSEroYZ6FKld2ML59DVIgNdqjiHrpYZ6FLFOXS1zECXKonLFtUuA12qpPuhq2EGulTxwiK1zECXKm7OpZYZ6FLFEbpaZqBLFW9woZYZ6FLl5IVFncnWIa2GgS5V5vdDd4Su9hjoUqXjskU1zECXKs6hq2UGulRxLxe1zECXKu7lopYZ6FLFdehqmYEuVdzLRS0bKNAjYndE3BsRRyPi6iXavCYibouIIxFx43DLlMbDEbpatnG5BhExBVwDvA44BtwcEYcy8+6qzdnAx4HdmflARDxvRPVKIzW/Dn3ChUirMMgI/ULgaGben5lPANcBe/vavAW4PjMfAMjMR4ZbpjQe8+vQTXS1Z5BA3wI8WD0/Vo7VXgg8OyL+IyJuiYi3LfaDImJfRMxGxOzc3NzqKpZGyHXoatkggR6LHOv/tG8ELgB+A3gD8GcR8cIF35R5IDNnMnNmenp6xcVKo+aVomrZsnPodEfk26rnW4Hji7R5NDN/BPwoIr4KvAS4byhVSmPiXi5q2SAj9JuBnRGxIyLOAC4BDvW1+SJwUURsjIinA68A7hluqdLoOUJXy5YdoWfmiYi4ErgBmAIOZuaRiLi8vL4/M++JiH8G7gA6wCcz865RFi6NgnPoatkgUy5k5mHgcN+x/X3PPwR8aHilSeNXL1vMTCIWO4UkrU1eKSpV6htbOEhXawx0qVJPtZjnao2BLlXqUbnz6GqNgS5V6hA30NUaA12qnDLlYp6rMQa6VHHKRS0z0KXKqVMuEyxEWgUDXap0HKGrYQa6VDllDr1zmobSGmSgSxXn0NUyA12quGxRLTPQpYonRdUyA12q1CHunuhqjYEuVTKTDWWDRUfoao2BLlU6CRs3bCiPTXS1xUCXKp1MpsoQ3UBXawx0qdLpJBtLoJvnao2BLlUyYWrKEbraZKBLlYSTI3RPiqo1BrpUcQ5dLTPQpUon8+QqF9ehqzUGulTpJNUIfcLFSCtkoEuVzKzm0E10tcVAlyqnjNDdPleNMdCliidF1TIDXSoyk0zYOOWFRWqTgS4VvQCfci8XNcpAl4pegHtSVK0y0KWic3KE7rJFtclAl4r+EboXFqk1BrpUpCN0NW6gQI+I3RFxb0QcjYirT9Pu5RHxZET87vBKlMajN0J32aJatWygR8QUcA2wB9gFXBoRu5Zo90HghmEXKY2DJ0XVukFG6BcCRzPz/sx8ArgO2LtIuz8GPg88MsT6pLHpPylqnqs1gwT6FuDB6vmxcuykiNgCvAnYf7ofFBH7ImI2Imbn5uZWWqs0UnlyhN7bbXGS1UgrN0igxyLH+j/qfwVclZlPnu4HZeaBzJzJzJnp6ekBS5TGY+GyRRNdbdk4QJtjwLbq+VbgeF+bGeC6iADYBFwcEScy8x+HUaQ0Ds6hq3WDBPrNwM6I2AF8F7gEeEvdIDN39B5HxLXAlwxztaZ/lYt5rtYsG+iZeSIirqS7emUKOJiZRyLi8vL6aefNpVb0AnyjN4lWowYZoZOZh4HDfccWDfLM/MOnXpY0fr0A3xBeWKQ2eaWoVPQC3Dl0tcpAl4pOpzeH7k2i1SYDXSoWzqFPsBhpFQx0qXAvF7XOQJeKhevQJ1mNtHIGulQs3MvFRFdbDHSpSK8UVeMMdKmYH6GXm0R3JliMtAoGulS4l4taZ6BLhXu5qHUGulS4l4taZ6BLxcJ16JOsRlo5A10q3MtFrTPQpWJ+hO5eLmqTgS4VC9ehT7IaaeUMdKnwnqJqnYEuFb3tcx2hq1UGulS4l4taZ6BLxck5dNehq1EGulT04vvkXi7muRpjoEvFyWWL4QhdbTLQpWLhHPoEi5FWwUCXik7/HLpzLmqMgS4V6V4uapyBLhW9G1q4l4taZaBLRS/AN4Tr0NUmA10qelMsGyLYEE65qD0GulT0RuQbNnRD3SkXtcZAl4pTR+iBca7WGOhSMT+HDhGeFFV7Bgr0iNgdEfdGxNGIuHqR198aEXeUr69FxEuGX6o0Wr0Aj94I3TxXY5YN9IiYAq4B9gC7gEsjYldfs28Bv5qZLwY+ABwYdqHSqGX/SVHPiqoxg4zQLwSOZub9mfkEcB2wt26QmV/LzO+XpzcBW4dbpjR69ZRL96TohAuSVmiQQN8CPFg9P1aOLeUdwJcXeyEi9kXEbETMzs3NDV6lNAb1SVHn0NWiQQI9Fjm26Cc9In6NbqBftdjrmXkgM2cyc2Z6enrwKqUxmJ9Dhw0bwguL1JyNA7Q5Bmyrnm8Fjvc3iogXA58E9mTm94ZTnjQ+WV0p6pSLWjTICP1mYGdE7IiIM4BLgEN1g4g4F7ge+IPMvG/4ZUqjt/BKURNdbVl2hJ6ZJyLiSuAGYAo4mJlHIuLy8vp+4H3Ac4GPR3cfjBOZOTO6sqXhO3UduiN0tWeQKRcy8zBwuO/Y/urxO4F3Drc0abx6AR5lhO4culrjlaJSkfVJUfdyUYMMdKnoXUjkSVG1ykCXivmTou7lojYZ6FLhXi5qnYEuFVmN0F22qBYZ6FLR8cIiNc5Alwr3clHrDHSp6PQtW3QdulpjoEvFgr1cOhMuSFohA10qXLao1hnoUlGfFHUvF7XIQJeK+b1ccC8XNclAl4rMJGL+wiKnXNQaA10qOpls6G7/XC4smnBB0goZ6FLRyW6QQ28/dBNdbTHQpaKTSVQjdPNcrTHQpSKrEbpz6GqRgS4VnU49h26gqz0GulR059C7gR6eFFWDDHSp6JRli+BeLmqTgS5VTk65bHCErvYY6FLRXYfefewculpkoEtFfWFReAs6NchAl4pO0rcO3URXWwx0qcgFUy6TrUdaKQNdKjod+vZyMdHVFgNdKuqTou6HrhYZ6FLhHLpaZ6BLRWayofwX4bJFtchAl4pT90N3ykXtMdClYuFeLia62jJQoEfE7oi4NyKORsTVi7weEfHR8vodEXH+8EuVRmvhXi6TrUdaqWUDPSKmgGuAPcAu4NKI2NXXbA+ws3ztAz4x5Dqlkct02aLatnGANhcCRzPzfoCIuA7YC9xdtdkLfDq7ywJuioizI2JzZj407IJvvG+Ov/jS3cs3lFbooR/8hM3P+jmgG+zf/f6Ped1HbpxwVVqPfv/l23jnRecN/ecOEuhbgAer58eAVwzQZgtwSqBHxD66I3jOPffcldYKwFlnbmTnOWet6nul09l5zlm8euc0AG++YCs/OfHkhCvSerXprDNH8nMHCfRY5Fj/36KDtCEzDwAHAGZmZlb19+wFz382Fzz/gtV8qzSwV71gE696waZJlyGtyCAnRY8B26rnW4Hjq2gjSRqhQQL9ZmBnROyIiDOAS4BDfW0OAW8rq11eCfxgFPPnkqSlLTvlkpknIuJK4AZgCjiYmUci4vLy+n7gMHAxcBR4HLhsdCVLkhYzyBw6mXmYbmjXx/ZXjxO4YrilSZJWwitFJWmdMNAlaZ0w0CVpnTDQJWmdiElt4h8Rc8B3Vvntm4BHh1jOsKzVumDt1mZdK2NdK7Me63p+Zk4v9sLEAv2piIjZzJyZdB391mpdsHZrs66Vsa6V+VmryykXSVonDHRJWidaDfQDky5gCWu1Lli7tVnXyljXyvxM1dXkHLokaaFWR+iSpD4GuiStE80F+nI3rB5jHdsi4t8j4p6IOBIR7yrH3x8R342I28rXxROo7dsRcWf5/bPl2HMi4l8i4hvl32ePuaYXVX1yW0T8MCLePYn+ioiDEfFIRNxVHVuyfyLiT8vn7d6IeMOY6/pQRHy93Hz9CxFxdjm+PSJ+XPXb/iV/8GjqWvJ9m3B/fa6q6dsRcVs5Ps7+WiobRv8Zy8xmvuhu3/tN4DzgDOB2YNeEatkMnF8ePwO4j+5NtN8P/MmE++nbwKa+Y38JXF0eXw18cMLv4/8Az59EfwGvBs4H7lquf8p7ejtwJrCjfP6mxljX64GN5fEHq7q21+0m0F+Lvm+T7q++1z8MvG8C/bVUNoz8M9baCP3kDasz8wmgd8PqscvMhzLz1vL4f4F76N5Hda3aC3yqPP4U8NuTK4XXAt/MzNVeKfyUZOZXgcf6Di/VP3uB6zLz/zLzW3T3/L9wXHVl5lcy80R5ehPdu4GN1RL9tZSJ9ldPRATwe8BnR/G7T+c02TDyz1hrgb7UzagnKiK2Ay8D/rMcurL8iXxw3FMbRQJfiYhbyo25Ac7Jchep8u/zJlBXzyWc+h/apPsLlu6ftfSZ+yPgy9XzHRHx3xFxY0RcNIF6Fnvf1kp/XQQ8nJnfqI6Nvb/6smHkn7HWAn2gm1GPU0ScBXweeHdm/hD4BPCLwEuBh+j+2Tdur8rM84E9wBUR8eoJ1LCo6N7G8I3AP5RDa6G/TmdNfOYi4r3ACeAz5dBDwLmZ+TLgPcDfRcQzx1jSUu/bmugv4FJOHTSMvb8WyYYlmy5ybFV91lqgr6mbUUfE0+i+YZ/JzOsBMvPhzHwyMzvAXzOiPzdPJzOPl38fAb5Qang4IjaXujcDj4y7rmIPcGtmPlxqnHh/FUv1z8Q/cxHxduA3gbdmmXQtf55/rzy+he686wvHVdNp3re10F8bgd8BPtc7Nu7+WiwbGMNnrLVAH+SG1WNR5uj+BrgnMz9SHd9cNXsTcFf/9464rp+PiGf0HtM9qXYX3X56e2n2duCL46yrcsrIadL9VVmqfw4Bl0TEmRGxA9gJ/Ne4ioqI3cBVwBsz8/Hq+HRETJXH55W67h9jXUu9bxPtr+LXga9n5rHegXH211LZwDg+Y+M46zvkM8gX0z1r/E3gvROs41fo/ll0B3Bb+boY+FvgznL8ELB5zHWdR/eM+e3AkV4fAc8F/g34Rvn3ORPos6cD3wOeVR0be3/R/R/KQ8BP6Y6O3nG6/gHeWz5v9wJ7xlzXUbrzq73P2P7S9s3l/b0duBX4rTHXteT7Nsn+KsevBS7vazvO/loqG0b+GfPSf0laJ1qbcpEkLcFAl6R1wkCXpHXCQJekdcJAl6R1wkCXpHXCQJekdeL/Aa4zBiEvDLjEAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# construct initial condition: 1 unit of heat at midpoint. \n",
    "v0 = np.zeros(n)\n",
    "v0[int(n/2)] = 1.0\n",
    "plt.plot(v0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "5ceb6d76",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7fab3a70ab20>"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABgO0lEQVR4nO3dd3hUVfrA8e+ZkkwmfdJ7CIQSBAKE3pFQLNhXLGtBV7Gs234q6666bHN1V13dVVnXhm1d14qKgEhHSgIEAoRACJBMeu/JtPP7YyYhQAIB0gjn8zx5MnPuufeeezN5c3Luue8VUkoURVGUvkvT0w1QFEVRupYK9IqiKH2cCvSKoih9nAr0iqIofZwK9IqiKH2crqcb0JbAwEAZGxvb081QFEW5aOzcubNUShnU1rJeGehjY2NJTU3t6WYoiqJcNIQQx9tbpoZuFEVR+jgV6BVFUfo4FegVRVH6uF45Rq8oyqXLarViNptpbGzs6ab0SgaDgcjISPR6fYfXUYFeUZRexWw24+3tTWxsLEKInm5OryKlpKysDLPZTL9+/Tq8nhq6URSlV2lsbCQgIEAF+TYIIQgICDjn/3ZUoFcUpddRQb5953NuVKBXlA5KOVbOntzKnm6GopwzFegVpQPyKhu4480d/PjN7ZTWNvV0c5QuVFlZyauvvtop23rnnXcICgoiMTGRxMRE3njjjZZly5YtIz4+nvj4eJYtW9Yp+2uPCvSK0gFLlu9HImmw2vnzioyebo7ShToz0APcfPPNpKWlkZaWxr333gtAeXk5S5YsYfv27ezYsYMlS5ZQUVHRafs8lQr0inIWKcfKWX2giJ9dPpD7p/bns115HMiv7ulmKV1k8eLFHDlyhMTERB599NEu2ceqVatITk7GZDLh7+9PcnIyK1eu7JJ9gZpeqShnte1IGQC3jY+mvsnOP9dlsf1oGQnhPj3csr5vyVf7O/2PakK4D09fPbTd5X/5y1/Yt28faWlpbS6fMmUKNTU1p5X/7W9/Y9asWaeVf/rpp2zcuJGBAwfy4osvEhUVRV5eHlFRUS11IiMjycvLO/eD6SAV6BXlLPaYK+kf5ImPQY+PQU+Ij7u6KHsJ27RpU4frXn311dxyyy24u7uzdOlS7rzzTtauXUtbz+ruyplGKtAryhlIKUnLrWLqwMCWshGRfuwxV/Vgqy4dZ+p595Rz6dEHBAS0vP7JT37C448/Djh78OvXr29ZZjabmT59epe0F1SgV5Qzyq9qpLS2icQov5ayEVF+rD5QRGW9BT+jW881TukS3t7ebQbyZufSoy8oKCAsLAyA5cuXM2TIEADmzJnDE0880XIBdvXq1TzzzDMX0OozU4FeUc6geYhmRKRfS1lz0N9rrmLqwDaf86BcxAICApg0aRKXXXYZ8+bN469//et5b+vll19m+fLl6HQ6TCYT77zzDgAmk4knn3ySMWPGAPDUU09hMpk6o/lt6lCgF0LMBV4CtMAbUsq/nLJcuJZfAdQDd0kpd7mW+QFvAJcBElgopdzaWQegKF1pT24lbloNg8O8W8qGRfoihHOZCvR904cfftgp23nmmWfa7akvXLiQhQsXdsp+zuas0yuFEFrgFWAekADcIoRIOKXaPCDe9XUf8FqrZS8BK6WUg4ERgJqErFw00nIrSQj3wV2nbSnzMejpH+TFHnNlzzVMUc5BR+bRjwWypJTZUkoL8BFwzSl1rgHelU7bAD8hRJgQwgeYCrwJIKW0SCkrO6/5itK1jpTUMTDE67TyQSHeHCmp64EWKcq560igjwByW703u8o6UicOKAHeFkLsFkK8IYTwbGsnQoj7hBCpQojUkpKSDh+AonSVeouN0tomYgJO/8hGmYyYK+qxO06fJqcovU1HAn1bkztP/XS3V0cHjAJek1KOBOqAxW3tREr5upQySUqZFBSkxj2Vnpdb3gA4g/qpYgKMWO2Swmr1cAyl9+tIoDcDUa3eRwL5HaxjBsxSyu2u8k9wBn5F6fVyyusBiG4j0DeX5ZTVd2ubFOV8dCTQpwDxQoh+Qgg3YAGw/JQ6y4E7hNN4oEpKWSClLARyhRCDXPUuBw50VuMVpSs1B/qYMwT63HIV6JXe76yBXkppAx4GVuGcMfOxlHK/EGKREGKRq9oKIBvIAv4NPNhqEz8FPhBC7AUSgT93XvMVpevklNXh7a7Dz3j6sznDfA1oNYLj5eqCbF/TmdkrN27cyKhRo9DpdHzyyScnLWsvTfHRo0cZN24c8fHx3HzzzVgslgtuR4eyV0opV0gpB0op+0sp/+QqWyqlXOp6LaWUD7mWD5NSprZaN8019j5cSnmtlLLrcnEqSifKKa8nymRsMweJTqshws+DHNc4vtJ3dGagj46O5p133uHWW289qfxMaYoff/xxfvGLX3D48GH8/f158803L7gdKk2xorQjp7y+zfH5ZtEmY8vwjtJ3dGaa4tjYWIYPH45Gc3KobS9NsZSStWvXcuONNwJw55138sUXX1xQG0ClQFCUNjkcktyKBi4fEtJunegAIyv3FXZjqy5B3y6GwvTO3WboMJj3l3YXd3aa4ra0l6a4rKwMPz8/dDrdSeUXSgV6RWlDcU0TFpujzamVzaJNRsrrLNQ0WvE2nD6Or/RN55LUrD3tpSnuqvTFKtArShvONLWy2YmZNw0khKtA3yXO0PPuKZ3Ro28vTXFgYCCVlZXYbDZ0Oh1ms5nw8PALbrMK9IrSBnOFM9BH+nu0WyfK39hSVz1tqu/ozDTF7WkvTbEQghkzZvDJJ5+wYMECli1bxjXXnJpx5typi7GK0oai6iYAQn0M7dYJ8XV31q1p6pY2Kd2jdZriC70Ym5KSQmRkJP/73/+4//77GTrU+SCV1mmKx4wZc1Ka4meffZYXXniBAQMGUFZWxj333HPBx6R69IrShqLqRrzddXi6t/8rEuDpjlYjKKpSaRD6ms5KUzxmzBjMZnOby9pLUxwXF8eOHTs6Zf/NVI9eUdpQVN1IsI/7GetoNYIgL3eKVL4bpZdTgV5R2lBU3UjIGYZtmoX4uKuhG6XXU4FeUdpQVN10xvH5ZsE+BopVj17p5VSgV5RTOByS4ppGgjsQ6EN9DCpVsdLrqUCvKKeoqLdgtUtCzjJGD86hm8p6K41Weze0TFHOjwr0inKKjkytbNbc6y9R4/RKL6YCvaKcoqjGORTTkaGb5gu2auZN39FdaYq1Wi2JiYkkJiYyf/78lvIeS1OsKJeS5ourHRm6CW0J9KpH31d0R5piAA8PD9LS0khLS2P58hPPclJpihWlGxRWOYN2sHfHplcC6oJsH9IdaYrbo9IUK0o3KappJMDTDTfd2X85fT30uOk0aoplF3l2x7McLD/YqdscbBrM42Mfb3d5d6QpBmhsbCQpKQmdTsfixYu59tprVZpiRekuxdUdm1oJzhSyIT7q7thLSWckNQPIyckhPDyc7OxsZs6cybBhw/DxOT05nkpTrChdoKi6qUPj881CvA1qjL6LnKnn3VM6q0ffnH44Li6O6dOns3v3bm644QaVplhRukNxTSNDwrw7XD/Yx52Dhe2ntVUuLt2RpriiogKj0Yi7uzulpaVs2bKFxx57TKUpVpTuIKWkvM6CybPjPfpAL3fKai98CpzSO3RHmuKMjAySkpIYMWIEM2bMYPHixSQkJAAqTbGidLmaJhtWuyTQy63D6wR6uVPVYMVic3ToAq7S+3V1muKJEyeSnt72s3B7LE2xEGKuECJTCJElhFjcxnIhhHjZtXyvEGJUq2XHhBDpQog0IURqZzZeUTpbc8/c5NnxQB/g+qNQXqd69UrvdNYevRBCC7wCJANmIEUIsVxKeaBVtXlAvOtrHPCa63uzGVLK0k5rtaJ0kfI650XVcwn0gV7OYZ7S2iZCfTs2W0dRulNHevRjgSwpZbaU0gJ8BJx6deAa4F3ptA3wE0KEdXJbFaXLNffoA85pjN75R6G0Vs28UXqnjgT6CCC31Xuzq6yjdSSwWgixUwhxX3s7EULcJ4RIFUKklpSUdKBZitL5modfAs5xjB6gVF2QVXqpjgT6tmbry3OoM0lKOQrn8M5DQoipbe1ESvm6lDJJSpkUFBTUgWYpSucrqzufMXpnoC9TPXqll+pIoDcDUa3eRwL5Ha0jpWz+Xgx8jnMoSFF6pbJaC55uWgx6bYfXcdbXqKEbpdfqSKBPAeKFEP2EEG7AAmD5KXWWA3e4Zt+MB6qklAVCCE8hhDeAEMITmA3s68T2K0qnKq9rwnQOwzbgvEU9wFPNpe8rOjN75QsvvEBCQgLDhw/n8ssv5/jx4y3Lli1bRnx8PPHx8SxbtqylvEfSFEspbcDDwCogA/hYSrlfCLFICLHIVW0FkA1kAf8GHnSVhwCbhRB7gB3AN1LKlRfcakXpImXneLNUs0Bvd0pUj75P6MxAP3LkSFJTU9m7dy833ngjjz32GADl5eUsWbKE7du3s2PHDpYsWUJFRQXQg2mKpZQrpJQDpZT9pZR/cpUtlVIudb2WUsqHXMuHSSlTXeXZUsoRrq+hzesqSm9VXmch8BzG55sFerqpHn0f0ZlpimfMmIHRaARg/PjxLTdPrVq1iuTkZEwmE/7+/iQnJ7Ny5UqVplhRukNZrYWEsNMzCJ5NoJc76XlVXdCiS1vhn/9MU0bnpil2HzKY0CeeaHd5V6UpfvPNN5k3bx4AeXl5REWduKzZnI5YpSlWlC7WkufmHMfowTkds7zOgsMh0WguPK2s0nudT1Kz999/n9TUVDZs2AA4P2unEkK0W36hVKBXFJfaJhsWu4OA8xm68XLH5pBUNVjxP4/1lbadqefdU861R79mzRr+9Kc/sWHDBtzdndd/IiMjWb9+fUsds9nM9OnTCQwMVGmKFaUrtdwsdR4XY5tvsCqra1KB/iLXmWmKd+/ezf3338/KlSsJDg5uKZ8zZw5PPPFEywXY1atX88wzz6g0xYrS1VpuljqPoZsg101TJTXqguzFrjPTFD/66KPU1tZy0003kZiYyPz58wEwmUw8+eSTjBkzhjFjxvDUU09hMpkAlaZYUbrUiTw35zNGfyKxmXLx66w0xWvWrGl32cKFC1m4cOFp5T2WplhRLgXnk7myWfM6FfWqR6/0PirQK4pL2QWM0fsb9c5tqLn0Si+kAr2iuJTXWjC6afFw63iem2Y6rQY/o149fETplVSgVxQXZ/qDMwzbNFRCY3W7i01GNxXolV5JBXpFcSmrs7R9IdZugw3PwQtD4IUE2PISOBynVTN5ulFWpy7GKr2PmnWjKC7ldU0t0yRPsuXvsO5PMGQ+2C3w3VOg84BxJz9Hx+TpxrGyuu5prKKcA9WjVxSX8lpLyzTJFlVm2PQ8DLkabn4PbvkI4qbDuj9C3cmPQXamQbB2X4OVLtFdaYq1Wi2JiYknza+HHkpTrCiXAill20M3a34H0gGzXYlXhYB5z4GlDtb+8aSqJk83Kuqd+W6Ui1d3pCkG8PDwIC0tjbS0NJYvP/GIjx5LU6wofV2dxU6TzXHyxdi6Utj/OYy5F/xjTpQHDYLEW2Hvf6GptqXY5OmO3SGpblS9+otZd6Qpbo9KU6woXai8to1nxe77DBw2Z1A/VeJtsOtdyPgKEm9xreuaS19nwc+o8t10hk0fH6I0t/bsFc9BYJQXU340sN3l3ZGmGKCxsZGkpCR0Oh2LFy/m2muvVWmKFaUrNc+WCWid52bvRxAyDEKGnr5C1Djwj4U9/2kV6J3j++V1Fvqr59v3WZ2RphggJyeH8PBwsrOzmTlzJsOGDcPH5/RnIag0xYrSSU7LXFmaBXk7YfYf215BCBh+s3PaZVUe+Ea0jO+ru2M7z5l63j2lM9IUAy3ph+Pi4pg+fTq7d+/mhhtuUGmKFaWrlJ06dJP5jfP70OvbX+myG2HDs3BoJYy5R+W76SO6I01xRUUFRqMRd3d3SktL2bJlC4899phKU6woXaklz03z0M2RdRA0BHwj2l8pMB58IiF7HXDij4S6O/bi1h1pijMyMkhKSmLEiBHMmDGDxYsXk5CQAKg0xYrSZcrrmjDoNRjddGBthJytkHR6CtmTCAH9pzsvyDrsGPRaPN20auimD+jqNMUTJ04kPT29zWUqTbGidBHnHHrX+GnOVrA1QtyMs68YNwMaqyB/NwD+nm4t6Y4VpbfoUKAXQswVQmQKIbKEEIvbWC6EEC+7lu8VQow6ZblWCLFbCPF1ZzVcUTpTeZ3lxLBN9jrQ6CF20tlXjJt+Yh2cDy0pU0M3Si9z1kAvhNACrwDzgATgFiFEwinV5gHxrq/7gNdOWf4zIOOCW6soXaSstlXmyuz1zumTbp5nX9EzEEKHw5H1gHOcXo3RK71NR3r0Y4EsKWW2lNICfAScehn4GuBd6bQN8BNChAEIISKBK4E3OrHditKpyptTFDfVQmF6x3rzzWInO6di2q2YPN2pUIFe6WU6EugjgNxW782uso7W+TvwGHB6XtdWhBD3CSFShRCpJSUlHWiWonSesrom5zz4/N3O3DaRYzq+cmQS2BqgaB8BXs6hGylVvhul9+hIoG/rtqxTP8Vt1hFCXAUUSyl3nm0nUsrXpZRJUsqkoCB1W6HSfeotNhqtDmfmSrNrtkPE6I5voPmPgjkVk6cbTTYH9RZ75zdUUc5TRwK9GYhq9T4SyO9gnUnAfCHEMZxDPjOFEO+fd2sVpQucdLOUORUC4sFo6vgGfKPAKxRyd2Ayqrn0F7vOzF65dOlShg0bRmJiIpMnT+bAgQMty5YtW0Z8fDzx8fEsW7aspbyn0hSnAPFCiH5CCDdgAbD8lDrLgTtcs2/GA1VSygIp5a+llJFSyljXemullLdfcKsVpRO13Cxl1EPujnMbtgHnfPrIJDCntFzQVTNvLl6dGehvvfVW0tPTSUtL47HHHuOXv/wlAOXl5SxZsoTt27ezY8cOlixZQkVFBdBDaYqllDbgYWAVzpkzH0sp9wshFgkhFrmqrQCygSzg38CDF9wyRekmzfPeQx1FUF8KUecY6AGixkLFUYK1zmfKqguyF6/OTFPcOklZXV1dS4KyVatWkZycjMlkwt/fn+TkZFauXNmzaYqllCtwBvPWZUtbvZbAQ2fZxnpg/Tm3UFG6WPPQTWiN607FiKRz34jrv4DQmn2AmkvfWda98zrFx7M7dZvBMXHMuOu+dpd3dpriV155hRdeeAGLxcLatWsByMvLIyrqxGh3czpilaZYUbpI83i6T8UB0Bkg+NTbRDogdDgg8K08ACSqu2P7sHNNU/zQQw/x0EMP8eGHH/LHP/6RZcuWtTkrSwjRbvmFUoFeueSV1Vlw12nQlaQ7g7z2PH4t3L0gYABuJftw045SPfpOcqaed0853wePLFiwgAceeABw9tTXr1/fssxsNjN9+nQCAwNVmmJF6QpltRYCjHpEwR4Yet35byhsOCJ3B/6e+pYnVikXn85MU3z48GHi4+MB+Oabb1pez5kzhyeeeKLlAuzq1at55plnVJpiRekq5XVNDDZWOpOThQ4//w2FDoeqXGI8mtT0yotYZ6Yp/uc//8nQoUNJTEzkhRdeaJlGaTKZePLJJxkzZgxjxozhqaeewmRyTulVaYoVpQuU11m4XJfjfBM2os06hXWFrM9dj0ZomB41nWBj8OmVwpx/JEbqc9lR79c1jVW6RWelKX7ppZfaXbZw4UIWLjw9FXZXpClWgV655JXWWhjsdRSEps0Lscv2L+P51OeRrhvC/7z9zzwx7gl+NOhHJ1cMdf6RGCKOsrJuUJe3W1E6SgV65ZJXXmch1v0IBA4EN+NJy15Le41X97xKckwyPx35UxzSwQs7X+AP2/5Ag62BO4feeaKyZwD4RBBny1Zj9EqvosbolUtag8VOg9VOeGPWaePzqYWpvLrnVeb3n89fp/6Vfr796O/Xn7/P+DvJMcm8uPNF9pftP3mDocOJbDxMTZONJpvKd6P0DirQK5e0sromfKjDq6kIQk4M2zTYGnj6h6eJ9IrkN+N+g1ajbVmm1+h5esLTmAwmfrv5t1jt1hMbDBmKX2MOblipqLOiKL2BCvTKJa28zsIA4brzMGhwS/kHGR+QU5PDkolLMOqNp63n6+7Lk+OfJKsyi08Of3JiQfAQNNJOP1FAmbppSuklVKBXLmlltRYGaE4O9PXWet7d/y6TIyYzNmxsu+tOj5rOyOCRvL3v7RO9+uAhAAwSZtWjV3oNFeiVS1pZnYWBwoxDZwC/GAA+PfwpFU0V3D/8/jOuK4TgvuH3UVBXwNfZrschBwxACi3xGrPq0V+kuitNsVarJTExkcTERObPn99S3lNpihWlzyqvayJe5CEDB4FGg91h590D7zImdAyJwYlnXX9S+CQSAhJ4e//bzjwlOnccpv4MErnqpqmLVHekKQbw8PAgLS2NtLQ0li8/kfm9R9IUK0pfVlZnIV6ThybYOe/9h/wfKKwr5JbBt3RofSEECwYt4GjVUdJK0gDQBA9hoMasAv1FqjvSFLenR9MUK0pfVVtVQbgoaxmf/+zwZ5gMJqZHTu/wNubEzuHZlGf59NCnjAweiQhJIDpjOdU11V3U6ktH5VdHsOTXdeo23cI98bu6f7vLuyNNMUBjYyNJSUnodDoWL17Mtddeq9IUK0pX8KjKcr4IGkJZQxnrc9dz25Db0Gv1Hd6GUW9kbuxcVhxdweKxi/EKGowGiXtFFnAeDzFRerXOSFMMkJOTQ3h4ONnZ2cycOZNhw4ad9B9AM5WmWFEukE/tEeeLoEF8e/RbbNLGdfHnnsHy+vjr+fTwp3x3/DuuCx7q3HZNVmc29ZJ0pp53T+mMNMVAS/rhuLg4pk+fzu7du7nhhhtUmmJF6WxBDUexCDfc/GNZvWMJA/0H0t/v3IPLsMBhRHhFsPr4aq6LuxoreoIajnRBi5Wu1h1piisqKjAajbi7u1NaWsqWLVt47LHHVJpiRekKEdbjlBtiKGooZXfxbmbHzD6v7QghmB0zm20F26iy1VFqiCbccryTW6t0h+5IU5yRkUFSUhIjRoxgxowZLF68mIQE553ZKk2xonSiRqudOMxUeyeRkrMGgNmxbQd6abPRlJ2N0Ghw69cPodWeVmd27Gze3v8263LXMchrAHENO7E7JFrNhY+xKt2rq9MUT5w4kfT09DaXqTTFitKJyivKiRSl7PWLZ/Wx1cT7x9PPt99JdRyNjZT9+w3K334bR309ABpvbwJ+8hMC7roT4ebWUndowFDCPcNZfWw1UX7xDCldRXlFGaaAwG49LkU5lRq6US5Z9XkZAFSZotldvJtZ0SdfSLOVlnLsph9R+soreE6ZQvhzzxL2l2cwjhlDyQsvcOzW27BXVbXUF0JweczlbCvYRq1/HAB15n3dd0CK0o4OBXohxFwhRKYQIksIsbiN5UII8bJr+V4hxChXuUEIsUMIsUcIsV8IsaSzD0BRzpet0JlieJ+hCYlkWuS0E8vKyzl+111YzGaiXv8XkS/9Hd/58/G79lqiXn2FiJdeoikzk5y7F2KvPjFfflrkNKwOK5meDgCshQdQlJ521kAvhNACrwDzgATgFiHEqY/hmQfEu77uA15zlTcBM6WUI4BEYK4QYnznNF1RLowozaRJ6ki3HCHAEMCQAGdCMulwkP/Y41hzcon611K8pk49bV2fObOJfOWfNB46RMFvn3SmPwBGBY/CS+/FAdtRGqQbouRgtx6TorSlIz36sUCWlDJbSmkBPgJOne9zDfCudNoG+Akhwlzva1119K4v2VmNV5QL4VGZxWEZTlrpDqZGTkUjnL8O5W+/Td3mzYQ88Ws8x7afvdJr6lSCf/FzalavpvK//wVAr9UzIXwCaWXbOCTDca841C3Hoihn0pFAHwHktnpvdpV1qI4QQiuESAOKge+klNvb2okQ4j4hRKoQIrWkpKSDzVeU8+dTm8U6QyA11hqmRjp77ZbcXEpeehnv5GT8br75rNsw3X03npMmUfzcX7EWFQMwNXIqZY0lbNKH4FNzuEuPQVE6oiOBvq25Yaf2ytutI6W0SykTgUhgrBDisrZ2IqV8XUqZJKVMCgoK6kCzFOUCWOrwayog1csdndAxIXwCAEXPPgs6HSG//W2Hbj0XGg2hTz+FtFopeeF5ACZHTAYg1csdL0spNFR22WEona8zs1c2++STTxBCkJqa2lK2bNky4uPjiY+Pb5lfDz2XptgMRLV6Hwnkn2sdKWUlsB6Ye66NVJROV+ocUjlqbGB40HA89Z7Ubd1K7ZrvCVy0CH1IcIc35RYdjenuu6n6cjkNe/YQ6BHIQP+BHPdsPGlfysWhswN9TU0NL7/8MuPGjWspKy8vZ8mSJWzfvp0dO3awZMkSKioqgJ5LU5wCxAsh+gkh3IAFwPJT6iwH7nDNvhkPVEkpC4QQQUIIPwAhhAcwC1BXp5SeV3yQCo2GCl05E8InIKWk5OV/oAsLw3TXnee8ucD770Pr70/JP/4JwLiwcZS5ldAoBJRkdnbrlS7UmWmKAZ588kkee+wxDAZDS9mqVatITk7GZDLh7+9PcnIyK1eu7Lk0xVJKmxDiYWAVoAXeklLuF0Isci1fCqwArgCygHrgbtfqYcAy18wdDfCxlPLrC261olyokoNs8TCCgInhE6nfupWG3bsJffopNK1uguoojacnpoV3U/L8CzTs2cP4sPG8d+A9Utw9maJm3py3b7/9lsLCwk7dZmhoKPPmzWt3eWemKd69eze5ublcddVV/O1vf2spz8vLIyrqxCBIczriHk1TLKVcgTOYty5b2uq1BB5qY729wMgLbKOidL6STNYY/NELT4YGDCX3F3eiCwnB94YbTqomHZLcjHIOpxRRUVSPEOAf6snAsSFEDPI/aRzfdOutlL/5FqWvvsbofz6PQMsaDxNTVI++T+loUjOHw8EvfvEL3nnnndOWNU/HbU0I0W75hVIpEJRLkr34ADu9dEQahmHZn0FD6k5Cfr34pN58RWEd3711gJKcGgyeegIivQBJdloJGT8UEBrny6y7h+AbZAScvXr/O35M6cv/IDi3iED9AHYYspGlmW3OVlDO7kw9757S0R59TU0N+/btY/r06QAUFhYyf/58li9fTmRkJOvXr2+pazabmT59OoGBgSpNsaJ0Cks9+bV5VPqHM953JOXvvYvG0/Ok3nx2WgnfvbkfnZuWy+8cQvyYELQ65yUtm8VO5vZCfvjsCP/9Ywpz77+M6IQAAPxvvpmy15ZS/t57xI0cyQ5LJjVFufhY6sDNs0cOVzk3nZWm2NfXl9LS0pb306dP529/+xtJSUnExcXxxBNPtFyAXb16Nc8884xKU6wonabsMKkGdwBG6+Ko/nYlvjdcj9bLC4Aju4pZ9fo+TBFeLHhqLIMnhLUEeQCdm5ahUyJY8ORYfAI9WPFqOsf3lTmXBQTgc/XVVH3xJYnuQ5ACUgwGKFXz6S8WnZmmuD0mk4knn3ySMWPGMGbMGJ566ilMJhOg0hQrSucoPkiqwR1hM5CwfT9YrZhuuw2AwuwqVr+1n+BYb67+aSJuHu3/inibDFz7i5F8+dJuVr6ezg2PJREY6YXpjh9T9dlnjNubx78CdWw3GLi8JBPCE7vpAJUL1VlpiltrPVQDsHDhQhYuXHhava5IU6x69Mqlp+QgqQYD9rp+uH2/EuOE8bjFxFBX1cS3/0rHy8+dKx8cccYg38zgpeeqh0fg7qFjxWt7aay1Yhg8GMOwYXitX429LpbtHgZQM2+UHqQCvXLJKShOJ1+vI/6oLzLPjN911znnL7+bgaXexhUPDMfg1fGHg3v6ujNv0XDqKpvY8B/nDBvf666FI1lE5QSR7aanuGR/Fx2NopydCvTKJSe1yjlePudQBRovL7yTk8nYUkDO/nIm3jCAgAiv09Zx2O3UlJVSU16Kw2E/bXlIPx/GXt2PrJ3FHE4twvfKKxFublye4byot71K3R2r9Bw1Rq9cWqwNpNpr8G/0YcLxw/hcdw0NFg1bPjlMxCB/Lpt6Il+fdDjI3LqJ/RvXkpO+B4fdBoBWrydmWCKXTU9mwNgJLfOcRyZHc3RPKRs/OkTUkvF4z7qc6Ws28KVVS6qtmqttTaBz75HDVi5tKtArl5ayLFINbsxM98LdVobv9dex9fMj2GwOpt86COF6vmtRdhar/vUyJcey8Q0OYeTcK/EPiwQkZXm5ZKVsY/kLfyZ0wEDm3P8IgdGxaLQapt82iI//lMKOr48y6rrrqV7xLdMPBZEysBHKjkDIqY9yUJSupwK9ckkpzkshR6/nwb12ygPC8PPpx8Ftuxg1Nwa/EOeNT3u+W8G6d17Hw9ePKx55lMETpiA0J49yTr/jXg5sXMemD9/hg9/8iln3PsjQaZcTGOnNZVMj2LfeTMKvR1Pt5c+YvZLlQ/UU5m0nVAV6pQeoMXrlkpKav5XQcklcQSXHk2aw9fMjGH3cGD03BiklWz/5D2veeJXoy0Zwx7MvM2TStNOCPIBGo+Wy6bO447l/EBY/iJWvvkjq158DMHZ+HG4eOrZ9dYwjI6cSf7wM/xpJSv627j5c5Tx0V5pirVZLYmIiiYmJzJ8/v6W8p9IUK0qfkVp5iMv3ObAjqBoyk4KsKsZc1Q83g47tn3/MD//7gKHTLufax57Cw9vnrNvz9PPnhieWED9uIhvee5NdK77E4Kln1JwYjqeXUTJsFhopuXyfZGeluiB7MeiONMUAHh4epKWlkZaWxvLlJxIC91SaYkXpM1Kt5Uw+KEgPjIMciU+QB0MmhbF/w/ds+e97DJkygzmLfoZGqz1pPUejDWtJPdbSBhxNtpOWaXV6rvrZ4wwYM551777Boe1bGDYjEk9fN3TFHmT6RTLtoCDFWtadh6qcp+5IU9yeHktTrCh9RWlNPpYKSUCZg69HX05UuYWxdydQcjSL717/B9GXDWfOokdahmps5Y3Ubi+gMaMMW3HDSdvShxoxDAnAc3wYOl93NFotVzzyKP/7w2/49pUXuO1PLzB6XiwbPzpESuxMbk97l6YqKKzOI9Tn1CdxKu05dOgP1NRmdOo2vb2GMHDgk+0u7440xQCNjY0kJSWh0+lYvHgx1157bc+mKVaUvmDnkRVMyJA4BHj4JODu60bUYA8++M0TGP38uerni9Hq9DjqrVStPk7d9gIQ4N7fD+PIELR+zqmR9vJGmrIrqVmfS80GM14Tw/GZFY3e4M78Xz7Be48/wvIXnmHB7/7KluVaHGIEAOMPSlKPruSqEReeu0TpOZ2RphggJyeH8PBwsrOzmTlzJsOGDcPH5/ThQpWmWFHOQWr+D0w46KBo6OUESj2DZkTw/dtLqS0vZ8Hvn8XD24fGI5WUf3QQR60Vz/FheE+LQufXxtz3y6OxlTdSsy6X2i15NKSXYrptMF7RJq762WP87w+/ZfNHbxM9bjb2dflUDhzDlAMpbM37QQX6c3CmnndP6Yw0xUlJSS3ph+Pi4pg+fTq7d+/mhhtu6JI0xWqMXrlkHM/KIKoUssOSqRYOvNyOcmjrJibedCthAwZRuzWf0jfS0Rh0BD88Ev9rBrQd5F10JgP+N8QT9MAI0EDJ0r3UpRYRNXQ4Y+ZfT/ra1QQGFdOE5HC/+UQXw7EjKhVCb9eRNMXNF1Fbf506bNOcpvjYsWMcO3aM8ePHtwT5iooKmpqaACgtLWXLli0kJCSclKYYUGmKFeVcVDRWELqvmhqvSJqs3mQY6vjh/X8TNmAQSVdfT9V3x6n88giGIQEEP5yIW6s0CLW1tRw5coT09HTS09PJzs6mrq6uZbl7tA8hj4zCvb8vFZ8comajmQk33UZgVAy7/vsGe/X11DSZqPcIIia9hqK6op44BUoHdUea4oyMDJKSkhgxYgQzZsxg8eLFJCQ477HoijTFoq1HV/W0pKQk2Xq+qaJcqO+Pf4/lxw9TF3k3BX5JHOE7oiv2c8ezL+N2WEP1d8cxJoXgf108Qiuora1l586dpKenn/TwiNaCg4MZPnw4o0aNwmg0Im0Oyj/OpGFvKb5XxVEbUsuHv/0Vh02juExMJ7xiK57m97Et+xtXxl3ZvSfgIpKRkcGQIUN6uhm9WlvnSAixU0qZ1FZ9NUavXBIydn/H5EpftiaMpNhYTGRuOqOuvAZDvhuV3x3BODIY/+vjsVgtbPh+Azt27MBmsxETE8PIkSMJDw/Hy/VgkpqaGvLz88nMzGTNmjVs3LiR8ePHM3nyZEw3D6bcnkHV19n4/2ggl82YjWP9GspjEtE6xjJx3xes3blGBXqlW6lAr1wSHOs2Y46Yil0KNKXfYXf3YnTiVVR+lIVhsAn/Gwdy9NhRvvjiC6qrqxkxYgRTpkwhMDDwtG0FBQURFxfH5MmTKS4uZsOGDWzcuJG9e/dy/fXXE3XLYErf2kfFp4cZf+uNpG3eiL3sO+zaG8kLm4R27Va4qQdOgnLJUoFe6fOqmqrot7cGc//JmAIO0nSkCP2426n+7Cj6YE/8Fwxkw6YNrF+/nsDAQO655x6ioqIAkNJOQ0MOdXVHsNmqANDp/fA0DsDDI4rg4GBuuukmkpKSWL58OW+//TazZs1i/G1jKHltL3Wf5cDQK/BJ+wy/2MOYo6YwcN93FNcXE2wM7snTolxCOhTohRBzgZcALfCGlPIvpywXruVXAPXAXVLKXUKIKOBdIBRwAK9LKV/qxPYrylnt3bMaDzEWm1ZPWf4migyR3NkQAwJ8bx3I519/yb59+xgxYgRXXnkler2OsvLNFBZ8RmnZBmy2yja3q9ebCAyYTljYDcTGjuP+++9n+fLlfPfdd5SVlTH71umUvZbOLH0C/3PbhKFkE1J/Bx6MYE/KNyRPu7t7T4RyyTproBdCaIFXgGTADKQIIZZLKQ+0qjYPiHd9jQNec323Ab9yBX1vYKcQ4rtT1lWULlX8zVfUhCejt6yltqGO/tE/xrfGhu+PB/Ppd19y+PBhZs2axcSJEykpXUV29t+pr89Cp/MjMHAG/n7j8fSKx03vj5QSq7Wc2rpDVFRso6T0OwoKP8PLawhxcb/gxhtvZP369WzcuJH6+nquuGoK8vNsAoKvptH8FkbjFnIjpuO2aiWoQK90k4706McCWVLKbAAhxEfANUDrYH0N8K50TuHZJoTwE0KESSkLgAIAKWWNECIDiDhlXUXpUrpd9dSG+WGrySQ+YS6jGnypHOrH5vTvOXz4MFdddRUJCQHsTrudysrtGI0DSEh4npDgeWg0J+bR26VEAEZjDL6+I4kIvxm7vZGioq85dvxV9u69D5NpChMm/AGj0cjKlSuRUnJZzEB+dDyQ3YOmczxrE5U+44ndVd1zJ0S55HRkHn0EkNvqvdlVdk51hBCxwEhg+zm3UlHOU+Wxw1h0Y7E1bEYn3BjqSCQbGyns4+DBg8yZM4ewsMNs33EVtbUZDB70R8aPW4Gb6Sr+W1TLwweOM3PHQfpv3EvE+j2Er9/DgI17mZWSyc8ycvi0pB5j0LWMH7eKgfFPUlW1i+07riQmtpB58+aRmZnJTu9s8nEwXDMeHFrsjduwuE+kMDOtp0+P0obOzF75zjvvEBQU1JKO+I033mhZtmzZMuLj44mPj2fZsmUt5T2VpritRAunTr4/Yx0hhBfwKfBzKWWbXRkhxH1CiFQhRGpJSUkHmqUoZ7f/k/9SZIrFbtnPjPgb0FskH+gOcvxIJrNmzcDX70sOZv4WP98kxo1dyXHj1dy97zgjftjHLw7msq68hjB3N24NM/F/saH8KjaEm0NNBLnpWF1axSMZOQzbsp/7Dpgp9LmJ8eNW4uU1hAMHfoWf/yqmTJnEsaz9vKnNQt8omTbgauxNuykIHMKhzz7u6dOjtKGz0xTffPPNLXfP3nvvvQCUl5ezZMkStm/fzo4dO1iyZAkVFRVA16Qp7sjQjRmIavU+EsjvaB0hhB5nkP9ASvlZezuRUr4OvA7OG6Y60C5FOauinY1YdTuI8hyAb1MkqwOLCastYsKERDy93qagYBOxMQ9SHXg/t2QUsLWyAJNey/1RwVwX7MdQL492k0o5pGRPTQNfFFXw38Jyvi6pYrq/N7+Jf4OoopfIzX2bwKB8hg2bDukH2WQKZGrJAILcwymz7KJs94X31JTO1zpNcXJyMn/96187fR+rVq0iOTkZk8kEQHJyMitXrmTBggWsXbuWDz/8EHCmKf7d737HAw88cEH760igTwHihRD9gDxgAXDrKXWWAw+7xu/HAVVSygLXbJw3gQwp5QsX1FJFOUdNefmU6GPQWL5lbNhCSgIbOV67jzqDG35+/6ai4jBR8c/yev043t2VRYBex5/iI7g1LAAP7dn/2dUIwUgfIyN9jDzWL5T38st46XgRc3cd4ZHoe7hxQCTHsv5IVHQxmw8mkVm7hwSfcUyyz2Z57vuUuV+BxWzGLTKyG87GxenJw2b21TacveI5uMzLgz/Et3/OOzNNMcCnn37Kxo0bGThwIC+++CJRUVHk5eW1TOGFE+mIeyxNsZTSJoR4GFiFc3rlW1LK/UKIRa7lS4EVOKdWZuGcXtk8nWAS8GMgXQiR5ip7Qkq54oJbrihnsf+j9dQ7DjLCNIUGjZ7Vtl1Ig2D8iBU0NFbjFv8Wd5kDONZQxr2RgTzeLwwvnfbsG26Dp07LouhgFoSZeCorjxePF7HKczxLBvyLuuyfMjaxnMy9M1mr3cd87UgS/CZyoC6b9I++Y/T/qdk3F5OOpikGuPrqq7nllltwd3dn6dKl3Hnnnaxdu5a2Us8IIdotv1AdmkfvCswrTilb2uq1BB5qY73NtD1+ryhdbl96Cb6ikv7eV7HCcytC08SQhLV4uZWRGfEefz6ix1/n4LORAxjv53X2DXaAn17Hy0NiuDLIj0czc7n9aCC/j3yDsON3M+iy9RxNn8RGzxQul+M4XruMA3sDGd0pe+6bztTz7inn0qMPCAhoef2Tn/yExx9/HHD21NevX9+yzGw2M336dAIDA7skTbG6M1bpk8oyzJTZzcwMnMNO/UFKqWXq+J3UWkr5u/1V0nJ1TPf34h8J0QS56Tt9/3MCfUny8eSefUd5PFcyzvoiDxh/TtikXWzcMJI9uixGBcxiU/l+Kg/l4jcw6uwbVbpFR9IUd1RBQQFhYWEALF++vCUR2Zw5c3jiiSdaLsCuXr2aZ5555qQ0xQsWLFBpihXlTHZ8uI44dwMNHm6k6wsYl7SdBnsev6r5O2l6f34aHcwHI+K6JMg3C3DT8XFif24LM7HdPZxfVb9Eg/0oY0fvYJcuB2n0IUxnJeWDDV3WBuXcdWaa4pdffpmhQ4cyYsQIXn755ZanTZlMJp588knGjBnDmDFjeOqpp1ouzKo0xYrSAXarg9cX/oYpYZNYbtzFoPjv0YdV8JLba+y3GLnF05sXxw3otvZIKXngh8N8YalntHsVDzUuoiY3kuNHZzCvNoFNxdtZ9Paf0HTgAvClQKUpPrtzTVOsPllKn7PnmxQGegWx3f0owZHp6MPKed7tn2RaPdHvqeDHYQFn30gnEkJwW7A/un0V7Gry5WX3l/GKysc3LJ09HoXEefixb+Wubm2TcmlRgV7pczJWrEbjH0FdcBoh/Q7ygmMJ2TZf7nDzQlvUQLivR7e3KczXA11ePQuN3qRbgnje8TSRA/ZSGrAbd1Mse5Z/2+1tUi4dKtArfUrOgeNEe4ST5pNC9OBUnnc8QZYmln8PjcWvxoZOIwjybv85sF0l3M/g/F4v+dfQWA5qBvCyfJT+Q35gp/duotxDyD984fOl+4reOKTcW5zPuVGBXulT9r79Ndn+1cQO3cA/tP/HYe1AlvoWMzfIl/zKBkJ9DWg13T/j1+imw9+op6CqgauC/XjBw8webSL/0j9Mv4T1ZPvVs+vNL7u9Xb2RwWCgrKxMBfs2SCkpKyvDYDCc03pqeqXSZ5TnFaEXRtyGfMoHnreyXwzjxYN/4aobfgdAfmVjjwzbNAvz9SC/shGABVHhVH/7D54a8FM8fWq5Mn412vSbqCouxTf49KdaXUoiIyMxm82onFdtMxgMRJ7j3dQq0Ct9Rtq/V5E/YDU7giawWczgl7Z93FL6PQS8D0B+VQNJMf491r5wPw/MFfXON8FDuC/vU0oHzOVlcTlBEcUMr97Eztd0zHz6xz3Wxt5Ar9fTr1+/nm5Gn6KGbpQ+oaaghEKPQ2T19+FzcROxMpNH8z6CkMtAq8PukBRWNRLu13M9+gg/A/mVrrwtbp4QGM+vCz8jynGQT8QtFAyS5LtlUVtU1mNtVPomFeiVPmHzW59xfPhx3hb34Vu9h/uDGhEFeyB8JAAlNU3YHJKwHgz0YX4eVDfaqGm0OgvCRyLyd3NPUB2+NQd5Q/MAJSMOsPHNT3usjUrfpAK9ctEr3n+MzEGbWer+ILEVBeiq/slEz3Boqm4J9PlVzp50hN+5XcTqTM3/TRRUOcfpCUuEmgLG+/RDV/l3AmsqedXwMEcHrKUoI6fH2qn0PSrQKxe9FT+8yJum2zFYrdy+92MC3DyIr610LgxPBGgZMunpoRuAvObhG9cfoaENDXhqHTyw/b/Y7W68GXg7Kzf9raeaqfRBKtArF7Vtqz7mnf5jqcDEb/+1jDUhR5gUMQlNQRroDBA0GIAC12yXsB6eddO6LYQOAwT6wr2MCxtHSvARfv3GuxQRyrL+SWxf9b8ea6vSt6hAr1y06mrK+bvI46AYyg3bMxkpizjsWc2UiCmQn+YMpFpn0rLcinq83XX4GHpuolmIjwG9VpBT7pp54+4FQYOgII2pkVPZ5lfKpPJDXJeayQHNcF7iOLXVpT3WXqXvUIFeuWj9buO/WK+fwYzcvVyzcg1HRoeiERomhU2Agj3OMXCXY2X1xAZ6dspDHM6XViOI8jeSU153ojAsEfJ3O/84CYF5Yj+u/XYt03L3stZtFn/Y9K8ea6/Sd6hAr1yU3tv+Hz70SGZoXRZTf9ARXLKLr2LKSAxKxLemCCw1LWPgAMfL6ogJMPZgi52iA4wcL6s/URAxCmqLCLXZGOg/kG/j6wgtSmXqZkioy+Z9jzm8t/39nmuw0ieoQK9cdPYVHeT3dREEOkq55qt0YiuO4D5yKNsdWUyJnALmHc6KUeMAsNod5FU09IpAHxvgyfGy+hO390eOcX7P3cGUiCmsc2TgMWwQMfUlXPvVLgIdZfy+Lpp9RRk912jloqcCvXJRqbU2cvf+I9jQ8sCWPWgcQwjL+Ibcyf0BnEMgudvBwwQBzrL8ygZsDklMgGdPNh2AmAAjtU02yuoszoLQYaDzgNwdTI2cil3aKZw8iND9yxG2gTzwQzo2tNy9P5tqS/2ZN64o7VCBXrloOKTkjm2rMRPK3dnbcGSXEmV0YNBZ+Ta2ihBjCAP9B0JuirOn7BqPP+YaKontJYEeODF8o9U7h2/MOxgeNBwfNx/WDmjE01pOhK8X8oiZe45sx0wod2xfg0Ml+lLOgwr0ykXj9/u28IMtmvlVq0nYXoFdn0TYro/xmjObTRUpTImcgmishNJMiBrTst7xMufFz9heMHTT/F9Fc5sA5x+lgr3o7DYmhU9iTfUOPKdMJuLAF9j0SQxNqWJ+5Xdss0XzVPrGHmq5cjFTgV65KHyRf4ylpV6Mt/7Aj7YFUyitBAX4412SScG0wdRZ65gaMRXMrkdQRo5tWfd4WT0eem2P5KE/VaS/Bxpx4r8MAKLGgsMKBWlMiZxCeWM51TNG4X0shcCgcPLtNfxoRygTLZt5o8yXj3MP99wBKBclFeiVXu9ATR0/zywmTh5mwc468gp3IsVYYos3oY+OZp2pEL1Gz7iwcZC7A4QGIka3rN8846Ynp1Y2c9dpCffzIOekHr3rj1LudiZHTEYg2BTXiC4ggNiqHTjEOAoLd3F9agPxjoP8X1YFaVWVPdJ+5eLUoUAvhJgrhMgUQmQJIRa3sVwIIV52Ld8rhBjVatlbQohiIcS+zmy4cmkos9i4bfdeDLKWm7M2EZ7vQ51HIKbgYHy2forvtdeyNncd48LGYdQb4fgPzguc7l4t2zhWVt8rZtw0iwkwntyj9woCU384vhV/gz8jg0eyJn89ftdfh8/GD/E2hVFr9CSuNJDrMrfiLav4cVoGxU3WnjsI5aJy1kAvhNACrwDzgATgFiFEwinV5gHxrq/7gNdaLXsHmNsZjVUuLVaHZOGevZTatNxZ/j4Tjs1kf+VG7I7RxLsfQwgonjaEvNo8ZsfMBmuDc2pl7JSWbTgckpzy+l5xIbZZTIDnyWP0AP2mOP9IOewkxyRzqOIQtXMnoLFZ6W/Mx2JNIr10LWPzpvLj4o+otgtu351Go93RMwehXFQ60qMfC2RJKbOllBbgI+CaU+pcA7wrnbYBfkKIMAAp5UagvDMbrVwafpuZxfZaDbdblhGfPoaCkl3oTUPxDgjEd81beE6cyHeNaWiFlhlRM8CcAnbLSYE+t6Iei81Bv8DeE+jjAj2pqLdSVtt0ojB2CjRVQeFeZsXMAmCt4wCeEycSuOkdPLyD0AUMoKbsINH7R3Fnw9vsbdDzywOZ6pF7yll1JNBHALmt3ptdZeda54yEEPcJIVKFEKnqEWLKO+YilhXWMdfxNdF7vRhsCSfXlkVj40gGh9fgKC7E//bb+e74d4wNHYufwQ+ObnKOz8dMaNlOZmENAINCvXvoSE7X3JbMopoThbGTnd+PbiLUM5ThQcNZfWw1fjffDHnHGdLfTkP9KLLq95JojyV4XyDX2f/HZ6VNvJpT2ANHoVxMOhLo27qCdWoXoiN1zkhK+bqUMklKmRQUFHQuqyp9zJqyap44nE+iTCXp4BHGVUxlh3k5vmHT8fDyImDz++hjojEPCyanJofk2GTnisc2OXPHGHxbttUc6AeG9MJAX9gq0HuHQuBA5zEAs2Nmk1GeQcWYAWgDAwlL/xyDly++EZNJzfuKiZUTGJqZx1j5A3/MLmRNWXVPHIpykehIoDcDUa3eRwL551FHUc4qvaaen6RnES2PckPeN/TLn0m9xYzDT09VWRxDh+mx7tmJ6bbb+froN+g0OpKjk8FS75xa2dwzdsksqiHK5IGne+95PHKQlzv+Rv3JgR6cbT++Few25sTOQSBYkbMKv+uvp2nj91w2xo/q8kE0eDRhtZcRXTCDawq/JVoeZdG+Ixyqa+yZA1J6vY4E+hQgXgjRTwjhBiwAlp9SZzlwh2v2zXigSkpZ0MltVfq4/EYLt+85hFGWc3/NG8hj0xiijyQ1byUGvxl4+hoI2/8FGk9PvK6dz4qjK5gSMcU5bHNss3Muety0k7Z5qKiGQSE+PXNA7RBCMCjU++ShG4C46c5kbOYUQj1DGRs6lq+zv8bv5h+BEEQcW4O70YBfeDLbc5czXBuN9eh0HmhYitZRzY/3HKLcauuRY1J6t7MGeimlDXgYWAVkAB9LKfcLIRYJIRa5qq0AsoEs4N/Ag83rCyH+A2wFBgkhzEKIezr5GJQ+oMZm57Y9WVRbG/mF9UXKMkYxrW4oW81fENQ/gaqSIEZODaB+5Qp8r7+e1Jr9lDaUcnX/q50bOLwK9EaIOdGjt9gcZJfUMSjUq5299pxBId4cKqzB4Wg1whk3HTQ657EAV8ZdSU5NDhlupfjMmUP9p/9h2ORgSvJC8I/pR0rxKqbWDqN0fyKP2F4kv6mJ2/ccps5u75mDUnqtDs2jl1KukFIOlFL2l1L+yVW2VEq51PVaSikfci0fJqVMbbXuLVLKMCmlXkoZKaV8s2sORblYWRwO7t2XTWZ9I484/kbVvhgmVI+i1lBKaVMu6KbiE2gg9Mj3YLNhuu1Wvsr+Cm+9N1Mjp4KUcHg19JsG+hPPhM0urcXmkAwK7V09eoBBoT7UWewnHisIzmsL0RPg8HcAJMck465156sjX2G6+y4ctbXEVGzHw9sNd69ZmKsOYjPWklQ1Bsv+EB50/J20mgbu3XcUi0NNu1ROUHfGKj3KLiUPHzjOhoo67pGv4X7Ig/51E4jyCGLTwf8SP34+lUV6kpLDqf7PB3hNm0ZDmD+rj63mirgrcNe6Q+khqMyBgbNP2nbLjJtedCG2WfN/GaeN08cnQ9E+qMrDy82LWTGz+Cb7G+Tg/hiTkqh9723GzIumxOxG/Li5rD/wAf3dwoisG48py8FCuZR15bX8/GCuSoCmtFCBXukxUkoezzSzvKSKW+U7DMgtxa16AqMaotla9CX+UZGUFgzAFO5J4OG12CsrCVx0P8uPLMfisHDTwJucGzrkHOpgQPJJ2z9YWINeK3rVHPpmzbOADhaeMlsm3vXH6vBqAG4aeBO11lpWHVtFwKJF2IqKiCjejn+okYriBDyDTKRWrmJsbRyOygkkFOVxs3yfz4oqeCorT82xVwAV6JUeIqXkj9kFvF9QxjXyUyZVZlCSP4FpFYMoNRVjLs0kcugN1JRamHRtDBXvvI3nxAkYRozgk0OfMCxwGINMg5wby1wBwUPBL+qkfaTlVDIkzAc3Xe/7mHsb9MQFepKWW3XygqDB4BcNB78BYFTwKOJ84/jk0Cd4TpqIYfhwKt54nYnX9qO6zEa/0TeRnb+bmuAqZlQMIffoJKbXpDFPfs0b5lJePF7UA0en9Da97zdA6fOklDyTXcArOcXMYjVXN24kM2MCM+uG4h7mwbqUZQydNpfDqYL4pGA8U1ZgLy0l8IEHSC1KJbsq+0RvvioPcrbC0OtO2ofN7iAtt5JR0f49cIQdMzLan905FSf3uoWAhGshex3UlyOE4KaBN7G3dC8Hyg8Q+OADWPPy8M3cQHSCiew9BgZPmsF329/EK9iLmY3D2Z8+gRstq5nCRp47WshLx1Swv9SpQK90Kyklf8ou4OWcYi7ne263fsXunZOYph1FiMHEmoy3CYiMpqlpLBqdYPycUEpffx3PyZMxjhnDsv3L8Hf3Z16/ec4NHvjS+f2UQH+wsIYGq52R0X7de4DnYFSMH2V1lpOfIQvOY3HYWnr11wy4Bk+9J+/ufxevadPwGDWK0n/+kwlXRWFttKM1TMc7OIh1Rz8kXGtiitto0nZO4F77h0wWW3nmaAF/P6bunr2UqUCvdBspJX84UsA/c4pJFuu4S37CzpSJjPEaSUyNiQNsp7q2lOGz78F8sIZxV8fR+L/3cFRVEfzLX5Bdmc0G8wZuGXwLBp1rds3+z5zZKgMHnLSv3TkVAIyO6b09+ua27XK1tUX4SPCPdR4b4O3mzQ3xN7Dq2CoK6woJ/tUvsZWUwJrPGD0vhiO7qhg5715KSo9zxGMf/SsCGO45kt2pk1jkeJspYit/OVrI80cL1Zj9JUoFeqVb2BySX2Xm8mpuMbPFOu6UH5K6YyLxPpcxtCiEmgH17E1fw4Qbf8zetQ0ERHgxaICk/J138LnqKgwJCbx74F3cte7cPPhm50YrjjsTmQ29/rT97TxeQbC3OxF+Ht18pB0XH+yNl7vu9EAvhPOYsjdArTPv0+1Dbgfg/Yz3MY4ejdfMmZQtXcrw0UYCIrzYu97GmPk3k7JzOY2D7SQWhhPleRk7UyZzv3yTaWIrfz1WyNNZ+Wo2ziVIBXqlyzXYHdy7/ygfFpRzvfiKO+V/SN0+hWjfYYzPi0YkeLBqzavEJo6msnQgDTVWLr9zCCXPPQc6HcGP/h/mGjNfZn3JtQOuxWQwOTe8611nErNhN522z105zvH53vCwkfZoNYLEKD92Hq88feGIBSDtkPYBAGFeYVzR7wo+zvyY0oZSQh5/DGm1Uvb3F5h5x2Aaaqw0NY0gYnAC33z/DzSDPJmcF0eI91B2pU7hJ/INrtCs43VzCQ9n5Kh59pcYFeiVLlVmsXHLniOsKq3iLvE+N7GCHdunEOZ7GZNz+6GP9ebrzS/jHRDEoIk/Jiu1mDFXxWLI3knt998T+MAi9CEhLN2zFI3Q8JNhP3Fu2G6F3e85p1SeMtumsKqRnPL6Xj1s02xUjD+ZhdVUNZzyEJGgQRA9EXa+A66gvGjEIqwOK2+kv4FbTAymhQupXv4VnoUHGTk7msxtxQyf/RMMXt6sSHkV90gvpucNwGQcQtrOqdzu+JBbxCd8VlTBHXuPUqXSJVwyVKBXusyB2gbm7jzErupafsrLzGEb27dNJthnMNPM/XEP92btkXexWpqYufD/2PJZHmEDfBkxMYDCp3+He3w8AXfeSXZVNl9lf8XNg28mxDPEufFDK6G2CJLuPm2/6zKLAZg6sPdnQZ02MBCHhI2H2kjNnXQ3VByFY84Hgkf7RHPNgGv4OPNjCmoLCLz/PvRRURQ8+SRJl4cSFO3ND58WkHzvo9TXVLKx4H8Ygry4vHAQXu4DSdkxnWvEOhaJ19lcUc2Vuw5zpF4lQrsUqECvdIkVJZVctfMwjdZafut4nAkUsGXzZEL9h3B5wUAMwV6k1K2iKDebuQ/8H1u/rESn1zD7nsso+euz2EpKCPvzn0Gv57kdz2HUGbnnslZpkna8Dj4Rp90kBfB9RjERfh4MDOl9OW5OlRjlj79Rz9qDxacvHDIfPEyw/fWWokXDF6EVWp7f+Twao5GwP/4Ra04OZf94iTk/GYp0SHauqmf2op+Tm7WPPXIjHr6ezC5NIMAnni2bJzONDH7NEkqb6pi38xBrVYrjPk8FeqVTNTkcPHnYzMJ9x4jRlvC07QHipB+bNiUxIHwYM/PiMfh5stu+nkO7tjD9rvvITHWnsqie2fcOxbH1e6o+/YyAn/wEj2GXsT53PVvyt/Bg4oMEeAQ4d5K7A45uhPEPgvbk9MONVjtbskq5fEhwrx6fb6bVCGYMCmZdZjF2xykXSfUGGHsfZH4DRfsB51j9wmELWXVsFTsKduA5biz+t91Gxbvvodm/g1l3J1CcU0Nuhj+TF9zBvm3fs9+wAw+jB8lFQ4gOTmDzpgkMkYLf2R4miHJu25vNn4/kYz11/0qfoQK90mmO1Ddy1c7D/NtcylW6H3jM8hButSPYvGkIw6KGMzErGjeTkT2azRzYvo6pt91NbUU8x9PLmPKjeEI8aih48ik8EhMJevghai21PLPjGfr79mfB4AUndrTxb86ebhvDNluzy2iw2pk5OLgbj/zCzBwSTGW9tWVK6EnG3Q9uXrDp+Zaiu4feTYRXBH/Y9gcabY0EP/Yo7oMHU7D410QG25hwbX8OpxbjYDTjb1hA2qYVHPTahbu3B9OO92NQ+GVs3DgM74ahPGF9kNn6NF7OKeba3Yc53tB0ehuUi54K9MoFs0vJv3KLmZVyiNyGWh4VL3K7498UHr+C3btjmDFgAkmZoRhifdjNBvZtW8OUW+/CwUjSN+SRmBxNwigfch94EKHXE/H83xB6Pc+lPEdRfRFLJi1Br9E7d2ZOdabxnfAQuJ2ew2b1/kI89FrGxwV081k4f1MHBqHTCFbua+OmJqMJxtwL+z6DogMAGHQGnp7wNMeqj/HSrpfQuLsT8cILSJuN3AcfYsTkQIZOCWfXquMYfCcz7rofsWvj1+w37sA93IcJWZFMiR9PSkoUFfmzuMv2N34uXuNQbS2Xp2TyTl6pmoLZx6hAr1yQ/bUNXLHzEE9n5TNce5Q/2BYxWhSRmjKXkpJQrgufRv99nrgn+LHO/BEHtq5l8i134GAkqSuOMWRiGOOvjMT8s59hyc0l8h8vo4+IYPWx1Xye9TkLL1vIiKARzp05HLDi/8Ar1NnTPUVtk43laflcMSwMg17bzWfi/PkY9MwaEsKnu8w02drIJT/pZ+DhB98+5kzJDEwIn8CCQQt4P+N9Nudtxj2uHxEvPE9TZib5jz7GlBvjGJAUzLbPs/Hwm8q4635E2voVbCn/HPcBvgxK92Z+1DTM5nDSdl/NBHGMP9gfZKA2n8WHzFy7O0s9saoPUYFeOS8lFiuPZeYyOyWT3PoafiaW8oh1Mbrq0axfN5YAn3iuc4zHdFSLbqIvX+54iYLsg8x7+Fc0Noxg57fHSZgUxrSb+5P/y19Rv3UbYX/4PcYxYzhUcYjfbvktw4OG8+CIB0/sdPd7kL8bZv8R3E9PPfz57jzqLHZuHx/djWeic9w+PoaKeivfprfTq5/5pPN5sq67ZQF+mfRLBvkP4rGNj5FTnYPX1KmE/PY31K5dS+ETTzDrjsHEjwlh2xfZaPQTSb7vpxw/kMbX+5aiH+9PcJaO67UT8Db0Y/36cRjrhvLLpkd4UPM2B2uquTwlkycPm9VTq/oA0RtviU5KSpKpqalnr6h0uxqbnTfMJfwzp5gmu505uu3Mt76GjwgjLW0EdXXeTIkcTdxhL3RebpQPqOD7r97AzcODeQ8vZu96Gzn7y0hMjmb83HDyf/lLajdsIOSpJzHdeiuFdYXc+e2dWBwW/nvVfwk2usbay4/Cv6Y60x3c9Y3z7tFWpJTMe2kTGiH45pHJF8WF2NYcDsnM59cT5O3O/xZNbKOCHf49A6rMsGgz+IQDYK4xs+CbBfi7+/PO3HcI8Aig9N//puT5F/CeO5fQv/yFLZ8dY9/GPOJGBjF4rJ0VrzyHtNuZdfX9+BwwYmu0cji+mh9yduHrW8Ow4XupcBTxpdvP+M46Ai+dlp/FhHJXRACe2ovnP6VLjRBip5Qyqc1lKtArHVFqsfGGuYS38kqotjmYqM/iesvfidTaOHZsNMeOBtIvNIZxlXF4VWrRD/Xhh2NfkL0vhehhiYy6YiFbPimkrqqJqQsGMnCgDvMjj9CYvo/Qp5/G/+YfUVJfwsJVCylpKOHN2W8yNHCoc+c2C7w1B8qPwP2bwD/mtPZ9n1HEPctSeeb6Ydwy9uLr0QO8sSmbP36TwUf3jW/7GkPJIXh9GoSPgju+bJlxtKtoF/d/dz8xPjG8MfsN/Ax+lL39DsXPPotx7FjCX3ie/Wn1/PDZEXwCDUy4NoiUL18n7+ABBiZOYmzYFVgPVlNtsrHN+wg5RbkMiC8iIiKFHIc3n7j9nBRrLCa9lnsiglgYGYi/vvc8bF1xUoFeOS9SSnZW1/NefhlfFlfQ5HAwUZfJXOtb9BclFBeP4FBmJP7ewYzXDCKswIg2wJ08Uw4/bPgIgMm33E1dbTx715rxNhmYc+9leJr3kr94MY6GBiKeexbvWbPIrszmgTUPUNFUwWuzXmN0yGhnIxx2+PRe55DFj96DhPmntdNiczD7xQ1oNYKVP5+KXntxjkg2WOzMemEDPh56vv7pZLSaNv4rSfsQvngARt4O8//Z8p/ND3k/8NO1PyXcK5xXZ71KlHcUVV99TcFvf4vW35+Ivz5Hpd8AVr95gLqqJi6bFo6H4QBbP3kfjU7HlGm3EVIUhq2ykfyIerZaMqiuK2Xw4AICg3aSKSP4Vn8PO2z9MWo0XB/iz+3hAST6GLv5LCntUYFeOSdH65v4qqSSz4oqOFjXiIewMVlsZbb9f4TTQH7eZRw9GoG/ZxAjNHHEFvmg83KnMqicDSkfUl9TSfzYSYQOuoL9m2qor7IwdEo4Y2cGUPnKS1R98ilu/fsT+dLfcR8wgBXZK/j9tt/jrnXn1ctfPdGTt1vhq59D2vuQ/AeY9Eib7X1lXRZ/XZXJO3ePYfqgi2daZVu+3pvPwx/uZsn8odw5MbbtSmv/BBufc86xn/sX0DiHU3YV7eKRdY8gpeT3k37P5dGX03jgAOaf/wJrTg7+t92Gz08eJHV9KekbzHj6ujN0shHz/q/I3rUDb79Apo6+Fe9Cb2z1Fo6GVZJmzaamoYT+/XMJCc0gB2++0y7gB5lEo9QyzMuD60L8uSrIl2gP9+47UcppVKBXzsjqkOysrmNDeQ3flVWxr9Y522KQJofJ9m+YwBZEYyTZ2RGUl0UQ4RVGQl04UfX+CC8dhe45bEv/nKamOiITRhAyYA5H92ior7YQ0s+HCclBuG3+kvJ338PR0EDA3XcR+NOfUmgt47mU5/g+53sSgxJ5duqzhHs5x56pLXb25I9ugGmLYcav22z7uoPF3LMshbmXhfLqbaO765R1GSkld76dwtYjpbx3z7i2h3CkhNW/ha3/hPg5cN1S5wVb4Hj1cR7d8CgZ5RnM6zePX43+FUHCm+IX/07FBx+g9fbGdNedWCZdzeavCyjJqcHL353Yy6zkH1xB/qEDeBi9GT/0eoLrw3HU28nxLGe/MZ/C2kICA3PpF5eL3b2YLUxms/YqshzOn1mitweXB/gww+RDorcRXVv/kShdRgV65STFTVZ2Vdezq7qOXdW17Kquo94h0OAgniMkyS2MkdvwanQnPz+U0tJoDDKYAfYwBtQH4YORWq9aMku2c6RgJ1o3PeEDxyLchlOS44HDLokc7MfQfo14/PAVNSu+QdpseM2cSfAvf0GuSfLBwQ/4IusL9Bo99w2/j7uG3oVOo3P24tM+gO+eBms9XP0SJN7a5nGsyyzm4Q92ERvoyf8WTcDo1jfGjasarFz/6hZKay0svX00E/q3c09Ayhvw7ePg4e+ciTTsJtBosdqtvJ7+Om+lv4UQguvjr+fWwbcSWthE8fPPU7dxE8LdHZ+r51M35kr2ZenIP1yFRicIia7D2riH/MwUpM1BfORYBpnGYKwxUinqyPIs4YgooEmUEBR8jIjwfKoMku2MZ6dmCkdkDA40eGthlI8no3y8GOljZKSPkSA3ffeeyEvMBQd6IcRc4CVAC7whpfzLKcuFa/kVQD1wl5RyV0fWbYsK9BdGSkmVzU6hxUpeo5Uj9Y0cqq0hq66GIw12SmzOf/W12ImSxxlAJpexl3jLUZoqvKmsCKWqIhI/eyhRTSaiHYH44kGlKCG7NA1zzSFswopPcDwabX/qa6NBGDB6aYk21RFWvgvdpq+xV1UhPDzwue4aiq4cwxZNNpvzNrO/bD96jZ7r46/n3mH3EmoMgfJs2P857FoGlTkQNR7mv+zM4niKstomXlt/hDe3HGVQiDfv3D2WUF9Dd5/mLpVbXs+db+/geFk9i6bFcd+U/vga2wiUhemw/KfOaaem/jD6Tki4BvxiMNfm8e/0f7M8azk2aWN40HAmR0xmYlM0QV9to3r5V8imJrQmE9bJV1PgN4KcUg8a6h0IGvDwOo7dkkV1yRHccCfKZzBxASPwlQFUUM9xTQlmQwWVmkL8THmY/PPQ+TdwSN+ffQwji0GYRTQO1yzuEJ2d/h46Bnh5Eu/pwwCjgUiDG2Huerx0ajbPhbqgQC+E0AKHgGTADKQAt0gpD7SqcwXwU5yBfhzwkpRyXEfWbculGOgdUmJxSCyu71bpwGJ30Gi3YHFYabJbqbFZqbJYqLJaqLI0UmVtotpqpcZmo9LmoMSmpdSup1y6Y+Xk3q1R1hJGPuHkEcMxoi05mGqqsdT6UFtrgupwfBpDCHD4ECL98LRCbVMZpY1mShrNVFiK0boFIkU4QhuORh+NRujxF+X41hzB73gKfiUHEUjs3h6UjuzHwQQvNkfUktmUQ6O9EQ0ahvkPJNk0jKuNMZjqypwBKmcb1LmSesVMgomPwMA5IASNVjtldRYKqxo5WFjN5sOlrMsspsnm4Eejo/jd/KF4uPXNIFHTaOXJL/bxRVo+nm5aZgwOZvKAQAaFehPiY8Dk6ea8MczhgINfww//APMO58reYRA1DsITKfH0Z3ndMdaU7WV/ZRYSiYfOgyFuMUzKNTL4QDWm3cfR1jXgEBoqg4dSFZ1EpVccldIfh6MJhy0HHHngKMBuKcXkFkqQIYogQySehkBqdDaKRCXlmhpqPYrQ+OZg9KxA61lPqbcfuW5R5BBDAeEUEEG9OPmuZg8sBGgaCdTY8NEJ/HQa/HRa/N30zi93A546A0adG546dwxaPUatDg+tDqNOh5sQaIVAK0AnBJqLbHptZ7jQQD8B+J2Uco7r/a8BpJTPtKrzL2C9lPI/rveZwHQg9mzrtuV8A/3YNV9iFc5ej0QgOfHDbuu188hFqzIJrdaTzka3WffEWRMt60qE61Vby0/dhvO1Aw02dDjE+QUrvbTgQT2e1OFHBf6U4+uowstWi5e1Hm9LPf4NtRjrtdBgQjQE4l4fhKfNC3ebRGezQVMjVnsdtdYaqm31VNtt2IUBofFGow1AJ70wWuwYmqoxNpTgWZeHV10+7g0FlPlYKfIT5JsgO0yQFSYoNIEUgiCbg34WG/FWC6MbGxjb2IjvKYmz8ghiN4NIYzCbSSSPExdT7Q5JveXkO0WDvd1JTgjh7kmxDAg+/aapvuhAfjXvbj3GmowiSmstJy3zdNOi12nQaTToNIIYCphAGsMdBxlqzyBElp5Uv1KjYZuHgd3uBg67uZGt11Gm0yAckvByGJAviSuUhJVDSKXEVONGk0cYtV4R1HqG02gIoMHNmwa9wKptwGEvRzpqcMOCj1aPp9YNo84bg84bh8EDq05Dk1bSpK+jyViEzaMMDFU0GG1UeHpSo/eiRudFpdaHKq0vlfhThyf1eFGLFzZx/sM9WmlD4ECLA43rS4sdDQ7nb6hs/RvJSb+1J793xgXR8lvbXr0T708t6ygvez2bZt94XuueKdB3ZFAzAsht9d6Ms9d+tjoRHVy3uZH3AfcBREef3zzosKYyrELb5g/O+YM9ubx16AeBkKf/IE8Ky/LkdU5e3/lNnPKhOH3fze+dfxI0ONBJOzppR+t6rXE40Doc6BwSrUM6X9slBrsNg82GwWrDYAGDBbRNGqRVj7TqwKpHWHUIhw6Nwxfh8AGHxGG3YXVYsWHF6qimwVFOJQ6EkGikRCMkAjsa6UAjbPhqG5CaOhzaBhy6BuxaO/U+UG7Q0uShocGgweGhRbiFoNfocJca3NEwDA0TpR7fOj2+Uo8BLQ6hwaI3YnH3YIvJE4vWg1pdAFVuIVS5BWPRnujZjefk+6C0QuDv6UaApxtB3u4MDPEm0t/jorsZ6kIlhPvwlxuG43BIzBUNZBbVUFrbRFltE+V1VmwOBzaHxG6X2ByBHHMM5bBD8ilgsNfhbyvGZC3C216Bu6MBg6OeUbKBCY0NaBrsWLFRLSxU6y3U9LNhibNzDDuZwo7FYUPbVI2+sRK3hnTcmyT6Rgcam8TNoUUjfNDaDQiHAbvNQJXNQGVTKQ40yMpWnR6pQa9xQ6/zQ0cgRqEhSKMFrRap1SI1AoemBuFWhnBvQupsSJ0Nq5uk3l1LvZsWm05g04FVp8Gq0WLRarBqdFg0OuxCIIXALgQOIZAa52uJwCE02NEgXd8daJAndeBOkIiWX/JTwzfINtZrrxN3coToKA+75eyVzkNHAn1brT31/LRXpyPrOgulfB14HZw9+g606zRfXrnwfFZTlIuCRiOIDjASHaDmrivnpiOB3gy0flZbJJDfwTpuHVhXURRF6UIduYUwBYgXQvQTQrgBC4Dlp9RZDtwhnMYDVVLKgg6uqyiKonShs/bopZQ2IcTDwCqcUyTfklLuF0Isci1fCqzAOeMmC+f0yrvPtG6XHImiKIrSJnXDlKIoSh9wplk3F2f2J0VRFKXDVKBXFEXp41SgVxRF6eNUoFcURenjeuXFWCFECXD8PFcPBErPWqv7qXadG9Wuc6PadW76YrtipJRBbS3olYH+QgghUtu78tyTVLvOjWrXuVHtOjeXWrvU0I2iKEofpwK9oihKH9cXA/3rPd2Adqh2nRvVrnOj2nVuLql29bkxekVRFOVkfbFHryiKorSiAr2iKEof12cCvRBirhAiUwiRJYRY3IPtiBJCrBNCZAgh9gshfuYq/50QIk8Ikeb6uqIH2nZMCJHu2n+qq8wkhPhOCHHY9d2/m9s0qNU5SRNCVAshft5T50sI8ZYQolgIsa9VWbvnSAjxa9dnLlMIMaeb2/VXIcRBIcReIcTnQgg/V3msEKKh1blb2s3tavdn18Pn67+t2nRMCJHmKu+W83WG2ND1ny8p5UX/hTMF8hEgDufDTvYACT3UljBglOu1N86HoycAvwP+r4fP0zEg8JSy54DFrteLgWd7+OdYCMT01PkCpgKjgH1nO0eun+sewB3o5/oMaruxXbMBnev1s63aFdu6Xg+crzZ/dj19vk5Z/jzwVHeerzPEhi7/fPWVHv1YIEtKmS2ltAAfAdf0REOklAVSyl2u1zVABs5n5/ZW1wDLXK+XAdf2XFO4HDgipTzfu6IvmJRyI1B+SnF75+ga4CMpZZOU8ijO5zGM7a52SSlXSyltrrfbcD7BrVu1c77a06Pnq5lwPnT4R8B/umLfZ2hTe7Ghyz9ffSXQt/dw8h4lhIgFRgLbXUUPu/7Nfqu7h0hcJLBaCLFTOB/GDhAinU8Dw/U9uAfa1WwBJ//y9fT5atbeOepNn7uFwLet3vcTQuwWQmwQQkzpgfa09bPrLedrClAkpTzcqqxbz9cpsaHLP199JdB3+CHk3UUI4QV8CvxcSlkNvAb0BxKBApz/Ona3SVLKUcA84CEhxNQeaEObhPNRk/OB/7mKesP5Opte8bkTQvwGsAEfuIoKgGgp5Ujgl8CHQgifbmxSez+7XnG+gFs4uUPRreerjdjQbtU2ys7rfPWVQN+RB5h3GyGEHucP8gMp5WcAUsoiKaVdSukA/k0X/ct6JlLKfNf3YuBzVxuKhBBhrnaHAcXd3S6XecAuKWWRq409fr5aae8c9fjnTghxJ3AVcJt0Dey6/tUvc73eiXNsd2B3tekMP7vecL50wPXAf5vLuvN8tRUb6IbPV18J9L3mIeSu8b83gQwp5QutysNaVbsO2Hfqul3cLk8hhHfza5wX8vbhPE93uqrdCXzZne1q5aReVk+fr1O0d46WAwuEEO5CiH5APLCjuxolhJgLPA7Ml1LWtyoPEkJoXa/jXO3K7sZ2tfez69Hz5TILOCilNDcXdNf5ai820B2fr66+0txdXzgfTn4I51/j3/RgOybj/PdqL5Dm+roCeA9Id5UvB8K6uV1xOK/g7wH2N58jIAD4Hjjs+m7qgXNmBMoA31ZlPXK+cP6xKQCsOHtU95zpHAG/cX3mMoF53dyuLJxjuM2fs6Wuuje4fsZ7gF3A1d3crnZ/dj15vlzl7wCLTqnbLefrDLGhyz9fKgWCoihKH9dXhm4URVGUdqhAryiK0sepQK8oitLHqUCvKIrSx6lAryiK0sepQK8oitLHqUCvKIrSx/0/9PsKGQ7SuzcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# simulate diffusion and plot, time intervals of 50\n",
    "\n",
    "epsilon = 0.4\n",
    "v = v0\n",
    "for t in range(1, 501):\n",
    "    if t % 50 == 0:\n",
    "        plt.plot(np.arange(n), v, label = f\"t = {t}\")\n",
    "    v = v + epsilon * (A @ v)\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c6f83a19",
   "metadata": {},
   "source": [
    "We observe the bell-shaped curve (Gaussian distribution) characteristic of 1d diffusion, just as we'd expect. Note that once we constructed the discretized approximation, we were able to perform the simulation in Python purely via linear algebra! "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7b039ac6",
   "metadata": {},
   "source": [
    "## Solving Linear Equations\n",
    "\n",
    "One of the most fundamental tasks in applied mathematics is solving linear systems of the form \n",
    "\n",
    "$$\\mathbf{A}\\mathbf{x} = \\mathbf{b}\\;,$$\n",
    "\n",
    "where $\\mathbf{A} \\in \\mathbb{R}^{n \\times m}$, $\\mathbf{x} \\in \\mathbb{R}^{m}$, and $\\mathbf{b} \\in \\mathbb{R}^{n}$. This equation represents a set of linear relationships between variables, a single one of which looks like this: \n",
    "\n",
    "$$\n",
    "a_{i1}x_1 + a_{i2}x_2 + \\cdots + a_{im}x_m = b_i\\;.\n",
    "$$\n",
    "\n",
    "Collectively, the equations in a linear system define a \"flat space\" called an *affine subspace* of $\\mathbb{R}^m$.  \n",
    "\n",
    "> 1. When $\\mathbf{A}$ is square and of full rank (determinant nonzero), this equation has a unique solution. \n",
    "> 2. When $\\mathbf{A}$ is not square or not of full rank, then this equation may have either 0 or infinitely many solutions. \n",
    "\n",
    "In Case 1 (\"the good case\"), we can use a simple built-in `numpy` method: `np.linalg.solve`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "5b2a7a41",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.73684211,  0.21052632, -1.61403509,  1.77192982, -0.66666667])"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# solve A@x = b for x\n",
    "\n",
    "n = 5\n",
    "\n",
    "A = np.random.randint(1, 5, size = (n, n))\n",
    "b = np.random.randint(1, 5, size = n)\n",
    "\n",
    "x = np.linalg.solve(A, b)\n",
    "x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "079caf0f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.allclose(A @ x, b)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "31ae98b2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# manual approach (not as efficient)\n",
    "# compute the inverse explicitly and \n",
    "# premultiply by it\n",
    "\n",
    "x_ = np.linalg.inv(A) @ b\n",
    "np.allclose(x_, x)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "907f803d",
   "metadata": {},
   "source": [
    "In Case 2 (\"the bad case\"), in which the matrix is either not of full rank or not square, we need to resort to subtler means. Suppose that the matrix $\\mathbf{A}$ has more rows than columns: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "3e2fc1f1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[2, 3, 2, 3],\n",
       "       [2, 2, 2, 4],\n",
       "       [1, 4, 4, 2],\n",
       "       [3, 1, 1, 2],\n",
       "       [1, 2, 3, 2]])"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A_ = np.random.randint(1, 5, size = (n, n-1))\n",
    "A_"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6b78edbe",
   "metadata": {},
   "source": [
    "In this case, there usually are no exact solutions to the equation $\\mathbf{A}\\mathbf{x} = \\mathbf{b}$. If we try the method from before, `numpy` will complain at us: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "efa6e65e",
   "metadata": {},
   "outputs": [
    {
     "ename": "LinAlgError",
     "evalue": "Last 2 dimensions of the array must be square",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mLinAlgError\u001b[0m                               Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-22-9379ccf0981e>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mx\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlinalg\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msolve\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mA_\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mb\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m<__array_function__ internals>\u001b[0m in \u001b[0;36msolve\u001b[0;34m(*args, **kwargs)\u001b[0m\n",
      "\u001b[0;32m~/opt/anaconda3/lib/python3.8/site-packages/numpy/linalg/linalg.py\u001b[0m in \u001b[0;36msolve\u001b[0;34m(a, b)\u001b[0m\n\u001b[1;32m    379\u001b[0m     \u001b[0ma\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0m_\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0m_makearray\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    380\u001b[0m     \u001b[0m_assert_stacked_2d\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 381\u001b[0;31m     \u001b[0m_assert_stacked_square\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    382\u001b[0m     \u001b[0mb\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mwrap\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0m_makearray\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mb\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    383\u001b[0m     \u001b[0mt\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mresult_t\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0m_commonType\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mb\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/opt/anaconda3/lib/python3.8/site-packages/numpy/linalg/linalg.py\u001b[0m in \u001b[0;36m_assert_stacked_square\u001b[0;34m(*arrays)\u001b[0m\n\u001b[1;32m    202\u001b[0m         \u001b[0mm\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mn\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0ma\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    203\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mm\u001b[0m \u001b[0;34m!=\u001b[0m \u001b[0mn\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 204\u001b[0;31m             \u001b[0;32mraise\u001b[0m \u001b[0mLinAlgError\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'Last 2 dimensions of the array must be square'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    205\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    206\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0m_assert_finite\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0marrays\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mLinAlgError\u001b[0m: Last 2 dimensions of the array must be square"
     ]
    }
   ],
   "source": [
    "x = np.linalg.solve(A_, b)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8389002b",
   "metadata": {},
   "source": [
    "One of the most common ways to approach this kind of problem is to compute the *least-squares approximation*, which is the minimizer $\\mathbf{x}$ of the function \n",
    "\n",
    "$$f(\\mathbf{x}) = \\lVert \\mathbf{A}\\mathbf{x} - \\mathbf{b} \\rVert^2\\; = \\sum_i \\left(b_i - \\sum_j a_{ij} x_{j}\\right)^2.$$\n",
    "\n",
    "Note that, if $\\mathbf{b} \\in \\text{range}(\\mathbf{A})$; that is, if $\\mathbf{b}$ is one of those lucky values such that $\\mathbf{A}\\mathbf{x} = \\mathbf{b}$ does indeed have an exact solution, then we can choose $\\mathbf{x}$ such that $f(\\mathbf{x}) = 0$ by finding the exact solution. \n",
    "\n",
    "Otherwise, we need to satisfy ourself with an approximation, i.e. a value $\\mathbf{x}$ such that $f(\\mathbf{x}) > 0$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "f1b76991",
   "metadata": {},
   "outputs": [],
   "source": [
    "# compute the solution x, error f(x), rank of A, and singular values of A\n",
    "x, fx, rank, s = np.linalg.lstsq(A_, b, rcond = None)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "4645d531",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.48002131, -1.35801811,  1.41928609,  0.2972829 ])"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x # approximate solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "65f0e02e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.61640916 2.27171018 1.31965903 2.09589771 2.61640916]\n",
      "[1 2 1 2 3]\n"
     ]
    }
   ],
   "source": [
    "print(A_ @ x) # should be close to the below\n",
    "print(b)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "131f5156",
   "metadata": {},
   "source": [
    "## Application: Linear Regression, Several Ways\n",
    "\n",
    "Actually, the problem of minimizing $f(\\mathbf{x}) = \\lVert \\mathbf{A}\\mathbf{x} - \\mathbf{b} \\rVert^2$ has a special name in statistics -- it's linear regression! Specifically, it's *orderinary least-squares multvariate linear regression*.  It's usually written like this: \n",
    "\n",
    "$$f(\\beta) = \\lVert \\mathbf{X}\\beta - \\mathbf{y} \\rVert^2\\;,$$\n",
    "\n",
    "where $\\mathbf{X}$ is the matrix of observations of the dependent variables, and $\\mathbf{y}$ is the vector of observations of the dependent variable. $\\beta$ is the vector of coefficients, and it's the thing that we want to estimate. We do this by finding an estimate $\\hat{\\beta}$ that makes $f(\\hat{\\beta})$ small. \n",
    "\n",
    "By the way, if you are familiar with the topic of *loss functions* in machine learning, this function $f$ is called the *square-error loss* for estimating $\\mathbf{y}$, and is probably the most important of all loss functions for regression tasks. \n",
    "\n",
    "Let's use least-squares approximation to perform 1d linear regression \"by hand\": "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "998e44dd",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x7fab3a8f2d90>"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD4CAYAAADvsV2wAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAav0lEQVR4nO3dfYwdV3kG8OfJxkEG2i6tF0g2WZwiYwhNEyfbxNQtStxSxwbJLolIKIU2QrLSJlWLqIVBFVT9EK4imhJCY1k0gkgoCSLp4haD22LagME0a2zng+DWTYDs2iLmww6QrbA3b/+4d9PruzNzZ+6cM3POzPOTrOzeO7kzkxu/c+Y977yHZgYREWm+s+o+ABERqYYCvohISyjgi4i0hAK+iEhLKOCLiLTE2XUfQJZly5bZ8uXL6z4MEZFo7N+//3tmNpb0XtABf/ny5Zienq77MEREokHy22nvKaUjItISCvgiIi2hgC8i0hIK+CIiLaGALyLSEkFX6YiItMnUgVncuvswjp6Yw3mjS7Fl3UpsWjXu7PMV8EVEAjB1YBbvfeARzJ2aBwDMnpjDex94BACcBX2ldEREAnDr7sPPB/sFc6fmcevuw872oYAvIhKAoyfmCr0+DAV8EZEAnDe6tNDrw1DAFxEJwJZ1K7F0ycgZry1dMoIt61Y624cmbUVEArAwMasqHRGRFti0atxpgO+nlI6ISEso4IuItIQCvohISyjgi4i0hCZtRUQK8N3vxicFfBGRnMr0uwnhQqGUjohITsP2u1m4UMyemIPh/y8UUwdmPR7tYhrhi4h0ZY3Cpw7MYrZAv5vezzqLxLzZGe8vXCiqHOUr4IuIIDtdA+CMn/v197vp/6z+YL/AZWO0PJTSERFBdrom6b0FSf1usrbvZQDWbNtTWWpHI3wREQzfnviDb754UVqmyMjdx0InaTTCFxFBdnvitPfGR5cmBum07UfIxNddL3SSRgFfRASd9sRLRs4MyEtGiC3rVhZuXZy2/YfecgmSQ341+XwFfBGRBf1zq93fN60axwfffDHGR5eC6Izsk1I5C7K2r2KhkzS0lNnjEExOTtr09HTdhyEiLbBm257Essvx0aXYu3Wts/30V/AAndF/1gWkCJL7zWwy6T1N2opI5UJ46rRfFWvKAtUsdJKmdMAneQGAuwG8HMBzAHaY2Yf7tiGADwPYAOBZAL9vZl8vu28RiU+Z9gQ+nTe6NHGEn1RjXzZY+17oJI2LHP5pAO82s9cAWA3gZpIX9W2zHsCK7p/NAO50sF8RidCw7Ql8yzMxG0qLhGGVDvhmdmxhtG5mPwLwOID+S9dGAHdbxz4AoyTPLbtvEYlPWopk9sQcLtz62UofROqVZ2I21ItVXk5z+CSXA1gF4Gt9b40DeKrn95nua8cSPmMzOncBmJiYcHl4IhKAtNQJgDNGzYCbFE+RFMygVEtVeX5fnJVlknwxgPsB/ImZPdP/dsK/klgeZGY7zGzSzCbHxsZcHZ6IBCIpddLP1ajZdQqmzpJKF5wEfJJL0An2nzSzBxI2mQFwQc/v5wM46mLfIhKX/tRJGhejZtcpmKIPYIXGRZUOAfwDgMfN7G9TNtsJ4BaS9wK4EsBJM1uUzhGRduhNnaTVv7sYNbtOwdRZUumCixz+GgBvB/AIyYPd194HYAIAzGw7gF3olGQeQacs80YH+xWRBtiybmXig0guRs15Sy2LqKuk0oXSAd/MvozkHH3vNgbg5rL7EpH6+HpYysWoOe3Y8lxMQnwIzBc9aSsiA/l+WKrMqDnPsWWtYhXiQ2C+qJeOiAxUVZ+ZvAYtH5j32EI+r2HvNtRLR0RKCan+3OXygSGfl4+7DbVHFpGBQqo/z7t8YJ5jC/28XD/Fq4AvIgPlrT+fOjCLNdv2eG2RkGf0nbfKJ6S6+iruNhTwRWSghYelXvLCJc+/9oKzzwwfVTUWy1o+MM/iJL2KLmziUxV3G8rhi7RY0UnC/z313PM/n5g7dUaOOSsl4TKAppVaDhuosyqEqizZ9Pk8wgIFfJGWKjpJOCigN20BkapLNqs4LwV8kZYqOiIfFNB9PNWapoqnXX3esaTdOfg+L+XwRVqq6Ih8UI45pAlQF3zdsdS5iIoCvkhLFZ0kHBTQQ5oAdaHsJGpaxVKdi6gopSPSUkUnCfPkmMu2SOj97KtfPYYvfvN4bT1uykyiZuX/63zYS60VRFoslMZh/QEySZlKnDLHNcx/n6yWDQC8tnNQawURSRRKq988T8/6KPEcZNj/Plmj+Nuuv9R7+WUa5fBFpHZ50xmxrB2blf+vc65DI3wRqV3Wwub928VgUP6/rjsrjfBFpHZ5FjaPqcRz06pxXHv5OEbYWRtqhMS1l9efPlPAF5HaJaU5fnf1RLQlnlMHZnH//tnnWzfPm+H+/bOV1NpnUUpHRGrVXwlz2/WXRhPY01TVV6goBXwRqU1TlxgMaWGVXgr4Ii1Wdx1+qCPhsvL0Farjv70CvkgEfASHvKNrn4Ep1JFwWYOqdOq6s9GkrUjgfDXbytPTxXejr5CWGHRpUK19Xf10NMIXqUGRUbOvtEee0bXvlEsVi37UJavWvq47G43wRSpWdNTsKzjkGV37DkxN67CZV113Ngr4IhUrejvvKzjk6V9fRWDatGoce7euxZPb3oi9W9c2PtgD9a0d4CTgk7yL5NMkH015/yqSJ0ke7P55v4v9isSo6KjZV3DIM7pu2qImoajrzsZVDv/jAO4AcHfGNl8yszc52p9ItIouBehzrdNBPV2qWj+2jerop+Mk4JvZgySXu/gskaYbZqKyzjbGobRQlvKqrNJ5HclDAI4C+FMzeyxpI5KbAWwGgImJiQoPT6QavkbNdT9EJeFztuJVd4T/z2b2Swnv/SyA58zsxyQ3APiwma0Y9Jla8Uokn6QVo+pYIUrql7XiVSVVOmb2jJn9uPvzLgBLSC6rYt8ibVCk8idtcW1pvkpSOiRfDuC7ZmYkr0DnQvP9KvYt0gZ5K3+a2qxM8nFVlnkPgK8CWElyhuQ7Sd5E8qbuJtcBeLSbw78dwA0W8urpIpHJWy9f1yP9EgZXVTpvHfD+HeiUbYqIB3krf5rarEzyUS8dEU+qrJrJW/lT9BkAaRYFfBEP6siV56mXL9OsbJgLmEpFw6KAL+JBqAt7DPsMwDAXsNgniJt4sVLAF/HAVa7cR9AZ5snZQRewpOMM9aKXR+wXqzTqlinigYsuk74XHyki6wKWdpxJcwVZnxWSplYzKeCLeOCiy2RIQSfrApZ2nCNkoc8KSdrFKu31WCjgi3jgov1tSCWUWRewtOOZN4u2tXLaxSrt9Vgohy/iSW+ufCHH/a77DubOxYdUQpk12Xvr7sOJxznek8sPbeJz0NzIfMpzoWmvx0IBX8SzYScAQ1vvNW2yN+s4Q2ytnOf7GE+52I5HkI7KopSOiGfD5uKHTQtV3RytjtWbypxjnu+jqSt9aYQv4lmZXHzREXJd5YRVjuTLnmOe76OpK30p4It4lpaLNwBrtu0ZKpCk5aBjrn3Pq+w55p0bCTEdVZZSOiKeJaUHFgxTW59Vnx9SZU8SF+mmsufY1HRNHgr4Ip715riTzJ2ax7s/dSh3EMwa4bp44MsXVw+SlT3HOuYcQqGUjognSWmXd913EEmFfQvlfnny0Vkj3NuuvzSoyp5ertJNLqqXmpiuyUMjfBEP0kazoy9cMvDfHVTBkzXCDXn06irdFPI5hk4jfBEP0kazLzj7LCxdMrLovX5ZQXDQCDfU0avLB8lCPcfQaYQv4kFawD45d+qM0WmefjP9E50AohzhtnmyNBQa4Yt4kDWa7W+5kDVaT6s5/+CbL8berWsrOBN3mlrbHhMFfBEP8k4sDgqCTaurVyqmXgr4Ih4UGc1mBcHQ6+olLgr4Ip64GM2G1DFT4qdJW5GaZT19qolOcUkjfJEaDWoEpolOcUkBX5zyseh2k+WZlHU50anvp90U8MWZulrzxizPpKyrIO3z+9GFJA5Ocvgk7yL5NMlHU94nydtJHiH5MMnLXOxXqpOny2FIi27HYlAjMFcNxwB/34/LYxS/XE3afhzANRnvrwewovtnM4A7He1XKpD3L3SIJYRVr/5U1KBJWZdB2tf3owt9PJwEfDN7EMAPMjbZCOBu69gHYJTkuS72Lf7l/QsdWmveGEaegxqBuQzSvr6fEC/0kqyqssxxAE/1/D7TfW0RkptJTpOcPn78eCUHJ9ny/oUOrYQwlpHnplXj2Lt1LZ7c9kbs3br2jNy3yyDt6/sJ7UIv6aoK+EkdopLagsPMdpjZpJlNjo2NeT4sySPvX+jQ2tY2YeTpMkj7+n5Cu9BLuqqqdGYAXNDz+/kAjla0bympyIITIfVKacJTqq7r8H18P3pWIB5VBfydAG4heS+AKwGcNLNjFe1bSor1L/SgC1UspYQhXUTTxHCM4ijgk7wHwFUAlpGcAfABAEsAwMy2A9gFYAOAIwCeBXCji/2KP0nBsEntePXMgLQRzRJT6UGYnJy06enpug+jddJ6tMewyEZea7btSUz3jI8uje7CJtKL5H4zm0x6T0/ayiJN68Heb+rAbGKwB/JN6MaSChLpp4AvizShuiXNwt1LmkETukoFSczUHlkWaXJdddLdy4I8pYSx1PaLJFHAl0WaXFeddZeSZ46iyXc/0nwK+LJIaA9QuZR2lzLeXVx82H+/CXc/0nzK4UuiQXXVsU5cFnmIzMe/L1InBXwpLOaJy7IPkcX6EJoIoDp8GYJq2EXCpTp8cUoTl27Fmh6T+CjgS2G+mpK1MfDFnB6T+KhKRwrzUbYZw2IleRVZZUt1/VIljfClMB8Tlz7bOVR551B0xK70mFRJAV+G4rodbtnAlxbUq06ZFL1wNaFnv8RDKZ0IhL4QtwtlHmjKSgdVnTIpeuFq8lPNEh4F/MA1KbedpUzgywrqVadMil64mvxUs4RHKZ3ANb1V8YIy8wJZQb3qlMkwT+JqtSipigJ+4NKC2eyJOazZtifKEsa0fPuwgS8rqFfdCkFP4krIFPADlxbMCDz/eky12z4mUbOCeh0BOOvC1cZnDSQcaq0QuKTlBgkg6VuLobWBr7YMMQTSNiwdKfVTa4WIJY1QyyzPVzdfk6gx5MHbMh8j4VLAj0B/MEsbJcdQux1L3bmPOwY9ZCV1U1lmhFzXbldZ5x9D3bmvUlgtniJ1U8CPkMva7arr/GOoO/f1sFYMFztpNqV0IuUqZ+0jrzwoHRJ6vt3nPAOgkk2pjwJ+5Mrmml0Ht6Syyy2fPoQ/3/kYTs6diiLI+ZxnCP1iJ82mlE7EXKRjXOeVk+4YTs0bTsydiqY1hFIv0lROAj7Ja0geJnmE5NaE968ieZLkwe6f97vYb+zKTpampWPe/alDuT/TdXDLc2cQer/3GOYZRIZROqVDcgTARwG8AcAMgIdI7jSzb/Rt+iUze1PZ/TWFiydO04LrfPdhujyf6TqvnPWcQK/QSxGVepEmcpHDvwLAETN7AgBI3gtgI4D+gC89XEyW5gmueT7TZXBLanOQRKWIItVzkdIZB/BUz+8z3df6vY7kIZKfI/natA8juZnkNMnp48ePOzi8MLmYLE1Kx5T9zLL60yEvOmfx8SkfLlIPFwGfCa/1t3r5OoBXmNklAD4CYCrtw8xsh5lNmtnk2NiYg8MLU9oI9ywydy6/P7iOMOmrqH40vWnVOPZuXYvbrr8Uz/X9n0AA116udIlIHVwE/BkAF/T8fj6Ao70bmNkzZvbj7s+7ACwhuczBvqOVNjqfNytUxbIQXJ/c9kZ86C2XeK8uKbtAtwH44jebe+cmEjIXAf8hACtIXkjyHAA3ANjZuwHJl5Od4SfJK7r7/b6DfTtVZYuBTavGce3lyaPcYatYfFeXFC0DVe8YkbCUnrQ1s9MkbwGwG8AIgLvM7DGSN3Xf3w7gOgB/QPI0gDkAN1hgfZmrXux66sAs7vnaU6nvDxsUfVaXaIFukbg5qcM3s11m9ioze6WZ/XX3te3dYA8zu8PMXmtml5jZajP7iov9ulTlYtcLF5f5jGteiEFRC3SLxE2tFbqqTD8kXVx6EagsKBZpzVB0xK7eMSJhUcDvqjL9kHURIYC3rZ6oJCgWTWNpgW6RuKmXTleV6Ye0i8gIiduuvxR/teli5/tMUjSNpZYDInHTCL+ryvRD2ki56uA5TBpLI3aReCng96gqmIWS21YVjUi7KODXJISR8jA5eRGJV+MCvo/Fp5sqlDsNEalGowJ+1Q9PVcnXhSyEOw0RqUajqnSqfHiqSlUvNC4izdSogN/U3i1NvZCJSLUaldJpatVJjBcyzaWIhKdRI/ym9W5Z6N6Z1nEn1AvZn009gnfdd1ApKJHANCrgN+lJ0N68fZJQL2RTB2bxyX3fWXSRUgpKpH6NSukAzak6yWqwNh5wiuTW3YdT70hCTkGJtEHjAn5TpAVHAti7dW21B1NAVlAPNQUl0haNSuk0SVpwDD1oph1flS2fRSSZAn6gYp2ATjruKls+i0g6pXQCFWvbg1iPW6QNGNjSsmeYnJy06enpug9DRCQaJPeb2WTSe0rpiIi0hAK+iEhLKOCLiLSEJm0rpP4yIlInBfyKNLlXv4jEQQG/h88ReFaLYwV8EamCAn6X7xF4jC2ORaRZnEzakryG5GGSR0huTXifJG/vvv8wyctc7Ncll4uMLLQ1vnDrZ7Fm2x5MHZiNtlWCiDRH6YBPcgTARwGsB3ARgLeSvKhvs/UAVnT/bAZwZ9n9uuZqBJ62HOHVrx6LslWCiDSHixH+FQCOmNkTZvZTAPcC2Ni3zUYAd1vHPgCjJM91sG9nXI3A0+4UvvjN443p1S8icXKRwx8H8FTP7zMArsyxzTiAY/0fRnIzOncBmJiYcHB4+WxZt/KMHD4w3Ag8606hKb36RSROLkb4THitv0FPnm06L5rtMLNJM5scGxsrfXB5uVotS7l6EQmVixH+DIALen4/H8DRIbapnYsRuKs7BRER11yM8B8CsILkhSTPAXADgJ192+wE8I5utc5qACfNbFE6pwmatK6uiDRL6RG+mZ0meQuA3QBGANxlZo+RvKn7/nYAuwBsAHAEwLMAbiy735ApVy8iIXLy4JWZ7UInqPe+tr3nZwNws4t9iYjIcPSk7QBqeCYiTaGAn0ENz0SkSdQPP4PLdgsiInXTCD9DFQ3PlDISkapohJ/B90NUaX13pg7MOvl8EZFeCvgZtqxb6bXhmVJGIlIlpXQyLKRWfKVc1CNfRKqkgD9A/0NUC73uXVwAzhtditmE4K6+OyLig1I6BbjOuftOGYmI9NIIv4BB69IWrbjxnTISEemlgF9AVs592Ie01HdHRKqilE4BWWWaqrgRkdAp4BeQlXNXxY2IhE4Bv4CsXvda6UpEQqccfkFpOXetdCUioVPAd0QVNyISOgV8h1RxIyIhUw5fRKQlFPBFRFqiVSkd9Z4XkTZrTcDXcoUi0natSenoSVgRabvWBHw9CSsibdeagK8nYUWk7VoT8NV7XkTarpGTtlnVOKrSEZG2KhXwSf48gPsALAfwLQBvMbMfJmz3LQA/AjAP4LSZTZbZb5ZB1TgK8CLSVmVTOlsBfMHMVgD4Qvf3NFeb2aU+gz2gahwRkTRlA/5GAJ/o/vwJAJtKfl5pqsYREUlWNuC/zMyOAUD3ny9N2c4A/AvJ/SQ3Z30gyc0kp0lOHz9+vPABqRpHRCTZwIBP8t9IPprwZ2OB/awxs8sArAdwM8nXp21oZjvMbNLMJsfGxgrsokPVOCIiyQZO2prZb6a9R/K7JM81s2MkzwXwdMpnHO3+82mS/wjgCgAPDnnMmVSNIyKSrGxZ5k4AvwdgW/efn+nfgOSLAJxlZj/q/vxbAP6i5H4zqRpHRGSxsjn8bQDeQPK/Abyh+ztInkdyV3eblwH4MslDAP4TwGfN7PMl9ysiIgWVGuGb2fcB/EbC60cBbOj+/ASAS8rsR0REymtNawURkbZTwBcRaQkFfBGRlqCZ1X0MqUgeB/Dtgv/aMgDf83A4IWvjOQM67zZp4zkDw533K8ws8SGmoAP+MEhO++7XE5o2njOg8677OKrUxnMG3J+3UjoiIi2hgC8i0hJNDPg76j6AGrTxnAGdd5u08ZwBx+fduBy+iIgka+IIX0REEijgi4i0RLQBn+Q1JA+TPEJy0dKK7Li9+/7DJC+r4zhdynHOb+ue68Mkv0KyET2MBp13z3a/QnKe5HVVHp8Pec6Z5FUkD5J8jOR/VH2MPuT4f/znSP4TyUPd876xjuN0ieRdJJ8m+WjK++5imZlF9wfACID/AfCLAM4BcAjARX3bbADwOQAEsBrA1+o+7grO+VcBvKT78/rYzznvefdstwfALgDX1X3cFXzXowC+AWCi+/tL6z7uis77fQD+pvvzGIAfADin7mMved6vB3AZgEdT3ncWy2Id4V8B4IiZPWFmPwVwLzrr6/baCOBu69gHYLS7SEusBp6zmX3FzH7Y/XUfgPMrPkYf8nzXAPBHAO5HyiI8kclzzr8D4AEz+w7QWVyo4mP0Ic95G4CfIUkAL0Yn4J+u9jDdMrMH0TmPNM5iWawBfxzAUz2/z3RfK7pNTIqezzvRGRXEbuB5kxwH8NsAtld4XD7l+a5fBeAlJP+9u1b0Oyo7On/ynPcdAF4D4CiARwD8sZk9V83h1cZZLCu74lVdmPBaf31pnm1ikvt8SF6NTsD/Na9HVI085/13AN5jZvOdgV/08pzz2QAuR2c9iqUAvkpyn5n9l++D8yjPea8DcBDAWgCvBPCvJL9kZs94PrY6OYtlsQb8GQAX9Px+PjpX/KLbxCTX+ZD8ZQAfA7DeOgvUxC7PeU8CuLcb7JcB2EDytJlNVXKE7uX9//t7ZvYTAD8h+SA6Cw3FHPDznPeNALZZJ7l9hOSTAF6Nzmp6TeUslsWa0nkIwAqSF5I8B8AN6Kyv22sngHd0Z7hXAzhpZseqPlCHBp4zyQkADwB4e+QjvV4Dz9vMLjSz5Wa2HMCnAfxhxMEeyPf/92cA/DrJs0m+EMCVAB6v+Dhdy3Pe30F3lT2SLwOwEsATlR5l9ZzFsihH+GZ2muQtAHajM7N/l5k9RvKm7vvb0anW2ADgCIBn0RkZRCvnOb8fwC8A+PvuaPe0Rd5hMOd5N0qeczazx0l+HsDDAJ4D8DEzSyzri0XO7/ovAXyc5CPopDreY2ZRt00meQ+AqwAsIzkD4AMAlgDuY5laK4iItESsKR0RESlIAV9EpCUU8EVEWkIBX0SkJRTwRURaQgFfRKQlFPBFRFri/wDsmy3UkiWgAQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "n = 100\n",
    "x = np.random.rand(n)\n",
    "beta = 2\n",
    "y = beta*x + np.random.rand(n) - 0.5\n",
    "plt.scatter(x, y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "bff0138b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# formally, x needs to be 2d for this to work\n",
    "# so we give it an extra dimension using reshape\n",
    "x = x.reshape(n, 1)\n",
    "\n",
    "beta_estimate = np.linalg.lstsq(x, y, rcond = None)[0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "67701179",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7fab3a957f40>"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD4CAYAAADvsV2wAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAA9QElEQVR4nO3dd3hU1dbA4d8mJBAIEKS3kKBIh1CjItKrSBMBFSleRXoPYLnIletHDGikB0QQRAEBRa6gqCAKKFW6dAiQgBB6C6Tt748JMYSZzCRzpiXrfZ48ZGbOnLN3EtbZs84+ayutNUIIIbK/XK5ugBBCCOeQgC+EEDmEBHwhhMghJOALIUQOIQFfCCFyiNyubkBGihYtqgMDA13dDCGE8Bi7du26pLUuZu41tw74gYGB7Ny509XNEEIIj6GUOm3pNUnpCCFEDiEBXwghcggJ+EIIkUO4dQ7fnISEBKKjo7l7966rmyI8XN68eSlbtize3t6ubooQTuFxAT86OpoCBQoQGBiIUsrVzREeSmvN5cuXiY6OJigoyNXNEcIpPC7g3717V4K9sJtSiiJFihAbG+vqpgjBqt0xTF53hHPX4ijt70to60p0ql3G8ON4XMAHJNgLQ8jfkXAHq3bH8ObX+4lLSAIg5locb369H8DwoC8XbYUQwoUmrzuSGuzvi0tIYvK6I4YfSwK+g3322WecO3cu9fFrr73GX3/9Zfd+o6Ki+PLLLzP9vj59+rBixYosH9eW9zuqz2n3P3jwYLOvtWvXjmvXrhl2LCEc7dy1uEw9bw8J+A6WPvjNmzePqlWr2r3frAZ8WyQmJtr1fkf12RZr167F39/fKccSwgil/X0z9bw9JOBnweLFi2nQoAHBwcG88cYbJCUlkZSURJ8+fahevTo1atQgIiKCFStWsHPnTl5++WWCg4OJi4ujSZMmqeUi/Pz8GDt2LHXr1qVFixZs376dJk2aUKFCBVavXg2YAnujRo2oU6cOderU4ffffwdg3LhxbNq0ieDgYCIiIkhKSiI0NJT69etTs2ZN5syZA5hmowwePJiqVavy7LPPcvHiRbN9atKkCW+99RaNGzdm6tSp7Nq1i8aNG1O3bl1at27N+fPnH3rPe++9R/369alevTr9+vVDa221z0uWLKFGjRpUr16dsWPHpu7Lz8+Pt99+m1q1avHEE09w4cIFAJYvX0716tWpVasWzzzzTOr2586do02bNlSsWJExY8akPh8YGMilS5eIioqicuXK9O7dm5o1a9K1a1fu3LmT5d+5EI4S2roSvt5eDzzn6+1FaOtKxh9Ma+22X3Xr1tXp/fXXX/88GDZM68aNjf0aNuyhY6Y/fvv27XV8fLzWWusBAwbohQsX6p07d+oWLVqkbnf16lWttdaNGzfWO3bsSH0+7WNAr127VmutdadOnXTLli11fHy83rNnj65Vq5bWWuvbt2/ruLg4rbXWR48e1fd/Jr/88ot+9tlnU/c7Z84cPXHiRK211nfv3tV169bVJ0+e1CtXrtQtWrTQiYmJOiYmRhcqVEgvX778oX41btxYDxgwQGutdXx8vH7yySf1xYsXtdZaL126VPft21drrXXv3r1T33/58uXU9/fs2VOvXr06wz7HxMTocuXK6YsXL+qEhATdtGlT/c0336T+LO6/PzQ0NLUv1atX19HR0Q/8TBcsWKCDgoL0tWvXdFxcnA4ICNBnzpzRWmtdvnx5HRsbq0+dOqUBvXnzZq211n379tWTJ09+qN8P/D0J4SLf/Bmtn5q0XgeO/U4/NWm9/ubP6CzvC9ipLcRUj5yl40rr169n165d1K9fH4C4uDiKFy/Oc889x8mTJxkyZAjPPvssrVq1srovHx8f2rRpA0CNGjXIkycP3t7e1KhRg6ioKMB0o9ngwYPZs2cPXl5eHD161Oy+fvzxR/bt25eaX79+/TrHjh3jt99+48UXX8TLy4vSpUvTrFkzi+3p3r07AEeOHOHAgQO0bNkSgKSkJEqVKvXQ9r/88gvh4eHcuXOHK1euUK1aNZ577jmL+9+xYwdNmjShWDFTIb+XX36Z3377jU6dOuHj40P79u0BqFu3Lj/99BMADRs2pE+fPnTr1o0uXbqk7qt58+YUKlQIgKpVq3L69GnKlSv3wPHKlStHw4YNAejZsyfTpk1j9OjRFtsnhKt0ql3GIdMw0/PsgP/xx04/pNaa3r17M2nSpIde27t3L+vWrWPmzJl89dVXzJ8/P8N9eXt7p04NzJUrF3ny5En9/n4ePSIighIlSrB3716Sk5PJmzevxXZNnz6d1q1bP/D82rVrbZ5+mD9//tR9VatWjT/++MPitnfv3mXgwIHs3LmTcuXKMWHCBKt3P5sGH+al/Vl4eXml9j8yMpJt27axZs0agoOD2bNnD0Dqzyr99mml77dMwxQ5neTwM6l58+asWLEiNRd+5coVTp8+zaVLl0hOTub5559n4sSJ/PnnnwAUKFCAmzdvZvl4169fp1SpUuTKlYvPP/+cpKQks/tt3bo1s2fPJiEhAYCjR49y+/ZtnnnmGZYuXUpSUhLnz5/nl19+sXrMSpUqERsbmxrwExISOHjw4APb3A/uRYsW5datWw/M3LHU55CQEH799VcuXbpEUlISS5YsoXHjxhm25cSJE4SEhPDee+9RtGhRzp49a7X99505cya1D0uWLOHpp5+2+b1CZEeePcJ3gapVq/Lf//6XVq1akZycjLe3NzNnzsTX15e+ffuSnJwMkPoJoE+fPvTv3x9fX98MR8yWDBw4kOeff57ly5fTtGnT1FF4zZo1yZ07N7Vq1aJPnz4MGzaMqKgo6tSpg9aaYsWKsWrVKjp37syGDRuoUaMGjz/+uNUAC6ZU04oVKxg6dCjXr18nMTGR4cOHU61atdRt/P39ef3116lRowaBgYGpKa6M+lyqVCkmTZpE06ZN0VrTrl07OnbsmGFbQkNDOXbsGFprmjdvTq1atVJH+dZUqVKFhQsX8sYbb1CxYkUGDBhg0/uEyK5URh+zXa1evXo6/QIohw4dokqVKi5qkfAUUVFRtG/fngMHDmS4nfw9iexGKbVLa13P3GsywhdCCBs4q96NI0nAF9lSYGCg1dG9ELayt96Nu5ws5KKtEEJYYU+9m/sni5hrcWj+OVms2h3joNZaJiN8IYQg41F4TCbq3aTfz534RIsnC2eP8iXgCyFyvIxSNgAKMDe9JX29G3P7scQRxdGskYAvhMjxrKVszAV7BQ/VuzG3H0tyKcWq3TFOHeVLDj+Trl27xqxZs1zdjAxLBN+3cePG1GJrYLprddGiRYa1ISoqiurVq5t9bfz48fz888+GHUsIR8qoRLGl1zQPX7DNzKg9SWun5/Il4GdSRgH//l2wRrG3THH6gN+/f3969eplb7Ns8t5779GiRQunHEsIe2VUotjSa2XMPG9pW39fb7zMlPZw1EInlkjAz6Rx48Zx4sQJgoODCQ0NZePGjTRt2pSXXnoptehZ2lHvlClTmDBhAmAqE9CmTRvq1q1Lo0aNOHz48EP7nzBhAv369aNVq1b06tWL2NhYnn/+eerXr0/9+vXZsmXLQ+/53//+R0hICLVr16ZFixZcuHCBqKgoIiMjiYiIIDg4mE2bNjFhwgSmTJkCwJ49e3jiiSeoWbMmnTt35urVq4CpTPLYsWNp0KABjz/+OJs2bQLg4MGDqSWha9asybFjxwDTSe7111+nWrVqtGrVirg40wgn7UIpgYGBqfts0KABx48fN+i3IYQxQltXwtvrwYDs7aUIbV0pU+WLLW07oUM1ki3c5OrMXL5H5/CHDx9u8232tgoODubjDIqyhYWFceDAgdTjbty4ke3bt3PgwAGCgoJSq1ya069fPyIjI6lYsSLbtm1j4MCBbNiw4aHtdu3axebNm/H19eWll15ixIgRPP3005w5c4bWrVtz6NChB7Z/+umn2bp1K0op5s2bR3h4OB9++CH9+/fHz88vtULk+vXrU9/Tq1cvpk+fTuPGjRk/fjz/+c9/UvudmJjI9u3bWbt2Lf/5z3/4+eefiYyMZNiwYbz88svEx8eTlJTEhQsXOHbsGEuWLOGTTz6hW7durFy5kp49ez7Up4IFC7J9+3YWLVrE8OHD+e677yz+nIRwifTxOOXx/bSNLfPoM9p28rojZi/iOmKhE0s8OuC7iwYNGhAUFJThNrdu3eL333/nhRdeSH3u3r17Zrft0KEDvr6mP4Kff/75geUBb9y48VBhsujoaLp378758+eJj4+32pbr169z7dq11Lo6vXv3fqBd98sQ161bN/UE9uSTT/L+++8THR1Nly5dqFixIgBBQUEEBwc/tH16L774Yuq/I0aMyLB9Qjjb5HVHSEh+MOInJOvUqZOZKV9sadvQ1pUemMEDDlzoxAKPDvgZjcSd6X5BM4DcuXOnFlCDf6pKJicn4+/vb9MnkrT7S05O5o8//kg9AZgzZMgQRo4cSYcOHdi4cWNqCimr7pceTlt2+KWXXiIkJIQ1a9bQunVr5s2bR4UKFR4qU3w/pZNe2tLEUqY4Z3OXu07Tcsa6sjZ9UtAa1q2DvXshzYpwRrE7h6+UKqeU+kUpdUgpdVApNczMNkopNU0pdVwptU8pVcfe47qKtXLHJUqU4OLFi1y+fJl79+6lpi4KFixIUFAQy5cvB0y14ffu3Wv1eK1atWLGjBmpj82dMK5fv06ZMqY/moULF1pta6FChShcuHBqfv7zzz+3WkXz5MmTVKhQgaFDh9KhQwf27dtnte1pLVu2LPXfJ598MlPvFdmHO911mpat68qu2h1Dw7ANBI1bQ8OwDZlud6faZdgyrhmnwp5ly7hm/wT7xERYuhRq14a2bWHOHLCyvkRWGHHRNhEYpbWuAjwBDFJKpV+xui1QMeWrHzDbgOO6RJEiRWjYsCHVq1cnNDT0ode9vb0ZP348ISEhtG/fnsqVK6e+9sUXX/Dpp59Sq1YtqlWrxrfffmv1eNOmTWPnzp3UrFmTqlWrEhkZ+dA2EyZM4IUXXqBRo0YULVo09fnnnnuOb775JvWibVoLFy4kNDSUmjVrsmfPHsaPH59hO5YtW0b16tUJDg7m8OHDmZ7tc+/ePUJCQpg6dSoRERGZeq/IPuwpUeBItlyYdcjJ6u5dU3CvVAlefBHu3YMFC+DwYbCw2JE9DC+PrJT6Fpihtf4pzXNzgI1a6yUpj48ATbTWD6+MnYaUR84eAgMD2blz5wMnI3chf0/OFTRujdmbmMB0I5MrUzzWUk0NwzaYvehaxt+XLeMsLx1q1o0bMHs2RETAhQtQvz68+SZ07MjW7dvZvXt3ltdvcFp5ZKVUIFAb2JbupTJA2qWKolOeeyjgK6X6YfoUQEBAgJHNE0K4WGl/X4vlBtKOmsG2KpTWZOZ6gbULs4bk+S9cgKlTYdYsuH4dWrSAL78kuXFj1qxdy+QmTdi0aRMlS5akb9++Fpc0zSrD5uErpfyAlcBwrfWN9C+beYvZE73Weq7Wup7Wut79xa6FZ4uKinLL0b1wPnOpk/SMSvEYnYKxNc9v1qlTMGgQBAZCWBi0bAk7dnDvu+9YcPo01WvUoEOHDpw+fZqPP/6YY8eOGR7swaCAr5TyxhTsv9Baf21mk2igXJrHZYFzWT2eO6/SJTyH/B05X6faZZjUpQZl/H3NjgLvM2J2jNHXCzJzA1aq/fuhZ0+oWBE++QRefhkOHeL6vHmEb9hAhQoVePXVV/Hx8WHx4sUcP36cYcOG4efnl6U2WmN3SkeZ5th9ChzSWn9kYbPVwGCl1FIgBLhuLX9vSd68ebl8+TJFihSR6X0iy7TWXL582SGjKJGxtKkTS3lxI25GMnqqZWZuwGLLFpg0Cdasgfz5YehQGDmSGKWYOnUqkZGR3Lx5k+bNm7NgwQJatmzplHhmRA6/IfAKsF8ptSflubeAAACtdSSwFmgHHAfuAH2zerCyZcsSHR1NbGysPW0Wgrx581K2bFlXNyNHc+TNSJauF9hzMskwz681fP+9KdBv3gxFisB778GgQfz1999MGT+exYsXk5SUxAsvvEBoaCh169bNcluywu6Ar7XejPkcfdptNDDI3mOBadqjtTtJhRDGc8QNU5kaNWeyXdZOJob1JzERli835eb37YNy5WDqVPSrr7J5927Ce/fmu+++w9fXlzfeeIMRI0ZQoUKFzB/HAB59p60QwjnsXdM1I5kpW5CVdpkL6ob05+5d+OwzmDwZTp6EKlVg4UKSunVj9fffE96yJVu3bqVIkSJMmDCBQYMGuXzyggR8IYRVGV0Adfac+bQj81xKkZTu4nvadlk6mdjVn+vXITLynzn0ISHw0UfcbdmSRYsX82GtWhw9epSgoCBmzpxJnz59yJcvX6b65aj7ESTgCyGsckatGVukH5mnD/b3WWtXlvpzfw79zJmmG6dat4Zx47hasyazIyOZ9sYbXLhwgbp167Js2TK6dOlC7ty2hVhHfoJKS+rhCyGssmsOuoFsXULQWrsy1Z+TJ2HgQChf3pSnb90adu3izNy5jFy9mnIBAbz99tsEBwezfv16duzYQbdu3WwO9uC8khMS8IUQVtlaa8aewmK2sOUThS2zfGyaU79vn2nefMWK8Omn0KsXHDnCvnfe4ZWICB599FGmTZtG586d2bNnDz/88APNmjXL0vRKZ32CkpSOEMKqTrXLsPP0FZZsO0uS1ngpxfN1/8mPOyslYWmqpZdSJGttc+47w9lBmzebplauXQt+fjByJHr4cDYePUr40KH88MMP5M+fn8GDBzN8+HDKly/vsH4Z/QlKAr4QOZitFwpX7Y5h5a6Y1Jx5ktas3BVDvfKPpK7m5IyLupamWk7qUiPTx3nggq7WppukhnSHLVu45/8IC1v2ZUblFtyLOUxSy3acOLSP4sWL8/7779O/f38eeeQRh/fL6MVRJOALkUNlZlRuLaA7KyVhxLz9ByQmwrJlptz8gQMQEMDeMRPplViR6P2/cuPzUSRe+xvvR8ow4O0wPnpnmEPuzja8XxZIwBcih8rMqNxaQHdWSgLsm7efKi7OVHd+8mSIioKqVWHRIi61aMFzr7xJzO+TSY67gU+pShRr8iq+FUPY6+dnd7DP6BOVIf2yQi7aCpFDZWZUbm1WS5YKi7nCtWum/HxgoKl6ZcmSsGoVp1avZsj27QQ8+ihn1y8kT+lKlHgpjJKvTCFfpadQubzs/rTiDqt9yQhfiBwqM6NyazlmZ6Uksuzvv003Ss2eDTdvQps2MG4cu/LnZ/KUKSxfvhwvLy969uzJ7kJPczVPiYd2YeunFUujeHe4eU0CvhA5VGYuFNoS0O0tkZB2300rF+OXw7H2nzxOnDClbT77DBISoFs39Jgx/BQbS/jEiaxfv54CBQowatQohg0bRpkyD5ddANs/rWR0XcQdbl4zfIlDI5lb4lAIYRxn3M5vSxvSB9j0Mj0TZ88e+OAD+OoryJ0b+vQhYfhwlu/eTXh4OHv37qVUqVIMHz6cN954g0KFCj3Upqz8XDJaBhEwbonEDDhtiUMhhGdxxoVCa2y5e9am1IfWsGmTacbN999DgQIwejS3Xn+dT9es4aM2bThz5gxVqlRh/vz5vPTSS+TJk8fsrrL6c8loFB/RPdgpUy8zIhdthRAuZWtKw+J2ycnwv//B009D48awcye8/z4Xd+3i3z4+BDRowPDhwwkICGD16tUcOHCAvn37Wgz29sjo4nb61b7K+Ptm6f4Be8gIXwjhUhktbJ5+uwckJMDSpabUzcGDplo306dz7Jln+HDWLD6rUYP4+Hg6depEaGgoTz75pIN68A9bLm678hOVjPCFEC5ly8LmD6Q+7tyBGTNMNW569TI99/nnbP/yS7pu3Eil4GAWLFhAr169OHToEF9//bVTgj38s2avv6936nN5vd0nzMoIXwjhUuZmAJmdpROYD95/31SiODYWnnqK5GnT+F4pJn/4Ib/++iv+/v6MGzeOoUOHUrJkSZf16V5icur3V+8kOKSuUFZIwBdCuEz62TAR3YMfDornzpnm0EdGwq1b0LYt8aNGseTsWSa/9RYHDx6kbNmyfPTRR7z22msUKFDANZ1J4Q7z7S2RgC+EcAmrtXyOH4fwcFi40FTzpnt3bgwaxCdbtxLRuzcxMTHUqFGDRYsW0aNHD7y9vTM6nNO4w3x7SyTgC5FDuXoOvqWR8KoF39EpbAOsWAHe3vDqq5zv1Yup335L5LPPcv36dZo2bcq8efNo3bp1lurPO5K1O5hd+XOXgC+EBzA6SNhaKdORwemBEa/WhJw9wMCty2l86k/THPrQUA63a8eURYv4vEkTEhMTef755wkNDaV+/fqGtMERMpqp46x1AyyRgC+Em3NEkLAlz+zo4FTa35dzV2/T/PgOBm79ijrnjhCbz5/ZrV+j5vDnmRwZybcffEDevHl57bXXGDlyJI8++qjdx3W0jMpQNAzb4NL8vgR8IVwgMyNnR1wEtCXP7NCLjwkJTEvYT6H5H/HYpTOcKVSCt1sO4MsC/uQ98ROH287jkUceYfz48QwePJhixYrZdzwnszTf3tX5fQn4QjhZZkfOjggStlTKdEhwunPHtD7slCnUPXOG649V5p0mY5h78za3/lxNXOxZAgMDmTZtGq+++ir58+fP+rHckDPXDTDHfe4IECKHyGjkbI61WvRZYUv9ekOPe/Uq/Pe/prthhw6FgACuLV3K7Fd78enmRcSum0nlskVZsmQJx44dY8iQIdku2IPr1w0wJOArpeYrpS4qpQ5YeL2JUuq6UmpPytd4I44rhCfK7MjZEUHClrouhhz33DkIDYWAAPj3vyEkhOiVKxkdEkK5117jzbfeokaNGvz000/s2rWLHj16kDt39k08uLqejlE/2c+AGcCiDLbZpLVub9DxhPBYmf1Y76jFRazVdbHruMeOmebQL1oESUnQowcHOndm8urVfNm9O1prunfvzujRo6ldu7Zd/fA0rqynY0jA11r/ppQKNGJfQmR3mVl45D5XBYlMH3fXLlMxsxUrwMcH/a9/8VujRoQvXszarl3Jly8fAwcOZMSIEQQGBjqs3cI8Z352elIptRc4B4zWWh80t5FSqh/QDyAgIMCJzRPCORwxYnfpTVRaw8aNprVif/oJChYkKTSUbypWJHzuXHbMnk2xYsWYOHEiAwYMoEiRIs5pl3iIYStepYzwv9NaVzfzWkEgWWt9SynVDpiqta5obZ+y4pUQ1llaks/hueHkZFi92hTot2+HEiWIGzyYhfnzM2XmTE6cOMGjjz7K6NGj6d27N76+zpmJktO5fMUrrfWNNN+vVUrNUkoV1VpfcsbxhcjOMjNf3pBPAvHx8OWXptTN4cMQFMSVKVOYdeMG06ZNIzY2lvr16/PBBx/QqVMnvLwyLn0snMcpAV8pVRK4oLXWSqkGmGYHXXbGsYXI7myd9WP3nbO3b8O8efDhh3D2LNSsSdTUqUQcO8a88eO5c+cO7dq1Y8yYMTzzzDNuV+NGGBTwlVJLgCZAUaVUNPAu4A2gtY4EugIDlFKJQBzQQ7vz6ulCeBBbZ/1k+c7ZK1dMC45MmwaXL0OjRuwZM4bJf/zBspEjUUrx0ksvMXr0aGrUqGFIn4RjGDVL50Urr8/ANG1TCGEwW2f9ZPrO2ehoUx36OXPg9m30s8+yvkULwteu5achQ/Dz82P48OEMGzaMcuXKGdYf4TjZ9w4HIVzMWTNnbJ31Y/P8/yNHYPJk0xz65GQSu3VjRXAw4UuXsnvECEqWLMmkSZPo378//v7+hvdHOI4EfCEcwNllcG2ZL2/1k8CuXaYZN19/DXnycLtvX+aXKsVHCxcStWQJZQIf5dEuo0gMepr/UZDKp27Tqba/xeO5ut6+eJgEfCEcwB2XuTP7SaDV43S6egRa9oGff4ZChYgdOpQZuXMzY8ECrly5wlNPPUWPof/mm8sluZtouvRm7QTm6rrv9squJysJ+EI4gFGVJo0OPKmfBJKTYdUq6N8FduyAkiU5MWYMH16+zII5c7h79y4dO3YkNDSUhg0b0jBsA3cTH2x72oJv6dvojic8W3n6ySojEvCFcAAjyuA6JPDEx8MXX5jm0B85Ao8+yo4332TykSOsnDKF3Llz88orrzBq1CiqVKmS+jZLJ6r7bUrfxvTB3tp+3Iknn6yskfLIQjiAEZUmM1tGOUO3bsHHH8Ojj8Krr6Lz5uWHceNoVq4cDSZN4sf16xkzZgxRUVHMmzfvgWAPlk9UXkqZbaOXhTn4zqr7bg9zJ+qMnvckEvCFcAAjyuAakha6fBkmTDDVoR8xgoTAQD4fPZpaycm0DQvj6LFjTJkyhTNnzjBp0iRKlSpldjeWTmBJFm6nSdLapXXf7WHpZGXpeU8iKR0hHCTtzJn7ufgRy/bYnIu3Ky0UHW26I3buXLhzh5vt2jHvsceI+OYbzm7eTLVq1fjss8948cUX8fHxsakvYD5Xb66NZdK87m4XPq1dF8noJObpJOAL4WBZzcVnpYwyhw+b6tAvXgzJyfzduTPT/P2ZvWIF19aupXHjxkRGRtK2bdtMlz6wNPXTUhtdWffdElt+F2UsnGjLeEA6yhoJ+EI4WFYvAmaqjPKOHRAWBt98Q5KPDzOqPcOkewlc/OZbSE6kS5cuhIaGEhISYmjfHLU4S0bsmblky+8iSydaDyEBXwgHsycXn+EoWWtYv94U6NevB39/lrZ/gUF/neTKng3glRu/Gi0o/tTz9HytDSEOCsLOHMnbO3PJlt+FK05iziIBXwgHs5SL10DDsA2ZDyZJSWyP+JQCH39ElZgjXPR7hK86dGNZ7Fk2/+8rcuX1o9CT3ShQtz1e+QuTBNliSiHYP2XS1usi7piOMoLM0hHCwczNcLnv/gh11e4Y6zuKj4f587n5WCUahL6Bun2ddtWaUc4nH0NWf8WRU2d4pPnrlBmwAP9nXsErf+HUt7rL/PdVu2NoGLaBoHFraBi2wbZ+p2HvzCVHLAjvSWSEL4SDpU0RmBtdxiUkMeqrvZZn8Ny6BZ98Ypp1ExPDkZJBDKnahB2n95J0cAPexYMo+txoHg1pSS6v3Hbf8OUoRtxIZu8Nbdk5XWMLCfhCOIi5i4sjlu3B3OS++1P+HgiC5fKY6tBPnw5XrhDz1FNMfeoppnz7HfrvU+QtX4siz44gb2BtlFL8fTOBiO7V3PaCoxF3sBpxQTW7pmtsIQFfCAewNJr1z+fN1TsJGb638OXzxA/+BPb8AHfucLBZM6b4+vLFjz+StHUrRao3JnftjuQp+dgD7yvt7+vWI1gjbiRz5/55Agn4QjiApdFsnty58PX2Mltr5tFLZxmwbQUd/9qI1ppNrVsRHhfHdxs24OvryxtvvMHIkSPZe80nw1Guu45gjagvBO7bP08gF22FcABLo9brcQkPlFzwUorgc0eY8/V/Wf/pANoc2sSIwNqULFWBZ374ga0HD9LjjZFUH/UF3+VvQ89lpwDsLtvgCjn9gqk7kBG+EA6Q0Wi2U+0ydAouDT/9ROzb/6HYzt+54JOP7hXq8e2VaO6d3EmJsuWZOXMmRYJbMeH74w+lhiZ1qcGWcc2c3S27SDrG9ZQ7ryVer149vXPnTlc3Q4hMS5/DB9NodlLHqnQ6tc10s9Sff3KlZEkmVKjMnL17iL99jfxlHmfAkBFMGvUauXPnpmHYBou3+XtawBfOoZTapbWuZ+41GeEL4QDpR7Pl/byIiN9P7ZeGwrFjnAkMJKJFCz754w9u/76RNm3aMGbMGJo0afJAjRujFlIRAiTgC+EwnWqXodNjBU0VKz/6CM6dY1/lykxu1Iglv/+Oio6mR48ehIaGUrNmTbP7MOpCpxAgF22FcIzYWBg/HsqXR48ezYbixWlTrx61Dh/mmz//ZOjQoZw4cYLPP/+ck0lFLN59Khc6hZFkhC+Ekc6cgSlTYN48EuPi+Lp+fcJv3mTXnj0UL16c999/nwEDBlC4sKnsgbW7T+VCpzCSBHxhKKMX3fYYf/1lWif2yy+5ozUL6tfno5gYTu7YQcWKFZk7dy6vvPIKefPmfeBtttx9auS88xz7+xGABHxhIIcsuu3utm2DSZPg22+55OvLzHr1mH7kCJe3biUkJIQpU6fSoUMHvLzMF0+zdlHWyADtyN+PnEg8gyE5fKXUfKXURaXUAQuvK6XUNKXUcaXUPqVUHSOOK5zLWqVDQxfddmdaw7p10LQpPPEEpzZuZEj9+gQAE7Zu5amnn+a3337jjz/+oHPnzhaDPVi++Fra3zc1QMdci0OTycqaZjjq92N0O4XjGHXR9jOgTQavtwUqpnz1A2YbdFzhJLb8p3bXKYT2luRNlZQEX30FdetCmzbsOniQHsHBPHbzJnP27KFHjx4cPHiQ1atX06hRI5uWEMzooqzRAdpRv58cc6LPBgwJ+Frr34ArGWzSEVikTbYC/kqpUkYcWziHLf+pMxqtuooho89790xTKytXRnfvzrqLF2lRpQr1YmNZe+IEo0aN4tSpU8yfP5+qVatmqn2dapexWCbB6ADtqN+Pu57oxcOclcMvA5xN8zg65bnz6TdUSvXD9CmAgIAApzROWGfLf2p3XAvUrpK8N27AnDkQEUHC+fN8FRjI5MBA9kZFUVprwsPD6devH4UKFbKrjZYuyho9B99Rvx+5V8BzOGsevrnPtmZrOmit52qt62mt6xUrVszBzRK2smV0mNFo1VWyNPq8eBHeeQfKl+fWmDFM9fPjseLF6RkVRbyvL/Pnz+fkyZOEhobaHewzYvQcfEf9fuReAc/hrBF+NFAuzeOywDknHVsYwNbRobuVrs3U6DMqyrSq1KefciEujumVKzMrOZmrx47RqFEjZn76Ke3atSNXLueMkxwxB98Rvx+5V8BzGFY8TSkVCHynta5u5rVngcFAOyAEmKa1bmBtn1I8zb144tQ7i0XMutQATEHK79ghhu/+ltb7NnAC+PCxx/js1CniExLo1KkToaGhPPnkky7qgRCZ4/DiaUqpJUAToKhSKhp4F/AG0FpHAmsxBfvjwB2grxHHFY5lLsB7WoVGS6NPgGXTvmLC5qW0PL6djV4+hBQswa5r5/GJiqJ3nz6MGjWKxx9/3JXNF8JQUh5ZmJXRyNjdR/UZSplDv3vgWGqd2sdXPr68nbcAJ29cJFee/JR6siM7l3xIyZIlXd1SIbJEyiOLTDNiwWm3kpgIK1ZAWBjxe/fya14/WucvTOztq3jl8aNws9fwq9kKrzz5rAZ7T0xtCQES8IUF2WZu9d27sHAhTJ7MjRMnmF64CB/k9ePm3Vt4FwukSNNXyV+5EcrL9F/B2lTCHFk+QmQbEvCFWR4/t/rGDYiMhIgIzv39N9NKlmR2vnzcuHqZPAE1KR7yPHmD6jxwN6wtUwmz3ScfkaNIwBdmueNNVDa5cAGmToVZszh0/TpTSpfm89y5Sbp4ka5du7Ip75P4lKpo9q22XJ/INp98RI4kC6AIs9zxJqoMnToFgwZBYCBbJk2iY758VAW+vHKF1/v14+jRoyxbtoygKuZXliqTsri4Ne5YPkIIW8kIX1hk7SYdt7h4uX8/fPAByUuWsFopJhcpwu937/LIvXu8++67DBo0iLR3bNv7ycVjP/kIgQR8kUUuv3i5ZQtMmsTdNWtY7OPDlEKFOHL1KoH58jF9+nT69u1L/vz5H3qbvXeFyl2lwpPJPHyRJQ3DNpi9qFvG39dxN2dpDd9/D5MmcW3zZiLz5WNqrlz8fesWtWvXZsyYMXTt2pXcuWUcI3IumYcvDOfUi5eJibB8OYSFcXbfPj7282NunjzcunOHVq1asXjMGJo1a2ZT/Xl35BapMZEjSMAXWeKIaZvpA9/YJuXpsPtHmDyZ/adOMaVQIb708kLHxdG9e3dCQ0MJDg62oxeu5/LUmMhRJOCLLDH64mXawFfg3m06/LCcJyZ+y8Y71wgvVIjvgXwJCQwaPJgRI0ZQvnx5g3pivMyM2GVev3AmCfgiS4y+eDl53RHyX7vEkJ3f8uKuNaxPiKOxTz4OA8V8fJg4cSIDBgygSJEimd63M1MmmR2xy7x+4UwS8EWWGVZb/dQp+n/1Ic/t/ZEvkpOo5p2H80Du/IV4pEkfTq+agq+v9VSRucAOODVlktkRu8ff0Sw8igR8D5EtL+zt3w9hYVxeupQTGsp75eZGciI+RctRNKQr+So+QdlH/GwO9uYCe17vXE5NmWR2xC7z+oUzScD3ANnuwt7mzRAWRtSaNUR4ezMvVy7uJCaSP7AWJep3IU+56iilMhX4LI2s0z93n6NSJpkdscu8fuFMEvA9QLa4sKc1rF0Lkyaxe8sWJufJw1e5cqG05uWXX2b06NEcTyic5cCX2QDuqJRJVkbs7rYspMi+JOB7AEvBLOZaHA3DNrj3yDAxEZYtQ4eFsf7AAcLz5uUnwNcrNyUadCBX9Wc5FlCO4wmF7Qp8lkbW/r7e3EtMdlrKREbswp1JwPcAloKZgtTn3S7NExcHCxaQGB7O8tOnCc+blz1ASX9/XunWl60+dYj3Mo2yjWi7pZH1hA7VAOcG4IxOXNnyWozwGBLwPYC5YKaA9EUx3CLNc+0azJ7N7YgI5sfG8lGePEQBlcqXZ15oKD179qRZxBbi053A7G27tZG1OwTVbHctRngcCfgewFwwMzfiBxfO3/77b/j4Y2JnzmTGrVvM8PbmCtCwXj2mjhlD+/btyZUrV4ZttLft7p4LzxbXYoRHk4DvIdIHM0vFy5w+f/vECZg8meMLFvBRfDwLvLy4C3Rs147Q0FAaNmxoto1u0fYMOCL1IjdZCVeTBVA8VGjrSvh6ez3wnL2lDRqGbSBo3Boahm1g1e6YjN+wZw+8+CI7Klak2yefUCkhgU+9venZty+HDh1i1apVZoO9I9putPupl5hrcWj+Sb1Y/ZlYIYunCFeTgO+hjFyRyuYApzX89hu6bVu+r12bpsuX00BrfvTzY8zYsUSdPs0nn3xC5cqVndZ2R8go9WIPdz/RiexPUjoezKictdXccnIyrFlD/P/9H0u3bmWKlxf7gTLFizNl1Chef/11ChYs+MD7raVE3Dnf7shrDCBTNoXrSMDPBuzNN1sKZBcv34TFi7n5f//HJ4cOEeHlRTRQ7fHHWThuHD169MDHx8dse9LPRgldsZcJqw9yPS7B7QOdI68xuPOJTmR/EvA9nBFT/dIHuLwJd3lh/888v30lb12PZVauXFwHGjdsyJyxY2nbtm2Gi42Y+8SQkKS5FpeQ5TY6k9S3EdmVITl8pVQbpdQRpdRxpdQ4M683UUpdV0rtSfkab8Rxs4tMXzBNw1I6ZtRXe23e3/3ccsG7txj0+zI+m9WHv3+KJOTGJcKUomXnzmzbto2Nv/5Ku3btrK4sZUvqw4icuKO4+zUGIbLK7hG+UsoLmAm0BKKBHUqp1Vrrv9Jtuklr3d7e42U39o7QLQXXpJS1im3ZX6eSuah25lvOLlnA9IR7jAW8c3vz6r9eZeSoUVSsWDFTfcroPgFb2u4OJPUisiMjRvgNgONa65Na63hgKdDRgP3mCPbOCLElr2xxf8ePk9yvH/8LCODVRZE0T7jHrwUL8vY773Am+iyzIyMzHezB/GyUrLZdCGEcIwJ+GeBsmsfRKc+l96RSaq9S6nulVDVLO1NK9VNK7VRK7YyNjTWgee7N3hkhtgbXB/a3ezf3XniB+Y8/TrV58+iQmMiZ0qWZOnUqZ2JimDhxIiVKlLDp+OakT4nk93m4fZITF8L5jAj45hK66cu8/AmU11rXAqYDqyztTGs9V2tdT2tdr1ixYgY0z71ZGuX65/O26f3pg6uXhfx66UJ54ddfud68OeF16hC0ciX/0po8VavyxRdfcDwqiqFDh+Ln55fVrjzUri3jmhHRPZjkdH8NCni+rqRMhHA2I2bpRAPl0jwuC5xLu4HW+kaa79cqpWYppYpqrS8ZcHyPFtq6EqEr9pKQ9GBUvH4ngVW7Y2wKimnzzemvCSidTNuoXQzct5LQtw4wRyluAi0aN+azN9+kZcuWVi/CppXZKaDmUlYa+OVw9v/0JoS7MSLg7wAqKqWCgBigB/BS2g2UUiWBC1prrZRqgOmTxWUDju0Qzixh26l2Gd78et9DAT8ZmLD6YKaPe3/7j9YepP4fP9Dqj6UsvXKeOkCSUnTr2pXQceOoU6dOptualQvMUj9GCPdhd8DXWicqpQYD6wAvYL7W+qBSqn/K65FAV2CAUioRiAN6aK3Tp33cgjNL2N4/scQlJJt9/f689Uy5c4eOm1byyLT3Cb94kY8AXx8f+r/+OiNGjSIoKCjL7c1KtUdPKJQmRE5hyI1XWuu1wNp0z0Wm+X4GMMOIYzmas0rYpj+x2O3qVZKmT+fbKVMIv3mTbUDRggX5z8iRDBw0iKJFi9p9iKyM1uUmJiHch9xpm46zUhDmTizpFbblwu25c8SFh7MoMpIP793jGFChdGlmvfMOvXv3Jl++fFZ3YWsKKyujdakfI4T7kICfjrNSENZOIN5einefszh7FY4d48p77zF7yRKmJSVxEahbtSpfTZhAly5d8PKyPlUTMpfCyupoXW5iEsI9SHnkdJxVwjajE0gZf18md61lPkj++Sdn2rdnxOOPE7B4Me8kJVHnmWfYsGEDOw4c4IUXXrA52EPmbvySkgNCeDYZ4afjrBSEpdGy2QCqNfz6K3vfeovJf/zBUkDlysWLXbow+t//pmbNmlluR2ZTWDJaF8JzScA3wxlBzaYTS3IyevVqfnnzTcIPH2YdkN/Hh6GvvcbwsWMJCAiwux0yi0aInEMCvgtZPLEkJJC4eDEr//1vwmNi+BMoUbAg/zd6NP0HD6Zw4cKGtUFm0QiRc2TbgO/Mm6cMc/s2d2bPZsH77/PhtWucAh4vVYq548fzSp8+5M2b1/BDyiwaIXKObBnwnXnzlCGuXOHSBx8wY/p0ZsTFcRl4onJlPnz/fTp07Jh6EdZRJzHJywuRM2TLWTqOWoTacNHRnPzXvxhcogQB4eH8Jy6Opxo2ZNOmTfz+1190TjO90uaFxoUQwoJsGfDdvn7LkSPs6tSJHgEBVJw/n7nJybzYuTMHDx5k9ebNPP300w8VNPOYk5gQwm1ly5SOu8480Tt38uPw4YRv2cIGoKCPD6P79GHo+PGUKZNxSsXtT2JmeOR1FCGysWw5wnfWzVM20ZqEH3/ki+rVCa5fnzZbtnDYz4/wf/+bMxcv8sGcORkG+/vr3VqqNOfqk5gl76zaz4hleyQFJYQbyZYjfLeYeZKczK0lS5j35ptEnD3LGaBq8eIsePddXnrtNXx8fKzuwlqBNXedPrlqdwxfbD3z0EnKEUXohBC2y5YBH1w48yQ+nguzZjF94kRmXbnCVaDRY48x84MPaNepE7ly2f6hKqMCa2XcOEUyed0Ri59I3DkFJUR2l20DvtPdvs3R99/nw2nTWHj7NvFA5wYNCP3wQ554+uks7dJScFTAlnHNst5WB8soqLtrCkqInCBb5vCd6vJltr3+Os8XLkzlSZNYeOcOvdu25fDhw6zcti3LwR4sB0d3D5qW2qfALVNQQuQUEvCzKPnMGb7r1InGxYvzxLx5bADe7NWLqHPnmLN2LY9Xsj+wudXF50ww124FvPxEgFumoITIKSSlk0nx+/fz5aBBTN60ib+AcvnyETFkCP96+20KFChg6LHc4uJzFnhqu4XI7pSbLi0LQL169fTOnTtd3QwAbvzyC3OHDuXjAweIAWoUKcKYd96h+6BBeHvbsDKVEEI4gVJql9a6nrnXZISfEa05t2wZU8eNI/L0aW4AzQID+TQsjFbduj10N6wQQrgzCfjmJCVxaMYMpvz3v3x+6RJJQNdatQj9+GPqNWni6tYJIUSWSMBPKz6eLe++S/iMGay+dQtfpejXpAkjZ82iQpUqrm6dEELYRQI+kHzjBquHDyf8iy/4Iz6eIl5evNu1K4OmT6dYyZKGH09qzAghXCFHB/y70dEsHjCAyWvXcjQ5maC8eZk+cCB9P/iA/H5+Djmmx9XqF0JkGzky4F/dv5/I/v2Z+vvvXADqFCrE0lGjeP7NN8mdO3fKCHy7Q0bgGZU5loAvhHCkHBXwz27YwMdDhzL34EFuAa1Kl2bMxIk069s3dcaNo0fgnljmWAiRPRhyp61Sqo1S6ohS6rhSapyZ15VSalrK6/uUUnWMOK6t9n/5Jb0CAqjQvDlTDx6kY6VK7F6zhnUxMTR/9dUHplcaudDI/dLGQePW0DBsA6t2x3hsuQQhhOezO+ArpbyAmUBboCrwolKqarrN2gIVU776AbPtPa41OjmZjVOm0K5IEWq+/DIrz55lUIMGnNi1i8WHDxPcrp3Z9xk1Are0JGHTysU8slyCEMLzGTHCbwAc11qf1FrHA0uBjum26Qgs0iZbAX+lVCkDjv2QpPh4VowaRUiBAjQNDWXntWv8t00bzkZF8fG2bZSvk/GHC6NG4JY+KfxyOJZJXWpQxt8XhanM8aQuNSR/L4RwOCNy+GWAs2keRwMhNmxTBjiffmdKqX6YPgUQEBCQ6cbEXblCv4gIiuTOTWTPnvSaPh1ff3+b3x/autJDi45kZQSe0ScFl9XqF0LkaEYEfHP1BdIX6LFlG9OTWs8F5oKplk5mG+NXsiSbV62iUps2eNmwqlR6RhX+ctd1dYUQOZcRAT8aKJfmcVngXBa2MUzVDh3ser8RI3CjPikIIYRRjMjh7wAqKqWClFI+QA9gdbptVgO9UmbrPAFc11o/lM7JTjrVLiO5eiGEW7F7hK+1TlRKDQbWAV7AfK31QaVU/5TXI4G1QDvgOHAH6GvvcT2B5OqFEO7EkBuvtNZrMQX1tM9FpvleA4OMOJYQQoisyVF32tpDCp4JITydBHwbSMEzIUR2IIuY28DIcgtCCOEqMsK3gaMLnkm6SAjhDDLCt4EjC55ZqrmzaneM3fsWQoi0JODbILR1JYcVPJN0kRDCWSSlYwOjyi2YI/XxhRDOIgHfRulvorpf697eE4DU3BFCOIukdLLAyLy7I9NFQgiRlozws8DaurSZmXXjyHSREEKkJQE/CzLKu2flJi2puSOEcAZJ6WRBRtM0ZdaNEMJdScDPgozy7jLrRgjhriTgZ0FGte4deZOWEELYQ3L4WWQp7y4rXQkh3JUEfIPJrBshhLuSgO8AMutGCOGOJIcvhBA5hAR8IYTIIXJkSkfqzwshcqIcF/BluUIhRE6V41I6ciesECKnynEBX+6EFULkVDku4MudsEKInCrHBXypPy+EyKmy9UXbjGbjyCwdIUROY1fAV0o9AiwDAoEooJvW+qqZ7aKAm0ASkKi1rmfPcW1hbTaOBHghRE5jb0pnHLBea10RWJ/y2JKmWutgZwR7kNk4QgiRnr0BvyOwMOX7hUAnO/dnGJmNI4QQD7I34JfQWp8HSPm3uIXtNPCjUmqXUqpfRjtUSvVTSu1USu2MjY3NcsNkNo4QQjzIasBXSv2slDpg5qtjJo7TUGtdB2gLDFJKPWNpQ631XK11Pa11vWLFimXiEA+S2ThCCPEgqxdttdYtLL2mlLqglCqltT6vlCoFXLSwj3Mp/15USn0DNAB+y2KbbSKzcYQQ4kH2TstcDfQGwlL+/Tb9Bkqp/EAurfXNlO9bAe/ZeVybyGwcIYT4h705/DCgpVLqGNAy5TFKqdJKqbUp25QANiul9gLbgTVa6x/sPK4QQohMsmuEr7W+DDQ38/w5oF3K9yeBWvYcRwghhP1yXGkFIYTIqSTgCyFEDiEBXwghcgiltXZ1GyxSSsUCp7P49qLAJQOb4wmkz9lfTusvSJ8zq7zW2uxNTG4d8O2hlNrprLo97kL6nP3ltP6C9NlIktIRQogcQgK+EELkENk54M91dQNcQPqc/eW0/oL02TDZNocvhBDiQdl5hC+EECINCfhCCJFDeHTAV0q1UUodUUodV0o9tLyiMpmW8vo+pVQdV7TTSDb0+eWUvu5TSv2ulPL4OkbW+pxmu/pKqSSlVFdnts8RbOmzUqqJUmqPUuqgUupXZ7fRaDb8bRdSSv1PKbU3pc99XdFOoyil5iulLiqlDlh43fj4pbX2yC/ACzgBVAB8gL1A1XTbtAO+BxTwBLDN1e12Qp+fAgqnfN82J/Q5zXYbgLVAV1e32wm/Z3/gLyAg5XFxV7fbCX1+C/gg5ftiwBXAx9Vtt6PPzwB1gAMWXjc8fnnyCL8BcFxrfVJrHQ8sxbTGblodgUXaZCvgn7JQi6ey2met9e9a66spD7cCZZ3cRqPZ8nsGGAKsxMIiPB7Glj6/BHyttT4DpsWFnNxGo9nSZw0UUEopwA9TwE90bjONo7X+DVMfLDE8fnlywC8DnE3zODrlucxu40ky259/YRoheDKrfVZKlQE6A5FObJcj2fJ7fhworJTamLJWdC+ntc4xbOnzDKAKcA7YDwzTWic7p3kuYXj8snfFK1dSZp5LP8fUlm08ic39UUo1xRTwn3ZoixzPlj5/DIzVWieZBn8ez5Y+5wbqYlqPwhf4Qym1VWt91NGNcxBb+twa2AM0Ax4FflJKbdJa33Bw21zF8PjlyQE/GiiX5nFZTGf+zG7jSWzqj1KqJjAPaKtNi9R4Mlv6XA9YmhLsiwLtlFKJWutVTmmh8Wz9276ktb4N3FZK/YZpoSFPDfi29LkvEKZNCe7jSqlTQGVMK+llR4bHL09O6ewAKiqlgpRSPkAPTGvsprUa6JVytfsJ4LrW+ryzG2ogq31WSgUAXwOvePBoLy2rfdZaB2mtA7XWgcAKYKAHB3uw7W/7W6CRUiq3UiofEAIccnI7jWRLn8+QssKeUqoEUAk46dRWOpfh8ctjR/ha60Sl1GBgHaYr/PO11geVUv1TXo/ENGOjHXAcuINphOCxbOzzeKAIMCtlxJuoPbjSoI19zlZs6bPW+pBS6gdgH5AMzNNam53e5wls/D1PBD5TSu3HlO4Yq7X22LLJSqklQBOgqFIqGngX8AbHxS8prSCEEDmEJ6d0hBBCZIIEfCGEyCEk4AshRA4hAV8IIXIICfhCCJFDSMAXQogcQgK+EELkEP8P1qQVZmlTVgQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(x, y)\n",
    "plt.gca().plot(np.linspace(0, 1), beta_estimate*np.linspace(0, 1), color = \"red\", label = \"estimated relationship\")\n",
    "plt.gca().plot([0,1], [0,2], color = \"black\", label = \"true relationship\")\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "79ce780c",
   "metadata": {},
   "source": [
    "This works in any number of dimensions! \n",
    "\n",
    "In fact, this least-squares problem has an analytic solution in terms of matrix inverses as well, given by \n",
    "\n",
    "$$\n",
    "\\hat{\\beta} = (\\mathbf{X}^T\\mathbf{X})^{-1}\\mathbf{X}^T \\mathbf{y}\\;.\n",
    "$$\n",
    "\n",
    "This matrix computation is easy to implement in Python using tools from the previous section. Rather than creating a plot, let's show how we perform multidimensional regression and recover an estimate of the coefficient vector $\\beta$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "81c8e93d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([0, 1, 2]), array([0.27025418, 1.24981514, 2.19446116]))"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def least_squares(X, y):\n",
    "    return np.linalg.inv(X.T @ X) @ (X.T) @ y\n",
    "\n",
    "# multidimensional data\n",
    "X = np.random.rand(100, 3)\n",
    "\n",
    "# true value of beta is [0,1,2]\n",
    "beta = np.arange(3)\n",
    "y = X @ beta + np.random.rand() - 0.5\n",
    "\n",
    "beta_hat = least_squares(X, y)\n",
    "beta, beta_hat"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "67a66791",
   "metadata": {},
   "source": [
    "The reason that `numpy` implements a specialized least-squares function like `lstsq` is that the analytic approach can be impractical to compute when the number of columns in `X` is large, as matrix inversion is a very slow operation. "
   ]
  }
 ],
 "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"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}