diff --git a/README.md b/README.md index 4dedbd87..19a72a7d 100644 --- a/README.md +++ b/README.md @@ -81,6 +81,7 @@ Example [IPython Notebooks](examples/) on how to use the library can be found in - [Example 0.0](https://mybinder.org/v2/gh/neurolib-dev/neurolib/master?filepath=examples%2Fexample-0-aln-minimal.ipynb) - Basic use of the `aln` model - [Example 0.3](https://mybinder.org/v2/gh/neurolib-dev/neurolib/master?filepath=examples%2Fexample-0.3-fhn-minimal.ipynb) - Fitz-Hugh Nagumo model `fhn` on a brain network +- [Example 0.6](https://mybinder.org/v2/gh/neurolib-dev/neurolib/master?filepath=examples%2Fexample-0.6-custom-model.ipynb) - Minimal example of how to implement your own model in `neurolib` - [Example 1.2](https://mybinder.org/v2/gh/neurolib-dev/neurolib/master?filepath=examples%2Fexample-1.2-brain-network-exploration.ipynb) - Parameter exploration of a brain network and fitting to BOLD data - [Example 2.0](https://mybinder.org/v2/gh/neurolib-dev/neurolib/master?filepath=examples%2Fexample-2-evolutionary-optimization-minimal.ipynb) - A simple example of the evolutionary optimization framework diff --git a/examples/example-0.6-custom-model.ipynb b/examples/example-0.6-custom-model.ipynb new file mode 100644 index 00000000..48ea17a5 --- /dev/null +++ b/examples/example-0.6-custom-model.ipynb @@ -0,0 +1,254 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Minimal model implementation" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "This notebook demonstrates how to implement your own model in `neurolib`. There are two main parts of each model: its class that inherits from the `Model` base class and its `timeIntegration()` function. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# change to the root directory of the project\n", + "import os\n", + "if os.getcwd().split(\"/\")[-2] == \"neurolib\":\n", + " os.chdir('..')\n", + "\n", + "%load_ext autoreload\n", + "%autoreload 2" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Model equations" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In this example we will implement a linear model with the following equation:\n", + "\n", + "$\\frac{d}{dt} x_i(t) = - \\frac{x_i(t)}{\\tau} + \\sum_{j=0}^{N} K G_{ij} x_j(t)$.\n", + "\n", + "Here, we simulate $N$ nodes that are coupled in a network. $x_i$ are the elements of an $N$-dimensional state vector, $\\tau$ is the decay time constant, $G$ is the adjacency matrix and $K$ is the global coupling strength." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Implementation" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We first create a class for the model called `LinearModel` which inherits lots of functionality from the `Model` base class. We define `state_vars` and `default_output` so that `neurolib` knows how to handle the variables of the system. Next, we define `init_vars` in order to use the autochunk integration scheme, so we can save a lot of RAM when we run very long simulations. \n", + "```\n", + "class LinearModel(Model):\n", + " state_vars = [\"x\"]\n", + " default_output = \"x\"\n", + " init_vars = [\"x_init\"]\n", + "```\n", + "Next we define a simple parameter dictionary called `params`. In here, we can define all the necessary parameters of the model and change their values later. In this example, we set the timescale $\\tau$, the coupling strength $K$, the integration time step `dt` (in ms) and the duration to 100 ms.\n", + "```\n", + "params = dict(tau=10, K=1e-2, dt=1e-1, duration=100)\n", + "```\n", + "We are now ready to set up the constructor of our model! This method is supposed to set up the model and prepare it for integration. All the magic happens in the background! We pass the `timeIntegration` function and the parameter dictionary `self.params` to the base class using `super().__init__()`.\n", + "```\n", + "def __init__(self, Cmat=np.zeros((1,1))):\n", + " self.params['Cmat'] = Cmat\n", + " super().__init__(timeIntegration, self.params)\n", + "```\n", + "That wasn't too bad, was it? We are finally ready to define the time integration function that crunches the numbers. Remember to put this function outside of the class definition, so we can use use `numba` acceleration to greatly increase the performance of our code. \n", + "\n", + "```\n", + "def timeIntegration(p):\n", + " N = p['Cmat'].shape[0]\n", + " t = np.arange(1, p['duration']/p['dt']) # holds time steps\n", + " x = np.ndarray((N, len(t)+1)) # holds variable x\n", + "```\n", + "Here we prepare the numpy arrays that will hold the simulation results. We have to prepare them before we can execute the numba code. \n", + "```\n", + " # either use predefined initial conditions or random ones\n", + " x[:, :1] = p.get('x_init') if p.get('x_init') is not None else rand((N, 1))\n", + "```\n", + "Next, we make use of a neurolib convention to prepare the initial conditions of our model. If you remember, we defined `init_vars` above in order to use the autochunk feature. The autochunk feature will automatically fill this parameter with the last state of the last simulated chunk so the model integration can be continued without having to remember the entire output and state variables of the model indefinitely. In this line, we check whether x_init is set or not (which it will be, when we use chunkwise integration). If it is not set, we simply use random initial conditions using `rand((N, 1))`. Remember that the convention for array dimensions is `array[space, time]`, meaning that we only fill in the first time step with the initial condition. \n", + "```\n", + "return njit_integrate(x, t, p['tau'], p['K'], N, p['Cmat'], p['dt'])\n", + "```\n", + "We're ready to call our accelerated integration part 🚀!" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Numba time integration " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "```\n", + "@numba.njit\n", + "def njit_integrate(x, t, tau, K, N, Cmat, dt):\n", + "```\n", + "We first have to let `numba` know which part of the code to precompile. We do this by simply placing the decorator `@numba.njit` in the line above the integration function. Easy way of getting 100x faster code! ❤️ `numba`!\n", + "\n", + "```\n", + " for i in range(1, 1 + len(t)): # loop over time\n", + " inp = Cmat.dot(x[:, i-1]) # input vector\n", + " for n in range(N): # loop over nodes\n", + "```\n", + "Next, we do some simple math. We first loop over all time steps. If you have prepared the array `t` as described above, you can simply loop over its length. In the next line, we calculate the coupling term from the model equation above. However, instead of looping over the sum, we use a little trick here and simply compute the dot product between the coupling matrix `G` and the state vector `x`. This results in a `N`-dimensional vector that carries the amount of input each node receives at each time step. Finally, we loop over all nodes so we can finally add up everything.\n", + "```\n", + " x[n, i] = x[n, i-1] + (- x[n, i-1] / tau + K * inp[n]) * dt # model equations\n", + "```\n", + "In this line, we integrate the equation that we have shown above. This integration scheme is called Euler integration and is the most simple way of solving an ODE. The idea is easy and is best expressed as `x_next = x_before + f(x) * dt` where `f(x)` is simply the time derivative $\\frac{d}{dt} x_i(t)$ shown above.\n", + "```\n", + " return t, x\n", + "```\n", + "We're done! The only thing left to do is to return the data so that neurolib can take over from here on. The outputs of this simulation will be available in the `model.outputs` attribute. You can see an example time series below. " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Code" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "import numba\n", + "import numpy as np\n", + "from numpy.random import random as rand\n", + "\n", + "from neurolib.models.model import Model\n", + "\n", + "class LinearModel(Model):\n", + " state_vars = [\"x\"]\n", + " default_output = \"x\"\n", + " init_vars = [\"x_init\"]\n", + " params = dict(tau=10, K=1e-2, dt=1e-1, duration=100)\n", + " def __init__(self, Cmat=np.zeros((1,1))):\n", + " self.params['Cmat'] = Cmat\n", + " super().__init__(timeIntegration, self.params)\n", + " \n", + "def timeIntegration(p):\n", + " N = p['Cmat'].shape[0]\n", + " t = np.arange(1, p['duration']/p['dt']) # holds time steps\n", + " x = np.ndarray((N, len(t)+1)) # holds variable x\n", + " # either use predefined initial conditions or random ones\n", + " x[:, :1] = p.get('x_init') if p.get('x_init') is not None else rand((N, 1))\n", + " return njit_integrate(x, t, p['tau'], p['K'], N, p['Cmat'], p['dt'])\n", + "\n", + "@numba.njit\n", + "def njit_integrate(x, t, tau, K, N, Cmat, dt):\n", + " for i in range(1, 1 + len(t)): # loop over time\n", + " inp = Cmat.dot(x[:, i-1]) # input vector\n", + " for n in range(N): # loop over nodes\n", + " x[n, i] = x[n, i-1] + (- x[n, i-1] / tau + K * inp[n]) * dt # model equations\n", + " return t, x" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [], + "source": [ + "Cmat = rand((12, 12)) # use a random connectivity matrix\n", + "model = LinearModel(Cmat)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "model.run()" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "Text(0, 0.5, 'Activity $x$')" + ] + }, + "execution_count": 5, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAEGCAYAAABhMDI9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAABbvklEQVR4nO3dd3gVVf7H8fe5vaX33ggtdCIiYAMpil1EseDay+rqurq6q79dy7quZe2911WsiEhRBER6kRJ6SYH0ntxez++PG5UoIihJgJzX89wnuZNzZ74z190PM2fmHCGlRFEURVG+p+nqAhRFUZTDiwoGRVEUpR0VDIqiKEo7KhgURVGUdlQwKIqiKO3ourqA3ys+Pl5mZ2d3dRmKoihHlDVr1tRLKRP29bcjPhiys7NZvXp1V5ehKIpyRBFClP3S39SlJEVRFKUdFQyKoihKOyoYFEVRlHZUMCiKoijtqGBQFEVR2lHBoCiKorSjgkFRFEVpp9sGQ0vLd+zc9Qhq2HFFUZT2um0wtNo3UVb2Am737q4uRVEU5bDSbYMhNmYkAI1NS7q4EkVRlMNLtw0GiyUHozGZpqZlXV2KoijKYaXbBoMQgtiYETQ2LkXKUFeXoyiKctjotsEAEBM7kkCgGbtjc1eXoiiKctjo1sEQGzMCgKbGpV1ciaIoyuGjWweD0ZiI1ZpPY5MKBkVRlO9162AAiIkZQXPzKoJBb1eXoiiKcljo9sEQGzuSUMhDS+t3XV2KoijKYaHbB0NM9DCE0Kp+BkVRlDadGgxCiAlCiG1CiJ1CiDv38fdMIcQCIcRaIcQGIcRpHV2TThdBZORA1c+gKIrSptOCQQihBZ4FTgX6AlOEEH1/0uxu4AMp5WDgQuC5zqgtNmYkra0b8PtbO2NziqIoh7XOPGMYBuyUUhZLKX3A+8BZP2kjgci236OAys4oLCZmBBCiuXl5Z2xOURTlsNaZwZAG7NnrfXnbsr3dA1wihCgHZgE37WtFQohrhBCrhRCr6+rqfndhUVGD0GqtNDQs+t3rUhRFOdIdbp3PU4A3pJTpwGnA20KIn9UopXxJSlkopSxMSEj43RvVaAzExoygoeEbNQy3oijdXmcGQwWQsdf79LZle7sS+ABASrkMMAHxnVFcXNyJeLyVOF07O2NziqIoh63ODIZVQL4QIkcIYSDcuTzjJ212A2MAhBB9CAfD779WdADi4k4EoKHhm87YnKIoymGr04JBShkAbgTmAlsI3320SQhxnxDizLZmfwGuFkKsB94D/iA76dqOyZSK1ZqvgkFRlG5P15kbk1LOItypvPeyf+z1+2ZgZGfWtLe4uBPZs+ctAgEnOp21q8pQFEXpUodb53OXios9ESl9avIeRVG6NRUMe4mOLmy7bXVhV5eiKIrSZVQw7EXdtqooiqKC4WfUbauKonR3Khh+Qt22qihKd6eC4Sd+vG11YVeXoiiK0iVUMOxDfPwYmptXqdFWFUXpllQw7ENC/ClIGVBnDYqidEsqGPYhMnIgBkM8dfVfdXUpiqIonU4Fwz4IoSE+bjQNDYsIhbxdXY6iKEqnUsHwCxISxhIMOmhqWtnVpSiKonQqFQy/ICZmBBqNmbr6eV1diqIoSqdSwfALtFoTcXHHU18/Tz0FrShKt6KCYT8S4k/B663Gbt/Y1aUoiqJ0GhUM+xEXdzKgUZeTFEXpVlQw7IfBEEt0dCH1KhgURelGVDD8ioT4U3A4tuJ27+7qUhRFUTpFtw0Gz5Yt1D39zK+2S0gYD0Bt7eyOLklRFOWw0G2DwbV6DfXPPou7qGi/7czmdCIjBlCjgkFRlG6i2wZD1Dlno7FaaXz77V9tm5h4KnZ7EW73nk6oTFEUpWt122DQ2mxEnXcurbPn4K+t3W/bxMRTAXU5SVGU7qHbBgNA7MUXQyBA87QP9tvObM4gIqK/CgZFUbqFbh0MhqwsbCeeSNP77xPy+fbbNinxVFrtG3C7yzupOkVRlK7RbYPB6w5QsqGemEsvIdjQgH32/s8GEhNPA6C2Tp01KIpydOu2wbB+3m5mPb+BUO8hGPLyaHzr7f2OiRS+nNRPXU5SFOWo122DodfwFJCwfUUNsZdegmfTJtxr1+33M4mJp9Haul5dTlIU5ajWbYMhKsFMQp6JrcuqiDzjDDSRkTS+9dZ+P5P0w91JszqjREVRlC7RbYPh22+/ZZtnPk31dmoq/cRcMBn7l1/i2/3LQ1+YzZlERg6kumZGJ1aqKIrSubptMGRlZREMBgjaGtm6rIqYSy5FaLU0vvHGfj+XnHQWDscWHI5tnVOooihKJ+u2wZCRkUFMTAzEN7FzTS1ExxF51pk0f/wJgcbGX/xcUtJEhNCqswZFUY5a3TYYhBAMGDCAZncNHr+LXd/VEnfFFUivl6Z33v3FzxkM8cTGjqKmegZShjqxYkVRlM7RbYMBYMCAAQCI+Ca2LK3CmJuLbcwYmt59l5DL9YufS046G4+3kubm1Z1VqqIoSqfp1sEQFxdHWloaPmsdlTuaaalzE3fllQRbWmj++JNf/FxCwilotRaqaz7rxGoVRVE6R7cOBoCBAwdidzcR0DvYsqQSy5DBmAcPpvH115GBwD4/o9VaSIgfR23tLEIhbydXrCiK0rG6fTAUFBSg0WgwpreyeWkVwUCIuKuuxF9ZSet+hslITj6TQKCVhoZvOrFaRVGUjtftg8FqtdKjRw9aZCWuVi8l6+uxnXwyxvwe1L/wIjIY3OfnYmJGotfHUVU9vXMLVhRF6WDdPhggfDnJ7XGhi3OwcVEFQqMh/vrr8e3ahf3LL/f5GY1GR3LymdTXz8fn++XbWxVFUY403TYYpJQUu8L9A7169cJisUBiIxXbmmiqdhIxfjyG3Fzqn3seGdr3bampKecjpV91QiuKclTp1GAQQkwQQmwTQuwUQtz5C20mCyE2CyE2CSH+11G1PFZaw5hV26j0+NDpdAwcOJDalj1IrZ9N31YitFrir78O744d2OfN2+c6bLZeREYMoKryw/2OzKooinIk6bRgEEJogWeBU4G+wBQhRN+ftMkH/gaMlFIWALd0VD2TkmOQSO7bVQnAkCFDCIVCWPIcbF1WRcAXJPK00zBkZ4fPGn7h//hTUs/H4dyG3V7UUaUqiqJ0qs48YxgG7JRSFkspfcD7wFk/aXM18KyUsglASrn/yZh/hyyzkRsyE5le28zSJgcJCQlkZGTQIsvxuPzsXFOL0GqJu+5avFu34liwYJ/rSU46A43GRGXVhx1VqqIoSqfqzGBIA/bs9b68bdneegI9hRBLhBDLhRAT9rUiIcQ1QojVQojVdXV1v7mgGzOTSDPquWtHOYGQZMiQIbTYmzAl+Sj6pgKAqNNPR5+RQf2zz+3zrEGniyAxcQLV1TMIBt2/uRZFUZTDxeHW+awD8oGTgCnAy0KI6J82klK+JKUslFIWJiQk/OaNWbQa7u2Rxhanhzcr6ykoKMBgMKBNbqS2tJXq4haETkf8ddfi2bQJx8KF+1xPaspkgkEHtbVzfnMtiqIoh4vODIYKIGOv9+lty/ZWDsyQUvqllCXAdsJB0SE8/iATE6I4PsbGwyXVtKKhf//+VDSUoDPD+q/DJzhRZ56JPjOTuief2ucdStHRwzCbs9TlJEVRjgqdGQyrgHwhRI4QwgBcCPx07OrphM8WEELEE760VNwRxby9vIwx//0Gly/Iv/LTcQSD/Ke4iiFDhhAIBIjq42bX2jrsjR6EXk/CTTfh3bqV1lk/fxpaCEFqyiSam1fgcpV0RLmKoiidptOCQUoZAG4E5gJbgA+klJuEEPcJIc5sazYXaBBCbAYWALdLKRs6op6C1Egqmt28tKiYXlYTV6Yl8G5VA7W2aFJSUqjzlCCRFC0Mz+8cOfE0jL16Uff0U0i//2frS0mZhBA6Kire64hyFUVROk2n9jFIKWdJKXtKKfOklA+0LfuHlHJG2+9SSnmrlLKvlLK/lPL9jqplSGYMEwek8NKiYmpaPdyWk0ycXsfdOysoPGYYDY31xPeBzYsr8XuDCI2GhFtuxl+2m+ZPPv3Z+ozGRBISxlFZ9RHBoKejylYURelwh1vnc6e6Y3xvAqEQj3+1nUidlrvyUljT6mJrYhpmsxmPtRKvK8DWZVUA2E46CfOgQdQ/+ywhz8//zz897RICgRZqamZ29q4oiqIcMt06GDLjLEw9LpsPVu9ha3UrFyTHMjTSwgOlNeQPGUrpnl3EZOnZsKAcGZIIIUi49c8EamtpevfnD2VHRw/Das2nvOKdLtgbRVGUQ6NbBwPATaN7YDPqeHDWVjRC8EivDFoCQeYnZiOlxJDRTHONi5IN9QBYhw3DOmoUDS+9RNBub7cuIQRpaRdjtxfR2rqhK3ZHURTld+v2wRBtMXDT6Hy+2V7Htzvq6Gszc016Ah81OtD3G0RJ1VYi4g18N7fshwfcEv58C8GWFhpefuVn60tJPhut1kJ5uTprUBTlyNTtgwFg6ogs0mPMPPDFFoIhyW05yaQZ9cxMzqXV5SJ+QJCaklYqtjcDYC4oIOqsM2l84w185e0fxdDpIkhOPpua2pn4/U1dsDeKoii/jwoGwKjT8tcJvdlabefD1XuwarU82DOdEn+InT0HUN68DVOknu/mlv3wmYRbbgGNhrrHHvvZ+tLSLiYU8lJZqR54UxTlyKOCoc0ZA1IYmhXDI3O30eL2My4+ilPjo1iSnM3WxmYyjzGyZ3MjtWWtAOhTUoi74nJaZ83CtXZtu3VF2HoTHX0se8rfIhTa97zRiqIohysVDG2EENx7ZgGNLh+Pf7UdgH/lp6HXaFnWazA13p0YzLp2Zw1xV16JNiGe2v889LMB9jIzrsDrraKuTo2fpCjKkUUFw176pUUxZVgmby8vY1u1nTSTgb/mJlMSncDcZgc5x9rYtbaOpmonABqrlcRbbsG9fj322e2HyoiPH43ZnM3u3a+qSXwURTmiqGD4idvH9cJm1HHPjE1IKbkyLYE+FgNL8gdQbShHq9OwZs6PZw1RZ5+NsU8fah/9LyGv94flQmjIzLiCVvsGWlrWdMWuKIqi/CYqGH4ixmrgtnE9WVbcwKyianQawZN9s/Hojbzig/yRsWxfUU1zjQsAodWSdMcd+CsraXj11XbrSkk5B50umt17XuuKXVEURflNVDDsw0XHZtEnJZIHvtiMyxdgQISFKxIj2ZqUydqYRrQ6Datm/TiKqnX4sURMmEDDiy/hKy//YblWayE9bQp1dV/icpXta1OKoiiHnQMOBiHEk0II0ZHFHC60mnBHdGWLh+cW7ALg7j45JAa8vCiNZJ+YxI6VNT/0NQAk3XkHaLXUPPifdutKT78UIXTsKX+jM3dBURTlNzuYMwY7MEMIYQUQQowXQizpmLK63rCcWM4dnMaLi3axo8aOSavh31kJ2A0m/hdnR2vQsuqL0h/a65OTSbjhehxff419r5nejMYkkpJOp7LyQ3y+xs7fEUVRlIN0wMEgpbwbeA9Y2BYItwJ3dlRhh4O7JvbBatTx90+LCIUkp/fMZaSjnrkYEKMT2LG6hoZKxw/tY6dOxZCbS80D/27XEZ2VdS2hkJs95W92xW4oiqIclIO5lDQGuBpwAvHAn6SU33ZUYYeDOJuRv5/Wh1WlTXywOjzN5/39ehDhdvJ8pAdMWlbNLP2hvTAYSL77Lvx79rTriLZZ80lIGEd5+VsEAvafbkZRFOWwcjCXku4C/k9KeRIwCZgmhBjdIVUdRs4fms6xObH8e9YW6uxe+ubmMLmlkio0FI2PZ9d3tdSUtv7Q3jpixI8d0Xv2/LA8O+t6AoFWyit+Ply3oijK4eRgLiWNllIubvu9CDgV+FdHFXa4EELwwDn98fhD/OuLzQBcMWwIfSpLmWH0UZthZNknO9s9xJZ05x0InY7qf97zw/LIyAHExh7P7t2vqhneFEU5rP3m21WllFXAmENYS6fy7m6l6ZMdyNCvP5XcI9HG9Sfl8dm6Sr7ZXkdubi7neRqJ8HqYcVwEJcXNlG38cWpqfXIyCbf+GefSpbR89tkPy7OzbsDvb6Cy6oMO2SdFUZRD4Xc9xyCldB+qQjpboNqFc2U1ztXVB9T+hpPzyE2wcvf0Ily+IONPOJ4Tt6ymWgT5Zngkyz7dRWivkImZMgXzoEHUPvgfAg3h0IiJGUZUVCFlZS8RCvk6ZL8URVF+r277gJvlmCSMWWZaZpUQbPX+anujTstD5w2gvMnNQ3O20qNHD46xmSis2c3yNB0rQz62Lf8xZIRGQ8q/7ifoclHz7wd/WJ6TfQNebxVVVR93yH4piqL8XgdzV9JNQoiYjiymM4mlTxPTfCMyEKRp+q4DGujumOxYLh+Rw1vLylhe3MiYMWMYuG09mSLEF8fZWDCnmIAv+EN7Y48exF97La1ffIHjm28AiI09gcjIwZSUPkMw+OuBpCiK0tkO5owhCVglhPhACDHhiH8KuvdEdLp6ooyf4NncgHuvPoL9uX18L7LjLPz14/UkpWWSn53FCUUrcBkEH+bpWDd/d7v2cddcjSEvj6p77iXocCKEIC/3z3i91VRWvt8Re6YoivK7HOwDbvnAq8AfgB1CiH8LIfI6qLaOFZcHl3yMTfsJen05zZ/tIOTy/+rHzAYtj5w/8IdLSmPGjMHWUMMk3GzJNPJ6URWOph/PBDQGAyn330+guprahx4CICZmBNHRx1Ja9jzB4BHbTaMoylHqoPoYZPh6S3XbKwDEAB8JIR7ugNo6XspAxEXvEqN9gpDDS/Pn2w/oY3tfUtrtNdG7d2/ily1koMnAzIFmPvt8R7v2liGDib3icpo//BDHokUIIcjN/TM+Xx3lFe90xJ4piqL8ZgfTx3CzEGIN8DCwBOgvpbweGAqc10H1dbzsURgu/CcRuum41jbi2Vp7QB+7fXwvcuKt3PHxBoaPOpGAz8ulrRVodYLHzB5272hq1z7hT3/CmN+DqrvuJtjcTEz0McTGHk9Z2YsEAo5f2IqiKErnO5gzhljgXCnleCnlh1JKP4CUMgSc3iHVdZbepxF51jB0Yg9N764h5Pz1TmGzQcsjkwZQ0eTm2WW1DBw4kJKVy7k3PY6KeB1/X17c7vZVjdFIyn/+Q6Cpier7w88F5ub+Gb+/iT17Xu+wXVMURTlYBxMMJillu0kFhBAPAUgptxzSqrqAOOYiYke5CfrNNL/wIRzAXUqF2bH88eQefLimHJHaDyEEpvXLmag3My9Nw7vflrZrby4oIP7662j94gta58whKnIgCfFjKdv9Ml5ffQftmaIoysE5mGAYu49lpx6qQg4HholXEpFThqsuC/f7Lx7QZ/40Jp+BGdH8c04JAwuPZdOmTdyeZiTZC/e6mqhobt+5HH/NNZj69aP6nnsJ1NXRo8cdhEJeSkqe6ohdUhRFOWi/GgxCiOuFEEVALyHEhr1eJcCGji+xc0VeeQl6ayNN69MJLnz5V9vrtRqeuGAQ/mCI98vMREREsOirL3kmPwO3XnDl0h2E9jr7EHo9qQ/9h5DbTeXf78JsyiItdQqVle/jdO7qyF1TFEU5IAdyxvA/4AxgRtvP719DpZSXdGBtXULotMReeTIhYaXpy2bkul9/1iAn3so/z+jLkpIWtOn9qaysJMJZweVeI+vMIR75rv2zDca8PJLuvAPnt9/S+Mab5OTchEZjZueuI/PmLkVRji6/GgxSyhYpZamUcoqUsmyv11E7HZk+NZKo8Tl4QsNxffwxbP7sVz8zuTCDCQXJPLvBT0x8EvPmzePOE7IpqA3wZHMjyxvaz8MQfeGFRIw9hdrHHye4rYrsrOuor59HU9OKjtotRVGUA3Igl5IWt/20CyFaf/qz40vsGrYTsjBm22gOXIf/g/+D7XP3214IwYPn9ifOauJrRwp2u50VK5fxVN8sopwhrlhXTJ3P3659yv33o4uLo+K2v5AaOxmjMZkdOx8kfKOXoihK1ziQM4ZRbb/eA/SRUkZIKSO//9mh1XWwFvcvP+ksNILYKX0RJjONobuQ718OxQv3u74Yq4FnLhrMplYD7oh0lixZQmIy/KXFSKsMcdXaYgJ73cKqjY4m7dFH8O8pp/6Bh8nL/Qt2exHV1dMP0R4qiqIcvIO5K8kGfCmE+FYIcaMQIqmjiuoMbywpYdzj37C7wfWLbbRRRmIm98LvS6VFcxO8NwXKlu13vYXZsfx1fC8+q4sHjY5Zs2Yx5cx8ztzgYYXLzUPFVe3aWwoLib/hBlo+m4FpBURGDmLnrofUFKCKonSZgxkr6V4pZQHwRyAF+EYIMa/DKutgI7PMjPAv56JXllPV8svjFZn7xGEbmYrDMQq3fiy8ez5UrNnvuq8+PpdRfdJY7k2hpKSE4t3buWl4NoN3eXh6Ty1z61vatY+//joshYXU3Psvsk2X4/M1UFzy5CHZT0VRlIP1W+ZjqCU8VlIDkHhoy+k8+dte4nH5CMNd33DJKyuod/zy085Rp+agT7XS5LqWgDEX3j4Xqjf+YnuNRvDo+QNpsWTSImzMnjOH7MHRXOszkdIU4MZNZZS6f9ye0GpJfey/aGxWWm97lpSE8ygvfwuHY9sh3WdFUZQDcTBjJd0ghFgIfA3EAVdLKQcczMbahuveJoTYKYS4cz/tzhNCSCFE4cGs/2A0DLuKuvgePKx9jrTm1Ux9deUv9jkInYbYKb2RIWg0PorURcDbZ0PdLw+6F20x8OwlQ1nqy8LpdLJgwQLGXtyHKWvcBH1BLi8qwRn4ce4GfWIi6Y8/jm/PHqxvtKDVRrBt+z0HNE+EoijKoXQwZwwZwC1SygIp5T1Sys0HsyEhhBZ4lvDT0n2BKUKIvvtoFwHcDHTofZvTdk3nbIubFmsMr5meQFO7mT+8vhKnN7DP9voECzHn5OMr99KS/RYg4M0zoH7nL25jUEY0159WyNZAIitXrqSxtZazz8rn7MV2tjk83LRld7uH3yyFhSTedhuuL74huXo4zc0rqan5/FDvuqIoyn4dTB/D36SU637HtoYBO6WUxVJKH/A+cNY+2t0PPAR4fse2ftW1A65lZN7pTI7W4tNp+Sjyv9SVF3P1W6vx+IP7/IxlcCLW41JwrHbiGv4BhALwxkSo37HP9gCXjcgmvWAYDmngfx98TM7gOMamxjB2g4tZ9S08Wtp+zunYP1xGxPjxBO9dgFWTw46dD+L3H7V3BSuKchj6Lc8xtP7G5xjSgD17vS9vW7b3toYAGVLKLw5ivb+JVqPlgeMfoHf2KVwaa0IbaOWLuCfYWLx7v+EQPTEXQ2YETV+58Z8xHWQQ3jj9F8NBCMEDkwZTHVWA297Mp7O+4sQpvTihIsiwmiCPldbweW1zu/YpDzyAMSsH27N2fL56du56qAOOgKIoyr4d8HMMez2/ENkRzzEIITTAY8BfDqDtNUKI1UKI1XV1db95m3qNnkdPepSYzOO5IT6KCGcp85OfY83Ocq56czVu38/DQeg0xF7cB2HQ0DDLS2jKjLZwmPiLfQ4mvZZHrhzHbpHEprUr2VNfxZipfRizqIWePg1/2rKbjfYfb5vV2qykP/M0+pIQUWtiqax8Xz0RrShKpzmYzuef/bN1X8v2o4JwP8X30tuWfS8C6AcsFEKUAsOBGfvqgJZSviSlLJRSFiYkJBxECT9n1Bp5avRTuDIKuTMhjrjm9SxIe4VVu6q46q1V+wwHXZSR2Cm9CdS7afpGIKd+Hh6m+83TfzEcUqLMXH3h2bilnjff+4jU3tEMPTmDibMbsEm4rKik3ZPRxtxc0h5/DPM7dvQOC1u2/o1gsEOvrimKogCdO+z2KiBfCJEjhDAAFxIemA/4YUymeClltpQyG1gOnCmlXH0Q2/hNLHoLz415jl1p/bkvIYGk+qXMz3qLlbtquPLNfYeDKS+aqAk5uIvqsW+ywh9mhsPhjYlQt+/bTI/rlUJe4YkY/Hb+8/onHHd2HtkJViYtaqXeF+APRSW4gz8Oh2E7/nhSbr2DyFd9uN1llJQ+3WHHQFEU5XudNuy2lDIA3AjMBbYAH0gpNwkh7hNCnPlbd+BQiTJG8eLYF1md0ov/JiSTVv01X+d9wIriOq54YxUu38/vVrKdkIZlUAKtX5bhrokLhwOEw6F26z63c+0ZowjFZBGo2MzrX69m/FX9SG4I8oddQb5rdXHjlrJ2dyrFTJ1K0oDJmJdqKCt9Cbt9U4fsv6IoyvfEr90nL4SIAmKAB4G9nz2wHw4jrBYWFsrVqw/dSUW1s5rLZl/GedWlXFNfTUn2ZMZsO4vC7DhevayQCJO+XXvpD1H38gb8VU4SrhuIwVAZvqQUCsCln0LKwJ9tw+ly8+/HnsLrD3Lq5MtIbdEw/60tVJ6dwqtGL9dmJHBvjx/75aXPR/F1U9l95mpMMZkce/xsNBrjIdtnRVG6HyHEGinlPp8VU8Nu/0SyNZmXx73MtKR03olNJKf0A77qO4fvyhq56OUVNDp97doLvYa4S/uisehpeGsTQVM2XD4b9BZ44wzY/fNOY6vFzNQpk7EKH29/MB2ZbaHX8GRSPqvifJONF/fU8Wr5j53qwmAg69HniJuTiCtYxo5193X0YVAUpRs7mM7nN4UQ0Xu9jxFCvNYhVXWxzMhMXhv/Gq8nZfBBbAJ5u95mQcFstte0MvnFZVS3tO8E1kYYiLusLyFXgIa3NiMjc8LhYI0PPyG9a8HPttErL4ehx44gW9Rx12uz6XN6FvFpNvp/WMnoCBv/t6OCL/caU0kXG0uvv76LdaWJ8qb3qd/9dUcfBkVRuqmD6XweIKVs/v6NlLIJGHzIKzpMZEVm8eqE13ghKYNpsQlk7HibbwpmUd3iZtILSylrcLZrb0i1EXtBL3x77DR+tB0ZmR4Oh5gc+N9k2DrrZ9uYOG4MsYnJ9PHv4KZ3lnDC5X3QABO+aqS/zcy1m8r4ruXH7RhzcygY9zLaBsGm727EZ//tt+oqiqL8koMJBo0QIub7N0KIWEB36Es6fGRHZbeFQzrvxSaSvO1tvunzOS6Pj0kvLGNbdfuhsc394omckI17fR0tc0ohIincIZ3cH6ZdAkUftWuv1Wq55MLJmPUaUhvX89dZGzlxam/sZQ6u2Rog0aDj4g3FbHP+eIYSMWQEPWNvI2DzsX7aecjAvofwUBRF+a0OJhj+CywTQtwvhPgXsBQ46icpzonK4dUJr/FSUhrvxiURt/UdFvSZgZYQk19cxpqy9l0tESemh4fNWFSOfXEFWGJh6meQNQI+vgrWvNGufWxsLOedew7xGie+srU8samcoadmUb2khvs9VvQawYXrd7HH82PfRuop15HkPoHW3Aq2P3e1GmhPUZRD6mDGSnoLOBeoAaqAawg/hHbUy43K5dXxr/FyQgpvxycTtfldvs6bRoJFw0Uvr2DOxh/HOxJCEH1GHqaCOFq+KMa1oQ6MEXDxh5A/Fj6/GRY9Gn7moU3fvn059thj6aurZeOmTXwedJLVP46dHxbzhC0eVzDEBet2tXsAru/EFzE7E6nIXUTFs6ozWlGUQ+dg52MwAlnA7cC9hJ9H6BZyo3N5bcLrvJaQwsvxyVi3fsSs5JcZmGLi+nfX8ObS0h/aCo0g7sJeGDIjaZy2DW9xM+jNcMG70H8yzL8fZt0OoR8fnBs7dizp6emcbCrj81U72J5jJDbFwq43tvNcWgpVXh8XrS+mtW2obo3GwKAx0xAGPcXmt6l98flOPiKKohytDuQBt55CiH8KIbYCTwO7CT//cLKU8pkOr/Awkhudy5unvsWHyVk8kpiCftdc3jM/yum9bPxzxiYenL2FUNuczkKvJW5qX3SxJurf2oyvygk6A5zzIoy4CVa9DB9dDv5w/4FOp2PSpEmYjXrOjizjuUXbcQ6LQasTVL6+g+d7ZLDF6WbqhmJcbU9HWyyZ9B3wGP4cSXHpYzS+/U6XHRtFUY4eB3LGsBUYDZwupRwlpXwa2PfQo91AZmQmb0x4g4XJefwjOQVRvoynvP/gmsIIXvymmFumrcPb9q96rVVP/BX90Bi01L9ahL/OBRoNjPsXjHsANn8G75wH7mYAoqOjOe+889D67EyKq+T+BdvRHJ+Is8mL+/1SnuqZyYoWJ5cXleBpC4ek5NNIS7kI59gQZZ/dT/PHn3TVoVEU5ShxIMFwLuE+hQVCiJeFEGMA0bFlHd5Sbam8MeENipJ78eekJEK1m/lb1a3cf1IUM9ZXctlrK2l2hTuLdTEm4q/qDxLqXyki0Nh2h9GIG+G8V2HPCnj9NGitBKBHjx6MGzcOk7Oac5JauGfpTswjE6jc0UzE3Coe75XBoiY7l28swRsKh0N+z7uxWfvQfCXseeIuWmfP7pLjoijK0eFAnnyeLqW8EOgNLABuARKFEM8LIcZ1cH2HrQRLAq9PeJ2q1P5clRRPwF7BpZuu4vUJRr4ra+bsZ5ews9YBgD7RQvyV/Qh5Q9S9WkSwtW2+5/6Twp3SzWXw6rgfxlcaPnw4gwYNIqplB6emB7hnQynWwji2r6ghY1kTj/bKYEGjnSs3luINhdBqjQwY8Bwai4Xmm3Xs+ftttM76+XMTiqIoB+Jg7kpySin/J6U8g/CQ2WuBOzqssiNAjCmGV8a/gj99KBcmRuMKejh56WV8cZoHuyfAOc8tYdH28ENohlQb8VcUELL7qXtlI0Fn2x1GeSfDH76AgDccDrsWIITg9NNPJz09nfSWjZycYeDeneVYC6JZ99VuCjY7ebhnOvMaWrl2Uym+UAizOZN+/Z7CF+vB/kcr5bfdRsvnM7vw6CiKcqQ62LuSgPBTz21zIow51AUdaSINkbw09iUSsk/k9Dgz9eZIesy7kq9OKiUt2szlb6zizaWlSCkxZkYS/4e+BBo91L9S9GM4pA6Cq7+GqPRwn8Pq19HpdFxwwQWYzWYKPEWMzLJwf0UV5lwbSz7aybHlAf6dn8ac+lau21SGLxQiLu54euTdhjO3Ad8fUqi84w5aPvusS4+PoihHnt8UDEp7Fr2Fp0c/zfCe5zAxWkNxXCaxX9/GjL7zOblnAv+csYm7p2/EHwxhzI0mfmpf/HUu6l8uIuhoe3AtOhOumAN5o2HmLTD3LiKsFqZMmYLX4+E4uYVjsiL5d2Md+hQz89/cwiktGu7vkcas+hau2FiKJxgiM/MaEhMn0lBYhjy7B5V3/o3mTz7t0uOjKMqRRQXDIaLX6PnXqH9xQf8rOc8WYFlKLwxLH+dl6/PcdHwa767YzaWvrqDe4cXUM4b4ywrw17up2zscTJEw5X0Ydg0sewamXUpKXCSTJ0+mob6e8eZijusRy6OuRojSM/uFIk5z63io7bLSJW23svbt8x9stl7UjS9DN2EAVXfdRdO0D7r2ACmKcsRQwXAIaYSGWwtv5c/H3M41RhcfZhTApk/4S/nNPH9mEmt3N3PG04tZu7sJU34M8X8oINjooe6lDQTtbeGg1cFpj8Cpj8D22fDaBHokWjjzzDMpKy3hrJhKTuqbyNPBFgJmLTOf3cA4t46n+mSytNnBheuLcYQMDOj/EhqtidpJFRjHDaP6n/+k/qWX1fAZiqL8KhUMHeCygst48IT/8G+9m3/n9ifUsINTl0xh9rkGdFrBBS8u590VZRjzooi/vIBgs5e6FzcQbPH+uJJjr4Ep06CxGF48kUGxHkaPHs3GoiImpTQzdkgqz9GKTwefP7OeE906XirIZp3dxfnrduHSJjFwwMv4/E3UX9yE7awJ1D32GLUPP6LCQVGU/VLB0EFOzz2dZ095lpm6AJenZ+DV6sn94kLmnlDGcXlx3PXpRm7/aAMyI4L4K/oRtPuofX49/nr3jyvpOQ6umgdGG7xxOsebd1BYWMjSpUs4L9XBmcMzeFHrwC0kM55ax7EeLa/3z2G7y8M5a3diN/SiX78nsTs203iRi+iLp9D4+utU/f0uNSqroii/SAVDBxqROoJ3TnuHGmsMp8YZaUjshWXOLbye8jE3j87hozXlTHphKbWRehKu7o/0B6l7fj2+CsePK0nsA1fPh9yTEF/cymmhr+hf0Jf587/m7DQ3F5+Uy6s6F45gkM+eWMsAB7w7IJcqr4/Tv9tBvWkEPXv+H/UN82idLIm76Y+0fPop5X+6mZDH88vFK4rSbalg6GB50Xm8e9q7pCYUMMbYTFGPE9CseIE/l/+FtydnUtbg4vSnF7OgJTxntNBrqHtpA55dzT+uxBwDF02DUbei+e4Nzm5+hV552cyePZsJKR5uOasPb5rcNPsDTH98LbkNQaYPySckJWet3cEe67lkZlxJecVbOE6TJP3j/3AsWMDuy68g0HhUztCqKMrvoIKhE8SZ43h1/KuMyz2Ni4KlfFgwDlm1juO/PpevztGQGWvh2rfX8MCyEiKv7o82ykj96xtxb6r/cSUaLZzyTzj/DbS1RZxf/RC5qfHMmDGDYdFOHrhkMP+zemkKBvnsyXVElbuZObQnSQY9F67fxdbI60hJPo+SkidwHu8j7fHH8WzeTOmFU/AWl3TdwVEU5bCjgqGTGLVGHjr+Ia4feD33ubZyV+/hBIw2kqdP5tMBK7hiRBZvLC1l8rurcUzKw5Bio+GdLTiWV7VfUcE5cOVX6PR6Lqy6n4xoHR9//DEZ1PHcNcP4ONpPPUFmPL2e4I5WZgzJZ2CEhWs272ax7S8kJIxj+477cQ50kvXmG4QcDkqnTMG5cmXXHBhFUQ47Khg6kRCCGwbdwMMnPMxX7nLOToyiOe8kdAvu4x+O+3njwnwqmt2c9cpyFh8Th6lXLM3Td9L8RTEytNedRMn94NpvMPQcw8VNj5NpcvLpp59iat3D238cwbwkqBFBvnh+A7Xf1fPBoDwmxEfxj11VvG+8k6joUWzZeieOtFqyP5iGLi6O3VdeRcuMGV13cBRFOWyoYOgCp+acyjunvUNAb2J0qJT1hRfDzq85acG5zJtkoCA1ils/2cB/LH4YloTj2woa3tlCyLfXaOfmGLjwXYzj7+Viz5vkaKuZPn06jortfHDzSDb1MrFbBJn3+ma2zNvDywVZXJuRwGsVTTwi7kLYjqFo459oMW4h+73/YRkyhMq/3kHt408gg912VHVFUVDB0GV6xfbi/YnvMzS5kEsavuW1YZORGi3xH57D+z0XcPPJ2XyytoIpO8vZNSIRz5aG8INwrT/O/YwQcNwfMVzxOVPMi+lBGZ9//jmlm9fzznXH4RsVz1Z9gGWf7uLbaTv4Z24qj/XKYEmLi/8L/h27dRRFG2+kwbeczJdfIvr882l48UX2XH89wdbWrjs4iqJ0KXGkP+xUWFgoV69e3dVl/GaBUICnvnuK1ze9znGx/XjKrce0+TPIGM6GYx/hxlkN7GlycVVBKhdvdWK0Goib2hdDmq39ilyNBD6+jg93GdhGD04eNZzjR4/j+W92sXp6CYVeHWn94zj9mn6sdLi5alMJISm5Xf8mOe6Z9Ct4ksTECTS9P43qBx7AkJpK+rPPYOzRo2sOjKIoHUoIsUZKWbjPv6lgODzMKZnDP5b+A7POzMtJY+m59AUQGjwT/su9Jb14b+UeesdbudujJ8cLMeflYxmU2H4loRDBxU8wY/4K1tOHwl5pnHbBlczbUsvrr29glEOHNc3C5JuHUKuXXLqhmBK3lyuNcznB8yr9Cp4gKfE0XGvWUH7zLUiXi9SHHyLilFO65qAoitJhVDAcIXY17+K2b25jV/Mubs09l8u2LkJUrIEBF7Aw9zZum1lGq9vP9VGRnNsYImpUOlGn5iC07SfUk+VrmPfOYyzx5NMnTnLu1XdQ3BzkXy+sprA6hM6iY/Kfh2BMsfDHzWV81dDKSfqNXOr7D4N730Nq6iT81dWU3/QnPEVFxF17LQk33YjQ6broyCiKcqipYDiCuANuHlzxIJ/u/JRj4gfxlC4D2/IXwJZI69j/ctu6JL7cXMMAm4nbHVp658YSe1FvtDZD+xX5nCx7617mllvIMjRx4dSrCcT34h9vrCVpvR2LEIz+Qx/6HpPMU2U1PFRSTZamnpuC9zAqbypZWdcQ8nqpvv9+Wj76GHPhUNL++1/0SUldc2AURTmkVDAcBkKhIM6mJlrrarE31uNzuwl4PQT8fjRaLTq9AZ3BgCUqGmt0DEtb1/Dgxv9i0pl4svflDF7yItRtQQ66hC9Sb+L/5uzG4QlwGQYui4gg6ZK+GDIifrbdojlv8unyncRg56KTehN7wjW8/NVOds/YTUpQQ97oNMaf35Nvmuxcv6kMb9DNtaHHODdzID3y7kQIQcuMGVTdcy8ao5HUhx/GdvyoLjiCiqIcSioYuoDX5aKsaC0VWzZRvWsHtaXFBHzeX//gXrR6PS22ALUmB7k5vfhDfIjYnR8iIlJpGfcYd21IYOaGKvI1Wu7EzLDTemAbmYoQ7S8tlW3+jvc/+gQZCjA5ZQ+5FzzIshozH764gTyXQJdp4bJbhlKnkVy1sYT1djcT5Wf8KclF/z4PoNHo8BYXU3HzLXh37FCXlhTlKKCCoZMEfD52rFjCpkXz2bOpiFAwgM5gJDEnj+TcHsSmpRMZn0hEXDwGixW9yYROrycUDBLw+fB7vbhamnE2N+JoaqSlpoq6PWWUlW5BtIZDJdlkZ2JmMdHaVhpTTmF5wd+4Z2EzDQ4vF2Lg+l7JpF7QG41F3662xsZG3nvteeodXiZql1J46sXU5l/Io8+sIWOPF79Jy9l/HEBqbjT/3FnBm5UN5Mqd/D1qMeMH3IdeH0nI7abm3/+m+cOPMA8eTOrDD2HIyOiKQ60oyu+kgqGDuVpbWD3zU4rmzcHjdBCVmET+sSPJGzKM1F590Gi1v3sbi3bN5+nZ/8ZY5+W4UA9GuNYx0FaMO6hnqW84b0Vdyrd2G2lo+IvFxsSpAzBmR7Vbh8fj4eP332FHaTnDWMu4HIHmjCd5ZYGTpnlVmIDcCRmccVZPZtU1c8vmXfhCPq7Vf8YtQ2/AbM4EoOXzmVTfdx8EgyTd9Xeizj33Z2cpiqIc3lQwdBCvy8WKT6exdu5Mgj4/+ceOYODYU8no2x+hOfTPDjZ7mrlv+X18VfYVQxOH8H9R40le/BBWVymlrjjeajiOGTFnUq+zcSI67hieSZ8zere7aykUCvHVl1+ybPly0kUN5+u+Jmr8nXwXcTqfvbCRRDf4My1cfcsQmjRw7YaNrHFqOF4s4/H+x5AeF/7vyF9ZSeXf/o5rxQpsp4wh5b770MXGHvJ9VhSlY6hgOMSklGxfvpiFb76Mo7mJPiNPZPh5FxKbmt4p256xawYPrnwQgL8MvoXzWpoQ8+9HBgNst5zEo7sHscjUH63QcJEIcfZIMz1Hj8Bosfywno0bNzLjs+nogm4mhT4jNyMNx9hHeOp/TqJK3dgNcOLUPhwzJJlHdm7l6Qo3cdTzQIafM3tMDNcSCtH45lvUPfYYmshIUv51PxEnn9zhx0BRlN9PBcMh5HbY+fKFp9i5ahmJOXmMveqPJPfo2Wnb/16Fo4J/Lv0nK6pWUJhUyL/6X0/aosdh2yxkYh+2513P3cuiWeXRkS0FZzk20yMvQO/jTyRnUCE6vZ66ujo++OAD6uvqOFn3HaOCS9CM/BOf+yaz5YtKTEGgfzTXXTOI1a11/HHjFqpCMZxn3clDgydi05sB8GzbTuVf/4p32zaizj6bpDvvQBsd3enHRFGUA6eC4RCp3L6FmU88jLO5iVFTpjJ04lloNL+//+C3klLy6c5PeXTVo/hDfm4c9Ecu0cSgnXs3tOyGgnNYmHQDd85toDoYYkTQR7/qmUTqHeQPG0nBiaNJyM3n888/Z+PGjeRFeDnb/iYR0fHUj3yIV2fYsNX5qbcIzry2P/m5kfx13Ww+cWSRKhp4uqAHIxOyAAj5fNQ/+xwNr76KNiqK5P/7PyLGj1N9D4pymFLBcAhsWbyQuc8/QURcAhNv/ivJefkdvs0DVeuq5f7l97Nwz0L6x/fn3mF3kr9pNix+DBD4jruFZ8pO5KXtjYSA00wu0ks/QbjtRCenUHDiKQTiklj47WIMWsFZhiX0si9B9j2X6f7rKFvsIAhoh8Zy7dT+fFW+iL+XhWgglisSgtzdZyhmbbhPxbN1K1V33Y1n0yYixp5C0v/9H/rExP1UryhKVzhsgkEIMQF4EtACr0gp//OTv98KXAUEgDrgCill2f7W2dHBIKVkxacfsGTa22T07c8Zf/k7ZtvPHyTralJK5pbO5d8r/k2rr5VL+17KDdmnY57/AGz+DKKzqBv8d/6xKJHZTjdxWg1T8zQk7l5A5ZaNCKEhecAQ6owRNNntFKZqGVf9PAatpKrv7fxv2RBMrSGqLHDK1N70znVx27r5zAsOI13n4rG+fTkhLjpcSyBA45tvUvfU0wiDgcTbbyN60qQO6ZBXFOW3OSyCQQihBbYDY4FyYBUwRUq5ea82JwMrpJQuIcT1wElSygv2t96ODAYpJUs/eIfln0yjz/EnM/66P6HV6X/9g12o2dPME989wcc7PibRksgdx9zB2JABMfvO8JPTmSPYGHEtf/vOxEYZpHekiT+PScNaspJNC+dhb24ilJaDKyKWmAgrk+K2klb6EaGoLOaZ72XL+iiCUtLa28YNl+czp/hV/tvUjxqRwrnxRu7vlU+cIfzgm6+0lKp//BPXypWYBgwg+R//wNyvoIuPkKIocPgEw3HAPVLK8W3v/wYgpXzwF9oPBp6RUo7c33o7MhgWv/82Kz6dRv/R4xh79Y1H1L9419Wu44EVD7C1cSsjUkfw92P+StaOBbDwP+CsI5R/JjPrz+WBKiM1SEZlRPPXM/sQUV/MxgVfsXXzJlxJmUidnp7xZs4Ofo6leSONyWcyrWwqoSYt5cYQA87M4djeG3l4+3o+5zQitIJ783OYnByDEAIpJa0zZ1Lz8MME6xuIvmAyibfcojqnFaWLHS7BMAmYIKW8qu39pcCxUsobf6H9M0C1lPJf+/jbNcA1AJmZmUPLyvZ7tek3+W7WZyx482X6jxnP2Kv+eESFwvcCoQDTtk3jmbXP4A16uazgMq7qeSHWVa/B0qcg6MObPYVXd43lBa+eViQT+yZx22l9SDIE2LBgHotXrMRhMKPxeelrcjLS/yVJmjrWRP+dZTsGIUNQmqhjysUx7G5+jiedJ7FD9GZ4pIl/98qiry1851LQbqf+mWdofOddtBERJPzlVqLPO++IPK6KcjQ44oJBCHEJcCNwopRyvwMMdcQZw/bli/n8iYfoUTicM269s0vvPDoU6lx1PL7mcT4v/pxYUyw3Db6Jc5KOQ7voEfjuLaTeTGvcxTxXOpK3pA6fgAsKM7h5bE8SI4ysWrSQed98iy8YRN9YS6q/noHGzWRGCRbIv1Fbl0KTJoSnn4VJJy3hw5oSpnERLqxMTUvgrznJxOrDl5c827ZRfd/9uNeswVRQQOIdf8U6bFgXHyFF6X4Ol2A4oEtJQohTgKcJh0Ltr633UAdDTfFO3vvH7STl5jPp7vvRG4yHbN1draiuiEdWP8La2rXkx+Rze+HtHKePha/vg60zkcZo6g2TeLzuBD5Ah1ar4eLhWVx3Yi5RRsHcOXP4bu1adDKErqIYo7OZfFstydGJfOe7Fr/Xyi5jiL7j3GTGvs7b7mP4WownUqfjjtxULk2JQ6f58fJS7WOPE6iqwjZ6NIm33YYxN6erD5GidBuHSzDoCHc+jwEqCHc+XySl3LRXm8HAR4TPLHYcyHoPZTC4HXbeufMWpAxxyYNPYImM+vUPHWGklHxV9hWPrXmMCkcFJ6SfwK1DbyXP1QoLHoQdc5GGWCoC5/CI6yRmokWn1TDl2EyuPykPb3MtX3zxBTU1NcSYjLBzIwF7C1F6D9GR6VQHzyIgItkRHWDCacuoDyzmHc3VbJK96Gs1cU+PNE6IDd/VFfJ4aHzzLRpeeomQ10vMBRcQf+Mf0cXEdPFRUpSj32ERDG2FnAY8Qfh21deklA8IIe4DVkspZwgh5gH9gaq2j+yWUp65v3UeqmCQoRCfPnQvuzeu58J7H+6Sp5k7kzfo5X9b/sdLG17CFXBxeu7pXD/wetJbqmHBv2HX14T0sezxnsvjvpP4HA1arYYpwzK4+vgcyrcXMX/+fILBIL0z0wluWUFFSTkCidUcj1czCo8hi/q8KsYO+x+rgvG8r72W2lAkJ8VEcFdeCv0jwkN0BBoaqHvmGZo/+BCN2UzclVcQc+lUtDZrFx8lRTl6HTbB0BEOVTCs+WI6C996hVOuuoGBY087BJUdGZo8Tby28TXe2/oeQRlkUv4krhlwDQn1u8IBUfINIV00lZ7TeSp0Cp9IAxqN4PzCDC4emsSWVd9SVFREZGQkw4cMJrDuSzavXofTr0WjMSEMA2g198Jy3AYKMmczn5P5THMhrSE95ybFcEdOMlnm8OU6765d1D72OI6vv0YbHU3c1VcTc9EUNGZzFx8lRTn6qGD4FQ3lu3n7zpvJHjiEs267u1sO41DjrOGlDS/xyY5P0Gl0TOkzhSsKriC6Zgssfjx8iUlYqPGP5zkm8l7ISlBKTu2fwhk9zOxe9y2VlZUkJyczdsxo5PrPKZo3m+ImIyDQ6DJwxmaSOXo9MVGbmKO/jC+CYwii4bK0OP6UmUSiMfyMiLuoiLonn8K5eDHahHjir72O6MnnozEY9r8TiqIcMBUM+xEMBPjf3X/BXl/HZY8+izW6e1/f3tO6h+fWP8cXxV9g1pm5oPcFTO07lfiWalj8OHLTJyC1NARG85bmTF4PxWEPhDguN47xmYLGzUtpbW2hZ8+ejB0zGtOO2Wyc/iZrKyNwByUII/rcKLJGbcJvCvKF8Wa+9PXHoBFMTY3nj5mJPwSEa/Vq6p54Etfq1ehSUoi7+iqizz0XjcnUxUdJUY58Khj2Y8WnH7D4/bc489a/k3/siENY2ZFtZ9NOXip6ibmlc9Fr9EzqOYk/FPyBZK8LljyFXPsuhPw4gsOYpTmdx0QvajwBeifZOCklSKB4JUG/lwEDBnDiqJHEVC6gZMaLLCmOptYVQGiCRPf3kzq0nAZ9HF9ZbuUrdw76vQIiyahHSolr2TLqnnwK9/r1aOPjifvDZURfeCFam62rD5OiHLFUMPyCltoa3vjLDW2XkO46xJUdHUpbSnl146vM3DUTBJyVdxZX9ruSDKGHVa8gV7yK8DbiDeWwXHMm/9EexxZniDirnuHxASLrNmCSXgYPHswJxx9PdP1qmr98hgUbrJTatQhdHUlDG4kvaKRGZDI/4q/MdSaj1wguTY3j+oxEUk2GcECsXEXDiy/iXLoUTWQkMRdfROzUqeouJkX5DVQw/ILpj9zP7qL1/OGx54mMTzjElR1dKh2VvLbxNT7d8SkBGeCUzFOYWjCVgdE9kRs+QC56Fk3LNoIymu3yVJ6zTmBmkx6tRjAwTpJi30aCcFBYOJSRI0cSbd+Oe/5zLF4l2NIai8a0hZRjK4jOcVAZzObriNv5yp2MEHBOUgw3ZCTSp+0pandREQ0vvYT9q3kIs5no884jduqlGDIzu/goKcqRQwXDPuxcvYLPHrmfEy6+nGPOPK8DKjs61bnqeHvL23y0/SPsPjsDEwYyte9URmecjK5sCaEFT6PZMw8pdTTK4XwWOZHH7bnYvUHSrZJsfxm52kaGDOzPyJEjSZAN+Ja8zLqljaxpyAfbOpKHbCUy00mlL51FlpuZJ/NwhySjYyP4Y2YiI6JtCCHw7txJw8uv0PLFFxAMYhs9mtipU7EMO6Zb3kCgKAdDBcM+7Fi5lLVzZnLe3+9Dq9N1QGVHN5ffxfSd03lnyzvsse8hzZbGxX0u5pwe52Cz1xBc9AJi4zQ0wRZ8oVTWWCfyqDyeNc06zFrI0dSRJ2oY2TeLkSNHkhFnJbD6HbZ+vZFV1YPwR2wgceAyItIcNDpjmS+uZGHEcJpCMDDCzPUZiUxMiEavEfhraml6/z2a359GsKkJY+/exE6dSuTpE9WdTIryC1Qw/AIppfqX5e8UDAVZWL6Qtza9xXe132HWmZmYO5ELel1A74hsQms+IrTkZXT2dUipo1J3PB/YTuP5+jR8QUjQuughajkx28pJo4bTMy8Xts9l15wFrCvJpiWikoR+87ClNOHymPiy6RwWJ51OldZEkkHHpanxXJoaR5JRT8jjoXXmTBrffAvvjh1o4+KIPvdcoiefjyEjo6sPlaIcVlQwKJ1iU/0mpm2bxuyS2XiCHgYkDOCCXhcwLmscxpqdBL56AW3ZJ2ikEx+prIsYz9OB4XzbFIGOENnaBgZHeThrRD+GDh2KuWUntfM+Zv13GvYYo4jp9TmRmSUE/BqWVIxiUdQ5bI3NRAdMTIzmirR4hkWFn5Z2LV9O49vv4Fi4EEIhrCNHEn3BZCJOPhmhP7zn1FCUzqCCQelULd4WPt/1OdO2TaO0tZQoYxRn5Z3FufnnkmtJJrDoPVjzP/TuNQA06wYw1zaWh+sLaAiYsAkP+fpmJhYkcubJx5IYbcW5cjobv97FNlci1vwFRGWvASQ7yvP4MnAa67JG4NYZ6G02cGl6AuclxRCt1+Gvrqb5o49p/ugjAtXV6BISiDrvXKInnY8hPa1rD5SidCEVDEqXkFKyqnoV07ZNY/7u+QRkgP7x/Tkr7ywm5EzAVleHf97r6Mo+RScrCGGk1HIC72pG8WZ9DgF0xAsHQ+NDXDgin5FDB6Ct207xrC/ZvMuDN6OCqNxv0BndNDVGM79+NEviTqImPg29lEyItTE1M5mRMTZEMIhj0bc0T5uGY9EikBLLsGFEnX02EePGqXGZlG5HBYPS5RrcDcwqmcX0ndPZ3rQdvUbP6MzRnJV3FsOTjiW0ahFyxdsYmr9CI5wERDSbbKN4zTWYz529kAgydHZOzIngktED6JkWT+vKOWxctJ1yXTO2vOWYYvYQ8OpYXz6IRaHj2Zg5BI/JQpIMMCUljkty0kg3GfBXVtLy2Wc0T5+Ov2w3wmQiYuxYos4+C+vw4QjtkT3/hqIcCBUMymFla+NWPtv5GV8Uf0GTt4kEcwLjs8czIWcCBbps/PM/QWydjsG3DI3w4tPEsdI0gpcdQ1jky0cAOWYv4/rEc+noQSTpvZTMmcvWXdvxp24nIn0tQhOiqS6Bb5pOYKVtGGVpPQAoCHo4Py2R8/MyidVrca9bR8v0z2idPZtQayu6xEQizzidyFNPw1TQV92coBy1VDAohyV/0M+iikXM2DmDxRWL8YV8pNnSwiGRPYE8Rxz+bz9Bu2sGxtBKhPDj1iayTF/IO46BfBPoi0SQZwswoSCJS08egK2xjG1ff8Ue92Z0mZswRlYTCujYXtmXb33HUxTXj4bYRDShEIN9ds5LiWNyQS8soRCOBQtpmT4dx+LFEAigz8ggcsIEIk+dgLFPHxUSylFFBYNy2LP77CzYs4DZJbNZXrmcgAyQHZnN+OzxjMscS2aNmcDST9Du+QITaxHCj09EsFo3iA9cg5gbHIwHI3kRIcb0SWLyyD7EVG1my4q5NOg3Y0orQqv34nNEsqFmEEu1I9ic0BO7LQp9wM8QVyMTE6OZPKA/EX4f9q+/pnX2HJzLlkEwiD4rk8gJpxI5YTzG3r1VSChHPBUMyhGl2dPMvN3zmFMyh1U1qwjJEGm2NEZnjubklJPoU5+AXDUHbflXmFiBRjgIYqBI25fPvIOYGTiGOmJIMoU4oUcM5x6bT1rdFoq3Tsdp3ooxYSdCI/E0J7C68RiWGoazPSEbt8mKLuCnb2sdp1j1TOqTT0ZEJI6vv8Y+ezbO5SsgFEKfloZt9GgixozGMnSouv1VOSKpYFCOWPXuehbuWcj83fNZXrUcf8hPtDGaE9JPYEzaaIbaeyK++xZN6VxMcik6EZ4mvFyTyYJAX2b6C1kje2LUaTk208bpg7Pp7d5GXekMnObNGON2A+BuSGNt82CWG4awPT4XuyUCTTBITl05I4WPMzJTGZyRiX/ZMhzzF+Bctgzp9aKJjMR2wglEjD4Z6/HHo42I6MKjpSgHTgWDclRw+p0sqVjC/D3zWVS+CLvPjlFrpDC5kOOTRjHSP4SYLSWI4nkY/SsxiC0IEcSDmVWyF3MCQ/g6OJgaYsmN1nFi7ySG2pow136Cz7wOfVR4RllfSxJbmgazWDOYLXH5NNnCc38nNFTT317PydFWTsnJIqamDtc3i3AsXEiwsRH0eiyFQ7GNGoV11CiMPXuqS07KYUsFg3LU8Yf8rKlZw8I9C1lSsYTS1lIA0mxpjEwZyWjTSHrXxKDb8i3apsWYNGvQiXoAyklmabA3C4IDWRbqi1sbwZA0KyMydKQ7vsSsX4wuZjdCSHyuGPY0DGJFsD9FET0pj00ipNFg9LrJrSxmaMDF2MQY8rUGLFu34lq0GO+OHQDoEhKwjhyJddQorCNHqOHBlcOKCgblqFduL2dp5VIWVyxmRdUKXAEXOqFjUOIgRsUcx3DHAFJK69DuWYAxuA6DZiMa4SWEoESmsiRYwMLQAFaFeqM1WhmSbqaXcSep2m9ISliFVhcgFNDT0tCL9Z5BrNb3ZntcFk6TBYD4hmqya3ZTqA1yfKSVTKcb0+atuJcuI9jSAkJg6tsX68iRWIYNwzJ4EBqreqhO6ToqGJRuxR/0s65uHYsrFrOschlbG7cikZi0JgbHD+YUw/EMdOSSUFaCtnY5BtZj1GxBI/yEEGwLZbIs1IfloT6sCfUkYIymf6KPfMNmsiKXkh6/BSHA64ynpOkYNsheFFly2B2TTECrQ4SCJNdVklOzm6GaAMdqIaGpBcvGrQSKiiAQAJ0Oc79+WI45JhwUQwaroFA6lQoGpVtr8bawuno1K6tXsrJ6JTubdwJg0VkYFncMozmOvk1pxJfvRN+yCp1mLSaxC40IALBHJrA62Ivlsg+rQz2p1afSM9pBlmkrObFryI3diUnnw96cyU5nIRvJZ6M1i4qoBKRGgzbgJ7mugvTq3fQNuCkM+Ui3O7HtKkVs3IT4PigKCrAMOwbz4CGYBw9Sl56UDqWCQVH20uhpZFX1qh9exS3FAGiFlsERAxnNCAY4s0ipqcTQtB6NKMKk2YZBOABolRbWhvJZGerFBplLUSiHqAhBrrWUnNgN9IjZRYK5npaWXLa4j2WbJpdt5lQqIhMIaTRA+NJTWnUZea31DA54yWltIaK8CvPmreh9fgAMWVmYBw3CPHgQ5sGDMfbooYbrUA4ZFQyKsh8t3hbW161nXe061tetp6i+CHfADUCWMYNx+hMZ5OpFdoMHW8MmYANGzRasouaHdZTLBNaF8lgbymNDKI9SXTqJtgayIneQFV1GduQejEEj25zHUix7sM2QRnFkMj59eCIhq7OVlNoKkuvKyXO10s/RQkpTE9bi3dhq69CFJBqrFdOA/uGwKCjA1K8fuqQkdeeT8puoYFCUg+AP+dnetD0cFLXrWV+3nkpnJQA6qWWk9hiOYwh9nQmkNFSjdW0G7RYsohiraAYgKAUlMoV1oR5slllsllmUaZKIiWgiO7qErIg9pFurafL1YkegP2XaFHaZk6iz/Xj5KLqlnpTaClLqKsh3ttC7qZ6YunosuyuIdHkwBoJo4+Iw9e2LqaAvpoICzAUF6FJSVFgov0oFg6L8To2eRjY3bGZzw2Y21W9iU8MmalzhM4aIoJUTNcMpDPajp91KfPMeCGxGp9mBTZRibrsEBVAlY9kcCgfFllAmJSKZoC1AakQV6RFVxJgcOEUmFfSgTJtCiTmBFnPbQ3NSEttcT0JDNYkN1aS1NpDf3EB6dTXW6loiXR6sHh/62NhwWPTpjbFnL4w9e2LMyUaoaU6VvahgUJQOUO+ubxcU25u2U+WsAgnxgWj6BXpRSD/6uqNItNcjfDvQanZhFnuIog6tCAHglga2yQy2h9LZIdPYKdOoM1jRWr2kRVYTZXbi0SbRaEijWpdIuTGeBnPUD3WYPU4S6qtJaKgmqbGG7NZGsqsriamuxer2EOHxYQlKzDk54ZDo2RNjz3xMPXuiS01VZxfdlAoGRekkrb5WdjTtYEfTDrY3bWd703Z2NO3AFXBhC5rJ9WQwONSXfv4sspxuTN4ytOzCKMqwihpse51duKWeXTKVnTKNHaF0dskUGg0W/NYQ0bZWtAYdHl0cLaZEasxxVFni8Gt/HLcp0t5EXFMdcU21xDXXk9ncQFZVOQk1NVg9fiI8Pmx6A9bcPAy5uRhycjDm5mDIycGQkaHOMI5yKhgUpQuFZIgKR8UPIVHSUkJJSwmlraV4/R6S/PFkeVPoGcimvz+JNI8Ho7cCrSjDKCqwiVqi2/ouAAJSQwXxlIRSKJVJlMkkqjQxuE1a3DYd0mTAp4/EaYqm0RpDnTW2XWBYXHbim+qIba4juqWBxJYm0uprSKsoJ9rtxuL1Yw0EiY5PwvLTwMjMRBsfr84yjgIqGBTlMBSSIaqd1T8ERXFL8Q8/m93NpPjiyfKmkuZLpFcgkWxfkEh/HXpZgU5UYxR1RNOAWXh/WGdAaqiQ8ZTJJEpkMrtlEvWaCFqNRpotNhyGKLymSBzmCBojo/FYLfD9/8lLSaSjheiWBqJbG4lpaSCxuYHUuirSqquJcTox+wJYhJaohAQi0jIwZmagz8jEkJGOPiMTfXoaGnWmcURQwaAoR5gWbwt77HvavXa37qbcXo7TYSfNl0iaL5EUTwI9AhGk+wNEBBvRU42Wasyinljq24UGgF2aKJcJVMh4ymUCVcTQqLfRaIik2hhPvSkej9mG3WLDbbUhTVrQhoPD7HYSaW8m0tFMpL2ZaHsT8c2NJDXUklJfR1xrKxZ/EFtEJFGJSeHgSEtFl5KCPjkFfWoK+uRk9YT3YUIFg6IcRdwBNxX2ih8Dw76bSkclVY4q7K3N2Fxmkv1xJPviSPeZSPNDXMiNJdiMjjq0ohGLaCKGJiw/CQ6XNFItY6iRsVQT/tmgsdFgiKDeEEu1IYFKSxKt1ij8JiMYtUijBnQaDF4PkY5mouzNRDhbiHC0EtPaTGxLOEASmxqIsTuw6A3YoqKJSEgkMjUNU1oa+pRU9MlJ6BIS0CUkoLFYuujodh8qGBSlG2n1tVLlqKLKudfLEX45m1vR2SHJF0uCL5Zkv5Fkf5DYkIuIkB2jbEQjmjGIZiJEC1E40BNst/6QFNQT9WOAyBgaRCRNOhuNuihqDXFUGxOotCbiNNt+CA9p0KLDj83lwOZsJcLZitVpJ8reTGxrC7GtzcS1NhHb3EKs34fVYsUSFY01Lh5bYhKW5BT0iYnoEhPDAZKYgMZqVf0dv5EKBkVRfuAP+ql311PrrqXOVUetq5Y6d/hng6MeX6sb7AFMLj3x3ihS/SYS/UFiQj6sIQd6aUcrWjBqGrHSShR2InDvc1suaaReRlJPFPUyigYiadRE0KCNplYXQ50hhkZDDI3GKJpM0XhMJqRBi9BJzCE3Fp8Tq8uJxe3E6nYQ6bAT5Wglxt5KdGsz8U4H8cEAkUYj5ohIzFHRmGPjsCYkYEpIQhcXiy42Fm1sHLrYGHUmspf9BYOus4tRFKVr6bV6UmwppNhS9tvOE/BQ766nzl1HjauGSncDTZ4mGt2NOJyt+Fs9SEcAjUNic+tICuiI84eIDvmxhbwYpAstDnSaVizCTh92EyFdROBCG5LgI/zai1MaaSKCZmmjSdpoxkazxkaTJpJGXSSN+iiaEuPZnZZHkyGSFmMErcYIpBaM+DCHPJgCbsxeN+aqeqzFZdhcTiKcDiKddqLcbqKDfuJCIeI0YLPYMEVGYIqKxhwdiykuDmNcHLqoKDRRUWjbXhqTqeO+kMOQCgZFUfbJpDORHpFOekT6r7b1BX3h0PA0tns1uBtwO5z4HV6CTh+4gmidISLdfmL8AaKCIWyhAMaQB51wo9c40ePCgItUGulBJVbpxiY9aPwS/IDn59t3SSOtWLBLC61YaJUW7BoLLRobrRobrVEWWuNs7NEmskkfQYvWRqvBhltrwKfRERBatCKI0enF1FKB2bsLs8eD1evG4nZh87ix+bxEBoNEEiJGCKK0gmijAbM1AmNEJMbISIzR0RijYtBHRaKLiEATEYHGZkNrtR5Rz4WoYFAU5XczaA0kWZNIsiYd8Gc8AQ+tvlZavC3hl6+FWm8rLZ4WnA47frsHv8tD0BVEOAOYPB4sHh9WvwdLwIdJ+jEIDzrhRYcHPV700osZDzHSjgkfFunFGvSgDwbDobK/eqQeB2Zc0ogTMw5MOIUJp96M02DGqQm/GjVm9mjNOLQWXFoTbo0Jj8aAzx/AX99EsKGFAIIQGkQwhCEQwOT3YvF5MXu8WPx+LAE/VhnEJkNECIjQaonUaYk0GrAZTZgtFvRWK4YIG4aISAyRUegjItFarWgsZjQWCxqzGaHX73+nfiMVDIqidAmTzoRJZyLRknjQnw2EAjj9Tuw++w8/HX4HDp8Dl9uJ1+nG7/Tgd3nxO73g8mL0uDF6vRj9HvR+N/qQBy0+dNKHTvjREX7ppQ8dfizSR5R0YsSHSfowB31Ygl50hA6qVo/U48aIGwNujHgMBjwGA25hDL8XBtwaE7XCSJkw4AkZ8biMeN16vI16/EKHv+2sJiC0+IWWEBpCQE+3h7vuuOegj9+v6dRgEEJMAJ4EtMArUsr//OTvRuAtYCjQAFwgpSztzBoVRTn86TQ6ooxRRBmjfr3xrwiEAngCHlwBFy6/C1fAhTvgptHnwu1x4XN78Lm9+DwefE4PIbcHjcuF1utC5/Og9bvRBNyIkB9N0IeQAXTSj0YG0MkAWgJoCKCVQbQE0RFALwPEYCdJNmHAjzHox4gfEz7MP+142Y8XLOf+7v3fl04LBiGEFngWGAuUA6uEEDOklJv3anYl0CSl7CGEuBB4CLigs2pUFKX70Wl02Aw2bAZbh20jGAriDXrxBD14A20/g16cPjdNPi8+rw+/z4vf58Pn9iFdTqTbifS6wy+fm5DXgwx4IeBFhPyIkJ/cvKEdUm9nnjEMA3ZKKYsBhBDvA2cBewfDWcA9bb9/BDwjhBDySL+nVlGUbk2r0WLRWLDoj4zbZTWduK00YM9e78vblu2zjZQyALQAcT9dkRDiGiHEaiHE6rq6ug4qV1EUpXvqzGA4ZKSUL0kpC6WUhQkJCV1djqIoylGlM4OhAsjY631627J9thFC6IAowp3QiqIoSifpzGBYBeQLIXKEEAbgQmDGT9rMAC5r+30SMF/1LyiKonSuTut8llIGhBA3AnMJ3676mpRykxDiPmC1lHIG8CrwthBiJ9BIODwURVGUTtSpzzFIKWcBs36y7B97/e4Bzu/MmhRFUZT2jsjOZ0VRFKXjqGBQFEVR2jni52MQQtQBZb/x4/FA/SEs50ig9rl7UPvcPfyefc6SUu7zfv8jPhh+DyHE6l+aqOJopfa5e1D73D101D6rS0mKoihKOyoYFEVRlHa6ezC81NUFdAG1z92D2ufuoUP2uVv3MSiKoig/193PGBRFUZSfUMGgKIqitNMtg0EIMUEIsU0IsVMIcWdX13OoCCEyhBALhBCbhRCbhBA3ty2PFUJ8JYTY0fYzpm25EEI81XYcNgghhnTtHvx2QgitEGKtEGJm2/scIcSKtn2b1jZwI0IIY9v7nW1/z+7Swn8jIUS0EOIjIcRWIcQWIcRxR/v3LIT4c9t/1xuFEO8JIUxH4/cshHhNCFErhNi417KD/m6FEJe1td8hhLhsX9v6Jd0uGPaaYvRUoC8wRQjRt2urOmQCwF+klH2B4cAf2/btTuBrKWU+8HXbewgfg/y21zXA851f8iFzM7Blr/cPAY9LKXsATYSnjYW9po8FHm9rdyR6EpgjpewNDCS870ft9yyESAP+BBRKKfsRHojz++l/j7bv+Q1gwk+WHdR3K4SIBf4JHEt49sx/fh8mB0RK2a1ewHHA3L3e/w34W1fX1UH7+hnhOba3ASlty1KAbW2/vwhM2av9D+2OpBfhuT2+BkYDMwFB+GlQ3U+/c8Kj+x7X9ruurZ3o6n04yP2NAkp+WvfR/D3z4+yOsW3f20xg/NH6PQPZwMbf+t0CU4AX91rert2vvbrdGQMHNsXoEa/t1HkwsAJIklJWtf2pGkhq+/1oORZPAH8FQm3v44BmGZ4eFtrv1wFNH3uYywHqgNfbLp+9IoSwchR/z1LKCuBRYDdQRfh7W8PR/T3v7WC/29/1nXfHYDjqCSFswMfALVLK1r3/JsP/fDhq7lEWQpwO1Eop13R1LZ1IBwwBnpdSDgac/HhpATgqv+cY4CzCoZgKWPn55ZZuoTO+2+4YDAcyxegRSwihJxwK70opP2lbXCOESGn7ewpQ27b8aDgWI4EzhRClwPuELyc9CUS3TQ8L7ffraJg+thwol1KuaHv/EeGgOJq/51OAEillnZTSD3xC+Ls/mr/nvR3sd/u7vvPuGAwHMsXoEUkIIQjPgrdFSvnYXn/ae8rUywj3PXy/fGrbnQ3DgZa9TlePCFLKv0kp06WU2YS/y/lSyouBBYSnh4Wf7/MRPX2slLIa2COE6NW2aAywmaP4eyZ8CWm4EMLS9t/59/t81H7PP3Gw3+1cYJwQIqbtbGtc27ID09WdLF3UsXMasB3YBdzV1fUcwv0aRfgUcwOwru11GuFrq18DO4B5QGxbe0H4Dq1dQBHhOz66fD9+x/6fBMxs+z0XWAnsBD4EjG3LTW3vd7b9Pber6/6N+zoIWN32XU8HYo727xm4F9gKbATeBoxH4/cMvEe4H8VP+Ozwyt/y3QJXtO3/TuDyg6lBDYmhKIqitNMdLyUpiqIo+6GCQVEURWlHBYOiKIrSjgoGRVEUpR0VDIqiKEo7KhgURVGUdlQwKN2SECJOCLGu7VUthKho+90hhHiuA7b3hhCiRAhx3SFY1yNtNd92KGpTlJ/S/XoTRTn6SCkbCD8khhDiHsAhpXy0gzd7u5Tyo9+7Einl7UII56EoSFH2RZ0xKMpehBAniR8n+7lHCPGmEOJbIUSZEOJcIcTDQogiIcSctnGpEEIMFUJ8I4RYI4SY+/2YNr+ynTeEEM8LIZYLIYrbtvuaCE+680ZbG21bu41t2/xzh+68orRRwaAo+5dHeGC+M4F3gAVSyv6AG5jYFg5PA5OklEOB14AHDnDdMYTnEPgz4TFvHgcKgP5CiEGEz2jSpJT92rb5+qHaKUXZH3UpSVH2b7aU0i+EKCI8a9ictuVFhCdT6QX0A74Kj+2GlvA4NwficymlbFt3jZSyCEAIsalt3d8AuUKIp4EvgC8PyR4pyq9QwaAo++cFkFKGhBB++ePgYiHC//sRwCYp5XG/dd1t6/LutTxEeFayJiHEQMIzlV0HTCY8MJqidCh1KUlRfp9tQIIQ4jgIz4chhCg4FCsWQsQDGinlx8DdhOdcUJQOp84YFOV3kFL6hBCTgKeEEFGE/zf1BLDpEKw+jfD0nd//A+5vh2CdivKr1LDbitIJ2u40mnkobldtW989dM4ttko3pC4lKUrnaAHuP1QPuAGXEJ7rWVEOOXXGoCiKorSjzhgURVGUdlQwKIqiKO2oYFAURVHaUcGgKIqitPP/RIQS2OB2abAAAAAASUVORK5CYII=\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "import matplotlib.pyplot as plt\n", + "plt.plot(model.t, model.output.T);\n", + "plt.xlabel(\"Time [ms]\")\n", + "plt.ylabel(\"Activity $x$\")" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.7.3" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/neurolib/models/model.py b/neurolib/models/model.py index cbbe0b4a..a8828007 100644 --- a/neurolib/models/model.py +++ b/neurolib/models/model.py @@ -14,6 +14,8 @@ def __init__(self, integration, params): if hasattr(self, "name"): if self.name is not None: assert isinstance(self.name, str), f"Model name is not a string." + else: + self.name = "Noname" assert integration is not None, "Model integration function not given." self.integration = integration @@ -33,6 +35,10 @@ def __init__(self, integration, params): assert isinstance(self.default_output, str), "`default_output` must be a string." + # if no output_vars is set, it will be replaced by state_vars + if not hasattr(self, "output_vars"): + self.output_vars = self.state_vars + # create output and state dictionary self.outputs = dotdict({}) self.state = dotdict({}) @@ -301,7 +307,8 @@ def clearModelState(self): self.state = dotdict({}) self.outputs = dotdict({}) # reinitialize bold model - self.initializeBold() + if self.params.get("bold"): + self.initializeBold() def storeOutputsAndStates(self, t, variables, append=False): """Takes the simulated variables of the integration and stores it to the appropriate model output and state object. @@ -397,18 +404,24 @@ def getMaxDelay(self): :return: maxmimum delay of the model in units of dt :rtype: int """ - dt = self.params["dt"] - Dmat = self.params["lengthMat"] + dt = self.params.get("dt") + Dmat = self.params.get("lengthMat") - if "signalV" in self.params: - signalV = self.params["signalV"] + if Dmat is not None: + # divide Dmat by signalV + signalV = self.params.get("signalV") or 0 if signalV > 0: Dmat = Dmat / signalV else: + # if signalV is 0, eliminate delays Dmat = Dmat * 0.0 - Dmat_ndt = np.around(Dmat / dt) # delay matrix in multiples of dt - max_global_delay = int(np.amax(Dmat_ndt)) + # only if Dmat and dt exist, a global max delay can be computed + if Dmat is not None and dt is not None: + Dmat_ndt = np.around(Dmat / dt) # delay matrix in multiples of dt + max_global_delay = int(np.amax(Dmat_ndt)) + else: + max_global_delay = 0 return max_global_delay def setStateVariables(self, name, data):