{ "cells": [ { "cell_type": "code", "execution_count": 1, "id": "f07fa1f2-187e-4ce0-af95-31d6120977fe", "metadata": { "ExecuteTime": { "end_time": "2024-10-31T11:15:47.604112Z", "start_time": "2024-10-31T11:15:47.544336Z" }, "pycharm": { "name": "#%%\n" } }, "outputs": [], "source": [ "from mimic.utilities import *\n", "\n", "from mimic.model_infer.infer_gLV_bayes import *\n", "from mimic.model_infer import *\n", "from mimic.model_simulate import *\n", "from mimic.model_simulate.sim_gLV import *\n", "\n", "import pandas as pd\n", "import numpy as np\n", "import seaborn as sns\n", "import matplotlib.pyplot as plt\n", "\n", "import arviz as az\n", "import pymc as pm\n", "import pytensor\n", "import pytensor.tensor as at\n", "import pickle\n", "import cloudpickle\n" ] }, { "cell_type": "markdown", "id": "82eb9f01", "metadata": { "pycharm": { "name": "#%% md\n" } }, "source": [ "## Used Bayesian inference to infer the parameters of a (linearised) gLV model\n", "\n", "The generalized Lotka-Volterra equation takes the form\n", "\n", "$$ \\frac{dX_i}{dt} = \\mu_i X_i + X_i M_{ij} X_j + X_i \\epsilon_{il} u_l $$\n", "\n", "where:\n", "- $X_i$ is the concentration of a species\n", "- $\\mu_i$ is its specific growth rate\n", "- $M_{ij}$ is the effect of the interaction of species $i$ on species $j$\n", "- $\\epsilon_{il}$ is the susceptibility to the time-dependent perturbation $u_l$" ] }, { "cell_type": "markdown", "id": "b4324950", "metadata": {}, "source": [ "### Bayesian inference with no shrinkage " ] }, { "cell_type": "code", "execution_count": 2, "id": "ed7e4c8c", "metadata": { "ExecuteTime": { "end_time": "2024-10-31T11:17:37.557214Z", "start_time": "2024-10-31T11:16:01.339288Z" } }, "outputs": [ { "ename": "FileNotFoundError", "evalue": "[Errno 2] No such file or directory: 'params-s5.pkl'", "output_type": "error", "traceback": [ "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m", "\u001b[0;31mFileNotFoundError\u001b[0m Traceback (most recent call last)", "Cell \u001b[0;32mIn[2], line 3\u001b[0m\n\u001b[1;32m 1\u001b[0m \u001b[38;5;66;03m# read in pickled simulated parameters, mu, M, epsilon\u001b[39;00m\n\u001b[1;32m 2\u001b[0m num_species \u001b[38;5;241m=\u001b[39m \u001b[38;5;241m5\u001b[39m\n\u001b[0;32m----> 3\u001b[0m \u001b[38;5;28;01mwith\u001b[39;00m \u001b[38;5;28;43mopen\u001b[39;49m\u001b[43m(\u001b[49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[38;5;124;43mparams-s5.pkl\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[38;5;124;43mrb\u001b[39;49m\u001b[38;5;124;43m\"\u001b[39;49m\u001b[43m)\u001b[49m \u001b[38;5;28;01mas\u001b[39;00m f:\n\u001b[1;32m 4\u001b[0m params \u001b[38;5;241m=\u001b[39m pickle\u001b[38;5;241m.\u001b[39mload(f)\n\u001b[1;32m 5\u001b[0m M \u001b[38;5;241m=\u001b[39m params[\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mM\u001b[39m\u001b[38;5;124m\"\u001b[39m]\n", "File \u001b[0;32m/opt/anaconda3/envs/my_mimic_env/lib/python3.10/site-packages/IPython/core/interactiveshell.py:324\u001b[0m, in \u001b[0;36m_modified_open\u001b[0;34m(file, *args, **kwargs)\u001b[0m\n\u001b[1;32m 317\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m file \u001b[38;5;129;01min\u001b[39;00m {\u001b[38;5;241m0\u001b[39m, \u001b[38;5;241m1\u001b[39m, \u001b[38;5;241m2\u001b[39m}:\n\u001b[1;32m 318\u001b[0m \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mValueError\u001b[39;00m(\n\u001b[1;32m 319\u001b[0m \u001b[38;5;124mf\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mIPython won\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mt let you open fd=\u001b[39m\u001b[38;5;132;01m{\u001b[39;00mfile\u001b[38;5;132;01m}\u001b[39;00m\u001b[38;5;124m by default \u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m 320\u001b[0m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mas it is likely to crash IPython. If you know what you are doing, \u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m 321\u001b[0m \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124myou can use builtins\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124m open.\u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m 322\u001b[0m )\n\u001b[0;32m--> 324\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43mio_open\u001b[49m\u001b[43m(\u001b[49m\u001b[43mfile\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43margs\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mkwargs\u001b[49m\u001b[43m)\u001b[49m\n", "\u001b[0;31mFileNotFoundError\u001b[0m: [Errno 2] No such file or directory: 'params-s5.pkl'" ] } ], "source": [ "# read in pickled simulated parameters, mu, M, epsilon, created in examples-sim-gLV.ipynb\n", "num_species = 5\n", "with open(\"params-s5.pkl\", \"rb\") as f:\n", " params = pickle.load(f)\n", "M = params[\"M\"]\n", "mu = params[\"mu\"]\n", "epsilon = params[\"epsilon\"]\n", "\n", "# read in the data\n", "num_timecourses = 1\n", "data = pd.read_csv(\"data-s5-r1.csv\")\n", "times = data.iloc[:, 0].values\n", "\n", "yobs = data.iloc[:, 1:6].values\n", "\n", "X, F = linearize_time_course_16S(yobs, times)\n", "\n", "# Define priors\n", "prior_mu_mean = 1.0 \n", "prior_mu_sigma = 0.5\n", "\n", "## NB prior_Mii_mean is 0, so not defined as an argument\n", "prior_Mii_mean = 0.0\n", "prior_Mii_sigma = 0.1\n", "\n", "prior_Mij_sigma = 0.1\n", "\n", "\n", "# Sampling conditions\n", "draws = 500\n", "tune = 500\n", "chains = 4\n", "cores = 4\n", "\n", "inference = infergLVbayes()\n", "\n", "inference.set_parameters(X=X, F=F, prior_mu_mean=prior_mu_mean, prior_mu_sigma=prior_mu_sigma,\n", " prior_Mii_sigma=prior_Mii_sigma, prior_Mii_mean=prior_Mii_mean,\n", " prior_Mij_sigma=prior_Mij_sigma,\n", " draws=draws, tune=tune, chains=chains,cores=cores)\n", "\n", "idata = inference.run_inference()\n", "\n", "# To plot posterior distributions\n", "#inference.plot_posterior(idata)\n", "\n", "\n", "summary = az.summary(idata, var_names=[\"mu_hat\", \"M_ii_hat\", \"M_ij_hat\", \"M_hat\", \"sigma\"])\n", "print(summary[[\"mean\", \"sd\", \"r_hat\"]])\n", "\n", "# Save posterior samples to file\n", "az.to_netcdf(idata, 'model_posterior.nc')\n", "\n", "# get median mu_hat and M_hat \n", "mu_h = np.median(idata.posterior[\"mu_hat\"].values, axis=(0,1) ).reshape(-1)\n", "M_h = np.median(idata.posterior[\"M_hat\"].values, axis=(0,1) )\n", "\n", "# compare fitted with simulated parameters\n", "compare_params(mu=(mu, mu_h), M=(M, M_h))\n", "\n", "predictor = sim_gLV(num_species=num_species, M=M_h.T, mu=mu_h)\n", "yobs_h, _, _, _, _ = predictor.simulate(times=times, init_species=yobs[0])\n", "plot_fit_gLV(yobs, yobs_h, times)\n" ] }, { "cell_type": "code", "execution_count": 3, "id": "7fcb005f031b1ce7", "metadata": { "ExecuteTime": { "end_time": "2024-10-31T11:20:49.557644Z", "start_time": "2024-10-31T11:20:49.420375Z" }, "collapsed": false }, "outputs": [], "source": [ "# read in pickled simulated parameters, mu, M, epsilon\n", "num_species = 5\n", "with open(\"params-s5.pkl\", \"rb\") as f:\n", " params = pickle.load(f)\n", "M = params[\"M\"]\n", "mu = params[\"mu\"]\n", "epsilon = params[\"epsilon\"]\n", "\n", "# read in the data\n", "num_timecourses = 3\n", "data = pd.read_csv(\"data-s5-r3.csv\")\n", "times = data.iloc[:, 0].values\n", "\n", "yobs_1 = data.iloc[:, 1:(num_species+1)].values\n", "yobs_2 = data.iloc[:, (num_species+1):(2*num_species+1)].values\n", "yobs_3 = data.iloc[:, (2*num_species+1):(3*num_species+1)].values\n", "ryobs = np.array([yobs_1, yobs_2, yobs_3])\n", "\n", "\n", "X = np.array([], dtype=np.double).reshape(0, num_species+1)\n", "F = np.array([], dtype=np.double).reshape(0, num_species)\n", "\n", "\n", "\n", "for timecourse_idx in range(num_timecourses):\n", " Xs, Fs = linearize_time_course_16S(ryobs[timecourse_idx], times)\n", " X = np.vstack([X, Xs])\n", " F = np.vstack([F, Fs])\n", " \n", "init_species = ryobs[timecourse_idx,0,:] \n" ] }, { "cell_type": "code", "execution_count": null, "id": "d13d05a1", "metadata": { "ExecuteTime": { "end_time": "2024-10-31T11:21:52.005379Z", "start_time": "2024-10-31T11:20:50.711980Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "X shape: (297, 6)\n", "F shape: (297, 5)\n", "Number of species: 5\n", "AdvancedSetSubtensor.0\n" ] }, { "name": "stderr", "output_type": "stream", "text": [ "Initializing NUTS using jitter+adapt_diag...\n", "Multiprocess sampling (4 chains in 4 jobs)\n", "NUTS: [sigma, mu_hat, M_ii_hat_p, M_ij_hat]\n" ] }, { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "0b6029cf4bee4f868421af968189d4bd", "version_major": 2, "version_minor": 0 }, "text/plain": [ "Output()" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
\n" ], "text/plain": [] }, "metadata": {}, "output_type": "display_data" }, { "name": "stderr", "output_type": "stream", "text": [ "Sampling 4 chains for 500 tune and 500 draw iterations (2_000 + 2_000 draws total) took 37 seconds.\n", "/Users/chaniaclare/Documents/GitHub/MIMIC/venv/lib/python3.10/site-packages/arviz/stats/diagnostics.py:596: RuntimeWarning: invalid value encountered in scalar divide\n", " (between_chain_variance / within_chain_variance + num_samples - 1) / (num_samples)\n", "/Users/chaniaclare/Documents/GitHub/MIMIC/venv/lib/python3.10/site-packages/arviz/stats/diagnostics.py:596: RuntimeWarning: invalid value encountered in scalar divide\n", " (between_chain_variance / within_chain_variance + num_samples - 1) / (num_samples)\n", "/Users/chaniaclare/Documents/GitHub/MIMIC/venv/lib/python3.10/site-packages/arviz/stats/diagnostics.py:596: RuntimeWarning: invalid value encountered in scalar divide\n", " (between_chain_variance / within_chain_variance + num_samples - 1) / (num_samples)\n", "/Users/chaniaclare/Documents/GitHub/MIMIC/venv/lib/python3.10/site-packages/arviz/stats/diagnostics.py:596: RuntimeWarning: invalid value encountered in scalar divide\n", " (between_chain_variance / within_chain_variance + num_samples - 1) / (num_samples)\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "mu_hat/mu:\n", "[1.3832308 0.74120724 1.77738998 0.95705513 0.80500161]\n", "[1.27853844 0.55683415 2.06752757 0.86387608 0.70448068]\n", "\n", "M_hat/M:\n", "[[-0.06 -0. 0. -0. 0. ]\n", " [ 0. -0.09 0. 0. 0. ]\n", " [-0.03 -0.01 -0.13 -0.01 0. ]\n", " [-0. 0.04 0. -0.01 0. ]\n", " [ 0. 0.01 0. 0. -0.15]]\n", "\n", " [[-0.05 0. -0.025 0. 0. ]\n", " [ 0. -0.1 0. 0.05 0. ]\n", " [ 0. 0. -0.15 0. 0. ]\n", " [ 0. 0. 0. -0.01 0. ]\n", " [ 0.02 0. 0. 0. -0.2 ]]\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAGwCAYAAABVdURTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAArlElEQVR4nO3de3RU5b3/8c9MIAmxSbjmAkSIErFyCYiCwWqCgJRSC6enHPScH3Co0qOHeyosLi7A06OxqGCpVEAXTY83qBdiFyIaAwleqBRCloAtIlKghQTUksQIScns3x9ppkzIZWYymT3z5P1aa1ay93727O/j45iPez+zt8OyLEsAAACGcNpdAAAAQCARbgAAgFEINwAAwCiEGwAAYBTCDQAAMArhBgAAGIVwAwAAjNLB7gKCzeVy6fTp04qNjZXD4bC7HAAA4AXLslRZWamePXvK6Wz+3Ey7CzenT59WSkqK3WUAAAA/nDp1Sr179262TbsLN7GxsZLq/uHExcXZXA0AAPBGRUWFUlJS3H/Hm9Puwk39pai4uDjCDQAAYcabKSVMKAYAAEYh3AAAAKMQbgAAgFHa3ZwbAAAuV1tbq7///e92lwFJkZGRLX7N2xuEGwBAu2RZlkpLS3X+/Hm7S8E/OJ1OpaamKjIyslXvQ7gBALRL9cEmISFBMTEx3NjVZvU32T1z5oyuvvrqVo0H4QYA0O7U1ta6g023bt3sLgf/0KNHD50+fVqXLl1Sx44d/X4fJhQDANqd+jk2MTExNleCy9VfjqqtrW3V+xBuAADtVqAuRa0tOKrUxW9qbcHRgLxfexWo8eCyFAAArbC24KhW538qSe6fc0en2VlSu8eZGwDhr2iVtLKzVPR4g+VVdlaFduDyYFNvdf6nnMGxGeEGQHgrWiXtekSSJe36X+k3P7hs+RECDtpMY8GmnkkBZ+XKlRoyZIjdZfiEcAMgvO161HP5eFHz24EAaC7Y1DMp4LTWn//8ZzkcDpWUlATleIQbAOFt1NIWti8LTh1oN7wJNvUIOPYg3AAIb5mLpNTMxrddkyVlLgxqOTDfGi+Djb/tW5KVlaU5c+Zo/vz56tKlixITE/Xss8+qqqpKM2bMUGxsrPr166e33npLkpSbm6vOnTt7vEdeXp7P30x6/vnn1bdvX8XHx+vuu+9WZWWle9uOHTv0ne98R507d1a3bt30/e9/X8eOHXNvT01NlSQNHTpUDodDWVlZ/nXeS4QbAOGtaNWVl6LqfV74z0nGQIAsGHtdm7b3xm9+8xt1795de/fu1Zw5c/TAAw9o8uTJGjlypIqLi3XnnXdq6tSp+uabbwJyvGPHjikvL0/btm3Ttm3bVFRUpMcee8y9vaqqStnZ2dq3b58KCgrkdDr1L//yL3K5XJKkvXv3SpLeffddnTlzRq+//npA6moK4QZAeGtpTs2uR4JTB9qNuaPTlO1lYMkee12bfC08PT1dDz30kNLS0rRkyRJFR0ere/fumjlzptLS0rR8+XJ9+eWX+vjjjwNyPJfLpdzcXA0cOFC33Xabpk6dqoKCAvf2f/3Xf9UPf/hD9evXT0OGDNGmTZt08OBBffLJJ5Lq7jwsSd26dVNSUpK6du0akLqaQrgBEN5anHPTwnbAD94EnLYKNpI0ePBg9+8RERHq1q2bBg0a5F6XmJgoSTp79mxAjte3b1/Fxsa6l5OTkz3e++jRo7rnnnt0zTXXKC4uTn379pUknTx5MiDH9xXhBkB4y1z0j0nDl80fuH1h3fKoZXXbgTbQXMBpy2Aj6YrnLjkcDo919fNpXC6XnE6nLMvyaF//+InWHK/+kpMk3XXXXfrqq6/07LPP6qOPPtJHH30kSaqpqfHpOIHCHYoBhL/MRVLGLOnRnnXL31kg3fGQvTUh7FiWJevCBZ/2mZ3RSzUXq/X0eyfc6xZkpWp2Ri+5fJjv4ujUqc2eSt6jRw9VVlaqqqpKV111lSQF9CvZX375pY4cOaJnn31Wt912myTp/fff92gTqGdGeYtwAwCAJOvCBR25cZjP+02QVN5/jF64fpz+35/e1p157+qIj+/Rv3i/HG30EM8RI0YoJiZGS5cu1dy5c/XRRx8pNzc3YO/fpUsXdevWTRs3blRycrJOnjypxYsXe7RJSEhQp06dtGPHDvXu3VvR0dGKj48PWA0N2XpZKicnRzfffLNiY2OVkJCgSZMm6ciRlv+VeOWVV3T99dcrOjpagwYN0vbt24NQLQAAjfv3I+9q+xsL9e9H3rW7lCt07dpVL7zwgrZv365Bgwbp5Zdf1sqVKwP2/k6nU5s3b9b+/fs1cOBALViwQI8/7vktxQ4dOmjt2rXasGGDevbsqYkTJwbs+I1xWA0vxAXRd7/7Xd199926+eabdenSJS1dulSHDh3SJ5984j511tCHH36o22+/XTk5Ofr+97+vl156ST//+c9VXFysgQMHtnjMiooKxcfHq7y8XHFxcYHuEgC71FT987LU0tNSZOP/DQEk6eLFizp+/LhSU1MVHR0tyb/LUoHSlpelwklj41LPl7/ftoabhs6dO6eEhAQVFRXp9ttvb7TNlClTVFVVpW3btrnX3XLLLRoyZIjWr1/f4jEIN4ChCDfwQXN/RGGfQIWbkPq2VHl5uSQ1+/33PXv2aMyYMR7rxo0bpz179jTavrq6WhUVFR4vAABQZ8CAAfrWt77V6OvFF1+0uzy/hMyEYpfLpfnz5+vWW29t9vJSaWmp+/v79RITE1VaWtpo+5ycHD388MMBrRUAAFNs3769ya+GN/x7Gy5CJtzMmjVLhw4duuLrY621ZMkSZWdnu5crKiqUkpIS0GMAABCu+vTpY3cJARcS4Wb27Nnatm2bdu/erd69ezfbNikpSWVlZR7rysrKlJSU1Gj7qKgoRUVFBaxWAAAQ2mydc2NZlmbPnq2tW7dq586d7qeGNicjI8PjeRaSlJ+fr4yMjLYqEwAAhBFbz9zMmjVLL730kt544w3Fxsa6583Ex8erU6dOkqRp06apV69eysnJkSTNmzdPmZmZevLJJzVhwgRt3rxZ+/bt08aNG23rBwAACB22nrl55plnVF5erqysLCUnJ7tfW7Zscbc5efKkzpw5414eOXKkXnrpJW3cuFHp6el69dVXlZeX59U9bgAAgPlsPXPjzS12CgsLr1g3efJkTZ48uQ0qAgDAB0WrpF2P/uMhrQsvW17KQ1ttFBITigEACDtFq6Rdj9T9vut/pT+/Jx0v+sfyP9YTcGwRUjfxAwAgbOx61HO5Ptg0tR1BQ7gBAMAfo5a2sH1ZcOrAFQg3AAD4I3ORlJrZ+LZrsurm4LSBrKwszZkzR/Pnz1eXLl2UmJioZ599VlVVVZoxY4ZiY2PVr18/vfXWW5Kk3Nxcde7c2eM98vLyjH5QJ+EGAAB/FK268lJUvc8LpaLH2+zQv/nNb9S9e3ft3btXc+bM0QMPPKDJkydr5MiRKi4u1p133qmpU6fqm2++abMaQhnhBgAAf7Q0p6Z+UnEbSE9P10MPPaS0tDQtWbJE0dHR6t69u2bOnKm0tDQtX75cX375pT7++OM2qyGUEW4AAPBHi3NuWtjeCoMHD3b/HhERoW7dumnQoEHudfUPvDx79myb1RDKCDcAAPgjc9E/Jg1fNnfl9oV1y6OWtenXwDt27Oix7HA4PNbVz6dxuVxyOp1X3FeuqaeAm4L73AAAoLoby1oXLvi2082zpfQfS0/0q1u+8X7pluy6332Y7+Lo1KnNJvj26NFDlZWVqqqq0lVXXSVJKikpaZNjhQrCDQAAkqwLF3TkxmF+7p1c9+PV2/zau3/xfjliYvw8dvNGjBihmJgYLV26VHPnztVHH32k3NzcNjlWqOCyFAAABuvatateeOEFbd++XYMGDdLLL7+slStX2l1Wm3JY3jzgySAVFRWKj49XeXm54uLi7C4HQKDUVEmP9qz7felpKfIqe+tBSLt48aKOHz+u1NRURUdHS/LzslSAtOVlqXDS2LjU8+XvN5elAABQ3STctro0hODishQAADAK4QYAABiFcAMAaLfa2bTTkBeo8SDcAADanfob3rXXZy+FqpqaGkl1d11uDSYUAwDanYiICHXu3Nn9eIKYmBi+rWQzl8ulc+fOKSYmRh06tC6eEG4AAO1SUlKSpPb7/KVQ5HQ6dfXVV7c6aBJuAADtksPhUHJyshISEox/1lK4iIyMlNPZ+hkzhBsAQLsWERHR6jkeCC1MKAYAAEYh3AAAAKMQbgAY51eFx+wuAYCNCDcAjHB5oPnlzqNaW3DUxmoA2IlwAyDsrS04ql/u9Awzq/M/JeAA7RThBkBYW1twVKvzP210GwEHaJ8INwDCVnPBph4BB2h/CDcAwpI3waYeAQdoXwg3AMLSGi+Djb/tAYQvwg2AsLRg7HVt2h5A+CLcAAhLc0enKdvLwJI99jrNHZ3WxhUBCBWEGwBhy5uAQ7AB2h8enAkgrNUHl2fyP75iG8EGaJ8INwDC3tzRabJqqqSP/rmOYAO0X1yWAmCEmbelun+/P/Nagg3QjhFuABjn8qADoP0h3AAAAKMQbgAAgFEINwAAwCiEGwAAYBTCDQAAMArhBgAAGIVwAwAAjEK4AQAARiHcAAAAoxBuAACAUQg3AADAKIQbAABgFMINAAAwCuEGAAAYhXADAACMQrgBAABGIdwAAACjEG4AAIBRCDcAAMAohBsAAGAUwg0AADAK4QYAABiFcAMAAIxCuAEAAEYh3AAAAKMQbgAAgFEINwAAwCiEGwAAYBTCDQAAMArhBgAAGIVwAwAAjEK4AQAARiHcAAAAoxBuAACAUWwNN7t379Zdd92lnj17yuFwKC8vr9n2hYWFcjgcV7xKS0uDUzAAAAh5toabqqoqpaena926dT7td+TIEZ05c8b9SkhIaKMKAQBAuOlg58HHjx+v8ePH+7xfQkKCOnfu7FXb6upqVVdXu5crKip8Ph4AAAgfYTnnZsiQIUpOTtbYsWP1wQcfNNs2JydH8fHx7ldKSkqQqgQAAHYIq3CTnJys9evX67XXXtNrr72mlJQUZWVlqbi4uMl9lixZovLycvfr1KlTQawYAAAEm62XpXzVv39/9e/f3708cuRIHTt2TGvWrNHzzz/f6D5RUVGKiooKVokAAMBmYXXmpjHDhw/XZ599ZncZAAAgRIR9uCkpKVFycrLdZQAAgBBh62Wpr7/+2uOsy/Hjx1VSUqKuXbvq6quv1pIlS/TXv/5V//d//ydJeuqpp5SamqoBAwbo4sWLeu6557Rz50698847dnUBAACEGFvDzb59+zRq1Cj3cnZ2tiRp+vTpys3N1ZkzZ3Ty5En39pqaGv30pz/VX//6V8XExGjw4MF69913Pd4DAAC0bw7Lsiy7iwimiooKxcfHq7y8XHFxcXaXAyBAvvm6XDFPXF33+4MnFfOteJsrAhBIvvz9Dvs5NwAAAJcj3AAAAKMQbgAAgFEINwAAwCiEGwAAYBTCTQCtLTiq1MVvam3BUbtLAQCg3QqrZ0uFsrUFR7U6/1NJcv+cOzrNzpIAAGiXOHMTAJcHm3qr8z/lDA4AADYg3LRSY8GmHgEHAIDgI9y0QnPBph4BBwCA4CLc+MmbYFOPgAMAQPAQbvy0xstg4297AADgH8KNnxaMva5N2wMAAP8Qbvw0d3Sasr0MLNljr+Nr4QAABAnhphW8CTgEGwAAgotw00rNBRyCDQAAwUe4CYC5o9M0545+HusINgAA2INwEyAPZF3r/n3OHf0INgAA2IRw0wYuDzoAACC4CDcAAMAohBsAAGAUwg0AADAK4QYAABiFcAMAAIxCuAmEolXq9Gh3zY7YKknq8P4T0srOUtEqe+sCAKAd6mB3AWGvaJW06xE5JD3Y8RXd4vxEkbsP123b9Ujdz8xFtpUHAEB7w5mb1tr1qMfirc7DzW4HAABti3DTWqOWeiw6HA23LwteLQAAgHDTapmLpNRMWY1tuyZLylwY5IIAAGjfCDetVbRKOl6khidsJEmfF0pFjwe5IAAA2jfCTWu1NKemflIxAAAICsJNazWYc+PzdgAAEFCEm9bKXCSNWibrsgtTNbf+VJKjbjIxXwMHACCoCDeBkLlIFx484V68lDFPWnmeYAMAgA0INwAAwCiEGwAAYBTCDQAAMArhBgAAGIVwAwAAjEK4AQAARiHcAAAAoxBuAACAUQg3AADAKIQbAABgFMINAAAwCuEGAAAYhXADAACMQrgBAABGIdwAAACjEG4AAIBRCDcAAMAohBsAAGAUwg0AADAK4QYAABilg7cNu3bt6tMbOxwOFRcXq0+fPj4XBQAA4C+vw8358+f11FNPKT4+vsW2lmXpv//7v1VbW9uq4gAAAHzldbiRpLvvvlsJCQletZ0zZ45fBQEAALSG1+HG5XL59MaVlZU+FwMAANBaTCgGAABG8frMze9+9zuNHz9eHTt21O9+97tm2/7gBz9odWEAAAD+8DrcTJo0SaWlpUpISNCkSZOabOdwOJhIDAAAbOPXnBtf598AAAAEC3NuAACAUbwON2vXrtXFixe9fuP169fzjSkAABB0XoebBQsW+BRWFi1apHPnzvlVFAAAgL+8nnNjWZZGjx6tDh282+XChQt+FwUAAOAvr8PNihUrfHrjiRMn+vw8KgAAgNZqs3ADAABgB74tBQAAjOLTgzPrOZ1OORyOJrdzEz8AAGAXv87cbN26Va+//rr7tWXLFi1evFjJycnauHGj1++ze/du3XXXXerZs6ccDofy8vJa3KewsFA33nijoqKi1K9fP+Xm5vrTBQAAYCi/ztxMnDjxinU/+tGPNGDAAG3ZskX33nuvV+9TVVWl9PR0/fjHP9YPf/jDFtsfP35cEyZM0P33368XX3xRBQUFuu+++5ScnKxx48b53A8AAGAev8JNU2655Rb95Cc/8br9+PHjNX78eK/br1+/XqmpqXryySclSd/+9rf1/vvva82aNU2Gm+rqalVXV7uXKyoqvD4eAAAIPwGbUHzhwgWtXbtWvXr1CtRbXmHPnj0aM2aMx7px48Zpz549Te6Tk5Oj+Ph49yslJaXN6gMAAPbz68xNly5dPCYUW5alyspKderUSS+++GLAimuotLRUiYmJHusSExNVUVGhCxcuqFOnTlfss2TJEmVnZ7uXKyoqCDgAEKLWFhzVmvxPtWDsdZo7Os3uchCm/Ao3a9as8Qg3TqdTPXr00IgRI9SlS5eAFRcIUVFRioqKsrsMAEAL1hYc1er8TyXJ/ZOAA3/4FW7+8z//UxcvXtTHH3+ss2fPyuVyqaamRu+9954k6Qc/+EFAi6yXlJSksrIyj3VlZWWKi4tr9KwNACA8XB5s6hFw4C+/ws2OHTs0bdo0ffnll7Isy2Obw+Fos/vcZGRkaPv27R7r8vPzlZGR0SbHAwC0vcaCTT0CDvzh14TiOXPmaPLkyTp9+rRcLpfHy5dg8/XXX6ukpEQlJSWS6r7qXVJSopMnT0qqmy8zbdo0d/v7779fn3/+uRYtWqQ//elP+tWvfqXf/va3WrBggT/dAADYrLlgU291/qdaW3A0SBXBBH6Fm7KyMmVnZ18xuddX+/bt09ChQzV06FBJUnZ2toYOHarly5dLks6cOeMOOpKUmpqqN998U/n5+UpPT9eTTz6p5557jnvcAEAY8ibY1CPgwBd+XZb60Y9+pMLCQl177bWtOnhWVtYVl7Uu19jdh7OysnTgwIFWHRcAYL81Xgaby9tzeQre8CvcPP3005o8ebLee+89DRo0SB07dvTYPnfu3IAUBwAw14Kx13l95qa+PeANv8LNyy+/rHfeeUfR0dEqLCz0+Fq4w+Eg3AAAWlR/FsabgJPNfW/gA7/CzbJly/Twww9r8eLFcjoDdpNjAEA7403AIdjAV34lk5qaGk2ZMoVgAwBotbmj05TdxCUngg384Vc6mT59urZs2RLoWgAA7dTc0Wmac0c/j3UEG/jLr8tStbW1WrVqld5++20NHjz4ignFq1evDkhxAID244Gsa/XLnZ9Jkubc0Y9gA7/5FW4OHjzovjfNoUOHPLZdPrkYAAB/PJDVuluNoH3zK9zs2rUr0HUAAAAEBDOCAQCAUQg3AADAKIQbAABgFMINAMB+RavU6dHumh2xVZLU4f0npJWdpaJV9taFsOTXhGIAAAKmaJW06xE5JD3Y8RXd4vxEkbsP123b9Ujdz8xFtpWH8MOZGwCAvXY96rF4q/Nws9uBlhBuAAD2GrXUY/GK26WNWha8WmAEwg0AwF6Zi6TUTFmNbbsmS8pcGOSCEO4INwAAexWtko4XqdH7239eKBU9HuSCEO4INwAAe7U0p6Z+UjHgJcINAMBeDebc+LwdIWVtwVGlLn5TawuO2lYD4QYAYK/MRdKoZbIuuzBVc+tPJTnqJhPzNfCwsbbgqFbnfypL0ur8T20LOIQbAID9MhfpwoMn3IuXMuZJK88TbMJIfbC5nF0Bh3ADAABapbFgU8+OgEO4AQAAfmsu2NQLdsAh3AAAAL94E2zqBTPgEG4AAIBf1ngZbPxt7y/CDQAA8MuCsde1aXt/EW4AAIBf5o5OU7aXgSV77HWaOzqtjSuqQ7gBAAB+8ybgBDPYSIQbwEMo3FkTAMJNcwEn2MFGItwAbqFyZ00ACEdzR6dpzh39PNbZEWwkwg0gKbTurAkA4eqBrGvdv8+5o58twUYi3AAhd2dNADDB5UEn2Ag3aNdC8c6aAIDWIdyg3QrVO2sCAFqHcIN2K1TvrAkAaB3CDdqtUL2zJgCgdQg3aLdC9c6aABCWilap06PdNTtiqySpw/tPSCs7S0Wrgl5Kh6AfEQgh9YGlubk3BBsAaEHRKmnXI3JIerDjK7rF+Ykidx+u27brkbqfmYuCVg5nbtDuhdqdNQEg7Ox61GPxVufhZre3NcINoNC6syYAhJ1RSz0WHY6G25cFrxYRbgC3ULmzJgCEncxFUmqmrMa2XZMlZS4MajmEG6ARdt5ZEwDCTtEq6XiRGp6wkSR9XigVPR7Ucgg3AACgdVqaU1M/qThICDcAAKB1Gsy58Xl7gBFuAABA62QukkYtk3XZhamaW38qyVE3mTiIXwOXCDcAACAQMhfpwoMn3IuXMuZJK88HPdhIhBsAAGAYwg0AADAK4QYAABiFcAMAAIxCuAEAAEYh3AAAAKMQbgAAgFEINwAAwCiEGwAAYBTCDQAAMArhBgAAGIVwAwAAjEK4AQAARiHcAAAAoxBuAEkqWqVOj3bX7IitkqQO7z8hrewsFa2yty4AgM862F0AYLuiVdKuR+SQ9GDHV3SL8xNF7j5ct23XI3U/MxfZVh4AwDecuQF2PeqxeKvzcLPbAQChjXADjFrqsehwNNy+LHi1AABajXADZC6SUjNlNbbtmiwpc2GQCwIAtAbhBihaJR0vUsMTNpKkzwuloseDXBAAoDUIN0BLc2rqJxUDAMIC4QZoMOfG5+0AgJBCuAEyF0mjlsm67MJUza0/leSom0zM18ABIKwQbgBJylykCw+ecC9eypgnrTxPsAGAMBQS4WbdunXq27evoqOjNWLECO3du7fJtrm5uXI4HB6v6OjoIFYLAABCme3hZsuWLcrOztaKFStUXFys9PR0jRs3TmfPnm1yn7i4OJ05c8b9OnHiRJNtAQBA+2J7uFm9erVmzpypGTNm6IYbbtD69esVExOjTZs2NbmPw+FQUlKS+5WYmBjEigEAQCizNdzU1NRo//79GjNmjHud0+nUmDFjtGfPnib3+/rrr9WnTx+lpKRo4sSJOnz4cJNtq6urVVFR4fECAADmsjXcfPHFF6qtrb3izEtiYqJKS0sb3ad///7atGmT3njjDb3wwgtyuVwaOXKk/vKXvzTaPicnR/Hx8e5XSkpKwPsBAABCh+2XpXyVkZGhadOmaciQIcrMzNTrr7+uHj16aMOGDY22X7JkicrLy92vU6dOBbliAAAQTB3sPHj37t0VERGhsrIyj/VlZWVKSkry6j06duyooUOH6rPPPmt0e1RUlKKiolpdKwAACA+2nrmJjIzUsGHDVFBQ4F7ncrlUUFCgjIwMr96jtrZWBw8eVHJycluVCQAAwoitZ24kKTs7W9OnT9dNN92k4cOH66mnnlJVVZVmzJghSZo2bZp69eqlnJwcSdL//M//6JZbblG/fv10/vx5Pf744zpx4oTuu+8+O7sBAABChO3hZsqUKTp37pyWL1+u0tJSDRkyRDt27HBPMj558qSczn+eYPrb3/6mmTNnqrS0VF26dNGwYcP04Ycf6oYbbrCrCwAAIITYHm4kafbs2Zo9e3aj2woLCz2W16xZozVr1gShKgAAEI7C7ttSAAAAzSHcAAAAoxBuAACAUQg3AADAKIQbAABgFMINAAAwCuEGAAAYhXADAACMQrgBAABGIdwAAACjEG4AAIBRCDcAAMAohBsAAGAUwg0AADAK4QYAABiFcAMAAIxCuAEAAEYh3AAAAKMQbgAAgFEINwAAwCiEGwAAYBTCDQAAMArhBgAAGIVwAwAAjEK4AQAARiHcAAAAoxBuAACAUQg3AADAKIQbAABgFMINAAAwCuEGAAAYhXADAACMQrgBAABGIdwAAACjEG4AAIBRCDcAAMAohBsAAGAUwg0AADAK4QYAABiFcAMAAIxCuAEAAEYh3AAAAKMQbgAAgFEINwAAwCiEGwAAYBTCDQAAMArhBgAAGIVwAwAAjEK4AQAARiHcAAAAoxBuAACAUQg3AADAKIQbAABgFMINAAAwCuEGAAAYhXADAACMQrgBAABGIdwAAACjEG4AAIBRCDcAAMAohBsAAGAUwg0AADAK4QYAABiFcAMAAIxCuAEAAEYh3AAAAKMQbgAAgFEINwAAwCiEGwAAYBTCDQAAMEpIhJt169apb9++io6O1ogRI7R3795m27/yyiu6/vrrFR0drUGDBmn79u1BqhQAAIQ628PNli1blJ2drRUrVqi4uFjp6ekaN26czp4922j7Dz/8UPfcc4/uvfdeHThwQJMmTdKkSZN06NChIFcOAABCUQe7C1i9erVmzpypGTNmSJLWr1+vN998U5s2bdLixYuvaP+LX/xC3/3ud7Vw4UJJ0s9+9jPl5+fr6aef1vr164Na++VcLpdclxySpEtffaVL1ZdsqwX+uVRVwRiGMcYv/DGG4e/yMXS5XLbVYWu4qamp0f79+7VkyRL3OqfTqTFjxmjPnj2N7rNnzx5lZ2d7rBs3bpzy8vIabV9dXa3q6mr3ckVFResLb4Tr/HkdeTW5buHV77XJMRAMjGF4Y/zCH2MY/urGsNdPzktxXWypwNbLUl988YVqa2uVmJjosT4xMVGlpaWN7lNaWupT+5ycHMXHx7tfKSkpgSm+gZir4trkfQEACEd2/l20/bJUW1uyZInHmZ6Kioo2CTjOLl2U9sH7AX9fAADCkbOLPWdtJJvDTffu3RUREaGysjKP9WVlZUpKSmp0n6SkJJ/aR0VFKSoqKjAFN8PpdMrZrVubHwcAADTP1stSkZGRGjZsmAoKCtzrXC6XCgoKlJGR0eg+GRkZHu0lKT8/v8n2AACgfbH9slR2dramT5+um266ScOHD9dTTz2lqqoq97enpk2bpl69eiknJ0eSNG/ePGVmZurJJ5/UhAkTtHnzZu3bt08bN260sxsAACBE2B5upkyZonPnzmn58uUqLS3VkCFDtGPHDvek4ZMnT8rp/OcJppEjR+qll17SQw89pKVLlyotLU15eXkaOHCgXV0AAAAhxGFZlmV3EcFUUVGh+Ph4lZeXKy6ObzgBABAOfPn7bfsdigEAAAKJcAMAAIxCuAEAAEYh3AAAAKMQbgAAgFEINwAAwCiEGwAAYBTCDQAAMArhBgAAGMX2xy8EW/0NmSsqKmyuBAAAeKv+77Y3D1Zod+GmsrJSkpSSkmJzJQAAwFeVlZWKj49vtk27e7aUy+XS6dOnFRsbK4fDEdD3rqioUEpKik6dOmXkc6tM759kfh/pX/gzvY/0L/y1VR8ty1JlZaV69uzp8UDtxrS7MzdOp1O9e/du02PExcUZ+y+tZH7/JPP7SP/Cn+l9pH/hry362NIZm3pMKAYAAEYh3AAAAKMQbgIoKipKK1asUFRUlN2ltAnT+yeZ30f6F/5M7yP9C3+h0Md2N6EYAACYjTM3AADAKIQbAABgFMINAAAwCuEGAAAYhXDjo3Xr1qlv376Kjo7WiBEjtHfv3mbbv/LKK7r++usVHR2tQYMGafv27UGq1D++9C83N1cOh8PjFR0dHcRqfbN7927ddddd6tmzpxwOh/Ly8lrcp7CwUDfeeKOioqLUr18/5ebmtnmd/vK1f4WFhVeMn8PhUGlpaXAK9lFOTo5uvvlmxcbGKiEhQZMmTdKRI0da3C+cPoP+9DGcPofPPPOMBg8e7L65W0ZGht56661m9wmn8fO1f+E0do157LHH5HA4NH/+/Gbb2TGGhBsfbNmyRdnZ2VqxYoWKi4uVnp6ucePG6ezZs422//DDD3XPPffo3nvv1YEDBzRp0iRNmjRJhw4dCnLl3vG1f1LdHSjPnDnjfp04cSKIFfumqqpK6enpWrdunVftjx8/rgkTJmjUqFEqKSnR/Pnzdd999+ntt99u40r942v/6h05csRjDBMSEtqowtYpKirSrFmz9Pvf/175+fn6+9//rjvvvFNVVVVN7hNun0F/+iiFz+ewd+/eeuyxx7R//37t27dPd9xxhyZOnKjDhw832j7cxs/X/knhM3YN/eEPf9CGDRs0ePDgZtvZNoYWvDZ8+HBr1qxZ7uXa2lqrZ8+eVk5OTqPt/+3f/s2aMGGCx7oRI0ZY//Vf/9WmdfrL1/79+te/tuLj44NUXWBJsrZu3dpsm0WLFlkDBgzwWDdlyhRr3LhxbVhZYHjTv127dlmSrL/97W9BqSnQzp49a0myioqKmmwTbp/BhrzpYzh/Di3Lsrp06WI999xzjW4L9/GzrOb7F65jV1lZaaWlpVn5+flWZmamNW/evCbb2jWGnLnxUk1Njfbv368xY8a41zmdTo0ZM0Z79uxpdJ89e/Z4tJekcePGNdneTv70T5K+/vpr9enTRykpKS3+H0q4Cafxa40hQ4YoOTlZY8eO1QcffGB3OV4rLy+XJHXt2rXJNuE+ht70UQrPz2Ftba02b96sqqoqZWRkNNomnMfPm/5J4Tl2s2bN0oQJE64Ym8bYNYaEGy998cUXqq2tVWJiosf6xMTEJucolJaW+tTeTv70r3///tq0aZPeeOMNvfDCC3K5XBo5cqT+8pe/BKPkNtfU+FVUVOjChQs2VRU4ycnJWr9+vV577TW99tprSklJUVZWloqLi+0urUUul0vz58/XrbfeqoEDBzbZLpw+gw1528dw+xwePHhQ3/rWtxQVFaX7779fW7du1Q033NBo23AcP1/6F25jJ0mbN29WcXGxcnJyvGpv1xi2u6eCI3AyMjI8/o9k5MiR+va3v60NGzboZz/7mY2VwRv9+/dX//793csjR47UsWPHtGbNGj3//PM2VtayWbNm6dChQ3r//fftLqXNeNvHcPsc9u/fXyUlJSovL9err76q6dOnq6ioqMkAEG586V+4jd2pU6c0b9485efnh/zEZ8KNl7p3766IiAiVlZV5rC8rK1NSUlKj+yQlJfnU3k7+9K+hjh07aujQofrss8/aosSga2r84uLi1KlTJ5uqalvDhw8P+cAwe/Zsbdu2Tbt371bv3r2bbRtOn8HL+dLHhkL9cxgZGal+/fpJkoYNG6Y//OEP+sUvfqENGzZc0TYcx8+X/jUU6mO3f/9+nT17VjfeeKN7XW1trXbv3q2nn35a1dXVioiI8NjHrjHkspSXIiMjNWzYMBUUFLjXuVwuFRQUNHk9NSMjw6O9JOXn5zd7/dUu/vSvodraWh08eFDJycltVWZQhdP4BUpJSUnIjp9lWZo9e7a2bt2qnTt3KjU1tcV9wm0M/eljQ+H2OXS5XKqurm50W7iNX2Oa619DoT52o0eP1sGDB1VSUuJ+3XTTTfqP//gPlZSUXBFsJBvHsE2nKxtm8+bNVlRUlJWbm2t98skn1k9+8hOrc+fOVmlpqWVZljV16lRr8eLF7vYffPCB1aFDB+uJJ56w/vjHP1orVqywOnbsaB08eNCuLjTL1/49/PDD1ttvv20dO3bM2r9/v3X33Xdb0dHR1uHDh+3qQrMqKyutAwcOWAcOHLAkWatXr7YOHDhgnThxwrIsy1q8eLE1depUd/vPP//ciomJsRYuXGj98Y9/tNatW2dFRERYO3bssKsLzfK1f2vWrLHy8vKso0ePWgcPHrTmzZtnOZ1O691337WrC8164IEHrPj4eKuwsNA6c+aM+/XNN9+424T7Z9CfPobT53Dx4sVWUVGRdfz4cevjjz+2Fi9ebDkcDuudd96xLCv8x8/X/oXT2DWl4belQmUMCTc++uUvf2ldffXVVmRkpDV8+HDr97//vXtbZmamNX36dI/2v/3tb63rrrvOioyMtAYMGGC9+eabQa7YN770b/78+e62iYmJ1ve+9z2ruLjYhqq9U//V54av+j5Nnz7dyszMvGKfIUOGWJGRkdY111xj/frXvw563d7ytX8///nPrWuvvdaKjo62unbtamVlZVk7d+60p3gvNNY3SR5jEu6fQX/6GE6fwx//+MdWnz59rMjISKtHjx7W6NGj3X/4LSv8x8/X/oXT2DWlYbgJlTF0WJZlte25IQAAgOBhzg0AADAK4QYAABiFcAMAAIxCuAEAAEYh3AAAAKMQbgAAgFEINwAAwCiEGwAAYBTCDQCjZGVlaf78+XaXAcBG3KEYgFG++uordezYUbGxsXaXAsAmhBsAAGAULksBMAqXpQAQbgAAgFEINwAAwCiEGwAAYBTCDQAAMArhBgAAGIVwAwAAjEK4AQAARuEmfgAAwCicuQEAAEYh3AAAAKMQbgAAgFEINwAAwCiEGwAAYBTCDQAAMArhBgAAGIVwAwAAjEK4AQAARiHcAAAAoxBuAACAUf4/LERKtLh5JOQAAAAASUVORK5CYII=", "text/plain": [ "