Skip to content

Instantly share code, notes, and snippets.

@KOLANICH
Last active March 9, 2022 16:45
Show Gist options
  • Save KOLANICH/d9583e821f5e9592d9ffe9d26f63e44b to your computer and use it in GitHub Desktop.
Save KOLANICH/d9583e821f5e9592d9ffe9d26f63e44b to your computer and use it in GitHub Desktop.
Seaborn jointplot with fitted gaussian underlay
Display the source blob
Display the rendered blob
Raw
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"id": "broke-detector",
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import scipy.stats\n",
"import scipy.linalg\n",
"from scipy.spatial.transform import Rotation as R\n",
"from matplotlib import pyplot as plt\n",
"import seaborn\n",
"\n",
"%matplotlib inline"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "spoken-handy",
"metadata": {},
"outputs": [],
"source": [
"def plotFunc(f, xR, yR, ax, **kwargs):\n",
" \"\"\"Utility function to plot 2D functions\"\"\"\n",
" xs = np.linspace(*xR)\n",
" ys = np.linspace(*yR)\n",
" dd = f(np.dstack(np.meshgrid(xs, ys)))\n",
" m = ax.pcolormesh(xs, ys, dd, shading='auto', **kwargs)"
]
},
{
"cell_type": "markdown",
"id": "sharp-ottawa",
"metadata": {},
"source": [
"# Defining distribution"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "gothic-nepal",
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAD6CAYAAABEUDf/AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Z1A+gAAAACXBIWXMAAAsTAAALEwEAmpwYAAAeR0lEQVR4nO3dbaxlV1kH8P9z3u65L3M7TDtlSlugKP1QfKlYChoQkNqiiSmtlmCC8skxWAwYTAWaYEkgUaJUg0YZpQGDShAYSlJMdYgRNQqW2tKWhlClhbZD22mZzn0/L/vxwz1j52X9nzt3nXvuObPu//elc/e6a+999z2zZvdZ63mWuTtERKRMtXHfgIiIjI4GeRGRgmmQFxEpmAZ5EZGCaZAXESmYBnkRkYKNfJA3s3ebmZvZeaO+loiInKwxypOb2cUArgbw3TP5/pZNeRuzo7wlORuZ8aawX0avWtDG7iO4P1jwHhW9YpFzenTv4cPYPItyaFhT2Cdoq0hb7vnC/B92reh0k51PtIAfHHH3vam2kQ7yAG4FcBOA28/km9uYxSvtDaO9IxmvaNBjXer14HTByEb6WYN/7K3V4udrNTd3HACC8/kUvw9vpu/dg2eBOn8WTv7RiAZy61X8WqTNOj1+vk6Xn2+tkz7e5X08aENwLe/30w3seNQHgNN/oILnt8UO+WcfYW0jC9eY2bUAHnP3e0d1DRERiQ31Jm9mhwDsSzTdDOB9WA/VbHSO/QD2A0AbM8PcjoiInGKoQd7dr0odN7MfBXAJgHtt/X8TLwJwt5ld6e7fP+UcBwAcAIB52zPZgS8RkbPMSGLy7n4fgPOPf21mDwO4wt2PjOJ6UgYWX8+JuwNBLL8ZfOwbQcyb9Yvi7m0er/cWv4+qRWLyDR5h9eg50U78vaoWnM9YjL8K4tBV8Gx75OfqZ0aUo8nwHUbr5EVECjbq1TUAAHd/8XZcR0RETqY3eRGRgm3Lm7zsMMFa+DC+zvpFa+uj2CtZDx+tk8dUsE5+aip5OIq7V9NBW7BOvmqlf+YqWgsfxdBJ6N3YGm8gfLbsN+JRDJ3F3QEY/d1H8wLB/AS/ix1Hb/IiIgXTIC8iUjAN8iIiBVNMXiYGjddHseYmj3kbW9ce9InaWOw9irv32/yvWL/N143TmHwjr0CZkeXr1s+MXpP19daN8hZ4/RfUNx+TlzOjN3kRkYJpkBcRKZgGeRGRgmmQFxEpmCZeJc9WJzwFbdGmIXTCDqDJUFHCU05iUz+YeO1N83vvt/m996fIzlBBMlSEJT3VyF4dg9ZNny8qoIZaZlIbPZ8mZc+E3uRFRAqmQV5EpGAa5EVECqaYvGyreAMQ8s4RbbwdFRvL2Hjb2zxez2LvvRkedw/bpvmz6JNbzNoYBECN5CFFIX6WQAXwQmm1IFmLbSYO8E1IlAw1PL3Ji4gUTIO8iEjBNMiLiBRMMXmJsbXruWvhMzbeDtfJR8XGyAbbYUGxmWDNO4mvd2f5/XVng7h7e/Mx+fC1LIih17ubj23XerzNm+kbCecMogmACVjz7tEGKmcxvcmLiBRMg7yISME0yIuIFEyDvIhIwTTxKlsunpSNJubIBGaQvJRTbKwKipD1ZvhfCTbB2p3jP1N3JphcnaJNqNgtRo822HjJ2dyw8xPWuvx89TXWKfN3LyOjN3kRkYJpkBcRKZgGeRGRgo00Jm9mvwXgRgB9AHe4+02jvJ5kytkAZIsTngDAmuTjyI4DYbyeJT31Znmf7iz/uVjsvRPE5HsztAn9Nm+jMfSABclLLA8pTHgK7oElPUVFyCZGoUlPzMgGeTN7PYBrAfy4u6+Z2fmjupaIiKSNMlzzdgC/7+5rAODuT47wWiIikjDKQf5SAK8xs6+a2b+Y2StS32Rm+83sLjO7qwu2LktERHIMFa4xs0MA9iWabh6cew+AVwF4BYDPmNlL3P2kgJi7HwBwAADmbc/OCpadxbLXwkfxdVZsjBQaA4BqKig21k5fqzcTxN2DmDyLvXfnaJcwJl9N8Y87X9fOz1frBJtykH5VsBa+Cn5Vzh5T8Nq45fH6KLbuGkqOG2qQd/erWJuZvR3A5weD+tfMrAJwHoCnhrmmiIicuVGGa74A4PUAYGaXAmgBODLC64mIyClGuYTyNgC3mdn9ADoA3nZqqEZEREZrZIO8u3cAvHVU5xcRkY2pQNlOkZPwFPXLTXhqRBOv6Taf4n2qcCendL+chKf1NnJ8F+2C/izfrqlqRrOo7IS8i9f4z1XrpX+uKkh4qoKdnLYtGWpC/uf/bN41SmUNREQKpkFeRKRgGuRFRAqmmLyEeIGyzISnMCafjq97lPAUbPLRI5t8dOaihCfaRGPvvV087u6zQRC9wfsZyV7yLg+iV8EGIFWTxNCjESB6BcwIvVsUX8+Iebvz55cl93xbfR9bTG/yIiIF0yAvIlIwDfIiIgXTIC8iUjBNvJaGJC9lJTxF58tMhoomZb2dnmCtyHEA6M3wa9GqkbO0S5jYxCZYfZ6XcmzN8LZ6PZiwJZOonTU+QdnvBxOvq2TiNUp4moRNnjKTobwKJkMnfKJ0q+lNXkSkYBrkRUQKpkFeRKRgislLHK+vb75AGd3hCQh3eWJJT6zQGAB0w12eSEx+nnZBbz6Ik5PY+/Quvm3lXJu3NYKYfK+f/rkW0KZ9VoNEKW+k2+gOT5nChKecti1OoNqJ9CYvIlIwDfIiIgXTIC8iUjDF5M9G0br2nD7RhiIk9m5REbJoLXywAUh/Ot3Wy4i7A9EmHzyWW833aNvM/Gry+Hm7lmifc1rpPgBQI0XIAGCxm5676Aa7fKyt8LkQr7FrZS6GJ6ezfhAnj9rYuvZwvXte/P9s3gAkh97kRUQKpkFeRKRgGuRFRAqmQV5EpGCaeC0MS2wKE57CZCgy0ZexwxOwwS5PZOK1OxtMvJIiZOtt5Dq7+G5NU/M8een8XYvJ4xfMHqN9djX4+apg0rNRm0keX+pM0T4L9WBCmbZwwbwwjE1eRn1yEpvCPsFPlbPTVKETsnqTFxEpmAZ5EZGCaZAXESnYyGLyZnY5gL8A0AbQA/Cb7v61UV1PNpCR8BS2BX2ihKcqaOtNp+8xLkJGm2ixsXqwyUeU2HTx3NHk8Qvaz9I+M/UObVut+PwES3pqNXiyltUyIu9RGDo4nZE2GqsHNoihk3mS3IJnW+0s3mhklG/yHwbwAXe/HMD7B1+LiMg2GuUg7wCOF3U9B8DjI7yWiIgkjHIJ5bsA3Glmf4j1f0x+eoTXEhGRhKEGeTM7BGBfoulmAG8A8Nvu/jkzezOAjwO4KnGO/QD2A0Ab6bXBO1IUQ8/ZlDtnLTwAkM0m0ArWwreiImT8WqwQWS+KuwdtFVkPv3vXCu1z0RyPr79w+pnk8QtaR2mfpvE1+c/2+ef9WCO9OUidBcMzRaerBQXFWJv1ghOGBcpIW+7a9SiGfhbH13MMNci7+2mD9nFm9tcA3jn48u8B/BU5xwEABwBg3vaUmY0gIjImo4zJPw7gtYM//yyAb4/wWiIikjDKmPyvA/gTM2sAWMUgJCMiIttnZIO8u/8bgJ8c1flFRGRjKlBWGDopa9HEaxC1Y4XIgonXKihC1gsmXrvT6XvsBvPx3V18Eq21K10c7PlzC7TPi2aepm2XtJ9KHt/b4Ofre97OS1O1c5LHo92kPLiWVaRwHZ8Xjtt6ZOK1z38fFiZDpdu8z2/Co/PJ/1NZAxGRgmmQFxEpmAZ5EZGCKSZ/NgoSpVhbVhEygCZDRQlP1VSQ8NTm996bSceNo4Qnn+Mx23Nm00lPF88epX1e3D5C2y5upuP1u2s8uWrV+fzEQjVN2+ok9l4FcXfvBzF5UtesxuudhW3GkqGihKdeEORnsffMAmWlbgCSQ2/yIiIF0yAvIlIwDfIiIgVTTH7cWAw9KigWYf1qGWvhAbopdxUVIWtHRciCjbfJevjeHI+vNuf4phz7yHr4F7bThcYAHncHgBc00ht2t4MqX0eD2HBUvIytr2ebiQBA1eO/42Y3fb44Jh8UKOuStrBAWVQ0jBUo29rNuoEgXl9o4TK9yYuIFEyDvIhIwTTIi4gUTIO8iEjBNPF6Fop3hiJtURGyIBnKm+mPiLd4n36U8ESKkAEA2yipmuETlLtnV2nbBdPpidIXTfGEpwsbR2nb3trmJ+aWESQABXpkgrXT439lPZh4rXU3dxwA6p1g4pVMsFovmMmtomQo8mzDZCjt/nQm9CYvIlIwDfIiIgXTIC8iUjDF5CdVRhEygBcbsyjhKWhjhcj6URGyIO7e4zW50JtJx19rszxwfO7MEm17Qfto+njjB7TP3jq/1q5aK3m8mxn/Xa148bLFfvpaa72g0FwnIybPc8nCZCjrsph88CyCAmVscxBtDDI8vcmLiBRMg7yISME0yIuIFEwx+e0QxdBzCpFFfdia94yNQQDAm+m2/lSwFr4dxOTb/Db60+kY8PQMDxyfP71I2y5oHk33qfM+5xiPk0/bVPJ41/mmIR3wZ7tcpePuALDYTT+o1Q6/v1qHP/daek/zeC38Go+H17okvp6zMQjA18NrY5Ch6U1eRKRgGuRFRAqmQV5EpGAa5EVECqaJ1wmVVYQM4JOyGUXIAF6ILCpC1p8KipAFyVCsENn8NJk1BHBei0+i7iU7OZ0TbIc0U5ulbUwXfIJyueI/8DO9Odp2tJOeeO2u8t9VbZX/Tupk7jqaeK13+M9ldOI1KlAW7fJE2rZ696f1xqxznq2GepM3sxvM7AEzq8zsilPa3mtmD5nZt8zsmuFuU0REcgz7Jn8/gOsBfOzEg2Z2GYC3AHgZgBcAOGRml7p7Xt1VERHJMtSbvLs/6O7fSjRdC+DT7r7m7t8B8BCAK4e5loiIbN6oYvIXAvjPE75+dHDsNGa2H8B+AGiD7BpRMpYolVGELGyLNg0hCU8AULXS/aK4e5jwlM4nAgDU2+l47q4pvjHIuU1eoOzcerptthYU8greeyoSe18KYs1P93mM/+kub1tYSz/EfhCTb/KpCzTII6yvBclQQUweLCbPNv8AL0IWtmljkKFtOMib2SEA+xJNN7v77cPegLsfAHAAAOZtj1LYRES20IaDvLtflXHexwBcfMLXFw2OiYjINhrVOvkvAniLmU2Z2SUAXgrgayO6loiIEEPF5M3sOgAfBbAXwB1mdo+7X+PuD5jZZwB8E0APwI3Fr6zZziJkQUyZtkUbgwQx+T6LybeCtfC87haqNo+jNklMfneLx+T3NPg6+RmyI0Yz891muUqf79mKP9unevO07UiHr5M/tpKevKgt899VnddJo7H3+mqwFr7D/8pal6yHj4qQhWvX020qQja8oQZ5dz8I4CBp+xCADw1zfhERGY7KGoiIFEyDvIhIwTTIi4gUTAXKxoxOykZFyKI2lgwV7P5URclQZAeoaOK1ChKeqhafSGu3usnjc00+8comVwGghc3P9a84zyh61tP392R/N+1zuMvbnlzhE6+ry+nZ6/oKf+4s4Wm9jUy8Brs/0SJkAN8BKpp4zZ2UJVSE7MzoTV5EpGAa5EVECqZBXkSkYIrJTygLEp6iAmUs9u4Nfj5vBkW5GukYcBUlPDV5m7d4rHSqkY7ZzrEdLwC0a+k4eWQ1yMuL2p7qp3+wx3vPo30eX91N255e5gXKqqX0tVrLQUw+SIZqrLCYfJDw1Ak2ACGbg3iwaUhUoIzG0BVbH5re5EVECqZBXkSkYBrkRUQKppj8dgiKl/FNQ6K18FGBMhaT57/qKojXs/XwUdw9arMGj7G2Gul4bqOWV9uug/SzWAhj8vxZPNbbnTz+yNp5tM/hFV6gbGGJ767CCpHVl2kXGncHeCGyrCJkAF/zHmwawoqQASpENkp6kxcRKZgGeRGRgmmQFxEpmAZ5EZGCaeJ1s8hEadbuTwDfASra/SlnUjaYXK2aQbExlgwVfHKqJp9EqwUTr00ywVo3fr6u88SwZZKxVTn/eY9W07Tte91zk8e/u7qH9vn+4i7a1l3kGWWtxfQ9NqOJ1+Vo4jX9bGtRwlM3SDRjSU9VMPGqImRjoTd5EZGCaZAXESmYBnkRkYIpJr8Nwng9i6/nbAwC8AJldX6+qI3F3qOYvAdtVg/i9ST23g9i6KtB5tXT/XQBsCjG/1SPx9AfXk0nPX1vaTft8+wSj/HbIv89NpbI8SDu3iAJTwBQI4XI4iJkPFGKFRvLKkK2UZsMRW/yIiIF0yAvIlIwDfIiIgXTIC8iUjBNvI4Z2wEq3P0pSJRy0uZB5UonCU/r/TZ3fP0eNj+5CvAkpbVglvfZ/gy/ESJKoDrc2U3bHllOJ0M9ESQ8dRamaFtrgf9OWNJTM6g02Vjhk561NZLY1MlIeAJ4tUlVmpw4Q73Jm9kNZvaAmVVmdsUJx3/OzL5uZvcN/vuzw9+qiIhs1rBv8vcDuB7Ax045fgTAL7r742b2IwDuBHDhkNcSEZFNGmqQd/cHAcBOWdPt7v99wpcPAJg2syl3XxvmeiIisjnbEZP/JQB3Fz/A5+z+FIkSqKI2ktjkUYGyKBmKnW8EU/Z9ctKlHo9rP9NLJzwBvEDZcp8XBju8eg5te3Qx3bawwBOeagubT3iK2hrLPGmoHsTkadJTGJMPEptY0lNmQTEar1eS1NA2HOTN7BCAfYmmm9399g36vgzAHwC4Ovie/QD2A0Abm59EExERbsNB3t2vyjmxmV0E4CCAX3P3/wnOfwDAAQCYtz2afhcR2UIjWSdvZrsB3AHgPe7+76O4hoiIbGyomLyZXQfgowD2ArjDzO5x92sAvAPADwN4v5m9f/DtV7v7k0Pd7XbJiaFnXyujQFlt82voPYjjh/F11pa5R4pX/GKdXvrjuBjE5I905mhbDen/MVwIznd4aZ5fayEd/+8f40XSWseCtfBBTL5JYu+NlagIGV/XbqwtKCgWFRujbYqhT5xhV9ccxHpI5tTjHwTwwWHOLSIiw1NZAxGRgmmQFxEpmAZ5EZGCqUDZJoW7PDFhYhP5dzaceA0mUelELj9dOCmbOcHKVP1glycy8Xp0rU379IKJXFbw7NkOPx+bXAWAtWPpCdt6kPDUXKRNaC7yFcPNpfQEZn0lmlwNEpu6pC0sQhYkQ5HkpagImQqUjYfe5EVECqZBXkSkYBrkRUQKppj8Fglj9VF8nbUFG4OEbaygWFCELDexiZ6uCgqe9fm9r66lk4qO1XgMfbnLi431yLUWV3ky1EqwyUftWPqvS3OB/7xhTH4p2ACEJEPFMfkgvt5Ntzk5DsTJUDTpKTcZSklUI6M3eRGRgmmQFxEpmAZ5EZGCKSY/biyWnxPHD9ro+vlMFoRQLQgNeycoULaW/jgugMfkLdgYvE9i8p1lXlDMgmJjDVJsrHmMdgnj7qwIGQA0ltMPMVwL3+nwNrZOPoq7s826oTXvZxO9yYuIFEyDvIhIwTTIi4gUTIO8iEjBNPG6VYLdpCxKXqKdtjhDKROb17Rgvq7W5fdua0FBsXr647gWFDXzqIJaL30tW+IFxRoLwU5OC+R4sMNTKyhC1ljmD7G2mp4oNXIcAE14AnhiU1bCU9AWTsgq4Wks9CYvIlIwDfIiIgXTIC8iUrCdG5MPYugj6cdOl3O+nE1DIlFeCwmjRglPtSBsXOsExcvIs4gSqNj9rd9H+lr1pc3H3QGgRYqNtRaDpKYlHvOuk4QnIEh6YklNQLwBCIvXuzb5KJ3e5EVECqZBXkSkYBrkRUQKtnNj8oGszbq3/Ca2uKBYEHu1IPZa66fbasHa9SgmX18N1tCzTbmjemzRev21dMfGMu+Ts/E223R7/Vo8Tl5b5g+KbgASrYUPYvJ0PXzmOnkar9da+ImjN3kRkYJpkBcRKdhQg7yZ3WBmD5hZZWZXJNpfaGaLZvY7w1xHRETyDPsmfz+A6wF8hbR/BMA/DHkNERHJNNTEq7s/CACWmCQ0szcB+A6AoHzTDhFNom7xJC+dYI12cooSishcXpjwtMbb6tGj6G3+WUQTr/XV9PFo4jUqKMYmWBvBDk/11WDiNWeXpygZKpiUZUlPSngq30hi8mY2B+B3AXzgDL53v5ndZWZ3dRGMDiIismkbvsmb2SEA+xJNN7v77aTbLQBudffF1Fv+idz9AIADADBve/RaISKyhTYc5N39qozzvhLAL5vZhwHsBlCZ2aq7/2nGuUREJNNIkqHc/TXH/2xmtwBYLH6An4QEqkBuMhSLeYcJTyScvH4fvM3ZXh5Bn/A+SPSvuRzE3XNi8kv8JmorUdx98/F1z9gYBABPespJeNqgn0yWYZdQXmdmjwL4KQB3mNmdW3NbIiKyFYZdXXMQwMENvueWYa4hIiL5lPEqIlIwFSgrDV0nHxUh46erd9L9qmZwD8GKqihsTF85gvBvPSyGxgqKRXF3foOs2BjbdBvYYOPttWDygsXrw4Ji2gBETqc3eRGRgmmQFxEpmAZ5EZGCaZAXESmYJl43aSJ2jYom0chOTmHCU5fPbNa66Z+XJRptxHKKkAX3Hk28NlbS/RrL0eQqb6uTnZzihKdgcjUoNsYSm7ISngCavKSEp/LpTV5EpGAa5EVECqZBXkSkYIrJT6ogsSVqY4XIrBclQwVtpF8Ukw+LkGW8Vlhwf/W1oG2VFBRbCeLuQXzdVkl8PSo0FrR5L9jkg7X1MwuKyY6lN3kRkYJpkBcRKZgGeRGRgikmP245cdQoXl+lY7bW47HcaJ18neYF8PeDqOCZ16KAPTlfNJ/QiTbRTt9IuLl2GJMnbVGhsSjuHsXXMzb5yNoARGvhi6c3eRGRgmmQFxEpmAZ5EZGCaZAXESmYJl4nVTS5GkzY0QJlQfGqWjf4tz7Y5YnxoIibB6djhcjCZK1g4rW2RnZyyplcBXhBsajQWJdPvHrQj03IK+FJNktv8iIiBdMgLyJSMA3yIiIFU0x+zJwko7BCY+udgsQmkgyFzGQo+hYQ3V9GHB/gMfkokavW4XMN1knHw8O4+1pQeY0kPWUVGgM2KEJHfuachKcN+knZ9CYvIlIwDfIiIgUbapA3sxvM7AEzq8zsilPafszM/mPQfp+ZtYe7VRER2axhY/L3A7gewMdOPGhmDQCfAvCr7n6vmZ0LIAiEFi6KvTIstg7E6+R7JEZd4/+eW1BRjG5CEsR/PYjJh3MNGWv8rRu0kXXy8ebaUXydXCtaC5+58TaLr2udvGzWUIO8uz8IAHb6X+qrAXzD3e8dfN/Tw1xHRETyjComfykAN7M7zexuM7tpRNcREZHAhm/yZnYIwL5E083ufntw3lcDeAWAZQBfNrOvu/uXE+ffD2A/ALQxc6b3LSIiZ2DDQd7dr8o476MAvuLuRwDAzL4E4OUAThvk3f0AgAMAMG97FHAUEdlCo0qGuhPATWY2A6AD4LUAbt2o0wJ+cOSQf/aREd3TyU7/5+Q8AEcAAFHeSDAvV5jnnoccp2dyMj2Pk43zebyINZjnrPw43tnsOgAfBbAXwFEA97j7NYO2twJ4L9aH0y+5+0TH5c3sLne/YuPv3Bn0PE6nZ3IyPY+TTerzGHZ1zUEAB0nbp7C+jFJERMZEGa8iIgXTIP+cA+O+gQmj53E6PZOT6XmcbCKfx1AxeRERmWx6kxcRKZgGeQBm9kYz+5aZPWRm7xn3/YybmT08KCp3j5ndNe77GQczu83MnjSz+084tsfM/snMvj347/PGeY/biTyPW8zsscHn5B4z+4Vx3uN2MrOLzeyfzeybgyKM7xwcn7jPyI4f5M2sDuDPAPw8gMsA/IqZXTbeu5oIr3f3yydxSdg2+QSAN55y7D0AvuzuL8V6Yt9OeiH4BE5/HgBw6+Bzcrm7f2mb72mcegDe7e6XAXgVgBsH48bEfUZ2/CAP4EoAD7n7/7p7B8CnAVw75nuSMXP3rwB45pTD1wL45ODPnwTwpu28p3Eiz2PHcvfD7n734M8LAB4EcCEm8DOiQX79F/O9E75+dHBsJ3MA/2hmXx/UFpJ1z3f3w4M/fx/A88d5MxPiHWb2jUE4Z+yhiXEwsxcD+AkAX8UEfkY0yEvKq9395VgPYd1oZj8z7huaNL6+LG2nL037cwA/BOByAIcB/NFY72YMzGwOwOcAvMvdj53YNimfEQ3ywGMALj7h64sGx3Ysd39s8N8nsZ7RfOV472hiPGFmFwDA4L9Pjvl+xsrdn3D3vq/vRv+X2GGfEzNrYn2A/xt3//zg8MR9RjTIA/8F4KVmdomZtQC8BcAXx3xPY2Nms2a26/ifsb4BzP1xrx3jiwDeNvjz2wCwUts7wvHBbOA67KDPia3vlPRxAA+6+0dOaJq4z4iSoQAMln79MYA6gNvc/UPjvaPxMbOX4Ll6RA0Af7sTn4eZ/R2A12G9suATAH4PwBcAfAbACwE8AuDN7r4jJiPJ83gd1kM1DuBhAL9xQjy6aGb2agD/CuA+PFe39n1Yj8tP1GdEg7yISMEUrhERKZgGeRGRgmmQFxEpmAZ5EZGCaZAXESmYBnkRkYJpkBcRKZgGeRGRgv0fS9OPFQvVEhIAAAAASUVORK5CYII=\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"def getRotMatrix(a):\n",
" return R.from_euler(\"Z\", np.array([a])).as_matrix()[0, :-1, :-1]\n",
"\n",
"sigmas = np.array([4, 2])\n",
"r = getRotMatrix(np.pi/180 * 45)\n",
"c = r @ np.diag(sigmas**2) @ r.T\n",
"m = [10, -10]\n",
"\n",
"d = scipy.stats.multivariate_normal(mean=m, cov=c)\n",
"s = d.rvs(1000)\n",
"\n",
"plotFunc(d.pdf, (m[0] - 3*sigmas[0], m[0] + 3*sigmas[0]), (m[1] - 3*sigmas[1], m[1] + 3*sigmas[1]), plt, zorder=-1)"
]
},
{
"cell_type": "markdown",
"id": "encouraging-paper",
"metadata": {},
"source": [
"# Defining functions to plot fitted gaussian under seaborn KDE plot"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "fifty-chassis",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"<seaborn.axisgrid.JointGrid at 0x7fd48c5fb730>"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAa4AAAGoCAYAAAAerAGHAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Z1A+gAAAACXBIWXMAAAsTAAALEwEAmpwYAACAC0lEQVR4nO29ebhsV1nn/33X3lW1q850h9zkJnfKnJCJhIRMIIJEiIoigoAKYtMKtCLabXerrd1tt2IPv8epu+1WUJS2RVQwgATUgAMigRBCyAiZyXyTO517zqlp773W74+11t5r79o1nFNVp06dej/PU6lz9riqTm699X3X930XKaXAMAzDMNOCmPQAGIZhGGY9cOBiGIZhpgoOXAzDMMxUwYGLYRiGmSo4cDEMwzBThT/pAeRgiyPDMLMMTXoA0wArLoZhGGaq4MDFMAzDTBVbLVXIMEMjpcJXHj+OW+47jDu+eRzPLDdRb0eoljycthTgwr0LuPqsXXjJOafg1MVg0sNlGGad0BbrnLGlBsNMF2EscdMdT+F//t2DEER48Zm7cN6p8zh1IUBQEmhFEkdWW3j8WB0PHF7FvU8v4+CuGr77hWfguy49HQd21Sb9EhiG57gGgAMXsy346uPH8bMfvQsVX+B1V+zHhXsXQNT7MyCSEvc/s4LbHj2KLz92HId21/DaF56B77zsdJy+VN2kkTNMBg5cA8CBi5lqYqnw23/7EP7gC4/ih645hOvP2d03YBURSYl7nlrGbY8ew1e+eRwHdtVwwwtOw8sv2INL9y3B93g6mNkUOHANAAcuZmpZaYZ494e+iiOrLfz4y8/FrrnySK4bSYmvP7OCrz15Avc8vYznV1p44f4duPrMXbj84A5ctn/HyO7FMDk4cA0ABy5mKnn6RANv/f0v4dw983jLdYfgi/EpopPNEA8cXsFDh1fx6NE1PPzcKpaqJVy6fwkvOrgTVx7aiUv3L6Hie2MbAzMzcOAaAA5czNTx0HOreMvvfQk3vOBUfOelp28oNTgMUik8u9zEw8+v4tEja3jouVU8ebyBy/Yv4RUXnopXX7wXZ50yt6ljYrYNHLgGgAMXM1Xc/8xJvPX3v4Tvv+oAXnbenkkPJ6HejnDfMydx15PLuP2xY9i7FOAHrj6I112xDwtBadLDY6YHDlwDwIGLmRrueWoZP/yB2/CWaw7hunN2T3o4XZFS4Z6nl/H3DzyPe59axg9ecxDv+tZzsKPG82JMXzhwDQAHLmYquPfpZbz192/DD193CNectXWDVp7nV5r4y7uewW2PHsNPftu5+JHrz2SHItMLDlwDwIGL2fJ8/dmT+MH3fwk/fO0hXHP29AQtl6dPNPDBWx9DO5L4rTdfgQv2Lkx6SMzWhAPXAHDgYrY0Dz23ije/71b8wNUHcf05p0x6OEOhlMLfP/A8/vTLT+Bfv+p8vOXaQ5tuLGG2PPw/xABw4GK2LI8dWcObfvdWvP7K/fiWLWTEGJZnlhv4rc8+iCsP7sR/ef2lbKNnXDhwDQAHLmZL8sSxOr7/d27Fa154Ol554WmTHs7IaYYxfvdzDyOKFX7vbVexcYOxcOAaAA5czJbj6RMNfP/v3IpXXXwaXnXR3kkPZ2xIqfCh2x7HNw6v4EM/eg13qmcADlwDwfYmZkvx7HITb3rfrfi2C0/d1kELAIQg/NA1B/Gigzvw+t/5Ap5Zbkx6SAwzFbDiYrYMh0828cbfvRUvOWc3vvuF+yY9nE3l5ruexj888Dz+7F3XcWf62YYV1wCw4mK2BDZoXXf27AUtAPiuy87Ay87fgze/74t4fqU16eEwzJaGAxczcZ5ZbuANv/MFXH/Obrz28tkLWpbXXHYGrj5zF37w/V/EiXp70sNhmC0LBy5molj34LeetwffM4NKK8/rrtiHC09fwNs+cBvWWtGkh8MwWxIOXMzEeOT5Vbzhd76AG15wGr7rsjMmPZwtARHhB158ELvny/jRD96OVhRPekgMs+XgwMVMhHufXsYbf/dWvPbyfXj1xdvbPbheiAhvf8nZUFB4z598FbFkzxLDuHDgYjad2x49hrf83pfwQ9ccwisuOHXSw9mSeILw4y8/F88uN/GzH70LW8z9yzAThQMXs6ncct9h/Nj/vR3v+tZzcO2UNszdLEqewE/fcD7ueWoZ/+kv7+PgxTAGruNiNo3/98XH8Ou3PIh/9e3n45w985MeztSw2orwXz51P1598V782xsv4Ma82xv+4w4ABy5m7Eip8F//6uv45Neexr+98UKcxq2N1s3JZohf/dT9+O7LzsDPvOp8Dl7bF/7DDgAHLmas1NsRfurDd+Kp4w389A3n8TL2Q7DcCPFfP30/brxkL372xgs5eG1P+I86ABy4mLHx5PE6fvSDt2PvUoC3v+QslHjl36E52Qzx3//q67j27N345ddeAiH4c26bwX/QAeDAxYyFLzx8BD/5oa/iOy89Hd9xyV5WByOk3o7w67c8gIO7aviNN12OoMTreW0j+B/KAHDgYkaKlAr/5x8exu/94yP4Fy8/F5fuW5r0kLYl7Ujidz/3MJphjPf/8FXYPV+Z9JCY0cCBawA4cDEj48hqC//qT+/EcystvPsV5/KH6ZiRSuEjX3kSX3r0KN7/w1fh4jP4S8I2gAPXAHDgYkbC33/jOfybP78LLzn3FLz+yn3wBc9nbRa3PnwEH7z1m/jF73oBvv+qA5MeDjMcHLgGgAMXMxT1doT33nw//ua+w3jny87mb/0T4snjdfyPzz6IFx3aife+7lLMV/xJD4nZGBy4BoADF7Nhbn34KP7tR76Gc/bM4y3XHsIcf1hOlGYY4/998Zt44LkV/Nabr8CLDu6c9JCY9cOBawA4cDHr5mQzxK/efD8+c/9h/Mj1Z+HKQ/wBuZW47dFj+MMvPIofuuYQ3vPK81D2OW07RXDgGgAOXMzAKKXw6XuexS994l68cP8S3vTig6yytijH6238/ucfxVorwm+++XJcuHdx0kNiBoMD1wBw4GIG4oljdfzix+7BI8+v4u0vOQsXns4fhFsdpRT+/oHn8adffgLveNnZeNe3ngOPC5a3OvwHGgAOXExP2pHE+z73MN7/j4/iOy7Zi++69HT43AFjqnh+pYX3/+MjEAT85puvwFmnzE16SEx3OHANAAcupiuff/AIfvFjd+OU+Qreeu0hnMrNcacWqRRuue8wbvrqU/jpG87D2647k9tFbU34jzIAHLiYDp4+0cB//sv78NUnjuMt1x7CVYd2TXpIzIh45kQDv/u5R7CjVsKvvfGF2L+zNukhMVk4cA0ABy4moRnGeP/nHsH7//ERfPtFp+G7X3gGKj73wdtuxFLh5rufxqfvfhY/9x0X4k0vPsC9JLcO/IcYAA5cDJRS+Kt7nsWv3Hw/9u+s4gevPshpwRng8WN1vO9zD+OMpSr+6xsuw74d1UkPieHANRAcuGac2x87hl+5+X4sN9r4wasP4RJuijtTRFLik197Bn9977P4l99+Ht5y7ZnsPJws/OYPAAeuGeWuJ0/g1/7mAdz/zEl834v24VvO3cOT9TPMk8fr+IN/egy+R/jV113KX2AmB/8jHAAOXDOEUgpffOQYfvvvHsLXnz2J11x2Br7twlN5gUcGgHYe/sMDz+PPb38Cr754L37mVRdgzwJ3+N9kOHANAAeuGWCtFeEvv/Y0/uALj2GtFeHGS/biW87dw62AmEJWWxE+dudT+PyDR/DPrj8TP/qys7lp7+bBgWsAOHBtU2Kp8MVHjuIv7ngSf3PvYbzgjEW84oJTcdn+JQh2kDEDcPhkEx+940nc89Qy3v7Ss/DD152JpWpp0sPa7vA/zgHgwLWNaIYxbn3kKP7qnmdxy32HsWuujGvP3oXrzzkFO2vlSQ+PmVKeOt7AJ+56Gnc+fhxvuHI/fuT6s3BwN9d/jQkOXAPAgWuKUUrhsaN1fP7B5/HZ+5/Dl795DId2zeGKgztw1aFd2LvElnZmdBxZbeGW+w7j7x94Di/cvwM/dM0hfNuFp3LKebRw4BoADlxThFIKjx5Zw+2PHccXHj6CWx85ilgqXLpvCZfsW8Kl+5awEHAqhxkvrSjGlx45hs89+DyePtHAd156Ol57+T5ceWgnW+mHh9/AAeDAtUVRSuHZk03c/8xJ3P3kMu54/ATuevIESp7ABXsXcP5pC7j49EXsXQq46wEzMQ6fbOLWh4/iS48dxclGhFdeeCpe+YLTcP25u7HIX6I2Av9jHgAOXBNESoWja208faKBp0408M2jdTz8/Kp+PLcKTxDO3D2Hg7tqOGvPHM7dM4/d82xPZrYmzy43ccfjx3H3U8v4xrMrOGvPHK45axcuP7ADF5+xhDN313hlgf5w4BoADlxDoJRCK5JohjHqbfuIsNaKsdqKsNoKsdKMcLIR4ng9xNG1Fo6thTi62sLR1TaOrrVQLXvYM1/BKfMV7Fmo4LTFAKcvBdi/s8YOLmZqaUcSjxxZxQOHV/DokTq+eXQNx9ba2L+zikO753BgZxWn76jilPkKds+VsVQrYTEoYSHwUSt7qJX9WU07zuSLXi9TE7j+4YHn8Wt/8w3EUsEOWUEHDwBQShdQSqWSn2OlIKVuaxNJhThWiKRCGJvf5fhfuyBgruxjruJjPvCxUPGxWC1hR62EHdUSF/8yM0Mrknj2ZBOHl5s4stbC8XqIRjue9LD6QgT4guALAd8j+ILgCWGe9UMQIARBkP6ZQCBCksYncx0A2DVXxq987yU4tLtwXTQOXAOwpQIXEf0VgFOK9nlzO0/1FnYfGPlN9Rugw52CMj8qAOmz3iaVUsnPUErCnCPba4EozzVHPrYtxCy8RmA2XudUvEYCABIgIhAJgAgEYcMBCGQPMk8dH/iyseKJ6kLsXjF3E3K2JtcjHW1MxBnPBHJ47OlvqHZ9tWDXEaXUjeO453ZiSwWuaYWIbldKXTXpcYyTWXiNwGy8zll4jcDsvM5ZhPNUDMMwzFTBgYthGIaZKjhwjYb3TXoAm8AsvEZgNl7nLLxGYHZe58zBc1wMwzDMVMGKi2EYhpkqOHAxDMMwUwUHLoZhGGaq4MDFMAzDTBUcuBiGYZipYksFrhtvvNG0WOIHP/jBj5l8DMwMfF52ZUsFriNHjkx6CAzDMFPBLH9ebqnAxTAMwzD94MDFMAzDTBUcuBiGYZipggMXwzAMM1Vw4GIYhmGmCg5cDMMwU8iR1RY+8PlHJz2MicCBi2EYZgp5ZrmJ//Lp+yc9jInAgYthGGZKmdVVqThwMQzDTCkzGrc4cDEMwzDTBQcuhmGYKWVWV7DnwMUwDDOlzGbY4sDFMAwztcyo4OLAxTAMw0wXHLgYhmGYqYIDF8MwDDNVcOBiGIaZQmjSA5ggHLgYhmGYqYIDF8MwzBRCNLuaiwMXwzAMM1Vw4GIYhplCZlhwceBiGIaZRmY4bnHgYhiGmUZ4jothGIaZKmY3bHHgYhiGmUpmWHDBn/QApoFvF98/6SEwzOSgPt9vldyccYyYW+SfT3oIQ0EzrLlYcTEMw0whrLgYhpkN+qmncV5zSpXZVoUDF8Mw24NxBKZR0WtsHNTWDacKxwgR3UhE3yCih4jo58Z9P4ZhmFlglhXXWAMXEXkAfhvAdwC4CMAPENFF47wnw8wMJDof08p2eA2bzAzHrbErrqsBPKSUekQp1QbwYQCvHfM9GYZhtj1cgDw+9gF4wvn9SbMtgYjeQUS3E9Htzz///JiHwzBTynZSV/2Yhde4QdzPy+bK8UkPZ2JM/P8MpdT7lFJXKaWu2rNnz6SHwzAMs2VxPy+rizsnPZyJMW5X4VMADji/7zfbGIbpxiYpDRLDp5qUVCMYSQ/y7wW7DxNmN1E4/sD1ZQDnEdFZ0AHrzQB+cMz3ZJjpYQxBahQBaVT3GnlgI8HBixlv4FJKRUT0bgB/DcAD8AGl1L3jvCfDMMxMMMOSa+wFyEqpTwH41LjvwzBTw4hU1mYqq41SNMahVZh9/2ZceW39v/74mLg5g2EYhlk/Y55d3NJwyyeGGTcjUFgjV1cbHdMIVI77WoZSX6y8ZhZWXAzDMNPIDEsuVlwMMw6GUFlDqatxW+kHuf46FNBI1Bc7DWcODlwMMyo2EDQ2FKQ2ydwxVCDJXmhd49nQfWcwbTjDgotThQzDMMx0wYqLYTbKBpXPQCprndcehzV+0Gv2VUjrVGBDK68ZUV2suBiGYZjpYoYjFysuhlkP45rHGuC6G1ZVw86JDaiQsqf0+FR1x9Pj2iRo/L0QmamEFRfDMAwzVbDiYph+rEOxDKyK+lxzHPNgG6bffQpU08A29z5uwA3Nd82Kw3CGez5x4GKYIjY5WHW9xqDj2Ky+hUUBpI/5Iv/aCoNQH1PFUIaNbcoMxy1OFTIMw0wjsxy4WHExjMuACqevyhq1uhpY1Y3h40w5KqdoHHkV1Md80dV0MUCKb12GjW2eMhzHn3paYMXFMAwzhdAMRy5WXAwzqvmsLtdZ9zldFVn/DyoSo/kuqqSjUoru20uFuYqoi/rqOWc1DsPGNkTMcOBixcUwDDOF2LgVxdszFdoLVlzM7DLsfNZ657EGVVddvkn3VFMj/vbd7V6JEuulwtzXVKS+1qO8RjE/tU3bQNl3ea0dY6k6WxqEAxczW2y1YFUQAAqDRrfANKLUYAeyW5qu836FwcwNYkXmjQLL/HrThpwy1JxshFiqliY9jE1ltsI0wzDMNmO5EU56CJsOKy5mNhiFzT13jYEUVp9U4MDqaj0qbBQU3U/Kznsq1fEalHucUsXmjfWmDZmuHFtrT3oImw4rLoZhmCnm6Fpr0kPYdFhxMduXMcxn9VVZ+f395rDc/f1U1QBFyLQOa79SEsHOxY7tzeMngaIprm4qLL2gHoMQ2XkvlVNQ7rzXoMqrz1zXLKu0505y4GIYZpsT7FwAAKh8QEn2LwJKoXns5GYOi9kAtbKHZ5Ybkx7GpsOBi9l+DNN5fViVlVNYA6urPsqqUEmto7GuDVaZa1pDdaGbjxDsXtLqy+KKncQ1KFLV5cxr2dfdMd/VjXFb1rdp+6fd8xU8eZwDF8NMJyNOC240WHWYLYqCVZcgRb1Sjvnzem1zKApYALJBxCuwsQMgCVR3LaF5fEXvEs6Hvvv5n3/NUnZPG9p7iB5B096/KAW4TWuyNsqeGQ1cbM5gmG1IsHOhe9Bitg2nLurA1S3tu11hxcVMN+NMC/aytndTWXZ74bYu6qrHtfuaNwqOKwxYvT7YuhQbJ19rzdhIpvcvVF9F6cN82tCONd9lo4tFflBm1Uo/X/FBBBxZbWPPQmXSw9k0WHExzDZhXCqLldvWhQDs31HFw8+vTnoomworLmZ62YjaWo/K6mFtJyE655cKt1Gxusof101ZdVNfhmDHQtd9CUp12W/Uied1KrK8qYIo8zXXqi+lpKOanPvZ15Cb7wLWYdjA7Cqp9bBvZxUPHF7BtWfvnvRQNg1WXAwzhQQ7FtKgtRn3Y9W1ZTmws4Z7npqt0gVWXMz0MIxzcCPzWf2s7T1cgonKKpq7cs919/VSV+Y6HcFqEEe86rZqcU5R6V/0k6vC7NyUEM7clLl90byXRNf5LiDnNLT0KUrOsE2t7RvlzFPm8OHbHp/0MDYVDlzMdDCBtGChtb1fsBKU3d/NDl+4v2AbCQQ75rPjWEf9FoAudVpux4uCZUiUSoOi5wGArus6upwbS3qODdZKyE7Dhj5gsG3Mujhz9xwePbqGZhgjKHmTHs6mwKlChtmCBDvmk8dWglOGW4+yL3BwVw1fe+LEpIeyabDiYrY241Ra/Tpe5K3tOXXUkQ50ldR61ZWrrPqZN9aLQKeiyassm3ZL3hPKqi+zCUIg2L0ESJkUJidffx3llUkbukYNe/W8Rd69t5syNOMamUljm6YZzz9tAV969BiumRGDBisuhpkgW1VZDQIXOW8dXrB3EZ9/8PlJD2PTYMXFbD02WlTcq21T0XxWURFxnwLijMrKKylXpWWMGNltwc7F3nNhKFBmHccMQNG8UaJm3HkmlSqwRGUV2dydoTkqKtilO8wnTXm7qaICi3w6BhrJPNesdop/wemL+F9/9xDWWhHmKtv/Y50VF8OMGWtdn1ZlNSiswCZHtezhvNPm8U8PHZn0UDaF7R+amelhk+ezerZqcpRSV5WVV2lEgKe3JR/gXZ2GqYrLj0Fl5sLcF9bxUotJBEdOMQEgys1bAabdUk6JKdeebi4QO9d2lVfO2h7sXkTzqFZfXR2G9jJFc12Dsk3nqzbKC/fvwN/cexivunjvpIcydjhwMZNngPqsfkFr6LRgrp9gV2t7JrBlg2Kwc6EzUAqBjiBFlAYnO9RMppCQj1JqwDQhuTZ2i1BJ0FCJKYOSgEKU7s8EocR0YX7wZJpSjO3JKDRfBLsXtXkjLhhkQW1Xdn8ubckMxIvP3In/8PF7EcUSvre9k2nb+9UxzJhJ04CcIsvDacPNZc9CgFMXK/jCw0cnPZSxw4qLmSxDpgcH7X6x7rSg18UOb34Odi12FOh2s8MXqitylA+Mokr2b9yQkVneQjl5PSuQkm1pelApAllllKgnJ31oT5YCjkxLLtMrhdlhjc+Mq+Pgzv2cDlwX15y1Gzd99Um87Pw9kx7KWGHFxTDrgNXV+mDVtblcf85ufOa+51BvR5MeylhhxcVMhh5Ka8MmjF5tm9Y7n+V5mbmtYMeCUWEo2N9p3lB2VWFHXankOGSurfdR9mtkh9LqJmsK1EsyX2XOVK66MsdIAiE1ZCRXkfY9Uek8VqJWZTrflelLiOxxGaPFaOepZtXuPig7amVcsHcBn777Wbz+yv2THs7YYMXFMD0IdrLCYqaLl523B3/8pW9OehhjhRUXs7kMobQ69ufdgwO2bSI7J1XUqsk2lN25mFjbu6os24ZIOErKbZkEo7yo4LjcHFd23gtQA3+lNAopI3KUuwsKqlOFCZVa0CWBRCLP9Dkxgezrj515LTuuRGUREKt0P6Df6zhO9mea7zpjcM8p7BjPbIgrDu3AB299DA8cXsH5p23PL10cuJjNo0vQGiQ1OIjdfd1pwSSQiExg02lBSu+TBLBcsLKNuN1UoEfZbU7gcp+TwGSDSy5w2e1+rZZ/ZxDV6+kvtuzKjkUpkKLMPpLkBLM0ZUjJ2iQKKhbmtiZ4eAoqtkHF2tNdc4aTPqTs/bTxo8Bo4VKYVtwEZsDk4QuBV1ywB3/4T4/hV7/v0kkPZyxw4GIYQ7BjYd0uvlHh12qFgau52ug4Nph3gpkTFzIBjZlpXnHhafi5v7gLP3vjhViqlSY9nJHDgYsZL+Owu3frgtHL5l5kWbdpwV1LnSqswA6fVVnCMWC4KqxAXZlz0tpfghIEv1Z10nlAo940+533ouSkEk2Qqrea6X6zrVoL4C3UAAXEa/XMPiWUVl1A8gypoHxzQKwSpZVRXp5RZ3C22fNdFWZTgPmUYf5nZtPYNVfGFQd34kO3fRP/4uXnTno4I4cDFzOzTMKqbZWVEkBztbmxlk4FNNdMMJNAdc4oMhNHorW1jV+YmVpuvHgvfvMzD+Cfv/RslP3t5cPjwMWMh1EqLbuvX/umnOmiaD4r2L2U3idRXu656bxXYqYoUlnCVVypErRzV8kcliB4JpAooXVNs97QQaqUjkHPcbmvt0cUyxQZ559VoshIAsFcAG9hLjFvkGMlj0+uJa8vo6qglVcy3yWc/oXI/w2cFlPuvJZ97+I4/RsUtX/qhSBu+zQEZ50yhzN2VHHTV5/Em158cNLDGSkcuJhtT7BzIbtMySbh12o6mEGrq8RAsYnZs+ZaU/snTLAjZ3HIVJmpJM0Xr6xu3uCYsfPdl52O3/67h/GGKw/A67Yc0BTCgYsZLSPo8E6COq9T1L6po+mt3l/dtaS3ufNUyRxXgavQEx0OQuWl10sKh/2cyrL77RyWSOeuEnXVaOgxlF0nIcHtdavvgcHng3JLaQGOadAtCLY7BdK5KecajYY2flCskoBWWdDLrsQnV4CYMtdRRI7D0Jnr6lBSzoScq46tLX7d0ovZKC84fRFzFR9/+bWn8b1X7Jv0cEbG2AIXEf0SgB8DYJfl/HdKqU+N637MFmA9QatXJ4wN2N2rp+zo0v1iwLSgn6b7gKy1PdnnBCvpO8eaZ2/ezF8BaDR1UFAlkQlY+tpunZZTu9UnbqVrP6ZpOnIDFvR1k5out8u6E3wAQETFtWKttToqCzX4C/OIjPoiG8CEs+CkmzJMUqU5k0Y33NovZqwQEV53xT78xi0P4DWXnb5tusaP+1X8hlLqcvPgoMWMnGDXYrIC76Twa1XtEISev2rWOy3s00RrRbsS/fnOGjJm+rjkjEUsBD7+4o4nJz2UkcGpQmZ4hlBag66jlbe7B7sWASFAboFxkvpzJpO8VFXpfWlaUPnmOA/ONtf6nlVheZXlLegPdkmkg5UgyHJWhYFSlZMxbrgpQud19sTOUyVGDEpNF05Wz97Hpv/ITdeZFJ70CSJKx5MqNn1yc62BoBYgszYYqeS19B3thGzw3MewEyLC9191AL92ywP4nsv3ISh5/U/a4oxbcb2biO4iog8Q0c6iA4joHUR0OxHd/vzzzxcdwjAZJq+waqnCWtuYwqJKteuDYbrhfl6uLh8f+LzzT1vAod1z+IN/enSMo9s8hlJcRPQZAEXrRP8CgP8D4JehZ2l/GcCvAXh7/kCl1PsAvA8ArrrqKv66NE2sd+XijRYWm9+DPTuz+8hRVK4dPjFaDDif5aUtn9J9IjVdmLkusTAHCN30qNFs6pZIpdSwoUSBuvK0+qJKNTOH1WqbObCMBV4/leec4KUA1W4YlWWu3bmkVqqu4lSFpbKoc9EsgkrGTXFqMEnW6zJmGH9+DtGymetyTReuScP+nJhCKLemlmOdZ4bC/bw8cP4l63pD33TVAfzyJ+/DG686gN3zlbGMb7MYKnAppW4Y5Dgiej+ATw5zL2a2maTKsv0CJbTCApB2yOgBVao6DSkIYbORllo5qcIi2q1UwZXLVVC5qlOA7c2dO2ut1lHhea5twxk7qrj+3N3473/9Dfy311826eEMxThdhacrpZ4xv74OwD3juhczATY4rzVwYbGjlPRqw44dvpvK6jGfldjc/dQFl8xn+Y7isurJIyhPW9ulR2iuNfS2sp0DI0jHBq/P6VRX7XZDi49Sqq6UoI7A1e2rc1PqYuJyuQrytRJTJrC51nc7BhLa3g6kbkABFAfKpGM8pdLNc9STXb3Z/vliyhUjj5mC+arMHNZ6GubOQHPdQfi+K/bj33zka3jLNYdw6f6lSQ9nw4zTnPHfiehy6H8CjwF45xjvxWwWo1qWpGDxxaL6rGD3UrKN8oFJ5AKT3e6nQS/pfuGnaUFlg53TBUOafwnK3MNbqEF5hHq9ASUJKAtIj9IgJeCcY8ZfqUJ6WjElx5XTIOV2hO/43M//nm3mjmas04WVShXwjAJsNZwUYcHFzMkyIogkwtldlAY9z+kon5hFTEcMgWyqr+NPPBkTBrMx5io+3njVAfz8X9yFj7/7pVNblDy2wKWUeuu4rs3MBpNID1rThQIGNl1QWZ/TbjUSFTYuWq1GsrJIpVLNBK7NTiUy08nLzt+Df3zwCP7fFx/D264/a9LD2RBsh2cGYxxKC0C3wmKrtJJ9RXZ3d/0rz+vYr7tf9EgLJiqMID2RpgWNylIlfZwNRtKnjBEDVR2wWqEOGLKcptUy1ndHaQFG6PRIFVJ+g/k98TuYwNWMGokRo1KqAn5NpwmttV3ptKKAggrte2v3qVR9ydRUYlUa2dSsk6JVRGmnDtcWzx3gpwpBhH/2kjPx3pvvx6su3ovTl6bPybo9yqiZbcVmKy1rbx/U2k4VbZhotxoZI8UksWNxH2GjwTZ7ppD9O2u44aLT8LMfuVt/iZkyWHEx/RnA9t7t2A61NUhhsavC3L6CRXZ330+2ZQqKASjhAX5WcUEISLPGVVJA7BEarZaZx9KHST+dz7LHI6hCeTpdJz2CKqWmjOTZKhq3wDg3x1Vkgdc7nM25ru8kHeu7vY5E2icwcadTYs4QSgc0IbVhozRXBUUKqtWAcoqXIdBhulBCmX6ElKor4Q6+4MMu+duPsBchmyrGxmtfeAb+/cfvwUe/8iTecNWBSQ9nXXDgYrYMm6m0/FpNz2Otma7tPdYrsmpFAQMrLOUHyTlF5oz04u5JAGJnkcgREzYbqPiBfj1T3paKGR7fE3jnt56DX7n5flx/7ik4Y8f0qHIOXEx31tvpfQgHoW3hpDflCouBbMPcRFmlKkwJkXEO6mcvtcEnKmsO8Al1Y29HiaB8gjTnykRlAaiadbSMU1D6pF2CMPNeibKDOSdw1uEC4rCRbe/kvp0F00KeXwVKga7ZskHMcRema2o5b429dpzcNlVNdgVjhUQJtlpNlAO98nIyHpnuT9yJ5BQbTwHc6mljnLl7Dq++eC/+1Z/diQ/96LUQU+Iy5MDFdDKqjhh2X78O78bKnvTUcwJSX7u72+kinxb0vMSy7i3OAXDSgiWRBiuPdKACkm2oaqNGu2XqsMqkr2VTieSkCMtGXQkgko3kZ1QACOfjv89ngkSqgnwTxKzRgqJmElTICTL2kpkYlFvqRHqAULYGTkF5AFWrUGu6ma7rrUhNFwrBQg2ROaYrhcua9GGzgwynG3vy3S88A79y8334/c8/ih972dmTHs5AsDmDmTjBzvGmCG3ni+ZaM+180QUqD2a8UH4A5ekHoNVVHI0u/RbHjeQBAPCC9DEkYZPThEyKJwjv+tZz8L/+7iHc+/TypIczEKy4mHWx7k7v+fSgq7QABKfsyCotIJsWdIuOkxShSQ/6pA0YgF7kMae4bFpQAmi0mp1pQc8xXxhru/RNWrBEiVqz/0qkAKRVVx4QqUbyMzxk19lyjRhJUW/u2UUhlWZuWpCACA1jOwc8rwohTPCKAbSbyen20onSMipLKJVNC9ojEws8HDeIM27LoJ6LbspmANeakrLncRvumMEMxGmLAX7omoN494e+ipvf81LUyls7NLDiYibGOM0Y/pxVWY2eKsu1ivdSWVZhAYCMRquu1kMc63vb+7uqj2GG4VvO24ODu6r4Dx+/d9JD6cvWDqvM5tNlfmtDRgygc17LUVp6E2WVlr2G3ebMZ7lKC8jb3T2oEmm3oKdXnaq3mokpI1FZvkjms5RNIcYNSCKgTFC2jVNJKyzlBcm2SDUSdSU9x95uu2U4tvK0LyE6FVaR78Hp9A6ZXsCKC4pTgSTNc+QUIFNJB1/RbqZGC1dlJePRP7RaDQTVKlSzAUVw1jVT6TmDsI75KlWklIpUFhstJsbbrjsL//7j9+BjX30S33vF/kkPpyscuJhNZ1xKK5nLWjVqqFSUk0vt7W0711PqPEb5QRJ8EnW1hdffk2aM5GvHIEXjs9Uz25dq2cNPftu5+KVP3IdL9+/AOXvmJz2kQjhwMZpBWzqts9N7x7zWIA5C39kmUut7orR829JJKy5/Tqss6ROaq83Usp5zDlKlClkmKKHTgnHFqjAkqiquFMxh2TkuL1tsbO3mcJSXsoqFAKhyl/5NmTfNnJyKD1JhMqeTsafn5pi0W9D8bDdGDfiiClUOoExneSUIysg5d8ms5PYF8T2YqwLRGJWPHEDSDTiX1dUKz3NhG+LQ7jm84cr9eNcffQWfePdLUS1vvW9sHLhmnQ0ErK77+nR6hydQ3bWU/KyfvWzAstuETQta8wWl9VmlNH0oluaM+aJlurhT0ulC+oTYdrcwc16tqIGY7HEw9zDBywsAXzsEpZf2KEyCle/URvlpkFLJSy3lZo0VJDW6mzGQ80TYHKCsOkYOAqitg0i+J6BKx2ZThkpohShK6dIqyhmjvW6pVoWKVLbDB7JD7WuFT8ahOn+324pSgVPYYmjW+LYLT8U3Dq/gFz92N37tjZdPejgdsDmD2TRGbXt3DRjdcM0X3Wzg1nghowbicB2mC1UGVBlkHgpNKDQyj42gVBOwD3sflNd/Hb+3aUNtQp/F5vGVsd+DGT1EhLe/5Czc9ugxfPi2xyc9nA5Ycc0ygxoxuu0bsNM7PIFg52J3I4af2ya8jNICjPnCOc6fq0H6hHqrCZR0GjFVWkaZzadFxFIAKGsVZpVWkhb0teUcZaSFyF6aIlS+UVZeKVUunoISWpEkistz1qtKuqx3eSMT63vquEjSg0JBxXaCzQaXKpLJOC805zjCxmZzlL5kbFKG2jCCtGN8L0v+sHQTUkOYLTrSgJz+2zSCkoefuuF8/PIn78Ml+5Zwyb6ts/AkKy5m7IxSaflzNfhzNTRX612VVt7inkf5QaJGZN/C4bLzQKKEFDbX/KDQ1EosGRPDjJ99O6p423Vn4p1/9BWcqLcnPZwEVlyzyBBra2W29TNioKDAeAAjhvK9jNLSzwL+4gIAJCpLmoJiQBsxYj81YITNBmIJSMeAAQCyGiD2tAtPmrZMspQqlnTeq2SUloQUJmD4KpnjIk8lzggh0gmr5K1KFFfB8vOKOuzuUKT99wAQp375ZO4pgrH5NyFigEQA5YWp4rKqL/kPuhov+hHMBVljhityzA0pmadSxYqqyHxRcJyyx6ku1+kBmzI2h+vO2Y2Hj6ziPX/yVfzhP7t6S/QzZMXFTAXJfNZqd8NAr7ksq7LidiOxjndeIKuulHJVzvBQVOr6WC/aMTg+5RXVBzRmMDPBm198AEfX2vj1Wx6Y9FAAsOKaPdYzr1VkfV+Hg7BwXmtQB6FVWiU9n6V8gUYjVVp6n0BcIlBgVi42S4Ikc1g+Ia6l81gybCAuO9b3RF2VIUu655L0TaASCqpk57aMyvAVhCeT/fZnq7hIKAg4c1thCa7kURUzJ6ZIqy5AOyEByChItkkvRvKdMrJjSB2IdhVmAqBEGVDt4oJno7i0wLMbi1VKgTDspJ+QkY4KgzFmKNVZeNynvRMrpq2HLwTe/Ypz8R8+fg9eeGAHvv2i0yY7nonenZk4g6QGSVCnEQPoasQAkLW9b8SIYSzv3tI8FIBGowlZtvv1NeISQc2ZdbVkE7ERIDawxdVA29qjhnaal6EDl00LlsvmeAUlTM2WDVa+ApVsKlB/kJIv4Zlg5QmZxGi7X8QlCHLShrVVfS0bpOwzAGkCVhQL89Y0km1oVaEgTCrQ7IcyKcQ0LajiBojMGkpOsFIFf1KXcqUKVa8XZxGlHiA5vRPJDTJJbHLs7oMGmmFs8BzMJs6OWhnveeV5+Dcf+Ro+8q7rce6pkytO5lQhM3JGZcbolx6kINcBwyExXxSlBals0oI65TZUOjAsAWFZPwDAawJeE+QNcU2vpYcZl0DxAN8tic0azOZw7qkLeNNVB/DPP/hlLDfCiY2DFdesMMAaW/njOtKDmeMIeSMGiHoXGLvqq5cRo+RllBY8AVkmyJK+jixpE0ZcIrSiJlAmxGWjwPwAqgw9l0WAtCrMPpu0oFJNqCRVKAFjeUfZqCxPgnz9Ld8v6ZYVHin4Rq55HkGQgqitwLOKS6SqQDi5N9ulXRrzRSwFYvPe2jRjFCvEcfo+xyaYkgpA8CApBmzKMrG2k3Y5UrAhEwa5mT2jhoL5AFG9YXwn2SJiUipNF8rsvszPSqF57GT2uG4M2hF+gO1mZ+/7MUPz8gtOxWNH1/CTH7oDf/DProY3AbMGKy5my9FXaXUxYViVFbd7qKx1Gi4oKoFC/UBYAvlN/fCaEMOoqgFR1Br7PRhmvbzl2kM4Vm/jv336/oncnxXXdmcUPQiBYut7bo6resqO/p3e12PEMEoLMH0HHaUVNhuIS4S47MxnlbQBIykiLgFxMo8FQDYhyyotKDbqCiUFGHXlmU66nq/gCwJqqyiZea2SF8P3tPoqiRiekSyJ4uricIgaekDSqItYSYTSA1UiRLG+n/uOK0WJISSWACRAyke6BDKSZ7veV3LnLl9+7TxVuVKFaje6mjFI5ea0MnSzn+d+z5svOoqIVWqDd49htTQ1+ELgJ7/tPPzSJ+7FBXsX8forN7eTPAcuZiSMouN7Rml5nQG3q9IqpcXE2RNsfrDZ0xGn55HMh6sv9RyTiAF/vWvS29t12tu9SgvKpAqFFNqA0fKB2AbwHikz0YTunMEwW4fFoIR/ecP5+OWb78OZp9Rw5aFdm3ZvDlzbmUGs7xttnut5WaWVOy4zr+X3dxBapdVcrUPllBYAqPka4nKqtAAgLhPiqg5aETUSe7ssAaCyVlwmnScrqVsQJQmKfZBvi4lj+IFOS3omWJV9iZJVV34EAKiION0mYvhWaQGIm46jsNSGVzad2RUhMjbAyLx3oRQg0vdrxz5kWIYKPaiSvp6UBCns+23fOwXEHpQXIZFVPdo3Cb/ac8ViMu5BIOsgtCsx6/25eSylsoXHgFZKRTZ4u8VVUYN0hAd4teMp4sCuGt7xLWfjnX/0Fdz04y/BgV21TbkvB67tyHo7Y+S3D9KDkAgQpB2E9jhnscdEMbkpRU+kAct2evdF8uhpxCibOi0nPRhVA8SmkW1cMgYMKusAJpuQJZWYMmxQgPBAPgG+gjDByvdjlEsmOPlpsAo8vc0+l7wYJZOu80kCbf1aBRTgx6gEqeKTJppE0kOk0p/1uV42rehFaDbmELc9oBJDCJXsl/ZZtEAwTXN7zYU7l6XQGjyy25Pg5FrbMz+bR2LEcM/LBzMnsLj7MsEnqyaVm0rkbhlTzxUHd+K7Vpp42x/cho/9xEuwGKy/oH69sDmD2TCjtL2v24hRCnqnBvPXkT7I9H1SomnSb+snbpYQN0uImvpefqWJUqAfQ1EyJozW+NY+Klc43ciMh1dftBfnnTqPd/3RVxDG4/8SwYprRhh1D0IIArmqCsha3+0230tUmCpYU8tfXEC9lxGj3GnEiKpBYsJIumQE5cSAEZdT84UqKz2HVVJQogUqycTeXiprJVX2YwQlXZNi1VXVD1ExaUGvre9bVjFKpTZKlQZKJBNVVGh9t89EaJtWHe3keJnuN/eIlYAMmlDtCmKhECcKF9nnInIihFRn+i/Z1mzohhqOeiKpfw/mAlCoknPJpPaS9KBEpw3eVWEoWMakh/LKvoaCDztWUVMDEeGt156J37jlG/iFm+7Gf3v9ZWmP0jHAiovZEKNQW9aMUcS6lJYoVlokfVDsQ1EbSqzPVi6NsoqN0aJklVWl/xpWYd3PPKLkMZiaovbkVpzlHoXMRvEE4SdecR5uf+w4fvvvHhrrvVhxbSdGUWTcraVTUQ/CpOu729IpN8fleUkHd7eVk9vpPek9mKyjVctY3gFtxHCVFmDs7oGxu3tNxNaAUVIg6UOVJVRFBzNRTouJK2Wtrux8VrXURjnS96kJvW1ufhWBWfcqEPq5LKJ0jkvIpC8hALTWjOUdhLIXwwv0OZHy0DJzW/V6DWgBUeyhXE6VFgB4VsFVmkAjDeh9v7RSAKCdCjIFCK+azGPZbekj3UjJ3JUCxSpj2ICSmWMBq7xyc1My3dY8upxcLzFlKLUxU8Y69rEy2zpUyx5+5lUX4D/95b3Yt6OK171oPDZ5DlzMuhil0mqu1oFSLtgGet/ASss01kWyhL2ffKgq6rN+UMvMeUkJeDG8ShO+N3gbGxusLJVqE8qEkFB1qiY/MONZKyFuePCqG7PbDwLF4y+OZpgids2V8a9fdQH+8yfvw6mLAV5y7ikjvwcHru3CoF3f+xUaA1nnIKDVFhGCXYvFa2vZ7/terhM8oNeQskXHngd/QTfmzKypVUqb5lIJaMbNjNJSpQChcpRWGYAoI/aaiWswLiut7jwFaY0SZQkvUVrGIQiFwLY3qrTglVuoehHmjDmi5qgsq7gqwsyFiRAlitFcK8MDUPGaqNYcJ6ESiK2b0ASutvRBlA1wtVoD8RrB82K0jRrTRczKfTe7k5u7crdRxgFothl7u33Y4/LnQgEUu/NZyQsz+1WnmzCzrcgi7wy7aO0tVkvbkgO7anj3K87FT3zoDvzJj12LF5w+usVkAQ5c0884re9OqpA8T2/vt0QJpTb35BrmZ39pHtInNNcaaX2WT0mQUnM1tJoNyHK2IwagO0QkRoyKVlqyjMSIAV/XQSnRBkwazqtIlGzAssGqHGJxThsI5kxgqpVamBNaDdV8/VwVISo2RUgRmmtlVCiGIIWl+eMQUBAmCkiT7pNIjRiRMmNAaitXgpJ9lWoLUaMM3wRWIRV6tnwrypQ5KUByPv+tESMJUipNCybNNxSS4KEDmgJJ5QQ7lTVlwBzvmjIsUiZLmAAFtVsDdIXfUO0WB70tzUVnLOGt1x7C2z5wG276iZdg347RuVrZnMEMRLBzYSTXaa4Vmxtsp3eXwo4YbnrQnpvY3DtTg9T2ErMDldsQpfUtP95aK6O5pmVdtdZArTad5oVStQrV6m8sYZhRcv05p+DVF+/FW3/vSzhRX9+/vV6w4tqGjMz6bp6D3UsdxyG/QCSQ7YwhUtu78nVnjEazmaQFrWFDlkRmTa2k07vTEQNlXWAcB0ZpVfQt4rLSSgsKcdAEGZu7H5juF0qhWo5A5RZqJauuQsyXdOCbN+qq5rWxYBaQrJrgR3VCWcSoeA3smNOdzksUQ8D2NEy/7cfmvY3gwzdKq2XloUjnu3xKVZggCRpg9caMWFF2VUjzJ5DkWNftNiR29+w2nZ6jxAJvFBZMejCfApQy+TnTLSNjyjDbkoLlYut7JkWYObc/PU0ZzFTwnZeejuP1Nt7+h1/Gh37sWgSl4V2zrLiYnoyyB2ERRWtqrUdpKWp3mDAo9EChVVmD2+BbayW01sporZUTddVPYa2tllFPHiU0co9elIIWZGNAi7yN1u42EXRsE/7mFxl31G4xTI4fuPogqmUf7/mTryIewZcRVlzTygas7/rX3GrG+R6EjvU92L2kN3kiW2Rsr5ufC3NbOpl5LW9pHsonNFfrkL7nKC0CBTXEntPKqaSDVuwTIjSAklZaEGVI0XTW1NJjkJUWlJkjolIMP4hBoafbN/khquU25so6qM2bFOFcqZWoqwVfB7VSU6ESRah4EXbPa0t3VYSoQJ/jG+lSpiixwK+uVCBBqFKI6ry+XqR8tEzasqlKaNQraKyVUTLuwZb550bOhBTptY4zfyPtWqf0FwBQBEWRWaHY7JMAZDv9GUZdtZ1WT0hVGHLmjEx7qFhb4jNtnpJ7W6WUqrAOQ4aLdCzw3ea3uOB4phBEeOfLzsZ//6uv41duvg//8bsvHu56IxoXw3TQa10tMrb3gZRWlCotUva7VsF8Vmgb+Pa3tEd1H+16Ce16CZWqdge6DsEiVlcqyWN+vpE8ulGt8VpaDGMpeQI/dcP5+Mx9h/GH//ToUNdixTWNjNL6nhzQw/ruFiAnXd8LrO9OSye7gnGz7sxrlbTSQlCD9AjtZgOqQpBGaUmfdIFxWTsIlTBtnAQQlwFZ8gAoSJP+U2UJlGOIUMDzgHJVB5FaxcxdldpYNMcuGLt7pRWhJiPAj3CKUVc1r4V507uwZhZurIgIFdIBsLmiU35zIsLS/BoA7SC0xcOhXcMLKvkqGJslTBbn6zh8cgnBfASPOtftUiBIEKSibJP1ZMrJKFgJwPycmbtyFRTM/Fayzc5hGfUVp9sgVea4YK6KaHXNafMkHWu8TM5B7A7SFB3n3YR5lbWBVY577mNlNrXMV3z861ddgP/0yftw6JQ5vOKCUzd0HVZcTAfDzmv581pNtVa6zw+1m51KxVVayrZx6iikzSotEdrg2duxFNW9pOVSudpCuTqYGlpd0XNLC/N1LMyP1lEYNisQYyxCZpityKmLAX7qlefhX/7pnXjk+dUNXYMV1zZhoJotvaFnzVawe8lZc6ugZssz5/p+psgYAJTwoDwPSnho1k3jXFNkDKQOwnajCWWdgT4QVauQYQOxjVWVMpQPszSJ3qZKXlZpQb8Mr6KgvDbK5RjVcqq0AGCx3MJCqYWo7mHBb6JUbWGH38CCWaNr3jzPiRaqxuBRMz0Nw1UPigh7F5dRNgtbuS5ApQhtMkrLKIAmkBQge84/LQ9Kq7EC3Ka8UgndZFeKZFkU5SqqgmVGChWXrdVKaraMm1BlVZZVVBQbJRU7iipTs+U4DgvW6Op0Gqb7Mysd9yo6ZhU1U5x/2gLecOV+/OgHb8cn3/NS1MrrC0UcuKaJDRoyMtvyhox8lwyzzlbHtrz1nQRU3pzhC/iL89qIUTZFuT6lvQpND0IZA7HZFtequss7AOVrpaV8IDYGiris57VkuQlVMh98ZanTgxWJUlWrINeIYdOD1TBCVYUoLaxih6/V3JJfx7wJTotCn1sTbVRNWjBaE2ZbiFPmVzNvgweVBCYFYzHXLxwAUFIxIpM2tEFOOGkym1qUihApD1KKNHBJgVgRlAJiX0K1bcQm7SiUlJgyKCaQCHSq0AYnaRyFYWcwq/hViEhBOalC6kgbSv1I7PBIU4R2mQqV/CdTdNyRIhygNyH3JWReeeFpeOi5VfziTffg1990+brO5VQhkzB0irDHulpFPQhlzoyhCrq8p2aMFJseVD3Sg9ZmXhowJbiyUsGKSQsuztexOOK04LAos9xKdmOBQSUq7lGoClKzlmB+c1atZZg8b7vuTHzp0WP4q3ueWdd5rLimgUHbOvWyvheeXGB979b1vY/13Z+r6ZWMkwLjtNBYzc8BAJrtJmTFKi1TYCwbzirF0AXGJUD6RmmVFIA2pK/06sChgChLKC9EqRwjqLiFxaZdU6jNF7sWljNKCwAWvQaWhN42ZwwZ8arAgtfC4nwdFSNTSlC2by+EsyiWNIojJCTf/pVJJUaUFgdb1tYCxCDExoQB6D6GYUO/aKvCIiW0oUMJSCUgY/PexqTVVkyJurJtDTPmjBgQHiBiQJh0n7D7yFVa6XNivlAS0eqa+dkWKLtpQUd5uSlCmybMKySlBu9LyApqpglKHn70pWfh33/sXlx/7ikDr57MiosZGmt772nGaBSYMdo5M0Y8mO1d9ejgLhse/KCddmLvw6qjssZFMF+glgB4tZwxIyyPbQyjhouOmVFx4emLuGTfIn7n7x8e+BxWXFNKL3t74b5uhoz1dn0XjpLyBPx5rbRaK3W93pa5j13VQ83PoRk2AZ+gymmRcah00bEsAbJi1Ic5R69sbJSWMWKoikQJAEoSXsU0zi2HSSunuXIbtTCCt1DHYkkHxB1+I6O0AGBJNDAnmlhdqUAJbb6oUZgorYp5zSUSjtIy44JK3eBQiCn5sQNbQBwrAWVUVWy+J0bSQyh9yJjQNrb5WApISVDlGLLlaaUFQEQVe3OQWTOMVAAoQERaYQGA51WN3d01Zxj15Kls13c9wLTJrqu+bCwtbKjbWVisXMPGkOtu8dzW7PJ9L9qPf3fT3Xjny87BUq2/6mLFxQxNV6UVdM6d2CLj5Hc7r9VRZJyzvUc2YBarrUFbJwGuxX28TWfX6sXtl6Lm+pTVMPNbpaDad34rWtta83nM7HHKfAWX7d+Bj97xxEDHs+La6gzrJOzTRNfObWWOI2ceCyheroQE/IV5tFbryarGrvVdzet1t1rt1PoeVQNILy0yBvRclhRNwBQdk/IRlxSkWa7EOgk9oSCq+oO5Ypcq8SPMl9qQDYH5UhO7FnQz3AXPFB17zYzSAgCxJlETbSzO11EjfZ0qKfhGXVWMxT2vtgCtuOykUpE4iJ31uHQ7XkJ5PkbDqCa79lZbeogrQCw9REZmtpsBwsiDEoCMRFJsjNi4CWPqaNdEMTLLlOjfVTq3FVlVlioqEaXPZFs4uQ11E0UlgdhcPHZUllI6TdjPAu821B1mbouV1szwknN241N3P4u3v/Tsvsdy4NqqbMCQ0bc7BpDtgiEo2x0DMEHK6aJhn730HMD2IBRQnkiNGB6lNV0+0G5oo4W1vkufEKkG4Get78kCkSX9AaxKLaiS7UGore+i2kTFdH+vmsA1V26h5rcR+x6W5lcxn/Qg1EFq0Wtgnmytlk4P1kQbe02AqxrLeoUEyrb7hXntRYErVNLZ7tjczTYJgcjMzZ1cMy2tZNq/sC19hI0yQiUQmvRhywSuWArEngJiY8yICCTLIBO4hAlCEIF+ViZNGANeuQrRNgaN2EkRSlPD5fYgdAJKUKsCNqglKULXiGFfodurUHVa4PMpwnyneLt5oxZ4ZiZ4wemL+K3PPoh2JFH2e39h51ThDBPs3Jj93TcuwfWkCGWps5O53tHb+t6PuOHBD/r3JtwME4ZlbVW/1up88bi8fLeMsLPzu0WJAdOEYWeasFwNeqYJAXCakNkyBCUPO2slPLPcP4XPimuKGEhR6QP1c1H397whI19YDEq7YziGjKSI2BNQHqHZaCZrb7mGDNf6Dg9QJYI0CiuSDSibIiybtbV8nS4EoC3v1vpuTBk+ABW04JckKsaIEZh5rjk/xBxClPwWFvwW5kyKcMlL7e7zooXVlQrmvDYW5huoUZQorapJC/oQKJn3zKfOuTK7mrGLhEoNGOb7X6g8tJWHEB5ETa/J1VJ+sipyvW7cl3EJrUhvC6UHFQu0BRDb9lWxAFlTRuwUHSMwJg3zl7IKK0pVlv3d/CUT44WbIgRMSjCWQKTThJTv5q5kNkVonptHT5o3IKfCkOuSYVGyU01xmpDpwkJQwol6iEO7ex/HiotZF/4AxaphzvquSkHG+l5EkSGjH+sxZIzbiGGpr3VRlga/llVQqt3dqDGI2hKlaqEpo1yworRLZaGGqFuxOMNMiDCWKHn9w9JQiouIvh/ALwF4AYCrlVK3O/t+HsA/h86ev0cp9dfD3GumGMSQ0W+dLSBrfSdHrVkl5igyFNnhrSFDCG3GmK9BkUgMGdKZuwK0IUN5gPRIFxND9yKUHumHUVfKrK8FYcwZvu6lZ40YqiwhPGXaOsVQpRjlUoTAzG1ZC3wQhlic0/VENdFKehBWhd6/IJpQa8DexeOomZZOZVKoWHVlvreVyYNntgnnu5xE5zf+0KiAEEDbHNsypoujKwtoKR+1+RCtWP/TasRlNOMS2o0yWrGHKNbHt83+WAqocoy46UGGRl21K9r6LvTcFgAIGUCESOe2YBSXfU5UmNJFyFKB1taS4yhnpqDIae8Uo7O9U6yQmdsCits7AeAVjplRoJTC4ZMt7NvRfzHUYRXXPQC+D8Dn3I1EdBGANwO4GMCNAP43UUEOhpkIG53bAoBWn2/p7Ub22//gaqsYVVqfCnOx81qbwcqqVqK1gnmttumSkS84lj0KjgdVW0WUq71V30ZbPHHRMTNODp9sYb7iD1THNZTiUkrdDyB1pqW8FsCHlVItAI8S0UMArgZw6zD32/YMus7WgPsyOMqr5zpbwl3t2DbPJXhLc3qZKM9xEDrFyAhqWm2Z20kvVWOJoc0zc1tw57UA5SuA2voZAPkKXjkGEVDx9Yd94EfJ3FbVdM6oqRbmTL/Cea+VdHafM05CQRJ7FvSHrVtgbB2Edl7LI5FRWnlsm6fI0WCRIrRMwD26or8I+PMSq7EOJg3j92+qEuqyDL8WoWHSgs2whFbsQUUeWkRA6COOBCCFbqgbU1JsLCICKMha360KU4Bn0rK2zROF0P+qm/Vse6fIuDRjPXcVnVxL20BJ6ay55agna4c321QcZ+e27HFF2HZYGbu8dH7kYmMmy91PLeO6c/pMbhnGZc7YB+CLzu9Pmm0dENE7ALwDAA4ePDim4WwzNlizBZiehEU1W3abSINZYm0nASU8nSJMrO8iSQcqTxs6mk7NlvS12lIC6Urzpgch4iasyFIlBeUrKEoNGcKX8DwJigklkx6seBGqphC36rcR1X0sLq4lhoyqaCdLkwQUYnWlggWvjYr5tHc7YljLu5/UbHW+nxIyMWXE5sM5VBItE8RaykdLelhZraGlfMwttLEqgyRgrcX6jTi5OoewLBBGZTRMerAV+wgbVZ0ibJmO8ZEA2ubNkjHIGDUoAiAA0W7DZEB1arBUBYVOitA8B+UAFCptc0+Ol0mQCmoVff38QpH5uqxMRwyndiu3rElHzZY914UD0chwPy93nHrGhEczWm5/7Bj+xcvPGejYvqlCIvoMEd1T8Hjt0CMFoJR6n1LqKqXUVXv27BnFJZkR09eQUWB/t7iLQ2424zZk2PTg3EJxOrPdKE5VqnbxdlLlzi4ZojPtZ1OERRZ4AECzdzqXTRnTi/t5Ob+0c9LDGRnPrzTx2NE1vOLCwVZE7qu4lFI3bGAcTwE44Py+32xjihim2HjA7u8QhGDnIshzU4HFxcZKdKqwZkMvDFlUbAwPyeKQbopQmTRhoq48/QBZ67tJE3pK571MywfPlyiXIpBSqHgmVVgKUfV1gJgTbYRCIqAQVWE6wos2KkZ2tFc9lClChSL4dl0so7h8dPYgdLHJwFjJxIjRNtW5LSi0TOHw8ysLaKoS5hdaOC51AKvLClaN0jqxuqDPKftohFqaNqMSZKhNGkmKsG2mftuBTg9KAoWUKzZuQ0RpqlAoHbREpNJO8JGe26IQINfyDmuJl6jM1xCfWNGvPlKpEcNNFSpnW9KPUKU9CYfsR9hzHyuzmeUz9z+HN1y5H0FpMCvEuOzwnwDwZiKqENFZAM4DcNuY7sWMEVtsvF5UKYBs5RQPlTMFx+NkYYyFxidX9Xsyv1C81leroVVSKWd9T8wY5axJg6QOdoMYMrxyH0NGF7VVGWLNrebxkxs+l2H6UW9H+PsHnsOPXH/WwOcMa4d/HYD/CWAPgJuJ6E6l1KuVUvcS0Z8BuA9ABOAnlCqo4mS6MrDpwlJUbAyk20ika20l+4uLje1CVIoEFBm1VVRs7GlTRlJsbCzvACCN1V3b4802O8eF9Fn5CuQDEArC9iX0Y/ieBAShbFJnFREjMD8HXgjhAVUvVVkVCrNzWyJGGaajPFIjRhGu7d3Oa4VKIjLbW0YJPL8yjzXlY2GhBYp9nJQ6WDRMz6qVOMDxNW3UCAOBUJbRiEuohyXIdgWtdgmoRGi1PUSheT/telsyBkwjYRERhAqAGMm8lm35RAR4phuGFyGdx7KGjEh2FhvHusg4Wq0nFnhy+hMijguLjfv2I8wf18eQ0RNWWzPLZ+9/Di899xQc3D34l6thXYU3Abipy773AnjvMNdnRkewc2Hd5wxSbDwJ2vXedtlxzG2dWNMqa6GLygKAplFa5VqIepSOUdo5rUqnoiJZ0fNaoftlw6in9RQbR91Tc5WFGhBt7HsjW+CZcdKKYnz6nmfw4Xdct67zuOXTpBm22DizPVdsLERmnirZ5xYbU/YcW2xs92snYUGxsSkqJi/txSqtCoOe15LSKDOrrgSgPHO0++zp+S3PM6sP+zHKXgx4hIpvXYUhSkZdlUWEStBChSJUhC0sjlEyD88oJaLiZrkWt5WTtbxbB2EbcaK0nl3Rne5L8yGWjcpalRWsmZ+X4yqajQBxhdD2BNpxBaum/+DJtXm0Iqu09BsRRR5kKECqopdStvNapgBZxNpFCKSuQY+qWnFF6TYKFYJyAEgFWu0sNqZYorJQA4Ux4uVV3QLKBjjb8km/aCR/RacjfGZuy5xTWGw8oFriYmMmz2fuP4wXn7kLF+xd3xdrDlxbkG5pwn4LRKbb3JotodVWxoRhz0HHApE25efP10ytVq7DhjlMeQRUa2g1G1AVE9SctKAS+lrSc9KDzs9JAPMBEhLkKQgTuHxPouxFUB5QMvVXZSFRMS0jKhSjIiKURZTs9ynG2koZi/Nr8ApWdrSBSUIlpoui/S0TzCIoNBRheXUODeVhcaGZCVarKsDJuIpmo4rlOEC5FmIlDLBqegzWozLiVgX1qIymByDy0TaBK257QCuAgrGqx8IELmhDhmt9L6jZ8pJUoPlL9qjZokgiXl5N03mFhgynZssxZDSPLqfHOvu6wTVbzHpohjFuvusZ/Mk7rl33udyrkOlK1+7vM8KyMWEsLhQbSpoNbZQo1zq7ZcQtHcCoS3oQwEDWd0AbMgpThH06ZAxjyGCYcfM39z2L687ejQv3rr+TDyuuSTHOBSItghDsXswch3xHeLcfYXJtnS6026RPmRQhgKT3oPK1EcPuU47isvZ3ZV+ChzRF6KepQuEpwJMomS4ZJSHhCwkpJMpGZZWNwmrXS/C9NCUoTFqwsVLCnGjBo+w3+DQFaN872ZE+1CrMKC4oLK/Ooa0E/HkdWGx6cE0GWJY6WB1e2YWVuILKXIjVUO8/GVWxEgaQrQrW4hKoEqHZKqUpwtAD4gpkKHQRsEkRAoAXVVPru1NY7FMVCLXRwktShHpeiyIFWl3Vf+rQqLBQJilCRHFiyLAqrNiQ0VlE3Dy6vLFi427bNnIMs21Za0X41N3P4i9+/PoNnc+Ki+nAn6v17Uk4SYJqsQIahQW+n8oC0uVJKnOdSkv2UFqIi5UW9TBjAChUWwDGVmjMhgxm3Nx89zO44QWn4pw98xs6nxXXZjNosfEA2zvImTOIHBXmWuDt5WyHeCC1wJvC4mRuKzku3Q9oW3xybG5eSx+nv5RbS7zephL1pZKpNQnhESAUPGHs8F6MsoghSaBkPuQ9o7A8SJBRVR6liqsIpbT3wfwGwAqFdL4L0M3RQ5BRWmTms0qom55Wdl7rmdXdWIkDBHNttGIfy6EOLCcj/XxibRFrkVZadduXsO0jbPtAXIGKBJQXgdoiaekkzLkiasO0XdSKK9RvuW/ntUKA2qZtUynQLZ0ilfYvdAwZQTUAIq28AGPI6NX93VFfRcXGhdZ3933mlY2ZdXCi3sZn7z+MT/3Ut2z4Gqy4tjHjtMBTpYqwz+q6nbRTdbHF6Ke0Vut6fzDX2d4pbAUIW/p1rUdppXRes986W6pbofGC/vttdGVjLjZmxs3H7nwKr79yP/bv3PgcLCuurUq3uS2guNg4v95W7rhuFnjVcT6hudbUrZ+8dF7LVVr62alPTea9nPku4TwcFQahv+mTmesiT0FA2+GFsA5BCWp78IRMvlmVSCJs+CBSiWtQmOf6ahmnLZxAbGRc2zTPJaRW7TgRlgqxkZx2BeMjq/OIlMDiQhMnY1NUbLq6A8Czq6foa1SB5Vj/Y1sOqzgZVRG2AqxGZXhBiHosMkoLAMLmXKK00DaNc0OCaBNIBBChsb2302Jj37W+JwXIWeu7tbYnDsJQIpivgaLU+p7UdsVF3d9lxkmolNJBq2Bl40K4+zuzAZ450cCXHjmGv/vXLx/qOhy4NpNxL1vipAgTC7ybIrRPrgXeBhrXpCFyQcpLf066YPiALYNyzRdJKtDTn3s6eJnJfaHM/UxvQmgrvCChU4Wmr6BHCkQKfrkFz5gzBOkA585vSZCxZwjEIIQmT+klH4w+lOkzSI5F3gasNjysrFYRKYHSfIxlGaBuHH8NVcaarKBer6EuKwjm2jgZB5n0oO2SsUYe0PKw1i4nfQnbDZ27j6WxmscEaps6rTbBiwJAZdODFGkHoWgBXj1nfW8DVILpjpF2xhBhbrmS5dW0R2FRkHI7ZGRqsVRH9/eBa7a4OwYzIB++/Qm842VnY+dc97XoBoEDF5Pgz9fQXK0DAza6nGaWV2uIILCw0EKkBJq5xSxtwKrOtSBl55cHmxr0ghBo5d6vSO9Tfgi0O/+JkQg6jBhA2oewq/W9W6YRw1vfm8c4RciMl3ufXsaTx+p4+0sH70nYDQ5cW5G8MsuvmYXeKcKMBT6fKnQNGxlzhvlZUCZFCBgjRqLMnG3mtnkLvD7OqC/7QPq7AiCM4hJCQZACCQlP2I4Xepug4nRT7NjZteIiRMpPUoj2Y7+kYkRmvEJZQwbhhOmEESyEWJVltFQJLVMZXVcV1OtVrMkAsgq04hpWTCBaiQIcr+8AADR8wlqslVbdNM9thCW0G3ouLJYSaHtAi0BmzS1qE4gCUCtM7O42FVgirbQoasILi63vYqW79Z2s9V3G6TpbVmVFUVoPUGCHbxw5kemOoQ/rbX1nQwazHmKp8P+++E38+9dcNHAH+F6wOYMBoC3w252V1Vqyhla3noP1ulY91bnO/e2W3ucHxetvZZRWHgqMMaXLueiitIwZo2vXd2vG2MLlCwzz2a8fxqmLAW68ZO9IrseKazMYtQU+X2hstvVs72S7wOdWOHbntuwKx64Bw57jznfp+6HjOEVIVjvOznHZc1T6Vlg1RVphkXnWt1NGqKXf3FtrJczNNxArAWVNGMqHr2JESmjFZK5tFVkEL1Faq6s6ANQW2pBKYFnquS0AaKGMhiyhXq+hIUsI5tqox+naWitRBWGripUoQMsHEAZYbVewZowYjYZZe0tJxLEAYh9IVJZILe/tEADBC/WcFWBs76UqRL3pKKnU+i58BdWsQ4TZFY3tM0US0Vo9a33Pz2fFztyUM2+VGDLcdbbW09ZJb+i+r+AYZvY42Qhx0x1P4c/edZ1eoWIEcOBidF/CTbgPqTIUshZ6ikpAKbduVdOHP7e+Iti1ehWV+c5zbMCya2fJgqa7a/UamrKE6lyr8MM3dJRWK8pOKkvjQkQp1KnBPIn9vw0U3NsWGRcxUNf3cOOrBbH1ndkM/vT2x/G6F+3D+aetvzynGxy4thK9LPAd23OFxYI667Yyc1yOkzDZhmQOqrVW1yrKcz5crXpyVJP7nHyRzs1hJecm11fpcY7S0k8Kwu9MkeUDjFSUbAuN7NNKqwRRA9p1D0dW9dyeb5yEAhJz81ra2FZNEiI5v6181Os1tKUPVVNoxTU0jAV+JQpwYm0JANAyK5SstivJciX1dhn1un6/29DzWXHLhzKWd9EWAAUQbUqUVmJtb+uCYlGqGjehmdey+0OFqi0yXkmt78JRVRu2vktnXqvA+s7rbDGj5MHDK7jryWX81puvGOl1OXBNiHUtFFnUlzC5UEH9llOT1bGsCSgxWBT1JQQV1WxRxpRBlaoOXOaKKh/gRIAobuhASCrZnwSz3LB7YdOCEgpRErD0h3BblpJA6dfSjhrKbFQEnJQ1c76+YawEWsZBeHLNBPqqQsMGJKOolhuLqEcVeEGINWNxXw3LaNg6rfocwsiH8kNI6xxsCR2wAIgwTQ/atKBnps28EPBVFWgDXkMHba8NCKczBrVVxvpOUdpvsBpUgChGfGIlCURk19uScV/re/P4ykit72zGYIqQUuGDtz6GX/jOF2Ah6L2G3nrhwLVNCHauv8My4PQl9Nfn9Gm3GkBpNPnqzaZhDBgV0wWjJbP/qKwJwwsKTBZOalB1S7AOkh5s9+hBCPTvQ7hW77HSGMNMnr/9xnPYUSvje6/YN/Jrc+AaJxvsAN9XjWUKi+0mJy2YUWTOOYAxZ+gfFVHySNOClDFl6OOypozkXMeUYfflLfAdn649XpotDJYgRKYFR6Ssnd1Dy1Q/E5lAI1Mjhk8x2uZ/Z+F0iJdmkPa41foc2tJHZS5EK9LHN+MSmiZ4Ldf1F4CWDzTaRoVZu3tzAa22B5RChO0SZGjelJZRwq1OI4YIke1BWKqCQldpOcXEka7XEivZImN9rlZclXmzmnGkU4KFXd97dIJvHDlh3+zxWt85RTjTrDYjfPQrT+JDP3btyAwZLhy4tgHBrg2qrU0yZQAAVBnI3a3ImNGLsFFBqVpsY+9Hs6FVUGWunTSkdWk3s1b3DhNG6JgwCiBUzZeDdrIlT79u74Our8Vd35mtzkfueALfddnpuOiMjX029YMD11alV19Ch2D3UsdxmWenL2HynMx79TFldDNnWGMGFRzvBZlz3XmtvKMQAOC10rkyRVrq1FaTvoNhrP8X9YMQYcuDkun/skpQYrQoUZxp62RpmMUeYyVQroUIwwra5hpWZa3UF9CWvp7PMiqrGfloGKXVbOqC5RYkIjOfpUIPMC2cqK0DighDCKO+rBHDazmW98SIkSotz1jgAy/Q5zTroFBmiowBIKgFmSJjALrQuMj67nZ9N/sax5btm6yfXZXVxQK/Iet77jhm9njiWB23PXoMf/szLx/bPThwMeuCytW+81tx3BhLabtfixDVyyjXequ0ViOrXCpzYZJ6zGOVVtF8Vmp1NyoqLGjfBGNnp+L5LKC/0ioFVSBUvL4Wsy348JefwLtfce7Q/Qh7wYFrXGxCQ91g12KaPy5qqOuuqWU7opOWQLovYQOq5KfnwyivxEGI9FzXQu8+kH1WlE8IwqQJc5scYZBOr+gLhO0AJW8NANA281qejAGUgQqwumYUEEVJWyjdIiq9fqkWagu9IjRDH6FVcMpHO/YRtoJEZbVigYbxvDeNwms2F9AO9XxW27gK47anlRYAhAQKc/NZbcooLQDwkVVagFZbtnlu4AVAqECr9WS+SoRO89woRmWhhvjESuIcpMhRVFZpWVdhrut78/iKXl/LHucWIEtHGbH1nRkBdz+1jOdWmnjrdWeO9T4cuCbNRmq3gGxqsNtx5PzspAcTU4YQnWlBcgKWDWCU258LXEla0HOUTn5oyqiNJC+oH0oSYmvEkAooSXixhzD2zLBNYIpLkMYG71f0syQfvu0e74RLCUIj8hArQiythd7UgEkfzYZO7bVLCmhX0Ix9hMao0Yo8qKiKMBIIIYGWn5gwVFuAQgIhAIUCItJKjEx60O2I4TlGDD9nxPBCBWECm+34LkKZ2OF1Rwz9uoKqtr7rbfkehHG2VsvukzkjhuzsjNFhyGDrOzMkUin86Zcfx89/xwtQ9sfbTZB7FU4pG7W/A4MvFpmHKiZN2IM4Wu/ikptH1AoQtSsQQQhRZHUHoIwrsLsJwwRn6t5zEOifHgSMGaNHepD7EDLTxK0PH0W15OE7Lx1NP8JesOIaNRvoS1iIU0zcYcroVnScV2D5TvBA0s2itWJqt4r2O0or3eb8DEd1Ze7X5bU4aUH7TGEZymvpdJ5ZNkSaG0YxodmsQZTS4CA9SnoL+kZF+SRB1Pm/cN5WH7crCM05cUmiZYqIrcoKYw/tpg4SsRRQfoio5UNFJi1o19GKaqBIOwd1Rwy9O+mIETodMUydlnCNGNba3gIqJb1MiU0JirZK+g2KUCKoBkAks50x8uk+6aQA3c4YbpGxfkMG6oyhf2TrO7N+wljiI195Er/xpsvHYn/Pw4prCtmo/X0YqNy9px4AKC/oorY2MEFb0eYLGQ43uRu1K4jbulFuT5XlFhUXdXaHY8Lo0d0dWIfSAvp3fF9jMwYzHXzm/sM4/7QFXHfO7k25HyuuSeEos3WZMuw5hd9q8vNZqSxSGWVFiUJT+a8ulFNaSOe8tJtwsKFm57LMh7ik5Mu5Iv2zkgLSfJOPzJpRBID8GGj5ifJSihCZ9bp881y0XlfcClIhUZaIzXsRGpt75MyfhbEAogCRFIiUBFolyMi82EhoAwYCILSKKwTZn02/QSBVXK4RQ7hGDGttb+ugRRHStbWceS0RmrW1whjx8kqmByFFsW7npF+EebEyVTtGUTWPnoRyzBn6uLhvSye2vjMbZbUV4RN3Po0/fed1m3ZPDlxTxrBqy5+fG3nRsfJ6F852g2QAFNV2WSoRyIshW2VIESfLnsS2J2FB7BZBmOyIu1jgAaOyTJBSfgiEnS2vOuezutjdfaOymk1QD6d+P6Vl2Wg7J1ZazCT42J1P4dUX78UFe0fX/b0fHLi2Et1WOu6Yu6LOn4vaPHWhVTDZr3rkpSmoomVMGd2Ok5FWY4lYc63urvqyX/xLESj2IWNKFGccp6/BnhIrAjwgJpGskCzM8a7isvNjCMuQJmDFUujzkQaxMBKQxsYemsGotg9l722LiqMqKBIAtVOVFVK6cnHbFBb7VXixVlnUTguKkzZPkcoorSIHob62RBBUEFnbu6u0AK22onzz3NT63jyiC4xVHGeLjNGlpZNDN+s7Ky2mH08db+CfHjyCW/7Vt27qfTlwjYoNmjK67us3wZnvBF9oznDquJz+hKn5oss5eZt7D5SfU1tukKIyIBvJhSgG4BmTQExATIAAZO69U45hwzcdPWKhQLHtVagyz/acxJSRrPAhUuNHLABZQRwLqEgbMGwXDhWSMV1ABykKQDFBhLY+C8n43b6DPlUBp++gcJcmSXoQAoEfACFAq2l60A1YABAElUx6MBOwgGzNVq4HYVKrpd+IzkUhizpjKMm1WsxQKKXwf7/4GN7zynOxZ6GyqffmwDVFTMSUUan2TS3KISzwpPqkC0dBVEkm6pQfJsuedIwFgQnctgtGceT2SlWtKOM+JoygChErYKV3arCyUNNKa5DXwjBbhH988AjakcQPj7nYuAgOXJvNoB3jC88d8KPN7U/o2twNytmWqKoewwqbDcDvoRpzRcXJzwAgKWlqqwS0ygJMX8EIiH0oVKGomQxCeQrKKiWT4iNShUoreU2O4krOjQKttAAok+NTsQ9l1VWcpgBJWjegDloipGS+yqooioAS9JIkotlMi42TDu5IegwGZiFIsaJ7D+r9WSMGgO7pQVdpAbrju6u0gGwPQtcqn7e+A+vujMGFxkwvjtfb+PCXH8f/ffs18L3NN6dz4JoSJqG2+tGRJtzINbwIBE8rr1imvZI2SlzJ5DgTe7ssDryJygL6FhV7fhWIeqssQCstxAqq2VtJWqW1UdiMwUwCpRQ+8PlH8YPXHMKl+5cmMgYOXFuBXirM9iDMW+E71BdlTB1d6fPlyJ3XokoV7UYzU4pF5ts8maaEcbuR/F9EMlVfmWfrzhak57nc+0QKyrfzOR4gq0CkEFt1lbSfUl2Vlr2gkiZQueLCFB6rOFVZdj5LSB14KQrNtrTXIMVZmzsA+K0myM5hRamC8owyo7ZCtRToOaqVNXNt1dnpvVoBeTJVWuZ9KFzFOOo0YmTmtXr1ICywvg+ltHh+iwFwy/2HsdaO8FOvPG9iY+DAxWwZlI0AwomATuDCAIFrEIgcpdijqzvgqKy4t8oqBVUIoaBM0XC3KwbzNSCON1xcDLDSYibHY0fXcNMdT+Evfvz6sfcj7AUHrlGwgU7whXRp85SkCYvmuArW6HI7wruFx/78HJq21dMoUM7D3lG6z20QBVBSf+hTrJLx2PkjvZyWnZOxDkgPEGZbIt16jEEfmCotq64kkhShiIxbUNq5rbZRWUIXBBsx46oskgDVm5l9wpnPsnNbFa8KtBRobS2dx3JWMHaVFuJYt3FylJY+Pjev1WUV44yDMJb9m+e618i8bzL3KystpjdrrQj/47MP4j+/9mKcvWd+omPhwLWZbKRbRhFdLfQbv2QHQU2nCQHAFTv2s9ALACuMbGyRSNKCiT9EpqJJxpQcrIzZg1QasMizKkslL0YVmEsS8iJMpunAZDzKdHO344pDEzTTtKVw0oIeVdNtTWtzN/sildZxRUovSQITpJp1iFAlAStdokSCYtMRwy4EGcviOq1RGDG4ByEzBqRS+J1/eBg3vOA0fM/l+yY9HA5czMaR4YA2dgqQtH3aTCjQc3EA0h6DPdKCpTQt2KsDBpAu/qhajXS5kS6Moss7pweZSfLRO56EVAr/4bsvmvRQAHDg2vIEuxaLU4QdXd0LumlsAFJp9o3cL+6uqoKpcsp9sSeZlEulCofauhAZAUTcTDtcWKXgUaqa7Bd8osS8UWTISFCUUYKkqpkxQrZT9ZU8p2lBG5x8MnVZbcBrNZN9IrHDm7RfDFSEVlnCLv6I1HRBISBiuwCk1PNZnkRsnIOZhSAT1dQlPTgCI0byNnFnDGYIvvjIUXzhoaP4y598KUoTsL4XwYGLGT/KBC8KQALJnNcwZGzsyX3Q19Lu4pkeg/0KiS2loAq0tcpSseqZmQ3m0w7vw2RwWWkxk+TRI2v4wy88hj/+0Ws2vTtGLzhwDcMo1t7qV5BcZL4YN5UawrpTdKzQ0YOQpILnVSHDRvoSZKpm7Pd1QXbqqW1UWxmEIPOFPvnmn7SlGlA0qnaqABXMjSi5OcWp+kpUVmxaNUGrKoqb2sToFBnr49P2TZWSUVnGgAHoOa7EgGENGVKBQq20KNTOQYokSDpKC8isUlzY6X1ERoxe1nee12L6sdwI8VufeQDvfd0luGTfZOq1usGBawuzFYuOLRQ2Aa/3Gl2FqCJFNECXhtTXsWE8r6oDYgwgbibpzF6UbDFxqwEV9va/BHNVUKRMwBrg4j1gpcVMkiiW+B+ffRBvuOoAXnPZGZMeTgccuLYCReqsoMh43SuLOtf15+fQOlnv3rrJftFWKimZyhQRJ6ooVWG227tVNa5SSoYaOa5Cr3NeTAnnmu65XVyE+Wd3fiyjvmAUV6xTghQDiAHR1grLjg0wxcRx2hQXAColU5e1Unf2ua5BlVjjKTYqK7NisZnrknFWaQFaZUlHSdk3xOn0nlFawHAOwkHmtVhpMQ5/fNvjOHWxgp/59vMnPZRCOHBtBoPa4G2NU7/0YLcANmgfRLfuSuU2KaQfiAU1Um6dlpAKwg8QmeJcSg/Lxh2nrljZEjJr4lCAiotO6j9+AtK0oBu4YuMQtKnCNuCFZg4rToNTaodXTsBy0oJ17QIUjsU9MZ1EUlviTUExQluflTVQUCQ704JSom+n9800YjCMwz89dAT3PrWMT77nW5IlhLYaHLiYDUNhE6o0fL/CUWHNFkkN2YDpQMt60oJA1oABDF9Gx+lBZtI8s9zAH33xm/jQj12Dpeqgy51vPhy4NsowXd4HIDO/tVF7u1S9exNmlJdJQ4GSNKCbKrSqiIS1hhMkAdRqwguqkK2GNmLkbkFAZlmoTIrQHpB2bUq39RqrTLf5dp5NpmaLTPcOp8jYbktTfHpbpRRAEKBO1iFkQVrQpgpDBZJKz2UJBVgDhknxZQqLZVZJAUgVl2t3V0oHrFhuvhGjy/HMbBJLXWT80zech4vP2FpmjDwcuGYAf35u7PcQlSrieMzrasHMVzmBCwraaOGkCgelXLFBT6usQQjm9DnRWj0NSEPAKovZKnzyrqexe66CH7n+zEkPpS8cuEbMUK2cxkjrZJeuDcnquG6bJGvIMHuEW9SbHp7MEbUaUKXANKRtZC4rlZO6E6mqchWXynozMiS1VtDjQhsQYWcrqsx4Het7aodPVValbNo6mZoskqnRQjhzWOl19L5aJQAimdjcXQMGAFMP5qgqAIgiIM4pJGOHbx5fSdSTimPnTRuzEaPgWGa2ee5kE5+6+xnc/J5vWb8JbAJw4GI6adZRqtXQbg9eKExhE8qrQpSqg7eC6oJXygUra/4Y0l9Qrhh3oYKZxxrsglZlwVjdRwErLWYr8aHbHsePvexsHNhVm/RQBoID16QYxxxZ8g264NrJt3lK5qn0t3jTzFak+7VbLf2GT9Kx4tsefqRARiMlKqXRQBxU4XlVxEZ5CeXMcckCu7yZ4xKlaupIbDQ71ZdTBJ2urpx1OVrFJdxxxyZgeWbs9XqnuopVoqrcru7JfFaoVZaUMnENum5Be26msDhp2aSySgtA8+hyVmkB5v0ejYNwIFhtMYYHD6/g0SNr+ODbr570UAaGA9cWZOjCY/thR9p0QSr9gCMFqIxbwnxgOz0ElbG/l4Oq7qAh4fQMNIEuAoTvGDqgP1u9hk4bkhdAthtQMu07CCdVmAQrE7iomVN3OTHU0YXeHON28nANGABQLge6ZqutgEbDOS5ruqBIJR0zrKlCRApBzaQGT5pFIeNYGzAAvfCjDS6u3b2oPstsax7V3d1VHDtpvy7pwQ0aMTr2d+7svo+ZST7+tafxU688D0FpRMsdbQIcuLY5/vwcWitrwHqbYzbqwNwcSrXqulKGgE4bQgKiXIUSlMYgQodzkKLm8D5yh3KlmgQ41ahD2YC0zusE8zVASk4NMtuaJ4/X8diRNbz+yv2THsq64MC1rVHFPyeKy9nsbLOd1FV9DajWgFilqUK71mPs+DTMOloClLG+I2rq2mUBrcKiZsf6Wt20QWGq0KYAHeWVqKtKFQSl7eh1V10h+RnIpQUzqcI07RfM1xCdXNEGDSAtKlbS6YLh9hZ0lFdBYXHz2Em9yaYFbSrWnuPus28eGzGYTeBvv/4cfvDqg1OltoDeVT59IaLvJ6J7iUgS0VXO9jOJqEFEd5rH7ww/VGZSlGob6EmYg8LRr8dVrlRRCvRDtRrJY6ME8zUdtEaksgBWWszWJYwl/umhI3jjiw9MeijrZljFdQ+A7wPwuwX7HlZKXT7k9bce6zFV5I4ttMrnto28sa6dx5JW6aQ/k+ucSNouUcbIgbU6MFdDuRzo+a5k3Fph2Vvo45U2QQB6ngxG1LnFxnYF5EGLjVW6LVFXge3u7qqrtIg6bf2kMkpLj9EpQHb22fmspNeglOl8lj3Obd8kZdaAAWQKi63KyhQWW1XkWN87TBj2uFG2cmK1xRTwtSdP4NxT56fGSegyVOBSSt0PYCp8/8wQNOtAbWPzXaMiKRYGoJomgMrea2INijufNYrrscpipoGvPHYc3/PCrdf5fRDGOcd1FhF9FcBJAL+olPrHooOI6B0A3gEABw8eHONwthFKdbaBUiqVMfZbuLYQojJfQ7PeTI6zDkLlzqXYBr+xSguCrVObCFRfA4JaOsfju+uMpMXLqSPfqBSBZDJMOwgHDA1GgSQBK3KdgfblKcdpmM5npTZ+1VFEDJUWG0OaJrlhjHhltWM+C0kxstuqqWA+C6mSss7BxO6+noa5yWvnVk5MMe7n5Y5TNx50Iilxx+PH8V9ef+mohrap9A1cRPQZAHsLdv2CUurjXU57BsBBpdRRIroSwMeI6GKl1Mn8gUqp9wF4HwBcddVVQ5aYTo6JdcxQ2cCTJ1pZg78w5/QidM5RaXotMQmQSJcpMSk+EQHStxb5AGGjAUQAeSYA2g9pkdrqVSZYmSECUANUEVcq1dS8sWoa2Eo3bejUaTk1W27dGWCClUrrsvTrTG3sQTUAwgjRyTVQLp3XkRYETN9Bt5QgV5917GTf+iz93hTY3fUO81TwHrERg0H28/LA+Zds+PPyG8+u4MCuGk5fGn7+ehL0DVxKqRvWe1GlVAtAy/z8FSJ6GMD5AG5f9whniHEuHFlZMLb4YWjUgVoNpWoV7Xzd1Qiwc1e2s8U4qSzUkvZNo4DTg8w0ccfjx/Gqi06b9DA2zFhShUS0B8AxpVRMRGcDOA/AI+O410zQR1UlJGk/R+6QQrSyCn9xQacM1xppGi8mqKTBn2dO7fwSp0S6hpVaqwPVGoRUaSGzTTMSkvW27FCVO+Qu4y9XzdIosdLzV0rZIaapQIWM0tL7VKY/YWrOcAwZiQqTyb5KNQDaaXoQUjndL9xUYU5xyZzKMvdrHjGFxUr2LyxOxt1pbR+qsLjHeQzjopTCHd88jp/8tvMmPZQNM6wd/nVE9CSA6wDcTER/bXa9DMBdRHQngI8AeJdS6thQI2WGIlrVaquyMJpO8f4ILPLlapAELdVspKaLMVKxa2itjs7yzjDTxGNH6yh5AhfuXZj0UDbMsK7CmwDcVLD9owA+Osy1mU6swiE4lvV0Lzq85VKlE0NSITq5Cn9pAZWFmk4bOjZ5MjJGJSb3HO7G1TWgNodypYqoXk/mwhRlxZ6+rnOeAEpVJ+DZ4KFy7kBn2EXF0onhwm1ZJZXTwT5VWYmClBKVhRooihFZpZUUJw8wnwXAbd8EqdVVsmKxva87twVHZbnX2UDfQZ7XYkbFFx85itdcdsZUu8HHuxoiMzDjnN9ycZVXZXGI+o2GDjp+bbBrlGrVNGg16/qxiVQWRq+0eF6LmTakVLj14SP43iv2TXooQ8EtnyaNVB1FyAOd41nru/mm7XmOMnBtfO4kkf7dNo31l+YQzFfROll3jO0SynyfyXRjt19x7PUEaeUFoFytmdtSclxHgbE5Nvm9SEDkhqpfn93mqKtEeWUVV34+C0ohsCnNSOqgHTv73TmnfvNZ9jhzv8aRE4PPZyWvxVFu9iWz0mI2kbueOoFTFwNcMMVpQoAD1/SSN2y4tV2uSUPaVRwdO7iJBtHyGvylOVQWqmitmAUggWR/En2E0x1XOTVb1kGx5gQla9QY8rUVWd8zQcoe5waPON1fWazpbZFCtLqm04bm9Xf0GFSyf1rQnKsy6b5sKcGG67Ny+wr39ziWYQblb7/+HN5y7aFJD2NoOHDNOEnqcHG++yrJU4JNB1p3oVVYo6Z5vKMckWG2PM8uN/Hg4VV87+XTnSYEOHANzjgWfhwEKXUasBsZk4bKPGl1lFNPSFUYkYIyneCjlVVT69VIs3jWNg+RZh9dpef2P4Rt/WVSZetp6ZgXEG5HjMzaYfrHoBYkrzexsQNAbNKBNr0oZWd6TrkGivQ8ON0v0vtlVZOS0umMUZAWtOc552Tunfy6wfQgKy1mCD59zzP4gWsOolqerk7wRXDg2gJsljFjECqLNTRXx29LHxSrotyarWSdLOUGroLaKIZhAADH1tq49ZGj+NXvm84WT3k4cG1BlJQgYeemOq3vCo5KcU0a7jlA1vhhP9QF0vmuWILMApMqJj3ntZgaNpC0b5Kp4EzW0yLn59w+ZPscFuGuyuy+vspCTa+4bLeFcRKkUqOGma9ylVncqZCgnNZKrsrK9xjMqLD0GpmVomFWMO41n5W/t71LnxZOPY/rcTzDDMpf3vU03njVAZwyX5n0UEYCBy4mQ7S6Bn9xHpXFWmLYGDeJqrJjWKtnU4RbgGSZEoaZMp5faeELDx3BZ3/m5ZMeysjgwLWJKKlG34w3823fWuOd5YXtF3Ur0aTTrl2IRH0k41KkC5Xn51CZD9A6WddzV1Z9WbUmnKJhV1XlFFa3V2s7WGjnX5r6s0NNFVnnHFeRRd6dc8qorOScAhVWYIfPqCz3OKfFk75M70a56a9sd2cmy8fufBJvufYQ9ixsD7UFcOCaHPYDyTV99OtJ6H5wmmcqCg3JdZwDC8wZOidnU4nW+i5AMRAvr8JbmEdlrobW2lqakoNzfDJOZ1xdur/nVVV8wi3eVWnqL/8CC4OUO69VUHflBiuVS+cVBbiCtKDb1d0aMjpMGD0s7RywmK3AE8fq+OrjJ/Cbb75i0kMZKRy4mK5Eq2taeZn+huvtLt+RAnS6VkxLsxlOETLTzJ9/5Qn8+MvPwWJQmvRQRgoHrq2A/XJdZCF313MSBQe4potkm3O8jRDWkEAEeE7aMH9JKWFPIkGIT2pl5C/MI6jVErNFa2Utc25lYa4zGkmZ1InBXjUvNBxhUTiv5aqrAqWUUVlFpou8gcJRSl3Tgrnju5owBikaZrXFTIivP3sSTx5v4IevP3PSQxk5HLi2AM1jywh2LU16GD2xAUgRpSpMFB+zXWC1xUwrSin8+e1P4me+/XxU/Omv28rDgWucFM1jDXyqNU0UnFs01yVRrJ4Si7xzgu0mIaRzAbvcMXXOezkXcFVYumsdbZ7yLsGM+iiwyHdTV8m2AkXmKqR8EbG7rcjY4vYizI9lkPZNPdQSKy1ms7jryWU0wxive9H+SQ9lLHB3eIbJwUqLmWaUUrjpzqfwL7/9fHijdjFvEVhxjQj7bXqjdvfmsWUEp5h0odsw1/4O9J7rAvrPdwFZp6FEWqDsqhlXfVnseOLc7xulqAC5SLnk7eyZc/vY4Yvcgpn9ncXESdDajKLiLucwzDDc+/RJtMIY33Xp6ZMeytjgwDUoQ6T9spfJBTglk2s2j5i5Ltch73bRsLgpQHsckC4MN6hhQykgzq34KKgzaOgBd76YjSzHksf94F5PkMrv72a6yB/nBrhkCMVBr2iMHLCYrc6n73kG/+Ll50BsU7UFcKpwy9E8tjzpIcwsnCJkpp2nTjTw6JE1vHYbdIDvBSuuLYj9AA12LRanDIu2wfFZZHoV5i7uqpW8+gJSBZa/jyhQCH1ERc9jClWVc1KHiaNAFTnXyZgu3OsWWNkzhcRmX/NI7gvDOFRWj3MZZhT8/Teew5tefABBafs5CV1YcW1hmsdOJg9mfPD7y2wHYqnw+YeO4M0vPjjpoYwdVlybQb/5sfx+t6s7tEJoHDmRnesyaig4ZYc5qFN5AU4XeXuqa/zIqy97truqsr5KukKwSy+DRtfmuAWqyFI0d+X+nlM0HSqrxxxY4ZpZMEGL57OYbcA9Ty3jwM4azjxlbtJDGTscuKacjFog2lJre211WGkx24mvPH4cr7ls+zoJXThwTYiB7fOOQsoUJRdZ5JXSDWFtACNKHHaZVYaL1Jddjyozr+WosDwF64R1jLuIgvmlwp+L1FX+uLxaK3ILFl47Z4LZgMrqeuyA5zLMKFFK4c7Hj+Nnb7xw0kPZFDhwjZih6rnclGFiaXcChBPAkrShG8CcbZkAhi7pQ6D3LKcNZs69101Ryg9A87jTfcMJAMHOheIg5f7cJSAVmS7Se6Q/No8e7xxmj07vPY/rBgcsZhN5drkJIsI5e7Z/mhDgwLXtSRyKu5cmOw4bqHqpNHOccj70gx0LoxsDlxow25T7n13B1WftSh3F2xwOXJuJU2ycbuqzuKRr1HCUV0cvQ9f6bXHMHM2jy/3VV36vKLbdD0oSrApUiipSJAUWebsWVrBjvnAsXdOBwPpVlt7Yua3X8QOeyzDj5NEjq3jpuadMehibBgeuGWIz1FcmBTjK67q1bes+l5UWs7154lgDl+ybbFZlM+HAtV5G1Pqp7z16Xb+fYcO1u+fnwpBTX71WXTYW+H7pBzdYFSopl465pAIV02MOq/F8qpyCXYtd57GA4VUWz2cx04BSCo8fq+O800aXVt/qcOCaUYZRMMD4lNW6xnDsZE8nIsPMAsfrIYKSwFJ1e61y3AsOXGOi69xVgWIrdCJ2cxgm+wdwGgKd815Ax9wX0Dt9aFXZhuue+qmqHsf1nMMCuhYPp5tGPI81wDUYZjN5bqWJAztrkx7GpsKBawvRNYC5HTUsPQwblkz6MNlIhUHDBrCRUxSggP5Bqui4UXSy4LQgs804utrGvp3VSQ9jU+FehQzDMFPM8Xobe5eCSQ9jU2HFNUZ6FiP3MHl0pBmLjs2nD11lkjNuuHSs7ZVcb8RKoot9vlBV9TpnWJXFaUFmm3OyEeKiM2bHUQiw4mIYhplqVlsRds+VJz2MTYUV10bpZ1kf8jp9DRuWXvNeLj1UmL7kaL7D9FRUQO9C5lEZKUahsvpch2G2CmvtGEu12XEUAqy4GIZhppq1VoTFYLYCFyuuTaBv4931zHe5x+fPyauJfi2bnMLivkppUAZtDTVq518fdcTzWcx2pdGOsRDM1kf5bL3aUbPOLhp9+xJ2ud5AJo+icRR9WA/Zf3BghgwUwwSpga6xwWsyzFajHs5e4OJUIcMwzBTTbMeolWcrcM3Wq90CDLReVx/DRtfz86qhSAlOoi3SKNJ4AyiidSmsAa/JMFudRhhjruJNehibCisuhmGYKUUqhXYkEfizFbhYcY2CDVjjhzFsuOdbBlJgeYa1829AsYzSJLFuhTXgdRlmWmhFEkHJg9jIiutTDCsuhmGYKUUHrtn7GGfFNWEGdhoCPRVSkfroed38tUfMOBx9rLAYJksrjFEtz1aaEODANTqGWGByIMOGe48B7zPIB33few54nYEYZ2pxyPswzDTSiiSqpdkLXLOnMRmGYbYJrUjOnBUeYMU1ejZDebn3cRniniNjCLUz1FhYZTEzSDOMsThDKx9bWHExDMNMKa1IYo7nuJiRMUT3+L6Fxr3u2Y+NWuBHrGhGovRYZTEzTjOMMV+ZvY9xVlwMwzBTiu6aMXuBa/Ze8WYyxHxXeokNqq/uFxz+Guu+5eRciQyznWmGMeZnrMEuwIFrcxhBANOX2UCt1iYycqOHvujor8kw24RmKGduLS6AU4UMwzBTzSymCocKXET0/xHR14noLiK6iYh2OPt+nogeIqJvENGrhx7pdmAM6kFJ1fOxGfcY2b2U7HwwDNOTWXQVDqu4bgFwiVLqMgAPAPh5ACCiiwC8GcDFAG4E8L+JaPbeXYZhmDEzi3NcQwUupdTfKKUi8+sXAew3P78WwIeVUi2l1KMAHgJw9TD32jZssqIYVC2NW7UVDIzVFcOMAE4VDsfbAXza/LwPwBPOvifNtg6I6B1EdDsR3f7888+PcDgMwzDbC/fzcnX5OABwHVcRRPQZIrqn4PFa55hfABAB+OP1DkAp9T6l1FVKqav27Nmz3tOnn1lQHKyuGGYkuJ+X80s7AQBz3KuwE6XUDb32E9GPAHgNgFcqpWxe6SkAB5zD9pttTC/W2f19S8JBiWE2lbnK7NkHhnUV3gjg3wL4HqVU3dn1CQBvJqIKEZ0F4DwAtw1zL4ZhGKYT7g6/fv4XgAqAW4gIAL6olHqXUupeIvozAPdBpxB/QikVD3mv2aKXcpmkGmNFxTBbitoM2uGHClxKqXN77HsvgPcOc32GYRimGDsvM4sLSc6extwODKp6+ikzVk8MM7VYR4HYQm3fNospdQAwDMPMNlKNqcZyCmDFNQC3yD+f9BAYhmEyqBkOXKy4GIZhppBxNbWZBjhwMQzDTCEzLLg4cDEMw0wjszzHxYGLYRhmCpnhuMWBi2EYZhpRmN3IxYGLYRhmCmHFxTAMw0wVHLgYhmGYqYJThQzDMMxUwYqLYRiGmSpmOG5x4GIYhplGuOUTwzAMM1XMbtjiwMUwDDOdzHDk4sDFMAwzhcxw3OLAxTAMM42wHZ5hGIaZLmY3bnHgYhiGYaYLDlwMwzBTyAwLLg5cDMMwzHTBgYthGIaZKjhwMQzDMFMFBy6GYZgphWjSI5gMHLgYhmGmlBmNWxy4GIZhmOmCAxfDMMyUQjOaK+TAxTAMM6XMZtjiwMUwDMNMGRy4GIZhphBfEPYuBZMexkTgwMUwDDOFnHnKHN731qsmPYyJwIGLYRhmCqmWPFx0xuKkhzEROHAxDMMwUwUHLoZhGGaq4MDFMAzDTBUcuBiGYZipggMXwzAMM1WQUltnHU0ieh7ANyc9jg1wCoAjkx7EmJmF1wjMxuuchdcITOfrPKKUunGQA4norwY9druxpQLXtEJEtyultnVBxSy8RmA2XucsvEZgdl7nLMKpQoZhGGaq4MDFMAzDTBUcuEbD+yY9gE1gFl4jMBuvcxZeIzA7r3Pm4DkuhmEYZqpgxcUwDMNMFRy4GIZhmKmCA9cQENGNRPQNInqIiH5u0uMZF0T0GBHdTUR3EtHtkx7PqCCiDxDRc0R0j7NtFxHdQkQPmuedkxzjsHR5jb9ERE+Zv+edRPSdkxzjsBDRASL6OyK6j4juJaKfMtu31d+SSeHAtUGIyAPw2wC+A8BFAH6AiC6a7KjGyiuUUpdvs7qYPwSQL+D8OQCfVUqdB+Cz5vdp5g/R+RoB4DfM3/NypdSnNnlMoyYC8DNKqYsAXAvgJ8y/xe32t2QMHLg2ztUAHlJKPaKUagP4MIDXTnhMzDpQSn0OwLHc5tcC+KD5+YMAvnczxzRqurzGbYVS6hml1B3m5xUA9wPYh232t2RSOHBtnH0AnnB+f9Js244oAH9DRF8hondMejBj5jSl1DPm52cBnDbJwYyRdxPRXSaVuG1SaER0JoArAHwJs/O3nDk4cDGD8FKl1Iug06I/QUQvm/SANgOla0W2Y73I/wFwDoDLATwD4NcmOpoRQUTzAD4K4KeVUifdfdv4bzmTcODaOE8BOOD8vt9s23YopZ4yz88BuAk6TbpdOUxEpwOAeX5uwuMZOUqpw0qpWCklAbwf2+DvSUQl6KD1x0qpvzCbt/3fclbhwLVxvgzgPCI6i4jKAN4M4BMTHtPIIaI5IlqwPwN4FYB7ep811XwCwNvMz28D8PEJjmUs2A9zw+sw5X9PIiIAvw/gfqXUrzu7tv3fclbhzhlDYGzEvwnAA/ABpdR7Jzui0UNEZ0OrLADwAXxou7xOIvoTAC+HXv7iMID/COBjAP4MwEHoJXbeqJSaWnNDl9f4cug0oQLwGIB3OnNBUwcRvRTAPwK4G4A0m/8d9DzXtvlbMikcuBiGYZipglOFDMMwzFTBgYthGIaZKjhwMQzDMFMFBy6GYRhmquDAxTAMw0wVHLgYhmGYqYIDF8MwDDNV/P9AefpXG5DfPwAAAABJRU5ErkJggg==\n",
"text/plain": [
"<Figure size 432x432 with 3 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"def fitGaussian(points):\n",
" m = np.mean(points, axis=0)\n",
" ms = points - m\n",
" c = ms.T @ ms / (len(points) - 1)\n",
" return scipy.stats.multivariate_normal(mean=m, cov=c), c\n",
"\n",
"def decompose(c):\n",
" S, e = scipy.linalg.eig(c)\n",
" return np.real(S), np.real(e)\n",
"\n",
"def computeBounds(S, e, n = 2.6):\n",
" s = np.sqrt(S)\n",
" ns = s * n\n",
" bounds = np.array([\n",
" [-ns[0], -ns[1]],\n",
" [ns[0], -ns[1]],\n",
" [-ns[0], ns[1]],\n",
" [ns[0], ns[1]]\n",
" ])\n",
" #print(bounds)\n",
" bounds = (bounds @ e).T\n",
" bounds[0] += m[0]\n",
" bounds[1] += m[1]\n",
" xbs = np.min(bounds[0]), np.max(bounds[0])\n",
" ybs = np.min(bounds[1]), np.max(bounds[1])\n",
" bounds = (xbs, ybs)\n",
" return bounds\n",
"\n",
"def jointPlotWithGaussian(points):\n",
" g, c = fitGaussian(points)\n",
" variances, rotations = decompose(c)\n",
" bounds = computeBounds(variances, rotations)\n",
" jp = seaborn.jointplot(x=points[:, 0], y=points[:, 1], kind=\"kde\", joint_kws={\"alpha\":0.1}, fill=True)\n",
" xr, yr = bounds\n",
" plotFunc(d.pdf, (*xr, 150), (*yr, 150), jp.ax_joint, zorder=-1)\n",
" return jp\n",
"\n",
"jointPlotWithGaussian(s)"
]
}
],
"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"
}
},
"nbformat": 4,
"nbformat_minor": 5
}
This is free and unencumbered software released into the public domain.
Anyone is free to copy, modify, publish, use, compile, sell, or distribute this software, either in source code form or as a compiled binary, for any purpose, commercial or non-commercial, and by any means.
In jurisdictions that recognize copyright laws, the author or authors of this software dedicate any and all copyright interest in the software to the public domain. We make this dedication for the benefit of the public at large and to the detriment of our heirs and successors. We intend this dedication to be an overt act of relinquishment in perpetuity of all present and future rights to this software under copyright law.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
For more information, please refer to <https://unlicense.org/>
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment