{"cells": [{"cell_type": "markdown", "metadata": {}, "source": ["# 2A.ml - Bayesian models with Python\n", "\n", "Mod\u00e8les de m\u00e9langes de lois. Statistiques bay\u00e9siennes. *bayespy*, *scikit-learn*."]}, {"cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [{"data": {"text/html": ["
run previous cell, wait for 2 seconds
\n", ""], "text/plain": [""]}, "execution_count": 2, "metadata": {}, "output_type": "execute_result"}], "source": ["from jyquickhelper import add_notebook_menu\n", "add_notebook_menu()"]}, {"cell_type": "markdown", "metadata": {}, "source": ["You can read [Probabilistic Programming and Bayesian Methods for Hackers](http://nbviewer.jupyter.org/github/CamDavidsonPilon/Probabilistic-Programming-and-Bayesian-Methods-for-Hackers/blob/master/Prologue/Prologue.ipynb). Results might be different between examples. The example used is the same but the default parameters the optimisation uses are different."]}, {"cell_type": "markdown", "metadata": {}, "source": ["We try different python model to deal with a Bayesian problem: a Gaussian Mixture. We will use the following example."]}, {"cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": ["%matplotlib inline"]}, {"cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": ["import matplotlib.pyplot as plt"]}, {"cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [{"data": {"image/png": "iVBORw0KGgoAAAANSUhEUgAAAlIAAAGfCAYAAACUQTyEAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAAsTAAALEwEAmpwYAAAsoklEQVR4nO3db2wd13nn8d8j+qa5SrtmFlVR6NpuVKCQUVfrMCaSdgksEKUbZevEIeRukwAt0PaFUWDbTbIGs1TTbWygC3HBbZsALRYI+udNgsSurXKTqlslXbovVoC9oUKrWsfWooib2NcJqmDNtI2YhqLPviCvdEnO3D/nnpk5Z+b7AQxYl+TM0Yj3zjPPec5zzDknAAAAjO9Q1QMAAABIFYEUAACAJwIpAAAATwRSAAAAngikAAAAPBFIAQAAeAoSSJnZtJk9YWYvmNnzZvZTIY4LAAAQs9sCHecTkv7COfezZvY6SYcDHRcAACBaNmlDTjO7XdKzkn7U0d0TAAA0SIiM1DFJ1yT9sZndK+mSpA86577T/01m9pCkhyTpDW94w3133313gFMDAAAU69KlS99yzh3J+lqIjNSspKclzTnnnjGzT0j6e+fcf8r7mdnZWbe2tjbReQEAAMpgZpecc7NZXwtRbP6ypJedc8/s/vkJSW8JcFwAAICoTRxIOee+KeklMzu++9I7JH1l0uMCAADELtSqvV+T9OndFXtflfRLgY4LAAAQrSCBlHPuWUmZc4cAAAB1RWdzAAAATwRSAAAAngikAAAAPBFIAQAAeCKQAgAA8EQgBQAA4IlACgAAwBOBFAAAgCcCKQAAAE+htogBGmFlvavlC1f1ysamjk63tXDquOZnOlUPCwBQEQIpYEQr612dOXdFm1vbkqTuxqbOnLsiSQRTANBQTO0BI1q+cPVmENWzubWt5QtXKxoRAKBqBFLAiF7Z2BzrdQBA/RFIASM6Ot0e63UAQP0RSAEjWjh1XO3W1J7X2q0pLZw6XtGIgOFW1ruaW1rVscXzmlta1cp6t+ohAbVCsTkwol5BOav2kAoWSADFI5ACxjA/0+EGhGQMWiDB7zEQBlN7AFBTLJAAikcgBQA1xQIJoHgEUgBQUyyQAIpHjRQA1BQLJIDiEUgBQI2xQAIoFlN7AAAAngikAAAAPBFIAQAAeCKQAgAA8EQgBQAA4IlACgAAwBOBFAAAgCcCKQAAAE8EUgAAAJ4IpAAAADyxRQwARGxlvcteeUDECKQAIFIr612dOXdFm1vbkqTuxqbOnLsiSQRTQCQIpAAERxYljOULV28GUT2bW9tavnCV6wlEgkAKQFBkUcJ5ZWNzrNcBlI9icwBBDcqiYDxHp9tjvQ6gfARSAIIiixLOwqnjarem9rzWbk1p4dTxikYEYD8CKQBBkUUJZ36mo7OnT6gz3ZZJ6ky3dfb0CaZIgYhQIwUgqIVTx/fUSElkUSYxP9MhcAIiRiAFIKjeTZ9VewCagEAKQHBkUQA0BTVSAAAAngikAAAAPBFIAQAAeCKQAgAA8EQgBQAA4IlACgAAwBOBFAAAgCcCKQAAAE8EUgAAAJ4IpAAAADyxRQwANNDKepf9EIEACKQAoGFW1rs6c+6KNre2JUndjU2dOXdFkgimgDExtQcADbN84erNIKpnc2tbyxeuVjQiIF0EUgDQMK9sbI71OoB8BFIA0DBHp9tjvQ4gH4EUADTMwqnjarem9rzWbk1p4dTxikYEpIticwBomF5BOav2gMkRSAFAA83PdAicgACY2gMAAPAULJAysykzWzezPwt1TAAAgJiFzEh9UNLzAY8HAAAQtSA1UmZ2h6T7Jf1nSf8hxDEBIHZsswIgVLH5xyV9RNIPBDoeAESNbVYASAGm9szs3ZL+zjl3acj3PWRma2a2du3atUlPC6BhVta7mlta1bHF85pbWtXKerfS8bDNCgApTI3UnKQHzOxvJX1W0kkz+9T+b3LOfdI5N+ucmz1y5EiA0wJoil72p7uxKadb2Z8qgym2WQEgBQiknHNnnHN3OOfeJOn9kladcz8/8cgAYFeM2R+2WQEg0UcKQAJizP6wzQoAKXBnc+fcX0n6q5DHBICj0211M4KmKrM/bLMCQGKLGAAJWDh1fM8KOSmO7A/brAAgkAIQPbI/O+hbBcSHQApAEpqe/aFvFRAnis0BIAExrlwEQCAFAEmIceUiAAIpAEgCfauAOBFIAUAC6FsFxIlicwBIACsXgTgRSAFAIpq+chGIEVN7AAAAngikAAAAPBFIAQAAeCKQAgAA8EQgBQAA4IlACgAAwBOBFAAAgCcCKQAAAE8EUgAAAJ4IpAAAADwRSAEAAHgikAIAAPDEpsVAjpX1rpYvXNUrG5s6Ot3WwqnjbBgLANiDQArIsLLe1ZlzV7S5tS1J6m5s6sy5K5JEMAUAuImpPSDD8oWrN4Oons2tbS1fuFrRiAAAMSIjBWR4ZWNzrNfRDD7TvUwRj49rhpQQSAEZjk631c0Imo5OtysYDWLgM93LFPH4uGZIDVN7QIaFU8fVbk3tea3dmtLCqeMVjQhV85nuZYp4fMOu2cp6V3NLqzq2eF5zS6taWe9WMUzgJjJSQIbeky/TC+jxme5liniwrCm8QdeMbBViRCAF5Jif6fDhjJt8pnuLmCIuun6orPqkvKDo9nZLG5tbB77/6HR7YLaK9yqqwtQeAIzAZ7o39BRxL/jobmzK6VbwEWp6q4zj96blHn78cmZQZKbca0aGDzEikAKAEczPdHT29Al1ptsySZ3pts6ePjEwE+LzM4MUXXNV5PH3B2nbzmV+38b1rdxrlpfJYxEIqsTUHgCMyGe6N+QUcV7mpbuxqbml1Ymn4YrM+GQFaVmOTrdzr9nCqeN7pgMlFoGgemSkACARgzIvIabhisz4jBKMDQuK9mf4ptstvb51SB9+7FlW8KEyBFIAkIismqt+k07DFdn2Iy8YmzIba9pzfqaji4sn9bvve7P+6cZrevX6ViH1XMComNpDY9E9Ganpb8uRtRpQmmwarsi2H3nTcr41Y6zgQywIpNBI9KNBqnr1Q3NLq4V038+qTwrx0BE6SGMFH2JBIIVG4mkWqSur8DrkQ0fIwnu2cUIsqJFCI/E0i9SFbq2QJ9ZtbtjGCbEgI4VG4mkWdVBG9/1YHzrYxgmxIJBCI9GPBhhNzA8dbOOEGDC1h0Yqa1oE8NW/nUqVPZKYQgMGIyOFxuJpFrGKaVUpU2jAYARSwAD0mkIVYltVykMHkI9ACsgRU1YAzRJrgTeAgwikgH16WaisAlt6TaEMMRd4A9iLYnOgTy8Llbf9hkRWAMWjwBtIBxkpoE9Wbcp+ZAVQNAq8gXQQSKFRhhWPD8s2kRVAWepS4M2CDdQdgRQaY5Ti8bzaFGmn1xQ3AWB0LNhAE1AjhcYYZc+wvNqUj7/vzbq4eJIPf2AMse7TB4RERgqNMcqScmpTkIoUpsxo44AmIJBCY4y6pLwutSmor1SmzGjjgCZgag+NwZJy1EURU2ZF7O3Hew5NQEYKjcG0Heoi9JSZT4ZrlKlF3nNoAgIpNArTdqiD0FNm4+7tN07gxXsOdcfUHqJVxFQDUAehp8zGzXCxGg+4hYwUopRKMS1QhdBTZuNmuFiNB9xCIIUojTvVADRNyCmzhVPH9zy4SIMzXKzGA25hag9R4okXKM/8TEdnT59QZ7ot004X/7OnT+QGaqzGA24hI4Uo8cQLlGucDBer8YBbCKQQpXGnGgCUi9V4wI6Jp/bM7E4ze8rMvmJmz5nZB0MMDM027lQDAABVCJGRuiHpYefcl83sByRdMrMvOue+EuDYaDCeeAEAsZs4kHLOfUPSN3b//x/M7HlJHUkEUgAaJ4XNhAGEE7RGyszeJGlG0jMhjwsAKaD/GdA8wdofmNn3S3pS0oecc3+f8fWHzGzNzNauXbsW6rQAEA06fgPNEySQMrOWdoKoTzvnzmV9j3Puk865Wefc7JEjR0KcFgCiQv8zoHlCrNozSX8o6Xnn3O9MPiQASFNenzP6nwH1FSIjNSfpFySdNLNnd//7mQDHBYCk0PEbaJ4Qq/b+lyQLMBYASBodv4HmobM5AARE/zM/tI1AqgikAACVom0EUmbOudJPOjs769bW1ko/LwA0UezZnrml1cxNyqfM9JpzUY4ZzWJml5xzs1lfIyMFADWWQrYnrz3E9u6DfoxjBnqCNeQEAMQnhSaho7SHiG3MQA+BFADUWApNQrPaRmSJacxAD1N7wD6x15MA4zg63c6sPzpkpmOL56P4Hd/fNuKQ2c1pvX40NkWMCKSAPinUkwDjWDh1fM/vdE9s9Uf9bSP2vw+l8hqb8iCFcTG1B/QJWU+yst7V3NKqji2e19zSqlbWu6GGCYxsfqajs6dPqDPdlmlnJdx+sdUf7R9zZ7qts6dPFB7Q9AK47samnG4Fmbx3MQgZKaBPqHoSMluISX+259ji+czvia3+qIrGpoMepHjfIg8ZKdTeOJmhUJvOprBSCs3Exsr5UijMR3wIpFBr46bqQ206ywcyYsXGyvkIMuGDQAq1Nm5mKFRtBh/IiMX+jKykSuqPUkCQCR/USKHWfDJDIWozslZK7f9AZnUQipZXq3f29AldXDxZ8ejisP99+OB9HT31wjXelxgZgRRqLa+HTtGZof19cfZ/IFOMjjJQPD1Y1vvwyUtdMnQYC4EUam2UzFBRBmW2uMGhDNTqDcb7ECFQI4Vaq6ofzTDc4FAGavUG432IEMhIofaq6EczTFVTjohPkbVyVWZkU8D7ECEQSAEl6b9h3t5uqTVl2tq+tZ8YN7jmKbpWblitXtMRaCIEAimgBPtvmBubW2odMr3xcEsb17e4wTVUGTU6MWZkY0GgiRAIpAAV34og64a59ZrT4dfdpvXffGew8yAt1OhUj0ATk6LYHI1Xxkal3DCRhWJwIH0EUmi8MvbF44aJLHTSBtJHIIXGKyNbxA0TWWJtzwFgdNRIofHKWAJNUSvyUKMDpI1ACo1X1hJobpgAUD8EUmg8skUAAF8EUoD8s0VFt00AAMSNQArwVHRXaiBlPGSgKVi1B3gqo20CkKIyerMBsSAjBXiiySaQrYytb+qCzF36yEgBnmiyCWTjIWM0ZO7qgUAK8ESTTSAbDxmjoTygHgikAE90pQay8ZAxGjJ39UCNFDABmmwCB9GbbTRl7KqA4hFIAUDJmlBgzEPGcGXtqoBiEUgBQInoP4YeMnf1QCAFACWiNQD6kblLH4EUAJSIAmM0YWq3SQikgBx82KEIFBg3Wx2mdvls3Iv2B0AGGuWhKLQGaLbUe0fx2XgQgRQqsbLe1dzSqo4tntfc0mp0b8LUP+wQL/qPNVvqU7t8Nh7E1B5Kl0JqO/UPO8SNAuPmSn1ql8/Gg8hIoXQpPNGwxQWAIqQ+tctn40FkpFA6nyeasosbF04d18ITl7W17W6+1pqygR92FGAiBan9nqY23mFS7x1FE9GDCKRQunFT25VNBbohf+6zst7Vwp9c1tZrO9/U3djUwp9clhTPdCWQwrR6v9TGO6qUp3ZTDwSLYM4NuDsUZHZ21q2trZV+XsRh/4ejtPNEk1dwO7e0mhl4dabburh4spAx5p2zd979HxxvfvQL2tjcOvC90+2Wnv3YOwsZ46Tq9qSP4ap4L00itfGivszsknNuNutr1EihdOOuWsqb8utubBa22m/QNGPWct+sIGrQ61VjCXMzpVYonNp40UxM7aES46S286YCJRWW5h90Tin9LT3YpqSZylox5pPtzPqZ1Fe4oRnISCF6Watceopa7TfonD39T8VvPNzK/J6816s2zpN+7D2/MLoyVoz5ZDvzfubtdx9JeoUbmoFACtHrTQXm6d38Q97w+6cf8/Q/FX/sPfeoNWV7vt6aMn3sPfd4j6FIoy5hZgqwXspoBurT3iTvZ5564RrNSxE9pvZqqm6FxPMzHS1fuJqb5h+0ukcab4XJ/mv38z95l5681B243De1lSyjLmFmCrB+il4x5lPXNOhnUl7hhmYgkKqhui4ZHnTzz7vhP/r55/TdrddGvhZZ1+7JS109eF9HT71wbWCQlNIH/qiBH8W+GJdPXRO1UEgZgVQCxs0u1TWLMOjm/+HHns38mVevH1w1N+haDJpiiGG5dchM4yiBX11vcHXL2MbEp2EjTR6RMgKpyPlkl+qcRci7+Q9bZbffuNcohmtXRaaxjje4umZsY+EzzV3F1DjBNEIhkAos9JvTJ7tU1yzCIHk3/O+77VBmL6dBxdZZ1+6QmVbWu5V+0FaRaUyt9msUdc3YxsRnmrvMqXGCaYREIBVQEW/OQRmSvKCtqCxCzE9weTd8SWNdi6xrJ0nbzlX+QVtVtizEDS6m352Ys44oB8E0QiKQCqiIN2dehmT6cGto0BbyxpXCE9ygG37etci6wZ89fUIPP35Z2/u2T6r6gzbVTGNsvzupXkeEQzCNkOgjFVARb868BnrOaWCvlvmZji4untSLS/fr4uLJiW9YPr1hYpF3LfJ6JEnSazl7UFb5QVtGM8UixPa7k+p1RDij9lEDRkEgFVARb868BnrfztnDragbfR2f4B79/HO5LRMOmWX+TJUftGU0UyxCbL87qV5HhEMwjZCY2guoqNqkrCmrQc0pi1C36ZCV9W5mawQpu2WCJJmkt999pMBRDZdSr6qeGH93UryOCKeOiyhQnSCBlJm9S9InJE1J+gPn3FKI46amzDdn2cvS67YM/pHPPTf2zzhJT17qavZH/nklH7i9eq7uxqamzLTtnDoJ3ADq9ruD+I2yuIFgGqFMHEiZ2ZSk35f0ryW9LOlLZvY559xXJj12isp6c5b9RDXO+YpYoRX6mFktEUZRVcH5/oLtXiF81YXbo+DpH2WKbXED6i9ERuqtkv7GOfdVSTKzz0p6r6RGBlJl2n+D6i80L+p8w45dxIdYbB+MVdT2ZBVs92xubevhxy9LivdGwdM/ykJrA5QtRLF5R9JLfX9+efc1FCxv1dnKejfIseeWVnVs8bzmllZHPmYRK7Tyjvnw45fHHl/PGw+3vMdTRW3PsOCt1+cqxL89kLLYFjeg/kpbtWdmD5nZmpmtXbt2razT1lpRy8onCdCK+BDL+9lt57wDyI+95x61prJX5g0zam2PbzCaZZTgreiWAiH/PkBRaG2AsoUIpLqS7uz78x27r+3hnPukc27WOTd75Ei1K5/qoqgnr0kCtCI+xIoIIuZnOlr+2Xs1ldPmIO/16XZrpOmB0NnCrOXaWYp66i4y+7n/PARrmAStDVC2EIHUlyT9mJkdM7PXSXq/pM8FOC6GKOrJa5IArYgPsaKCiPmZjn775+7NHO8H3nanWof2BlOtQ6ZHHrhnpGOHzhb29z4apKin7jKaapYVrKHe6BOGsk1cbO6cu2FmvyrpgnbaH/yRc278teUYW1HLyifp+1PECq39xzy0u/TfZ3yjjleSHvvSS3u+d9s5Pfr55/Thx54d+vcKmS3cv2Lx4+97s6Tx9hCcVBl1JxQJIxQWN9wS0z6XdRWkj5Rz7s8l/XmIY2F0RS0rnzRAK+JDrP+Y+1fxjTO+vA+V/eOdW1rV1vbeYO01d6tZ57CVgyGaUK6sd/Xo55/b0yC0d96zp0/o7OkTpX1AltFUkyLhcnBjbY7YVjzXFZ3NE1dU0LL2tf+nzzzzkrad05SZHrwvnic83wBynA+VUW7eg7IlkwajWcHi/vOG2ENxVGU01YyxA3rdcGNtFrK85SCQwgEr6109eal7c/ps27mBHb2reMLNCyAHjWWcD5Xpw63crWL65QVck2YLB/WNGnTeopTRVJMO6MXjxtosZHnLQSCFA8b5sI3pCXfYWEb9UFlZ7+ofv3tjpHMOypZMki0c9kFXRZam6LoTOqAXjxtrs5DlLQeBFA4Y58M29BOuT3arfw+6/frHMuqHyvKFq9p67WAxe5br37uhlfVu8Jt93lilemdpKBIuFjfWZiHLW47SGnIiHdM5Xb+zXs8Lurobm2P3A/JZ/t7/M3l6Yxy1NcOgp/Pp9t5r8Or1rUKW6Oe1fJhut/TgfR0tX7hKryWMjR5LzUIriHKQkcIBGZ0Fcl/Pe8I16ebro073+WS3htUS9cbYf+5hGa+8v1Ovh9P+DY+LqDEZ1JYhlqlUpIfp0+Yhy1s8AqkCpbrM+Nub2UXWWa9npY5N0v6Ya5Rgw6d+Y1htR7s1pbfffURzS6sj/zvk/Z3efvcRffrpr2f+zKCMmK+8tgwUC2MS3FiBsJjaK0gKXZrztuMYp2N6Vuo4r7rIt4B6UP3GoK91ptt68L6OnrzUHevfYX6mowfv66i/r7mT9OSlbu60p0ml/NtSLAwAcSGQKkgZW2pMYlCgN24dxfxMRxcXT+rFpft1cfFk7jYmwwpafeo38n7m4+97sy4untRTL1zz+nd46oVrmVk156SsXfictOeYRe0Zx4asABAXAqmCxJ45GFaPNEmBom9Bq895h/2M779D3te/vbk1NONWZDaSYmEAiAs1UgWJfZnxsABjkjqKSQpafc476Gd8/x2G/dygr/m2hBilpo5iYQCIC4FUQWLv31F0oDdKQFRGMb7vv8Ownxv0NZ8s2DiNTSkWBoB4MLVXkNj7d1Q9RVRWMb7vv8Ognxt2TJ86pthr6gAA2czlNQ0q0OzsrFtbWyv9vNiryvYMc0urub2aLi6eLGUMg0xybbI2HG63pgYGcMcWz2fWXpmkF5fu9/gbAABCMbNLzrnZrK8xtddgVU4RxVyMP+n+gT51TLHX1AEAshFIoRIxBw4h9g8cN0iNvaYOQHFSbd6MHdRIoRJV12gN4pstm6R3VOw1dQCKkULzZgxGRgqViHkZv0+2bNh04KitDWL4+wMoT4gMOKpFIIXKVBk4DApsfKbZhq26C73RMFMBQD3EXC+K0TC1h8YZlkr3mWYb9GEYurUBUwFAfbDtU/rISCFpPpmZUVLp42bLBk0Hhn7iZCoAqA8WmtySaqadjBSS5ZuZGRbY+BSNDyqeD/3EyVQAUB8sNNmRcqadjBSS5ZuZGZQ9yioaX3jish753HP69uaW9x54IZ84Y24dAWB8LDRJO9NOILUr1ZRik/lmZgal0rPezFvbThubW5L89sALvUKRqQAAdZNypp1ASpN3skY1fDMzgwKbDz/27NDzjvKU9BsrV/SZZ17StnOaMtMH3nZnsK1vYm4dAQA+Us60E0gp7ZRik02SmcnLHk0fbunV61tDf37QU9JvrFzRp57++s0/bzt388+/NX9i6LFHwVQAgDpJOdNOIKV0UopMP+4VOjOzst7VP373xkjfO+gp6TPPvJT7eqhACgDqJOVMO4GU0kgppjj9WEbgFyozs7Le1cOPX9a2c0O/d9hTUt4xRjk2ADRVqpl22h8o7n3fekI3dSxaSktZe2MdFOiMszR5ymys1wEA6SIjpTRSiqlMP/akVHeWNdZ+nen2WIXiH3jbnXtqpPpfBwDUC4HUrthTiilMP/ZLKfAbNCafzGSvDmr/qj3qowCgfmoXSNW1IDu1FQ0pBX55Y50y8+4w/FvzJwicAKABalUjlVJdzrhS20Yghbqznryx/vbP3Rvt9QUAxKFWGamU6nJ8xD792C+FurOelMYKAIhLrQKplOpymiC1wC+VsQIA4lGrqb28+psY63IAAED6apWRSqEgu67F8LHg+gIAylSrQCr2WpcUu5OnhOsLAChbrQIpKe5al7oXw1eN6wsAKFutaqRiRzF8sbi+AICy1S4jFbNJm1RS/zNYSk1AAQD1QCBVokmK4an/GS7EYgPfYLXsIJegGgiP9xV8EEiVaJJi+CLqf+r2oTHpYgPfYLXsIJegGgiP91V6YrmHmXOu9JPOzs66tbW10s+bmv5fkrx/JZP04tL9XsfOyt7EvO1M0eaWVjOnBjvTbV1cPBn853yVfT6gCXhfpaXse5iZXXLOzWZ9jWLzSO3fNzCPb/3PoAxXU/kWq5dd5E5RPRAe76u0xHQPI5CKVNYvyX6TNBvlQ+Mg3874ZXfUp4M/EB7vq7TEdA8jkIrUoF8G0066eZIUJh8aBy2cOq52a2rPa6MEq74/56vs8wFNwPsqLTHdwyg2j1TeUv5Q8/UpbKdTNt9i9bI76sfewR9IEe+rtMR0D6PYfEyhVwnkHa+MQrpQf5dYVk4AAJqjzHvPoGJzAqkxhA5uhh0vtgAlazySgl+TmP7OAAAQSAUSenls3vF6x4wpiMgL+l7fOqRXr28d+H6fa5JSSwYCPgBoDtofBBJ6lcCgn+s1g1tZ73odO7S8paZZQZTkd01iWs46yP7WFLH9WwEAykMgNYbQqwSG/dzm1rYefvxyFDfocQMjn2sS03LWQVIJ+AA0x8p6V3NLqzq2eF5zS6tR3DeagkBqDKGXx2Ydb79t56LIduQFRtPtVrBrEtNy1kFSCfgANANZ8moRSI1hfqajs6dPqDPdDtLLqf94g8SQ7cgLIh954J5g1ySVPi6pBHwAmoEsebXoIzWm+ZlObpDgU4DcO15WoXW/rGxHmQXPeT1Wsl6bJLD0OV7Zhd8x9S8BALLk1SKQCmTSncN73/Pw45e1nbGScn+2o4qdyvcHkUWMYVCgmqWq6yDRuA9AHPIaOJMlLwdTe4GESK3Oz3T02z9370jTWzGkcps8hvmZjhZOHdfR6bZe2djU8oWr1CMAqEQqZRF1RUYqkLwUandjU3NLqyNnLkbNdgw7XxkZkhjSyVWNoYpMGABkIUteLQKpQPJSqybdfH3Um+0o01t55xvnPJOKIZ1c1RgGZcL48AJQtnHLIhAOU3uBZKVWTdL+aqdQ007DWieUMb0VQzq5qjHEkI0DAFRvokDKzJbN7AUz+2sz+1Mzmw40ruRktUbI23wnxM12lNYJRd/UQ7eDSGkMtEAAAEgT7rVnZu+UtOqcu2Fm/0WSnHP/cdjPpbrX3rhG2ZsvxNL90HsAYriU9gUEAEymsL32nHNfcM7d2P3j05LumOR4dTNs2ilUN9oYpthSM+l2CjFk4wAA1QtZbP7Lkh4LeLzkDVtJEapguegVGyvrXT3yuee0sbmzQfEbD7f0sffck2zQEGrFHcWdAIChU3tm9peSfjjjSx91zv333e/5qKRZSaddzgHN7CFJD0nSXXfddd/Xvva1ScZdC8cWz2fWUZmkF5fuL3s4mVbWu1r4k8vaem3vSFtTpuWfvTfJQIKpUADAOAZN7Q3NSDnnfnrIwX9R0rslvSMviNo9ziclfVLaqZEadt4miKF9wDDLF64eCKIkaWvbJbvUnxV3AIBQJl219y5JH5H0gHPuepghNUcKtU2DgotUA4+8QNVJXvVSAIDmmrRG6vckfZ+kL5qZJD3tnPuViUfVEKFqm0Jt2pt1nEGNP2PKnI0ja9PhHjqUAwDGMVH7A19NaX9QhlDL8POO8+B9HT32v1+qVY2UdCtozAsSqZcCAPQU1v4A1Qu1aW/ecZ564ZqW/+29mm63br7+xsOtpIMoaSfbdHHxpCzn64OmLSdtnQAAqA/22ovcsGm7QYXT40z5DTpOnZf5j1vwz2bFAIB+ZKQGqDLzsLLe1Zsf/YI+9NizAxt25t3wb2+3tPDE5T0/u/DE5dy/Q1O3PBm34D9UBhAAUA8EUjlCdR2f5Ny9Bpj99t+08wKBre3XtLW9t65pa9vp0c8/l3nOrOO0Dpmuf+9Graewxu1QTusEAEA/pvZyhOo6Hurc/fpv2nkr/z702LOZP/vq9YPBWdZxbm+39J3v3bj5/XWewhpn6jKF3l8AgPIQSOWoMvMw7Bz7b9pZgUBeIDVI/3HmllYPZMTKCiSlcC0dQstqnRBb7y8AQHkIpHJUmXkY1Ltp1Jv2dLuVOTXYv/qu3/7AJe/8ZQSSMRd0F72vIQAgLdRI5aiy63jWuaWdtgOj9od65IF71Dq0d3F/65DpkQfuOfC9WfVgeW0ByggkYy/o7rVOeHHpfl1cPEkQBQANRkYqR5WZhxDnHucYWYGL087myf3l6mUFkhR0AwBSQSA1QJX9k0Y597A6olHHnxegOO2sYis7kKSgGwCQCqb2EhWyPcPtOXVTnb799l7Z2NTyhaultEBIYTNnAAAkAqlkhaojWlnv6jvfu3Hg9dYh09vvPlJKL639jU8ljdXbCQCAqjC1F7m86btQdUTLF64eaNwpSd//+tv01AvXCu+llbdC7+zpE2waDACIHoFUxAa1AQhVR5QXeG1c39JGTvPOkEXfVTY+BQBgUkztFSTEPn2DgoxQdUTTh7Pro6YPt0rZf48VegCAlBFIeRoUKIUqBB8UZIy7R1wed3BW7+brZRR9N3WzZABAPTC152FY5+1Q01XDpu9CtGf4dkb3897rZfTSYssVAEDKCKQ8DAuUQk1XjRNk+O5NV0awNghbrgAAUkYg5WFYoBSqEHzUIGOSveliyAhV2fgUAIBJEEh5yAuUpg+3NLe0enOvuhDbq4wSZEwylUhGCAAAfwRSHrKyOK0p0z9+94Ze3W0Z0L9XXafg4GTSqUQyQgAA+CGQ8pCVxfnOP93Qxr7C7V4QVXRjybL2pvOtwwIAoK4IpDztz+IcWzyf+X1l9EMqo85pkjosAADqij5SgVTZDylUT6lBQu3tBwBAnZCRCqTq1W9F1znRgRwAgIPISAVSRlaoSnQgBwDgIDJSAdV59VvVGTcAAGJEINUQk664o98UAAAHmcvbtbZAs7Ozbm1trfTzNtX+FXfSTjbJd+qRNggAgCYxs0vOudmsr1Ej1QAhV9z1grLuxqacbrVBWFnvBhotAADpIJBqgJAr7miDAADALQRSDRByxV1e8JXVWR0AgLojkGqAhVPH1W5N7XnNd8VdXvBlEtN7AIDGIZBqgJA9rhZOHZdlvO4kpvcAAI1D+4OGCNXjan6mow899mzm1+hyDgBoGjJSGFuHLucAAEgikIKHkDVXAACkjKk9jI0u5wAA7CCQgpc67ysIAMComNoDAADwRCAFAADgiUAKAADAE4EUAACAJ4rNa2ZlvRtsNV3IYwEAUEcEUjWyst7VmXNXtLm1LWlnI+Ez565I0tgBUMhjAQBQV0zt1cjyhas3A5+eza1trz3wQh4LAIC6IpCqkby97nz2wAt5LAAA6opAqkby9rrz2QMv5LEAAKgrAqkaCbkHHvvp7dSJzS2t6tjiec0trWplvVv1kAAAkaHYvEZC7oHX9P30KLYHAIzCnHOln3R2dtatra2Vfl5gVHNLq+pm1IN1ptu6uHiyghEBAKpiZpecc7NZX2NqD8hAsT0AYBQEUkAGiu0BAKMgkAIyUGwPABgFxeZAhqYX2wMARkMgBeSYn+kQOAEABmJqDwAAwBOBFAAAgCcCKQAAAE8EUgAAAJ4IpAAAADwFCaTM7GEzc2b2gyGOBwAAkIKJAykzu1PSOyV9ffLhAAAApCNERup3JX1EUvm7HwMAAFRookDKzN4rqeucuxxoPAAAAMkY2tnczP5S0g9nfOmjkn5dO9N6Q5nZQ5IekqS77rprjCECAADEyZzzm5EzsxOS/qek67sv3SHpFUlvdc59c9DPzs7OurW1Na/zAgAAlMnMLjnnZrO+5r3XnnPuiqQf6jvJ30qadc59y/eYAAAAKaGPFAAAgCfvjNR+zrk3hToWAABACshIAQAAeCKQAgAA8EQgBQAA4IlACgAAwBOBFAAAgKdgq/aA0FbWu1q+cFWvbGzq6HRbC6eOa36mU/WwAAC4iUAKUVpZ7+rMuSva3NqWJHU3NnXm3BVJIpgCAESDqT1EafnC1ZtBVM/m1raWL1ytaEQAABxEIIUovbKxOdbrAABUgUAKUTo63R7rdQAAqkAghSgtnDqudmtqz2vt1pQWTh2vaEQAABxEsTmi1CsoZ9UeACBmBFKI1vxMh8AJABA1pvYAAAA8EUgBAAB4IpACAADwRCAFAADgiUAKAADAE4EUAACAJwIpAAAATwRSAAAAngikAAAAPBFIAQAAeCKQAgAA8EQgBQAA4Mmcc+Wf1OyapK/t/vEHJX2r9EHUF9czPK5pWFzPsLieYXE9w6rL9fwR59yRrC9UEkjtGYDZmnNuttJB1AjXMzyuaVhcz7C4nmFxPcNqwvVkag8AAMATgRQAAICnGAKpT1Y9gJrheobHNQ2L6xkW1zMsrmdYtb+elddIAQAApCqGjBQAAECSCKQAAAA8RRVImdnDZubM7AerHkvKzGzZzF4ws782sz81s+mqx5QiM3uXmV01s78xs8Wqx5MyM7vTzJ4ys6+Y2XNm9sGqx1QHZjZlZutm9mdVj6UOzGzazJ7Y/fx83sx+quoxpczMPrz7fv8/ZvYZM3t91WMqQjSBlJndKemdkr5e9Vhq4IuSfsI59y8k/V9JZyoeT3LMbErS70v6N5J+XNIHzOzHqx1V0m5Ietg59+OSflLSv+N6BvFBSc9XPYga+YSkv3DO3S3pXnFtvZlZR9K/lzTrnPsJSVOS3l/tqIoRTSAl6XclfUQS1e8Tcs59wTl3Y/ePT0u6o8rxJOqtkv7GOfdV59z3JH1W0nsrHlOynHPfcM59eff//0E7N6hOtaNKm5ndIel+SX9Q9VjqwMxul/SvJP2hJDnnvuec26h0UOm7TVLbzG6TdFjSKxWPpxBRBFJm9l5JXefc5arHUkO/LOl/VD2IBHUkvdT355fFjT8IM3uTpBlJz1Q8lNR9XDsPn69VPI66OCbpmqQ/3p0u/QMze0PVg0qVc64r6b9qZ5bpG5K+7Zz7QrWjKkZpgZSZ/eXuPOn+/94r6dcl/WZZY6mDIdez9z0f1c6UyqerGylwi5l9v6QnJX3IOff3VY8nVWb2bkl/55y7VPVYauQ2SW+R9N+cczOSviOJ2khPZvZG7WTxj0k6KukNZvbz1Y6qGLeVdSLn3E9nvW5mJ7RzoS+bmbQzDfVlM3urc+6bZY0vNXnXs8fMflHSuyW9w9EszEdX0p19f75j9zV4MrOWdoKoTzvnzlU9nsTNSXrAzH5G0usl/TMz+5RzrpY3qpK8LOll51wvU/qECKQm8dOSXnTOXZMkMzsn6V9K+lSloypA5VN7zrkrzrkfcs69yTn3Ju38Mr+FIMqfmb1LOyn/B5xz16seT6K+JOnHzOyYmb1OO0WSn6t4TMmynaekP5T0vHPud6oeT+qcc2ecc3fsfma+X9IqQdRkdu85L5nZ8d2X3iHpKxUOKXVfl/STZnZ49/3/DtW0eL+0jBRK9XuSvk/SF3ezfE87536l2iGlxTl3w8x+VdIF7aw2+SPn3HMVDytlc5J+QdIVM3t297Vfd879eXVDAg74NUmf3n14+qqkX6p4PMlyzj1jZk9I+rJ2SkzWVdPtYtgiBgAAwFPlU3sAAACpIpACAADwRCAFAADgiUAKAADAE4EUAACAJwIpAAAATwRSAAAAnv4/rWmM3i+PqEUAAAAASUVORK5CYII=\n", "text/plain": ["
"]}, "metadata": {"needs_background": "light"}, "output_type": "display_data"}], "source": ["import numpy as np\n", "y0 = np.random.multivariate_normal([0, 0], [[2, 0], [0, 0.1]], size=50)\n", "y1 = np.random.multivariate_normal([0, 0], [[0.1, 0], [0, 2]], size=50)\n", "y2 = np.random.multivariate_normal([5, 2], [[2, -1.5], [-1.5, 2]], size=50)\n", "y3 = np.random.multivariate_normal([-2, -2], [[0.5, 0], [0, 0.5]], size=50)\n", "y = np.vstack([y0, y1, y2, y3])\n", "X=y\n", "\n", "fig, ax = plt.subplots(1, 1, figsize=(10, 7))\n", "ax.plot(y[:,0], y[:,1], \"o\");"]}, {"cell_type": "markdown", "metadata": {}, "source": ["

bayespy

"]}, {"cell_type": "markdown", "metadata": {}, "source": ["The module [bayespy](http://www.bayespy.org/) allows to build and estimate simple bayesian models. I just replicate the example on the [Gaussian mixture model](http://www.bayespy.org/examples/gmm.html)."]}, {"cell_type": "markdown", "metadata": {}, "source": ["We define the model:"]}, {"cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [{"data": {"text/plain": ["('0.5.22', '1.22.1')"]}, "execution_count": 6, "metadata": {}, "output_type": "execute_result"}], "source": ["from bayespy import __version__ as v1\n", "from numpy import __version__ as v2\n", "v1, v2"]}, {"cell_type": "markdown", "metadata": {}, "source": ["Si cela ne marche pas avec la version 1.14 de numpy, vous devriez essayez la version 1.13 (a priori, cette [exception](https://github.com/numpy/numpy/blob/master/numpy/core/einsumfunc.py#L710) pose probl\u00e8me)."]}, {"cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": ["N = 200 # number of data vectors\n", "D = 2 # dimension\n", "K = 10 # maximum number of clusters"]}, {"cell_type": "code", "execution_count": 7, "metadata": {"scrolled": false}, "outputs": [], "source": ["from bayespy.nodes import Dirichlet, Categorical, Gaussian, Wishart, Mixture\n", "alpha = Dirichlet(1e-5*np.ones(K), name='alpha')\n", "Z = Categorical(alpha, plates=(N,), name='z')\n", "mu = Gaussian(np.zeros(D), 1e-5*np.identity(D), plates=(K,), name='mu')\n", "sigma = Wishart(D, 1e-5*np.identity(D), plates=(K,), name='Lambda')\n", "Y = Mixture(Z, Gaussian, mu, sigma, name='Y')"]}, {"cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": ["Z.initialize_from_random()"]}, {"cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": ["from bayespy.inference import VB\n", "Q = VB(Y, mu, sigma, Z, alpha)"]}, {"cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": ["Y.observe(y)"]}, {"cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": ["import time\n", "if not hasattr(time, 'clock'):\n", " # bayespy still clock and it was removed in python 3.8\n", " setattr(time, 'clock', time.perf_counter)"]}, {"cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [{"name": "stdout", "output_type": "stream", "text": ["Iteration 1: loglike=-1.497701e+03 (0.012 seconds)\n", "Iteration 2: loglike=-1.359319e+03 (0.009 seconds)\n", "Iteration 3: loglike=-1.345008e+03 (0.007 seconds)\n", "Iteration 4: loglike=-1.333128e+03 (0.007 seconds)\n", "Iteration 5: loglike=-1.321888e+03 (0.008 seconds)\n", "Iteration 6: loglike=-1.308002e+03 (0.007 seconds)\n", "Iteration 7: loglike=-1.290391e+03 (0.007 seconds)\n", "Iteration 8: loglike=-1.275647e+03 (0.006 seconds)\n", "Iteration 9: loglike=-1.264646e+03 (0.008 seconds)\n", "Iteration 10: loglike=-1.251051e+03 (0.006 seconds)\n", "Iteration 11: loglike=-1.229172e+03 (0.006 seconds)\n", "Iteration 12: loglike=-1.182531e+03 (0.006 seconds)\n", "Iteration 13: loglike=-1.162046e+03 (0.007 seconds)\n", "Iteration 14: loglike=-1.129390e+03 (0.008 seconds)\n", "Iteration 15: loglike=-1.114487e+03 (0.007 seconds)\n", "Iteration 16: loglike=-1.089347e+03 (0.008 seconds)\n", "Iteration 17: loglike=-1.086340e+03 (0.007 seconds)\n", "Iteration 18: loglike=-1.083643e+03 (0.008 seconds)\n", "Iteration 19: loglike=-1.080921e+03 (0.007 seconds)\n", "Iteration 20: loglike=-1.077116e+03 (0.008 seconds)\n", "Iteration 21: loglike=-1.064775e+03 (0.007 seconds)\n", "Iteration 22: loglike=-1.052069e+03 (0.006 seconds)\n", "Iteration 23: loglike=-1.031870e+03 (0.006 seconds)\n", "Iteration 24: loglike=-1.020058e+03 (0.006 seconds)\n", "Iteration 25: loglike=-1.008205e+03 (0.011 seconds)\n", "Iteration 26: loglike=-9.924005e+02 (0.011 seconds)\n", "Iteration 27: loglike=-9.920968e+02 (0.016 seconds)\n", "Iteration 28: loglike=-9.917380e+02 (0.018 seconds)\n", "Iteration 29: loglike=-9.912025e+02 (0.008 seconds)\n", "Iteration 30: loglike=-9.903257e+02 (0.009 seconds)\n", "Iteration 31: loglike=-9.892435e+02 (0.008 seconds)\n", "Iteration 32: loglike=-9.877160e+02 (0.007 seconds)\n", "Iteration 33: loglike=-9.859192e+02 (0.006 seconds)\n", "Iteration 34: loglike=-9.836574e+02 (0.009 seconds)\n", "Iteration 35: loglike=-9.814675e+02 (0.007 seconds)\n", "Iteration 36: loglike=-9.798386e+02 (0.006 seconds)\n", "Iteration 37: loglike=-9.793723e+02 (0.006 seconds)\n", "Iteration 38: loglike=-9.792857e+02 (0.006 seconds)\n", "Iteration 39: loglike=-9.792329e+02 (0.008 seconds)\n", "Iteration 40: loglike=-9.791737e+02 (0.006 seconds)\n", "Iteration 41: loglike=-9.791036e+02 (0.008 seconds)\n", "Iteration 42: loglike=-9.790199e+02 (0.007 seconds)\n", "Iteration 43: loglike=-9.789201e+02 (0.006 seconds)\n", "Iteration 44: loglike=-9.788011e+02 (0.006 seconds)\n", "Iteration 45: loglike=-9.786598e+02 (0.006 seconds)\n", "Iteration 46: loglike=-9.784926e+02 (0.007 seconds)\n", "Iteration 47: loglike=-9.782935e+02 (0.007 seconds)\n", "Iteration 48: loglike=-9.780498e+02 (0.007 seconds)\n", "Iteration 49: loglike=-9.777315e+02 (0.006 seconds)\n", "Iteration 50: loglike=-9.772782e+02 (0.005 seconds)\n", "Iteration 51: loglike=-9.766198e+02 (0.006 seconds)\n", "Iteration 52: loglike=-9.757762e+02 (0.010 seconds)\n", "Iteration 53: loglike=-9.748480e+02 (0.011 seconds)\n", "Iteration 54: loglike=-9.738679e+02 (0.011 seconds)\n", "Iteration 55: loglike=-9.728080e+02 (0.007 seconds)\n", "Iteration 56: loglike=-9.716088e+02 (0.006 seconds)\n", "Iteration 57: loglike=-9.701531e+02 (0.005 seconds)\n", "Iteration 58: loglike=-9.684097e+02 (0.006 seconds)\n", "Iteration 59: loglike=-9.668369e+02 (0.007 seconds)\n", "Iteration 60: loglike=-9.658912e+02 (0.006 seconds)\n", "Iteration 61: loglike=-9.654709e+02 (0.004 seconds)\n", "Iteration 62: loglike=-9.652997e+02 (0.000 seconds)\n", "Iteration 63: loglike=-9.652227e+02 (0.000 seconds)\n", "Iteration 64: loglike=-9.651824e+02 (0.023 seconds)\n", "Iteration 65: loglike=-9.651584e+02 (0.006 seconds)\n", "Iteration 66: loglike=-9.651426e+02 (0.008 seconds)\n", "Iteration 67: loglike=-9.651314e+02 (0.007 seconds)\n", "Iteration 68: loglike=-9.651230e+02 (0.007 seconds)\n", "Converged at iteration 68.\n"]}], "source": ["Q.update(repeat=1000)"]}, {"cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [{"data": {"image/png": "iVBORw0KGgoAAAANSUhEUgAAAlIAAAGfCAYAAACUQTyEAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAAsTAAALEwEAmpwYAABwKElEQVR4nO3deXicZb0//vf9zJLZMjPZ96ZtulFKaaG0hEUKRQVEBQT0uNbvERDUL+pBfwcX9DrqqV5yfQFRjogo5ygKCKKI4EKhHpbS0tLS0o02bbZm3zOTTGZ57t8fycQ0nSSzPDPPM5P367p6pU1mnvuTSTPzmfv+3J9bSClBRERERIlT9A6AiIiIKFsxkSIiIiJKEhMpIiIioiQxkSIiIiJKEhMpIiIioiQxkSIiIiJKkiaJlBDCK4R4UghxWAhxSAhRr8V1iYiIiIzMrNF17gPwFynl9UIIKwCHRtclIiIiMiyRakNOIYQHwF4AiyW7exIREdE8osWM1CIA3QB+KYQ4G8BuALdLKf1TbySEuBnAzQDgdDrPXbFihQZDExEREaXX7t27e6SUJbG+psWM1DoArwO4UEq5QwhxH4AhKeU3Z7rPunXr5K5du1Ial4iIiCgThBC7pZTrYn1Ni2LzVgCtUsodE/9+EsA5GlyXiIiIyNBSTqSklB0AWoQQyyc+tQnAwVSvS0RERGR0Wu3a+wKARyd27B0H8GmNrktERERkWJokUlLKvQBirh0SERER5Sp2NiciIiJKEhMpIiIioiQxkSIiIiJKEhMpIiIioiQxkSIiIiJKEhMpIiIioiQxkSIiIiJKEhMpIiIioiQxkSIiIiJKEhMpIiIioiQxkSLSw/btwJYt4x+JiChraXVoMRHFa/t2YNMmIBgErFZg61agvl7vqIiIKAmckSLKtG3bxpOoSGT847ZtekdERERJYiJFlGkbN47PRJlM4x83btQ7IiIiShKX9ogyrb5+fDlv27bxJIrLekREWYuJFJEe6uuZQFHytm9nIk5kEEykiIiyCTcrEBkKa6SIiLIJNysQGQoTKSKibMLNCkSGwqU9IqJsws0KRIbCRIqIKNtwswKRYXBpj4iIiChJTKSIiIiIksREioiIiChJTKSIiIiIksRicyKiGFRVRSgUQjgcnvyoKMrkH5PJNPl3q9UKIYTeIRORDphIEREBCIfD8Pl86O7uhs/nQzAYnEyOpJSTt5spYXK73fB4PHA6nbDZbLBarRmJm4j0xUSKiOatSCQymTz19vYCAPLy8mCz2eB0OuO+jpQSY2NjaG1tnUy6zGYzvF4viouL4XK5oCispCDKRUykiGjeUVUVfX19aGpqQjgchtVqhdfrTXp5TggBm80Gm802+blIJIKBgQF0dXXBbDajtLQUhYWFcDqdXAYkyiFMpIhoXvH5fDh+/DhGRkbgdrthNqfnadBkMsHlcgEYT6q6u7vR3t4Oi8UymVQ5HI60jE1EmcNEiojmBVVV0dbWhpaWFjidThQWFmZsbJPJhPz8fADjtVjt7e1obW2F1+tFVVXV5NeIKPswkSKinBeJRHDs2DEMDAygoKBA13ols9kMj8cDABgZGcHbb78Nt9uN6upquN1uLvsRZRkmUkSU06YnUUbicDjgcDgQCARw+PBh2Gw21NTUwOv1sjidKEswkSKinGXkJCrKtmcPCnbuxOj69Rg680wcPXoUVqsVtbW1KCgo4AwVkcExkSKinKSqalYkUdWbN0MEg5BWK1ofeQTWtWsRDAZx5MgReL1e1NbWsiidyMA4d0xEmbN9O7Bly/jHNGtrazN0EgUA9p07IYJBCFWFCIVg37kTAGC1WlFUVIRAIIB9+/ahqakJoVBI52iJKBbOSBFRZmzfDmzaBASDgNUKbN0K1NenZSi/3z+5K87IRtevh7RagVAI0mLB6Pr1p3zd6XTC4XCgq6sLXV1dWLRoEYqKirjcR2QgTKSIKDO2bRtPoiKR8Y/btqUlkVJVFQ0NDbDb7YYv2A6sXYvWRx6BfaJGKrB27Wm3EULA4/EgHA7j2LFjaGtrw6JFi9gygcggmEgRUWZs3Dg+ExWdkdq4MS3DdHZ2YnR01NBLelMF1q6NmUBNZzabUVhYiEAggLfffhvV1dWorKyEyWTKQJRENBMmUkSUGfX148t527aNJ1FpmI0aHR1Fc3PzZJ+mXGSz2ZCXl4e2tjb09fWhrq5usoM6EWUeEykiypz6+rTVRQFAR0cHzGaz4Zf0UiWEQEFBAQKBAPbv34+amhpUVlbm/PdNZET8rSOinBAKhdDd3T2vZmdsNhsKCgpw8uRJHDx4EKOjo3qHRDTvMJEiopzQ19cHKeW829GmKAoKCgoQCoWwb98+9PT0QEqpd1hE8wYTKSLKelJKtLW1zavZqOmcTify8/Nx9OhRHDt2DOFwWO+QiOYFJlJElPWGh4cRDAZhsVj0DkVXZrMZRUVFGBgYwIEDB7jUR5QBTKSIKOt1dHQgLy9P7zAMw+PxQEqJ/fv3Y2BgQO9wiHIaEykiymqqqmJgYAB2u13vUAzF4XDA4XDg0KFDaGtrY90UUZowkSKirDY6Ojovi8zjYbVaUVBQgKamJtZNEaUJEykiymo+n49J1CwURWHdFFEaMZEioqw2MDAAm82mdxiGx7opovRgIkVEWUtVVQwODrLQPE7RuqnDhw+jq6tL73CIcgITKSLKWtFlKi7txc9qtcLj8aChoYFF6EQaYCJFRFkrEAgwEUiCyWRCQUEBmpub0dzczMeQKAVMpIgoa42OjsJs5tnrAGDbswcFDz4I2549cd0+erRMR0cHGhoaEIlE0hwhUW7iMxARZS0mUuNse/agevNmiGAQ0mpF6yOPILB27Zz3E0KgoKAAvb29CIVCWLp0KR9PogRxRoqIslYgEIDJZNI7DN3Zd+6ECAYhVBUiFIJ9586E7l9QUAC/34+DBw8iGAymKUqi3MREioiyViAQmPfn6wHA6Pr1kFYrpMkEabFgdP36hK/hdrsRDofZa4ooQZzDJaKsFA6H2dF8QmDtWrQ+8gjsO3didP36uJb1YnG5XBgdHcWBAwdw5pln8tgdojgwkSKirMTjTk4VWLs26QRqqmjyxGSKKD5c2iOirMREKn3sdjssFguX+YjiwESKiLKSqqrsf5RG0WTq4MGDTKaIZqFZIiWEMAkh9gghntXqmkREs2F9VHrZ7XaYzWYmU0Sz0HJG6nYAhzS8HhHRjKSUnJHKACZTRLPTJJESQlQDeB+An2txPSKieHBGKjOYTBHNTKsZqXsBfBWAqtH1iIhmpceMVKLHsOQSJlNEsaWcSAkhrgbQJaXcPcftbhZC7BJC7Oru7k51WCKijM5IRY9hKb73XlRv3jxvkymTyYTDhw9jbGxM73CIDEGLGakLAXxACNEI4DEAlwkhfj39RlLKn0kp10kp15WUlGgwLBHNS9u3A1u2QNmxI6MzUqkew5IrHA4HpJQ4cuQIQqGQ3uEQ6S7lREpKeaeUslpKuRDARwC8KKX8eMqRERFNt307sGkT8M1vIv/aa+Havz9jQ2txDEuucLlcCAaDOHbsGCKRiN7hEOmKfaSIKHts2wYEg0AkAgSDcL/5ZsaGjh7D0nP77Wh95BFNuohnM7fbjeHhYRw/fhyqyvJYmr80PSJGSrkNwDYtr0lENGnjRsBqHU+mrFYMn3tuRs+50uoYllzh9XrR29sLi8WC2tpa7qKkeYln7RFR9qivB7ZuBbZtw9j558PndMKrd0zzXEFBAdrb22GxWFBVVaV3OEQZx0SKiLJLfT1QXw8xNga89Zbe0cx7QggUFBSgubkZFosFpaWleodElFGskSKirGQymXKqs3k296hSFAVerxcNDQ3o7+/XOxyijOKMFBFlJZPJpHcImon2qBLBIKTVmpXF7CaTCW63G0eOHMGqVavgcrn0DokoIzgjRURZSQgBs9mcEzvGcqVHlcVigcvlYsNOmleYSBFR1rJYLDnRxyiXelRZrVYoioKjR4/mxM+GaC5c2iOirGW1WnOiu3a0R5V9506Mrl+fdct607lcLgwMDODEiROoq6tjWwTKaUykiChr2e12jIyMIC8vT+9QUpZrPaq8Xi96enrgcDhQWVmpdzhEacOlPSLKWg6HIydmpHKV1+tFU1MTd/JRTmMiRURZy2q1ctnIwBRFgdvtxtGjRzEyMqJ3OERpwUSKiLKW1WrVOwSag8Vigc1mw5EjRxAMBvUOh0hzTKSIKGtZLBa9Q6A42Gw2SCnR0NCQE+0qiKZiIkVEWctsNudMC4Rc53K5MDw8jObmZr1DIdIUEykiympOp5NLRlnC4/Ggvb0dvb29eodCpBkmUkSU1VwuF3fuZQkhBDweD44dO4bR0VG9wyHSBBMpojQLBoPo7+/HwMAAl6DSwOFwsO4mi5jNZuTl5bHzOeUMNuQkShOfz4c33ngDr7zyyuQLRn5+Pi6//HKsXLmShdIasVqtkFLqHQYlwOFwYGBgAE1NTVi8eLHe4RClhIkUURocO3YMjz32GCKRCMrKyiaTJr/fjyeffBKlpaX4xCc+Aa/Xq2+gOcBmswEApJTsKZVFPB4POjs74fF4UFRUpHc4REnj0h6Rxo4ePYr/+Z//gdvtRnV19SkzT06nEwsXLoTP58Mvf/lLdnzWQLTpIwvOs0u0XqqhoYH1UpTVmEgRaaixsRG/+tWvUFJSAqfTOePtSktLMTY2hl/84hcYHh7OYIS5yev1IhAI6B0GJchsNsNqteKdd95hvRRlLSZSRBoZHR3F7373OxQWFsLhcMx5+5KSEvj9fjz//POs8UmR0+nkY5ilHA4HxsbG0NTUpHcoRElhIkWkkRdeeAEjIyPIz8+P+z6VlZV46623cODAgTRGlvvsdjsTqSwWrZdifynKRkykiDRw/PhxvP7666isrIz59dKGBqx+7jmUNjSc8nkhBMrLy/HHP/6RS3wpsFgsyMvLQzgc1jsUSgLrpSibMZEiSpGqqnj22WdRVFQERTn9V6q0oQFX3n03zn36aVx5992nJVMOhwPhcBivvPJKpkLOSR6Ph3VSWSx63A/P46Nsw0SKKEVHjx5Fd3c3PB5PzK+XHzkCJRyGIiWUSATlR46cfpvycuzYsQODg4PpDjdneTwedjjPck6nEz6fDx0dHXqHQhQ3JlJEKVBVFS+88AIKCgpmvE3H8uVQzWaoigLVZELH8uWn3cZsNkMIgR07dqQz3Jxmt9vTcl3bnj0oePBB2Pbs0fS2uUiL79/r9aK5uRl+v1/DyIjShw05iVJw/PhxtLe3Y9GiRTPepquuDs/fcQfKjxxBx/Ll6Kqri3m78vJyvPbaa6ivr0+oYJ3G2Ww2mEwmqKoac4k1qWvu2YPqzZshgkFIqxWtjzyCwNq1Kd82F2n1/SuKArvdjmPHjmHVqlUwmUxpiJZIO5yRIkrBSy+9NOtsVFRXXR32XXXVjEkUMD4rBQC7d+/WLL75RAiBwsJCjIyMaHZN+86dEMEghKpChEKw79ypyW1z0Wzff6IzVXa7HYFAAG1tbekKl0gzTKSIktTd3Y2WlhZNj3kpKSnB66+/zuaESSosLNS0Tmp0/XpIqxXSZIK0WDC6fr0mt812sRKjmb7/6ExV8b33onrz5riTKa/Xi9bWVu5mJcPj0h5Rkvbv36/5soPNZkNHRwcaGxtRN8vsFcU2Wzf5ZATWrkXrI4/AvnMnRtevn3WpKpHbzsS2Z09K90/39aLXjLWEN9P3P3WmChMzVfHEIoSAy+XCsWPHcNZZZ03O2BIZDf9nEiUhEolg586dKCkp0fzaTqcTu3fvZiKVBIvFAqfTiWAwCKvVqsk1o0mC1redTusaq3Rcz75zJ8xtbTMmRrG+/+hMFUKhhGfq8vLyMDo6ipaWllnrEIn0xESKKAmNjY3w+/1pSaSKiopw4MAB+Hw+uFwuza+f60pKStDc3KxZIpUpyc7cZOJ6pyRlZjOk2QxEInElRqnO1Hk8HrS3t6OwsHDGFiNEemIiRZSE3bt3a76MFKUoCqSUOHz4MNatW5eWMXKZy+XKyoaO02duIl4vCh58MOlluVRmgqY7JSmLRDB4ww0IVVbGHVsqM3VCCLjdbhw7dgyrV6+GxWJJ6jpE6cJEiihBwWAQhw4dQkVFRdrG8Hq92Lt3LxOpJDgcDphMJkQikazaOj915ibi9aL0P/8zpWU5LWq2oqYnZUPXXJPRGjCr1YrR0VG0trZyiY8Mh4kUUYJOnjwJVVXT+iKdn5+PlpYWjIyMwOFwpG2cXBRtgzAwMJB1S6PRmZuCBx/UZFkulZmg6dfRopA+lZott9uNjo4OFBUVwe12Jzw+Ubqw/QFRgt55552YO4hmOpg4GdHlvebm5pSvNR9p3QYh04zYSiGwdi36b7kl6cQs1T5b0V18x48fZ3sQMhTOSBElQEqJ/fv3o7Cw8JTPRw8mVsJhqGYznr/jjlmbb8bDbrfj4MGDWLFiRUrXmY+i9WtSSgghdI4mcVouy02VjnYI8dKiZisvLw/9/f1ob29HdXV1GqIkShwTKaIE9PT0YHh4+LRu5lMPJsbEwcSpJlIFBQU4dOhQ1tX6GIHFYoHX60UgEEjbGXzpptWyXJTeR9holRx6PB60traisLCQy95kCFzaI0pAY2NjzM/HczBxoiwWC4LBINrb21O+1nxUWlqKQCCgdxiGYYQjbFJdHgTGl71tNhtOnDgBKaWG0RElhzNSRAnYv39/zF428R5MnChFUdDQ0MBljCTk5+dDCJG1y3ta07Idgt4cDgf6+vrQ29uL4uJivcOheY6JFFGcRkZG0NzcjKqqqphf76qr0yyBiiooKMBbb72FSy65RNPrzgdmsxkFBQXw+Xxp6/mVTdJVd6UXt9uNEydOwO12Z13zVcotXNojilNrayuklFCUzP3aOJ1O9PT0YGBgIGNj5pLS0lIEg0G9w4gp1sG/6abF0ppRmM1mCCHQ2tqqdyg0z3FGiihOR44cQV5eni5jNzc3w+v16jJ2NnO5XBBCQFXVjCbAc9G78DtX5Ofno7OzEyUlJcjPz9c7HJqnjPPMQmRwR48eTemsr2T7TDkcDjRo0JtqPjKZTCguLsbIyIjeoZzCCIXfuUAIAafTiePHj2flsUCUG5hIEcXB7/djcHAQNpstqftH+0yd+/TTuPLuuxNKpvLz85lIpaC4uNhwy3tGbLiZrWw2GwKBAHp6evQOheYpLu0RxaGrqyup+5U2NKD8yBE4e3uT7jNls9nQ1dUFn8+XdUeeGIHL5TLc2Xu5Vvitt/z8fDQ3N6OwsDDmqQNE6cT/cURxaG9vT3gL/dRu59JkgqoogJRJ95nq7u5mIpUERVFQWlqKnp4eQ9XRaN1wcz4zm81QVRXt7e2oqanROxyaZ5hIEcWhoaEh4Rfhqd3OVVXFOxdfDF9RUVJ9phRFQVtbG0++T1JRURE6Ojr0DoPSyO12o62tDSUlJUkvwRMlgzVSRHNQVRVNTU0zzgbNVEQ+vdv5sQsuwL6rrkqq15TL5cLRo0eTip/G20jYbDbD1UoZnR4tGpKlKAosFgsP+qaM44wU0Rz6+/sRCoVi1l7Mdlixlt3OXS4XWlpaDLeNP1sIIVBVVYXjx4+zeWOcsrFFg8vlQl9fH4aHhw21jEu5jc/IRHPo7u6e8WtTl++UiSLyqbrq6pKehZrKbDYjHA6jt7c3pevMZ9E+XDyfLT7Z2qLB4XDwHD7KKCZSRHNoa2ubcSdQOg4rng0TqeSZzWaUlJTA7/frHcqsjLKclq0tGmw2G0ZHR9HX16d3KDRPcGmPaA6NjY0zntWWrsOKYzGbzWhra8OKFSvSNkauKy0tRWdnp95hzMhIy2nZ3KLB5XKhsbERHo+H7RAo7fg/jGgWqqpO7gSaSToOK47F6XSiqakp7ePkMqfTCbvdjmAwaMhaqanLaZhYTksmgbHt2aNJApStLRosFgv8fj86OjpQXV2tdziU45hIEc1iaGhoxkLzTHM6nTh58iQLzlNk5KLz6HIaQqGkl9MSmdXSKuEyIrfbjZMnT6KkpES3MzJpftD/1YHIwIxUk2Q2mxEKhTA4OIiCggK9w8la0aJzIyakWiynxTurZaRlxHRQFAUmkwkdHR2ora3VOxzKYcZ6FiHKtO3bgS1bxj/G0NnZabgXWyMld9nIbDajtLTUsEXngbVr0X/LLUknNfEWiWfrrrxE5Ofno6OjA2NjY3qHQjmMM1I0f23fDmzaBASDgNUKbN0K1NefcpOmpiY4HA6dAjydoijo7OzEkiVL9A4lq5WUlORsp/N4Z7W0WEY0OiEETCYTTwWgtGIiRfPXtm3jSVQkMv5x27bTEqnm5mZDNfZzOBxobGzEhRdeqHcoWc3pdMLtdiMQCOTkcSLxFIln8668ROTn56OzsxMVFRU5+bMm/RlrzYIokzZuHJ+JMpnGP27ceMqXR0dHMTIyYqhCVafTifb2dr3DyAlVVVUYGRnROwxdpbqMmA2EELBYLGhra9M7FMpRTKRo/qqvH1/O+853Yi7rDQ8PQwihU3CxWa1WDA8PIxQK6R1K1nO73Tx/b55wuVzo7u7G6Oio3qFQDko5kRJC1AghXhJCHBRCHBBC3K5FYEQZUV8P3HnnaUkUMJ5IGe2YiWhiNzw8rHMk2U8IgZqaGsMWnZN2hBCTDW2JtKbFjFQYwL9JKVcCOB/A54QQKzW4LpGuBgcHDZdIAeMvCkyktOH1eifPMaTcFp2Vmu/LuaS9lBMpKWW7lPLNib8PAzgEoCrV6xLpraury1D1UVGqqmJoaEjvMHKCyWRCVVUVfD6f3qFQmgkhYLVacfLkSb1DoRyjaY2UEGIhgLUAdmh5XSI9dHd3G3KXj8ViQXd3t95h5IyioiIA4wlqooxywDDFx+Vyobe3l8u5pCnN2h8IIVwAngLwRSnlaW+XhRA3A7gZABYsWKDVsERpY9REymazoaurS+8wcobFYkFFRQU6Ozvhdrvjvl+udwbPVVarFW1tbVi6dKneoVCO0GRGSghhwXgS9aiU8vexbiOl/JmUcp2Uct1sB8ASGUEkEsHg4KAhl/aYSGmvpKQE4XA4oZq4+dAZPBc5nU709fUhEAjoHQrlCC127QkADwM4JKX8f6mHRKS/aOsDo7U/AMYTqb6+vqSWoig2m82GkpKShJZ84j2KhYwl2u28s7NT71AoR2gxI3UhgE8AuEwIsXfiz1UaXJdIN0Yu5o6e/ccCaW2Vl5cndCZbtDN4z+23c1kvy7hcLnR2drKHGGki5RopKeUrAIz3tp0oBcPDw4af8RkeHk6opodm53Q6UVRUBJ/PB5fLFdd94jmKhYxHURRIKdHT04PKykq9w6Esx87mRDH09fXBbDbuUZRSSkPPmmWrqqoqzlLME263G21tbYhEInqHQlmOiRRRDF1dXYbcsRclhEB/f7/eYeScqbNSFJ9sbQFhMpkQiUT4e0QpM+5bbiIdGbUZZ5TNZmOxbJpUVVVh3759kFIacrOBkWR7Cwin04nW1lYUFhZO1h4SJYr/c4hi6O3tNfSMlN1uZwuENHE4HCgpKTH8rJQRZoKmt4DI/8MfdI8pEVarFYFAgMvklBLOSBFNMzIyglAoZOgaKZvNhp6eHr3DyFmVlZXo6ekx7KyUUWaCoi0gEApBmkzw/P73EOFwVs1OORwOtLa2wuPxGPJnTcbHGSmiaaI9pIzMbDYjGAxidHRU71Bykt1uR2lpqWEPhzZKM9CpLSCGrrsOIhzWPaZE2Ww2+Hw+w89AknExkSKaxqgvntMJIbgkkUYVFRWIRCIJdTvPFCM1Aw2sXYv+W27B0DXXGCamRFmtVnR0dOgdBmUp465dEOlkcHBQ0+uVNjSg/MgRdCxfjq66Os2uK6XE8PAwysrKNLsm/ZPNZkNpaSl6enrg8Xj0DucU0Zkg+86dkwlLwYMPYnT9et2W06bHlA3LelFOpxP9/f0IBoOwWq16h0NZhokU0TTd3d2wWCyaXKu0oQFX3n03lHAYqtmM5++4Q7NkKppIUfpUVlaiq6vLkLVS0WagRqmXmhpTOtn27NE8WRNCQEqJ/v5+vjGhhHFpj2iaoaEhzRKp8iNHoITDUKSEEomg/MiRuO9b2tCA1c89h9KGhphfN5vNXNpLs7y8PFRUVBj6cTZKvVQmRJPG4nvvRfXmzZruDnS5XGhrazPkUi4ZGxMpommGh4c1m97vWL4cqtkMVVGgmkzoWL48rvtFZ7LOffppXHn33TGTKYvFYugX+FxRUVEBKSXC4bDeocRkpHqpdEtn0mixWBAMBjnLSwnj0h7RNMPDw3HNSMVT+9RVV4fn77gj4RqpqTNZmJjJmn5fq9XKRCoDLBYLamtr0djYiIKCAr3DOU021yYl6pR2C2lIGq1WKzo7O3mGJSWEiRTRND6fDyUlJbPeJpHap666uoTroqIzWYhEZpzJMpvNfPecISUlJWhvbzdsMfJ8OTw53Umjw+FAX1+fYX/OZExMpIimCIfDcTXjjGfGKBXxzGRZrVb2vskQRVGwcOFCHD58GIWFhXqHk5aC62yRzqQxuqGgr68P5eXlaRmDcg8TKaIpAoFAXLuz4pkxStVsM1mlDQ0oO3wY+woKDLmjLBd5PB54PB74/X44nU7d4jDSLr1MyWTi6HK50N7ejrKyMv5eUVyYSBFNEW+n8GRrn7QwdVlxrcmE4Cc+gbyNGzM2/nwlhMCCBQuwf/9+OBwO3V5kpxZcY6LgOpcTqUwnjhaLBcPDwxgeHmatFMWFu/aIpggEAnHftquuDvuuuiqjSRRweksF9cUXMzr+fOZ0OlFWVqZrbdp82qUH6NPeIS8vj53OKW6ckSKaIpFESi+nLCsqCvznnQe73kHNI1VVVeju7kYkEoHJZMr4+PNplx6Q/p16sTgcDvT39yMUCmnWU45yFxMpoikCgYDhG/JNXVZ8y+vFplWrUKx3UPOI1WpFTU0NWlpa4PV647qP1jU+82WXHqBP4hhdth0aGkJRUVHax6PsxkSKaAqfz2fYAtPpfau66urQ0twcd10Xaae0tDTudgjzsThca3okjna7HV1dXUykaE6skSKaYnBw0JBT+bN1OmcilXkmkwkLFy6Mq1ZqPh3hkktsNhuGhoYQDAb1DoUMjokU0RTxdjWfy1zn5CVqpjP7zGYzBgcHNRmDElNQUIDCwsI5e3nNt+LwXMPfL5oLl/aIptAikUqk63m8ZupbxfP29COEQG1tLd566y2oqgpFif2+dL4Vh+cSh8OBzs7OOU86oPmNiRTRFFokUunoej5T3yqr1cp3zDqy2WyoqalBa2vrrIXn86k4PJfk5eWhr68PgUAANptN73DIoJhIEU0xPDyc8rvPdHU9j9Xp3GKx8JgYnZWVlaGrqwtjY2PIy8vTOxzSmBACg4ODTKRoRkykiCZEIpG4ztmbSya7nlssFvT396ft+jQ3k8mExYsX48CBA7BarYbd9UnJcTqd6OjoQFlZmd6hkEExkSKaMDo6qtmL4Gzn5GnJbDYjEAjMWqND6ed2u1FWVoa+vj4eK5JjrFYr+vv7MTo6CrudrW/pdHzmJZqQDV3Np4smftkYe66pqamBlBKhUEjvUEhjiqKwFpFmxESKaIJeyUiqrRKEEOwlZQAWiwWLFy/W9Rw+Sg+73Y7u7m69wyCD4tIe0QQ9EimtWiVwRsoYCgsLUVBQAJ/PB5fLpXc4KdH6WJtsFl3ei6eTPc0/nJEimhAOhzN+zt5MjTYTIYRAOBxOQ3SUqGhvqXA4jEgkonc4SYsea1N8772o3rwZtj179A5Jd0IIjIyMpH2cSCSCSy65BL29vWkfi7TBRIpogh7JSLRVgqooSbdKkFIykTIQm82G2trarK6pma/H2tj27EHBgw/GTBwtFgv6+vrSHoOiKPjf//1ffO5zn0v7WKQNLu0RTdAjGdGiVQITKeMpLS1FX19f1i7xRY+1QSg0b461metwabvdjr6+PixatCitLS6EEDjvvPPw+OOP45e//CV3CmYBzkgRTQgGg7r0AOqqq8O+q65KqV0CEyljEUJg0aJFCIfDWfmziR5r03P77aclFLlqrlk4RVEQDoczsrx31113AQB++9vfpn0sSh0TKaIJY2NjWdmLiVvujclms2Hx4sVZexZiYO1a9N9yy7xIooD4DpdWFCUjP893v/vdAICvfOUrGa/bpMRl36sGUZpkayIlhEAwGNQ7DIqhqKgIhYWFbImQBeKZhXM4HOjp6Ul7LHl5ebj22mvR19eH1157Le3jUWqy71WDKE2yNZEymUwYGxvTOwyKQbz+OhY/8QTse/dy1jALzDULZ7VaMTIykpE3Lps3bwYAfP/730/7WJQaFpsTTRgbG4PJZNI7jIQpisJEyoi2bwc2bYI5GMQZVivevvdeYONGvaMiDfj9/rT3k3rPe94DIQSeffZZtLW1obKyMq3jUfKy7+03UZqEQiHOSJF2tm0DgkEgEoEIBlH5zjtZ3RKBxlkslozUSdlsNnzgAx8AAPzkJz9J+3iUvOx71SBKk2ydkWIiZVAbNwJWK2AyAVYrPB/8IEwmE+vZDGi2/lHT5eXlYWBgIP1BYXx5z2w244EHHsjK3Z/zBRMpognBYDBtM1Kpnqc3G0VReESMEdXXA1u3At/5DrB1K8wXX4y6ujoMDw9zJ5aBJNrF3WKxIBAIZCSxueKKK2C1WnGWz4d3Pv3p8eViMhwmUkQTgsFgWmakoufpnfv007jy7rs1T6YURWEhs1HV1wN33jn+EYDH40FlZWXGZjRobsl0cc/UQeE2mw13XHgh/hIOY9mvfw1s2sRkyoCYSFFu274d2LIlricfrWakpJRQVRWqqiISiaD04MGUz9Objclk4oxUFqmurobT6YTf79c7FEJ8/aNiyURjTgD4ZG0trBjfGSaDwfHaOzIU7tqj3DWxawrB4HitytatkzMDwHjCEwwGMTIygpGREZw8eRJWqxWqqmJsbAyjo6Pw+/0IBAIIhUKTXaojkcgpf8LhMFRVhZRy8g+AyS7pZ/l8OBvjv2xhKfHTQ4dwqL0dZrMZJpMJZrMZZrMZFovltL9P/WixWJCXlwebzQabzTb5d0VRWHeTRUwmE5YsWYJ9+/YhFArBYrHoHdK8Fu0fZd+5E6Pr18fVgNRms2FgYABlZWVpj2/hpz6FsZ//HBKAFAJW7vw0HCZSlHNCodD44aJPPIHSsTEIVYU6Nobt3/setq5fj8HBQQwNDcHn8yESiUBRFAghcODAAeTl5QEYT7IURYHJZIKiKJO3iX4EMJkE2Wy2yc9N/wgAbUVFuNPjweq+PrxVUIAWjweOKYlXOBxGKBSC3++fnM2ampRF/x019drRz+fl5aG7uxtFRUUoLi5GSUkJ8vPzT/vDF21jsNlsWLJkCY4cOYLCwkJdjiaifwqsXZtQB/e8vDwMDg5CSpn2n53poovw0+uvR89TT2G33Y4/n3ceX7gNhj8PyjqqqmJoaAi9vb3o6+tDd3c3WlpacPLkSbS3t2NgYAAmkwnLenvxdSFgFgJhIfCiqqKlpQVmsxl2ux1ut3tyKa+2rQ3vGhjAkYoKHC4oSFvs0WRspiXEMwYGsLqvD/sKC3HI653zetFEzO/3o729Hc3NzQgGgwiHw6ckfdGkrLi4GLW1tViyZAkqKytRVlaG0tJSOJ1OLb9NikNhYSEqKyvR2dkJbxw/65nY9uxJaDYlE4wYk5aiv1eBQCAjhwpf+rWv4eLnn4cC4Lnnnptsi0DGwESKDMvv96O1tRXt7e04efIkWlpa0N7ejq6uLqiqCkVRJhMEq9UKm80Gu90Oj8czXgxaU4OHSkuxuKUFx2tqMFRZCW+McWrb2nDzE09ACYcRPnkSd553XlxJTLzOGBjAljfegEVVEVKUGa9/xsAAvr9zJyxSIiQE/n39+jnjEEJMLhHO9WIspcTo6CgOHz6M3bt3AxgvVI9EIsjPz8eCBQuwePFi1NbWorS0FGVlZfB6vZwtSaPq6moMDw/D7/cnlcxGd5yJYBDSajXEAcNGjCkdpJQZS6TWrFmDwsJCtLS04J577mEiZTBMpMgQoklTS0sLDh48iMOHD6O3txdCCKiqOrmEZrPZUFlZGXdReFNlJZrm6Ai8uKUFpkgEJgBQVazu69M0kVrd1weLqsIEQKoqPnbsGB5dsuS0MTadPAmrlBAArFJi08mTccUhhJhMKGdLeoQQcDgccDgcKCoqmvx89NDj1tZWHDlyBOFwGCaTCaqqYkV/P+qDQSiXXori978fixYtYnKloVTrpabuOMPEjjO9kxYjxpQOFosFg4ODKEjjDHaUEAI33XQTvve972H79u04efIkqqqq0j4uxYeJFGVcNGlqbm7GoUOHcPjwYfT09EBRFKiqCrvdDpfLhZqamoy8YB+vqUHEZIIMh6EqCkpGR3HGwIBmydS+wkKEFAVyIpla29uLVf39ms58RWupkmnfIISA1WqF1Wo9ZVartq0NN7/0EkyRCMKvvILvvvwy3ikqQkFBAVavXo2zzjoLdXV1qKioYGKVglTqpaI7zhAKJbTjLJ2MGFM6WK1W+Hy+jI33yU9+Ev/5n/8JKSUefvhh3HXXXRkbm2Yn9GgMt27dOrlr166Mj0uZJ6VEd3c3Dh8+jL17956SNEUiETgcDrhcLtjtdl1fjCsaG7Hw5Zfxvq4umKWcdQkuGWcMDOBjx45hbW8vTADCAH61dCkeX7z4lNv84I03YFZVhBUF/18C4w8NDWH9+vWaFpNfumMH3vvKKzBJiYgQ+OtFF+HF9esRCAQwODiIsbExCCHgdrtx3nnnYe3atVi6dGlK9T7zWVNTU1L1UumqR4r3urFul+s1UsD4c9vg4CDOO++8jD13nXPOOdizZw9KS0vR3t6elUdaZSshxG4p5bpYX+OMFGnO7/fj6NGj2Lt3L3bs2IHe3l4A4++8o7U4RpvBOF5WBnteHsxSTi7BabnEd8jrxaNLlmBVfz/kRKK0r7DwtNv8f+edl1Cx+VSqqmoSa1R0pg6RCCImE45PzBDa7fZT6kLKT5yA54EH8MeiIrxTVITKykpceumlOP/88zOyPTxX1NTUYHh4GD6fDy6XK+77JbrjLB7x1jnNdLt0xGQ00eewYDA4uds33W677TZ88YtfxMjICF566SVs2rQpI+PS7JhIUcqklOjs7MT+/fvxv//7v3jnnXcAjNd/FBQUGDJxmk5VVexxuyeX4KYnOonupovlkNeLO+dIlA55vUldP1pLpqWmykr87MYbJ4v1Y9Wa1ba14eY//hGmSATXmkx48IYbcCAQwG9+8xs8+uijWLZsGd773vdi7dq1CSUH85GiKFi6dCn279+PYDAIq9WqWyzx1jnNl3qomUR70WUqkbr++uvx+c9/HsFgEPfddx8TKYNgIkVJUVUVjY2N2Lt3L/7xj3+gvb0dAOB2u1FdXQ1FUVDb1obF+/fP+CJsJFJKvJ2fHzPRibXrDkDcidX0JEzLQvap8adjmX6uYv3JQn0pgUgEda2taN6wAfn5+ZBSor29Hffffz9MJhPOP/98XHrppVi5ciXMZj71xJKXl4fly5fjwIED8Hg8uh2iHW+d03yph5qJEAJjY2PIz8/PyHherxeXXXYZnn/+efztb39DT08PiouLMzI2zYzPZpSQ4eFhbN++HX/605/Q1dUFADFnnaItBUwTy0I/u/FGwydTQOwZoem77jadPInL29rmbGcAxN/6QAt61DvGWv6LEkKgsLAQhYWFCIfD2LVrF1599VUUFBTghhtuwEUXXQSbzZbxmI0uPz8fixcvRkNDg27NOuPt9p1MV/BcYjabM3ZUTNRnP/tZvPLKK4hEIvjVr36FL33pSxkdn07HRIrmJKVEY2MjXnjhBbT+7ndY0dGB6sWLYVu6dMb7TJ+pWNzSYuhEarYkZOquu/BEcefUxGq2WqrpSZjWrRXiUdvWNuvyXCriWf4Dxl9wysvLAQA+nw8/+9nP8Jvf/AbXX389Nm7cCIfDoWlcs9q+ffy8so0bTzkyyEhKS0sxOjqKjo6OjGyvjyXeOqf5UA81E4vFkvFE6oorroCUEiMjI7jvvvvwxS9+0fClE7mOiRTNaGxsDLt27cIzzzyDEydO4IyBAXx9Yjt85NChWWeZZpupMKqZnoym1zYBwOVtbTMWjU81PQmb7bbpkImZwXh6dU3lcrngcrkwMjKCRx55BE899RRuuukmbNiwIf0vCHOcv2gkNTU1CAQCGB4eztjSkdEYffef2WzG6OhoRse0Wq248cYb8d///d/o6enBjh07cP7552c0BjoVE6n5IoF34aFQCC+//DIeffRRDA8Pw+12Y8GCBdjQ0RH3LFO8MxVGMtus1PQlv7mKxqfe787zzsOmkyc1jDR+Rp4ZdDgcqK2thc/nw913341zzjkH//qv/5renX7bto0nUZHI+Mdt2wybSCmKgsWLF+PAgQMYHR3NSAdtI8mGDukmkwnBYBCRSCSj9Ww33XQTnnjiCYyMjOD+++9nIqUzJlLzQZzvwlVVxc6dO/GrX/0KXV1dKCkpQeGUGZSZZplmWjpKdKYilnQuS02VaH1RokXj0Zqqy9va0lYnFet7MPLM4NSfbePChTh48CD+7d/+DV/60pdw7rnnpmfQjRvHfweivwsbN6ZnHI1YLBYsX74c+/fvh8lk0nUnX6Zly45AIcTkaQCZsmHDBrhcLvh8Pvz+97+f17OWRsBEaj6I4114W1sbHnjgARw+fBiFhYVYuHDhaZeJNcuUzqWjTBesa7GsFKtNgp51UkadGYz1sxWVlfD7/diyZQs+/elP46qrrtJ+qa++fvyNhMFrpKay2+1YtmwZDh06hIKCgnnThDGbdgRGIpGMjieEwGc+8xn88Ic/hMlkwmOPPYabbropozHQPzGRmg9meReuqipeeOEF/PKXv4TFYsHChQtnffGaPsuUzqUjIy9LxTLTDr1M1UnNNKumxcyg1mb62TqdTlRVVeHhhx/GyMgIbrjhBu0Hr6/PigRqKq/Xi4ULF6KxsVG3nXyZlk07AjOdSAHA5s2bcffdd8Pv9+Oee+5hIqUjJlLzwQzvwkdGRnDvvffizTffREVFRVJb0dO5dJTpZalU2wdsOnkSVlWFgtNnnl6YSGS2VlWlZTYq215YZ/vZWq1WLFiwAI8//jiWLl2KNWvW6BeogZSXl2N0dBTd3d267eTLtGzYESil1CWRqqurw5IlS/D222+jqakJ+/btw+rVqzMeB2mUSAkhrgBwHwATgJ9LKb+vxXVJQ9PehQ8ODuL73/8+Tpw4Mecs1GzSuXSUyWWpVJOoMwYG8J7WVggAEoA6MfM0dZYqPOUx1jqZ0qOHVCrm+tlaLBYUFxfjnnvuwX333cfz+zCeLC9cuBDhcBgDAwN8TAxCr0QKAG699VZ85StfwdjYGB544AH89Kc/1SWO+S7lxXYhhAnATwBcCWAlgH8RQqxM9bqUPoODg7jrrrvQ1NSEmonz01LRVFmJlzZsSEuik85ra2nTyZOwAIg+kjuLinDI6z2lPsoqJa5qbcWWN97AGQMDmo4vhIiZTNW2teHSHTuw4a23cOmOHahta9N03FTM9bPNz89HIBDAq6++muHIjEtRFNTV1cHpdMLn8+kdTspse/ag4MEHYduzR+9QkhbduaeHD3/4w4hEIohEIvj1r3+NQCCgSxzznRZVi+sBHJNSHpdSBgE8BuCDGlyX0kBVVTz44IPo6upCVVWV3uGkLJoopJogLO7sxCdOntQswSkaG8MZAwOT9VHR96sKAPPEsp+WZkqibn7iCbz35Zfxob//HVe8/DJufuIJQyVTcyksLMSzzz6r+TmC2cxkMmHZsmUwmUzw+/16h5O0aHuD4nvvRfXmzVmbTJlMJoTDYV3GLioqwsUXXwxgPMl+5plndIljvtMikaoC0DLl360TnyMDevHFF7Fz507NkyitEppEx7z5iSfw3ldeSSlBqG1rw5f//Gfc1Nyc9GzR1qoqhBQF0Zf7ZUND2PLGGwCAn65YgWiaM3XZT0uxZhUnC7on/q0AMIXDePdrr2VNMuVyudDX14fBwUG9QzEUi8WCFStWQFXVrJ2FmNreQEy0N8hGiqLoNiMFAJ///OeRn5+P4eFh/OhHP9ItjvksY/tohRA3CyF2CSF2dXd3Z2pYmiIcDuPXv/41ysvLNS1O1iqhiV4r3oRs6s4v08TOr2QsbmmBSVVhxj9ni84YGMCHjx+PO6k65PXi/zvvPLxZVAQV48WC0Wt5QiEIjC/7qQD+VlmZkfYH0YLu6GyYivFf+KWNjVk1MyWE0PWFyqjy8vJwxhlnIBAIZOXjE21vIE0mw7c3mI3JZEIoFNJt/CuvvHLy77t27cJJnZr/zmdaJFInAUzdTlU98blTSCl/JqVcJ6VcV1JSosGwlKjDhw9jZGRE8w7JWiU0iSZkk4mCECnt6jteU4OIoiAMIKwoGLRYsOWNN/DJo0cTmqE65PXi0SVLEJpyrX2FhZPLe2EAIUXB1jhnAxNJ5mIt7UULuv968cV46t3vxtGFCyExMTOVws9pJumYlVRVFVJKOJ1Oza6ZSxwOB8444wz4fD5dX8yTEW1v0HP77YbsWh4vIYSuS89WqxUf+9jHYDKZIITAL3/5S91ima+02LX3BoClQohFGE+gPgLgoxpclzS2Z8+etHTf1apNQaJ9o7Ta1ddUWYl7rr4anj17cKSiIqUGmtPP5YveL94jZaJm6kk1m1izjFN7SHWUlGBRa+tkE0wt20mko3lqbVsbyo8cwdDatXC5XBpFmnvy8/OxbNkyHDlyBF6vN6MdtlOVDe0NssEtt9yC//mf/8HIyAgeeOABfP3rX8+6lijZLOVESkoZFkJ8HsBfMb6i8Qsp5YGUIyPNDQ8Pp+WICa0SmmQSMq2aTZ4oL8e+qiq43W4ASKmBZqLHx8SSaDIXz5NmOttJaN08dTIxC4eBt98G/s//ybommplUWFiIuro6NDQ0zKvu50YgpdT98V6zZg3Ky8tx/PhxDA8P47XXXsOFF16oa0zziSZ9pKSUzwF4TotrUfrk5+enbfpfi4RG7+NMosnITLNKyZo+u/TTFSvgCYVmvXai3dCnL+2l8/zDWLRunrq4pQWmcHg8kQyHDX24sFGUlpYiEongxIkTKCws1P3FnSYkcGB8Km677TZ885vfnJyVYiKVOexsPo+sXr0af/rTn/QOY1Z6HmcyNRlJZFYp1vl6U02fXfr8oUMQUs66ZJdMMieEQG1bG849cADr3n4biqpm5IxCQPsk+EBxMTaZTFAAiCw4XFhXU16oKyZeqJlMZc6sM1JxHhivhU984hP4+te/DlVV8fTTT8Pv97O2MEOYSM0jZ555JvLz8zEyMgKHw6F3ODHNNJOSbsnWE8RTyzR1dkkKAUXKuJbsEl0iXNTRgZuffRbmcHhyl2AmzyjUKgkeGhpCs8WCEw89hBUdHVlzuLAuYrxQM5kykDgOjNdKaWkpLrzwQrz44oswmUx48skn8alPfSotY9Gp+Bs2j5jNZnzkIx9BZ2enIRscxrtrT4+eVTPtoJs62zRbo80XKivxl+pq/OSMM07b1aeVZW1tMEUik7/UKpCRMwq11NfXB5/Ph7vuugsrPv1p4M47mUTNJtYLNYCKigosXLgQfX19hvxdzyVSypnfiEUPjDeZTjswPh2iPaV8Ph/uu+++tI5F/8QZqXnm8ssvx4EDB7B9+3YsWLBA73BOEU/Bcjp2hwGzz0jNNus0Vy3T9PturarStP5qqqNVVYi8+SYQiUBVFLyxahXePPNMwx+vA4y3OTh58iQ8Hg++853vYPHixXqHlB2iL9TRGakpL9SVEz/3pqYmFqDrZYYD49Plfe9732SJwqFDh3D8+HH+LmUAE6l5RgiBm2++GW1tbWhubkZ1dbVhtsker6mBqigQE4lArJkULXaHJbJ8eMbAAD527NiMO+jmqmWKtfvu8cWLccjrnZzl0iqhOlFermuxfrKGhobQ09ODd7/73fjEJz7Buo5EzPFCzWQq/SKRCMzmWV5Kpx0Yn05WqxUf/ehH8fDDD0NVVTz00EPYsmVLRsaez5hIzUMOhwPf/va3cf/99+ONN97AggULsqb3TKzdYYkkRonMaE2dTVIARBB7OW62WqaZZqyS6RMVj5nqlPSqPZvN6Ogourq64PF48I1vfANr1qwxTFKfVeZ4oWYylV6RSMRQNae33HILHn30Ufj9fjz00EP43ve+x595mjGRmqccDgfuuOMOPProo3jmmWfg9XpRUFCga0yLW1qgTCQtUlVjzjZN3x0GIKGlvplmtBRFOe1FfOpsUhjA3qIiPLpkSUIJz0wzVqk0/ZzJTE+W6VoOTVYwGER7ezvsdjs2b96Myy67DDabTbd45oNoMtXY2IiCgoKseeOUDVRVTUt/vmStXbsWJSUl8Pv9CAaDePHFF3H55ZfrHVZOY5o6H23fDmzZAtPOnfjkJz+J//iP/4DdbkdjY6OuZ3bFe+RLU2UlXtqwAU2VlQkfTzPTGIqinNaLaerRLmFFOSWJSuT4lkNeL/YVFk6e4Tf12hEAEAKDFsuc15nLTImUVkf4pGpoaAhNTU3o6enBDTfcgAceeABXXXUVk6gMqaysRF1dHQYGBhAOh/UOJ2cIIWDR4PdXK0II3HbbbbDb7fD5fPjxj3+sd0g5jzNS802M7dIr6+tx9913489//jOefPJJRCIRlJeXp/VdVqylpmR6ESXaCHKmMWIlITPNJiW6LDfT7X+6YgU+d+gQFCnx2cOH0ZSfn9Ks1EzLYlo3y0xEJBJBV1cXgsEgKioqcMstt+D8889nHZROysrKYLFY8M4778DlchlqJiVbSSlnr5HSwSc+8Ql885vfhJQSf/3rXzEwMABvBg5Kn6+M9dOn9Juhr0leXh6uu+46bNy4EX/961/x7LPPIhwOo7S0VPMZg9mWmhLtRZRM8hVrjFgzUkDs+qdEl+Vmur0nFJqxp9RcTT5jmWlGSo+O8SMjI+jp6QEA1NfX44orrsDy5ctZA2UAhYWFOPPMM3Ho0CGoqsoZwRQZbUYKAMrLy3H++efjH//4B0wmEx577DF89rOf1TusnMVEar6ZZbs0MP4k+y//8i943/veh61bt+IPf/gDurq6YLPZUFRUpElthdbnsmnRCDKRYsxEj2+Z6fZaFqLPdd5XujvGSynh8/nQ398PAPB4PPjYxz6Giy++WPfaOzpdfn4+Vq1ahUOHDiESiXCGMEmqqkIIYciasy984Qt48803MTw8jB/96EdMpNKIidR8E2dfE7fbjWuvvRZXX3013n77bbzwwgvYvXs3pJRwu93weDxJzy7oudQ0EyHE5KzUXN9Xose3zHR7rQrRV/T3Y3l7OyIdHRktIldVFYODgxgaGoIQAuXl5fjwhz+MNWvWYOHChdwpZHAOhwNnnnkmjhw5Ap/PB5fLpXdIWSccDhtqx95UV1999WQz1sbGRhw8eBArV67UOarcJGItZ6TbunXr5K5duzI+LqVmcHAQu3fvxl/+8hc0NjZCCAGv14v8/PyEk6pMbsePd6wdO3bAZrPpngBEZ6TMEzNVs81ITb2tajandUeelBIjIyPo7++fXAZdvnw5LrroIqxevRqlpaVcustCoVAI77zzDvx+P+toEjQ0NISKigpUVVXpHUpMn/nMZ/DII49AURTceuut7HaeAiHEbinluphfYyJFiZJSoq2tDTt37sSrr76K5uZmKIoCk8kEr9cLu91umBfURLb+v/HGGzCbzYYoHI23RurDx4/jk0ePwgQgIgT+etFFeGnDBk1iUFUVfr8fw8PDCIfDkFKitLQU9fX1WL16NZYsWQK73a7JWKSvSCSCY8eOTRYlG+X31+gGBgawYsUKuN1uvUOJadeuXdi4cSP8fj88Hg96enoM8fyWjWZLpPiIUsKEEKiqqsK1116La6+9FkNDQzh27Bjeeust7NixAy0T2+vNZjPcbjecTqduT8yJ1GOZTKaYBeepSqZwPN4Di6fWWakpLJNOT5qiP6+amhps2LABy5cvx7Jly1BUVMQX2RxkMpmwdOlStLa2orW1FV6vly+4cZBSGnZpDwDOPfdcFBcXw+/3Q1VV/O1vf8NVV12ld1g5h78plDK3241zzjkH55xzDjZv3ozu7m6cOHECBw8exFtvvTU5YyWlhNPphMPhgM1my8gLciL1WIqiIBKJaDp+ujqYRx3yevGVtWuxZmAAoxO9teYyV9K0bNkyVFdXo7Ky0nC7kSh9FEXBggUL4HQ6cezYMdjtdu7om0UwGITD4TB0wimEwOc+9zl861vfwvDwMH784x8zkUoDLu1R2g0PD6OpqQlHjx7F4cOH0dTUhL6+vskZICkl7HY7HA4H7Ha75jVK8dZIHThwAKOjo8jLy9Ns7KlLb2EAv1q6FI9rfIhoMBiExWLB6tWrJz8npcTY2BgCgQACgQDGxsZOeVxramqwcuVKJk0Uk9/vx5EjRwBg3hWh2/bsgX3nToyuX4/A2rUz3s7o9VFRnZ2dqK2txdjYGPLy8tDW1obCOXYa0+m4tEe6im61XrVq1eTnxsbG0NXVha6uLrS1teH48eNobGzEyZMnAYy/kwqHw7DZbLDb7bDb7bBYLEnNYsW79V8IofnSXqKtEhJZBgyHwwiHw/D7/VAU5ZSZPyklvF4vqqqqUF1djerqahQXF6OkpIRJE83J6XRi1apVaGhoQH9//7ypm7Lt2YPqzZshgkFIqxWtjzwyYzKlqiry8/MzHGHiysrKcOGFF+LFF1+EyWTCb3/7W3zuc5/TO6ycwkSKdJGXl4eamhrU1NTg3HPPnfx8JBJBb28vOjs70dnZiRMnTuDEiRNob2+H3++fbFMAjM+6RE9et1qtsFqtsFgssFqtMJvNCT/xm81mzROpRFolrOjvx5ZduyaXAW8/80zsczoRiURw5tAQ1vl8eNPtxn6Xa7KbssfjgcfjwZIlS3DdddehqKgIhYWFKCgoYLJEKbFarVi+fPm8qpuy79wJEQxCqCoQCsG+c2fMRCraViBb+m994QtfwBtvvDG5vMdESlu5/VtBWcdkMqG0tBSlpaU466yzTvlaOByGz+eDz+fD8PDw5J+enh709vaip6cHfX19GBgYwMjIyClJF/DPJ7/ox+jXTSYTFEWZvF60yV70IOPon6jojM9UZwwMYM3gIPa43Xg7P3/yNlJKvC4EthcWQlVVyN7eybGnklJiaXs7zFP6R10iJQLLluHMoSF8cc8emCIRqJ2deOKWWzC4ciWsViuEEOjp6cGyZctw4YUXzvzAbt8+Z+8wounmW93U6Pr1kFYrEApBWiwYXb8+5u18Ph/Ky8sN2Ygzlve9732Tf29qamJPKY0xkaKsYTab4fV64+p1E026oslRIBBAKBRCOBxGKBRCKBSarB2K1hHt27cPR48ehd1uRzgcRiQSQSgUQiQSmTzkNZpUTU2yVg4O4gf79sEysXz3rYsuwrGSEphMJlgsFphMpsm2CjabDTabDVarFXl5eZMzaXl5eahsagL+3/+DGokAJhPy3/9+vKuuDqufew7mSASKlBCqimXt7dg35V1yJBKZ/cy0GOcrMpmiRBQVFcFms+Ho0aMYGBhIqSGvkQXWrkXrI4/MWSMVDodRUlKS4eiSZ7FY8MlPfhIPPvggwuEwHnroIdxzzz16h5UzmEhRTkok6Yp6+eWX8dJLL00Wj5Y2NKD8yBF0LF+Orrq6Ge+3+rnnYN21CwoAAeBDxcXYd+WVCcfcs3Qpnr/jjtPG7Fi+HKrZDEQiUE0mdCxffsr95jwvbYbzFYkSEa2bam5uRmdnJ9xud04uHwfWrp21yHx0dBQejyfreqjdcsstOPDzn+P8sTHs/vnPEf7hD3N+qTZT+CgSgPGlpWAwiFAoNPnRZDJN7qTLxXef0+Xl5U22PyhtaMCVd98NJRyGajbj+TvumDGZmivRiTchA4CuurrTbtNVVxczwYpSVXX2Gak5zlckipfZbMbixYvh9XrR0NAAk8k073b1BQIBLFy4UO8wEnaWz4fngkFYAAR9Puy6/36c/6Uv6R1WTmAiNQ9FIhGcOHECDQ0NOH78OE6cOIHW1laMjY2dVhMUrSdyOBwoKSnBkiVLsHTpUlRWVqKqqgoej0fPb0VTVqt1svap/MgRKOEwlIlGnuVHjsyYBM2W6ExPyF7/yEeQ5/fHlVRNH2Om20cikdlbNsR5viJRvAoLC+F0OtHY2Ije3l54PJ55MbsRCoVgsViyYrfeabZtQx4ABYAEcORnP2MipZHc/5+fa5IsGg6FQti7dy+2b9+OnTt3YmxsDKqqIi8vDw6HA8XFxTM+EUopEQ6HMTQ0hJdffhkvvPDC5Db7iooKbNiwAWvXrsXSpUtnnxkxuKntFeaaZZpupkTnlIQsHEb9o49CSDnnLFcipJRzP+719UygSFN5eXlYtmwZuru70djYCLPZnPOzU0NDQ1i2bJnu53EmZeNGCJsNodFRhAD8oqEBH+jvR0FBgd6RZT0mUtkkiaLhoaEhvPTSS/jDH/6AoaEh5OXlobCwMKGERwgBi8Vy2jsxKSX8fj+eeeYZ/PGPf4TJZMLFF1+Myy67DMuXL8+6J5upieRcy2nxmpqQSSEgIhEowJyzXKUNDVjy2muQABouuGDO8efDbAAZjxACpaWlcLvdOHHiBPr6+pCfn5+TtVM+nw+FhYXZ28yyvh5i61b85tOfxk+PHMF+qxWPPfYYbr31Vr0jy3p89s0mCRQNj42N4dlnn8UTTzyBSCSC4uJi1NbWahqOEAIul2vyXWg4HMYrr7yCl156CQUFBbjyyivxrne9C8XFxZqOmy7Tk5HZltMS8c4FF0AA6F2wAOc/9tics1ylDQ246oc/hDKxU3DZK6/g+a9+dcZYhBBMpEhXNpsNK1asQF9fHxobG+H3++F2u7PuzdRMos1vFy5cmN31ovX1KL/vPhy44Qb4h4dx//33M5HSAJ99s0mcRcNvv/02HnjgAXR2dqKioiIty221bW1Y0tqKY9XVk13DzWYzKioqAAAjIyP4zW9+g9/+9rd4z3vegw984AMoKyvTPA4tJdPEczbT66OOXXBBXLNc0eXAaCSKqs46e8VEioxACIGioiJ4PB60t7fj5MmTsFqtObHcNzQ0hMWLF2t6fJReLr/88skZw8bGRhw+fBgrVqzQOarslhtvF+aLaNHwd74Tc1kvFArhF7/4Be666y4EAgHU1tamLYm69amncMVrr+HWp55CbVvbabdxOBxYsGABKioq8MILL+Bzn/sc/uu//gvt7e2axwNgfNlzy5bxj0nS8kmytKEBa595ZrI+SpmylLfvqqtmnemKLgdKjBeFqooyZ41WNtemUW4xm82oqanB2WefDafTid7eXgSDQb3DStrw8DA8Hk9W9Y2ajclkwr/+67/CarUiHA7j5z//ud4hZT2+jc02MxQNDw8P4+6778b+/ftRU1OT1o67S1pbYYpEYJrY0baktXXGs+zMZjOqqqoQDoexbds2bN26FZdccgmuu+467Q771KjhpM1m0+SImMmZqFAIAoAqRFwF61FddXV47itfSahGKpe7TVN2stvtWL58OQYHB3HixAmMjIzA5XJl1eypz+eDxWLBkiVLsntJb5qbbroJ999/PwKBAB5++GH84Ac/yJou7UaUPf+jaUb9/f341re+hY6ODtTW1qb9F/5YdTUiJhMQiSBiMuFYdfWc94kmVJFIBK+++ir+8Y9/4MMf/jCuueaa1AtTNWo4qVWDvcmdehhPotpWrsSeD3xAs3YH00kpmUiRYXk8Hpx11lno6elBa2srwuEwXC6X4QvSfT4fzGYzzjjjDMPHmqilS5di2bJl2LdvHyKRCF544QW8973v1TusrMWlvSzn9/vxve99D11dXaiurs7Iu6amykr814c+hL9ccAH+60MfmnE2KhaTyYTKykpUVFTgt7/9Le688040NzenFlC0dsxkSqnhZPTsumjvrGRFl+ZURYFqNiecRCUj27os0/xiMplQVlaGNWvWYPHixQgGg+jr6zPskl8uJ1FRt99+O1wuF4aHh/HAAw/oHU5WE1qfdh+PdevWyV27dmV8XENK4TDZYDCIH/zgB9i3bx9qamrSEl66dXV1YWxsDB//+Mfxvve9L/lpf40O5b377rvhcDhSrpdKpJt5KiKRCDo7O/HNb34zbWMQaU1VVQwMDKClpQWjo6NwOByGmFWVUmJgYAAOhwMrVqzI2SQKGC8HKS0tRSAQgM1mQ3t7e0JHas03QojdUsp1sb7GpT09pVjb86tf/Qp79uzRvK1BJpWWliIYDOKRRx7B9u3b8YUvfCG52imNGk66XC6MjY2lnEhp1TphLqFQKCd2RdH8oigKCgsLUVBQgKGhIbS0tKCvrw8mkwlOp1OXOqpAIAC/34/q6mpUVlbmfM1Qfn4+3v/+9+PJJ5+EyWTCE088gZtvvlnvsLISl/b0FKu2J0779+/Hs88+i5qamqwvgrRarVi4cCFaWlrwpS99Cc8999zkmXeZ5na7EQqFdBk7GUykKJsJIeDxeLBq1SqsXr0alZWVGBsbQ19fH4aHhzPyPBAKhdDX1wcAWLVqVdo36xjJ5z73OTidTvj9fvz4xz/WO5ysxRkpPSV5mKzf78ePfvQjFBYW5swvvBACZWVlCAQCeOihh7B9+3Z8+ctfzvjxBW63G01NTRkdMxWhUAhut1vvMIhS5nA44HA4UFlZCb/fj76+PnR1dSEcDsNisSAvL0+zNh+RSAR+vx/hcBh5eXlYsmQJCgsLc6aBaLze9a53weVywefz4ejRozh+/DgWL16sd1hZZ379rzGaOfpCzeR3v/sd+vv7c+rA4CibzYaFCxfi6NGj+NrXvpa+vlMz8Hg8WTcjlZUHqBLNIHpiwoIFC3DOOedg5cqVKC0tBTC+Qzn6x+fzIRwOz9myREqJUCiE0dFRDA0Noa+vD36/H6WlpTjrrLOwZs0aFBcXz7skChh/rD//+c9Ptn757//+b71DykosNs8yvb29uPXWW1FeXq5ZHUGsLuVG0NXVBbPZjLvuuitj75L27t2L3//+96fVnWWqeDxRra2tuPzyy3HhhRfqHQpR2kUiEYyNjWF0dBSDg4MYHh5GIBAAgMkSBynlaeUONpsNdrsddrsdHo8HLpcr60sitNLe3o5FixZhbGwMFRUVOHnyJB+bGFhsnkP+/Oc/A9DukNpol3LTRE+oRNsZpFNpaSn6+vrwta99Dd/4xjewatWqtI9pt9tPexKZftTL83fcYZhkKhKJsEaK5g2TyTS5BFhUVDT5eVVVJ/9EIhGoqgpFUWCxWGAymZgYzKKiogIXXXQRtm7diqGhIezcuRMbNmzQO6ysMv/mMrNYf38//vznP2t6Zt3ULuWmiS7lRlJYWIj8/Hx8+9vfxmuvvZb28Ww222lPupMNNqcc9WIUiqIYYts4kZ4URYHZbIbVaoXdbofT6YTdbtf8/Mxc9cUvfhH5+fkIBAI8MiYJTKSyyM6dOycLL7US7VIeESLuLuWZlp+fj+LiYvzwhz/EX/7yF02OcZlJrMaWpzTYTOCol0xhIkVEqbjiiitgNpsRiUTw2GOPZVWdqBEwkcoif//73zVvmJZMl/LatjZs2rkz5mHF6RLdzfPTn/4Ujz32WNqSqVhJSVddHZ6/4w7svuYaQy3rRTGRIqJUmM1m3HzzzZOnO/zlL3/RO6SswkQqS7S3t6OxsTEtW92bKiuxdf36uJOoW596Cle89hpufeqpjCZTeXl5WLBgAR5//HH8/Oc/T0syZbfbY163q64O+666KqEkqrShAaufew6lDQ1puX1U3MfDbN8ObNky/pEoW/H/cVrccsstUBQFw8PD+OlPf6p3OFmFxeZZYu/evTF3o2Ta1JoqTNRUJVOcnuxOQbPZjIULF+K5555DQUEBrr/++oTHnk20ODUSiaTUoyvRAvVkC9rjPrA4xS76RIbA/8dps2jRIpx99tnYsWMH/C+8gMC3vgXbFVfw8Y0DZ6SyxN69e+F0OnUbP7qc57PZUq6pSnVWS1EUVFdX49e//jW2JdANPl75+fkp1wgkWqCeTEF7dGdSXE0KU+iiT2QY/H+cVl/60pewyeHAc8EgrN/5znjSypm/OXFGKgtIKXHw4EHdGnBOb5Hw9CWXwBUIJN13SotZLYvFgsrKStx///0oKirCWWedlXAcM8nPz4fP50up9ihaoI5IJK4C9URvDyTYjDPJLvpEhsL/x2l1zTXX4NAnPwkrAEXKfyarnJWaFROpLNDV1YVAIICSkhJdxp+e+LgCAWxdvz7p60V3CmIiMUt2p6DNZkNhYSG2bNmC73//+1iwYEHSMU2Vn5+P/v7+lK4RLVCPt4lnorcHEkykol30t20bf/HhEyNlI/4/Tqu8vDy4rr4awd//HhKAyWKBwmR1TkykskBra6uutVFaJT5R0Z2CWnRTz8/Px9jYGL7//e/jBz/4gSbHpWh1cHFXXV1CxemJ3j4YDCbWU6y+ni88lP34/zitrvyP/8BVzz6LiyMR1H3yk9jMx3pOrJHKAgMDAxk5BX0mybRIiOea8e4UnEtxcTF6enpw//33a/I4eTweBIPBlK+Tbjxnj4i0duaZZ6J32TJ8NxLBXc8/n9a+fbmCiVQW6OrqimsHmRb9nWa6hpaJT6pixVhVVYU33ngDv/vd71K+vtvthqqqKV8n3YLBIAoKCvQOg4hyzJe//GU4nU709fVhz549eodjeEykskBHRwfy8vJmvY0W/Z307BE1W0xTk6aZYhRCoKamBo8//jjmPBB7jj40brc7K06Cl1KisLBQ7zCIKMfccMMNUFWVR8bEyfivFoTe3t45t7hrcWbe9GusO3gw4x3Mp4qVNM32fZrNZpSUlODHP/4xfD5f7ItG+9B885szbu3Nz8/PyHR2sg04o4QQXNojIs25XC5cf/31kFLiN7/5ja6lJdmAiVQWCIfDcxaba3Fm3tRrqIqC9QcP6jo7FStpmuv7dLlc8Pv9ePzxx2NfNI4+NNHkJJ3JVLQB57lPP40r7747qWRKSslEiojS4vbbb4fdboeqqnjxxRf1DsfQmEhlgXhe0LUoCJ96jZ0rV0JR1fEkJhzGe19/PePJVKykKZ7vs6KiAn/+85/xzjvvnH7RaB8ak2nGPjRmsxlutzutBefJNOCcSlVVCCHgcrnSFCERJSyHjq8599xzUVlZieHhYS7vzYHtD7JAvDMjTZWVKReDR69R29aG8w4dAsJhKACWNTVh8cmTmu3aizeWWG0S5vo+zWYz8vPz8cADD+CHP/whLBbLP78YZx+akpISdHd3z1mblqxkGnBONTY2hoKCgqyo5SKaF3Lw+Jovf/nL+NKXvoQ//elPGB0djf9cz3mGz8IUUzSJOVpbC4nx/yjJ1l6lGkcyuwWLiorQ3NyM55577vQv1tcDd94565NcaWkpAoFAouHGLdqAc/c118R9rt5UejZoJaIYcvD4mo9+9KMAxt/MP/vsszpHY1xMpLKAw+FIaTt+sm0Rmior8dfzz0fYbJ619kqLtgvxSmSsiooKPProo+jo6Eh4nNLSUoyNjaVcED6brro67LvqqoSTKGA8kSotLdU8JiJKUhxlA9nG7XbjQx/6EAKBAB588EG9wzEsLu1lAa/Xi4YkX8inn5OX6NLcXF3IU71+IhIdKy8vDyaTCQ899BC+8Y1vJNQd3uPxoKa1FVf++tdQwmGoZnNSM0fpEgwGOSNFZCQ5enzN7bffjt/97nd4+eWX0dfXx5YrMXBGKgsUFxdjbGwsqftq0RZhtuU1La4fr2TGKi8vx5tvvonXX389obHcbjcWNTWlVBAer2RmvRRFgdvtTks8RJSkOMoGss26deuwYMEChEIh/PGPf9Q7HENiIpUFysrKku7jEWvnm5ZLcdOv77PZ0rbMl0yLByEESkpK8OCDDyaUjObn5+NEbS1UsxmqoiRVEB6PVNogsPUBEaWbEAL//u//DiklHn74Yb3DMSQu7WWBwsLCpHdnTV+aA6DpUtzU6/tsNlz7j3+kbZkv2cOOXS4XmpqasH37dmyMs27BZrOhe8kS/OmLX0R1QwM6li9Py7Le1DYImJj1imccVVWZSBFRRnz0ox/FrbfeildffRX9/f08mmqalGakhBA/FEIcFkLsE0I8LYTwahQXTVFRUZHS/acuzaVjKS56fVcgkPZlvmR38RUUFODJJ59MqGi/pKQETZWVSReExyPaBiGRWa9QKASn05m21gxERFPZ7XbcdtttAIBnnnlG52iMJ9Wlvb8DWCWlXA3gHQB3ph4STVdWVgYhhCZt+rXogK7HtVPldrvR3t6Ot99+O+77pLsFApBcG4RAIIDi4uK0xkVENNUXv/hFAMBPfvITfQMxoJSW9qSUf5vyz9cBXJ9aOBSLyWTCggULMDAwkHKBcTy78BJdOov32npzOBx4+umnsXr16rhuX1ZWlpGTz7vq6hKa8RodHcXixYvTGBER0akWLlyIDRs2YMeOHRgYGIDX69U7JMPQstj8/wB4XsPr0RTLly+H3+/X5FozLY/FOiRYq2sbQXFxMfbt24empqa4bu/1euPuKp/OflPTBQIBlJWVpX0cIqKpvvvd7wLg8t50cyZSQogXhBBvx/jzwSm3+TqAMIBHZ7nOzUKIXUKIXd3d3dpEP4+sWrUKoVAorWNkspVBPDbs24ebf/97bNi3T5PrCSFgNptjdzuPwe12x9V7SosDiBPFd4NElGmbNm0CgMl6KRo3ZyIlpbxcSrkqxp8/AoAQYjOAqwF8TM7y9l1K+TMp5Top5To2EkzcGWecASD+c/eSYaQapw379uHGrVuxvKkJN27dqlkyVVZWhpdeegn9/f1z3jbeZdRUDyBOlBCCO/aIKOOEEPje976n2epIrkh1194VAL4K4ANSyhFtQqJYvF4vamtrMTw8nLYxojVOf7nggoweThzL2ceOAQDEtH+nymw2Q1VVbN26dc7bOp1OKIoy506/ZHbepUJKyWacRKSLO++8M+kG0bkq1RqpHwPIB/B3IcReIcRPNYiJZnDRRRdhcHAwrWMYpcbprSVLAABy2r+1UFpaimeeeQbhcHjW2wkhUFxcjJGR2d8jTN95ByBt9VLhcBgmkwkOh0PzaxMRzUUIAavVqncYhpLqrj3tXt1oTmeddZbeIcwold1+se67Y2Jn3dnHjuGtJUsm/60Fm82G7u5uHD9+HMuWLZv1tgsWLMC+ffvgcrlmvV105120Xipd5/P5/X5UVVUl3aCViIi0xc7mWWTRokXIy8tDMBg01DuCVA4unu2+O1av1jSBmm7v3r1xJVI7d+6M+5rJdiqPl9/vj7t9AxERpR/f1mYRk8mEDRs2oLe3V+9QTpHKbj+9dgp6vV68+uqrc96uqKgorp17UYnUSyXTMiEcDqPSgK0liIjmK85IZZnLLrsML730EqSUCb3Aa2Gm5bvobj9MzCpFd/vFs9w3033TzeVyobm5Gb29vSgqKprxdkVFRQntlIzWS5UfOTLr+XzJLgEKIWaNl4iIMouJVJZZuXIlqqurMTg4iLNHRjLSRby2rQ3rDh7E+oMHoajqaUtwsTqaR5fszOEwpKLgyUsvjblMp1c3dCEEhBA4dOgQLrroohlvZ7PZ4PV6EQgEYLPZ4rp2PJ3Kk1kCVFUVUkoUFhbGFQcREaUfE6ksI4TAjTfeiD9/4xu49dVXk6pLSsTUhEhgoh3BxBLc1PGaKitP+feS1laYw2EoAKSq4kMvvYSO4uKYMU6/byqF64mw2Wx4/fXXZ02kgPE6qYaGhrgTqXhElwARicTdMiEQCKCkpAQWi0WzOIiIKDVMpLLQ+vXr0T48PFlbFCux0Uq0hknBeCsCFYhrCe5YdTWkokCq6ngCpqpxxZhK4XqiCgsLsXv3bqiqOusuuNraWuzfv1/TseNdApzK5/PhzDPP1DQOIiJKDYvNs5DVakXpDTcgrChp70I+tdt52GTC9rPOiiu5aaqsxJOXXoqIoiACIGI2xxVjporPa9va8N4338Tizk7MdWRRSUlJWurRuurqsO+qq+Le1RcIBLBgwQLN4yAiouRxRipLrbn1VvzH3/+O80dHcXzBgrTN2qRSw7Rj9Wp0FBfPet/py3iZKD6fOut1uaKg8fnnUbZ584y3LykpSevRPPESQoDHKxERGQsTqSzl8Xiw8CMfwa9feAFVaS7Qnl7DpNV9Z1rGS3fx+dRZLxmJILJ1KzBLIuV0OuHxeBIqONdatNCciRQRkbFwaS+LXXnllYhEInOeBWdUMy3jpfuYmumHM++ao2s5ANTV1WFoaCgt8cRjZGQE5eXlhmrESkRETKSyWnV1NdauXYuuri5d46hta8OmnTtR29aW0P2mJzSZ6iE19XDme6++Gq/HUf+0ePFijI6OZiC62IaHh7F48WLdxicioti4tJflPvrRj+KOO+5AOByG2Zz5H2cqu+yaKivx9CWXTJ6nl8kWCNElx0gkgt6OjjkbnJaUlOh6vl0oFEJtba1u4xMRUWyckcpyixYtwgc/+EG0t7frMn4qu+xq29pw7T/+gaXNzbj2H/+YnNGKJmdXvPYabn3qqYRnumYaK9asmclkQjgcnnO2qbi4GEIIXZdRWR9FRGQ8TKRywA033ACPx5P2Gp5YyUgqy3MzJWFat0CYKzFTFAV+v3/Wa5jNZlRVVcHn86UUSzLC4TAsFgsKCgoyPjYREc2OiVQOsNvtuO2229Db25u2GZOZkpGp9UaJNs/02WyAEFCnJWFa107NlZgpioKRkZG5r7NkiS4F58PDw6itrdV1aZGIiGLjM3OOWLt2LS688MI5l/iSLQyfLRlJZpdddFkPqgpVCDx9ySWnnd2XTHIWSzyJWTx9ohYtWqTL0t7w8DBWrlyZ8XGJiGhuLDbPEUIIfPrTn8abb76JisZGrOzqOq1QO5XCcK0bZU4mZgAiUsIVCJzy9VR6V00XT2+qeDqXV1ZWTtZUZbqwn4XmRETGxEQqhxQVFeH29eux5o47YFbV05KlqbNKiZ7Pp3WjzEx0MJ9qtsQs3q7lFosFy5cvR1NTE4qKirQMb0ajo6Nwu90ZG4+IiBLDRCrHnOf3Q6pqzGQp1eQlnlmieNsWZKKDebwSOUdv1apVOHDgQMYSm76+PlxwwQVpOeuPiIhSx0QqxyiXXQb1u99FOBCAOi1ZSnfykujSoZbLd6lQVRUOhyOu29bU1ADAnH2ntBIOh7FkyZK0j0NERMlhIpVr6uuhvPgijv3sZ3i4oQFj5eVQcOpM0db169MydCpLh3qRUiISicDj8cR1e7fbjfLycvh8PuTn56c1tmgtVlVVVVrHISKi5HHXXi6qr8fSX/wCCz78YbS0tGDByZOaN7iMRa8jX1IRDofhcDgSOsNuzZo1GBgYSF9QEwYHB7Fs2TJYLJa0j0VERMnhjFSOEkLgM5/5DBobG1H+t79lZKYoHUuH6T4qZmxsDMXFxQndZ9GiRXEXqKfC5/PhzDPPTPs4RESUPCZSOSwvLw9f/epX8dMDBxA+cACY2MmXzpkiLeueUmnXEC+/359wj6aysjLYbDYEg8GEZrKSsWDBgrRen4iIUsOlvRxXUlKCj9x3H/7jkkvwzLp1aUlG0kXro2JiCQQCWL16dUL3URQFZ511Fnp7ezWPJ2p4eBhlZWVx124REZE+mEjNA0uXLsU13/8+flVdjXeyqB9RJmquFEXBwoULE77f6tWrEQwGNY8nqr+/Hxs2bEjb9YmISBtMpOaJ8847DzfddBNOnjyJcDisdzhx0fqomOkikQiEEKhOIkGrrq6Gx+OJ64y+REWPoVmxYoXm1yYiIm0xkZpHrrjiClx33XVobm7OqmQq0XP84tXX14fVq1cnVeekKArq6+vTsrzX29uLFStWwOVyaX5tIiLSFhOpeUQIgY9//OP40Ic+hObm5rQuTU2X7GHJ6TQyMoLLL7886fuvXLkSUkrNDzL2+/0477zzNL0mERGlB3ftzTOKouDjH/84vF4vHn74YVRUVMBut6d1TC1232ndBiHa7DLRQvOpvF4vlixZgvb2ds2OjBkbG4PD4UiqbouIiDKPM1LzkBAC73//+/HlL38ZXV1d8Pl8aR0v1d130URMy4ai3d3dOP/88+M+GmYmF1xwAYaGhlKOJ6qzsxMXXXQRzGa+xyEiygZMpOaxd73rXfjmN7+JwcHBtHbqTnX3ndZtEFRVxdjYGK655pqUrgOMN+esqqrS5PELBoMwm80499xzU74WERFlBhOpeW7NmjX43ve+h3A4jJ6enrSMkeruu8lEDIAUAj6bLaV4uru7sW7dOixatCil6wDjS6WXX365JolUZ2cnLr744pRnyYiIKHOYSBGWLFmCLVu2wGazob29PS1jpLL7rqmyEk9fcgmgKFBUFdf+4x9JL++pqorR0VF8+MMfTur+sdTV1aGsrAyDg4NJXyMUCkEIgXXr1mkWFxERpR8TKQIAVFVVYcuWLSgvL0dTU5Ph2iO4AgFASihASst7J0+exKWXXoq6ujrNYlMUBe9973vR29ub9A6+trY2XHzxxWx5QESUZZhI0aSioiJs2bIFH/jAB9Da2prWuqlEadHlfGhoCC6XC5s3b9Y8vqVLl2LNmjXo6OhI+L4DAwMoKirCRRddpHlcRESUXkyk6BQ2mw2f+tSn8N3vfhcWiwXNzc2IRCJ6h5VynVU4HEZvby/+7//9v8jPz09LjFdeeSXMZjP8fn9CcQ0MDOC6665L+wHIRESkPSZSFNPKlStxzz334IorrkBLS0tK9T9aSbbOSlVVNDc347rrrsOaNWvSExwAl8uFa665Bp2dnQiFQnPeXkqJlpYWXHzxxaipqUlbXERElD5MpGhGDocDn/nMZ/Dtb38bQgg0Nzdr3sU73aSUaG5uxmWXXYaPfexjEEKkdbwzzjgDV111FVpaWmZNplRVRVNTE84++2xcdtllaY2JiIjSh4kUzWn16tW49957cemll6KpqclQtVOziSZRa9aswWc/+1koSmb+u1944YW4+uqr0dzcjOHh4dO+HggEcOLECaxZswbXXXcdm28SEWUxIaXM+KDr1q2Tu3btyvi4lBopJfbs2YNf/OIXaGtrg9frhcfjSfssTzLC4TBaWlpwwQUX4Atf+AJsKfaeSsbhw4fx17/+FT09PcjLy4MQAoFAAHa7HZdeeinWr18Pk8mU8biIiCgxQojdUsqY/WmYSFHCIpEI3nzzTTz66KNobm6Gx+OB1+s1TEI1PDyMnp4e3Hjjjbjxxht1TVZUVcWxY8fQ3NyMcDiMsrIyrFy5Enl5ebrFREREiWEiRWmhqir27t2LRx99FCdOnIDL5UJRUZFuCVU4HEZbWxs8Hg9uu+02HrVCRESamC2RYnEGJU1RFJxzzjlYs2YN9u/fj9/97nc4dOgQFEVBUVFRxo46CYfD6OrqQigUwgc/+EFcf/31PGaFiIgygokUpUxRFJx99tk4++yz0dHRgcO//CV6f/977Ha5cLS4GIWFhXA4HHHNVNW2tWFJayuOVVfP2eYgEAigu7sbAHDJJZfg/e9/P2prazX5noiIiOLBRIo0VX7iBMq/8x3IYBDXWiz421e/ij90dqKlpQVSSlitVjidTjgcjtN2q9W2teHWp56CKRJBxGQ6rfFmOByGz+fD0NAQAMDtduOGG27A5ZdfjsLCwox+n0RERAATKdLatm1AMAgRiUAAuMJmw3v/67/Q3d2No0eP4uDBgzh69Ohk8bWiKJBSQlVVnHPkCEzhMEwAZDiMov378fLEbVRVhcViQV1dHa6//nqsWrUKNTU1hilwJyKi+YmJFGlr40bAagWCwfGPGzdCCIHS0lKUlpbiwgsvBDBeqN7T04Pu7m74/X6MjIzA/MYbkN/4BtRwGMJsxvJbbsHta9eipKQEZWVlhtoZSEREBHDXHqXD9u3jM1MbNwL19Zm7LxERURqw/QGR3pggEhFlLbY/INLT9u3Apk3/XO7cupXJFBFRjuBZe0TpNlGAj0hk/OO2bXpHREREGmEiRZRu0QJ8k2myAJ+IiHIDl/aI0q2+fnw5jzVSREQ5h4kUUSbU1zOBIiLKQVzaIyIiIkoSEykiIiKiJDGRIiIiIkoSEykiIiKiJDGRIiIiIkoSEykiIiKiJGmSSAkh/k0IIYUQxVpcj4iIiCgbpJxICSFqALwHQHPq4RARERFlDy1mpO4B8FUAUoNrEREREWWNlBIpIcQHAZyUUr6lUTxEREREWWPOI2KEEC8AKI/xpa8D+BrGl/XmJIS4GcDNALBgwYIEQiQiIiIyJiFlcityQoizAGwFMDLxqWoAbQDWSyk7ZrvvunXr5K5du5Ial4iIiCiThBC7pZTrYn0t6UOLpZT7AZROGaQRwDopZU+y1yQiIiLKJuwjRURERJSkpGekppNSLtTqWkRERETZgDNSREREREliIkVERESUJCZSREREREliIkVERESUJCZSREREREliIkVERESUJCZSNH9t3w5s2TL+kYiIKAma9ZEiyirbtwObNgHBIGC1Alu3AvX1ekdFRERZhjNSND9t2zaeREUi4x+3bdM7IiIiykJMpGh+2rhxfCbKZBr/uHGj3hEREVEW4tIezU/19ePLedu2jSdRXNYjIqIkMJGi+au+ngkUERGlhEt7REREREliIkVERESUJCZSREREREliIkVERESUJCZSREREREliIkVERESUJCZSREREREliIkVERESUJCZSREREREliIkVERESUJCZSREREREliIkVERESUJCZSREREREkSUsrMDypEN4CmjA+svWIAPXoHYUB8XGLj4xIbH5fY+LjExsclNj4usWn1uNRKKUtifUGXRCpXCCF2SSnX6R2H0fBxiY2PS2x8XGLj4xIbH5fY+LjElonHhUt7REREREliIkVERESUJCZSqfmZ3gEYFB+X2Pi4xMbHJTY+LrHxcYmNj0tsaX9cWCNFRERElCTOSBEREREliYkUERERUZKYSGlACPFvQggphCjWOxYjEEL8UAhxWAixTwjxtBDCq3dMehJCXCGEOCKEOCaE+He94zECIUSNEOIlIcRBIcQBIcTtesdkJEIIkxBijxDiWb1jMQohhFcI8eTEc8shIUS93jEZgRDiSxO/Q28LIX4rhLDpHZNehBC/EEJ0CSHenvK5QiHE34UQRyc+Fmg9LhOpFAkhagC8B0Cz3rEYyN8BrJJSrgbwDoA7dY5HN0IIE4CfALgSwEoA/yKEWKlvVIYQBvBvUsqVAM4H8Dk+Lqe4HcAhvYMwmPsA/EVKuQLA2eDjAyFEFYD/C2CdlHIVABOAj+gbla4eAXDFtM/9O4CtUsqlALZO/FtTTKRSdw+ArwJg1f4EKeXfpJThiX++DqBaz3h0th7AMSnlcSllEMBjAD6oc0y6k1K2SynfnPj7MMZfFKv0jcoYhBDVAN4H4Od6x2IUQggPgHcBeBgApJRBKeWArkEZhxmAXQhhBuAA0KZzPLqRUv4vgL5pn/4ggP+e+Pt/A7hG63GZSKVACPFBACellG/pHYuB/R8Az+sdhI6qALRM+XcrmDCcQgixEMBaADt0DsUo7sX4mzNV5ziMZBGAbgC/nFjy/LkQwql3UHqTUp4EcDfGV0TaAQxKKf+mb1SGUyalbJ/4eweAMq0HYCI1ByHECxNrz9P/fBDA1wDcpXeMepjjcYne5usYX8J5VL9IyciEEC4ATwH4opRySO949CaEuBpAl5Ryt96xGIwZwDkA/ktKuRaAH2lYosk2E/U+H8R4olkJwCmE+Li+URmXHO/3pPnqkVnrC+YaKeXlsT4vhDgL4/953xJCAOPLV28KIdZLKTsyGKIuZnpcooQQmwFcDWCTnN/Nyk4CqJny7+qJz817QggLxpOoR6WUv9c7HoO4EMAHhBBXAbABcAshfi2lnO8vjq0AWqWU0VnLJ8FECgAuB3BCStkNAEKI3wO4AMCvdY3KWDqFEBVSynYhRAWALq0H4IxUkqSU+6WUpVLKhVLKhRj/RT9nPiRRcxFCXIHxpYkPSClH9I5HZ28AWCqEWCSEsGK8EPQZnWPSnRh/9/EwgENSyv+ndzxGIaW8U0pZPfGc8hEALzKJAiaeV1uEEMsnPrUJwEEdQzKKZgDnCyEcE79Tm8Ai/OmeAfCpib9/CsAftR6AM1KUDj8GkAfg7xOzda9LKT+rb0j6kFKGhRCfB/BXjO+o+YWU8oDOYRnBhQA+AWC/EGLvxOe+JqV8Tr+QyOC+AODRiTckxwF8Wud4dCel3CGEeBLAmxgvo9iDeXxUjBDitwA2AigWQrQC+BaA7wN4QgjxrwCaANyo+bjze9WFiIiIKHlc2iMiIiJKEhMpIiIioiQxkSIiIiJKEhMpIiIioiQxkSIiIiJKEhMpIiIioiQxkSIiIiJK0v8Psrgb7RZdOb4AAAAASUVORK5CYII=\n", "text/plain": ["
"]}, "metadata": {"needs_background": "light"}, "output_type": "display_data"}], "source": ["import bayespy.plot as bpplt\n", "fig, ax = plt.subplots(1, 1, figsize=(10, 7))\n", "bpplt.gaussian_mixture_2d(Y, alpha=alpha, scale=2, color=\"black\", fill=True, axes=ax)"]}, {"cell_type": "markdown", "metadata": {}, "source": ["We get the result of the optimization:"]}, {"cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [{"data": {"text/plain": ["(array([[-1.36058471e-01, -5.54785591e-01],\n", " [ 0.00000000e+00, 0.00000000e+00],\n", " [ 0.00000000e+00, 0.00000000e+00],\n", " [ 7.32146904e+00, 4.25483005e+00],\n", " [ 0.00000000e+00, 0.00000000e+00],\n", " [ 8.75857358e-02, 9.11669540e-03],\n", " [ 3.72250304e+04, 1.40121275e+04],\n", " [ 0.00000000e+00, 0.00000000e+00],\n", " [ 0.00000000e+00, 0.00000000e+00],\n", " [-7.71374280e+00, -7.35395097e+00]]),\n", " array([[[ 6.52316512e-01, -1.63828746e-01],\n", " [-1.63828746e-01, 8.96234523e+00]],\n", " \n", " [[ 2.00000000e+05, 0.00000000e+00],\n", " [ 0.00000000e+00, 2.00000000e+05]],\n", " \n", " [[ 2.00000000e+05, 0.00000000e+00],\n", " [ 0.00000000e+00, 2.00000000e+05]],\n", " \n", " [[ 1.28998747e+00, 6.13731151e-01],\n", " [ 6.13731151e-01, 6.90986759e-01]],\n", " \n", " [[ 2.00000000e+05, 0.00000000e+00],\n", " [ 0.00000000e+00, 2.00000000e+05]],\n", " \n", " [[ 8.65297348e+00, 9.29871380e-03],\n", " [ 9.29871380e-03, 3.55806853e-01]],\n", " \n", " [[ 4.64596327e+03, 1.74893727e+03],\n", " [ 1.74893727e+03, 6.60113981e+02]],\n", " \n", " [[ 2.00000000e+05, 0.00000000e+00],\n", " [ 0.00000000e+00, 2.00000000e+05]],\n", " \n", " [[ 2.00000000e+05, 0.00000000e+00],\n", " [ 0.00000000e+00, 2.00000000e+05]],\n", " \n", " [[ 2.88152054e+00, 9.46475241e-01],\n", " [ 9.46475241e-01, 2.54402981e+00]]]))"]}, "execution_count": 15, "metadata": {}, "output_type": "execute_result"}], "source": ["from bayespy.inference.vmp.nodes.gaussian import GaussianWishartMoments\n", "parent = Y.parents[1]\n", "(mu, _, sigma, _) = parent.get_moments()\n", "mu, sigma"]}, {"cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [{"data": {"text/plain": ["array([[-0.22515766, -0.06601764],\n", " [ 0. , 0. ],\n", " [ 0. , 0. ],\n", " [ 4.75563304, 1.93368381],\n", " [ 0. , 0. ],\n", " [ 0.01009479, 0.02535878],\n", " [ 8.21785907, -0.54595493],\n", " [ 0. , 0. ],\n", " [ 0. , 0. ],\n", " [-1.96797779, -2.158508 ]])"]}, "execution_count": 16, "metadata": {}, "output_type": "execute_result"}], "source": ["import numpy as np\n", "mu2 = np.linalg.solve(sigma, mu)\n", "mu2"]}, {"cell_type": "markdown", "metadata": {}, "source": ["The way you can build your model is quite nice but it still needs some development. scikit-learn proposes a better interface."]}, {"cell_type": "markdown", "metadata": {}, "source": ["

scikit-learn

"]}, {"cell_type": "markdown", "metadata": {}, "source": ["We try to solve the same problem with another module: [scikit-learn](http://scikit-learn.org/stable/auto_examples/mixture/plot_gmm.html#example-mixture-plot-gmm-py)."]}, {"cell_type": "code", "execution_count": 16, "metadata": {"scrolled": false}, "outputs": [{"data": {"text/plain": ["GaussianMixture(n_components=10)"]}, "execution_count": 17, "metadata": {}, "output_type": "execute_result"}], "source": ["from sklearn import mixture\n", "gmm = mixture.GaussianMixture (n_components=10, covariance_type='full')\n", "gmm.fit(X)"]}, {"cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [{"data": {"text/plain": ["BayesianGaussianMixture(n_components=10)"]}, "execution_count": 18, "metadata": {}, "output_type": "execute_result"}], "source": ["dpgmm = mixture.BayesianGaussianMixture(n_components=10, covariance_type='full')\n", "dpgmm.fit(X)"]}, {"cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [{"data": {"image/png": "iVBORw0KGgoAAAANSUhEUgAAAzEAAAGrCAYAAADw9BmMAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAAsTAAALEwEAmpwYAABi5ElEQVR4nO3deXxcV33//9eZGWm0jmTZkm3JWzaTBIUk2ISwhSVhcwIBWkAQKAmlNqQLgVBooPRL6Y+2FMpWoFgFUiApooUAwYQtLSRNE0jikMXZnMV2bMu2bMvyyFpGmpnz++PM2JKsZZY7y515Px8PPcaeuXPnjKPMmc/9nM/nGGstIiIiIiIifhEo9QBERERERESyoSBGRERERER8RUGMiIiIiIj4ioIYERERERHxFQUxIiIiIiLiKwpiRERERETEVxTEiIiISNUxxnzNGPPxeR5/iTHm8QzOc6Ux5o55Hv+NMeY9uY5TRGanIEYEMMb0GGN+Z4wZMcYMpP58tXH+3RhjjTGXz3jO51P3X5n6+5Wpv39+xnGXp+7/9+K9IxGR6mWM2WmMGTPGDBtjhowxdxpj3muMOf69x1r7Xmvt3811Dmvt/1prn1WcETupueL0eR6/0hiTMMYcM8ZEjTH3G2MuK+YYvZCaW//MGPOgMWbUGLM/Fez1TDnmN6l/j3NnPPeHqftflvr7J1J/f/+M496fuv8TRXhLUgIKYqTqGWOuBb4IfAZYBiwF3gu8CKhNHbYd+KMpzwkBbwGemnG6p4C3pB5Pe1fq+SIiUjyvs9Y2A6uBfwQ+AnwjkyfO+AwvN3dZa5uAVtz7+U9jzKKZB5X5e/gScA1wLbAY6AL+GnjNjONmzr2LgRcAB+c7LkVzb4VTECNVzRjTAnwSuNpa+31r7bB1fm+tvcJaG0sd+hPgxVMmitcADwL7Z5xyP/AQ8OrU+duAFwI3F/q9iIjIyay1R621NwNvBd5ljOkGSGXZ/7/Un19mjNljjPmIMWY/cH36vvR5jDErjTE3GWMOGmMOG2O+PPV1jDGfNcYcMcbsMMa8dq7xGGPebYx5NHXsL4wxq1P335465IFUpuWtC7yvJPBNoB44LZWR+L4x5gZjTBS40hjTaYy52RgzaIx50hjzJ1PGETTGfNQY81QqY7XVGLMy9diZxphfpZ73uDHmLVOet8EY80jqOXuNMR9K3b/EGLMllfkaNMb879TM15TnrwWuBnqstb+y1o5ZaxPW2justVfOOPxG4K3GmGDq728DfghMzDjuHqDBGPPs1Gs8G6hL3S8VSkGMVLsXAGHgxwscN546Jp3q/iPg23Mc+21OXBHqST0vNsexIiJSBNbau4E9wEvmOGQZ0IbL3Gyc+kDqS/QWYBewBpc56JtyyPOBx4ElwD8B3zDGmJkvkFqW/FHgTUA78L/Ad1Pjuyh12LnW2iZr7ffmez+pTMt7gGPAE6m7Lwe+j8vS3Jga4x6gE/hD4O+NMa9IHftBXFCwAYgA7wZGjTGNwK+A/wA6cPPYV40xZ6ee9w1gUyrL1Q38T+r+a1Ov1Y5b0fBRwM4y9FcAu6219873/lL6gUeAV6X+Pt/c+x1OzL3vSv1dKpiCGKl2S4BD1tp4+o7U2umh1Hrqi6Yc+23gj4wxrcBLgR/Ncc4fAi9LZXnm+8AVEZHi6scFKrNJAv/PWhuz1o7NeOwCXCDwl9baEWvtuLV2ajH/Lmvtv1lrE8C3gOW4L/IzvRf4B2vto6l55++B89LZmAxdaIwZwmX+3wa80Vp7NPXYXdbaH6WyNEtwy6I/khrv/cDXOfFF/z3AX1trH0+tQHjAWnsYuAzYaa293lobt9b+HvgB8ObU8yaBs40xEWvtEWvtfVPuXw6sttZOpmqKZgtiljBjFUMqCzZkjBmf5d8iPfeeCbRaa++a49/lBuBtxpgaXOB1w5z/glIRFMRItTsMLJm6dtha+0JrbWvqsalFoHfgrjB9DNgyyySXPm4M+Clufe9ia+3/FW74IiKShS5gcI7HDlprx+d4bCUuUInP8fjxL+XW2tHUH5tmOW418MXUF/ah1FhMalyZ+q21ttVau8Rae6G19tYpj+2e8udOYNBaOzzlvl1TXmslJ9d1psf4/PQYU+O8ApepAvgDXPZmlzHmNmPMC1L3fwZ4EvilMeZpY8xfzTH+w7hg5zhr7QpccBPG/XtMdRMue/NnzJNdsdY+k3r9vweesNbunutYqQwKYqTa3YVb6nX5Qgem3IBLmS+UXfl26jhdCRIRKQPGmOfhvsDP1Q55tqxB2m5glQfF8rtxS7Fap/zUW2vvzPO8aVPfQz/QZoxpnnLfKmDvlLGcNscYb5sxxiZr7fsArLX3WGsvxy01+xHwn6n7h62111prTwVeD3zQGHPxLOf/H2CFMWZ9Rm/IBYU/A97HwkvE0nOvVkBUAQUxUtWstUPA3+LW+/6hMabZGBMwxpwHNM7ylC8BrwRun+WxqW5LHfcvHg5XRESyZIyJGNeGuA+4wVr7UA6nuRvYB/yjMabRGFNnjHlRDuf5GnDdlAL0FmPMm6c8fgA4NYfzniSVibgT+IfUeJ8D/DEnLq59Hfg7Y8wZxnlOqvvXFmCtMeadxpia1M/zjDFnGWNqjTFXGGNarLWTQBS3DA9jzGXGmNNTtUBHgUT6sRnjehzYDPQZY15pjKlP1Ry9cJ6381HgpdbanQu87e/h6mf+M5N/I/G3cm6/J1IU1tp/MsbsBT6Mu3ozAjyNa8d5J3DllGMHgf/O4Jw2k+NERKRgfmKMieO+SD8CfA4XRGTNWpswxrwOdyHrGVzG4z+ArJYLW2t/aIxpwn2BX437sv8r4L9Sh3wC+JYxph7YaK3N98v423DvuR84gqv5SS8/+xxu+dYvcUu5HsPV1xw2xrwq9fjncBe8H8A1AgB4J/DlVODxOG6pGcAZwJdxy66PAF+11v56jnH9KfDnqfOfDgzh2iG/FffvO421tj/1HuaVWs5960LHSWUws9dciYiIiIiIlCctJxMREREREV/xJIgxxrSmNlh6zLgNnF6w8LNERESKQ/OUiEhl8aom5ovAz621f2iMqQUaPDqviIiIFzRPiYhUkLxrYlIb+t0PnDrHpkYiIiIlo3lKRKTyeJGJOQU4CFxvjDkX2Aq831o7MvUgY8xGYCNAY2PjujPPPNODlxYRkVxt3br1kLW2vdTjKALNUyIiPjTfPOVFJmY98FvgRdba3xljvghErbUfn+s569evt/fee29erysiIvkxxmy11ma04ZyfaZ4SEfGn+eYpLwr79wB7rLW/S/39+8BzPTiviIiIFzRPiYhUmLyDGGvtfmC3MeZZqbsuxm0qJSIiUnKap0REKo9X3cn+HLgx1fHlaeAqj84rIiLiBc1TIiIVxJMgxlp7P1Dx66pFRMSfNE+JiFQWTza7FBERERERKRYFMSIiIiIi4isKYkRERERExFcUxIiIiIiIiK8oiBEREREREV9RECMiIiIiIr6iIEZERERERHxFQYyIiIiIiPiKghgREREREfEVBTEiIiIiIuIrCmJERERERMRXFMRI1qKxKL1be4nGoqUeioiIyEmiUejtdbciUpkUxEjW+rb1sWnLJvq29ZV6KCIiIifp64NNm9ytiFSmUKkHIP7T090z7VZERKSc9PRMvxWRyqMgRrIWCUfYuG5jqYchIiIyq0gENmqaEqloWk4mIiIiIiK+oiBGpEqo0FVERMqZ5inJhoIYkSqhQlcRESlnmqckG6qJEakSKnQVEZFypnlKsqEgRqRKqNBVRETKmeYpyYaWk4mIiIiIiK8oiBEREREREV9RECMiIiIiIr6iIEZERERERHxFQYyIiIiIiPiKghgREREREfEVBTEiIiIiIuIrCmJERERERMRXFMSIiIiIiIivKIgRERERERFfURAjIiIiIiK+oiBGRERERER8RUGMiIiIiIj4ioIYERERERHxFQUxIj4WjUJvr7sVEREpR5qrpBAUxIj4WF8fbNrkbkVERMqR5iophFCpByBSTNGo+xDt6YFIpNSjyV9Pz/RbERHxt0qbp0BzlRSGMjFSVSrtalAkAhs3Vs5EJyJS7SptngLNVVIYysRIVdHVIBERKWeap0QyoyBGqkr6apCIiEg50jwlkhktJxMREREREV9RECMiIiIiIr6iIEakjFkLk5MwMVHqkYiIiIiUD9XEiJSBRAJ274annoJHHoGdO2F0FGIxSCbBGGhpgbVrYf16OP98qKkp9ahFRERESkNBjEgJHTwIt94Kt9/uAhZrobERmppg8WIIBNyPte7xbdvgt7+F006D974Xli4t9TsQERERKT4FMSJFZq3LuvzsZ3DXXRAKQUcH1NbO/RxjoK7O/XR0wL598LnPwd/93fzPExEREalECmJEimhsDP7zP+F//gfCYVi5EoLB7M+zbBns2gU7dsCznuX9OEVERETKmYIYkSIZGIAvfAH274dVq3ILXqYyxv2IiIiIVBsFMSJFcPgwfOpTrsvYqlX5n+/YMVc748W5RERERPxGLZZFCsxa+M533FIyLwrxEwmX1bnySlcjIyIiIlJtFMSIFNjTT8P998Py5fmfKxZztTCvfS0897n5n09ERETEj7ScTKTA9u51t/nUr1gLhw65bM6VV8LLX656GBEREaleCmJECiwScUGItdkHHsmk20tmdNTtDXPllaqDEREREVEQI1Jg3d2wdi088URmLZWthZEROHLE1b+cey5cdhmcfrqyLyIiIiKgIEak4EIhuPZauOkmuPVWF6QEg1Bf7zItiYT7icddkGKtawBw6aVwwQWwYoX3Y4pGoa8PenpcpkhERKScaJ6ShSiIESmCujp4+9thwwZ46ilXnL9nD9TWumCmsRGam6Gry/20thY269LXB5s2uT9v3Fi41xEREcmF5ilZiIIYkSJqbYV169xPKfX0TL8VEREpJ5qnZCEKYkSqUCSiK1siIlK+NE/JQrRPjIiIiIiI+IqCGBERERER8RUFMSIiIiIi4iueBTHGmKAx5vfGmC1enVNERMQrmqdERCqHl5mY9wOPeng+ERERL2meEhGpEJ4EMcaYFcClwNe9OJ+IiIiXNE+JiFQWrzIxXwA+DCQ9Op+IiIiXvoDmKRGRipF3EGOMuQwYsNZuXeC4jcaYe40x9x48eDDflxWpGNEo9Pa6WxHxnuYpkfxonpJy5EUm5kXA640xO4E+4BXGmBtmHmSt7bXWrrfWrm9vb/fgZUUqQ18fbNrkbkWkIDRPieRB85SUo7yDGGvtddbaFdbaNUAP8D/W2nfkPTKRMlHoK1A9PbB5s7sVEe9pnpJKp3lKqpH2iRFZQKGvQEUisHGjuxUREcmW5impRiEvT2at/Q3wGy/PKVJq6StPugIl4n+ap6QSaZ6SauRpECNSidJXoERERMqR5impRlpOJiIiIiIivqIgRkREREREfEVBjEgFUA9/EREpd5qrxEsKYkQqgHr4i4hIudNcJV5SYb9IBVBnGhERKXeaq8RLCmJEKoA604iISLnTXCVe0nIyERERERHxFQUxIiIiIiLiKwpiRERERETEVxTEiIiIiIiIryiIERERERERX1F3MhERERHJWDIJR46AtdDYCHV1YEypRyXVRkGMiIiIiCzo8GH40Y/g7rthctIFLsmkC2TWr4fzzoOzzoJwuNQjlWqgIEZERERE5vXQQ/Av/+KClmXLIDTlG2QsBnfdBbfdBvX18NrXwkteAq2tJRuuVAEFMSIiIiIyp3374EtfgpYWaGo6+fFwGJYvd38eH4cf/hB+/GPYsMEFNPX1xR2vVAcV9ouIiIjInH78Y7d0bLYAZqa6Oli1CpYuhZ/8BK67Du67z9XPiHhJQYyIiIiIzGpoyNXALF2a3fNqamD1aggG4QtfgM9+Fvr7CzFCqVYKYkRERERkVk884bIowWBuz29uhlNOgaeegr/+a/jv/1ZWRryhIEZEREREZrV7d+4BTJoxrhnA0qXwrW/B174GIyPejE+ql4IYEREREZnV00+7FspeCIddVubee+Fv/xaeecab80p1UhAjIiIiIrMaG8s/EzOVMbBypTvvJz/pAhqRXCiIEREREZFZJRIu8PDa4sXQ1gZf/KLrfpZMev8aUtkUxEhGorEovVt7icaipR6KiIjISaLxOL39/UTj8VIPpaJYW5ggBqChwXUw+8EP4F//1WVnRDKlIEYy0retj01bNtG3ra/UQxERETlJ38AAm7Zvp29goNRDqSgtLTAxUbjzh0In6mQ+9zkV/EvmQqUegPhDT3fPtFsREZFy0tPRMe1WvLF0KTz+eGFfwxi3QeaOHfCZz8AHPwiRSGFfU/xPmRjJSCQcYeO6jUTC+lTxRDQKvb3uNp9jyl0lvAcR8YVIKMTGzk4iIV2f9UL64zsSmTsTE4vB1q3uNl/GwIoVsHevC2SGhvI/ZyY0TfmXghiRUujrg02b3G0+x5S7SngPIiJVKP3x/bvfzV0Ts20bbNnibr3S1QUHDsCnPw2Dg96ddy6apvxLlytESqGnZ/ptrseUu5nvIRp1M0VPj9YKiIiUsfTH9vOf7zqIzaa7e/qtVzo7Yf9++Id/gA9/GNrbvT3/VJqm/EtBjEgpRCKwcWP+x2Rrxw74t39zl9UaG91rNDW5P8/8aWhwO5PV1rrbcNhVYM5ltk/+me8hfckLvH9vIiLimfTH9/Cw61A2W5eycBjWrSvM6y9b5jIy//AP8JGPuNocL8ycqjRN+ZeCGJFKkcnlo3AYWlvh6FH3MzjomvOPj8P27bB2rQta4MRsZYz7CQTcTzAI7343nH769HNn8slfCdklEZEq0tzsalWGh/PPTMRibulZd7ebjhaydCkcPAh///dw9dVw2235Z0gWmqo0TfmHghiRSpFJENHZ6XLzM/X2uib9b3nLiedaC5OT7icedz/j467/5cqVJ56bDp42bIDNm+f/5C9EdklERApq/Xq4+eb8g5h0DQ1knsFpb4dDh9z0dttt7r5sp5Gp1/gWClI0TfmHghiRSpHP5aPZnmuMy8qkMzNzSQdPmzfrk19EpAKdeSb8+Mf5nyfXGpolS9xCgVAIXvWq7F935jU+TVWVQd3JRCrF1MtH2faLTD83l8tsPT0LZ2BERMS3Vq9217USifzOM7WGJtvWzJ2dbhXzv/0bHDuW3etqmqpMCmJEKk2x+0XmEgCpMb+IiG/U1cEZZ3j3kZ1ra+Zly1zXsi9/ObsASNNUZVIQI1Jp/HDJSY35RUR8Zd06777Qd3fDZZfl1pp5xQp47DH45jfzzwzNR9NU+VNNjEwTjUXp29ZHT3cPkbAapPuSH6oS1f5FRHIUjcfpGxigp6ODyHxt38VTZ5wx96aX2cqnNbMxbnnbXXe56e7tb/duXFNpmip/ysTINH3b+ti0ZRN923TpoZpF43F6+/uJxuOFeYF8anBEpKr1DQywaft2+gYGSj2UqrJihevzMjlZ6pG4oGXZMviXf4Hvf78wr6FpqvzpEoZM09PdM+1WqlP6SwLAxs7OEo9GROSEno6OabdSHKEQvPCFcMcdrsi+1B59FH73O/jnf4auLjc2qS4KYmSaSDjCxnVlvhRJCk5fEkSkXEVCIV1cKZEXvAB+/etSj8JJ19OcfrorwI9EcquxEf/ScjI5ycjECNbaUg9DSij9JUHrzUVEJO2006ClBUZHSz2SE3U1LS1uH5kvfhF27Sr1qKSYFMTINIlkgr+7/e/Yd2xfqYciIiIiZSQQgFe8Ag4dKvVIpmtqgvp6+MIX4OjRUo9GikVBjEyzd3gvR8ePUheqK/VQZAEFL74XERGZYf16sNb9LCQWy35Ty1y1tblNMDdvBk2L1UFBjEyza2gX0ViUgNGvRrlThx4RESm25cvh1FNhaGjhY3Pd1DJXnZ3w8MPwwx8W5/WktPRNVabZMbSDgAnQWNM47f5oLErv1l6iMW1dWy56OjrYvHatiu9FRFB2uliMgVe9KrNlW/lsapkLY2DlSrj5ZrjnnuK8ppSOghiZZvvh7YRDYWqCNdPu1/4x5UfF9yIiJyg7XTznnusK6ycm5j8uXXwfDhdnXOBaQS9f7jqW7dlTvNeV4lMQI8dZa9l/bD9t9W3H70tnYDacsYHNl23W/jEiIlJW0hmYDW1tyk4XSV0dvOxlUK7xYkODC5y+9CVXJyOVSUGMHDc6OcrY5BirWlYdvy+dgbnliVvYuG4jkbC2rhURkfKRzsDcMjio7HQRvfjFroC+XHdkWLIEBgfhm9+ERKLUo5FC0P/pctyh0UOMTo7S0XjiKlY686IMjIiIlCNtzlsaK1bAWWfBM89Ae3upRzO7ri6491645RZ43etKPRrxmjIxArhlY/92379xbOIYHQ0nJoJIOKIMjMxKRbQiUmyzfe6oPrA0jIHXv94t1yrXbMzEhNvT5j/+Ax58sNSjEa8piBHALRv71P9+iv7hflrqWko9HPEBFdGKSLHpc6e8nHkmrF6dWbvlUti2DX72Mzh8GL7yFdi/v9QjEi/psoUAbrnYbTtv48j4EQUxkhEt4RCRYtPnTnkxBt70Jvj852HRolKP5mTp1s7d3a4l9Fe/Cn/911BbW9pxiTeUiRHALRtb07qGptomLR3zsWIu8dISDhEpNn3ulJ9zzoGlSyGa4TZysRhs3epuC21qi+eODle/86MfFf51pTgUxAgASZvkwMgB6kJ1NNc2l3o4kiMttRARkWIKBuENb3CdwDKxbRts2eJui23FCvfajzxS/NcW7ymIEQCGxocYj4/TVNtEOFTEXanEUz0dHVnvk6ACfRGpBJPJJI+MjPDoyAhj6qlbVOvXQ0sLjIwsfGx3N1x22YmlXpnwKnsTCrnWy//6r255mfibghgBXHvlWDzG4obFpR6K5CGXpRbzZW8U4IiIHzwxOsrf7NjBZ3bv5jO7d/M3O3awe3y81MOqGrW1cPnlcPDgwsdOXeKVqYWyN9kEOZEIjI3B9ddDMpn5GKT8KIgRAAbHBhmPj9PeUKbN3qVg5sveaHmaiJS7nWNj/OMzzzCaSLA6HGZVOMx4Msn/t2sXexTIFM0LXgD19VCIf/KFsjfZLlHr6nJBz+23ezdGKT4FMRUoGovSu7WXaCzDKjugf7ifhE2wrGlZAUcm5Wi+7E0uy9NERBbiVZZ3OB7nS3v30hgIsKim5vj9i2tqCBjDtw8cIFmum5hUmPp6uPRSOHDA+3MvlL3JdomaMS6QueEG2LfPu3FKcSmIqUB92/rYtGUTfdv6Mn7O7qO7CQVCLG1cWsCRid+oE5CIFIJXWd4fHDzI0XictikBTFpHKMRjo6P8LtO2WZK3l77U1Z1MTBT3dXNZolZXBzU10NsLWjHtT/pmUoF6unum3Waif7ifhpoGFtWXYaN3ERGpKF7s93JwYoLbh4bomuObqzGGJaEQPzl8mAsjEYwxOb+WZKa52dXG/OAHbhPMctfRAU8/DT//ucvkiL8oE1OBIuEIG9dtzHi/l0QywcHRg9SH6mkJa6NLEREpLC+yvL8YHMQYQ3Ce4KQ5GKQ/FmN3MTYlEQAuvhgaGzPrVFZqxri2y9//PuzcWerRSLbyDmKMMSuNMb82xjxijHnYGPN+LwYmxTM0PgS4q1YtdQpiRKSyaJ6qPEfjcX49NMTyBbZeTwc5v9WSsqKpr4e3vjWzTmXloKbGZZA2by7OBpziHS8yMXHgWmvt2cCFwJ8aY8724LxSJIfHDmNThY/a6FJEKpDmqQrzwLFjJK0llMESsY6aGm4bGlKBfxG94AWwbBkMDZV6JJlZvNgV+P/4x6UeiWQj7yDGWrvPWntf6s/DwKNAV77nleJJt1deVLeIYCBY6uGIiHhK81Tl+fWRI7RkuBQtHAgwlkiwv9jV5lUsFIIrroDBQfBL7LhiBdxyC+zaVeqRSKY8rYkxxqwBzgd+5+V5pbD2De9jMjlJR6Pa6FaLpLVsjUb5jwMH2HLoEIc0uUuV0Dzlf4OTk+wYH6clmPlFN2sMO8bGCjgqmam7G84+G/yyzVgo5Gp5/v3fIZEo9WgkE54FMcaYJuAHwDXW2pMWnxpjNhpj7jXG3HvQLwslq8Tu6G4CJqD2ymXAq70TFvKfAwN8cc8efnPkCDcdOsRHnn6aXw0OHl9WKFKJNE9VhsdGRjCQVbexOmN4dHS0cIOqItGoa0u8UJmRMdDTA6OjkEwWZ2z5WrIEnnpKm2D6hSdBjDGmBjcx3GitvWm2Y6y1vdba9dba9e3t2hW+nOwd3kvABFjevLzUQ6l6Xu2dMJ8nRkf52eAgq+vq6Eztbr20poZv79/PjQcOKJCRiqR5qnI8ODJCfSC7ry+RYJDHFMR4oq8PNm1ytwtZswZe/GL/bChpDHR2wne/C4cPl3o0shAvupMZ4BvAo9baz+U/JCmmRDLBwZGD1IXqaKtvK/Vwql5PRweb167Na++Ehdx19CjhGW1JawMBTqmr4xeDg/xicLBgry1SCpqnKoe1lm0jIxnXw6TVBQIcnpxkwi8pgTLW0+M6efVkuBXdG97glmdNThZ0WJ6pq3N1PN/9rn/qeaqVF5mYFwHvBF5hjLk/9bPBg/NKEaTbKwdNUO2Vy4AXeyfMJ2EtvxseZvEsu1sHjGFlOMx3BwbYriuWUlk0T1WIg5OTjCYShLPMxBhjCBjDoF++SZexSAQ2bnS3mWhvh0sv9U82Blw25u674cEHSz0SmY8X3cnusNYaa+1zrLXnpX5u8WJwUniDY4MY3BV5bXRZ+fZPTDCWTM75BaA2EKAtFOKre/cS0xVLqRCapyrHrvFxbBa1MDMNFrjeUGb36ldDOAx+6a1gjKuPuf56V9Mj5cnT7mTiP4NjgyRJkrRJWutaSz0cKbCxRIKFpv+WUIiheJw7/dLgX0SqxmOjo8y/veXcEtZyWJmYkmhqchtg7tvnnyVakYhrXrBlS6lHInNREFPl+of7SSaT1NfUEw6FSz0cKbDxDLMrHTU1/ODQIcbUZ1JEysiO8XEas2itPFWtMfRrS/aSeclL4PTTwU+N/zo7tXdMOVMQU+X2Du/FGEN7ozrxVIMkkMlFsIZgkGOJBLfnmI0pVqtoEake1lr6YzEacg1iAgGO6DOpZIJBuOoqtzyrHP4zxGKwdau7nUt675hvfrM8xizTKYipcvuG92EwLGtcVuqhSBEEYMHlZGnLamv54aFDjOSQjSlGq2gRqS7DiQQT1hLKsSYmZAzHlF0uqZUr4XWvg717Sz0S2LbNLRXbtm3+45YsgR074LbbijMuyZyCmCpmrWVgxH3J7GzuLPFopNCi8Tg3HTqUcYvRukCAmLX8Xw7ZmHSr6A1tbcrIiIgnDk1OZnwRZjYKYsrDhg2waBEcPXryY5lkR7zS3Q2XXQZnnDH/a6b3jvmv/4Lh4cKPSzKnIKaKDU8Mk7DuA31p09ISj0YKrW9ggP+3cydPj49n/JzFoRC3Dg1lvQFmulX0LYODysiIiCcOT07mtRlvyJicMsvirfp6ePe73WaSM6+pZZod8UI4DOvWwRNPLPyadXUwMQE/+1nhxyWZK8xmFOILR8aOYDAYY9SZrAr0dHSQsJa7olGS1hLIYElGYyDAM7EYu2MxVtXV5fSaU29FRHLVH4th8mivHDJGWeEy8exnwwteAPfdB11dJ+7v7p5+WwyZvmZnJ/z85/Cyl4GmtPKgTEwVOzJ+BJsq815Ut6jEo5FCi4RCvK+ri5XhMKMZLilLbxB3b4459EJv3iki1WNnLEZDlptcThUExpLJvLI54g1joKcHAoHpe8eksyPhIjZLzfQ1QyHXnOCmm4ozLlmYgpgqdnj0MEmbJGETc2ZiorEovVt7icaixR2cFMwpdXUMTU6ydXg4ow0t20Ih7jp6VBO/iJTU3nmCmFgyueBnmjFua2dt5FseFi2Ct70N9u+ffe+YYtbHZGr5crjrLlfoL6WnIKaK7R3eS8AEaKptmnOPmL5tfWzasom+bX1FHp0Uyil1dTwyNsaWwUG2jYwseHxDIMChyUkOTEwUYXTziEaht9fdikjVGYrHqZ0jiNk2MpLRZ5oxJuP9sqTwXvISV1g/294xxayPyVQg4Fou9/XNHnhpmiouBTFVrH+4H4CljXMX9fd097D5ss30dPcseD7tDeIPy8JhTq2r47K2NrobGxc83hiDBbZPzfmXQl8fbNrkbkWkqiSsZSKZnPNLS3djY0afaRPJJN/cv1/zVJkIBuFd73JLymb+J0l3DytmfUwm2tvhscfgoYdOfkzTVHFpoXoVO3DsAAbD8qblcx4TCUfYuG5jRudL7w0CsLFTLZvL1Sl1ddQGApzf1JRRcT+4dstPjY1xUWtrYQc3n56e6bciUjXGk0m3HGyOz6xwIMC65uYFz/P0+Dg3DAzQFAxqnioT6b1jbr4Z1qw5cX+6VqXcGANtbXDjjXD22a5WJk3TVHEpE1OlJhOTDI0PkbRJFtUv8qTuJb03iDpRlbeGYJCzGhqy2rm6MRBgZxatmQsiEoGNG92tiFSVgxMTbB8dzbue5dS6Or54+umap8rMpZe6LmWHDpV6JJlpaXG1PHfeOf1+TVPFpSCmSg2NDx2/qrV131ZP6l7Uico/LoxEstr0rSEYZG8sRlLF/SJSAt8bGODOaDSjOr751AQCvGf5cs1TZSYcdsuwRkbcfix+sHQp/Od/wuhoqUdSvRTEVKkj40cAMBjefPabM657Ef+aWrOUXjeeaVASNIYkMDg5WcARiojM7tVtbbwwEsmojm8+Fvd5JuVn1SqXkfn1r6HUif9MNDS4oOvWW0s9kuqlIKZKHRlze8RYLCtbVrJx3UYi4cLkP1XwXx7SNUt9AwMsqqnhrMbGrJaUGeBgkYMY/e6ICLi6vDMbGgjnsU8MgLWW4Cz367OmPBw8CL/97cnLtMrV8uXwgx/AF7+ojmSloHxqlRoYGcBgsNbOuUeMV1TwXx7Sa8DTty9raeGrIyMsrqnJ6PlJaxmYmOCsPK+EZkO/OyICkPni17klrSVgzKyZGH3WlIcrrnDZjYcectmYurpSj2h+tbXw5JPw7W9Dfb2rh5HiURBTpfZE9xAwAZrDzYQChf01mPnlWUojXbOUdl5zM62hEMcSCZqCs12bPFmsyDUx+t0REYC4B589CWtpDARm7XCmz5ryEInAhz4Et98OX/86nHKK6wZWzp7/fNeh7PWvL/VIqo+Wk1Wpfcf2YTAsbZp7jxivqOC/PNUGAry5vZ1DGS4RM7gvAcWk3x0RAW8+e+LW0jDHBRt91pSXl7wELrgA9u4t9UgWFom4YOvhh0s9kuqjIKYKWWs5cOwAAMualpV4NJINr9dtXxCJsKSmhmGtAxeRMpawNu/uiJPWZpx1ltx5sWu9MfBHfwSNjXD0qHdjK5SODvjRj/zTWa1SKIipQiOTI0wkJphMTtLV3FWQ11CRZGFMLc73Qk0gwFs7OjiUwX+nJNCoLwAiUgIJa+fc6DJT8RlBjOapwvBq1/pIBN73Pjh8GMr9P1FDAwwPwz33lHok1UV50yp0ZOwIARPAYlncsLggr6EiycIoxLrt5zY301lby5HJSRbNU+RvcB2CRESKrcaDz56ZQYzmqcLwctf6M8+Eyy+Hm2+GNWvKuz5m8WL44Q9P1MhI4ekbSRU6Mu7aKwNzdiaLxqL0bu0lGsstH9zT0cHmtWtVJOmxQqzbDhrDHy9fztFEYsHdsLvCYc9eV0QkU3WBADO/v8aSSbYODy/4uXX8eGtpn3KhRvNUYXi9a/3rXw+nngoHDnhzvkJpbnYtoh98sNQjqR4KYqrQkbEj2NTa4kV1i2Y9pm9bH5u2bKJvW275YBVJ+svpDQ28Y+lS9sZis647H0kkWFxTQ2dtbQlGJyLVbrYgZtvICFsGB9k2MpLROZLW0jnlQozmKX+oqXHLyoxxS7bKWWury8ZkGFdLnvR/bhXqH+4nGAhiraWlrmXWY3q6e6bdSuW7eNEi+icm+NXgIKvC4WnLNw5OTvKGJUvyXpMuIpKL2Zaydqf2rOrOcO+qgDEZ74sl5aWjA/70T+Ezn3F7x5Trf8bWVti5Ex57DM4+u9SjqXzKxFShdHvlJQ1LCJjZfwUi4Qgb120kEvYoHyxlzxjDO5Yu5c3t7eyZmGDfxATHEgn2xGJ01Nbyqra2Ug9RRKpUXSDAzBxxOBBgXXMz4QzrZZLAknL99isL6u6Gt7wFnnkGitztP2PGuGVlP/pR+Y6xkigTU4XUXlnmEjCG17e3093UxK8GB+mfmODMhgbe2N4+5/4KIiKFNlsQk424tdQaQ0SfY7722te6TMd998HKlaUezewWL4bt22H3bli1qtSjqWwKYqpM0iY5NHoIYwxdkcK0Vxb/O7W+nk1d+v0QkfIQNgaL2+csl2Wto4kEneGwlsT6XCAAV10Fe/a4Ivr29lKP6GTGuO5kd9wBb397qUdT2bScrMpEY1EslngyTmezWkqKiEj5CwUC1BhDrvXSo8kkq9RdsSI0NMCf/7nbWDLDng5F19EBt90GY2OlHkllUxBTZY6MHcFgCJrgnJ3JMqWNwkREpFhaQiEmsmz7FEvE2LrvPo7FYwpiKkhnJ1x9NezfX54bYdbWQiymdsuFpiCmygyNDx3fI2ZRfX5BjNe7x4t/KIAVkWLrrK1lLMsgZtvAw2zZ/hN2HNnBMgUxFeX88+GNb5y70D8Wg61b3W0ptLTAz3+uAv9CUk1MlRkcG8Rai8XmnYkpxO7x4g/a6VpEim1lOMwjWa4f6u54NtZCS+spysRUoMsvd4X+jzwCK1ZMf2zbNtiyxf153bqiD43WVtixA/btc5kj8Z4yMVVm7/BeAiZAQ00DdaG6vM6ljcKql693uo5GobfX3YqIb3SFw8SzvKwdDoY5a+m5rKxrpFlzVcUJBuFP/sRlPQ4fnv5Ydzdcdpm7LQVj3PjuuSf752qayoyCmCqzb3gfAEublqpLi+TM1wFsXx9s2uRuRcQ32mpqCOQwbx2Nx3lOU1MBRiTloLkZrrnGLRsbHj5xfzjsMjClTMAtWQL/8z+QSGT3PE1TmfHhNxDJx/6R/Vgsy5uWl3ooIqXR0zP9VkR8ob2mhmQOBQZx4KyGBu8HJGVj5Ur44Afh05927Y3r60s9Iqe+HgYG4Omn4YwzMn+epqnMKBNTReLJuCvst5YVkRULP0GkEkUisHGjuxUR31gUCtEcDBLLorjfWou1ljV1+S2flvJ35pnwvve5GpSJiVKP5oSaGrjzzuyeo2kqMwpiqsjQ+BAGl4pvbyzDHaJERETmYIzhrMZGjmbRFfFoIsEp9fUsqqkp4MikXFxwAVxxBezenf0SrkJpb4e77irPVtB+pyCmiqSDGGNM3p3JpLjU0lhEBJ7T2MhYFkvKjsbjvLilpYAjkrRyKUZ/1atcQf/OnZBlR+6CqKlxmaE9e0o9ksqjIKaKpJeSQf57xEhxaU8eERE4rb4eUkvEFpK0ble081TUXxTlUoxuDLz5zXDRRbBrV3ns02ItPPlkqUdReVTYX0UOjR7CYjEYWsK6MuUn2pNHRASW1taysq6OaDxOywLdEQ9NTtLd2EiblpIVRTkVowcC8K53wdAQPPaYK/wvpUgE7r0XLrmktOOoNMrEVJG9w3sBWFy/mGAgWOLRSDZ83dK4XJTLWgcRyZkxhksWLeLIAktrrbWMJJO8fsmSIo1Myq0YvbYWrr4aurpcsX8pRSLwxBOuDfR8NE1lR0FMFemP9gNuj5hsqSZDfK+Uax00M4l45rymJoLGnLTxZSwRY+u++4glYhyOx3lWfT1nlEuvXSmJxkb4wAegqQkOHizdOIJBt6Rs5875j9M0lR0FMVVkYGQAi6WruSvr56omQ3yvpwc2by7NWodyWSwuUgEioRAXt7bSP+Oy9raBh9my/Sc8eOBhhhMJ3tTerk2dhUWL4EMfcvvHlDKQMcYtbZuPpqnsKIjxsWgsSu/WXqKxhcPmicQEwxPDGAydzZ1Zv1ZPRweb165VTYb4VynXOpRyZhIpoWzmqWxc3t5OYyjEsSl9dLs7ns2lZ7yOlpYzeE1bG8/SBpeSsnw5fPSjEA7DgQOlGUNLC2zdOv8xmqayoyDGx/q29bFpyyb6ti0cNkdjUYImSMAEWNKQ/Rph1WSI5KHcFouLFEk281Q2GoNB3r1sGYcmJ4/vGxMM1LKk7Wye3RThD5WFkRmWLoXrrnNLy/bvL/7rNze7/WuGh4v/2pnw4zSlb6Q+1tPdM+12PkPjQ8f/3FbfVqghiRTX8DDs3QsNDdDa6q9PX5EqkM08la3zm5v56OrVfHnPHnbHYlhgw+LFXLZ4MTUBXaOVk7W3u0DmM5+B/n7ozH5hSs6McV3TDh50AY3kT/+X+1gkHGHjuo1EwrN/cZuaxj86fhSLJWETWQcxCxX1q+jfX4r936sgrzc5Cd/8Jrz//fBP/wR/+7fuz5/9LNx/v3tcREoum3kqF2sbGvjM6afzoc4lnHPsDl7TUkdD8OTum4Va1iaFUcgi87Y2+MhHoKPDXQMD1zVs69aFu4fly1o4erSwr1FNFMRUsKlp/CPjRxifHKcl3EI4FM7uPAsU9V+/bx+btm/n+lL3MJSMFLtJQ0Fe7yc/gd/8xvXOXLXKbQKwciXs2AFf+IK71Pbgg+Wxy5mIzMmL5WbhQID/feImPvjTP5nzPNf//no2bdnE9b+/PufXkeIpdJF5ayt8+MNuCtm9Gx56CLZsgW3bCvN6ackkHDlS2NeoJlpOVsGmpvFvfvxmkiRzLuoH2NDWRm9/Pz0dHdNrY9Lrjou0/jgadR9sPT1aPZSLYm+c6fnrxePws5+5oGXqFddAwK0VaG93l7o++1m44AK44grXnkZEyo5Xy83Sz99wxgZ6t/bS090zPftjZtwWWDQWpW9b38njkIwUY+PMSMR1LfvCF1wG5tJLobu7cK8Hbu+aUtTjVCplYipYJByhp7uHvm197Dyyk6RNsiKyIvvzpIr6f3r48KxX1K9atozNa9dy1bJlXg19Xn5sA1hOit2kwfPXO3gQEgmYbxfulhY45RS3tOxv/sbtMiYiZSe93AzIa7lX+jy3PHHLrJmdq867is2Xbeaq867Ke8yZKFRDg2pRrCLzpib44AfhnHNgyZL5pxUv1NWVrjtaJVIQU+HSH6T/u/t/McbkFMSkrWtu5pWtrSddUS/2l2I/tgGsJCWvkRoczCzrZwysWOE2B/jUp+C//1vLy0TKlFdf+nu6e9h82eaTMjsL1eZ4ba5xSPHMV1cz9bGGBrjmGnjRi9yK5EKWVNbUwNBQ4c5fbbScrML1dPeQtEn+75n/IxgIsrhhcc7n+r9olGc1NJS8zXL6Cs1stNSs8NI1LgAbZ2ntstDjeaury+74lhb3nG99C3btgj/6IxfYiEjZ8GpZ2dTMTinNNw4tNSuO9KoNOPk7w8zHwmH44z923cq+9z3XjrmxsTDjUuM872gmr3CRcIS3n/N27t57N+DaK+fyAWqt5ZGREc5rairkcPM234eWeGOhGpeC19zksoFdOAxr1sBtt8HEhJutCr1uQEQyNvNLfyV/0U9nnYCyCLgq1Xx1NbM9Fgi4upgVK+ArX4HRUVdi6SVrp5dySn4UDxZIObVzPDp+FIMhaZO01bdNS9svNM700qAnxsYYSSS4oMzTG1pqVngLLR8s+PLClhbX4iXbpWGBgAtk7rrLrSNQG2apcuU0T800c3nZgnNVGb+XmbTUrDjmq6uZ77Fzz4VPfMJNNbt2uenGK/F44TI81UiZmAIppystR2NHicVjLG5YTF2oblrafqFxppcGfXDFChLWsiKcXXvmYptvqZlUiIYG11Z5eDj7NYPGuEDm7rtdc4D3vte1ixGpQuU0T800c3nZgnNVGb+XmcplyZvMrbMTPv5xuPFGl8Dv6sp+JfNsRkbg1FPzP484CmJytFCqu5C7FGfr6PhRRiZH6G52vQOnfoAuNM70kqBnxsdpCgbpKMASnEqrY6m091OWnvtc+PGPc/sHTgcy993n1gxcfbVbbiZSYfw0T80084v+gnNVgd9LpS1vq7T3Uwj19W7l8dq18O1vu2T+8uX57SZhrdsdQLyh5WQ5WqiTSrE7ocznwMgBJpOTs3YmW2ickVCIdy9bxt5YjOc0NWEKsBdMpbVMrrT3UzLztZZ51rPy6zRmDKxe7TbE/OpXXY5fpML4aZ5ayIJzVYHfS6W1TK6091MoxsBFF8Hf/z08+9mue9nRoycej8Vg61Z3u5B43J3v9NMLN95qU/WZmFyvRpTzFayZ9h/bj8Gwojm39sp7YjGG4nH2xmJE43HPax2KsalVMVXa+ymZ+bo0nHaay+3HYrlnUYxxy9Luvx9+8AN4y1uKtmGrSDaqYZ7yQiGzC5X2b1lp76fQOjrgL/4CHn4Y/v3fXTDT1QXbtsGWLe6YdevmP8f+/a6Ns1ZoeKfqMzG5Xo2YedXHi6LCQhUmDowMEA6Fc26v/NT4OHtiMT69e/dJG116oVibWk0130X+fJXi/VSk+bo01NTAS1/qNr7MRzqQ2bIF7rwzv3OJFEg1zFNeKGR2oRRZq0L+W/spC1cujIHubrft2Jvf7KafSARe9Sp3/3yiUdfZ/7LLijPWalH1mRivrkZ4UVRYqMLEg6MHqQvVZRzERONx+gYG6OnoIBIKcfvQEN2Njbxn+fLCtc0tMrVi9oGFujRcdBH8/OduWVk+GZRg0F1S+/rX3eYAyvVLmamGeSpbs2VdKi27UC7/1jJdOOyCkZe+1DW73LLFZVlCIVi0yNXSGOOmpokJOHDA9Y/5wAdcRke8U/VBjFddQrz48CzEB/BEYoKRiRFCgRCL6zMLYqZuVvie5cvZHYvxopYW3lOIjQtLREu+KkBnJ5x1FuzeDUuW5HeuujpobYUvfhE++Uk3E4mUiUqfp3Ix2xf8Suv6VS7/1jK75maXhXn5y+GJJ1yJ5T33wDPPnAhimpvh4otd0NPSUuoRVx5j8ymOzdH69evtvffeW/TXrUaHRg/x57f8Ocual/H5V39+weOj8TjX79sHxnDVsmXsn5jgXY8+yudOP50X6P9AKTcPPgif/7wr0vdCf7/LxFx7bVXsSGaM2WqtXV/qcZQjzVPlKxqLcv3vrwcDV513FZFwRN22pGzE466VciAATU0qtczXfPOUJzUxxpjXGGMeN8Y8aYz5Ky/OKd6IxqKMTo7S1dyV0fHX79vHNU89BdYSCYX45r59/HZ4mDuGhujt7yeqLk5STs4+G9ravCtuWr7cVWr+5CfenE/KhuapynH976/nml9cA5bjAct3H/oum7Zs4uqfXl2W9TpSPUIhl3VpblYAU2h5LyczxgSBrwCvBPYA9xhjbrbWPpLvuSV/0ViUWCLGqpZVmT0h/X9c6rY1FOKS1lZqg8HjS8w2erisLL2nyoYNcMst2lulKuWzsU4oBFdcAV/4gjczhjGuif9NN7k2zmedld/5pCxonvK/sckx+of7OTBygIcGHgLgV0//iv3H9nNs4hhHxo9wauup3PjQjeyJ7uHcpedSE6yhNljL0salnNJ6CitaVrCqZRWL6xdntV1AOsuz4YwN3PLELcr2VCHt/1aevKiJuQB40lr7NIAxpg+4HNDk4JFc0+TRWJQbHriBWCLGykhmuytdtWyZ+4O1RONx/qyri02dnYwkEtwTjbKhrS2XtzCndIH9FVe4nXFBhfZVJ98uC+efD+ecA089Benf33zU1Lgam69+1W0O0Nyc/zml1DRPFUEuc9Vsz7HWsnd4L7uP7ubxw4/z+KHH3VYBxmCtJZFM8KKVL2I4Nszo5CjN4WYW1S+ivaGdW5++lecsfQ4NNQ1YLJOJSZ4cfJIHDzyIxS2f72js4M1nv5nzl59PwCy8ICVdf3PFOVdw40Nuoqqk2htZmJoBlScvgpguYPeUv+8Bnj/zIGPMRmAjwKpVGWYFfCSTD+/f7vktzxx9hucsfQ6nLTqNmmBNRufOtUNJ37Y+/vm3/8w5Hedk3JksEgpRHwiwaft26oPB41mX/zp4kBsHBriotdXTTEy6sH7DBnje82BszF3x0JWOKpJvlwVj4O1vh7/+a7cY2Yt9jCIR1zDgJz9x5xa/0zxF4Xdpz2Wumvqcy591OXfvvZvbdt3G4OggFks4FKa5tplVLaumZU+OjB9hyxNbiIQjrOt0G3Q8Ofgk2w5uY03rmuP3ATTWNrKYE3Pg0fGjfOnuL9Hd3s01F16z4FycLqzfcMYGntf5PMbiY0RjUWVjqoiaAZWnonUns9b2Ar3gCiaL9brFksmHd3+0n75tfdzw4A10NXfxytNeyXnLzuPWp2+dd1LJtUNJT3cPv9n5G8bj49QGa+nd2pvR5JVuozy1nfJs93lhahfd+np3paO+Xlc6qspCrZQz0dUFr32ta7ns1ZfPzk745S/hJS9xS8yk4mmeml2mwU8uc1VPdw9Hxo5waPQQH/zFB9kxtIPndT6PVa3z/3/c3dE97Xau+2bTUtdCJBzhgQMPsHXfVi5cceG8x0/telZfU8+mLZuoD9UrG1NFvJimxHteBDF7gakz/IrUfVUlkw/vN539Js7uOJvvPPAdnhx8ku888B0+e+dnuWvPXRwZO8JHXvyRWZ+Xa9vISDjCaYtOY3d0Nz9/8ue876fvAxaevCKh0EnZltnum00+60Znu9IRi7kmUdlcXNfa1Sp16aVw++2uLUxjY/7nCwZdRH3DDfCRj7hWM+JXmqfI/YJYpsFPtnOVtZb79t3H44cfpzZYy9Hxo9y15y4W1y9mXf3825+HQ+Fp2Za57ptNLB7joQMP0RRuIhwMZzxe8K7tsbqpieTPiyDmHuAMY8wpuEmhB6i69ReZfnifueRMPvnyT/LQwEPc/PjNPHrwUYZjw/xu7+/4jwf/gzee9Ubqa+o9G9f+Y/tpq2/j7ee8nYAJFLzffD7rRiMR+OM/hiefdKt47rkHDh92G0q95z3FGQMoCPKtxka39GvzZlizxpu2MB0d8OijsHWrW+8ofqV5itwviBVqv5K7997NN+77Bl2RLupCdSyqW4QxZsFMSr7u338/P3/q57z12W/l3GXnZvVcr/aiyXcjSwVBIh4EMdbauDHmz4BfAEHgm9bah/MeWQULBoKct+w8zl16Ls8cfYY7nrmD23fdzo8f/zENtQ284cw3ePI6k4lJYokY57edn9cHbzQep29ggJ6ODiILpERyXTd64IALWm69FY4edZmXYNDtcnvh/Jl+z8aQpgI+H3v+890v0YED0N6e//mMcYHMd74D3d0uMyO+o3kqP4XaRPLWp2+lrb6NulAdkHkmZTaxeIxtA9vo7ugmHJo9uxKLxzgwcoAlDUv4iwv+gr99+d9mVNhfCPkGhvkGQSKVwJOaGGvtLcAtXpyrmhhjWN26mtWtq+np7mH/sf2eXlEZnhgmFo+xpmVNXufpGxg43l65p6Nj3oAmm3WjR4/Cl7/sApVHHnHfF5csgRUrYO9eaGiAv/xLOO207Mab79pVFfD5WDAI73wnfOITsGiRN0X+TU2waxf87GdwySVK0/mU5qnyc87Sc/ivh/+LxtpGaoO1eZ1r28A2tjyxBXA1MemABmBwbJCJxAR1oTouW3sZL171YjoaM6vvLFTGI9/AsFDZMfG/alpNUrTCfplfMBCkK5LZhpSZGo4NY4xhadPSBY+dL9sytah/akCTa5eyZNJlXTZtggcecHXTL3uZKzk4fNhdRH/FK+AP/sB9fyy2Si/gyyaz5kunnOLqY265BVav9mZZWWcnbNni2jh/+MPuvkr+JREpgkvPuBRrLVu2byGejBMJR4iEIwQDwTmfM1fGpbujm0QyQVdzF7ftuo279tzF4bHDrO9czwVdF3DhigtZu3htxl1BwQUwV//06rJsq1yo7Fi5qKYv4l6rptUkFfgNxl8Kua41YRM01zazsmXhzkrzBSdTi/rz6VJmrcu4/Nd/wS9+4QKY7m540YtgfNwFL6tWwTXXZJ99kcx5EYiWvTe+EbZvhz17vNs7JhRyLZy/9jWl6aSqFGqeCgaCXH7m5bxszcu4a/ddPHDgAZ4cfJKETWBwe8Kk93ZJ3z5x+Anu3HMng2ODnLnkTLd3DJZkMsnqRas5o+0MXnnaK3lo4CHeff67WRlZmdXGllP1bevjxodu5IpzrlDGo8iq6Yu416ppNYmCmBKbb11rvhPHkbEjtNS10N4wvTZgtivx8wUnM4/P9ouvtfDYYy54eeopaGlxxfpLlsDpp8O+fS7jctVV8OIXe7MCSOZWqHbZZaWmBt77Xvj4x+HYMW9SesuXw+OPwyc/qUuDUlUKOU+Ba3n8mjNew2vOeA2JZIIj40fYObSTmx69iUtOvYRwMEzCJphMTDI0PsQFz1zAZWsvY3HDYhprGo+/7o8e+xFvO+dtRMIR/uDsP8jvTTN9yZaK54urmr6Ie63SV5NMpa+LJTbfutZ8C/dWt67m0jMuPekq1GxX4ucLTvK5cn/gANx4o8u6RCInmkaNjcHixTAxAe94hwte6uqyfIOSk1wCUV9assQFMv/8z+6XK9/o2Bh3jl//Gt71Lm/GKOIDhZynZgoGgixpWMJNj97EP9zxD6xpXXPSeWcLUHq39vLen74XY4xny6wqfclWOaumL+KSOwUxZSzbwr2ZV8SWNCzh5ae8/KTHN5z5h2xeuzbjK/Ezr9xnslY1FnP7BP7oR+57Xzp4OXYMDh1yF8bf9ja46CI1fJICOvdceMMb3C/iKafkXx+zbJnbi+ZNb4LmZi9GKOJr+c5Tcz2+4YwNbL5sc8bnnTkOtSAWqXzava3E0lex+rb1nfRY+ipQJBwhGovSu7WXaCya07mmPn7LY99nY2dnxkXd6Sv36ePTa1X7Zn8ZduyA//f/4Ac/gKVL3fe+I0dg5063tOyP/9hdHH/1q08EMNEo9Pa6WxFPXX652+Nl9+78zxUKuc4Ud9+d/7lEfKIk89QTtxw/byamjiOT18lFJu9PRIpHmZgSy/QqViYp+4XO5VVLxvnWqt5zD/zrv7pMS3s77N/vApczznArcM49d/ZVPSrik4IJBl3k/KlPufWNSxfu1jevxYtdZ4pXvMKbzmciZc6X81QBWhBrbxaR8mKstUV/0fXr19t777236K87n3g0zkDfAB09HYQi5RfbFSo1nkm73dmOmW1J2Z498KEPuT/X17ttOi65BNavd/sFzjsOtVOUQhsYcPvHhMOuu0Q+du1y6cY1a7wYWckYY7Zaa9eXehzlqBznqXJfIlWweSqD8852jNfjKfd/f5FKNN88VX7f1ktkoG+A7Ztc8XrnxsIVPef6IVioAsNMivZnO2a2zMmSJfCe97jvh8uXQ1dX5heqVcQnBdfRAddeC5/+tEsPtrbmd77HHvN9ECP+UqxMQNnNUxm879mO8frfS4X+IuVFQUxKR0/HtNtCKbd0dCbtdmc7ZrYlZRMT8PTTyqZIGTvtNLjuOvjHf3SBzKJFuZ2npQV+/3t4zWu8HZ/IPIq1S3vZzVMZvO/ZjpntPmVTRCqHlpMVWSV/gPb2uuzM5s3KqkiZ27XLZWRCIWhry/75ySTs3Qtf/apbnuZTWk42N81TFTpPbe1l05ZNbL5sc1kEaCIyv/nmKXUnK7KZHVQqSU+PC2Dm2pxKHcikbKxe7TIyiQQcPpz98wOpj84DB7wdl0gZqOh5qrtn3tbN6kAm4h8KYnyunD5w03Utcy0lW6g1s0hRrVwJH/2o+3OuwcjIiHfjEalQZTVPLRCgFaI1s4gUhmpifK7c1i7PZ77WzCIl0dUFH/84/Mu/wDPPwIoVJ7IsC7EWRkcLOz6RCuCreapIdUcikj8FMT5XjA9cr9ofqwOZlKWODvjYx+A734Hbb3eBTKZ1LsFgYccmUgGKMk95VMejDmQi/qHlZD5XjLXLWgYmFa+uzvUHv/JKt59MepfWhSxbVvChifhdUeYpLQMTqTrKxMiCtAxMKs5s6UVj4OKL4eyz4dvfhocfdi2YW1tP3vAoGnX3t7cXe+QiMgstA5NKo03AF6ZMjCxooYJ9dR0T35kvvbh8OXz4w/CBD7ggZtcu2L3bdTGLRmHPHjh61GVutJxMpCwslO0pp+YCIpnQKpiFKRMjeUv/jwYu2NHVAyl7C6UXjYHzz4fzznP7wdxzjwtkhobg+c+HF77QdTcTEV+YrblAJe+HI/6nVTALUxAj00TjcfoGBujp6CASyuzXY+b/aDODmoxeV4GPFFOmXSaMcYX+K1YUfkwikpFcgo/Zlptl2zVNQY8Uk5ohLUxBjEzTNzDApu3bAdjY2ZnRc2b+j5bL1YNcAh8REak+ubRsnq3rWLZ1NH5qFS1SDRTEVJFMsh09HR3TbnORy9UDpU1FZqEUpVSZTLIdXhXxZ9tOWc0DRE5WymlKhf1VJJMisUgoxMbOzuNLyeLROP29/cSj8YKObaHmASJVSZWdUmUyaZU8s4i/WEX7xWgVLeI3pZymlImpIrlkOwb6Bti+yS0v6+jpYKBvgMWXL2bfN/bR+tJWWl/U6v1ARcRRilKqTC7ZjqnLvHq6e1S3IlJEpZymKjKIiUfjDPQN0NHTQShSkW8xJ7ks8+ro6Th+mw5ouh7uYuypMQI1AQUxIoWkys7KpaWCs8p2iRdMD3xUtyJSXKWcpipyOVn6y/ZA30Cph+J7oUiIzo2dhCIhOno6OOOrZ2AnLOHVYeJHpy8xG98zziPveITxPeMlGq2IiE9oqaBnpi7z6unuYfNlm+fM5OyJ7uEdN72DPdE9RR6liHitItMUU7MHkp+ZWa3Wl7USvTtK7fJaEsOJacc+/VdPM3CjCxzPvuHsUgxXRMQftFTQMzObAcyXgfmrW/+KGx+6EYAb3nRDsYYoIgVQkUFMOntQaNWwbG1qTUznxk6OPXiMQCgABmzcTjv21H88ddqtiIjMoVhrMKpg2Vo2S8j+8ZJ/nHYrIv5VkcvJiqXYy9aK1YEFTnQla9vQxtrNa+no6cBay/Dvhgm2BSEBgfD0X5+6FXWcfcPZ1K2oK/j4REQkA8VethaNQm+vuy3Sa/Ws2jDvErKpVkRWcMObbmBFRBvYivhdZaYPiqTYy9aKWbCYDtDWbl57PKs1cWiC+FCc8KowidEEJmwKOgYREclTsZetFXPn4tRrRdjMRjXAEKk6CmLyUKxla2nF3GhrtgBtfOc41rolZHbSUrOopuDjEBGRPBS7dVAxgybVFYlUNQUxPpJL68lczRagjTw8QqDWLSGzk5ZQq359RERkimIGTWpBLlLVfFcTU6wd5GU6ay0jD44QbA0CkJxIElqkIEZE5CTFrAsREalSvgtitAdMacSPxJkcnCRQdyITo+VkIiKz0B4wIiIF57tL6doDpjTGd41jjMGYE8X8wUiwhCMSESlTqtUQESk43wUxxS6mF2d8x/i0vJ0xhmCDghgRkZOoVkNEpOB8t5xMSmP0yVECjdN/XZSJEREREZFSUBAzi2poHpDNe7TWEtsVI9g0PWipWayaGBGRkqiG5gHV8B5FJGcKYmZRDc0DsnmPiWiC5GiSQE2qqD9uMTWGYLMyMSIiJVENzQOq4T2KSM58VxNTDNXQPCCb9zhxYGJauJsYSxDuDE8r8hcRkSKqhuYB1fAeRSRnCmJmUQ3NA7J5jxMHJrBJe/zvydEktWfXFmpoIiKykGpoHlAN71FEcqblZHmohtoZgLGnxwiET/yqJMeT1K2qK+GIREQkI6orEZEKpSAmD5VeO5MO0kYfm96ZzBhDbUfttGPyDeSqJSCUPOjLmEj2Kr2uJNPPBS8+P/QZJAvQr0hxKYjJQ0dPB2s3r824dsZvX9TTQVr0juj0PWEMhNpC047JN5Cr9IBQPFDpX8aypdlSMtHTA5s3Z15X4rffq0w/F7z4/NBnkCxAvyLTFfrjRDUxOYhH4wz0DdDR05FV7Uz6izrgi5qbjp4OEmMJhu8ZxgRdEb+1FpJQ015z/Jipt/m8lhfnkQqmIt/p0rMlqG5AThaNut+Rnp7sfj/89nuV6eeCF58f+gySBehXZLpCf5woiMlBrsGI376ohyIhlrxuCSMPjhy/LzmeJLQ4RLAuePwYLwKyamimIHlSke90mi1lPrl+e8jm9yoWg6NH3c+xYxAKQVMTLFkCzc3ZjzkXmX4uePH5oc8gWYB+RaYr9DSlICYHuQYjfvyiPnl40mVfUpLHkjSc3VDCEYkIoNlS5pfrt4eFfq/27YNHH4W774Ynnjhxv7WQbrtvLZx5Jrz0pbBunQtuRKTqFHqa0idLDvwYjORqYv8ETNkOJjmWpP70+tINSEREFublt4dk0gUuW7bAY4+5YKW5Gbq6IDBLaW0yCbt3w1e+AmvXwp/8CSxd6s1YRERSFMTIvMZ3jROsn17UH+4Ml25AIiJSHNbCI4/AjTdCfz80NsKqVScyLnMJBGDxYmhrg7174ROfgI9/HDqr4+KfiBRHVXUny6U7WCk7ipVDN7PYnhiBhim/Jpbj7ZVFRMRjubTzKUQLoCNHXCbl05+G0VFYswba208OYGIx2LrV3c5kjMvABALuXGNj3o1PRKpeVQUxubTxLWbr35lBS7avHY1F6d3aSzTmzURmk5bJg5ME6tyviU1YCEDNkhpPzi8iIjPk0qPV676uDz4I110H998Pp5wCra0nHpsZtGzb5paZbds29/na293ysttu818LZxEpW1W1nCyXgvz0sW0b2ujv7aejp4NQpDD/bDO7nmU73r5tfWza4rrRbFyX/1ro5FgSLJiAu/KWGEkQXhHGBM20NtOF+vcQEak6uRTkp4/dsMEFCD09riYmW8kk/Oxn8L3vQUeH6zQ2UzpoAVe0393t/py+ncvSpfDLX8KuXfDe97r7Ct2YYmqb6Vz+PUSkrFXVt89cCvLTz+nv7S/4Hi8zg5Zsx9vT3TPtNl+JkcS0ov7E0QStL24F/LfnjYiIL+RSkJ9+Tm9v7psyTE7CN78Jd9zh6l5q5si4zwxawmEXzCyksRGeeQYuvji7zTfz4bc9b0QkK1UVxOQjvfFjcixJPBovSPYh365nkXDEkwxMWmJ0ehCDhYZnufbKftvzRiRnuporftHT4+pOxsbc722mv6/xOHzjG3DXXXDqqfMX7mcatMzGGPdTrIBCeylJlajWaaqqamLyEYqECNYHefKaJ4tSH1MOkiNuORmAtRabtIRXu85k6YBLS8mk4nldbyBSKJEI1NfDNddk/vuaSMC3vgV33umK9xfqPJav0dHCnn+qdIaqmr7VSVWq1mlK30CzUG3Zh8RoApt0UUxyJEm4K0yoSb8yUmV0NVf8JNvf15/8BH7zG1fAX+gAxlptfClSANU6TVV1JibbFsbzZR+8bodcqPbK2Zw3MZI4nomJD8VpOm+WIk+RSqeruVJq2XT0mu/3deZ5Hn0UbrrJ1cDMtmnlQuZrrzwbY9z+MQtRBzORrFTrNFXVQYyX7ZO9bsVcqNbO2Zw3cSwB6X0uEyfqYUREpIi8Wisy9TxDQ/DVr8KSJXMX8S8kk/bKacmku21ry26cIiJzqOq8rpfLw7xealaopWvZnDc5lsQEjFtSZqD+tPppj6vNsohIEXi1VmTqeb73PdcAIJPMyFwyba8McPgwnHNOZsvJvFwbU60VzyJVoKozMVOXh2W7fOvIbUd4+uNPc+C7B4jtj03rLObFMrBCFc5nc1474Ta3jB+NU39aPcHG4LTHi7kRqIhI1Zq6ViSfpVbp8+zc6dop5xPAwIlOZeHwwseOjsJFF2V2Xi/XxiirI1KxdPk8JZt9TyYOTLD/3/cTag1x5NYjDP16iM4/66T5Oc0VtX9KsCWIjVsSRxO0vKnlpMerrdGBiEjJ5bv3ibXwyU/Cb3/rlpLl2i45GxMTbsnaWWcV/rVmqtaKZ5EqkFcQY4z5DPA6YAJ4CrjKWjvkwbiKLpsv5JNHJjEBQ6glRKglRGIkwZ5/3sOKa1dU1Bf7UEsIG7ME6gM0r28++fE897URESm0SpqngPy/lA8MuMzJpZdmtgzMC/v2wetf79o/F1sum4eKiC/ku5zsV0C3tfY5wHbguvyHVBrZLLNKjp7YPwUg2Bikpr2G/q/0kxxLVsz+KaFIiMRogkWXLDppKZmIiE9UzDwF5L/U6r77XBCzfn1my8DyNTbmXueSSwr/WiJSVfIKYqy1v7TWpos/fgusyH9I5WW2Wplg5OQv9MGmICZg2POVPSTGE8UcYsE0djey5m/W0P7G9lIPRUQkJ9UwT2VcJ2Mt/PrXmXUI84K10N/vskbNJ2fzRUTy4WVh/7uBn831oDFmozHmXmPMvQcPHvTwZQtrtuL12mW1bgd7a6cdW7O0htgzMQ7/9HCxh+mJmQGbCRqaz2smUFvV/R9EpHJU5DyVcfH6yAgcPAiNjcUZV38/PPe5mRf0Z0r7yIgIGdTEGGNuBZbN8tDHrLU/Th3zMSAO3DjXeay1vUAvwPr16+1cx5Wb2WpcQk0hwl1hEsOJk5aN1XbVMvjTQVpe2EJ4eeap+ng0zr7r92EwLLtqWUmWo1VSUwIRqR7VPk9lXCdz+LDb1NKY3F4nFoPf/979+fzz51+OdvSoq4F597tz20hzPvk2NxCRirDgN2Vr7bwLWY0xVwKXARfbmamJCjBX8XrrK1oZ+M7AScFGoCaAqTEc+uEhuq7uyvh1BvoGeOqap9w56gMlCSIqqSmBiFSPap+nMi5eP3LELfHK1bZt8ItfuD/X1Mzd2WxkxG2med110HJyZ8u8qeOYiJB/d7LXAB8GXmqtHfVmSP4QWRdh4IYBbMJigtOvatUsrSF6d5TFr1tM3cq6jM7X0dNBYiyBwZQsiFC3MRGpNNU8T50kGMw9CwOum9nk5Ik/z2Z01C1Z+8AHYO3a3F9rPuo4JiLkXxPzZaAZ+JUx5n5jzNc8GJMvhFpCNK9vZvLg5EmPmYDBhAxH7zia+fkiIVa+fyUr3r+i5J3Nst34U0SkjFXtPHWShob8nh8Ow4UXup/ZlpKNjMCBA/Bnfwbnnpvfa2VCtTEiVS3f7mSnW2tXWmvPS/2816uB+UHbq9tIjiWx1pKMJRneOkwylgSgpqOGod8MzdmprJwDhdmaGYiI+FG1z1PTtLW5upZ773W3C4nFYOvWzI49dMgtIfvAB1z75mLItKGBiFQk/29mUkJ1p9bRvK6ZkW0jxPbFGNwyCEDzOtfRKzmRZOShESLPO7mffzkX0as2RkSkAi1a5LIWP/2pW1Y2V01L2rZtsGWL+/Ncx1oLe/a4AOm666Ar81rQvKk2RqSqKYjJgzGG9re0M3z/MPXPqqeNNhq7T7SuDNQFOPb7Y7MGMeUcKCxUGxOPxhnoG6Cjp6PkS99ERCQL114LAwNw9tkLH5uue5mr/mV42NW/XHABXHklNDV5NsyMzFcbE426DE1PT+4bg4pIWdMGIDmYuhQsvCzMksuXED8Yp+m5TQTC7p80GUsytmOMY78/Nus50oGCH4MALTcTESlzc9WLXHQRvOc9rnbF2vmXjIXDLgMzs/5lbAx27oR4HK65Bv70T4sfwCxES81EKp7/vkFnyeusQTwaZ/vV2xm40X2B79zYyeJLFzO2fYzR7aOEu9yH/ci2EYZ+McTYKWPE9scIL8t8z5hyV85ZJBER3/E6axCNwtVXw42pLXGmZiuMgSuucDUsjzziNqT87/92XccuvHDuc1rrznvkiAtY3vpWeNnL3F4w5UhLzUQqXsVnYrzOGgz0DTBw4wAdV3Qc/xIfCLl9XYKNQSYPu25ljd2NNJzTQGxHjIPfn3/n53Iu8p+Nn7NIIiJlx+usQV+fC2CuuGL2L/H19fCXf+keGxlx942MwMSEC1asdUHNsWNu6dmuXfDkk/D44/DOd8LnPgevfW35BjBwYqmZlpKJVKyK/xbqddZg6vmmfokPtYRYee1Kdv/Tbib2T1CztIa217QRag6x9J1L5z1nORf5i4hIgXmdNZh6vrm+xAeDsGGDa4X8pS+57MrAgKtzMcYFKIsXw+mnw3nnwZ13unqaN7wBLr7Ym3GKiOTBlGLz4vXr19t777236K9bDBMHJ+j/Wj/jO8ZJxpO0XdLGsj9aNu9zvFzypqJ7EcmUMWartbZI/XD9pZLnqXklk5BIQE3N9Pu9XPKmonsRydB885S+5Xqstr2W1R9bzcSBCeKDcRrOPHlzsZmBxkLdwLKhrI6IiOQsEHA/cHKwMVcnsGyll8+Bd+cUkaqjIKYATMAQXh4mvHz2Yv5CBhoquhcREU8UKthQ0b2IeEBBTAkUMtDwMqsjIiJVrFDBhpdZHRGpWhXfnawYsu0upu5eIiJSdHPtHTMXdfgSkTKmIMYD2bZxnhn0+K3FsoiI+FA2rZxnBjzZBkAiIgWmVIAHsl0eNrMmJt8aGXUkExGRBWWzPGxmPUy+9THqSCYiHtM33hzk211sZtCTb42MOpKJiMg0swUN2dSizAx48q2PUUcyEfGYgpgceB005FuMr45kIiIyjddBQ77F+OpIJiIeUxBD9suxyi1zoo5kIiIVLtvlWOWWOVFHMhHxWNUFMbMFLNkGFcqciIhIwcwWsGQbVChzIiIVruqCmNkClmIHFcqciIjInGYLWIodVChzIiJlruqCmNkCFi+DCnUKExGRvMwWsHgZVKhTmIhUgKrbJ6bQG01mu2eMiIjINIXeZDKb/WJERMpU1acKvM6cqN5FREQ85XXmRPUuIlIBqi4TM5PXmZNCZ3pERKTKeJ05KXSmR0SkCKr+m3Y5ZU5UTyMiIicpp8yJ6mlEpExUfSamnDInqqcREZGTlFPmRPU0IlImSv/NXY4rp6yQiIjIScopKyQiVU1BTBnR/jEiIlLWtH+MiJSJql9OJiIiIiIi/qIgRkREREREfEVBTIHFo3H6e/uJR+OlHoqIiMjJolHo7XW3IiI+oSCmwLLtOKagR0REiiqXjmMKfESkxFTYX2DZdhxLBz2AivxFRKTwcuk4lg58QIX+IlISCmIKLNuOY2qzLCIiRZVLxzG1WhaREtNysjITioTo6OlgoG9AS8pERERERGahIKYMZVtHIyIiUlS51NGIiHhIy8nKkJaUiYhIWdNyMhEpMQUxZSjbOhoREZGiyqWORkTEQ1pOJiIiIiIivqIgRkREREREfEVBjIiIiIiI+IqCGBERERER8RUFMSIiIiIi4isKYkRERERExFcUxIiIiIiIiK8oiBEREREREV9REOOxeDROf28/8Wi8IMeLiIjkJRqF3l53W8jniIgUkIIYjw30DbB903YG+gYKcryIiEhe+vpg0yZ3W8jniIgUUKjUA6g0HT0d0269Pl5ERCQvPT3Tbwv1HBGRAlIQ47FQJETnxs6CHS8iIpKXSAQ2biz8c0RECkjLyURERERExFcUxEhG1IBARETKmpoPiFQVBTGSETUgEBGRsqbmAyJVRTUxkhE1IBARkbKm5gMiVUVBjGREDQhERKSsqfmASFXRcjIREREREfEVBTEiIiIiIuIrCmJERERERMRXFMSIiIiIiIivKIgRERERERFfURAjIiIiIiK+oiBGRERERER8RUGMiIiIiIj4iidBjDHmWmOMNcYs8eJ8IiIiXtNcJSJSOfIOYowxK4FXAc/kPxwRERHvaa4SEaksXmRiPg98GLAenEtERKQQNFeJiFSQvIIYY8zlwF5r7QMZHLvRGHOvMebegwcP5vOyIiIiGct0rtI8JSLiH6GFDjDG3Aosm+WhjwEfxaXnF2St7QV6AdavX68rYSIi4hkv5irNUyIi/rFgEGOtvWS2+40x5wCnAA8YYwBWAPcZYy6w1u73dJQiIiLz0FwlIlJdFgxi5mKtfQjoSP/dGLMTWG+tPeTBuERERPKmuUpEpDJpnxgREREREfGVnDMxM1lr13h1LhERkULQXCUiUhmUiREREREREV9RECMiIiIiIr6iIEZERERERHxFQYyIiIiIiPiKghgREREREfEVBTEiIiIiIuIrCmIqUDwap7+3n3g0XuqhiIiInCwahd5edysikgMFMRVooG+A7Zu2M9A3UOqhiIiInKyvDzZtcrciIjnwbLNLKR8dPR3TbkVERMpKT8/0WxGRLCmIqUChSIjOjZ2lHoaIiMjsIhHYuLHUoxARH9NyMhERERER8RUFMSIiIiIi4isKYkRERERExFcUxIiIiIiIiK8oiBEREREREV9RECMiIiIiIr6iIEZERERERHxFQYyIiIiIiPiKghgREREREfEVBTEiIiIiIuIrCmJERERERMRXFMSIiIiIiIivKIgRERERERFfURAjIiIiIiK+Yqy1xX9RYw4Cu4AlwKGiD6DwKvF9VeJ7gsp8X5X4nkDvqxBWW2vbS/TaZU3zlC9V4nuCynxflfieQO+rEOacp0oSxBx/cWPutdauL9kACqQS31clvieozPdVie8J9L6kNCr1v08lvq9KfE9Qme+rEt8T6H0Vm5aTiYiIiIiIryiIERERERERXyl1ENNb4tcvlEp8X5X4nqAy31clvifQ+5LSqNT/PpX4virxPUFlvq9KfE+g91VUJa2JERERERERyVapMzEiIiIiIiJZURAjIiIiIiK+UjZBjDHmWmOMNcYsKfVY8mWM+Ywx5jFjzIPGmB8aY1pLPaZ8GGNeY4x53BjzpDHmr0o9nnwZY1YaY35tjHnEGPOwMeb9pR6Tl4wxQWPM740xW0o9Fq8YY1qNMd9P/X/1qDHmBaUeU76MMR9I/f5tM8Z81xhTV+oxyfw0T5UvzVP+onnKH8p9niqLIMYYsxJ4FfBMqcfikV8B3dba5wDbgetKPJ6cGWOCwFeA1wJnA28zxpxd2lHlLQ5ca609G7gQ+NMKeE9TvR94tNSD8NgXgZ9ba88EzsXn788Y0wX8BbDeWtsNBIGe0o5K5qN5qnxpnvIlzVNlzg/zVFkEMcDngQ8DFdFlwFr7S2ttPPXX3wIrSjmePF0APGmtfdpaOwH0AZeXeEx5sdbus9bel/rzMO6Dpqu0o/KGMWYFcCnw9VKPxSvGmBbgIuAbANbaCWvtUEkH5Y0QUG+MCQENQH+JxyPz0zxVvjRP+YjmKV8p63mq5EGMMeZyYK+19oFSj6VA3g38rNSDyEMXsHvK3/dQIR+kAMaYNcD5wO9KPBSvfAH3RStZ4nF46RTgIHB9avnB140xjaUeVD6stXuBz+Ku6u8Djlprf1naUclcNE+VPc1T/vIFNE+VPT/MU0UJYowxt6bW0838uRz4KPA3xRiHlxZ4T+ljPoZLCd9YupHKXIwxTcAPgGustdFSjydfxpjLgAFr7dZSj8VjIeC5wL9aa88HRgBfr3k3xizCXSk+BegEGo0x7yjtqKqb5inNU+VI85RvaJ4qgVAxXsRae8ls9xtjzsH94zxgjAGXzr7PGHOBtXZ/McaWq7neU5ox5krgMuBi6+/NePYCK6f8fUXqPl8zxtTgJoYbrbU3lXo8HnkR8HpjzAagDogYY26w1pbVh04O9gB7rLXpq5Dfx+eTA3AJsMNaexDAGHMT8ELghpKOqoppntI8VW40T/mK5qkSKOlyMmvtQ9baDmvtGmvtGtwvwXPLfWJYiDHmNbhU6euttaOlHk+e7gHOMMacYoypxRV13VziMeXFuG8i3wAetdZ+rtTj8Yq19jpr7YrU/0s9wP9UwMRA6vNgtzHmWam7LgYeKeGQvPAMcKExpiH1+3gxPi8CrVSap3xB85RPaJ7ylbKfp4qSialCXwbCwK9SV+5+a619b2mHlBtrbdwY82fAL3CdKb5prX24xMPK14uAdwIPGWPuT933UWvtLaUbkizgz4EbU19QngauKvF48mKt/Z0x5vvAfbilPL8Heks7KqkymqfKm+Yp/9E8VWTG3xlkERERERGpNiXvTiYiIiIiIpINBTEiIiIiIuIrCmJERERERMRXFMSIiIiIiIivKIgRERERERFfURAjIiIiIiK+oiBGRERERER85f8HAG72ymlClmoAAAAASUVORK5CYII=\n", "text/plain": ["
"]}, "metadata": {"needs_background": "light"}, "output_type": "display_data"}], "source": ["import itertools\n", "import matplotlib as mpl\n", "color_iter = itertools.cycle(['r', 'g', 'b', 'c', 'm'])\n", "f, axarr = plt.subplots(1, 2, figsize=(14,7))\n", "\n", "for i, (clf, title) in enumerate([(gmm, 'GMM'),\n", " (dpgmm, 'Dirichlet Process GMM')]):\n", " splot = axarr[i]\n", " Y_ = clf.predict(X)\n", " for i, (mean, covar, color) in enumerate(zip(\n", " clf.means_, clf.covariances_, color_iter)):\n", " v, w = np.linalg.eigh(covar)\n", " u = w[0] / np.linalg.norm(w[0])\n", " # as the DP will not use every component it has access to\n", " # unless it needs it, we shouldn't plot the redundant\n", " # components.\n", " if not np.any(Y_ == i):\n", " continue\n", " splot.scatter(X[Y_ == i, 0], X[Y_ == i, 1], 0.8, color=color)\n", "\n", " # Plot an ellipse to show the Gaussian component\n", " angle = np.arctan(u[1] / u[0])\n", " angle = 180 * angle / np.pi # convert to degrees\n", " ell = mpl.patches.Ellipse(mean, v[0], v[1], 180 + angle, color=color)\n", " ell.set_clip_box(splot.bbox)\n", " ell.set_alpha(0.5)\n", " splot.add_artist(ell)\n", " splot.set_title(title)"]}, {"cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [], "source": []}], "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.9.5"}}, "nbformat": 4, "nbformat_minor": 2}