Skip to content

Instantly share code, notes, and snippets.

@michaelosthege
Created June 11, 2021 07:53
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save michaelosthege/730b1a20b1267791916c98ba4661762a to your computer and use it in GitHub Desktop.
Save michaelosthege/730b1a20b1267791916c98ba4661762a to your computer and use it in GitHub Desktop.
Display the source blob
Display the rendered blob
Raw
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"id": "b6c4287e",
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"from matplotlib import cm, pyplot\n",
"import pymc3 as pm"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "9d511476",
"metadata": {},
"outputs": [],
"source": [
"D = 2\n",
"BOUNDS = [\n",
" (-2.5, 2.5), # vertical\n",
" (-3.5, 3.5), # horizontal\n",
"]\n",
"SIZES = (6, 8)\n",
"N = np.prod(SIZES)\n",
"\n",
"# create grid coordinates corresponding to well positions\n",
"TRUTH = {\n",
" \"optimum_x\": np.array([1.5, 2])\n",
"}\n",
"X = np.empty((48, 2))\n",
"Y = np.empty(N)\n",
"for ir, r in enumerate(np.linspace(*BOUNDS[0], SIZES[0])):\n",
" for ic, c in enumerate(np.linspace(*BOUNDS[1], SIZES[1])):\n",
" i = ir * SIZES[1] + ic\n",
" X[i] = (r, c)\n",
"Y[:] = 100 * 1 / np.sum(np.exp(np.abs(X - TRUTH[\"optimum_x\"])), axis=1)\n",
"\n",
"# create coordinates for cross-sections\n",
"X_denses = {}\n",
"for d, (lb, ub) in enumerate(BOUNDS):\n",
" x_dense = np.zeros((30, D))\n",
" x_dense[:,d] = np.linspace(lb, ub, 30)\n",
" X_denses[f'D{d}'] = x_dense\n",
"del x_dense"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "8563c180",
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAfMAAAF3CAYAAACxLqKFAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAACaj0lEQVR4nOydeXwdVfn/32dm7pY9TdMtKd3pCmUpa9lENkF2BFF2EFAQARUQ9y9uPxVFAVHcRVSU1cpWQbayt+xQum9pmzZJm/1uM3N+f0xuSZub5Ca5ZyZJz/v1qpLcyTnz3DlzPs/znE1IKdFoNBqNRjN0MYK+AY1Go9FoNANDi7lGo9FoNEMcLeYajUaj0QxxtJhrNBqNRjPE0WKu0Wg0Gs0QR4u5RqPRaDRDnMDEXAgRFUK8JoR4WwjxvhDiu0Hdi0aj0Wg0QxkR1DpzIYQACqWUrUKIELAI+JKU8pVAbkij0Wg0miGKFVTF0vMiWjt+DHX80zvYaDQajUbTRwIdMxdCmEKIt4CtwH+llK8GeT8ajUaj0QxFAovMAaSUDrCPEKIMeEgIMUdK+V7na4QQlwOXAxQWFu4/Y8YM/29Uo9FoNJoAWLJkSb2UsrK36wIbM98VIcS3gTYp5U+7u2bevHly8eLFPt6VRqPRaDTBIYRYIqWc19t1Qc5mr+yIyBFCxIBjgA+Duh+NRqPRaIYqQabZxwJ/FkKYeE7FP6WU/wnwfjQajUajGZIEOZv9HWDfoOrXaDQajWa4oHeA02g0Go1miKPFXKPRaDSaIY4Wc41Go9FohjhazDUajUajGeJoMddoNBqNZogT6A5wgw0pJa4Ed5d9dAzh/fPOhhmaSCmREhwJUn60Cb4ATGN42Oe4kHYlbqcHaBoCyxCYxtC2z5WSRNolkXZxXO8BCiEIW4JYyMAyxZC2L2W7NLSnaIynSXfYZ5mC0miIisIQUcsM+hYHxLa2FCvq2tjYlCBpuxgCCsImkyoKmTKygFho6NonpWRzc5KVdW1sbUmSdiWmISiNWkwZWciEETFC5tCNG10paYqnaWhPE087uFJiCkFB2KSiIExJ1BoU795uL+aulKQdT+R6QyCxDLCGiDBIKXEkpGxP6HrDEJKwJYaUfUlb0p5ysV1J5o67e5QhU1AQNggPEeFzXElTe5qGNpuULTHEzo4YsON3QkBx1KSiKEQsPDSEIZ52WLa1leX1rSRtF9MQOK7c4Uwbgh2/C5kGUyoKmDmqmKLI0Oi2NmyP89A7tTy/soGU4xIyBClH4rgSIcAyPCcsabuMLAxx8pwxHDezkuIhYJ+Ukvc3t/DA25t5s6YZ8J5X2nFxXO+/Q6aBEJByJBPKY5y5z1gOmzyCsDX4hd2Vkg2Ncd6vbaGhPYXR0V84rkSSCYLEjmtHFUWYM6aYcSXRwPqWQbOday7kcztXV0pSTtcoPFdMAWFz8Iqe7UgSadnvY+jCJoStwSl6UkriKZe2lNsv+wwBhWGDaMgYlPa5rmRLc4rtbTbQt6MEhYCwKRhbFqEwMjhFvT3l8Or6bdQ0JRDk5kiD99wARhdFOHjCCEqig1P0Vte38cvn1rCmoR3HlTnbFzENJJIjplZw+fwJg1bUX1u3nV+9sJamuE3Szv0djIU8ET9j7lg+vd84rEEYrUsp+XBrK29tasKVYPdBIKyOLOC86jImVxTkrW/JdTvX3U7MpZSkXbBziFRzIWJ+5KENBqT0RDwf9gkBsZAYVPbZrqQ5bufFvpApKImag8q+tqRDzbbkjgigvwgBZQUWY0rCGIPEPiklqxvaeGV944DsE4BhCParKmXmqKJB45DZjsvflmzkgbdqSTv9czQBQoYgEjL46sencOCE8rze40BoSdrc/twaXlvXSHIAL2DEMhhZGObrx09jUkVBHu9wYDQnbJ5bXU9zwu6TiO+KZQgqC8McNqmCgjxkybSYZ0FKScLO/6HploDQIIjSHVfSnsr/84xYDIrUWNJ2aYo7eS1TAGUF5qAY06trSVHXkiZfr6TAG3eeVBkN3D5XSl5Y3UBNU2JAHWVnLENQURjmmKkjA4/yWhI2Nz7yAZuakwMSus5ELIMTZlZyxfwJgfct67fH+erDHxBPOd6chgEi8PqUq4+YyDHTez0QTDk1jXGeW90wYCc6g8Brn8fuWUllUWRgZQ32g1b8RpWQA9gSUo5XR1CoEnKApO1NUAqSZDr/Qg5ee9je7gRu39bm/Ao5eLalHcmqrXHSAdrnSsn/VtazoTF/Qg5elqa+NckTy7Ziu8HZ15Kw+dID77GhMZE3IQfPeX1iaR2/eHZNoH3Lum3tXPfg+7Qk7LwIOXhtM2m73P78Wh7/YGteyuwv67fHeXZVA3aehBw63j1XsnB5HVtbknkqtWd2CzFXKeQZHAnpgPoTV6oT8gxJ2xuHD4K049KUyL+Qd6Yp7uRVaPrCttY09a35FfLOOC6sqU94s+AD4KW126htTuIoMNCR0BhP88zK+kAEz3Zcbvr3UupaU0raT9J2eXZlA399vSbvZefCtvYUX314KfGUo6T/TNkuv160jlfXbldQeu/UtSZ5fk2DkrYJnsP53xV1NCXSSsrvzG4h5mlXrZBnsF1PWP1ESklCsZBniKel7x2mlFJJRN6lHqApbvtuX8p2qW1OKRPyDLYjqW1Kqa0kCzVNcdZuiyvrLMET9K2tKVbWtymrozv+8cYmNuZx6CAbSdvl/rc2s6LOX/uklNz69GriaVtp/5lyXH7y9CqafRC8ztiuy7OrGpQ7ubYreW5Vg3JtGPZi7sr8TAbLlaTtb7o97eQ+WzYfxNP+il1L0un3ioO+4rjQnvKvsUgp2bAtqVzIocNZabdpS6p3jDKkbJcXVm9TKuQZbFfy2oZG2lK28royrGlo519vbsprar07Uo7kB0+uIJXLGtM88cyKej6obfGl/0zaLrc9s0Z9RZ1YvKHJl2cH3uTB92tblNYx7MU85V/fBXidpl/Og7fO2p+6MjguvqVrHdebme8nbSnXt+xKS8LxrTMBr21u3J70zdl8Z3OTr2PZjitZvKHRt/rueH4NaR896e3xNE/4NL5sOy6/emEdCZ/ap+1K3qhp4sMtrb7U15K0WVnf6oujCZ59b29qUjo3Z1iLuSulb1FdZ2zXn+g8qDlNKdufL7Xdb0+sg7hP0Xl9nie85YLtSl/sc1zJ8ro2X98/CWxojJO01bebjY0JVta1+TJ8lyGTbvejb3lpzXbfhwxTjsv9b23ypa6lW1p8fXYACMHKBnXOyrAWcx/e6axI+r8ZTV/wS1R3xY+5AZn18kEQT7vKO8xkx9asfiMl1LeqH5tct73d/84SAMEKH8bOH3l3s+9iB15E+e5mtelagH++uYm4z+1TSnhtXSNNcbXt03ElK+r9dTQz9b5f26qsbxmcWwzliVy0buGTT3DD9dfiuA4XXnwpX7nhpqzXLVn8Okcddgh/ufcfnH7mWb3X7Xp7nqtC5ph1+MIVl/LE449SWTmKV5e80+Xz5cs+5POXX8rbb73Bt77zPa657ss51e+4YCjcYCyX9OXGmg1cfcUl1G2pxTAMzrvoMi7/whd3uqa5qYkvfO5CNtZswLFtPn/N9Zx73oU9lutKzz6V24E3x3ObVNTc1Mi3v3oVK5d9AEJwy613sc/+B+34/A933cajD90HgOPYrF6xjBfeXktp+Yhuy2xNOEgpla5dXtXQlnVS2O//7yu8tehpSsor+P59TwHwwF0/5c3nFyKEQcmICi779q2UV47p8rcL//57nnv470gpOfK0czn+M5d1ucaRklX1bcwZU5J/ozrx/MptXTJj8W1bePfP3yXZ3IAwDKrnn8bEo8/hwwdvp+7dRQjToqCymr3O/wahguIuZa79333UvPgIIKmefyoTj/50l2uSaZdnV9Sz9zh19jUn0qzbFu/yez/sM4XgtXWNHDtD3drzLS1JsrV81W0TvOxDYyJNeSycT5OAYRyZ5+L9OI7D9V+6mocWPMaSt9/nX/f9g6UffJD1um/cfBPHHHd8zvWr9vpynQfz2fMv5MFHHuv28/LyEfz41tu45trcRPyj+tUamHZ6X/NpWRbf/f6PWbT4XR57ehF//O1dLPtw5+f3h9/exfQZM3nmpSU8+NhTfOfmG0ilep7VLSBv62m7oy3HIYQfffsG5h91LAuee5MHF77C5KnTd/r8ks9fywMLX+aBhS9z7U3fZd7Bh/Uo5ODtDpdUnNVpaM8eXR32yU/x5V/+ZaffnXj+FXzv7wu55W9PsM9hH+eR3/2iy9/VrFzGcw//nW/9eQG3/O1J3l70NLXrs0+Yak7aSqPm5kQ66/MTpsn0M6/h8G/fx8Ff/R3rn7+f1s1rGDnjQOZ/414O+8a9FI4az+on/9zlb1s2raLmxUc45MY/cOjN91D37iLatq7vcp0EPtisdlx5RV1b1k2i/LAvYbss3aI281DflszqaPrRNgEa2tRkHoatmOfSFy9+/TUmT5nKpMmTCYfDnHX2OfxnwSNdrrvrzts57fQzqKwclXP9ErXj5rmK6fzDjqB8RPede+WoUew/7wCsUKiP9ffp8j6TyiEyHz1mLHvvsy8ARcXFTJs+g9pNO4+5CSFobfFSW22trZSVj8Cyek5IeZutqDUwlxR7a0szS159kTPP9TIJoXCYktKybq9/7OF/ceKpn8qxfnVjUAnb6fb7m77fQRSWlO30u1jRR1FcMt6eNWOwae0Kpuy1H5FoDNOymL7fwbzx7BNZ6zCEoCmhbmZod2IXLR1J6R4zALCihRSNmUiicSsjZx2EYXptrmzSHBKNXSextdWupWzSbMxwFMO0KJ+2H1veei5r/Rub1O4ZsGJrG6ksY5R+2afaWaltTWYNFPxom7YrqWtVs4nMbi3mmzZupLq6esfPVVXVbN60scs1Cx55mMsuv7LP96Ay9glo/w/f6u9rZ7V+3Vree+dt9pt34E6/v/TyL7B8+YfsvecEjjpkP773/27FMHpv9io3yMkc19obNevXUj5iJN+4/krOOv5QvvWVq2hvzz4eHI+3s+jZpzj2xFN7LdeVXrpWFS0JG7OPKfz7f/Vjrj/pIF5+4mFOv6Jrlqh6ynSWvfkqrY3bSSbivPPSMzRs2Zy1LAFK1yxvakpg9/IA2xs20bxhOWUT5+z0+5qXFlA565Au1xeNncy2lW+Ram3CSSWoe/8lEtu3ZC3bNAQNber2DFjd0N7r5FqV9m1RvGNacx8dvXy2TYBtiuYEDFsxz4VskfOuntcNX76OW37wI0xzcJ1ANXR21FdPW2srl55/Drf86KcUl+w8lvjM0wuZs9dc3lm+jv8tep2vffVaWpqbey1zMDhitm2z9L23OOf8y7j/yZeIFRTw+ztvzXrts/99jH0POLjXFHsGlSuq+rOByllfuIGfPfoqh5xwGk//809dPh83aRonXvB5fnL1Z7n1mvMZP21mt++kRO0wUNKWPabx7UQ7b939NWacdS1WrHDH71c9/keEaTH2wBO6/E3R2ElMPvZ8Ft/+RRbfcS0lVdMQZvYMkiFQut68t6yNavtUD+H1tfx8ts3+1J8ru7WYV1VXU1Pz0TaJGzfWMGbsuJ2ueeONxVx43rnMnDaJhx+8n2uvuYoFjzzs851quiOdTnPJeedw5tnnctIpp3f5/B9//QsnnXIaQggmTZnKHhMmsmL5sl7LVXmsRa5ljxlbxeixVey93wEAHHfSaXzw7ttZr338kftzTrH35R76gxD9r+DgE05j8f8ez/rZkad+mu/+9TFuvvt+ikrKGD1+Uvb6UXvokSkytXTFdWze/O3XGHvg8YzZ92M7fr/xlUfZ+t6LzL34u93eW/X8Uzj0a3/hoOt/TaiwhMLK6qzXSdhxvrYKejpF0A/7VJ8p09/y89E2B1J/bwxbMc/lC9t/3gGsWrmCtWvWkEqluP+f93HSJ0/Z6ZoPlq9m6Yo1LF2xhtPOOIvbfnknJ596Wm730I/7zpVh++A6yOW7k1Jy3VWXM236DK68+tqs11SNH88Lz/4PgK1bt7BqxXImTOr+RdtRv8IeJdeiR44azZhxVaxZtRyAVxY9y5RpM7pc19LcxOJXXuRjx5+U8z2oPPY1bBp9Sm10niz05vP/ZezEKVmva95WD0BD7UYWP/MEBx9/StbrvHtQZ19B2MTK8v1JKXnvnu9TNGYikz7+mR2/r3v/ZVYvvIf9r/wJZjjabbnJlm0AxLfVsuWtZxl7wHFZr3McmZejNbuju3Pi/bIvonIZCfTpBEE1bVNN7z1sl6bl0ldZlsWtt93OqSedgOM6XHDhxcyaPZvf3f1rgH6Nk3dGpSAYhsgpX3vxBZ9h0QvP0VBfz4wpe3DzN79NOu2N2Vz6uSvZUlvLkfMPpKWlGcMw+NUdv+C1N9+jpKTnpS8K+0rAO7rT6WXG9WuvvMS//nEvM2fP4ej53gmBN3/rFjbWbADgwksv5/obbuaaKy/jyIP3RUrJN7/7fSoqRvZaf0ihgUIIQqbIafndzbfcyo1fvJR0KsX4CZO45da7uO+e3wFwzvne8penn1jAoUceTUFBYU9F7cAQEAupcwdLoyHsbtLQd339aj5c8jKtjdu57qQDOe3y63nnxWeoXbcKYRhUjKnioq/9EIDtdbX88Xs3cv0vvNnRd9x4Ba1N2zGtEBfccEuXyUoZHCmVLP3JMHlkYVaHrHHV22x67XGKxk3hxR+cD8Cep3yepf/6GW46xeu3XwNA2cQ5zP7MjSQa63jv3h8w76qfA/DW3V8j1daEYVrMOucrhAqyv4Nhy6As1rcJq31h+qginlu5rcvuhH7ZN3FETJltACMLwlnHzf1omwIYNcAjUbtj2J5nLqUk7vNWp50xBEQtdYKg8sjTXAibEFEoCPGUS4uP+4h3RgAlMZOIwjPca7YlaQqogQpg2piY0jPOH3x3U2DPL2wKzt03ewo3H6QclzN+tziwU+j2GlvMj0+bpaz8FXVt3PjIB75vGgNev/np/ao4/0B1z295XSuvrW/0bSvXzoQMweGTKxhflrvDstufZy6EUJrm7g2FWUxfyu8NlWlaUBsZ94bEe+lUUhgxAmufhtG3VGN/GFMcCcy+ikJ1UTl4adJxpWqiq94IGYJ9q0uV1jFhRCyw44CjlsGsMUVK6xhZGFY+Lt8djpSMVNQ+h62YAygMrAKvWwgRqH2KtQDTUJ/K746QKbxhDIWUxIIZ4RJAeaG6FG2G6ZXFyr/DbFiGYOaorruP5ZtT9xpDNIgXUMAxM3ofJhoIYdPg8CkjAgkYTEOwj2JnpTwWIhYKZnXSqKKIsrq1mCvAEGpnm2YIK0zj91ivqXY+AHjlq5zk0xMFYfUNxzQEpQEJ+ohC9fVWFIYpCuD5mYagqrT7SVj54ug9RwayPHTO2GIqFY25duaMuWOVZ292JWwKTtlrjPKsnxCCvcYUZ53EqBLLEMwZo87RHNZiLoQgCL1TOJS8E6YhAvGeQz59qdGQ/0MlhlA7E7ozI4tDvttXFDV966Tnjiv1tcO0DMGc0cW+ONKxkMkJMyt9aysAEcvg3P2rfKlryshCJoyI+ZqONoTgpNmjfalr0ogC3/vOiGUwrkSdozmsxRzA72yKKdSPJ3cmFvK3RUYsf7IO4DljJVF/H2BpzFSedcgQCRlUFFm+CbohYFyZf2O9E8tjVBT4Nz5ZGDaZpTDy2ZWLDhpPYdif7ErIFBw6qZy9FB6wsitf/fgUwjnslpgPIpbBJYfsQXmB+iEgAMs0mD+xwjdn0xSCIydXqF3hpKzkQYIQgoiPeuB3ZtEwhG91GsL/iWmRkOFb9BMLCd9Ti5UlYSwf7BMdQu5HXR/VKTh88og+b+3aH0whOHLKSN8cTYBoyOSmY6cqXfWQIRYyuerwicrr6Ux1WYzzDqhSbp8hYFJFjE/Ozv3si3ywR3mMcSVR5RG6KWBaZaHy4ZFhL+bgRcp+jJ9HfBhLzkbYEsonpAkgFhaB2FcSM5W/cJYBRX56fR0YQjChQm2HIoDSmEVJzH/7CsMWR0yuUCrophAcuEcZ5QrXXnfH3lUlnLXPWKWCF7EMvvOJPSmM+D/H4vS5Y9lrXLEyh1rgbVLz9eP2DKRvmT9xBAUhU92ubALKYiHmVZepqaBzXcprGCSEDJSOn0dMf9PrnRFCEAsJpbO/C8LC16inM4YQlBdYygTPNKCswAqkMwEv+zBxpBpBFwKKYybjysKB2Te+LMahE8uVCLopBPtWlbJnpdrlTD3x2XlVnDR7lBJBj1gG3/7Ensz0cfigM6Yh+OYJXv35ts8QnpDfevpsRhapXU7YHWHL4MSZoykImXkXQ0N4GygdN32UL9qw24i5EIKwJZRMTgtSyDMIIYiF85+BMAQURtQv1eoN0xCMKLTyPsYVNgUjCqzAHJUMsbDJpMoYITN/k/4E3sz16vJIYEKeYXJFIUdNrfCW/eXhVgTehLdDJpQzOyCh23EvQvC5QydwwYHVRKz87B8QNgWlUYsfnjxD+bry3u/F4HsnTeeIKSPyJuhRy2B8WYxfnjWHcT6sPuiJWMjkk7NGU1EYzlv/YhmCcSVRPjFjlLLtW3dl2O4A1xOulCTtgZ+MZQp/lmn1FduRxNMDf65h00vhDyb7pJTE0y6tyYHtTiWA4qhJ1K+lBzniSsmWphTb2+x+t0+Btx1u9YhIYMv7uiORdnhx7TZqW5L93pjEMgTlsRBHTK6gKIDUc0+s29bODxauZGtLkkRv54hmQeBFi4dNHsHnD5/g2wS7XHljQxP/76mVJNIOqX4cvefNuzH49H7j+NS+4wIPgjojpWRZXStLappwXUl/ehijYwL0oRNGMHFEQV7uK9cd4HZLMYeOM6UlpJ2+i7ohvLT9YGqIuyKlJGXLftkXMjwRDzoa7wlvO1unz06LwJvoVhA2B7V98ZRDfUualoS3JWouVgrhCV1FkUV5YSjwbENPbGyK887mFhrakkh6P2ZA4A23lMYs9hpTwoTy2KByMjvjuJL/La/nn29uoq41Rcp2e31+mYmlc6tK+PR+VcweG2y2oSfaUjb/eW8rD7+zmaTt5rTta9QycCUcNW0En9p3HNV92M7Ub9pSNu/XtrCivg3I7Uhfq2OZ8PTKImaNLiaax2VUWsz7gONKHNc747m7b8PbCMabKDWYO8ldkdKzze6wMVu7FHTsuGYIQoMw09ATUkoStiRlu6QdmdU+s2MWftgyiAyyTENv2I6kqT1Na9IlkXaxXfnRZJ0OW8OWoCBiUhqzKAgbQ8q+lqTNqoY2apuTbI+nsF254/1ypfff5bEQY4ojTKooDGSS20BYtrWV51c28N6mFtZtj+NKiSkEEu+dLAqbTBlZyD7VJRy950jlW9HmE8eVLNnQxCtrt/H+5lY2NSV2bJgl8dpueUGIPUcVMm+PMo6cWjHoMkU9Ybsu67bHqWlMUN+WpC3tYCAQAqQEiaQobFFZFGZ8WYzxpTElAYIW836S+T4y34pgaIlbb0gpPdskIIanfZ0F3RDDyz7HlbhSIqXXaZrG8LIvabvYrouU3lBBxBxazklPSClpStgkbRdDQEHIDGSGuiocV9KUSJO2JZYpKAybeY1Qg8ZxJUnb3eFkRizDl+xsrmI+fFpSnsh0HMOj++jKjgNohqmBQqid1R80piEwh+vDw5u9HRmm83KFEEqPLg0a0xCMKBg6mYW+YhrBbTGdC8PzrdFoNBqNZjdCi7lGo9FoNEMcLeYajUaj0QxxtJhrNBqNRjPE0WKu0Wg0Gs0QR4u5RqPRaDRDnMDEXAgxXgjxjBBiqRDifSHEl4K6F41Go9FohjJBrjO3gS9LKd8QQhQDS4QQ/5VSfhDgPWk0Go1GM+QILDKXUm6WUr7R8d8twFKgKqj70Wg0Go1mqDIoxsyFEBOBfYFXs3x2uRBisRBicV1dne/3ptFoNBrNYCdwMRdCFAEPANdKKZt3/VxKebeUcp6Ucl5lZaX/N6jRaDQazSAnUDEXQoTwhPxeKeWDQd6LRqPRaDRDlSBnswvg98BSKeXPgroPjUaj0WiGOkFG5vOB84GjhRBvdfw7McD70Wg0Go1mSBLY0jQp5SKG7UGcGo1Go9H4R+AT4DQajUaj0QwMLeYajUaj0QxxtJhrNBqNRjPE0WKu0Wg0Gs0QJ8i92QclUkrv/zt+FoC3im54IKX0bJOAGJ72ufKjnw0xvOxzXInTYaAQYBliWNnXnnJI2i4SCJuCwrA5bOyTUlK3vZV4Io1hCEoKo5QWx4K+rbxh2w5121tJpW0sy6S8uICCWDjo28obacelPeVgS4llGBSGTSxj8LRNLeZkOkhw5EcivjMSQ3jCYBlgDKHORUrPNtuR2BJkFgMFEsPwhCFkDi3xk1KSSEuStkva2VnIM5gGhAxBJGQQsYaW+NmOy/a2NC0Jh3jKxXYlosMJyzhlkZBBYcSkrMCiMDK0xK+hLcWSmmZW1LezuTlJynExOzpIR4IlYExJhCkVBexfXcKY4kjAd9w3Fr+/jn89sZhFb6xi6arNSCkxTC8hmk47lBRFmTu9mqMPnsG5Jx7AuFFlwd5wH3Acl/++tJR/P/M2L7+1mtU1dZiGgWEIpIRU2mZ0RQn7z96D4+bP4lPH709J0dBxXtKOy9ubWli6tY112+M0JWxMQ3jvngRHSkYUhJhYHmPW6CLmjCna0XaDQMhsvfsgZd68eXLx4sV5KUtKiSMh7XQn4N1jCAgZBPrgekNKScqW/bLPMiBiCYxBbJ/jStqSDom07JN9AoiFvYhvMNsXTzlsaUrRHLdBZHfCdsUQXpusLA5RURwe1E7nh1vbeGpFAxsaEzvexZ4w8GwbVRTm49Mq2Hts0aB1Wmzb4e+Pvc5P/rCQjVsbSSTSuL08wEjYi6uOPGBPbrzseA7dZ4oft9ovmlvj/Pofz3P7354hkUzT2p7s9W8KYmGkK/nU8fvzlUuOZdqE0T7caf/YHk/z7KptvLa+CYBUb40TiJhekDB/YhlHTC6nKJK/OFkIsURKOa/X63ZHMXelJGn3XeR2xRQQHoSRrO1I4umBP9ewCeFBFslKKYmnXFqS7oDKEUBJzCQaGlzTRlwp2bw9SUNbOicBz0Ym/T5xZIyCiJnfGxwgLUmb+96qZVVDe06dZDbCpmBsSYTz9hvHiIJQnu9wYCxdtZnP3vgH1m9qoC2e6lcZsWiIM47Zl5/d+KlBF8k+9fJSLrr5T7TFUySS6T7/vWkahC2TGy87ni9fdCyWNXjap5SSl9Y18p8P6rxsbT+ap2UITENwztwxzB1XnJf70mLeDWlHkh6YDnQhYg6OKD2TcrbzaJ8hvEh2MER5ritpbHew3b5F4z0RtgRlscGRmk6kHFbXxbGd/NgnBIwsCjG2LDIo7Fu6pZV73tiM7bj96ig7Y3Q4LGfuNZp540vzc4MD5Jd/fZrv3PEfEimbgfar0bBFYUGE+2+7goPnTs7THfafdNrhC7f8jQf++wbxRN9FfFcKYmEmjKvg4ds/zx5jR+ThDgdGW8rhD6/VsKk52W8nszMhUzBtZAHn7TeOiDWwgCFXMR9cYYlCdqSd8yzkAEmHHZOSgsKLWPMr5ACuhLakxA3YPseVNLTZpPMo5AApW7Ktze41Daqa9pTDii3tpPMk5OCl5hta06xvSAxYXAbKGzVN/HnxJpL2wIUcvHaZciT3v7uF51ZtG3iBA0BKyU0/e5Dv3vko8WQ6L991ImXT0NjGSVfewf9e+TAPd9l/kqk0p1x1Z96EHKA9nmL5mi3M/8yPWbU+2KOtmxM2P39+LRsaE3kRcvCCxuV17dzx4noSaScvZfbGbiPmaRdshf1ZkIIupZdWz1M7zEp7SgYmeK6UbG+zs05uywe2C9vbBh5N9ZdE2mHVlnYl9rkSmuI2NduSgdn3fm0r/3x7C2kFBqYdyeMf1vPKusa8l50r3//NY/z2X4toT/Qvrd4T7YkUn7rubl55e3Xey84Fx3E55/rf8uo7a/Im5DvKdl0amtr42MW3snFrY17LzpX2lMMdL66nKWHnvf+0XcnW1iS/fnkDaUdBFLkLu4WYO27+I9ZsJB0C6TBTtjdjXSUSiKdkIPY1xx2ljgp4gt6S8MeD7owrJWu2xpU5KuBF6Nvb0zS22+oq6YameJq/vrFJiZBnSLuSh9/byubm3idi5ZvnFy/nZ396SomQZ2hPpDjrS7+hqSWurI7u+OU9/+OFJSuJ92N8PBeklGxrauczX/ldIH3LfW/X0phIKw0UaltSPLpUffZh2Iu5lJKkj310ymc9cF3pW52u9CIhP0mkXZIqUyqdiKclKT+8vk7UNiaVCl0GKaFmW8KXCOGjOiX3vrEZ2wf7bFfylyWbfM2OtcWTXHDTH5UJXWda40mu/eF9yuvpzIp1W7jl148qdVTAi/7fW7mJu//5gtJ6duXdzS0sq2tTHgilXckr65pYt12tMzbsxdyn4YodONLfdHs+Zq33haSNb+l2KSXNcX8fYFPc8S1CSKQd6lv7P2u9r7gdgu4Xb21qYUNTQmnWIYMEGuNpnl+9XX1lHXzr9n/T1OrP95lM2fz7mXdYtGSlL/UBXPz1v5BM+5PNaY+nuPkXD7OlodmX+lK2y31v1/oWnKQ7nE2VfeewFnMppdJx8u5QMckuG46bfZMU1aR9+lL7uoY8H2QmVvnB1uaUb0KeoSXukPYp+7BweYNv3yV4z+2ZVdt8cTZb25P86aGX+rU8q7+0J1L88LeP+1LX28tqWLpqs68TX13H5Xf/WuRLXW9uavZ9jlM85bC8rl1Z+cNazH3OmO7Alf5Er6kgPBW8oQTV0auU3qYwQdA2wDXsueC4MpAxbID6VvUCtKExwfa4f0KXwXZcPtzapryevz/6WiDL/V58cxUbatVnH35xz9MkU/62z0TK5s5/PIttq33vpZT8b+U2Xx1NgKQj+d/KBmXlazEfonVL6c+kvu5QPc5kd2yvGwRpRyr32psCEnIJNLSqHQMFeGltI3YADzDpSF5Yo17sfvX3Z/u9KcyAkPDXBa8orSKZSvPgf9/Ecf3vYGzb5ZnXlimtY3NLkuZEMO/fuu0JWpJq6h62e7PvOFCkFxY++QQ3XH8tjutw4cWX8pUbbsp63ZLFr3PUYYfwl3v/welnntVruaozOLmW/4UrLuWJxx+lsnIUry55p8vny5d9yOcvv5S333qDb33ne1xz3ZdzKtdxJZapLjLJZSxrY80GvnjlJdRtqUUYBudfdBmf+/wXd7qmcft2rrv6c6xds5pIJMrP77ybmbPm9FiuwNtFT+VGQC0Ju9cU+9pVK7jh6ot2/Lxx/Vo+f/3NnHfpVTtd9/rLL/CT/7sJO52mfEQFv/9nz6lY1/X2nQ6Z6nz5VQ3tWd+/J39xM6sXP0tBaQUX3rEAgK2rl/LUr76Dk05imCZHX/ltxu65905/11K3mcdvu5H27fUIYbDX8Wez3ykXZK17w3a149jJVJqVG7rOTpauTWrlQ+A6gItROoXQ2INIrX0SmfAcDOmkEGaYyIxPd/l7e+tbONs+AAQiWkFoj6MRxs5ddDJt88yry/ja5z6hwjQA3l+5mXDI6hKZ+2FfeyLFK2+v5thDZ6kyj3XbE4G1TcsQbGhMMGt0Ub7NGsZinsM1juNw/ZeuZsFjC6mqrubwQw7kpE+ewsxZs7pc942bb+KY447PuX7VYp5rZPzZ8y/k8iuv4orLLsr6eXn5CH586208uuARJfX3l1zGdS3L4jvf+zF777MvrS0tHHfkQRzxsY8zfcZHz+8Xt/4/Zu81lz/eez8rln/I1778Je5f8GSP5Uo8ZyKicKfQ9hyWIEycMo1/Pv4i4LXB4w6aztHHn7zTNc1NjfzwG9dz518eZGzVeLbV974ERgiIp1xCMTVibjtutyn22R8/nX0++Vme+PlHTvMLf/oJh5x7FZP2P4LVi5/jhT/9hLN/cM/O92yaHHnJjYyeMptUeyt/vf5MJuxzKBV7TO1SR9qVNCdsSqJqurf3V24mFgnTYu/iNAiT8JRTEWYYKR1SKx7ELZlAeOJH/UZ64yKE2fWwGJlqxal/h/CMzyAMi9TaJ3C2r8CqmNnl2neWbcy7TZ1544P1OE6W9umDfY7j8oLiSX6rG9qzBgt+tM2k7bJ+uxoxH7Zp9lzEdPHrrzF5ylQmTZ5MOBzmrLPP4T9ZRO2uO2/ntNPPoLJyVJ/uQeW4cq4TU+YfdgTlI7rfLrFy1Cj2n3cAVqhvyqU6g5rLcq3RY8ay9z77AlBUXMy06TOo3bRpp2uWL1vK4UceDcC0PWewYf066rZu6b1+hQZmdiPsC6+++CzVe0xiXPUeO/3+8Uf+xdEnnMzYqvEAjBhZ2WtZrszNmegvW1tT3Ub91XMOIFq0y/arQpBqbwUg1dZC4Yiu71nRiFGMnjIbgHBBERXVU2htyP4cLUOwSeGa83eXb8yaghZCIMyOIz+l6/3rhJQSp3EVRvm0rOVKKcG1kdIF10aECrNeF0+mqdvWMjAjemDxe2tpz7JBjF/2vb9y88AM6IWaxuxtw4+2KYG1ipaoDd/IPIe+ctPGjVRXV+/4uaqqmsWvv9rlmgWPPMxjC59myeLL+nYPeClbFQQ4XO4LffWD1q9by3vvvM1+8w7c6fez5+zFYwse5qBD5vPGktep2bCOTRs3Ujmq51ObVE5g7E/W5sl/P8AnTuk6vLNuzUrsdJpLzzmR9tZWPnPJlZx85md6Lc9R6KzEbbdP7f6oy27mwW9fxnN//DHSdTn3x3/v8fqmLTVsXb2UMdPnZv1cAnGFa1KbWuPdTtKS0iW17J/IVBPmyL0wCsd89FnbZoQVw4iUdfk7ES7CGrUPyQ/+DMLCKBmPWbJHl+sAQpZJS1uCyhH5OchjV+obu59A6Id9qte1J/owwS7fbROgXVHbHLaReS5ki5x3naF6w5ev45Yf/AjTHDyn+2h2pq21lcvOP4f/++FPKS4p2emzL153A42N2/n4YfP4w2/uZM7e+wyKk5r6InbpVIrnnnqMY086vctnjm2z9L23uOOP/+JX9zzE3b/8MetWr+i1TJWJlb76QW8//neOvOwmLv/Dsxx12ddYePs3ur02FW9jwY+u4ajLvkakoJtUpZRKl/w5rttt+UIYRGZ8msisi5DtW3HjH81edrYvx+wuarUTuE1riMy6gMici8CxcbZlnwgmBEonp/VUth/2qV4O15fS89426fv7kSu7tZhXVVdTU1Oz4+eNG2sYM3bcTte88cZiLjzvXGZOm8TDD97PtddcxYJHHvb5TrsS/BlYg4N0Os2l55/DGWefy0mndBW74pISfvGr3/H0osXc/ps/0tBQzx4TJvVarspVR4boW4ey6Nn/MmPOXCqyDPOMHlvFoUceQ6ygkPIRFex/4HyWLX2v1zIVzn0j1MeJkR/872GmHXIcAHvOP4Ha5V0nagI4dpoFP7qGmUeezLRDj+u2PCFEn++hL8QiIcxevkBhRTCKxuG2rAe8iNZpWo1Zll3s3NYaRLgEYcUQwsQsm4zbVpv9WtclElY3oaMw1nXMe1dU2hdS7GxbfZjYmu+2CRBSNLF22Ip5Lt/X/vMOYNXKFaxds4ZUKsX9/7yPkz55yk7XfLB8NUtXrGHpijWcdsZZ3PbLOzn51NNyugeVghv0iauq68/lhZNSct3VlzNt+gyuvPrarNc0NTaSSnlpu3v//AcOPvSwLtF7f+vvL0KIPs2Uf+Lf/+KEUz6V9bOjjj2JN197Gdu2icfbefetxUyeOr3H8gwBEYUd5sjCcJ+2jS0aMYqa914DYMM7r1A2bkKXa6SULLz9G4yonsL+p13cY3lSSiqLwn276T4wdY9RhENdvz9px5G2Nx4rXRunpQYRKQfAbdmAiJQjwtkjNhEqwm2vRbreqWtOSw0iWp71WttxqRpVlh9jsrDXtHFZBdUv+yaMU3skamVh7m0j320TYGxJ785Sfxi2Y+a59JWWZXHrbbdz6kkn4LgOF1x4MbNmz+Z3d/8agMsuv7Lf9Qu6puzziWmInGahXXzBZ1j0wnM01NczY8oe3PzNb5NOe5NbLv3clWypreXI+QfS0tKMYRj86o5f8Nqb71HSi+CpjOzAi+5629ThtVde4v5/3MvM2XP4+GHecb9f+9YtbNywAYALL72cFcs/5ItXXIJpGuw5fSY/u+PuXusWoHTZFkAsbNCaw8Eu8Xg7r7zwDN/4wS92/O5ff/09AJ8671ImT5vOoUcew9nHH4IwDE7/9AVMnd77sp6CsDr7CsMm0ZBBW6qroD/6k+upee914s3bufviIznk3C9y7NW38Mxvv4/rOFjhCMde9X8AtDZsYeEd3+SMb9/NpqVvsPSZRxg5YU/u+dJpAMw//zomzzuySx22KxmlUMz3nTk+6wliMt1Gev3THXlUiVk2FbN0IgDO9pVdUtDe9f8jPOVkjMIxGKVTSC37JwgDERuJWTE7a/1T9xjVa2ZgIOw/ewKxaIh0687t0y/7DtlH7fntU0YWsKqhvUv36UfbDJuCiSNiSuwSQZ9z3BfmzZsnFy9enPP17T7vW94ZU0DEUifmUkpak8HZFw2pTWUmbZemdsf37VwzVBRaStfRb2lKUtsUwKYjeM7KXuOLlDqbv3llg9KtK3tiTHGYrx7V+1DKQBh/9E3Ub29VWkc2hIBLzpjPHd84V1kdDY2tTDru66T9PtgCKIyFue2msznvlIOV1bGsro0/L95EMoBdt8Km4JrDJvQpOhdCLJFSzuvtumGbZgdQqKW91634mxVCBJpqVx2ZhxUKaW8YQr19pQUhpePyPVEcM5VvRXrA+FIiATzDkCGYN7609wsHyFnH7ad8bDcbBdEI55zQa78+ICrKipg9ZVzvFyrAdlxOOLznTZ0GyuQRMXWz0HohFjIZU6wmazS8xTygScsCf8a0wwF5K5YBhmIxEEIQDQdjX0HYUC520ZBBNOT/62cIGFWiLgWdYe8xRWpnEXaDBA7yQcyvOvcopanu7hhRVsBh+3fdjCTffPniYykqUDO22x2GIfjE4XMYWZ7/DVU6EzINDtqjDL99zZApOGpKubK+ZViLuRFQ9GoZasfLO9cTBH45EQXhYLyxmMLx5M6MKgn73j4tU/jyvVqmwSETSpVOJNwVA5gzpsgX+6ZOGMW+M8b76q8URMNcf8ExvvQtp35sLpbPzko0HOK6C4/xpa7DJpX5flCOlF7GShXDWswB/NYDgX8iK4Qg4vMURtNA6Z7lnbEMQSzk7wtXFDGUZx0ylMYsIj5G50LA+BFR3zqxY6ZV+DpcYpmCk2f1vgNevvjl188hqnCJ2K6Mqijm4jMO9aWuUMjk5zedTUFMfRYHIBIyOeaQGRy410Rf6qsoDHPQHqXKlontStgUHD99JLEsqyDyxbAXc0MIXyPYiOVPVJ4hZApf00V+i2tx1PQterUMtbO8d0UIwYSRMV+iOwGUF4QoUrRfeTZiIZPP7DvWlw4zbApOmzOKsph/4jpnWhXXX3QsBVH1gheNhLj3x5cqXV++K+d8Yh6H7jPZl7kB0WiYX32r950L88nJsyp9yeIIPOfhqCnZl+Lli2Ev5gAhw59NVvwYS94VP8eWYyHhe2pKCEFpzJ8XrjRm+W5fxDKoKosoF3TLFIwr93cMFGDm6CL2GVesVNAtw1vuc6APY+W7ctNlxzNlj8qs687zRUE0zHUXHMN+s7Jvf6oKIQS//b8LKC6MKm2fsUiIP3zvAirK1I6V70rINLhw3jilq3K8egQX7j9OuTbsFmIuhCBqqRV0U3hOQxAYQlCgWNAjFkqXavVE2DKUCroAygrMwOyrKA4zqjisrMO0DMHU0QW+DY/syqfmjmHqyAIlnaZlCMaWRLjkgCrfHTEAyzJ54u5rqB5TrkTQC6Jhzj5hf775+RPzXnYujBlZwlO/v1aZoMciIX761bM48Yi98l94Dkwoj3HefuqyR2FTcPnB1Uo3McqwW4g5qBV0S3hj80F0JhlMQ52gRyxPUIMkGjIoUyDoXvrZDNy+MWURRpfkV9CF8KKCaWMKArXPNAQXH1DFrFFFeR1DD5uCCeVRvnDIeOWb/PTEiNJCXrjnq+w5cXReU+4F0TCXnHEov/rWZwLtW2ZOGcvzf/kKFWVFRML5GaYRwhPyX379HC45c35eyuwvc8YUc+G8cYRNkTd9MARELYMrDxnPpBEFeSq1Z4b1pjHZkFKSdiFf+wVETP8mhOWClJJEWubFPiG81Ppgss92JU3tdl7sC5leCn8w2deasFnXkMBxB3ZYiBAwojDEuLIIxiCxT0rJGzXNPPDuFtKu7NfpceDZZhmCE2eM5LBJ5b4PbXVHOu3ww98+zs//8jTJVLrfzy8StiiIhvn99y7gE4rXXPeF7c3tfPF7/+DxF97NekRqrhREw4wbVcbffnIpe+1Zlcc7HBj1bSnuWbKJra2pXnef7ImwKZhYHuPcfcdSkoc5KrluGrPbiXkGV0pSTv+OowQvrR50NN4TtuOJen+fbtj0lqANRvuklLSnXNqSnqL3xUZvm11v1no0pH49eX9wXcnmpiQNrV6H2ZdXVAgvi1JdHvF1sltfaE7YPPDuFj7c6h21aef4EprCe98mjYhx1t6jGdmHPbb95N3lG7nqlr/x7opN2LaDneM+9bFoGCklnzp+P/7fl8+kvMSfiK6vPP7Ce1z/o39Rt72F9kQ66+mT2Sgq8OaGXHv+x/nqJccTUjjPoL+4UvLimkaeXF6P60qSfRD1iCkIWwYnz6xkv+qSvPUtWsxzxJWStJPTNuc7lp35tY58oEgpcSSkbEku/YkhPAEfSvYlbU/Y047ckSLb9VGKjt+FTEFhxPDSaUPAPseVbG9LU9+SImnLHaetdX5lM78TQEnMorIkHNj6/L7SnLB5ZV0jL69vojVpEzYNbFfidIi7aYBlGKQdl1jIZN74EuZPLGdEgX8zugfCsjW13PG3Z7l/4RvEEykiYYtEMo1tuwhDELZMQiGTeCJF1ahyPv/pIzj/1EMGrYh3RkrJS2+t5rY/P8X/Xv0QEJim4dnnOJiGQSRsYRgG8WSKWVPGcu0Fx3DGMfv4OiO/vziu5IMtrTyzahs1jQksQyDxHE9Xeu9dqON3roRJI2J8bMoIplUW5D1TpMW8H0jpPahdAwVDeP+GggB0h+w449mRnhhkTBR4neZwsM9xIe1IXCl3CJxhCEKGwBwiDkp3uFKSSLnE086OFLy3z4AgFjYJDREHpTviaYeNTUk2tyRJ2S4SCJsGo4vDVJdGKRwiDkp31NY38+bS9axYt5V4IoVhGJQURdlrzyrmTq/O6djRwYqUkrUbG3hj6XrWb95GMmkTCplUlBayz8zxzJoylnBocGaJcsFxJVtbU9Q0JWhJ2jiuN8xTGrWoLosysjCkdKhHi7lGo9FoNEMcfdCKRqPRaDS7CVrMNRqNRqMZ4mgx12g0Go1miKPFXKPRaDSaIY4Wc41Go9FohjhazDUajUajGeIEKuZCiD8IIbYKId4L8j40Go1GoxnKBB2Z/wk4IeB70Gg0Go1mSBOomEspnwe2BXkPGo1Go9EMdYKOzHtFCHG5EGKxEGJxXV1d0Lej0Wg0Gs2gY9CLuZTybinlPCnlvMrKyqBvR6PRaDSaQcegF3ONRqPRaDQ9o8Vco9FoNJohTtBL0/4OvAxMF0LUCCEuDfJ+NBqNRqMZigR6yKyU8twg69doNBqNZjig0+wajUaj0QxxtJhrNBqNRjPE0WKu0Wg0Gs0QR4u5RqPRaDRDHC3mGo1Go9EMcQKdzT5YkFIiAVeC43r/n0EIz+MxDDAEGEIEdZv9Rkq5wzbHlTvbB5gGmIbANEAMQftcKbEdSbrjn9vxPAXe8wpbAssQhEwxJO2zHUk85RBPu7QnHeyOByiEIGIZFEQMYiGDaMgYkva1Jm3q21LUtSWpb02RciUgsQyDioIQo4oiVBSGKYlYQ84+KSVrt8VZWdfG+7UtrGloJ2W7CCEoDJvMGlPEnqOKmD6qiJFF4aBvt8+kHZdV9e2sqGvlvU0tbGpKkHYkpikoi4WYM7aYPSsL2XN0EcWRoSc3acelvi1FfVuKra1J2lMOrpQYhqAobDG6OEJFQZiKwhCWEWxsPPS+3TwipcR2Ie32dA24AI73syEklgGmGPzC50pJypaknZ6vcxzA8QTCNCQRS2AMAfvSjqQ95ZC0JQKQWa+SpJyPPo+GBAVhE8sY3LZJKWlLutS3pmhPugixs5PZcRWJtEtzwnNchICKwhDlhSEsc3Db57iSddvbea+2heaEjSEg3dVA6ttSrGpoRwIFIZM5Y4qZNKKAkDm4k4qtSZv/fljHA29vpi3pvYAJu2tH8/7mFmIhg7QrmTqykE/tO5YDJ5RjDvL2ubUlyYL3tvDY+1sBsF2XlNP1+b1d00QkZJJ2XA7Yo4wz9xnLzNFFg75v2dae4v3aFtZtb8cwRJcgyLsmTU1THFMIJDClooCZo4spjYYCuWchZfYucDAyb948uXjx4gGXk4uI94YAwiaD8qWTUpJMywHZZwiIhQTGILTPcSXNCYd0ls4jV8KmoCRmDspMSzzlsGFbMmsH0hsZayqKLCpLwoPSvvXb47y4dpuXUemjgRkn7MDxZUwdWTjoRMFxJQ+9vZl7Xt+IEJDMIuA9EQsZxEImNx47lb3HlSi6y/4TTzvc/eI6nl5eDzK7A9YdAohYBlVlUb527DSqyqLqbrSftKVsFq3ZRl1bCteV3QQI2fEygVBVGuPQieVELDMv9ySEWCKlnNfrdbubmLtSkrS7i+L6jik8UR8snYrtShKpvjXCngibELYGR3paSi8SbU26ebFPACUxk4g1OKI8V0rqmlM0tNoDtk8AlikYPyJCLJyfTmWgJG2Hl9ZuZ2NzAqevXsouWIagojDMEZMqKBgk9m1sTPCDhSvY1JTIGoX3hbApOHrPkVwxfwLR0OCw751Nzfxo4UraUnbWKDxXDAEh0+D8A6o4fe7YQeFwSilZ1dDGq+sbcWXfnejOGMJrn/MnVbBHWWzA96bFPAuOK0n2knLuDwKIWsELesp2Sdr5L9cyvPR0kPZJKWlNOsTT+W+vhWGDwkiwHabrStY1JIin8uOoZBACxo+IUBwNdkStLWXz+IdbiaedAXWUnfGyYwafmDGK0lgwqc0MS2tb+Pp/PiRpu3mzL2wKxpVG+fFpswIfb35qWR13PL+2z5mGnohaBgdMKOPGY6YGmuGUUvLahkZW1LcN2MnsjGkI5o4tYa+xA8uw5CrmgyMk8QFVQg5elJ+wvUYRFKqEHMB2IZ6WgdknpaRFkZADtKXcHeOaQeBKyVoFQg7enI8N25K0JhQ1jhyIpx0eXbqlY/JQ/sqVQNJxeezDLTQn0vkruI8s39rKzQs+JJ7On5ADpBxJTWOCLz/4Pm2p4J7f08vquD3PQg7eHILX1jXyw4UrcAPsW15dn38hB09z3t7czLubm/NabnfsFmIupToh31EHkHSCEXTHlcqE/KM6IGUH88LF0y4JRUKeoS3l5r2zypVN21MkFAh5Bilh/bZkIPa5UrJwWR0JW519KUfyxLKtpB3/7WtsT3Pzgg8HnFbvDtuV1DYnueWJFYH0LR9uaeX259eSUmRf0nZZvKGJe16rUVJ+byyra2VVQ/6FPENG0Dc0xpWU35ndQsz9Crpc6UWxfiKlJJ7y5yVPOShr9N3huJLWpD9fanPc8T1CaE3YNCcGPkbeG1JCzbaE74Lw7uZmWlM2qqtN2ZLFGxrVVpKFW59ZRdJW28GkXcmyLa08vaxeaT27krJdfrBwhXInMGm7PPh2Lavq25TWsyutSZslNU19noTZVxxXsmhNg/J2MuzF3O7HjOCBkHbxVRCSdv4mu+WCn+l2KSVNcf/S3xJoSfhXn+NKarYnlQtdhqQtaWj1Lx3dGE/z7uYW5Z0lgCMlqxraqW1JKK8rwwurGnhvU4svDnzCdrlz0Voa2lLqK+vgj69soDnuT3tJOS7ff3IFtk/ZFSklz61u8C04sV3Jy+u2K61jWIu5lL2vsVaBYgdsB0HYJ33MPqSdvi9dGihJ2786t7elfXU0pYS6lrRvzuabG5twfHRsHSlZUtPkS11SSn770npl6fVs2I7kgbc2+1JXcyLNYx9sITmAWet9pTGe5uW1agUvw9bWFI3xtG+BkCuhpjGudG7HsBZzV+ZvCVpfsKU/Y+dBjWH7VW97Kpgx7HhKvYckpfSWoAXwCFt8yHYk0g4bm9SPE+7K9vY0TT5Ek+9saqZF9USVXbBdyRNLtyobv+7Mk0vrPtq0wCfiaZf73tjkS13v1Tb7Hii4EpZubVVW/rAW84FsmjJQVDu0QWUdILPtrVoDXVcOaC3rQPBjKKEt6QYyg9eVUNeqPlW7os7f8c8MrpQs3dqivJ4H3tpMMqAOZtHqbUrLl1Ly0Nu1gQQLG7YnWLetXWkd8bTDpmb/hmMySGClglnzGYbtdq6Z/ch7Y+GTT3DD9dfiuA4XXnwpX7nhpqzXLVn8Okcddgh/ufcfnH7mWb2Wa7ve+mxVyByzDl+44lKeePxRKitH8eqSd7p8vnzZh3z+8kt5+603+NZ3vsc11305p/od19vTXRWdt2Dtjo01G7j6ikuo21KLYRicd9FlXP6FL+50TXNTE1/43IVsrNmAY9t8/prrOfe8C3usW+BFQSGFW6K2Ju2c2mdzUyPf/upVrFz2AQjBLbfexT77H7Tj8z/cdRuPPnQfAI5js3rFMl54ey2l5SO6LTOZlriuVLq737rG9qwO7e//7yu8tehpSsor+P59TwHwwF0/5c3nFyKEQcmICi779q2UV47p8rcL//57nnv470gpOfK0czn+M5d1uUYCNU1qO2opJe9saunSNuPbtvDun79LsrkBYRhUzz+NiUefw4cP3k7du4sQpkVBZTV7nf8NQgXFXcpd+7/7qHnxEUBSPf9UJh796S7XxNMuL6/ZxtF7jlRjHLClJUl7lkjBD/tA8tbGZiaMKMi/YR1saUliCtHFmVbdNsHrWxraU4wqiuTdrmEbmecidI7jcP2XruahBY+x5O33+dd9/2DpBx9kve4bN9/EMccdn3v9ip3aXIPWz55/IQ8+8li3n5eXj+DHt97GNdfmJuI76lccmaed3pcyWZbFd7//YxYtfpfHnl7EH397F8s+3Pn5/eG3dzF9xkyeeWkJDz72FN+5+QZSqZ4jUwkD2io2F9pynKH/o2/fwPyjjmXBc2/y4MJXmDx1+k6fX/L5a3lg4cs8sPBlrr3pu8w7+LAehRy8HaoSCqNKKSWN3axrP+yTn+LLv/zLTr878fwr+N7fF3LL355gn8M+ziO/+0WXv6tZuYznHv473/rzAm7525O8vehpatevyVpHe9pROpFqS0sy6++FaTL9zGs4/Nv3cfBXf8f65++ndfMaRs44kPnfuJfDvnEvhaPGs/rJP3f525ZNq6h58REOufEPHHrzPdS9u4i2reuz1rNsq9qsx4q6NswsG0T5YV/Kkby3Se267Pq2VNZtaP1om66UyiYxDlsxz0VrFr/+GpOnTGXS5MmEw2HOOvsc/rPgkS7X3XXn7Zx2+hlUVo7KuX6J2nHzXMV0/mFHUD6i+869ctQo9p93AFaobztoqZ50mouYjh4zlr332ReAouJipk2fQe2mncfchBC0trR6B5e0tlJWPgLL6j0hpVLMpZQ5LfdpbWlmyasvcua5XiYhFA5TUlrW7fWPPfwvTjz1UznU70V4qmhNOd1609P3O4jCkrKdfhcr+iiKS8bbs+40uGntCqbstR+RaAzTspi+38G88ewTWeuwhGCbwnHz7sQuWjqS0j1mePcQLaRozEQSjVsZOesgDNNrc2WT5pBo3Nrlb9tq11I2aTZmOIphWpRP248tbz2Xtf5tbWkSCsfYlm1pJZ5lFq9f9ql2VrZ0s+LBj7bpSNjSmt0ZHCjDVsxz0dFNGzdSXV294+eqqmo2b9rY5ZoFjzzMZZdf2fd76PNf9KHsgHfhVV19XwP/9evW8t47b7PfvAN3+v2ll3+B5cs/ZO89J3DUIfvxvf93K0YORxWqzDx4jl7v19WsX0v5iJF84/orOev4Q/nWV66ivT17RxePt7Po2ac49sRTc6pf5QYr8ZRDXzP49//qx1x/0kG8/MTDnH5F1yxR9ZTpLHvzVVobt5NMxHnnpWdo2NL9zO52hZMYt7Wle/3+2hs20bxhOWUT5+z0+5qXFlA565Au1xeNncy2lW+Ram3CSSWoe/8lEtu3ZC07ZAqa4uom321u7n25pEr7mhXvVtjeR0c2322zVdHEyWEr5rmQLXLe1fO64cvXccsPfoRpDo7DDjIMnR311dPW2sql55/DLT/6KcUlO++D/MzTC5mz11zeWb6O/y16na999VpamntP46l2xHLROtu2WfreW5xz/mXc/+RLxAoK+P2dt2a99tn/Psa+Bxzca4o9g8pREidXAztx1hdu4GePvsohJ5zG0//8U5fPx02axokXfJ6fXP1Zbr3mfMZPm9ntOylRu9dD2ul58qKdaOetu7/GjLOuxYoV7vj9qsf/iDAtxh54Qpe/KRo7icnHns/i27/I4juupaRqGsLMnkESQq0z1lvZqu1TPrm2j20jn23Tq7+vd5wbu7WYV1VXU1Pz0TaCGzfWMGbsuJ2ueeONxVx43rnMnDaJhx+8n2uvuYoFjzzs8512JfhzhgYH6XSaS847hzPPPpeTTjm9y+f/+OtfOOmU0xBCMGnKVPaYMJEVy5f1Wq7K71eI3JyFMWOrGD22ir33OwCA4046jQ/efTvrtY8/cn9OKfYMKs+1MHI1MAsHn3Aai//3eNbPjjz103z3r49x8933U1RSxujxk7JeJzL3oAjLNLqdPOg6Nm/+9muMPfB4xuz7sR2/3/jKo2x970XmXvzdbg8sqp5/Cod+7S8cdP2vCRWWUFhZnfU6KVF6Xn1PZ8X7YZ/qQ1f62zby0Ta9+vtVfa8MWzHP5XntP+8AVq1cwdo1a0ilUtz/z/s46ZOn7HTNB8tXs3TFGpauWMNpZ5zFbb+8k5NPPS23e+jHfedK0KcGqq4+lwYvpeS6qy5n2vQZXHn1tVmvqRo/nhee/R8AW7duYdWK5UyY1P2L9lH96iwU5Pb9jRw1mjHjqlizajkAryx6linTZnS5rqW5icWvvMjHjj8p5/othUsRYiGjT1reebLQm8//l7ETp2S9rnmbt51pQ+1GFj/zBAcff0rW6xAQDamzr7wgtONc9c5IKXnvnu9TNGYikz7+mR2/r3v/ZVYvvIf9r/wJZrj7M7yTLd6Ss/i2Wra89SxjDzgu63VpV1ISVXdKXGVROGv79Mu+IsUnGPalbeS9bYKyI3uH7dK0XMTAsixuve12Tj3pBBzX4YILL2bW7Nn87u5fA/RrnLwzKo8MNQ2R0yStiy/4DIteeI6G+npmTNmDm7/5bdJpb3LQpZ+7ki21tRw5/0BaWpoxDINf3fELXnvzPUp2SVfvSg7DzgMiZIpeN3V47ZWX+Nc/7mXm7DkcPd87IfDmb93CxpoNAFx46eVcf8PNXHPlZRx58L5IKfnmd79PRUXvy3rClkIxF4JIyMhpRvnNt9zKjV+8lHQqxfgJk7jl1ru4757fAXDO+d7yl6efWMChRx5NQUFhT0V1qt8TXFUUR6xuJ3/e9fWr+XDJy7Q2bue6kw7ktMuv550Xn6F23SqEYVAxpoqLvvZDALbX1fLH793I9b/wZkffceMVtDZtx7RCXHDDLV0mK2WwXUlFQViJbQBTRxZmHVNuXPU2m157nKJxU3jxB+cDsOcpn2fpv36Gm07x+u3XAFA2cQ6zP3MjicY63rv3B8y76ucAvHX310i1NWGYFrPO+QqhguzvYFnMUnqG+4zRxSz8sK7LJEm/7JtWmVs77i+jiyJsa+86QdKPtmkIGF3UvcMzEIbteeZSShTOEekVQ0BUoSC4rqTNpwNWshGxhFLBS6Rdmn3cJ70zAiiNmYQVbhSwuTHJtrbgGuiMsQVK05kL3q9VOqO8JwpCJp+aO673C/uJKyWn//b1wDY1OmRiOd/6xJ7Kyt/YmODqf73r61a1GUKG4OKDx3P63LHK6lizrZ2X1m7zfQc48IKUo6eOZExx7oK+259nLoRQOi7YGyo3jIHg0+wqN4wBlG7Y0hvSh/qLomZg7TNsCeXjktVlsUDsE8DYkvxvyNEZQwhmjum6KYofRC2DgyaWKa1jbGkksPfPMAR7V/WcFRwoo4rCgZ2f7rooyxoNWzEHUJhJ7BXV74IQAoWZtl7qVj9JxTREYB1KNCSUDpGANy4YhENmCBhZpG68NcP0yqJAJmkaQjBrtHqhPXPuWKVDFd0hgSOnViitwxCCU/YaQziA9290cZgpI9Wm2QvDFqMV7MDWGwKYXBHrcYLhQBjWYh5U5GMKtePlGUIK09w9EfGp3sKwEYggFITUe0lCCCoKQ77bJ4HSmPqpMgVhs0+pxHxRErUYoXC8PMP+e5QSsfz1pk0BR+85kqgP7fPE2aN8X/4atQw+ta+64ZHOzBlTknUSo0oMIZg5Wl3WYViLuRAikOjch3cN8BqH3/YJ1A8hZAiZQvlEu10Jm0Lpsp/OlBeFfI3ORUdUrnJP9s7sW1WqPEPVGdMQ7F9d6ktdhhBcdFA1Ub9eBrwVCOf4JHYjCsJ8bFqFr9F5YcTkiClqsw4ZxpZEKIr4N/9b4GUdymPqsmLDWszBEx4/HbCQoXZZ065EQsLX6C4aVp+CziCEoDTq7wtXHPUv2rIMQVV5xDdBDxmCymL1KfYMFYVhZowq8kXQDQHjS6NUlcbUV9bBcTMqmTKywBf7opbBRQeNZ7Ti+QCduWL+RGI+jeWFLYOvHTtN6aTTzgghOGpyRdZteVVgGoLDJikeHlFa+iDAz7FlQ/gXtWYQQhAN+dMgQya+p6YsU1AY9udLLYoYyucC7EpJzPLGzxXXI4DxFVHfHLEM+1aVEQupty9kGBwyIbfd7/KFEIIbj51KSPFLbxmwR3mMU/YarbSeXSkIm9x4zFQiiu2LWAbHTx/J7LH+TiosjYWYO67El/k/B08oJ6Y4ZTvsxRy8SFnxPgQIIGL6M1a+K5YpCCsOYE3h31j5rhSEDeV1x0LCtyhkV6rKI4rXtUNVeVjpRirdYRqC46aPUjbpBzwH87jplb5FdZ2pLIrwnU/sqUzwTAHlBWG+e9J0XzN+GfatLuXCA6uV2Rc2BdNHFXL5/AlKyu+NOWOKGV8aVSboliGYUVnIlAq1k/pgNxFz8DoVVYIugIgVjJBniFiGsgyEKSDmY3p9V4QQlERNZYIeCwnlu071hGkIJlXGvFn0eS5bCBhXGqa0wL/0+q4URyxOnDmKiJX/CY0hQ3D89EpfJr11x9yqUr55wrS8C55lCEYWRfj5GbMpUzjW2hunzx3L+QdU5d2+iGUwa0wx/3fSDKU7EvaEEILDJ1coEXTTEOxZWcj+1WV5Lbc7hu2mMd3hSknSzt9BGoYILiLPRtqRJNL5e6Yh04vIB4N9UkraUy5tqfxsZiHwUutBReS74rqSzU0pmuL2gE/FE3i79I0fEaUwQEelM20pm+dXN7CtPT3gDTssQ1AcsThqSoXSrU37wsq6Nr735Aoa4+mcjrjtiYhlsG91Cdd/bArFPs4b6YmX12zj1v+tJmm7A35+YdPglL1Gc+GB1YEJeWeklLxX28Lbm5q9g4IGgBBgCsGB48uYVlk04HvLddOY3U7MwXtwaQfsAZoeDmAMORdc6Qn6QA5WEnjRuN9jyLlgu5KmuI3r9t8pC5letD8Y7WtNOtRsSyBl/05YEkBpgcWY0vCgs09KyYr6Nl7f0Iik7ydkeZPNBHPHlTB7THEgqeeeSDku97xWw7/frUVKbx/1vhC1vHkb131sMvMn+zsHIBea4mlue3Y1b25oJu26fW6fsZBBaTTEzcdPU75ta39oSqR5blUDLUm7Xw6LZQgqCkIcPrmCwjyNfWoxzwFXSuw+irrAi1b9Wks+EBxXkrIlfQkSjI6xcdMY3PZ5DpmkLeXmtEe9wBP+sOVNqLOMwZFt6A5XSlriDnWtKVK2RNCzsBvCs6+swKKiKKR80tJASdkuKxvaeL+2hZTjIqWku8doCG/ei2kIZo0qYs/KIl/WWg+EutYk/3lvK/95fwtIScqR3YpD1PIOpqkoDHH2vuM4cmrFoLdvdX0bD71dy/OrGrAMQcLOLuwC72ATx5VMrSzkU/uO44A9ygadk9kZKSW1LUneq21mS0sSQ4genbKQIXClpKo0xpwxxYwsDOe1b9Fi3geklLgSnI5ISEqvY8w8DkN425caouPEq0EsAtmQ0ovSHdf7/87tUnTYZhoCU+DbGuR84riesKcdT9g722d0rGYImQYhSwy6SC4XkmmX9pRDe8olnnJ2tFEhMs6JSSxsUBA2h9zzk1LS0J6ivi3FlpYU29pTO0TPNATlsRCjiyOMLAgzsig85J5f2nH5oLaFFVvbeHdzC+u3x0k7LgJBLGywZ2URs8cWM2N0ERNHxIZc39KWsnl/cwvLt7bx7qZmaluS2I7ENAQlUYtZYzzbZo8p9nVZXb5oTzlsbU1S15ZkS0uyw2mRGEJQEDIZXRyhsjDMqOIIUUWbCGkx12g0Go1miLPbH7Si0Wg0Gs3ughZzjUaj0WiGOFrMNRqNRqMZ4mgx12g0Go1miKPFXKPRaDSaIU6vq9qFt1biQKAKb8XWJuA1OZSmwWs0Go1GM4zpUcyFEMcBvwJWABs7fl0NTBVCfEFKuXAglQshTgB+AZjA76SUPxpIeRqNRqPR7I70Fpn/AjhGSrm28y+FEJOAx4CZ/a1YCGECdwLHAjXA60KIf0spP+hvmRqNRqPR7I70NmZu4QntrmwEBnq6wYHASinlaillCvgHcOoAy9RoNBqNZrejt8j8D3gR8z+ADR2/Gw98Gvj9AOuu6lQmeE7DQQMsU6PRaDSa3Y4exVxK+UMhxMN4EfMheFuT1wCfzUM6PNsmxF0m1QkhLgcuB9hjjz0GWKVGo9FoNMOPXmezSymXAkuFEJUdP9flqe4avCg/QzXeTPld678buBu8vdnzVLdGo9FoNMOGHsfMhcd3hBB1wIfAMiFEnRDiW3mo+3VgmhBikhAijJe6/3ceytVoNBqNZreitwlw1wLzgQOllBVSyhF449rzhRDXDaRiKaUNXA08CSwF/imlfH8gZWo0Go1GszvSW5r9AuBYKWV95hdSytVCiPOAhcDPB1K5lPIxvCVuGo1Go9Fo+klvkXmos5Bn6Bg3H+jSNI1Go9FoNHmgNzFP9fMzjUaj0Wg0PtFbmn2uEKI5y+8FEFVwPxqNRqPRaPpIb+vMTb9uRKPRaDQaTf/QR6BqNBqNRjPE6XXTmN0JKSWu9LahkxJExx51hvDGFYTItmnd0EFKiSM92zK77wjANIaHfY4rsV3vGYJnkyHAMgXGELdNSs+2pC1xXYnEa58hQxAJGcPCvnjaJZ5ycTseoGFANGQSCw99+xxX0pKwaUs5uFJ2vHeC4qhFYdgc8u9e2nFpTtgk0t7zEwLClkFx1CJqGUPevpTt0p5ysB2vfzEEhEyDgohByBwcMfFuL+ZSStIuOG6WvWR3QSAJmWCKoSN8jitJOxI7B/tMIQlbYsjYl3l28ZRL2ukQOLLbKYCwKYiFDULm4LcNPPvaUi5N7TZJ281+kfCcM9OAoohJacwaNJ1LbziupKElRV1LmkTa3eE8d3Y0wbMvbBlUFocYWRzGGiLPL5l22LA9waamBIm0i2kIpJQf2dfRWCVQFLGYMCLG6JIIpjE07GtN2qxriLO1JUnakTvZlwkOXOkJe1ksxISKAkYWhoZM39KadKhtStGasHHdj4K7na/z3r2SmMXo0giFkeBGpoWUQ2eH1Hnz5snFixfnpSxXSlIOO6K4vmIJCJmDV/QcV5KwZb/s84QPQqYYlPZJKUnaLm2p/tlnCiiKGIStwSl6Ukoa220a4/ZOWZRcEEAkZDCyKERkkNpnO5KN2xPUt6YR5P4OZppieYHF+IrooHVa2lMOH9a20tDmLfjJ1T7TU0DGl0eZWlk4aEV9e3uaD2tbaU32rX2awstGTK0spLo8Omj7lm1taTZuT+6IwnPF6MhGjB8RpbQgf3GyEGKJlHJer9ftbmLupSsh3U2g0xcyojeYXjpP6GRe7DMExEKDK0XtupLmpEPaGXhZEVNQFB1cKdyk7bKlOYXtyD6J+K4IoKzAorzAGlSdZlN7mtV1iR1DBf3FEDChIsqIosET6UkpWb8tzvKtbX12wjqTSeHOrS6hvGDwbOfhuJLlW1qpaUz0OwiCjDNtsXd1CQXhwTPHOm27rKmL05p0BmSfIbx3b4+RMaw8aIMW8yxIKUkOIBrvDsvwUrhB40pJe2pgnWQ2YpYYFKnNtCNpijt5tc8TPTMvL91AaY6nqW+182afwMuujCsLPnXrCV2ChpZ03t4/Q0BpzGLSqFjgDpntSt5Y30RzPI2TR/smjyxgSmVhfgocAPG0w2trG0nZbl6eX0cSgn2qS6gsjgy8wAHSmrBZXtuet7aZmRMxfWwBsQE6LLmK+eDMUylASknCzr+QA9gupOxgnSJXStoUCDlA3JaknTyE+gMg5Uga8yzk4EVPje3exJYgaWrPr5CDZ1vKkdRsT+KoaPi53oeUrKmL51XIwXuXm+I2K2vbcQMMSmxX8traRpryKOTg2be6vp3lW1rzV2g/iKccXl693Zvclif7JJ59b9U0U9uUyE+h/aQlnl8hB88+25Us3dRGeyoPacQc2C3EPBORq3zdbelFjkEgOyJylSRsr3EGgd0RkatCAo1xJzDBa0nYNLTlV8g7Y7uSjY3JwASvZluSxnZbiSPtSmhNOqzZGs9/4TkgpeSNdY20JdXZt25bnHUN7fkvPAfSjsuraxuV9W2uhHc3tbCtLZgNRdtTDiu25FfIO+NKWLa5rfsJrHlktxBz21UTke9K2iWQDjNhq4nIdyWelvg9LCOlpCmh3rOVQHPC8d0+25HUtaaVPz/b8Sb2+E1LwmZrS0rp+5eJ0IOwb922OE0JNUKewZWwfGsbrUlbXSXd8P6mFlKKhSgTods+Z/9cKVmlUMgzOC6s3tquvG8Z9mLuyvxMBsuVpI2vgmB3LDvzi4TPwwntqfyl9nrDdiHhY2ORUrKlJYUfzUUCzXHHV/scV7J6a9wX+1wJ6+rivg4HtaUcVmxt86V9ZgTPz2Bha0uSutaUL4GC40qW1vo7nLB5e9K3bGo85VLXojb7MOzF3Kfhih1IyOu4WY91Sem7uNouvqWjHVfSnvbXvtaUf9mH9pRL0k/nAdiquEPpTG1T0tehGVfCxm1J3+r7YHOLb44mQCLlsKnRn/FlKSXvb/LPPldCbXOS5oQ/2YeU7VLbrDZj1BlXesNNKvvOYS3mruzfOuSBknb8ic5z2ehGBSm/vFk/Uyqd8Ct6bYyrGyfvDtuRvtgnpWRrsz9Zhx11Ag1taV+czXjaobHd37S+I2FNg/p0LUBda8r3OSSuxLe5AVubU4F0ng2t6pzpYS3mts9ReYbMTE3VJAOacGe76p0VKSUJn6PyDO0+zA1IO/5G5Rm8yX7qo5/GdttXIc8ggHofsg/rt/kzfLArybRLkw/Pb019u28Zxs7UNieVD5W4UlLX4s/wwc71Qm1TSlnfMrzFPIfvbOGTT7DP7BnsNXMaP/3xj7q9bsni1ymOWjz0wP051a166E7mmHX4wuWXMnn8GA7ab+9uy/nq9V9i7qw9OWTePrz15hs51a96nD6XsayNNRs4/cRjmL//Xhx+wFzu/tUvu1zz4gvPMaWqgo8duj8fO3R/fvqj7/VarivVD5W0JnJbZnfk/jM48cgDOPljB3HasfO7fL5qxTLO+sRRzKwu43d33pZT3e1J9RP96nNchva339/J2ccfxDknHMzXr7mEZDJ7Gvn9t5dw0NRynn7s4R7LcyXUt6qPmDc1JbM+v59940ucc/gsrjj1iB2/e/7Jf3P5KYfziTmjWf7eW92Wme1vd8WRsFnxUq60073DoNo+IbysgEpaE05Ojli+2yZ4mTFVTvywFfNcOivHcbj+S1fz0ILHWPL2+/zrvn+w9IMPsl73jZtv4pjjjs+5ftVikGv5nz3/Qh7892Pdfr7wycdZtXIFb72/jF/c+Wuuu+aq3OpXnHpI57ADmmVZfPcHP+bFJe/y+P8W8Ye7f82yD7s+v4MPOYxnXlrCMy8t4Ss3fSOn+lWvO+/LEMJfH3ycBc+8ysP/fbHLZ2Vl5XzrBz/lsi98qU/1q57405bsPS22tXYT9/351/zlkWe574lXcF2HhQse6HKd4zjc8eNvc/DhH8+p7kTKVeqspB2XdDfe7LGnfZrv/eYfO/1u4tQZfPMXf2TOvEN6LDfb32Zju+L0fnPC7naTIdX2OS7Khy/akr2Luaq2CdCW0mLeJ3LRmsWvv8bkKVOZNHky4XCYs84+h/8seKTLdXfdeTunnX4GlZWjcq5fojYVnauYzj/8CMrLR3T7+WML/s25nz0fIQQHHnQwTY2N1G7e3Hv9ip2VXLRu9Jix7L3PfgAUFRez5/QZbN60KT/1KzYwX8t9KipHsfe+87Csvm37qXLdq+3InNun7TgkE3Fs2yYRj1M5ekyXa+7782/42PGnUj6yMqcyhUDpEEZPYrfXvEMoLi3b6Xd7TNmT8ZOm9lputr/NRltKbWalKd79vAM/7FMt5q2J3OaqqGibrvTqV8FuLeabNm6kurp6x89VVdVs3rSxyzULHnmYyy6/ss/3oFIO8hUYb9q0kerq8Tt+rqqqZtMu34HK+rujr5H/+nVrefedt9h/3oFdPlv82iscdch+fPqMT/Lh0vdzKk/lLOzMUbS5IITgorNP5tRjDuUff/l9fuonf85ENhK2m/WEqV0ZNWYc5132RU4+bA6fOHhPCotLukQ4W2s38ezC/3DmZy/p0z2onDzZlnQC3XEO1DpjLTmKnSraFc8lyaVtqGyb7Toyzz/ZvNtdD2244cvXccsPfoRpDp4DAYC8Tb7J5TsIgr6Y19rayiXnnc0tP7qV4pKSnT7be+6+LPlgFc++/AaXXXEVF557Vt7r7yt9eXb3/edp/v30y/zh7w/z1z/czWsvL8rLPSjd5MSV5NKCmpu28/xTj/LIc+/w+MvLSMTbeezh+3a65me33MQXb/xun98/tZu4yEAmv2UQQigd5gpy619Qv/FWLuYpbZuKvt/d+jzzqupqampqdvy8cWMNY8aO2+maN95YzIXnnQtAQ309Tz7xGJZlcfKpp/l5q13p7uDuPlJVVU1NzYYdP2/cWMPYXb6DwUw6neaS887mzLPP5ZOnnt7l887ifszxn+DG679IQ309FSNH+nmbO9MHX2n0GO9ZVFSO4tgTT+adNxZz4CGHDfwWFPpruRb92ovPMq56AuUV3rP42PEn886SVznxtHN2XLP03Tf5+jVe5NO4vYGXnl2IaVkcddwn83IPQxOp9gEG/O2prj2X8pW2TUUGDlsxz+UL23/eAaxauYK1a9YwrqqK+/95H3/8y707XfPB8tU7/vvySy/mEyeelLOQD4XX7ROfPJm777qTs87+NK+/9iolpaWMGTs2T6X3n1x8FSkl1171OfacPoPPf/G6rNds2VLLqFGjEULwxuLXcF2XERUVvdavMmWV67Nrb2vDlS5FRcW0t7Wx6NmnuforX8vLPag8Zcw0RE5+5phx43n3rcUk4u1EojFef+k5Zu61707XPPL8uzv++ztf/TyHf+z4nITcUPgALUMgRP6yY31FSpSe8hcK+IRE1Sfg5fLVqWqbAIaiZzd8xTyHayzL4tbbbufUk07AcR0uuPBiZs2eze/u/jVAv8bJ/cI0RE7juhef/xkWvfAcDfX1zJiyBzd/49ukbW+CyaWfu5LjTziRhU88ztxZe1JQUMCv7s5tXFb1iZqWIXrdnObVl1/kX3+/l5mz5/CxQ/cH4Ovf/h41NesBuOjSK/jPww/wp9/djWmZxKIxfvPHv+Y0jKDyyFchBFYOz6++bitfuOjTANiOzSlnnM2RRx/H3/70WwA+c9HnqNtSy2nHHUZrSwuGYfDHu+/giUVvUFxc0m25hoCIpU7toiEjp1TmnH3m8fETTuW8k4/AtCymz9qb0z99EQ/c67XBMz97ab/qd0HpOdlFEQtDiKzp4B9+5Qreef1Fmhu3cd7RcznvqhsoLi3jrh/cTNO2Br71hc8wefocfvDbf9KwtZbbvnUdt/z6793+7QlnfrZLHUIIpUcul8ZC1DYnsz5DP+wriqgd0iyImCTtniehqWqboM6+YXueuZQSH/ZW6BZDQNRS98IFsdVpZ0KG12mrIp5yaVU0UaQ3BFAcNZQKXm1Tija/9xruQAB7jIhgmerse2dDS2DHAhsC9p1QrGzuh+NKnv6wPrBJYqUxi4MnlSsrv7E9zZL1TYGdkjipIsaeo4uUlb+1OcWGbYlAMiuGgEmVMcoLc199stufZy6ECHTkR/UXqzoy7g2VaT5QGxn3hkS9fbGwEVj7FIJul1blC9XRVU8UhE2lkzhNQyh1ZHtCABV9EIL+UBy1cAIK8kwDygrU2lcYMQN792RH/SoYtmIOoDCw6r1uxX2Zl6pVW0dPKAzqAO/ZBeWwWMbwFrviqFqxA6gsCQfy/AwBo0rDyuvZY0QsEPuEgKqymNI6TEMwqkj9d5gNgWCk4roLwoby97vbukMGYUUdtxZzBRhC/SQOQOm4WU+EDPXL14QQxEIBvXBh9Q3HNITScd3uEHhpWtUURcxAsitCQFmBevuqyqLK68hGSdTypd1MrCjA78cngOryqPK+UwjBmFL/nU1DwJiyiLrylZU8CBBC+N4gwRM7P/CcBn/q6kxY4VyAzkRD/qeihfDPSSovsHy3LxoyCKlOq+C9e2NLI762TyFgVHHYF0c6ZBqMLfHXPlPA5MpCX+oqjVlEQ/46m0J4GQ8/GFnsf+bBEEKpozmsxRzA7+DH8GE8MoMQQukku2yETX+yDuDVUxj2176SiOHbpjmRkEFR1L/xOwFUFqsdj+zMyOKQr2PLIVMwtlxd5LMr08cU+fYuCGBEYZhKn9LfQgj2ri7xzVkxBUypLCDmkwNhGoI9KqK+2ScETB4VU9q3DHsxF0Lg5/Ck30OhpiF8c1gMH6PWDF4k6U9dEUsoG8/qjpFFIaVrojMIoKLI8iUq31GnEEweVeBLhykETBlV4Ju4ghed711V7It9piGYM65YfUWdKIlaTKzw5/nFwiaTKgrUV9SJiqKQsslonfEcMYsSxcNbw17MwXsR/AhgI2YwW6GGTeHPCxcSvtsnhKAkYqrd8AovMiiK+P86GEIwpiSifIOhWNigJOr/thLRkMEeFVGlz88QMK4s4kvHvCuVxRGqytRGeIaAudUlvjua4EXLxVFLqX2WIdhnfGkgfcvkUTHlczvClmCPCvXDB7uFmAOETJSOn4d9mAHdHUIICkJqBd0rPxj7DENQHlMn6IaAsgIzMPuiIYOxpWElgi46yh9TEg5sz/2RxWGqyyNKnp8hYFRJmLEKJxb1xswxRYwpiSjpXwwBe40rVj7Du/v6BfMmlHVslJP/8i1DcMDEMgoDmAwKXnZl5thCJUtRBZ6QzxhX6Is27DZi7u2ahJIIPWIGuy4aPhL0fGdRBXSUG6x9Zoeg5/trtgwoD1DIM8TCJuPK8juhSuDNKh9bGpyQZxhdGmFinscohYBx5RGqRwQzs/yj+/BS4OPzuFxN4LX5fcaXMqY0WPssQ3DgxDJGFIbz9v5ldiE8aFJZIBmjzkRCBrOqComGjLw9P0NAQcRg1rgi34a2hu0OcD3huJKUM/BzSgwRXGq9O6SUpB1JMg+bi4UMbxx5sNnXnnLzsvtdYVgQC/k34S0XHFdS35qmLen0u30KPKEbXRIOZPlbTyTTLqvr2omn3H6fbGYICFsGU0bFiA0y+xrb07xV00za6b99poDyghBzqkqU7kLYV6SUbG5K8kFtK64r+90+DQHVZVH2HF0UeJDQGSklmxuTbG5KDWh3OENAdXmEyjxlw3LdAW63FHPoED0X+nMssOCjtP1gEoLOuFKSsj0b+4ppQMQMPhrvCduVtCXdXvdvz0bEEhQGuHFELsRTDtvabZJpN+dOMyPiJTGLspg1aO2TUrKtzWZzY5KUnbvoGcKLEseUhvPWUarAcSXrtrWzriGO4+Z+dr0hvH3fp1QWUFk0eO1L2i6r69rY2JgAyMm+TFMcURhmamUBpTH/VlX0lXjKYdP2JI3tNojcDtQxhBccVhSGGFsWIZLHVRxazHNESokrPVF3ZffRusATb8vHpVn5QEqJ7ULalbhu9/Z5HaW3vGco2edKSSLtkrI9OyU7H7KT+dkyvYmCXipt6NiXdlya4zbxDhvBE+zMc5TSS8dGLEFx1KIwPLgyDb3RlnSoa0nRmnBIpt0u4+pSelF4UdRgZFHYW8o3ROyTUlLfmmJTU4KmuE0i7WIYH20zLaUX3RaGTcoLw4wvj1IUGTpnXzmupLY5SW1zkuZ4mrQjd9gn8ewTQlAUMaksClNVFvV97fpAsB2X+pY0jXGbeMrBdXc+jTPz7hWEDcoKLCqKw0ocaC3m/STzgnXGi3iGRgfSGxnnJYMQw8s+15U7nDKBdxTmUBLvnvAcs47n13GktTXEnK+ekFKStN0dTqchvPHM4WKf40riaQfHlTtOzhtswzwDIe24JNIuboeIh01BxBpe9qUdiZQfbS6l8rCiDLmK+dBxA30i6ANaVBPUrnh+YRhi2M7qFEIEfta0SoQQQypy6yumIYZU5N1XQqY/uwsGhWdf0HfRPcP3m9doNBqNZjdBi7lGo9FoNEMcLeYajUaj0QxxtJhrNBqNRjPE0WKu0Wg0Gs0QR4u5RqPRaDRDnEDEXAjxKSHE+0IIVwjR6/o5jUaj0Wg03RNUZP4ecAbwfED1azQajUYzbAhkBwMp5VIYPruOaTQajUYTJIN+zFwIcbkQYrEQYnFdXV3Qt6PRaDQazaBDWWQuhHgKGJPlo69LKR/JtRwp5d3A3eDtzZ6n29NoNBqNZtigTMyllMeoKluj0Wg0Gs1HDPo0u0aj0Wg0mp4Jamna6UKIGuAQ4FEhxJNB3IdGo9FoNMOBoGazPwQ8FETdGo1Go9EMN3SaXaPRaDSaIY4Wc41Go9FohjhazDUajUajGeJoMddoNBqNZogTyAS4wY6UO+9NM9y2nR3O9u1qGwxP+yQgGF62gWdf5glq+4YeGfsyVg0n+3Z992Bw2afFHHClxHbBld6/XRFIDAGmAaYYXA+wN6SUuBJsV+K44GSxzxASU4BpCCxj6Nlnu5B2XGwnu32mAMuEkGkMSfsSaUnKdkk53jPsjBAQMgRhSxAJGYTMoWMbgONK4imXhO2Sst0u9pkCwpZBJCSIhU0sY2jZl3Zc2pIOybQkabtd+hfLEERCgmjIoDBsYgwh+6SUpGxJe8ohkfba566+dMgURC2DaNigIGwMwXfPJZF2SdoS2/nIEQNP0ENm5vmZRCwRqH27rZhnRC7tZhfwna7FEwnH8X42hSRkgjGIG2ZG5JJZXrBdyTgx6Y4vImRIwpYY1Pa5UnZ0kDu/YNlwJDg2JG0XISBqicBfvN6wXUlb0iGR+igayIaUkHIkKUfSmnSxDCiMmERDg9u+lO3SHHeIp10E3dvnSIinXeJpaGx3iFqC4phFNDR4RwillMTTLk3tNqle2qftSuykpD3psg2bgohBacwibA1u+1qTDk3tDo7bs31pR5J2HFqTDggojpqUxKxB7ZQ5rmdfa8JByh7ePTq9ewkXo8O+wqgZSN+5W4q5KyUpp3cR746MOFgdoj7YOk3HlcTTvYtcd6RdSKckEVMSMgefKKRsSVvK7f3CLEgJ8bQX7RZGBl8kK6WkLenSmuyffbYLzXGHtiSUFVhYg8w+V0q2t9nEU+6O9tmXdpqwJcmWNJGQwYhCC3OQiYLtSOpaUr2K+K5krm1LurQnUxTHTMoLrEH37iXTLltbUrhu356b7Pif5rhDS9yhvNCiOGoOOvvaUw7b2uxeA6BdyQR8jXGH5oRDRVHId4dztxPztCNJ96+f7ILdIeoRSw6KKFZKL1LNl31Jx4scoqHBkYXICF0+7JNAa9IlbDJo0n+O6wndrqnmviLxRL2+1aY4YlAQGRz2JdMu9a3pfjvRGSSQSLtsbkpRUWgRC5t5ub+B0pKw2dZq99uJziCBlrhDe9JhVEl4UETpUkq2t9u0xJ282Le9zaatw77B4JC5UtLQapNMuwO2z5VQ35ImFvYcTr/eveBbiY+k8ijkGSSQsL2OOEiklCTyKOQZHAntKYnbV1c1z7hS0pzIj5B3JuV4op5t4pyf2I6kvtXG7mPE0xstSZeWhBO4ffGUS13LwIW8M1JCQ6tNW8LOX6H9pKk9nRchz5BxyDY3pkjmu9H39V6kpK4lnRch31EmkLQlmxqT2NkmuviI60q2NqdJ5EHIM0igvaPN+/Xu7TZinna8MWRVJB0CE7yMkKuyz2uYwQm6lJKWRNfJQ/nCdoMVdMeVNPQjtZcr7SlJS8JRU3gOJNIuDa3pvDopGSSwvd2hLRmcfc1xm8b2/AldZyRQ25QipbLz6ql+KalvTe80LJJPHBc2NyUDC4ZcKdnakiatyKFI2dI3Qd8txNxx8x+xZiNpZ18apRrVjgpkUpsyEPvaU+qEPIPtevb5jewYQ1b9tbanZCARnuN6YqDSvEzaNogIL5l2veensA4JbGlOBeJMtyQc2pNqhDyD40JdSyqQvmV7m61MyKFjkpwtaYqrdzaHvZjLjsluvtQFvtWVwZUSv4ISR6K04Wcj7fj3/Lzshr/2tSW7LsdSRWPc8V0Q+jOZqD9I8JwGH+2TUrK1JaVU6DK4Lmxv9Xc4wXakckclQzItfc+uJNIu8X5OpO0LEs8pUp1dGfZins7zGGRvONLf8XO/o8mk41/2ITPhzU/afEy32663nMyvJyg7ZhP7RTzl+JoNsB1vSZFfNLbbuD6Z503Y9Pf7rG/1x1EBz76GVtu3vlNKqWzopzsaWm2lfcuwFvPMWmu/8et9c1yZdZMU1aR8qjTl9H95XX/xNtjxp64gxnkTaYnrU4fZlMcJU7kg8ZwVP5wxV0qvLuU1fYTEcyD8IGW7JAMYdmr1aTJje8r1JWPUGcf1VhupYliLeVCTJF3pT/Tqd8r7o3r9sS+IMWyvXh9Sb1Lu2BDGb9p9SC2mbDeQMWwpUdphZmj3OWOUIZ52fYle/XZUwHNWmnxyxvI5Mz9XMul2VQxrMU8HN8FVeXQupT+T+rLWjXpHyXGl8klv3eFt7au28qAcFfBHzFuT/neWoL7DzNAU92csORstiqNXb3grmM5TSvXOdNpxfZ8bkyGh0BkbtmLe+UCDnlj45BPsM3sGe82cxk9//KNur1uy+HWKoxYPPXB/TvWrbiu5lv/fhU+w314zmTtrT372k//X5fPt27fzmbPP4JB5+3DUYQfzwfvv5VSu6ugg15ft6isvY88JYzl03tysnzdu3875nz6Tww7cl2OO6It9Od9qv0jZvY+VJxIJPnH0oXx8/v4cefBcfvKD73a5ZsXyD/nksYczYVQRd93+s5zqdiXKU+25pmgdx+ETRx3EReee3uWzlcuXcdrxRzJ1bAm/uePnOdeteqKRlDLnrJjjOJz0sYO59DNnZP38lRef58SjDuK4w/bjnFOOzalM9WKXe9vIt30SlM8L6EvmJt/tUwhvdrsKhu0OcLl8XY7jcP2XrmbBYwupqq7m8EMO5KRPnsLMWbO6XPeNm2/imOOOz7l+1WKey/vmOA5f/tIXeeTRJ6mqruao+Qdx4idPZsbMj+y79cc/ZK+99+Fv/3yQ5cs+5Mtf+iILnvhv72UrFrtcU7SfOe8CPnfFF/j85y7O+vnPfvJD5uw9l3v+8QDLl33IDdd9kYcf690+25FKt3rNZd5BJBLh/n8vpLCoiHQ6zaknHMXRx57A/gcctOOa8vIRfO///ZzHH30k57oFXocdUbTzljdXJbfn94ff3MHUPafT0tLS5bOy8nK++8NbefKxf/epfrdjEqqqncXSjuxxP/nO/PFuz77WLPY1NzXyzRu+xJ/ue4Sq6j2or9uaU/2qxCBDsg/OkAr7Eoqdsb7s8pbv9iml52zGwvmPo4dtZJ5LX7L49deYPGUqkyZPJhwOc9bZ5/CfBV07xbvuvJ3TTj+DyspRfboHlWM/uURWnn1Tdth35qfO4dEFOze8D5d+wFEfOxqAPafPYN26tWzdsqX3+n1Y950Lhx52BOUjRnT7+bIPl3LkUR/Zt379upzsU5mGk7Lr6WfZEEJQWFQEQDqdJp1Od9kacmTlKPbZbx4hK5R7/aidb5ERu97YvLGGpxc+zqfPy+6Ijawcxdz95mGFcrcNvOhHpX25iunmTTU8898nOKcb+x554D6OP+lUqqr3ADx7cyHjrKgiV7FTZZ9qZyXn56eoffbFWeoLw1bMc9HRTRs3Ul1dvePnqqpqNm/a2OWaBY88zGWXX9n3e+jzX+ROLs1h86aNVFeP3/HzuKoqNu1i3157zeXfjzwEeOK/Yf06Nm6s6bVs1SNO+fKD5uy1Nws67Fuy2LNv06Yc7AtuSHsnHMfhmMPmsde0Ko782MfZb96BeSlX5ZyAXIv+zte/ys3f+QGGkeduSPEwgpPjEN7/ff2r3PTt73dr35pVK2hqbOTTpx7HyR8/lAfuuzen+gVqn1+ufpAq+1S/e7kGWarap6qmOWzFPBeyPdRdI58bvnwdt/zgR5jm4DjMYQc5NIhc7LvuqzfSuH078w/cj9/86g723mdfLGv4jL586cs30tjYyBEH789v77qTvefui2X2bp/K/qQvZZumyVOLFvPG+2t4c8liPvwgtzH/fN6DirKfevIxRo6sZO999gvsHlQW/vTCxxhZOYq95nZvn23bvPfOG/zhbw/x53/+mztu/SGrV63ovfBcc/z9JWD7lAcKOVyjun2qYPj02v2gqrqampqPorSNG2sYM3bcTte88cZiLjzvXAAa6ut58onHsCyLk089zc9b7YLI4YUeV1VNTc2GHT9v2riRsbvYV1JSwl2//QPgif9e06cwYeKkfN9u38lTh1VSUsKdv/k94Nm3z6yp7JGDfSrPOepP2aVlZRx62BE88/RCZsyaM+B7UOnFi8z/9PD8Fr/6Ev994lGeeeoJkskkLS3NfOmKi/jFb/6Un3tQ+ABzKXvJqy/z1BP/8exLJGltbebaz1/MbXf9ccc1Y8dVMaJiJAWFhRQUFnLgIYex9L13mDxlWs+FS7XHLgdtn+ozxnLpO1W2T1X2DdvIPJe5L/vPO4BVK1ewds0aUqkU9//zPk765Ck7XfPB8tUsXbGGpSvWcNoZZ3HbL+/MWchVNspc7Vu9cuUO+x74132c+MmTd7qmsbGRVCoFwJ//8DsOPexwSkpKei1b9al++Zp71tTJvr/86fccOj83+0yFb4YQIqfvr76+jqbGRgDi8TjPP/c/pk6bPvD6AVOhgZYpek2V3vSt7/Hae6t46a3l3PHbv3Do4UflTcgBLIXHalpm78/vhm/ewsvvrGLRG8u4/bd/4dDDjtpJ6ACO/cTJvP7Ki9i2Tby9nbfeeJ2pe87otX6J2vaZy8RPlfapPhI1l7ahsn1aiibWDtvIPJf2YFkWt952O6eedAKO63DBhRcza/Zsfnf3rwH6NU6eQaDWezZzcC8ty+Int/2S00/+BI7jcP6FFzNz1mx+/1vPvks/dyXLPlzKFZdehGmazJg5kzt+/bsc6x+oBT1jGSKnSWiXXfhZXnzhORoa6pk9bQI3fePb2Ok0ABdfdgXLli3lC5+7GNM0mT5jJr/81W9zq1+xgZYhep2ktbV2M1/6/KU4joMrXU457SyOPeEk/vyHuwG48JLL2bqllhM+dggtLc0YwuC3d93Oc6+8TXEvDovKmfqm4Yldf8Y+7/mj93zOv/hzbN1Syyc/Pp/WlmYMw+D3v76Dp196s1fbJGrti1hGv8d17/2TZ99nL/ocU/ecwRFHH8snjjwAwzA457MXMX3m7F7LsAyhtG+JhAxEwumXjfmwLxJS++5FLIOk3b919ANtnwLv+1WBCPqc474wb948uXjx4pyvbw9wYw5TQMRS1yhdKWkLaAcx8GwLK+ww0463b3lQFEcNpdFdS8Lxfd/5zowusZQKwtbmlC87sWXDMgRjy8JK61jfkAhsU6PCiEFlsTr7HFeyYVtSWfk9IYARRRbFUXVxZjzlKD1yuCcEMKokRNjKXdCFEEuklPN6u27YptlBffTYE314Vv3CyDFVqwrV9qkuvycE6ttOVJF3ngthU21kB1AQMZWPfXZbt4I1vLuiYp1wLgigMKJ2Mq5pCKWZjZ6QQCyk1r5IyPD39K1OCKEuazSsxVxxm+gRxcM+AEQCeuFMw3MmVCKEIBzQ84tY6sUuZIpAHBYBFEbVV+yHoHZHUVR9wymNWYE4K0JAzAdHMCj7YiFD+RCXIURg7bMoairrW4a1mBtCBNIgLUPteHnneoJAZXq9M0FFr6rH7DKojrCyIYQ/z88QgoKI/88vGjKUT6ACCFtGINFraUydGHSmIGKon1a+CwIoifnzThTHgskcFSl854e1mAOBRHd+aZAQAr/1wBT+DV946T5/6soQsdRnHTJEQ0LprORslPgkBuB/dCeAsgL/5vRWFIV8tc8UUBzzxz5DCMoL/H1+kZDwzYEPmYavQyUCbx6OSkdz2Iu5aQgUzkPrQsT0JyrPEDKFLyn9DNGQ+hR0ZwrDhm8diuFTCjOD6Ogw/SJi+ddZgvfujSj0RxAyUZ2f0XIkZPiS0gfPvsqSsG+OJkBx1CTkU+cpgJFFYV/7lvJCy7e+0zQ851Ylw17MwRs79+OZWUL9GsldEUIQ8yktHDH9i1ozCCEo9CldWxgxfO1MwFsCV+SDfUJ4KVq/KYiY3lInxfWETEGxT8LamfJCS+mqB/D6rqKo6fuwkxCCyuKQ8om2AqgospSPle+KIYQv2RXPvpDyvmW3EHMhBApXOgBeVBfUhDtDCAoUC3rYpE/LKfJJyFRvX1FE7VK0niiMGMTC6uZ3CKCi0MIIyL6RRRYhS519lkGH6PhvnyEEY8rCyoaeBBANG4woDGZLkJBpMKYkrLRtlhaYFKnuoLshGjIoK1A3fp4Rcj/6zt1CzOEjQVfx0Ezhf3q9yz0Y6gQvbHobLQRJJGQos68oEsxkpgxCCEqiJrFw/u9BiGCinp3vwYvwwnkWdIHn6I0uCQfmqEBmXXvE28wlj+UKvFUBowJyVDJEQgZjSsN5j9C9OQ4mZQV9O3Us3xRFLcoL8y/oAhhZHPJtbH63EXPwvOiold8JXCHDE7sgX7YMpiEoDIu82SeAWEgELuQZIiGD4qiRt3EuU0BJNFghzyCEoCRm5TVKiFiCyoCFPIPRIegleZpF7KWeDUaXhAIV8gyWKRhXHqYwT+vrM5unjAxYyDNEQgZVZREioYE7LJl9HEaXhikNWMgzFEYsRpWEMPMwiX+Hk1ka9nVoZNhu59odQggilrfLUcrp/94BhvBE3O8x5N4whCAW8nZQS/Zvx0LAc1L8WG/dVyxDUBI1iKddknb/y/GclMFnXzRkEC4WNCccEv3YwVDgReMlMf/HWHsj47DEwgbb2mzSdm5Hie5UBp5wjii0Ahv26Q5DCEYWhyiKmtS3pHDc/vUv0ZBgZFF4UDhhnbFMwZiSMK0Jh+3t3g5qfbEvY01hxGREkTXo+s6wZTC2NExT3KY14e3O2B/7SmImxQrXk3fHbifmGUxDEBUSV0Lazf2MWUuANQhFvDNCCMKWIGRKbBdSjszJPs+j9GeHsIEghKAgbBILSZK29y8X+0zhdZShQW6fYQjKCiwcVxJPubSnXFzZ/UFkmd+HTW+yYHgQOimdCZkGo0vCpB2XlrhDe+qjbW13tU90+n1B2KA4ag46Ed+VaMigqjxC0pY0xW3iKXcnOzqTOWJBCG+SW0l0cGRSukMIQXHMoihq0p5yaYrbpGzZvX18dDBMSdQbG/d7knBfEEJQVhCiJCZpSzq0Jhxst4d3T3hnEGQmYBaE/Z9Em2G3FXPwHpwpvIYmpRcluBLcDo9a4P2P0fFPMDjS6bkihLdOO2QKpPQEz5HgdlI+0TED3xRDyzbomAcREkRD3vOzXS/j4kqJ7OggDSGwDG8991CzzzQERVGToqiJ40rSjvcv01a9KNXYsZvcULMvZBqMKDIolxLHhZTjkrLdHRGfN6nUc05UHy6Sbz5qm2Gk9J5bypakHBfpsqNfCVsGEcsYcu3TW2ViUhgxcaUkbUuStusdHtTROE3Dm2sTtvzZyCefGEJQHPX2iHddScqRpGwXx/2ob7FMQdg0CFliUAR3u7WYd0Z07BZnCIblTIIdjgsEu2m9Ijo7LsMR0xBeNmlwDDHmFSEElgmWaVIQ1B6+CslkysIWdLyBwwpDCCIhoew0sKAxDEHU8HePhv4wuO9Oo9FoNBpNr2gx12g0Go1miKPFXKPRaDSaIY4Wc41Go9FohjhazDUajUajGeJoMddoNBqNZogTiJgLIX4ihPhQCPGOEOIhIURZEPeh0Wg0Gs1wIKjI/L/AHCnl3sBy4GsB3YdGo9FoNEOeQMRcSrlQSpnZWfsVoDqI+9BoNBqNZjgwGMbMLwEe7+5DIcTlQojFQojFdXV1Pt6WRqPRaDRDA2XbuQohngLGZPno61LKRzqu+TpgA/d2V46U8m7gboB58+b195AzjUaj0WiGLcrEXEp5TE+fCyEuBD4JfFxKqUVao9FoNJp+EshBK0KIE4AbgSOllO1B3INGo9FoNMOFoMbM7wCKgf8KId4SQvw6oPvQaDQajWbIE0hkLqWcGkS9Go1Go9EMRwbDbHaNRqPRaDQDQIu5RqPRaDRDHC3mGo1Go9EMcbSYazQajUYzxNFirtFoNBrNECeQ2eyDESklrgRHgttpCxsBGOKjf0KIwO6xv0gpkYDrevZJCRkTO9s2lO1zZcY+77+RQIdNpiGGvH2OBMeV2K6k8xZLpgDTNLAEiCFqnysljgu2I3FcuaNtCrxnZ5kC0wBjCNoGnn22C47jPb8MQoBlCM9GY2g+ux39ZkfbdNyPPjMEHc9OYA7Rtpmxz3YltiN30gbv3RNYHf1L0Pbt1mKe6STTzkfilg2n04eGkISMoSEMUkrSDthu9/a5u3xgGZKQOTQ6TldKUrZnY1YkOEC60wMMm5KQJYaEfbYrSdpu9/bh7YWc6UENARFTELbEkGibtguJtIvtdn9d2pUI2xN404CoZRAyh8q7J0mk5U79x654bdO7IGRANGRgDgFhz7x7ibTsse9MdbIvYkIkZGAag9s2ANeVJGyXpN39Nd6759kngIgliISC61t2SzHPdCTpHjqR7nAlJB06Hp4clKKQEfH+2Ge73j9DSCLW4BR1Kb1OpCcR6I6U43UwIVMSGaSi57iS9pTbowhkw5UQtyVxWxK1xKC1z3YkbSm3iyPZHZnLHBfaUi4CKIwYhMzBZ5vsELn2dN93qE67kE66GMKzzxqEopd59xJ23+1LOpB0XCzDs28w9i2u9N69nhzobEggYXvfS9iEgrDh+7u324m5lJKk0zUi7XM5QMLuiGQHkSftuJKk3XOmIRdcCfG0F8kOphSg7Uji/egodyXteGVFwwyaTlNKSdLuX0e5KwnbiwwLwoMnEpIdHWWqjx1ll3KA1qRLqKPTHCyi4EpJW7LnTENu5UBLwiVqCaKhweOQ2a5n30D7TtuFprhLQWhwZZHSjmffQN++lAPphEtR2MDy0eHcrSbAuVIStwcu5J2xXc/jHAxnxdiuJJEHIe+MF8kODvuSaTcvQp5BAvGUJDXQ3jcf9yK9aDUfQp7BkdCSdLH7GuIrwJWS5sTAhbwzaQeaEy5uPl/ofuK4kub4wIW8Mwlb0pJ0B8W7l3IkLYmBC3ln2tOSeGpw2BdPubTmQcgzyI53L9Gf9Gg/2W3E3JWe0KkpO3hBtzsicjVlBy/oyXR+hWCnsm0CFfSMkKu6hdaUu9PEK79xZf6FIIOUwQu643qOioo7cFwCF/RMxKqCpAPtAQt6PM9O9E5lp6Vvgr5biLlUKOQZMoIeBK5CIc9gu/R5HClfpGx1Qp4haRNYBBvvZRJYPmhNurgBdJhSyo66FdZBcIInOxwVlTguysS097q956eSlIMyMe2NZFqdkGeIp+VOk3BVsVuIuWohyJBZwuAnfjgqGdIuvguCK9U7KhkSaem7INiO9K19tgUQASXtnZcrqcKb4+G/4LWl1ETku5J2IeWz4GUcMT9IpL1liX7iyv5NVOwPbT44m8NezB2356Uh+cbvdHRvy+ryTSLtr32JlH91SSDp08sNH6XX/cLxWRAcNz+TFXMlafvrTKedHpZFKqA95W92JZGWSjMqu9Lqc3bFz2yHxHt+KhnWYi6lf1FPZ/x6waWU/Vp+NqA6wTfnyG9HDPzNPiTtntfoqiBh+5d9CCJSjvvkHGVm5vuJBN/GX72Mn9+Rsn99p+30b2nrQEg5KJ3bMazF3JX+Rq0ZbOlP9BrUGLZf9fqdVvSzXs/RDGCMl/7tP9BXXOlv1JrBdvElXeu4+V0Vkysp25++JRnQu5fwSWH9qsfPeoe1mAcQGOxAdT+d2fgmCFypPnoN0r60D0MlQYkBeJN+VBOUIwaQ9KHh+LnkqDMS9c50Zr+DIHB8cMaCcjTBGwpS1bcMWzHP7KnbGwuffIJ9Zs9gr5nT+OmPf9Tl8//8+xEO3G8uB8/bl8MOPoCXXlyUU/2qJ/18tEliz/z3ySfYd84M9p45jVt/0tW+2279CYccsC+HHLAvB+y7FyUxi23btvVarmr7cu2Pv3DFpUzeYwwH7b931s+XL/uQjx85n5GlMX7581tzrl+10KZz+AKvvuIypk0YyyHz5mb9fNHzz7LHmBEcftD+HH7Q/vz4B7fkVLfjQ+Yol6zD1Vdexp4TxnJoN/b98uc/5YiD9+eIg/fn0HlzGVkcZnsObdMPsctFy3uzr7mpiXPPOpXDD9qPQ+btzb1/+VNO9acUv3xS5tb+VdmnelWJnWP7UNU+VT0+MRgW7OfKvHnz5OLFi3O6Npd15Y7jMHf2dBY8tpCq6moOP+RA/nTP35g5a9aOa1pbWyksLEQIwbvvvMMFnzmHN99b2mv9AoiF1O3+k8u6csdx2Gf2dP7dYd8Rhx7IH+/5GzNnzsp6/WP/WcAdt9/GY08+3Wv9poCoQvsS6dy2VHxx0fMUFhZxxWUX8eqSd7p8Xrd1K+vXr+PRBY9QVlbONdd9Oaf6I5a3O5UqWpJOry/1i4uep6iwiCs/dzEvL367y+eLnn+W22/7Gfc9+O8+11+kcLtQKSWN8d57rJc6nt3nP3cxL2WxrzNPPLaAu27/BY88/lRO91AWU7edputKmnJYjtabfT/7yQ9pbmriO9/7EfV1dRy47yw+XL2RcDjcY7mGgNKY2e/7742UndvETFX2hUwoiqizrz3l5LRCRlX7LAgJIqHc42ghxBIp5bzerhu2kXkunuXi119j8pSpTJo8mXA4zFlnn8N/Fjyy0zVFRUU7OoX29racOwiJ2ugnF+8um32P7mJfZ/71z3/wqbM/nVP9qiPXXL3X+YcdQfmIEd1+XjlqFPvPOwArFOpj/eoMlDK35Vq92TYQ1NqX23WH9sG+B/55H2fk2DZBbeYo16xRb/YJBK2trd6qhrZWystHYFm977DtKs6s2G5uBqqyb7Bk/VS1z7Sid2/YinkubX3Txo1UV1fv+LmqqprNmzZ2ue7fDz/EvnNmcuapn+Su3/4+93vI+cq+k5N9mzZSPX5n+zZt7GofQHt7O08tfIJTTz8zt/pzuqr/BJ0wGgQ7hObE66+9wmEH7cdZp57E0g/ez/nvVM6qzXfR7e3tPP3Uk5xy6hmB3cPOZeen8MuuvIrly5Yya8p4DjtwH374k59hGLl1yUrty5OY9tc+1e9e0O0zX9/vrgxbMc+FbN5ttsj7lNNO5833lvKP+x/i/77zLT9uLS/kah/AY48u4OBD5jNCUSSoyT9777Mf73y4mkWvvsHln7+K887JzREDxY5mnst74rH/cNDBh/YpS+H/or++87+nFjJnr7l8sGoDz728hBuu/xLNzc1B31bevrnBal++m0Zf26eqlrlbi3lVdTU1NTU7ft64sYYxY8d1e/1hhx/BmtWrqK+v9+P2BkxVVTU1G3a2b+y47Pbd/8/7+NQ5uacxNcFTUlJCUVERAMedcCLpdJqGHNumyrOc8l32Q/ffx5mf6lvbFEotzA9/u+dPnHzq6QghmDxlKhMmTGTF8g+Dvq28fXOD1b5809f2qaplDlsxz2Voe/95B7Bq5QrWrllDKpXi/n/ex0mfPGWna1atXLkjwn3zzTdIpVJUVFTkdg99vuvc6a99J+5iH0BTUxMvvvAcJ518au719+Vm+0HQpyIGXX8ubKmt3dE2l7z+Gq7rMiLXtqnQwHwW3dzUxIuLnucTWdqtX/egquzq8Xvw3LP/A2Drli2sXLGciRMn5/S3Kk+1zVfZ/bVvKPUt/WmfqtrmsD3PPJcGaVkWt952O6eedAKO63DBhRcza/Zsfnf3rwG47PIrefihB/j7X+/BCoWIxWL85d5/5NwRquwwTQG9TcjM2HfaJ0/AcRzOv+hiZs3a2T6ABY88xNHHHEdhYWHO9as+Its0wM1hNvvFF3yGRS88R0N9PTOm7MHN3/w26XQagEs/dyVbams5cv6BtLQ0YxgGv7rjF7z25nuUlJT0XL9SMRAYovexu0sv/CwvPv8cDQ31zJ46gZu+8ZFtl3zuCh556AH++LvfYFoWsWiU3//l3pzbnMozznMt+rILP8uLL3TYN82zz+6w7+LLrgDgP/9+mI99/Ng+tU3w2o8qvFUAvSdLe7PvKzd9nasuv4T5B+yDlJJv3/JDKkaO7LVcgeK+xRQ5bZShyj6Vzw4gZOR2KJaq9qlqFcmwXZomO84uDwpDQFTh0iZXSuJpZcX3SshA6dKttCNJ+Liv967EQgJLoaK3p5xAthrOUBpVt3QLoCnuBDaJUABlBeqWNuW69E4VlgHFUXX2ZY50DYqoJYiF1Sl60nZp9/HMh10pihiE+tC37PZL0zLRT1CojOxAfSqqN1R7z6q/v17rV2yfpbqCHjCE2sgOvLXCQWEprlsIobx99ERfhKA/BNlvgnr7VEXGudevptxhK+ag7ksbDHULIQKzT6D+hTeM4Jwxy/BB7AJsm2EfPKVIgC9f1Ie6/aijO1RmxMBr+5GABmAN4UOgYATnjIVMdX3LsBbzoKI704fIB4KLflQ2yM5EFHda3aG6s4SODjOg5+fH92oawTibfogBdLwD6qvpQtgEw4d3LyhnJWIJX/qWoOyLKvTih7WYCyECiYD8EllDiEAcFr/eA9Pwf1a5Jwb+VNqXLR3zVqdPjhhALAD7YiG1cwEyCCGUbmfcHSrFoDOG4X/fKfDPgQ+Z/g8nmIbaFP+wFnPoSJn6XJ8fnnMGv9NhEcs/MRBCKN3fPht+1mcIQczH7IPAPzEAsP5/e3cfYlldx3H8/Tnn3IeZddcglcqVNBBJNrMIMfynzGQr2agIjJ6gQIIEo6KypSJCEATrDwOJiqCsCEoKJdTIkKAHy3xsNUQqtweWCCsfdmfv3G9/nHN3ruvMzszOued3z5nPCxa9u3fH79c7cz6/8zu/8zu56Dc4+1Bkzc5WDYpmB9PDQo0NNAEWB1mjx84dg2YGYlAeW3YMmo2/HTNc1AfbIMybvP4jmr8W2mR/uZpfPJJnauwAPSjKM5Im9Yvmrt81ebCcWOw3FwhN99dkIGQzfrDR6v9NsTjjAJro5bNf+Ha8Imvu2LnQm/1ArPNhDuU3ZRNnCE2etU4rGpgSy2h+FmBi0EDglWd1zX92ktjRz2Y+5bdQKMkqXknsHM7+MHPKIGt0RmwizzTzMy4BOxMMxKAcbM7yFlsoTxJm/f9wLQu9bOaXDft5M5cPtkWYw+wDb1g0O71+vF4+21sehr00AxVYmW6fVdYWWXnWk6q/TKrCaDZff1hs7pGLdcuz2Qb6Zu/brVu/EIszOmsWsHOYNT5jNG3Ym90ZbK6yv5THlvJxwLP5+r28mp1qoL9tE+ZQnnnVvYK43BwmbZDDZLq9/hmIMujSBfmEJBb69U+59/O0QT6RzeigstjPGr1OvpYiE7uG9Q5YVAVByiCfGPQyTqn5GnOewa5h1uh18tVIYrGf176epJ+nDfKJSaDXPWAZFuWsTVP9dXY717XkmVhQsLS8oR0LT6iXNXNP8mb0cpFnwZHR1h/1NyyaW9m9EZMVxL08eG5pa8/FmlyDnKf+ykDPWVoOnl3a2g5cRVYGeepB5rS8CvTDR4PDo619cw6K5laub1QvF6cuZDy7NN7y7n4LPTV2m9ZGDXsZvTx4Zmm8pWeOi3J9wzwMwiYmA5Z+Hjy9NN7SI5gzlbNFTR9btl2Yw8qisXEER08i1IusDPJ5+kGblkkMi2AcsLS8uVAXK1P289pfnokdAxiNYWkUm+ov08o1+Hntr5+L3jDjyCg4MtrcoKWXlWeJqXe5WstkhmVQBIdHY45scsvlQVFuSDNPg7Bp5aK4nME4OHJ0c6E+uTVr0NNcDcKm5ZnYOcgYjeHw0TGjTYT6ZADdz+drkDKtyMWpw+zYdtKbyYY8K+9fb2ofjuNtyzCfyKpQjyg/tPG4DPbjP79M1f3HamYrzDqougd9ISsHLcvjMtSXx8/vb7KbW56t9NmW/iYrYJfHZX+Tfz6vv+pzm2xikvLa42ZMZiEGRdnTaByMVukvUznwKrJyL/l5DYHjZdlk6jYYjWG0/ML+BGTZpL90B8mTUWSiGOQsVicMo/GYUTWwnu4vn/rs5nkAPW3lZy9nPA6OjqP6/CCm+xMUKgOyyNXYZlpbJYl+IfpFeUw5OvW9OTlxEFV/WdVfln6Wb1uH+YQkCtHZFQSZRJZwr+xZW9mecf4PFJsliSJnpg99SWl6UNZFZTBAn27+AGaZGDR4i1fT8jkI6Y3qaHyZmZltHw5zMzOzlnOYm5mZtZzD3MzMrOUc5mZmZi3nMDczM2u5JGEu6UuSHpR0v6Q7Jb0sRR1mZmZdkOrM/IaIuCAiLgRuAz6fqA4zM7PWSxLmEfHfqZc7eOGma2ZmZrZByfbtkXQd8AHgP8AbU9VhZmbWdoqtPB7mRF9Y+hnwklX+aH9E/HjqfdcCw4j4whpf5yrgqurlHuDhumudI6cB/0pdxAx1ub8u9wbur+3cX3udFxE713vTzMJ8oyS9HLg9IvZs4L2/i4jXNVBWEu6vvbrcG7i/tnN/7bXR3lKtZj936uU+4NEUdZiZmXVBqmvm10s6DxgDfwE+kqgOMzOz1ksS5hHxrpP8q1+rtZD54/7aq8u9gftrO/fXXhvqLfk1czMzM9sab+dqZmbWcq0Nc0mflBSSTktdS126vs2tpBskPVr1eKukF6WuqU6S3i3pEUljSZ1ZWStpr6THJD0u6TOp66mTpG9KOiSpc7e8SjpL0t2SDlTfl9ekrqlOkoaSfivpgaq/L6auaRYk5ZL+IOm2E72vlWEu6SzgzcBfU9dSs65vc3sXsCciLgD+BFybuJ66PQy8E7gndSF1kZQDXwXeApwPvEfS+WmrqtW3gL2pi5iREfCJiHglcDHw0Y59dkeASyPi1cCFwF5JF6ctaSauAQ6s96ZWhjnwZeBTdGwb2K5vcxsRd0bEqHr5a2B3ynrqFhEHIuKx1HXU7CLg8Yh4IiKWgO8Db09cU20i4h7g36nrmIWI+EdE3Ff9+/8oA+HMtFXVJ0pPVy971a9OHTMl7QbeBnx9vfe2Lswl7QP+FhEPpK5lFiRdJ+lJ4L1078x82oeAn6YuwtZ1JvDk1OuDdCgQtgtJZwOvAX6TuJRaVVPQ9wOHgLsiolP9AV+hPHEdr/fGZHuzn8iJtoIFPgtc3mxF9Vlvm9uI2A/sr7a5vRpYdZvbebWRbXwl7aecArylydrqsNFtijtEq/xep85+uk7SKcAPgY8dN/vXehGxDFxYrb+5VdKeiOjE+gdJVwCHIuL3kt6w3vvnMswj4rLVfl/Sq4BzgAckQTlNe5+kiyLinw2WeNLW6m0V3wVup2Vhvl5/kj4IXAG8KVp4X+QmPr+uOAicNfV6N/D3RLXYJknqUQb5LRHxo9T1zEpEPCXpF5TrHzoR5sAlwD5JbwWGwC5J34mI96325lZNs0fEQxFxRkScHRFnUx5oXtuWIF9P17e5lbQX+DSwLyKeTV2Pbci9wLmSzpHUB64EfpK4JtsAlWc83wAORMSNqeupm6TTJ3fESFoALqNDx8yIuDYidldZdyXw87WCHFoW5tvA9ZIelvQg5aWETt1KAtwE7ATuqm6/uzl1QXWS9A5JB4HXA7dLuiN1TVtVLVi8GriDcgHVDyLikbRV1UfS94BfAedJOijpw6lrqtElwPuBS6uft/urs7yueClwd3W8vJfymvkJb9/qMu8AZ2Zm1nI+MzczM2s5h7mZmVnLOczNzMxazmFuZmbWcg5zMzOzlnOYm9kxkparW5geqZ5G9XFJWfVnL66ewvW0pJtS12pmK+ZyBzgzS+a56ql9SDqDcifCUyl3IjwMfA7YU/0ysznhM3MzW1VEHAKuAq6WpIh4JiJ+SRnqZjZHHOZmtqaIeILyOHFG6lrMbG0OczNbz2pPTjOzOeIwN7M1SXoFsEz5vGgzm1MOczNblaTTgZuBm9r4uFqz7cQPWjGzYyQtAw8BPWAEfBu4MSLG1Z//GdgF9IGngMsj4o9JijWzYxzmZmZmLedpdjMzs5ZzmJuZmbWcw9zMzKzlHOZmZmYt5zA3MzNrOYe5mZlZyznMzczMWs5hbmZm1nL/B+gcPzoIy8E0AAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 576x432 with 1 Axes>"
]
},
"metadata": {
"needs_background": "light"
},
"output_type": "display_data"
}
],
"source": [
"fig, ax = pyplot.subplots(figsize=(8, 6))\n",
"\n",
"ax.scatter(X[:, 1], X[:, 0], color=cm.Blues(Y / Y.max()), s=1000)\n",
"for (y, x), v in zip(X, Y):\n",
" ax.text(x, y, f\"{v:.1f}\", ha=\"center\", va=\"center\")\n",
"ax.set(\n",
" ylabel=\"D0\",\n",
" xlabel=\"D1\",\n",
" ylim=(-3, 3),\n",
" xlim=(-4, 4),\n",
")\n",
"pyplot.show()"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "5efac6a7",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"<?xml version=\"1.0\" encoding=\"UTF-8\" standalone=\"no\"?>\r\n",
"<!DOCTYPE svg PUBLIC \"-//W3C//DTD SVG 1.1//EN\"\r\n",
" \"http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd\">\r\n",
"<!-- Generated by graphviz version 2.47.2 (0)\r\n",
" -->\r\n",
"<!-- Pages: 1 -->\r\n",
"<svg width=\"384pt\" height=\"363pt\"\r\n",
" viewBox=\"0.00 0.00 384.45 362.86\" xmlns=\"http://www.w3.org/2000/svg\" xmlns:xlink=\"http://www.w3.org/1999/xlink\">\r\n",
"<g id=\"graph0\" class=\"graph\" transform=\"scale(1 1) rotate(0) translate(4 358.86)\">\r\n",
"<polygon fill=\"white\" stroke=\"transparent\" points=\"-4,4 -4,-358.86 380.45,-358.86 380.45,4 -4,4\"/>\r\n",
"<g id=\"clust2\" class=\"cluster\">\r\n",
"<title>cluster48</title>\r\n",
"<path fill=\"none\" stroke=\"black\" d=\"M152,-8C152,-8 246,-8 246,-8 252,-8 258,-14 258,-20 258,-20 258,-334.86 258,-334.86 258,-340.86 252,-346.86 246,-346.86 246,-346.86 152,-346.86 152,-346.86 146,-346.86 140,-340.86 140,-334.86 140,-334.86 140,-20 140,-20 140,-14 146,-8 152,-8\"/>\r\n",
"<text text-anchor=\"middle\" x=\"243\" y=\"-15.8\" font-family=\"Times New Roman,serif\" font-size=\"14.00\">48</text>\r\n",
"</g>\r\n",
"<g id=\"clust1\" class=\"cluster\">\r\n",
"<title>cluster2</title>\r\n",
"<path fill=\"none\" stroke=\"black\" d=\"M20,-232.91C20,-232.91 120,-232.91 120,-232.91 126,-232.91 132,-238.91 132,-244.91 132,-244.91 132,-334.86 132,-334.86 132,-340.86 126,-346.86 120,-346.86 120,-346.86 20,-346.86 20,-346.86 14,-346.86 8,-340.86 8,-334.86 8,-334.86 8,-244.91 8,-244.91 8,-238.91 14,-232.91 20,-232.91\"/>\r\n",
"<text text-anchor=\"middle\" x=\"120.5\" y=\"-240.71\" font-family=\"Times New Roman,serif\" font-size=\"14.00\">2</text>\r\n",
"</g>\r\n",
"<!-- ls -->\r\n",
"<g id=\"node1\" class=\"node\">\r\n",
"<title>ls</title>\r\n",
"<ellipse fill=\"none\" stroke=\"black\" cx=\"70\" cy=\"-301.38\" rx=\"54.39\" ry=\"37.45\"/>\r\n",
"<text text-anchor=\"middle\" x=\"70\" y=\"-312.68\" font-family=\"Times New Roman,serif\" font-size=\"14.00\">ls</text>\r\n",
"<text text-anchor=\"middle\" x=\"70\" y=\"-297.68\" font-family=\"Times New Roman,serif\" font-size=\"14.00\">~</text>\r\n",
"<text text-anchor=\"middle\" x=\"70\" y=\"-282.68\" font-family=\"Times New Roman,serif\" font-size=\"14.00\">Lognormal</text>\r\n",
"</g>\r\n",
"<!-- f -->\r\n",
"<g id=\"node5\" class=\"node\">\r\n",
"<title>f</title>\r\n",
"<polygon fill=\"none\" stroke=\"black\" points=\"245.5,-213.93 154.5,-213.93 154.5,-160.93 245.5,-160.93 245.5,-213.93\"/>\r\n",
"<text text-anchor=\"middle\" x=\"200\" y=\"-198.73\" font-family=\"Times New Roman,serif\" font-size=\"14.00\">f</text>\r\n",
"<text text-anchor=\"middle\" x=\"200\" y=\"-183.73\" font-family=\"Times New Roman,serif\" font-size=\"14.00\">~</text>\r\n",
"<text text-anchor=\"middle\" x=\"200\" y=\"-168.73\" font-family=\"Times New Roman,serif\" font-size=\"14.00\">Deterministic</text>\r\n",
"</g>\r\n",
"<!-- ls&#45;&gt;f -->\r\n",
"<g id=\"edge3\" class=\"edge\">\r\n",
"<title>ls&#45;&gt;f</title>\r\n",
"<path fill=\"none\" stroke=\"black\" d=\"M98.09,-269.24C109.4,-257.36 122.86,-244.02 136,-232.91 141.14,-228.56 146.69,-224.23 152.31,-220.06\"/>\r\n",
"<polygon fill=\"black\" stroke=\"black\" points=\"154.5,-222.8 160.54,-214.1 150.39,-217.13 154.5,-222.8\"/>\r\n",
"</g>\r\n",
"<!-- scaling -->\r\n",
"<g id=\"node2\" class=\"node\">\r\n",
"<title>scaling</title>\r\n",
"<ellipse fill=\"none\" stroke=\"black\" cx=\"322\" cy=\"-301.38\" rx=\"54.39\" ry=\"37.45\"/>\r\n",
"<text text-anchor=\"middle\" x=\"322\" y=\"-312.68\" font-family=\"Times New Roman,serif\" font-size=\"14.00\">scaling</text>\r\n",
"<text text-anchor=\"middle\" x=\"322\" y=\"-297.68\" font-family=\"Times New Roman,serif\" font-size=\"14.00\">~</text>\r\n",
"<text text-anchor=\"middle\" x=\"322\" y=\"-282.68\" font-family=\"Times New Roman,serif\" font-size=\"14.00\">Lognormal</text>\r\n",
"</g>\r\n",
"<!-- scaling&#45;&gt;f -->\r\n",
"<g id=\"edge2\" class=\"edge\">\r\n",
"<title>scaling&#45;&gt;f</title>\r\n",
"<path fill=\"none\" stroke=\"black\" d=\"M295.88,-268.12C285.82,-256.49 273.88,-243.61 262,-232.91 257.22,-228.6 252.02,-224.33 246.74,-220.23\"/>\r\n",
"<polygon fill=\"black\" stroke=\"black\" points=\"248.5,-217.18 238.41,-213.96 244.29,-222.77 248.5,-217.18\"/>\r\n",
"</g>\r\n",
"<!-- sigma -->\r\n",
"<g id=\"node3\" class=\"node\">\r\n",
"<title>sigma</title>\r\n",
"<ellipse fill=\"none\" stroke=\"black\" cx=\"320\" cy=\"-187.43\" rx=\"54.39\" ry=\"37.45\"/>\r\n",
"<text text-anchor=\"middle\" x=\"320\" y=\"-198.73\" font-family=\"Times New Roman,serif\" font-size=\"14.00\">sigma</text>\r\n",
"<text text-anchor=\"middle\" x=\"320\" y=\"-183.73\" font-family=\"Times New Roman,serif\" font-size=\"14.00\">~</text>\r\n",
"<text text-anchor=\"middle\" x=\"320\" y=\"-168.73\" font-family=\"Times New Roman,serif\" font-size=\"14.00\">Lognormal</text>\r\n",
"</g>\r\n",
"<!-- L -->\r\n",
"<g id=\"node6\" class=\"node\">\r\n",
"<title>L</title>\r\n",
"<ellipse fill=\"lightgrey\" stroke=\"black\" cx=\"204\" cy=\"-76.48\" rx=\"41.94\" ry=\"37.45\"/>\r\n",
"<text text-anchor=\"middle\" x=\"204\" y=\"-87.78\" font-family=\"Times New Roman,serif\" font-size=\"14.00\">L</text>\r\n",
"<text text-anchor=\"middle\" x=\"204\" y=\"-72.78\" font-family=\"Times New Roman,serif\" font-size=\"14.00\">~</text>\r\n",
"<text text-anchor=\"middle\" x=\"204\" y=\"-57.78\" font-family=\"Times New Roman,serif\" font-size=\"14.00\">Normal</text>\r\n",
"</g>\r\n",
"<!-- sigma&#45;&gt;L -->\r\n",
"<g id=\"edge5\" class=\"edge\">\r\n",
"<title>sigma&#45;&gt;L</title>\r\n",
"<path fill=\"none\" stroke=\"black\" d=\"M288.3,-156.66C273.31,-142.58 255.35,-125.71 239.87,-111.17\"/>\r\n",
"<polygon fill=\"black\" stroke=\"black\" points=\"242.17,-108.53 232.49,-104.23 237.38,-113.63 242.17,-108.53\"/>\r\n",
"</g>\r\n",
"<!-- f_rotated_ -->\r\n",
"<g id=\"node4\" class=\"node\">\r\n",
"<title>f_rotated_</title>\r\n",
"<ellipse fill=\"none\" stroke=\"black\" cx=\"199\" cy=\"-301.38\" rx=\"50.82\" ry=\"37.45\"/>\r\n",
"<text text-anchor=\"middle\" x=\"199\" y=\"-312.68\" font-family=\"Times New Roman,serif\" font-size=\"14.00\">f_rotated_</text>\r\n",
"<text text-anchor=\"middle\" x=\"199\" y=\"-297.68\" font-family=\"Times New Roman,serif\" font-size=\"14.00\">~</text>\r\n",
"<text text-anchor=\"middle\" x=\"199\" y=\"-282.68\" font-family=\"Times New Roman,serif\" font-size=\"14.00\">Normal</text>\r\n",
"</g>\r\n",
"<!-- f_rotated_&#45;&gt;f -->\r\n",
"<g id=\"edge1\" class=\"edge\">\r\n",
"<title>f_rotated_&#45;&gt;f</title>\r\n",
"<path fill=\"none\" stroke=\"black\" d=\"M199.33,-263.73C199.44,-251.07 199.57,-236.92 199.68,-224.32\"/>\r\n",
"<polygon fill=\"black\" stroke=\"black\" points=\"203.18,-224.21 199.77,-214.18 196.18,-224.15 203.18,-224.21\"/>\r\n",
"</g>\r\n",
"<!-- f&#45;&gt;L -->\r\n",
"<g id=\"edge4\" class=\"edge\">\r\n",
"<title>f&#45;&gt;L</title>\r\n",
"<path fill=\"none\" stroke=\"black\" d=\"M200.94,-160.89C201.34,-149.98 201.82,-136.89 202.28,-124.35\"/>\r\n",
"<polygon fill=\"black\" stroke=\"black\" points=\"205.79,-124.12 202.66,-114 198.79,-123.87 205.79,-124.12\"/>\r\n",
"</g>\r\n",
"</g>\r\n",
"</svg>\r\n"
],
"text/plain": [
"<graphviz.dot.Digraph at 0x222872b7b48>"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"with pm.Model() as pmodel:\n",
" # GP hyperparameters\n",
" ls = pm.Lognormal('ls', mu=np.log(4), sd=0.1, shape=(2,))\n",
" scaling = pm.Lognormal('scaling', mu=np.log(10), sd=0.1)\n",
"\n",
" # construct the GP\n",
" mean_func = pm.gp.mean.Zero()\n",
" cov_func = scaling**2 * pm.gp.cov.ExpQuad(\n",
" input_dim=2,\n",
" ls=ls\n",
" )\n",
" gp = pm.gp.Latent(mean_func=mean_func, cov_func=cov_func)\n",
" \n",
" # Add observations\n",
" f = gp.prior('f', X=X) \n",
" sigma = pm.Lognormal('sigma', mu=np.log(1))\n",
" pm.Normal(\n",
" 'L',\n",
" mu=f, sigma=sigma,\n",
" observed=Y - np.mean(Y)\n",
" )\n",
"pm.model_to_graphviz(pmodel)"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "570ca92a",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
" <div>\n",
" <style>\n",
" /* Turns off some styling */\n",
" progress {\n",
" /* gets rid of default border in Firefox and Opera. */\n",
" border: none;\n",
" /* Needs to be in here for Safari polyfill so background images work as expected. */\n",
" background-size: auto;\n",
" }\n",
" .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n",
" background: #F44336;\n",
" }\n",
" </style>\n",
" <progress value='5001' class='' max='5001' style='width:300px; height:20px; vertical-align: middle;'></progress>\n",
" 100.00% [5001/5001 00:16<00:00 logp = 16.433, ||grad|| = 1,377.2]\n",
" </div>\n",
" "
],
"text/plain": [
"<IPython.core.display.HTML object>"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n"
]
}
],
"source": [
"with pmodel:\n",
" map_ = pm.find_MAP()"
]
},
{
"cell_type": "markdown",
"id": "51c11498",
"metadata": {},
"source": [
"### Works with 5 points"
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "6568d873",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
" <div>\n",
" <style>\n",
" /* Turns off some styling */\n",
" progress {\n",
" /* gets rid of default border in Firefox and Opera. */\n",
" border: none;\n",
" /* Needs to be in here for Safari polyfill so background images work as expected. */\n",
" background-size: auto;\n",
" }\n",
" .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n",
" background: #F44336;\n",
" }\n",
" </style>\n",
" <progress value='500' class='' max='500' style='width:300px; height:20px; vertical-align: middle;'></progress>\n",
" 100.00% [500/500 00:01<00:00]\n",
" </div>\n",
" "
],
"text/plain": [
"<IPython.core.display.HTML object>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"name = f\"x_{np.random.randint(300)}\"\n",
"x_dense = X_denses[\"D1\"][:5]\n",
"\n",
"with pmodel:\n",
" gp.conditional(name, x_dense)\n",
" \n",
" mapp = pm.sample_posterior_predictive(\n",
" [map_]*500,\n",
" var_names=[name],\n",
" random_seed=None\n",
" )"
]
},
{
"cell_type": "markdown",
"id": "e683ffbf",
"metadata": {},
"source": [
"### Crashes with all of them"
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "5eb98c09",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
" <div>\n",
" <style>\n",
" /* Turns off some styling */\n",
" progress {\n",
" /* gets rid of default border in Firefox and Opera. */\n",
" border: none;\n",
" /* Needs to be in here for Safari polyfill so background images work as expected. */\n",
" background-size: auto;\n",
" }\n",
" .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n",
" background: #F44336;\n",
" }\n",
" </style>\n",
" <progress value='0' class='' max='500' style='width:300px; height:20px; vertical-align: middle;'></progress>\n",
" 0.00% [0/500 00:00<00:00]\n",
" </div>\n",
" "
],
"text/plain": [
"<IPython.core.display.HTML object>"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"ename": "LinAlgError",
"evalue": "Matrix is not positive definite",
"output_type": "error",
"traceback": [
"\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[1;31mLinAlgError\u001b[0m Traceback (most recent call last)",
"\u001b[1;32m<ipython-input-8-6a2976ea130b>\u001b[0m in \u001b[0;36m<module>\u001b[1;34m\u001b[0m\n\u001b[0;32m 8\u001b[0m \u001b[1;33m[\u001b[0m\u001b[0mmap_\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m*\u001b[0m\u001b[1;36m500\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 9\u001b[0m \u001b[0mvar_names\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mname\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m,\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 10\u001b[1;33m \u001b[0mrandom_seed\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;32mNone\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 11\u001b[0m )\n",
"\u001b[1;32m~\\AppData\\Local\\Continuum\\miniconda3\\envs\\murefi_env\\lib\\site-packages\\pymc3\\sampling.py\u001b[0m in \u001b[0;36msample_posterior_predictive\u001b[1;34m(trace, samples, model, var_names, size, keep_size, random_seed, progressbar)\u001b[0m\n\u001b[0;32m 1731\u001b[0m \u001b[0mparam\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0m_trace\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0midx\u001b[0m \u001b[1;33m%\u001b[0m \u001b[0mlen_trace\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 1732\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m-> 1733\u001b[1;33m \u001b[0mvalues\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mdraw_values\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mvars_\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mpoint\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mparam\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0msize\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0msize\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 1734\u001b[0m \u001b[1;32mfor\u001b[0m \u001b[0mk\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mv\u001b[0m \u001b[1;32min\u001b[0m \u001b[0mzip\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mvars_\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mvalues\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 1735\u001b[0m \u001b[0mppc_trace_t\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0minsert\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mk\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mname\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mv\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0midx\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
"\u001b[1;32m~\\AppData\\Local\\Continuum\\miniconda3\\envs\\murefi_env\\lib\\site-packages\\pymc3\\distributions\\distribution.py\u001b[0m in \u001b[0;36mdraw_values\u001b[1;34m(params, point, size)\u001b[0m\n\u001b[0;32m 789\u001b[0m \u001b[1;31m# This may fail for autotransformed RVs, which don't\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 790\u001b[0m \u001b[1;31m# have the random method\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 791\u001b[1;33m \u001b[0mvalue\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0m_draw_value\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mnext_\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mpoint\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mpoint\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mgivens\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mtemp_givens\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0msize\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0msize\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 792\u001b[0m \u001b[0mgivens\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mnext_\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mname\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;33m(\u001b[0m\u001b[0mnext_\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mvalue\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 793\u001b[0m \u001b[0mdrawn\u001b[0m\u001b[1;33m[\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mnext_\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0msize\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m]\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mvalue\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
"\u001b[1;32m~\\AppData\\Local\\Continuum\\miniconda3\\envs\\murefi_env\\lib\\site-packages\\pymc3\\distributions\\distribution.py\u001b[0m in \u001b[0;36m_draw_value\u001b[1;34m(param, point, givens, size)\u001b[0m\n\u001b[0;32m 953\u001b[0m \u001b[1;32mreturn\u001b[0m \u001b[0mpoint\u001b[0m\u001b[1;33m[\u001b[0m\u001b[0mparam\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mname\u001b[0m\u001b[1;33m]\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 954\u001b[0m \u001b[1;32melif\u001b[0m \u001b[0mhasattr\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mparam\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;34m\"random\"\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;32mand\u001b[0m \u001b[0mparam\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mrandom\u001b[0m \u001b[1;32mis\u001b[0m \u001b[1;32mnot\u001b[0m \u001b[1;32mNone\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 955\u001b[1;33m \u001b[1;32mreturn\u001b[0m \u001b[0mparam\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mrandom\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mpoint\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mpoint\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0msize\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0msize\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 956\u001b[0m elif (\n\u001b[0;32m 957\u001b[0m \u001b[0mhasattr\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mparam\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;34m\"distribution\"\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
"\u001b[1;32m~\\AppData\\Local\\Continuum\\miniconda3\\envs\\murefi_env\\lib\\site-packages\\pymc3\\model.py\u001b[0m in \u001b[0;36m__call__\u001b[1;34m(self, *args, **kwargs)\u001b[0m\n\u001b[0;32m 103\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 104\u001b[0m \u001b[1;32mdef\u001b[0m \u001b[0m__call__\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m*\u001b[0m\u001b[0margs\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 105\u001b[1;33m \u001b[1;32mreturn\u001b[0m \u001b[0mgetattr\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mobj\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mmethod_name\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;33m*\u001b[0m\u001b[0margs\u001b[0m\u001b[1;33m,\u001b[0m \u001b[1;33m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 106\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 107\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
"\u001b[1;32m~\\AppData\\Local\\Continuum\\miniconda3\\envs\\murefi_env\\lib\\site-packages\\pymc3\\distributions\\multivariate.py\u001b[0m in \u001b[0;36mrandom\u001b[1;34m(self, point, size)\u001b[0m\n\u001b[0;32m 280\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 281\u001b[0m \u001b[1;32mif\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_cov_type\u001b[0m \u001b[1;33m==\u001b[0m \u001b[1;34m\"cov\"\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 282\u001b[1;33m \u001b[0mchol\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mlinalg\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mcholesky\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mparam\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 283\u001b[0m \u001b[1;32melif\u001b[0m \u001b[0mself\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0m_cov_type\u001b[0m \u001b[1;33m==\u001b[0m \u001b[1;34m\"chol\"\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 284\u001b[0m \u001b[0mchol\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mparam\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
"\u001b[1;32m<__array_function__ internals>\u001b[0m in \u001b[0;36mcholesky\u001b[1;34m(*args, **kwargs)\u001b[0m\n",
"\u001b[1;32m~\\AppData\\Local\\Continuum\\miniconda3\\envs\\murefi_env\\lib\\site-packages\\numpy\\linalg\\linalg.py\u001b[0m in \u001b[0;36mcholesky\u001b[1;34m(a)\u001b[0m\n\u001b[0;32m 762\u001b[0m \u001b[0mt\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mresult_t\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0m_commonType\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0ma\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 763\u001b[0m \u001b[0msignature\u001b[0m \u001b[1;33m=\u001b[0m \u001b[1;34m'D->D'\u001b[0m \u001b[1;32mif\u001b[0m \u001b[0misComplexType\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mt\u001b[0m\u001b[1;33m)\u001b[0m \u001b[1;32melse\u001b[0m \u001b[1;34m'd->d'\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m--> 764\u001b[1;33m \u001b[0mr\u001b[0m \u001b[1;33m=\u001b[0m \u001b[0mgufunc\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0ma\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0msignature\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0msignature\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mextobj\u001b[0m\u001b[1;33m=\u001b[0m\u001b[0mextobj\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 765\u001b[0m \u001b[1;32mreturn\u001b[0m \u001b[0mwrap\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mr\u001b[0m\u001b[1;33m.\u001b[0m\u001b[0mastype\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0mresult_t\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mcopy\u001b[0m\u001b[1;33m=\u001b[0m\u001b[1;32mFalse\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 766\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n",
"\u001b[1;32m~\\AppData\\Local\\Continuum\\miniconda3\\envs\\murefi_env\\lib\\site-packages\\numpy\\linalg\\linalg.py\u001b[0m in \u001b[0;36m_raise_linalgerror_nonposdef\u001b[1;34m(err, flag)\u001b[0m\n\u001b[0;32m 89\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 90\u001b[0m \u001b[1;32mdef\u001b[0m \u001b[0m_raise_linalgerror_nonposdef\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0merr\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mflag\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[1;32m---> 91\u001b[1;33m \u001b[1;32mraise\u001b[0m \u001b[0mLinAlgError\u001b[0m\u001b[1;33m(\u001b[0m\u001b[1;34m\"Matrix is not positive definite\"\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0m\u001b[0;32m 92\u001b[0m \u001b[1;33m\u001b[0m\u001b[0m\n\u001b[0;32m 93\u001b[0m \u001b[1;32mdef\u001b[0m \u001b[0m_raise_linalgerror_eigenvalues_nonconvergence\u001b[0m\u001b[1;33m(\u001b[0m\u001b[0merr\u001b[0m\u001b[1;33m,\u001b[0m \u001b[0mflag\u001b[0m\u001b[1;33m)\u001b[0m\u001b[1;33m:\u001b[0m\u001b[1;33m\u001b[0m\u001b[1;33m\u001b[0m\u001b[0m\n",
"\u001b[1;31mLinAlgError\u001b[0m: Matrix is not positive definite"
]
}
],
"source": [
"name = f\"x_{np.random.randint(300)}\"\n",
"x_dense = X_denses[\"D1\"]\n",
"\n",
"with pmodel:\n",
" gp.conditional(name, x_dense)\n",
" \n",
" mapp = pm.sample_posterior_predictive(\n",
" [map_]*500,\n",
" var_names=[name],\n",
" random_seed=None\n",
" )"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "65840ac2",
"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.7.9"
}
},
"nbformat": 4,
"nbformat_minor": 5
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment