{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# EE16A Discussion 12B" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Question 3: 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": 1, "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\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": [ "### 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": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "200\n" ] }, { "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'])\n", "print(len(x_data))" ] }, { "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": 3, "metadata": { "collapsed": true }, "outputs": [], "source": [ "degree=4\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": 4, "metadata": { "collapsed": true }, "outputs": [], "source": [ "D = DataMat\n", "\n", "# Least Squares\n", "params = np.linalg.solve(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": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "38.5961287395\n" ] }, { "data": { "text/plain": [ "" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEKCAYAAAARnO4WAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsnXecFEX2wL81s8sGyUhSWAGRaABRWdQTUURR7zxPRVA5\ns+dPPTFgQJFgREU5vDPcyalnwpzOLCh64IIuiBIkg4hKEFAQWMLO+/1RPbu9s909PbOTdqe+n09/\ndqe7urqmp+rVq1evXikRwWAwGAx1n0C6C2AwGAyG1GAEvsFgMGQJRuAbDAZDlmAEvsFgMGQJRuAb\nDAZDlmAEvsFgMGQJGSvwlVIXKKXEdmxVSn2tlLpKKZUTY17trDwuSFJxU4JSaoxSKmY/2kR/f6VU\nQCn1N6XUT0qpkFLqDY+0q2y/4R6l1Eal1Cyl1DilVLtElKe2oJTqoJTabr2LjukuTzSUUtcopf6U\nwPyeVEp9q5TaopT6zWrPf1VKBSPSBZVS1yql5iultln17HWl1ME+nhEpN8LHXIe0/ZRS05VSO5RS\nm5RSzyilWrrkW6yUel8p9YtVpnlKqcERae5WSn1o1XHH9qaU6qSU+rtSaqH1Dn5SSr2llDrEIW2h\nUmqCUuoHpdRO65nnRnsHXsQkONPEWcAaoKH1/9+BFsCodBYqTUwC3k93IYAzgWHA9UAJsDFK+g+A\nMYACGgOHApcCVymlhorI68krakbxCPArUJDugvjkGmA68FqC8itAt9/lgAAnAhOBjuj6FOYO4Cbg\nHuBjYG/gVuATpdQhIrLGx7PCciPMNvtFpdTvgA/R7ekMoBlwJzBVKdVLRHba0p4CvA48D5wD7AK6\nAfkRz/wrMBd4G/izS7kGAMcBTwGlQCPgRmCmUupoEZltS/sa0AcYCSwG/gQ8q5RSIvJs9FfggIhk\n5AFcgK4UHSPOfwz8GmNe7ay8Lkj390rTu0zo9wdGW/kFfKRdBTzrcL4+MAPYDrRJ8fvIS8NvcA6w\nDi1Eq9XrTDzcfrsEP2MysDXi3I/A5IhzXaz39pco+TnKDYd0U4BlQI7t3GHWvVfYzjUA1gN/8/Fd\nAtbfjm7tDd15qYhzjYDNwNO2c0c75YHuTH4EgvG874w16XhQCjRUSrUAUErlKqXutEwHu6y/dyql\nct0yUEpdbw2RmkecV0qpFUqpF6zPYVPIX5RSt1vDr1+UUv9VSrWJuDdqOWz5Xa6UukcptdYyVT1r\nDd86KqU+sIZ6y5RS50c8o5pJR2kTV4k1JP1FKTXT0kjiQil1kpXfDqXUr0qpN5RSnW3XV6G1dYBy\nt6FrNETkN+AKtNb3l4gy9FVKTbXezTbrnRwYkSZovd+flDaTfKyU6mKVZ4wt3Rjr3IHhdwu8ZLv+\nJ+udbbfe38tKqSKH93KZ0iaIMqXUz0qpfyulmvr5rkqpJsCDwHDgF/9vyTW/Q5Q2cWy0fqfFSqkR\ntutKaZPIYqsu/qSU+odSqmFEPsOUNrHsUEptVkqVKqVOt66tAvYDzlWVZpGnalp2BzYCeyLO1aP6\newp/TpTMKgY+EpGKZ4tIqVWe023pzgKaAw9Ey1BEQj7S/CyW5Lad+xVYAuwbUT6A9yKyeB9obbse\nE7VR4HcAyoHfrM//AW4GngZORQ+VbrLOu/EkEAIujDg/AGgPPBZxfgS6174IPfTsA0QOqWIpxwhg\nH+B8tGnqbOuZrwPvoCvcN8CTSqnuHt8Dq7xPAYOsfEqBt5VSJ0W5rxrWPe+g3+3ZwP8BBwLTlVLh\nyni69TzQ76GPdU/MiMjXaG3lKFsZTgGmWmU4D60ZNwD+p5Rqa7t9LHAL+n2fhh6ev+XxuDeBT4E/\nABOsZ10OvAosRJup/mJ930+VUg1sZRoHPIzWCv8A3ACcBLynIuzPLtwHLBKRZ3yk9UQpdQTajLY/\ncC1wCrozsSsgd1nnPgJ+bz3/AuAdpVTAyudctBCbDJwMnAu8AoQ7sdOBtWhzXPh3vsO6Vymlcnwc\n1d6N7d7GSqkz0G3gwYhkjwDnKaVOU0o1VEp1sM6twdZZR2G6Uqrc6uwec+icy9GmmUh2outAmKOB\nTcBBStvQ9yilvldKjfb520fFKtuBwLcR5cOhjGFT04HEQ7qHjT6GZp3Rcw1N0A2yHHjDSnOglWZM\nxL0jrfMHW5/bETE8QgutZdiGV2ib2be2z+H7pkXkP9w6v0+c5fg4It1r1vnzbOeaoDWf0bZzY/RP\n5j6ktN7Vh8CbDt+j2hAz4v5SYClVh7ntgd3Ag7Zzd3qVIyLPVXiYBdDCy/7OlwFTI9I0BH7GGlZb\n7+Y34JGIdNdF/g7hdwYMi0hbH21PfyLifHt0I7vG9u7KgVER6Y6y8v1jlO//O3Qj7RZRr+My6QCf\nAd8DhS7Xm1rPeyri/HnWc/9gff4HMCee3w441sor2jHN4d5TbddDwN0uzx5pvfdw2sXA/j7ez4lo\nZeBkoJ+Vz1ZgHpBvS/cFMCvi3v2sMu20nXsf2IEeYVxvffc70W1zgksZXE06LumfQ5s2O9rOnWzl\nMTAi7RPW+RFx1Z94bkrFYWsY9qMcrdE1tdJc4dR4qBRwf434fIEtzRHWuf7W59ZowXatQz43OlQq\nAYrjLMcNEenuts43jzj/IzDJ9nkMEYIW6IW2662zKmv4XS1yKIdrBQT2su6/0+HaNGC27XMiBf5M\nYKH1/wFWOS9Cd1z2479YAgo4xkrXLyKvItwFflFE2hOs88c7POsb4DUr3aVWuv0d0m3B1hE6fLd6\n6NHDOId6HbPABwrRbWCcR5qwoOgfcT4HXb8fsD6fb/3efwf649CBuP126BHXYT6Ozg73NrKuHY+u\n97uAuyLS/B9ayI5FC9gz0crIciwlK8b3dpr1Ti62nTvXOncn2gmkC7oz3QPssKX70Ep3XUSej1pl\nb+TwPN8CHz3aF+Aih99rIVoB64NWci623osAN8X6HkSkVnjpnI4eym0FvhORMtu18DDtp4h71kZc\nr4aIfKGUmg1cjh6qX4L+sZ1MMJsiPoeHVeFZ+ljLsTni8y6P85GeABVYJo6p6IrxV2A1+jvcAXR1\nu8+FJmgvmsjvAPp77Bdjfn5pCyyy/m9h/f23dUSy2vrb2vq7PuL6Oo/nRH6v8LOmuKTfHJFumUu6\nZh7PvAb9Xh9SSjW2zhVafxsopRqIyFaP+yNpgh7FeXmpONZFEdmjlNpou/40um5djFZYdiul3kUL\ntlVRyvEb2hslGlLthLZXl1ofpyqldgG3KaUeEZEfLPPGBOB+ERkdvk8p9TG6A7oBbcqKhbfQXjpH\nYNUrEXlOKdUFPVq/1Srri8C7VDWXhD3QPorI80O07OiGHqXGjGVSvBsYKSJP2K9Zv9eZaM+gz63T\n69AdxASc22lUaoPAny8ibo0tLIhboXt/bJ/t1914BPinZZ++BHhZRKLdk4xyxMtJaI1pkNhc1ZRS\nhe63uLIZXelbOVxrRRK+g1KqB3ouY5J1Kty4RuAsiMMdY7iytwAW2K47+lBbRAqf8LMuiMgjzNaI\ndAOo3iHbrzvRDf3ufnC4Ngf4GujhcX8km9Fa+b4eaex1seJ7Kb12pVn4umg18p/o+t8E/f0eQAu9\n3lHK0Rf4xEd5P0Vr6F6Uojux9uj31AnIo7JTwCrvJqXUcmJXZFwRkdus+ZkOwHoRWaeU+hbtihrG\nqW7UGKXUULT8eUBE7nIp30Kgh9LrVfZCT+yG10XMiOe5tUHge/GZ9XcweqIqTHhxwrQo908GxqN7\n0SKqT9amqhzxEhbsu8MnlFKd0PZlP77KFYjINmvEc5ZSaoyIlFv57QcciR76JwylVH30ROh2tOAB\nbaddBXQXkXEet89Da2xnUVXwnBVDET5HC/WOIuI1wf8RWsgWiUiklheNcVROcIc5CT2Zfx76+/pG\nRLYrpaajJzRvF5EdDslmojvGwejRX5iz0e19mkO+m4EXlVK9qeoxtRPnNQOzgcN9FNnP6KUvujNe\nYX0Oj4oPQ0+0AxUTmx3RHWWs/BEtMGdFXhCRbej6FHZa6IIe9YR5Az1iPjGczuIkoCzinC8sT6gn\n0eba4dHSh0dcSnv8XQV8KCLLPW9yoVYLfBGZr5SaDIyxNJjP0fau29B+vJ4/hojssFzNrgXmicjn\nXumTVY4aMAVtwnlaKfUA2tQxFm36iMcD6za0x83bSqlH0BObY9GTm1Hd0jzYWylVjDYZNaJy4VVz\nYIiI/Aha61RKXQm8qZSqh/bI+BmtuR8JrBaRB0Vks1Lqb8AtSqmt6PdwKJUN1Y973Bal1A3Aw0q7\n575nfc990UJomog8LyLLlVL3Av9Q2j31U3RDb4ueB5gkIo7arogsotJcBWjXXOvfWfaRq9KurU+i\n5yWmeRR9uFWGEus3X4PWUHuIyF8tTfgBYIRSahvaRNEVbauejuVRpZT6F1ogl6BNY52AoWhTRZiF\nwO+UUqeiBfHPIrLKMkNV0cCjYXlfXYiei1mNngcYCFwG/NNWB1Yppd4GblTaBflT9MjkRrTm/6gt\nz1FoL7f9ReQ769wHVJo5d6KVn+Ho0dRztnt7Ws8PdyBHo81F99nlgNW2nwJutzyc5qDnPC4B7hDt\nXhzOsy+6TodHyYcp7QaMiLxipTkGrWh+DTxltYswO0XkK1t+I4Dv0HN5RcCV1t+jiJd4DP+pOPC/\ngKIeujJ/h9Z0v7M+59rStMN9IUQf69qVDtfC910Scf5Y6/yxcZYjMr8x1vmciPOrsE2a4TxpOwgt\nVMrQw8/BaK1ylZ/v7/CdT0ILgR1oAfgmEZNvxD5pa5903wx8idZ+93O5pw96Inqz9b1WAS8AfWxp\ngujR1FqrrNPQnYJg88hxe7e26yejRwlb0KONpWhPiG4R6YaitedtaBv2t2hPl5gWjeG+oPBK63xX\nH3n0RAvOX6zvvgjbJB66Y70WPYLYhTaBPQw0tKU533pn69GCcSXaNmxP0wX4n/VehAjPnxi/dxe0\nN9r31vPWoTugc4lYwIceud6GFtrbrPK/Axzh0m7a2c79zfpttlrffTl6FN8o4t7u1vPD73AOcGEU\nGfO9lecSIry+rHTTcPFYciiz07EqIr87rd8l/L7+A7SN9zcQEe2SmM0ope5C+9bvIyJb0l0eQ/xY\nk1wvA8eIyP/SXZ5YUEo9DzQWkZPTXRZD3aVWm3RqgjWk64wW9v8ywr52YdmbT0HbZcvQ7qk3o7Xw\n6R63ZirHoEdrBkPSyFoN31o63hK9knCoxOYeZ0gz1grkh4GD0Auz1qPNHCNET0IaDIYIslbgGwwG\nQ7ZRG2PpGAwGgyEOMsqGv/fee0u7du1iv/Gbb6BhQ4jnXoMB+PFHWLcOQjaHzkAAWraEffZJX7kM\nhips3gwrVkC3blBQuURi9uzZP4tIc487gQwT+O3ataO0NCb3Xk3fvrql/q9WOWYYMoSZM+H446sK\ne9Cff/0VXnsNiuMKRmswJJj77oObboKSEq3kWiilvvNze41NOkqptkqpT5TesmuBUmqYdb6pUuoj\npdRS62+Tmj7LlQ4ddK+HbrxDhkCvXvrvzJlJe6qhFmOvJ2efDTuc1qwCZWUwcWJqy2YwuLJqFTRt\nWkXYx0IiNPw9wPUiMkfpGOKzlVIfoReYTBWRcUqpm9Euczcl4HnV6dABfvyRO27ZwbiJBezYASIw\ndy689RYMHw5jxyblyYZayOjRMH48FfXEi1AIli5NTbkMhqisWlUj03WNNXwR+UlE5lj/b0WvctsX\nHZI0HKPkP+h4FsmhQwcAXrx3Fdu3VzbiUAi2b9eN22j6BtD1YPx4qtQTLwIB6NSp8l4zejSklXQL\nfDtWnJCe6MUwLUUkHNVwLd6RDOMi3ABPvVoL/P1CKxzTmWG5IczEie7mGyfy8+Hqq/Wo4Pjj4cUX\nYc4ceOkl/Xn06Oh5GAwJQSRzBL4V/fBV9E5BVVatinb2d9SnlN4rtFQpVbphwwbfzws3wBdegC82\n7Q9AB5wFvhmWG8IsWeJfsy8s1OZAqD4qMKNHQ8pZv15rK+3bx51FQgS+FbbzVeA5EXnNOr1OKdXa\nut6a6ptVACAi/xKRw0TksObNo3oVAVWH5QAbaM5v7OUq8O3DckN206mTrg9OKAVFRdpkM2gQTJ2q\n5368RgVm9GhIGatW6b/p1PCVUgq9i8y3ImLfjPgtdEQ+rL9vRt4bL9UboGIFHdgf5xDR4WG5wTBs\nmK4PThQUaJNNaSlMnlzpiuk1KgiF4IMPjF3fkAIyQeCjYzMPBY5TSs21jpPR4W9PUEotRceP9trQ\nIiacGuAKOjhq+MGgHpYbP2oD6HowfLg214Q1fbv5xqmeeI0KQK+FMXZ9Q9IJC/z94t9ttMZumSIy\nHR1/24nja5q/E506aZdL+0KZFXRgQMV+w7o4wSA89hhcckkySmGorYwdCwMH6pHi0qVwwAFa83dT\nCoYN0+69YROiG3a7/sCBRskwJJiVK7UPfoMGcWdRK2PpOA3LV9CBQnbQknUopTW2W281wt7gTHGx\nNttEmm/c0kaOCrwwdn1DUlixAvbfv0ZZ1EqB79QAV6BfxBHNVnD22ZUTbgZDIhg7VtepQYO0vb6J\nx7rxUAi++sr9usEQF8uXZ6fAh+oNsMvJ2hf/rQnLo2psBkM82EcFJ57ore0vW2Zs+YYEsns3fPdd\nxSLTeMmo4GmxUlxsE+w720FhQPeCBkOSiWbXLy83tnxDAlm9WleqbNXwq5GXp52ozQorQwoImxWD\nQfc0xpZvqCnhaAJXnqQV2QVlRuBXcsAB2mfTYEgBY8dCx47u180Kb0NNsIfzCC3TAv+P1+9fI1Nh\n3RP4S5f6WztvMCSAnj3dbflmhbchXiKD/O3PcsrIY3nZPtx7b/wL/OqewP/1V/j553SXxFDHcIuU\n6bVy16zwNsRLZDSB/VnOCjogBNi5M35387on8MGMow0JxStSppuPvlI1inFlyHIiownsz3KWU2m/\n//bb+LR8I/ANBg+c4udHRsoMuwh37Vop9EV0o3QLteA3tr6JwZ+dVDUFCh1YUUXgh0LxOQTULYHf\nvr12mzAC35AgRo1yd72M9MJZubJquA+3EMp+Y+ubGPzZy7BhepQI0JJ11GdbFYEP8Ym5Wi/wq2hA\nf85lR+v2RuAbEsLo0TBlivt1uxeO3xDKfkYMsaQz1E2Ki6FbN/1/OAqwXeDH6xBQqwW+kwb02Y8H\n8NNnxjXTUDPCAtfL4cve6PyGUPa7YbqJwW+YNElP/DsJ/HgdAmqtwHfTgBaFDqDB2qXMLDGumYb4\n8bMVYm6udgjr1Uv/VW4xY6kMobx6tXfHEB4xROtAzCC27lNcDDfeCF1ylhNCsYp2UUN5R6PWhlZw\na5BLOYD6bOM/966l+I3WqS+YodYxc6auT0uWaI192DB/WyGGQtrkI6KFfU2Xf9hHDE4hwJ3SGeo2\nY8fChlnL2fi/thzUNS9qKO9o1FqB79Ygl6I9dfZ8uxQwAt/gzejReqS4Y4euT3Pnwuuv60gdXgQC\nOp5VGHtdjFf424fpXrF6jH9/dtH81+XQe39KP655XrXWpOO2C1FY4B/e2Ix5Dd64mQV37oQtW9zv\nCwbdBbpS0LZt9BDKdpyG6fHszGWooyQgLHKYhGj4SqkngFOB9SJyoHVuDHApsMFKdouIvJuI54G7\nBrSaInaRyymdjcDPdpxMNXZB6cdOH2QPBzGPTiyhAytpnbOB+vV2sWN7iE00ZQPNWUInvuFgfqI1\nIormzXUI5SFDtCOBk1kmjFLaf3/SpOpCPNaduQx1kK1bYcOGzBL4wFPAP4CnI85PEJHxCXpGFcIa\n0Pjx2mshFNIaUF5+Dr8WdGDf34ynTjbjZKp56y1dZ8Ib47iZBVvzI2fyCn/gLYqZSX22VVwrzylk\nx548dqJozC8EqZTm39OG9xnIlsJTYNdAhg2rF3VrRBHtv+9GlRDghuwjHO49QQI/ISYdEfkM2JSI\nvGIhchOUQYP05+ZHdjJuDFnMpElw113RfdjtZkFFiIG8ywcMYA1teIhhtGItT3Ihg5nMQXzD73ps\nJbhjG/M/3URR4UbqsYvmrOcYPuVqJjKL3pzNC1z/vz/CvvtS/Mpw7rx4JYWF3uU1bpYGV5Yt038T\nJPARkYQcQDtgvu3zGOA74BvgCaCJy32XAaVAaVFRkSSE664Tyc8XKS9PTH6GWsOoUSLBoIgW9c5H\n//46bUmJSGFBSP7EK7KQLiIga9hHRjFGOvNtlXsCAZEhQ6o+p7BQnw9fLywUGXPrLpF33xU54wyR\nnByRYFDWnXKh9Gy4zLNMvXql530ZMpy77tIVZMsWz2RAqfiR034S+cqousBvCQTRo4i7gCei5dEr\nUbX+0Uf1V/vuu8TkZ6gVlJRooeslWMPCedQoEfniC/l+nyNEQObTTYbwnOSwy/GevDydf+TzBg/W\nwnrw4OrXZc0akWHDRPLzZVegntzNCCnkN8fy2DsTg6GC888Xad06ajK/Aj9pXjoisk5EykUkBDwO\nHJGsZ1WjSxf9d9GilD3SkH78TMIC5Ie2sffd1yHFxbRhDctvfYK7Bn3Dkl7n0Ll7ruMuVqEQvPde\n1XP2PW6d9lGe+f2+DFn3N048YAWftBzCCO5hEV3oz0dVyxOjm6UJqJZFLFkCnTsnLj8/vYKfg+oa\nfmvb/9cCL0TLI2Ea/tq1WnWaODEx+RlqBYceGl2778WXsowOIiBv7Xu5XPinX+TQQys19JISbQ10\nurew0EGLdyFs8lGqUos/JjhdFtJVBGQCw6RQbZfCQmu04ROnfGPNw1CLaNZM5LLLoiYjlSYdYDLw\nE7AbWANcDDwDzEPb8N+ydwBuR8IEfigk0rixyOWXJyY/Q0oJm0rsgtgPgwdX2tSrHyG5kr9LGfXk\nO9rK7/hUAoHqgrN798pzTqaXwYP9ld/NtNSo3nZ5ve1fRUBWNT5Y5ry6Iqb34pZvLJ2RoZbw88/6\nxx0/PmrSlAr8RB0JE/giIn36iPTtm7j8DCmhJhqsm0DMY4c8yzlaq+dUacJG1xGAm7CPZXJ18GAf\nnca774o0aSLStKnIlCm+3o2vfA11h5IS/eO+9VbUpH4Ffq1daRuVLl2MDb+WUdOQwE6rU5vxM1Po\nz7k8zy3cxWm8yWaauuYRfq4TfmPY+ImcObPJQPjiC2jVCk48Ef75zxrnazyR6xhLrLVECQycVHcF\nfteusG6dDlNoqBUkIiSwfW3GHw5cwex6fTiMUs5WL3EPt6ACAc+oluC+KbnfyVW3sB9hNm+2NjJ5\npqPuxU46CS6/HO6807PHiZbvhg1mArdOsWSJjuORwL0y67bAB6Pl1yISpcEWF8PkMYt5c/Mx7Ndg\nM8v+9Qmbjj+LJk2gUSNo0cI9lHEgoDeeqEkMG6+NzcNUjFoWNNDR2oYOhdtug2uucY3FEC3f7783\nO2LVKZYsgQ4doF69hGVZdwV+2DXz22+NG1uGE/59vEIMBALQtKnP33HBAujbV4eznDaNl9f04fPP\n4ZdftHa9fr17x5KfD48/7ryCOxySIRp205IXFaOW3Fx46im49lp46CE9jHAooD1fpw5LxOyIVacI\nB4FKJH4M/ak6Ejppu2ePSL168r8+Nxg3tgwmcpLW7cjJ0e6SUX/HefNE9t5bL1ZZuDDqYqxk1ouS\nEj0v63sSOBQSGT5cX7juOv3ZJd+iIvc8zQRuHaC8XFfIa6/1lZysn7QNBtnephNbZn1r9gXNUJwm\naSMJBHRseqW0Ruz5O65cCQMG6CHwp59C166e8wL2UMaxavF+KC7W87FudvfISeCZsxRDvr+PF5pf\nBQ8+yA8XjnTNd++93Z9rJnBrJ3ZLxJWn/6gruNHw/VPS9ixZQkejBWUoXm6GoLXjwYNFTjjBhzvi\n2rUi+++vb5o/v+IZ0RZjJTuGjV/feftIR1Eu/+JSEZB3B0xwfXduaw5MqIbaR+RI93g1VQTkyaFT\nfd1P1mv4wLzdXejACvIoq3bNaEHpJ9o2gi1a6InK0lLvydwfF/2qVemffoJ334Xu3Suue3m2xLtV\nYCxzQn42Mokc6QgBLudRXuN0TvzwOhaPe61avl4TuGZHrNqF00j3AFkMwF0vd0qsJcJPr5CqI9Ea\n/kN9nhcB6c48owVlIN4rY7W2k5vrraHnqt0yt9WJ2sj//vvVnpHo1anxLgzzCrTmNtIpYJt8TrHs\nDOY7FtQtYqeZn6pdOP3+D3KN/EahBFW5L0sEWb/SVkReu22OCMiZvJSQxm5ILH6jW3odD+dcrf+Z\nNMn1OYkSjMkKbeBldtqb9bI6b389Eb2iehiGqBE7DRmP0+//NifLHHpUmDaj/a5ZL/BHjRJpVrBN\nBGQkt1fRGo0WlDn4iV/vdvw11wqDfd11UePvJEIwJiu0QTR7/LWnLtGxoQ45RGTbtvgeYshYnH7/\nJXSUFxhURaHwkllZLfDtmthK9pPnGFLx4oJBkccfT8hjDAmic+fYhX0/psoeFRQ5+WQZPXJPNddO\npXQgtERqvH4ngP0Gfwun69zZvdOrGDm8+67+Uuec4+quaaidRI4c61EmuwnK7Yz0PYrMaoFv18Te\nYaB8xSEJ0cQMySGaLT/yaMNq2UAzWaC6yVMP/eppFsrPT9xozo9njF8bv9v6A7f7SkpEJh+kdz96\nuscDxnRTxwjXBxDpzjwRkME873sUmdUC366J3csNUkY9CbK7miZmyAxiseXnslNm0Ed+pYF0Voul\nqCj6oi27ZhRv6OVo5Sws1CNHPzZ+r3yCQZEuXaqWraJzICSv8CfZQ0D6531mzJJ1jPBCvbN4UQTk\nEL5yHUVGktUC366JDeU/IiBdWFhFEzNkFk4Tqzk51QXig1wjQuVEvB/7f1gzSsTmIV4TwH5t/LHM\nBUR2DvXZIkvoKKtpI20Kfnadq4inQzOkn8GDRUarMVKOkny2O44inchqgW9vJD2o6qljvHMyF6eJ\n1RNOqKzwf+IVEZC/cbVv80/46Nw5cR42bhPAfm38sSwGc+ocejJbdpIrb/F7GXx2pT3f7IZV+ykp\nEXk5OEiWsn9M9TSrBb5IZeUvUDtkDwG5ndtM5a+FhDvvdqyQX2goJfSWXHbGJOwDAfE0/SRqXsfv\n6tdYVslYQ8W7AAAgAElEQVS6dQ5X8zcRkPvaTqzynhLRoRnSy9oWB8rbwd/H5EacUoEPPAGsp+qe\ntk2Bj4Cl1t8m0fJJeGgFSxNbmddZZrU53VT6WsqYkbtlRuAo2UwjKWJVzNp9YWF0T6BEVD2/QtdP\nunDddQ++FpL/cqrsCtQTmT3b7IZVV9i9WyQ3V34478aY3Ij9CvxEhVZ4Cjgp4tzNwFQROQCYan1O\nKcXFMHkytPv9QRyRP89XLHND5jE6bxxHhmbwTPEjNO+1Hw0a+LvPHsKgZ8/Eh1iIxE8YBT/p3ntP\nx7V/8UWv/XsUVxQ8iTTdG4YOZcHsMrSeVR0TRqQWsWIF7N7NPv27MXmyDisyebK/fRh84adX8HMA\n7aiq4S/G2rgcaA0sjpZHojX8CsaM0eqPWbRSawhruEM7z5I9KigbTtB2jmgLtZTS5ptIzSiVJg+/\ni7zsfvhFRfrv4YdHd1GtMsR/7z0RkPvVDZ7pjaNCLeH11/WP9sUXMd1Gqm34DgL/F9v/yv454r7L\ngFKgtKioKI435INXX9Vf9csvk5O/IaGE51/qs1WW0FFWUSStCzbLRRdFd9/0Et6ZGHvG734A4SMc\nQdTekU3KuUzKUXIk0x3vycvTk9/Gc6cWcJdeayFbtsR0W0YJfOvz5mh5JE3DX7xYf9Unn0xO/oaE\nYdfEH0MLsmOYJqA1ey/BGAzWLIhZqoknllCTJlUF9+DBIg3YIitoJ0vZXwr5rdqIJzfXeO7UGs47\nT6Rt25hvywSBnzkmnT179JLL665LTv6GhBGefDyej0RA7mO4b2HYuXO6Sx8b0fYD8GPWad1afz6G\naSIgf+fKagI/FWYsQ4I49FCRE0+M+Ta/Aj+Z8fDfAs63/j8feDOJz/ImGIRu3fhlxnyzt22Gs2QJ\nFMpvPM6lLKYTo7jd132BABx6aJILl2Ci7QfgRXjHr/Xr9c5dn9GXCVzDVTxMX6ZVpHPLv2I/XUPm\nEArBt99C165Je0RCBL5SajJQAnRWSq1RSl0MjANOUEotBfpbn9PGV+UHsf2L+bz4IsyZAy+9pD0h\nRo9OZ6myj2ibh3TqBOMYwX58x0U8QRkFFdeU0n23E7Vx0w+vzVn8EgpV5nErd7GcDvyTv5BHmeNG\n5/b7Ij13YtnYxZAEvvtO78fZrVvynuFnGJCqI5EmHfsS8/79RW7OuV8EpAkbzdA2TfhZCTr/0c9E\ncF5NW1goFRO3mTTxGi+J2A8ARFq1qnwn/flQBOTunNuke3f/C7zMKt0M4J139MufPj3mW8nmlbZO\nng8DeF8E5Hd8Ws3GaRalJB9fbpHbtol07CgbG7eXvQt+cxXq9onX/v1rtwdKpOeQk53ej+C2v5NP\n2w2V8pxcmfvc/IQtBDOkgPu1UiobN8Z8a9YKfLfKuw9rRECu5O+OGpIhufhaCTp8uD4xdaovb5q6\nopU6dWD27x2zQF6/XqRZM5Ejj5TRt5VHHRF5/TZGIUohF14oO5u0iCv4XdYKfPfKG5L17C2Pc3G1\na8Gg0WKSTbSAYYO7ztU/xKWX+sov27TSmNcQ/EdHiZVHH43aeUb7bYxClBrW7HOYfBw4Pi4FJmsF\nvlfl/YATpJRDjRaTZJxC9HoFDAuqclnSrFikeXPfw9lsjB0T0xqCUEjk+ONFGjYU+fFHz3yjuYca\nhSj5lEzfI9vJlwe4Ni4Fxq/AT6ZbZlrw8nz4ip4cyHxy2VXlvIiJNZIoRo+ujANj94YqLNSeNE5c\nnjuJAzbOhAcegKZNfT3Hy6WxrsaOCceG8hVfRSl47DHtf3njjZ75Dhvm7S0UChkXzmTz4t3LKaCM\nbzi42rVEutDWOYE/bJi7YPmKnuSxi24srHI+UcGzsp2ZM2H8eO0fHhbGYX/xF16AwYOrBwzbr2A9\n9wdugmOPhfPO8/0sr47d/J4WHTvCDTfAs8/C9OmuyYqLoUUL92yMQpR8cr/9BoB5HFTtWiIVmDon\n8N0iEeblwfycngD05Ksq9yiltVHje1wzJk7UbsROlJVpwT91KgwapH29Bw2CL/sOp6B8Gzz6KJ6O\n4xF4dey10Sc/Hnz5zY8YAW3bwlVXQXm5a159+7q/ftOBJp/ehd9QToCFVPfBT+j792P3SdWRDD98\nu71z9G3lspW95CH+6mr7rY1eHplCLDs5iYjIxx/rC7feGtfzMjEYWqqIyUPp5Zd1oocfds0v2ybB\nM42Nx/xRFqnOcb9/snXSNhpbDjpSvt37aOnc2T3Mrqng8RHLTk6yc6cOftOhg8j27XE/M5OCoaWK\nmIVzeAK3cWPtsinOE+vZ3IGmnQ4dZF63s+J+/0bgu3HllSL168uQs8uzzssj2cQkiMKLTN55J23l\nra3E5aG0YIHeFf7SSz1HB9nYgaadrVv1D3H77XG/fyPw3Zg0SQTktG5LYjM/GHzhS0v86SeRBg1E\nTjklbeWszcRsOgtz/fUSUkqOzvvCjGwziZIS/QO88UbcWfgV+HVu0jYqPfXEbb/GXxkvjyQwdmz1\nidmpU/X5CkaMILSjjGuZYAJ1xUHcHkqjRvFLXkvG7bwGqO7TaiJopolvtIcOB1d3yUw4fnqFVB0p\n0fDLykRycmTN0JvNJFU6mDVLBOSBnBtrfUiEdFGTCdaxRXqEeyYvxTyydbL7GxLAVVeJ1K8vUl4e\ndxYYDd+FvDzo3p19183xteG0IYGEQmy96GrW0orRe0YiEb7648cbTd8PfjdLd2JR8QV8zcHcy03k\nUVblmtfowG1BnQkvngC++QYOOqjmsbL94KdXSNWREg1fROSCC/Qy/lDITFKlkqeeEgE5n6fMZHkC\ncKu7Xpp4SYnIKXl6N7Hh3OdrdGBcNpNIKKT3rbzsshplg5m09eChh/RX//771DwvyykpEbngjC2y\nIaeVfBnsLYpyM1meJJxCgyulQyyHBfOoUSLvBk6RX2goe7M+qkktG+MWpYxVq/SLfOSRGmXjV+Bn\nn0kH4PDD9d8vv0xvObKAsCmg66t3sveetVxR/hBu1c5MltcMp9AWoP9fsAD69dO/x9ix0Pq5+6mv\ntvFw87HOE+s2sjFuUcqYM0f/TdH+nEkX+EqpVUqpeUqpuUqp0mQ/zxc9ekBODnzxRbpLUqcJC6B9\nti/lGibwJBfwJUe4ps+WkAjJwiu0BWgvnPA8SY/BXQlecTmDNj3G5FHfetr9TdyiJDJnjt63MxUe\nOqQulk4/EekhIoel6Hne5OfDIYcYgZ9kwgLoPm5kJ3mM4B7HdGayPDH42RS9iuvl6NFQv74OsOaB\niVuUHGbOhDlPfMWy3C4MuaggJQ4L2WnSATjiCG3SCYXSXZI6y5IlcLR8xum8wThuZh2tqlzPz/fw\n1TfEjJ9N0auYYJo3h5Ej4Z134KOPXO+piVeQwZmwqbPVj3P4vOzQ1Hk9+TH01+QAVgJzgdnAZQ7X\nLwNKgdKioqIaTVzExJNP6smShQtT98wsY8jZ5TKLw+V79pUCtlWb8DMbayQWP5uiV4tpVFYm0q6d\nSM+eUf3AjUdbYgj/Tq34UQRkGBNq7PVEpnjpAPtaf1sAXwPHuKVNmZeOiI4tAtpV0JAUlox5TgRk\nKP9xFD5mp7HEM2qUSF6eu8DPy3PY9P3ZZ/XF559Pd/GzgrDX00DeEQH5HZ9W+Y369489T78CP+km\nHRH5wfq7HngdPGbtUkmXLtCggbHjR8FXzHUnyso44IkRzMvpybM4b2wiYjw8Es3YsTBtGnTvXtW8\nEwhoP4VQCKZMqbp4asziIXpOa+RI2LXLNW9DYgjPtYT35ZhLjyrXP/44eaadpAp8pdReSqkG4f+B\nAcD8ZD7TN4GAds80At+VGq2ufOghWL2aN383HpS7G2bTpnF2KAZXioth/nyYMUPvMtarFxx3nBb4\nu3dXTuyGVziPuy/ALYFxsGIFT/b5FzNn1qCjN0QlPNdyKHNYSke20rDK9VAoiavO/QwD4j2ADmgz\nztfAAuBWr/QpNemIiNx8s0hursyctsPECImgRqsr16/Xm2efeqpnPjk5Ivn5PjfxMNSIaBuVQ0g+\n5lhZSwtpHNwiubnmd0kW4TaxgnbyAoMSsqCNTLHhx3KkXOC/9poIyDF5M03ljqBGqyuvukrPyFoT\n4k4hk/PyRHJz4+xQDDETLaQyiByODmw3ijHmd0kyd9+wSQTkJu5JyKpzvwI/e90ygdlBPZ3QY+dM\nJGKYm+2BvOJeXbl4MTz2GFx2GXTtCjiHTD7mGNizxzkLE6Y38fhx2fySI3iZMxnOeJqzvtp187sk\njhEnafv9V/R0vJ6sBW1ZLfDHT96X7yjiSD6vdi3bK3fcqytvugkKCmDMmCqni4th8mQoLdV/N240\ny/VTidfiKTu3chcF7GAkd1a7Zn6XBGLNHS4scF6LmqwFbVkt8JcsgRkcxdFMh4gNIbK9cse1uvLT\nT+HNN+Hmm6FFC8/8zXL91NO+fXQtfymdmMQlXM5jtGdFtetNm1b+byZ2a8AXX0DHjlx0Q7PULmjz\nY/dJ1ZFqG/7gwSJXqn+IgOzHSu8FKllITJtal5eLHHaYlLVoI0PP3B51AtyE3E0dThE0vY7W/CDb\nKJBnOafatfx8nZ/XvriGqjiGq95nH5FzzqlyvSYL2jCTttEpKRHpnT9XBOQcnjVCxwHfldFavHNx\nvad9C4GYOhRDXPhZfRt5BIMid6tbREB6MMdx8VZ+vvO9SumFXabtaJw6xv3z1+gPEycm7DlG4Ech\nLMj2ablHfqGhPMrlRujEy/btUtaqSOaoQx1j3Xt1nma5fnKJ7o5Z/ejfX+S0Y3+Rn2kq73JSTPfa\nf/Nsb0OPP647z8h380e0d+BfDy9JWH03At+DyF73fQbI1xwkrVsboeOF605K48aJgPTjY8fGbzbJ\nSB9+3DGd3AEPPVTkOsaLgPTlk7iFfra2pVGjnIU9iNzDTbKTXMljR8I6RiPwXXAa4o7kdilHSeuC\nzVlXQf1uTO1msx13vV5k9WnDU6MKEUPqGTy40mTm5wjPXQ0eLFKotstq2sjnFAuEYhb42drRRzOj\nTaWffMFhCe0YjcB3wWmI24+pIiAnq3ezqoL6nXjzqsAP51wtoUBALu6zIKoQMaSeWG34YeETvu8i\nJomA/IE3KtJ42fBNRx9l0SJ7ZAv15e9cmdCO0a/Azzq3TKcFRbPozR6C9JEZWeOK6bQdXnjR2R13\nwIABlW52bjsp7c8yLt3zCC83uoRn53RzfZbZJCN9uMWyz8mB3Fx3d8DwfS8XnM8iOnM3t5Cjyiks\n1EstbrxRp/ciW91rvRYtdmERDfiNL2wxJFPpAp51At/J/3s7e/EVPfkd/8uaCuq1HZ6I3g8jHCjN\nrQLfwwh2ksc1v4xl507nvPLzzSYZ6cZppfP//geffVb1XOQmNGPHwocf5/DeUXfRnYU8dMSzTJ0K\nAwfqOtG2LbRsCUo5PzdbO3qvNSZHoBdczaJ3xbmUdox+hgGpOtJlwweR+7leyqgnMz/ZnvQyZAJ+\nJ/MKC7WbXaQduDclIiCjGe15fzyxvQ0ZRigkcthhIkVFMvaWsmpmwJwcHRfJuNdqvMxoj/IX2Uyj\nKt5sqbThZ52G7zbE/TzvOPLYRe/y6mEW6iLNmvlLV1amq2W9evazwniG8xOtGM9wz/tnzzYrMWs9\nSsG4cbB6NVvvf6yaGXDPHt2Gjj/ebFkJzjImzFHMoIQ+hEWvUnoFdMrw0yuk6kiHH37Y/3vWlC0S\nCgbltW631vkwydF2RYo8WrXSWlz4c9iP+FL+GZPHRjZrfXWBb1r2l/XsLfXZYjxyfFBSoke44ZFP\nY3SEzFu4M+FtA+OlExujRomUBPrIDPrUaQEVq9eGUlX9iXPYJYvoJAvoKkF2+84nkcNXQ3o4r/MX\nIriHT85Gj5xo2D12TuZtEZBjmJbwtuFX4GedSceJsMfK1FA/juAL6rO1zoZJ9pqsdUIEyssrP1/K\n43RmCTdxL+XkxPz8bI9CWpvZ0/NwXrHCJ+/NhirXstUjJxp2h4ejmc4ucvmSwx3TpqJtGIFPpRD8\nmOPIodyKnqmpawLKy2UsGg3YwhjG8CnH8DanuqYLBt3zyPYopLWZYcPgzvw7KWAHt3JXlWvZ6pET\nDbvHztFMZza92IGzP2sq2kbSBb5S6iSl1GKl1DKl1M3Jfl48hIXg5xzJTupxHB9XXKtrAsrLZUwp\n7Wbndv0G7qcFGxjOeMDZFy8QgH33NaGP6yLFxXDajZ15Jngh/8ej7McqlNIdfKNGWjGqS6PhRBAO\nM55HGUfwBdM52jVtKtpGsjcxDwIPAwOBbsAQpZT7Cp00ERaCZRTwOUdWEfh1TUB5xbkvKICDD3Ye\nAezDD1zPA0xmMKUuQ1LQed92Wxyx9A21grFj4ZDXRqMCinvzxxAIaKXop59i3OQ+Swh77ByVN5s8\ndnkK/FS0jWRr+EcAy0RkhYjsAl4ATkvyM2PGLgQ/oR89+YombALqnoByc0sNr7J024nqdkaRwx5u\n5a6oKzUvucT7GWYRVu3m0D+04efBf+WssqfpUj6/or6E573uvVev1DbuuJqxY+HJi7WZeNvBR9K9\nu5YraWkbfmZ24z2AM4FJts9DgX9EpLkMKAVKi4qK4puiTgDhuDJ9lF5QNFi9UCe9dMK4hSV2CrZ1\nIN/IHgLyANdK06aV6aOFNjahj+sul5z+s/xCQ3md04w7rh9OPVWkc+eKj4luG2SCW6YfgW8/0umW\nKaJf+pBBe2RzsJlMa3d+VgooJ7fNdxgom2gs+xZszMp3YqjOoYeK3MKdIiDFfG7ccb0oLxdp0kTk\n4ouT9gi/Aj/ZJp0fgLa2z22scxlJcTE8/2KQxoMG0HfH+xQfEUp3kVJOpMnnOKZyMu8xPvcWLr6h\nqTHHGAA9r/V3NYy1tGQcN0PEntCR1DVvt0g89/f9+mvYvBmOPTZdxavET68Q7wHkACuA9kA94Gug\nu1v6dGv4FTz9tFZLSkt9x4uva4RHO4sKesj6wiKZOW2Ha7psfD/ZTngkeAV6T+gTeS+qlp8pzTvR\nRA0zPl5vJCM//JC0MpAJJh1dDk4GlgDLgVu90maMwF+3TgRkyrF3ZPdGzf/+t/7ikyc7XjYbWWc3\no0aJNCrYKctpL3Po4bi9pX3Fdl3cE8Fr5XqFGWvgQJEuXZJajowR+LEcGSPwRWRr18Pk88CR2WuP\n3LJFpGVLkT59dLTECHxVdEOdp6RE5O/FegP7s5nsKvCDwbpZJzw3OwmInDtol8hee4lccUVSy+FX\n4JuVti58EBjIEaGZFe6ZdmqbPdLTvujGPffAunUwYYJjwHOvEA217f0Y4qe4GK6aMYQlBQdzJyPJ\nYbdjuubN66Y7rtfK9VAI6n39JWzbphcoZABG4Lvw390DCRJiAB9Wu1abVt+OHq3r2osvwpw5PhfH\nrFoFDz4I550HvXs7JolW0WvL+zEkgECANw6/m44s52L+Xe2yUtCvXxrKlQK8Vq4HAjCw3lT9Avr2\nTW3BXDAC34VdPY5gPc05jTerXastq2+9tjH0DAp30036S95zj2ve0Sp6bXg/hsRxzLiTmRE4mtGM\npYDtVa4VFNStxYt2vFau5+fDgNyPoUcP/xtQJBkj8F24+togbwf/yKm8TR5lVa7VltW3cZldZszQ\nw4Abb4Q2bVzzjlbRa8P7MSSO4j6KRX++h9asZZh6CKi6ghTiMCvWArxWrt88bAeN5n+eMeYcwEza\nevH0Oe+JgJyq3q6VXijRtjGs9rrLy/VWdvvsI/Lbb1HzD3vpmK3tDGE2HXWqbM1tLMcesqnCTTcb\nvLkcV85+9JH+wu++m/TnY7x0EsDOnbJ7r4bycfuLamV4AKcwCXYPgrCbXLiyjtxPrz9YOuppz3zt\nvvf9++s9b2vj+zEkga+/1pL9pptEJMu9ua67TqRePV/KU00xAj9RnHOOSLNmIrt3p7skMeOnsYW1\nr/pslTXsI7M4XPYqKHfVvrJBWzNUJebFdeedJ5KfL7JmTVS3xTq9LWKXLiIDBqTkUUbgJ4pXXtGv\n6ZNP0l2SuPAyu9g7hHu4SQSkDzNcta+s1taylLg6+OXLRXJzRf7yl9jNinWFlSv1F5wwISWPMwK/\nhoS1miMP+U12BvPlpzOvSneRHPGjfXlFxlRKpDPfyi5y5Aku8NS+slpby0Jq1MFfdZVIMCjXnLLE\nl1mxzvHII/pLLlqUkscZgV8DIrWa1zhdfqKVjB65J91Fq0JNzSta+wrJh/SXzTSS5qzz1L6yVlvL\nUmrUwa9dK7LXXrKh/9nZOSo89VSRDh0cV6knA78C37hlRuDku/4s59KKtcy+/+OMcSeL28feRqdO\ncJZ6lROYwkjuZAMtqlxfsaKqC53xvc8uarS4rmVLuO469p7yIg+cOye7NsMpK4OPP4aBAx1XqacV\nP71Cqo5M0PCdtJo8dshmGsl/+HPGmC0SYV6ZNXWrfK/ayFccIkF2u+blZPPPKm0tS/Hr5eXKL79o\nh4cTT/S14Uedibz6/vv6Jb39dsoeiTHpxIeb2eJfXCJbqC9H9Ui+i5UfEmJeuflmEZDj8qa7Nmw3\nrx7je1/3SUgH/8AD+oaPP/ZMVqe8vy6/XAdM2749ZY80Aj9O3LSaY5gmAvJQn+czQhOpsfa1aJH2\npDj//Irv06SJu8C3jxrM1oXZQ407+B07RNq0Eend29WeXadGjuXlIq1bi5xxRkofawR+nJSUaBfi\nyIqnKJfvaCsf5Z4kwaA/TSSZHUONGkl5uUjfviKNG+vJNQszKWtwoqYd/LJb9L4K13V43fH+OuX9\n9fnnuuDPPpvSxxqBHyejRonk5DhXvtsZKeUoKWKVq5ANN45WrcR3x1CTssalff3rX/qGSZOqnK7x\nqMFgiGDUKJEGBbtlIV3kWzpLPbWrWh2tU4rGDTdoAbJ5c0ofawR+HHhpzSBSxCopR8lYbnMUiN27\ne99fkyGq22ghZu3rhx9EGjUS6dev2hC7Tg2tDWnHXp9O5S0RkKv5W7X65KVogEhRUS2pe6GQSMeO\nKVtdayftAh8Yg96wfK51nBztnnQLfK+hZfh4m5NlDfs4erVEm/iMd4ia0AmtP/1J26yWLPF8lpmU\nNdSUqu0pJO8zQDbTSPZmfbU5IS9FSalaUgfnzdMFfvTRlD86UwT+8FjuSbfAjza0BJHf86YIyGm8\nHjVtIoaoCdW6X3tNBOS5Q8bFtTLXYIiFyPbUhYWym6A8yl+qtYNIpaZWjjJHjtQa0k8/pfzRRuDH\nQbShJYgE2S3fs6+8z4C4hH2stvBoow7fw91Nm+TX+q3la3WI5LLLaO+GpOPUniYwTMpR0lN9Va0d\nlJTo+pzI0XHKCIX0ytoTTkjL4zNF4H8HfAM8ATRxSXcZUAqUFhUVJfetRCHa0DJ83MKdIiDdmRez\nwPerpYS17IICf3lGE9obThgiu8iRnsyufZqToVbi1J4as0k20Ew+CxwjJZ9Xd9OsbRO44Xb6584z\nRUCWjXwyLeVIicAHpgDzHY7TgJZAEL2r1l3AE9HyS7eGL+Jsww4Gq1a6pvwsv1FYJdhYtCMWO6Sf\n4W1MQvuFF0RARnJH7dOcDLUap/Z0Ve5j+sOLL1ZLn4meYm4OE/Z2+jeulh3kSauCX9IyYk67hl/l\nIdAOmB8tXSYIfJHqNuwTTqgufB/iKtlJruzDGk8hHwzqdRh+l5M//ri/UYZvob1mjUiTJvJNYW/X\n8AmZqDkZ6g7V5oSm75HfOh4sGwrbypE9tlXzOsskTzE3h4mLLqosZ5Dd8hMt5SXOTFs50y7wgda2\n/68FXoh2T6YI/EichpntWS57CMg4bnSsnA0bRp/wdKpMkaOJWI5qry8UEjnxRJGCArnm5MUZpzkZ\nspNRo0ROyPtUBOQuRlSbS8oUTzGvzsfelgagY+eczqtpGzFngsB/Bphn2fDfsncAbkemCny3YeZk\nzpYt1Je9WV/tWn6+dwX1O18Qi4ZfTWhPnKgvPvxwxmlOhuzEXg+f5HzZRY50Y37FiPiEE6ouYEyn\np5gfN20QeYkzZQPNpB5laRsxp13gx3NkqsB3E5ZdWCh7CMj9XB+zIPVbmexH2ETk61mzZulYOb//\nfcUCq0zRnAzZgZPt217vm7FBfqapfMbRoiivUpczoU76cdNuzjrZSa6M57q0jpiNwE8wkcIyfDzF\nn2U7+dKaHxy1brehnZ/K5CTUw7ZDT6G9caPIfvvpY+PGKs/NBM3JUPdxs323bl21Tl+IjrNzIf9O\n26jTbVLWj5v2cO4TAenCwopzwaBI586pbV9G4CcBp6iS7Vkuu8iRh/k/xwrh9pWiVaZg0F2oewrt\nUEg2Hv172R3IlaGdZxmhbkg5XuZDe3wpfYTkU34nP9O0imk0VXbwyI7J7mjRv79zIEV72RfRSf7H\nUY7XUzmCNgI/iURW6L9zpewhIAczt9oP7ja0i2ZTf/zx+DTxKf30GoGrmWjMNoa04GWudDJLdmWB\n7CRXnuLPvpSlROEnpENOjraMOilnfflEBOTPPOU5CsjLS77SZQR+krFrBk3YKOvZ2+rpQ76HpYm2\nqS++62URkGc4t0o5Uj1ENmQ30cyVrVpVF7S3M1IE5GTejqos1QS7+aaoyN88Wl6enkzu1Utr/fXq\n6fOvc5psoJkUsC1qHt27J/672DECPwXYQyFfrLQtcij/iUlwJ8ym/uWXUhYskBn0kTx2OA4vzeIq\nQyrws3iqpEQLz3C6epTJ1xwkP9BaGrMpKQpKPAsaI9tOOHz6/iyVcpTczkjfeSRT4TICP8WUzCiX\nxc2KZXOwmVx+2o+p1aZXrhTZZx/5sV6RtGCta6Wrxa/XUIuIxQXYPsrtyWzZRY48Fxxabb6qppsI\n1dQNulevqnk8xFVSRj1pyU8VAj1aHslUuIzATwfffitSUCCb+gyUwWeHYq6kcVXuH38U2X9/kcaN\nZfhJ81wrnlJ6CFvrN4g21ApiMVfaR7kvdx+tb3jjjYo87PU4GNSearESjxt05KgknEdjNslvFMqT\nnKZ9cXQAABflSURBVF+lffnpNJKFEfhp4u2B/xABuYYJUSu6nbhi3v/8s8iBB+oNk0tKfE1CxVIm\ng6EmxGWu3LlTpEcP2dm0pRQVVF/QGD68hL6T4hSPG7S9k3n88co8RjFGBKo5aXh59CTbN98I/DRQ\nUiJSWBCS1/ij7CYoxzHFdShb7b5YV8H+8IOeCcrLE5kypeJ0pGZVq+OLG7KTb76RnYE8eYeBVRZk\nRQphp7rrpjh17+7P7OKmKIXzaKo2yS80lFc5vZpAP+EEd6Gf7LZmBH4aCA/56rNF5tFdfqZpxbJx\nr0nTaMPN/v0jbli+XMferl9fZOpUEamq1fTvX+lV4OWJYCZyDZnKPW0fFgG5nvtd20Vk3fVSnPLy\novnURz/y8kTuCd4qAnIQXzsK9HStZjcCPw3Yh43tWS4/0Fp+pJV0ZImAuw0v2nAzELBVmE8+EWne\nXKRpUx0+QbzNQbUtvrjBICIy+OyQvMwZsoscOYKZvuqul+Lkd89pr6O52iC/BerLy8FBngI9HavZ\njcBPA5HuaF1ZIOvZW36ipfRWs1xteH6WcNcv2COrrrpfj2W7dBFZtEhEopuDTjgh8+KLG7IXv44J\nJSUiTQObZSX7yQraSWM2VTOzRNZdP8pNpDto5OE10g4vsJz7/IKMC09iBH4acBK+XVgoy2kv28mX\nZSMerwhkFu0++9GJRfI/jtIf/vQnkV9/rbg3mlbTv7+JkmnIDGJ1TLjoIpEjmCk7yZUP6V9lPwen\nuhvL5ilOppe8PB3W3On+A/lG9hCQDw64InkvqAYYgZ8mnCpSUcF6Wd7uOH1iwACRr75yvC9ScLdh\ntTzGZbKboGyisdzW7ulqHYYfrcZEyTSkGy+lxh4WOZKLLqpc1Pg3rvbcOS5W5we76aV7dy8bf0im\ncJxspIl8+f7PyXlBNcQI/DTiaMMrLxf5xz9EGjWqFPyTJmnTzO7dIiJy8vFl0olFcgn/kv9yiuwh\nIDvJlYe4SlqptY7mFy+txu57b5/IzZRhqCF78OMH7yXI3+l0jQjIPw9/POHhSqKNsC+yInr+9+RH\nav4ikoQR+JnK5s0iY8eKtG9fVTLn50vI1iJW00bu4FYpYlVFYwgHVLNvh9i/v3dDMr73hkzArx+8\nq5lx9269e1swKPLf/3o+y8+kqd+YOvvyvfyqGsqvPftqpS1D8SvwlU6bGRx22GFSWlqa7mKkBhGY\nNw+++gqWL4edO6GwkNe+asedHxzB1zu7EBJFIAD5+XDYYVBaCjt26FvjpbAQpk6F4uLEfRWDIRpD\nhsBLL0Eo5J0uEIBBg2DyZIeLW7fCccfB/PnwwQdwzDFxlWX0aBg/PnpbClDOhwzgyMBMCpbOgw4d\n4npeKlBKzRaRw6Im9NMruB3AWcACIAQcFnFtBLAMWAyc6Ce/rNDwfRCpocS7sbnb5JXxvTekmlhi\n2XiKgQ0bZPt+XWRbbkM5v9PnUT19Ij2CYinHndwiAvJY7yeS8k4SCakw6QBdgc7ANLvAB7oBXwN5\nQHtgORCMlp8R+M7UJA5IzA3KYEgSTrFxnBQSL1fhUaNEDshfLUvoKFvZS/qrKY6mSq8Vt37a0hno\nUONPBi+uFfNdfgV+oCbDCBH5VkQWO1w6DXhBRHaKyEpL0z+iJs+qzcycqYe0vXrpvzNnxnb/kiU1\nM+PYCQSgU6fE5GUwxMLYsdqc2L+/rodO5OfD1Vc7X5s5U5tilpa15Rg+YyXt+a+cwsDtrzB+fGW7\nCqfbvr2y3YRC+vPChdHbUj8+5jnOpSRwJN/f+Pe6Zf700ytEO6iu4f8DOM/2+d/AmS73XgaUAqVF\nRUVJ7APTQ1xB0SLwszDL72F87w2ZQDzeNJEj3ab8LJ9TLAJyByNlyNnljun8HkqJnNv8A9kW2EtW\nN+ouX3ywKUVvo+aQKJMOMAWY73CcZksTt8C3H3XNpBNXULQY84mlMhsvHUMmEWsIAidPn3qUySQu\nEgH5okE/kZUrfYUqqX4+JJfX+7eUB3NEDj5Yhx13KW8mhhhPmMD3lUl1gT8CGGH7/AHQJ1o+dU3g\nR9M0ior8V5pIjcguyMOVOHL/TfuGzJlWQQ2GWHEf6YbkEjVJtuU0kD0Fe8ldje913XYwENARxe1t\nqYVaLy8Gh+gPxx0n8ssv1Z6diJF6Mkm3wO9O1UnbFWThpK0f3+NYKo2T906khpSOwE0GQyqINmIe\nPug7eTdwigjID7SWEdwl+/K948i6pETkqt+vkn+1vk225TSQUDAocuedInv2xPzcTGhjKRH4wOnA\nGmAnsA74wHbtVrR3zmJgoJ/86prA92t7z5RKYzBkOm62/4suqhTKR/OZTOG4igY2j+7yFH+W+3JG\nyGdH3Sxy8cUiPXpUNsDTT9e71bkQLV5VJrg5+xX4ZuFVEpk5E44/XnsHeOG52MRgyCJmzoSJE7Vn\nWqdOMGxY9UWC4TRLl8IBB+g0EyfCiy9qMRymA8sZxEscyzQOCi6kFWsJKKBJE+jRQy/cOvdcaN/e\ns0y9esGcOd7X0y22UrLwKtFHXdPwRfz5HoO7b3wmTxQZDIkkVju5vW00aeKjfTlEqvVDLFE40wUm\nlk7mUFKiJ2i9vAacKk2mTxQZDIkiVjt5ZNuI5qEWDiIYj9JkbPhG4MdMPKFbM72SGQyJwo+dPKzR\nd+6svc+iCfpIoV8TpSnTQ4z7Ffg1Wmlr8E9xMQwfroOXhVcZBgL68/Dh1e2UEyfq4E5OlJXp6wZD\nXcFrNXkoBNOm6fmwF1+ExYuhvDx6nkpV/h/OO7zi1r4y1w/hVcKDBmmb/aBB+vPYsf7zyARy0l2A\nbGLsWBg4sPqEk9PS7WgNYOnS5JbVYEglnTrB3LnO0TSVgg0b/Al50HOyHTroe77/3rkdhZWmWMIm\nFBfX/iizRuCnGL+VxqsBmHg4hrrGsGHw1lvOHm2BQPSwyva0J50Ezz+vNfHVq53TZavSZEw6KcQe\nRO2EE2DAAPeAasOG6UBSTngFmDIYaiNeJs8WLdxHu5HY20anTu5B2rJWafJj6E/VUZcnbb28Ctwm\ngDJ9oshgSDROK8X9LGB0ahvZ5PiAWXiVOfhdgOW0G5XTIpPabkc0GGLBq/0Eg7pd9Ojh3DbCu1uV\nlWkzTngHueHDa9+Eqxd+F14ZgZ8ChgypvgrQCbPi1mBwpiaCuyZKk5+Vv5mAEfgZRLSl2XaaNIF3\n383MSmUwpJNUj3Yj977N5NGBEfgZhN8NnMOEffMzrVIZDNmClxnJyfSabvwKfOOlkwK8PG6ciGdh\niMFgcCfWbUbr6sJH44efAsIuZ3YbZDTiWRhSG9i9ezdr1qyhrKws3UXJGvLz82nTpg25ubnpLkpa\niDTNzJ2rff69RtF1deGjEfgpInKVbZMmMGsWbN3qnL42Vyov1qxZQ4MGDWjXrh3KvvbdkBREhI0b\nN7JmzRraRwkDXBexb2gexh5eYeBAZ6Wqri58NCadFFJcrD1wSkvho4/glFOyb2FIWVkZzZo1M8I+\nRSilaNasWZ0dUUUz1cRrmqmrCx+NwE8jdbVSRcMI+9RSV9/36NGVAdXmzIEXXoAjj4QDD6wU/PGa\nZmINdlhbqJHAV0qdpZRaoJQKKaUOs51vp5TaoZSaax2P1byodY+6WqkSSayTbYbswG6qsQt0EViw\nAPr10x1CTcIr1JUImVXwsxzX7QC6Ap2pvol5O2B+rPnV5dAKXmTTxuMLFy70nTZZG8AEAgE55JBD\npFu3bnLwwQfL+PHjpby83POelStXynPPPVezB/vg4osvlgULFnimef3116OmiSSW914b8Iqfbw+f\n8Pjj2RFegVTEwxeRb0VkcY16HEMV2/7kyUazB2cNLt5Y5pEUFBQwd+5cFixYwEcffcR7773H2Chq\n26pVq3j++efjf6hPJk2aRLdu3TzTvPHGGyxcuDDpZclkvEw1YcrKtEZuRtE2/PQK0Q6cNfxtwFzg\nU+B3HvdeBpQCpUVFRcntBg1px6+m6WcHpHjZa6+9qnxevny5NG3aVEKhkKxcuVKOPvpo6dmzp/Ts\n2VNmzJghIiK9e/eWhg0byiGHHCIPPvigazo7K1eulM6dO8s555wjXbp0kTPOOEO2bdsmIiJTpkyR\nHj16yIEHHigXXnihlJWViYhI37595csvv6wo5y233CIHH3yw9O7dW9auXSszZsyQJk2aSLt27eSQ\nQw6RZcuWycSJE6Vr165y0EEHydlnn+34neuihh8toFrFXrZS90fRJGqLQ2AKMN/hOM2WJlLg5wHN\nrP97Ad8DDaM9K1tNOtmEX8Fz6KH+GnI8RAp8EZFGjRrJ2rVrZdu2bbJjxw4REVmyZImE6+Qnn3wi\np5xySkV6t3R2Vq5cKYBMnz5dREQuvPBCuf/++2XHjh3Spk0bWbx4sYiIDB06VCZMmCAiVQU+IG+9\n9ZaIiNxwww1yxx13iIjI+eefLy+//HLFc1q3bl3RYWzevNnxO9c1ge8VCTPaXtF1Eb8CP6pJR0T6\ni8iBDsebHvfsFJGN1v+zgeVAHXQyNCSLdMUy3717N5deeikHHXQQZ511lqvpxG+6tm3bctRRRwFw\n3nnnMX36dBYvXkz79u3pZH2J888/n88++6zavfXq1ePUU08FoFevXqxatcrxGQcffDDnnnsuzz77\nLDk52bG0JuzwkJfnnqYue7rFS1LcMpVSzZVSQev/DsABwIpkPMtQN0mly+qKFSsIBoO0aNGCCRMm\n0LJlS77++mtKS0vZtWuX4z1+00W6RMbiIpmbm1uRPhgMsmfPHsd077zzDldeeSVz5szh8MMPd01X\n1xg7Vu912717VeUgq230UaipW+bpSqk1QB/gHaXUB9alY4BvlFJzgVeAy0VkU82KasgmUuWyumHD\nBi6//HKuuuoqlFL8+uuvtG7dmkAgwDPPPEO5tZFqgwYN2GpbFu2WLpLVq1dTUlICwPPPP8/RRx9N\n586dWbVqFcuWLQPgmWeeoW/fvr7LbC9LKBTi+++/p1+/ftx77738+uuv/Pbbb3G9i9pIcTHMnw8z\nZsDgwXXIfTJJ1Gj8JyKvA687nH8VeLUmeRsMsWz6Hgs7duygR48e7N69m5ycHIYOHcp1110HwBVX\nXMEZZ5zB008/zUknncRee+31/+3df2hV5xnA8e9DanuhKbpO1+pSZ5oVMR0anYqROgrWzQapc8QR\n/5miWCLbnIExUoRaqAXrmIrMbritSTeKprjVhWKZRjIGit2PGE008Ue2iFHrjwztVNC5PfvjnLjj\n5d7cc73nR5LzfCB4PO97733uc3KenHPuue8LOJdNioqKmDZtGitWrMjaL93kyZPZsWMHK1eupLy8\nnDVr1pBKpWhoaGDp0qXcu3ePWbNmUVtb6zv+mpoaVq9ezfbt29m9ezerVq3ixo0bqCpr165lzJgx\nhSVoGBoJE4xHwYZHNpHq6upiypQpcYcRid7eXhYtWkRnZ2fcoSQq70lkwyMbY4x5gBV8Y0IyadKk\nIXF0b8wAK/jGGJMQVvCNMSYhrOAbY0xCWME3xpiEsIJvEqe4uBiAixcvUl1dHchzdnd3U1FRwfTp\n0+np6WHu3LlAdKNsGuOHFXyTWBMmTGDPnj2BPNfevXuprq7m6NGjlJWVcfjwYcAKvhlakjHSkhma\n1q1zZooOUkUFbNvmq6v3i1GNjY00Nzdz+/Ztenp6WLJkCZs3bwZg//79bNiwgTt37lBWVkZDQ8P9\nswSAffv2sW3bNoqKijh48CCtra0UFxdz8+ZN6uvr6erqoqKiguXLl1NXVxfs+zUmD3aEb4yrvb2d\npqYmOjo6aGpq4vz581y7do2NGzfS0tJCW1sbM2fOZMuWLQ88rqqqitraWurq6mhtbX2gbdOmTcyb\nN4/29nYr9iZ2doRv4uPzSDwq8+fPZ/To0QCUl5dz7tw5rl+/zsmTJ+8PcXz37l0qKyvjDNOYh2YF\n3xjXY57B1QeGI1ZVFixYwK5du2KMzJhg2CUdYwYxZ84cDh06dH8o41u3bnH69Gnfj08fVtmYOFnB\nN2YQ48aNo7GxkWXLljF16lQqKyvp7u72/XjvsMpbt24NMVJjcrPhkU2kbJjeeFjeR7ZIhkcWkR+L\nSLeIHBeRD0VkjKftNRE5KyKnROQbhbyOMcaYwhV6SecA8BVVnQqcBl4DEJFyoAZ4HlgIvDMwx60x\nxph4FFTwVXW/qg7MmHwEKHGXFwO7VfWOqv4DOAvMLuS1zMgxlC4jJoHl2wwI8kPblcDH7vIXgfOe\ntj53nUm4VCpFf3+/FaGIqCr9/f2kUqm4QzFDQM778EWkBXg6Q9N6Vf2922c9cA94P98ARORV4FWA\niRMn5vtwM8yUlJTQ19fH1atX4w4lMVKpFCUlJbk7mhEvZ8FX1ZcGaxeRFcAiYL7+/7DtAvCMp1uJ\nuy7T8+8EdoJzl07ukM1wNmrUKEpLS+MOw5hEKvQunYXAj4BXVPW2p6kZqBGRx0SkFHgO+HMhr2WM\nMaYwhQ6t8FPgMeCAiAAcUdVaVT0hIh8AJ3Eu9XxXVf9T4GsZY4wpQEEFX1W/PEjbW8BbhTy/McaY\n4Aypb9qKyFXgXAFPMRa4FlA4QbK48mNx5cfiys9IjOtLqjouV6chVfALJSJ/9fP14qhZXPmxuPJj\nceUnyXHZ4GnGGJMQVvCNMSYhRlrB3xl3AFlYXPmxuPJjceUnsXGNqGv4xhhjshtpR/jGGGOysIJv\njDEJMawKvogsFZETIvJfEZmZ1pZzwhUReVJEDojIGfffz4UUZ5OItLs/vSLSnqVfr4h0uP1Cn+pL\nRN4QkQue2Kqy9Fvo5vGsiNRHEFfWiXTS+oWer1zvXRzb3fbjIjIjjDgyvO4zItIqIifdfeAHGfq8\nKCI3PNv39YhiG3S7xJEzEZnsyUO7iHwmIuvS+kSSLxF5V0SuiEinZ52vWhT4vqiqw+YHmAJMBv4I\nzPSsLweO4QzzUAr0AEUZHr8ZqHeX64G3I4j5J8DrWdp6gbER5u8N4Ic5+hS5+XsWeNTNa3nIcX0d\neMRdfjvbdgk7X37eO1CFMwy4AHOATyLaduOBGe7yEzgTDqXH9iLwUVS/T363S1w5S9uun+J8OSny\nfAFfA2YAnZ51OWtRGPvisDrCV9UuVT2VocnvhCuLgffc5feAb4YTqUOcAYa+DewK83UCNhs4q6p/\nV9W7wG6cvIVGs0+kEzU/730x8Gt1HAHGiMj4sANT1Uuq2uYu/wvoYvjMMRFLzjzmAz2qWsi3+B+a\nqv4J+Gfaaj+1KPB9cVgV/EH4nXDlKVW95C5/CjwVclzzgMuqeiZLuwItIvI3d16AKHzfPa1+N8tp\nZNyT13gn0kkXdr78vPe484OITAKmA59kaJ7rbt+PReT5iELKtV3izlkN2Q+64sgX+KtFgeet0NEy\nAyc+JlwJgqqqiDz0Pak+41zG4Ef3L6jqBRH5As6Io93u0cBDGywu4GfAmzg76Js4l5tWFvJ6QcSl\n/ifSCTxfw42IFAO/Bdap6mdpzW3ARFW96X4+sxdnaPKwDdntIiKPAq/gzredJq58PaDQWpSPIVfw\nNceEK1n4nXDlsoiMV9VL7inllYeJEXxNDPMI8C3gq4M8xwX33ysi8iHOKVxBO4rf/InIL4CPMjT5\nnrwmyLgk80Q66c8ReL7S+HnvoeTHDxEZhVPs31fV36W3e/8AqOo+EXlHRMaqaqgDhfnYLrHlDHgZ\naFPVy+kNceXL5acWBZ63kXJJx++EK83Acnd5ORDYGUMGLwHdqtqXqVFEHheRJwaWcT647MzUNyhp\n102XZHm9vwDPiUipe3RUg5O3MOPKNpGOt08U+fLz3puB77h3nswBbnhOzUPjfh70K6BLVbdk6fO0\n2w8RmY2zf/eHHJef7RJLzlxZz7LjyJeHn1oU/L4Y9ifUQf7gFKk+4A5wGfiDp209zifap4CXPet/\niXtHD/B54CBwBmgBngwx1kagNm3dBGCfu/wszqfux4ATOJc2ws7fb4AO4Lj7izM+PS73/1U4d4H0\nRBTXWZxrle3uz8/jylem9w7UDmxLnDtNdrjtHXjuFgs5Ry/gXIo77slTVVps33Nzcwznw++5EcSV\ncbsMkZw9jlPAR3vWRZ4vnD84l4B/u/VrVbZaFPa+aEMrGGNMQoyUSzrGGGNysIJvjDEJYQXfGGMS\nwgq+McYkhBV8Y4xJCCv4xhiTEFbwjTEmIf4HTS6c5dUVe58AAAAASUVORK5CYII=\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)" ] } ], "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.1" } }, "nbformat": 4, "nbformat_minor": 1 }