{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# EE16A Discussion 7B" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Question 2: 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": "code", "execution_count": 8, "metadata": { "collapsed": true }, "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\n", "\n", "\"\"\"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=[x_input[1], x_input[-1]]\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", " \n", "\"\"\"Function that defines a data matrix for some input data. 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", " \n", "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" ] }, { "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": 9, "metadata": { "collapsed": false }, "outputs": [ { "data": { "text/plain": [ "" ] }, "execution_count": 9, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD8CAYAAAB0IB+mAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJztnX+UVdWV57/7PYoqyyaG31SUarAjBEoErEqsGEeDFCjq\nip1klCLRoaUNQ2ISMmuMMTFCiMs1odW4nDE/lpCkHRMhZjqJrCQsA6jLkYZMChoNP+SHQiM2SomI\nRigtqD1/nHety+Pee879+e6P/Vnrrnrv3fvuO3V/7LvPPt+zNzEzBEEQhPxTqnUDBEEQhGQQgy8I\nglAQxOALgiAUBDH4giAIBUEMviAIQkEQgy8IglAQxOALgiAUBDH4giAIBUEMviAIQkEYUOsG2Bk2\nbBiPGTOm1s0QCsh//Afw2mtAX1//Z6USMHIk8KEP1a5dgmDCpk2bXmfm4brtUmXwx4wZg66urlo3\nQygYGzcC06efauwB9f7oUeDXvwba22vTNkEwgYj+3WS70CEdIhpNRE8R0XYi2kZECyufDyGiNUS0\nu/J3cNjfMmHjRmDOHKC1Vf3duDGJXxWyhv06mT0bOH7cebueHuCBB5JtmyDERRQe/gkA/52ZNxPR\nIACbiGgNgH8AsI6Zv0dEtwO4HcA3Ivg9VxYvBu69V928zMCWLcCqVcCttwJLlsT5y0KWqL5OvOjr\nA3bvTqZdghA3oT18Zj7IzJsrr98GsAPA2QCuBfBwZbOHAfx92N/yYvly4O67gWPH+m/ivj71/t57\nxdMXFBs3quvBfp14USoB48b1f1d6j0KWiVSlQ0RjAEwF8CcAI5n5YGXVqwBGRvlbQP8N2NQEfOEL\nwMmTzttJt1yweOAB9/CNEw0NwFe/qnoF06cDv/wlsHkz8Nhj6v3ixfG1VRCiJjKDT0R/A+BfAHyN\nmd+yr2OVdN/RnyKi+UTURURd3d3dxr9n3YArVwKvvuq9rXTLBYtdu8w9+8ZGFQ4ETu8VSO9RyCKR\nGHwiqoMy9r9g5l9XPn6NiJoq65sAHHL6LjM/xMxtzNw2fLhWVQTg1G65CfZuuVBsxo1T14MTREBz\nswrZXH89sG6dGvvx6hVI71HIElGodAjATwDsYObv21atAjC38nougMfD/pZF0G65ICxcqK4HJ844\nQ4VsurqAFSv6pZhevYK+PuCJJySuL2SDKDz8TwC4EcDlRLSlslwF4HsAZhDRbgAdlfeRYNotB4By\nWXXLRUctAOo6uPVWFa6xPH17+MbpOvHqFQDAkSMS1xeyAaWppm1bWxubTLyaM0fdXNUTZaopl4Ef\n/xi4+eaIGijkho0bVU9x927gvPOU5+/mFFgTs0xDiI2NKhwkToaQFES0iZnbdNtlMpeOV7ccULHY\nxkbgjjvE2AvOtLersE11+MZt2+pegRcS1xfSSiYNvtMNSKQ8+qYmNXPSGnAThChYskRdU9dfr+L1\ngz3mjff1Af/2b8m1TRBMyWRIx8JPt1wQokQXViyXVQ9TnA4hCUxDOpk2+IJQK0zi+hLLF5Ii1zF8\nQag1VlixXHbfRmL5QliiTuchBl8QArJkCfDhD7uvlxneQhjiSOchBl8QQjB1qrtyR2Z4C0FxSvJn\npfNYujS4py8GXxAMcOtae0mEZYa3EBSvbALvvhtcbi4GXxA0eHWt3TT6RMDYsbVrs5BtdNkEduwI\n5uWLwRcED7y61lamTEujP2FCv9FnVjelW8zVdDBOcvAXE10osK8voCCAmVOztLa2siCkiRkzmJX5\nPn0plZg7O9V2GzYwNzY6b9fYqNZbLFqkPiPq309jo/rcjul2Qv7YsKH/vLstdnMJoIsNbGzmPXzx\ngIS4WLwYWLvWfb1dhWOaQtmkx+BnOyGftLcDEye6rw8qCMi0wZcqREJcWAbXK45qv+lMUyibFkyX\nHPzC8uXRCwIya/DFAxLixKTmQl0d8PrrypC//roaqHXDSqG8f7/3g8HqMegeIKLvzz/t7cBtt/lL\n5a1jQLRNTA4TD0imtAsmWDmZdu1SHvvChWY1F/r6VMiHWRn7sFlK7D2GceOALVucc/WIvr84LFkC\nzJoVXc6wzBp88YCEKFi8WPUIjx9X19OWLcBvfgPU13t/r1QCenv739uvxaDG395NX7gQWLXKOVeP\n6PuLRXt7dM5rZkM6XlWIxAMSTHALC777LvDWW+7fK5fdDToRMHq0PoWyHaduepDKXIKgI6oi5j8l\nokNEtNX22XeI6JWqsoeRITMcBR06BZff2shWYZ0RI9wNPjMwfLgqrHLFFfqCKURKv+9Uv6E6B7+9\nsLogBCEqD/+fAVzp8Pn9zDylsvwhot8CIB6Q4I2JgstPbWRAee7r1gGXXWbWu9RVZgPU7+/d677e\nT2UuQdARicFn5mcAvBHFvvwgHpDgxPLlwN136xVcuuLk1Qwfrgyuae/S7pR4ITJLISnijuF/hYie\nr4R8HCOaRDSfiLqIqKu7u9v3D4gHJNhZvBhYsAA4edJ5/bFjwJ13qtcmHriF3XP307u0nBJdSUQR\nGQhJEKfB/xGAcwFMAXAQwH1OGzHzQ8zcxsxtw4cPj7E5Qt6xBmHdjL3Fk096Jz5zoq7u1HEhP73L\n9nbveL6IDISkiE2WycyvWa+JaBmA38X1W4IAmA/C9vWpB8OsWafrnHt6gBdeOP2h0dcHrF59qveu\nk8vZ9f1DhgADB6r9V+NXZOA0b0B6toIJsRl8Impi5oOVt58GsNVre0EIi59B2J4eFdoZNqzfcD74\noFo3bdrpBr+3t/8hYWJcq/X9pZJa6urUvvv61PuGBn8iA6d5A6tWqX3I2JWgI5Ii5kS0AsAnAQwD\n8BqAxZX3UwAwgH0A/qvtAeCIFDEXgOAe7Jw5So3jNDvViVKpP/egZXzHjgW2b3d+cJRKKnSzYoW+\n/W4FzuvrgUsvBd54w/+sSa/9SsH0YmNaxDwSD5+Z5zh8/JMo9i0UizAerNfsVCfsDwZLxeNm7K1t\nTAZXvUJLvb3A0KHAH/9o1kbT/Uo6EcGEzM60FfJH2IR4fgZh3TDNjumFSebMIMn9JJ2IEBYx+EJq\niCIlcLV6pqVFhWrs8kmvrJbWNk6YDq7q9P1HjgRL463bb3e3ZIkVvBGDL6SGqDxY+9yMrVuBp54C\nLr9caeHPOkulRnAz+qWSKjwRZga3ib4/SBpv3X5fflnqQQjeFMLgS1WsdGOdH68UA6WSkjYGOY+r\nVwP/+q/Am28q7/rQIfcHS0MDsGxZuBnccc2wte/X6YHFLPUgBA0mdRCTWuKoaSt1QdNN9flxWwYM\nYG5o8H8evWrNAvFeFxs2MA8ebF6X1M9+m5vNau0KxQBFqWnrhVTFSjdO56eaUklJGYmUR+z3PHqN\nC9hTGceRh8nvDFvTnmh7u5o/4IYM4GaTRCIRJk+FpJaoPfzOTnfPUbyg2uN1fgDlHXd2Ms+YEfw8\nXnhh9B62H7x6GI2Naj2z/55oZ6faxu2YzJkT7/8lREvYSATEwxcZW9rRzYwdMUINVHZ1BT+PcRTK\n8eOJmSRaC9ITlXoQ+SHRSITJUyGpJQ4PX7yg9OJ1fqz4el2dt4euO4+mHrYpQT2xDRvU/9vaqv7a\nfzdoT9Rqi3UMZXwqm0QRiYChh19zI29fojb4y5Yxl8vR3exCtOgGVE0Wk/MYlWGM+uFhESbs5PUg\nEbKB7vwPHqw/r4U3+NZN7uQ1iheUHhYtcn8omxh76zxahu/CC50NXxSGMa4xIemJFhtdT7f6Wnei\n0AbfyxMrl5XnL6SH8eODGfx589T3naSdRMwtLdF6vKaeuO7hY2FtN3689ESLjGlP1+taKLTBF3VO\ntjDxcNxugGXLvG+WhoboenMmnrhpjN9t/oHb90wfIkI2cYtImNquQhv8WkvxBH8EjeWXSmoCkm7S\nlt0zCmM4dTF8r4dPdRu8eqAf+cipbZPJg8UgzEQ9U4OfS1lmHFI8IT7cpIsDNMm7+/qAV15Rt4IX\nVgqDxYtVrplf/hLYvFnlzveTe0YnsVy3ziz5m9dkMGZgypT++sx+JXuSRiS7JFIK0+SpkNSSRAxf\nYqLpxWlgdcYM/56/0zJ+fHTXhNsAsGnP0k8P1E94UnoC2Seo7UKRQzrMolHOC1FIN3Whn6jGdUzV\nNn5UOX4GisXJyQdBbFeiBh/ATwEcArDV9tkQAGsA7K78HazbT9Q6fNEo54PqG8Dv0tioVwJFcemZ\nGl2T7axr1yumW/0QEaFCfvBru5I2+JcCuLDK4P8TgNsrr28HsFS3nziyZQr5wH4DDBpkZujtnlFS\nWndT78xrO9MMovaHSBIPNCG9JB7SATCmyuDvBNBUed0EYKduH2LwBQs3NY1uohaRCt9Ue0ZJhjxM\nvTO7Dr+5Wf396Ef1PZnqh4jumMjkrfyTBoP/pu012d9XfW8+gC4AXc3NzfEeFSETuA0+zpunj+d7\nGe80juuYevPWYmUQNXmQWUt9vRr8Fg1/fkmVwa+8P6Lbh3j4gk6j7mUYy+VwScySJsiA9ODBpxpu\nXYppKwGdKHfyjanB1yidQ/EaETUx80EiaoIa1BUET7w06idPen/3wx/WFzBpbzerS5sEXv+rG0eO\nqGXLFmDVKlWjV/lT7vT29r+2a/hnzUrPsRCSIc6JV6sAzK28ngvg8Rh/ywiZlJJ+dDny3SiVgAsv\njL49cRL0fwX6DfehQ+4F2QH3/futpyvkg0gMPhGtALABwHgiOkBE/wjgewBmENFuAB2V9zUj7CxL\nIRp0D12vWdJEQLnsvC6LRT+8/ldT+vq8j5fX96oLx4hDVABM4j5JLVHG8O0qj44OlUQrCYWG4I7J\nTFCdmsYauE3TwGtQophUBjCPGuV8TFpazKWoMks326DIM239KB+IZFJKEviRRerUNPaB146ObCtQ\nvCaV+THcToPRUU4EE9JNYQ1+EK9p1KjQPyto8DsT1ERNkxev1OkBFsRwO2EiRfU6N+IQJUvQbK6F\nNfg6mZrTUi6LFxM3UaesLppXGmYOge7hqTs34hAlQxgHprAGX3fxihcTP05eStSpDYqYOyauOQQ6\nJ0kcovgJ68CYGvw4dfg1Ydw4pVHu6zP/DvPpigUhGIsXK4338ePquFp68c5OpaQ5duz07wRR2HhJ\nGp0UKHkgrjkECxcCv/qV+zyHvj4l4RTNfnx4zcmwJLRRHP/cFUBZuFAZED9IUZRo8CrWsXKlMvpu\nxUP8XsxS5CY62tuBESPc14tDFD9JOTC5M/huVYnq64G6OufvECltvmiPw6HzUo4dU1Whrr9eab2v\nv169182OdcLrwZ5FTX4QotTNX3aZu25fHqDxk5gDYxL3SWqJQ4dvj3fq8qpnVeWRFpKuJZzGZGhJ\nEbVCqWiD4GkjqRh+zY28fUkieZo9Ja1bSlm5wIORVM55O2lKhpYUURhnp4H1Ij9A00CY4y8GX0MR\nVR5xI15iMoS9dr16B0V8gKaJoMff1ODnTqVjShFVHnFjjZ/ce6+K2Vt5Xhoagg3MCs6EuXbtA+v2\n79gzaK5YEW17BXPizuaau0FbU0TlEQ9LlpgNzEqiruCEuXZN5H9CjjHpBiS1JBnSkfBD7chLSoRa\nEebaDTOwHnTavxA/MAzpFNbDd5NvBtWFC2Z4afXvvVc8fRPCXLtBeweSXjwnmDwVklpqUeJQBqmS\nRQbLo8Pt2vXyxIP0DqQ3nH4gKh0hLdgN0ODBwUMKgh6n1OBEKsWyZZj9yv/kIZ1+TA1+YVU6QjJU\n59bxQgbLw+GkwAHUcd+2DZg2DbjtNjWAPmuWGqDdvRs47zw1c9ktFCSKtvwQu8Enon0A3gZwEsAJ\nZm6L+zeFdOBmgNwoSkqEuNAVRe/pObV4uek4lVdCQnlIZ4ukBm2nMfMUMfbFQmeALGSwPBpMiqIH\nkV5K3qJ4qIU0ubAqHSF+dAaooSF8EjWhH5Oi6EFCMKJoi55aqZ6SMPgMYC0RbSKi+dUriWg+EXUR\nUVd3d3cCzRGSYtw49wyMANDbCzz4oJrZKUYjPCapwYOGYEwn1Al6ailNJtb1AcP+ANHZzPwKEY0A\nsAbAV5j5Gadt29rauKurK9b2CMmxcSNwySXuhTWIgNmzZSp/lCxeDCxdCrz7rvP6+nrg0kuBw4eV\n4fcarBXiYc4c5dm7md6ODmDNGn/7JKJNJiHz2D18Zn6l8vcQgN8A+FjcvylER5g4oxTWSJ4lS4Cn\nnwZaWk4N75RKwIABypNcu1YmT9USXajzySfjOyexGnwiOpOIBlmvAcwEsDXO3xSiI4o4o66wxpAh\nklMnatrbga1bgfXrVZWx1lbg8suVwe/tPT2MsHQpMHPmqedAch3Fh26spa8vxtCOiVg/6ALgXADP\nVZZtAO7w2r5WE68kR8jpRDW70ms/AwYwNzRITp0k0BUqt0+kGjCAua5OzktceN0TQSe0QWbamiGJ\nvJyJcnal08zO+nplVGS6fjLokqaZLHJeomPRIv0D2I85NDX4hZZlSiIvd6KcXemk8Lj0UuDECeft\nJU1v9JhINnXIeYmOJUvU4KwbcU1oK7TBl9zg7kRdL6C9XalxurrU38OHZbp+kphINnXIeYmW735X\nzWVwIq4JbYU2+JIjxJ24Z1dKAZrkGTs2vJc/ZEj/axnYDUdNJrSZxH2SWpKO4dei6HaWCFJU2XQA\nXFLuJodTBs2gS0OD2p+MfZmjuyeiSNEOGbTVI0ZHj5+L0a8RCPJAEfxhogipXspltbitr69Xht9p\nHRHzjBly71gk9WAUg6/BMmSjRqmLWzyVcAR9eEoBmngxlWPal44OZbTDKnqKfg8tW+b+4Iz6wSgG\n34Pqpy6ROjFNTWJ0vPDqmkqRjHQSRI7Z2ioyzrAsWuTdS4r6wWhq8As3aOskxWRW+V6OHi1ebhHT\ngTfdrFsZAE8nfuWY1oC5yDiDY9kYtxxSdpKWgBfO4IsUsx/T1Akm8xWGDnX/HVHd1A6/ckxLgeX1\nvfp6s30W9UFvWgfCIkm7UziDL56owsuI33WXyq1ieR26h+TNNwPPOOY/VUiRjNrhJv0bMACoq3OX\nA3pJBr/xDVUq0U1DblHUB71JIRo7SdqdwtW0lXJtCi8jzqzSs65fr2563UNyxw7n4wkoYy9FMmqL\nWw1bwLuure57o0cDb74JHDrkfH0U9UHvZWOcSNTumAT6k1qSGLQVKabCdFCusVGpCdzmK5goPoT8\n4CQztJKtibxW4VcKG4XdgQzaOiPl2hReMXc7PT3qshw4MNjvbNokMzHzglsY8MQJdQ9Nny7VsABn\nG+MGkZoBnRgmT4Wkllro8O3676KkSV60SE2eMfVARo1SXlwYiV7Rvb48INJbf2zYoHq4ut5xFPcG\nRIfvj6JMFffb3bTmKITVZBcxbJY3dGHAGt6+qcXPxLcw94apwS9cSMeJIqVJ9isZYzbTE5tSNOlr\nnpCEd/7xo9hJ4t6I3eAT0ZVEtJOI9hDR7XH/XhCKpM33KxkLQrnsvq5I0te8EXcG1TziZwJbEvdG\n3DVtywB+AGAWgIkA5hDRxDh/MwhF0uZ7XYBEwMiR4WZYlkrA2WeLJ5hHnAYjidQD/qyzlGOUp95w\nFPiZ+JbEvRG3h/8xAHuY+SVmfg/ASgDXxvybvilSV9XrAjzjDOCCC8L1ABoagDvvFE8wr9irl40a\npe6Pvj7g4MFgRe7zjh/FThL3RtwG/2wAL9veH6h8liqK1FXVyVK9KlFZ6GZq3nyzSF/zTHu7umfe\nekuN71SPey1dqmZqixxXUV3is6VF2ZWa3BsmI7tBFwD/GcBy2/sbATxYtc18AF0Aupqbm4MNUUdA\n0XKzu6Ul9ioKAzAPGXK6jNUttbGkPs4vpuqTvN9HQYn63oChSoc4xhE8Ivo4gO8w8xWV99+sPGT+\nh9P2bW1t3NXVFVt7dGzc6D3VvAhs3Ki65ceOnb6usVF5KkU7JsLptLaqhHumyLUTL0S0iZnbdNvF\nHdL5M4DziGgsEQ0E0AlgVcy/GZjqQttFvDhlJrJggt/0yXlTu1WTlfq+sRp8Zj4B4MsAngCwA8Bj\nzLwtzt+MmqycyCipjjl6TZUv4vER/KddzpvazY5pmvFUYBL3SWqp5UxbJ4oy+zYocnyKTfW4l27G\n9pw5tW5x9KQlGSNkpm04ijT7NghyfAR7T3DwYO9tS6V8qd0ssjZpUwy+C1k7kV7EEXbJ0/ERgmON\ne+kyPg4fns/xn6xN2hSD70LWTqQbccUX83J8hGjQzeCeNi3Z9iRF1iZtisF3IWsn0ok4wy55OD5C\ndOhmcOcxnANkb9KmGHwXsnYinYgz7JKH4yNEh07OC+RTzZU1GbMYfBeydiKdiDPskofjI0SLm5wX\nyJBsMQB+ZMy1JtaZtn6p9UxbJ7I8+3bOHHVzuRVsnz0bePTR/v9x1y4VitH9j/bthwxRMdo33sje\n8RHiR2ZuJ4PpTNuaa+/tS9p0+FnHRCPsV0sv2vviEab0p5RFTAaIDl8wiav6GdQV7X3xCKvyEjVX\nuhCD70JWUgbo2ukVX/Q7qCva+2IRxQNe1Fwpw6QbkNSSlpBOVsIWYdvptyi1FLEuFlGEY9KSeiDv\nQEI6wchK2CJu7wsAXnrp1F6DeGvFIopwjKi50oUY/CqyEraIop26jIdHjpwasxXtfbGI6gFvKlvM\nShg105h0A5Ja0hDSyUrYIqp2mmY8rFb1FKUyWJFJMhyTlTBqWoGEdIJh4tWkwROJw/vyynho9Rqy\nNMlECEdS4ZishFFzgclTIaklDR7+hg3MDQ3OXk25zDx0qPpr4omE0S+btDNq7ysrvRshWcLWX9Xd\nB6LVDw8MPfyaG3n7kgaDv2gR84AB3obPy8haF/eoUeYPhjBtjTK84lXAvFTKZwELIV5MQjXiaISn\n5gYfwHcAvAJgS2W5SvedWht8L69Zt5RKzC0t+u8H9b7dvKSw3pfp/y8SOsEvpteTl6MBMDc3y7Wn\nIy0G/1Y/36m1wffqWpoafZNt/HZRkxzQkkFZISpMQzU6R4tIrkEdpgZfBm1teOmOTXBKUua0jZ/p\n5FEPaIWZmSsIfjDV8dsHh4lO35ZZBnCjIm6D/xUiep6IfkpEmqqXtUc3ESkK/E5Q8tLbHzumMl6a\n3gSmeVGssnVdXeqvTI4RguBHSWY5GqNHu+8vTfNgskoo80ZEa4loq8NyLYAfATgXwBQABwHc57KP\n+UTURURd3d3dYZoTGt1EpCgwnaBkeeKPP+7d69i/3yyZlUjfhKTxO1GvvR0YNsx9f2lMtpYGibYv\nTOI+YRcAYwBs1W1X6xg+s3MMu1wOHtcPEoesjtn7UQm5IdI3oRb4HRNKo1LMTTCRpsliSMGgbZPt\n9X8DsFL3nTQYfObTlS8zZgQbzCVSD4umJncFTfXFtGxZMKWQzmiL9E2oFU5KMi/VWZqUYm5Gfd68\ndLUzDQb/EQB/AfA8gFX2B4DbkhaDX43OWDotH/iAXibpdDGF6U14Hb40ek5CMdF5xmlRink9fLwU\nebXoMZsa/NiGKJn5RmaexMwXMPOnmPlgXL8VN0EGc997T33PbcDTLaZ+8mSwNuoGgyXxmZAGvMaS\n7roLmDkTmDUrHUoxL8GElyIvjWMNFiLLNCDIYG5Pj/dgqNfF5AYRUC47r9MZbUlTKySN04Cm13XP\nDKxZo0QIq1fXXikWVKad6lThJt2ApJa0hnSYzbNKmnbtgoSJ7LHDoN3dKGfmCoIbbmGbpibzaz2p\na9NtPEE3A9htKZeZx49P9v5CrWP4QZY0G3zm/gtj8GDzk+/2L+kupnLZ3ajrjHacSdsEQYdX7Nue\nXyqosxQl1Q8mu9Cio8M9kaKpw5fU2IMY/BgxzbnjNRiqUyMsWxbME0+TVEwoJl4SYMughnGWosIk\npcOAAcx1dcE8fWupr4/f6RKDHzMmWnldtzRqNULaJG1CMdGFK0eN0jtMcSnH7L3f5maz3kZ9vZJm\nt7Yqr3/gQP9Gv6Ul+v/FjqnBl0HbgFhTwWfPBkaNUoOpVh4Q08HQqPPWZKU8o5BvdCkVpk1T13lH\nh/t2cSjHqlOL7N+vzLGO3l5g6FA1gHzxxWY5s6rZsSMls3BNngpJLVny8KtJw2CoTK4S0oCfnqau\nlxvVeFSY1OfWvRNUl28tcY5JQEI6tSPoRRr24vYaCCZSXVgZyBWSwE+40s1ZsvZhv47LZaVU80uY\n1OdWeEk3NmHy0IgLMfg1IuigaRSDrSaDUFGMFQiCCWF6vbpr2cvoOzlOQWTQ9ofMsmX6fXgpeuKe\nzS4GvwYEHTSNcrC12rMKM6gsCLWis1NvhJ2uXTfHqaUluNLGSnzotY9SSQ3suhn9uO81U4Mvg7YR\nEnTQVJfz/s479b9tzWr83e/UwNL06WogePRo56ISujYJQi3Ztct7/cmTp1+7Xmkb9uwBBg4M1hZm\n/T4aGoDvfhe47bZ0z2YXgx8hphV+/HwPAJ580jvffbX64MkngfXrgauvVvnFg7RJEGqJSWqC6mvX\ny3Hq7QX+7u+U8Q2KfR9uBj3tFePE4EeInwo/pt8DlGF2y8ujK2wydGiwNglCHJgWDFm40D1vFKB6\nrdXXrs7hamjQy0HdesPV+/Ay6KmuGGcS90lqkRi+++I21VxX2KSjQyZjCenArzBh3jx/409+UoA7\nqYjq61Va81oNvIYBMmhbG4LOnl20SC/tcjo8Jtr7tOQXF4qLl1NDpAY8nZyPefNOTcXgVTnOr8Nl\nVxG1tOjz5qTZQRKDX0OCytFmzHC/2Ny8C1PtfUdH//Rw0eELSWOig/er03ciiHNjImdOu4NkavBJ\nbZsO2trauKurq9bNqBkbN6rB12PHTl/X2KgGpdatU7HKceP6B2rXrXOPXRKpdaWSij/eemt6BpCE\n4tDaqgQFOhob1fUcJu5t5d3fvRs47zw1HlC9P2ubXbuA118HXn7Z/R5qblb3Wapi8VUQ0SZmbtNu\naPJUcFsAXAdgG4A+AG1V674JYA+AnQCuMNlfXjz8MLh5KJde6r+weda6pUJ+Mc0tn0RaZJPEh7pQ\natpAQjr8rQA+A+CZqqfNRACdAFoAXAngh0TkMeYuWDjJuh54QI3425U4QRHtvVALTKvGmUiFTZU+\nTts5qdoLH9rwAAAQfUlEQVS8yJ2SzeSpoFsAPA2bhw/l3X/T9v4JAB/X7Uc8fGfC5AHJqsci5A+n\n3Dh+lTCmSh+vGbd+7qWs9IhR45m2ZwN42fb+QOWzQmLqkbgRtLamE7nzWITMYPVeg6ZF1s05se4r\nr+22bzf37NM0QzYqtAafiNYS0VaH5dooGkBE84moi4i6uru7o9hlqqieBfvYY+q918zZanQTs/wQ\nR55xQTClvV0VKv/2t/2nIDBNXaIrlO4GkRqgTeMM2cgw6QboFkhIx5GokqKFzeWdFWmZUCz8ypdN\n6z3otnMbPNbdk2muFY0kdfgOBr8FwHMA6gGMBfASgLJuP3kz+LrYe3Oz/1q11RerPUZZXX/TXpA5\nbReoIPjFZCbthg3qvvIy9uef71+rn/Za0aYGP5QOn4g+DeB/ARgO4E0AW5j5isq6OwDMA3ACwNeY\nebVuf3nT4Ztoj61urEnXsVpfPH266nba9caAXoMsCFlEN0+lsxNYudJ5vX27devUa9P7RPe7YecN\nREEiOvyolzx6+Cba46woAQSh1rjNU5k3L77Zsrp8VXHPGzABkg+/9phqj0UbLwgKnaLNLf3wsWPu\nA7WAqgsRdBA2aNrzNCIGP0ba21W4RpeD2+uiCSvpFISsYKpos9IPP/igen/LLcATT3grcIYPDx52\nCZr2PJWYdAOSWvIW0rEwGUhymmyS9oEiQYgKv4o2P+kR7EkEg4gXoixBGheQbJnpIkjq1rRfZIIQ\nFSZxcksWOX78qSmTTWXJYZymtKcYNzX4EtJJCHt4x2SySdD6uIKQRXRx8qef7g/37NypatrqsFev\nsvbtNDPXhLSXLjRlQK0bUCSWLAFmzTKTg+VpoEgQdIwbB2zZoq7taoiA7m4zIw8AgwcD556rvuOW\n9thymvzE9dvbay+/DIsY/IQxvWi8boDMDRQJgoaFC4FVq5y17qWS833gRKkEXHkl8OijyhPfv995\nu6I6TRLSSRC74mbGDGDmTHf1jZekU/LhCHnDK+Q5YoS3AseO/d7IlbomKkwC/UkteR609VIV6FK8\npnWgSBCixim/jskERqd7o0jCB0iJw/TgNTXbjtM0bZNybYKQZ7zun3JZ3RdTpjjfG4sXqwHanh4V\nxslrqU/T1Api8BNgzhylLtAd6lJJjf6vWJFMuwQhK4Qx3GGcJnvt23Hj0utwicFPEaYFnAGlMPjD\nH9J5UQlCLUm6t2s9ZI4fV85amnsHYvBTxJw5apq4qdLATwZNQRCiJwsZMu2YGnxR6SSAaRI1iyAT\nQwRBcMdvTqq8TnwUg58ATpIzHVm+qAQhTQQpM5rXiY9i8BOiemp2RwcwaJD79lm+qAQhLZgWPq8m\nrxp+MfgJYqV17epShZyvvjqfF5UgJIUuVBM0NJPXiY+hDD4RXUdE24ioj4jabJ+PIaLjRLSlsvw4\nfFPzR14vKkFIgupQzcqVwMUXA+ef32/4g4Zm/CY7zAphPfytAD4D4BmHdS8y85TKsiDk7+SSvF5U\nghA3TqEaQL3etg2YNk09EMKEZvKSIdNOqORpzLwDAMieh1TwhZ8MmoIgKLxCNYAK19x7r9rOLSmb\nSS86Dxky7cSZLXMsEW0BcBTAt5n5/8b4W5kmbxeVIMSNV6jGoqdHeeS33uo+S7do953W4BPRWgCj\nHFbdwcyPu3ztIIBmZj5MRK0AfktELcz8lsP+5wOYDwDNzc3mLRcEobB4pQ+3sGL0K1ZIL9pCa/CZ\nucPvTpn5XQDvVl5vIqIXAYwDcNo0WmZ+CMBDgJpp6/e3BEEoHl758y3sMXrpRStikWUS0XAiKlde\nnwvgPAAvxfFbgiAUD0vwUF/vvo0o3U4nrCzz00R0AMDHAfyeiJ6orLoUwPOVGP7/AbCAmd8I11RB\nEIR+lixRtW5bWk5V4ojSzR1JniYIQuYpet0I0+RpUtNWEITMIzF6MyS1giAIQkEQgy8IglAQUh/S\n6e3txYEDB9DT01PrphSGhoYGnHPOOairq6t1UwRBiJDUG/wDBw5g0KBBGDNmjKRwSABmxuHDh3Hg\nwAGMHTu21s0RBCFCUh/S6enpwdChQ8XYJwQRYejQodKjEoQcknqDD0hytqSR4y0I+SQTBt8PfmtX\nmlAulzFlyhS0tLRg8uTJuO+++9CnqUi+b98+PProo+F/XMPNN9+M7du3e27z29/+VruNIAj5J1cG\nP0jtShPOOOMMbNmyBdu2bcOaNWuwevVqLNEkxU7K4C9fvhwTJ0703EYMviAIANQgXVqW1tZWrmb7\n9u2nfebEhg3MjY3MKmnqqUtjo1oflDPPPPOU9y+++CIPGTKE+/r6eO/evXzJJZfw1KlTeerUqbx+\n/XpmZr7ooov4Ax/4AE+ePJm///3vu25nZ+/evTx+/Hj+3Oc+xx/5yEf4s5/9LL/zzjvMzLx27Vqe\nMmUKn3/++XzTTTdxT08PMzNfdtll/Oc///n9dn7rW9/iCy64gC+66CJ+9dVXef369Tx48GAeM2YM\nT548mffs2cMPPPAAT5gwgSdNmsSzZ892/J9Nj7sgCLUHQBcb2NiaG3n7Esbgd3YyEzkb/FJJrQ9K\ntcFnZj7rrLP41Vdf5XfeeYePHz/OzMy7du1i63946qmn+Oqrr35/e7ft7Ozdu5cB8LPPPsvMzDfd\ndBPfc889fPz4cT7nnHN4586dzMx844038v3338/Mpxp8ALxq1SpmZv7617/Od911FzMzz507l3/1\nq1+9/ztNTU3vPzCOHDni+D+LwReE7GBq8HMT0glauzIsvb29+MIXvoBJkybhuuuucw2dmG43evRo\nfOITnwAA3HDDDXj22Wexc+dOjB07FuMquV7nzp2LZ545varkwIEDcc011wAAWltbsW/fPsffuOCC\nC/D5z38eP//5zzFgQOqVuYIgRERuDH6Y2pV+eemll1AulzFixAjcf//9GDlyJJ577jl0dXXhvffe\nc/yO6XbVChk/ipm6urr3ty+Xyzhx4oTjdr///e9xyy23YPPmzfjoRz/qup0gCPkiNwZ/4UKV/9qJ\nKPNid3d3Y8GCBfjyl78MIsLRo0fR1NSEUqmERx55BCdPngQADBo0CG+//fb733Pbrpr9+/djw4YN\nAIBHH30Ul1xyCcaPH499+/Zhz549AIBHHnkEl112mXGb7W3p6+vDyy+/jGnTpmHp0qU4evQo/vrX\nvwY6FoIgZIvcGHyrIEJjY7+nH1Ve7OPHj78vy+zo6MDMmTOxuCL9+dKXvoSHH34YkydPxgsvvIAz\nzzwTgAqblMtlTJ48Gffff7/rdtWMHz8eP/jBDzBhwgQcOXIEX/ziF9HQ0ICf/exnuO666zBp0iSU\nSiUsWLDAuP2dnZ245557MHXqVOzevRs33HADJk2ahKlTp+KrX/0qPvjBDwY/OIIgZIbU58PfsWMH\nJkyYYLyPLOfF3rdvH6655hps3bq11k3xfdwFQagdhc2HL3mxBUEQnAlb4vAeInqBiJ4not8Q0Qdt\n675JRHuIaCcRXRG+qflnzJgxqfDuBUHIJ2Fj+GsAnM/MFwDYBeCbAEBEEwF0AmgBcCWAH1pFzQVB\nEITaEMrgM/MfmdnS9G0EcE7l9bUAVjLzu8y8F8AeAB8L8Tthmin4RI63IOSTKFU68wCsrrw+G8DL\ntnUHKp/5pqGhAYcPHxYjlBBcyYff4KZxFQQhs2gHbYloLYBRDqvuYObHK9vcAeAEgF/4bQARzQcw\nHwCam5tPW3/OOefgwIED6O7u9rtrISBWxStBEPKF1uAzc4fXeiL6BwDXAJjO/W74KwBG2zY7p/KZ\n0/4fAvAQoGSZ1evr6uqk8pIgCEIEhFXpXAngNgCfYuZjtlWrAHQSUT0RjQVwHoD/F+a3BEEQhHCE\n1eE/CKAewJpKDpeNzLyAmbcR0WMAtkOFem5hZudcAoIgCEIihDL4zPxhj3V3A7g7zP4FQRCE6EhV\nagUi6gbw7yF2MQzA6xE1J0qkXf6QdvlD2uWPPLbrb5l5uG6jVBn8sBBRl0k+iaSRdvlD2uUPaZc/\nityu3GTLFARBELwRgy8IglAQ8mbwH6p1A1yQdvlD2uUPaZc/CtuuXMXwBUEQBHfy5uELgiAILmTK\n4BPRdUS0jYj6iKitap02/z4RDSGiNUS0u/J3cEzt/CURbaks+4hoi8t2+4joL5Xtupy2ibhd3yGi\nV2xtu8pluysrx3EPEd2eQLtc6ypUbRf78dL976T4n5X1zxPRhXG0w+F3RxPRU0S0vXIPLHTY5pNE\ndNR2fhcl1DbP81KLY0ZE423HYQsRvUVEX6vaJpHjRUQ/JaJDRLTV9pmRLYr8XmTmzCwAJgAYD+Bp\nAG22zycCeA5q1u9YAC8CKDt8/58A3F55fTuApQm0+T4Ai1zW7QMwLMHj9x0At2q2KVeO37kABlaO\n68SY2zUTwIDK66Vu5yXu42XyvwO4CiorLAFoB/CnhM5dE4ALK68HQdWfqG7bJwH8LqnryfS81OqY\nVZ3XV6G06okfLwCXArgQwFbbZ1pbFMe9mCkPn5l3MPNOh1Wm+fevBfBw5fXDAP4+npYqSOWbuB7A\nijh/J2I+BmAPM7/EzO8BWAl13GKD3esqJI3J/34tgP/Nio0APkhETXE3jJkPMvPmyuu3AexAwJTj\nNaAmx8zGdAAvMnOYSZ2BYeZnALxR9bGJLYr8XsyUwffANP/+SGY+WHn9KoCRMbfrPwF4jZl3u6xn\nAGuJaFMlTXQSfKXSrf6pSzcysloGAbHXVagm7uNl8r/X+viAiMYAmArgTw6rL66c39VE1JJQk3Tn\npdbHrBPuTlctjhdgZosiP26pK2JOBvn3o4CZmYgCS5QM2zkH3t79Jcz8ChGNgEpA90LFGwiMV7sA\n/AjAXVA36F1Q4aZ5YX4vinaxeV2FyI9X1iCivwHwLwC+xsxvVa3eDKCZmf9aGZ/5LVSm2rhJ7Xkh\nooEAPoVK+dUqanW8TiGsLfJD6gw+a/Lvu2Caf/81Impi5oOVLuWhIG0EjOoEDADwGQCtHvt4pfL3\nEBH9BqoLF+pGMT1+RLQMwO8cVhnXMoiyXeRcV6F6H5EfrypM/vdYjo8JRFQHZex/wcy/rl5vfwAw\n8x+I6IdENIyZY80bY3BeanbMAMwCsJmZX6teUavjVcHEFkV+3PIS0jHNv78KwNzK67kAIusxONAB\n4AVmPuC0kojOJKJB1muogcutTttGRVXc9NMuv/dnAOcR0diKd9QJddzibJdbXQX7NkkcL5P/fRWA\n/1JRnrQDOGrrmsdGZTzoJwB2MPP3XbYZVdkORPQxqPv7cMztMjkvNTlmFVx72bU4XjZMbFH092Lc\nI9RRLlBG6gCAdwG8BuAJ27o7oEa0dwKYZft8OSqKHgBDAawDsBvAWgBDYmzrPwNYUPXZhwD8ofL6\nXKhR9+cAbIMKbcR9/B4B8BcAz1cunKbqdlXeXwWlAnkxoXbtgYpVbqksP67V8XL63wEssM4llNLk\nB5X1f4FNLRbzMboEKhT3vO04XVXVti9Xjs1zUIPfFyfQLsfzkpJjdiaUAT/L9lnixwvqgXMQQG/F\nfv2jmy2K+16UmbaCIAgFIS8hHUEQBEGDGHxBEISCIAZfEAShIIjBFwRBKAhi8AVBEAqCGHxBEISC\nIAZfEAShIIjBFwRBKAj/HwSavK75jtIOAAAAAElFTkSuQmCC\n", "text/plain": [ "" ] }, "metadata": {}, "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'])" ] }, { "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": "code", "execution_count": 82, "metadata": { "collapsed": true }, "outputs": [], "source": [ "degree=23\n", "\n", "DataMat=data_matrix(x_data,degree)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Least Squares\n", "Recall the system of equations we are trying to solve here:\n", "\n", "$$\\texttt{DataMat}\\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", "Why does this make sense? Think about it.\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 a $(d+1)$ dimensional numpy array called params." ] }, { "cell_type": "code", "execution_count": 83, "metadata": { "collapsed": true }, "outputs": [], "source": [ "D = DataMat\n", "\n", "# Least Squares\n", "params = np.dot(np.linalg.inv(np.dot(np.transpose(D), D)), np.dot(np.transpose(D),y_data))" ] }, { "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. \n", "\n", "Once you're done, re-run all the blocks of code for different degrees $d$! What do you observe?" ] }, { "cell_type": "code", "execution_count": 84, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "38.4588390567\n" ] }, { "data": { "text/plain": [ "" ] }, "execution_count": 84, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEKCAYAAAARnO4WAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsnXecFEX2wL81u7BLliBBcAUDIDkpaxYFFb0zBwycCfXn\n6YkB71TORczeoR5nFsyeiFk8zwSKcRddEQVUchIkKDkssDvv90f17Pb29vT0zE7amfp+Pv3Zne7q\n6uru6levXr16pUQEg8FgMGQ+gVQXwGAwGAzJwQh8g8FgyBKMwDcYDIYswQh8g8FgyBKMwDcYDIYs\nwQh8g8FgyBLSQuArpS5SSolt26KU+l4pdbVSKjfKvDpaeVyUoOImBaXUbUqpqH1m433/SqmAUupf\nSqlflVJBpdRbHmmX2t5huVLqd6XUDKXUvUqpjvEoT7qilGpqvbMZSqn1SqmNSqmvlFKnuqS9Tin1\njfV8ypRSC5VS9yulWqai7NFgfauXxDG/e5VSP1jPa4dS6melVJFSqqFL2guVUt8qpTYrpdYppT5S\nSh0RwzXft+ronY79oW/HbdvDkbZAKfWcUmq5Ve75Sqk7lVKNHOlaKqXGK6UWW+mWKKUeVkrt6Uj3\njFLqJ+vetlry7y9KqRxHuhyl1K1WPjuVUguUUtf6vfeohGkSOAv4BWhq/f8Q0BooSmWhUsRE4P1U\nFwI4ExgJ3AAUA79HSP8BcBuggD2AfsBlwNVKqeEi8mbiippSCoA/A88BY4EK4FzgTaXU1SLyiC1t\nC+ANYA6wBeiLruODlFIDRCSY1JJHx0VoufF0nPJrCjwDzAN2AocCo4H+wCmhREqpy4EngMeBm4CG\nwPXAR0qpQ0TkOz8XU0qdC/SOkOweYIpj3xZbHo2AqUA94FZgOXAQ+r0fAJxjpVNWPp3R7/cnoBtw\nOzDAKndIqWuAlneLAAGOB8YD+6O/vxCPot/BHcAMYBAwTinVWESqNWCuiEjKN+sGBNjfsf9jYFOU\neXW08roo1feVomcZ1/sHxlj5BXykXQq86LK/MfAlsB3okOTnkZek6zQCGrrsnwYs93H+FdZz7p/q\nOhShnNOBLxJ8jXusZ9HKtu8roNiRrgmwC7jHZ77NgdXohliAOx3HQ9/OiAj5HBcSyo799wLloXqA\nFvQCXOFI93/W/i4RrjMJ2GL7XYBWJG5zpHsY2AG0iPQM0sKk40Ep0FQp1RpAKVXP6jYtVUrtsv7e\nqZSqFy4DpdQNVtfH2YVSVjfrZet3qDt3hVLqdsuEsVEp9Y5SqoPj3IjlsOX3f0qpe5RSq5U2Vb2o\nlGqolNpfKfWB1X1bqJS60HGNGiYdpU1cxTaTQYlS6qRYH65S6gQrvx1KqU1KqbeUUl1sx5eitXWA\nChWjqUhEtqK13wZowWYvw1FKqWnWs9lmPZMejjQ51vP9VSm1XSn1sVKqq1We22zpbrP29Qg9W+AV\n2/HTrWe23Xp+ryqlClyey+VWl7pMKfWbUuoppVSLCPe4TUS2uxwqBfby8ZhCPadyH2lroJTqpJR6\nwapnO626Pd6R5gLHfb2glGrnSHOeUuo7q15uVkrNVkpdYR2bDhwFHKaqTB3TYylvBNyeRX1goyPd\ndmA3/k3T9wFzRGRS7YpHfeuvszwbrbIoH+kgcrl/p/ozONg65z1HuveBfGBohPzSXsN/jeot5kvW\n79vRrext6Bf+kksrfZH1uwW69furI+/jrXRHO85bal1nKHAh8Bsw3XFuNOVYhu7mHw9cZ6V7HpgN\nXAMMQXfvg0B32/m36ddT7br3owXmYCu/h61rnBDu/j2e+QlobeEj4GTgPGAhsA5ob6Xpi+5uC1Bo\nbXt65LkUFw3fdnwl8LHt90nWc3wb3X0/Ba3JbQD2tqW703o+91nP6yZgvlWu25zPDN0tvgU4xvZ+\nQ1rV08CJ6G73T8ASoIlDS9ttPevjgIutcs8AcmKo28VoIeN2LBdtmigE5gJTY/x+OlnvbRlwObqb\nfyHwH1uay637f9m6/xHAWus5NrbSHG49539Zdew4q47+zTreDZgJfG+rD90c9xNpc32G1rHG1nVX\nAU85jl9i1ZVL0abC9sBjaOHZ2cczOhwow9Kq8dbw11nX2oQ2yfR0pMu3ntun1jNpbNW1X4FHbemU\nlWYOMMBKdzDwI/A/lzIq6znsAZyBNiPdajt+mlW+vo7zhlj7I/Z0UiLgXW70IqvAXawbbo4WbBXA\nW1aaHjg+cGv/3639vRwv7SJbmmfRwkzZ9r0B/OTysqc78h9l7d8rxnJ87Ej3hrX/Atu+5lYFG2Pb\ndxsOge/IJ2A9qw+Bt13u46Jw51rpSoEFQK5DcOwGHrDtu9OrHI48l+It8Isdz3whMM2Rpim6kf2X\n7dlstX9I1v7rne+BKoE/0pG2MfrjfdqxvxPaJHCt7dlVAEWOdIdZ+Z4aZb0OCdnzXY41to6Ftvex\nNTxRXud56xntFeZ4DrAG+MSx/3Dr2tfY6vr6CNeajotJx1bvIm1LXc7t4UjzHC4NA1rol9nS/Qoc\n5OP51Ec3qHfa9rkJ/HboMYLTgSPQY09L0IK3qyNta+BzR7kn4DB9ok19bzjS/Rdo4FLOP9jSBIG7\nHce7WceudOwvsvY/EfFZxFLB4r1RJfDtW4VVkVtYaf6Mey8gVNH+4vh9kS3Nwda+wbYXuxu4ziWf\ncD2BwhjLcaMj3d3W/j0d+1cBE22/b6Omht/fqixrrAoRelY/u5TjIvu5LpUw6Kzwtg/6W9vveAr8\nEuBH6/8DrHJeQk0t8B1gppXuSCvdIEdeBYQX+AWOtCEN6FiXa/0AvGGlu8xKt59Lus3YGkIfz+Jo\ntHB6LszxAFrrOxz4C9q2PB1bAxzFtVYDL3scDwmKGrZp6529bv1/lJXuRbTw2SNM/XAT+PWt+4m0\n9XQ5N986dhRws/Ws/+NIcwq6pz7eeo9/QPdO1wE9IjyfvwOLsQlZXAR+mHP3tsrzgqO8n6C1/Aus\nOjrKSveY4/yX0N/2FVa6K6z39S41G4dm1nM4Fi0ndgF3OdJ8hO6ZHY/uCZwGrLfu57FI95NuXjqn\nob10tgDLRKTMdixkQ/3Vcc5qx/EaiMjXSqlv0d36qejubDlak3Cy3vF7p/U3P8ZybHD83uWxP58w\nKKX2Rg8A/ogWEMvR93AHcGC488LQHN19dN4D6PvYJ8r8/LI38LP1f2vr71PW5mS59TdkY17rOL7G\n4zrO+wpda2qY9Bsc6RaGSefLbVIpdRDaFPAxuq7VQLQnTqn18wul1Gy0EDkTbXaJhpbo7yYc4eos\n6PfdwirTp0qps9D1600ApdSnwPUi8oNXAURkl1Jqlo+yisu5ZVQ9i0+VUr8CzyilHhKREqWUAp4E\nXhORSo8VpdSH6Pp0B1p21MAaoxmNfg95Sqk82+E8pd0tt4hIRZj7WqGU+gKtNIa4FN2gHyAiobry\nmVJqE/CkUupxEfneGl87F61oTrOlW4zumf8Rbc4MXWuT7TlMU0rtAm5VSj0qIiut/RcB/6HKg28z\n8Fd0z8Tt/VYj3QT+HNsDdBISxG3RNlpsv+3Hw/Eo8IRSqj365b8qIpHOSUQ5YuUEtAZwtohUftzK\nxV/ZBxvQH15bl2NtScA9KKX6oAcvJ1q7QgNzN+MuiEMNY6gSt0Z3y0O08bicU6iErnWRI48QWxzp\njqNmg2w/HhalVE+0a+os4AwR2R3pHIvQh76/z/R2fkPbtMNhr7NO2gLfhn6IyGvAa0qpxmihdh/w\nvlKqg3i4iyo9z2KJj7IuQ/dCvbA/ixL0u25t2x8q6y6l1Pd4Kzz7ohWpF12OjbK2vuj35ZeewEYX\nWfW19fdA9DhHT+t3qUe6twlPKbon2Ak9joQl+I9WSu2FbqgXAb2s9F9EKni6CXwvPrP+DgPusu0/\n3/o7PcL5k4Bx6C5WAbpFTEU5YiUk2CsFiFKqM9q+7KXd1UBEtlk9nrOUUreFtBul1D5oP+iH4lPk\nynI2Bh5Be1U8Ye2ehzYndBeRez1Onw1sQ8/L+MS2/6woivAVWqjvLyJuvboQH6FNXQUi8lEU+QOg\nlDrAymMx8AcR2RHF6UdZfxd5pnLnQ+B0pVQ7EXHT8uahe0TDsPWmlFKHontz9ztPEO1Z9V+l1L5o\nM0pLtPlkJ9od0skqtC96JHZGTlLjWWywzhtgT6SUqg/0QT/vcMxCD2I7+QTdCDxF+B5dqIdwOGCf\ncLga2EMptb9D6A+0/q60pcMq9zSPdOEImdhq3J+IrAJWWb2fa9E9nekR8qs7Al9E5iilJgG3KT37\n9ivgEPTEh0kiMjvC+TuUUs+iPWVmi8hXqShHLZiKNuE8r5S6H23qGIs2fcTiXnsr2o74X6XUo+hB\nxLHowc0aAiAKWimlCtEmo2ZUTbzaEzjXqqiIiCilrgLetj7cV9Caaht0o7NcRB4QkQ1KqX8Btyil\ntqCfQz90txq0gPZERDYrpW4EHlHaPfc96z7boz+q6SLykogsUkrdBzxsuad+irbD740eB5goIp+4\nXUNp1+GP0LbsMUA3/S1W8p2I7FRKNUN3x/+DHjQXtLngerRW+IYtz6PRguliEXnW4xbHoD1vvlJK\n3Y0WYO3R3lsXiEiFUqoI3cN9ES3o2qMVlgVYk6iUUrejn/8naAHeAe2lM0tE1lnX+hH4s1LqHLRA\n3iIi80RkFzU1WU+UUr3QStiraKGWh7ZzjwTeE5FiAOu5TUBP3tuANegJXI3uLVxry/NP1v0cKyKf\nishGXASh9W6Wich027770fWpBN0r6oLugQaprtw9i35f/1NK3YX+Bgegv6lv0XNOQL/Lu4AXlFJ3\noIVyV/T7WkGV2ewktDfYO1ZeTdBegpejB2JX2cp4JbpOLkH3zi5EN0jHevXAKolk5E/GRhi3TJd0\n9dGDiMvQmu4y63c9W5qOhBm0RAtmAa5yORY6b4Rj/9HY3DdjKIczv9us/bmO/UuxDXjiPmh7tlVp\nytCmiWHoyrfUz/273PMJaM+ZHWgB+DaOySBEP2hrH3TfAHyDdnXcJ8w5h6A/4A3WfS1F27APsaXJ\nQX84q62yTkc3CoLNIyfcs7UdPxEtzDajexshYdfNkW44+qPfhvZ++QntAht20pitnoTbOlrp8qxr\nzrfy34QW9KNxeOmg3VYFm9utx/X3Q/dif7Oe4yIcg8zoAcbv0dry78ALQDvH9T5Am9F2ooXSU9i8\nf9BC5n/oHpPg8GqL8rtvg+5xL7He6+9WfbkKx4Q5tHJ6NVpj34LubUwHjgsjS46OcO0ag7ZoB4Jv\nrLq426pvL+EyQQo9EP6K9Yx2WO9zHNDckW5v6xkuoUpQT8ByfbbSdEU3Dius574GbZ45n5oDu1ej\ne2xl6EbpDWzu3JE2ZWWSFVit8Uh0Bd6c6vIYYkcpdSZaMzxSRD5PdXkSgaWtn4z2bMmeD9WQMOqM\nSac2KKX6ortnI4EnjbCvWyilBqK1zxlozaY/evJVCT4GquowR6F9sY2wN8SFrNDwlQ4R0AbdXR0u\nIlu8zzCkE0qp7uhB357oiVlr0fbOm0XEzZvGYDC4kBUC32AwGAxpEg/fYDAYDIknrWz4rVq1ko4d\nO0Z30rZt8PPPsP/+0KxZQsplyHxWrYI1ayBoc2wLBKBNG9jLT6xLgyHRbNgAixdDt27QoEG1Q99+\n++1vIrJnmDMrSSuB37FjR0pLo3LlhbVr9Vf5l7/ANdckpmCGjKakBI49trqwB/170yZ44w0oLExN\n2QyGSh59FK66CqZNg7bVJ00rpZb5yaLWJh2l1N5KqU+UUj8qpeYqpUZa+1sovQTZAutv89pey5U9\n94SGDWGJntVdUgLnngv9++u/JSUJuaqhjmOvJ+ecAzvCzIktK4Px492PGQxJZc0aUApatYo5i3ho\n+OXADSIyUynVBPhWKfURegLENBG5Vyl1E9qN7m9xuF51lIJOnWDJEsaMgXHj9McrArNmwZQpMGoU\njB0b9ysb6ijOeuJFMAgLFiSnXAaDJ2vXQsuWkBu72K61hi8iv4rITOv/LehZie3R4UxDcUueA2os\n5hw3OnVifeli7roLtm+v+oiDQf173Dij6Rs0JSW6PtjriReBAHTuXHWu6T0aUsaaNdC6deR0HsTV\nS8eKmNcXPUGmjVQFclqNd3TDmAh9gBM/7kS9lUuoqHD/gk233BBi/Pjw5hs38vP10NCYMdrOP3ky\nzJwJr7yif48Zk7iyGgzVWLNGj1fWgrgJfCsi4uvo1YOqzWS1Zgq6SmOl1w8tVUqVrlu3zi2JK6EP\n8OWXYc72fWnCVlrxm2ta0y03hJg/379m37ChNgdCzV6B6T0akk66CHylF+9+Hb1KTSja35rQAsnW\nX+cCFgCIyJMiMkBEBuy5Z0SvIqB6txxgoRVCfP8wUU7t3XJDdtO5s64PbigFBQXaZHP22doZYuxY\n716B6T0akkY6CHwrHvNT6LVKH7AdmoIO3Yn11yvQf1Q4P8D5aGnemfmu6UPdcoNh5EhdH9xo0ECb\nbEpLYdKkKldMr15BMAgffGDs+oYEs307bN2aeoGPXoBjOHCMUmqWtZ2IDok7RCm1AL0SvdciF1Hh\n/ACX0IlyclwFfk6O7pYbP2oD6HowapQ214Q0fbv5xq2eePUKQM+HMXZ9Q0JZY63oWUuBX2u3TBH5\nAr3YhRvH1jZ/Nzp31i6XoYky5dRjMftyANUN9Tk58PjjMMJ1VVFDtjJ2LAwdqnuKCxbAAQdozT+c\nUjBypHbvDZkQw2G36w8dapQMQxxZa1nE00DDTzpu3fIFHFCp4SulNbbRo42wN7hTWKjNNk7zTbi0\nzl6BF8aub4g7cdLw66TAd/sAF9CZA1jAXm2DnHNO1YCbwRAPxo7Vderss7W9vrnHvPFgEL77Lnll\nM2QB2SzwoeYH2KR/ZxqxnZXfrIqosRkMsWDvFRx/vLe2v3ChseUb4khI4KfTxKtkY/8AL73P8rs0\nDveGJODl7QNQUWF89A1xZM0a2GMPyMurVTZ1WuBXI+RoP9/dNdNgiCchs2JOTvg0xpZvqC2haAIf\nvbiGVRWta61AZI7Ab99eq1xG4BuSxNixehmGcJgZ3obaYA/nUX/jGhZuaVNrt9/MEfiBgPavM1+Y\nIYn07Rvelm9meBtixRnkrw1rWEMbtm+H++6L3VSYOQIf9NdlNHxDAggXKdPLlm9meBtixRlNICTw\nAXbujN3dPPME/qJFUF6e6pIYMgivSJnhfPRDyzQYDLFgjyZQj100Z2OlwAf46afYtPzME/jl5Xrd\nR4MhDrjFz3dGygy5CB94YJXQF9EfZTibq9/Y+iYGf3ZiNwW2tuJO2gV+MBibQ0BmCfxu3fTfH39M\nbTkMGUNRUfiQCk4vnCVLqq+LGy6Est/Y+iYGf/YycqTuJYI250B1gQ+xDVfWeYFv14Au/seBeqcR\n+IY4MGYMTJ0a/rjdC8dvCGU/PYZo0hkyk8LCKv3VTeDH6hBQpwW+UwN6/s0mLFcFfP+yEfiG2hES\nuO5hkYV8dlBPlVeb/uEnhLLfBdNNDH7DxIl64H8vVgGwir0qj8XqEFBnBX44DWiudIM5c40GZKgV\ndoHbmXlcxwNM5mzmcwDbaMQOGrJL6vHs5HwWNujB6J+HM5wXaMZG1/xCIZSXL/duGEI9hkgNiPE+\nznwKC+Gvf4V96mmBv5q2EUN5R6LW4ZFTRTgNaC7dOVqmc+mDFRRO9pgGaTBYlJTo+jR/vu4mjxwJ\nv/y0havlGS5jAj2ZA8ASOlLKAN7mFH6nJXmBcpoEN3FA2XyO4H1O50W204CXGcbd3MIiPGZluWDv\npjtDgIdLZ8hsxo6FNd+vZNP7renVo37EUN6RqLMCP5wG9CPdaEAZ235cCuyX7GIZ6hhjxuie4o4d\nuj4t/W4D/V79B/+teJRmbKaEgVzDeN7kNH5h78rzAoHqwlgR5CC+4RKeZjgvMJwXeJz/4xbuZitN\nfJXF3k33isFv/PuzizYVq+DAvSgtrX1eddakE24Voh/RIx1HtDB2fIM3drNgQMq5lgdZIPtxQ8V9\nvMdQDmYGh1DCQ1xTTdjn5NRUNoQAXzOQK9UTHNl+MW/teRlX8Qg/0IujmO5ZDrdueiwrcxkylJUr\nYa+9Iqfzg4jUegOeRi9SPse27zZgJTDL2k6MlE///v3FL8XFIg0biuhPr2prwiYRkGVX3uM7L0Nm\nUlwsMmyYSL9++m9xcfXjw4aJKCXSm++klH4iIO9znPRiVo16BTptw4Yi7drVPGbfQtW46NgvZB4H\nSDkBGcmDAkHXPLt3r1k25z307+9+D4YsoE0bkcsu80wClIofWe0nUcRM4Eign4vAHxVNPtEIfBGR\noiL9AQYC+k4CAf17Y9MOIsOHR5WXIbMI1Q2lqteNoqKqNAP6lksRt8lucuRX2sjpvOYpyAsKqgRw\nqM45t0BA5Nxzdf7FxSKtG2yWNzhVBORJRkiA8hrnNGxoBLkhDLt26Uo8ZoxnsqQKfH09OiZb4IuE\n0YCOO65KzTJkHRMmiOTkuAvkSuG6erXMbn2MCMhzDJfm/O4p7O2ae7jepZvwLioSadSgQu5gtAjI\nJM6RXHbVaCSGDUvJozKkO8uX60ryxBOeydJF4C8DfrBMPs3DnHc5UAqUFhQUxOUZybXX6i+voiI+\n+RnqDEVF4YV9aPvLgK9E2raV8rwGckX9pyMKeqfmHrqOW+/S3oMIUVws0ry5yA38UwTkdU6THHa7\nNiYGQzVKSnQFeecdz2R+BX4iB20fA/YF+gC/Ave7JRKRJ0VkgIgM2HPPPeNz5e7d9UjckiXxyc9Q\nJwgNwlZUhE9zBq/xj9JB/L6rMTnfzKDNTRf7Wpy8Xr3qnjHOJTbPPjv8OsqFhXpJxAcDoxjJvzid\nN3mMKwE98mvcLA1hWaV98GnfPi7ZJcwtU0TWhP5XSk0A/puoa9WgVy/9d/Zs2M+4ZmYLXrNTAW5g\nHOO4kS84jPN2vMUr21oxdiwMHarPXbBAz2L9+eeajUYwCO+9V907prDQ21vG7t/fogXUrw//LhtJ\nK37jVu5kDW24lTujdrN0mzdgvHYylJUr9d908tIRd5NOO9v/1wEvR8ojFhu+K1u36oGOsWPjk5+h\nTtCvXziTTFDGcqsIyMucLXnskEBAZPDgml48xcUi+fn+7PNeuA0a5+aK1KsnElBBeYLLREAur/+M\nqxkomnzDmZIMGcDNN+uKE8E8TZK9dCahzTa7gV+AS4EXgNloG/4UewMQboubwBcR6dxZ5PTT45ef\nIWlEcqcMh7v3TFDu4mYRkAlcKoqKanZ5p+Ds3r1qn5sd38/gqtegbl6eyJAhIgf32y2zWx8jFfXz\nRL7+2vdz8TtYbMgQ/vQnkQ4dIiZLqsCP1xZXgX/mmSL77x+//AxJoTYarJtADAn7x7iimrAPt4UT\n9tEMrob8+yM2GuvWieyzj0j79iKrV8cvX0PmcPTRIoceGjGZX4FfZ2faRqRXL7361datqS6JwSe1\nDQnsnJ16HQ9wC/fwOFdwJY/hp7qHruuG38FVP5EzS0qAVq3grbdg/Xo96hthpTYTUC0LWbYM9tkn\nbtlltsAXgblzU10Sg0/iERI45D3z74Nf5AFu4IOmZ3JD3iMEAno1iUCgamGJcITz2PE7uBou7EeI\nDRtsC5n06QNPPgmffQZ33lmrfNetM3HyM4qKCvjlFyPwfdG7t/77/fepLYfBN/HSYAu3TuWq0oth\n0CCOX/si06bncMwx0Lw5NGsGrVuHF/qBgF54ojYxbLwWNg9RrddywQVw4YVwxx3w6acx57tihVkR\nK6NYvRp274aCgrhlmbkCf599oEkT+OEHsy5omhN6P17TJgIB7doY8T0uXKjNI127anNJXh7vvQdf\nfQUbN2rteu3a8A1Lfj5MmODfx94Nu2nJi2q9locfhv33h/PPh99/j5ivW4MlYlbEyiiWLdN/46jh\np3yg1r7FddBWROSww2RJwRHGjS2NcQ7Shttyc7W7pOd73LxZpFs3kRYtRBYtEhFvzxb7IG0i6kVo\nhq3vQeCZM0Xq1xc5+WSRYNAz34KC8HmaAdwMYdIk/ULnzImYlKwftAVWt+nFHst/YPt2iWkQ0JBY\n3AZpnQQCkJenNdqyMo/B3GAQhg+HefP0at/77gt4jwsoBXvvHZsW74fQDNtwdnfnIHDJzr483/0+\nmDKFxwufDVs/Cwv1eG84zABu3cRpiVj+uaXhx9Gkk3Kt3r7FW8OfMOBxEZAClhotKA3xcjMErR0P\nG6b91iO6I95/v97x4IPVrhF+MpaLhp0A/PrOV8bmoUI+5mjZRBPpmr8kbI/Db8ROQ93AzR358dw/\ny7b85r7Ox2j48MVmHWKhFz/UOGa0oNTjNUgLenB15EgoLfUezM2dVQo33QSnnqpPsOHl2RJrDJto\nxoT8LGRi7+kECXAxzwDwaNnF3P/PoGv+XgO4ZkWsukU4d+T25ctYuLMgvpYIP61CsrZ4a/gXnbFZ\nBOQW7jRaUBripaWG7Ov16nlr6M3UJvm18X4ie+8t8vvvNa4R79mpsU4M81rIxK2nczFPiYCMVP8K\n2xONJmKnIX0J19P9gR7yNif7skSQ9TNtRcdFX8i+Mpmz4vKxG+JLpAFVP9uLOX+SYCAg8vnnYa8T\nL8GYqNAG7manoEzhD7KdfDm9+0+eZTIrYtVtwr3/TTSRf3O1NG8e+b1mvcAPfeSvc5rM44BqWqPR\ngtIHP/Hrw21n1J+i/7n11ojxd+IhGBMV2iBcT6cNv8o6WsrCFgP0ykeGjMTt/bdknQhYS2NGlllZ\nLfDtmtho7hABacxmAS1cJkyIy2UMcaJLl+iF/R6slw0N2on07CljR++s4doZaa3YWPA7AOw3+Fso\nXZcu4Ru98+u/ov8xkV8zFree40CKRUBO4h1fvUi/Aj8jB23trnjf0ReA3ugZtyLa/c6QPvTtG3kB\nEicPch2Nd6zlzVOe5b4H69dw7RTRUTUGDYrfzFM/A8BjxujZrpMnw8yZ2kPUbfarPd28eVXx90MT\nqkIDu/vddBacdx7B2+9g9PGlZvJgBuI2UW9/FgKwiKr1PPyGF/HET6uQrC1eGr5dE2vHShGQvzC+\nhiZmSA/7lvt3AAAgAElEQVSiteWfyH9FQO5ktBQURJ60ZdeMYg29HKmcDRvqnqMfG79XPjk5Il27\nVi/b3Teul19Ue5nLgZLPdjM4m6HYJ+oVcZtUoCSPHa69SCdks0mnuk0sKKtpLU9zUaWt1XjnpB9u\nA6u5uS6Ck62yjL1lNt2lPmW+7P8h+3o8Fg/xGgD2a+OPZiwg1DgM5kMRkPu5zrOLX5sGzZB6QrLr\neS6QpRTUqBvhZFdWC3ynBvU/TpDv6B3RDmZILW4Dq0OGVBeIofj2h/NZVDb/Ll3i52ETbgDYr40/\nmslg9sbhIa4SATmaj10bB7MaVt0nJLu+olCmcozveprVAl+kuiZ2FzfLLnJljwZlpvLXMeyNd2d+\nlp3Uk2f5U1TCPhAQT9NPvGZd+539Gs0sWXvj0JCtMo8DZCkF0pSN1RoHsxpW5lBUJLKWVvIEl/tu\nuJMq8IGngbVUX9O2BfARsMD62zxSPvH2ww9pYn/dV3s6/PBMaVzzNySHoiKRhg2C8hGDZQPNpDWr\noxL4DRtG9gSKR9XzK3T9pAvVXWfwtYMpkXIC8jQX1WhEzGpYGcLGjSIg/+l9n2834mQL/COBfg6B\n/w/gJuv/m4D7IuUT92iZIRYu1Ldq/DHrLPPu1I32fXs/JE2a+BP0Tvt6MmLP+J3k5ZUuUgTR2/m7\nCMjZ9d+sFATJaNAMSeLbb/VLe/1136ck3aQDdHQI/HlYC5cD7YB5kfJImMCvqBBp2lTkyisTk78h\n7tgHHy88c6uU7dlepG9fGfP3cs+BWqW0+capGSXT5OF3kpfdD7+gQP896CDvcBMgkqd2yneqr2xp\nuKfImjURJ68ZR4U6Rigs8g8/+D4lHQT+Rtv/yv7bcd7lQClQWlBQEOMT8sGRR4oUFiYuf0PccGq4\nRWqsCMjdJ30R0X3TS3inY+wZv+sBhLZQBNFZ/5kjkpcn6w8ZKo0aeC/OnpenB7+N504d4dZbdeUs\nK/N9SloJfOv3hkh5JEzDFxG59lqRBg1EyssTdw1DrXFq4m1ZJVtoJK9wpuTkeAvGnJzaBTFLNrHE\nEmrevEpwLx71iAjIjfzDs8dTr57x3KlTnHGGSOfOUZ2SDgI/fUw6IiLPPadvd+7cxF3DUGucg49P\nMkJ2Uk/2Y0FEYdilS6pLHx2R1gOIOD7RIChT8s+UXeRKIV+FFfjJMGMZ4siBB4qcckpUp/gV+IkM\nrTAFuND6/0Lg7QReKzJ9dYiFhy/9zkxPT2PsMfJ7MJtLeJpHuIpF7O95XiAA/foloYBxJNJ6AF4E\ng7B9h+KiXRNYwd68zDD2YEONdOHyj8s0fUP82bVLL9TRrVtCso+LwFdKTQKKgS5KqV+UUpcC9wJD\nlFILgMHW75QxdnJXyshjZ8l3njFODIkl0uIh9ng1/+RGNtOUO7gV0HFmcnLc862Li354xebxywbZ\ng/MCk2nHrzzDxUCVhHdb6DyE2wJA0SzsYkgQCxdCeXnCBH7cTDrx2OJp0rF7eQwerBfAnsFBMo1B\npmubIvzMBA3ZtY/jfRGQ67i/2ru65JL0G3iNlXisBwAibduK/LXeAyIgf+OeymfSvbt/V1QzSzdN\nePVV/QK+/Taq08jmmbbhPB8e53JZzx4CwWo2TjMpJfFE4xY55tYK+V71lkV0kvqU1RA+9oHXwYPr\ntgeK03Oohp3ep+Au/iooXxYMkwqU3Hvku5WTt+I1EcyQJMaO1UJp27aoTstage9Vea/gMRGQfVhS\nQ0MyJJaoZoJOniwC8nDhC57eNJmilbo1YPb79i2Qt20T6dNHpFkzkXnzRMSfK6rXuzEKUXJZN/gc\nWd2oU9QKTNYKfK/KexAzREBO4/Vq+3NyjBaTaHwHDNu9W7vbdOvm6UKbbVqp7zkES5eKtGqlPT02\nbRKRyK6okd6NUYiSQ1GRyGzVU97hpKgVmKwV+F6VN5/tspscuZ2/Gy0mgbiF6PUd2uDZZ/XO117z\nvEY2xo7xPYfgk0+0FnP88b6WRozkHmoUosRTXCzStMEuKaO+3MeNUSswWSvwvQQLiMyme2UL6qph\nGmpFODNLaLDVUyPfuVOkY0fdUgSDnteJJsRwVjJhgn4QI0ZEfJbFxd6hGYxClHiGDRPpxfciIOfy\nn6gVGL8CP+OWOBw5UrvohWMm/ejLd9X2hZanM9SOkhIYN45qyw0Gg/r3yy/DsGF6GbeQK2JoGb9R\no/Qybzz1FCxdCnfe6e1TiL/lBrOaESNg9GiYOBHuucczaWEhtG4d/rhITRdOQ3yZPx/6WHJpFn2q\nHXNzoY2VjBP49vUh7YIlLw/q1dNr3LZnFW1YXXmOUnr9UeN7XDvsawk7KSvTgn/aNDj7bO3rffbZ\n+vfYsegT77wTDjsMTjgh4rW8Gva66JMfCxH95u+4Ay64QAv+F1/0zOuoo8K3saYBTTydO0M/vmM7\nDZhHl2rH4vr8/XQDkrUlwg/fbu8sKhIZkvepCMhQ3nW1/dZFL490oVZmlvvv14mmT/d9vXQMhpYs\nfHso7dwpMmiQDqgzdWrY/LJtEDzdKC4W+TRwlBQzMKbnT7ba8CMxY+pmqUDJ+JZjw9otTQWPjZhj\nzm/erD1LBg+O+prpFAwtWUQtnDdsEOnRQ6RJE5HvvqvMwzmwns0NaMoJBmV7XjN5Mvf/Ynr+RuB7\n0bWrfNP+5Kzz8kg0MWuJd9yhE82YkdTy1lVi8lBasUKkQweRtm3lgWuWhO0dZGMDmhYsWiQCsuhv\nT8T0/I3A9+L882V1vfbGyyMBRK0lrl+vJwqdfHJSy1mXidl0NmeO7G6yh/ysukgLfjM923Titdf0\nS/j665hO9yvwM27Q1hcDBtBm90raqdWuh80gVeyMHesxMGsjNOD41IHjCG7azA9n3ZGaAtdBYvZQ\n6t6dOw+awj6ylHf4I/lUH2E3ETRTyKxZOjJgz56JvY6fViFZW9I0/M8+EwE5Le+/xoafAkK9gNas\nkS00kkkMM7biKKjNAGu/fiKn85pUoORNTpEcdvvu2brZ/Q1x4qST9DhLjGBMOh5s2SKilHx81G1m\nkCrJ2IXV/Vwn5QSkMz+bhjZKYh1gDQ2sX8VDIiCPcYWEggl6DaxnStyitGWvvUQuuCDm043Aj8SB\nB4r88Y9mkCrJhAYc9+IX2UGePMXFZrA8RsLVXS9N3N7g3s1NIiCjucOzwTUumwlm5Ur9MB98MOYs\njMCPxPDhIu3aJe96WYxdADVvrmvdI1wpu8iVjiw2g+VxxC00uFI6xHJIMFf2DlRQnmO4CMj/1X8q\nrLaejXGLksqbb+qH+dVXMWfhV+DnJnaEII3p3x9eeAF+/RXatUt1aTKWMWN0uIUdO7SIANiHpYxg\nIhMZwVI6VaY1g+W1wx7awo4IzJ0LgwbBX/+qB9CHDoXx4xWPzH+K3qvW8Oi6y1ED2wIn1sjXaynG\neE77z1pmzNBhAKxlWBNJwr10lFJLlVKzlVKzlFKlib6eb/r313+//Ta15chg3GLrANzKHQQJcBej\nq6XPlpAIicIrtAVoL5xx4/R7KSyESZNgxrf16D3/NVTv3nDWWfD11zXOM3GLEsyMGdC7t3cQsDiR\nLLfMQSLSR0QGJOl6kenTR9fW0vRpgzINNwF0APO5kOd4jCtZSQfAJYiaISb8LIru6nrZpAm8+y60\naQMnn6x7vTZM3KLEUFIC551TwdZPS/lw48FJieOVnX74AI0bQ9euRsNPIG4CaAxj2Uke93IT+fne\nvvqG6PCzKHpYE0zbtvDOO7BlC5xzjl5I2yJcQELTSMfOmDFw7LHwwys/0zi4hZcWDeTYY/X+hOLH\n0F+bDVgCzAK+BS53OX45UAqUFhQUxDxoERN/+pNZzieBOAf7ujNbKlByD38TMAtrxBs/i6J7xjQS\nEXnxRZ3wxhtd8zcebbXH/p4u5ikRqLVrMunipQO0t/62Br4HjgyXNqleOiIi48frR7ByZXKvmyU4\nF9Z4lTNkE00qp/WbhTXiT1GRSF5eeIGfl+dj0fcrr9SJ33476eXPBuyK0GNcIRtoJoqKyncUQwzB\n9AmtICIrrb9rgTeBgxN9Td8MsIYUjB0/LBFjrntgX1ijLzM5k9d5gOtZT0tAV2/j4RFfxo6F6dOh\ne/fq5p1AAHJztUln6lS9/sMrr+BuRnjwQejXDy69FNasSWbxswK7qXMgM/iag7GL4o8/TpxpJ6EC\nXynVSCnVJPQ/cBwwJ5HXjIrQwK2x47sSsjNOnhxBQHgQWljjdopYT3Me5LrKY4EAtGgRe4NicKew\nEObMgS+/1KuM9e8PxxyjBf7u3VXCJrQa2X33wXHH2d7Bd3l8P+oFdq3fwqddLuPcYWLeSxwJjbU0\nYDs9mc3XDh04GKzypoo7froBsW7AvmgzzvfAXGC0V/qkm3RERHr0kPWHnmRihDiI1+zK4mKRo/O+\nEgG5ibur5ZObK5Kfb6brJ4NIC5Xbbfy5uXq9lOt4QATkUvWUeS9xJPRtHcbnIiB/YEqtJ7SRLjb8\naLZUCPyZvS+U1bQRZYsnYip3HGdXBoOypOAIWUNraaK2VJ6fl6eFipmunxwihVR22xQVMo1BspGm\n0pZV5r3EkaIikRsYJwLSmtWuzz8acehX4GevWya6y/TCjwNowxrasQqo6uYmrEtVR4jb7Mp336Xj\n8s/ZNmoMJ53TuNIN88gjq3n+VcOE6Y0/flw2nQgBruAJ8tjJg1xn3kscGTsWTm5dwlL2YS1tahxP\n1IS2rBb448dDyW4943YA1Qdus71yx2V2ZUUF3HQT7L8/ne6+jEmT9Pj4pEnw++9mun4y8Zo85cVC\nDuBubmEYkxkc/MC8l3ghQmHwS2bkHOZ6OFET2rJa4M+fD7PoTTk59Kf6wG22C524zK58/nkdxOXu\nu3WsEBtmun7y6dQpei0f4D7+xjw68yh/pm2zqqnTtfHgynqWLaP+b7/S+LhDkzuhzY/dJ1lbsm34\nodjg39NT3mVodBNUsoBYYq6HJucc0me7/Nagg2zpdrBIMOiazoTcTQ5uETSj3Y5hqgjI6Nx7pajI\nxMePBtdw1aEJbt99F5cJbZhB28iEhM4ELpV1tJTQQhBG6FQRTWW0C4FR/EME5Pi8T8IKgVgX8TD4\nx8/sW+eWk1N9wlxoe5s/ykaaSvv6ayU/3/1cpfTELvPtaMI1jDMGXCnSpIlIeXlcrmMEfgRCgqxt\nW5HL1ZMiIPuxwAidGLELlub8LuvZo7LX5NV4mun6icWvO6Z9GzxYC23n/q78KLvJkX9zdcQ8zDck\nMmGCe8MJIrPoJaUthsStvhuB74Gz1e3F9yIgf97jRSN0PPBaSckuWP7BKKlASU++j8mn2BA/YnHH\n7N8//HmhhWtCsV8iCf1s/ZaKisIL+yZskgqUjGFM3BpGI/DD4NbFDVAuW2gkj+ZenXUV1O/C1JFs\ntiEBsS8LpYz68gwX1hAihuQTGqfyK+xDY1fhztuTNbKJJvIyZ/vKKxsb+khmtCF8IAIymA/j1jAa\ngR+GcF3cjzlaZnBQVlVQvwNvfgZYQyaANzhVttBI2rGyhhAxJJ9obfih9+l13j05o6UCJd2Z7au3\nkG1EMqPdRpGUE5DGbI5bw+hX4GedW2a4CUUlFNKHWSybV5b8QqUAt9WoQpPO7rhDx1YJudl5raRU\nVgYjRsBnn8EgPuY03uJubuFX9qpMYxbJSB3hYtnn5mpP2XDugF4x8Ln+enbVb8ztOd4LGGSre22k\nhWgO40t+oBdbaQIk1wU86wR+OP/vGQykPrsZ0uq75BcqBXgJcRH46KOqQGmRZt3+9BOU7yznX1zL\nEjryANdXHs/PN4tkpJqxY/UCM2efXbXgzOef60bavs+5CI3bedOmwdGnt+B/+4/k9IrXOKrFbJRy\nv262NvRec0xyKGcgM/iKQyv3JbVh9NMNSNaWKhs+iLRllQjIkpEPJrwM6YDfwbyGDbW5JpId+Aoe\nEwE5g1er7Y8ltrchfQmZAZuzXjbSVF7jjMpga8a9VuNlDuvDTBGQc/lPDTNabcCYdNwJ11Xd3LAd\nG5sV0HF1dkwXbNnSX7qyMl0t69cPn6YZG7mDW/mUI3mdM6od+/ZbMxMzU7CbATfQnIf4C6fxBp3K\n5xMI6B6hWbLSXcaEOIwvAfgSHVJBKT0DOmn4aRWStaXCD9/u//3bMWfJmkYdMz5McqRVkWr0ftrq\nkLnhjt/PdVKBkj7M9PTYyGatLxNwDka2ZrXsIE8e5/Ks9cjxorhY93DtveOXGCYraC/2SZ7x+DYw\nXjrRUVQk8rd694ugw5VmqoCK1mtDqfD+xCByIHNlF7nyJCOi8gIx1D3czICPc7nsIE9aszorPXIi\n4Wwkl1Igkzkr7t+GX4GfdSYdN0Jd1c93DwT0smOZGibZa7DWDREd9DLMUR7lz2yhCbdwt6/8sj0K\naV3GbTDyfm6gPru4hoey0iMnEnaHhwKWsQ/L+YLDXdMm49swAp8qITiTfuwml0KqJHymCahILmPR\ncAEvcjSfchP38ht7Vu7PyQl/TrZHIa3LuEVQXUBn3uJUruRRrh2xNTUFS2PsjeRRfArApxzlmjYZ\n30bCBb5S6gSl1Dyl1EKl1E2Jvl4shIRgGQ34nt7VBH6mCSgvlzGloE0bfyF092AD4xhFCQOZyIjK\n/YEAtG9vQh9nIm6DkUrBA4EbacEGfrzx6YzqDccDeyN5FJ+ynubMpqdr2mR8G4lexDwHeAQYCnQD\nzlVKdUvkNWPBLgS/5DAGMoNcdgOZJ6C84tw3aAC9evnrAdzFaFrxG1fyGHbLYH4+3HprHGLpG9IS\nu29+27b6+/hKDqGYQg6Z+QiDjwlGtch9pmNvJI/kMz7nCMJZ0pPybfgx9Me6AYcAH9h+3wzcHC59\nqgZt7QOZZ/CqCMjBlGTsIKNXWGI//vkHq6+lAiX/DlwT1vfahD7ObJyD/+eh47sP4QPJy9NzNzLd\n2y0aSqesFAF5oMP90r27SH5+fL8N0sFLBzgTmGj7PRx42JHmcqAUKC0oKIj9jmtJSEC1U7+KgIxS\n/8xoARUuLHGkYFutmpfLoub9ZWerdvL11E2eoY1N6OPMxel9Up8yWU1reZs/GndcN156ST+Q0lIR\nif+3UWcEvn1LpVumSNVLWJa3v3zT/uSsFFCRAqUtvuFh/WPSpFQX1ZBC3HqCt/N3qUBJRxYbd1wn\nV1wR1wVPnPgV+IketF0J7G373cHal5YUFuoFtgvOO4IBO76g8OBgqouUdLyCZt32f6vpNHG0nlJ5\nzjmpLaghpbgN/j/BFQQJ8H88XiN9pnm7OYm4vu+nn8Lhh3u7sCWBRAv8b4ADlFKdlFL1gWHAlARf\ns/YccQSsXw8//ZSVCzWHC5p145pR2n/10UcJRczKxudjcB/8X0kH3uJURjCRfKpP9sg0bzc7Y8Zo\nHWjyZJg5E155pSrwIABr18LPP8NR7u6YScVPN6A2G3AiMB9YBIz2Sptqk04lCxaIgEw58TGzUHOI\njz/WD+Hvf6/cZRayzm6cA/MgchSfiIBcxNM1Zmxn4poIftaKkFe1I0gibVqkgw0/2i1tBH4wKDtb\ntpWXc84z4QFERHbuFOnaVaRTJ5Ht20XEZ0U3ZDyhca/mzUPvPyhz6CZfUVitTuTkZGad8FrspDK+\n0NVXizRqJLJrV8LK4Vfgm5m2bijFzIZHcEjF566H65o9stZml3HjdJf04Ye1sz6RF0WpS8/HEDuh\nca+qiI+KiYzgEErozpzKdHvumZlrIkRaK2LBAmDqVG2/r1cvqWVzwwj8MHwmR1DACjqypMaxumSP\njGhfjMTixXoJrDPOgBNPrNztq6Ibsgb7IO4LDGcn9bmMCYAe7hk0KIWFSyBeM9cDASjs8ItWloYM\nSW7BwmAEfhjW9TwGgGP4uMaxujL71msZQ19B4UTg6qv1enj/+le1Q5Eqel14Pob4YR/E/Z1WvMHp\nDOcF8iijQYPMnV3tNXM9Px+u6zlV/zACP70549ZurKYtg5la41hdCQ9Qa7PLG2/Ae+/B7bdDhw7V\nDkWq6HXh+Rjih9OddwKX0YINnFf/dUaN0mky0ZvLy4151CjYb8lUaN0aevRIbUFD+DH0J2tLm0Fb\ni1k9L5C17Ck5qqJOeqFECpPg+bg3bxZp316kd2+R3btdk5jwCQYnoUHcAf0qZHWjfWVT36OywpvL\ndeZsMCjSpo3Ieecl/Pr4HLTNTXWDk870vv5YuPhFbjh+DtPW9eKAA7RmW1cGnzp3hlmztBnHid3s\nUlKitf358/W+kSOh8NXbYNUqeO01bdKxYU9/6KHaRrt+PXXu+RjiT2Fh6P0H4J4RcMstvPXTfLaX\nVdn47GbFoUMzo75U3beNH2bDmjVpY84BjIbvyfLlWiV54IFUlyQm/LhOumlfA/NnSbnKEbn88hp5\nZoO2ZqhOSHuNOhjaqlVSrnLkPv7q7baYqdyvV9CTFSsSfimMH36c6NxZ5MQTU12KmPEyu7g1CIoK\n+YpCWcue8vX7v1fLy/jeZx+1beA/bnaqrKa11GNn9GbFus4JJ+j5K0nAr8A3g7ZhCPmuv7J+MGUf\nfsqMz3elukiuRPKxDxcmYexY90HdEUzkEEq4UY3jgWdbVDtmfO+zi1p7eQFf976MNqzlj7xT41hG\ne3OVlcFnn6WXOQeMhu+GXas5mbdEQE7I+zjtzBa11b6cg7otWSe/01w+4SiBYA3tq1aDwIY6h69Z\npBEo/qJcVqgO8j9OyK5e4bvv6pt8772kXA6j4ceGU6uZxrHspD7H7nw3rRY0j4f25fSlv50imrKZ\nq3gEUCxeXL3XYHzvs4t4TK4rPCyHhUdcwvF8QEe1DKjutpgJA7auvPMONGoERx+d6pJUx0+rkKwt\nHTR8N63mPY6Xn+iSVoNMcdG+bDb5Hvwg5QTk31xdIy8vm39WaGtZitdiOIFAFMHQli6VoFLyavcx\nngt+xDw4nG4EgyIdOoicdlrSLokZtI0NN7PFVTwkArI/89PGbBEv80pRkUjDBkGZyjHyGy2kOb/7\n8uoxvveZT1wb+OOPF9l777ALgGSU99fMmfomnn46aZf0K/CNSceBm9niXU4C4A/qXTp3To8Y8PEy\nr4wdCzPHvM2xfMy9DW5nAy1c04UGZb0GgQ2ZRaRZpFGZY0aMgBUr4MMPaxyKh3kyrfjvf/XkFFvs\nqbTBT6uQrC0dNPziYr3AsFOjmU13mcqx0rKlDvXqRxNJZBc1btpXWZnIvvuKdO8uB/XdbQZlDTWo\n7fqrxcUi55+1U37P3VNKOpxe4/x4mCfTioMOEhk4MKmXxJh0YqOoSCQ3t2bFu4e/yS5ypSkbPYVs\n6ONo29Z/w1CbstbavHLvvfrkjz6Kn83WYLCwm2r+wSjZRa50bLC6Wh3NKO+vZct0oe++O6mXNQI/\nBry05oEUi4BcwPNhBWL37uHPr+3gZrjeQq20r1WrRBo3Fjn55Ij3bwZlDdHirE9d+EkE5Ebuq1af\nvBQNECkoqEN1b9w4XegFC5J62ZQLfOA29ILls6ztxEjnpFrge3UtIShLKZAp/CFsxfSqtLXpoiZs\nQOuii0Tq1atWOc2grCFeuH1Pn3KEzOMACahg5XfgpWiAzqPO1MGDD9ZaWZJJF4E/KppzUi3wI3Ut\n/8kNspN6sgfrIwr2eHVR4611h3oEF3T5WgRk5QV/DZsmVputwSDi/j0N5zkRkCOZXu07cCo1dbKX\nuXixLuh99yX90n4FvvHSseHl+QIwmXOoz25O4e2YrxHtBCWvcAbbt8M55/j3Yqhc/epl4f/mXctq\n2tDv9dE1Vr8KLVtXWqr/ZuzkGENCcfueXuNMNtKMy5hY7TsIeX/tvXf4/NI+fMcrr+i/Z5+d2nJ4\n4adViGVDa/jLgB+Ap4HmYdJdDpQCpQUFBYltBiMQqWsJQVlMR3mXoTFr9361lJCW3aCBvzwjdXft\n9zaMl0RALuapuqE5Geok4b6nh/mzbCdfvv5gfY1z6toAbuXYWt+grGjaTTb3KExJOUiGSQeYCsxx\n2U4B2gA56FW17gKejpRfqk06Iu427Jycqgp3L3+VXeRKS9ZFJeijsUP66d5G25CE7KkN2SrL6SDf\n0F8UVQu71DnXN0OdwO17OiTfmpj00EM10qejp1g4hwn7dxpy6vhz/SdTMtaQFIHvdwM6AnMipUsH\ngS9S04Y9ZEiV8O3ObBGQkTzoS8jn5Ii0a+d/OvmECZE9fcJ9DF5CO6Q53UaRCMhhfJ7WmpMhc3Ab\nE9rSpZ8s3aOX9OsbrOF1lk6eYuEcJi65pHo5n+Ay2UpDacKmlJQz5QIfaGf7/zrg5UjnpIvAd+Ls\nZpZwsMymu0AwrABu2jTygKdbZbL3JqLdvB7fsGEi+6hlsp18eYlhaaE5GbKToiKRa+o9KgLSn29q\neIKli6eYV+Nj74U0ZKtsook8zUUp6zGng8B/AZht2fCn2BuAcFu6CnxnN3MET4qAHExJWOGbn+9d\nQSOPF0S3RRLaxcUir+ScI9toIHuzLOWakyE7CdX7pmyUbTSQx7m8Wo94yJDqExhT6Snm7aZdUx4c\nzmcp6zGnXODHsqWrwHcK58Zsli00kicZ4VkRvASp38rkZiKKqbv72WciIHfmjkm55mTIDtxs3/Z6\n/wwXymYaSyO21KjL6VAnIw0g6y0oczlQvqVvZY8/FT1mI/DjjLOb+RQXyxYaSWM2e2rd4bp2/ipT\nTaEesh1GJbQrKvQFO3SQGZ9sS7nmZMh8wtm+27Wrqs+H8KUIyBU8Fr0CE0fCDcpGmgEMIsfzngjI\n+bxQuS8nR6RLl+R+X0bgJ4BQxWjeXGQAeuLSNfzLs0KEu6VIlSknJ7xQj9TddVbghbc8pTN66aXE\nPiCDQbzNlfb4UhCUUvrJHLqJczwsWXZwZ8Nkd7QYPNg9kKJ9+5DB8gt7ua7Zm8wetBH4CSRUoT/j\ncFlMR8nBPcqkV9cukjfChAmx2TCdFXgPtVHW0FqWdThUL8xgMCQYL3Ol0yz5J54VATmGqb6VpXjh\nJ4WOy24AABnLSURBVKRDbq6OPuKmnBXylQjIX7nXs1HIy0u8pm8EfoIpKhI5p/4bIiBn8GpYwe3H\nSydeNnW3CjyO66UCJYfmf2vMN4akEMlc2bZtVT3NY4espZW8ySm+laXaYO/9FhT4G0fLy9ODyf37\na62/fn29fxqDZDWtpSFbI+bRvXv878WOEfhJoPiLcvm18X4ys95BkhMIxhTcLJ7eCE7Nqis/yi5y\n5QkuM5OrDEnDz+Sp4mItPAMBkTsYLRUo6chi38pSLMQyodFpXgqFTz+GqSIgf2G87zwSqXAZgZ8s\nJk4UAfnpH1NSPhhaXbMKyvscJxtoJq1Ym5QussEgEt3kqaIikf3yf5Hd5Mg/GBV2vKq2iwjV1g26\nf/+qPHLYLd/RW5ZSIPUpqxTokfJIpMJlBH6y2LVLZL/9RHr31t4wEnslrW3ltmtWJ/OW2DUQpXQX\nts4vEG2oE0RjriwuFvmy4BzZnNNMLj59Y43QBXahmZOjPdWiJRY3aGevJJTHSB4UATmN1yvT+Mk7\nkeLNCPxk8sIL+lG++mrMsevjEfM+pIHksUMW0Ulm011y2VWjUhr/e0MyiMpc+e23unLee2/luV4a\nuZfQd1OcYnGDtjcyEyboPNqzQjbT2AqgWN2zyMujJ9G++UbgJ5PycpEDD5Ttex8gezQo89WVtRPP\n+CFFRSJ31LtNBO35UKfjixuyi+OOE2nTRmT7dhk2LLIQdqu74RSn7t39mV3ctlDgw57dymUag2Qr\nDaUTi2oI9CFDwgv9RH9rRuAnm/f0BIybuCdsCx/Ohhepuzl4cOTLh7Sa0w/8UXaq+jKt7XnSv7+3\nJ4IZyDWkFZ98oivmI4/40siddddLccrLi+xTH2m7JUev/xwKK+4m0FMVB8gI/BTwcbNTZSsNa8Sq\nCW3hbi9S5Q4EvCtMZSWjQj7lCPmd5rJPgzVSVFT34osbsphgUKSwUKRjRzn/7F2e9dat7nopTn7X\nnA6rdPGh7CZH3m96pjRsEPQU6KmIA2QEfgq4+g9LZBsN5A1OFbeZg+FseH6mcIfrEtq1mlAQp4t4\nuvKcIUPSL764IXuJ6JgwZYoIyMKbJ3hGjlWqZt31o9zY3UHD5evc15UfZQPN5Ad6yJF9NqVFYDcn\nRuCngOJikdH17hMBGc5zvgR26LxImkc480tIq2nPCtlAM/mYo8UexGnw4PSKL27IXnw5JoS0/Pbt\n5Yo/bY9KAYpm8RQ300teng5rbj+vA8tlEZ3kV9pIR7U0bRUkI/BTxJi/l8vngSNkE01kH5ZE5aUT\nybXL7fH06yeiqJCpHCNbaCT7saDGOekSX9yQvXgpNfawyCIiMn26PvDPf8oll1QPxeC1cly0zg92\nTb1795o2/g4slwXsJxtpKgP4Oq0VJCPwU8jMN5bIttwmMqvRoXL+WTt9V5IhQ8IL+3Dml2HDRK5X\n94uAXMLEGh9SyPd+8OCq6eHp0g01ZA9+/OCrCfITTtBRCtevj8qEEoty49ZQFLC0UtgPpCTtFSQj\n8FPNyy/rx3vFFb5P8RtQzb4c4hUDv5My6ruOG9htkkarN6QSv37wlVr0rFm60l55ZdTX8tNAeMXU\n6c838ittZAPN5GBKpKAg/RUkI/DTgZtu0o/4kUd8nxJOQznyyJpxQFqyThbTUVbQvjJ8gu8PymBI\nIn4cE2qMVY0cqSv8jBlxLYtXTJ2TeUu20lAW01G68mNYU2q6kRSBD5wFzAWCwADHsZuBhcA84Hg/\n+WWcwC8vFznpJF2LJ0/2fZpTQ3Fb2DyH3TKVY2QHeTKAr30J+xoflMGQJKKJZVMpBjZtEtlrL5E+\nfUR2766Wl58QJG7pwpVDUSE3c5dUoGQGB0lrVnuaUtONZAn8A4EuwHS7wAe6Ad8DeUAnYBGQEym/\njBP4IiLbtokccYQOsTdlSkxZ1LR/BmUCl4qAXMgzvoV9jQ/KYEgibrFxIo5VvfqqPjBmTLU8Ipkq\nvWbcOjX7PVgvU/iDCMhLDJMGbKtzPeKkmnRcBP7NwM223x8Ah0TKJyMFvoh8PXWTLGwxQHaTI48d\n/HTUFcgZBfMBrhUBGcutUQv7uqKxGDKTSH7wrgJ2+HCRQEDmPPaZLy+cSB5B9t99mCmL6CQ7qSdX\n8ZDYXZrr0phXqgX+w8AFtt9PAWeGOfdyoBQoLSgoSPBjST4hTaMpm+QDtBvOvbmjZczfy33nEbJ/\n5rJLJnKJCMi/uEbcBmmNDd9QF4jKm2bzZpH995f1+e1kb5ZHNFX68QhSVMj1jJMy6styOshAiiu9\n2uqiJ1vcBD4wFZjjsp1iSxOzwLdvmabhOzWNXHZVmmI+CxwpM19b5Duf/fNXyFSOsWn20Ql7L/9l\ngyEVRDVj9YcfZEugqfxAD2nKRk9TZSSPoAK1XKYxSATkDU6tdHiIpAzFKzZ/Iki1hm9MOhJe0zif\nF2QzjaWM+vLL8L+J/Ppr+EzKykQeeki25+8hW2gkF6tnXLuogUDN9TftCzKnWwU1GKLlzqM/kl3k\nSgkHSwt+czVVFhdrLd21F0C5/EX9WzYHmspmGsul6imBoC/zTTzClyeSVAv87o5B28XZOGjrpWm0\nY6U8x3ARkN2BeiJ//KPIgw+K/Pe/Ih9+KPL889oHuVUrfcLRR8vMVxbU8N5xakjpGOfDYIgHxcUi\nZ9Z/W3aQJz/TWXoxq5qp8pJLwtnug3IMU6UU/UFuGHiczHxtke/vJJ7hyxNFsrx0TgN+AXYCa4AP\nbMdGW94584ChfvLLNIHvx/d4f+bLv3Ovkx3t93OvTWedJfLRRzrGiMGQ5RQViQzO+0xW0k52kCf3\n8jfZu8E6V2Ffj51yOq/JdI4UAVmhOsjkMyZH/S1FisKZDm7OfgW+0mnTgwEDBkhpaWmqixE3Skrg\n2GNh+3bvdIEAnH02TLp/FaxYAbt2QatWcMABkJubnMIaDGlASQmMHw/z50PnzjByJBQW1kzzzH1r\nOfmzUQxd/yLk5jK71SDeX92XDexBMzbRk9kcwec0YzMr6MCE5n/jpDdHMPCo/KjL1L8/zJzpfTzV\nYksp9a2IDIiY0E+rkKwt0zR8EX++x/YBJyfpPFBkMMSTaO3kxcUiNwydK8+3vkHm5PSU3egoa7vI\nlR/pKo9zuZzA/yRAea3mnkQThTNVYEIrpA9eA0lelSbdB4oMhngRrZ3cLTyCokIaslUUFTUcG0JB\nBGNRmowN3wj8qIkldGu6VzKDIV74sZOHertduojn4ihuW22VpnQPMe5X4AcSbVsyaAoLYdQoaNhQ\n2+xB/23YUO932inHj4cdO9zzKivTxw2GTGH+fC2a3QgGYfp0PR42eTLMmwcVFZHzVKrq/1DewaAe\nUxs3To8F+GXsWJg2TY+19e+v/06bpvfXJcyIYBIZOxaGDtXCesECPSbrNigFkT+ABQsSW1aDIZl0\n7gyzZum67UQpWLfOn5AHaN4c9t1Xn7Nihft3FFKa3L69cBQWRpc+HTECP8n4rTReH0AgoI8bDJnC\nyJEwZYq7R1sg4P4duBEIwAknwEsvaU18+XL3dNmqNBmTThIpKYFzz9UVccgQOO44/f+559bsXo4c\nCflhPMjy8+GaaxJfXoMhWXiZPFu3Dt/bdWL/Njp3rsrLSdYqTX4M/cnaMnnQ1mvRhUghXtN1oMhg\niDduM8X9TGB0+zayyfEBM/EqffA7AathQz0QZDf5hCaiRLL5GwyZitf3k5Ojv4s+fdy/jTFj9ABt\nWZk24wQCuhcwalTdG3D1wu/EKyPwk8C552rvgkiPunLG7aTklMtgqCvURnDXRmnyM/M3HTACP42I\nNDXbTvPm8L//pWelMhhSSbJ7u6FGZscOraylc+/ACPw04txz4ZVX/HsahHzz061SGQzZgpcZyc30\nmmr8CnzjpZMEvDxu3IhlYojBYAiP3UPOzSvOSaZOfDR++Ekg5HJmt0FGIpaJIXWB3bt388svv1BW\nVpbqomQN+fn5dOjQgXr16qW6KCnBaZqZNUv7/Hv1ojN14qMR+EnCOcu2eXOYMQO2bHFPX5crlRe/\n/PILTZo0oWPHjij73HdDQhARfv/9d3755Rc6deqU6uIknZISLeztphl7eIWhQ92Vqkyd+GhMOkmk\nsFB74JSWwkcfwUknZd/EkLKyMlq2bGmEfZJQStGyZcuM7VFFMtXEaprJ1ImPRuCnkEytVJEwwj65\nZOrzHjOmKqDazJnw8stw6KHQo0eV4I/VNBNtsMO6Qq0EvlLqLKXUXKVUUCk1wLa/o1Jqh1JqlrU9\nXvuiZh6ZWqniSbSDbYbswG6qsQt0EZg7FwYN0g1CbcIrZEqEzGr4mY4bbgMOBLpQcxHzjsCcaPPL\n5NAKXmTTwuM//vij77SJWgAmEAhI7969pVu3btKrVy8ZN26cVFRUeJ6zZMkS+c9//lO7C/vg0ksv\nlblz53qmefPNNyOmcRLNc68LeMXPt4dPmDAhO8IrkIx4+CLyk4jMq1WLY6hm2580yWj24K7BxRrL\n3EmDBg2YNWsWc+fO5aOPPuK9995jbAS1benSpbz00kuxX9QnEydOpFu3bp5p3nrrLX788ceElyWd\n8TLVhCgr0xq56UXb8NMqRNpw1/C3AbOAT4EjPM69HCgFSgsKChLbDBpSjl9N088KSLHSqFGjar8X\nLVokLVq0kGAwKEuWLJHDDz9c+vbtK3379pUvv/xSREQGDhwoTZs2ld69e8sDDzwQNp2dJUuWSJcu\nXeS8886Trl27yhlnnCHbtm0TEZGpU6dKnz59pEePHnLxxRdLWVmZiIgcddRR8s0331SW85ZbbpFe\nvXrJwIEDZfXq1fLll19K8+bNpWPHjtK7d29ZuHChjB8/Xg488EDp2bOnnHPOOa73nIkafqSAava1\nojO9F028ljgEpgJzXLZTbGmcAj8PaGn93x9YATSNdK1sNelkE34FT79+/j7kWHAKfBGRZs2ayerV\nq2Xbtm2yY8cOERGZP3++hOrkJ598IieddFJl+nDp7CxZskQA+eKLL0RE5OKLL5Z//vOfsmPHDunw\n/+3df2yU9R3A8ffHo9CEEpgTsKwwKnMNRWsRMIDgj6AbNkTHhKVEWQkEU4U5MHPBmPgjaCIso4SM\nbXGbIFWBiQMbwYgQFiL+2LAWqLTFlhULIpaqVcTwQz7743nKjttd77nePXft3eeVNDw83+/d83k+\nT+/T55577vvNy9OGhgZVVZ09e7ZWVFSo6qUFH9CqqipVVX344Yd16dKlqqpaVlamL7/88sXt5Obm\nXvyD8cUXX4Td53Qr+J2NhBltruh05LXgR72ko6q3qeo1YX5e7eQxZ1S1zV1+H2gC0vAmQ+OXVI1l\nfu7cOebPn8+1117LzJkzI1468dpv6NCh3HjjjQDce++9vPXWWzQ0NJCfn8+P3Z0oKytj9+7d//fY\n3r17M23aNADGjBlDc3Nz2G0UFRVxzz338MILL9CrV2Z8tabjhoc+fSL3Sec73brKl9syRWSgiATc\n5auAq4HDfmzLpKdk3rJ6+PBhAoEAgwYNoqKigsGDB7Nv3z727t3L2bNnwz7Ga7/QWyJjuUUyKyvr\nYv9AIMD58+fD9tu6dSsLFiygurqacePGReyXbp580pnrdtSoS08OMvoafRTx3pY5XUSOAhOArSLy\nhtt0E7BfRGqATUC5qn4eX6gmkyTrltXW1lbKy8tZuHAhIkJ7ezu5ublcdtllVFZW8p07kWq/fv34\nOuhr0ZH6hfr444955513AHjppZeYNGkSBQUFNDc309jYCEBlZSU333yz55iDY7lw4QItLS3ceuut\nLFu2jPb2dk6dOtWlXPRE48dDbS3s2QOlpWl0+6RP4nr/p6qbgc1h1r8CvBLPcxsTy6Tvsfj2228p\nLi7m3Llz9OrVi9mzZ/PQQw8B8MADD3D33Xezbt06pk6dSt++fQHnskkgEOC6665jzpw5EfuFKigo\nYPXq1cydO5fCwkLuv/9+srOzWbNmDTNnzuT8+fOMGzeO8vJyz/GXlpYyf/58Vq1axYYNG5g3bx7t\n7e2oKg8++CADBgyIL0E9UDpMMJ4MNjyySaq6ujpGjhyZ6jCSorm5mWnTplFbW5vqUDIq75nIhkc2\nxhhzCSv4xvhk+PDh3eLs3pgOVvCNMSZDWME3xpgMYQXfGGMyhBV8Y4zJEFbwTcbJyckB4JNPPmHG\njBkJec76+nqKi4sZPXo0TU1NTJw4EUjeKJvGeGEF32SsIUOGsGnTpoQ815YtW5gxYwYffPABI0aM\n4O233was4JvuJTNGWjLd06JFzkzRiVRcDCtXeuoa/MWotWvXUlVVxenTp2lqamL69OksX74cgO3b\nt/P4449z5swZRowYwZo1ay6+SwDYtm0bK1euJBAIsHPnTnbt2kVOTg6nTp1iyZIl1NXVUVxcTFlZ\nGYsXL07s/hoTAzvDN8ZVU1PDxo0bOXDgABs3bqSlpYWTJ0/y1FNPsWPHDqqrqxk7diwrVqy45HEl\nJSWUl5ezePFidu3adUnbM888w+TJk6mpqbFib1LOzvBN6ng8E0+WKVOm0L9/fwAKCws5cuQIX375\nJQcPHrw4xPHZs2eZMGFCKsM0psus4Bvj6hM0uHrHcMSqyu2338769etTGJkxiWGXdIzpxPjx49mz\nZ8/FoYy/+eYbDh065PnxocMqG5NKVvCN6cTAgQNZu3Yts2bNoqioiAkTJlBfX+/58cHDKldUVPgY\nqTHR2fDIJqlsmN7UsLynt6QMjywivxORehHZLyKbRWRAUNsjItIoIg0i8tN4tmOMMSZ+8V7SeRO4\nRlWLgEPAIwAiUgiUAqOAqcAfO+a4NcYYkxpxFXxV3a6qHTMmvwvkuct3ARtU9Yyq/gdoBG6IZ1sm\nfXSny4iZwPJtOiTyQ9u5wOvu8g+AlqC2o+46k+Gys7Npa2uzIpQkqkpbWxvZ2dmpDsV0A1HvwxeR\nHcCVYZoeVdVX3T6PAueBF2MNQETuA+4DGDZsWKwPNz1MXl4eR48epbW1NdWhZIzs7Gzy8vKidzRp\nL2rBV9XbOmsXkTnANGCK/u+07RgwNKhbnrsu3PM/CzwLzl060UM2PVlWVhb5+fmpDsOYjBTvXTpT\ngd8Cd6rq6aCmKqBURPqISD5wNfCveLZljDEmPvEOrfAHoA/wpogAvKuq5ar6oYj8HTiIc6lngap+\nF+e2jDHGxCGugq+qP+qk7Wng6Xie3xhjTOJ0q2/aikgrcCSOp7gCOJmgcBLJ4oqNxRUbiys26RjX\nD1V1YLRO3argx0tE9nr5enGyWVyxsbhiY3HFJpPjssHTjDEmQ1jBN8aYDJFuBf/ZVAcQgcUVG4sr\nNhZXbDI2rrS6hm+MMSaydDvDN8YYE4EVfGOMyRA9quCLyEwR+VBELojI2JC2qBOuiMjlIvKmiHzk\n/vs9n+LcKCI17k+ziNRE6NcsIgfcfr5P9SUiT4jIsaDYSiL0m+rmsVFEliQhrogT6YT08z1f0fZd\nHKvc9v0icr0fcYTZ7lAR2SUiB93XwK/D9LlFRNqDju9jSYqt0+OSipyJSEFQHmpE5CsRWRTSJyn5\nEpHnROQzEakNWuepFiX8taiqPeYHGAkUAP8ExgatLwT24QzzkA80AYEwj18OLHGXlwDLkhDz74HH\nIrQ1A1ckMX9PAL+J0ifg5u8qoLeb10Kf4/oJ0MtdXhbpuPidLy/7DpTgDAMuwHjgvSQdu1zgene5\nH86EQ6Gx3QK8lqzfJ6/HJVU5Czmun+J8OSnp+QJuAq4HaoPWRa1FfrwWe9QZvqrWqWpDmCavE67c\nBTzvLj8P/MyfSB3iDDD0C2C9n9tJsBuARlU9rKpngQ04efONRp5IJ9m87PtdwDp1vAsMEJFcvwNT\n1eOqWu0ufw3U0XPmmEhJzoJMAZpUNZ5v8XeZqu4GPg9Z7aUWJfy12KMKfie8TrgyWFWPu8ufAoN9\njmsycEJVP4rQrsAOEXnfnRcgGX7lvq1+LsLbyFRPXhM8kU4ov/PlZd9TnR9EZDgwGngvTPNE9/i+\nLiKjkhRStOOS6pyVEvmkKxX5Am+1KOF5i3e0zIQTDxOuJIKqqoh0+Z5Uj3HOovOz+0mqekxEBuGM\nOFrvng10WWdxAX8CluK8QJfiXG6aG8/2EhGXep9IJ+H56mlEJAd4BVikql+FNFcDw1T1lPv5zBac\nocn91m2Pi4j0Bu7EnW87RKrydYl4a1Esul3B1ygTrkTgdcKVEyKSq6rH3beUn3UlRvA0MUwv4OfA\nmE6e45j772cishnnLVxcLxSv+RORvwCvhWnyPHlNIuOS8BPphD5HwvMVwsu++5IfL0QkC6fYv6iq\n/whtD/4DoKrbROSPInKFqvo6UJiH45KynAF3ANWqeiK0IVX5cnmpRQnPW7pc0vE64UoVUOYulwEJ\ne8cQxm1AvaoeDdcoIn1FpF/HMs4Hl7Xh+iZKyHXT6RG292/gahHJd8+OSnHy5mdckSbSCe6TjHx5\n2fcq4JfunSfjgfagt+a+cT8P+htQp6orIvS50u2HiNyA8/pu8zkuL8clJTlzRXyXnYp8BfFSixL/\nWvT7E+pE/uAUqaPAGeAE8EZQ26M4n2g3AHcErf8r7h09wPeBncBHwA7gch9jXQuUh6wbAmxzl6/C\n+dR9H/AhzqUNv/NXCRwA9ru/OLmhcbn/L8G5C6QpSXE14lyrrHF//pyqfIXbd6C841ji3Gmy2m0/\nQNDdYj7naBLOpbj9QXkqCYltoZubfTgffk9MQlxhj0s3yVlfnALeP2hd0vOF8wfnOHDOrV/zItUi\nv1+LNrSCMcZkiHS5pGOMMSYKK/jGGJMhrOAbY0yGsIJvjDEZwgq+McZkCCv4xhiTIazgG2NMhvgv\nAu6q/XJn2FsAAAAASUVORK5CYII=\n", "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "cost=np.linalg.norm(y_data-np.dot(DataMat,params))\n", "print(cost)\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, cost=%f' %(len(params)-1,cost),fontsize=16)" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [ "x" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.6.0" } }, "nbformat": 4, "nbformat_minor": 1 }