{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# EECS 16A Discussion 12B" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Problem 1: Polynomial Fitting" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this discussion, we are trying to fit data (observations) of the form $\\{(x_i,y_i),i=1,2,...,n\\}$ to a polynomial that we know looks like this:\n", "\n", "$$y = f(x) = a_0 + a_1x + a_2x^2 + a_3x^3 + a_4x^4$$\n", "\n", "In other words, we want to find $a_0$, $a_1$, $a_2$, $a_3$, and $a_4$ that best fit the data.\n", "\n", "More generally, we might want to fit the data to a polynomial of different degree (for instance, if we do not know that the polynomial looks like as above), so we could try to solve for some $a_0,a_1,\\ldots,a_d$ that define a $d$-degree polynomial.\n", "\n", "Note that the observations are not perfect -- they are *noisy*, which means that $y_i \\neq f(x_i)$ in general! That is what makes this problem interesting.\n", "\n", "This first block of code contains functions that will help us set up some useful objects -- the polynomial curve for a set of parameters $a_0$, $a_1$, $a_2$, $a_3$, and $a_4$ and a \"data\" matrix that we will use to compute the least squares estimate later." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Initialization" ] }, { "cell_type": "code", "execution_count": 55, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import numpy.matlib\n", "import matplotlib.pyplot as plt\n", "from mpl_toolkits.mplot3d import Axes3D\n", "from numpy import pi, cos, exp\n", "\n", "%matplotlib inline" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Construct Data Matrix\n", "Use this block to make the data matrix. The input data is given in the array x_data, and the output data is given in the array y_data.\n", "\n", "To construct the data matrix, which we call DataMat, you only need to specify what degree polynomial you will use to fit the data.\n", "\n", "If x_data has the form x_{data}$=\\begin{bmatrix}x_1\\\\x_2\\\\ \\vdots \\\\x_n \\end{bmatrix}$, then DataMat has the form DataMat $=\\begin{bmatrix}1 &x_1&x_1^2& \\cdots &x_1^d\\\\1 &x_2&x_2^2& \\cdots &x_2^d\\\\ \\vdots & \\vdots & \\vdots & \\ddots & \\vdots \\\\1 &x_n&x_n^2& \\cdots &x_n^d \\end{bmatrix}$, where $d$ is the degree of the polynomial. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Function that defines a data matrix for some input data" ] }, { "cell_type": "code", "execution_count": 56, "metadata": {}, "outputs": [], "source": [ "\"\"\"You'll understand why it's constructed like this after\n", "doing the worksheet!\"\"\"\n", "def data_matrix(input_data,degree): \n", " # degree is the degree of the polynomial you plan to fit the data with \n", " Data=np.zeros((len(input_data),degree+1))\n", " \n", " for k in range(0,degree+1):\n", " Data[:,k]=(list(map(lambda x:x**k ,input_data)))\n", " \n", " return Data\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Function that helps with plotting " ] }, { "cell_type": "code", "execution_count": 57, "metadata": {}, "outputs": [], "source": [ "\"\"\"Function that defines the polynomial curve for a set of parameters and a range. The set of parameters defines the degree of the\n", "polynomial.\"\"\"\n", "def poly_curve(params,x_input):\n", " # params contains the coefficients that multiply the polynomial terms, in degree of lowest degree to highest degree\n", " degree=len(params)-1\n", " x_range=[np.min(x_input), np.max(x_input)]\n", " x=np.linspace(x_range[0],x_range[1],1000)\n", " y=x*0\n", " \n", " for k in range(0,degree+1):\n", " coeff=params[k]\n", " y=y+list(map(lambda z:coeff*z**k,x)) \n", " return x,y\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Printing the data matrix" ] }, { "cell_type": "code", "execution_count": 58, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[1.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00]\n", " [1.000000e+00 5.000000e-01 2.500000e-01 1.250000e-01 6.250000e-02]\n", " [1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00 1.000000e+00]\n", " [1.000000e+00 1.500000e+00 2.250000e+00 3.375000e+00 5.062500e+00]\n", " [1.000000e+00 2.000000e+00 4.000000e+00 8.000000e+00 1.600000e+01]\n", " [1.000000e+00 2.500000e+00 6.250000e+00 1.562500e+01 3.906250e+01]\n", " [1.000000e+00 3.000000e+00 9.000000e+00 2.700000e+01 8.100000e+01]\n", " [1.000000e+00 3.500000e+00 1.225000e+01 4.287500e+01 1.500625e+02]\n", " [1.000000e+00 4.000000e+00 1.600000e+01 6.400000e+01 2.560000e+02]\n", " [1.000000e+00 4.500000e+00 2.025000e+01 9.112500e+01 4.100625e+02]]\n" ] } ], "source": [ "x_data = [0.0, 0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4, 4.5]\n", "y_data = [24.0, 6.61, 0, -0.95, 0.07, 0.73, -0.12, -0.83, -0.04, 6.42]\n", "\n", "D = data_matrix(x_data, 4)\n", "print(D)\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Least Squares\n", "Recall the system of equations we are trying to solve here:\n", "\n", "$$\\texttt{D}\\vec{a} = \\vec{y},$$ where $\\vec{a} = \\begin{bmatrix}a_0\\\\a_1\\\\:\\\\a_d\\end{bmatrix}$ and $\\vec{y} = \\begin{bmatrix}y_1\\\\y_2\\\\ : \\\\y_n \\end{bmatrix}$.\n", "\n", "In the next block, you will implement code to compute $\\vec{a}$, the set of optimal polynomial coefficients (optimal in a least squares sense) to fit the data. Put the coefficients in an array called params." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Finding the least squares solution" ] }, { "cell_type": "code", "execution_count": 59, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[ 24.00958042 -49.99515152 35.0039627 -9.99561772 0.99841492]\n" ] } ], "source": [ "# params = np.linalg.solve(np.dot(np.transpose(D), D), np.dot(np.transpose(D), y_data))\n", "params = np.linalg.inv(np.transpose(D) @ D) @ np.transpose(D) @ y_data\n", "print(params)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Plot Curve Fit\n", "Use the next block to compare your fitted polynomial to the data. All you need to do is enter the polynomial coefficients (in degree of lowest degree to highest degree) into an array called params. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Plotting the fitted polynomial and the original datapoints" ] }, { "cell_type": "code", "execution_count": 60, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Text(0.5, 1.0, 'Polynomial of Degree 4')" ] }, "execution_count": 60, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXYAAAEKCAYAAAAGvn7fAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAsw0lEQVR4nO3deXxU5dn/8c+VEBLCoogRomChsq9BogRRlIIW0F8VKxVc6lrrjra4PxV9tFYLijx9rNZqRVGsuJbHlUVcQKKGHWQJSpQgyCL7Tub6/XGfwBATMklm5sycXO/Xa16TmTlnznVOku/cc5/7nCOqijHGmOBI8bsAY4wx0WXBbowxAWPBbowxAWPBbowxAWPBbowxAWPBbowxAWPBnuBE5HIR0bDbNhGZLyI3ikidKr7XRyLyUYxKjQsRaelth8urMW9U119E/p+ILBSR3V5NR1Yw3X3l/A4LRWSCiPwyWvUkAxH5uYjs9LZDa7/rCaoqBYPx1RCgGGjk/fw34BjgXj+L8sEaoBfwtZ9FeB+qLwGfATcAe4Ftlcx2KlACZAKtgAuA90XkReAyVQ3FruKE8XdgC1DP70KCzII9ecxT1RXez5O91s5walmwq+oeIN/vOoDjgIbARFX9JMJ5PlfV/WGPnxWRW4HHgHnAo9EtsWIikgbs1zgeoSgiFwHdgb8AY+K13NrIumKS15dAIxE5BkBEBojILBHZJSJbROQtEWlX0cwi0kxE9orI8HJeu8/7utzYe/yRiMwQkf4iMsd7bZGIDC5n3krrCHu/ASIyz5t2roj0FJE6IvKQiKwRkR9FZJyI1A+b9yddMSJykoi8JiLF3nst896jWq1CEckWkRdEZIOI7BGRBSJySfj2AYq8h8969XxUnWWp6hhgLu5DOryGLBF5SkRWezUsFZFryqm1v7ftdovIChG52ttmRWHTlG6z60XkryLyPbAHONJ7/XwRyfd+r5tF5FUROb6cZV3jdQPu9rbNsyJyVCTr6f0tPQaMADZHun1M9ViwJ69WuK/120VkAPAOsB24ELgO6AzMEJHjyptZVdcCbwGHhIWIpAJX4Vqim8JeOgEYi/vnPB/XJfJqeD9pFetoDYwCHsZ1LaUDk4AngWzgcuC/gYuBkZVsi+NxLd5rgQFenVcCz1Uy3094HyIfAwOBu4HzgIXA+LBgfcarGeBBXNfQ9VVdVpj3gBalYSoijYAZwCDgPuBs4P+AJ0XkprBaO3Jwew/16h0O/KKC5dwDtMX9zgcDu0XkWuB14Ctc19Dvcb+zj0WkYdiyHgaeAKYCvwJuw23r97y/mcr8FViqquMjmNbUlKraLYFvuIBToB2u66wx7p+vBHjLm6YAKATqhM3XCtgHPBb23EfAR2GPz/De+7Sw537lPZdXZr59QJuw547xarg77Lmq1LEP+Hk5y51aZv3fAFaGPW7pTXd5BdtLvO10CRACmlS0/hXMf6P3/meUeX4qsA5I9R63PlwdZea9z5u2TgWv/957vaf3+E/A7vDt7T3/T2BD6fsAE4D1QGbYNNnevEXlbLM5gIQ93wDX3/2vMstphdtncEvY/CXAvWWm6+2973mVrP9puG8IHcv8Tbf2+/8rqDdrsSePpbgw/BG3A+ol4EqvhXki8IqG9d+q6kpgJnB6RW+oqh/hWmq/D3v698ACVS3bj12oqoVh867DBV1pK7OqdSxX1W/KrB/AB+Wsd3MRkYrWQ0QaicgjIvI1LkD2AeNxId+movkq0AdY7W2bcC8CWUDHKr5fJErXrbS/ewDwObDS65qq4+2s/QBoElZDHvCuqu4sfSNVXYPboVuet9RLVk8v3M74l8osZxVuu/fxpjsT9+2+7HSf43YY96ECIlIX+AcwRlW/qnRLmKiwnafJYzBuVMw24FtV3Q0gIs1xwbCmnHnWAj+r5H2fBEZ7fe0NcKFyYznT/VjOc3uADO/nxlWsY1OZx3sP83wdIBXYT/meA/rjdiTPA3YAJ+O6DjIqmKciR1HxOpS+Hm0tvPvS5R6D+0awr4Lpm3j32bgP17J+AH5ezvNl1+sY735qBcvZVGa6FRVM16SC5wFuwf1t/I8cHA6a6d03FJGGqlrZaCJTRRbsyWORHhwVE24TrqXXrJzXmlF+IId7ATdK4XLcP+BO3LeBqqppHdUiIhnAucB9qjo27Pku1XzLH3HdXmU1C3s92gYB36nqKu/xRlxg/2THtmeZd7+Gg6EbrmkF85UdAbPRu78cWFzO9NvKTHcWP/3gDX+9PB1x2251Oa/NAeYDOYeZ31SDBXuSU9UdIjIbGCIi96lqCYCI/Aw4BTfe/XDzbxWRl3BdMA2Al1V1a7zrqIF0XGu+bOv28mq+38e4deitqjPDnr8IF7ZR7U4QN9wxB9eyLfU+cBMu7MtrkZfKBwaJSGZpd4yIZOP6vsv71lHWZ7jwbq2qzx9muim4/RXHq+qUCN433MPAuDLPDQDuwO0HWVZ2BlNzFuzB8Cfc6Ii3ReTvuIC+H7djLJKx0X/nYD/7Uz7WUWWqukVE8oE/isga3M7FK3HjzKtjHK6l/IaI3IPr/roY18/8+9IPrGrqKSIluO6hn+NGoQwEngf+J2y6MbhRRZ+KyBhc+NUH2uN2dJ/rTfeg9x4fiMho3Ifcn3BdMZUe7OR9qN8GPCEiWbjROVtw2+503I7mCar6tYg8AvyvN3T1Y9wO2ha47fKMqk6vYBlLObj/BHDDL70fP6/gW6ipIQv2AFDV90XkbNywwIm4fumPgNtV9fsI5l8gIsuBrao6x686amAYbl/BE8Aub9nDgber+kbeN4/TccPzHsYdhLQMuFRVX6xhnTO8+524FvUXwABVPWSHsfdhdQpun8EduKDd7NXxeth0X3nbexRunVcDj+BaxC0jKUhV/yEiq3DDFy/CZcJq4FPc/orS6e4WkSW4o2xvwHXrrAKm4UZCmQQih+4kN7WR1wpbAvxOVZ/1ux5TfSLSALeT8x1Vvcrveow/LNhrMW9ETWtcd0lrXF/rLn+rMlUhIn/D9ZV/DxyL+6bSHThJVRf4WZvxT8Tj2EWkhYhMF5GvRGSxNzyu9PDz1eIODZ8nIoNiV66JsquBD3GjKC6yUE9KGbjul8nA07ihnv0t1Gu3iFvs3t72bFWd4x1qPBt3uPVvgO2qOjpmVRpjjIlYxDtPvSPa1ng/b/N2pFR35IExxpgYqVYfuzdc6RPcyYL+gBszvBV3rpA/6qEnjyqd5xq8E07Vr1+/R/v27atdtDHG1EazZ8/eoKpZlU1X5WD39rp/DPxZVd8Qkaa4scMKPIDrrrnycO+Rm5urBQUFVVquMcbUdiIyW1VzK5uuSicBE3dy/teBl1T1DQBV/UFVS9Rd/eWfuHN0GGOM8UlVRsUI8CywRFUfC3s+O2yywcCi6JVnjDGmqqpy5Glv4FJgoYjM8567GxgmIjm4rpgiDj0FbK2Snw9jx8Ly5dC2LQwfDnl5fldljKltqjIqZgYHzxsd7t3olZO8Ro6E0aNh1y5QhXnzYNIkGDEC7r/f7+qMqb59+/ZRXFzM7t27/S6l1sjIyKB58+akpaVVa347V0wU5Oe7UN+9s4QuLGITjVkVOp6dO93zAwday90kr+LiYho2bEjLli05zPVOTJSoKhs3bqS4uJhWrVpV6z3sCkpRMHasa6mnsY85nMhVHDzdyu7d7nVjktXu3btp0qSJhXqciAhNmjSp0TckC/YoWL7cdb/sIYOVtKIDSw68FgpBoZ37ziQ5C/X4qun2tmCPgrZtIcXbkkvocEiwp6S4140xJl4s2KNg+HDI8K6suYQOtGU5qd7lOTMy4OabfSzOmDjLz4dhw6BHD3efX/ay6NWQmppKTk4OnTp1olu3bjz66KOEQoe/lkhRURETJkyo+cIrcfXVV/PVV4e/sNZbb71V6TTRZMEeBXl5bvRLZiYskw6ks5cTZCWZme5523FqaouRI6FfP3jlFZgzByZOdI9HjqzZ+9arV4958+axePFipkyZwnvvvcf9lQw3i1ewP/PMM3Ts2PGw08Q72FHVuN969OihQTRrlup/9Z+lCjrqtP/orFl+V2RMzX311VcRTTdrlmpmpqrb43ToLTNTa/T/UL9+/UMef/3113rUUUdpKBTSlStX6qmnnqrdu3fX7t2768yZM1VVtWfPntqoUSPt1q2bPvbYYxVOF27lypXarl07veiii7R9+/b661//Wnfs2KGqqlOnTtWcnBzt3LmzXnHFFbp7925VVT399NP1yy+/PFDn3XffrV27dtWePXvq2rVrdebMmdq4cWNt2bKlduvWTVesWKFjx47VDh06aJcuXfTCCy8sd53L2+5AgUaQsRbs0bZ5s9usDz/sdyXGREWkwT50qKpI+cGekuJer66ywa6qesQRR+jatWt1x44dumvXLlVVXb58uZbmy/Tp0/Xss88+MH1F04VbuXKlAjpjxgxVVb3iiit01KhRumvXLm3evLkuW7ZMVVUvvfRSHTNmjKoeGuyATpo0SVVVb7vtNn3ggQdUVfWyyy7TV1999cBysrOzD3wwbNq0qdx1rkmwW1dMtB1xBGRnw9KllU9rTICUjg4rTyxHh+3bt4/f/e53dOnShSFDhlTY5RHpdC1atKB3794AXHLJJcyYMYNly5bRqlUr2nojIS677DI++eSTn8xbt25dzjnnHAB69OhBUVFRucvo2rUrF198MS+++CJ16kT/cCIL9ljo0AGWLKl8OmMCJHx0WFnRHh32zTffkJqayjHHHMOYMWNo2rQp8+fPp6CggL1795Y7T6TTlR1qWJWhh2lpaQemT01NZf/+/eVO984773DDDTcwZ84cTjrppAqnqy4L9lgoDXa7nqypRcJHh5UVzdFh69ev59prr+XGG29ERNiyZQvZ2dmkpKQwfvx4SkpKAGjYsCHbtm07MF9F05X13XffMWvWLAAmTJjAqaeeSrt27SgqKmLFihUAjB8/ntNPPz3imsNrCYVCrFq1ir59+/LII4+wZcsWtm/fXq1tUREL9ljo0AG2boU1a/yuxJi4CR8dVtpyT0khKqPDdu3adWC4Y//+/TnrrLMY6Q21uf7663n++efp1q0bS5cupX79+oDr7khNTaVbt26MGTOmwunKateuHU888QQdOnRg06ZNXHfddWRkZPDcc88xZMgQunTpQkpKCtdee23E9Q8dOpRRo0bRvXt3CgsLueSSS+jSpQvdu3fn5ptv5sgjj6z+xilHta6gVFOBv9DGtGnQvz9MnerGehmTxJYsWUKHDh0inr70LKeFhdCmTXKd5bSoqIhzzjmHRYv8P/t4eds90gtt2EnAYqH0l7FkiQW7qXXy8pInyIPKumJiITsbGjWyHajGJJmWLVsmRGu9pizYY0HERsYYY3xjwR4rFuzGGJ9YsMdKhw6wdi1s3ux3JcaYWsaCPVbCd6AaY0wcWbDHigW7MVHToEEDAL7//nsuuOCCqLzn0qVLycnJoXv37nz99deccsopQPzOChlLFuyx0qoVpKdDPE/VaUzAHXvssbz22mtRea+33nqLCy64gLlz53LCCSfw2WefARbs5nBSU6FjRwjA0CljEkVRURGdO3cGYNy4cZx//vkMGDCANm3acPvttx+YbvLkyfTq1YsTTzyRIUOG/OSQ/XfffZfHH3+cJ598kr59+wIHvxXceeedfPrpp+Tk5DBmzJg4rVl02QFKsdS5M3z4od9VGBM9t9wC8+ZF9z1zcuDxx6s167x585g7dy7p6em0a9eOm266iXr16vHggw8ydepU6tevzyOPPMJjjz3Gvffee2C+QYMGce2119KgQQNGjBhxyHs+/PDDjB49mrfffrsGK+UvC/ZY6twZxo+HTZugcWO/qzEmcPr168cRRxwBQMeOHfn222/ZvHkzX3311YFT7+7du5devXr5WWbcWbDHkveVkUWL4LTT/K3FmGioZss6VtLT0w/8XHqaXFXlzDPP5OWXX/axMn9ZH3ssdeni7q2f3Zi4ycvLY+bMmQdOsbtjxw6WL18e8fxlT/ebjCzYY6l5c3fOGAt2Y+ImKyuLcePGMWzYMLp27UqvXr1YWoUrmpU93W8yivi0vSLSAngBaAoo8LSqjhWRo4BXgJZAEfAbVd10uPcK/Gl7w/XuDXXqwMcf+12JMdVS1dP2muioyWl7q9Ji3w/8UVU7AnnADSLSEbgTmKaqbYBp3mNTqnNn12K3qykZY+Ik4mBX1TWqOsf7eRuwBDgOOBd43pvseeC8KNeY3Lp0gR9/tKspGWPiplp97CLSEugOfA40VdXS1FqL66oxpcJHxhiTpPy40lptVtPtXeVgF5EGwOvALaq6tUwxiut/L2++a0SkQEQK1q9fX61ik1KnTu7egt0kqYyMDDZu3GjhHieqysaNG8mo6MrgEajSOHYRScOF+kuq+ob39A8ikq2qa0QkG1hXQbFPA0+D23la7YqTTVYWNG1qwW6SVvPmzSkuLqZWNch8lpGRQfPmzas9f8TBLiICPAssUdXHwl6aBFwGPOzd/6fa1QRVly6wcKHfVRhTLWlpabRq1crvMkwVVKUrpjdwKfALEZnn3QbhAv1MESkE+nuPTbjOnWHxYgiF/K7EGFMLRNxiV9UZgFTwcr/olBNQnTvDrl2wciWccILf1RhjAs6OPI0HGxljjIkjC/Z4KB0ZY/3sxpg4sGCPhwYNXBfM/Pl+V2KMqQUs2OMlJyf6FygwxphyWLDHS04OrFgBSX46UGNM4rNgj5ecHHe/YIGvZRhjgs+CPV5Kg926Y4wxMWbBHi/HHQdNmliwG2NizoI9XkRsB6oxJi4s2OMpJ8eNZd+/3+9KjDEBZsEeTzk5sGcPLFvmdyXGmACzYI8n24FqjIkDC/Z4atcO0tMt2I0xMWXBHk9pae6EYBbsxpgYsmCPt9KRMXaZMWNMjFiwx1tODmzYAN9/73clxpiAsmCPN9uBaoyJMQv2eOva1d3PnetvHcaYwLJgj7dGjaB1a5gzx+9KjDEBZcHuh5NOgoICv6swxgSUBbsfcnNh1Sr44Qe/KzHGBJAFux9yc9397Nn+1mGMCSQLdj907+7O9mjdMcaYGLBg90PDhtC+vQW7MSYmLNj9kptrwW6MiQkLdr+cdBKsWWNHoBpjos6C3S+lO1Ct1W6MiTILdr906wapqfDll35XYowJmIiDXUT+JSLrRGRR2HP3ichqEZnn3QbFpswAysyETp2sxW6MibqqtNjHAQPKeX6MquZ4t3ejU1YtUboD1U7ha4yJooiDXVU/AX6MYS21T26uO4Xvd9/5XYkxJkCi0cd+o4gs8LpqGlc0kYhcIyIFIlKwfv36KCw2AGwHqjEmBmoa7E8CJwA5wBrg0YomVNWnVTVXVXOzsrJquNiA6NrVXS7PdqAaY6KoRsGuqj+oaomqhoB/AidHp6xaIj3dXXgjP9/vSowxAVKjYBeR7LCHg4FFFU1rKtCrl2ux79/vdyXGmICoynDHl4FZQDsRKRaRq4C/ishCEVkA9AVujVGdwZWXBzt3wsKFfldijAmIOpFOqKrDynn62SjWUjv16uXuZ81yZ300xpgasiNP/fazn0GzZtbPboyJGgt2v4m47phZs/yuxBgTEBbsiaBXL1ixAmx8vzEmCizYE0Fenrv//HN/6zDGBIIFeyLIzXVnerTuGGNMFFiwJ4LMTHcaXwt2Y0wUWLAnil694IsvoKTE70qMMUnOgj1R9OoFO3bAIjt41xhTMxbsiaL0QKWZM/2twxiT9CzYE0WrVnDssfDpp35XYoxJchbsiUIETjvNBbtdUckYUwMW7ImkTx9YvRpWrvS7EmNMErNgTySnneburTvGGFMDFuyJpFMnaNzYgt0YUyMW7IkkJQVOPRU++cTvSowxScyCPdH06QOFhbB2rd+VGGOSlAV7ointZ58xw986jDFJy4I90Zx4ojt3jHXHGGOqyYI90aSluaNQbQeqMaaaLNgT0Wmnwfz5sHmz35UYY5KQBXsi6tvXHX368cd+V2KMSUIW7ImoZ0+oVw+mTfO7EmNMErJgT0Tp6a475sMP/a7EGJOELNgTVb9+sHixjWc3xlSZBXui6tfP3Vur3RhTRRbsiSonx503xvrZjTFVZMGeqFJT4YwzXLDb+dmNMVUQcbCLyL9EZJ2ILAp77igRmSIihd5949iUWUv16wfffmvnZzfGVElVWuzjgAFlnrsTmKaqbYBp3mMTLaX97NYdY4ypgoiDXVU/AX4s8/S5wPPez88D50WnLANAu3buOqgW7MaYKqhpH3tTVV3j/bwWaFrRhCJyjYgUiEjB+vXra7jYWkLEtdqnTYOSEr+rMcYkiajtPFVVBSrcy6eqT6tqrqrmZmVlRWuxwTdgAGzYALNn+12JMSZJ1DTYfxCRbADvfl3NSzKHOOss13J/7z2/KzHGJImaBvsk4DLv58uA/9Tw/UxZRx8NJ59swW6MiVhVhju+DMwC2olIsYhcBTwMnCkihUB/77GJtoED4YsvXJeMMcZUoiqjYoaparaqpqlqc1V9VlU3qmo/VW2jqv1VteyoGRMNAwe6g5QmT/a7EmNMErAjT5NBbq7rkrHuGGNMBCzYk0FKituJ+sEHEAr5XY0xJsFZsCeLgQNh/Xob9miMqZQFe7L45S9t2KMxJiIW7MkiK8sNe3z7bb8rMcYkOAv2ZHLeefDll1Bc7HclxpgEZsGeTM49191PmuRvHcaYhGbBnkzat4e2beGtt/yuxBiTwCzYk4mI646ZPh02b/a7GmNMgrJgTzbnnQf799voGGNMhSzYk03PntC0qXXHGGMqZMGebFJS4Fe/gnffhT17/K7GGJOALNiT0Xnnwfbt8OGHfldijElAFuzJqF8/aNQIXn3V70qMMQnIgj0ZpafD4MHwxhvWHWOM+QkL9mQ1dChs2cLoMz+gRw8YNgzy8/0uyhiTCCzYk9T9M/qxgSZkf/oKc+bAxImuh2bkSL8rMyb48vNdYypRG1UW7EkoPx/+OiaN1/k15/If6rGTUAh27oTRoxPvj8yYIBk50jWiXnmFhG1UWbAnobFjYdcu+DdDacAOzuadA6/t3u1eN8ZEX36+azzt3OmuVgkkZKPKgj0JLV/u/qg+oQ9raMaFvHLgtVAICgt9LM6YACttVAG8zvn8gUcPvJZIjSoL9iTUtq07TilEKq8yhHN4m4ZsBdzzbdv6XKAxAVXaqGrFN5zPm6Sx78BridSosmBPQsOHQ0aG+3kCF5HBHi7gNcA9f/PNPhZnTICVNqp+ywuEEF7i4gOvJVKjyoI9CeXlwYgRkJkJX0pPltKOyxlHZqZ7Pi/P7wqNCabhw6FeeojLeJ5p9KOYFgdeS6RGlQV7krr/fpg2DX5zoTDluMvpw6fMGLeC++/3uzJjgisvD/5nyKe0oogX5HLAtdQTrVFVx+8CTPXl5Xl/SKsvhePvofuC52HIA36XZUygXZkyjv2ZDUkZNJgeK6FNG9eST5RQBxAtHbMTR7m5uVpQUBD35QbawIGweDEUFbkmhDEm+rZvh2bN3JHfzzwT98WLyGxVza1sOkuAoLj8cli1ys74aEwsvfEG7Njh/t8SWFSCXUSKRGShiMwTEWuK++Hcc+HII+Ff//K7EmOC65ln4IQToHdvvys5rGi22Puqak4kXxNMDGRkwKWXwuuvw7p1fldjTPAsWgSffgrXXOOuP5zArCsmSK67DvbutVa7MbHw1FNQty5ccYXflVQqWsGuwGQRmS0i15Q3gYhcIyIFIlKwfv36KC3WHKJDBzjjDPcHWFLidzXGBMf27fDCCzBkCGRl+V1NpaIV7Keq6onAQOAGEelTdgJVfVpVc1U1NysJNkzSuv56+PZbeP99vysxJjgmTIBt29y34iQQlWBX1dXe/TrgTeDkaLyvqYbzznPDsf7+d78rMSYYVOHJJ6FLFzjlFL+riUiNg11E6otIw9KfgbOARTV9X1NNaWnwu9/Be+/B11/7XY0xyW/WLJg3z7XWE3ynaalotNibAjNEZD7wBfCOqlo/gJ+uvRbq1IHHH/e7EmOS36OPQuPG8Nvf+l1JxGoc7Kr6jap2826dVPXP0SjM1MCxx8LFF7vRMRs3+l2NMclrxQp4803XWq9f3+9qImbDHYPqj390l3V56im/KzEmeY0Z47o3b7zR70qqxII9qDp3dueP+dvf3KVdjDFVs3EjPPccXHIJZGf7XU2VWLAH2YgR8MMP8OKLfldiTPJ58kl3Hbw//MHvSqrMgj3I+vaFHj3gkUdg/36/qzEmeWzb5gYfnH02dOrkdzVVZsEeZCLwpz+5HUAvveR3NcYkjyeecF0x997rdyXVYudjDzpVyM2FLVtg6VI3DNIYU7Ht26FlSzj5ZHj3Xb+rOYSdj904InDffe5gpfHj/a7GmMRX2lofOdLvSqrNWuy1gSqcdBL8+CMsW+aGbxljfiqBW+tgLXYTTsRd/XrlSvjnP/2uxpjENWqUa63fd5/fldSIBXttMWgQnH66+3q5ZYvf1RiTeFavdsF+4YWuxZ7ELNhrCxF3zosNG+Avf/G7GmMSz733uusYBOD/w4K9NunRw10+7/HHoajI72qMSRzz57ujTG+6CVq18ruaGrNgr23+/GfXer/jDr8rMSYxqMKtt7qLwd9zj9/VRIUFe23TogXceSdMnAgffOB3Ncb4b8IEmD4dHnzQnZ43AGy4Y220Zw907Qr79rkrr2dm+l2RMf7YtAnat3dDHD/7DFJT/a7osGy4o6lYejr84x9u+OODD/pdjTH+ueceN6DgqacSPtSrwoK9tjrjDLjiCje8a/58v6sxJv4++sidwfGmm6B7d7+riSrriqnNNm50521v0gQKCiAjw++KjImPrVtdd2RamrueaZJcHcm6YkzlmjSBceNg8WK46y6/qzEmfm69FVatghdeSJpQrwoL9trul790X0UffxymTPG7GmNi7/XX3fWA77gDevXyu5qYsK4Y464Sk5vrdiLNng3Nm/tdkTGxUVjoDtTr2BE++QTq1vW7oiqxrhgTuXr14LXX3MWvf/1ru0ZqHOTnw7BhLmOGDXOPTYzt3AkXXOD61SdOTLpQrwoLduN06OD6G7/4Am64wR2NZ2Ji5Ejo1w9eeQXmzHEZ069fUp/+O/GFQnD11bBwobua2PHH+11RTFmwm4MGD4b/+i/X//joo35XE0j5+TB6tGs8ln52hkLu8ejR1nKPmZEj4eWX3Sk1Bgzwu5qYs2A3h7r/fvjNb+C22+yKSzEwdqzbpVGe3bvd6ybKnnvOHYh31VXudBq1gF0A0xwqJcV1yWzYAFdeCUcd5a7UbiIXCsH69fDDD+62di2sWwfbtjHwwx2cqjupzw5SKSFECiFSKCGVfaE0Mqc3goeOgCOOcOctyc52O7OPO85O/VAdEye6Lpj+/d3BSCJ+VxQXNirGlG/rVujb151L5tVX4Ve/8ruimMrPd63l5cuhbVsYPhzy8iqZaft212e7YAGsWOFuhYXu+rIV7IDek1qPrSX12Ukm+6lDCiFSKSGFEOns4ciUraSF9pa/vKOOcuc0ad/e7RNp397d2rWzyx2W5803YcgQOOUUeO+9QIxXj3RUTFSCXUQGAGOBVOAZVX34cNNbsCeJTZtcf+ScOTBhAvkthlQ9/JLAyJGuf3vXLtfvnZLiDsIdMcL1TAGuE/yLL9wnwNy57mjFwsKDHeXp6XDCCdC6tbu1bOla202bHrw1bEj+Fyn06+ferqzMTJg2DfJydrurXG3aBN9/D8XFB2/ffANLl8K33x6cMT0dunVzQ2xyc92tY0eoc/ALebU+uJJI2fV7oN2LtP7zFW5bTJ4MDRv6XWJURBrsqGqNbrgw/xr4OVAXmA90PNw8PXr0UJMktmxR7d1bSxC9K22UCiEF1ZQU1cxM1Xvv9bvAmpk1y62HS+iDt6as0WF1X9PvL7xF9aSTVOvUOfhiq1aqgwer3n+/6n/+o7pypWpJScTLvPdet8yUFK3+tty+XXXOHNWXXlIdMUL1jDNUGzU6WGPDhqoDBqg+9JD+84qZekS9PSpSg+UlsNLt6dYvpCNklCro1y37qm7e7Hd5UQUUaCS5HMlEh30D6AV8EPb4LuCuw81jwZ5c8qfv1NdSh6iCPs3Vms6uA/mRmenCMVkNHeoCIZPtOoB39VFu1QV0PhCQe1IzVPv0Ub37btV33lH98ceoLHfWLLfsHj3cfVS2YUmJ6rJlqi++qHrddaqdOh1Yjx3U02n01bv4s55IgQolSf+7Uz30gzmDnfoMV6qC/pvf6JH1dif9+pUVz2C/ANf9Uvr4UuB/DzePBXtyGTpUNYUSfYB7VEHn0k3bsvRA62/oUL8rrIaSEtWCAv3bsQ/pNPrqbuqqgu4iXafQT2/nYe3JLO3ZfY/fldbI785bp4N5Qx/jFp1DzoGgX8fR+jJD9e8nP6daXOx3mdVW+sH8c1bobLqrgj7APSqUJO/f5mFEGuxxGxUjItcA1wAcH/CDA4Jm+XIIkcKfeJBZ9GIclzOHE7mbh3gidAOFhUkyuKq42F01asoUmDoVNm7kRmA+XfkbNzGZs5jBqezCjT5JSYEL2/tbck3N/i6LOQzmTQYDkMU6zmQKZzGZs5hM9hf/hua4s3yedZa79enjjkZOAoXLQlynT/IId7CPNM7h/3iHcwDQkNsNUitFkv6Hu2FdMYE3dOjB/mBQzWa1vsNAVdDZdNc/nZmg33d37lR97z3VW25R7djx4Ao0a6b629+qjh+vBW+vKbePPQjdTKo//d2F31IkpLcPmK86apRq//6q6enuhfR01TPPdM/Pn68aCvm9GuX78ktdcvSpqqDv8UttzneHrl+K6rBhfhcZXcSxK6YO8A3QioM7Tzsdbh4L9uRS/g7GkP6aV3W1HOueOPts1YKCqC5z6FDVE0+sQh90KKS6cKHq6NEumMoG1ejRqgsW/CSoorIzM0FVtHO43A+uHTtU339f9dZbD+mf12bNVC+9VHX8eNW1ayNebpV/f5FavFj14otVQfc2ztJr6z6reDv1g/bBXFbcgt0ti0HActzomHsqm96CPflUFH4P3rFV9aGHVBs3di/06eNGauzaVeNlVTqKY/9+NzLk8cfdKJWjjz74X92xowuo9993gVWJmOzMTBDV/uAqLlZ97jm3QZo0Obhtc3JUb7tNddIk1Q0bKlxeVEfh7NvnfpfnnOPetF491TvvVN2yJdAfzGVFGux2gJKJWOlY4cJCaNOmzFjorVvddSP/8Q831rpRIxg40B3YdMYZcOyxES+j/HHeSqeMb3jjntm03T7HnV74yy/deG+AVq3g9NPdrV8/aNEiSmsdDIf93UUiFHLj9ydPdreZM93F0MEdINW7N/TqxcLUHPre0JGNu356lOyBcfqRLnf7dndq3fffd2dMW7cOjj7aXT/g+uvdz9FavyQR1wOUqsqCPcBCIfjwQ/ePOGmS+2cEF7Q9ehw8gKdZM3fY/BFHuANp9u+H/fv58+1bWPTRBo5mPcexmjYU0poVtGYF9fHSPi0NunRxB5/06eNuFuTxtXOnu5ziZ58dvG3cCEAIYQWtWURnVtCaVbRgFS1YK8fSe2AjHn26oTsKtKQE9u6FPXvcqRfWrIGiooNH886f7z480tPdaS0uuQQGDXKPaykLduO/UMj988+a5W4LFrjW/J49Ec2+lzS+4eesoDWFtOErOrKzQw9emtupVv9zJyRVWLGCEb9cSP2VC+nCQjqziJYUkUFkv+8DsrIOHknbv7/7NpAko3RiLdJgT5JxaiYppaTAySe72/Dh7rlQCFavdifJ2rLF3UpKXKs9NZX/HtOQidOz+EGz+JGjCJF6yNtdmANYpiceEWjThtU92zDx2/MJhUpfULJYz/F8x3Gyhv5527npiu2um6VOHXexi7p14ZhjDp7wLCur1pysK1Ys2E18paS4bpMKuk7OOhoeqeBcKhkZcPPNMa7P1Mjw4a4H7uDvT1jPMaznGJbUg7seAwLY951o7HzsJqHk5bmTb2Vmus8AcPeZme75IO4QCxL7/SUG62M3Cam2jHIIKvv9xYbtPDXGmICJNNitK8YYYwLGgt0YYwLGgt0YYwLGgt0YYwLGgt0YYwLGgt0YYwLGgt0YYwLGgt0YYwLGgt0YYwLGgt0YYwLGgt0YYwLGgt0YYwLGgt0YYwLGgt0YYwLGgt0YYwLGgt0YYwLGgt0YYwLGgt0YYwLGgt0YYwLGgt0YYwLGgt0YYwKmRsEuIveJyGoRmefdBkWrMGOMMdVTJwrvMUZVR0fhfYwxxkSBdcUYY0zARKPFfqOI/BYoAP6oqpvKm0hErgGu8R7uEZFFUVh2ojoa2OB3ETEU5PUL8rqBrV+yaxfJRKKqh59AZCrQrJyX7gHycRtRgQeAbFW9stKFihSoam4kBSYjW7/kFeR1A1u/ZBfp+lXaYlfV/hEu8J/A25FMa4wxJnZqOiomO+zhYCDI3SvGGJMUatrH/lcRycF1xRQBv49wvqdruNxEZ+uXvIK8bmDrl+wiWr9K+9iNMcYkFxvuaIwxAWPBbowxAeNbsIvIEBFZLCIhEQnE8CQRGSAiy0RkhYjc6Xc90SQi/xKRdUE9/kBEWojIdBH5yvu7HO53TdEkIhki8oWIzPfW736/a4o2EUkVkbkiErjReSJSJCILvVO3FFQ2vZ8t9kXA+cAnPtYQNSKSCjwBDAQ6AsNEpKO/VUXVOGCA30XE0H7cAXYdgTzghoD9/vYAv1DVbkAOMEBE8vwtKeqGA0v8LiKG+qpqTiTj2H0LdlVdoqrL/Fp+DJwMrFDVb1R1L/Bv4Fyfa4oaVf0E+NHvOmJFVdeo6hzv5224gDjO36qiR53t3sM07xaYkRMi0hw4G3jG71oSgfWxR89xwKqwx8UEKBhqExFpCXQHPve5lKjyuirmAeuAKaoapPV7HLgdCPlcR6woMFlEZnunZzmsaJwrpkKHOx2Bqv4nlss2pjpEpAHwOnCLqm71u55oUtUSIEdEjgTeFJHOqpr0+0xE5BxgnarOFpEzfC4nVk5V1dUicgwwRUSWet+iyxXTYI/0dAQBsRpoEfa4ufecSRIikoYL9ZdU9Q2/64kVVd0sItNx+0ySPtiB3sCvvOtBZACNRORFVb3E57qiRlVXe/frRORNXNdvhcFuXTHR8yXQRkRaiUhdYCgwyeeaTIRERIBngSWq+pjf9USbiGR5LXVEpB5wJrDU16KiRFXvUtXmqtoS93/3YZBCXUTqi0jD0p+Bs6jkA9nP4Y6DRaQY6AW8IyIf+FVLNKjqfuBG4APcjreJqrrY36qiR0ReBmYB7USkWESu8rumKOsNXAr8IqBXBMsGpovIAlwjZIqqBm5YYEA1BWaIyHzgC+AdVX3/cDPYKQWMMSZgrCvGGGMCxoLdGGMCxoLdGGMCxoLdGGMCxoLdGGMCxoLdGGMCxoLdGGMC5v8DBGVUBKd477wAAAAASUVORK5CYII=\n", "text/plain": [ "