{ "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": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "fig=plt.figure()\n", "ax=fig.add_subplot(111,xlim=[-1,5],ylim=[-5,25])\n", "x,y=poly_curve(params,x_data) ###### Plotting the fitted data\n", "\n", "ax.plot(x_data,y_data,'.b',markersize=15) ###### Plotting the original data\n", "ax.plot(x,y,'r') \n", "ax.legend(['Data points','line fit'])\n", "plt.title('Polynomial of Degree %d' %(len(params)-1),fontsize=16)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## [Practice] Another Polynomial Fitting Problem" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Input and output data\n", "\n", "`x_data` is the input (set of $x_i$'s), and `y_data` is the output (set of $y_i$'s). Here we plot our test points. Does it look like $y$ should be a quartic polynomial in $x$?" ] }, { "cell_type": "code", "execution_count": 61, "metadata": {}, "outputs": [], "source": [ "np.random.seed(10)\n", "\n", "# Parameters corresponding to the polynomial we want: $y = 1600 + 40x -102x^2 -x^3 + x^4$\n", "smarap=np.array([1600,40,-102,-1,1])/100.0\n", "x_data=np.linspace(-11,11,200)\n", "y_data=np.dot(data_matrix(x_data,4),smarap)+(np.random.rand(len(x_data))-0.5)*10\n", "\n", "y_data_training = y_data[0::2]\n", "x_data_training = x_data[0::2]\n", "y_data_test = y_data[1::2]\n", "x_data_test = x_data[1::2]\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Plotting the original datapoints" ] }, { "cell_type": "code", "execution_count": 62, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "200\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAD4CAYAAAAJmJb0AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAxRklEQVR4nO2de5AV9bXvv2vvGWYcwzED8pgIc8FECEyQx0y8cxMjIgOKWseTU4pDqZcTjuGS6A25ddGYIkKIZVWI5lDmmnO8YszxmgjqPTFSiZZB1DIqk5uBoOEhD5WjeBBGRDTh4cis+8ev22n29OPX797d61PVNXvv/u3u3/T+9er1W7/1IGaGIAiCkE9KaXdAEARBiA8R8oIgCDlGhLwgCEKOESEvCIKQY0TIC4Ig5JiatDtg5cwzz+QxY8ak3Q2hgPzHfwAHDgB9ff2flUrAiBHAZz6TXr8EQYdNmza9y8zD7PZlSsiPGTMG3d3daXdDKBhdXcDMmacKeEC9P3IE+NWvgPb2dPomCDoQ0b877QttriGi0UT0LBFtJ6JtRLTY+HwIEa0not3G38aw59KhqwuYNw9obVV/u7qSOKtQbVjHydVXA8eO2bc7fhy4665k+yYIURKFJv8xgP/JzJuJaDCATUS0HsA/ANjAzD8kolsA3ALgOxGcz5Hly4E771Q3LDOwZQuwbh2wZAmwYkWcZxaqicpx4kZfH7B7dzL9EoQ4CK3JM/N+Zt5svP4QwA4AZwG4AsADRrMHAPxd2HO5cd99wO23A0eP9t+4fX3q/Z13ikYvKLq61HiwjhM3SiVg3Lj+78osUag2IvWuIaIxAKYC+AOAEcy839j1DoARUZ4L6L/pmpqAr38dOHnSvp1MuQWTu+5yNs3YUV8PfOtbSvufORN4+GFg82bgkUfU++XL4+urIERBZEKeiD4F4N8AfJuZP7DuY5Ugx1ZvIqKFRNRNRN09PT3a5zNvurVrgXfecW8rU27BZNcufQ2+oUGZ+oCB2r/MEoVqIRIhT0S1UAL+l8z8K+PjA0TUZOxvAnDQ7rvMfC8ztzFz27Bhth5AA7BOuXWwTrmFYjNunBoPdhABzc3KHDN3LrBhg1rLcdP+ZZYoZJ0ovGsIwM8A7GDmf7LsWgdgvvF6PoDHw57LJOiUWxAWL1bjwY7TTlPmmO5uYM2afrdJN+2/rw946imx0wvZJQpN/ssArgNwERFtMbZLAfwQwCwi2g2gw3gfCbpTbgAol9WUW/ycBUCNgyVLlCnG1Oitphm7ceKm/QPA4cNipxeyC2Upn3xbWxvrBEPNm6duqMrglUrKZeCee4Drr4+og0Ju6OpSM8Ldu4FzzlEavpMiYAZL6ZoHGxqUqUcUCyEpiGgTM7fZ7avK3DVuU25A2VYbGoClS0XAC/a0tyuTTKVpxqltpfbvhtjphSxRlULe7qYjUpp7U5OKYDQXzQQhClasUGNq7lxlf290id/u6wP+9Kfk+iYIblSlucbEz5RbEKLEy2RYLquZpCgaQhK4mWuqWsgLQlro2OnFNi8kRe5s8oKQNqbJsFx2biO2eSEsUaTSECEvCAFZsQL43Oec90uktRCGqFJpiJAXhBBMnerscSOR1kJQ7BLpmak0Vq70p9GLkBcEDZymzW7uvBJpLQTFLar/xAl/ruEi5AXBA7dps5MPPREwdmx6fRaqG6+o/h079LV5EfKC4ILbtNnMQGn60E+Y0C/omdWN6GRD1V1Qkxz2xcTLzNfX52NRn5kzs7W2trIgZIlZs5iVyB64lUrMnZ2q3caNzA0N9u0aGtR+k2XL1GdE/cdpaFCfW9FtJ+SPjRv7f3enzSouAXSzg1ytek1eNB0hLpYvB55+2nm/1XtGNx2xzszATzshn7S3AxMnOu/3s6hf1UJeqvUIcWEKWTe7qPVG001HrFs0XHLYC/fdF82iftUKedF0hDjRqVlQWwu8+64S3u++qxZbnTDTEb/5pvvDwJwZeD00xP8+/7S3Azff7C8tth018XUxXnQ0HQknF3QwcyDt2qU088WL9WoW9PUpcw6zEvBhM4RYZwbjxgFbttjnxhH/++KwYgUwZ064HF1VK+RF0xGiYPlyNfM7dkyNpy1bgMceA+rq3L9XKgG9vf3vrWMxqMC3TsEXLwbWrbPPjSP+98WivT2cwlq15hq3aj2i6Qg6OJn8TpwAPvjA+XvlsrMQJwJGj/ZOR2zFbgoepIKVINgRVSHv+4noIBFttXz2fSJ6u6IkYGRIpKHghZfnld9awWYxmuHDnYU8MzBsmCpGcvHF3kVGiJR/vV39g8oc9tbi4oKgS1Sa/L8CuMTm81XMPMXYnojoXABE0xHc0fG88lMrGFAa+oYNwPTperNIrwpmgDr/G2847/dTwUoQ7IhEyDPz8wDei+JYfhBNR7DjvvuA22/39rzyKtBdybBhSsjqziKtiogb4hIpxEncNvkbiegVw5xja6EkooVE1E1E3T09Pb5PIJqOYGX5cmDRIuDkSfv9R48Ct96qXuto2iZWDd3PLNJURLzKBYqjgBAXcQr5fwHwWQBTAOwH8GO7Rsx8LzO3MXPbsGHDYuyOkHfMhVQnAW/yzDPuycXsqK09dZ3Hzyyyvd3dPi+OAkKcxOZCycwHzNdEtBrAb+I6lyAA+gupfX3qYTBnzkA/5OPHgVdfHfig6OsDnnzyVC3dy7XN6n8/ZAgwaJA6fiV+HQXs/PplBis4EZuQJ6ImZt5vvP0qgK1u7QUhLH4WUo8fV2abM8/sF5Z33632zZgxUMj39vY/GHQEaqX/famkttpadey+PvW+vt6fo4CdX/+6deoYshYl2BFJIW8iWgPgQgBnAjgAYLnxfgoABrAXwH+zCH1bpJC3AATXVOfNU140dlGidpRK/Tn9TIE7diywfbv9w6JUUmaZNWu8++9U5LuuDrjgAuC99/xHL7odV4qGFxu3Qt6RaPLMPM/m459FcWyhWITRVN2iRO2wPgxM7xsnAW+20VkgdTMb9fYCQ4cCv/udXh91jyupPAQnqjbiVcgfYZPO+VlIdUI366QbOhkpgyTQk1QeQhBEyAuZIYr0upVeLy0tygxjdXV0yxZptrFDd4HUy//+8OFgKbG9jtvTI9lXhYGIkBcyQ1SaqjV2YutW4NlngYsuUr7qZ5yh0hI4CfpSSRVrCBNJreN/HyQlttdx33pL6ikIAymEkJfqUdnG/H3cwvtLJeWGGOR3fPJJ4KWXgPffV1r0wYPOD5P6emD16nCR1HFFulqPa/eQYpZ6CoINTnUB09jiqPEqdTKzTeXv47TV1DDX1/v/Hd1qrwLxjouNG5kbG/XrdPo5bnOzXu1ZoRggzzVe3ZDqUdnG7veppFRSbodESvP1+zu62fmtaYHjyHvkN9JVd8bZ3q78+52QRdjqJDaLg5P0T2OLWpPv7HTWEEXbSR+33wdQWnBnJ/OsWcF/x2nTotek/eA2k2hoUPuZ/c84OztVG6drMm9evP+XEC1hLQ4oqiYvLmfZxitCdfhwtdjY3R38d4yjuIwfjUsnmVmQGafUU8gPsVscnKR/GlscmrxoO9nF7fcx7eW1te6auNfvqKtJ6xJU49q4Uf2/ra3qr/W8QWecZl/MayjrTdVJFBYHuGjyqQt26xa1kF+9mrlcju4GF6LFa1FUZ9P5HaMShlE/MEzCmJTcHh5CdeD1+zc2ev+uhRTy5o1tpx2KtpMdli1zfhDrCHjzdzSF3bRp9sIuCmEY1xqPzDiLjdeMtnKs21E4Ie+mcZXLSsMXssP48cGE/IIF6vt2bphEzC0t0Wq2uhq31wPHxGw3frzMOIuM7ozWbSwUTsiLV011oaPJOA361avdb5D6+uhmbToat67N3ik+wOl7ug8OoTpxsjzoyq7CCfm03eYEfwS1zZdKKijIK5DKqgGFEZZeNnm3B05lH9xmmp///Kl9k4C+YhAmeK5wQl5snNWH3eJoTY23oNex55saUBTC0m0RV3cG6Wem6XexVzT+6iao7CqckI/LC0KIF7vF0VmzvIW4zjZ+fHRjwmkRV3cG6Wem6eeBIBp/9RNUdhVOyDOLD3FeiMLN0susE9U6ja4W5kdb87PYK4pNPggiu2IX8gDuB3AQwFbLZ0MArAew2/jb6HWcqP3kxYc4H1QOer9bQ4O3B08UQ09X0Oq0M8eum4228sEhzgb5wa/sSkLIXwBgWoWQ/xGAW4zXtwBY6XWcOLJQCvnAOugHD9YT7pX28iTWaXS1MLd2upk5rQ+OJB5iQnZJxFwDYEyFkN8JoMl43QRgp9cxRMgLJk4LiF7BU0TKNFOpASVpztDVwqx+8s3N6u8Xv+g9Y6l8cHhdE3E2yD9pCfn3La/J+r7iewsBdAPobm5ujv1iCNnHaQFxwQJv+7ybwM7iOo2u1m5uZmZOnYeXudXVqQVs8bjJL6kLeeP9Ya9jiCYvePmQuwnDcjlcorCkCbKo3Nh4qrD2StdsJnkTj5t84ybka0ImsXTjABE1MfN+ImqCWpgVBFfcinycPOn+3c99zrvoR3u7Xp3WJHD7X504fFhtW7YA69apmrVKh3Kmt7f/tTWF7Zw52bkWQnzEmU9+HYD5xuv5AB6P8VxaSK3X7OOVY96JUgmYNi36/sRJ0P8V6BfWBw86FyUHnI/vt76sUL1EIuSJaA2AjQDGE9E+IvpHAD8EMIuIdgPoMN6nxvLlqpL9ww8DmzcDjzwile3TwOtB61bkgwgol+33VWOhDLf/VZe+Pvfr5fa9ymIrogTlFCc7ThpblDZ5q3dGR4dKVJWEZ4XgjE5EppcXjLn4mqXF06BEEegFMI8caX9NWlr03UYlWra6QdEiXv14LBBJoEgS+HFh9PKCsS6ednRUt+eIW6CXH2Ftt6AcZXCWkG0KJeSDaEcjR4Y+reCB34hMHS+YvGifdg+tIMLaDh23UbffRpSgZAmaYK5QQt7LpcxuK5dFW4mbqNM/F037DOPj7/XA9PptRAlKhjBKS6GEvNeAFW0lfuy0kajTChQxV0tcPv5eipEoQfETVmlxE/Jx+smnwrhxyoe4r0//O8wDPQ2EYCxfrnywjx1T19X05+7sVB4wR48O/E4Qzxg390M7z5E8EJeP/+LFwKOPOsch9PUpd0vxqY8Pt5gJ09016PWP008+FRYvVkLDD6WSejgI4ejqUgL+6NF+AWz6c69dqwR9Q0O/y1+ppN4vWeJ/ALu5H8rv6Y/2dmD4cOf9ogTFT5xKS+6EfHu7EhqVwqSuDqittf8OkfKdF9/gcHhpI0ePAhs2AHPnKl/suXPVe68oVTvcHubV6DMfhCj92qdPd/arl4dm/MSqtDjZcdLY4vCTt9ovvfKSV6t3RlZIurZuFhOOJUXUnkVFW8jOGnHa5FMX7NYtiQRl1vSuTulZZVAHI43aullKOJYUUQhku8XxIj80s0CY6y9C3oYiemfEjWiDyRB27LrNAor40MwSQa+/m5DPnXeNLkX0zogbcz3kzjuVDd7Mq1JfH2xxVbAnzNi1Lo5bv2PNTLlmTbT9FfSJw4Mqdwuvuoh3RjysWKG3uCrJsIITZuzquOoJOcNJxU9jS9JcI6aF9MhLOoK0CDN2wyyOBw25F+IHLuaawmryTq6WQf22BT3cfOnvvFM0eh3CjN2gswBJ1V3FOEn/NLY0yv/JQlOyyIJ3dDiNXTeNO8gsQGa92QfiXSOkiVXoNDYGNxcI3til2SZS6YpNYezXVU8ezNnHTcgX1rtGSIbKXDZuyIJ3OOw8ZwB13bdtA2bMAG6+WS2Cz5mjFll37wbOOUdFEDuZecQTrbqJXcgT0V4AHwI4CeBjZm6L+5xCNnASOk4UJR1BXHgVBj9+/NQC3rrrTm5J/+TBnH2SWnidwcxTRMAXCy+hYyIL3tGgUxg8iJuk5AmKh6TciAvrXSPEj5fQqa8Pn6hM6EenMHgQ84p4okVPkt5KSQh5BvA7ItpERAsrdxLRQiLqJqLunp6eBLojJMW4cc6ZDQGgtxe4+24VYSmCIjw6abaDmld0g9wEb5J2Iyb2mt+FPQHRWcz8NhENB7AewH9n5uft2ra1tXF3d3es/RGSo6sLOP9852IURMDVV0sYfZQsXw6sXAmcOGG/v64OuOAC4NAhJezdFlyFeJg3T2nwTqK3owNYv97fMYlok5M5PHZNnpnfNv4eBPAYgPPiPqcQHWHshlKMInlWrACeew5oaTnVdFMqATU1SmN8+mkJaEoTLzPmM89E+5vEKuSJ6HQiGmy+BjAbwNY4zylERxR2Q69iFEOGSA6bqGlvB7ZuBV58UVXjam0FLrpICfne3oEmgpUrgdmzT/0NJLdQfHitnfT1RWy2cXKgj2IDcDaAl41tG4Clbu3TCoaSnBwDiSrK0e04NTXM9fWSwyYJvIp1W4ObamqYa2vld4kLt3siaJAZJOLVGUmWZU+UUY52EZZ1dUqQSKh8MnglJtPZ5HeJjmXLvB+6fsShm5AvtAulJMtyJsooRzvPjAsuAD7+2L69pLyNHh33Si/kd4mOFSvUAqsTUQaZFVrIS25tZ6LOt9/errxourvV30OHJFQ+SXTcK72Q3yVafvADFWtgR5RBZoUW8pKTw5m4oxylaEvyjB0bXpsfMqT/tSzOhiOxIDMnO04aW9I2+TQKT1cTQQoL6y5iS/ra5LDLTBl0q69Xx5O1LH287oko0p1DFl7tEUHjjZ8B6PfGD1OdXtBDx5OjciuX1ea0v65OCXu7fUTMs2bJvWOS1MNQhLwNpvAaOVINaNFIwhH0gSlFW+JF13XSunV0KEEd1hOn6PfQ6tXOD8uoH4Yi5CuofLoSqR+jqUkEjRtu004pLJFNgrhOtraKy2VYli1znw1F/TB0E/KFW3i1c5tkVvlVjhwpXi4P3cUzr+hXWcTOJn5dJ81Fb3G5DI4pY5xyNllJwl27cEJe3Cb70U1boBNPMHSo83nEWyY9/LpOmp5Tbt+rq9M7ZlEf7rp1FEziljuFE/KicSrcBPdtt6lcJqZ24fVgvP564HnbvKIKKSyRHk5uejU1QG2ts+uem3vfd76jygg6+XibFPXhrlO8xUrccqdwNV6llJnCTXAzq1SnL76obnSvB+OOHfbXE1ACXgpLpItTTVfAvc6r1/dGjwbefx84eNB+fBT14e4mY+yIXe44GevT2JJYeBW3SYXuwlpDg/ICcIon0PHUEPKDnUugmdBMXGEVft1Wo5A7kIXXfqSUmcLNhm7l+HE1FAcNCnaeTZskIjIvOJn4Pv5Y3UMzZ0rVKMBexjhBpCKRY8VJ+qexpeEnb/XPLkrK4WXLVECLrqYxcqTS1sK40xVdu8sD4ibrj40b1UzWaxYcxb0B8ZP3pihh2n6nkmYMQVif6SKaxPKGl4kvxds3s/gJRgtzb7gJ+cKZa+woUsphv+5dzHr+vroUzU01T0hSOf/48bSJ694QIY9i+c77de8KQrnsvK9Ibqp5I+7MpHnET1BZXPdG7EKeiC4hop1EtIeIbon7fEEoku+826AjAkaMCBfpWCoBZ50lGl8esVtQJFIP9TPOUMpQnma9UeAnGC2ueyPuQt5lAD8FMAfARADziGhinOcMQpGmoW6D7rTTgHPPDafp19cDt94qGl9esVb5GjlS3R99fcD+/cEKvecdP542cd0bcWvy5wHYw8yvM/NHANYCuCLmc/qmSNNQLxdSt4pNJl4Rk9dfL26qeaa9Xd0zH3yg1msq17FWrlQR0+I6q6gsf9nSouRKYveG04psFBuAKwHcZ3l/HYC7K9osBNANoLu5uTnY0nIEFC23uVOKX7dCKgDzkCEDXU6d0gRLGuH8ous1kvf7KChR3xtw8a4hjnEVjoiuBHAJM19vvL8OwH9m5hvt2re1tXF3d3ds/fGiq8s9zLsIdHWpKffRowP3NTQojaRo10QYSGurSmqni4ydeCGiTczcZrcvbnPN2wBGW96PMj7LJJXFpos4ICUiWNDBbyrivHmpVZLlerdxC/k/AjiHiMYS0SAAnQDWxXzOSMnyjxcXlTZEtzD1Il4fwX8K47x5qVnRTdmdGk52nKg2AJcC2AXgNQBL3dqmGfFqR1GiYIMi16fYVK5jeUVOz5uXdo+jJysJDyFpDfyTlR8vq8j1EZj7FxAbG92FfLmczzGRlXw+bkJeIl4dyFMUbBwmlTxdHyE45jqWVybFYcPyuZ5TDYGUIuQdqIYfT4e47IV5uT5CNHhFUs+YkWx/kqIaAilFyDtQDT+eF3EmXsvD9RGiwyuSOk8BhVaqIZBShLwD1fDjeRGnSSUP10eIDi/XWyCfXljV4HIsQt6BavjxvIjTpJKH6yNEi5PrLZBxF8OQ+HE5ToNYI179knbEqx3VHAU7b566oZyKll99NfDQQ/3/465dyszi9T9a2w8Zomyu771XfddHiB+JoE4Gt4jX1N0mrVuWXCjzgI6bo19fd/GNLx5hymJmxcUw70BcKIuJjp3Uz8JskSpoCYqw3lnihZU+IuQdqJZwfa9+utkL/S7Mim98sYjioS5eWBnAScVPY8uKuaZaTBJh++m3MLMUci4WUZhaJDI6GSDmGn2qxSQRt5YFAK+/fursQLSyYhGFqUW8sNJHhHwF1WKSiKKfXpkEDx8+1QYrvvHFIqqHuq6LYbWYSKsOJxU/jS0L5ppqMUlE1U/dTIKV3jhFqaBVZJI0tVSLiTSrQMw1+uhoL1nQOOLQshobnduZs4OsB34I0ZGUqaVaTKRVi5P0T2PLgia/cSNzfb299lIuMw8dqv7qaBxh/It1+hm1llUtsxghWcLWI/W6D8SXPjyQfPL6LFvGXFPjLuzcBKs5oEeO1H8YhOlrlKYTtyLepVI+iz4I8aJjhhHlIjwi5DVx0469tlKJuaXF+/tBtWwnbSjKqu/i7iZEie54clMuAObmZhl7XqQi5AF8H6po9xZju9TrO2kLebdpo66g12njd/qZ5KKULKwKUaFrhvFSrohkDHrhJuTjXnhdxcxTjO2JmM8VGje/YB3sEoHZtfETyh31olSYCFlB8IOun711gZdoYFtmWYQNg3jXWPAKDooCv0FDbv7wR4+qTJK6A183D4lZ0q27W/2VgBUhCH48wEzlYvRo5+NlKU6lmohbyN9IRK8Q0f1EZOugR0QLiaibiLp7enpi7o47XsFBUaAbNGRq3I8/7j67ePNNvYRR4qYmJI3f4Ln2duDMM52Pl8WEZllwp/bEyY6jswF4GsBWm+0KACMAlKEeJLcDuN/reGnb5JntbdLlcnA7fRC7YqUN3o93jxPipiakgd81nix6eDk5PWQpgAtpe9cAGANgq1e7LAh55oEeK7NmBVuQJVIPiKYmZ8+XygG0enUwDx8vQS1uakJa2HmAuXmLZcnDy0mQL1iQrX6mIuQBNFle/w8Aa72+kxUhX4mXgLTb/uZvvF0a7QZQmFmD2+XLooYkFBMvDTgrHl5uDxw3T7o0ZsZuQj5Om/yPiOjPRPQKgBmGoK9KgizIfvSR+p7ToqWTjfzkyWB99FrQleRiQhZwWxu67TZg9mxgzpxseHi5OT24edJlbu3ASfqnsWVVkw8aJOU2bQvik2+af4JOEbOiIQnFwM4kozPuszImg8zg05oZI22bvO6WVSHPrJ+tUXfaFmQAWW2BQQV1lBGyguCEk0mmqSm8ghQ1TusDXpG4Tlu5zDx+fLL3lwj5iDAHQ2Oj/g/u9C95DaBy2VmQewnqOBOjCYIXbjNfaz6noApSlFQ+jKzOEh0dzskKdZW8pGYlIuQjRtd84zZt8/IiWL06mMadJbcuoZi4mWTcTI66ClJU6KRTqKlhrq0NptGbW11d/IqWCPkY0PFl95pyRm0jz5r7mVBMvEyRI0d6K0lx2bWts9zmZr1ZRV2dcqNubVXa/aBB/gV9S0v0/4sVNyEvaQ0CYoZhX301MHIkUC73593QLawQdZ6YaildKOQbr3QGM2aocd7R4dwuDo+vyrQeb76pRLAXvb3A0KEqzceXvqSXo6qSHTtSjIZ1kv5pbNWkyVeShQVNCXgSsoCfGaXXbDaq9aUwacTNeyeo37y5xbnGADHXJEvQgRl2QLst5hKp6aksxgpJ4McU6aQgmcewjuNyWXmY+SVMGnHTdOS11qDzoIgLEfIJEnThM4oFU52FpChs/4KgQ5jZrddYdhP0dspSUJ9388GyerX3Mdw8ceL2nRchnxBBFz6jXDCt1KDCLAwLQlp0dnoLXrux66QstbQE95Axkwu6HaNUUouzToI+7nvNTcjLwmuEBF349MoZf+ut3uc2U57+5jdqcWjmTLWYO3q0fSEGrz4JQprs2uW+/+TJgWPXLWXCnj3AoEHB+sLsfYz6euAHPwBuvlk5XZgLyrpOGHEiQj5CdCvh+PkeADzzjHu++EqvgWeeAV58EbjsMpWfO0ifBCFNdArrVI5dN2Wptxf47GeVwA2K9RhOQjyLldVEyEeIn0o4ut8DlDB2KuzhVQxk6NBgfRKEONAtsrF4sXJLdoJo4Nj1UrLq671dN51mvZXHcBPimaus5mTHSWMTm7zz5hTm7VUMpKNDAqSEbODXuWDBAn/rSX7Sadt5/9TVqRThaS2ehgGy8JocQaNYly3zdsOyuzw6vvGSfVJIGzdFhkgtWtopHAsWnJoGwa3Cml8ly+r909Linacmy0qRCPmECeo6NmuWfy1C1ze+o6M/NFv85IWkCZNi2M/9FESh0XE9zrpS5CbkSe3PBm1tbdzd3Z12N1Kjq0stoB49OnBfQ4NaWNqwQdkex43rX2zdsMHZFkmk9pVKyp64ZEm6i0BCMWltVU4BXjQ0qPEcxo7d1aXuld27gXPOUfb9yuOZbXbtAt59F3jrLed7qLlZ3Wep29ZdIKJNzNxmu0+EfLZYvlwtmB4/rhZ6TOHc1qYWco4dc/fE8SKKm0gQ/DJvHvDII955X0oltZi5Zk18fTHvMd17qbVV3XtZxk3Ih/KuIaKriGgbEfURUVvFvu8S0R4i2klEF4c5T5Gwc8G66y41yKweNEER33ghDdzKT1rRcevV9dCxa2fnjeZGLjzQnOw4OhuACQDGA3gOQJvl84kAXgZQB2AsgNcAlL2OlxebfNSEybuhu4ArCHFjl4vGrweLroeOW+Srn3spy4utVhBXxCsz72DmnTa7rgCwlplPMPMbAPYAOC/MuaoZXc3DCa9gKT/kQjMRqhJzlho0xbBXTIh5X7m1275dX4NPO1I1KuIKhjoLwFuW9/uMzwZARAuJqJuIunt6emLqTnpURqM+8oh67xbBWolXsJQf4sjTLQi6tLcD69cD3/ue//B/3bQhbu3cBDyRWmTNSqRqVHiKDiJ6moi22mxXRNEBZr6XmduYuW3YsGFRHDIz6GoeXujaM90gyo9mIlQ/QcL/ddOGeM18nRSm005TyphTpGrYGXla1Hg1YOaOAMd9G8Boy/tRxmeFwivx2NVX67lmtbcr4Wz1ujGxukiWSur9yZOqDZH6bPhwYPp0e1cyQUiL9nZ/43HcOGDLFnsPHdMM2dWlXCKdKJWAiROB118f6MHmpgBVeuRs2QKsW1cdLsmRuFAS0XMAljBzt/G+BcBDUHb4zwDYAOAcZj7pdpy8uVDq+Aab2rXOQKn0/505U2k/Vn9gwNtHWBCqEa84ks5OYO1a+/3Wdhs2qNe694nXebPgkuzmQhnWu+arUPb2EwAOAHjKsm8plFfNTgBzdI6XN+8at2jUalzBF4S0cYpoXbAgvqhVr/xQcZb10wUxetc8xsyjmLmOmUcw88WWfbcz82eZeTwzPxnmPNWKri1dfNcFQeFl93ay5R896mwaBVRdhaALqUFTiGcFSTUcI6Yt3SuHtdtAqdbFHkHwi64nmpnK9+671fsbbgCeesp9sXXYsOAmlaApxDODk4qfxpY3c43Jxo0qUZjfAJAo6r4KQjXgN4Nk5b3htoUtYh9lec64gGShTJ8gaVCzPrAEISp07N5mNsrx409NP6yzhVWUsp6u203Ii7kmIaymG50AkKD1YgWhGvGyez/3XL8pZ+dO5SbshbXKk3nsIHEqQDbL+uni6ScvRMeKFcCcOXquW9W+2CMIfnDzgScCenr0BDsANDYCZ5+tvuOUQthUlPzY6f369WcFEfIJoztQdAI/BCEvLF6sgovsfNFLJe8Uxda2l1wCPPSQ0rjffNO+XZEUJTHXJIjVU2bWLGD2bGevGTf3S8k/I+QNN3Pm8OH6Cfqs90bVe8VEhZOxPo0tzwuvbt4AXulSs7rYIwhRY1fqTyeo0O7eKJLzAqT8X7q4hUVbsQuR1illJgh5xu3+KZfVfTFliv294VRprRpyzvhByv+lzLx5yivA61InUfpMEKqRMMI6jKJkrQU7blx2lSwR8imjW8QYUJ4BTzyRzYEkCGmS9Ky2MvNklmcBIuRTRreIsYmfzJSCIERPNWSetBJbIW9BD79FP4IEawiC4IzfHFB5CkYUIZ8Adu5hXlTbQBKErBKkBGeeghFFyCdEZVh0RwcweLBz+2obSIKQRYKW4MyTj70I+QQxU6R2d6tixpddlp+BJAhp4GWGCWp2yVMwogj5FMnTQBKEpKk0w6xdC3zpS8AXvtAv7IOaXfwmFMwyoYQ8EV1FRNuIqI+I2iyfjyGiY0S0xdjuCd/V/JGngSQISWJnhgHU623bgBkz1EMgjNmlmjNPWgmboGwrgL8H8L9t9r3GzFNCHj/3+MlMKQiCws0MAyhTzJ13qnZOic90ZsvVmnnSSighz8w7AICsiZsF3+RhIAlCkriZYUyOH1ea95IlztGyRbjv4kw1PJaI/gTgAwDfY+bf2zUiooUAFgJAc3NzjN0RBCEvuKXiNjFt7mvWFHu27CnkiehpACNtdi1l5scdvrYfQDMzHyKiVgC/JqIWZv6gsiEz3wvgXkBFvOp3XRCEouKWf97EanMv8mzZU8gzc4ffgzLzCQAnjNebiOg1AOMA5C9ngSAIiWM6LaxcCZw4Yd9GPNQUsbhQEtEwIiobr88GcA6A1+M4lyAIxWTFClX7taXlVA8a8VA7lVA2eSL6KoD/BWAYgN8S0RZmvhjABQB+QES9APoALGLm90L3VhAEwUJ7O7B1q9RdcEOyUAqCIFQ5koVSEAShoIiQFwRByDFx+slHQm9vL/bt24fjx4+n3ZXCUF9fj1GjRqG2tjbtrgiCEJLMC/l9+/Zh8ODBGDNmjETWJgAz49ChQ9i3bx/Gjh2bdncEQQhJ5s01x48fx9ChQ0XAJwQRYejQoTJzEoSckHkhD0hunKSR6y0I+aEqhLwf/NZyFARByDO5EvJBajnqUC6XMWXKFLS0tGDy5Mn48Y9/jD63zEgA9u7di4ceeijciTW4/vrrsX37dtc2v/71rz3bCIKQT3Ij5IPWctThtNNOw5YtW7Bt2zasX78eTz75JFZ4VA5ISsjfd999mDhxomsbEfKCUFxyI+SD1nL0y/Dhw3Hvvffi7rvvBjNj7969+MpXvoJp06Zh2rRpeOmllwAAt9xyC37/+99jypQpWLVqlWM7K3v37sXnP/95XHPNNZgwYQKuvPJKHDXS7G3YsAFTp07FpEmTsGDBApwwsjJdeOGFMKOEP/WpT2Hp0qWYPHky2tvbceDAAbz00ktYt24dbrrpJkyZMgWvvfYafvKTn2DixIk499xz0dnZGc2FEQQhmzBzZrbW1lauZPv27QM+s2PaNGalw9tvNofW5vTTTx/w2RlnnMHvvPMO//Wvf+Vjx44xM/OuXbvY/B+effZZvuyyyz5p79TOyhtvvMEA+IUXXmBm5q997Wt8xx138LFjx3jUqFG8c+dOZma+7rrreNWqVczMPH36dP7jH//IzMwAeN26dczMfNNNN/Ftt93GzMzz58/nRx999JPzNDU18fHjx5mZ+fDhw7b/s+51FwQhfQB0s4NczY0mH6aWYxh6e3vx9a9/HZMmTcJVV13laBbRbTd69Gh8+ctfBgBce+21eOGFF7Bz506MHTsW44x/Yv78+Xj++ecHfHfQoEG4/PLLAQCtra3Yu3ev7TnOPfdcXHPNNfjFL36BmprMh0oIghCC3Aj5xYtV/mg7os4r/frrr6NcLmP48OFYtWoVRowYgZdffhnd3d346KOPbL+j267SfdGPO2Ntbe0n7cvlMj7++GPbdr/97W9xww03YPPmzfjiF7/o2E4QhOonN0LeLCLQ0NCv0ceRV7qnpweLFi3CjTfeCCLCkSNH0NTUhFKphAcffBAnT54EAAwePBgffvjhJ99zalfJm2++iY0bNwIAHnroIZx//vkYP3489u7diz179gAAHnzwQUyfPl27z9a+9PX14a233sKMGTOwcuVKHDlyBH/5y18CXQtBELJPboQ8oIoIbNgAzJ2r/OTnzlXvPRxhPDl27NgnLpQdHR2YPXs2lht+md/85jfxwAMPYPLkyXj11Vdx+umnA1AmkXK5jMmTJ2PVqlWO7SoZP348fvrTn2LChAk4fPgwvvGNb6C+vh4///nPcdVVV2HSpEkolUpYtGiRdv87Oztxxx13YOrUqdi9ezeuvfZaTJo0CVOnTsW3vvUtfPrTnw53gQRByCyZzye/Y8cOTJgwIaUeJcvevXtx+eWXY+vWrWl3pVDXXRCqndjyyRPRHUT0KhG9QkSPEdGnLfu+S0R7iGgnEV0c5jyCIAhCMMKaa9YD+AIznwtgF4DvAgARTQTQCaAFwCUA/tms+So4M2bMmExo8YIg5IdQQp6Zf8fMpmtGF4BRxusrAKxl5hPM/AaAPQDOC3GeMN0UfCLXWxDyQ5QLrwsAPGm8PgvAW5Z9+4zPfFNfX49Dhw6J4EkINvLJ1zv5owqCUFV4RsIQ0dMARtrsWsrMjxttlgL4GMAv/XaAiBYCWAgAzc3NA/aPGjUK+/btQ09Pj99DCwExK0MJglD9eAp5Zu5w209E/wDgcgAzuV/dfhvAaEuzUcZndse/F8C9gPKuqdxfW1srFYoEQRACEta75hIANwP4W2Y+atm1DkAnEdUR0VgA5wD4f2HOJQiCIPgnbOKSuwHUAVhvhNN3MfMiZt5GRI8A2A5lxrmBme1DPAVBEITYCCXkmflzLvtuB3B7mOMLgiAI4chUxCsR9QD49xCHOBPAuxF1J0qkX/6QfvlD+uWPPPbrPzHzMLsdmRLyYSGibqfQ3jSRfvlD+uUP6Zc/itavXCUoEwRBEE5FhLwgCEKOyZuQvzftDjgg/fKH9Msf0i9/FKpfubLJC4IgCKeSN01eEARBsCBCXhAEIcdUlZAnoquIaBsR9RFRW8U+zyIlRDSWiP5gtHuYiAbF1M+HiWiLse0loi0O7fYS0Z+Ndt12bSLu1/eJ6G1L3y51aHeJcR33ENEtCfTLsfhMRbvYr5fX/26k6njY2P8HIhoTRz9szjuaiJ4lou3GPbDYps2FRHTE8vsuS6hvrr8LKX5iXLNXiGhaAn0ab7kOW4joAyL6dkWbRK4XEd1PRAeJaKvlsyFEtJ6Idht/Gx2+O99os5uI5gfqADNXzQZgAoDxAJ4D0Gb5fCKAl6FSLIwF8BqAss33HwHQaby+B8A3EujzjwEsc9i3F8CZCV6/7wNY4tGmbFy/swEMMq7rxJj7NRtAjfF6JYCVaVwvnf8dwDcB3GO87gTwcEK/XROAacbrwVBFeir7diGA3yQ1nnR/FwCXQqUhJwDtAP6QcP/KAN6BChhK/HoBuADANABbLZ/9CMAtxutb7MY8gCEAXjf+NhqvG/2ev6o0eWbewcw7bXZ5FikhlVznIgD/1/joAQB/F2N3zXPOBbAmzvNEzHkA9jDz68z8EYC1UNc3Nti5+EzS6PzvV0CNHUCNpZnG7xwrzLyfmTcbrz8EsAMBazSkwBUA/g8rugB8moiaEjz/TACvMXOYaPrAMPPzAN6r+Ng6jpxk0cUA1jPze8x8GKoS3yV+z19VQt4FnSIlQwG8bxEmgQuZ+OArAA4w826H/Qzgd0S0ycirnwQ3GlPm+x2miJEVfAmItfhMJXFfL53//ZM2xlg6AjW2EsMwEU0F8Aeb3f+FiF4moieJqCWhLnn9LmmPqU44K1ppXC8AGMHM+43X7wAYYdMmkusWNgtl5JBGkZIsoNnPeXDX4s9n5reJaDhUJs9Xjad+LP0C8C8AboO6KW+DMiUtCHO+KPrF+sVnIr9e1QYRfQrAvwH4NjN/ULF7M5RJ4i/GesuvodJ8x01mfxdj3e1vYdSfriCt63UKzMxEFJsve+aEPHsUKXFAp0jJIahpYo2hgTkWMtHBq59EVAPg7wG0uhzjbePvQSJ6DMpcEOrm0L1+RLQawG9sdmkXfImyX2RffKbyGJFfrwp0/nezzT7jNz4DamzFDhHVQgn4XzLzryr3W4U+Mz9BRP9MRGcyc6zJuDR+l1jGlCZzAGxm5gOVO9K6XgYHiKiJmfcbpquDNm3ehlo3MBkFtR7pi7yYazyLlBiC41kAVxofzQcQ58ygA8CrzLzPbicRnU5Eg83XUIuPW+3aRkWFHfSrDuf7I4BzSHkiDYKa6q6LuV9OxWesbZK4Xjr/+zqosQOosfSM00MpSgy7/88A7GDmf3JoM9JcHyCi86Du71gfQJq/yzoA/9XwsmkHcMRiqogbx9l0GtfLgnUcOcmipwDMJqJGw7Q62/jMH3GvLEe5QQmmfQBOADgA4CnLvqVQnhE7AcyxfP4EgM8Yr8+GEv57ADwKoC7Gvv4rgEUVn30GwBOWvrxsbNugzBZxX78HAfwZwCvGIGuq7Jfx/lIo743XEurXHijb4xZju6eyX0ldL7v/HcAPoB5AAFBvjJ09xlg6O+7rY5z3fCgz2yuW63QpgEXmOANwo3FtXoZawP5SAv2y/V0q+kUAfmpc0z/D4hkXc99OhxLaZ1g+S/x6QT1k9gPoNeTXP0Kt42wAsBvA0wCGGG3bANxn+e4CY6ztAfC1IOeXtAaCIAg5Ji/mGkEQBMEGEfKCIAg5RoS8IAhCjhEhLwiCkGNEyAuCIOQYEfKCIAg5RoS8IAhCjvn/v4U39LR7Uv4AAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# x_data and y_data have already been defined\n", "\n", "fig=plt.figure()\n", "ax=fig.add_subplot(111,xlim=[-11,11],ylim=[-21,21])\n", "x,y=poly_curve(smarap,x_data)\n", "ax.plot(x_data,y_data,'.b',markersize=15)\n", "ax.legend(['Data points'])\n", "print(len(x_data))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Generating the matrix D" ] }, { "cell_type": "code", "execution_count": 63, "metadata": {}, "outputs": [], "source": [ "degree=4\n", "D = data_matrix(x_data,degree)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Finding the least squares solution" ] }, { "cell_type": "code", "execution_count": 64, "metadata": {}, "outputs": [], "source": [ "\n", "params = np.linalg.inv(np.transpose(D) @ D) @ np.transpose(D) @ y_data" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Plot Curve Fit\n", "Use the next block to plot the fitted polynomial and find the magnitute of error. \n", "\n", "Once you're done, re-run all the blocks of code for this problem for $\\textbf{different degrees}$! What do you observe?" ] }, { "cell_type": "code", "execution_count": 65, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "38.596128739466714\n" ] }, { "data": { "text/plain": [ "Text(0.5, 1.0, 'Polynomial of Degree 4, ||error||=38.596129')" ] }, "execution_count": 65, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYQAAAEKCAYAAAASByJ7AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAABY5klEQVR4nO2dd5wURfbAv29nYRdEJUhSRDICBhDUNZwRFNE7zwzeeSqG09MT/YlZCaeeiTN7egb0TBhPjzMLZlxQRCRnUEGQjBJ2gd33+6N62N7Z7pme3Z6wO/X9fPqzO93V1dXd1fXqvXr1SlQVi8VisVjyMl0Ai8VisWQHViBYLBaLBbACwWKxWCwOViBYLBaLBbACwWKxWCwOViBYLBaLBUiTQBCR80REXduvIvKdiFwuIvlJ5vWJiHySoqKmBRFp5zyH86pxbqj3LyK/FZHpIlLilKmxT7oRHu9wvoi8KCLHh1We2oCIdBCRzc5z6FSDfKrUARF5prbV71Teh4jsJSL/FZHvRWSLiKwWkU9FZIBH2rYi8m8R+cFJO09EbhORnQJc55mY+h3d7vdIe6mIzBGRUudat4pIPZ98/yQiXzv1Zb2IfCEi+7qOtxGRh0Sk2FWn2nnk00dEHneuu9m57gsi0t4j7W4iMlpEVjnPYVLQbzSpxjgEzgCWArs4/z8EtACGpbkcmWY5cAiwMJOFcITxC8CXwGXAVuDXBKcdDpQBDYH2wOnAeyLyPHCuqpanrsRZwz+BDUCDTBckB2gErAZupqLtuAh4W0ROU9X/ADiN/jigHnAL8ANwIDAS6AycFeBaq4Dfxexb7v4hIjcAtwP3Ae8BPZ1rtAYujEn7d+BK4G7gWsw3c5DzN0on4EzgG+Bz4Difsg0EegAPAjOBPZz7nCwiPVX1R+eaBcBHwG7ONVcAFwBviUg/Vf0k7hNQ1ZRvwHmAAp1i9n8MbEgyr0+AT9JR7mzcwrx/YC/nvQwOkHaEkzbf49hVzrGr0/ws6gGS5mueDfzsfOhV6nSSeSlwXsy+Z0J8vwLU9zlWUMO8C1z/p/Q+PK6dD/wI/M+17zinHMfFpL0T2A40TJDnM8DSBGkKMR2mZ2L2DwXKgR6ufYc4+36fIM881/8XOvfQziNdc499eznX+Jtr3x+dPI6KqQfTgK8SPdtMjyF8DewiIi0ARKS/ozptEZENIvKmiHT1O1lEWonIVhEZ4nFshKNaNXF+f+Koa31FZIpzbIaInOJxbsJyuPLrLyJTnbTfisjBIpIvIn8XkeUistZRR3dynVvFZCQiB4rIayKy1MlrrpNHtXqhItJaRJ51VOxSEZkmIn90Px9gifPzKac8n1TnWqp6H/AtUOk9iEhzEXlMRJY5ZZgjIhd7lLWv8+xKRGSBiFzoPLMlrjTRZ/YXEblbRH4CSoHGzvFTRWSiSzV/VUTaelzrYjHmyhLn2TwlIk2D3KdTl+7FNADrgz6fMAhyfyKyRESeF5HBIjIHo/GdKBUm2yOc89YDk5xzdhGRh0XkJ+cdzRWRq0REXPke5Zx/qog8ISKrMEIxI6jqdoyGtt21u77z95eY5OsxpnGh5uyD0Vjejdn/npP/7137LgUWq+qb8TLUgBq1qq7y2Pc9RqvZw7W7CNiiLk1AjVT4ADhQRPYgDpkWCO0x5oeNItIfeBvYiFHvLsW8gC/8bkJVVwBvApUaGRGJYNSkV1R1netQR+ABzEd9KkYdfFVcduAky9EJuAfTCzkDKADGAo9iVMjzgL8BfwCGJ3gWbYGpwCVAf6ecg4GnE5xXBUf4fAqcANyIqajTgedcDfKTTpkBbsP0aP6S7LVcvAvsGW2kRGQX4AtgAEa7OBH4H/CoiPzVVdbuVDzvgU55hwDH+FznJqAL5p2fApSIyCXA68AsjAnrz5h39qmI7Oy61p3AIxjTwu+AazDP+l2nziTibmCOqj4XIG1oBL0/h6OB/8OYMfpjeoZRXgAWO3lcLyJ5mGd/PvAP4LeYxu1ejFkklocwDd85mLqd7H3kOZ2lRFuVdsl1bisRGYapAw+7kowD5gN3iUh3EWkkIsdg6tJjqropQBFbOJ2E7WLGH66LqRdlzt+tMeeVOn/3ce07HPhORK51OkTbxXRAzyAkRKQbxuQ+O6aM2zySe5WxKqlQ6TxUm/MwakxXjLrXBFOpy4A3nTSTMS8033Vee+fm7nXt+wSXKgoc5eT9G9e+3zn7imLO2wZ0du1r4ZThRte+ZMqxDejgcd1xMff/H0xvIfq7HR5qdoyKl49R/8qBZn7373P+5cSojc7+ccBKIOL87hSvHDHnjsDHZOQc/7Nz/GDn9y1Aift5O/ufwNiE853fL2J6OQ1daVo75y7xeGZTcJmJMD22DcDomOu0x3y4V7rOLwOGxaQ7zMk3kWr/G8xH1T2mTqfUZBT0/px9S4DNQCuf7+++mP0n+ZThSeded4v5xt6o7n3E1KFE2wiPa4xyHf8VONUjTQuMHd6d1xO4zDJx3sWVwF8xHZEBznnlwJMx76IMuCvm3D8513rfta8Eo60sxpgZ+wGvOulO9imDr8nII20+ptO3Emji2v8XJ49uMek/cvYPipdvujWEOZhGdC1mYO4FYLDToz0AeFmNOgiAqi4GJgBH+mWoRjWahWmQovwZmKaqE2OSz1fV+a5zV2IeaLRXm2w55qnqopj7A3jf477buNXwWBzV/S4RWYj5GLcBz2GEQ2e/83w4AlimVQeQngeaA92TzC8I0XtT529/jFlisbv3h3k2zVxlKALeUdXN0YxUdTlmoNuLN9Wp4Q6HYAYaX4i5zo+Y536Ek64fRiOOTTcJ08AcgQ8iUh/4F6ZBnZXwSYRL0PuLMlGN5uzFGzG/j8A0ei/G7H8eY4I5JMH5yfI4ZqA30fa4x7n3O8d+i9FGXxSRk6IHRaQQeBkjFM7BfKvXYLT8RxIVTFXvV9WHVPUjVX1HVS/CaOkXiEhnJ81GYDRwuYgMFJHGInI08HeMoHCbf/KAnYHTVPVFVf3QKctMjBZcUx4GDgX+qJWtIC9iOlz/FpF9xXgc3UhFPYlrokq3l9EpGE+BX4HvVbUEjOsVpkFZ7nHOCszgSTweBUaJGUtohGmMLvdIt9ZjXylmsAiM5pJMOdbF/N4aZ38+EKGy3dPN00BfjMfVVGATxiPhEVf5gtIU/3uIHg+bPZ2/0eu2wGggXuorGKEARhtY6XH8Z6CDx/7Y+2rh/B3nc511MekWJCiPF1di6saDUuGWG/UU2VlEdlbVRN5Z1SXo/UXxeu9+x5oCa1U11gTiV0/i5R2EFXi/61iqNFqquhTTdoDxmPkEozW85ey7AKPJdFLVqPfeZyKyAXhcRB5T1e+SLO8YzLvvg7EaAFyNqSsvYtqKEsw3ey2Vn88azID+FNc9lIvIeIxZuNo4ps+LMV59H7iPqep6ETkV+DcV5sKFGO3sVhK8w3QLhBmq6vVBrsP0LFt5HGuFd0Pu5lngDoxq3ASjNr9QjfLVtBzVwundnIxRlR9w7d/X/6y4rMWY52Jp5ToeNgOAH9Rxf8N8ECuJGWh2Mdf5u5yKRs9NS5/zNOb3GufveZjeVyy/xqQ7jqqNqPu4F90xz26Zx7EpwHcY98NUEPT+osQ+n3jH1gJNRaR+jFDwqyfx8g7CMBKPpYEZ/xiRIM1kTGMdZV9gnUsYRPnK+dsN856qw477VtVfgFNFpDnmOS3BdA7uxoyZRZkJ9EqUX7KIyE3AdcBf1WcsS1U/F5GOmA5ZBJiH0Za2YNxbfUm3QPBEVTeJyDfAGSIyQlXLwExKwahFDyU4/xcReQFjKmoEjHFeXFrLUQMKMC8utjd9XjXz+xRzD4ep6gTX/rMxjXSoZg8RuQrTIF7p2v0exib7g2Oa82MiMEBEGkbNRiLSGmPbD9Ij/RLTKHZS1X/HSfchpufZ1lHfk+FOjE3cTX/Mh/lHKoRbKgh6f9XhU0xDcQaVO1B/wGi1xSFf73EqevTx+CneQWfQ+XAqz+NZATQRkU4xnc6Dnb9ewjwRf8A03l/HHlDj9bPKKc9NGDPNq64kbwDHikgfVZ3sKnc/r/yCICJXYBxAblLVh+Oldcyq853zGmHmbjynCQbXs0IgONyC8Xh4S0T+iWnYR2IG1P4R4Px/UjGO8FgGy5E0qrpBRCYCV4vIckzlGkxld7JkeAbTM/+PU1mXYip3P+DPUUFXTQ4WkTKMGasDxmPlBIyK+qAr3X0Ym+nnInIfptHcCdgb4wBwspPuNieP90VkFEY43oIxGSV0yXM6A9cAjzi9tncx72oPjB35E8eGu1BE7gIeFuNC/ClG3d8T81yeVNWPfa4xh4rxIcC4wTr/TnI3QGJciZ8GjvYYw0maoPdXzezfxfRqH3PynonR9C4E7lDV1TUtvxtV/YkEjX0sYtyjm2LG8FZgeuUXYMypZ7uSPoPxrnpHRG7HTEzrg6lL3zjnR/McD+ylqp2c33thxutewpgUCzDm7fOAf7m1DhE5yynPXIw14lTMxLLTYsyGT2Eme74uIjdjvumLMZp7pclnInK6829v5+8JYlx7V6nqp06agZhxlPeAj0SkyJXFL+6xLRG5w7nn1Rgt4RpMZ/MGEpFoNDuMjYAeGZheVzFGtdkA/BfoGpPmE3y8bJyX9LXPsU+ALzz2L6HqRJOg5fgiZl875z4vjNk/ApeHDh5eRs6+dzG9wZWYQaMTqTrJxPf+Y67ZGlPJV2PGSaZhBqDcaarjZRTdNmE+nheB433OaYIRDIsxPc6VGC+QK2PS9cOMm5QCizCC/Q3g20TP1nV8AGai4y8Yk+F8zABg95h052C0kk0YV9fZzrNuE0adxjQCVbw8fPII5J0T9P6cuvx8Mt8fZsD6YYw2thVjXriKyp5cRznn963pfSS7YTz3PnLqTinwPca1+zCPtN2BVzAD7lucexmFywvH9Q0tcf1uinFf/x7TSdiMMQVeToyHEqbxn+6k+QXj31+lLK5v8HmM6a0E06Yc55FOfbZPXGmeCZLOSTsa0wnc6vx9CGga5HmLk0Gtx+n1zQYuUtWnMl0eS/VxVNwFwNuqekGmy5MMIvIi0FhVq8Ta8UirwPmq+oxr3zMYt8OjUlXGsKkr92HJLpNRtXA8lDphzDrLqepCZ8lyROQhjK38J2B3jLmrCcbtr7ZxBKYXabHUOmq9QMDYO4dh1MOzVXVLhstjSZ5C4C6MZ9FWjGdIX1WdFvesLERV22S6DBZLdan1AkFVR5DYRc2SxaiZBGSxWDJMnRlDsFgsFkvNyCoNYbfddtN27dolf+K0abDLLlCdcy0W4Kef4Oefodzl6JqXBy1bwu67Z65cFksl1q2DRYuge3doUBEI+Ztvvlmtqs1rmn1WCYR27doxefLk5E888kjzJX/+efiFstR5Jk6EY4+tLAzA/N6wAf7zHygq8j7XYkkrd98N110HxcWmE+wgIt+HkX2Ng9uJyJ4i8rGIzBKRmU48IUSkqYh8KGaZxQ/FWZcgJXToYKQm5uMeNAh69zZ/J8aGt7NYqFxPzjoLtvi4IpSUwAO10dfJUjdZsgSaNq0kDMIkDA1hO2alrCliYrN/IyIfYibDjFfVO0XkeuB6zFT/8OnQAX76iVtv3MKdDzRgyxZQhalTYexYGDoURo5MyZUttZDhw2HUKHbUk3iUl8P8+fHTWCxpY/HilJrGa6whqOpydSL6qZm6PRszrf5kTDgDnL+/r+m1fOlggmK+fNcSNm+u+MjLy2HzZvPxW03BAqYejBpFpXoSj7w86NKl4lyrfVoyypIl0L59yrIPdT0EJ75LL0yM+ZZq4tqDiUHiF72y2kQ/0BOv6AjAXuWLPNNZtd8S5YEH/M1DXhQWwhVXGK3i2GPh5ZdhyhR45RXze3iQ2J0WSxioGoGQzRpCFCfcwOuYODWVIo2q8W317I+JWeN2sohMXrWqyrKhvkQ/0Jdegq/XGg2hI7GRbw1W7bdEmTcvuGbQsKExN0JVrcJqn5a08/PPpneb7RqCiNTDCIMXVPU/zu6fnTDG0XDGniGQVfVxVe2jqn2aNw/mNeVW+wFW0ZyN7EQHvDUEt9pvyW26dDH1wQsRaNvWmITOPBPGjzdjT/G0Cqt9WtLG4sXmbzYLBGdZyKeA2ap6r+vQWOBc5/9zMRFDQ6HqByosooOvQIiq/RbLkCGmPnjRoIExCU2eDGPGVLiaxtMqysvh/fftuIIlDSxZYv5mucnoMExI4WNEZKqzDcAsKtJPROZjloa8M4RrAd4fqJ9AiESM2m/9yC1g6sHQocYcFNUU3OYhr3oST6sAM1fIjitYUk5UQ0ihQKix26mqfkHFAuuxHFvT/L3o0sW4lLonEi2kI8fxAWaowhQnEoHHHoMLL0xFKSy1lZEj4YQTjKY5fz507mw0B79Ow5Ahxn05aqL0wz2ucMIJthNiCZnFi6FFC9N7SRGhehmlCy+1fxEdaMgWWvIzIuaZ3XSTFQYWb4qKjFko1jzklzZWq4iHHVewpIQUu5xCLRUIXh/oYoyn0UHNFnHWWRUDghZLGIwcaerUmWea8YImcebdl5fDt9+mr2yWHGHxYisQ/Ij9QPceYATC2PsWJuzxWSzVwa1VHH98fG1hwQI7lmAJkbIy+OGHlAfwzKrgdslSVORq+EvbQQOBhd5zESyWMEk0rlBWZscSLCGydCls2wYdO6b0MrVWQ6hCQYFxIl+wINMlseQAUbNlJOKfxo4lWGpKNBrDJX1Nuzar1AqE4HTubHxSLZY0MHIkdOrkf9zOkLfUBHe4lPIFxvLx+6GdUmqKrHsCYf78YLEJLJYQ6NXLfyzBzpC3VJfYIIydWEAJBSwo2YO77krdBMi6JRC6dIH162HNmkyXxFLH8It0Gm/ms50hb6kusdEYOrKQxbRHyaO0NHXu9HVLIHTubP5aPd0SIvEinfrNURBJuYegpQ4TG42hIwtZQIV9cvbs1GgJdVMg2HEES0h4rZ8QG+k06gLdrVuFUFA1H61fKIugayvYNRhyk8qmRqUTC1hIxYByeXlqHBbqlkBo3964fVgNwRISw4b5u5bGehEtXlw5nIpfiOygayvYNRhylyFDjJYJ0IKVNGJTJYEAqWnmar1AqNSD+lM9trRubwWCJRSGD4dx4/yPu72IgobIDqJxJJPOUjcpKoLu3c3/nTAup26TUaocFmq1QPDqQX32U2d++swKBEvNiDbI8RzW3B9l0BDZZ50VTHDYNRgsTz5pHBOiC3+5NYRUOSzUWoHg14OaU96ZnVfMZ2KxdT21VJ8gS23WqwerV5uGfvXqChXfi2iI7B9+iC84ohpHIgFjleC6T1ERXHst7J2/kDLyWEK7hKHaa0qtDV3h98HOpzM7s5F/37WCojdbp79gllrHxImmPs2bZ3r8Q4YEW2qzvNyYlFSNMKjp9Be3xuEV4t0rnaVuM3IkrJ64gDUT2rLf3vUThmqvKbVWIPh9sPMwX8r22fMBKxAs8Rk+3GiaW7aY+jR1KrzxhomEEo+8PBNaJoq7LlZXOLjNAPFiJdn5DbnFbhsWQlFHJscZzwqLWmsy8lvFaj7G9fTAxlantsTHz+xYWgq//OJ/XiTi3+CLwJ57Jg6R7cbLDFCdld0sdZSFC1Me1C5KKBqCiIwGTgJWquo+zr4RwEXAKifZjar6ThjXA/8e1A+0pZT6nNjVCoRcx8sU5G5Ig4wTRNjOvkynK3Npx/fsnr+Snepvo2RzGWtpykpasJCOfEsvVtAaVWje3ITIHjTIODp4mX2iiJj5C08+WbWRT3ZlN0sdZMMGM0AVL2hWiIRlMnoGeBh4Nmb/fao6KqRrVCLagxo1ynhdlJebHlRhYYRfGnZgj412clou42UKGjvW1Jnowkl+ZsfdWcZpvM7J/JeDmUQjNu04VlZvJ7Zsr08pQmPWE6GitV/KHrzLADY0PAlKj2fIkIKES2+qViyV60WlEO+W3CMazj9NGkIoJiNV/QxYG0ZeyRC7SM6ZZ5rfzQ/tYt0wcpgnn4Tbb0/sw+82OwrlDOBt3uc4ltGGBxlCK1bwNOczkDHsw3SO6PkLkc0bmfHpWto2XEN9ttKCnzmCTxnC/UykiLN4iaGfnwxt2lD0+jXcfsGihEvgWjdSiy/RcP5p0hBQ1VA2oB0ww/V7BLAEmAaMBpr4nHcxMBmY3LZtWw2Fq69WLSxULSsLJz9LrWHYMNVIRNWIAu+tb1+TtrhYtWGDcj2V13QWe6uCLmV3HcYI7crsSufk5akOGlT5Og0bmv3R4w0bqo68qVT13XdVTztNNT9fNRLRn08arD13WRi3TL17Z+Z5WbKcv//dVJBff42bDJisYbTjYWSi3gKhJRDBaCG3A6MT5dE7rK/iscfMrX3/fTj5WWoFxcWmUY7X8EYb72HDVPXrr/XH3Q9SBZ1Bdx3Ii5rPVs9zCgpM/rHXGzjQNOYDB1Y9rsuWqQ4ZolpQoNvy6unfuUEbstGzPG5hY7Hs4IILVFu2TJgsLIGQMi8jVf1ZVctUtRx4AjgoVdeqQteu5u/cuWm7pCXzBBkkBigs30Szv1+NHnwwbVjKwpuf5vYzpzG/9yC69qjnuQpaeTm8+27lfe41lr3W8Z74w+4M+vl+ju+8iI9ans0N3MFsunEc71cuT5JupDbgXQ4xb176zEWQUg2htev/q4CXEuURmoawfLnpej34YDj5WWoFBxyQWDs4gMk6n46qoGP3uETPP3W9HnBARQ+/uNhYG73ObdjQQwvwIWpSEqnQAo6MfK4z6a4Kei9XagPZog0bOtpKQLzyTTYPSy2iZUvVwYMTJiObTEbAGGA5sA1YClwAPAdMx4whjHULCL8tNIFQXq7auLHqpZeGk58lrURNMe6GOggDB1bY9Ktu5XoZD2kJ9fV79tQj+ETz8qo2rD16VOzzMu0MHBis/H6mq13rb9Y39vyrKuj3u+6rU15dkNRz8cs3GWFlqSWsX29e7l13JUyaVQIhrC00gaCqesghqkcdFV5+lrRQkx6wX4NZwBZ9nrNVQf/HidqU1b4ahJ8wSGbwd+DAAELl7bdVmzY12/jxgZ5NoHwtdYevvjIv9803EyYNSyDU2pnKCdl7b5gzJ9OlsCRBTUM+e83ubcZqxtGXP/AiN3I7v2Msa2nmm0f0ul4EjSEUJPLpxKYD4KuvoFUrOO44eOyxGudrPa3rGNEx0DQGrqq7AqFbN1ixwqyxbKkVhBHy2T035Xf7LGJK/SL6MJmz5BXu4EYkLy9uVFLwDokCwQd//cKqRFm3zlno5tmOUFwM/fvDpZfCbbfFlUiJ8l21yg4w1ynmzjUvPE2T0qCuCwQw6xhaagVh9YCLimDMiLn8d90RtG20jgWPf8zaY8+gSRPYdVdo0cI/VHVenlmYpCYxhIYMMcIjHju0nlm7wJtvwp/+BLfcAldf7fsQEuX74492RbU6xbx5ZhXI+vXTdsm6KxD23tv8nTPHuullOdH3Ey+EQ14eNG0a8D3OnAlHHglbt8Inn/Dq0kP48kujLK5bBytX+guewkJ44gnvGfDRkBeJcJuu4rFD68nPh6efhr/+Fe67Dy67zLOA7ny9BJqqXVGtTjF3boULfboIYyAirC3UQeXt21Xr19fPD7nGuullMbGDyH5bfr5xB034HqdPV91tN9XWrVVnzUo4WS2V9aK4WLVJkyQGqcvLVa+91hy46irz2yfftm3987QDzHWAsjJTIa+6KlBy7KByAiIRNrfpwi+TZtt1abMUr0HkWPLyzNoEIqZHHfc9Ll5sBmjr14dPP4Vu3eKOS7hDVSerBQShqAiOP97f7h87SD1xkjDo+zsZ0/wKuO8+lp1/s2++u+3mf107wFw7cVsyLvv9MlPB070SUhhSJawtVA1BVYv3PEPn0cn2orKUeG6UYHrXAweq9usXwN1y+XLVjh3NSdOn77hGoslqqY4hFHTuQGVNqVz/xcWqoO8cf7/vs/Obc2FDYdQ+YjXlfvKhKujocz4KdD5WQ0jMtG3d6MAiCiipcsz2ojJPomUqW7QwA6mTJ8cfbP5pzgbjqbN8ObzzDuyzz47j8TxzqrsUZTJjUkEWuqmqKQmX8k/+wykc//5VzL3zjSr5xhtgtiuq1S68NOVOasL33/Zq1/RaMsKQKmFtYWsIDx7yoiroPkyzvagsJP7MYtNbqlcvfg+/nmzT71odZwYZ3nuvyjXCnt1b3Ylz8QLh+WlKDdikX1KkpZFCz4L6RVy142O1C6/3fz9X6C800jwpD2TJwM5UTszrt3yrCnoGL4fSGFjCJWh00njbI/lXmH+eeML3OmE1nKkKHRHPrLUbK/WHgo5moHzRIs8yxY24asl6vN7/O/TXyRyww3Sa6L1agZCAYcNUmzXYpGWI3sLISr1O24vKHoKsX+C3/bXeo+afq65KGP8ojIYzVaEjEo0HXHXSPBOba//9VTdtqt5FLFmL1/tfSHt9kYGVOhzx2iwrEOLg7sktol2lBxuJxO1MWjJA167JC4OjGa/bJaI6YIAOv3l7FddVEROoLswec9AB6qDB+aLpunb1F4o7NI933jE3NWiQrzuqpXYSq3kWsEXLEB3O8MBaqBUIcXD35N7mBP2W/UPpyVlSQ6KxhNitDT/oKprpTOmuzzy4Ia7ZqbAwPG0wiGdP0DEGv/kXfucVF6uO2fd2VdBne/7DmobqGNH6AKrdmaEKOpAXA2uhViDEwd2Tu5uhuoUCjbCtSk/Okh0kM5aQz1adwCH6C410b5mjbdsmntTm7llVN7R2onI2bGg0zyBjDPHyiURU9967ctl2CA/K9TVO1e3kad+Cz6zZs44Rnch4Cq+rgh7AZF8tNBYrEOLg7sn9iWdUYccauda7KDvxGvjNz6/aYP6Dq1SpcBQIMv4Q7VmFsbhMvAHqoGMMyYxFxAqPRvyi8+ikP9BG2zRY7TtWUh2BZ8k8Aweq3iS3qYLuxK+eWqgXViDEwf0R9eIbVdDTeSWhHc6SWbwGfvv1q/ggTuU1VdD7uSKweSm6de0anoeQ3wB10DGGZCbLeQmPXnyjpdTTsfxWB55VMZ5gV1Or/RQXq74UOVsXs1dS9dQKhAREP46Gslm3k6cjGWY/jlpIVLi3Y5GuZxct5mCtR2lSwiAvT+OalsIaVwo6eziZWcZ+wuMK7lcFvXvPByo9pzAEniWz/NSqp76Xd0JSbtJZJRCA0cBKKq+p3BT4EJjv/G2SKJ/QQ1c4PbklBV10UptT7EdRSxlx8zadkHeYrmcXbcuSpLWDhg0TezKFUfWCNspB0kXrrn9wvHL9Hyfp1rz6qt98Y1dTqyts365aWKjLBl2dlJt0WAIhrNAVzwD9Y/ZdD4xX1c7AeOd3WikqgjFjYK+T9uWgwumBYtlbso/hBXdyaPkEniv6J81778XOOwc7zx0iolev8ENYxBIkTEWQdO++a9Y1ePllE67bG+EvDZ5Gm+4G55zDzG9KMP2wqtgwLbWIxYuhpITd+3ZnzBgTtmXMmGDrcIRCGFLFCCjaUVlDmAu0dv5vDcxNlEfYGsIORoww3Sc7qafWEO0hn9N1km6XiK7qZ+woiSayiRjzUGzPKp0mlaCT4NzzENq2NX8PPDCxC24lE8K776qC3iPXxE1vHSlqCf/9r3lpSVZIsslkpN4CYb3rf3H/jjnvYmAyMLlt27ZJPYTAvGYGI/Xrr1OTvyVUouM/jfhV59FJv2dPbd1gnQ4enNg9NV7jno2xf4KuBxHdohFg3YLuyfyLtQzRQ/nC85yCAjM4bz2PagF33GFe2vr1SZ1WqwSC83tdojxSpiHMnWtudfTo1ORvCQ13T/4xTEN3BJ8oGM0gXsMZidQsyFy6qU4spyZNKjfsAweq7swvuoh2Op+O2pCNVTSmevWs51Gt4ZxzVPfYI+nTaoNAyB6TkTNQo//3f6nJ3xIa0cHRYzHx4O/h6sCNZdeumS59ciRaDyKI2ah1a/P7CD5RBX2Iy6oIhHSYySwh0bu36nHHJX1aWAIhleshjAXOdf4/F/hvCq8Vn0gEundn/RfT7drKWc68edBQN/IEFzGPztzCrYHOy8uDAw5IceFCJtF6EPGIrhi3cqVZ+e0zjuQ+ruRyHuEoPt6Rzi//Hes5W7KH8nKYPRu6d89YEUIRCCIyBigGuorIUhG5ALgT6Cci84G+zu+M8W3ZvpR8NZ2XX4YpU+CVV4wnx/DhmSxV7pFocZkuXeBObmAvvmcwoymhwY5jIka2e1EbF4WJt3hPUMrLK/K4idtZSAf+xZ8poASR+OfFeh4ls/CPJQX88IOR8hkUCKGZjMLYwjQZuafw9+2rel3+KFXQZqyyqnOGCDKTdsajn6niPRu5YUPdMbCcTQPD1SWM9SBAtVWrimfSlw9UQe/Iv1l79Ag+Ac7Ocs4C3n7bPPwvvkj6VLJtDCGMLSyB4OW5cRzvqYIeycdVbKx20k7qCeT2uWmTaqdOuqZxe92twUbfRt89MNy3b+32oIn1fPIaJwjSsLufyaftztGySL5OfWFGaBPlLGngnnvMQ1+zJulTrUDwwa9yt2aZKujlPOjZw7KklkAzaYcONTvGjQvkDVRXerVeAs5930k32CtXqjZrpnrooTr8lrKEGlW8d2M7TGnkvPO0tFmragUntALBB//KXa6raKZPcEGVY5GI7QWlmkQB3QZ2m2pexIUXBsov13q1Sc+h+Pe/TcJHH00oXBO9G9thSg9Ld++jH+cdU60OjhUIPsSr3B9yrH5Nb9sLSjFeIZjjBXSLSJnOa1ak2rx5YHU5F2P3JDWHorxc9dhjVXfZRfWnn+Lmm8j91XaYUk/xF9t1Ew30Xq6sVgcnLIGQSrfTjBDPc+NberEv08lnW6X9qjbWS1gMH14Rh8ftzdWwofEE8uKSek/Sec1EGDUKmjYNdJ14Lpt1NXZPNDZXoPg2IvDYY8a/9Npr4+Y7ZEh8b6fycuuimmpe/vtCGrKFaexX5Vg6XYTrnEAYMsS/4fmWXhSwlW7MrrQ/rOBmuc5Ep03fvLmisY76y7/0EgwcWDWg214NVnJP3nVw5JFwzjmBrxVP8Nv36dCpE1xzDTz/PHzxhW+yoiJo0cI/G9thSj31Zk8D8BQI6ezg1DmB4BdJsqAAZuT3AqAX31Y6R8T0Zq3vdc144AHYssX7WEmJEQzjx8OZZxpf9zPPhK+PHEqDsk3w6KPEdZyPIZ7gr41zEqpDoHkDN9wAe+4Jl18OZWW+eR15pP/jtwI29RzccDpl5DGLqnMQ0vr8w7A7hbWlYh6C2946/ObtupGGej9DfG3PtdFLJVtIZiUwVVX96CNz4KabqnW9bAxWly6S8rB69VWT6JFHfPPLtUH6bGPNEb/XOdK12s8fO6hcPX7Zp0hnNT9Cu3b1D6NsP4DqkcxKYFpaalaT79BBdfPmal8zm4LVpYukG+/oAHPjxsYlVb0H/nNZwGacDh10evczqv38rUCoLpdeqrrLLjrorLKc81JJNUk1VNFJOG+/nbHy1laq5WE1c6Zqfr7qRRfF1S5yUcBmnF9/NS/i1lur/fytQKgujz+uCvq7HguSM29YAhGol7l8uerOO6ueeGLGylmbSdo0F+Xqq7VcRA8v+MpqxtlEcbF5AW++We0swhIIdW5QOSG9zMDyUY2/tV4qKWDkyKoDx+PHm/07uOEGyreUcBX32UBq1aDaHlbDhrG+oCV3ll4JVPXZtRFQM8Q042HEflU9jNJOGFIlrC0tGsKWLaqRiC4990Y7iJYJJk1SBf1H/rW1PuREpqjJAPDItk+qgp7OK0lrxl7jDpYQuOwyozGXlVU7C6yGUE0KC6F7d/ZY+W2gBdEtIVJezsbz/8oKWjF8+81ozFyFUaOsphAEP9fqIHV3TtF5fMd+3MV1FFBS6Vg87cJvwqENHx8C06fDPvvUPBZ6GIQhVcLa0qIhqKr+6U+qLVuqlpfbQbR08swzqqDn8owdzA8Bv7obrydfXKx6YoFZjW4odwfSLqxLagopLzfeX3/+c42ywQ4q14CHHjK3/sMP6blejlNcrHreab/oqvxW+nXkYBXK7GB+ivAK/S5iQmhHG+5hw1TfyTtR17OL7sbKhCa7XIwblTa+/948yDhzRIIQlkDIAh0lAxx0kPn71VeZLUcOEDU1dHv9NnbbvoK/lD2IX7Wzg/k1wyt0CJj/Z86Eo48272PkSGj9wj00kk080nyk98C/i1yMG5U2pkwxf7Nk/deUCwQRWSIi00VkqohMTvX1ArH//lCvHnz9daZLUqeJNlC7b57PldzH05zH1xzkmz5XQk6kinihQ8B4EUXHaXoO7EbkL5dw5trHGDNsdtxxBxs3KoVMmWIeYjZ4GJG+WEZHq2pPVe2TpuvFp6DACAWrIaSUaAN1N9dSSgE3cIdnOjuYHw7xevJRKrmWDh8OjRoFioaa63GjUsHEiTDlqSksqN+NQRc0zAqHitw0GYExG02ebHReS0qYNw8O0885hTe5k+v5mVaVjhcWxpmrYEmaeD35KJVMPM2bw803w1tvwbhxvufUxKvJ4k3UlNrqpyl8WXJA9nhthTEQEW8DFgNTgG+Aiz2OXwxMBia3bdu2RgMrSfH002YwZ9as9F0zxxh0VplO4kD9kT20AZuqDEjahVfCJZ43kG9MqZIS1fbtVXv2TOgHbz3ywiH6nlrxkyroEO6rsdcWtWhQ+XBVPQA4AbhMRI6IEUiPq2ofVe3TvHnzNBTHwQ4sp5yRPV7hIL7mZm5jCw2rHLcLr4RLtCdfUOCfpl49WL3aFTL72wK47TaYOtVMMkiQf+AFeiy+RE2p0TD8U6gYUN68GW65JVMlS4PJSFWXOX9XAm9AnFHFdNK1q7Gf2oHluASKue9FaSmdR9/AzPz9eQ7vhW9UrYdK2IwcCZ98Aj16VDYf5eVBfr4RwuPGVZ5cNmLOQOjZ05iPtm7NVNFzhuhYT1QgTKVnpeMffZQ501FKBYKI7CQiO0f/B44DZqTymoGJRKBPH6shxKFGs1MffhiWLOG/h9+DSsQzSV6eWTGzWgLH4ktREcyYARMmmFXqeveGY44xAmHbtoqB5+gM8TvvzuNGuQMWLeLpQx5n4sQadAQsCYmO9RzAFObTiV/ZpdLx8vIMztoPw+7ktwEdgO+cbSZwU7z0aZuYFuXaa1Xr19eJn5bYGC0x1Gh26po1ZvZl//5x88nPVy0sDLjIi6VGxJtcZrZyHc/RuoIW2jjyi9arZ99Lqoh+E4topy9xZigT/rAzlUPgtddUQX9TMMlW/hhqNDv1qqtMomnTVNU7JHZBgWq9etUUOJakSRQyG1QPxAQeHMYI+15SzB3XrFEFvZY7Q5m1H5ZAyF23U+CbemZU7IDSL9EYNTrXA61Ve3bqokXGXHT++bDvvoB3SOwjjoDt272zsGGYwyeIS+rXHMSrnM5QRtGclVWO2/cSHtcfX3VA2U2mJvzltEAYNWYPlrAXhzGhyrFcr/zVnp16ww3GleVvf6u0O9ZDZc0aGw4hncSbXObmJm6nAVu4mduqHLPvJUQmm6ANcxv08jycqQl/OS0Q5s2DCRzmCITKrVOuV/5qzU6dONGMPA8dCrvvHjd/Gw4h/bRvn1hLmE8XnuICLuEx2rOoyvGmTSv+twPPNWDSJOjYkfOv2S27JvyFYXcKa0v3GMLAgap/kUdUQduxKP4EnhwkqUXXy8tVDztMS5u21HNP+zXhAL0NqZw+vCKgxttas0w30UCf4w9VjhUWmvzirctsqYxnOPI99tjRwIQx4Q87qFxziotVDy6cqgr6B56zjZIHgSvrf/6jCnpZ/X8FbiSSEjiWahFk9rLXDPI75AZV0P35tsrxggIjGLzOFVHt189+O1G8BGfHwqXmx/33h3YdKxBqSLSh273ldl3PLvool9hGqbps3aqb9+yss6SbRtiWVI/fhkNILYndTatuffuq/v6odbqapvoO/ZM61/3Oc/0beuIJI1xjn83vMZ2nvx5YHFp9twKhBsRK7fc4Tr9jX23d2jZK8fBdictZcOhE3vJsHOwiKpkjiLupl7vjAQeoXs09qqBH8Em1hUKufkvDhnkLA1C9g+u0lHpawJbQBKcVCNXES4UexggtQ7R1g3U5V4GDLpzuZzO+/dr1qs2a6VeNjlYoj9vIWNLPwIEVJrkgW3TsbOBA1YayWX9kD53AIXHfbby8crEjkMhM9xFH6Vf0CVVwhiUQcs7LyGsRkQkcRh5Kz5KJOeVqGjQ0hddKXNH5GvXvvRPWrOHpHqMA8byO9RrKHEHdTaNEPciGDAEaNGAEIziUYn7L/3akKSgIlmeueurFW6gojzL6MJlJHLxjXza5uOecQPCacDWJgykjj0N0Qs5U4HiN/K23wnHHVbgR+lXwPfmBy7ffx/92/SOjp/ovAWgXUckcfmsZ5Oeb6SJ+7o7R815pcB5z6cLfuZF8KaNhQ7juOrOmTsOqAWwrkasdgXiTOrsxm53ZWEkgZJPgzDmB4OX/vpGdmcIBHMUnOVOB4/ViVOHDDyu0Bb8Kfhs3A/DXX26ntNQ7r8JCu4hKpvGaKf755/DZZ5X3xS5SNHIkfPBRPu8cdjv7MJMHDnqB8ePhhBNMndhzT2jZEsRbMczZjkC8OTYHYYJpfuUK+pxVgjMMu1NYW6bGEED1Lq7RUurppI83pbwM2UDQwcaGDY0bYawduhffqILewXVxz+/bN9N3aqkx5eVmEGivvXTkjSVVxpLy801cKus+bIg3hvAYF+taGqtQZscQsgE/Fbq44Gjqs42DtlUNY1EXadYsWLqSElNt69d371VGMZTVNOMOboh7/jff2JmstR4RuOMO+P57frnnX1XMjNu3m2/o2GPtkqjg3cZEOYRivuIgok2viJlBnjWEIVXC2jIxDyHq/z5p3C9aHonoG91uqPNhsIcNM5OLgnqLtGpleoHR3wN4SxX0ch5MyuMkl3uNtZ7ycp3e4hj9mebaiF+sR1EAiouNhhzVnBqzVssQvYlbQ/82sG6n4TJsmGpx3iH6JUV1ugFLduaqSGV/6gjbdAbddS6dtR6lSbsi5rJvem3nnK4mPPYtjPR8t9a1uCruiYEn8LYq6JF8HPq3EZZAyDmTkRdRj5vx5UdzIF/TiF/rbBjseIPJXqhCWVnF78GMpgezuJ472UZ9/xN9yCYXO0tybOt1EP/hVIYyimasrnQsqwZGswi3Q8Zv+Jxt5FcaUHaTDd+GFQhUNJIfczT5lHE4X+w4lg0vKUziucQlohG/8jeG8TmH8wan+KaLeK+YCWSXi50lOYYMgdsKb2MnNnEDd1Q6lqseRYlwexwdzhd8Q2+24O2vmw3fRsoFgoj0F5G5IrJARK5P9fWqQ7SR/JJD2Uo9jubjHcey4SWFSTyXOBHjRuh3/BruoRU/M5T4k9D22MOGtq6LFBXBb6/txvOR87iMR9iTHxAxHYBddzUdp7qkTYdBdGJgASUcxFd8zm9802bDt5FSgSAiEeAR4ASgOzBIRLqn8prVIdpIbqEhEyniGD7acSwbXlKYxJu52qAB7LeftwaxO8sYyijGMJCvXJNqYikshFtuqcZaCpZawciRsO9rw8nLgzsbjCQvz3Sali/3n+mey0Q9jg4vmEwBW/mCw33TZsO3kWoN4SBggaouUtWtwEvAySm+ZtK4G8lx9OUApuywkWbDSwoTP7fb6CxVv5XMbuUWIpRxI39PONP1wgvjX8NOUqvdHPD7tqw+6zLO2vIMnctm76gv0XG3u+4yM92tu7Fh5EgYPdiYoTftdyg9eph2JSu/jTBGpv024HTgSdfvc4CHY9JcDEwGJrdt27Z6Q+whEA3edrAYT4qz5cU66WUUxS/stFcwtP2YqmWI3s1Qbdq0In2i0NU2tHXd5aJTVukGdtbXONW6GwfhxBNV9957x8+wvw1qg9tpEIHg3jLpdqpqXsqgM7frukgz/bTdOTnZgFV1Sy3XD+irq2mquzdYm5PPxFKVAw5QvYWRqqAHMsm6G8ejrEy1cWPVCy9M2SXCEgipNhktA/Z0/W7j7MtKiorgxZcjND7reI7Y/D5FB5VnukhpJ9akdDzv049x3FlvGBde0yTzKq0lK+jSBR6Qq1hJc+7keohZkzyWuuatF0vc9aWnTYP16+E3/gPKWUMYUsVvA/KBRUB7oD7wHdDDL32mNYQdPPus6dZMnhx4vYC6RnGx6tlnbtMFhT10eaOOOvGzUt90ufh8cp2oJvlXHlAF7csHCbWEbPm8wybh+tKjRpkDP/6YsjJQG0xGppwMAOYBC4Gb4qXNGoHw88+qoOOOvi23FxJ/4glz46+95nnYLrSe2wwbptq4QYkuZi/9mt4abxEdkR1rytcp4s3832EmGzBAtUuXlJaj1giEZLasEQiq+uvevXVC3mG5aw/99VcTxOjQQ020yxgCfQiWOk9xseojB/9bFfQ0XvUVCJFI3awT8daszstT/cOZW1UbNVK99NKUliMsgWBnKvvwQd4JHFxeTGPWVTlW2+yhce2bfowaBStWwD/+4RnwPl4IjNr2fCzVp6gI/jLhDywo7MHt3ESE7Z7pmjfPApfKFBBv5n95ORR89xVs3GgmaNQCrEDw4X/bTyBCOcfzfpVjtWn2ctBlMiuxdCncc4+JY+zzFSf6EGrL87GEQCTCmwf+na7M4zyeqXJYBI4+Ov3FSgfxZv7n5UH/+uNr1QOwAsGH0p4Hs4KWnMIbVY7VltnL8ZbJjBu077rrTMK77vLNO9GHUBuejyU8Dr/rt0zMO4QRjKCQyqpjgwZ1a3Knm3gz/wsL4fj8j6BXL2jaNL0FqyZWIPhwxVUR3s4/mQG8QwEllY7VltnL1TLrTJgAL74I11wD7dr55p3oQ6gNz8cSHkWHCLP+eAdtWMbl8ghQeQYuVMNsWQuIN/P/hiGb2WVmMRxzTGYLmQxhDESEtWXToLKq6rNnv6sK+lv5X630okm0TGaVx11WZnbusYfqxo0J8496GdmlEy1R1hX111/rN9Ej91+3ww05F7zRPGcev/++ueF330359bFeRmmgtFS37bSLftR+cK0Mv+AVhsLtARF1A4xW5uF7jVYFnTfihbj5uuce9O1r1lyujc/HkgKmTDEV7KabVDXHvdGuvlq1fv1AnauaYgVCujj7bNVmzVS3bct0SZImyMcY7b3twgZdTkv9gkO1YYNy395bLvT2LJVJevLhwIGmUixfntAts04vu9mtm+ktpQErENLFq6+ax/Txx5kuSbWIZ9ZxC4w7uVYVtDdf+/becrq3l6NUqwMwf75ZhPuyy5I3W9YVFi0yN3jffWm5nBUIKSbaKzps/1+1NFKoy0+/PNNF8iRI7y1eZFMR1U7M01Lq6VOcH7f3ltO9vRykRh2ASy5RrVdPrzhpYSCzZZ3jkUfMTc6dm5bLWYGQQmJ7Ra9xmv5MCx1xc3aZjWpqvjG9t3J9iwH6C420Jcvj9t5ytreXo9SoA7BsmWqDBrqy/x9zU6s86STVjh09Z/mngrAEgnU7jcHLd/95/kALVjLlnvFZ4y5X7TkGLrp0gVPlTU7kHYYzkp9pVen4okWVXQTt3IPcokaTD3ffHa64gubvv8A950zLrcWSSkpg/Hg44QTPWf5ZTRhSJawtGzQEr15RfUp0LY31Wc7JGrNIGOabSeN/1R+ljU5lP42wzTcvrzGHnOrt5ShBvdR8WbvWrANw0kmBFoSpM5Fz33vPPKR33knbJbEmo9TgZxZ5nAv1Fxrpob02ZbqIqhqS+eaaa1RBjyn4wvfD9/NKsnMP6j6hdADuuMOc8PnncZPVKe+1Sy81hd+8OW2XtAIhRfj1io7kY1XQBw8dkxU9mRr33qZPN54gF1yw436aNPEXCG6twy6NmTvUuAOwaZNq69aqhx3ma0+vU5pnWZm539NOS+tlrUBIEcXFqoWFVSumUKY/0EbH1TteI5FgPZlUCo4afUTl5aqHH67atKnqqlU7dttBY4sXNe0ALLzmUVXQIR3Gep5fp7zXvvzSFPz559N6WSsQUsSwYabj7FU5RzBMyxDdi8W+jXD042nVSgMLjpqUtVq9t6efNic8+WSl3TXWOiyWGIYNU92lwVadQxedQxctkNIqdbROdUSGDlWtV0913bq0XtYKhBQQr9cNqnvyvZYh+jdu9mwwe/SIf35NVGA/bSPp3tvy5cY2dNhhRr0NeP+1TnW3ZBx3fTqBt1VBr+TeKvUpXkcEVNu2rSV1r7zcuJr275/2S2e9QABGAMuAqc42INE5mRYI8VTX6PYWA3QZrT29chINzFZXBQ51wO2001QLClRnz457LTtobKkpsd/TO/TXdeyqu7GyyphUvI6USC2pg999Zwr8r3+l/dK1RSAMTeacTAuERKorqP6ON1VBf8ebCdOGoQKH2mt3wnC8uN8d1ZrZbLEkQ+z3tDezdCv5+ih/rvIdxHZ6aqWWesstpge1YkXaL20FQgpIpLqCaoRtupTd9X36VUsYJGuLT6S1BFanV6/WX3dqoVPkAM13tBvb+7ekEq/v6T6G6HbytKdMrfIdFBeb+hymdp02ystVO3Qw4X8zQG0RCEuAacBooIlPuouBycDktm3bpuZpBSSR6hrdbuB2VdB9+S5pgRC0lxPtpTdoECzPRI36yv7n6FbydT+m1r6el6VW4vU9NWatrqKZfpJ3lBZ/WdUNtbYNMEe/0z91LVYFXXDT6IyUIysEAjAOmOGxnQy0BCKYVdluB0Ynyi/TGoKqtw09EqlcKZuwRjfSUJ/m3MCCIBk7aBD1OalG/fXXVUFHMqz29bwstRqv7+mKev80P157rUr6bPR083PocH+nD/BX3UKBtmqwPiMad1YIhMAXgXbAjETpskEgqFa1offrV7VxfpDLtZR62pplcYVAJGLmqQSdrv/EE8G0lMCN+rJlqk2b6oyGfTSfrbWm52WpO1QZk/p8m27suK+ubNhWD91/YxWvuWzydPNz6Bg8uKKcEbbpClroK5yesXJmvUAAWrv+vwp4KdE52SIQYvFSYzuwQMsQvYPrPCvvLrskHpD1qmyx2kgym+eSmP36qTZsqFcOmJt1PS9LbjJsmGrfgs9UQe/guipjWdni6RZPOLm/peMxS+2ezBsZ07hrg0B4DpjujCGMdQsIvy1bBYKfGjuGs/RXdtLdWFnlWGFh/AocdLwiGQ2hSqN+//3m4GOPZV3Py5KbuOvhU5yvW8nX7szYoVH361d5gmcmPd2CuKGD6kucqatpqvUpyZjGnfUCoTpbtgoEv8a0K7N1O3l6D1cn3dAGrWzuLWqCCnStyZPNfIOTTtoRQyZbel6W3MDL9u6u981Ypatpqp9xuApllepyNtTJIG7ozflZS6mn93JlRjVuKxDSTGxjGt2e4U+6mUJtxU+evXY/1TFIZfNq9KO2y7iN+tq1qu3aqe65Z6VYRarZ0fOy1H38bO+tW1eu0+fzlCroeYzOmNbqN2gcxA19KHergnZj5o59kYhq167p/b6sQMgAXlFBO7BAt5Kvj3CpZ4Xxu6VElS0S8W/04zbqZWW69rCTdFtePf1T12Lb6FvSTjzzpDu+F5igkZ9zmK6imTZldaDOVJjECi63I0jfvt6BLiu2cp1LZ/2Mwz2Pp1MDtwIhg8RW+Ae5XLeTV8XHP57qmMim/8QT1evJf3i0mSPxVx60ZiFLRohnDvUye/Zgum4lX//NOYE6U2ERJGRGfr6JVRcvJP45/DuuFlFQkPpOmRUIGcbds2jMWl3Jbk5PoTyw2hu2TX/OHWa+wQsMqlSOdKvgltwmkTm0VauqDfEIhqmCnsTYhJ2pmuA2D7VtG2wcr6DADHb37m20hvr1zf7XOUXX0EQbsClhHj16hH8vbqxAyALcoa4vlCdVQf/Is0k17KHZ9L/+WksiDfRLirSALZ7qq518ZkkHQSaXFRebxjWarh6lOpX9dBmttTFrU9KBqc6Ez9hvJxoevyPztQzR27gxcB6p7JBZgZBlFE8o07nNinRdpKlecvJP6e2Nf/+9aqtWuqz+XtqCFb6VshY/XkstIhkXZ7eW3JMpuo2IPh/5U5XxspouMlVTN+/evSvn8RCXaQn1dziTBIl0nMoOmRUI2cicOaoNGujaQ07QgWeVJ12Jq1X5f/pJtVMn1V131aH9p/tWTBGjItf6BcwttYJkzKFuLfm1HreYE8aO3ZGHux5HIsbTLlmq4+Ydq9VE82jKat1IQx3NeZW+ryBCJVVYgZClvNX/IVXQq5yFQIKaj6q15sHq1cY4udNOqhMmBBokS6ZMFktNqJY5tLRUdd99tbRpS23XwF/bjScUvDpW1XHzdguhJ56oyOMWRqqC9mB6pXTxPJJSPTfBCoQspLhYdacGZfo6p+g2InoM43xV5djzkp5F/NNPqvvtZ0a8xo/fsTu2Z1ar48tbcpNp07Q0UqjvcnylCWuxjbRX3fXrWPXoEcys49eRiubRVNbqOnbVN/ldlQa/Xz9/oZDqb80KhCwkqlI24hedQXddTdMdvYh4g7qJ1NkqIdbnz1dt395oBh98oKqVe0V9+1Z4RcTzpLADzZZs5fY9H1UFvZp7fL+L2Lobr2NVUJBoTkHiraBA9Y7ITapUDX0fbfAzFQ3ACoQsxK2WtmehLqO1Lqeldmaugr8NMZE6m5fnqlAff6zavLlqs2aqX32lqvHNTbUtvrzFoqo68KxyfZXTdCv52oevAtXdeB2roGuex9tayErdlLeTvhI5K26Dn4loAFYgZCGx7nZ7M0t/prmuoIUeLJN8bYhBpsjv1KBMl1x2d8W8+DlzVDWxualfv+yLL2/JXYI6ThQXqzbLW6tLaKuLaFdpFnPUjBNbd4N0fmLdXWO3eJr6/Vyh28nTb8fMzrrwL1YgZCFejXNXZutC2utmCnXh9Y/vCDSX6LzYPD7nMPPjtNNUN2zYcW6iXlHfvjbKqSU7SNZxYvBg1QOZpFso0HEcoxFn6Ve/upvM4jpepp2CAhO23uv8HkzXbUT0g06Xpu4B1QArELIUr4rWtsFKXdjuGLPjuONUv/3W87zYhr0NP+gjXKql1NPVNNVh7f5dRaAE6RXZKKeWTBOv0+MOex3L4MGq58kzqqD3cmXclQeTdc5wm3Z69Ig3xlCuH3GUrqapfv3e6tQ8oBpiBUIW42lDLCtTffhh1V13NY+9Xz/Vp54ypp/t21VVdcCxJdqV2Xohj+tbDNDt5O0InNdKVniad+L1itxzD9wDzdmi5lpyhyDzAOI19O90GaIK+tiBT4YeDiaRhn4+o1VBxw54tMbPIVVYgVBbWbdOdeRI4yXkbrljuic/0EZv5SZty5IdH0s04J17uc2+feN/aHbugSUbCDoPwNeMuW2b0a4jEdWxY+NeK8igbtCYRm34QTfILrqh5xGmU5elhCUQxOSVHfTp00cnT56c6WKkB1WYPh2+/RYWLoSSEmjUiP9M2Ytb3z+YaaVdKVchLw8KC6FPH5g8GbZsMadWl4YNYfx4KCoK71YslkQMGgSvvALl5fHT5eXBmWfCmDEeBzduhGOOMd/Nhx/C4YdXqyzDh8OoUYm/pTzKeJ/jOTRvIg3nT4MOHap1vXQgIt+oap8aZ1QTaQKcAcwEyoE+McduABYAc4Hjg+SXExpCAGJ7OE88Ed5ym3bugSUTJBNLKG4zsHKlbm7bRTfl76znd/kioadSrEdTMuW4jRtVQR876KmUPJMwIRtMRkA3oCvwiVsgAN2B74ACoD2wEIgkys8KBG9qEocl6Q/OYkkRXrGJvDos8Vyhhw1T7VT4o86ls/7KTnqMfORpCo03YznIt3Qar6qCjo5cWCvG28ISCHk11C5mq+pcj0MnAy+paqmqLnY0hYNqcq3azMSJRmXu3dv8nTgxufPnzauZmchNXh506RJOXhZLMowcacyVffuaeuhFYSFccYX3sYkTjalnQUkbjuAzltCOt3QAAza/yqhRFd9VNN3mzRXfTXm5+T1rVuJv6VjG8QJ/YGLeISy97qHcMq+GIVWoqiE8DPzR9fsp4HSfcy8GJgOT27ZtmwLZmVmqFbQuhiAT14Judu6BJRuojjdQrKbcjFX6BYeqgo5kmA46q8wzXdBNRHVQiw91U95O+v2u++hX761J09OoOaTLZASMA2Z4bCe70lRbILi3umYyqlbQuiTzSaayWy8jSzaRbIgHL0+l+pToaM5TBZ3c6EjVxYsDhYLx2v/n+qO1LJKvuu++JnikT3mzMYR82gRCoEyqCoQbgBtcv98HDkmUT10TCIl6Km3bBq9UsT0qd0MfreSx67+6FwzPtgpssSSLv6ZcroNltG7K31m3NdxZ79j1Di1ks68w2Gefyt9Sc1mlr0TOMj+OPVZ1/foq1w5D008l2S4QelB5UHkROTioHMT3OplK5eV9FNvDykRgLYslHSTSuK85Y7G+nXeSKuiP7KE3cpvuzlJPzby4WPWy336v/2o9TDfl72w0g9tvN/MdkrxuNnxjWSEQgFOApUAp8DPwvuvYTRjvornACUHyq2sCIajtP1sqlcWS7fiNPQweXNFo/4ZPdRzH7PjAprGP/ptz9O78G/Szw643iXv1qvgATz9ddeZM32smiheWDW7cYQkEOzEthUycCMcea7wb4hF3Mo7FkkNMnAgPPGA867p0gSFDqk6ijKaZPx86dzZpHngAXn7ZNNNROjOP3/MmxzKeHpE57MEyJC8PmjSBffc1k9wGDkw44ax3b5gyJf7xTDdbWTExLeytrmkIqsF8r8F/bkA2D2RZLGGSrJ3e/W00aVK97ysIyURRzRRkg8ko7K0uCgRVU3HbtvWvrH6VKtsHsiyWsEjWTh/7bSTysIsGeaxOp8qOIViBEDrVCc2b7ZXQYgmLIHb6qEbQtavxnkskCGKFQk06VdkeQj4sgVCjmcqW4BQVwdChJrhcdJZmXp75PXRoVTvpAw+Y4FtelJSY4xZLXSHebPzycvjkEzMe9/LLMHculJUlzlOk4v9o3tEZy+6ZzUGIzrI+80wzZnDmmeb3yJHB86gN5Ge6ALnEyJFwwglVB8S8psYn+kDmz09tWS2WdNKlC0yd6h0NVQRWrQomBMCMGXfoYM758Ufv7yjaqUomLEVRUd2PEmwFQpoJWqnifSA2HpGlrjFkCIwd6+2Rl5eXOGy2O23//vDii6Yn/8MP3ulsp8obazJKI+4gd/36wXHH+Qe8GzLEBPryIl4AMIulNhLPpNqihb+2HIv72+jSxT+Inu1U+RDGQERYW10eVI7nFeE3QJXtA1kWS9h4zbQPMsHT69vIJccM7MS02kPQCWpeq5l5TcKp63ZMi8VNvO8nEjHfRc+e3t9GdHW0khJjJoquQDh0aN0aEA5rYpoVCGlg0KCqsyi9sDOWLRZvatKw16RTFWTmdDZgBUItItHUdzdNmsA772RnpbNYMkm6teXYtZezWbuwAqEWEXSB8SjRuQnZVuksllwhnpnKy7SbacISCNbLKA3E8xjyojoTZywWiz/JLmObqxND7TyENBB1qXPbQBNRnYkztYFt27axdOlSSkpKMl2UnKGwsJA2bdpQr169TBclI8SafqZONXMe4mnhuTox1AqENBE7S7lJE5g0CX791Tt9Xa10S5cuZeedd6Zdu3aIO7aAJSWoKmvWrGHp0qW0b98+08VJOxMnGmHgNv24w1eccIJ3pytXJ4Zak1EaKSoyHkSTJ8OHH8KJJ+bexJmSkhKaNWtmhUGaEBGaNWtWZzWyRKag6pp+cnViqBUIGSRXK50VBumlrj7v4cMrAt5NmQIvvQSHHgr77FMhGKpr+kk2GGVdoUYCQUTOEJGZIlIuIn1c+9uJyBYRmepsj9W8qHWPXK10yZDsYKAlN3CbgtwNvirMnAlHH20ERk3CV+RKhNNK1GSaM9AN6Ap8AvRx7W8HzEg2v7ocuiIeXtP16yqzZs0KnDZVCwTl5eXp/vvvr927d9f99ttPR40apWVlZXHPWbx4sb7wwgs1u3AALrjgAp0ZZ31fVdU33ngjYZpYknnutYF46ye4w1M88URuhK8gG9ZDUNXZqjq3RhLJUmlsYcwYqxmAdw+wurHsY2nQoAFTp05l5syZfPjhh7z77ruMTNDtW7JkCS+++GL1LxqQJ598ku7du8dN8+abbzJr1qyUlyWbiWcKilJSYnr0VgtPgjCkCt4awibgW+BT4Ddxzr0YmAxMbtu2bWrEpyVrCNpTDbKCVnXZaaedKv1euHChNm3aVMvLy3Xx4sV6+OGHa69evbRXr146YcIEVVU9+OCDdZdddtH9999f7733Xt90bhYvXqxdu3bVs88+W/fee2897bTTdNOmTaqqOm7cOO3Zs6fus88+ev7552tJSYmqqh555JH69ddf7yjnjTfeqPvtt58efPDBumLFCp0wYYI2adJE27Vrp/vvv78uWLBAH3jgAe3WrZvuu+++etZZZ3nec13UEBIFvHOvpVzXtXDStYQmMA6Y4bGd7EoTKxAKgGbO/72BH4FdEl0rV01GuUTQhumAA4J96NUhViCoqu666666YsUK3bRpk27ZskVVVefNm6fROvnxxx/riSeeuCO9Xzo3ixcvVkC/+OILVVU9//zz9Z577tEtW7ZomzZtdO7cuaqqes455+h9992nqpUFAqBjx45VVdVrrrlGb731VlVVPffcc/XVV1/dcZ3WrVvvECjr1q3zvOe6JhDiRTJNtFZ5XSQsgZDQZKSqfVV1H4/tv3HOKVXVNc7/3wALgTroRGlJFZmKZb9t2zYuuugi9t13X8444wxf00zQdHvuuSeHHXYYAH/84x/54osvmDt3Lu3bt6eLcxPnnnsun332WZVz69evz0knnQRA7969WbJkiec19ttvP/7whz/w/PPPk5+fG1OLog4ZBQX+aeqyp16qSInbqYg0F5GI838HoDOwKBXXstRN0umSu2jRIiKRCC1atOC+++6jZcuWfPfdd0yePJmtW7d6nhM0XazLZzIuoPXq1duRPhKJsH37ds90b7/9NpdddhlTpkzhwAMP9E1X1xg50qy13KNH5c6DHSOoPjV1Oz1FRJYChwBvi8j7zqEjgGkiMhV4DbhEVdfWqKSWnCJdLrmrVq3ikksu4fLLL0dE2LBhA61btyYvL4/nnnuOMmch35133plfXdPK/dLF8sMPP1BcXAzAiy++yOGHH07Xrl1ZsmQJCxYsAOC5557jyCOPDFxmd1nKy8v58ccfOfroo7nrrrvYsGEDGzdurNazqI0UFcGMGTBhAgwcmEPuoSmiRvqlqr4BvOGx/3Xg9ZrkbbHEhvsIK+Txli1b6NmzJ9u2bSM/P59zzjmH//u//wPgL3/5C6eddhrPPvss/fv3Z6eddgKMWSYSibD//vtz3nnn+aaLpWvXrjzyyCMMHjyY7t27c+mll1JYWMjTTz/NGWecwfbt2znwwAO55JJLApd/4MCBXHTRRTz44IO89NJLXHDBBWzYsAFV5YorrqBx48Y1e0C1kKBrlVviY8NfW9LK7Nmz6datW6aLkRaWLFnCSSedxIwZMzJdlJx67rmIDX9tsVgsllCxAsFiSRHt2rXLCu3AYgmKFQgWi8ViAaxAsFgsFouDFQgWi8ViAaxAsFgsFouDFQiWnKNRo0YA/PTTT5x++umh5Dlnzhx69uxJr169WLhwIYceeiiQviipFksYWIFgyVl23313XnvttVDyevPNNzn99NP59ttv6dixI19++SVgBYKldpEbkbAs2cmVV5qVzMOkZ0+4//5ASd0Tx5555hnGjh3L5s2bWbhwIaeccgp33303AB988AHDhw+ntLSUjh078vTTT+/QMgDeeecd7r//fiKRCOPHj+fjjz+mUaNGbNy4keuvv57Zs2fTs2dPzj33XK666qpw79diCRGrIVgsDlOnTuXll19m+vTpvPzyy/z444+sXr2a2267jXHjxjFlyhT69OnDvffeW+m8AQMGcMkll3DVVVfx8ccfVzp255138pvf/IapU6daYWDJeqyGYMkcAXvy6eLYY49l1113BaB79+58//33rF+/nlmzZu0IYb1161YOOeSQTBbTYkkZViBYLA4FruD60XDTqkq/fv0YM2ZMBktmsaQHazKyWOJQVFTEhAkTdoSq3rRpE/PmzQt8fmzYbIslm7ECwWKJQ/PmzXnmmWcYNGgQ++23H4cccghz5swJfL47bPZ9992XwpJaLDXHhr+2pBUbhjkz2Odet8mK8Ncico+IzBGRaSLyhog0dh27QUQWiMhcETm+pgW1WCwWS2qpqcnoQ2AfVd0PmAfcACAi3YGBQA+gP/DP6BrLFovFYslOaiQQVPUDVY2u6D0RaOP8fzLwkqqWqupiYAFwUE2uZak7ZJOZMhewz9sSlDAHlQcD7zr/7wH86Dq21NlnyXEKCwtZs2aNbaTShKqyZs0aCgsLM10USy0g4TwEERkHtPI4dJOq/tdJcxOwHXgh2QKIyMXAxQBt27ZN9nRLLaNNmzYsXbqUVatWZbooOUNhYSFt2rRJnNCS8yQUCKraN95xETkPOAk4Viu6fcuAPV3J2jj7vPJ/HHgcjJdR4iJbajP16tWjffv2mS6GxWLxoKZeRv2Ba4Hfqepm16GxwEARKRCR9kBn4KuaXMtisVgsqaWmoSseBgqAD0UEYKKqXqKqM0XkFWAWxpR0maqW1fBaFovFYkkhNRIIqtopzrHbgdtrkr/FYrFY0kdWzVQWkVXA9zXIYjdgdUjFCRNbruSw5UoOW67kqIvl2ktVm9e0AFklEGqKiEwOY/p22NhyJYctV3LYciWHLZc/NridxWKxWAArECwWi8XiUNcEwuOZLoAPtlzJYcuVHLZcyWHL5UOdGkOwWCwWS/WpaxqCxWKxWKqJFQgWi8ViAWqZQBCRM0RkpoiUi0ifmGMJF+QRkfYiMslJ97KI1E9ROV8WkanOtkREpvqkWyIi0510KV8qTkRGiMgyV9kG+KTr7zzHBSJyfRrK5bvQUky6lD+vRPfuhGN52Tk+SUTapaIcHtfdU0Q+FpFZzjcwxCPNUSKywfV+h6WpbHHfixgedJ7ZNBE5IA1l6up6DlNF5BcRuTImTVqel4iMFpGVIjLDta+piHwoIvOdv018zj3XSTNfRM5NRfkqoaq1ZgO6AV2BT4A+rv3dge8wYTTaAwuBiMf5rwADnf8fAy5NQ5n/AQzzObYE2C2Nz28EMDRBmojz/DoA9Z3n2j3F5ToOyHf+vwu4KxPPK8i9A38BHnP+Hwi8nKZ31xo4wPl/Z8yCVLFlOwp4K131Keh7AQZgQuMLUARMSnP5IsAKzOSttD8v4AjgAGCGa9/dwPXO/9d71XmgKbDI+dvE+b9JKstaqzQEVZ2tqnM9DiVckEdMsKVjgNecXf8Gfp/C4kaveSYwJpXXCZmDgAWqukhVtwIvYZ5vylD/hZbSTZB7PxlTd8DUpWOd95xSVHW5qk5x/v8VmE3tWWPkZOBZNUwEGotI6zRe/1hgoarWJApCtVHVz4C1Mbvd9civLToe+FBV16rqOswKlf1TVU6oZSajOARZkKcZsN7V8KRj0Z7fAD+r6nyf4wp8ICLfOOtCpIPLHbV9tI+amunFjdwLLcWS6ucV5N53pHHq0gZM3UobjpmqFzDJ4/AhIvKdiLwrIj3SVKRE7yXTdWog/p2yTDwvgJaqutz5fwXQ0iNN2p9bTaOdho4EWJAnGwhYzkHE1w4OV9VlItICEzF2jtObSEm5gEeBWzEf8K0Yc9bgmlwvjHJp8IWWQn9etQ0RaQS8Dlypqr/EHJ6CMYtsdMaH3sSEnk81WftenHHC3+Gs9x5Dpp5XJVRVRSQr/P+zTiBoggV5fAiyIM8ajKqa7/TsfBftCUKicopIPnAq0DtOHsucvytF5A2MyaJGH1LQ5yciTwBveRwKvLhRmOUS74WWYvMI/XnFEOTeo2mWOu94V0zdSjkiUg8jDF5Q1f/EHncLCFV9R0T+KSK7qWpKA7kFeC8pqVMBOQGYoqo/xx7I1PNy+FlEWqvqcsd8ttIjzTLMOEeUNpjx05RRV0xGCRfkcRqZj4HTnV3nAqnUOPoCc1R1qddBEdlJRHaO/o8ZWJ3hlTYsYuy2p/hc72ugsxiPrPoYdXtsisvlt9CSO006nleQex+LqTtg6tJHfgIsTJxxiqeA2ap6r0+aVtHxDBE5CPN9p1RYBXwvY4E/Od5GRcAGl7kk1fhq6Zl4Xi7c9civLXofOE5Emjjm3eOcfakj1SPsYW6YRmwpUAr8DLzvOnYTxkNkLnCCa/87wO7O/x0wgmIB8CpQkMKyPgNcErNvd+AdV1m+c7aZGNNJqp/fc8B0YBqmQraOLZfzewDGi2Vhmsq1AGMrnepsj8WWK13Py+vegb9hhBVAoVN3Fjh1qUOqn49z3cMxpr5pruc0ALgkWs+Ay51n8x1mcP7QNJTL873ElEuAR5xnOh2Xh2CKy7YTpoHf1bUv7c8LI5CWA9uc9usCzLjTeGA+MA5o6qTtAzzpOnewU9cWAOen+pnZ0BUWi8ViAeqOychisVgsNcQKBIvFYrEAViBYLBaLxcEKBIvFYrEAViBYLBaLxcEKBIvFYrEAViBYLBaLxeH/AZCr3qw8TYohAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "error = np.linalg.norm(y_data-np.dot(D,params))\n", "print(error)\n", "\n", "fig=plt.figure()\n", "ax=fig.add_subplot(111,xlim=[-11,11],ylim=[-21,21])\n", "x,y=poly_curve(params,x_data)\n", "\n", "ax.plot(x_data,y_data,'.b',markersize=15)\n", "ax.plot(x,y,'r') \n", "ax.legend(['Data points','line fit'])\n", "plt.title('Polynomial of Degree %d, ||error||=%f' %(len(params)-1,error),fontsize=16)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "anaconda-cloud": {}, "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.7.3" } }, "nbformat": 4, "nbformat_minor": 1 }