diff --git a/docs/source/toy/index.rst b/docs/source/toy/index.rst index 3503ec66d..69a63255b 100644 --- a/docs/source/toy/index.rst +++ b/docs/source/toy/index.rst @@ -38,4 +38,5 @@ Some toy classes provide extra functionality defined in the simple_harmonic_oscillator_model sir_model stochastic_degradation_model + stochastic_logistic_model twisted_gaussian_logpdf diff --git a/docs/source/toy/stochastic_logistic_model.rst b/docs/source/toy/stochastic_logistic_model.rst new file mode 100644 index 000000000..a5fb8eec2 --- /dev/null +++ b/docs/source/toy/stochastic_logistic_model.rst @@ -0,0 +1,7 @@ +************************* +Stochastic Logistic Model +************************* + +.. currentmodule:: pints.toy + +.. autoclass:: StochasticLogisticModel diff --git a/examples/README.md b/examples/README.md index e4a7a9415..9b668c988 100644 --- a/examples/README.md +++ b/examples/README.md @@ -117,6 +117,7 @@ relevant code. - [Simple Harmonic Oscillator model](./toy/model-simple-harmonic-oscillator.ipynb) - [SIR Epidemiology model](./toy/model-sir.ipynb) - [Stochastic Degradation model](./toy/model-stochastic-degradation.ipynb) +- [Stochastic Logistic model](./toy/model-stochastic-logistic-growth.ipynb) ### Distributions - [Annulus](./toy/distribution-annulus.ipynb) diff --git a/examples/toy/model-stochastic-logistic-growth.ipynb b/examples/toy/model-stochastic-logistic-growth.ipynb new file mode 100644 index 000000000..1ccd30e0e --- /dev/null +++ b/examples/toy/model-stochastic-logistic-growth.ipynb @@ -0,0 +1,177 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Stochastic Logistic Growth model\n", + "\n", + "This example shows how the stochastic logistic growth model can be used.\n", + "This model describes the growth of a population of individuals, where the birth rate per capita, initially $b_0$, decreases to 0 as the population size, $\\mathcal{C}(t)$, approaches a carrying capacity, $k$.\n", + "\n", + "The population grows starting from an initial population size, $n_0$, to a carrying capacity $k$ following a rate according to the following model [(Simpson et al., 2019)](https://doi.org/10.1101/533182):\n", + " $$A \\xrightarrow{b_0(1-\\frac{\\mathcal{C}(t)}{k})} 2A$$\n", + "\n", + "The model is simulated according to the Gillespie stochastic simulation algorithm [(Gillespie, 1976)](https://doi.org/10.1016/0021-9991(76)90041-3), [(Erban et al., 2007)](https://arxiv.org/abs/0704.1908):\n", + " 1. Sample a random value r from a uniform distribution: $r \\sim \\text{Uniform}(0,1)$\n", + " 2. Calculate the time ($\\tau$) until the next single reaction as follows:\n", + " $$ \\tau = \\frac{-\\ln{r}}{\\mathcal{C}(t)b_{0} (1-\\frac{\\mathcal{C}(t)}{k})} $$\n", + " 3. Update the population size at time t + $\\tau$ as: $ \\mathcal{C}(t + \\tau) = \\mathcal{C}(t) + 1 $\n", + " 4. Return to step (1) until population size reaches $k$\n", + " " + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "import pints\n", + "import pints.toy\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Specify initial population size, time points at which to record population values, initial birth rate $b_0$, and carrying capacity, $k$" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYgAAAEGCAYAAAB/+QKOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAX1UlEQVR4nO3de9RddX3n8fdHELRojECkkZAJrSw7EStqFKx0FoLOAKJxHAreWmSwzJpaxZnOVO2sNVJrV2Xa8V4ZUlFDvSBLRZh6mbIw3ipkTASlgq6mSiCZQFAuIiqKfuePvR84xhNyTp5zP+/XWs96zv6dfZ7zPWvD+eb7++7926kqJEna1UPGHYAkaTKZICRJXZkgJEldmSAkSV2ZICRJXe077gAG5eCDD65Vq1aNOwxJmiqbN2/+blUt6/bczCSIVatWsWnTpnGHIUlTJcnW3T3nFJMkqSsThCSpKxOEJKkrE4QkqSsThCSpq5EliCQ3JrkuybVJNrVjBya5Isk/tb8f3Y4nyTuSbEny9SRPGVWckqTGqCuIZ1XVUVW1pt1+HXBlVR0BXNluA5wEHNH+nA2cP+I4JWnujfs6iLXAce3j9cDngNe24xdVsxb51UmWJlleVTvGEqUkLdKHNt7EZdduH8rfXv3YJbzheU8Y+N8dZYIo4O+TFHBBVa0DDun40r8FOKR9fChwc8drt7Vjv5AgkpxNU2GwcuXKIYYuaVYM84v6wWz8zu0AHH34gSN/7701ygRxbFVtT/IY4Iok3+x8sqqqTR49a5PMOoA1a9Z45yOpT+P6shyncX1RH334gaw96lBecvT0/GN2ZAmiqra3v3cmuRR4OnDrwtRRkuXAznb37cBhHS9f0Y5JehD9fuFP479qF2sav6jHZSQJIskBwEOq6u728b8G3ghcDpwBvLn9fVn7ksuBP0xyMXA0cJf9B6m7zqTQ7xe+X5Z6MKOqIA4BLk2y8J4fqqrPJPkKcEmSs4CtwGnt/p8CTga2AD8EzhxRnNJU2F1S8AtfgzSSBFFV3wae1GX8e8AJXcYLeOUIQpOm0mXXbuf6Hd9n9fIlJgUNzbhPc5XUo86qYSE5fOQ/PGPMUWmWmSCkCba7qaTVy5ew9qhDxxma5oAJQppgTiVpnEwQ0oRzKknjYoKQJky3XoM0DiYIaQLYa9AkMkFIY+K1DJp0JghphEwKmiYmCGmEPCtJ08QEIQ3RrovneYGbpon3pJaGaKFiWGDTWdPECkIaMisGTSsThDRgXsegWeEUkzRgndNKTilpmllBSAPgSquaRSYIaS959bNmnQlC2kte06BZZ4KQ+uBUkuaJTWqpDzagNU+sIKQ+WTVoXpggpD3wugbNKxOE1IVnKEkmCKkrz1CSTBDSbtlr0LwzQUgtew3SL/I0V6nlKazSL7KC0Fzzwjdp96wgNNesGqTds4LQ3LNqkLqzgpAkdWUFobnj2UpSb6wgNHfsO0i9sYLQXLLvIO2ZCUJzwWklqX9OMWkuOK0k9W+kFUSSfYBNwPaqOiXJ4cDFwEHAZuB3q+onSfYHLgKeCnwPOL2qbhxlrJo9TitJ/Rl1BXEOcEPH9nnAW6vqccAdwFnt+FnAHe34W9v9pL58aONNnH7BVZx+wVX3Vw+SejeyBJFkBfBc4D3tdoDjgY+2u6wHXtA+Xttu0z5/Qru/1DOnlaTFGeUU09uAPwYe2W4fBNxZVfe129uAhf+DDwVuBqiq+5Lc1e7/3c4/mORs4GyAlStdq1+/zGklae+NpIJIcgqws6o2D/LvVtW6qlpTVWuWLVs2yD8tSXNvVBXEM4HnJzkZeBiwBHg7sDTJvm0VsQLY3u6/HTgM2JZkX+BRNM1q6UF5Oqs0OCOpIKrq9VW1oqpWAS8CPltVLwU2AKe2u50BXNY+vrzdpn3+s1VVo4hV082+gzQ4475Q7rXAxUneBFwDXNiOXwj8bZItwO00SUXqyns6SMMx8gRRVZ8DPtc+/jbw9C77/Bj4nZEGpqm1UDWsXr7EqkEaoHFXENJAWDVIg9d3DyLJAe0V0ZKkGbbHCiLJQ2h6AC8FngbcC+yf5LvAJ4ELqmrLUKOUduHZStLw9VJBbAB+HXg98KtVdVhVPQY4FrgaOC/Jy4YYo/RLPFtJGr5eehDPrqqf7jpYVbcDHwM+luShA49M2gP7DtJw7bGC2DU5dOtBdEsgkqTptscEkeQhSV6S5JNJdgLfBHYkuT7JXyZ53PDDlCSNmj0ISVJXPfcgkqyqqp8vDNqDkKTZtscE0dFf+DjwlM7nkhxTVVfbg9AoeGqrNFq99CBOS/Jm4JFJ/mV7XcSCdcMLTfpFntoqjVYvU0z/QLNE9yuAtwCPT3In8P+AHw0xNumXeGqrNDq9TDFtBy5K8s9V9Q8ASQ4CVtGc0SQNjdNK0vj0MsUUgIXk0D7+XlVtrqp7OveRBs1pJWl8epli2pDkY8BlVXXTwmCS/WhOdT2D5lTY9w8lQs09p5Wk8eglQZwI/Hvgw0kOB+6k6UnsA/w98LaqumZ4IWreOK0kTYZeehA/Bt4NvLu93uFg4EdVdeewg9N88gZA0mTo64ZB7fUOOxa2kyw1UWgYnFaSxq+nBJHkAOAJwJEdv48EDgCWDi06SdLY9HLDoBuBhwLX05zWegPwYuCoqto51OgkSWPTSwXxv4HjgL+pqksAkvxXk4MGyca0NHl6uR/Eq4BTgJOTfCXJSUANPTLNFa93kCZPTz2IqtoKvDzJE4A/A341ybOqasNQo9NcsTEtTZZe7gdxv6r6RlW9EHgW8N+SfH44YUmSxq2XJnWq6hemlKpqI/DsJCfsbh9pT+w7SJOtpzvKJXlVkpWdg+1SGyRZT7PchtQX+w7SZHOpDY2VfQdpcrnUhiSpq0UttSH1y76DND36OotJWiz7DtL06KuCkAbBvoM0HXpOEEn2B/4dza1G739dVb1x8GFJksatnwriMuAuYDNw73DCkSRNin4SxIqqOnFokUiSJko/TeovJ3ni0CKRJE2UfhLEscDmJN9K8vUk1yX5ei8vTPKwJP83ydeSfCPJn7bjhyfZmGRLko90XJ29f7u9pX1+Vb8fTJK0OP1MMZ20iPe5Fzi+qn7QXmz3pSSfBv4z8NaqujjJ/wLOAs5vf99RVY9L8iLgPOD0Rby/JKlPPSeIqtqa5EnAb7dDX6yqr/X42gJ+0G4+tP0p4HjgJe34euBcmgSxtn0M8FHgXS4IOL28OE6aTj1PMSU5B/gg8Jj25wNJXtXH6/dJci2wE7gC+Gfgzqq6r91lG7Bw1dShwM0A7fN3AQd1+ZtnJ9mUZNNtt93WaygaMS+Ok6ZTP1NMZwFHV9U9AEnOA64C3tnLi6vqZ8BRSZYClwK/0Wes3f7mOmAdwJo1a6wuJpgXx0nTp58mdYCfdWz/rB3rS7vI3wbgGcDSJAtJagWwvX28HTgMoH3+UcD3+n0vSdLe6ydBvA/YmOTcJOcCVwMX9vLCJMvayoEkDweeA9xAkyhObXc7g+ZiPIDLeeAeE6cCn7X/IEmj1U+T+i3tLUaf2Q6d2cd9IJYD65PsQ5OULqmqv0tyPXBxkjcB1/BAwrkQ+NskW4DbgRf1GqckaTD6Xe57M81SG32pqq8DT+4y/m3g6V3Gfwz8Tr/vI0kanF7uSf2lqjo2yd00p6be/xTNGayesyhJM6iXO8od2/5+5PDDkSRNin6W+z6vql67pzEJvDhOmgX9nMX0nC5ji1l+QzPMi+Ok6ddLD+I/An8A/Noui/M9EvjysALT9PPiOGm69TLF9CHg08BfAK/rGL+7qm4fSlSSpLHrpUl9F81aSC9O8mjgCOBhAEmoqi8MN0RJ0jj006R+BXAOzZIY1wLH0KzFdPxwQtM06WxKg41paRb006Q+B3gasLWqnkVz4dudQ4lKU6ezKQ02pqVZ0M+V1D+uqh8nIcn+VfXNJI8fWmSaOjalpdnST4LY1i649wngiiR3AFuHE5Ykadx6ShBJAry6Xar73CQbaJbg/swwg5MkjU9PCaKqKsmngCe2258falSSpLHrp0n91SRPG1okkqSJ0k8P4mjgpUm2AvfwwGquvzmUyCRJY9VPgvg3Q4tCU8kF+aTZ1s8U0x9U1dbOH5o1mjSnXJBPmm39VBDPAXZd2vukLmOaI177IM2uvV3NNcAjcDXXueO0kjQ/XM1VfVmYVlq9fInTStKM63k11yRnAi8EVi28rl3N9Y1DjVATx2klaT7004P4BM2y35uBe4cTjiRpUvSTIFZU1YlDi0QTy76DNJ/6Oc31y0meOLRINLE8nVWaT/1UEMcCL0/yHZopJq+kniP2HaT500+COGloUUiSJk7PCaK9clpzwr6DpJ57EGm8LMl/b7dXJnn68ELTONl3kNTPFNO7gZ8DxwNvBO4GPkZzn2rNIPsO0nzra7nvqnpKkmsAquqOJPsNKS5J0pj1kyB+mmQfoACSLKOpKDQj7DtI6tTPdRDvAC4FHpPkz4Ev0azPpBlh30FSp37OYvpgks3ACTTXQLygqm4YWmQaC/sOkhb0cxbTeuCWqvrrqnoXcEuS9w4vNEnSOPUzxfSbVXXnwkZV3QE8efAhSZImQT9N6ockeXSbGEhyYK+vT3IYcBFwCE2Te11Vvb39Gx+hWUL8RuC09uyoAG8HTgZ+CLy8qr7aR6zqkY1pSbvTTwXxP4Grk/xZkjfR3E3uL3t87X3AH1XVauAY4JVJVtPcgOjKqjoCuJIHbkh0EnBE+3M2cH4fcaoPNqYl7U4/TeqLkmyiuVCugBdW1fU9vnYHsKN9fHeSG4BDgbXAce1u64HP0dzjei1wUVUVTVJammR5+3c0YDamJXXTT5N6f+AoYAlwEHDqwrIb/UiyiqZ3sRE4pONL/xaaKShoksfNHS/b1o7t+rfOTrIpyabbbrut31AkSQ+inymmy2j+ZX8fcE/HT8+SPIJmeY7XVNX3O59rq4Xq5+9V1bqqWlNVa5YtW9bPSyVJezCyO8oleShNcvhgVX28Hb51YeooyXJgZzu+HTis873bMQ2AjWlJvRjJHeXas5IuBG6oqrd0PHU5cEb7+AyaKmVh/PfaFWSPAe6y/zA4NqYl9WJUd5R7JvC7wHVJrm3H/gR4M3BJkrOArcBp7XOfojnFdQvNaa5n9hGnemBjWtKejOSOclX1JZqE0s0JXfYv4JV7+36SpMXzjnJzwr6DpH71U0GQ5EnAb7ebX6yqrw0+JA3DQt9h9fIl9h0k9aTnBJHkHOD3gYUzkD6QZF1VvXMokWng7DtI6kc/FcRZNHeVuwcgyXnAVYAJQpJmUD+nuQb4Wcf2z9h941mSNOX6qSDeB2xMcmm7/QKaaxs0oWxMS1qMPSaIJI+jWTPpLUk+R3M9BMCr8ermiWZjWtJi9FJBvA14PUB7T4avArRXVb8NeN7QotOi2ZiWtLd66UEcUlXX7TrYjq0aeESSpInQS4JY+iDPPXxQgUiSJksvU0ybkvx+Vf1N52CSVwCbhxOW9paNaUmD0kuCeA1waZKX8kBCWAPsB/zbYQWmvWNjWtKg7DFBVNWtwG8leRZwZDv8yar67FAj016zMS1pEPpZrG8DsGGIsWgvOa0kaRj6uZJaE8obAEkahr5Wc9Xk6FY1OK0kaZCsIKaUVYOkYbOCmGJWDZKGyQpCktSVFcQU8WwlSaNkBTFF7DtIGiUriClj30HSqFhBSJK6MkFIkrpyimnC2ZiWNC5WEBPOxrSkcbGCmAI2piWNgwliAjmtJGkSOMU0gZxWkjQJrCAmlNNKksbNBDEhnFaSNGmcYpoQTitJmjRWEBPEaSVJk8QEMUZOK0maZE4xjZHTSpIm2UgqiCTvBU4BdlbVke3YgcBHgFXAjcBpVXVHkgBvB04Gfgi8vKq+Ooo4x8FpJUmTalQVxPuBE3cZex1wZVUdAVzZbgOcBBzR/pwNnD+iGEfiQxtv4vQLruL0C666v3qQpEk0kgRRVV8Abt9leC2wvn28HnhBx/hF1bgaWJpk+SjiHAWnlSRNi3E2qQ+pqh3t41uAQ9rHhwI3d+y3rR3bwS6SnE1TZbBy5crhRTpgTitJmgYTcRZTVVWS2ovXrQPWAaxZs6bv14+KZytJmkbjPIvp1oWpo/b3znZ8O3BYx34r2rGp5bSSpGk0zgricuAM4M3t78s6xv8wycXA0cBdHVNRU8tpJUnTZlSnuX4YOA44OMk24A00ieGSJGcBW4HT2t0/RXOK6xaa01zPHEWMg+a0kqRpN5IEUVUv3s1TJ3TZt4BXDjei4VuYVlq9fInTSpKm0kQ0qWeV00qSpplLbUiSurKCGCD7DpJmiRXEAHk6q6RZYgWxSN2qBvsOkmaBFcQiWTVImlVWEH3qrBjAqkHS7LKC6FNnxQBWDZJmlxXEXrBikDQPTBA98PRVSfPIKaYe2IiWNI+sIDrs2oBeYCNa0jyyguiwawN6gVWDpHlkBbELKwVJasx9grABLUndzf0U02XXbmfjd24HnEqSpE5zX0GsfuwSVj92CW943hPGHYokTZS5TxAmBknqbu6nmCRJ3ZkgJEldmSAkSV2ZICRJXZkgJEldmSAkSV2ZICRJXZkgJEldparGHcNAJLkN2LqXLz8Y+O4Aw5kGfub54GeeD4v5zP+iqpZ1e2JmEsRiJNlUVWvGHcco+Znng595PgzrMzvFJEnqygQhSerKBNFYN+4AxsDPPB/8zPNhKJ/ZHoQkqSsrCElSVyYISVJXc58gkpyY5FtJtiR53bjjGYYkhyXZkOT6JN9Ick47fmCSK5L8U/v70eOOdZCS7JPkmiR/124fnmRje6w/kmS/ccc4SEmWJvlokm8muSHJM+bgGP+n9r/pf0zy4SQPm7XjnOS9SXYm+ceOsa7HNY13tJ/960mespj3nusEkWQf4K+Bk4DVwIuTrB5vVENxH/BHVbUaOAZ4Zfs5XwdcWVVHAFe227PkHOCGju3zgLdW1eOAO4CzxhLV8Lwd+ExV/QbwJJrPPrPHOMmhwKuBNVV1JLAP8CJm7zi/Hzhxl7HdHdeTgCPan7OB8xfzxnOdIICnA1uq6ttV9RPgYmDtmGMauKraUVVfbR/fTfPFcSjNZ13f7rYeeMF4Ihy8JCuA5wLvabcDHA98tN1l1j7vo4B/BVwIUFU/qao7meFj3NoXeHiSfYFfAXYwY8e5qr4A3L7L8O6O61rgompcDSxNsnxv33veE8ShwM0d29vasZmVZBXwZGAjcEhV7WifugU4ZExhDcPbgD8Gft5uHwTcWVX3tduzdqwPB24D3tdOq70nyQHM8DGuqu3AXwE30SSGu4DNzPZxXrC74zrQ77R5TxBzJckjgI8Br6mq73c+V835zjNxznOSU4CdVbV53LGM0L7AU4Dzq+rJwD3sMp00S8cYoJ13X0uTHB8LHMAvT8XMvGEe13lPENuBwzq2V7RjMyfJQ2mSwwer6uPt8K0L5Wf7e+e44huwZwLPT3IjzbTh8TTz80vbqQiYvWO9DdhWVRvb7Y/SJIxZPcYAzwa+U1W3VdVPgY/THPtZPs4LdndcB/qdNu8J4ivAEe1ZD/vRNLguH3NMA9fOv18I3FBVb+l46nLgjPbxGcBlo45tGKrq9VW1oqpW0RzTz1bVS4ENwKntbjPzeQGq6hbg5iSPb4dOAK5nRo9x6ybgmCS/0v43vvCZZ/Y4d9jdcb0c+L32bKZjgLs6pqL6NvdXUic5mWa+eh/gvVX152MOaeCSHAt8EbiOB+bk/4SmD3EJsJJmqfTTqmrXZthUS3Ic8F+q6pQkv0ZTURwIXAO8rKruHWd8g5TkKJqm/H7At4Ezaf4ROLPHOMmfAqfTnKl3DfAKmjn3mTnOST4MHEezpPetwBuAT9DluLaJ8l00U20/BM6sqk17/d7zniAkSd3N+xSTJGk3TBCSpK5MEJKkrkwQkqSuTBCSpK5MEFKfkhyU5Nr255Yk29vHP0jy7nHHJw2Kp7lKi5DkXOAHVfVX445FGjQrCGlAkhzXce+Jc5OsT/LFJFuTvDDJ/0hyXZLPtEufkOSpST6fZHOS/7OYlTelQTNBSMPz6zTrQD0f+ACwoaqeCPwIeG6bJN4JnFpVTwXeC8zclfyaXvvueRdJe+nTVfXTJNfRLOXymXb8OmAV8HjgSOCKZoUE9qFZtlqaCCYIaXjuBaiqnyf5aT3Q8Ps5zf97Ab5RVc8YV4DSg3GKSRqfbwHLkjwDmiXZkzxhzDFJ9zNBSGPS3ub2VOC8JF8DrgV+a7xRSQ/wNFdJUldWEJKkrkwQkqSuTBCSpK5MEJKkrkwQkqSuTBCSpK5MEJKkrv4/rlfX08Vi+jgAAAAASUVORK5CYII=\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "n_0 = 50\n", + "model = pints.toy.StochasticLogisticModel(n_0)\n", + "\n", + "times = np.linspace(0, 100, 101)\n", + "\n", + "# $b_0$ = 0.1, $k$ = 500\n", + "params = [0.1, 500]\n", + "\n", + "values = model.simulate(params, times)\n", + "\n", + "plt.step(times, values)\n", + "plt.xlabel('Time')\n", + "plt.ylabel(r'Concentration ($A(t)$)')\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Given the stochastic nature of this model, every iteration returns a different result. However, averaging the population values at each time step, produces a reproducible result which tends towards a deterministic function as the the number of iterations tends to infinity: $$ \\mathcal{C}(t) = \\frac{k\\mathcal{C}(0)}{\\mathcal{C}(0) + e^{-b_{0}t}(k-\\mathcal{C}(0))} $$\n" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAEWCAYAAACT7WsrAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3dd3gc9bXw8e/ZJmnVm7vcwJYLxdgGY2M6BIMNhtBJaKGFhLyESwIkNyEJN7kJ9yYhPYEEDAZCzXUgYFro1cbGxt244KIuW7K6tJL2vH/MyF4LrSXZklblfJ5nn90pu3NmZ2fO/srMiKpijDHGtMUT6wCMMcb0XpYkjDHGRGVJwhhjTFSWJIwxxkRlScIYY0xUliSMMcZEZUmiB4jINhE5oweW8xcR+eFBvK9L4juE5Y8UkWoR8R5qDLEgIqNFREXEF+tYupuIvCUi17uvvyIir0ZMO0FENrnb8nwRGSwi74hIlYj8KnZRt01EXhKRq2McQ7WIjI1lDO0Z0ElCRGaLyAciUiEiZSLyvogc6067RkTei3WM0bQVn6p+XVX/K1YxdXT5rZOSqu5Q1SRVbe7eCLtGTyX93k5VH1fVL0WMugf4g7st/wncCOwCUlT19p6MrSP7r6qeraqPdHT+Lohpb4KNiCFJVbd253IPVb//5xONiKQALwA3A08DAeBEoCGWcZnuISICiKqGYx1LVxERby9LrKOAta2G1+lBnLErIj5VbeqyyLpZX4u3U1R1QD6A6cCeKNMmAvVAM1DdMh+QCiwESoHtwA8AT8T7bgDWA1XAOmCqO34b8B1gFVABPAXEu9PScZJVKVDuvh4R8ZnXAFvdz/wc+MoB4nsY+GnEe+cDK4FKYAswJ8r6bgPOcF/HAb8BCtzHb4C4iHnvAArdadcDChzeevlAlrsue4Ay4F2ckuujQBioc2O/Axjtfo7PfW8GsMBdRjnwzyhxe4Ff4fxb/Ry4pdXnvAX8DHjfXd7hwCzgY3c7fAzMcuc9FVgd8dmvAR9HDL8LnN9O/FcDO9x4/vMAv725wAp3u+wEftxq+mzgA/e72wlcE/H9/hlYDNQAZ7i/hbfcedcC50V8zjk4v8MqIB/4zoG2TZRYzwQ2uN/XH4C3gesjfpvvua+3tPpengAagZA7fIa7/e9y592N8+csw31/y3d4nfsdvuOO/xrOPlUOvAKMiohNga8Dm9x1+SMgRNk/2li3t3B+w9H2pzjgl248xcBfgAR32ilAHnAnUOT+LqLuyzi/w2Z3OdU4Ja6WdWjZf6IeX1q+azeecpzf+9kHOk502bEy1gfrWD2AFPeH+ghwNpDeavreHSBi3ELgOSDZ/VF/BlznTrsYZ0c81v2hHt7yg8Y5CC8FhuEcANcDX3enZQIXAkH3c5/BPSgCiTgHklx3eCgw+QDxPcy+g/RxODv2mTg753BgQpTvYhv7ksQ9wEfAICAb52D1X+60Oe4OMdmN9zGiJ4mf4+xUfvdxIs4/+f2W1+oA0XJwfxEnkaa77z05StxfxzkIjnDn/TdfTBI73Hh9wGCcHexKd/hydzgTSMDZgbPcZRa72zPZnVYHZLYT/1/deY/GKZFOjBL3KcCR7nY5yl3W+e60UTg7+uVuHJnAlIjvtwI4wX1vMrAZ+D5OSfg0970tv5dC4ET3dTr7/rRE3Tat4sxyP+8id77bgCbaSBJRvpe9vwd3+Fac39YInAPw/cATrb7DhTi/+wScPzmbcQ7iPpyD5gcRn6c4B+I0YCTOwXVOtP2jjfV7K9q6uOPuA57H2WeTgX8BP4/Yhk3Ave66JHCAfbn18lqtQ8v+c6DjyzU4SfcGnD9HN+P8iRIOcJzokmNlTxyQe+vD/fE9jPOPoMn9QQyOsgN4cf4VTYoYdxPwlvv6FeDWKMvZBnw1Yvh/gL9EmXcKUO6+TsT5h3Qh7j+YiPna+lE/zL6D9P3AfR38HraxL0lsAc6JmHYWsM19/VDLTuIOH070JHGP+4M//EDLc4dHu5/jc3/gYVol7ShxvwHcFDF8Bl9MEvdETL8SWNrqMz5k3z/1d4EvA8cDr+L8052DU8pY1YH4I0uAS4HLOvj9/6ZlWwHfAxZFme9hYGHE8Ik4STuyNPsEbskEJ0HehNMmEPk5UbdNq/muAj6KGBacfeVgk8R64PSI4aE4Bz5fxHc4NmL6S7gHSXfYA9Sy78+XArMjpj8N3BVt/2hj/d46wLoITmntsIhxM4HP3den4BwP4g/w+Xv35dbLixinOPtRe8eXa4DNEdOC7nuHcIDjRFc8BnTDtaquV9VrVHUEcATOP/3fRJm95R/m9ohx23H+oQPk4BxgoymKeF0LJAGISFBE7heR7SJSCbwDpLn1zTXApTj/mAtF5EURmdDB1WsvnmiG8cV1HBYxbWfEtMjXrf0vzr/AV0Vkq4jc1cHl5wBlqlrewVjbiydyXOt1g/234ds4O/9J7uu3gJPdx9sdiKfNbdyaiMwQkTdFpFREKnC2b5Y7ub3t1np9dur+7SyR63MhTpXTdhF5W0RmuuM7um32+37VOTodaJu3ZxSwSET2iMgenKTRjFPCa7Gz1fy/jZi/DOfgPTxing595wchG+dAvDxi+S+741uUqmp9y8CB9uUOLK+94wtErKuq1rovkw7xONGuAZ0kIqnqBpx/Pke0jGo1yy6cfz2jIsaNxKmSAOfHfdhBLPp2IBeYoaopOAcocHYGVPUVVT0T51/XBpwqjbbia+1g4yngi+tY4L4uxKkqaJET7UNUtUpVb1fVscB5wH+IyOktkw+w/J1AhoikdSDWjsQTuazW6wb7b8PWSeJt2k4S7X337fk7Tqk1R1VTcap+xJ3W3nZrvT45IhK5H+9dH1X9WFXn41Qd/hPnn3Z72yZSIRHfqdv4H3Wbd8BOnHr0tIhHvKrmR8yjrea/qdX8Car6QQeW1dlt1Nb+XodTbdOy7FRVTTrAew64L7cTU3vHlwMHH/04ccgGbJIQkQkicruIjHCHc3DqgT9yZykGRohIAECdXiRPAz8TkWQRGQX8B069PMDfgO+IyDRxHO7O055knB/jHhHJAH4UEeNgEZkvIok4ddzVOFUxX4ivDQ8C14rI6SLiEZHhHfx38QTwAxHJFpEs4O6IdXza/cyJIhIEop4TISLz3O9AcOrRm1vF3mbfcFUtxKlm+JOIpIuIX0ROamteN55b3XVLw2lEPJDFwHgRuUJEfCJyKTAJp14bnPaXXJz2nKWquhZnp52B86+wRdT4OygZp7RULyLHAVdETHscOENELnFjzBSRKVE+ZwnOv+c73O/pFOBc4EkRCYhzHkOqqjbi1FmHod1tE+lFYLKIfFmcc0D+H071xsH6C87+M8qNI1tE5rcz//dEZLI7f6qIXNzBZbW3fxxwfrd09lfgPhEZ5C5/uIicdYDPiLovRywj2u++veNLVO0cJw7ZgE0SOA1yM4AlIlKDkxzW4PwbAKe+ey1QJCK73HHfwqmn3IrT0+DvOPX0qOozOD0Y/u5+9j9xGrza8xucRq9dbgwvR0zz4PxQCnCK2ifjNFhFi28vVV0KXIvT+FaB80+4I0nrp8AynJ5Yq4FP3HGo6kvA74A3caorWhJqW92Gx+E0JFfj1Pv/SVXfdKf9HCcR7RGR77Tx3itx/lVtAEqAb0eJ9a84bQercHoLLcZpW2qzW6iq7gbm4Wzj3Tg9k+ap6i53eo27vmtVNeS+7UNgu6qWRHxUe/G35xvAPSJShZOEn46IcQdOFdHtONt8JU5DeFvrE8JJCmfj/H7+BFzllorB+R63uVUfX8fpGQcH3jaRn78Lp0PGL3C+r3E4PcUO1m9xSlCvuuv+Ec4+2CZVXYTTMPykuw5r3HXtiAPuHx2c/07c37m7/H/j/ImI5kD7Mjjrf5GIlIvI79p4f9TjSzsOdJw4ZC29TYzpNBGZiLPjxmkv6CMuImfjdAjoSDI0xnTAQC5JmIMgIheISJyIpOP8y/tXrBKEiCSIyDlutcxwnOL9oljEYkx/ZUnCdNZNOFVAW3CqdbqsWHsQBPgJzrkOK3B6y9wdw3iM6XesuskYY0xUVpIwxhgTVb+6wF9WVpaOHj061mEYY0yfsnz58l2qmt3WtH6VJEaPHs2yZctiHYYxxvQpItL6SgR7WXWTMcaYqCxJGGOMicqShDHGmKgsSRhjjInKkoQxxpioejRJiHMD+dUislJElrnjMkTkNRHZ5D6nu+NFRH4nIptFZJWITO3JWI0xxsSmJHGqqk5R1enu8F3A66o6DnjdHQbnao/j3MeNOPf2NcYY04N6w3kS83Fu9ALO/abfwrlE73ycWzUqzqV600RkqHu/AWPMAPHMZ8+weOvirvkwVVCQMIg6wxLGGadKxo4xpBcMd+bDfajS1BymKdzsjmfveFD3jkJKUl0Twfpmd6jl/XToed/FkTTiqfV8uMvbf7iF3wc3PPKLTn0dHdHTSUJxriWvwP2q+gDOPaVbDvxF7LuV4XD2v5VhnjtuvyQhIjfilDQYOXJkN4ZujOmovy/ZwXMr276pmoSb8DaH8DQ1Os/NITzNjdSxgnr9DG8YfM3gbVa8YdBwGd5mIS6cxvDSENl7GlFxD5fiHDKdZ/fwud84YO98su8ecVGVAcu/MNbfwfWubX+WbiMyuls+t6eTxGxVzXfv9PSaiGyInKiq6iaQDnMTzQMA06dPt6sVGtNJ5U89TeULL3xhfHFVA7urv3g/qQb20OipodnjISyefc8ihMVDWIRmEU4SD7jD4YhpSLtH6v1pKogXwYsQoC7Og3PbaA/iPiJfOzfcaxkve6c3BDIBDwmN5e54Afche2vehdT6HaTXbEZaSgkKguLzCH6PZ+94ccfDvtyTODqB1HGJe1dRWvKS+ywRz3vfF5G7IsdHfkbkMmTiPJhymTtt3/cZiE/o3PfaQT2aJFruZauqJSKyCOc2kcUt1UgiMhTnMtTg3Ns18n66I+jg/V6NMY5oCSBS7ccfo8CWnPE0eqDJozR5oIY6wsEmwl6vc+D3eGgWL2FPEpDU5meJhvFoGK+G8aji1Wb84X3DnnCzMx2ojxtLQ8JhCD7Ai4gPP178eJG9ScFLeXA8ABl1mwHwBb34kg7u0DX+6CCTv9rRO6Aa6MEk4d5/1aOqVe7rLwH34NzO8GqcWyReDTznvuV54BYReRLnFocV1h5hTPsiE0Ptxx8DEDz2WArrSilqqCAkXho9HkIeLyHx0nDESEJeH/KFu77GURfwIYEwTf4QTf5mmhK8NAfCTPD7OTEQR3wAEvwQ74c4P/i9Anjdh2NtyZF8tmvivo9NzIbkIZRv2gPAsHFpB1yfBGD8cYOZfOJph/jNmIPRkyWJwcAipyiID/i7qr4sIh8DT4vIdcB24BJ3/sU49/rdjFPVd20PxmpMn/Lv391J8yvObaqHbKuhKj7AjhEpNEwcSllaAo3swe/1QnDfbdfV40e9QfAFUW8C+N1nbwLqjQdvHIPET3bLAd89uLdY0/Ki0X1EUfB528lg2Lg09+A//FBX33SjfnXToenTp6tdBdb0N/tVGVUVQU0pAGuSj6UgaSJN4QqkcQ+NVNBIJcq+dgQVH+pPQn1J4Esk0ZtA0JNOTdNExBNPatznADQnZJGSnfOFZXcVSwa9m4gsjzgtYT+9oQusMaaVtqqMth2WTGn8BKoGHUVTuJwm3QMNr7vvEDyeZOICg/D4UvB5UwkEUvB44t2GXEdLMkjDqnBMx1iSMCaGIpNBaV0pu+t2A+CtPYrCQcdS4y8jNDWVEGXQXIlQBk1lNPuSaQ5kU5aexu4hx/D2Lj/NHh8zxjjVSfOnDOeKGdYl3Bw6SxLG9LDy/7mVytfeAWBT0zSKB8+iPk5oCoYJJ1QSDpfQkFqJNr8KjYqKB18gg+TEXPCmU5oYz4rxs/b7zOmplhhM97AkYUw3OFDX0015AYoH30R9nFAXfxjhpjzqWY6nrhBprgcgLpDO8LFHMO2yy3m3IpHn1zg9w9cVVjJpaApP3TSzx9bFDGyWJIzpBlueWYB3805KhgehOQTN+7r/bMu5hTq/EPKsx1P5KhJuIC4Qx+ipx3DYtBmMOWY6iWnpe+d//v4P9yaHSUNTmD/FGoBNz7EkYUwXieyGmplXxZbB8OLFYbLyppJZfiw01eCp2Ym3+k1oqCU+Lp7DZh7P+ONPYPTRU/HHxe/9rMjLWljpwcSSJQljDkHkxefmvvARo4th9yAPW8fOonrQNOatjiOvrIrmhqVocymIl2DaeMbPOJGTrvhS1MSw5PMyAGaMybDSg4kpSxLGdEC0K5HmvOBjZtUMvOKjLnMG64ZA8uAESvak0lS3goqqtYSb60nOyuHYc28id9ZJBFNS21zGcyvz95YaZozJsIZo0ytYkjCmAxZvXczGso3kZuTuN35U9bE0+oeR3FgAPgjH1VBcvZ2GyrWIB8bPOIFj5pzLsNyJe89XiHaFVKtWMr2RJQljoogsPTSvTea8slvIzZiw31nPJf5MUut3ctIxL/DRFti2SwgkJJA47VTej5vEkvgUHn27At7+aO/nRlYlRbJqJdMbWZIwJorI0sOk0qmkVKVBaDXUVzgzxKcSV7+aqtAKnlwiJKSkMvvy85nypXO46tFVbsngi59rVUmmL7EkYcwB5GbksmDOAha9/gj4tnHBpBcpX1nJjh1B1sQF2a2NxCOcctUNHHX6Wfjj9zVEW9WR6Q8sSRgT4ZlXb2NxwbsANBdPY9Ku6Sx6/RF2VWeSlQTV8xfy+qs3sD0Qwq9NHOlN4Ji5F5A9d36b3VaN6essSZgBL7LtYVnxMhCYrnFM2jWdlNrhECwlVQsI7fyQv938DKpKrjeBM/72MPGJ+26+E9k7ydoXTH9hScIMeJFtD9M1jnM0kYuvfZtFv/oEgKlfGsbLP72basIMEz9H+YMMO3c+8YlJdtKb6fcsSRjDvraHtffexWe7JrLoV59QvCWPcOWrPLt0G8HGJmYnpDHj8Sf3e5+VHkx/Z0nCDHhZGwaRmX84i15/hIKqLwGQElhLXdk/0XADE73xTAgmkH7uuW2+30oPpj+zJGEGpMh2iJz8mSTUDoXgbjJ1I3V73qWkvJLU+kZmJGdx5GOP7/dea6A2A4klCTMgLf70ITbW5JNLgKDOICFYwrHXHsfzP/k+9drMZF+Q3OR00tzSg11XyQxUliTMwLFsAax+FoCs/KHklp9HridIaW0mXlnG0/f8k0TgNH8qUx7/+35vtesqmYHKkoQZMJ5Z/TCLG0sgkEhOxXwSGkbQPCKJUMnfaajPY5QnwMSt+STlpgF2uW5jADyxDsCYnrJYatgY8MOQIwkGkgiklLJ7x0M01udxdEklx/qSSMrNJWXePGBf6QHsukpm4LKShOnfIqqYaKwhN5DIgjkLWPjW3yjd9i9SsrI41ZdC+ogMRj26cF/pIeJucFZ6MAOZlSRM/7b6WShaDUBW+Wnkrv8mf/v2Lyj9/J8EPNmcvKeBhM82753dSg/G7M9KEqbf2e8GQVIMQwfBkEHkrDuKQNEKKkJbCfoPZ9zuRuKyQhQPH8NzKRNZbqUHY77AkoTpdyK7t2YVTCGzfAZs8OEveA9tKuaUq26g6a+PU9YY4o7Zd+7r0oqVHoxpzZKE6X9qSskNNbLAn8Oi8hMpqU2hofp5mprLmXLOzUybO5eX7ltATUMTYPd3MOZALEmY/imQCNe8SOinr1BXtBCvt5E5s08j8d9Psf3fTzGkdAdF2SOtWsmYdliSMP1CZDvERkLkEmDXjm0UrH8IAS655+fU3/0T6jdsIH7CBIqyR7I6dwZnxzZsY3o9SxKmX4hshzih4DjGFk/i0dfuQMMehk3+GoPHHMZ2IH7CBEY9upA77v8w1iEb0ydYkjD9QtbWseSWzidXEsgrj6O2+lm8/jiGT7qGySdNAqC4qoHd1Q3cEdGLyRhzYJYkTJ+135Vcd80koW44oSEemvIfIj4pka/+972kDhqyd/7d1Q17G6utF5MxHWNJwvRZkXeUCyLE+TdTtvNjgqlBLvvJ/gmiRWKczxqrjekESxKm76oqIjcUYkFhCU/XVlBQvoRAYoCLfvDTvQmi/KmnqXzhBYC9PZqMMR1nScL0WVlbx5JZOp9nw0J+2TtAmIv+87/wvPch21/4TwBqP/4YgG0jcqlJGcZW69FkTKf0eJIQES+wDMhX1XkiMgZ4EsgElgNXqmpIROKAhcA0YDdwqapu6+l4TS8TccG+zNJzia8dQmHdYpRajrvgNgaNHsv2H/54b1fXbSNyeWXI0RSeOAfA2iGM6aRYlCRuBdYDLV1L7gXuU9UnReQvwHXAn93nclU9XEQuc+e7NAbxml5k4eurKCg9F8RLoHYwjbXPEW4o5fQZp5D2r7+x/V9/o3LtOvIyRvDw7JvtWkzGHKIevQqsiIwA5gJ/c4cFOA1wr+XMI8D57uv57jDu9NPd+c1As2wBLJgLC+ZSUJpLXF0O+IM0h17FU1fMmTfcQtryldRv2ABAXsYIXhlyNGC9mIw5VD1dkvgNcAeQ7A5nAntUtckdzgNa9ujhwE4AVW0SkQp3/l2RHygiNwI3AowcaY2S/VLL5b6HHAnipSG5lJOnxfHB09uY5E0gZcFjVnowppv0WJIQkXlAiaouF5FTuupzVfUB4AGA6dOna1d9ruk91pYcyWdll0DgSBLqSqj3rOCDp99npCfAmE3bYcIEKz0Y0016siRxAnCeiJwDxOO0SfwWSBMRn1uaGAHku/PnAzlAnoj4gFScBmwzwLy/azx1tVnUlW2gwVdAYvGHDBs/kWmfF+FtdZkNKz0Y07V6rE1CVb+nqiNUdTRwGfCGqn4FeBO4yJ3tauA59/Xz7jDu9DdU1UoKA1AZzewO5rNl2ssklbxJQnMzUz8vpHHjxliHZky/1xvOk7gTeFJEfgqsAB50xz8IPCoim4EynMRiBoiF999HweYkAOJqh0JCPqetGExxcxnH5e8mfmwmTJhAyrx5MY7UmP4tJklCVd8C3nJfbwWOa2OeeuDiHg3M9BoFm5OIqx1KQ7CQhmAhyc2fUPBZPjO8iQwam8GoRxfy9yU7eG5lPtgF+4zpNj3aBdaYzmgIFnLX/97AvHNGU5+fz/Rzv0yON27v9OdW5rOusBKwxmpjuktvqG4yBoC17+bz2dJiABJqh1EXLKD48y38+4HfM0h8jHz5Teo3biR+woS977GursZ0L0sSptd4/5X3qStPpC5YQG0wTH3aSv513xICYeXobYV4xqdQMXwMz6VMZLlVMRnTIyxJmF6jrKGC2mA5Oyf/FRSOXX8YVbvqOcmXROr49L1dXdcVVjIJq2IypidYkjC9ShAPC65ZxrIXFvH2yw9y8pXXkfXUov3msSomY3qOJQkTUwuffZ6CFdUAxNUOoyFYQOGmjbzz2EMMEz9ZT/7fF9ohjDE9x3o3mZgqWFFN3J5UABqCBQxJX80Lv/0fElQ48vMCRIR4Ox/CmJixkoSJuYa0Cu762RXoQ3N54VPYtNvLKb5Eksfn2vkQxsSYlSRMr7Hkwzo+KxImSxyJn23ZO97OhzAmdqwkYXrcfudDVKRRl5BP2e/n8lFFPFm1deSmpSOtqpissdqY2Oh0khCRRKBeVZu7IR4zAHy2tJhdedVkjUiiLiGf3Wkf8dJq8CrM8HkZ/dijVsVkTC/RbpIQEQ/OxfW+AhwLNABxIrILeBG4X1U3d2uUpl8prSulLLibDya/xsbM5czenE5RRTIz/EkkxDuX3WipYpo0NMWqmIyJoY6UJN4E/g18D1ijqmEAEckATgXuFZFFqvpY94Vp+rrIKqa6IqU2sQ6AoyuCDN6cRO6skwi8uZx11ZX7TpizKiZjYq4jDddnqOp/qeqqlgQBoKplqvoPVb0QeKr7QjT9QUsVE0BdUhF1GZ/w1/wSjv8kkcQAnH7dzeyubqCmwbmTrZUejOkd2i1JqGpj5HBbbRKt5zGmLVkjkrjg9qlc+/CNEKphyZax7KpP4ITGRkpu/iZDSndQlD3SSg/G9CLtliRExCMiV4jIiyJSAmwACkVknYj8r4gc3v1hmv4mvj6Vpdt8DGr2EbellHWFlWxJGcbq3BmxDs0YE8HaJEyPKK0rpaxuN9e+/Hs+C4c4e3Um8UnJpH9exZaUYTxz8Z0AVsVkTC/TkSRxhqo2isjo1m0SwD+Af4iIv9siNP1CWd1uapucxuoTPk8nodLPGbd/g83f/x8S43xWxWRML9VudVNEe8P/tZ4mIse3mseYqIK+BH555D0M3ZTEkNQQP1ghexuqjTG9U0fOk7gEmAoki8hEYGNEieIB4KhujM/0F80haArx2n99HZ82M2FjERlb72VIZQGNY8bFOjpjTBQdqW56H4gHrgd+DeSKyB6gAKjrxthMH7ffZcCrsmlsfpf8cuHIij0EymHS0Skw9Ai7wqsxvVhHusDmAwtFZIuqvg8gIpnAaJyeTsa0qeUy4A1pFTQEthAoWkfOEUfT9PYa8gYJZz+6MNYhGmPa0ZEusALQkiDc17tVdbmq1kTOY0xrDWkV3PnTyxlZ8i+kuYnAu2sYVLwj1mEZYzqoI2dcvyki3xKRkZEjRSQgIqeJyCPA1d0TnukPNi/9kMKmOMaVlBMIC6WDR+E5Y06swzLGdEBH2iTmAF8DnhCRMcAenDYKL/Aq8BtVXdF9IZo+LdzIGw/fT3JjiNFazRH//jDWERljOqEjbRL1wJ+AP7nnQ2QBdaq6p7uDM31LZEM1QOLuVKTyFaqbdjErrxhPQgyDM8YclI60SfxQRL4JzvkQqlpoCcK0JfJ+1QCBmrU0NG0lp66a5LhmNo+33tLG9DUdqW66BOc+EvsRkeuBbFX9eZdHZfqsvferDod58MKzCHh9/POsH/EPf4JdcsOYPqgjSaLRrXJq7VHgE8CSxAD2zGfPsHjrYgByGmcRVGDBXFbvhIpAHEft2cW3bjkttkEaYw5aR5JESESGqmph5EhVbRARuxzHALf0jU3kfD6boC+BpMpUUmp38tmnRbydOpj0unoSQ3ZZL2P6so50gf0V8JyIjIocKSKDAO2WqEyfkblzNJk1w8nNmEBmXT7DC5ezJnM0jR4P6Y1JvKvRVfUAABwmSURBVDfpnFiHaIw5BB3p3fSMiASB5SLyEbASJ7lcDPy4e8MzfUFdahkX3D6HTaddS6nHx4rwUIpyZrBg4nFMGpoS6/CMMYegQ2dcq+ojwBic25T6gXrgClV9vGWebo3S9Amh5jBrkzNoDCSRP+ZEuwWpMf1Ah246JCL/AJ5T1UdbRraccY1ztvWbwMPdE6LpzRIrGwlWNbL9yqvYGQ5SGYjj/Ju+zvdmnxLr0IwxXcDOuDaHJFjVQKAhTGj3GjYNGUS6J8SEE06OdVjGmC7SY2dci0g88A4Q5y73WVX9kZt4ngQygeXAlaoaEpE4YCEwDdgNXKqq2zqzTNM9Is+sTvEPA81jx9FDCe0A3+TJWO2jMf1HR3o37XWIZ1w3AKep6tHAFGCOe2e7e4H7VPVwoBy4zp3/OqDcHX+fO5/pBSLPrI4L5ZNcuYQVOzysTz+Kt7IvinF0xpiu1KkkcSjU0XJhH7/7UOA04Fl3/CPA+e7r+e4w7vTTrYG892g5s/r4bX+gWnfS7Iun7ogzrKHamH6mI20SXUZEvDhVSocDfwS2AHtUteVGx3lAy1FmOLATQFWbRKQCp0pqV6vPvBG4EWDkyP2uZm66SeKeOoLVTWw/8xiKm+Ioj4tnzjXXcueZp8c6NGNMF+uxkgSAqjar6hRgBHAcMKELPvMBVZ2uqtOzs7MPOUbTvmB1E4GQ0gxsGJpFwOPhyNPPinVYxphu0OGShNuQfCHObUv3vk9V7+nsQlV1j4i8CcwE0kTE55YmRgD57mz5QA6QJyI+IBWnAdv0AqGAUPqN71H35ELWH305Ho831iEZY7pBZ6qbngMqcKqLGjq7IBHJxrlY4B4RSQDOxGmMfhO4CKeH09XucgCed4c/dKe/oap2GZAYibyQ3yyOB61nyaKnKcsaT1X66NgGZ4zpNp1JEiNU9VDuOTkUeMRtl/AAT6vqCyKyDnhSRH4KrAAedOd/EHhURDYDZcBlh7Bsc4iWLl5GTv5MgghNgaE0V79Ig4Z4M/V4cmIdnDGm23QmSXwgIkeq6uqDWZCqrgKOaWP8Vpz2idbj63GuD2V6gZEbD8PXNIy0UAENDWvZrTspHnk8OaNyrEeTMf1YZ5LEbOAaEfkcp7pJcHq22u3GBoBgrRII5THT+z4vhcrRgJ+f/eQ24oKJsQ7NGNONOpMkzu62KEyfEAoIPx93OuNXP8u28adbgjBmAOhwklDV7SJyNHCiO+pdVf20e8IyvcEzr97G4oJ3AZjF9UhYGbn5dWoTs5h59twYR2eM6QkdPk9CRG4FHgcGuY/HRORb3RWYib3FBe+y0b1zrQch1Lyd+LpyvvKtb/KVmWNiHJ0xpid05mS664AZqnq3qt4NHA/c0D1hmd4iV+JZcM0ywuqhvvkzRh11DGOmTI91WMaYHtKZNgkBmiOGm91xpp+a/EmISeuaeOmx82nMHIkS4uSvfs2u8mrMANKZJLEAWCIii9zh89l3ToPpJyJPmjs97zh2Z0+nzhemLvQygYQJZI+yaiZjBpLONFz/WkTeBk5wR11rNxvqfxZ/+hAba/LJJUBD8vWEAiMQ33vIHi+zLr0m1uEZY3pYp64Cq6rLcS7LYfqrmlJyQ40s8OfweFjQunVUNK7h+AsvZ9rZR8Q6OmNMD2s3SYjIe6o6W0SqcO7/sHcSzsl0Kd0WnekRkVVMzcXTmLRrOouGTKMirohQ9QcEU9M49rwvxzhKY0wstNu7SVVnu8/JqpoS8Ui2BNE/LN66mI1lGwGYtGs6KbXOZTakfiXNupsTLvkqgfiEWIZojImRzpwn8YXbh7Y1zvRNuRm5LJizgMOqPaRX53H0ivsI17xDXBMcceqZsQ7PGBMjnTlPoq0jhV2qo59prGkmHFLe31NBfcBHbcYkPF67V4QxA1VH2iRuBr4BjBWRVRGTkoEPuisw04OqiqCmFBbMRXQeTd5GClP8VKeNYOLVt8Q6OmNMDHWkd9PfgZeAnwN3RYyvUtWybonK9KyaUgjVgB/CeKhq3oSfRr79/e+QOcLuFmHMQNZuklDVCpw70l0uIunAOCAeQERQ1Xe6N0TTHaL1aNrjX0d9aBtTzjrHEoQxplP3uL4euBXnPtQrca7d9CFwWveEZrpTS4+m3Ixcpm+bhr9pKGUlawjXv42oMPOiK2IdojGmF+hMw/WtwLHAdlU9Fecuc3u6JSrTI1p6NGVWKynVeaQWP0uIYgYNH08wJTXW4RljeoHOnHFdr6r1IoKIxKnqBhHJ7bbITLfK2jCIzPzDWfT6I1THDyexLo/6Y8aSXFHBpff+LNbhGWN6ic6UJPJEJA34J/CaiDwHbO+esEx3y8w/nITaoQAk1ecjtR9TvHUzsy+7En8gLsbRGWN6i85c4O8C9+WPReRNIBV4uVuiMj2iLljIBf97Ay+eMZ8tac1kjx7LxNmnxDosY0wv0qkL/LVQ1be7OhDTA5YtgNXPApBYdS7BOmH7lVfR1FxGozeVk7/yNcTTmcKlMaa/68jJdC0X9ou800zLsF3grw95ZvXDLG4sgUAis+qEQAhW7apga1Yq6k9j1FFTYh2iMaaX6ch5Esk9EYjpfktLjiBn19EEA0k0BlJRLWbVUVMZnLeM7MvtzGpjzBd15jyJu9sar6r3dF04pjtl7ppCQu0wRg4ZSlneGuKq1zGsaRNHnH4WZ849PtbhGWN6oc60SdREvI4H5gHruzYc09Uiz6w+ovp4kmrzmLryKcpWreHjscPx+uOZdclXYhylMaa36kzvpl9FDovIL4FXujwi06Uib0eaVjuDuJCAD7YMzaEq0MwJ8y8iMS091mEaY3qpQ+nKEsS5RIfpzVpuR6qDCeLBmxDHyEceZt2IbBrikpk2d36sIzTG9GKdaZNYzb7bl3qBbMDaI3q5rOIZZO6awqIh06iMLyElXMbad94gsaqILZPm44+Lj3WIxpherDNtEvMiXjcBxara1MXxmC7W0lgNkBIuY0jjJt574iWqU4aze/DkGEdnjOntOlPdVAxcCNwH/A64RUTsb2gfUBcs4ILbp3JY6XPs2L2Mmj3lvJ0xC0Taf7MxZkDrTJJYCEwGfg/8wX39aHcEZbpOYrWSXRJ2zqzeuZVdCcquwUeQPXY886cMj3V4xpherjPVTUeo6qSI4TdFZF1XB2S6VrBWCYQUfLA2Zwgqwvfuvp2UrOxYh2aM6QM6U5L4RET2nnElIjOAZV0fkjlkyxbAgrnOAyUUELz/eSeVcUrB6BMsQRhjOqwzJYlpwAcissMdHglsbOn1pKpHdXl05qAsfH0VBaXngnhJCQzF35DPQ7+8D3xJFI6aGevwjDF9SGeSxJxDWZCI5OC0awzG6Ur7gKr+VkQygKeA0cA24BJVLRcRAX4LnAPUAteo6ieHEsNAUbFtIilNw4hrLMLfkE99w1IS/cVsmTSfc6eNjnV4xpg+pDNnXG8XkaOBE91R76rqp51YVhNwu6p+IiLJwHIReQ24BnhdVX8hIncBdwF3AmcD49zHDODP7rNph9MOkcdJvg9YU1rB+iwPww7P5T/uvh6xHk3GmE7ocJuEiNwKPA4Mch+Pici3Ovp+VS1sKQmoahXOdZ+GA/OBR9zZHgHOd1/PBxaq4yMgTUSGdnR5A10oIIx6dCFvTJ2BaCOnXn2jJQhjTKd1prrpOmCGqtYAiMi9wIc4XWI7RURGA8cAS4DBqlroTirCqY4CJ4HsjHhbnjuuMGIcInIjcCPAyJEjOxtKv3bVfYs5YsdH7Bp8BEPH2e3IjTGd15kkIUBzxHAz+9+IqGMfIpIE/AP4tqpWRv67VVUVEY365jao6gPAAwDTp0/v1Hv7k4XPPk/BimoAUgLDiWvIZ+Tmf4PHy6T5V8Q4OmNMX9WZJLEAWCIii9zh84EHO7MwEfHjJIjHVfX/3NHFIjJUVQvd6qQSd3w+kBPx9hHuONOGirfLSanPJK6xiEBDJb76D0hvKOPEK67huNOt45kx5uB0puH61yLyFjDbHXWtqq7o6Pvd3koPAutV9dcRk54HrgZ+4T4/FzH+FhF5EqfBuiKiWsoAa9/N57OlxQD46rNIqt7JtII/4g/X8/rww0gbnMPUc+wqr8aYg9eRe1zHA18HDgdWA386yAv7nQBcCawWkZXuuO/jJIenReQ6YDtwiTttMU731804XWCvPYhl9mufLS1mV141WSOSCITyiKv/GD0rlc9KB1FT4OP8q67H5/fHOkxjTB/WkZLEI0Aj8C5Ot9SJwLc7uyBVfY/obRintzG/At/s7HIGmqwRSVxw+1QWn/FVAEbd/B5v3PZ1Rk+ZyNipx8U4OmNMX9eRJDFJVY8EEJEHgaXdG5I5kMgqpsLtlZQH4Mn7P+RqBQTeeewhmhtDnHaNdXk1xhy6jiSJxpYXqtpkB57YiqxiamgsJGfnO0z6dCmDdynbhgXZ9t5bHP/lS0kfald4NcYcuo4kiaNFpNJ9LUCCOyw4tUIp3RadaVNLFdOiM28gp6SW3SOSKR3koSgjm5T0wRx3/sWxDtEY00+0e8a1qnpVNcV9JKuqL+K1JYgY2zkoyFkvLiVr5mjqm3yces2NdktSY0yX6cylwk0v0FRaSv2GDWy/8ipGldQRp3VU/mkuH2wKMzZbOXy6Xd7KGNN1OnMynekFqosL8NaH+GTXWnSQsm28EloPKh5Ou+DsWIdnjOlnLEn0AZE9mup8g/AF8vjjFaMZ1biFU4tS2bJWOOmrXyP19C/HOFJjTH9j1U19QEuPJgB/qIBg1acsufYfPCbD8WxOInvkaKaefV6MozTG9EdWkugDmkpLSarezdSVL7Nr1XJ2DHIapt/fBNUNcO4Nt+D12aY0xnQ9K0n0AU27dxOurQWgINvLZ7khCu+byyfb4eiRMGz8hBhHaIzpr+zvZy8V2Q5RIen4fU3cMetmmvXrJGqIV5dBUhyceL5dwM8Y030sSfRSkWdW+xoK8Vd+zLbAc0CY2ZsHsatamP/dHxJnXV6NMd3IkkQv1bodYtsgYdLQYwlsDjJsSzLjZ55o50QYY7qdtUn0Uq3bITbnCg8VlHDqJ0HifXDaNTfGOEJjzEBgJYleJLIdotKTQUoQRv11Id99YArxWs+K7ZBfk8ics6eTmJYe42iNMQOBlSR6kcjzIVLCZQxv2rpvYl0S726NZ+zUY5l09Y9iFKExZqCxkkQvEtkOUbP6U0gNs/a/XyY+vZ6x63Lw+f2cecMtdp8IY0yPsSTRi+xth/BAeXozSyY18VFWNXF5GSRXxHHaLTeRlJEZ6zCNMQOIJYlexhMM7m2H2O4Xjok/ksM2l5MwIYcJs0+JdXjGmAHGkkSMRTZW7yGVQCiPSx6Ywk5fI6Mb/JywIYOqxDBX3/bfVs1kjOlx1nAdY5GN1YFQHnFVHwOQ0+Rnev4ESrZt4cwbb7HeTMaYmLAk0Qu03I50aP5vyaj8gKdvXMlvTnkC76ZaJp98BuOOnRnrEI0xA5QliRhr2rmZ+lXL2H7mMQwtDQMQqqvlpT/8muSsLE61k+aMMTFkbRIxVle+BwmF2UCI2kHCqglJyEN/oaKkmEvu/m/igsFYh2iMGcAsScRAZGN1TWAYPvJ58vbjWFdYyYS80SS98wbHX3g5IyYdEeNIjTEDnVU3xcD+jdX5JFcto3b7jXg2Xcrk9Z8zLHcSMy+8LMZRGmOMlSRiorSulLLgbj6Y/BrnvfwRcWFBwlczp/Tf+Hxe5n7rO3i83liHaYwxVpKIhYbCPDLzqrjsV0sZXaKkhj3clLCOpMoC5n3jVlKyB8U6RGOMAawk0WMi2yF8ocHEh/KYQICatCbWjZpI0UvPc8zZ5zL++NkxjtQYY/axJNFD3n97DXVFSl1qGZmhKvwNy7njoj+xc/sOzt/2FCMOz+Xkr34t1mEaY8x+LEn0kIbiIoJV9Uz4/EGy86qozfDgab6SOcWv4PP7mXfbnXh9/liHaYwx+7E2iR4SrGwgPqTkVNeRmtJE2sgEvtq0hMSaUi749ndJybJ2CGNM72MliW4U2Q4RCgwDzUPPSqWeVDbqNLa9/zazL7uKMcdMj3GkxhjTNitJdKPW50MkVi9n8vffI2neH9i++nPGzZjFcedfHOMojTEmOitJdKOmnZtJKt/D1DV/pqykjsJsD+VFBbzw23vJGD6COd+4zS7/bYzp1XqsJCEiD4lIiYisiRiXISKvicgm9zndHS8i8jsR2Swiq0Rkak/F2ZXqyvfQ5F6Xadsg4ePcZBb94icgwvzv/oBAfEKsQzTGmAPqyeqmh4E5rcbdBbyuquOA191hgLOBce7jRuDPPRRjl2pEqQ/A3VdM4McXjCPUfDiVpcXMv/37pA8ZFuvwjDGmXT1W3aSq74jI6Faj5wOnuK8fAd4C7nTHL1RVBT4SkTQRGaqqhT0T7cH7+Ae/ZvPOJAAa/cPxh/IZ3XA7Yz57gYzyVXzpltsZMdEu3GeM6Rti3XA9OOLAXwQMdl8PB3ZGzJfnjvsCEblRRJaJyLLS0tLui7SDNu9MotLvlBL8jfmUJq/m1szPyS5cxfEXXs6kE0+NcYTGGNNxsU4Se7mlBj2I9z2gqtNVdXp2dnY3RNZ5KY0FXP7IjTxwyqt8OiiPD599giNOPZNZF18R69CMMaZTYt27qbilGklEhgIl7vh8ICdivhHuuF7p4/ueY/O6GgAq/MPwN+YzY8GFZBUVM3VdGmOnHsuZN9xiPZmMMX1OrEsSzwNXu6+vBp6LGH+V28vpeKCiN7dHbF5XQ6WmAk4VU37SMgbtDnPSqjSastOYd+uddulvY0yf1GMlCRF5AqeROktE8oAfAb8AnhaR64DtwCXu7IuBc4DNQC1wbU/FeTAam+rwhat49dS/skHrmFiayBmrckgdNpJLfvRz/PHxsQ7RGGMOSk/2bro8yqTT25hXgW92b0QHL/JyGwA1vsH4QnnUhJrILQ8y8dNBJA/O5uIf/oxgSmoMIzXGmEMT6zaJPmndi2soK1dSwmUAJFVX4q9fTlLFD8ld9QT+pGQu/uFPSUxLj3GkxhhzaCxJHISmkkKS6hqYtesvAKzXOj49LJGp6yqJz0jlkrt/TnJGVoyjNMaYQ2dJ4mA0h/D4lFFXOOdD/KS0lNyV2SRkpXDx3T+zy34bY/oNSxIdFNnNtdI/jKSGfC4N/YDk8u2MX/UYdXHCJT+2EoQxpn+xJNFBLd1cU6SCuFA+RUnLaSx/ifErm6hMbGLtSWmWIIwx/Y4liY5qDpHSnM/lpzzHVaEthErSmbUyk7osH0VnZnLFxHmxjtAYY7qcJYkD2K+KyZNFSuNO1hRUkFGZTc62VEZPmcZ5t33PzoMwxvRbliQOoHUVU0HiJ3xYmsGowjCFo/x8+7s/wOvzxzpMY4zpNpYkWok8Ua7Sk0FKuIzL/3oFX/n9VEauzmRUYRxF0xI5bv6FliCMMf2eJYlWIk+US6qpZnj9CvJ//ThHrhiEp8nDuf/xPcbPOCHWYRpjTI+wJNFK0+7dJNXWMsvzDgXFK/hgVJAlS7Ooj1fWTK3nDksQxpgBxJJEK41NdTT6mrln9noS1mYxpjCRomwPS46KY+aos2IdnjHG9ChLEq00hpsIhXdz2HuZ+Bu8lI9J4Jf//RTiifVV1Y0xpudZkgCe/e7dVJQMRzVMrRbRGPoEjc9i/dTzOePk4yxBGGMGLEsSQEXJcBo8CYTqXqFZK/AGUvnOn+8nLpgY69CMMSamLEkADU07qGtaTaM/nm25F3LSl063BGGMMViSAMDnDeHXVG7+4+8JpqbFOhxjjOk1LEkANy7831iHYIwxvZK1yBpjjInKkoQxxpioLEkYY4yJypKEMcaYqCxJGGOMicqShDHGmKgsSRhjjInKkoQxxpioRFVjHUOXEZFSYPtBvj0L2NWF4fQFts4Dg63zwHAo6zxKVbPbmtCvksShEJFlqjo91nH0JFvngcHWeWDornW26iZjjDFRWZIwxhgTlSWJfR6IdQAxYOs8MNg6Dwzdss7WJmGMMSYqK0kYY4yJypKEMcaYqCxJACIyR0Q2ishmEbkr1vF0BxHJEZE3RWSdiKwVkVvd8Rki8pqIbHKf02Mda1cSEa+IrBCRF9zhMSKyxN3WT4lIINYxdiURSRORZ0Vkg4isF5GZA2Ab3+b+pteIyBMiEt/ftrOIPCQiJSKyJmJcm9tVHL9z132ViEw9lGUP+CQhIl7gj8DZwCTgchGZFNuoukUTcLuqTgKOB77pruddwOuqOg543R3uT24F1kcM3wvcp6qHA+XAdTGJqvv8FnhZVScAR+Ose7/dxiIyHPh/wHRVPQLwApfR/7bzw8CcVuOibdezgXHu40bgz4ey4AGfJIDjgM2qulVVQ8CTwPwYx9TlVLVQVT9xX1fhHDyG46zrI+5sjwDnxybCriciI4C5wN/cYQFOA551Z+lv65sKnAQ8CKCqIVXdQz/exi4fkCAiPiAIFNLPtrOqvgOUtRodbbvOBxaq4yMgTUSGHuyyLUk4B8qdEcN57rh+S0RGA8cAS4DBqlroTioCBscorO7wG+AOIOwOZwJ7VLXJHe5v23oMUAoscKvY/iYiifTjbayq+cAvgR04yaECWE7/3s4tom3XLj2mWZIYYEQkCfgH8G1VrYycpk5/6H7RJ1pE5gElqro81rH0IB8wFfizqh4D1NCqaqk/bWMAtx5+Pk6CHAYk8sVqmX6vO7erJQnIB3Iihke44/odEfHjJIjHVfX/3NHFLUVR97kkVvF1sROA80RkG04V4mk49fVpbrUE9L9tnQfkqeoSd/hZnKTRX7cxwBnA56paqqqNwP/hbPv+vJ1bRNuuXXpMsyQBHwPj3N4QAZxGr+djHFOXc+vjHwTWq+qvIyY9D1ztvr4aeK6nY+sOqvo9VR2hqqNxtukbqvoV4E3gIne2frO+AKpaBOwUkVx31OnAOvrpNnbtAI4XkaD7G29Z5367nSNE267PA1e5vZyOByoiqqU6zc64BkTkHJz6ay/wkKr+LMYhdTkRmQ28C6xmXx3993HaJZ4GRuJcZv0SVW3dQNanicgpwHdUdZ6IjMUpWWQAK4CvqmpDLOPrSiIyBaehPgBsBa7F+TPYb7exiPwEuBSnB98K4HqcOvh+s51F5AngFJzLgRcDPwL+SRvb1U2Wf8CpdqsFrlXVZQe9bEsSxhhjorHqJmOMMVFZkjDGGBOVJQljjDFRWZIwxhgTlSUJY4wxUVmSMOYgiEimiKx0H0Uiku++rhaRP8U6PmO6inWBNeYQiciPgWpV/WWsYzGmq1lJwpguJCKnRNy74sci8oiIvCsi20XkyyLyPyKyWkRedi+TgohME5G3RWS5iLxyKFfsNKarWZIwpnsdhnPdqPOAx4A3VfVIoA6Y6yaK3wMXqeo04CGg353xb/ouX/uzGGMOwUuq2igiq3Eu+/KyO341MBrIBY4AXnOupoAX55LXxvQKliSM6V4NAKoaFpFG3dcIGMbZ/wRYq6ozYxWgMQdi1U3GxNZGIFtEZoJzOXcRmRzjmIzZy5KEMTHk3jL3IuBeEfkUWAnMim1UxuxjXWCNMcZEZSUJY4wxUVmSMMYYE5UlCWOMMVFZkjDGGBOVJQljjDFRWZIwxhgTlSUJY4wxUf1/S/7kW7SAD6AAAAAASUVORK5CYII=\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "def simulate_n(n):\n", + " values = np.zeros(len(times))\n", + " for i in range(n):\n", + " values += model.simulate(params,times) / n\n", + " plt.plot(times, values, label=r'$n=%s$' % n)\n", + " \n", + "for i in range(5):\n", + " values = model.simulate(params, times)\n", + " plt.step(times, values)\n", + "\n", + "simulate_n(1000)\n", + "plt.title('Stochastic logistic growth across different iterations')\n", + "plt.xlabel('Time')\n", + "plt.ylabel(r'Population ($C(t)$)')\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "By running the model 1000 times and averaging the results we can see that indeed we converge on the deterministic mean:" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYkAAAEGCAYAAACQO2mwAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nO3dd3wUdf7H8ddnNwkhJCSUgEhAOtJbQBBBFBEBBQ4VESwoWO5QOfXuJ5Y7LFfksJcT8UARCyqKYEOKoKggEKoCAiJg6DWEkLr7+f2xAwZIIIEks9l8no/HPrLznZnd9zJhP5n5znxHVBVjjDEmLx63AxhjjAleViSMMcbky4qEMcaYfFmRMMYYky8rEsYYY/IV5naAolS1alWtU6eO2zGMMaZUSUpK2quq8XnNC6kiUadOHZYuXep2DGOMKVVEZEt+8+xwkzHGmHxZkTDGGJMvKxLGGGPyZUXCGGNMvqxIGGOMyVeJFgkR2Swiq0VkhYgsddoqi8hsEdng/KzktIuIvCAiG0VklYi0Lcmsxhhj3NmTuERVW6tqojM9Cpirqg2Buc40QC+gofO4HXilxJMaY0wZFwzXSfQDujnPJwHzgQec9jc1MJb5IhGJE5EaqrrDlZTGmFLL51eyMw6Tk5WJLzsTny+HnOxsstVDVlR1fH4/nv2b0KzDqM+H35+D3+cjJ7wCh2Mb41cleucSJDsN9fvw+3yo+sksV5n9ldvgV6V68iw8vgxU/ajfD+rncPma7K7SHr9C3a1TEX82qIL6UVUOVGjA9sodUL+fFlsngyqq/kBoVXZVbEpyXAe8/kxab30zsG5gJqiyNbYd22LboUD386vRqlZckf/blXSRUGCWiCjwqqqOB6rn+uLfCVR3ntcEfsu1brLTdlyREJHbCexpULt27WKMbow5KzmZ+NNTSE9LISPtEOlph8hKT2NX1U4cycohcttCovavQbPT0ax0yMkg26/MrPEnMnN8dNz5LvUOL8Prz8KjWYT5s0iVaB6JfpysHD9/OfIM7XOW48VHONmEqY+tVKdH1lhU4b2Ix7nAs+64SCv99eiX9Q8APot4kGae468pW+hryvXZjwAwP+Je6nh2HTd/tq8tf8r+CwBLyv2deEk5bv7Hvgt5IDsSgDXlniRKMo+b/3ZOd57IiUbw82vksyf9k72a04exObHEcITVkeNOmr9w036e91UEID6mXLEUCSnJmw6JSE1V3SYi1YDZwN3ADFWNy7XMAVWtJCKfAk+q6rdO+1zgAVXN95LqxMREtSuujSkmqpCZiu/wXg4f2Enawd1kHNrLL1W7sz/LS8zWudTePhNv1iHCcw4TkZNKpC+NWyu8yJ7MMP6Y+TpD5dOTXrZexlv48fCPsAncEDb3WHu6RpBCND1kHJHhXkbqW1zgX0m2RODzROCTcFLDKvFGtVFEhHnofuhjamZvRj3h4AlHveGkR1RmZcIQwr0eGu+bTUz2fsTjBW8Y4gknu3wV9p57CV6Ph+p7FxGecwTxehFP4OGPrEx61RZ4PEL0/p8I02zE48XjLKPlYvHF1sIrQsShzQh+PJ4wRLyIRyAiCqKqIgJhR/YgIogIHo8XEQ+ERSAR0YAi2YfxeDwIgngk8BreMPCEISii/mPri0iu5eSsN62IJOXqAjhOie5JqOo25+duEZkGdAB2HT2MJCI1gN3O4tuAWrlWT3DajDHFwJe6m5QNC0nds5XsA9vxpe7Am7abKVXuYm1mZTrt/ZC7MsbjBWKdB8CtmU+zWWtwo3cFLcKSSJNoUr0VyPJWJysymkbx5WlSPg7NuZK52U3wlIvGGxlDWGQ04ZHRvF2zAxUiw4nWNmwLU8qVjyYyMoryEWGc4xFWH0t4WZ65uxx71i7P+b+v1fA0/wIDTj27XtdTz6/e4tTzY093pKPyaea7o8SKhIhUADyqmuo8vxx4HJgB3Aw86fyc7qwyA7hLRKYAFwAp1h9hzJnTtH2krJ1L6vaNZO3dhDdlK1Hp2/lv9F3MPtKIlqnf8Er4s1QG/CrsoyJ7iGPDkR2kx8VyoGoiX/jCkApVCYuuSkTFeCIrxvNitbrERUdRqUJPKkQ8d9Jftr9/ibc6TcKiP1Rizl5J7klUB6Y5v0BhwDuqOlNElgDvi8gwYAsw0Fn+c6A3sBE4AtxSglmNKZ3S9pK9/UcObFlF+rY1ePavZ1aFvkzLbEfUnpW873mYOGCfxpCs8WwKSyCTSDrUrUz9Cr2ZGd6BmPgEKlVL4JxK0ZxfPpxJntxf+te59cmMS0qsSKjqJvL4U0JV9wHd82hXYEQJRDOm9FGFg1vI3LqMTRkV+D67EVu3bOKxDQMIB6oBhzSKjXoum46kUblGORrVvpBPIt+jUkIjap1TnaZx5Qn3erjU7c9iglownAJrjCkA9fs5MGsMGZu+p+K+VUT7DlIOWOfrzBPZI6gWHUHDSiMIr96Y2DqtSEioQ5NqMfwzwut2dFOKWZEwJhhlp6NbfyBlzVy2793Py+WGsfjX/byeOYVyZLPY05oDlVoSXjuRGg3asvi8eKpVjAR6uJ3chBgrEsYEkfSVH5H2/QRidy8mXLOIVg97/M1ZFnWAzvWr8ON5H9GuXjX6xkfj8Zz9qY/GnI4VCWPc4vfDtiRSV07n44qDmbn+EB22fEof2cQXchn7a1xEtWbd6NikDt9XiSqS8+GNKSwrEsaUJFXYvpxDS99D1kwjJnMXkeplelYlDsYnkt7xPvacX4PrzqtERJgN0mzcZ0XCmJLg95Oa5WPBN3Po/f0gItXLN/6W/Bh3A7Gt+vKfVg2oFx/tdkpjTmJFwpji4vehG+eQsmA8aw5XYNjewaRn+7iz0v3Et72Ky9s15rLKUW6nNOaUrEgYU9TSD5K55A2yv3+V6IztZGssP2hv+rc5l+va16ZVQh/rXzClhhUJY4rQjpR0tr9zP+12TWWZrynzYm+ibueB3Nb2PKLL2X83U/rYb60xZ2vfL6TOGcMbmZfw/LpYamgXeje4nJ49LufBWnG212BKNSsSxpypA1s4/OUTRK37kHD1stMfyw0dhzLsom7Usr4GEyKsSBhzBtJnjibih5cI8wtvaC9S2vyJkZclUi0m0u1oxhQpKxLGFFROFpnq4fXvt5C5cCfx/q5saX4Xw3p1dobEMCb0WJEwpgB0/Zccmf4X/p09iLcOtebS84fzYK/zGVw9xu1oxhQrKxLGnMqh7aRN/wsVfvmMHf5zyYyJY/KwDnRpGO92MmNKhBUJY/KRk/Q2/s//ijcni+cZRJWef+HfHesT5rXhMkzZYUXCmDys23mIj7/6lc5ZdZlTbxR/uroH1a3fwZRBViSMOUoV//K3+XbdNob/1JKKkR1oe90wHmtew+1kxrjGioQxABkppE8bSfmfp+HztaJbo0v599UtqRJdzu1kxrjKioQxyUmkv3sT4Wnbec5/HTWveohX259nV0obgxUJU8b5Dm7DP6EX+/wxPB39JCNuHkyDanZaqzFHWZEwZZMqKRk5jJy2jYqZw6nSshf/vvpCIsO9biczJqhYkTBlT+ou0t8ewuiD/fj2UH0e63cHQy44z+1UxgQlKxKmbNn5IxlvXose2Qee7rx7e0fa16nsdipjgpYVCVN2rJ9F9ns3cyAnkn/EjOGh4YOoGVfe7VTGBDUrEqZM0OQk9J1B/OyvxcTaTzLmxh52EyBjCsD+l5iQ5/MrDy/0UiF7ENLuZsb274DXY6e3GlMQViRM6FIlZ94YHt/anCnrhHsuHcm9PRrZ9Q/GFIIVCROa/H5yPvsrYUn/o2L2tTzc+yFu61rP7VTGlDpWJEzo8fvJmTGSsBVv8mpOH2pc9QhDOtZxO5UxpZIVCRNa/H5yPh5B2Kp3eCmnP1WvepxBdg2EMWfMioQJKRlHDpG8bhmf5gzg3H6PM7B9LbcjGVOq2d1TTGhQJTsrk7s+3EifQ6M4t99jViCMKQK2J2FCgv+rf/JL0lcs2H83D/drw8D2td2OZExIsD0JU+rponF4FoxlxaEY/tyzOTd1quN2JGNCRokXCRHxishyEfnUma4rIj+IyEYReU9EIpz2cs70Rmd+nZLOakqBn6bBzFHM8rUj+aJ/8sdLGrqdyJiQ4saexEhgba7pMcCzqtoAOAAMc9qHAQec9med5Yz53a8L8E29jaX+hixo9R/u79nM7UTGhJwSLRIikgD0Af7nTAtwKTDVWWQS0N953s+ZxpnfXexSWZPLV1tz+C7nfN6t9x8eHdDOrqQ2phiUdMf1c8D/AUdv/VUFOKiqOc50MlDTeV4T+A1AVXNEJMVZfm/uFxSR24HbAWrXts7KMiE7nSXb0rljVjpta49l0hAbi8mY4lJiexIiciWwW1WTivJ1VXW8qiaqamJ8fHxRvrQJRtkZZE64kjWT/kytSlGMvzHR7iZnTDEqyT2JzkBfEekNRAIVgeeBOBEJc/YmEoBtzvLbgFpAsoiEAbHAvhLMa4KNKtnT/kS5nUtZRVcmDG1PbFS426mMCWkltiehqg+qaoKq1gEGAV+p6hBgHnCNs9jNwHTn+QxnGmf+V6qqJZXXBB//dy8QvuZDns4ZyNU3jKBu1QpuRzIm5AXDdRIPAPeJyEYCfQ4TnPYJQBWn/T5glEv5TDD45SuY8yif+i6gWp+HubBBVbcTGVMmuHLFtarOB+Y7zzcBHfJYJgO4tkSDmaC1bOM21F+fH5o/zuMdbcA+Y0qKDcthgt7G3anc+F08DeKf5b0B7e1UV2NKkBUJE9Syp49kxprylI/ozit2JpMxJS4Y+iSMyZMum0z48jcIT9vBi9e35dy48m5HMqbMsT0JE5x2/ojv0/tZ7GtKWPeH6FS/ituJjCmTbE/CBJ/Mw2S8eyP7feX54LxHuePiRm4nMqbMsiJhgk7aT58TlrKZxyPuZ/TgS/DYkBvGuMYON5mgoqrc91M9NmQ/w7PD/kBcVITbkYwp06xImOBxYAszFy7jy58ieKTPxbSqFed2ImPKPCsSJjj4cjgy5RY67VxHr4bvcGvnum4nMsZgfRImSGR98zRRu5J4Ouw2/nFdJ+uHMCZIFLpIiEgFEbErmkzR2bES79djmOHrRK/r76ZKdDm3ExljHKctEiLiEZHBIvKZiOwG1gE7RGSNiIwVkQbFH9OErOwMDk8Zxh6tyPp2o23gPmOCTEH6JOYBc4AHgR9V1Q8gIpWBS4AxIjJNVd8qvpgmVB3IUCandmRHhfqM7nPSOI/GGJcVpEhcpqrZJzaq6n7gQ+BDEbE7v5hCU7+fh2esYXZmbz4e3tnGZTImCJ32cNOJBSKvPom8iogxp5R1hP0vdefIT19wb49GNDs31u1Expg8WJ+EcUX6l49RZf8yEqpW4vYu9dyOY4zJR0HObpoH1CfQJ3GOqtZS1WrARcAiAn0SNxRjRhNqfltMuaRXecd3GTcNuYkwr52JbUywKnCfhIjUOdppDdYnYc5QdgaH37+Dg1qF1C5/o1H1GLcTGWNOoTB9Eh+dOE9EOp6wjDGnlLbsA6JTNzGu4j3c2r2l23GMMadx2j0JERkItAViRKQJ8HOuPYrxgP1PNwX2t83N2ZL9BE9cP5RwO8xkTNAryP/S74A1QCXgGWCjiCwTkU+B9OIMZ0KI38ei5Sv4aPl2OnXtSdNzK7qdyBhTAKfdk1DVbcCbIvKLqn4HICJVgDoEznQy5rQyF75K69mjuaTyM9x1qZ0QZ0xpUZDDTaIB3x1tU9V9wL4TlymmjKa0S9kGc5/gB//5jLi2t100Z0wpUqBTYEXkbhGpnbtRRCJE5FIRmQTcXDzxTCg4OO0+1JfNsuaPkFjX7lVtTGlSkFNgrwBuBd4VkbrAQSAS8AKzgOdUdXnxRTSlmW/dF8Rtnsl/w4YwvO+lbscxxhRSQfokMoD/Av91roeoCqSr6sHiDmdKv+WLFxDtr0X9AQ8QE2mX0xhT2hSkT+JvwH5Vfdm5HmJH8ccyoWDbwXRu2tiVi+r04dUWtU+/gjEm6BTkcNNAoP2JjSIyHIhX1X8XeSpT+u3dyOSPvsOv1fhb/zaI2J3mjCmNClIksp1DTieaDCwDrEiY46my/4O7uGPnaqp0+5JalaPcTmSMOUMFObspS0RqnNioqpmADcdhTpK1+iMq71rIW+WHMLRbU7fjGGPOQkGKxNPAdBE5L3ejiFQD7NoIc7zMw2R+Ooqf/OfR7ur7begNY0q5gpzd9IGIRAFJIrIIWEGguFwLPFq88UxpkzLr38Rm7ebLOi9xX8Pqbscxxpylgl5xPUlEPgL6A82BNGCwqi7JtYztVRg+3+RDtSeDr77W7SjGmCJQkI7reSLyITBdVScfbTx6xTWBq63nAW8UT0RTWsxdu4sHd3Thod63cU5spNtxjDFFwK64NkUia8NXfPfxAupV7crQC+u6HccYU0RK7IprEYkEvgHKOe87VVVHO4VnClAFSAJuVNUsESkHvAm0IzCY4HWqurkw72lKSE4WRz4ayZD0bLr1H05EmHVWGxMqCvW/WVWzVXXHGQ7JkQlcqqqtgNbAFc6d7cYAz6pqA+AAMMxZfhhwwGl/1lnOBKHUr18kLn0rM2rcTdcmNd2OY4wpQiX2J58z3PhhZzLceShwKTDVaZ9EoHMcoJ8zjTO/u9hlu8EndRfh343lK39b/nDtULfTGGOKWIkeFxARr4isAHYDs4FfgIOqmuMskgwc/VO0JvAbgDM/hcAhqRNf83YRWSoiS/fs2VPcH8GcYO+Mh/H4stjY5kHqVK3gdhxjTBEr0SKhqj5VbQ0kAB2A84vgNceraqKqJsbHx591RlNwfr8yYVdDXvEOYXBvGwbcmFBUkLObAHA6kq8mcNvSY+up6uOFfVNVPSgi84BOQJyIhDl7CwnANmexbUAtIFlEwoBYct0Nz7hvxsrtvLK7OWOvGUx0uQL/KhljSpHC7ElMJ9BPkEPgYrqjjwIRkXgRiXOelwd6AGsJXGNxjbPYzc77AMzg9zveXQN8ZRfsBY/M1dPZ8ek/aHNuFFe3TXA7jjGmmBTmz78EVb3iLN6rBjBJRLwEitP7qvqpiKwBpojIP4DlwARn+QnAZBHZCOwHBp3Fe5uilJ1BxqcP0C07nMQrH8fjsfMJjAlVhSkS34tIC1VdfSZvpKqrgDZ5tG8i0D9xYnsGgfGhTJA5NO95YjN3MKn209xTz/qBjAllhSkSFwFDReRXAtc8CIEzW1sWSzITnA7vJmLRs8zxJzLg6sFupzHGFLPCFIlexZbClBp7P/k7sb4sNrd9gMsq2c2EjAl1BS4SqrpFRFoBXZymBaq6snhimWCkqozd3Z7K3mhG9LJTXo0pCwp8dpOIjATeBqo5j7dE5O7iCmaCzyerdvDejnOo22uknfJqTBlRmP/pw4ALVDUNQETGAAuBF4sjmAkuWeu+hOnjaX/ObVzdzk55NaasKEyREMCXa9rntJlQ58shdcYoWuakcX/vVnjtlFdjyozCFInXgR9EZJoz3Z/fr2kwIezQwolUObKJD2o8yp2NargdxxhTggrTcf2MiHwNdHaabrGbDZUBGYeQef9iib8xV1xzm9tpjDElrFC9j6qaRODGQKaM2PPlWOJ9B1jddCy3xke7HccYU8JOWyRE5FtVvUhEUgnc/+HYLAIX01UstnTGVarK49vbU0OGM6Jv/9OvYIwJOac9BVZVL3J+xqhqxVyPGCsQoW3u2t18ssVLzR53ExsV7nYcY4wLCnOdxEm3D82rzYSG7OQVxH04kE5V0hh8QW234xhjXFKYocJ75NFmQ3WEIlX2ffQX6uVs5PYebQn3lui9qYwxQeS0//tF5I8ishpoLCKrcj1+Bc5oRFgT3I789Dnn7F/C9Lib6NaqvttxjDEuKsjZTe8AXwD/Bkblak9V1f3Fksq4x5fDkc8eYqf/HBKvvh8Ru3DOmLKsIB3XKaq6WVWvBw4B1YHzgOYi0rW4A5qSdWDhJKqmb+brWiNoUbuq23GMMS4rzD2uhwMjCdyHegXQkcDYTTYcaAj595bzifQN585rhrsdxRgTBArTIzkSaA9sUdVLCNxl7mCxpDKuWL5lP++vPkjFi27nXLtXhDGGwhWJDOeWoohIOVVdBzQunlimpGnKNipO7sHFFbZyZzfrrDbGBBRmWI5kEYkDPgZmi8gBYEvxxDIlLfmjR0jI3sw1l7S2e0UYY44pzAB/f3CePioi84BYYGaxpDIlKjN5JTW3TOPDyD8woGtHt+MYY4LIGf3JqKpfF3UQ4xJV9nz4VypoBRL6/s3uFWGMOU5BBvg7OrBf7m+Po9M2wF8pl7J2LgkHfuCdKiMY3Kye23GMMUHmtEVCVWNKIohxx1M/VyU15y7uuvZ+t6MYY4JQYa6T+Hte7ar6eNHFMSVpw84U3lmynRsuuJ4GNSq7HccYE4QKcwpsWq6Hj8DgfnWKIZMpCRmHiJpwEVdFLGXkZY3cTmOMCVKFObvp6dzTIvIU8GWRJzIlYusn/6J29lYuTmxN5QoRbscxxgSpsxkDOorAEB2mlMnZv4XqP01gtrcrva/o43YcY0wQK0yfxGp+v32pF4gHrD+iFNrywShqqhJxxWOUC/O6HccYE8QKc53Elbme5wC7VDWniPOYYpby2xrq7/icjysOol9iG7fjGGOCXGGKxC7gT8BFBPYoFojIuKPjOZnS4ellfjZmPcyj195o94owxpxWYfok3gSaAS8CLznPJxdHKFM81m3by1uLttDggt40ql3D7TjGmFKgMHsSzVW1aa7peSKypqgDmeKh2elUfP1i7ih3GbdfNtbtOMaYUqIwexLLROTY6G8icgGwtOgjmeKwYcZYzs1JpnW7C6lkp7waYwqoMHsS7YDvRWSrM10b+PnoWU+q2rLI05kikXFgOwmrX+b7sAvo3vtat+MYY0qRwhSJK87mjUSkFoF+jeoEOr7Hq+rzIlIZeI/A1dubgYGqekACvarPA72BI8BQVV12NhnKqk3vjaKBZhPR65+Eec/m0hhjTFlT4G8MVd0CxAFXOY84Vd1y9FGAl8gB7nf6NToCI0SkKTAKmKuqDYG5zjQEhv1o6DxuB14paFbzu707tlJ/x+d8FXc1ie3aux3HGFPKFLhIiMhI4G2gmvN4S0TuLuj6qrrj6J6AqqYCa4GaQD9gkrPYJKC/87wf8KYGLALiRMROySmkJ789wJU5Y2gy0K57NMYUXmEONw0DLlDVNAARGQMsJHBKbKGISB2gDfADUF1VdzizdhI4HAWBAvJbrtWSnbYdudoQkdsJ7GlQu3btwkYJaSvXb2JqUjJ3XHwh59W0+mqMKbzCHKAWAqO/HuXj+BsRFexFRKKBD4E/q+qh3PNUVfl96I8CUdXxqpqoqonx8fGFjROy/JlpnDvlckZHTeXuSxu6HccYU0oVZk/ideAHEZnmTPcHJhTmzUQknECBeFtVP3Kad4lIDVXd4RxO2u20bwNq5Vo9wWkzBbBm6uM09++hwUX9iS53RnepNcaYQnVcPwPcAux3Hreo6nMFXd85W2kCsNZ5raNmADc7z28Gpudqv0kCOgIpuQ5LmVM4tHMTDTZM5NtyF3NR975uxzHGlGIFucd1JHAn0ABYDfz3DAf26wzcCKwWkRVO20PAk8D7IjIM2AIMdOZ9TuD0140EToG95Qzes0z6bcp91FOIHzDGxmcyxpyVghyHmARkAwsInJbaBPhzYd9IVb8l/z6M7nksr8CIwr5PWbf+l42ce2AJ39a4mR6Nm7gdxxhTyhWkSDRV1RYAIjIBWFy8kcyZUlUenr2bPZ4X+HhwD7fjGGNCQEH6JLKPPrH7RwS3r+bPJWnzPv7Yuz1xFSu6HccYEwIKsifRSkSOnqoqQHlnWggcFbJvoyCQumcrnb4ewjOVetO33ZWnX8EYYwrgtEVCVe3+lqXA5nfupZH6OP+qe/F4rLPaGFM0bLS3EPDrD5/Q4sAcvj3nRs5vaoPxGmOKjhWJUs6XlU65L/+PrZxD4g02PpMxpmhZkSjlZnyzhPQcJfnCfxAbE+N2HGNMiLHxGkqxXYcy+PuCdNrUmsikHp3cjmOMCUG2J1FaqfLVW0/i8R3hsT+0siurjTHFwopEKfXjzNe4fvezPNNkI3WrVnA7jjEmRFmRKIUO799JzR+eYI33fLoMvNftOMaYEGZFohTaOHkkFTQN6fs8EeHWrWSMKT5WJEqZn7+fTusDM1lU40aatOrodhxjTIizIlGKZGT7+Nd3qXzuvYTEm/7pdhxjTBlgRaIUeXbOer7eV4mYQa8RFRXtdhxjTBlgRaKU2LhkFq2/v4db21SkS0O7l7cxpmRYr2cpkJmeSvkv7qGl10/nXs3djmOMKUNsT6IUWD35AWr6d7D7krFUrFjJ7TjGmDLEikSQ27BkFm23vcP3lfrS5uJ+bscxxpQxViSCWHpmDllfPMJ2qUazoS+4HccYUwZZkQhiY778mRuP3MuePhOJjbXDTMaYkmcd10FqadISJn+/gxsvbEmb9s3cjmOMKaOsSAShlL3bqfvJNbwUk8glvd53O44xpgyzw01BRv1+tk68hWhNo/6VfyEy3G4xboxxjxWJILN86hhaHFnEkkb30qiV3UjIGOMuKxJB5Le1i2n209MsK9eBCwc96HYcY4yxIhEsMnN8PDVzLT9KA2oOfR2P1zaNMcZ99k0UJP712Vqm76rK3ms/pnqNBLfjGGMMYEUiKKyc8SINlozmtgsT6NnsHLfjGGPMMXYKrMuS1y6mcdJjUL4Z113R1O04xhhzHCsSLjpyaD98MJRDEk21oZOJiAh3O5IJEdnZ2SQnJ5ORkeF2FBNEIiMjSUhIIDy84N81ViRcon4fG18dQhPfTlZf9iZta9Z2O5IJIcnJycTExFCnTh1ExO04JgioKvv27SM5OZm6desWeD3rk3DJu5/Npv7hJBY1up+2Xa50O44JMRkZGVSpUsUKhDlGRKhSpUqh9y5tT8IFs9fs4qHvfKxvNpnR1/dwO44JUVYgzInO5HfC9iRK2KY1S/h+yn9omRDLqEGXIR7bBMaY4FVi31AiMlFEdovIj7naKovIbBHZ4Pys5LSLiLwgIhtFZJWItC2pnMVpz67fiPxgMCM8Uxl/bT7B4VwAAA/RSURBVH0bl8kYE/RK8s/YN4ArTmgbBcxV1YbAXGcaoBfQ0HncDrxSQhmLTfqRNPa+di2V/AdJ6T+Zc6rXcDuSMcacVokVCVX9Bth/QnM/YJLzfBLQP1f7mxqwCIgTkVL7rer3+Vn9yk00yVnL+gvHUr91V7cjGVOq3HrrrVSrVo3mzZsf1z5z5kwaN25MgwYNePLJJws071TrmJO5fUC8uqrucJ7vBKo7z2sCv+VaLtlpO4mI3C4iS0Vk6Z49e4ov6RlSVSa9/x4dUuewpP5dtOo51O1IxpQ6Q4cOZebMmce1+Xw+RowYwRdffMGaNWt49913WbNmzSnnnWodk7egObtJVVVE9AzWGw+MB0hMTCz0+sXttQWb+NfKiuS0/h/DB17tdhxTBj32yU+s2X6oSF+z6bkVGX3V6e+YOGDAAJo2bco333zD5s2bmThxIpdddlmh369r165s3rz5uLbFixfToEED6tWrB8CgQYOYPn06TZs2zXdet27d8l3H5M3tPYldRw8jOT93O+3bgFq5lktw2kqVxZ+9zvyZU+nTsgbDBl5jZzKZMmf16tXExcXxzTff8Pzzz/P2228fN79Lly60bt36pMecOXNO+9rbtm2jVq3fvyYSEhLYtm3bKeedah2TN7f3JGYANwNPOj+n52q/S0SmABcAKbkOS5UKy+d/ROvFf+Gh6OY0vPZ+PB47Z924oyB/8ReHI0eOkJKSwr333gsEhgqJi4s7bpkFCxa4Ec0UQokVCRF5F+gGVBWRZGA0geLwvogMA7YAA53FPwd6AxuBI8AtJZWzKKz+7nPOn3cH28JqUfuPH1Au3O1abEzJW7NmDe3atcPrDZzqvWrVqpM6nrt06UJqaupJ6z711FOnPSxVs2ZNfvvt967L5ORkatasecp5p1rH5K3Evr1U9fp8ZnXPY1kFRhRvouKxZslX1J11K3u81ah0x2dUrFTN7UjGuGL16tW0bt362PSqVavo16/fccuczZ5E+/bt2bBhA7/++is1a9ZkypQpvPPOO6ec17hx43zXMXmzg+RFKGnLAZI+fY1DnlgqDP+MuGr2F4opu04sEj/++ONJexIFdf3119OpUyd+/vlnEhISmDBhAmFhYbz00kv07NmTJk2aMHDgQJo1Cxxay2/eqdYxeZPAH+2hITExUZcuXerKey/etJdb3lhKtehwptzYkOo1bFRX4561a9fSpEkTt2OYIJTX74aIJKlqYl7L255EEVi98EsqTrqEljGHmHJnZysQxpiQYUXiLC3/Zgb1Zt5EtDeHl4a0p3rFSLcjGWNMkbEicRYWf/Y6zebewp6wakTdPpMq5xb8Rh7GGFMaWJE4Q/Onv07i4nvZFNGIynfNpfI557kdyRhjipwViULy+5WnvvyZexZWYE7sAOrcO8tOczXGhCy7yqsQMo4c5qv/jeK17d0Z0OF8Lul3DeFeq7PGmNBl33AFtG/HVrY+cwlX7HuL5zsc5F9/aGEFwphT8Hq9tG7dmmbNmtGqVSuefvpp/H7/KdfZvHlzsV7cNm7cON58881TLrN06VLuueeefOefmPF0y5d2tidRAOsWz6bq58NJ0HRWXPgSV/S8we1IxgS98uXLs2LFCgB2797N4MGDOXToEI899li+6xz9Ah48eHCB3ycnJ4ewsIJ9ld15552nXSYxMZHExDwvGcgz4+mWL+3sT+FTUFUWTBtH/c+uI0PKs3Pgp7S1AmFKo9f7nPxY/FpgXtaRvOcvd0ZsTdt38rxCqlatGuPHj+ell15CVfH5fPz1r3+lffv2tGzZkldffRWAUaNGsWDBAlq3bs2zzz6b73Lz58+nS5cu9O3bl6ZNmzJ//nwuvvhi+vXrR7169Rg1ahRvv/02HTp0oEWLFvzyyy8APProozz11FMAdOvWjQceeIAOHTrQqFGjY0OEzJ8/nyuvvBKAr7/++tjItG3atCE1NfWkjLmXP3z4MLfccgstWrSgZcuWfPjhhyf9W9SpU4cHH3yQ1q1bk5iYyLJly+jZsyf169dn3Lhxx5YbO3bssc89evToY+39+/enXbt2NGvWjPHjxx9rj46O5uGHH6ZVq1Z07NiRXbt2FXo75cWKRD5SM7K5970V/G1xGEsrdCXm7m+p16yD27GMKbXq1auHz+dj9+7dTJgwgdjYWJYsWcKSJUt47bXX+PXXX3nyySfp0qULK1as4N577813OYBly5bx/PPPs379egBWrlzJuHHjWLt2LZMnT2b9+vUsXryY4cOH8+KLL+aZKScnh8WLF/Pcc8/luYfz1FNP8fLLL7NixQoWLFhA+fLlT8qY2xNPPEFsbCyrV69m1apVXHrppXm+b+3atVmxYgVdunRh6NChTJ06lUWLFh0rBrNmzWLDhg0sXryYFStWkJSUxDfffAPAxIkTSUpKYunSpbzwwgvs27cPgLS0NDp27MjKlSvp2rUrr7322hlspZPZ4aY8bEiax/IvJjAjbRB/vuxiOlwyDK8N9W1Ks1s+y39eRNSp51eocur5Z2DWrFmsWrWKqVOnApCSksKGDRuIiIgo8HIdOnSgbt3fr01q3749NWoE7nJcv359Lr/8cgBatGjBvHnz8swxYMAAANq1a3fSTY0AOnfuzH333ceQIUMYMGAACQkJp/xcc+bMYcqUKcemK1WqlOdyffv2PZbt8OHDxMTEEBMTQ7ly5Th48CCzZs1i1qxZtGnTBgjsoWzYsIGuXbvywgsvMG3aNAB+++03NmzYQJUqVYiIiDi2R9OuXTtmz559yqwFZUUil5zsLJa+/Xfa/TqeaKnMhzf9jTZNGrody5iQsGnTJrxeL9WqVUNVefHFF+nZs+dxy8yfP/+46VMtV6FChePaypUrd+y5x+M5Nu3xeMjJyckz09FlvF5vnsuMGjWKPn368Pnnn9O5c2e+/PLLgn3Y08id7cTcOTk5qCoPPvggd9xxx3HrzZ8/nzlz5rBw4UKioqLo1q0bGRkZAISHhyMip/w8Z8IONzm2rE3i1zEX0nHzK6yseDFR9yykTZPGbscyJiTs2bOHO++8k7vuugsRoWfPnrzyyitkZ2cDsH79etLS0oiJiTnu/hL5LVdSfvnlF1q0aMEDDzxA+/btWbdu3UkZc+vRowcvv/zysekDBw6c0fv27NmTiRMncvjwYSBwp73du3eTkpJCpUqViIqKYt26dSxatOiMXr8wbE8C+GDxJi78bCAVJZOkDs+S2PtWtyMZU+qlp6fTunVrsrOzCQsL48Ybb+S+++4DYPjw4WzevJm2bduiqsTHx/Pxxx/TsmVLvF4vrVq1YujQoYwcOTLP5UrKc889x7x58/B4PDRr1oxevXrh8XiOy3j0kBDAI488wogRI2jevDler5fRo0cfO6RVGJdffjlr166lU6dOQKBT+q233uKKK65g3LhxNGnShMaNG9OxY8ci+6z5saHCgaWb9zNv9ifcctWlVD2n1ulXMCbI2VDhJj+FHSrc9iSAxDqVSbztZrdjGGNM0LE+CWOMMfmyImFMiAqlQ8mmaJzJ74QVCWNCUGRkJPv27bNCYY5RVfbt20dkZOFujGZ9EsaEoISEBJKTk9mzZ4/bUUwQiYyMPO0FgSeyImFMCAoPDz/uamRjzpQdbjLGGJMvKxLGGGPyZUXCGGNMvkLqimsR2QNsOcPVqwJ7izBOaWCfuWywz1w2nM1nPk9V4/OaEVJF4myIyNL8LksPVfaZywb7zGVDcX1mO9xkjDEmX1YkjDHG5MuKxO/Gn36RkGOfuWywz1w2FMtntj4JY4wx+bI9CWOMMfmyImGMMSZfViQAEblCRH4WkY0iMsrtPMVBRGqJyDwRWSMiP4nISKe9sojMFpENzs9KbmctSiLiFZHlIvKpM11XRH5wtvV7IhLhdsaiJCJxIjJVRNaJyFoR6VQGtvG9zu/0jyLyrohEhtp2FpGJIrJbRH7M1ZbndpWAF5zPvkpE2p7Ne5f5IiEiXuBloBfQFLheRJq6m6pY5AD3q2pToCMwwvmco4C5qtoQmOtMh5KRwNpc02OAZ1W1AXAAGOZKquLzPDBTVc8HWhH47CG7jUWkJnAPkKiqzQEvMIjQ285vAFec0Jbfdu0FNHQetwOvnM0bl/kiAXQANqrqJlXNAqYA/VzOVORUdYeqLnOepxL48qhJ4LNOchabBPR3J2HRE5EEoA/wP2dagEuBqc4iofZ5Y4GuwAQAVc1S1YOE8DZ2hAHlRSQMiAJ2EGLbWVW/Afaf0Jzfdu0HvKkBi4A4Ealxpu9tRSLwRflbrulkpy1kiUgdoA3wA1BdVXc4s3YC1V2KVRyeA/4P8DvTVYCDqprjTIfatq4L7AFedw6x/U9EKhDC21hVtwFPAVsJFIcUIInQ3s5H5bddi/Q7zYpEGSMi0cCHwJ9V9VDueRo4HzokzokWkSuB3aqa5HaWEhQGtAVeUdU2QBonHFoKpW0M4ByH70egQJ4LVODkwzIhrzi3qxUJ2AbUyjWd4LSFHBEJJ1Ag3lbVj5zmXUd3RZ2fu93KV8Q6A31FZDOBQ4iXEjheH+ccloDQ29bJQLKq/uBMTyVQNEJ1GwNcBvyqqntUNRv4iMC2D+XtfFR+27VIv9OsSMASoKFzNkQEgU6vGS5nKnLO8fgJwFpVfSbXrBnAzc7zm4HpJZ2tOKjqg6qaoKp1CGzTr1R1CDAPuMZZLGQ+L4Cq7gR+E5HGTlN3YA0huo0dW4GOIhLl/I4f/cwhu51zyW+7zgBucs5y6gik5DosVWh2xTUgIr0JHL/2AhNV9Z8uRypyInIRsABYze/H6B8i0C/xPlCbwDDrA1X1xA6yUk1EugF/UdUrRaQegT2LysBy4AZVzXQzX1ESkdYEOuojgE3ALQT+GAzZbSwijwHXETiDbzkwnMAx+JDZziLyLtCNwHDgu4DRwMfksV2dYvkSgcNuR4BbVHXpGb+3FQljjDH5scNNxhhj8mVFwhhjTL6sSBhjjMmXFQljjDH5siJhjDEmX1YkjDkDIlJFRFY4j50iss15flhE/ut2PmOKip0Ca8xZEpFHgcOq+pTbWYwparYnYUwREpFuue5d8aiITBKRBSKyRUQGiMh/RGS1iMx0hklBRNqJyNcikiQiX57NiJ3GFDUrEsYUr/oExo3qC7wFzFPVFkA60McpFC8C16hqO2AiEHJX/JvSK+z0ixhjzsIXqpotIqsJDPsy02lfDdQBGgPNgdmB0RTwEhjy2pigYEXCmOKVCaCqfhHJ1t87Af0E/v8J8JOqdnIroDGnYoebjHHXz0C8iHSCwHDuItLM5UzGHGNFwhgXObfMvQYYIyIrgRXAhe6mMuZ3dgqsMcaYfNmehDHGmHxZkTDGGJMvKxLGGGPyZUXCGGNMvqxIGGOMyZcVCWOMMfmyImGMMSZf/w98n9AK8WnPCAAAAABJRU5ErkJggg==\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "simulate_n(1000)\n", + "plt.plot(times, model.mean(params, times), '--', label=\"Deterministic mean\")\n", + "plt.legend()\n", + "plt.xlabel('Time')\n", + "plt.ylabel(r'Population ($C(t)$)')\n", + "plt.show()" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.8.5" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/pints/tests/test_toy_logistic_model.py b/pints/tests/test_toy_logistic_model.py index 11ae30630..0648048d5 100755 --- a/pints/tests/test_toy_logistic_model.py +++ b/pints/tests/test_toy_logistic_model.py @@ -12,7 +12,7 @@ import pints.toy -class TestLogistic(unittest.TestCase): +class TestLogisticModel(unittest.TestCase): """ Tests if the logistic (toy) model works. """ diff --git a/pints/tests/test_toy_stochastic_degradation_model.py b/pints/tests/test_toy_stochastic_degradation_model.py index a7262c90b..d5d34a5ba 100755 --- a/pints/tests/test_toy_stochastic_degradation_model.py +++ b/pints/tests/test_toy_stochastic_degradation_model.py @@ -13,7 +13,7 @@ from pints.toy import StochasticDegradationModel -class TestStochasticDegradation(unittest.TestCase): +class TestStochasticDegradationModel(unittest.TestCase): """ Tests if the stochastic degradation (toy) model works. """ diff --git a/pints/tests/test_toy_stochastic_logistic_model.py b/pints/tests/test_toy_stochastic_logistic_model.py new file mode 100755 index 000000000..67fc04025 --- /dev/null +++ b/pints/tests/test_toy_stochastic_logistic_model.py @@ -0,0 +1,116 @@ +#!/usr/bin/env python3 +# +# Tests if the stochastic logistic growth (toy) model works. +# +# This file is part of PINTS (https://github.com/pints-team/pints/) which is +# released under the BSD 3-clause license. See accompanying LICENSE.md for +# copyright notice and full license details. +# +import unittest +import numpy as np +import pints +import pints.toy + + +class TestStochasticLogisticModel(unittest.TestCase): + """ + Tests if the stochastic logistic growth (toy) model works. + """ + def test_start_with_zero(self): + # Test the special case where the initial population count is zero + + # Set seed for random generator + np.random.seed(1) + + model = pints.toy.StochasticLogisticModel(0) + times = [0, 1, 2, 100, 1000] + parameters = [0.1, 50] + values = model.simulate(parameters, times) + self.assertEqual(len(values), len(times)) + self.assertTrue(np.all(values == np.zeros(5))) + + def test_start_with_one(self): + # Run a small simulation and check it runs properly + + # Set seed for random generator + np.random.seed(1) + + model = pints.toy.StochasticLogisticModel(1) + times = [0, 1, 2, 100, 1000] + parameters = [0.1, 50] + values = model.simulate(parameters, times) + self.assertEqual(len(values), len(times)) + self.assertEqual(values[0], 1) + self.assertEqual(values[-1], 50) + self.assertTrue(np.all(values[1:] >= values[:-1])) + + def test_suggested(self): + # Check suggested values + model = pints.toy.StochasticLogisticModel(1) + times = model.suggested_times() + parameters = model.suggested_parameters() + self.assertTrue(len(times) == 101) + self.assertTrue(np.all(parameters > 0)) + + def test_simulate(self): + # Check each step in the simulation process + np.random.seed(1) + model = pints.toy.StochasticLogisticModel(1) + times = np.linspace(0, 100, 101) + params = [0.1, 50] + time, raw_values = model._simulate_raw([0.1, 50]) + values = model._interpolate_values(time, raw_values, times, params) + self.assertTrue(len(time), len(raw_values)) + + # Test output of Gillespie algorithm + self.assertTrue(np.all(raw_values == np.array(range(1, 51)))) + + # Check simulate function returns expected values + self.assertTrue(np.all(values[np.where(times < time[1])] == 1)) + + # Check interpolation function works as expected + temp_time = np.array([np.random.uniform(time[0], time[1])]) + self.assertTrue(model._interpolate_values(time, raw_values, temp_time, + params)[0] == 1) + temp_time = np.array([np.random.uniform(time[1], time[2])]) + self.assertTrue(model._interpolate_values(time, raw_values, temp_time, + params)[0] == 2) + + # Check parameters, times cannot be negative + parameters_0 = [-0.1, 50] + self.assertRaises(ValueError, model.simulate, parameters_0, times) + self.assertRaises(ValueError, model.mean, parameters_0, times) + + parameters_1 = [0.1, -50] + self.assertRaises(ValueError, model.simulate, parameters_1, times) + self.assertRaises(ValueError, model.mean, parameters_1, times) + + times_2 = np.linspace(-10, 10, 21) + parameters_2 = [0.1, 50] + self.assertRaises(ValueError, model.simulate, parameters_2, times_2) + self.assertRaises(ValueError, model.mean, parameters_2, times_2) + + # Check this model takes 2 parameters + parameters_3 = [0.1] + self.assertRaises(ValueError, model.simulate, parameters_3, times) + self.assertRaises(ValueError, model.mean, parameters_3, times) + + # Check initial value cannot be negative + self.assertRaises(ValueError, pints.toy.StochasticLogisticModel, -1) + + def test_mean_variance(self): + # Check the mean is what we expected + model = pints.toy.StochasticLogisticModel(1) + v_mean = model.mean([1, 10], [5, 10]) + self.assertEqual(v_mean[0], 10 / (1 + 9 * np.exp(-5))) + self.assertEqual(v_mean[1], 10 / (1 + 9 * np.exp(-10))) + + # Check model variance is not implemented + times = np.linspace(0, 100, 101) + parameters = [0.1, 50] + self.assertRaises(NotImplementedError, model.variance, + parameters, times) + + +if __name__ == '__main__': + unittest.main() diff --git a/pints/toy/__init__.py b/pints/toy/__init__.py index 30b2db076..a6aee7fcf 100644 --- a/pints/toy/__init__.py +++ b/pints/toy/__init__.py @@ -36,3 +36,4 @@ from ._sir_model import SIRModel from ._twisted_gaussian_banana import TwistedGaussianLogPDF from ._stochastic_degradation_model import StochasticDegradationModel +from ._stochastic_logistic_model import StochasticLogisticModel diff --git a/pints/toy/_hes1_michaelis_menten.py b/pints/toy/_hes1_michaelis_menten.py index 719a6cf89..fc7b6643b 100644 --- a/pints/toy/_hes1_michaelis_menten.py +++ b/pints/toy/_hes1_michaelis_menten.py @@ -34,13 +34,6 @@ class Hes1Model(pints.ForwardModel, ToyModel): Extends :class:`pints.ForwardModel`, :class:`pints.toy.ToyModel`. - References - ---------- - .. [1] Silk, D., el al. 2011. Designing attractive models via automated - identification of chaotic and oscillatory dynamical regimes. Nature - communications, 2, p.489. - https://doi.org/10.1038/ncomms1496 - Parameters ---------- y0 : float @@ -48,6 +41,13 @@ class Hes1Model(pints.ForwardModel, ToyModel): implicit_parameters The implicit parameter of the model that is not inferred, given as a vector ``[p1_0, p2_0, k_deg]`` with ``p1_0, p2_0, k_deg >= 0``. + + References + ---------- + .. [1] Silk, D., el al. 2011. Designing attractive models via automated + identification of chaotic and oscillatory dynamical regimes. Nature + communications, 2, p.489. + https://doi.org/10.1038/ncomms1496 """ def __init__(self, y0=None, implicit_parameters=None): if y0 is None: diff --git a/pints/toy/_stochastic_logistic_model.py b/pints/toy/_stochastic_logistic_model.py new file mode 100644 index 000000000..bf276e3f2 --- /dev/null +++ b/pints/toy/_stochastic_logistic_model.py @@ -0,0 +1,163 @@ +# +# Stochastic logistic model. +# +# This file is part of PINTS (https://github.com/pints-team/pints/) which is +# released under the BSD 3-clause license. See accompanying LICENSE.md for +# copyright notice and full license details. +# +from __future__ import absolute_import, division +from __future__ import print_function, unicode_literals +import numpy as np +from scipy.interpolate import interp1d +import pints + +from . import ToyModel + + +class StochasticLogisticModel(pints.ForwardModel, ToyModel): + r""" + This model describes the growth of a population of individuals, where the + birth rate per capita, initially :math:`b_0`, decreases to :math:`0` as the + population size, :math:`\mathcal{C}(t)`, starting from an initial + population size, :math:`n_0`, approaches a carrying capacity, :math:`k`. + This process follows a rate according to [1]_ + + .. math:: + A \xrightarrow{b_0(1-\frac{\mathcal{C}(t)}{k})} 2A. + + The model is simulated using the Gillespie stochastic simulation algorithm + [2]_, [3]_. + + *Extends:* :class:`pints.ForwardModel`, :class:`pints.toy.ToyModel`. + + Parameters + ---------- + initial_molecule_count : float + Sets the initial population size :math:`n_0`. + + References + ---------- + .. [1] Simpson, M. et al. 2019. Process noise distinguishes between + indistinguishable population dynamics. bioRxiv. + https://doi.org/10.1101/533182 + .. [2] Gillespie, D. 1976. A General Method for Numerically Simulating the + Stochastic Time Evolution of Coupled Chemical Reactions. + Journal of Computational Physics. 22 (4): 403-434. + https://doi.org/10.1016/0021-9991(76)90041-3 + .. [3] Erban R. et al. 2007. A practical guide to stochastic simulations + of reaction-diffusion processes. arXiv. + https://arxiv.org/abs/0704.1908v2 + """ + + def __init__(self, initial_molecule_count=50): + super(StochasticLogisticModel, self).__init__() + self._n0 = float(initial_molecule_count) + if self._n0 < 0: + raise ValueError('Initial molecule count cannot be negative.') + + def n_parameters(self): + """ See :meth:`pints.ForwardModel.n_parameters()`. """ + return 2 + + def _simulate_raw(self, parameters): + """ + Returns tuple (raw times, population sizes) when reactions occur. + """ + parameters = np.asarray(parameters) + if len(parameters) != self.n_parameters(): + raise ValueError('This model should have only 2 parameters.') + b = parameters[0] + k = parameters[1] + if b <= 0: + raise ValueError('Rate constant must be positive.') + + # Initial time and count + t = 0 + a = self._n0 + + # Run stochastic logistic birth-only algorithm, calculating time until + # next reaction and increasing population count by 1 at that time + mol_count = [a] + time = [t] + while a < k: + r = np.random.uniform(0, 1) + t += np.log(1 / r) / (a * b * (1 - a / k)) + a = a + 1 + time.append(t) + mol_count.append(a) + return time, mol_count + + def _interpolate_values(self, time, pop_size, output_times, parameters): + """ + Takes raw times and population size values as inputs and outputs + interpolated values at output_times. + """ + # Interpolate as step function, increasing pop_size by 1 at each + # event time point + interp_func = interp1d(time, pop_size, kind='previous') + + # Compute population size values at given time points using f1 + # at any time beyond the last event, pop_size = k + values = interp_func(output_times[np.where(output_times <= time[-1])]) + zero_vector = np.full( + len(output_times[np.where(output_times > time[-1])]), + parameters[1]) + values = np.concatenate((values, zero_vector)) + return values + + def simulate(self, parameters, times): + """ See :meth:`pints.ForwardModel.simulate()`. """ + times = np.asarray(times) + if np.any(times < 0): + raise ValueError('Negative times are not allowed.') + if self._n0 == 0: + return np.zeros(times.shape) + + # run Gillespie + time, pop_size = self._simulate_raw(parameters) + + # interpolate + values = self._interpolate_values(time, pop_size, times, parameters) + return values + + def mean(self, parameters, times): + r""" + Computes the deterministic mean of infinitely many stochastic + simulations with times :math:`t` and parameters (:math:`b`, :math:`k`), + which follows: + :math:`\frac{kC(0)}{C(0) + (k - C(0)) \exp(-bt)}`. + + Returns an array with the same length as `times`. + """ + parameters = np.asarray(parameters) + if len(parameters) != self.n_parameters(): + raise ValueError('This model should have only 2 parameters.') + + b = parameters[0] + if b <= 0: + raise ValueError('Rate constant must be positive.') + + k = parameters[1] + if k <= 0: + raise ValueError("Carrying capacity must be positive") + + times = np.asarray(times) + if np.any(times < 0): + raise ValueError('Negative times are not allowed.') + c0 = self._n0 + return (c0 * k) / (c0 + np.exp(-b * times) * (k - c0)) + + def variance(self, parameters, times): + r""" + Returns the deterministic variance of infinitely many stochastic + simulations. + """ + raise NotImplementedError + + def suggested_parameters(self): + """ See :meth:`pints.toy.ToyModel.suggested_parameters()`. """ + return np.array([0.1, 500]) + + def suggested_times(self): + """ See :meth:`pints.toy.ToyModel.suggested_times()`.""" + return np.linspace(0, 100, 101)