Ordinary Least Squares Optimization.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h1 style=\"background-color:#0071BD;color:white;text-align:center;padding-top:0.8em;padding-bottom: 0.8em\">\n",
    "Ordinary Least Squares Optimization\n",
    "</h1>\n",
    "\n",
    "Accompanying notebook for the recipe:\n",
    "\n",
    "  * Christian Bauckhage: \"NumPy / SciPy Recipes for Data Science: Ordinary Least Squares Optimization\",<br/> Technical Report, March 2015 ([Download from ReseachGate](https://www.researchgate.net/publication/273133972_NumPy_SciPy_Recipes_for_Data_Science_Ordinary_Least_Squares_Optimization))\n",
    "\n",
    "  * Abstract of the paper: *In this note, we study least squares optimization for parameter estimation. By means of the basic example of a linear regression task, we explore different formulations of the ordinary least squares problem, show how to solve it using **NumPy** or **SciPy**, and provide suggestions for practical applications.*\n",
    "\n",
    "In this notebook you find now in addition to the technical report as  `lsq_solution_V0(X, y)` a direct implementation of the formula $(X^TX)^{-1}X^Ty$ based on numpy matrices reading `(X.T * X).I * X.T * y`. It comes with a slight performance penality, especially for small matrices: 200% for n = 100, but only 3% for n = 1000000. `X` and `y` in this case have to be instantiated with `np.matrix(...)` not `np.array(...)`.\n",
    "\n",
    "As an experiment we write `ŷ` ([Unicode 0177](https://www.fileformat.info/info/unicode/char/0177/index.htm)) instead of `yhat` to be closer to $\\hat y$. This is not yet a recommendation. Try for yourself. There are [different strongly](https://stackoverflow.com/questions/2649544/unicode-identifiers-in-python) [diverging opions](https://softwareengineering.stackexchange.com/questions/16010/is-it-bad-to-use-unicode-characters-in-variable-names) on this topic on the net. The benefit of close resemblence to the mathematical formula competes with the difficulty of typing the character.\n",
    "\n",
    "<p style=\"background-color:#66A5D1;padding-top:0.2em;padding-bottom: 0.2em\" />"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import numpy.linalg as la\n",
    "import numpy.random as rnd\n",
    "\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import warnings\n",
    "warnings.simplefilter(action='ignore', category=FutureWarning)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def create_data(n, xmin=-2, xmax=12, a=1.1, b=2.0):\n",
    "    x = rnd.random(n) * (xmax - xmin) + xmin\n",
    "    y = a * x + b + rnd.randn(n) * 0.5\n",
    "    return x, y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "x, y = create_data(25)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(7.83, 10.79), (3.86, 5.80), (9.66, 13.30), (3.65, 6.79), (10.35, 13.58), (4.50, 7.27), (6.26, 9.46), (2.06, 4.08), (-1.14, 0.59), (9.98, 13.83), (10.82, 13.97), (-0.29, 1.40), (10.10, 13.07), (11.35, 13.58), (4.62, 6.38), (3.00, 4.97), (2.52, 4.86), (-0.91, 0.77), (11.86, 15.56), (6.54, 9.61), (9.84, 13.14), (-0.52, 2.01), (2.70, 4.18), (9.31, 12.38), (-0.38, 1.14), "
     ]
    }
   ],
   "source": [
    "for xi, yi in zip(x, y):\n",
    "    print('({:.2f}, {:.2f})'.format(xi, yi), end=', ')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAEwRJREFUeJzt3X9sXXd5x/H3QxKGyw+5KIYRt1nCBNlYUxZkpkI0xigsHVQ0iphEt6JuIEWaBhTEUhLxR7V/lkhBDCQmUNSGdlrVqSrBMAqkVQOrJkGHUxfSEgKIHyVOIUYlgKi3pumzP3zdOo7t++Oc43vvue+XVNn3+Pqep1X60TfP+f6IzESS1P+e0+0CJEnlMNAlqSYMdEmqCQNdkmrCQJekmjDQJakmDHRJqgkDXZJqwkCXpJpYvZI3W7t2bW7YsGElbylJfe/o0aO/yMyRZu9rGugRcRC4GjidmZfNu/4+4L3AU8DdmXljs8/asGEDExMTzd4mSZonIn7SyvtaabncCly14MP/HLgGuDwz/wj4aLsFSpLK1TTQM/N+4PEFl/8e2JeZ/9d4z+kKapMktaHTh6KvBP40Ih6IiP+KiNeWWZQkqX2dPhRdDVwMXAG8FrgzIl6ei+zFGxE7gZ0A69ev77ROSVITnY7QTwKHctb/AE8Daxd7Y2YeyMyxzBwbGWn6kFaS1KFOR+jjwJuAr0XEK4HnAr8orSpJqoHxySn2Hz7BqTMzrBseYte2TWzfMlrZ/VqZtngH8EZgbUScBG4CDgIHI+Jh4Eng+sXaLZI0qMYnp9hz6BgzZ88BMHVmhj2HjgFUFupNAz0zr13iR9eVXIsk1cb+wyeeCfM5M2fPsf/wicoC3aX/klSBU2dm2rpeBgNdkiqwbnioretlMNAlqQK7tm1iaM2q864NrVnFrm2bKrvnim7OJUmDYq5P3lOzXCRJndm+ZbTSAF/Ilosk1YSBLkk1YaBLUk0Y6JJUEwa6JNWEgS5JNWGgS1JNGOiSVBMGuiTVhIEuSTVhoEtSTRjoklQTTQM9Ig5GxOnGcXMLf/aPEZERsegB0ZKkldPKbou3Ap8E/m3+xYi4FHgL8Gj5ZUlS9Vb6EOeqNR2hZ+b9wOOL/OhfgBsBD4eW1HfmDnGeOjND8uwhzuOTU90urWMd7YceEW8HpjLzWxFRckmStLgyR9TdOMS5am0HekRcBHwE+IsW378T2Amwfv36dm8nScCzI+q5EJ4bUQMdBXA3DnGuWiezXH4f2Ah8KyJ+DFwCPBgRv7vYmzPzQGaOZebYyMhI55VKGmjLjag70Y1DnKvWdqBn5rHMfElmbsjMDcBJ4DWZ+bPSq5OkhrJH1J0c4jw+OcXWfUfYuPtutu470nP99qYtl4i4A3gjsDYiTgI3ZeYtVRcmSfOtGx5iapHwXjiibrXP3u4hzmW3fKrQNNAz89omP99QWjWStIRd2zadF6hw4Yi63dBt5xDnfniI6kpRSX1h+5ZR9u7YzOjwEAGMDg+xd8fm88K07D77fP3wELWjaYuSVJXlWibNRtRVhm6rLZ9ucoQuqWcUXexT5cyVTh6irjQDXVLPWKpl8qE7v9XSzJIqQ7eVlk+32XKR1DOWao2cy9kdRuY/5ISlZ6hUtT9LOw9Ru8FAl9QzlupTzzdz9hz/9J+P8L9nn15yNksvh26VbLlI6hmLtUwW88snzlY2m6WfOUKX1DMWtkyeE/FMu6UVvTSFsBsMdEk9ZX7LZOFCIZh9yPk7q5/DmZmzF/xuL00h7AYDXVLPWuohJ9B01eggMtAl9bTlHnLW6bShMhjokvrSIM9mWYqzXCSpJhyhS1pRdTuYuZcY6JJWTD/sKd7PbLlIWjFVbm8rA13SCuqHPcX7mYEuacXU8WDmXtI00CPiYEScjoiH513bHxHfjYhvR8TnImK42jIl1UE/7Cnez1oZod8KXLXg2r3AZZl5OfA9YE/JdUmqoX7YU7yftXJI9P0RsWHBtXvmvfwG8I5yy5JUVy4Iqk4ZPfR3A19e6ocRsTMiJiJiYnp6uoTbSZIWUyjQI+IjwFPA7Uu9JzMPZOZYZo6NjIwUuZ0kaRkdLyyKiOuBq4ErM9vYsFiSVImOAj0irgI+DPxZZj5RbkmSpE60Mm3xDuDrwKaIOBkR7wE+CbwQuDciHoqIT1dcpySpiVZmuVy7yOVbKqhFklSAK0UlqSYMdEmqCQNdkmrCQJekmjDQJakmDHRJqgkDXZJqwjNFpQHgwcyDwUCXas6DmQeHLRep5jyYeXAY6FLNeTDz4DDQpZrzYObBYaBLNefBzIPDh6JSzc09+HSWS/0Z6NIA8GDmwWDLRZJqwkCXpJpo5Qi6gxFxOiIennftxRFxb0R8v/H14mrLlCQ108oI/VbgqgXXdgP3ZeYrgPsaryVJXdQ00DPzfuDxBZevAW5rfH8bsL3kuiRJbeq0h/7SzHwMoPH1JeWVJEnqROUPRSNiZ0RMRMTE9PR01beTpIHVaaD/PCJeBtD4enqpN2bmgcwcy8yxkZGRDm8nSWqm04VFXwCuB/Y1vn6+tIqkGnNfclWpaaBHxB3AG4G1EXESuInZIL8zIt4DPAr8VZVFSnXgvuSqWtNAz8xrl/jRlSXXItXCUqPw5fYlN9BVBvdykUq03CjcfclVNZf+SyVabhTuvuSqmoEulWi5Ubj7kqtqBrpUouVG4du3jLJ3x2ZGh4cIYHR4iL07Nts/V2nsoUsl2rVt03k9dDh/FO6+5KqSgS6VyNOB1E0GulQyR+HqFnvoklQTBrok1YSBLkk1YaBLUk0Y6JJUEwa6JNWEgS5JNWGgS1JNGOiSVBMGuiTVRKFAj4gPRsQjEfFwRNwREc8rqzBJUns6DvSIGAXeD4xl5mXAKuCdZRUmSWpP0ZbLamAoIlYDFwGnipckSepEx4GemVPAR4FHgceAX2XmPWUVJklqT5GWy8XANcBGYB3w/Ii4bpH37YyIiYiYmJ6e7rxSSdKyirRc3gz8KDOnM/MscAh4/cI3ZeaBzBzLzLGRkZECt5MkLadIoD8KXBERF0VEAFcCx8spS5LUriI99AeAu4AHgWONzzpQUl2SpDYVOoIuM28CbiqpFklSAa4UlaSaMNAlqSYKtVyklTY+OcX+wyc4dWaGdcND7Nq2ie1bRrtdltQTDHT1jfHJKfYcOsbM2XMATJ2ZYc+hYwCGuoQtF/WR/YdPPBPmc2bOnmP/4RNdqkjqLQa6+sapMzNtXZcGjYGuvrFueKit69KgMdDVN3Zt28TQmlXnXRtas4pd2zZ1qSKpt/hQVJUqc1bK3O85y0VanIGuylQxK2X7llEDXFqCLRdVxlkp0soy0FUZZ6VIK8tAV2WclSKtLANdlal6Vsr45BRb9x1h4+672brvCOOTU6V8rtSvfCiqylQ5K8VtAKQLGeiqVFWzUpZ74Gqga1DZclFf8oGrdCEDXX3JB67ShQoFekQMR8RdEfHdiDgeEa8rqzBpOW4DIF2oaA/9E8BXMvMdEfFc4KISapKachsA6UIdB3pEvAh4A/C3AJn5JPBkOWVJzbkNgHS+IiP0lwPTwGci4tXAUeCGzPzt/DdFxE5gJ8D69esL3E514TFyUjWK9NBXA68BPpWZW4DfArsXvikzD2TmWGaOjYyMFLid6mBu/vjUmRmSZ+ePuyhIKq5IoJ8ETmbmA43XdzEb8NKS3LBLqk7HgZ6ZPwN+GhFz0wquBL5TSlWqLeePS9UpOsvlfcDtjRkuPwT+rnhJqrN1w0NMLRLezh+Xiis0Dz0zH2r0xy/PzO2Z+cuyClM9OX9cqo57uWhFOX9cqo6BrhXn/HGpGu7lIkk1YaBLUk0Y6JJUE/bQ1ZRL9aX+YKAPqFZD2qPepP5hy2UAtbOfikv1pf5hoA+gdkLapfpS/zDQB1A7Ie1Rb1L/MNAHUDsh7VJ9qX8Y6AOonZDevmWUvTs2Mzo8RACjw0Ps3bHZB6JSD3KWywBqdz8Vl+pL/cFAH1CGtFQ/tlwkqSYcoQtwNahUBwa6XA0q1UThlktErIqIyYj4YhkFqXPjk1Ns3XeEjbvvZuu+I4uu/FyMq0GleihjhH4DcBx4UQmfpQ4VGWW7GlSqh0Ij9Ii4BHgbcHM55ahTRUbZrgaV6qFoy+XjwI3A0yXUogKKjLJdDSrVQ8eBHhFXA6cz82iT9+2MiImImJienu70dmqiyCjb1aBSPURmdvaLEXuBdwFPAc9jtod+KDOvW+p3xsbGcmJioqP7aXkLe+gwO8o2mKX+FxFHM3Os2fs6HqFn5p7MvCQzNwDvBI4sF+aqlqNsSc5DrxGX80uDrZRAz8yvAV8r47MkSZ1xLxdJqgkDXZJqwkCXpJow0CWpJgx0SaoJA12SasJAl6SaMNAlqSYMdEmqCQNdkmrCQJekmjDQJakmDHRJqgkDXZJqwkCXpJow0CWpJgx0SaqJjgM9Ii6NiK9GxPGIeCQibiizMElSe4ocQfcU8KHMfDAiXggcjYh7M/M7JdVW2PjkFPsPn+DUmRnWDQ+xa9smz9yUVFsdB3pmPgY81vj+NxFxHBgFeiLQxyen2HPoGDNnzwEwdWaGPYeOARjqkmqplB56RGwAtgAPlPF5Zdh/+MQzYT5n5uw59h8+0aWKJKlahQM9Il4AfBb4QGb+epGf74yIiYiYmJ6eLnq7lp06M9PWdUnqd4UCPSLWMBvmt2fmocXek5kHMnMsM8dGRkaK3K4t64aH2rouSf2uyCyXAG4Bjmfmx8orqRy7tm1iaM2q864NrVnFrm2bulSRJFWryAh9K/Au4E0R8VDjn7eWVFdh27eMsnfHZkaHhwhgdHiIvTs2+0BUUm0VmeXy30CUWEvptm8ZNcAlDYwi89D7knPTJdXVQAW6c9Ml1dlA7eXi3HRJdTYwgT4+OcWUc9Ml1dhABPpcq2Upzk2XVAcDEeiLtVrmODddUl0MRKAv11JxbrqkuhiIQF+qpTI6PGSYS6qNvg308ckptu47wsbdd7N13xHGJ6eWfK/bAEgaBH05D73d+eRz11xQJKnO+jLQl5tPvlRIuw2ApLrr+UBfbKm+e51L0oV6OtCXaq0MX7SGXz5x9oL3O59c0iDr6YeiS7VWMvEhpyQt0NOBvlQL5VczZ93rXJIW6OmWy7rhoUX3X1nXmD9ugEvSs3p6hO78cUlqXU+P0J0/LkmtKxToEXEV8AlgFXBzZu4rpap5bK1IUms6brlExCrgX4G/BF4FXBsRryqrMElSe4r00P8E+EFm/jAznwT+A7imnLIkSe0qEuijwE/nvT7ZuHaeiNgZERMRMTE9PV3gdpKk5RQJ9FjkWl5wIfNAZo5l5tjIyEiB20mSllMk0E8Cl857fQlwqlg5kqROReYFg+rWfjFiNfA94EpgCvgm8NeZ+cgyvzMN/KSjG66stcAvul1Eh6x95fVr3WDt3dBJ3b+XmU1bHB1PW8zMpyLivcBhZqctHlwuzBu/0xc9l4iYyMyxbtfRCWtfef1aN1h7N1RZd6F56Jn5JeBLJdUiSSqgp5f+S5JaZ6Av7kC3CyjA2ldev9YN1t4NldXd8UNRSVJvcYQuSTVhoC8QEVdFxImI+EFE7O52Pa2IiEsj4qsRcTwiHomIG7pdU7siYlVETEbEF7tdSzsiYjgi7oqI7zb++7+u2zW1KiI+2Pjz8nBE3BERz+t2TUuJiIMRcToiHp537cURcW9EfL/x9eJu1riYJere3/jz8u2I+FxEDJd1PwN9nj7ecOwp4EOZ+YfAFcA/9End890AHO92ER34BPCVzPwD4NX0yb9DRIwC7wfGMvMyZqcev7O7VS3rVuCqBdd2A/dl5iuA+xqve82tXFj3vcBlmXk5s2t59pR1MwP9fH254VhmPpaZDza+/w2zodI3ew5HxCXA24Cbu11LOyLiRcAbgFsAMvPJzDzT3arashoYaiwSvIgeXumdmfcDjy+4fA1wW+P724DtK1pUCxarOzPvycynGi+/wewq+1IY6OdracOxXhYRG4AtwAPdraQtHwduBJ7udiFtejkwDXym0S66OSKe3+2iWpGZU8BHgUeBx4BfZeY93a2qbS/NzMdgdlADvKTL9XTi3cCXy/owA/18LW041qsi4gXAZ4EPZOavu11PKyLiauB0Zh7tdi0dWA28BvhUZm4Bfktv/rX/Ao1+8zXARmAd8PyIuK67VQ2WiPgIs+3S28v6TAP9fH274VhErGE2zG/PzEPdrqcNW4G3R8SPmW1xvSki/r27JbXsJHAyM+f+NnQXswHfD94M/CgzpzPzLHAIeH2Xa2rXzyPiZQCNr6e7XE/LIuJ64Grgb7LEueMG+vm+CbwiIjZGxHOZfUj0hS7X1FREBLN93OOZ+bFu19OOzNyTmZdk5gZm/3sfycy+GClm5s+An0bE3KnlVwLf6WJJ7XgUuCIiLmr8+bmSPnmgO88XgOsb318PfL6LtbSscXTnh4G3Z+YTZX62gT5P40HF3IZjx4E7m2041iO2Au9idnT7UOOft3a7qAHxPuD2iPg28MfAP3e5npY0/lZxF/AgcIzZLOjZlZcRcQfwdWBTRJyMiPcA+4C3RMT3gbc0XveUJer+JPBC4N7G/6ufLu1+rhSVpHpwhC5JNWGgS1JNGOiSVBMGuiTVhIEuSTVhoEtSTRjoklQTBrok1cT/AwD7AA9HFqPmAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(x, y)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Ordinary Least Squares"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAHCRJREFUeJzt3Xt4XWWZ9/HvTVpoUNqADUgCIYXBytliFLCoWNAW6AuFd1A5iYDUQQYQoUBAIcBocdqBwbc4vBVqUWnlFFKcMhQGRBSRoSVIkVIdKEObAi2WlINJmyb3/LF3kr32zmEf1s7ah9/nuria/ey197rhCr/r7rOeZy1zd0REpPhtF3UBIiISDgW6iEiJUKCLiJQIBbqISIlQoIuIlAgFuohIiVCgi4iUCAW6iEiJUKCLiJSIUSN5svHjx3t9ff1InlJEpOitWLHibXevHu64YQPdzBYA04EN7n5gwviFwD8C24Cl7n75cN9VX1/P8uXLhztMREQSmNn/pHNcOlMuC4FpSV/+BeBE4GB3PwCYm2mBIiISrmED3d2fBDYlDZ8P3OjuW+LHbMhDbSIikoFsL4p+DPismT1jZr8xs08NdqCZzTSz5Wa2fOPGjVmeTkREhpNtoI8CdgYOB2YB95iZDXSgu8939wZ3b6iuHnZOX0REspRtoK8Dmj3mv4AeYHx4ZYmISKayXbbYAkwBnjCzjwHbA2+HVpWISIloaW1jzrLVrG/voKaqkllTJzJjUm1ezpXOssXFwFHAeDNbB1wLLAAWmNmLwFbgLNejj0REAlpa22hsXklHVzcAbe0dNDavBMhLqA8b6O5+6iBvnRFyLSIiJWXOstV9Yd6ro6ubOctW5yXQtfVfRCRP1rd3ZDSeKwW6iEie1FRVZjSeKwW6iEiezJo6kcrR27Er7/SNVY6uYNbUiXk534jenEtEpJzMqH6DGRVfhQo4astNdI2bEO0qFxERyVD3NrjtSNi4KvZ6XB1PfPscGHj/ZWgU6CIiYfpTC9x7Vv/rry2BvY8akVMr0EVEwrDlPZi9R//rfabAGc1578oTKdBFRHL11C3w6DX9r7/1DOz68REvQ4EuIpKtzW1w8/79rw87H469MbJyFOgiItlYcgG0/qL/9WV/gQ/vGl09KNBFRDLTehcs+Vb/6+PmwqfPi66eBAp0EZF0dG+DGz4SHLvqDdh+x2jqGYB2ioqIDGfpZcEwbzgXmjYXVJiDOnQRkcH9bRP884Tg2Pf+ChWFGZ2FWZWISNRuPQw2vtz/+sRbYVJh3zVcgS4ikujNF+G2ycGxps3R1JIhBbqISK+mccHX5/4n7PmpaGrJwrAXRc1sgZltiD9uLvm9y8zMzUwPiBaR4vXSkmCY7/iRWFdeRGEO6XXoC4F5wM8SB81sT+CLwOvhlyUikn8tz73OjAcPCg5+52UYu3s0BeVo2A7d3Z8ENg3w1s3A5YAeDi0iRWfTLZ8NhPmvug9nv+67aXmlJ8KqcpPVHLqZnQC0ufsfbQTvJCYi5a2ltY05y1azvr2DmqrK7B4W8cHbMGcfdkkY2q9zAR2MAfL3AOeRkHGgm9mOwNXAl9I8fiYwE6Curi7T04mIALEwb2xeSUdXNwBt7R00Nq8ESD+Aky56Pth9BBd1XRgYy9cDnEdCNjtF9wEmAH80s9eAPYDnzOyjAx3s7vPdvcHdG6qrq7OvVETK2pxlq/vCvFdHV6yjHta65SlhPnmH5pQwh/w9wHkkZNyhu/tKoO+WYvFQb3D3t0OsS0QkYLDOOXk8eVrmqc6Tgh/4P7fAJ7/OrKSOH4Z/gHMoUz55NGygm9li4ChgvJmtA6519zvyXZiISKKaqkraBgj1xI46cVrmzIpHuKFzYfDghA1CvUGcbkCHMuWTZ8MGurufOsz79aFVIyIyiFlTJw7bUcemZbbx2pjTA589Z/u5LLgq9Ra3MybVph3GQ035FE2gi4gUgnQ66iUdZzF+zLuBz9V3LsI6cz9/ulM+UVKgi0hBGWqeetCOOn5XxPEJq6g/1fljNlIFhHOhM50pn6jpfugiUjB656nb2jtw+uepW1rbBv9Q07iUW9zWdy7qC/PhLnSma9bUiVSOrgiMhfXdYVGgi0jBGGye+tJ7/siEK5cy+cbH+8N97X+l3kzrmk20nPgStVWVGFBbVcnskw8KZY57xqRaZp98UF6+OyzmPnI79xsaGnz58uUjdj4RKS4Trlw67L1EKkdXsKriK4GxpRVT6Jo+r6DCNUxmtsLdG4Y7TnPoIlIwBpun7nV+xYNcUfHLwFh95yIAKgtsCWEUNOUiIgVjoHnqXq+NOY0rRveH+QVbL+oLc8hg12gJU4cuIgUjeWnidma8skPqVpjEIE9USEsIo6BAF5GC0rc0cct7MHuPwHuf23IzG0fVUFW5He0dXSmfLaQlhFFQoItI4UlevQJM6FxETVUls+PLBDO9D0s5UKCLSOFY81u4c3pw7Oq3YPQY1gxweCHfKCsKCnQRKQwDdOWJN9NKlsl9WMqFAl1EorX0Mnj2J8GxIYJcBqdAF5ERlXivljVjTgu+OekMOPHWaAorAQp0ERkxvfdqWVXxFRiT9Ka68pxpY5GIjJhbHn4hZdv+17ZeweQxD0RUUWlRhy4iI6NpHL9OGurdIGRlviEoLMN26Ga2wMw2mNmLCWNzzOxlM3vBzB4ws6r8likiRWvdipQVLId23hbY7VnuG4LCks6Uy0JgWtLYo8CB7n4w8GegMeS6RKQUNI2D26cEhvbrvptNjO17rQ1B4Rk20N39SWBT0tgj7r4t/vIPwB4pHxSR8vXId1PXlV/bDk2bC/6e4sUsjDn0c4C7Q/geESkFyUG+8wS4+Pm+l9oQlD85BbqZXQ1sA+4a4piZwEyAurq6XE4nIoUsw52eEr6sly2a2VnAdOB0H+KxR+4+390b3L2huro629OJSKHq6U4N82PnKMwjkFWHbmbTgCuAz7v738ItSUSKhrrygjJsoJvZYuAoYLyZrQOuJbaqZQfgUTMD+IO7/0Me6xSRQvLXV+D/HRocu6gVdtk7mnoESCPQ3T31cSFwRx5qEZFioK68YGmnqIik5w+3wcNXBMeu2QTbDfwMUBl5CnQRGZ668qKgQBeRwd18IGxeGxxTkBcs3W1RRFK5x7ryxDCffLHCvMCpQxeRIE2vFC0FuojEvL8B5u4bHDvvcaj9ZDT1SMYU6CKirrxEKNBFytmfHoB7vx4cu/otGJ38fDgpBgp0kTKQ+GDmmqpKZk2dyIwl+6ceqK68qCnQRUpc74OZO7q6AZj9wTV8bsnK4EEK8pKgQBcpcXOWre4L89fGnBZ884CT4JSFI1+U5IUCXaTErW/vSA1yYELnItaccnwEFUm+KNBFStmW91mTFObnbf0Oj/Y0UKsHM5ccBbpIqRpgKWJ95yJAD2YuVdr6L1JqXn0iJcyXHvt7Jo95QA9mLnHq0EVKySAbhI4Hjj/sgBEvR0aWAl2kFDTPhBfuDo5pKWLZUaCLFLvkrnyXvWOPg5Oyk84zRRcA04EN7n5gfGwX4G6gHngN+LK7v5O/MkUkhe6/IknSuSi6EJiWNHYl8Ji77ws8Fn8tIiOhuys1zI//F4W5pPWQ6CfNrD5p+ETgqPjPdwJPAEkPGxSR0KkrlyFkO4e+m7u/AeDub5jZriHWJCLJ3lwJtx0ZHLv4Bdh5r2jqkYKU94uiZjYTmAlQV1eX79OJlB515ZKmbDcWvWVmuwPE/9ww2IHuPt/dG9y9obq6OsvTiZShX/8gNcyveUdhLoPKtkN/EDgLuDH+55LQKhIpYQPel3ygHZvqyiUL6SxbXEzsAuh4M1sHXEssyO8xs3OB14FT8lmkSClIvi95W3sHjc2x+5L3hbqCXHKQziqXUwd56+iQaxEpCYN14Yn3Je/V0dXNnGWrmfGJGriuKvhFR14CxzSNWN1S/LRTVCREQ3Xh69s7BvzMU50nwXVJg+rKJQu626JIiIbqwmuS7j++O39NffDEN59UmEvWFOgiIRqsC1/f3sGsqROpHF0BxB4F9/SYC4MHNW2G3Q/Jd4lSwjTlIhKimqpK2gYI9ZqqSmZMqqV23UN8asVlwTe/uxFGbT9CFUopU4cuEqLELrxX39OBmsalhnnTZoW5hEYdukiIepcfJq5yad5pLrsteSp4oObJJQ8U6CIhmzGpNriuvDPhzcP+AY79YSR1SelToIvkgzYISQQU6CJh6ngHflgfHDv7YdjriEjKkfKiQBcJi7pyiZgCXSRXqx+GxV8Jjl25FsaMjaYeKVsKdJFcqCuXAqJAF8nG4tNg9dLgmIJcIqZAF8lUcle++yGxe7CIREyBLpIuTa9IgVOgiwynqxO+v1twbMZt8InBHhUgEg0FushQ1JVLEckp0M3sEuAbgAMrgbPdvXPoT4kUgbYV8JMpwbHvvAxjd4+mHpE0ZB3oZlYLXATs7+4dZnYP8FVgYUi1iURDXbkUqVynXEYBlWbWBewIrM+9JJGILLsanp4XHLu2HcyiqUckQ1kHuru3mdlc4HWgA3jE3R8JrTKRkZTclY+qhO++GU0tIlnKZcplZ+BEYALQDtxrZme4+y+SjpsJzASoq6vLoVSRPND0ipSQXJ5YdAywxt03unsX0Ax8Jvkgd5/v7g3u3lBdXZ3D6URC1NOTGuZTvqcwl6KWyxz668DhZrYjsSmXo4HloVQlkk/qyqVE5TKH/oyZ3Qc8B2wDWoH5YRUmErpNr8KPJgXHvvUM7PrxaOoRCVlOq1zc/Vrg2pBqEckfdeVSBrRTVErb0z+GZY3Bse/9FSr0qy+lR7/VUlRaWtuYs2w169s7qKmqZNbUif0PZE6mrlzKjAJdikZLaxuNzSvp6OoGoK29g8bmlQDBUJ87Ed5PWkOuIJcykMuyRZERNWfZ6r4w79XR1c2cZav7B5rGBcP80K8pzKVsqEOXorG+vWPwcU2viKhDl+JRU1WZMvYRNrNmzGnBwa8/pDCXsqQOXYrGrKkTA3PoryUHOSjIpawp0KVo9F74fPahn/L9rjnBN69aD9t/KIKqRAqHAl3yKqNlhmmYsWR/ZiQPqisXARTokkdpLzNMx/3fgJX3BscU5CIBuigqeZPWMsN0NI0LhvmBf68wFxmAOnTJmyGXGaZDSxFFMqIOXfJmoGWGQ4336RpgXfnp9ynMRYahDl3yJnmZIUDl6ApmTZ04+Icy6MrDvuAqUuwU6JI3veGaVui2rYCfTAmOXb4GdtxlwO8O9YKrSIlQoEtezZhUO3zAZjFXPtQFVwW6lCsFukTnP5vgdzcHx9KcJ8/5gqtICVKgSzSSu/LaBjjvsbQ/XlNVSdsA4T3sBVeREpbTKhczqzKz+8zsZTNbZWZHhFWYlKh5n04N86bNGYU5xC64Vo6uCIwNe8FVpMTl2qHfAjzs7n9vZtsDO4ZQk5Sinm64PukC55d/DvufkNXXZXTBVaRMZB3oZjYW+BzwdQB33wpsDacsKSl52iCU1gVXkTKSS4e+N7AR+KmZHQKsAC529w8SDzKzmcBMgLq6uhxOJ0XnvTfhX5KmQL7zMi2v9DDnxsfVWYuELJc59FHAocC/ufsk4APgyuSD3H2+uze4e0N1dXUOp5Oi0jQuNcybNtPySg+NzStpa+/A6V8/3tLaFkmZIqUkl0BfB6xz92fir+8jFvBSzv78SOoUyzXv9E2xhHbDLhFJkfWUi7u/aWZrzWyiu68GjgZeCq80KTrJQX7Ql+H//iQwpPXjIvmT6yqXC4G74itcXgXOzr0kKTpLL4Nng8E92EVPrR8XyZ+cAt3dnwcaQqpFio07XFcVHDtlIRxw0qAfyeqGXSKSFu0UlezcsCt0bwmOpbEUUevHRfJHgS6Z6WiHH+4VHLv4Bdh5r4GPH4DWj4vkhwJd0qcnCIkUNAW6DG/ts3DHMcGx770NFaOjqUdEBqRAl6Eld+V7HQlnL42mFhEZkp4pKgP7zZyUMJ/QuYjJb12qXZ0iBUodepka8nmcSUF+fc85LNgam3LRo95ECpcCvQwN9jzOox+bzk7vvxo4dvKYB1I2AulRbyKFSVMuZSj5fio7sJVVFV8Jhvn5v4emzdqqL1JE1KGXocQwfm3MaakHJCxF1FZ9keKhDr0M1VRV8ne2LiXMp+zwy5R15XrUm0jxUIdehp7qPAl26H+9tqeaL/k8Zk87KOVYbdUXKR4K9HKyYiH86uLA0ITORdRUVTJ7iJDWVn2R4qBALxfJG4Q+fyV8oZE10VQjInmgQC91Pz8JXnk8OKb7r4iUJAV6qereBjd8JDh29sOw1xEDHj7kRiMRKQoK9FKU4V0RB9toBNoNKlJMcl62aGYVZtZqZv8eRkGSvWW/X54a5pevGXaKRQ9uFikNYaxDvxhYFcL3SC6axjH1kaMDQ/t1303L6uF3dGo3qEhpyCnQzWwP4Hjg9nDKkYyt+lVKV17feRf1nYvS7rIH2/Wp3aAixSXXDv1fgcuBnhBqkUw1jYO7z+h7uWjbFOo7FwHWN5ZOl63doCKlIeuLomY2Hdjg7ivM7KghjpsJzASoq6vL9nSSaMkF0PqLwNBAd0WE9Lps7QYVKQ3m7tl90Gw2cCawDRgDjAWa3f2MwT7T0NDgy5cvz+p8ArjDdVXBsa8uho8fl7JSBWJd9uyTD1IwixQ5M1vh7g3DHZd1h+7ujUBj/GRHAZcNFeaSo2GWIqrLFhGtQy90f9sE/zwhOHbJSzAuNah1zxWR8hZKoLv7E8ATYXyXJMhwg5CIlDd16IXotd/BwuODY9dsgu0qBj5eRAQFeuFJ7sr3/RKcfm80tYhIUVGgF4rHboDfzg2OaXpFRDKgQC8EyV35CfPg0DOjqUVEipYCPUo3HQDvrguOqSsXkSwp0KOw9QP4QU1w7IJnofpj0dQjIiVBgT7StBRRRPJEgT5S3lwJtx0ZHLv6LRg9Jpp6RKTkKNBHQnJXvusB8K3fR1OLiJQsBXo+PfP/4T8uD45pekVE8kSBni/JXfkxTXDkJVFUIiJlQoEetgXT4PWng2PqykVkBCjQw9LdBTeMD45943HY45PR1CMiZUeBHgYtRRSRAqBAz8WmNfCjTwTHrlwLY8ZGU4+IlDUFeraSu/IdxkLj2mhqERFBgZ65lffB/ecGx65tB7No6hERics60M1sT+BnwEeBHmC+u98SVmEFKbkrP+x8OPbGaGoREUmSS4e+DbjU3Z8zs52AFWb2qLu/FFJtOWtpbQvnocn3nQMv3h8c00VPESkwWQe6u78BvBH/+T0zWwXUAgUR6C2tbTQ2r6SjqxuAtvYOGptXAqQf6j09cP3OwbEz7oe/OybMUkVEQhHKHLqZ1QOTgGfC+L4wzFm2ui/Me3V0dTNn2er0Al1LEUWkyOQc6Gb2YeB+4Nvu/u4A788EZgLU1dXlerq0rW/vyGi8z/sbYO6+wbFL/ww77RZSZSIi+ZFToJvZaGJhfpe7Nw90jLvPB+YDNDQ0eC7ny0RNVSVtA4R3TVXl4B9SVy4iRWy7bD9oZgbcAaxy95vCKykcs6ZOpHJ0RWCscnQFs6ZOTD34lcdTw/yadxTmIlJUcunQJwNnAivN7Pn42FXu/lDuZeWud5582FUuyUG+/wz48p0jVKWISHhyWeXyO6Cgd9PMmFQ7+AXQZVfD0/OCY+rIRaSIld1O0ZbWNmYs2T84ePLtcPAp0RQkIhKSrOfQi9Fbtx6XEub7dd9NS/dnIqpIRCQ85dGhd3XA9z9K4sLDwzrn8Ra7ABmsTRcRKWClH+iz94Qt/cvj3/NKDtpyR+CQYdemi4gUgdIN9AE2CO3T+XO6qUg5dMi16SIiRaI0A/2mA+DddX0vf1pxCtd9cNKAhw66Nl1EpMiU1kXRtudi68oTwpymzVw/SJgDzD75IM2fi0hJKJ0OPXmD0DefhN0PAQa/DUBtVaXCXERKRtF26C2tbUy+8XEuuaoxGOY7T4htEIqHOWR4GwARkSJVlB16S2sbVzX/kZcqToXt+8cfmvZbjjv84JTj074NgIhIESvKQP/NQ7/kpYrr+14v3vYFGredR+0Tb3Pc4QN/ZsjbAIiIlICCD/TEx8jVj9uOZT3f5Oau2D1XNvpYjtgyj23xfw2tJxeRclbQgZ74GLlTKx5jdsKGoOlb/okXfe/A8VpPLiLlrKADvfcxcjeN/jEnV/wOgPu7P8sNoy5iy6geSHjEnC5yiki5K+hVLr1TKM/37EO3G5/p/BGXdp3P5o4uZp98ELVVlRix5YdaTy4i5a6gO/Te9eM/657Kz7qnBsZ1kVNEJKigO3StHxcRSV9Bd+haPy4ikr6cAt3MpgG3ABXA7e5+YyhVJdDUiohIerKecjGzCuBW4Fhgf+BUM9t/6E+JiEi+5DKH/mngv939VXffCvwSODGcskREJFO5BHotsDbh9br4WICZzTSz5Wa2fOPGjTmcTkREhpJLoNsAY54y4D7f3RvcvaG6ujqH04mIyFByCfR1wJ4Jr/cA1udWjoiIZMvcU5rq9D5oNgr4M3A00AY8C5zm7n8a4jMbgf/J6oQjazzwdtRFZEm1j7xirRtUe1QyrX0vdx92iiPrZYvuvs3M/hFYRmzZ4oKhwjz+maKYczGz5e7eEHUd2VDtI69Y6wbVHpV81Z7TOnR3fwh4KKRaREQkBwW99V9ERNKnQB/Y/KgLyIFqH3nFWjeo9qjkpfasL4qKiEhhUYcuIlIiFOgJzGyama02s/82syujriddZranmf3azFaZ2Z/M7OKoa8qUmVWYWauZ/XvUtWTCzKrM7D4zezn+3/+IqGtKl5ldEv99edHMFpvZmKhrGoyZLTCzDWb2YsLYLmb2qJn9Jf7nzlHWOJBB6p4T/315wcweMLOqsM6nQI8r8puNbQMudff9gMOBC4qo9l4XA6uiLiILtwAPu/vHgUMokn8HM6sFLgIa3P1AYkuPvxptVUNaCExLGrsSeMzd9wUei78uNAtJrftR4EB3P5jYXp7GsE6mQO9XtDcbc/c33P25+M/vEQuVornnsJntARwP3B51LZkws7HA54A7ANx9q7u3R1tVRkYBlfFNgjtSwDu93f1JYFPS8InAnfGf7wRmjGhRaRiobnd/xN23xV/+gdgu+1Ao0PuldbOxQmdm9cAk4JloK8nIvwKXAz1RF5KhvYGNwE/j00W3m9mHoi4qHe7eBswFXgfeADa7+yPRVpWx3dz9DYg1NcCuEdeTjXOA/wjryxTo/dK62VghM7MPA/cD33b3d6OuJx1mNh3Y4O4roq4lC6OAQ4F/c/dJwAcU5l/7U8Tnm08EJgA1wIfM7IxoqyovZnY1senSu8L6TgV6v6K+2ZiZjSYW5ne5e3PU9WRgMnCCmb1GbJpripn9ItqS0rYOWOfuvX8buo9YwBeDY4A17r7R3buAZuAzEdeUqbfMbHeA+J8bIq4nbWZ2FjAdON1DXDuuQO/3LLCvmU0ws+2JXSB6MOKa0mJmRmwed5W73xR1PZlw90Z338Pd64n9N3/c3YuiU3T3N4G1Ztb71PKjgZciLCkTrwOHm9mO8d+foymSC7oJHgTOiv98FrAkwlrSFn905xXACe7+tzC/W4EeF79I0XuzsVXAPcPdbKyATAbOJNbdPh//57ioiyoTFwJ3mdkLwCeAH0RcT1rif6u4D3gOWEksCwp256WZLQaeBiaa2TozOxe4Efiimf0F+GL8dUEZpO55wE7Ao/H/V28L7XzaKSoiUhrUoYuIlAgFuohIiVCgi4iUCAW6iEiJUKCLiJQIBbqISIlQoIuIlAgFuohIifhfsabYVJqo8lkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "X = np.vander(x, 2)\n",
    "w = la.lstsq(X, y)[0]\n",
    "ŷ = X.dot(w)\n",
    "\n",
    "plt.plot(x, y, 'o')\n",
    "plt.plot(x, ŷ, '-')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Comparison of Implementation Variants"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Computing a data matrix for linear regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def data_matrix_V1(x):\n",
    "    n = len(x)\n",
    "    return np.vstack((x, np.ones(n))).T\n",
    "\n",
    "def data_matrix_V2(x):\n",
    "    return np.vstack((x, np.ones_like(x))).T\n",
    "    \n",
    "def data_matrix_V3(x):\n",
    "    return np.vander(x, 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### All three variants create the same matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[4.7 1. ]\n",
      " [1.1 1. ]]\n",
      "\n",
      "[[4.7 1. ]\n",
      " [1.1 1. ]]\n",
      "\n",
      "[[4.7 1. ]\n",
      " [1.1 1. ]]\n"
     ]
    }
   ],
   "source": [
    "x = [4.7, 1.1]\n",
    "\n",
    "print(data_matrix_V1(x)); print()\n",
    "print(data_matrix_V2(x)); print()\n",
    "print(data_matrix_V3(x))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Solving Least Squares for Linear Regression\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def lsq_solution_V0(X, y):\n",
    "    w = (X.T * X).I * X.T * y\n",
    "    return w\n",
    "\n",
    "def lsq_solution_V1(X, y):\n",
    "    w = la.inv(X.T.dot(X)).dot(X.T).dot(y)\n",
    "    return w\n",
    "\n",
    "def lsq_solution_V2(X, y):\n",
    "    w = np.dot(la.pinv(X), y)\n",
    "    return w\n",
    "\n",
    "def lsq_solution_V3(X, y):\n",
    "    w, residual, rank, svalues = la.lstsq(X, y)\n",
    "    return w"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### All three variants create the same regression line"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABE8AAAFACAYAAABNzFJlAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3Xl8XFX5x/Hv0zRtwxqWIrRSWrayFGyhqIDsYFGQRlBQkEVBcGPVCAVsWYRWqgI/UbAglB1BSkC2AiI7ioUghUJlh6YsZUmhNGnT9Pz+mMncM5OZyWTmZm7uzOf9es2Le+7czDxJ2y/PnNx7rjnnBAAAAAAAgOwGRF0AAAAAAABAf8bkCQAAAAAAQB5MngAAAAAAAOTB5AkAAAAAAEAeTJ4AAAAAAADkweQJAAAAAABAHkyeoChm9oKZ7R51HQDijzwBECYyBUCYyBR0YfKkQpjZbDM7J8v+iWb2rpkNDPP9nHNbO+ceCvM1JcnMBpvZlWb2SbLuU8J+DwD5VVCeHGxmT5jZUjML/fUBFKaCMuW3ZvaymX1qZi+Z2RFhvweAnlVQplxgZm8nP/e8aWZnhP0eCBeTJ5VjpqTDzcwy9h8u6Xrn3IrevFjYodMLZ0naTNJGkvaQ9Esz2zeiWoBqNVOVkScfSbpI0rSI3h9AwkxVRqZ8JukbktaUdKSki81sp4hqAarZTFVGpvxF0hbOuTUk7STpUDM7MKJaUAAmTypHk6S1Je3StcPM1pK0v6RrkuPByd+avGVm75nZZWZWl3xudzNbYGanmtm7kq4ys3XN7E4zazWzj8zsUTMbkDz+DTPb23vdi8xsYfJxkZkNznjdn5vZ+2b2jpl9P8/3cYSkc51zHzvnXpR0uaSjwv5hAcirIvLEOfeAc+5mSQv75scEoECVkilTnHMvOedWOuf+LelRSTv2yU8MQD6VkinznXOfebtWSto01J8UQsXkSYVwzrVJulmJyYcuB0t6yTn33+T4N5I2lzRWiX+YwyVN9o5fX4kg2kjSsZJ+LmmBpKGSPifpdEkuy9ufIenLydf9gqQvSjoz43XXTL7f0ZL+mAy4NMl9wyT919v9X0lb5/3mAYSqEvIEQP9RiZmS/BC2g6QXejoWQLgqKVPM7DQzW5J871Ul3ZD/u0eUmDypLFdL+nbXrKoSgXK1JJmZSfqhpJOdcx855z6VdL6k73hfv1LSFOfcsmQodUjaQNJGzrkO59yjzrlsIXKYpHOcc+875xZJOluJ0+a6dCSf73DO3S1piaTRWV5nteR/F3v7FktavdAfAIDQxD1PAPQvlZYplynxC57ZBX33AMJWEZninJumxGed7SRdq/TPQehnmDypIM65xyQtkjTRzDZW4jciXbOXQyWtIunp5OlorZLuTe7vssg51+6Np0t6RdJ9ZvaamZ2W462HSXrTG7+Z3Nflw4xrD5cqmCjxLUn+dw1v3xqSPs3xvgD6SAXkCYB+pJIyxcymSxoj6eAcH64A9LFKyhSX0CypTYnJGPRTTJ5UnmuUmHk9XNJ9zrn3kvs/UOIf5NbOufrkY03nnP+POa0BcM596pz7uXNuYyUWSDvFzPbK8p4LlTjlrcsIFbHGgHPuY0nvKHEKXJcviFNigajENk8A9EuxzxQzO1vS1yR91Tn3SbGvAyAUsc+UDAMlbRLSa6EPMHlSea6RtLcSp6pd3bXTObdSicVXLzSz9STJzIab2YRcL2Rm+5vZpslT3z6R1Jl8ZLpR0plmNtTM1lXiesLrSqj/TDNby8y2SH4fM4t8LQCliXWemFmNmQ1RohkZYGZDzKy2mNcCEIq4Z8okSYdK2sc592ExrwEgVLHNFDMbYGbHJT/zmJl9UdJPJf2jt6+F8mHypMI4596Q9IQSCw7dkfH0qUqcjvYvM/tE0gPKf13vZsljlkh6UtKfXPZ7nP9a0hxJz0maK+mZ5L5iTJH0qhKnwD0sabpz7t4iXwtACSogTw5X4jdPlyqxIn+bEs0UgAhUQKacr8RvmV82syXJx+lFvhaAElVApnxTic89nyoxAfOH5AP9lHGpJgAAAAAAQG6ceQIAAAAAAJAHkycAAAAAAAB5MHkCAAAAAACQB5MnAAAAAAAAeTB5AgAAAAAAkAeTJ71kZrub2crk7en2jboeVB4z2zz596vTzI6Juh70HfIEfY08qS5kCvoamVI9yBP0tTjmCZMnxVnonFvNOXevJJnZBmZ2h5ktNDNnZiPzfbGZjTSzf5rZUjN7ycz29p4bbGYXJl/rYzP7k5nVes9vaWYPmtliM3vFzL6Z8doHm9mLZvapmc0zs4YcNTyYrHVgcryemd2YfN/FZva4mX3JO97M7Awze8vMPjGzm8xsDe/54WZ2u5l9ZGYLzOxHeb7/L5vZ/cljF5nZLWa2QZ7j1zaz28zsMzN708wOzXHcVcnvaVNv35KMR6eZdbt/uplNSX6t/2dxsJk9kfxzeijj+F2yvLYzs4O8Y042s3eTP88rzWyw99xYM3s0+dwCM5vc9Zxz7n/OudUkPZrrZ4KKQp6QJ+QJwkSmkClkCsJCnpAn5InPOcejFw9Ju0takLHvc5J+ImlHSU7SyB5e40lJv5dUJ+kgSa2Shiafm6LEX6C1JQ2V9C9JZyefGyjpf5JOkVQjaU9Jn0naPPn8cEnLJX1NkknaT9JSSetlvP9hkh5J1jowuW/j5OtukHztYyV9IGm15PNHSnpJ0oaSVpN0u6Srvdf8p6SLJNVK+oKkjyTtkeP7/5qkb0taQ9Iqkq6UdG+en9eNkv6afN+vSFosaeuMY77ifU+b5nidVSUtkbRrxv5NJM2VtFDS3t7+vSUdLGmypIcK+HvxqaRVk+MJkt6TtLWktSQ9JGmad/w8Seclf9abSHpH0gEZr/mQpGOi/jvPo+8e5Al5kufvBXnCo9cPMoVMyfP3gkzh0asHeUKe5Pl7UbV5EnkBoX9D0vcl/d0bvyLpZm/8tqSxJbz+7soIEu+5geohSCRtLmmZpNW9fY9K+lFye46kb3vPHSrp7eT2mOQ/BPOev0/SucntL0l6P+P9Fkna0RuvqUQYfdkPkhy1fiJp++T23yQ1es/tJKk9GQSrJV9rqPf8DEnXFvgz3U7SpzmeW1WJcNzc23dtxj/KgZKaJW3bQ5AcKek1/+eX3H+PpK9LesMPEu/5YwoIkqskXeWNb5B0vjfeS9K73nippK288S2SJmW8ZmyCpFIf5Al5Qp7wCPNBppApZAqPsB7kCXlCnpT/UYmX7TwsaRczG5A8LapW0s6SZGYbK/GX/rlsX2hmrXkep4VU39aSXnPOfert+29yv5SYPTW/LEmfN7M1M/b7z49Jbs+R9KKZHWBmNcnT15Yp/fs9X9Klkt7NV6SZjZU0SIkgzlXXYEmbefsznx+jwuwq6YUcz20uqdM59z9vn//zkqSTJT3inMv65+o5UtI1LvmvVJLM7NuSljvn7i6w1m7MbBVJ35J0tbd762Sdfs2fM7N1kuOLJB1hZrVmNlqJGfwHiq0BfYY8IU9yIU9QDDKFTMmFTEFvkSfkSS7kSR+puMkT59xrSpxKNFbSbpJmS2oxsy2S40edcytzfG19nse0kEpcTYlTsHyLJa2e3L5H0olmNtTM1pd0QnL/KkqcQva+pMbkX8CvJr+nVZL1d0q6RokZwGXJ/x7nnPtMksxsvBKh2u36N58lruu7VolT57pqvUfSMZa4dnFNSad21ZUMxccl/crMhpjZdkqcmrdKTz8MM9tWiVPEGnMckvfnZWYbSjou+Rr53meEEj+rq719qykRrCf1VGcPDlLidL+H89Tdtd3153ynEuHTpsSf61+cc/8psQ6EjDwhT3K8D3mCopApZEqO9yFT0GvkCXmS433Ikz5UcZMnSQ8rcarZrsnth5T4S7Sb0v+wo7BEievefGsoEX5S4pqwZknPSnpCUpOkDiVOTeuQ1KDEdX3vSvq5pJslLZAkSyz8c4ES3/sgJb7fKyyxUM8ASX+SdKJzbkWu4sysTtLfJf3LOTfVe+pKJa7De0iJGdN/JvcvSP73MEmjlDhF8FJJ13vP5XqvTZUMTudcroWCevp5XSTpHC/wcjlC0mPOude9fWcrcZrd6zm+plDdZnfVve6u7U/NbG1J90o6R9IQJa6pnGBmPymxDvQN8oQ8yUSeoBRkCpmSiUxBscgT8iQTedKXerquJ44PST+UdIcSC+J8Xol/eDdIel3S+DxftyTP4/TkMbur9Ov/2pV+/d8jSl7/l+X4YyU9mef1nlBiplWSfiHptoznm5L76yWtVCKA3lXiukCX3N4leexgJWatb5A0oIef8VeVCIqsxyVfY2qer99Iievtsn7f3nFd1/9t5u27Rsnr/5RYeOo97/tyye/t0IzX+Z+kH2Tse1aJ2dOur+1UYtGnUzOOy3n9nxIhsELSJlm+//O88Z5KXv8nabykjzOOP0nSnRn7HlJMrv+r5Ad5kvY8eeLIEx6lPciUtOfJFEem8Cj+QZ6kPU+eOPKkrx+RF9An31TiH+unkl5JjtdI/gX5RFJNia+dNUiUmE1bNfkXebSkIXle41+Sfpv8mm8qfeXp4ZKGKXH93JeVmNX8qve12ya/bpVkQLwuaXDyud2S/zDGJsfjJH2Y/Edvktb3Hjskax2uxIxtrRKzr03KsqCSEithb5J8na0kPS/pWO/5LZU4PWuQpO8l6xia4/sfLulVeYsx9fAzv0mJGeBVlTgFL7XytKT1Mr4vl/y51Xlfv5MSK3SvnvG662R87dtKrIjdtdp2TfJn/SMlwn6IpNqM1zhdiWsPM2veV4lw2kqJlacfVBB+ayT/zA9V4uyv9ZVYjfy8jNd4SDEJkkp+iDwhT8gTHiE+RKaQKWQKj5AeIk/IE/KkvP/moi6gz76xxG2QrvLGcyTdE8Lr7q7sQeIyH95zl0m6zBuPTP4laZM0X+m3itpVidnJpcnnDst4n+mSPlZiZvgeZayyLOlnSix49KkSqyz/PMf3MVLpt+3aLTleqvTZ564Z2s2T9SyV9KakUzJe7yQlZj8/k/SY8s92T0m+V9pMt/f86f6flRIh1pR87beUMcOa5c8h82fyZxWwCrYyVp6WdFSWP9eZGV/zkqSjc7zeKUrMEH+ixMrUg73n9pT0HyVC8V1JlytxLaX/9bEJkkp/kCfkibePPOFR8oNMIVO8fWQKj5Ie5Al54u0jT/r4YcmCUSAz21WJ07yWSTrEOTc74pJQYcxsMyVCZpCknzjnZkZbEfoKeYK+Rp5UFzIFfY1MqR7kCfpaHPOkx8kTM7tS0v5KLN4zxtt/vBIzfisk3eWc+2VfFgqgMpApAMJCngAIE5kCIJ9C7rYzU4lrmVLMbA9JEyVt65zbWolr2QCgEDNFpgAIx0yRJwDCM1NkCoAcepw8cc49osTCQ74fK7EQzLLkMe/3QW0AKhCZAiAs5AmAMJEpAPIp5MyTbDaXtIuZ/dvMHjazHcIsCkDVIVMAhIU8ARAmMgWApMQ9uov9urWUuD3SDpJuNrONXZYFVMzsWCXu261VV111+y222KLYWgH0gaeffvoD59zQiMsgU4AK0Q8yhTwBKgiZAiAspeZJsZMnCyTNSobGU2a2UtK6Sty2KY1zboakGZI0fvx4N2fOnGJrBdAHzOzNqGsQmQJUjH6QKeQJUEHIFABhKTVPir1sp0mJezbLzDZX4vZCH5RSCICqRqYACAt5AiBMZAoASQWceWJmN0raXdK6ZrZA0hRJV0q60syel7Rc0pHZTl0DgExkCoCwkCcAwkSmAMinx8kT59x3czz1vZBrAVAFyBQAYSFPAISJTAGQT7GX7QAAAAAAAFQFJk8AAAAAAADyYPIEAAAAAAAgj2JvVQygn2lqbtH02fO1sLVNw+rr1DhhtBrGDY+6LAAxRaYACBOZAiAsUeUJkydABWhqbtGkWXPV1tEpSWppbdOkWXMlicYEQK+RKQDCRKYACEuUecJlO0AFmD57fipAurR1dGr67PkRVQQgzsgUAGEiUwCEJco8YfIEqAALW9vSxhvae9pAH3bbDwCFSM8Op8/boiz7AaAwfnYM0TKto8Xd9gNAITJzY5g+kGllWfKEyROgAgyrr0tt7zNgjh4dfLKeHHJ82n4AKJSfHa8N/p4eG3yiRtk7ZAqAonRlxybWopeGfF+zB5+ath8ACuXnxlW1v9ETQ07QVwc8XZY8YfIEqACNE0arrrZGv6u9VJcP+r0k6djO09Q4YXTElQGIo65MeWPIoRpgTpL03sBhZAqAojROGK19aufqH4MbJUlXrPi66mpryBQAveb3KHvU/FeS9EzNtmXJEyZPgArQMG64bhz3gg6qeVSS9PjKMbqvY1tNnz1fTc0tEVcHIG4axg3XizWHpMZbLrtaSzscmQKgKA0bLdPlNVNT48s6D9DggXwMAdB7mT3KHssu0qKOwWXpUUgtoBK8+aTGPndOanjY8tMlBatP82EHQK+ctWZq8+DO89TmaiWRKQCKsGyJ9H/jUsOR7TdIklrbOsgTAL137tDU5p9XNuh1t56k8vQoTJ4AcffJQumqfVPDrqakC6vZA+gVb+LkvIE/01Mdo9KeJlMAFMw5aWpw61B6FAAluekwqXO5JOl521xTlx+c9nRfZwqTJ0CcdbRLv98yNcxsSrqwmj2Aglw4JrX5+kbf1uVLdsp6GJkCoCBn16c2R7Zfn/UQ8gRAQZ64RHrpztRw/7azsh7Wl5nC5AkQV85J530uNWyaOE+W41BWswfQozuOlxa/LUlaWjdMX3/tWzkPJVMA9Mg7i+2O/Z+V5ehSyBMAPXrtIem+M1LDLTv/mvPQvswUJk+AuPJ+m6MprZo+e75clsNMYjV7APk9c630zDWp4T7uj2rr6Mx6KHfIANAjb+JEp7yo3zzwOj0KgOJ8/KZ0zcTUcOcht0XWozB5AsSR35RMapHMcp6i5pRYlRoAslrwtHTHz4LxWYvznvI69cBtyBQAufk9yg9mS2sMo0cBUJzlS6WLtw3GEfcoTJ4AceM3Jcc/Iw1eTVLuU9SGczosgFyWvC9dsWcwPmuxpPx5wgcdADn5Pcp+v5NGfFkSPQqAIjgnnb9BMO4HPQqTJ0Cc+E3JYbdK62ySGjZOGK262pq0wzm9HkBOnR3SbzcLxsmmRCJPABThj18Ktrc9RNrhmNSQTAHQaxlLFHSJMk8G9vk7AAiHP3Gy56+kzfZOe7prpnX67Pla2NqmYfV1apwwmt8SA8ju3HWDbW/iRCJPAPTSPadKi15KbK+yrnTgjLSnyRQAveJ/7jnjXcmCBaejzBMmT4A4uGq/YHuTvaRdf5H1sIZxw2lEAPTMb0omf5z1EPIEQEGeu1n692XB+JevZj2MTAFQEL9HOel5qbb7ZTpR5QmX7QD93T+nSm8+FowPnxVdLQDiz29KTntbGkArAKBI7/xXmvXDYJxxFhsA9Irfoxx5p1S/YXS1ZEHHBPRn8++RHp4WjGlKAJTCb0p++h9pyBrR1QIg3j77UPrzrsGYHgVAKfweZd9p0qhdoqslByZPgP7qg5elG78TjGlKAJTCb0q+c6M0dPPoagEQb50rpOkbB2N6FAClmLF7sL3VROnLP46slHx6nDwxsyvN7H0zez7Lc78wM2dm62b7WgBFav9EumR8MK6gpoRMASLgT5zsdqq0xdejqyVkZAoQgXPXCbbpUQCU4v7J0sLmxHbtqtLB10RbTx6FnHkyU9K+mTvNbENJ+0h6K+SagOq2cqU0zbu+r4KakqSZIlOA8rn+28H2RjtLe5weXS19Y6bIFKB80hac/ii6OvrGTJEnQPm80CQ9fnEwPmNhdLUUoMfJE+fcI5KyJeOFkn4pyYVdFFCtmppbpHPWCsYT50VYTd8gU4DyWfjnb0kv35caN429PMJq+gaZApRHU3NL2sTJXV97UhpQE2FF4SNPgPJ57L5bpVuOTI3j8LmnqFsVm9kBklqcc/81757LOY49VtKxkjRixIhi3g6oCk3NLWq4favUeLP2azRw1lxJqvhb+5EpQPiev/kcjXnn/tR4ZPsNqiNTMo8jT4ACZPYoX132G71951vqGLQmeZJ+LJkCFODeJ57Wvk/8IDWOS4/S6wVjzWwVSWdImlzI8c65Gc658c658UOHDu3t2wFVw29K9lo2XR0aqLaOTk2fPT/CqvoemQL0gbef0ph5v0sNR7bfIElkSgbyBCiM36PMWLGf/uc2JE+yIFOAAnSu0L737ZkaxqlHKeZuO5tIGiXpv2b2hqTPS3rGzNYPszCgqninwZ68/Md61QUzri2tbVFUVE5kChCmpR9Jf9knNexqSrosJFMA9IbXo7zn6nX+isNSY3oUAL3mLTgdtx6l15ftOOfmSlqva5wMkvHOuQ9CrAuoHl5TckfnjrptZfo9zWt6OEU07sgUIETOSReMSg0zmxJJGlZfV86Kyo5MAULkLw4r6UvL/pQ2pkcB0CtepmzW3v2uOv29RynkVsU3SnpS0mgzW2BmR/d9WUCVyGhKTug4vtshna6y1iYjU4A+dHZ9anNk+/VZD2mcMLpc1ZQFmQL0kWsa0obZJmPpUQAUzPvc85VlF6kjy3kc/b1H6fHME+fcd3t4fmRo1QDV5Oy10oY7D7lNau9+qtrwfj4D21tkCtBHvKZkn8E3SO3dD1lrldp+vRBbMcgUoA88dbn02j9TQ3qU1PMjy1QKUFm8HmVS7ala0L5et0Pi0KMUs+YJgFLdcpTkVgbjsxarccJo1dWm3/Kvrram38/AAugH/LPYfvS4frrv2Kx5MuUbW5e5MACxs+Bp6e5fBGN6FACl8HuU8T/Ql75+ZGx7lKJuVQygBE9dLr1wWzA+a7Gk4LZc02fP18LWNg2rr1PjhNH9fgYWQMT8puQbF0vrj1FDcilD8gRAryz9SLoiuAsGPQqAkvg9yqDVpf0vVNcFgXHMEyZPgHJ674Vuv83xNYwbHovgANBP+E3J0C2l7Y9KDckTAL2SseA0PQqAkmSs7ajTF6Q245onXLYDlMvypdKlOwXjjKYEAHolsyn56b+iqQNAZfAWnNaU1ujqABB/f9opfVwhn3uYPAHK5fwNUps7D7lNo067SztPe1BNzS0RFgUgljImTsgUACXxMuVrg6/WqEl3kycAivPYRdL7L6SGldSjcNkOUA5eUzKm8wYtaU2sWN/S2qZJs+ZKUixPXQMQgRu+kzbcsvOvaiNTABTL61GO6JyiFxfXSiJPABRhwdPSA1NSw0rrUTjzBOhrXlNy4OAZWtKR/nRbR6emz55f5qIAxNKzN0j/uyc13HnIbWrr6Ew7hEwBUDCvR7m65iA90pF+9xzyBEDB2j9JW3C6EnsUJk+AvuSfWn/wtWpevFrWwxYmZ2QBIKcPXpGafhyMz1qcMzvIFAA98nuU1YfprM8OynoYeQKgR85J0zYMxhXaozB5AvQVvykZe5i01QEaVl+X9dBc+wFAkrRimXTJ9sE4ufAamQKgKJkLTv/8RfIEQPH8BacnfyypMnsUJk+AvpDZlDT8SZLUOGG06mpr0p6qq61R44T002QBIM2v1wu2vRXryRQAvZbZoyQzhTwBUBQ/U37xsjQgMcVQiZnCgrFA2HI0JVKwONL02fO1sLVNw+rr1DhhdGwXTQJQBn6mnLko7SkyBUCv3HpM+pgeBUAp/B7le7Ok1YJf9lRipjB5AoTpgk3Sx1nuad4wbnisQwNAGflNyQnN0sBB3Q4hUwAU5IXbpLm3BGN6FACl8HuUL/1Y2nSvbodUWqYweQL0QlNzS+7Z0/t+JS39IDg4S1MCAL68meI3JQdeIa29cTRFAoiNnJny8ZvSLUcFB9KjAOhBwT3KoNWkr02LpsgyY/IEKFBTc4smzZqbuuVW2r3K13xFeuL/goNpSgD0IG+m3L5VcOCW35C2/XYUJQKIkVyZYis7NPHvXwgOpEcB0IOCexRJOr2l3OVFhgVjgQJNnz0/673KL7/3KemaA4KdNCUACpArU7o1JYdcV8aqAMRVrkxh4gRAb+XKkwNu3zr9wCrLFCZPgAJlvye5013LjgqGVRYgAIqXLVPeGHJo+g4yBUCBesyUM94rYzUA4ixbnpw/8AoNkAt2VGGPwuQJUKBs9yR/Y8hhwWBKaxmrARB3mZny0KCT0w+owqYEQPEyMyVt4uSnT0m1Q8pcEYC4ysyTvQY8rUMHPhjsqNIehckToECZ9ypPa0pOfVMyi6AqAHHlZ8qJNbdq5ADvt8JV2pQAKJ6fKX6P0jz2HGno6KjKAhBDfp6srw/1l0G/C56s4h6FBWOBAvn3Kn+8/ZvBE0ffL9XVR1QVgLjqypR77rldJy+/NXiiipsSAMXryhR/3aT3hu6scQ0nRlUSgJjqypPf3TtPjy47PniiynsUzjwBeqFh3PD0iZM9zpQ2/GJ0BQGItYbNh+jPyycFO6q8KQFQmswFpz/307sjqgRA3DWMG65Hl30r2EGPwpknqD5571neE/+e5muOkHZr7JsiAcRG0ZninDR942BMUwJUvdB6FIlMARBeppy+sG8KjBkmT1BV8t6zvKcgyWxKTp7bFyUCiJGSMuVs73I/FpwGql5JeTJtRPqYiROg6oX2uefo+6VBq/ZVmbHS42U7Znalmb1vZs97+6ab2Utm9pyZ3WZmLPiAWMh1z/Lps+fn/0J+mxMaMgWVJJRM+cXLLDhdJPIElaToPHngLKnd60voUYpGpqCShNKj7HYqSxR4ClnzZKakfTP23S9pjHNuW0n/kzQp84uA/ijbPcvz7Zckzdg9fUxTUqqZIlNQIYrKFL8pOfRmabX1Qq6qqswUeYIKUVSevPUv6bELgzE9SqlmikxBhSi5R1l1PWmP00OuKt56nDxxzj0i6aOMffc551Ykh/+S9Pk+qA0IXeY9y3var3/PkBY2B2OakpKRKagkvc4UvykZ+z1p8wl9UFX1IE9QSXqdJ+2LpSu9DKFHKRmZgkpSUo8iSY0vh1xR/IVxt50fSLonhNcB+px/z/IudbU1apwwuvvB778k3eMtCEtTUi5kCmKjV5mS2ZQ0/LEPK0MSeYLY6FWeSOnrnNCjlAuZgtgoqUchU7IqacFYMztD0gpJ1+efFmtVAAAgAElEQVQ55lhJx0rSiBEjch0GlEXX4kg9rjq9Ypn0py8FYwKkLMgUxE3BmUJTUnbkCeKm4DyR0jNl8kfdn0foyBTETcGZcuE26WN6lJzMOdfzQWYjJd3pnBvj7TtS0o8k7eWcW1rIm40fP97NmTOnuEqBcvKbkgoPEDN72jk3vszvOVJkCqpFlU2clDtTyBNUHT9TTn5BWrOyryIhU4A+9Mh06cFfB2N6lLyKOvPEzPaVdKqk3QoNECA2/KbkjPeiq6OKkCmoWFcfkD6u8KakPyBPUNH8HuVbV1X8xEl/QKagYrU8U1UTJ2Eo5FbFN0p6UtJoM1tgZkdLukTS6pLuN7NnzeyyPq4TKA+/KfnJv6TaIdHVUqHIFFSNZ66VXn84GNOUhI48QVXxe5QtD5DGHBhdLRWKTEHVWP6ZdPkewZgepSA9nnninPtult1/6YNagGj5TcnXLpDW2zK6WioYmYKq8OGr0h0/C8Y0JX2CPEHVyLz875Bro6mjwpEpqBrnDwu26VEKFsbddoD485uSDb4gfem46GoBEG+dHdIftgvGNCUASlFl6yYB6GN+pvzqg+jqiCEmT4DMpuS4R6KpA0BlOHfdYJsPOQBKwcQJgDD5mXJCs1RTG10tMcTkCaobTQmAMPmZcvrC6OoAEH83HZY+pkcBUAq/R5n4J2ntjaOrJaaYPEH1Omed9DFNCYBS+E3JsQ9Jg1aNqhIAcff8rdJLdwZjehQApfB7lI33kMYdlvtY5MTkCarTHSdIK1cEY5oSAKXwm5K9JkvDxkVXC4B4+/hN6W8/CMb0KABKkXmm/RFN0dRRAZg8QfV55QHpmauDMU0JgFL4TclaI6Vdfh5ZKQBibmWndPG2wZgeBUApWKIgVEyeoLos/Ui67qBgTIAAKEVmU3Lif6OpA0BlOGftYJseBUApmDgJHZMnqB7OSReMCsYECIBS0JQACJOfKae+GV0dAOLv5iPSx/QooWDyBNXj7Ppge/LH0dUBIP5YcBpAmPyJk6Pukurqcx8LAPnMu0Oad3swpkcJDZMnqA5+U3LKS9IA/uoDKNLsM1hwGkB4/B7lyz+RRn4luloAxNsn70g3Hx6M6VFCxSdIVD6/Kfn21dIaG0RXC4B4e+Nx6clLgjFNCYBSpF3+Z9K+UyMrBUDMrVwp/X6LYEyPEjomT1DZ/KZky29IWzdEVwuAeGtfLM38ejCmKQFQim7rJrVGUweAynDOWsH2FPKkLzB5gsqV2ZQccl00dQCoDNNGBNtMnAAoBQtOAwiTnymNr0lm0dVSwZg8QWWiKQEQJj9TJn8UXR0A4u+CTdLH9CgASuH3KN+7VVp1ndzHoiRMnqDyMHECIEx+ppz4X2lATXS1AIi3h6ZJSz8IxvQoAErh9yjjDpc23Tu6WqoAkyeoLFfskz6mKQFQCr8pabhUWmtkZKUAiLmWp6WHvAVh6VEAlCLzF8YTL8l+HELD5Akqx5yrpAVPBWOaEgCl8JuSUbtKYw+NrhYA8bb8M+nyPYMxPQqAUnCmfSSYPEFl+PBV6c6TgjEBAqAUmU3JkX+Ppg4AleH8YcE2PQqAUjBxEhkmTxB/nSukP2wXjAkQAKWgKQEQJj9TzlwUXR0A4u8P26eP6VHKiskTxN+53orSBAiAUjBxAiBMfqb8bI40cFB0tQCItyf+IH34SjCmRyk7Jk8Qb35TMqklujoAxN91B6WPaUoAlMLvUfb7nbTuZtHVAiDe3ntBuu/MYEyPEgkmTxBfflNy9APS4NWiqwVAvD13i/TKA8GYpgRAKfweZYOx0g7HRFcLgHjraJcu3SkY06NEZmDUBQCZmppbNH32fC1sbdOw+jo1ThithnHD0w/ym5Jdfi5tuEN5iwQQCwXlSevb0izvgw1NCYAcet2jSNJxD5evQACxUVCeSNJ5nwu26VEi1eOZJ2Z2pZm9b2bPe/vWNrP7zezl5H/X6tsyUS2amls0adZctbS2yUlqaW3TpFlz1dTsXZLjNyWDVpf2mlz2OlE8MgXlUlCerOyULhoTjGlKYodMQbn0ukeRyJSYIU9QLgXliZSeKWe8V9Ya0V0hl+3MlLRvxr7TJP3DObeZpH8kx0DJps+er7aOzrR9bR2dmj57fmKQ2ZScvqBMlSFEM0WmoAx6zBNJOmftYHtKa5kqQ8hmikxBGfS6R2HiJI5mijxBGRTUo/iZ8qPHpdohZaoOufQ4eeKce0TSRxm7J0q6Orl9taSGkOtClVrY2pZ7P01JRSBTUC5580RKz5RT35DM+r4ohI5MQbnkzZRbvp++kx4llsgTlEuvepR9zpHWH5P1eJRXsWuefM45944kOefeMbP1ch1oZsdKOlaSRowYUeTbodLkusZvWH2dWrKEyUtDjkrfQVNSacgUFK23eTKsvi69KTny71IdZ2FXmIIyhTxBNr3NlINXf056YVawgx6l0tCjoGgl9yhrjZJ2PrGMFSOfPr/bjnNuhnNuvHNu/NChQ/v67RAD+a7xa5wwWnW1NWnHnznoJg3W8mAHTUlVI1Pg622e1NXW6PH2bwY7vnicNGrX8haNfoM8QabeZsqGtZ/qNx3Tgh30KFWNTIGv5B5Fkk58tnwFo0fFTp68Z2YbSFLyv++HVxIqXb5r/BrGDdfUA7fR8Po6maTdh7yqYwbckTquaeK8MleLMiFTUJTe5Mnw+jq9WHNI2rFNG/DbnApFpqAovcqUNYfo0ZrjUsftPOS27os9ohKQJyhKyT0Kn3v6nWIv27lD0pGSpiX/e3toFaHi9XSNX8O44WoYN1x3PvU/7X/3oannR7bfoLpZc1PHoKKQKShKoXkiqdu6SWRKRSNTUJRiM2Vk+/VSe+K3yl3HoWKQJygKPUrlKeRWxTdKelLSaDNbYGZHKxEe+5jZy5L2SY6Bggyrryto//5375DaHtl+g6Qsq1AjdsgUhKnQPMnWlEhkSiUgUxCmYjJl+/ZLJSUWnCZT4o08QZgKzpO/p58FS4/Sf/V45olz7rs5ntor5FpQJRonjNakWXPTTmOrq61R44TRwUFeU7JJ+7VpX59rFhfxQKYgTAXlyXkbpH1NV1PShUyJNzIFYeptj3LU8kZ9qPTJWTIlvsgThKmgPHnpbunpmakhPUr/VuxlO0DRuk49y7bytKS0pmTXZReqU+mLKeWaxQVQfXrMk9lnSB1LU8dnNiUSmQIg0Jse5e4Be+ihleO6vQaZAkAqIE+WLJJuCubr6FH6PyZPEIm0a/x8XlPy3JhJWvTfYVK+2VoAVS9nnrz9lPTkJalh08R5quvpN0AAql4hPYokLf/GH8kUAHnlzBPnpN9umhrSo8QDkyfoP/ymZO2Nte23TtPUTbLfGx0A8lr+mfSXfYLxWYvVkNwkUwD0WsbECZkCoCRn1wfbkz9Ww4DEUqTkSf/G5AnKqqk5x2RIZlNyQrOkPLO1AKA8mXL+sOCgsxanNskUAPlkzZTbt0o/iEwBUICCPvec9LyUnDghT/q/Hu+2A4SlqblFk2bNVUtrm5ykltbkLf2y/DYHAHpSUKac+X5k9QGIl2yZslfT9ukH0aMAKEBBPUrDpVL9hpHViN5j8gRlM332/LTr+CTpxZpD0g+iKQFQoB4z5UePSQMHl7kqAHGVmSkn1MzS6ubd6YIeBUCBeuxRNvyyNPbQMleFUjF5grLJvNVW06Az0w+gKQHQC5mZ8sYQrwnZ80xp/W3KXBGAOPMzZWt7Q6fU/i14kh4FQC/k7VEk6ejZZawGYWHyBGXj32rrqJp7NXbAa8GTNCUAesnPFL8pWazVpF0boygJQIx1ZcpgLdddg09P7d95yG1RlQQgpnL1KJL43BNjTJ6gbBonjFZdbY02snd1Vu01qf1NE+dFWBWAuOrKlMym5J8Tn4qoIgBx1pUp84ccldq3ZedfuVUogF7L1aPwuSfemDxB2TSMG66p39xaDw8+JbWvaeI8VpUGUJSGccO7rZtEpgAoVmam7DX4Rk09cBsyBUCvNYwbrudrv5e2jx4l/rhVMcqq4Y4xwWBKqxrMoisGQLxluVNXQzSVAKgEfqZ8/x79Y6OdoqsFQLz9c6pqVnYEY3qUisCZJygfvyn5xcsSEycAinXZV9LHXD8MoBR+j7L9URITJwCK9e7z0sPTgjE9SsVg8gTl4Tcl37pKWm296GoBEG/N10nvzg3GNCUASpF5Fts3Lo6mDgDx19khXbZzMKZHqShMnqDv+U3JyF2kMQdGVwuAeFu8QLr9p8GYpgRAKbJc/gcARTt33WCbPKk4TJ6gb2U2JUfdGU0dAOLPOenCrYMxTQmAUjBxAiBMfqac+mZ0daDPMHmCvkNTAiBMZ9cH21Nao6sDQPzRowAIk58ph90q1dXnPhaxxeQJ+gZNCYAw+Zly8gssOA2geDcfmT6mRwFQCr9H2apB2mzv6GpBn2LyBOE7K32mdecht6mpuSWiYgDEnteUnDfwZ2p6jYkTAEWaf680ryk1pEcBUJKMXxg3bXZ+RIWgHAZGXQAqzP1TJLnUcGT7DVJ7mybNStwZo2Hc8IgKAxBLXlPySOc2urx9J9WRJwCK8dmH0o2HpIb0KABKkjFxMrL9BnqUCseZJwjPgjnS4xelhiPbb0htt3V0avrs+VFUBSCuMpqSIzomSSJPABTBOWn6xqkhPQqAkmSZOJHIk0rH5AnCsXypdMVeqaHflHRZ2NpWzooAxFmOpqQLeQKgV7wFp0e1X9ftaTIFQMHoUaoWkycIx/kbpDZ3HnJb1kOG1deVqxoAcdZDUyKRJwB6wcuUgwb9WS5L+0umACjIHcenDelRqktJa56Y2cmSjlFikYu5kr7vnGsPozD0P03NLZo+e74WtrZpWH2dGieMTlzP53/QOeM9NT7/oSbNmqu2js7U7rraGjVOGB1B1YgTMqW6ZM2U+3dJP2biPNWRJygCeVJdCupRDviDDrddNI9MQRHIlOqSNVPWfEV65prgGHqUqlP0mSdmNlzSCZLGO+fGSKqR9J2wCkP/0tTcokmz5qqltU1OUktrcoE1vyk57hGpdogaxg3X1AO30fD6Opmk4fV1mnrgNiychLzIlOqSLVPeuO1saemHwUFnLSZPUBTypLoU1KNsMFba7ggyBUUhU6pLtkw5b9a/pWsOCA6iR6lKpd5tZ6CkOjPrkLSKpIWll4T+aPrs+WmzqpL0Yk2wYr12P13a4AupYcO44QQHikGmVInMTNnC3tJJA/4aHHDW4tQmeYIikSdVosceRZKOezi1SaagSGRKlciWKf+p+X4woEepWkWfeeKca5H0W0lvSXpH0mLn3H1hFYb+JXPhozeGHBoMBq8h7X5qmStCpSFTqoufKYPUoXsHnxY86TUlQDHIk+qSt0eRyBSUjEypLnkz5VcfCtWrlMt21pI0UdIoScMkrWpm38ty3LFmNsfM5ixatKj4ShEpf+Gjbk3JpLfLXA0qEZlSXfxM+d+QI1PbuRacBnqDPKkueXsUJk4QAjKluuTKlO8O+oNUU+qFG4izUu62s7ek151zi5xzHZJmSdop8yDn3Azn3Hjn3PihQ4eW8HaIUuOE0aqrrenWlDRNnBdRRahAZEoVyZYp23VezSJrCAt5UkW68uTBQaek7adHQYjIlCqSrUc5f+WROuRre0VYFfqDUiZP3pL0ZTNbxcxM0l6SXgynLPQ3DeOGd7t+uGniPK7xQ5jIlCqSmSnHDpqqyQfuQKYgLORJFWkYN1x//cKz2njAu6l99CgIGZlSRTJ7lLdtA231zVPJFBS/YKxz7t9m9jdJz0haIalZ0oywCkM/c/3B6eOzFqshmkpQociUKuPfBWPHn2nGhJ9EVwsqDnlSZT54Wds+PzUY06MgZGRKlfF7FEkbTnlJG0ZUCvqXki7acs5NkTQlpFpQZlnvX55tRvWFJunl2cGY64fRR8iUeCs4UzKaEk04rzwFoqqQJ/FWcJ50rpAuGR+M6VHQR8iUeCu6RyFT4GHFmyrVdf/yrttwtbS2adKsuZKUHiSfvifdEizmSIAAyKbgTKEpAdCDgvNEks5dJ9gmTwBkQY+CsJSy5gliLNv9y9s6OjV99vxgh3PS7zYPxgQIgBwKyhSaEgAFKChPpPRMOe2tMlQGII4KypTLMxaDpUdBFkyeVKnM+5dn3X92fbA9pbWPKwIQZz1mChMnAApUUI/iZ8rhTdKQNbt/AQCogExpvk5qmRM8QY+CHJg8qVL+/cuz7vebklNelMzKUBWAuMqbKbcek76TpgRAHr3qUcYdLm2yRxmqAhBXeTOl9S3p9p8GO+lRkAeTJ1Wq6/7lvrraGjVOGJ3elBx4hbTGsDJXByBucmXKb8e+K829JdhJUwKgBwX3KJI08ZIyVgYgjnJmylc3ky7aJthJj4IesGBslepaHKnbqtO3bxUcNHIXadtvR1QhgDjJlimn7/E57XjPIcFBNCUAClBQjyKRKQAKUlCmsEQBCsDkSRVrGDc8/wrTR91Z3oIAxFreTOFDDoBe6LFHIVMA9ELeTGl8lSUKUBAu20ECTQmAMPmZMvmj6OoAEH/0KADC5GfKd2+SVl03uloQK0yegKYEQLj8TDnhWWlATe5jASCfO09JH9OjACiF36Ns1SCN/lp0tSB2mDypdr8dnT6mKQFQCr8p2f8iae1R0dUCIN7eeEya85dgTI8CoBSZvzA++Opo6kBsMXlSzf5zhbTk3WBMUwKgFH5Tsv420vjvR1cLgHhbtkSauV8wpkcBUArOtEcImDypVh+8LN3182BMgAAoRWZT8qPHoqkDQGWY6i/sSI8CoARMnCAkTJ5Uo84O6ZLxwZgAAVAKmhIAYfIz5VcfRFcHgPi7aJv0MT0KSsDkSTU611tRmgABUAomTgCEyc+U45+RamqjqwVAvD12odT6VjCmR0GJmDypNn5TcvrC6OoAEH9XH5A+pikBUIrMBafX2SS6WgDE27tzpQfOCsb0KAgBkyfVxG9KfvigNGjV6GoBEG/z7pBefzgY05QAKIXfo3z+iyw4DaB4He3SZV8JxvQoCAmTJ9XCb0p2O00avn10tQCIt8Ut0s2HB2OaEgClyLz875j7o6kDQGU473PBNj0KQsTkSTXwm5INxkp7TIquFgDxtrJTunCrYExTAqAUrJsEIEx+ppy5KLo6UJGYPKl0mU3JcQ9nPw4ACnHO2sE2H3IAlIKJEwBh8jPlZ09LAwdFVwsqEpMnlYymBECY/EyZtCC6OgDE310/Tx/TowAohd+jHHCJtO6m0dWCisXkSaW6eGz6mKYEQCn8puSYB6XBq0dXC4B4e/kB6T9XBGN6FACl8HuUjfeQtjs897FACZg8qUSPTJc+fj0Y05QAKIXflOx5pvR5FpwGUKQl70vXHxSM6VEAlCLzTPsjmqKpA1WByZNKs+Bp6cFfB2OaEgCl8JuStTeRdm2MrhYA8eac9NvNgjE9CoBSsEQByqykyRMzqzezv5nZS2b2opntGFZhKMKyT6Ur9gzGBAhihkzpZzKbkhOeiaYOoAjkST90dn2wPaU1ujqAIpAp/cy0jdLHfO5BGQws8esvlnSvc+5bZjZI0ioh1IRiTf18sE2AIJ7IlP6C3+Yg/siT/sTPlFPfkMwiKwUoEpnSX9w/WWr3JmDpUVAmRU+emNkaknaVdJQkOeeWS1oeTlnoNb8pmfyRJKmpuUXTZ8/XwtY21a9SK+ekxW0dGlZfp8YJo9UwbnhExQLdkSn9yJ93Sx8nmxIyBXFBnvQzfo/y/XukurUkBZnS0tqmGjN1Oqfh5An6ITKlH3njMenxi4MxPQrKqJTLdjaWtEjSVWbWbGZXmNmqmQeZ2bFmNsfM5ixatKiEt0NOflNy8jxpQI2amls0adZctbS2yUn6eGmHWts65CS1tLZp0qy5ampuiapiIBsypT/416XSO88GY68pIVMQI+RJf+H3KF85RdpoJ0npmSJJnc5JIk/Qb5Ep/cHSj6SZ+wVjehSUWSmTJwMlbSfpUufcOEmfSTot8yDn3Azn3Hjn3PihQ4eW8HbIym9KDrlOWjMxqzp99ny1dXTm/LK2jk5Nnz2/r6sDeoNMidq7z0v3ej9y7zRYMgUxQ570B36Pssq60t5TUsN8mUKeoB8iU6LmnHTBqGBMj4IIlDJ5skDSAufcv5PjvykRKigXvynZ9hBpy2+khguTv8nJp5BjgDIiU6K0fKl02c7BOOP6YTIFMUOeRC1z3aRfvpo27CkvyBP0M2RK1PIsOE2PgnIpevLEOfeupLfNbHRy116S5oVSFXqW2ZQcOCNtOKy+rseXKOQYoFzIlIidv0GwnWXhNTIFcUKeRKyABad7ygvyBP0JmRIxP1MaX+224DQ9CsqlpFsVSzpe0vVm9pyksZLOL70k9KiApqRxwmjV1dbkfIm62ho1Thid83kgImRKFPxM+dUHWQ8hUxBD5EkUrv5G+jjHXTDyZQp5gn6KTImC36N8b5a06rrdDqFHQbmUdKti59yzksaHVAsKUeDtQ7tWlGbVacQJmRIBP1NOeFaqqc16GJmCuCFPIvD0TOn1R4JxntuH+pnC3XYQB2RKBPwe5YvHSZvulfUwehSUS0mTJyizpp+kj3u4p3nDuOEEBYDc/KbkwCuktUflPlZkCoA8PnhZ+vuJwbiHHkUiUwDk4fcoA2qlr1+Q93DyBOVQ6mU7KJeX7paevT4YF9CUAEBOflOyxf7Stt+OrhYA8bZimXSJ9wt5ehQApcg8035y9kuKgXJj8iQOPnlHuum7wZimBEApMpuS71yf/TgAKMSv1wu26VEAlKLAJQqAKDB50t+tXCn9fotgTIAAKAVNCYAw+Zly5vvR1QEg/m46LH1Mj4J+hsmT/u6ctYJtAgRAKc5ZJ31MpgAohT9x8rM50sDB0dUCIN6eu1l66c5gTI+CfojJk/7Mb0pOezu6OgDE392/lFauCMY0JQBK4fco37hYWnez6GoBEG8fvyHN+mEwpkdBP8XkSX/lNyVHPyANWSO6WgDE2yv/kJ76czCmKQFQCr9HGbmLtP1RkZUCIOY6O6SLvxCM6VHQjzF50h/5Tcnup0sb7hBdLQDi7bMPpOsODMY0JQBKkblu0lF3Zj8OAApx7rrBNj0K+jkmT/obvympHyHtfmp0tQCIN+ek6ZsEY5oSAKVgwWkAYfIz5fR3oqsDKBCTJ/1JZlNy0txo6gBQGc6uD7antEZXB4D4m75p+piJEwCl8D/3/OhxadAq0dUCFIjJk/6C3+YACJOfKb98XTKLrhYA8fbA2dJni4IxPQqAUvg9yoSp0vpjoqsF6AUmT/qDy/dKH9OUACiF35QcdZe0ytrR1QIg3t54XHrs98GYHgVAKfweZYOx0o4/ia4WoJeYPInav2dILXOCMU0JgFL4TcnOJ0ojvxJdLQDire1jaebXgzE9CoBSZJ5pf9zD0dQBFInJkyi9N0+6pzEY05QAKIXflAxZU9rnnOhqARBvzkm/GRmM6VEAlIIlClABmDyJSkebdOmOwZgAAVCKzKbktLeiqQNAZfAXnJ78cXR1AIi//9sufcznHsQUkydROW/9YJsAAVAKfpsDIEx+pvziZWkA7SKAIj3yW+mjV4MxPQpijP8bRsFvSn71QXR1AIi/a7+ZPqYpAVAKv0c57FZptfWiqwVAvLU8LT14bjCmR0HMMXlSbn5TcvwzUk1tdLUAiLdnrpVefTAY05QAKIXfo+xwjLTZ3tHVAiDeln0qXb5nMKZHQQVg8qSc/Kbkm3+W1tkkuloAxNsHr0h3/CwY05QAKEXa5X8m7fe7yEoBUAGmfj7YpkdBhWDypFz8pmSzCdIXvhNdLQDibcVy6ZLtgzFNCYBSdFs3qTWaOgBUBj9TJn8UXR1AyJg8KYfMpuSwm6OpA0Bl+PXQYJuJEwClYMFpAGHyM+XkedKAmuhqAULG5ElfoykBECY/U854L7o6AMTfzUekj+lRAJTC71EOvkZac3h0tQB9oOTJEzOrMbNmM7szjIIqynkbpI9pSoAekSl5+E3JT5+SaodEVwsQA+RJHnP/Js27PRjTowA9IlPy8HuUbQ+RtpoYXS1AHwnjzJMTJb0YwutUlntPlzqWBmOaEqBQZEo2flOy3++loaOjqwWID/Ikm4/flG49OhjTowCFIlOyyTzT/sAZ0dQB9LGSJk/M7POS9pN0RTjlVIaH/3GX9K8/psZNE+dFWA0QH2RKd03NLWlNyQfrjJd2ODrPVwCQyJNcbn/6TenibVNjehSgMGRKd5k9iiQmY1HRSj3z5CJJv5S0MoRaKsJd/56n3R49NDUe2X6DJs2amwgXAD0hUzxNzS1quH2rtH27vN9IngCFIU8yNDW3aOLfg4kTehSgV8gUT1Nzi77YtGvavi07/0qeoKIVPXliZvtLet8593QPxx1rZnPMbM6iRYuKfbt4cE773bNjajiy/QZJUltHp6bPnh9VVUAskCndZU6cjGy/gTwBCkCeZOdnyuj2mZLoUYBCkCndfXLnrzTMPkiN6VFQDUo582RnSQeY2RuSbpK0p5ldl3mQc26Gc268c2780KFDM5+uLGfXpzZHtl+f9tTC1rZyVwPEDZniu3LftGHXZKxEngAFIE8yeafW77VsupZpUGpMpgA9IlN8rz2sIzpvTQ3pUVAtip48cc5Ncs593jk3UtJ3JD3onPteaJXFjdeUbNs+Q5KlPT2svq7MBQHxQqZ4/nWZ9NaTqaHflEjkCdAT8iSD16M0dhyrV1367UPJFCA/MsWzZJF0zQGpIT0KqkkYd9uB15Q8vMuN6qhNXziprrZGjRO4MwaAAix8Vrr31NRwy86/pj1NngDoFa9HeWf9PXXngL3SniZTABRs5Urpt5umhvQoqDYDw3gR59xDkh4K47Vix2tKLh34PV1wv1P9KgM0eOAALW7r0LD6OjVOGK2GccPzvAgAX9VmSvsn0ozdUsNR7TeQJ0CJqjZPpHOTue0AABeLSURBVG53wdjpjWPIFKBEVZ0p56yV2hzVfj15gqoTyuRJ1fKaknlupH6z5OuSpI+XdqiutkYXHjKWAAFQGOekaRumhl2nwZInAIqSMXFCpgAoiZcpW7VfKScjT1B1uGynWBlNydeXnZ82ZrVpAL2SZ8Fp8gRAr1z/7bRh5poEZAqAXvE+93xt2VQt1ZDUmDxBNWHypBjnrJs2HJXRlHRhtWkABfGaku3bL1PmgtMSeQKgQE/PlF6+LzWkRwFQEq9HmdxxlF50G3U7hDxBtWDypLfuOEFa2RGMz1qcc1VpVpsG0CP/LLaj7tKQ+s9lPYw8AdCj9+ZJfz8xGNOjACiF36OM2k3/WH1i1sPIE1QL1jzpQVNzi6bPnq+FrW06ZPXnNK3j6uDJsxZLkhonjNakWXPV1tGZeorVpgFk42fK60MODZ7Y9ZfSyK+ocUILeQKgIH6ejFpzgB5c9p3gSXoUAL2Us0eRpCPvUGMzPQqqG5MneTR5ATFUH2tax7TgyWRTIim1QFJX2LDaNIBs/Ex5w2tKlqw6QqvteYYk8gRAYZoyPsRkmziRyBQAhfEz5ZFBJ6Y/mcwU8gTVjsmTPKbPnq+2jk7VqFP/GfLT1P6dh9ymxzOObRg3nOAAkFdXpryR8ducCZ0Xp2UKeQKgJ115IiktU3YZ/Dc9mnEsmQKgJ12ZctrAGzViwKLU/szPPeQJqhmTJ3l0LX706pDDU/tGtt8ga2dRJAC9t7C1TX8ddE7aPjIFQDG6ehR/4mS79sv0cfvyqEoCEGMLW9u0+4Bm/Wjg31P76FGAdCwYm8ew+rq0pmSz9mtS+wGgt3682iP60oCXUuOu24eSKQB6K7NHOXDZWfpIa5AnAIqy7ZpLNXPQ9NSYHgXojsmTPB5v/2Zq+0vtl6hDA1kUCUBxWp7RL1dclhp2NSVkCoBi+D3K6R1H6xm3OXkCoDidHbp92TGpIT0KkB2X7eTi3Zrr+Nqz9X772hrOokgAirH0I+nyPVLDnYfcJmtnoTUARfJ6lH8M2Ek3du5FjwKgeOeum9qkRwFyY/IkG/+e5nv+Sn/Y9ST9IbpqAMTZypXSBaOC8VmLuy04DQAFO29Y2nCvyffo9YhKAVAB/M89v/pQj9fw8RDIhX8dmfwAGbGjtOsvJKXf95yZWAAFO2etYNu7fSiZAqDXZh0ndXwWjJOZQp4AKIr/uecXr0jJiRMyBciOyRPfH7ZPH//gXknp9z2XpJbWNk2aNVeSCBIAuflNyZnvpzbJFAC99uwN0nM3BWNv4oQ8AdBrfo/yg9nSakMlkSlAPiwY2+W+M6UPXwnG3m+Iu+577mvr6NT02fPLVR2AuPGbkpPmSgMHp4ZkCoBeeXeu1PTjYEyPAqAUfo+y7zRpxJdTQzIFyI3JE0l66W7pCW9VE68pkRL3Pc8m134AVc5vSg67VaofkfY0mQKgYO2Lpcu+EozpUQCUwu9RNvuq9OUfpz1NpgC5MXny0WvSTd8NxhlNiZT7/ubc9xxAN35TsmujtNne3Q4hUwAUxDlpmjf5So8CoBS/3Tx9fNgt3Q4hU4DcqnvypKNN+r9xwThLUyJJjRNGq662Jm0f9z0H0I0/cbL+ttKeZ2Y9jEwBUJCz64PtKa1ZDyFPABTkzlOkJe8FYz73AL1W3QvGnrd+sJ0jQKRgcSRWnQaQ0593Sx//6NGch5IpAHrkT8ae/o5klvUw8gRAj56/VZrzl2DM5x6gKNU7eeI3JTl+m+NrGDec0ACQ3YPnSe88G4zzNCVdyBQAOfk9yvHPSINWyXs4eQIgp0Xzpb/9IBjTowBFq87Ldvym5NQ3c/42BwB69MoD0iMXBOMCmhIAyMnvUQ65Xlpnk+hqARBvy5ZIf/xiMKZHAUpSfZMnflNy3KNSXX3uYwEgn9a3pesOCsY0JQBK4fcoO/5M2nL/6GoBEG/OSVO9s0foUYCSVdfkid+UHHCJtMG20dUCIN5WLJMuGhOMaUoAlMLvUdbeWJpwXnS1AIi/AhacBtA7RU+emNmGZvZPM3vRzF4wsxPDLCx0flOyzcHSdodHVwuAbmKXKb9eL9hm4gToV2KXJ1d9PX18QnM0dQDIKnaZ4n/umdTCEgVASEpZMHaFpJ87554xs9UlPW1m9zvn5oVUW3jOXjvYHlgnHXR5dLUAyCU+meI3JZM/jq4OALnEJ08e/b305uPBmMlYoD+KT6b4PcpPn5IGrxZdLUCFKXryxDn3jqR3ktufmtmLkoZL6tMQaWpu6d2ts275vuQ6g/GZ7/ZleQCKFEWm9DpPpPSmpPE1aUB1Xf0IxEFsepTXH5X+cXYwZuIE6Jdi2aN860pp6Oi+Kg+oSqHcqtjMRkoaJ+nfWZ47VtKxkjRixIiS3qepuUWTZs1VW0diMqSltU2TZs2VpOxhMucq6YVZwZimBIiFcmRKr/NESm9KfvigtOo6Rb8/gPLotz3KJ+9IV3sLwtKjALEQix5l/A+kMQdlPw5A0Ur+lamZrSbpVkknOec+yXzeOTfDOTfeOTd+6NChJb3X9NnzUyHSpa2jU9Nnz+9+8MJm6c6TgjFNCRAL5cqUXuWJlN6U7Pc7afj2Rb83gPLotz1KZ4f0+y2CMT0KEAux6FFWXU/a/8Ki3xtAbiVNnphZrRIBcr1zblZPx5dqYWtbYfuXfiTN2D0Y05QAsVDOTCk4T6T0pmTLA6QdjumjqgCEpd/2KJJ07rrBNj0KEAv9tke57lvp48aX+6AiAFJpd9sxSX+R9KJz7vfhlZTbsPq6nvevXCldMCoY05QAsVDuTCkoTyRpasZpt4dc20cVAQhLv+1RpP9v795j5CrvM45/f6xtvFxim1uCDQ04BSuh/AFyCJeIRqTFxqqA0iiCpEBCJBQIFaEC1Q6BNrQNENqojUTSQrBEEwrkAq7VQg0qRVWbmmJuMQYMDgWCTcEFbKDYsbHf/jFnd2aXmdnLnD3zzu73I6085zI6j85Onrz7MvPO0MnYq16fwESSypLtGGX192DD/fVt/+6RJlQn7zw5CTgXOCUiHi9+loz0pE5csWgB/dP7huzrn97HFYsaFkO6Zk79sQUi9ZJKO2VUfbLiK/Crhh6xU6RekecYpXHi5PLnoK+UpeckTbz8xigvPQT/vLS+7RhFmnCdfNvOvwOVfmn4wAJJLVeebhyUXPlqldEkdajqThmxT564Ex7/Yf0JDkqknpH9GOWCVbDPQVXGk9SB7MYo72yG5afWn+AYRapEz/0njzOPmdd8lenGQcmlT8D0mdWFktSTWvbJq+vg7gvr2w5KJI3CqMYoi74Jv3Z8daEk9aSWfbJ7F/zFr9e3HaNIlen423ay0Dgo+dyPYM5hXYsiqcdtfwu+d2J920GJpE40jlE+8mk44SvdyyKp912zX/2xYxSpUr0/edI4KPnkZXDkou5lkdTbUoLrDq1vOyiR1Im//OjQ7XMn/Et/JE1mjX/3fH1z93JIU1TWH9tZ8djG1p/1g6EFctBR8Ft/UnVEST1ixD4B+Mbs+uM/3lJtQEk9ZcRO+afL4e1N9W0nYyW1MKoxSuPfPX/4NEybUW1ISflOnqx4bCPL7lrLtp27ANi4ZRvL7loLFIso3XzK0Cdc/LOqI0rqESP2CQwdlHztFYhK15qU1ENG7JR1d8PDN9ef4MSJpBbGPEY5byV8YG7VMSWR8cd2bli1frBEBmzbuYsbVq2HB6+DjY/UDzgokdRG2z6BoYOSP3gUZuxVYTpJvaZtp2x+Fn78hfoBxyiS2hjTGOWUq2D+b1aYTlKjbCdPNm3Z1nT//lvXwYPX1nc4KJE0glZ9smnLtqGDks/+APb/SEWpJPWqVp3y+patcOPH6zsco0gawajHKId+Ak6+vKJUkprJdvJk7uz+9+3bn62s3PPr9R0OSiSNQrM+AfjvmZ+rbxx/MXzs9IoSSeplzTsl8czML9Q3HaNIGoVWY5Tv7nXT0B1fuq+CNJLayXby5IpFC+if3je4PY33eGTmRfUTHJRIGqXhfQJw757L6hu/8Xuw+FokaTSadcoLMz9f33DBaUmj1KxPzp/xAKftfrC+w797pCxku2DswAJJAytPb5h5Xv2gBSJpDIb3yQX7ruajO1+sHZw5Cz6zvIvpJPWa4Z0y5F1sLjgtaQyG98nxH3iDb+z4fv0E/+6RspHt5AnUyuTMY+YN/bzf1W92L5CknjXYJy+thuXfqe2csQ8sfam7wST1pKZjlMuecsFpSWM22CfvvgHfOrx+wIkTKSvZfmxn0JCvD90Ee+QfWVKmtrwEyxfVt7+2sXtZJPW+Pz+4/viin8Gsed3LIqm37drpxImUuazeebLisY2Db1mbO7ufH3zoR8wfOHjZUzBj727Gk9RjGjvlw7Om8eCvPls/6KBE0hg1dsoX932Iq3e+Wztwzh3wwaO6G05STxn+d89/bP/d+kHHKFKWspk8WfHYRpbdtXbwe84Xv/1T5m+/vXbwkjX+1xxJYzK0U5ITJ5I60tgpJ+7xJFfv/GsA1hx7PQsXnNbldJJ6yfC/e1ZsOx8GlkpywWkpW9l8BuaGVesHC2TJHqu5avoPATh3xl/BAUd0M5qkHtTYKY3fgnHSzLu7FUlSDxvolAXxEn8/45sAXLTjUi596sguJ5PUaxrHKLdP/zMOjLcAOHnPn7jgtJSxbN55smnLNgAOYCvfnVFbzPGcHVeyevtB3YwlqUcNdMrfTv/24L7Dtt9GbN/WrUiSetimLdvYg92s2nMpAH+68/e5d/cniC12iqSxGRijnNe3ihP6ngLgqO238O72Hd2MJWkE2bzzZO7sfgD6YzsAl+24iP/cfdTgfkkai4HuWLu7tvjakdtvBcJOkTQuc2f3M5MdbE6zuPm9Jdyya8ngfkkai4HeeDPtC8DHt9/I/9Fvn0iZi5RSZRdbuHBhWrNmTdNjwz/7B9A/vY9rzzoaYMiCSlcsWjD4neiSOhMRj6SUFnY7x3jYKVJ+erVT2vUJ2ClSt0zGTrFPpO7otE+y+djOQCkMLwtgSLls3LKNZXetHfIcSRrOTpFUJjtFUlnsE6k3ZTN5ArVSGF4MJ133wJBZWYBtO3dxw6r1loiktuwUSWWyUySVxT6Rek82a560sqnFQmyt9ktSO3aKpDLZKZLKYp9Iect+8qTVwkkuqCRpPOwUSWWyUySVxT6R8pb95MkVixbQP71vyL7+6X2DnwuUpLGwUySVyU6RVBb7RMpbR5MnEbE4ItZHxIaIWFpWqEZnHjOPa886mnmz+wlg3ux+rj3raD/3J01CdoqkslTRJ2CnSFOFYxRJ414wNiL6gBuB3wZeBh6OiJUppafKCjeg2YJKkiYXO0VSWarsE7BTpMnOMYok6OydJ8cBG1JKz6eUdgB3AGeUE0vSFGSnSCqLfSKpTHaKpI4mT+YBv2zYfrnYN0REXBgRayJizebNmzu4nKRJzk6RVBb7RFKZ7BRJHU2eRJN96X07UroppbQwpbTwwAMP7OBykiY5O0VSWewTSWWyUyR1NHnyMnBow/YhwKbO4kiawuwUSWWxTySVyU6R1NHkycPAERFxeETMAM4GVpYTS9IUZKdIKot9IqlMdoqk8X/bTkrpvYi4BFgF9AHLU0rrSksmaUqxUySVxT6RVCY7RRJ0MHkCkFK6B7inpCySpjg7RVJZ7BNJZbJTJEVK71vraOIuFrEZeLGyCw51APC/Xbr2SHLOBubrRM7ZoJZv75RST65qZqe0lHM2yDtfztkg73wD2T7ci51in7SVc76cs0He+XLOBnZKJ3L+3eacDfLOl3M2yDtfKX1S6eRJN0XEmpTSwm7naCbnbGC+TuScDfLPl7Oc713O2SDvfDlng7zz5Zwtd7nfu5zz5ZwN8s6XczbIP1/Ocr53OWeDvPPlnA3yzldWtk4WjJUkSZIkSZr0nDyRJEmSJElqYypNntzU7QBt5JwNzNeJnLNB/vlylvO9yzkb5J0v52yQd76cs+Uu93uXc76cs0He+XLOBvnny1nO9y7nbJB3vpyzQd75Ssk2ZdY8kSRJkiRJGo+p9M4TSZIkSZKkMZt0kycRsTgi1kfEhohY2uT4nhFxZ3H8oYg4rKJch0bEv0bE0xGxLiIubXLOpyJia0Q8XvxcXUW2huu/EBFri2uvaXI8IuI7xb37eUQcW1GuBQ335PGIeCsivjrsnErvXUQsj4jXIuLJhn37RcT9EfFc8e+cFs89vzjnuYg4v8J8N0TEM8Xv7u6ImN3iuW1fB1NJrn1SXDvrTsm1T4pr2ymdZ7NPxiHXTsm9T4rr2ymjz5Ntn7TJZ6eMUa59UlzbThl/rqz6pLhetp1SeZ+klCbND9AH/AKYD8wAngA+Nuyci4G/KR6fDdxZUbaDgWOLx/sCzzbJ9ingH7t4/14ADmhzfAlwLxDA8cBDXfod/w+17+ju2r0DTgaOBZ5s2PctYGnxeClwfZPn7Qc8X/w7p3g8p6J8pwLTisfXN8s3mtfBVPnJuU+K62XdKb3QJw2/Zztl7Nnsk/G91rLslNz7ZDSvJTtlyLWy7ZM2+eyUsb/OsuyT4np2Snm/Z8coY882YX0y2d55chywIaX0fEppB3AHcMawc84Abi0e/wT4dETERAdLKb2SUnq0ePw28DQwb6KvW7IzgL9LNauB2RFxcMUZPg38IqX0YsXXHSKl9G/AG8N2N762bgXObPLURcD9KaU3UkpvAvcDi6vIl1K6L6X0XrG5Gjik7OtOMtn2CUyKTsmhT8BOGVc2+2Rcsu2USdAnYKcMyrlPWuWzU8Ys2z4BO6VEXe8TyLtTqu6TyTZ5Mg/4ZcP2y7z/f6iD5xQ3dSuwfyXpCsXb5o4BHmpy+ISIeCIi7o2Io6rMBSTgvoh4JCIubHJ8NPd3op0N3N7iWDfvHcAHU0qvQO3/NICDmpyTwz0EuIDabHozI70Opoqe6BPItlN6oU/ATimDfTI6PdEpmfYJ2Cmd6pU+ATtlNHqiT8BO6VCufQK90yml9sm00mLlodls6vCvExrNORMmIvYBfgp8NaX01rDDj1J7W9Y7EbEEWAEcUVU24KSU0qaIOAi4PyKeKWbzBnT73s0ATgeWNTnc7Xs3Wl29hwARcSXwHnBbi1NGeh1MFdn3CWTdKVn3CdgppVzcPhmL7Dsl4z4BO6UKOdxDO2V0su8TsFM6MQn6BLp/D0vvk8n2zpOXgUMbtg8BNrU6JyKmAbN4/9uQJkRETKdWILellO4afjyl9FZK6Z3i8T3A9Ig4oIpsxTU3Ff++BtxN7S2BjUZzfyfSacCjKaVXhx/o9r0rvDrwdr7i39eanNPVe1gs1PQ7wOdTSk3LaxSvg6ki6z4prpltp/RAn4Cd0hH7ZMyy7pSc+6S4pp3Smaz7BOyUMcq6T4pr2imdyblPIPNOmag+mWyTJw8DR0TE4cVs3dnAymHnrAQGVvr9DPBAqxtapuIzhrcAT6eUvt3inA8NfBYxIo6j9vt5faKzFdfbOyL2HXhMbaGdJ4edthI4L2qOB7YOvF2rIufQ4q1r3bx3DRpfW+cD/9DknFXAqRExJ2qrUp9a7JtwEbEY+CPg9JTSuy3OGc3rYKrItk8g707pkT4BO2Xc7JNxybZTcu6T4np2Suey7ROwU8Yh2z4BO6UkOfcJZNwpE9onqeJVgyf6h9rKyM9SW4H6ymLfNcXNA5gJ/BjYAPwXML+iXJ+k9jalnwOPFz9LgC8DXy7OuQRYR23F7NXAiRXet/nFdZ8oMgzcu8Z8AdxY3Nu1wMIK8+1FrRRmNezr2r2jVmavADupzap+idrnSP8FeK74d7/i3IXA9xuee0Hx+tsAfLHCfBuofe5w4PU3sAL7XOCedq+DqfqTa58U1862U3Lvk+L6dkpn2eyT8d3LLDsl5z5p91qyU1pmybZP2uSzU8Z+H7Psk+Ladkpn+bLpk+J62XZK1X0SxZMlSZIkSZLUxGT72I4kSZIkSVKpnDyRJEmSJElqw8kTSZIkSZKkNpw8kSRJkiRJasPJE0mSJEmSpDacPJEkSZIkSWrDyRNJkiRJkqQ2nDyRJEmSJElq4/8BZkVGy6s2Q+0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1368x342 with 4 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "n = 25\n",
    "x, y = create_data(n)\n",
    "X    = np.vander(x, 2)\n",
    "\n",
    "X0 = np.matrix(np.vander(x, 2)) # (n, 2)\n",
    "y0 = np.matrix(y).T             # (n, 1)\n",
    "\n",
    "w0 = lsq_solution_V0(X0, y0)\n",
    "ŷ0 = X0 * w0\n",
    "\n",
    "w1 = lsq_solution_V1(X, y)\n",
    "ŷ1 = X.dot(w1)\n",
    "\n",
    "w2 = lsq_solution_V2(X, y)\n",
    "ŷ2 = X.dot(w2)\n",
    "\n",
    "w3 = lsq_solution_V3(X, y)\n",
    "ŷ3 = X.dot(w3)\n",
    "\n",
    "fig, axs = plt.subplots(1, 4, figsize=(19, 4.75))\n",
    "for a in axs: a.plot(x, y, 'o')\n",
    "\n",
    "axs[0].set_title('Version 0\\nw =' + str(np.array(w0).T[0])) # convert (2, 1) to (2,)\n",
    "axs[0].plot(x, ŷ0, '-')\n",
    "\n",
    "axs[1].set_title('Version 1\\nw =' + str(w1))\n",
    "axs[1].plot(x, ŷ1, '-')\n",
    "\n",
    "axs[2].set_title('Version 2\\nw =' + str(w2))\n",
    "axs[2].plot(x, ŷ2, '-')\n",
    "\n",
    "axs[3].set_title('Version 3\\nw =' + str(w3))\n",
    "axs[3].plot(x, ŷ3, '-')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Performance comparison"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "         n        V0        V1        V2        V3\n",
      "       100     0.005     0.002     0.008     0.008\n",
      "      1000     0.004     0.002     0.009     0.008\n",
      "     10000     0.010     0.009     0.031     0.016\n",
      "    100000     0.050     0.047     0.344     0.088\n",
      "   1000000     0.557     0.550     1.616     1.111\n"
     ]
    }
   ],
   "source": [
    "import timeit\n",
    "\n",
    "print(9*' ' + 'n', end='')\n",
    "for v in [0, 1, 2, 3]: \n",
    "    print(8*' ' + 'V{}'.format(v), end='')\n",
    "print()\n",
    "\n",
    "for n in [100, 1_000, 10_000, 100_000, 1_000_000]:\n",
    "    print('{:10}'.format(n), end='')\n",
    "    for v in [0, 1, 2, 3]:\n",
    "\n",
    "        if v == 0:\n",
    "            setup = 'x, y = create_data(n); X = np.matrix(np.vander(x, 2)); y = np.matrix(y).T' \n",
    "        else: \n",
    "            setup = 'x, y = create_data(n); X = np.vander(x, 2)'\n",
    "        \n",
    "        t = timeit.timeit(\n",
    "            stmt    = 'lsq_solution_V{}(X, y)'.format(v), \n",
    "            setup   = setup, \n",
    "            number  = 100, \n",
    "            globals = globals()\n",
    "        )\n",
    "        \n",
    "        print('{:10.3f}'.format(t), end='')\n",
    "    print()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<table style=\"width:100%\">\n",
    "  <tr>\n",
    "      <td colspan=\"1\" style=\"text-align:left;background-color:#0071BD;color:white\">\n",
    "        <a rel=\"license\" href=\"http://creativecommons.org/licenses/by-nc/4.0/\">\n",
    "            <img alt=\"Creative Commons License\" style=\"border-width:0;float:left;padding-right:10pt\"\n",
    "                 src=\"https://i.creativecommons.org/l/by-nc/4.0/88x31.png\" />\n",
    "        </a>\n",
    "        &copy; C. Bauckhage and O. Cremers<br/>\n",
    "        Licensed under a \n",
    "        <a rel=\"license\" href=\"http://creativecommons.org/licenses/by-nc/4.0/\" style=\"color:white\">\n",
    "            CC BY-NC 4.0\n",
    "        </a>.\n",
    "      </td>\n",
    "      <td colspan=\"2\" style=\"text-align:left;background-color:#66A5D1\">\n",
    "          <b>Acknowledgments:</b>\n",
    "          This material was prepared within the project\n",
    "          <a href=\"http://www.b-it-center.de/b-it-programmes/teaching-material/p3ml/\" style=\"color:black\">\n",
    "              P3ML\n",
    "          </a> \n",
    "          which is funded by the Ministry of Education and Research of Germany (BMBF)\n",
    "          under grant number 01/S17064. The authors gratefully acknowledge this support.\n",
    "      </td>\n",
    "  </tr>\n",
    "</table>"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}