{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "## Optimization with `gradpy`" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "At the core of therapy planning for IMRT is an optimization problem which answers the following question: What are the optimal parameters, i.e. beam intensities and collimator positions, for which a radiation dose can be delivered as close to the desired dose for a given region, respecting any additional constraints, such as maximum and minimum local doses? Such questions are commonly framed as minimizations of a cost or objective function which penalizes deviations from the target, as well as failure to satisfy imposed constraints.\n", "\n", "This framework for solving optimization problems is quite general, and has been addressed by a variety of numerical optimization methods. As typical cost functions may be nonlinear with respect to the parameters to optimize, we implement the [Broyden-Fletcher-Goldfarb-Shanno (BFGS) algorithm](https://en.wikipedia.org/wiki/Broyden%E2%80%93Fletcher%E2%80%93Goldfarb%E2%80%93Shanno_algorithm), an iterative method used to solve nonlinear optimization problems. The BFGS algorithm relies on an available gradient of the objective function, to be provided by automatic differentiation using the `gradpy` package, and builds an improved estimate of the Hessian matrix at each iteration using the current solution and gradient information.\n", "\n", "In the following example, we illustrate the use of our BFGS implementation with `gradpy` to find the optimal parameters of a polynomial fit to data. We start by importing the relevant modules and define a set of random $(x,y)$ data points to fit by a polynomial:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": true }, "outputs": [], "source": [ "from gradpy.autodiff import Var\n", "from therapy_planner.bfgs import BFGS\n", "import numpy as np\n", "\n", "data = np.array([[-2.9902373, 2.98118975],\n", " [-2.29029546, 10.58188901],\n", " [-1.64611399, 7.63332972],\n", " [-0.99102336, 5.20796568],\n", " [-0.34860237, 4.60892615],\n", " [ 0.36251216, 6.58063876],\n", " [ 0.98751744, 6.74456028],\n", " [ 1.74502127, 4.1763449 ],\n", " [ 2.42606589, 1.80673973],\n", " [ 2.9766883, 9.29002312]])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next we define functions that will comprise the overall cost to minimize. `MS_error` computes the mean squared error between the target $y$ and the polynomial fit `poly`. `L2_reg` is an L-2 regularization on the model parameters with weight alpha." ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": true }, "outputs": [], "source": [ "def poly(x,params):\n", " return params[0] + sum([params[i]*x**i for i in range(1,len(params))])\n", "\n", "def MS_error(x,y,params):\n", " return sum([(poly(xi,params) - yi)**2 for xi,yi in zip(x,y)])\n", "\n", "def L2_reg(alpha,params):\n", " return alpha*sum([p**2 for p in params])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We choose a polynomial order and define a list of the corresponding number of parameters needed. The cost function is defined as a sum of the mean squared error and regularization contributions. As the cost is a function of `Var()` objects in the `params` list, it is a differentiable `autodiff` object." ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "collapsed": true }, "outputs": [], "source": [ "order = 5\n", "params = [Var() for _ in range(order+1)]\n", "\n", "alpha = 0.1\n", "cost = MS_error(data[:,0],data[:,1],params) + L2_reg(alpha,params)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Using the implemented BFGS algorithm, we pass in the function to minimize (`cost`), parameters on which it depends (`params`), and an initial guess for each parameter. As `params` is a list of `Var` objects, the values are updated in place. The optimization returns the final step size and number of iterations performed, reporting `Minimum found` if the minimum is successfully obtained within the specified minimum step size tolerance." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Minimum found.\n" ] } ], "source": [ "step, Niter, found = BFGS(cost,params,np.ones(len(params)))\n", "fit_params = [p.value for p in params] # extract fitted parameters" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally, we visualize the data and the fit." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlgAAAFtCAYAAADbORRfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4xLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvDW2N/gAAIABJREFUeJzs3Xd4lFX6xvHvmZBAICSUUEKvIkhT\noqISFRs2QBcrNhRULLu6q+4uuvpbV1dddV1dGypYVnTtsqCIYkGJSgkgBEWKNElCCSWQAgmZ9/fH\nISASSEJm5ky5P9c11yTzTmbuoJl55rznPMd4noeIiIiIBI7PdQARERGRaKMCS0RERCTAVGCJiIiI\nBJgKLBEREZEAU4ElIiIiEmAqsEREREQCTAWWiIiISICpwBIREREJMBVYIiIiIgGmAktEREQkwOq4\nDpCamup16NDBdQwRERGRKs2dOzff87xmVd3PeYHVoUMHsrKyXMcQERERqZIxZnV17qdThCIiIiIB\npgJLREREJMBUYImIiIgEmAosERERkQBTgSUiIiISYM5XEYpI5PH7/eTn57N161bKy8tdx4k4cXFx\nNGrUiNTUVHw+fc4ViUYqsESkxtauXYsxhg4dOhAfH48xxnWkiOF5HmVlZaxfv561a9fSrl0715FE\nJAj00UlEaqyoqIjWrVuTkJCg4qqGjDEkJCTQunVrioqKXMcRkSBRgSUih0SntmpH/34i0U1/4SIi\nIiIBpgIrSvj9HhPn5zD4yUzS75/G4CczmTg/B7/fcx1NJGz4/X6uv/56mjZtumcO2bnnnus6lohE\nIU1yjwL+km2Mfv07MlcXU1xqV3TlF5Yy5r1spmTnMfbyfvh8micjMmXKFF566SWmT59Op06dSExM\nxPP2fgg5+eST6dmzJ0899ZTDlCISDTSCFYn85bDyK9i6BoBJn39J5rKNe4qrCiVl5cxYls/khbku\nUoqEneXLl5OWlsbxxx9Py5YtSUlJoVGjRq5jiUgUUoEVSUq2wMd3wb+OgFcGw7z/ADB+WQOKqVf5\nj5SVM27GylCmFAlLI0aM4Pe//z1r1qzZc3pwxIgRe04Rjhgxgi+//JKnn34aYwzGGFatWuU2tIhE\nLBVYkaJwI7w8GGaNhVZHwQUvwYA/AJBXePBGj3n5m6EgJxQpRcLWE088wT333EObNm3Iy8tjzpw5\n+x0/7rjjuPrqq8nLyyMvL4+2bds6SisikU5zsCLF14/DpuUw/C3ocuo+h9JSEskvLD3gj6aVrYFn\nboLT/gr9rgYtD5dgeOmc/W874jw45looLYbXLtz/eN/hcORlULQJ3rpy/+NHXwM9h0HBWnjv+n2P\nXf1hjeKlpKTQsGFD4uLiaNmyZaXHExISqF+/fqXHRURqQu+0keLUe2Dkx/sVVwAjB3QkMT6u0h9L\njI9j1NknQKsj4cM/wLsj7RwuERERCRqNYIWzTT/B1D/DeWOhQVNI61Pp3Yb0acWU7DxmLMunpGxv\n8ZQYH0dG11QGH98XTvgfZD4Gn/0NEhrAUK2SkgA72IhSQv2DH2/Q9ODHU9rUeMRKRMQlFVjhatNP\n8NLZ4C+DwnX2DegAfD7D2Mv7MXlhLuNmrCSvoIS0lERGZXRkcO9We1s0ZNwGngfNDg/RLyESWRIS\nErR5tYgEhAqscOR58L+boXwnXP0RNO9e5Y/4fIahfVsztG/rg9/xxNv3fr0uG1r0BO0lJwJAhw4d\nmD17NqtWrSIpKYkmTZpoSxsROSR65QhHC96ANd/A6X+rVnF1SHLmwnMnwef3BefxRSLQ7bffTkJC\nAj169KBZs2asWbPGdSQRiVDml12MXUhPT/eysrKcZggrnrf71OAuuObj4K348zyYfAvMe8UWcifc\nEpznkai0ePFiuncPUvEfQ/TvKBJ5jDFzPc9Lr+p+OkUYboyBKydC8abgtlMwBs59HHYUwKd/hbbH\nQrv+wXs+ERGRGKJThOFkyyrYWQh16kJyq+A/n88HQ56ERu3g3VG22BIREZFa0whWuPCXw9sjAAPX\nfh66ief1kuGCF2H1t5DQMDTPKSIiEuVUYIWLuS9D7nz4zbjQr+pr3c9eAHbttCNoIiIicsh0ijAc\nFOXbBqAdMqDXBe5yrPoanugDGxa7yyAiIhIFVGCFg68egdJCOPtRtz2pmnaxqxffuQbKStzlEBER\niXAqsFzzPNi8EnpeAM0dd1hv2MJuy7PhB5h2j9ssIiIiEaxaBZYx5kRjzCRjTI4xxjPGjPjVcWOM\n+asxJtcYU2KMmW6MOSIoiaONMXDZW3Y1XzjoehocOxpmvwBr1Z9MRETkUFR3BCsJWATcAlR27uiP\nwG3Ab4GjgQ3ANGOMlqUdjOdB0Sb7dZ0Et1l+aeBdkNQClnzkOomIiEhEqlaB5XneFM/z7vQ87x3A\n/8tjxhgD3Ao85Hneu57nLQKuAhoCwwMdOKr8PAv+2Q1WTHedZF/1kuGGr+HUu10nEQm6c889lxEj\nRriOISJRJhBzsDoCLYFPKm7wPK8E+Ao4PgCPH73mjIP4RGhztOsk+2uQaq83LoXCjW6ziISJ6dOn\nY4whPz/fdRQRCXOBKLBa7r5e/6vb1//i2D6MMdcZY7KMMVkbN8bom3fhRvh+IvQdDgkNXKepXMkW\neP5kmKaRLBERkZoI5CrCX+8abSq5zd7R8573PC/d87z0Zs2aBTBCBJn/H/CXQfpI10kOLLExHHs9\nLPgvrP7GdRqRWisuLmbEiBEkJSXRokULHnjggX2OT5gwgaOPPpqGDRvSvHlzLrzwQnJycgBYtWoV\nAwcOBKBZs2YYY/acWpw6dSoZGRk0btyYJk2aMGjQIBYvVj85kVgWiAJr3e7rX49WNWf/US0BO7l9\n7ivQ8URodpjrNAd34u2Q0hY+vA3Ky1ynkSjh93tMnJ/D4CczSb9/GoOfzGTi/Bz8/ko/kwXM7bff\nzrRp03j33Xf57LPPmD9/Pl999dWe46Wlpdx7770sWLCADz74gPz8fC699FIA2rZty7vvvgvA999/\nT15eHk888QQARUVF3HrrrcyePZvp06eTkpLC4MGDKS0tDervIyLhKxBb5azEFlmnA3MAjDH1gAzg\njgA8fvQxBq78n20uGu4SGsCZD8Gbl8Gs5+D4m10nkgjn93uMnjCXzOX5FJeWA5BfWMqY97KZkp3H\n2Mv74fMFvuFuYWEh48eP58UXX2TQoEEAvPTSS7Rp02bPfa655po9X3fq1Ilnn32W7t27s3btWtq0\naUOTJk0AaN68OampqXvuO2zYsH2e66WXXiI5OZnZs2czYMCAgP8uIhL+qtsHK8kY09cY03f3z7Tb\n/X07z/M84HHgz8aY3xhjegIvA4XA68EKHvGadISWvVynqJ7Dz4Hug2HHVtdJJApMWpC7T3FVoaSs\nnBnL8pm8MDcoz/vTTz9RWlrKcccdt+e2pKQkevXa+3c4b948hg4dSvv27WnYsCHp6ekArFmzpsrH\nHj58OJ07dyY5OZkWLVrg9/ur/DkRiV7VPUWYDszffUkE7t399d92H38YeAx4GsgC0oAzPM/bHtC0\n0WDzCnj9Eshf7jpJ9RkDF70Kp/zFdRKJAuMzV+5XXFUoKStn3IyVQXle+1nwwIqKihg0aBD169fn\n1VdfZc6cOUydOhWgylN9gwcPZuPGjTz33HPMmjWL+fPnU6dOHZ0iFIlh1e2DNd3zPFPJZcTu457n\neX/1PC/N87x6nuedtLsflvza3Jdh2Sfhu3LwQCr2SFzxJWz92W0WiWh5BQff57Kq44eqS5cuxMfH\nM3PmzD23FRUVsWiRfan68ccfyc/P54EHHuDEE0/k8MMPZ8OGDfs8RkKCbQhcXr63QNy0aROLFy/m\nzjvv5LTTTqN79+5s376dXbt2BeX3EJHIoL0IQ8nzYNH70PV0SE5znabmCjfCaxfCFw9UfV+RA0hL\nSazV8UOVlJTEyJEj+dOf/sS0adP4/vvvueaaa/YUS+3ataNu3bo89dRTrFixgg8//JC77963RUn7\n9u0xxvDhhx+yceNGCgsLady4MampqbzwwgssX76cL7/8ktGjR1OnTiCmuIpItfw8G94fDduCM8Xg\nUKjACqX130PBGjunKRIlNYNjrrVtG9b/4DqNRKiRAzqSGB9X6bHE+DhGZXQM2nM/+uijDBw4kPPP\nP5+BAwfSs2dPTjzxRMC2XnjllVeYOHEiPXr04N577+Wxxx7b5+dbt27Nvffey1133UWLFi24+eab\n8fl8vPnmmyxcuJCePXty0003cd9991G3bt2g/R4i8iuNO0D3IVA3fHboM1XNSwi29PR0LysrRjYV\n/vIR+OLvcPtSSGruOs2hKd4MT/SBDgPg0v+6TiOOLF68mO7dux/Sz1asIpyxLJ+Ssr2n2hLj48jo\nmhq0VYThqDb/jiLihjFmrud56VXdTyNYoZTUDPpcErnFFUD9JnDCLbBkCqyZWfX9RX7F5zOMvbwf\nDw3rRa/WKaQmJdCrdQoPDesVU8WViATQmpmw8quq7xdCmiQQSv1G2Euk638DLHoPCjdUfV+RSvh8\nhqF9WzO0b2vXUUQkGmT+C7blwOhM10n2UIEVKttyoX5TqBMF8zISGsANX+9dWSgiIuLS9jxIqnT7\nY2d0ijBUJt8Kzw90nSJwjIHyXXYky+93nUZERGLZ9vXQsIXrFPtQgRUKpUWwYrrdezCaLP0I3rka\nFk9ynURERGKVvxyKNmgEKyb99AWU74RuZ7lOEljdzoamXeGrR22PL4kprlcgRzr9+4kESFE+eH5o\nqAIr9iz5COqmQPvjXScJLF8cZNwG67Pt7ygxIz4+npKS4HRcjxUlJSXEx8e7jiES+RIbw/Vf2T5Y\nYUQFVrD5y2HpVNu9PS4KX0x7XQCN2sNXj2gUK4Y0b96cnJwciouLNRJTQ57nUVxcTE5ODs2bR3DL\nFpFwUScB0vqE3RwsrSIMOgPD34yO1YOViYuHjD/AzGftMG1SM9eJJASSk5MByM3NpayszHGayBMf\nH0+LFi32/DuKSC3kLYCcedDnUoiv5zrNHiqwgs3ngzZVNnyNbH0vhyOvtL+rxIzk5GQVCCLi3tKP\n7S4pfYe7TrIPvSMG2/SHIGeu6xTBFVfHFlc7tsHmla7TiIhIlPP7PSbOz2Hwk5mkf9qFwWUPMjE7\nH78/fKYsaAQrmDb9BNMfhHqNoHU/12mCy/Ng/BnQIBVGfOA6jYiIRKmK/Uwzl+dTXFoO1CWf9ox5\nL5sp2Xlhs+WWRrCCqWJlXbS1Z6iMMXDUFbBqBqz+1nUaERGJUpMW5P6iuNqrpKycGcvymbww11Gy\nfanACqalU6FFT2jc3nWS0Oh3NdRPha8edp1ERESi1PjMlfsVVxVKysoZNyM8pqqowAqWsh3w8yzo\ndLLrJKGTUB+Ovxl++tyu6BAREQmwvIKD9+Cr6nioqMAKlq1roF4UNhetSvpI21R12TTXSUREJAql\npSTW6nioaJJ7sDQ7DG5fZtv3x5J6yXDznLBr+CYiItFh5ICOjHkvm5Ky/U8TJsbHMSqjo4NU+9MI\nVjAZY7eTiTUVxVVpkdscIiISdYb0aUVG11QS4/d9f02MjyOjayqDe7dylGxfKrCCwe+HsRkwf4Lr\nJO7Mfw3+2d12dxcREQkQn88w9vJ+PDSsF70a7yKVrfRqUZeHhvUKmxYNoFOEwbHhB1i3EEwMjl5V\naJMOOwtg9gswcIzrNCIiEkV8PsPQvq0Zuv1N+OxeuCkXEhq4jrUPjWAFw5rdfaDaH+c2h0vNukG3\ns2H281Ba7DqNiIhEo8L1kNAw7IorUIEVHGtmQsNW0ChG+l8dyPG/g5LN8N1rrpOIiEg0KtwQtouq\nVGAFmufZEax2/e0k91jWrj+0OQZmPmv/XURERALpghfhuumuU1RKc7ACbddO6HIqdDzJdRL3jIFz\nHrX9wGK92BQRkcAzBuo2dJ2iUiqwAi2+Hgx50nWK8JHWx3UCERGJRp4HH9wK3c6Bw85wnWY/OkUY\naIUbbJsG2WtbLrx2Iayc4TqJiIhEi53bYO7LsPFH10kqpQIr0F4ZDO+McJ0ivCQ2sXsTfvu06yQi\nIhIttq+31w1bus1xACqwAql4s62kW/Z2nSS8xNeDo0fC0qmQv9x1GhERiQaF6+x1klYRRr81M+11\nrG3wXB1Hj4K4eJj1rOskIiISDTSCFUPWfAtxCdDqKNdJwk9Sc+h1IXz3uh3pExERqY2yIqiTGLYj\nWFpFGEhrZkKrI+0pMdnfcTdD446xuQG2iIgEVr8RcNRVrlMckAqsQDrpj2qoeTAtetiLiIhIIIRx\nj0UVWIHU9XTXCcKfvxx+mAiJjaHzKa7TiIhIpPrkbtvI+sTbXSeplAqsQFkzC/Ds9jByEAa+eADq\nJkOngWH96UNERMLY0qnQvLvrFAekSe6B8uU/4IM/uE4R/nw+6H8D5M6Dn2e5TiMiIpFq+3pICs8V\nhKACKzD8flg7B9od6zpJZOhzKdRrBN8+5TqJiIhEotJi2FkADcNzBSGowAqMLStty/60vq6TRIaE\nBpB+Nfz4IWxZ5TqNiIhEmj1NRjWCFd3WZdvrNHVwr7ZjroNm3WFbnuskIiISacpKoElnSGnjOskB\naZJ7IKzLBhNnCwapnuRWcMPXmuQuIiI11+II+N081ykOSiNYgZBxG1z/lRqM1pQxUFoE639wnURE\nRCSgVGAFQkJ9aNnTdYrI9MZl8NYVdqGAiIhIdXz7DEwY5jrFQanAqq3izTDtHti41HWSyNR3OGxa\nDj995jqJiIhEinULYcOPrlMclAqs2sqdD18/sXdFg9RMj/PsKpCZz7pOIiIikWL7urBu0QAqsGqv\nYgVhC50iPCR1EuDokXYES6OAIiJSHYXh3WQUVGDV3rqFkNIO6jdxnSRy9RsBcQmw5EPXSUREJBJE\nwAiW2jTU1rpsaNnLdYrIltQcbs6Cxu1dJxERkXDn99v33TA/c6QCqzZ2ldpJ7j3VYLTWKoqr8l0Q\np/8tRUTkAHw+uGqS6xRV0jtZbdRJgDuWQ3mZ6yTR4ZsnYd6rcOO34ItznUZEROSQaQ5WbRljCy2p\nvZS2kL8ElnzkOomIiISrpZ/Av4+C/OWukxyUCqza+OoRmHKH6xTR4/BzIbkNzBrrOomIiISrrath\n809QN8l1koNSgVUbSz6CDYtdp4gecXXgmFGwagas/951GhERCUfbcsAXDw2au05yUCqwDlX5LlsE\ntNQE94A66iqokwiznnOdREREwlFBDiSn2cnuYUyT3A/VpuWwawekqcAKqPpNYOhTkNbXdRIREQlH\n23LsdJIwF5DyzxgTZ4y5zxiz0hizY/f1/caY6C3gKjq4qwdW4PW6AFK7uE4hIiLhqE06dD3NdYoq\nBaoA+hNwE3AVkA30Bl4BdgL3Beg5wosx9vRg6mGuk0Sn3O8g60U45zH1xRIRkb1O/5vrBNUSqHeu\n44HJnudN3v39KmPMJODYAD1++Ol1gb1IcBSshXmvQJfToMcQ12lERCQceJ69NsZtjmoI1AyxTGCg\nMeZwAGNMD+AUYEqAHj+8eN7e/8gSHN3OgkbtNNldRET2ypkHf0+DFdNdJ6lSoAqsfwCvAj8YY8qA\n74FXPM97JkCPH1625cLDnWDxB66TRC9fHBxzHazO3DvfTUREYlvBz7CrBBKbuE5SpUAVWBcDVwLD\ngaN2f32jMWZkZXc2xlxnjMkyxmRt3LgxQBFCaF02lGyGBqmuk0S3Iy+H+PpqPCoiIta2HHudEv6r\nCAM1B+sR4FHP897Y/X22MaY9MAYY/+s7e573PPA8QHp6euSda1uXDRhocYTrJNEtsTH0vxHiE10n\nERGRcFCQYz94JzZ2naRKgSqw6gPlv7qtnGhtZLpuATTpBHUbuk4S/U6923UCEREJF9vWQnLriJjk\nHqgCazLwZ2PMSuz8qyOBPwD/CdDjh5d12WqEGUr+clj6MXQ9HeLiXacRERFXOp8CrY5ynaJaAlVg\n/Rbb7+oZoDmQB7wAREaziprw+6HHUG2RE0o/fQFvXAoXvAg9h7lOIyIirvQb4TpBtRnPcbuB9PR0\nLysry2kGCXN+PzzVD+qnwqhprtOIiIgLfj/s2GrnXzk8RWiMmet5XnpV94vOOVLBtKMAykpcp4gt\nPh8ccz2snQ05c12nERERFwp+hoc7wvwJrpNUiwqsmvr6CXioHZTvcp0ktvQdDgkNYaZaNoiIxKSC\ntfY6uZXbHNWkTd5qKn8pNGqv/fFCyO/3mLR4O+P9j5I315CWM4ORGZ0Y0qcVPl/4ryQREZEA2NMD\nq63bHNWkKqGm8pdpg+cQ8vs9Rk+YS+byfIpLbVuM/NxtjHkvmynZeYy9vJ+KLBGRWFAxgpXS2m2O\natIpwpoo3wWbV0BqF9dJYsakBbm7i6t926yVlJUzY1k+kxfmOkomIiIhtS0H6jWChAauk1SLCqya\n2Loayks1ghVC4zNX7ldcVSgpK2fcjJUhTiQiIk50OwsG3uk6RbXpFGFN1E2GQQ9Cu+NcJ4kZeQUH\nX7FZ1XEREYkSXU6zlwihEayaSGoGx90ITTu7ThIz0lIOvg9hVcdFRCRKrFsEO7a5TlFtKrBqYv33\nsGW16xQxZeSAjiTGx1V6LNFXxqiMjiFOJCIiIVdaBGNPgDkvuE5SbSqwauLD2+D90a5TxJQhfVqR\n0TV1vyIr0VdOBt8xuIPfUTIREQmZgt0tGpLbuM1RA5qDVRP5S+Hwc12niCk+n2Hs5f2YvDCXcTNW\nkldQQlpKIqPSmzA4LRVfSuT8sYmIyCHaFlktGkAFVvUVb4biTZDa1XWSmOPzGYb2bc3QvpHzhyUi\nIgG0ZwQrct4HdIqwuvKX2Wu1aAgf5WXwwR9gzjjXSUREJJgqurhHyDY5oBGs6stfaq+bqslo2IiL\nhw2LYdk0OGqEti8SEYlWh59rt8ipU9d1kmrTCFZ1dTkNLnzF7kMo4eO4G6FgDfz4geskIiISLC17\nwpGXuU5RIyqwqis5DY44T6Mk4abb2dC4A8x8xnUSEREJlpVfwebI2rlDBVZ1LXwb1v/gOoX8mi8O\njr0Bfp4Fa7NcpxERkUDzPHj9Epj9vOskNaICqzp2lcL718Oid1wnkcoceRkcfS3Ub+o6iYiIBFrJ\nFigriqgVhKBJ7tWzZRV45VpBGK7qNoRzHnWdQkREgqFiBWEE9cACjWBVT8UKQvXACm8/z4GFb7lO\nISIigRSBXdxBBVb17GnRoAIrrM18Gj68HXZud51EREQCJQK7uIMKrOrJXwZJLaFesuskcjDH/xZ2\nFsC8/7hOIiIigXL4uTD8bUhq4TpJjWgOVnWc+SBsX+c6hVSldT9oPwC+fQaOuc42IhURkcjWsKW9\nRBiNYFVHYiNofrjrFFIdJ/zODid//77rJCIiEgiLP4A1s1ynqDEVWFUp3gyf3w8bl7pOItXR5XRo\n2x9KC10nERGRQPj4TpjzgusUNaZThFXZsBi+egTa9YdmatMQ9nw+uGYqGOM6iYiI1JbfD9tyI64H\nFmgEq2p7WjSouIoYxtg/yjUzXScREZHaKNoA/jJIiawWDaACq2qblkOdxIjrvxHzssbDi4Mgb4Hr\nJCIicqj29MDSCFb0yV8KTbvYU08SOXpdCAlJ8M2TrpOIiMihitAeWKACq2rbctXBPRIlNoKjroJF\n78HWNa7TiIjIoehyOtzwDaR2c52kxlRgVWV0Jgx9ynUKORT9b7DXM8e6zSEiIocmoT60OALi67lO\nUmMqsKpiDCQ0cJ1CDkWjttBzGKyYbie9i4hIZPnudfjhf65THBK1aTiYNbNg3itwyt2QnOY6jRyK\ns/4BdZM1h06khvx+j0kLchmfuZK8ghLSUhIZOaAjQ/q0wudTGxQJkW+ehMYdoMdQ10lqTO86B5M7\nD757TVuuRLL6TSCuDpTtsBcRqZLf7zF6wlzufD+b7JwC8gtLyc4pYMx72YyeMBe/33MdUWKB58GW\nVbbAikAqsA5my2q7Eq1+U9dJpDYKN8ITvW3rBhGp0qQFuWQuz6e4tHyf20vKypmxLJ/JC3MdJZOY\nUrgByoojtsDSKcKD2boaGrVXV/BIl9TMNor95kk4ehTUqes6kUhglGy1I+0FOVC0EYrybUPG4260\nxz++C0qLoF6KvTRqB22OhsbtD/qw4zNX7ldc7XnKsnLGzVjJ0L6Rt2xeIsyWVfa6cUenMQ6VCqyD\n2bK6yhciiRAZt8Gr59kJk+lXu04jcmjyl8PWVdDlNPv9i2fCxsV7j8c3gK6n7S2w1mbB5p9gRwGU\nl9rbepwHF71iv/7639D6KGh3/D7zFPMKSg4ao6rjIgFR0WJHI1hRKC5ePbCiRaeTodVR8PXjcOQV\ndl6WSIAEdUL41p9h0Tuw8G3Y8D3UawR/XGkLotP+apevN+4IDVL3X/E88uO9X5eVQP4yMLsLqaJ8\nmHYP4EFKW+h9MfQdDk07k5aSSH5h6QEjpaUk1u53EqmO3hfaDwx1k10nOSTG89xOVkxPT/eysrKc\nZpAYsfgDePMyuPRN6Ham6zQSJSomhP96zlJifBwZXVMZe3m/Qy+yvn0GPh5jv25zDPS6ADoNtB/8\nAjF1Ycc2WPaJHdld8QV4frj4NSbuOJIx72VTUrb/acLE+DgeGtZLpwglZhlj5nqel17V/fQxXmJH\nt7NhxBRof7zrJBJFqjMhvNrFSPFmmP0CHDYIWvWFjhlwyl+g5wXQJAjzUOol26Kt1wWwLQ8Wvgkd\nT2RIQkOmfDOfGblQUr63kKsoGgf3bhX4LCK/9sWDdk7hUVe4TnJIVGAdyNJP4Nsn4fzn1QMrWvh8\n0OEE+7XnafGCBERAJoRvy4Nvn4K5L0NpoZ2e0KovtOxlL6GQnAYDbgXs8vKxjSYweV0B4+KGkedL\nI61pMqMyOjG4t/pgSYjMfdnON1SBFWXWL4KVX0Hdhq6TSKBlPg7LP4WrJqvIklqr9YTwLx6AzH+B\nf5fdeWDA7+3WII75Ln2doYveZegX99vVXA0zoM3joOJKQqG0GArXRewEd1AfrAPbsgrqp0LdJNdJ\nJNDqJcOqGfDT566TSBSoasJ3pcf9fjuKCmDioNeF8Nt5MGxcWBRXgB3x7X0h3DQHzn4U1mXDuoWu\nU0ms2LraXgfj1HiIqMA6kK1hBhXhAAAgAElEQVRq0RC1+l4GyW1g+oN73+REDtHIAR1JjI+r9Fhi\nfByjMn71BvHzHBh3KiyZYr8/6Y9w3jPh+0ZSJwGOuRZ+Nx96/sbe9uOH9rSmSLDs6YHVwWWKWlGB\ndSBbdjcZlehTpy6ceBusnQPLP3OdRiLckD6tyOiaul+Rtd+E8J3bYfItMP402JYL7D7VFimnqes3\nsdc7t8PEG+HZ42DRe24zSfTaUWD7ukVwgaU2DZXxPHjtQts76fibXaeRYNhVCk/2g4YtYNSnrtNI\nhPP7PSYvzGXcjL19sEZldNw7IXzV1zBxtO1pddxNcPKfI3t+Z/4yeO8620W+/01wxn3gq3wUT+SQ\n7TmNHl4fQtSmoTaMgcvfcZ1CgqlOAgx9EpJauE4iUcDnMwzt2/rAqwW3rgZfHbhmKrTrH9pwwZDa\nFUZOg4/vhJlPw5aVcPEEFVkSWGFWWNWUCiyJXZ1Odp1AotnauXZbm57DoM+lcMT5EB9FHdDj6sDZ\nD9tiq3iTiisJrPeug9bpcOx1rpMcMs3BqszCt+zpo+3rXCeRYCveDG+PgKUfV3lXkWrxPJg5Fl4c\nBNP/AeW77CfxaCqufumYa+0pT7AT+PMWuM0jkc/vh+8nQsEa10lqRSNYlclfBptXQP2mrpNIsNVt\nCLnzYfNK6HpGxA9Ji2M7tsGk38IPE+Gws+D8Z2Nn30vPgym32Q2pL5kAnU9xnUgi1fY8KN8Z0RPc\nQSNYldu6GpJb227KEt3i4uHEOyDvO1jykes0Esl2FsILA2HxZDjtXrjkdUhs7DpV6BgDw9+y7Sb+\neymsmO46kUSqPS0awrR1STWpwKqMWjTElt6X2D9k9cWS2qibBH0usTsEDLjVNuqMNQ1bwpX/gyad\n4PVL7G4YIjUVBT2wQAVW5baujvj/sFIDcXVss8d1C+3og0h1+cth2v/B2t2tZk68Y+9+l7GqQSpc\nOck2ap77ius0Eol8cZDaDVLauk5SKzEyOaAGPA+6nAodMlwnkVDqdREU5ED7GH9zlOrbsQ3eHQXL\nPrbNa9tU2RYndiQ1gxFT9vb60ubqUhN9LrGXCKcC69eMgaFPu04hoRZXB066w3UKiRSbV9p5RvlL\n4ZzH4OiRrhOFnwa7FwkV5cMbl8EZ90Pbo91mEgkhnSL8tfJdmocTy3Lm2rkjpUWuk0i42vST3Utw\nex5c8b6Kq6r4d0HRBvjvxXZ1tkhVnjsJvnnSdYpaC1iBZYxJM8a8YozZaIzZYYz5wRhzUqAeP2Tm\nvgQPtLKfuiT27CqFpR/BzGddJ5Fw1bgD9LwArv0cOkXeS1zINWwJl70Dnh9eu8j2nhM5kJ3b7aru\n8jLXSWotIAWWMaYR8DV299JzgO7Ab4ENgXj8kNq62r4QJDZxnURcaH8cdDsbvn4Cija5TiPhJOtF\n2JZnJ+Ce/TA07ew6UeRo2tm2rdi6Gt68wn6QEalMxQrCJpHdogECN4L1RyDP87wrPc+b7XneSs/z\nPvM8b3GAHj90tqyGRu1ic4m1WKfeA6WFMOOfrpNIOPD7YeoY+OD3MPt512kiV/vjYegzULTRbq0j\nUpkoadEAgSuwzgNmGWPeNMZsMMZ8Z4y52ZgIXDayVT2wYl7z7tB3OMx5wRbcErvKdsA7V8PMZ6D/\njXDK3a4TRbbeF8LoTEhO01xXqZwKrP10Am4EVgCDgCeAh4CbAvT4obNlle3fIrHt5DttN+6GLV0n\nEVd2FMCE39htb874O5z5oEa2A6FOgi1c3x0F2e+4TiPhpmEadDsnKnZBCFSbBh+Q5XnemN3fzzfG\ndMUWWE/9+s7GmOuA6wDatWsXoAgB4C+HY66HNlpKHPNSWsNxN7pOIS75y+2p4mHjodcFrtNEF+OD\nbbl238YWR9hRYxGwf2tR8vdmvAAM0xpjVgPTPM8b9YvbrgDGep7X4GA/m56e7mVlZdU6g0hQfD8R\nfvwAfvOCGiXGioIcu9F7fD1bZPniXCeKTtvXwdgMSGwE135htxoS8fvDfqTYGDPX87wqOwsH6rf4\nGuj2q9sOAyJrAkvJVrtyTHMDpEJxPmS/DT/8z3USCYUNP8K40+DDP9jvVVwFT8OWcMF42LQcJt+i\n112xH2gebAOZj7tOEhCBKrD+BfQ3xtxljOlijLkQ+B0QWS3R5/0HHukEO7a6TiLhot/V0KIXfPIX\nKC12nUaCKWcuvHSWbYzZX6eHQ6LjiTDwLlj+KRT87DqNuFawFsqK7KhmFAhIgeV53hzsSsKLgEXA\n34G7gWcC8fghs3U11EuJisl1EiC+ODjrH/bF/5t/u04jwbIqE14ZYvfOG/kxtOzpOlHsGPAHuHGm\nbY8jsW3PCsLI74EFAezk7nneh57n9fE8r57neYd5nvdvLxATvEJpi1o0SCU6nABH/AYy/2XnjUh0\nqVjRltwarvkYmnRynSi2+Hy2bYPfD7NfgJItrhOJK1HUogG02fO+tq6GZr+eSiYCnHEf9BwGSS1c\nJ5FAi68Hl74BKW2gQarrNLFr0zLb0PWnL+CS17SoJBZtWQm+OvbDThQI76n6oeR5sHWNRrCkcilt\noPu59kXf73edRgJh4Vt7u/W36qviyrVm3eC0v8KSD2H+BNdpxIXW/eC4myAuOsZ+VGBV8O+yc216\nnOc6iYSzOeNg/OlQvst1EqmNuS/De9fZ0ZIo2FQ2avS/ETpkwNQ/7z1dJLGj+2A4/W+uUwSMCqwK\ncfHQbwS0VZNROYgGzSEnC7LGu04ih2rms7YtQNfT4bK37d++hAefD857xjYinXiTWjfEEn+5bT4b\nRf/NVWBV2JYLeQs0MiEH130wdD4FPvubXVIskSXzcTs60n0wXPwaxCe6TiS/1qgdDH0KBo7RPKxY\nsmUVPNYdvnvddZKAUYFVYeFb8NyJsKvEdRIJZ8bAuf8Czw8f3hZVn7ZiQv2m0OsiuOBluyeehKce\nQ6HDAPv1rp1us0ho5C+116mHuc0RQCqwKmzLhbrJtg+OyME07gCn/AWWfQLrsl2nkap4Hmzc/eJ9\n1BXwm+ejZhJt1JvxT9tZX0VW9Nu4xF43U4EVfbblQHIr1ykkUhw7Gq77EtJ6u04iB+P3w0d/gucy\n9r6A67RT5Gh+BKxbCNMfdJ1Egi1/KSS1tM2+o4QKrArbclVgSfX54vYWV/nL3GaRyvnL4YNbYPZz\ncPSoqDr1EDO6nQlHXgFfPwG537lOI8G0cUlUjV6BCqy9VGDJofjxQ3jqaPjpc9dJ5JfKd8H719v9\nRTNuhzPu18hVpDrjfmjQDCb9VouQotkJv4u6PUBVYFUY+rT9lCtSE51Phaad7bL/0iLXaaTCdxMg\n+2049R449W4VV5EssRGc9bAd4cid7zqNBEuPodDtLNcpAsq43i4wPT3dy8rKcppBpFZWfQ0vnw3H\n3gBnPeQ6jYCde7Xic+hymuskEgieZ88ypETHFiryK9vybNubtD4RsbrXGDPX87z0qu6nESywG/j+\nOAV2FLhOIpGowwlwzHUw61lY/qnrNLGrtAjeH223vPL5VFxFE2P2FlervlZ7lGizeBKMPw1KNrtO\nElAqsABWfwNvXGo/IYkcitP/Bq2OgqJ810liU8lWePV8WPgm5MxznUaCZenHdrR44Zuuk0ggbVwC\ndVMgqYXrJAGlZjCwt7DSJHc5VPGJMOozO3IioVW40RZXG3+EC1+2czkkOnU5HdocA1PH2BFKbdAd\nHfKX2hWEUTZXUu8GYAushCTbaFTkUFUUVwvehLmvuM0SKwpy4KUzYdNyGP6Giqto5/PB4Cdg53b4\n+E7XaSRQNi6B1G6uUwScRrAAtq21o1dRVj2LA54Hi96BFV9Cm6OhRY+D3t3v95i0IJfxmSvJKygh\nLSWRkQM6MqRPK3w+/f9YpbpJ9m936NPQrr/rNBIKLXrACbfAjEfhqCv3bqkjkalkCxRtiLoeWKAR\nLEs9sCRQjIGhz0C9ZHh3JJQdeG9Lv99j9IS53Pl+Ntk5BeQXlpKdU8CY97IZPWEufr8m8h7QxiVQ\nWmy7Pl85ScVVrMm4DVqnw45trpNIbcXXhxEfQo/zXCcJOBVYYN8QBz3gOoVEi6RmcN6zsOEH+OQv\nB7zbpAW5ZC7Pp7i0fJ/bS8rKmbEsn8kLteiiUitn2P3pPrrDfq+R59iTUB9GfQqHn+06idRWnbp2\nFLJxe9dJAk4FFtihyRZHuE4h0aTr6XDczTBn3AFXtY3PXLlfcVWhpKyccTNWBjNhZPrhfzDhN9Aw\nDU76s+s04pIxtrP7rOegcIPrNHKofvoClnzkOkVQqMAq2QKzX4Atq10nkWhz2r1w+XvQ+qhKD+cV\nHPj0YXWOx5w54+CtqyCtL1wzFRq1dZ1IXNu6Gj6+Cz79q+skcqhmPgOf3+86RVCowMpfDlNut0u8\nRQIprg50OdV+vXbufn3W0lISD/rjVR2PKSVb4IsH4bBBcOX/oH4T14kkHDTtDMfdBN+9Bj/Pdp1G\nDsXGJVG7EbsKrG059lqT3CVYSovg9Qvhjcv2mfQ+ckBHEuPjKv2RxPg4RmV0DFXC8FVeZre9SWwM\nIz+Bi1+z829EKpx4hz1lPOV28Fd+yl3CVFmJ3XmhWfS1aAAVWL9oMqo9riRIEhrAkKcgd57dFHr3\nNh9D+rQio2vqfkVWYnwcGV1TGdw7xov+kq3w2gXwxe7TB00721FBkV+qmwRn3A95C2Ce+s9FlPxl\ngBe1I1h6tdqWA3US7SdkkWA5/GwY+BdbLLToCSf8Dp/PMPbyfkxemMu4GXv7YI3K6Mjg3jHeB2vL\nKnjtIti8Anpf7DqNhLuew2BVZlQ2q4xq+UvtdZSOYKnAquiBpaXeEmwn3g7rF8G0u6F1P+hwAj6f\nYWjf1gztqxHUPX6eDf+9FPy74Ir3oWOG60QS7oyBwY+7TiE1dcT5dhFQSnQuWFGBNfhxO4FWJNiM\ngfOfg/bHQ7vjXKcJTyVbYMIwu8fc8LchtYvrRBJJSrbCV49Av6v1/04k8MVBk06uUwSN5mDVS4HG\nHVynkFgRXw+Ovd7uqbZ1DayY7jpRePD77XViYxg2HkZ+qjdIqbnyUpj7Mky7x3USqY4vHoQlU12n\nCJrYLrD85bZ/ipb3igsf/QlevxhWfe06iVuFG+E/Q2DRu/b7w86ABk3dZpLIlNQcMv4ASz6ElV+5\nTiMHU74LZvwT1nzrOknQxHaBVbQRMv9lV5+IhNqQJ6FRO1tk5cx1ncaNtXPh+ZNg7RwtsZfA6H+j\nndPz8Z36fyqcbVkJ/rKoneAOsV5g7emBpQnG4kCD1L1NM189H1ZH7ye5/XgeZL0EL51p52GM/AR6\nX+Q6lUSD+EQ47a+wLhsWvOE6jRzIxiX2OopXfsZ4gVXRAyvG+w2JO8mt4KrJ0KA5fPH3PT2yot7P\ns+CDW6FDBlz3JaT1cZ1IoknPYXDCLdD2GNdJ5EDyKwqsrm5zBFFsryJUk1EJB43b2xEcz7MrDXft\nhLiE6Gwdsi0PktOgXX8Y/hZ0Od1O+BcJJGPg9L+5TiEHU5Rv33vrJbtOEjSx/cpWuN6+kdXXhFpx\nrH4TO7F7Vym8diFM/XN0zR8p2wFT74R/94UNi+1thw1ScSXBtXUNvD0CCnJcJ5FfO/NBuCW65z/H\n9qvbKXfDH1foRV7Ch6+O7fQ+ayy8dSXs2OY6Ue2tzYIXToGZT8NRV6otioSO54cfP4TpD7pOIpWJ\ni3edIKhiu7IwBuo2dJ1CZC+fD858AM58CJZMgbEDYM1M16kOjefBlDtg3Gm2gejwt+HsR+wkZJFQ\naNwBjh4F370GG350nUYqbPjRjtSv/951kqCK7QJr6hjIfsd1CpH99b8Brt7dgG/iDbZnTKTwPPx+\nj4nf5TJ4/tGkl7/E4DpjmVh0BH5/jEzil/CRcTskJMFn97pOIhVy5sKyT+wUnSgWu5Pc/X6YM84O\nUfa6wHUakf21OxZu+NouxoirY+cxbcuBpp1dJzuwn+fgn/Z/jN75WzLzDMWlDQDIzytkzHvZTMnO\nY+zl/WJ7I2sJrQZN7YrCz++zrVDaa5sq5/IW2KK3SRi/lgVA7I5gFW+y2ypoBaGEs7oN9zbi++ph\nePYE+OKB8JublTMXJlwA409jUk4SmTnlFJfuO0m/pKycGcvymbww11FIiVn9b7AjWamHuU4iYAus\nlr2ifv5zdP92B7Onyah6YEmEOPpau/Luy3/AE33gm6fsqJZrk35nJ7HnZMGp/8f45Jso3lX5CFVJ\nWTnjZqwMcUCJeQkN4NS7tQVTOPCX2yawLXu7ThJ0MVxgqcmoRJjkNLjoFbhuOrTqC5/cZZt1hlpp\nEXz3Xygrsd+3ORpO+QvcshAy/kDe9tKD/nheQUkIQopUYs0seOeayJrTGG1KttrXrxhoAhu7c7BK\nC+05YJ0ilEjT6ki44n27mW2DZva2dYsg60XofbF94Qp0k1LPg5x5MP8/kP0ulG63qwGPOA+OumKf\nu6alJJJfeOAiKy1FqwjFkaINdlPxjidCvxGu08SmBk3h6imuU4RE7BZYvS+yl1jZmkSiT8cT9369\n/nv47nXIGm+Xpve+GDqfYkeXfHGH9vhlJbaIKt4Mzx4P2/Mgvj702F1Utat8svDIAR0Z8142JWX7\nN0pNjI9jVEbHQ8sjUluHnwttj4UvHoReF9pThxJaFTtWxIDYLbAqxMh/aIlyfS6Gw8+GxR/Awjfg\ny4fh26fhzz/b4/P+Y4fmU1pDnUSIrweJTexQPdhmjPnLoGCtnZ+4bpEdCbtgPCQ2hsPPgbS+0GNo\nlVtbDOnTiinZecxYlr9PkZUYH0dG11QG99ZpeXGkYgudFwfZZr4Zt7lOFHteu9DunvKb51wnCbrY\nLbCmjrFvHCf90XUSkcCo2xD6XmovhRtg0097V+n8+CEsnbrv/Vv2gtGZ9usvH4a876BeCiS3sYVX\n19PtMWPgnH9WO4bPZxh7eT8mL8xl3IyV5BWUkJaSyKiMjgzu3UotGsStdv3hsDPh6ycg/Rr7PiCh\n4Xmwdo6dWhADYrfAWvrx3k/vItEmqbm9VBj+pj3VV7gBdpXY1Ye/bPJ3yet2ZCpAOxv4fIahfVsz\ntK/mOEoYOvUeWJVpT3lL6GxdAzu2Qlof10lCIjYLLM+zqwgPP9t1EpHQqd/EXiqTokJIYkiLI+xF\nQitv9+bOMVJgxWabhpIt9lO8VhCKiMSu+a/Z0+MSGnkLwMRB89gobmOzwFKTURERWTvHFlhbVrlO\nEhvSetuu+vH1XCcJidgssHbthNRu0Kid6yQiIuLKSX+0bUymP+Q6SWzoMRQG/d11ipCJzQKrTTrc\nPNs2bBQRkdiU3AqOuRYWvAEbFrtOE91Ki+0imxgSmwWWiIgIwAm/t7t6fH7/Pjf7/R4T5+cw+MlM\n0u+fxuAnM5k4Pwe/X82pD8mqGfBoV/h5tuskIRObqwi/eMA2Urz0dddJRETEpQZN7WmrBql7bvL7\nPUZPmEvm8nyKS22z3PzCUsa8l82U7DzGXt5P/dxqKm8BYKB5d9dJQiY2R7DWZcPW1a5TiIhIOOh3\nld2tYLdJC3L3Ka4qlJSVM2NZPpMX5oY6YeTLWwBNOwes114kiM0Cq3D9vk0YRUQktpUW2dOEq75m\nfObK/YqrCiVl5YybsTLE4aJA3sKY6X9VIUYLrA2Q1MJ1ChERCRfGB/MnwGf3kldQctC7VnVcfqV4\nMxSsUYEVCMaYO40xnjHmqWA8fq143u4RLBVYIiKyW3winHgH/DyLtHplB71rWkpiiEJFibh4OG8s\nHHaW6yQhFfACyxjTH7gWWBjoxw6IXTugwwBtkyAiIvs68gpo1J6RZhKJ8XGV3iUxPo5RGR1DHCzC\nVWxE3+ww10lCKqAFljEmBXgNGAlsCeRjB0x8IlzxPvS+yHUSEREJJ3US4OQxDNn2Bhkty/YrshLj\n48jomsrg3toFpEZWfgUbfnSdIuQC3abheeAdz/M+N8bcE+DHFhERCa7eF+FbO5ux/ToweX0Txs1Y\nSV5BCWkpiYzK6Mjg3q3UoqGmJt9q2zNc8prrJCEVsALLGHMt0AW4IlCPGRTfvw/T/g9GfKCtckRE\nZF++ODj3X/iAoWkwtG9r14kiW+FG2PwTHHm56yQhF5BThMaYbsADwGWe55VW4/7XGWOyjDFZGzdu\nDESE6ivIsT2w6iaH9nlFRCRybFkFU++E8oNPeJcqrPnWXncY4DaHA4Gag3UckAosMsbsMsbsAk4C\nbtz9fd1f3tnzvOc9z0v3PC+9WbNmAYpQTYXrIa4u1EsJ7fOKiEjk2LgEZj4N8191nSSyrf4G6iRC\nWl/XSUIuUAXWRKAX0PcXlyzgjd1fVzmqFTIVPbCMzqGLiMgBdD0D2h4LXz4CZep7dcjWfAtt0u0C\nghgTkDlYnudtBbb+8jZjTBGw2fO8RYF4joApXA8N1QNLREQOwhg45W545VyYMx6Ov9l1osh0xfu2\n0WgMir1O7m2PtZ9MREREDqZjBnQ6GTIfg53bXaeJTPWbQGoX1ymcCHSbhj08zzs5WI9dKwPHuE4g\nIiKR4pR74LsJmux+KL57HYry4YTfuU7iRNAKrLDkefZa869ERKQ62vSzF6m5ua+AVx6zBVZsnSLc\nngf3NYMFb7hOIiIikWTNLLsZtFRPWQnkzIX2x7tO4kxsFViF68FfZvdFEhERqa7Zz8OUO2D7etdJ\nIsPaLPt+204FVmwo3GCvk7SKUEREamDgnbBrp53wLlVb8y1goN2xrpM4E2MF1u5PHknN3eYQEZHI\n0rQz9B0OWS/C1p9dpwl/OwqgdT9IbOw6iTMxWmBpBEtERGropD/Z6y//4TZHJBj0dxg5zXUKp2Jr\nFWHL3nDsaKhTt+r7ioiI/FKjtjDg95CQ5DpJZPDF1hjOr8VWgXXYIHsRERE5FAPvdJ0g/M18Fha9\nC1dNhvhE12mcia3yckcB+P2uU4iISCTz+2Hh25C3wHWS8LTiSyjZEtPFFcRagfXcSfDeta5TiIhI\nJCsrgql/gmn3uE4Sfvx+WPNNTPe/qhBbBVbhBk1wFxGR2qnbEDJuhxXT7UX22vCDPVsUw/2vKsRO\ngbWz0H7qUIsGERGprfRrIKUtfHrv3m3YBFZ/Y681ghVDBVZFi4aGLd3mEBGRyBdfD04eA7nzYPFk\n12nCR6N20Pcyex3jYmcVoZqMiohIIPW5xBZXMT6Zex/dzrQXiaECq2EaDLwLUru5TiIiItHAFwfD\n33CdInwU5YO/HBpqrjPE0inCJh3hpD9CSmvXSUREJJqUFsM3T9nrWDZnPDzWHYo3u04SFmKnwNq+\nzl5EREQCKe87+OQumPmM6yRuLZli9x+s38R1krAQOwXW5/fB8ye7TiEiItGm/fHQ7WzIfBwKN7pO\n40ZBji00Dz/bdZKwETsFVuEGTXAXEZHgOO1eKCuO3Y2gl0yx191UYFWIoQJrvZqMiohIcDQ7DPqN\ngLkvQf5y12lCb8kUaNIZUg9znSRsxM4qwu3roWVv1ylERCRanTwGtq4Gf5nrJKE35El7mtAY10nC\nRmwUWP5yKNqoESwREQmepGZw+buuU7iR0sZeZI/YOEXo+WHIv+Hwc1wnERGRaLctD756JHa20Pn2\naVj0nusUYSc2Cqy4eDjycmh9lOskIiIS7X76HD6/H75/33WS4Csvg+n/gOWfuU4SdqL6FKHf7zFp\nQS7jv1xKXsEO0ho3YGRGZ4b0aYXPp/PEIiISBH0usT2xPv2rXVUXX891ouBZ/TXsLFB7hkpE7QiW\n3+8xesJc7nw/m+x1xeSX+MnO3c6Y97IZPWEufn+MDN2KiEho+eJg0N/thPdvn3SdJrh+nAJ16kGn\nk10nCTtRW2BNWpBL5vJ8ikvL97m9pKycGcvymbww11EyERGJep1Ohu6DYcZjULDWdZrg8DzbnqHT\nQEho4DpN2InaAmt85sr9iqsKJWXljJuxMsSJREQkppzxd+hxHvjiXScJjuJNUK+RTg8eQNTOwcor\nKKnVcRERkVpp3B7Of9Z1iuBpkAo3ZMbOaskaitoRrLSUxFodFxERCYj1P8Ck30H5LtdJAmtXqb1W\nc9FKRW2BNXJARxLj4yo9lhgfx6iMjiFOJCIiMWnTMpj3it1GJ1psWQ3/6AA/THKdJGxFbYE1pE8r\nMrqm7ldkJcbHkdE1lcG9WzlKJiIiMaX7EOh4ku2NVbTJdZrAmPeK3dy6VV/XScJW1BZYPp9h7OX9\neOi87vRq6pFavw69Wqfw0LBejL28n/pgiYhIaBgDZz0MO7fD539znab2dpXCvP/AYYOgUTvXacJW\n1BZYYIusoZ0Nk4suI+vc9Uz+7QCG9m2t4kpEREKr+eFw7PUw9xXImes6Te38ONnu73v0KNdJwlrU\nriLco3CDvdZGzyIi4tLJYyCxCTTv4TpJ7WS9BI3aQ+dTXScJa9FfYG1fZ6+TmrvNISIisa1eMpx0\nh/3a7wdfhJ5EOvsR2J4XuflDJPr/dQrX2+uklm5ziIiIAOQthKePse0bIlHz7tD5FNcpwl4MFFgb\nAGMboomIiLiW3BpKNsOk34K/8h1HwlJpEbw/OnILwxCL/gKr31Uw4gOIi9KtCkREJLI0aGpXFeZk\nwaznXKepvux3YMF/YUeB6yQRIfoLrORW0GGA6xQiIiJ79RwGXQfB5/fB5gjYG9fzIGu8naDfrr/r\nNBEh+gssERGRcGMMnPsYmDhbuIS7nHmQtwDSr9HWONUU/asIRUREwlFKGxg1DVK7uU5StazxkJAE\nvS92nSRiaARLRETElebdbbuDghxY/73rNAfWuAMcO9q2mpBq0QiWiIiIS54H/70Edm6D67+Ceimu\nE+3vpD+6ThBxNIIlIiLikjG2eefWn2HyLbbgChf5y+H798MrU4RQgSUiIuJau/5wyl9sMTP3Jddp\n9pr6Z5j0Oyje7DpJxFGBJSIiEg5OuNXu7/fRn2HdItdpYOnHsHwanPQn27tLakQFloiISDjw+eD8\n5+DIy6FRW7dZdu20ox004pQAAAk8SURBVFdNu8Ix17nNEqE0yV1ERCRcJDWz/bEAdmyDOvWgTkLo\nc8x8FjavgMvedfP8UUAjWCIiIuFm1054+RyYOBr8/tA/f5OO0O9q6Hpa6J87SqjAEhERCTd16kLP\n38Cid2Hqn0K/iq/HUBj8eGifM8roFKGIiEg4OuFWKMqHb5+CBs1C04tqVSasnQP9b7RFnhwyFVgi\nIiLhyBg4/T4o3gRf/B2SW8ORlwXv+dZlw3+HQ4NUOPpaFVi1pAJLREQkXPl8MORJqJsMHU8M3vNs\n+glePR/qJsGVE+211IrmYImIiISzuHg4+2HbusFfblf47doZuMcvWAv/GQqeH66YCI3aBe6xY5gK\nLBERkUixYrrtT/XimbYwCoTc76C0EC5/D5odFpjHFBVYIiIiEaPLqXDxBMhfBs+dCCu+PPTH8pfb\n6+7nwi0LoFXfwGQUIEAFljFmjDFmjjFmmzFmozFmsjGmZyAeW0RERH6h+2C47gu7svDV82DeqzX7\n+dIiyHwc/tUTFk+2t9VLCXzOGBeoEayTgWeA44FTgF3Ap8aYJgF6fBEREamQ2hVGfQZ9h0PjDva2\nLatgydS9I1O/VloEX/8bHu8Nn/4fNO8Oya1ClTjmBGQVoed5g375vTHmCqAAOAGYHIjnEBERkV+o\nmwRDn977/bz/wIx/2knqHTLspPU6dWHwE7ZR6YtnwrqF0PkUOHkMtD3GXfYYEKw2DQ2xo2NbgvT4\nIiIi8ksnj4GWvWDOeFj5le2jldDQHjMGTv4z1G8K7fq7zRkjjBeE9vvGmLeArkC653n7jVUaY64D\nrgNo165dv9WrVwc8g4iIiEigGWPmep6XXtX9Ar6K0BjzGDAAGFZZcQXged7znuele56X3qxZs0BH\nEBEREXEqoKcIjTH/Ai4BBnqetyKQjy0iIiISKQJWYBljnsAWVyd7nvdjoB5XREREgsPv95i0IJfx\nmSvJKyghLSWRkQM6MuT/27vXGDvGOI7j39+ukro2lFTr1gRBqNJKRCwtEXcqXhCppOISKYlIBH2l\nQkQiQSOoW0IQL1yCLaJxi7prE7ZeoChFG+mm1WZjtdTfi5my9Lq7z55nZ+b3SU6yZ2b3nF/+e86Z\n/3nmmZmjx9LWptzxKi1JgyXpfuBSYBqwWtKYclVPRPSkeA4zMzNL56+/gqufWsR733Tz2/piRk93\nz3pmvbCYVxevYO70SW6yBiHVHKyZFEcOvgms6HO7IdHjm5mZWUIvf778P83VRr1/bGDBkm46u5Zn\nSlYPSRqsiNAWbrNTPL6ZmZml9dh7Szdprjbq/WMDjy5Y2uJE9eJrEZqZmTXQijW9g1pvW+cGy8zM\nrIH23WPkoNbb1rnBMjMza6DLTxzPyBHtm103ckQ7V3SMb3GienGDZWZm1kDnHT2WjkNGb9JkjRzR\nTschozl3gi8EPRhDdS1CMzMzG8ba2sTc6ZPo7FrOowv+PQ/WFR3jOXeCz4M1WG6wzMzMGqqtTZw/\ncRznTxyXO0rteBehmZmZWWJusMzMzMwSc4NlZmZmlpgbLDMzM7PE3GCZmZmZJeYGy8zMzCwxN1hm\nZmZmibnBMjMzM0vMDZaZmZlZYoqIvAGklcAPQ/w0o4HuIX6OpnFN03I903NN03I903NN02pVPQ+M\niL239UvZG6xWkLQwIibnzlEnrmlarmd6rmlarmd6rmlaw62e3kVoZmZmlpgbLDMzM7PEmtJgPZw7\nQA25pmm5num5pmm5num5pmkNq3o2Yg6WmZmZWSs1ZQTLzMzMrGXcYJmZmZkl1qgGS9Ijkr6V1Ctp\npaSXJB2eO1dVSdpT0n2Svixr+qOkByXtlTtbVUm6StLbkn6VFJIOyp2paiTNlLRU0u+SFknqyJ2p\nqiSdJOllST+Xr8cZuTNVmaRZkj6VtLbcBnVKOjJ3riqTdI2krrKmayV9KOns3LmgYQ0WsBCYARwO\nnA4IeEPSiJyhKmwsMA64ETgKmA6cBDyTM1TF7QzMB2ZnzlFJki4C5gB3AMcAHwCvSToga7Dq2hX4\nArgO6M2cpQ6mAA8AJwCnAH9SbIP2zBmq4n4CbgKOBSYDbwEvSpqQNRUNn+Re/gM+Bw6LiK9y56kD\nSWcB84BREbE2d56qkjQZ+BQYHxHfZ45TGZI+Broi4so+y5YAz0XErHzJqk9SD3BtRDyeO0tdSNoV\nWANMi4jO3HnqQtIqYFZEPJQzR9NGsP4haRfgMmAZ8H3eNLWyO7AO+C13EGsWSTsCkyhGAPuaTzFi\nYDbc7EaxHV6dO0gdSGqXdDHFyOsHufM0rsEq52f0AD3AmcCpEbEuc6xakDQKuA14JCL+zJ3HGmc0\n0A788r/lvwBjWh/HbJvmAJ8BH+YOUmWSjiq36+uAucAFEbE4c6zqN1iSbi8nX27tNqXPnzxNMTfj\nZOBr4FlJO+fIPlwNoKYbRwQ7gZ8p5mRZaSD1tEH5/7wHbWaZWVaS7gZOBC6MiA2581TcV8BE4Hjg\nQeCJ4XDwwA65AyRwL/DUNn5n2cYfImINxT7vJZI+ohiavRB4csgSVk+/alrOI3i1vHtORPw+VMEq\nql/1tAHrBjaw6WjVPmw6qmWWjaR7gIuBqRHxXe48VRcR64FvyrsLJR0HXA9cni9VDRqsiOim+GAd\nCJW3ndIlqr7+1FTSbsBrFHU8IyJ6hjJbFQ3yNWrbKSLWS1oEnAY822fVacDzeVKZ/ZekORTN1ZSI\n+DJ3nppqYxhs1yvfYG0vSQdTjFS9AawE9gNupthnOy9jtMoqm6v5FBPbpwG7lLsKAVaV3yqsHySN\noRiBObRcdEQ5t21ZRKzKl6wy7gaelPQJ8D5wNcXpROZmTVVR5ej0weXdNuAASRMp3t8ede0nSfcD\nl1J8Xq4u3+8APf5yOjCS7gReAX6kOGjgEorTYWQ/F1ZjTtMgaX+KC0FOAkZR7DJ4F7jN3yIGppw3\n9PYWVk+NiHdal6YeJM0GbtnMqst8ePz2kTSTYh7gvhTncLo+It7Nm6qatvIefyIiZrQ2TfVJ2tIG\n99aImN3KLHUh6XFgKsUX0zVAF3BXRLyeMxc0qMEyMzMza5XKH0VoZmZmNty4wTIzMzNLzA2WmZmZ\nWWJusMzMzMwSc4NlZmZmlpgbLDMzM7PE3GCZmZmZJeYGy8zMzCwxN1hmZmZmif0N9JN5cUtVADwA\nAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "%matplotlib inline\n", "import matplotlib.pyplot as plt\n", "\n", "fig, ax = plt.subplots(1,1,figsize=(10,6))\n", "ax.scatter(data[:,0],data[:,1],label='data',zorder=3,s=60)\n", "x = np.linspace(-3,3,100)\n", "ax.plot(x,poly(x,fit_params),linestyle='dashed',color='C1',label='fit')\n", "plt.xticks(size=14)\n", "plt.yticks(size=14)\n", "plt.legend(loc='upper center',fontsize=14)\n", "plt.show()" ] } ], "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.6" } }, "nbformat": 4, "nbformat_minor": 2 }