Skip to content

Instantly share code, notes, and snippets.

@marufeuille
Last active March 14, 2020 10:49
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 marufeuille/eb5435e81fadc2ce9c0bc5a8912d2757 to your computer and use it in GitHub Desktop.
Save marufeuille/eb5435e81fadc2ce9c0bc5a8912d2757 to your computer and use it in GitHub Desktop.
Display the source blob
Display the rendered blob
Raw
{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Collecting tensorflow-probability\n",
" Downloading tensorflow_probability-0.9.0-py2.py3-none-any.whl (3.2 MB)\n",
"\u001b[K |████████████████████████████████| 3.2 MB 4.0 MB/s eta 0:00:01\n",
"\u001b[?25hRequirement already satisfied, skipping upgrade: numpy>=1.13.3 in /opt/conda/lib/python3.7/site-packages (from tensorflow-probability) (1.18.1)\n",
"Requirement already satisfied, skipping upgrade: six>=1.10.0 in /opt/conda/lib/python3.7/site-packages (from tensorflow-probability) (1.14.0)\n",
"Requirement already satisfied, skipping upgrade: gast>=0.2 in /opt/conda/lib/python3.7/site-packages (from tensorflow-probability) (0.2.2)\n",
"Requirement already satisfied, skipping upgrade: cloudpickle>=1.2.2 in /opt/conda/lib/python3.7/site-packages (from tensorflow-probability) (1.3.0)\n",
"Requirement already satisfied, skipping upgrade: decorator in /opt/conda/lib/python3.7/site-packages (from tensorflow-probability) (4.4.2)\n",
"Installing collected packages: tensorflow-probability\n",
"Successfully installed tensorflow-probability-0.9.0\n"
]
}
],
"source": [
"!pip install --upgrade tensorflow-probability"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"import tensorflow as tf\n",
"import tensorflow_probability as tfp\n",
"import matplotlib.pyplot as plt\n",
"plt.style.use('seaborn-whitegrid')"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [],
"source": [
"tfd = tfp.distributions\n",
"tfb = tfp.bijectors\n",
"psd_kernels = tfp.math.psd_kernels"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"WARNING:tensorflow:From /opt/conda/lib/python3.7/site-packages/tensorflow_core/python/ops/linalg/linear_operator_lower_triangular.py:158: calling LinearOperator.__init__ (from tensorflow.python.ops.linalg.linear_operator) with graph_parents is deprecated and will be removed in a future version.\n",
"Instructions for updating:\n",
"Do not pass `graph_parents`. They will no longer be used.\n",
"Step 0: NLL = 692.2976799192386\n",
"Step 100: NLL = -17.5254993250389\n",
"Step 200: NLL = -25.008495231716566\n",
"Step 300: NLL = -28.797935669302518\n",
"Step 400: NLL = -29.543987900853736\n",
"Step 500: NLL = -29.560702537329917\n",
"Step 600: NLL = -29.560715783945888\n",
"Step 700: NLL = -29.56071578399623\n",
"Step 800: NLL = -29.5607157839962\n",
"Step 900: NLL = -29.560715783996045\n",
"Final NLL = -29.560715783996102\n"
]
},
{
"data": {
"text/plain": [
"[<matplotlib.lines.Line2D at 0x7f08c36d5310>,\n",
" <matplotlib.lines.Line2D at 0x7f08c03e2510>,\n",
" <matplotlib.lines.Line2D at 0x7f08c03ee950>,\n",
" <matplotlib.lines.Line2D at 0x7f08c03eeb10>,\n",
" <matplotlib.lines.Line2D at 0x7f08c03eecd0>,\n",
" <matplotlib.lines.Line2D at 0x7f08c03eee50>,\n",
" <matplotlib.lines.Line2D at 0x7f08c03eefd0>,\n",
" <matplotlib.lines.Line2D at 0x7f08c040a290>,\n",
" <matplotlib.lines.Line2D at 0x7f08c040a450>,\n",
" <matplotlib.lines.Line2D at 0x7f08c03eee90>]"
]
},
"execution_count": 3,
"metadata": {},
"output_type": "execute_result"
},
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD1CAYAAABA+A6aAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOy9eZhcd3nn+z1LLaf2rau6em+ppVZrs7VZWLKRdwWBbeFAbLiETO4YknCTAGO44GTuxTN48BBDJhkuM9mehInJgDA4gtiJZWMLL1i2ZEnW2lKr97W69n2vOvePt49Oq9WSWt1VXUufz/PU00ttvzp1zvf3/t7fuzCiKIpQUFBQUKh52EoPQEFBQUGhNCiCrqCgoFAnKIKuoKCgUCcogq6goKBQJyiCrqCgoFAnKIKuoKCgUCfwlXzz48ePV/LtFRQUFGqSbdu2zfv/igo6cO2BLYTe3l709PSUcDSlQRnXwqnGMQHKuG6GahwTUL/jup4hrLhcFBQUFOoERdAVFBQU6gRF0BUUFBTqBEXQFRQUFOoERdAVFBQU6gRF0BUUFBTqBEXQFRQUFOoERdAVFMpNJgMUi1f//9IlwONZ/vEo1C2KoCsolBNRBHp7gYGBK/+fTgPRKBCJVGZcCnXJkgS9r68P9913H374wx9edd8777yDT3ziE3j00Ufx/e9/fylvo6BQu2SzQKFA4u31yv8Ph+lnKlWZcSnUJYsW9GQyiW9+85u4/fbb573/6aefxve+9z386Ec/wltvvYX+/v5FD1JBoWZJpUjQNRpgfFwW8PFxstqHh4FksqJDVKgfFi3oarUaf/u3fwun03nVfWNjYzCbzXC73WBZFnv27MGRI0eWNFAFhZokkSDhFgSAYYA33wR++UvgvfcAnqf7T55UXC8KJWHRxbl4ngfPz/90n88Hm812+W+Hw4GxsbF5H9vb27vYISCdTi/p+eVCGdfCqcYxAaUbl+bsWejOnEEmlUJRp4N6cBBMLgcmk0HaaAQfi6EwOAhxagq5xkYU7PZlGVcpqcYxAStzXGWptiiK4lX/Yxhm3scupepYvVZTKxfVOK5qHBOwxHGJInDxIqDTkYvFagWmp4G2NqCnB1CpAIcDCIWAWAy47TaKgmEYYN268o2rTFTjmID6Hdf1qi2WRdBdLhf8fv/lv6enp9HQ0FCOt1JQqD5yOXKlTE0Bx48DLEu3qSlg7VrA5wN27gT0eop2OXmShFwU6XYN40dB4UaUJWyxpaUF8Xgc4+PjyOfzOHz4MHbv3l2Ot1JQqD5yOfrp9QLBIAm4Wk2bn0NDwOgoWe5tbYDbTb8HAhSTnslUduwKNc2iLfSzZ8/i29/+NiYmJsDzPA4dOoR77rkHLS0tuP/++/HUU0/hiSeeAADs27cPnZ2dJRu0gkJVk82SK8XrBbRacr243ST0DEPW+tmzFPmSyQAmE2A0AmNjJO5abaU/gUKNsmhB37hxI5577rlr3r9jxw4cOHBgsS+voFC7ZLNAfz/FnrMs0NaGo0kVXhvIo2m4D7xOwK5gGp0qFc5MJ/DyUByvW1dhR9aHXYZ2/Mbe7ZX+BAo1SsVb0Cko1B1+P22CJpOAWo1Tnjj+P08eUX0LtlliiKu0mIyksOGiH/0DUzBkM3BobPDlCvjLgyeRdrqxf0tzpT+FQg2ipP4rKJSa06eBeJwSikwm/GIogUv6BuR5Fd5t24RhezNONnTi38IqhHgdhFwGa71D0GTTUCViePbQxUp/AoUaRRF0BYVS4vXSxifP002nw0SegyGbRp7lcL5xNSYsjWDBwKszI6XRQp9JYtP0INb5hiHk0pgOxiv9KRRqFMXloiAjiuT/1WgqPZKa5OdHh/GrvzqANWePI2t34N5VFmyyWqE3A5piDmNmF0SGhU9ngSMegjmbREytgzGbAFcowsKrYE9E0KlX7CyFxaGcOQoUaTExAZw5Q9EXShr6TXPw+Bhe/s7fwzrUhwyvgjfP4qULAZyajOE3t7agoNMjIhgBABmVBkmTFXf1uJA02dBva4HIMtDmsuiMTePLd7RW+NMo1CqKoK9kolGqyX32LMVA63QULz0xQda6woL5nz8/AUfQC3siDEYUkVKpwWYzePP8JHZtXYU/+J170GwRwABotgj4o0/fgd99cBs+u+9W9K/fgbRKAzNbwMNuHvucymWpsDgUl8tKJJUiP28qRWnoTU2Uiq5SUTr64CAlxNygroiCTNLrhyseQJFhwIoitNks7MkIBrUGYONGPGS14qGdq6580tgY7m5pwd0P3QF8bQzo6wPyUcoc3bCBQh4VFG4C5YxZaRSLJNj5PNDRAWzaREkvKhXdb7WSpT45qVjpN8GGYgzGTALjZheKLIfGRAAFhkGocw0d0/lobQVcLkoquvVW2rtgWYqSuahEuijcPIqgrzTGxqh+SGcnWeDz1Q1pbqbNUZ9v+cdXo/z7VgYqhkGaV2PaYEWWUyNkcWL/Q/P3C7iK226jDFG1msIdFUFXWASKoK8g2EiEkl4aG8kqvBZSKvrUFImLwvURRexg49i5yg6nwCKj0iLjasTd996KB3Z2Lew11q0jaz0SoZXR4KBS10XhplF86CuFbBaqyUkq39rUdOPHNzcDFy6Qld7YWP7x1TKpFDA9jdUtNvzp1lvIrZVOA+vXAwbDwl5DqyVRHx+n1ZHXS6uprgVOCAoKUCz0lcP0NBhRBFatWlh5Vr2eLPXp6fk71ivIjI9TMS6tlo5ZJEJ+c5Pp5gpt9fTInY0CAYo2Uo69wk2gCPpKQBSBYBAFo/HmkoYaG2nzdFZte4V5GBykui1Op9wftLGRJsWbqW3e1UXZpYKAIV8c//l7L2Hbl36M3f/1dRw8OVGesSvUFYrLZSUQiQD5PAoWy809z2gkl8H0NNDQoDReuBaDg/TTZqPoIIaRVzg3g8kEOBy4OBXBexMJGNVjsJk7cClsxpMvnAEAdCuVdRWug2KhrwT8fkClQnGh/tzZNDaSTzcQKP246oFkkvYZGIaO09AQbW62tJDFfjMIAuB2o7dvEh6NEdZUFLpMAsZMAqlcQSnapXBDFEGvd3I5ygi9VojijTCbKS7d41Hi0ufD4yFB1+tpE9NsBu68kybCmz3eWi3Q3IxsOoOYSgBXKKA55oMtSaUYJsOpMnwAhXpCEfR6JxgkIV5K1mdjI4XQxWKlG1e94PHQhFko0HHu7l58VBDLAm43GL0eYFioiwW0Bj3Q5ih8sckilHDgCvWIIuj1jt9P1uNS2ppJvmBpw0+ByOfJKo/PlLs1m0nMpazbxeB0YlNPC6DRIKLRoTM8iYZ4EALP4qt7u0szboW6RRH0eiaRoHhoh2Npr8NxJFLpdGnGVS/EYsD779MqSK+n42y1Lq0Gi8OBdasa8fDOTiRcTdAU8rgjNIhvP9itdDFSuCFKlEs9Ew6TH/datURuBkGgBBoFmWCQYtATCRJzs/n6GbgLYSY6ZrvJhO2dHwLeY2gfpMtcmjEr1DWKhV7PJJMkxBy39NfSahULfS4eD4V0CjO+bZOJarEsBUEgt00qRe4yhiEf/fDwkoerUP8ogl7PJJMUoVIKBIGyFpX6IjJDQ3SM3W45K3Spgq5WA2vXUuijyUSbrcEg0N9fmjEr1DWKoNcr2Sxt2pVK0KVNVcVKlzl3jgTX4aC6LYKwdEEH5IQup5Pi2fN5KgOQSCz9tVcQnN9PXbhWULitIuj1ihSRUkoLHVAEXSIcpggXnidLev16YM2apW9AA0B7O2Xm+nwk7DxPvnolueumUE1MACMjK+qcVQS9XkkkyP8qlCh2WYp0UTZGyeK7dIlcISoVuUekYlyl2K9QqahjkcVC76HVkj89GFSKdS2URAK8z0eT4uQkraTC4bo/fkqUS72STJIQlLKNmSCsKGvnWvzba6fwzt/9Mz56agg5vRE8TLh9qdEtc9FogJ07ya1TLFKIZCIBthQTxkogGCSXSyxGk68gUH1/nqdVlNO5tHyBKkWx0OuVZJJC4EqJVrviLfSD74/ixz94GYbxEXBiEaMaM/7zhTwOfjBZ+jdrbibrX6UiP7rXCy4cLv371BuiSBFIxSLtR4yNAaOj8orK45H3P+Y+p7eXwkRrFEXQ65FSb4hKaLV0kWSzpX3dGuKvfnEClogfqmIRIhiMm12Y0hjLUziLZSmEUaej6KJAAKyyQrox0SgQiYBNpeh8DQSADz4gYySVouMai9EtGARee402TycmyBCKRiv9CRaNIuj1SKk3RCUkf/wKttLjviAAQJdJgQGDSbMTGU5dvsJZzc200ioUgIkJMLmc0hbwRgSDgNcLNpejCTGdprBPq5W6dTEMWewnTgAvvgi88w5Z8atXk0umhgVd8aHXI8lkaTdEJWaHLppXZuZih4Y21RrjfiRUGgT0VmR5VfkKZ7W1ydmnwSCYdJqs9VJP1vXCzObn6eN9eLUvjV9GvfjUwHnsdqiwmufxi9EU/v7gBTReOIdbEr/EvS411jZbaTNbyvSt4SJ0ioVejyQSpd8QBch6WeGRLp/b0gA1irCloohqdAjqzVBr1OUrnGW1UuMMlgXCYbCxmLIxfT3CYbx9pBevnRjCOZ0Tw9YmpFI5fHDJgwMHDuO1P/s7qEeGIWSTyMXiODDN4HRGRdEwoRBFKuVyNXuMF22hf+tb38KpU6fAMAz+5E/+BJs3b7583/79+2Gctev/ne98By6Xa2kjVVg4yWT5LOgVXgJgT7Mehm4zVB/k0Gt0wGw14ZlHNpWvcJZaTW4CngfSabCh0Io+/vOSSJCR4XAAHg9eP3wKHDgM2VoAAHGNgMaYH4eHIvDrLMizHJqifhQYDhctTfiJJ4fNqRT50NeupdeUesTWGIsS9KNHj2JkZAQHDhxAf38/nnzySTz//PNXPOa5554ryQAVbhJpQ7TUES4SgrCyE1zicWzjU4DbiM2/fT9+72N7gM4yV0Fsa6MwxmQSnFRBU0FmcpL83sUiMDSEaDwNUTAhoDPDmo4ix3BI8xqkWR6nG7vQEvUipDWBF/PgxSLOcyY6vn4/JXWp1fR6DQ2V/mQ3zaLW5EeOHMF9990HAOjq6kI0GkVcqgkNIKGkKFeOcm2ISmi15KdciZEuuRz5r8fGyPW0evXyWHE2G624Uqmaj8IoOaIo16M/cgQYGoKbKyCt0kDIpbFpsg9FlsWoxQVDPoue6QFYUzGcc3ehwHBQ57PgXU5aAWUyVGzNZCILvQZLBizKQvf7/diwYcPlv+12O3w+HwwzPSvD4TCeeOIJTExMYOfOnfjSl74E5hrtuHp7exczBABAOp1e0vPLRSXHxU9Pgw8EkJ7Hh16KcTGJBDTDw8gWCiiWIJmmlr5DNpGAYWAAhpERsAD8iQSyw8Moljk2XOX3wwxAE40CgQCGT51CSqMp63veDJX8DplUCprBQeQNBlh/8QsUtVrc1sgjM+SHOHEB7aEJxNUChpztuKVFD2PvKAqFIsbMTjTG/GiPTqOrqYB3/zWIw2EV/NDgwrqt+J22IjYXixBLHViA8h6vRQm6OGfmEkXxCsH+8pe/jIceeggajQZf+MIX8Morr2Dv3r3zvlZPT89ihgCAJoOlPL9cVHRcUuzyunVX3VWScRUKFBFgt5MrYInU1Hc4NUW1VYpFYO1aGFetAnp6yh9x4nYDx48DExPIh8NwO51AV1fVZDpW9DucnqY6Nx98QBEqH/oQnOEwLLpRvH7Rh6woQmWz4dMP7cTdXBQXA4M4MeCHMxFESyGB9T0dCK7pxts/fwvWuA+CSoM+XwsORNQwbmDxkftK/7mWeryOHz9+zfsWJegulwt+v//y316vF45ZRYk+/elPX/79rrvuwsWLF68p6AolRBRpg+hmu83fDBwn1xhpaSl9JE01E4uRgORyQEcH/W85LGWjkdwuajW4ZJJqkqTTVSPoFSUWA44eBQYGyMCIx4GTJ7Fp+3bod23GKp2O9nxyXmBoHN1b1qH7Y83A+fNUuOvWbnx8Sod2wYL1sRBsyRhaI15ccrTjBy+fxkfuu7XSn/CmWNTVuHv3bhw6dAgAcP78eTidzsvulmAwiM997nPIzaTPHjt2DGvWrCnRcBWuSyJBoj7zXZQNh0MudrSSiMVoA47jSDx4vjTFuG4Ex9EkLQhgZkoAKHXpQee610sWutNJ38n779O52dGBzMaN5A8fGSFRv+UW4I/+CLj/fopm6egAIhHYLp1DVGeEx2iDqpBDa3ASaV6NeCBcc370RVnoW7duxYYNG/DYY4+BYRh84xvfwAsvvACj0Yj7778fO3fuxKOPPgq1Wo3169cr1vlyIW0OlVvQjUaKBAgEyHJcKfh8VO9Dp6MMzuUMa7PbAb0e/qEJvPDP7+KvL1rAtrXiq3tXcK/RdBro66Pz/s47qe5NNEouR5sNqrExEny1mjawH3yQEu6MRhL0hgbg7bfRkwvjPX0TRIZDnmHREvODK+bhMqjJSCr39VRCFh2H/pWvfOWKv9fN8tk+/vjjePzxxxc/KoXFEYtRWCG/DAnAdjv5lLPZ0jR1qHYKBRLzYJD811br8rhbJCwW9GdYjIYyYMQIhGwKg+EUnnzhDACsTFGPxYBTpyhMN5UiA+ORR0i8338fmvPnafLduFEWe72eLHqGATZvBjwePBDpw/hQGCKADK+CKx6GOxPDb9+1kaKKakjQV5ADtM6R/OfLdfLZ7fQzGFye96s0mQz19cxmgc5OcoMsp4VuNuOwv4gcw8GcisGaotDFVK5QnsJgtUB/P30nOh25/1iWJlueB0QRKr+fBHzrVjJ0xseBs2dJ+BsaSORXr8bm1U483sZCMOqR4VVwFhP4uj2Guze4a65LlFLLpV5IpciKXC5B12ho6er30xK33onHgcFBEgspQmE5LXSbDcOsAS5eA0s6joZ4kCZxhilfYbBqJpcD3n6b3C4OB/nKOzpoFaVWA9PTKJhMdG42NZHYT01RPH9zs1znaPNmgGWxcWMKG+Nxssifew6YnpksOjsr+SlvGkXQ64Xl8p/Pxm6nkz4er6ll6aLwemlD1GCg6B5geS10gwG83YbItAENiRBcsQDUhXx5C4NVM/39FKlSKNAk29VFE2x/P1ngoojErl3Atm0k5m43na9zJ2GDAfjQh2ileeIEuWK6umjy9nrpMRs3Ls/mdwlQXC71QixGJ+ty+rMtFroAIpHle89K0dtLE5fUbg5YXgudYbB391pE9GawKMId9UFTyEJQceUrDFbNvPMO4PNhrMjhvw8X8PCPL+J3//EE3jvWR6vVnTtRcDplQ4Nhrv992WxkzScS8oozHqframio/J+nRCiCXi9UwkrmOPJfzir7ULf09pI16HaTMKhUyx6Dv2vPrdiy1gGe49AS86JNx5a3MFg1ksuRwJ44gcloBq8HOfRzRkwbbHhD14Q/8ljwr+6N8irqZli9mowUgJLHJPfNwEDNlNRVBL0eSKdpp78Sbg+Dgayaem6+GwzShppKRUIhipWpxLdxI9pu6cK2FjM+5mDw0qd7VpaY+3zUWejMGSCTwdlIDlG1Fn69FQmNgCLLwcvr8F+OLrJ4nCBQ+KPkc+/ro/M6GgUOHaLyulWOIuj1QCX85xIGAwmcVBSsHpmepotZraaNt0xmed0tEiyL5Pbt5Av2+0lwVhJTU7QinCm1EMszCAtGxDU6JFTyBLukTWKdjkS9o4Os8ldfpffNZIB33636RCNF0GudfF5eGlbCapQmkXp2u4yOkqBLCUX5fMVqZRcbGmgTL5sF3nyzImOoCIkEuVtsNqpbnslAq+YwbnQhqSbrXGLJm8RtbcCtt9LEOTICXLpEmaiBAG3EVjFKlEstI4rkT8xmge4KbYzxPC1VY7H6DV8cH6fP53TS8tvhkNvCVYIHHwT+8R8pKiOVKn2rwWokHKaNzUKBIlD8fqxZ146E2o6kWp5cS7JJzHFU3K63F2PBJN44MYnTl/4FVoHH3guT2PpwkiJfqvC4166Fns2CqbGg/5IjFfZva1taQ4t0emk+cMmPXuXL0UUzOEguJYuFYtA3b65sT8/mZpo8fT5yB6wEwmE6zy5cAM6dAwCsfeg+/PsH1sPUYAMDoNkilG6TuK0N7zrX4FBMjXiugLhGh8GiBmcOvY3+P/sehetWIbVroQcC0IyM0NJoOVLdq41wmFwtDgfdFkMiQc0apIlRrSahstvl3f6FYDCQuKRS9de8OJkk0WQYqv/R2lrpEdHk3d5OrqCTJ4FVqyo9ovKSTtMtlwNef52+k64uoKsLd3Mc7n50c+krT+r1+LvhHNyGBmzy9qMoAilBB6/GiJHeYXS99pqcYFZF1K6FrtORhT6rjO+KYmqK/LiLqUleKJBv8MIFcte0ttLOvl5P4j4wQCnSPt/CXq+e/eheL0W58DyJaDWg1QK33Ua/v/VWfUcYAUA4jMMXvPjGMwfwP18+h1dH4ziltcv1dMpRRthsxkiGxZTJDrYoYnVoHK3haUwb7RhXGan+ehVWvKxdQU+lwIVCZKXU61L/WmQyZKU4HGQ5zsbrpYJFY2Pz954sFuVsOpcL2LCBfMNuN1l6mzbRT56nY7uQ/qFqNd3qUdCnpihxShAWF9tcLnbtIiGTEp7qmF8e6cOfvzkCzucDX8ghXmTxZx4DXj8xXL69DJZFvr0TQb0VMY2AlvA07PEgimDgae2iSf7YsfK89xKoSUE/eHICjz/1PH524B386V++hJferO6d55IjFcSSMhZnI1nVPh/5Gi9dkl0qokjWdzxONSpaWq5OaWYYet1160jEpPTnG2E01kzyxU0xOSlvPFbTpm9jI03ok5MUVlmv5HL436+fRyaThSGTgCUdR5LX4LTFjR++PVDWUN0vfmwTJho7MGRthi0VgSsehI4p4I49G+kBv/pV1TXsrjnn88GTE3jyhTPoiafAikXk/X788B8OIWeyrJwki2BQrkk+m1SKTrC2NvKB+/0kyBcuABYL1Yd2OMh1MN9kMBenk1wzsdiNLSGDgaz5dLpiIX0lJ5uVY5Db2xe/V1EOTCbKbHznHbIU67WJTCQCXywDczYFczoOVTGPPlsHVAwDXzRTVkEnPdmFn6f8GIh5odUb8OCdPdjZ04z+oyb0/fRX+MnJLDwbt+H392+vCv2pOQv92UMXkcoVMGhrQUalBVcsoiEwhf9x8Np99uqKZJJEcz5BljLZrFZajrvdFF7V1ATEYuBiMbLKFypMNhu5XhZipdejHz0WI0HP5WifoZomKr0eWL+eMhrfe6/qLMWSEYnAYtGjOeqFOR1HjuUwZXDAFQ/AYtGXPcFr/5Zm/MP3fg9fePbL+D8/9WHsbDbi2PF+/DiggjoaQo93CBFfCE++cAYHT06UdSwLoeYEXcoCUxVy0GQzYIsFiAAs/SvE7RIMym6R+e4zGGjTU4LjLgt7tqOD/OYLhWWpbnQ4fOMNIK2WVgx1VKjr1TfP4eC/vY9z0zE8cymHg5eq6LMxDEV4abWUizAyUukRlYdkEp/98BqsD45ByCbBi0Ce51EU9Pit37precbAssDu3XQdpdM4dHIU2SKDPMuiPTIJV8xfNXXpa07QpSwwZzwIezKMjvAUYmodegpRyuCrd0IhWm7PDdVMJsmSnp6mCJVTp8hfLm1q8jyKi4lVb2gg8ViIlW42k1VbB5vUrw/G8Pc/PwYhHESBYXFU58aTB89VhRV2mbVryRUWCAAXKy8mJadYBDIZ3K+OYa8+C41ahQlTAya7b8EfPL4XD+7qWr6x8DywfTtgMEDl98KSjsOrt0OTzaE9OAmIYlXUpa85Qf/q3m4IKg5ZTo2kWouGRBir4tN4aI21vjeHAHJnZLNX9/HMZKjz+cQELb3jcRL+wUHaGM3lgFwO7GI2LVUqej+//0rLfz5MJnpMHbhd/umoF+p4FOZ0HHG1DkGjrWqssMs4nbTiSibJSq+38MV0ms7ft99Ge2wa9922Bn/4hw/hF8/8VmX81Y2NQEcHHCpAn0ti2OYGAxHNER80+WxV1KWvOUHfv6UZzzyyCWJ7BybMLuh54JEOHbYJOUrRrmeCQVr+zU36mZigEEOep6Sg5ma62HU6+v9MzLl6dHRxLbWcThKLG7WbM5nqpj56PJKAPRmBkE/DrzMjJJgALLHwU6nR62lPpFC4MkGsXkinyUgbHCTXYXc3ieoyly2+DMMAq1Zh+/a1YDkeEa0JeZZHYzwAG/JVUZe+5qJcABL1/VuaMfDTFFb/LA5waYBh8P6//RpfPJzAZDiFJotQfx3RpU1PlqUsUa+XltzHj5PY7thBJ720QZnPk+V2+DCwYwdEhqFwRr2e7ltohq1OR5tP4TC5YK4Fy9J4IpHqitleBM3qIlrDHrCiiKDOgqiWjmk1WGGXYVnKI3jjDdq8DQQqW2Om1EQidI77fBRme++9lc/OdDhwywO3g4+E4fUWkVJr4Sqm8J92OfFAFWhNTQq6RG7NGhKwkydxirfgjTN9SG1qhKgzY6IeO6IPDZGwZjKyb3xsjPynnZ2UPShlzeVyJMAcRxeGxYKCxULPSyRoQrjjjoVHCVhnXFo3mgjMZhpTpUrMloj/o5MD/0oIIhiMWBsBhqnO7kAbNtBmtM9HbrGOjkqPqHR4PGSdF4v0OW22yrc65HngQx/ChnPnsGEjC/QItGfFVket9JpzuVxBoUAilsmg793TYDIprPUNX96Uqzqf51JIp0mI1WpKFkqnyX/K87T7vncviXkiQQJ/+jS5W5qa6CL3+1HU6cj9cuoUWXQ346KyWum43sidYjbTzxp3u9yhS2G3nQOn5jFoay1t4adSsmYNiVwotPBSDbVAKkXhmF4vrSi7u6tn1afX0zWXz9P5XizSNVcFwQC1a6F7PNAMDpILwG6Hs/cseMEEkWUxbmrAuLUJumwKzMgCMx2rHb+fTpiGBrnhAsOQKK9dS7HlUjndYpGE3GKhk+7cOSCZhFq64KXuRlNTlJyyEHQ6mkxCIfLTXwuNhkLpIhHyvdcoXCSCViYLtNrxv/7LpyhzthppaKCJ/cyZ+gldnJqifaH+fvq7uZkmrmpa8bW30/HW68moGhujhiOtrRUtUFe7FnpDA7KtrcDWrcBHP4qozQUGIlrDHkWeibwAACAASURBVOztexdMsQBnPIi1fIYiQ2odycXS3U1iOTBAVkEwSMI9OEhCn8mQRe52U7q6Xk9CbDZDlJ7PsmTNR6M3t5FmtdJzbhTtIoUv1mrURS4HJhajz2CxXH8CqzQcR2KXz5PAVGHBqJvG5yPxzuXo87lc1VdRsrGRxlgokIAPDdHKeWCgouHTtSvoHIeiyUSbQBs2oGtrN/JqLZhCEe3hKXQEJ2AvZvDbH2qvjyy6QICs6nSahMZqpRPH56NN0SNHgPffJzE1meTnsSyJukaDTE8PHS+epwslk7lx5MpsLJaFu11EkcS/FkmnwcXjdHxcruqyDOdj0yb6ngcHyVVRy8yE2CKVonNerSardyn1/suB1UoTPcPQsY9EaMxSE+sKuV9qV9Bn09mJnp2bcGt3EwQNh6aoD/cG+/GHd3fh7nXO2rdaRJFcHXo9uU/0etrQlMoHZzJ0Ik1MXF1KVBTp1t8P9eAgCXhrK/1PEvSFioDBQK8fDt/4cSxbu/Ho6TTYUIhWdk5n9Qv6li1yiGqt+9Glc3FsjAwCk4k2/KsNnY7cXSqVXG3U46E6StFoxRqP1K4PfTYmE7B7N7oiEXR1NQOnTuE+RwzY1ikXrKplEgn6DCoVLfE0GvKd+/3kcurupuW2IND/pJojySRZ86EQEI1CFYvRibdlCy0Px8cpIen0aWDfvoW5FqRImWLx2vHADENjqFVrMZUC7/XSSqa19eoiaNWG200ridFRsg5ruVBXMknn+PAwTahW6+Jq/pcbjiNB9/vpeI+P07E3meg6mpqi1fAyh5HWh4UO0And2EgHWhI8qfJfrVvofj9Zu6kU7fpnMiTU7e3Aww/LiUTt7ST64+N0i8fp83d1Aa2t1BAkHKaLZNMm2jg9c4Ya3544sbCxWK0k5jdyuwhC7Qp6PE619lmWrMO5NeerDaORRC+bpcm5lkml6PyemKBVpM1WHV2i5sNkImNq7VrSnOlp2UrnuJtzZ5aI+hF0o5EsFa2WDm46TVao9Hst4/eT5VIskrtlwwY5XNFkohNeanhgs5HFcMst9HssRifZ9DS4RIJONKmDPc/T8ZIiYRYSxmgwkMV6o05RgkBuoFqsrxMKgQuH6djUirW7fj19txcu1Pb5LlUT9Xjo2nW7q89/LqHX06TT3Ewr12CQInOkfasKZO7Wj6BrNLQE0utJ5ASBwoiSSbJcqiBGdFGIIrlTJDHX68l6mZ6mSUxye2g09LnPn6fNMZ+PLgqbjaJeursharUk4h4Pve74OFnpt9xCjRJefJGWutc7VgxDS8po9PrRQ1JH9Fqz0otFIBwGn0zK5RNqAclKnJykpX8tUizKdYhSKTIcqiX2fD6kc8NgIA9BNksBClJYcCp144iwErNoQf/Wt76FRx99FI899hhOz1nmvfPOO/jEJz6BRx99FN///veXPMgFYzaTS0AQ6KfPJ+8416rbJZGgi9RopJPD7ycxVatJnGdbylotuWQ8HvKNG40k5nY70NVFYZ5mM1000mTX0gLccw9ZePE4nZBSzPu1kOqpX89Kr1VBz2RoApV6rVb7hqjEqlV0zkciNKnXIlIEl99P543FQgZHtSIZWMkkne+CQC6vgwfJWs9mlz0wYFGCfvToUYyMjODAgQN4+umn8c1vfvOK+59++ml873vfw49+9CO89dZb6JcSBMqN0UiJMqtWkSUpinKLrloV9IsXabMrFqP+kQaDvBpRq8nKllYgkYi8IZpOk3Uu+X8FAen160nQpegYq5UeY7XSMWtqotc7doy6yff2zm+Fq9W0GggEri38KhVNOLUm6Ok0fW6A/KO1IuitrfSd5HLkaqzFFWkyScaK10vWejX7zwFaHXd30xjdbjku/YMP6PoKh5fd7bIoQT9y5Ajuu+8+AEBXVxei0SjiMzPR2NgYzGYz3G43WJbFnj17cOTIkdKN+HoYjbQM6umhg1sskqh88AGFQdUaokjpz8Uizf7NzcC2bSSSDgctswFyk3g8NGlt3EjLvzVraKd99pJPpaLnSQ2dW1pky9/tJvGWarA0NNAFdq3Nz4YGevz1Ys1rMdIlnQaGh6mQ2bp1la8dslB0Ovrei0USxBoMGX3lyCV880fv4q9f+gDvDgRxRjTIpSSqFYahYIy776afJhNdQ8PDsgtpGVlU2KLf78eGDRsu/2232+Hz+WAwGODz+WCbVa/b4XBg7Dpi2itZQ4sgnU5f9XzNxASYRALGTAaqS5cQczjAx+MQf/ITRMNhFJch62++cS0G1u+H7exZ5CwWwOtFpqcHhVOnIAwNIdHYiMLgILhkEqqZFVDBZEKhWAQbCqFgt0Nz8SLyPh/ybrc8rlwOWp8PEEUUGAZcMIhcJgNRFGEcHgabSCDb24uEywXN+DiKgQBy81lJogjN6CiKfj9y1wgr4ycnwUUiyFzHWizVsSoVqqEhNPX3I6/RoE+vR+HSpUoP6Qqud7xsggBDMolMfz8Cx45d/t4rOaaF8vpgDG/+9NfojMZhSscRFhn881AKHYfP4K7uxV2zy3puZTKwGAxQjYwA8TjyuRwKdjsy8Tgy+fwVkVLlHNeiBF2cc4GKoghmZsBz7wNw+b756FlCOcze3t6rn6/TkdugqwsIBGBOpYA77wTGxuD2+6mYV5njiucd12J4/XXyfzc20ue67z6y2Ht6KKpF2hAdGCBLYP162jNgWbLQGxpoH8HtBiwWeVxWK1lxa9fSBqrFQlb62rVycpKUnReLkRtreJjCsWb31TSbyZ3V1XV1QhNAq4HRUXr+NY55yY5VqTh3DsjnEWprw9pdu6ouZPG6x2vbNuDoUehZFjZBWLZSs6X4Dh//+evYEA3DGQtBm8sgpDfDqzHi7d4M/mD/4l572c+tyUlq2p3JyFmj8TitYtetu9zHYKnjOn782v2TF+Vycblc8M/aEPN6vXDMbJTNvW96ehoN16uhXWqamkiYdu2iuGyvl27t7XRwayWTLh6Xm1bk8yTKDEObLe3tNGlJ1uOqVeRqAUiAAfrMLS0kykND5D6RaGqi+GrJRRUOU2RBVxdNIOEwuWsMBjop+/vpdefG1drt5BaS6szMpdY2RjMZ2j8AkO3qqjoxvx4HT07g998O4M3JJN47PYKjvzpZ6SHdFBFPAA3xAMzpOAosj7BghF9vwWSkhkIwN2+ma0YQ6GexSIZMKrVsurMoQd+9ezcOHToEADh//jycTicMM77GlpYWxONxjI+PI5/P4/Dhw9i9e3fpRnwjVCqyTGcn2pw7R77RhfbGrAamp0koGUbOWBwaor81GjpBolESWoahGGQpJd9iId93LkfWMc+TKEsbnBwnt7ETBJow1GqyIhoa6EQcGiLBl9KYWfbqlH+tlh5zrZO11gQ9GKQIF60W+eYqK5N7HQ6enMCTL5zBB7wVKZUWYiaDX73yPg6eqJ0OXpuYOJqiPgAiEmoNIlojolpDdTUUuRE2GxlKUtBANkt6w3F0nS5DCOOiXC5bt27Fhg0b8Nhjj4FhGHzjG9/ACy+8AKPRiPvvvx9PPfUUnnjiCQDAvn370FmJWgw2G914Xq4PvnFjbQh6NkuuD6+X3CPt7XSSXLhAIhqLkTtDqoEtpReHQnR/WxsJus9HVnpXF3DxItTj42RFzLY8paSN1lZ6HSmOf3qarHopGsbtpjFls1e6TxoayG0j1d2YDcfJFkotMDVFoq7VIlfN8c9zePbQRaRyBeR1ZsQ0OlhSUVjiQfzli2ewf2ttfI4vuTMYyKUQ5rXIiSoEdWYwglB9DUWuhyAAra344P2LOHJqENlYHHphCpszGtx27w66Juf2Ay4xi67l8pWvfOWKv9fNqhe9Y8cOHDhwYPGjKgUmE92MRrLOz5yRSwOkUrL1WI34fGRRRyLUadztlmugNzWR71pKL/Z6yRJnGBJ6t5tWJVYriXFT0+UTje3vJ1fO7PoSFgv5WqUkCbtddrWcOXO59C4sFhL0SOTKNnQWC02aPt/Vgg7UVgmA4WE6PmYzCjVUy13qc5rj1QgLJrSFpmFIJ5H01IDxAgCFAnaycdhX2/Hr6QzimSIKzkb83w/fin3V1lDkejAMXvGL+NeBJExFHk35HFLJIn727iDyBiN2zUSOMWUMZayP4lzzodeTmEQiwJ495G8+eZJ865FIdQt6IEBtrYxG2ugaHSULWEozluLspWYXfr/cFs5qpZ9OJ1mbgQA9zmqlULxw+EpBZ5grsyEbGuh/gkAxtFJlRqlGzlxBZxhaLUxP0yQwd3NUEOSxV7tPenKS/OhNTcjXUG/OJouAiXAKYBj49WbkOQ6GTBLruBqZSMNhwONBl8uErhYVGSUf/jCwsQqrLN6AvzzqgYMTkDDaYUin0Bz1IhoO4rU3z2JXhxVwOMCUscl1/aT+z4XnybLM50mIbruNTpzR0RvXIakkxSLFzUciwM6dNP5AgP52OMj9Mjfd3+eT3S3SRKXXk1B7vSSmLIuiwXDj0rcmE72uTkfWtxTrHovR8YxGr25cIXVLmu+4CgLdVwv1RaSOP2YzijMRCbXAV/d2Q1BxAIAJkxN5hoOxmMHvt9eIvTY5See4VKe/pYUMg9kRVTXCYAqIawRkeA0G7G4I+RR6fIMQPJNUAG9oCEwZG7/Ur6ADcgu0WIyqC5rNZKlLldyqkVSK0vL1eqqxMrsuh8t1daJFQwNZxlLTi9m43SSkM7WZi9JGzfWWfNLKJpcjYddo5NrmUuMKKZJGQppY5isZIE0w1S7o0SitaKSVT7WXzJ3F/i3NeOaRTWi2CLhkbwUjqNFjVeF2pkb6ug4O0n6NWk0i7naTsF+vGXmVYmuwIqkWkFTrILI8RJEFI4poYPO0MrZaUSzj56p/QRcEvHesD5/78Wl8Y4jFq6dGcOalN6o3k254mKzq5mYSVcmdwnF0QsxN0jKbZfGZK+hS+7SpKSAWQ8FgIMG6npXOMLRxk8mQJR6LkU9d8r3PjqaZjcMxf+aoVkuvWe1+9GCQPhfP08RZY+zf0oxff/0evPzfPos7d6xFIy/Syu1GK7JqoL+fDIFikc4jrbYmrXOAVksFQY+oVo+4WkBCrYMWBdzPhOi6jkZRTsdjfQu61Yp3/Hm8+GYvwqEYRqxuTHB6XHjjGF7713crPbqryeWA8+cxNjCBP39zFP/uGwfwJy/14b0+D4l1Lnd1LDjD0ManVJRsLlKBqeFh+nshbheHg143n6cJxWiUK8eZTOT+yWZpcjl1ikTcYiGxn1sqQAqzrHZB9/lo8tJq5eJjtYjFQhvmxSLta1R73kU6Ta6uQoFWjgwjr5JqkP1bmvGlR7ahRcchJpgwvnodute0oMk502c3ny9r+GLtrWluBo7D9wezWJXLQpdNI6EWMGZpwibPJfz6+Vdw72/eXV3LuslJnH3jfQxOxRFzJtGa8WDQ1oSf/rofCXcr7nGKdDLM3Xy026/dbYjjaBP1wgWofD7qcDQ6Kjf/mI/ZF1MsJr92MEgrglCIImCkAmiJBAm9wXC1Owag96l2l4vUIUerrVkxASBnCR8+TJPx9HR113T3eMgwkCp/dnbSOVpDLq+57Nu1FvssM3t3ZjNw4IDcl8HjAVfG/Zn6ttABHGcsSKk0sGTiKLIs/Dozkiot7OND1bUcHRwETp3CuWO9EMUi2iI+NEe8MGVSiIPH373vkaNEbrZypE4HNDeDm53ccK3PnsuRy4Xj5OYUmQwJ9sQECZ5eT+6snh45wQmQQ0RzuStfU+oaVa37FrkciUouR5PSMtVAKRtr19J3kUzKlSOrEVGkiUfacN+4kaoX1rCYA6CVsstFhlRzM30foRBNsAYD8nNdoyWk7gXd7rRi0tgAZyKEIsMgIpgQ0hrQnE8urEPPcnH2LHDhAoRwEDlwyHA8EhotQloTTrtWIxaOyZu8i7F2XS4UdTqyiNTq+QXd76cELKkWdTwul8lta6MJZWKCMkqbmuS4bSm6RQr1m2ula7XVXZM+HqdIi3yePncNxaDPS3c3fW/pNPmn506w1cLwME04okiTqNtd3eHEN0NLCxlSLEv9ejdvxvlQFj/48Zv4wg/OYPd/fR0HT06U/G3rXtC/urcbPkczRACWdBwxrR4xowVbGjRyLZRKI3VpyWSgUXMQ8mmkNAJGLU2YtDRA5DjYTVpKjGKYRbsvcs3NdPHE4+Qmmb0xnMvJE1w8Tn7kcJgETiqr29JCYj01RZE4qRQtJSUBFwQ5zXk2s2u0VyOxGG1YZbMUFtrYWOkRLY2GBnkfxOerSG/LGyKKdMwjETl5zW6vnQ5RNwPP45DKje+HjDhpdCOhpryBJ184U3JRr3tB37+lGf/Xv7sHKosZzngIDQY1Hti9Ae0WLVUpXOYWUVdQLMqt4Px+QKXC6jUtSAsGXGhox/GmdTBmU3Dlkvj03lvkHqCLtHRFtZosa46j1xgYkOu7jI/TeDiOluoOBx2bQkHujORw0IU3OUlivno1bcZKEwPDkJVeY4L+q8Mn8caRizg9EcEfD6px8IPJSg9p6bS1kXUoTcDVRipF/v1olM6PxkYS9nqx0Ofw349MYEJrxpi5ETmW9u1SuQKePXSxpO9T94IOAPt3tONrX/44vv1AB/5qVQZbdTkSrqmp6zdoKDd9fbQknpggUVSp0GXksXljOzzrtiBgtMFlUOOP71qFB/bMVFNcatNrp5OW45KPfGCALPFgkC4qo5Gsd4eDhHxsjHzq4TA9vr2dnr96NYn7XL+50UiTxexORxxHm7hVKOgHj4/hp//yPgzRIGJqAW/rm/HkC2fw+uA8m7u1hFQaIh6vLteiRHLG5Sm1b5Nq7tepoA8lReQ5HnG1Dlle3iOQyjaUihUh6ABoA2/jRhKlvj4SsL4+ulUCUaSTeXRUTh4SRSAcxprtG/DzZ34Lvd/9Tfz9H92Du7d2yJ1zlrrByDB08fA8WdfJJIm6ZCXpdCS8UgecwUF6nCiS/12loqgJKcFJqt8iTYzSOOez0qtQ0L//4ikYI37ocml4TA6E9GakcgX8rxOhSg9tabS0kLsskZAzYKuJZJKMhUyGjAyzmc6RMqbFVxKbwwwRQFgwIq6R3UqlriZZn0dvPhoaqKbL7/0e1XPR62mz79VXK5NkJImyVLdFryc/Os9T/RYpomX1ajn9HiALXRTn7/W5UPR6ep18nlwwDCMv0aXqi8kkTYKSb93hIJ/nXFGWhF2KP9fp6DPUiKDHvAF0BifBiUWcd3RAZOiS8CXyFR7ZEmlro8lVFEk4l3K+lAOpTwHLkuGgVsvnXh3ylY/0gNFqIeRld6mg4kpeTbKKgrCXASmVfetWslqmpzH86xN4euIHeEvfAofDhK/u7cb+Zajw9uLRIbz0T2/DMDqIO/yX0H6XDbeGvSScHR3yAznuyifO9kcvpYGxzUYup85OspCk95E2pRIJstgtFlrFPPAATThSNygJvZ4uytmuq/ni0bXa+WPoK0wPkmiJeJDlVeh1rb78/wZ9jV8as6MspFLIVRIOePDEOE7+2c+w/WQvBDUH2xYGW3n+2rkUdcD+Lc3QPXgLfvjKGYyCLPNyaE2Nn7WLxO0GWlsxzmgw0DcGTjMBfZsNE2EVnnzhDACUVdQPnpzAd/75A2ybGIUqn8UIZ8Dge/2wGBPouHPH/GVoJSRBX2oIoCToweCVqe48TxOF1OGou5va3nm9dNzGx8kalyxztfrKSBeA/Ojh8JW102dPRFUk6J9frQWfiiKm1mPYRt+5oOLwO1vLFyu8LGg0ZBywLH1fweDlFmiV5ODJCTx14H18buQS+GIBE2oLXj7rR2hrDPdur50Kl4vhgR2r8ECjCr0qFXo2by7Le6wcl8tsHA6guRnHEhw02TS6/KPQZWlzohw7z3N59tBFMMkEGuNBZHk1PCY7mFwGJ7xJ8m9fL3RLquuyVPeFVkvvM19Im14vF/BqaSGB7u8nIWYYEnXJh88wV1vk88WjV2OkS7GI29UJrNazSFosCOtNaLYIeOaRTbhnVR2Ii9NJ31kqdXUNoArx7KGLsAam0B6aREwjIMOrEeB1+POTNb5nsRBmrmu2jNfAyhR0tRpoacEYb0Ke5dAemoQjLifalHrneS6T4RQaoz6IYhFJlRZFhkNapcE4b7iy2Na1KJU/2mYjS3zua+l0ZF3nchR10NZGLpWhIRLpU6eurAIpRbZIqwZBICEJzbpI1WqyFqtJ0NNpYGwMNoHHvfdsxflnH8Gvv37PsrjcloWuLvou4nHa+K4CJsMp7Bo+DUM2hSljA5IqLRIaHS5kV4CzYCaCh1EEvQysXo2C3YGkSgtbKgpXzA91nkLvyt3HsNmkQWvUizzHgxWLEHJpxNU6FGwOcrfcqBHEEmLRr0BqhzXXSp+9MQqQq6WlhTZoN2wgAT9xQi64ZTTSBuvsAlw2G00C+Vmbi9W2MRqPk39ZFGllNHe/otZpbSUrXcoYrYJM3fVcCrdOXkBQMCKot6DAcpgy2tFoM1R6aOVnptZ7sYybvytX0I1G7NpzC5J6A1S5LG6d6IUukyjLzvMVpFL4T6tFNCdDUOULMKWisKZiSOuNuONDPVfXO58PjYYs6KUWylepSIznCvrsjVGAxFna+HQ4gB07SKhPnaKfJhNZ87MF3W4noZz92tUm6OEw3VSq2q/fMh8zrkWwLFnolcy5mOEpvQe2bBIXHe0wZxMIaQ2IWZ211Tt0KbhcEMuYDbtyBR3Azt+4HRtvXYui3oBN3kF8auwY/nKrDvs3lyn1O58Hzp/HvelJfKxVi0Y2i67ABAxaFR58eBd2dNquvyEqUaqNUYCEN5MhS1WCZeUWdAC5SxoaKI08laJN1M5Oes7gIN0vJbFICALdZgu6INBEVMns3NkEAuRCEoTar98yHy4XTVQaDbm/Llyo7HhSKeyYvIButxF2NQNrMopkSzueenR7/bi5KswKcFxdh64urFvfjnVuE3DuHO7YbALy07Q8bW2lC72UfTClRs8NDei+ZQ26WwJ0kW3pBravpsJYC6nEJoUrptNLz6yz2SgKYnychFwqHavXX1nAy+0mARwfp8Sitja5MBfD0CQz1wK02+nxUqne2RNRNdTs8HjIrWQy1WfIHMdRpy5pZXThAuVgVMq1NDgIjI+jpc2Fz2/eQJP7479bn6ujCrGyBX2mrCxEkazM8XE5ikQKu9u+vTTpyIUC1UDJZkmQAwGyXmdHjnCcnGl5PUoZMcIwZG0Xi5S1ynEk8jodiXUmQ+PleUpCGhuT+5tarRTOaLHQmOJxeh0p289mo88VDNJzZ4+70oKeyVDJBcllVMtNLa7DS1kTeE8ODu8U3j/wNprW78GDu9fe+InloK+PzpFNm+icam+vye5Q1cyKdrkAoJPK6STLdGiIrMxYjE683t7ShXtNTspp9NksiQnHAbfeStZuMChnWd4IlqUJqFSbXAxDtZuNRippmkrJoYder/y4hga6EKWJT2pVJ4n+7I3RRILGaTLJbheNpnra0aVScslct3thK6Ma4+DJCXzt1RH06hzgCwVovdN49vmjZSnbuiAuXSLjyeGg66urq25T/SuFcjSlhsDt7XRRB4Mk4qJI7pFS1HrJZoFjx8hC3bBBtlTsdrKOf/M35Q2shbLUIl1zYVmKYmEYckVIrdgkvzlA97W00PsGgzQRqtW02lCrSRyTSdogvXiRXEyS6Mfjsmvmek2ql4toVI5waWurmizKUvLsoYuIF4BBewvyHAd33A8hEip7nsW8SK3mtFo6T1QqMiIUSooi6FLd6FWrKOVeiiBhGHIlSOK+WIpFakF16RJNGuPjwAcfkBunqYlEXa8H7rmHrPWFIggktEuNdJkNx9HxCIXoGEgRErNXKRYLvbfXSysKg4HGwXF0zEIhEnupNZ3VSq8hxaSbTLJrppIEg3SbCSWr1abE10PKpxixNiGm1sOYSaIt5Cl7nsW8jI2RcWA207llt8thswolQxF0u5068LS0kD+vUCCXi1QiNhRaWrjX2bPkm25rI4H7l3+hk1uq/SxtxqnVN9ff1GQiUSx1YTGnk8TY56PxNDfT8ZidJNTQQJZ4IkHjTybpuOn19FiPhx6XTMoFv6RxGo1yk41KEgjQZ9JoaJVUys3vKkHKp/CaHAgLRqiKObSHJ8qeZzEvkjvTZCILfXY9IIWSoQg6w8i77B0dchOJbJZOukKBol4Ww9QU8PrrJGIPPUTikUqRIG7eTD8Xu+w0GmnsUpXDUqHRkFXt89Fndzho4pEaYAByz1Gvl34vFMitotfTT8llk8/TcTQYZNGXxj1fM+nlIp+nz5dIkMVYy02hr8NX93ZDUHHw6yyI6IzgisAm/xC+ek8FXB1SMxlpg3316hs/R+GmUQQdIKtBKmbU1kb+vkCAhE2vvzLNHbixC0YUaRI4eZJE8P775SaxAL1PdzeJ+WILJrEsiWM5kkWklUogINdPz2blzU2OI8tcsnBn+8ULBYoQamqiv1MpOXJH2ijV6yub5CJ1y5ldi7sO2b+lGc88sglOhxmjFjeyOh12cQnsv/AmbQgvF8kkbbZLEWQOh+JuKRMrO2xRQq+nE43nSbzGxsgizeVI3KRwQ7WaLMuBAVqmX0MI+Kkpui8ex3sFA7753BloRoexPTGFR4xpdHfPlDY1LrEAlNksNwlYSinduej1JMLT02S9Go0k2sGgHN7X0EAWeipFFryU5s8wtOElxTonk3QMGUZuOm000uqlUglGUrccaUO00iGUZWT/lmZK2unwAq+CjvvEBHD4MG3GL8fewdQUnUssS8fa7a7rY15JFAsdILHR68mabmkhMRscJHdGeztZlq+/Tq4Eaek4O5xvNoEA+FAI4Dgcu+TFn5+JYtt7r+G2kdNQez14cSyJ4+mZkrNLFfS53YJKictFk5iUXGS10mQmtZrTamn8gQBNLOm0PLFYLCT+Wq3sR5eKRM0ed6XcLpEIiRrLkqDX6TvEhgAAIABJREFU4YboVTQ0yF2BtFo5ims5GBmh80iq8NnRUZd7FtWAIugSUlEsrZbcBZEILRNbWuQl41tvkRXqcJCIzuoCc/DkBD76pz/DZ77wP/DsSwN4+53z+NWxS2j1jaI9PAVLOg5HMoJRvQP/esFH1v5Ckoiuh1ZLAlpqPzpAF79GI5cEkJbIs0XA6ZRXLvk8iTrH0QUbCsmCDsi9SkVx/qYYy4nPRxOOWi2HatY7bW1y0bWpKfoeluv4j47KkVA2W13G/FcLi3K55HI5fP3rX8fk5CQ4jsMzzzyDVqnJ6wx33HEHOjs7L//9gx/8AFw1V7MzGumET6XI7XL2LHD+PMWN8zwJlNVKF7/LRf7wwUGgWMShoRj+7NV+fPSDXyIiGDDMuvHGz99CiNPAmEmi39YMRhRhyKYR1RqRzGdlUVsqZjONZXaGZimQPufoKFnWBoNcP12qe2I203uyLB2jaJTE3OmUQx+zWRJ7yYWTTNJnNxrJQl/uxBJRlD+TTkdCtxLo6KDjLm2kZ7NlF/SDJyfw3146g4f+7RXs9njR2G1AR0tL3W5CVwOLEvQXX3wRJpMJ3/3ud/HGG2/gu9/9Lv7iL/7i8v2iKMLpdOK5554r2UDLjpSlKYUTOhzkOz9zRk6k2bePBOncOXK9TE4Cra145RcX8ZHRftiSEYyYXHAlwghxaniMDljScQR1JmiKBRzpuAWqYgHrkFi6u0XCZCL3j+SfLiV2O31Gj4cifqxWclVIrhXJVZVK0U9RpLZ1ej0dR2kFk0xe2TxaEvRIZPm7F0kJLtksiflKsRZVKgpBPXWK3B+ZTFldXgdPTuDJF85AF/bDkQihkM7gTW8GfXkzHqjDJK5qYVHm0ZEjR3D//fcDIEv8+PHjV9yfTCZRqJaKeguFYUhkVCpAEHAqp8E/vDuC7/0/f4MXfnkafYEU8O67tFw/dIgEXRCAVauQCEfREZqCLpfG9qmLgFhEQDCjPTSJtM4ARmSQ4VQYtjUjaHdj/72bSyfo5QpfBMh6djrptdPp+d0uBgNZ4Dodiby02SUlkIgiCboUDjrHj84td9ao30+TEsOQpVjHjYmvYtUqOSEtGi1rLsCzhy4ilSvAlozCkQhDhQLGDA14+lwVlH2oYxZlofv9fthmLm6O48CyLLLZLNQzM28ymUQgEMAf//Efw+v1Yt++ffjsZz8772v19vYucuhAOp1e0vPnwvn9UE1P4/S5CQwcHQAKIkyZFGJiAa8ORuDR9qJzYwqi2YyiwQA2FkNuZAROLo80r0ZKpQZfLILP52AUgY5cBC3bNuLt/ghOGRrhMGrwhTYGHdokekdHS+ZuUE1Pgx0aAo4fBwoFiGo1svP4hhd1vPJ5aIaHUQiFkG9uhtrjASYnkZ1JDGHjcahHRsAkEuCDQcSbmyFqNGADAainpoB8HoVQCLlQCPzUFLhoFJmZjVXN6ChyanVJv8MbIRw5AsvgILhiEUmOQ3BwcF4feqnPrVKxlHHx8TgsHAet349UXx/yhQIiq1cv+Tycb0xSNqo76oUtGUaGVeO8sxNjsdyyHdd6/A5vxA0F/fnnn8fzzz9/xf9OnTp1xd+iKIKZdVEIgoAvfvGLePjhh5HL5fCZz3wGW7duxcaNG696/Z6ensWOHb29vUt6/lVkMsDAAL796jjawcKYzyGt1iArauDXGvDBRBIfvt0IPPggXQRHjgBnz+JTDSJe8woIaAwwphPoCk9h0NWB3d0urL9tHT6yOQ3s3UsujN5eem53CQv6t7bSyoFlyVoOhyk6Z471uejjJfnpu7roM4yNUdimVksRPyxL1vjAAFm9uZwc7aLVkmXf00PW/vAw+XNn2tQN9PZidSm/wxtx6BCtIsxmGG+7Da716+d9WMnPrRKxpHE1NVEdoYkJeMMZvPLmIL7vOwtzY8OSOtDPN6YmyxQmwim0haZgyKYxaGtGv70NTRZh2Y5rXX6HwFUekdnccGr+5Cc/iZ/85CdX3D7+8Y/D5/MBoA1SURShmuULNRgM+OQnPwm1Wg29Xo/bb78dFy9WoCDQzaLRAOvX4xXLagw62hDT6mHIZzBsbgQrAoHczEbh1BS5DNasAQD0bFmDOz7UA5XZCpVYhGjQ4bFtzVjf3kDLWquVbpkMuR9KnchiMJDAtrfTDSitC0YqcerxyBvDM98/OE52t/A88P77tN/A8/IEk06T8Et+dGmpr9OBLUXnpYUSj9OEwjAra0NUQqcDWlsxHs/jRL8X8VQWXYExTIRTePKFMyWtwvjVvd0QeBbdvhFwYgEjNjcStoaV05moQixqrbV79268/PLLAIDDhw9j586dV9x/8eJFfO1rX4Moisjn8zhx4gTWzIhfLeC26uE1WBASTEjxGkT0JgR1ZmhMBrI683kSN7Ua+PCHgeZm3NrTgv/4wGp88dFdePgjm7HDNdMFJxAg4WBZ+h0ob5Ycz5NwllLQNRqyzH0+8onbbGSxS/1CDQb5/8UihQKazbQfwbI0kaVS9Doq1RWCDlFctnK6v/rZYfzirV789IwHL45n8Etfje3zLBWVCrBacSLOQpNJIQcWrngAjkQIqVyhpFUY929pxnfubkZX0o8Mp0KwZRWe/sStSmeiMrMoQd+3bx+KxSI+9alP4Z/+6Z/wxBNPAAD+5m/+BidPnkR3dzcsFgs++clP4tOf/jT27NmDzZs3l3Tg5eSre7sRtLkQ1+jR52gHwGDU1YHfWO8it4JKReIWCJBISeVwLRZg82Zw0ShFFDgcJHqShRsMkmVf7l1+s1kuY1sqpHo3k5P0eYpFuZSBZHnfdhutEKS2bmq1nF0rxaPr9XKZAGkDVbqvjBw8OYHXfvIKCsk0GBEY53T4j0d8lasNXikcDgxpTNAUctDnUiiCgSsWAF/Il7wK40dVYdxuErFlfRue+g/7FTFfBha1KSrFns/l85///OXfn3zyycWPqsLs39IMJrcbP/lRFqPJIj4cH8Nj963BVnMRF154Gf94aADTaREdbAYf7bZjqzZLXVi2bweGhsDkciSqU1MUhSL128xklqfdltlMkRxSZ6FSoFbTimN6mkITpXBJl0sW9EyGJrVAgBKypG5MUj10p1NubZfPA2o1RI5bFkH/6xfew96gD5pcBgWWxaCtBR6VAc8euriyhMZiQd7WgFxgHK5YEGddAANAyGdgspco8krivffoHFy/norRKZQdpZbLNXj4tk48fFsnidHRo0A6jTfGkzjdF0RzzoNVuQx0uTR+EV2Nwp712NHaSiKXy6FgsZAroVgkwUsk5M3D5Yh7ntlwRDRa2tZqjY1y2J/LRTXeg0Fyx2g0smiHQnLzDbVarocuZYkCl3t5Fmdnk5YR1cgInIkAVMUc/IYG9LpWAQxTmdrglUQQcOeezRge74cpEYUr5seUxQWrmMd/KKV/2+Oh64ZhqNb/SgoPrSBK6v+NkKowNjbi/53Q4pcdWzFmciKm1cGvt2DY0ID/fSlBVurUFGCxINvSQtarwUAWeSQiZ5ouV2ak2UyCvpTmHHPheRLycFh+bak0gMFAgm400oQSidDFrFaTayUSIeGW3CwzbhdRiosu5TjnYXtqGvZ4BAWWw5ClCR4j1aGvSG3wSqLT4bZtXdi8ZTUENYe28DTcAof/eG9H6VYqogi88QZN/FYrcNddpXldhRuiWOgLYabZsTcUR8DeAncsgJRGQESth6aQQyCRJp+6IADt7SgGg/S7RkNWrSR6y9lZXgo1lES2VLhcJMxeL4n0wIBcaCwQoBVNQwOl14siibrkAopGyVIThMuCXpQKRaVS5avAl8ngk+YU/MUcJnU2jNqakFQLEFTcyou6mKl3v66nHetyEdwN4PP3uYFVJYy8ikapLIaUjauk+i8bioW+EGw2gOexM+cHJ4oYsrlRZFjkeBV4sYAOPk+PWbcOcDhQkKramUxyRItaXVphvRFSsbFSZ5BKvUdvuYVuWi25XkZH6b2kXqk8LxfBslrJ/STVDpm1MSoKMxZyOd0uExPoKUTR5TIgZ7bCZ7DD2mDFM49sWln+c0AO2XQ45PNxYIA2rksVPurxUK6CFPGklMpdNhQLfSGoVMD69fjM3RHEXjqFtCjivLMTBZaFCXns+fidlDgzk1xV1OtlV41OR7flLugvNcCIRGiDstRIjQpWr6YLl+dJ2AcHgR07gLVryc3U30+fXaMhcZda1fn9QCYDUa2WywOUi/PngUAAbrsRn/nYPfjMPffQ5LtSEQTavNbraa8jkaD68Js3///tnXlwU/e1x7+SbK2WLMnybhZjAnbYg/PAYTMBAmkSQhuY0JS8QGgyKS5lMiSTpCFt/miTtKX0nzDzJg3Qhr7JvMLQPPpKA1nwmyTsmPIgGIPBGLxLsmRZkmVrue+Pw5Vs8CILSVas85nRWNLVvTr+Sfrec8/vnPOLjvhevhy6Yi0piX8DtiSGBT1cUlOxdNUC+FQq/PUfZ3FOmg51bja+t/BhLCmb2Ld8XCLpu8TWSFWr6XT0Q/X5hrde6XBQqUgUSkootFRXR4tdq1SU+WOxhJavczrpJk6QiemLKlXsBD0QAKqq6OqgoICqVO+1bfF3HXEBk7Q0+nxSUsijdrnuXdB7eoCaGvo8VSqgn+pwJnawoA8HiQQrVjyIFUtmxr9LYCT0zvOOdidGEaUyFNYpKqI89fp6ahOg1VJeutcbiq+LwiqV9s1Ht1hCMfdo0tJC7RZ8PrpqyMhgQRdbMuj19JnI5fRZ3Lhxb/Fuj4euysTVqLTa2FwdMgPC10KR8F0Qc4A8JCC2lZgqFQmxx0NioNdTPN3jIaFITQ0tNN3eHsp+0Wjox+/3k6AHAiQu0ebECTrJ5ORQt0GxkjaZEbOVxJbRLhd9PpcvR37Mzk5qK/3553QSFfvpR7pmLhMRLOijGbG/eyzj0+LybR4PiUJODt1vaKDnRcEfP55OLN98Q5NwNhtw4QJSGxpiVzHa3Q188QVdDYwfT965Uhm78NN3iby8UMaR2Uwn1uvXI+uRHgjQVZlMRr18GhrouEVFIaeCiQss6KMdMc87VoiCLr6HwRDKQxdbAHR1USMzr5cu60+dIq8wNRUyqzXUyCvagn7sGKVXZmVRPUBqKnvnIjod3cTaCJmMPOyPPx52tktKWxt9nvX1tICGx0MnjAkTkmN5vwSCBX20o1bTDyxWHQ2lUspgEStDxcm2QIBCLCoVeX82G4VeJBL68U+aBOTnQyrGXXvlpkeNixdDPXbEsA8LeoicHKoRMBqp95BCARw5Qn31w8XlQorVSieEY8doriQ/n8R8woTY2c70Cwv6aKd3jDuW7yF66FptKE3N6Qxdcjc00ATpAw+QZ+j3k8fs85HwazTkoUerYtTrpZiw00n/u8NBgsWCHiI3lzx0sU3F0qX0/D/+Ed7VkiAA9fUQxJqDEycorLVuHfCDHyR3augIwYI+2olHR0OlkuLVgkACLsaoPR6639FBAjtlSmhx6NOngZYWyOx2El2ZjEQlWieeujqanFWrKX4+eTJ5jApFdI4/GkhLoyyUvDw6CZeUkAjX1dHcg9geeSDMZjqR+3yhk8Ds2dSkTlzEhIkrLOijHYWCfqzxyHQRs1S0WhJnQaDCIrebPPL0dNrm94faAIhFReK+0Qi7BAJAZSVd/ufl0W3y5PgXdyU6SiV9JuK6sXY7Vf+aTBR2OXs2lH1040Zfp8Dno+whux2q06epeGvMGOCxx7gR1wjC0/3JgFodn0yXri66L/Z1aWigE8rUqRRWqa8PFZrY7YDBAIm4AIbTGUqhu5cOkT09VLH6r3/R/52fTzH0aK8SNRpQKslL7+4mQW9uBvLzcUadjRPHr8D61Ye4MOsKNj2Yg4cnGmmuY/x4CtM0NdE+9fWQ19fTcebNo3HmUv8RgwU9GRBzwGMVbugt6AYD/bilUvLEdToSAbWaqhHFDpQOBxAIICCXU5ilpyck6PdCUxMJk8dDJxbRHvYa70Yioc8uNZWuXtra8PXlZlSeqkOHRIVcjxVzThzG3pYp6P735Uh1OvDxhx+jw9aJHLUM60oMmNvTRt56bm7oZM2hlhGDBT0ZUKsBsxmSnp7YHP/OTBeVirxsuZxi54JAl/VOJwluRga9prsbAbF1gN9Pl/c+H92XySKzxeGgq4OODvLMjUb6+10pBos3SmWoK6jdjk++rkcgRQGDuxMuuRLabje6vX4c3v3faFTpYbCZ8YC9GbKAHxcuulCgdkGjkFCLXJ2Ojse9W0YMHvlk4LbHJIllpotSGYrT9/SEPG2fL3R/3Djywru6go2hBJWKXu/zhTJcIg0PibH45ma6n5VFdkVzkY/Rhuihy+WAwYAuhxMXsyfCptLiQnYR6oz58EkkaJYqIevqQqrgR5PWBJtCgwJrA641O+DLyAjWFbB3PrKwh54M3P6RSWOduuhwkDDX1pKw9vTQZJrRSL1dZLLQEnW3i48EuZwu/R0Oep0gkDhE0mpYbN/b0kLvlZ4e6izI9I9aTR51ZydgMMCkkqLL3gqfNAUeuQouuQqTrA24ahqDtB432jR6QACMbjvalTocyypEjykHJ09bcbL6GyA/HxtXK5KvLXGCwIKeDEilgFIJibiocywQF6q4coVCKBMn0vs2N1OBj05HHrNOR+luMlmoGZdaTZOk4gIgkcbRHQ6gsxN139aivrYd/+NvgaNBiyeypuAJcaFupi8GA/VcsVoBtxvlKx7Epwe+htpNqzulux0ocjSj25CBTEs9HAoNnHIVprZeQ5M2E5cnzkDX9Tqcy5yANp0JTp8Cbxy4AAAs6iMAh1ySBbU69h46QOGUceNIuNPSKA1OJqOMF9HzlkjoplJRGEirpVi7202vjSTk4vMBTifOVNWi7tINmGVKmDUGNHYJeO3QVXxyrjG6/+9ooqCAagRcLiweb8CjTy2EQpeGAkcb8lP8mL1gBlYunQ671gSj24Ex9lak+v3o0BkhlQJeAWjTZqBToYEgkaLL68fvDteM9H+VlLCHniyoVJB4vfc24TgYSiXFYU2mvkvtabXkAQYC1Pzp/vspfu73AyoVpGJrX4mEWujq9RSq8XqHN5HpcADt7fj2i5OY2ONBY2YRuuRKmNV6uH0Cfne4hj3GwRDL9ZubUT6zCOUZS+kzslopdKXXAyvn4l+ffAFHlxdOUzZWLJ6OG2da0ZWqQKeibxZR0i2+nSCwoCcLvStGY7EUnlRKC1rcSVpaaHUjq5Vi6jpdMNtF4vVS/NZopPCM2NDL5Rpe7Pv2akkaawvsCg2atRlwy1WwaOgYLDBhMHYsrf3a0kJXWh4PVY7eugX4/VigBxY8W04hs7Q0YNIk/PPcTZxXGu9qwpV0i28nCBxySRbuXCUoXohZFH4/eYF2O3nfAKDVwq/Xk4hkZJDwNzb2zYwJh6Ym8v7lchgFL9rVeriVajTpshCQ0tUIC0wYZGTQ1VR7O52gnU666jIaQ+uO5uTQ5zR1KpCejsem58Fi6Ds/kZSLbycILOjJgkxGGSWxrBgdCL2eMlAyM+nE0tlJ6YsOB7rFpfqsVhIUj4da3ooLSg+FIFDPltv57MVGFVxp6XDLVWhX0SpNLDBholQCDz1EbRJmzaLeLpcv0+dhtdKVXVMTbc/MBLxePLhuJdYun4x8vQoSAPl6VXIuvp0gcMgliQgolfH30AHKpDCbSaRNJmoBoFQCDgcC6enUa8XpJNuUytCCCYWFoSrUgXA4yOtXKICaGuRlaDBn8jQc68mBRCJBnl6FV5dPZoEJF5mMPq+ODmDBAlqQ5OTJ0GpPkybRrbaWwnhjx+JhtxsVj88ZacsZsKAnFYHeRTzxXLUnLY3ez2ajOO2tW2TDtWtI1WioqVNqKnnn16/TItMWC8Vyx4/v9Q/cbvglk4Xy3C0W8hrdborPp6djVukkfLRoEYUHmOFjMlHYxW6n+Y6xY0nkbTYKl928SSG08eN5AYsEg0MuSYQgphbGO+wikYTCLlIp3bdYAJcLqU1NJPhqNVBWRq9vbibBvnatb3/0a9doRXmAYu3V1STmHg89bmsjERdbCzCRodXSFY8o3HPmkMgvXEjxdLeberfwGCcc7KEnEQGxb7nLRZ5XPBFFXKwIPXUKwO3ly2w28r71esqqsFgo1NLYSNuMRrJZjKuL5f3Xr1M2hkRCMV6plIRH7CnCRI7JRONvNJInPnZsqEdLVxeLeYLCgp5MyGQkdCMxMarT0fvbbKE+6BMmINDZGeqOmJcXyoTJygIuXULloeN4s0GF1LprGK8Q8NycMVica6W4rujJ63QUIigooApVsTcJEzkmE4XFxLBV74ZbLOYJS8Qhl1OnTqGsrAxHjx7td/vBgwfx1FNPYc2aNdi/f3/EBjJRRq0emYlRiYQKVOx2io3fvmT3ZmbStitXgL/+lSZDW1uB7m5UtXXhy78cAq5fQ1q3G1e9KXj/y6v4Zt9ndJzMTArjiAsTT5hAE3nFxRzbvVdSUugEGc+5FuaeiejTunnzJvbs2YPZs2f3u93tdmPnzp3Yv38/UlNTsWrVKixduhR6bpI08mg05M0OtxIzGhgM9N5+P1WM3ryJVLOZPME5c0jUHQ7y4i9dwv9W3cAYSxPKAwHUGfMhEQQYO8y4cuAY5hXfzo7pXVU6ZQp59gyTpETkoWdmZuL9999H2gAL7p4/fx7Tpk2DVquFUqlEaWkpqqqq7slQJkrEY43RgRDDLloteddGI3yZmZTlUloKzJ1LQl9SAlit0Lc0wOSy48FbFzG1qRaqbjemtdWisL4GtTYPfvnlDfzH9R68pZmOiwXFFG5hmCQmIg9dNUQMzWKxwNhr/UaTyQSz2RzJWzHRRhR0lyv+y7JJpVSUIl4ZFBbCI7bcdTopfl5XR/nnABqu+uBqrEdhewMyujqw6PpZ6Ls60aVOw186NLCneiFT63FSNwae5ltYYpNhRXz/I4ZJKIYU9H379mHfvn19ntu8eTMWLFgw4D5C71Sz248lA8Q0q6urw7GzXzwezz3tHysS2q6aGsibmiCYzfB2dIy0SfAEArjW2Ej2jBsHeXMzpA4H/Ho9Ji+cho/+bzzyrE1Q+LoxpfUaWg1ZSOm0wGRvg0RjQJdcAX23E43KdLz7dQPGjYvOUnMJ/RkmmF2JaBOQnHYNKehr1qzBmjVrhnXQ7OxsVFZWBh+3tbVh5syZ/b62pKRkWMfuTXV19T3tHysS3q7bi0sgAWysrq5G0cyZlHZYXEzhl5oawGTCxHwr0id247/+KYPxykX4MzIxb/FsXPnPT6DtdsEvlaJWMxYuOV113PLKojbuCf8ZJhCJaBMweu06e/bsgNtiUlg0Y8YMXLhwAQ6HAy6XC1VVVSgtLY3FWzGRoNFQmKO7e6QtIXQ6ykN3uajIyGikXPTubiwvUGP36hK899gkVHxvGuZkK9GePw42lRYKnw82JeXTN+sykWvghaCZ5CYiQa+srMSzzz6Lr776Cjt27MDzzz8PAPjggw9w7tw5KJVKbN26FRs3bsSGDRtQUVEBbSxatjKRIcbO7faRtUNEq6VilWvXKKVRoaBYv9dLMfWrVynnXKEA5HLMePoxNJsKYE3Tw6w1ojZjDLrSjdyAi0l6IpoULS8vR3l5+V3Pv/jii8H7K1aswIoVPEWVkMjlJJg2G7VLHWlkMkpnvHmTJkfFlrvjxgGnT1Njr5QUYN48QKdD+cWLkK58CHvPteCKdjzyDGpuwMUw4ErR5EWvpz4oI5GPficuF+WPp6dTyGXMGOrm53BQbrnFQp66zRZcBGPh8jIs3FjQt4KRYZIc/jUkKwYD/U2EsEtHB51UZs2iWHprK3nsjY1UcPRv/0Yi39lJIj9xYt/eIgzDAGAPPXlRKulms1GRz0hit5NnrtNR+X5tLU3cKhShZmJFRbS9u5vb4jLMALCgJzN6PXnD8e6P3huxP3tBAT1OTycPXCKhCdKmJmrklZND3rnYYIxhmLvga9ZkxmAgsRzBAiOZ00l3eletpqeTN56XR/ZJJH0nb7nbH8P0Cwt6MqNWU8aLzTZiJkg7Oym00p/XnZYWEnetNhQzZw+dYfqFQy7Jjl5PaYHxCrv4fCTMUikQCEA2VE+ZoqJQK1yNhkI0PBnKMP3Cgp7smEwk6PX1JJ6xpKuLyvplMspSEQS6DdZWuXcPoLFjqfUuwzD9woKe7KhU1OWwoYEmSGNVaOT1UvaK6J3X1gIpKRCkUgqthAOHWhhmUFjQGRJxp5PyvjWa8AU2XPx+Kt/3+ah9rkpFJ4/mZgR0Ol5diGGiBAcjGWL8eJogvX6dhDdadHdTjxaPh0I6ajUJeE4OMH06vLm50XsvhklyWNAZQiajoh6vl7zne8XnA27dAr79lgqDxo2jbJU735MnOBkmanDIhQmhVlPr2rY26q0SSY8XsXS/tZXum0y0IPRI94thmCSABZ3pS14e5aW3tFD/lOFgtVIc3uulzJX8fJ7IZJg4woLO9EWhADIyKJUxO5vi6uHgcgE3btCk6oQJ0Z9YZRhmSDiAydyNOFHZ3Bz+Pi0tFBO/7z4Wc4YZIVjQmbuRy6kDo9VKlZlD0dVFHROzskjUGYYZEVjQmf7JyqIqznD6pbe0ULZKVlbs7WIYZkBY0Jn+USjo5nAM/rrublqMIjNz5FrwMgwDgAWdGQydjipIBWHg17S03N3elmGYEYEFnRkYrZbK9l2u/rd7PBRnN5k4z5xhEgAWdGZgtFr629nZ//b6eoqdc/k+wyQELOjMwKSkUPVof3F0i4XCMQUF7J0zTILAgs4Mjk5HIZdAIPSc10vtdrVaCrcwDJMQsKAzg6PT0aRo77DLzZv03LhxI2cXwzB3wYLODI5GQ3FyMezS2Ei56bm5lNbIMEzCwInDzOCIKwp1doZWNTKZqJ85wzAJBXvozNBotVTe39pK1aAcamGYhIQFnRma9PRQ8dBwW+oyDBM3OOTCDI2tV3CDAAAJAElEQVRKBcyYwY23GCbBYQ+dCQ8Wc4ZJeCIW9FOnTqGsrAxHjx7td/v8+fPx7LPPBm9+vz9iIxmGYZihiSjkcvPmTezZswezZ8/ud7sgCMjKysLevXvvyTiGYRgmfCLy0DMzM/H+++8jbYCVadxuN3vkDMMwcSYiQVepVJANElN1u92wWq342c9+hrVr1+Kjjz6K2ECGYRgmPIYMuezbtw/79u3r89zmzZuxYMGCAfdRqVTYsmULnnzySXi9Xqxbtw4PPPAApk6detdrq6urIzCb8Hg897R/rGC7wicRbQLYruGQiDYByWnXkIK+Zs0arFmzZlgHTUtLC+4jl8tRVlaGmpqafgW9pKRkWMfuTXV19T3tHyvYrvBJRJsAtms4JKJNwOi16+zZswNui0naYk1NDV577TUIggCfz4eqqircd999sXgrhmEY5jYSQRhsfbH+qaysxK5du3D9+nUYjUZkZmZi9+7d+OCDD/Dggw9i1qxZePfdd3H27FlIpVIsXrwYP/nJT+46zmBnGoZhGKZ/BsowjEjQGYZhmMSDK0UZhmFGCSzoDMMwo4TvRHOuU6dOYcuWLXjnnXewePHiu7YfPHgQf/7znyGVSvH0009j9erV8Hq9eP3119HU1ASZTIZ3330XY6LYKXCo41+8eBG/+c1vgo9ra2uxc+dOmM1mbN++HTm3+4k/9NBD/c4vxMImgFoyFBYWBh//6U9/QiAQGNGxAoBDhw5h9+7dkEqlKCsrw8svv4zDhw/HbKzeeecdnD9/HhKJBD//+c8xffr04LZjx45hx44dkMlkWLhwISoqKobcJ9Y2nThxAjt27IBUKkVhYSF+/etf49KlS9i0aRPG3W5nPGnSJLz11ltRtWkou1atWgWtuJg4gO3btyM7OzvmYzWYXa2trXjllVeCr7t16xa2bt2KwsLCuIzXlStXsGnTJqxfvx7r1q3rsy3m3y0hwamvrxdeeukloaKiQvjyyy/v2u5yuYRHHnlEcDgcQldXl7B8+XLBZrMJBw4cEN5++21BEAShsrJS2LJlS1TtGs7xOzo6hGeeeUbw+/3CgQMHhD179kTVlnBtCgQCwve///1h7xdru9xut7B48WKhs7NTCAQCwurVq4WrV6/GbKxOnjwpvPjii4IgCMLVq1eF1atX99n+6KOPCk1NTYLf7xeefvpp4erVq0PuE2ubli1bJjQ3NwuCIAibN28WKisrhZMnTwq/+tWvomrHcO168sknh71PPOwS8Xq9wtq1awWn0xmX8XK5XMK6deuEbdu2CXv37r1re6y/WwkfchmqzcD58+cxbdo0aLVaKJVKlJaWoqqqCsePH8eyZcsAkFca7Yya4Rx/165dWL9+PaRSKVwuV1TtGI5NA7VkGOmxUqlUOHjwINLS0iCRSKDX62G322M2VsePH8fSpUsBABMnToTD4YDT6QRA3lx6ejpyc3MhlUqxaNEiHD9+fNB9Ym0TABw4cCB4pWI0GmGz2WL6XQrXrv5siPVYDec9/va3v2H58uXQaDRxGS+5XI4//vGPyMrKumtbPL5bCS/oQ7UZsFgsMBqNwccmkwlms7nP8zKZDFKpFD09PVGzK9zjezwefP3111iyZAkAEtXPPvsMzz//PDZs2IDLly/HzaaBWjIkwliJJ+wrV66gsbERM2bMiNlYWSwWGAyG4OOMjAyYzWYAgNlsHvD7NNA+sbYJCI1PW1sbjh07hkWLFsHtduPs2bP48Y9/jB/96Ec4ceJE1OwJ1y673Y6tW7di7dq1+MMf/gBBEGI+VuHYJbJv3z6sXr0aAOIyXikpKVAqlf1ui8d3K6Fi6JG0GRDuyLoUBAESiWTA56Nl1/nz58M6/ueff47y8nJIpXTunDt3LqZPn465c+fizJkzePXVV/H3v/89LjYN1JIhUcbqxo0b2Lp1K37/+98jNTU1amN1J4P9v3duAxD179NwbRKxWq146aWX8Itf/AIGgwHFxcWoqKjAkiVLUFdXhw0bNuDIkSOQy+Vxs+vll1/GypUroVAosGnTJhw5ciTmYxWOXQBw7tw5TJgwIXgyjMd4DcdmIPrfrYQS9EjaDGRnZ6OysjL4uK2tDTNnzkR2djbMZjOKi4vh9XohCAJSU1OjZtfrr78e1vGPHj2KH/7wh8HHvSc7SktL0d7eDr/fP+hVSLRsGqglQyKMVUtLCyoqKvDb3/42WBYdrbG6k+zsbFgsluDjtrY2mEymfre1trYiMzMTKSkpA+4TDQazCQCcTideeOEFbNmyBfPnzwcAFBUVoaioCABQWFgIk8mE1tbWqE5oD2XXM888E7xfXl4e/D7FcqzCsQugAsiysrLg43iM13BsjsV3K+FDLkMxY8YMXLhwAQ6HAy6XC1VVVSgtLcW8efPw6aefAiBRnTNnTlTfN9zjX7x4EcXFxcHHO3fuxOHDhwFQeMFoNN6zQIVr00AtGRJhrN588028/fbbmDJlSvC5WI3VvHnzgse9dOkSsrKygl5cQUEBnE4nGhoa4PP5cPToUcybN2/QfaLBUMd/77338Nxzz2HRokXB5/bv3x8Mm5nNZlitVmRnZ0fNpqHsam9vxwsvvACv1wsAOH36dPD7FMuxGsoukQsXLvT57cVjvAYjHt+thK8UDafNwKeffopdu3ZBIpFg3bp1WLlyJfx+P7Zt24YbN25ALpfjvffeQ25ubtTsGuj4ve0CgLKyMhw/fjy4X0NDA954442gqEYzpSscm/pryTDSY6XX67Fq1ao+47B+/XpMnjw5ZmO1fft2nDlzBhKJBL/85S9x6dIlaLVaLFu2DKdPn8b27dsBAI888gg2btzY7z69xSKWNs2fP7/PdwoAHn/8caxYsQKvvPIK3G43enp68NOf/rSP4MfarmXLluHDDz/EoUOHIJfLcf/992Pbtm2QSqUxH6uh7AKAJ554Anv27Al6ux0dHTEfLzFdubGxESkpKcjOzsbDDz+MgoKCuHy3El7QGYZhmPD4zodcGIZhGIIFnWEYZpTAgs4wDDNKYEFnGIYZJbCgMwzDjBJY0BmGYUYJLOgMwzCjBBZ0hmGYUcL/A11WfkfU5TRNAAAAAElFTkSuQmCC\n",
"text/plain": [
"<Figure size 432x288 with 1 Axes>"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# データ生成関数\n",
"f = lambda x: np.sin(10*x[..., 0]) * np.exp(-x[..., 0]**2)\n",
"\n",
"# -1から1の範囲で50個の点を打つ(x軸)\n",
"observation_index_points = np.linspace(-1., 1., 50)[..., np.newaxis]\n",
"# 関数fを使ってデータを生成(各点それぞれに正規分布に従うノイズを乗せる)\n",
"observations = f(observation_index_points) + np.random.normal(0., .05, 50)\n",
"\n",
"# 10番 ~ 19番のデータをあえて欠損させる\n",
"ind = np.ones(50, dtype=bool)\n",
"ind[10:20] = False\n",
"observation_index_points = observation_index_points[ind]\n",
"observations =observations[ind]\n",
"\n",
"# パラメータの指定\n",
"amplitude = tfp.util.TransformedVariable(\n",
" 1., tfb.Exp(), dtype=tf.float64, name='amplitude')\n",
"length_scale = tfp.util.TransformedVariable(\n",
" 1., tfb.Exp(), dtype=tf.float64, name='length_scale')\n",
"kernel = psd_kernels.ExponentiatedQuadratic(amplitude, length_scale)\n",
"\n",
"observation_noise_variance = tfp.util.TransformedVariable(\n",
" np.exp(-5), tfb.Exp(), name='observation_noise_variance')\n",
"\n",
"# 事前分布\n",
"gp = tfd.GaussianProcess(\n",
" kernel=kernel,\n",
" index_points=observation_index_points,\n",
" observation_noise_variance=observation_noise_variance)\n",
"\n",
"# OptimizerはAdam\n",
"optimizer = tf.optimizers.Adam(learning_rate=.05, beta_1=.5, beta_2=.99)\n",
"\n",
"@tf.function\n",
"def optimize():\n",
" with tf.GradientTape() as tape:\n",
" loss = -gp.log_prob(observations)\n",
" grads = tape.gradient(loss, gp.trainable_variables)\n",
" optimizer.apply_gradients(zip(grads, gp.trainable_variables))\n",
" return loss\n",
"\n",
"# -1から1の間を等間隔に100点分割し、indexとして指定する(その点のサンプリングが行える)\n",
"index_points = np.linspace(-1., 1., 100)[..., np.newaxis]\n",
"gprm = tfd.GaussianProcessRegressionModel(\n",
" kernel=kernel,\n",
" index_points=index_points,\n",
" observation_index_points=observation_index_points,\n",
" observations=observations,\n",
" observation_noise_variance=observation_noise_variance)\n",
"\n",
"# 学習\n",
"for i in range(1000):\n",
" neg_log_likelihood_ = optimize()\n",
" if i % 100 == 0:\n",
" print(\"Step {}: NLL = {}\".format(i, neg_log_likelihood_))\n",
"\n",
"print(\"Final NLL = {}\".format(neg_log_likelihood_))\n",
"\n",
"# 指定したindex点に従って、10系列サンプリングを実施\n",
"samples = gprm.sample(10).numpy()\n",
"# ==> 10 independently drawn, joint samples at `index_points`.\n",
"\n",
"# お絵かき\n",
"import matplotlib.pyplot as plt\n",
"plt.scatter(np.squeeze(observation_index_points), observations)\n",
"plt.plot(np.stack([index_points[:, 0]]*10).T, samples.T, c='r', alpha=.2)"
]
},
{
"cell_type": "code",
"execution_count": null,
"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.6"
}
},
"nbformat": 4,
"nbformat_minor": 4
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment