{ "cells": [ { "cell_type": "code", "execution_count": 2, "id": "db13b080", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAA9AAAAH6CAYAAADvBqSRAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjAsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvlHJYcgAAAAlwSFlzAAAPYQAAD2EBqD+naQAARftJREFUeJzt3XvU7XddH/jPb+/9PM+5X3OSkxvJSUiChlQh3IKBGBEwI+WigpdaJUh1BnUN02Fc6IgKFu20rnEqlLZLLQwSsMigLWgGbQUEWiiUEIOAJJDLSULOybnfntve+zd/sGQmQ8x5f+PZ5CS8XmuxliTvfPPZv9/39tk7hq7v+74AAACAhzR4pAsAAACARwMNNAAAAAQ00AAAABDQQAMAAEBAAw0AAAABDTQAAAAENNAAAAAQ0EADAABAQAMNAAAAAQ00AAAABDTQ8Ld429veVl3XPeh/Xvva1z7S5Z0Sx44dq1/+5V+u7/me76lt27ZV13X1tre97UGzv/3bv13XXHNNnXXWWbWwsFC7du2q66+/vu64444H5B7quXVdVzfccMPsPxgA3zSc1w+UntdVVYcPH66f+7mfq0suuaTWrl1bF1xwQf3ET/xE3XXXXbP9QPAoNnqkC4DT3Rve8IbatWvXA/7YE5/4xEeomlNr37599YY3vKEe97jH1bd927fVhz70ob81e9NNN9WuXbvqhS98YW3durVuv/32+u3f/u16//vfXzfffHOdc845VVX17Gc/u37v937v6/763/zN36ybb765nvOc58zq4wDwTcx5/VXpeT2dTuu5z31ufe5zn6tXvepVdemll9Ztt91Wb3nLW+oDH/hAff7zn6+NGzd+gz4hPHpooOEkrrvuunrKU54SZZeWlmp+fr4Gg0fHP9xx9tln11e+8pXauXNnfepTn6qnPvWpf2v2LW95y9f9sRe/+MX1lKc8pd7+9rd/7Vv+iy66qC666KIH5BYXF+tVr3pVfdd3fVft3Lnz1H4IACjn9d9Iz+uPf/zj9clPfrLe/OY310//9E9/LXvZZZfVK17xivqP//E/1kte8pJT/2HgUe7RsWvAaehDH/pQdV1Xv//7v1+/+Iu/WOeee26tW7eujhw5UgcOHKjXvOY1dcUVV9SGDRtq06ZNdd1119XNN9/8oGO8+93vrte//vV17rnn1saNG+sHfuAH6vDhw7W8vFyvfvWr68wzz6wNGzbU9ddfX8vLy19Xyzve8Y668sora+3atbVt27b6oR/6odq9e/dJP8PCwsLfqaG98MILq6rq0KFDD5l73/veV0ePHq1/8A/+wcP+ewHAw+G8fvDz+siRI1VVddZZZz0ge/bZZ1dV1dq1ax/23w8ey/wCDSdx+PDh2rdv3wP+2BlnnPG1//tXf/VXa35+vl7zmtfU8vJyzc/P1+c+97n6oz/6o3rpS19au3btqj179tS/+Tf/pq655pr63Oc+97V/fOpv/Pqv/3qtXbu2Xvva19Ztt91Wb3rTm2pubq4Gg0EdPHiwfuVXfqU+/vGP19ve9rbatWtX/dIv/dLX/to3vvGN9brXva5e9rKX1Stf+cq6//77601velM9+9nPrptuuqm2bNlySp/H/v37azKZ1F133VVveMMbqqpO+o9l33DDDbV27dr6vu/7vlNaCwD8Def1A53svH7KU55S69evr9e97nW1bdu2uuyyy+q2226rn/u5n6unPvWp9d3f/d2ntB54zOiBB/XWt761r6oH/U/f9/0HP/jBvqr6iy66qD9x4sQD/tqlpaV+Mpk84I/dfvvt/cLCQv+GN7zha3/sb8Z44hOf2K+srHztj//wD/9w33Vdf9111z1gjKuuuqq/4IILvvbf77jjjn44HPZvfOMbH5C75ZZb+tFo9HV//KF88pOf7Kuqf+tb3/qQuYWFha89h+3bt/e/9Vu/9ZD5/fv39/Pz8/3LXvayuBYASDmvH1xyXr///e/vzz777Ac8s+c///n90aNH43rgm41foOEk/uW//Jd16aWX/q1//sd//Me/7h9zWlhY+Nr/PZlM6tChQ7Vhw4a67LLL6tOf/vTXjfFjP/ZjNTc397X//vSnP73e9a531Ste8YoH5J7+9KfXb/3Wb9V4PK7RaFTvfe97azqd1ste9rIHfOu+c+fOuuSSS+qDH/xg/cIv/ELzZ34oN954Yy0tLdXnP//5esc73lHHjx9/yPx73vOeWllZ8Y9vAzBTzusHSs7rHTt21JOe9KT6mZ/5mbr88svrM5/5TP2zf/bP6vrrr68/+IM/OKX1wGOFBhpO4mlPe9pD/ktJ/v//xs+qr/6bLf/Fv/gX9Za3vKVuv/32mkwmX/tz27dv/7r84x73uAf8982bN1dV1fnnn/91f3w6ndbhw4dr+/btdeutt1bf93XJJZc8aG3/30P+VLn22mur6qv/spYXvehF9cQnPrE2bNhQP/MzP/Og+RtuuKG2bdtW11133SmvBQD+hvP6gU52Xn/5y1+ua6+9tt7+9rfX93//91dV1Yte9KK68MIL6+Uvf3ndeOONzm54EBpo+Dt6sH/Jxq/92q/V6173unrFK15Rv/qrv1rbtm2rwWBQr371q2s6nX5dfjgcPujYf9sf7/u+qr568HddVzfeeOODZjds2NDyUZpdfPHF9aQnPaluuOGGB22g77rrrvrIRz5SP/mTPzmTywEApJzXDzyv3/a2t9XS0lK94AUveED2hS98YVVVfexjH9NAw4PQQMMMvOc976lrr722fvd3f/cBf/zQoUMP+Bea/F1dfPHF1fd97dq16yH/sbVZWlxcfNB/02hV1bve9a7q+94/vg3Aaemb+bzes2dP9X3/gF/dq6pWV1erqmo8Hn9D64NHC/8zVjADw+Hwa986/40/+IM/qHvuueeU/n2+7/u+r4bDYb3+9a//ur9f3/e1f//+U/L3GY/HdfDgwa/74//1v/7XuuWWW/7Wf2Tune98Zz3ucY+rq6+++pTUAQCn0jfzeX3ppZdW3/f17ne/+wHZd73rXVVV9aQnPemU1ASPNX6Bhhl4wQteUG94wxvq+uuvr2c+85l1yy231A033FAXXXTRKf37XHzxxfVP/sk/qZ//+Z+vO+64o1784hfXxo0b6/bbb68//MM/rJ/8yZ+s17zmNQ85xpvf/OY6dOhQ3XvvvVX11f/N5rvvvruqqn72Z3+2Nm/eXMeOHavzzz+/fvAHf7Auv/zyWr9+fd1yyy311re+tTZv3lyve93rvm7cz372s/WXf/mX9drXvra6rjulnxsAToVv5vP65S9/ef3Gb/xG/dRP/VTddNNNdfnll9enP/3p+p3f+Z26/PLL6yUveckpfQbwWKGBhhn4hV/4hTp+/Hi9853vrH/37/5dPfnJT64//uM/rte+9rWn/O/12te+ti699NL6zd/8zXr9619fVV/9l5k873nP+9r/H9ND+Y3f+I268847v/bf3/ve99Z73/veqqr60R/90dq8eXOtW7euXvnKV9YHP/jBes973lOLi4t1zjnn1A//8A/XL/7iL9aFF174dePecMMNVVX1Iz/yI6fgUwLAqffNfF5v3769PvWpT9Uv/dIv1fve97761//6X9f27dvrFa94Rf3ar/1azc/Pn9oHAI8RXf///+dIAAAAgK/j/wcaAAAAAhpoAAAACGigAQAAIKCBBgAAgIAGGgAAAAIaaAAAAAhooAEAACAwSoMr/Wo86KCbTV8+7adxdlb/89b+Z7P/X6OG9zwcDOPspCZxdtDwHdCxxWN5DZO8hi996Utx9p577omzq8fzeo8dy7NVVbXvaBy9/fbb42y/uhRnL7jggjh7xhlnxNkdO3bE2Z0X7IqzGzZsiLPdmvk4u+PCi+Ps6WDcsD6r2vbtlrOj5ZSZThtqGMzm/BpUvgeeCp/8b/8tzrbsH4973Llx9swzz4yzC+vWxtnT4Y5RNbv7wKBh2JYauoZsy7iD6vIaujzbcg63rPHxeBxnjx8/HmePHDkUZ/fu3Rtnq6oWl/I6Ws7s/ffcF2fvv//+ODtqmBN79uyJsxvXrY+zT3rSk+LspU/99ryGjRvj7CWXPiHOtqyN4XAuzs7P5/eRqqqlyYz2iYaNraVnmNUdY3SSM9sv0AAAABDQQAMAAEBAAw0AAAABDTQAAAAENNAAAAAQ0EADAABAQAMNAAAAAQ00AAAABDTQAAAAEBg90gWcDvq+j7Nd182wkkeX6XQaZ1ue8ep0Nc4uLS3F2f37989k3L1798bZgwcPxtnxieNxtqXeqqpzN2+Os+edd16c3XffPXl23744e+DAgTh75513xtl1X7o9zo5G+Xa5acf2OPvtR4/G2Y3nnBNnW/aqhXVr42zLc6iqGg6HTXlOjfn5+Ti7devWOLthw4Y4Ozc3F2cfjWZ2H5jRnaSl2pZxp+NJnJ1M8mzLHaPF5z73uTh7tGF/3rcvvwvcfffdcbaqanllMc4uLubZyfH87rC6mt/Ntm0/YybjjpdX4mzL3ezz//7fx9kzzsg/2569+T3n/PPPj7Pz82vi7MUXXxxnq6pGw/yMb7kN5Lta1Wqfz4lHil+gAQAAIKCBBgAAgIAGGgAAAAIaaAAAAAhooAEAACCggQYAAICABhoAAAACGmgAAAAIaKABAAAgoIEGAACAwGgWg0776SyGrb7vZzJu13UzGXdW9bYaNJTR8izm+obn1lDDnj174+yBAwfi7Je+9KU4e/DgwTh75NDhOHv06NE4u25+Ls6+5ud/Oc5WVf3eb/7vcXZ5eTnOLi7m2clkEmePHTsSZ1dWVuLsap9/hzgej+PsxrPOiLNf/vKX4+yTrr46zm7fvj3OXn7FFXG2eb9sWPuDlqGn+TnT8k3xN+O3ytu2bYuzCwsLcXY0yq8Ys7k1zFY3PT3O+NSoZXYP8uxkdTXOLp9YzGtoeL5LS0txds/tu+PskSP52XPw4P44+5rX/kKcrar6X/+nn42zw+Ewzq7p82zDEViH79mThxsGHq3m94a//uRn4uyexfw9b9myJc7edvMX4mzLmf2il700zp575s44W1W1ZuPGhnR+aLcc74Mu339m1XOezDfjXQEAAACaaaABAAAgoIEGAACAgAYaAAAAAhpoAAAACGigAQAAIKCBBgAAgIAGGgAAAAIaaAAAAAhooAEAACAweqQL6Pv+kS6hSdd1cbb1s02n0zg7GDx2v/s4duxYnP22v/eUfOC/9zCKeZT4J7/88035f/j6fzyjSh5dXnrtd8bZyWQSZ49OVuLsnj174uy+xcU4e+WVV8bZiy+5JM6uX78+zrZaWcmf2/zoET++Tmst50mLkef+sLTcB1qyLXeSGs3n2QYt95FNZ501kxqa7HqkC6j6uZ9+VVP+jf/yTTOq5NHlZ1/8fXF2fj6f71sWtsTZ8XgcZ7/yla/E2d27d8fZM88/L8623Kmrqi77e0+Ms2vXro2z82vX5NmGdzccDuPsqfTY7cIAAADgFNJAAwAAQEADDQAAAAENNAAAAAQ00AAAABDQQAMAAEBAAw0AAAABDTQAAAAENNAAAAAQ0EADAABAoOv7vk+Cy9OVmRQQ/u2bdV03k2xLva2frZvO5lkMGsZdWcnf84mjx+LsmeecF2db/PRP/aM4u7iynGcXF+Ps8vJsxl08djzOTqfTOFtVNWyIt8zjccMzbrFx48Y42/IsPvDRjz2cck7qmu+4aibjrs7oK88Xv/jFcfY5z3lO09g7d+6Ms2vWrImzGzZsiLOj0SjOzsqght/Qv9+nP/OZOLtt27Y4e+bZO+LscNjwmQf5OXzamE7yaMO+NDej/XluJa+327g5L6LBX/3pn8XZW2+9Nc5+8YtfjLPdeBxnV1dX4+xkkj/f9rtv2xmfGk3yOlrmcMuVpOUZ/+q735UP3OAfPve5Mxl3sRrWXEMvstSQPe+8tvv3c1/wvXH2oosuirOXX3FFnB0MHvnfd092Zj/yFQIAAMCjgAYaAAAAAhpoAAAACGigAQAAIKCBBgAAgIAGGgAAAAIaaAAAAAhooAEAACCggQYAAICABhoAAAACo0e6gBZd181k3L7vZzJua72DQZ5fWVnJx53mn2/Dxi1xtjbm0Q/88fvi7Lvf/e44e+DEgTi7uLIcZ1ue7/LybMZdXcrHbZ3Dw2mebZnHg4bscDiMs0tLS3kNg/x7wWu/45lxtuU5/PnHPhZnW1x9zXfMZNwbb7xxJuNWVT3nOc+Jszt37oyz69evfzjlfNOY1bnWsm75f7XsH32fb9DzCxvyIhby6J2f+XSc3b/73ji7e+/uOLt37944u7q6Gmfnm95Fvo7a3nHb+pzVep5O87nWcrbOatz/+QUvibMtd7M3/9mfxNkWL3v+d8fZyWQSZ6cNc+3WW2+Ns1VVq3+SP4tnPOMZcfbc88+Ps5s3b46zj9SZ5BdoAAAACGigAQAAIKCBBgAAgIAGGgAAAAIaaAAAAAhooAEAACCggQYAAICABhoAAAACGmgAAAAIaKABAAAgMEqDXdfFg/Z9/7CKObWGcbKl3q6b5hUM8hqqqoaVP+O9h/bG2fPOPj/O/un73xdnj37p7ji75/hfx9nfffvb4+wPvfDFcXawMomzc3m0+mn+3ro+/85qMIiXZ43H4zj71ULy6GCQ1zwc5nO+ZdyW/aclOxzO5jvE77z6mTMZ90Mf/VicfdbVV8XZPXvvjbN/+ic3xtmqqmHDcfCc5zwnzq6bX4iz85s350XMSsPaOBXGfX5WTSp/SX3D3tFi2lDDrHQtG2NVVUPNk2l+qKxd2Bhn/+ovPxlnD33prjj76b2fiLNLBw7nNRw6FGe7lZU4uyFOVm06nr+L8bhhbUzzelvvycOG9dwy9mrXsi/lNSxN8xrG43zcEw1rbtiwnK+/7vl5uMFbP/CBOPuy5+bn33R1Nc4uHT8RZ6uqvviFz8fZ1eWlvI6VPPu85z0vzrbcOy+68JI4ezJ+gQYAAICABhoAAAACGmgAAAAIaKABAAAgoIEGAACAgAYaAAAAAhpoAAAACGigAQAAIKCBBgAAgIAGGgAAAAKjR7qA08FgkH+P0PfTODuZTprqOH78RJw97+zz4+yfvv99cXb37t1xdum+PXH2p//pr8TZFr//H/4ozr7k+dfF2a7r4mzL/BkOh3F2PB7H2dFodku55fO1ZGf1jFv0DTW0aPlsfd/H2au/4xlx9iMf+y9x9knffkWcvffee+NsVdWNN97YlE+99KUvjbMbNmyIs03zvWE9f6NNJvn5M53m51qLlrldM1qLs9Ty+TbObYyzd37583H2r+78qzh74iv3xdnDhw/H2cnR/O6ysrISZ4cNZ+Cs5nDLfjDtZziHm5ZSXkdLtsWjbTm3PIeWufZDz39unP39P/uzOHvddz47zq5ZsybOVlUdXV6Ks7fffnuc/cr9e+Pszp074+yzn50/i1PJL9AAAAAQ0EADAABAQAMNAAAAAQ00AAAABDTQAAAAENBAAwAAQEADDQAAAAENNAAAAAQ00AAAABDQQAMAAEBAAw0AAACBru/7Pgmu9KvxoOGQp42Weoc1irNzw2FTHV11cfatb/8/4+zy8RNxdu/evXH2l37ll+Lso82LX/ySOLu6mq+N8XgcZ5eXl+PsZDKJs1Vtc34wyL9na8m26Lp8bbR8tvF0OpNxJzWbPXDSz6beT330v8TZq7/72jhbVXXgwIE42/Ken//858fZn/qpn4qzl16wK84u9fm6W9O1nQd/Vx/++Efi7M6dO+PsBRdcEGeb1syMvs7vGtZiyxlcVbW+1sTZu++/M87e+fkvxNmbbropzh67d0+cXVpairPDlTzbckaMGo61acNevnOl7T2nRuP8LtBSb1VV17DXtFjp5mYy7nLDEbiyshJn907yZ9wy7tFR/j5a7nwnZnRmv/sDfx5nn3/NM+NsVdVKwznc8iwOHDkcZ5/5zLzmV7/61XH28m95Ypw92W/MfoEGAACAgAYaAAAAAhpoAAAACGigAQAAIKCBBgAAgIAGGgAAAAIaaAAAAAhooAEAACCggQYAAICABhoAAAACo0e6gFnp+34m484Nh3G2q65p7A995C/i7Cf2fiLOHj98JM6eOHEizv4PP/U/xNl/9W/+VZydlRe+8EVxtu/zdzeruTZsmGutptNpnB0M8u/ZPtjlz+274uSjz6zmRMu4Le/4e7//xXH2j//jH8XZ1rHvvffeOLt79+44+9nPfjbObty4Mc5u2rQpztZoduv5wbS8/65h3bZoG3dGa6Zh3A21tmnsQ4v74+z+5Ty7b9++OLt37944Ozl6NM9OJnF2uLIUZ0ejhmtnP6vfePK1OKu9fDB4R9tfMP2RhrHz59Y13HVa1vOgf1vDuD88mxoankNDtG1fm830qe99zrPj7B9/OO8tqqq+69qr42zLs2g5Wz/2sY/F2SuvvDLObtmyJc6ee/bjHvLP+wUaAAAAAhpoAAAACGigAQAAIKCBBgAAgIAGGgAAAAIaaAAAAAhooAEAACCggQYAAICABhoAAAACGmgAAAAIjGYxaNd1cbbv+1mUUMPKa2jRNYx71333No39n770gTi7+8474+xwnD/j8XgcZyerq3H2x37wH8TZt/+7G+Ls33/O98TZ0cpcnB0M8u+WJg1z4o+7fNznNiyNQeN87xrq+POWsad59M/zaD23IduypXQN9c5oq2rTUG9N8oInSytx9id+8h81FFH1u//Xb8fZ73/ZS+Ps7V+8Lc5+/CMfi7OjPp/v1157bZyt0XyePQWOHj0RZ8cNZ8SstJytLTbUujh7tPJnVlW1uHgwzt76hb+Osx/+8Ifj7IkTec0b+3zfb7mb9X2+MXXVsC9N8ztGy5m9PJnE2Zb77LR7e5wd/N6PxdmqqsGPvzPOTupH84H/z2EcbdklxnV9nv2hd8XZQfeDcXbYUPFwkM+JFi3raNIwL1tcd+3VTfkbP/jROHvNs6+Ks13DPrGwsBBnP/GJT8TZyy67LM6ee/bjHvLP+wUaAAAAAhpoAAAACGigAQAAIKCBBgAAgIAGGgAAAAIaaAAAAAhooAEAACCggQYAAICABhoAAAACGmgAAAAIjB7pAk4H88OFODutPs7evP+zTXXs2bMnzo7H4zg7XZ3OZNzJZDKT7Eu/98Vxdnm5i7N9n7+7lnrfN8jHbfFnDdnnNI7ddflzaxn7PzVkn9uQfbRpmWstptN8Lbdo2Xvm5+ebxn7FP3plnP237/6dOPvkpz4lzn70ox+NswsL+XnwjGc8I85uXL8+zp4KBw8ejLPLy8szqaFpHTTsSRsqf0fHajHOLi7l2aqqAwcOxNlbb701zrasx40bN8bZWe0fLVpqmFW9LcMObxg2jPyjefTH39EwbtvYLed7d32eHQ5/N882POOGcps+22CQ/0bYdXnBTeNOGz5cg5a10XofuebZV8XZD//Ff4mzV139tDg7GuXt6X333Rdnb7755jh77TXf/ZB/3i/QAAAAENBAAwAAQEADDQAAAAENNAAAAAQ00AAAABDQQAMAAEBAAw0AAAABDTQAAAAENNAAAAAQ0EADAABAYJQGh90wHnQ8HcfZvu/jbIvlvNwm9x/eF2c/+d8+0TT2F2/76zjbdV2cnV+ZxNnp6mpew2r+nte0vOd+Gkfn5uIpPLO59gPD/Huo9zTU8P3r18bZzY2fbdzwjFue29+f5OMuNIw7mTTM4WnDsxjm66jlEQ/zcms6zZ/ZYJAXMR7n63PY8Bzuv/++OFtVNZ3me0qLrsuf2/79e+Psrbd+Ic5++tOfjLPf87zr4uyp0PKZjx8/Emdb9oOWc6qrPNti2LBujx441DT27tvviLP77tsTZ0d9/iyG+TKoaribtWipt8Z5wd2M7g3jhhrmXj6XZwcNF8/h9Xm2qgaVjz0Y5HeSQcPvaHNzr4qzx/9VwyF4/WIcHa42POOGEkaVr41Rw7jdasMCneYDt+zDs7r7tmo5D4bD/D3fd19+J7n55pvj7Mn4BRoAAAACGmgAAAAIaKABAAAgoIEGAACAgAYaAAAAAhpoAAAACGigAQAAIKCBBgAAgIAGGgAAAAIaaAAAAAiM0mBXXTxo180mO51O4+zGWhNnj9VynP1PH//TOHvbbbfF2aq2z9diMpnMJDvo+zjbN2SHw+FMxm0xNzcXZ4ejfA7/w4Ya1ubDNq2jqqppQ7zlGXfT2cyJlux4PM6zfb7mWtbG0kq+p6yursbZaUN2MMi/H215vq371P79++PsC170vXH2/f/1j+PsFd/+xDh78ODBOHvLLbfE2e953nVx9lS444474uzFF18cZ2d1vq+v+Ti7VPla7Gslzh47dizOVlUdOHAgzrbsHwsLCzMZd1bnZcuwTfv+w6gl0XTPadhHB4O2c7jFqKGOlnU3HOR3nZZ70Xzl97jxXH5mzzX87tdyVg26fNyWO2rL/JmVWa37qqprnn1VnP3wR/5LnL3q6qfF2ZZ5ed9998XZk3nk3ywAAAA8CmigAQAAIKCBBgAAgIAGGgAAAAIaaAAAAAhooAEAACCggQYAAICABhoAAAACGmgAAAAIaKABAAAgMHqkC2jR9/1Mxp1MJ3H281/8fJzdfe/upjr6pXFTPjW/nGeH02mc7cZ5tsWalmnZ5dHBIP++aDTIa1ieH86khulcPm7XNTyIqpoO83zLuuumebZl3MkkX6NN49Zs6j1y7GicPXbsWJzthvmcWFlZibMtz7fVtGFPOXjw4ExqWL9+fZw9ejR/d5/97GcfTjnfEHfdcWecPXzwUJydNmw1wxl9R9+y2y0ePh5nv3LH3U11fOE/35TXceRInD1ndS7Orq6uxtkNK7O5Y6xp2PdbzA8azp7uvXF2MnlhnF1Ym++N69bMx9nWM3u1y+sYDPKxJ4N83K7h/jJoWPtzC/l9a67hetgP8/BoNT+nWs7Llvfccj+shr6lVVMdM9Ly3IYN96KW8/1kHvmnBAAAAI8CGmgAAAAIaKABAAAgoIEGAACAgAYaAAAAAhpoAAAACGigAQAAIKCBBgAAgIAGGgAAAAIaaAAAAAiMHukCBl3ewy+M5uPs/pUjcXbSTeLsiRMn4uzKykqcraoaTPqZjL2+z5/bdDqNs90kz7YYNsyJFl3Xxdmm59Dly2gwaPhsDdnRqG0pT2b01dkgn8JNz3g4HMbZvs+LGPezmcMbNmyYybiLy8txtmW+Ly0txdmW99aab6njqc94Spz9+Mc/EWef+ayr4uzBgwfj7DfacsNcOXz48ExqWNNwxTheeb2Dhu/+W/aD1ufQ8v6PHTsWZ+cbPl/LOm95Fk376Di/Q7XUO57+UT7un7wozk6v+w/5uN3LG7L5Z5ullvOyZTdv2/tnM4dn9Yxb5vujzekyL5997XfE2b/4yMfi7FVXPy3Ott5fHopfoAEAACCggQYAAICABhoAAAACGmgAAAAIaKABAAAgoIEGAACAgAYaAAAAAhpoAAAACGigAQAAIKCBBgAAgMAoDfbVx4MOu2GcnfSTOFtdHl1dXomz4/E4zk4meb1d11BwVU27fOxJNWT7aZ6d5uPOTfM50WKwmtfb4s1dPu5rRg2fbZh/DzUY5DVMpw3PYdL2Lvph29xMrTR8Jdf3ec0tz6JpPTesjZZ6VxpqmFtYiLPLDetz2LBn12q+X7acBVVVLdN43PA+BrPZfprW83Q2y+iUOHLkSJw9ePDgDCvJdA0H/OJkMc6OF5fzIlYb7iNVNV3N13nXcF4ur+Y1z8/Px9mWNTNtqHfQtSzGhnH/5AVxtvvv3h9nV+r5cfb48qY4O54sxdnRKL5+V1VVP8jXx9JSXse44d213GkPVT4vG66ztdpwB5/0+frsJ6txdrySP9/pOD9bu4Z+aNCwX04b7i5VVV3DXGu5F83KYJCf2aeyXr9AAwAAQEADDQAAAAENNAAAAAQ00AAAABDQQAMAAEBAAw0AAAABDTQAAAAENNAAAAAQ0EADAABAQAMNAAAAgVEa7KuPB+2qi7ODbjY9/GQyibPj8TjOHjt2bCY1VLV9mzEcDuNs1+fvo+vybDXMib7Ps4NB/tneFCerqmGu/UZDvT/dMH9anu9y5eO2zOGqqtUZfXU2aRi3ZU60rKXpdJqP2zdkW/aUhhpa6l1ZWclraJwTqbY9omowyCfFrLItWubl6ezw4cNx9o477phdIaFh5fv+pMvX4tGjR+Psvn374mxV1eLiYpxtWeeta2wWNbRkG658M/T342TL8z1x4kScXR3M6q5VNW2IN+2jC3NNdaSWW3qGlfzD9Q2frWUOt5yXLXeBWfUi1eXvbVb7SevYsypjNIpb2ea+7KH4BRoAAAACGmgAAAAIaKABAAAgoIEGAACAgAYaAAAAAhpoAAAACGigAQAAIKCBBgAAgIAGGgAAAAIaaAAAAAiM0uB4Mo4HnRvOxdlhdXG2xWDSzyTbr+TPYTiNo1VVNVd5HQtz8aur+ekkzk4rL3phmtfbYs04H/d/aRj3nze8559vmJYLK/nzbTFtmD9dl8/LqqrJMP/ubDLJP9+4y5/xtOEDvrNh3O9vqHelz7N93/DZxqszGbca1vKg4fmO+jzbN+xTX9WwmLp8Xg4ast+MVhvO7M994fNx9gMf+EBexPPzaMvb7BrOnrk1C3F2peGZVVUtLy/H2YWFvI6uYQ/rJw3rvGEtdsP8jjHX8Ny6rmE/uK7lEMyjc31DeJzXsDiYzXlSVdUPGt5dyzMe53N4bi6/24/z9qIarp21srISZ1dX83N4Om2Zww139YX8mQ0aNsHlhnty33C+V1V1Xf7uWufxY4kbCAAAAAQ00AAAABDQQAMAAEBAAw0AAAABDTQAAAAENNAAAAAQ0EADAABAQAMNAAAAAQ00AAAABDTQAAAAEBilweXl5XjQuXVzD6uYR0rf93F2bm52n63ruji7sLAQZwer+btrMRjkz63ls3WTh1PNyf18Sw0zyrbMtdNFy+eryj/fO9tLifxfDdkXNLyPlnf3aHzPs9IyfwaDR/473TVr1sTZTZs2zbCSv5udO3fG2dEovgo03QVarE5X42zLnDpy5MhMxq2q2rJlS5ydTqdxdr4bN9WR6iYzOlxPA7Pac1ve26TPn2/rXtfy+cbjhvkzyOf8pGH+HM8rqOHSMM6uTvL30cL5Pnut++vp7pG/rQAAAMCjgAYaAAAAAhpoAAAACGigAQAAIKCBBgAAgIAGGgAAAAIaaAAAAAhooAEAACCggQYAAICABhoAAAACozQ4XVyOBx2uXZ9X0M2mhx/1eXbakN2+YVOcvXeaj1tVtXmSZ7esXxtn16yciLN9nz+MTcuNHzB0ZNDNZNy+y+vturyGUcNjaHm+86uz+35reRAv/aaaB8P8ub2yzz/f7wzyxfEDDet5rmHxtzyH6XQ2a6Mbr86khkE1PLTG5TkY5O+5bzgPVifjtkJCGzfne/zOc86eSQ2nwuXf8q1xdnFxMc4eP3L04ZRzUtOGud0yp+YW8r1u/cZ1cbaqan7Lxji7vJzfocaT43G2ZV+azOZorXHTUdWw1zRkW/a7lu25JTtuuDc0XEe+OnbTc8tfyKSlkEl+Dq++JN9TWp7buKHclnvctOV8H+fPYdBQ76ClGWk9iBu07CktWt7Ho4FfoAEAACCggQYAAICABhoAAAACGmgAAAAIaKABAAAgoIEGAACAgAYaAAAAAhpoAAAACGigAQAAIKCBBgAAgMAoDZ44cSIedOu2bXG2a+jhl1fyGo7UkTjb4vjx43G267qmsRcW5uPsaBS/uhoM8mfc932c7bo82zZu23NrGHg24z7Gtb2P/D23+MGWChrm2qyyjzYt77j1ObTsPy2WV1fi7Hdc/bQ4+x8++kdx9lf+6Rvi7Dfa9u3b4+zKSv4sX/mjr4iz//eH/jTOnvj2/HxfWFiIs9sa7iOXXnppnK2qOnDXV+Ls7t2742zLmplMJnH2sWw6ncbZvp/Nfjft8xpa7znThrO1peZJ5TW3GI9n8ywahm0ybsi2zLWGaOPZ+ti+z17zXVfH2Q//+Ufj7NXXXPVwynlQfoEGAACAgAYaAAAAAhpoAAAACGigAQAAIKCBBgAAgIAGGgAAAAIaaAAAAAhooAEAACCggQYAAICABhoAAAACGmgAAAAIjNLgiRMn4kH76TTPDvNsN5/3+4eG++Jsjfo4uqaGcXbT/Nq8hqraPOri7HCymtfR5Z9vMp3E2a7huVWfZ9fO6GudSf54q2+ot8W0z+f7sM8fRGu9XcOcaLGcf7yWKVFdNYQneRHTlr2qpeDJeDbjtqzPhnFbllzrXOuH+ejDYb6/Htt3f5z9q89+Mc6+6h//TJw9+9vOjrPfaKP5uTg7bphX//xN/zzOfs/PPi/O7r5/d5zdtmV7nD12+EicXVpZibNVVWs2b8jHvivfExYa9udpw+JdGeTvuUVfsxl30LA/jyvPnmjYkxqO4RpPW/bGtn10tSHfN2QHDedli+l4NneM2YxaNcjboeoneRWrDc93PMnXUd90h2s8s5vSuZZn/ME//4s4+8KX/P04e+4V58bZk/ELNAAAAAQ00AAAABDQQAMAAEBAAw0AAAABDTQAAAAENNAAAAAQ0EADAABAQAMNAAAAAQ00AAAABDTQAAAAEBilwb7v40FXV1fj7FycrKphS3g21q9fH2c3bdrUNPZw+URrOdm4w/zBdV2XZ2saZ6fTPNtSb8u87Cv/bC3jtmh5vpNudt9vtdTRNu5Mhm3SNCcaspPJJM52DfN9VvW2rY1cy1quaptr69ati7Ot+2vqnHPOibPnnXfeTGo4FXbu3Blnjxw5Emd37Djj4ZRzUktLS3F22HAZ2LhxY5zdsWNHnK2q2t+Qb5mv3eL+ODurPaFl3baM27J/9NN8zz0d9v1p37Y3tpg27NJte/Rs7jrT6Wzm2nRGc3iW7+6xrGWuDWbUw7W857PPPvuU/X39Ag0AAAABDTQAAAAENNAAAAAQ0EADAABAQAMNAAAAAQ00AAAABDTQAAAAENBAAwAAQEADDQAAAAENNAAAAARGcXAUR+vgwYNxduvWrXG2hnl03bp1ebjBpu15vffsva9p7L7Ls9PpNM52XT7wYJB/p9KNGsZteHnTlgfRYNr3Mxl3VlreW9/42Vrec8tca9FS87TPa2ipdzKZxNmWeict7y5OVvXVsDZaamjZe/JoVTVt2zUej+PsgUNHGivJXHLJJXH2iiuumEkNp8I555wTZ3fu3Blnu242++ji0vE4+5WD98TZNWvWxNltO7fH2aqq4W3zeR1b1sfZ4/v3xdnlSb5mlrp89Q6H+crtBy1nVcMdYyXfn1v2pT6/+jbtudOGpdF6Zk8b1l3f8NvYZJA/45Y7yXCaz8uGYWvY8D66huyk4RBseXVrpvk6mjSMuzhouD+13uEa1nPX5XNtteG+1WJuzUKcffxll56yv69foAEAACCggQYAAICABhoAAAACGmgAAAAIaKABAAAgoIEGAACAgAYaAAAAAhpoAAAACGigAQAAIKCBBgAAgMAoDQ6Hw3jQyWQSZ48fPx5na00e7bpuJtnLLrsszu7bty/OVlWt3H1XnO37fibZFi3PraWG02Hc08Gjrd5ZmtV8b8lOp9M82zdkG8atasnm+oap1rqftJwHR48ejbObNm1qqiN19tlnx9k1axoOpW+wM844I862vdPZzMH5+fk4u7KyMpMatmzZ0pR/4hOfGGcPHDgQZ//qi1+Os4NB/jtIS7bpDlWzOatmdQaeDvenWZrV/bfpbhYn2zza7kWzqneWz2FWd/BZ7duPf/zj42xLD3cyfoEGAACAgAYaAAAAAhpoAAAACGigAQAAIKCBBgAAgIAGGgAAAAIaaAAAAAhooAEAACCggQYAAICABhoAAAACozTYzeW99rgmcfbIscNxtrbn0UOHDsXZ0Sh+DLV58+Y4u3Xr1jhbVbXvK/fE2ZWVlTg76JvKiHUt4zZk5/rZFDzXUHDfUnCDaT/Nww0l9I3PbLUh2zU8t67yz9eS7acN873L95/qGt5HQ3Y8yWuYThtqGM7mO8++y7PdoK2GyTR/FouLi3H28suviLOfrE/F2UsvfXycHQ4bHtw32Pat2+Jsy/4xt27+4ZRzUqOFfNxxwz7area73fHjx+NsVdWZZ5wRZzdv3Bhnz7ng/Dh7zz35vaHlmBi3hFuOtZajapjfzfrBMM+Ox3G25fzrJ/l+0HrNGXQtY+eDDyoft2vI1jB/H4OGM6XpvGwwbDoE8+yo4bN1+SOrcZfva5PGa0NfDe+uYY3WXMvNM/ekK58cZ7duz8/Fk/ELNAAAAAQ00AAAABDQQAMAAEBAAw0AAAABDTQAAAAENNAAAAAQ0EADAABAQAMNAAAAAQ00AAAABDTQAAAAEBilwbm5uXjQ6XSaV9B1ebbByspKnF1eXo6zGzdujLObNm2Ks1VVexueW8sz7vs+zg4G+Xcqs3lzVV03m+91Wp5DS7ZF1zDfB32ebVpzfE3Le255xrPKDoaP/HeeLXtEVdV4kn++lv11y5YtTXWktm/ZHmePLR6bSQ2nwpo1a2Yz8Gg2O/9wOIyzszrTJpNJnK2qmjb8BnHeeefF2aN33R1n77vvvji7tJTP15ZnvH4ufw4t487ujpGP3FTvjO6zVVV9zaaOmWVb7pIt96LG8yc1bbhvzep+2HRvqDw7Ho+b6uiH+edr2TMf//jHx9mP13+Os09/+tPj7IYNG+LsyTzytzEAAAB4FNBAAwAAQEADDQAAAAENNAAAAAQ00AAAABDQQAMAAEBAAw0AAAABDTQAAAAENNAAAAAQ0EADAABAYBQHF+bjQfu+n0n27r274+y5550fZ++59444e/Y5Z8XZpzz1yXG2qurJF14YZz/zmc/E2Ts//KE4Oxp2cXahn8TZFt10OpNxa5rPtS5/DE1zuGXcQcNj6AZ5DVVVw65hjVbD2E1zouU9j/NRp6szqaFreGazMm15Fw36fhhnu77te9fVaT4nnnDpZXH2Tf/Hm+PssZWjcXa5luLssDt9v4Oen5/RmT2fb2L3HLkvzp67Y2c+7uF74mzXsO9PJrM506qqduzYEWfvOzPPbtt3Zpxd7lfi7MGDB+Psmi6+Sla1nK0tB2ZDdNjnZ0TXcEYMGoqYNt5zJg1rtMWw4c7X8jomTedlnl0dzOYZj1vmZUO9iw3PYTJsuWvl0eZ9reFZnHFWvv+843feHme/3NCX7Tgr3y/H4/wueTKn7+kPAAAApxENNAAAAAQ00AAAABDQQAMAAEBAAw0AAAABDTQAAAAENNAAAAAQ0EADAABAQAMNAAAAAQ00AAAABEZpcH5+fiYFDAZ5Dz+ZTGZSw9q1a2cy7sLCQlN+sH59nD3zzDPj7MrZZ8fZw4cPx9nlI4fi7HA4jLNzg7k423VdnK2G6OlgMMgLntXaeKzr+34m47bsay1zuG8YdzqdxtmW9bmyshJnq6qe8IRvibMveMEL4uy/rd9tqmMWZjV/ToXRKD7e28znc3BWz2duLj8japKvg1bT1XzfXd9wvl944YVxdjwex9n9g/xZtOwfk6NH4myLrs/3paZxG/bclr182jDdW8adpa7Li57Vczsd9tFZ1dDyzFrO4WFDvU335KraddFFcfbKK6+Ms++tP4yzmzdvjrOrq6txtvVZPJTTYwUDAADAaU4DDQAAAAENNAAAAAQ00AAAABDQQAMAAEBAAw0AAAABDTQAAAAENNAAAAAQ0EADAABAQAMNAAAAgVEanI76eNC5ubk4O51O42wNuzh6/NjBOLt+69Y4u3zsQJxdt2ZtnK2qmqxbibNnnnt2nF28/OI4e+DWW+PsoeP5Mx4M8ve8tSE7HA7j7Gg8ibOTSZ6dlXG+5GpaDeuoqrqm787yQoZdnp1OV+PsoBvPJDvs8uc27WczL6vyfW0yyN9by966dm2+V1122WVxtqrqWc+6Js5+3wteEmdPLB+Nsy3Poq987Q+7/N19o43m83O4xaRhzXQNz2ff0r44e8a6M+Ls3qP3xdmWequqGpZjU/ZbvjVfY5u3bIyzB3ZsjrNbbtseZ+/+3Gfj7H335e9jy2o+1wYND7ib5vtz3+dn2mDScGg3Gg9mcydZGjU8t4b1MW2Y7y3PeLlh4EnLuA33nKZ6W7aUhuyhg/n5d/bZeb9QVfWd3/ldcfblP/Ljcfbw+HicHaxteM8r+V3yVPILNAAAAAQ00AAAABDQQAMAAEBAAw0AAAABDTQAAAAENNAAAAAQ0EADAABAQAMNAAAAAQ00AAAABDTQAAAAEBilwfvvvz8edP369XF23bp1cXY4HMbZ0Wg23w10XRdnWz5bVdVcw+ebn5+Ps+OD98XZlZWVOLt2PI2zu3fvjrPDuTVxlq9qmZezHLulisGgYY3mU21mmp5xP5saWp5Zyz78Ld/yLXH2Gc94RpytqrrssnzsFtNpPin6Pn8hfcPLaxn3G21ubm4m4w6H+WduWTPj1eWHU85JtZzD4/G4aeyu4bwcjeLrVq1puAucddZZcXbdJK+3xeDYkTi7tLQUZ/u9++Nsy34wN8jXRsu4LRv/6bJ3NO2NDdlpw3E5mUzi7OokH3h1dTXPNvye2DInpnMLcbbFrl274mzrmf2UpzyltZxIy3tu0dIbnsq7sl+gAQAAIKCBBgAAgIAGGgAAAAIaaAAAAAhooAEAACCggQYAAICABhoAAAACGmgAAAAIaKABAAAgoIEGAACAwCgN3nnnnfGg5557bpxdu3ZtnO26Ls7WsI+jy6tH4uz8+k1xdun4gThbVVXD+HXU2oWNcfayK78tzm4658w4u/ecs+Lswhe25ON++Y44u7S0FGcHtZpn52bz3dJ0Oo2z44bvt/o+n+9VVat9PnZLzX3la3Q8HsfZwSivdzgc5uM2PLaWZ7y6lM+1c845J86ecfbZcXbHjh1x9lu+9Yo4e+WVV8bZqqqLL7w0zh48si/ODubyGiYNc3g6Xo6zTetuXR49FVZWVuJsy+c4fCI/L1ts2ZyfaUf7o3F2wzAf9+D0YJytqppvOLPnGrL9NN8bR/P5QjhjZ74nTLp8zdTR/K4zN5/vz0e/tDvOHjmSz8vV5XyNLy3l66i6/L1VtZ3Z08ls9qVJN5u7zuI0f88td4Ej3XycXR7n72655YxoyG5cuyHOtrjmu66Ns894xjOaxn7ipZfH2UMnDsfZbpjfD0cN87KhM2zrI0/CL9AAAAAQ0EADAABAQAMNAAAAAQ00AAAABDTQAAAAENBAAwAAQEADDQAAAAENNAAAAAQ00AAAABDQQAMAAEBglAZXV1fjQVuyg0Hew/d935CNo9V1XZxdXTocZ+fWb86LqKqVo1+Jsy01z2/ZEmfPaRh3YTyJsy3OWLMuzt5///1x9tDePXH2+PHjcXY8HsfZ0Sheck1a1kar6XQaZycNZSwtLcXZuw7ka6PJaC6Otqy5cx63K84+/vGPj7MXXXZZnL3iiivi7IW78hp2bNsZZ6uqji3me2Y/yidQy9nRMocfK/bu3Rtnjxw5Emc/d+vn4+zOnflc2XXh4+Ls+vXr4+yBdQfi7Nb5rXG2qmp5ejDOtuwfw4aNdH5+Ps6OurVx9qyzzoqzwyOH4uymTZvi7MFRfhe455574uzRw/metHv37jjbTWZzJ6qq6qvl/ptnZ7U3rk7ycVvuUJPh8OGUc1LDhnFb7nEt6+iCCy6Is8973vPi7K5dl8fZqqpDJ/L10XSnHeR7YMt+2dLwncq7sl+gAQAAIKCBBgAAgIAGGgAAAAIaaAAAAAhooAEAACCggQYAAICABhoAAAACGmgAAAAIaKABAAAgoIEGAACAwCgNbtu0LR60JTuY5j38oKHf76uLs9WQ7Yf5qMsrxxpqqJrfuCHOri4dzgfuVvMatp0ZZ7dMGp7x+k1xdOOadXF285l5vYd37Iiz999/f5zdu3dvnF1aWoqz4+kkz/Z5tqpq0uXvblzTPDsZx9nV8XKcXTu/kGfXro2z6zZvmcm4T3vms+LsJZdckme/9Vvj7K6LLo6zXeXPdzJejLNVVZP4lKmaTvO51vX5uINq2Lj7hoL7hiK+wW677bY4+4UvfCHOfvQTH4uzT3jCE+Ls4Npr4uy5554bZ1vW7ZHpkThbVbWwkJ9ri6uH4mw/yudrP8zvReNBw11nPn9u535rvodtPXEizl76+HwPu+uuu+Ls7Z//YpztN+bPYf/9e+Ls8ePH4+xX8/ndoW/Yl47XbPawlUE+h1cqv7/suT+/b3UN95zzLrgwzu7cuTPOvugHfijOPvnJT46zW3ZcEGf7SVsvsrowF2dbzuw82Xa+d30+csvaOBm/QAMAAEBAAw0AAAABDTQAAAAENNAAAAAQ0EADAABAQAMNAAAAAQ00AAAABDTQAAAAENBAAwAAQEADDQAAAIFRGtyxY0c8aEt2ZWUlzk6n0zhb3TCPdl2ezStoNl09FmcHazbE2X5l/8Mp56TWbtmS19D3cXb9uefG2Y0bN8bZpYZ6zzjjjJnUcODAgTh78ODhOHv8+PE4W1W1PMmzLetudXU1zrbMiUsvvTTOrlu3Ls5uPfOsONvynq+44oo4e84558TZs87K6+1qIc4eO5rPy35N/t7+ppJZaDoPvgndeeedcfYLX/hCnG1Z46dDdjDIfydo2ZOqqk6sHIyza+e3xNmlGZ3ZLXedJsP8vrVmzZo4Oxjme9jmzZvj7NOe9rQ4e/7558fZ/ffvibNHjx6Ns1VV73v/H8bZpjvtXP7uWvbcwXA2627nzp1xtuW8/O7nf0+c3bp1a5x98pOfHGe37Lggzvb9iTjbsj6rqrppPn9a9teWq0DrXvxI8As0AAAABDTQAAAAENBAAwAAQEADDQAAAAENNAAAAAQ00AAAABDQQAMAAEBAAw0AAAABDTQAAAAENNAAAAAQ0EADAABAYJQGx+NxPOiXv/zlOHv+uefF2cEg7/e7bhhnq7o42Xd9nm2ooKpq2vBXTCbH42w3vz7O9itH42zNzcXRdedszsc9c2cc3bS8nI97PH9mi4cP5zXceXuc3bNnT5z9zGf+Ms7evvf+OFtVdWJ5Kc52Xb4+lqaTOLvSTePsk6/89ji7bdu2OHvFt10ZZ3ft2hVn792zP87u2LEjzq5fvz3OLh89FGen0/xddOM8W1U16PJ9ouvzuTZtyPZ9vrcOurVxdtIw37/RVldX42zL2bpty9Y4u3Xzlji7ZctssoOWg3jSNren+RSs4+NDcXbN/JY4u7Sa7/3TQf75JpN8bq827OVzG/L1VZP85a3dtjHPrtsUZ3esia/JNbc5/2ybTpyIs1VVa/5zvu6ON9x1pl2+9sfTvA9o2lMazuzHXXRZnL344ovj7OUXPz7ObtqUz58tOy6Is5OD98TZ5VFDD9Bwh6uqqtFCHB0O815r0PCbbcuZ3XAVaO7LHopfoAEAACCggQYAAICABhoAAAACGmgAAAAIaKABAAAgoIEGAACAgAYaAAAAAhpoAAAACGigAQAAIKCBBgAAgMAoDS4sLMSDjsfjh1XMqdR13SNdwmmjXzkaZ7v5jbOpYXw8D08meXZlJY6OT5yIs0eP5s+sZb4PBvl3Vhs35u9ix44dcbaq6sDhQ3H22LFjcXY6ncbZljXa8txaati/f3+c7fs+zj7hymfH2RZLRw7G2cm0YR01aHkOVVV9teUfTU7nc2bdunVxdvv27XF2ebwcZ1v2pabsGXl23HBGtGrZa1r2sJVJvi/Nz+XvrsXRE7fG2eFwmA88iq+dVfN5du2mTXH26J78+bacEfvv3xtnl5aW4mxV1bZt2+Jsy339yFK+nlca1tKkYb5v2bIlzracPyca7nxPeOZ1cbbF9NC9cbbbnM/h1UNH4uyk5U5dVfPr8nXXsq/VaXAXOJVntl+gAQAAIKCBBgAAgIAGGgAAAAIaaAAAAAhooAEAACCggQYAAICABhoAAAACGmgAAAAIaKABAAAgoIEGAACAwCgNbtqwOR50Mpk0lJD38F03jLP9oGuoIddXw7iNJQwnDX9B3zDwZJwPu3i4YdyVONptWJ+PG8/KqsP33JRnD+yPs/v359lDBw/E2YMHD8bZcZ+/t3Ub18XZqqp+kE+gjRvzd9e3TOE+r+H+e++NsycO53P4ZT/xj+Nsi7++6RNxdmFhIc6ubF+Ks2vWrImzo1HDomvZe6pq2hY/DTw2vlc+a8cZcfbo4UNxdtLnb/Rx554XZ9fOrY2zyyeW4+zc3FycnU7bZuugYap0Xb45Thr2xhOr98fZfnIkzq5fd3GcbbHvzo/G2T233hFnb7vttjj7lbvuirMrK/k9p2VeLi/n2aqq48cXG9L5xJxbmJ9JdjLJ5/DS8RNx9tff9vY42+J/+tHr4mxLj/Pru/63OPusZz0rzq7bsjXOXnnVNXG2qmpl+Xge7vI9s5/mc6Jlv6xpnm0a9yQeGzcFAAAAmDENNAAAAAQ00AAAABDQQAMAAEBAAw0AAAABDTQAAAAENNAAAAAQ0EADAABAQAMNAAAAAQ00AAAABEZpcGFhIR50Op0+rGJOpu/7mWRbdF03k3FbtXy+bkbvoxYX4+jBu2+Ks+PJSpzdfPmT4uxj2a/+9z/Slv+375lRJY8u/+P1L42zo1G8Xda/P/ffx9kLL7wwzl79rGvj7MatW+PsZHk5zrbqRvmeeTrs26fJFv93tmbNmji7du3aODu/Jr8L7NixI85u2LAhzg4Gs/nuv33cGZ2tDSaTSR5uuAvcf+CWODs6cSzObr3gqjhbF+TR+u6G7GngF364reBf//0PzqiSR5f/9SdeFGdbzuz9q4fj7NLSUpzds2dPnL3llnzNPeHvfVucPbhnd5ytqtqy46w429LvjRv2quFwGGerHpk7hl+gAQAAIKCBBgAAgIAGGgAAAAIaaAAAAAhooAEAACCggQYAAICABhoAAAACGmgAAAAIaKABAAAgoIEGAACAwCgNdl0XDzo3NxdnB4O8h+/7fibZls9WNW3INupaxs4/X63f0jDu8Tg52XMwzh687+44e+Tw/jj7f3/0n8bZ3XffGWf37t0bZ4fT/F2srKzE2X4yjrOve8c742xV1Rtf/oI4u7q6Gmcn3Wy+kxsP8z2lZT2PBw1rfzKJo3vvuy/O7r///jh7x533xNnLLrsszn7nd35nnN20dWucraqaTBr2qgYte3zbuDMZ9huuG+QfZN36NXF27Zr1cXZhbj7OTqf5+TcaxVeXJm13gaqu2vKpvs/3muk0PydqNT9/VpcW4+w9t98RZ9/3F/8szi7tO5Bnl5bi7M2f+nScbbmjzs0txNk3vusDcbaq6pev/96mfKplzs/qDt6y9ufn83GHDff1rWvyd7fY5+NOlvI79d1//bk4Oxzm723LuvyzVVU9efuOvI5Rfh6cOJLfdebWrYuz44Z7Z+se/1D8Ag0AAAABDTQAAAAENNAAAAAQ0EADAABAQAMNAAAAAQ00AAAABDTQAAAAENBAAwAAQEADDQAAAAENNAAAAARGj3QBgy7v4ac1nWElmeFgGGenfVu9fd/H2a7rGkY+HienX/lKnL399tvj7O7du+Ps4UP74uyePXvi7JGjh+Jsy/MdDPLswsJCnK1pvjz/t1e+OB+3qmnlDwb5Gp20rOdpvj76Ll93LetoVsbj8Uyy99xzT5zdu3dvnG15x09/+tPjbFXV5m1nxtmmd9ew7r4ZtbzTNWvWxNmNGzfG2fn5+Tg7K21nZZtZndkte8Lq6mqePXEizt55551x9u5bb42zhw4dirN1LK93MpnE2W3btsXZlve2bt2GOPum/+XH42xV1ZH5I3G25Ww9HbS8u5b30ZLdtGlTnF27dm2cXezn4mxLvTfddFOcPXbsWJytqjq6nM+fq666Ks4Oh/k9ruX8mk7yelvGPelYp2wkAAAAeAzTQAMAAEBAAw0AAAABDTQAAAAENNAAAAAQ0EADAABAQAMNAAAAAQ00AAAABDTQAAAAENBAAwAAQKDr+75/pIsAAACA051foAEAACCggQYAAICABhoAAAACGmgAAAAIaKABAAAgoIEGAACAgAYaAAAAAhpoAAAACGigAQAAIPD/AASojFlrB7jZAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "Failed to read frame 155.\n", "Failed to read frame 155.\n", "Failed to read frame 156.\n", "Failed to read frame 157.\n", "Failed to read frame 158.\n", "Failed to read frame 159.\n", "Failed to read frame 160.\n", "Failed to read frame 161.\n", "Failed to read frame 162.\n", "Failed to read frame 163.\n", "Failed to read frame 164.\n", "Failed to read frame 165.\n", "Failed to read frame 166.\n", "Failed to read frame 167.\n", "Failed to read frame 168.\n", "Failed to read frame 169.\n", "Failed to read frame 170.\n", "Failed to read frame 171.\n", "Failed to read frame 172.\n", "Failed to read frame 173.\n", "Failed to read frame 174.\n", "Failed to read frame 175.\n", "Failed to read frame 176.\n", "Failed to read frame 177.\n", "Failed to read frame 178.\n", "Failed to read frame 179.\n", "Failed to read frame 180.\n", "Failed to read frame 181.\n", "Failed to read frame 182.\n", "Failed to read frame 183.\n", "Failed to read frame 184.\n", "Failed to read frame 185.\n", "Failed to read frame 186.\n", "Failed to read frame 187.\n", "Failed to read frame 188.\n", "Failed to read frame 189.\n", "Failed to read frame 190.\n", "Failed to read frame 191.\n", "Failed to read frame 192.\n", "Failed to read frame 193.\n", "Failed to read frame 194.\n", "Failed to read frame 195.\n", "Failed to read frame 196.\n", "Failed to read frame 197.\n", "Failed to read frame 198.\n", "Failed to read frame 199.\n", "Failed to read frame 200.\n", "Failed to read frame 201.\n", "Failed to read frame 202.\n", "Failed to read frame 203.\n", "Failed to read frame 204.\n", "Failed to read frame 205.\n", "Failed to read frame 206.\n", "Failed to read frame 207.\n", "Failed to read frame 208.\n", "Failed to read frame 209.\n", "Failed to read frame 210.\n", "Failed to read frame 211.\n", "Failed to read frame 212.\n", "Failed to read frame 213.\n", "Failed to read frame 214.\n", "Failed to read frame 215.\n", "Failed to read frame 216.\n", "Failed to read frame 217.\n", "Failed to read frame 218.\n", "Failed to read frame 219.\n", "Failed to read frame 220.\n", "Failed to read frame 221.\n", "Failed to read frame 222.\n", "Failed to read frame 223.\n", "Failed to read frame 224.\n", "Failed to read frame 225.\n", "Failed to read frame 226.\n", "Failed to read frame 227.\n", "Failed to read frame 228.\n", "Failed to read frame 229.\n", "Failed to read frame 230.\n", "Failed to read frame 231.\n", "Failed to read frame 232.\n", "Failed to read frame 233.\n", "Failed to read frame 234.\n", "Failed to read frame 235.\n", "Failed to read frame 236.\n", "Failed to read frame 237.\n", "Failed to read frame 238.\n", "Failed to read frame 239.\n", "Failed to read frame 240.\n", "Failed to read frame 241.\n", "Failed to read frame 242.\n", "Failed to read frame 243.\n", "Failed to read frame 244.\n", "Failed to read frame 245.\n", "Failed to read frame 246.\n", "Failed to read frame 247.\n", "Failed to read frame 248.\n", "Failed to read frame 249.\n", "Failed to read frame 250.\n", "Failed to read frame 251.\n", "Failed to read frame 252.\n", "Failed to read frame 253.\n", "Failed to read frame 254.\n", "Failed to read frame 255.\n", "Failed to read frame 256.\n", "Saved rotation data to /home/davic/projects/IACV_project/bowling-analysis/notebook/spin/intermediate_data/Rotation_data_TEST_5.csv\n" ] } ], "source": [ "from pathlib import Path\n", "import numpy as np\n", "import cv2\n", "import pandas as pd\n", "from matplotlib import pyplot as plt\n", "import time\n", "from IPython.display import clear_output, display\n", "\n", "# --- Configurations ---\n", "VIDEO_NUMBER = \"5\"\n", "PROJECT_ROOT = Path().resolve().parent.parent\n", "INPUT_VIDEO_PATH = str(\n", " PROJECT_ROOT\n", " / \"data\"\n", " / f\"recording_{VIDEO_NUMBER}\"\n", " / f\"Recording_{VIDEO_NUMBER}.mp4\"\n", ")\n", "INPUT_CSV_PATH = str(\n", " PROJECT_ROOT\n", " / \"data\"\n", " / \"auxiliary_data\"\n", " / \"circle_positions\"\n", " / f\"Adjusted_positions_new_{VIDEO_NUMBER}.csv\"\n", ")\n", "\n", "# --- Load video and CSV ---\n", "cap = cv2.VideoCapture(INPUT_VIDEO_PATH)\n", "fps = cap.get(cv2.CAP_PROP_FPS)\n", "\n", "if not cap.isOpened():\n", " raise IOError(\"Error: Could not open video.\")\n", "\n", "df = pd.read_csv(INPUT_CSV_PATH)\n", "if df.empty:\n", " raise ValueError(\"Error: The CSV file is empty.\")\n", "\n", "frame_count = int(cap.get(cv2.CAP_PROP_FRAME_COUNT))\n", "offset = 2\n", "output_data = []\n", "\n", "\n", "# --- Process frame pairs ---\n", "for frame_number in range(frame_count - 1):\n", " cap.set(cv2.CAP_PROP_POS_FRAMES, frame_number)\n", " ret, frame1 = cap.read()\n", " if not ret:\n", " print(f\"Failed to read frame {frame_number}.\")\n", " continue\n", "\n", " ret, frame2 = cap.read()\n", " if not ret:\n", " print(f\"Failed to read frame {frame_number + 1}.\")\n", " continue\n", "\n", " # --- Get CSV data for current frame pair ---\n", " try:\n", " df_row1 = df.iloc[frame_number]\n", " df_row2 = df.iloc[frame_number + 1]\n", "\n", " if (\n", " df_row1[[\"x\", \"y\", \"radius\"]].isnull().any()\n", " or df_row2[[\"x\", \"y\", \"radius\"]].isnull().any()\n", " ):\n", " continue\n", " except IndexError:\n", " print(f\"No CSV data for frame {frame_number}.\")\n", " continue\n", "\n", " # --- Compute ball centers and ROI bounds ---\n", " ball_center1 = np.array([int(df_row1[\"x\"]), int(df_row1[\"y\"])])\n", " ball_center2 = np.array([int(df_row2[\"x\"]), int(df_row2[\"y\"])])\n", " ball_radius = int(df_row1[\"radius\"])\n", "\n", " def roi_bounds(center, radius, frame_shape):\n", " x_min = max(center[0] - radius - offset, 0)\n", " x_max = min(center[0] + radius + offset, frame_shape[1])\n", " y_min = max(center[1] - radius - offset, 0)\n", " y_max = min(center[1] + radius + offset, frame_shape[0])\n", " return x_min, x_max, y_min, y_max\n", "\n", " x_min1, x_max1, y_min1, y_max1 = roi_bounds(ball_center1, ball_radius, frame1.shape)\n", " x_min2, x_max2, y_min2, y_max2 = roi_bounds(ball_center2, ball_radius, frame2.shape)\n", "\n", " roi1 = frame1[y_min1:y_max1, x_min1:x_max1]\n", " roi2 = frame2[y_min2:y_max2, x_min2:x_max2]\n", " gray1 = cv2.cvtColor(roi1, cv2.COLOR_BGR2GRAY)\n", " gray2 = cv2.cvtColor(roi2, cv2.COLOR_BGR2GRAY)\n", "\n", " center_roi1 = ball_center1 - [x_min1, y_min1]\n", " center_roi2 = ball_center2 - [x_min2, y_min2]\n", "\n", " # --- Mask for feature detection ---\n", " mask = np.zeros_like(gray1)\n", " cv2.circle(mask, tuple(center_roi1.astype(int)), int(ball_radius * 0.8), 255, -1)\n", "\n", " feature_params = dict(\n", " maxCorners=100, qualityLevel=0.001, minDistance=0, blockSize=3\n", " )\n", "\n", " p0 = cv2.goodFeaturesToTrack(gray1, mask=mask, **feature_params)\n", " if p0 is None:\n", " raise ValueError(\"No features detected. Try relaxing parameters.\")\n", "\n", " # --- Optical flow tracking (forward + backward for validation) ---\n", " lk_params = dict(\n", " winSize=(21, 21),\n", " maxLevel=3,\n", " criteria=(cv2.TERM_CRITERIA_EPS | cv2.TERM_CRITERIA_COUNT, 15, 0.01),\n", " )\n", "\n", " # Forward tracking\n", " p1, status_forward, _ = cv2.calcOpticalFlowPyrLK(\n", " gray1, gray2, p0, None, **lk_params\n", " )\n", " # Backward tracking\n", " p0r, status_backward, _ = cv2.calcOpticalFlowPyrLK(\n", " gray2, gray1, p1, None, **lk_params\n", " )\n", "\n", " # Calculate forward-backward error\n", " fb_error = np.linalg.norm(p0 - p0r, axis=2)\n", "\n", " # Filtering criteria\n", " fb_threshold = 10.0 # Pixel threshold for forward-backward error\n", " movement_threshold = ball_radius # Reject points moving too far (could be tuned)\n", " low_movement_threshold = ball_radius / 5 # Reject points with moving too little\n", "\n", " good_pts = []\n", " old3d, new3d = [], []\n", "\n", " for i, (old_pt, new_pt, back_pt, s1, s2, err) in enumerate(\n", " zip(\n", " p0.reshape(-1, 2),\n", " p1.reshape(-1, 2),\n", " p0r.reshape(-1, 2),\n", " status_forward.ravel(),\n", " status_backward.ravel(),\n", " fb_error.ravel(),\n", " )\n", " ):\n", " if s1 and s2 and err < fb_threshold:\n", " displacement = np.linalg.norm(new_pt - old_pt)\n", " if (\n", " displacement < movement_threshold\n", " and displacement > low_movement_threshold\n", " ):\n", " ox, oy = old_pt - center_roi1\n", " nx, ny = new_pt - center_roi2\n", " oz = np.sqrt(max(ball_radius**2 - ox**2 - oy**2, 0))\n", " nz = np.sqrt(max(ball_radius**2 - nx**2 - ny**2, 0))\n", " old3d.append([ox, oy, oz])\n", " new3d.append([nx, ny, nz])\n", " good_pts.append((old_pt, new_pt)) # Save for visualization\n", "\n", " # --- At this point, old3d and new3d are already filtered; just convert them to arrays ---\n", " old3d = np.array(old3d)\n", " new3d = np.array(new3d)\n", "\n", " if old3d.shape[0] < 3:\n", " low_movement_threshold = (\n", " ball_radius / 10\n", " ) # Reject points with moving too little\n", "\n", " good_pts = []\n", " old3d, new3d = [], []\n", "\n", " for i, (old_pt, new_pt, back_pt, s1, s2, err) in enumerate(\n", " zip(\n", " p0.reshape(-1, 2),\n", " p1.reshape(-1, 2),\n", " p0r.reshape(-1, 2),\n", " status_forward.ravel(),\n", " status_backward.ravel(),\n", " fb_error.ravel(),\n", " )\n", " ):\n", " if s1 and s2 and err < fb_threshold:\n", " displacement = np.linalg.norm(new_pt - old_pt)\n", " if (\n", " displacement < movement_threshold\n", " and displacement > low_movement_threshold\n", " ):\n", " ox, oy = old_pt - center_roi1\n", " nx, ny = new_pt - center_roi2\n", " oz = np.sqrt(max(ball_radius**2 - ox**2 - oy**2, 0))\n", " nz = np.sqrt(max(ball_radius**2 - nx**2 - ny**2, 0))\n", " old3d.append([ox, oy, oz])\n", " new3d.append([nx, ny, nz])\n", " good_pts.append((old_pt, new_pt)) # Save for visualization\n", "\n", " # --- At this point, old3d and new3d are already filtered; just convert them to arrays ---\n", " old3d = np.array(old3d)\n", " new3d = np.array(new3d)\n", "\n", " if old3d.shape[0] < 3:\n", " print(\n", " f\"Frame {frame_number}: only {old3d.shape[0]} points after filtering, skipping rotation\"\n", " )\n", " continue\n", "\n", " # --- Kabsch algorithm for rotation ---\n", " P = old3d # - old3d.mean(axis=0)\n", " Q = new3d # - new3d.mean(axis=0)\n", " H = P.T @ Q\n", " U, _, Vt = np.linalg.svd(H)\n", " R = Vt.T @ U.T\n", " if np.linalg.det(R) < 0:\n", " Vt[-1, :] *= -1\n", " R = Vt.T @ U.T\n", "\n", " cos_theta = (np.trace(R) - 1) / 2\n", " cos_theta = np.clip(cos_theta, -1.0, 1.0)\n", " theta = np.arccos(cos_theta)\n", " if np.isclose(theta, 0):\n", " axis = np.array([0.0, 0.0, 1.0]) # default/fallback axis\n", " else:\n", " axis = np.array([R[2, 1] - R[1, 2], R[0, 2] - R[2, 0], R[1, 0] - R[0, 1]]) / (\n", " 2 * np.sin(theta)\n", " )\n", " axis /= np.linalg.norm(axis)\n", " print(f\"Rotation matrix:\\n{R}\")\n", " print(f\"θ = {theta:.3f} rad, which is {theta * 180 / np.pi:.1f}°\")\n", "\n", " # Save data for output\n", " output_data.append(\n", " {\n", " \"frame\": int(df_row1[\"frame\"]),\n", " \"x\": df_row1[\"x\"],\n", " \"y\": df_row1[\"y\"],\n", " \"radius\": df_row1[\"radius\"],\n", " \"x_axis\": axis[0],\n", " \"y_axis\": axis[1],\n", " \"z_axis\": axis[2],\n", " \"angle\": theta, # in radians; multiply by 180/np.pi if degrees preferred\n", " }\n", " )\n", "\n", " # --- Compute rotation poles ---\n", " pole3d = axis * ball_radius\n", " opp_pole3d = -pole3d\n", "\n", " pole2d = center_roi1[:2] - pole3d[:2]\n", " opp_pole2d = center_roi1[:2] - opp_pole3d[:2]\n", "\n", " # --- Visualization overlay ---\n", " mask = np.zeros_like(roi1)\n", " mask2 = np.zeros_like(roi2)\n", "\n", " for old_pt, new_pt in good_pts:\n", " c = tuple(old_pt.astype(int))\n", " a = tuple(new_pt.astype(int))\n", " cv2.circle(mask, c, 1, (0, 0, 255), -1)\n", " cv2.circle(mask2, a, 1, (255, 0, 0), -1)\n", " cv2.line(mask2, c, a, (0, 255, 0), 1)\n", "\n", " output_roi1 = cv2.addWeighted(roi1, 1, mask, 1, 0)\n", " output_rgb1 = cv2.cvtColor(output_roi1, cv2.COLOR_BGR2RGB)\n", "\n", " output_roi2 = cv2.addWeighted(roi2, 1, mask2, 1, 0)\n", " output_rgb2 = cv2.cvtColor(output_roi2, cv2.COLOR_BGR2RGB)\n", "\n", " cv2.circle(output_rgb1, tuple(center_roi1), ball_radius, (0, 255, 0), 1)\n", " cv2.circle(output_rgb2, tuple(center_roi2), ball_radius, (0, 255, 0), 1)\n", "\n", " cv2.circle(output_rgb1, tuple(pole2d.astype(int)), 2, (0, 0, 255), -1)\n", " # cv2.circle(output_rgb1, tuple(opp_pole2d.astype(int)), 2, (0, 0, 255), -1)\n", "\n", " print(f\"Frame {frame_number}: {len(good_pts)} points used\")\n", "\n", " # --- Display results ---\n", " clear_output(wait=True)\n", " fig, axes = plt.subplots(1, 2, figsize=(10, 5))\n", " axes[0].imshow(output_rgb1)\n", " axes[0].axis(\"off\")\n", " axes[0].set_title(f\"Frame {frame_number}\")\n", " axes[1].imshow(output_rgb2)\n", " axes[1].axis(\"off\")\n", " axes[1].set_title(f\"Frame {frame_number + 1}\")\n", " plt.tight_layout()\n", " plt.show()\n", " display(plt.gcf())\n", " plt.close()\n", " time.sleep(1 / fps)\n", "\n", "output_df = pd.DataFrame(output_data)\n", "OUTPUT_CSV_PATH = str(\n", " PROJECT_ROOT\n", " / \"notebook\"\n", " / \"spin\"\n", " / \"intermediate_data\"\n", " / f\"Rotation_data_TEST_{VIDEO_NUMBER}.csv\"\n", ")\n", "output_df.to_csv(OUTPUT_CSV_PATH, index=False)\n", "print(f\"Saved rotation data to {OUTPUT_CSV_PATH}\")" ] } ], "metadata": { "kernelspec": { "display_name": ".venv", "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.12.3" } }, "nbformat": 4, "nbformat_minor": 5 }