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.08183424 -0.43254686 -0.18893     0.60685865  0.96228237 -0.49706461\n",
      " -0.91780677 -0.17740619  0.29815188 -0.84552489 -0.08772226]\n",
      "Fitted coeffs =  [-0.19930902 -0.43254686  6.56500165  0.60685868 -1.85188364 -0.49706461\n",
      " -0.48267869 -0.17740619  0.2523617  -0.84552489 -0.09221752]\n",
      "frac error =  [ 1.43552093e+00  4.29743339e-09 -3.57483279e+01  4.46997345e-08\n",
      " -2.92447009e+00  4.75474501e-09 -4.74095528e-01  1.42266096e-08\n",
      " -1.53580061e-01 -4.64360780e-10  5.12442337e-02]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjMAAAGdCAYAAADnrPLBAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA7RElEQVR4nO3deXwU9f3H8ffmTiAJhFwgCYRDUIkg4QpWjoIERSvFKlKqaCkoDf2JeAFarPZXoWillZ8VrRWtVqG0SCvVauSySkANoIAE5ZJADgzCJkDItfP7Y5INgSQkIbuzu3k9H495zOzMdzafL0vYN9+5bIZhGAIAAPBSflYXAAAAcDEIMwAAwKsRZgAAgFcjzAAAAK9GmAEAAF6NMAMAALwaYQYAAHg1wgwAAPBqAVYX4A4Oh0O5ubkKDw+XzWazuhwAANAIhmGouLhYnTp1kp9f/eMvrSLM5ObmKiEhweoyAABAM+Tk5Khz5871bm8VYSY8PFyS+YcRERFhcTUAAKAxioqKlJCQ4Pwer0+rCDPVh5YiIiIIMwAAeJkLnSLCCcAAAMCrEWYAAIBXI8wAAACvRpgBAABejTADAAC8GmEGAAB4NcIMAADwaoQZAADg1QgzAADAqxFmAACAVyPMXIQ8e4k27StUnr3E6lIAAGi1WsWzmVxhxaeHNHfVDjkMyc8mLZiQrIkDE60uCwCAVoeRmWbIs5c4g4wkOQxp3qqdjNAAAGABt4WZhQsXymazadasWc51Z86cUXp6ujp06KC2bdvq5ptvVkFBQa39Dh06pHHjxiksLEyxsbF68MEHVVFR4a6y63Sg8JQzyFSrNAwdLDxtTUEAALRibgkzn376qV544QVdeeWVtdbfd999evvtt7Vy5Upt3LhRubm5mjBhgnN7ZWWlxo0bp7KyMm3atEmvvvqqXnnlFc2fP98dZdcrKbqN/M55Grm/zaau0WHWFAQAQCvm8jBz8uRJTZ48WX/605/Uvn1753q73a4///nPeuaZZ/T9739fKSkpWrZsmTZt2qTNmzdLkt5//319+eWXev3119WvXz9dd911+vWvf63nnntOZWVlri69Xh0jQ7VgQrIibGc00JatKNspPTmhjzpGhlpWEwAArZXLw0x6errGjRun0aNH11qflZWl8vLyWut79+6txMREZWZmSpIyMzOVnJysuLg4Z5u0tDQVFRVp165d9f7M0tJSFRUV1Zpa2sSBifos8VmtDH5CayeIk38BALCIS8PM8uXLtXXrVi1YsOC8bfn5+QoKClK7du1qrY+Li1N+fr6zzdlBpnp79bb6LFiwQJGRkc4pISHhIntSt6BL+kqS2tt3u+T9AQDAhbkszOTk5Ojee+/VX//6V4WEhLjqx9Rp7ty5stvtziknJ8c1P6hj1TlA+V+45v0BAMAFuSzMZGVl6ejRo+rfv78CAgIUEBCgjRs36tlnn1VAQIDi4uJUVlamEydO1NqvoKBA8fHxkqT4+Pjzrm6qfl3dpi7BwcGKiIioNblEx37mPI8wAwCAVVwWZkaNGqUdO3Zo+/btzmnAgAGaPHmyczkwMFBr16517rNnzx4dOnRIqampkqTU1FTt2LFDR48edbbJyMhQRESELr/8cleV3nixl0s2P+nUUam4/sNeAADAdVx2B+Dw8HD16dOn1ro2bdqoQ4cOzvVTp07V7NmzFRUVpYiICP3iF79QamqqhgwZIkkaM2aMLr/8ct1+++1atGiR8vPz9eijjyo9PV3BwcGuKr3xgsKk6Eulb7OlvM+l8PpHiwAAgGtYegfgxYsX64YbbtDNN9+sYcOGKT4+XqtWrXJu9/f315o1a+Tv76/U1FT95Cc/0R133KEnnnjCwqrPEV913gyHmgAAsITNMAzjws28W1FRkSIjI2W321v+/Jm9a6Vj+6SkYVJs75Z9bwAAWrHGfn/zoMmL1WOUOQEAAEvwoEkAAODVCDMt4Wi2tO11qfBrqysBAKDVIcy0hHW/lv6ZLn31ntWVAADQ6hBmWkJH87EG3AkYAAD3I8y0BC7PBgDAMoSZltCpnzkv3COVnrS0FAAAWhvCTEsIj5ciLpEMh5S33epqAABoVQgzLeWS/ub8SJa1dQAA0MoQZlrKJSnmnDADAIBbcQfglnL5eCkuWep0ldWVAADQqhBmWkpUkjkBAAC34jATAADwaozMtKTc7dLuf0kdekr9JlldDQAArQIjMy3pyGfSf38n7VhpdSUAALQahJmWdPYVTYZhbS0AALQShJmWFHuF5B8snTkhfbff6moAAGgVCDMtKSCo5qGTR7ZaWwsAAK0EYaalOQ81fWZtHQAAtBKEmZbWeYA5P7TZ2joAAGglCDMtLXGIOS86IlWUWVsLAACtAPeZaWmRnaWZn0lR3SU/siIAAK5GmHGF6J5WVwAAQKvB0AEAAPBqhBlXKDku/f2n0nODJUel1dUAAODTCDOuEBwhfZ0hfZstFeyyuhoAAHwaYcYV/PxrLtHO2WJtLQAA+DjCjKskVF2izf1mAABwKcKMqyQONueEGQAAXIow4yqXDJBs/lLRYelEjtXVAADgswgzrhLcVopPNpcPZVpbCwAAPoyb5rlS0jDJcEj+gVZXAgCAzyLMuNK1T0g2m9VVAADg0zjM5EoEGQAAXI4w4w5lp6XiAqurAADAJxFmXO3TP0sLE6UPfmV1JQAA+CTCjKtFJUmOcungfyXDsLoaAAB8DmHG1RKGSH6Bkj1HOn7Q6moAAPA5hBlXCwqreU7TgQ+trQUAAB9EmHGHrteY84P/tbYOAAB8EGHGHZKqwswBzpsBAKClEWbcofMgKSBUOpkvHd1tdTUAAPgU7gDsDoEh0tX/I4V1kNrGWl0NAAA+hTDjLiPnWV0BAAA+icNMAADAqxFm3Ml+RMp6RTqy1epKAADwGYQZd/rv09Lb90qfv2l1JQAA+AzCjDv1GG3O935gbR0AAPgQwow7JQ2T/AKk7/ZLx/ZZXQ0AAD6BMONOweFSYqq5vG+dtbUAAOAjCDPu1mOUOf/6fWvrAADARxBm3K1nmjnfv1EqPWltLQAA+ADCjLvFXia17yo5KqS87VZXAwCA1+MOwO5ms0m3vCK16yKFRVldDQAAXo8wY4VOV1ldAQAAPoPDTFYzDKsrAADAqxFmrPJ1hvTStdL7j1pdCQAAXo0wY5WKUunwJ9LufzE6AwDARSDMWKX7SCkgRDpxSCrYZXU1AAB4LcKMVYLaSN2/by7vftvaWgAA8GIuDTMLFizQwIEDFR4ertjYWI0fP1579uyp1ebMmTNKT09Xhw4d1LZtW918880qKCio1ebQoUMaN26cwsLCFBsbqwcffFAVFRWuLN09Lr/JnO9axaEmAACayaVhZuPGjUpPT9fmzZuVkZGh8vJyjRkzRqdOnXK2ue+++/T2229r5cqV2rhxo3JzczVhwgTn9srKSo0bN05lZWXatGmTXn31Vb3yyiuaP3++K0t3j17XS/7BUuFXUsFOq6sBAMAr2QzDfUMC3377rWJjY7Vx40YNGzZMdrtdMTExeuONN/SjH/1IkpSdna3LLrtMmZmZGjJkiN59913dcMMNys3NVVxcnCRp6dKlevjhh/Xtt98qKCjogj+3qKhIkZGRstvtioiIcGkfm2z5ZCl7jfS92dLox6yuBgAAj9HY72+3njNjt9slSVFR5p1vs7KyVF5ertGjRzvb9O7dW4mJicrMzJQkZWZmKjk52RlkJCktLU1FRUXatavuE2dLS0tVVFRUa/JYfW+Teo2reZo2AABoErfdAdjhcGjWrFm6+uqr1adPH0lSfn6+goKC1K5du1pt4+LilJ+f72xzdpCp3l69rS4LFizQ448/3sI9cJHLbjQnAADQLG4bmUlPT9fOnTu1fPlyl/+suXPnym63O6ecnByX/0wAAGANt4SZmTNnas2aNVq/fr06d+7sXB8fH6+ysjKdOHGiVvuCggLFx8c725x7dVP16+o25woODlZEREStyeN9t1/a9H+Sw2F1JQAAeBWXhhnDMDRz5ky99dZbWrdunZKSkmptT0lJUWBgoNauXetct2fPHh06dEipqeY5JKmpqdqxY4eOHj3qbJORkaGIiAhdfvnlrizffSrLpRdGSO8/In3zsdXVAADgVVwaZtLT0/X666/rjTfeUHh4uPLz85Wfn6+SkhJJUmRkpKZOnarZs2dr/fr1ysrK0l133aXU1FQNGTJEkjRmzBhdfvnluv322/X555/rvffe06OPPqr09HQFBwe7snz38Q+U+vzQXN7+V2trAQDAy7j00mybzVbn+mXLlunOO++UZN407/7779ebb76p0tJSpaWl6Y9//GOtQ0jffPONZsyYoQ0bNqhNmzaaMmWKFi5cqICAxp2/7NGXZlfL+UT687VSYJj0wFdScLjVFQEAYKnGfn+79T4zVvGKMGMY0v8NlI59Lf1gidT/DqsrAgDAUh55nxk0wGaTrppsLm/jUBMAAI1FmPEkV94m2fyknM3SsX1WVwMAgFcgzHiSiI5S91FSUFupoO67GwMAgNrcdgdgNNINz0hhHaSgNlZXAgCAVyDMeJp2iVZXAACAV+Ewk6cyDCl/h9VVAADg8QgznqiyXHpxhLT0e1LBl1ZXAwCARyPMeCL/QKldgrn86UvW1gIAgIcjzHiqgdPM+efLpTN2a2sBAMCDEWY8VdIwKbqXVH5K2v6m1dUAAOCxCDOeymaTBk83lzc/J1VWWFsPAAAeijDjyfr+2LznzIlD0perra4GAACPRJjxZEFh0uB7zOXdb1tbCwAAHoqb5nm6gT+T4vpIl461uhIAADwSYcbThUVJva+3ugoAADwWh5m8SelJ6USO1VUAAOBRCDPeYu9a6Q99pX/9wupKAADwKIQZb9Ghh3nzvP3rpYMfWV0NAAAegzDjLdp3kfrfYS6v+435IEoAAECY8SrDHpD8g6VDm6R966yuBgAAj0CY8SYRnaSBU83ldb+WHA5r6wEAwAMQZrzN9+6TAttIuduknX+3uhoAACxHmPE2bWOla2aby0e/tLYWAAA8ADfN80apM6Ueo6ROV1ldCQAAlmNkxhsFhhBkAACoQpjxdsf2SR8/a3UVAABYhsNM3uxUobT0e1L5afOmejzDCQDQCjEy483aREuDppvL/75fOlNkbT0AAFiAMOPtRsyR2idJxbnS2iesrgYAALcjzHi7wFDpxj+Yy5++JH2zydp6AABwM8KML+g2XLrqdkmGtOpuqeSE1RUBAOA2hBlfkfak1L6rZD8kbf6j1dUAAOA2XM3kK0IipAkvSV+/Jw170OpqAABwG8KML0kYaE4AALQiHGbyVRVl0sZFXK4NAPB5hBlf9dZ0af1vpNUzJIfD6moAAHAZwoyvSp0p+QdJ2Wuk/z5tdTUAALgMYcZXdR4gjfudubz+N9L2N62tBwAAFyHM+LL+d0hD/8dc/tdMae9aa+sBAMAFCDO+bvTjUvItkqNC+tsd0pEsqysCAKBFEWZ8nZ+fdNMfpaThks2Pk4EBAD6H+8y0BgFB0m1/lb47IHW80upqAABoUYzMtBbB4bWDzKHN5gQAgJcjzLRGhXulNyZKf7lJ2v221dUAAHBRCDOtUUQnKTFVqjgjrbhd2vKiZBhWVwUAQLMQZlqjoDBp4utS/ymSDOndB1Xw+jTlHTtudWUAADQZYaa18g+QbvyDPu99nyoNm+L2rdTRP3xfb3+4xerKAABoEsJMK5ZXdEY//HygppTP0XGjrfr67deR9/9PefYSq0sDAKDRCDOt2IHCU3IY0keOZN1Y9hv9o/IaLa6YoIOFp60uDQCARiPMtGJJ0W3kZzOXDxsxur98hipsweoaHSY5KqXXb5a2v8HJwQAAj0aYacU6RoZqwYRk+dvMRONvs+nJCX3UMTJU+ny5tPcDafUMadl1Uv5Oi6sFAKBuNsPw/f92FxUVKTIyUna7XREREVaX43Hy7CU6WHhaXaPDzCAjSRVl0ubnpI2LpPLTkmzmM55GzpWiullaLwCgdWjs9zdhBg07kSO9/6j05WrztV+AdNVPpOuflvwDLS0NAODbGvv9zWEmNKxdgnTrq9L0DVKP0ebTt4/tI8gAADwGD5pE43S6SvrJP6RvNkn+QTXri/KkN26V+k6SrrxVahNtXY0AgFaJMIOm6TK09uvP35TyvzCnjPnSpWnS5eOlS8dIIZGWlAgAaF0IM7g4KXeaT+Te9rqUt13KXmNOfoFS0jDphmek9l0tLhIA4MsIM7g4YVHSoGnmVLBL2vkPafcaqXCPdPAjKeysw047V5nn3HQeaAacqkvCAQC4GIQZtJy4K8xp1Hyp8GupYKcU3LZm+3+fkQp2mMuhUdIlKebUsa8Ue5kUlWRN3QAAr0aYgWtE9zSnaoYh9fi+eRVU/g6p5Dtpb4Y5SVKHntIvPqtp/9nLUlBbqV0XqX0XqW0cIzkAgDoRZuAeNpt07RPmckWpeUfhI1nmVLBLiu1d09YwpPcelcpP1awLCJHaJUrhHaXEIdLIeTXbDm0xTzYObS+FtpMCgt3SJQCAZ/CaMPPcc8/pqaeeUn5+vvr27aslS5Zo0KBBVpeF5ggIljqnmFNdKsukK8ZLx7+RTnwjFR2RKs5IhV+ZU2BoTVvDkF4bX3WX4qpVAaGyhbaTQtpJXa+Wxv2upv3GpySbzFGfwDApqE3N1CbGPNxVreyUeRm6XwCjQgDgwbwizKxYsUKzZ8/W0qVLNXjwYP3+979XWlqa9uzZo9jYWKvLQ0sLCJbG/7HmdWW5ZD9sBpvi/Nr3sqkoldolqvREngLLiuRnM2SrKJGKS6TivPOvpPrv02YwqkuXq6W73ql5/ftk6fQxyeYn+QebdQWESAFB5n13bv1LTdtVd0ulRWe1CTaDkH+QFNFJSk2vabv9Dan0pOQfUBWWAs3Db/5BUkiE8qIG6kDhKSVFt1HH0m8kR/lZbQJr9gkIMq8k82B59pKavkSGXngHD+Xz/XBUmuHdUSEZDvO1o0Iyquah7c1JMtt9u8dsY1Se09ZhnvvWobvZtuSE+Yw3R0XtdjZ/8+9x7GVSp35m2/Iz0qFNNX/XA4LN/3QEtamZ+/m7848LXsQrwswzzzyjadOm6a677pIkLV26VP/+97/18ssva86cORZXB5fzDzT/gazrBOHAEOX9ZIOuXrhOhuFQuEoUYTulKNtpLZvUUx3ad6hpaxhSyl1S2UlzJKfsVM1Ufvr84FNRWrWfQ6ooMadqEZfUbrt/vXSyoO7645Jrh5kPn5K+219n0+KwRF19fKEchuRnk7JiHlf7oj11v2/beOmBs7a9eqN0ZJsZkvwCzRElvwDzCyC0nXkX52rvPCjlfV6z3dm2KmDd+mpN2y0vmidzn93GuY+/NHyO+TMl6av3pWN7JZufsnJO6J/b81UpmwzZ9MP+nTXwpplmCJOknE+k7w6YYdFmq5r8zEk2qee1NaNwR3ebgba6jaraVz+NJWGwFBRmLhfuNYOvDMn5sBajpm2X1JoQWLhXOvZ11baq7c4nvBhS1+9pxc5izV21Q911WH399uvWlE4a1KVdzRd59Zf/5T+QIjubux7ZKn39fs22c9umTKkZBTy0Wdr6Wk1wcH7pO8z51bPMmiXpwIfSut+cEyIcNfuNfky67Eaz7b710qppVe0qVV5RrvYV5eoghwLkUFafeUq55SGz7Tcfm39/6nPtE9LV99Z8Fi+Nqr/tiLnSiKp/l4vzpH9Mrb9t6syaMHPqqPTaD+tv2/8O6QdLJEl5hccU/uZNCgiPUUhknBTWwRxZbRNtnl8X3ZNbQrQyHh9mysrKlJWVpblz5zrX+fn5afTo0crMzKxzn9LSUpWWljpfFxUVubxOWOdA4Sk5DEnyU5HaqMhoo8OG9FVYilI7nxVmbDbpuoWNf+MHvpYqS81QUz1Vvz73vJy0J6XSYvMQWcUZ80GdFWfMUZW2cbXb9kyTTuabI06VZVXzcpWVlerDw35VfZEchvTliUANaRsrf6Ncqqww36+yzPwCO/eREqUnpbLiuvtyJqr266O7pZwtdbcNCKn9et9a6av/1P/nNKLmd1OfvyntWiVJSpGUcnaJO6T8YbcrPibGfJ31qrT99frf9/6zDil+9rL0yYv1t/2f7TVhd9tfpI//UH/bGZvMq+4k81YCG56st2nhbe9o7qoTchjSCP/teiTwDekLmdO54q6oCTO5W6UNC+qvodvwmjDz3YGG/xySb6lZPv2dlLO5/rYlJ2qWHRXSqW+dLwMlBZ51tHTN9sPqNKbEHKGxnTPiYfMz150dXKsFhEgRnSU/v6pDsGe38zdDRbWgtlLS8Jr3sFW1cVSaf5ejL639M+P61PxeVJSa582VnjTDW5AZQFd8ekhLVq3TR8E7pGP1/DlceZs04QVzuaJMWvUz8yKD2MukmN5m2OHcOp/i8WGmsLBQlZWViour/YUQFxen7OzsOvdZsGCBHn/8cXeUBw+QFN1GfjY5Q4Ak+dts6hoddnFvHBQmqZHvkfyjxr9vPYHqs32FSv9T7YAxuWye3rx5iFK7d6jduPp/5Web9KY5ylRZftb/8qv+x65zzvn5/i/NL7qz21Qvn+vKiVLnAbXbVO9TWV41UlIlYbBks6mwuESb9x+TnwzZZMivago/Xq746u+62N5St5EyR00c5oiI4aiZzg5rEZdIHfvVbiejpl9ntw3vKMUnV72wnXW+U9Xy2WEtoqN5e4Dq7VJVe3P58Mmav1c5Rqw2VPZVpfx0VZcoRbUNNfvu529+SZ/9JR5zmTkKePZ2P/+a1+3PGmXs2Fca9dj5X/jV+13S/6w/30HSra+d0/asUNGhR+3PYsYmyeavrUeKNetvO+SQnyoMP1XKT8UK05jC02aYSRwiPXr0rDobOEcsvo80e1f928/WLkGa8q/GtY3sLM34+Pz1hmEGG8OhPHuJ5q7aoWAjXD8te0AdbEWKthUrfXCk2lbYpVOF5qHosy8oKDoiffnP2u9p8zcPhXXqL13xQ6nX2MbVCI/l8U/Nzs3N1SWXXKJNmzYpNTXVuf6hhx7Sxo0btWXL+f+7rGtkJiEhgadm+7AVnx7SvFU7VWkY8rfZ9OSEPpo4MNHqspokz16iqxeuOy+UfTRnpNedp+ErfaEfnmXTvkL9+E/n/5v/5rQ6An+1U8ekHX+Tvs2Wjmabo5Kl9prtwx6Uvv+ouVxaLH22zHyobuxlnPjvARr71GyPH5mJjo6Wv7+/Cgpqn49QUFCg+Pj4OvcJDg5WcDBDiK3JxIGJGnZpjA4WnlbX6DCv+ge6WsfIUC2YkHxeKKMv1qEfnqVZo7BtOkhDZtS8NgzzXJ6CXdLhT6WeY2q2HfivlPFLcwrvJPUcLV0xQep6Tc15YfBIHj8yI0mDBw/WoEGDtGSJefKXw+FQYmKiZs6c2agTgBub7ABPkGcv8epQdjZf6Qv98BwuHYU98KF5vtXBj2pf9dgmxnyAburPpahuLfOz0CiN/f72ijCzYsUKTZkyRS+88IIGDRqk3//+9/rb3/6m7Ozs886lqQthBgB8h8tDWXmJeYXX7rfN821Kjpvrf7659r2o4HI+c5hJkiZOnKhvv/1W8+fPV35+vvr166f//Oc/jQoyAADf0jEy1LUjS4Gh5nkzPUZL1z8t7d8g5W6rHWTeedC8JHzw3TX34IFlvGJk5mIxMgMAaDGnCqVnLjMvIQ+OkAZNk4b+glDjAo39/vardwsAADhfaHvphy9IsZebd/7+7++kZ/tLn/zJvB8U3I4wAwBAU/j5S30mSPd8LE183bwRX8l30jsPSEuvNu+uDbcizAAA0Bx+fubjI+752Dy3JjRKKsqT2vDMQHfzihOAAQDwWP4B5nkzyT8yH8IZ0bFm2/GDPCfKDRiZAQCgJYS2Nx8NUW3Pu+a5NBufqnqsCFyFMAMAgCvs32A+J239/0qv31xzvxq0OMIMAACucN1vpZv+KAWESvvXSy+Nlo7ts7oqn0SYAQDAVa6aLP3sAymis3Rsr/TSKOlgHU8Hx0UhzAAA4ErxfaRp66RLUsxDTa9PkAq/troqn8LVTAAAuFp4nHTnv6W//9R8WGWHHlZX5FMIMwAAuENgqHTra+ZN92w25dlLdKDwlJKi23jtU8w9BYeZAABwF/8AyWbTik8P6fsL/6MTr0zS3N8u1opPD1ldmVcjzAAA4EZ59hLNXbVDd/m9o+v9P9HzgYv197dWKc9eYnVpXoswAwCAGx0oPCWHIb1YeYPWV/ZVqK1MLwX+Vvn7dlhdmtcizAAA4EZJ0W3kZ5MqFKAZ5bP0meNSRdpOK/nDu7mxXjMRZgAAcKOOkaFaMCFZ/jabzihYM8rv16nQTgo4sV9aeZdUWWF1iV6Hq5kAAHCziQMTNezSGB0sPK2u0WFqU3KF9Oc0807BGfOlsU9aXaJXYWQGAAALdIwMVWr3DuZl2fHJ0oQXpIhLpF7XWV2a12FkBgAAT3DZjVKPa6XAEKsr8TqMzAAA4CnODjLHD3L+TCMRZgAA8DQ7/i79caj0399ZXYlXIMwAAOBpDIdUfkr6cJGUv9PqajweYQYAAE+TfIvU+wbJUSH98+ccbroAwgwAAJ7GZpPG/U4KiZTyPpc2/cHqijwaYQYAAE8UHi+NXWgub1xknhCMOhFmAADwVH0nSUnDpIoz0n/mWV2NxyLMAADgqWw26bqnpIBQKSpJclRaXZFH4qZ5AAB4stje0n27pDYdrK7EYzEyAwCApyPINIgwAwCAt8jdLq24XTpjt7oSj0KYAQDAGxiG9NY90u5/SR/93upqPAphBgAAb2CzSaN+aS5vfl4qyrW2Hg9CmAEAwFv0ul5KTJUqSnhu01kIMwAAeAubTRr5iLm89S+MzlQhzAAA4E26fk9KHCpVlnHuTBXCDAAA3sRmk0Y8bC5nvSIV51tajifgpnkAAHibpOHSFROkLkOlkHZWV2M5wgwAAN7GZpNuWWZ1FR6Dw0wAAMCrEWYAAPBWleXSttelv4yXKsqsrsYyhBkAALyV4ZDWPiHtXy/tesvqaixDmAEAwFsFBEuDppvLmUvMRx60QoQZAAC82YCfSgEhUv4OKecTq6uxBGEGAABvFhYl9fmRufzZn62txSKEGQAAvN3An5rzXW9Jp45ZW4sFCDMAAHi7S1Kkjv3MRxxsf93qatyOMAMAgC8YNE3qMdoMNa0MdwAGAMAXXPUTc2qFGJkBAABejTADAIAvOXFI2rhIOmO3uhK34TATAAC+5K+3St/ultrGSSlTrK7GLRiZAQDAl/S9zZxvf8PaOtyIMAMAgC/pe5tk85NyNkvH9lldjVsQZgAA8CXh8eYl2pL0+XJra3ETwgwAAL4m+VZzvmtVq3j4JGEGAABf0+s6KSBUOrZXyv/C6mpcjjADAICvCW4rXTpGCmwjffuV1dW4HJdmAwDgi8YulMa3k4LCrK7E5QgzAAD4oohOVlfgNhxmAgDAlxmGdPJbq6twKZeEmYMHD2rq1KlKSkpSaGiounfvrscee0xlZWW12n3xxRe65pprFBISooSEBC1atOi891q5cqV69+6tkJAQJScn65133nFFyQAA+J5j+6Q/DpFeHOHTVzW5JMxkZ2fL4XDohRde0K5du7R48WItXbpU8+bNc7YpKirSmDFj1KVLF2VlZempp57Sr371K7344ovONps2bdKkSZM0depUbdu2TePHj9f48eO1c+dOV5QNAIBviegknciRig5LedutrsZlbIbhnqj21FNP6fnnn9f+/fslSc8//7weeeQR5efnKygoSJI0Z84crV69WtnZ2ZKkiRMn6tSpU1qzZo3zfYYMGaJ+/fpp6dKljf7ZRUVFioyMlN1uV0RERAv2CgAAD7fidmn3v6RhD0rff9Tqapqksd/fbjtnxm63Kyoqyvk6MzNTw4YNcwYZSUpLS9OePXt0/PhxZ5vRo0fXep+0tDRlZmY2+LNKS0tVVFRUawIAoFXqfYM5372m4XZezC1hZu/evVqyZInuvvtu57r8/HzFxcXValf9Oj8/v8E21dvrs2DBAkVGRjqnhISElugGAADe59Ixkl+A+SRtH31WU5PCzJw5c2Sz2Rqcqg8RVTty5IjGjh2rW265RdOmTWvR4uszd+5c2e1255STk+OWnwsAgMcJbS91/Z65nP1va2txkSbdZ+b+++/XnXfe2WCbbt26OZdzc3M1cuRIDR06tNaJvZIUHx+vgoKCWuuqX8fHxzfYpnp7fYKDgxUcHNxgGwAAWo3eN0j7N0jZa6Sr/8fqalpck8JMTEyMYmJiGtX2yJEjGjlypFJSUrRs2TL5+dUeBEpNTdUjjzyi8vJyBQYGSpIyMjLUq1cvtW/f3tlm7dq1mjVrlnO/jIwMpaamNqVsAABat97jzGc09UyzuhKXcMnVTEeOHNGIESPUpUsXvfrqq/L393duqx5Vsdvt6tWrl8aMGaOHH35YO3fu1E9/+lMtXrxY06dPl2Remj18+HAtXLhQ48aN0/Lly/Xkk09q69at6tOnT6Pr4WomAAC8T2O/v13yOIOMjAzt3btXe/fuVefOnWttq85OkZGRev/995Wenq6UlBRFR0dr/vz5ziAjSUOHDtUbb7yhRx99VPPmzVPPnj21evXqJgUZAADg29x2nxkrMTIDAGj1HA4pd6u0b710zWzJz//C+1jM0pEZAADgYYxK6bUJUqld6jZCShhodUUthgdNAgDQGvgHSt1Hmst7M6ytpYURZgAAaC16XmvOv37f2jpaGGEGAIDWokfVI4Jyt0knj1pbSwsizAAA0FqEx0vxV5rL+zdaW0sLIswAANCadBthzvdvsLKKFkWYAQCgNek23JznbrW2jhbEpdkAALQmXa6Wpn4gdbrK6kpaDGEGAIDWJDDUp+4xI3GYCQAAeDnCDAAArU3Jcente6UXR5iPOfByhBkAAFqboLbSjr+b95vJ/8Lqai4aYQYAgNbGP1Dq+j1z+cCH1tbSAggzAAC0Rl2vMeffbLK2jhZAmAEAoDXqMtScH9rk9efNEGYAAGiN4q80z505Y5eOfml1NReFMAMAQGvkHyAlDDaXvfxQEzfNAwCgteoyVLLnmCcEezHCDAAArdX3ZkvDHrC6iovGYSYAAForP9+IAb7RCwAA0HyV5dLp76yuotkIMwAAtGbbXpcWdpHem2d1Jc1GmAEAoDULj5fKT0mHNltdSbMRZgAAaM0uGWDOjx+QThVaW0szEWYAAGjNQttJMb3N5cOfWlpKcxFmAABo7TpXjc7kfGJtHc1EmAEAoLXrPMicMzIDAAC8UkJVmDmyVaqssLaWZuAOwAAAtHbRvaTLfmA+fLKy1HxukxfxrmoBAEDL8/OTJr5mdRXNxmEmAADg1QgzAADAVJQrffWe1VU0GYeZAACAVHJCeuYyc/mhA1JYlKXlNAUjMwAAwLx5XlQ3czlvu5WVNBlhBgAAmDpdZc5zt1lbRxMRZgAAgIkwAwAAvFp1mDlCmAEAAN4o/kpJNqnosHTyqNXVNBphBgAAmEIipOie5nLudktLaQouzQYAADWGPSjJVnPIyQsQZgAAQI0rb7W6gibjMBMAAPBqhBkAAFDbwY+lj5+VTn5rdSWNwmEmAABQ2zsPSEe/lDr0kHpfb3U1F8TIDAAAqC3+SnNesNPaOhqJMAMAAGqLTzbn+V9YW0cjEWYAAEBtzjCzw9o6GokwAwAAaqsOM8cPSmfslpbSGIQZAABQW1iUFNHZXC7YZW0tjUCYAQAA5/OiQ01cmg0AAM43/CHz0Qaxl1ldyQURZgAAwPku6W91BY3GYSYAAODVCDMAAKBuny+X1syWju2zupIGcZgJAADULesV6VCmlDBY6tDd6mrqxcgMAACoW/UVTQWefUUTYQYAANQt9nJzfnS3tXVcAGEGAADUjTADAAC8Wmxvc150RCo5YWkpDSHMAACAuoVE1jzWwINHZ1weZkpLS9WvXz/ZbDZt37691rYvvvhC11xzjUJCQpSQkKBFixadt//KlSvVu3dvhYSEKDk5We+8846rSwYAANWq7wB8bK+1dTTA5WHmoYceUqdOnc5bX1RUpDFjxqhLly7KysrSU089pV/96ld68cUXnW02bdqkSZMmaerUqdq2bZvGjx+v8ePHa+fOna4uGwAASNK4p6WHDkj9b7e6knrZDMMwXPXm7777rmbPnq1//OMfuuKKK7Rt2zb169dPkvT888/rkUceUX5+voKCgiRJc+bM0erVq5WdnS1Jmjhxok6dOqU1a9Y433PIkCHq16+fli5d2ug6ioqKFBkZKbvdroiIiJbrIAAAcJnGfn+7bGSmoKBA06ZN02uvvaawsLDztmdmZmrYsGHOICNJaWlp2rNnj44fP+5sM3r06Fr7paWlKTMzs8GfXVpaqqKioloTAADwTS4JM4Zh6M4779Q999yjAQMG1NkmPz9fcXFxtdZVv87Pz2+wTfX2+ixYsECRkZHOKSEhobldAQCgdTMMKWO+9NoPpeICq6upU5PCzJw5c2Sz2RqcsrOztWTJEhUXF2vu3LmuqrtBc+fOld1ud045OTmW1AEAgNez2aTsf0v71klHv7S6mjo16dlM999/v+68884G23Tr1k3r1q1TZmamgoODa20bMGCAJk+erFdffVXx8fEqKKid8Kpfx8fHO+d1taneXp/g4ODzfjYAAGim2MvMq5mOfil1H2l1NedpUpiJiYlRTEzMBds9++yz+t///V/n69zcXKWlpWnFihUaPHiwJCk1NVWPPPKIysvLFRgYKEnKyMhQr1691L59e2ebtWvXatasWc73ysjIUGpqalPKBgAAFyP2cmn3274xMtNYiYmJtV63bdtWktS9e3d17mzefOfHP/6xHn/8cU2dOlUPP/ywdu7cqT/84Q9avHixc797771Xw4cP1+9+9zuNGzdOy5cv12effVbr8m0AAOBiMVV3AvbQG+dZdgfgyMhIvf/++zpw4IBSUlJ0//33a/78+Zo+fbqzzdChQ/XGG2/oxRdfVN++ffX3v/9dq1evVp8+fawqGwCA1qc6zBR+bZ4Q7GFcep8ZT8F9ZgAAuAgVpdJv4iXDIc3OliI6uuXHWn6fGQAA4CMCgqX2Xc1nNRXnWl3NeVxyzgwAAPAx0zdIwRHmpdoehjADAAAuLCTS6grqxWEmAADg1QgzAADgwuyHpTduk16+zupKzsNhJgAAcGFBbaWv3jWXzxRJIZ5zdTAjMwAA4MJC20ltqx7+fOxrS0s5F2EGAAA0TvSl5vzbr6yt4xyEGQAA0DjVYaZwj7V1nIMwAwAAGiemlzkv5DATAADwRs7DTIzMAAAAbxR9qRQUbt5Az+GwuhonLs0GAACNE9FJmpvjcY80IMwAAIDG8bAQU43DTAAAwKsRZgAAQON9vkJ6brD0n7lWV+LEYSYAANB4jgrp2+yauwF7AEZmAABA43XoYc6/229tHWchzAAAgMbr0N2c23Ok8hJra6lCmAEAAI0X1sG8z4wkfXfA2lqqEGYAAEDj2WxSVNXozLG91tZShTADAACaxnnezD5r66jC1UwAAKBp4q6QjiZLweFWVyKJMAMAAJrqe7PMyUNwmAkAAHg1wgwAAGgeh8Mjnp5NmAEAAE336o3Sk52k/M+troQwAwAAmqGiTKoo8Yg7ARNmAABA07Xvas6PH7SyCkmEGQAA0BxRSebcA+4CTJgBAABN174qzDAyAwAAvBIjMwAAwKtVnzNTdESqKLW0FO4ADAAAmq5NjHRJihTeUSo9KQUEW1YKYQYAADSdzSZNW2d1FZI4zAQAAC5Cnr1Em/YVKs9eYlkNjMwAAIBmWfHpIc1d9YWCjDKV2YK1YEKyJg5MdHsdjMwAAIAmy7OX6KPVL+qLoKn6v8Bn5TCkeat2WjJCQ5gBAABNdqDwlE4YbdTWdkaJtqOSpErD0MHC026vhTADAACaLCm6jXKMOElSou2obHLI32ZT1+gwt9fCOTMAAKDJOkaG6ufjh+ult6/XQSNOwTZDj0/oo46RoW6vxWYYhuH2n+pmRUVFioyMlN1uV0REhNXlAADgM/LsJTpYeFpdo8NaPMg09vubkRkAANBsHSNDLRmNORvnzAAAAK9GmAEAAF6NMAMAALwaYQYAAHg1wgwAAPBqhBkAAODVCDMAAMCrEWYAAIBXI8wAAACvRpgBAABejTADAAC8GmEGAAB4NcIMAADwaq3iqdmGYUgyHyUOAAC8Q/X3dvX3eH1aRZgpLi6WJCUkJFhcCQAAaKri4mJFRkbWu91mXCju+ACHw6Hc3FyFh4fLZrO12PsWFRUpISFBOTk5ioiIaLH39SS+3kf65/18vY/0z/v5eh9d2T/DMFRcXKxOnTrJz6/+M2NaxciMn5+fOnfu7LL3j4iI8Mm/oGfz9T7SP+/n632kf97P1/voqv41NCJTjROAAQCAVyPMAAAAr0aYuQjBwcF67LHHFBwcbHUpLuPrfaR/3s/X+0j/vJ+v99ET+tcqTgAGAAC+i5EZAADg1QgzAADAqxFmAACAVyPMAAAAr0aYaYKDBw9q6tSpSkpKUmhoqLp3767HHntMZWVlDe535swZpaenq0OHDmrbtq1uvvlmFRQUuKnqpvnNb36joUOHKiwsTO3atWvUPnfeeadsNlutaezYsa4t9CI0p4+GYWj+/Pnq2LGjQkNDNXr0aH399deuLbSZvvvuO02ePFkRERFq166dpk6dqpMnTza4z4gRI877DO+55x43VXxhzz33nLp27aqQkBANHjxYn3zySYPtV65cqd69eyskJETJycl655133FRp8zSlf6+88sp5n1VISIgbq22aDz/8UDfeeKM6deokm82m1atXX3CfDRs2qH///goODlaPHj30yiuvuLzO5mpq/zZs2HDe52ez2ZSfn++egptowYIFGjhwoMLDwxUbG6vx48drz549F9zP3b+DhJkmyM7OlsPh0AsvvKBdu3Zp8eLFWrp0qebNm9fgfvfdd5/efvttrVy5Uhs3blRubq4mTJjgpqqbpqysTLfccotmzJjRpP3Gjh2rvLw85/Tmm2+6qMKL15w+Llq0SM8++6yWLl2qLVu2qE2bNkpLS9OZM2dcWGnzTJ48Wbt27VJGRobWrFmjDz/8UNOnT7/gftOmTav1GS5atMgN1V7YihUrNHv2bD322GPaunWr+vbtq7S0NB09erTO9ps2bdKkSZM0depUbdu2TePHj9f48eO1c+dON1feOE3tn2TeafXsz+qbb75xY8VNc+rUKfXt21fPPfdco9ofOHBA48aN08iRI7V9+3bNmjVLP/vZz/Tee++5uNLmaWr/qu3Zs6fWZxgbG+uiCi/Oxo0blZ6ers2bNysjI0Pl5eUaM2aMTp06Ve8+lvwOGrgoixYtMpKSkurdfuLECSMwMNBYuXKlc93u3bsNSUZmZqY7SmyWZcuWGZGRkY1qO2XKFOOmm25yaT2u0Ng+OhwOIz4+3njqqaec606cOGEEBwcbb775pgsrbLovv/zSkGR8+umnznXvvvuuYbPZjCNHjtS73/Dhw417773XDRU23aBBg4z09HTn68rKSqNTp07GggUL6mx/6623GuPGjau1bvDgwcbdd9/t0jqbq6n9a8rvpqeRZLz11lsNtnnooYeMK664ota6iRMnGmlpaS6srGU0pn/r1683JBnHjx93S00t7ejRo4YkY+PGjfW2seJ3kJGZi2S32xUVFVXv9qysLJWXl2v06NHOdb1791ZiYqIyMzPdUaJbbNiwQbGxserVq5dmzJihY8eOWV1Sizlw4IDy8/NrfYaRkZEaPHiwx32GmZmZateunQYMGOBcN3r0aPn5+WnLli0N7vvXv/5V0dHR6tOnj+bOnavTp0+7utwLKisrU1ZWVq0/ez8/P40ePbreP/vMzMxa7SUpLS3N4z4rqXn9k6STJ0+qS5cuSkhI0E033aRdu3a5o1y38KbP72L069dPHTt21LXXXquPP/7Y6nIazW63S1KD33tWfIat4kGTrrJ3714tWbJETz/9dL1t8vPzFRQUdN65GXFxcR57jLSpxo4dqwkTJigpKUn79u3TvHnzdN111ykzM1P+/v5Wl3fRqj+nuLi4Wus98TPMz88/b7g6ICBAUVFRDdb64x//WF26dFGnTp30xRdf6OGHH9aePXu0atUqV5fcoMLCQlVWVtb5Z5+dnV3nPvn5+V7xWUnN61+vXr308ssv68orr5TdbtfTTz+toUOHateuXS59oK671Pf5FRUVqaSkRKGhoRZV1jI6duyopUuXasCAASotLdVLL72kESNGaMuWLerfv7/V5TXI4XBo1qxZuvrqq9WnT59621nxO8jIjKQ5c+bUeULW2dO5/7AcOXJEY8eO1S233KJp06ZZVHnjNKd/TXHbbbfpBz/4gZKTkzV+/HitWbNGn376qTZs2NBynbgAV/fRaq7u3/Tp05WWlqbk5GRNnjxZf/nLX/TWW29p3759LdgLtITU1FTdcccd6tevn4YPH65Vq1YpJiZGL7zwgtWloRF69eqlu+++WykpKRo6dKhefvllDR06VIsXL7a6tAtKT0/Xzp07tXz5cqtLOQ8jM5Luv/9+3XnnnQ226datm3M5NzdXI0eO1NChQ/Xiiy82uF98fLzKysp04sSJWqMzBQUFio+Pv5iyG62p/btY3bp1U3R0tPbu3atRo0a12Ps2xJV9rP6cCgoK1LFjR+f6goIC9evXr1nv2VSN7V98fPx5J45WVFTou+++a9Lft8GDB0syRx+7d+/e5HpbSnR0tPz9/c+7+q+h35/4+PgmtbdSc/p3rsDAQF111VXau3evK0p0u/o+v4iICK8flanPoEGD9NFHH1ldRoNmzpzpvKDgQiOAVvwOEmYkxcTEKCYmplFtjxw5opEjRyolJUXLli2Tn1/Dg1spKSkKDAzU2rVrdfPNN0syz2I/dOiQUlNTL7r2xmhK/1rC4cOHdezYsVpf/K7myj4mJSUpPj5ea9eudYaXoqIibdmypclXfTVXY/uXmpqqEydOKCsrSykpKZKkdevWyeFwOANKY2zfvl2S3PoZ1iUoKEgpKSlau3atxo8fL8kc6l67dq1mzpxZ5z6pqalau3atZs2a5VyXkZHhtt+3pmhO/85VWVmpHTt26Prrr3dhpe6Tmpp63mW8nvr5tZTt27db/rtWH8Mw9Itf/EJvvfWWNmzYoKSkpAvuY8nvoMtOLfZBhw8fNnr06GGMGjXKOHz4sJGXl+eczm7Tq1cvY8uWLc5199xzj5GYmGisW7fO+Oyzz4zU1FQjNTXVii5c0DfffGNs27bNePzxx422bdsa27ZtM7Zt22YUFxc72/Tq1ctYtWqVYRiGUVxcbDzwwANGZmamceDAAeODDz4w+vfvb/Ts2dM4c+aMVd1oUFP7aBiGsXDhQqNdu3bGP//5T+OLL74wbrrpJiMpKckoKSmxogsNGjt2rHHVVVcZW7ZsMT766COjZ8+exqRJk5zbz/07unfvXuOJJ54wPvvsM+PAgQPGP//5T6Nbt27GsGHDrOpCLcuXLzeCg4ONV155xfjyyy+N6dOnG+3atTPy8/MNwzCM22+/3ZgzZ46z/ccff2wEBAQYTz/9tLF7927jscceMwIDA40dO3ZY1YUGNbV/jz/+uPHee+8Z+/btM7KysozbbrvNCAkJMXbt2mVVFxpUXFzs/B2TZDzzzDPGtm3bjG+++cYwDMOYM2eOcfvttzvb79+/3wgLCzMefPBBY/fu3cZzzz1n+Pv7G//5z3+s6kKDmtq/xYsXG6tXrza+/vprY8eOHca9995r+Pn5GR988IFVXWjQjBkzjMjISGPDhg21vvNOnz7tbOMJv4OEmSZYtmyZIanOqdqBAwcMScb69eud60pKSoyf//znRvv27Y2wsDDjhz/8Ya0A5EmmTJlSZ//O7o8kY9myZYZhGMbp06eNMWPGGDExMUZgYKDRpUsXY9q0ac5/iD1RU/toGObl2b/85S+NuLg4Izg42Bg1apSxZ88e9xffCMeOHTMmTZpktG3b1oiIiDDuuuuuWkHt3L+jhw4dMoYNG2ZERUUZwcHBRo8ePYwHH3zQsNvtFvXgfEuWLDESExONoKAgY9CgQcbmzZud24YPH25MmTKlVvu//e1vxqWXXmoEBQUZV1xxhfHvf//bzRU3TVP6N2vWLGfbuLg44/rrrze2bt1qQdWNU30p8rlTdZ+mTJliDB8+/Lx9+vXrZwQFBRndunWr9bvoaZrav9/+9rdG9+7djZCQECMqKsoYMWKEsW7dOmuKb4T6vvPO/kw84XfQVlUsAACAV+JqJgAA4NUIMwAAwKsRZgAAgFcjzAAAAK9GmAEAAF6NMAMAALwaYQYAAHg1wgwAAPBqhBkAAODVCDMAAMCrEWYAAIBXI8wAAACv9v/K+AiILOQdOAAAAABJRU5ErkJggg==",
      "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 = 10   # order of polynomial\n",
    "n = 10   # 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",
    "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 = 37411\n",
      "np.linalg.cond gives 37411\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": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Fitted coeffs =  [-0.09322416 -0.43254686  0.08412339  0.60685865  0.05847538 -0.49706461\n",
      "  0.01186025 -0.17740619 -0.05571574 -0.84552489 -0.04429305]\n",
      "frac error =  [ 1.39182716e-01 -1.05235175e-14 -1.44526222e+00  4.20775569e-15\n",
      " -9.39232614e-01  2.98180099e-14 -1.01292238e+00 -7.32195955e-14\n",
      " -1.18687032e+00  3.67656176e-15 -4.95076272e-01]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjMAAAGdCAYAAADnrPLBAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA6j0lEQVR4nO3deXwU9eH/8ffm2hCSbAjkAEkgHIJIBAlX0IIHJSgeFGtRqYWWrwcGK0pVUIvaC4sH3gWrglZRxJ/aar2Qy4NAIYBAEJQbEhLubICQhGR+fwzZJJCEBLI7e7yej8c8dnb3s5v3sIR985nZWZthGIYAAAB8VJDVAQAAAM4FZQYAAPg0ygwAAPBplBkAAODTKDMAAMCnUWYAAIBPo8wAAACfRpkBAAA+LcTqAJ5QUVGhvLw8RUVFyWazWR0HAAA0gGEYKioqUps2bRQUVPf8S0CUmby8PCUlJVkdAwAAnIVdu3apbdu2dd4fEGUmKipKkvmHER0dbXEaAADQEE6nU0lJSa738boERJmp3LUUHR1NmQEAwMec6RARDgAGAAA+jTIDAAB8GmUGAAD4NMoMAADwaZQZAADg0ygzAADAp1FmAACAT6PMAAAAn0aZAQAAPo0yAwAAfBpl5hzsKSzW0i37taew2OooAAAErID4biZ3mLtipyZ/sE4VhhRkk6aOSNXIPslWxwIAIOAwM3MW9hQWu4qMJFUY0kMfrGeGBgAAC3iszDzxxBOy2WyaMGGC67bjx48rMzNTLVu2VGRkpG644QYVFBTUeNzOnTs1bNgwRUREKD4+Xvfff79OnDjhqdi12rb/qCoMaXjQt/pLyGvqZNutcsPQ9v3HLM0FAEAg8kiZWbFihWbOnKmLLrqoxu333nuvPv74Y82bN09LlixRXl6eRowY4bq/vLxcw4YNU2lpqZYuXao33nhDs2fP1pQpUzwRu04prZoryCb9MniJfh2yQGlBPynYZlP7VhGW5gIAIBC5vcwcOXJEo0aN0j//+U+1aNHCdXthYaFee+01PfPMM7riiiuUlpamWbNmaenSpVq2bJkk6csvv9SGDRv01ltvqWfPnrrqqqv05z//WS+99JJKS0vdHb1OrR3NNHVEqjYa7SVJFwbt0N9GdFdrRzPLMgEAEKjcXmYyMzM1bNgwDR48uMbt2dnZKisrq3F7165dlZycrKysLElSVlaWUlNTlZCQ4BqTkZEhp9OpnJycOn9mSUmJnE5njaWpjeyTrF8OGypJuimpkIN/AQCwiFs/zfTuu+9q1apVWrFixWn35efnKywsTDExMTVuT0hIUH5+vmtM9SJTeX/lfXWZOnWqHn/88XNMf2YxKb0kSWEHfpAMQ7LZ3P4zAQBATW6bmdm1a5fuuecevf322woPD3fXj6nV5MmTVVhY6Fp27drlnh/U6nwpKFQqcUqHd7jnZwAAgHq5rcxkZ2dr79696tWrl0JCQhQSEqIlS5bo+eefV0hIiBISElRaWqrDhw/XeFxBQYESExMlSYmJiad9uqnyeuWY2tjtdkVHR9dY3CIkTIrvaq7nr3fPzwAAAPVyW5m58sortW7dOq1Zs8a19O7dW6NGjXKth4aGasGCBa7HbNq0STt37lR6erokKT09XevWrdPevXtdY+bPn6/o6Gh169bNXdEbJyHVvDy41docAAAEKLcdMxMVFaXu3bvXuK158+Zq2bKl6/axY8fqvvvuU2xsrKKjo3X33XcrPT1d/fv3lyQNGTJE3bp106233qpp06YpPz9fjzzyiDIzM2W3290VvXGu/KOU8VcpItbqJAAABCRLv85g+vTpCgoK0g033KCSkhJlZGTo5Zdfdt0fHBysTz75ROPGjVN6erqaN2+u0aNH609/+pOFqU8R3cbqBAAABDSbYRiG1SHczel0yuFwqLCw0H3HzwAAgCbV0PdvvpupKXzztPTm9dLubKuTAAAQcCgzTWHncmnrYimXMgMAgKdRZppC4skDnQvWWZsDAIAARJlpCoknP57NuWYAAPA4ykxTSDz5beAFOVJ5mbVZAAAIMJSZptAiRbI7pPISad9Gq9MAABBQKDNNIShIatPDXM9bbW0WAAACDGWmqbS5WGoWK5UeszoJAAABhZPmNZWy41KIXbLZ3PP8AAAEmIa+f1v6dQZ+JTTc6gQAAAQkdjO5g/9PdgEA4DUoM01pyZPS9FQpe7bVSQAACBiUmaZUdkwq3MknmgAA8CDKTFNqc7F5SZkBAMBjKDNNqU1P83LvBvPTTQAAwO0oM03JkSRFtJQqTphfbQAAANyOMtOUbLZqu5pWWZsFAIAAQZlpapVlJpcyAwCAJ1BmmlpSf7PQtOpsdRIAAAICZwBuap0HmwsAAPAIZmYAAIBPo8y4S+kx6fAuq1MAAOD3KDPukPOR9ESS9O+7rE4CAIDfo8y4Q8tO5rlmcldJFeVWpwEAwK9RZtwh/gIpLEoqPWKeDRgAALgNZcYdgoKl83qZ67v+Z20WAAD8HGXGXZL6mpe7V1ibAwAAP0eZcZe2J8sMMzMAALgVZcZd2vY2Lw9ukY4esDYLAAB+jDMAu0tErNRvnBSTLAXRGQEAcBfKjDtd9YTVCQAA8HtMGQAAAJ9GmXG3A1ukVW9KJUVWJwEAwC+xm8nd/jVcOrxTim4jdeLbtAEAaGrMzLhbu0vNy+3fWZsDAAA/RZlxt/aXmJc7KDMAALgDZcbd2p0sM7mrpNJj1mYBAMAPUWbcrUV7Kfo8qaJM2s3ZgAEAaGqUGXez2aR2A8z1HUutzQIAgB+izHhC5a4mDgIGAKDJ8dFsT+hylRQZLyWnW50EAAC/Q5nxhKhEqeswq1MAAOCX2M0EAAB8GmXGU47skxZNlT4cZ3USAAD8CmXGk5Y8IX0/Rzqy1+okAAD4DcqMp0TGSYmp5vrWJdZmAQDAj1BmPKnD5ebl1kXW5gAAwI9QZjyp48kys2WRZBjWZgEAwE9QZjwpOV0KtktFedL+H61OAwCAX6DMeFJoM6ndyRPnbVlobRYAAPwEZcbTOl4h2R3SiRKrkwAA4Bc4A7Cn9fk/qf9dUnCo1UkAAPALlBlPC2tudQIAAPwKu5msYhhSUYHVKQAA8HmUGSvs/0l67iJp5kCposLqNAAA+DTKjBVikqVjB6Uj+VL+91anAQDAp1FmrBBirzqB3qbPrc0CAICPo8xY5fyh5uWPlBkAAM4FZcYqnYdIskl71kjOPVanAQDAZ1FmrBIZL52XZq7/+Jm1WQAA8GGUGSt1HWZe/vCxtTkAAPBhbi0zU6dOVZ8+fRQVFaX4+HgNHz5cmzZtqjHm+PHjyszMVMuWLRUZGakbbrhBBQU1z7+yc+dODRs2TBEREYqPj9f999+vEydOuDO6Z1w4XEr7rXTpvVYnAQDAZ7m1zCxZskSZmZlatmyZ5s+fr7KyMg0ZMkRHjx51jbn33nv18ccfa968eVqyZIny8vI0YsQI1/3l5eUaNmyYSktLtXTpUr3xxhuaPXu2pkyZ4s7onhHbQbr2WSlloNVJAADwWTbDMAxP/bB9+/YpPj5eS5Ys0cCBA1VYWKi4uDjNmTNHv/zlLyVJGzdu1AUXXKCsrCz1799fn332ma655hrl5eUpISFBkjRjxgw9+OCD2rdvn8LCws74c51OpxwOhwoLCxUdHe3WbQQAAE2joe/fHj1mprCwUJIUGxsrScrOzlZZWZkGDx7sGtO1a1clJycrKytLkpSVlaXU1FRXkZGkjIwMOZ1O5eTk1PpzSkpK5HQ6ayxeyzCkncukzyebJ9IDAACN4rEyU1FRoQkTJuiSSy5R9+7dJUn5+fkKCwtTTExMjbEJCQnKz893jaleZCrvr7yvNlOnTpXD4XAtSUlJTbw1Tchmkz65T1r2srSJTzUBANBYHiszmZmZWr9+vd599123/6zJkyersLDQtezatcvtP/OcdLvevNzwb2tzAADggzxSZsaPH69PPvlEixYtUtu2bV23JyYmqrS0VIcPH64xvqCgQImJia4xp366qfJ65ZhT2e12RUdH11i82oXDzcstC6Xiw1YmAQDA57i1zBiGofHjx+vDDz/UwoULlZKSUuP+tLQ0hYaGasGCBa7bNm3apJ07dyo9PV2SlJ6ernXr1mnv3r2uMfPnz1d0dLS6devmzvieE9dFiu8mVZQxOwMAQCO5tcxkZmbqrbfe0pw5cxQVFaX8/Hzl5+eruLhYkuRwODR27Fjdd999WrRokbKzs/Xb3/5W6enp6t+/vyRpyJAh6tatm2699VZ9//33+uKLL/TII48oMzNTdrvdnfE966KR5uX37t8NBwCAP3HrR7NtNlutt8+aNUtjxoyRZJ40b+LEiXrnnXdUUlKijIwMvfzyyzV2Ie3YsUPjxo3T4sWL1bx5c40ePVpPPPGEQkJCGpTDJz6a7cyTnukmyZB+v0aKTTnTIwAA8GsNff/26HlmrOITZUaS3rxe2vejNPxlqePlVqcBAMBSDX3/btjUBjxjxKtSRKwUFGx1EgAAfAZlxptExlmdAAAAn8O3Znuj8hPS3o1WpwAAwCdQZrzNgS3SMxdIrw2RSo9ZnQYAAK9HmfE2LVKk0HCppFDK+dDqNAAAeD3KjLcJCpLSxpjr2bMsjQIAgC+gzHiji2+VgkKk3Suk/HVWpwEAwKtRZrxRZLx0wbXm+kpmZwAAqA9lxlul/da8XPueVHLE2iwAAHgxyoy3ShkotewklRZJWxdbnQYAAK/FSfO8lc0mDf271LyV1Kan1WkAAPBalBlv1nmw1QkAAPB67GbyFcWHrU4AAIBXosx4O8OQPp4gPXW+tGet1WkAAPA6lBlvZ7NJxwul8hJp2ctWpwEAwOtQZnzBgPHm5br3JWeetVkAAPAylBlfcF6a1O4SqaJM+u45q9MAAOBVKDO+YuD95uXKWZJzj7VZAADwIpQZX9HhMimpv3nsDLMzAAC4UGZ8hc0mXTbJXF/3nlRWbG0eAAC8BCfN8yUdLpOuelLqPkIKbWZ1GgAAvAJlxpfYbFK/261OAQCAV2E3k68yDCl/ndUpAACwHGXGF5WXSW+NkGZcKuWttjoNAACWosz4ouBQqXm8uf7lH81ZGgAAAhRlxldd8YgUbJe2fyNt+szqNAAAWIYy46tikqT0u8z1Lx7io9oAgIBFmfFlP5soRbWRDm2TvnnG6jQAAFiCMuPL7FHSVX8317+dLu370do8AABYgDLj6y64VuqcIbVoJ5UUWZ0GAACP46R5vs5mk4a/bM7ShNitTgMAgMdRZvxB81Y1rxuGWXIAAAgA7GbyJxUV0tIXzRPqVZRbnQYAAI+gzPiTojxp0d+kLQulZS9bnQYAAI+gzPgTR1tp6N/M9QV/4qsOAAABgTLjb3qNlroMk8pLpfd+Ix07aHUiAADcijLjbyo/3dSivXR4p/ThneaxNAAA+CnKjD9qFiP96k3zu5t++kL65mmrEwEA4DaUGX/Vuoc07CnJFiw1b2l1GgAA3IbzzPizXr+RkvpJcV2sTgIAgNswM+PvqheZowekQ9stiwIAgDtQZgLFgS3Sa4Ol2ddKh3dZnQYAgCZDmQkUYZGSbFLhTunN6yTnHqsTAQDQJCgzgSIqQRr9HykmWTq4lUIDAPAblJlA4mgrjf5Yij5P2v+j9HqGufsJAAAfRpkJNC3aS7/9VGqRIh3eIb0+VPu2rNLSLfu1p7DY6nQAADQaH80ORC3aS7/7QnrrBjmPFOnqVzdqn7FHQTZp6ohUjeyTbHVCAAAajJmZQBWVoPwR7+uqgxO1z3BIkioMQw99sI4ZGgCAT2FmJoBtPRKiXKPq7MC/Dv5Klwat1+7crmrt6GhhMgAAGo4yE8BSWjVXkE2qMKRWKtRDIXMUYStR+afXSBH/lNpfanVEAADOiN1MAay1o5mmjkhVsM2m/XLo5rIpKmreTsFH8qTZ10jzH5VKj1kdEwCAetkMwzCsDuFuTqdTDodDhYWFio6OtjqO19lTWKzt+4+pfasItQ4vlz57UFrzlnlni/bSsKelToMtzQgACDwNff9mZgZq7Wim9I4t1drRTLJHSsNfkm56xzwfzaHt0ju3SEUFVscEAKBWHDOD2nW9Wkr5mbTob1JES/MMwpWOO6VwZrgAAN6BMoO62aOkoVNr3rYjS3r7l1LaGKn/XZLjPEuiAQBQid1MaJzv35FKj0hZL0rP9ZA+ypRyV1mdCgAQwCgzaJxrn5NGvS+1u0SqKDMPFP7n5dLMgdLKWVL5CasTAgACDGUGjWOzSZ1/bn6/0+++lLr/UgoOk/Z8Ly2fIQUFV40tL7MuJwAgYHDMDM5ecj9zOXpA+n6OFJlolh1JKjkiPdtdOi/N/Fh3xyulVp2r7gcAoIlQZnDumreUBtxd87ati6TiQ9Lmr8xFkhzJUrt0qW0fc3anRXuPRwUA+B92M8E9ul4jjcuShvxF6nCZuSuqcKe0dq706R+krUuqxu7fbB5vs+t/UkmRZZEBAL7JZ2ZmXnrpJT355JPKz89Xjx499MILL6hv375Wx0JdbDYpoZu5DLhbKj0q7VxmFpZdy6Xk9KqxWxZKn91fdT26rTlrE9vevEz9ldSinXmfYdS5q2pPYbG27T+qlFbNzRMAAgACgk+Umblz5+q+++7TjBkz1K9fPz377LPKyMjQpk2bFB8fb3U8NERYc6nTleZyqqgE85iaghzpSL7k3G0uO741708ZVFVmls+QlkyTohKlyARziWipdQel93OO6NPyfjpgi9HUEakameowS1RohLkEh/rEMTv+VMr8ZVvYDsC7+cR3M/Xr1099+vTRiy++KEmqqKhQUlKS7r77bk2aNOmMj+e7mXzIsYPSgS3SoW3mVykc2i79/M/mcTmSNH+K9N1zdT782pK/aJ3RQcE2m1YP3qDob/5cdact2CxVoc3M5ZevmwcoS9IPn0ir3pCCQs3SExxac73vHVLc+ebY3FXST19KtqCT5ch2cv3k0nWY1LKjOXbfj9K2Jea46mMqH9P+kqpjhw7v0rJF/9HclbtdE1Aj+ySpf4dW5v1te0uxHcx15x5px3en/wFUlrXWPasyHN0vbf+m+qCaY+O7mQdnS+ZxTtu/PfVJq1ZbdZbiupjrJUX1jv0s167M+cdUYUjNbKWacekxDTo/vtrPr/aY6POk+K7m1ROl0s6ltWeVTYqMr8pQUW7O9NU1NiK2atskaffKU8ZWyxDuqPozk6T8dZJRoS9yCvT8wp9UYdhks0m/v7Kzhl7cseq1kMzXueJEzZ9duR4SXlXGJfPvdOUpDKqXa5vN/DsXk1R1mzOv6lOBpz63LViKbl019uj+esYGae6GY5r8wTpVGJLDdkyPXXO+ftGz7ekF3xYkNYupul56zNy2U5+38npoRNXtJ0olo6Lm2KBQKahpj2iglAWOhr5/e/3MTGlpqbKzszV58mTXbUFBQRo8eLCysrJqfUxJSYlKSkpc151Op9tzoolExJpLUp/a7//ZROmikVJRvnRkr3QkX7l78rTk+5/ksB3VXiNGklRuGCp0FinaFiwZ5eZjjXKpxGkuklS9xh/cahaUulxwXVWZyVstLZ5a99hWnaveFHOzzWOE6vLL111l5tDmZer//UPqH1rt/jUnF0m67sWqN9CC9dL/G1v38141rVqh2iTNG1P32MGPSZfea64f2i7N/XXdYwc+IF3xsLnuzJPeuanOobknrlaFYT5XC8OpQSt+L62oY3Dv30nXTDfXS4qkN6+vO8NFN0kjZprr5aXSrKvqHnvBddLIf1Vdf3Wwar7w1XQaLP36/1Vdfy1DKjuqDEkZYdXGfSOVbO0v+21fVN02+2rp6L7an7d1T+mOaseIvXGddHhH7WNbnS+Nr/aH9K8R0r4fah8b3Va6L6fq+pxfmX/falERHqvJhS+q4uSm/yPkGQ2Yv0GaX8vgkHDpkWrfxTZvjPTTF7UMPGnKoari8uHtUs6Hp48JCjGPm7vvh6qitPCv0sZPzNtD7FWXYc2lsChpyJ/NfwskadcK83fUHqmF24/ria/3a39FtJy2SP11RA+N7JNcdz4EBK8vM/v371d5ebkSEhJq3J6QkKCNGzfW+pipU6fq8ccf90Q8eFq4w1wSLnTdFFRYrEdWLXT9Qy1JwTabQq6YJP3icfN/q6VHpbJiqeyYuZQeq/rfvWS+kUXEmmPLy8wTAlZfr/7Jq7iuUu+x5v9AXYtRtR7dpmqso63U7frTx1Rej6z6e727LErrylNrbK5NhlLPcygmIlSKqva/8GaxUsrAquunTrBWHxseLbW7tHLg6eOjq30lRWhzKal/3WOrf31FiF06r3etY4tKTigvv6VraJlCtLYiRR1aNVekPbjac568rP5nFhQkxV9Yy88/eRmVWG1DbVLLznWPjaz574ZatKuaOXBt1smViFY1x0YlqqT4iA4dK5Xt5JjKeYkQW5Ts1cdGtDz5s43TM9ijaj6vPUqyO2ofG3rKLEOIXQppdsr2nVwPPuWfb1uwuVQ+X7XnLjeMGr8fHldxwlyCq7XCwl3S3g11P+bn1f4N//4daeVrkqQrJF1x8mnKDZsOfxKlvfHzFd/u5Mze9u+k/ZskR5K5xCSZBQl+zet3M+Xl5em8887T0qVLlZ5eddDoAw88oCVLlmj58uWnPaa2mZmkpCR2M/mxuSt26qEP1qvcMBRss+lvI7r73P/W9hQW65InTi9l30663Oem0v1lW/x3O8zfk28fvEyto2spS8HVpgfLy6rK98nHVo2VWcAqZ2bKis3SUn1s+QmpvMScRYtpVzV234+SM/fkfxpKpBMnl9KjUmmRlH63FHKytSz7h/Tj5ypyHtbevQVqaXMqxnbUFXHFr7LVp1sn88p/J0orXq35BxDVxpxZjesq/ewPUmTcWf9ZwrP8ZjdTq1atFBwcrIKCghq3FxQUKDExsdbH2O122e32Wu+DfxrZJ1kDz4/T9v3H1L5VhE+90VRq7WimqSNSTytlbIt1/Hc7gsztiGnAjEX1YnMmp84s1Sfu/Kpdt2fSf5zUf5yOFBbr5ydLWYhOqIWKFBdUpNfaVJuJjO8mnT9UKtwtHd4llRRKRXnmsnWxdFm14yyXzTCPjTqvl5TU15wRbOLje+AZXj8zI5kHAPft21cvvPCCJPMA4OTkZI0fP54DgOF39hQW+3Qpq85ftoXt8B6NnoUtPmSey2rfRunwzqpjviRp9jU1D46PaCV1GGSeG6vD5TUPxoYlGvr+7RNlZu7cuRo9erRmzpypvn376tlnn9V7772njRs3nnYsTW0oMwDgP5qslG3+yjz/1e6V5jmwyqp2XSksUnpgq3ncEizjN7uZJGnkyJHat2+fpkyZovz8fPXs2VOff/55g4oMAMC/tHY0a5qZpU6DzUUyP1aeu9LcFbVlkRSTXLPIfDhOapsmpd5ofggBXsUnZmbOFTMzAIBGqaioOn5mz/fSzJOfHgxpJl34C6n/nVLrHtblCxANff/mSCcAAE5V/UDgmGQpY6r5aagTxdL3c8xy8+ZwcxbH/+cEvB4zMwAANIRhmMfX/G+mtP6DqhNy/upN83xSaHJ+dcwMAACWs9nMs5Mn9ZGu+KOU9ZK0M0vqMqxqTHlZ4z7OjiZBmQEAoLFatJOunmZ+P1jQyTMvlx2X/nmFdME10qX3SaHh1mYMIBwzAwDA2aosMpK04SNpb4605O/SjEtq+SJWuAtlBgCApnDRSOnG2eZ3gh3YbJ6U76vHqr7NHG5DmQEAoCnYbObHtsevkHr9RpIhfTtdem2IdHCb1en8GmUGAICmFO6QrnvB/JRTeIyUt8r8Aky4DQcAAwDgDt2ul85Lkz59QBr2lNVp/BozMwAAuIujrXTzHCm6TdVt2781zzCMJkOZAQDAU75/V5o9TProTu056NTSLfu1p7DY6lQ+j91MAAB4jE2yBUtr5ypn9U+6q+wenbCFauqIVI3sk2x1OJ/FzAwAAJ7SY6QOXjdLx41QDQ5epedDX5TNKNdDH6xnhuYcUGYAAPCgjdGX6P/K/qASI0RDg1fo6dB/yDDKtX3/Mauj+SzKDAAAHpTSqrmWGqm6q+welRnBGh68VH8JnaX2LZtZHc1nUWYAAPCg1o5mmjoiVYuN3rqnLFPlhk1pnc5T62i+y+lscQAwAAAeNrJPsgaeH6ft+/vpoH6hLh17WR3Jp1FmAACwQGtHM7V2NJPUsurGE6WSM1eKTbEsly9iNxMAAN7g2EHpX8PN5Xih1Wl8CmUGAABvYLNJhbukQ9ulf2dKhmF1Ip9BmQEAwBs0ayHdOFsKCpV++FhaPtPqRD6DMgMAgLc4L00a8hdz/ctHpD1rrc3jIygzAAB4k353SF2vkSrKpI/GmQcFo16UGQAAvInNJl3zrBTRUipYL339pNWJvB5lBgAAbxMZJw17WkpMlbpdb3Uar8d5ZgAA8EYX/kLqeq0UzFv1mTAzAwCAt6peZEqOWJfDy1FmAADwZidKpQV/kp5NlZx5VqfxSpQZAAC8WVCItO0bqfig9MXDVqfxSpQZAAC8WVCQeTCwLUjK+UDa/p3VibwOZQYAAG/X+iIpbYy5/uUjfNXBKSgzAAD4gssmS2GRUt4qc4YGLpQZAAB8QWS8dMk95vpXj0snSqzN40UoMwAA+Ir0TCkyUTpSIO1eaXUar8GZeAAA8BVhzaVfvSG1aC9FJVqdxmtQZgAA8CXJ/a1O4HXYzQQAgK/anS2VFVudwnKUGQAAfNF/fi+9eoWUPdvqJJajzAAA4Iva9DQvl74olZdZGsVqlBkAAHxRj1uk5vGSc7e0PrDPO0OZAQDAF4WGS/3uMNeXPh/QZwWmzAAA4Kt6/04KbS4VrJe2LLQ6jWUoMwAA+KqIWKnXreb60uetzWIhygwAAL6s/13mN2rv3SgVH7Y6jSU4aR4AAL6sRTtpzH+ltn2k4FCr01iCMgMAgK9rN8DqBJZiNxMAAP6iolw6vMvqFB5HmQEAwB/kZkvP9ZDmjAy4j2lTZgAA8AexHaSj+6W9OdKu5Van8SjKDAAA/qBZC+nCX5jrq9+yNouHUWYAAPAXF//avMz5UCo9am0WD6LMAADgL9oNkFqkSKVHpA3/sTqNx1BmAADwFzabdPEocz2AdjVRZgAA8Cc9bpZkk3Z8Kx3aYXUaj+CkeQAA+BNHWynjb+YZgWOSrU7jEZQZAAD8TfpdVifwKHYzAQAAn0aZAQDAHxXkSB/fI339pNVJ3I4yAwCAPzqwRcqeLa2cJVVUWJ3GrSgzAAD4o85DJHu05MyVdi2zOo1buaXMbN++XWPHjlVKSoqaNWumjh076tFHH1VpaWmNcWvXrtXPfvYzhYeHKykpSdOmTTvtuebNm6euXbsqPDxcqamp+vTTT90RGQAA/xIaLl1wrbm+7n1rs7iZW8rMxo0bVVFRoZkzZyonJ0fTp0/XjBkz9NBDD7nGOJ1ODRkyRO3atVN2draefPJJPfbYY3rllVdcY5YuXaqbb75ZY8eO1erVqzV8+HANHz5c69evd0dsAAD8S/cR5uXGT/x6V5PNMDzzPeFPPvmk/vGPf2jr1q2SpH/84x96+OGHlZ+fr7CwMEnSpEmT9NFHH2njxo2SpJEjR+ro0aP65JNPXM/Tv39/9ezZUzNmzGjwz3Y6nXI4HCosLFR0dHQTbhUAAF7sRKn0ZCeppFD63ZdScj+rEzVKQ9+/PXbMTGFhoWJjY13Xs7KyNHDgQFeRkaSMjAxt2rRJhw4dco0ZPHhwjefJyMhQVlaWZ0IDAODLQsKk8zPM9R/897uaPFJmNm/erBdeeEF33HGH67b8/HwlJCTUGFd5PT8/v94xlffXpaSkRE6ns8YCAEBAuuAaqUV7KSrR6iRu06gyM2nSJNlstnqXyl1ElXJzczV06FDdeOONuu2225o0fF2mTp0qh8PhWpKSkjzycwEA8Dpdr5V+v0YacLfVSdymUV9nMHHiRI0ZM6beMR06dHCt5+Xl6fLLL9eAAQNqHNgrSYmJiSooKKhxW+X1xMTEesdU3l+XyZMn67777nNddzqdFBoAQGAK8v+zsDSqzMTFxSkuLq5BY3Nzc3X55ZcrLS1Ns2bNUtApf5jp6el6+OGHVVZWptDQUEnS/Pnz1aVLF7Vo0cI1ZsGCBZowYYLrcfPnz1d6enq9P9tut8tutzdiywAA8HMnSqXd/5PaX2p1kibnlrqWm5uryy67TMnJyXrqqae0b98+5efn1zjW5ZZbblFYWJjGjh2rnJwczZ07V88991yNGZV77rlHn3/+uZ5++mlt3LhRjz32mFauXKnx48e7IzYAAP6p7Lj09PnS7GHS4Z1Wp2lybikz8+fP1+bNm7VgwQK1bdtWrVu3di2VHA6HvvzyS23btk1paWmaOHGipkyZottvv901ZsCAAZozZ45eeeUV9ejRQ++//74++ugjde/e3R2xAQDwT6HhUlxXc/2nL63N4gYeO8+MlTjPDAAg4H3ztLTgT1LnDGnUe1anaRCvO88MAACwUOeT55vZ9rVUVmxtliZGmQEAIBAkXChFnyedKJa2f2t1miZFmQEAIBDYbFLnn5vrfnbcDGUGAIBAUbmr6ccvJD86ZLZR55kBAAA+rMMg6YpHqkqNn6DMAAAQKMKaSwPvtzpFk2M3EwAA8GmUGQAAAkl5mbT2Pek/vzfX/QBlBgCAQGILlj57UFr1hpS7yuo0TYIyAwBAIAkKMg8ElqSti6zN0kQoMwAABJoOl5mXWxdbmaLJUGYAAAg0lWVm9wqppMjSKE2BMgMAQKBp0V5qkSJVnJC2f2d1mnNGmQEAIBBVzs5s+9rSGE2BMgMAQCBqf6l56dxtbY4mwBmAAQAIROcPlf6wWYqMszrJOaPMAAAQiOyR5uIH2M0EAECg8/Fv0KbMAAAQqApypDevNxcfxm4mAAAClT3KPHFeUIhUcsRndzsxMwMAQKCKSZYcSeb5Znb/z+o0Z40yAwBAIGt3iXm5Y6m1Oc4BZQYAgEDWnjIDAAB8WVJ/8zI3WyovszbLWaLMAAAQyFp2kpq1kE4cl/LXWp3mrPBpJgAAAllQkNRpsFR8SPLR081QZgAACHQ3vGp1gnPCbiYAAODTKDMAAMBUVCCVHrU6RaNRZgAAgDTnJunp86Wf5ludpNEoMwAAQHKcZ17u8r0zAVNmAACAlNTPvNy13NocZ4EyAwAApKS+5uWe76Wy49ZmaSTKDAAAkGLaSREtpYoyqSDH6jSNQpkBAACSzSa16WWu562yNksjUWYAAIDpvJNlJte3ygxnAAYAAKaOV0jHDpqXPoQyAwAATMn9zcXHsJsJAAD4NGZmAABAlbJiac9aKcQutelpdZoGYWYGAABUWfay9PoQ6btnrU7SYJQZAABQpY3vfaKJMgMAAKq0udi8PLxDOnrA2iwNRJkBAABVmsVIsR3N9T2rLY3SUJQZAABQU+se5mX+OmtzNBBlBgAA1NT6IvOSMgMAAHxSYqp5uWettTkaiPPMAACAmtr0koY+UbW7yctRZgAAQE0RsVL/cVanaDB2MwEAAJ/GzAwAADidc4+0bYkUHCZ1H2F1mnoxMwMAAE63a7n04R3Sd89ZneSMKDMAAOB0lR/P3rtBKi+zNssZUGYAAMDpYtpL9mipvFTa/6PVaepFmQEAAKcLCvKZ881QZgAAQO0qy4yXnwmYMgMAAGoX38283LvB2hxnQJkBAAC1S7jQvPTyMsN5ZgAAQO3iu0mj3q+aofFSlBkAAFC7sAip88+tTnFG7GYCAAA+ze1lpqSkRD179pTNZtOaNWtq3Ld27Vr97Gc/U3h4uJKSkjRt2rTTHj9v3jx17dpV4eHhSk1N1aeffuruyAAAoFLBBmnhX6UVr1qdpE5uLzMPPPCA2rRpc9rtTqdTQ4YMUbt27ZSdna0nn3xSjz32mF555RXXmKVLl+rmm2/W2LFjtXr1ag0fPlzDhw/X+vXr3R0bAABI0r4fpK+nSd+/a3WSOrm1zHz22Wf68ssv9dRTT51239tvv63S0lK9/vrruvDCC3XTTTfp97//vZ555hnXmOeee05Dhw7V/fffrwsuuEB//vOf1atXL7344ovujA0AACrFV36i6QeposLaLHVwW5kpKCjQbbfdpn/961+KiIg47f6srCwNHDhQYWFhrtsyMjK0adMmHTp0yDVm8ODBNR6XkZGhrKysen92SUmJnE5njQUAAJyFlh3Nb84uPSIV7rQ6Ta3cUmYMw9CYMWN05513qnfv3rWOyc/PV0JCQo3bKq/n5+fXO6by/rpMnTpVDofDtSQlJZ3tpgAAENiCQ6VWXcz1Au8830yjysykSZNks9nqXTZu3KgXXnhBRUVFmjx5srty12vy5MkqLCx0Lbt27bIkBwAAfiGh8kzAOdbmqEOjzjMzceJEjRkzpt4xHTp00MKFC5WVlSW73V7jvt69e2vUqFF64403lJiYqIKCghr3V15PTEx0XdY2pvL+utjt9tN+NgAAOEuVJ83z0pmZRpWZuLg4xcXFnXHc888/r7/85S+u63l5ecrIyNDcuXPVr18/SVJ6eroefvhhlZWVKTQ0VJI0f/58denSRS1atHCNWbBggSZMmOB6rvnz5ys9Pb0xsQEAwLmoLDMHNlubow5uOQNwcnJyjeuRkZGSpI4dO6pt27aSpFtuuUWPP/64xo4dqwcffFDr16/Xc889p+nTp7sed88992jQoEF6+umnNWzYML377rtauXJljY9vAwAAN2s3QLp7ldSivdVJamXZGYAdDoe+/PJLbdu2TWlpaZo4caKmTJmi22+/3TVmwIABmjNnjl555RX16NFD77//vj766CN1797dqtgAAAQee6T5qaagYKuT1MpmGIZhdQh3czqdcjgcKiwsVHR0tNVxAABAAzT0/ZvvZgIAAGe28b/S+2OlNXOsTnIaygwAADizvRuk9e9L2762OslpKDMAAODMKk+ct2+TtTlqQZkBAABnFneyzOz/UfKyw20pMwAA4MxiO0hBIeZ3NDlzrU5TA2UGAACcWXCoFNvRXPeyXU2UGQAA0DBx55uXlBkAAOCTWnUxdzUVH7Q6SQ1u+ToDAADghy6dIA16UAoJszpJDZQZAADQMPYoqxPUit1MAADAp1FmAABAw33xsPTPK6XcbKuTuFBmAABAw+Wvk3JXSns3Wp3EhTIDAAAarmUn8/LAZmtzVEOZAQAADUeZAQAAPs1VZrZYm6MaygwAAGi4lie/0uDgFqmiwtosJ1FmAABAw8W0M88CfOK413zhJGUGAAA0XHCIuaupRXvp2AGr00jiDMAAAKCxxi2VgoKtTuHCzAwAAGgcLyoyEmUGAAD4OMoMAABonP2bpdcypH9eYXUSSRwzAwAAGsseKe1aJtmCpBMlUojd0jjMzAAAgMaJTJDCIiWjQjq03eo0lBkAANBINlvVyfO84GsNKDMAAKDxYjuYlwe3WZtDlBkAAHA2WqSYl4coMwAAwBe5Zma2WptDlBkAAHA2WnaSHMnmwcAW46PZAACg8dqlS/euszqFJGZmAACAj6PMAACAc2MYlv54ygwAADg7Xz0uPdVFWvmapTEoMwAA4OyUl0pH8i0/1wxlBgAAnJ0W7SVJB3dt1J7CYstiUGYAAMBZWbI/UpK0b+dGXfLEQs1dsdOSHJQZAADQaHsKi/XoN+ZsTLJtrwyjQg99sN6SGRrOMwMAABpt2/6j2m201KsnrtJ2I1EhqlCZEaTt+4+ptaOZR7NQZgAAQKOltGquCluI/nLiVtdtwTab2reK8HgWdjMBAIBGa+1opqkjUhVss0kyi8zfRnT3+KyMxMwMAAA4SyP7JGvg+XHavv+Y2reKsKTISJQZAABwDlo7mllWYiqxmwkAAPg0ygwAAPBplBkAAODTKDMAAMCnUWYAAIBPo8wAAACfRpkBAAA+jTIDAAB8GmUGAAD4NMoMAADwaZQZAADg0wLiu5kMw5AkOZ1Oi5MAAICGqnzfrnwfr0tAlJmioiJJUlJSksVJAABAYxUVFcnhcNR5v804U93xAxUVFcrLy1NUVJRsNluTPa/T6VRSUpJ27dql6OjoJnteb+Lv28j2+T5/30a2z/f5+za6c/sMw1BRUZHatGmjoKC6j4wJiJmZoKAgtW3b1m3PHx0d7Zd/Qavz921k+3yfv28j2+f7/H0b3bV99c3IVOIAYAAA4NMoMwAAwKdRZs6B3W7Xo48+KrvdbnUUt/H3bWT7fJ+/byPb5/v8fRu9YfsC4gBgAADgv5iZAQAAPo0yAwAAfBplBgAA+DTKDAAA8GmUmUbYvn27xo4dq5SUFDVr1kwdO3bUo48+qtLS0nofd/z4cWVmZqply5aKjIzUDTfcoIKCAg+lbpy//vWvGjBggCIiIhQTE9Ogx4wZM0Y2m63GMnToUPcGPQdns42GYWjKlClq3bq1mjVrpsGDB+unn35yb9CzdPDgQY0aNUrR0dGKiYnR2LFjdeTIkXofc9lll532Gt55550eSnxmL730ktq3b6/w8HD169dP//vf/+odP2/ePHXt2lXh4eFKTU3Vp59+6qGkZ6cx2zd79uzTXqvw8HAPpm2cr7/+Wtdee63atGkjm82mjz766IyPWbx4sXr16iW73a5OnTpp9uzZbs95thq7fYsXLz7t9bPZbMrPz/dM4EaaOnWq+vTpo6ioKMXHx2v48OHatGnTGR/n6d9BykwjbNy4URUVFZo5c6ZycnI0ffp0zZgxQw899FC9j7v33nv18ccfa968eVqyZIny8vI0YsQID6VunNLSUt14440aN25cox43dOhQ7dmzx7W88847bkp47s5mG6dNm6bnn39eM2bM0PLly9W8eXNlZGTo+PHjbkx6dkaNGqWcnBzNnz9fn3zyib7++mvdfvvtZ3zcbbfdVuM1nDZtmgfSntncuXN133336dFHH9WqVavUo0cPZWRkaO/evbWOX7p0qW6++WaNHTtWq1ev1vDhwzV8+HCtX7/ew8kbprHbJ5lnWq3+Wu3YscODiRvn6NGj6tGjh1566aUGjd+2bZuGDRumyy+/XGvWrNGECRP0f//3f/riiy/cnPTsNHb7Km3atKnGaxgfH++mhOdmyZIlyszM1LJlyzR//nyVlZVpyJAhOnr0aJ2PseR30MA5mTZtmpGSklLn/YcPHzZCQ0ONefPmuW774YcfDElGVlaWJyKelVmzZhkOh6NBY0ePHm1cf/31bs3jDg3dxoqKCiMxMdF48sknXbcdPnzYsNvtxjvvvOPGhI23YcMGQ5KxYsUK122fffaZYbPZjNzc3DofN2jQIOOee+7xQMLG69u3r5GZmem6Xl5ebrRp08aYOnVqreN/9atfGcOGDatxW79+/Yw77rjDrTnPVmO3rzG/m95GkvHhhx/WO+aBBx4wLrzwwhq3jRw50sjIyHBjsqbRkO1btGiRIck4dOiQRzI1tb179xqSjCVLltQ5xorfQWZmzlFhYaFiY2PrvD87O1tlZWUaPHiw67auXbsqOTlZWVlZnojoEYsXL1Z8fLy6dOmicePG6cCBA1ZHajLbtm1Tfn5+jdfQ4XCoX79+XvcaZmVlKSYmRr1793bdNnjwYAUFBWn58uX1Pvbtt99Wq1at1L17d02ePFnHjh1zd9wzKi0tVXZ2do0/+6CgIA0ePLjOP/usrKwa4yUpIyPD614r6ey2T5KOHDmidu3aKSkpSddff71ycnI8EdcjfOn1Oxc9e/ZU69at9fOf/1zfffed1XEarLCwUJLqfd+z4jUMiC+adJfNmzfrhRde0FNPPVXnmPz8fIWFhZ12bEZCQoLX7iNtrKFDh2rEiBFKSUnRli1b9NBDD+mqq65SVlaWgoODrY53zipfp4SEhBq3e+NrmJ+ff9p0dUhIiGJjY+vNesstt6hdu3Zq06aN1q5dqwcffFCbNm3SBx984O7I9dq/f7/Ky8tr/bPfuHFjrY/Jz8/3iddKOrvt69Kli15//XVddNFFKiws1FNPPaUBAwYoJyfHrV+o6yl1vX5Op1PFxcVq1qyZRcmaRuvWrTVjxgz17t1bJSUlevXVV3XZZZdp+fLl6tWrl9Xx6lVRUaEJEybokksuUffu3escZ8XvIDMzkiZNmlTrAVnVl1P/YcnNzdXQoUN144036rbbbrMoecOczfY1xk033aTrrrtOqampGj58uD755BOtWLFCixcvbrqNOAN3b6PV3L19t99+uzIyMpSamqpRo0bpzTff1IcffqgtW7Y04VagKaSnp+s3v/mNevbsqUGDBumDDz5QXFycZs6caXU0NECXLl10xx13KC0tTQMGDNDrr7+uAQMGaPr06VZHO6PMzEytX79e7777rtVRTsPMjKSJEydqzJgx9Y7p0KGDaz0vL0+XX365BgwYoFdeeaXexyUmJqq0tFSHDx+uMTtTUFCgxMTEc4ndYI3dvnPVoUMHtWrVSps3b9aVV17ZZM9bH3duY+XrVFBQoNatW7tuLygoUM+ePc/qORuroduXmJh42oGjJ06c0MGDBxv1961fv36SzNnHjh07NjpvU2nVqpWCg4NP+/Rffb8/iYmJjRpvpbPZvlOFhobq4osv1ubNm90R0ePqev2io6N9flamLn379tW3335rdYx6jR8/3vWBgjPNAFrxO0iZkRQXF6e4uLgGjc3NzdXll1+utLQ0zZo1S0FB9U9upaWlKTQ0VAsWLNANN9wgyTyKfefOnUpPTz/n7A3RmO1rCrt379aBAwdqvPG7mzu3MSUlRYmJiVqwYIGrvDidTi1fvrzRn/o6Ww3dvvT0dB0+fFjZ2dlKS0uTJC1cuFAVFRWugtIQa9askSSPvoa1CQsLU1pamhYsWKDhw4dLMqe6FyxYoPHjx9f6mPT0dC1YsEATJkxw3TZ//nyP/b41xtls36nKy8u1bt06XX311W5M6jnp6emnfYzXW1+/prJmzRrLf9fqYhiG7r77bn344YdavHixUlJSzvgYS34H3XZosR/avXu30alTJ+PKK680du/ebezZs8e1VB/TpUsXY/ny5a7b7rzzTiM5OdlYuHChsXLlSiM9Pd1IT0+3YhPOaMeOHcbq1auNxx9/3IiMjDRWr15trF692igqKnKN6dKli/HBBx8YhmEYRUVFxh/+8AcjKyvL2LZtm/HVV18ZvXr1Mjp37mwcP37cqs2oV2O30TAM44knnjBiYmKMf//738batWuN66+/3khJSTGKi4ut2IR6DR061Lj44ouN5cuXG99++63RuXNn4+abb3bdf+rf0c2bNxt/+tOfjJUrVxrbtm0z/v3vfxsdOnQwBg4caNUm1PDuu+8adrvdmD17trFhwwbj9ttvN2JiYoz8/HzDMAzj1ltvNSZNmuQa/9133xkhISHGU089Zfzwww/Go48+aoSGhhrr1q2zahPq1djte/zxx40vvvjC2LJli5GdnW3cdNNNRnh4uJGTk2PVJtSrqKjI9TsmyXjmmWeM1atXGzt27DAMwzAmTZpk3Hrrra7xW7duNSIiIoz777/f+OGHH4yXXnrJCA4ONj7//HOrNqFejd2+6dOnGx999JHx008/GevWrTPuueceIygoyPjqq6+s2oR6jRs3znA4HMbixYtrvOcdO3bMNcYbfgcpM40wa9YsQ1KtS6Vt27YZkoxFixa5bisuLjbuuusuo0WLFkZERITxi1/8okYB8iajR4+udfuqb48kY9asWYZhGMaxY8eMIUOGGHFxcUZoaKjRrl0747bbbnP9Q+yNGruNhmF+PPuPf/yjkZCQYNjtduPKK680Nm3a5PnwDXDgwAHj5ptvNiIjI43o6Gjjt7/9bY2idurf0Z07dxoDBw40YmNjDbvdbnTq1Mm4//77jcLCQou24HQvvPCCkZycbISFhRl9+/Y1li1b5rpv0KBBxujRo2uMf++994zzzz/fCAsLMy688ELjv//9r4cTN05jtm/ChAmusQkJCcbVV19trFq1yoLUDVP5UeRTl8ptGj16tDFo0KDTHtOzZ08jLCzM6NChQ43fRW/T2O37+9//bnTs2NEIDw83YmNjjcsuu8xYuHChNeEboK73vOqviTf8DtpOhgUAAPBJfJoJAAD4NMoMAADwaZQZAADg0ygzAADAp1FmAACAT6PMAAAAn0aZAQAAPo0yAwAAfBplBgAA+DTKDAAA8GmUGQAA4NMoMwAAwKf9f/vMiT1MPal4AAAAAElFTkSuQmCC",
      "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",
    "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.\n",
    "\n",
    "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}