07-Interpolation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Interpolation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 07.01 Interpolation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Given input data such as:\n",
    "$$\n",
    "(t_1, y_1), (t_2, y_2), \\cdots, (t_m, y_m) \\quad t_1 < t_2 < \\cdots < t_m\n",
    "$$\n",
    "\n",
    "determine a function $f : \\mathbb{R} \\rightarrow \\mathbb{R}$ such that:\n",
    "$$\n",
    "f(t_i) = y_i \\quad i = 1, \\cdots, m\n",
    "$$\n",
    "\n",
    "$f$ is referred to as **interpolating function**\n",
    "\n",
    "Additional constraints might be imposed:\n",
    "* smoothness\n",
    "* monotonicity\n",
    "* convexity\n",
    "\n",
    "Use cases:\n",
    "* graaphics or plotting\n",
    "* differentiating or integrating \n",
    "* replacing complex function with simpler one\n",
    "\n",
    "#### Interpolation vs. Approximation\n",
    "Interpolation fits given input data points **exactly** whereas approximation does not.\n",
    "\n",
    "Example of approximation: Smooth noisy data using least squares.\n",
    "\n",
    "#### Basis Functions\n",
    "Family of functions for interpolating data points spanned by a set of basis functions $\\phi_1(t), \\cdots, \\phi_n(t)$.\n",
    "\n",
    "Interpolating function $f$ is linear combination of basis functions:\n",
    "$$\n",
    "f(t_i) = \\sum_{j=1}^{n} x_j \\phi_j (t_i) = y_i\n",
    "$$\n",
    "where\n",
    "* $x$ is an n-vector of parameters\n",
    "* $\\phi_j (t_i)$ is the value of *jth* basis function at $t_i$ forms an $m \\times n$ matrix $A$\n",
    "* $y$ is the product of the linear system $Ax$\n",
    "\n",
    "#### Existence, Uniqueness, and Conditioning\n",
    "Existence depends on the number of data points $m$ and basis functions $n$.\n",
    "* If $m > n$, then interpolant usually doesn't exist.\n",
    "* If $m < n$, then interpolant is not unique.\n",
    "* If $m = n$, then basis matrix $A$ is nonsingular and data can be fit exactly.\n",
    "\n",
    "Sensitivity of $x$ depends on $\\text{cond}(A)$ which depends on basis functions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 07.02 Polynomial Interpolation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Unique polynomial of degree $n-1$ passsing through $n$ distinct data points.\n",
    "\n",
    "#### Monomial Basis\n",
    "Basis functions given by sequential powers:\n",
    "$$\n",
    "\\phi_j (t) = t^{j-1} \\quad j=1, \\cdots, n\n",
    "$$\n",
    "\n",
    "Polynomial evaluated at $t_i$ with this basis:\n",
    "$$\n",
    "p_{n-1} (t_i) = x_1 + x_2 t_i + x_3 t_i^2 + \\cdots + x_n t_i^{n-1}\n",
    "$$\n",
    "\n",
    "The $n \\times n$ matrix is known as **Vandermonde matrix**:\n",
    "$$\n",
    "Ax = \\begin{bmatrix}\n",
    "1 & t_1 & t_1^2 & \\cdots & t_1^{n-1} \\\\\n",
    "\\vdots & \\vdots & \\vdots & \\ddots & \\vdots \\\\\n",
    "1 & t_n & t_n^2 & \\cdots & t_n^{n-1} \\\\\n",
    "\\end{bmatrix}\n",
    "\\begin{bmatrix}\n",
    "x_1 \\\\\n",
    "x_2 \\\\\n",
    "\\vdots \\\\\n",
    "x_n\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "For monomial basis, matrix $\\text{cond}(A)$ increases with degree of polynomial.\n",
    "* Ill-conditioning does not prevent us from fitting points.\n",
    "* Residuals at fitted data points is zero.\n",
    "* Coefficients will be poorly determined.\n",
    "* Relatively high cost to solve linear system $O(n^3)$.\n",
    "* Scaling can help reduce the growth in condition number.\n",
    "* Cost of evaluation using Horner's rule (see below) is $2n$.\n",
    "\n",
    "#### Horner's Rule\n",
    "Rather than explicitly compute powers, performed nested multiplications of $n$ additions and $n$ multiplications.\n",
    "$$\n",
    "p_{n-1} (t_i) = x_1 + t_i (x_2 + t_i (x_3 + t_i (\\cdots + (x_{n-1} + tx_n) \\cdots)))\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Determine the polynomial of degree 2 used to interpolate 3 ($m=3$) data points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x:  [-1.  5. -4.]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "# Input data $(t_i, y_i)$.\n",
    "m = 3\n",
    "t = np.array([-2,0,1], dtype='d')\n",
    "y = np.array([-27,-1,0], dtype='d')\n",
    "\n",
    "# Solve Ax = y for x.\n",
    "A = np.column_stack((np.ones(m), t, t*t))\n",
    "x = np.linalg.solve(A, y)\n",
    "print(\"x: \", x)\n",
    "\n",
    "# Evaluate the polynomial at the data points.\n",
    "vft = np.vectorize(lambda t: x[0] + t*(x[1] + x[2]*t))  # Horner's rule.\n",
    "yhat = vft(t)\n",
    "np.testing.assert_almost_equal(y - yhat, np.zeros(m))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the observations against the fit using polynomial interpolation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAEWCAYAAACNJFuYAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAw4ElEQVR4nO3defxWc/7/8cervZRKaKiUEKVJKyYyRUNMKmnslJB9GWOyjWUG30F+zGSZDJK9TKbFMlH0sYcktCChVKKS6tP6+Xx6/f54n08un30/13V9nvfb7dyu5ZzrXK/3Odd1va73+33OeZu7IyIikqhG3AGIiEjyUXIQEZF8lBxERCQfJQcREclHyUFERPJRchARkXyUHKqImd1sZk/GHUdpmdn/zGxo3HFUBjPrZWafx/C+SfdZMLPTzeyVEi6bdPGXhZm1MTM3s1plfP11ZvZwRceVLJQcKoiZDTOzT81sk5mtNLN/mVmTuOMqjYK+9O5+rLs/FldMFSHhRyAzYfrY3d909/0TlvvGzPqWYr03m1lWtL6fzOwdM/tN5ZSicrn7U+5+dHnXY2a9zWxZCZct149zVSqoXO7+f+5+blwxVTYlhwpgZn8C7gD+DDQGDgVaA9PNrE4VxpH0X7KYNXH3htF0UAWtc4K7NwR2A94C/mtmVkHrliLo8165lBzKycx2Bv4KXOru09w9y92/AU4C2gBnJCxez8wmmNkGM5tjZgclrOdqM1sezfvczI6Knq9hZteY2WIzW2Nmz5rZLtG83H9e55jZUuC1qBnokjwxfmxmg6P7/zSzb81svZl9aGa9ouf7AdcBJ+f+s46ezzCzcxNi+YuZLTGzH8zscTNrnCeWoWa21MxWm9n1CTEcbGazo/f93szuLmR7LjSz/gmPa5nZKjPramb1zOzJaDv8ZGYfmFnzMuy23HXv+DdoZk8AewHPR+UfWZp1uXsW8BjwK6CZme1pZlPN7Ecz+9LMziskhhfN7NI8z31iZidE993MLjCzRVGZ789NPiXcH2dH+3tttJ4e0fp/MrP7Et5zmJm9lfC4wM9JaUWfn1vM7O3os/2Kme0azX4juv0p2ua/iV4zPPocrDWzl82sdcL63MwuNrNFwKKE5y4zs6+iz90oM6tR3DYqINazo/fdEK3r/Oj5nYD/AXvazzXPPS1PTdvMBpjZ/GjbZphZ+4R535jZVdG2X2fhd6BeWbZplXF3TeWYgH5ANlCrgHmPAc9E928GsoAhQG3gKuDr6P7+wLfAntGybYB9ovuXA7OAlkBd4MGEdbYBHHgc2AmoD5wFvJ0QQwfgJ6Bu9PgMoBlQC/gTsBKolxDjk3nKkAGcG90fDnwJtAUaAv8FnsgTy0NRHAcBW4H20fx3gTOj+w2BQwvZnjcCTyU8/j2wMLp/PvA80ACoCXQDdi7BPsqNrVae53sDyxIefwP0LcW+37G9on0zClgaPX4DeACoB3QGVgFHFvC6k4D3EtZ5ELAGqBM9duAFoAkhea0C+pVif4yJYjga2AJMBnYHWgA/AL+Nlh8GvJUQR6k+J0Vs0wxgMdAu+lxkALcXtl+AgVGZ2kfv/RfgnYT5DkwHdgHqJzw3M3puL+ALSveZrZXwWdsHMOC3wCaga0HlKmA/tgM2Ar8jfKdHRu+bux+/Ad4H9oziXAhcEPfvV5Gf77gDSPUp+hKtLGTe7cD0hA/SrIR5NYDvgF7AvtEXtS9QO886FgJHJTzeg5BkaiV8uNsmzG8UfUhbR49vA8YWEf9a4KCEGItKDq8CFyXM27+AWFomzH8fOCW6/wahhrVrMdtzX2AD0CB6/BRwY3R/OPAO0KmU+yg3tp8SpqvyfuEpW3LYFq3vB+A1QsJqBeQAjRKW/TswLu92JvxwrwX2ix7fBTyQ8DoHDk94/CxwTSn2R4uE+WuAkxMePwdcEd0fRkJyKO3nJGG5vNs0A/hLwuOLgGl59kticvgfcE6e78kmfv48O1GSzbON+uV5j1dLsY3y/bGLlp0MXF5QuQrYjzcAz+aJeznQO+GzdUbC/DuBMaX5HFf1pGal8lsN7GoFt3/uEc3P9W3uHXffDiwj1Ba+BK4gfNh+MLPxZrZntGhrYFJUVf2JkCxygOaFrHcD8CJwSvTUqYQfWACiqu3CqGr7E6GPZFdKZk9gScLjJYQvWWIsKxPubyL8WwM4h/Dv6rOoOag/BYi2xULgeDNrAAwAno5mPwG8DIw3sxVmdqeZ1S5h7BASU5NouqsUryvKs9H6dnf3I939Q8J2+jHaF7mWEP6t/4K7bwEmAGdETSGnEsqZqLBtWpL98X3C/c0FPG5IAcr5OcmrsPgL0hr4Z8Ln/UfCP/nEbfdtAa9LfG4JYdtAybYRAGZ2rJnNipoCfwKOo4zfjej7/W2euEuzHWKn5FB+7xKaTwYnPmlmDYFjCf9ccrVKmF+D0FS0AsDdn3b3wwlfDid0cEP4gB2b8KPWxN3rufvyhPV6npieAU6N2nDrEarcRO3GIwlNGU3dvQmwjvDlK2g9ea2I4su1F6FJ7fuCF08I0H2Ru59KaNK4A5gYteUW5BnCj+RAYEGUMPDQn/NXd+8A9AT6E5rRKkpx5S+pFcAuZtYo4bm9CP8kC/IYcDpwFLDJ3d8txfuUaX8UpQSfk4pS0Pb+Fjg/z+e9vru/U8zrWiXc34voe0UJt5GZ1SXUpO4CmkdlfokyfjeifqFWFL7Pk56SQzm5+zpCc8m9ZtbPzGqbWRtC9X8Zv/wX2M3MBke1jCsISWWWme1vZkdGH9AthH9026PXjAFuy+2UM7PdzGxgMWG9RPig/o1wNE3uuhoRvhirgFpmdiOwc8Lrvgfa5HbmFeAZ4I9mtneU/P4vWn92MfFgZmeY2W5RLD9FT28vZPHxhDbyC/m51oCZ9TGzX5tZTWA9oXmgsHWUxfeEtunEuL8xs2GlWYm7f0to/vq7hU70ToSaU4HnBkTJYDvw/8hfayhKmfdHMYr7nFSUVYRyJ27zMcC1ZnYggJk1NrM/lGBdfzazpmbWitBPNyF6vqTbqA6h32gVkG1mxxI+g7m+Jxxo0LiQ938W+L2ZHRXVZv9E+H6/U8jySU/JoQK4+52EI33uIvxovUf4B3SUu29NWHQKcDKh/fZMYLCHo1zqEvonVhOqnrsD10av+ScwFXjFzDYQOqcPKSaerYSOt74k/LgSmmSmETrslhASUWJ1/D/R7Rozm1PAqscSfrzeIHSmbwEuLWC5gvQD5ptZZlSmU9x9cyHxf0eokfXk5y85hCOBJhK28ULg9SgezGyMmY0pYSyF+Tvwl6hJ4yoLhyE3I2zz0jqV0Ka9ApgE3OTuM4pY/nHg1xSSQApRnv1RlOI+JxXC3TcR+sTejrb5oe4+iVCzHG9m64F5hBp4caYAHwJzCc2qj0TPl2gbRU2AlxF+5NcCpxG+d7nzPyMkmq+iWPfM8/rPCf2P9xK+x8cDx7v7thLEnpQs6hwRkTzM7HDg4qg5rLLf6yxgRNS0KKVgZk7o0P8y7ljSiU4iESmEu79FOLGtUkUd7xcRDn0VSQpqVhKJkZkdQ2jn/p5fNgGKxErNSiIiko9qDiIikk9a9Dnsuuuu3qZNmzK/fuPGjey0U2GH3KeOdCkHqCzJKF3KASpLrg8//HC1u+9W0Ly0SA5t2rRh9uzZZX59RkYGvXv3rriAYpIu5QCVJRmlSzlAZcllZksKm6dmJRERyUfJQURE8lFyEBGRfNKiz6EgWVlZLFu2jC1bthS7bOPGjVm4cGEVRFW5krUc9erVo2XLltSuXZoLqIpInNI2OSxbtoxGjRrRpk0brJhRGzds2ECjRo2KXCYVJGM53J01a9awbNky9t5777jDEZESSttmpS1bttCsWbNiE4NULjOjWbNmJarBiUjySNvkACgxJAntB5HUk7bNSiIiBdq2jdo//QSLFsG6dfDTT+F2yxbYvDnc5k45OeD+8wRQqxbUqQN16/5826AB7Lzzz1OjRtCkCTRtCjVS8z+4kkMlGz16NP/6179YuXIlV199Nddccw2TJ0+mXbt2dOjQIe7wRNJHdjasWAFLl8KSJeF22TL4/vtfTuvXc1hVxVSzJuy6K+y2W5iaN4eWLaFVq59v99oLdt8dkqyGnbTJwcz6EQaFqQk87O63xxxSmTzwwAPMmDGDli1b7nhu8uTJ9O/fX8lBpLTcw4/+Z5/BF1+Ef/9ffBGmpUvDP/1ETZuGH+TmzaFLl3C7224sWr2a/bp3D//umzQJ//YbNIB69cJUv36oEdSsGX60E3+4c3Jg27Ywbd0apk2bYMMGWL/+52ntWli1Cn74IdyuWgUffACTJoXXJNppJ9h3X9hnn3C7777QoQMceGCILwZJmRyiYSDvB35HGGrzAzOb6u4L4o2sdC644AK++uorjj32WIYPH87ixYs57bTTmDp1Kq+//jq33norzz33HPvss0/coYoknw0b4KOP4JNP4NNPwzRvXng+V6NG0K4dHHoonHYatG4d/om3bh3+lRdyzaHlGRnsV9bLZ9SsGZJH/fple707rF4dajXffhuS2uLF8OWXsGABvPBCSDy59twzJImOHeGgg6BbNzjgAJZurMXwqTCoLtz6BIwdAHsVNohpGSRlcgAOBr50968AzGw80WDzZVrbFVfA3LmFzq6fkxN2eGl07gz/+EeRi4wZM4Zp06Yxc+ZMXnjhBQB69uzJgAED6N+/P0OGDCnde4qkq23bQiL44IMwzZ4NCxf+3M7ftCn8+tdw1lnhtkMH2G+/UBNIsuaYYpn93MzUpUv++Tk5IWksWADz54eEOH8+jBkT+kQA6tcnc8+DGLpHV5oeugeLW/Vm+FSYcWbFhZmsyaEFvxyzdhl5xk02sxHACIDmzZuTkZHxixU0btyYDdE/jLrbtlEjb3UzkTvZRc0vwPZt29ia+A+m0FU7mZmZbNmyhW3btrFhwwaysrLYvHnzjvgqSk5OToWvs6Js2bIl3z4qSmZmZqmWT2bpUpaKLEetzEx2nj+fxp9+SuNPP6XRZ59RM/q3vK1pU9YfcAAbhg4ls107Nuy7L9t23TV/085nn4WpDFJinzRoAD16hAkgJ4cGy5bR6IsvaPjFFzSe/wUnfziO5bUP5coeh2NARRYpWZNDsdz938C/Abp37+55r0q4cOHCn08Ie6Do0RfLevJYnRIsY2Y0bNiQevXqUadOHRo1akTt2rWpX79+hZ+wlownweWqV68eXQr6l1QIXTUz+ZSrHJs2wVtvwYwZ8OqroZbgHo786dIFLr4YDjsMDjmEOi1asKsZu1Zo9L+UDvuk7xPw1ZrtjGw6jbtX9mafpnBJ74pbf7Imh+VAq4THLaPn0kKjRo2S9h++SIVwD30EL74IL78M774bmo5q14bf/AZuugmOOAIOPrjQfgEp2tgBMHxqDbLrNmCfGuFxRUrW5PABsJ+Z7U1ICqcAp8UbUsU55ZRTOO+88xg9ejQTJ05Uh7Skh40bQ83gpZfCtGxZeL5zZ7jsMjjqKOjVS8mgguzVOPQxZGRUbI0hV1ImB3fPNrNLgJcJh7KOdff5MYdVJt988w0Aw4YNY9iwYQAcdthhLFiQUgdeiRRs7Vp4/nn4739DDWHLlnAE0dFHw1//CsceC3vsEXeUUgZJmRwA3P0l4KW44xCp7pau4xeHTI47fA0tX5kIzz0HM2eGk89atoTzzoNBg+Dww8OZw5LSkjY5iEhyGD4VVn63gXbfTue610bxq89fgZzscCjpVVfB4MHQvXvqHVIqRVJyEJGC5eTA9Olcfs84+s6bSv2szSxv0opHel/J+XedGk7IUkJIW0oOIvJLn38O48bB44/DihUc0bAZEw85m/q/O4CRO19M22Y1OL9z3EFKZVNyEJFw5u2zz8KDD4bDTmvWDJ3Jo0ezrld/Hnu5LoPqZtB2a40KP2RSkpOSg0h1tnBhSAiPPRYuXb3//nDnnXDGGTuOMtqLyj1kUpJTal5oPIV98803dOzYMe4wmDt3Li+99PPBYFOnTuX221PywrdSWtnZ4UijPn3CNYoeeAD69Qu//gsXwp//rMNPRTWHdJCdnU2tWqXblXPnzmX27Nkcd9xxAAwYMIABA9RekNbWroWHH4b77gtXAm3TBm6/Hc4+O4wnIJJANYdKdvfdd9OxY0c6duzIP6KruGZnZ3P66afTvn17hgwZwqZNmwC45ppr6NChA506deKqq64CYNWqVZx44on06NGDHj168PbbbwNw8803c+aZZ3LYYYdx5plncuihh7Jw4cId79u7d29mz57N+++/z29+8xu6dOlCz549+fzzz9m2bRs33ngjEyZMoHPnzkyYMIFx48ZxySWXAKF2c+SRR9KpUyeOOuooli5dCoQT+S677DJ69uxJ27ZtmThxIgDfffcdRxxxBJ07d6Zjx468+eabVbJtpYS++AIuuiicizByJLRtG8YU+PJLuPpqJQYpmLun/NStWzfPa8GCBfmeK8z69et9yU/uRz3uvvc/w+2Sn0r88kLNnj3bO3bs6JmZmb5hwwbv0KGDz5kzxwF/66233N397LPP9lGjRvnq1au9Xbt2vn37dnd3X7t2rbu7n3rqqf7mm2+6u/uSJUv8gAMOcHf3m266ybt27eqbNm1yd/e7777br776and3X7Fihbdr187d3detW+dZWVnu7j59+nQfPHiwu7s/+uijfvHFF++INfFx//79fdy4ce7u/sgjj/jAgQPd3X3o0KE+ZMgQz8nJ8fnz5/s+++zj7u533XWX33rrre7unp2d7evXr8+3LUqzP9zdZ86cWarlk1lsZXnvPffBg93N3OvUcT/7bPePPirz6rRPklN5ygLM9kJ+V9WsFBk+FRavhe0ebivi2uhvvfUWJ5xwAjtF15IZPHgwb775Jq1ateKww8JAhWeccQajR4/miiuuoF69epxzzjn079+f/v37AzBjxoxfXGpj/fr1ZGZmAqEpqH404MhJJ51E3759uf3223n22Wd3jBWxbt06hg4dyqJFizAzsrKyio373Xff5b///S8AZ555JiNHjtwxb9CgQdSoUYMOHTrw/fffA9CjRw+GDx9OVlYWgwYNonPnzuXZbFIe7jBtWuhUzsgIo4hde224tlHz5nFHJylEzUqRr6LEAOH2q7WV916W58QhM6NWrVq8//77DBkyhBdeeIF+/fqFWLZvZ9asWcydO5e5c+eyfPlyGjZsCLAj6QC0aNGCXXbZhU8++YQJEyZw8sknA3DDDTfQp08f5s2bx/PPP8+WLVvKFXvdunV33PdoIJYjjjiCN954gxYtWjBs2DAef/zxcr2HlMH27aGpqFs3OO64MHzm//t/oW/httuUGKTUlBwibZtCjeg3u4aFx+XVq1cvJk+ezKZNm9i4cSOTJk2iV69eLF26lHfffReAp59+msMPP5zMzEzWrVvHcccdxz333MPHH38MwNFHH8299967Y51zixjRbvDgwdx5552sW7eOTp06AaHm0KJFCwDGjRu3Y9miLhves2dPxo8fD8BTTz1Fr169iiznkiVLaN68Oeeddx7nnnsuc+bMKXrDSMXZvh3+859w5dPBgyEzE8aOha++giuvDBfBEykDJYfI2AGwT1OoaeG2Ik706dq1K8OGDePggw/mkEMO4dxzz6Vp06bsv//+3H///bRv3561a9dy4YUXsmHDBvr370+nTp04/PDDufvuuwEYPXo0s2fPplOnTnTo0IExY8YU+n6DBg1i/PjxnHTSSTueGzlyJNdeey1dunQhOzt7x/N9+vRhwYIFOzqkE9177708+uijdOrUiSeeeIJ//vOfRZYzIyODgw46iC5dujBhwgQuv/zysmwuKY3t22H8+DBk5kknQVYWPPlkGFry7LN14Tspv8I6I1JpqogO6XSQzOVQh3QF2b7dfcoU906d3MH9wAPdx493z86uuPcohPZJcqqsDmnVHERSxauvhlHUBg4Mw24+/TR88gmcfHK43IVIBVJyEEl2H34YRlHr2xdWrICHHgrNR6eeCjX0FZbKkdafLI+OppF4aT+U0TffwOmnh7ESPvkE/vGPcELbueeGsZhFKlHaJod69eqxZs0a/TDFzN1Zs2YN9erVizuU1LF2bbi+0f77h+E3r702nM18+eWg7ShVJG1PgmvZsiXLli1j1apVxS67ZcuWtPjxStZy1KtXj5YtW8YdRvLLzg5XSL3xxpAghg6FW24Jl70QqWJpmxxq167N3nvvXaJlMzIy6NKlSyVHVPnSpRzV0quvhprB/Plw5JFw991hpDWRmKRts5JISli8GE44IXQ2b9oUznKeMUOJQWKn5CASh82b4YYbwngK06fD3/8ejkAaNEjjMktSSNtmJZGk9cIL4UJ4X38dRly74w7Yc8+4oxL5BdUcRKrKkiWhZnD88eGoo5kz4YknlBgkKSk5iFQyy8kJl9Bu3z40Id1xB8ydC717xx2aSKHUrCRSmWbPpusFF4TzFAYOhNGjYa+94o5KpFiqOYhUho0bwyWzDzmEOmvXwnPPweTJSgySMpIuOZjZzWa23MzmRtNxccckUiqvvAIdO8I998B55/HBuHFhrAWRFJJ0ySFyj7t3jqaX4g5GpETWrYPzzoNjjoG6deGNN2DMGLKjkftEUkmyJgeR1DJtWqgtjB0LV18dOpyLGUFPJJlZsl2YzsxuBoYB64HZwJ/cPd+IzmY2AhgB0Lx58265w1qWRWZm5o5xmVNZupQDUqcstTIz2eeBB9jjf/9jY+vWfHb11Wxo3/4Xy6RKWYqTLuUAlSVXnz59PnT37gXOLGwUoMqcgBnAvAKmgUBzoCahVnMbMLa49RU0ElxppMuoUOlSDvcUKcurr7q3auVeo4b7Nde4b95c4GIpUZYSSJdyuKssuShiJLhYDmV1974lWc7MHgJeqORwREpn82a47rowvkK7dvDuu3DwwXFHJVKhkq7Pwcz2SHh4AqFGIZIc5syBbt1CYrjkEvjoIyUGSUvJeBLcnWbWGXDgG+D8WKMRAcjJCWc233QT7LZb6IA+5pi4oxKpNEmXHNz9zLhjEPmFZcvCBfJefx1OOgn+9S/YZZe4oxKpVEmXHESSyqRJcM45sG0bjBsHZ52lS2pLtZB0fQ4iSWHTJrjwwnBmc9u2oW9h6FAlBqk2lBxE8lqwIHQyjxkDf/4zvPMO7Ldf3FGJVCk1K4kkeuwxuOgiaNgQXn4Zjj467ohEYqGagwiEq6iefTYMGwaHHBIuf6HEINWYkoNIbjPSY4/BjTeGAXn22KP414mkMTUrSfX21FMwYkRoRnrlFehbopP3RdKeag5SPW3bBpdeGs5f6NYtNCMpMYjsoOQg1c/y5WH85vvuC6O1vfqqmpFE8lCzklQvGRlw8smhA3rChHDGs4jko5qDVA/uYdjOvn2haVN4/30lBpEiKDlI+tu8OZzdfOWVcPzxITF06BB3VCJJTclB0tu338IRR8ATT8Bf/wrPPQc77xx3VCJJT30Okr7eegtOPDFcJ2nyZBg4MO6IRFKGag6Snh56CI48MtQS3ntPiUGklJQcJL3k5MAf/xhObOvTR/0LImWk5CDpY/16GDAgDOF52WXw4ovhyCQRKTX1OUh6+PrrcCTSZ5+FkdouuCDuiERSmpKDpL6334ZBgyA7O1xm+6ij4o5IJOWpWUlS24QJoeO5adPQ8azEIFIhlBwkNbnD7bfDKaeE8RfefRfatYs7KpG0oWYlST1ZWWG0tocfhlNPhUcfhbp1445KJK2o5iCpZf166N8/JIbrr4cnn1RiEKkEqjlI6lixAo49FubPD8nhnHPijkgkbSk5SGpYuBD69YMff4SXXtL4ziKVTMlBkt8774RzGGrXhtdfh65d445IJO2pz0GS25Qp4fDUZs1CklBiEKkSsSQHM/uDmc03s+1m1j3PvGvN7Esz+9zMjokjPkkSDz4IgwfDr38dTnRr2zbuiESqjbhqDvOAwcAbiU+aWQfgFOBAoB/wgJnVrPrwJFbutH7iiXAJjH79YOZM2G23uKMSqVZiSQ7uvtDdPy9g1kBgvLtvdfevgS+Bg6s2OonV9u1w5ZXsPXYsnHFGGIdhp53ijkqk2km2DukWwKyEx8ui5/IxsxHACIDmzZuTkZFR5jfNzMws1+uTRaqXw7Kz2X/UKH71yit8NWAAS88+OzQnpbhU3y+50qUcoLKURKUlBzObAfyqgFnXu/uU8q7f3f8N/Buge/fu3rt37zKvKyMjg/K8PlmkdDk2b4aTT4ZXXoFbbmHpYYfRu0+fuKOqECm9XxKkSzlAZSmJSksO7t63DC9bDrRKeNwyek7S2fr14VDVN9+EBx6ACy+ENPlXJ5Kqku1Q1qnAKWZW18z2BvYD3o85JqlMa9aEQ1XfeQeefjokBhGJXSx9DmZ2AnAvsBvwopnNdfdj3H2+mT0LLACygYvdPSeOGKUKrFwJv/sdLFoUOp5///u4IxKRSCzJwd0nAZMKmXcbcFvVRiRVbulS6Ns3XC/ppZfCmAwikjSS7WglqQ6+/DI0Ja1bFzqge/aMOyIRyUPJQarWggWhxrBtG7z2mi6HIZKklByk6nz6aagx1KwZLqB34IFxRyQihUi2o5UkXc2dC336hCurZmQoMYgkOSUHqXxz5oQO5wYNQo1h//3jjkhEiqHkIJXrgw9CU9LOO4fEsO++cUckIiWg5CCVZ9as0Pm8yy4hMey9d9wRiUgJKTlI5XjvvTCU5+67h8TQunXcEYlIKSg5SMX74IOfE8PMmdCyZdwRiUgpKTlIxfrww5AYmjVTYhBJYUoOUnHmzAnXSmrSJCSGVq2KfYmIJCclB6kYc+eGzudGjUJiUB+DSEpTcpDymz8/JIaGDcMJbm3axB2RiJSTkoOUzxdfhPMY6tQJ10rS4aoiaUHXVpKy+/rrkBi2bw9NSTrBTSRtKDlI2SxbFhLDxo0hMbRvH3dEIlKBim1WMrNLzaxpVQQjKWLlypAYVq+Gl1+Ggw6KOyIRqWAl6XNoDnxgZs+aWT8zs8oOSpLYjz+Gw1WXLQsjuPXoEXdEIlIJik0O7v4XYD/gEWAYsMjM/s/M9qnk2CTZbNgAxx4bOqGnToXDD487IhGpJCU6WsndHVgZTdlAU2Cimd1ZibFJMtmyBQYNCmdAT5gQmpVEJG0V2yFtZpcDZwGrgYeBP7t7lpnVABYBIys3RIldVhaccko4VPXxx0OSEJG0VpKjlXYBBrv7ksQn3X27mfWvnLAkaWzfDsOHw5QpcO+9cOaZcUckIlWg2OTg7jcVMW9hxYYjScUdLr8cnnwSbr0VLrkk7ohEpIroDGkp3C23wH33wZ/+BNddF3c0IlKFlBykYGPGwE03wdChMGoU6AhmkWpFyUHymzgRLroI+veHhx5SYhCphpQc5Jdeew1OPx169gyHrNauHXdEIhKDWJKDmf3BzOab2XYz657wfBsz22xmc6NpTBzxVVtz5sDAgdCuHTz/PDRoEHdEIhKTuC68Nw8YDDxYwLzF7t65asMRFi8OZz83awbTpkFTXU5LpDqLJTnkHgKryzQliR9+gGOOgZyccCG9Fi3ijkhEYmbhyhgxvblZBnCVu8+OHrcB5gNfAOuBv7j7m4W8dgQwAqB58+bdxo8fX+Y4MjMzadiwYZlfnyzKUo4amzfT+cor2enrr/n47rtZ36FDJUVXOumyTyB9ypIu5QCVJVefPn0+dPfuBc5090qZgBmE5qO808CEZTKA7gmP6wLNovvdgG+BnYt7r27dunl5zJw5s1yvTxalLkdWlvvvf+9eo4b7lCmVElNZpcs+cU+fsqRLOdxVllzAbC/kd7XSmpXcvW8ZXrMV2Brd/9DMFgPtgNkVHJ64w4UXwosvhnMaBgyIOyIRSSJJdSirme1mZjWj+20Jlwr/Kt6o0tTf/gYPPwzXXw/nnx93NCKSZOI6lPUEM1sG/AZ40cxejmYdAXxiZnOBicAF7v5jHDGmtUcfhZtvDmc/33JL3NGISBKK62ilScCkAp5/Dniu6iOqRmbMgBEjwmhuOvtZRAqRVM1KUsnmzYMTT4T27eE//9HZzyJSKCWH6mLFCjjuONhpp9AJ3bhx3BGJSBKL6wxpqUqZmXD88fDjj/Dmm9CqVdwRiUiSU3JId9nZYYjPuXNh6lTo0iXuiEQkBSg5pLsrrwzNSPffD7//fdzRiEiKUJ9DOrv//jDu8x//GMZnEBEpISWHdDVtGlx2WehrGDUq7mhEJMUoOaSj+fPh5JPh17+Gp5+GmjXjjkhEUoySQ7r54YcwvGeDBmHAnjS58qSIVC11SKeRGtu2wQknwMqV8MYbOmRVRMpMySFduNPurrvgnXfC2c89esQdkYikMDUrpYs77+RX06eHC+kNGRJ3NCKS4pQc0sHzz8O11/L9kUeGS3CLiJSTkkOqmzcPTjsNunbl8z//WVdZFZEKoeSQylavDiO4NWoEU6awvV69uCMSkTSh5JCqtm0LfQsrVsDkydCiRdwRiUga0dFKqeryy+H11+HJJ+Hgg+OORkTSjGoOqejBB2HMGLj6ajj99LijEZE0pOSQat5+Gy69FI49Fm67Le5oRCRNKTmkkmXLwjCfbdromkkiUqnU55AqtmyBwYNh40Z47TVo0iTuiEQkjSk5pAJ3uOAC+OCDcGRShw5xRyQiaU7NSqng3nvhscfg5pth4MC4oxGRakDJIdm9/noY6nPgQLjhhrijEZFqQskhmS1fDiedBPvuC48/DjW0u0SkaqjPIVlt3RrOgN60CTIyYOed445IRKoRJYdkdcUVMGsWTJwI7dvHHY2IVDOxtFOY2Sgz+8zMPjGzSWbWJGHetWb2pZl9bmbHxBFf7MaO/fkM6BNPjDsaEamG4mrEng50dPdOwBfAtQBm1gE4BTgQ6Ac8YGbV60yv2bPhoougb1+49da4oxGRaiqW5ODur7h7dvRwFtAyuj8QGO/uW939a+BLoPpcVW716lBTaN4cnnkGaqnVT0TiYe4ebwBmzwMT3P1JM7sPmOXuT0bzHgH+5+4TC3jdCGAEQPPmzbuNHz++zDFkZmbSsGHDMr++QuTk0Onaa2kydy4f3XsvG/bfv9SrSIpyVBCVJfmkSzlAZcnVp0+fD929e4Ez3b1SJmAGMK+AaWDCMtcDk/g5Sd0HnJEw/xFgSHHv1a1bNy+PmTNnluv1FeLmm93B/cEHy7yKpChHBVFZkk+6lMNdZckFzPZCflcrrd3C3fsWNd/MhgH9gaOiIAGWA60SFmsZPZfepk2Dv/4VzjoLzjsv7mhERGI7WqkfMBIY4O6bEmZNBU4xs7pmtjewH/B+HDFWmaVLw5gMHTvCv/6lMaBFJCnE1eN5H1AXmG7hx3CWu1/g7vPN7FlgAZANXOzuOTHFWPm2boU//AGyssL5DA0axB2RiAgQU3Jw932LmHcbUD1GsbnySnj/fXjuOWjXLu5oRER20MV64vLMM/DAA/CnP4VxGkREkoiSQxw+/xxGjICePeHvf487GhGRfJQcqtqmTaGfoW5dmDABateOOyIRkXx0Cm5Vu+wy+PRTeOklaNmy+OVFRGKgmkNVeuIJeOQRuO46OPbYuKMRESmUkkNVWbAgjAP929+GE95ERJKYkkNV2Lgx9DPstBM8/bQuqCciSU+/UlXhsstg4UJ4+WXYc8+4oxERKZZqDpXt6afD4D3XXw+/+13c0YiIlIiSQ2VatAjOPx8OPxxuuinuaERESkzJobJs3QqnnAJ16qifQURSjn6xKsvIkTBnDkyZAq1aFb+8iEgSUc2hMkyZAqNHw+WXw4ABcUcjIlJqSg4V7dtv4eyzoWtXuOOOuKMRESkTJYeKlJMTBu7JygrXTapbN+6IRETKRH0OFem22+DNN8NlMvYtdMgKEZGkp5pDRXnrrXBZjDPOCJOISApTcqgIa9eG5qQ2beD+++OORkSk3NSsVF7u4US3FSvg7bdh553jjkhEpNyUHMpr7Fj4z3/CiG4HHxx3NCIiFULNSuXx2WfhonpHHhlOehMRSRNKDmW1bRucdhrUrx+OTqqhTSki6UPNSmV1ww3w0UcwaZIuwy0iaUd/d8siIwNGjYLzzoNBg+KORkSkwik5lNbatXDmmeEkt3vuiTsaEZFKoWal0nAP40CvXAnvvBOG/RQRSUOx1BzMbJSZfWZmn5jZJDNrEj3fxsw2m9ncaBoTR3yFeuIJePbZcCZ0jx5xRyMiUmnialaaDnR0907AF8C1CfMWu3vnaLognvAK8PXXcMkl0KsXXH113NGIiFSqWJKDu7/i7tnRw1lAyzjiKLHs7HC9JLNQe6hZM+6IREQqVTJ0SA8H/pfweG8z+8jMXjezXnEF9Qt33hn6GO6/H1q3jjsaEZFKZ+5eOSs2mwH8qoBZ17v7lGiZ64HuwGB3dzOrCzR09zVm1g2YDBzo7usLWP8IYARA8+bNu40fP77MsWZmZtKwYcMC5zVctIiuF17I6l69WHDjjaH2kKSKKkeqUVmST7qUA1SWXH369PnQ3bsXONPdY5mAYcC7QIMilskAuhe3rm7dunl5zJw5s+AZmza5d+jgvsce7mvWlOs9qkKh5UhBKkvySZdyuKssuYDZXsjvaiyHsppZP2Ak8Ft335Tw/G7Aj+6eY2Ztgf2Ar+KIEYDrroMFC2DaNNhll9jCEBGpanGd53AfUBeYbqGZZpaHI5OOAP5mZlnAduACd/8xlghffRX+8Q+4+GI45phYQhARiUssycHdCxxD092fA56r4nDy++knGDYM2rULndEiItWMzpAuyCWXwHffhSOUGjSIOxoRkSqXDIeyJpeJE+Gpp+Avf9HgPSJSbSk5JPr+e7jwQujWDa6/Pu5oRERio+SQK/eiehs2wGOPQe3acUckIhIb9TnkevJJmDw5jNNw4IFxRyMiEivVHIC6q1bBpZfCYYfBH/8YdzgiIrFTcnBn/zvvhKwsGDdOF9UTEaGaJ4el6+Afwx9kl9mzGT1kFEt3K/D0CxGRaqdaJ4frH1rMeU9fxdJO3bmn64UMnxp3RCIiyaFad0h/ux5m730YX51/DtuzjK/Wxh2RiEhyqNY1h1r77cOwi14mc9fdqWHQtmncEYmIJIdqnRzGDoB9moIRbscOiDsiEZHkUK2Tw16NYcaZ0HH3cLtX47gjEhFJDtU6OYiISMGUHEREJB8lBxERyUfJQURE8lFyEBGRfJQcREQkH3P3uGMoNzNbBSwpxyp2BVZXUDhxSpdygMqSjNKlHKCy5Grt7rsVNCMtkkN5mdlsd+8edxzllS7lAJUlGaVLOUBlKQk1K4mISD5KDiIiko+SQ/DvuAOoIOlSDlBZklG6lANUlmKpz0FERPJRzUFERPJRchARkXyqZXIws1Fm9pmZfWJmk8ysSSHL9TOzz83sSzO7porDLJaZ/cHM5pvZdjMr9FA2M/vGzD41s7lmNrsqYyypUpQlqfcJgJntYmbTzWxRdFvgMFJmlhPtk7lmljSD1Ba3jc2srplNiOa/Z2ZtYgizREpQlmFmtiphP5wbR5zFMbOxZvaDmc0rZL6Z2eionJ+YWddyv6m7V7sJOBqoFd2/A7ijgGVqAouBtkAd4GOgQ9yx54mxPbA/kAF0L2K5b4Bd4463vGVJhX0SxXkncE10/5qCPl/RvMy4Yy3LNgYuAsZE908BJsQddznKMgy4L+5YS1CWI4CuwLxC5h8H/I8wdtmhwHvlfc9qWXNw91fcPTt6OAtoWcBiBwNfuvtX7r4NGA8MrKoYS8LdF7r753HHURFKWJak3yeRgcBj0f3HgEHxhVJqJdnGieWbCBxlZlaFMZZUqnxeiuXubwA/FrHIQOBxD2YBTcxsj/K8Z7VMDnkMJ2TcvFoA3yY8XhY9l4oceMXMPjSzEXEHUw6psk+au/t30f2VQPNClqtnZrPNbJaZDaqa0IpVkm28Y5noT9Y6oFmVRFc6Jf28nBg1xUw0s1ZVE1qFq/DvRq1yhZPEzGwG8KsCZl3v7lOiZa4HsoGnqjK20ihJOUrgcHdfbma7A9PN7LPon0iVqqCyJIWiypL4wN3dzAo7Xrx1tF/aAq+Z2afuvriiY5UiPQ884+5bzex8Qo3oyJhjSgppmxzcvW9R881sGNAfOMqjRrs8lgOJ/yJaRs9VqeLKUcJ1LI9ufzCzSYTqdpUnhwooS1LsEyi6LGb2vZnt4e7fRVX7HwpZR+5++crMMoAuhDbyOJVkG+cus8zMagGNgTVVE16pFFsWd0+M+2FCf1EqqvDvRrVsVjKzfsBIYIC7bypksQ+A/cxsbzOrQ+h4S5ojSkrKzHYys0a59wmd8QUe8ZACUmWfTAWGRveHAvlqRWbW1MzqRvd3BQ4DFlRZhIUryTZOLN8Q4LVC/mDFrdiy5GmXHwAsrML4KtJU4KzoqKVDgXUJTZtlE3cvfBwT8CWhfW5uNOUeebEn8FLCcscBXxD+zV0fd9wFlOMEQtviVuB74OW85SAcqfFxNM1PxnKUtCypsE+iGJsBrwKLgBnALtHz3YGHo/s9gU+j/fIpcE7ccRe1jYG/Ef5MAdQD/hN9j94H2sYdcznK8vfoe/ExMBM4IO6YCynHM8B3QFb0PTkHuAC4IJpvwP1ROT+liKMXSzrp8hkiIpJPtWxWEhGRoik5iIhIPkoOIiKSj5KDiIjko+QgIiL5KDmIVBIza2JmF8Udh0hZKDmIVJ4mhCuYiqQcJQeRynM7sE80TsCouIMRKQ2dBCdSSaJBcF5w945xxyJSWqo5iIhIPkoOIiKSj5KDSOXZADSKOwiRslByEKkkHsYKeNvM5qlDWlKNOqRFRCQf1RxERCQfJQcREclHyUFERPJRchARkXyUHEREJB8lBxERyUfJQURE8vn/HUnvJQUNaUMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "# Fit the model to a dense 1D-grid of points from [t_min, t_max].\n",
    "fitn = m*20\n",
    "fitx = np.linspace(np.min(t), np.max(t), fitn)  # Actually t.\n",
    "fity = vft(fitx)\n",
    "\n",
    "# Plot the fitted model against observations.\n",
    "plt.scatter(t, y, s=16, c='dodgerblue', label='observations')\n",
    "plt.plot(fitx, fity, c='red', label='fit')\n",
    "plt.title('Observations vs. Fit, Polynomial Interpolation')\n",
    "plt.ylabel('y')\n",
    "plt.xlabel('t')\n",
    "plt.grid()\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 07.03 Lagrange and Newton Interpolation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Lagrange Basis\n",
    "The Vandermonde matrix for the Lagrange basis is the identity matrix.\n",
    "$$\n",
    "l_j(t) = \\prod_{k=1,k \\neq j}^n (t - t_k) / \\prod_{k=1,k \\neq j}^n (t_j - t_k) \\quad j = 1,\\cdots,n\n",
    "$$\n",
    "\n",
    "Polynomial evaluated at $t_i$ with this basis:\n",
    "$$\n",
    "p_{n-1}(t_i) = y_1 l_1(t_i) + y_2 l_2(t_i) + \\cdots + y_n l_n(t_i)\n",
    "$$\n",
    "\n",
    "Lagrange form of a polynomial obtained for $m$ data points will give same result as $m-1$ degree polynomial obtained using monomial basis.\n",
    "\n",
    "Cost of computing the weights is $O(n^2)$ (less than solving linear system as is the case for monomial basis) and cost of evaluating the interpolant is $O(n)$, but the constant is larger $5n$ in comparison to methods using Horner's rule.\n",
    "\n",
    "#### Newton Interpolation\n",
    "Newton basis functions given by:\n",
    "$$\n",
    "\\pi_j(t) = \\prod_{k=1}^{j-1} (t - t_k) \\quad j = 1,\\cdots,n\n",
    "$$\n",
    "\n",
    "Polynomial evaluated at $t_i$ with this basis:\n",
    "$$\n",
    "p_{n-1}(t_i) = x_1 + x_2(t - t_1) + x_3(t - t_1)(t - t_2) + \\cdots + x_n(t - t_1)(t - t_2)\\cdots(t - t_{n-1})\n",
    "$$\n",
    "* Avoid repeated computations using Horner's rule, cost of evaluation will be $2n$.\n",
    "\n",
    "The basis matrix $A$ formed from the linear combination of these basis functions is lower triangular and can be solved using forward substitution (cost = $O(n^2)$).\n",
    "\n",
    "Newton basis functions can be constructed incrementally.\n",
    "* Given $m$ data points and an interpolating polynomial of degree $m-1$, then a new interpolating polynomial of degree $m + n - 1$ can be constructed for any additional $n$ data points using Newton interpolation "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use Lagrange interpolation to interpolate 3 ($m = 3$) data points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.interpolate as interp\n",
    "\n",
    "# Input data $(t_i, y_i)$.\n",
    "m = 3\n",
    "t = np.array([-2,0,1], dtype='d')\n",
    "y = np.array([-27,-1,0], dtype='d')\n",
    "\n",
    "# Implement L1, L2, and L3.\n",
    "# L1 = \\frac{(x-t_2)(x-t_3)}{(t_1-t_2)(t_1-t_3)}\n",
    "L1 = lambda x: ((x-t[1])*(x-t[2]))/((t[0]-t[1])*(t[0]-t[2]))\n",
    "# L2 = \\frac{(x-t_1)(x-t_3)}{(t_2-t_1)(t_2-t_3)}\n",
    "L2 = lambda x: ((x-t[0])*(x-t[2]))/((t[1]-t[0])*(t[1]-t[2]))\n",
    "# L3 = \\frac{(x-t_1)(x-t_2)}{(t_3-t_1)(t_3-t_2)}\n",
    "L3 = lambda x: ((x-t[0])*(x-t[1]))/((t[2]-t[0])*(t[2]-t[1]))\n",
    "\n",
    "# Implement the polynomial p_2(x).\n",
    "p2x = lambda x: np.dot(y, np.array([L1(x), L2(x), L3(x)]))\n",
    "\n",
    "# Evaluate the polynomial at the data points.\n",
    "vp2x = np.vectorize(p2x)\n",
    "yhat = vp2x(t)\n",
    "np.testing.assert_almost_equal(y - yhat, np.zeros(m))\n",
    "\n",
    "# Compare to polynomial returned from scipy.interpolate.lagrange.\n",
    "poly = interp.lagrange(t, y)\n",
    "x = np.linspace(np.min(t), np.max(t), num=100)\n",
    "np.testing.assert_almost_equal(vp2x(x), poly(x))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the observations against the fit using Lagrange interpolation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAEWCAYAAACNJFuYAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAxgElEQVR4nO3deZgU1fX/8fdhF0FANEQBQVQQJMiqBsUvKFEwBBCJO7KouGtiDG6JmkR/cYsaXIJREcUFDIZFY1BQxh0VlaiAGygIiAIiMKwzcH5/3BpsZ996qrvn83qeenqp6qpzq7r7VN1bVdfcHRERkUQ14g5ARERSj5KDiIgUoOQgIiIFKDmIiEgBSg4iIlKAkoOIiBSg5BAjM7vBzB6LO46yMrP/mtnwuONIBjPrZWafxB2HJIeZTTCzGyvw+Wwza1OZMaUqJYckMrMRZvahmW02s1Vm9g8zaxx3XGVRWAJz9/7u/khcMVUGM2ttZh792POG/7n7q+7eLmG6L82sbxnmm5YJvyLMLMvMzinltBX6c65KhZXL3Ru4+5K4YqpKSg5JYma/A24Bfg80Ao4AWgGzzKxOFcZRq6qWlaYaRz/4Bu5+aNzBlIYF1fK3q+9z1amWX7BkM7M9gD8Bl7j7THfPcfcvgZOB1sCZCZPXM7PJZrbRzN4zs0MT5nOlma2Ixn1iZsdG79cws6vMbLGZrTWzp8xsz2hc3h7x2Wa2DHgpqga6OF+M/zOzIdHzv5vZV2a2wczeNbNe0fv9gGuAU/L2rKP3d+1RRbH8wcyWmtm3ZvaomTXKF8twM1tmZmvM7NqEGA4zs3nRcr8xszuKWJ+LzGxAwutaZrbazLqaWT0zeyxaD9+b2Ttm1qwcmy1v3r3NbHn0fCKwH/BMVP4x5Z1vNL+8bbbRzBaa2YkJ42qa2d+idfSFmV0crbta0fgsM7vJzF4HNgNtzGxktG42mtkSMzsvfznM7HfRdvnazEYmjG9qZs9E6/4dM7vRzF5LGH+wmc0ys++i797JZVl/hS3XzEYDZwBjovX5TPT+vmb2dLRNvzCzSxPmd4OZTYm28QZgRMJ7Rf1u2kfr63szW2BmA4uItYmZPRstd130vEU07iagF3BPFOs90ftuZgdGzxtF3/fV0ff/DxYlbQu1Bq+Z2e3RvL8ws/6lWYcpw901VPIA9ANygVqFjHsEeDJ6fgOQAwwFagNXAF9Ez9sBXwH7RtO2Bg6Inl8GzAVaAHWB+xPm2Rpw4FFgd2A34Czg9YQYOgDfA3Wj12cCTYFawO+AVUC9hBgfy1eGLOCc6Pko4HOgDdAA+DcwMV8sD0RxHApsA9pH498EhkXPGwBHFLE+rwMeT3j9S2BR9Pw84BmgPlAT6AbsUYptlBdbrXzv9waWJ7z+Euhbhm1fYH0ljPs1sC9hp+wUYBOwTzTufGBhtE2bALMT44vW+TLgkGg71Y7WwwGAAf9HSBpdE8qRC/w5mvaEaHyTaPykaKgffR++Al6Lxu0evR4ZLasLsAboUES5Er8PJS13AnBjwmdrAO9G27hO9D1aAhyf7zcyOJp2N4r/3dQmfB+vieZ3DLARaJd/+YTv/EnROmgI/AuYVli5Et5z4MDo+aPA9OizrYFPgbOjcSOiGM8lfC8vAFYCFvf/U6m/y3EHkIkD4c92VRHjbgZmRc9vAOYmjKsBfE3YYzkQ+BboC9TON49FwLEJr/eJvoi1+OFPr03C+IaEP6JW0eubgPHFxL8OODQhxuKSw4vAhQnj2hUSS4uE8W8Dp0bPXyEcYe1Vwvo8MPqB149ePw5cFz0fBbwBdCrjNsqL7fuE4QqSmBwKmXY+MCh6/hJwXsK4vhRMDn8uYX7TgMui572BLSQkv+j7dAThzyqH6A8zGncjPySHU4BX8837fuD6Ipab+H0ocrnR8wn8ODkcDizLN7+rgYcT1ucrhazjon43vQg7NzUSxj8J3FDY8vPNtzOwrrByJbzn0fexJrCdhIRJ2FHJip6PAD5PGFc/+uxPy/I9jXNQtVJyrAH2ssLrR/eJxuf5Ku+Ju+8ElhOOFj4HfkP4IXxrZpPMbN9o0lbA1Oiw+XtCstgBNCtivhuB/wCnRm+dRviDBcDMroiqJ9ZH82sE7FXKsu4LLE14vZSQGBJjWZXwfDPhKAHgbKAt8HFUtTGAQkTrYhHwKzOrDwwEnohGTwSeByaZ2Uozu9XMapcydgiJqXE03F6Gz5WZmZ1lZvMTtltHfljP+5KwzfI9L/Q9M+tvZnOjqp/vCXvpidttrbvnJrzOW/d7E7ZRUctrBRyeF2c07zOAn5aupEUutzCtgH3zLesaivguF/Ze4u8mGr6K3suzFGiefwZmVt/M7o+qhDYQdlYam1nNEksY1nNtCn73E5ez63vv7pujp0Wth5Sj5JAcbxKqT4YkvmlmDYD+hL3tPC0TxtcgVCusBHD3J9z9KMIPyAkN3BB+GP0T/tQau3s9d1+RMN/8t9t9EjjNzH4O1APmRMvsBYwhtIc0cffGwHpCVUVh88lvZRRfnv0I1QrflPA53P0zdz8N+ElUtilmtnsRkz9JSGqDgIVRwsBDe86f3L0D0BMYQKhGqyyVcttiM2tFqF67GGgareeP+GE9f03Y9nlaUtCuWMysLvA0cDvQLJrfcwnzK85qwjYqanlfAS/n+341cPcLSjHvkuRfn18BX+RbVkN3P6GYz/wo3ny/m5VAS/txg/1+wAoK+h3hSPdwd98DODpvlsUsN88awtFX/u9+YctJS0oOSeDu6wnVJXebWT8zq21mrYGnCHs4ExMm72ZmQ6KjjN8QkspcM2tnZsdEfwJbCYfqeXtD44Cboj8czGxvMxtUQljPEb7IfwYmJ+xZNST8UawGapnZdcAeCZ/7BmhtRZ8d8yTwWzPbP0p+/y+af24R0+9iZmea2d5RLN9Hb+8sYvJJwHGEutu8owbMrI+Z/Sza29tA+MEWNY/y+IZQD54Y95dmNqKYz9Sw0FCeN9Ql1OM7YT0TNdJ2TPjMU8BlZtbcwunOV5YQVx1Ce9NqIDdq7DyuNAVy9x2EtqEbor3ng/lxQn0WaGtmw6Lvbm0z62Fm7Usz/xLkX59vAxstnHyxm4WG+Y5m1qOE+RT6uwHeIhypjIni7g38ivD9ya8h4Xf1vYUTOq4vIdZdonX4FOF32DD6LV4OZMxpzEoOSeLutxIOj28n/Gm9RdhLOtbdtyVMOp1Qx7sOGAYMcfccwg//ZsIeyirC3vXV0Wf+DswAXjCzjYQfxeElxLON8IfQl4Q/V0KVzExCY9pSQiJKPIz/V/S41szeK2TW4wnJ7hVCo+BW4JLiYknQD1hgZtlRmU519y1FxP814YisJzA5YdRPgSmEdbwIeDmKBzMbZ2bjShlLUf4K/CGq8rjCwmnITQnrvCinEf508obF7r4Q+FtUhm+AnwGvJ3zmAeAF4APgfUIyzyVUFxYQVRVeSviDWgecTvhOlNbFhOrDVYT19SThDzZv3scRqiFXRtPcQvhOVtRDQIdofU6L/mQHEOr7vyB83x+MYitOob8bd99OSAb9o3ndB5zl7h8XMo+7CA3cawjbc2a+8X8HhkZnG40t5POXENrylgCvEX5X40uIO21Y1FgiIqVgZkcBF0XVYclcTn9gnLu3KnHiylneLYTG0uFVsbyKMLMbCGcMnVnStFJ+OnIQKQN3fy0ZiSGqUjnBwjUczQlVHFMrezkJyzvYzDpZcBjh5ICkLU/Sj5KDSGowQjvVOkK10iLCuf/J0pBQzbiJUE33N0JVjQigaiURESmEjhxERKSAjLiJ1V577eWtW7cu9+c3bdrE7rsXdXp9+siUcoDKkooypRygsuR5991317j73oWNy4jk0Lp1a+bNm1fuz2dlZdG7d+/KCygmmVIOUFlSUaaUA1SWPGa2tKhxqlYSEZEClBxERKQAJQcRESkgI9ocCpOTk8Py5cvZunVridM2atSIRYsWVUFUyZWq5ahXrx4tWrSgdu2y3CxVROKUsclh+fLlNGzYkNatW2NW/I0qN27cSMOGDasosuRJxXK4O2vXrmX58uXsv//+cYcjIqWUsdVKW7dupWnTpiUmBkkuM6Np06alOoITkdSRsckBUGJIEdoOIuknY6uVREQKtX07tb//Hj77DNavh++/D49bt8KWLeExb9ixA9x/GABq1YI6daBu3R8e69eHPfb4YWjYEBo3hiZNoEZ67oMrOSTZ2LFj+cc//sGqVau48sorueqqq5g2bRpt27alQ4cOcYcnkjlyc2HlSli2DJYuDY/Ll8M33/x42LCBI6sqppo1Ya+9YO+9w9CsGbRoAS1b/vC4337wk59Aih1hp2xyMLN+hM42agIPuvvNMYdULvfddx+zZ8+mRYsfemScNm0aAwYMUHIQKSv38Kf/8cfw6adh7//TT8OwbFnY00/UpEn4Q27WDLp0CY97781na9ZwUPfuYe++ceOwt1+/PtSrF4bddgtHBDVrhj/txD/uHTtg+/YwbNsWhs2bYeNG2LDhh2HdOli9Gr79NjyuXg3vvANTp4bPJNp9dzjwQDjggPB44IHQoQMcckiILwYpmRyiLh/vBX5B6FbzHTObEfWmlTbOP/98lixZQv/+/Rk1ahSLFy/m9NNPZ8aMGbz88svceOONPP300xxwwAFxhyqSejZuhPffhw8+gA8/DMNHH4X38zRsCG3bwhFHwOmnQ6tWYU+8VauwV17EPYdWZGVxUHlvn1GzZkgeu+1Wvs+7w5o14ajmq69CUlu8GD7/HBYuhGefDYknz777hiTRsSMceih06wYHH8yyTbUYNQMG14UbJ8L4gbBfSf3nlUFKJgfgMOBzd18CYGaTiDqWL9fcfvMbmD+/yNG77dgRNnhZdO4Md91V7CTjxo1j5syZzJkzh2effRaAnj17MnDgQAYMGMDQoUPLtkyRTLV9e0gE77wThnnzYNGiH+r5mzSBn/0MzjorPHboAAcdFI4EUqw6pkRmP1QzdelScPyOHSFpLFwICxaEhLhgAYwbF9pEAHbbjex9D2X4Pl1pcsQ+LG7Zm1EzYPawygszVZNDc37cj/Fy8vWRbGajgdEAzZo1Iysr60czaNSoERujPYy627dTI//hZiJ3cosbX4id27ezLXEPpshZO9nZ2WzdupXt27ezceNGcnJy2LJly674KsuOHTsqfZ6VZevWrQW2UXGys7PLNH0qy5SyVGY5amVns8eCBTT68EMaffghDT/+mJrR3vL2Jk3YcPDBbBw+nOy2bdl44IFs32uvglU7H38chnJIi21Svz706BEGgB07qL98OQ0//ZQGn35KowWfcsq7E1hR+wgu73EUBlRmkVI1OZTI3f8J/BOge/funv+uhIsWLfrhgrD77it2XuW9eKxOKaYxMxo0aEC9evWoU6cODRs2pHbt2uy2226VfsFaKl4El6devXp0KWwvqQi6a2bqqVA5Nm+G116D2bPhxRfDUYJ7OPOnSxe46CI48kg4/HDqNG/OXmbsVanR/1gmbJO+E2HJ2p2MaTKTO1b15oAmcHHvypt/qiaHFUDLhNctovcyQsOGDVN2D1+kUriHNoL//Aeefx7efDNUHdWuDT//OVx/PRx9NBx2WJHtAlK88QNh1Iwa5NatzwE1wuvKlKrJ4R3gIDPbn5AUTgVOjzekynPqqady7rnnMnbsWKZMmaIGackMmzaFI4PnngvD8uXh/c6d4dJL4dhjoVcvJYNKsl+j0MaQlVW5Rwx5UjI5uHuumV0MPE84lXW8uy+IOaxy+fLLLwEYMWIEI0aMAODII49k4cK0OvFKpHDr1sEzz8C//x2OELZuDWcQHXcc/OlP0L8/7LNP3FFKOaRkcgBw9+eA5+KOQ6S6W7aeH50yOeGotbR4YQo8/TTMmRMuPmvRAs49FwYPhqOOClcOS1pL2eQgIqlh1AxY9fVG2n41i2teuo2ffvIC7MgNp5JecQUMGQLdu6ffKaVSLCUHESncjh0waxaX3TmBvh/NYLecLaxo3JKHel/OebefFi7IUkLIWEoOIvJjn3wCEybAo4/CypUc3aApUw4fyW6/OJgxe1xEm6Y1OK9z3EFKsik5iEi48vapp+D++8NppzVrhsbksWNZ32sAjzxfl8F1s2izrUalnzIpqUnJQaQ6W7QoJIRHHgm3rm7XDm69Fc48c9dZRvuR3FMmJTWl543G09iXX35Jx44d4w6D+fPn89xzP5wMNmPGDG6+OS1vfCtllZsbzjTq0yfco+i++6Bfv/Dvv2gR/P73Ov1UdOSQCXJzc6lVq2ybcv78+cybN48TTjgBgIEDBzJwoOoLMtq6dfDgg3DPPeFOoK1bw803w8iRoT8BkQQ6ckiyO+64g44dO9KxY0fuiu7impubyxlnnEH79u0ZOnQomzdvBuCqq66iQ4cOdOrUiSuuuAKA1atXc9JJJ9GjRw969OjB66+/DsANN9zAsGHDOPLIIxk2bBhHHHEEixYt2rXc3r17M2/ePN5++21+/vOf06VLF3r27Mknn3zC9u3bue6665g8eTKdO3dm8uTJTJgwgYsvvhgIRzfHHHMMnTp14thjj2XZsmVAuJDv0ksvpWfPnrRp04YpU6YA8PXXX3P00UfTuXNnOnbsyKuvvlol61ZK6dNP4cILw7UIY8ZAmzahT4HPP4crr1RikMK5e9oP3bp18/wWLlxY4L2ibNiwwZd+737so+77/z08Lv2+1B8v0rx587xjx46enZ3tGzdu9A4dOvh7773ngL/22mvu7j5y5Ei/7bbbfM2aNd62bVvfuXOnu7uvW7fO3d1PO+00f/XVV93dfenSpX7wwQe7u/v111/vXbt29c2bN7u7+x133OFXXnmlu7uvXLnS27Zt6+7u69ev95ycHHd3nzVrlg8ZMsTd3R9++GG/6KKLdsWa+HrAgAE+YcIEd3d/6KGHfNCgQe7uPnz4cB86dKjv2LHDFyxY4AcccIC7u99+++1+4403urt7bm6ub9iwocC6KMv2cHefM2dOmaZPZbGV5a233IcMcTdzr1PHfeRI9/ffL/fstE1SU0XKAszzIv5XVa0UGTUDFq+DnR4eK+Pe6K+99honnngiu0f3khkyZAivvvoqLVu25MgjQ0eFZ555JmPHjuU3v/kN9erV4+yzz2bAgAEMGDAAgNmzZ//oVhsbNmwgOzsbCFVBu0Udjpx88sn07duXm2++maeeempXXxHr169n+PDhfPbZZ5gZOTk5Jcb95ptv8u9//xuAYcOGMWbMmF3jBg8eTI0aNejQoQPffPMNAD169GDUqFHk5OQwePBgOnfuXJHVJhXhDjNnhkblrKzQi9jVV4d7GzVrFnd0kkZUrRRZEiUGCI9L1iVvWZbvwiEzo1atWrz99tsMHTqUZ599ln79+oVYdu5k7ty5zJ8/n/nz57NixQoaNGgAsCvpADRv3pw999yTDz74gMmTJ3PKKacA8Mc//pE+ffrw0Ucf8cwzz7B169YKxV63bt1dzz3qiOXoo4/mlVdeoXnz5owYMYJHH320QsuQcti5M1QVdesGJ5wQus/8299C28JNNykxSJkpOUTaNIEa0X92DQuvK6pXr15MmzaNzZs3s2nTJqZOnUqvXr1YtmwZb775JgBPPPEERx11FNnZ2axfv54TTjiBO++8k//9738AHHfccdx999275jm/mB7thgwZwq233sr69evp1KkTEI4cmjdvDsCECRN2TVvcbcN79uzJpEmTAHj88cfp1atXseVcunQpzZo149xzz+Wcc87hvffeK37FSOXZuRP+9a9w59MhQyA7G8aPhyVL4PLLw03wRMpBySEyfiAc0ARqWnisjAt9unbtyogRIzjssMM4/PDDOeecc2jSpAnt2rXj3nvvpX379qxbt44LLriAjRs3MmDAADp16sRRRx3FHXfcAcDYsWOZN28enTp1okOHDowbN67I5Q0ePJhJkyZx8skn73pvzJgxXH311XTp0oXc3Nxd7/fp04eFCxfuapBOdPfdd/Pwww/TqVMnJk6cyN///vdiy5mVlcWhhx5Kly5dmDx5Mpdddll5VpeUxc6dMGlS6DLz5JMhJwceeyx0LTlypG58JxVXVGNEOg2V0SCdCVK5HGqQriQ7d7pPn+7eqZM7uB9yiPukSe65uZW3jCJom6SmZDVI68hBJF28+GLoRW3QoNDt5hNPwAcfwCmnhNtdiFQiJQeRVPfuu6EXtb59YeVKeOCBUH102mlQQz9hSY6M/mZ5dDaNxEvboZy+/BLOOCP0lfDBB3DXXeGCtnPOCX0xiyRRxiaHevXqsXbtWv0xxczdWbt2LfXq1Ys7lPSxbl24v1G7dqH7zauvDlczX3YZaD1KFcnYi+BatGjB8uXLWb16dYnTbt26NSP+vFK1HPXq1aNFixZxh5H6cnPDHVKvuy4kiOHD4S9/Cbe9EKliGZscateuzf7771+qabOysujSpUuSI0q+TClHtfTii+HIYMECOOYYuOOO0NOaSEwytlpJJC0sXgwnnhgamzdvDlc5z56txCCxU3IQicOWLfDHP4b+FGbNgr/+NZyBNHiw+mWWlJCx1UoiKevZZ8ON8L74IvS4dsstsO++cUcl8iM6chCpKkuXhiODX/0qnHU0Zw5MnKjEIClJyUEkyWzHjnAL7fbtQxXSLbfA/PnQu3fcoYkUSdVKIsk0bx5dzz8/XKcwaBCMHQv77Rd3VCIl0pGDSDJs2hRumX344dRZtw6efhqmTVNikLSRcsnBzG4wsxVmNj8aTog7JpEyeeEF6NgR7rwTzj2XdyZMCH0tiKSRlEsOkTvdvXM0PBd3MCKlsn49nHsuHH881K0Lr7wC48aRG/XcJ5JOUjU5iKSXmTPD0cL48XDllaHBuYQe9ERSmaXajenM7AZgBLABmAf8zt0L9OhsZqOB0QDNmjXrltetZXlkZ2fv6pc5nWVKOSB9ylIrO5sD7ruPff77Xza1asXHV17JxvbtfzRNupSlJJlSDlBZ8vTp0+ddd+9e6MiiegFK5gDMBj4qZBgENANqEo5qbgLGlzS/wnqCK4tM6RUqU8rhniZlefFF95Yt3WvUcL/qKvctWwqdLC3KUgqZUg53lSUPxfQEF8uprO7etzTTmdkDwLNJDkekbLZsgWuuCf0rtG0Lb74Jhx0Wd1QilSrl2hzMbJ+ElycSjihEUsN770G3biExXHwxvP++EoNkpFS8CO5WM+sMOPAlcF6s0YgA7NgRrmy+/nrYe+/QAH388XFHJZI0KZcc3H1Y3DGI/Mjy5eEGeS+/DCefDP/4B+y5Z9xRiSRVyiUHkZQydSqcfTZs3w4TJsBZZ+mW2lItpFybg0hK2LwZLrggXNncpk1oWxg+XIlBqg0lB5H8Fi4MjczjxsHvfw9vvAEHHRR3VCJVStVKIokeeQQuvBAaNIDnn4fjjos7IpFY6MhBBMJdVEeOhBEj4PDDw+0vlBikGlNyEMmrRnrkEbjuutAhzz77lPw5kQymaiWp3h5/HEaPDtVIL7wAfUt18b5IxtORg1RP27fDJZeE6xe6dQvVSEoMIrsoOUj1s2JF6L/5nntCb20vvqhqJJF8VK0k1UtWFpxySmiAnjw5XPEsIgXoyEGqB/fQbWffvtCkCbz9thKDSDGUHCTzbdkSrm6+/HL41a9CYujQIe6oRFKakoNktq++gqOPhokT4U9/gqefhj32iDsqkZSnNgfJXK+9BiedFO6TNG0aDBoUd0QiaUNHDpKZHngAjjkmHCW89ZYSg0gZKTlIZtmxA37723BhW58+al8QKSclB8kcGzbAwIGhC89LL4X//CecmSQiZaY2B8kMX3wRzkT6+OPQU9v558cdkUhaU3KQ9Pf66zB4MOTmhttsH3ts3BGJpD1VK0l6mzw5NDw3aRIanpUYRCqFkoOkJ3e4+WY49dTQ/8Kbb0LbtnFHJZIxVK0k6ScnJ/TW9uCDcNpp8PDDULdu3FGJZBQdOUh62bABBgwIieHaa+Gxx5QYRJJARw6SPlauhP79YcGCkBzOPjvuiEQylpKDpIdFi6BfP/juO3juOfXvLJJkSg6S+t54I1zDULs2vPwydO0ad0QiGU9tDpLapk8Pp6c2bRqShBKDSJWIJTmY2a/NbIGZ7TSz7vnGXW1mn5vZJ2Z2fBzxSYq4/34YMgR+9rNwoVubNnFHJFJtxHXk8BEwBHgl8U0z6wCcChwC9APuM7OaVR+exMqdVhMnhltg9OsHc+bA3nvHHZVItRJLcnD3Re7+SSGjBgGT3H2bu38BfA4cVrXRSax27oTLL2f/8ePhzDNDPwy77x53VCLVTqo1SDcH5ia8Xh69V4CZjQZGAzRr1oysrKxyLzQ7O7tCn08V6V4Oy82l3W238dMXXmDJwIEsGzkyVCeluXTfLnkypRygspRG0pKDmc0GflrIqGvdfXpF5+/u/wT+CdC9e3fv3bt3ueeVlZVFRT6fKtK6HFu2wCmnwAsvwF/+wrIjj6R3nz5xR1Up0nq7JMiUcoDKUhpJSw7u3rccH1sBtEx43SJ6TzLZhg3hVNVXX4X77oMLLoAM2asTSVepdirrDOBUM6trZvsDBwFvxxyTJNPateFU1TfegCeeCIlBRGIXS5uDmZ0I3A3sDfzHzOa7+/HuvsDMngIWArnARe6+I44YpQqsWgW/+AV89lloeP7lL+OOSEQisSQHd58KTC1i3E3ATVUbkVS5Zcugb99wv6Tnngt9MohIyki1s5WkOvj881CVtH59aIDu2TPuiEQkHyUHqVoLF4Yjhu3b4aWXdDsMkRSl5CBV58MPwxFDzZrhBnqHHBJ3RCJShFQ7W0ky1fz50KdPuLNqVpYSg0iKU3KQ5HvvvdDgXL9+OGJo1y7uiESkBEoOklzvvBOqkvbYIySGAw+MOyIRKQUlB0meuXND4/Oee4bEsP/+cUckIqWk5CDJ8dZboSvPn/wkJIZWreKOSETKQMlBKt877/yQGObMgRYt4o5IRMpIyUEq17vvhsTQtKkSg0gaU3KQyvPee+FeSY0bh8TQsmWJHxGR1KTkIJVj/vzQ+NywYUgMamMQSWtKDlJxCxaExNCgQbjArXXruCMSkQpScpCK+fTTcB1DnTrhXkk6XVUkI+jeSlJ+X3wREsPOnaEqSRe4iWQMJQcpn+XLQ2LYtCkkhvbt445IRCpRidVKZnaJmTWpimAkTaxaFRLDmjXw/PNw6KFxRyQilaw0bQ7NgHfM7Ckz62dmluygJIV99104XXX58tCDW48ecUckIklQYnJw9z8ABwEPASOAz8zs/5nZAUmOTVLNxo3Qv39ohJ4xA446Ku6IRCRJSnW2krs7sCoacoEmwBQzuzWJsUkq2boVBg8OV0BPnhyqlUQkY5XYIG1mlwFnAWuAB4Hfu3uOmdUAPgPGJDdEiV1ODpx6ajhV9dFHQ5IQkYxWmrOV9gSGuPvSxDfdfaeZDUhOWJIydu6EUaNg+nS4+24YNizuiESkCpSYHNz9+mLGLarccCSluMNll8Fjj8GNN8LFF8cdkYhUEV0hLUX7y1/gnnvgd7+Da66JOxoRqUJKDlK4cePg+uth+HC47TbQGcwi1YqSgxQ0ZQpceCEMGAAPPKDEIFINKTnIj730EpxxBvTsGU5ZrV077ohEJAaxJAcz+7WZLTCznWbWPeH91ma2xczmR8O4OOKrtt57DwYNgrZt4ZlnoH79uCMSkZjEdeO9j4AhwP2FjFvs7p2rNhxh8eJw9XPTpjBzJjTR7bREqrNYkkPeKbC6TVOK+PZbOP542LEj3EivefO4IxKRmFm4M0ZMCzfLAq5w93nR69bAAuBTYAPwB3d/tYjPjgZGAzRr1qzbpEmTyh1HdnY2DRo0KPfnU0V5ylFjyxY6X345u3/xBf+74w42dOiQpOjKJlO2CWROWTKlHKCy5OnTp8+77t690JHunpQBmE2oPso/DEqYJgvonvC6LtA0et4N+ArYo6RldevWzStizpw5Ffp8qihzOXJy3H/5S/caNdynT09KTOWVKdvEPXPKkinlcFdZ8gDzvIj/1aRVK7l733J8ZhuwLXr+rpktBtoC8yo5PHGHCy6A//wnXNMwcGDcEYlICkmpU1nNbG8zqxk9b0O4VfiSeKPKUH/+Mzz4IFx7LZx3XtzRiEiKietU1hPNbDnwc+A/ZvZ8NOpo4AMzmw9MAc539+/iiDGjPfww3HBDuPr5L3+JOxoRSUFxna00FZhayPtPA09XfUTVyOzZMHp06M1NVz+LSBFSqlpJkuyjj+Ckk6B9e/jXv3T1s4gUScmhuli5Ek44AXbfPTRCN2oUd0QiksLiukJaqlJ2NvzqV/Ddd/Dqq9CyZdwRiUiKU3LIdLm5oYvP+fNhxgzo0iXuiEQkDSg5ZLrLLw/VSPfeC7/8ZdzRiEiaUJtDJrv33tDv829/G/pnEBEpJSWHTDVzJlx6aWhruO22uKMRkTSj5JCJFiyAU06Bn/0MnngCataMOyIRSTNKDpnm229D957164cOezLkzpMiUrXUIJ1BamzfDieeCKtWwSuv6JRVESk3JYdM4U7b22+HN94IVz/36BF3RCKSxlStlCluvZWfzpoVbqQ3dGjc0YhImlNyyATPPANXX803xxwTbsEtIlJBSg7p7qOP4PTToWtXPvn973WXVRGpFEoO6WzNmtCDW8OGMH06O+vVizsiEckQSg7pavv20LawciVMmwbNm8cdkYhkEJ2tlK4uuwxefhkeewwOOyzuaEQkw+jIIR3dfz+MGwdXXglnnBF3NCKSgZQc0s3rr8Mll0D//nDTTXFHIyIZSskhnSxfHrr5bN1a90wSkaRSm0O62LoVhgyBTZvgpZegceO4IxKRDKbkkA7c4fzz4Z13wplJHTrEHZGIZDhVK6WDu++GRx6BG26AQYPijkZEqgElh1T38suhq89Bg+CPf4w7GhGpJpQcUtmKFXDyyXDggfDoo1BDm0tEqobaHFLVtm3hCujNmyErC/bYI+6IRKQaUXJIVb/5DcydC1OmQPv2cUcjItVMLPUUZnabmX1sZh+Y2VQza5ww7moz+9zMPjGz4+OIL3bjx/9wBfRJJ8UdjYhUQ3FVYs8COrp7J+BT4GoAM+sAnAocAvQD7jOz6nWl17x5cOGF0Lcv3Hhj3NGISDUVS3Jw9xfcPTd6ORdoET0fBExy923u/gXwOVB97iq3Zk04UmjWDJ58Emqp1k9E4mHuHm8AZs8Ak939MTO7B5jr7o9F4x4C/uvuUwr53GhgNECzZs26TZo0qdwxZGdn06BBg3J/vlLs2EGnq6+m8fz5vH/33Wxs167Ms0iJclQSlSX1ZEo5QGXJ06dPn3fdvXuhI909KQMwG/iokGFQwjTXAlP5IUndA5yZMP4hYGhJy+rWrZtXxJw5cyr0+Upxww3u4H7//eWeRUqUo5KoLKknU8rhrrLkAeZ5Ef+rSau3cPe+xY03sxHAAODYKEiAFUDLhMlaRO9ltpkz4U9/grPOgnPPjTsaEZHYzlbqB4wBBrr75oRRM4BTzayume0PHAS8HUeMVWbZstAnQ8eO8I9/qA9oEUkJcbV43gPUBWZZ+DOc6+7nu/sCM3sKWAjkAhe5+46YYky+bdvg17+GnJxwPUP9+nFHJCICxJQc3P3AYsbdBFSPXmwuvxzefhuefhrato07GhGRXXSznrg8+STcdx/87nehnwYRkRSi5BCHTz6B0aOhZ0/461/jjkZEpAAlh6q2eXNoZ6hbFyZPhtq1445IRKQAXYJb1S69FD78EJ57Dlq0KHl6EZEY6MihKk2cCA89BNdcA/37xx2NiEiRlByqysKFoR/o//u/cMGbiEgKU3KoCps2hXaG3XeHJ57QDfVEJOXpX6oqXHopLFoEzz8P++4bdzQiIiXSkUOyPfFE6Lzn2mvhF7+IOxoRkVJRckimzz6D886Do46C66+POxoRkVJTckiWbdvg1FOhTh21M4hI2tE/VrKMGQPvvQfTp0PLliVPLyKSQnTkkAzTp8PYsXDZZTBwYNzRiIiUmZJDZfvqKxg5Erp2hVtuiTsaEZFyUXKoTDt2hI57cnLCfZPq1o07IhGRclGbQ2W66SZ49dVwm4wDi+yyQkQk5enIobK89lq4LcaZZ4ZBRCSNKTlUhnXrQnVS69Zw771xRyMiUmGqVqoo93Ch28qV8PrrsMcecUckIlJhSg4VNX48/OtfoUe3ww6LOxoRkUqhaqWK+PjjcFO9Y44JF72JiGQIJYfy2r4dTj8ddtstnJ1UQ6tSRDKHqpXK649/hPffh6lTdRtuEck42t0tj6wsuO02OPdcGDw47mhERCqdkkNZrVsHw4aFi9zuvDPuaEREkkLVSmXhHvqBXrUK3ngjdPspIpKBYjlyMLPbzOxjM/vAzKaaWePo/dZmtsXM5kfDuDjiK9LEifDUU+FK6B494o5GRCRp4qpWmgV0dPdOwKfA1QnjFrt752g4P57wCvHFF3DxxdCrF1x5ZdzRiIgkVSzJwd1fcPfc6OVcoEUccZRabm64X5JZOHqoWTPuiEREkioVGqRHAf9NeL2/mb1vZi+bWa+4gvqRW28NbQz33gutWsUdjYhI0pm7J2fGZrOBnxYy6lp3nx5Ncy3QHRji7m5mdYEG7r7WzLoB04BD3H1DIfMfDYwGaNasWbdJkyaVO9bs7GwaNGhQ6LgGn31G1wsuYE2vXiy87rpw9JCiiitHulFZUk+mlANUljx9+vR51927FzrS3WMZgBHAm0D9YqbJArqXNK9u3bp5RcyZM6fwEZs3u3fo4L7PPu5r11ZoGVWhyHKkIZUl9WRKOdxVljzAPC/ifzWWU1nNrB8wBvg/d9+c8P7ewHfuvsPM2gAHAUviiBGAa66BhQth5kzYc8/YwhARqWpxXedwD1AXmGWhmmauhzOTjgb+bGY5wE7gfHf/LpYIX3wR7roLLroIjj8+lhBEROISS3Jw90L70HT3p4Gnqzicgr7/HkaMgLZtQ2O0iEg1oyukC3PxxfD11+EMpfr1445GRKTKpcKprKllyhR4/HH4wx/UeY+IVFtKDom++QYuuAC6dYNrr407GhGR2Cg55Mm7qd7GjfDII1C7dtwRiYjERm0OeR57DKZNC/00HHJI3NGIiMRKRw5A3dWr4ZJL4Mgj4be/jTscEZHYKTm40+7WWyEnByZM0E31RESo5slh2Xq4a9T97DlvHmOH3sayvQu9/EJEpNqp1snh2gcWc+4TV7CsU3fu7HoBo2bEHZGISGqo1g3SX22AefsfyZLzzmZnjrFkXdwRiYikhmp95FDroAMYceHzZO/1E2oYtGkSd0QiIqmhWieH8QPhgCZghMfxA+OOSEQkNVTr5LBfI5g9DDr+JDzu1yjuiEREUkO1Tg4iIlI4JQcRESlAyUFERApQchARkQKUHEREpAAlBxERKcDcPe4YKszMVgNLKzCLvYA1lRROnDKlHKCypKJMKQeoLHlaufvehY3IiORQUWY2z927xx1HRWVKOUBlSUWZUg5QWUpD1UoiIlKAkoOIiBSg5BD8M+4AKkmmlANUllSUKeUAlaVEanMQEZECdOQgIiIFKDmIiEgB1TI5mNltZvaxmX1gZlPNrHER0/Uzs0/M7HMzu6qKwyyRmf3azBaY2U4zK/JUNjP70sw+NLP5ZjavKmMsrTKUJaW3CYCZ7Wlms8zss+ix0G6kzGxHtE3mm1nKdFJb0jo2s7pmNjka/5aZtY4hzFIpRVlGmNnqhO1wThxxlsTMxpvZt2b2URHjzczGRuX8wMy6Vnih7l7tBuA4oFb0/BbglkKmqQksBtoAdYD/AR3ijj1fjO2BdkAW0L2Y6b4E9oo73oqWJR22SRTnrcBV0fOrCvt+ReOy4461POsYuBAYFz0/FZgcd9wVKMsI4J64Yy1FWY4GugIfFTH+BOC/hL7LjgDequgyq+WRg7u/4O650cu5QItCJjsM+Nzdl7j7dmASMKiqYiwNd1/k7p/EHUdlKGVZUn6bRAYBj0TPHwEGxxdKmZVmHSeWbwpwrJlZFcZYWunyfSmRu78CfFfMJIOARz2YCzQ2s30qssxqmRzyGUXIuPk1B75KeL08ei8dOfCCmb1rZqPjDqYC0mWbNHP3r6Pnq4BmRUxXz8zmmdlcMxtcNaGVqDTreNc00U7WeqBplURXNqX9vpwUVcVMMbOWVRNapav030atCoWTwsxsNvDTQkZd6+7To2muBXKBx6sytrIoTTlK4Sh3X2FmPwFmmdnH0Z5IlaqksqSE4sqS+MLd3cyKOl+8VbRd2gAvmdmH7r64smOVYj0DPOnu28zsPMIR0TExx5QSMjY5uHvf4sab2QhgAHCsR5V2+awAEvciWkTvVamSylHKeayIHr81s6mEw+0qTw6VUJaU2CZQfFnM7Bsz28fdv44O7b8tYh5522WJmWUBXQh15HEqzTrOm2a5mdUCGgFrqya8MimxLO6eGPeDhPaidFTpv41qWa1kZv2AMcBAd99cxGTvAAeZ2f5mVofQ8JYyZ5SUlpntbmYN854TGuMLPeMhDaTLNpkBDI+eDwcKHBWZWRMzqxs93ws4ElhYZREWrTTrOLF8Q4GXitjBiluJZclXLz8QWFSF8VWmGcBZ0VlLRwDrE6o2yyfuVvg4BuBzQv3c/GjIO/NiX+C5hOlOAD4l7M1dG3fchZTjRELd4jbgG+D5/OUgnKnxv2hYkIrlKG1Z0mGbRDE2BV4EPgNmA3tG73cHHoye9wQ+jLbLh8DZccdd3DoG/kzYmQKoB/wr+h29DbSJO+YKlOWv0e/if8Ac4OC4Yy6iHE8CXwM50e/kbOB84PxovAH3RuX8kGLOXiztoNtniIhIAdWyWklERIqn5CAiIgUoOYiISAFKDiIiUoCSg4iIFKDkIJIkZtbYzC6MOw6R8lByEEmexoQ7mIqkHSUHkeS5GTgg6ifgtriDESkLXQQnkiRRJzjPunvHuGMRKSsdOYiISAFKDiIiUoCSg0jybAQaxh2ESHkoOYgkiYe+Al43s4/UIC3pRg3SIiJSgI4cRESkACUHEREpQMlBREQKUHIQEZEClBxERKQAJQcRESlAyUFERAr4/950EZ527K43AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "# Fit the model to a dense 1D-grid of points from [t_min, t_max].\n",
    "fitn = m*20\n",
    "fitx = np.linspace(np.min(t), np.max(t), fitn)  # Actually t.\n",
    "fity = vp2x(fitx)\n",
    "\n",
    "# Plot the fitted model against observations.\n",
    "plt.scatter(t, y, s=16, c='dodgerblue', label='observations')\n",
    "plt.plot(fitx, fity, c='red', label='fit')\n",
    "plt.title('Observations vs. Fit, Lagrange Interpolation')\n",
    "plt.ylabel('y')\n",
    "plt.xlabel('t')\n",
    "plt.grid()\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use Newton interpolation to interpolate 3 ($m = 3$) data points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "x:  [-27.  13.  -4.]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import scipy.linalg as la\n",
    "\n",
    "# Input data $(t_i, y_i)$.\n",
    "m = 3\n",
    "t = np.array([-2,0,1], dtype='d')\n",
    "y = np.array([-27,-1,0], dtype='d')\n",
    "\n",
    "\n",
    "def newton_basis(t):\n",
    "    \"\"\"\n",
    "    Return the lower triangular basis matrix used by Newton interpolation.\n",
    "    \"\"\"\n",
    "    A = np.zeros((t.size, t.size))\n",
    "    for i in range(t.size):\n",
    "        aij = 1.\n",
    "        for j in range(i+1):\n",
    "            A[i,j] = aij\n",
    "            # a_{ij} = \\pi_{j=1}^n (t_i - t_j)\n",
    "            aij = aij * (t[i]-t[j])\n",
    "    return A\n",
    "\n",
    "\n",
    "# Solve Lx = y for x using forward substitution.\n",
    "A = newton_basis(t)\n",
    "x = la.solve_triangular(A, y, lower=True, unit_diagonal=False)\n",
    "print(\"x: \", x)\n",
    "\n",
    "# Implement the polynomial p_2(t_i).\n",
    "p2ti = lambda ti: x[0] + x[1]*(ti-t[0]) + x[2]*(ti-t[0])*(ti-t[1])\n",
    "\n",
    "# Evaluate the polynomial at the data points.\n",
    "vp2ti = np.vectorize(p2ti)\n",
    "yhat = vp2ti(t)\n",
    "np.testing.assert_almost_equal(y - yhat, np.zeros(m))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the observations against the fit using Newton interpolation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAEWCAYAAACNJFuYAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAw50lEQVR4nO3debxW4/7/8denQUmphI1KKaJ00ohT8hUd4nQq6WtOiTIPx3FMOeTga/wZMuUgkaGcnJThRNEmQ0gyVKYMqYSSajfuvfv8/rjWzm0P7Xmv+773+/l4rMc9rOH+XPe67/VZ67rWWpe5OyIiIolqxB2AiIgkHyUHEREpQMlBREQKUHIQEZEClBxERKQAJQcRESlAySHJmNkoM3si7jhKy8z+a2ZD4o6jMphZTzP7PO44pOTMLNPMzizjvHuaWZaZ1azouFKJkkMVM7OhZvaJma03s+Vm9oCZNYo7rtIoLIG5+9Hu/lhcMVUEM2tpZh5tGPKGj9x9lrvvmzDdt2bWuxTLHRUt9/iE92pF77WswPiHmtmbFbW8QpZf4nKXZ+Nc1fKXy90Xu3t9d8+NM664KTlUITP7G3AL8HegIXAw0AKYbmbbVWEctarqs1JUo2jjUN/dD6igZf4CXFfd90ZLSt9TEnB3DVUwADsCWcDx+d6vD/wMDItejwImAROBtcBc4ICE6S8HlkbjPgeOiN6vAVwBLAJWAs8AO0XjWgIOnAEsBt4A/gucny+Wj4CB0fO7ge+BNcAHQM/o/T7AZiA7Ks9H0fuZwJkJsVwNfAf8BDwONMwXy5AolhXAyIQYDgTmRJ/7I3BHEd/nQqBvwuta0ffYGagLPBF9D78C7wMZJVhHebHVyvf+YcCS6Pl4YAuwISr/ZSVY7ijgyej7HZIQrwMto9d1gNuj7+RHYAywfTTudeC46HmPaL4/R6+PAOYBbYGNQG4U16/R+IbR9/9ztD6uBmpE44YCb0afuwr4Bjh6G+X4Fuhd3LzAjVEcG6NY7o3e3w+YTkiUn5PwXwDGAQ8ALwHrgN7Re2OiedZG30OLhHm6R+t2dfTYPWFcJr/9HlsDr0W/hxXRumhU1PrM/zsA9gCmRnF/BQzPt26fib7jtcB8oGvc25sK2WbFHUB1GQgb1RzybXiicY8BT0fPRxE2vIOA2sCl0R+vNrAvYYO9RzRtS6B19PwiYDbQjLCheTBhmXk/9seBHYDtgdOAtxJiaEfYkNaJXp8KNCFsxP4GLAfqJsT4RL4yJP4Zh0V/olaE5PcfYHy+WB6K4jgA2AS0jca/AwyOntcHDi7i+7wGeDLh9Z+BhdHzs4DngXpATaALsGMJ1tHvNgoJ7x9GlByi198SbSRLuO5HEZJVP+DraF3mTw53EjZAOwENovhvisb9E7gnen4VYQfgloRxd0fPhwJv5vvsx4Ep0TJbAl8AZyRMnw0Mj76nc4BlgBVRjq3lLm7exN9D9HoHwm/39KjsnQgb6nbR+HGEjXwPws5F3ei9tcChhN/03Xnli76nVcDgaHknRa+bFPJ73Bv4U7SMXQg7R3cVtT7z/w6i6e+PYupISLSHJ6zbjcAx0fdwEzA77u1NhWyz4g6gugyEje3yIsbdDEyPno9K/HFFf5QfgJ7Rj/wnwl5V7XzLWEh0FBG93j3689ZK+LG3ShjfgLCH1iJ6fSMwdhvxryI6gqH45PAqcG7CuH0LiaVZwvj3gBOj528A1wE7F/N97h1tOOpFr58EromeDwPeBjqUch3lxfZrwnApFZQcoufvEjakW5MDYNG6aJ0wzx+Bb6LnRwAfR8+nAWfm/UYIe9N5R3tDSUgOhI3VZqINcPTeWUBmwvRfJYyrF8W0WxHl2Fru4ualYHI4AZiVb3kPAtdGz8cBj+cbPw6YkPC6PuGIpDkhKbyXb/p3gKGFfX6+6QYAHxa1PhN+B7Wiz8oFGiSMvwkYl7BuZySMawdsKM3vLlkHtTlUnRXAzkXU9+8ejc/zfd4Td98CLCEcLXwFXEz4Qf5kZhPMbI9o0hbAZDP71cx+JSSLXCCjiOWuBV4ETozeOomwgQXAzC41s4VmtjpaXkNg5xKWdQ9CFUae7wh/tMRYlic8X0/440Oo+moDfGZm75tZ38I+IPouFgJ/MbN6hL3yp6LR44GXgQlmtszMbjWz2iWMHUJiahQNt5divpK4GhhJ2AvNswth4/pBwvqbFr0PYaPXxswyCHuujwPNzWxnQjXcG0WVg3CUkn9dNE14vXU9uPv66Gl9SqY087YADsorX1TGU4DdEqb5vpD5En+zWYSqnT0o+BuDgmUDwMwyov/KUjNbQziKK81v+Zfo/1LU5+T/LddNh3Y9JYeq8w6h+mRg4ptmVh84mrC3nad5wvgahKqiZQDu/pS7H0L4szmhgRvCn+johI1aI3ev6+5LE5br+WJ6GjjJzP5I2FjNjD6zJ6Hu9Xigsbs3IhzyWxHLyW9ZFF+ePQlVaj8WMx/u/qW7nwTsGpVtkpntUMTkTxOSWn9gQZQwcPdsd7/O3dsR6qX7EqrRKkpx5S96RvfphCq3cxPeXkGo894/Yd01dPf60TzrCe0+FwGfuvtmwpHRJcAid8/bscgf1wrCEVv+dbGUypc/lu+B1/P9Puu7+znbmAd+/1+oT6hOWkbB3xgUXbb/i5b9B3ffkXAUbwnjt7U+lwE7mVmDEnxOWlFyqCLuvppQXXKPmfUxs9rRaYzPEI4MxidM3sXMBkZ7HxcTkspsM9vXzA43szqEes4NhMY0CA13N5pZCwAz28XM+hcT1kuEP9g/gYnRUQqEKqccQt1qLTO7htCgnudHoGWUuArzNPBXM9sr+kP/X7T8nGLiwcxONbNdolh+jd7eUsTkE4AjCdU0eUcNmFkvM/tDdMbLGsIGsqhllMWPhPaUxLi/NbOhJZx/JCH5AluPDh8C7jSzXaPlNTWzoxLmeR04P3qEUG2S+DovrmZ5Z755OBXzGcLvokH027iEsOdc2fJ/Ry8Qjn4GR7/92mbWzczaFrOcY8zskKhM1xOq074n/HbbmNnJ0WnBJxCqdF4oZBkNCI3Nq82sKeFswW3FulX0WW8DN5lZXTPrQDi6TblrkUpLyaEKufuthAbF2wkbrXcJe1RHuPumhEmnEOpo8xrcBrp7NqFB7WbCHuFywt71ldE8dxMaNF8xs7WExumDiolnE6GxuDcJG1dClcw0QuPld4RElHjI/+/ocaWZzS1k0WMJye4NQmP6RuCCbcWSoA8w38yyojKd6O4bioj/B8IRWXfC2V15diOc8bWGUPX0ehQPZjbGzMaUMJai3ARcHVWPXBptuJoQvvNiuftbhHaWRJcTjihmR1UfMwhtNXleJ2zk3ijiNYQzcuYDy80s72jiAkJ7xteEs4ueIqyfynY3MMjMVpnZ6Kha5khCNeYywu/3FsJvelueAq4lVCd1Iez14+4rCUeEfyOchXQZ4ey1FYUs4zrCWWyrCVWp/8k3/nfrs5D5TyK0QywDJhPaSWYUE3fKyzuzQETKyMwOAc6LqsOkgpjZOMKJAFfHHUt1lPKNJiJxc/c3CXvlImlD1UoiIlKAqpVERKQAHTmIiEgBadHmsPPOO3vLli3LPP+6devYYYeiTqVPHelSDlBZklG6lANUljwffPDBCnffpbBxaZEcWrZsyZw5c8o8f2ZmJocddljFBRSTdCkHqCzJKF3KASpLHjPLf5X5VqpWEhGRApQcRESkACUHEREpIC3aHAqTnZ3NkiVL2LhxY7HTNmzYkIULF1ZBVJUrWctRt25dmjVrRu3apbkxqojEKW2Tw5IlS2jQoAEtW7bEzLY57dq1a2nQoME2p0kFyVgOd2flypUsWbKEvfbaK+5wRKSE0rZaaePGjTRp0qTYxCCVy8xo0qRJiY7gRCR5pG1yAJQYkoTWg0jqSdtqJRGRQm3eTO1ff4Uvv4TVq+HXX8Pjxo2wYUN4zBtyc8H9twGgVi3YbjuoU+e3x3r1YMcdfxsaNIBGjaBxY6iRmvvgSg6VbPTo0TzwwAMsX76cyy+/nCuuuILnnnuONm3a0K5du7jDE0kfOTmwbBksXgzffRcelyyBH3/8/bBmDT2qKqaaNWHnnWGXXcKQkQHNmkHz5r897rkn7LorJNkRdtImBzPrQ+gwpCbwsLvfHHNIZXL//fczY8YMmjVrtvW95557jr59+yo5iJSWe9jof/YZfPFF2Pv/4oswLF4c9vQTNW4cNsgZGdCpU3jcZRe+XLGCfbp2DXv3jRqFvf169aBu3TBsv304IqhZM2y0EzfcubmweXMYNm0Kw/r1sHYtrFnz27BqFfz8M/z0U3j8+Wd4/32YPDnMk2iHHWDvvaF16/C4997Qrh3sv3+ILwZJmRyi7h3vA/5E6ELzfTOb6u4L4o2sdM4++2y+/vprjj76aIYNG8aiRYs4+eSTmTp1Kq+//jo33HADzz77LK1bt447VJHks3YtfPghfPwxfPJJGD79NLyfp0EDaNMGDj4YTj4ZWrQIe+ItWoS98iLuObQ0M5N9ynr7jJo1Q/LYfvuyze8OK1aEo5rvvw9JbdEi+OorWLAAXnghJJ48e+wRkkT79nDAAdClC+y3H4vX1WLYVBhQB24YD2P7wZ4NyxZSYZIyOQAHAl+5+9cAZjaBqBP5Mi3t4oth3rwiR2+fmxtWeGl07Ah33bXNScaMGcO0adOYOXMmL7wQurbt3r07/fr1o2/fvgwaNKh0nymSrjZvDong/ffDMGcOLFz4Wz1/48bwhz/AaaeFx3btYJ99wpFAklXHFMvst2qmTp0Kjs/NDUljwQKYPz8kxPnzYcyY0CYCsP32ZO1xAEN270zjg3dnUfPDGDYVZgyuuDCTNTk05fd9Fi8hX3/IZjYCGAGQkZFBZmbm7xbQsGFD1kZ7GHU2b6ZG/sPNRO7kbGt8IbZs3symxD2YIhftZGVlsXHjRjZv3szatWvJzs5mw4YNW+OrKLm5uRW+zIqycePGAutoW7Kysko1fTJLl7JUZDlqZWWx4/z5NPzkExp+8gkNPvuMmtHe8ubGjVmz336sHTKErDZtWLv33mzeeeeCVTuffRaGMkiJdVKvHnTrFgaA3FzqLVlCgy++oP4XX9Bw/hec8ME4ltY+mEu6HYIBFVmkZE0OxXL3fwH/AujatavnvyvhwoULf7sg7P77t7mssl48tl0JpjEz6tevT926ddluu+1o0KABtWvXZvvtt6/wC9aS8SK4PHXr1qVTYXtJRdBdM5NPucqxfj28+SbMmAGvvhqOEtzDmT+dOsF550GPHnDQQWzXtCk7m7FzhUb/e+mwTnqPh69XbuGyxtO4Y/lhtG4M5x9WcctP1uSwFGie8LpZ9F5aaNCgQdLu4YtUCPfQRvDii/Dyy/DOO6HqqHZt+OMf4dpr4dBD4cADi2wXkG0b2w+GTa1BTp16tK4RXlekZE0O7wP7mNlehKRwInByvCFVnBNPPJHhw4czevRoJk2apAZpSQ/r1oUjg5deCsOSJeH9jh3hwgvhiCOgZ08lgwqyZ8PQxpCZWbFHDHmSMjm4e46ZnQ+8TDiVday7z485rDL59ttvARg6dChDhw4FoEePHixYkFInXokUbtUqeP55+M9/whHCxo3hDKIjj4TrroOjj4bdd487SimDpEwOAO7+EvBS3HGIVHeLV/O7UybHHbKSZq9MgmefhZkzw8VnzZrB8OEwYAAccki4clhSWtImBxFJDsOmwvIf1tLm++lc9dpt7Pb5K5CbE04lvfRSGDgQunZNvVNKZZuUHESkcLm5MH06F905jt6fTmX77A0sbdScRw67hLNuPylckKWEkLaUHETk9z7/HMaNg8cfh2XLOLR+EyYddDrb/2k/LtvxPFo1qcFZHeMOUiqbkoOIhCtvn3kGHnwwnHZas2ZoTB49mtU9+/LYy3UYUCeTVptqVPgpk5KclBxEqrOFC0NCeOyxcOvqffeFW2+FU0/depbRnlTuKZOSnFLzRuMp7Ntvv6V9+/Zxh8G8efN46aXfTgabOnUqN9+ckje+ldLKyQlnGvXqFe5RdP/90KdP2PovXAh//7tOPxUdOaSDnJwcatUq3aqcN28ec+bM4ZhjjgGgX79+9Oun+oK0tmoVPPww3HtvuBNoy5Zw881w+umhPwGRBDpyqGR33HEH7du3p3379twV3cU1JyeHU045hbZt2zJo0CDWr18PwBVXXEG7du3o0KEDl156KQA///wzxx13HN26daNbt2689dZbAIwaNYrBgwfTo0cPBg8ezMEHH8zChQu3fu5hhx3GnDlzeO+99/jjH/9Ip06d6N69O59//jmbN2/mmmuuYeLEiXTs2JGJEycybtw4zj//fCAc3Rx++OF06NCBI444gsWLFwPhQr4LL7yQ7t2706pVKyZNmgTADz/8wKGHHkrHjh1p3749s2bNqpLvVkroiy/g3HPDtQiXXQatWoU+Bb76Ci6/XIlBCufuKT906dLF81uwYEGB94qyZs0a/+5X9yMed9/r7vD43a8lnr1Ic+bM8fbt23tWVpavXbvW27Vr53PnznXA33zzTXd3P/300/22227zFStWeJs2bXzLli3u7r5q1Sp3dz/ppJN81qxZ7u7+3Xff+X777efu7tdee6137tzZ169f7+7ud9xxh19++eXu7r5s2TJv06aNu7uvXr3as7Oz3d19+vTpPnDgQHd3f/TRR/28887bGmvi6759+/q4cePc3f2RRx7x/v37u7v7kCFDfNCgQZ6bm+vz58/31q1bu7v77bff7jfccIO7u+fk5PiaNWsKfBelWR/u7jNnzizV9MkstrK8+677wIHuZu7bbed++unuH35Y5sVpnSSn8pQFmONFbFdVrRQZNhUWrYItHh4r4t7ob775Jsceeyw7RPeSGThwILNmzaJ58+b06BE6Kjz11FMZPXo0F198MXXr1uWMM86gb9++9O3bF4AZM2b87lYba9asISsrCwhVQdtHHY4cf/zx9O7dm5tvvplnnnlma18Rq1evZsiQIXz55ZeYGdnZ2cXG/c477/Cf//wHgMGDB3PZZZdtHTdgwABq1KhBu3bt+PHHHwHo1q0bw4YNIzs7mwEDBtCxY8fyfG1SHu4wbVpoVM7MDL2IXXlluLdRRkbc0UkKUbVS5OsoMUB4/HpV5X2W5btwyMyoVasW7733HoMGDeKFF16gT58+IZYtW5g9ezbz5s1j3rx5LF26lPr16wNsTToATZs2ZaedduLjjz9m4sSJnHDCCQD84x//oFevXnz66ac8//zzbNy4sVyx16lTZ+tzjzpiOfTQQ3njjTdo2rQpQ4cO5fHHHy/XZ0gZbNkSqoq6dIFjjgndZ/6//xfaFm68UYlBSk3JIdKqMdSIttk1LLwur549e/Lcc8+xfv161q1bx+TJk+nZsyeLFy/mnXfeAeCpp57ikEMOISsri9WrV3PMMcdw55138tFHHwFw5JFHcs8992xd5rxt9Gg3cOBAbr31VlavXk2HDh2AcOTQtGlTAMaNG7d12m3dNrx79+5MmDABgCeffJKePXtus5zfffcdGRkZDB8+nDPPPJO5c+du+4uRirNlC/z73+HOpwMHQlYWjB0LX38Nl1wSboInUgZKDpGx/aB1Y6hp4bEiLvTp3LkzQ4cO5cADD+Sggw7izDPPpHHjxuy7777cd999tG3bllWrVnHOOeewdu1a+vbtS4cOHTjkkEO44447ABg9ejRz5syhQ4cOtGvXjjFjxhT5eQMGDGDChAkcf/zxW9+77LLLuPLKK+nUqRM5OTlb3+/VqxcLFizY2iCd6J577uHRRx+lQ4cOjB8/nrvvvnub5czMzOSAAw6gU6dOTJw4kYsuuqgsX5eUxpYtMGFC6DLz+OMhOxueeCJ0LXn66brxnZRfUY0RqTRURIN0OkjmcqhBuoJs2eI+ZYp7hw7u4L7//u4TJrjn5FTcZxRB6yQ5VVaDtI4cRFLFq6+GXtT69w/dbj71FHz8MZxwQrjdhUgFUnIQSXYffBB6UevdG5Ytg4ceCtVHJ50ENfQXlsqR1r8sj86mkXhpPZTRt9/CKaeEvhI+/hjuuitc0HbmmaEvZpFKlLbJoW7duqxcuVIbppi5OytXrqRu3bpxh5I6Vq0K9zfad9/Q/eaVV4armS+6CPQ9ShVJ24vgmjVrxpIlS/j555+LnXbjxo1psfFK1nLUrVuXZs2axR1G8svJCXdIveaakCCGDIHrrw+3vRCpYmmbHGrXrs1ee+1VomkzMzPp1KlTJUdU+dKlHNXSq6+GI4P58+Hww+GOO0JPayIxSdtqJZGUsGgRHHtsaGxevz5c5TxjhhKDxE7JQSQOGzbAP/4R+lOYPh1uuimcgTRggPpllqSQttVKIknrhRfCjfC++Sb0uHbLLbDHHnFHJfI7OnIQqSrffReODP7yl3DW0cyZMH68EoMkJSUHkUpmubnhFtpt24YqpFtugXnz4LDD4g5NpEiqVhKpTHPm0Pnss8N1Cv37w+jRsOeecUclUiwdOYhUhnXrwi2zDzqI7VatgmefheeeU2KQlJF0ycHMRpnZUjObFw3HxB2TSKm88gq0bw933gnDh/P+uHGhrwWRFJJ0ySFyp7t3jIaX4g5GpERWr4bhw+Goo6BOHXjjDRgzhpyo5z6RVJKsyUEktUybFo4Wxo6Fyy8PDc7F9KAnksws2W5MZ2ajgKHAGmAO8Dd3L9Cjs5mNAEYAZGRkdMnr1rIssrKytvbLnMrSpRyQOmWplZVF6/vvZ/f//pd1LVrw2eWXs7Zt299NkyplKU66lANUljy9evX6wN27FjqyqF6AKnMAZgCfFjL0BzKAmoSjmhuBscUtr7Ce4EojXXqFSpdyuKdIWV591b15c/caNdyvuMJ9w4ZCJ0uJspRAupTDXWXJwzZ6govlVFZ3712S6czsIeCFSg5HpHQ2bICrrgr9K7RpA++8AwceGHdUIhUq6doczGz3hJfHEo4oRJLD3LnQpUtIDOefDx9+qMQgaSkZL4K71cw6Ag58C5wVazQiALm54crma6+FXXYJDdBHHRV3VCKVJumSg7sPjjsGkd9ZsiTcIO/11+H44+GBB2CnneKOSqRSJV1yEEkqkyfDGWfA5s0wbhycdppuqS3VQtK1OYgkhfXr4ZxzwpXNrVqFtoUhQ5QYpNpQchDJb8GC0Mg8Zgz8/e/w9tuwzz5xRyVSpVStJJLoscfg3HOhfn14+WU48si4IxKJhY4cRCDcRfX002HoUDjooHD7CyUGqcaUHETyqpEeewyuuSZ0yLP77sXPJ5LGVK0k1duTT8KIEaEa6ZVXoHeJLt4XSXs6cpDqafNmuOCCcP1Cly6hGkmJQWQrJQepfpYuDf0333tv6K3t1VdVjSSSj6qVpHrJzIQTTggN0BMnhiueRaQAHTlI9eAeuu3s3RsaN4b33lNiENkGJQdJfxs2hKubL7kE/vKXkBjatYs7KpGkpuQg6e377+HQQ2H8eLjuOnj2Wdhxx7ijEkl6anOQ9PXmm3DcceE+Sc89B/37xx2RSMrQkYOkp4cegsMPD0cJ776rxCBSSkoOkl5yc+Gvfw0XtvXqpfYFkTJScpD0sWYN9OsXuvC88EJ48cVwZpKIlJraHCQ9fPNNOBPps89CT21nnx13RCIpTclBUt9bb8GAAZCTE26zfcQRcUckkvJUrSSpbeLE0PDcuHFoeFZiEKkQSg6Smtzh5pvhxBND/wvvvANt2sQdlUjaULWSpJ7s7NBb28MPw0knwaOPQp06cUclklZ05CCpZc0a6Ns3JIaRI+GJJ5QYRCqBjhwkdSxbBkcfDfPnh+RwxhlxRySStpQcJDUsXAh9+sAvv8BLL6l/Z5FKpuQgye/tt8M1DLVrw+uvQ+fOcUckkvbU5iDJbcqUcHpqkyYhSSgxiFSJWJKDmf2vmc03sy1m1jXfuCvN7Csz+9zMjoojPkkSDz4IAwfCH/4QLnRr1SruiESqjbiOHD4FBgJvJL5pZu2AE4H9gT7A/WZWs+rDk1i502L8+HALjD59YOZM2GWXuKMSqVZiSQ7uvtDdPy9kVH9ggrtvcvdvgK+AA6s2OonVli1wySXsNXYsnHpq6Idhhx3ijkqk2km2BummwOyE10ui9wowsxHACICMjAwyMzPL/KFZWVnlmj9ZpHo5LCeHfW+7jd1eeYWv+/Vj8emnh+qkFJfq6yVPupQDVJaSqLTkYGYzgN0KGTXS3aeUd/nu/i/gXwBdu3b1ww47rMzLyszMpDzzJ4uULseGDXDCCfDKK3D99Szu0YPDevWKO6oKkdLrJUG6lANUlpKotOTg7r3LMNtSoHnC62bRe5LO1qwJp6rOmgX33w/nnANpslcnkqqS7VTWqcCJZlbHzPYC9gHeizkmqUwrV4ZTVd9+G556KiQGEYldLG0OZnYscA+wC/Cimc1z96Pcfb6ZPQMsAHKA89w9N44YpQosXw5/+hN8+WVoeP7zn+OOSEQisSQHd58MTC5i3I3AjVUbkVS5xYuhd+9wv6SXXgp9MohI0ki2s5WkOvjqq1CVtHp1aIDu3j3uiEQkHyUHqVoLFoQjhs2b4bXXdDsMkSSl5CBV55NPwhFDzZrhBnr77x93RCJShGQ7W0nS1bx50KtXuLNqZqYSg0iSU3KQyjd3bmhwrlcvHDHsu2/cEYlIMZQcpHK9/36oStpxx5AY9t477ohEpASUHKTyzJ4dGp932ikkhr32ijsiESkhJQepHO++G7ry3HXXkBhatIg7IhEpBSUHqXjvv/9bYpg5E5o1izsiESklJQepWB98EBJDkyZKDCIpTMlBKs7cueFeSY0ahcTQvHmxs4hIclJykIoxb15ofG7QICQGtTGIpDQlBym/+fNDYqhfP1zg1rJl3BGJSDkpOUj5fPFFuI5hu+3CvZJ0uqpIWtC9laTsvvkmJIYtW0JVki5wE0kbSg5SNkuWhMSwbl1IDG3bxh2RiFSgYquVzOwCM2tcFcFIili+PCSGFSvg5ZfhgAPijkhEKlhJ2hwygPfN7Bkz62NmVtlBSRL75ZdwuuqSJaEHt27d4o5IRCpBscnB3a8G9gEeAYYCX5rZ/5lZ60qOTZLN2rVw9NGhEXrqVDjkkLgjEpFKUqKzldzdgeXRkAM0BiaZ2a2VGJskk40bYcCAcAX0xImhWklE0laxDdJmdhFwGrACeBj4u7tnm1kN4EvgssoNUWKXnQ0nnhhOVX388ZAkRCStleRspZ2Age7+XeKb7r7FzPpWTliSNLZsgWHDYMoUuOceGDw47ohEpAoUmxzc/dptjFtYseFIUnGHiy6CJ56AG26A88+POyIRqSK6QlqKdv31cO+98Le/wVVXxR2NiFQhJQcp3JgxcO21MGQI3HYb6AxmkWpFyUEKmjQJzj0X+vaFhx5SYhCphpQc5Pdeew1OOQW6dw+nrNauHXdEIhKDWJKDmf2vmc03sy1m1jXh/ZZmtsHM5kXDmDjiq7bmzoX+/aFNG3j+eahXL+6IRCQmcd1471NgIPBgIeMWuXvHqg1HWLQoXP3cpAlMmwaNdTstkeosluSQdwqsbtOUJH76CY46CnJzw430mjaNOyIRiZmFO2PE9OFmmcCl7j4net0SmA98AawBrnb3WUXMOwIYAZCRkdFlwoQJZY4jKyuL+vXrl3n+ZFGWctTYsIGOl1zCDt98w0d33MGadu0qKbrSSZd1AulTlnQpB6gseXr16vWBu3ctdKS7V8oAzCBUH+Uf+idMkwl0TXhdB2gSPe8CfA/sWNxndenSxctj5syZ5Zo/WZS6HNnZ7n/+s3uNGu5TplRKTGWVLuvEPX3Kki7lcFdZ8gBzvIjtaqVVK7l77zLMswnYFD3/wMwWAW2AORUcnrjDOefAiy+Gaxr69Ys7IhFJIkl1KquZ7WJmNaPnrQi3Cv863qjS1D//CQ8/DCNHwllnxR2NiCSZuE5lPdbMlgB/BF40s5ejUYcCH5vZPGAScLa7/xJHjGnt0Udh1Khw9fP118cdjYgkobjOVpoMTC7k/WeBZ6s+ompkxgwYMSL05qarn0WkCElVrSSV7NNP4bjjoG1b+Pe/dfWziBRJyaG6WLYMjjkGdtghNEI3bBh3RCKSxOK6QlqqUlYW/OUv8MsvMGsWNG8ed0QikuSUHNJdTk7o4nPePJg6FTp1ijsiEUkBSg7p7pJLQjXSfffBn/8cdzQikiLU5pDO7rsv9Pv817+G/hlEREpIySFdTZsGF14Y2hpuuy3uaEQkxSg5pKP58+GEE+APf4CnnoKaNeOOSERSjJJDuvnpp9C9Z716ocOeNLnzpIhULTVIp5EamzfDscfC8uXwxhs6ZVVEykzJIV240+b22+Htt8PVz926xR2RiKQwVSuli1tvZbfp08ON9AYNijsaEUlxSg7p4Pnn4cor+fHww8MtuEVEyknJIdV9+imcfDJ07sznf/+77rIqIhVCySGVrVgRenBr0ACmTGFL3bpxRyQiaULJIVVt3hzaFpYtg+eeg6ZN445IRNKIzlZKVRddBK+/Dk88AQceGHc0IpJmdOSQih58EMaMgcsvh1NOiTsaEUlDSg6p5q234IIL4Oij4cYb445GRNKUkkMqWbIkdPPZsqXumSQilUptDqli40YYOBDWrYPXXoNGjeKOSETSmJJDKnCHs8+G998PZya1axd3RCKS5lStlAruuQceewxGjYL+/eOORkSqASWHZPf666Grz/794R//iDsaEakmlByS2dKlcPzxsPfe8PjjUEOrS0SqhtocktWmTeEK6PXrITMTdtwx7ohEpBpRckhWF18Ms2fDpEnQtm3c0YhINRNLPYWZ3WZmn5nZx2Y22cwaJYy70sy+MrPPzeyoOOKL3dixv10BfdxxcUcjItVQXJXY04H27t4B+AK4EsDM2gEnAvsDfYD7zax6Xek1Zw6cey707g033BB3NCJSTcWSHNz9FXfPiV7OBppFz/sDE9x9k7t/A3wFVJ+7yq1YEY4UMjLg6aehlmr9RCQe5u7xBmD2PDDR3Z8ws3uB2e7+RDTuEeC/7j6pkPlGACMAMjIyukyYMKHMMWRlZVG/fv0yz18hcnPpcOWVNJo3jw/vuYe1++5b6kUkRTkqiMqSfNKlHKCy5OnVq9cH7t610JHuXikDMAP4tJChf8I0I4HJ/Jak7gVOTRj/CDCouM/q0qWLl8fMmTPLNX+FGDXKHdwffLDMi0iKclQQlSX5pEs53FWWPMAcL2K7Wmn1Fu7ee1vjzWwo0Bc4IgoSYCnQPGGyZtF76W3aNLjuOjjtNBg+PO5oRERiO1upD3AZ0M/d1yeMmgqcaGZ1zGwvYB/gvThirDKLF4c+Gdq3hwceUB/QIpIU4mrxvBeoA0y3sDGc7e5nu/t8M3sGWADkAOe5e25MMVa+TZvgf/8XsrPD9Qz16sUdkYgIEFNycPe9tzHuRqB69GJzySXw3nvw7LPQpk3c0YiIbKWb9cTl6afh/vvhb38L/TSIiCQRJYc4fP45jBgB3bvDTTfFHY2ISAFKDlVt/frQzlCnDkycCLVrxx2RiEgBugS3ql14IXzyCbz0EjRrVvz0IiIx0JFDVRo/Hh55BK66Co4+Ou5oRESKpORQVRYsCP1A/8//hAveRESSmJJDVVi3LrQz7LADPPWUbqgnIklPW6mqcOGFsHAhvPwy7LFH3NGIiBRLRw6V7amnQuc9I0fCn/4UdzQiIiWi5FCZvvwSzjoLDjkErr027mhEREpMyaGybNoEJ54I222ndgYRSTnaYlWWyy6DuXNhyhRo3rz46UVEkoiOHCrDlCkwejRcdBH06xd3NCIipabkUNG+/x5OPx06d4Zbbok7GhGRMlFyqEi5uaHjnuzscN+kOnXijkhEpEzU5lCRbrwRZs0Kt8nYu8guK0REkp6OHCrKm2+G22KcemoYRERSmJJDRVi1KlQntWwJ990XdzQiIuWmaqXycg8Xui1bBm+9BTvuGHdEIiLlpuRQXmPHwr//HXp0O/DAuKMREakQqlYqj88+CzfVO/zwcNGbiEiaUHIoq82b4eSTYfvtw9lJNfRVikj6ULVSWf3jH/DhhzB5sm7DLSJpR7u7ZZGZCbfdBsOHw4ABcUcjIlLhlBxKa9UqGDw4XOR2551xRyMiUilUrVQa7qEf6OXL4e23Q7efIiJpKJYjBzO7zcw+M7OPzWyymTWK3m9pZhvMbF40jIkjviKNHw/PPBOuhO7WLe5oREQqTVzVStOB9u7eAfgCuDJh3CJ37xgNZ8cTXiG++QbOPx969oTLL487GhGRShVLcnD3V9w9J3o5G2gWRxwllpMT7pdkFo4eataMOyIRkUqVDA3Sw4D/Jrzey8w+NLPXzaxnXEH9zq23hjaG++6DFi3ijkZEpNKZu1fOgs1mALsVMmqku0+JphkJdAUGurubWR2gvruvNLMuwHPA/u6+ppDljwBGAGRkZHSZMGFCmWPNysqifv36hY6r/+WXdD7nHFb07MmCa64JRw9JalvlSDUqS/JJl3KAypKnV69eH7h710JHunssAzAUeAeot41pMoGuxS2rS5cuXh4zZ84sfMT69e7t2rnvvrv7ypXl+oyqUGQ5UpDKknzSpRzuKkseYI4XsV2N5VRWM+sDXAb8j7uvT3h/F+AXd881s1bAPsDXccQIwFVXwYIFMG0a7LRTbGGIiFS1uK5zuBeoA0y3UE0z28OZSYcC/zSzbGALcLa7/xJLhK++CnfdBeedB0cdFUsIIiJxiSU5uHuhfWi6+7PAs1UcTkG//gpDh0KbNqExWkSkmtEV0oU5/3z44YdwhlK9enFHIyJS5ZLhVNbkMmkSPPkkXH21Ou8RkWpLySHRjz/COedAly4wcmTc0YiIxEbJIU/eTfXWroXHHoPateOOSEQkNmpzyPPEE/Dcc6Gfhv33jzsaEZFY6cgBqPPzz3DBBdCjB/z1r3GHIyISOyUHd/a99VbIzoZx43RTPRERqnlyWLwa7hr2IDvNmcPoQbexeJdCL78QEal2qnVyGPnQIoY/dSmLO3Tlzs7nMGxq3BGJiCSHat0g/f0amLNXD74+6wy2ZBtfr4o7IhGR5FCtjxxq7dOaoee+TNbOu1LDoFXjuCMSEUkO1To5jO0HrRuDER7H9os7IhGR5FCtk8OeDWHGYGi/a3jcs2HcEYmIJIdqnRxERKRwSg4iIlKAkoOIiBSg5CAiIgUoOYiISAFKDiIiUoC5e9wxlJuZ/Qx8V45F7AysqKBw4pQu5QCVJRmlSzlAZcnTwt13KWxEWiSH8jKzOe7eNe44yitdygEqSzJKl3KAylISqlYSEZEClBxERKQAJYfgX3EHUEHSpRygsiSjdCkHqCzFUpuDiIgUoCMHEREpQMlBREQKqJbJwcxuM7PPzOxjM5tsZo2KmK6PmX1uZl+Z2RVVHGaxzOx/zWy+mW0xsyJPZTOzb83sEzObZ2ZzqjLGkipFWZJ6nQCY2U5mNt3MvoweC+1Gysxyo3Uyz8ySppPa4r5jM6tjZhOj8e+aWcsYwiyREpRlqJn9nLAezowjzuKY2Vgz+8nMPi1ivJnZ6KicH5tZ53J/qLtXuwE4EqgVPb8FuKWQaWoCi4BWwHbAR0C7uGPPF2NbYF8gE+i6jem+BXaOO97yliUV1kkU563AFdHzKwr7fUXjsuKOtSzfMXAuMCZ6fiIwMe64y1GWocC9ccdagrIcCnQGPi1i/DHAfwl9lx0MvFvez6yWRw7u/oq750QvZwPNCpnsQOArd//a3TcDE4D+VRVjSbj7Qnf/PO44KkIJy5L06yTSH3gsev4YMCC+UEqtJN9xYvkmAUeYmVVhjCWVKr+XYrn7G8Av25ikP/C4B7OBRma2e3k+s1omh3yGETJufk2B7xNeL4neS0UOvGJmH5jZiLiDKYdUWScZ7v5D9Hw5kFHEdHXNbI6ZzTazAVUTWrFK8h1vnSbayVoNNKmS6EqnpL+X46KqmElm1rxqQqtwFf7fqFWucJKYmc0Aditk1Eh3nxJNMxLIAZ6sythKoyTlKIFD3H2pme0KTDezz6I9kSpVQWVJCtsqS+ILd3czK+p88RbRemkFvGZmn7j7ooqOVbbpeeBpd99kZmcRjogOjzmmpJC2ycHde29rvJkNBfoCR3hUaZfPUiBxL6JZ9F6VKq4cJVzG0ujxJzObTDjcrvLkUAFlSYp1Atsui5n9aGa7u/sP0aH9T0UsI2+9fG1mmUAnQh15nEryHedNs8TMagENgZVVE16pFFsWd0+M+2FCe1EqqvD/RrWsVjKzPsBlQD93X1/EZO8D+5jZXma2HaHhLWnOKCkpM9vBzBrkPSc0xhd6xkMKSJV1MhUYEj0fAhQ4KjKzxmZWJ3q+M9ADWFBlERatJN9xYvkGAa8VsYMVt2LLkq9evh+wsArjq0hTgdOis5YOBlYnVG2WTdyt8HEMwFeE+rl50ZB35sUewEsJ0x0DfEHYmxsZd9yFlONYQt3iJuBH4OX85SCcqfFRNMxPxnKUtCypsE6iGJsArwJfAjOAnaL3uwIPR8+7A59E6+UT4Iy4497Wdwz8k7AzBVAX+Hf0P3oPaBV3zOUoy03R/+IjYCawX9wxF1GOp4EfgOzof3IGcDZwdjTegPuicn7CNs5eLOmg22eIiEgB1bJaSUREtk3JQUREClByEBGRApQcRESkACUHEREpQMlBpJKYWSMzOzfuOETKQslBpPI0ItzBVCTlKDmIVJ6bgdZRPwG3xR2MSGnoIjiRShJ1gvOCu7ePOxaR0tKRg4iIFKDkICIiBSg5iFSetUCDuIMQKQslB5FK4qGvgLfM7FM1SEuqUYO0iIgUoCMHEREpQMlBREQKUHIQEZEClBxERKQAJQcRESlAyUFERApQchARkQL+P1v/fksrbs1gAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "# Fit the model to a dense 1D-grid of points from [t_min, t_max].\n",
    "fitn = m*20\n",
    "fitx = np.linspace(np.min(t), np.max(t), fitn)  # Actually t.\n",
    "fity = vp2ti(fitx)\n",
    "\n",
    "# Plot the fitted model against observations.\n",
    "plt.scatter(t, y, s=16, c='dodgerblue', label='observations')\n",
    "plt.plot(fitx, fity, c='red', label='fit')\n",
    "plt.title('Observations vs. Fit, Newton Interpolation')\n",
    "plt.ylabel('y')\n",
    "plt.xlabel('t')\n",
    "plt.grid()\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 07.04 Orthogonal Polynomials"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Inner Product\n",
    "\n",
    "Inner product of polynomials on interval $[a, b]$ given by:\n",
    "$$\n",
    "\\langle p, q \\rangle = \\int_a^b p(t) q(t) w(t) dt\n",
    "$$\n",
    "where\n",
    "* $w(t)$ is a **weight function**\n",
    "\n",
    "Polynomials $p$ and $q$ are orthogonal if $\\langle p, q \\rangle = 0$.\n",
    "\n",
    "Set of polynomials $\\{p_i\\}$ are orthonormal if they are mutally orthogonal and the magnitude of each polynomial is 1.\n",
    "* Similar to vectors, but object is a polynomial.\n",
    "\n",
    "If we apply Gram-Schmidt process with the weight function defined below to a set of polynomials $\\{p_i\\}$ then we obtain a family of orthogonal polynomials bearing the names listed below.\n",
    "\n",
    "| Name | Symbol | Interval | Weight Function, $f(t)$ |\n",
    "|------|--------|----------|-------------------------|\n",
    "| Legendre | $P_k$ | $[-1, 1]$ | 1 |\n",
    "| Chebyshev, 1st | $T_k$ | $[-1, 1]$ | $(1 - t^2)^{-1/2}$ |\n",
    "| Chebyshev, 2nd | $U_k$ | $[-1, 1]$ | $(1 - t^2)^{1/2}$ |\n",
    "| Jacobi | $J_k$ | $[-1, 1]$ | $(1 - t)^{\\infty} (1 + t)^{\\beta}$ |\n",
    "| Laguerre | $L_k$ | $[0, \\infty]$ | $e^{-t}$ |\n",
    "| Hermite | $H_k$ | $[-\\infty, \\infty]$ | $e^{-t^2}$ |\n",
    "\n",
    "#### Chebyshev Polynomials\n",
    "**Equi-oscillation property:** successive extrema of $T_k$ are equal in magnitude and opposite in sign.\n",
    "* Distributes error uniformly when approximating continuous functions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 07.05 Convergence"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Question: When you interpolate a continous function, $f(t)$, by a polynomial, $p_{n-1}(t)$, at a given set of sample points, how much can the polynomial differ from the continuous function between the data points?\n",
    "\n",
    "1. Tight Bounds, to be applied at some point $\\theta$ in $[t_1, t_n]$\n",
    "$$\n",
    "f(t) - p_{n-1}(t) = \\frac{f^{(n)}(\\theta)}{n!} (t-t_1)(t-t_2)\\cdots(t-t_n)\n",
    "$$\n",
    "where\n",
    "  * $f^{(n)}(\\theta)$ is the nth derivative of $f(t)$ at some point $\\theta$ in $[t_1, t_n]$\n",
    "  \n",
    "2. Loose Bounds, to be applied across the entire interval $[t_1, t_n]$\n",
    "$$\n",
    "\\max_{t \\in [t_1, t_n]} |f(t) - p_{n-1}(t)| \\leq \\frac{Mh^n}{4n}\n",
    "$$\n",
    "where\n",
    "  * $M$ is an upper bound on $f^{(n)}(\\theta)$\n",
    "  * $h$ is the maximum spacing between samples\n",
    "  * $n$ is the number of samples\n",
    "\n",
    "#### Expected Accuracy\n",
    "* Equally spaced sample points are expected to have inflated errors at the extrema of the interval $[t_1, t_n]$.\n",
    "* Sample points spaced at Chebyshev points distribute the error more evenly, resulting in smaller errors at the extrema.\n",
    "* Avoid high degree polynomial interpolants.\n",
    "  * A polynomial of degree $n$ has $n-1$ inflection points and these oscillations may bear no relation to the data to be fit. \n",
    "  * The cost to determine and evaluate a polynomial increases with $n$.\n",
    "\n",
    "\n",
    "#### Truncated Taylor Series\n",
    "For completeness, consider Taylor series approximation of $f(t)$ when $t$ near $a$ as a polynomial interpolant of exactly one sample $(a, y)$.\n",
    "$$\n",
    "p_n(t) = f(a) + f'(a)(t-a) + \\frac{f''(a)}{2}(t-a)^2 + \\cdots + \\frac{f^{(n)}(a)}{n!}(t-a)^n\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compare the Taylor approximation of $f(t) = ln(t)$ at $t = 1$ for Taylor polynomial of degree 0 to 4."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYAAAAEWCAYAAABv+EDhAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABSvElEQVR4nO3deXxU1d348c+ZyUz2fSUEEvZ93xREQREBEVyLdaut1Yra5+mjdrN9qvaxrW391WpbtVjXFsR9VxZlE1DZwxKWJEAgIQlZJ5lss53fH3cIISQQIMkkM9/36zWvzNz1e+5Mzvfec5ejtNYIIYQIPCZfByCEEMI3JAEIIUSAkgQghBABShKAEEIEKEkAQggRoCQBCCFEgJIE0EUppZKVUuuUUtVKqf/n63iaU0rdqZRa3wXiOKyUmuHrOJpSSn2ulPpeG6dt9/iVUlOUUtlKKbtS6tr2XHZnU0pNVUrt93Uc/koSQNd1D1AKRGmtH7qQBSmlHvFWBnalVL1Syt3k8572CffCeCvCOm9MxUqpV5VSEb6O63xorWdrrV/zYQi/Bf6utY7QWn/QfGRXTJqt0Vp/pbUe1BHLVkqt8f4/VCulqpRSW5VSv1BKBXfE+roiSQBdVzqQpc/jTj2lVFDTz1rr33srgwjgXuDrE5+11sPaKd62xqaUUq397q7xxjgWGA/8uvMi8yvpQKcm9ua/uW7kAa11JNADeAi4GfhMKaXae0VdcRsFTAJQSv1UKfVus2HPKqWeaWHanyul3mk27Bml1LPe93cqpQ569xwOKaVubcP645RS+Uqpa7yfI5RSOUqpO1qY9lXge8DPvHvEM5RSwUqpvyqljnlffz2xp6KUmuZd9s+VUkXAK+ewXZ5RSh1tsgc01Ts8RSlVq5SKbzLtWKVUiVLK0sJyJiulNiulbN6/k5uMW6OU+p1SagNQC/Q9U0xa6wLgc2C4d/55Sqk9SqlK77KGtLD+M8Z7oslKKfWUUqrC+73NbjJtqlLqI6VUufd7ubvJuMeUUm8rpf7j/c53KaUGKqV+qZQ67t1+M5uV94fe9/2UUquUUmVKqVKl1GKlVEwr38VEpdQW73dRrJT6S2vbSCl1tzfOcm/cqd7hud7t+7H3t3PGvdk2bJc4pdQr3t9chVLqA+/w035zSimTMvagc73lfUspFddkWW8rpYq8v5F1SqlhTcbNUUplebdvgVLq4abraTLdYaXUw0qpnd7lvKmUCmky/mdKqUJvvD9USmmlVP8zbQMArXWN1noNMA+4GLjau7yzlekOpVSed9z/qiZHV97fzTve300VcKdSKlop9ZI3xgKl1BNKKXOT5f1AKbXXu62XK6XSzxb7BdFaB8QLI8PXADHez0HAcWBcC9OmY1RUkd7PZqAQuAgIB6qAQU2WO6yNMcwEioAk4EXgnTNM+yrwRJPPvwW+8c6bCGwE/s87bhrgAv4IBAOhZ1juncD6Jp9vA+K92+Mhb3wh3nGfAQubTPs08LfmywHigArgdu9yvuv9HO8dvwY4Agzzjre0ENdhYIb3fS+MPdj/AwZ6v7crAQvwMyAHsLYw39nidQJ3e7/PhcAxQHnHrwOeA0KA0UAJcLl33GNAPXCVN/7XgUPAr7wx3Q0carLeNcAPve/7e2MP9n5v64C/tlLur4Hbve8jgIta+Q4vx2geHOtd7t+AdS0ts5X5m67zbNvlU+BNINZb1sta+80B/43xG03zDvsn8EaT9f4AiPSO+yuwo8m4QmCq930sMLbJevKbxb4JSMX43e0F7vWOm4Xx+x0GhAH/ATTQv5Xt0Pg9NRu+Dvij932rZQKGAnbgEsAKPOXdlie27WPez9di7GyHAu97lxGO8b+8CfiRd/r5GL/tIRi/s18DGzu0XuzIhXe1F8Ze5d3e93Mxmlham3Y9cIf3/ZVArvd9OFAJ3MAZKtozLPdvwC6gAG8F2cp0r3JqAsgF5jT5fBVw2Pt+GuDAW3GfZf130iQBtDC+Ahjlfb8A2OB9b/b+c01svhyMin9Ts+V8Ddzpfb8G+O1Z4jrs/WeqBPIwKuNQ4H+Bt5pMZ/Juu2lN5pvRxnhzmiwnDKNySMFIOG68Cd87/g/Aq973jwErm4y7xhur2fs50rusmCblPa1i8Y67FtjerNwn4l8HPA4knGVbvQT8qcnnCIyKJqP5Ms+wrZsmgNa2Sw/AA8S2sIzTfnMYlfEVTT738MYV1ML8Md71RHs/HwF+hHHOq/l6mieA25p8/hPwgvf9y8Afmozrz/klgKXAi2crE/AbTk1wYd5t0jQBNE3MyUADTeoNjJ2l1d73nwN3Nfut1wLpZ/u/Pt9XwDQBeb2GsceL9++/zzDtEowvB+AW72e01jUYFc29QKFS6lOl1OBziGERRtPGq1rrsnOYLxWjYjwhzzvshBKtdf05LA8A7+H0Xu/hdCUQDSR4R38IDFVK9cFIgjat9aY2xHYivp5NPh9tQzjXaq1jtNbpWuv7tNZ1zZettfZ4l9WzhfnPFm9Rk+XUet9GeNdRrrWuPkP8xU3e1wGlWmt3k88nlnUKZVzNtdR7uF+FsVea0Hw6r7swjnj2KaMZbW4r0zXfJnagjJa3SVu0tl16YWyXilbma/6bSwfeV0ZTXSVG5ekGkpVSZqXUk96mlCqMihxObosbgDlAnlJqrVLq4rbEi1FBntjuqZz6O2vLb64lPYHys5Wp+fq82675/3TTGNIxjqIKmyzvnxhHAifGP9NkXDmgOP/v9awCLQF8AIxUSg3HOAJYfIZp3wamKaXSgOvwJgAArfVyrfWVGHsD+zCac87K29a3CKMJ4b62tE02cQzjB3JCb++wxrDOYVkn4pmK0aTyHYy9vBjAhvGjw/vP/RZGsryd1hNm89hOxFdwIfG1tGyllMKomAqaT3gO8ba0jjilVGSTYc3jP1+/xyj7CK11lDe2Fk8waq2ztdbfxagQ/gi8o5QKbyXeptskHKMZrz3ibeooxnaJaWV88+/0KDDbm8RPvEK0cU7nFowmjhkYOxkZJ8IH0Fpv1lrPxyj7Bxjf47kqxGiqOaHXuS5AKdULGAd85R10pjKdsj6lVCjG99BU0210FOMIIKHJsqL0yQsxjmI0BzVdV6jWeuO5lqOtAioBeCuIdzAq801a6yNnmLYE4xDxFYz23b3QuEc33/tP14DRFOBpYwiPYPwgfgD8GXi96Qmgs3gD+LVSKlEplYBx+PmfNs7bmkiMdtwSIEgp9Rsgqtk0r2M0E8yj9Qr1M2CgUuoWpVSQUmoBRvvoJxcYHxgVwdVKqSuUcfL5IYzt3to/RVviPYXW+qh3eX9QSoUopUZi7I1f6PYFYxvbAZtSqifw09YmVErdppRK9B7lVHoHt/TbegP4vlJqtDJO8v4e+FZrfbgd4m2ktS7EaJZ4TikVq4yT6ZeeYZYXgN+dOHHp/a3O946LxPjeyjCaSn5/YiallFUpdatSKlpr7cQ4x9bW/6mm3sLYLkOUUmEYzYdtopQKU0pdhnEUuQnjN322Mr0DXKOMCyCsGE0+rV495N2eK4D/p5SK8p5g7udd74l1/VJ5T457Txjf1NYynI+ASgBerwEjaFvlsARjj2VJk2Em4EGMvbBy4DKME2cnblqxt7QgpdQ473x3eJsO/oiRDH7RxrifALYAOzHOIWzzDrsQy4FlwAGMJoV6mh02a603YPwzbtNaN2/mOTFNGcYR1UMY/+A/A+ZqrUsvMD601vsx9pr/hnHi8xqMy0UdrUx/1nhb8V2MvdJjGCfqHtVaf3EBoZ/wOMbJWhvGCdX3zjDtLGCP9zf0DHCztxnsFN64/hd4F2MvtB/G5Ysd4XaMNu99GBdN/OQM0z4DfASsUEpVY5w8neQd9zrGb6wAyPKOa76ew97moXuBs15Z15zW+nPgWWA1xsnUE+toOMNsf/fGWoxxYvpdYJY3CZ+xTFrrPcCPMc4ZFGIk+uNnWd8dGCeMszDOt72D0ZKA1vp9jHphqXc77AZmt7KcdnHiTH/AUEr1xvgxp2itq3wdT3eglFoFLNFa/8vXsbRFd4tXdAxlXC68GwjWWrs6YX0RGEduA7TWhzp6fe0hoI4AlHED0oPAUqn820YpNQFjD/ZNX8fSFt0tXtG+lFLXKeOemViMvemPO7LyV0pd420+Cse4DHQXJ09wd3kBkwC8X1AVxtUhj/o4nG5BKfUa8AXwk2ZXyHRJ3S1e0SF+hNEMk4txtc7CDl7ffIxmw2PAAIxmu27TrBJwTUBCCCEMAXMEIIQQ4lRd7uFEZ5KQkKAzMjJOGVZTU0N4eEuXSndP/lYe8L8y+Vt5wP/K5G/lgQsr09atW0u11onNh3erBJCRkcGWLVtOGbZmzRqmTZvmm4A6gL+VB/yvTP5WHvC/MvlbeeDCyqSUavGSaGkCEkKIACUJQAghApQkACGECFDd6hxAS5RSHDp0iPr6c34QZpcUHR3N3r17WxwXEhJCWloaFstp/bEIIcQ56/YJIDw8nMjISDIyMlDt34tbp6uuriYyMvK04VprysrKyM/Pp0+fPj6ITAjhb3zeBOR9Tvh2pdR5PTnSbDYTHx/vF5X/mSiliI+P95sjHSGE7/k8AWB0udZym0cb+Xvlf0KglFMI0Tl8mgC8na1cDchTG4UQogUuh5uv3jyAs7b9H9vj02cBKaXeweh7NRJ4WGt9Whd4Sql7gHsAkpOTxy1duvSU8VFRUQwYMKATou0cbrcbs7n1PmJycnKw2WydGNGFs9vtRESc1ltit+Vv5QH/K5M/lef4Lg8leyDponoSM8LOaxnTp0/fqrUef9qIjups+GwvjA5EntMnO37+5GzzjBs3Tje3bdu204Z1tmeeeUYPHjxY33LLLfr999/Xjz/+uNZa6/fff1/v2bOncbqHHnpIf/nll2dcVlVV1RnHZ2VlXXjAnWz16tW+DqFd+Vt5tPa/MvlLeSqP1+jn71+tl7+464LKBGzRXaxT+CnAPKXUYYwedS5XSrVHF3yd7rnnnmPlypUsXryYP/3pT9x3330AfPDBB2RlZTVO9+Mf/5gnn3zSV2EKIboRrTVfvZmNyayYfEPHtHL47DJQrfUvgV8CKKWmYTQB3Xahy/2H2t3quGn/TGXYPXEA7FlUzpofHWt12vv18Dat79577+XgwYPMnj2b2267jeDgYBISEti4cSMfffQRa9eu5YknnuDdd9+lX79+lJWVUVRUREpKyrkVTAgRUA7vLCVvdxmTb+hPRGxwh6yjK1wF1K298MILpKamsnr1apKSkhg7diwAkydPZt68efz5z39mx44d9OvXD4CxY8eyYcMGX4YshOjinA43X72VTWyPcEZentZh6+kSN4JprdcAa9pjWW3dcx92T1zj0UB7KSwsJDHxtCeuniIpKYljx1o/8hBCiK2fHaa6rJ5r/2cMZnPH7afLEUA7Cg0NPeuNWvX19YSGhnZSREKI7qbsmJ3tK44w6KIUeg6K7dB1SQJoR0OGDCEnJ6fxc2RkJNXVp3ZNe+DAAYYPb9tRihAisGiPZu2S/VhCzUy5oX+Hr08SQDu69NJL2b59+4nLXLn55pv585//zJgxY8jNzcXpdJKTk8P48adfjiuEEHs3FlKYY2PKDf0JjbR2+Pq6xDmA7u7w4cON72fMmMGXX37JjBkzmDJlyimXgb7//vvceOONBAXJZhdCnKq2ysHG93JIHRDD4It7dMo65QignT3yyCPU1ta2OM7lcvHQQw91ckRCiO5g/VsHcDa4ueyWQZ323C/ZFW1nycnJzJs3r8VxN910UydHI4ToDg7tLCV7y3EmXtOHuB6d15m9HAEIIYQPNdS5WLtkP/E9wxl7VXqnrlsSgBBC+NDX7+VQa2tg+u1DMAd1bpUsCUAIIXykYH8Fe746xqgrepGcEdXp65cEIIQQPuCod7Hq33uJSgxl4ry+PolBEoAQQvjA1+/nUlVWzxV3DMFibb0PkI4kCcAH6urquOyyy3C73QDk5+fz5ptv4nA4mDVrFi6Xy8cRCiE6Uv6+cnavLWDU5b1IHRDjszgkAfjAyy+/zPXXX9/Y89eXX37Jtm3bsFqtXHbZZbz55ps+jlAI0VEc9S5Wvb6PmOQwJs33TdPPCZIA2sl3v/tdFixYwMSJE0lPT+fTTz9tddrFixczf/58ANavX8+DDz7IO++8w+jRo5k7dy6LFy/urLCFEJ1sw9vZ2CvqueJ7vmv6OcGvbgR7/OM9ZB2ratdlDk2N4tFrhp11uszMTObPn8+bb77ZWKlfffXVp03ncDg4ePAgGRkZAFxyySVMmDCBp556iuHDh1NZWcnmzZvbtQxCiK7hUGYJWRsKGXtVOil9o30djn8lAF+pr6+npKSERx99FIChQ4dSUVFBTU0N9913H1arlWnTpnHrrbdSWlpKTEzMKfPv37+fwYMHA2A2m7FarVRXVxMZGdnZRRFCdJDaKger/7OPhF4RTLymj6/DAfwsAbRlT70j7N69mwEDBhASEgLAtm3bGDVqFO+99x433ngj11xzDQsWLODWW289rc+A0tJSoqOjT3lAXENDQ+OyhBDdn9aa1f/ei6POzfz/GdrpN3y1pmtE0c1lZmZy5MgR6uvrqamp4dFHH+V//ud/yM/Pp1evXgCNJ3xjY2Nxu92NSeDw4cOkpqY2LqusrIyEhAQsFkvnF0QI0SH2fHWMw7vKuPi6fsSnRvg6nEaSANpBZmYm119/PZMmTWLChAksXLiQKVOmkJaWRn5+PgAej6dx+pkzZ7J+/XoABg8eTGlpKcOHD2fjxo189dVXLZ47EEJ0T+XHatjwdja9hsQycnrH9e97PvyqCchXMjMzWbRoEc8+++wpw6+//noeeOABPv30U6655prG4ffffz9PP/00M2bMICIigk2bNjWOe/LJJ3nqqac6LXYhRMdxOd2seGk3lhAzV9w5FGXqnMc8t5UkgHaQm5vLgAEDThseHh7OK6+8ctrwsWPHMn36dNxud2PTEBhXCF199dUMHDiwQ+MVQnSOje/lUlZQw9X3jyQ8OtjX4ZxGEkA7ONHMcy5+8IMfnDbMarVyyy23tEdIQggfO7SzlF2r8xl1eS8yRiT4OpwW+ewcgFIqRCm1SSmVqZTao5R63FexCCFEe6our+fL17JI6BXBxdf183U4rfLlEUADcLnW2q6UsgDrlVKfa62/8WFMQghxQdxuDyv+tRuPW3PV3cMxW7rutTY+i0wb7N6PFu9L+yoeIYRoD9+8n0vRwSouv30IMUlhvg7njHyampRSZqXUDuA4sFJr/a0v4xFCiAtxKLOEHV8cZcRlPek/LsnX4ZyV0tr3O91KqRjgfeDHWuvdzcbdA9wDkJycPG7p0qWnzBsVFdXiFTjdVfMrg5rLycnBZrN1YkQXzm63ExHRdW5+uVD+Vh7wvzL5ojwN1ZqDKzTWCOgzQ2Eyt+8lnxdSpunTp2/VWo9vPrxLXAWkta5USq0GZgG7m41bBCwCGD9+vJ42bdop827fvt2vnplztmcAhYSEMGbMmE6M6MKtWbOG5t9bd+Zv5QH/K1Nnl8fpcPPuH7disdZz00MTiEoIbfd1dESZfHkVUKJ3zx+lVChwJbDPV/EIIcT50FqzZvE+yo7ZmfmDYR1S+XcUXx4B9ABeU0qZMRLRW1rrT3wYjxBCnLNdawo48G0xk+b1ofeweF+Hc058eRXQTq31GK31SK31cK31b30VS2eTLiGF8A8F+yvY8HY2GSPiGTcrw9fhnLOue4GqH5MuIYXo/qpK61i2aDfRSaHM+MGwLvecn7aQBNBOpEtIIQKHo97FZ8/vQmvNnIUjCQ7tEtfTnLPuGXVrPv8FFO1q32WmjIDZT551MukSUojAoD2aL1/bS/kxO3N/PIqY5K59s9eZyBFAO2itS8iDBw9y1113ceONNzZOey5dQgohup5vPjrIwe0lTL6hP72Hdq+Tvs351xFAG/bUO0JrXUL27duXl1566ZQEIF1CCtF97d1YyLZleQydmsqoK3p1yjod+QUU/vKXmObPa/dl+1cC8JGmXUK63W4effRR/vSnP7U4bdMuIUNCQqRLSCG6iYIDFaxZvI+0wbFcevNAlOr4k751u/dw9N570Q4H5qqqdl++NAG1g9a6hGyNdAkpRPdSUVTD5//cRXRiKLPuGY7Z3PFVp33dOvLuuAOT1UrGG0tw9u3b7uuQBNAOMjMzuf/++8nMzCQrK4vbbrsNMPbm7733XrZv384f/vCHxunvv/9+XnvtNYDGLiF3797N5MmTefvtt/nRj37kk3IIIU5XY2vg479lYjIprr5/FMFhHX90XvHWWxxdeB/WjHTSl75BcL+O6VNAmoDaQWtdQsbHx/PCCy+cNly6hBSie3DUu/j0Hzupq3Zw3UNjiU7s2Mc8aI+Hkr/8hbJ/vUT41Kn0fPppzBHhHbY+SQDtQLqEFML/uN0elr+4m9J8O3MWjiApPapD1+epq+PYz39B9YoVxHz3ZlJ+9StUUMdW0ZIAhBCiGe3RrHptL0f2lDP9tsEd3qevs7iY/PsfoH7PHpJ+8XPivve9TjnJLAlACCGa0Fqz/u1sDmwqZtL8vgy9JPXsM12Aul27yb//fjx2O2n/+DuRl1/eoetrSk4CCyFEE1s/P8zO1fmMuqIX42ald+i6bJ9+St5tt6GCgkh/441OrfxBjgCEEKJR5qqjfPvRIQZNSmHKDf07rBlGu92UPP00Zf96idBx40h79hmC4jv/rmJJAEIIAWStP8b6t7LpOzqR6XcM7rCne7ptNgoe/ik1X31FzM0LSHnkEZTV2iHrOhtJAEKIgLf/2yJWL95H72FxzLxrWIfd6FW/fz/5D/wYZ1ERKY8/TuyC73TIetpKzgEIIQJa9pZivnxtLz0HxjD7RyMwWzqmWrR9/DGHF9yMbmgg/bXXfF75gxwBCCECWPaWYla+nEVK3yjmLBxJkNV89pnOkcfh4PiTf6RiyRLCxo+n59N/ISgxsd3Xcz7kCMAHpEtIIXyvaeU/94FRWEPaf3/YkZ9P3i23UrFkCXHf/z69X3m5y1T+IAnAJ6RLSCF868Dmog6v/KtXreLQ9TfgyMsj7e9/I/nnP0N1saf8SgJoJ9IlpBDdQ9aGY6x8OYse/aI7pPL3OBwU/e735N93P9a0NPq89y6RM2a06zrai1+dA/jjpj+yr3xfuy5zcNxgfj7x52edTrqEFKLr27k6n6/ePECvoXHMvncElnZu83ccPkz+gw/SkLWX2NtvJ+mnD2Py0SWebeFXCcBXWusS8oMPPuDTTz+lqqqKu+66i5kzZ55Tl5CRkZGdXRQh/JLWmq3L8vj2w4P0GZXAVT8c3q5X+2itsb33HkW/+z0mi4W05/7R6Xf1ng+fJQClVC/gdSAZ0MAirfUzF7LMtuypd4TWuoS89tprufbaa6moqODhhx9m5syZ0iWkEJ1MezQb3s0h88ujDJiQzBV3DmnX6/zdNhuFjz5G9bJlhE2cSOofn8TSo0e7Lb8j+fIIwAU8pLXeppSKBLYqpVZqrbN8GNN5OVuXkE888QT3338/IF1CCtGZ3G4Pq1/fx/5vixg5PY1LbhrQrnf41mzcyLFHfoWrtJTEBx8k/q4foMztfylpR/FZAtBaFwKF3vfVSqm9QE+gWyaAE11COp1OHnnkEaZMmYLWml/84hfMnj2bsWPHNk5/okvIGTNmnNIl5KJFi8jNzZUuIYVoB456F8tf3M2RPeVMmteHcbMz2u3ZPp66Oo7/v79Q8Z//YO3bl4w33iB0xPB2WfbJlXigcDscWA77Pyc87Yftu3xAaa3bfaHnHIRSGcA6YLjWuqrZuHuAewCSk5PHLV269JR5o6KiWuyNqzPNnj2bZ5999rQ4nn/+ed544w3Gjh3LiBEjuOuuuwDYsWMH//jHP3jxxRdPW9Ytt9zC448/3mqZcnJysNls7V+IDmS324mIiPB1GO3G38oD/lemylI7ZVvDqK+EHuMVcf3ab6/fcvAgUa+9TlBxMbXTp1N93bXQTid6Te564sp3EF+2mbjyrQQ7KtCYsEUPZneP7+BKGXNey50+ffpWrfX400ZorX36AiKArcD1Z5t23Lhxurlt27adNqyz9ezZU7vd7nOa56WXXtIul+uUYQ0NDfqFF14443xZWVnnHJ+vrV692tchtCt/K4/W/lWmsmN2/c//+VK/8F9r9KGdJe22XHd9vS7+85911pCh+sC06dq+YUP7LLgiT+tvF2n97+u1/m2i1o9Gaf37NK3f+p7WO5ZqbS/VWl/YdwRs0S3UqT69CkgpZQHeBRZrrd/zZSwXQrqEFKJrOJJVxvIX96A1XPfgmHbrxrF2+3YKf/2/OHJzibnpRpJ+/nPM53vE5HFDwTY48LnRvFO82xge1xcm3AUDZ0H6ZDB3/HlAX14FpICXgL1a67/4Kg4hhH/YvTafdW9mE9cjjLix7nap/D01NRx/5hkq/v0fglJS6PXii0RMveTcF9RQDbmrjAr/wHKoLQVlht4Xw5X/B4NmQ0LnN2X78ghgCnA7sEsptcM77BGt9We+C0kI0d243R42vJXNrrUFpI+IZ+Zdw9j4zfoLXq79q68oeuxxnAUFxN5yC4kPPog5IrztC6jIgwPLjNfh9eB2QEg09L/SqPD7XwGhsRcc54Xw5VVA64GO7/VYCOG36qodLFu0m2PZlYye0YuLr++P6QIv83SVllL8+z9Q9dlnWPv2Jf0//yZs/OnnT0/jcUP+ZqPC378MSvYaw+MHwMR7jEq/10Vg7jr333adSIQQ4hyUHKnmsxd2UlftZMb3hzJoUsoFLU+73VS+9RbHn/4ruq6OhB8/QPzdd5/5UQ71Nsj50mjWyV4BdeVgCjLa8MfebrTnx/e7oLg6kiQAIUS3k7XhGOveOEBopIXrHx57we39dbt2UfT4b6nfvZuwiy4i5Te/Ibhvn5YnLsv1tuV/DnkbweMymnIGzDQq/H6XQ2jMBcXTWSQBCCG6DZfDzbqlB9i7sZBeQ2K58gfDCI08/2vwXeXllPz1GSrffhtzQjypTz1F1NVzTr1hzO2Co9+ebM8vPWAMTxwMFz9gVPppE7pU005bdb+IhRABqaKohuX/2kNZvp3xczKYMLfPebf3a5eLiqVvUvLss3hqaoi743YSHngA84kHMNZVeJt2lkH2SqivBJMFMqbA+Ltg4FUQ18oRQjciCUAI0eXt+7qQtUsPEGQxMfeBUaQPjz/vZdnXb6D4yT/gyMkl7OKLSPnVrwju3x9Ks2HXq8YJ3CNfg3ZDWDwMmmNU+P0uh5D2ua+gq5AE4AN1dXXMmjWLVatWYTabyc/PZ8OGDVx33XXMmjWLtWvXnvJ0UCEClaPOxbqlB9j/bRGpA2K48gfDiIgNPq9lNeTmcvxPf8a+di2W3r3p+czTRA4IRR14GT5bBuW5xoRJw+CSnxhNOz3Hgan7PNztXEkt4wMtdQmZlZXFggULGruEvPXWW30cpRC+VZhTycpXsrCX1zNhbh/Gz8k4ryYfV0kJJX/7O5XvvIMpLJSk22YS268S07Z74GsbmK2QMRUuWmjs6cf07oDSdE2SANrJd7/7XTweD4cOHaK4uJjnnnuu1ad6Ll68mCVLlgAnu4SMiYlh+fLlPPfcczzxxBOSAETAcrs8bP70ENuW5REZH8L1Px1HSt/oc1+O3U75yy9T9soraIeD2NHhJGTkEOTKhvwkGHoNDJwNfadBsP88CO9c+FUCKPr972nY275dQgYPGUzKI4+cdTrpElKIC1eaX80Xr+6lLN/O4ItTmLpg4Dn32euptRP/8evkPvTfuGscRPaqI2lkFdYBQ2Hgg0alnzoGTNIlul8lAF9prUvIvXv38swzz1BaWsoVV1zBwoULpUtIIVrgdnvYvjyPzZ8eJjjcwpyFI+gzKrHtC6gpRe/9nMq3llD6xUGCak0EpzhJunkwodOuM9rzo3t2XAG6Kb9KAG3ZU+8IrXUJOWTIEF544QU8Hg933HEHCxculC4hhWjmeF4Vq17fR1mBnf7jk7js5kGERJzlSZhaw/EsOLAMvfdzbBv2ULo7AmdNECG9Y2m4fBLpD/4BrGGdU4huyq8SgK+cqUvIjz76iOeff57bb78dkC4hhTjB2eBm0yeHyPziCGFRVmbfO4K+o8+w1+9qgMNfGZdpHliOrjiCLS+U0n0JOG0xBPdPJ+3BnxMxfRpr166Vyr8NJAG0g9a6hASYN28e8+bN4+qrr2581r90CSkC3aHMEta9eQB7eQNDL0ll8vX9CA5rYafHftx4xs7+zyF3NThr8KhQbNUjKNsUgfN4JcFDB5L2+/uJuPzyduvyMVBIAmgHmZmZLFq0iGefffaU4WvWrOG9996joaGBOXPmNA6///77efrpp5kxYwYRERFs2rSpcdyTTz7JU0891WmxC9GZqkrrWP92NocyS4lLDef6h4fRo3/MyQm0hqJdJ5+1U7DVGB7VE8/gG6g8GEHZxxtxHc8nZMQIkh/7HRHTp0vFf54kAbSD3NzcFvvwnTZtGtOmTTtt+NixY5k+fTput7vxXgAwrhC6+uqrGThwYEeGK0SnczncbFuex7YVR1AKLr6uH6Nm9MJsNoGzDg59dbKHrKoCY6ae42D6r3ElXUzFym1U/GkJbpuNsPHj6fGH3xM+ebJU/BdIEkA7kC4hhWiZ1pqcrcf5+r1cqsvrGTA+ick39CciqBJ2vG5U+LmrwVUHlnDoNx2m/RIGzKShtJbyV1/D9sED6IYGIq64gvgf3kXYmPPrGF2cThKAEKJDFB20seGdbIoOVhGfFsG1cyPp6fgE3l4Gx7YbE0X3hjG3waBZkDEVbbZS++0myn/+W+yrV6MsFqLnzyPuzjsJ7td1n6vfXUkCEEK0q8riWr75MJfcbSWEhWmmj9zB4Np/YVpxDFDGo5Ov+I1xQ1bSEFAKT309VR9+Qvnr/6Zh/37MsbEkLLyX2FtuISghwddF8luSAIQQ7aKmsoHNH+wm69tKzMrFhMgPGR36LtbKIONJmoNmG52mhJ+s0B1Hj1LxxlJs776L22YjeOBAevzuCaKuvhqT3AvT4SQBCCHOn8dDbfZWtn26n93ZiWitGB62nPGpXxM2bCoMXArpUyDoZKct2uXCvnYtFUvfpGb9ejCZiJwxg9hbbyFswgQ5sduJJAEIIc6NowZyV1O7axU7tlnYZbsUt05mUEIW46cEET3ubkh8CppV5M6CAirffZfKd9/DVVxMUFISCffdR8xNN2JJubD+fMX5kQQghDi7yqONXSLW5Oxme9Vs9tRdhVtb6D/AwYQbhhGbceVps3kcDuxffknlu+9Rs2EDAOGXXELK//6aiGnTUNLvBWBcLeWs9lBzzEXNMecpfwd8N5qUSR1zV7NPt75S6mVgLnBcaz3cl7EIIZrwuKFg28lr84t3U+lKZbvrdvZVLURrE4MmJjFuTl9ikk+tnLTW1O/eg+2DD7B98gkem42glBQSFi4k5obrsfQMrIeyuRuMit1e4K3YC5yYrYoR9xu9mrmdmpfi9uK0e1qcP2aA1T8TAPAq8HfgdR/HIYSor4KDq41n7WSvgNpSUGYKY29gh/VhDhZHYw4yMXRqD0bP6E10YugpszsLC7F98gm2Dz7EkZuLslqJnDGD6OuvJ/zii1Bm/+pZS2uNo8qDPd9JTb4Te4GL3ldFENHTeKTF1j+UsOMvpdSXuk+bN6qvpTEBmC0Kc4hCexThqRbCewYRnmohrIfxt8clHfdMI58mAK31OqVUhi9j8AXpElJ0GRWHYf8yRma+Aev2gMcJITG4+83koPkaMvcnUZxVS3BYEOOu6snIy3sRFnXyhK67spKqFSuo+vgTar39WISOG0fKbx8natYszFHdsw9drTUNlR7sR51ojyZxtJHs6kpdrLj5KDUFLuz5ztP22me/37sxAQDUl7pRZgjvEUR4T4v3FURUH+sp891xeBBBYarTT4ArrXWnrvC0AIwE8ElrTUBKqXuAewCSk5PHLV269JTxUVFRLT6GoStbtGgRLpeL++67DzB6CNu/fz+//e1v+f3vf0+/fv1YsGBBi/Pm5ORgs9k6M9wLZrfbiYjwnx6XunV5tJuoqgPEl20moXQz4bVHALCH9KA8cRKFERdzpGQA5bkKVx1YIyBukCK2D5iCjMpJ1dcTnLmTkC1bsGZlodxuXMnJ1E+cSP3ECbgTz+E5/h3kbN+Ruw5MVlDegxLbKiv2zVacJSacx004S8zoeqO8YSOd9HnG+J/zOGDvVScvY1UhGkuiG0uCh6BED3Hz6gkb5gLAVaXQTgiK0Y3r6cgyncn06dO3aq3HNx/e5RNAU+PHj9dbtmw5Zdj27dsZ0wVuDT+XLiEnT57MkiVLyMjIYP369cyfP5+YmBgiIyMbu4T87LPPWpx37969DBkypCOL0u7WrFnT4jORuqtuV556G+R8abTlZ6+AunIwBUH6ZBg4Gz1gJp+tOoqlOoncbcfxuDW9hsYxcnoa6cPiUSaF216Dfe0aqpctw752HdrhICglhajZs4maM4eQ4cO61OWba9asYfLYqeR9bqf6iAP7ESfVR5yNfxvK3dy0uR9J4409+3UPHGPXP8pPWYYlwkRELwvJE0O54tW0xuFHv7ATlhxERC8L1mhTp5X7Qn53SqkWE4BftTN89dYBSo/a23WZCb0imPqdsz+cTbqEFF1KWe7JJ2rmbQSPC0JjjRuxBs6C/ldQ5wpl/zdFZP2jkIpCjTW0jOGX9mT4ZT2JTQnHbbNh++hDqlespGb9eqPST0wk5jvfIWrWVYSOHYvyQbeKbocHe76L6sMOqvOcVOed+OskYUwIl/y/HgDUlbpZcfPRFpdhsirqSl2Nn/sviCZ+ZAgRvSyNL2tUy5V7rxnd9OivBX6VAHyltS4hAWpqarjssst47LHHmDt3rnQJKTqG2wVHvz151U7pAWN44mC4+AGj0u81EQ8mju4tZ9+/j3IwswSPS5PcJ4rUCYq5t02BsmKqv/iAvC+/oHbTZnC7CerRg5ibFxA1c2anVPpup8Z+1EnVIYdRyR92MurBBEJijXaU5d85yqEPq1uc1+M42aIRkRZE3+ujiOxtITLdQkRvo2KP7G0lNNGMMp2s3FOnhpM6NbxDy9UVnTUBKKV+DPxHa13R3itXSr0BTAMSlFL5wKNa65fOd3lt2VPvCK11CQnwxz/+ke985zuN00qXkKLd1FV4m3aWQfZKqK8EkwUypsD4u2DgVRDXB4CyAjv7PzjEgW+LqLE5CA4PYvjUngyZnEK47Si7X1nK0Zv/RsPevQBY+/Uj/q67iJxxBSEjRrRrM4f2aGqLXHhcmsjexsnQin0NrLm3gKpDxhU1utkVkRnXRBIy0bgaJqqP1ajI072vDGvj++h+J0+umq0mZr/bu93i9kdtOQJIBjYrpbYBLwPLdTudONBaf7c9luNrrXUJuXLlSoYOHXpKhS9dQooLUprtvSFrudG0o90QlgCD5hhP1Ox3OQQbR47V5fVkr8jjwKZiyvLtKJMifVgcU8bGkWDbR93616n4+zpKSkoIVwrT2LEk/fRhIqZfTnDfPu0SbsHaGkq312E76KDqoJOqgw6qDjlw12sGfDeamUt6AWAOURxbW2vMpIy998gMK1F9rET2sRCaeLKqmvKXFC55uke7xBfozpoAtNa/Vkr9LzAT+D7wd6XUW8BLWuvcjg6wO2itS8hf/epX1NTUkJWVRWhoKHPmzMFkMkmXkKLt3E448rW3H9xlUO79l0seDpf8xHiiZs9x4G2WqbE1cPDrfLK3FFOYY1y5kpQRycXTY0iu2Inrm9epfWErRU4npshIwi+ZQuS0aewwmbjsmmvOKbT6Cje2nAZsOQ6qch3Ych3YchzMfq93Y4Wd+ddSDn1wenNNSIIZc8jJo4qINAvzVmQQ2cdCZG8LZmvrzUxd6WRzd9emcwBaa62UKgKKABcQC7yjlFqptf5ZRwbYHbTWJeTvfvc7AF599VUSEhIwef9JpUtIcUa15UaTzoHPjSaehiowW6HPpXDRQqNpJ+Zk04a9op6DO0rI3VbCsZxK0BCbFMKoQU6SijajPluJ6/hxqoDgAQOIv/N7hE+dStiYMSjvkaZes6bFUOrLXVRmOzAHq8Zr4Y9vreOjmYdpKD/9BicAW46jMQFkXB1JeKqFqL5WovtaieprIaqPFWvUqddFmoIUva70n5Or3UVbzgH8N3AHUAr8C/ip1tqplDIB2UDAJ4DWuoQ84c477zzls3QJKU6hNZTsP3kC9+i3oD0QngRD5xsncPtOg+CTFWRFUQ2HMks5uKOE4kNVAMREK4bGHiM+ey2WNeuNRUdFETZ5MhFTLyF8ypQzPnTtyIpqir+to/JAA7ZsB5XZjsZKvt9NUcx6y0g6YclBNJS7CQo3Ed3fSnQ/a+PfqH5W4oYFNy5z6A/j2ntriXbUliOAOOB6rXVe04Faa49Sam7HhNW9SJeQ4py5HJC33qjw938Old5/r5SRMPVho9JPHdPYtONxeyjKruTwzlIO7SylsthoL48NrWOwcy8xu5cTZssHi4WwUaMI/+//InzyZEKGD0eZzXjcmurDDio+q6Zyf4PxOuBg5hu9CEs2qoGsRRXkvlt1SphB4SZiBhh77yeEpwZx57FBhKUESXNMN9eWcwCPnmHc3vYNRwg/VlNq3Ih1YBnkrAJHNQSFGHv3l/wEBlwF0ScflFZX7eBIVjl5u8s4sqeMhloXJqWJcxcx8Mg3JBRuIaShkuAhQwi//krCJ1+MZcgYXM5gwnsYTTvlWfUsX3AQW7YDd8Pp125U7m9oTAB9rosiqp+VmAFWYgYGEz3QSljy6ZW8MqnG5YvuTe4DEKKjaA3Hs4wKf/8yyN8MaIhIgeHXGz1k9bkMrMbljR63h+JcG0eyyjiyu5TjecZNjVZdT3zpLuKLdxBXsY+wjJ4Ej52AO/kX2BnBsbxQKlY3UP5cA/Yjh+l3QxSz3jGaa0Ligyjf3QAYe+4xg4KJGRRM7CArMYOCiR958nLjQbfGwK2duoWEj/lFAtBaB8ShqK8f2yHawFlvNO3s916qaTOetUOP0TDtF8YJ3B6jQSm01thK6sjfV8CRXccp2F+BwwFoD1H2I/Qp3UV8eRbRkUG4E0aT/N+3ETNtIkGJiSy76Qi5f6gCar0vg8mi8LhO/k5Ck8zctKUfMQOtWCP962mc4sJ1+wTgdrspKysjPj7er5OA1pqysjK5Qawrqi4+2bSTuxqcNRAUCv2mw6UPG49fiDKuW6+xNVCwuZgjOwrJ31dOTa3xmw2uLye+Yi9x5fuIcjmod/XjePEoMo/dhNNtPFFz/r0ZBCUaJ4KTxodSddBB3LAQYocGEzskmLghwUT1tTY+tA2MSyaTxoUiREu6fQKoqamhurqakpISX4fSLk7cINaSkJAQ0tLSWhwnOpHWULizsYcsCrYaw6PSYNQC49r8PlPBEkpdtYOCAxUc3byJ/AOVVNUa/3JBzlpiKw+QasuhZ6qJhDH9CRowlbfnX4Xbc/L570FhisQhwcQPDyE45uQe/NifJzL2575/6qbo3rp9AtBa06dP+9y12BWsWbOmSzzdVDTjrINDX8GBz7lo14ewtgxQxk1Y039t3IWbPJzaaicFu4s48tlmCg/XYKs3Lok0uxuIrswho/wI5kIXDcVJVNUOojT9KmZ+MaxxNYPuKiAyw0LcsBDiR4QQlWFBmYzmonp3Pcdrj2N32LE7jVeds45aVy21zlrq3fXUueqod9XT4G7A6XE2/nV5XDjdTtzajVu78WgPHu1BYzQXVVRUsGTFEpRSmJQJszIbL5MZi8mCxWTBarYSbA5u/BsSFEJoUCihQaGEBYURbglvfEVaI4mwRBBuCffrI/PurtsnACE6TFUhZC832vIPrgFnLVjCqY4eTshVj8PAq6hqiCR/0yHynj9CcdFhapTxQDGTu4Fo20EGmssIrTdTtDqKqtrBlDunEZ5qIXaSleQxDtQwOxuPbaSsrozy+nLKf1ROeX05lfWV2PbbqNxZia3BRpWjCpfHdeZ4vYJUEMFBwQSbgxsr7iBTkPFSQY0VvFIKhdEJiVM7qXXVorVuTBBu7W5MHi6Py0gmbiOpODyONsViVmYirZFEB0cTHRxNTHAMMcExxIXEERsSS1xIHPEh8SSEJpAQmkBsSCxBJqmWOotsaSFO8HigKPPktfmFO4zh0b1g9K3ogbMoU8PZ+PY6TO+EUFz8LQ0Wo8IPcrmIqswjqsQGRUEkz+lFxC9HU2iqZsf+oxyalU9N7MfYgkspcRynrK4Ml3aBE1h5MoQgUxBxwUblGBMcQ/+Y/kQHRxNljSLKGkWkNfKUvetwSzhhQWGEWow98WBz8HlVoOf6rHmP9lDvMo446lx11DhrqHPVNR6Z2B3Gq8pRRZWjCluDDVuDjZLaErIrsimvL6fB3XDack3KRHxIPIlhiSSFJZEclkxKeAop4Sn0CO9BangqiWGJkiTaiWxFEdgctcbe/YkHrNmLAAW9JuK67DfklY4hd3MdZTvsVFkacAVlAQlYHVVEVB4lrLqcck8l+xJrODKqDltSKeUhx3HigK9OriYsIowUawpJoUn0jZ9EUlhS415vQmgC8SHxxIXGEWmJ7BZNJiZlIswSRpjl/Pqr1VpT56qjrK6MsvoyyurKKK0rpaSuhJK6EopriymwF7C1eCvVjlOfJRSkgkgOTyYtIo2ekT3pFdmLtMg0ekf2Jj0qnXBL4D3W+XxJAhCBx1ZwssI/tBZc9WCNpL7XTHI8Ezi4N4rK7CBqvonBYwYIJVRXE1qZhSOxkKyYfWxKO4qtyaNroqxR9IzoycjIIfSMmEHPiJ70CO9Bj4ge9AjvQaRV+nZoSinVmEB6RfU647S1zlqKaosotBdyrOYYhfZC8u35FNgLWHN0DeX1p/bkFR8ST4yOYc3GNfSJ7kOf6D70je5LakQqJtX5Hdh0ZZIAhP/zeODYdm+l/zkU7UJrKPOMJKf0bko86ZR7YrDnxYIyoSxuQp1HUXW7OR6Vw97UwxxIrSbYHEGfuAxCG5K5dcA1pEemkx6VTlpkGtHB0b4upd8Ks4TRN7ovfaP7tji+xllDfnU+R6qPkFeVR15VHjuP7mTVkVVUNJzsxiTEHEKf6D70j+lP/9j+9I/pz8DYgSSHJXeLo66OIAlA+KcGu7dp53M4sAJPdQlFVX3JrZxKfsMCqkN64gyOAcDsrsPScAin5WuOxB4iu2cZ4QkZDEruz+De85gTY1Q+CaEJKKWM9vJR03xZOtFEuCWcQXGDGBQ3qHHYiXMalfWVHKo6xMHKg+TacsmtzOXbom/5+ODHjdNGB0czMHYgg2IHMThuMIPjBtM3pi8Wk/8/7kISgPAflUfgwHL0/mW4cjaQX9qX7LqxlHh+QnVQOu6gMAgBi6kSXLlUeQ5REm1DD0lhzKRRTIq9idtiB5AQmuDrkoh2EhMSw5iQMYxJOvXSaluDjZzKHLIrstlfsZ8D5Qd458A71LuNzpusJisDYwcyNH4owxKGMSx+GP1i+vndyWf/Ko0ILB43FGzFk/UZtVu+4GCBlYO1I6j0XEld8N14zFawgrW+EI9rB9WqhIbkUFJHj2Ti+DkM7TUQi9n/9/LE6aKDoxmXPI5xyeMah7k9bvKq89hbtpe9ZXvJKs/is0Of8daBtwAIDQplSNwQRiaOZETCCEYljiI5PNlXRWgXkgBE99JQDbmrqN22jB3fFpJf3R+7Hkx98KNokxmCPVjr89H1m6CHkx6z+zNs1EX0i1ngd3tvon2ZTebGcw1X9zV65fNoD0eqjrC7bDd7Svews3Qni/cuxulxApASnsKoxFGMSRrD6KTRDIod1K1+Z90nUhGwnMcOUfLxe2RlZlJSm0SdHkJdyPWgzCirG0t9HtR+AxFmeo4fzZSrriexd5SvwxZ+wKRMZERnkBGdwdy+RvcnDreD/eX7ySzJJLMkkx0lO1h+eDlgHCWMShzVeHQxMnEkwebgM63CpyQBiC7H3eCi+JMv2LXmc0qrLTR4BlAfMgJtGoOyurE25BFm2UrMsBh69p3IiAnTCY2Tn7LoHFazlRGJIxiROILbuA2Aopoith/fzrbibWw7vo3ndjyHRmM1WRmZOJIJKROYkDKBUYmjsJqtZ1lD55H/GuFzWmsaysrY++nrZG8qoqY6jfrgAXjM8yHYQ3DdUSz2zQRFRTH88umMve4KzMFyPbfoOlLCU5jdZzaz+8wGjJPM249vZ0vRFjYXb+afO//J85nPE2IOYUzSGCb1mMRFqRcxJG6IT+9NkAQgfKK+ws2eJSs4sG019Y5YnJZBOK2jAbCairDU7iAoxMqQqVcwat40QmQPX3Qj0cHRTOs1jWm9pgFQ5ahia9FWNhVt4pvCb/jrtr/CNmO6SSmTmJw6mSk9p5AS3nqfzR3Bp/9VSqlZwDOAGfiX1vpJX8Yj2pfHrTnyeTXH3wvlsL2K4BQHm/7zChVlNlymvjSE9AbLLMzYsTiyiU9tYNj8Kxg8QfpFFv4lyhrF9N7Tmd57OgCldaV8U/gNXx/7mm+OfcOKvBUA9I3uyyU9L2Fq2lTGJo3t8OYinyUApZQZ+AdwJZAPbFZKfaS1zvJVTKL9eNyaD2ccouTwboJGb2TVm2txWgfhsoyDMA8hdYcJq11HTFIqkxbcTI/x16BMgXk3pgg8CaEJzO07l7l956K1Jqcyh43HNrK+YD1v7HuD17NeJzQolIt7XMylaZcyNW1qh8ShfNXNoFLqYuAxrfVV3s+/BNBa/6G1ecaPH6+3bNlyyrC2PMXw8Y/3kHWs6kJD7hSVlZXExMT4Oozz5mnwkH7kG/qVHcfq7kt9aF9QZoKc1Zid+6i2lJGZOobaHukoc/es8Lv7d9QSfytTdy6Ph3pqTPuxm3ZTbd6JSxmPsxjYcCfv3vPQeS1TKbVVaz2++XBfNgH1BI42+ZwPTGo+kVLqHuAegOTkZNasWXPKeLvdftqw5vLzG6is8lxYtJ3E7XZTWVnp6zDOiap1clHp16TaTZjUEBpCxuAJBl1fgKVuHUfMmvX1F+NW47AkuLFEeKDa5uuwz1t3/I7Oxt/K1P3Lk04o6YQwB1dQEfXB+wjRPc9a150rXx4B3AjM0lr/0Pv5dmCS1vqB1uY53yOA7qQ7lEdrjaOykI2vv8zRnW7qPcNwBseBdhNSl4uFo8RGDOfY6yNx1Zz8fVkiFDPf6EXG3O59jX53+I7Olb+Vyd/KAxdWpq54BFAANH0ObJp3mOiKtKb4qy188/ablNsicASNwmWdjDI7CXHtI7Q+k8FT5jL6lnuwhJvwuDUfHTjM8W9rcdZ4sISbSJoURu/Z8lhkIboKXyaAzcAApVQfjIr/ZkAu/+hKXA4qN3zJhvc/pKQiEYd5DE7rHEzBDoIb9tGjr5Opd91KbI+rTpvVZFbMW57Bkc+r2fz+XiZcN4TesyMxddN2fyH8kc8SgNbapZR6AFiOcRnoy1rrPb6KR3jZS9DZK/j2i2Xsz0rBwQQcwd9BBTsJrt1HiCeX8fNuYtA1V2G2nLkyN5kVGXOjOBxRR8a07t3sI4Q/8ul9AFrrz4DPfBlDwNMajmfhyfqMXWtXsiu/D/XOi2kI/T5YPYTUZBPuymfcnHkMveFKzFa5A1cIfyG3VwYiZz3krYf9yyjcuZJ1+f2w2y+jPuSnEGQmxHWUlIQsJsy7kh7DLscSLpW+EP5IEkCgqC6G7BVwYBkNuatZZUujoGQmbtP/4bJEYDFXElKxkwETxnDpT7/n62iFEJ1AEoC/0hqKdnn7wV2GLtjKLiLYUjILt/0JHKEZKIuT4Kps4kKiuXThfJImRgZs36hCBCJJAP7EWQeH1nkr/eVQVUClycwnDZOxF/0UrcbgCQolRBURWbKbkdNnMvyuKwkKlSYeIQKRJIDurqoQspcbFX7uanDVoS3hfJ0ygU2HZhBWNQJnaD9MJgcxpkOMmjWYQdcswGw2+zpyIYSPSQLobjweKMqE/UbTDoU7jOHRvSkfvYB3S0Kp3mgm/OAELJYoTJQQfXQfI2ddycgfz/Jp6EKIrkUSQHfgqIWDa0427diLAAVpE9CX/y/bEvuwcsNWIt+NRlnGEGwyEWw7SEq9gzH/dRW9roiWtn0hxGkkAXRVtnyjsj+wzGjXd9WDNRL6Xw4DZ2HPmMJHBRvIensDaYejiAyZgdlcT3hBNj2TRjH2d98jZkDX7YtUCOF7kgC6Co8Hjm0n49Bi2Pdr4woegNgMGPd9GHgVpE/hoD2fJZlLqHv8K3rbLqWH9QaCVDnDUgoZd89MgiKvwhop7ftCiLOTBOBLDXY4uNrbtLMCao6Tjgl6XwRX/hYGzoKEgbi1h7X5a3n7kx+T8kUoyXXTcFsmY63PJ2bPEWYtuYHofvKQNSHEuZEE0Nkqj5zatON2QHA09L8CBs1mQ3EIl1w5D4BqRzXvZ/2bj7a+z9ivBzPSNQ+POZRw20GCD9Yy+AezGPpyAkEhchmnEOLcSQLoaB43FGyF/Z8bFf9x7/Pu4vrChLth0CzofTGYLQC41qzhaNVRFu9bzKrtX3D59ou5wvMjPCqIqOIczDkJDL33Wob+J04qfiHEBZEE0BHqqyB3lVHhZ6+A2lJQZkifDFf+HwyaDQkDTplFa82Okh28ePxFjvwnj7m7ZnCt+2eAoqc+THJMXyz9bmDE0ngsYVLxCyEunCSA9lJx+OS1+YfXg8cJITEwYKZxArf/FRAae9psbo+bL498yWt7XuPg0UNcu3MWU1zfAxQp1TlM+dEEkqdd2dmlEUIEAEkA58vjhqOb4IC3aadknzE8YRBctNA4gdtrEphb3sT1rno+zPmQ17Jeo7S0nBuz5jKl9g7ARFz+PhybexB343ySp6V2XpmEEAFFEsC5qLdBzpfGXn72CqirAFMQpE+Bsd8z9vTj+51xEbYGG0v3LWXx3sXU2GtZkHsdkWVj8SgLcYVZuDYlYkqbxZUfpNJjSngnFUwIEYgkAZxNWa5R4e//HI58DR4XhMbBgKtONu2ERJ91Mcdrj/P6ntd5+8Db1DsauOXYAmLyRuIyhRJzPAu9OQJ73Qzi7qzi2j/3l64ThRAdThJAc24XHP3m5FU7ZdnG8MQhMPnHRtNO2gQwte1mq6PVR3ll9yt8kPMBbo+H79q/S+KuodQTTlRNDuOnxrDj5SuJmxnO5D+l8O2e9VL5CyE6hSQAgNryk007OSuNph6zFTIugYn3wMCZxh255+Bg5UH+tetffHboM0zKxALrbaRu6E+VI4Igez6TMwoZ+dStmCPC6fM9jzySWQjR6QIzAWgNpdmNnaVw5BvQbghLgMFzjb38ftMh+Nzvrt1fvp9FOxexMm8lIUEh3Nbj+/RZN4CikhAa6stJP/wNhRsu4tjsfoyJMNr4pfIXQvhC4CQAtxPyNnrvwv0cyg8aw5OHwyU/gYGzoec4MJ1fZby3bC//3PlPvjzyJeGWcO4aeDdDto8iZ52DEncDfcs2UL52KHkV84jqY2HMwwntVzYhhDgPgZEA1v4ZNj4LDVVgDoY+U+Gi+4w9/ZheF7TofeX7eH7H86w6uopISyT3jriXiUWXsePfBWS7HfQs20p0ZSJZX84HTAz5QQyX/LWHPLBNCOFzPkkASqmbgMeAIcBErfWWDl1hZAoMnW/cgdvnMgiOuOBF7i/fz/OZz/PlkS+JtERy36j7uDJkHptfzmJT1XGibXlMTq9mx4oZFBwOxhplYtqingxYcPYrhoQQojP46ghgN3A98M9OWdvY241XOzhYeZDnMp9j+eHlRFgiWDhqITf0WkDmG7ks370fa0M1o51bGPvrWwkdNgy7pZgjy+3MXNqL6L7WdolBCCHag08SgNZ6L9Cteqk6Wn2U53c8z6eHPiXEHMLdI+7m9iF3kP+1jQ9f3IbLqUkv/orx1w0leNKvCR0SBsCEx5IY/+tEzMFyolcI0bUorbXvVq7UGuDhMzUBKaXuAe4BSE5OHrd06dJTxtvtdiIiLrxJpzU2l41ltmVstG/ErMxcGnkpM6JmYLaFU7ihnrqaYGIr9pNu2YVjzizy/pJGw2Ezff9ZiSX+3LdtR5fHF/ytTP5WHvC/MvlbeeDCyjR9+vStWuvxzYd32BGAUuoLIKWFUb/SWn/Y1uVorRcBiwDGjx+vp02bdsr4NWvW0HxYe7A12Hhp10ss2bcEt3Zz06CbuGfkPcSY4tj0QTaZa45hdTQw4vgnjP7J9dD/Lj65Oo+arAZCk8yM7jeRxNGh57zejiqPL/lbmfytPOB/ZfK38kDHlKnDEoDWekZHLbsj1TprWbx3Ma/sfgW7087cvnO5b/R9pEWmcXhXKZ+9tp4au4fUY+sZO0KR9syTlOcE8fGkXOqOu4kdGszcT9OJypD2fiFE1xYYl4G2gdPj5P3s93kh8wVK6kqYljaN/xr7XwyIHUBdtYMVL+4ke2sp4TWFTCxbwdBf/4jwyZM59lUNn849hKPKQ9qMcGa905vgaLnEUwjR9fnqMtDrgL8BicCnSqkdWuurfBGL1ppVR1fx161/5XDVYcYkjeGpy55ibPJYtNYc2FzEusV7cdS56HN4GSMnRtLj+ecwR4RTc8zJx7MO46rV9Lspiiv/k4bZKid7hRDdg6+uAnofeN8X624qsyST/7fl/7H9+Hb6RPfhmenPML3XdJRS1NgaWLt4H4d2lhFVdZhxxz9lwG/+m4hLL22cPzzVwqQnkqnIauCyF1LlIW5CiG4lIJuAjlYf5Zltz7D88HISQhP4zcW/4br+1xFkCkJrzf5vClm3dD+uOgf9cz9i6HALqc+9hDkmBgC3U2O2GJX96P9JQGvdrS5pFUIICLAEYGuwsWjnIpbsW4LFZGHhqIXcOexOwizGNfu1VQ7WLN7HocxSoqsPM/Tw2/T7+UKirrmmsYI/+qWdNT86xtxP0okdHAx0r/sZhBDihIBIAE6Pk7f3v83zmc9ja7Bxbf9reWDMAySFJTVOk7v9OGv+vQ9HbQP9cz9kYEIFPd/8F9a0no3TFH1dy2fX5OGq0+x7rYKL/9DSVa5CCNE9BEQCeGzjY3yU+xGTUibx0wk/ZVDcoMZxDXUuvnrzAPu/KSKqoYhRu14i/Y5rSbhvISro5OapzG7gU2/lP/jOGC76XbIviiKEEO0mIBLAbUNu48r0K7ks7bJTmmuOZVew8pUsasrr6XNkOf1qtpD23JOET5x4yvx1JS4+np1HfZmb9DkRTH+xJ8okzT5CiO4tIBLAkPghDIkf0vjZ7faw+eNDbF2eRzh2xm57ntTRvUj947sExcaeMq+rzsOn8/KoynWQMCaEmW/2whQklb8QovsLiATQVGVxLStf3sPxvGrSqjPpu/Pf9PivhcTfdReqhc5gDn1YRfE3dUT0tjD3k3SsEXKTlxDCPwRUAtj3TSFr3ziAyeNixIHXSHEeIu3lfxI2YUKr8wy4OQYVpIgdFEx4qqUToxVCiI4VEAnAUe9i7Rv7OfBtMQlWGwPX/YnYEf3o+Zd3sSQltThP02v7+98onbgIIfxPQDy3YM3i/WRvKmZgw1ZGrPgVqTfPJ/2VV1qt/OtKXbx/6SGKN9V2cqRCCNF5AuIIYPQgBzFvv0JU0W56/PlPRF8zt9VptdasuecYhetr2fjTIq5d00du9BJC+KWASACeVZ8Q5ywk7Y0lhAwZcsZpDyyxcfD9KiyRJq54PU0qfyGE3wqIBJD8i5+T8OMHTrvEszl7gZOvHjgGwCV/7UFUujzTXwjhvwIiASirlSDrmStzrTWr7y6godJD+pwIhnw/pnOCE0IIHwmIk8Btse/VSo58bic41mzc6StNP0IIPycJwMscrAhJMDP12R5yvb8QIiAERBNQWwy8JYaMuZFYIiUnCiECQ8DXdlrrxvfWKLM0/QghAkbAJ4DlC46y+f+O46rz+DoUIYToVAGdAA5/Vk3u21XseKoUR7UkACFEYAnYBKC15ptHigEY/5skwpLkdIgQIrAEbALI+7Sassx6wnoEMeL+OF+HI4QQnc4nCUAp9Wel1D6l1E6l1PtKqZjOXL/Wmi2/KwFgzMMJBIUEbB4UQgQwX9V8K4HhWuuRwAHgl5258oLVNRR/U0dIvJlhP5K9fyFEYPJJAtBar9Bau7wfvwHSOnP9B9+vAmDUT+KxhMvevxAiMKmm18H7JAClPgbe1Fr/p5Xx9wD3ACQnJ49bunTpKePtdjsRERHntE6toWaLhdAhLswRvi1/c+dTnq7O38rkb+UB/yuTv5UHLqxM06dP36q1Hn/aCK11h7yAL4DdLbzmN5nmV8D7eBPR2V7jxo3Tza1evfq0Yd2Zv5VHa/8rk7+VR2v/K5O/lUfrCysTsEW3UKd22LWPWusZZxqvlLoTmAtc4Q2ww9WXu3BUeYjKkMc8CyGEr64CmgX8DJinte60fhezXqrg330O8O2jxZ21SiGE6LJ8dQb070AksFIptUMp9UJHr1Brzb6XKwFIGhfa0asTQoguzye3v2qt+3f2Oou/raNiXwOhyUH0nh3Z2asXQoguJ2Cugdz7cgUAg26PwWyRJ34KIURAJABnjYfspTYA6epRCCG8AiIBHHzPhrPaQ/JFocQNDfF1OEII0SUERAKoznNisiqG/CDW16EIIUSXERDPQB7/6ySGL4zDLA99E0KIRgGRAABC4gOmqEII0SaySyyEEAFKEoAQQgQoSQBCCBGgJAEIIUSAkgQghBABShKAEEIEKEkAQggRoCQBCCFEgPJ5n8DnQilVAuQ1G5wAlPognI7ib+UB/yuTv5UH/K9M/lYeuLAypWutE5sP7FYJoCVKqS26pc6Ouyl/Kw/4X5n8rTzgf2Xyt/JAx5RJmoCEECJASQIQQogA5Q8JYJGvA2hn/lYe8L8y+Vt5wP/K5G/lgQ4oU7c/ByCEEOL8+MMRgBBCiPMgCUAIIQJUt0gASqlZSqn9SqkcpdQvWhgfrJR60zv+W6VUhg/CPCdtKNOdSqkSpdQO7+uHvoizrZRSLyuljiuldrcyXimlnvWWd6dSamxnx3gu2lCeaUopW5Pv5zedHeO5UEr1UkqtVkplKaX2KKX+u4Vputt31JYydZvvSSkVopTapJTK9Jbn8Ramad+6TmvdpV+AGcgF+gJWIBMY2mya+4AXvO9vBt70ddztUKY7gb/7OtZzKNOlwFhgdyvj5wCfAwq4CPjW1zFfYHmmAZ/4Os5zKE8PYKz3fSRwoIXfXHf7jtpSpm7zPXm3e4T3vQX4Frio2TTtWtd1hyOAiUCO1vqg1toBLAXmN5tmPvCa9/07wBVKKdWJMZ6rtpSpW9FarwPKzzDJfOB1bfgGiFFK9eic6M5dG8rTrWitC7XW27zvq4G9QM9mk3W376gtZeo2vNvd7v1o8b6aX6XTrnVdd0gAPYGjTT7nc/qX3DiN1toF2ID4Tonu/LSlTAA3eA/F31FK9eqc0DpMW8vcnVzsPVz/XCk1zNfBtJW32WAMxh5mU932OzpDmaAbfU9KKbNSagdwHFiptW71O2qPuq47JIBA9TGQobUeCazkZNYXXcM2jOerjAL+Bnzg23DaRikVAbwL/ERrXeXreNrDWcrUrb4nrbVbaz0aSAMmKqWGd+T6ukMCKACa7v2meYe1OI1SKgiIBso6Jbrzc9Yyaa3LtNYN3o//AsZ1UmwdpS3fY7ehta46cbiutf4MsCilEnwc1hkppSwYFeVirfV7LUzS7b6js5WpO35PAFrrSmA1MKvZqHat67pDAtgMDFBK9VFKWTFOfHzUbJqPgO95398IrNLesyRd1FnL1KztdR5G+2Z39hFwh/dKk4sAm9a60NdBnS+lVMqJtlel1ESM/6Uuu9PhjfUlYK/W+i+tTNatvqO2lKk7fU9KqUSlVIz3fShwJbCv2WTtWtcFne+MnUVr7VJKPQAsx7h65mWt9R6l1G+BLVrrjzB+BP9WSuVgnLi72XcRn10by/RfSql5gAujTHf6LOA2UEq9gXHFRYJSKh94FOMkFlrrF4DPMK4yyQFqge/7JtK2aUN5bgQWKqVcQB1wcxff6ZgC3A7s8rYxAzwC9Ibu+R3RtjJ1p++pB/CaUsqMkaje0lp/0pF1nTwKQgghAlR3aAISQgjRASQBCCFEgJIEIIQQAUoSgBBCBChJAEIIEaAkAQghRICSBCCEEAFKEoAQF0ApNcH7wL4QpVS49znuHfr8FiHai9wIJsQFUko9AYQAoUC+1voPPg5JiDaRBCDEBfI+z2kzUA9M1lq7fRySEG0iTUBCXLh4IAKjV6oQH8ciRJvJEYAQF0gp9RFGr259gB5a6wd8HJIQbdLlnwYqRFemlLoDcGqtl3if4rhRKXW51nqVr2MT4mzkCEAIIQKUnAMQQogAJQlACCEClCQAIYQIUJIAhBAiQEkCEEKIACUJQAghApQkACGECFD/H9XRka2vrEupAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import math\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "# Define f(t), f'(t), f''(t), ..., f^{(4)}(t).\n",
    "ft = lambda t: math.log(t)\n",
    "d1ft = lambda t: 1./t\n",
    "d2ft = lambda t: -1./math.pow(t, 2)\n",
    "d3ft = lambda t: 2./math.pow(t, 3)\n",
    "d4ft = lambda t: -6./math.pow(t, 4)\n",
    "\n",
    "# Define Taylor series approximation of f(t) for t near a.\n",
    "a = 1.\n",
    "p0n = lambda t: ft(a)\n",
    "p1n = lambda t: p0n(t) + d1ft(a)*(t-a)\n",
    "p2n = lambda t: p1n(t) + (d2ft(a)/2.)*math.pow(t-a, 2)\n",
    "p3n = lambda t: p2n(t) + (d3ft(a)/6.)*math.pow(t-a, 3)\n",
    "p4n = lambda t: p3n(t) + (d3ft(a)/24.)*math.pow(t-a, 4)\n",
    "\n",
    "# Define the points to sample along the x axis.\n",
    "numt = 100\n",
    "t0, tN = 0.1, 3.\n",
    "t = np.linspace(t0, tN, num=numt)\n",
    "\n",
    "# Define a list of functions to plot and their labels.\n",
    "funcs = [p0n, p1n, p2n, p3n, p4n]\n",
    "labels = [\"$p_0(t)$\", \"$p_1(t)$\", \"$p_2(t)$\", \"$p_3(t)$\", \"$p_4(t)$\"]\n",
    "\n",
    "# Plot the true function and a marker at (a, y).\n",
    "vft = np.vectorize(ft)\n",
    "plt.plot(t, vft(t), \n",
    "         linewidth=2., linestyle='--', c='darkviolet', label='f(t)')\n",
    "plt.plot(a, vft(a), marker='o', markersize=5.0, c='darkviolet')\n",
    "\n",
    "# Plot the Taylor series approximations of increasing degree.\n",
    "for func, label in zip(funcs, labels):\n",
    "    vft = np.vectorize(func)\n",
    "    plt.plot(t, vft(t), label=label)\n",
    "\n",
    "plt.title('y vs. x for Taylor Polynomials of Increasing Degree')\n",
    "plt.ylabel('y')\n",
    "plt.xlabel('x')\n",
    "plt.grid()\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 07.06 Piecewise Polynomial Interpolation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Fit large number of samples with low-degree polynomials to avoid excessive oscillations in the interpolant.\n",
    "* Break the interval $[t_1, t_n]$ into $k$ subintervals.\n",
    "* Each point at which interpolant changes referred to as **knot** or **breakpoint**.\n",
    "* At degree of 1, sacrifices smoothness to eliminate oscillation.\n",
    "\n",
    "#### Hermite Cubic Interpolation\n",
    "Hermite cubic interpolant is piecewise cubic polynomial interpolant.\n",
    "* Assume $n$ knots with $n-1$ piecewise subintervals.\n",
    "* For each knot, there are 4 parameters, $x$, in the cubic eg $p_{3}(t) = x_0 + x_1 t + x_2 t^2 + x_3 t^3$.\n",
    "* As a result, there are $4(n-1)$ parameters to be determined.\n",
    "  * Input samples give $2(n-1)$ equations at each subinterval.\n",
    "  * First derivatives give $n-2$ equations at each subinterval.\n",
    "  * This still leaves $n$ free parameters to be chosen.\n",
    "* To preserve monotonic data, then Hermite cubic is appropriate.\n",
    "\n",
    "#### Cubic Spline Interpolation\n",
    "Cubic spline is piecewise cubic polynomial interpolant that is twice differentiable.\n",
    "* Assume $n$ knots with $n-1$ piecewise subintervals.\n",
    "* For each knot, there are 4 parameters, $x$, in the cubic eg $p_{3}(t) = x_0 + x_1 t + x_2 t^2 + x_3 t^3$.\n",
    "* As a result, there are $4(n-1)$ parameters to be determined.\n",
    "  * Input samples give $2(n-1)$ equations at each subinterval.\n",
    "  * First derivatives give $n-2$ equations at each subinterval.\n",
    "  * Second derivatives give $n-2$ equations at each subinterval.\n",
    "  * This still leaves 2 free parameters to be chosen.\n",
    "* For maximum smoothness, then a cubic spline is appropriate."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Example: Cubic Spline Interpolation\n",
    "Determine the system of equations required to solve a cubic spline interpolation for 3 data points $(t_1, y_1), (t_2, y_2), (t_3, y_3)$.\n",
    "\n",
    "For $3$ knots there are $3-1=2$ subintervals and thus $4*2=8$ parameters which we refer to as $\\alpha$ and $\\beta$.\n",
    "\n",
    "The polynomial for the first interval $[t_1, t_2]$ given by coefficients $\\alpha$ and for the second interval $[t_2, t_3]$ given by coefficients $\\beta$:\n",
    "$$\n",
    "p_{[t_1, t_2]}(t) = \\alpha_1 + \\alpha_2 t + \\alpha_3 t^2 + \\alpha_4 t^3 \\\\\n",
    "p_{[t_2, t_3]}(t) = \\beta_1 + \\beta_2 t + \\beta_3 t^2 + \\beta_4 t^3\n",
    "$$\n",
    "\n",
    "The system of equations is given by:\n",
    "$$\n",
    "\\begin{bmatrix}\n",
    "1 & t_1 & t_1^2 & t_1^3 & 0 & 0 & 0 & 0 \\\\\n",
    "1 & t_2 & t_2^2 & t_2^3 & 0 & 0 & 0 & 0 \\\\\n",
    "0 & 0 & 0 & 0 & 1 & t_2 & t_2^2 & t_2^3 \\\\\n",
    "0 & 0 & 0 & 0 & 1 & t_3 & t_3^2 & t_3^3 \\\\\n",
    "0 & 1 & 2 t_2 & 3 t_2^2 & 0 & -1 & -2 t_2 & -3 t_2^2 \\\\\n",
    "0 & 0 & 2 & 6 t_2 & 0 & 0 & -2 & -6 t_2 \\\\\n",
    "0 & 0 & 2 & 6 t_1 & 0 & 0 & 0 & 0 \\\\\n",
    "0 & 0 & 0 & 0 & 0 & 0 & 2 & 6 t_3 \\\\\n",
    "\\end{bmatrix}\n",
    "\\begin{bmatrix}\n",
    "\\alpha_1 \\\\\n",
    "\\alpha_2 \\\\ \n",
    "\\alpha_3 \\\\ \n",
    "\\alpha_4 \\\\\n",
    "\\beta_1 \\\\\n",
    "\\beta_2 \\\\ \n",
    "\\beta_3 \\\\ \n",
    "\\beta_4\n",
    "\\end{bmatrix} =\n",
    "\\begin{bmatrix}\n",
    "y_1 \\\\\n",
    "y_2 \\\\\n",
    "y_2 \\\\\n",
    "y_3 \\\\\n",
    "0 \\\\\n",
    "0 \\\\\n",
    "0 \\\\\n",
    "0\n",
    "\\end{bmatrix}\n",
    "$$\n",
    "\n",
    "In this system:\n",
    "* Rows 1-4 ensure the interpolant is an exact fit at each data point.\n",
    "* Row 5 gives the first derivative at $t_2$.\n",
    "* Row 6 gives the second derivative at $t_2$.\n",
    "* Row 7 ensures the second derivative at $t_1$ is 0.\n",
    "* Row 8 ensures the second derivative at $t_3$ is 0."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use cubic spline interpolation to interpolate 3 ($m=3$) data points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.interpolate as interp\n",
    "\n",
    "# Input data $(t_i, y_i)$.\n",
    "m = 3\n",
    "t = np.array([-2,0,1], dtype='d')\n",
    "y = np.array([-27,-1,0], dtype='d')\n",
    "\n",
    "# Construct a cubic spline interpolator.\n",
    "cs = interp.CubicSpline(t, y)\n",
    "\n",
    "# Evaluate the interpolator at the data points.\n",
    "np.testing.assert_almost_equal(y - cs(t), np.zeros(m))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the observations against the fit using cubic spline interpolation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYcAAAEWCAYAAACNJFuYAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAAxW0lEQVR4nO3de7xWY/7/8deng5JSCQ0VCaUknRwmZYqGmKYSQ44dkONgjLNxmMF3kB8mh8kgOReZksOgaDuHJIcKCaVSlFS74967z++Pa+3c9r73+bDu+97v5+OxHvdhrfu+P9da616fta5rrXWZuyMiIpKoVtwBiIhI6lFyEBGRQpQcRESkECUHEREpRMlBREQKUXIQEZFClByqiJldb2aPxR1HWZnZ/8xsaNxxVAUz62VmX8T4+8WuE2Y2x8x6V19EYGZuZntFz8eY2TXV+fupxsx6m9niCnw+Y+ahkkM5mdkwM/vUzNab2TIz+7eZNYk7rrJItrFy96Pc/eG4YqoMZtY62uhlJwwfu/ub7t4uYbpvzaxvGb+7rZk9bWYrzGy1mX1iZhebWe2Kxu3u+7p7Vlk/Z2anm9nnZrbWzJab2Ytm1qgcv3+2u99Q1s+VIr5hZvZWKaet0Ma5OiUrV1XNwzgoOZSDmf0VuAW4FGgMHAzsDkw1s22qMY461fVbaaqJuzeMhv0r+mVmtifwHvAdsJ+7Nwb+BHQHyrwxrgxm9jvg/4AT3b0R0B6YEEcsqUD/iUrk7hrKMADbA9nA8QXebwj8CIyIXl8PTCT8UdcCs4D9E6a/HFgSjfsCODx6vxZwBbAAWAk8BewQjWsNOHA6sAh4A/gfcH6BWD4GBkfP/0XYmK0BPgR6Re/3AzYDOVF5Po7ezwLOSIjlb8BC4AfgEaBxgViGRrGsAK5OiOFAYGb0u8uB24uYn/OA/gmv60TzsStQH3gsmg8/Ax8AzUuxjPJjq1Pg/d7A4uj5o8AWYENU/stK8b2PAS8UM37r9ye89y3Qt5TrROK0tYGrovVgbbTsWiX5zUuAycXENA4YA0yNvud1YPeE8Q7slTDtjYllAf4aLfvvgeEJn6sH3BYt++XRb2xbRAzDgLcKlPMS4BNgdTQ/6gPbRctjS7RMsoFdKft/Iv+9kcDSKPZLCsR+ZzRuafS8XrJlmPC7a4G5wDHR++2BjUBeFOfPBedh9PpM4CvgJ2AKsGuBeX82MJ+wft8DWJzbt18tt7gDSLeBsFHNpcCGJxr3MPBk9Px6wob3OKBu9Gf4JnrejrDB3jVhBd8zen4hMANoGa3E9yV8Z/5K/0j0R9oWOA14OyGGDtGKlr+ynwI0I2x0/wosA+onxPhYgTJk8UtyGBGt2G0Iye+/wKMFYrk/imN/YBPQPhr/LnBq9LwhcHAR8/Na4PGE138A5kXPzwKeAxoQNpbdgO1LsYzyYysyOUSvvyXaGJdy2S8jYQOZZPyvvr/gbxS3TiSZ9lLg02hdsWj+Nkvym70IG9S/A4fkL/eE8eMIG7ZDo/XpX/x6Q11ccsgF/hHFejSwHmgajb+DsLHbgXDU9BzwzyLmyzAKJ4f3CRv+HQg7CGcXMw/L+p/If+/J6L39CDsc+fP2H9H37QzsBLwD3FDEOvInfklQJwDrgF2SlSvJPDyMsNPUNYr7LuCNAvP+eaAJsFsUY784t2+/KkvcAaTbQNjYLiti3M3A1Oj59cCMhHG1CHswvYC9CHtjfYk2DAnTzSM6iohe70LYoNRJWOnbJIxvFK2wu0evbwLGFhP/KqK9VUpODq8C5yaMa5cklpYJ498HhkTP3yBssHYsYX7uRdh4NYhePw5cGz0fEf1xO5VxGeXH9nPCcEmSP/63lC055BT35y34/QV/o7h1Ism0XwADSxnXUYSN88+EvdjbgdrRuHHA+IRpGxL2dltFr4tLDhtISLCEdfZgQrJaR7RDE437LfBNEfENo3ByOCXh9a3AmGLmYVn/E/nv7VPgNx6Mni8Ajk4YdyTwbVG/XyCW2fnLpWC5kszDB4FbC8z7HKB1wrzvmTD+KeCKsqzrVTmozaHsVgA7FlG3uUs0Pt93+U/cfQvhMH1Xd/8KuIiwsfjBzMab2a7RpLsDk8zsZzP7mfDHyAOaF/G9a4EXgCHRWycSNrAAmNklZjYvajz9mdBGsmMpy7oroUop30LCHzIxlmUJz9cT/gAQDvPbAp+b2Qdm1j/ZD0TzYh7wRzNrAAwAnohGPwq8DIw3s6VmdquZ1S1l7BASU5NouK0MnyvKSsIyroik60SS6VoRNmIlcvf/ufsfCXvhAwkbrTOK+M1sQhVHst8saKW75ya8zl++OxGO5j5MWE9fit4vraLWm2TK9J8o4r2F/FLmZOt10vlhZqeZ2eyE3+5IOf8/0bxfCbRImKYs86FaKTmU3buE6pPBiW+aWUPCHtyrCW+3Shhfi3BYvBTA3Z9w956EFd8JDdwQVuijEjZqTdy9vrsvSfheLxDTk8CJZvZbQt3t9Og3ewGXAccTqgOaEOp4rYjvKWhpFF++3QhVDctL+BzuPt/dTyQcut8CTDSz7YqY/ElCUhsIzI0SBu6e4+5/d/cOQA+gP6EarbKUVP6CpgHHFjN+HWGjCUB0BlPBDWaR60QB3wF7liU4d9/i7q8CrxE2Ysl+syEhiST7zdJaQTiq2DdhHW3s7pWxYUu2TMrzn4CEchPW3fwyJ1uvC80PM9udUG16PqFKrwnwGeX8/0TrfzNCW2PKU3IoI3dfTaguucvM+plZXTNrTTgkXEzY283XzcwGR0cZFxGSygwza2dmh5lZPUKjVn4jHISGvZuiFRMz28nMBpYQ1ouElfAfwIRojxRClVMuoS6zjpldS2hQz7ccaB1tpJJ5EviLme0RbVT+L/r+3CKm38rMTjGznaJYfo7e3lLE5OOBI4Bz+OWoATPrY2b7RRvZNYRD8qK+ozyWE9pTEuP+1syGFTH9dUAPMxtlZr+Jpt/LzB6LTmP+EqhvZn+IjnD+RqhrTpR0nUjyWw8AN5jZ3hZ0MrNmBScys4FmNsTMmkbTHQj8rsB3Hm1mPaMz6W4gVG0l29MulWiZ3g/cYWY7R3G0MLMjy/udCZYDzcysccJ75flPAFxjZg3MbF9gOL+cxfUk8Lfoe3YktHslu/5kO0IC+DH63eH8OukuB1oWc4bik8BwM+sc/df/D3jP3b8tReyxU3IoB3e/lXAmyW2EjVb+6Y2Hu/umhEmfJTRirQJOJZxBlEPYYNxM2ANbRti7vjL6zL8IDX2vmNlawp/8oBLi2URoLO5LwsaVUCXzEmGjtZCQiBI3Ck9HjyvNbFaSrx5LSHZvEBpONwJ/Li6WBP2AOWaWHZVpiLtvKCL+7wlHZD349WmYvyGc3bOGUJXwehRP/sVGY0oZS1H+SdhI/BxVv21D2LNLtrHG3RcQ6tZbR2VbDTxDOCtrbbTjcC5hw76EcCRR8Jz9otaJgm4n7HC8Qij/g4TG1oJWEc6ImR9N9xgwyt0fT5jmCUJi+4nQqH9KEfOjLC4nnKwww8zWEI6q2hX/kZK5++eEjerX0XLZlXL8JyKvRzG+Ctzm7q9E799IWGafEBr9Z0XvFYxlLvD/COvmckLD9tsJk7wGzAGWmdmKJJ+fBlxDWEe+JxwJDik4XaqyqCFEpMYzs57AeVF1WEYws3GEBta/xR1LdYmO5PPPAivxKFeS0wUjIhF3fwso1ZW8IplO1UoiIlKIqpVERKQQHTmIiEghGdHmsOOOO3rr1q3L/fl169ax3XZFnYKfPjKlHKCypKJMKQeoLPk+/PDDFe6e9OLFjEgOrVu3ZubMmeX+fFZWFr179668gGKSKeUAlSUVZUo5QGXJZ2YLixqnaiURESlEyUFERApRchARkUIyos0hmZycHBYvXszGjRtLnLZx48bMmzevGqKqWqlajvr169OyZUvq1i3LDVVFJE4ZmxwWL15Mo0aNaN26NWZW7LRr166lUaNYenmsVKlYDndn5cqVLF68mD322CPucESklDK2Wmnjxo00a9asxMQgVcvMaNasWamO4EQkdWRscgCUGFKEloNI+snYaiURkaQ2b6buzz/D/PmwejX8/HN43LgRNmwIj/lDXh64/zIA1KkD22wD9er98tigAWy//S9Do0bQpAk0bQq10nMfXMmhio0ePZp///vfLFu2jMsvv5wrrriCyZMn07ZtWzp06BB3eCKZIzcXli6FRYtg4cLwuHgxLF/+62HNGg6prphq14Ydd4SddgpD8+bQsiW0avXL4267wc47Q4odYadscjCzfoROPmoDD7j7zTGHVC733nsv06ZNo2XLllvfmzx5Mv3791dyECkr97DR//xz+PLLsPf/5ZdhWLQo7Oknato0bJCbN4cuXcLjTjsxf8UK9u7ePezdN2kS9vYbNID69cOw7bbhiKB27bDRTtxw5+XB5s1h2LQpDOvXw9q1sGbNL8OqVfDjj/DDD+Hxxx/hgw9g0qTwmUTbbQd77QV77hke99oLOnSAffcN8cUgJZND1C3kPcDvCT1pfWBmU6KemdLG2Wefzddff81RRx3FiBEjWLBgASeddBJTpkzh9ddf58Ybb+SZZ55hzz3L1FWwSM2wdi189BF88gl8+mkYPvssvJ+vUSNo2xYOPhhOOgl23z3sie++e9grL+KeQ0uysti7vLfPqF07JI9tk3XMVwrusGJFOKr57ruQ1BYsgK++grlz4fnnQ+LJt+uuIUl07Aj77w/dusE++7BoXR1GTIFB9eDGR2HsANitcdE/W1YpmRyAA4Gv3P1rADMbT9T5fLm+7aKLYPbsIkdvm5cXFnhZdO4Md95Z7CRjxozhpZdeYvr06Tz//PMA9OjRgwEDBtC/f3+OO+64sv2mSKbavDkkgg8+CMPMmTBv3i/1/E2bwn77wWmnhccOHWDvvcORQIpVx5TI7Jdqpi5dCo/PywtJY+5cmDMnJMQ5c2DMmNAmArDttmTvuj9Dd+lK04N3YUGr3oyYAtNOrbwwUzU5tODXfR0vpkCfsWY2EhgJ0Lx5c7Kysn71BY0bN2ZttIdRb/NmahU83EzkTm5x45PYsnkzmxL3YIr8aic7O5uNGzeyefNm1q5dS05ODhs2bNgaX2XJy8ur9O+sLBs3biy0jIqTnZ1dpulTWaaUpTLLUSc7m+3nzKHxp5/S+NNPafT559SO9pY3N23Kmn32Ye3QoWS3bcvavfZi8447Fq7a+fzzMJRDWiyTBg3ggAPCAJCXR4PFi2n05Zc0/PJLGs/5khM+HMeSugdz8QE9MaAyi5SqyaFE7v4f4D8A3bt394J3JZw3b94vF4Tde2+x31Xei8e2KcU0ZkbDhg2pX78+22yzDY0aNaJu3bpsu+22lX7BWipeBJevfv36dEm2l1QE3TUz9VSoHOvXw1tvwbRp8Oqr4SjBPZz506ULnHceHHIIHHQQ27RowY5m7Fip0f9aJiyTvo/C1yu3cFnTl7h9WW/2bArn966870/V5LAEaJXwumX0XkZo1KhRyu7hi1QK99BG8MIL8PLL8O67oeqobl347W/huuvg0EPhwAOLbBeQ4o0dACOm1CK3XgP2rBVeV6ZUTQ4fAHub2R6EpDAEOCnekCrPkCFDOPPMMxk9ejQTJ05Ug7RkhnXrwpHBiy+GYfHi8H7nznDBBXD44dCrl5JBJdmtcWhjyMqq3COGfCmZHNw918zOB14mnMo61t3nxBxWuXz77bcADBs2jGHDhgFwyCGHMHduWp14JZLcqlXw3HPw3/+GI4SNG8MZREccAX//Oxx1FOyyS9xRSjmkZHIAcPcXgRfjjkOkplu0ml+dMjmu50pavjIRnnkGpk8PF5+1bAlnngmDBkHPnuHKYUlrKZscRCQ1jJgCy75fS9vvpnLVa6P4zRevQF5uOJX0kktg8GDo3j39TimVYik5iEhyeXkwdSoX3jGOvp9NYducDSxp0ooHe1/MWbedGC7IUkLIWEoOIvJrX3wB48bBI4/A0qUc2rAZEw8azra/34fLtj+PNs1qcVbnuIOUqqbkICLhytunnoL77gunndauHRqTR49mda/+PPxyPQbVy6LNplqVfsqkpCYlB5GabN68kBAefjjcurpdO7j1VjjllK1nGe1G1Z4yKakpPW80nsa+/fZbOnbsGHcYzJ49mxdf/OVksClTpnDzzWl541spq9zccKZRnz7hHkX33gv9+oWt/7x5cOmlOv1UdOSQCXJzc6lTp2yLcvbs2cycOZOjjz4agAEDBjBggOoLMtqqVfDAA3D33eFOoK1bw803w/DhoT8BkQQ6cqhit99+Ox07dqRjx47cGd3FNTc3l5NPPpn27dtz3HHHsX79egCuuOIKOnToQKdOnbjkkksA+PHHHzn22GM54IADOOCAA3j77bcBuP766zn11FM55JBDOPXUUzn44IOZN2/e1t/t3bs3M2fO5P333+e3v/0tXbp0oUePHnzxxRds3ryZa6+9lgkTJtC5c2cmTJjAuHHjOP/884FwdHPYYYfRqVMnDj/8cBYtWgSEC/kuuOACevToQZs2bZg4cSIA33//PYceeiidO3emY8eOvPnmm9Uyb6WUvvwSzj03XItw2WXQpk3oU+Crr+Dyy5UYJDl3T/uhW7duXtDcuXMLvVeUNWvW+MKf3Q9/xH2Pf4XHhT+X+uNFmjlzpnfs2NGzs7N97dq13qFDB581a5YD/tZbb7m7+/Dhw33UqFG+YsUKb9u2rW/ZssXd3VetWuXu7ieeeKK/+eab7u6+cOFC32effdzd/brrrvOuXbv6+vXr3d399ttv98svv9zd3ZcuXept27Z1d/fVq1d7Tk6Ou7tPnTrVBw8e7O7uDz30kJ933nlbY0183b9/fx83bpy7uz/44IM+cOBAd3cfOnSoH3fccZ6Xl+dz5szxPffc093db7vtNr/xxhvd3T03N9fXrFlTaF6UZXm4u0+fPr1M06ey2Mry3nvugwe7m7lvs4378OHuH31U7q/TMklNFSkLMNOL2K6qWikyYgosWAVbPDxWxr3R33rrLY455hi2i+4lM3jwYN58801atWrFIYeEjgpPOeUURo8ezUUXXUT9+vU5/fTT6d+/P/379wdg2rRpv7rVxpo1a8jOzgZCVdC2UYcjxx9/PH379uXmm2/mqaee2tpXxOrVqxk6dCjz58/HzMjJySkx7nfffZf//ve/AJx66qlcdtllW8cNGjSIWrVq0aFDB5YvXw7AAQccwIgRI8jJyWHQoEF07ty5IrNNKsIdXnopNCpnZYVexK68MtzbqHnzuKOTNKJqpcjXUWKA8Pj1qqr7LStw4ZCZUadOHd5//32OO+44nn/+efr16xdi2bKFGTNmMHv2bGbPns2SJUto2LAhwNakA9CiRQt22GEHPvnkEyZMmMAJJ5wAwDXXXEOfPn347LPPeO6559i4cWOFYq9Xr97W5x51xHLooYfyxhtv0KJFC4YNG8YjjzxSod+QctiyJVQVdesGRx8dus/8f/8vtC3cdJMSg5SZkkOkTVOoFW2za1l4XVG9evVi8uTJrF+/nnXr1jFp0iR69erFokWLePfddwF44okn6NmzJ9nZ2axevZqjjz6aO+64g48//hiAI444grvuumvrd84upke7wYMHc+utt7J69Wo6deoEhCOHFi1aADBu3Lit0xZ32/AePXowfvx4AB5//HF69epVbDkXLlxI8+bNOfPMMznjjDOYNWtW8TNGKs+WLfD00+HOp4MHQ3Y2jB0LX38NF18cboInUg5KDpGxA2DPplDbwmNlXOjTtWtXhg0bxoEHHshBBx3EGWecQdOmTWnXrh333HMP7du3Z9WqVZxzzjmsXbuW/v3706lTJ3r27Mntt98OwOjRo5k5cyadOnWiQ4cOjBkzpsjfGzRoEOPHj+f444/f+t5ll13GlVdeSZcuXcjNzd36fp8+fZg7d+7WBulEd911Fw899BCdOnXi0Ucf5V//+lex5czKymL//fenS5cuTJgwgQsvvLA8s0vKYssWGD8+dJl5/PGQkwOPPRa6lhw+XDe+k4orqjEinYbKaJDOBKlcDjVIV5ItW9yffda9Uyd3cN93X/fx491zcyvvN4qgZZKaqqpBWkcOIuni1VdDL2oDB4ZuN594Aj75BE44IdzuQqQSKTmIpLoPPwy9qPXtC0uXwv33h+qjE0+EWvoLS9XI6DXLo7NpJF5aDuX07bdw8smhr4RPPoE77wwXtJ1xRuiLWaQKZWxyqF+/PitXrtSGKWbuzsqVK6lfv37coaSPVavC/Y3atQvdb155Zbia+cILQfNRqknGXgTXsmVLFi9ezI8//ljitBs3bsyIjVeqlqN+/fq0bNky7jBSX25uuEPqtdeGBDF0KNxwQ7jthUg1y9jkULduXfbYY49STZuVlUWXLl2qOKKqlynlqJFefTUcGcyZA4cdBrffHnpaE4lJxlYriaSFBQvgmGNCY/P69eEq52nTlBgkdkoOInHYsAGuuSb0pzB1Kvzzn+EMpEGD1C+zpISMrVYSSVnPPx9uhPfNN6HHtVtugV13jTsqkV/RkYNIdVm4MBwZ/PGP4ayj6dPh0UeVGCQlKTmIVDHLywu30G7fPlQh3XILzJ4NvXvHHZpIkVStJFKVZs6k69lnh+sUBg6E0aNht93ijkqkRDpyEKkK69aFW2YfdBDbrFoFzzwDkycrMUjaSLnkYGbXm9kSM5sdDUfHHZNImbzyCnTsCHfcAWeeyQfjxoW+FkTSSMolh8gd7t45Gl6MOxiRUlm9Gs48E448EurVgzfegDFjyI167hNJJ6maHETSy0svhaOFsWPh8stDg3MJPeiJpDJLtRvTmdn1wDBgDTAT+Ku7F+rR2cxGAiMBmjdv3i2/W8vyyM7O3tovczrLlHJA+pSlTnY2e957L7v873+s2313Pr/8cta2b/+radKlLCXJlHKAypKvT58+H7p796Qji+oFqCoHYBrwWZJhINAcqE04qrkJGFvS9yXrCa4sMqVXqEwph3ualOXVV91btXKvVcv9iivcN2xIOllalKUUMqUc7ipLPorpCS6WU1ndvW9ppjOz+4HnqzgckbLZsAGuuir0r9C2Lbz7Lhx4YNxRiVSqlGtzMLNdEl4eQziiEEkNs2ZBt24hMZx/Pnz0kRKDZKRUvAjuVjPrDDjwLXBWrNGIAOTlhSubr7sOdtopNEAfeWTcUYlUmZRLDu5+atwxiPzK4sXhBnmvvw7HHw///jfssEPcUYlUqZRLDiIpZdIkOP102LwZxo2D007TLbWlRki5NgeRlLB+PZxzTriyuU2b0LYwdKgSg9QYSg4iBc2dGxqZx4yBSy+Fd96BvfeOOyqRaqVqJZFEDz8M554LDRvCyy/DEUfEHZFILHTkIALhLqrDh8OwYXDQQeH2F0oMUoMpOYjkVyM9/DBce23okGeXXUr+nEgGU7WS1GyPPw4jR4ZqpFdegb6lunhfJOPpyEFqps2b4c9/DtcvdOsWqpGUGES2UnKQmmfJktB/8913h97aXn1V1UgiBahaSWqWrCw44YTQAD1hQrjiWUQK0ZGD1AzuodvOvn2haVN4/30lBpFiKDlI5tuwIVzdfPHF8Mc/hsTQoUPcUYmkNCUHyWzffQeHHgqPPgp//zs88wxsv33cUYmkPLU5SOZ66y049thwn6TJk2HgwLgjEkkbOnKQzHT//XDYYeEo4b33lBhEykjJQTJLXh785S/hwrY+fdS+IFJOSg6SOdasgQEDQheeF1wAL7wQzkwSkTJTm4Nkhm++CWciff556Knt7LPjjkgkrSk5SPp7+20YNAhyc8Nttg8/PO6IRNKeqpUkvU2YEBqemzYNDc9KDCKVQslB0pM73HwzDBkS+l94911o2zbuqEQyhqqVJP3k5ITe2h54AE48ER56COrVizsqkYyiIwdJL2vWQP/+ITFcfTU89pgSg0gV0JGDpI+lS+Goo2DOnJAcTj897ohEMpaSg6SHefOgXz/46Sd48UX17yxSxZQcJPW98064hqFuXXj9dejaNe6IRDKe2hwktT37bDg9tVmzkCSUGESqRSzJwcz+ZGZzzGyLmXUvMO5KM/vKzL4wsyPjiE9SxH33weDBsN9+4UK3Nm3ijkikxojryOEzYDDwRuKbZtYBGALsC/QD7jWz2tUfnsTKnd0ffTTcAqNfP5g+HXbaKe6oRGqUWJKDu89z9y+SjBoIjHf3Te7+DfAVcGD1Riex2rIFLr6YPcaOhVNOCf0wbLdd3FGJ1Dip1iDdApiR8Hpx9F4hZjYSGAnQvHlzsrKyyv2j2dnZFfp8qkj3clhuLu1GjeI3r7zC1wMGsGj48FCdlObSfbnky5RygMpSGlWWHMxsGvCbJKOudvdnK/r97v4f4D8A3bt39969e5f7u7KysqjI51NFWpdjwwY44QR45RW44QYWHXIIvfv0iTuqSpHWyyVBppQDVJbSqLLk4O59y/GxJUCrhNcto/ckk61ZE05VffNNuPdeOOccyJC9OpF0lWqnsk4BhphZPTPbA9gbeD/mmKQqrVwZTlV95x144omQGEQkdrG0OZjZMcBdwE7AC2Y2292PdPc5ZvYUMBfIBc5z97w4YpRqsGwZ/P73MH9+aHj+wx/ijkhEIrEkB3efBEwqYtxNwE3VG5FUu0WLoG/fcL+kF18MfTKISMpItbOVpCb46qtQlbR6dWiA7tEj7ohEpAAlB6lec+eGI4bNm+G113Q7DJEUpeQg1efTT8MRQ+3a4QZ6++4bd0QiUoRUO1tJMtXs2dCnT7izalaWEoNIilNykKo3a1ZocG7QIBwxtGsXd0QiUgIlB6laH3wQqpK23z4khr32ijsiESkFJQepOjNmhMbnHXYIiWGPPeKOSERKSclBqsZ774WuPHfeOSSG3XePOyIRKQMlB6l8H3zwS2KYPh1atow7IhEpIyUHqVwffhgSQ7NmSgwiaUzJQSrPrFnhXklNmoTE0KpViR8RkdSk5CCVY/bs0PjcqFFIDGpjEElrSg5ScXPmhMTQsGG4wK1167gjEpEKUnKQivnyy3AdwzbbhHsl6XRVkYygeytJ+X3zTUgMW7aEqiRd4CaSMZQcpHwWLw6JYd26kBjat487IhGpRCVWK5nZn82saXUEI2li2bKQGFasgJdfhv33jzsiEalkpWlzaA58YGZPmVk/M7OqDkpS2E8/hdNVFy8OPbgdcEDcEYlIFSgxObj734C9gQeBYcB8M/s/M9uzimOTVLN2LRx1VGiEnjIFevaMOyIRqSKlOlvJ3R1YFg25QFNgopndWoWxSSrZuBEGDQpXQE+YEKqVRCRjldggbWYXAqcBK4AHgEvdPcfMagHzgcuqNkSJXU4ODBkSTlV95JGQJEQko5XmbKUdgMHuvjDxTXffYmb9qyYsSRlbtsCIEfDss3DXXXDqqXFHJCLVoMTk4O7XFTNuXuWGIynFHS68EB57DG68Ec4/P+6IRKSa6AppKdoNN8Ddd8Nf/wpXXRV3NCJSjZQcJLkxY+C662DoUBg1CnQGs0iNouQghU2cCOeeC/37w/33KzGI1EBKDvJrr70GJ58MPXqEU1br1o07IhGJQSzJwcz+ZGZzzGyLmXVPeL+1mW0ws9nRMCaO+GqsWbNg4EBo2xaeew4aNIg7IhGJSVw33vsMGAzcl2TcAnfvXL3hCAsWhKufmzWDl16CprqdlkhNFktyyD8FVrdpShE//ABHHgl5eeFGei1axB2RiMTMwp0xYvpxsyzgEnefGb1uDcwBvgTWAH9z9zeL+OxIYCRA8+bNu40fP77ccWRnZ9OwYcNyfz5VlKcctTZsoPPFF7PdN9/w8e23s6ZDhyqKrmwyZZlA5pQlU8oBKku+Pn36fOju3ZOOdPcqGYBphOqjgsPAhGmygO4Jr+sBzaLn3YDvgO1L+q1u3bp5RUyfPr1Cn08VZS5HTo77H/7gXquW+7PPVklM5ZUpy8Q9c8qSKeVwV1nyATO9iO1qlVUruXvfcnxmE7Apev6hmS0A2gIzKzk8cYdzzoEXXgjXNAwYEHdEIpJCUupUVjPbycxqR8/bEG4V/nW8UWWof/wDHngArr4azjor7mhEJMXEdSrrMWa2GPgt8IKZvRyNOhT4xMxmAxOBs939pzhizGgPPQTXXx+ufr7hhrijEZEUFNfZSpOASUnefwZ4pvojqkGmTYORI0Nvbrr6WUSKkFLVSlLFPvsMjj0W2reHp5/W1c8iUiQlh5pi6VI4+mjYbrvQCN24cdwRiUgKi+sKaalO2dnwxz/CTz/Bm29Cq1ZxRyQiKU7JIdPl5oYuPmfPhilToEuXuCMSkTSg5JDpLr44VCPdcw/84Q9xRyMiaUJtDpnsnntCv89/+Uvon0FEpJSUHDLVSy/BBReEtoZRo+KORkTSjJJDJpozB044AfbbD554AmrXjjsiEUkzSg6Z5ocfQveeDRqEDnsy5M6TIlK91CCdQWpt3gzHHAPLlsEbb+iUVREpNyWHTOFO29tug3feCVc/H3BA3BGJSBpTtVKmuPVWfjN1ariR3nHHxR2NiKQ5JYdM8NxzcOWVLD/ssHALbhGRClJySHeffQYnnQRdu/LFpZfqLqsiUimUHNLZihWhB7dGjeDZZ9lSv37cEYlIhlBySFebN4e2haVLYfJkaNEi7ohEJIPobKV0deGF8Prr8NhjcOCBcUcjIhlGRw7p6L77YMwYuPxyOPnkuKMRkQyk5JBu3n4b/vxnOOoouOmmuKMRkQyl5JBOFi8O3Xy2bq17JolIlVKbQ7rYuBEGD4Z16+C116BJk7gjEpEMpuSQDtzh7LPhgw/CmUkdOsQdkYhkOFUrpYO77oKHH4brr4eBA+OORkRqACWHVPf666Grz4ED4Zpr4o5GRGoIJYdUtmQJHH887LUXPPII1NLiEpHqoTaHVLVpU7gCev16yMqC7bePOyIRqUGUHFLVRRfBjBkwcSK0bx93NCJSw8RST2Fmo8zsczP7xMwmmVmThHFXmtlXZvaFmR0ZR3yxGzv2lyugjz027mhEpAaKqxJ7KtDR3TsBXwJXAphZB2AIsC/QD7jXzGrWlV4zZ8K550LfvnDjjXFHIyI1VCzJwd1fcffc6OUMoGX0fCAw3t03ufs3wFdAzbmr3IoV4UiheXN48kmoo1o/EYmHuXu8AZg9B0xw98fM7G5ghrs/Fo17EPifu09M8rmRwEiA5s2bdxs/fny5Y8jOzqZhw4bl/nylyMuj05VX0mT2bD666y7WtmtX5q9IiXJUEpUl9WRKOUBlydenT58P3b170pHuXiUDMA34LMkwMGGaq4FJ/JKk7gZOSRj/IHBcSb/VrVs3r4jp06dX6POV4vrr3cH9vvvK/RUpUY5KorKknkwph7vKkg+Y6UVsV6us3sLd+xY33syGAf2Bw6MgAZYArRImaxm9l9leegn+/nc47TQ488y4oxERie1spX7AZcAAd1+fMGoKMMTM6pnZHsDewPtxxFhtFi0KfTJ07Aj//rf6gBaRlBBXi+fdQD1gqoWN4Qx3P9vd55jZU8BcIBc4z93zYoqx6m3aBH/6E+TkhOsZGjSIOyIRESCm5ODuexUz7iagZvRic/HF8P778Mwz0LZt3NGIiGylm/XE5ckn4d574a9/Df00iIikECWHOHzxBYwcCT16wD//GXc0IiKFKDlUt/XrQztDvXowYQLUrRt3RCIihegS3Op2wQXw6afw4ovQsmXJ04uIxEBHDtXp0UfhwQfhqqvgqKPijkZEpEhKDtVl7tzQD/TvfhcueBMRSWFKDtVh3brQzrDddvDEE7qhnoikPG2lqsMFF8C8efDyy7DrrnFHIyJSIh05VLUnngid91x9Nfz+93FHIyJSKkoOVWn+fDjrLOjZE667Lu5oRERKTcmhqmzaBEOGwDbbqJ1BRNKOtlhV5bLLYNYsePZZaNWq5OlFRFKIjhyqwrPPwujRcOGFMGBA3NGIiJSZkkNl++47GD4cunaFW26JOxoRkXJRcqhMeXmh456cnHDfpHr14o5IRKRc1OZQmW66Cd58M9wmY68iu6wQEUl5OnKoLG+9FW6LccopYRARSWNKDpVh1apQndS6NdxzT9zRiIhUmKqVKso9XOi2dCm8/TZsv33cEYmIVJiSQ0WNHQtPPx16dDvwwLijERGpFKpWqojPPw831TvssHDRm4hIhlByKK/Nm+Gkk2DbbcPZSbU0K0Ukc6haqbyuuQY++ggmTdJtuEUk42h3tzyysmDUKDjzTBg0KO5oREQqnZJDWa1aBaeeGi5yu+OOuKMREakSqlYqC/fQD/SyZfDOO6HbTxGRDBTLkYOZjTKzz83sEzObZGZNovdbm9kGM5sdDWPiiK9Ijz4KTz0VroQ+4IC4oxERqTJxVStNBTq6eyfgS+DKhHEL3L1zNJwdT3hJfPMNnH8+9OoFl18edzQiIlUqluTg7q+4e270cgbQMo44Si03N9wvySwcPdSuHXdEIiJVKhUapEcA/0t4vYeZfWRmr5tZr7iC+pVbbw1tDPfcA7vvHnc0IiJVzty9ar7YbBrwmySjrnb3Z6Nprga6A4Pd3c2sHtDQ3VeaWTdgMrCvu69J8v0jgZEAzZs37zZ+/Phyx5qdnU3Dhg2Tjms4fz5dzzmHFb16Mffaa8PRQ4oqrhzpRmVJPZlSDlBZ8vXp0+dDd++edKS7xzIAw4B3gQbFTJMFdC/pu7p16+YVMX369OQj1q9379DBfZdd3FeurNBvVIciy5GGVJbUkynlcFdZ8gEzvYjtaiynsppZP+Ay4Hfuvj7h/Z2An9w9z8zaAHsDX8cRIwBXXQVz58JLL8EOO8QWhohIdYvrOoe7gXrAVAvVNDM8nJl0KPAPM8sBtgBnu/tPsUT46qtw551w3nlw5JGxhCAiEpdYkoO7J+1D092fAZ6p5nAK+/lnGDYM2rYNjdEiIjWMrpBO5vzz4fvvwxlKDRrEHY2ISLVLhVNZU8vEifD44/C3v6nzHhGpsZQcEi1fDuecA926wdVXxx2NiEhslBzy5d9Ub+1aePhhqFs37ohERGKjNod8jz0GkyeHfhr23TfuaEREYqUjB6Dejz/Cn/8MhxwCf/lL3OGIiMROycGddrfeCjk5MG6cbqonIkINTw6LVsOdI+5jh5kzGX3cKBbtlPTyCxGRGqdGJ4er71/AmU9cwqJO3bmj6zmMmBJ3RCIiqaFGN0h/twZm7nEIX591OltyjK9XxR2RiEhqqNFHDnX23pNh575M9o47U8ugTdO4IxIRSQ01OjmMHQB7NgUjPI4dEHdEIiKpoUYnh90aw7RToePO4XG3xnFHJCKSGmp0chARkeSUHEREpBAlBxERKUTJQUREClFyEBGRQpQcRESkEHP3uGOoMDP7EVhYga/YEVhRSeHEKVPKASpLKsqUcoDKkm93d98p2YiMSA4VZWYz3b173HFUVKaUA1SWVJQp5QCVpTRUrSQiIoUoOYiISCFKDsF/4g6gkmRKOUBlSUWZUg5QWUqkNgcRESlERw4iIlKIkoOIiBRSI5ODmY0ys8/N7BMzm2RmTYqYrp+ZfWFmX5nZFdUcZonM7E9mNsfMtphZkaeymdm3Zvapmc02s5nVGWNplaEsKb1MAMxsBzObambzo8ek3UiZWV60TGabWcp0UlvSPDazemY2IRr/npm1jiHMUilFWYaZ2Y8Jy+GMOOIsiZmNNbMfzOyzIsabmY2OyvmJmXWt8I+6e40bgCOAOtHzW4BbkkxTG1gAtAG2AT4GOsQde4EY2wPtgCygezHTfQvsGHe8FS1LOiyTKM5bgSui51ckW7+icdlxx1qeeQycC4yJng8BJsQddwXKMgy4O+5YS1GWQ4GuwGdFjD8a+B+h77KDgfcq+ps18sjB3V9x99zo5QygZZLJDgS+cvev3X0zMB4YWF0xloa7z3P3L+KOozKUsiwpv0wiA4GHo+cPA4PiC6XMSjOPE8s3ETjczKwaYyytdFlfSuTubwA/FTPJQOARD2YATcxsl4r8Zo1MDgWMIGTcgloA3yW8Xhy9l44ceMXMPjSzkXEHUwHpskyau/v30fNlQPMipqtvZjPNbIaZDaqe0EpUmnm8dZpoJ2s10Kxaoiub0q4vx0ZVMRPNrFX1hFbpKv2/UadC4aQwM5sG/CbJqKvd/dlomquBXODx6oytLEpTjlLo6e5LzGxnYKqZfR7tiVSrSipLSiiuLIkv3N3NrKjzxXePlksb4DUz+9TdF1R2rFKs54An3X2TmZ1FOCI6LOaYUkLGJgd371vceDMbBvQHDveo0q6AJUDiXkTL6L1qVVI5SvkdS6LHH8xsEuFwu9qTQyWUJSWWCRRfFjNbbma7uPv30aH9D0V8R/5y+drMsoAuhDryOJVmHudPs9jM6gCNgZXVE16ZlFgWd0+M+wFCe1E6qvT/Ro2sVjKzfsBlwAB3X1/EZB8Ae5vZHma2DaHhLWXOKCktM9vOzBrlPyc0xic94yENpMsymQIMjZ4PBQodFZlZUzOrFz3fETgEmFttERatNPM4sXzHAa8VsYMVtxLLUqBefgAwrxrjq0xTgNOis5YOBlYnVG2WT9yt8HEMwFeE+rnZ0ZB/5sWuwIsJ0x0NfEnYm7s67riTlOMYQt3iJmA58HLBchDO1Pg4GuakYjlKW5Z0WCZRjM2AV4H5wDRgh+j97sAD0fMewKfRcvkUOD3uuIubx8A/CDtTAPWBp6P/0ftAm7hjrkBZ/hn9Lz4GpgP7xB1zEeV4EvgeyIn+J6cDZwNnR+MNuCcq56cUc/ZiaQfdPkNERAqpkdVKIiJSPCUHEREpRMlBREQKUXIQEZFClBxERKQQJQeRKmJmTczs3LjjECkPJQeRqtOEcAdTkbSj5CBSdW4G9oz6CRgVdzAiZaGL4ESqSNQJzvPu3jHuWETKSkcOIiJSiJKDiIgUouQgUnXWAo3iDkKkPJQcRKqIh74C3jazz9QgLelGDdIiIlKIjhxERKQQJQcRESlEyUFERApRchARkUKUHEREpBAlBxERKUTJQURECvn/FIR2+9eMnvEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "%matplotlib inline\n",
    "\n",
    "# Fit the model to a dense 1D-grid of points from [t_min, t_max].\n",
    "fitn = m*20\n",
    "fitx = np.linspace(np.min(t), np.max(t), fitn)  # Actually t.\n",
    "fity = cs(fitx)\n",
    "\n",
    "# Plot the fitted model against observations.\n",
    "plt.scatter(t, y, s=16, c='dodgerblue', label='observations')\n",
    "plt.plot(fitx, fity, c='red', label='fit')\n",
    "plt.title('Observations vs. Fit, Cubic Spline Interpolation')\n",
    "plt.ylabel('y')\n",
    "plt.xlabel('t')\n",
    "plt.grid()\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 07.08 B-splines"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "B-splines are basis functions for spline functions of any degree.\n",
    "* Can be defined in various ways: recursion, convolution, and divided differences.\n",
    "* Linear function used in recursive definition.\n",
    "$$\n",
    "v_i^k(t) = \\frac{t - t_i}{t_{i+k} - t_i}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Summary: Interpolation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* Interpolating functions are **exact** fit at sample points.\n",
    "* Interpolating functions are linear combinations of basis functions with coefficients determined depending on choice of technique.\n",
    "* Polynomial interpolation.\n",
    "  * Coefficients for Lagrange interpolation are cheaper to compute and better conditioned than monomial interpolation, but cost of evaluation for Lagrange interpolation is a little higher $5n$ vs. $2n$.\n",
    "  * Coefficients for Newton interpolation can be constructed incrementally.\n",
    "  * Convergence depends on location of sample points.\n",
    "    * Equally spaced sample points are expected to have inflated errors at the extrema.\n",
    "    * Sample points spaced at Chebyshev points distribute the error more evenly.\n",
    "* Piecewise polynomial interpolation.\n",
    "  * Used to fit large number of samples with low-degree polynomials.\n",
    "  * Cubic spline is preferred when smoothness is important."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}