{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# EE16A Discussion 13A" ] }, { "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+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJztnX+UFdW157/7Xppu20dMgyAdoQdM\nhGCLgN2JncTRIA2KuuJL8tQm0eHJMwyJTsisMcbEJ4S4XBOeGpczJnEJSZ5jIqjvJZGVyDKAuhx5\nkElD0PBDfig8xYfSIqIRWht6zx+nyi4uVXVO/bx1q/ZnrVp9761zq07XPbVrn32+Zx9iZgiCIAj5\np1TtCgiCIAjpIAZfEAShIIjBFwRBKAhi8AVBEAqCGHxBEISCIAZfEAShIIjBFwRBKAhi8AVBEAqC\nGHxBEISCMKjaFXBy6qmn8pgxY6pdDaGA/Md/AG+8AfT3D3xWKgGnnQZ87GPVq5cgmLBhw4Y3mXm4\nrlymDP6YMWPQ3d1d7WoIBWP9emDatOONPaDeHzoE/PrXQEdHdeomCCYQ0b+blIsc0iGi0UT0NBFt\nJaItRDTf+nwoEa0iop3W36ao5zJh/Xpg1iygrU39Xb8+jbMKtYaznVx9NXDkiHu53l7g3nvTrZsg\nJEUcHv5RAP+DmTcS0RAAG4hoFYC/B7CGmX9IRLcAuAXAd2I4nycLFwJ33aVuXmZg0yZgxQrgppuA\nRYuSPLNQS1S2Ez/6+4GdO9OplyAkTWQPn5n3MfNG6/W7ALYBOB3AFQAetIo9COBvo57Lj6VLgTvu\nAA4fHriJ+/vV+7vuEk9fUKxfr9qDs534USoB48YNfFd6j0ItE6tKh4jGAJgC4I8ATmPmfdau1wGc\nFue5gIEbsLkZ+NrXgGPH3MtJt1ywufde7/CNGw0NwDe/qXoF06YBjzwCbNwIPPqoer9wYXJ1FYS4\nic3gE9HfAPhXAN9i5nec+1gl3Xf1p4hoLhF1E1F3T0+P8fnsG3D5cuD11/3LSrdcsNmxw9yzb2xU\n4UDgxF6B9B6FWiQWg09EdVDG/lfM/Gvr4zeIqNna3wxgv9t3mfkBZm5n5vbhw7WqIgDHd8tNcHbL\nhWIzbpxqD24QAS0tKmRz1VXAmjVq7MevVyC9R6GWiEOlQwB+BmAbM//IsWsFgNnW69kAHo96Lpuw\n3XJBmD9ftQc3TjpJhWy6u4FlywakmH69gv5+4MknJa4v1AZxePifA3AtgIuIaJO1XQrghwCmE9FO\nAJ3W+1gw7ZYDQLmsuuWioxYA1Q5uukmFa2xP3xm+cWsnfr0CADh4UOL6Qm1AWVrTtr29nU0mXs2a\npW6uyokylZTLwP33A9dfH1MFhdywfr3qKe7cCZx5pvL8vZwCe2KWaQixsVGFg8TJENKCiDYwc7uu\nXE3m0vHrlgMqFtvYCNx6qxh7wZ2ODhW2qQzfeJWt7BX4IXF9IavUpMF3uwGJlEff3KxmTtoDboIQ\nB4sWqTZ11VUqXt/kM2+8vx/485/Tq5sgmFKTIR2bIN1yQYgTXVixXFY9THE6hDQwDenUtMEXhGph\nEteXWL6QFrmO4QtCtbHDiuWydxmJ5QtRiTudhxh8QQjJokXAJz7hvV9meAtRSCKdhxh8QYjAlCne\nyh2Z4S2ExS3Jn53OY/Hi8J6+GHxBMMCra+0nEZYZ3kJY/LIJvP9+eLm5GHxB0ODXtfbS6BMBY8dW\nr85CbaPLJrBtWzgvXwy+IPjg17W2M2XaGv0JEwaMPrO6Kb1irqaDcZKDv5joQoH9/SEFAcycma2t\nrY0FIUtMn86szPeJW6nE3NWlyq1bx9zY6F6usVHtt1mwQH1GNHCcxkb1uRPTckL+WLdu4Hf32pzm\nEkA3G9jYmvfwxQMSkmLhQmD1au/9ThWOaQplkx5DkHJCPunoAM46y3t/WEFATRt8WYVISArb4PrF\nUZ03nWkKZdMF0yUHv7B0afyCgJo1+OIBCUlisuZCXR3w5pvKkL/5phqo9cJOofzKK/4PBrvHoHuA\niL4//3R0ADffHCyVt45B8VYxPUw8IJnSLphg52TasUN57PPnm6250N+vQj7MythHzVLi7DGMGwds\n2uSeq0f0/cVh0SJg5sz4cobVrMEXD0iIg4ULVY/wyBHVnjZtAn7zG6C+3v97pRLQ1zfw3tkWwxp/\nZzd9/nxgxQr3XD2i7y8WHR3xOa81G9LxW4VIPCDBBK+w4PvvA++84/29ctnboBMBo0frUyg7ceum\nh1mZSxB0xLWI+c+JaD8RbXZ89n0ieq1i2cPYkBmOgg6dgivo2sj2wjojRngbfGZg+HC1sMrFF+sX\nTCFS+n239Rsqc/A7F1YXhDDE5eH/M4BLXD6/h5knW9sTMZ0LgHhAgj8mCq4gayMDynNfswa48EKz\n3qVuZTZAnX/3bu/9QVbmEgQdsRh8Zn4WwFtxHCsI4gEJbixdCtxxh17BpVucvJLhw5XBNe1dOp0S\nP0RmKaRF0jH8G4noBSvk4xrRJKK5RNRNRN09PT2BTyAekOBk4UJg3jzg2DH3/YcPA7fdpl6beOA2\nTs89SO/Sdkp0SyKKyEBIgyQN/k8BfBzAZAD7ANztVoiZH2DmdmZuHz58eILVEfKOPQjrZextnnrK\nP/GZG3V1x48LBelddnT4x/NFZCCkRWKyTGZ+w35NREsA/C6pcwkCYD4I29+vHgwzZ56oc+7tBV58\n8cSHRn8/sHLl8d67Ti7n1PcPHQoMHqyOX0lQkYHbvAHp2QomJGbwiaiZmfdZb78IYLNfeUGISpBB\n2N5eFdo59dQBw3nffWrf1KknGvy+voGHhIlxrdT3l0pqq6tTx+7vV+8bGoKJDNzmDaxYoY4hY1eC\njlgWMSeiZQA+D+BUAG8AWGi9nwyAAewB8F8dDwBXZBFzAQjvwc6apdQ4brNT3SiVBnIP2sZ37Fhg\n61b3B0eppEI3y5bp6++1wHl9PXDBBcBbbwWfNel3XFkwvdiYLmIei4fPzLNcPv5ZHMcWikUUD9Zv\ndqobzgeDreLxMvZ2GZPBVb/QUl8fMGwY8Ic/mNXR9LiSTkQwoWZn2gr5I2pCvCCDsF6YZsf0wyRz\nZpjkfpJORIiKGHwhM8SRErhSPdPaqkI1TvmkX1ZLu4wbpoOrOn3/wYPh0njrjtvTI1liBX/E4AuZ\nIS4P1jk3Y/Nm4OmngYsuUlr4U05RqRG8jH6ppBaeiDKD20TfHyaNt+64r74q60EI/hTC4MuqWNnG\n/n38UgyUSkraGOZ3XLkS+Ld/A95+W3nX+/d7P1gaGoAlS6LN4E5qhq3zuG4PLGZZD0LQYLIOYlpb\nEmvayrqg2aby9/HaBg1ibmgI/jv6rTULJNsu1q1jbmoyX5c0yHFbWszW2hWKAYqypq0fsipWtnH7\nfSoplZSUkUh5xEF/R79xAWcq4yTyMAWdYWvaE+3oUPMHvJAB3NoklUiEyVMhrS1uD7+ry9tzFC+o\n+vj9PoDyjru6mKdPD/87nntu/B52EPx6GI2Naj9z8J5oV5cq43VNZs1K9v8S4iVqJALi4YuMLevo\nZsaOGKEGKru7w/+OSSyUE8QTM0m0FqYnKutB5IdUIxEmT4W0tiQ8fPGCsovf72PH1+vq/D103e9o\n6mGbEtYTW7dO/b9tbeqv87xhe6J2XexrKONTtUkckQgYevhVN/LOLW6Dv2QJc7kc380uxItuQNVk\nM/kd4zKMcT88bKKEnfweJEJtoPv9m5r0v2vhDb59k7t5jeIFZYcFC7wfyibG3v4dbcN37rnuhi8O\nw5jUmJD0RIuNrqdb2dbdKLTB9/PEymXl+QvZYfz4cAZ/zhz1fTdpJxFza2u8Hq+pJ657+NjY5caP\nl55okTHt6fq1hUIbfFHn1BYmHo7XDbBkif/N0tAQX2/OxBM3jfF7zT/w+p7pQ0SoTbwiEqa2q9AG\nv9pSPCEYYWP5pZKagKSbtOX0jKIYTl0M3+/hU1kHvx7oJz95fN1k8mAxiDJRz9Tg51KWmYQUT0gO\nL+niIE3y7v5+4LXX1K3gh53CYOFClWvmkUeAjRtV7vwguWd0Ess1a8ySv/lNBmMGJk8eWJ85qGRP\n0ojULqkshWnyVEhrSyOGLzHR7OI2sDp9enDP320bPz6+NuE1AGzaswzSAw0SnpSeQO0T1nahyCEd\nZtEo54U4pJu60E9c4zqmapsgqpwgA8Xi5OSDMLYrVYMP4OcA9gPY7PhsKIBVAHZaf5t0x4lbhy8a\n5XxQeQME3Rob9UqgOJqeqdE1KWe3Xb+YbuVDRIQK+SGo7Urb4F8A4NwKg/9PAG6xXt8CYLHuOElk\nyxTygfMGGDLEzNA7PaO0tO6m3plfOdMMos6HSBoPNCG7pB7SATCmwuBvB9BsvW4GsF13DDH4go2X\nmkY3UYtIhW8qPaM0Qx6m3plTh9/Sov5+6lP6nkzlQ0R3TWTyVv7JgsF/2/GanO8rvjcXQDeA7paW\nlmSvilATeA0+zpmjj+f7Ge8sjuuYevP2ZmcQNXmQ2Vt9vRr8Fg1/fsmUwbfeH9QdQzx8QadR9zOM\n5XK0JGZpE2ZAuqnpeMOtSzFtJ6AT5U6+MTX4GqVzJN4gomZm3kdEzVCDuoLgi59G/dgx/+9+4hP6\nBUw6OszWpU0Dv//Vi4MH1bZpE7BihVqjV/lT3vT1Dbx2avhnzszOtRDSIcmJVysAzLZezwbweILn\nMkImpWQfXY58L0ol4Nxz469PkoT9X4EBw71/v/eC7ID38YOupyvkg1gMPhEtA7AOwHgi2ktE/wDg\nhwCmE9FOAJ3W+6oRdZalEA+6h67fLGkioFx231eLi374/a+m9Pf7Xy+/71UuHCMOUQEwifuktcUZ\nw3eqPDo7VRKtNBQagjcmM0F1ahp74DZLA69hiWNSGcA8cqT7NWltNZeiyizd2gZFnmkbRPlAJJNS\n0iCILFKnpnEOvHZ21rYCxW9SWRDD7TYYHedEMCHbFNbgh/GaRo6MfFpBQ9CZoCZqmrx4pW4PsDCG\n2w0TKarfbyMOUbqEzeZaWIOvk6m5beWyeDFJE3fK6qJ5pVHmEOgenrrfRhyidIjiwBTW4Osar3gx\nyePmpcSd2qCIuWOSmkOgc5LEIUqeqA6MqcFPUodfFcaNUxrl/n7z7zCfqFgQwrFwodJ4Hzmirqut\nF+/qUkqaw4dP/E4YhY2fpNFNgZIHkppDMH8+8Nhj3vMc+vuVhFM0+8nhNyfDltDGcf1ztwDK/PnK\ngARBFkWJB7/FOpYvV0bfa/GQoI1ZFrmJj44OYMQI7/3iECVPWg5M7gy+16pE9fVAXZ37d4iUNl+0\nx9HQeSmHD6tVoa66Smm9r7pKvdfNjnXD78Fei5r8MMSpm7/wQm/dvjxAkyc1B8Yk7pPWloQO3xnv\n1OVVr1WVR1ZIey3hLCZDS4u4FUpFGwTPGmnF8Ktu5J1bGsnTnClpvVLKSgMPR1o5551kKRlaWsRh\nnN0G1ov8AM0CUa6/GHwNRVR5JI14iekQte369Q6K+ADNEmGvv6nBz51Kx5QiqjySxh4/uesuFbO3\n87w0NIQbmBXcidJ2nQPrzu84M2guWxZvfQVzks7mmrtBW1NE5ZEMixaZDcxKoq7wRGm7JvI/IceY\ndAPS2tIM6Uj4oXrkJSVCtYjSdqMMrIed9i8kDwxDOoX18L3km2F14YIZflr9u+4ST9+EKG03bO9A\n0ovnBJOnQlpbNZY4lEGqdJHB8vjwart+nniY3oH0hrMPRKUjZAWnAWpqCh9SEPS4pQYnUimWbcMc\nVP4nD+nsY2rwC6vSEdKhMreOHzJYHg03BQ6grvuWLcDUqcDNN6sB9Jkz1QDtzp3AmWeqmcteoSBR\ntOWHxA0+Ee0B8C6AYwCOMnN70ucUsoGXAfKiKCkRkkK3KHpv7/GLl5uOU/klJJSHdG2R1qDtVGae\nLMa+WOgMkI0MlseDyaLoYaSXkrcoGaohTS6sSkdIHp0BamiInkRNGMBkUfQwIRhRtMVPtVRPaRh8\nBvAHItpARHMrdxLRXCLqJqLunp6eFKojpMW4cd4ZGAGgrw+47z41s1OMRnRMUoOHDcGYTqgT9FRT\nmkys6wNGPQHR6cz8GhGNALAKwH9j5mfdyra3t3N3d3ei9RHSY/164PzzvRfWIAKuvlqm8sfJwoXA\n4sXA+++776+vBy64ADhwQBl+v8FaIRlmzVKevZfp7ewEVq0Kdkwi2mASMk/cw2fm16y/+wH8BsCn\nkz6nEB9R4oyysEb6LFoEPPMM0Np6fHinVAIGDVKe5OrVMnmqmuhCnU89ldxvkqjBJ6KTiWiI/RrA\nDACbkzynEB9xxBl1C2sMHSo5deKmowPYvBlYu1atMtbWBlx0kTL4fX0nhhEWLwZmzDj+N5BcR8mh\nG2vp708wtGMi1g+7ATgDwPPWtgXArX7lqzXxSnKEnEhcsyv9jjNoEHNDg+TUSQPdQuXOiVSDBjHX\n1cnvkhR+90TYCW2QmbZmSCIvd+KcXek2s7O+XhkVma6fDrqkaSab/C7xsWCB/gEcxByaGvxCyzIl\nkZc3cc6udFN4XHABcPSoe3lJ0xs/JpJNHfK7xMeiRWpw1oukJrQV2uBLbnBv4l4voKNDqXG6u9Xf\nAwdkun6amEg2dcjvEi8/+IGay+BGUhPaCm3wJUeIN0nPrpQFaNJn7NjoXv7QoQOvZWA3GlWZ0GYS\n90lrSzuGX41Ft2uJMIsqmw6AS8rd9HDLoBl2a2hQx5OxL3N090QcKdohg7Z6xOjoCdIYgxqBMA8U\nIRgmipDKrVxWm9f++npl+N32ETFPny73jk1aD0Yx+BpsQzZypGrc4qlEI+zDUxagSRZTOaZz6+xU\nRjuqoqfo99CSJd4PzrgfjGLwfah86hKpH6a5WYyOH35dU1kkI5uEkWO2tYmMMyoLFvj3kuJ+MJoa\n/MIN2rpJMZlVvpdDh4qXW8R04E0361YGwLNJUDmmPWAuMs7w2DbGK4eUk7Ql4IUz+CLFHMA0dYLJ\nfIVhw7zPI6qb6hFUjmkrsPy+V19vdsyiPuhN14GwSdPuFM7giyeq8DPit9+ucqvYXofuIXn99cCz\nrvlPFbJIRvXwkv4NGgTU1XnLAf0kg9/5jloq0UtDblPUB73JQjRO0rQ7hVvTVpZrU/gZcWaVnnXt\nWnXT6x6S27a5X09AGXtZJKO6eK1hC/iva6v73ujRwNtvA/v3u7ePoj7o/WyMG6naHZNAf1pbGoO2\nIsVUmA7KNTYqNYHXfAUTxYeQH9xkhnayNZHXKoJKYeOwO5BBW3dkuTaFX8zdSW+vapaDB4c7z4YN\nMhMzL3iFAY8eVffQtGmyGhbgbmO8IFIzoFPD5KmQ1lYNHb5T/12UNMkLFqjJM6YeyMiRyouLItEr\nuteXB0R6G4x161QPV9c7juPegOjwg1GUqeJBu5v2HIWomuwihs3yhi4MWMXbN7MEmfgW5d4wNfiF\nC+m4UaQ0yUElY8xmemJTiiZ9zROS8C44QRQ7adwbiRt8IrqEiLYT0S4iuiXp84WhSNr8oJKxMJTL\n3vuKJH3NG0lnUM0jQSawpXFvJL2mbRnAjwHMBHAWgFlEdFaS5wxDkbT5fg2QCDjttGgzLEsl4PTT\nxRPMI26DkUTqAX/KKcoxylNvOA6CTHxL495I2sP/NIBdzPwyM38AYDmAKxI+Z2CK1FX1a4AnnQSc\nc060HkBDA3DbbeIJ5hXn6mUjR6r7o78f2Lcv3CL3eSeIYieNeyNpg386gFcd7/dan2WKInVVdbJU\nv5WobHQzNa+/XqSveaajQ90z77yjxncqx70WL1YztUWOq6hc4rO1VdmVqtwbJiO7YTcAfwdgqeP9\ntQDuqygzF0A3gO6WlpZwQ9QxULTc7F5pif0WhQGYhw49UcbqldpYUh/nF1P1Sd7vo7DEfW/AUKVD\nnOAIHhF9BsD3mfli6/13rYfM/3Qr397ezt3d3YnVR8f69f5TzYvA+vWqW3748In7GhuVp1K0ayKc\nSFubSrhnirSdZCGiDczcriuXdEjnTwDOJKKxRDQYQBeAFQmfMzSVC20XsXHKTGTBhKDpk/Omdquk\nVtb3TdTgM/NRADcCeBLANgCPMvOWJM8ZN7XyQ8ZJZczRb6p8Ea+PEDztct7Ubk5M04xnApO4T1pb\nNWfaulGU2bdhketTbCrHvXQztmfNqnaN4ycryRghM22jUaTZt2GQ6yM4e4JNTf5lS6V8qd1sam3S\nphh8D2rth/QjibBLnq6PEB573EuX8XH48HyO/9TapE0x+B7U2g/pRVLxxbxcHyEedDO4p05Ntz5p\nUWuTNsXge1BrP6QbSYZd8nB9hPjQzeDOYzgHqL1Jm2LwPai1H9KNJMMuebg+Qnzo5LxAPtVctSZj\nFoPvQa39kG4kGXbJw/UR4sVLzgvUkGwxBEFkzNUm0Zm2Qan2TFs3ann27axZ6ubyWrD96quBhx8e\n+B937FChGN3/6Cw/dKiK0b71Vu1dHyF5ZOZ2OpjOtK269t65ZU2HX+uYaISDaulFe188oiz9Kcsi\npgNEhy+YxFWDDOqK9r54RFV5iZorW4jB96BWUgbo6ukXXww6qCva+2IRxwNe1FwZw6QbkNaWlZBO\nrYQtotYz6KLUsoh1sYgjHJOV1AN5BxLSCUethC2S9r4A4OWXj+81iLdWLOIIx4iaK1uIwa+gVsIW\ncdRTl/Hw4MHjY7aivS8WcT3gTWWLtRJGrWlMugFpbVkI6dRK2CKueppmPKxU9RRlZbAik2Y4plbC\nqFkFEtIJh4lXkwVPJAnvyy/jod1rqKVJJkI00grH1EoYNReYPBXS2rLg4a9bx9zQ4O7VlMvMw4ap\nvyaeSBT9skk94/a+aqV3I6RL1PVXdfeBaPWjA0MPv+pG3rllweAvWMA8aJC/4fMzsnbjHjnS/MEQ\npa5xhlf8FjAvlfK5gIWQLCahGnE0olN1gw/g+wBeA7DJ2i7VfafaBt/Pa9ZtpRJza6v++2G9by8v\nKar3Zfr/i4ROCIppe/JzNADmlhZpezqyYvBvCvKdaht8v66lqdE3KRO0i5rmgJYMygpxYRqq0Tla\nRNIGdZgafBm0deCnOzbBLUmZW5kg08njHtCKMjNXEIJgquN3Dg4TnViWWQZw4yJpg38jEb1ARD8n\nIs2ql9VHNxEpDoJOUPLT2x8+rDJemt4EpnlR7GXrurvVX5kcI4QhiJLMdjRGj/Y+XpbmwdQqkcwb\nEa0mos0u2xUAfgrg4wAmA9gH4G6PY8wlom4i6u7p6YlSncjoJiLFgekEJdsTf/xx/17HK6+YJbMS\n6ZuQNkEn6nV0AKee6n28LCZby4JEOxAmcZ+oG4AxADbrylU7hs/sHsMul8PH9cPEIStj9kFUQl6I\n9E2oBkHHhLKoFPMSTGRpshgyMGjb7Hj93wEs130nCwaf+UTly/Tp4QZzidTDornZW0FT2ZiWLAmn\nFNIZbZG+CdXCTUnmpzrLklLMy6jPmZOtembB4D8E4C8AXgCwwvkA8NqyYvAr0RlLt+0jH9HLJN0a\nU5TehN/ly6LnJBQTnWecFaWY38PHT5FXjR6zqcFPbIiSma9l5onMfA4zf4GZ9yV1rqQJM5j7wQfq\ne14Dnl4x9WPHwtVRNxgsic+ELOA3lnT77cCMGcDMmdlQivkJJvwUeVkca7ARWaYBYQZze3v9B0P9\nGpMXREC57L5PZ7QlTa2QNm4Dmn7tnhlYtUqJEFaurL5SLKxMO9Opwk26AWltWQ3pMJtnlTTt2oUJ\nEzljh2G7u3HOzBUEL7zCNs3N5m09rbbpNZ6gmwHstZXLzOPHp3t/odox/DBblg0+80DDaGoy//G9\n/iVdYyqXvY26zmgnmbRNEHT4xb6d+aXCOktxUvlgcgotOju9EymaOnxpjT2IwU8Q05w7foOhOjXC\nkiXhPPEsScWEYuInAbYNahRnKS5MUjoMGsRcVxfO07e3+vrknS4x+AljopXXdUvjViNkTdImFBNd\nuHLkSL3DlJRyzNn7bWkx623U1ytpdlub8voHDw5u9Ftb4/9fnJgafBm0DYk9Ffzqq4GRI9Vgqp0H\nxHQwNO68NbWyPKOQb3QpFaZOVe28s9O7XBLKscrUIq+8osyxjr4+YNgwNYD82c+a5cyqZNu2jMzC\nNXkqpLXVkodfSRYGQ2VylZAFgvQ0db3cuMajoqQ+t++dsLp8e0tyTAIS0qkeYRtp1MbtNxBMpLqw\nMpArpEGQcKWXs2Qfw9mOy2WlVAtKlNTndnhJNzZh8tBICjH4VSLsoGkcg60mg1BxjBUIgglRer26\ntuxn9N0cpzAyaOdDZskS/TH8FD1Jz2YXg18Fwg6axjnYWulZRRlUFoRq0dWlN8JubdfLcWptDa+0\nsRMf+h2jVFIDu15GP+l7zdTgy6BtjIQdNNXlvL/tNv257VmNv/udGliaNk0NBI8e7b6ohK5OglBN\nduzw33/s2Ilt1y9tw65dwODB4erCrD9GQwPwgx8AN9+c7dnsYvBjxHSFnyDfA4CnnvLPd1+pPnjq\nKWDtWuCyy1R+8TB1EoRqYpKaoLLt+jlOfX3Axz+ujG9YnMfwMuhZXzFODH6MBFnhx/R7gDLMXnl5\ndAubDBsWrk6CkASmC4bMn++dNwpQvdbKtqtzuBoa9HJQr95w5TH8DHqmV4wzifuktUkM33vzmmqu\nW9iks1MmYwnZIKgwYc6cYONPQVKAu6mI6utVWvNqDbxGATJoWx3Czp5dsEAv7XK7PCba+6zkFxeK\ni59TQ6QGPN2cjzlzjk/F4LdyXFCHy6kiam3V583JsoMkBr+KhJWjTZ/u3di8vAtT7X1n58D0cNHh\nC2ljooMPqtN3I4xzYyJnzrqDZGrwSZXNBu3t7dzd3V3talSN9evV4Ovhwyfua2xUg1Jr1qhY5bhx\nAwO1a9Z4xy6J1L5SScUfb7opOwNIQnFoa1OCAh2Njao9R4l723n3d+4EzjxTjQdUHs8us2MH8Oab\nwKuvet9DLS3qPstULL4CItrAzO3agiZPBa8NwJUAtgDoB9Bese+7AHYB2A7gYpPj5cXDj4KXh3LB\nBcEXNq+1bqmQX0xzy6eRFtkk8aEulJo1kJIOfzOALwF4tuJpcxaALgCtAC4B8BMi8hlzF2zcZF33\n3qtG/J1KnLCI9l6oBqarxpkpudmvAAARK0lEQVRIhU2VPm7l3FRtfuROyWbyVNBtAJ6Bw8OH8u6/\n63j/JIDP6I4jHr47UfKA1KrHIuQPt9w4QZUwpkofvxm3Qe6lWukRo8ozbU8H8Krj/V7rs0Ji6pF4\nEXZtTTdy57EINYPdew2bFlk358S+r/zKbd1q7tlnaYZsXGgNPhGtJqLNLtsVcVSAiOYSUTcRdff0\n9MRxyExROQv20UfVe7+Zs5XoJmYFIYk844JgSkeHWqj8H/8xeAoC09QluoXSvSBSA7RZnCEbGybd\nAN0GCem4EldStKi5vGtFWiYUi6DyZdP1HnTlvAaPdfdklteKRpo6fBeD3wrgeQD1AMYCeBlAWXec\nvBl8Xey9pSX4WrWVjdUZo6xcf9O5IHPWGqggBMVkJu26deq+8jP2Z58dXKuf9bWiTQ1+JB0+EX0R\nwP8GMBzA2wA2MfPF1r5bAcwBcBTAt5h5pe54edPhm2iP7W6sSdexUl88bZrqdjr1xoBegywItYhu\nnkpXF7B8uft+Z7k1a9Rr0/tEd96o8wbiIBUdftxbHj18E+1xrSgBBKHaeM1TmTMnudmyunxVSc8b\nMAGSD7/6mGqPRRsvCAqdos0r/fDhw94DtYBaFyLsIGzYtOdZRAx+gnR0qHCNLge3X6OJKukUhFrB\nVNFmpx++7z71/oYbgCef9FfgDB8ePuwSNu15JjHpBqS15S2kY2MykOQ22STrA0WCEBdBFW1B0iM4\nkwiGES/EuQRpUkCyZWaLMKlbs97IBCEuTOLktixy/PjjUyabypKjOE1ZTzFuavAlpJMSzvCOyWST\nsOvjCkItoouTP/PMQLhn+3a1pq0O5+pV9rHdZuaakPWlC00ZVO0KFIlFi4CZM83kYHkaKBIEHePG\nAZs2qbZdCRHQ02Nm5AGgqQk44wz1Ha+0x7bTFCSu39FRffllVMTgp4xpo/G7AWpuoEgQNMyfD6xY\n4a51L5Xc7wM3SiXgkkuAhx9Wnvgrr7iXK6rTJCGdFHEqbqZPB2bM8Fbf+Ek6JR+OkDf8Qp4jRvgr\ncJw4741cqWviwiTQn9aW50FbP1WBLsVrVgeKBCFu3PLrmExgdLs3iiR8gCxxmB38pmY7cZumbbJc\nmyDkGb/7p1xW98Xkye73xsKFaoC2t1eFcfK61KdpagUx+Ckwa5ZSF+gudamkRv+XLUunXoJQK0Qx\n3FGcJufat+PGZdfhEoOfIUwXcAaUwuCJJ7LZqAShmqTd27UfMkeOKGcty70DMfgZYtYsNU3cVGkQ\nJIOmIAjxUwsZMp2YGnxR6aSAaRI1mzATQwRB8CZoTqq8TnwUg58CbpIzHbXcqAQhS4RZZjSvEx/F\n4KdE5dTszk5gyBDv8rXcqAQhK5gufF5JXjX8YvBTxE7r2t2tFnK+7LJ8NipBSAtdqCZsaCavEx8j\nGXwiupKIthBRPxG1Oz4fQ0RHiGiTtd0fvar5I6+NShDSoDJUs3w58NnPAmefPWD4w4ZmgiY7rBWi\nevibAXwJwLMu+15i5snWNi/ieXJJXhuVICSNW6gGUK+3bAGmTlUPhCihmbxkyHQSKXkaM28DAHLm\nIRUCESSDpiAICr9QDaDCNXfdpcp5JWUz6UXnIUOmkySzZY4loj8DeAfAPzLz/03wXDVN3hqVICSN\nX6jGprdXeeQ33eQ9S7do953W4BPRagAjXXbdysyPe3xtH4AWZj5ARG0AfktErcz8jsvx5wKYCwAt\nLS3mNRcEobD4pQ+3sWP0y5ZJL9pGa/CZuTPoQZn5fQDvW683ENFLAMYBOGEaLTM/AOABQM20DXou\nQRCKh1/+fBtnjF560YpEZJlENJyIytbrMwCcCeDlJM4lCELxsAUP9fXeZUTpdiJRZZlfJKK9AD4D\n4PdE9KS16wIALxDRJgD/AmAeM78VraqCIAgDLFqk1rptbT1eiSNKN28keZogCDVP0deNME2eJmva\nCoJQ80iM3gxJrSAIglAQxOALgiAUhMyHdPr6+rB371709vZWuyqFoaGhAaNGjUJdXV21qyIIQoxk\n3uDv3bsXQ4YMwZgxYySFQwowMw4cOIC9e/di7Nix1a6OIAgxkvmQTm9vL4YNGybGPiWICMOGDZMe\nlSDkkMwbfECSs6WNXG9ByCc1YfCDEHTtShPK5TImT56M1tZWTJo0CXfffTf6NSuS79mzBw8//HD0\nk2u4/vrrsXXrVt8yv/3tb7VlBEHIP7ky+GHWrjThpJNOwqZNm7BlyxasWrUKK1euxCJNUuy0DP7S\npUtx1lln+ZYRgy8IAgA1SJeVra2tjSvZunXrCZ+5sW4dc2Mjs0qaevzW2Kj2h+Xkk08+7v1LL73E\nQ4cO5f7+ft69ezeff/75PGXKFJ4yZQqvXbuWmZnPO+88/shHPsKTJk3iH/3oR57lnOzevZvHjx/P\nX/nKV/iTn/wkf/nLX+b33nuPmZlXr17NkydP5rPPPpuvu+467u3tZWbmCy+8kP/0pz99WM/vfe97\nfM455/B5553Hr7/+Oq9du5abmpp4zJgxPGnSJN61axffe++9PGHCBJ44cSJfffXVrv+z6XUXBKH6\nAOhmAxtbdSPv3KIY/K4uZiJ3g18qqf1hqTT4zMynnHIKv/766/zee+/xkSNHmJl5x44dbP8PTz/9\nNF922WUflvcq52T37t0MgJ977jlmZr7uuuv4zjvv5CNHjvCoUaN4+/btzMx87bXX8j333MPMxxt8\nALxixQpmZv72t7/Nt99+OzMzz549mx977LEPz9Pc3PzhA+PgwYOu/7MYfEGoHUwNfm5COmHXroxK\nX18fvva1r2HixIm48sorPUMnpuVGjx6Nz33ucwCAa665Bs899xy2b9+OsWPHYpyV63X27Nl49tkT\nV5UcPHgwLr/8cgBAW1sb9uzZ43qOc845B1/96lfxy1/+EoMGZV6ZKwhCTOTG4EdZuzIoL7/8Msrl\nMkaMGIF77rkHp512Gp5//nl0d3fjgw8+cP2OablKhUwQxUxdXd2H5cvlMo4ePepa7ve//z1uuOEG\nbNy4EZ/61Kc8ywmCkC9yY/Dnz1f5r92IMy92T08P5s2bhxtvvBFEhEOHDqG5uRmlUgkPPfQQjh07\nBgAYMmQI3n333Q+/51WukldeeQXr1q0DADz88MM4//zzMX78eOzZswe7du0CADz00EO48MILjevs\nrEt/fz9effVVTJ06FYsXL8ahQ4fw17/+NdS1EAShtsiNwbcXRGhsHPD048qLfeTIkQ9lmZ2dnZgx\nYwYWWtKfb3zjG3jwwQcxadIkvPjiizj55JMBqLBJuVzGpEmTcM8993iWq2T8+PH48Y9/jAkTJuDg\nwYP4+te/joaGBvziF7/AlVdeiYkTJ6JUKmHevHnG9e/q6sKdd96JKVOmYOfOnbjmmmswceJETJky\nBd/85jfx0Y9+NPzFEQShZsh8Pvxt27ZhwoQJxseo5bzYe/bsweWXX47NmzdXuyqBr7sgCNWjsPnw\nJS+2IAiCO1GXOLyTiF4koheI6DdE9FHHvu8S0S4i2k5EF0evav4ZM2ZMJrx7QRDySdQY/ioAZzPz\nOQB2APguABDRWQC6ALQCuATAT+xFzQVBEITqEMngM/MfmNnW9K0HMMp6fQWA5cz8PjPvBrALwKcj\nnCdKNYWAyPUWhHwSp0pnDoCV1uvTAbzq2LfX+iwwDQ0NOHDggBihlGArH36Dl8ZVEISaRTtoS0Sr\nAYx02XUrMz9ulbkVwFEAvwpaASKaC2AuALS0tJywf9SoUdi7dy96enqCHloIib3ilSAI+UJr8Jm5\n028/Ef09gMsBTOMBN/w1AKMdxUZZn7kd/wEADwBKllm5v66uTlZeEgRBiIGoKp1LANwM4AvMfNix\nawWALiKqJ6KxAM4E8P+inEsQBEGIRlQd/n0A6gGssnK4rGfmecy8hYgeBbAVKtRzAzO75xIQBEEQ\nUiGSwWfmT/jsuwPAHVGOLwiCIMRHplIrEFEPgH+PcIhTAbwZU3XiROoVDKlXMKRewchjvf4TMw/X\nFcqUwY8KEXWb5JNIG6lXMKRewZB6BaPI9cpNtkxBEATBHzH4giAIBSFvBv+BalfAA6lXMKRewZB6\nBaOw9cpVDF8QBEHwJm8eviAIguBBTRl8IrqSiLYQUT8RtVfs0+bfJ6KxRPRHq9wjRDQ4oXo+QkSb\nrG0PEW3yKLeHiP5ilet2KxNzvb5PRK856napR7lLrOu4i4huSaFenusqVJRL/Hrp/ndr9vgj1v4/\nEtGYJOrhct7RRPQ0EW217oH5LmU+T0SHHL/vgpTq5vu7kOJ/WdfsBSI6N4U6jXdch01E9A4Rfaui\nTCrXi4h+TkT7iWiz47OhRLSKiHZaf5s8vjvbKrOTiGZHrgwz18wGYAKA8QCeAdDu+PwsAM9Dzfod\nC+AlAGWX7z8KoMt6fT+Ar6dQ57sBLPDYtwfAqSlev+8DuElTpmxdvzMADLau61kJ12sGgEHW68UA\nFlfjepn87wC+AeB+63UXgEdS+u2aAZxrvR4Ctf5EZd0+D+B3abUn098FwKVQmXQJQAeAP6ZcvzKA\n16G06qlfLwAXADgXwGbHZ/8E4Bbr9S1ubR7AUAAvW3+brNdNUepSUx4+M29j5u0uu7T590nlfrgI\nwL9YHz0I4G+TrK91zqsALEvyPDHzaQC7mPllZv4AwHKo65sY7L2uQtqY/O9XQLUdQLWladbvnCjM\nvI+ZN1qv3wWwDSFTjleBKwD8H1asB/BRImpO8fzTALzEzFEmdYaGmZ8F8FbFx8525GWLLgawipnf\nYuaDUAtOXRKlLjVl8H0wyb8/DMDbDsMSOkd/AP4zgDeYeafHfgbwByLaYKWJToMbrW71zz26kbGt\nZRAS57oKlSR9vUz+9w/LWG3pEFTbSg0rjDQFwB9ddn+GiJ4nopVE1JpSlXS/S7XbVBe8na5qXC8A\nOI2Z91mvXwdwmkuZ2K9b5hYxJ4P8+1nAsJ6z4O/dn8/MrxHRCKgEdC9a3kAi9QLwUwC3Q92gt0OF\nm+ZEOV8c9WLzdRViv161BhH9DYB/BfAtZn6nYvdGqLDFX63xmd9CZapNmsz+LtY43RdgLb9aQbWu\n13EwMxNRKnLJzBl81uTf98Ak//4BqK7kIMsz88zRb4KunkQ0CMCXALT5HOM16+9+IvoNVEgh0o1i\nev2IaAmA37nsMl7LIM56kfu6CpXHiP16VWDyv9tl9lq/8SlQbStxiKgOytj/ipl/Xbnf+QBg5ieI\n6CdEdCozJ5o3xuB3SaRNGTITwEZmfqNyR7Wul8UbRNTMzPus8NZ+lzKvQY0z2IyCGr8MTV5COtr8\n+5YReRrA31kfzQaQZI+hE8CLzLzXbScRnUxEQ+zXUAOXm93KxkVF3PSLHuf7E4AzSSmaBkN1h1ck\nXC+vdRWcZdK4Xib/+wqotgOotvSU1wMqTqxxgp8B2MbMP/IoM9IeTyCiT0Pd34k+jAx/lxUA/oul\n1ukAcMgRzkgaz152Na6XA2c78rJFTwKYQURNVvh1hvVZeJIeoY5zgzJSewG8D+ANAE869t0KpbDY\nDmCm4/MnAHzMen0G1INgF4DHANQnWNd/BjCv4rOPAXjCUZfnrW0LVGgj6ev3EIC/AHjBanDNlfWy\n3l8KpQJ5KaV67YKKVW6ytvsr65XW9XL73wH8AOphBAANVtvZZbWlM5K+PtZ5z4cKxb3guE6XAphn\ntzMAN1rX5nmowe/PplAv19+lol4E4MfWNf0LHAq7hOt2MpQBP8XxWerXC+qBsw9An2W//gFq3GcN\ngJ0AVgMYapVtB7DU8d05VlvbBeC6qHWRmbaCIAgFIS8hHUEQBEGDGHxBEISCIAZfEAShIIjBFwRB\nKAhi8AVBEAqCGHxBEISCIAZfEAShIIjBFwRBKAj/H/UrvZYf/kE0AAAAAElFTkSuQmCC\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.59612873946671\n" ] }, { "data": { "text/plain": [ "Text(0.5,1,'Polynomial of Degree 4, cost=38.596129')" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAEKCAYAAAARnO4WAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzsnXmcFMX1wL81s8seciOXwgqILIcK\niMqiRkQRRU2MMSKoxNv4UyMeeCOHR8SIGkw8EokaL7yveEZQNOCCrohyyA0iKoIcgsACu/N+f1TP\nbu9sd0/P7Fy7U9/Ppz+7011dXdNT9erVq1evlIhgMBgMhoZPIN0FMBgMBkNqMALfYDAYsgQj8A0G\ngyFLMALfYDAYsgQj8A0GgyFLMALfYDAYsoSMFfhKqfOUUmI7timlvlRKXaGUyokxr05WHuclqbgp\nQSk1XikVsx9tor+/UiqglPqrUuoHpVRIKfWaR9rVtt+wQim1USk1Ryk1USnVKRHlqS8opboopXZY\n76JrussTDaXUVUqp3yUwv8eVUl8rpbYqpX6x2vOflFLBiHSFSqkJSqmlSqmdSqlvlVJP+qkvDnIj\nfMxzSDtIKTXTesYmpdRTSqm2LvmWKKXeVUptUUptV0rNV0oNj0jzZ6XUf6067tjelFLtlVJ3KaXK\nrLw2KKWmK6WOdkgbVErdqpRapZTapZRappS6Kto78CImwZkmzgDWAk2t//8GtAHGprNQaWIK8G66\nCwH8HhgFXAuUAhujpH8PGA8ooDlwCHAxcIVSaqSIvJq8omYUDwE/AwXpLohPrgJmAq8kKL8CdPtd\nAQhwAjAZ6IquT2GmAL8FxgFlQBEwAZiulOotIr/4eFZYboTZbr+olPoV8F903TwdaAXcYT2jn4js\nsqU9GXgVeBY4C9gN9ATyI575J2Ae8CbwB5dy9QPOBB4HZgONgMuAGUqp34jIm7a0DwHnAbcDc4BB\nwCSlVGMRuSP6K3BARDLysL6oAF0jzn8I/BxjXp2svM5L9/dK07tM6PdHN0QBAj7SrgaedjjfGJgF\n7AA6pPh95KXhNzgL+BEtRGvV60w83H67BD9jKrDN9rkQqAD+HJHuROu9nRAlP0e54ZBuGrAcyLGd\nO9S69zLbuSbAeuCvPr5LwPrb1a29oRWenIhzOcAS4GPbuSKgEhgfkfbvwE6gZTzvO2NNOh58BjRV\nSrUBUErlKqXusEwHu62/dyilct0yUEpdaw2RWkecV0qplUqp56zPYVPIH5VSt1kmjC1Kqf8opTpE\n3Bu1HLb8LrWGdessU9XT1jC2q1LqPWu4u1wpdW7EM2qZdJQ2cZVaQ9ItSqnZlkYSF0qpE638diql\nflZKvaaUKrZdX43W1gEq3Yau0RCtpV2G1vr+GFGGgdYwd5s1fH5PKXVgRJqg9X5/UNpM8oFSqrtV\nnvG2dOOtcweG3y3wgu3676x3tsN6fy8qpYoc3sslSpsgypVSPyml/qWUaunnuyqlWgD3AaOBLf7f\nkmt+vZVSr1qmg51KqSVKqZts15VS6mrr/G7rHf1dKdU0Ip9RSptYdiqlNittZjjNurYa2A84W1Wb\nRZ6oa9kd2IgW8GGC1rE1Il34vSVKZpUA74tI1bNFpMwqz2m2dGcArYF7o2UoIiEfabbYn2mdq0CP\nDPa1nT4c/V3ficjiXfTIYmi0ZzlRHwV+Z3TPFx7W/Ru4EXgSOAV4ArjBOu/G40AIOD/i/BAr/0ci\nzt+E7rUvQA89BwBPR6SJpRw3AfsA56JNU2daz3wVeAtd4b4CHldK9fL4HqC19ynoinkmegj8plLq\nxCj31cK65y30uz0T+D/gQGCmUipcGU+zvhvo9zDAuidmRORL4HvgSFsZTgamW2U4B60ZNwH+p5Tq\naLt9AnAz+n2fih6ev+HxuNeBj4DfAPdbz7oUeBlYhDZT/dH6vh8ppZrYyjQReBCtFf4GuA6tcb6j\nIuzPLvwFWCwiT/lI64lS6nC0GW1/4GrgZHRnYldA7rTOvQ/82nr+ecBbSqmAlc/ZaCE2FTgJOBt4\nCQh3YqcB69Amj/DvfLt1r1JK5fg4ar0b273NlVKno9vAfeHrIrINeAq4Umkbe2OrDdwDfImuG36Y\nqZSqtDq7Rxw650q0aSaSXeg6EOYoYBNwkNJ2+wql5xTG+fzto6KUaoR+v19HlA+HMoZNTQcSD+ke\nNvoYmhWjhzwt0A2yEnjNSnOglSZy2DPGOn+w9bkTEUMstNBaDijbuVeAr22fw/fNiMh/tHV+nzjL\n8UFEules8+fYzrVAaz7jbOfG65/MfUhpvav/Aq87fI9aQ8yI+8uAZdQc5nYG9gD32c7d4VWOiDxX\n42EWQAsv+ztfDkyPSNMU+AlrWG29m1+AhyLSXRP5O4TfGTAqIm1jtD39sYjzndGN7Crbu6sExkak\nO9LK97dRvv+v0I20Z0S9jsukA3wMfAsUulxvaT3viYjz51jP/Y31+e/A3Hh+O+AYK69oxwyHe0+x\nXQ8RYbqx0gTRHaw9r9lAax/v5wS0MnAS2uY9BtgGzAfybek+BeZE3LufVaZdtnPvok0oW9BzVsdY\n9b8CuN+lDK4mHZf0f7ae+yvbuZ5WHv8XkXasdf4fcdWfeG5KxWFrGPajEq3RtbTSXObUeKgWcH+K\n+HyeLc3h1rnB1uf2aMF2tUM+1ztUKgFK4izHdQ4/uERWaLT2O8X2eTwRghY9CfQm2j4csr2rxQ7l\ncK2AwF7W/Xc4XJsBfG77nEiBPxtYZP1/gFXOC9Adl/34D5aAAo620g2KyKsId4FfFJH2eOv8cQ7P\n+gp4xUp3sZVuf4d0W7F1hA7frRF69DDRoV7HLPDR9u1Ke34OaU7CVq9t53PQ9fte6/O51u/9N2Aw\nDh2I22+HHnEd6uModri3mXXtOHS93w3cGZHmLnRnfK31W5+DVkTKgL3ieG+nWu/kItu5s61zd6Cd\nQLqjO9MKYKct3X+tdNdE5PmwVfZmDs/zLfDRI9gQMMHh2vvo+YMT0Lb/09CjDQEejvU9iEi98NI5\nDT3bvg34RkTKbdfCw7QfIu5ZF3G9FiLyqVLqc+BS9FD9IvSP7WSC2RTxOTysCs/Sx1qOzRGfd3uc\nj/QEqMIycUxHC5U/AWvQ3+F2oIfbfS60QHvRRH4H0N9jvxjz80tHYLH1fxvr77+sI5I11t/21t/1\nEdd/9HhO5PcKP2uaS/rNEemWu6Rr5fHMq9Dv9QGlVHPrXKH1t4lSqoloE4ZfWqBHcWs90jjWRRGp\nUEpttF1/El23LkQrLHuUUm+jBdvqKOX4BW1zjobUOiHyM1pwg/aI2Q3cqpR6SES+s8w3N6KFc1Ud\nUErNAZai2+lkH8+28wbaS+cwtPkTEXlGKdUdPVq/xSrr88Db1DSXhD3Q3o/I879o2dEL+CTG8gCg\nlPo12tLwLxEZ55DkPOAZqj3ztgLXo82/Tu00KvVB4C8QEbfGFhbE7dCuXtg+26+78RDwD8s+fRHw\noohEuycZ5YiXE9Ea0zARqRICSqlC91tc2Yyu9O0crrUjCd9BKdUHPZcxxToVblw34SyIwx1juLK3\nARbarjv6UFtECp/ws86LyCPMtoh0Q6jdIduvO9ET/e6+c7g2F22T7uNxfySb0drgvh5p7HWx6nsp\nvXalVfi6aBXyH+j63wL9/e5FC73+UcoxEO0tF42P0CYQL8rQnVhn9Hs6yDr/mT2RiCxTSm0hdkWm\nRjYRed5qzc90AdaLyI9Kqa/RrqhhnOqGnagTtU4opY4DXkTP2/3RKY2IfAcco5TaB91RrwAOti7P\ndLonGvVB4HvxsfV3OHqiKszZ1t8ZUe6fCkxC+9cWUXuyNlXliJewYN8TPqGU6oa2L3tpgbUQke3W\niOcMpdR4Eam08tsPOAI99E8YSqnGaDvtDrTgAe2athroJSITPW6fj9bYzqCm4DkjhiJ8ghbqXUXE\na4L/fXSjLhKRSC0vGhOpnuAOcyJ6Mv8c9Pf1jYjsUErNBM5RSt0mIjsdks1Gd4zDqTnBeSa6vc9w\nyHcz8LxSqj81hc8unNcMfI7WlqPhZ/QyEC2IV1qfw6Piw9GmNaCqXjfHufOMxm/RJstPIy+IyHZ0\nfQo7LXRHj3rCvIYeMZ8QTmdxIlAOLIi1MEqpAWgngunoeTvPTkNEvge+V0op9KhxMXHKlHot8EVk\ngVJqKjDe0mA+Qc923wpMFZH5Ue7fabmaXQ3MF5G4hmZ1LUcdmIY24TyplLoXbeqYgDZ9xOOBdSva\n4+ZNpdRD6InNCWh7alS3NA/2VkqVoE1GzaheeNUaGGFVaERElFKXA69bngsvoCdr26I7nTUicp+I\nbFZK/RW4WSm1Df0eDqG6ofpxj9uqlLoOeFBp99x3rO+5L1oIzRCRZ0VkhVLqbuDvSrunfoRu6B3R\n8wBTRMRR2xWRxVSbqwDtmmv9O8c+clXatfVx9LzEDI+ij7bKUGr95mvRGmofEfmTiGyyzt+klNqO\nNlH0QNuqZ2J5VCml/okWyKVo01g3YCTaVBFmEfArpdQpaEH8k4istsxQZcSA5X11PnouZg16HmAo\ncAl6AvJ7K+n/0COfe62RR3jh1Rj07/NvW55j0ZOY+4vIN9a599FKwAL0ZOuR1jv7Em0eCd/b13r+\nXOvUUWjvq7/Y5YDVtp8AbrM8nOai5zwuAm4X2yIwpdRAdJ0Oj5IPVdoNGBF5yUrTHf0b/IT2POqn\n5XjV82bb8vs/dF1bZeV5rlXO46J1Eq7EY/hPxYH/BRSN0JX5G7Sm+431OdeWphPuCyEGWNcud7gW\nvu+iiPPHWOePibMckfmNt85HLshYjW3SDOdJ22FooVKOHn4OR2uVq/18f4fvfCJaCOxEN7DXiZh8\nI/ZJW/uk+2b0cH0isJ/LPQPQE9Gbre+1GngOGGBLE0SPptZZZZ2B7hQEm0eO27u1XT8JLSC2okcb\ny4DHsLxqbOlGorXn7Wgb9tdoT5eYFo3hvqDwcut8Dx959EULzi3Wd18M3GC7rtBKzBK0tv8DejTV\n1JbmXOudrUdr8qvQ7qr2NN3RAniHVbYn6tCeu6O90b61nvcjugM6m4gFfGjT073Wb7HTuud5h3oY\n/m072c791fpttlnffQV6FN8s4t5e1vPD73AucH4UGfOtledSIry+rHQzcPFYcvj9PdNZaa+wfsNy\ntCnuFfToN67fQES0S2I2o5S6E+1bv4+IRC72MNQjlFK/R9tFjxaR/6W7PLGglHoWaC4iJ6W7LIaG\nS7026dQFa0hXjBb2/zTCvn5h2ZtPRscYKUe7p96I1sLjmtBKM0ejR2sGQ9LIWg3fWjreFr2ScKTE\n5h5nSDOW+96DaK+OpmjTxH+Am0RPQhoMhgiyVuAbDAZDtlEfY+kYDAaDIQ4yyoa/9957S6dOnWK/\n8auvoGlTiOdegwH4/nv48UcI2ZzdAgFo2xb22Sd95TIYarB5M6xcCT17QkH1EonPP//8JxFp7XEn\nkGECv1OnTpSVxeTeqxk4ULfU/9UrxwxDhjB7Nhx3XE1hD/rzzz/DK69ASUl6ymYw1OAvf4EbboDS\nUq3kWiilvvFze51NOkqpjkqpD5VSi5RSC5VSo6zzLZVS7yu9Ldf71iKK5NCli+710I13xAjo10//\nnT07yr2GrMReT848E3Y6rVkFysthcqyRWwyGZLF6NbRsWUPYx0IiNPwK4FoRmat0DPHPrdVu56HD\n3E5USt2Idpm7IQHPq02XLvD999x+804mTi5g504QgXnz4I03YPRomDAhKU821EPGjYNJk6iqJ16E\nQrBsWWrKZTBEZfXqOpmu66zhi8gPIjLX+n8bepXbvuiQpOFl0P9Gx7NIDl26APD83avZsaO6EYdC\nsGOHbtxG0zeArgeTJlGjnngRCEC3btX3mtGjIa2kW+DbseKE9EUvhmkrIuGohuvwjmQYF+EGeMqV\nWuDvF1rpmM4Myw1hJk92N984kZ8PV16pRwXHHQfPPw9z58ILL+jP45yC2hoMyUAkcwS+Ff3wZfRO\nQTVWrYp29nfUp5TeK7RMKVW2YcMG388LN8DnnoNPN+0PQBecBb4ZlhvCLF3qX7MvLNTmQKg9KjCj\nR0PKWb9eayudO8edRUIEvtIbdb8MPCMir1inf1RKtbeut6f2ZhUAiMg/ReRQETm0deuoXkVAzWE5\nwAZa8wt7uQp8+7DckN1066brgxNKQVGRNtkMGwbTp+u5H69RgRk9GlLG6tX6bzo1fCtG87/Q+5Le\nZ7v0BjoiH9bf1+v6rDC1G6BiJV3Yv8beI9WEh+UGw6hRuj44UVCgTTZlZTB1arUrpteoIBSC994z\ndn1DCsgEgY+ONz0SOFYpNc86TkKHvz1eKbUMHT/aa0OLmHBqgCvp4qjhB4N6WG78qA2g68Ho0dpc\nE9b07eYbp3riNSoAvRbG2PUNSScs8PeLf7fROrtlishMdPxtJ46ra/5OdOumXS7tC2VW0oUhVfsN\n6+IEg/DII3DRRckohaG+MmECDB2qR4rLlsEBB2jN300pGDVKu/eGTYhu2O36Q4caJcOQYFat0j74\nTZrEnUW9jKXjNCxfSRcK2UlbfkQprbHdcosR9gZnSkq02SbSfOOWNnJU4IWx6xuSwsqVsP/+dcqi\nXgp8pwa4Ev0iDm+1kjPPrJ5wMxgSwYQJuk4NG6bt9S081o2HQvDFF6krmyFLWLEiOwU+1G6A3U/S\nvvhv3L8iqsZmMMSDfVRwwgne2v7y5caWb0gge/bAN99ULTKNl4wKnhYrJSU2wb6rExQGdC9oMCSZ\naHb9ykpjyzckkDVrdKXKVg2/Fnl52onarLAypICwWTEYdE9jbPmGuhKOJnD5iVqRXVhuBH41Bxyg\nfTYNhhQwYQJ07ep+3azwNtQFeziP0HIt8H977f51MhU2PIG/bJm/tfMGQwLo29fdlm9WeBviJTLI\n3/6soJw8VpTvw913x7/Ar+EJ/J9/hp9+SndJDA0Mt0iZXit3zQpvQ7xERhPYnxWspAtCgF274nc3\nb3gCH8w42pBQvCJluvnoK1WnGFeGLCcymsD+rGAF1fb7r7+OT8s3At9g8MApfn5kpMywi3CPHtVC\nX0Q3SrdQC35j65sY/NlJTVOg0IWVNQR+KBSfQ0DDEvidO2u3CSPwDQli7Fh318tIL5xVq2qG+3AL\noew3tr6JwZ+9jBqlR4kAbfmRxmyvIfAhPjFX7wV+DQ3oD7nsbN/ZCHxDQhg3DqZNc79u98LxG0LZ\nz4ghlnSGhklJCfTsqf8PRwG2C/x4HQLqtcB30oA+/v4AfvjYuGYa6kZY4Ho5fNkbnd8Qyn43TDcx\n+A1TpuiJfyeBH69DQL0V+G4a0OLQATRZt4zZpcY10xA/frZCzM3VDmH9+um/yi1mLNUhlNes8e4Y\nwiOGaB2IGcQ2fEpK4PrroXvOCkIoVtMpaijvaNTb0ApuDXIZB9CY7fz77nWUvNY+9QUz1Dtmz9b1\naelSrbGPGuVvK8RQSJt8RLSwr+vyD/uIwSkEuFM6Q8NmwgTYMGcFG//XkYN65EUN5R2Neivw3Rrk\nMrSnTsXXywAj8A3ejBunR4o7d+r6NG8evPqqjtThRSCg41mFsdfFeIW/fZjuFavH+PdnF61/XgH9\n96fsg7rnVW9NOm67EIUF/mHNzZjX4I2bWXDXLti61f2+YNBdoCsFHTtGD6Fsx2mYHs/OXIYGSgLC\nIodJiIavlHoMOAVYLyIHWufGAxcDG6xkN4vI24l4HrhrQGsoYje5nFxsBH6242SqsQtKP3b6IBUc\nxHy6sZQurKJ9zgYaN9rNzh0hNtGSDbRmKd34ioP5gfaIKFq31iGUR4zQjgROZpkwSmn//SlTagvx\nWHfmMjRAtm2DDRsyS+ADTwB/B56MOH+/iExK0DNqENaAJk3SXguhkNaA8vJz+LmgC/v+Yjx1shkn\nU80bb+g6E94Yx80s2J7v+T0v8RveoITZNGZ71bXKnEJ2VuSxC0VzthCkWpp/SwfeZShbC0+G3UMZ\nNapR1K0RRbT/vhs1QoAbso9wuPcECfyEmHRE5GNgUyLyioXITVCGDdOfWx/RzbgxZDFTpsCdd0b3\nYbebBRUhhvI27zGEtXTgAUbRjnU8zvkMZyoH8RW/6rON4M7tLPhoE0WFG2nEblqznqP5iCuZzBz6\ncybPce3/fgv77kvJS6O548JVFBZ6l9e4WRpcWb5c/02QwEdEEnIAnYAFts/jgdXAV8BjQAuX+y4B\nyoCyoqIiSQjXXCOSny9SWZmY/Az1hrFjRYJBES3qnY/Bg3Xa0lKRwoKQ/I6XZBHdRUDWso+MZbwU\n83WNewIBkREjaj6nsFCfD18vLBQZf8tukbffFjn9dJGcHJFgUH48+Xzp23S5Z5n69UvP+zJkOHfe\nqSvI1q2eyYAy8SOn/STylVFtgd8WCKJHEXcCj0XLo1+iav3DD+uv9s03icnPUC8oLdVC10uwhoXz\n2LEi8umn8u0+h4uALKCnjOAZyWG34z15eTr/yOcNH66F9fDhta/L2rUio0aJ5OfL7kAj+TM3SSG/\nOJbH3pkYDFWce65I+/ZRk/kV+Enz0hGRH0WkUkRCwKPA4cl6Vi26d9d/Fy9O2SMN6cfPJCxAfmg7\ne//5GqSkhA6sZcUtj3HnsK9Y2u8sinvlOu5iFQrBO+/UPGff49ZpH+XZ3+7LiB//ygkHrOTDtiO4\nibtYTHcG837N8sToZmkCqmURS5dCcXHi8vPTK/g5qK3ht7f9fzXwXLQ8Eqbhr1unVafJkxOTn6Fe\ncMgh0bX7fnwmy+kiAvLGvpfK+b/bIoccUq2hl5Zqa6DTvYWFDlq8C2GTj1LVWvzRwZmyiB4iIPcz\nSgrVDikstEYbPnHKN9Y8DPWIVq1ELrkkajJSadIBpgI/AHuAtcCFwFPAfLQN/w17B+B2JEzgh0Ii\nzZuLXHppYvIzpJSwqcQuiP0wfHi1Tb32EZLL+ZuU00i+oaP8io8kEKgtOHv1qj7nZHoZPtxf+d1M\nS80a7ZBXO/5JBGR184Nl7ssrY3ovbvnG0hkZ6gk//aR/3EmToiZNqcBP1JEwgS8iMmCAyMCBicvP\nkBLqosG6CcQ8dsrTnKW1ek6RFmx0HQG4CftYJleHD/fRabz9tkiLFiItW4pMm+br3fjK19BwKC3V\nP+4bb0RN6lfg19uVtlHp3t3Y8OsZdQ0J7LQ6tRU/MY3BnM2z3MydnMrrbKalax7h5zrhN4aNn8iZ\ns1sMhU8/hXbt4IQT4B//qHO+xhO5gbHUWkuUwMBJDVfg9+gBP/6owxQa6gWJCAlsX5vxmwNX8nmj\nARxKGWeqF7iLm1GBgGdUS3DflNzv5Kpb2I8wmzdbG5k81VX3YieeCJdeCnfc4dnjRMt3wwYzgdug\nWLpUx/FI4F6ZDVvgg9Hy6xGJ0mBLSmDq+CW8vvlo9muymeX//JBNx51BixbQrBm0aeMeyjgQ0BtP\n1CWGjdfG5mGqRi0Lm+hobSNHwq23wlVXucZiiJbvt9+aHbEaFEuXQpcu0KhRwrJsuAI/7Jr59dfG\njS3DCf8+XiEGAgFo2dLn77hwIQwcqMNZzpjBi2sH8MknsGWL1q7Xr3fvWPLz4dFHnVdwh0MyRMNu\nWvKiatSSmwtPPAFXXw0PPKCHEQ4FtOfr1GGJmB2xGhThIFCJxI+hP1VHQidtKypEGjWS/w24zrix\nZTCRk7RuR06OdpeM+jvOny+y9956scqiRVEXYyWzXpSW6nlZ35PAoZDI6NH6wjXX6M8u+RYVuedp\nJnAbAJWVukJefbWv5GT9pG0wyI4O3dg652uzL2iG4jRJG0kgoGPTK6U1Ys/fcdUqGDJED4E/+gh6\n9PCcF7CHMo5Vi/dDSYmej3Wzu0dOAs+eoxjx7V94rvUVcN99fHf+GNd8997b/blmArd+YrdEXH7a\n97qCGw3fP6Udz5CldDVaUIbi5WYIWjsePlzk+ON9uCOuWyey//76pgULqp4RbTFWsmPY+PWdt490\nFJXyTy4WAXl7yP2u785tzYEJ1VD/iBzpHqemi4A8PnK6r/vJeg0fmL+nO11YSR7lta4ZLSj9RNtG\nsE0bPVFZVuY9mfv94p+1Kv3DD/D229CrV9V1L8+WeLcKjGVOyM9GJpEjHSHApTzMK5zGCf+9hiUT\nX6mVr9cErtkRq37hNNI9QJYAcOeL3RJrifDTK6TqSLSG/8CAZ0VAejHfaEEZiPfKWK3t5OZ6a+i5\nao/Ma3eCNvK/+26tZyR6dWq8C8O8Aq25jXQK2C6fUCK7gvmOBXWL2Gnmp+oXTr//fVwlv1AoQVXp\nyxJB1q+0FZFXbp0rAvJ7XkhIYzckFr/RLb2OB3Ou1P9MmeL6nEQJxmSFNvAyO+3NelmTt7+eiF5Z\nOwxD1IidhozH6fd/k5NkLn2qTJvRftesF/hjx4q0KtguAjKG22pojUYLyhz8xK93O/6Ua4XBvuaa\nqPF3EiEYkxXaIJo9/upTlurYUL17i2zfHt9DDBmL0++/lK7yHMNqKBReMiurBb5dE1vFfvIMI6pe\nXDAo8uijCXmMIUEUF8cu7AcxXSpUUOSkk2TcmIparp1K6UBoidR4/U4A+w3+Fk5XXOze6VWNHN5+\nW3+ps85yddc01E8iR46NKJc9BOU2xvgeRWa1wLdrYm8xVL6gd0I0MUNyiGbLjzw6sEY20EoWqp7y\nxAM/e5qF8vMTN5rz4xnj18bvtv7A7b7SUpGpB+ndj57sc68x3TQwwvUBRHoxXwRkOM/6HkVmtcC3\na2J3c52U00iC7KmliRkyg1hs+bnsklkMkJ9pIsVqiRQVRV+0ZdeM4g29HK2chYV65OjHxu+VTzAo\n0r17zbJVdQ6E5CV+JxUEZHDex8Ys2cAIL9Q7g+dFQHrzhesoMpKsFvh2TWwk/xYB6c6iGpqYIbNw\nmljNyaktEO/jKhGqJ+L92P/DmlEiNg/xmgD2a+OPZS4gsnNozFZZSldZQwfpUPCT61xFPB2aIf0M\nHy4yTo2XSpTks8NxFOlEVgt8eyPpQ01PHeOdk7k4Tawef3x1hf8dL4mA/JUrfZt/wkdxceI8bNwm\ngP3a+GNZDObUOfTlc9lFrrzBr2X4mdX2fLMbVv2ntFTkxeAwWcb+MdXTrBb4ItWVv0DtlAoCchu3\nmspfDwl33p1YKVtoKqX0l1xtC4aJAAAgAElEQVR2xSTsAwHxNP0kal7H7+rXWFbJunUOV/JXEZC/\ndJxc4z0lokMzpJd1bQ6UN4O/jsmNOKUCH3gMWE/NPW1bAu8Dy6y/LaLlk/DQCpYmtiqvWOZ0OM1U\n+nrK+DF7ZFbgSNlMMylidczafWFhdE+gRFQ9v0LXT7pw3XUPvhaS/3CK7A40Evn8c7MbVkNhzx6R\n3Fz57pzrY3Ij9ivwExVa4QngxIhzNwLTReQAYLr1OaWUlMDUqdDp1wdxeP58X7HMDZnHuLyJHBGa\nxVMlD9G63340aeLvPnsIg759Ex9iIRI/YRT8pHvnHR3X/vnnvfbvUVxW8DjScm8YOZKFn5ej9aza\nmDAi9YiVK2HPHvYZ3JOpU3VYkalT/e3D4As/vYKfA+hETQ1/CdbG5UB7YEm0PBKt4VcxfrxWf8yi\nlXpDWMMdWTxHKlRQNhyv7RzRFmoppc03kZpRKk0efhd52f3wi4r038MOi+6iWmOI/847IiD3qOs8\n0xtHhXrCq6/qH+3TT2O6jVTb8B0E/hbb/8r+OeK+S4AyoKyoqCiON+SDl1/WX/Wzz5KTvyGhhOdf\nGrNNltJVVlMk7Qs2ywUXRHff9BLemRh7xu9+AOEjHEHU3pFNyblEKlFyBDMd78nL05PfxnOnHnCn\nXmshW7fGdFtGCXzr8+ZoeSRNw1+yRH/Vxx9PTv6GhGHXxB9BC7KjmSGgNXsvwRgM1i2IWaqJJ5ZQ\nixY1Bffw4SJN2Cor6STL2F8K+aXWiCc313ju1BvOOUekY8eYb8sEgZ85Jp2KCr3k8pprkpO/IWGE\nJx+P430RkL8w2rcwLC5Od+ljI9p+AH7MOu3b689HM0ME5G9cXkvgp8KMZUgQhxwicsIJMd/mV+An\nMx7+G8C51v/nAq8n8VneBIPQsydbZi0we9tmOEuXQqH8wqNczBK6MZbbfN0XCMAhhyS5cAkm2n4A\nXoR3/Fq/Xu/c9TEDuZ+ruIIHGciMqnRu+Vftp2vIHEIh+Ppr6NEjaY9IiMBXSk0FSoFipdRapdSF\nwETgeKXUMmCw9TltfFF5EDs+XcDzz8PcufDCC9oTYty4dJYq+4i2eUi3bjCRm9iPb7iAxyinoOqa\nUrrvdqI+bvrhtTmLX0Kh6jxu4U5W0IV/8EfyKHfc6Nx+X6TnTiwbuxiSwDff6P04e/ZM3jP8DANS\ndSTSpGNfYj54sMiNOfeIgLRgoxnapgk/K0EXPPyxCM6raQsLpWriNpMmXuMlEfsBgEi7dtXvZDD/\nFQH5c86t0quX/wVeZpVuBvDWW/rlz5wZ861k80pbJ8+HIbwrAvIrPqpl4zSLUpKPL7fI7dtFunaV\njc07y94Fv7gKdfvE6+DB9dsDJdJzyMlO70dw29/JR51GSmVOrsx7ZkHCFoIZUsA9WimVjRtjvjVr\nBb5b5d2HtSIgl/M3Rw3JkFx8rQQdPVqfmD7dlzdNQ9FKnTow+/eOWSCvXy/SqpXIEUfIuFsro46I\nvH4boxClkPPPl10t2sQV/C5rBb575Q3JevaWR7mw1rVg0GgxySZawLDhPebpH+Lii33ll21aacxr\nCP6to8TKww9H7Tyj/TZGIUoNa/c5VD4IHBeXApO1At+r8r7H8VLGIUaLSTJOIXq9AoYFVaUsbVUi\n0rq17+FsNsaOiWkNQSgkctxxIk2binz/vWe+0dxDjUKUfEpnVsgO8uVero5LgfEr8JPplpkWvDwf\nvqAvB7KAXHbXOC9iYo0kinHjquPA2L2hCgu1J40Tl+ZO4YCNs+Hee6FlS1/P8XJpbKixY8KxoXzF\nV1EKHnlE+19ef71nvqNGeXsLhULGhTPZPP/nFRRQzlccXOtaIl1oG5zAHzXKXbB8QV/y2E1PFtU4\nn6jgWdnO7NkwaZL2Dw8L47C/+HPPwfDhtQOG7VewnnsCN8Axx8A55/h+llfHbn5Pi65d4brr4Omn\nYeZM12QlJdCmjXs2RiFKPrlffwXAfA6qdS2RCkyDE/hukQjz8mBBTl8A+vJFjXuU0tqo8T2uG5Mn\nazdiJ8rLteCfPh2GDdO+3sOGwWcDR1NQuR0efhhPx/EIvDr2+uiTHw++/OZvugk6doQrroDKSte8\nBg50f/2mA00+/Qu/opIAi6jtg5/Q9+/H7pOqIxl++HZ757hbK2Ube8kD/MnV9lsfvTwyhVh2chIR\nkQ8+0BduuSWu52ViMLRUEZOH0osv6kQPPuiaX7ZNgmcaG4/+rSxWxXG/f7J10jYaWw86Qr7e+ygp\nLnYPs2sqeHzEspOT7Nqlg9906SKyY0fcz8ykYGipImbhHJ7Abd5cu2yK88R6NnegaadLF5nf84y4\n378R+G5cfrlI48Yy4szKrPPySDYxCaLwIpO33kpbeesrcXkoLVyod4W/+GLP0UE2dqBpZ9s2/UPc\ndlvc798IfDemTBEBObXn0tjMDwZf+NISf/hBpEkTkZNPTls56zMxm87CXHuthJSSo/I+NSPbTKK0\nVP8Ar70WdxZ+BX6Dm7SNSl89cTuo+RfGyyMJTJhQe2J2+nR9voqbbiK0s5yrud8E6oqDuD2Uxo5l\nS15bJu66Cqjt02oiaKaJr7SHDgfXdslMOH56hVQdKdHwy8tFcnJk7cgbzSRVOpgzRwTk3pzr631I\nhHRRlwnWCUV6hPt7Xoh5ZOtk9zckgCuuEGncWKSyMu4sMBq+C3l50KsX+/4419eG04YEEgqx7YIr\nWUc7xlWMQSJ89SdNMpq+H/xulu7E4pLz+JKDuZsbyKO8xjWv0YHbgjoTXjwBfPUVHHRQ3WNl+8FP\nr5CqIyUavojIeefpZfyhkJmkSiVPPCECci5PmMnyBOBWd7008dJSkZPz9G5io/mLr9GBcdlMIqGQ\n3rfykkvqlA1m0taDBx7QX/3bb1PzvCyntFTkvNO3yoacdvJZsL8oKs1keZJwCg2ulA6xHBbMY8eK\nvB04WbbQVPZmfVSTWjbGLUoZq1frF/nQQ3XKxq/Azz6TDsBhh+m/n32W3nJkAWFTQI+X72DvinVc\nVvkAbtXOTJbXDafQFqD/X7gQBg3Sv8eECdD+mXtorLbzYOsJzhPrNrIxblHKmDtX/03R/pxJF/hK\nqdVKqflKqXlKqbJkP88XffpATg58+mm6S9KgCQugfXYs4yru53HO4zMOd02fLSERkoVXaAvQXjjh\neZI+w3sQvOxShm16hKljv/a0+5u4RUlk7ly9b2cqPHRIXSydQSLSR0QOTdHzvMnPh969jcBPMmEB\n9BeuZxd53MRdjunMZHli8LMpeg3Xy3HjoHFjHWDNAxO3KDnMng1zH/uC5bndGXFBQUocFrLTpANw\n+OHapBMKpbskDZalS+Eo+ZjTeI2J3MiPtKtxPT/fw1ffEDN+NkWvYYJp3RrGjIG33oL333e9py5e\nQQZnwqbOdt/P5ZPyQ1Ln9eTH0F+XA1gFzAU+By5xuH4JUAaUFRUV1WniIiYef1xPlixalLpnZhkj\nzqyUORwm37KvFLC91oSf2VgjsfjZFL1WTKPycpFOnUT69o3qB2482hJD+Hdqx/ciIKO4v85eT2SK\nlw6wr/W3DfAlcLRb2pR56Yjo2CKgXQUNSWHp+GdEQEbyb0fhY3YaSzxjx4rk5bkL/Lw8h03fn35a\nX3z22XQXPysIez0N5S0RkF/xUY3faPDg2PP0K/CTbtIRke+sv+uBV8Fj1i6VdO8OTZoYO34UfMVc\nd6K8nAMeu4n5OX15GueNTUSMh0eimTABZsyAXr1qmncCAe2nEArBtGk1F0+NXzJCz2mNGQO7d7vm\nbUgM4bmW8L4c8+hT4/oHHyTPtJNUga+U2ksp1ST8PzAEWJDMZ/omENDumUbgu1Kn1ZUPPABr1vD6\nryaBcnfDbNkyzg7F4EpJCSxYALNm6V3G+vWDY4/VAn/PnuqJ3fAK54l/CXBzYCKsXMnjA/7J7Nl1\n6OgNUQnPtRzCXJbRlW00rXE9FEriqnM/w4B4D6AL2ozzJbAQuMUrfUpNOiIiN94okpsrs2fsNDFC\nIqjT6sr16/Xm2aec4plPTo5Ifr7PTTwMdSLaRuUQkg84RtbRRpoHt0purvldkkW4TaykkzzHsIQs\naCNTbPixHCkX+K+8IgJydN5sU7kjqNPqyiuu0DOy1oS4U8jkvDyR3Nw4OxRDzEQLqQwih6ED241l\nvPldksyfr9skAnIDdyVk1blfgZ+9bpnA50E9ndBn12wkYpib7YG84l5duWQJPPIIXHIJ9OgBOIdM\nPvpoqKhwzsKE6U08flw2P+NwXuT3jGYSrVlf67r5XRLHTSdq+/0X9HW8nqwFbVkt8CdN3ZdvKOII\nPql1Ldsrd9yrK2+4AQoKYPz4GqdLSmDqVCgr0383bjTL9VOJ1+IpO7dwJwXsZAx31LpmfpcEYs0d\nLipwXouarAVtWS3wly6FWRzJUcyEiA0hsr1yx7W68qOP4PXX4cYboU0bz/zNcv3U07lzdC1/Gd2Y\nwkVcyiN0ZmWt6y1bVv9vJnbrwKefQteuXHBdq9QuaPNj90nVkWob/vDhIperv4uA7Mcq7wUqWUhM\nm1pXVooceqiUt+kgI3+/I+oEuAm5mzqcImh6He35TrZTIE9zVq1r+fk6P699cQ01cQxXvc8+Imed\nVeN6XRa0YSZto1NaKtI/f54IyFk8bYSOA74ro7V458JGT/oWAjF1KIa48LP6NvIIBkX+rG4WAenD\nXMfFW/n5zvcqpRd2mbajceoY989fqz9Mnpyw5xiBH4WwINunbYVsoak8zKVG6MTLjh1S3q5I5qpD\nHGPde3WeZrl+conujln7GDxY5NRjtshPtJS3OTGme+2/eba3oUcf1Z1n5Lv5Ldo78E+HlSasvhuB\n70Fkr/suQ+RLDpL27Y3Q8cJ1J6WJE0VABvGBY+M3m2SkDz/umE7ugIccInINk0RABvJh3EI/W9vS\n2LHOwh5E7uIG2UWu5LEzYR2jEfguOA1xx3CbVKKkfcHmrKugfjemdrPZTrxWL7L6qOkpUYWIIfUM\nH15tMvNzhOeuhg8XKVQ7ZA0d5BNKBEIxC/xs7eijmdGmM0g+5dCEdoxG4LvgNMQdxHQRkJPU21lV\nQf1OvHlV4AdzrpRQICAXDlgYVYgYUk+sNvyw8AnfdwFTREB+w2tVabxs+Kajj7JokQrZSmP5G5cn\ntGP0K/Czzi3TaUHRHPpTQZABMitrXDGdtsMLLzq7/XYYMqTazc5tJ6X9Wc7FFQ/xYrOLeHpuT9dn\nmU0y0odbLPucHMjNdXcHDN/3YsG5LKaYP3MzOaqSwkK91OL663V6L7LVvdZr0WJ3FtOEX/jUFkMy\nlS7gWSfwnfy/d7AXX9CXX/G/rKmgXtvhiej9MMKB0twq8F3cxC7yuGrLBHbtcs4rP99skpFunFY6\n/+9/8PHHNc9FbkIzYQL894Mc3jnyTnqxiAcOf5rp02HoUF0nOnaEtm1BKefnZmtH77XG5HD0gqs5\n9K86l9KO0c8wIFVHumz4IHIP10o5jWT2hzuSXoZMwO9kXmGhdrOLtAP3p1QEZBzjPO+PJ7a3IcMI\nhUQOPVSkqEgm3FxeywyYk6PjIhn3Wo2XGe1h/iibaVbDmy2VNvys0/Ddhrif5B1LHrvpX1k7zEJD\npFUrf+nKy3W1bNTIflaYxGh+oB2TGO15/+efm5WY9R6lYOJEWLOGbfc8UssMWFGh29Bxx5ktK8FZ\nxoQ5klmUMoCw6FVKr4BOGX56hVQd6fDDD/t/z5m2VULBoLzS85YGHyY52q5IkUe7dlqLC38O+xFf\nzD9i8tjIZq2vIfBV28Gynr2lMVuNR44PSkv1CDc88mmOjpB5M3ckvG1gvHRiY+xYkdLAAJnFgAYt\noGL12lCqpj9xDrtlMd1kIT0kyB7f+SRy+GpID+cUfyqCe/jkbPTIiYbdY+ck3hQBOZoZCW8bfgV+\n1pl0nAh7rEwPDeJwPqUx2xpsmGSvyVonRKCysvrzxTxKMUu5gbupJCfm52d7FNL6TEXfw3jJCp+8\nNxtqXMtWj5xo2B0ejmImu8nlMw5zTJuKtmEEPtVC8AOOJYdKK3qmpqEJKC+XsWg0YSvjGc9HHM2b\nnOKaLhh0zyPbo5DWZ0aNgjvy76CAndzCnTWuZatHTjTsHjtHMZPP6cdOnP1ZU9E2ki7wlVInKqWW\nKKWWK6VuTPbz4iEsBD/hCHbRiGP5oOpaQxNQXi5jSmk3O7fr13EPbdjAaCYBzr54gQDsu68JfdwQ\nKSmBU68v5qng+fwfD7Mfq1FKd/DNmmnFqCGNhhNBOMx4HuUczqfM5CjXtKloG8nexDwIPAgMBXoC\nI5RS7it00kRYCJZTwCccUUPgNzQB5RXnvqAADj7YeQSwD99xLfcyleGUuQxJQed9661xxNI31Asm\nTIDer4xDBRR3548nENBK0Q8/xLjJfZYQ9tg5Mu9z8tjtKfBT0TaSreEfDiwXkZUisht4Djg1yc+M\nGbsQ/JBB9OULWrAJaHgCys0tNbzK0m0nqtsYSw4V3MKdUVdqXnSR9zPMIqz6zSG/6cBPw//EGeVP\n0r1yQVV9Cc973X23Xqlt3HE1EybA4xdqM/H2g4+gVy8tV9LSNvzM7MZ7AL8Hptg+jwT+HpHmEqAM\nKCsqKopvijoBhOPKDFB6QdFw9VyD9NIJ4xaW2CnY1oF8JRUE5F6ulpYtq9NHC21sQh83XC467SfZ\nQlN5lVONO64fTjlFpLi46mOi2waZ4JbpR+Dbj3S6ZYrolz5iWIVsDraSGZ3OzUoB5eS2+RZDZRPN\nZd+CjVn5Tgy1OeQQkZu5QwSkhE+MO64XlZUiLVqIXHhh0h7hV+An26TzHdDR9rmDdS4jKSmBZ58P\n0nzYEAbufJeSw0PpLlLKiTT5HMt0TuIdJuXezIXXtTTmGAOg57X+pkaxjrZM5EaI2BM6kobm7RaJ\n5/6+X34JmzfDMcekq3jV+OkV4j2AHGAl0BloBHwJ9HJLn24Nv4onn9RqSVmZ73jxDY3waGdxQR9Z\nX1gks2fsdE2Xje8n2wmPBC9D7wl9Au9E1fIzpXknmqhhxifpjWTku++SVgYywaSjy8FJwFJgBXCL\nV9qMEfg//igCMu2Y27N7o+Z//Ut/8alTHS+bjayzm7FjRZoV7JIVdJa59HHc3tK+Yrsh7ongtXK9\nyow1dKhI9+5JLUfGCPxYjowR+CKyrceh8kngiOy1R27dKtK2rciAATpaYgS+KrqhwVNaKvK3Er2B\n/ZlMdRX4wWDDrBOem50ERM4etltkr71ELrssqeXwK/DNSlsX3gsM5fDQ7Cr3TDv1zR7paV904667\n4Mcf4f77HQOee4VoqG/vxxA/JSVwxawRLC04mDsYQw57HNO1bt0w3XG9Vq6HQtDoy89g+3a9QCED\nMALfhf/sGUqQEEP4b61r9Wn17bhxuq49/zzMnetzcczq1XDffXDOOdC/v2OSaBW9vrwfQwIIBHjt\nsD/TlRVcyL9qXVYKBg1KQ7lSgNfK9UAAhjaarl/AwIGpLZgLRuC7sLvP4aynNafyeq1r9WX1rdc2\nhp5B4W64QX/Ju+5yzTtaRa8P78eQOI6eeBKzAkcxjgkUsKPGtYKChrV40Y7XyvX8fBiS+wH06eN/\nA4okYwS+C1deHeTN4G85hTfJo7zGtfqy+jYus8usWXoYcP310KGDa97RKnp9eD+GxFEyQLH4D3fR\nnnWMUg8ANVeQQhxmxXqA18r1G0ftpNmCTzLGnAOYSVsvnjzrHRGQU9Sb9dILJdo2hrVed2Wl3spu\nn31Efvklav5hLx2ztZ0hzKYjT5Ftuc3lmN6bqtx0s8Gby3Hl7Pvv6y/89ttJfz7GSycB7Nole/Zq\nKh90vqBehgdwCpNg9yAIu8mFK+uY/fT6g2Vjn/TM1+57P3iw3vO2Pr4fQxL48kst2W+4QUSy3Jvr\nmmtEGjXypTzVFSPwE8VZZ4m0aiWyZ0+6SxIzfhpbWPtqzDZZyz4yh8Nkr4JKV+0rG7Q1Q01iXlx3\nzjki+fkia9dGdVts0Nsidu8uMmRISh5lBH6ieOkl/Zo+/DDdJYkLL7OLvUO4ixtEQAYwy1X7ympt\nLUuJq4NfsUIkN1fkj3+M3azYUFi1Sn/B++9PyeOMwK8jYa3miN6/yK5gvvzw+yvSXSRH/GhfXpEx\nlRIp5mvZTY48xnme2ldWa2tZSJ06+CuuEAkG5aqTl/oyKzY4HnpIf8nFi1PyOCPw60CkVvMKp8kP\ntJNxYyrSXbQa1NW8orWvkPyXwbKZZtKaHz21r6zV1rKUOnXw69aJ7LWXbBh8ZnaOCk85RaRLF8dV\n6snAr8A3bpkROPmuP83ZtGMdn9/zQca4k8XtY2+jWzc4Q73M8UxjDHewgTY1rq9cWdOFzvjeZxd1\nWlzXti1ccw17T3uee8+em12b4ZSXwwcfwNChjqvU04qfXiFVRyZo+E5aTR47ZTPN5N/8IWPMFokw\nr8yZvk2+VR3kC3pLkD2ueTnZ/LNKW8tS/Hp5ubJli3Z4OOEEXxt+NJjIq+++q1/Sm2+m7JEYk058\nuJkt/slFspXGcmSf5LtY+SEh5pUbbxQBOTZvpmvDdvPqMb73DZ+EdPD33qtv+OADz2QNyvvr0kt1\nwLQdO1L2SCPw48RNqzmaGSIgDwx4NiM0kTprX4sXa0+Kc8+t+j4tWrgLfPuowWxdmD3UuYPfuVOk\nQweR/v1d7dkNauRYWSnSvr3I6aen9LFG4MdJaal2IY6seIpK+YaO8n7uiRIM+tNEktkx1KmRVFaK\nDBwo0ry5nlyzMJOyBifq2sEvv1nvq3BNl1cd729Q3l+ffKIL/vTTKX2sEfhxMnasSE6Oc+W7jTFS\niZIiVrsK2XDjaNdOfHcMdSlrXNrXP/+pb5gypcbpOo8aDIYIxo4VaVKwRxbRXb6mWBqp3bXqaINS\nNK67TguQzZtT+lgj8OPAS2sGkSJWSyVKJnCro0Ds1cv7/roMUd1GCzFrX999J9KsmcigQbWG2A1q\naG1IO/b6dApviIBcyV9r1ScvRQNEiorqSd0LhUS6dk3Z6lo7aRf4wHj0huXzrOOkaPekW+B7DS3D\nx5ucJGvZx9GrJdrEZ7xD1IROaP3ud9pmtXSp57PMpKyhrtRsTyF5lyGymWayN+trzQl5KUpK1ZM6\nOH++LvDDD6f80Zki8EfHck+6BX60oSWI/JrXRUBO5dWoaRMxRE2o1v3KKyIgz/SeGNfKXIMhFiLb\nU3cWyR6C8jB/rNUOIpWaejnKHDNGa0g//JDyRxuBHwfRhpYgEmSPfMu+8i5D4hL2sdrCo406fA93\nN22Snxu3ly9Vb8llt9HeDUnHqT3dzyipRElf9UWtdlBaqutzIkfHKSMU0itrjz8+LY/PFIG/GvgK\neAxo4ZLuEqAMKCsqKkruW4lCtKFl+LiZO0RAejE/ZoHvV0sJa9kFBf7yjCa0Nxw/QnaTI335vP5p\nToZ6iVN7as4m2UAr+ThwtJR+UttNs75N4Ibb6R+KZ4uALB/zeFrKkRKBD0wDFjgcpwJtgSB6V607\ngcei5ZduDV/E2YYdDNasdC35SX6hsEawsWhHLHZIP8PbmIT2c8+JgIzh9vqnORnqNU7t6YrcR/SH\n55+vlT4TPcXcHCbs7fSvXCk7yZN2BVvSMmJOu4Zf4yHQCVgQLV0mCHyR2jbs44+vLXwf4ArZRa7s\nw1pPIR8M6nUYfpeTP/qov1GGb6G9dq1IixbyVWF/1/AJmag5GRoOteaEZlbIL10Plg2FHeWIPttr\neZ1lkqeYm8PEBRdUlzPIHvmBtvICv09bOdMu8IH2tv+vBp6Ldk+mCPxInIaZnVkhFQRkItc7Vs6m\nTaNPeDpVpsjRRCxHrdcXComccIJIQYFcddKSjNOcDNnJ2LEix+d9JAJyJzfVmkvKFE8xr87H3paG\noGPnnMbLaRsxZ4LAfwqYb9nw37B3AG5Hpgp8t2HmVM6UrTSWvVlf61p+vncF9TtfEIuGX0toT56s\nLz74YMZpTobsxF4PH+dc2U2O9GRB1Yj4+ONrLmBMp6eYHzdtEHmB38sGWkkjytM2Yk67wI/nyFSB\n7yYsu7NIKgjIPVwbsyD1W5nsR9hE5OtZc+boWDm//nXVAqtM0ZwM2YGT7dte71uxQX6ipXzMUaKo\nrFGXM6FO+nHTbs2PsotcmcQ1aR0xG4GfYCKFZfh4gj/IDvKlPd85at1uQzs/lclJqIdth55Ce+NG\nkf3208fGjTWemwmak6Hh42b7bt++Zp0+Hx1n53z+lbZRp9ukrB837dH8RQSkO4uqzgWDIsXFqW1f\nRuAnAaeokp1ZIbvJkQf5P8cK4faVolWmYNBdqHsK7VBINh71a9kTyJWRxXOMUDekHC/zoT2+lD5C\n8hG/kp9oWcM0mio7eGTHZHe0GDzYOZCiveyL6Sb/40jH66kcQRuBn0QiK/TfuFwqCMjBzKv1g7sN\n7aLZ1B99ND5NfNogvUbgSiYbs40hLXiZK53Mkj1YKLvIlSf4gy9lKVH4CemQk6Mto07K2UA+FAH5\nA094jgLy8pKvdBmBn2TsmkELNsp69rZ6+pDvYWmibepL7nxRBOQpzq5RjlQPkQ3ZTTRzZbt2tQXt\nbYwRATmJN6MqS3XBbr4pKvI3j5aXpyeT+/XTWn+jRvr8q5wqG2glBWyPmkevXon/LnaMwE8B9lDI\nFyptixzJv2MS3AmzqX/2mZQHC2QWAySPnY7DS7O4ypAK/CyeKi3VwjOcrhHl8iUHyXe0l+ZsSoqC\nEs+Cxsi2Ew6fvj/LpBIltzHGdx7JVLiMwE8xpbMqZUmrEtkcbCWXnvp9arXpVatE9tlHvm9UJG1Y\n51rp6vHrNdQjYnEBto9y+/K57CZHngmOrDVfVddNhOrqBt2vX808HuAKKaeRtOWHKoEeLY9kKlxG\n4KeDr78WKSiQTQOGyvAzQzFX0rgq9/ffi+y/v0jz5jL6xPmuFU8pPYSt9xtEG+oFsZgr7aPcF3uN\n0ze89lpVHvZ6HAxqTyhjVO4AABibSURBVLVYiccNOnJUEs6jOZvkFwrlcc6t0b78dBrJwgj8NPHm\n0L+LgFzF/VErup24Yt7/9JPIgQfqDZNLS31NQsVSJoOhLsRlrty1S6RPH9nVsq0UFdRe0Bg+vIS+\nk+IUjxu0vZN59NHqPMYyXgRqOWl4efQk2zffCPw0UFoqUlgQklf4rewhKMcyzXUoW+u+WFfBfved\nngnKyxOZNq3qdKRmVa/jixuyk6++kl2BPHmLoTUWZEUKYae666Y49erlz+zipiiF82ipNskWmsrL\nnFZLoB9/vLvQT3ZbMwI/DYSHfI3ZKvPpJT/RsmrZuNekabTh5uDBETesWKFjbzduLDJ9uojU1GoG\nD672KvDyRDATuYZM5a6OD4qAXMs9ru0isu56KU55edF86qMfeXkidwVvEQE5iC8dBXq6VrMbgZ8G\n7MPGzqyQ72gv39NOurJUwN2GF224GQjYKsyHH4q0bi3SsqUOnyDe5qD6Fl/cYBARGX5mSF7kdNlN\njhzObF9110tx8rvntNfRWm2QXwKN5cXgME+Bno7V7Ebgp4FId7QeLJT17C0/0Fb6qzmuNjw/S7gb\nF1TI6ivu0WPZ7t1FFi8WkejmoOOPz7z44obsxa9jQmmpSMvAZlnFfrKSTtKcTbXMLJF1149yE+kO\nGnl4jbTDCyznPbsw48KTGIGfBpyEb3cWyQo6yw7yZflNj1YFMot2n/3oxmL5H0fqD7/7ncjPP1fd\nG02rGTzYRMk0ZAaxOiZccIHI4cyWXeTKfxlcYz8Hp7oby+YpTqaXvDwd1tzp/gP5SioIyHsHXJa8\nF1QHjMBPE04VqahgvazodKw+MWSIyBdfON4XKbg7sEYe4RLZQ1A20Vxu7fRkrQ7Dj1ZjomQa0o2X\nUmMPixzJBRdUL2r8K1d67hwXq/OD3fTSq5eXjT8k0zhWNtJCPnv3p+S8oDpiBH4acbThVVaK/P3v\nIs2aVQv+KVO0aWbPHhEROem4cunGYrmIf8p/OFkqCMgucuUBrpB2ap2j+cVLq7H73tsncjNlGGrI\nHvz4wXsJ8re6XSUC8o/DHk14uJJoI+wLrIie/znpobq/iCRhBH6msnmzyIQJIp0715TM+fkSsrWI\nNXSQ27lFilhd1RjCAdXs2yEOHuzdkIzvvSET8OsH72pm3LNH794WDIr85z+ez/Izaeo3ps6+fCs/\nq6byc9+BWmnLUPwKfKXTZgaHHnqolJWVpbsYqUEE5s+HL76AFStg1y4oLOSVLzpxx3uH8+Wu7oRE\nEQhAfj4ceiiUlcHOnfrWeCkshOnToaQkcV/FYIjGiBHwwgsQCnmnCwRg2DCYOtXh4rZtcOyxsGAB\nvPceHH10XGUZNw4mTYrelgJU8l+GcERgNgXL5kOXLnE9LxUopT4XkUOjJvTTK7gdwBnAQiAEHBpx\n7SZgObAEOMFPflmh4fsgUkOJd2Nzt8kr43tvSDWxxLLxFAMbNsiO/brL9tymcm63T6J6+kR6BMVS\njju4WQTkkf6PJeWdJBJSYdIBegDFwAy7wAd6Al8CeUBnYAUQjJafEfjO1CUOSMwNymBIEk6xcZwU\nEi9X4bFjRQ7IXyNL6Srb2EsGq2mOpkqvFbd+2tLp6FDjjwcvrBfzXX4FfqAuwwgR+VpEljhcOhV4\nTkR2icgqS9M/vC7Pqs/Mnq2HtP366b+zZ8d2/9KldTPj2AkEoFu3xORlMMTChAnanDh4sK6HTuTn\nw5VXOl+bPVubYpaVd+RoPmYVnfmPnMzQHS8xaVJ1uwqn27Gjut2EQvrzokXR29IgPuAZzqY0cATf\nXv+3hmX+9NMrRDuoreH/HTjH9vlfwO9d7r0EKAPKioqKktgHpoe4gqJF4Gdhlt/D+N4bMoF4vGki\nR7ot+Uk+oUQE5HbGyIgzKx3T+T2UEjm79XuyPbCXrGnWSz59b1OK3kbdIVEmHWAasMDhONWWJm6B\nbz8amkknrqBoMeYTS2U2XjqGTCLWEAROnj6NKJcpXCAC8mmTQSKrVvkKVVL7fEgubfQvqQzmiBx8\nsA477lLeTAwxnjCB7yuT2gL/JuAm2+f3gAHR8mloAj+aplFU5L/SRGpEdkEersSR+2/aN2TOtApq\nMMSK+0g3JBepKbI9p4lUFOwldza/23XbwUBARxS3t6U2ar08HxyhPxx7rMiWLbWenYiRejJJt8Dv\nRc1J25Vk4aStH9/jWCqNk/dOpIaUjsBNBkMqiDZiHj3sG3k7cLIIyHe0l5u4U/blW8eRdWmpyBW/\nXi3/bH+rbM9pIqFgUOSOO0QqKmJ+bia0sZQIfOA0YC2wC/gReM927Ra0d84SYKif/BqawPdre8+U\nSmMwZDputv8LLqgWykfxsUzj2KoGNp9e8gR/kL/k3CQfH3mjyIUXivTpU90ATztN71bnQrR4VZng\n5uxX4JuFV0lk9mw47jjtHeCF52ITgyGLmD0bJk/WnmndusGoUbUXCYbTLFsGBxyg00yeDM8/r8Vw\nmC6sYBgvcAwzOCi4iHasI6CAFi2gTx+9cOvss6FzZ88y9esHc+d6X0+32ErJwqtEHw1Nwxfx53sM\n7r7xmTxRZDAkkljt5Pa20aKFj/blEKnWD7FE4UwXmFg6mUNpqZ6g9fIacKo0mT5RZDAkiljt5JFt\nI5qHWjiIYDxKk7HhG4EfM/GEbs30SmYwJAo/dvKwRl9crL3Pogn6SKFfF6Up00OM+xX4dVppa/BP\nSQmMHq2Dl4VXGQYC+vPo0bXtlJMn6+BOTpSX6+sGQ0PBazV5KAQzZuj5sOefhyVLoLIyep5KVf8f\nzju84ta+MtcP4VXCw4Zpm/2wYfrzhAn+88gEctJdgGxiwgQYOrT2hJPT0u1oDWDZsuSW1WBIJd26\nwbx5ztE0lYING/wJedBzsl266Hu+/da5HYWVpljCJpSU1P8os0bgpxi/lcarAZh4OIaGxqhR8MYb\nzh5tgUD0sMr2tCeeCM8+qzXxNWuc02Wr0mRMOinEHkTt+ONhyBD3gGqjRulAUk54BZgyGOojXibP\nNm3cR7uR2NtGt27uQdqyVmnyY+hP1dGQJ229vArcJoAyfaLIYEg0TivF/SxgdGob2eT4gFl4lTn4\nXYDltBuV0yKT+m5HNBhiwav9BIO6XfTp49w2wrtblZdrM054B7nRo+vfhKsXfhdeGYGfAkaMqL0K\n0Amz4tZgcKYugrsuSpOflb+ZgBH4GUS0pdl2WrSAt9/OzEplMKSTVI92I/e+zeTRgRH4GYTfDZzD\nhH3zM61SGQzZgpcZycn0mm78CnzjpZMCvDxunIhnYYjBYHAn1m1GG+rCR+OHnwLCLmd2G2Q04lkY\nUh/Ys2cPa9eupby8PN1FyRry8/Pp0KEDubm56S5KWog0zcybp33+vUbRDXXhoxH4KSJylW2LFjBn\nDmzb5py+PlcqL9auXUuTJk3o1KkTyr723ZAURISNGzeydu1aOkcJA9wQsW9oHsYeXmHoUGelqqEu\nfDQmnRRSUqI9cMrK4P334eSTs29hSHl5Oa1atTLCPkUopWjVqlWDHVFFM9XEa5ppqAsfjcBPIw21\nUkXDCPvU0lDf97hx1QHV5s6F556DI46AAw+sFvzxmmZiDXZYX6iTwFdKnaGUWqiUCimlDrWd76SU\n2qmUmmcdj9S9qA2PhlqpEkmsk22G7MBuqrELdBFYuBAGDdIdQl3CKzSUCJk18LMc1+0AegDF1N7E\nvBOwINb8GnJoBS+yaePxRYsW+U6brA1gAoGA9O7dW3r27CkHH3ywTJo0SSorKz3vWbVqlTzzzDN1\ne7APLrzwQlm4cKFnmldffTVqmkhiee/1Aa/4+fbwCY8+mh3hFUhFPHwR+VpEltSpxzHUsO1PnWo0\ne3DW4OKNZR5JQUEB8+bNY+HChbz//vu88847TIiitq1evZpnn302/of6ZMqUKfTs2dMzzWuvvcai\nRYuSXpZMxstUE6a8XGvkZhRtw0+vEO3AWcPfDnwBfAT8yuPeS4AyoKyoqCi53aAh7fjVNP3sgBQv\ne+21V43PK1askJYtW0ooFJJVq1bJUUcdJX379pW+ffvKrFmzRESkf//+0rRpU+ndu7fcd999runs\nrFq1SoqLi+Wss86S7t27y+mnny7bt28XEZFp06ZJnz595MADD5Tzzz9fysvLRURk4MCB8tlnn1WV\n8+abb5aDDz5Y+vfvL+vWrZNZs2ZJixYtpFOnTtK7d29Zvny5TJ48WXr06CEHHXSQnHnmmY7fuSFq\n+NECqlXtZSsNfxRNorY4BKYBCxyOU21pIgV+HtDK+r8f8C3QNNqzstWkk034FTyHHOKvIcdDpMAX\nEWnWrJmsW7dOtm/fLjt37hQRkaVLl0q4Tn744Ydy8sknV6V3S2dn1apVAsjMmTNFROT888+Xe+65\nR3bu3CkdOnSQJUuWiIjIyJEj5f777xeRmgIfkDfeeENERK677jq5/fbbRUTk3HPPlRdffLHqOe3b\nt6/qMDZv3uz4nRuawPeKhBltr+iGiF+BH9WkIyKDReRAh+N1j3t2ichG6//PgRVAA3QyNCSLdMUy\n37NnDxdffDEHHXQQZ5xxhqvpxG+6jh07cuSRRwJwzjnnMHPmTJYsWULnzp3pZn2Jc889l48//rjW\nvY0aNeKUU04BoF+/fqxevdrxGQcffDBnn302Tz/9NDk52bG0JuzwkJfnnqYhe7rFS1LcMpVSrZVS\nQev/LsABwMpkPMvQMEmly+rKlSsJBoO0adOG+++/n7Zt2/Lll19SVlbG7t27He/xmy7SJTIWF8nc\n3Nyq9MFgkIqKCsd0b731Fpdffjlz587lsMMOc03X0JgwQe9126tXTeUgq230UairW+ZpSqm1wADg\nLaXUe9alo4GvlFLzgJeAS0VkU92KasgmUuWyumHDBi699FKuuOIKlFL8/PPPtG/fnkAgwFNPPUWl\ntZFqkyZN2GZbFu2WLpI1a9ZQWloKwLPPPstRRx1FcXExq1evZvny5QA89dRTDBw40HeZ7WUJhUJ8\n++23DBo0iLvvvpuff/6ZX375Ja53UR8pKYEFC2DWLBg+vAG5TyaJOo3/RORV4FWH8y8DL9clb4Mh\nlk3fY2Hnzp306dOHPXv2kJOTw8iRI7nmmmsAuOyyyzj99NN58sknOfHEE9lrr70AbTYJBoP07t2b\n8847zzVdJMXFxTz44INccMEF9OzZ8//bu/sYqaozjuPfXxCZxDWAlSp0tSA1BmzIYtGwRBpTfKEb\noqVZmuWf0ko0S2ssJE2DMVFTNVGaIiGtNbaVbRsDa2mlG4ORl27TRCN9weVFdoVdC2EREdaCRRIo\n7dM/7ll6mczszDJz7+zuPJ9kwuWeM3OfOffsM3fuvXMOy5YtI5PJsG7dOhYtWsT58+e59dZbaW5u\nLjr+pqYmHnjgAdauXcuGDRtYunQpp06dwsx4+OGHGTduXGkNNAyNhAnG0+DDI7tUdXZ2Mm3atEqH\nkYqDBw+yYMEC9u7dW+lQqqrdq5EPj+ycc+4invCdS8jkyZOHxNG9c/084TvnXJXwhO+cc1XCE75z\nzlUJT/jOOVclPOG7qlNTUwPABx98QGNjY1les6uri7q6OmbOnElPTw9z5swB0htl07lieMJ3VWvS\npEls3LixLK+1adMmGhsbeeedd5g6dSpvvfUW4AnfDS3VMdKSG5qWL49mii6nujpYs6aoqvEfRrW0\ntNDW1saZM2fo6elh4cKFrFq1CoAtW7bw+OOPc/bsWaZOncq6desufEsA2Lx5M2vWrGHUqFFs376d\n9vZ2ampqOH36NCtXrqSzs5O6ujqWLFnCihUryvt+nRsEP8J3Lujo6KC1tZU9e/bQ2trK4cOHOXHi\nBE899RTbtm1j586dzJo1i9WrV1/0vIaGBpqbm1mxYgXt7e0XlT3zzDPMnTuXjo4OT/au4vwI31VO\nkUfiaZk3bx5jx44FYPr06Rw6dIiTJ0+yb9++C0Mcnzt3jvr6+kqG6dwl84TvXDAmNrh6/3DEZsZd\nd93F+vXrKxiZc+Xhp3ScG8Ds2bN58803Lwxl/Omnn7J///6in589rLJzleQJ37kBTJgwgZaWFhYv\nXsyMGTOor6+nq6ur6OfHh1V+7rnnEozUucJ8eGSXKh+mtzK83Ue2VIZHlvQjSV2Sdkt6VdK4WNkj\nkrolvSfpnlK245xzrnSlntLZCnzRzGYA+4FHACRNB5qAm4H5wPP9c9w655yrjJISvpltMbP+GZPf\nBmrD8n3ABjM7a2b/ALqB20rZlhs5htJpxGrg7e36lfOi7f3A62H5c8DhWFlvWOeqXCaToa+vz5NQ\nSsyMvr4+MplMpUNxQ0DB+/AlbQOuzVH0qJn9IdR5FDgPvDzYACQ9CDwIcP311w/26W6Yqa2tpbe3\nl+PHj1c6lKqRyWSora0tXNGNeAUTvpndOVC5pG8BC4B59v/DtiPAdbFqtWFdrtd/EXgRort0Cofs\nhrPRo0czZcqUSofhXFUq9S6d+cAPgHvN7EysqA1okjRG0hTgRuAvpWzLOedcaUodWuEnwBhgqySA\nt82s2czelfQKsI/oVM93zew/JW7LOedcCUpK+Gb2hQHKngaeLuX1nXPOlc+Q+qWtpOPAoRJe4mrg\nRJnCKSePa3A8rsHxuAZnJMb1eTObUKjSkEr4pZL0t2J+Xpw2j2twPK7B8bgGp5rj8sHTnHOuSnjC\nd865KjHSEv6LlQ4gD49rcDyuwfG4Bqdq4xpR5/Cdc87lN9KO8J1zzuXhCd8556rEsEr4khZJelfS\nfyXNyiorOOGKpCmSdoR6rZIuTyjOVkkd4XFQUkeeegcl7Qn1Ep/qS9ITko7EYmvIU29+aMduSStT\niCvvRDpZ9RJvr0LvPQwX0hrKd0ianEQcObZ7naR2SfvC38D3ctS5Q9Kp2P59LKXYBtwviqwNbbZb\n0i0pxHRTrB06JH0iaXlWnVTaS9JLkj6StDe27ipJWyUdCP+Oz/PcJaHOAUlLSg7GzIbNA5gG3AT8\nCZgVWz8d2EU0zMMUoAcYleP5rwBNYfkFYFkKMf8YeCxP2UHg6hTb7wng+wXqjArtdwNweWjX6QnH\ndTdwWVh+Fni2Eu1VzHsHvgO8EJabgNaU9t1E4JawfCXRhEPZsd0BvJZWfyp2vwANREOnC5gN7Eg5\nvlHAh0Q/Tkq9vYAvA7cAe2PrVgErw/LKXH0euAp4P/w7PiyPLyWWYXWEb2adZvZejqKCE64oGuzn\nK8DGsOpXwNeSjDds8xvA+iS3U2a3Ad1m9r6ZnQM2ELVvYiz/RDppK+a930fUdyDqS/PCfk6UmR01\ns51h+V9AJ8Nnjon7gF9b5G1gnKSJKW5/HtBjZqX8iv+SmdmfgY+zVsf7Ub5cdA+w1cw+NrN/Es0w\nOL+UWIZVwh9AMROufAY4GUssaUzKMhc4ZmYH8pQbsEXS38O8AGl4KHytfinP18hKT14Tn0gnW9Lt\nVcx7v1An9KVTRH0rNeE00kxgR47iekm7JL0u6eaUQiq0Xyrdp5rIf9BVifYCuMbMjoblD4FrctQp\ne7uVOlpm2amICVeGgiLjXMzAR/e3m9kRSZ8lGnG0KxwNJBIX8DPgSaI/0CeJTjfdX8r2yhGXFT+R\nTtnba7iRVAP8DlhuZp9kFe8kOm1xOlyf2UQ0NHnShux+Cdfp7iXMt52lUu11ETMzSancHz/kEr4V\nmHAlj2ImXOkj+ip5WTgyyzspSzEKxSnpMuDrwJcGeI0j4d+PJL1KdEqhpD+UYttP0s+B13IUFT15\nTTnjUu6JdLJfo+ztlaWY995fpzfs47FEfStxkkYTJfuXzez32eXxDwAz2yzpeUlXm1miA4UVsV8S\n6VNF+iqw08yOZRdUqr2CY5ImmtnRcHrroxx1jhBdZ+hXS3T98pKNlFM6BSdcCUmkHWgMq5YASX5j\nuBPoMrPeXIWSrpB0Zf8y0YXLvbnqlkvWedOFebb3V+BGRXc0XU70dbgt4bjyTaQTr5NGexXz3tuI\n+g5EfemP+T6gyilcJ/gl0Glmq/PUubb/eoKk24j+vhP9MCpyv7QB3wx368wGTsVOZyQt77fsSrRX\nTLwf5ctFbwB3SxofTr/eHdZduqSvUJfzQZSkeoGzwDHgjVjZo0R3WLwHfDW2fjMwKSzfQPRB0A38\nFhiTYKwtQHPWuknA5lgsu8LjXaJTG0m332+APcDu0OEmZscV/t9AdBdIT0pxdROdq+wIjxey40qr\nvXK9d+CHRB9GAJnQd7pDX7oh6fYJ272d6FTc7lg7NQDN/f0MeCi0zS6ii99zUogr537JikvAT0Ob\n7iF2h13CsV1BlMDHxtal3l5EHzhHgX+H/LWU6LrPduAAsA24KtSdBfwi9tz7Q1/rBr5daiw+tIJz\nzlWJkXJKxznnXAGe8J1zrkp4wnfOuSrhCd8556qEJ3znnKsSnvCdc65KeMJ3zrkq8T9awpj7DI8B\n6wAAAABJRU5ErkJggg==\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": "markdown", "metadata": {}, "source": [ "(d) Solve for $a_0$, $a_1$, $a_2$, $a_3$, and $a_4$." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[ 24.00958042 -49.99515152 35.0039627 -9.99561772 0.99841492]\n" ] } ], "source": [ "x_data = [0.0, 0.5, 1, 1.5, 2, 2.5, 3, 3.5, 4, 4.5]\n", "y_data = [24, 6.61, 0, -0.95, 0.07, 0.73, -0.12, -0.83, -0.04, 6.42]\n", "\n", "D = data_matrix(x_data, 4)\n", "params = np.linalg.solve(np.dot(np.transpose(D), D), np.dot(np.transpose(D), y_data))\n", "\n", "print(params)" ] } ], "metadata": { "anaconda-cloud": {}, "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.3" } }, "nbformat": 4, "nbformat_minor": 1 }