diff --git a/demo/notebooks/prototype_interface.ipynb b/demo/notebooks/prototype_interface.ipynb index 854cd484..0d7a36cd 100644 --- a/demo/notebooks/prototype_interface.ipynb +++ b/demo/notebooks/prototype_interface.ipynb @@ -766,7 +766,7 @@ ], "metadata": { "kernelspec": { - "display_name": "stochtree-dev", + "display_name": "venv", "language": "python", "name": "python3" }, @@ -780,7 +780,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.10.14" + "version": "3.8.17" } }, "nbformat": 4, diff --git a/demo/notebooks/tree_inspection.ipynb b/demo/notebooks/tree_inspection.ipynb new file mode 100644 index 00000000..d87e4eb2 --- /dev/null +++ b/demo/notebooks/tree_inspection.ipynb @@ -0,0 +1,449 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Deeper Dive on Fitted Forests in StochTree\n", + "\n", + "While out of sample evaluation and MCMC diagnostics on parametric BART components (i.e. $\\sigma^2$, the global error variance) are helpful, it's important to be able to inspect the trees in a BART / BCF model (or a custom tree ensemble model). This vignette walks through some of the features `stochtree` provides to query and understand the forests / trees in a model." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Load necessary libraries" + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import pandas as pd\n", + "import seaborn as sns\n", + "import matplotlib.pyplot as plt\n", + "from stochtree import BARTModel\n", + "from sklearn.model_selection import train_test_split" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Demo 1: Supervised Learning" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Generate sample data where feature 1 is the only \"important\" feature." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "# RNG\n", + "random_seed = 1234\n", + "rng = np.random.default_rng(random_seed)\n", + "\n", + "# Generate covariates and basis\n", + "n = 1000\n", + "p_X = 10\n", + "X = rng.uniform(0, 1, (n, p_X))\n", + "\n", + "# Define the outcome mean function\n", + "def outcome_mean(X):\n", + " return np.where(\n", + " (X[:,9] >= 0.0) & (X[:,9] < 0.25), -7.5, \n", + " np.where(\n", + " (X[:,9] >= 0.25) & (X[:,9] < 0.5), -2.5, \n", + " np.where(\n", + " (X[:,9] >= 0.5) & (X[:,9] < 0.75), 2.5, \n", + " 7.5\n", + " )\n", + " )\n", + " )\n", + "\n", + "# Generate outcome\n", + "epsilon = rng.normal(0, 1, n)\n", + "y = outcome_mean(X) + epsilon\n", + "\n", + "# Standardize outcome\n", + "y_bar = np.mean(y)\n", + "y_std = np.std(y)\n", + "resid = (y-y_bar)/y_std" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Test-train split" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "sample_inds = np.arange(n)\n", + "train_inds, test_inds = train_test_split(sample_inds, test_size=0.5)\n", + "X_train = X[train_inds,:]\n", + "X_test = X[test_inds,:]\n", + "y_train = y[train_inds]\n", + "y_test = y[test_inds]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Run BART" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [], + "source": [ + "bart_model = BARTModel()\n", + "param_dict = {\"keep_gfr\": True}\n", + "bart_model.sample(X_train=X_train, y_train=y_train, X_test=X_test, num_gfr=10, num_mcmc=100, params=param_dict)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Inspect the MCMC (BART) samples" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlUAAAGwCAYAAACAZ5AeAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuNSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/xnp5ZAAAACXBIWXMAAA9hAAAPYQGoP6dpAACwrklEQVR4nOzde1zT9fcH8NcYuzAuAzdQVG46zBsoipIMvJdWaprdsBJE7eKF0m5aea/U0i7esouIXaS+pWla2S+1FNDMy7wrOSQhURFkQxjs/vsDNzd2YzhgzPN8PHw8ZPt8PnvD91uczvu8z2Ho9Xo9CCGEEELIHfFq6QUQQgghhHgCCqoIIYQQQlyAgipCCCGEEBegoIoQQgghxAUoqCKEEEIIcQEKqgghhBBCXICCKkIIIYQQF/Bu6QV4Ap1Oh5KSEvj7+4PBYLT0cgghhBDSAHq9Hjdv3kT79u3h5XXneSYKqlygpKQEYWFhLb0MQgghhDRCcXExOnbseMfPoaDKBfz9/QHU/Y8SEBDQwqshhBBCSENUVlYiLCzM+Hv8TlFQ5QKGLb+AgAAKqgghhJBWxlWlO1SoTgghhBDiAhRUEUIIIYS4AAVVhBBCCCEuQEEVIYQQQogLUFBFCCGEEOICFFQRQgghhLgABVWEEEIIIS5AQRUhhBBCiAtQUEUIIYQQ4gIUVBFCCCGEuACNqSGEEELucnKFCmVVKlTWqhHgw4LQlw0+j93Sy2p1KKgihBBCWiFXBUIlshq8vuUkci6UGV8bGC3EsvGxaB/o48olezwKqgghhJBWxlWBkFyhsngOAOy/UIY5W05idUqcx2asqqqqUFVV5dJnUk0VIYQQ0oo4CoTkClWDn1VWpbJ4junzyqoa/qzW5ODBg+jduzfeeOMNlz63VQVV+/fvx+jRo9G+fXswGAxs27bN7H29Xo/58+cjNDQUPj4+GD58OC5cuODwuWvXrkVkZCS4XC4SEhLw999/N9F3QAghhNwZVwZClbVqu+/fdPB+a6PT6TBv3jwkJSWhoKAAu3fvdunzW1VQVV1djV69emHt2rVW33/vvfewatUqrF+/HocOHYKvry9GjBiB2tpam8/87rvvMHv2bCxYsADHjh1Dr169MGLECJSWljbVt0EIIYQ0misDoQAuy+77/g7eb20YDAby8/Oh0+nw9NNP48CBA659vl6v17v0ic2EwWDgxx9/xNixYwHUZanat2+Pl19+Ga+88goAQC6Xo23btsjKysKTTz5p9TkJCQno168f1qxZA6Auig0LC8PMmTMxZ84cq/colUoolUrj15WVlQgLC4NcLkdAQIALv0tCCCHEXEFpFYZ9sM/m+3tmD0LnEL8GPUuuUGFmtgT7b2W+eGwm0pOiEBcWCAAIb8NDiD/Ho+qqysvLsW/fPjzyyCOorKwEn8932e/vVpWpsqewsBBXr17F8OHDja/x+XwkJCTg4MGDVu9RqVQ4evSo2T1eXl4YPny4zXsAYOnSpeDz+cY/YWFhrvtGCCGEEBvkChV0ej02pMYjM60fZgwVgcdmGt8fGC2E0K/hARCfx8ay8bEYGC0Ej83EqpQ4SIoqMHnTEUzedAT3fbgfM7MlKJHVNMW30yIEAgEeeeSRJnm2xwRVV69eBQC0bdvW7PW2bdsa36uvrKwMWq3WqXsAYO7cuZDL5cY/xcXFd7h6QgghxL4SWQ1mZEtw34f7MXnTEaRnHYakqAKrUuLAYzMxMFqI5eNjnc4qtQ/0weqUOPyakYxNeYXIk5abvd+YAnh38PPPP0Otbt6aMGqp0AgcDgccDqell0EIIeQuYevEX560HF4MBn7NSEYgj9XobTo+j11XAF8voDIwFMC3hm3AGzduYPr06fj2228xf/58LFq0qNk+22MyVe3atQMAXLt2zez1a9euGd+rTygUgslkOnUPIYQQ0tzsnfjLuVAGjU5/xwGPJ5wE/P333xETE4Nvv/0WTCYTLFbzFtp7TFAVFRWFdu3aYc+ePcbXKisrcejQIQwYMMDqPWw2G3379jW7R6fTYc+ePTbvIYQQQppbcwQ8rfkkoEKhQEZGBu6//36UlJSgS5cuOHDgAN56661mXUerCqqqqqpw/PhxHD9+HEBdcfrx48dRVFQEBoOBl156CW+//TZ++uknnDp1ChMnTkT79u2NJwQBYNiwYcaTfgAwe/ZsfP7559i0aRPOnTuHF154AdXV1Zg0aVIzf3eEEEKIdc0R8Aj92BgYLbT6nrMF8M1NLpfjm2++AQBMnz4dEokE/fv3b/Z1tKqaqiNHjmDIkCHGr2fPng0ASE1NRVZWFl577TVUV1fj2WefhUwmQ1JSEnbt2gUul2u8p6CgAGVlt1OoTzzxBK5fv4758+fj6tWr6N27N3bt2mVRvE4IIYS0FEPAs9/KFqCrAh7DScA5W06afU5jC+CbU2hoKLKyssBmszFixIgWW0er7VPlTlzd54IQQgipz9q8P7FIgJlDoxHRhodQFw0/Ngxqvlmrhj+XBaFf4wY1twau/v3dqjJVhBBCyN3Kl83EgzGhSEuMhFKjA8fbC5JiGdKzDiM+Ishlw4/5PPcNovR6PU6dOoXY2NiWXopVFFQRQgghrUBZlQpzt56y+l5rannQWFeuXMGUKVOwe/duHDlyBDExMS29JAutqlCdEEIIuVt5QsuDxtq6dStiYmLwyy+/gMFg4NQp68FlS6NMFSGEENIKNOUJQEMdVWWtGgE+LAh93WMLUC6X48UXX8SmTZsAAL1798bXX3+NHj16tPDKrKOgihBCCGkFmuoEoLUC+IHRQiwbH4v2Lip+b6ydO3di06ZN8PLywuuvv46FCxeCzW75YM8WOv3nAnT6jxBCSHMokdVgwfbTuCc0AHFhgVBqdAjisRDehocOQTynnydXqDBjswQ5UuuBmquK3xtLr9fjpZdewuOPPw6xWOzy57v69zcFVS5AQRUhhJDm8t8NBeZuPWk2p6+xmaX8q5UY8VGOzff3zB6EziF+Dp/jrtuHjlBLBUIIIeQuJVeoMPfHUxaDj/dfKMOcLSedyizJFSr8V1Fj95qGFL+7YvtQq9WioqICQqH1ju6tBZ3+I4QQQloJe4OVDW0VnHmWI46K3+UKlUVAZVjLnC0nIVc4/ozCwkIMGTIEo0ePhkajcXi9O6OgihBCCGklXNlWobJWDUmxDGKRwOr7ydFCaPV6FFyvshocyRUqXJHXNjjIkytUKCitgqSoAgXXqyCrVmLjxo2IjY1FTk4OTp8+7batEhqKtv8IIYSQVsKVbRUCuCxk5hZiVUocACDPZEtRLBJg3qjuGLs2DwqV1mI774qsBn/+cx1RAl+7n2EI8upvEWqrZUDuZ/jv+H4AQFJSEjZt2oROnTo1eP3uiIIqQgghpJVwtq1C/QJyP443qpUayGvUaOPLRnxEEDKyJUhPikK6OMo4/uZaZS1+PX0FCpUWgHnNFgBcuqHAzpMlSBdH2V2vP5dldYvwxp7PoDi3H15Mb8xfsAhvvfE6mEzmnf54WhwFVYQQQkgrweexsWx8LOZsOWkWWA2MFmL5+FizInVrBeRJIgHSxFHIyJYAADLT+mHtH1Ks2Ss1XiMWCTBJHIW3fz5n9tmG7TxvLwZW772APGk54sKDIBYJzLJcpmsS+rGt1oEFDU6H9mY52gx/Dk8/N8kjAiqAgipCCCGkxTnTkqB9oA9Wp8ShrEqFm7Vq+HNZEPqZX2+rgDxXWg49gPSkKGTmFuKvi+V4+f57kDEsGn5sJjjeTOw4VYKMbIkxS2XqZq0abG8vYxBla/vQNMi7WFZt8RzvACHaPbXc+ExPQUEVIYQQ0oIa05KAz7PfB8reKcE8aTmmJnVCnwmByMwtxEe7LxjfS44W4pUR9+Crg5egUGnBYzORnhRlbDTKZddllHhsJhQqLRQqrcX2YYSAh46BPuDz2FCpVPBj289C3cl4HXdDQRUhhBDSQhy1JGhsR3NHpwSDfNlY8dt55Nbbtsu5UAbogaxJ/fHC10exekIfrPztvNn2YHK0EKtS4ozZLIVKa/b+77MGgs9j4+zZs3j66afxeMpTGBid5PLxOu6IWioQQgghLcSVfadMOTolyGJ6WTQQNciRlkGr0+PrKQl4/7fzFtflXChDVl4h0pMsi9STo4UQ+rLw8ccfo0+fPpBIJFjz8YdY+FAXDIw2b+xprQ6staNMFSGEENJCXNl3ypS9U4JikQBVSvvPVWq0uFyhsVqADtTVZk0bLDLLUA2MFiIjoQ0eGzsKe/bsAQCMHDkSmZmZCG0X5LAOzBNQUEUIIYS0EFf2nTKQK1Qor1ZhwZgeWPjTGaun/7gs+3VOvmxvXLqhsHsNl8XEntmDjEGSP0uPPjHdcPnyZfB4PKxcuRLPPfccGAwGAMd1YJ6AgipCCCGkhTjbd8oR06J3HpuJZwd2wqsj7gHb2wt6PVCt1MCf6w0208tmKwSxSAA2ywscb/sVQnwflsWw5fnz52PDhg346quv0KVLF6fW7gkoqCKEEEJaiDN9pxyxVvTeq2MgqpUaLN8lNQug3n80FjOHRgOw7KQ+Y0g0NFqdcYSNtcAr2UbA9/hTqRg06nFUa+rG29hrDeGJKKgihBBCWlBD+k41RP2i9/SkKFyR1+DnU1csAqMFP51BVlo8RsW2N+ukXlpZi5AANvQACkpvYnJSFLwAs2L1ZJEA747tabG+ElkNXv/hJHKk5sHhu+NioNLqIK9x3IOrtaOgihBCCGlhrqg3ql/0HhcWCABWM00KlRbTNkvw3bMD6gIehRp+3LqQ4IlP/8I97fyxYHQPyBVKPBATijSTwOvAX4cw+r6B+HHLD4iKqjsBKFeoLAIq4FZriK0n0Ts8yFjUbtqDy5mmp60BBVWEEEKIB6hf9K7U6MD0YmDGUNHt5p0sJo4VVeDbv4uwbHwsFv502iwLJRYJ8EVqP8iqVSi+oYDQj4MSeS3e/vkcqmtVkB/8H+R52YBehxdnv4IPP92Eylo1fNhMi4DKIFdajkkmMwL3XyjD61tOYukjMVj40xnsPldqfM9R01N3R0EVIYQQ4gHqF73zWEy0D/TBlwf/tZjtl5nWD8t3nbfIYuVJy8HAebPMklgkwGsDAvDClHTUluTXPbtrMqa9uQyj1+QiPSkKA6OD7a5NqdGZfZ1zoQyXyhV4sn84DhSUWx3c3BozVtT8kxBCCGlhcoUKBaVVkBRVoOB6FeQK55t+GoreDU02tXo9luw8YzVwullrvweVYesQAP7IO4wpY4eitiQfDI4vhKNfRfDDr0PP9cOqlDhIiirAZdkPJ6ydJJTXqLHRShPRO2l62tIoU0UIIYS0oMbM/jOwVpNkKHpXarQ2u6bLaxw1/7ydWWIFR8C7bTRiOwShNC4d3gF1WalgPw4+3P2PMTiz16JBUiyzeJ1zazBzutiyM3trHbJMmSpCCCGkhTia/WcvY3W5QoEzJZXIv3YTFQo1fpRcxivfn0C1SovOIX4WW26mHPWgMn2f4cVEyCNvYfH6bGNAJRYJwL4VFAHAH/mlmDFEBLFIYPacJJEAk8RRyMwtNHvdNNCyts7WOmSZMlWEEEJIC2nI7D9rtUX/3VDg9a0nLXpMTRJHYcH201jxWC+73dolxTIkRwutfra1zJIX1w8+7LqQITlaiMlJUbgirzW+/+m+i4jtEIiHYkKNLRq4LCZC+Vys/L98Y82U6TozsiUALAO81jxkmYIqQgghpIU0ZvafXKHC3HoBFXC7dUJceBCuVtaiXQDXZrf28yVyLBrTAwu23z79V3vpJPz+O4i0iZ/hxW+Pm12fLBKijS8bG1LjUVpZCz8OE2182chM6wdvLwba+LKh0ekQyGOB6cXAFXkt2vO5COVz8c64GKSLqyCrUYPj7QVJsQwZ2RIoVFok1QvgWvuQZQqqCCGEkBZiyCbx2EykJ0VZtD4I8LHMNpVVqWzWShlqlP6rqEGQjW7tySIhJid3xuOfHsST/cMxIT4UX61ZjuxvP8E1AItXroUi+F7j9WKRAKniSDz1xSEoVFrw2Ez8PDMJ87efrtcUVIiXR9yDT/6UIjUxCqF8bl3/LQBand5iHWKRAK+O6IqKahU2pMajY5AP2gVwW21ABXhYUBUZGYlLly5ZvD5t2jSsXbvW4vWsrCxMmjTJ7DUOh4Pa2lqLawkhhBBXE/qxcV+3EDzRPxwb8wrNWh8kiQR4tE8HXCqvhkqjQ5VSgwAfFuQ19k/GGZp0SkurECX0NRauVyhU8GEzodXpUVGtxpP9w7Huh93478f3oC4rAgAExo3E6HGPYHz/aFSrNFAotWB5eyHnwnXj89OToiwCKgDIkZZBDz1eG9kVq/ZcwMrHehnfM3SNL72pRNGtQc2SYhlSPv/LuDU4MFqI1Slxd/YDbWEeFVQdPnwYWu3tfdvTp0/jvvvuw2OPPWbznoCAAOTn5xu/NkzTJoQQcvdp7g7ffB4bC8f0wGtbLLfzjhXJcEVei7V/SJFr8t7mKQn2n+nDwsGL5egeGoDXb/V84rGZ+LuwCiEBHGMm7OKBX1Gc9Ro0ajWYvoFY+P4q/OvTFX06t7faFHRVShzmbDmJwV2CzYI/U7nScrxQq0HX0ACLejA+j42yKhUmbzpi9V57NWSthUcFVcHB5s3Hli1bhs6dO2PQoEE272EwGGjXrl1TL40QQoibu5PWBrY0JEirVeustiJIT4rCmj+kFu8duFiOJJHALNAySBIJEOLPQXx4ELgsJtLEUbhZo8Z/8lrsPFVifBaPzUR6z1h4M70xYNBwvLRoBYoU3ng9JhSLd5yxyELlScvhBeC7ZwegWqUBj800Kz43+55r1IgLC7RaD9aYGrLWxGNbKqhUKnz99ddIT0+3m32qqqpCREQEwsLC8PDDD+PMmTMOn61UKlFZWWn2hxBCSOt1J60NbCmR1WBGtgTDPtiHcesOYNjKfZiZLUGJrMbsOluBRlxYoNVgKzO3EGniKCTfavJpkCwSYvHDPfHB7/lIyzqMJz//C+lZh1Gp1GD13gtmAdWqlDhIKnlo8/SHuNRnGmb/VIij/96AVqe3Wa+VIy3Hv+XVeH9XPlbdyn5Zw/H2gi/HG2qdHhevV+FaZa2xsakPm4kZQ0U2722trRQMPCpTZWrbtm2QyWRIS0uzec0999yDzMxMxMbGQi6XY8WKFUhMTMSZM2fQsWNHm/ctXboUixYtaoJVE0IIaQmNbW1gi6MgzXQMi63WB7b6TClUWmRkS7B5agJeG3kPSiuVxvfe/vks9p6/bna9Rqc3C87Sk6KwMa8QedJysAS3f9flSMtxuV7AZ21Nhtqp9KQoi21AQzuGuLBAzMyWYFVKnPGzDJJubSUaTgAatOZWCgYem6nasGEDHnjgAbRv397mNQMGDMDEiRPRu3dvDBo0CFu3bkVwcDA+/fRTu8+eO3cu5HK58U9xcbGrl08IIaQZuXpbqiFBmoFhZl999hp0KlRalFep8MSnfxlbEgT7cywCKgD47eedqD6fa/zaVgasIUL8OeCxmciVlmNAJ/NGn8m3+k+dLZFDUiwzC95M5UrLkVVvPE1rb6Vg4JGZqkuXLmH37t3YunWrU/exWCzExcVBKrVegGfA4XDA4XDuZImEEELciL1GmYDz21LOBGl8G60PrlXWIlkkRI7UdoNOhUprzBate6qP2TU6pQIVe7/Aayf/Dww2D5z298A7INhqBszQ0iHYn2P3M/+5dtOYZdLq9diQGm88bRjsz8FHu//BhIQIZGRLsDolzm5B+1sPdcfwriHw57Ig9GvaAwHNxSODqo0bNyIkJAQPPfSQU/dptVqcOnUKDz74YBOtjBBCiDvy43pj85QEyGrUxh5RmbmFUKi0jdqWcjZIM7QcKKtSQVajglKtw+FLNzApKRI66K12Tp+z5SRmDBUZe1uFteEZr6n97yzKd66ERn4NDAYD3YaMQxUvEMDtDJghkOoTHgQeiwkd9Nh7vhTpSVE2PzMjW4K48ECkJ0VBrdEZT/IlRwuxaHQPdG/PN27r2RuTAwC1ai16hwc5/mG2Ih4XVOl0OmzcuBGpqanw9jb/9iZOnIgOHTpg6dKlAIDFixfj3nvvhUgkgkwmw/vvv49Lly5hypQpLbF0QgghLcDaqT9DC4Hv/i7C4od7Op1FMWzpWetmbitI4/NuZ2vkChXaBnAhr1Hh9ZFdodLocL1KiWA/Dv65dhNztpzEsvGxZr2tZgwVIUkkwG9frUF5bjb0Oh2CQzvgjWVrMPK+oVjxWz52ny+FpFiGoV2DMSEhwqI3llgkQGxHPvpHtTEOOg7ksSBTqMFkMLDisV7gsphoF8DB7+euGb+fd8fFoFqlQffQAKyZ0AfHiirAZVkvRjewlf1r7rYWruRxQdXu3btRVFSE9PR0i/eKiorg5XV7j7qiogJTp07F1atXERQUhL59++LAgQPo3r17cy6ZEEJIC7FVUJ4nLYfXrSCibQDX6WeWV6uwYEwPLPzpjEWLhobUDhkCLElRBcasyTO+bji5N/u+Lhb1Spm5hVg7oQ/8z0di/X4dfHsOA3f4s/joLAtH1ecwb1R3vDyiC6qVGoyP64C3tp+2O+pm8qYj4LGZ2DZNjA9+/8fs2uRoIRaN6YHRse3B8fbCgp/OYPe5UuP7YpEA93dvC7FIYLV+y1ZgaRrgGjJpiZ3qhjcH+bLdPsBi6PV6fUsvorWrrKwEn8+HXC5HQEBASy+HEEJIAxWUVmHYB/tsvr9n9iB0DvFr8PNsBQUclhcCfdhO1w5ZWx+PzcSX6f3x6PqDFte/NDwaf1+8jj17/4RPZG+z9+paLvSARqdHcUUN0rMO2/zcDanxmLzpCGYMFeFEUYXVNgsDo4V4/7FeeOX7E1aL8od2DUbGsC5Y8dt5s55ahsAytF7vL7lChRnZEuPPztrJwTvtG1afq39/e1ymihBCCGkoV576q5/1MhSRr9krNY5gccU2okKlRelNpdXrYzrw8dHuCxYBFVA3RqZKpYFeDzjKpxjqoeLCAm0Wm++/UIaKatunHPeev47XR3bFysd7o6pWg5u1artF6aYnJm2dHLTWksKdUFBFCCHkruXKU3+u7nUF2D4ZyOd6o/rsn+BG9ALT93axt73icB6bCY1WjwqFGu0cbGnaa+dgqrJWjU2T+kOj04HBYKBWrTUr9K9RaXFPOy7aNiAJZBrgOgrm3HWcDQVVhBBC7lqNKSi3palGsJieDLxZq4ZGUYlFrz6Hsh0/wkfUH8GPzDNODrEVDBm20z74v3zkSMsxY6jIZr2ToV2DWCRAKN/+NluVUgudXm+RVTIU+gf4NDwoNQ1wHZ0cdNdxNh7b/JMQQghxxJAJqt98szHNKF3d68oUn8dG5xA/XD17CGOGDsDO7T/C29sbDw1NwufP9MG6p/ogM60ffDneGN4txOJ+w3aaoTYqM7cQk8RREIvqN/AUYv6o7hB3FuD1kV1RerMWSfWuMRCLBNBbCaiAuoL3rLxChycATZk2QXWUKXPXcTaUqSKEEHJXq58JamwzSldmveqrqalBxkuz8MVndRM/Okd3war1G7C1iINnv5aYfc6SsT0xontbtPHjQKnRGVsgZOYWGq8zjLpJT4pCujgKvhxvVCs1kBTLMG7dAShUWiSLhHjjoa5IE0dBD1jtW8VkMGx2Z8+VlqOyVo32aFhRuelWpyFT5szJQXdAp/9cgE7/EUIIAepO/9Wvf7J12s0Z/5bK0Se+PyqK/4F/39EIHJSGgd3aI+1WQ07TGXrJ0UI82DMUc388dfs1kRCp4kiLaw3WPdUH0745ZvH6988NQOrGv5GeFIW4sEDw2N5QqOqCr8zcQqx4rJfV+wz+99wA9I9qA8B2/ynT1/k+LGOApwdstqS4k5+lKTr9RwghhLgpV2W9TMkVKszbcR7c+15ESNUN+ETFAajLBOkBi8HGORfKkJYYafaMHGkZdDaGIAN2arE4TLNROIZWC47uMwjg1oUZ1hqsDowW4u2xPbF451mzHlembRPWuPhn2dQoqCKEEEJcyLQzuisYThWygyOA4Aiz9/Kk5cbO56asFXrbutZQmG4NAzCbBVh/W87wtaRIZsxmGbYcr8lr0Ma3LhM1f/tp9AoLRFpipPH9Y0UVWLTjDLq155sFVfXbJrhzEFUfBVWEEEKIm9Dr9fi///s/jBgxwviao1OF1gKohrZESBYJMG1INCZvsmwEmhwtxO/nruHlEfcAqMt2ZeYWYlVKHBioy5QZurhzWV5Y84fULAuWHC3EoHtCUF2twpP9w62OxJkkjoL3rZOLpty5bYI9FFQRQgghbqCkpATp6en47bff8O233+KJJ54A4PhUYf0Ayl7mKSSAg2+mJEBeowbH2wunLstRq9YiLjzQonP5gjE9MHp1Lr46eAlfpSegRF4DpUYHbwYDk8RRSE/qhFq1FsF+HCzfdc6iqDznQhle33ISix/uYfOEIAC8cv89Vtfqrm0T7KGgihBCCGlh33//PZ5//nncuHEDXC4XlZWVxvfsnSpMqhdADYwWYtoQkdURNGKRAHvOleLeTgJ8+3eR8Xk8NhPzRnXH/FHdUaPSGmuXrlbWQqHSQqHS4kplrVktlakNqfFWx9gAdYFVrVpn84RgnrQccx6wzFQB7ts2wR4KqgghhJAWIpPJMHPmTHz99dcAgD59+uDrr79Gt27djNfY6qo+MFqId8fFQKXVYXjXEGMwpFBpER8RZHZtskiIV0bcg8/3F+DJ+LAGFdNfu6k01kvp9Xqz2ipnOMo4VSgs33fntgn2UFBFCCGEtBCpVIrs7Gx4eXnhjTfewLx588BmWwYTzpwq5POApY/E4FK5ArJb23ySYhk+21+AuQ90Q/WtzFOADwtRQl+bdUveXgxMSeoELssLn+dcRKo4EjroLbYJOwTZb2/gKOPk7WWeqWpM41V3QUEVIYQQ0kLi4+OxZs0a9OrVCwMGDLB7bUNPwskVKszZesqsfYFhTM3rW09aBEWG9gX1CXzZOHqpAtuPX0autBx/XbxhbBaq1OgQ6MNC5xA/1Kq1dkfe+LGZSI4WWp2LKBYJUKPSYteLyahVa1tF2wR7qPmnC1DzT0IIIe6ioLQKwz7YZ/bajKEiSIoqbHYoN7QvqC//aiVGfJRj87P2zB4EoR8b567exOq9Fyy6rs8cGo1u7fxxs1aDOVtPItdKV/bv/i7Cisd6tUggRc0/CSGEkFZGo9Hg4sWL6NKlS5N/lrUWDHFhgVabfgL22xdY675u6matGp1D/BDRhoexvTtg1vAu0Oj04LGZYHkxwPdhGTNsy8fHWmxJfvd3ERY/3LPVZqbqo6CKEEIIaUIFBQWYOHEiCgoKcOrUKQQHBzfp51lrwWCtl5UpW8XkDR0SHRrogwGdBJi79aTZSUDT7cUOQTz4cbyNdWHjeneAMCnKYwIqgIIqQgghpEno9Xps2LABL730Eqqrq+Hv749Tp05h6NChLnm+rVl69Vsw8NhMdAzywYbUeLNu5pm5hcZMlLVicrlCBZ1ejw2p8WAwGBb3mJ7QkytUmPvjKYvWCta6o3tSEFUfBVWEEEKIi127dg1Tp07Fjh07AAADBw7Epk2bEBkZ6ZLn25qlZ8gKGVowHLlUgVUpcXhv13mLeqZVKXHIyJYgPiLIon2BtefXv8f0hJ5hlI41rbU7emNQUEUIIYS42Nq1a7Fjxw6w2Wy88847mDVrFphMpkueLVeoLAIewDwrZGjBIFOo8da2U2YBFXC7m/m8Ud0xuEuwWcBj6/l50nJ4MRj4NSMZgTyW2T2ORum0xu7ojUFBFSGEEOJib7zxBi5cuIC5c+ciNjbWpc9uaFaIz2PXXWunm/nC0T0QWq+dgr3n51wog0ant8g6BXBZ4LGZFkOVDVuGrbE7emNQUEUIIYS4GJfLRXZ2dpM825mskKNrq5UaAOb1WRqd/U5L1rJOQj82MtP6YfXeCxZDkzPT+rXK7uiNQUEVIYQQ0khKpRJqtRp+fn7N9pkNPZHX0Gvr109tSI1v8PNNrd0rtTo02YvBwJqUOLvP9BReji8hhBBCSH2nT59GQkICpk+f3qyfazjdZ039mXmOrvXjelvUT0mKZRCLBFbvSRIJwGVZhg5124y2twzLqlQ2vx9PQkEVIYSQRpMrVCgorYKkqAIF16sgV3j+L0+dToeVK1eib9++OHHiBH755RdcvXq12T7fMGC5frBkbWaeo2urlRqL+qnM3EK8OqIrkuoFVmKRAGniKCz86YzF/85UqF6Htv8IIYQ0iqNj/Z7o0qVLSEtLw59//gkAGDVqFL744gu0bdu2WdfhzIBle9dKiiosrleotJBVq9A7PAiTbs35M3RAz8iWQKHSWrRIcGZL0pNRUEUIIcRpDTnW72l9ifR6PcaPH4+jR4/C19cXH374IaZMmQIGg9Ei63Gmkaata20FQwq11uZYG8Ay81S/4aip+luSnoy2/wghhDitIcf6PQ2DwcDq1auRlJSE48ePY+rUqS0WULmKrZorjrf98KB+5smZLUlPRpkqQgghTrtba2gGDBiA/fv3t/pgysAQDM3ZctIsy1R6U4nkaKHVwNlW5smZLUlPRUEVIYQQp3l6DY1arQaLZf178JSAysBWMDSoS7BFsOUo8+Tps/0coaCKEEKI0zy5hubQoUN45plnsGbNGtx///0tvZxmYS0Y4vPgMPNka6jz3Yqh1+vtt04lDlVWVoLP50MulyMgIKCll0MIIc2iRFZjM5NRf/RJa6BWq/H222/jnXfegVarRUJCAg4ePOhxmSlX8YTTn67+/e1RheoLFy4Eg8Ew+9O1a1e793z//ffo2rUruFwuYmJi8MsvvzTTagkhpHUzbBvtmT0I26YlYs/sQVidEtcqA6rz588jMTERixcvhlarxVNPPYVdu3ZRQGWDo9Ofd0O/Mms8bvuvR48e2L17t/Frb2/b3+KBAweQkpKCpUuXYtSoUdi8eTPGjh2LY8eOoWfPns2xXEIIadU8oYbmypUr6Nu3LxQKBYKCgvDJJ5/giSeeaOllubWGDnW+23hcUOXt7Y127do16NqPP/4YI0eOxKuvvgoAWLJkCX7//XesWbMG69evt3mfUqmEUqk0fl1ZWXlniyaEENJiQkNDMWXKFJw7dw4bN25Ehw4dWnpJbu9uPf3piEdt/wHAhQsX0L59e3Tq1AlPPfUUioqKbF578OBBDB8+3Oy1ESNG4ODBg3Y/Y+nSpeDz+cY/YWFhLlk7IYSQlvH+++9j165dFFA1kKef/mwsjwqqEhISkJWVhV27duGTTz5BYWEhkpOTcfPmTavXX7161WK0QNu2bR3OcJo7dy7kcrnxT3Fxscu+B0IIccbdOHuvKbDZbHh5edSvxCblzFDnu4lHbf898MADxr/HxsYiISEBERER+N///ofJkye77HM4HA44HI7LnkcIIY3hCaevmtOePXuwePFi7Ny5E/7+/i29HLfQ2JYItpqG3m0d1OvzqKCqvsDAQHTp0gVSqfX5Re3atcO1a9fMXrt27VqDa7IIIaSl3I2z9xqrpqYGc+fOxccffwygroTj3XffbeFVtbw7Dcqpg7olj851VlVVoaCgAKGhoVbfHzBgAPbs2WP22u+//44BAwY0x/IIIaTR7sbZe41x7Ngx9O3b1xhQvfDCC3jzzTdbeFUtz1UtEfg8NjqH+KF3eBA6h/jd1QEV4GGZqldeeQWjR49GREQESkpKsGDBAjCZTKSkpAAAJk6ciA4dOmDp0qUAgBdffBGDBg3CypUr8dBDD+Hbb7/FkSNH8Nlnn7Xkt0EIIQ411+mr1twx++eff8bYsWOh0WjQrl07ZGZmmpWJ3M2oJULT8Kig6r///kNKSgrKy8sRHByMpKQk/PXXXwgODgYAFBUVmRUiJiYmYvPmzXjrrbfwxhtvIDo6Gtu2baMeVYQQt9ccp69ae81WcnIyOnbsiL59+2L9+vUQCq0XVt+NqCVC06AxNS5AY2oIIc1NrlBhZrbE5uy9O62pkitUmJEtsZrNcMXzm0tpaSmCg4OpM3o9BaVVGPbBPpvv75k9CJ1D/JpxRS2DxtQQQggxnr6qf6zdVaevPKVmKyQkhAIqK6glQtPwqO0/Qgi5m9zp6St79VKtaXvoxx9/xKlTpzB//vyWXkqrQS0RmgYFVYQQ0oo1dvaeo3qp1tAxu7KyEi+99BI2btwIABg2bBjEYnELr6r1oJYIrkfbf4QQcpdpyHF6d98e2r9/P3r16oWNGzeCwWBgzpw5iI+Pb9E1tUbUEsG1KKgihJC7jKN6qauVtU1es3Unli9fjsGDB+Pff/9FZGQk9u/fj6VLl9KkC9LiaPuPEELuMo7qpf6rqEG7AK7bbg/16tULer0e6enp+PDDD+nUNXEbFFQRQshdxlG9FABj88fG1mw1pZEjR+LEiROIjY1t6aUQYoa2/wghxMPIFSoUlFZBUlSBgutVFiNHhH5sJNuolxKLBJAUy9zqdJ81FFARd0RBFSGEeJASWQ1mZEsw7IN9GLfuAIat3IeZ2RKUyGqM1/B5bCx5uCfEIoHZvWKRAJPEUcjMLWzx0316vR5ZWVk4dOhQi66DEGfQ9h8hhHgIR6f6TLug89hMjIptj3RxFJQaHTjeXpAUy5CRLUF8RFCLnu67fv06nnvuOfz444+Ijo6GRCKBr69vi62HkIaiTBUhhHiIhnZBL5HV4I0fTyHYn4NvDl3C2SuVAIDuoQHYkBqPpY/EtFgd1c8//4yYmBj8+OOP8Pb2xqRJk8DlcltkLYQ4izJVhBDiIRrSBd00m3W8WIYvUvth5W/nsWav1HhdcrQQC8f0AAOAL8cb1UoN5DWWXdftdWR3llarxfTp0/Hpp58CALp3746vv/4acXFxjXoeIS2BgipCCGnlDMGNRqdHZlo/HCuqQGZuIRQqrdl1/lyWWTbryf7heP+388iTlptdl3OhDPO3n0ZceBCOF1UgTRyFjGwJFCqtses6A8Brt4IzHpuJ9KQoJHYSgMvygsCPA5VGhyqlpsHBFpPJRHV1NQBg1qxZePfddylDRVod2v4jhJBW7HKFAmdKKpF/7SZu1mpwrKgCZ0vkWJUSBx6babzO0AW9slYNHpuJGUNFGNGjrUVAZZAnLUdcWCBypeXYmFeI9KQoALfrs/7857oxoFqVEgdJUQWmfHkENxRqvPHjKdz34X6bhfK2rF69Gn/88Qc++OADCqhIq0SZKkIIaaX+u6HA61tPmgVGySIBFo7piRvVSnz/3ADo9MD5q5VIFgnrek75qLAqJQ4b8wrRPdR+00ylRgegLsBKF0cZX99/oQypiZEAgPSkKGzMK0SetBwzhoqMfzdlrVDemsDAQAwePNjJnwIh7oOCKkIIcWNyhQqlN5WQ1ajhy2bCl+ONQJ+6dgdz6wVUAJAjLcf8n+q27tbslUIsEmDm0Gjob73vy/E2Bj6mgZI1HO/bmxmGAKv+13FhgcZ6LNO/G7YE48ICodTowGUxIVOo4c/1xldffYUnn3ySxsoQj0NBFSGEuKkSWQ1e/+EkcqS3T/QZgqQgHgs5drbuDAGTpEiGi9erwGUxcUVeAx7HG3HhQZAUySAplkEsEljdAjQ0ATUwDbBMvzYNtgx/N2wJbswrNCuA7yPQoeyXj5Dz516cPXsWy5cvd/In4llcWehP3AMFVYQQ4obkCpVFQAXAGABNGyyye79SozMLbt748bTxPbFIgFUpcZiz5SSWjY81e67h/Um3itMNX5sGWAOjhSi9qQRgHmwZ/m66JWhQfXYftv/fOuiU1fDx8UFkZGRDfxQeqURWY9FTzHAIoH2gTwuujNwJCqoIIcQNlVWpLAIqgzxpOeY80NXu/RxvL6vBjeF+oO70X0a2BOlJUZic1Ak8NhM6nR4HLpYbT/sliQTG039A3S/+5bcCsYHRQrNsl+HvptuA2toq3Pi/T6A4tw8AwA6NxvffbkbXbl0hV6juysyMM01aSetCQRUhhLghRz2nNFo9kqOFVpt9GjJLpsFNfYYtQoVKa7xmaNdgPJUQgbiwQKyb0AeRQh5UGj3ktWpseX4AeGxvBPJYxl/4q1PiUF6twri4Dlj40xlk5hZiVUq9vlJaDWovnQAYXuAnPgH+gCfACOqAYSv33bWZmYY0aaWgqnWioIoQQtxQgIPZe3KFGkse7on5209j/wXzmivD1t2Kx3rZfYZpPZRYJMCEhAjMzJZA3FmAeaO6481tpy22p94dF4PyapWxGajfreags4ZHg8fxhrcXAxqd3ngP0zcQwtGvwIvtA077ewDc3ia8WzMzDWnSSlonCqoIIcQNCf3YGBgtNAuYgLoi8Lce6oZAXzZkCiWWjO0JlUaHaqUGPmwmjhXJjFt39YvL6+sk9MW2aYnw5XiDzfSCvEaFHTOS4Mf1xivfn7C+PbX1JHrfOlkIwLg9OGfLSTzZPxyJnQQQ+LHNsmg+kb2NzzCtz+KxmYgNC8QVeS0ullXfNcXajgLmlh5mTRqPgipCCHFDfB4by8bHmtXe8NhMbEiNx7o/pGaF56bbaO0CuIiPCIK8Rg2BLxvJIoHVU4LJorru5zVqLby8GAjksRAhrBtaXFBaZXN7KldajkkmrRhypeXQA/js6Ti8/skWrNnb0Vggr9frkWujAN7WCcHGbgm2ppN0tgJm4HaTVtI6MfR6vd7xZcSeyspK8Pl8yOVyBATYb6ZHCCHOMPSpulmrhsCXg8uyGshq1OCymGbjaO7rFoK3x8WgqlZjDCxYXgzIatR4f1e+WdF7skiIaUNEmLzpsHGUjWkwIymqwLh1B2yuad1TfTDtm2PGr9U3LiPw789w9tRxvLjqO4wemgSNTo/2fB9o9XqotTpU1qghKZYZ1ztjqAiSogqr7RwGRgud2hJsjSfpSmQ1mLPlpFlgZTgEEOqma/ZErv79TZkqQghxE7ayLXwe+1bPqhNmWSfT1ghP9A/HK/87bvZ+crQQM4aIEB8VhDRxJJQaHQJ9WAjisfBM5t9mswFN65scbU8ZthX1ej2qTuxCxd4vUKJWwtc/AJLzF7Gt+PaIGbFIgHfHxuCzfQX4/Vyp8XV7RfTOFGu31pN07QN9sDolDmVVKtysVcOfy4LQz32za6RhKKgihBA3YC/b4stm3upZZb01wvLxsci00joh50IZdHo94sKDMHnTEePrYpEAT/YPtwhq9l8oQ+lNJdjeXhjeLQRdQwPMOqIb5gpKimXQVlWg/NePUXOx7rkJ4mS0Hz0bxyrMA7I8aTne2n4ayx6JgVKjM2Zm6ndor6+hxdqt+SSdIWAmnuOOgqra2loaekkIIXfIUbZl8cM9HPassjcYuf44GkmRDK+P7GoRMGXmFuJyRQ2YXgzMeaArFvx0xizwShIJMH90D0z4/C/oSk6j5uIRsNgcTJ71BmbMyMCvZ6/h/K3tPVM5F8pQq9aZZWa4LCbsaWixNp2kI+7E6aBKp9PhnXfewfr163Ht2jX8888/6NSpE+bNm4fIyEhMnjy5KdZJCCEey1G2pbpekFJfVa319w3z94L9OVj3VB9wWUyc/E+GXh0D8f6u81a3EoP9OTj5nwzr9xdYBGq50nIs2XEGT98bgR6PvITZ6lKUh8TjV0Tg1zV5xmcYTh+aulmrRucQP2NmRq5QuaRYm07SEXdi/7ytFW+//TaysrLw3nvvgc2+/X/6nj174osvvnDp4ggh5G7gKNuiUNoPqgK4lv99bDhdJymqwJg1eZj2zTGkZx3G0X9vgMvywtEimdn1edJyZOUVgsPyQkgA12bmK0dajpE92uHLA//iZo/xYAdHmD1jY14h0pMsBzXXD24MpxsHRgvNXjcUazd0W8xwks4aOklHmpvTmaovv/wSn332GYYNG4bnn3/e+HqvXr1w/vx5ly6OEELuBo6yLSyml53WCEJo9XqLwci2RtTkSMuhu/V+/ZqqXGk5lGqdWb2TXqOGViGHd8DtwEVWo27QMGcDW8GNK4q1DcGZrZN0VLNEmpPTQdXly5chElkO8tTpdFCrae+aEEKcZa9vUZJIgD//KcWkpCgADORIy4zbegM6CcBkMFBRrcLrI7vio93/YO/56wDsn66zFvgYVJs0DVWVFqJs50owvJho98wKMJh1wR+Pbb8eyjQoS3YQ3LiiWJtO0hF34XRQ1b17d+Tk5CAiIsLs9R9++AFxcXE27moeS5cuxdatW3H+/Hn4+PggMTERy5cvxz333GPznqysLEyaNMnsNQ6Hg9ra2qZeLiGEALCdbUmOFiI1MRJztpzE0/dG4KX7ojFzmAgBPixILlVg6pdHjLVLYpEAr47oion3RkKr1zusJbJ1+s6P441/rsgQJP0Vx3/8FNBp4MXjQ32jBOzgCCRHC8Fj2//VEdbGB+ue6oNAHxY6h/ihbUDTH2iik3TEHTgdVM2fPx+pqam4fPkydDodtm7divz8fHz55ZfYuXNnU6yxwfbt24fp06ejX79+0Gg0eOONN3D//ffj7Nmz8PX1tXlfQEAA8vPzjV8zGIzmWC4hhBgZsi0yhRrVKg2qVVr4c72x59w1vPdoLDbkFuKj3ReM19cvCq/b5juPUbHtEezPgUKlsft51kbYJEcLceW/S1j50rM4fvgvAICPKAGCkTPB9A2EWCRAamIkfj5VgiSRwKxbuum6fjtzDSeLZVg+PrZZAipC3IXTQdXDDz+MHTt2YPHixfD19cX8+fPRp08f7NixA/fdd19TrLHBdu3aZfZ1VlYWQkJCcPToUQwcONDmfQwGA+3atWvq5RFCiF3VKi3e2m4+xPjdcT2xMdeyNsrwtWltVJ60HPMe6o5HPjmA9KQoizorg2SREKWV5tl4Q8A0NfVpFB7+C1yeL+YsWoZHU57G1UolAEBSXDdXEABWpdTtTJgGVskiId4a1Q1XZLV4rE9H6gxO7jqN6lOVnJyM33//3dVrcTm5XA4AaNOmjd3rqqqqEBERAZ1Ohz59+uDdd99Fjx49bF6vVCqhVCqNX1dWVrpmwYSQu5atXlVtA7h2i8JnDe9i1m/KkGjPzC00Bj559bqsLx7TA3oAm6ckQFajBsfbyxgwqfunoqNaia8/Wwthh3Do9cDkTUeMdVyrU+Kg1OjA8vLCJHEU5jzQDUU3FMZnjFt3AAqVFtumJSICtncICPFEd9T8s6qqCjqd+b68u8y+0+l0eOmllyAWi9GzZ0+b191zzz3IzMxEbGws5HI5VqxYgcTERJw5cwYdO3a0es/SpUuxaNGiplo6IeQuZKtXlaPO4/IatVm39ORoITLT+iE96zAysiVIT4pCujgKKq0O4UE8nPxPhodW52J1SpzZfQZMvyBgxOvQ+rfFQ6tysSE1HkI/NjLT+uFmrQbyW3MHD1+6gbMlcrw+site+f6ERV8q6g9F7kZOB1WFhYWYMWMG/vzzT7Nibr1eDwaDAa3Wfj+V5jJ9+nScPn0aubm5dq8bMGAABgwYYPw6MTER3bp1w6effoolS5ZYvWfu3LmYPXu28evKykqEhYW5ZuGEkLuSrV5V1mqf7Dl6qQLFNxT4/rkBuKnUwJ/jDS6LCZVWi8oaNS7L6/69Lau8Cb1OC4aX9ZN8PDYTG1Lj0TaAgy9S+2H5rvNmGS+xSIBJ4ii8vyvfoj0D9Ycidyung6qnn34aer0emZmZaNu2rVsWdc+YMQM7d+7E/v37bWabbGGxWIiLi4NUav0oMlB3OpDD4dzpMgkhxMhWrypJscxObZTArFv6qct13dI35FzEqz+cNF4nFgkwY4gItWodzpbI8VxXDWY+Nhw1nYeAn/Co1c/153jjyc/+wobUeGRZ6Xdl+DouPAhxYYHG16k/FLmbOR1UnThxAkePHrXbpqCl6PV6zJw5Ez/++CP+/PNPREVZ78Nij1arxalTp/Dggw82wQoJIeQ2uUKFsioVKmvVaONrvVdVZm4hMtP6gclgWLRbmDZYhCc/+8u49fbuuJ7YkHPR5uDlEV0FKPq/jcj69jPo9Tr4KXdD3/dhMLzNAzqxSACVtm7bkcFgOGz0GeDDwg/PDwCPzYQvx9thHytCPJXTQVW/fv1QXFzslkHV9OnTsXnzZmzfvh3+/v64evUqAIDP58PHp+4UysSJE9GhQwcsXboUALB48WLce++9EIlEkMlkeP/993Hp0iVMmTKlxb4PQojnK5HVmBWm89hMZKb1gx4wq63qGxGEsHrNLX053jjy7w1M3nTYrJbJXlH7H39JcPjD9bhwti6Ddd+YxzD9zXfwv5PlVrf15Iq67chatf2SDqVGB7VGhwlfHDK+NjBaiGXjY9GeTv+Ru4zTQdUXX3yB559/HpcvX0bPnj3BYpn/F05sbKzLFuesTz75BAAwePBgs9c3btyItLQ0AEBRURG8vG7XKFRUVGDq1Km4evUqgoKC0LdvXxw4cADdu3dvrmUTQu4y1k76KVRapGcdxrxR3fHy/fegRFZjPFG3cMcZLHq4JzqH+AEACkqrMPfH0xbPtVXUrlMrcS37DVxRyODPDwJn8HN4aUkGXvn+BNKTojBtsAhavR5qjQ6SYhk2H7qE7u35ABzXdAX6sHDgonkgt/9CGeZsOYnVKXG0DUjuKk4HVdevX0dBQYFZF3IGg+EWhep6vd7hNX/++afZ1x9++CE+/PDDJloRIYRYsnXST6HSYu7WU9iQGo9p3xwze0+puR2kOFvU7sXiIGjwJHSqPI5X3vkIr/1SDI63FxQqLY4XVaB3WKCxiWiySIhUcaSxH5W9mq4kkQDtA7nIzC20eG//hTKUVamMQZXpVmeADwtCX+qATjyP00FVeno64uLikJ2d7baF6oQQ4s5sBUUG1jJOpkGKraL2U5flSBYJkSO1DNhGjn0co3u9iCJZLZJFQrTxrWuTwGMzEeTLwsrHe4HNrAvKvjl0ybitaLPflUiIuQ92RZXSduf2m7e+z/pbnQBtERLP5HRQdenSJfz0009WhyoTQghxzFZQZGAr42QIUmwNYGYwgOcGRkIHvUWd1PQh0fDnMvH72WuYP7o7FCoNgv3ZqKzRgKFnQKZQY8nOswDquqUrNTrkScuhUGmRkS3BvIe6Y84DXfFfRQ3YzLptyUfXH8S9ndpg9a2gy9CA9FhRBTJzC+HPZdlsamq6RQiAsljEIzgdVA0dOhQnTpygoIoQQhrJVlAE1AVAkmKZ2WuGbuZcFhOSogoE+LCw9JEYLPzpDH4/V2q8Tld0AsMen4WX3v8CG1LjodTojHVZkzcdxvfPDUDPjoH4r6IGX+RetOi0viE1HpM3HTE2DZ02WASmFwMKlRZCPzae+PQvs8J4HpuJCQkRyMorNCuQF4sEyEzrB6Ef2+ZWJwAcuVSBCoUab/98Fl1DAxAXFogr8lqU8lgIb8NDhyBeY3/EhLQIp4Oq0aNHY9asWTh16hRiYmIsCtXHjBnjssURQogn4vPYWDY+FnO2nDQLrAZGCzFtiAjpWYeNr/HYTKxKicPGvEKLBpvvjovB3Ae74doNOda+txizl38GAFizchm+G/2qxedeuqGARqe3CKiAuhOHHCYD22eIUaPSQqHUguXthZwL13H6Pxm6tudbdE1PT4rCRhs9rLwYDKxJibO71ZmeFIW3d57BkwkRFt9fkkiAZY/EomMbCqxI68HQN6S624TpyTmLh7lRR/XmVFlZCT6fD7lc7jZjeggh7s9QvH2zVg1/LgtCPzYUKi1eNwm2ZgwVQVJUYbVQfGC0EOn36PD8lEnIz88HAPj3GYXAwWnwYnEtrt+QGg8AVsfTGIK3TfWyTkkiAZaM7YnH1h9EWZXK4nnWnmWw68Vk6PR6PLjK+mSLDanxkBTLbH5/ydFCrKEThKQJufr3t9OZqvqz/gghhDQOn2dZO8TnwawnFZfFNMvgmPp9/wFkP/8aNBoNQkND8f6qT/CrrC1ybQQokmIZuoda/8VhK+uUKy3H/G1n8PJ9XSzaOBgahNpysawaZ69U2jw9CABxYYE2v7+ceicICXF3zg2VIoQQ0uT4PDY6h/ihd3gQauw032S3EyF+QBIef/xxHPj7GP6obo80cRTEIoHZdUkiAd4ZF4P8K5U2i+DjwgJtBj450jL07MjH0K7BmDFUhA2p8fjk6T6IFPhixlCRzQ7qHG8vZOYWYpKVNQ2MFqJjkI/DgdE3bWwfyhUqFJRWQVJUgYLrVZArVFavI6Q5OZ2pAoB9+/ZhxYoVOHfuHACge/fuePXVV5GcnOzSxRFCyN3O3klBBsMLn3/5HXpEBOPi9WrsPleKAwXlSE+KQro4yqxQXafVYcHoHsiRllnNHDkKbq7Ia5ExrAtW/nbeovZpVUqcsc+VQfKtgnvD6UHDmvg+LATx2MaBy+VV9oMhfyvfP7VoIO7K6UzV119/jeHDh4PH4yEjIwMZGRnw8fHBsGHDsHnz5qZYIyGE3LUMJwWtGRgtRFhIIBgMhrEgXKHSYs1eKSZvOoJp3xzD5E1HkJlbCC2AhTvOoPRmLRaM6oHkepmjQB/7bR7aBXDx/m/nLcbg5ErLkZVXiPSk27NWxSIB5o3uYWwKaromby8GOof4Gbc+IwQ8JNVbi+n3Zwi+DBy1aKCMFWlJTmeq3nnnHbz33nuYNWuW8bWMjAx88MEHWLJkCSZMmODSBRJCyN2Mz2MjmXUR+3Z/AQx7CQxG3X8LJ0cLsWBMD5RX1wURfDtBUXpSFBZuP40caTl2nyvFp/suIj0pCmniukAovA0PAT4sm20eBkYLodPD5vZgrrQcrz/QFd1DA4yZscKyaovTgoBl5qlDEA/LHonF3B9PWWSelo+PtainsteioX4Xd0Kam9NB1cWLFzF69GiL18eMGYM33njDJYsihBACyOVyZGRk4MsvvwQALJ/4GIY9/DiUah0OXCzH6NW5UKi0xvYK93ULMetbZZDYSWC2ZWfIHBnsmT0IbQO4Nts8LB8fixJZjd21Ft+oMRutYzhpaMpa5gkAOrbhYY1Jcb7hJKS14MhRN3pbNViENAeng6qwsDDs2bPHovnn7t27ERYW5rKFEULI3ezPP/9EamqqcQj83LlzMWni05i95YzVra83fjyFpY/EQKnRWQRFjoYiGwIRXzYTSx7uiWqVBgqVFnwfFkL8OeDz2FazTqZMPyNZJETpTaWxaWlcWCAAIOxWM09bcwAbkmFy1I3eWg0WIc3F6aDq5ZdfRkZGBo4fP47ExEQAQF5eHrKysvDxxx+7fIGEEHK3Wb9+PaZNmwa9Xo9OnTrhq6++QmJiIgpKq+xufdWqdWbtGAwZn/r9peoL8GHZLf7m8+pqu5KjhVY/37QLvFgkwKSkSEQL/fBLRjLmbTtllhVLjhZi+q0Gp4ZAzZkic3vd6G1lwghpLk4HVS+88ALatWuHlStX4n//+x8AoFu3bvjuu+/w8MMPu3yBhBDi6epnbu5NGgRfX1+kpKTggw8+gJ+fH4CGbX0ZisDrMw1E6meQvBjA6z+ctBjEfORSBfb9cx3xEUGoVqqxYHR3LPrpjFmxerJIgAVjeuBqZS3iwuqaec7YLMH26WIs3nnWorA950IZdHo90pOijMGW6RxAR9kqe93ordVgEdKcGtVSYdy4cRg3bpyr10IIIXcdWxmiPw9J0Le7eZlFY7e+TAORI5cqLMbebEiNtwioTMfjzN16yvjaWw91w2v1BiuPWZOH1SlxZt3VtTq9zaxanrQc6eIos9ecKTJvH+hjNSNHARVpaU4HVYcPH4ZOp0NCQoLZ64cOHQKTyUR8vGVxIiGEEEv22gMAwOrIcLNA4U62vgyBiEyhxlvbTpmd5LPWo8pah3WFSos3fjwNsUiAuPAgs20902cMjBZCodLY+9atfqYzReYNrcEipDk53adq+vTpKC4utnj98uXLmD59uksWRQghnk6v1+OTz77Arz9+Z/V9Q+bGlCHjVL9vVUO3vvg8NjQ6vcWWnLVCdnsd1vOk5catw/rPMK7Fx/5arH0mFZmT1s7pTNXZs2fRp08fi9fj4uJw9uxZlyyKEEI8WWlpKZ599lls374dDBYX3LAYePNDLK6zlrm5060va3VZkmKZRZd1Rx3WTd9PjhYivA0Pe2YPMq5FrlDZzKqZFrYbUJE58QROZ6o4HA6uXbtm8fqVK1fg7d2oEi1CCLlr7NixAzExMdi+fTtYLBb44ifB9LfeUdxejZRhNqCtwnRbrNVlWZvP56gNg2lm6r3xsYhu62+2FltZteRoIWYOjTZ2Wzc8g4rMiSdg6PV6vTM3pKSk4MqVK9i+fTv4fD4AQCaTYezYsQgJCTGeCLybVFZWgs/nQy6XIyDA+gR4QgiZOXMm1qxZAwDo2bMnPvl8Iz47o7VZI9WQ03DOkitUmJktsfhMHpuJeaO6o294EIorFAj25+D9XZZjaYC6wGj+qO7wYjAcZskMJxtNs2oAqMicuAVX//52OrW0YsUKDBw4EBEREYiLiwMAHD9+HG3btsVXX311xwsihBBPFRUVBQaDgZdffhlLliwBl8tFp641VtsDvDc+FgBQUFpl0STzTthqSRAfEYTBXYJRo9Zi8qYjxtN/OpiPp0m+tbbQBg4utlVQTkEU8UROZ6oAoLq6Gt988w1OnDgBHx8fxMbGIiUlBSzW3VlkSJkqQkhD6HQ6HDt2zOKUtLVsTrVKa7MZZ0OaZDpi7TP5PDYkRRUYt+4AAPN+VkqNDhxvL4S34SG6rf8dfz4h7sDVv7+dDqr279+PxMREi/opjUaDAwcOYODAgXe8qNaGgipCiCvJFSrMyJZY7fPUVNuCBgWlVRj2wT6b7++ZPQidQ/xsvm9rBA0h7qjFt/+GDBmCK1euICTE/KSKXC7HkCFDoNXanw9FCCGeTKfTYe3atXjkkUfQoUOHRj2jrEpldxxNQ5tkNsad9MKyN+rGFdk1Qtyd06f/9Ho9GAyGxevl5eXw9fV1yaIIIaQ1Ki4uxn333YeMjAxMmjQJOp39tgS2NGQcjSm5QoWC0ipIiipQcL0KcoX9WX/2NLYXlr1GpnO2nLyjNRHSWjQ4U/XII48AABgMBtLS0sDhcIzvabVanDx50jhgmRBC7iZ6vR6bN2/G9OnTIZfLwePxMH78eKv/AdoQzoyjaYrsUGN6YdnKrvHYTMSGBeKKvBYXy6ppS5B4tAYHVYb2CXq9Hv7+/vDxuf0PK5vNxr333oupU6e6foWEEOLGbt68iSlTphjbySQkJOCrr75CdHR0o5/Z0C04R9mhO6m9qn9qz5ANM62VAmCsn+J4e2HGUBEycwuhUNWVgZjODzQdaUNbgsRTNTio2rhxIwAgMjISr7zyCm31EUIIAC6Xi3///RdMJhMLFizA3Llz77gRsq22B/W34Jqr9spaNiw5WojpQ0RIzzpsDKLEIgFWpcQhI1sChUprdX6gYW13GvQR4o6c/id/wYIFTbEOQghplVgsFr7++mvIZDL069fPZc9tyBacs7VXzjCc4tPq9Viy44xFE9CcC2XQ6fVIT4oyZqEMwZPhtbiwQLMMlammLrgnpCU4HVQZmtfZcvHixTtaECGEtDZ3stVnj63GmQbO1F45wzQztSE13mpXdaAuiEoXR9l8zdH8wDsJ+ghxR04HVS+99JLZ12q1GhKJBLt27cKrr77qqnURQohbUavV2LJlC5544olGF6C72p20P7Clfp2WM4OVDfg+LGyblggui2n33sYGfYS4K6eDqhdffNHq62vXrsWRI0fueEGEEOJu8vPz8cwzz+Dw4cNQqVSYOHFiSy8JQMNrr5xRv06roYOVTQXdGvgsV6hcHvQR4s6c7lNlywMPPIAtW7a46nF3ZO3atYiMjASXy0VCQgL+/vtvu9d///336Nq1K7hcLmJiYvDLL78000oJIe5Mr9dj3bp1iIuLw+HDhxEYGAgej9fSyzJjqL3aM3sQtk1LxJ7Zg7A6Ja7Bs/nqq1+nJSmWQSwSWL1WLBJAUiwze800WGpszytCWqs7O6Ji4ocffkCbNm1c9bhG++677zB79mysX78eCQkJ+OijjzBixAjk5+dbdIEHgAMHDiAlJQVLly7FqFGjsHnzZowdOxbHjh1Dz549W+A7IMSztNaxJSUlJZg8eTJ27doFABg2bBiysrLQsWPHFl6ZJUe1V86oX6eVmVuIVSlxACwHKxtO/xlYC5Ya0/OKkNbK6dl/cXFxZvUEer0eV69exfXr17Fu3To8++yzLl+kMxISEtCvXz+sWbMGQN3IiLCwMMycORNz5syxuP6JJ55AdXU1du7caXzt3nvvRe/evbF+/Xqrn6FUKqFUKo1fV1ZWIiwsjGb/EVJPax5bcuLECfTr1w9MJhPLly/HjBkz4OXlsuS+25IrVJiZLTHbsjMMVk7sJACXxQTfh2XMRlGwRFqzFp/9N3bsWLOvvby8EBwcjMGDB6Nr1653vKA7oVKpcPToUcydO9f4mpeXF4YPH46DBw9avefgwYOYPXu22WsjRozAtm3bbH7O0qVLsWjRIpesmRBP1ZSNKZtDr169sGHDBvTt2xfdu3dv6eU0G2t1WgqVFieLZXiqf7jFtqI7/29ISHPzqD5VZWVl0Gq1aNu2rdnrbdu2xfnz563ec/XqVavXX7161ebnzJ071ywQM2SqCCG3teRQYFd55plnWnoJLuHsFixt2RHSOI2qqdJqtdi2bRvOnTsHAOjRowfGjBkDJtP+8VlPweFwzGYfEuKOWrqWqSkbU7pSbW0tzp8/j969e7f0UppEY7dgXVmnRcjdwumgSiqV4sEHH8Tly5dxzz33AKjbDgsLC8PPP/+Mzp07u3yRDSUUCsFkMnHt2jWz169du4Z27dpZvaddu3ZOXU9Ia+AOtUxN1ZjSlSQSCZ555hlcvXoVp0+f9rh/7lv7FiwhrY3TVZcZGRno3LkziouLcezYMRw7dgxFRUWIiopCRkZGU6yxwdhsNvr27Ys9e/YYX9PpdNizZw8GDBhg9Z4BAwaYXQ8Av//+u83rCXF3jn6RyhWqZlmHoTGlNS3do0ir1WLp0qVISEjAmTNn4O3tjcLCwhZbT1NpyBYsIcR1nM5U7du3D3/99ZdZ+wSBQIBly5ZBLBa7dHGNMXv2bKSmpiI+Ph79+/fHRx99hOrqakyaNAkAMHHiRHTo0AFLly4FUNfMdNCgQVi5ciUeeughfPvttzhy5Ag+++yzlvw2CGk0d6llaorGlK5w8eJFTJw4EXl5eQCAcePG4dNPP0VwcHCLrKcptZYtWEI8hdNBFYfDwc2bNy1er6qqApvd8mnkJ554AtevX8f8+fNx9epV9O7dG7t27TIWoxcVFZkdi05MTMTmzZvx1ltv4Y033kB0dDS2bdtGPapIq+XoF2mFQgW5onkCK3cseN60aRPy8vLg7++PVatWITU11W3Gzrhaa9iCJcSTON2nauLEiTh27Bg2bNiA/v37AwAOHTqEqVOnom/fvsjKymqKdbo1V/e5IOROFJRWYdgH+2y+vyE1HpsO/NsqekU1BZVKhRdffBGvvfYaoqKiHN/QilnrOWUwMFpINVXkrufq399O11StWrUKnTt3xoABA8DlcsHlciEWiyESifDxxx/f8YIIIXfGXi2TYaxIc9dXuRM2m41PPvnE4wMqgMbEENLcnM5UGUilUmNLhW7dukEkErl0Ya0JZaqIuymR1VjUMolFAkwSRyEjWwKFSgsA2DN7EDqH+LXUMpvUzZs3odVqERgY2NJLaXGG9hrusgVLiLto8Y7qBiKR6K4OpAhxZ4ZapsuyGvxbrgDH2wuSYplZQAV4bqFybm4uJk6ciISEBGRnZ7f0cloc9ZwipHm4bKAyIcR5jWnQ2dB7+Dw2yqpUmPbNMZvP8rRCZZVKhQULFmD58uXQ6/XQ6XS4fv26R57sI4S4HwqqCGkhjWnQ6ew9hvoqW4XKLdkrytVOnz6NZ555BsePHwcApKWl4eOPP6YteUJIs/H8keuEuKHGNOhszD2NKVSWK1QoKK2CpKgCBderWkUxu16vNwZUQqEQW7duxcaNGymgIoQ0K8pUEdICGtOgs7FNPZ3pFeVsJqyl5wsaMBgMfPHFF1iyZAnWr1/vceNmCCGtQ6OCqpycHHz66acoKCjADz/8gA4dOuCrr75CVFQUkpKSXL1GQjxOYzpd30l37IYUKjs7J64p5gveSZDWt29fbNu2rVGfSwghruD09t+WLVswYsQI+Pj4QCKRQKlUAgDkcjneffddly+QEE/UmE7XTd0d25k5ca6YL2i6zXjxehX+u6HAjGwJhn2wD+PWHcCwlfswM1uCElmN2X3V1dWN+O4IIaTpOR1Uvf3221i/fj0+//xzsFi3/yUuFotx7JjtU0aEkNucHTYsV6ig0+uxITUemWn9MGOoCDw20+499e83BDAXrt3EpbJqnCg2r5mqcBAImWbC7nRQb4msxiyA2iq5jNe3Og7Sdu3ahejoaOzYscPu8wkhpCU4vf2Xn5+PgQMHWrzO5/Mhk8lcsSZCPJ4zw4atbbOJRQKsSolDRrYE8RFBdrtj27p/SlIn/JF/Hf0j24Dt7QUfFhMzhoqQmVto1svKwDQTdidbkdayXHFhgVizV2r1+v0XylB0rQKfvL8In3zyCQDggw8+wKhRozx2Zh8hpHVyOqhq164dpFIpIiMjzV7Pzc1Fp06dXLUuQjxe+0AfvP9YL1RUq1BZq0GAjzeCeGy0DeAar7G1zZYnLYcXg4FfM5IRyGMZA6r6NUl+HG/M337a4n5JkQw8thfa+nNQo9ZCVqMGl8VEBz4Xayf0wfTNx8wCq+R6mbA72Yq0luVSanQ2r1eW5GPUsBdRVFgAAHjxxRexdOlSCqgIIW7H6aBq6tSpePHFF5GZmQkGg4GSkhIcPHgQr7zyCubNm9cUayTEYxiCniqlGnwfNuZtO40cqe1Cb3vbbDkXyqDR6Y0BlbWMVHK0EKmJkThQUG4WJD03qBN0emDnqSvIk5YbXxeLBJgxRITnBnXCh79fML62aEwPs0zYnfS/spbl4nhbr0TQVJbh6jevAzoNOnTogKysLAwfPtzmswkhpCU5HVTNmTMHOp0Ow4YNg0KhwMCBA8HhcPDKK69g5syZTbFGQjyCadAzY6gIkqIKs4AGsDxp19BtNlsZrZwLZdDp9UhPijLbXhtyTwiW7zpv8fmGrxeN6YHoEH/jeJv6nNm+rM9alktSLINYJLBYj3eAEN3vewLdAzT47NNPEBQUZPO5hBDS0pwOqhgMBt588028+uqrkEqlqKqqQvfu3eHn55lDWQlxhfpBj6MaIkPPKUMAwmMzkZ4UhbiwQCg1OnBZTBwrqkCAT9379jJaedJypIujrL5u63qFSmscbzMwWoipSZb3O9P/ypS1LFdmbiFWpcSBASDXZF0Do4V499X16Cigf78QQtxfo5t/stlsdO/e3ZVrIcRj1Q967NUQAbczUEI/Nu7rFoIn+odjY16hWSCWJBLgyfgwAI4Lx+t/nkJpWYhuqlatQ2ZaP1yrrMWQLsE2A6XGDOq1luVSqLT49q9/sXx8LGrVOqeCNEIIcRdOB1VDhgyxWyC6d+/eO1oQIZ6oftBjq4bIwFDozeexsXBMD7y25aRFZilXWo65P57C22N7wo9j/x/lQB/zLTcui2njyjpeDCA96zCSo4UY0EmAE8UV8OO6rmN6/SzX/t924PMP34DvuH3oEBJ4x88nhJCW4HRQ1bt3b7Ov1Wo1jh8/jtOnTyM1NdVV6yLEo9SvI7JVQwRYFnrXqnU2t+pyLpRBWlqF0ptKJEcLrW4BJkcLERbkg99nDUS1UgMfNhPyGrXt60VC5N4qns+5UIY3t51CXHgQ1uyV3nHHdFN8Hht6lQILZs/EN998AwD48MMPMfv1N91i9A0hhDjL6aDqww8/tPr6woULUVVVdccLIsQT1a8jMtQQAea1TdYKveU19htpKjU6LNl5Fplp/QDAoh9VamIkRnycg/iIICx9JAZztp7C0UsV2DZdjCU7zpqdPkwWCTApKRIzNkuMr5nWZJkW0gO4o+Bn7969SE1NxX///QcvLy+8+eabSJ/xMmZkS1w6+oYQQpoLQ6/X613xIKlUiv79++PGjRuueFyrUllZCT6fD7lcjoCAgJZeDnFTJbIaszoiHpuJeaO6o094IGpUWqs1RCWyGvxbVo0JXxyy+dwNqfGYvOkIeGwmfnh+AJheDFyW1UCvr8uImTbz3DwlwfisT5/pi1A+F9dvKsFje0Oh0lhcb7DuqT7GwnUem4lfMpIxr17/K2eCny1btuDRRx8FAIhEInz11VfoFtvHIqAyfXb92YP1uctwZ0JI6+Hq39+NLlSv7+DBg+ByuY4vJOQu5expOcOJwV5hgTa3CsUigbHlgUKlxb/lCmT/XYReNk4XymrUxpOE7QO5qKrVgsFgwJfDxORNh612UgfMa8DSk6Iwb9sp5DhoB2HPyJEjER0djWHDhmHFihXw9fVFQWmVw9E3znSNpwwXIaS5OR1UPfLII2Zf6/V6XLlyBUeOHKHmn4Q44MxpOcOJwaOXKqxuFYpFAkwSRyEj+/ZWHcfbCzkXypCWGGn1mTwWE6tS4ixOEiaLBMaMV/3AyjRwAxreDsIeX19fHD16FP7+/sbXqpRqzBgqsmgbYcic2Rp942i4c0OCPEIIcQWngyo+n2/2tZeXF+655x4sXrwY999/v8sWRsjdznBiUKHSIiNbgvSkKMwa3gXymrrXJcUyZGRLjEGQafBjq2WDVq/HxrxCi6xXXdaJgbce6oY3fjxtfD05Woj5o7rjckUNMtP6oVatdXjS0N7cP1OmARUA8H3YkBRVmAVspjMObY2+achwZwqqCCHNwamgSqvVYtKkSYiJiaHOxoQ0MdMTgwqVFmv2So0F7pvy/jUrMK+ftQr0YVltGBoawIGkSGb183KkZXhxeDQ2pMYDAEL5Pvgj/xr+q6jBF7kXjYGY4X1buCwmJEUV8Od645fvv8Z/lwrx3nvv2b1HrlBh3rbTVju8e6Guw7ut0Td3MtyZEEJcyamgislk4v7778e5c+coqCKkiVnrPG7IWmWm9cMLgztDVqM2jpIxZK0GRgsRKeAhM60fVu+9YL7NFy00Zn6s1U+ptXpwWUz4c70xMfMQnuwfbhZQAfbbQSRHC1Gr0aKg6DLmvzIT5w79CQAYO3YsEhMTbX6vZVUqsyDRVI60HG+N6m4z23Qnw50JIcSV7HcgtKJnz564ePFiU6yFEGLC0Hl8YLTQ7PX4iCBEtuEhUuiLb/8uwuRNR7Bmr9QYUC1+uCcUai3W7r1guc13oQwb8wqRbmXsDAAE8lg4eLEcT372F8qqVIgLC7R4xrd/F2H+qB5IEgnMXk8SCTBvVHc8+vrHeOJ+Mc4d+hMsFhsJT76IbrF97H6vjrJN/1XUQK6w3lrCEHxa42i4MyGEuJLTNVVvv/02XnnlFSxZsgR9+/aFr6+v2fvUUoAQ17F3YlCuUGHeqO6Q1ajhx2aCy2LiyKUKPLgqB6tT4ixO5xlYmwVoaO/AYABJnYUY0b0tAAZ0eh14bCYUKq1xO/H+7m3x/q5z6B0ehEniKCg1OmO27LEpGSjYmQUAYAVH4v5pS7B0ymhU1moRZGd8n6NsEwCbtVF3MtyZEEJcqcFB1eLFi/Hyyy/jwQcfBACMGTPGbFyNXq8Hg8GAVmt/phghxDnWTgxaayGQJBIg7Vaw5Gi2oCkem4nMtH5Yu/cC5m49ZXxdLBJg5pBorJ3QB6/+cALLxsdiY14h4sICsfv8dew+f93iWTW+IngxmfCLH4fApKdwupaF0pu1eO+3url+ttobCP3YNju8i29lxJQaLSS3hkjX70HV2OHOhBDiSg0OqhYtWoTnn38ef/zxR1OuhxDigK0WArnScjAAfPfsACg19v/jRuDLwYbUeCg1OnQM8sH7u85bZLbqisQZeHXkPcia1B8n/5NBUiSDMsF2wOYT2Ruf78jD4n3mw6NzHLQ34PPYWPxwD7xVr1hdLBJgwegeWPbrOew1CeKs9aBqzHBnQghxpQYHVYbG64MGDWqyxRBCHLPXQiBHWo60m7Uoray1mflJEgmw+/w1rNkrBY/NxJfp/W1uFeZIy5B2MxKTNx0xtjfwtjNQHQAiIiMBk6DK0DjUUXsDJhh4KCYU6SZbigAsAirDs6gHFSHE3ThVU8Vw8C9TQkjTMxR1W2uZcKyoAhqdHkt+PlfX+kCvNwuYkqOFSE2MREa2BDx2XSNQQ98rWwxbiYYMUtqACLS5+BtKWB3ADetpdm39RqH1v7bX3kBeq0LPDnzcrNVAqdGBwWAggOttEVAZUA8qQoi7cSqo6tKli8PAqqVm//37779YsmQJ9u7di6tXr6J9+/Z4+umn8eabb4LNtv0v3cGDB2Pfvn1mrz333HNYv359Uy+ZkEYJ4LKMAZGhM7ohwBrQSQA20wtrJ/TBocIbeP2BrkirVEKp0SHQh4VQPhcPrc6FQqXFjKGiupOAYusnAQ1MR9TsO3oOZz57GScOH4CvIBTsiavhxa4bT2Wo6TL0yrLW8d1eewN/Lhtvbjtltv3nqCcW9aAihLgTp4KqRYsWWXRUdxfnz5+HTqfDp59+CpFIhNOnT2Pq1Kmorq7GihUr7N47depULF682Pg1j8dr6uWSu0RTDPkV+rExb1R3Y2f0+gGWgVgkwAM922GmSU+qX19MNv7dMG4mLjwIySKB1S1AQ6ZJr9dDff5PVOz+FJcVVeD68DDr1dfx6IShqKhWo20AFyXyGrTxZWH1rZE69Tu+J9tpb2Cr+acj1IOKEOJOnAqqnnzySYSEhDTVWu7IyJEjMXLkSOPXnTp1Qn5+Pj755BOHQRWPx0O7du2aeonkLtNUQ375PDb6hAcaT+qlJ0VZHT2TJy3H2zvPIT0pyhhsKVQaJIkEyJWWG7f1MnMLsWZCHACG1S7tM78+jIqfluHm+TwAAKd9V7QZ9TK+qghF4W/5mDeqB4orFJiy6Qhm3ReN/pFtsOYPqUXB+fQhIpvfk63mn/YajVIPKkKIu2lw88/WWE8ll8vRpk0bh9d98803EAqF6NmzJ+bOnQuFQmH3eqVSicrKSrM/hJhyNOTXViPLhqqs1Rj/bq1Bp0GOtAxxYYHGrwN92Hh3XAySRALjtp5CpcWMzRL0jQzCD88PwP+eG4DNUxIQFx6EjGwJarQMdIsMBYPJRODAiWj71HKwgkIB1J04XLzzDEL8OQCAnu35mLzpCOLCg7AhNR7rnuqDDanxiAsPQnrWYZRVWf++bTX/zMwtxCRxFJLrNfekHlSEEHfk9Om/1kIqlWL16tUOs1QTJkxAREQE2rdvj5MnT+L1119Hfn4+tm7davOepUuXYtGiRa5eMvEgTTHk13Qr0XSosaOeVIb3DZkdPo+NlY/3hkKpMZ4QVKi0+Gj3BXy0+wJ4bCY2pvXDiSKZcevurSXL8Fy7ZLDbdrZ4fp60HEwvBtY91Qch/hykJ0UhM7fQ6hgcWzVQ1pp/GurEvBgMvDgsGm8+1A1MBgNMLwYELthGJYQQV2twUKXTNbyZoCvNmTMHy5cvt3vNuXPn0LVrV+PXly9fxsiRI/HYY49h6tSpdu999tlnjX+PiYlBaGgohg0bhoKCAnTubPkLBADmzp2L2bNnG7+urKxEWFhYQ74dcpdwdsivo9qr+luJM4aKjNt4poXk1nC8vZAsqhtfY9A2oK64fLmVTuT3dmqDIF82HohphzRxJJQaHdrw2FYDKoNL5QpM++YYACBZJMC2aWIUlleDxfTCsaIKY5Blqwaq/pxDW3Vihu1TCqgIIe7I6TE1ze3ll19GWlqa3Ws6depk/HtJSQmGDBmCxMREfPbZZ05/XkJCAoC6TJetoIrD4YDD4Tj9bHL3aMiQX7lChfJqFfQAFm4/bVYoblp7ZW0rMTO3EKtMCsJtDjgWCRHsz0Gv8EA8uCoH8RFBZjVdhk7kV+S1+OeKDLu+34RRvZ7Boh1nrJ7Cs9XGgcdiGq/NkZZj0c4ziAsPwpq9UmN/q+/+LrJaA2UIKDOGReOFwZ2RV1AOFpNhtU6M+lMRQtyZ2wdVwcHBCA4ObtC1ly9fxpAhQ9C3b19s3LgRXl5Oz4vG8ePHAQChoaFO30uIQf3Mi6mB0UJwWV6YkS1Br7BAnC2Ro1d4ENJuNb00BCoLtp/Gisd6Wd1KVKi0yMiWID0pCg/FhmJMbHss+fmsxdiaVHEknvzsL+NWnLWgRF6jxsFjJ/H8lElQXZXi2uViHG3/sNnnSYplGNo1GBMSIiyyR0kiAYZ3a2ucEQiYzxfMu9Xp3VoNlLVi/uRoIRaP6YFP91kf3E79qQgh7srtg6qGunz5MgYPHoyIiAisWLEC16/fbhhoONl3+fJlDBs2DF9++SX69++PgoICbN68GQ8++CAEAgFOnjyJWbNmYeDAgYiNjW2pb4V4AHtDft8dF4MFP51BzoUyTBZHoXdYoNV2CJPEUSivVtncSlSotMjMLcTo2FDooMebD3aDWqdHhUKFNjw2dp25atbSwMA0KLkiU+CRmfNxaus6qJS18OL6oUtsXxytFwtm5hbi22fvxfJd5y2yR7nScuhx3uyUIWBe65UrLUet2ryEwFYxf86FMsz/6YzF80xRfypCiDvymKDq999/h1QqhVQqRceOHc3eMxTZq9Vq5OfnG0/3sdls7N69Gx999BGqq6sRFhaG8ePH46233mr29RPPY2vIb3m1CrvPlQIA+DwWVvxfvtV2CACwcHQPm1uJhrqjt3eeNds6FIsEmD+qu81icaAuKLl69Soee/IpHN23FwDQrnsCmINfwPCHRiJ70xGz6xUqLcqqVDZPGZpmpgzq13rVD4Tsjtu5UIa0xEir7wHUn4oQ4p48JqhKS0tzWHsVGRlpdooxLCzMops6Ia5kbcjvxbJq49+9mQy7gYpWp0eAL8vqHD9n+lMB5vVQap0elytV+Of8WTC8OQgZPhlfrpyHrAP/2qzRqlXbH9JsmpmqP54GsAyEHBXz20L9qQgh7spjgipCWgvTzFNFtf3AolqpwYr/y0dqYiR0er15Q83OApvbYznSMrwwuLPxfVun6eLSFuKJ5B5YJ6nBmZJKzBreBRqdHqNiQyG5VIElP58zZrv4PvazQ4bMlLXxNEkiAfy4t/91I1eo4GNS3G5NxyAfi7o06k9FCHFnFFQR0sxMi9g1Ovv931RaHXafK8WBgnKkJ0Uh/VYxexseC95M+wcx2N5exs+xldXKR0fkXWfji9SeeP+38/ho9wXje8nRQuycmYSyKiXUWj1q1Vq7pwwFfmzsmCnGqf/kZrVc4lszAauVdQ1LDcXpvcIC7XZLbxfAtbp9SgEVIcRdUVBFSDMzLWI/VlRhO1CJFiKvoO51hUprlmGaMVSExE4Cu58TxGNjdUoc/rsuR87+P5EntT7TMqZjIFb+ZlmAnnOhDPO2n8bDvdqjV1gQlu86h0kmJ/oMxLdOGU74/BCAum3J7567F8U3asDx9jLOANw8JcGsOP3opQpjW4i8eu0kTLNRFEQRQloLhr61tUp3Q5WVleDz+ZDL5QgICGjp5ZBWwqxP1a3TgAYDo4VYMKYHRq/OtVpsnpnWD8eKKiApqrAZkK1JiYP03Ck888wzyM/PR8iE5eB06GZ2HY/NxHfP3YvRq/NsrnNDajy+OXQJTyVEgMlggM9jgcti4ka1CkqNzqy5J1AXZBl6VJnWcQX4sODP9cbIj3LMPt+071UnoS9C+VwKpAghzcLVv78pU0VICzEtYl9j45SgrdN73l4Mswag9TNHbz14D1a8twzL3lkCjUYDYUhb6DSWc/fSk6JQWqm0u06lRoe956/jqYQIpGUdhlgkQP+oNujRnm+xpZh0a6svI1titY5r3VN9zJ5dPwO39YVElFWpcLGs2mpneUIIcWcUVBHiBqydEgRgs4GowJdt1gDUUGvF8fbC73+dwJNjRuD4kb8BADFJ9+ODVWvB4PqDwWCYZZZMhy3bYihAV2p0ZkETALz1UDfMeaAr/quoQUQbHvg+LCz86UzdkOahIougy9FInVq1Fo98csDs+zftAE8IIe6Mtv9cgLb/SFMpkdVYbSD65kPdsHjnWatbf2Hyk8hd/wa8ODy8smg5itv0xYGCG8b3TU/nrU6Jg6RYZnsbUSRAfFQbqLV6PBQTCrlCDT+uN9RaHeQKNY6YBGh7Zg9C5xA/49gZpUaLB1flmj1vxlCRzc9KEgnQ+9a2oamB0UIaS0MIaRK0/UfIXcReA1FbReOTxOnY4q9Em9hB+LfWBwes9LHyYjDwa0YyVFodZmZLbG4jzh/TA8U3FNiQW2hWIzWgkwBcFhN9w4OAJOD8lUpj7yhD1u3YpRuoz9aWZXK0EKmJkWZtGAxoLA0hpLWgoIoQFzJkaSpr1S6rCbK1Nbj0l3OICw8y2/qTFMvw7aFL6Do8BXFhgZhcrzO6Qc6tdg4h/hzERwRZ3Ua8VlmL/24FVHnScpu9rpJEAiwc09Ps+SWyGouxNID5zMJ5D3VHrVoLfy4LWr0eY9fm2e0ATwgh7o6CKkJcxNpw4KaqCeLz2HhlSDjmfJOLNXu5xtfFIgEWjO6BsWvzsOKxXnafcbNWjc4hfsb2DqaBUnK0EIvG9IBKqzNmlGz1usqVlmPBT6fx7tgY8HlsY9sEW32oFCotThbLMDUpyhgsFpRWAajbHjScBDQMls7MLaSxNISQVoGCKkJcwNZw4P0XyjBny0mbNUGNzWzt378fEydOhK+fP3759U8otAzIa9SQFMvw6+kriAsPdFgUbghUbG0x8nlsSIoqjNfHhQXa7OCeJy1HjVqLC9duQqZQ47UR9wBgYGSPtvjg93+w9/ztAefJVrqiC/3YyEzrh8/2Fxg/q1atRWJnAR6KCaWxNISQVoGCKkJcwN5wYFs1QY3JbCmVSsybNw8rVqyAXq9HVFQUWLU3EN0hAjOzJdh/ocy4TXetstZux3LTQMXWFqPpSB3T2X7WlMhrkJ51e7tRLBJgxhARXhreBRPvjYRCrQXH2wvhbXgItfL9fbH/IiYkRFjdXlwytifKq1UQUIsFQogbs/+fsoSQBnE0HLh+TZCjzJZcYdlT6uTJk+jfvz/ef/996PV6TJ48GSdOnEB0dLSxS/vAaKGxbul6lRKLRvdAcrTQ7DmO5ufJFSoUlFZBUlQBbyYDA2/d7yjz1S7AxyxQy5OWY80fUpy+LMeRogpM++YYJm86Ai8Gw+LesioVurYPsLm9OH/bGWyVXMbMbAlKZDV210EIIS2FMlWEuIAfx/4/Sr713m9MZuvdd9/FyZMnERwcjM8//xwPP/yw2fu2tvGsNRa1FVDVz57x2ExkpvWDHoCkWIZkkQA5VjJfYpEAJ4or8EVqP0z4/C9jwXmetBzp4ii0Dair+0qulyEzqKxV291ezJGWIU0ciTV7pXa3UwkhpCVRUEWIC7CZXja32sQiAdj1hh87m9kCgNWrV8PHxwfLli1D27Ztrd5naxuvIQGIteyZQqVFetZhzBvVHX0iAjE+rgPmbz9tFliZ9r3qGx6E9KQos+BIqdHBj+uNoV2D8ezAzlY/O4DLwhV5rd31GbYfqcUCIcRdUVBFiAvIauz1jYqCvEYFwNf4eoCD02zWTrsFBwdj48aNrlmwFbayZwqVFnO3nsKG1HgwGQy8OrIr0m4qzdo4ZGRLoFBpkSMtw4vDo81O8An92FAoNcgY1gUzNx/D15MTLIrzhX5sXKu0/zMx3X6kFguEEHdEQRUhLuDHYSHl80MWvZ4MAceOGUlm1wv92DZH0PRv541AbvOXOzrKnhm+p/8qajDtm2M2r5PXqM36YyWLBHh5RFd8uk+K1RP64K1tp8wyXYbi/EgBD0kiAXJtZPskxTLj19RigRDijqhQnRAXEPqxER9RN2Jl8qYjxqLsNXuliI8IsqgjMi0sN9WpJh9/Lp2IdR+vbM7lA3CcPTMEiXwf5wKaHGk53v/tPB7tG4aVv523qMkyFOf7cryx7JFYi8J6Q7YvM7cQgOXJRUIIcRcUVBHiAraCJHsn7QyF5XtmD8Lm1FgMKNmCP1a9jOulpdiyZQvU6ubd4jJkz6wxZIrqGnF6I0kksHtdfXnScoQEcKwWuQO366Q6tuFhza2fyZYXBmDzlATEhQcZtxcdnVwkhJCWRNt/hLiIvSaatvB5bJw9cRQTn3kGBQV1jS9nz56Nd955ByyW7YxQU43DMXRXN92WNC1ENxSub3k+EW9tP21Wg2Vvfh8AVNVaH0FjYKiTMi22lytUaBvAxfCuIQ36eRJCSEuioIoQF7J1+s6W6upqjBkzBmVlZQgLC0NWVhaGDh1q956mHIdjCAxLbypRdEMBAGaF6ADQPTQAgTyWRasGby8GHliVY3N+XwDX/r9urNVJOfvzJISQlkTbf4S0IF9fX6xevRpPP/00Tp486TCgakzTUGfxeWxEt/VHt9AAbDrwL9bslRoDJdPtNz6Pjc4hfugdHoTOIX4I5LEQHxFk9ZlJIgH4PJbN7UWqkyKEeAKGXq/Xt/QiWrvKykrw+XzI5XIEBAS09HKIBysorcKwD/bZfH/P7EHoHOLnss8zbDM2dDuzRFZjsX2YHC3E0nEx6NiGZ/V9Q6BmbXQNIYQ0JVf//qbtP0KaSW1tLbhcrtP3mdZPaXT2/xuoIf2bnKnHcnb7zVFdWWPqzgghpLWgoIqQZpCdnY3Zs2dj9+7d6NGjR4Pvq18/tSE13u71jvo3NWU9loGjQIzqpAghnopqqghpQhUVFUhJScGECRNw9epVfPDBBw2+11r9lKRYBrGNdgaO6pKaox6LEELuZhRUEdJEfv/9d8TExODbb78Fk8nEwoULsX79+gbfb21sTGZuISaJoywCq4b0b2rIEGdb5AoVCkqrICmqQMH1KgrACCHECtr+I6QJHDhwAPfffz8AoEuXLvjqq6/Qv39/p55hbWyMQqVFRrYE6UlRePPBblBpdA2uS2rMEGegebYMCSHEE1CmipAmMGDAAIwePRrTp0+HRCJxOqACbI+NUai0WLNXCo4309jOoCE1So0Z4kxbhoQQ0nCUqSKkCTAYDGzduhXe3o3/R8ze0OXG9HVqzPMasmVIReeEEFKHMlWE3CFbrd7uJKACGjdP0NXPa+yWISGE3I0oU0VII+n1enz22Wf48ccf8fPPP4PJZLr8M1zd18nZ5zVmy5AQQu5WHpWpioyMBIPBMPuzbNkyu/fU1tZi+vTpEAgE8PPzw/jx43Ht2rVmWjFpra5cuYJRo0bh+eefx2+//Ybs7Owm+6z642BcMTi5oc8zbBlaQ6NlCCHEnEcFVQCwePFiXLlyxfhn5syZdq+fNWsWduzYge+//x779u1DSUkJHnnkkWZaLWlOrmoLsHXrVsTExOCXX34Bh8PBhx9+iAkTJrh4te7RxsDVW5CEEOLJPG77z9/fH+3atWvQtXK5HBs2bMDmzZuNg2w3btyIbt264a+//sK9997blEslzchVbQFWrlyJV155BQDQu3dvfP311051SG/u9boCjZYhhJCG8bhM1bJlyyAQCBAXF4f3338fGo3G5rVHjx6FWq3G8OHDja917doV4eHhOHjwoM37lEolKisrzf4Q9+XKtgCPPvoogoKCMHfuXBw6dKhJAip3bGPg6i1IQgjxRB6VqcrIyECfPn3Qpk0bHDhwAHPnzsWVK1dsjga5evUq2Gw2AgMDzV5v27Ytrl69avNzli5dikWLFrly6aQJubItQEREBAoKChAUFOTKJZqhNgaEENI6uX2mas6cORbF5/X/nD9/HgAwe/ZsDB48GLGxsXj++eexcuVKrF69Gkql0qVrmjt3LuRyufFPcXGxS59PXMvVbQGaMqACqI0BIYS0Vm6fqXr55ZeRlpZm95pOnTpZfT0hIQEajQb//vsv7rnnHov327VrB5VKBZlMZpatunbtmt26LA6HAw6H06D1t2ZyhQplVSpU1qoR4MOC0Ld11tE42xZAq9Vi5cqV4PP5eO6555pyaVZRGwNCCGmd3D6oCg4ORnBwcKPuPX78OLy8vBASEmL1/b59+4LFYmHPnj0YP348ACA/Px9FRUUYMGBAo9fsCZqjULq5gjZnOokXFhYiNTUVOTk58PHxwQMPPIDw8HCXr8keV3dSr89TgmVCCHE3DL2tdtCtzMGDB3Ho0CEMGTIE/v7+OHjwIGbNmoUHHngAmzZtAgBcvnwZw4YNw5dffmmcxfbCCy/gl19+QVZWFgICAowtGA4cONDgz66srASfz4dcLkdAQIDrv7lmJleoMCNbYrWuZ2C0EKtT4u74l3Bzn24rkdVgzpaTZoGKoS1AaKAP9Ho9srKykJGRgaqqKvj5+eHjjz/GpEmTwGAwXL6eO13vnTzXXU4VEkJIS3P172+3z1Q1FIfDwbfffouFCxdCqVQiKioKs2bNwuzZs43XqNVq5OfnQ6FQGF/78MMP4eXlhfHjx0OpVGLEiBFYt25dS3wLbqOpC6UdnW5zRdBWn6O2AE8//TQ2b94MAEhKSsKmTZtsbis3h6ZoY9ASP3dCCLmbeExQ1adPH/z11192r4mMjLSY08blcrF27VqsXbu2KZfXqjR1oXRLnW7j82wHJYmJifj++++xZMkSvPLKK00ycsZZ9tbbGHSqkBBCmpbHBFWepiXrXpq6UNodT7dNmzYN9913H7p06dLsn91c3PHnTgghnoSCKjfU0nUvTV0o7Y6n2xgMhkcHVIB7/twJIcSTuH2fqruNO3TTbup5by01pFelUmHhwoW4dOlSkzzf3dFwZEIIaVoec/qvJbny9EBBaRWGfbDP5vt7Zg9C5xC/O/qMhjJsQZoWSgNwybZkU51us+XMmTN45plnIJFIMGjQIPzxxx8tcqqvpTX3z50QQtwZnf7zcO5U91K/ULqh25INqQdrriG9Op0Oq1atwpw5c6BUKiEQCDBz5sy7MqACaDgyIYQ0JQqq3Iy71r009Di+aeDFYzORnhSFxE4CsL29EOTLthpg6QGgCWKciooKPProo9i7dy8A4IEHHsCGDRsQGhrq+g9rRVx9qpAQQkgdCqrcTFMXiTdWQ47jAzALqFalxGFjXiHW7JUarzVsNelNrjV9z5XF+AEBAVCpVODxeFi5ciWee+65uzZDRQghpOlRobqbaeoi8caqcFAgf7NWbRZ4pSdFYWNeIfKk5WbX7b9Qhj//uY7Xf2j6Ynwmk4mvvvoKEokEzz//PAVUhBBCmhRlqtyQu9W9yBUqqDQ6u9f4c1lm9WBxYYFmGSpTIf4c5EibpwllZGSkS55DCCGEOEKZKjfF57HROcQPvcOD0DnEr8kDKrlChYLSKkiKKlBwvcosW1RWpcKBi+UQiwRW702+tS1pWg+mtBOE2XsPcL4YX6FQYP369Rbd8gkhhJDmRJkq4vBUX5VSDRaTgQWjemDJzjPIMdnSE4sEWDSmhzHoM9SDcbxtx+v23gOcK8Y/fPgwnnnmGeTn58PLywvPPvtsg+8lhBBCXIkyVXe5hjQb5fuw8XfhDYxdl4de4UHYkBqPDanx2DFTjH6RbYz3mNaDSYplNjNbpTeVd9yEUq1WY9GiRRgwYADy8/PRvn17REVFOfGdE0IIIa5Fmaq7nKNTfaU3lVi846yx4Ny0TkosEmBUbHsIfG8HQYZ6sPJqFcbFdcDCn85YZMCGdAnGoC7BNptQ8nlsu72upFIpnnrqKfz9998AgMcffxyffPIJ2rS5HeARQgghzY2Cqruco2aj8hq1zaLyPGk55o/qblHvZdoHaY2dgntbxfiOtiNVKhVOnDgBPp+PdevWISUlhU72EUIIaXEUVN3lHDUb9WEzkZnWD8eKKpCZWwiFSmv2/s1aTd0WoY1CenuNJq2915Amo927d0d2djbi4+MRFhbm6FskhBBCmgUFVR6qIaNiAPvNRsUiAX49fRVr9kohFgmwZkIcTl2WI6ZDILy9GGjjywaTwcD5q5XoGMRDhyDeHa+7IU1G+Tw2xo0bd8efRQghhLgSBVUeqKEz+oC6bNG742Iw98dTZteLRQJMEkchI1sCoG6rzwsMPBDTDulZhy2uW7D9NN4c1R2VNWr4cRs/aNl0O1KnrIbqehG4HbsZX2vO2YeEEEKIMyio8jANndFnUCKrwcIdZ9ArLBBpiZHw5XijWqmBpFiGjGyJ2XZfjrQMaeJIs+caCtjjwoPw1rbTiAsPwpq90kaPnDFsR9ZeOomynz+EXqVAaPpaeAfUnRZsqdmHhBBCiCMUVHmY8mqVMUBSanTgspjGeqj63cpNA7Dd50oBAOue6oNp3xyz+XxrjTvzpOVIF0dhzV4p0sV1bQ1MgzgADdqKBAA/bx14x77Bpd+zAQDegaHQ1VYCAUKLdgsN3eIkhBBCmgMFVR7oeFGFWeuD5Fv1UDM2S4zbZ3KFClfktRYZLUeNOQ3v89hMpCdFIS4sEEqNDiH+HMwYKoJGd7ur+f4LZbhaWYu3fz7XoK3I48eP4+mnn8a5M2cAAH69RiJo6GR4sX0sZh86s8VJCCGENAcKqjyIXKHCgu2nkVtviHFdB3QGnh3YCQE+LGNAktI/3OIZp0vkeHdcT7QN4FpkuuLCAyEploHHZmJVShw25hVa9K0aE9sePDbTuG34X0VNg7Yi5QoVvvx2K86cOQOBMBgfrl0P8ZD7UVlj2YrB2S1OQgghpDlQUOVBSm8qzUbImMqRluGNh7pBq9NhX3450hIjEezPMbuGx2YipgMfG3MLLUbRZKbGw4vBQFrWYaQnRWFjXqGxnsogT1qOJTvPID0pyuYwZQPTrUhDkLdfF4+Aex+DT7+x+P1mBwxjMdEp2M/i3oaeECSEEEKaE42p8SCyGvsn4ypr1NDrGdh56gombzqC/zt7zWyUTHpSFDbUC6iAumBp7R9SBAdw0DciCHFhgRYBlUGOtBxxYYEA6gYtS4plNtdzs1ZtlnVieDERNCgVTB7fbEyOxffh4AQgnRAkhBDSEiio8iC+bKbd9/253liy84wxIMrMLcQkcZQxsHIULNWqdcgYFo1QPhczhorAs/F5So0OA6OFWPJwT2TmFlq8r62WQauQw5/LalDWqT5HDUtb6wlBuUKFgtIqSIoqUHC9ympASQghxH3R9p8H8WV7QywSWA2MkkVC8NhMHC2SGV9TqLTIyJYgPSkK6eIo+HG8LQrQTWuqCsuqjScDxSIBVqXEWbRdAIBOQl/jqb/4iCCzxqKKC4dQvmsVOtzTGwLf8SgsV9j9nqxlnew1LG3oQGZ3Q4X3hBDS+lGmyoME8liYOTTabEsPqAuAUsWRePvnc1iVEmeWYVKotFizV4qZ2RIE+3Hw3XP3Ii4sEAwGA2evVGLG5mOQFFXU3ce6fV+etBwb8wqRnhRl9lkDo4UI5XONI2iWjY/FwGghdEoFyn9dhetbl0CnkMNPWQ6dsrpRWSfT59b/bNMTgq2Fo8J7ylgRQkjrQJkqD8LnsRHRhoeZQ6MxbbAI8ho1ON5eZo08a9Rai0Jyw2m+BT+dtihQN2SjGKjbKjRl6E9lYC2oaR/ogwkRCvy5+DVUFf0LBoOBGRkv4b1l74LL5YKpUDUq69Q+0MfmQObWhgrvCSHEM1BQ5WFCA31QpdTgvg/3W33fNBAybPXd370t3t913mqBOgBjEJae1Mnieb4cb6x7qg843l4ovam0eF+n0+GN115GcdG/CA8Px5dffolBgwYZ3zdkneZsOWkWWDUk62RvWHNrQoX3hBDiGSio8kBVSo3Da0x7TcWFBdpsxWAahPmwmMhM62essVKotKhWasw6sA+MFpr1ifLy8kJWVhY+/vhjfPDBB+Dz+Raf0disk6d0VPfUwntCCLnbUFDlxhobNDj6Jd2Oz8XWFxLxzs9nkSctx1MJEXavN4ymUag0mLzpiHFbcPOhSxYtE/ZfKEN5dV0NkEyhRrVKA2VAR8x59yOAxYEtzmadPKmw2xML7wkh5G5EQZWbupOgQejHRnK00GqdTpJIgF9PXzXLTjVkNI1YJDAGUHnScjAAvDayK5787C/jNmInfx18/AOh1ulwWVaDd345Z3YSMfnWlt6dBj2e1lH9TrZACSGEuA86/eeG7vQ0GJ/HxtJxMUiqdwowSSRAmjgKmbmFZoORJcUyixODBmKRANcqazHp1n0GudJyYw+pj5/sja3/y8ZjQ/oh+7v/QanWYWm9gAoAci6U4XUXnGZrTG8rd2fYAt33ymD8kpGE758fgHmjutvsBUYIIcT9eExQ9eeff4LBYFj9c/jwYZv3DR482OL6559/vhlXbskVQUPHNjwsHx+LzVMSsO6pPtgxU4ze4UHGU4Cm2an6TUANkkVCLBjdA9erlFb7UfHYTLz9QASmTHwKBzcsgk5ZjT3b/4fr9sblXCizWtDuDE8t7K5WafHW9tN4cFUuHlt/EPd9uB8zsyUokdW09NIIIYQ0gMds/yUmJuLKlStmr82bNw979uxBfHy83XunTp2KxYsXG7/m8XhNssaGclXQ0CGIBz+ON65U1kJWrTZro2DITuVJyy2agAJ1Pa/4PiyMWZNnEUwZ/LH7/7By3iyUl16DF9MbPUelQ9blIbMsmDUyhRqSoopGF5d7YmG3p21pEkLI3chjgio2m4127doZv1ar1di+fTtmzpwJBoNh914ej2d2ryNKpRJK5e1sS2VlpfMLtsOVQQOfx8a/5dXwZpr/DDJzC7HqVtdzQ2C1Zq8UybfGy2h0OjAZDPSNCLKaNesdpMJb0yZCp9XAu01HfLz+C7x3VAMGHNdoVdaqMXnTEQCWdWINKc73xMJu6lVFCCGtn8ds/9X3008/oby8HJMmTXJ47TfffAOhUIiePXti7ty5UCjsj05ZunQp+Hy+8U9YWJirlg3gdtBgTWOCBj8OC3/+c91se8+QnYoLD8LmKQn433MDsGNmEt4d2xNsJgOllUr8U1qF+aO6I1lkvpZkkQAZDyeiTeLjGPvUZISmfYTIrrHG9+3VaBmeZagVMmRiLpVXo+B6FV7+33H8ePwyblSrkH/1Js5eqcTlCvP/PTytozrguVuahBByN/GYTFV9GzZswIgRI9CxY0e7102YMAERERFo3749Tp48iddffx35+fnYunWrzXvmzp2L2bNnG7+urKx0aWDlitNgphkfP443Itrw0Dc8CF6Asd5JodLieFEFeocFIuPLusxRZlo/rN17ATnScvDYTKyd0AcPxrRDmjgSSo2urslnpRJCPw7C70uFOCECvKIKixot0yyYgWFczjeHLpnNDdx/oQzXbyqx7k8pJiREYGNeodlWZZJIgGWPxKJjm9vbsp7UUR3wzC1NQgi52zD0er2+pRdhz5w5c7B8+XK715w7dw5du3Y1fv3ff/8hIiIC//vf/zB+/HinPm/v3r0YNmwYpFIpOnfu3KB7KisrwefzIZfLERAQ4NTn2WMIjJwNGqy1Y0iOFmLGEBF8OUxcq1SC5e0FJoOBgxfLjY08ZwwV4WyJHN3b8xEXFog2vmx88H/52J9/DQymefydJBLgtZFdkZ51GMvGx+L6TSV2niwxBlE8NhOrU+LAZTHNxuUYPkssEiAuPMgYPH377L3IlZZBUlRhfSB0tBBrPLiuSK5QYWa2xOaWJtVUEUKI67n697fbZ6pefvllpKWl2b2mUyfz8SkbN26EQCDAmDFjnP68hIQEAHAqqGoqjRnDYqvgOedCGRgA3hnbE8t35SMtMRITb9U1GcSHB6F3WKAxU/TFxL745X+bUHlkO9o9sxJMn9v/h8uVluOFWg2e7B+O7/4uwjvjYjCgkwDztp9GzoUyY3H7U18csrrO+nMDeWwm4sICzTJU9dfvyXVF1KuKEEJaP7cPqoKDgxEcHNzg6/V6PTZu3IiJEyeCxXJ+y+T48eMAgNDQUKfvdQeOCp5VWj2Wj4/FuSuWxfV8Hgsr/i8fedJyaG6W4/WpKbhx4E8AQNXxXQgdlIL0pCjEhQVCqdEhgOuNh3u1h75XKC7LasD3YWHeQ91QXFEDpUYHX7b9/3sZTgka6q8cnRr09LoiT9vSJISQu43bB1XO2rt3LwoLCzFlyhSL9y5fvoxhw4bhyy+/RP/+/VFQUIDNmzfjwQcfhEAgwMmTJzFr1iwMHDgQsbGxVp7evBozpqYhBc+dQ/yg0mixITUeSo0OXBYTx4oqwGJ6IU9ajupzObjxf2txubYKDG82AgdPQtt7xxhnBZpmk5JFQqSKI431UZunJBhP9m1Itd/KwtCpfZI4Cn/klyI+oo3d6/25LI+Z92eLpwyJJoSQu5HHBVUbNmxAYmKiWY2VgVqtRn5+vvF0H5vNxu7du/HRRx+huroaYWFhGD9+PN56663mXraFxo6p8ePY/5/Un8tCiawG87adQY7UpOZKJMSg6GBUnfkD5TtXAgA6iLqj98R5OF3tj8nJnbExr9CyS7q0DCxvBlbfKkzX6PTYPDUBBwrKceqy3NgLq75kkQDhbXgYFdseGdkSAMCQZ0OQJBIg18r1A6OF4LK8MCNb4hHz/gghhHgejwuqNm/ebPO9yMhImNblh4WFYd++fc2xLKc0thFkiawGRy5V2AxkBkYL4cf1xivfnzALqIC64Oj5wZ3B65KIyuAt4IkS4D/kKbz0aH9szCu0We/EYzMxISECWXmFZl3Uk0QCTE7qhNiOfLMTh4DhFGAUUj7/C2sm9EF8RBD2XyhDetZhfJHaDwzkm61vYLQQ746LwYKfzlBzTEIIIW7L44IqT9CYRpCGQOzopQqr7QwMBc/VSo3NZx+8WI5B3Ttg/8QPwfBmoUYLY6d1Wxmw9KQoqxmsXGk5GGDgzVHdsGBMD1wqVxhbMkiKZcbtwvSsw/g1IxkanR43a9UI4HpjxeO9UFWrMasrKq9WYfe5Uqd+Js3N07cmCSGE2EdBlRtqTCNI00DMdOSMIZARBfshNNAHkqIKm8/NzC3E9uliLNLrjVtwCpUWkqIKPNjTeuG+3RN70jKotTpcvF6Dad8cs3qNQqVFhUKF3uFBZq+3rXey9WJZtc11A4CspmWHKDd2u5YQQojnoKDKDTWmEaRpIGYYOWNq27RERMAXAVwW9DotKg9tgabyOgQjppvdVyKrxdtje0Kh1uHfsmpjZun/zl61uq3o6MReraouqOOxmWYnBw3F8Zm5hQ1qbOnoZ6JU6yBXtEy2iub2EUIIASiockuNmW3nKOjgspmQK1S4WVqMmq1vQVZwCgDg13MYOB3qivrFIgEOX7oBoX9beHt54ZtDl8yaeVrbVgz0sf+5vhwmCsuqsCE1Hmv+kJoFe2KRAJlp/Ro0dkfox0ZytNDq1qVYJMCBi+VoG8BtkeCF5vYRQggBPHj2X2vW0Nl2coUKBaVVkBRVQKfXY+kjMcaZeqaSRUL8euoKRr6wAOJ7++F6wSl4c3kQPDgL7Pb3AICxtcH5K5XQ6oBdZ65ixhCRsYeUYVbgQzGh2DFTjE+e7oMNqfHw8qrrdm5NkkgAnQ7oFsrHuj+kFlmuPGk51v5xO8gy/X4KrldBrri9pcfnsbFwTA+LmYKGdWfmFrZYHyua20cIIQSgTJXbctQI0uooGpEQG1LjMXnTEWNHc7FIgGlDRHjt5Rfx1/ZvAACJSclYueZT6HyDITMZIfPd30VY8nBPVKs0+HTfRcR2CMRDMaGY80BXFN+oMV739s/njM/nsZnYNk2MRTvPmAVNydFCzH2gK3R6PapqNWan/0wZOqVXq7QOa5IYAOLCg8xqxUyL3ltqPh7N7SOEEAJQUOXWbDWCtDmKRloGQI9vn70X/1XcDoImbzqMZx99HHk//wB+0lPYuPkDdAnlm80WHNe7A4RJUeDz2LhcoUCf8EBM33wM6UlR6CT0s1toXlhebRHshLfhIcSfg/JqFWrU9uuu5DVqLNxR1y6hfu3VpfJqML0YaBvAhcCXjZPFMquF8ba2RZtDY7ZrCSGEeB4KqlohezU8OdJypN1UWgRB3eL6Y/43fyA5phNktVoUXK+C0JeNziF+ZtfJFSos+OkMJid1woPyGrQN4ILHsdxSNOXtxbAIdPbMHgQA0OsBXytbkqZ4bCamJHVC34gg9OoYiC9yL5p3bb+17dk+0Mct5+PR3D5CCCEABVWtkqManvon8nhsJqIEvviykoEskyHK1o78l1WpcKCgHBMSwvHLqSvIlZZjxlCRzYaiYpEAkmKZ2Wv1u5/PGCpCskhgdQtQLBLgWFEFSuS1OFFUgYSoNpAUmT8vx+QUnbvOx3PXdRFCCGk+VKjeCtmq4dFr1Kj4MwtnDueZvf7WQ92wuF7NE3D7yL9coTIWiZdXq/DV5ARcldfi2K3gJjO3EJPEURZF4snRQswYEo3M3ELja0kiAd4e29Os+3lmbiHmjbJdZP72z+cQFxaIHGk51vwhRXpSlMX3ZjhFB9RlhjqH+KF3eBA6h/i5TeDirusihBDSPChT1QpZq+FRXf8XZTtXQl1aiHcv5uLVL35FQnRdw872gT4ouRUkGQrMDY5cqkCFQo1520+bbSmKRQKsSokzFoGbNhTl+7DA92HheLEMf/9bjtUpcWaF4yWyGrPu57bqrkyLzA3ZtTxpOdLFlkEVQKfoCCGEuDcKqloh0xqeff+U4ubhbajY/yWg1aCNQIh33/8AuaW1+OzA7a2++kGSQXpSFOZtO2WxNWfIaqUnRWHNXqlZQ9Ft0xLB9GLg0g2F1Wae3UPrtUOH9borUxzv20lTWw1FbZ2io/EwhBBC3AEFVa1U+0AfrBzfA2NGvYa/8vYDAEY88CBWrfkES/64YjEwuX6QZGBvzIytrFGADwt6AJKiCotmnqtS4uDNYFjcIymWNbguyzTAMrB1io7GwxBCCHEXVFPVigUH+qFv7xj4+vris88+w68/7wTTr43Vk4E8NhNx4UEY0aMt1j3VB5lp/TBjqAhW4h8z9bNGA6OF8OV4Y+H201abeW7MK4RWr7doCGqoy0qy07wTqKvTKr2ptPhMa6foHI2HMW0eSgghhDQ1ylS1cu+99x5mzZqFzp07A7B+MtAwYmZjXqFZZilJJMDoWOuDkg1Ms0aG4MZeM888aTlmDonG0nExeOPHU8a6L4VKi+/+LsLy8bGoVesgq1FBqdbhwMVy45ak4fk8NhP9I9s4PEVH42EIIYS4EwqqWjkej2cMqADrJwPTk6KwMa/QIrOUKy3HsSKZ3caVomA/bJuWaBbcSIoq7K7Jm8lAxzY8hy0G5AoV2gZwMbxriMX7DQmGaDwMIYQQd0Lbf61AdXU15s+fj+rqaofXGk4GmooLC7RaywQAS3aexeKHe9qcMxgh9LVoEeBoLItKo4NcoXLYYuBOWxDQeBhCCCHuhDJVbu7QoUN45plncOHCBZSWlmL9+vV2r7fW3dvWaTqgbluuskblVONKoR8bydFCq1tvYpEABy6Wo20At8m33mg8DCGEEHdCmSo3pVarMX/+fIjFYly4cAEdO3bEo48+2qB7Dd2998wehG3TEtFJ6Gv3el8Oy6msEZ/HxsIxtpt5ZuYWNsvWmyGAtJVlo3oqQgghzYkyVW7o0qVLGD9+PI4ePQoAmDBhAtasWYOgoKAGP8N0GLNcoXJ5RocB2G3m2VxbbzQehhBCiLugoMoNBQYG4vr16wgKCsInn3yCJ5544o6edycDf2011hT4snGyWGa1x1Vzb72ZBpCEEEJIS2Ho9Xp9Sy+itausrASfz4dcLkdAgGU38cY4evQo2rVrhw4dOjTqfmvBEACnMjqOGmuWyGpsBmqh1HiTEEKIm3P1728KqlygKYKqO+GKLuNyhQozsiVWi9EHRguxOiUOfB7bGLzR1hshhJDWxtW/v6lQvQVVVFTgu+++c+kzXdVl3FFjzSuVtQ1qm0AIIYTcLSioaiG7d+9GTEwMUlJSsH//fpc9tyFdxhvCUWPNi9erMTNbghJZjdNrJIQQQjwRBVXNrKamBi+99BLuu+8+XL58GZ07d4aPj+vqj1zVZdxRY02OtxfN2COEEEJMUFDVjCQSCfr27YuPP/4YAPDCCy/g+PHj6Nevn8s+w1Vdxq11ZjcQiwSQFMsAOJf9IoQQQjwZBVXNqKSkBOfOnUO7du3wyy+/YN26dfD1td+Y01n2giFnWh3Yaqxp2uDTgGbsEUIIIXT6zyWcOT3wxRdfYOzYsRAKrQc+ruDKVgdyhQpX5LW4WFZtbPCZmVsIhUprvGbP7EHoHOLnsvUTQgghzcHVp/+o+WczmzJlSpN/hiu7jBvuefeXczRjjxBCCLGDtv+awNWrV3H27NkWXYMrWx3QjD1CCCHEsVYTVL3zzjtITEwEj8dDYGCg1WuKiorw0EMP/X979x4WVbX/D/w9CMyM4IDI3WBAUQQBA00ESypQVOxIecyQk5cIlYOShX7NMhG70DHLyjoHew4HtLTUTmqagmjiBRBQgQqQ2xlBdIDCUBCVy3x+f/hz18gMCo4I8nk9zzxPe+219v6stXDPpz1rz6B///6wtLTE8uXL0dra2uFxL126hNDQUMhkMpiamiIsLAyNjY1djnPXrl1wd3fHjBkzcO3aw/N1A7f/SPPh1/ywMcSTvzmdMcYY+/96zcd/zc3NmDlzJnx8fJCQkNBuf1tbG4KCgmBtbY2MjAwolUrMmTMHBgYGeO+997QeNzQ0FEqlEqmpqWhpacH8+fOxYMECbNu2rdMxRkRECO0GDx6M2tpayOXyTh+np+Lf2GOMMca063UL1ZOSkrB06VLU19erlR84cADTpk3DxYsXYWVlBQCIj4/HihUr8Ouvv8LQsH0yUFRUBFdXV+Tk5GDMmDEAgOTkZEydOhVVVVWwtbW9q5huLXQDAJFIhBUrVmDNmjUQi8X30FPGGGOM3U/8MzVaZGZmwt3dXUioACAwMBBXrlxBQUGB1jampqZCQgUAAQEB0NPTQ1ZWltZz3bhxA1euXFF7AYC9vT2OHTuGuLg4TqgYY4yxPuahSaqqq6vVEioAwnZ1dbXWNpaWlmpl+vr6MDMz09oGAOLi4mBiYiK87OzsAADp6el4/PHH76UbjDHGGOulHmhS9frrr0MkEnX4Onv27IMMUaOVK1fi8uXLwuv8+fMAoJNbh4wxxhjrnR7oQvXo6GjMmzevwzpDhgy5q2NZW1sjOztbraympkbYp61NbW2tWllraysuXbqktQ0AiMVi/niPMcYYY2oeaFJlYWEBCwsLnRzLx8cH7777Lmpra4WP9FJTUyGTyeDq6qq1TX19PU6fPo3Ro0cDAH788UeoVCp4e3vrJC7GGGOM9Q29Zk1VZWUl8vLyUFlZiba2NuTl5SEvL0/4TqlJkybB1dUVL774IvLz85GSkoJVq1YhMjJSuKuUnZ2NESNG4MKFCwAAFxcXTJ48GeHh4cjOzkZ6ejoWL16MF1544a6f/GOMMcYYA3rR91StXr0amzdvFrY9PT0BAEeOHMGTTz6Jfv36Yd++fYiIiICPjw+MjIwwd+5crF27VmjT1NSE4uJitLT88QPAW7duxeLFi+Hv7w89PT3MmDEDn376afd1jDHGGGMPhV73PVU9ka6/54Ixxhhj9x9/TxVjjDHGWA/Uaz7+68lu3ey79SWgjDHGGOv5br1v6+pDO06qdKCurg4AhC8BZYwxxljv0dDQIPzc3L3gpEoHzMzMANx8QlEXk9IbXLlyBXZ2djh//nyfWkfWF/vNfeY+P8z6Yr+5z3/0mYjQ0NCgsyf+OanSAT29m0vTTExM+swf6C0ymazP9Rnom/3mPvcNfbHPQN/sN/f5Jl3eDOGF6owxxhhjOsBJFWOMMcaYDnBSpQNisRgxMTF96vcA+2Kfgb7Zb+5z39AX+wz0zX5zn+8f/vJPxhhjjDEd4DtVjDHGGGM6wEkVY4wxxpgOcFLFGGOMMaYDnFQxxhhjjOkAJ1V34d1334Wvry/69+8PU1NTjXUqKysRFBSE/v37w9LSEsuXL0dra2uHx7106RJCQ0Mhk8lgamqKsLAwNDY23oce3Lu0tDSIRCKNr5ycHK3tnnzyyXb1Fy1a1I2R3xsHB4d28b///vsdtrl+/ToiIyMxaNAgGBsbY8aMGaipqemmiO/duXPnEBYWBkdHR0ilUgwdOhQxMTFobm7usF1vm+vPP/8cDg4OkEgk8Pb2RnZ2dof1d+7ciREjRkAikcDd3R379+/vpkh1Iy4uDo899hgGDBgAS0tLBAcHo7i4uMM2SUlJ7eZUIpF0U8T3bs2aNe3iHzFiRIdtevs8a7pmiUQiREZGaqzfW+f42LFjeOaZZ2BrawuRSITdu3er7ScirF69GjY2NpBKpQgICEBpaekdj9vZ68LtOKm6C83NzZg5cyYiIiI07m9ra0NQUBCam5uRkZGBzZs3IykpCatXr+7wuKGhoSgoKEBqair27duHY8eOYcGCBfejC/fM19cXSqVS7fXyyy/D0dERY8aM6bBteHi4Wrt169Z1U9S6sXbtWrX4lyxZ0mH9V199FXv37sXOnTtx9OhRXLx4Ec8991w3RXvvzp49C5VKhU2bNqGgoAAbNmxAfHw83njjjTu27S1zvX37drz22muIiYnBmTNnMGrUKAQGBqK2tlZj/YyMDISEhCAsLAy5ubkIDg5GcHAwfvnll26OvOuOHj2KyMhInDx5EqmpqWhpacGkSZNw9erVDtvJZDK1Oa2oqOimiHVj5MiRavGfOHFCa92HYZ5zcnLU+puamgoAmDlzptY2vXGOr169ilGjRuHzzz/XuH/dunX49NNPER8fj6ysLBgZGSEwMBDXr1/XeszOXhc0InbXEhMTycTEpF35/v37SU9Pj6qrq4Wyf/3rXySTyejGjRsaj1VYWEgAKCcnRyg7cOAAiUQiunDhgs5j17Xm5maysLCgtWvXdljPz8+PXnnlle4J6j6Qy+W0YcOGu65fX19PBgYGtHPnTqGsqKiIAFBmZuZ9iLB7rFu3jhwdHTus05vmeuzYsRQZGSlst7W1ka2tLcXFxWms//zzz1NQUJBambe3Ny1cuPC+xnk/1dbWEgA6evSo1jrarnm9RUxMDI0aNequ6z+M8/zKK6/Q0KFDSaVSadzf2+eYiAgA7dq1S9hWqVRkbW1NH3zwgVBWX19PYrGYvv76a63H6ex1QRO+U6UDmZmZcHd3h5WVlVAWGBiIK1euoKCgQGsbU1NTtbs8AQEB0NPTQ1ZW1n2P+V59//33qKurw/z58+9Yd+vWrTA3N4ebmxtWrlyJpqambohQd95//30MGjQInp6e+OCDDzr8WPf06dNoaWlBQECAUDZixAjY29sjMzOzO8K9Ly5fviz8cHhHesNcNzc34/Tp02pzpKenh4CAAK1zlJmZqVYfuPlvvLfPKYA7zmtjYyPkcjns7Owwffp0rde0nqq0tBS2trYYMmQIQkNDUVlZqbXuwzbPzc3N+Oqrr/DSSy9BJBJprdfb5/h2CoUC1dXVanNpYmICb29vrXPZleuCJvyDyjpQXV2tllABELarq6u1trG0tFQr09fXh5mZmdY2PUlCQgICAwPxyCOPdFhv9uzZkMvlsLW1xU8//YQVK1aguLgY3333XTdFem+ioqLg5eUFMzMzZGRkYOXKlVAqlfjoo4801q+uroahoWG7tXdWVla9Yl41KSsrw8aNG7F+/foO6/WWuf7tt9/Q1tam8d/s2bNnNbbR9m+8t86pSqXC0qVLMX78eLi5uWmt5+zsjP/85z/w8PDA5cuXsX79evj6+qKgoOCO//Z7Am9vbyQlJcHZ2RlKpRKxsbF44okn8Msvv2DAgAHt6j9s87x7927U19dj3rx5Wuv09jnW5NZ8dWYuu3Jd0KTPJlWvv/46/vGPf3RYp6io6I6LGnu7roxDVVUVUlJSsGPHjjse/89rxNzd3WFjYwN/f3+Ul5dj6NChXQ/8HnSmz6+99ppQ5uHhAUNDQyxcuBBxcXG97iceujLXFy5cwOTJkzFz5kyEh4d32LYnzjXTLDIyEr/88kuH64sAwMfHBz4+PsK2r68vXFxcsGnTJrz99tv3O8x7NmXKFOG/PTw84O3tDblcjh07diAsLOwBRtY9EhISMGXKFNja2mqt09vnuKfps0lVdHR0h9k7AAwZMuSujmVtbd3uCYFbT3tZW1trbXP74rfW1lZcunRJa5v7oSvjkJiYiEGDBuEvf/lLp8/n7e0N4Obdjwf1Rnsvc+/t7Y3W1lacO3cOzs7O7fZbW1ujubkZ9fX1anerampqunVeNelsvy9evIinnnoKvr6++OKLLzp9vp4w15qYm5ujX79+7Z7I7GiOrK2tO1W/J1u8eLHwYExn70QYGBjA09MTZWVl9ym6+8vU1BTDhw/XGv/DNM8VFRU4dOhQp+8U9/Y5Bv54362pqYGNjY1QXlNTg0cffVRjm65cFzTq3HKwvu1OC9VramqEsk2bNpFMJqPr169rPNatheqnTp0SylJSUnr8QnWVSkWOjo4UHR3dpfYnTpwgAJSfn6/jyLrHV199RXp6enTp0iWN+28tVP/222+FsrNnz/a6hepVVVU0bNgweuGFF6i1tbVLx+jJcz127FhavHixsN3W1kaDBw/ucKH6tGnT1Mp8fHx61QJmlUpFkZGRZGtrSyUlJV06RmtrKzk7O9Orr76q4+i6R0NDAw0cOJA++eQTjfsfhnm+JSYmhqytramlpaVT7XrjHEPLQvX169cLZZcvX76rheqduS5ojKVzofdNFRUVlJubS7GxsWRsbEy5ubmUm5tLDQ0NRHTzj9DNzY0mTZpEeXl5lJycTBYWFrRy5UrhGFlZWeTs7ExVVVVC2eTJk8nT05OysrLoxIkTNGzYMAoJCen2/nXGoUOHCAAVFRW121dVVUXOzs6UlZVFRERlZWW0du1aOnXqFCkUCtqzZw8NGTKEJkyY0N1hd0lGRgZt2LCB8vLyqLy8nL766iuysLCgOXPmCHVu7zMR0aJFi8je3p5+/PFHOnXqFPn4+JCPj8+D6EKXVFVVkZOTE/n7+1NVVRUplUrh9ec6vXmuv/nmGxKLxZSUlESFhYW0YMECMjU1FZ7gffHFF+n1118X6qenp5O+vj6tX7+eioqKKCYmhgwMDOjnn39+UF3otIiICDIxMaG0tDS1OW1qahLq3N7v2NhYSklJofLycjp9+jS98MILJJFIqKCg4EF0odOio6MpLS2NFAoFpaenU0BAAJmbm1NtbS0RPZzzTHQzGbC3t6cVK1a02/ewzHFDQ4PwXgyAPvroI8rNzaWKigoiInr//ffJ1NSU9uzZQz/99BNNnz6dHB0d6dq1a8Ixnn76adq4caOwfafrwt3gpOouzJ07lwC0ex05ckSoc+7cOZoyZQpJpVIyNzen6Ohotf9DOHLkCAEghUIhlNXV1VFISAgZGxuTTCaj+fPnC4laTxUSEkK+vr4a9ykUCrVxqayspAkTJpCZmRmJxWJycnKi5cuX0+XLl7sx4q47ffo0eXt7k4mJCUkkEnJxcaH33ntP7e7j7X0mIrp27Rr9/e9/p4EDB1L//v3p2WefVUtIerrExESNf+9/vrH9MMz1xo0byd7engwNDWns2LF08uRJYZ+fnx/NnTtXrf6OHTto+PDhZGhoSCNHjqQffvihmyO+N9rmNDExUahze7+XLl0qjJGVlRVNnTqVzpw50/3Bd9GsWbPIxsaGDA0NafDgwTRr1iwqKysT9j+M80x081MPAFRcXNxu38Myx7feU29/3eqbSqWit956i6ysrEgsFpO/v3+78ZDL5RQTE6NW1tF14W6IiIju/sNCxhhjjDGmCX9PFWOMMcaYDnBSxRhjjDGmA5xUMcYYY4zpACdVjDHGGGM6wEkVY4wxxpgOcFLFGGOMMaYDnFQxxhhjjOkAJ1WMMcYYYzrASRVjrE8QiUTYvXv3gw6jU3pjzIz1ZZxUMdYLZGZmol+/fggKCnrQofR4a9as0fhL9EqlElOmTLnv5+9tiZC28WKMdR4nVYz1AgkJCViyZAmOHTuGixcv3tdzERFaW1vv6zkeBGtra4jF4gcdBmPsIcZJFWM9XGNjI7Zv346IiAgEBQUhKSlJ2Dd79mzMmjVLrX5LSwvMzc2xZcsWAIBKpUJcXBwcHR0hlUoxatQofPvtt0L9tLQ0iEQiHDhwAKNHj4ZYLMaJEydQXl6O6dOnw8rKCsbGxnjsscdw6NAhtXMplUoEBQVBKpXC0dER27Ztg4ODAz7++GOhTn19PV5++WVYWFhAJpPh6aefRn5+fod9Pn/+PJ5//nmYmprCzMwM06dPx7lz59RiHjt2LIyMjGBqaorx48ejoqICSUlJiI2NRX5+PkQiEUQikTBef76DdO7cOYhEIuzYsQNPPPEEpFIpHnvsMZSUlCAnJwdjxoyBsbExpkyZgl9//VU4b05ODiZOnAhzc3OYmJjAz88PZ86cEfY7ODgAAJ599lmIRCJhGwD27NkDLy8vSCQSDBkyBLGxsWrJa2lpKSZMmACJRAJXV1ekpqZ2OEYAcOPGDURFRcHS0hISiQSPP/44cnJyhP1JSUkwNTVVa7N7926IRCJhv7bxqq+vx8KFC2FlZQWJRAI3Nzfs27dPOM5///tfjBw5EmKxGA4ODvjwww/VzuPg4IB33nkHc+bMgbGxMeRyOb7//nv8+uuvmD59OoyNjeHh4YFTp06ptTtx4oQwJ3Z2doiKisLVq1fvOBaM9Qhd/oloxli3SEhIoDFjxhAR0d69e2no0KGkUqmIiGjfvn0klUqpoaFBqL93716SSqV05coVIiJ65513aMSIEZScnEzl5eWUmJhIYrGY0tLSiOiPX3v38PCggwcPUllZGdXV1VFeXh7Fx8fTzz//TCUlJbRq1SqSSCRUUVEhnCsgIIAeffRROnnyJJ0+fZr8/PxIKpXShg0b1Oo888wzlJOTQyUlJRQdHU2DBg2iuro6jf1tbm4mFxcXeumll+inn36iwsJCmj17Njk7O9ONGzeopaWFTExMaNmyZVRWVkaFhYWUlJREFRUV1NTURNHR0TRy5EhSKpWkVCqpqamJiIgA0K5du4iISKFQEABhXAoLC2ncuHE0evRoevLJJ+nEiRN05swZcnJyokWLFgmxHT58mL788ksqKiqiwsJCCgsLIysrK2Gsa2trCQAlJiaSUqmk2tpaIiI6duwYyWQySkpKovLycjp48CA5ODjQmjVriIiora2N3NzcyN/fn/Ly8ujo0aPk6empFrMmUVFRZGtrS/v376eCggKaO3cuDRw4UBjbxMREMjExUWuza9cuunXp1zZebW1tNG7cOBo5ciQdPHiQysvLae/evbR//34iIjp16hTp6enR2rVrqbi4mBITE0kqlVJiYqJwHrlcTmZmZhQfH08lJSUUERFBMpmMJk+eTDt27KDi4mIKDg4mFxcX4e+5rKyMjIyMaMOGDVRSUkLp6enk6elJ8+bN0zoGjPUknFQx1sP5+vrSxx9/TERELS0tZG5uTkeOHFHb3rJli1A/JCSEZs2aRURE169fp/79+1NGRobaMcPCwigkJISI/kiqdu/efcdYRo4cSRs3biQioqKiIgJAOTk5wv7S0lICICRVx48fJ5lMRtevX1c7ztChQ2nTpk0az/Hll1+Ss7Oz8EZLRHTjxg2SSqWUkpJCdXV1BEBICm8XExNDo0aNaleuKan697//Lez/+uuvCQAdPnxYKIuLiyNnZ2et49HW1kYDBgygvXv3ajzPLf7+/vTee++166eNjQ0REaWkpJC+vj5duHBB2H/gwIEOk6rGxkYyMDCgrVu3CmXNzc1ka2tL69atI6I7J1VEmscrJSWF9PT0qLi4WOO5Z8+eTRMnTlQrW758Obm6ugrbcrmc/va3vwnbSqWSANBbb70llGVmZhIAUiqVRHTz73LBggVqxz1+/Djp6enRtWvXNMbCWE/CH/8x1oMVFxcjOzsbISEhAAB9fX3MmjULCQkJwvbzzz+PrVu3AgCuXr2KPXv2IDQ0FABQVlaGpqYmTJw4EcbGxsJry5YtKC8vVzvXmDFj1LYbGxuxbNkyuLi4wNTUFMbGxigqKkJlZaUQm76+Pry8vIQ2Tk5OGDhwoLCdn5+PxsZGDBo0SO38CoWi3fn/3KasrAwDBgwQ6puZmeH69esoLy+HmZkZ5s2bh8DAQDzzzDP45JNPoFQquzS+Hh4ewn9bWVkBANzd3dXKamtrhe2amhqEh4dj2LBhMDExgUwmQ2NjozAm2uTn52Pt2rVqYxAeHg6lUommpiYUFRXBzs4Otra2QhsfH58Oj1leXo6WlhaMHz9eKDMwMMDYsWNRVFR0dwOgRV5eHh555BEMHz5c4/6ioiK18wLA+PHjUVpaira2NqHsbsYXgDDG+fn5SEpKUhunwMBAqFQqKBSKe+oTY91B/0EHwBjTLiEhAa2trWpvtkQEsViMzz77DCYmJggNDYWfnx9qa2uRmpoKqVSKyZMnA7iZGAHADz/8gMGDB6sd+/ZF20ZGRmrby5YtQ2pqKtavXw8nJydIpVL89a9/RXNz813H39jYCBsbG6SlpbXbd/tanz+3GT16tJAo/pmFhQUAIDExEVFRUUhOTsb27duxatUqpKamYty4cXcdG3AzCbnl1jqj28tUKpWwPXfuXNTV1eGTTz6BXC6HWCyGj4/PHceksbERsbGxeO6559rtk0gknYq5M/T09EBEamUtLS13bCeVSnVy/rsZXwDCGDc2NmLhwoWIiopqdyx7e3udxMTY/cRJFWM9VGtrK7Zs2YIPP/wQkyZNUtsXHByMr7/+GosWLYKvry/s7Oywfft2HDhwADNnzhTeuFxdXSEWi1FZWQk/P79OnT89PR3z5s3Ds88+C+DmG96fF4s7OzujtbUVubm5GD16NICbd8Z+//13oY6Xlxeqq6uhr6+vtmi7I15eXti+fTssLS0hk8m01vP09ISnpydWrlwJHx8fbNu2DePGjYOhoaHa3RJdSk9Pxz//+U9MnToVwM0F9b/99ptaHQMDg3bn9/LyQnFxMZycnDQe18XFBefPn4dSqYSNjQ0A4OTJkx3GMnToUBgaGiI9PR1yuRzAzYQpJycHS5cuBXAzCW1oaMDVq1eFpDkvL0/tOJrGy8PDA1VVVSgpKdF4t8rFxQXp6elqZenp6Rg+fDj69evXYdwd8fLyQmFhodZxYqyn44//GOuh9u3bh99//x1hYWFwc3NTe82YMUP4CBC4+RRgfHw8UlNThY/+AGDAgAFYtmwZXn31VWzevBnl5eU4c+YMNm7ciM2bN3d4/mHDhuG7775DXl4e8vPzMXv2bLW7NiNGjEBAQAAWLFiA7Oxs5ObmYsGCBZBKpcIdiICAAPj4+CA4OBgHDx7EuXPnkJGRgTfffLPdU1+3hIaGwtzcHNOnT8fx48ehUCiQlpaGqKgoVFVVQaFQYOXKlcjMzERFRQUOHjyI0tJSuLi4ALj51JlCoUBeXh5+++033Lhxo8tzoGlMvvzySxQVFSErKwuhoaHt7uo4ODjg8OHDqK6uFhLM1atXY8uWLYiNjUVBQQGKiorwzTffYNWqVcI4DR8+HHPnzkV+fj6OHz+ON998s8NYjIyMEBERgeXLlyM5ORmFhYUIDw9HU1MTwsLCAADe3t7o378/3njjDZSXl2Pbtm1qT4/eivf28fLz88OECRMwY8YMpKamQqFQ4MCBA0hOTgYAREdH4/Dhw3j77bdRUlKCzZs347PPPsOyZcvuaXxXrFiBjIwMLF68GHl5eSgtLcWePXuwePHiezouY93mQS/qYoxpNm3aNJo6darGfVlZWQSA8vPziYiosLCQAJBcLldb4E1EpFKp6OOPPyZnZ2cyMDAgCwsLCgwMpKNHjxLRHwvVf//9d7V2CoWCnnrqKZJKpWRnZ0efffYZ+fn50SuvvCLUuXjxIk2ZMoXEYjHJ5XLatm0bWVpaUnx8vFDnypUrtGTJErK1tSUDAwOys7Oj0NBQqqys1Np3pVJJc+bMIXNzcxKLxTRkyBAKDw+ny5cvU3V1NQUHB5ONjQ0ZGhqSXC6n1atXU1tbGxHdXJw/Y8YMMjU1FZ7EI9K8UD03N1c4p6ZxuH2h95kzZ2jMmDEkkUho2LBhtHPnTpLL5WpPO37//ffk5ORE+vr6JJfLhfLk5GTy9fUlqVRKMpmMxo4dS1988YWwv7i4mB5//HEyNDSk4cOHU3Jy8h2f/rt27RotWbJEGKfx48dTdna2Wp1du3aRk5MTSaVSmjZtGn3xxRdqC9W1jVddXR3Nnz+fBg0aRBKJhNzc3Gjfvn1Cu2+//ZZcXV3JwMCA7O3t6YMPPlA77+3jcvscEGmeh+zsbJo4cSIZGxuTkZEReXh40Lvvvqt1DBjrSUREt33gzhhjXVRVVQU7OzscOnQI/v7+DzocxhjrVpxUMca67Mcff0RjYyPc3d2hVCrxf//3f7hw4QJKSkrUFiQzxlhfwAvVGWNd1tLSgjfeeAP/+9//MGDAAPj6+mLr1q2cUDHG+iS+U8UYY4wxpgP89B9jjDHGmA5wUsUYY4wxpgOcVDHGGGOM6QAnVYwxxhhjOsBJFWOMMcaYDnBSxRhjjDGmA5xUMcYYY4zpACdVjDHGGGM68P8AM91nAKgwPLUAAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "forest_preds_y_mcmc = bart_model.y_hat_test[:,bart_model.num_gfr:]\n", + "y_avg_mcmc = np.squeeze(forest_preds_y_mcmc).mean(axis = 1, keepdims = True)\n", + "y_df_mcmc = pd.DataFrame(np.concatenate((np.expand_dims(y_test,1), y_avg_mcmc), axis = 1), columns=[\"True outcome\", \"Average estimated outcome\"])\n", + "sns.scatterplot(data=y_df_mcmc, x=\"Average estimated outcome\", y=\"True outcome\")\n", + "plt.axline((0, 0), slope=1, color=\"black\", linestyle=(0, (3,3)))\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAGwCAYAAABVdURTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuNSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/xnp5ZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA8VElEQVR4nO3de3hU1aH//88EciWZCZAm4ZJAFL4IahRBIGhBKxbxBpb2eDieghU9VaFC8XuUWG2/1Z8NSrUqckRrFWulKCpoqUopyE3iBQS52CIgBQoEDMLkRjIx2b8/ODPN7EzmlslcNu/X88zzMHv2zKy9wuz92WutvbbNMAxDAAAAFpEU6wIAAABEEuEGAABYCuEGAABYCuEGAABYCuEGAABYCuEGAABYCuEGAABYSudYFyDampubdfjwYWVlZclms8W6OAAAIAiGYai6ulo9e/ZUUpL/tpkzLtwcPnxYBQUFsS4GAAAIw8GDB9W7d2+/65xx4SYrK0vS6cqx2+0xLg0AAAhGVVWVCgoKPMdxf864cOPuirLb7YQbAAASTDBDShhQDAAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALIVwAwAALOWMu/1CR3HWuVRZ41JVfaPs6cnK6ZIiR0ZKrIsFAMAZh3ATAYdPntK9b2zT+t2VnmWj+udozsRi9cxOj2HJAAA489At1U7OOlerYCNJ63ZXavYb2+Ssc8WoZAAAnJkIN+1UWeNqFWzc1u2uVGUN4QYAgGgi3LRTVX2j39erA7wOAAAii3DTTva0ZL+vZwV4HQAARBbhpp1yMlM0qn+Oz9dG9c9RTiZXTAEAEE2Em3ZyZKRozsTiVgFnVP8cPTKxmMvBAQCIMi4Fj4Ce2emaN2mwKmtcqq5vVFZasnIymecGAIBYINxEiCODMAMAQDygWwoAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFhKTMPNM888o+LiYtntdtntdpWUlOjdd99tc/2FCxfKZrN5PdLS0qJYYgAAEO86x/LLe/furTlz5qh///4yDEMvvfSSxo8fry1btujcc8/1+R673a5du3Z5nttstmgVFwAAJICYhpvrrrvO6/nDDz+sZ555Rh9++GGb4cZmsyk/Pz/o72hoaFBDQ4PneVVVVXiFBQAACSFuxtw0NTVp8eLFqq2tVUlJSZvr1dTUqE+fPiooKND48eO1c+dOv59bVlYmh8PheRQUFES66AAAII7YDMMwYlmA7du3q6SkRPX19crMzNSiRYt09dVX+1y3vLxcu3fvVnFxsZxOp379619r3bp12rlzp3r37u3zPb5abgoKCuR0OmW32ztkmwAAQGRVVVXJ4XAEdfyOebhxuVw6cOCAnE6nXn/9dT3//PNau3atBg0aFPC9jY2NGjhwoCZNmqSHHnooqO8LpXIAAEB8COX4HdMxN5KUkpKifv36SZKGDBmiTz75RE8++aSeffbZgO9NTk7W4MGDtWfPno4uJgAASBBxM+bGrbm52asbyZ+mpiZt375dPXr06OBSAQCARBHTlpvS0lKNGzdOhYWFqq6u1qJFi7RmzRqtWLFCkjR58mT16tVLZWVlkqQHH3xQI0aMUL9+/XTy5EnNnTtX+/fv16233hrLzQAAAHEkpuHm2LFjmjx5so4cOSKHw6Hi4mKtWLFCV155pSTpwIEDSkr6V+PSiRMndNttt6miokJdu3bVkCFDtHHjxqDG5wAAgDNDzAcURxsDigEASDyhHL/jbswNAABAexBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApRBuAACApcQ03DzzzDMqLi6W3W6X3W5XSUmJ3n33Xb/vWbJkic455xylpaXp/PPP1zvvvBOl0gIAgEQQ03DTu3dvzZkzR5s3b9amTZv0ne98R+PHj9fOnTt9rr9x40ZNmjRJU6dO1ZYtWzRhwgRNmDBBO3bsiHLJAQBAvLIZhmHEuhAtdevWTXPnztXUqVNbvXbjjTeqtrZWy5cv9ywbMWKELrzwQi1YsCCoz6+qqpLD4ZDT6ZTdbo9YuQEAQMcJ5fgdN2NumpqatHjxYtXW1qqkpMTnOuXl5RozZozXsrFjx6q8vLzNz21oaFBVVZXXAwAAWFfMw8327duVmZmp1NRU3X777Vq6dKkGDRrkc92Kigrl5eV5LcvLy1NFRUWbn19WViaHw+F5FBQURLT8AAAgvsQ83AwYMEBbt27VRx99pDvuuENTpkzR559/HrHPLy0tldPp9DwOHjwYsc8GAADxp3OsC5CSkqJ+/fpJkoYMGaJPPvlETz75pJ599tlW6+bn5+vo0aNey44ePar8/Pw2Pz81NVWpqamRLTQAAIhbMW+5MWtublZDQ4PP10pKSrRq1SqvZStXrmxzjA4AADjzxLTlprS0VOPGjVNhYaGqq6u1aNEirVmzRitWrJAkTZ48Wb169VJZWZkkacaMGRo9erQee+wxXXPNNVq8eLE2bdqk5557LpabAQAA4khMw82xY8c0efJkHTlyRA6HQ8XFxVqxYoWuvPJKSdKBAweUlPSvxqWRI0dq0aJFuv/++3Xfffepf//+WrZsmc4777xYbQIAAIgzcTfPTUdjnhsAABJPQs5zAwAAEAmEGwAAYCmEGwAAYCmEGwAAYCmEGwAAYCmEGwAAYCmEGwAAYCmEGwAAYCmEGwAAYCmEGwAAYCmEGwAAYCmEGwAAYCmEGwAAYCmEGwAAYCmEGwAAYCmEGwAAYCmEGwAAYCmEGwAAYCmEGwAAYCmEGwAAYCmEGwAAYCmEGwAAYCmEGwAAYCmEGwAAYCmEGwAAYCmEGwAAYCmEGwAAYCmEGwAAYCmEGwAAYCmEGwAAYCmEGwAAYCmEGwAAYCmEGwAAYCmEGwAAYCmEGwAAYCmEGwAAYCmEGwAAYCmEGwAAYCmEGwAAYCmEGwAAYCmEGwAAYCkxDTdlZWW6+OKLlZWVpdzcXE2YMEG7du3y+56FCxfKZrN5PdLS0qJUYgAAEO9iGm7Wrl2radOm6cMPP9TKlSvV2Nio7373u6qtrfX7PrvdriNHjnge+/fvj1KJAQBAvOscyy9/7733vJ4vXLhQubm52rx5s0aNGtXm+2w2m/Lz84P6joaGBjU0NHieV1VVhVdYAACQEOJqzI3T6ZQkdevWze96NTU16tOnjwoKCjR+/Hjt3LmzzXXLysrkcDg8j4KCgoiWGQAAxBebYRhGrAshSc3Nzbr++ut18uRJbdiwoc31ysvLtXv3bhUXF8vpdOrXv/611q1bp507d6p3796t1vfVclNQUCCn0ym73d4h2wIAACKrqqpKDocjqON33ISbO+64Q++++642bNjgM6S0pbGxUQMHDtSkSZP00EMPBVw/lMoBAADxIZTjd0zH3LhNnz5dy5cv17p160IKNpKUnJyswYMHa8+ePR1UOgAAkEhiOubGMAxNnz5dS5cu1erVq1VUVBTyZzQ1NWn79u3q0aNHB5QQAAAkmpi23EybNk2LFi3SW2+9paysLFVUVEiSHA6H0tPTJUmTJ09Wr169VFZWJkl68MEHNWLECPXr108nT57U3LlztX//ft16660x2w4AABA/YhpunnnmGUnSZZdd5rX8xRdf1M033yxJOnDggJKS/tXAdOLECd12222qqKhQ165dNWTIEG3cuFGDBg2KVrEBAEAci5sBxdHCgGIAABJPKMfvuJrnBgAAoL0INwAAwFIINwAAwFIINwAAwFIINwAAwFIINwAAwFIINwAAwFIINwAAwFLaNUPx559/rgMHDsjlcnktv/7669tVKAAAgHCFFW6+/PJL3XDDDdq+fbtsNpvckxzbbDZJp29mCQAAEAthdUvNmDFDRUVFOnbsmDIyMrRz506tW7dOQ4cO1Zo1ayJcRAAAgOCF1XJTXl6u1atXKycnR0lJSUpKStKll16qsrIy3XXXXdqyZUukywkAABCUsFpumpqalJWVJUnKycnR4cOHJUl9+vTRrl27Ilc6AACAEIXVcnPeeefps88+U1FRkYYPH65HH31UKSkpeu6553TWWWdFuowAAABBCyvc3H///aqtrZUkPfjgg7r22mv17W9/W927d9err74a0QICAACEwma4L3Vqp6+//lpdu3b1XDEVr6qqquRwOOR0OmW322NdHAAAEIRQjt/tmuempW7dukXqowAAAMIWVripr6/XvHnz9P777+vYsWNqbm72ev3TTz+NSOEAAABCFVa4mTp1qv7yl7/o+9//voYNGxb3XVEAAODMEVa4Wb58ud555x1dcsklkS4PAABAu4Q1z02vXr0889wAAADEk7DCzWOPPaZ7771X+/fvj3R5AAAA2iWsbqmhQ4eqvr5eZ511ljIyMpScnOz1+tdffx2RwgEAAIQqrHAzadIkHTp0SL/61a+Ul5fHgGIAABA3wgo3GzduVHl5uS644IJIlwcAAKBdwhpzc8455+jUqVORLgsAAEC7hRVu5syZo7vvvltr1qzR8ePHVVVV5fUAAACIlbDuLZWUdDoTmcfaGIYhm82mpqamyJSuA3BvKQAAEk+H31vq/fffD6tgAAAAHS2scDN69OhIlwMAACAiwgo327Zt87ncZrMpLS1NhYWFSk1NbVfBAAAAwhFWuLnwwgv9zm2TnJysG2+8Uc8++6zS0tLCLhwAAECowrpaaunSperfv7+ee+45bd26VVu3btVzzz2nAQMGaNGiRfrd736n1atX6/777490eQEAAPwKq+Xm4Ycf1pNPPqmxY8d6lp1//vnq3bu3HnjgAX388cfq0qWL7r77bv3617+OWGEBAAACCavlZvv27erTp0+r5X369NH27dslne66OnLkSPtKBwAAEKKwZyieM2eOXC6XZ1ljY6PmzJmjc845R5J06NAh5eXlRaaUAAAAQQqrW2r+/Pm6/vrr1bt3bxUXF0s63ZrT1NSk5cuXS5K+/PJL3XnnnZErKQAAQBDCmqFYkqqrq/XKK6/oiy++kCQNGDBA//Ef/6GsrKyIFjDSmKEYAIDE0+EzFEtSVlaWbr/99nDfDgAA0CGCDjdvv/22xo0bp+TkZL399tt+173++uvbXTAAAIBwBN0tlZSUpIqKCuXm5npunOnzA0O4cWZZWZnefPNN/f3vf1d6erpGjhypRx55RAMGDPD7viVLluiBBx7QP/7xD/Xv31+PPPKIrr766qC+k24pAAASTyjH76CvlmpublZubq7n3209Qrkj+Nq1azVt2jR9+OGHWrlypRobG/Xd735XtbW1bb5n48aNmjRpkqZOnaotW7ZowoQJmjBhgnbs2BH09wIAAOsKaUBxeXm5jh8/rmuvvdaz7Pe//71+8YtfqLa2VhMmTNC8efPCvq/UV199pdzcXK1du1ajRo3yuc6NN96o2tpaz1VZkjRixAhdeOGFWrBgQcDvoOUGAIDE0yEtN5L04IMPaufOnZ7n27dv19SpUzVmzBjNnj1bf/rTn1RWVhZeqSU5nU5JUrdu3dpcp7y8XGPGjPFaNnbsWJWXl/tcv6GhQVVVVV4PAABgXSGFm61bt+qKK67wPF+8eLGGDx+u3/72t5o1a5aeeuopvfbaa2EVpLm5WTNnztQll1yi8847r831KioqWk0OmJeXp4qKCp/rl5WVyeFweB4FBQVhlQ8AACSGkMLNiRMnvILF2rVrNW7cOM/ziy++WAcPHgyrINOmTdOOHTu0ePHisN7fltLSUjmdTs8j3PIBAIDEEFK4ycvL0759+yRJLpdLn376qUaMGOF5vbq6WsnJySEXYvr06Vq+fLnef/999e7d2++6+fn5Onr0qNeyo0ePKj8/3+f6qampstvtXg8AAGBdIYWbq6++WrNnz9b69etVWlqqjIwMffvb3/a8vm3bNp199tlBf55hGJo+fbqWLl2q1atXq6ioKOB7SkpKtGrVKq9lK1euVElJSfAbAgAALCukGYofeughfe9739Po0aOVmZmpl156SSkpKZ7XX3jhBX33u98N+vOmTZumRYsW6a233lJWVpZn3IzD4VB6erokafLkyerVq5dnoPKMGTM0evRoPfbYY7rmmmu0ePFibdq0Sc8991womwIAACwqrHtLOZ1OZWZmqlOnTl7Lv/76a2VmZnoFHr9fbrP5XP7iiy/q5ptvliRddtll6tu3rxYuXOh5fcmSJbr//vs9k/g9+uijTOIHAICFhXL8DvvGmYmKcAMAQOLpsHluAAAA4h3hBgAAWArhBgAAWArhBgAAWArhBgAAWArhBgAAWArhBgAAWArhBgAAWArhBgAAWArhBgAAWArhBgAAWArhBgAAWArhBgAAWArhBgAAWErnWBfAqpx1LlXWuFRV3yh7erJyuqTIkZES62IBAGB5hJsOcPjkKd37xjat313pWTaqf47mTCxWz+z0GJYMAADro1sqwpx1rlbBRpLW7a7U7De2yVnnilHJAAA4MxBuIqyyxtUq2Lit212pyhrCDQAAHYlwE2FV9Y1+X68O8DoAAGgfwk2E2dOS/b6eFeB1AADQPoSbCMvJTNGo/jk+XxvVP0c5mVwxBQBARyLcRJgjI0VzJha3Cjij+ufokYnFXA4OAEAH41LwDtAzO13zJg1WZY1L1fWNykpLVk4m89wAABANhJsO4sggzAAAEAt0SwEAAEsh3AAAAEsh3AAAAEsh3AAAAEsh3AAAAEsh3AAAAEsh3AAAAEsh3AAAAEsh3AAAAEsh3AAAAEsh3AAAAEsh3AAAAEsh3AAAAEsh3AAAAEsh3AAAAEsh3AAAAEsh3AAAAEuJabhZt26drrvuOvXs2VM2m03Lli3zu/6aNWtks9laPSoqKqJTYAAAEPdiGm5qa2t1wQUXaP78+SG9b9euXTpy5IjnkZub20ElBAAAiaZzLL983LhxGjduXMjvy83NVXZ2dlDrNjQ0qKGhwfO8qqoq5O8DAACJIyHH3Fx44YXq0aOHrrzySn3wwQd+1y0rK5PD4fA8CgoKolRKAAAQCwkVbnr06KEFCxbojTfe0BtvvKGCggJddtll+vTTT9t8T2lpqZxOp+dx8ODBKJYYAABEW0y7pUI1YMAADRgwwPN85MiR2rt3r37zm9/o5Zdf9vme1NRUpaamRquIAAAgxhKq5caXYcOGac+ePbEuBgAAiBMJH262bt2qHj16xLoYAAAgTsS0W6qmpsar1WXfvn3aunWrunXrpsLCQpWWlurQoUP6/e9/L0l64oknVFRUpHPPPVf19fV6/vnntXr1av3lL3+J1SYAAIA4E9Nws2nTJl1++eWe57NmzZIkTZkyRQsXLtSRI0d04MABz+sul0t33323Dh06pIyMDBUXF+uvf/2r12cAAIAzm80wDCPWhYimqqoqORwOOZ1O2e32WBcHAAAEIZTjd8KPuQEAAGiJcAMAACyFcAMAACyFcAMAACwloWYoTmTOOpcqa1yqqm+UPT1ZOV1S5MhIiXWxAACwHMJNFBw+eUr3vrFN63dXepaN6p+jOROL1TM7PYYlAwDAeuiW6mDOOlerYCNJ63ZXavYb2+Ssc8WoZAAAWBPhpoNV1rhaBRu3dbsrVVlDuAEAIJIINx2sqr7R7+vVAV4HAAChIdx0MHtast/XswK8DgAAQkO46WA5mSka1T/H52uj+ucoJ5MrpgAAiCTCTQdzZKRozsTiVgFnVP8cPTKxmMvBAQCIMC4Fj4Ke2emaN2mwKmtcqq5vVFZasnIymecGAICOQLiJEkcGYQYAgGigWwoAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK4QYAAFgK95aKEWedS5U1LlXVN8qenqycLtx7CgCASCDcxMDhk6d07xvbtH53pWfZqP45mjOxWD2z02NYMgAAEh/dUlHmrHO1CjaStG53pWa/sU3OOleMSgYAgDUQbqKsssbVKti4rdtdqcoawg0AAO1BuImyqvpGv69XB3gdAAD4R7iJMntast/XswK8DgAA/CPcRFlOZopG9c/x+dqo/jnKyeSKKQAA2oNwE2WOjBTNmVjcKuCM6p+jRyYWczk4AADtxKXgMdAzO13zJg1WZY1L1fWNykpLVk4m89wAABAJhJsYcWQQZgAA6Ah0SwEAAEsh3AAAAEsh3AAAAEsh3AAAAEsh3AAAAEuJabhZt26drrvuOvXs2VM2m03Lli0L+J41a9booosuUmpqqvr166eFCxd2eDkBAEDiiGm4qa2t1QUXXKD58+cHtf6+fft0zTXX6PLLL9fWrVs1c+ZM3XrrrVqxYkUHlxQAACSKmM5zM27cOI0bNy7o9RcsWKCioiI99thjkqSBAwdqw4YN+s1vfqOxY8d2VDEBAEACSagxN+Xl5RozZozXsrFjx6q8vLzN9zQ0NKiqqsrrAQAArCuhwk1FRYXy8vK8luXl5amqqkqnTp3y+Z6ysjI5HA7Po6CgIBpFBQAAMZJQ4SYcpaWlcjqdnsfBgwdjXSQAANCBEureUvn5+Tp69KjXsqNHj8putys9Pd3ne1JTU5WamhqN4gEAgDiQUOGmpKRE77zzjteylStXqqSkJEYlihxnnUuVNS5V1TfKnp6snC7cWBMAgHDENNzU1NRoz549nuf79u3T1q1b1a1bNxUWFqq0tFSHDh3S73//e0nS7bffrqefflr33HOPbrnlFq1evVqvvfaa/vznP8dqEyLi8MlTuveNbVq/u9KzbFT/HM2ZWKye2b5bpAAAgG8xDTebNm3S5Zdf7nk+a9YsSdKUKVO0cOFCHTlyRAcOHPC8XlRUpD//+c/66U9/qieffFK9e/fW888/n9CXgTvrXK2CjSRt2n9Ca7/4SkP7dFVNwze05gAAECSbYRhGrAsRTVVVVXI4HHI6nbLb7bEujvYeq9EVj6/1WpaR0klPTRqsFz/Ypw/2HPcspzUHAHCmCuX4bfmrpeJdVX1jq2W3XFrUKthI0rrdlZr9xjY561zRKh4AAAmHcBNj9rTkVssGF2S3CjZu63ZXqrKGcAMAQFsINzGWk5miUf1zvJY1fNPs9z3VPlp7AADAaYSbGHNkpGjOxGKvgJPa2f+fJctHaw8AADgtoea5saqe2emaN2mwKmtcqq5vVNeM060560xXUEmnBxXnZHLFFAAAbaHlJk44MlJ0dm6mLizsqj45XVq15king80jE4u5HBwAAD9ouYlT5tacrLRk5WQyzw0AAIEQbuKYI4MwAwBAqOiWAgAAlkLLDdrFKjf8tMp2AAAIN2gHq9zw0yrbAQA4jXtLJbBArQ0d2RrhrHNp+h+3tLrhp3Q6GMybNDghWj6ssh1IPLQW+ka9oC2hHL9puUlQ5taGjJROeuDaQbqoMFv1jU1ypKfogWU7tH5Px7RGVNa4fAYC6V+3iEiEHZJVtgOJhdZC36gXRAoDihOQs87VKtg8NWmwlm87rLFPrNfKvx3Tz5Zt9wo2UmRvvOnrhp8ZKZ00/Tv99LspQ3W81qW9X9XE/U0+fW1HS9zqApFm/v26JcKNcZ11Lu09VqMtB05E/PedyPWC+EPLTQIytzaY7yI+uCBbT6/e4/O97W2NcDcZm28R4Q5YL36wz+u7I3HW1ZHN1L5uXNoSt7pAOPz9n03U1kJ/rcV1rqZ2/zYTtV6CRXdbdBFuEpC5tcEcZoK98WaoP7aWO7fp3+mnS/p19wQqc8Byc591tTV2xVcZpNM7upqGxlbda5HeobpvXNryVhcZKZ10y6VFGnlWdzlPnW6BYkeEYAXqWknE1sK2Wotf/GCfSt/c7lmvPScziVgvwaK7LfoINwnI3NpgDjOBbrxpT08O+cdm3rm9sGGfnpo0WJL0wZ7jYbUW+ToTfOHmizV/9R6t33M6QG05cMITmDpih+q+censN7Zp3e7KDm2BgvUF6lqZN2lwQrYWBmotdgt0MuNPItZLMIL5P8GJU+Qx5iYBuVsb3MxhZsvBk7qkX3ef7x3VP0ddUjuH3Ldt3rnVuZp01x+3aHBhV/1uylBlpfnPyS1bi/Yeq9EXR6t17+uftdphzlu929NKM7gg22vnGWiHGm6fvPtWF6tmjdZb0y7RSx3wHTgzBNO1Yv79tuTvxrgdOd4lEF+txebfiJt7O0MVbr3Eu2D+TyDyCDcJxL1z+8fxWj04/jzPjsAcZl7YsE8/uqSoVcBx33izpv6bkH9svpqM61xNenr1Hk19aZMyU/2HG3dr0fQ/btEVj6/Vwa/rtN60czTvMM0tUuHuUFseFL78qkZHq+pbHSTcNy5Nstlalcs9UHrKyL764lhNQgyURmwE07Xibi0M5ca4LX87N/zPRl3x2Fr95I9bdPjkqYiWvy2BWovNwulCCqdeEkGw3W2xDK9WRLdUgmhrMN/PrhmohsYmff+i3vr5Wzu0bnelp1XlgWsH6efXDtIpV5PXjTe3HDjh97t87ZgCNRl3SencauyKm7u16P8u+VdLja+dY6DutXB2qC3rrWWXU8uQ5G88BN1UiS+aAzmD7VoJ5sa47nI3GYYe+tPOVqE7mt0a5rFpgbq+w+1CsuINgwP9nwhnmAACI9wkAF99tnWuJpW+ud1rorlgdpZ7j9WEtGNquYP9dv+cNie7y85I9hq70vI1X61FvsrQVveaO4iEukM111sw4wTMO6KOGFsQK2fi1RrRPmj4GqDe8ntbdq34uzFuy3L/bsrQVsHGLVpXEZnHppl/my21twvJajcMDvR/wnzi5xaP+5hE2ocQbhJAsJdIBruzNF/pJLW+QujLr2qU0ilJpUu3e7V6GIahDaZWD3eTsSNDbQYsc2uRr52jeZl50HKoO1RzvQUz6Nm8I+rIy+qjyddB/sqBufp/15+r+sbmhNhZhao9Azlb7sQd6cnqktpZNfXfBKwncwhwc9d1ZY1LX1bW+v0Mc7l9tVi6f6+DC7J1vNYlReGKvpatKrUNjV6txW6J3oXUEdr6PxHKMIF4qM9Ea10i3CSA9l4iGehKJ19dL+YrldxdXbdcWqQ7L+untOROcqS3bh0yB6y2WovMZXAve+Hmi0+Pe/HRvWbufnNra4dqrrdgurXOzs302hF1xNiCcIV61uSvWyMjpZNuHFaoe97Y1mYXXaILd96UULsyzcxdK/b0ZKV0StLsN7cHdWAwlzuac0oFYv59x6ILKZjfQbi/lUiGfHNAnvuDC1RT/03AEz+zeLgEPhGv+CLcJID2XiLZ1pVOt1xapFsuKVLvrun6/5Z/7rXz9tVi4R5A/PTqPVo1a7TOzs30+73+WotalmGaKSw9HWCHGewO1VxvwXZrtTw4NXzTFNR7QtWeOYbc/B3MAnVrBNvdFsv7l7VXKCcFbQVBX/WUkdJJxQXZ+kdlrSqcp+TISGm13S1DQFv3L2urro/Xeg8kNbdYxlNXabS7kIL5HbTntxLM+u0tp3m/Gez+PZa/tUScYJFwkwBC6cf3xd+VTpL0zl2Xtjr4tbfFIlBrkbsM2w6e1E3DCtXDtBPx90MJdodqrrdQurXc3+Gsc4U1yZ+/HVF75xhya+tgFky3RqDutuO1LtW6mkK6f1mkJ1gMVlt1HexBw18QNNdTOK0mgQ4M5rr+3ZShXusEM6dUy26qL47VqFuX1oEr0fn6HZiDZvfMVN2/dEebt54J9FsJtH7L97X1+w71M4PZv8eqS6itwG0WD61LZoSbBBCozzbQDizQTr7W1bp1or1XQwRqLXKkJ6trRkqHNmOb6819kLBJbY4bCvQZwRzc/E1Tb7PZQr7yJdSzpkDdGlLg8NrUbPidkTYaEywGw99OP5iDRqAgaH4eTqtJoBYkc12bQ7i5lbNzp/jppgpHuC0Q5v/Xvrb7dDgN/7cSaH0p8MlJqJ8ZaP8uKSZdQubQ7088TrBIuEkQ7blEMtBOPju99X/M9l4NEai1aNmdIwN2awXL387S1/iHx/7tQp99321p+RnNhqEH/7SzzYPb3B9c4DcUtHXli/sM9IizXl9W1noNYg32rCnYbg3Jd+BpefbvamoO6f5l0ewqCeUS6UCDexu+afIbBM3Pw2k1CXRy0dRseJUhUCtnnelkJJ66qQJpTwuEeZ/ia7uDGXxt+6om5N+WW1utMpv2n9DaL77S0D5ddaIu9DGS/vbve4/V+AxL5n1GoG7jUAbGm7ezI6+O6yiEmwQSbv92oDODjJROrcJPOK0cLUVrKvVgdpa+6i3PHtr3uD9j77Eav5flnqj1P019Wzvflmeg5kGsgc6azPNkBOrWkE7vrC7t193ztzWX4X9uusjrMwLdv6w9V5X52wG3fG6+11gwl0ifnZvpd3CveTvNO3Hzc/N2B9NqEujkos71jdcycytnVlqyunf5Vyunuau0rboP5uAXTcF0K/kav+Rm3qf42u5Ag6+D+W21DEONzYZXt7OvVhlfJzD+tLXva2v/7utEMdQW5FAHxpu309c+xP3+eL06jnBzhgjU8mMOP3WuJr368QE9MrFY9Y3NEW8tikTSb+/OMhyBuhiq6r0PVOYdsK8WE3MAMj/3ddbUctzPN03NuvetHZ4DfaBujbTkTsrOSNa/Dy3QfUu3a93uylbfGWgCxWAmWGzrcmXJ941RzTtg83NzV1go3+luJTQP7g10FZ855JvXD7bVxN/JhbklRvJu5TQP3jd/XjCBueV3Rruryh1eza1kwZaxrbm2fG13oMHXgX5bgcrk6/cfzO+15ee4933Bds/5OlEM9P/O3IIcauueeTv9BW5J2nusJu4uKiDcnEH8tfxEembQ9o4TCkYwffDu74zUDj3gbKOme2yZd8C+dnzmAGR+HujSfXMLRiiDt1teFdbyOwNNoBjodV9/i0A3Rg10IAoUFIP5+5v/zwQbBN1dmc2mA2ywLVb+fl++Bq23LLuvE4FAV/QFc5VXdkZK0N0U4WrZemBuJQvmgNtyoLV5ri1fJwqBBl8H+m0FKtMD1w5q9Z2BPtOt5b4vmBZnfxOoBvp/Z25BDrVl1dd+zlfgjue5bwg38Ij0ZZ0dPZV6MH3wUmTHHgRqkeraxf809b52fIEGsZrPmnplp+vhP3/e5vsDdWu05P6bm+faCDSBYqDXff0t3DdGbSusBHoeKCgG8/c3/58JJQi6uzIfCXMepLZ+X+GeCPi7oi/QVV7hzN8TjJatEd26pHhduRTM+CW3dbsrdbKuUfe/tcNz4DTPtZWTmdrqoN9ybqz/d925rca/BPqt5DvS/JYppVNSq7oO9fcXzNVU/kKdr+80M7cghzoWKdCM9L4G4/vajli24BBu0KE6ch6MYPrg3SI1F0OgA1GePc3vNPW+WgYyUjp5fYevM9KWZ03vzfi2V0tNoPWDmZPIXJfmHXTXjGS/9y8zT7Do628RKKwEeh4oKAbz9w+0ncFcxdcR8yC150TA1//JQFd5dUTLjvks3nzlUqDxS2a1rtYz97aca2v13aO9gqbb0D5dddn/+ZZ6ZKdr77Ear/cH+q0sub3Eb5mcp1yt6jrU31+gq6kChTpf+wwzcwtyqGORgpmRvq2Bzu7tiPXcN4QbJCxzK0q0ZhMOdCAKNE29uWXAfOYdqM/ePAA1Elcy+GqRcu+gQ7l/mft1X1ehBAorgZ6Heom0WXV9o4pyugS1nYH4azVxC3VsWXtOBMz/J9OSvQ9+gVrF2tuy4+ss3vz3NofRQNNN+JqioqWqU40661uZfv9PhjrXla8rR1vqkprcqq67ZoQ2vjDQuL1AoW7VrNEhtyCHOhYpmBnp2ztzfkfz/78LiGPuM9ZR/XMkddyditv67rNzM3VhYVfPWdneYzXacuCE9n51+mzx7NxMFRd0Vd+cLpo3abBWzRqtZXeO1KpZozVv0mBPl4d5O17YsE8/uqRIl/br7vWdnrOmdO+dpXv9S9paP4Szf3cZ2voM83b76uY6OzdT3bu0/s62wkqwz31tpzsoFuV08fmdLWWlJQe9ncGK9Oe1R8u/TQ9HmleZQp2/J1AXn7POO7z6ao0w/73dB8zBhV31uylD1btrur5tqje3YIKG+/fs7/9kqL+t3KzUVn/Lluu4g0rL7+yT0yWk/wPhzDvWUnV9Y8D/d+4W5La2e3BBdqtxf+a/tTtQ/cfzH8mRntyqbqN1RWy4aLlBQjPPQROon7gjhHs5ekuhzMdjbi3wdyuLUA6ukRwj5evMMtC4HfOVSebn5q6wU66msAbmRnosWEePLQuHuasq1Pl7Qu3i9XUW76uVxNxK5qtbyd8UFS3XCfb3HOpcV+GMfwrl/0A484615Os2Me7tcncjbjlwotX9rFput7llNZxW72hcEdseNsMwjJiWIMqqqqrkcDjkdDplt4c42Qni3uGTp9rcMZlv8RAJbd03yP29HTl7aDS3MxzmMnqulnp/j8/Zm0/97y0b3Dto8w47mOCQCPUSTe4Bvs2GoQeXf+6pd/OVav9z00W685VPPe8zPzdbdudIXVjY1fN877EaXfH4Wq91/HVttfx7uMvo6+8bq79nyzLZQ5j8Llj+tisjpZN+8sctbYaGtvYpoVy5ZP57/W7KUE19aVOb5W1r3F60/z6hHL8JN7AcfzvLSPO1U28pmMG84YrmdobLVxkldWi5E6FeYqHlgcgdPBZ+sE8b/ncyu5YHt1APds46l88Dsju8Du3TVbUN34T194j137OjLneOZKgL9STL/Pcyh91A729rOzoiCLZEuPGDcINI2nLghG74n41tvm4+wwViqa0DUaCWnZbaOthZsdUsVi2z7u8ONtSFc5LlL+y6hfr36+h5b0I5fjPmBmiHeB9UB7Tk71YkLce/hHP7lXgce9Re4dxYM1JCuXounCuXInHfvZbibd4bwg3QDvE+qA4IViQOdh05r1UsxPvlzm7hnmT5C7uhimUQ9IVwA7RDNG4zAURLJA92VpAoLbPxcJIVb0EwLua5mT9/vvr27au0tDQNHz5cH3/8cZvrLly4UDabzeuRlpYWxdIC3txnvG3NYwMgMblDgy/x1DIbD/MtxVsQjHnLzauvvqpZs2ZpwYIFGj58uJ544gmNHTtWu3btUm5urs/32O127dq1y/PcZrNFq7iAT1ZrjgeQWC2zsR7zFA+tRy3F/Gqp4cOH6+KLL9bTTz8tSWpublZBQYF+8pOfaPbs2a3WX7hwoWbOnKmTJ08G9fkNDQ1qaGjwPK+qqlJBQQFXSwEAghLry9ETRUdfMZcwV0u5XC5t3rxZpaWlnmVJSUkaM2aMysvL23xfTU2N+vTpo+bmZl100UX61a9+pXPPPdfnumVlZfrlL38Z8bIDAM4MtMwGJ9atRy3FdMxNZWWlmpqalJeX57U8Ly9PFRUVPt8zYMAAvfDCC3rrrbf0hz/8Qc3NzRo5cqT++c9/+ly/tLRUTqfT8zh48GDEtwMAAAS+/1y0xHzMTahKSkpUUvKv29KPHDlSAwcO1LPPPquHHnqo1fqpqalKTU2NZhEBAEAMxbTlJicnR506ddLRo0e9lh89elT5+flBfUZycrIGDx6sPXt83+QNAACcWWIablJSUjRkyBCtWrXKs6y5uVmrVq3yap3xp6mpSdu3b1ePHj06qpgAACCBxLxbatasWZoyZYqGDh2qYcOG6YknnlBtba1+9KMfSZImT56sXr16qaysTJL04IMPasSIEerXr59OnjypuXPnav/+/br11ltjuRkAACBOxDzc3Hjjjfrqq6/085//XBUVFbrwwgv13nvveQYZHzhwQElJ/2pgOnHihG677TZVVFSoa9euGjJkiDZu3KhBgwbFahMAAEAcifk8N9HGXcEBAEg8oRy/4+L2CwAAAJFCuAEAAJZCuAEAAJZCuAEAAJYS86ulos09frqqqirGJQEAAMFyH7eDuQ7qjAs31dXVkqSCgoIYlwQAAISqurpaDofD7zpn3KXgzc3NOnz4sLKysmSz2SL62VVVVSooKNDBgwe5zLydqMvIoB4jh7qMHOoycs6kujQMQ9XV1erZs6fX/He+nHEtN0lJSerdu3eHfofdbrf8f7JooS4jg3qMHOoycqjLyDlT6jJQi40bA4oBAIClEG4AAIClEG4iKDU1Vb/4xS+Umpoa66IkPOoyMqjHyKEuI4e6jBzq0rczbkAxAACwNlpuAACApRBuAACApRBuAACApRBuAACApRBuImT+/Pnq27ev0tLSNHz4cH388cexLlLcKysr08UXX6ysrCzl5uZqwoQJ2rVrl9c69fX1mjZtmrp3767MzExNnDhRR48ejVGJE8OcOXNks9k0c+ZMzzLqMXiHDh3Sf/7nf6p79+5KT0/X+eefr02bNnleNwxDP//5z9WjRw+lp6drzJgx2r17dwxLHJ+ampr0wAMPqKioSOnp6Tr77LP10EMPed0XiLr0bd26dbruuuvUs2dP2Ww2LVu2zOv1YOrt66+/1k033SS73a7s7GxNnTpVNTU1UdyKGDPQbosXLzZSUlKMF154wdi5c6dx2223GdnZ2cbRo0djXbS4NnbsWOPFF180duzYYWzdutW4+uqrjcLCQqOmpsazzu23324UFBQYq1atMjZt2mSMGDHCGDlyZAxLHd8+/vhjo2/fvkZxcbExY8YMz3LqMThff/210adPH+Pmm282PvroI+PLL780VqxYYezZs8ezzpw5cwyHw2EsW7bM+Oyzz4zrr7/eKCoqMk6dOhXDksefhx9+2OjevbuxfPlyY9++fcaSJUuMzMxM48knn/SsQ1369s477xg/+9nPjDfffNOQZCxdutTr9WDq7aqrrjIuuOAC48MPPzTWr19v9OvXz5g0aVKUtyR2CDcRMGzYMGPatGme501NTUbPnj2NsrKyGJYq8Rw7dsyQZKxdu9YwDMM4efKkkZycbCxZssSzzt/+9jdDklFeXh6rYsat6upqo3///sbKlSuN0aNHe8IN9Ri8e++917j00kvbfL25udnIz8835s6d61l28uRJIzU11fjjH/8YjSImjGuuuca45ZZbvJZ973vfM2666SbDMKjLYJnDTTD19vnnnxuSjE8++cSzzrvvvmvYbDbj0KFDUSt7LNEt1U4ul0ubN2/WmDFjPMuSkpI0ZswYlZeXx7BkicfpdEqSunXrJknavHmzGhsbver2nHPOUWFhIXXrw7Rp03TNNdd41ZdEPYbi7bff1tChQ/WDH/xAubm5Gjx4sH772996Xt+3b58qKiq86tLhcGj48OHUpcnIkSO1atUqffHFF5Kkzz77TBs2bNC4ceMkUZfhCqbeysvLlZ2draFDh3rWGTNmjJKSkvTRRx9FvcyxcMbdODPSKisr1dTUpLy8PK/leXl5+vvf/x6jUiWe5uZmzZw5U5dcconOO+88SVJFRYVSUlKUnZ3ttW5eXp4qKipiUMr4tXjxYn366af65JNPWr1GPQbvyy+/1DPPPKNZs2bpvvvu0yeffKK77rpLKSkpmjJliqe+fP3eqUtvs2fPVlVVlc455xx16tRJTU1Nevjhh3XTTTdJEnUZpmDqraKiQrm5uV6vd+7cWd26dTtj6pZwg7gwbdo07dixQxs2bIh1URLOwYMHNWPGDK1cuVJpaWmxLk5Ca25u1tChQ/WrX/1KkjR48GDt2LFDCxYs0JQpU2JcusTy2muv6ZVXXtGiRYt07rnnauvWrZo5c6Z69uxJXaLD0S3VTjk5OerUqVOrK0+OHj2q/Pz8GJUqsUyfPl3Lly/X+++/r969e3uW5+fny+Vy6eTJk17rU7feNm/erGPHjumiiy5S586d1blzZ61du1ZPPfWUOnfurLy8POoxSD169NCgQYO8lg0cOFAHDhyQJE998XsP7L//+781e/Zs/fu//7vOP/98/fCHP9RPf/pTlZWVSaIuwxVMveXn5+vYsWNer3/zzTf6+uuvz5i6Jdy0U0pKioYMGaJVq1Z5ljU3N2vVqlUqKSmJYcnin2EYmj59upYuXarVq1erqKjI6/UhQ4YoOTnZq2537dqlAwcOULctXHHFFdq+fbu2bt3qeQwdOlQ33XST59/UY3AuueSSVtMRfPHFF+rTp48kqaioSPn5+V51WVVVpY8++oi6NKmrq1NSkvchplOnTmpubpZEXYYrmHorKSnRyZMntXnzZs86q1evVnNzs4YPHx71MsdErEc0W8HixYuN1NRUY+HChcbnn39u/Nd//ZeRnZ1tVFRUxLpoce2OO+4wHA6HsWbNGuPIkSOeR11dnWed22+/3SgsLDRWr15tbNq0ySgpKTFKSkpiWOrE0PJqKcOgHoP18ccfG507dzYefvhhY/fu3cYrr7xiZGRkGH/4wx8868yZM8fIzs423nrrLWPbtm3G+PHjuXzZhylTphi9evXyXAr+5ptvGjk5OcY999zjWYe69K26utrYsmWLsWXLFkOS8fjjjxtbtmwx9u/fbxhGcPV21VVXGYMHDzY++ugjY8OGDUb//v25FByhmzdvnlFYWGikpKQYw4YNMz788MNYFynuSfL5ePHFFz3rnDp1yrjzzjuNrl27GhkZGcYNN9xgHDlyJHaFThDmcEM9Bu9Pf/qTcd555xmpqanGOeecYzz33HNerzc3NxsPPPCAkZeXZ6SmphpXXHGFsWvXrhiVNn5VVVUZM2bMMAoLC420tDTjrLPOMn72s58ZDQ0NnnWoS9/ef/99n/vGKVOmGIYRXL0dP37cmDRpkpGZmWnY7XbjRz/6kVFdXR2DrYkNm2G0mC4SAAAgwTHmBgAAWArhBgAAWArhBgAAWArhBgAAWArhBgAAWArhBgAAWArhBgAAWArhBgAAWArhBgD+l81m07Jly2JdDADtRLgBEFVfffWV7rjjDhUWFio1NVX5+fkaO3asPvjgg1gXDYBFdI51AQCcWSZOnCiXy6WXXnpJZ511lo4ePapVq1bp+PHjsS4aAIug5QZA1Jw8eVLr16/XI488ossvv1x9+vTRsGHDVFpaquuvv16S9Pjjj+v8889Xly5dVFBQoDvvvFM1NTWez1i4cKGys7O1fPlyDRgwQBkZGfr+97+vuro6vfTSS+rbt6+6du2qu+66S01NTZ739e3bVw899JAmTZqkLl26qFevXpo/f77f8h48eFD/9m//puzsbHXr1k3jx4/XP/7xjw6pGwCRQ7gBEDWZmZnKzMzUsmXL1NDQ4HOdpKQkPfXUU9q5c6deeuklrV69Wvfcc4/XOnV1dXrqqae0ePFivffee1qzZo1uuOEGvfPOO3rnnXf08ssv69lnn9Xrr7/u9b65c+fqggsu0JYtWzR79mzNmDFDK1eu9FmOxsZGjR07VllZWVq/fr0++OADZWZm6qqrrpLL5YpMhQDoGLG+LTmAM8vrr79udO3a1UhLSzNGjhxplJaWGp999lmb6y9ZssTo3r275/mLL75oSDL27NnjWfbjH//YyMjIMKqrqz3Lxo4da/z4xz/2PO/Tp49x1VVXeX32jTfeaIwbN87zXJKxdOlSwzAM4+WXXzYGDBhgNDc3e15vaGgw0tPTjRUrVoS+4QCihpYbAFE1ceJEHT58WG+//bauuuoqrVmzRhdddJEWLlwoSfrrX/+qK664Qr169VJWVpZ++MMf6vjx46qrq/N8RkZGhs4++2zP87y8PPXt21eZmZley44dO+b13SUlJa2e/+1vf/NZzs8++0x79uxRVlaWp8WpW7duqq+v1969e9tbDQA6EAOKAURdWlqarrzySl155ZV64IEHdOutt+oXv/iFLrvsMl177bW644479PDDD6tbt27asGGDpk6dKpfLpYyMDElScnKy1+fZbDafy5qbm8MuY01NjYYMGaJXXnml1Wvf+ta3wv5cAB2PcAMg5gYNGqRly5Zp8+bNam5u1mOPPaakpNMNy6+99lrEvufDDz9s9XzgwIE+173ooov06quvKjc3V3a7PWJlANDx6JYCEDXHjx/Xd77zHf3hD3/Qtm3btG/fPi1ZskSPPvqoxo8fr379+qmxsVHz5s3Tl19+qZdfflkLFiyI2Pd/8MEHevTRR/XFF19o/vz5WrJkiWbMmOFz3Ztuukk5OTkaP3681q9fr3379mnNmjW666679M9//jNiZQIQebTcAIiazMxMDR8+XL/5zW+0d+9eNTY2qqCgQLfddpvuu+8+paen6/HHH9cjjzyi0tJSjRo1SmVlZZo8eXJEvv/uu+/Wpk2b9Mtf/lJ2u12PP/64xo4d63PdjIwMrVu3Tvfee6++973vqbq6Wr169dIVV1xBSw4Q52yGYRixLgQAdLS+fftq5syZmjlzZqyLAqCD0S0FAAAshXADAAAshW4pAABgKbTcAAAASyHcAAAASyHcAAAASyHcAAAASyHcAAAASyHcAAAASyHcAAAASyHcAAAAS/n/AWfQ4zzjcOCeAAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "sigma_df_mcmc = pd.DataFrame(np.concatenate((np.expand_dims(np.arange(bart_model.num_samples),axis=1), np.expand_dims(bart_model.global_var_samples,axis=1)), axis = 1), columns=[\"Sample\", \"Sigma\"])\n", + "sns.scatterplot(data=sigma_df_mcmc, x=\"Sample\", y=\"Sigma\")\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Compute the test set RMSE" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "1.2873862776376062" + ] + }, + "execution_count": 7, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "np.sqrt(np.mean(np.power(y_test - np.squeeze(y_avg_mcmc),2)))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Check the variable split count in the last \"GFR\" sample" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([28, 23, 23, 19, 18, 28, 35, 34, 40, 29], dtype=int32)" + ] + }, + "execution_count": 8, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "bart_model.forest_container_mean.get_forest_split_counts(9, p_X)" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([3199, 2738, 2626, 2031, 1976, 3262, 2024, 2764, 2891, 3670],\n", + " dtype=int32)" + ] + }, + "execution_count": 9, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "bart_model.forest_container_mean.get_overall_split_counts(p_X)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The split counts appear relatively uniform across features, so let's dig deeper and look at individual trees, starting with the first tree in the last \"grow-from-root\" sample." + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [], + "source": [ + "splits = bart_model.forest_container_mean.get_granular_split_counts(p_X)" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([0, 0, 0, 0, 0, 0, 0, 0, 0, 1], dtype=int32)" + ] + }, + "execution_count": 11, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "splits[9,0,:]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This tree has a single split on the only \"important\" feature. Now, let's look at the second tree." + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([0, 0, 0, 0, 0, 0, 0, 0, 0, 1], dtype=int32)" + ] + }, + "execution_count": 12, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "splits[9,1,:]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This tree also only splits on the important feature." + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([0, 1, 0, 0, 0, 0, 0, 0, 0, 0], dtype=int32)" + ] + }, + "execution_count": 13, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "splits[9,20,:]" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "array([0, 0, 0, 0, 0, 0, 0, 0, 1, 1], dtype=int32)" + ] + }, + "execution_count": 14, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "splits[9,30,:]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We see that \"later\" trees are splitting on other features, but we also note that these trees are fitting an outcome that is already residualized many \"relevant splits\" made by trees 1 and 2." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now, let's inspect the first tree for this last GFR sample in more depth, following [this scikit-learn vignette](https://scikit-learn.org/stable/auto_examples/tree/plot_unveil_tree_structure.html)." + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": {}, + "outputs": [], + "source": [ + "forest_num = 9\n", + "tree_num = 0" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "node=0 is a split node, which tells us to go to node 1 if X[:, 9] <= 0.49971201341971494 else to node 2.\n", + "\tnode=1 is a leaf node with value=[-0.313].\n", + "\tnode=2 is a leaf node with value=[0.406].\n" + ] + } + ], + "source": [ + "nodes = np.sort(bart_model.forest_container_mean.nodes(forest_num,tree_num))\n", + "for nid in nodes:\n", + " if bart_model.forest_container_mean.is_leaf_node(forest_num,tree_num,nid):\n", + " print(\n", + " \"{space}node={node} is a leaf node with value={value}.\".format(\n", + " space=bart_model.forest_container_mean.node_depth(forest_num,tree_num,nid) * \"\\t\", \n", + " node=nid, value=np.around(bart_model.forest_container_mean.node_leaf_values(forest_num,tree_num,nid), 3)\n", + " )\n", + " )\n", + " else:\n", + " print(\n", + " \"{space}node={node} is a split node, which tells us to \"\n", + " \"go to node {left} if X[:, {feature}] <= {threshold} \"\n", + " \"else to node {right}.\".format(\n", + " space=bart_model.forest_container_mean.node_depth(forest_num,tree_num,nid) * \"\\t\",\n", + " node=nid,\n", + " left=bart_model.forest_container_mean.left_child_node(forest_num,tree_num,nid),\n", + " feature=bart_model.forest_container_mean.node_split_index(forest_num,tree_num,nid),\n", + " threshold=bart_model.forest_container_mean.node_split_threshold(forest_num,tree_num,nid),\n", + " right=bart_model.forest_container_mean.right_child_node(forest_num,tree_num,nid),\n", + " )\n", + " )" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "venv", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.8.17" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/include/stochtree/tree.h b/include/stochtree/tree.h index d027c4e6..b252f9cb 100644 --- a/include/stochtree/tree.h +++ b/include/stochtree/tree.h @@ -461,6 +461,22 @@ class Tree { return node_type_[nid]; } + /*! + * \brief Whether the node is a numeric split node + * \param nid ID of node being queried + */ + bool IsNumericSplitNode(std::int32_t nid) const { + return node_type_[nid] == TreeNodeType::kNumericalSplitNode; + } + + /*! + * \brief Whether the node is a numeric split node + * \param nid ID of node being queried + */ + bool IsCategoricalSplitNode(std::int32_t nid) const { + return node_type_[nid] == TreeNodeType::kCategoricalSplitNode; + } + /*! * \brief Query whether this tree contains any categorical splits */ @@ -500,18 +516,35 @@ class Tree { [[nodiscard]] std::vector const& GetInternalNodes() const { return internal_nodes_; } + /*! * \brief Get indices of all leaf nodes. */ [[nodiscard]] std::vector const& GetLeaves() const { return leaves_; } + /*! * \brief Get indices of all leaf parent nodes. */ [[nodiscard]] std::vector const& GetLeafParents() const { return leaf_parents_; } + + /*! + * \brief Get indices of all valid (non-deleted) nodes. + */ + [[nodiscard]] std::vector GetNodes() { + std::vector output; + auto const& self = *this; + this->WalkTree([&output, &self](std::int32_t nidx) { + if (!self.IsDeleted(nidx)) { + output.push_back(nidx); + } + return true; + }); + return output; + } /*! * \brief Get the depth of a node diff --git a/src/forest.cpp b/src/forest.cpp index b0f47c96..d510a16f 100644 --- a/src/forest.cpp +++ b/src/forest.cpp @@ -227,7 +227,8 @@ cpp11::writable::integers get_tree_split_counts_forest_container_cpp(cpp11::exte StochTree::Tree* tree = ensemble->GetTree(tree_num); std::vector split_nodes = tree->GetInternalNodes(); for (int i = 0; i < split_nodes.size(); i++) { - auto split_feature = split_nodes.at(i); + auto node_id = split_nodes.at(i); + auto split_feature = tree->SplitIndex(node_id); output.at(split_feature)++; } return output; @@ -243,7 +244,8 @@ cpp11::writable::integers get_forest_split_counts_forest_container_cpp(cpp11::ex StochTree::Tree* tree = ensemble->GetTree(i); std::vector split_nodes = tree->GetInternalNodes(); for (int j = 0; j < split_nodes.size(); j++) { - auto split_feature = split_nodes.at(j); + auto node_id = split_nodes.at(j); + auto split_feature = tree->SplitIndex(node_id); output.at(split_feature)++; } } @@ -262,7 +264,8 @@ cpp11::writable::integers get_overall_split_counts_forest_container_cpp(cpp11::e StochTree::Tree* tree = ensemble->GetTree(j); std::vector split_nodes = tree->GetInternalNodes(); for (int k = 0; k < split_nodes.size(); k++) { - auto split_feature = split_nodes.at(k); + auto node_id = split_nodes.at(k); + auto split_feature = tree->SplitIndex(node_id); output.at(split_feature)++; } } @@ -282,8 +285,9 @@ cpp11::writable::integers get_granular_split_count_array_forest_container_cpp(cp StochTree::Tree* tree = ensemble->GetTree(j); std::vector split_nodes = tree->GetInternalNodes(); for (int k = 0; k < split_nodes.size(); k++) { - auto split_feature = split_nodes.at(k); - output.at(num_features*num_trees*i + split_feature*num_trees + j)++; + auto node_id = split_nodes.at(k); + auto split_feature = tree->SplitIndex(node_id); + output.at(split_feature*num_samples*num_trees + j*num_samples + i)++; } } } diff --git a/src/py_stochtree.cpp b/src/py_stochtree.cpp index 0420e396..45ff4127 100644 --- a/src/py_stochtree.cpp +++ b/src/py_stochtree.cpp @@ -164,7 +164,7 @@ class ForestContainerCpp { return forest_samples_->NumSamples(); } - int NumLeaves(int forest_num) { + int NumLeavesForest(int forest_num) { StochTree::TreeEnsemble* forest = forest_samples_->GetEnsemble(forest_num); return forest->NumLeaves(); } @@ -428,7 +428,8 @@ class ForestContainerCpp { StochTree::Tree* tree = ensemble->GetTree(i); std::vector split_nodes = tree->GetInternalNodes(); for (int j = 0; j < split_nodes.size(); j++) { - auto split_feature = split_nodes.at(j); + auto node_id = split_nodes.at(j); + auto split_feature = tree->SplitIndex(node_id); accessor(split_feature)++; } } @@ -449,7 +450,8 @@ class ForestContainerCpp { StochTree::Tree* tree = ensemble->GetTree(j); std::vector split_nodes = tree->GetInternalNodes(); for (int k = 0; k < split_nodes.size(); k++) { - auto split_feature = split_nodes.at(k); + auto node_id = split_nodes.at(k); + auto split_feature = tree->SplitIndex(node_id); accessor(split_feature)++; } } @@ -460,11 +462,11 @@ class ForestContainerCpp { py::array_t GetGranularSplitCounts(int num_features) { int num_samples = forest_samples_->NumSamples(); int num_trees = forest_samples_->NumTrees(); - auto result = py::array_t(py::detail::any_container({num_trees,num_features,num_samples})); + auto result = py::array_t(py::detail::any_container({num_samples,num_trees,num_features})); auto accessor = result.mutable_unchecked<3>(); - for (int i = 0; i < num_trees; i++) { - for (int j = 0; j < num_features; j++) { - for (int k = 0; k < num_samples; k++) { + for (int i = 0; i < num_samples; i++) { + for (int j = 0; j < num_trees; j++) { + for (int k = 0; k < num_features; k++) { accessor(i,j,k) = 0; } } @@ -475,14 +477,144 @@ class ForestContainerCpp { StochTree::Tree* tree = ensemble->GetTree(j); std::vector split_nodes = tree->GetInternalNodes(); for (int k = 0; k < split_nodes.size(); k++) { - auto split_feature = split_nodes.at(k); - accessor(j,split_feature,i)++; + auto node_id = split_nodes.at(k); + auto split_feature = tree->SplitIndex(node_id); + accessor(i,j,split_feature)++; } } } return result; } + bool IsLeafNode(int forest_id, int tree_id, int node_id) { + StochTree::TreeEnsemble* ensemble = forest_samples_->GetEnsemble(forest_id); + StochTree::Tree* tree = ensemble->GetTree(tree_id); + return tree->IsLeaf(node_id); + } + + bool IsNumericSplitNode(int forest_id, int tree_id, int node_id) { + StochTree::TreeEnsemble* ensemble = forest_samples_->GetEnsemble(forest_id); + StochTree::Tree* tree = ensemble->GetTree(tree_id); + return tree->IsNumericSplitNode(node_id); + } + + bool IsCategoricalSplitNode(int forest_id, int tree_id, int node_id) { + StochTree::TreeEnsemble* ensemble = forest_samples_->GetEnsemble(forest_id); + StochTree::Tree* tree = ensemble->GetTree(tree_id); + return tree->IsCategoricalSplitNode(node_id); + } + + int ParentNode(int forest_id, int tree_id, int node_id) { + StochTree::TreeEnsemble* ensemble = forest_samples_->GetEnsemble(forest_id); + StochTree::Tree* tree = ensemble->GetTree(tree_id); + return tree->Parent(node_id); + } + + int LeftChildNode(int forest_id, int tree_id, int node_id) { + StochTree::TreeEnsemble* ensemble = forest_samples_->GetEnsemble(forest_id); + StochTree::Tree* tree = ensemble->GetTree(tree_id); + return tree->LeftChild(node_id); + } + + int RightChildNode(int forest_id, int tree_id, int node_id) { + StochTree::TreeEnsemble* ensemble = forest_samples_->GetEnsemble(forest_id); + StochTree::Tree* tree = ensemble->GetTree(tree_id); + return tree->RightChild(node_id); + } + + int SplitIndex(int forest_id, int tree_id, int node_id) { + StochTree::TreeEnsemble* ensemble = forest_samples_->GetEnsemble(forest_id); + StochTree::Tree* tree = ensemble->GetTree(tree_id); + return tree->SplitIndex(node_id); + } + + int NodeDepth(int forest_id, int tree_id, int node_id) { + StochTree::TreeEnsemble* ensemble = forest_samples_->GetEnsemble(forest_id); + StochTree::Tree* tree = ensemble->GetTree(tree_id); + return tree->GetDepth(node_id); + } + + double SplitThreshold(int forest_id, int tree_id, int node_id) { + StochTree::TreeEnsemble* ensemble = forest_samples_->GetEnsemble(forest_id); + StochTree::Tree* tree = ensemble->GetTree(tree_id); + return tree->Threshold(node_id); + } + + py::array_t SplitCategories(int forest_id, int tree_id, int node_id) { + StochTree::TreeEnsemble* ensemble = forest_samples_->GetEnsemble(forest_id); + StochTree::Tree* tree = ensemble->GetTree(tree_id); + std::vector raw_categories = tree->CategoryList(node_id); + int num_categories = raw_categories.size(); + auto result = py::array_t(py::detail::any_container({num_categories})); + auto accessor = result.mutable_unchecked<1>(); + for (int i = 0; i < num_categories; i++) { + accessor(i) = raw_categories.at(i); + } + return result; + } + + py::array_t NodeLeafValues(int forest_id, int tree_id, int node_id) { + StochTree::TreeEnsemble* ensemble = forest_samples_->GetEnsemble(forest_id); + StochTree::Tree* tree = ensemble->GetTree(tree_id); + int num_outputs = tree->OutputDimension(); + auto result = py::array_t(py::detail::any_container({num_outputs})); + auto accessor = result.mutable_unchecked<1>(); + for (int i = 0; i < num_outputs; i++) { + accessor(i) = tree->LeafValue(node_id, i); + } + return result; + } + + int NumNodes(int forest_id, int tree_id) { + StochTree::TreeEnsemble* ensemble = forest_samples_->GetEnsemble(forest_id); + StochTree::Tree* tree = ensemble->GetTree(tree_id); + return tree->NumValidNodes(); + } + + int NumLeaves(int forest_id, int tree_id) { + StochTree::TreeEnsemble* ensemble = forest_samples_->GetEnsemble(forest_id); + StochTree::Tree* tree = ensemble->GetTree(tree_id); + return tree->NumLeaves(); + } + + int NumLeafParents(int forest_id, int tree_id) { + StochTree::TreeEnsemble* ensemble = forest_samples_->GetEnsemble(forest_id); + StochTree::Tree* tree = ensemble->GetTree(tree_id); + return tree->NumLeafParents(); + } + + int NumSplitNodes(int forest_id, int tree_id) { + StochTree::TreeEnsemble* ensemble = forest_samples_->GetEnsemble(forest_id); + StochTree::Tree* tree = ensemble->GetTree(tree_id); + return tree->NumSplitNodes(); + } + + py::array_t Nodes(int forest_id, int tree_id) { + StochTree::TreeEnsemble* ensemble = forest_samples_->GetEnsemble(forest_id); + StochTree::Tree* tree = ensemble->GetTree(tree_id); + std::vector nodes = tree->GetNodes(); + int num_nodes = nodes.size(); + auto result = py::array_t(py::detail::any_container({num_nodes})); + auto accessor = result.mutable_unchecked<1>(); + for (int i = 0; i < num_nodes; i++) { + accessor(i) = nodes.at(i); + } + return result; + } + + py::array_t Leaves(int forest_id, int tree_id) { + StochTree::TreeEnsemble* ensemble = forest_samples_->GetEnsemble(forest_id); + StochTree::Tree* tree = ensemble->GetTree(tree_id); + std::vector leaves = tree->GetLeaves(); + int num_leaves = leaves.size(); + auto result = py::array_t(py::detail::any_container({num_leaves})); + auto accessor = result.mutable_unchecked<1>(); + for (int i = 0; i < num_leaves; i++) { + accessor(i) = leaves.at(i); + } + return result; + } + private: std::unique_ptr forest_samples_; int num_trees_; @@ -1044,8 +1176,25 @@ PYBIND11_MODULE(stochtree_cpp, m) { .def("GetForestSplitCounts", &ForestContainerCpp::GetForestSplitCounts) .def("GetOverallSplitCounts", &ForestContainerCpp::GetOverallSplitCounts) .def("GetGranularSplitCounts", &ForestContainerCpp::GetGranularSplitCounts) + .def("NumLeavesForest", &ForestContainerCpp::NumLeavesForest) + .def("SumLeafSquared", &ForestContainerCpp::SumLeafSquared) + .def("IsLeafNode", &ForestContainerCpp::IsLeafNode) + .def("IsNumericSplitNode", &ForestContainerCpp::IsNumericSplitNode) + .def("IsCategoricalSplitNode", &ForestContainerCpp::IsCategoricalSplitNode) + .def("ParentNode", &ForestContainerCpp::ParentNode) + .def("LeftChildNode", &ForestContainerCpp::LeftChildNode) + .def("RightChildNode", &ForestContainerCpp::RightChildNode) + .def("SplitIndex", &ForestContainerCpp::SplitIndex) + .def("NodeDepth", &ForestContainerCpp::NodeDepth) + .def("SplitThreshold", &ForestContainerCpp::SplitThreshold) + .def("SplitCategories", &ForestContainerCpp::SplitCategories) + .def("NodeLeafValues", &ForestContainerCpp::NodeLeafValues) + .def("NumNodes", &ForestContainerCpp::NumNodes) .def("NumLeaves", &ForestContainerCpp::NumLeaves) - .def("SumLeafSquared", &ForestContainerCpp::SumLeafSquared); + .def("NumLeafParents", &ForestContainerCpp::NumLeafParents) + .def("NumSplitNodes", &ForestContainerCpp::NumSplitNodes) + .def("Nodes", &ForestContainerCpp::Nodes) + .def("Leaves", &ForestContainerCpp::Leaves); py::class_(m, "ForestSamplerCpp") .def(py::init, int, data_size_t, double, double, int, int>()) diff --git a/stochtree/forest.py b/stochtree/forest.py index 866208ee..4eefae79 100644 --- a/stochtree/forest.py +++ b/stochtree/forest.py @@ -11,6 +11,10 @@ class ForestContainer: def __init__(self, num_trees: int, output_dimension: int, leaf_constant: bool, is_exponentiated: bool) -> None: # Initialize a ForestContainerCpp object self.forest_container_cpp = ForestContainerCpp(num_trees, output_dimension, leaf_constant, is_exponentiated) + self.num_trees = num_trees + self.output_dimension = output_dimension + self.leaf_constant = leaf_constant + self.is_exponentiated = is_exponentiated def predict(self, dataset: Dataset) -> np.array: # Predict samples from Dataset @@ -149,14 +153,14 @@ def get_granular_split_counts(self, num_features: int) -> np.array: """ return self.forest_container_cpp.GetGranularSplitCounts(num_features) - def num_leaves(self, forest_num: int) -> int: + def num_forest_leaves(self, forest_num: int) -> int: """ Return the total number of leaves for a given forest in the ``ForestContainer`` forest_num : :obj:`int` Index of the forest to be queried """ - return self.forest_container_cpp.NumLeaves(forest_num) + return self.forest_container_cpp.NumLeavesForest(forest_num) def sum_leaves_squared(self, forest_num: int) -> float: """ @@ -166,4 +170,227 @@ def sum_leaves_squared(self, forest_num: int) -> float: Index of the forest to be queried """ return self.forest_container_cpp.SumLeafSquared(forest_num) + + def is_leaf_node(self, forest_num: int, tree_num: int, node_id: int) -> bool: + """ + Whether or not a given node of a given tree in a given forest in the ``ForestContainer`` is a leaf + + forest_num : :obj:`int` + Index of the forest to be queried + tree_num : :obj:`int` + Index of the tree to be queried + node_id : :obj:`int` + Index of the node to be queried + """ + return self.forest_container_cpp.IsLeafNode(forest_num, tree_num, node_id) + + def is_numeric_split_node(self, forest_num: int, tree_num: int, node_id: int) -> bool: + """ + Whether or not a given node of a given tree in a given forest in the ``ForestContainer`` is a numeric split node + + forest_num : :obj:`int` + Index of the forest to be queried + tree_num : :obj:`int` + Index of the tree to be queried + node_id : :obj:`int` + Index of the node to be queried + """ + return self.forest_container_cpp.IsNumericSplitNode(forest_num, tree_num, node_id) + + def is_categorical_split_node(self, forest_num: int, tree_num: int, node_id: int) -> bool: + """ + Whether or not a given node of a given tree in a given forest in the ``ForestContainer`` is a categorical split node + + forest_num : :obj:`int` + Index of the forest to be queried + tree_num : :obj:`int` + Index of the tree to be queried + node_id : :obj:`int` + Index of the node to be queried + """ + return self.forest_container_cpp.IsCategoricalSplitNode(forest_num, tree_num, node_id) + + def parent_node(self, forest_num: int, tree_num: int, node_id: int) -> int: + """ + Parent node of given node of a given tree in a given forest in the ``ForestContainer`` + + forest_num : :obj:`int` + Index of the forest to be queried + tree_num : :obj:`int` + Index of the tree to be queried + node_id : :obj:`int` + Index of the node to be queried + """ + return self.forest_container_cpp.ParentNode(forest_num, tree_num, node_id) + + def left_child_node(self, forest_num: int, tree_num: int, node_id: int) -> int: + """ + Left child node of given node of a given tree in a given forest in the ``ForestContainer`` + + forest_num : :obj:`int` + Index of the forest to be queried + tree_num : :obj:`int` + Index of the tree to be queried + node_id : :obj:`int` + Index of the node to be queried + """ + return self.forest_container_cpp.LeftChildNode(forest_num, tree_num, node_id) + + def right_child_node(self, forest_num: int, tree_num: int, node_id: int) -> int: + """ + Right child node of given node of a given tree in a given forest in the ``ForestContainer`` + + forest_num : :obj:`int` + Index of the forest to be queried + tree_num : :obj:`int` + Index of the tree to be queried + node_id : :obj:`int` + Index of the node to be queried + """ + return self.forest_container_cpp.RightChildNode(forest_num, tree_num, node_id) + + def node_depth(self, forest_num: int, tree_num: int, node_id: int) -> int: + """ + Depth of given node of a given tree in a given forest in the ``ForestContainer``. + Returns ``-1`` if the node is a leaf. + + forest_num : :obj:`int` + Index of the forest to be queried + tree_num : :obj:`int` + Index of the tree to be queried + node_id : :obj:`int` + Index of the node to be queried + """ + return self.forest_container_cpp.NodeDepth(forest_num, tree_num, node_id) + + def node_split_index(self, forest_num: int, tree_num: int, node_id: int) -> int: + """ + Split index of given node of a given tree in a given forest in the ``ForestContainer``. + Returns ``-1`` if the node is a leaf. + + forest_num : :obj:`int` + Index of the forest to be queried + tree_num : :obj:`int` + Index of the tree to be queried + node_id : :obj:`int` + Index of the node to be queried + """ + if self.is_leaf_node(forest_num, tree_num, node_id): + return -1 + else: + return self.forest_container_cpp.SplitIndex(forest_num, tree_num, node_id) + + def node_split_threshold(self, forest_num: int, tree_num: int, node_id: int) -> float: + """ + Threshold that defines a numeric split for a given node of a given tree in a given forest in the ``ForestContainer``. + Returns ``np.Inf`` if the node is a leaf or a categorical split node. + + forest_num : :obj:`int` + Index of the forest to be queried + tree_num : :obj:`int` + Index of the tree to be queried + node_id : :obj:`int` + Index of the node to be queried + """ + if self.is_leaf_node(forest_num, tree_num, node_id) or self.is_categorical_split_node(forest_num, tree_num, node_id): + return np.Inf + else: + return self.forest_container_cpp.SplitThreshold(forest_num, tree_num, node_id) + + def node_split_categories(self, forest_num: int, tree_num: int, node_id: int) -> np.array: + """ + Array of category indices that define a categorical split for a given node of a given tree in a given forest in the ``ForestContainer``. + Returns ``np.array([np.Inf])`` if the node is a leaf or a numeric split node. + + forest_num : :obj:`int` + Index of the forest to be queried + tree_num : :obj:`int` + Index of the tree to be queried + node_id : :obj:`int` + Index of the node to be queried + """ + if self.is_leaf_node(forest_num, tree_num, node_id) or self.is_numeric_split_node(forest_num, tree_num, node_id): + return np.array([np.Inf]) + else: + return self.forest_container_cpp.SplitCategories(forest_num, tree_num, node_id) + + def node_leaf_values(self, forest_num: int, tree_num: int, node_id: int) -> np.array: + """ + Leaf node value(s) for a given node of a given tree in a given forest in the ``ForestContainer``. + Values are stale if the node is a split node. + + forest_num : :obj:`int` + Index of the forest to be queried + tree_num : :obj:`int` + Index of the tree to be queried + node_id : :obj:`int` + Index of the node to be queried + """ + return self.forest_container_cpp.NodeLeafValues(forest_num, tree_num, node_id) + + def num_nodes(self, forest_num: int, tree_num: int) -> int: + """ + Number of nodes in a given tree in a given forest in the ``ForestContainer``. + + forest_num : :obj:`int` + Index of the forest to be queried + tree_num : :obj:`int` + Index of the tree to be queried + """ + return self.forest_container_cpp.NumNodes(forest_num, tree_num) + + def num_leaves(self, forest_num: int, tree_num: int) -> int: + """ + Number of leaves in a given tree in a given forest in the ``ForestContainer``. + + forest_num : :obj:`int` + Index of the forest to be queried + tree_num : :obj:`int` + Index of the tree to be queried + """ + return self.forest_container_cpp.NumLeaves(forest_num, tree_num) + + def num_leaf_parents(self, forest_num: int, tree_num: int) -> int: + """ + Number of leaf parents in a given tree in a given forest in the ``ForestContainer``. + + forest_num : :obj:`int` + Index of the forest to be queried + tree_num : :obj:`int` + Index of the tree to be queried + """ + return self.forest_container_cpp.NumLeafParents(forest_num, tree_num) + + def num_split_nodes(self, forest_num: int, tree_num: int) -> int: + """ + Number of split_nodes in a given tree in a given forest in the ``ForestContainer``. + + forest_num : :obj:`int` + Index of the forest to be queried + tree_num : :obj:`int` + Index of the tree to be queried + """ + return self.forest_container_cpp.NumSplitNodes(forest_num, tree_num) + + def nodes(self, forest_num: int, tree_num: int) -> np.array: + """ + Array of node indices in a given tree in a given forest in the ``ForestContainer``. + + forest_num : :obj:`int` + Index of the forest to be queried + tree_num : :obj:`int` + Index of the tree to be queried + """ + return self.forest_container_cpp.Nodes(forest_num, tree_num) + + def leaves(self, forest_num: int, tree_num: int) -> np.array: + """ + Array of leaf indices in a given tree in a given forest in the ``ForestContainer``. + + forest_num : :obj:`int` + Index of the forest to be queried + tree_num : :obj:`int` + Index of the tree to be queried + """ + return self.forest_container_cpp.Leaves(forest_num, tree_num) \ No newline at end of file