Skip to content

Instantly share code, notes, and snippets.

@qpleple
Created May 9, 2018 23:35
Show Gist options
  • Save qpleple/be060c04e65777d4a67020d8377be4d9 to your computer and use it in GitHub Desktop.
Save qpleple/be060c04e65777d4a67020d8377be4d9 to your computer and use it in GitHub Desktop.
Display the source blob
Display the rendered blob
Raw
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"%matplotlib inline\n",
"import numpy as np\n",
"from matplotlib.pyplot import *\n",
"def color(i): return get_cmap(\"tab20\").colors[i % 20]"
]
},
{
"cell_type": "code",
"execution_count": 51,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"E = 0.383330853793\n",
"E = 1.60033993287\n"
]
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAecAAAHVCAYAAADLvzPyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xl83NV97//XmVWj0WjfvW8YbDYvbGGHhCWQkEJoydYm\nNw3pbRa6JM0vXRJuknvTlJaGNm0TAiVLE9JACCELIaEsMZvBGDAYYzAG25IlWfs2mv38/hiNkW3Z\nGkmzz/v5eOhhaeY73+/HsjxvneV7jrHWIiIiIoXDke8CRERE5FAKZxERkQKjcBYRESkwCmcREZEC\no3AWEREpMApnERGRAqNwFhERKTAKZxERkQKjcBYRESkwrnxduLGx0S5dujRflxcREcm5Z599ts9a\n2zTTcXkL56VLl7Jly5Z8XV5ERCTnjDF70jlO3doiIiIFRuEsIiJSYBTOIiIiBUbhLCIiUmAUziIi\nIgVG4SwiIlJgZgxnY0yFMeZpY8wLxpjtxpj/M80xXmPMfxtjdhljNhtjlmajWBERkXKQTss5DFxk\nrT0FOBW4zBhz5mHHfBQYtNauBP4Z+FpmyxQRESkfM4azTRqb/NI9+WEPO+wq4LuTn98NXGyMMRmr\nUkREpIykNeZsjHEaY54HDgC/tdZuPuyQBcA+AGttDBgGGqY5z/XGmC3GmC29vb3zq1xERKREpRXO\n1tq4tfZUYCFwujHmxLlczFp7q7V2o7V2Y1PTjEuLioiIlKVZzda21g4BDwOXHfZUJ7AIwBjjAmqA\n/kwUKCIiUm7Sma3dZIypnfzcB7wDeOWww+4D/mjy8/cCD1lrDx+XFhERkTSksytVG/BdY4yTZJj/\n2Fr7C2PMl4At1tr7gNuB7xtjdgEDwHVZq1hERKTEzRjO1tptwLppHv/ClM9DwLWZLU1ERKQ8aYUw\nERGRAqNwFhGRWUvYBOOxsZkPlDlJZ8xZRETkEL/s+TEJa3FNVHPL1lvoHu+m1d/KDetv4IrlV+S7\nvKKncBYRkVlr9rTx/PDTfOepnzAWHQega7yLG5+4EUABPU/q1hYRkVlr9y0GA3VV1dSPxGDy7tlQ\nPMQtW2/Jc3XFT+EsIiKz1updQDwRZ4m/mR9/YRd//62OgwHdPd6d5+qKn8JZRERmze3wMDA+wlW/\n6WZ3m5el3WE++6NusJZWf2u+yyt6GnMWEZE5WR1v57zbN/Phv1xEZ43h1pve4KO/GWbVV7Rr8Hwp\nnEVEZE7OWfh2Nv3fZ4guHGCk703+7otn8r/Xfox3tF6Y7OLWzsFzpnAWEZE5aapfyb5LT+HkgQPU\nOhu547I7kk9cfz0EAvCP/6iAniONOYuIyJy4HC5aKtqp9lUe+sTXvgaPPAKf+czBSWIyOwpnERGZ\nszbvIvwVFTgdU+Kkrg4efBB+9zt46qn8FVfE1K0tIiJz1l6xCGMMgcrDWs91dfDEE+B2w5YtsGGD\nurhnQS1nERGZs2ZvO4lE4siubUgGczwOn/gE/OVfqot7FhTOIiIyZy6Hi7FQ6MiWc4rTCb/+NTz6\nqAJ6FhTOIiIyL6MTQfwVPqKJyPQHpMagBwZgTDtZpUPhLCIi8zIaDOIwhgPhrqMfVFcH3/lOspv7\nllvUgp6BwllEROZldCKItZaucMfMB1sL3/++urhnoHAWEZF5iScSBMNhukOdMx9cVwe//W3yNqvP\nfCb7xRUphbOIiMzb6ESQA+H9JGxi5oNTAf22t2W/sCKlcBYRkXkbnQgStVEGo33pvaCuDq65Bu67\nTyuJTUPhLCIi8zY2MQGQXtf2VOeeCw8/rIA+jMJZRETmLRyNUun00x2eZTinbrN65BH4+tezUlsx\n0vKdIiKSES3edg6E98/+hamATiRgfBwqK8t+qU+1nEVEJCOave2MxIaZiI/P/sV1ddDQADfcoNus\nUDiLiEiGtHjbAeg51mIkM7nppuRtVmUe0ApnERHJiEZPC4YZVgqbSeo2q8ceg6efzlxxRUbhLCIi\nGeF2uGnwNNMzl3HnqVLbTZ5xRjKgy7AFrXAWEZGMafa20Rfuxs43UF0uiMXgk58sy9usFM4iIpIx\nTZ5WIjbCUHRg/idzueCBB5K3WZVZQCucRUQkY5q9rQD0Rrozc8LUbVZDQxAMZuacRUDhLCIiGVPr\nbsBl3PSGMxTOkAzo22+HcDi5UEkZtKAVziIikjEO46DJ08KByDxmbB+NMfBf/1UWt1kpnEVEJKMa\nvS0MRHrT26FqNqZuN/nZz2b23AVG4SwiIhnV5GklZmMMRfszf/JUQJ9zTubPXUAUziIiklGNnhYA\neiM92blAXR285z3w05+W7CxuhbOIiGRUjbsOl3HTF85SOKdccEHyNqsSHINWOIuISEY5jIMGTxN9\n2Wo5p0wdg77lluxeK8e0ZaSIiGRcg6eZ18ZexlqLyeb2j6mABhgbA7+/JLabVMtZREQyrtHTTNRG\nGI0NZ/9idXXJj099qmTGoBXOIiKScQ2eZgD6Igdyd9Gbby6ZpT4VziIiknF17kYMhv5chnNqqc9N\nm2DLltxdNwsUziIiknFuh5sadx0Dkd7cXji13eRpp8FTTxVtC1rhLCIiWVHvbmIg2pf7C6e2m/z2\nt2F8PPfXzwCFs4iIZEW9p5GR2BDRRCSn1+0KdfC7of8hcdu3oaoqp9fOFIWziIhkRb27EYDBbCzj\neRRdoQ5+1XM3XaEOwolQzq6baQpnERHJinpPEwADkdx0baeCucoZ4MrWP8DnrMzJdbNB4SwiIlkR\ncNXgNC4GczDufHgw+13F2Z2dohXCREQkKxzGQa2rjqHoQNavZbHUuuu5rPnqog9mUMtZRESyqNbT\nkNUx54l4cjZ2e8Uirm77UEkEMyicRUQki+rcDYzGhokmohk/d1eogzs7buPVse0A2V3DO8cUziIi\nkjV17gYAhjPctd0V6uD+nrvxO6tYULEko+cuBApnERHJmhp3PQBDscGMnTM1+ctfIpO/pqNwFhGR\nrKlx1QKZazmPx8a4v4RmZR+NZmuLiEjWuBxuqpzVGZux7XdV8bb6i1nkW1aywQwKZxERybIadx3D\n0fl1a3eFOnDgoKWineMDJ2WossKlbm0REcmqGncdw/MYc06NMT8x8BC2SHeZmi2Fs4iIZFW1q5ZI\nIkwoPjHr105d+euS5veU1O1SxzJjOBtjFhljHjbGvGyM2W6MuWGaYy4wxgwbY56f/PhCdsoVEZFi\nUz05KWwkNjSr15Xakpyzkc6Ycwz4S2vtVmNMAHjWGPNba+3Lhx23yVp7ZeZLFBGRYlbtngzn6BDN\n3ra0X/fK6LayDGZII5yttV1A1+Tno8aYHcAC4PBwFhEROUK1qwZIv+VsrcUYw3mNlxJJhIt6d6m5\nmtWYszFmKbAO2DzN02cZY14wxtxvjFl7lNdfb4zZYozZ0tvbO+tiRUSk+LgdHioclYzGRmY8tivU\nwb1dPyAYH8dpnGUZzDCLcDbGVAE/Af7MWnv4d3grsMRaewrwr8C9053DWnurtXajtXZjU1PTXGsW\nEZEiE3BVMxYbPuYxqTHmSCJcNrOyjyatcDbGuEkG8w+stfcc/ry1dsRaOzb5+a8AtzGmMaOViohI\n0Qq4qo/Zci7nyV/TSWe2tgFuB3ZYa28+yjGtk8dhjDl98rzZ2yNMRESKSpWrhrHYyLQt4p7QfgXz\nYdKZrX028CHgRWPM85OP/TWwGMBa+03gvcD/NsbEgAngOlvufRIiInJQwFVNnDihRBCf03/Ecwt9\nSzin/h0K5knpzNZ+DDjmXd/W2m8A38hUUSIiUlr8zmTojsXGDobzQKSXWncDla4qLm3+vXyWV3C0\nQpiIiGSd3xUAYDw+CiTHmH/a9QM2D/4un2UVLIWziIhknd85Gc6x0UMmf51cvTHPlRUm7UolIiJZ\n53NWYjD0hPfz1OCjmvw1A7WcRUQk6xzGQYWjkt3BVxXMaVDLWUREcsLvqiJgqrmk6T0K5hkonEVE\nJKu6Qh30Rw5Q6axiPD6qYE6DwllERLLm4OQvV4BmTzt9ke58l1QUNOYsIiJZcciSnC3JMeZQfKLs\n181Oh8JZREQybrq1siscPhIkiCTC+S6v4CmcRUQk4wYivVS5Dp2VXeH0ATCRCOaztKKgMWcREcmY\nWCKKy+FmbfU6VledhMvxVsyk9mYOxYPgrs9XiUVBLWcREcmIrlAHd3Z+m+5QJ8AhwQzgdVQAEEqE\ncl5bsVE4i4jIvAV8Pn7Vczce4yXgqpn2mFQ4RxTOM1I4i4jIvAR8Po5ftGTGlb9S4RzWhLAZKZxF\nRGTOBiJ9HL9oCZFodMYlOT0OLwDhuFrOM9GEMBERmbNadz3dgwN0DwzgX33slb8cxoHbuIlYtZxn\nopaziIjMWneok/HYGA7jYF/vAaLxWFqvczs8RBPRLFdX/BTOIiIyK12hDn7Zcxeb+n8769e6jYdo\nIpKFqkqLwllERNI2deWvcxveMevXux0edWunQeEsIiJpmW5JztlyGzcxdWvPSOEsIiIzstayefDR\neQUzgNO4iNn0xqfLmWZri4jIjIwxXNr8eyRsYl77MbscLoJpTh4rZ2o5i4jIUXWFOvif3l8Qt3F8\nzsp5BTOAy7jVck6DwllERKaVGmPuC/dkbJtHh3GSsPGMnKuUKZxFROQIXaEONv/P11l774tc2foH\nB3eUmi8HDhIkMnKuUqZwFhGRQ6RazG5/Daf/84P4b/53sDYj53YYBwmrcJ6JwllERA7hNE7q3Y1c\ncM6f43j8CXj8cRgdzci5HTjUrZ0GhbOIiAAwHksGcLO3jfe0fSA5+WvBAvjZz8Dthq98BSLzW93L\nGENm2uClTeEsIiJ0hTr4787b2TH6ApAM0SM8+yxccUUGWtGK55konEVEylxqjNnvDLDYt2L6g3w+\nuOsuWL4cLrwQwnObvW1QyzkdWoRERKSMzWpJTpcLvvlN2LQJvF4IhaCiInfFlhG1nEVEytREPMj9\nPT+Z3ZKcxsB550FHBxx/PJt+/i9s693Glp4tXHL3Jfxy9y+P+XKLxcE0XeZyCIWziEiZ8jkrOa/h\nkrmtlb1wIc/+1ftZ+4G/5NSXBgDoGu/ixiduPGZAJ2wCYxQ9M1G3tohImekKdRC3cRb6lrCy6oQ5\nn+fzzVtY8IlF3PTve+lsdDMUcDPsdxK+80/hjI9CfX3yo6Hh4OdOhnE4M/iXKVEKZxGRMpIaY651\n17Gg4g+nn5Wdpu7xbrqO9/Ohv1lO40iMmrE4NeNxasajEI3Cjh0wMHDIx2l9PbRcvAbu+VQG/1al\nR+EsIlImpk7+uqz5mnkFM0Crv5Wu8S46Wrx0tHgPPt7mb+Ov3vvlaV/zaO8vOBDqZPm8rlz61PEv\nIlIGZjUrO003rL+BCuehs7UrnBXcsP6Go74mbuM4He55X7vUqeUsIlIGdo3vyGgwA1yx/AoAvvD4\nF4gkIrT527hh/Q0HH59O3MZwGkXPTPQdEhEpYdZajDGcXX8xkUSYCqcvo+e/YvkV3P3q3QDccdkd\nMx4fs1FcCucZqVtbRKREdYU6+Mn+7zIWG8FhHBkP5rmIJCJ4HN6ZDyxzCmcRkRKUGmOO2zimgN7q\no4kIbocn32UUvML5FxMRkYzIxuSvTIkkwriNwnkmCmcRkRJyINxVsMFsrSWSCON1aD3umWhUXkSk\nhFS7alniW8FZ9RcWVDBDcjJYnPgRt1/JkdRyFhEpAf2RA8RtjAqnj7c3v6vgghkgnAgB4HXkf2Ja\noVM4i4gUua5QB/d2/ZAnBx7JdynHFIpPAKhbOw0KZxGRIjZ18te6mjPzXc4xBeNBACqd/jxXUvgU\nziIiRaqQZ2VPZyIxDiS3qpRjUziLiBShuI3zUO8viyaYASbiyXBWy3lmmq0tIlKEnMbJZS2/R4Wj\nsiiCGWA8NobHeLQISRrUchYRKSJdoQ6eG9oMQIOnuWiCGWAsPoLfVZ3vMoqCwllEpEikxphfHXuJ\naCKS73JmbSw2SpUrkO8yioLCWUSkCBw++asYu4bHY6P4nQrndCicRUQKXLHNyp5OJBFhIhGk2lWb\n71KKgiaEiYgUuNHYMAFXNVe0/H5RBjPAaGwIgGp3TZ4rKQ4KZxGRApXaXvG4qrWs8B+P0zjzXdKc\nDUcnw1kt57SoW1tEpAB1hTr4YcetdE7sASjqYAYYjg4CUOOuy3MlxUHhLCJSYFJjzBUOH7XuhnyX\nkxFD0X4qnVV4HN58l1IUZgxnY8wiY8zDxpiXjTHbjTE3THOMMcb8izFmlzFmmzFmfXbKFREpbaUw\n+Ws6Q9EB6tz1+S6jaKTTco4Bf2mtXQOcCXzCGLPmsGMuB1ZNflwP/EdGqxQRKQND0YGSDGZrLYPR\nvpLpBciFGSeEWWu7gK7Jz0eNMTuABcDLUw67CvietdYCTxljao0xbZOvFRGRNNS46ji15nSOrzq5\nZIIZYCQ2RNRGafA057uUojGrMWdjzFJgHbD5sKcWAPumfN0x+ZiIiMygO9TJSHQIYwwbat9WUsEM\n0B/pBaDB05TnSopH2uFsjKkCfgL8mbV2ZC4XM8Zcb4zZYozZ0tvbO5dTiIiUlK5QB7/suYtN/b/N\ndylZ0x85gMFQ527MdylFI61wNsa4SQbzD6y190xzSCewaMrXCycfO4S19lZr7UZr7camJv0GJSLl\nberkrwsaL893OVnTG+mmzt2A2+HOdylFI53Z2ga4Hdhhrb35KIfdB/zh5KztM4FhjTeLiBxdqc7K\nPpy1lt5wN03etnyXUlTSWSHsbOBDwIvGmOcnH/trYDGAtfabwK+AdwK7gCDwkcyXKiJSGqy1bBl6\nvOSDGZJLj4YSEzR5WvJdSlFJZ7b2Y4CZ4RgLfCJTRYmIlDJjDJc0XUXcxqgs4WAG6A4nRzhbKzRH\neDa0QpiISI50hTp44MC9xBJRvM6Kkg9mgJ5wJx7j0WSwWdLGFyIiOTB1jDliI7goj8lRXaFOmr3t\nOIzagrOh75aISJYdPvmr0unPd0k5MREfZzDaR3vF4nyXUnQUziIiWVQus7Kn0xlKrk21wKdwni11\na4uIZJHH4aXR08Lbm95VVsEM0DmxB4/x0KiZ2rOmcBYRyYLR6DBV23bR4HLx7pOvI7lkRPmw1tIx\n8QbtviUab54DfcdERDKsK9TBj/ffwd5XHoa///uyC2aAwWg/Y/FRFvuW57uUoqSWs4hIBk0dY278\ngz+A9/uhsxP27YMzz8x3eTmzd2I3AIt8y/JcSXFSy1lEJEOmnfxlDLz2Grz73fDII/kuMWfeDL5G\no6eZKlcg36UUJYWziEgGhOIT/LrnnulnZV9wAfz4x3DttfDAA3mrMVeCsTF6wvtZWrkq36UULXVr\ni4hkQIXTx4VNl9PkaZt+VvYFF8C990J/f85ry7U3g7sAWKZwnjOFs4jIPHSFOggnQiytXDlzS/Hs\ns5N//uu/QlsbvPe92S8wD3aNv0Ktu15Lds6DwllEZI5SY8zVrloW+5anf8vQuefC5ZdDOAwf+EB2\ni8yx8dgoXeF9bKw9uyxnqWeKwllEZA6mTv56Z8t7Z3cv76mnwoMPwiWXQF0dvPOd2Ss0x3aNvwLA\nCv/xea6kuCmcRURmKSNLcq5dC5s2QWsrjI+Dv/jX27bWsnPsRZq9bdS66/NdTlHTbG0RkVnaE9yV\nmbWyly+Hykq4+mq4+ebMFZgnvZFuBqP9rK46Kd+lFD21nEVE0pSwCRzGwRl157Ou5ky8zorMnPi2\n2+DiiyEUgr/+68ycMw92jG7DZVzq0s4AtZxFRNLQFergrv13MBwdxBiTuWAGWLQIHn0U7roruWBJ\nEXI6HOwa38FK/wl4Hd58l1P0FM4iIjPoCnVwf8/dYMFl3Nm5SFsbbNkCq1YlJ4tZm53rZElTTS0x\nG2VNYF2+SykJCmcRkWNITf7y52I/ZqczeXvV5z8Pn/oUJBLZu1aGtdbV0+xto8mr7SEzQeEsInIU\nveHu+c/Kni2vN9lyfu45+PjHIR7P/jXnqa4qQIXHw8nVp+W7lJKhcBYROYoadz3LK4/LXTAfvHBN\ncg3uxkaIRnN33Tmw1tJe30AoEtFynRmkcBYROUxvuIdoIoLH4eHCpnfmNphTqqrgq1+FgQH4sz8r\n2JDeH9pHoLKSroH+2S3EIsek76SIyBRdoQ7u676TJwYeyncpSfX18PrryR2twuF8V3OErcNPEolF\nOTA8lO9SSorCWURk0tSVvzbWnpPvcpIqKuAnPwGXC667Lt/VHGJ/aB/7Q3vZ39+PLbLZ5YVOi5CI\niJChJTmzxeOBH/0IXngh+XUkknwsj6y1PD34OyqdVfQMDea1llKklrOIlL2ETfBI3/2FGcwpLhds\n2AC//CVcdBEMD+e1nD0Tr9MT3s+G2rPUas4ChbOIlD2HcXBZ89WFG8xTXX45nHwyvOMdyclieRC3\ncTYPPEqNq07raGeJwllEylZXqINnBh/DWkudp6HwgxnA4YB/+7fkntBf/nJeSnh59HmGYgOcWX8B\nTuPMSw2lTmPOIlKWpo4xn1xzGl5TROtBGwP/+I/J26v27EkuXNLampNLB+PjbBl8nAUVS1jiW5GT\na5YjtZxFpOwcPvmrKDdqMCY5KeznP4fzz4eOjpxcdvPAo8RslLPrL8YYk5NrliO1nEWkrBT0rOy5\n+OQnk1tNnncePPQQLF2atUt1Tuzl1fHtrKs5kzpPQ9auIwpnESkzE/Eg1a4a3tlybfEHc8pnPgM+\nX3KxkiyFczQR4Xf9v6baVcu6mjOzcg15i8JZRMpCOB7C66xguf84llauLL2lJj/xieSf//zPcOml\nsGZNRk//9OBjjMSGeVfrdbgdWdo2Uw4qsZ9OEZEjdYU6+GHnrewJvg5QesE8VVMTvP3tby1YkgEd\nE2/y0uizrA2so71iUcbOK0enlrOIlLSpY8yNnjLYa/iDH0zO3r70Unj4YTjhhHmdLhSf4OG++6l1\nN3BG3fkZKlJmonAWkZJVcpO/0nXttdDWBsuWJfeDds7tXmRrLQ/33U8oHuSy5qvVnZ1DJdy3IyLl\nbCQ6VJ7BnHLOOcklP888Ex59dE6neGHkafZOvM5Z9RfS5C2DXocConAWkZIUcNWwofZt5RnMKS4X\nfO1ryZb0b34zq5d2TOzh6cFNLK88jrWBdVkqUI5G4SwiJaU71MFgpB9jDKfWnF6+wZxy0UXw05/C\n9denvVnGSHSIB3vvo9bdwPmNl2uxkTxQOItIyegKdfDLnrt5bOC3+S6lsJx9Nrz8MlRXw0svHfPQ\ncDzErw/cA8Clze/B48jv1pTlSuEsIiVh6uSvixqvzHc5haeyEnp6krtZ/fCH0x4St3F+0/szhqOD\nvKPpKmrcdTkuUlI0W1tEil7ZzsqerdZW+O1vk7dZhcPwkY8cfMpayyN997M/tJcLG9/JAt/iPBYq\najmLSNF7fvhpBXO6TjwxuQb3k0+CtUAymB/r+w27xndweu15HFe1Ns9FilrOIlL03t50JdFEhEoF\nc3pWr4Zbb4Xdu7G/+Q1bz67n+A/8Bb6f/Bvras/Id3WCWs4iUqS6Qh38svsuIokIbodHwTwH1ukk\ndNNXiH3/DkIXncuGP/oaBIP5LktQOItIEUqNMY/FRoglIvkupyhZa3m86jXu/uGHOene7SxsWINZ\nsQLe//7kqmKSVwpnESkqh0/+Uot59uI2zsN9v2L76HOsXHUxvk3PYC64AG67DcbG4FOfOjgeLfmh\ncBaRotGtWdnzFkmE+XXPPbw2/jKn1Z7DmXUXYFpa4Nxz4T//E1auhCeegK9+Nd+lljWFs4gUDZ/T\nT4u3TcE8R2OxEX7WdSedoT2c33Ap62vPOnT1r/e9D7ZtS87ovvVW+O5381dsmVM4i0jBG44OYq2l\nxl2nYJ6jrlAH9+z/PqOxYS5vuYbjAycfeVBtLTzwAOzbl9x68q/+atZrcktm6FYqESloXaEOHnvy\n31lx4jtZ33ROvsspOtZato8+x5MDDxNw1fCu5uuo8zQc/QWBANx/f3KbyfXr4QMfSAb2+vW5K1oU\nziJSuFKTv95+y//QHngTvnc2aBOGtEUSEX7X9wCvB19hsW85FzVegddZMfMLKyuTf/7617B8OVx5\nJTz+eHJ/aMkJdWuLSEGaOiu78bv34d71RnL7Q0lLd6iTn+z/DruDOzmj7jwua746vWCe6hvfgIUL\nk93dl1wCfX3ZKVaOoJaziBScSCLMAwd+euis7F/8IrkedCIBDrUrjiZu42wdepLnhp/C7wzwrtbr\naKtYOLeTeTzw3/8NH/84xGLw7nfDgw++1bKWrFE4i0jB8Ti8XNx0JfXuprcmfzVMjpN+8INwxRXJ\nmcVyiN5wN4/03c9AtI9V/jWc3fB2vA7v/E7qcsHttycXJnnb2+Daa+G++5Jj0pI1+vVTRApGV6iD\nXeOvALDIt2z6Wdmf+xz82Z8ld1cSINnT8MTAQ/y0678IJSa4tPn3uKjpivkH81QOB5xyCjz2GFx/\nvRYpyTK1nEWkIKTGmAOuGpZVrsJpjtIyO+kkuPtuuOaa5D25ra25LbSAWGt5ffwVnhx8hGB8jBOq\nTuGMuvNmP7acDmPgW98CtxvuuAMWLIAvfSnz1xFA4SwyL8H4OM8Pb6a9YjHtFYvwZLKlUkamTv66\nouXaowdzyrnnwjPPJIM5Gk0GRpnpCe3nicGHOBDuotHTwiXNV9Hibc/uRY1JThJbvRr+6Z+Sq4n9\n4R9m95plasZwNsb8J3AlcMBae+I0z18A/Ax4Y/Khe6y1+nVKysJApI+XR1/gxZFnSdgEB0b7GRof\n49z6s3nXqt/H7Si/0Jitw9fKTnuBkSVLYPduxi45n49+bjk7PIO0+lu5Yf0NXLH8iuwWnUcDkT6e\nGdrEm8FdVDr9nN9wGcdVrcVhcjRKaQx8+tNw+ulw/vlsHXiRbQ3biCQiXHL3JSX//c+VdFrO3wG+\nAXzvGMdsstZemZGKRIrIQt8SWuMruXX7f9Bc3cDCmhYu3zzImf/0Ke65/TkqTlpPe8Vi2ioW0uJd\noLCexv7Q3jmvlf1LdtBxSoIbv/I0H/n8Mrro4sYnbgQouYAYjPSxdfgpdo3vwGM8bKw9m5OrN+J2\nePJT0JkQL1HuAAAgAElEQVRn8uZ7384pf/GPnHHDYjadWk3XeOl+/3PN2DQG9Y0xS4FfHKPl/JnZ\nhvPGjRvtli1bZvMSkYJ0yd2X0DXexbL9YW76973E3Q6+87G1jJ5zEisbl+KvqMAYQ8JaxicmGJ0I\nMhIMMjoRJJ5I5Lv8vDFA6t3H6XDM6XuxrXcbkXiYv/vufk7dFeSar6wCwOPwcHLTNMtTFqGqCh9t\n9Q3UBwIkrKVncID9/f3EEvnf1nFb7zb+5tu7effjQ7z7q6vY15Ic1mnzt/Gb92rZz+kYY5611m6c\n6bhMjTmfZYx5AdhPMqi3H6Wo64HrARYvXpyhS4vkV/d4NwDjXkPNeJyfnF/N/e0h2P0MoXGLw+Eg\n4PNRXemn2ldJa30D7Q2NWGsJhsOMTiSDemxignA0mue/TW4EfD5WtC1gZ8c+JiLhOf+SEklEwBj+\n74faeOb6l/FEEkQ8juTjRa4+EKC1roHqykpi8Tid/X10Dw4QK6C9liOJCF/86EJeXFFJsOKteQKp\n/xMyd5kI563AEmvtmDHmncC9wKrpDrTW3grcCsmWcwauLZJ3rf5Wusa7ONDg4bovruA/v/YmGMPP\n3n8Kd1x2xxHHRxNResNddIU76A510OPrImqTYVLp9NPiXUCrt51mbzuNnhZcjtKatzl1jPmPz7th\nXptYpHotEk4Hw1UuqoNx+jwO2vxt037vC91YbJRXxrbxyug2xuNjBFw1nFS9gdVVJ+FZkafu62NI\nff/vvqD+kMdb/eU7gz5T5v2/3lo7MuXzXxlj/t0Y02it1TpvUhZuWH8DNz5xI6F4iP5aN//rc0v5\nmzsP8OcnfHza490ON+2+xbT7kr1HCZtgMNpHd6iT7nAn3aFO3gi+CoADBw2eZlq87TR5W2n2tlHj\nqjt0m78iMufJX0cx9Xs/7HdSMx5nrCHADetvyFDF2RdLRHlz4nV2jr5IZ2gPFssi3zLOCbyDxb7l\nuZvoNQdTv/8pFc6Kovr+F6p5h7MxphXosdZaY8zpJBc26Z93ZSJFIjXx5Zatt9A93o1nwWIiP/wn\nrmg4J7my0kc/eszXO0wygBs8zaxlHQDB2BgHIl30hLs4EN7PK2Mv8tLoVgA8xkujt4VmTytN3lYa\nPS0EXDUFH9h94Z6MBjMc+r0fqdrNilgtH3vbjQU/GSlu43RO7GHX+A7eDO4iaiNUOQOsqzmT46rW\nUuOuy3eJaTn8Z78cZsvnyowTwowxdwIXAI1AD/BFwA1grf2mMeaTwP8GYsAE8BfW2idmurAmhEnJ\n6+mBCy+E666DL3xhXqdK2ARD0X4OhLvpjXRxINzNQKSXBMmxWo/DS6OnhcbJkG/0NFPjrp/5fuEc\niiWiPDHwMBtq35ad/Zivugo+8hF4z3syf+4MiCTCdEy8yRvB19gb3E3EhvEYL8v8x7HSfzztFYsL\nupUsmZGxCWHW2mMuYGut/QbJW61EZKqWFnjoIbjoouS9oX/3d3M+lcM4qPc0Ue9p4nhOAiBuYwxE\n+uiL9NAb7qEv0sP2keeIk5ww5MBJvaeRek8TDe5G6jyN1LubqHT6c9rK7g13U+2qxeus4LzGS7J3\nofp6GBjI3vlnKWET9EcO0Bnaw97gG/SEO0mQoMLhY5l/FUsrV7HItxSnKa05BZIZ+qkQyabW1mRA\nP/xwxk/tNC6avMmu7RMCyceSLewB+iMHDn7sC+7m1cRLB1/ndVRQ526gzt1ArbuBOk/yc78zkPHQ\n7gp1cH/P3SypXMnFTVleCiHP4Ry3cfoiPZNzBzrYH9pHJBEGoMHdxMk1p7HYt5wWb7tayDIjhbNI\ntrW2JndQ+ulPYft2+Nu/zdqlki3sRuo9jaxizcHHJ+JBBiN9DET76I/0MhTtZ3fwVcKJtybyuIyb\nGlctNe56atx11LjqqHbXUu2qnVNre+rkrzPrLsjUX/HochjOqV+Ckr0W3fRGuumPHCBmYwBUu2pY\nXnkcCyqW0F6xiMpsdONLSVM4i+TKWWfB3/xNcj/ieY5Bz5bPWYlvygxxSG6aEEoEGYwOMBTpZyg2\nwHB0kP7IAd4IvorlrfkoLuMi4Kol4KqmylVNwFVDwFVDlTOA3xWg0uk/pDWY6VnZaWlogH37MnrK\nuI0xEh1mKDow+dHPQLSXwUj/weEDl3HT6GnmhMAptHoX0OJdkJu/r5Q0hbNIrqS6uC+6KPn59dfn\ntRxjDD6nH5/TT3vFokOei9s4Y7ERRmJDjESHDv45Gh+hO9x5sLv24LkwVDr9+J0BfI5KuiOdOI2L\n1YET6QnvxxerpMLhw+esxOPwZqdbd5Yt57iNMxEPMhEfJxgfZzw2ylh8lLHYyOTffZjx+Oghr/E7\nq6hzN3Ji9ZLkWL6nmTp3g7qpJeMUziK5lApolwtCIajIwtZ+GeA0zmTXtrsOfEc+H06EGYsNMx4b\nYyw+ejDYgrFRRuPDWGuJ2BCbB3837fk9xoPHUYHXWYHHeHA7kh8e48HlcOM0LlyTHw7jxIEDh3Hg\nwMlbvevJT6y1JEhQ6R2g4cAedg1vJm7jxGyMWCJKxEaIJsJEExFCiRDhRIhwfIKIPXIVMYPB7wxQ\n5QqwoGIxAXcN1a5aat311LrrteuY5IzCWSTXUvsPf+hDcNxx85rFnS9ehxfv5G1bKV2hDt4Ivso7\n667FGEPcxgnFg0zEg4QSEwf/DMcnAzIRIjwZmsH4GJFohGgiQszGiNvYwdvE0lVX0c2Gqrd+IXDg\nxOVw4THeyfB3U+n0U+duwOuooMLpw+f0U+n043NW4nce2T0vki8KZ5F8uemm5H3QUJQBPVVqVrbf\nGWB9zVlUOH04jRO/KzkmPRcJmyBmYyRsPPmBJWGT47xTx8MdODDGgWOhwbz9y/wv48JpXApZKWoK\nZ5F8aW1N3mJ18cVw5ZWwbl2+K5qTwyd/VTin6QefA4dx4DGFt560SC7oV0uRfGptha1bk8G8bVu+\nq5m1vMzKFikDCmeRfPN6IRKB978fvvSlfFczK9FElBpXrYJZJMPUrS1SCDweePDB5G1WkPP7oGcr\nFJ+gwuljceUyFvqWaHxXJMP0P0qkUKRus9qxAyYm8l3NUXWFOriz41Z2j+8EUDCLZIH+V4kUktZW\nuPPO5D3Qt92W72qOkBpjrnRW0eJdkO9yREqWwlmkEEUicPPNBTUGrclfIrmjMWeRQpTB7SYzYSw2\nomAWySGFs0ihSo1Bb9qU70qoclVzRt15LKs8TsEskgPq1hYpZK2tcO21cPfd8OUv5/zyXaEO+sI9\nAJxYvV7BLJIjCmeRYnDOOcmJYjkcg06NMT828CDW2plfICIZo3AWKQapLu4f/Qi+/e2sX27q5K93\nNF2FeWsrKBHJAY05ixSLVEB7PMn7oH2ZWcP6cJqVLZJ/ajmLFJPWVqivh49+NGtd3C+NbFUwi+SZ\nWs4ixejmmzO+1Ke1FmMMFza+k6gN43P6M3JeEZk9tZxFilGqi/vHP4YXXpj36bpCHfy8+0eE4hO4\nHC4Fs0ieKZxFilVqu8lTToHnn5/zaVJjzBPxIHEbz2CBIjJXCmeRYubxJJf6/OAH5zQGrclfIoVJ\nY84ixS613eSFFya/TnMMujvUqWAWKVBqOYuUgtZWePhh2LkzuaNVGvyuAO0VixTMIgVI4SxSKlpb\n4Qc/gPFxuPXWox42GOknYRMEXNVc3nKNglmkACmcRUpNNApf//q0Y9BdoQ7u6fo+zw49kYfCRCRd\nGnMWKTWp26wuuggcDvjbvwUOnfy1JnBqnosUkWNRy1mkFKUCes0aQLOyRYqNwlmkVLW2wtVXE/vR\nD3njlr9SMIsUEXVri5Q41wUXsSx+Iqe0LFcwixQJhbNIqWttpY3WfFchIrOgbm0REZECo3AWEREp\nMApnERGRAqNwFhERKTAKZxERkQKjcBYRESkwCmcREZECo3AWEREpMApnERGRAqNwFhERKTAKZxER\nkQKjcBYRESkwCmcREZECo3AWEREpMApnERGRAqNwFhERKTAKZxERkQKjcBYRESkwCmcREZECo3AW\nEREpMApnERGRAqNwFhERKTAKZ2A4GGVwLJrvMkSkSMTils7+ENbafJciJarsw9lay+adw2zdPUJC\n/9FEJA2v7h/n6deGGRqP5bsUKVEzhrMx5j+NMQeMMS8d5XljjPkXY8wuY8w2Y8z6zJeZPT97fj8/\nfPYNRoIxPv1fL3Dvc535LklECthPt3ayfd8oW/f1c+U3Nuk9Q7IinZbzd4DLjvH85cCqyY/rgf+Y\nf1m5ce9znXz+nhd5ZNcBdh4Y4awlTXzlFzv0n01EpnXvc528vHeMeMJy3/YOOocm+Pw9L+o9QzJu\nxnC21v4OGDjGIVcB37NJTwG1xpi2TBWYTTc9sJOJaJwFwwd48KmX8LocXLSqlZse2Jnv0kSkAN31\ndAdrWmt5aPs+mvftAmAiGtd7hmScKwPnWADsm/J1x+RjXYcfaIy5nmTrmsWLF2fg0vOzf2gCgMt3\nPsZnH/0e++9awwkXXMrLrW/Lc2UiUmgSCcsf79rE0u9/gYu3PM6w18+5f3I7GHPwvUQkU3I6Icxa\ne6u1dqO1dmNTU1MuLz2t9lofALedfjXeRIwvnX4djAzze6cs0SxMETnEru4gzWOD3LXkdC782LdY\nNHIAjAHeei8RyZRMhHMnsGjK1wsnHyt4n710NT638+DXDy86mS+c+2HqF7azp1e/CYtIUjAc55WO\nMX51xYf590Vn0u+vPficz+3ks5euzmN1UooyEc73AX84OWv7TGDYWntEl3Yhes+6BXz16pNYMPlb\nb2t1BfuGx3hjYIwX3hglHE3kuUIRKQTP7R4hHEtw9/N7uf7cZQffMxbU+vjq1SfxnnUL8lyhlBoz\nU/etMeZO4AKgEegBvgi4Aay13zTGGOAbJGd0B4GPWGu3zHThjRs32i1bZjwsd559FjZs4MBoiE/+\n13P8/inLCFQ6uPTU5nxXJiJ59EZPkOffGOWBV/bzB2cu4MLjJ98TJt8zRGbDGPOstXbjTMfNOCHM\nWvu+GZ63wCdmUVtBaw5U8I0PruPffrObk9vr2bSzj3NXN+a7LBHJg8GxCI+83E8snuD3z2h/K5hF\nsqzsVwg7aONbv8g0Byr4+MXLCEZibHtzlKde78tjYSKSD8PBKP/069eo9XlY0lLBRSe0HHrAxhkb\nPyJzpnA+ivZaHxtWVNNeU8m3H9nD5t39+S5JRHJkOBjlT/7rWU5uq8fthktOapn5RSIZpHA+hhMW\nBKitcnHJ6jY+fedzCmiRMjAcjPKB25/ilNZ63E4HF5+kYS3JPYVzyhe/eMRDxhhOW1mD1+3kmlOW\n8JHvPKOAFilhqWD2Gherm2s4eUkAn8c5/cHTvGeIZIrCOeXGG6d9uKrCxfEL/KxqrOaMJQ0KaJES\nlQrmvX1BPnTacmoqXSxrPcbiIkd5zxDJBIVzSnv7UZ9a1e6nqsLJtacuYVGdTwEtUmJSwfxq9xhf\nvepUrIV1y6txTK4ANq1jvGeIzJfCOaXr6OumOB2GdcurCUctX77yZNpqKhTQIiViajD/2/vWE4nA\nitZK6qrcx37hMd4zROZL4ZymxmoPS5oq2N8f4bY/PE0BLVICpgbzNz+4HkfChc/j4IRF/nyXJmVO\n4Zyyfv2Mh5y4JIDH7WBPd4gffuwMBbRIEZsazN/60AYWVFcxEoxxytIAbmcab41pvGeIzJXCOeXZ\nZ2c8xONycPKSAIPjMUbHE9x5/ZkKaJEidHgwb1xSzysdY7TXe2mrr0jvJGm8Z4jMlcI55frr0zps\nQYOX1loPL+8bw+92K6BFiszhwXzB6iae2z2C02E4ZWkg/ROl+Z4hMhcK55Rvfzutw4wxnLKsGmMM\nz78xQlOVVwEtUiQOD+YLj2/mzQMT9I9GOXFJgIqj3dM8nTTfM0TmQuE8B5VeJ2sXVXFgOMLe3hDN\ngQoFtEiBmy6YJ8JxXto7RtPkhE+RQqFwnqNlLT4aAm5e3DNKKBJXQIsUsOmC2VrLc2+MYK3l1OUB\nzLHuaRbJMYVzSmfnrA43Jnnvczxhef6NUay1CmiRAjRdMAPs6wvRMxRhzaIqqipm3D33SLN8zxCZ\nDYVzyhxmXgZ8LtYsqqJrMExnfxhAAS1SQI4WzKFInG17RqmvcrOitXJuJ9dsbckihXPKu989p5et\nbKukzu/ihTdHCEcTgAJapBAcLZittbzw5ijxuGX9iuq5d2fP8T1DJB0K53kyxrB+RQ2xuOX5yfEr\nUECL5NPRghmgsz/M/oEwJyysIuCbQ3e2SA4onDOgutLFCQur2D/wVvc2KKBF8uFYwRyKxHn+zRHq\nqtysbJ9jd7ZIDiicU771rXm9fGV7cqH8598cIRSJH3xcAS2SO8cKZmuTkzfjccuGFTPsOJWOeb5n\niByLwjllnqv9OIxhw4pq4vG3Zm+nKKBFsu9YwQzJ2dldg2HWLMpQd7ZWCJMsUjinZOAex4DPxZrF\nydnbe3tDhzyngBbJnpmCORiOs+3NUeoDbla2Zag7W/dFSxYpnDNsZWsljQE32/aMEgzHD3lOAS2S\neTMFs7WWra+PkLCwYT6zs0VySOGcYanZ21h49vXhQ7q3QQEtkkkzBTPA7p4JekcinLRkjouNiOSB\nwjnlyiszdip/hZOTlgboG4myqzt4xPMKaJH5SyeYRydibN87Skuth6XNvswWkMH3DJHDKZxTfv7z\njJ5uSVMFbXVeXt47xnAwesTzCmiRuUsnmBMJy5ZdwzgdhvXLs9CdneH3DJGpFM4p73pXRk+XWnvb\n7XKwZdcI8YQ94hgFtMjspRPMADs6xhgaj7FuefXstoJMV4bfM0SmUjin/OIXGT+l1+1g/fJqRoIx\nXt43Nu0xCmiR9KUbzH0jEV7dH2RJUwXt9VnaCjIL7xkiKQrnLGut87K8xceuriAHhsLTHqOAFplZ\nusEciSXYsmsYv9fJyUsDOa5SJDMUzjlw4pIAAZ+TLa+/tTnG4RTQIkeXbjBba3l+9wihaIKNK2tw\nOfUWJ8VJP7kp9sgx4UxxOgynraohGktMe3tVigJa5EjpBjPA3t4QnZObWtQH3NktLIvvGSIK55Rb\nb83q6Wsq3Zy4JEDPUITd3RNHPU4BLfKW2QTz6ESMF94cpbHazXG52NQiy+8ZUt4Uzikf/3jWL7G8\nxUdrnYeX9o4yNH7k7VUpCmiR2QVzPGF55rVhHA7YuLImN6uA5eA9Q8qXwjmHjDGsX16Dx+3g6deG\nicanH38GBbSUt9kEM8BLe0YZDsbYsKIGXzZumxLJMYVzjnndDk5bWcN4KM4Lh+1edTgFtJSj2Qbz\n/oEQu3smWNlWSVudN0dVimSXwjnlvvtydqnGag8nLPSzry/EnsN2rzqcAlrKyWyDeTwUZ+vrI9T5\nXaxdVJWjKifl8D1Dyo/COWXDhpxebvUCP03VHl54Y4ThY4w/gwJaysNsgzk5zjwEwGmranA4crzb\nVI7fM6S8KJxTFizI6eWMMWxcWY3HNTn+HDv6+DMooKW0zTaYITnOPDgeY/2Kavz52G0qx+8ZUl4U\nznlU4XFy2qoaxkJxnntj5Jjjz6CAltI0l2Du6H9rnDlry3OK5JHCOc8aqz2sWVRFZ3+Y3T1Hv/85\nRQEtpWQuwTw6EeO53SPUV7lzP84skiMK55SPfSxvlz6uvZLWOg8v7hmlfzQy4/EKaCkFcwnmWDzB\n5leHcJg8jTNPlcf3DCl9CueUPK72Y4xhw4oaKj1Onn51mFAkPuNrFNBSzOYSzNZatu4eYXQizumr\naqn05vl+Zq0QJlmkcE7J88xLj8vBGcfVEI0nePq1YRLT7P98OAW0FKO5BDPA691BOvvDrF1URVON\nJ8tVpkGztSWLFM4pW7fmuwJq/G7WLa+mfzTKS3tH03qNAlqKyVyDuXc4wkt7xmir87IqF+tmp6MA\n3jOkdCmcC8yiRh8rWit5vXuCvb0zTxADBbQUh7kGczAc5+nXhqjyOdmwsjo362aL5JnCOaWtLd8V\nHHTikioaq908t3uEwbFjL1CSooCWQjbXYI7FLU/tHMJaOOO4WtyFtD9zAb1nSOkpoJ/0PNu/P98V\nHOQwhtNX1eJ1O3jq1aG0JoiBAloK01yD2VrLc7uHGQ7G2LiyhoAvDwuNHEsBvWdI6TEzLXyRLRs3\nbrRbtmzJy7Wns+dPP8feT3zukMcWNFSwvLWSWNzy5CuDR7xmcZOPJc0+wtEET786dMTzy1oqWdhY\nQTAc59ldw0c8v7LdT1udl9GJGM/vHjni+QWNFby0Z5RKrwuPy3B4Z96axVU0BDz0j0Z4ee/Ywccj\n8QQ7ukb46bZ9/L9rTmRZQxU7O8aPOP+py6sJ+Fx0DYbZtf/I5zesrKHS66SjL8QbPcEjnj/9uOQv\nEHsOTN8Ff9bxdbicht3dQTr7j1xD/Ny19QC8tn+c7sHwIc85HYa3nVAHwCsdY/QOH3qLmcft4Izj\nagHYvneUgdFDexh8XicbV9YAsO3N0SOWSK3yuVi3vBqA53aPMDYRO+T5Gr+bk5cGANiya5iJ8KG/\nINUH3KxdnHx+86tDRKKHrvDWVOPh+IXJe3Cf2DFI/LAJfq11Xla1+wHYtH2Aw+X7Z2/1Qj/NNV6G\nxqO8+OaR8x+O9rOXctLSALV+N2/0BLn/hR6CkTjHtQSoq3QD6f3s7eub4OV941R6HVQettNUQfzs\n3Xhj8kNkFowxz1prN850nFrOBayqwsmGFTWMTsQYm4iR7q9RHqeDE9qqaQx4+Mh3nuGVriPffI+l\n/Xv/QfPPfjT7gkuI9/XX4Lrr8l1GwfD07GfNn75vVq8ZDkb56q93HBHM6egZCvPyvnEaAm5tASll\nSS3nIvBKxxg7OsZZu7iK4yZbW+k4MBrifbc+RddwiDs+fBpnLG9I74V/+qewdi184hNzrLgEdHXB\nnXfCX/xFvispDD09cNJJcOBAWofPtSsbYHg8yqPbB6n2OTl3bT3OfC40IpJhajmXkNUL/Cxo8LJ9\n7xj7B469xeRUcx6DHhmBQGCO1ZaItjYF81SBAIymd3vffII5FInz5M4h3C7DGatrFcxSthTORSC1\nglhdlZstu4bTnsENcwzo//gPuPbaeVRcAvbsgYsuyncVhcPng717ZzxsPsEci1ue3DlEJGY5a3Wt\nurOlrCmci4TTYTjzuBq8LgdP7RwiGE5vBjfMIaC3boXgkRPAyoq1sHt3vqsoHMbA5s0QOnrPzXyC\n2VrLs68PMzQe47RVNdT60x+fFilFCuciUuFxctbxdcQSlidfGZpxD+ipZhXQn/0s7NqVgYqLmMsF\nsdjMx5WTP/mTo445zyeYAV7aO8b+gTAnLamirc6biWpFiprCuchUV7o447gaRkMxNr+a3hrcKWkH\n9OioxpwrKmDp0nxXUViOMu4832B+vSvIrq4gK1orWdmW/oRHkVKmcC5CzTVe1i+vpnckwtbdI8xm\nxn1aAa1whsZGeOyxfFdRWKYJ5/kG8/6BENv2jNJW5+WkJdqbWSRF4VykFjf5OGGhn319IbbvO3IR\niGOZMaBvugmamjJYbREKheD/+//yXUVhufHGQ3oT5hvMfSMRnnltmLoqNxtX1mjNbJEpFM5FbPUC\nP8uafby2P9ktOBtHDWhr4eqrk9265Sweh3/913xXUVguvTTZo8D8g3kkGOOpnUNUep2ctboWl1PB\nLDJVWuFsjLnMGLPTGLPLGHNEc8IY82FjTK8x5vnJjz/OfKlyOGMMpywL0Fbn5cU9o3T0pX8PNBwl\noINBqK3NUsVFxOWCaPq3rJWFj3wEfvjDeQdzMBzn8VcGcToMZ59Qh9etNoLI4WZcSd4Y4wT+DXgH\n0AE8Y4y5z1r78mGH/re19pNZqFGOwRjDaatqeHzHIFteH8btMrTUpj/bNRXQ77v1KT54+2aWx0b5\ngcPLVX//EJ+9dDXvWbcgi9UXMLdbs7WnuPe5ThKvDvN895PcubORRMJy2x+dNutgDkcTPLZjkHjc\ncu7aOiq9updZZDrp/Mp6OrDLWrvbWhsBfgRcld2yZDacDsNZq2up9rnY/OoQfSORmV80RXOggj96\n2xJicUt4cJhxj4/OoQk+f8+L3PtcZ5aqLnAOB+zcmezmL3P3PtfJ5+95kQO48YeDROMWp8PB8MTs\nehYisQSP7xgkFIlz1vG11MxirW2RcpNOOC8A9k35umPyscNdY4zZZoy52xizKCPVSdrcLgdnn1CH\nz+PkqZ1DDI3P7o3zW4++gQXW9ryeDCRrmYjGuemBndkpuBh0dkIi/XvJS9VND+xkIhonEBpjfecr\nQHLns9n8bKT2ZR6ZiHHGcbU0BDzZKlekJGRqsOfnwFJr7cnAb4HvTneQMeZ6Y8wWY8yW3t7eDF1a\nUrzuZEC7nIbHdwwyEky/W3b/0AQr+vZx44O38teXfSq5ItTk42XriiuOuSJWuUj9DPzD+R9m+eB+\nPrj1l4c8PpN4wrL51SH6R6NsXFkzq2EXkXKVTjh3AlNbwgsnHzvIWttvrU1tinobsGG6E1lrb7XW\nbrTWbmwq91t1sqTS6+ScNXU4jOGxHYNH7FN8NGudE9xx9438w/l/xONLTz34eENVGbdwtEoYw8Ho\nwZnUw74AH772Rj79xI+48PVnaK/1zfj6RMLy9KtDHBiOsH5FNQsbyvwuAJE0pRPOzwCrjDHLjDEe\n4DrgvqkHGGPapnz5bmBH5kqU2aqqcHH2CXVYa3lsxyDjoRkCJhjke/f9X35x0kXcdfI7Dj5sgOGJ\naPq7WZWaMp+xnZqVnUhYPM7kW8W+2lY+/nt/wz/96p/5yuJjz21IJCzP7BqmeyjCqcsCLGmaOcxF\nJGnGcLbWxoBPAg+QDN0fW2u3G2O+ZIx59+RhnzbGbDfGvAB8GvhwtgqW9FRXujjnhOQ63JteHmQ8\ndJSNMuJx+OAHqT/1RNpv+QcW1PowwIJaH//nqjUsrq+c3XaTpeT//T+orMx3FXkx9Xap2/7oNP7h\nvSiO1PMAABQBSURBVCcf/Nk4sHYdu/7PTVz4V3981J2qEgnLll3DB9fLXtZSnt9Hkbkys1n6MZM2\nbtxot2zZkpdrl5Oh8SiPvTyIy2k4d009/orDbl358z+H55+HX/8avEeOBR4YDfG+W5+iazjEHR8+\njTOWN+SocsmXtO9jvvlmuOOO5DKnNTUHH05Yy5bXhumcDGatly3yFmPMs9bajTMdp7v/S1yt3805\na+qIxS2bXh5gbGoX97/8SzKU77ln2mCGOe4HXSrWrIE338x3FTk1qwVG/vzP4YIL4JprIJLs4k4k\nLM9MBvOJCmaROVM4l4FUQMcTlk3bBxmdiMHPfgZ///fwq19BXd0xX1+2AR2NltWY86xX/jIGvv71\nZNf/n/wJ8XiCza8NHezKXqVgFpkzhXOZqPW7OXdNPRbYdtdDJP74j+Hee2HZsrReX5YBXUaztee8\nJKfTCXfeid22jc6/+Du6ByOcsjSgFrPIPCmcy0h1pYtz19Sx8Bc/4tnP3cTAmv+/vTsPjvOu7zj+\n/u4lre7TR+TbsR1fCTYmmAQSJ0BCMy0J4ZjQGRooTE+g/MNMpkzpDDPlaJmStrSlQDIJDBOultY0\nYQLBJgkhcWKn8W3HsmLHsi1bl1fHStrr1z92E8m2jpW0935eMxo92n1296vfPquPnuf5Pb/fllk9\nvuwC+qabpjzcX0rmO1Z2tCLInq9/j+Czu9m6wFi1SJ2/ROZLHcLK0PBojOeOXmIsmmD7ugZa62d3\nLbM6iZWO+QbzWDTB7471EwrH2La6niUtuo5ZZDrqECZTqq70cUtq0oHfHevnbF8GZrMqRV/4Apw4\nke8qsiYTs0s9c7iPgXCM7WsbFMwiGaRwLlOVAS/v2thIfbWfF18N0dGVofmgS8nu3XDxYr6ryIr5\nBnMoHOXpQ32MRRPcvL6RRY2lf/hfJJcUzmUs4PPwzvWNLGoIsP/UIEfODDGb0xwlH9AlOm3kfIO5\nZyDCs4f7weCWjU201JXxEK8iWaJwLnM+r/H2dQ0sb63k+Nlh9p0cIJFQQAMl2Vt7vsF8pmeE5472\nU+H3cOvGJuqqZpwSXkTmQOEseMzYsqqO9UuqOdMzynPH+onE0p8qsWQD+vHHk4NslIj5BLNzjmOd\nQ+xtH6Cpxs+tm5qoqvDO/EARmROFswBgZly3pIa3rq6jdzB5PnFopgkzJijJgH7hBTh3Lt9VZMR8\ngjmecLx8coCjncMsbankpvWNBHz60yGSTfqEyWWWtQZ55/pGIrEEvznYR3do+pmHJiq5gH7wQdi3\nL99VzNt8gnksmuC3R/p5vWeU65ZU89bVdXg9lsVqRQQUzjKJlroAOzY1URnw8Nyxfjq6wml3FCup\ngC6BKSPnE8yXhqPsPthLKBzlxjX1rF9Sg5mCWSQXFM4yqepKH7dubGJBfbIn98sdA8TT7ChWMgFd\n5L215xPMr3eP8PShPgDetaGJtmZdwyySSwpnmZLf5+Ed6xq4rq2a17tHeeZwH+GxKeaFvkJJBPRn\nP5scwrMIzTWYEwnHgVMD7DuZ7Ph12+ZmGmv8Wa5WRK6kcJZpmRnrl9awfW09Q6Nxdh3opat/LK3H\nFn1A33ADLFqU7ypmba7BHB6L88yRfk52jbB6URU3r2+kwq8/ESL5oE+epGVxUyW3bU5ePvP88Usc\nfn0wreuhizqgP/c5ePTRfFcxK3MN5q7+MXYd7GVwJMaNa+q5fkUtHnX8EskbhbOkrabSx62bmlix\nIMir58I8c6Sf4dGZD3MXbUAX2SAkcwnmeMJx4NQgzx+/RDDgZccmnV8WKQQKZ5kVryc5YMnbrq1n\ncCTGroO9nOkZmfFxRRnQRRTOcwnmwZEYTx/q42RXmFULg+zY1ERtUCN+iRQChbPMyZKWSm7f3Exd\n0Mfe9gFePHFpxlHFii6g3/52WLs231XMaLbB7Jzj5Pkwuw70MhKJs31tPTes1PXLIoVE8znLvCSc\n48S5YY52DlPh87BlVd2MMxRpPujMmW0wh8fivNwxQHcowsKGAFtX1VEZ0DCcIrmi+ZwlJzxmrGur\nYcemJvw+4/njl9jbHpp2L7po9qC/852C7hA2m2B2ztHRFebX+3vpH4zylpW1vGNdg4JZpEApnCUj\nGqqT18Sua6ums2eUp/b30tk7OuXIYkUR0J2dcOpUvquY1GyCeWAkxrNH+tl/apCmWj/vvqGZlQur\nNNqXSAFTOEvGeD3GhqU17NjcRDDg4aUTIZ4/donhKSbQKPiALtAOYekGczzhOHJmiF0HehkIx9iy\nqo6brmvQbFIiRUDhLBnXUJ2cUnDz8lp6B6M8tb+Xo2eGiMWv3osu6ICurk4GdAFJJ5idc5zrG+Wp\n/T0cPzvMkuZK3vuWFlYsCGpvWaRIqEOYZNVIJM6h00N09o4SDHjYtKyWtuaKq0JCncRmlk4wD4Rj\nHDw9yMVQhNqglxtW1NFaH8hDtSIyGXUIk4IQDHh525p63rUhOQfwS+0hnj7UR8/A5VNRFuQe9L59\n8NRT+a4CmDmYRyPJXti/PtBL/1CU65fXcvv1zQpmkSKlcJacaKkLcNvmJrauqmMkkuDZI/387lg/\nl4bHp2QsuIB+/nn42c/yWwPTB/NYNMGh04P88pUeXu9Ojol9x5YWVi+uwqND2CJFS+EsOWNmLF8Q\n5L1vaWHD0hr6BqPsPtjHnlcvEUqFdEEFdAFMGTlVMEdiCY6cGeKXr/Rw4nyYa5oqec8NzVy/opaA\nTx9rkWJXWL1dpCz4vMa6tmpWLgzSfj7Mya4w5/rGWNRYwbprqt4M6I9++wU+8chL+TsH7fNBNDrz\nelkyWTCPRuK0nw/z2oURYgnH4sYKNiytoa5KH2WRUqIOYZJ3kViCjq4w7efDROOOpho/115Thc/n\n+MPv7MlfJ7GuLhgYyMsQnlcG89ZljbSfD9PZM0rCwZLmCta2VVNfpbmWRYpJuh3CFM5SMGLxBKcv\njtLeNUx4LEFVwENrg5+/+flBXusdzn1Ah0LJcF66NHevyXgwt18Y4p8+spUqr5+ewSheDyxvDbJ6\ncRU1ldpTFilGCmcpWgnnON83RseFMD0DUczgRPcAz73WzQN3rWX76pbcFPLjH8NPfpL8ypFQOMpf\n/uBlWoJBdqxZiHNQVeFh5YIqli8IUuHX+WSRYpZuOOvfbyk4HjPamitpa65kIBzj1MURvB7j2pY6\nXj0zysDQRd66qoGmGn92B9XIYYew8Fick13DvNDezwc3rwBgUUMFyxcEWdQQ0OAhImVG4SwFra7K\nx/Uratm0rIbj54d4/P+6qPR7eeZwP8GAh8WNFSxuqqC5NpD5KQ+zOHync47BkTjn+8c43zdK/3Dy\ndUYicSqb4faNrdpLFiljCmcpCh6Psb6tluY6P/c/9CKNwUo+9rYVnO4eoePCCF5P8lrqBfUVtNT5\nqa/yzX9vc80auPfezPwCJIO3ZyBCdyjCxVCEkUhy5q66oJe9nT3sOn6BL39w04zTPopI6dM5Zyk6\nE4f6fOj+baxoquXipTEuhCIMj8YB8HuNpho/DTV+GmuSYR0MeHJ2eDgaSxAKx7g0HKV/KEb/UJTh\nsfHaWusDLGwIUFXh5Y+/92La8zGLSHFThzApaVONxR0eS+6d9g5G6RuKMhAePyzt9xq1QR/VlV5q\nKr1UVXgJBrxUBjxU+D34vXZ5eO/ZAw8+CI89dtlrO+eIJxyRmGMkEmc0kiAciTM8mvwaHIm9uVcM\nEAx4aKzx01zrp6Uu8OZe/WymfRSR0qAOYVLSphqopKrCy7LWIMtagwDE4o5QOEpoOMbASIzBcIye\ngQhnehKTPq/fa3i9htdjNB7tZs3RDl77959gY2N03XIn8YQjGndM9j+t32tUV3ppqQtQF/RRV+Wj\nodpHZeDqKRoVzCIyHYWzFK10RhLzeY3m2gDNtZdPABFPOMJjyb3ekUicSCxBJOaIxhLEE454AoLV\nlfhcnAUH9uD8AeK//wd4PYbfZ/h9HgI+o9Kf3PMOBrwEfJbWYXMFs4jMROEsRW2uQ316PclD3LXB\naVaKLIa1K6nxjcHqpSy5tn7e9SqYRSQdulZDil7WJsvYuDE5AEkoBPUKZhHJHYWzlISsBPSFC/Dl\nL2cknBXMIjIbCmcpGRkP6FAIHnlk3uGsYBaR2VI4S0nJaEC/MWXkPMJZwSwic6FwlpKTsYB+Y2zt\nOYazgllE5krhLCUpIwG9eDH89rfJcG5omNVDFcwiMh8KZylZ8w5o56CjIzmn8yz2nBXMIjJfCmcp\nafMK6OFh+MAHkoe3/f60HqJgFpFMUDhLyZtzQL/RISzNvWYFs4hkisJZysKcAnoW4axgFpFMUjhL\n2Zh1QAcC8JnPzBjOCmYRyTSFs5SVWQW0xwN33DFtOCuYRSQbFM5SdmYV0PfcA7W1k96lYBaRbFE4\nS1lKO6Dj8UnDWcEsItmkcJaylVZAm10VzgpmEck2hbOUtRkDesWKy845K5hFJBcUzlL2pg3oO+6A\nBckAVjCLSK4onEWYJqCffBK8XgWziORUWuFsZu8zs+Nm1m5mD0xyf4WZ/Sh1/x4zW5HpQkWy7cqA\n/vqTxxh97TQP/PwY2/7uVxw9N6BgFpGcmDGczcwL/Cvwe8AG4KNmtuGK1T4J9DvnrgW+AXwt04WK\n5MIbAV1T4eWbu0+Cc4z4KojGHV6Ph9BINN8likgZSGfP+Uag3TnX4ZyLAD8E7r5inbuBR1PLPwXe\nbWaWuTJFcmdBbSUeS340DBjxBwCIxBP8w5PH81iZiJSLdMK5DTgz4efO1G2TruOciwEhoPnKJzKz\nPzGzvWa2t7u7e24Vi+TAhYFRAL5x00fZ1zZ+oOjcpZF8lSQiZSSnHcKcc992zm1zzm1rbW3N5UuL\nzMo1DUEAvnXTR+itabzqdhGRbEonnM8CSyf8vCR126TrmJkPqAdmObO9SOH4/J3rCPq9l90W9Hv5\n/J3r8lSRiJSTdML5JWCNma00swBwH7DzinV2Avenlj8E7HLOucyVKZJb92xp4yv3bqatIYgBbQ1B\nvnLvZu7ZcuUZHRGRzPPNtIJzLmZmnwaeBLzAw865w2b2JWCvc24n8BDwfTNrB/pIBrhIUbtnS5vC\nWETyYsZwBnDOPQE8ccVtX5ywPAp8OLOliYiIlCeNECYiIlJgFM4iIiIFRuEsIiJSYBTOIiIiBUbh\nLCIiUmAUziIiIgVG4SwiIlJgFM4iIiIFRuEsIiJSYBTOIiIiBUbhLCIiUmAUziIiIgVG4SwiIlJg\nFM4iIiIFxpxz+Xlhs27gdF5efHItQE++iygQaovLqT3GqS3GqS0up/YYN11bLHfOtc70BHkL50Jj\nZnudc9vyXUchUFtcTu0xTm0xTm1xObXHuEy0hQ5ri4iIFBiFs4iISIFROI/7dr4LKCBqi8upPcap\nLcapLS6n9hg377bQOWcREZECoz1nERGRAqNwFhERKTBlG85m9mEzO2xmCTObssu7mb3PzI6bWbuZ\nPZDLGnPFzJrM7FdmdiL1vXGK9eJm9krqa2eu68ymmd5nM6swsx+l7t9jZityX2XupNEeHzez7gnb\nw6fyUWcumNnDZnbRzA5Ncb+Z2T+n2uqAmW3NdY25kkZb7DCz0ITt4ou5rjFXzGypme02syOpLPmr\nSdaZ+7bhnCvLL2A9sA74DbBtinW8wElgFRAA9gMb8l17Ftri74EHUssPAF+bYr2hfNeapd9/xvcZ\n+AvgW6nl+4Af5bvuPLfHx4Fv5rvWHLXHLcBW4NAU998F/AIwYDuwJ98157EtdgD/m+86c9QWi4Gt\nqeVa4NVJPidz3jbKds/ZOXfUOXd8htVuBNqdcx3OuQjwQ+Du7FeXc3cDj6aWHwXuyWMt+ZDO+zyx\njX4KvNvMLIc15lK5bPdpcc49A/RNs8rdwPdc0gtAg5ktzk11uZVGW5QN59x559zLqeVB4CjQdsVq\nc942yjac09QGnJnwcydXN34pWOicO59a7gIWTrFepZntNbMXzKyUAjyd9/nNdZxzMSAENOekutxL\nd7v/YOpQ3U/NbGluSitI5fJ3Il3vMLP9ZvYLM9uY72JyIXWaawuw54q75rxt+DJRWKEys6eARZPc\n9QXn3P/kup58mq4tJv7gnHNmNtX1dcudc2fNbBWwy8wOOudOZrpWKQo/Bx5zzo2Z2Z+SPKpwe55r\nkvx7meTfiSEzuwv4b2BNnmvKKjOrAf4T+JxzbiBTz1vS4eyce888n+IsMHGPYEnqtqIzXVuY2QUz\nW+ycO5865HJxiuc4m/reYWa/IfmfYimEczrv8xvrdJqZD6gHenNTXs7N2B7OuYm/+3dJ9lsoVyXz\nd2K+JoaTc+4JM/s3M2txzpXkhBhm5icZzD9wzv3XJKvMedvQYe3pvQSsMbOVZhYg2RGopHopp+wE\n7k8t3w9cdVTBzBrNrCK13ALcDBzJWYXZlc77PLGNPgTscqkeHyVoxva44rzZ+0mebytXO4E/SvXM\n3Q6EJpwmKitmtuiNvhhmdiPJjCnJf2JTv+dDwFHn3D9Osdqct42S3nOejpl9APgXoBV43Mxecc7d\naWbXAN91zt3lnIuZ2aeBJ0n2YH3YOXc4j2Vny1eBH5vZJ0lO4/kRgNQlZn/mnPsUyd7t/2FmCZIf\nuK8650oinKd6n83sS8Be59xOkh/C75tZO8kOMfflr+LsSrM9Pmtm7wdiJNvj43krOMvM7DGSvZBb\nzKwT+FvAD+Cc+xbwBMleue1AGPhEfirNvjTa4kPAn5tZDBgB7ivhf2JvBj4GHDSzV1K3/TWwDOa/\nbWj4ThERkQKjw9oiIiIFRuEsIiJSYBTOIiIiBUbhLCIiUmAUziIiIgVG4SwiIlJgFM4iIiIF5v8B\ndmkUYKNNVj0AAAAASUVORK5CYII=\n",
"text/plain": [
"<matplotlib.figure.Figure at 0x10a54ba90>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"figure(figsize=(8, 8))\n",
"u = np.linspace(0, len(x)-1, 100)\n",
"\n",
"# 5 points in V shape\n",
"x = [-1, -0.5, 0, 0.5, 1]\n",
"y = [1, 0.5, 0, 0.5, 1]\n",
"plot(x, y, 'o-', color=color(0))\n",
"\n",
"# Fitting a parabola\n",
"t = range(len(x))\n",
"px = np.poly1d(np.polyfit(t, x, 2))\n",
"py = np.poly1d(np.polyfit(t, y, 2))\n",
"plot(px(u), py(u), color=color(1))\n",
"\n",
"# Error of the parabola fitting\n",
"for t in range(len(x)):\n",
" t_closest = np.argmin(np.sqrt((x[t]-px(u)) ** 2 + (y[t]-py(u)) ** 2))\n",
" plot([x[t], px(u[t_closest])], [y[t], py(u[t_closest])], 'r', lw=1)\n",
"\n",
"# Fitting a line\n",
"t = range(len(x))\n",
"px = np.poly1d(np.polyfit(t, x, 1))\n",
"py = np.poly1d(np.polyfit(t, y, 1))\n",
"plot(px(u), py(u), '--', color=color(1))\n",
"\n",
"# Error of the line fitting\n",
"for t in range(len(x)):\n",
" t_closest = np.argmin(np.sqrt((x[t]-px(u)) ** 2 + (y[t]-py(u)) ** 2))\n",
" plot([x[t], px(u[t_closest])], [y[t], py(u[t_closest])], 'r--', lw=1)\n",
"\n",
"# Rotation of the 5 points by an angle theta and translation of [1.5, 1.5]\n",
"theta = np.pi / 4\n",
"R = np.array([[np.cos(theta), -np.sin(theta)], [np.sin(theta), np.cos(theta)]])\n",
"x, y = np.dot(R, np.array([x, y])) + 1.5\n",
"plot(x, y, 'o-', color=color(4))\n",
"\n",
"# Fitting a parabola\n",
"t = range(len(x))\n",
"px = np.poly1d(np.polyfit(t, x, 2))\n",
"py = np.poly1d(np.polyfit(t, y, 2))\n",
"plot(px(u), py(u), color=color(5))\n",
"\n",
"# Error of the parabola fitting\n",
"print \"E =\", sum([np.min(np.sqrt((x[t]-px(u)) ** 2 + (y[t]-py(u)) ** 2)) for t in range(len(x))])\n",
"for t in range(len(x)):\n",
" t_closest = np.argmin(np.sqrt((x[t]-px(u)) ** 2 + (y[t]-py(u)) ** 2))\n",
" plot([x[t], px(u[t_closest])], [y[t], py(u[t_closest])], 'r', lw=1)\n",
"\n",
"\n",
"# Fitting a line\n",
"t = range(len(x))\n",
"px = np.poly1d(np.polyfit(t, x, 1))\n",
"py = np.poly1d(np.polyfit(t, y, 1))\n",
"plot(px(u), py(u), '--', color=color(5))\n",
"\n",
"# Error of the line fitting\n",
"print \"E =\", sum([np.min(np.sqrt((x[t]-px(u)) ** 2 + (y[t]-py(u)) ** 2)) for t in range(len(x))])\n",
"for t in range(len(x)):\n",
" t_closest = np.argmin(np.sqrt((x[t]-px(u)) ** 2 + (y[t]-py(u)) ** 2))\n",
" plot([x[t], px(u[t_closest])], [y[t], py(u[t_closest])], 'r--', lw=1)\n",
"\n",
"\n",
"axis('equal')\n",
"show()"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 2",
"language": "python",
"name": "python2"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 2
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython2",
"version": "2.7.13"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment