phys512/polynomial_fit.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "d6186669",
   "metadata": {},
   "source": [
    "# Polynomial fitting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "c9b00a3e",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f0c00302",
   "metadata": {},
   "source": [
    "Let's consider the example of fitting a polynomial to some measurements of $y_i(x_i)$ (ie. $y$ measured at a set of points $x_i$). An obvious way to solve the normal equations\n",
    "\n",
    "$$\\mathbf{A^T} \\mathbf{d} = \\mathbf{A^T}\\mathbf{A}\\mathbf{a}$$ \n",
    "\n",
    "is to calculate the inverse of $\\mathbf{A^T}\\mathbf{A}$ and write\n",
    "\n",
    "$$\\mathbf{a} = (\\mathbf{A^T}\\mathbf{A})^{-1} \\mathbf{A^T} \\mathbf{d}.$$\n",
    "\n",
    "For simplicity, let's take a constant error $\\sigma_i$ for each point, so the $\\sigma_i$'s drop out of the equations and we'll just ignore them. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "5473e41a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Polynomial coeffs =  [-0.06810457 -0.91610815  0.89798733 -0.29154606  0.26086863  0.93971567\n",
      "  0.24332917 -0.44333463 -0.39585528  0.974724    0.94526284  0.89132267\n",
      "  0.58953405  0.61203039 -0.01513055  0.06471341 -0.17141171 -0.9525837\n",
      "  0.97697325  0.870457    0.30180934 -0.3201125  -0.95648576 -0.84369759\n",
      " -0.30176928  0.95358891  0.25481124 -0.17775    -0.42285418 -0.83911108\n",
      "  0.15120805]\n",
      "Fitted coeffs =  [-4.32457335e+09 -1.67263774e+09  5.02037577e+09  1.62012140e+09\n",
      " -2.39298194e+09 -5.89549081e+08  5.35570374e+08  7.41558693e+07\n",
      " -2.05309226e+06  2.19871111e+06 -3.29462755e+07  3.11737081e+06\n",
      "  1.03325380e+07 -2.25757037e+06 -2.39449501e+06  6.85293308e+05\n",
      "  6.23022355e+05 -2.52531237e+05 -1.55019220e+05  8.31031818e+04\n",
      "  2.93902694e+04 -1.08301980e+04 -4.34265523e+03 -2.27619395e+03\n",
      "  7.54748993e+02  1.06871481e+03 -1.52381191e+02 -1.53604104e+02\n",
      "  1.92676516e+01  7.22158292e+00 -8.90406676e-01]\n",
      "frac error =  [ 6.34990204e+10  1.82580816e+09  5.59069779e+09 -5.55699986e+09\n",
      " -9.17313022e+09 -6.27369641e+08  2.20101180e+09 -1.67268390e+08\n",
      "  5.18647078e+06  2.25572585e+06 -3.48540902e+07  3.49746511e+06\n",
      "  1.75266168e+07 -3.68865829e+06  1.58255669e+08  1.05896643e+07\n",
      " -3.63465547e+06  2.65100362e+05 -1.58673941e+05  9.54697496e+04\n",
      "  9.73792518e+04  3.38314747e+04  4.53921941e+03  2.69687893e+03\n",
      " -2.50207958e+03  1.11972907e+03 -5.99015965e+02  8.63158088e+02\n",
      " -4.65657119e+01 -9.60622992e+00 -6.88861959e+00]\n",
      "rms deviation / max(y) = 4.58505\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiIAAAGsCAYAAADg5swfAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABFKklEQVR4nO3dd3wUdeLG8c9ueoCEkgABQglSlN4NKsJZwI4ooqLAidjv9M4G3u/O84pY8PQO8cRTwYoCgh0UQUClSIsUCdJ7AqEkgYS0nd8fA4EIhGzI7nd293m/XvvaYZndfYaQ7JOZ73zHZVmWhYiIiIgBbtMBREREJHSpiIiIiIgxKiIiIiJijIqIiIiIGKMiIiIiIsaoiIiIiIgxKiIiIiJijIqIiIiIGKMiIiIiIsaoiIiIiIgxAVNE5s+fzzXXXEODBg1wuVx8/PHHXr/G5MmT6dixI7GxsTRp0oTnn3++6oOKiIhIhQVMETl8+DAdOnRg3LhxlXr+jBkzGDx4MPfccw+rV6/mlVde4cUXX+Tll1+u4qQiIiJSUa5AvOidy+Vi+vTp9O/fv/SxgoIC/vSnPzFp0iQOHjxI27ZtefbZZ+nduzcAt956K0VFRUyZMqX0OWPHjuW5555j27ZtuFwuP2+FiIiIBMwekTN54IEHWLhwIR988AErV65k4MCB9OvXj/Xr1wN2UYmOji7znJiYGHbs2MHWrVtNRBYREQl5QVFEtm3bxoQJE5gyZQoXXXQRzZs355FHHuHCCy9kwoQJAPTt25dp06Yxe/ZsPB4Pv/zyCy+88AIAu3fvNhlfREQkZIWbDlAVVq1aRUlJCS1btizzeEFBAXXq1AFgxIgRbNy4kauvvpqioiLi4uJ48MEH+etf/4rbHRR9TEREJOAERRE5dOgQYWFhLFu2jLCwsDJ/V716dcAeV/Lss8/y9NNPk5GRQWJiIrNnzwYgJSXF75lFREQkSIpIp06dKCkpYc+ePVx00UXlrhsWFkbDhg0BmDRpEqmpqSQmJvojpoiIiPxKwBSRQ4cOsWHDhtI/b968mbS0NGrXrk3Lli0ZPHgwQ4YM4YUXXqBTp07s3buX2bNn0759e6666iqysrKYOnUqvXv35siRI6VjSubNm2dwq0REREJbwJy+O3fuXPr06XPS40OHDmXixIkUFRXxj3/8g7fffpudO3eSkJDA+eefz1NPPUW7du3IysrimmuuYdWqVViWRWpqKv/85z/p0aOHga0RERERCKAiIiIiIsFHp4uIiIiIMSoiIiIiYoyjB6t6PB527dpFjRo1NAW7iIhIgLAsi9zcXBo0aHDGubocXUR27dpFcnKy6RgiIiJSCdu3b6dRo0blruPoIlKjRg3A3pC4uDjDaURERKQicnJySE5OLv0cL4+ji8ixwzFxcXEqIiIiIgGmIsMqNFhVREREjFEREREREWNURERERMQYFRERERExRkVEREREjFEREREREWNURERERMQYFRERERExRkVEREREjFEREREREWNURERERMQYFRERERExJmSLyO7sfBZszGJ3dr7pKCIiIiHL0Vff9ZXpP6xi3YxxJLOH20qGM3pAOwZ1a2w6loiISMgJuT0iu7Pz+dvnPzMyfBKDw2dT08rhiWmrtWdERETEgJArIpuzDnPAqs4vnoYAdHH/QollsSUrz3AyERGR0BNyRaRZQjXcLljqaQVAV/c6wlwumibEGk4mIiISekKuiCTFxzB6QDuWWa0B6Ob+hacHtCUpPsZwMhERkdATckUEYFC3xjx+91AAOoVvYVCHBMOJREREQlNIFhGAusmtoEYDXJ4i2L7IdBwREZGQFLJFBJcLzvkNhMfAga2m04iIiISkkJxHpNQlf4UrX4CIaNNJREREQlJoF5HqiaYTiIiIhLTQPTTzawWHTCcQEREJOSoiWxfCuB4w+XbTSUREREJOaB+aAaheF/amQ9Z6OJwF1XQqr4iIiL9oj0id5pDUAawSWPup6TQiIiIhRUUEoM0AALKXfKCL34mIiPiRigjwmScVj+UiPnMxg599jw+XbDMdSUREJCSEfBHZnZ3PgzOy+NbTEYBb3LN5Ytpq7RkRERHxg5AvIpuzDuOx4N2SSwEYGDaPcKuALVl5hpOJiIgEv5A/a6ZZQjXcLpjn6cC0kgv5oLgPxa4omibEmo4mIiIS9EJ+j0hSfAyjB7TD5Qrjj0X3sYzzeHpAW5LiY0xHExERCXohv0cEYFC3xvRqmciWrDyaJsTaJaTwMERWMx1NREQkqIX8HpFjkuJjSG1eh6RqYTDnH/BSO8jZZTqWiIhIUFMR+TV3OGyaC3n7YMZjYFmmE4mIiAQtFZFfc7vh6hfBFQZrP2PjzLE6lVdERMRHVEROpX470lo9CECjRU/x0LPjNMmZiIiID6iInMLu7HwG/NSFmSXdiHIV83rEGCZPn649IyIiIlVMReQU7EnOXDxYdD8LSs6jhiufiRGj2bFzp+loIiIiQUVF5BSOTXJWQCR3Fj3CNyWdeLb4ViJr1GHBxiztGREREakiLsty7mkhOTk5xMfHk52dTVxcnF/f+8Ml23hi2mpKLItwl8V1nRoxfcVOPBZcEfYjTySvIbnX7dD8NxBVo2re1LKg+AgU5oHLZZ/Bc+wWFmE/JiIi4nDefH6riJRjd3Y+W7LyiI10c/0rC/BYABZfRD5BG/fWo2u5ILEV1G4OcQ2gRn3o9cjxF1k6AfZtOF4wCg9BUZ49YZo7HIZ9fnzdt66BzfNPHcblhj9lQHiU/ecf/wd70yE+GWo2hlpNIbE1RGpqehERMcubz2/NrFqOpPgYkuJjWLAx62gJAXDxcNG9XB/2HUNrriL60Da7EOxNt/86pnbZIrJqKmz9/tRvEB5d9s8R5ZSIsKjjJQTgl69gw6yy67jcdiGq3xaue0WlREREHE9FpAKOjRk5VkbSrcY8V3IbdG/FGzMWcZ5rM8muLG5pHcZ5SXHszs5nc9ZhmiVUI6lNf2jY2S4dkdXK3iKq2Ydjjh1yuX48uMMg/Oh1bjzF9s0qgaIjZUN1vh2SOsDBbZC9HbLWQ14W7Dt6H3HCtXLmPw8eD7Tsaz9Hh3hERMQh/HJoZty4cTz//PNkZGTQoUMHxo4dS/fu3c/4PNOHZk504piRMJeLx/q14tmZ6SfsKcF+/IpWPDvDftztgtED2tGrZeLxYnL0YnplykpVXWAvNxMyV0H+QWh3o/2YZcG/zoPco9PV12hgF5K2A6DJhfYEbiIiIlXIUWNEPvzwQ4YMGcKrr75Kjx49eOmll5gyZQrr1q2jbt265T7XSUUEjo8ZaZoQy+asw9z6v8UnreNylZ0V3nX0sROLCcCoaau8LiunKi9nXMdTwsEf3qDkl6+ptft73MXHz/gpiK1PQZe7iLvk4cq9tg/XMf3+yuisdUy/vzIqoz9fOxg4qoj06NGDbt268fLLLwPg8XhITk7md7/7HSNHjiz3uU4rIifanZ3PBc/MKbNHxA14zvA8N8AJh3mO/vGMZeX6Tg1Lz9qp7DrRrkLGX5hHo4zZJG6bQZwrj9eKr2Jd+8eOPs8iwuXhHwM6+uT9K7qOL19bGZ2VMVi2w/T7K2NwZPTrHnQfc0wRKSwsJDY2lqlTp9K/f//Sx4cOHcrBgwf55JNPyqxfUFBAQUFB6Z9zcnJITk52ZBGBih2uqYxTlZWqWufYYxFWIX3caaRbyWyxkgBIda/h+YjxTCzux/ueS8izok7xqv7L6IvXVkZnZazIOqbfXxmV0V8ZXVR+D7rTOOasmaysLEpKSqhXr16Zx+vVq0d6evpJ648ePZqnnnrKl5Gq1KBujenVMrH0cE1SfAw1YyNKy4kbsI7ejqnIf1APv3pSFa5z7LECIpnpKTtOZ2DYPBq5svi/iHe52/qU14qv5t2SS8nnV2f3+CmjL167qtYx/f4VWcf0+1fVOqbfvyLrmH7/iqxj+v0rso7p96/IOr58bQv70H4X1zrqug7y3fRFWC4X3V012EUddlt1GPnRqpPKyqBujct/M4fz6R6RXbt20bBhQxYsWEBqamrp44899hjz5s1j8eKyYywCbY/I6Zw4lmT+L3vL7DV5ekBbgLMuK75o81EUcn3Y99wf9gnJ7r0AZFlxvFR8A5NKfkMJYcYzOmUd0+8fLBkrso7p91dGZfRNRoum7gx6ulaT4MrmpeIbS/9uVuSjtHCffEmRAiucNVZTbit8gryjvyCGuVx8P7IPgKP2knizR8Snp0wkJCQQFhZGZmZmmcczMzOpX7/+SetHRUURFxdX5haIkuJjSG1eh6T4GAZ1a8z3I/swacT5fD+yD4O6NS7z2A+jfsMzN7Qj7OgptWEuF6NvaMfoAWUfu6FzwypZp7znFRDJFM8lLLzqK35s/ze2WvVIcOXwj4gJfFnn347I6JR1TL9/sGQMlu0w/f7KGDgZz3Ht4LHwyays/ThzIh/mHxET+G3YTNxY2GvDSiuFJZ6WLPa0ZomnFRs9SeRbkUS5igmnpLSEALRiMxO+28QFz8zh1v8t5oJn5gTc1eL9Mli1e/fujB07FgCPx0Pjxo154IEHAnqwalU7cS/KiQOUTnysqtap8PMO5JK/8A2arvw37r7/YHfKAOdlNLiO6fcPlozBsh2m318ZHZ5x9xwK5z5PZMYKSrkjKGjQjcyEVCJ6PcT8jQfL2YPuoalrDzEUsNayD8XUJJcFUb/nF6sR/1f0W1ZbKcDxvSQm94w4ZrAq2KfvDh06lPHjx9O9e3deeuklJk+eTHp6+kljR34tlIqIo+UfhKg4SuccSf8CcnZB1+FoHhIRkQqY8w97ckl3OJxzGbQfCC36QlT1MqudqdCceLj/fHc6b0WPIcqTR4nl4rnimxlfcjXgYtKI80ltXsfAhtocVUQAXn755dIJzTp27Mh//vMfevToccbnqYg4UP5BGNcdDmVC8vnQ/xWo09x0KhER5ygpghXvQL12kNzNfiw3E5ZNhK53QPXEs3r5E4uJ+1AGS169h6vDFgHwZUl3Hi6+nzkj+wLmxo04rohUloqIA3k8sPQN+Oav9gX8IqrBVWOgwy2aOl5EZMNs+PJR2L/R/mXtjpk+/9n44Y9bWfXJv/lL+EQiXSVk1OnB913/w2OfbjR2do1jBqtKEHK7ofsIuG8RNL0Iig7Dx/fCtBFwJMd0OhERM3IzYeod8O4Au4TEJkCb68E60zSXZ29Q9ybc/9g/WX/5O3giYqm/bzF7v/hn6Vk7HsseZ7I7O7/8FzJERUQqp2YyDPkEfvN/4AqDVVNg/EVweJ/pZCIi/mNZ8NMH8HI3WP2RfRX0HvfC71fA+feAO+zMr1EFkuJjaHPBVbhvn87eptfw7+Lry/x9iWWxJSvPL1m8pSIilecOg16P2rse4xvbe0hia5tOJSLiP2s/g+l3Q0E2NOgEI+bAFc9AtKHhBI3Pp/j6/1HoKjszdpjLRdOEWDOZzkBFRM5ecne4ex5c9cLxY6GFefZ4EhGRYNb6KmjWC/r8Hwz/xi4jhiXFxzB6QDvCXRZ3hn1BI9d+nh7Q1hETnZ2KBqtK1fOUwPs3QUQsDPgfRJw8RbyISECyLFgzDVpfDeFH9zp4Svx2CMYbhz5/gupLx3Gk2WVED5ni1xMKNFhVzNqdBpvnw9pP7YFb+QdNJxIROXtF+fbg/Kl3wMwTJuR0YAkBqN59CLgjiN48C9Z+yu7sfBZszHLcoFUVEal6DbvAbR/Zk6Bt/QEmXGlPgCYiEqgOZ8HEq+GnSfaA1Nop9t4RJ6vbGi76IwC5n/+J3s987chp4FVExDea9YLffgnV68OeNfBGXziwxXQqERHv7dsIb1wGO5dCdE24fTr0/F1gzJ10wYOUxCZSI287g9yzAeedzqsiIr5Tvx0M/xpqN4fsbTDhKti/yXQqEZGK27EM3rjc/tkV3xiGz4KU3qZTVVxkNba0+z0Avwv/mCgKAWedzqsiIr5Vq4m9ZyShJeQfgEN7TScSEamYwjyYdDPkZUFSB7jzG0hsaTqV12J7DGOnlUCiK5sBYd8BzjqdN9x0AAkBNerDsC/sQzPJ3U2nERGpmMhYuOF/sOi/cMMbJ12gLlAk1Y5jeZvhbF39ORs8DUuv7OuU03l1+q6YsXcdRNWAuAamk4iIlFV0JPimHbAsduccOenKvr6i03fF2XatgDf7wTvXQ95+02lERI7bNBde7mr/shRMXC6S4mNIbV7HMXtCjlEREf+LqWVPBLQ3Hd67EQoOmU4kIgI7lsKkWyF7OywcZzqNb2TvhLnP2vcOoSIi/lerqX36W0wt2LkMPrrTnplQRMSUjNXw7g32FcVT+sCVz5tO5BvT74a5T0Pae6aTlFIRETPqngu3ToawKPhlBnz9Z9OJRCRUHdhqHyo+chAadYeb3zs+fXuw6Xirfb9qqmMmZFMREXOSu8P1r9rLi8bBktfN5hGR0JN/0L421uE9UK8tDJ4MkdVMp/Kd1ldBWCRkrYM9Pzti2ncVETGr7QD4zdG9Iaum6hCNiPjX7L/Z49VqNIDBU+xDxsEsOh5aXA7Az19P4IJn5hif9l1FRMy76GG4+iV73IhDLx4lIkHq0ifh3Gvg1g9DZzqBtgMAiF3/KZ6jh2dMTvuuIiLmuVzQ9bcQccIpZR6PuTwiEjqi42HQu5DU3nQS/2nRF487kqbuTJq7jl+Q1NS07yoi4iyeEpj9d5h+l2MGUolIkNkwGxaPD92fMVHVKWqUSr4VSQvX8dN4TU37rinexVky18APL4GnGBp2gfPvNZ1IRILJga3w0XD72leR1aDTbaYTGRE1YBxT0/OZ9ckGwDI67bumeBfnWfRfmDkS3OEw9HNokmo6kYgEg6J8+0q6GSuhQWe4Y2bwnqZbQbuz830y7bumeJfA1uMeaDfQ3ivy0Z32by4iImfDsuDzP9olJLYODHon5EsIYE/73qym0WnfVUTEeVwuuPpFqJ0COTvgswdD91iuiFSN5W/DT++Dyw03ToD4RqYTmffzp/BKKsz6i9EYKiLiTFE14IbX7cMzP38CK94xnUhEAtXeX2DG4/byJX+BlIvN5nGKem3sSSUv+7vRGCoi4lwNu9iTnUXE2oVERKQyti+CkgL7GjI9HzSdxjnqNIekDuA2WwX0012crefvoU1/+0J5IiKV0XkIJLS0f44Y/tCVk6mIiLO53WVLiKdEs6+KiPcan286gZyGqqEEjs3fwbjusCvNdBIRcbpDe+Ht/rB3nekkcgYqIhI4lrwO+zbAJ/dDcaHpNCLiZF8+DJu+hY/v01l3DqciIoHjyjH2+f+Zq2H+86bTiIhTrZ5mn23nCoOrXrCnBBDHUhGRwFE90S4jAN+9oEM0InKyQ3vhy0fs5V6PQIOORuPImamISGBpOwDOuw6sEvj0ASgpNp1IRJzCsuCLP0LePqjXDi56xHQiqQAVEQk8V74A0TUhYxX8ON50GhFxijXTYO2n9rxD/V+B8EjTiaQCVEQk8FRPhMuespd3LjObRUScI22SfX/RI5DU3mwWqTDNIyKBqdMQqNkEUnqbTiIiTnHLB/b1ZNrfbDqJeEFFRAKT2w3N+5hOISJOEhZuz6IqAUWHZiTwHd4HX/0JCg6ZTiIi/lZSBAtfgaJ800mkklREJLBZFrw7ABa+DN//y3QaEfG3BWPhq1Hw1rWauCxAqYhIYHO5oNej9vKCl2H/ZrN5RMR/DmyFec/Zy13v0MRlAUpFRAJf66vsQaslBTDrz6bTiIi/fP1/UJwPTS+CDhqgGqhURCTwuVzQd7Q9nfPaz2DTPNOJRMTXNs+35wxxueGKZ7U3JICpiEhwqHcedBtuL88cpRlXRYJZSTHMGGkvdx0O9dqYzSNnRUVEgkfvURBTC/asgeVvmU4jIr6yfKL9fR5dE/o8YTqNnCXNIyLBI7Y2XPIXyPwZ2lxvOo2I+EpKH2h1lT2XUGxt02nkLLksy7nnO+Xk5BAfH092djZxcXGm44iIiJNYlsaGOJQ3n986NCPBrbjQdAIRqSqekrJ/VgkJCioiEpz2pMPb/eGLP5hOIiJVZfIQ+Ph+yNltOolUIRURCU4FubDpW0h7H/asNZ1GRM7W1gWQ/jn8NAkKdTmHYKIiIsEpuRucey1YHvjmr6bTiMjZsCyY9Rd7ufMQSGhhNo9UKRURCV6XPGlPcvbLTNj8nek0IlJZaz+DHUsgIhZ6jzSdRqqYiogEr4RzoMswe3nO33VBLJFAVFIMs5+yl1Pvhxr1zeaRKqciIsHt4scgPBq2L4YNs02nERFvrXgb9m2A2DrQ8/em04gPqIhIcKtRH7rdaS8vm2A2i4h4x7Jg8Wv2cq/HIFrzSQUjzawqwe+Ch6B2CnS63XQSEfGGywV3zISlb0DXO0ynER/x2R6Rf/7zn/Ts2ZPY2Fhq1qzpq7cRObPqifYF8cIjTScREW/F1ISLHtb3bxDzWREpLCxk4MCB3Hvvvb56CxHvlRTBga2mU4jImRzYogHmIcJnh2aeesoe5Txx4kRfvYWId3b/BJOHQkQM3PMDuDVESsSR8g/C+F6Q0BIGvQc16plOJD7kqDEiBQUFFBQUlP45JyfHYBoJOjWbQN5+KMiGNdOg3Y2mE4nIqSz6LxzJtmdIrpZgOo34mKN+JRw9ejTx8fGlt+TkZNORJJjE1ISeD9jL858Hj8doHBE5hbz9sOgVe7n3SHCHmc0jPudVERk5ciQul6vcW3p6eqXDjBo1iuzs7NLb9u3bK/1aIqfU426Iioe96fZ1K0TEWRaOg4IcqNsGzr3OdBrxA68OzTz88MMMGzas3HVSUlIqHSYqKoqoqKhKP1/kjKLjofsI+G6MfTv3Gl1KXMQpDu+Dxa/ay31GaRxXiPCqiCQmJpKYmOirLCL+cf599q7f3T/Zs622uNR0IhEBWDTOvrJu/fbQ+mrTacRPfFY3t23bRlpaGtu2baOkpIS0tDTS0tI4dEiXbxbDqtWBLr+1lzd8YzaLiNgsCzbNs5cvfkx7KkOIy7J8c6L2sGHDeOutt056/Ntvv6V3794Veo2cnBzi4+PJzs4mLk5T+0oVys2A/ZuhSarpJCJyjKcE1s+CFpfrsEyA8+bz22dFpCqoiIiIiAQebz6/VTlFDu+Dg9tMpxAJXZlroOiI6RRiiIqIhLbVH8FLbWHmKNNJREJTcSG8NxD+3R4yVplOIwaoiEhoq9sGivIg/QvYt9F0GpHQs/IDyNkJuOwp3SXkqIhIaKvbGlr0BSxY+LLpNCKhxVMC379oL/f8HYRrHqlQpCIicsHv7fu09+FwltksIqFkzXTYvwliakGXYabTiCEqIiJNLoAGnaD4CCx53XQakdBgWcf3hvS4F6Kqm80jxqiIiLhc9m5hgB9fg8I8s3lEQsHGOZC5GiKqQY+7TKcRg1RERMC+uFbNxnAkB3b8aDqNSPDbttC+7zzEPjQjIcura82IBK2wcLj+NajVBOIamE4jEvx+83/QZgDE1DSdRAxTERE5RtO9i/hXvfNMJxAH0KEZkVPRnCIivnF4H+TsMp1CHERFRORElgWTboWxnWG7xoqIVLkF/4GX2sP3L5lOIg6hIiJyIpfr+DHrxa8ajSISdApyYdkE8BRBQgvTacQhVEREfq3H3fb9z59oF7JIVVrxLhzJhtrNoeUVptOIQ6iIiPxaUgdo3BM8xbD0TdNpRIJDSTEsesVeTr0f3Pr4EZv+J4icyvn32PdLJ+jy5CJVYe2ncHAbxNaBDreYTiMOoiIiciqtroL4ZMjLgtUfmU4jEvgWj7fvuw6HyFizWcRRVERETiUsHLrdaS+v/NBsFpFAl70Tdv8E7nDoNtx0GnEYTWgmcjqdh0C1BGh7o+kkIoEtviH88WfYtghq1DedRhxGRUTkdGJrQ6fbTKcQCQ6xtaH1laZTiAPp0IxIRXg8uiqvSGUczjKdQBxORUTkTNK/hHHdYN4zppOIBJaSInj1Qnj9Mjiw1XQacSgVEZEzsmDfBlj+jk7lFfHG2s8gdzcc2KyxIXJaKiIiZ9KyH8Q1gvz98PPHptOIBI7SU3bvgPAos1nEsVRERM7EHQZdh9nLS94wGkUkYOxKg+2L7FN2u95hOo04mIqISEV0GmL/QN3xI+xeaTqNiPP9+Jp9f15/HZaRcqmIiFREjXpw7jX28lLtFREp16G9sGqKvdzjHrNZxPFUREQqqnSm1SlwJMdsFhEnW/kBlBRCg87QqKvpNOJwmtBMpKKaXADn3wetroSoGqbTiDhXj3ugZmOIrA4ul+k04nAqIiIV5XJBv9GmU4g4X1gEnHed6RQSIHRoRkREqo5lmU4gAUZFRMRb2Ttg5ij4/A+mk4g4y8Ht8O8OMO95FRKpMBUREW8dzoJFr9gzrR7aYzqNiHMsfwsOboXN8zQ2RCpMRUTEWw06QsOu4CmC5W+bTiPiDCUnfD90G242iwQUFRGRyjj2g3b52/aVeUVCXfoXcCgTqteD1lebTiMBREVEpDLO6w9R8Ud3Q881nUbEvGMT/XW63T5rRqSCVEREKiMyFtrfZC8ve8tsFhHTstbD5vngckOXYabTSIBRERGprC5D7fv0L+wBrCKhaukE+77F5VAz2WwWCTia0Eyksuq3s3/wJrQES+NEJIS1uxGOZEPbAaaTSABSERE5G4OnmE4gYl7DzvZNpBJ0aEZERESMUREROVseD2z8Fn78n+kkIv61Zy189hDsSjOdRAKYDs2InK3dafBOfwiPto+Vx9QynUjEP5a/DcsmwOG9cPN7ptNIgNIeEZGz1aAT1G0DxUdg5WTTaUT8o7gAfppkL3ceajaLBDQVEZGz5XIdP5V32Vu62JeEhvTPIf8A1GgA51xiOo0EMBURkarQ/ib70MyeNbBzmek0Ir537LoynW4Dd5jZLBLQVEREqkJMLTjvOnt52USjUUR87sAW2DQXcNlFROQsqIiIVJVjx8lXT4OCXLNZRHxpxbv2fUpvqNXEaBQJfCoiIlWlSU+ocw7EN4SD20ynEfGdmNpQIwk6DzGdRIKAy7KcO7IuJyeH+Ph4srOziYuLMx1H5MxyM6F6XXsAq0gwKym278M0C4SczJvPb/0PEqlKNeqZTiDiHyogUkV0aEbEFwrzIGOV6RQiVevQHlg34/jeEJEqoCIiUtV2pcELreC9m8BTYjqNSNVJew8m3QyTNTZEqo6KiEhVq3uuPa9C7i77GjQiwcCyjs8d0vpKs1kkqKiIiFS18Chod5O9vOIds1lEqsrWH2D/JoisAef1N51GgoiKiIgvHJvkad2XkLffbBaRqnBsb0i7GyCqutksElRURER8Iak91G8HJYWwaqrpNCJnJ/8A/PyJvay5Q6SK+ayIbNmyheHDh9OsWTNiYmJo3rw5Tz75JIWFhb56SxFn6Xh0r0jau2ZziJyt1dPsq0vXbQMNOptOI0HGZ0UkPT0dj8fD+PHjWbNmDS+++CKvvvoqTzzxhK/eUsRZ2g0EdwTs/gn2bTSdRqTytv9o33e8VZP1SZXz68yqzz//PP/973/ZtGlThdbXzKoS8FZNhYadoXaK6SQiZ2f3SohvBLG1TSeRAODYmVWzs7OpXfv0/4kLCgooKCgo/XNOTo4/Yon4TrsbTScQqRpJ7U0nkCDlt8GqGzZsYOzYsdx9992nXWf06NHEx8eX3pKTk/0VT8T3PB7TCUS84ymBgkOmU0iQ87qIjBw5EpfLVe4tPT29zHN27txJv379GDhwICNGjDjta48aNYrs7OzS2/bt273fIhGn2bMWPhgMU3S2gQSYzfNgTEuYMdJ0EgliXh+aefjhhxk2bFi566SkHD8evmvXLvr06UPPnj157bXXyn1eVFQUUVFR3kYScTZ3OKR/Dq4wyM2AGvVNJxKpmLT3oegweIpMJ5Eg5nURSUxMJDExsULr7ty5kz59+tClSxcmTJiA261pSyQEJbSA5B6wfTH89AFc+JDpRCJndiQb1n5uL3e41WwWCWo+awY7d+6kd+/eNG7cmDFjxrB3714yMjLIyMjw1VuKOFfHwfZ92nv2NTtEnG7Nx1CcDwmt7DO/RHzEZ2fNzJo1iw0bNrBhwwYaNWpU5u/8eMawiDO0uR5mPA5Zv8Cu5dCwi+lEIuX7aZJ93/EWzR0iPuWzPSLDhg3DsqxT3kRCTnQcnHu1vZw2yWwWkTPZtxG2LQSXG9oPMp1GgpwGbYj4S4db7PvVU6FYlzoQB/vpA/s+pQ/ENTCbRYKeXyc0EwlpKb3hnMvgnEvBKjGdRuT0Og+B8Eio38F0EgkBfp3i3Vua4l1ERCTwePP5rUMzIiIiYoyKiIi/HcmGZW/B2s9MJxEpq+CQPQvw6mn29O4ifqAiIuJvaZPgs9/Ddy+YTiJS1tpP7VmA5/zDPmNGxA/0P03E39rdaE/7vmsF7Ek/8/oi/pL2vn2vuUPEj1RERPytWgK0uNxe/klziohDHNgKW74DXND+ZtNpJISoiIiY0OHoD/qVH+pYvDjDsblDmvWCmslms0hIURERMaFlP4iuCbm77Uuti5hkWfDTscMyg81mkZCjIiJiQniUPVYEjv8mKmLKtoVwYAtEVj9+KQIRP1ERETGlwy32mQmFh3VFXjHL8kCj7tCmP0RWM51GQoxmVhUxxbLgUCbUqG86iYitpAjCIkynkCCgmVVFAoHLpRIizqISIgaoiIg4QfZOOLzPdAoJRWs/h7z9plNICFMRETFt1pPwYhtY+qbpJBJqsnfAh7fBC61VhMUYFRER0xJbAZY9uZlzh2xJMFr5IWBBo25QrY7pNBKiVERETDv3WoiIhf0bYccS02kkVFhW2SndRQxRERExLaq6XUbg+AeDiK/tWAr7Ntgl+LzrTKeREKYiIuIEx34jXTMNio6YzSKhIe09+/7cayGqhtksEtJUREScoGkviGsER7Lhlxmm00iwK8qH1dPs5Y63ms0iIU9FRMQJ3G5of5O9vGqq2SwS/HYuh8JDEJ8MTS8ynUZCXLjpACJyVOfboXaKjteL7zW9AB5Oh/2b7RIsYpCKiIhT1E6xbyL+UL2ufRMxTFVYRCSUFOaZTiBShoqIiNMseR3G94KM1aaTSLCxLHj9EphwFWRtMJ1GBFAREXGeTXNh90/2TKsiVWnXCtjzM+xcCtUTTacRAVRERJynw9HTKVdOhpJis1kkuBybMK/11RAdbzaLyFEqIiJO0+IyiK0Dh/fAxjmm00iwKC6A1UdPDdeU7uIgKiIiThMWAe0G2ss/acp3qSK/zIT8A1AjCVL6mE4jUkpFRMSJOhz9jTX9S/vDQ+RspR0dc9R+ELjDzGYROYGKiIgTJXWAuudBSQGsmW46jQS6Q3tg/df2sqZ0F4fRhGYiTuRyQechsPUHSGxtOo0Euqga0P+/9tkyia1MpxEpw2VZlmU6xOnk5OQQHx9PdnY2cXFxpuOIiIhIBXjz+a1DMyIiImKMioiI0+3bCN+/CB6P6SQSiBaOgx/+bY8TEXEgjRERcbLiAnitNxTkQMOu0EyXbBcvlBTBd/+CvCx7rFHLvqYTiZxEe0REnCw8Cs67zl7+6QOzWSTwrJ9ll5BqdaH5JabTiJySioiI0x073fLnj6HwsNEoEmDS3rPv298EYdoBLs6kIiLidI1ToVZTKDwEaz83nUYCxeF98MtX9rLmDhEHUxERcTqX6/hMq7oir1TU6qngKbInx6vXxnQakdNSEREJBO0H2feb5kL2TqNRJEAcOyzTcbDZHCJnoCIiEghqN4PGPSEqDvasNZ1GnK7oCNQ5x/7/cuwCiiIOpdFLIoHi+v9C9XoQEWM6iThdRDTc+CYU5ev/izieiohIoKjV1HQCCTQqIRIAdGhGJNBYFhzcbjqFOFXmz7An3XQKkQpTEREJJAe2wstdYXwvKC40nUacaO7T8EoPWDDWdBKRClEREQkk8Y2g4BDk74f1X5lOI05zeB+sm2kvp/Qxm0WkglRERAKJO8yeJRMgTXOKyK+cOHdI/bam04hUiIqISKA5Nkvm+q/s34BFjtHcIRKAVEREAk3dcyGpI3iK7d+ARQAyVsPun8AdAW1vNJ1GpMJUREQC0bG9Imnvm80hznFs+v9W/aBaHbNZRLygIiISiNreAO5w2J2mmVbFPqV7wzf2sg7LSIDRhGYigahaAlz2d/swTUJL02nENJcL7v4O1n8N51xqOo2IV1RERAJV6n2mE4iThEfCuVebTiHiNR2aEREJZMWF4PGYTiFSaT4tItdeey2NGzcmOjqapKQkbr/9dnbt2uXLtxQJLQe2wswnYMbjppOIKUvfhP90hOXvmE4iUik+LSJ9+vRh8uTJrFu3jo8++oiNGzdy4406rUykyuRlwaJxsGwiHMk2nUZMSHsPDm61r7QrEoB8OkbkD3/4Q+lykyZNGDlyJP3796eoqIiIiAhfvrVIaGjQGRLPhb1rYfVH0PUO04nEnzJWQ8ZKe+6QdvolTwKT38aI7N+/n/fee4+ePXuetoQUFBSQk5NT5iYi5XC5oNNt9rJ2zYee0rlDroDY2maziFSSz4vI448/TrVq1ahTpw7btm3jk08+Oe26o0ePJj4+vvSWnJzs63giga/9IHtOkV3LIXON6TTiL8WF8NMH9vKxCe5EApDXRWTkyJG4XK5yb+np6aXrP/roo6xYsYKvv/6asLAwhgwZgmVZp3ztUaNGkZ2dXXrbvn175bdMJFRUT7R/IwZY8a7ZLOI/v8y0xwhVrw/nXGY6jUiluazTtYLT2Lt3L/v2lX+hrZSUFCIjI096fMeOHSQnJ7NgwQJSU1PP+F45OTnEx8eTnZ1NXFycNzFFQssvX8H7N0FMbXh4nT2nhAS3d2+EDbPgwj/ApX81nUakDG8+v70erJqYmEhiYmKlgnmOnuteUFBQqeeLyGk0vwTqtYOUi6EoT0UkFPR6BKrXhU63m04icla83iNSUYsXL2bJkiVceOGF1KpVi40bN/LnP/+ZzMxM1qxZQ1RU1BlfQ3tERLxgWfbgVRERw7z5/PbZYNXY2FimTZvGJZdcQqtWrRg+fDjt27dn3rx5FSohIuIllRARCUA+m0ekXbt2zJkzx1cvLyKn4vHAlvlQmAetrzSdRnxhyw/2nDFdhkFSe9NpRM6aLnonEkzWTIOPhkOtZtCyH7h1Oamgs/QNu4gAXP0vs1lEqoB+SokEk1ZXQGQNOLAZti0wnUaqWt5+WPu5vdxZg1QlOKiIiASTyGrQdoC9rDlFgs+qKVBSYJ8hldTRdBqRKqEiIhJsjp3OueZjXQgvmFjW8Wn8Ow/R4GQJGioiIsGmUVdIaAXF+bB6muk0UlV2p0HmKgiLgvYDTacRqTIqIiLBxuU6Pn5ghS6EFzSWv23fn3sNxNQym0WkCqmIiASj9jfbF8I7kgP5B02nkaoQn2xfV0aDVCXI+Gxm1aqgmVVFzsLeXyChhcYSBJOSYnC5dVq2OJ5PrzUjIgEisaXpBFLVwvQjW4KParVIsCvMg/2bTKeQytq/CdbNBE+J6SQiPqEiIhLMNsyGf7WGaXebTiKVtXg8TBoEnz1oOomIT6iIiASzem2h8DDs+BEy15hOI94qzIOfJtnLbfobjSLiKyoiIsGsRj1odfTid0snmM0i3lsz3Z6UrmYTSPmN6TQiPqEiIhLsuv7Wvl/5of0btgSOZUfLY5dhOlNGgpb+Z4sEu2a9oVZTKMixr84rgWH3StixBNwR0Ok202lEfEZFRCTYud3Qeai9rMMzgePY3pBzr4bqdc1mEfEhFRGRUNDpNnum1Z1LdSpvILAs2LvOXu56h9ksIj6m2XFEQkH1unDdK/YF8WqnmE4jZ+JywbAv7AvdJXU0nUbEp1REREJFh0GmE4g3XC5o0Ml0ChGf06EZkVBUUmw6gZxObqZ9sUKREKEiIhJK9m+GDwbDm5fb4xDEeeb8DV5oDcvfMZ1ExC9URERCSXQ8bPgGdi6zTw0VZ8nbD6umQtFh+8rJIiFARUQklMTWhrY32ss/vmY2i5xs+VtQfATqt4fkHqbTiPiFiohIqOk+wr5f87E9HkGcoaQYlrxhL/e42x6sKhICVEREQk2DjvZv254iWDbRdBo55peZkL0dYmpD2xtMpxHxGxURkVDU/S77fumbUFxoNovYfhxv33cZChExZrOI+JGKiEgoOvdaqF4PDmXA2k9Np5GcXbDlB3C5NZOqhBxNaCYSisIjofdIKCmCFpebTiNxDeChlXYZqdnYdBoRv1IREQlV+s3bWeIbafZbCUk6NCMiYlLhYdMJRIxSEREJZR6PPYPnG5fbk2mJf3k8ML4XvHsDHNxmOo2IESoiIqHM5bInNtu+GJZNMJ0m9Kz/CvZtgO1L7NN2RUKQiohIKHO5IPV+e3nxazqV198WvGzfdxkKUdXNZhExREVEJNS1GQA1kuxTeVdPNZ0mdOxaAVu/B3c49LjHdBoRY1REREJdeOTxCc4WjtNVef3l2N6QNgMgvqHZLCIGqYiICHT9LURUg8zVsGmu6TTB7+B2WDPdXu75gNksIoapiIgIxNSCTrfZywtfNpslFCybAFYJNL0IkjqYTiNilCY0ExHb+fdA1i8ar+APvR6F+GRIaGE6iYhxLsty7gHhnJwc4uPjyc7OJi4uznQcERERqQBvPr91aEZExF88JfYkZiJSSkVERMo6nAWz/w7zx5hOEnxWTYVx3WHlFNNJRBxDRUREytq+GL4bA9+/BPkHTacJHh4PfP8v2LceDm41nUbEMVRERKSslldA4rlQmAtL/mc6TfBY9yXsTYeoOOh2p+k0Io6hIiIiZbndcNHD9vKi/+rqsFXBsuC7F+zl7iMgpqbROCJOoiIiIidrcz3Uagp5+2D526bTBL5Nc2HXcgiPgR73mk4j4igqIiJysrBwuOAhe/mH/0BxgdE4Ae/Y3pAuQ6F6otksIg6jIiIip9bxVvtieLm74KcPTKcJXBmrYMt39sXtev7OdBoRx9HMqiJyauFR9liRzNXQrJfpNIGrXlsY+rldSOIbmU4j4jgqIiJyet1HmE4Q+FwuaHaRfRORk+jQjIiIrxzJNp1AxPFURETkzDJ/hg8Gw4p3TScJHNsWwQvnwrdPm04i4mgqIiJyZhtnQ/rnMPdZKC40nSYwfPs0FB2G3N2mk4g4moqIiJxZ1+FQvT5kb4MV75hO43xbvofN88AdAb0eNZ1GxNFURETkzCJjj8+2Ov95KMo3m8fJLOv44ZjOQ6BmY7N5RBxORUREKqbLUIhrZB9qWDrBdBrn2jQXtv4AYVHQ6xHTaUQcT0VERComPAouPnqY4ft/QUGu2TxO5PHArL/Yy13vgLgGZvOIBAAVERGpuI6DoXYKHN4LS143ncZ5Mlcdv8KuxoaIVIhfikhBQQEdO3bE5XKRlpbmj7cUEV8Ii4C+T8Nlf9PF204lqQM8sBQG/A+q1TGdRiQg+KWIPPbYYzRooF2UIkGh1RVwwYMQEW06iTPVagKt+plOIRIwfF5EZsyYwddff82YMWN8/VYi4m8lxXA4y3QK8/L2w45lplOIBCSfFpHMzExGjBjBO++8Q2xs7BnXLygoICcnp8xNRBxq53J49QKYpuvRMO85eP03MPvvppOIBByfFRHLshg2bBj33HMPXbt2rdBzRo8eTXx8fOktOTnZV/FE5GzF1IJ9G2HjHFg303QaczJ/hh9fs5ebXmA2i0gA8rqIjBw5EpfLVe4tPT2dsWPHkpuby6hRoyr82qNGjSI7O7v0tn37dm/jiYi/1G4GqffZyzMfh6IjZvOYYFkw4zGwSqD11dD8N6YTiQQcl2VZljdP2Lt3L/v27St3nZSUFG666SY+++wzXC5X6eMlJSWEhYUxePBg3nrrrTO+V05ODvHx8WRnZxMXF+dNTBHxh4JD8HI3yN0FvZ+A3o+bTuRfq6fB1N9CeDTc/6M9UFVEvPr89rqIVNS2bdvKjPHYtWsXffv2ZerUqfTo0YNGjRqd8TVUREQCwOqPYOodRz+MF0OtpqYT+UfBIRjXHXJ2hmYJEymHN5/f4b4K0bhx2esrVK9eHYDmzZtXqISISIBoMwCWTYTN82HmKLhlkulE/vHdGLuE1GwMF/zedBqRgKWZVUXk7LhccOUYcIdDUV7oXBCvfnuITYB+z0BEjOk0IgHLZ4dmqoIOzYgEkJ3LoUEnu5iEioJciKphOoWI43jz+a09IiJSNRp2Do0SUlJ8fFklROSsqYiISNU6kg2fPQgbZptOUvX2b4KxnewBus7dmSwSUFRERKRqLRhrD1797EH70EWwKCmG6ffAwW2wdILpNCJBQ0VERKrWhX+Amk0gezvMHGk6TdX54UXYvhgia8B140LjMJSIH6iIiEjViqwG/f8LuGDFu/akX4Fu53KY+4y9fNUYTVwmUoVURESk6jW9AHo9Yi9/9hAc2Go0zlkpOATT7gJPMZzXH9oPMp1IJKioiIiIb1w8Ehp1h4Js+OhOKCkynch7lgWf/R72rYcaDeDqF3VIRqSKqYiIiG+EhcMNr0NUHBzcag/yDDSeEnvSMnc4DJwAsbVNJxIJOprQTER8a/N8qNMC4pJMJ6m8fRuhTnPTKUQChiOuNSMiAkCzXmX/XFJs7y1xstwMiKkF4VH2n1VCRHxGh2ZExH/WTLevWOvkwav5B+Ht6+Cd6yFvv+k0IkFPRURE/KO4EOY9B/s3wrsD4PA+04lOVlwIk2+Hven2LKpFeaYTiQQ9FRER8Y/wSBg8FeIawb4N8P5AKDxsOtVxJcUw/W57TEtkdbh1MsQ3Mp1KJOipiIiI/8Q3hNun2eMvdi6D9wc5o4x4SuCT+2DNNHBHwMCJkNTedCqRkKAiIiL+ldgKbp1iT5W+5Tt490az16QpKYJP7oeVH4IrzD5Nt8Vl5vKIhBgVERHxv+RucPt0iIqHbQvsEmBKbgasnwUutz3vybnXmMsiEoIcfg6diASt5G4w9BP7ejRdh5vLUTMZBk+BQ5nQ6gpzOURClIqIiJjToJN9O+ZIDmycDW2u9+37bp5vn6Z73rX2nxt29u37ichpqYiIiDNYFnzxR1g1BX75CvqNtge1VqWifJg7GhaMhfBoe7xKYquqfQ8R8YrGiIiIM1gWJLS0x2r8NAnGdoW09+3Hq+K1138D/70Afvg3WB57r0t88tm/toicFV1rRkScZetC+OxByFpn/7leW7joj3Bef3CHefdalgWb59kTqW39wX6sRpJ9FV2NBxHxGW8+v7VHREScpUkq3PM9XPpXe2KxzNUw9Q57L4m3CnJh0i12CQmLgtQH4L5FKiEiDqIxIiLiPOGRcOEfoPNQ+PF/sPQNaHXl8b//4d+w5mOo1QSi4sAdDkcOQs4u+/bgT+ByQXSc/RqeIvv1NFOqiOPo0IyIOJ+npOxhmbf7w6ZvT7/+H9dCXAOfxxKRU/Pm81t7RETE+X49NuTK5+0L0x3cZk8R7ymG6HioVhfqnWffi0hAUBERkcCT0MK+iUjA02BVERERMUZFRERERIxRERERERFjVERERETEGBURERERMUZFRERERIxRERERERFjVERERETEGBURERERMUZFRERERIxRERERERFjVERERETEGBURERERMcbRV9+1LAuAnJwcw0lERESkoo59bh/7HC+Po4tIbm4uAMnJyYaTiIiIiLdyc3OJj48vdx2XVZG6YojH42HXrl3UqFEDl8tVpa+dk5NDcnIy27dvJy4urkpf2wm0fYEv2Lcx2LcPgn8btX2Bz1fbaFkWubm5NGjQALe7/FEgjt4j4na7adSokU/fIy4uLmj/g4G2LxgE+zYG+/ZB8G+jti/w+WIbz7Qn5BgNVhURERFjVERERETEmJAtIlFRUTz55JNERUWZjuIT2r7AF+zbGOzbB8G/jdq+wOeEbXT0YFUREREJbiG7R0RERETMUxERERERY1RERERExBgVERERETEmJIrIli1bGD58OM2aNSMmJobmzZvz5JNPUlhYWO7zjhw5wv3330+dOnWoXr06N9xwA5mZmX5K7b1//vOf9OzZk9jYWGrWrFmh5wwbNgyXy1Xm1q9fP98GraTKbJ9lWfzlL38hKSmJmJgYLr30UtavX+/boGdh//79DB48mLi4OGrWrMnw4cM5dOhQuc/p3bv3SV/De+65x0+Jyzdu3DiaNm1KdHQ0PXr04Mcffyx3/SlTptC6dWuio6Np164dX375pZ+SVp432zhx4sSTvlbR0dF+TOud+fPnc80119CgQQNcLhcff/zxGZ8zd+5cOnfuTFRUFOeccw4TJ070ec7K8nb75s6de9LXz+VykZGR4Z/AXho9ejTdunWjRo0a1K1bl/79+7Nu3bozPs/f34chUUTS09PxeDyMHz+eNWvW8OKLL/Lqq6/yxBNPlPu8P/zhD3z22WdMmTKFefPmsWvXLgYMGOCn1N4rLCxk4MCB3HvvvV49r1+/fuzevbv0NmnSJB8lPDuV2b7nnnuO//znP7z66qssXryYatWq0bdvX44cOeLDpJU3ePBg1qxZw6xZs/j888+ZP38+d9111xmfN2LEiDJfw+eee84Pacv34Ycf8sc//pEnn3yS5cuX06FDB/r27cuePXtOuf6CBQu45ZZbGD58OCtWrKB///7079+f1atX+zl5xXm7jWDPYHni12rr1q1+TOydw4cP06FDB8aNG1eh9Tdv3sxVV11Fnz59SEtL46GHHuLOO+/kq6++8nHSyvF2+45Zt25dma9h3bp1fZTw7MybN4/777+fRYsWMWvWLIqKirj88ss5fPjwaZ9j5PvQClHPPfec1axZs9P+/cGDB62IiAhrypQppY+tXbvWAqyFCxf6I2KlTZgwwYqPj6/QukOHDrWuu+46n+apahXdPo/HY9WvX996/vnnSx87ePCgFRUVZU2aNMmHCSvn559/tgBryZIlpY/NmDHDcrlc1s6dO0/7vIsvvth68MEH/ZDQO927d7fuv//+0j+XlJRYDRo0sEaPHn3K9W+66SbrqquuKvNYjx49rLvvvtunOc+Gt9vozfem0wDW9OnTy13nscces9q0aVPmsUGDBll9+/b1YbKqUZHt+/bbby3AOnDggF8yVbU9e/ZYgDVv3rzTrmPi+zAk9oicSnZ2NrVr1z7t3y9btoyioiIuvfTS0sdat25N48aNWbhwoT8i+s3cuXOpW7curVq14t5772Xfvn2mI1WJzZs3k5GRUeZrGB8fT48ePRz5NVy4cCE1a9aka9eupY9deumluN1uFi9eXO5z33vvPRISEmjbti2jRo0iLy/P13HLVVhYyLJly8r827vdbi699NLT/tsvXLiwzPoAffv2deTXCiq3jQCHDh2iSZMmJCcnc91117FmzRp/xPWLQPsaVlbHjh1JSkrisssu44cffjAdp8Kys7MByv3sM/E1dPRF73xlw4YNjB07ljFjxpx2nYyMDCIjI08ai1CvXj3HHg+sjH79+jFgwACaNWvGxo0beeKJJ7jiiitYuHAhYWFhpuOdlWNfp3r16pV53Klfw4yMjJN28YaHh1O7du1y89566600adKEBg0asHLlSh5//HHWrVvHtGnTfB35tLKysigpKTnlv316evopn5ORkREwXyuo3Da2atWKN998k/bt25Odnc2YMWPo2bMna9as8fkFPv3hdF/DnJwc8vPziYmJMZSsaiQlJfHqq6/StWtXCgoKeP311+nduzeLFy+mc+fOpuOVy+Px8NBDD3HBBRfQtm3b065n4vswoPeIjBw58pQDh068/foHws6dO+nXrx8DBw5kxIgRhpJXXGW20Rs333wz1157Le3ataN///58/vnnLFmyhLlz51bdRpTD19vnBL7exrvuuou+ffvSrl07Bg8ezNtvv8306dPZuHFjFW6FVIXU1FSGDBlCx44dufjii5k2bRqJiYmMHz/edDSpgFatWnH33XfTpUsXevbsyZtvvknPnj158cUXTUc7o/vvv5/Vq1fzwQcfmI5ykoDeI/Lwww8zbNiwctdJSUkpXd61axd9+vShZ8+evPbaa+U+r379+hQWFnLw4MEye0UyMzOpX7/+2cT2irfbeLZSUlJISEhgw4YNXHLJJVX2uqfjy+079nXKzMwkKSmp9PHMzEw6duxYqdesjIpuY/369U8a5FhcXMz+/fu9+j/Xo0cPwN7z17x5c6/zVoWEhATCwsJOOsusvO+f+vXre7W+aZXZxl+LiIigU6dObNiwwRcR/e50X8O4uLiA3xtyOt27d+f77783HaNcDzzwQOng9zPteTPxfRjQRSQxMZHExMQKrbtz50769OlDly5dmDBhAm53+TuDunTpQkREBLNnz+aGG24A7JHS27ZtIzU19ayzV5Q321gVduzYwb59+8p8cPuSL7evWbNm1K9fn9mzZ5cWj5ycHBYvXuz1mUVno6LbmJqaysGDB1m2bBldunQBYM6cOXg8ntJyURFpaWkAfvsankpkZCRdunRh9uzZ9O/fH7B3Dc+ePZsHHnjglM9JTU1l9uzZPPTQQ6WPzZo1y6/fb96ozDb+WklJCatWreLKK6/0YVL/SU1NPelUTyd/DatCWlqa0e+18liWxe9+9zumT5/O3Llzadas2RmfY+T70GfDYB1kx44d1jnnnGNdcskl1o4dO6zdu3eX3k5cp1WrVtbixYtLH7vnnnusxo0bW3PmzLGWLl1qpaamWqmpqSY2oUK2bt1qrVixwnrqqaes6tWrWytWrLBWrFhh5ebmlq7TqlUra9q0aZZlWVZubq71yCOPWAsXLrQ2b95sffPNN1bnzp2tFi1aWEeOHDG1Gafl7fZZlmU988wzVs2aNa1PPvnEWrlypXXddddZzZo1s/Lz801swhn169fP6tSpk7V48WLr+++/t1q0aGHdcsstpX//6/+nGzZssP72t79ZS5cutTZv3mx98sknVkpKitWrVy9Tm1Dqgw8+sKKioqyJEydaP//8s3XXXXdZNWvWtDIyMizLsqzbb7/dGjlyZOn6P/zwgxUeHm6NGTPGWrt2rfXkk09aERER1qpVq0xtwhl5u41PPfWU9dVXX1kbN260li1bZt18881WdHS0tWbNGlObUK7c3NzS7zPA+te//mWtWLHC2rp1q2VZljVy5Ejr9ttvL11/06ZNVmxsrPXoo49aa9eutcaNG2eFhYVZM2fONLUJ5fJ2+1588UXr448/ttavX2+tWrXKevDBBy2322198803pjahXPfee68VHx9vzZ07t8znXl5eXuk6Tvg+DIkiMmHCBAs45e2YzZs3W4D17bfflj6Wn59v3XfffVatWrWs2NhY6/rrry9TXpxm6NChp9zGE7cJsCZMmGBZlmXl5eVZl19+uZWYmGhFRERYTZo0sUaMGFH6Q9RpvN0+y7JP4f3zn/9s1atXz4qKirIuueQSa926df4PX0H79u2zbrnlFqt69epWXFyc9dvf/rZM0fr1/9Nt27ZZvXr1smrXrm1FRUVZ55xzjvXoo49a2dnZhragrLFjx1qNGze2IiMjre7du1uLFi0q/buLL77YGjp0aJn1J0+ebLVs2dKKjIy02rRpY33xxRd+Tuw9b7bxoYceKl23Xr161pVXXmktX77cQOqKOXa66q9vx7Zp6NCh1sUXX3zSczp27GhFRkZaKSkpZb4fncbb7Xv22Wet5s2bW9HR0Vbt2rWt3r17W3PmzDETvgJO97l34tfECd+HrqNhRURERPwuoM+aERERkcCmIiIiIiLGqIiIiIiIMSoiIiIiYoyKiIiIiBijIiIiIiLGqIiIiIiIMSoiIiIiYoyKiIiIiBijIiIiIiLGqIiIiIiIMSoiIiIiYsz/A/Ffp+a/ogrYAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# first generate some data from a random polynomial\n",
    "\n",
    "k = 30   # order of polynomial\n",
    "n = 100   # number of data points\n",
    "\n",
    "x = np.linspace(-2,2,n)\n",
    "a0 = -1.0 + 2*np.random.rand(k+1)\n",
    "print(\"Polynomial coeffs = \", a0)\n",
    "poly = np.polynomial.Polynomial(a0)\n",
    "y = poly(x) + np.random.normal(scale=0.0, size=n)\n",
    "plt.plot(x,y,'.')\n",
    "\n",
    "# Compute the design matrix\n",
    "# A_{ij} = f_j(x_i)\n",
    "A = np.zeros((n, k+1))\n",
    "for kk in range(k+1):\n",
    "    A[:, kk] = x**kk\n",
    "\n",
    "# Now do the linear algebra part \n",
    "rhs = np.transpose(A)@y\n",
    "lhs = np.transpose(A)@A\n",
    "a = np.linalg.inv(lhs)@rhs\n",
    "\n",
    "print(\"Fitted coeffs = \", a)\n",
    "print(\"frac error = \", (a-a0)/a0)\n",
    "\n",
    "xx = np.linspace(-2,2,1000)\n",
    "poly = np.polynomial.Polynomial(a)\n",
    "yy = poly(xx)\n",
    "plt.plot(xx,yy,'--')\n",
    "\n",
    "dev = np.sqrt(np.mean((y-poly(x))**2))\n",
    "print('rms deviation / max(y) = %g' % (dev/np.max(np.abs(y)),))\n",
    "\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "13fd394a",
   "metadata": {},
   "source": [
    "If you run this with no noise (`scale=0` on line 10), you should see that for low polynomial order the coefficients are fit to machine precision, but the errors become large pretty quickly as you increase the polynomial order. Something is going horribly wrong with the matrix inversion as the polynomial order increases!"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1aaa27d6",
   "metadata": {},
   "source": [
    "## Using SVD"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cac127b3",
   "metadata": {},
   "source": [
    "Singular value decomposition (SVD) is extremely useful in these situations where you are dealing with a close-to-singular matrix. The ratio of the largest to smallest singular values is known as the **condition number** of the matrix, and measures how close to singular it is."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "e8d9067e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Range of singular values = 1.16147e+14\n",
      "np.linalg.cond gives 1.16147e+14\n"
     ]
    }
   ],
   "source": [
    "U, Sdiag, VT = np.linalg.svd(A,0)\n",
    "\n",
    "print(\"Range of singular values = %g\" % (max(abs(Sdiag))/min(abs(Sdiag))))\n",
    "\n",
    "print(\"np.linalg.cond gives %g\" % (np.linalg.cond(A)))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "762f72b0",
   "metadata": {},
   "source": [
    "Let's use SVD to rewrite the normal equations:\n",
    "\n",
    "$$\\mathbf{A^T} \\mathbf{d} = \\mathbf{A^T}\\mathbf{A}\\mathbf{a}$$ \n",
    "\n",
    "$$ \\mathbf{VSU^T} \\mathbf{d} = \\mathbf{VSU^T}\\mathbf{USV^T} \\mathbf{a}$$ \n",
    "\n",
    "$$ \\mathbf{VSU^T} \\mathbf{d} = \\mathbf{VS^2V^T} \\mathbf{a}$$ \n",
    "\n",
    "$$ \\mathbf{V^TVS}\\mathbf{U^T} \\mathbf{d} = \\mathbf{V^T}\\mathbf{VS^2V^T} \\mathbf{a}$$ \n",
    "\n",
    "$$ \\mathbf{SU^T} \\mathbf{d} = \\mathbf{S^2V^T} \\mathbf{a}$$ \n",
    "\n",
    "$$ \\mathbf{U^T} \\mathbf{d} = \\mathbf{SV^T} \\mathbf{a}$$ \n",
    "\n",
    "$$\\Rightarrow \\mathbf{VS^{-1}}\\mathbf{U^T} \\mathbf{d} =  \\mathbf{a}$$ \n",
    "\n",
    "The matrix $\\mathbf{VS^{-1}}\\mathbf{U^T}$ is called the **pseudo-inverse**. We can use it to map from the data vector $\\mathbf{d}$ to the model parameters $\\mathbf{a}$.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "1595e75e",
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Fitted coeffs =  [-0.06810442 -0.91610793  0.89798456 -0.29154781  0.2609005   0.93970728\n",
      "  0.24313354 -0.44324493 -0.3951416   0.97437286  0.94372559  0.89199829\n",
      "  0.59155273  0.61116028 -0.01757812  0.06556702 -0.17004395 -0.95316315\n",
      "  0.97619629  0.87065506  0.30215454 -0.32018661 -0.95658875 -0.84368467\n",
      " -0.30175591  0.95358771  0.25480914 -0.17774992 -0.42285404 -0.83911109\n",
      "  0.15120804]\n",
      "frac error =  [-2.11980533e-06 -2.44679418e-07 -3.07735337e-06  5.99586571e-06\n",
      "  1.22147670e-04 -8.92834498e-06 -8.03952611e-04 -2.02320286e-04\n",
      " -1.80288818e-03 -3.60245783e-04 -1.62627498e-03  7.58003564e-04\n",
      "  3.42420505e-03 -1.42168454e-03  1.61763938e-01  1.31906354e-02\n",
      " -7.97943866e-03  6.08285240e-04 -7.95278498e-04  2.27540718e-04\n",
      "  1.14377625e-03  2.31537019e-04  1.07666287e-04 -1.53112627e-05\n",
      " -4.43319898e-05 -1.25424394e-06 -8.24456806e-06 -4.50115938e-07\n",
      " -3.31363628e-07  4.82297004e-09 -4.12189271e-08]\n",
      "rms deviation / max(y) = 3.81606e-08\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiIAAAGsCAYAAADg5swfAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA420lEQVR4nO3deXxU1f3/8fedyb4NSwIkrAEEVJBNwKBFqVRQqmD4KrV8K/qluPzAShcVaJViW4Pot9pavtXWiq21YrUB6la1IlgVEVmUoKAgyJKwBGESQjJJZs7vjyFjAtll5mZmXs/HYxzm5ty5n5OZybw9995zLWOMEQAAgA0cdhcAAACiF0EEAADYhiACAABsQxABAAC2IYgAAADbEEQAAIBtCCIAAMA2BBEAAGAbgggAALANQQQAANgmbILIW2+9pSuvvFJZWVmyLEsrVqxo8XO8+uqruuCCC5SamqqMjAxNmTJFu3fvPuO1AgCA5gmbIFJWVqbBgwdryZIlrVp/165dmjRpkr75zW9q8+bNevXVV1VcXKzc3NwzXCkAAGguKxwvemdZlpYvX67JkycHlnk8Hv30pz/VM888o2PHjmngwIG6//77dckll0iSnn/+eV133XXyeDxyOPz564UXXtCkSZPk8XgUGxtrQ08AAIhuYTMi0pTZs2dr7dq1WrZsmT766CNdc801mjBhgj777DNJ0vDhw+VwOLR06VJ5vV653W499dRTGjduHCEEAACbRMSIyJ49e9S7d2/t2bNHWVlZgXbjxo3TyJEjdd9990mS1qxZo2uvvVZHjhyR1+tVTk6OXn75ZbVr186GXgAAgIgYEdmyZYu8Xq/69eunlJSUwG3NmjXauXOnJOnAgQOaOXOmpk+frvXr12vNmjWKi4vTf/3XfykMsxgAABEhxu4CzoTjx4/L6XRqw4YNcjqddX6WkpIiSVqyZIlcLpcWL14c+Nlf//pXde/eXevWrdMFF1wQ0poBAECEBJGhQ4fK6/Xq0KFD+sY3vlFvmxMnTgQOUq1RE1p8Pl/QawQAAKcLm10zx48f1+bNm7V582ZJ/tNxN2/erD179qhfv36aNm2arr/+euXn52vXrl16//33lZeXp5deekmSNHHiRK1fv1733nuvPvvsM23cuFE33nijevbsqaFDh9rYMwAAolfYHKy6evVqjR079rTl06dP15NPPqmqqir98pe/1F/+8hft379f6enpuuCCC7Rw4UINGjRIkrRs2TItXrxYn376qZKSkpSTk6P7779fAwYMCHV3AACAwiiIAACAyBM2u2YAAEDkIYgAAADbtOmzZnw+nwoLC5WamirLsuwuBwAANIMxRqWlpcrKyjrtjNVTtekgUlhYqO7du9tdBgAAaIW9e/eqW7dujbZp00EkNTVVkr8jaWlpNlcDAACao6SkRN27dw98jzemTQeRmt0xaWlpBBEAAMJMcw6r4GBVAABgG4IIAACwDUEEAADYhiACAABsQxABAAC2IYgAAADbEEQAAIBtCCIAAMA2BBEAAGAbgggAALBN1AaRIne53t1ZrCJ3ud2lAAAQtdr0tWaC5dn1ezQvf4t8RnJYUl7uIE0d0cPusgAAiDpRNyJS5C4PhBBJ8hlpfn4BIyMAANgg6oLIruIyna3d+knMs7rW+aYkyWuMdhefsLkyAACiT9QFkez0ZPVz7NPsmJW60rFWkuS0LPVKT7K5MgAAok/UBZFMV6KmXNBPkpRkeeS0LN2XO1CZrkSbKwMAIPpE5cGqF53dQ9ogDejo1Ns3jCWEAABgk6gMIopLkSQly6NkQggAALaJul0zkqS4k8eDVJbZWwcAAFEuOoNI7MkgUsWZMgAA2Ck6g0hcsv++skwyxt5aAACIYtF5jEhSunTTaik22e5KAACIatEZRJwxUtZQu6sAACDqReeuGQAA0CZE54iIJK17TCo7LI2YKaV2trsaAACiUvQGkXd/J7n3SP0uJ4gAAGCT6N01UzOXSBVziQAAYJfoDSI1c4lUMpcIAAB2id4gUjOXCCMiAADYJnqDSCzTvAMAYLfoDSJx7JoBAMBuURxE2DUDAIDdovf03Yt+JJ3/P5Kru92VAAAQtaI3iHTsY3cFAABEvejdNQMAAGwXvSMiB7ZIO1dJ7XtJ50yyuxoAAKJS9I6I7PtAev0e6cNldlcCAEDUit4gUnPWDPOIAABgG4JIFfOIAABgl+gNIsysCgCA7aI3iLBrBgAA20VvEKkZEWHXDAAAtoneIMK1ZgAAsF30ziOSmiVd/08pLsXuSgAAiFrRG0RiE6TeF9tdBQAAUS16d80AAADbRXcQ2fw36Z3fSmVH7K4EAICoFL27ZiTpjV9IpYVS9jek5I52VwMAQNSJ7hGR+JMHqnqO21sHAABRKsqDSKr/vpIgAgCAHUIWRBYtWiTLsjRnzpxQbbJpNafuekrtrQMAgCgVkiCyfv16PfbYYzrvvPNCsbnmqxkRIYgAAGCLoAeR48ePa9q0afrjH/+o9u3bB3tzLcOuGQAAbBX0IDJr1ixNnDhR48aNa7Ktx+NRSUlJnVtQxXGwKgAAdgrq6bvLli3Txo0btX79+ma1z8vL08KFC4NZUl0jZkgDrpA69AndNgEAQEDQRkT27t2r22+/XU8//bQSEhKatc68efPkdrsDt7179warPL9OZ0t9vim17xnc7QAAgHpZxhgTjCdesWKFrr76ajmdzsAyr9cry7LkcDjk8Xjq/Kw+JSUlcrlccrvdSktLC0aZAADgDGvJ93fQds1ceuml2rJlS51lN954owYMGKC77rqryRASEkd3S5+vlhLbS+dMsrsaAACiTtCCSGpqqgYOHFhnWXJysjp27Hjactsc2CK9cLvUbSRBBAAAG0T3zKo1Z81w+i4AALYI6UXvVq9eHcrNNS0woRlBBAAAO0T3iEhgQjNmVgUAwA7RHURqX2smOCcPAQCARkR3EIk/GUR81VK1x95aAACIQtEdRGpGRCQOWAUAwAYhPVi1zXE4palPS7GJdUMJAAAIiegOIpJ09rftrgAAgKgV3btmAACArRgR2blKOrZXyh4jdci2uxoAAKIKIyL/+bX0wg+k/RvsrgQAgKhDEAnMrsqkZgAAhBpBJMHlv/eU2FsHAABRiCASn+a/r3DbWwcAAFGIIFIzIkIQAQAg5AgiBBEAAGxDEElg1wwAAHZhHpHsMdI1f5badbe7EgAAog5BpH0v/w0AAIQcu2YAAIBtCCJV5dLWFdKmv9pdCQAAUYddM1Xl0nPT/f8+7zuSk18JAAChwohIzYRmErOrAgAQYgQRZ4wUl+L/d8UxW0sBACDaEEQkJjUDAMAmBBGJ680AAGATgohUa0SEY0QAAAglgojErhkAAGzCuaqSNHq2NOS7UtdhdlcCAEBUIYhI/uvNAACAkGPXDAAAsA0jIpJ09AupcJOU1FHK/obd1QAAEDUYEZGkz9/0T/O+dondlQAAEFUIIlKts2aO2VoGAADRhiAiSYnt/fcnvrS3DgAAogxBRJISO/jvywkiAACEEkFEkpJqgshRyRh7awEAIIoQRKSvRkR81Tpw+LC9tQAAEEUIIpKe/bBYFSZWknTtQy/p2fV7bK4IAIDoEPXziBS5yzUvf4ves76vCsXpiEnV/PwCjemXoUxXot3lAQAQ0aI+iOwqLpPPSMtNrYnMjNHu4hMEEQAAgizqd81kpyfLYdVd5rQs9UpPsqcgAACiSNQHkUxXovJyB+lsx15d4XhPZzv26b7cgYyGAAAQAlG/a0aSpo7ooW/vKVDylqdUesFPlDqih90lAQAQFaJ+RKRGsitDkpTqK7W5EgAAogdBpEYSs6sCABBqBJEaibVmVwUAACFBEKlRMyLChe8AAAgZgkiNmivwsmsGAICQIYjUqNk1c4JdMwAAhAqn79ZIy5K+/bCU1NHuSgAAiBoEkRrxKdL5N9pdBQAAUYVdMwAAwDYEkdr2b5AK8iX3PrsrAQAgKhBEanvtHun5G6U979ldCQAAUYEgUluKf5p3lR22tw4AAKIEQaS25E7+++OH7K0DAIAoQRCpLZkREQAAQokgUltg10yxvXUAABAlCCK1BUZE2DUDAEAoEERqqzlGhF0zAACEBDOr1taxj/Tth6TULLsrAQAgKhBEakvqIJ3/P3ZXAQBA1GDXDAAAsA1B5FR710sF/5BKiuyuBACAiBfUIJKXl6cRI0YoNTVVnTp10uTJk7V9+/ZgbvLr+9dc6fn/8V93BgAABFVQg8iaNWs0a9Ysvffee3r99ddVVVWlyy67TGVlZcHc7NeTwpkzAACESlAPVv3Xv/5V5/GTTz6pTp06acOGDRozZkwwN916NXOJHD9obx0AAESBkJ4143a7JUkdOnSo9+cej0cejyfwuKSkJCR11ZF28tTdksLQbxsAgCgTsoNVfT6f5syZowsvvFADBw6st01eXp5cLlfg1r1791CV95XUTP99KQerAgAQbCELIrNmzVJBQYGWLVvWYJt58+bJ7XYHbnv37g1VeV+pCSKcNQMAQNCFZNfM7Nmz9eKLL+qtt95St27dGmwXHx+v+Pj4UJTUsLSaERF2zQAAEGxBDSLGGN12221avny5Vq9erezs7GBu7sxo30ua+GspravdlQAAEPGCGkRmzZqlv/3tb1q5cqVSU1N14MABSZLL5VJiYmIwN9168anSiBl2VwEAQFSwjDEmaE9uWfUuX7p0qW644YYm1y8pKZHL5ZLb7VZaWtoZrg4AAARDS76/g75rJiwVfSQVfyplDpbSz7K7GgAAIhbXmqnP2w9J/5ghffaa3ZUAABDRCCL1YVIzAABCgiBSn8CkZgfsrQMAgAhHEKlPahf/PbOrAgAQVASR+rhOTrrmtmFmVwAAoghBpD7tevjv3fslb7W9tQAAEMEIIvVJ6SI54yTjZap3AACCKCTXmgk7Dod01SNSYgcpqaPd1QAAELEIIg0Z/B27KwAAIOKxa6YRRe5yvbuzWEXucrtLAQAgIjEi0oB/vvW+Vr26QiUmUavNMOXlDtLUET3sLgsAgIjCiEg9itzlWvXqCj0cu0QznS/LZ6T5+QWMjAAAcIYRROqxq7hMe3wZkqRu1mFJktcY7S4+YWdZAABEHIJIPbLTk1UofxDJtI7IKa+clqVe6Uk2VwYAQGQhiNQj05WoH159kTwmVjGWT12to7ovd6AyXYl2lwYAQEQhiDRg6shecnTsJUlafl1nDlQFACAICCKNiM3oJ0nqWME1ZwAACAaCSGM69vHfF39mbx0AAEQo5hFpzODrpF4XSZ3PtbsSAAAiEkGkMZ3P8d8AAEBQsGsGAADYhiDSlG0vS2sekI5xwCoAAGcau2aa8tYDUuFGKaO/1K673dUAABBRGBFpSvpZ/vsjnDkDAMCZRhBpSkZ///2hbfbWAQBABCKINKXzQP/9wa321gEAQAQiiDSl08nTd4u3S9WV9tYCAECEIYg0xdVNindJvmqp+FO7qwEAIKIQRJpiWV/NrHqwwN5aAACIMJy+2xwT8qTYJKlDb7srAQAgohBEmiNriN0VAAAQkdg1AwAAbEMQaa73/yg9dyNTvQMAcAYRRJpr89PS1nxp3/t2VwIAQMQgiDRX1/P99/s+sLcOAAAiCEGkubr5g0jJzrUqcpfbXAwAAJGBINJML32ZJUmKP1SgSxa9pmfX77G5IgAAwh9BpBmK3OW67bUSHTPJireqdLZ2a35+ASMjAAB8TQSRZthVXCafsbTeN0CSdIHjY3mN0e7iEzZXBgBAeCOINEN2erIclvSu7xxVGafSLbeclqVe6Ul2lwYAQFhjZtVmyHQlKi93kO7LL9ffPZeowkrSfbkDlelKtLs0AADCGkGkmaaO6KEx/TK0u/iEeqUnEUIAADgDCCItkOlK/CqAVFVIsQn2FgQAQJjjGJGWOlAg/fGb0hPj7a4EAICwx4hIS6V2kQo3ScYnufdJrm52VwQAQNhiRKSlktOl7qP8/972sr21AAAQ5ggirTFgov9+2wv21gEAQJgjiLRGTRDZ/Y5UesDeWgAACGMEkdbo0FvqNlIyXumjv9tdDQAAYYsg0lpDvitJqtr4V7274zDXnQEAoBUIIq117tXa1/lS3X7gck17/D1duGgVV+QFAKCFOH23lYoq4zVmzwz5jP+xMdL8/AKN6ZfBrKsAADQTIyKt5L8ib91lXJEXAICWYUSklWquyNveuHV9zGuSLP3Wew1X5AUAoAUYEWmlmivynuPYq9tjlutW50r9ZkI7dssAANACBJGvYeqIHlp81w90LPMixVleffvwn+wuCQCAsEIQ+ZoyXYlqd9V9/gcFz0s737S3IAAAwghB5EzIHCyNmClJql45W+9t2828IgAANANB5EwZ93MdT+yqmJJ9OvL0TF246A3mFQEAoAkEkTOkqMKpG9wzVWmcGu34WFk6rPn5BYyMAADQCE7fPUN2FZfpA18/zamapY9Mb+0znST55xXhTBoAAOpHEDlDauYVedl3QWCZ07J0VtU2FR09T7u+9Cg7PZlQAgBALeyaOUNq5hVxWpYkfwj5wyVVav/sVdr90GWa+/g/uR4NAACnYETkDJo6oofG9MvQ7uIT6pWepPgdr6jS51CO42O9GneXHvdeocX5xzWmX4Yk/+4cRkkAANEs6CMiS5YsUa9evZSQkKBRo0bp/fffD/YmbZXpSlROn47KdCVqW7sxmlC5SO94z1WiVanbYlborbjbdPCpmZpz/xJ994/rAqMkRe5yvbuzuM7BrfUtAwAgkgQ1iDz77LP60Y9+pAULFmjjxo0aPHiwxo8fr0OHDgVzs21Gdnqy9qqLplXN182Vc7TV11PJlkdDil/Q3TFPSZJ8Rpr7jy36wf2/16w/vh4IJs+u36MLF61qcVg5U22C+dzUGH5t7N4+NVJjW9p+a2tE/SxjjGm6WeuMGjVKI0aM0O9+9ztJks/nU/fu3XXbbbdp7ty5Ta5fUlIil8slt9uttLS0YJUZVM+u36P5+QXyGiOnJf188DHFFDynQtNRj3hzJUlJqlBB/Aw5LKNSk6gi00FF6qhDvnZyK1nrff31mm+kLEtymmpd6tyk7+b0lc8Ro9+/vU+VJkbVVowuHpCp/E/KVGg6ymFJeVefq9Sy3frf1z+Tz1iyLGncOV306tZD8spShRJ0R+6FkqR5+R+ps/lSDku66/IBkqT7X9kmn5EclnTpoB7665aywOOHL884rc2EgV30r4IDqjAxOmq5lJc7SJL0cP4ayZg6bWrWuWPiQHni0zUvf4t8RupiHdXEgZ3qtKmp51evfKqDpr2/b7mDlFBxWItf3nra9n1GMpZDP8wdc7JvW9TRHFOc5T1t+3ddPsDf9uWDgWXfG5Sk17fsDTyee3L7i17ZJq+xdNjqoKuHdtXyTfvlMiVKtKrqtKm9Xnlil0DfOlilmjywQ53t116vyLSXZTl09dCuWrVpu+KN57Q2NeuNHnKu8jf7n6eddVz3Xp59Wpuavh407WQsp/JyBym2qkQPvrCx3t+Zw5LmTL5QxhGreflblGROKMXyaN4V/u3nvfzVc18+sIv+VlAmj4mRw5IeuKqPYqrLdN/Ln9Rqk6lXCorkM1KplaKFucMkSQvzP1CSqfC3GdRFr2z5avvzJw5QVUya7lq5XT4jJVqVunZQO728pSjQ5qcTz/a/H176RCUmUVVWnPJyB8nhrdSilesD7SYOytRLtdabO+l8+Zzxmpe/RU5TrVSr/LQ2d3/7HP975sWdqjBxcljSlCGd9e/NOwNt7jnZ5t4XP5bPSB4rThOHZmv5pv2yjFcpVoUWXOlvs/CFj2WMZFnSgivPkdcRr7tWfiqfkWIsn6ae114vfvTV9muvV2FiVGXF6eqhXbVi014lnHw//PyqcyVJP//nV+/9KwZ313MfFstnJKfl0+KrzpIkLajV5srBmXrhwyJVGUfgdyZjtHD5hkCbeyf5n/uelf71ZDk0cWgvLd+0Xz5jlGhV6Rcn29y98qvnvmpwllZ8eCDwfvC/Hh7dvbIg0P9Jg7O08sNCmZPPe2/u0MDn02mq5bCkSUOytHJzYeB5fzl5oIykuSu+eu8tmux/P/5sRUFg2aQhXbVy8355jRV4r0vST/M/DLSZPCRLK2o9931XD5SxnIHPp9Py6eohXbV88/5abfzPM3/5FnmNQw5Lgc9+TZuabdU8T2PLxvTLiJpd8i35/g5aEKmsrFRSUpKef/55TZ48ObB8+vTpOnbsmFauXHnaOh6PRx6PJ/C4pKRE3bt3D+sgIvlTcc1xI5J04aJV/g/5Sb2sIj0Zu1i9HAfrXf+v1ZfqZ9UzJEkdVKKNCbc0uK3nqsfojmr/z1NVri0JMxps+6J3lH5QdbtkST5jtDthWoNt3/AO1YyqOwKPt8VPV4JVVW/btd5zdF3Vz/zDbZb0QdxN6mAdr7ftZl8f5Vb9IvD7eDv+B+pmFdfbdruvm8ZXLpbkH8p7Le4n6usorLftXl+GLq78zcm+Sf+M+6nOc+yqt22xSdP5nkcDj5+Nu1ejHNvqbVtm4nWuZ2ng8ZOx9+sS54f1tpWk3p6/Bfr2f7EP6wpnw7sm+1c8KY/iJEn/G/t7TXH+p8G2Qyse1VH5PxO/jPmT/jvmjQbbXljxG+1XhhyS5sU8rZkxLzXY9lueB7RTXeUz0g9jntftMfkNtr3S80ttMb0lSbc6/6m7Ypc12Haq526tN2dLlvTfjld1b+yfG2x7Y9WdetM7RJJ0jXO1Hoj9Q4Ntb628Xa/4RskhaaJzrR6JfaTBtj+pvEX5Zox8RvqmY6OeiHuwwbY/q7pRf/V+S5KU49iqZ+J+1WDbX1V9V3/0fluSNMTaoRXx9zTY9uHqKXq4eookqZ+1V6/F39Vg20err9Si6uskSd2sQ3o7fk6Dbf9c/S0tqL5RkpQutz5IuLXBtjV/IxySkqwKFcT/T4NtX/ReoNlVPzj5qPl/IxySPm7ib8S0qp8FPp8b4xv/GzG58heBx+/E/0Bdm/gbUfO357XYJv5GVP0m8Pls7G/EYZOmEaf8jTjf2q4KxcmjWHkUJ4+JlUexOq5EXVO5MNC3Kx3vqrN1VAdMR+1RJ+32ddZxK1l5uYM0dUQPFbnLIzKctCSIBO1g1eLiYnm9XnXu3LnO8s6dO2vbtvr/yOfl5WnhwoXBKsk2ma7EOm+wvNxBgVESh6QvTKYuqXxIiapQpvWlulpH1MX6Uhk6plTrhDb5+gbWNZLW+/opTtWKU7ViT97irGrFyKtSJQXa+iQdM8myZGSdXNuSZMnIISOPYuWreVJJHhPbYB+q5azzuFKxshqIsFUn29Y8d6ViG3zuKjnrhLIq45Sngbdl7Rp8J9f1mPrbViqmTt+qm2hbpybTcNuqU9s28ryS6vTNK0ejbWurbqSG09o2UUOglpM1VBpno21qavYZS1WNtK27niWvsZrcvmo9d0Na879GvmasZ1T39bBD8MafW672ax2M525Wm1ZsvzmrNPe5v07/nZZRsjxK1sn/eT75li4xiXW2f61ztb7hLKiz7hGTqoIXeuuDraP039svVIWJDYyaTB3Ro/VFhamgjYgUFhaqa9euevfdd5WTkxNYfuedd2rNmjVat27daetE6ohIfWqPkrz16eFau28s3Zc7UJLqhBWjup+rmsTf2AeptW2C+dzU2LZqbE4bu7ffeBsjyTq5zDRjPSuwnlXPN1Xt5zYnD6Gz5JOjVttTazQnW9S0ddY8quf36LMcqjZftY2Vt8Ht+2Sp+mTwteRTnKpP/tv/n9p/ub1yBNpKRkmqbPD36JNDlYoN/M7izVd/c0+t2ydHYJTO/7yeBl8jb622Dkkp1okGXw9/24TA86ToRM33+Gm/D59lqdR89T9yLp2QZfnqfW4jh0qVFKgxxZTV+zqf/LFKrJTA86SpTM4GXg+vsXRMqYHlaSpTvCoVb1UpUVVKsCoVa6oUb1XJKZ/e9Q0K9G2m80UNdOxWV6tYPa2DyrDcgec5YlI13POoalJMlnVUv7/1CpVV+cJ+hKQlIyJBO1g1PT1dTqdTBw/W3d1w8OBBdenSpd514uPjlZaWVucWqWqfXTN1RA+9PXesnpl5gd6eO1ZTR/Sos+yded/Uoil15yjJmzLotHlLpgzrekbaBPO5qbFt1Rj+/bBqLTtPTstxcplDU4Z1q/M4b8p5ddZzWE7lDusuh+WUkUMOy6n7pgzWfbmD5bCcge3lDusuWTHyyilZMfrllCH6Ze6QwDLLcgZqNHLIWLH6xZSh+kXuUBkrVtWKCSz7Ve7gQF8cllNXDuslrxWnSsXKa8Vp4ZRhWpg7TF4rTtWKCfTfYTnlUZyqrXgtnDJcC3OHq9qKDyybNKxnrd+RQwumnK8FueerykpQheJVZSVo4rDeqrISVKnYOr+zSitR5UpQpZWoe6aM0D25IwLLqq34Wr9/Sx4rUXdPGaG7c0fIYyXqhBLksRJ1+bA+qrbi67xGP80dpXIrWceVpHIrWeOHnRV4XGkl1nkdjytJJ6xkXTbsLJ2wklV68vH8KaP0s9yRdV7/+VNGaV5ujsqsVJUoRWVWqr41rL/KrFSVKqnOe6TMSpFbKTpupWrcsP46bqUGHs+bklPnvVZmpeibw85WieXSUaWpxHLpzikX6s7cC1VqpdV5P5ZZKTqs9ipSZ30/d4L+++qr9KH6a63vXK0159Xp2x+939YPq2brmsqfa4Tn9zqn4gld5fmFFlRN1++rr1JNCHHIp2WxP1fS4xfp2T/9ry5a9O+omXcq6Aerjhw5Uo884t9v6/P51KNHD82ePTtqDlY9k2qPotQk5VOXnak2wXxuagy/NnZvnxqpsS1tvzU1njryfeeE/rr/X9sCIzJnWfu0Iu5uJVv+EaqPfNmaV32zHr/rhrAcGWkTB6tK/tN3p0+frscee0wjR47Uww8/rL///e/atm3baceO1IcgAgCIFKeGldpnVTokJemErne+pltjXlCqVS6PiVHR+Xcp/huztevIibDaXdNmgogk/e53v9MDDzygAwcOaMiQIfrtb3+rUaNGNWtdgggAIJLVhJOkOIeu/r935TP+M5/ui31clzk3SJKe816seVUz5LNiwuaA1jYVRL4OgggAIFqcOu/Unwd+qAs+fUDrfQM0vequwPE9b88d2+ZHRtrE6bsAAKD5Tr1e2a7iUbp+q1MFvmxVyj8FgtcY7S4+0eaDSEsQRAAAaCNOnXfqPTOwzrwsvaxDgckxI0XQL3oHAABaLtOVGDgN2CGffhbztN5I+Ikyj39sd2lnFCMiAAC0UYHdNYfLNGTds3J+Vq2q/Fv1wfjl6tW5Q0TsomFEBACANizTlaicvulKnPwbVcR1UOyR7Vr753sCV2YPdwQRAADCQFF1ku4s+64k6ZaYF9TRHNX8/AIVucttruzrIYgAABAGdhWX6Z/eHG309VWS5dEPY/4ROIsmnBFEAAAIA9npyXJYln5VNU2S/8q+PazisD+LhiACAEAYqDmLZrMG6D/egSpVkn55kTPsD1jlrBkAAMJEzVk0hbu7ydM1S2PS0+0u6WsjiAAAEEYyXYnKHHye3WWcMeyaAQAgXBkj7X3ffx+mCCIAAIQjn0969CLpT9+S9m+wu5pWI4gAABCOHA6p80D/vzc9ZW8tXwNBBACAcDX4O/77T16QvNX21tJKBBEAAMJVr4ukxPbSiSMqeO+VsJxllSACAEC4csbq844XS5I+eOUvYXn9GYIIAABhqshdrvt29ZUkTXCulzG+sLv+DEEEAIAwtau4TG95B6nMxKuLdVQDrL1hd/0ZJjQDACBMZacnq9qK1YLqG7TfpGunyZLTssLq+jMEEQAAwlTN9Wfm51vyGiOnZem+3IFhdf0ZgggAAGGs5vozu4tPqFd6UliFEIljRAAACHuZrkTlOLYq892f+6d8DyOMiAAAEAk2/0368BkpPk3qPtLuapqNEREAACJB91H++73v2VtHCxFEAACIBDVBZN+GsJrunSACAEAkyBjg3y1TVSYd2mp3Nc1GEAEAIBI4HFK3Ef5/h9EBqwQRAAAiRY8L/Pd7wuc4EYIIAACRoutw//3RXfbW0QKcvgsAQKToOVr64cdSWpbdlTQbQQQAgEgRmyi5utpdRYuwawYAANiGIAIAQCTZ9R/pmeuk1++xu5JmIYgAABBJKtzS9pelnavsrqRZCCIAAESSzPMkSb5D21T0pdvmYppGEAEAIII8+6lRiUmSw1el/3nwb3p2/R67S2oUQQQAgAhR5C7XvOUF+sz4z5zprULNzy9Qkbvc5soaRhABACBC7Couk89IO3z+IHKWY5+8xmh38QmbK2sYQQQAgAiRnZ4shyXtMP4JzfpahXJalnqlJ9lcWcMIIgAARIhMV6Lycgfpc3VTmYmXkaX7cgcq05Vod2kNsowxxu4iGlJSUiKXyyW32620tDS7ywEAICwUHT2u3UfK1SsjxZYQ0pLvb6Z4BwAgwmS2T1Fm+xS7y2gWds0AAADbEEQAAIhE7z0qLblAWvt/dlfSKIIIAACRqPxL6fAn0uFtdlfSKIIIAACRqH22//7oblvLaApBBACASNShJojssreOJhBEAACIRO17+e/d+6TqSltLaQxBBACASJTSWYpNkoxPcu+1u5oGEUQAAIhElvXVqEgb3j3DhGYAAESqzudKsqS2O4k6QQQAgIg15XG7K2gSu2YAAIBtCCIAAMA2BBEAACKVe5/0fznSQ4PsrqRBHCMCAECkSnBJhz72/9tTKsWn2ltPPRgRAQAgUsWn+sOIJLn321tLAwgiAABEsrRu/nv3PnvraABBBACASObq6r8vIYgAAIBQc9WMiLBrBgAAhFpazYgIQQQAAIRax75Sp3Ol1C52V1IvTt8FACCSnTvZf2ujgjIisnv3bs2YMUPZ2dlKTExUnz59tGDBAlVWVgZjcwAAIEwFZURk27Zt8vl8euyxx9S3b18VFBRo5syZKisr04MPPhiMTQIAgKYYI1mW3VXUYRkTmmsDP/DAA/r973+vzz//vNnrlJSUyOVyye12Ky0tLYjVAQAQwZ78tlT0oXTDi1Lm4KBvriXf3yE7RsTtdqtDhw6NtvF4PPJ4PIHHJSUlwS4LAIDIV3FM8pRIpQelTLuLqSskZ83s2LFDjzzyiG6++eZG2+Xl5cnlcgVu3bt3D0V5AABEtpSTZ8wcP2BvHfVoURCZO3euLMtq9LZt27Y66+zfv18TJkzQNddco5kzZzb6/PPmzZPb7Q7c9u7d2/IeAQCAulI7++9LD9pbRz1atGvmxz/+sW644YZG2/Tu3Tvw78LCQo0dO1ajR4/WH/7whyafPz4+XvHx8S0pCQAANCUwIhLmQSQjI0MZGRnNart//36NHTtWw4cP19KlS+VwMHcaAAC2SG27u2aCcrDq/v37dckll6hnz5568MEHdfjw4cDPunRpmzO7AQAQsVL8u2ZKivepzF2uTFeizQV9JShB5PXXX9eOHTu0Y8cOdevWrc7PQnS2MAAAOOm1wlh18/VUQVGa5i5apbzcQZo6oofdZUkK4TwircE8IgAAfD1F7nJduGiVfLW+7Z2Wpbfnjg3ayEhLvr85cAMAgAi2q7isTgiRJK8x2l18wp6CTkEQAQAggmWnJ8tRZ1Z3I6dlqVd6kl0l1UEQAQAggmW6EpWXO0h/in1QH8V/Xxc5PtZ9uQPbzAGrIZviHQAA2GPqiB7yfJSo+L0n9MhVWWrfRg5UlRgRAQAgKsS3y5Iktfd+aXMldRFEAACIBifnEmlrs6sSRAAAiAbJ6f77E4yIAACAUEuqCSLF9tZxCoIIAADRoGZEpIwgAgAAQi0tS+pynpQxwO5K6uD0XQAAokHmYOmW/9hdxWkYEQEAALYhiAAAANsQRAAAiBZPXC7l9ZAKN9ldSQBBBACAaFFZKnncUtkRuysJIIgAABAt2uBcIgQRAACiRRucS4QgAgBAtEjq6L9nRAQAAIRcEiMiAADALsk1IyIcrAoAAEKtXQ//DKvts+2uJIAp3gEAiBZ9x/lvbQgjIgAAwDYEEQAAYBuCCAAA0cJbLf1msLSoh1R+1O5qJBFEAACIHs4Y6fhhqcJNEAEAADZIbOe/J4gAAICQS2zvvy8/ZmsZNQgiAABEk0AQYUQEAACEWoLLf08QAQAAIceuGQAAYJuOfaXMIVJSB7srkcQU7wAARJeL5vhvbQQjIgAAwDYEEQAAYBuCCAAA0WTfBv80749/y+5KJHGMCAAA0cUZIx3dLVVV2F2JJEZEAACILgnt/PflRyVjbC1FIogAABBdauYR8XqkqnJ7axFBBACA6BKfKllO/7/bwOyqBBEAAKKJZX01KlJxzNZSJIIIAADRJ7Gd/74NjIhw1gwAANGm0zlSXIrkiLW7EoIIAABRZ+pTdlcQwK4ZAABgG4IIAACwDUEEAIBo8+7vpIfPk1bfb3clBBEAAKJOZZl07AuptNDuSggiAABEnYQ0/31Fib11iCACAED0iT8ZRDwEEQAAEGqMiAAAANswIgIAAGzThkZEmFkVAIBok9heap8tpWbaXQlBBACAqNO+l3T7ZrurkMSuGQAAYCOCCAAAsA1BBACAaPTU1dJvBkuHPrG1DIIIAADR6Nge6ehu6cSXtpZBEAEAIBq1kblECCIAAESjNjKXCEEEAIBodHJE5PP9hSpyl9tWBkEEAIAo9HmpU5L0/Dsf68JFq/Ts+j221BH0IOLxeDRkyBBZlqXNmzcHe3MAAKAJRe5yvbm7QpKUZpXLZ6T5+QW2jIwEfWbVO++8U1lZWfrwww+DvSkAANAMu4rLtMo3RF9WpWqjOUuS5DVGu4tPKNOVGNJaghpEXnnlFb322mv6xz/+oVdeeSWYmwIAAM2UnZ6stWaQ3vENCixzWpZ6pSeFvJag7Zo5ePCgZs6cqaeeekpJSc3rmMfjUUlJSZ0bAAA4szJdicrLHSSnZUnyh5D7cgeGfDRECtKIiDFGN9xwg2655Radf/752r17d7PWy8vL08KFC4NREgAAqGXqiB4a0y9Du4tPqFd6ki0hRGrhiMjcuXNlWVajt23btumRRx5RaWmp5s2b16Ji5s2bJ7fbHbjt3bu3ResDAIDmy3QlKqdPR9tCiCRZxhjT3MaHDx/WkSNHGm3Tu3dvXXvttXrhhRdknRzykSSv1yun06lp06bpz3/+c7O2V1JSIpfLJbfbrbS0tOaWCQAAbNSS7+8WBZHm2rNnT53jOwoLCzV+/Hg9//zzGjVqlLp169as5yGIAAAQflry/R2UY0R69OhR53FKSookqU+fPs0OIQAAIPIxsyoAALBN0Cc0k6RevXopCHuAAABAmGNEBAAA2IYgAgAAbEMQAQAAtiGIAAAA2xBEAACAbQgiAADANiE5fbe1ak755Sq8AACEj5rv7eZM3dGmg0hpaakkqXv37jZXAgAAWqq0tFQul6vRNkG51syZ4vP5VFhYqNTU1DoX0DsTSkpK1L17d+3duzcir2ND/8JfpPeR/oW/SO9jpPdPCl4fjTEqLS1VVlaWHI7GjwJp0yMiDocj6NemSUtLi9g3mET/IkGk95H+hb9I72Ok908KTh+bGgmpwcGqAADANgQRAABgm6gNIvHx8VqwYIHi4+PtLiUo6F/4i/Q+0r/wF+l9jPT+SW2jj236YFUAABDZonZEBAAA2I8gAgAAbEMQAQAAtiGIAAAA20RFENm9e7dmzJih7OxsJSYmqk+fPlqwYIEqKysbXa+iokKzZs1Sx44dlZKSoilTpujgwYMhqrrlfvWrX2n06NFKSkpSu3btmrXODTfcIMuy6twmTJgQ3EJbqTX9M8bonnvuUWZmphITEzVu3Dh99tlnwS20lb788ktNmzZNaWlpateunWbMmKHjx483us4ll1xy2ut3yy23hKjipi1ZskS9evVSQkKCRo0apffff7/R9s8995wGDBighIQEDRo0SC+//HKIKm2dlvTvySefPO21SkhICGG1LfPWW2/pyiuvVFZWlizL0ooVK5pcZ/Xq1Ro2bJji4+PVt29fPfnkk0Gv8+toaR9Xr1592mtoWZYOHDgQmoJbKC8vTyNGjFBqaqo6deqkyZMna/v27U2uF+rPYVQEkW3btsnn8+mxxx7T1q1b9dBDD+nRRx/V/PnzG13vhz/8oV544QU999xzWrNmjQoLC5WbmxuiqluusrJS11xzjW699dYWrTdhwgQVFRUFbs8880yQKvx6WtO/xYsX67e//a0effRRrVu3TsnJyRo/frwqKiqCWGnrTJs2TVu3btXrr7+uF198UW+99ZZuuummJtebOXNmnddv8eLFIai2ac8++6x+9KMfacGCBdq4caMGDx6s8ePH69ChQ/W2f/fdd3XddddpxowZ2rRpkyZPnqzJkyeroKAgxJU3T0v7J/lnr6z9Wn3xxRchrLhlysrKNHjwYC1ZsqRZ7Xft2qWJEydq7Nix2rx5s+bMmaPvf//7evXVV4Ncaeu1tI81tm/fXud17NSpU5Aq/HrWrFmjWbNm6b333tPrr7+uqqoqXXbZZSorK2twHVs+hyZKLV682GRnZzf482PHjpnY2Fjz3HPPBZZ98sknRpJZu3ZtKEpstaVLlxqXy9WsttOnTzeTJk0Kaj1nWnP75/P5TJcuXcwDDzwQWHbs2DETHx9vnnnmmSBW2HIff/yxkWTWr18fWPbKK68Yy7LM/v37G1zv4osvNrfffnsIKmy5kSNHmlmzZgUee71ek5WVZfLy8uptf+2115qJEyfWWTZq1Chz8803B7XO1mpp/1ryuWxrJJnly5c32ubOO+805557bp1lU6dONePHjw9iZWdOc/r45ptvGknm6NGjIanpTDt06JCRZNasWdNgGzs+h1ExIlIft9utDh06NPjzDRs2qKqqSuPGjQssGzBggHr06KG1a9eGosSQWb16tTp16qT+/fvr1ltv1ZEjR+wu6YzYtWuXDhw4UOc1dLlcGjVqVJt7DdeuXat27drp/PPPDywbN26cHA6H1q1b1+i6Tz/9tNLT0zVw4EDNmzdPJ06cCHa5TaqsrNSGDRvq/O4dDofGjRvX4O9+7dq1ddpL0vjx49vcayW1rn+SdPz4cfXs2VPdu3fXpEmTtHXr1lCUGxLh9Pp9XUOGDFFmZqa+9a1v6Z133rG7nGZzu92S1Oh3nx2vY5u+6F2w7NixQ4888ogefPDBBtscOHBAcXFxpx2L0Llz5za7P7A1JkyYoNzcXGVnZ2vnzp2aP3++Lr/8cq1du1ZOp9Pu8r6Wmtepc+fOdZa3xdfwwIEDpw3vxsTEqEOHDo3W+t3vflc9e/ZUVlaWPvroI911113avn278vPzg11yo4qLi+X1euv93W/btq3edQ4cOBAWr5XUuv71799fTzzxhM477zy53W49+OCDGj16tLZu3Rr0i3uGQkOvX0lJicrLy5WYmGhTZWdOZmamHn30UZ1//vnyeDx6/PHHdckll2jdunUaNmyY3eU1yufzac6cObrwwgs1cODABtvZ8TkM6xGRuXPn1nvgUO3bqX8U9u/frwkTJuiaa67RzJkzbaq8+VrTx5b4zne+o6uuukqDBg3S5MmT9eKLL2r9+vVavXr1metEI4LdP7sFu3833XSTxo8fr0GDBmnatGn6y1/+ouXLl2vnzp1nsBc4E3JycnT99ddryJAhuvjii5Wfn6+MjAw99thjdpeGZurfv79uvvlmDR8+XKNHj9YTTzyh0aNH66GHHrK7tCbNmjVLBQUFWrZsmd2lnCasR0R+/OMf64Ybbmi0Te/evQP/Liws1NixYzV69Gj94Q9/aHS9Ll26qLKyUseOHaszKnLw4EF16dLl65TdIi3t49fVu3dvpaena8eOHbr00kvP2PM2JJj9q3mdDh48qMzMzMDygwcPasiQIa16zpZqbv+6dOly2kGO1dXV+vLLL1v0fhs1apQk/6hfnz59WlzvmZKeni6n03naWWaNfX66dOnSovZ2ak3/ThUbG6uhQ4dqx44dwSgx5Bp6/dLS0iJiNKQhI0eO1Ntvv213GY2aPXt24AD4pkbf7PgchnUQycjIUEZGRrPa7t+/X2PHjtXw4cO1dOlSORyNDwYNHz5csbGxeuONNzRlyhRJ/iOl9+zZo5ycnK9de3O1pI9nwr59+3TkyJE6X9zBFMz+ZWdnq0uXLnrjjTcCwaOkpETr1q1r8ZlFrdXc/uXk5OjYsWPasGGDhg8fLklatWqVfD5fIFw0x+bNmyUpZK9fQ+Li4jR8+HC98cYbmjx5siT/0PAbb7yh2bNn17tOTk6O3njjDc2ZMyew7PXXXw/p5625WtO/U3m9Xm3ZskVXXHFFECsNnZycnNNO82yrr9+ZtHnzZts/bw0xxui2227T8uXLtXr1amVnZze5ji2fw6AdBtuG7Nu3z/Tt29dceumlZt++faaoqChwq92mf//+Zt26dYFlt9xyi+nRo4dZtWqV+eCDD0xOTo7JycmxowvN8sUXX5hNmzaZhQsXmpSUFLNp0yazadMmU1paGmjTv39/k5+fb4wxprS01PzkJz8xa9euNbt27TL//ve/zbBhw8xZZ51lKioq7OpGg1raP2OMWbRokWnXrp1ZuXKl+eijj8ykSZNMdna2KS8vt6MLjZowYYIZOnSoWbdunXn77bfNWWedZa677rrAz099j+7YscPce++95oMPPjC7du0yK1euNL179zZjxoyxqwt1LFu2zMTHx5snn3zSfPzxx+amm24y7dq1MwcOHDDGGPO9733PzJ07N9D+nXfeMTExMebBBx80n3zyiVmwYIGJjY01W7ZssasLjWpp/xYuXGheffVVs3PnTrNhwwbzne98xyQkJJitW7fa1YVGlZaWBj5jksyvf/1rs2nTJvPFF18YY4yZO3eu+d73vhdo//nnn5ukpCRzxx13mE8++cQsWbLEOJ1O869//cuuLjSppX186KGHzIoVK8xnn31mtmzZYm6//XbjcDjMv//9b7u60Khbb73VuFwus3r16jrfeydOnAi0aQufw6gIIkuXLjWS6r3V2LVrl5Fk3nzzzcCy8vJy8//+3/8z7du3N0lJSebqq6+uE17amunTp9fbx9p9kmSWLl1qjDHmxIkT5rLLLjMZGRkmNjbW9OzZ08ycOTPwh7StaWn/jPGfwnv33Xebzp07m/j4eHPppZea7du3h774Zjhy5Ii57rrrTEpKiklLSzM33nhjnZB16nt0z549ZsyYMaZDhw4mPj7e9O3b19xxxx3G7Xbb1IPTPfLII6ZHjx4mLi7OjBw50rz33nuBn1188cVm+vTpddr//e9/N/369TNxcXHm3HPPNS+99FKIK26ZlvRvzpw5gbadO3c2V1xxhdm4caMNVTdPzamqp95q+jR9+nRz8cUXn7bOkCFDTFxcnOndu3edz2Jb1NI+3n///aZPnz4mISHBdOjQwVxyySVm1apV9hTfDA1979V+XdrC59A6WSwAAEDIhfVZMwAAILwRRAAAgG0IIgAAwDYEEQAAYBuCCAAAsA1BBAAA2IYgAgAAbEMQAQAAtiGIAAAA2xBEAACAbQgiAADANgQRAABgm/8PkrlLH9E5QbQAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "aSVD = VT.T @ np.diag(1/Sdiag) @ U.T @ y\n",
    "\n",
    "print(\"Fitted coeffs = \", aSVD)\n",
    "print(\"frac error = \", (aSVD-a0)/a0)\n",
    "\n",
    "plt.plot(x,y,'.')\n",
    "\n",
    "xx = np.linspace(-2,2,1000)\n",
    "poly = np.polynomial.Polynomial(aSVD)\n",
    "yy = poly(xx)\n",
    "plt.plot(xx,yy,'--')\n",
    "\n",
    "dev = np.sqrt(np.mean((y-poly(x))**2))\n",
    "print('rms deviation / max(y) = %g' % (dev/np.max(np.abs(y)),))\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4fab5983",
   "metadata": {},
   "source": [
    "## Orthogonal polynomials"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7a3ee4b8",
   "metadata": {},
   "source": [
    "SVD does a lot better than straight inversion of the normal equations, but still fails for large enough polynomial order. One way that we can do better is to choose a set of basis functions that are orthogonal in the domain that we are interested in. For example, the Legendre polynomials are orthogonal on the domain $x=(-1,1)$ (which you can rescale to by rescaling your $x$ variable). "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b10f5bf9",
   "metadata": {},
   "source": [
    "> **Exercise**:  Try changing the basis functions to Legendre polynomials and see how that improves the fit, both for the inversion and SVD methods. How do the results change when you change the number of data points and the order of the polynomial? Check also how the condition number of the design matrix changes as a function of the polynomial degree when you change the polynomials you are using to fit. Hint: You can use [`numpy.polynomial.legendre.legvander`](https://numpy.org/devdocs/reference/generated/numpy.polynomial.legendre.legvander.html) to generate the design matrix $\\mathbf{A}$ for you."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}