{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "Markov jump process: Reaction network\n", "=====================================" ] }, { "cell_type": "markdown", "metadata": { "raw_mimetype": "text/restructuredtext" }, "source": [ "In the following, we fit stochastic chemical reaction kinetics with pyABC and show how to perform model selection between two competing models." ] }, { "cell_type": "raw", "metadata": { "raw_mimetype": "text/restructuredtext" }, "source": [ "\n", "This notebook can be downloaded here: :download:`Markov Jump Process: Reaction Network `." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We consider the two Markov jump process models $m_1$ and $m_2$ for conversion of (chemical) species $X$ to species $Y$:\n", "\n", "$$\n", " m_1: X + Y \\xrightarrow{k_1} 2Y\\\\ m_2: X \\xrightarrow{k_2} Y.\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Each model is equipped with a single rate parameter $k$.\n", "To simulate these models, we define a simple Gillespie simulator:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# install if not done yet\n", "!pip install pyabc --quiet" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "\n", "\n", "def h(x, pre, c):\n", " return (x**pre).prod(1) * c\n", "\n", "\n", "def gillespie(x, c, pre, post, max_t):\n", " \"\"\"\n", " Gillespie simulation\n", "\n", " Parameters\n", " ----------\n", "\n", " x: 1D array of size n_species\n", " The initial numbers.\n", "\n", " c: 1D array of size n_reactions\n", " The reaction rates.\n", "\n", " pre: array of size n_reactions x n_species\n", " What is to be consumed.\n", "\n", " post: array of size n_reactions x n_species\n", " What is to be produced\n", "\n", " max_t: int\n", " Timulate up to time max_t\n", "\n", " Returns\n", " -------\n", " t, X: 1d array, 2d array\n", " t: The time points.\n", " X: The history of the species.\n", " ``X.shape == (t.size, x.size)``\n", "\n", " \"\"\"\n", " t = 0\n", " t_store = [t]\n", " x_store = [x.copy()]\n", " S = post - pre\n", "\n", " while t < max_t:\n", " h_vec = h(x, pre, c)\n", " h0 = h_vec.sum()\n", " if h0 == 0:\n", " break\n", " delta_t = np.random.exponential(1 / h0)\n", " # no reaction can occur any more\n", " if not np.isfinite(delta_t):\n", " t_store.append(max_t)\n", " x_store.append(x)\n", " break\n", " reaction = np.random.choice(c.size, p=h_vec / h0)\n", " t = t + delta_t\n", " x = x + S[reaction]\n", "\n", " t_store.append(t)\n", " x_store.append(x)\n", "\n", " return np.array(t_store), np.array(x_store)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Next, we define the models in terms of ther initial molecule numbers $x_0$, an array ``pre`` which determines what is to be consumed (the left hand side of the reaction equations) and an array ``post`` which determines what is to be produced (the right hand side of the reaction equations).\n", "Moreover, we define that the simulation time should not exceed ``MAX_T`` seconds.\n", "\n", "Model 1 starts with initial concentrations $X=40$ and $Y=3$.\n", "The reaction $X + Y \\rightarrow 2Y$ is encoded in ``pre = [[1, 1]]`` and ``post = [[0, 2]]``." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "MAX_T = 0.1\n", "\n", "\n", "class Model1:\n", " __name__ = \"Model 1\"\n", " x0 = np.array([40, 3]) # Initial molecule numbers\n", " pre = np.array([[1, 1]], dtype=int)\n", " post = np.array([[0, 2]])\n", "\n", " def __call__(self, par):\n", " t, X = gillespie(\n", " self.x0, np.array([float(par[\"rate\"])]), self.pre, self.post, MAX_T\n", " )\n", " return {\"t\": t, \"X\": X}" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Model 2 inherits the initial concentration from model 1.\n", "The reaction $X \\rightarrow Y$ is incoded in ``pre = [[1, 0]]`` and ``post = [[0, 1]]``." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "class Model2(Model1):\n", " __name__ = \"Model 2\"\n", " pre = np.array([[1, 0]], dtype=int)\n", " post = np.array([[0, 1]])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We draw one stochastic simulation from model 1 (the \"Observation\") and and one from model 2 (the \"Competition\") and visualize both" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAs0AAAEWCAYAAACdXqrwAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3de7yVZZn/8c9XILA8ixpIsG2sRhDaThsYiinGA6ZmNNpvSDooU6O9QhObyVM2ouJoNZOl+YsoNZxUMLWxSPMnhJOaURtEURhH80AgxklCJo9w/f5Yz8LFZu+91t5rPetZh+/79dqvvdbzPGuti8Xe17r3fd/XfSsiMDMzMzOzru2WdQBmZmZmZrXOjWYzMzMzsyLcaDYzMzMzK8KNZjMzMzOzItxoNjMzMzMrwo1mMzMzM7Mi3Gi2miRphqQfZR1HT0m6W9KpWcdhZmbFc7KkWZK+Ws2YrH650WyZkXSapOWS/izpBUnflbRP1nGVqrOGfUQcFxFzsorJzCwtkqZIape0VdLapEE6Puu48orl5OQz54EO5z8fEZdVM06rX240WyYk/RPwNeDLwN7AXwPDgHslvaVKMfStxuuYmdU7SV8CvgX8K3AQMBT4v8CkLOMyqyY3mq3qJO0FXAKcFRG/iIjXI+JZ4O+BFuBTyaUDJM2T9JKkpZLeW/Ac50lak5x7QtJRyfHdJJ0v6feSNkq6VdJ+ybkWSSHps5JWAb9MekrO7BDfI5JOSm5/W9IfJG2RtETS3yTHPwxcCExOel0eSY7fJ+lzBbFcJOk5Sesk3Shp7w6xnCpplaQNkr6SyhtuZlaGJG9dCkyLiDsi4n+TvP2ziPiypP6SviXp+eTrW5L6J4+dIGm1pHOTPLhW0sckHS/pfyRtknRhwWvNkHRbN7l/sKTbJa2X9IykLybHu83Jkg4DZgHjkvObk/M/lDSz4Pn/UdJTSVw/lTS44FxI+rykJyVtlnStJKX53lttcaPZsvB+YABwR+HBiNgK3AUckxyaBPwY2A+4GfhPSf0kvQc4ExgdEXsCxwLPJo85C/gY8CFgMPAicG2H1/8QcFjyuFuAU/InJA0n1+P98+TQ74DWghh+LGlARPyCXI/LvIjYIyLey65OS77+FngnsAfwnQ7XjAfeAxwF/EuS2M3Mask4cjn7J12c/wq50cJW4L3AGOCigvNvTx5/MPAvwPfJdY68D/gb4KuSDim4vqvcvxvwM+CR5LmOAqZLOrZYTo6IlcDngYeS87tMBZR0JHAFuQ6cQcBzwNwOl30EGA2MSq47tov3xBqQG82WhYHAhoh4o5Nza5PzAEsi4raIeB34Jrmk+9fANqA/MFxSv4h4NiJ+nzzm88BXImJ1RLwKzAA+3mEqxoykp+Rlch8CrZKGJec+CdyRPJaI+FFEbIyINyLi35PXfU+J/85PAt+MiKeTPwguAD7RIZZLIuLliHiE3AdBZ41vM7Ms7U/XORtyue7SiFgXEevJjSR+uuD868DlSS6fSy7HfzsiXoqIx4EV7Jz7usr9o4EDIuLSiHgtIp4m1wD/RIX+nZ8Ero+IpclnwAXkeqZbCq65MiI2R8QqYBG5PxSsSbjRbFnYAAzsYk7xoOQ8wB/yByNiO7AaGBwRTwHTyTWI10maWzCENgz4STJ0thlYSa6RfVDBaxQ+70vkepXzSfcU4Kb8eUn/LGmlpD8lz7c3bzbqixlMrqci7zmgb4dYXii4/WdyvdFmZrVkI13nbOg81w0uuL8xIrYlt19Ovv+x4PzL7Jz7Os395PL74Hx+T3LyheycU8ux078j6ezYSK5XO885u4m50WxZeAh4FTip8KCkPYDjgIXJoXcUnNsNGAI8DxARN0fEeHJJNMgVFUIu2R4XEfsUfA2IiDUFLxUd4rkFOEVSfghyUfKafwOcS24Ibt9kOO9PgLp4no6eT+LLGwq8wc4fFmZmtS6fsz/WxfnOct3zZbxeV7n/D8AzHfL7nhFxfHJ5sZzco5wt6W3ketnXdPkIaypuNFvVRcSfyA3fXSPpw8lctRbgVnI9Cv+RXPo+SSclvRvTySXt30h6j6Qjk0KTV8j1UmxPHjMLuDw/3ULSAZKKVXffRS5RXkpuPlz+ufYk18hdD/SV9C/AXgWP+yPQkiT1ztwCnCPpkOQPgvx8u66GOM3Mak6Ss/8FuDYp4ntrkrePk/R1crnuoiTfDkyuLWed/U5zP/Bb4CXlCsF3l9RH0uGSRiePK5aT/wgMUdcrNN0CTJXUmny+/CuwOClUN3Oj2bIREV8nN6z2b8AWYDG5XoSj8vOJgTuByeSK+T4NnJTMcesPXEluGscLwIHk5p4BfBv4KfD/JL1ELtGOLRLLq+SKEo8mV3SSdw/wC+B/yA3ZvULBsCG5QhWAjZKWdvLU15P7A+BXwDPJ48/qLhYzs1qU1HR8iVyB33pyufBM4D+BmUA78CiwHFiaHOutTnN/MsXjI+TmET9D7jPgB+SmzUHxnPxL4HHgBUkbOp6MiAXAV4HbydXX/AWVmy9tDUARxUYrzMzMzNInaQZwaER8qti1ZtXmnmYzMzMzsyLcaDYzMzMzK8LTM8zMzMzMinBPs5mZmZlZEV0tVF5TBg4cGC0tLVmHYWbWY0uWLNkQEQdkHUc1OWebWb3qLmfXRaO5paWF9vb2rMMwM+sxSc8Vv6qxOGebWb3qLmd7eoaZmZmZWRFuNJuZmZmZFeFGs5mZmZlZEW40m5mZmZkV4UazmZmZmVkRbjSbmdkOkvpIeljS/OT+IZIWS3pK0jxJb8k6RjOzLLjRbGZmhc4GVhbc/xpwVUQcCrwIfDaTqMzMMlYX6zSbWYNqvwGW35Z1FKV5+0g47sqso0iVpCHACcDlwJckCTgSmJJcMgeYAXy30q99yc8eB+DiE0dU+qmtkdVTDrHqq3Dedk+zmWVn+W3wwvKso7A3fQs4F9ie3N8f2BwRbyT3VwMHd/ZASadLapfUvn79+h6/8Irnt7Di+S29CNmamnOIVZF7ms0sW28fCVN/nnUUTU/SR4B1EbFE0oSePj4iZgOzAdra2qLC4Zl1zTnEqsSNZjMzA/gA8FFJxwMDgL2AbwP7SOqb9DYPAdZkGKOZWWY8PcPMzIiICyJiSES0AJ8AfhkRnwQWAR9PLjsVuDOjEM3MMuWeZjMrLq1imxeW54ZWrZadB8yVNBN4GLgu43jMzDLhRrOZFZcvtql0A/ftI2Hkx4tfZ1UVEfcB9yW3nwbGVON1V6zdwuTvPcSk1oOZMnZoNV7SalFP/kj3H95WRW40m1lpXGxjKZrUmluUY8Xa3AoabjQ3sZ78ke4/vK2K3Gg2M7PMTRk7lCljhzL5ew9lHYrVAv+RbjXIjWYz21XH4VEPgZqZWZPz6hlmtquOGwZ4CNTMzJqce5rNrHMeHrWMLH5mEzcvXuV5zfWut6vueGTLalTqPc2S+kh6WNL85P4hkhZLekrSPElvSTsGMzOrD/mCwDuXeQ+VutfbLa49smU1qho9zWcDK8ntLgXwNeCqiJgraRbwWeC7VYjDzMxq3JSxQ91gbiQesbIGkmqjWdIQ4ATgcuBLkgQcCUxJLpkDzCCFRvMlP3scgItPHFHppzarfeVuRuLhUTMzs52kPT3jW8C5wPbk/v7A5oh4I7m/Gji4swdKOl1Su6T29evX9/iFVzy/hRXPb+lFyGYNoLfDonkeHjUzM9tJaj3Nkj4CrIuIJZIm9PTxETEbmA3Q1tYWFQ7PrPF5WNTqWH53wDzvElgHvFSlNbg0p2d8APiopOOBAeTmNH8b2EdS36S3eQjgyWtmZrZDvhgwz7sE1omOO/l5xMoaTGqN5oi4ALgAIOlp/ueI+KSkHwMfB+YCpwJ3phWDmZnVn/zugHneJbCOeITLGlgW6zSfB8yVNBN4GLgurRfqOLxXDR5CtKrprtjPw6JmZmYVVZVGc0TcB9yX3H4aGJP2a3Yc3qsGDyFaVXUcCi3kYVEzM7OKatgdATsO71WDhxCt6jwUamZmVhWp7whoZmZWrvx0u5sXr8o6FDNrUg3b05yVUudRe+6zlaX9BnjuARg2PutIzFKXn27nKXBmliX3NFfQpNaDGT5or6LXrVi7xdvEWnnyBYCet2xNYMrYocw7Y1xJ+dXMLC3uaa6gUudRe+6zVcSw8dA2NesorEFIGgD8CuhP7rPhtoi4WNIPgQ8Bf0ouPS0ilmUTpZlZdtxoNjMzgFeBIyNiq6R+wAOS7k7OfTkiuljf0MysObjRbGZmREQAW5O7/ZKvyC6izuXrRlwXUgXdrQXfGa8Pbw3Oc5rNzAwASX0kLQPWAfdGxOLk1OWSHpV0laT+XTz2dEntktrXr1+fSnz5uhHXhVRJfi34Unl9eGtw7mk2MzMAImIb0CppH+Ankg4HLgBeAN4CzCa3q+ulnTx2dnKetra2VHqo83UjrgupIq8Fb7aDe5rNzGwnEbEZWAR8OCLWRs6rwA1UYUdXM7Na5EazmZkh6YCkhxlJuwPHAP8taVByTMDHgMeyi9LMLDuenpERF7M0iZ4W0pTKBTdWeYOAOZL6kOtQuTUi5kv6paQDAAHLgM9nGWTe4mc2cfPiVc6fldQxXznPmO3EjeYMeHerJpIvpKn0B48LbqzCIuJR4IhOjh+ZQTjdmtR6MIuf2cSdy9Y4f1ZSx3zlPGO2EzeaM+BilibjQhqzipoydqhXz0iL85VZlzyn2czMzMysCDeazczMzMyK8PSMjOULAnvCxYMZ6U1RnwtpzMzMGoJ7mjOU392qJ7wTVoZ6ujsWuJDGLEX5ToebF6/KOpT61n4D3HBCz/ObWZNxT3OG8gWBPeHiwYy5SMasJngVogoqXDXDf+SbdcmNZjMzqztehajC3CFgVpSnZ5iZmZmZFeGe5jpUWDzoosAqab8BnnsAho3POhIzMzPLgHua60xh8aCLAqsov2qG5/uZ1Zz8ltpmZmlyo7nOTBk7lHlnjGPeGeN6vPKGlWnYeGibmnUUZlYgXxDoDgQzS5sbzWZmVremjB3K2EP2yzoMM2sCntNc5zy/uQJK2bTEm5SYmZk1Nfc01zHPb66QUjYt8fqlZmZmTc09zXWscHMUr1VaJq9RalbXCkfdwCNvXepsZM0jaWYlcU+zmZnVtcJRN/DIW7c6G1nzSJpZSdzTbGZmda1w1A088laUR9bMesWN5gbScXgSPES5k64K/jw0aWZmZkV4ekaD6Dg8CR6i3EVXBX8emjRD0gBJv5X0iKTHJV2SHD9E0mJJT0maJ+ktWcdqZpYF9zQ3iI7Dk+Ahyk55WNKsK68CR0bEVkn9gAck3Q18CbgqIuZKmgV8FvhuloGamWXBjWYzMyMiAtia3O2XfAVwJDAlOT4HmEEdNJrzW2s39fQ0r5RhVlFuNDe4zuY5QxPNdS780PCHhVm3JPUBlgCHAtcCvwc2R8QbySWrgYO7eOzpwOkAQ4dmm1smtR7M4mc2ceeyNc2R57qSn5JWmPc8Hc2s19xobmCTWjv9bGPF2i0AzfFhUvih4Q8Ls25FxDagVdI+wE+Av+zBY2cDswHa2toinQhLM2XsUNdz5HlKmlnFuNHcwDqb5wxNONfZHxpmPRIRmyUtAsYB+0jqm/Q2DwHcGjWzppTa6hmuxDYzqx+SDkh6mJG0O3AMsBJYBOSHaE4F7swmQjOzbKXZ0+xKbDOz+jEImJPMa94NuDUi5ktaAcyVNBN4GLguyyB7omm31s7XcriOw6yiUms0N1olttU4b1xiVpaIeBQ4opPjTwNjqh9ReTrWdDRtLYfrOMwqJtU5zY1SiW11oKteFX9omDWlpt9a27UcZhWXaqO5USqxrU74Q8LMzMxSUpVttCNiM7likh2V2MkpV2KbmZmZWc1LradZ0gHA68nSRflK7K/xZiX2XFyJbWZm1nve9c+satLsaR4ELJL0KPA74N6ImA+cB3xJ0lPA/tRRJXYjyW8xW/fab4AbTsh9SJiZNZt8PUch13KYpSLN1TMaqhK7kTTUFrOuEjezZud6DrOq8I6ATajhtpj1B4aZmZmlrKRGczI/+R+BlsLHRMQ/pBOWmZn1lnO2mVnlldrTfCdwP7AA2JZeOFZN+XnNNT1Fo6tNS/Jc8GLWGedsM7MKK7XR/NaIOC/VSKyq6mZec7GtYD2X2awzztmdyG+rXZfbaXvXU7PMldponi/p+Ii4K9VorGrqal6z5yyb9ZRzdgf5bbXrdjtt73pqlrlSG81nAxdKeg14PTkWEbFXOmGZmVkZnLM7yG+rXdfbabsDwSxTJTWaI2LPtAMxM7PKcM42M6u8kpeck/RR4IPJ3fuSjUqszuXn+JWiqvMA8/P3PF/PrFecs83MKqvUJeeuBEYDNyWHzpb0gYi4ILXILHX5OX6lqPo8QG9aYtZrztndq9mCwO5WC3IHglnmSu1pPh5ojYjtAJLmAA8DTsB1LD/HrxSZzAP0/D2z3nLO7kJNFwR2N7rmDgSzzPVkR8B9gE3J7b1TiMXMzCrHObsTNV8Q6M4Cs5pVaqP5CuBhSYsAkZsnd35qUVlNSmVI02uPmqWhxzlb0juAG4GDgABmR8S3Jc0gt7vg+uTSC72UnZk1o1JXz7hF0n3k5sgBnBcRL6QWldWc1IY0vfaoWcX1Mme/AfxTRCyVtCewRNK9ybmrIuLfUgrXzKwudNtolvSXEfHfkv4qObQ6+T5Y0uCIWJpueFYrUh3S9HCkWUWUk7MjYi2wNrn9kqSVQOnVwnVo8TObuHnxquznNXu1ILO6UKyn+UvA6cC/d3IugCMrHpGZmfVWRXK2pBbgCGAx8AHgTEmfAdrJ9Ua/2MljTk9em6FDa6i4rguTWg9m8TObuHPZmuwbzV4tyKwudNtojojTk5vHRcQrheckDUgtKjMz67FK5GxJewC3A9MjYouk7wKXkWt0X0auQf4Pnbz2bGA2QFtbW/T6H1ElU8YO5c5la7IO400ecTOrebuVeN2vSzxmVlz7DXDDCbmvF5ZnHY1ZI+pVzpbUj1yD+aaIuAMgIv4YEduS5eu+D4ypaKRmZnWi2Jzmt5Ob07a7pCPIVWED7AW8NeXYrFEVDkV6ODJTr7/+OqtXr+aVV14pfrF1a8CAAQwZMoR+/fplFkM5OVuSgOuAlRHxzYLjg5L5zgB/BzxW8cDNrCTO2ZXTm5xdbE7zscBpwBDgmwXHXwIu7GmAZjt4KLImrF69mj333JOWlhZybSbrjYhg48aNrF69mkMOOSTLUMrJ2R8APg0sl7QsOXYhcIqkVnLTM54FzqhgvGbWA87ZldHbnF1sTvMcYI6kkyPi9nKDNLPa8sorrzj5VoAk9t9/f9avX1/84hSVk7Mj4gHe7Jku1NBrMme6goZXzbAecs6ujN7m7FLXab5d0gnACGBAwfFLe/RqZlZznHwro5beR+fs0mS+goZXzbBeqKVcU8968z6WVAgoaRYwGTiLXE/E/wGG9fjVrCHkdwa8efGq0h/k4j/rwuWXX86IESMYNWoUra2tLF68uKLPf/zxx7N58+aynuOLX/wil176Znvz8ssvZ9q0aeWGlhrn7NJMGTuUsYfsl20Q+alqbVOzjcOsRM2cs0vdRvv9ETFK0qMRcYmkfwfuLvvVre70emdAF/9ZJx566CHmz5/P0qVL6d+/Pxs2bOC1116r6GvcdVf5swtmzpxJa2srn/rUpwD4wQ9+wMMPP1z286bIOdvMKq7Zc3apS87lyzT/LGkw8DowqOxXt7ozZexQ5p0xjuGD9ur5g/M9Ku5VscTatWsZOHAg/fv3B2DgwIEMHjwYgJaWFs4991xGjhzJmDFjeOqppwBYv349J598MqNHj2b06NE8+OCDAGzdupWpU6cycuRIRo0axe23377jeTZs2ADAj370I8aMGUNraytnnHEG27ZtY9u2bZx22mkcfvjhjBw5kquuumqXOPfaay8uv/xyzjzzTM4880wuvfRS9tlnn9TfnzI4Z5tZxTV7zi61p/lnkvYBvgEsJVdF/f2yX93MasYlP3ucFc9vqehzDh+8FxefOKLL8xMnTuTSSy/l3e9+N0cffTSTJ0/mQx/60I7ze++9N8uXL+fGG29k+vTpzJ8/n7PPPptzzjmH8ePHs2rVKo499lhWrlzJZZddtuN6gBdf3HnTupUrVzJv3jwefPBB+vXrxxe+8AVuuukmRowYwZo1a3jssdxKal0NC55yyilcffXV9OnTh09/+tPlvjVpc87ugfyUs0KTWg9Ob56zCwCtApyzq5+zizaaJe0GLIyIzcDtkuYDAyLiTxWJwOpWSVXn/nCwbuyxxx4sWbKE+++/n0WLFjF58mSuvPJKTjvtNCCX9PLfzznnHAAWLFjAihUrdjzHli1b2Lp1KwsWLGDu3Lk7ju+77747vdbChQtZsmQJo0ePBuDll1/mwAMP5MQTT+Tpp5/mrLPO4oQTTmDixImdxrp69WrWrl3LbrvtxtatW9ljjz0q9j5UknN2z+SnnBXq8fSznnIBoNWpZs/ZRRvNEbFd0rXAEcn9V4FXy35lq2slV537w6FudNe7kKY+ffowYcIEJkyYwMiRI5kzZ86OBFxY3Zy/vX37dn7zm98wYEBJu0LvEBGceuqpXHHFFbuce+SRR7jnnnuYNWsWt956K9dff/0u15x99tlccsklrFy5kksuuYRvfOMbPXr9anHO7pkpY4fuksM69jqnwmvVW5mcs6ufs0ud07xQ0snyOieW6FHVuavDrQtPPPEETz755I77y5YtY9iwNxd5mDdv3o7v48aNA3LDg9dcc81OjwE45phjuPbaa3cc7zjUd9RRR3Hbbbexbt06ADZt2sRzzz3Hhg0b2L59OyeffDIzZ85k6dKlu8R59913s27dOj7zmc/w1a9+lTvuuGOnnpMa5JxtZhXX7Dm71DnNZwBfAt6Q9Aq5JYwiInpRDWZmlrN161bOOussNm/eTN++fTn00EOZPXv2jvMvvvgio0aNon///txyyy0AXH311UybNo1Ro0bxxhtv8MEPfpBZs2Zx0UUXMW3aNA4//HD69OnDxRdfzEknnbTjuYYPH87MmTOZOHEi27dvp1+/flx77bXsvvvuTJ06le3btwPs0qvxyiuvMH36dG677TYk8ba3vY1vfOMbnHnmmfzyl7+swrvUK87ZZlZxzZ6zFRFlPUE1tLW1RXt7e9ZhWAf5Icx5Z4zr+qIbTsh99zBkTVq5ciWHHXZY1mF0qqWlhfb2dgYOHJh1KCXr7P2UtCQi2jIKKRP1nrMnf+8hVqzd0qtVgo76811M6vNrDtqzm6Ho/JQ150XrIefsyuppzi6pp1nSwog4qtgxaz7vXPVj/nj1RV1/QLgA0KzqnLPL01lxYKmO+NMC9tBzsOcRXV/kGg+zutRto1nSAOCtwEBJ+5Ib4gPYC+h9VrGGMKn1YN655tfs8eKqrj8g/OFgvfTss89mHULdcc6ujM6KA0v1+L/24VneyQj3IluTaYacXayn+QxgOjAYWMKbCXgL8J0U47I6MGXsUB5f2NcfEGa1wznbzCwl3TaaI+LbwLclnRUR13R3rZmZZcs528wsPSXNaY6IayS9H2gpfExE3JhSXFYP2m9gxGvL+R3Du13XNNWdtcxsF87ZZmaVV2oh4H8AfwEsA7YlhwNwAm5my28D4OG9j+7yktR31jKzXThnlym/k2kvtLz+NCtiWNENUtyZYFZ/St3cpA34QER8ISLOSr6+mGZgVieGjef0cy5j3hnjOv3qzZJN1lwuv/xyRowYwahRo2htbWXx4sUVff7jjz+ezZs39/rx9957L+PGjSO/POe2bds44ogj+PWvf12pENPgnF2O/E6mvbB138O67UiAXGfCncvW9Or5zbLWzDm71M1NHgPeDqwt9YklvYNcr8ZB5Ho4ZkfEtyXtB8wjN2z4LPD3EfFiV89jZo3roYceYv78+SxdupT+/fuzYcMGXnvttYq+xl133VXW44855hiuu+46rrvuOj73uc9xzTXX0NbWxvvf//4KRZiKHuds66CX6ygfBJyefHWlKtt0m6Wg2XN2qT3NA4EVku6R9NP8V5HHvAH8U0QMB/4amCZpOHA+sDAi3gUsTO6bWRNau3YtAwcOpH///gAMHDiQwYMHA7mF8s8991xGjhzJmDFjeOqppwBYv349J598MqNHj2b06NE8+OCDQG6nqqlTpzJy5EhGjRrF7bffvuN5NmzYAMCPfvQjxowZQ2trK2eccQbbtm1j27ZtnHbaaRx++OGMHDmSq666apc4r7rqKq644goef/xxvvOd7/C1r30t9femTD3O2ZLeIWmRpBWSHpd0dnJ8P0n3Snoy+b5vVf4FZlZzmj1nl9rTPKOnTxwRa0l6OSLiJUkrya0TOgmYkFw2B7gPOK+nz28Za78BnnsAho0veumKtVuY/L2HPIev1t19fq+HpLv09pFw3JVdnp44cSKXXnop7373uzn66KOZPHkyH/rQh3ac33vvvVm+fDk33ngj06dPZ/78+Zx99tmcc845jB8/nlWrVnHssceycuVKLrvssh3XQ24710IrV65k3rx5PPjgg/Tr148vfOEL3HTTTYwYMYI1a9bw2GOPAXQ6LDho0CCmT5/OuHHjuPrqq9lvv/0q8e6kaUYvHpPv6FgqaU9giaR7gdPIdXRcKel8ch0dztlmWXPOBqqbs0tdPeO/JA0D3hURCyS9FehT6otIagGOABYDByUNaoAXyI1mdfaYHSNcQ4e6oVVz8kUyRTYuye+s5YJA68wee+zBkiVLuP/++1m0aBGTJ0/myiuv5LTTTgPglFNO2fH9nHPOAWDBggWsWLFix3Ns2bKFrVu3smDBAubOnbvj+L777twhunDhQpYsWcLo0aMBePnllznwwAM58cQTefrppznrrLM44YQTmDhxYqexTps2jfPPP39HbLWsNzm7aTs6Oiv6q8JOpu5MsHrU7Dm71NUz/pFcA3Y/chXZBwOzgKJbskraA7gdmB4RWyTtOBcRISk6e1xEzAZmA7S1tXV6jWVs2Hhom9rtJfmdtTyHrw5007uQpj59+jBhwgQmTJjAyJEjmTNnzo4kV5gv8re3b9/Ob37zGwYM6GLr9i5EBKeeetw6pCoAABbfSURBVCpXXHHFLuceeeQR7rnnHmbNmsWtt97K9ddfv8s1u+22207x1LJycnby+BaapaMjX/RX2EhOeSdTdyZYRThnVz1nlzqneRrwAXK7ShERTwIHFnuQpH7kGsw3RcQdyeE/ShqUnB8ErOtp0GbWGJ544gmefPLJHfeXLVvGsGHDdtyfN2/eju/jxo0DcsOD11xzzU6PgVzxx7XXXrvjeMehvqOOOorbbruNdetyKWfTpk0899xzbNiwge3bt3PyySczc+ZMli5dWuF/ZSZ6lbNh146OwnORK0fvsqMjItoiou2AAw4oJ/bqyxf9FX4V6RAox5SxQ726kNWlZs/Zpc5pfjUiXsu32CX1pYvEmafcxdcBKyPimwWnfgqcClyZfL+zp0GbWWPYunUrZ511Fps3b6Zv374ceuihzJ49e8f5F198kVGjRtG/f39uueUWAK6++mqmTZvGqFGjeOONN/jgBz/IrFmzuOiii5g2bRqHH344ffr04eKLL+akk07a8VzDhw9n5syZTJw4ke3bt9OvXz+uvfZadt99d6ZOncr27dsBOu3VqEM9ztnJdV12dETEWnd0mDW3Zs/Zyq9j1+1F0teBzcBngLOALwArIuIr3TxmPHA/sBzYnhy+kNxw363AUOA5ckvOberu9dva2qK9vb1onFYF+fl/+eHMEpdkmvy9h1ixdsuOnhXP46sNK1eu5LDDDss6jE61tLTQ3t7OwIEDsw6lZJ29n5KWRERbNePoZc4WuTnLmyJiesHxbwAbCwoB94uIc7t7/brK2TeckPvei+XlypWftjbvjHFVf22rT87ZldXTnF1qT/P5wGfJNYDPAO4CftDdAyLiAaCrySQlzauzGlTYYO7BnL/8HD7wPD6zKuhxziY3nePTwHJJy5JjF5IbFbxV0mdJOjpSiTgLPVgFKC0uCDSrH6U2mncHro+I7wNI6pMc+3NagVkN68Wi//mCQPDC/laaZ599NusQ6lmPc3ZTdnSUuApQWlwQaI2kGXJ2qYWAC8kl3LzdgQWVD8fMzCrAObtUJawClBYXBJrVl1J7mgdExNb8nYjYmqz7adYr+SFJ8PzmrEVE3SylVstKqQ+pIudsswblnF0ZvcnZpfY0/6+kv8rfkfQ+4OUev5oZuUZyvmdlxdot3LlsTcYRNa8BAwawcePGWmvw1Z2IYOPGjT1ehzRFztlmDcg5uzJ6m7NL7WmeDvxY0vPk5ry9HZjcsxDNcjy/uXYMGTKE1atXs379+qxDqXsDBgxgyJAhWYeR55xt1oCcsyunNzm71G20fyfpL4H3JIeeiIjXexifmdWYfv36ccghh2QdhlWYc3YJamDljEJeRcNK4ZydrVJ7mgFGAy3JY/5KEhFxYypRmZlZuZyzu5PxyhmFvIqGWX0oqdEs6T+AvwCWAduSwwE4ATeLjpuaVFBhUWAh97iY9Y5zdokyXDmjUH7KmqermdW2Unua24Dh4ZnnzauXm5oUU7jpSSH3uJiVxTnbzKzCSm00P0aukGRtirFYrevFpibFFBYFFnKPi1lZnLPNzCqs1EbzQGCFpN8Cr+YPRsRHU4nKzMzK4ZzdlRSnmlVCV9PVwFPWzLJWaqN5RppBWEbyHx6lyOADxtXkZr02I+sAalZKU80qoavpauApa2a1oNQl5/5L0kHkqrEBfhsR69ILy6qiJ70tVf6AcTW5We85ZxeRwlSzSuhquhp4yppZLSh19Yy/B74B3EduofxrJH05IkrsprSaVeMfHv6gMOs552wzs8ordXrGV4DR+Z4KSQcACwAnYDOz2uOcbWZWYaU2mnfrMLS3EdgthXjMzKx8ztkd1XgBYCm8pr1ZtkptNP9C0j3ALcn9ycBd6YRkVVFjW8h2p7tqcvAHhlknnLM7quECwFJ4TXuz7HXbaJZ0KHBQRHxZ0klAvoX1EHBT2sFZimpoC9nudFdNDv7AMCvknF1EjdZwlMJr2ptlr1hP87eACwAi4g7gDgBJI5NzJ6YanaWrRraQ7U531eTgDwyzDpyzzcxSUmyO20ERsbzjweRYSyoRmZlZb/U6Z0u6XtI6SY8VHJshaY2kZcnX8ZUP2cysPhRrNO/TzbndKxmImZmVrZyc/UPgw50cvyoiWpOv5p4XbWZNrdj0jHZJ/xgR3y88KOlzwJL0wrJU1VERYCm8c6DZDr3O2RHxK0ktKcZmKemsWNr50KzyijWapwM/kfRJ3ky4bcBbgL9LMzBLUZ0UAZbCOwea7SSNnH2mpM8A7cA/RcSLnV0k6XTgdIChQ/17WC2dFUs7H5qlo9tGc0T8EXi/pL8FDk8O/zwifpl6ZJauOigCLIV3DjR7Uwo5+7vAZUAk3/8d+IcuXns2MBugra0tevl61kOdFUs7H5qlo6R1miNiEbAo5VjMzKwCKpWzk0Y4AJK+D8wv9znNzOpVqZubmNU875ZlVlmSBkXE2uTu3wGPdXe9mVkjc6PZGoJ3yzIrj6RbgAnAQEmrgYuBCZJayU3PeBY4I7MAeyq/bXZeHW+f3RsukDarPDearSF4tyyz8kTEKZ0cvq7qgVRK4bbZULfbZ/eGC6TN0uFGs5mZNaY63ja7HC6QNktHsc1NzMzMzMyannuaG1nHOX15TTq3D1wUaGZmZr3jRnMj6zinL68J5/aB5/eZNY0G2/W0HN4t0Kxy3GhudE06py+vsEDQ8/vMmkQD7XpaDu8WaFZZbjSbmVnjaZBdT8vh3QLNKsuN5kbR2fzlJpu7XApvgGJmZma94dUzGkV+/nKhJpq7XIpJrQczfNBeuxxfsXYLdy5bk0FEZmZmVi/c09xImnz+cjHeAMXMzMx6K7VGs6TrgY8A6yLi8OTYfsA8oIXclqx/HxEvphWDmZk1gSbfMrunvKKGWe+kOT3jh8CHOxw7H1gYEe8CFib3zczMeq/j9DRPTetSZ9PUPEXNrDSp9TRHxK8ktXQ4PAmYkNyeA9wHnJdWDA2vsHfFPStl6apAsJB7YsxqmKenlcQrapj1XrULAQ+KiLXJ7ReAg7q6UNLpktolta9fv7460dWbwt4V96z0WlcFgoXcE2NmZtbcMisEjIiQFN2cnw3MBmhra+vyuqbn3pWydVUgWMg9MWZmZs2t2j3Nf5Q0CCD5vq7Kr29mZo0kv2W2mVnKqt1o/ilwanL7VODOKr++mZk1Em+ZbWZVklqjWdItwEPAeyStlvRZ4ErgGElPAkcn960n2m+AG07IfXXczMRStfiZTdy8eFXWYZhZR94y28yqILVGc0ScEhGDIqJfRAyJiOsiYmNEHBUR74qIoyNiU1qv37Bc/JeJSa0HA7gY0BqapOslrZP0WMGx/STdK+nJ5Pu+WcZoZpYVb6Ndj/LFf1N/7t6VKpkydihjD9kv6zDM0vZDvL6+mVmn3Gg2MzMgt74+0HEEcBK5dfVJvn+sqkGZmdWIzJacM3bd+rUU3sQkU6VsgtKRN0WxOlfy+vpVlc+fzokVkc9tzldmXXOjOUu9Sfiex5yZ/LzmnlixdguAP4SsIXS3vr6k04HTAYYOrcLPe2H+dE4sSz63OV+Zdc+N5qx5c5K6UcomKB15UxRrAH+UNCgi1na3vn4mG1I5f1ZEPrc5X5l1z3OazcysO15f38wMN5rNzCzh9fXNzLrm6RlZcAFLU+lN8WApXLBjlRYRp3Rx6qiqBlJMfuvsYeOzjqThdMxXzjNmb3KjOQsuYGkavSkeLIULdqypeevsVHTMV84zZjtzozkrLmBpCr0pHiyFC3as6Xnr7IrrmK+cZ8x25jnNZmZmZmZFuNFsZmZmZlaEp2eYmVl9cBF11fW0kNmFg9bI3Gg2M7P64CLqquppIbMLB63RudFsZmb1w0XUVdPTQmYXDlqj85xmMzMzM7Mi3NNcSfn5dsV4Pp5VwOJnNnHz4lUeCjUzM6sC9zRXUn6+XTGej2dlys81vHPZmowjMTMzaw7uaa40z7ezKpgydqgbzNY8vGqGmdUA9zSbmVlt86oZZlYD3NNsZma1z6N4Zpaxxm00l1qUV0keOrQq627jAW8yYGZmVjmNOz2j1KK8SvLQoVXRpNaDGT5or07PrVi7xXOezczMKqhxe5rBw3nW0LrbeMCbDFhd6zhS6FG8utFx9MsjXtZIGren2czM6lPHkUKP4tWFjqNfHvGyRtPYPc1mTay7+c6lcA+RZcojhXWn4+iXR7ys0TRuo9lDedbE8puf9NaKtVsA3Gg2ACQ9C7wEbAPeiIi2bCMyM6u+xm00H3dl1hGYZaa7+c6lcA+RdeJvI2JD1kGYmWWlcRvNZmZWf9pvgOcegGHjs47EKqC308Q8PcxqkQsBzcysmAD+n6Qlkk7v7AJJp0tql9S+fv363r9SftUMF/7Vve6WxeyOCwitVrmn2cw65aWjrMD4iFgj6UDgXkn/HRG/KrwgImYDswHa2tqirFcbNh7appb1FJa93k4T8/Qwq1XuaTazXXjpKCsUEWuS7+uAnwBjso3IzKz63NNsZrvw0lGWJ+ltwG4R8VJyeyJwacZhmZlVnRvNZmbWnYOAn0iC3GfGzRHxi2xDMjOrPjeazawk5W6W0hXPla5tEfE08N7UXyi/dba3zDZ6lm+cQ6xa3Gg2s6LK3SylK95ExXYobDB75Yym1pN84xxi1eRGs5kVVe5mKV3xXGnbibfONnqWb5xDrJoyWT1D0oclPSHpKUnnZxGDmZmZmVmpqt5oltQHuBY4DhgOnCJpeLXjMDMzMzMrVRbTM8YATyXFJUiaC0wCVmQQi5llLK0Cw0obPngvLj5xRNZhmFkH9ZJDrPoqnbezaDQfDPyh4P5qYGzHi5KtWk8HGDrUE/zNGlFaBYZWh7xihvWCc4hVU80WAlZ0S1Yzq0lpFRhaHTruyqwjsDrkHGLVlEUh4BrgHQX3hyTHzMzMzMxqUhaN5t8B75J0iKS3AJ8AfppBHGZmZmZmJan69IyIeEPSmcA9QB/g+oh4vNpxmJmZmZmVKpM5zRFxF3BXFq9tZmZmZtZTmWxuYmZmZmZWT9xoNjMzMzMrwo1mMzMzM7Mi3Gg2MzMzMytCEbW/b4ik9cBzvXjoQGBDhcPprVqKBRxPMY6ne7UUTy3FArvGMywiDsgqmCw0SM4uxrGmw7Gmw7GWrsucXReN5t6S1B4RbVnHAbUVCzieYhxP92opnlqKBWovnnpST++dY02HY02HY60MT88wMzMzMyvCjWYzMzMzsyIavdE8O+sACtRSLOB4inE83auleGopFqi9eOpJPb13jjUdjjUdjrUCGnpOs5mZmZlZJTR6T7OZmZmZWdncaDYzMzMzK6IuG82SPizpCUlPSTq/k/P9Jc1Lzi+W1FJw7oLk+BOSjs0yHkn7S1okaauk71QiljLjOUbSEknLk+9HZhzPGEnLkq9HJP1dlvEUnB+a/J/9c1axSGqR9HLB+zOr3FjKiSc5N0rSQ5IeT36GBmQVj6RPFrw3yyRtl9SaYTz9JM1J3peVki4oN5Z6U2t5O41Y08qhacRacL5i+SzNWNPIL2nEmsXvegmxflDSUklvSPp4h3OnSnoy+Tq1VmOV1Frw//+opMlpx9qpiKirL6AP8HvgncBbgEeA4R2u+QIwK7n9CWBecnt4cn1/4JDkefpkGM/bgPHA54Hv1MD7cwQwOLl9OLAm43jeCvRNbg8C1uXvZxFPwfnbgB8D/5zhe9MCPFaJn5kKxdMXeBR4b3J//yx/tzpcMxL4fcbvzxRgbsHP9bNASyX//2r5q8z3ruJ5O8VYK55D04q14HxF8lnK72vF80uKsVb1d73EWFuAUcCNwMcLju8HPJ183ze5vW+Nxvpu4F3J7cHAWmCfNH9mO/uqx57mMcBTEfF0RLwGzAUmdbhmEjAnuX0bcJQkJcfnRsSrEfEM8FTyfJnEExH/GxEPAK+UGUOl4nk4Ip5Pjj8O7C6pf4bx/Dki3kiODwAqUbVazs8Pkj4GPEPu/ck0lhSUE89E4NGIeAQgIjZGxLYM4yl0SvLYcpUTTwBvk9QX2B14DdhSgZjqRa3l7VRiTSmHphIrVDyfpRlrGvklrVir/bteNNaIeDYiHgW2d3jsscC9EbEpIl4E7gU+XIuxRsT/RMSTye3nyXWiVX2n1XpsNB8M/KHg/urkWKfXJI2uP5H7y7SUx1YznjRUKp6TgaUR8WqW8UgaK+lxYDnw+YJGdNXjkbQHcB5wSZkxlB1Lcu4QSQ9L+i9Jf5NxPO8GQtI9ydDauRnHU2gycEvG8dwG/C+53pFVwL9FxKYKxFQvai1vd6fWcmh3aimfFVNr+SWtWKv9u17O70ct/m4VJWkMuZ7q31corpL1rfYLWu2TNAL4Grm/7jMVEYuBEZIOA+ZIujsiKtkz3xMzgKsiYmt6nb0lWwsMjYiNkt4H/KekERGRVe9lX3JTjUYDfwYWSloSEQszigfI/dEF/DkiHssyDnI9LNvIDSvuC9wvaUFEPJ1tWJaGWsqh3ZhB7eSzYmoyv3TBv+spkjQI+A/g1Ijo2HOeunrsaV4DvKPg/pDkWKfXJEMkewMbS3xsNeNJQ1nxSBoC/AT4TERU4q+4irw/EbES2EpunmBW8YwFvi7pWWA6cKGkM7OIJRmq3ggQEUvI/cX97jJiKSsecj0Gv4qIDRHxZ+Au4K8yjCfvE1Sml7nceKYAv4iI1yNiHfAg0FahuOpBreXt7tRaDk0r1krnszRjTSO/pBVrtX/Xy/n9qMXfrS5J2gv4OfCViPhNhWMrTVqTpdP6IvcX59PkCkLyE8lHdLhmGjtP0L81uT2CnQtKnqb8YqVex1Nw/jQqVwhYzvuzT3L9STXy/3UIbxYCDgOeBwZm/f+VHJ9B+YWA5bw3B+R/dskVVawB9sswnn2BpSTFm8AC4IQs/6/IdQqsAd5ZAz/L5wE3JLffBqwARlUirnr4KvO9q3jeTjHWiufQtGLtcM0M0i8ErKn8kmKsVf1dLyXWgmt/yK6FgM8k7+++ye2yPkdSjPUtwEJgepo/p0X/DVm+eBlv/PHA/5DrXftKcuxS4KPJ7QHkqoGfAn5LwYcm8JXkcU8Ax9VAPM8Cm8j1oq6mQyVpNeMBLiI3F2tZwdeBGcbzaXIFKsvIJcyPZf3/VfAcM6jAh0wZ783JHd6bE7N+b4BPJTE9Bny9BuKZAPymEnFU4P9rj+T44+Q+RL9cybjq4avM/8uK5+2U/p9TyaFpva8FzzGDlBvNFfgZqHh+SelnoOq/6yXEOppc++J/yfWGP17w2H9I/g1PAVNrNdbk///1Dr9brWnH2/HL22ibmZmZmRVRj3OazczMzMyqyo1mMzMzM7Mi3Gg2MzMzMyvCjWYzMzMzsyLcaDYzMzMzK8KNZms4kvaXtCz5ekHSmuT2Vkn/N+v4zMzsTc7ZVi+85Jw1NEkzgK0R8W9Zx2JmZt1zzrZa5p5maxqSJkian9yeIWmOpPslPSfpJElfl7Rc0i8k9Uuue5+k/5K0RNI9yb73ZmaWMudsqzVuNFsz+wvgSOCjwI+ARRExEngZOCFJwteQ28rzfcD1wOVZBWtm1uScsy1TfbMOwCxDd0fE65KWA32AXyTHlwMtwHuAw4F7JZFcszaDOM3MzDnbMuZGszWzVwEiYruk1+PNCf7byf1uiNy+9+OyCtDMzHZwzrZMeXqGWdeeAA6QNA5AUj9JIzKOyczMOuecbalyo9msCxHxGvBx4GuSHgGWAe/PNiozM+uMc7alzUvOmZmZmZkV4Z5mMzMzM7Mi3Gg2MzMzMyvCjWYzMzMzsyLcaDYzMzMzK8KNZjMzMzOzItxoNjMzMzMrwo1mMzMzM7Mi/j/jX4C8rd8uDgAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "%matplotlib inline\n", "import matplotlib.pyplot as plt\n", "\n", "true_rate = 2.3\n", "observations = [Model1()({\"rate\": true_rate}), Model2()({\"rate\": 30})]\n", "fig, axes = plt.subplots(ncols=2)\n", "fig.set_size_inches((12, 4))\n", "for ax, title, obs in zip(axes, [\"Observation\", \"Competition\"], observations):\n", " ax.step(obs[\"t\"], obs[\"X\"])\n", " ax.legend([\"Species X\", \"Species Y\"])\n", " ax.set_xlabel(\"Time\")\n", " ax.set_ylabel(\"Concentration\")\n", " ax.set_title(title);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We observe that species $X$ is converted into species $Y$ in both cases.\n", "The difference of the concentrations over time can be quite subtle.\n", "\n", "We define a distance function as $L_1$ norm of two trajectories, evaluated at 20 time points:\n", "\n", "$$\n", " \\mathrm{distance}(X_1, X_2) =\n", " \\frac{1}{N}\\sum_{n=1}^{N} \n", " \\left |X_1(t_n) -X_2(t_n) \n", " \\right|, \\quad t_n = \\frac{n}{N}T, \\quad N=20 \\,.\n", "$$" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that we only consider the concentration of species $Y$ for distance calculation. And in code:" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "N_TEST_TIMES = 20\n", "\n", "t_test_times = np.linspace(0, MAX_T, N_TEST_TIMES)\n", "\n", "\n", "def distance(x, y):\n", " xt_ind = np.searchsorted(x[\"t\"], t_test_times) - 1\n", " yt_ind = np.searchsorted(y[\"t\"], t_test_times) - 1\n", " error = (\n", " np.absolute(x[\"X\"][:, 1][xt_ind] - y[\"X\"][:, 1][yt_ind]).sum()\n", " / t_test_times.size\n", " )\n", " return error" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For ABC, we choose for both models a uniform prior over the interval $[0, 100]$ for their single rate parameters:" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "from pyabc import RV, Distribution\n", "\n", "prior = Distribution(rate=RV(\"uniform\", 0, 100))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We initialize the ABCSMC class passing the two models, their priors and the distance function." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "INFO:Sampler:Parallelizing the sampling on 4 cores.\n" ] } ], "source": [ "from pyabc import ABCSMC\n", "from pyabc.populationstrategy import AdaptivePopulationSize\n", "\n", "abc = ABCSMC(\n", " [Model1(), Model2()],\n", " [prior, prior],\n", " distance,\n", " population_size=AdaptivePopulationSize(500, 0.15),\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We initialize a new ABC run, taking as observed data the one generated by model 1.\n", "The ABC run is to be stored in the sqlite database located at ``/tmp/mjp.db``." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "INFO:History:Start \n" ] } ], "source": [ "abc_id = abc.new(\"sqlite:////tmp/mjp.db\", observations[0])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We start pyABC which automatically parallelizes across all available cores." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "INFO:ABC:Calibration sample before t=0.\n", "INFO:Epsilon:initial epsilon is 10.65\n", "INFO:ABC:t: 0, eps: 10.65.\n", "INFO:ABC:Acceptance rate: 500 / 1023 = 4.8876e-01, ESS=5.0000e+02.\n", "INFO:Adaptation:Change nr particles 500 -> 115\n", "INFO:ABC:t: 1, eps: 6.75.\n", "INFO:ABC:Acceptance rate: 115 / 272 = 4.2279e-01, ESS=1.0816e+02.\n", "INFO:Adaptation:Change nr particles 115 -> 89\n", "INFO:ABC:t: 2, eps: 5.316484503773134.\n", "INFO:ABC:Acceptance rate: 89 / 167 = 5.3293e-01, ESS=6.2536e+01.\n", "INFO:Adaptation:Change nr particles 89 -> 82\n", "INFO:ABC:t: 3, eps: 3.9215865820412024.\n", "INFO:ABC:Acceptance rate: 82 / 220 = 3.7273e-01, ESS=5.8737e+01.\n", "INFO:Adaptation:Change nr particles 82 -> 92\n", "INFO:ABC:t: 4, eps: 3.2.\n", "INFO:ABC:Acceptance rate: 92 / 330 = 2.7879e-01, ESS=2.1439e+01.\n", "INFO:Adaptation:Change nr particles 92 -> 75\n", "INFO:ABC:t: 5, eps: 2.45.\n", "INFO:ABC:Acceptance rate: 75 / 401 = 1.8703e-01, ESS=2.1633e+01.\n", "INFO:Adaptation:Change nr particles 75 -> 96\n", "INFO:ABC:t: 6, eps: 2.145736215699017.\n", "INFO:ABC:Acceptance rate: 96 / 549 = 1.7486e-01, ESS=2.9474e+01.\n", "INFO:Adaptation:Change nr particles 96 -> 102\n", "INFO:ABC:t: 7, eps: 1.75.\n", "INFO:ABC:Acceptance rate: 102 / 792 = 1.2879e-01, ESS=7.8885e+01.\n", "INFO:Adaptation:Change nr particles 102 -> 58\n", "INFO:ABC:t: 8, eps: 1.4.\n", "INFO:ABC:Acceptance rate: 58 / 357 = 1.6246e-01, ESS=4.9798e+01.\n", "INFO:Adaptation:Change nr particles 58 -> 57\n", "INFO:ABC:t: 9, eps: 1.25.\n", "INFO:ABC:Acceptance rate: 57 / 583 = 9.7770e-02, ESS=4.1186e+01.\n", "INFO:Adaptation:Change nr particles 57 -> 61\n", "INFO:ABC:t: 10, eps: 1.0627609318694404.\n", "INFO:ABC:Acceptance rate: 61 / 1185 = 5.1477e-02, ESS=5.7597e+01.\n", "INFO:Adaptation:Change nr particles 61 -> 46\n", "INFO:ABC:t: 11, eps: 0.95.\n", "INFO:ABC:Acceptance rate: 46 / 1461 = 3.1485e-02, ESS=4.2846e+01.\n", "INFO:Adaptation:Change nr particles 46 -> 52\n", "INFO:ABC:t: 12, eps: 0.8.\n", "INFO:ABC:Acceptance rate: 52 / 4317 = 1.2045e-02, ESS=3.9330e+01.\n", "INFO:Adaptation:Change nr particles 52 -> 53\n", "INFO:ABC:t: 13, eps: 0.75.\n", "INFO:ABC:Acceptance rate: 53 / 4867 = 1.0890e-02, ESS=4.8614e+01.\n", "INFO:Adaptation:Change nr particles 53 -> 51\n", "INFO:ABC:t: 14, eps: 0.6526353965182403.\n", "INFO:ABC:Acceptance rate: 51 / 15746 = 3.2389e-03, ESS=4.4594e+01.\n", "INFO:Adaptation:Change nr particles 51 -> 55\n", "INFO:ABC:Stopping: minimum epsilon.\n", "INFO:History:Done \n" ] } ], "source": [ "history = abc.run(minimum_epsilon=0.7, max_nr_populations=15)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We first inspect the model probabilities." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEvCAYAAABIeMa5AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAbt0lEQVR4nO3deZhddZ3n8fc3C5tASEgJSAiJWToGMWJiAk36ARt0ArQgW8uijSzyNCOK3dqYmWYgD9jTwW5kGgcakRbQbllHmDQgiw1o48iSICABghECJBAIW9iXUN/545zSa1GVqtStcyuV8349Tz11z3K/91s3qfrcs/1OZCaSpPoaMtANSJIGlkEgSTVnEEhSzRkEklRzBoEk1ZxBIEk1N2ygG5A6LFq06P3Dhg27EPgwfkhptXbggTVr1hw3ffr0Zwe6GbWWQaD1xrBhwy7cdtttP9TW1vbikCFDvMClhdrb22PVqlVTV65ceSGw/0D3o9byU5fWJx9ua2t72RBovSFDhmRbW9tqiq0x1YxBoPXJEENg4JTvvX8Tash/dKmTiJh+wAEHjO+Yfueddxg5cuS0T3ziExPXpc7222+/89NPP73W3a+9WUeqmkEgdbLpppu2L1myZNNXX301AK6++uott9lmm3cGui+pKgaB1IW999579ZVXXrkVwKWXXjrq4IMPfqFj2TPPPDN07733njB58uSp06ZNm3LnnXduCrBy5cqhu++++6SJEyfu9NnPfnbHxgEdzzvvvFE777zzh6ZMmTL1iCOO2HHNmjUt/5mk7hgEUhc+//nPv3D55ZePfP311+Ohhx7abLfddnutY9nJJ5/8gWnTpr3+yCOPPHjGGWesOOqoo8YDzJ079wO77bbbq0uXLl184IEHvvT0009vBHDPPfdsctVVV41auHDhww8//PCDQ4YMyfPPP3/rgfrZpM4MAqkLs2bNemP58uUbf+973xu19957r25cdtddd21x7LHHPg+w//77v/LSSy8Ne+GFF4bccccdWxxzzDHPAxx22GGrt9xyy3cBbrjhhi0eeOCBzaZNm/ahKVOmTL399tu3fPTRRzdu/U8ldc2DVFI35syZ89Jpp522w0033bTk2Wef7fPvSmbGoYce+vy55567oj/7k/qLWwRSN0444YTnvv71rz81c+bMNxrnz5o165WLLrpoa4Brr712i5EjR64ZNWpU+6677vrKxRdfvDXAFVdcseXLL788FGDOnDkvX3vttSNXrFgxDIpjDI888shGrf55pO64RSB1Y8KECe+ccsop7xlu4cwzz3zqyCOPHDd58uSpm266afvFF1/8GMD8+fOfOvjggz84ceLEnWbMmPHqdttt9zbA9OnT3zzllFNW7LXXXpPb29sZPnx4nnPOOU9Mnjz57Vb/TFJXwltVan1x3333LZs2bdpzA91Hnd13332jp02bNm6g+1BruWtIkmrOIJCkmjMIJKnmDAJJqjmDQJJqziCQpJozCKQGhx566LhRo0ZNmzRp0k4D3Ut/Wrp06fBZs2ZNnjBhwk4TJ07c6Ywzznj/QPek9YcXlGm9NW7uddP7s96y+fst6mmdY4455rmTTjrp2aOPPnp8T+v22bwR/fpzMW91jz/X8OHDOeuss5bPnj379RdffHHILrvsMnXfffd9efr06W/2ay8alNwikBrss88+r7a1tW1wY0TvuOOO78yePft1gJEjR7ZPmDDhjSeeeMJhLgQYBFLtLFmyZKMHH3xwsz322OPVge5F6weDQKqR1atXDznooIMmzJ8//8lRo0a1D3Q/Wj8YBFJNvPXWW7HffvtNOPTQQ1846qijXhrofrT+MAikGmhvb+ewww7bcfLkyW/OmzfvmYHuR+sXg0Bq8OlPf3r87Nmzpzz22GMbb7PNNh85++yzRw90T/3h5ptv3vyaa67Z+vbbb99iypQpU6dMmTL18ssvHzHQfWn94DDUWm84DPXAcxjqenKLQJJqziCQpJozCCSp5gwCrU/a29vbY6CbqKvyvffaghoyCLQ+eWDVqlUjDIPWa29vj1WrVo0AHhjoXtR6Djqn9caaNWuOW7ly5YUrV678MH5IabV24IE1a9YcN9CNqPU8fVSSas5PXZJUcwaBJNWcQSBJNTfoDhaPHj06x40bN9BtSNKgsmjRoucys62rZYMuCMaNG8fChQsHug1JGlQi4vHulrlrSJJqziCQpJozCCSp5gwCSao5g0CSaq6yIIiI70fEsxHR5SBWUTgnIpZGxP0R8bGqepEkda/KLYKLgTlrWb4PMKn8Oh745wp7kSR1o7IgyMyfAy+sZZUDgB9k4Q5gq4jYrqp+JEldG8gLyrYHnmyYXl7Oe7rzihFxPMVWA2PHjm1Jc9LajJt73XvmLZu/XyV1q6w92Or2R+0N5T3ur9owSA4WZ+YFmTkjM2e0tXV5hbQkqY8GMghWADs0TI8p50mSWmggg2AB8Bfl2UO7Aqsz8z27hSRJ1arsGEFEXArsCYyOiOXAacBwgMw8H7ge2BdYCrwOHF1VL5Kk7lUWBJl5eA/LE/hSVa8vSeqdQXGwWJJUHYNAkmrOIJCkmjMIJKnmDAJJqjmDQJJqziCQpJozCCSp5gwCSao5g0CSam4g70cwuM0b0cW81a3vQ5Ka5BaBJNWcQSBJNWcQSFLNeYxgfeOxB0kt5haBJNWcQSBJNWcQSFLNGQSSVHMGgSTVnEEgSTXn6aN10dVpqeCpqZLcIpCkujMIJKnmDAJJqjmDQJJqziCQpJozCCSp5gwCSao5g0CSam7DvqDMi6gkqUduEUhSzVUaBBExJyKWRMTSiJjbxfKxEXFrRPwqIu6PiH2r7EeS9F6VBUFEDAXOBfYBpgKHR8TUTqudAlyRmbsAhwHnVdWPJKlrVW4RzASWZuajmfk2cBlwQKd1EtiyfDwCeKrCfiRJXagyCLYHnmyYXl7OazQP+FxELAeuB77cVaGIOD4iFkbEwlWrVlXRqyTV1kAfLD4cuDgzxwD7Aj+MiPf0lJkXZOaMzJzR1tbW8iYlaUNWZRCsAHZomB5Tzmt0LHAFQGb+EtgEGF1hT5KkTqoMgruBSRExPiI2ojgYvKDTOk8AewFExIcogsB9P5LUQpVdUJaZayLiROBGYCjw/cxcHBGnAwszcwHwNeB7EfFXFAeOv5CZWVVPqkhXF+550Z40aFR6ZXFmXk9xELhx3qkNjx8Edq+yB0nS2g30wWJJ0gAzCCSp5gwCSao5g0CSas4gkKSaMwgkqeYMAkmqOYNAkmrOIJCkmjMIJKnmDAJJqjmDQJJqziCQpJozCCSp5gwCSao5g0CSas4gkKSaMwgkqeYMAkmqOYNAkmrOIJCkmjMIJKnmhg10A1KVxs297j3zls3fbwA6kdZfbhFIUs0ZBJJUcwaBJNWcQSBJNderIIiIT0eEoSFJG6De/nH/LPCbiPhWREypsiFJUmv1Kggy83PALsBvgYsj4pcRcXxEbFFpd5KkyvV6d09mvgxcBVwGbAccCNwTEV+uqDdJUgv09hjBARFxNXAbMByYmZn7ANOAr1XXniSpar3dIjgIODszd87Mf8jMZwEy83Xg2O6eFBFzImJJRCyNiLndrPPnEfFgRCyOiB+t808gSWpKb4NgZWb+vHFGRJwJkJn/0dUTImIocC6wDzAVODwipnZaZxLw34DdM3Mn4Kvr1r4kqVm9DYJPdjFvnx6eMxNYmpmPZubbFMcWDui0zheBczPzRYCOLQ1JUuusNQgi4oSI+DUwJSLub/h6DLi/h9rbA082TC8v5zWaDEyOiF9ExB0RMWddfwBJUnN6Gn30R8BPgL8HGvfxv5KZL/TT608C9gTGAD+PiJ0z86XGlSLieOB4gLFjx/bDy0qSOvS0aygzcxnwJeCVhi8iYlQPz10B7NAwPaac12g5sCAz38nMx4BHKIKhcxMXZOaMzJzR1tbWw8tKktZFT0HQcRbPImBh+X1Rw/Ta3A1MiojxEbERcBiwoNM611BsDRARoyl2FT3a2+YlSc1b666hzPyz8vv4dS2cmWsi4kTgRmAo8P3MXBwRpwMLM3NBuexTEfEg8C7wN5n5/Lq+liSp79YaBBHxsbUtz8x7elh+PXB9p3mnNjxO4K/LL0nSAOjpYPFZa1mWwJ/2Yy+SpAHQ066hT7SqEUnSwOhp19CfZuYtEXFQV8sz88fVtCVJapWedg3tAdwCfLqLZQkYBJI0yPW0a+i08vvRrWlHktRqvR2GeuuIOCci7omIRRHxTxGxddXNSZKq19tB5y4DVgEHA4eUjy+vqilJUuv0dIygw3aZeUbD9Dcj4rNVNCRJaq3ebhHcFBGHRcSQ8uvPKa4KliQNcj2dPvoKxdlBQXHTmH8tFw0BXgW+Xml3kqTK9XTW0BatakSSNDB6e4yAiBhJMUT0Jh3zOt++UpI0+PQqCCLiOOAkinsK3AvsCvwSxxpSPxg397ou5y+bv1+LO5HqqbcHi08CPg48Xo4/tAvw0tqfIkkaDHobBG9m5psAEbFxZj4M/FF1bUmSWqW3xwiWR8RWFHcUuzkiXgQer64tSVKr9CoIMvPA8uG8iLgVGAHcUFlXkqSWWZezhj4GzKa4ruAXmfl2ZV1Jklqmt4POnQpcAmwNjAYuiohTqmxMktQavd0iOBKY1nDAeD7FaaTfrKoxiXkjupi3uvV9SBu43p419BQNF5IBGwMr+r8dSVKr9TTW0HcojgmsBhZHxM3l9CeBu6pvT5JUtZ52DS0svy8Crm6Yf1sl3UiSWq6nQecu6XgcERsBk8vJJZn5TpWNSZJao7djDe1JcdbQMoohqXeIiKMcdE6SBr/enjV0FvCpzFwCEBGTgUuB6VU1Jklqjd6eNTS8IwQAMvMRYHg1LUmSWqm3WwSLIuJCfn+HsiP5/YFkSdIg1tsg+EvgS8BXyun/BM6rpCNJUkv1GAQRMRS4LzOnAN+uviVJUiv1eIwgM98FlkTE2Bb0I0lqsd7uGhpJcWXxXcBrHTMzc/9KupIktUxvg+B/VNqFJGnArHXXUERsEhFfBQ4FplDch+BnHV89FY+IORGxJCKWRsTctax3cERkRMxY559AktSUno4RXALMAH4N7ENxYVmvlAeZzy2fNxU4PCKmdrHeFsBJwJ29rS1J6j89BcHUzPxcZn4XOAT4k3WoPRNYmpmPlnczuww4oIv1zgDOBN5ch9qSpH7SUxD8bmC5zFyzjrW3B55smF5ezvud8vaXO2TmdWsrFBHHR8TCiFi4atWqdWxDkrQ2PR0snhYRL5ePA9i0nA4gM3PLvr5wRAyhuC7hCz2tm5kXABcAzJgxI/v6mpKk9+ppGOqhTdReAezQMD2GP7yr2RbAh4HbIgJgW2BBROyfmQ5fIUkt0ttB5/ribmBSRIwv72VwGLCgY2Fmrs7M0Zk5LjPHAXcAhoAktVhlQVAeUzgRuBF4CLgiMxdHxOkR4YVokrSe6O0FZX2SmdcD13ead2o36+5ZZS+SpK5VuWtIkjQIGASSVHMGgSTVnEEgSTVnEEhSzRkEklRzBoEk1ZxBIEk1V+kFZdJ6ad6Ibuavbm0f0nrCLQJJqjmDQJJqziCQpJozCCSp5gwCSao5g0CSas4gkKSaMwgkqeYMAkmqOYNAkmrOIJCkmjMIJKnmDAJJqjmDQJJqziCQpJozCCSp5gwCSao5g0CSas4gkKSaMwgkqeYMAkmqOYNAkmrOIJCkmhtWZfGImAP8EzAUuDAz53da/tfAccAaYBVwTGY+XmVP6rtxc6/rcv6yTVrciKR+VdkWQUQMBc4F9gGmAodHxNROq/0KmJGZHwGuAr5VVT+SpK5VuWtoJrA0Mx/NzLeBy4ADGlfIzFsz8/Vy8g5gTIX9SJK6UGUQbA882TC9vJzXnWOBn1TYjySpC5UeI+itiPgcMAPYo5vlxwPHA4wdO7aFnUnShq/KLYIVwA4N02PKeX8gIvYG/hbYPzPf6qpQZl6QmTMyc0ZbW1slzUpSXVUZBHcDkyJifERsBBwGLGhcISJ2Ab5LEQLPVtiLJKkblQVBZq4BTgRuBB4CrsjMxRFxekTsX672D8DmwJURcW9ELOimnCSpIpUeI8jM64HrO807teHx3lW+viSpZ+vFweL+0NXFTl7oJEk9c4gJSao5g0CSas4gkKSaMwgkqeYMAkmqOYNAkmrOIJCkmjMIJKnmNpgLyvR7XlwnaV24RSBJNWcQSFLNGQSSVHMGgSTVnEEgSTVnEEhSzXn66ADp6hRP8DRPSa3nFoEk1ZxBIEk1ZxBIUs15jKAH7suXtKFzi0CSas4gkKSaMwgkqeYMAkmqOYNAkmrOIJCkmjMIJKnmDAJJqjmDQJJqziCQpJozCCSp5gwCSaq5SoMgIuZExJKIWBoRc7tYvnFEXF4uvzMixlXZjyTpvSoLgogYCpwL7ANMBQ6PiKmdVjsWeDEzJwJnA2dW1Y8kqWtVbhHMBJZm5qOZ+TZwGXBAp3UOAC4pH18F7BURUWFPkqROIjOrKRxxCDAnM48rpz8PzMrMExvWeaBcZ3k5/dtynec61ToeOL6c/CNgSS/bGA081+NafVNV7cFWt8ra1q2+9mCrW2XtwVZ3XWvvmJltXS0YFDemycwLgAvW9XkRsTAzZ1TQUmW1B1vdKmtbt/rag61ulbUHW93+rF3lrqEVwA4N02PKeV2uExHDgBHA8xX2JEnqpMoguBuYFBHjI2Ij4DBgQad1FgBHlY8PAW7JqvZVSZK6VNmuocxcExEnAjcCQ4HvZ+biiDgdWJiZC4B/AX4YEUuBFyjCoj+t8+6k9aD2YKtbZW3rVl97sNWtsvZgq9tvtSs7WCxJGhy8sliSas4gkKSaMwgkqeYGxXUEvRURUyiuVt6+nLUCWJCZDw1cV2tX9rw9cGdmvtowf05m3tBE3ZlAZubd5dAec4CHM/P6ppv+w9f5QWb+RX/WLOvOprg6/YHMvKmJOrOAhzLz5YjYFJgLfAx4EPifmbm6j3W/AlydmU/2tbdu6nacYfdUZv40Io4A/hh4CLggM99psv4HgYMoTtt+F3gE+FFmvtxc5xrMNpiDxRHxDeBwiqEslpezx1D8Ul2WmfMret2jM/OiPj73K8CXKH7JPwqclJn/t1x2T2Z+rI91T6MY42kYcDMwC7gV+CRwY2b+XR/rdj79N4BPALcAZOb+falb1r4rM2eWj79I8b5cDXwK+Pe+/vtFxGJgWnkW2wXA65TDmZTzD+pj3dXAa8BvgUuBKzNzVV9qdar7bxT/bpsBLwGbAz8u+43MPGotT++p9leAPwN+DuwL/Kp8jQOB/5qZtzXVvFouIt6fmc82XSgzN4gvik82w7uYvxHwmwpf94kmnvtrYPPy8ThgIUUYAPyqybpDKf6YvAxsWc7fFLi/ibr3AP8K7AnsUX5/uny8R5Pv468aHt8NtJWP3wf8uom6DzX232nZvc30S7Fr9VMUp0GvAm6guC5miybq3l9+HwY8Awwtp6OZf7vG/xfl482A28rHY5v5/1bWGAHMBx6mOBX8eYoPOPOBrZqpvZbX/EkTz90S+Hvgh8ARnZad12Rf2wL/TDHo5tbAvPK9vwLYrom6ozp9bQ0sA0YCo5rpeUPaNdQOfAB4vNP87cplfRYR93e3CNimidJDstwdlJnLImJP4KqI2LGs3VdrMvNd4PWI+G2Wm/2Z+UZENPNezABOAv4W+JvMvDci3sjMnzVRs8OQiBhJ8cc1svx0nZmvRcSaJuo+0LDVdl9EzMjMhRExGWhmN0tmZjtwE3BTRAyn2Ao7HPhHoMsxXXphSLl76H0Uf6xHUPxh3RgY3kS/HYZR7BLamGJrg8x8ouy/GVdQbBnumZkrASJiW4pgvIIiMNdZRHS3VRwUW9F9dRHwG+D/AMdExMEUgfAWsGsTdQEuBq6j+De8Ffg3ii2wzwDn897BN3vrOd779217ig9oCXywj3U3qC2COcBS4CcUF1lcQPEJbSnFwHbN1H6G4j/djp2+xlHsy+1r3VuAj3aaNwz4AfBuE3XvBDYrHw9pmD+CTp+K+1h/DHAl8L9pYouoU81lwKPAY+X37cr5m9PcJ/cRFL+Yvy3fl3fK+j+j2DXU17rdfoLueO/7WPevyv4eB74C/AfwPYpPlKc1+R6fBNxf1nsYOLqc3wb8vMnaS/qyrBd13y1/T27t4uuNJure22n6b4FfUHzKbup3hD/cun1iba+7jnW/Vv5N27lh3mPN9Pq7Ov1RZH35ovg0uStwcPm1K+WmcJN1/wWY3c2yHzVRdwywbTfLdm+i7sbdzB/d+J+oH96X/SgOuFb5b7oZML4f6mwJTAOmA9v0Q73JFf7MHwA+UD7eimL4lZn9VHunst6Ufu75JuDkxveWYmv5G8BPm6j7ADCpm2VPNlH3IRo+JJXzvgAsBh5v8r24r+HxNzst6/NuzvL5HR/Cvg1sATzaH/9+G8zBYkkDp9ytN5dit8f7y9nPUIwnNj8zX+xj3UMo/ni+Z+j5iPhMZl7Tx7rfAm7KzJ92mj8H+E5mTupL3bLG6cC3suEswHL+RIr34pC+1m6otT/w34Fxmblt0/UMAklVaubMug2pbn/XLk+HnpCZDzRb1yCQVKmIeCIzx9a9bpW1m627IZ01JGmAVHVm3WCrW2XtKns2CCT1h22A/wJ0PhYQwP+rUd0qa1fWs0EgqT9cS3Fx5L2dF0TEbTWqW2Xtynr2GIEk1Zyjj0pSzRkEklRzBoE2eBGxTUT8KCIejYhFEfHLiDhwgHrZMyL+uGH6LyOi34fxltaFB4u1QYuIAK4BLsnMI8p5OwJ9HjK7F685LDO7GyhvT+BVyrM8MvP8qvqQesuDxdqgRcRewKmZuUcXy4ZSDJO8J8VonOdm5nfLUWDnUYz2+GFgEfC5zMyImE4xzsvm5fIvZObT5Vkb9wKzKe5P8AhwCsUw6M8DR1IMA34HxUBqq4AvU9xn4NXM/MeI+CjF6JSbUQySd0xmvljWvpPi3g9bAcdm5n/237ukunPXkDZ0O1EM09uVY4HVmflx4OPAFyNifLlsF+CrwFSK4X13L4dq/g5wSGZOB74PNN7kZ6PMnJGZZwG3A7tm5i4UN0s6OTOXUfyhPzszP9rFH/MfAN/IzI9QjjbasGxYFjfu+Wqn+VLT3DWkWomIcyk+tb9NMdTzR8qBzaAYsnpSueyuzFxePudeiiHHX6LYQri52OPEUIob83S4vOHxGODyiNiOYqvgsR76GkFxA5eOeztcQjHKZIcfl98Xlb1I/cYg0IZuMcWQ5ABk5pciYjTF3eCeAL6cmTc2PqHcNfRWw6x3KX5XAlicmbt181qvNTz+DvDtzFzQsKupGR39dPQi9Rt3DWlDdwuwSUSc0DBvs/L7jcAJHXfniojJEfG+tdRaArRFxG7l+sMjYqdu1h0BrCgfN95n+BWKceT/QGauBl6MiD8pZ32e4uY5UuX8ZKENWnmA9zPA2RFxMsVB2tcobphyJcVulnvKs4tWUdxOsLtab5e7kc4pd+UMA/4XxVZHZ/OAKyPiRYow6jj28O8UtyM9gOJgcaOjgPMjYjOKu5Qdve4/sbTuPGtIkmrOXUOSVHMGgSTVnEEgSTVnEEhSzRkEklRzBoEk1ZxBIEk1ZxBIUs39fzIt6F6q3DMjAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "ax = history.get_model_probabilities().plot.bar()\n", "ax.set_ylabel(\"Probability\")\n", "ax.set_xlabel(\"Generation\")\n", "ax.legend(\n", " [1, 2], title=\"Model\", ncol=2, loc=\"lower center\", bbox_to_anchor=(0.5, 1)\n", ");" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The mass at model 2 decreased, the mass at model 1 increased slowly.\n", "The correct model 1 is detected towards the later generations.\n", "We then inspect the distribution of the rate parameters:" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAGoCAYAAAATsnHAAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nOzdeXxU5dXA8d+ZPZOEhCwQtoDshE0Fd9+KKH3RKlrFhbpRFa3WpdXa2rq1llpbq69WqUsVF7RVW1sFtbii4oICVpF9h4AEspCQZbY787x/3EkIIXtCMoHz/XzGzL33uTdnJjJnnnuf+xwxxqCUUkolGkdnB6CUUkrVRxOUUkqphKQJSimlVELSBKWUUiohaYJSSimVkDRBKaWUSkiaoJSqh4gMEBEjIq5mtJ0uIh93RFxKHUo0QakuT0Q2i0hYRLLqrP9vPMkM6JzIauJ4QkTWiEhMRKZ3ZixKdSWaoNTBYhMwrXpBREYD/s4LZx9fA9cCX3Z2IEp1JZqg1MFiDnBpreXLgOdqNxCRNBF5TkQKRWSLiNwuIo74NqeI/ElEikRkI/C9evZ9SkR2iMh2EZkpIs7mBGaMmWWMeQ8ItukVKnWI0QSlDhaLgG4iMiKeOC4Enq/T5mEgDRgInISd0H4Y3zYDOAM4AhgPTK2z7zOABQyOt/kucGW7vwqlVA1NUOpgUt2LmgSsArZXb6iVtH5pjCk3xmwG7gcuiTc5H3jQGJNvjCkBfl9r357A6cBPjDGVxphdwP/Fj6eUOkCaHKGkVBcyB/gIOIw6p/eALMANbKm1bgvQJ/68N5BfZ1u1/vF9d4hI9TpHnfZKqXamCUodNIwxW0RkE3Zv54o6m4uACHayWRlfl8veXtYOoF+t9rm1nucDISDLGGO1d9xKqfrpKT51sLkCmGiMqay90hgTBV4GficiqSLSH7iJvdepXgZuEJG+ItIduLXWvjuAt4H7RaSbiDhEZJCInNScgETEIyI+QAC3iPiqB2copRqm/0jUQcUYs8EYs6SBzdcDlcBG4GPgb8Ds+La/Am9hDwn/EvhXnX0vBTzYva/dwD+BXs0M620gABwPPBF//p1m7qvUIUu0YKFSSqlEpD0opZRSCUkTlFJKqYSkCUoppVRC6vAEJSKzRWSXiCxvYPtFIrJMRL4RkU9FZGxHx6iUUqrzdfggCRH5DlABPGeMGVXP9uOBVcaY3SJyGvBrY8wxTR03KyvLDBgwoN3jbanCwkIAsrOzOzkSpVQiWrp0aZExRj8gmqHDb9Q1xnzUWPkDY8yntRYXAX2bc9wBAwawZElDo4s7zmmnnQbAf/7zn06ORCmViERkS9OtFCT+TBJXAA1+0ovIVcBVALm5uQ0161CamJRSqn0k7CAJETkZO0H9oqE2xpgnjDHjjTHj9ZSaUkodXBKyByUiY4AngdOMMcWdHU9LPPTQQwDceOONnRyJUkp1bQmXoEQkF3uamUuMMWs7O56Weu+99wBNUEqp5lu6dGkPl8v1JDCKBD6zdQDEgOWWZV05bty4XXU3dniCEpG/AxOALBHZBtyFXcoAY8xjwJ1AJvCXeGkDyxgzvqPjbIlwfjk4BU/vFObOndvZ4SiluhiXy/VkTk7OiOzs7N0Oh+OQmX8uFotJYWFhXkFBwZPAlLrbO2MU37Qmtl9JF6pUamKGwie/wYSieId2J+uyPMR5KH0BUkq1g1GHWnICcDgcJjs7u6ygoGC/W47g0OpKHhDRPSFMKIo7x09o7W6euu8x/vSnP3V2WEqprsVxqCWnavHXXW8u0gTVRtauAADJx9iVF3at3sZnn33WmSEppdRBQRNUG0V2VQGQlJcJTuHHF8zglVde6eSolFJdXX5+vuvMM888rG/fvqNHjhw54vDDDx/+3HPPpXdGLK+//nrqO++8k1y9/Mc//jH7kUceyTzQvzfhRvF1NVZhFZLkwtHNgzs7iUhBZdM7KaVUI2KxGGeeeebgH/zgB8Xz5s3bBLB27VrPP/7xjwOWoCKRCG63u95t77//fmpKSkp00qRJlQA///nPCw9UHLVpD6qNIruqcPfwIyK4cpIpWruDe++9t7PDUkp1YfPmzUt1u92mdiIYOnRo+LbbbttlWRZXX31131GjRo0YOnRo3n333ZcFdi/n6KOPHjZ58uSBhx122MgpU6YcFovFAFi4cKH/qKOOGjZy5MgRJ5544pAtW7a4AY4++uhhl19+eb9Ro0aNmDlzZs+//e1vaWPGjBk+YsSIvOOPP35ofn6+a82aNZ7nnnsu+7HHHus5fPjwvPnz56fcdNNNve+8886eAJ9++mnS2LFjhw8dOjRv0qRJgwoLC53Vx77mmmv6jB49esSAAQNGzZ8/P6Wl74MmqDayCgO4spMAcPdMJsX4WP31yk6OSinVlX3zzTdJY8aMqapv24MPPpiVlpYWXb58+aqvv/561bPPPpu9evVqD8CqVauSZs2alb9+/foVW7du9b7zzjspoVBIbrjhhtzXXnttw4oVK1ZddtllRT/72c/6VB8vHA7L8uXLV/3mN7/ZOWnSpIqvvvpq9apVq1ZOnTq15O67784ZNmxY+NJLLy380Y9+tHP16tUrJ0+eXFE7nunTpx92zz33bFu7du3KkSNHBn7xi1/0rt5mWZZ88803q/7whz/k33333b1pIT3F1waxqgixigjuHn4A3Dn2z8d//3BnhqWUOshccskluV988UWK2+02ffv2Da1evdo/d+7c7gDl5eXOlStX+jwejxk9enTloEGDIgAjR46s2rBhgycjI8Nat25d0sSJE4eCffowOzs7Un3sadOmlVQ/37Rpk+fss8/uW1hY6A6Hw45+/fqFGouruLjYWV5e7vze975XATBjxozi8847b2D19vPOO283wPHHH195yy23eFr6ujVBtUGk0B7B56pOUD3ta4iRnVV4B6R1WlxKqa5t9OjRgddee6179fKcOXO27tixwzV+/PgRffr0Cd9///1bzz333D2193n99ddTvV5vzVB1p9OJZVlijJHBgwcHvvrqq9X1/a7U1NRY9fPrrrsu98Ybbyy46KKLyl5//fXU1vR6avP5fAbA5XIRjUalpfvrKb42sOIj+NzZSUSjUb5Y/SUvez7l5dd1FJ9SqvXOPPPM8lAoJH/4wx9qZsGuqKhwAEyaNKns0UcfzQ6FQgKwbNky7549exr8LB8zZkywpKTE9e677yYDhEIhWbJkia++tuXl5c7c3NwIwDPPPFMzSi81NTVaXl7urNs+MzMz2q1bt2j19aWnnnoq87jjjquo2661NEG1QbTC7iU7u3lYunQpb7/9NiEi5DuK2LlzZydHp5TqqhwOB/PmzduwcOHC1D59+owePXr0iIsvvnjAr3/9620//elPi4YPHx4cPXr0iCFDhoycMWNG/0gk0mDvxOfzmRdffHHDrbfe2nfYsGF5I0eOzPvwww/rHbBw2223fTtt2rRBI0eOHJGZmWlVrz/33HNL33jjjfTqQRK193n66ac3/eIXv+g7dOjQvGXLliXde++937bX+9DhFXUPlPHjx5uOLlhY+p9NVHyynb4zT+SVV15h8+bNTPX9D38vfpfsPj25/PLLic8nqJRSAIjI0rrzi3799debx44dW9RZMXW2r7/+Omvs2LED6q7XHlQbmICFw2dfxisoKCAnJ4fkbikc7R1Ofn4+mzdv7twAlVKqC9ME1QaxgIUjyUUkEqGoqIhevXrx5aqv6bU7Ba/Xy5dfftnZISqlVJelCaoNYkE7Qe3cuRNjDDk5OezaU0Sy8TFmzBhWrlxJVVW9tzIopZRqgiaoNogFLMTnoqCgAIBevXpx9rRzcIhw+IgxRKNRli1b1slRKqVU19QpCUpEZovILhFZ3sB2EZE/i8h6EVkmIkd2dIzNYeKn+Hbs2IHP5yM9PR1Hin0vWo/kDNLT09m2bVsnR6mUUl1TZ/WgngEmN7L9NGBI/HEV8GgHxNRi1af4qgdIiAiz//4MYA9Bz8zMpKSkpPGDKKWUqlenJChjzEdAY5/cZwHPGdsiIF1EenVMdM1jjImf4nOy49udOC17NokdZbsAiFVEyMjIoLi4mINlKL9S6uBVVFTkvPfee7ObbmlbvXq1Z8yYMcNzc3NHfe973xsYDAbb/Z6aRL0G1QfIr7W8Lb5uHyJylYgsEZElhYUdMvt7DROOQQyqIhFixqJylz1byN33/Q6AWKWdoEKhEIFAoENjU0qpliouLnY+9dRTPZrb/qabbup73XXX7dy6devytLQ066GHHspq75gSNUE1izHmCWPMeGPM+OzsZif+dhEL2jdZFxSVAVBVbIhGYjj8bhCIVoTJyMgA0NN8SqmEd/PNN/fNz8/3Dh8+PO/qq6/u21jbWCzGZ599lvrDH/5wN8Dll19ePG/evHavVZWok8VuB/rVWu4bX5cwTCCeoArtBIXlYtfWcu5/9Ldc5ZxMcrwHBXaC6tu30b+3UkrVuOWfX/dbW1Dub89jDs1Jrbpv6tj8hrbff//9284444yk1atXr9y9e7dj+PDhefW1e+GFFzb27t3bSk1NjVYXOBwwYEB4586dLZ6tvCmJmqDmAteJyIvAMUCZMWZHJ8e0j1g8Qe0qKgM/OKIedqwvJRAIUOkPkFkRoXt3ezJi7UEppbqS7t27x1avXt1gYbsdO3Z0SO7olAQlIn8HJgBZIrINuAtwAxhjHgPeBE4H1gNVwA87I87GVCeoqqhdLiU9sxs7NpQxa9YsCp9YRrQygsvlIi0tTROUUqpFGuvpdITdu3c7jjvuuOH1bXvhhRc2HnHEEcHy8nJndZn4zZs3e3r27Blu7zg6JUEZY6Y1sd0AP+6gcFqlOkGFxU5QfQdls/nrEkzM4EhxE/m2EqBmJJ9SSiWytLS0aGVlpQOa7kEBHHvsseVPP/1096uuumr37NmzM88444zS9o6pSw+S6EzVgyQsp4XP5yNnQHdClRa3/PSXfPLloppSHBkZGdqDUkolvJycnOi4ceMqhgwZMrKpQRJgX7N6+OGHc3Jzc0ft3r3bdeONN7b7bOyJeg0q4VUPksBrkZycTHKafX3QabxUmRAmaGGsGBkZGQQCAYLBID5fvTXClFIqIcybN29Tc9vm5eWFv/nmm1UHMh5NUK0UC1hEBYzbIiUlBX+aF4AfXXk92VVhSl/dQCxg4ffbA3ECgYAmKKWUagE9xddKsWCUiIGYhElOTsbfze5BVZWFampExYIWXq+duEKhUKfFqpRSXZH2oFopFogQjhkisZDdg4onqJf/9i/S3GWc7z8JE4zW9JqCwWBnhquUUl2OJqhWsiojBGJRrJjdg6ooKcREPiI5WobTmwbYpwG9ydqDUkqp1tAE1UrRSotK7GH/KSkpfDjnKUIVS3GJg7TyXuCzT/H5MrUHpZRSraHXoFopFtiboKzKctZ98Snd+06g98gbcafGk1JpuV6DUkqpVtIE1VqhKAGxE9SGzxbi9SfTZ9gEdmyrZGmFXeZ9x4rVeg1KKdUltLTcxpQpUw4bMGDAqCFDhow877zzBoRCoUOm3EZCM8ZANEYwnqB2rFjGqInfJTUzDbf48GZlYjAUrtuAwyE4nU7tQSmlElpLy21cdNFFJRs3bly+Zs2aFcFgUB588EEtt5EQogYxYLmiAJhIiH55o/CneRCc3HXH3eASYgGLTf9dgtfr1R6UUiqhtaTcBsAFF1xQ5nA4cDgcjB8/vnLbtm2HzGzmCc2E7cQUdUcRgFiMXoOHEV1rJ6GqPWFcyR58FSmsX7wIn8+nPSilVPO9+uN+7FrZruU26JFXxdmz2qXcxrhx42q+cYdCIXnppZcyH3jggXaf4FYTVCvE4gkq5orhiEF6zxz8aekkd9sNwK9vm8ltY6aS7E9ny7fr8PYZpAlKKdVlNGey2GqXXXZZ7rHHHlsxefLkivaOQxNUK5iwXd496opCZYReg4cB4I/Px5eT1QfxufB5ktm9bRvJg0bqKT6lVPM10tPpCE2V26juQd188829ioqKXG+99daGAxGHJqhWqJnJnBDGCtNriP13rJ6P77RJU3BsL8fj9BEo30O6y0llZWWnxauUUk1pabmNBx54IOv9999PW7hw4Rqn03lAYuqUQRIiMllE1ojIehG5tZ7tuSKyQET+KyLLROT0zoizIaE99ui9cDSARKP0HmL3oDw+J04nVBSU4vA5cdk1GMGytAellEpoLS238fOf/7x/UVGRa/z48SOGDx+e97Of/axXe8fU4T0oEXECs4BJwDZgsYjMNcbUzta3Ay8bYx4VkTzsCrsDOjrWhuxNUGEkGiUrdwAApf/8J+4Kw/p5XzLqB2cjUTv/xyJhvQallEp4LSm3YVnW0gMZC3ROD+poYL0xZqMxJgy8CJxVp40BusWfpwHfdmB8TQqXx2eQMBZulwuXx0Plos8puPMu3FYFXpKACIQNDqeTaKCKUChELBbr3MCVUqoL6YwE1QeofQFwW3xdbb8GLhaRbdi9p+vrO5CIXCUiS0RkSWFh4YGItV7hiuoEFcPrs6877XnjdRwpKaSPGkLEnUJ4ywaIGbr36EO4otzeLxzusBiVUqqrS9QbdacBzxhj+gKnA3NEZL9YjTFPGGPGG2PGZ2c3e4aONotUWsQwxAT8/mQAKhd9jv/oo0npk0XE353giq8ByOrZj1BZKaDz8SmlVEt0RoLaDvSrtdw3vq62K4CXAYwxnwE+oN2n0WitaFWECPZIvuTUFMLbthPJzyf5mGN46/03CTmSCG1YC0D37L5UlRYDOh+fUkq1RGckqMXAEBE5TEQ8wIXA3DpttgKnAIjICOwE1XHn8JoQDUQJEgEgNa07VZ9/DoD/2GPo0TsDxEnU2DfzpqVlY+Kn9rQHpZRSzdfhCcoYYwHXAW8Bq7BH660QkbtFZEq82c3ADBH5Gvg7MN0YYzo61oZEQxaBeIJKy8ig8vNFODMy8A4ZwpnfPw2AsLEn9vV5kpGonay0B6WUUs3XKdegjDFvGmOGGmMGGWN+F193pzFmbvz5SmPMCcaYscaYw40xb3dGnA0xoShBsRNUenY2VZ9/gf+YoxERklLt2SQiPvvalNfth5idoLQHpZRKVC0tt1Ft+vTp/fx+/xEHIqZEHSSR0EwkRgg72aR5k7B27iRp7FgAbrvzFwDEevQGwOP0ITHtQSmlEltLy20AfPTRR/7S0tIDdj+tJqjWiMQIip2gXAUFAPhG2BP/jjlyJADhrJ72dvHUnOLTHpRSKlG1tNyGZVnccsstfR966KFtByqmVme++IwQK4wx9U4oeLAyxiDRGCFHPNlstG+89o2w34abfn4jj/14AVZaD0xlDKmyAIOgCUop1Tx3fHJHv/W717druY3B3QdX/faE37ZbuY3f//73PU4//fTS/v37R9ozztpanaCMMdH4fHq5xpit7RlUIosEoziBoLGTTWzNatx9+uDsZk984XAIvhQ3EX93KA0QLdxNUkoqQRE9xaeU6hKamix28+bN7ldffbX7okWL1hzIONp67rA7sEJEvgBqpus2xkxpeJeuLVARxilCCPtnZNVqfHkjarZPmTKF47MvpmefARhrN9ESISm1G6UYIpED9kVDKXUQaayn0xGaKrexbt0675YtW3wDBgwYDRAMBh25ubmjtm7durw942hrgrqjXaLoQgLlEVwYwiaC2+kkvGUL3c48o2b7KaecgmOLj2DYgbFCRPeAv1s3xMQ0QSmlElZLym2MGzcueOGFF35dvez3+49o7+QEbRwkYYz5EFgNpMYfq+LrDlqBighOgYhEcQtgTM0ACYDzrzifQSP6EyiPIC5DLBglKTUNolFNUEqphNXSchsdoU09KBE5H7gP+AAQ4GERucUY8892iC0hBcrDJIvYCSo+OXn1AImvdn3FlW9fyXd2ncfQPUfh8DiIRSCpWzdM8R6dLFYpldBaUm6jtqqqqv+2dyzQ9lN8twFHGWN2AYhINvAucNAmqOCeMN1EiDgMSWELR3Iyrpwc8svzuf7967FKLUrCu4iGwEpx4Ag7SUrtRsyKaA9KKaVaoK0JylGdnOKKOcjvrQqW2aP3LDF4qqrwDB6EiDB7+WyCVpALnBfg9tg3Yxc5qujp6olXKpBoVIeZK6VUC7Q1Qc0Xkbew58sDuAC7ftNBK1xu94KiYnCXleEdOIigFWT+pvlM6j+JX/7PL9m+Zjev/t9/yY+W0dPVH08wDCZGWBOUUko1W5sSlDHmFhE5FzghvuoJY8y/2x5W4gqXhzEYog7wlO3Be+yxvLf1PSoiFZw9+GwAUjN9AJRKFHE4cVWEkJiO4lNKqZZo8xxKxphXgFfaIZYuIVIVwSIKIngiYTwDB/Hq+r/RJ6UP43PGc+qppyI4+P6QW4m6kiEK7tIKiMWwLE1QSinVXK1KUCLysTHmRBEpB2qXwbAHXhvTrV2iS0BWpUU4XqzQHYkQ7pfNF599wRWjrsAhDi644AIAPJu8eF1ZEAJHSQViooStqD1VkkhnvgSllOoSWjWgwRhzYvxnqjGmW61H6sGcnACiQYuw2AnKEzN8IVuImRgT+k0AYMaMGcyYMYNumUl4LfutiO4qgpg9Jj0anzhWKaUSSUvLbcRiMa6//vo+AwYMGDVw4MCRM2fObNFM6M3R6hF3IuIUkdWt3HdyfB6/9SJyawNtzheRlSKyQkT+1to425MVjoIVIxTvQfkzM1m44xPSvemMzBy5T9vUDB+xkF0byuwux+mw32q9F0oplYhaWm7j4Ycfzty2bZt7w4YNyzdu3Ljihz/8YUl7x9TqBGWMiQJrRCS3JfvFZ0GfBZwG5AHTRCSvTpshwC+BE4wxI4GftDbO9hSoiOACIvEelD8nh4+3f8wJfU7A6XACMGHCBCZMmEBqpo/KcrudKyh4PPFChjpQQimVgFpabuPJJ5/s8dvf/naH02l/9vXp08dq75g6Y7LYo4H1xpiNACLyInAWUHvepxnALGPM7vjxdu13lE4QKA/jFFPTgwp397A7tJv/6fM/NW2mT58O2CP5IvGZJsSdhCeewDRBKaWa8u2vbusXWreuXctteIcMqep9z+/ardxGfn6+d86cOd3feOON7hkZGdasWbO2jh49ul3vpemMyWL7ALXfpG3AMXXaDAUQkU8AJ/BrY8z8ugcSkauAqwByc1vUkWuVQEUEJ5Gaa1D5SeUIwvG9j69pU52g8leXEIkPHxFXEhh7H01QSqlE19RksQDhcFh8Pp9Zvnz5qmeffTZ9+vTpA5YuXdqu5Tfaeh/UhyLSHxhijHlXRPzYCaU94hoCTAD6Ah+JyGhjTGmd3/8E8ATA+PHjTd2DtLdgeRgXUYLYSeYb13aGZwynu697TZvqBNQt00eU+BBHdxJE7E6gXoNSSjWlsZ5OR2iq3Ma4ceOCPXv2DE+bNm03wCWXXFJ63XXXDWjvONo6WewM7B5MBjAIu3f0GHBKI7ttB/rVWu4bX1fbNuBzY0wE2CQia7ET1uK2xNtWdg/KIkgYRzTKp44N/G+vS/ZpM2nSJADee/d9EIi5hKg3CYL2GVDtQSmlElFLym0AnHbaaaXz589PHT58ePGbb76Z2r9//3afKqetp/h+jH1N6XMAY8w6EWlqFMhiYIiIHIadmC4EflCnzavANOBpEcnCPuW3sY2xtlmgPIJLYgQJ47IiFCRbHNNr37OTV155JQBOl4NuWUlYJoqVkoKzyk5Q2oNSSiWi2uU2Jk6cWPb4449va6z93XffXTB16tTD/vKXv/T0+/2xv/71r5vbO6a2JqiQMSZcfeOpiLjY98bd/RhjLBG5DngL+3TgbGPMChG5G1hijJkb3/ZdEVkJRIFbjDHFbYy1zQIVYfxOCBPBFbVwOt0c2ePIfdpcfPHFNc97D04jtLIYR3I6vqIqAEKBQIfGrJRSzdWSchtZWVnRDz74YP2BjKetCepDEfkVkCQik4BrgXlN7WSMeZM6k8oaY+6s9dwAN8UfCSNQHqGbQwgRQUyEMdlj8Lv3HWhTVWUnIr/fT+8h6YS+KcTvyyA9PslsZXl5h8etlFJdUVtLY9wKFALfAFcDbxpjbmtzVAkqWBHGJUJYLMKEGJ8zfr82p59+OqeffjoAvYd0J2LA6UzFH7QTVLCqcr99lFJK7a+tPajrjTEPAX+tXiEiN8bXHXQC5RFcCBGihBwRjuhxxH5trrnmmprn3bJ8GLcDMW588WtPAU1QSinVLG1NUJcBdZPR9HrWHRSClREcTrvce9AdYXTW6P3aVE8WCyAieNO8SGkQd8QCYwjqNSillGqW1s5mPg175N1hIjK31qZuQLvPx5QIotEYoSoLZ6qDiCOGJPlI86bt166srAyAtDR7W3J2Eq7SIJY3C2IxTVBKKdVMre1BfQrsALKA+2utLweWtTWoRBSssK8hCUJMICWr/kl/zzrrLAA++OADALrnphJYX0ppv+MQE6UioIMklFKqOVpbbmOLMeYD4FRgoTHmQ+yE1Re7JtRBJ1BT6t1ezunZv952N9xwAzfccEPNsq+7XV23IuMIiMWoDOg1KKVU4mlpuY3XXnstNS8vb8Tw4cPzxo0bN2z58uXe9o6praP4PgJ8ItIHeBu4BHimrUElokBFGGMMMaf9lh2WPbDedueccw7nnHNOzbIj2Q1AzJ2BIxojFG73m62VUqrNWlpu48Ybb+z//PPPb1q9evXK8847r+Suu+7q1d4xtTVBiTGmCjgH+Isx5jxgZBP7dEnBcnuao+pqun271z8bfVFREUVFRTXLjhQ7QXkExBgsnepIKZWAWlpuA6C0tNQJUFZW5uzVq1e7f7i1dRSfiMhxwEXAFfF17TFZbMIJVIRxESIkdg2NJF9Sve2mTp0K7L0G5Yz3oNL9gBGMFTvgsSqlurb3nlvVr2R7RbuW28jok1J1yqUj2q3cxmOPPbb5nHPOGeL1emMpKSnRxYsXr2rPeKHtCeon2IUF/x2frmggsKDtYSWeQHkEB2GqvyL4fL5629188837LFf3oNLTvUiJAyFGwAqQ5Ko/wSmlVGdrzmSxDzzwQM9//etf6yZOnFh5xx139Lzmmmv6vfTSS1vaM442l9vAnu4oRURS4kUIb2hqv2zsKB4AACAASURBVK4oUBHB64kQis872FCCOvPMM/dZFo8TXEK3bklQbE9VuLF0IyOzDsozoUqpdtBYT6cjNFVuo1evXtaqVauSJk6cWAlw6aWX7p48efKQ9o6jreU2RgPPYZfbEBEpBC41xqxoj+ASSbA8jEtChONvWUMJqqCgAICcnBzAvlnXmewmyetEcIIjxtrdazVBKaUSSkvKbUQiESoqKpzLli3zjhkzJvT66693Gzx4cLC9Y2rrKb7HgZuMMQsARGQC9rRHxze2U1cUqIggJkBYkhDA4/HU2+7CCy8E9l6DAnsknxeQmBPjiLK2aLVd3UoppRJES8ptuN1uHnrooS1Tp04dJCKkpaVFn3nmmWbPhN5cbU1QydXJCcAY84GIJLfxmAkpUB7GHQsSxo3H5aa6xEhdt956637rHMluYgELhzhBHGz8dvWBDlcppVqsJeU2Lr300tJLL720tOmWrdfWBLVRRO4A5sSXLyYBCgseCMHKCJ5omJBYDZ7eA5g8efJ+65zJbqziID6Xm4AjQuX6jRhjGkxySiml2n4f1OVANvAv4BXsqY8ub2onEZksImtEZL2I7N/l2NvuXBExIrJ/XYsOZGKGYEUEVyxKGAtfUsOjP/Pz88nP3/f6piPZTawiQnJqCgDpBREKA4UHNGallOrqWjtZrA/4ETAYuxbUzcaYZt2kJSJOYBYwCdgGLBaRucaYlXXapQI3Ei8n35mCVRGMATcQFgtfUsNnMS+55BKgzjWoFA8mHCW9TxbbduwiuzydtbvX0sPf7Ju2lVLqkNPaU3zPAhFgIXAaMAL7nqjmOBpYHx+Sjoi8CJwF1B0x8lvgD8AtrYyx3VTPw+fGQRiLtKSGT/Hdfvvt+62ruVk3Mwt2QGoohbW713JinxMPTMBKKXUQaG2CyjPGjAYQkaeAL1qwbx+g9jmwbcAxtRuIyJFAP2PMGyLSYIISkauAqwByc3NbEELLBCvsYoNuXISINHoN6tRTT91vXfV8fGkpds/LY/ysLtaBEkop1ZjWXoOqOZ1njLHaKRYARMQBPADc3FRbY8wTxpjxxpjx2dnNnoS3xap7UE6Hyz7F10iC2rhxIxs37jtOpHo2iWSvPdmvERcriw+6W8WUUqpdtTZBjRWRPfFHOTCm+rmI7Gli3+1Av1rLfePrqqUCo4APRGQzcCwwtzMHSgTitaAcDjcRoo0mqMsvv5zLL993nEh1D8oTs9/uiMtB+fYtlIe1NpRSKjG0tNzGPffck52bmztKRMbt2LGj5mxcLBZj+vTp/XJzc0cNHTo07+OPP271nIKtrQflNMZ0iz9SjTGuWs+7NbH7YmCIiBwmIh7gQqCmKq8xpswYk2WMGWCMGQAsAqYYY5a0Jtb2ULUnjDFRYi43SMOzSAD85je/4Te/+c0+66qvQbki9rDyqFPoW2RYVdzucysqpVSrtLTcxkknnVTxzjvvrO3du3e49vp//OMfaRs3bvRt3rx5+aOPPrrl2muvbfX1l7YOM2+x+CnB64C3gFXAy/GJZu8WkSkdHU9zVJYFiVCCcdhvV2MJ6qSTTuKkk07aZ534nOAUXKHqBBWjbxGsLG50LkallOowLS23ccIJJwSGDRsWrrv+tddeS7/ooouKHQ4Hp5xySuWePXtcW7ZscbcmprbeqNsqxpg3gTfrrLuzgbYTOiKmxpTs3kPMFGMc9gzkXm/DhSPXrFkDwLBhw2rWiQjOdC+m1L5cZxwwdHcSK4pXQOFa2PAeFG+AlJ4weCL0GXcAX41SKtG99eiD/Yryt7RruY2sfv2r/vean7RbuY2GjrNjxw73gAEDahJXr169wlu2bHH379+/xfWiOiVBdTV7SishupuYwx6F11gP6uqrrwb2vQ8KwNM3leDmMgCMw0Gf3X7+teUDWPBXwIC3G4T2wIKZcNh34KxZkH7gRiYqpVRDmlNuoyNogmqGYLmF09pDLF6KsbEEdc8999S73tM3hcDXhTiSHBiHE2eFh/xoCWXH/Yi0Y6+DtL4QLIMvn4MP/wiP/Q+c/ywMnND+L0gpldAa6+l0hKbKbTTWg+rVq1dk8+bNNbNp79ixw9Oa3hNogmqWaJXgipQTbaIWFMDxx9c/kbunXyoAbhFiDicBdyo9Swv46tTTOCktfrrXlwbHXw/DToeXLoYXzodpf4fBp7TvC1JKqTpaUm6jMVOmTCn9y1/+0mPGjBklCxYsSE5NTY22NkF1+CCJrsYKR5GIE48VIBp/txpLUMuXL2f58uX7rXf3SgYxeK0IxuEg4EvisCIHS3ct3f8gmYPgstchawi8+APY/mV7vRyllKpX7XIbzRkkMXPmzB49e/Ycs3PnTs/YsWPzLrjggv4A559/fln//v1D/fv3H3XNNdf0nzVrVqur7GoPqglVe+xrfZ5ICEsM0Pggieuuuw6ocw0qEsTx1s9xMx6PI4lKxx4CHh/jyrNYuLOB5JOcCZe8Cn892e5NXfUBpOjcfUqpA6cl5TZuv/32XbfffvuuuusdDgdz5szZ2h7xaA+qCVV7QgC4TYywRHE7XTidzgbb33fffdx33317V+xaBU9Ngi+fxdPLi5MkcLoJ+9wMLfGyongFQauB07kp2XDhC1BVAi9dAtZ+IzqVUuqgpQmqCTuLiwFwOVyEieB1N9x7AjjqqKM46qijoOAbmHcjPHoClOXDtBfxnDABd9SJy+Un5LTI/rYSK2bxTdE3DR+w11g4exbkL4L/dPq8uUop1WH0FF8TdhTadZvE5SYsFl5PAwmqqgTWvUPx0ldJKfoSb1UBOL0w/ocw4VeQnIk/EsXzuhu3+LBi5Ti/LSQp5GLJziUclXNUw0GMOhd2LINPHoT+J8CY8w/AK1VKqcSiCaoJRSWlgAeH22cXK/Ql7dtg1yp4fyasnQ8xC3fEyeelKXznsj/ZicWfUdNU3E782anEdpaQ6vAScQjfCeSy6NtFXDP2msYDmXgHbF0Er98EfcdDxsD2f7FKKZVA9BRfE8pKK4lShcPl33cmc2Ng4f3w6PGwaSEcey3MWMDGqe/SbcZcOHrGPsmpmr9vOhGiDEodyx5/MieU92RZ4TIqwhWNB+J0wbl/BYcDXrkSoq0atamUUl2GJqgmBPaEwZTjdCfZPSi/z05Ob9wE790NeWfDDf+F7/4W+hzJ4UccyeGHH97g8XwpSUQkSndvDuU9+zCoACxj8XlBMwoHp+fCmX+G7Uthwe/a8VUqpVTi0QTVhEilwWOV4nYkERKLpGQ/fPx/sGQ2nHAjTJ1tDwmPW7x4MYsXL27weB6PfYN1iieLquxs/Ou/xe/y8+n2T5sX0Miz4cjL4OMHYfMnbXptSilVrb3KbTz66KMZQ4cOzRs6dGjeEUccMfyzzz5Lauw4jdEE1QQJuPBEynA7vYSI4A/vtHtOo6bCqb+B+OwS1W655RZuuaXh0XbV91CJ00XUl0Zk8xZOTD+ST779BGNM84L633uge3949RoINXFqUCmlmqG9ym0MHjw49Mknn6xZu3btyl/+8pffXn311f1bG5MmqEYErACekB9PqBRxekEged1cyB4GZz2yX3ICeOSRR3jkkUcaPGZ1DyoiUdxiTz57clV/tldsZ1NZM++R86bA2Y9C6VZ4p95J4JVSqkXaq9zGpEmTKrOzs6MAJ598cmVBQYFn/72bR0fxNWJz0Va8UT+eqmJw2V8C/MEdcNEj4K6/1zpq1KhGj1mdoMLGwm/8GGBUoReS4f389xmY3szRef2PtwdmLJoFI86AQROb/bqUUomt5J9r+0UKKtu13IY7J7kqY+rQA15uo7aHH3446+STTy5rbcydkqBEZDLwEOAEnjTG3Ftn+03AlYAFFAKXG2NaPZ9Ta23avg0QJFxG1Gm/VcnDJkC/oxvc59NP7WtJDU0aW32Kr5wgae5MIgNyca7YwMjTRrIgfwFXjr6y+QGecgesexteuw6u/cyebFYppdqoPcptzJs3L/X555/P+vTTT1e39hgdnqBExAnMAiYB24DFIjLXGFP7zfgvMN4YUyUi1wB/BC7o6Fi/LdgF9ARTRSxeTdd/7BWN7vOrX/0K2L8eVLXqHlSlK0oPT0/KRwyj6uNFTJxxOQ9/PYvCqkKy/c28TulOgu8/Zk+lNP9X9owTSqkur7GeTkdoS7kNgM8//zzp2muv7f/GG2+sy8nJibY2js7oQR0NrDfGbAQQkReBs4CaBGWMWVCr/SLg4g6NMK6oqIxUehIlSMxpX29Kzurd6D6PP/54o9ure1CRJDcpVjo7U3uSUl7OhKpcHgYW5C/g/GEtmCmi73g44Sfw8QMw4kwYNrn5+yqlVFx7ldtYt26d57zzzhs0e/bsTWPGjAm1JabOGCTRB6j97WBbfF1DrgD+U98GEblKRJaIyJLC+JRE7alidxCDIeKIEom/U35/46eFhw0btk+597qqe1DRZHvC2ZCxE1bGim/pl9qPBfkLGty3QRNuhR558PpP7aKHSinVQu1VbuP222/vVVpa6rr++uv7Dx8+PG/UqFEjWhtTQg+SEJGLgfHASfVtN8Y8ATwBMH78+GaO0W6+cJnBQSUhtwuvI4bX6Wl0JnOADz/8EICTTqo35JoEZfx2jyxUHMAzeBBVn3/OyVefzN9X/53KSCXJ7uTmB+rywpRH4KlT4Z274MwHm7+vUkrFtUe5jZdeemkL0C5jBjqjB7Ud6FdruW983T5E5FTgNmCKMaZN3cTWCFpBnFU+vFYJUV8KQQnj9zZcqLDaXXfdxV133dXg9uoEJcl2gpJKg/eoo6haupSJPU8kEovw8faPWx5w33H2qL6lT9tTLymlVBfXGQlqMTBERA4TEQ9wITC3dgMROQJ4HDs57ZehO8K28m0kh9JJqtoJnhSCREjyNj3qc/bs2cyePbvB7U6nE5fLhcMH4ViUZFca4dF5mECAgV8X0d3bnfe3vt+6oE/+FXQfAPNugEigdcdQSqkE0eEJyhhjAdcBbwGrgJeNMStE5G4RmRJvdh+QAvxDRL4SkbkNHO6A2bJnKynhdDzBQjwOLwGJkOxvesaOgQMHMnBg4/cy+Xw+YmJRGYuR4kqnIi0Vd79+7HnpZU7qdxILty0k0prJYD3J9lx9JRvhg9+3fH+llEognTKThDHmTWPMUGPMIGPM7+Lr7jTGzI0/P9UY09MYc3j8MaXxI7a/zVuX4Y55cQZLcTt99ik+f9PXhd59913efffdRtukpqZSFagk6HCR4k6nKH8L6eefR9WSJUxmFOWRcj75tpXz7A08CY68FD59GLY3UE5eKaW6AJ3qqAGbVn8BgLHKcDt8BImQnNJ0gpo5cyYzZ85stE1qairl5eXEfC78rjQKN28i/ZxzwO2m30sf092dxhsb32h98JN+C8k9YO71WpZDKdVlaYKqT9l2dpbaU0yZWCVOpw8jhuRuqU3uOmfOHObMmdNom+oERTcPTnFSsa0QZ0YGPX5yI5Vvv8tPl/Zkwdb3m64R1ZCkdPje/bBzuV2FVymluiBNUPUIf/wAldEsACwJ4U2ypxBK6d50gurXrx/9+vVrtE1qaiqVlZU4u9v3QLkjbipKism4/HLSzz+f4W+s5OfPB/joi5db/yJGnGHXqvrwj1C4tvXHUUodEtqr3Ea1Dz/80O9yucY9/fTT3VsbkyaousoL2LD8RboFe+IgguV14E1KB2jWKb758+czf/78RtukptqJznSPz07hSqdw6yZEhJy77qTnnXcyeKeQ+bOHCG/bbwR+851+H7j9MPc6iMVafxyl1EGvvcptAFiWxS9+8Yu+J5xwQptmDtAEVdenD7PWKWRU9SLFKiGcnITDZd//1NQsEgD33nsv9957b6NtqhOUdAdjDCmudAq3bLbXOZ1k/GAa62dehiMYZsP0S4hVVbXutaT0gMn3Qv7nsPjJ1h1DKXVIaK9yGwD33HNPj7POOmt3VlaW1ZaYEnomiQ5XWQRLZrPmsNFkbOlNyp417HA5iIk9e0RyctM9qBdffLHJNjU9KF+Eqhik+3tRuGXfG7i/e+pVXP/fv/GrF3aw+29/I/PKFsxyXtvYC+Gbf8C7v4bBp0DmoNYdRynVYV599dV+u3btatdyGz169Kg6++yzD3i5jU2bNrnnzZvXfdGiRWvOP//8FkyJsz9NULUtvB+sIBu9fRkdScG7ez1WlpMKQjjEQUpKSpOHyMnJabJNdYKKEmJP1JDu6cGqrZ/v0ybdl07uyd9j2aLXOPzJJ0m/cBrOZpxi3I8ITPkzPHoC/GM6XPmuPTWSUko1oC2TxV577bX97r333m1NTQvXHJqgqhWthy+eIHb4xRR+W2qvi+wgxZXNbqkkIzW9yXn4AObNmwfAmWee2WAbv9+PiFBZVYnb4yGHZMp2FGCFw7g8e4tPXjTiIn514muMebaM3c8/T9aPrm7da0vra1fgfXEavH0HnP7H1h1HKdUhGuvpdIS2lNtYtmxZ8qWXXjowfhzXggUL0lwul7nkkktKWxqHJqhq79wBLh+rjjgfzyr7NF1Y9pDsGsI2qaRvVm6zDnP//fcDjScoh8NRM9S8R49+yK4qUp0ZFKxfS9+8vRV58zLzyDnqRFZ+/Cmj//43Mq+8AnG18k82/PR4Bd6/wGHfsUf5KaVUXHuV29i+ffs31c/PPffcAWeccUZZa5ITaIKyrXgV1rwJp/6aT0tXk1HVC7eEqfK7yEzvzSopp0fvns061D//+c9mtatOUAMHpcOuKtK9Pchf9c0+CQrgqtFX8fDhC8l7ZRcVH3xA6qmntvTV7XXqb2DrZ/DatZA9HLIGt/5YByFjDNGYIWoMsRhYsVjNz/rWxYwh2oJ11b8DwACmZv79+DpTe6n28v7bzd6dEREcAk4RRASnw152OASHCE7Zd7n2c7dT8LoceJxOPC4HHpfDXnY5cDns46lDQ+1yGxMnTix7/PHHtzXWfubMmT0efvjhnOLiYvfYsWPzTj755LL4TObtRhNUVQm8+TPoNRaOu45P3pnB0PB3San8lvK0VHqlZQHl9Ozd9LUlgKysrGa1S01NpaSkhKxjMgh/up0eaYexbeVyOHffdkf2PBLXCcdS8s4ivC/MaXGCsqIxglaMYCRKyIphnfIX+vzjDKxnz2XZ//6DSnc60ajBiu39cI7GYlhRe9mKGWLG7Lccjbc3xhAz2B/gxmAMRGP7P4/FP+Bjxv4dxlBznNrPYyaeKEyt57F997ePQfx32x/+e5/XOl58n9rHi9WKP2bqJpIWvbUHPRHwOPcmLa/Lic/tIMXrIrn64XGS7HXtsy7F6yTF6ybd7yYtyf6Z7veQ7HFqwktw7VFuo7ZXXnllc1viObQTVCwG/74aArvhkn9TEQ2xfMcKjiq/CF/JAnYkR4g43RCF7Ozm3b/2r3/9C4Bzzjmn0Xapqals2bKFHgO6sSoGaZ5eLF07n6gVwelyx8MzlAYinDvkGt45fBEXfPQF//r3Qooye1MRtNgTtKgIWZQHI/GfFhVBi/KQRSAcJRiJYtXzqXukXM/fPb/D+/IFXBn+JXto/UAbEfZ+K5f9v6E7BJzxb+L1fssXibettb+jzrGq2zjAJY6a4znj2+zj1d6/nnjqrHc66jzi210NrKvZVs86h+zdb591Ttn7u0So+9ksAoLUPN/nZ9317Lsd9h6vOiHXTtD7fDGovVwryceMIWwZwtEYYav6YX+RCVuxmvWh+PNQJEYgYlEZilIZsiiprKIyvHc5ZDV+r53LIbWSlof0JDdpfjfpSR66+91kpnjJTPGQleIhM9l+nuJ1aVI7hB3aCWrB72Dd23D6nyBnNJ9vfY/eu4dC1ImvYjkkQ0UsjAMHGRkZzTrkn//8Z6DhBBWNGYorQwSMh0AgwKtf59PLKeTEUrDCIa596DXyXT0pqghRXBkmGk8w2YPGctZn/2X7Yw9x/xE/RARSPC5SfS5SfW5SfC4ykj3kZvhJ9bnwe1x4XQ58bvtbr8/txOdy4nU78LrGsWZXb0Z//GM+6/EIWyY/C/7Mmg/m2h/SrpqfDpzO6g9ocDkcOAT98FA1ItEYVaEoFWH7i1JZIEJpVZjSqgilgeqfEcriywV7gqwuKKcsYH/Bqo/H5SAr2VOTvDKTvXYCq5XEsuLbMpI9eF1tHzmmEsehm6A+/j9Y+Cc44mI4yr7H6IVVL5BXeiweE8RKs0j1ZFNqVZDRLa3REXxhK0ZxZYjiijC3PPA0xZURHv1gA0UVob2P8jBFFSFKqsIYA30dpZzqgf/79ydcEs6hv8dDN3cW7p0bycnrz6g+3chK8dqPVC8ezxAWrryAiYtXcOGsYWQOHITD0ZbkcAH06Ubyy5eR98b34QcvQ/bQNhxPHercTgdpfgdpfneL9w1ZUXZXRmq+mBVX2P+eiuL/rorj69ftrKCwIkS4gd5aqs9FZjyhZSR7yEy2E1dGsp3Iqp8nYEKLxWIxcTgch9yJ5lgsJkC9f9BDL0FZYfum1UWzYNRUOOMhEGFxwWL+u/0rZpRMI3PnJ2xJ9jDusFP5VLaTlNabF7/YGk82YQorQhSVh2qWywL1zxju9zhrvt3lZvo5sn93slM8ZKV6yUhy8dWbW7l6mJvxucOIzt3IyF4nYSo+4vKLforLvf8/cs/1v8Rcehdf/e5GTp39etvfi2GnwfTX4e/T4ImT4Lu/hXGXg0MnGFEdy+tykpPmJCet6arVxhgqw1GK4//+imsltaKKMMWVYUoqQ+SXVPFVfim7K8P1nuoGSPW6yIgnq8xkr53QUvYmtu5+D2nVpyWT3HRLcuN2HpB/H8sLCwvzsrOzyw6lJBWLxaSwsDANWF7f9k5JUCIyGXgIcAJPGmPurbPdCzwHjAOKgQuMMZtb+/uCkSh7ghHCGz4mY+Fd+IuXs67/D3gv4ycUz19LYUUliwK/Y0jhCRjLiVX1FRFvlFLLQ4UnyH82GzZstEdOpvpcZMd7NkN7pnL8oOpejv0N7evPF5LiMky/8Bz8nsbf3uiWPFasWMG0c7/PN//ZQq/IQBaXzGXFB+8wdtLp+7WfOP58XpnyCnmvLuPVh3/C2dc/2PZTbP2Ohh8thFevgTduhqXP2pV5h3wXHAnz7VKpGiJCSnxgRv/Mpq+fGmPYE7DssxyVYYorwpRU7k1sJfHHtt1VLNtWSkkjCQ0g2eMk3e+hW5KbtCQX6Uke0uLX09KS7Eeqz0Wyx0WKz44zNf4zxVf/Z4JlWVcWFBQ8WVBQMIpDawq6GLDcsqx6p8qR2sNVO4KIOIG1wCRgG3YJ+GnGmJW12lwLjDHG/EhELgS+b4y5oLHj9ho80pz32+fZE4xQFoiwJxChIhAiI5jP0WYZZzs/4QjHer41GdwduZT5saMB8Hmr8PV5kaxQkLNX3Ehy6TeUOT4ht+dxbEqNQYqXUy64lOy0JDKTPfjcjX9oT5gwAYAPPvigyfdiw4YNzJkzh6lTp5Kysxuud7awLrSctcFFnHfHPWT07rPfPlErwkfnTKT7xiI+nTqM7974JwZ3b4fh4rEYLH8F3rsbyrZCWi7kTYEBJ0LPUfbNvnq9SR0CjDHsCVoUV4Tsa2bx62Zl8eeltZ6XBcL7rG9qoAjAlj+csdQYM74DXkqX1xkJ6jjg18aY/40v/xLAGPP7Wm3eirf5TERcQAGQbRoJdkBOX3PHJdftXWGq/yM1Q6BMzXN7hTH2HSYiDgxOMBaChcftZ7crQqmjkoum/YAhw5p/baYqPrFrcyaWjUaj/PnPf6a8vJy8vDzCy8tJsiBoVRK0KnB4fDicLvaO4Yo/NYZooAqJxogBxh6tYG+LJ5E2pRJj4jfd1Pd21xlWppTax37/auqsuOL+WzVBNVNnnOLrA9SexmMbcExDbYwxloiUAZlAUe1GInIVcBVAr169yE+Oz75han9+7v+BLXU+XesuiZTTo1sWJ004hcFDh7TgpTUvMVVzOp1ceeWVfPTRR6xcuZKIJ4KFRcxV/S2swdlEoOlpAZVSqkvrjB7UVGCyMebK+PIlwDHGmOtqtVkeb7Mtvrwh3qaovmMCjBs3zixduvTABt8Mzz//PAAXX3xxq4/R2N+kZiYCYxpt11EaiiERYquRSLGoQ57X59MeVDN1Rg9qO1C75Gzf+Lr62myLn+JLwx4s0aBEuR/nySftukttSVCNvZZEeZ1KKXWgdUaCWgwMEZHDsBPRhcAP6rSZC1wGfAZMBd5v7PpTInnnnXc6OwSllDoodHiCil9Tug54C3uY+WxjzAoRuRtYYoyZCzwFzBGR9UAJdhLrEtz13L+klFKq5TrlPihjzJvAm3XW3VnreRA4r6Pjag/PPPMMANOnT+/UOJRSqqs7lG4I6xDPPPNMTZJSSinVeh0+iu9AEZFyYE1nx1GPLOoMj08AGlPzJWJcGlPzJWJcw4wxqZ0dRFdwMM3FtyYRh26KyJJEi0tjar5EjEtjar5EjEtElnR2DF2FnuJTSimVkDRBKaWUSkgHU4J6orMDaEAixqUxNV8ixqUxNV8ixpWIMSWkg2aQhFJKqYPLwdSDUkopdRDRBKWUUiohdbkEJSKTRWSNiKwXkVvr2e4VkZfi2z8XkQEHOJ5+IrJARFaKyAoRubGeNhNEpExEvoo/7qzvWAcgts0i8k38d+43tFVsf46/V8tE5MgDHM+wWu/BVyKyR0R+UqdNh7xXIjJbRHbFZ86vXpchIu+IyLr4z+4N7HtZvM06EbnsAMd0n4isjv99/i0i6Q3s2+jfup1j+rWIbK/1N9q//DNN/1tt55heqhXPZhH5qoF9D8j7FD92vZ8Fnf3/VZdWXbahKzyw5+7bAAwEPMDXQF6dNtcCj8WfXwi8dIBj6gX/z96dx0dVnQ0c/z2zZU/YCQHCvogCoiyiiFtVtG51t7Zqq1Lr1tbaqq1b1fatrVtrra1L61qVuhUrintRFARkkZ2A7ARCgITsszzvH/cmDCEkE8gkk+T5FK/PAwAAIABJREFUfj7DzNx77txn5g7z5Jx77jkc4T7OwJktuHZMxwP/bYHPay3QpZ71pwPv4EyJdRQwu5mPZT7QpyU+K2AicASwOGrZH4Bb3ce3AvfXsV0nYI1739F93DGOMZ0C+NzH99cVUyzHuoljuhu4OYbjW+//1aaMqdb6B4E7m/Nzcl+7zt+Clv5eteZba6tBjQXyVHWNqlYBLwNn1ypzNvCs+/hV4CSJ4xwVqrpFVb9yH+8GluFMuNganA08p45ZQAcR6dFM+z4JWK2q65ppf3tR1Rk4AxFHi/7uPAucU8empwLvq+oOVd0JvA9MildMqvqeqobcp7NwpqdpNvv5nGIRy//VJo/J/b9+IfBSU+yrMer5LWjR71Vr1toSVF2z8dZOBnvNxgtUz8Ybd25z4ihgdh2rx4vIQhF5R0QObY54cCabfk9E5okz+3BtsXye8XIx+/8RaYnPCqC7qm5xH+cD3eso05Kf2Q9xarx1aehYN7Xr3WbHf+ynyaqlPqdjga2qumo/65vlc6r1W5Do36uE1doSVMISkXTgNeCnqlpca/VXOE1ZI4FHgTebKawJqnoEcBpwnYhMbKb91ktEAsBZwL/rWN1Sn9Ve1Gl3SZhrMETk10AIeHE/RZrzWD8ODAAOB7bgNKklikuov/YU98+pvt+CRPteJbrWlqAaMxsvEuNsvAdLRPw4X8gXVfX12utVtVhVS9zH0wC/iHSJZ0zuvja599uAN3CaXaLF8nnGw2nAV6q6tfaKlvqsXFurmzjd+211lGn2z0xErgDOAC51f+D2EcOxbjKqulVVw6oaAZ7cz75a4nPyAecCr+yvTLw/p/38FiTk96o1aG0JqmY2Xvev8ItxZt+NVj0bLzTDbLxum/fTwDJVfWg/ZbKrz4OJyFiczz3eSTNNRDKqH+OcbF9cq9hU4DJxHAUURTVFxNN+/8ptic8qSvR353LgP3WUmQ6cIiId3aatU9xlcSEik4BfAmepatl+ysRyrJsypujzlN/Zz75i+b/a1L4FLFfVjXWtjPfnVM9vQcJ9r1qNlu6l0dgbTs+zlTg9hH7tLrsH5z8wQDJO01Ee8CXQP87xTMCpsi8CFri304FrgGvcMtcDS3B6Ms0Cjm6Gz6m/u7+F7r6rP6vouAR4zP0svwZGN0NcaTgJJytqWbN/VjgJcgsQxGnvvxLnXOWHwCrgA6CTW3Y08FTUtj90v195wA/iHFMezrmJ6u9WdQ/VHGBafcc6jjE9735fFuH8+PaoHZP7fJ//q/GKyV3+TPX3KKpss3xO7uvv77egRb9XrflmQx0ZY4xJSK2tic8YY0w7YQnKGGNMQrIEZYwxJiFZgjLGGJOQLEEZY4xJSJagjKmHiPxURFJbOg5j2iPrZm7aPfcCS1FnZITa69biXB+2vdkDM6adsxqUaZdEpK87V9FzOKMJPC0ic915fH7jlrkR50LPj0XkY3fZKSLyhYh8JSL/dsddM8bEgdWgTLvkjja9Bmekilki0klVd4iIF+eq/xtVdVF0DcodE/B14DRVLRWRW4AkVb2nhd6GMW2ar6UDMKYFrVNnHiyAC93pF3w4E88NwxmyJtpR7vKZ7nCBAeCLZorVmHbHEpRpz0oBRKQfcDMwRlV3isgzOGM61iY4k8pd0nwhGtN+2TkoYyATJ1kViUh3nOlAqu3Gmb4bnMFrjxGRgVAzOvbgZo3UmHbEalCm3VPVhSIyH1iOM3L4zKjVTwDvishmVT3BnZvpJRFJctffjjNitzGmiVknCWOMMQnJmviMMcYkJEtQxhhjEpIlKGOMMQnJEpQxxpiEZAnKGGNMQrIEZYwxJiFZgjLGGJOQLEEZY4xJSJagjDHGJCRLUMYYYxKSJShjjDEJyRKUMcaYhGQJypg6uFPCq4g0OOK/iFwhIp81R1zGtCeWoEyrJyJrRaTKnZI9evl8N8n0bZnIQEQGi8h/RKRARHaIyHQRGdJS8RjTmliCMm3FN0DNTLciMhxIbblwanQApgJDgO7Al8B/WjQiY1oJS1CmrXgeuCzq+eXAc9EFRCRLRJ5zazPrROR2EfG467wi8oCIbBeRNcC369j2aRHZIiKbROQ+EfE2FJSqfqmqT6vqDlUNAg8DQ0Sk88G+YWPaOktQpq2YBWSKyCFu4rgYeKFWmUeBLKA/cBxOQvuBu+5q4AxgFDAaOL/Wts8AIWCgW+YU4KoDiHMikK+qhQewrTHtiiUo05ZU16JOBpYBm6pXRCWt21R1t6quBR4Evu8WuRB4RFU3qOoO4P+itu0OnA78VFVLVXUbTk3o4sYEJyK9gMeAmw7s7RnTvjTYQ8mYVuR5YAbQj1rNe0AXwA+si1q2DujpPs4BNtRaV62Pu+0WEale5qlVvl4i0hV4D/irqr4U63bGtGeWoEyboarrROQbnNrOlbVWbweCOMlmqbsslz21rC1A76jyuVGPNwCVQBdVDTU2LhHpiJOcpqrqbxu7vTHtlTXxmbbmSuBEVS2NXqiqYWAK8FsRyRCRPjhNbdXnqaYAN4pILzeh3Bq17RacBPOgiGSKiEdEBojIcQ0FIyKZwHRgpqre2lB5Y8welqBMm6Kqq1V17n5W3wCUAmuAz4B/Af9w1z2Jk0gWAl8Br9fa9jIggFP72gm8CvSIIaTvAGOAH4hISdQtt6ENjWnvRFVbOgZjjDFmH1aDMsYYk5AsQRljjElIlqCMMcYkJEtQxhhjElKbuQ6qS5cu2rdv35YOwxhj6jVv3rztqtq1peNoDeKaoERkEvAnwAs8paq/r7V+IvAIMAK4WFVfjVp3OXC7+/Q+VX22vn317duXuXP317vYGGMSg4isa7iUgTg28bljnz0GnAYMAy4RkWG1iq0HrsC5HiV6207AXcA4YCxwl3vxpDHGmHYinuegxgJ5qrpGVauAl4Gzowuo6lpVXQREam17KvC+O0XBTuB9YFIcYzXGGJNg4pmgerL3YJob2TMwZ5NsKyKTRWSuiMwtKCg44ECNMcYknlbdi09Vn1DV0ao6umtXO+dojDFtSTwT1Cb2Hh26F1Hz88RxW5OgIqWlbL7lFna9/kZLh2KMaQXi2YtvDjBIRPrhJJeLge/GuO104HdRHSNOAW5r+hBNcwmXlLDhqqspX7CAov++jT8nh7SjxrV0WMaYBBa3GpQ7b871OMlmGTBFVZeIyD0ichaAiIwRkY3ABcDfRWSJu+0O4F6cJDcHuMddZlqpHf98hvKFC+nx+/8j0K8vm372MyKlpQ1uZ4xpv+J6HZSqTgOm1Vp2Z9TjOTjNd3Vt+w/2TIVgWjFVpXjaNFLHjqXDOefgz85m/RU/oHTWLDJOOqmlwzPGJKhW3UnCtA6Vy5dT9c03ZJ5+OgCpRxyBJy2NkhmftnBkxphEZgnKxF3xtGng85FxyskASCBA2tHjKZkxA5uPzBizP5agTNwVvzudtKPH4+u4ZzCQtGOPJbRlC1V5eS0YmTEmkVmCMnEVzM8nuGED6ROO3Wt5+sSJANbMZ4zZL0tQJq7KFywAIGXU4Xst92dn4++TW7PeGGNqswRl4qp8/nwkOZnkoUP3WZc8ZCgVK1e0QFTGmNbAEpSJq7IFC0g+7FDE799nXdKQwQTXb7DroYwxdbIEZeImUllJxdJlpI4aVef65CFDQJVK6yhhjKmDJSgTNxVLlkAwSMrhh9e5PmnIEKfcCmvmM8bsyxKUiZvyhYsASBk5ss71/p498aSlUbliZXOGZYxpJSxBmbipXLkSb5cu+Lp0qXO9eDwkDR5MpdWgjDF1sARl4qZy5UqSBw+qt0zSkMFUrFxpI0oYY/ZhCcrEhYbDVOblkTRocL3lkgYPJlJcTGjbtmaKzBjTWliCMnFRtX49WllJ0uD6E1Qgtw8AwfXrmyMsY0wrYgnKxEXlylUAMSQoZ+LkqvUb4h6TMaZ1sQRl4qJy5UoQIWnggHrL+Xv0AK+Xqg1WgzLG7M0SlImLypUrCeTm4klJqbec+P34c3Ksic8Ysw9LUCYuKleubLB5r1ogN9ea+Iwx+7AEZZpcpLycqvXrY05Q/tzeVG2wBGWM2VtcE5SITBKRFSKSJyK31rE+SURecdfPFpG+7nK/iDwrIl+LyDIRuS2ecZqmVZm3GlRjr0H1ziVSVER41644R2aMaU3ilqBExAs8BpwGDAMuEZFhtYpdCexU1YHAw8D97vILgCRVHQ4cCfyoOnmZxFe50hm6KKmBi3Sr1fTks1qUMSZKPGtQY4E8VV2jqlXAy8DZtcqcDTzrPn4VOElEBFAgTUR8QApQBRTHMVbThCpXrkSSkwnk5sZU3t/bKVdlHSWMMVHimaB6AtF/Em90l9VZRlVDQBHQGSdZlQJbgPXAA6q6o/YORGSyiMwVkbkFBQVN/w7MAalctZKkAQMQrzem8oHevQAIWg3KGBMlUTtJjAXCQA7QD/i5iPSvXUhVn1DV0ao6umvXrs0do9mPipWrYj7/BOBJTcXbtYs18Rlj9hLPBLUJ6B31vJe7rM4ybnNeFlAIfBd4V1WDqroNmAmMjmOspomEduwgvH17oxIUgD8nh9CWLXGKyhjTGsUzQc0BBolIPxEJABcDU2uVmQpc7j4+H/hInWGt1wMnAohIGnAUsDyOsZom0tgOEtX8PXIIbrYEZYzZI24Jyj2ndD0wHVgGTFHVJSJyj4ic5RZ7GugsInnATUB1V/THgHQRWYKT6P6pqoviFatpOtUJKrmxNagePQhu2WLTbhhjavji+eKqOg2YVmvZnVGPK3C6lNferqSu5SbxVaxcibdjR7z7maRwf/w9eqCVlYR37sTXqVOcojPGtCaJ2knCtFKVbgcJ52qB2Pl75gBYM58xpoYlKNNkNBJxJilsZPMeuKOaA8Etm5s6LGNMK2UJyjSZ4MaNaFlZoztIAPjcBGU9+Ywx1SxBmSZzoB0kALwdOiApKQQ3WQ3KGOOwBGWaTEV1F/OBAxu9rYjU9OQzxhiwBGWaUOXKVfh798aTlnZA2/tzcixBGWNqWIIyTaYxkxTWxWpQxpholqBMk4iUl1O1di3JQw4iQeX0ILx9O5HKyiaMzBjTWlmCMk2icuVKiERIOuSQA34NX7bbk2/r1qYKyxjTilmCMk2iYtkyAJIPqT0nZez82d0BCObnN0lMxpjWzRKUaRIVy5bjycysGRHiQPiyswGrQRljHJagTJOoWLaM5KFDGz3EUTR/d7cGtcVqUMYYS1CmCWgoROWKFSQfxPkncCYu9GRlEbImPmMMlqBME6hauxatrCR52MElKHBqUUFr4jPGYAnKNIHqDhIH04Ovmi+7u9WgjDGAJSjTBMoXLkJSUkjq3/+gX8vfPdtqUMYYwBKUaQLlCxaQMnw44jv4+S992d0Jb9+OVlU1QWTGmNbMEpQ5KJGKCiqWLydl5MgmeT2/29U8uK2gSV7PGNN6WYIyB6ViyRIIhUgZdXiTvN6ea6HsPJQx7V2DCUpEvCKy/EBeXEQmicgKEckTkVvrWJ8kIq+462eLSN+odSNE5AsRWSIiX4tI8oHEYOKrfMECgKavQdm1UMa0ew0mKFUNAytEJLcxLywiXuAx4DRgGHCJiNQeB+dKYKeqDgQeBu53t/UBLwDXqOqhwPFAsDH7N82jfMEC/Lm5+Dp3bpLX83W3GpQxxhHrWe2OwBIR+RIorV6oqmfVs81YIE9V1wCIyMvA2cDSqDJnA3e7j18F/iLOUASnAItUdaG7n8IY4zTNSFUpm7+AtPHjm+w1velpeNLTCeZbTz5j2rtYE9QdB/DaPYENUc83AuP2V0ZVQyJSBHQGBgMqItOBrsDLqvqHA4jBxFHlqlWEt28n7ajah/Xg2LVQxhiIMUGp6v9EpDswxl30papui19Y+IAJ7v7KgA9FZJ6qfhhdSEQmA5MBcnMb1QJpmkDpzM8BSDv66CZ9XbsWyhgDMfbiE5ELgS+BC4ALgdkicn4Dm20Cekc97+Uuq7OMe94pCyjEqW3NUNXtqloGTAOOqL0DVX1CVUer6uiuXbvG8lZMEyr9/HMC/fvj79GjSV/XalDGGIi9m/mvgTGqermqXoZzfqmhZr85wCAR6SciAeBiYGqtMlOBy93H5wMfqaoC04HhIpLqJq7j2PvclWlhkcpKyubMIe2YY5r8tf3dswkVFKBB6xdjTHsWa4Ly1GrSK2xoW1UNAdfjJJtlwBRVXSIi94hIdeeKp4HOIpIH3ATc6m67E3gIJ8ktAL5S1bdjjNU0g/KvvkIrKkg7pmmb9wB8PbJBlVCBXaxrTHsWayeJd90OCy+5zy/CaXarl6pOq11OVe+MelyB02xY17Yv4HQ1Nwlo94cfIUlJpI0Z03DhRqq5Fip/K/6cA58A0RjTusXaSeIXInIeUN2e84SqvhG/sEwi03CY4unvkn7ccXjS0pr89X3uxIV2LZQx7VvMo3uq6mvAa3GMxbQSZXPmEi7YTubpp8Xl9aNrUMaY9qveBCUin6nqBBHZDWj0KkBVNTOu0ZmEVDxtGpKaSvpxx8Xl9T0ZGUhqqvXkM6adqzdBqeoE9z6jecIxiS5SWcnu6dPJOOEEPCkpcdmHiNjMusaY+A4Wa9qe4mnvEC4qosMFdfZtaTJ2LZQxJm6DxZq2aeeLLxIYOIDUcWPjuh8bTcIYE8/BYk0bU75oERWLF9P9jttxxvSNH1+PbELbtqHhMOL1xnVfxpjEFM/BYk0bU/jPf+JJTyfr7LPjvi9/92wIhwlt347f7XZujGlfYhpJQlX/B6wF/O7jOcBXcYzLJJiqtWvZPf09Ol5yCd709Ljvz5ftXgtl56GMabdiHSz2apz5mv7uLuoJvBmvoEziKXz6H4jPR6fLL2uW/dm1UMaYWMfiuw5nFIliAFVdBXSLV1AmsQS3bqPozTfJOu9cfF26NMs+bTQJY0ysCapSVauqn7gjjGs95U0bsuO5Z9FwmM4//GGz7dPboQOSlGQ1KGPasVg7SfxPRH4FpIjIycC1wFvxC8skinBREbteepnM004j0Lt3wxsAwXCQ6eum878N/2P5juWUhcrontqdkV1HcsGQC+if1b/B1xARuxbKmHYu1gR1K3Al8DXwI2Caqj4Zt6hMwtg5ZQqRsjI6X31Vg2VVlXe+eYeHv3qY/NJ8uqZ05fBuh5PmT2NTySamrJjCC8te4KwBZ3Hr2FvJCNQ/QIk/uwdBS1DGtFuxJqgbVPVPQE1SEpGfuMtMG6XhMDtfeonUceNIHjq03rK7q3Zz+2e389GGjzi086HcedSdTOg5Ya/rpQrLC3lu6XM8u+RZ5ubP5fFvPU7/DvuvTfmzu1M2Z26TvR9jTOsS6zmoy+tYdkUTxmESUMknnxDavIWOl3633nIbijfw3be/y4yNM/j5kT/nxdNf5Nhex+5zMW/nlM787Mif8expz1IZruTydy9naeH+J0r2dc8muG0bGok0yfsxxrQu9SYoEblERN4C+onI1KjbJ8COZonQtJidL/4LX3Y2GSeeuN8ym0s2c+V7V7KrchdPnvIkVxx2BV5P/SM/jOw6kudOe44UXwo//uDHbNy9sc5yvuzuEAoRLiw8qPdhjGmdGqpBfQ48CCx376tvNwGnxjc005KCmzdT+vnndDj/fMRXd0twQVkBV793NSXBEp44+QlGZ4+O+fVzM3P528l/IxQJce2H11JSVbJPGbsWypj2rd4EparrVPUT4FvAp+4oEluAXjhzQpk2qujttwHIOrvu4RZ3VOzg6veuZnv5dh7/1uMc0vmQRu+jf1Z/HjnhEdYVr+PeWfeiuveVC3YtlDHtW6znoGYAySLSE3gP+D7wTLyCMi1LVSmeOpWUUaPq7FpeXFXMNe9fw8aSjfzlpL8wsuvIA97XmOwxXDvyWqZ9M4038/YenMRqUMa0b7EmKFHVMuBc4K+qegFwaIMbiUwSkRUikicit9axPklEXnHXzxaRvrXW54pIiYjcHGOcpglUrlhB5ao8ss46c591pcFSfvzBj1m1axWPnPAIY7LHHPT+rhp+FaO7j+YPc/7A1tI9ycjbsSPi9xPK33LQ+zDGtD4xJygRGQ9cCrztLqv3TLiIeIHHgNOAYcAlIjKsVrErgZ2qOhB4GLi/1vqHgHdijNE0keJ33gWvl4xJk/ZaXhGq4IaPbmDJ9iU8cNwDTOg5oUn25/V4+c3RvyEUCXHf7PtqmvrE48HXowfBzZagjGmPYk1QPwVuA95Q1SUi0h/4uIFtxgJ5qrrGHSbpZaD2PA1nA8+6j18FThK3b7KInAN8AyyJMUbTREo++pDU0aPxdexYs6wqXMVPP/kpc/Pn8rsJv+Ok3JOadJ+5mblcd/h1fLLhEz7d9GnNcn9ODsHNm5t0X8aY1iHm6TbcyQkfE5F0N+nc2MBmPYENUc83usvqLKOqIaAI6Cwi6cAtwG/q24GITBaRuSIyt6CgIJa3YhpQtX49lavyyDhpT9fyUCTEL2f8kpmbZnL30Xdzev/T47LvS4ddSt/Mvvxxzh8JRoIA+HvmENy0KS77M8Yktlin2xguIvNxajNLRWSeiDR4Duog3A08rKr79j2OoqpPqOpoVR3dtWvXOIbTfuz+8CMA0k90akjhSJhff/ZrPlz/IbeOvZVzB50bt337PX5+MeYXrC1eyyvLX3GW5eQQKiggUlXVwNbGmLYm1ia+vwM3qWofVc0Ffk7UsEf7sQmI7gLWy11WZxl3hPQsoBAYB/xBRNbiNC/+SkSujzFWcxBKPvyQpCFDCPTqSSgS4o6ZdzDtm2n85IifcOkhl8Z9/8f2PJajc47mrwv/yq6KXfhznEp3yJr5jGl3Yk1Qaapac87JvTYqrYFt5gCDRKSfiASAi4GptcpMZc8wSucDH6njWFXtq6p9gUeA36nqX2KM1RygcHExZfPnk37C8QQjQW6ZcQtvrXmLG0bdwFXDGx4stimICL8Y/QtKg6U8vvBx/D1zAOw8lDHtUKwJao2I3CEifd3b7cCa+jZwzyldD0wHlgFT3A4W94hI9dWfT+Occ8rDGZ1in67opvmUzp4N4TCBo8dx08c38d6697h59M1MHjG5WeMY2HEgFwy+gFdWvMK2TKdHX5WdhzKm3Yl1NPMf4nRYeB1nosJP3WX1UtVpwLRay+6MelwBXNDAa9wdY4zmIJXOnAmpKdyQ/ygLdn7N7eNu56KhF+1/g7IdsPJdWP8FbFsOZdtBI5CUAR37Qc4o6Hcc9DwCpHEDj1wz8hqmrp7KE/mvcZnXazUoY9qhehOUiCQD1wADceaC+rmqBpsjMNP8dsz4mK9zI6zYncdDxz/EyX1Orrtg/mKY+QgseQMiIUjuANnDIecIEA9U7IKti2GZ26KblQujvgdjr4bUTjHF0iWlCxcPvZhnFj/DZV07WU8+Y9qhhmpQzwJBnBrTacAhOJ0WTBsSjAR55YNHOHLzNlYf2ZkXT/8ngzoO2rdg2Q744C746nkIpMHYyTD8AqemVFcNqXQ75H0Ai6bAJ79zktoRl8MxN0JmToNx/eDQH/DK8lfYkh4ixWpQxrQ7DSWoYao6HEBEnga+jH9Ipjkt3r6Yuz+/m9wPl3EkcM1Vj9O5ruS0+mN488dQWgDjr4OJN0NKx33LRUvrAiMvdm5bl8Lnf4Y5T8JXz8HEn8P468GXtN/NOyZ35HvDvsfK5MfpsWH9wb1RY0yr01AniZrmPLfTg2kjNuzewB0z7+DSaZeys2In36s4HF92Np0GH7Z3QVWY+Wd4/juQlAlXfwSn/rbh5FRb92Hwnb/B9XNhwAnw4T3w2DinhlWPy4ZdRnHHJCLbCtCgtS4b0540lKBGikixe9sNjKh+LCLFzRGgaVqbSzZz9+d3c9YbZzFtzTQuPeRS3jj7DdIXryN1zJi9Z8GNhOG/P4X374BDz4HJn0CPAx+5HIBO/eDiF+F7r4PHBy+cB1NvhIq6v05ZSVkMPnQCHoWlS2cc3L6NMa1KvU18qlr/1Kim1cgvzeepr5/itVWvIQgXDLmAq4ZfRbfUblSuWUO4sJDUsVEjk4dD8PrVsOR1mHATnHRno3vi1WvgSXDNZ865qc8fhdUfwdl/gf7H71N0wtjzKfjbB7z96dMcOrJpxwA0xiSuWLuZm1aqpKqEJ75+gheXvkiECOcOPJerR1xNdlp2TZmyL51Ti2ljxzoLIhF46ydOcvrWb2BCnPrF+JPh5Htg6BnO+a3nzoYxVznLAnuuA88acAgFwNaVC1lauJRhnWsPim+MaYtivVDXtDIRjfDGqjc4440z+OfifzKp3yT++53/csf4O/ZKTuAkKF/37vhzc51zTu/dDgtegONuiV9yitZ7rFObOuo6mPM0PH4MrJ9Vs9rXrSuSnERusZ/HFz4e/3iMMQnBElQblF+az9XvXc2dn99Jr4xevPTtl/jthN/SM732YPLO7LmlX87Zc/5pxgMw6zEY+yM4/rbmC9qfApN+B1f8FzQM/5gE798JoUpEhEDvXI4I5vDJhk9YWri0+eIyxrQYS1BtzOwtszn/rfP5evvX3D3+bp4/7XkO63LYfstXfbOW8Pbtzvmnxa/Bx/fBiItg0u+b9pxTrPpOgB9/DkdeDjP/BE8cD1sW4u+TS/YuyAhkWC3KmHbCElQb8vqq1/nR+z+iS3IX/n3mvzlv8Hl798qrQ835p36Z8OZ10HscnPUoeFrwq5GUAWf+Cb77b+fi4CdPJCBbCW/YxGVDv88nGz5hWeGylovPGNMsLEG1Ec8vfZ67Pr+Lo3ocxQunv0CfzD4xbVf25Zf4unbBP+MmZxiii16o9+LZZjX4FLj2Cxh2DoGdn6FVVVzEAKtFGdNOWIJqA6asmMIf5vyBk/uczKMnPkp6ID2m7VSVsjlzSO1SjpQXOtcnpXeLc7SNlNoJzn+awLdvAiDw1GV8P6kXH2/42GpRxrRxlqBauffXvc99s+5jYq+J3D/xfvxef8zbVq1dS6iggNS0TU4LAn5JAAAfbUlEQVSTWs6oOEZ6cAITvwtAVdZRXLr4PTIiyuOf3e30OjTGtEmWoFqxZYXL+NWnv2JE1xE8eNyD+D2xJyeAsrdfACD1uFOd8fISmC87G/H7qep4DJmXv8P3Q0l8vGspy148y5nqwxjT5liCaqWKKou48eMb6ZDcgUdOeIRkX3LjXqBkG2XTXsCXJgQuS/zJisXrJdC3D1VrvoE+47n0e++T4Qnwt5KV8Ph4mHoDFNuI58a0JZagWiFV5e7P72Z7+XYeOf4RuqR0adwLRCLo65Mp3Qyp4ycgSbGds2ppgQEDqVy9GoDMlE58f8RVfJQSYNkRl8CCl+DPR8AHv4HyXS0cqTGmKViCaoVeXfUqH6z/gJ+M+gmHdjm08S/w+Z+pmv8p4QoPacef0vQBxknSwIEEN2wgUl4OwKWHXEqGP4O/pfnhhrlwyJnw2UPwyAj4+P+gfGcLR2yMORiWoFqZ1btW84cv/8DROUdz2aGXNf4FNsyBj+6lVEYDkHrUUU0cYfwkDRwAqlR98w0AmYFMvj/s+3y04SOWRcrhvCfhR59C/4nwv9/Dw8Phw3uda6mMMa1OXBOUiEwSkRUikicit9axPklEXnHXzxaRvu7yk0Vknoh87d6fGM84W4tgJMgtM24h1Z/Kbyf8Fo808vCV74JXfwiZOZRV9Mefk0OgV6/4BBsHSQMHAlCZl1ez7NJhbi1q4d+cBT1GONdyXTPTGTH90wfhkeHwwd3ODL/GmFYjbglKRLzAYzhTxQ8DLhGR2sNQXwnsVNWBwMPA/e7y7cCZ7my+lwPPxyvO1uTZJc+yYucK7hx/Z+PPO6k6HQl2b0a/8xRl8+aTOm5cfAKNk0BuLvh8VOatrlkWXYtaviOqN1/2YXDhs86FvoNPhc8ecZr+3rsDSgpaIHpjTGPFswY1FshT1TWqWgW8DJxdq8zZwLPu41eBk0REVHW+qlZ3yVoCpIhIggxv0DI2FG/gbwv/xkm5J3FS7gHMiTT3H7BsKpx4B5XlWYSLikg7qnUlKAkECPTts1cNCuqoRUXrdgic/w+4bjYM/TZ88Rf40whnIFqrURmT0OKZoHoCG6Keb3SX1VnGnVK+COhcq8x5wFeqWll7ByIyWUTmisjcgoK2+1exqnLvrHvxeXzcNvYARhjPXwzv3gYDToKjb6R0ljOVRWurQQEkDRhI5eq9E1R1LerD9R+yePviujfsOsQ5R3XtbGf+qZl/dmpU798FpYXNELkxprESupOEiByK0+z3o7rWq+oTqjpaVUd37dq1eYNrRm9/8zZfbPmCG0fdSPe07o3buLIE/n0FpHSA7/wdPB7KZn+Jv08u/uzsBjdPNEkDBxJcv4FIRcVeyy879DI6JXfigbkPoPWNLtF1sJOorvsShp7ujJhec47KEpUxiSSeCWoT0DvqeS93WZ1lRMQHZAGF7vNewBvAZaq6mnaqqLKIP875IyO6jOCiIRc1bmNVePsm2LEaznsK0ruioRBlc+aQNq719N6LljRkMKhSuWLFXsvT/Glcd/h1zNs6j483fNzwC3Ud7Hwm182GIac556j+NMK5jsp6/RmTEOKZoOYAg0Skn4gEgIuBqbXKTMXpBAFwPvCRqqqIdADeBm5V1ZlxjDHhPTj3QYoqi7hz/J14Pd7GbbzgRVj0ijMzbr+JAFQsW0akpITUcWPjEG38pQwfDkD54n2b8s4ddC79svrx8LyHCUaCsb1g1yFw/tNw7Sy3M8XDbo3qN1ajMqaFxS1BueeUrgemA8uAKaq6RETuEZGz3GJPA51FJA+4Cajuin49MBC4U0QWuLcEG2Y7/ubkz+GNvDe47NDLGNJpSOM23rYM3r7ZSUwTf1GzuPQL5/xT2tjWmaB82dl4O3em4ut9E5TP4+OmI29ibfFaXl35auNeuNtQpzPFtV/AoFOcRFVdo7JEZUyLkHrb61uR0aNH69y5c1s6jCZTFa7ivKnnEYwEeePsN0jxpTRi41J48kQoK3SuB8rYc95q7aXfQ8vL6ff6a3GIunls+NE1BDdvov9bb+2zTlW58r0ryduZx1vfeYuspKwD28m25TDjD7D4dfCnwrjJMP4GSKvdh8eYxhGReao6uqXjaA0SupNEe/bU10+xtngtdxx1R+OSkypM+yUUrIBzn9wrOYWLiiifP5+04ybGIeLmk3zYYVSuXkOktHSfdSLCL8f8kqKqIh6d/+iB76SmRjULhkxyr6M6zKmVFrbbU6LGNCtLUAloTdEanvr6KU7rdxrH9DymcRvPfRoWvADH/RIGnLDXqtKZMyESIX1iK09Qww+DSISKZXVPWDi001C+O/S7TFkxha8Lvj64nUUnqmHnwFfPwqNHwr8uhm9m2HxUxsSRJagEE9EI93xxDym+FH455peN23jd5/DOLTDoVKdjRC0l//sf3g4dSBkxoomibRkphx0G1N1Rotp1h19H15Su3DvrXkKR0MHvtNtQ+M7j8NPFzjm9jV/Cs2fCY+Oc2lXxloPfhzFmL5agEszrq15n3tZ5/Hz0zxs3nFHRRphyGXTsC+c+AbV6/GkoRMmMT0mbMAHxNrI3YILxdemCPyeH8q/m77dMeiCdW8bewrIdy3h5+ctNt/OM7nDir+FnS+Csv0BKR/jgLnh4mJOwZv8ddm1o+HWMMQ3ytXQAZo8tJVt4YO4DjMkewzkDz4l9w6pSePlSCFbAFW87F+XWUjZnDuGdO8k4+eQmjLjlpI4bR8lHH6GRCOKp+++sk/uczISeE3h0/qOckHsCPdNrD2RyEPwpcMT3nVvhalj4Eix7C975pXPrMgT6HgN9joG+EyCj9V0U3VTCEaUqFKEyFHbvnceV7uOaZcEwVeEIwXAEAEEQcV5DRBBAxFnu9UCSz0uy30uy3+Pee0kNeMlK8ZPsb91/hBmHJagEoar85ovfOE18R98T+0jl4RD8+weQvwgufsm5rqcOxe9OR1JSSJ94bBNG7cRdHgxTUhGiPBh2blVhKoIRKqKehyIRQhElXOsWvUwBj4DH/THyeKTmB6lmuTg/Vt06D2Bg0RtMffVjQgOG4Pd5CHgFv9dTcwv4hPP73Mi8/Cu46aNb+f3RfyXZ76tZn+Rz7r0eObgPofMAOPF257Y9D1ZMc85PLZrijIEI0GkA9BoDvUZDzyOh+2HgCxzsx18vVaUqHJUAapJBVKIIRqgKh6kM7r2+so5topNJVXWZYISqsLu++nHNvVMmFGn+83TJfg8dUgJ0SPWTleKnY2qAHh2SyclKoUeHZHpkpZDTIZluGckHf/xN3FiCShBv5r3JzM0z+dW4X9ErI8YpMKpHilg1Hc542OltVlexUIjd779PxgnH40mpv0egqrKrLMi23ZVs213B1mLnfltxJbvKqigqD+51Ky4PUeX+xducOpWn8CIw45V3eG1QRb1lfZlnUN5zCqf+8zdUFZ6wz3qvRwh4Pfi9QsDnJeAVAr7qJLfnPlDz3Cnn90pNkgt4PW6S9BDwnU4g9wz8fcJ03r2c7J3z6LZrPl2Xf0DqIqe5MeQJUJhxCAVZwynIHM7WzMPY4e9OMAyhiPMDHwwpoYhTowiG1b2P7JUIaieOyloJqCkEfB6SvB6S/B6SfF7nuXsL+Dyk+J1aS/SyJJ93r8c12/idzyjJH71+T/kknwef1/njTNX5o8Xph6KoUvM8HFEqQmEqgk5irAiGqQiFKasKU1QeZFdZkF1lVc59eZBV23YzY1UBZVXhvd6b3yv065LGoG4ZDOiWzqBu6Qzslk6/LmlWC0sAlqASQH5pPn+c80dGdx/duOGMZjzg9Co79ucw+of7LVY2Zw7hHTvIOHVSTQJav6OMdTvKWF9Y6jwuLGPjznIKdlfWmXAyknx0TAuQleL8RdojK4VM93Fmio/MZD8pfi8pAS8p/j1NLykBL8nuD5TXI3g9gs8jeNx7r0fwinMvIqgqEYWIOj9I1fdKreURpWDli1yTsYMbf3E8wXCEqtCeH/Gq6h/1UISq0Che+GYrS/iAa8aeRmf/gD1lQkpVOEww7DRDVdc4glH3le59WVWIonLdZ7nzOnv2uTcPMMa9KTkUcrgnj1GePA7fmcfwXa9wmLwAQIFmsiAykFk6jM84gk3envi9zg92wOvB59YQoxNBh9TAfhNB0n4SQWC/icQpH/DuSSQBrwdPG6lhqCrF5SE2F5Wzpaiczbsq2LCzjNXbSliyuYhpi7fUdMr0eYShPTIY2auDc+vdgYHd0q221czsQt0WFoqEuOq9q1hauJTXznyN3pm9G94InNG4378DRl4C5zxOTWO9S1XZWlzJiq274bd30WHBLO654o+sLAqxu2LvXm1dM5LI7ZRK744pdM9ymj26ZybRLSOZbhlJdMtMIjWQeH/L5N97H7tef53Bs2fhCdTfXFZUWcR5U88jxZfCS99+ifRAelxiqm5Wq06O0tDvWTiIt2Ap/vyv8OfPx7NpDlLojtbesR8ccgYMvxCyh+9zjE3TqgiGWVNQSl5BCcu3FLNoYxELN+6q+f+SFvAyvFcWR/XvzPj+nTk8twNJvsbXsuxC3dhZgmphj8x7hKcXP83vJvyOMwecGdtG1cnp0HPh3CcJ4mFF/m6+3lTE4k1FrNy6mxX5uymuCJFeVcaL797DjIHjmXXWD+nfJZ0+nVPJ7ZRKrnufiMknFiUzZrBh8o/o9dfHyDix4UmX5+TPYfJ7kzm659H8+YQ/N35sw+aycx2seg9WToc1H0Mk5HS6GHUpjPo+pHZq6QjbjUhE+aawlIUbdrFwwy7mrd/Jks3FqDrnuUb36cT4AZ0ZP6AzI3t1iKmGZQkqdpagWtCMjTO47sPrOH/w+dw1/q6Ytol89ic8H9zJ+pxJ/KPbbSzYVMrSLcVUuecbMpN9DM3OZHB2OoO7Z3DYrHdJ+fuf6PfmGyQPHRrPt9PsNBhk1bETSTvmGHo++EBM20xZMYV7Z93LZcMu4xdjftHwBi2tbAcsfRMWvgIbZoEvGYZfAGMnO9Pbm2ZXVBZk9jeFfL66kFlrClmevxuAjql+jhvclROGdmPioK50TKu7Vm8JKnaWoFrI5pLNXPDWBeSk5/DC6S+Q5K17wuCyqhALNuxi3jfbGbTg90wqeYP/ho/iJ8HrSA4EOKxnFiN7d2B4zyxG9upA704piNsUpJEIa848C09aGv2mvNKcb6/ZbLnrboqmTmXwzM/wpKbGtM3vv/w9Ly57kbvH3815g8+Lc4RNKH8xzHnS6R0YLIPc8TDhZ87gttb812IKSyqZubqQT5Zv45OVBeworcIjcHjvDpwwpBsnDO3GoTmZNf8vLUHFzhJUCyiuKubydy5na+lWXj7jZXIzc2vWbS+pZO7ancxdu4M563ayZFMR/kg5f/Y/xsneeczseiFbj7qdEb070b9Ler0nsIvff59NN9xIzh//QNaZMTYftjKlX37J+ssup+dDD5J5+ukxbROKhLj+w+uZvWU2Dx3/ECfk7tuzL6GV74QF/4JZj0PRBuf81LE/h0PO2ucCbdO8IhFl0aYiPl6+jU9WbGPhxiIAumcmcdIh3Tn5kO6ceEh3S1AxsgTVzCpCFVz74bXM3zafx096nOzAcOas3cHctTuZs24HawqcAVADPg+H9+rAyd13c/H6u0jftRyZdL8zqnYMVJVvzj0PLSuj/9v/RXyt8zxTQzQSIe+kbxHo3Zs+zz0b83YlVSVMfn8yy3cs55ETHmFir1Y4PmE46NSmPnsICvOgy2CYcBMMPx+8/paOzgAFuyv538oCPly2lRkrCyitCrPu/jMsQcXIElQzKq0q5+rp1/H1jrkMlsmsWz+U7SWVAGSl+BnTtyOj+3ZiTN+OHJaTSdKSKc7o2V6/M3zR4FNj3lfxe++x6caf0OP3/0eHcxoxKkUrVPjMM2z7/f30eelfpI4aFfN2RZVFTH5/Mit3rOS+Cffx7f7fjmOUcRQJO+epPn0Iti6GDrlwzE/h8EvBn9zS0RlXZSjMrDU7OH5IN0tQMbIEFUe7K4J8tX4X89buYPa6jSwJ/xlJWUP55vPJ9hzLmL6dGN23I2P6dmJg16jmut358O5tsOR16DPBSU5ZsQ/TEy4pZc0ZZ+DNSKffG2+02dpTtUhpKXknfYuUUaPo/fhfG7VtSVUJN358I3Py5/DjkT/mmpHXxD6KR6JRhZXvOtfHbZoLad1g3I9gzJXOmIEmIdg5qNhZgmoiqsqmXeXMW7fTOYe0bifL853uqN6kAjL7vEjYW8AFfW7mqlHn0SOrjhEdQlUw+2/wv/shXOVMmTHhpkafV8i/9z52/utf9H3pX6QcfngTvcPEVvDXv7L9z4+S+8wzpB01rlHbVoYrufeLe/nP6v8woecE7j3m3sYN1JtoVOGb/zmXI6z+EPxpzpiBR10LHfu0dHTtniWo2FmCOkAFuytZtHEXizYW8fWmIhZtLKpprksNeDkityNH5HYglPYlr639C0neJB447gHG9qhjqvVQJSx8GWY+AjvWwOBJcOrvnDHeGqnorbfY/Itf0vH73yf717862LfZakTKy1lzzjkQCtN/6n/wpKU1antV5ZUVr/DA3AdI9aVy85ibObP/mTU9r1qt/MXw+aOw+FXQiNPj78gfwKCTrUNFC7EEFTtLUA0IR5R1haWs3FrCyq27WbzJSUhbipzx30RgYNd0hvfKYkTPLEb37cTQ7AxWF63igbkPMGvLLI7odgT3T7yf7LRaI1oXbYJFL8PsJ6AkH3ocDif8GgafckCxln7+ORt+dA0phx9O7tNPIQ2MrtDWlH31Fesu/R4Zp5xCzwcfOKCmzbydedz1xV0sKljEYZ0PY/KIyRzf+/jWn6iKNjmTWc5/AUq2QkYOHPodOOw86HmEdVNvRpagYmcJylVUHmR9YRnrdpSyrrCMvG0lrMjfzeqCkr0G3ezXJY3hPbMY0SuLEb06cGhOJmlJzg+hqrKwYCHPL32e99e9T3ognRtG3cBFQy7ac15j13rI+xAWvwZrPwMU+h3nXM/S//gD+qFQVXa9MoX8++4jqV8/cp97Fl/H9nnOofCfz7Dt/vvJPP10evzut3iSG99JIKIR/pP3H55Y9AQbSzYyuONgLhl6Cd/K/RYdkvedyqRVCQed81TzX4S8DyAShKxcGPQtGHAS9JsIyZktHWWbZgkqdnFNUCIyCfgT4AWeUtXf11qfBDwHHAkUAhep6lp33W3AlUAYuFFVp9e3r/0lKFWlrCpMwe5K8osr2Fpzc55vdAdN3VUW3Gu7HlnJDOqewZDu6e59BgO7pdcko2plwTLmbZ3H7C2z+WTjJ6wrXkdGIIMLB1/ID4ZcQlbxZtiyEDYvgLWfOt2BATr1hxEXwYgLnccHQFUp+3IO2x99lLK5c0mbMIGeDz+ENyPjgF6vrdj+5JMUPPgQgX796P6r20g75pj9zhlVn1AkxDvfvMPTXz/N6qLV+MTHuJxxHNvzWEZ2HcmQTkPwe1pxd+7yXbD8bVj+X2d6kKoSEA90G+ZMC5JzhDN9S+dBkNa5paNtMyxBxS5uCUpEvMBK4GRgIzAHuERVl0aVuRYYoarXiMjFwHdU9SIRGQa8BIwFcoAPgMGqGq69n2o9Bx2ml/zuRXaVBSkqd4bZ3+k+3neEaUjxe8nOSqZXxxRyO6W649Ol0aezM2iq3x+hPFhORbiCsqoSSsoLKSzbRmHpVvJLt7CmeB3flGxkbVk+IQ3jFw9HBLrwbdI5tbSU1F0bnd54uPtOyoLccdD/BBhwAnQdGnNtSYNBImVlhHftomr9Bqo2rKdyxUpKZ84kuHEj3g4d6HrTz+hw/vkH9EPcFpV+/jlbbr+D4ObN+Pvkkj7xOJIPOQR/r54EevXCm5WFJCfHNLuwqrJ8x3LeXfsu7619j40lGwFI8ibRP6s/uZm55Gbk0j21O1nJWXRM6kiHpA5kBjIJeAN7bp5A4jYVhqpgw2wnUW2aCxvnQWXRnvUpHaHzQMjq5fQOTO/q3Kd1gUA6BNLAnwqBVKdThj8FPD735rUmxCiWoGIXzwQ1HrhbVU91n98GoKr/F1VmulvmCxHxAflAV+DW6LLR5fa3v0PSUvQfQ/s6r0tNWmCf/xb7ebtSz8dQe50AHhSPglfBi+Kr2aFnzw0BjwfEu28kUZ+77mc5AKEQGty7dgfgSUsjdexYMk45hczTJh1QU1Zbp1VVFE+fTtGb/6Fs3jy0oo55o3w+PIGAc74u+kd0z1Su+yyLaJhQJEwwEiSsYUIaJqLh/X219tFwkpKof1uaRt21jdMBLe3oecstQcUonhfI9AQ2RD3fCNTu/1tTRlVDIlIEdHaXz6q17T4XAonIZGAywIAOKRQOSXOmiQbUnYnVLVjzn12in0f9UAhS84vgzODqwSeCTzz4xIvP6yfZEyDFm0SSNxmvzw/eAPiS3Fuy87xWDWbvH6NaPzl1/SDWjssjeNLSnFtGJoHevfD3zsXXtYvVlhoggQBZZ55J1plnolVVBDdvpmrTJoIbNxEp2U2kshKtrEIrKtBgVc12NX+0Rf+xEP3bvNdy53FEw1SFq6gMV1EVrqQqXEVVJEhEw4Q1QkQjRDRMJBIhQsT9vdeol9aol4te4jxKuNSg6oyyrmHnsUbcW9Tj6LI176COz7bxO2/SYs1uXksH0Hq06is4VfUJ4AlwzkGd8a/EulDXJA4JBAj07Uugb9+WDsW0dy8nRt24NYjnn+CbgOjZ93q5y+os4zbxZeF0lohlW2OMMW1YPBPUHGCQiPQTkQBwMTC1VpmpwOXu4/OBj9RpX5kKXCwiSSLSDxgEfBnHWI0xxiSYuDXxueeUrgem43Qz/4eqLhGRe4C5qjoVeBp4XkTygB04SQy33BRgKRACrquvB58xxpi2xy7UNcaYZmTdzGNn3cCMMcYkpDZTgxKR3cCKlo6jGXUBtrd0EM2oPb3f9vReof293yGq2r6He4lRq+5mXsuK9lRtFpG59n7bpvb0XqF9vt+WjqG1sCY+Y4wxCckSlDHGmITUlhLUEy0dQDOz99t2taf3CvZ+zX60mU4Sxhhj2pa2VIMyxhjThliCMsYYk5DaRIISkUkiskJE8kTk1paOpymJSG8R+VhElorIEhH5ibu8k4i8LyKr3Ps2Nce7iHhFZL6I/Nd93k9EZrvH+BV3fMc2QUQ6iMirIrJcRJaJyPi2fHxF5Gfud3mxiLwkIslt6fiKyD9EZJuILI5aVufxFMef3fe9SESOaLnIE0+rT1DuzL2PAacBw4BL3Bl524oQ8HNVHQYcBVznvr9bgQ9VdRDwofu8LfkJsCzq+f3Aw6o6ENgJXNkiUcXHn4B3VXUoMBLnfbfJ4ysiPYEbgdGqehjOOJ0X07aO7zPApFrL9nc8T8MZDHsQztx2jzdTjK1Cq09QONPC56nqGlWtAl4Gzm7hmJqMqm5R1a/cx7txfrx64rzHZ91izwLntEyETU9EegHfBp5ynwtwIvCqW6TNvF8RyQIm4gycjKpWqeou2vDxxRkgIMWdYicV2EIbOr6qOgNn8Oto+zueZwPPqWMW0EFEejRPpImvLSSoumbu3Wf23bZARPoCo4DZQHdV3eKuyge6t1BY8fAI8EugelrWzsAuVQ25z9vSMe4HFAD/dJs0nxKRNNro8VXVTcADwHqcxFSEM8dsWz2+1fZ3PNvN79eBaAsJql0QkXTgNeCnqlocvc6dQ6tNXC8gImcA21S1vUyM7QOOAB5X1VFAKbWa89rY8e2IU2voB+QAaezbHNamtaXjGW9tIUG1+dl3RcSPk5xeVNXX3cVbq5sC3PttLRVfEzsGOEtE1uI0156Ic46mg9skBG3rGG8ENqrqbPf5qzgJq60e328B36hqgaoGgddxjnlbPb7V9nc82/zv18FoCwkqlpl7Wy33/MvTwDJVfShqVfRsxJcD/2nu2OJBVW9T1V6q2hfnWH6kqpcCH+PMugxt6/3mAxtEZIi76CSciTrb5PHFado7SkRS3e929fttk8c3yv6O51TgMrc331FAUVRTYLvXJkaSEJHTcc5bVM/c+9sWDqnJiMgE4FPga/ack/kVznmoKUAusA64UFVrn5ht1UTkeOBmVT1DRPrj1Kg6AfOB76lqZUvG11RE5HCcDiEBYA3wA5w/Htvk8RWR3wAX4fRQnQ9chXPepU0cXxF5CTgeZxqRrcBdwJvUcTzdJP0XnGbOMuAHqmqjnbvaRIIyxhjT9rSFJj5jjDFtkCUoY4wxCckSlDHGmIRkCcoYY0xCsgRljDEmIVmCMqYeIvJTEUlt6TiMaY+sm7lp99xrUURVI3WsW4sz8vb2Zg/MmHbOalCmXRKRvu4cYs8Bi4GnRWSuO0/Rb9wyN+KMF/exiHzsLjtFRL4Qka9E5N/uGInGmDiwGpRpl9yR4dcAR6vqLBHp5F7Z78WZr+dGVV0UXYMSkS44Y8edpqqlInILkKSq97TQ2zCmTfM1XMSYNmudOwcPwIUiMhnn/0QPnMkvF9Uqf5S7fKbTKkgA+KKZYjWm3bEEZdqzUnCmkwduBsao6k4ReQZIrqO8AO+r6iXNF6Ix7ZedgzIGMnGSVZGIdMeZhrvabiDDfTwLOEZEBgKISJqIDG7WSI1pR6wGZdo9VV0oIvOB5Tizm86MWv0E8K6IbFbVE0TkCuAlEUly198OrGzWgI1pJ6yThDHGmIRkTXzGGGMSkiUoY4wxCckSlDHGmIRkCcoY8//t1bEAAAAAwCB/62nsKIlgSVAALAkKgCVBAbAUGCrgqBryXUcAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "from pyabc.visualization import plot_kde_1d\n", "\n", "fig, axes = plt.subplots(2)\n", "fig.set_size_inches((6, 6))\n", "axes = axes.flatten()\n", "axes[0].axvline(true_rate, color=\"black\", linestyle=\"dotted\")\n", "for m, ax in enumerate(axes):\n", " for t in range(0, history.n_populations, 2):\n", " df, w = history.get_distribution(m=m, t=t)\n", " if len(w) > 0: # Particles in a model might die out\n", " plot_kde_1d(\n", " df,\n", " w,\n", " \"rate\",\n", " ax=ax,\n", " label=f\"t={t}\",\n", " xmin=0,\n", " xmax=20 if m == 0 else 100,\n", " numx=200,\n", " )\n", " ax.set_title(f\"Model {m+1}\")\n", "axes[0].legend(title=\"Generation\", loc=\"upper left\", bbox_to_anchor=(1, 1))\n", "\n", "fig.tight_layout()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The true rate is closely approximated by the posterior over the rate of model 1.\n", "It is a little harder to interpret the posterior over model 2.\n", "Apparently a rate between 20 and 40 yields data most similar to the observed data.\n", "\n", "Lastly, we visualize the evolution of the population sizes. The population sizes were automatically selected by pyABC and varied over the course of the generations. (We do not plot the size of th first generation, which was set to 500)" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXcAAAEGCAYAAACevtWaAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3deXzU5bX48c/JvicEQlYgbIaEXQIqiIIKUYuCXKy0rtUrreLS3paq9/60XlsrvfS2tyrWorVq3XBF1FpQxJ0t7EsM+5IQkhBICNmX5/fHTGiACVlm+c5Mzvv18pWZ73xnvgckJ0+e7/OcI8YYlFJK+ZcAqwNQSinleprclVLKD2lyV0opP6TJXSml/JAmd6WU8kNBVgcA0KtXL5Oenm51GEop5VPWr19/1BiT4Og1r0ju6enp5ObmWh2GUkr5FBE50NZrOi2jlFJ+SJO7Ukr5IU3uSinlh7xizl0p1T00NDRQUFBAbW2t1aH4lLCwMNLS0ggODu7wezS5K6U8pqCggOjoaNLT0xERq8PxCcYYysrKKCgooH///h1+n88m9yUbC1mwLJ/D5TWkxIUzLyeDGaNTrQ5LKXUOtbW1mtg7SUTo2bMnpaWlnXqfTyb3JRsLeejdrdQ0NAFQWF7DQ+9uBdAEr5SX08TeeV35O/PJG6oLluWfSuwtahqaWLAs36KIlFLKu/hkcj9cXtOp40op5SpLlixhx44dp54/8sgjfPrpp22e//nnnzNt2jRPhHYan5yWSYkLp9BBIk+JC7cgGqWUu3jbvbXGxkaWLFnCtGnTyMrKAuCxxx6zLJ5z8cmR+7ycDMKDA087Fh4cwLycDIsiUkq5Wsu9tcLyGgz/ure2ZGOhU5+7f/9+hgwZwo033khmZiazZs2iurqaxx57jLFjxzJs2DDmzJlDS5e6SZMm8dOf/pTs7Gx+97vfsXTpUubNm8eoUaPYs2cPt912G2+//TYA69atY/z48YwcOZJx48ZRWVl52rWrqqq4/fbbGTduHKNHj+b9998HYPv27YwbN45Ro0YxYsQIdu3a5dSfEXx05N7yk3vBsvxTI/g5lwzQm6lK+ZD//mA7Ow6faPP1jQfLqW9qPu1YTUMTv3x7C6+vPejwPVkpMfzqmqHtXjs/P5+//vWvTJgwgdtvv51nnnmGe+65h0ceeQSAm2++mQ8//JBrrrkGgPr6+lP1r3bt2sW0adOYNWvWaZ9ZX1/PDTfcwOLFixk7diwnTpwgPPz02YTHH3+cyy67jBdeeIHy8nLGjRvHFVdcwbPPPsv999/PjTfeSH19PU1Np99T7AqfTO5gS/AzRqdyoraBMb/+hJqG5vbfpJTyGWcm9vaOd0afPn2YMGECADfddBNPPvkk/fv353/+53+orq7m2LFjDB069FRyv+GGG9r9zPz8fJKTkxk7diwAMTExZ52zfPlyli5dyu9//3vAtjT04MGDXHTRRTz++OMUFBQwc+ZMBg8e7PSf0WeTe4uYsGAuGtiLZduP8NBVQ3SZlVI+or0R9oT5nzm8t5YaF87iH1/k1LXPzBMiwt13301ubi59+vTh0UcfPW0XbWRkpFPXa2GM4Z133iEj4/Qp5MzMTC644AI++ugjrr76av7yl79w2WWXOXUtn5xzP9PUrEQOlFWzs/ik1aEopVzE8b21QJfcWzt48CCrVq0C4LXXXuPiiy8GoFevXpw8efLUHLoj0dHRZ82lA2RkZFBUVMS6desAqKyspLGx8bRzcnJyeOqpp07N52/cuBGAvXv3MmDAAO677z6mT5/Oli1bnP4z+k1yB1i+/YjFkSilXGXG6FSemDmc1LhwBNuI/YmZw11yby0jI4OFCxeSmZnJ8ePHueuuu7jzzjsZNmwYOTk5p6ZWHJk9ezYLFixg9OjR7Nmz59TxkJAQFi9ezL333svIkSOZMmXKWTV0Hn74YRoaGhgxYgRDhw7l4YcfBuDNN99k2LBhjBo1im3btnHLLbc4/WeUlp8gVsrOzjbONuu47plvaGhq5sN7J7ooKqWUq+Xl5ZGZmWlpDPv372fatGls27bN0jg6y9HfnYisN8ZkOzrfL0buADlDk9hWeMLhHJ1SSnU3fpPcW6ZmPtGpGaXUOaSnp/vcqL0r2k3uIvKCiJSIyLZWx64Xke0i0iwi2Wec/5CI7BaRfBHJcUfQjgxIiGJw7yiWbS/21CWVUl3gDVPBvqYrf2cdGbm/CFx5xrFtwEzgy9YHRSQLmA0Mtb/nGREJxEOmDk1k7f5jHK+q99QllVKdEBYWRllZmSb4Tmip5x4WFtap97W7zt0Y86WIpJ9xLA8clqGcDrxhjKkD9onIbmAcsKpTUXVRztAkFq7cw4rvSpg1Js0Tl1RKdUJaWhoFBQWdrk3e3bV0YuoMV29iSgVWt3peYD/mEcNTY0mKCWP59iOa3JXyQsHBwZ3qJqS6zrIbqiIyR0RyRSTXVT/FRYSpQxP5clcpNfXO12ZQSilf5erkXgj0afU8zX7sLMaYRcaYbGNMdkJCgssCyBmaRG1DM1/u0l/7lFLdl6uT+1JgtoiEikh/YDCw1sXXOKdx/eOJDQ9mmS6JVEp1Y+3OuYvI68AkoJeIFAC/Ao4BTwEJwEcisskYk2OM2S4ibwI7gEZgrjHGo/MjwYEBXD6kNyvySmhsaiYo0G+W8iulVId1ZLXMD9p46b02zn8ceNyZoJw1dWgi724sZO2+Y4wf1MvKUJRSyhJ+Oay95LwEQoMCWL5DNzQppbonv0zuESFBTBycwPLtR3SzhFKqW/LL5A6QMzSRwxW1bCtsu42XUkr5K79N7pdnJhIgsHyHrppRSnU/fpvc4yNDGNc/XpdEKqW6Jb9N7gBTs5LYWXySfUerrA5FKaU8yr+T+1Btv6eU6p78Ormn9YhgaEqMLolUSnU7fp3cwVZrZsPB45RU1rZ/slJK+Qm/T+5ThyZiDHy6o8TqUJRSymP8PrlnJEbTr2eErppRSnUrfp/cRYSpWYl8u+colbUNVoejlFIe4ffJHWzz7g1NhpX5WuNdKdU9dIvkPrpvD3pFheiSSKVUt+HqHqpeKTBAmJKVyAebi6hrbCI0KNDqkJQFlmwsZMGyfA6X15ASF868nAxmjPZYi1+lPKpbjNzBtlv1ZF0j3+4pszoUZYElGwt56N2tFJbXYIDC8hoeencrSzY67AKplM/rNsl9/KCeRIYEsny7bmjqjhYsy6em4fSmYDUNTSxYlm9RREq5V7dJ7qFBgUwa0ptPdhTT1Kw13rubw+U1nTqulK/rNskdbKtmjp6sY9Oh41aHojzoeFU9oUGO/6mnxIV7OBqlPKNbJfdJGQkEBwrLdGqm21i77xhXP/kV9U3NBAfKaa+FBwcyLyfDosiUcq92k7uIvCAiJSKyrdWxeBH5RER22b/2sB8XEXlSRHaLyBYROd+dwXdWTFgwFw3sxTJtv+f3mpoNT67YxexFqwgNCuD9uRezYNZI4iNDAEiICuWJmcN1tYzyWx0Zub8IXHnGsQeBFcaYwcAK+3OAq4DB9v/mAH92TZiukzM0kQNl1ewsPml1KMpNik/UcuPzq/nDJzu5dmQKH943keFpscwYncr7cycA8NMpgzWxK7/WbnI3xnwJHDvj8HTgJfvjl4AZrY6/bGxWA3EikuyqYF1hSmYiIlrj3V+t/K6Eq/70FZsPVbBg1gj+eMMookL/tZ0jrUc40WFB7DisvXWVf+vqnHuiMabI/vgIkGh/nAocanVegf3YWURkjojkikhuaannygL0jgljdJ84lmlvVb9S39jMbz7cwY9eXEfv6FA+uPdirs/ug8jp8+wiQmZyDHlFmtyVf3P6hqqxTV53egLbGLPIGJNtjMlOSEhwNoxOmTo0iW2FJyjUZXB+Yf/RKmY9+y3Pf72PWy7qx5K5ExjUO6rN87OSY/juSCXNuiRW+bGuJvfilukW+9eWYumFQJ9W56XZj3mVnKFJAHyiUzM+7/1NhUx76msOlFXz7E1jeGz6MMKCz11eIis5hur6Jg4cq/ZQlEp5XleT+1LgVvvjW4H3Wx2/xb5q5kKgotX0jdfo3yuSwb2jdEmkD6uub+SXb2/m/jc2MSQpmn/cP5ErhyV16L1ZKTEAOu+u/FpHlkK+DqwCMkSkQETuAOYDU0RkF3CF/TnAP4C9wG7gOeBut0TtAjlDk1i7/xjHq+qtDkV1Ul7RCa556mveWl/APZMH8cacC0ntxGakQb2jCAwQnXdXfq3dqpDGmB+08dLlDs41wFxng/KEqUMTeXrlblZ8V8KsMWlWh6M6wBjDK6sP8OuP8ogND+bVOy5g/KBenf6csOBABiVEsUOTu/Jj3WqHamvDU2NJjg3TJZE+oqK6gbte2cDD72/nogE9+fj+iV1K7C0yk6N15K78WrdN7i3t977cVUpNfVP7b1CWWX/AVkLg07xi/vPqIfzttrH0igp16jOzUmIoqqjlmE7LKT/VbZM72Obdaxua+XKXtt/zRk3NhoUrd/P9v6wmIADevms8cy4ZSECAtP/mdmQm226q6uhd+atu0YmpLWP7xxMbHsyy7UdOLY9U1mndKSkxJozosEB2lVQxbUQyv505nJiwYJddq3Vyn+DE9I5S3qpbJ/fgwAAuz+zNirwSGpuaCQrs1r/IWKqlU1JLQ40jJ2o5cgJuGJvG/Jkjztpp6qxeUaH0jg7V5ZDKb3X7bDY1K4mKmgbW7juzfI7yJEedkgC+3lXm8sTeIislRlfMKL/V7ZP7peclEBYcwPIduqHJSlZ0SspKjmF3yUnqGvWGuvI/3T65h4cEMnFwAsu1xrul2uqI5M5OSZnJMTQ2G3aXaPln5X+6fXIH26qZwxW1bCvUX9GtMi8nw+OdkrQMgfJnmtyBy4f0JjBAWK5lgC0zfVQKvSJDCA4UBEiNC3d7p6T0npGEBQeQV1TptmsoZZVuvVqmRY/IEMalx7Ns+xF+PlV7alphS0EFRSfq+M2MYdx0YT+PXDMwQBiSFMOOogqPXE8pT9KRu93UoYnsLD7JvqNVVofSLS3OPURYcADXjkrx6HVtjTsq9X6L8jua3O2m2jcxaa0Zz6uub2TppsNcPTzZpRuVOiIrJYaKmgYOV9R69LpKuZsmd7vUuHCGpcbokkgL/GPrEU7WNXJDdp/2T3axrORoQG+qKv+jyb2VnKwkNhw8Tkmlb4/ilmwsZML8z+j/4EdMmP8ZSzZ6XTOs0yxed5D+vSIZ1z/e49fOSIpBRGvMKP+jyb2VqUOTMAY+3VHS/sleqmUbf2F5DQYoLK/hoXe3em2C311yknX7j3PD2LObWXtCVGgQ/eIjdOSu/I4m91bOS4yiX88IlvnwvLujbfw1DU0sWJZvUUTn9lbuIQIDhJnnu2/JY3uyUmLIO6LJXfkXTe6tiAg5Q5P4ds9RKmsbrA6nS6zYxt9V9Y3NvLOhgMuH9KZ3dJhlcWQlx3CgrNpn/58r5Ygm9zNMzUqkocmwMt83a7zHhDveuuDObfxd9dl3xRw9Wc/scZ6/kdpaS/nf/CO6mUn5D03uZxjdtwe9okJ9cknk17uOcqKmkTN7Wbh7G39XLV53iMSYUC4ZnGBpHKfKEOhNVeVHnEruInK/iGwTke0i8lP7sXgR+UREdtm/9nBNqJ4RGCBMyerN5/mlPlUtcP/RKua+toHBiVH89rrhpLYaqf8i5zy3buPviqKKGr7YWcr1Y/pYXkc/KSaMuIhgXTGj/EqXv6tEZBhwJzAOGAlME5FBwIPACmPMYGCF/blPmTo0iZN1jXy7p8zqUDqksraBO1/ORQSev2Uss8f15ZsHL2PVQ5cRHCheuev27dwCmg1834K17WcSEbKSY3TFjPIrzgyZMoE1xphqY0wj8AUwE5gOvGQ/5yVghnMhet74gT2JCg1i+Xbv39DU3Gz42eJN7D1axTM/PJ++PSNOvZYcG86sMWm8mVtAyQnvWbvf3GxYnHuI8QN7nhavlTKTY/juSCWNTc1Wh6KUSziT3LcBE0Wkp4hEAFcDfYBEY0yR/ZwjQKKjN4vIHBHJFZHc0lLvunkZGhTIoN6Rts01Xr4R6H8/yefTvBIe/l4m4x30Ar3r0kE0NRue+2qvBdE59u2eMgqO13DDWOtH7S2ykmOoa2xmf5n3/ZajVFd0ObkbY/KA3wHLgX8Cm4CmM84xgMOKTMaYRcaYbGNMdkKCtTfUzrRkYyHbD5+g2eDVG4E+2HyYhSv3MHtsH24dn+7wnL49I7h2ZAqvrD7Isap6zwbYhsW5h4gND/aqpuQtK2Z2aPlf5SecupNljPmrMWaMMeYS4DiwEygWkWQA+1ef2+65YFk+DU2n/0zyto1A2wormPf2ZrL79eCx6cPOubvz7kkDqWlo4m/f7PNghI4dr6pn2bYjXDc6lbDgQKvDOWVQ7yiCA0Xn3ZXfcHa1TG/7177Y5ttfA5YCt9pPuRV435lrWMHbNwKVVtYx5+Vc4iNC+PNNYwgJOvf/xsGJ0Vw5NIkXv93PCYs36izZVEh9U7NXTckAhAQFMLh3tC6HVH7D2TVo74jIDuADYK4xphyYD0wRkV3AFfbnPqWtDT8J0aEejuRsdY1N3PXKeo5V17PoluwOxzR38iAqaxv5+6oDbo6wbcYY3lh7iBFpsaemQbyJrba7JnflH5ydlplojMkyxow0xqywHyszxlxujBlsjLnCGHPMNaF6zrycDMIdTBlUVNfz5U7rbv4aY/jV+9vJPXCcBbNGMiw1tsPvHZ4Wy6XnJfDC1/uoqbdm/f7mggryiyu9btTeIislhtLKOkor66wORSmn6Q5VB2aMTuWJmbaNQC39PB+9NosBvaP50YvreHWNNaPfl1cd4I11h5g7eSDXjOx8x6J7LhtEWVU9r6896Ibo2rd43SHCgwO5tguxe0Kmvba7jt6VP9Aeqm2YMTr1rF2ds8b04b7XN/Jf721j/9EqHroqk4Az9/q7ybe7j/LYhzu4IrM3P5/StVICY9PjGdc/nkVf7uXGC/sSGuS5G5pVdY0s3VTI1cOTifZwt6WOykr+VxmCS87zrhVcSnWWjtw7ISo0iEU3j+G28ek899U+7np1vUemOA6WVXP3axsY0CuSP94wyqkfKPdMHsSRE7W8u8Gzyzo/2lpEVX2T5UXCziUuIoSU2DAduSu/oMm9k4ICA3j02qH86posPtlRzA2LVrl19+fJukb+/eV1GAPP3ZLt9Kh34uBejEiL5c+f7/Hobsw31x1iQEIk2f28u9RQVoqWIVD+QZN7F/1oQn+euyWb3SUnmbHwG75zQ7OHltICe0qrWPjD80nvFen0Z4oIcycP4uCxaj7cUtT+G1xgd0kluQeOc0O2Nd2WOiMzOYa9R6uobfCdonFKOaLJ3QmXZyby5o8voskYZv15FZ/nu3a/1v99upNPdhTzX1dncvHgs0sLdNWUzETOS4xi4crdNDc73EDsUm/mFhAUIMw8P83t13JWVnIMTc2GncW6U1X5Nk3uThqWGsuSuRPoGx/BHS/l8vfVrllJ89GWIp78bDfXj0njRxPSXfKZLQICbKP3XSUnWb7DvcXR6hubeWd9AZdn9vaKfQLtaVl/r/PuytdpcneB5Nhw3vrJRVx6XgIPL9nGbz7cQZMTI+Lthyv4xVubOb9vHL+57tylBbrqe8OTSe8ZwcKVu7GVAHKPFXnFlFXVM3tsX7ddw5X6xkcQGRKo8+7K52lyd5HI0CCeuyWb28an8/zX+7jrlfVU1zd2+nOOnqxjzsvriYsI5tmbx7htuWJQYAB3TRrI1sIKvtx11C3XAFuRsKSYMJ9ZWhgQIPadqjoto3ybJncXCgwQHr12KI9ek8WnecXc8JfVnVpJU9/YzN2vbODoyToW3Zzt9qbR141OIyU2jIWf7XbL5x8ut3dbyk4j0EP7AVwhMzmGHUUnPHI/Qil30eTuBrfZV9LsKbWtpOno/O2jH2xn7f5j/M+sEQxP63hpga4KCQpgziUDWLv/GGv2ur7r1Fu5BRgv6bbUGVkpMZysa6TguHcUilOqKzS5u0nrlTTXP9v+Spq/rz7Aa2sO8pNLBzJ9lOf6nc4e15deUSE8vdK1o/fmZsObuYe4eFAv+sR7R7eljspM1obZyvdpcnejYamxvD/3YvrGR3D7i+vaXEmzak8Z/710O5cN6c28nK6VFuiqsOBA7rh4AF/tOsrmQ+Uu+9xv9hylsLyG73tpkbBzyUiMJkA0uSvfpsndzZJiw3jrJxcxOaM3Dy/Zxq/PWElz6Fg1d7+6nn49I/i/2aMsmZu+6cK+xIQFsdCFo/fF6w4RFxHM1CyHXRa9WnhIIP17RepySOXTtHCYB0SGBrHolmx+/eEO/vr1PlbvLeNYVT1HKmoJDBCCAuD5uycQY1FBreiwYG6b0J8nV+wi/0glGUnRTn3esap6lm8v5sYL+3pVt6XOyEqJZcOB41aHoVSX6cjdQ1pW0sw8P5Xth09QVFGLARqbDc1GXDol0hU/Gp9OREggz3zu/Oj9vY3e2W2pMzKToyksr6GixtrOVUp1lSZ3D1uz9+zeJfVNzZb3Z+0RGcJNF/bjg82H2X+0qsufY4zhzXWHGNknjiFJ3tdtqaOydKeq8nGa3D3Mm/uz/vvF/QkKDODZL/Z0+TM2HSq3dVvyseWPZ8pK0eSufJsmdw9rqz9rW8c9qXdMGDdk9+GdDQVd/mHT0m3pmpHJLo7Os3pHh9ErKkTLECif5VRyF5Gfich2EdkmIq+LSJiI9BeRNSKyW0QWi0iIq4L1B476s4YHB3p8CWRbfnzpAIyBRV/u7fR7q+oa+WDzYaaN8N5uS52RmRxDnhtKOSvlCV1O7iKSCtwHZBtjhgGBwGzgd8AfjTGDgOPAHa4I1F846s/6xMzhZ7X0s0pajwiuG53KG+sOcvRk5xpFf7TF1m3Jl2+ktpaVHMPOIydp8GBTE6VcxdmlkEFAuIg0ABFAEXAZ8EP76y8BjwJ/dvI6fsVRf1ZvctekgbyzoYC/fr2PB64c0uH3Lc49xMCESMZ4ebeljspKiaG+qZk9pSd9+uaw6p66PHI3xhQCvwcOYkvqFcB6oNwY01IOsQDw3iymHBqQEMXVw5P5+6oDVFR3bCngruJK1h84zuyxfb2+21JHaW135cucmZbpAUwH+gMpQCRwZSfeP0dEckUkt7S0tKthKDeZO3kQJ+saeWnV/g6dv3jdIYIChOvO95+f5QN6RRISFKA3VZVPcuaG6hXAPmNMqTGmAXgXmADEiUjLdE8aUOjozcaYRcaYbGNMdkKCb9T67k4yk2O4IrM3L3yzj6q6c9elr29s5t2NhUzJSqRXlPd3W+qooMAAMhKjtba78knOJPeDwIUiEiG238MvB3YAK4FZ9nNuBd53LkRllbmTB1Fe3cCra87dOvDTvGKOVdX7zY3U1rLstd3d2a1KKXdwZs59DfA2sAHYav+sRcADwH+IyG6gJ/BXF8SpLDC6bw8mDOrJc1/to7ahqc3z3lh3iJTYMCYO9r/fwDKTozlWVU9JZedWDillNafWuRtjfmWMGWKMGWaMudkYU2eM2WuMGWeMGWSMud4Yo98VPmzu5EGUVtbxVu4hh68Xltfw1a5SZmX38aluSx2VlWJrmqLz7srX6A5VdU4XDejJ+X3jePaLvQ7Xe7ck/evHpHk6NI8YkmyrkKm13ZWv0eSuzklEuOeyQRSW17Bk4+n3xpuaDW/lFvhkt6WOigkLpk98uCZ35XM0uat2Tc7oTVZyDH/+fM9pjUa+2W3rtuSPN1Jby0qO0bXuyudoclftEhHmTh7E3qNVfLyt6NTxxesO0SMimCk+2G2pMzKTY9h3tIrq+nMvCVXKm2hyVx1y5bAkBiREsnDlHowxtm5LO45w3eg0QoN8s9tSR2Ulx2AM5B/R9e7Kd2hyVx0SGCDMnTSIvKITfPZdCe9uKKChyfj9lAz8qwyBzrsrX6LJXXXYtaNS6BERzE9eWc9vPsojOFC6xVx0Wo9wosOCdDmk8ima3FWHfbSliJN1jTQ02W6qNjQZHnp361mraPyNiNhqu3eDH2TKf2hyVx22YFn+qcTeoqahyfL+r56QlRzDd0cqaW7WMgTKN2hyVx3mzf1f3S0rJYbq+iYOHKu2OhSlOkSTu+owb+7/6m5ZLTdVdd5d+QhN7qrDvL3/qzsN6h1FUED3uIGs/IOzbfZUN9LSGnDBsnwOl9eQEhfOvJwMr24Z6CphwYEMTIjS5ZDKZ2hyV53i7f1f3SkrJYbVe8usDkOpDtFpGaU6KDM5mqKKWo5X1VsdilLt0uSuVAdlJdtqu+u8u/IFmtyV6qBMre2ufIgmd6U6qGdUKIkxobocUvkETe5KdUKmvWG2Ut5Ok7tSnZCVHMPukpPUNbbdMFwpb9Dl5C4iGSKyqdV/J0TkpyISLyKfiMgu+9cergxYKStlpcTQ2GzYXXLS6lCUOqcuJ3djTL4xZpQxZhQwBqgG3gMeBFYYYwYDK+zPlfILmVqGQPkIV03LXA7sMcYcAKYDL9mPvwTMcNE1lLJces9IwoMDySvSrkzKu7kquc8GXrc/TjTGtDTaPAI4bLApInNEJFdEcktLS10UhlLuFRggZCRFs6OowupQlDonp5O7iIQA1wJvnfmaMcYADgtgG2MWGWOyjTHZCQkJzoahlMdkpcSQV1SJ7Z+3Ut7JFSP3q4ANxphi+/NiEUkGsH8tccE1lPIamckxVNQ0cLii1upQlGqTK5L7D/jXlAzAUuBW++NbgfddcA2lvEZLbfc8vamqvJhTyV1EIoEpwLutDs8HpojILuAK+3Ol/MaQpGhEtAyB8m5Olfw1xlQBPc84VoZt9YxSfikyNIj0npFaQEx5Nd2hqlQXZGkZAuXlNLkr1QWZydEcKKumsrbB6lCUckiTu1JdkJViu6maf0Q3MynvpMldqS44VYZAp2aUl9LkrlQXJMWE0SMiWG+qKq+lyV2pLhARW213XeuuvJQmd6W6KCs5hu+OVNLY1Gx1KEqdRZO7Ul2UmRxDXWMz+8uqrA5FqbNocleqi1pWzOzQ8r/KC2lyV6qLBiZEERwoOu+uvJImd6W6KCQogMG9o3XFjPJKmtyVckJWipYhUN5Jk7tSTlFCpDkAABLLSURBVMhMjqG0so7SyjqrQ1HqNJrclXLCqdruOnpXXkaTu1JOyNIyBMpLaXJXygmxEcGkxoXryF15HU3uSjkpMzlal0Mqr+NUJyallG1qZmV+KbUNTYQFB1odTqct2VjIgmX5HC6vISUunHk5GcwYnWp1WMpJOnJXykmZyTE0NRt2FvveTtUlGwt56N2tFJbXYIDC8hoeencrSzYWWh2acpImd6Wc1FKGwBfn3Rcsy6emoem0YzUNTSxYlm9RRMpVnEruIhInIm+LyHcikiciF4lIvIh8IiK77F97uCpYpbxRnx4RRIUG+dy8+4GyKgrLaxy+driN48p3ODty/xPwT2PMEGAkkAc8CKwwxgwGVtifK+W3AgKEIUnR5PlIAbE9pSf5jzc3cdn/ftHmOSlx4R6MSLlDl2+oikgscAlwG4Axph6oF5HpwCT7aS8BnwMPOBOkUt4uKyWG9zYUYoxBRKwOx6GdxZU8/dluPtxymJCgAG4bn07f+HDmf3z61ExwoDAvJ8PCSJUrOLNapj9QCvxNREYC64H7gURjTJH9nCNAoqM3i8gcYA5A3759nQhDKetlJsfwct0BCo7X0Cc+wupwTrP9cAVPf7abj7cdISIkkDsvGcCdEwfQKyoUgNjwkFOrZYIChYiQQK4clmRx1MpZziT3IOB84F5jzBoR+RNnTMEYY4yIGEdvNsYsAhYBZGdnOzxHKV/RslN1++ETXpPctxSU8+SK3XyaV0x0aBD3TB7E7Rf3Jz4y5LTzZoxOPbX0cfXeMmYvWs3fvtnPXZMGWhG2chFnknsBUGCMWWN//ja25F4sIsnGmCIRSQZKnA1SKW+XkRRNgNjKEFg96l1/4DhPfbaLz/NLiQ0P5mdXnMdtE9KJDQ9u970XDujJ5UN688znu5k9tg89zvhBoHxHl2+oGmOOAIdEpGVy7nJgB7AUuNV+7FbgfaciVMoH/HPbEQJEeHLFLibM/8ySdeKr95Zx4/Or+bc/f8uWggrm5WTw9QOTuf+KwR1K7C0euGoIVXWNPL1ytxujVe7m7A7Ve4FXRSQE2Av8CNsPjDdF5A7gAPB9J6+hlFdr2QjU2GybXbRtBNoC4PadnsYYvtldxpOf7WLtvmP0igrlv67O5MYL+xIR0rVv7/MSo7l+TB9eXrWf28ane800k+ocMcb66e7s7GyTm5trdRhKdcmE+Z85XC8eFCBMHNyL1B7hpPWIIDUu3P44nF6RoQQEdHxVzZklAn4x9TziIkN4csUuNh4sJykmjB9fOoAfjOvrkhIIRypqmfT7leQMTeJPs0c7/XnKPURkvTEm29FrWltGKSe1teGnsdlQUlnHhoPlVNQ0nPZaSFAAqXG2RJ8aZ/svLT6c1LgIUnuEkxQTRqA9+bf8ZtCyXLGwvIb/eGszxkBqXDi/mTGM67PTCA1yXV2bpNgw7ri4PwtX7uHfLx7A8LRYl3228gxN7ko5KSUu3OHIPTUunI/umwjAybpGCo/XUHC8msLyGvvjGgrKa8jLK+HoydM7OQUFCEmxYaTGhbOloOKsEgHGQFx4MCt/MYmQIPdUEfnxpQN5fe0hnvg4j1f//QKvXb+vHNPkrpST5uVknDayBggPDjxtI1BUaBAZSdFkJEU7/IzahqbTkn5hefWpx2cm9hYVNQ1uS+wAMWHB3HvZIP77gx18sbOUSRm93XYt5Xqa3JVyUstNU2fK5oYFBzIwIYqBCVFnvdbWnL4nSgTceEE//vbNfuZ//B0TByecmipS3k+Tu1Iu0HojkKt15DcDdwkJCmBeTgb3vr6R9zYWMmtMmtuvqVxDS/4q5eVmjE7liZnDSY0LR7DN5T8xc7jHGmp8b3gyI9Ni+cPyfGrbmCJS3kdH7kr5AHf+ZtCegADhwasy+cFzq3nx2/385FItS+ALdOSulGrXRQN7ctmQ3ixcuZvjVfVWh6M6QJO7UqpDHrjSVpZgoZYl8Ama3JVSHZKRFM2sMWm8vOoAh45VWx2Oaocmd6VUh/1synmIwP8u1x6r3k6Tu1Kqw5Jjw7nj4v4s2XSYbYUVVofjVks2FjJh/mf0f/Ajyyp9OkOTu1KqU34yaSA9IoKZ//F3VofiNi31fArLazC0VPrc6lMJXpO7UqpTbGUJBvP17qN8ubPU6nBc4kRtA3lFJ/h0RzEvfbuf/7dk61llH2oamvjtP/Lwhkq6HaElf5VSnVbX2MQVf/iCqNBgPrr34k6VL3alM0shOyr7YIzheHUDhfaaPQWn6vfYvx6v5kRtY4ev2Ts6lHH947mgfzwXDOjJ4N5RlhVV05K/SimXCg0KZF7OEO57fSNLNhUy83zPlyVwVAp53tub+TSvmLiI4FZF2Gqorj99FB4ZEniqzn52vx6n6uy31NyfsfAbDpfXnnXNuPBgLhzQkzX7yvhwSxEA8ZEhjE3vwQX9e3LBgHiGJMV4RQ0eHbkrpbqkudkwfeE3HKuqZ8XPL3VJk5DOaKugGkBcRPC/6uT3sNXIb6mfn9YjnNjw4HOOts/8wQG2ej4tZR+MMRw8Vs2afcdYs/cYa/aVUXDcFktMWBBj0+O5YEA84/r3ZFhKDEGB7pkB15G7UsrlAgKEh64ewg+fW8PLq/Yz5xLPlSUwxrSZ2AXY9MhUpz6/vUqfIkK/npH06xnJ97P7ALbfHNbuK2PN3mOs3XeMFd+VALbfEsak26dx+sczIi2OkKCADk0pOUNH7kopp9z2t7VsOHCcL385mbiIELdfr6KmgQff2cLH2444fD01LpxvHrzM7XG0p+RELWv22RL9mn1l7Cw+CUBYcABpceHsL6s+1XcXTv/NoKPONXJ3KrmLyH6gEmgCGo0x2SISDywG0oH9wPeNMcfP9Tma3JXyXXlFJ7j6ya+4c+IA/vPqTLdea/2B49z3+kaKT9Ry1bAkPs0rpqah+dTrXUmQnlJ2so51+4+zZl8Zf1914LTE3qKzP5jOldxdMRE02RgzqtUFHgRWGGMGAyvsz5VSfiozOYZ/Oz+NF7/dT8Fx95QlaG42PPP5br7/l1WIwFs/uYinfng+T8wcYVkp5M7qGRXKlcOS+NU1Q2lykNih7X68XeGOOffpwCT745eAz4EH3HAdpZSX+I8p5/HB5sP8YflO/nDDKJd+dkllLT9/czNf7TrK94Yn89uZw4kNDwasLYXsjLb67rqyu5azI3cDLBeR9SIyx34s0RhTZH98BEh09EYRmSMiuSKSW1rqHxshlOquUuLC+dGE/ry3qZDth11XluDLnaVc/aevWLvvGE/MHM7TPxx9KrH7snk5GYSfsbrI1d21nE3uFxtjzgeuAuaKyCWtXzS2CX2Hv38YYxYZY7KNMdkJCQlOhqGUstpdkwYSG+6asgQNTc3M//g7bnlhLfGRISy952J+MK6vZZuFXM0T3bWcmpYxxhTav5aIyHvAOKBYRJKNMUUikgyUuCBOpZSXiw0P5p7Jg/jNR3l8tauUiYO7Nmg7dKya+97YyMaD5fxgXF8emZZFeIhn19B7grunlLo8cheRSBGJbnkMTAW2AUuBW+2n3Qq872yQSinfcPNF/UjrEc78j7+juY2bhufyj61FXP3kV+wuPsnTPxzNEzOH+2Vi9wRnpmUSga9FZDOwFvjIGPNPYD4wRUR2AVfYnyulugFbWYIMth8+wdLNhzv8vtqGJv7zva3c/eoGBiRE8Y/7JzJtRIobI/V/XZ6WMcbsBUY6OF4GXO5MUEop33XNiBSe+2ovC5blc+WwpHbLEuwsruTe1zaSX1zJjy8dwC+mZhDspu363Yn+DSqlXCogQHjoqkwKy2t4ZfWBNs8zxvDG2oNc+/TXlFXV8dLt43joqkxN7C6if4tKKZebMKgXl5yXwFOf7aaiuuGs10/UNnDv6xt58N2tZPeL5x/3T+TS83TVnCtpcldKucWDVw7hRG0Dz3yx+7Tjmw6V870nv+LjbUeYl5PBy7ePo3d0mEVR+i+tCqmUcouslBiuG53K81/uZcnGQkpO1BEdFkRlbSMpceG8+eMLGdMv3uow/ZaO3JVSbjM8NZYmA8Un6jDAidpGRGDu5IGa2N1Mk7tSym2e/2rfWceaDSxcuceCaLoXTe5KKbdpq8qhK6sfKsc0uSul3KatKoeurH6oHNPkrpRyG09UP1SO6WoZpZTbtNeLVLmPJnellFv5akMNX6fTMkop5Yc0uSullB/S5K6UUn5Ik7tSSvkhTe5KKeWHxNbD2uIgREqBtgs/W6sXcNTqILrIV2P31bhBY7dKd429nzHGYa1kr0ju3kxEco0x2VbH0RW+Gruvxg0au1U09rPptIxSSvkhTe5KKeWHNLm3b5HVATjBV2P31bhBY7eKxn4GnXNXSik/pCN3pZTyQ5rclVLKD2lyd0BE+ojIShHZISLbReR+q2PqLBEJFJGNIvKh1bF0hojEicjbIvKdiOSJyEVWx9RRIvIz+7+XbSLyuoiEWR1TW0TkBREpEZFtrY7Fi8gnIrLL/rWHlTG2pY3YF9j/zWwRkfdEJM7KGNviKPZWr/1cRIyI9HLFtTS5O9YI/NwYkwVcCMwVkSyLY+qs+4E8q4Pogj8B/zTGDAFG4iN/BhFJBe4Dso0xw4BAYLa1UZ3Ti8CVZxx7EFhhjBkMrLA/90YvcnbsnwDDjDEjgJ3AQ54OqoNe5OzYEZE+wFTgoKsupMndAWNMkTFmg/1xJbYE4zMFqUUkDfge8LzVsXSGiMQClwB/BTDG1Btjyq2NqlOCgHARCQIigMMWx9MmY8yXwLEzDk8HXrI/fgmY4dGgOshR7MaY5caYRvvT1UCaxwPrgDb+3gH+CPwScNkKF03u7RCRdGA0sMbaSDrl/7D9Q2m2OpBO6g+UAn+zTyk9LyKRVgfVEcaYQuD32EZeRUCFMWa5tVF1WqIxpsj++AiQaGUwTrgd+NjqIDpKRKYDhcaYza78XE3u5yAiUcA7wE+NMSesjqcjRGQaUGKMWW91LF0QBJwP/NkYMxqownunBk5jn5+eju0HVAoQKSI3WRtV1xnbGmmfWyctIv+FbVr1Vatj6QgRiQD+E3jE1Z+tyb0NIhKMLbG/aox51+p4OmECcK2I7AfeAC4TkVesDanDCoACY0zLb0lvY0v2vuAKYJ8xptQY0wC8C4y3OKbOKhaRZAD71xKL4+kUEbkNmAbcaHxnA89AbAOCzfbv2TRgg4gkOfvBmtwdEBHBNu+bZ4z5g9XxdIYx5iFjTJoxJh3bDb3PjDE+MYI0xhwBDolIhv3Q5cAOC0PqjIPAhSISYf/3czk+cjO4laXArfbHtwLvWxhLp4jIldimIq81xlRbHU9HGWO2GmN6G2PS7d+zBcD59u8Fp2hyd2wCcDO2Ue8m+39XWx1UN3Ev8KqIbAFGAb+1OJ4Osf+28TawAdiK7XvLa7fEi8jrwCogQ0QKROQOYD4wRUR2YftNZL6VMbaljdifBqKBT+zfr89aGmQb2ojdPdfynd9elFJKdZSO3JVSyg9pcldKKT+kyV0ppfyQJnellPJDmtyVUsoPaXJXPklEEkXkNRHZKyLrRWSViFxnUSyTRGR8q+c/EZFbrIhFqRZBVgegVGfZNwktAV4yxvzQfqwfcK0brxnUqjDVmSYBJ4FvAYwxXrnGWnUvus5d+RwRuRx4xBhzqYPXArFtvpkEhAILjTF/EZFJwKPAUWAYsB64yRhjRGQM8Acgyv76bcaYIhH5HNgEXAy8jq2U7P8DQoAy4EYgHFsVwiZsRc/uxbY79aQx5vciMgp4FluVyD3A7caY4/bPXgNMBuKAO4wxX7nub0l1dzoto3zRUGw7QR25A1tFxrHAWOBOEelvf2008FMgCxgATLDXEHoKmGWMGQO8ADze6vNCjDHZxpj/Bb4GLrQXNXsD+KUxZj+25P1HY8woBwn6ZeABe53xrcCvWr0WZIwZZ4/pVyjlQjoto3yeiCzENrquBw4AI0Rklv3lWGCw/bW1xpgC+3s2AelAObaR/Ce22R4CsZXsbbG41eM0YLG9qFYIsK+duGKBOGPMF/ZDLwFvtTqlpSDdenssSrmMJnfli7YD/9byxBgz196aLBdbAa97jTHLWr/BPi1T1+pQE7Z//wJsN8a01c6vqtXjp4A/GGOWtprmcUZLPC2xKOUyOi2jfNFnQJiI3NXqWIT96zLgLvt0CyJyXjsNP/KBhJZerSISLCJD2zg3Fii0P7611fFKbEWrTmOMqQCOi8hE+6GbgS/OPE8pd9DRgvI59pugM4A/isgvsd3IrAIewDbtkY6tJrbYX2uzXZwxpt4+hfOkfRolCFsnq+0OTn8UeEtEjmP7AdMyl/8B8La9o869Z7znVuBZe1OGvcCPOv8nVqrzdLWMUkr5IZ2WUUopP6TJXSml/JAmd6WU8kOa3JVSyg9pcldKKT+kyV0ppfyQJnellPJD/x80V+rNRbitCAAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "populations = history.get_all_populations()\n", "ax = populations[populations.t >= 1].plot(\"t\", \"particles\", style=\"o-\")\n", "ax.set_xlabel(\"Generation\");" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The initially chosen population size was adapted to the desired target accuracy.\n", "A larger population size was automatically selected by pyABC while both models were still alive. The population size decreased during the later populations thereby saving computational time." ] } ], "metadata": { "celltoolbar": "Raw Cell Format", "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.8" } }, "nbformat": 4, "nbformat_minor": 4 }