{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Custom observation models\n", "\n", "While *bayesloop* provides a number of observation models like `Poisson` or `AR1`, many applications call for different distributions, possibly with some parameters set to fixed values (e.g. with a mean value set to zero). The [sympy.stats](http://docs.sympy.org/dev/modules/stats.html) and the [scipy.stats](http://docs.scipy.org/doc/scipy/reference/stats.html) modules include a large number of continuous as well as discrete probability distributions. The observation model classes `SciPy` and `SymPy` allow to create observation models to be used in *bayesloop* studies on-the-fly, just by passing the desired `scipy.stats` distribution (and setting values for fixed parameters, if necessary), or by providing a `sympy.stats` random variable, respectively. Note that these classes can only be used to model statistically independent observations. \n", "\n", "In cases where neither `scipy.stats` nor `sympy.stats` provide the needed model, one can further define a custom observation model by stating a likelihood function in terms of arbitrary [NumPy](http://www.numpy.org/) functions, using the `NumPy` class.\n", "\n", "## SymPy.stats random variables\n", "The [SymPy](http://www.sympy.org/en/index.html) module introduces symbolic mathematics to Python. Its sub-module [sympy.stats](http://docs.sympy.org/dev/modules/stats.html) covers a wide range of discrete and continuous random variables. In the following, we re-define the observation model of the coal mining study `S` defined above, but this time use the `sympy.stats` version of the Poisson distribution:" ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "execution": { "iopub.execute_input": "2026-04-27T19:26:16.851518Z", "iopub.status.busy": "2026-04-27T19:26:16.851146Z", "iopub.status.idle": "2026-04-27T19:26:17.842801Z", "shell.execute_reply": "2026-04-27T19:26:17.842413Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ " + Trying to determine Jeffreys prior. This might take a moment...\n", " ! WARNING: Failed to determine Jeffreys prior. Will use flat prior instead.\n" ] } ], "source": [ "import bayesloop as bl\n", "import numpy as np\n", "import sympy.stats\n", "from sympy import Symbol\n", "\n", "rate = Symbol('lambda', positive=True)\n", "poisson = sympy.stats.Poisson('poisson', rate)\n", "\n", "L = bl.om.SymPy(poisson, 'lambda', bl.oint(0, 6, 1000))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "First, we specify the only parameter of the Poisson distribution (denoted $\\lambda$) symbolically as a positive real number. Note that providing the keyword argument `positive=True` is important for SymPy to define the Poisson distribution correctly (not setting the keyword argument correctly results in an error). Having defined the parameter, a random variable based on the Poisson distribution is defined. This random variable is then passed to the `SymPy` class of the *bayesloop* observation models. Just as for the built-in observation models of *bayesloop*, one has to specify the parameter names and values (in this case, `lambda` is the only parameter).\n", "\n", "When creating a `SymPy` observation model, *bayesloop* attempts to derive the [Jeffreys prior](https://en.wikipedia.org/wiki/Jeffreys_prior) symbolically. With the current SymPy stack, this derivation can fail for some distributions; in that case *bayesloop* falls back to a flat prior and reports this in the notebook output. This behavior can be turned off using the keyword argument `determine_jeffreys_prior`, in case one wants to use a flat parameter prior directly or avoid the symbolic calculation:\n", "```\n", "M = bl.om.SymPy(poisson, 'lambda', bl.oint(0, 6, 1000), determine_jeffreys_prior=False)\n", "```\n", "Alternatively, you can of course provide a custom prior via the keyword argument `prior`. This will switch off the automatic determination of the Jeffreys prior as well:\n", "```\n", "M = bl.om.SymPy(poisson, 'lambda', bl.oint(0, 6, 1000), prior=lambda x: 1/x)\n", "```\n", "See also [this tutorial](priordistributions.ipynb) for further information on prior distributions. Having defined the observation model, it can be used for any type of study introduced above. Here, we reproduce the result of the [regime-switching example](changepointstudy.ipynb#Exploring-possible-change-points) we discussed before. We find that the parameter distributions as well as the model evidence is identical - as expected:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "execution": { "iopub.execute_input": "2026-04-27T19:26:17.844177Z", "iopub.status.busy": "2026-04-27T19:26:17.844066Z", "iopub.status.idle": "2026-04-27T19:26:18.017256Z", "shell.execute_reply": "2026-04-27T19:26:18.016720Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "+ Created new study.\n", "+ Successfully imported example data.\n", "+ Observation model: poisson. Parameter(s): ('lambda',)\n", "+ Transition model: Regime-switching model. Hyper-Parameter(s): ['log10pMin']\n", "+ Started new fit:\n", " + Formatted data.\n", " + Cached observation likelihoods (0.8 MiB).\n", " + Set uniform prior with parameter boundaries.\n" ] }, { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "662661a372b64910af20a9ca71eee34b", "version_major": 2, "version_minor": 0 }, "text/plain": [ " 0%| | 0/110 [00:00" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "%matplotlib inline\n", "import matplotlib.pyplot as plt # plotting\n", "plt.style.use('seaborn-v0_8-whitegrid') # plot styling\n", "\n", "S = bl.Study()\n", "S.load_example_data()\n", "\n", "T = bl.tm.RegimeSwitch('log10p_min', -7)\n", "\n", "S.set(L, T)\n", "S.fit()\n", "\n", "plt.figure(figsize=(8, 4))\n", "plt.bar(S.raw_timestamps, S.raw_data, align='center', facecolor='r', alpha=.5)\n", "S.plot('lambda')\n", "plt.xlim([1851, 1962])\n", "plt.xlabel('year');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally, it is important to note that the `SymPy` module can also be used to create random variables for which some parameters have user-defined fixed values. The following example creates a normally distributed random variable with a fixed mean value $\\mu = 4$, leaving only the standard deviation as a free parameter of the resulting observation model (which is assigned the parameter interval ]0, 3[):\n", "```\n", "mu = 4\n", "std = Symbol('stdev', positive=True)\n", "\n", "normal = sympy.stats.Normal('normal', mu, std)\n", "L = bl.om.SymPy(normal, 'stdev', bl.oint(0, 3, 1000))\n", "```\n", "\n", "## SciPy.stats probability distributions\n", "We continue by describing the use of probability distributions of the `scipy.stats` module. Before we show some usage examples, it is important to note here that `scipy.stats` does not use the canonical parameter names for probability distributions. Instead, all continuous distributions have two parameters denoted `loc` (for shifting the distribution) and `scale` (for scaling the distribution). Discrete distributions only support `loc`. While some distributions may have additional parameters, `loc` and `scale` often take the role of known parameters, like *mean* and *standard deviation* in case of the normal distribution. In `scipy.stats`, you do not have to set `loc` or `scale`, as they have default values `loc=0` and `scale=1`. In *bayesloop*, however, you will have to provide values for these parameters, if you want either of them to be fixed and not treated as a variable.\n", "\n", "As a first example, we re-define the observation model of the coal mining study `S` defined above, but this time use the `scipy.stats` version of the Poisson distribution. First, we check the parameter names:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "execution": { "iopub.execute_input": "2026-04-27T19:26:18.018551Z", "iopub.status.busy": "2026-04-27T19:26:18.018448Z", "iopub.status.idle": "2026-04-27T19:26:18.020596Z", "shell.execute_reply": "2026-04-27T19:26:18.020273Z" } }, "outputs": [ { "data": { "text/plain": [ "'mu'" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "import scipy.stats\n", "\n", "scipy.stats.poisson.shapes" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In `scipy.stats`, the rate of events in one time interval of the Poisson distribution is called *mu*. Additionally, as a discrete distribution, `stats.poisson` has an additional parameter `loc` (which is **not** shown by `.shapes` attribute!). As we do not want to shift the distribution, we have to set this parameter to zero in *bayesloop* by passing a dictionary for fixed parameters when initializing the class instance. As for the SymPy model, we have to pass the names and values of all free parameters of the model (here only `mu`):" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "execution": { "iopub.execute_input": "2026-04-27T19:26:18.021750Z", "iopub.status.busy": "2026-04-27T19:26:18.021660Z", "iopub.status.idle": "2026-04-27T19:26:18.141894Z", "shell.execute_reply": "2026-04-27T19:26:18.141361Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "+ Observation model: poisson. Parameter(s): ('mu',)\n", "+ Started new fit:\n", " + Formatted data.\n", " + Cached observation likelihoods (0.8 MiB).\n", " + Set uniform prior with parameter boundaries.\n" ] }, { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "16c5a2ee547c4a5c8c7242be6719f7f1", "version_major": 2, "version_minor": 0 }, "text/plain": [ " 0%| | 0/110 [00:00" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "L = bl.om.SciPy(scipy.stats.poisson, 'mu', bl.oint(0, 6, 1000), fixed_parameters={'loc': 0})\n", "S.set(L)\n", "S.fit()\n", "\n", "plt.figure(figsize=(8, 4))\n", "plt.bar(S.raw_timestamps, S.raw_data, align='center', facecolor='r', alpha=.5)\n", "S.plot('mu')\n", "plt.xlim([1851, 1962])\n", "plt.xlabel('year');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Comparing this result with the SymPy example above, we find the same model evidence value. In the current dependency stack, the automatic Jeffreys-prior derivation for the SymPy model falls back to a flat prior, and the SciPy implementation also uses a flat prior by default. The built-in `Poisson` observation model still provides its analytic Jeffreys prior, so its evidence can differ from these generic model wrappers. Custom priors are still possible for SciPy and SymPy models using the keyword argument `prior`.\n", "\n", "## NumPy likelihood functions\n", "\n", "In some cases, the data at hand cannot be described by a common statistical distribution contained in either `scipy.stats` or `sympy.stats`. In the following example, we assume normally distributed data points with known standard deviation $\\sigma$, but unknown mean $\\mu$. Additionally, we suspect that the data points may be serially correlated and that the correlation coefficient $\\rho$ possibly changes over time. For this multivariate problem with the known standard deviation as \"extra\" data points, we need more flexibility than either the `SymPy` or the `SciPy` class of *bayesloop* can offer. Instead, we may define the likelihood function of the observation model directly, with the help of [NumPy](http://www.numpy.org/) functions.\n", "\n", "First, we simulate $1000$ random variates with $\\mu=3$, $\\sigma=1$, and a linearly varying correlation coefficient $\\rho$:" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "execution": { "iopub.execute_input": "2026-04-27T19:26:18.143063Z", "iopub.status.busy": "2026-04-27T19:26:18.142960Z", "iopub.status.idle": "2026-04-27T19:26:18.305506Z", "shell.execute_reply": "2026-04-27T19:26:18.305021Z" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAq8AAAFuCAYAAACvPfERAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjksIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvJkbTWQAAAAlwSFlzAAAPYQAAD2EBqD+naQAApmRJREFUeJztnQecHLX1x9+ejfudDW5gG1wx1WCHXkICgUCANCB/Qi8BE3oNYFrohN5rIAmht0ACCQklhNAhgMEG3LvPxr133/4/b/bmTquTNJJGM6PZfd/P5+y73RmNRtJIP715eioUi8UiEARBEARBEEQOqMk6AwRBEARBEAShC4lXgiAIgiAIIjeQeCUIgiAIgiByA4lXgiAIgiAIIjeQeCUIgiAIgiByA4lXgiAIgiAIIjeQeCUIgiAIgiByA4lXgiAIgiAIIje0hhywbt06WLx4MbRt2xZqakhvEwRBEARB+EZDQwOsXr0aOnfuDK1bt65u8YrCdcqUKVlngyAIgiAIgoigX79+0LVrV6hq8YoWV2SzzTaDjh07Zp0dImHWr18P48aNg8GDB0OrVq2ovCscqu/qguq7uqD6ri6WL18O06ZNa9JtVS1eQ1eBdu3aQYcOHbLODpFCZ4dgXZN4rXyovqsLqu/qguq7Ouu7JmEXz0wdSNesWQNXXXUV7LTTTrD77rvDbbfdBsViMcssEQRBEARBEB6TqeX12muvhY8++ggeeeSRwNR87rnnQq9eveCXv/xlltkiCIIgCIIgPCUzy+uiRYvghRdegGuuuQa222472G233eDEE0+EL774IqssEQRBEARBEJ6TmeX1008/hU6dOsHOO+/c9Nnw4cOzyg5BEARBEASRAzITr9OnT4fevXvDSy+9BA888ACsXbsWDjnkEDj11FOljr4YPyx0BiYql7COqa6rA6rv6oLqu7qg+q4uGhoaKlu8rlixAqZOnQpPP/003HDDDTB37ly44ooroH379oH7gIgJEyaknk8iO0aNGkXFX0VQfVcXVN/VBdU3URHiFXdeWLZsGdx6662BBRapr6+Hp556SipeBw0aFLgaEJU/U8eObsiQIRQqqwqg+q4uqL6rC6rv6mLZsmWpGBozE6/du3cPgtiGwhXp378/zJo1S3oOuhNQ3M/qAeua6rt6oPquLqi+qwuq7+qgJuH4rk3XgYzYfvvtg/1vJ0+e3PTZpEmTysQsQRAEQRAEQXghXgcMGADf//73YcSIETBmzBh455134KGHHoIjjjgiqywRabF8OcBf/wowdiyVedJgOV9/PYDijQZBEARB5IlMNym45ZZbgjivKFhxodZRRx0FxxxzTJZZItLgoYcA/v3v0u9PPw3QsSOVexLU1wM8/HDz7/fcQ+VMEARB5J5MxWttbS3cdNNNWWaByIJQuCJTpwJsvTXVQxKgYGXLmSAIgiAqgMzcBggioFikgiAIgiCINFm3DuDyywHOPhtg3rzclT2JVyJbSLwSBEHo95foCnT11QDz51OpEfa88grAyJG4Uh7g7rshb5B4JQhfBiWMjYezYZdpEkSeefZZgHPOAfj666xz4geffFJahIn/51BwEB4xfXrz7zl8vki8EtlCAqvEffcBnHsuwI03ZlwhBOEJS5cCPPYYwMSJABddlHVu/IAN/v7pp1nmhMg7xXwbN0i8EtmS8wfIGf/8Z+n/Dz/MOicE4QfLlmWdA/9Iad94osooFCBvkHglCIIgiDxA4pUgAki8EkSlQlZtIs9Q+20JiVeCCCDxSmQLdcbJQYM/kWeo/baE+st8tmPcTXL16qxzUlGQeCUIgiD8g8RrS0i85o8nnwS44AKAiy+mNu0QEq8EQRCEf5B4VYvXGhq+cwFugR5GiqBFiM6g1k9kCw1QBEGIICtjS9avb/6dxCtRxZB4JbKFxCtBEFFCzRVvvAHwpz8BLF+ezzInyytBBLQu/UcQGUHilSCINMQrboN5552l35csATjrrPyVO4lXggggyytBVCo0MSDyjGu3gS+/bP799dch94K+Vassc0IQmULilcgWElgEQaQhXn3fRUjnfsnyShBVIF4XLwZ4912AVauyzgkhg8RrcmVAZduSBQsA/v1vgKVL6ZmsNrcBnxc4vfkmwOGHAzz8sPo4Eq9EEhQ8n9gJ8PhpdsDllwPceGOznxNBVNOqahKvLbnsMoDbbwe4+eZkypzwV7z6PEDfcUfJyPLXv6qPI/FKEFUgXidPLv2P1leCqDbxSqGGWjJ9eun/zz9PpswJd1ST5VV30kniNd+QQcEZOXmaiYqFHuaWg7SrMiHxSuSZarK8VrN4XbsWqgYa75xRIa2fyC30MLcUma7KhMqWyDPVtmBL574rTbz+618A//d/APffn3VO8sOoUSW3p6+/bv5s9GiAu+4qhYOrkvGhAlo/kSpvvw1w7bVmD0kFP0BOIMsrQUQ/F3HJi9irJsvrPfcArFsH8I9/0JsiXS65BOC//wW46KLmz0aMKIV/O/tsqBZokwJCH+xkbrml9PunnwK8+GL80qNX28lZXqlsiTxDbgOVL15ZVq4E6NgRKhpfjTWFnLyVYKiw1k8kypo15ULWBSSwSLwSRBp9Q17Enuq+WUGfl/sxEa9JiEV23MoaX8VrDqmw1k/k7sGjhzk5twEqW6IaLK8ffABw/PEAjz2mPi4vYk/X57XSdthascJ9+znnHIBjjgEYOxa8gPpkZ+TkaSYq8jUeQpZXCpWVFjRwVGZ/c/31APPnAzz7bDJ9VNpUk89rkuIV12fg2gxMF2O+ExVFhbV+Bhqo3EPiNRlowRZB2PU3vPuSKuxSXsReNUUbSFK8srvoJeGSkHddUvQoLxZUWOuPQTXFmvNJvOb8AXICLdjKppyJ/NdXuOmETj+el0UpZHl1g4/insY7Z3hYuxk0kttuK+0r/Z//mF8HO8vXXgP48kuoeMjymgwU5zUdaOCovP5m6tTyv0m85hfX1lEfJyvUBzmDxOuyZQBvvVXq9G691bwEX3oJ4O67AS69FGDWLKhoXEUYYKGHmdwG0oIsr+bMnFla7JLFc6ojXletqrw3aOQ24AayvFY0JF7jCrI//7n59/feg4omicGfBAW5DaQFTZTMmDcP4LTTAC64AOCLL8BL8cr3HyrxmlX9f/UVwNNPAyxerHc8uQ0QRBWL1zy/WqgmyyuJV7K8JgGKmDlz1KKAxKya555rfj6vuw5ysWBLFdMzi/pGy/DFFwM88QTAvffqnUOWVzdUi+W1WJ3rRjysXUeQKHIPLdhKBlqw5b48zzwT4Fe/KoXLkZUz9RH6E/ssXsfrPBd8n6TKZxb1/e235fFodajWTQpcRxvw0TDlq3gteFhWEVRY6ycShRZsJUNSg2q1zcjRyjV3bmnxJPpqIuF2xqLycFXu6BN6++0Ao0dDRdGa2T08i/ipOiHkTMRrFs+DzTV13QYqbZMC123MR3Eftw2Kzi9WWT/fCNM7pc/rr78OZ5xxRtln+++/P9x1113xE6/SCk0UEq/plKsrUVVplkV8RYw/7dqJVyqfdFIptuPBB4vPT0q8ok8o8u9/A7z8MlQMWYsjkaWcFyT8s6NybcpiTMAg+aboug0Q+ROvcUlKvBbzp5cyFa8TJkyAvffeG6655pqmz9q2besmcZ8qI8xLDk3zZVC0gWRIakCqpIEOXymeemrpf7Sm9u1b/v3f/w6wZEnpd5mATEq8VipZi1edSZ3PbgPoc40W+aQsr9R+1fg43sbVJaI6b6jOfizTqcnEiRNh8ODB0L1796afuro6qCgwfBZahM4/X72YIA9QtIF0ytXVxCuPndry5QAPPgjw/PPl5fDiiwALFpRcA264QXxeFOTzWt3iNW1eeMHuPF3x6pOBRsUbb0BXnFBGPaOu76cSxWte6rwaxGu/fv2SSdyXSsYNEHAGPn48wN/+BrmGog3ky23Al2fAhD/9CeCVVwAefRTg88/FizlCf1ZT0UKWVzM22AAyRcfn1SRUVtqTOVH52YT/kn2Xh+d70iSouftu6PLvf0Phj39UH+v6fnx0G0hCvDZkaKRYuDCzONCZuQ0Ui0WYPHkyvPvuu/Dggw/C+vXr4YADDoCzzjoL2rRpIzynoaEhOE6LdeughinQBtl569frHSeh7FxsRNz5NRMnNlVscdYsKGax8MEVa9bEKquQFmlw6YR1rF3XeWft2sgysaGwbh0UHNRX0rD1XfPqq02fF0eOhOL22we/F1q3Vt5LYfXqsu9Zmo7lyxmFjoMycfFMeEmhkMi96T7fhbVry+sc64sbGwpr1pQfg5Z5Wbq6Y4IjCh06tGiTDbgpTqdOLY4tyxcaCSR5K6xf35QmjiXejyeffdZcBq+9Buu5NS7sfRcbGtzeT0ODF8+mbt3GGStqLO6z0NDQ3JYa25MR2OeeeioUli2DhrPOAvjBD0rXT0lMZyZe6+vrYeXKlYFQveOOO2DGjBlw7bXXwqpVq+Cyyy6T+sjqUrNyJfQPfeDQyjtypPi4Zcu0jpMxkDl3/uTJsIg7f8CiRYGIQJbMmAFzDdP3iY5jxsDGMcpKVGYLJk2ChZJ0Ro0aVfZ3zfLl0NChg5+vgxyW65Qvv4T1XbrETrfrlCnQxUF9pQXWN9s2Fs6YAQsa87zhrFmwkeJeuk+ZAnXM9yzhsa0WLYJ+bDmPHAnrO3eOne+BOSpjEzpPmwbdErw3/vnm6T5tWlmdTh45Eho6diw7ptuUKdCZOWbu2LGwpFs3YXp1EyZA9xTrqm7u3LLrIVM//hjWbbSRsg1NHz0a1qCLjIB+CxZAK1yUiGsUv/0W6j1vb10mTYKujfldsnRpizJn73vpzJkwx+H9dJw0ycl4FRf2Hqd99RWsnT/fOq3CqlUwgLunguAzHXrMnAm1jefhxBCfLxNqP/kEetTXl/649lqY2LUrpElm4rV3797w0UcfQefOnaFQKMBWW20VKPbf/OY3MGLECGgl8LcaNGgQdBLMWoUsWwY1jP/s0KFDxcctWqR3nAT23Nr+/TGB8u9RhDT6utb16QO9RennZUHX0qXlZYVWMYs8s2nUbbYZ9OXKBC0yOLANGTKkuR28/z7U3HILFAcMgCIu2PG9rExYsaKsTIZsuy2AZAA2oYBWjxhtOy3Y+t6AbRv9+8NmYZ6nTIGa99+X3kvhrbfK7pWl6dh588rLeZttnJRznP4jc0aPDl7nFvfaC+CnPy3/bvr0RO5N+HwLKLzzTlmdbofPBTfZKHz4YdkxtZtt1qIPbmLmzHTrauFCqHn99bKPth0wAEDgKlc2jmy1FQCOJQIKG20UjJdIXffu0MP39jZhAhRqawPhWldb26LMy8aCXr2gl8v7Wb7ci2eTzcM2W28NgG3U1T1tv32L8UP3Pgtvvw0FfOWPEqRTJ/PymT9feN1ly5YZGRpzGW2gC2ddGjhwIKxevRoWL14MGwlmpzU1NcrOjju4TOBIz9M9TgZ7LqbFn8+kXxB9jz4jl1xSeh2GC1HQsugrKLLZ+7Vd1MGkEXTEkjSwLprq46abSsejG8a4cQDYCVQSUe3IRbpZL8CJAPMXDsxIAX0Gwzy3b6++F3zlJZnQNB2L37NpKNpepZZxC8K3XPhcHXQQhnspj/Oa4L2VPd8youqL75PwLZcsTb7+k64rrvyCa6IhQ3Rd3eefuQdV3+kNhQIUmfzW8PnVHAusSLj9ZtK3i9pwjaWGYbWJTdnj8YLrok5Lg8w8mt955x3YZZddAteBkG+++SYQtCLhmhimjsa4P/Wbb+rvUx3FAw8AzJhRigf45JPgNUms1rZJI+9RG9KKNpCHBR0q2E4wKoSezYIt3/0F02b1ar+iVdhsUpDEolJbRPnV2UVKt13m4fk2aUOu7ydrYY/1yEdYcL1gq6HBPs2cRz7IzPI6bNiwIKYr+reefvrpMH36dLjpppvgJAwr5QLdgjWtAIxJi6b2zTcvRRKIm6epU9WrqH2CHxhcDP5ZD5A+kFQZ5F2csYNP1Mp3HdFCobLM+qeovhHbF2f1SSUKBxoO7ryz5ELAT1pMog1wVlvniMqPnyDonhfnOCJ9sD86/fTAVSnxZ7RYndvDZiZe0Xf1kUcegeuvvx4OPfRQ6NixI/zyl7/0X7w2+ogEoa90GkBUo8hTB5RESCff7h8XFzz+OK4oBDjiiHTcE9IKlZX0YJ2k5TWqnfgSKitvZWwr7nGije5O6PN2660togAk+lzcfz/AJ5+IzzHZHjYL8arT31WSeM3S8poluNteuJjJZZ8jmoAVK6jc8uLzuvnmm8Mfo2K/Ef7ADyaVaHl96SWAf/yj2UrS6GtbEZsU5E1YsZbXvIjXPMOXh+r5vuee0qYR+IMbsPTqBfB//4cLF5IXr++9Jz/HRLyGlmPfxKuqXbLn50G0VOszprNpik+W1xziYRTflB+aNCs+6lq+C4skLIS+dW5sGBNcxJIGSb3ez/tr8jyK17yVse3zPW1a8+9TpgTRQODcc5PLW1R+bNwGknarsQ0or2rreROvJnl0fT9Zlo+snl1vD9vQkF2fk3H7q1zx6surF1aQ5qGzSdvn1bcyYQe/tHZoScvymjdhZSJedXxe03IbqAbx2q5d8vdu81yYWl6TJAnLq2laWZNlHn0sH5s84Q6dN99c2nmQLK9+uA0QOSMJy4VvgioL8aqzqjqNBTh58nnlXSB0LK9pWN58a88m8BMAkdgP60QkXl1jM9nQfeWum55vllcftofFRcbdu+uFdaxWy6vs2jZ5uvHGUnjI//635J7Dp1fMKNpAxpDlNesKzPr6cQa3Slywxd5jpVleXaSL5YOLZWLsEuPE8sqv2ia3geQtr2wZZyFedSYbvotXlwu2spgoYZhI3OIVV9LbRPjwbSxIyjfVJShcQ6ZPL/+uSD6vlYePbgN5x5cFW0mWaSVZXpNwG3j+eYCrrwY477zkX7uqLK824jUNNwrfJmMuxSsbXzkL8RpX+OXF51V1TNaW/TvuKP2PIaBGj86n9TPk5ZdLEWVMQ17GubZPcV5zTuVaXnXxacGW71TDgq1K9nl1ke4TT5T+x1XmbIzitMUrb/Uhy2vyb1aytry6XOxUCT6vWfeduIOVz3mMKuuHHiod89Zb+RWvRbK8Vh6+CEXXVsJVq8A5uMvZc8+VVgynbXn1pZ5EA3RaVvM8WV7ThC3/KKuZL5sU+NaeTTCxTPoqXk2slr5aXk2sx1mCWzab3EvaMc+z9LdNi6JAvOrci0870VWl5RV3WokbkiJPjRatXhhL8bHH3KaL6f35zwA33KDe5asaLK/sQ51W3tLaHjZvCyJU5c+3RbK8Jm95Zf9OQ7zaTDYq3fLKnp9136mz/Sot2Conbp3pWF6LEW3s9dcBDj88vsWZQmVZgoGqjzmm5Htn20nkLc7r00+X0nj2WafZCnx/Qj7/3H+f1yRhRVCexSumYZIufmdan2kueIkSMrK8qELVJSFe8jQZNn0DwP7dti2kDkUb8KutmQrxKDesLO/HdV+W1L2I3kY0GI4fd91V7r+e07U5+bW8/u53pUrCQPJff53Ng5B155EEqg7Gl+1h01qwlbRlRnadOO0K6+TyywFOPrk8kHz4nez6F14IcOyxABMmQC7Q7bBZ61ASbgNpr2BPkqjnm73XNPo+m0ldJVhefQ+VZXJ9E7cB1/joNuDa51VE1u0iJfIlXvnBWPXKMA3LayWGilJ1MElsUuDbYM/eY1ri1WUZ/Oc/ACNHAnz7LcCsWXpt7Z13AMaMAViyBODKK/1puyqxpFs3qkVf5PNq7zaQxXPr2m3A1zivedmkIK7lNc3+w+WxUWAUBnxLmvR1bC2vJuDuebhzHi5u84xcidcaDNOBA6xL4lS0q9e7PpnuVXlIwnKVR7cBjA346actX734YHnF3VhkyO5n0aJyP3JdTPKJE89//MMsrqLKbcCFeGXTWLrU7N5l+fKtPadteXU5cNqUrcmCLV/Fq6+WV95I5NryWgkLtnBB9YgR8v4pjXssWlxDds5ll5XexqFrIf9WjnxeDVH5ZKZdsDodUVqCFCMGuLjntC2vWVsPbNwG0GUFLZR33unmmi4HVZPJh845LuoOY0L+5jcA998P8MgjLb8fPx464XMdtcOTTdtTbXQQlsfcuQDHH1/6mTFDL92oNPNIlE+7juXV5f3r+PdFnROVXpLoCgvbfKXdd65YIb8+TprffrtlNBz2Xvh+JukFpSa4urbKeOACUZk1JGh5ZSf0aWxMU6mW10TcAVy6DUSlldTD+e67AEceafbKV4ZKyIgGty+/BDjxRIAHHxS/jsadWHBXFhm+W175OkPLIb6WR3C7PpEV84orSoGvde/NpXuCqv5k7c+FeMXf//c/gLFjWx73l780D3z8CtdFi6DmggugJ0a7QMtskpZX2STh4YdLVnQUz7iYwQSfBuA0La8q/2lXiERd1NsOE8trFq+tXfrt+iJesVzPOQfgllsAfv976elF3m0gi8mDi2NVJG1dFvU3xQrqgwwg8erSbUCUlmq1sytw72MceD/7TB3qyrV4xc7n0ktLlqtXXin5WbLcdFPpdXG4K4vv4hXrJ8r6N368Og20LqIVEUUahiTRIa1BNUnLK05UrroK4IILSr62aG0VWSPC8sUZPfqFMb5UNSKrrAyXPq/soLxsmX4eVGlWerQBWRsNJ38u4kiKrs/vrBZ1TpZ1JUpf9JmJ+MjSbQDf7rGE10eXG/xBXnvNH8urSXppPbeuxWuDwMBicw0fXBYrXrzqEtWJ4etB7KzjzM5ciA7XD2xU5+4y2gD/t87A7/MsUdQe+HuMEq9ofQzBSBi613VVJmlaXtn2f/PNzb8PHw5wwgmlrWQR1loWti+00mDcYhS9OunbulbouA3E6bh9bs+uLa933928JajK8or+4Ljt5rXXxsuP6FnEzVRMzolKL210rLG6+fJFvOq6KkUt2KoEt4Gkr5Ol5bUQMflImfyJVxcF9swzAKeeWnq9GwfXCwBczH7ilk8cn1ebV2I+WapE1iL+Hk0WDOrWhUvLa5riNSwb2avcRx9t+X3otxi6Xqhw4fOq4zYQB5/bs+toA998U1qMoqobTAPdl9D38aOPAKZPt88Pfw2c6Pz1r24C/Ecd66vlNUu3AVk/pdvnpC1+bMs1SZKwLhc1/GArkMoVr6rjwv3Z0V/TZLtVG8urz7HtXEcbyLt4FYVcS2rr1kpwGwjLZtQo9XG8uDWJOJBktAGR5RWPefxxgNtvb34VqsIn61FcdOM4qwZHPo04b4L4ssQwcKbnsPB5/vhjuwgTcYgjXrNua7LnR/W2zifxim0TXZpwkx/dMIIsOBnDN0yTJ0Nm6AjVV14p7aD11FNur+2Za0FrqFR0HwyTcEdRIsMHIZam5dWFeLXJb1LlLBKvaQjLrBds2RKWDesqIYIXMLrWa5WVSbfM0JUFtz/eemu5AGbLDN2J8M1MOCiffbY6/WoKlcX2mSrLK0vrGEOMTVmaWF7RiottF/3UkxiYdVwERJ+pJg1RaclAv3NcmDh4MMBhh4ERuG4BhRC/q1p4fV2LZtZuA1iuGPUEJyx//7vZtfH70A3m/fcBXnzRH7cBHpx8I08+WXLfqVDyZ3l13UhMLANxfV7zYJUxjTaQhXhN020gDfGadaisqG0bZYT1HxVeip8g6lg0TcSrqk5wMRZaWvBVNn/dqHJ+4w2zPEblxXd0Q+Hh2ypZ2bkMp+d6DYHoO1zg6iJe8yeflNzQcNGs6npZWV5x8ewHH5RceSZNAiMuvrgU0YaPGqMjXn1asIV5CS3tCxaYXZv9Ps5iRNfitUEjfFwOraqVKV5dN/CkxatJo/ChAcURrzYBuG0GtzQtr0lsieub24AtUT6vSYpXnZijPHw0jKjzNtigpWi7916AP/xBPnDn2fKq2w5x4Y6u24DomcpKvOpaNG24+upSlJHf/tZPyysbwk62U6UM2fMquz4bXcTEbSAN8ap7LB833VXs3bj3KHpGiylNmH3QJxUnXkWFqtvYXIrXOAG0ZWCaGFIIX31m7TbAD0w2M1AXg31SD2tWPq9Zuw3EjTagEih4Tf4Z0/V51bW86rahdu3Mog3w4hXdCf75z9Irw9Aqi/5lojTzSJRbECvidd0G0havJgJF53O0OEa5xZjkxebtlOxz27YmSh/bNIZcNAmyL5vA4at50bWyjvOq6mfZvOBamKOPLrkMmeYx7UWAxRTFq2fkz+cVw7XI9g22Ic6CrTReGeKrmnCBmW8+ry4srz49eDo+r2lYXpPqAF0v2NIRr/gdL15169y1eOUxFa/sgiEMB7XDDiW/Mlme84bu86yyvGYtXk3dBhDZvbz3XknQhaHdttjCPD8619K1eLtqW3w6GBEC3yYgs2eXFiuapMOnx76SV7kNpI3uxAZjlyO4MAv9W/fcU7/sv/5a/X3cOhQZU4oZ9TkZ93X5tLzqzA6z8HlNwooo2sUpTnpRqPwf+fsz2c5TdozrRRlJh8qKukebDtrl/dhYXm3bTFg2KoGCQkfXomeSLxvxajrZatOm5DOLfoL4w/vJsRsxmObFB6IsUyrLq+o7lrT9A20srzJwgVPISy/Jj+PvEf+W+SLGcRtgYymrjrMRryH8/vU66ehOGPi+KWl3KT69qDdEqti2OnmbOrW0+M8kT6aINEcxxxPmqrK86qJboSaW17ihopIQPjbgIIUO/DgYozN+3uK8on/ZKacAHHwwOEX0WikLt4E417ARr7ZuC+F5KoEi2rjCRryqQmXppmf6vKJ4xS1rw80W+HLW9Wn0lajJqMryKrtP3j/SJ8ur7H50fExVzxXv0/2rX5WiLKBlXvdaUcfgwrIPPwQnYPpYh+3bx1uwqRKv+BmWGXsvUddJWryqDFUufFWjNtDQTcdXy2uBfF7TIQ3Lq84r8LgVbnq+zn2jNRdfh4wZA3Dbbe58XtMSr8iDD9qdp0KUl2qI8xr3tbtKoIjCYrkQJTYLtmSTT1mbR/GKWx+bkCfLa9Rk3Mbn1WfxavrmQfeVNy9e0SiAbwdffVWdpuwzUX7YLYxVaelw110ARx7ZLIbjbg8tym+YN7b/jLK8uobPl8pQFbUewFWfzL+9qSSf11Wrgq3gC6q3FFXtNuAavkH/5S+lYMRpuQ1EIZut2jZYDM5cX1/unO5KvNr4vMYpM9cPrcgCGaeOdfOXhkBWpWtbBzpuAyLxaiM2Xfi8ykSCSrzyW2KKzlddw2f4/Jss2JJZ63nxGsdtwLVLkWm0AV3LaxwDiIg02hDWy3XXubG8ihCJ17TjvLqwvOqKVwyV9vbb0XnCaCUPPQS5FK/FYsnXP9ykgb8uutm8+SYUTF0dLSG3Ab5B//GPpf8feACgd28z4eXaUhB1nqhDVaWHrzX+/OfkQmW58LMzTcvlq4ys3AZcWl5t4rzGdRtQiVdRmB0XvuFsnnXzb1qXuGBLZPFSne+LFUQHlSuG6Pu03QZscNkH677yjitebUWcK5GbleU1XAyXFCbiNaosRS42uGkBWlKHDy+5suny8ssA3/0uwFZbgTGiZ7QhpQkzLlgNBfqf/tTy+3/9C9Kkci2vcd0G2GDTsjR1rHJJuQ3YLL6JEq6q68W1vMo6ujgPnuuHVsfyajI46ta9C4EcZUVMwm1Ax+c1jtsAe1wSlledBVsqy6sI07rD4zGiCIbgSttqayteVW4DvI9zHPGaluVV9Izj6ntWiLuyvOos4nKxoNGEuNFGVIKcfbPJXgfLjF+Z7+p+Fi5sGdOZz0tcUKiFVsgRI8zPv/BC/XjXvrzZeZuxLGN/lTGVIV5NLZBJ+7yin9xXX7W0DLq0yuiKV3xgUYjr3mfabgNxLK9piNc4mxTYWlBM2wnGAMa9rF97TX2cTGTG8XnFvKrEaxzLa9LiNUrwo+XVVLyaluXHHwd+YsHrRIwpmiZRz6JNqCze/zPtaAMufF6xHvjtP5N0G7C1vLoaT2RWZRSB6B+rAicrogkKvrHAnffY6AVJjYUsaAk98USAk08u35ghSrxeconZxIeN0GCLSWQHWT4aUrS8evaGqTLcBnQ6BGxsm26afLQBfGhPP73UwZ97bnQ+k7a83nlnqTPecUeIjSgkjInlVbWZBG5ZiC4bO+8McOCBevlRlSf6IKH/MkYk2GMPvfRsF1O4DpVl2k5w+9MwBjK+wkpLvKLYibKsJWV5TWPBFg7qSfu8sosbXngBYK+9wHvLaxgKSke8StpHYe1aKFx/fSkdDGrfoYM/4nX06JbHqNwGTLaXddHHmB4Xhaz949iB1kUZaDnFV/+iOsaIC0OHmuXXxf08/nhzP8f7oKomGRj2Do/fe2+9snEhGG3uN6loA4UEFoYnTGVYXnUK8rTTxIO3bJ9nHWuuaDBEv5BwwNMN9myDiRUB0d0lRlWWcbaHDb+TCQgcwEaOLDmz675OwXPxWNy3/ne/K69f9EHCQQg/TzLaQNyHGDtU073GVXk2CXUW1/qN146yrMUNlYVtAhcVJml5VR1vKl5N2wM78PObIiRNVHmoFjjJ7lNzk4IN33gDCrgwFvsl2SYsNs+WiwVbGObKhLhuAzrHxHm7GIVMmKuEK4LGAdnkFdsBP+akIV5VE4moScb8+fYTPZdgFCBc/ITh0aLy0RAR5zVPC0gNqR7LK7JoEUC3bnZbVYrSjLvyD/fCZgOdizoo05kf+3mcV3Y8onuLivsqS0eUX7ZjEQkeEXgu+qaFnezgwQCHHKJ3LoayQV/DXXYBGDBAnn8d8RrHrxktBTwmbcjE6iMbaJK0vIq+172/L75o9tPedtuW1zZNT1aXsvrD45O2vGYpXqMmYqo+RvYdX9+SPqjdxIlqS6coPzqw5+Cz8dZbAL16AQwZoj/hF4lXV5ZXH90G0orfyeZXx/c37jXiTDKytLyG2+uiIYwfH0RjVFFxz+hO9utfA/zgB+AUD0SxN5bX4cOHw8V8wHzX4jXuw6HT2etawqZMAbjiimTN+YsXm53Lnx/10Ji4Dcic++MscgstcyHjxqmvzYLWWtzaE/evjuPzqrNITYUoJl4c8ao6Nwmf1yjxqrMITiVeZQLHpduADBvxmqblFSdgH3xgvyjK1vJqIl5lbgNsu2jVKhnLK7ph3HNPyZ8RjQSy9Ph7EYlXWrAVH5UbUBqYTDKi2sgbb0DiiMZvU8vrqlLs1UTI2HXAC/H697//Hd7WiZHmm3iN46v4+uvx8qJjeUWHe1N0Xweaug2E6eoMkLpxB3VFi2jGLXK4t4k2kMQDHEe8qspEJnLihMqKEk6iduNi1m4TKsu07lQhoWTXNb03tv4wukEUTz0FgL6iuJr6ggtKv4c+z6ZE5TVBt4Gy813GsmbPwckpuzAuKcur6wVbaUcbSEuE2Ew4XWJjeU27jHRcmfjjiwZ5030GbNOpJreBRYsWwU033QRD8LWOLbpC1eSBsfF5jWpI7He61kXTBVtxxauMuJZXXauHqeVVB3YrRBWuLa+2mHQMfIesOlfXbUDXFQLLJso1xfQ1ly42Pq+yupTda5RLkc7iG5eWV/SFCwUZWqJD3/CnnwY46ijw0vIqaR9llleZj6nNsyVrW1i2cXxeRW0Eyx/rRNfNSZY/nUmVjc8rvolC30l0jZKB/saDBkEqJO3z+te/lnaPTMNtwAWi+42aiJtaXpOCxCsuVrwRfvrTn8IcfAXmsiBVn8n2JXfR2fOfmTR+E59XnRkU+vi6QkfYqfwQw791xB9/z7IFXLqDW1REiVCspWl5RYubzK0ja7cBLAfRAC5q7z5YXuMu2HIpXpN0Gwh3tkFsYkTyRIlVlXiNaXktsGkn5TbAglZtXauTqB74NoLnXHaZ+WJLnTbD/y0TNFHlg25p2Ia/+UZ+DE58+vcHLyyvcQQR+lCjUE/K8pqWkcLGGNBgkDfZsTbuiRn7vWZqef3ggw/gf//7H7z88stwJfogRlAsFoMfngascIGQquGObcCONPw84kFpwPQEcU3ZNDE99u/i+vXBT9lnxSIUGv8Ov8cOpUbg61hsaCh9z1BgzpfdM5+npnwvWNCiDKIQliWyZk2LtIpr1pTlrezaDQ3l+QrT5ctHUGbrG9MI/h83DmrQZ01UBhhuR1S+fJlgB87dU4sywwGUq9/gOxSHzLloMWpxz8xgx7YtUX0GrFsHhVNPhQLGA9ZtzzJWrmxRzrI6L65eLcxPAfPDpoGdvMTSU1Z3a9dCcdUqZRvDY/j2i3lgPwufadGzLU133brme+GeQ9U5ZddtbC8opETPWHHlSuHnYV6x/bdoLyZ1h+0Fyzq8/9atxe2lKfHyZ6rsKxvXD77/4tor39abjmt8joXfcX2CqM0FzzWWfXjfYZocsnpRwd5D2XNRU9OinQv7rXnzoPDlly3bLNY356pSwy46s2m3Idwz23QM9m2XXhq4OBWPPFJY90F7kwiPGl2LMCfAw7ZkOnbIYJ/vsG8X9bXS/lKHUaMi84t9lao9NaAY48ai8DORnoiLUGOsXl1+fX685PvqtYymiboWIhkfWrTvRqTpYqQJCSb9eC7F6+rVq+G3v/0tXHHFFdCuXTutc5avWAHrBHEjZ40bBys4f7H248ZBL+7YaaNHw9q5c2GgxkKmuZMnwxJ2MRBOyOfMgc2YNGeNGQObMH8vmzEDVrdqBV2Zz9bV1EDrxr+Xz5oFs0eOhO5PPw11gvtYVl8P3zLXbL1gAfSVbLk2+YsvoKFTp+D3gUxas8eOheWN5bnhxImwkSjOpoLZ48Y1nc/SatEi6MeltbK+Htozn82dMAGW9OhR+mP9+rJ8Bflt3x5az5sHfZnPURxOHjmy7NgpaCno0gVGjRoF/a64AlpJrExTRo2CPgsWtChfvkxmjhwJq7i8s99PxKgPrVtDlwkTyuoO+Xb8eFjWuXPT35vMmgUdBPcVMmDRIig0WkOXzJgBc7k2hHT6/HPoqQhQPXv8eFjesaP0+3YTJkCPp5+GFdtsA8u33basnc+fNKnFPYQsmDgRFgry033q1LL2OPnzz8vuqYl168rKbfHUqbBs1CjorWhja2bPhjbc9wsnT4YNBecsMbAmLpk2rals206ZAn002vmy6dOhE5f/eSNHQs/6+rLPQ4orVgSCR5hWfT0sGzsWNubby7hxsKzxudRhILOQaPHs2UF+ZNRNngzdJfc5UXGeDL4/Y58fpN+8edBK1E/NmAFt5sxpUa/IukKh6XkM0pw+vSzNkM0aGprqe8XcuTBLcEx/nHwbbhKBYg77E4Rtq7MmT4Yus2eX9Vch0776CtbOnx88t/2uvFJ4zbCtyJ4FXZbNnFnWxyN8n7ho6lSYP3IkdPziC9gYF+QhN90kr3eJeLXJX1OaMc6XsXD+fJjamDa2K348WTV3btBX27DhmDGRY93yqVOho+KY+ZMnwyJuLPp2woSg/w/GYsflEegWzsqPE47+7NjE5Wcl14ZnjR0b9FF8PySrU4yvPEBw7PqGBpgiKHubNrC6rg4qWrzec889sO2228J3cY9fTTp26AAdBQVTiz47fEBkbAjcsdtsvTUU3nwTChqFW4uvUvg0Z8woSxOvy/5dt8kmUOzXr+yz4kYbNb0iw+83HjoUan77WwBBHup69YJNmGsWzjtPmtft0Ee4UVSV5WngwKZ8F775RutedcoSF4jw5VmH98Zeu2/f5nNxhsd8tx2GOqqtBZg1qzyddu2gM5YJ89nWW24Jo2bPDvygN8COWXIPQ7A++fJFP1+cvcry1Qj7/dDttiu9VpwwocU91vbrV3ZuoWdPKDDx95ruK0y3S5cmN4W63r2ht6gs581rcR2tOgjzcNNNQSeEkRYahg0rb4N9+0rrHPPTV5Bu4b33ys7ZbpttxGWOlkb2Wr16QZF7BniC9s+9fq/bdNOy6+FMHYVMXW0tFDRfX9VtvHGpbNGa8vHHWu0cz2GPq+vTB/oMHQqF3r2hwL6SD1G4FwX3PmBAy/aCYdcUdcdTw7SdID/bby9/hTd3rrSshxpcU9af1fXsGfRPIYUuXaAgKUd8xgoCkVesqyuzDNV161aWJoLWt8UNDU31je2yJ3sMno/p4/UNIzAUO3VqKouyPmDLLUv9oWDb0G1wj3nsIz78EGrweoJrYn1jW2kC3zhYDNJYdmwfH1BfX14PffrApngM97mIoO+SuF3Y5C9IU1B+cQif7w27dIENw3vHiQk/nnTrBt1l7TjCDz+wlkfklx+veMK+vqzdbL556XmeM8dZeSj7eW5swLooaxtdu5be2LH54/plZf/AjY0hxS5doAubl//9DwqjRxvrB2S5aMORShKvGGFg3rx5MGzYsODvNY3Wqn/961/wOVrCBGBHJxrcWuFn/AOMf3PHtsIHAF/XawyQrfB8Pk1cccqcG1yX+bspb/xnjX8X8HxBvpqODb9nX+dIjm3FHsvmif/c0JdFeN+SziMYpPjyCM9FQcl+F6bL33+xWPqOT2fduuBzlZgJ7pUtX5w5fvhh872Hx4VuAWWZF+RbcI/BWapz+Tpj2khTfeIWi7jABleHYyxarh1ptWcW3MWt8fxWKLok7U24SEaWLntPWA6i47h8B/XP1bPwmny7keRR9nwL0w3ziAtO/vEPrXYeHMHmH/8RtcnyTIk/5tqetD1EZorJz6uvlraMxF2LRKIN/ZBl/UF4Td3FduG12evzeZekFdSR7DvuNXYB3SIE5RG4JDTWdyBQ2fxjnEts1zgmGPZfTXUa3h/fV4vGj/C+cQIrq2/RM2GxmCcou4hxpekYVbvk8y6+mHH+gjQF5ecCTK1GMTYJywa57bZSLG9sFzJxiwaDiPwGE/6o8YQfi9h6cFwe0n6evz7fZ/HjT416PGlKB5EcW1b2+Ebk2mtb5EUX3T48t+L1scceg3XMK7lbbrkl+P8CHOCzWsnHouNQrhMqKqmK1Fk563K1rigtVbQB2WIQjUUhhXPPhf5oQY2KT8efK1uYFbVgK6xH200KVMya1RwW7fzzAV5+GZzC78Jis2BLN+SZqO5sFmy58IkK82xSnlF1aYLuAj9Txo8HQBH7k5+0/C6qbjEc1H33Aey3H8Bxx0VfK+pZdLFgS+J3KV2wNWpUyz3pTbDJc3iOysrrKsqIznkm7dKDVd/aqMYHGTNmlDaaQC6/XP6844Q+CtXCNRnh+J1EOeuMqTqao1h0u2BLFEbSQzITr7179y77u2Ojj19ffH1jim7Imrir8qI6MI2FYIljswpbB5EIUYlX2bX5jQRE4hUXtOFrOeysVNhGG9CNhCA6VtYGUGyLfB15cafTRkzaEHbuqvzZRBswiUpgszrWxdaKmBfTTThchzlzvWI6ZMEC8eeqiQK+wr/hhtLvzz9fCp0Vtc1pVF5V7UA3VJZs0RCeH8ZOZcWr7YYLISqBGjXhV5WXSUxr0/zFaZdZjzcmsHkV5Vv0ma7PpY54jUMS5cwsrIa//a20Cxy/A2hUuys6ivPKEvcZrJY4r6kSV7zyFS8Sb/wx3CtyZ4T5e/BBeR5tric7xzRUlii/WF74Coj/XIZs1XtUXnl43zxZh2AbKuuTTwCuuQYAF6tFlQmGhIvKd1Q7Rb/aMAwaL8xtNikQiVLd46LuJSnximkoFr0JkdWl7qYYUenJPjMVh7L8RIlX03zYxnlVDZj85/gKUuRyhPUnEq+ai3edWl7Dz1WWV5O+ziZ/omN0+rekxGsS6UaNTeEYwU4idLc4T0q8Jml5DdvQzTeX3MpwR7jrrhMfI/u7aCheddpfnJ3IqlG8/u53v0veAmLS4egEVzdtSEk8AK+8Ym79tEFkQTGxvOK9/+9/bstF9/74gV2Wbx0xwucVO82rry79LlgM0uJaeEzcdrDRRvIYvqq0dXfYkg0YInETVQc6FnsbsLNnFi7Echuwde2J6zYwbZqZeFUNKqI2HrVjV1TbtrG88mAamO+2beXXZsWrzi5jKsJ7MGmrsnPSsryauiLZHmtCEta3KJc2XONx/PEln+/wzaxuPgyjUkTmz+Q7W8IyYLe+RjczFv7+41peGxoqxvLqxfawiZCGz6trtwETK23aPq/hK0ndDl3UIfMuA7J0dDEZQFV/qyyvUQPLqae697OMOka1IMiF24Cu37OO5VVnW2EbMIJB3ADxcQekuG4DoggHtuKVtzzplLGqPMLFeCJU34kQhEAr22HL5RqB8B5Mtk3WsXQmaXlNym3AMFJD4gJGx+cVXYHQEmlazi76FBVJWl5Vn0W142ICbgM5sbxWhnjF17Y6AZlNLD46VhUdtwEWk0YWtchIhmqXKx1k54hey/AdXNQOW7riSRfd+8Pj8IEMZ7Um4jWun6TIAd+lz6vJuSbbwybpNuDC8mpDHAtXElvfTpliJl5V5WYjXlXlEeUCZHKfor5Zdi1XEwp+NyUda3HUPeseq0Ln/nzweU1CwOi6tLGLUHX7iqTEq45VPknxGmV5bTB4CyI6XzRpJMtrivznP6Wt8JJ2G4g6RmSRsLUk2GxNGeZB9HsSmFpeVf5oSVtezzgDYPjw0srVOOLVdJZr46MUp97StLzqiJikLK82pLFgSydNnJiir7pglz1ryyu/sEVn0FeJ1yhLZfg9vgXAWJOG4rXM8upykWl4D3zEAlVbVbkNifIYdWxca72JYFIdE6csk3Yb0O3jKlm8ohsd/2xEiVf++1GjSqEYdSHLq4dgiJmohyOueI2afYcO5y4exihLss6CCZcPNBusONxFyyRUVpbi9bPPmq2uuGBMNpuNqnMUD6adWB4sr7quDSK3GZu8ZmV5dRkqK454ffLJlr7qccXr/PnxJ9+6dcQ+y9gX4PbNKqJ2TnPppx+6J4V+6En5vLq0vEatpTBNz8UbrQQsr4UkxauLPuWxx0oLatMSrxjajje6mboNvPaa2TV17iMnlldvFmw5R1RJcd0GojoZkaWNt7yqHMvZY6Mc0GWNMCm3gfCh6dWrefEFf68Youd73wPAnUpMFkwk7TbAD54yEaOyitx5Z8nCb9pJuvazjEKVvq7lVbeebMWrb5ZX2/Zn6zYg2fI5VrQBDNXm0vKq6zagsymDbpxlUZ5sEFmiVC4uJpZX/B/7AX6y4KN4tS3Lv/4V4NlnwTk2kxTd/tbV4mTs50XpJtVvo9EtjtuAKTrn58TntbLEK1Z8GGYjruVV5/WOyKoWdQ3dVZE61rmojtCleA3vCweqsIxFHcuZZ5YCSYtEW1aWV37wNF2whQ/zG2/Y5dHG8hqng3LhNqCbtuligah8JI1sImE7MNlaXqP6CBvxyseG1cmHamKlK16jdoPTyYvr2NSi9qUz0dIRr2jpwo0gbNGZ8JhMqpIQVY88AomQB7cB2aLipI0Otm4DplSQ5bUyFmyFuHQbQKfx//63fBYSNUOOchvA83Ubho54jRo840YbwN/DcERhWjiwqiYIonRU+UVw6z8bdDsUfiZp4jagsyhJRRI+r7Lvw207ZfD3jWG70I3iyy+r0/KqsrhHgWWtM8F1KV5VFpE0La/sJF1HvEbdr84qdBNE5RfX5zX87ve/j5c3E8urzaI7n9EVr+F32N7vvz/bCXHSllceU7cBU8jy6ilxVpOKnKnx57DDmrdajIo2oOM2oJsHW/Gqu+hCBwyP9cEHAEceWS5edfZuNxGv/MYFutjen6nlNU7HZWN5jUJ2PobGMbG8opXcJD6iK/Hqi89r3IHJNs5rGpZXnTJ2YXkN936P42vsQrzihDq8Z9EiWdc+r76LVx9FranlFXed0sWV24Do7Spa201jStuStNtAsXIsr5XlNhD1cJisymP9OGX7hIsGc5vtMm3ivOpYXuO4DeD/KFzDBSbhAMVaXlVpmIhXW1z5KiYpXpOI82rr58Z2ShgjVSZcowZ59u9KsLym5TaA18GFg1HXsxGvfLSBNN0GonYoi1rs5MJtACduYd9rannVscCHeYwrCnWs9eExNhOQrJ8xFbpvBcPykIWSY0F3Lpy4udphi68L9E/n/VKzFK9x298tt4g/xzLE3b0uvJDEaya4dvwPwdeq775bWogUZVWLeihdWV5llgRXbgN8nESRz6sqb2mIV3zYokKKiTDdYStOvpOINqB6/ak6l70PfieXarC8ul6wZeo2gMHX33knOl0b8cp/Jyrj6dNLg/zgwWIXE/Zv3WgDmE5UfxAlXl1ZXkNklleVeEU//Yceyp/lVYaP4tXU8tq+vfr7r75qucDKNWkKV9Fz7NJtYOlS9f18+GHJysyP/Z5SuZZXl1x6qb5VLY7llU8r6vuowdOmPF59FWDgQIDOncXf67gNiPxEkxCvNsIVMQ2VlRefV12/PVG+dNIQiT+bOs1KvIomEuz/SUcb0BGuiOz5Uvm88t/xeUP/5tNPL+UPJ33bbRfP8sq6EUXFsuavk5TlNcTUbQDLTiVc+TzGoZotr+G9Rvnm64rX99+HiiNKvMYZi+bP19v1LyfitbIWbLl2/I9CZEEwWZwgAy0kUeF0dNwGbAYCDO6NvpCy2Iw6bgMi8RrXgukK3DvddJMC27Ykmsy4EK8qcakrbKMW5Ohc16XoSwPZKziXllcXbdxmwVZUGT/xRPP94t7xuq/Jo9pZKF5VE1oTtwGdNiXaFIHtk2RuA7L71fHxS8qnsposr6YTxnbt1N/H3UrYR/i26NJtYL6GeK2tzY14rVzLaxrilfczi1qwFeU2gMei2Dn/fL04r0mIV1ZAx7G8iq7Nr4jOAoxRa+I2ENdinITPq8onVdfyquNzFvW57YTEF/EaZ8GWLNoABjmfNw+gWzdwjsmgEvXan/2f/9zU5xXBPkF2TpTl1dToIOp/2M9MLa8uFtm5tLyaiFdZ3rN6uxGFyYSRt7ymIVazXuj2738nt2Br3jw98errxKeixWtSPq8y0A/W5YKtsIHprADH+xMd50rAywZKHCTYV3S694hbYY4eDZnDLuzw1W0A8zdyJMCgQQCdOpmJS9eWV1ywiFvqHnRQy7A1th1pVgOrzIrh0oL85pulAQi3I/76a4Af/Qhgiy3M0pXlJ47llbVGyiZruuKV/S5MV7Voi0+Lr39TtwHRtaIsr6qQdzqDta8LtmT4KkBMyo+fpPBv/CrR8srj0ud1noZ4jYpY4xGVJV7TtryKrs93GvX16u95wt2rosCO7de/Ts7yKhOvOpZXDOK9aFHLQd0HMP8y30eZeE17wRZuU4hCq29fgLvvbtlJ61pX44pX7DgffbT0uyjeog8+rz//OcCLL+odK7JiLF5sH2dYZQHHegvbPS4GSlu88mWclHhlLa86x0f9rdNvi8RpHMuriXj1zW0gvKcvvij1uXvuqbaCZ43JIkn+GL6NVaN4tdU1q1bphfuSvTX1kMoSr2n7vIquH9VpRH2vm28ccEVp4bZ+W24JsNNOyYnXKJ9Xnx3pscxM3QbStryGImvqVIBly0qvcvg0dK8Xx23AVVvlcTmwyhYW6lwX82+7WxKW3z//qX+8SVnJjjVxG+DrVxSGTyVeVRMMkeXVlXjV6bNcuw2kYXnFcmLf4uDC2FGjAI45Jp54xXOnTQO47LJmkbL//v67DeiUI19faBXEPvGzzwC+//1kxGvWbgM8zz3nJn/HHBO9TTOC7ca3MqgK8Zp1oetEG1B1ziYPo2qRwdVXl6w9ccpDFjdPd5MCX8H64esonN2iFc6124DMypvUoiCXlledxYdZW16jViSrwPu0nWjh7nsm1zER7KI2gmVmkoaO5dXW59W1eDV9WxRleTVdsBXHt1QXzF+YB7SAhZMmXFtw/PHxxCs7ibr33pJ49dXyqiNeZf0P1iu64yAffVQy0hB66AjXsN/IieWVog24RGeQchGNQFfoJiFedeK8+oxIOF11FcDEiaXFNq7dBmwsr1HWNtX5KmGYhHi1aWMuB9Y44jWtThpDuplcS1SmJi4Dtj6vqvPjuA3w15G9+Uh6wZaNz2uYlgvLa3g+bj0egoH4r7yyZV5NtvgVWdV9Fa8miyRVbgMY47Ua3AbSZn1+3AZIvPrmNuCSOI1QFkNVx20gjw+nbIvarLeHFYkW3UVZccVr1OBpWzY+WV7TAF0/8iBeddsHW386C7aiog249nk1tbyq2mO4ODVuvx3mT2eMUC0gjRKvIb67DdhEeOBDZ5F4dQ9ZXqt0wRa+yo8SBUn5EboWr5XqNiBzt0CBISKO5VVU36b1KxItLtwGXPi8kuXVX/HKixdT61ze3AbYCbXM1UaWNvpRpiledSPSxHFnqAS3AdXkCqkGn9e0WZcfn1eyvLoWRlEBr12J16zcBvJueZVNINh7Ou64ZMWrSXqm4lW3fSUhXqP2uM+L5RW3geZf5aYpXtnXyiGmgcNVlledBVtJRhtQhcpy4TZgGjVEFQkl7Bfi+r6zbgO6b+d0nhM8VhYaLO8LtqImPUmQE+GWGGR5zYisfTV0xGtUHnUfnr//PRvLKw4Seba8IqI6Yu+JjWMbx21AtIDPheXVNtqA7itJW7eBLNqFbmg5EbJyRDEQJ12eN94wexb/9rdSnN844lW1YCvENtoA+53rUFkuFmyJ8m4rPnmLrm1fEObPZFGvjljDSdbbb7f8/MknoeJ8XvnyILcB95B4rdJQWSiKona2cWV5lYlLFpfxSdmBI2qTAt95/HH19+xAGGfAEvm3ma48D8UrK7hN3Qb4+0lqwVYW4jXONWXl6No1BsWFzIdcBi4idGl5lS1iYmHrUscyr+s2gHm/9lqASy4p7UqochtIyvJquxCF7eviWP7CctJxG8C28tprpUWkUSxYUPrh+fJLqDifVxKvybOO3AayezDwFdANN6j9mCrBbUAHl/FJXQ/sRx0FXsEO7rzYi+M2ILK8mroN3HknwP/9X/PrTVPLq8gf0IXllS+bLMRr794AW21lZ4WRlSOmpesCoYsoDJsKWTg3F24DcS2vonRVdf/OO6XQRhjX9JFHkncbEPXBbFvdfHOAU08FLaJ8aZNwG8C6xk0uli6FiiNOnFf+b9fPKFLtbgPrKdpANmD8vDvuKMVuxFdvPsZTiwrL4urhMRVJPLLBy1WorCQ6Hlf54f0D47gNxPV5xQEMXztjfWDbDtNIUrzqvrZMYzCJAp+ZG28sWdN33tnsXFk9YJquhbipeHW9YEvHL1LX55UlnDTo1j1urqJ6HezCbSDK59VkcuLK8hrmD9vBwoVQtcTxec3aLbAaWLs2N+XsmYKIySefZJ0DteU1T9EGVG4DLgZ238Srr5ZX0WTI1G1AJF518qDj85q15TWsu7o6c+uryvLq+l7wtbkM2WSQtbyZug3orMyOs/CPT1e3vETPhAu3AbYMRe2Wj5Cg21bYdF2IV8R0u+BKwqXPq6/hwPLsk7ueLK/ZMHo0eE9aq0BR9LiMT1pNlldXPq8yy2vcTQpMRQdrPUrS8pq2eOUHB9PBQlaOOAHl76VNG0gMfvvfkFtvdbdgy9TnNUnxmsSCLfYzUd75RWa6/Y8L8YrX23XX5r8XLYJE8fnVdxyf16jNLnwjjwub15HPazbMmwdeo2N1c2WyRx/JceOSsbxWonhlB3c+oHuWllfTHbZ03QZsBg9eDGft88q3Q9M2JStHDG3Fp+Uy+sCgQXriFV+xu1qwJaqbtHxeo/LGTuhc+LxGWV5N8urCbQD7lj33hNTwWdS5jPOaxH36LPzTYB1tD5sNph18Fg3jscfy8Zoh6QVbvolXFvb+4vi8ioSvqXhduVKcruqarnxeVdZ3hC+XrMWrK7cBXO2dpOWVj9aBLg9RRC0EdW15NZ2wmdS9LJSV7nX5voO3pIrSt7W8it5amJKEG4oKX2O8xl2wxd+X6TORNj6N5xUoXnMcbT6H4H7MeZn5JW159e2VCvvAJrlgSxQ+y9TnVZUfF+I1RHYMu1uQrNzyIF5lnbRIvLq0vPJpd+rkXpDEjTagez2d7WH5a8hCWWG52IhX/NvEbSBtn1cSr8384Q8AAweKN+Lgidqd0GcLc57Fa9ETDRIBiVff8KXhqIRLpYtXfuVyHLcBUX2a1LHpgi2XPq+ywSGs/6w3KeAtmK4sryKB5NLyyj8/MrcBFtM2GDfOq64wCNPV7RNU4pXPA3sNVTvD79n7i7K8mohXV24DaU7sfLacYSQg/HHhF+u7eM0j6/JjefX43W2VUi2WV982OkjK8sqfa+pD68LyyrtByI7jj4myvOL3WYrXpCyvSJI+r3zb33DD6HNMB5S0Q2XF8Xll8ytq23zaUZZXl+KVLK/ZEdUGfXcbyCPrKdoAkXfxqto604VIcSGA07C8xomXKzs3js9rVPuIWrAVfh9nhy3W59Un8WqzYEsmYpJ0G+DF63bb+ec2oCsMTOO8ivLGfhbl9iK6lo4Y5cVrmnFeMW9p9nc++7zqoGt5NY19nDbt2kHuWLvWHw0SAVlefcP3hlOpPq+y1chxt4flz40bbSDqXJc7bEWJVz4/WYtXG7cBWZ6TXLDF5xvTHjEie7cBFlPxqlv3KCLjWl5FltioumevybsZpGF5TdNtIO/iNURnJzJf2W8/gG7doBotrw0HHQQVL16nTp0Kv/rVr2DYsGHw/e9/Hx5++GHwgj59ss6Bv7jyec2L20DcOK9xow1ELVrgEaUt8nnVcRuIijbA5y+PC7ZEeT7llHTFK+YhyrJrKkhUGwEkYXmN2jggZPVqKIjSDvMryqfI0sr/HdX22Gva+rxStIF0iSpvX6MLDR4McNZZ8jbmc6SddbRgK5KGhgYYPnw4DBkyBF588cVAyJ533nnQs2dP+PGPfwyZsvfeAB9+CDB+fPrX9t1Z2pXbQJ4sry4XbJmK10mTsrO8qvyeRfnJ44ItfiC57bZSHFbequPSbYAXdVhuUQOaqXjVaWNpiFdRG2J3DtNZsBXl8yqzpLLXt91hK2rb2TT7zDyMISedBNC9O8ANN8RPK6q8fRWvYV3L2hj2W77mfc4cyAuZTQHmzZsHW221FVx55ZXQr18/+N73vge77bYbfMoG5s4KbHRxG9fGG5uf8+abpQDpPiPriE2Fg28+rzILYloLtnTL7+23zTt80YItV24DPlleTa+P9cC35c03F4c2StLnFa8VVf9xF2yJ2rAL8SrbpEAVHmrGDPduA6K6l20Zm7Z45dOpdLcBri6KgwbBep1wcJVieY0Sr3kMoWVCSuNAZuK1R48ecMcdd0CnTp2gWCwGovWTTz6BnXfeGTIHG1dcfxpbZ+377gOv0Rko8ihe2Y6SteTEtbzqLtjSHdzuuit7twGZH2DWllebBVuyc9J0G9ARr3EXbCXlNqASlLL2MHGi2YIt02gDUeLVxAc1j5bXLMWrqGxraqDBZvLni+WV3dpXh/BZrlbxWkjn/qwUxLp162D+/PmwvrFxofhcs2YNfPPNN3DggQcap7fPPvtAfX097L333rD//vtLj8Pr4E/SYPdZQN+sGNcqtm1rd359PfhMMLQUi1DD3VuxVSsoGAj+hkKhRRpNaTV+nkZdN11z/fqm+kKXlhp8ABsags+L69ZJ8xqVJqbBtgNMCwdSvm0UW7cW+wI68Fkq1tQ03xvWEd7r2rXS9tnQuKimIMhnkF6h0PR5kUmH/VwIDmoCcWJb33hfQRk3gtc2eeaKjT6v7DkNTHpsnQf146g9Bs8Ke83gnwZhGwvzo6ov4TXWrCkrGz79BkH9YnmE5+j2f8Xw2Wks/7LvuLLl3QbY+m5AIYLpMM+hLJ0Grv/BdodtXnhe6ULl7bTxXnWeabZN47Nj1Q80jls259pgm08n126c7LeoH3xWDfIkap882Je6eiZVFHv3Nnv2sL4b27Gw/zTsp1xR3GijYDe/wpQpybcBH8XrG2+8AZdffjksWrSoxXfdu3e3Eq933XVX4EaALgQ33HADXHbZZcLjlq9YAeuWLIGkmT9lCrTp1g1qY1TyisWLoUMKeU2buVOmwIr27aEvd284kWmFOxNpMmP8eOgTUT5LRL5xCVFcsSLoLJEZY8ZA76VLg79Xz5kDC8aOhU0s6nLV7NlBmWzAnLtgyhRoO3MmdOTLr6EBWiXkMrJk9myoa7zejK++gtXLl0Pv2bOhneSe5k6cCEtGjoQNJ06EjQTHrJo3r+ncZdOnQ6fG35fPmdPivlgCAagILG5a38vnzoXZI0c2/d2tvh46G9QTDvLFNm2gFXPORCa9gWy9zZolLAsbFs6cCRsyaU395htovWgR9BakH+an+7RpTXWow+rZs2EGcy8bTppUln9Mt8f06VDLfIZtcErjOb1nzpS2D5bFM2bAvJEjoevMmdCFOR7rGtNrrUiDre+ZX30Fq1asgI0mTy4rG2Rtu3Zlz9DcqVOhO/P3upoaWDJtWov6WY+CEZ9rnICuWtU0uV46axYsGz9e65leMH16U7qzxoyx6geCch01qqw9JcmM0aMj+9akmDd5MqxdtqysnFbNnw81rVoZPd/YPvvMng1tFfdhOubYMm/+fOhmUJ4r586F+pEjoTfTT7aYXIjidydMcdUqqD/oIOj95ZeJXufb2bNTWfRuLF5vvfVW2G+//eD444+HI444Ah566KFAyF5zzTVw2mmnWWUCF20hq1evhgsuuAAuvPBCaCN4TdexQwfoqLMPeExqcfu6E0+EmmOPtU6jrl8/KMyaBZVGLS5m2W47qOHqAWd1BYNXY7VDhrRIoymtYjHo6Opqa6GQ1isWZnFH7VZbQU2XLoHrSLFrV+g2aJA0ryrqcOFCu3ZQYDqqus02CzYeKEyfXnZssXt3KCTkKxS0xdGjg99rt9gCYKutoNC1KxQWLhQeX9u/P8DQocE5BcF91/XqBYUFC0q/b7wxFCZPLv3ep0+L+2rhSiPotG3ru/aCC2Djvn2b/i58/DEURo3SPr/YoQNAhw4l61wjQ4cObfqdrXMstxrJzkANGJ2gXTuoufNOrevWDRxYVq7bYpzXOXOEbSzMT+Htt4V1ISNot8y9FMaPLzsf0y288UbZZ8XaWugSXg+fZ6zjQqFkeRL5qTa25z6Y1qhRUGDXK7RrB0W08kgs7Xx9B/3KkCGldPi+pWfPsmeodvPNy8qq2K0bbDRgQMvzunWDAooHFK3M81236aZQxPrUKM/aAQOajsNxwaYfKHbuDF2GDSv1KSlYpYK2msI4Kbz2gAEAvXqVPzs9esB8nHwZPN9B+8T6W7xYekyxY0ejMceW2sGDjcoT+8QemP8ePZr6yRb5TkF0t6CmBuoOPTSIhlBz+eWQFIVNNoEJAP6J1+nTp8ODDz4Im222GWy77bYwd+5c2HfffaGmpgZuuukmOOSQQ7TSQUvryJEjg3NDBg0aBGvXroVly5bBRmji5sCGn4aYaYUPBF5/n30A3nrLKo1C+/YV6dvSCv0M8Ye7t4LgM2U66AMVcXxa9c0HrG+FA12jH2LwesdkgQdDkHcuEH6TQOXLT6M8bCngRDC8t3AhjSJAf6vQRw8HWsExbF0H5RP+HtUGsFwV32vVN26levHFwQKqVjhQGqTf4nphnFfmnODZb85Q8+eK57nVXnuV/PM0r13g8tkK60eSd1l+Iq+B+eHPFd0n2zbZRVbh/bRpA4VbbwXAt1AvvADwySfl18Hj8Yeve3xlb1DfwVXDuuCfDb68uGsFeRBcqylv+DlTP8Hnmv0Ve60gjzb9QPg8YR5T2BWqlWrzDRNwIsNY77WujfcoqJ/AVcagP2+l6H+a0o343hV8e4siaK9hBBFR/2k5nsSmWCzVz7BhpYgQ8+Ylchm0sqeBsamnrq4OVjbu+tO/f38YM2ZM8PuAAQNghmR2LgKPPeOMM+Dbb79t+mz06NGBaBUJ11Sx2TXG9zimrpA9eKYLsHwLlcXHeQ3rPu4OW7qhslwuCOLhoycgqnuKCqclW8QSVacuLMvY9tBSiRZk0XeuNikwqR/TjTv4ckhjwZZttAHsxzp2BNhmG/E9ykJlmexiFRXnNWpTAlmfJKsXk3y5WrDFp5WHBVtduzpbsIX+x0agi0FUv5tWSDDTeovSDxVo1MoC49EEQ1pdddVVMGHCBNhll13gr3/9K3z11VfwzDPPBD6vJq4C22yzDVxyySVBWm+//TbcfPPN8Otf/xoyx0VnI3klm3tkMSl79jR7KH0X92y0gTg7bPEdrCy9tMSrTpzXEJ1NCkzivCoEHvqe4qvuQCipULUxl9EGTOrHdOMOG/Ead4ctUZtTHcOK1xBRHmUDtekKe5NoA3w+ZNED8DNRvZkIa1ehsvi0TDA1DLgSrzYRYUQxd00nMgi64ETdR1pRFUzznodoAwUP8pC2eL300kuhb9++gZUUX/lvv/32cNhhh8ETTzwBF+PrPE3wtcB9990H7du3h8MPPzxI95hjjoFjY/iZOsN0y0NRQ0D/pkpE1hGhtXz77fMbKks2GMYRr7oWr6TFvKl4tbW8RnXyiudpXW0tFH//e4A//hGssYnzarPHvei6ccSr6famScV5ZQkjh7D3LSorleWV/0xVhqFY1onzyotVWZ+E54lCNNnGebW19IV5s33zYNo3uLJI2uRXIl7RbcCIjz6y25AlCUzzHlVuPgjHQqH6QmUtWbIErrvuusDHFbnllluCKAFt27aFcePGGaWFu2ndc8894B2mllec3fMx5w4+GOD116HikA0U2GAxwLuuj5RvbgOywTGO24BI+MrixqZleQ07fJ2OP023gdDfMqrjV31vY3nVbYcqAZGk5TX0X0zDbYD9W2R5NRGv+Dkf6xoXK4pivLL5tXUbMLG8ph3nVZSWCXgPja56WrgSdbb5Fe2AZiOEs95sIa549dltoJB/8Wrcon7wgx+0CJOFGw3MmjULjjzySKgITC2vfAf5058G8dQqEtUOWyaBqH12G+A3KXC5w5bOFq6+WF51toc1Ea+qOo961cYfZ/pdXLcBVf2E4vOww/TSEu0YJcu7StSl5fOq4zYgEpjdupV/xkSGMNphiy973rqn8nmVuQ3oImvrabqhmQb4z1K8uvJ59WmrdFufV9vv06CQf/GqNWI+99xz8MADDzSFOTn00EObLK+sRXYghpiqBExf82Dnwsaww799fi0eB5V1zKSTzYvlNe4OWyLLq2iATnJLPdMFW6odbPjBycTnFSd0M2cmJ16z8nkN83TccQDPP+/W8mrioxxXvLLoWl5lwkwkXlVxPlUine9XeLGq8nmVuQ3kyfJqOtF3JV5duQ00RhvIrXg1LQdJRJmqEa8poaWwfvazn8EGG2wQ7JyAC6xOOOEEqMWQNY1g+Av0Xd3VdBu1SrG88h0kdjaVKl5lHT92mCavvn0uH1duAybbwybZmdhYXp97DkAU15QXDiaW186dpV8Fu/CkJV7DmJ98/TBxUbUFhMp/EsWbKByN6LWqa/Ea5fMqapvhMXit8Pco8Sr7Du+HF6877tgi1FaL/IpENt+/CsSR1OdVZnm18XnNyvJq6lIUx2/cxYItUVuwufe8ug34umDr0kurz/KKwhUFLNKnTx/4zne+A619Fh9pi1d+YKtW8WpieTVt4JtuikGGIVduAyY+r2mLV5Ugxzz/+c96g5OJeFUtYtS9fxfiFdupKDTTb35jJ15l4GYPuuJVlk5SbgONWxcLj2FjkUa5DciEGQouvj/43vcA7r9fnV9d8cqWl8ryKqo3W59X20ms6XgSV7y6it9p+zZItGDLpn/Lq+XVV/G6667VZ3ll2WmnneDNN9+E8ePHB9uzhaxZswa+/vprePjhhyH3mM6U+Q4yDDpeiaD1TNTwTcWrKbgjyPDhkKtoAzLrluh1fFpuA6+9VrIwqqwaqvvlJy8m0QYUltdEfF6/+12APfcEuOGG8s9xIVG4w02Yf3xmVX7qzJumMlR9BPp4iiyNJnXtym1AlK4L8SorB9zUYZNNyiefGAYNF3WKFvaq4ryK3Ab4v2XCWuY2oIsLy2va4tUVDn1eq8rymge3gSTxyfLKgtvAPv/887D11lvDl19+CcOGDYNp06YFO2bhdrEVgWlnIwoJ47NPZxwwlm8W4jXtBV6u3AZEQnD+fEgVdjDBLTwffdRevPKvuNl0ojotldsAb0mT5cHE8oqWXlG7YVfBh3UblXcMBffznwO8+GL55/xzfsEFAPfeiytb1WI9Ku9xLa9hu2UnYXHFq6hPC9PnJ+tYxrjt9957A4wf32zVlvWLppZX9jiZJVXlNhAlYAcPBsAFyGwZkXi1X7CVZ5/XSnEbyMJYkiDGV/nHP/4RhMd6+umngy1iMUzWW2+9BQcddFCwtWtF4CIun62fj+/IrFNoOUlSvKb0QDRdK6loAytWlH5E10wKPu1//COeeJUt2IoiC7cB0fFsO9UVr8iJJ0ZfE1+NP/MMwCmnqMNfubK8/uQn8nPZc3TcBkT+srqWV74cwu10zzuv5CqALhSI7I2UyvLKC1C8Fns/Mp9XvL6tz+vvfgewww5keTVF5L9tG23AF8ur6wVbhBOMW9SyZctg2223DX4fPHhwYH1F/9dTTjkl2CWrIoj7mifs8CtRvKoWluiKV1y4YUraZekqzit/rswfLS2f15Aon1cZKp/XKHTcBnTiwcrQfZ3MttMw/6Lj/u//Sv+jm4XsurKwcaL8sOV70kmlfiK8hq3lVVUeKvGKrisyv1gXPq98jFf+eFledSJx8PkwjTagWiDHX9Ol24DtBDWrccSl20A1hsry2fJa8CAPabsNbLrppoFva69evWDzzTcPxCuGzsIQWktVoVDyhCsfJez4w51qKh2TaANnn+235ZUd3NiV1zbw5+Jrex9e49iKV5XPqwosT13xarsoS3Se6DO2naosr0cfXbKk9u7tLj9h+WIs6IMOarZE2kYbUF2ftaDydbp4sXyTAvY8tg9UWbZ1xausT1VFG+jVq+U12eNk4ihOtAGReI27w5btOoisxKurBVuFgp3bgC+W10oUr8UYY1pefV5PPPFEuOCCC+D666+HAw88EA455JDA8vrZZ58FUQgqgrihTSrJ8or7zeOkBEX4AQfIj9P1eUWrq83WuWm7DbiyvCbZSehiE/9UBm+10i0bPAf9pWWX1O3wTMI1iUQKLrwSuT3IhC/uCqVC9YxHlTsrZmT3H7Ujmqrc/vpXgKOOEn8nEq8hvDBUXUsmXtFtwEa8snnCtHECsccearcBmSVV5TYQVTei+8rK57UCLK9W4tUXTPMeHu+zeK0AjMXrL37xC+jXrx907Ngx2JTg3nvvhWeffRa22247OPPMM6EicCVeKyHiAFpQMD4crhDeZ5/44tXmwc3iYWfFq+sdtlTXy4N4ZdOT7cIlykOHDrgdH/oeidNlj1VdX3WNqOMvuwzg2Wf13AZ0sLW86h7LhzYLY9RGnYc8/XSzeOWFKu6SKFuwxX7O3l+SlleRewRGrunRQ35NNk0Tt4Eoyyu7U5qtiwx/PT6tPODQ5zXXgs02VJbt92nQkKBLRkrtXEtdHXPMMcFGBCpwy9jzzjsP/iyLD1mNbgOVIF6xIWJ4G/xx4TZg07BlPm1JwXbArkNlqa6ZB/HKW61MO0EUI1HiNUnLK7+bVVzxquojogSpzrG85RUnxrriVYXM8op1Lwt/prK8ihZsxXUbkN2byPJqsmBL5fN63XUA220nzq/uRI0nTCNv1keyvIrLQbb5iK6Ps+4zi5N90eJeIkBrVNtll11g5513Dn7QzxVdBDbaaCP43ve+B/vuuy/07t0bvvjii6aFXLlHNVPWCdlUSZZXXeGju2CLfXDPOEM/D2nOVtlIEa63h1VdMylM07YNlaXTjkSWtDBd0e+ydKLSEOU1PJ8VNGHMXd0ywvBJuvkxmTTo+rzyfYruNfg6RcurrO5kllcTtwHTBVuvvgrw1FMAq1apr+fK51UG/x1b3nEtr3EmwVlgK7ZFz1zehLuqTZxwgjrKR1S0AR/6+mKV+LyewYiM448/Ptgi9kiuE8fNC57BEDGVgMryip0ZHxJMNGjKzs8bug0RO3YdYc92BPvvD7D11gCnnRadh7TFa3g9mQDlX98iF10EcOONdjvdJHl/pu0wqQVbCL9laHhJF+JV5z7xGHaSFT7LuuWPsayffDI9twH+dTr/jOnmm5+AqSyvSfq8qib0bLny1xblk82DLG8ytwHZhJS/pgvLa5gmThiSBjeCCDfgiIurt2S4YMtXl4mNNwZYuLDkToQb4Yjg2zaW8cknA/ztb/65DRQUMbJZfFiLERPjFjVy5EjYbbfdWny+/fbbw9ixY6EqxKuvgfWTQLfTwY7d5qHUKc+sxatooOvTp+VncSYrSXXuuLNRkj6v/AIb1XmIzDqfhOVVZLHHv0VWQZMyYs+3WbClEw5K1/Ia9VzI6nLJErl4NfF5NQ2VZfKMqAS9TpxX/CyML2srXl1YXsP7cCVe77qrtBnGTjuV3Byi+iVbqsFt4OabAZ54ohQOTwaf96hnLksf5xrLNzEuSWmsNi5d3FnroYcegtXMFpcY+/Wuu+6CoaoGUCkLtkysi74+sEk8DGHHfscdzbErddLTaeh4TNox/6Ksiz17tvwsrpsIbifqig03LC08ufLKZKMNmLoNyMrIhXjVWbCFx4iElUmHyz7XKj9vE8urDN7yajqQykJQyTYpUIlX1bXiug2YHCt68yWrewxzhpsN8Onq7uDm0vKK1m4b+LziRjEYwu2KKwC23FJ8LZ8WbMkmFz4gW9THH6P62wRcQJlkWRQ8cBtIiRqb7WE//vhj2GOPPYL4rhgqa6+99oKJEyfCdfwssJosr+HuOxiTMBQh1ebzigwciCv85MeJxEQU4eu/MN7mfvtBovCvxvkdh0aMEHd4cS2vmC4uFkHRyce2NGWrrQCOO67kY2raWaomCjphhlT1Lpv8pWV5lQ1WttYCV+JV1/LKt7GougjjTPPXVPlym4hXmZUJt9MVwef/wgvNy48XryqfV4RfbKqaDKssr3HFq6u4pWy5JCmEMDIID/YpGDIRy0XUz9u4DWS5+l7n2q7ivOIzgX17kgKz4IF49SnaAAuGx3r11Vfh/fffDwQrgou4dt999yDea0WgMvvL7vFnPwPABWv42iZsQJVQHroPg27HbiNeQxFyww0AY8YADBsG8PrrkBgq6+Ill5Ri1X7yiVvxGlqJwgngBx9ALNhyZd6SZNY5RVhetX3iTH1ek7C8su1B9SZGlleJ32+qlleZeOU/j3IbkK2slr1F4POvmqTJyo/va1Q+r6HrDJ+uTLzy+XO5YMsWvu7YezWdzOiCMZn5SR4Kr5//vHTN4cNLK+Efeyz+gi38PqsNCXTKy/SZk4nXXXcttacsFmwVPAjP5RgrddWmTRv4/ve/H/xUJCrLKwqn+nrxOfwMvxrdBqIcx23dBsJX4ehvncYrD1aUsJaeMC+yEDy22Ih6Fez5c+a4s7zid0lYXl2s4E/L8sqWj43ldZdd7C2veByKsenTS89CVL4xLBnubGZreY2yiMt8mWX3w/eJuruuqcSryufVVLzGWbAl6/PCNPDt3B/+AMbwaaqs4XFFChpgsM3g4iXeR3f77ZvvBZ9jWXswtbxmKV51yst0zJKJV/b5nzwZEqFgYXnFt3277w7wwAPp5iEmnjqieCxejz1WP51qF68ibDpbG8EbB7weK0rY8D1hXkQiLI7g1BGvOHGyaXv8DkVxwEHfpvyjxKsLy6uonSRheWWFj0qMi9I87zw7Ky/rNoDWeVywg9ssR6V1+uniiBl4DyLBwPu86m4Piz7gP/pRSSjefrs8Pybilb0e24ZFbgCqvOGKcj5d2W5vfHmaWF5lE5kwTdwOGOvs1lshFrZuNTrg26277y65G0X1STLxKvB5VS7YytIfNurNAsLn3Tb8Yfj3L34B8IMfQCLU1Oh9zt4Duofk0DJL4tXEbeCkk0qBg885R7N0K6B4dRs1b8nQfYh0fV7ThhUl7Gt3EwuiCToDw49/bCfy8NUsLtwK/bLjWl5l9YFbCcuIcKVxEipLNFgk7TZgannVfS5Cwk0yWLcBfAOBC3YwZA9/Deyf2J3wMK/o3uViwZYI9voY8u6++wAGDZIfzwsBXfGPaWP7P+uskkBm61AmXsPP+PrG41GcYX8eteDJxPIqawthmvj9vvsCDB4MsUirP7QVr6I+Pg8+rzI3P5XwU6Upuy9sB6gh8Pl1TUHTaKSzyYfn5C/HLpH5nsksr+HnOGvSaXhsg5CtvvUdW4uYzsIc0d865yQNXo99DcqKV3YgYkGrist8isodr4nXsTkfV1zrnqvqnFHYyKwoKJowbi8uTEAfcFF+dES/qhxV3+n4prlwG2DLRyVeVdZA3c/5Hd6iJn+YDgpYVZ7DdEXi9cMPAWbNEqdvcj8yRG1HFqWGvR6usEc/y3CxJnuP+HyqJgoyEfbTn5ZCJcmu6Vq82iKqOxlR18K2/8gjeteN6qs123JRx20gK1T+w1GfR6UZVV5ZuUqIxCtOhkX9hg0p1Wd1i1fZq6Ow8FWDoc5qcL4SMYg9+pbIgiH7iG5D5FcNY9BnV76dWYhXU8vrr38dL586Fmn87Je/lAeAZ5GFedPJo2pwxE5Plgamj5s0/OlPLeNrRlheI1e1i46zsbzyExOda6owtbyailcc4NhBLkqgY/mINjLQdRvAV9oPPyxO3+R+TNrl+eeLj1WlzT4DuHDIhbAWtX1MI2yzUeJVJnJcD+Yq0aNzzzgR0MGmrxZZWUW+56bp6nLbbWbHs9d2XX9xxKvOVusmRBkE8Hm6/343ro7k82oIzhxcERZ+XF9L9nzsFPfcsxQOSfVaLW2idrfSeXBx4Qi7HziCf6Ogi0rPxuc1aTBPoi1Ew+9kFsQ4D62OVQM/w4Hnj38Uh7FhkQ3mOlv42roNsNeR1XOSobJ0fF6xc47rNsBiGm3AxvKqeo2vI14xXBYvXnERoo7vXtSkwrTcRMd36SIOf6dKm7USycSrbCDGBUn8RBsFKr41UG0+EmUpk0X1EOUNXRZ0wXpi/d1V4jOqr8S0XIkxE9cwVV26XIBrGmJQx/KqeuWum7YIWXvCtn3tteCUQiF6USa+McMQizmhAmI5QSko/vz57tKTuQ2YPvTs8apXf1my996l2Kw21o+QAQPEn2+xRX7dBqIWbCU9M1a1EXRZQb9HfhDWOR/FK3s/IlSDtCragGogcOXzqrIMpOXzapsfW/Fqannl2yaKqlGjyj9DwaeDa7cB2fEqn1Vb8Sq7Fu4uxq/+fuqp0rMhOgfbLJahzPKK5Y0uM//5j/h7Ud7QhQfb4ezZAM8/D5Fg/4zhAdEooArB6NIYoOOiIrq+oE0q3QZcvj43tRyy+UKffZ2NJGx9XlUCkgXDj7keX3h0N+kIQVH7zTfi73TTcIxHKsoSDJR89NEl658rdNwGdD6Xide0xZgKUWfj6rW+zutZX90GbCyvSUcbMOlAZefr+F7bWF6j8h8OuLq+0DJMhK2JeLWtO9OFKHHdBqIEhcjyiqJn3rzyz6ImMKL0RXl0NRE3tery4jUqb2zb4MUrgm1Cdj2V5RXXTTzzTCmqg6xMRWWEz8P++4vDponAqAwYZ1W02AuFLRoQ0C0t6jlyaXmV1Znp2zWXuyfGGSuwHPHZiXo7lXSoxnDscXmdguY6Atlx7CJQT8i/eMUZARb4zjsDnHBCspbXuG4DtukkSVRe4gxOOrFQbS2vuM1eUqh8XpOyvOp0+GxnZitebdwG0DIfV7yy8SFFsOfbLu7QtbzqlIEupiGATMUrb/HjrW6i++M/++9/W6a7ciV4b3nV9XnFe4naoOLcc5t/P/hgMCIs83C3Mv4aUZvRxF1pH/WsY7z1O+8suaXpiFfZMXyfFvVMm0w4Vce6tLzGaY+bbQbw6KOlHx/cDl1SULg+6PS9/POl2mglpfvMv9tA2HFggR1ySMns/5e/6J0r6xRMLa8y8mB5NZmVmSISCToLk3TygAuX3nsPYMoU+/yprsfmnbWoqF5/u6xXUbmYiFeZqNIR3ax4xckh7pbUuJtebMurbKB3IV51fF7xb1HZ+Gp5xQ0m+O2JVedhfvjtU0XoHBOmp8qjqzZvOqk19Xnda69S2aHLDb9pQRRhm12wIDrfIuKutDcp4zhuWFFtS9fyKnIbUE3yXIrXuG8CamvdWTx9GucLmuJVd3KJm22gK9Kll0ZfKyEqR7y6fI2VtM+rTyS5EtSV5VWWB9exVmVpi8Sr6No4c7fFVNTbWl51fML4OKbsOTKf15huA8U0La8628j6YnlFn0hTy6vuVs06pLFgS/R5VDtlQxXKrIn8oIyRXtLwa3ZpecXvTTbGidNf60yMVN/LPpO5TYTtFDedGD8enJCGcNL1efXZ8spiI17xOH6Rdsrk320gCfEapqFKS0eMphU6JS5pilebB9p0wOfBgOS6AfpFr5bZ14Vh/kUWRPSbw9d3UREcREQNFK7cBnTEK+/iwp6j6+jPX8fEbcB2Nx7dOK8uxWvSltdvvy0Xo1ETal3Lqy6+Wl5/8pPmZ/Tii80XfJmgcgvQuYbtZAz9aDF80SabgDY6bgMyXLgNyJ45/vPjjgM444zSbmyu4ouK8mhDnDLMWrwWNPsXky2gTcs1JZHumYpKWbxGDcS6llcXFZ4VUQ0tzmstrJuoBzhO56/zkKArCQblNl3QFyWyZAMaLpzAbTJN/Sp1Bm+2vUYtcnApXnWswjEtr5CW5RXv34Xvpk1+XFheXbkN6JKVz2tUuhguCnfzwk0G0KLqsk5tLa8yy65te8aQir176+RQnB/R36p+w4XbgOwe+fvEiT4uWMN7dCl2fH5VnwdqJO2RxGvO3QayaAhJo/PaKg68iLPpAF1MDkwfvqjdcqIWaUR9L7qm6Doy0rK8iqwmcRZsyUJlyVaGx/V5jfKFZI+1wVSYmIpXXBn/5z+buQ3ssQfk3vKq87ziGxXc3tWknm3QDU2Fca1xsoyr/l24DdiME3HqI8oNS2eCIbGAl7kFyY6rNKLuUXeb+STyUpT4vMrih5N4TQC+Y3HhWypbsBXH59XXGZnIOuWyk+FFYNpuA2F7MBW6ccVrnHiDor/5HeGiFl4l5TaQlOU1Kn2d73QGyLC9u3qll3ScV+R//5NfTzSxwEV2rjZt8dVtQCcfaVhe2fLB2MsY8WbXXeXHqM43+c72HNX4mJTlVTTGJDkGYujMLAnvLeqeMcb6b39rtmGFK4oS8Xr88W7jYCdIptOfb7/9Fs466yzYeeed4bvf/S7ccMMNsFq2S4kMXkQsWpSceDWtQBKv0ZZXHbLww5GJwyi3AVvLq+w6ITvtVP4K8bzzmtso+o+l6TYgQmZp1Y02oGt5VdWjiUU17rNtKz7i9gm8wJBNejAIvo9uAzrXkV3L5HyXeVNNuETXMDF6mLbnKOLcMx/+yMZaKisPVb5cCyPctOi66yAxXPm84vc77thyK23d65jAX1smXnGXLdze++yz7dtipfu8FovFQLiuXLkSnnjiCbj99tvhrbfegjvuuMMsIX4wdLHTliuH67TF6w47uE8zabcBH/1w8FjZgKVreY07+PL5PeKI8r8xWPmDD5b2oe/Zs2V6ss4m7z6vJsJW5PIgS8ekfYTbHuM5qlf0JtZA3fYSle/we1exbNn0XVg3dS3PLvpaV6+m40aUsBWvLib6OmmgpQ3j5g4fnox4Fe2wFWXRj0PSK+F1RWUWC7biitcwkge65KjO94DMQmVNmjQJRo4cCe+99x50a5zxoZi98cYb4SLeZ0gFPxiaiNc4cV5N9wWXpeOC224rCY4xYwA+/dT8/DTdBlxaXpNyG1AN/uE1o15/x/V51bHehPt48x2N7HhRujbiNY7Pq+T6TkJlmVjf4ojXAw4obSWJq8AxLmSar7J1fF5Nw8ixYYt4okSG6XOVR/Ea9y1Kmm4DNu0Lg/LjtfjXxa5e9Ued56EwagrfhfGBVWHaTNwGoo5Pk2JE3Pk4Ric8NoXQoJlZXrt37w4PP/xwk3ANWabar12nY9Hdbs/G8moqRtMSrxhfdIst7DvruOJVdX6Uz2vS4tX0ulHiNPwerZ3hzlOiMFxxfV5NBg5R0PU41mo2zqvIX01HJMrcBnTq0tZtQCSYdcWrSbvEc9GvETdvUOFiwZapeA0xsbyqjmXLKclB1sbKpzrfJg0ZtgsIffB51akzDFWl4+eok5asHuK6haTNVVcBjBgBcOaZpb+vv760mAlf82+zjV4aWYjXfv3iLdgy2UY8YzGemeW1rq4u8HMNaWhogMcffxx25Z3dOVcD/GFp4Afbgw+GwsyZpdcgK1ZA4e235entuCMUxo5t8TnmJUyzhrleA7PPeKFYDH6CdBoaoCjYJQSrsMCf35hfNt24NGBajXkzTbch4rygzAX3VlYuTHnxFDbYoKwMROnx1w7ruOl//NEoX2H+168vndvQoFc2hUKpTGprhcc34IMZ5uV3vwOYNaskHrn8FWpqIvNWli4OkEwaLdpOWMci2raFQteuUGD2r2/xXBjkq7huXXPbDs/j2nGLOuPrqFgsOwYtq+H3wvouFGB9+Gwp6lXWFkoZK6/jBhR6XF7DZ5AvB9kzHAtBm1M9KzrtswEHIVU5h+VTU6PdFxSx/UiMBmXtjrtWi++1Lsb0m8Vic31wda6sZ0m6sfMmQfXMyPIprBMRin7JJv/suBRce9AgKOAbOTZdyTMRlTdRfmTtu+y8IGOFsrFbNpaKKP7wh1B47TXFXZffmyxvUce3EPShIQyP2XrrkoUa++nGN1NRY2aLNi0ZS1uUM9tPmo7nZ50FNSef3DJPXBuswXtgnsMW+dpggxZ5kuaRz0NKGzJ5s8PWzTffDF9//TU8//zz0mOWr1gB6zB0DMPEL75oeSCu4sO3qY8/DrXc8SHra2thav/+MEDw/ZSvvgq+RwYy38+eOBGWN4aS6DN3LrRt/G7lnDlQP3Jki3Q2nD4dNmLOnxge09BQlm5cgjJo3Rrqpk6F7ly6De3bw+I994QNX39dfO7IkbDB7NmwmSQ/K0aPhlmCe2PzP3/KFFgkOAbZeMEC6Mgcu2DaNFjIHSsriyVLlwb/r543D2YI0t/k22+hQ0Q5Lpo+HeaPHAk1K1dCf50yr6lpqqd+DQ3QihvUp4wZA+txcsQi2Dayz8KFTe1DhzkzZsBS5h57zpoFnZjzp3/zDaxRLEbs1bo1tGeOnztlCiwRlFn32bOhLiJfq779Fto1HrN81ixY1bo1dGWfg0mTYGMujfUNDTCFuR7fppbMng1zG7+v++EPocvbb8MGc+c2fY8+caNwu0Gs13nzpPW6cMYMWCBpa5gee81Z06bBmjVroK/gGey7eDG0Zj5fNWcOzJSka0ubGTNgU+4+Zo4dC6skA6ZOnzBv2jRYzOSz08SJ0JO9j7lzg/vQbu+4U+wGG8AGkmNnTZgAKxrfQvDXavrewErTZfJk6Nr4XOPzHdbHhtOmlfWVWEJse4qisHp1i7583pQpZWVlS885c8qeRRZZ3zQwos8LabVoEfSTpD3l669hfX29UV6719eXPd/T9tkHNvv447JjsMz7dO4MbadPh9V9+gjzL8obOy7K2mzQvtesKft81qRJUOzYsak/R76dNAmWYXQGfN5nz1b243M22AA69eoFHTgRLqJpjNV8ntjjTZGlv2Dq1KC++Ta9cNo0Yd/VbvJk6C3oowYsWgSFxpjNa7t1gw0Y44SMifX1UHvQQdDjqafKPuefpwGLF0OhsR9aPGMGzOPyVbNsWVn/MWvcONhEkEdRGUyfNg2gTx+oCvGKwvXRRx8NFm0NxoUoEjp26AAdMTg1w9ChQ6XHF956CwqSbecaLrkEuuy0E9Rw6SFD0Nm7c+fgd/b7WvSDabxeoXt3KDRWXF2PHtBDkA+8doE5vymvaCUQXNeWod/5Tum1QH19i3SLtbXQ+ZJLoFBfD4Xp01uei3maPr3svGL79lBYubJ0b2vXQk/BvZWVC76qkNRDoW/fsuvW9e8PfbljCyefDIVnnmm+frEYdHR1tbVQwBl7jx7QTVS+vXoF96Wirk8f2BTPXbFCr8xbt26qp8KWW0Jh3Liyr4dgWUf5POG5PXtCgemso6jFmT1zj4V//xsKEyc2f4+vqhTbzxYGDIACBrQPj8fg36Iye/99KHzzjTIvdWjFXbiw9Hvv3lAcOLDFc9CinXXpAl3Y682YUXZMXb9+0Dv8Hv8//XSo+elPS+cWi4CyfMiQIdCqVSsobLIJFNCiLcpbv36wmeyZnzWrPJ9YZr17l33WVLfdukGB8e2t694duiv6Eis6d25RTrUYxkry2lGnfbJ9UMCyZeXl3LNn6T7QmrXxxlBYsSIyzWKfPk2DZIvroTtSeL0lS1reD/bXJuVWVxf0y/h8d/rZz5rrY+zYsr6y2LVreXuKYs2alnlDtx4HdRqMI8yzyCLrm8rqZMCAFn1eE/PnS+s9GIdw9bdJXt99t+z53mbvvYNnqfDii02fBWV+990AH38MtTvtBN26dBEnNm9eWd7YcTFE2L633LLs846DB8NYFNWN/XmLsRSfd94gwKY5eDAUVq2K7Oub7k2St6jjTZGlj30U1jffpuv69hX3XW3bCvuoGqyXxuhLRTTIYb/85JPKPA3FcwVtin+egu8bxWvdZptBHz5fq1eX96WDB4vzKCgDHPMnuH6L5aN4veaaa+Cpp54KBOz+uOOGAmz4YeMPwcFOimK1cXCexCeuFVoawnSZ71t83vgdvlYS+kXha0v2/PAY2V7clrQKr8NdryxvIt/FME/ceQUM3dHYAQYCPSIEUSvZ/SNolWHTDsud5ZhjAPbZp3kVN1ffwnNKmY8sx6C94HFYdzplXlPTXE+4IIqb/LTC+9HxZ9W9XpguvqZi0+XbDv6tui4KavZ49DXWaJMiglddYdsWtI+m9saew7cB9KNk6519diQr2bHcg7JXlF1BVQ54z2w+wzIVPYPcNQoug9qHCMpJWY8a7aUVlquinZQ9Kyg8NSxLBRRIU6aIr8fWm+n9iNh8c1h/5pmw6MMPofakk5rrg6sn6TMvQ5Q3UZuzgWtXLNJ86t6Loq1b5V/Wb/DPANY5Ljw0SSvqGZZcr6ZxVzt2/C5rN6qoIOx1dZ4PNn8YsxQ3+MBJ8ksvRR9viqxNhH0h/2ya6gReY0SUk0rXCK/N9vH8d+h2yedJlkeOGt61KSEy9Zq+55574Omnn4bbbrsNDnIVl9AFOosqfIo2oFrdqLPykf8Od4kJF1rx4VN0zmfRddS3WbiQ5IItRBSCSncVt+kK5ajVnVH5x85G53idfPHRBnQWbvCf8fdjElrM1SYFnIBWphO13a4NScQejVqwxf6Nu09FgZY03QVJOvWuw377wfwf/7j8DUbcle1JxnnN04ItUTnalkPcRXSqHba4iasSTeHagsMOA3j2WYBf/QpSRTbmxl2sqXtcQ4O7WLU2pLQYLzPL68SJE+G+++6D4cOHww477ABzGf83jESQOMOGJbu6PYsVla46ayx/3DMc6wSD48dBd/WrTSzEJOK8ssfxQlXDatmE6YyeF3umq4Z1xWtacV75+4m4btng5mp7WIXFrEXnnYR4dSX2TLbDZtMXhVDjwZBfqntPa5OCuELJh+1hZdhGz7Ap27iTgLhpiY6JshhGpathcZQSZ9V81uKV76Nc1mUxItpADshMvL755pvB6uL7778/+GEZK4gA0IIDDwTYd1/7MBiq7TVdWUx9EK82oa7w7wEDSj9xiRtuTPWdSX3YiFc+ryaxMyU+hLkTr6IwNzaWV5OysxWvJpZXnkq0vPJtQgSGKlS11aj+IynxappuEhOFuJsMxD3fheU1jvCzmVDoilfT0HR5FFg2b/tcUCwmP0lSkVJdZSZe0eKKP1agqDr1VPuLR1l2XVj6VOkkha3bgE46PPjKb/ny0u+NERgSE/Bx6iN8kG2sMCJLni6NC56cuQ1E1QlvYXAV51V0jo5Q4Ms7qvxVk4aoa8uugWWq63uVRHiXNCyvLsTrjBn2lldXE3QXcUDxHH7i5QJby6lOPlTfueizEFfi1TYd0Q5bhudXhHhNy22gaCheZWCcW4z+1Liw1jc8jxScEFEPUlIdctLYDvo2ndSVV5YEU+/eagu4Lz6vuriyvEbt9MaLzSjLa1Y+r7aWV9PrstdQDdomQkIVfJ+3tObVbUD17OKCtSj23FPfbSBJAeFCKNkE6LdJV/c7nXwkbXl1WWe241nUG8GkfF6zwtRglLV4rZHU6+67l3bv/N73vCz/zKMNJIqtNUXH0mfr/5M2LvymROBikMceK1kjTTp3l/nRSStsA7rXVdWxiXhds0b9PVqrV63y321ANy86A5CCRLaHxTJtDPlWsW4DqkmOjuV1yBB1P5mV5dWFeHXl85qV5TWpvjvJOhG1pagJcByf10MPBXjhBfCKMK+uLa+6FB1ZXj032lW25VXWKdv6cbmKNpAUtoOlrT8SWg9Nrdg2LgCm57CDjekEhj0+juU1Ci7Yd2y3AV6oyGIc6tRtlOXVRsREWV7Z+3Pl84p/y47PasFW3D5Bd3tYHctr2GZ8sLyqrpu11TFJn1fX52Xl+2vqNmAiXlWW18aNDnJheY0rXnXLv2g45uXU57WyxasO11zTUki4ek2dF8trkvnU3UNe9cCZPPQqwaKDymfORLxisG4VvJ+wiUVNR7zKdjhxEW3Apr2YiFdXPq8qshKvcYmIs1lWdlHiNUwrjnj12fLqS6gsV9Ez0havDn1eY+VJZXl1HZvZ5wVbcdpzjeN25gH5zHVc2EaAO0U88YT8e53P82B5NfG7cplv3QHFZsWz6HO+M4uahZ50kvx4Pi2TBVtnnaX+np8wRdWBqeWVTz9p8Wr7NiNJy6vquny7SGvBVtxni59Aqco1ym3AVLwm2Z+56INcLPrSSdelcQD7lJ/8pLQAlu8zXAh423RE54nS2WEHvTzFmVioLK+m8bTTJKnJVBTFBPoyE8jy6qBydN0GXJv3bY+Pi+3CgCTRHVBUwtCkPrAzY7c7jBrAcSUluzWeK8tr1N7OqggNNlYPl5snsCv08bouLDDr1qXv82qS17R8Xl2nqfo7aqISfs/Wt8r3Om+W1zTEq841osark08GwL3o99oLYpOk24CIiy7S83lV9SFxfF5dunK5QsdglKSFtqEhH9ogJtVpebUla18Sm+tFvb4WkaTlVQZuV/jznwNssgk0nHGGvXjFARlDfKBIQ4H4i19EX5sd5F2J1yjY3YVcDOabbAKA+8wj55yjn64Ivgz41882Pq8R4lX7WYm7eCbvPq8uA5eHeWHT5CeRUSLDVR/nwvK6dGn8NLJ6NS2aJNqm4wqd8tNZFCjaYYsljs+rj5bXpBZsxanbQkzxanLtlHSPhzWfAraFaype03Yb4O8L9zaPu71rXExWAJ94Yuln5kx1GqrP8bPNNwf4059KViSdHVbYdNISr1FpmVqRsA5xV7TFi9WLGGTl/93vArzzjljI4avB/v0Bpk8HuOQS+fXjiFfdtPLk8xrHEnz66QD33msuXk36HJHbAD4zrAiMsrz6JF6j0vTdXzCpe3blNuBSmJvcq2pnQ599Xl2K0jgT87THf3IbSBBR4f7wh82CT4bpAJe1eL3lFvmqc9V5SXbIOp1NnIc+nInjnu26WwPKohMkKV7xmhdeWApVdPXVbl6j4jlRq2/58r/77tKGH+ymH7yAx2vfcQfAn/9c2i7YRpRFbRagu11hlP8h7meO/r5Ytib4tkmB7Dnhn+colw6VZUomXmX5cOlPGZWOTbr8VtZpWF518qnbtlzk19fFN3HcBlTi1Ue3gaTcEXWPKxaTe2OlA1leE0RUuKecArDrrgBbb20uXl35yMbFlV9Lkq+ebMRrnAVbOsjO4dM3WbAVBXbIaO3EH508uaoTXsz061f6YcUlv2ArLItwEZiNKIuwvBZcDfDHHQdw7LHN+ZE9s1m5DdiK1223LVnA0b3GZJJz/fUAV1xRHlOYP5ete9VObUlaXl34vB5xBMAnn8jTzCragG7bTsptIGvLKxZBnFBZKvGaJ7cBU/r2BRg/vnknvDji9ayIBcRRbLZZ8+/bbQc+4Ok0LWEriagRoDDBmbvKF5G97uGHN/9+yCF+Wl5tv0tSvOqUSZwZq4tVyrLPTWf5cV5xm8TzdH2v/IIt23L/8Y+bf2cXxSVtnbKxhvngNoCTCGTTTVu2jxtuKFmVo9Lj/0afd9wNT0R4DbYsePHKppem5dXmOeYnl75EG6gm8SoLkB9njMmb24Cpz6usfZ1/fmmyiiL2Zz8zv37IPfeUXOlk6NRF9+4A550HcOCBpXypIMtrgtgWLjvA7bYbwAUXlBqeLFxIGpVoMoPVSSNOOpVueTUVr9jpykKARVkM+Gulufgk6hW+yl2D5eijS4KiR4/SjmwqdMWjaTnga3b01508GeCXv4x//bionq2rrgL4+GOAnXcG+Ppru/RMJnSiaAMqAZimz6tNuqrFZnFI85VrXFzmJ8m0Klm8Yn/jovxw6/U//KHZdUs3vX33BXj88dLblt/8piR+Veimu/fepR9P8NDmngK2jYof1HHP3ySuYwJ7DZtA/3GPTUK88pjMWFWi8KabAO66C2DGjHTEK6YrE69R5WASz9M0TyLiWrdF6WKUguOP18qWM7eBFgkXAG68EWDKlHIB7aPlFa0sBxxg9jrUhXhl791UvCZl3cyL5dWl24ALXL7tS0u8xok24Jt4xeg26OLjKtoAe3+69dGxI8B99wHMmxdtNDBJV5eUJnT5dBvwwfIatyNBIfWd79jlQ3YNV52ky8anu8OW6hiTh17VmeEr1HPPjb6e7HNT/ypVXkwtr1lZkXQ/i7twIinLaxjOB+teNdFLQmDgdfnBQzf/uoNyHPEafs6WvWqntyQHpSQsr3kKleWKJK3jMSjGEa/Yt/goXjE8IQ/625u2addxXsPX/HyfJ4PEq4fIBqS0xKvqOtiwzj4bnGI6AIcrzH1zG+Bx6fNqInjjiledld6656bhP61bnqaTBh3SXJEtul5SbgMY8soG3fLUmeiZiFeVBTRNy2texKtvllcffV6TtLxmGV0Bfcl3373kz49CFsM85iHyg4eTGxsq223AtXh1Hec1au9x02uYug2g83W1+bzKYsPqHGsq0Hy0vKpW/eP98iIujs+rAWVuA6p7TaqjTUq82rx5MLlPF24DqvBw7N9pDsYu3AZctZW4C7bSxKV4jVPfuCL9yy9Lv3ftClBfb5+n0OdThMsIMKb06lXaECfJCZnL82WQ5TVH+GB5DVf28rtJxcF2hl9J0QaiRFQcy6tL8ZqVz6sq3qquldXVQrm03AaytLza5td2hXpcn9c8W151J72VbHlNMiKECbiYGRdIYpSMdu3UbgNRqBYr4YJM3FYX+31cDe8TeRGFhZzkkyMndu0UQmXJwPivodDEVcCur7P//pCZ5dXFsaYdv43l1eSVkY3lNSnxqhLSUYNeUqGyVCJNd+BLwPKq/Uy78nFLS1DYWl51xXSa4tX3UFk8eQuVVQluAyG4WcpRRzUvXorjNqA6BusGV9Q//bRXK+ETNwrF5Qc/aP4dN8qJ4zqR0X1VttuAC/DVOvqzYNgK3R2bXHWeOJO87Tb1MapFKLrk1W3Alc9rUm4DcQa9rNwGeLJwG1Ae6NEAkKTYtrVE24hX1hpvanlNSry6SNcXtwGTfhk3lXj9dYCDDtI/Jyo/nTo1/246hrlC1U7jilfRznA+4LPbwEknldwW0fVBtatoFBgmFMN5nXACpE0+La9xMRGVYRzXjTdOv7HpDHy+Rxtw4TZgclxUmaXpNsAKuoED1WmrznUJmw+d7TR1BYuNSJO13TQEqkt3nbxaXkULtlT5TdLy6iJUVlSaaboNsBt0mESVQSsYGi1sBYWoPtAiiTskoXAVbUUdl3DzDOyz+H5Ohmm78XHBVhSu8+ayX+zUCWD4cICDD46fFu7+lQHVaXlNy3oTt/GaCj1bt4G8Wl5Fn7v0ebVxeZAdj2GTTOo2qT27cfcmXAE/aRLAkUfaCVWbco9aIKYr1pYuBSfgJiMYvHvqVEgU221+0wiVpeM2EGUhI7cBcdkec0zp1XmfPqW3dmkheyuCuyytXp2M5RW35UXROmCAfAFVUm4DPr+JyakvaWxScpPJp3j99a/z0QjSeE3gwm3A5rppileTjitKFJr4vPLpxxGv/EKprCyvSBgI3zYSgyu3ASYdpdvAJZeUXqfiay4XsZHDa6OATVq82r46HDasZNHAQOO4+CXORM/E8orpXXcdwGOPlfvF5dFtIA3LqyyfOFnFgPVpo3KxSsplAAXrnnuanVONllef3AYqgFyJ14Yzzyw1iG22ySZUlilxHyyX4hVj0Y0Z0/w3rtCUXcdlhxD31bsoDdXnNm4DSfm84mv58eNLv6Po+uYbdT7SsLyqEN1fkj6vui4vu+5ashx17uwmvFya/YDtAIblef/9APPnqy13Os+qqeUVwxzdfHPL4/MmXl2R91BZPuLa59VHdNt0pfr7J0yuxGuwZ3BtLeSGuA+WqZVS9RCgFQBj7X32GcCPfgRwyCHiNFyTh1BZSYlX9AX79tuSOMDyfuKJ6GtmKV7T9nll0+Gtf/xxUftz25DGwGfrNoCgpSzqlbPOxFN2TQywHleoJbWi30W6qsWJSVtes8K3/LjKp8miXV9IOm9t2gCsWQPeQW4DCZLWbM31awIMO4IBn/EV9LvvmlmvsKFjSBGb68bBRgDqDmQ2IipOqCzTdoMC9Jxz9PORlttAEnFeY1petaMNVJN41UEnPT4Phx5aWrCBbhPIWWcB3HVX6Xec2GZRXklYXpMSr+gSsHJl6ffwf1/wWczFsbzK8Nnyauu3rnvck08CfPghwC23QDWSL8tr3nAtXjHgM3LjjfnxeY276EmFiX+fC8urS0GZZ8urCJt6RReAFSuCX4tZ7FGexkCftnjVeSaOP7787332KQkyXGCksvTmzfLq6nnl87LRRgAzZ5Z+X7AAvMK1mMPY5h9/DPDzn6cnXkUhIvlj8E3Mfvu5HU+ympB17w7Qrx/AlCkAxx2nn15bD8ODpUhli9esfV579Chf5e2qI5KFFcpDnFcbt4GsLK9pC+9qs7xefHHpbcAGG8CC/feH3jrhvFyShtUmaV/OOHFe2Xats+AmyVBZrsoJt+tEf90tt3TnasI/9z6LV9ftCxdL4r3ajF+24lUUsYA/Bt8U+Gx1NXV/Q+spljO6RuadIkUbSI60Gj1aMvC18ciRLcMSVat4TXqHLRsRlZTbgM01fbO8JhltAGNZ/v730LDBBrB+8uTmz/fYo/T6es6cZOOxsuJmq63yKV510nfVbvOwYAv9eNG3HC3JrvImsryG+CZeXY9t2E9hjFjXqOpGNI75vKDPxWJNtKJiqDEVebjnFPVVZVtefQDDzfAhZ3TRWZ3oYpOCvIpXG8uryTkuIiWwnHpqaQU57tgWFYTcF8urCFcLtpCePVuGEcP0TzsNEgc3H8Ed9GbMADj7bEgFEq965RJnAHQZkULUtvE1b0gWvtqVgOo5EMV8zqN4dR3FR3TPRU/aXzi2DR5cmtxhlJSEqU7xmoeGr2rsspXZu+wCcN99pd9xX2nX182D24BN3tPyeUVrIlr4ULza+OZWmuU1a/A+sNNN+5pJppemdTTJtH3qo3nxiv6fr7wCsGqVfEFmVvj+Kl3HbYCfzIqOyQN8XcR9m+ZzGRx4YMkYgAvKU1rEmMMRp4owdRvAGc9NN5V8Z9g4rlGkGSorbkglVdq+W17x2j77NKXt81qN+Ojzanst2fVcpO2TCOPzUlcH8MADAIsXR7/qJczFq8jy6lN7yEq8iih6YnkN36KlSGWPOD5VbFJuA/wxaNkz9d/zzW1A11fIlQVQd4etPHagacV5JfHqr+XV1ULDNC2vPj1roryghQl/iHTcBirhHuP2kT5bXjPAox6CsLK8+tTJ+7BgK0m3AZ/DsuRxk4JqhNwG3O5ElgbUtv3zec0Drt3O8lgGCeKF8lmzZg0cfPDB8NFHH2WdlfxZXn3HxaDkMlSWSfpJhsrSIdwBCbdHTQOyvGZTxnnxFRXlPam+yKdJuU95qQa3AZ1oA9VoeSXKyLw0V69eDeeffz6MD/eAJ9K3vCbZMbjwcU3a5zWLTQp0wGDdBxwAsPXW6VyPLK/5I03xKkpbtLjGBrK8Vhcq8cqGsAv7vjyKV7K8Vq54nTBhQiBci0nN3vNkoRQhe2Bd+wSlba1xlYaN72WXLgCdO5cWW4RbYeqK47QtMBj7b9iw9K6n64ZBC7bytWArqWu5DEvlc7QBsrymK15xzcbhhwNg7Odf/zq/dcDfVxLRBoo51zh5Fa8ff/wx7LLLLnDuuefC0KFDs8yKn5jGefWRJPNnuz3sHXcAjBnTvHsT+bzqlR1bhjz0SiwbshR9uJuVqy0qs54oVorPa17ETFQ7PfpoyD1Zv7mrcDItzSMNd51qaGiA9QavqQoNDVAQPMwNrl51OaJG0uHg/bLfhfkurF/fdF/FQgGKce9Hch0nFIvGaYd1HFrk8V/hPXL5Dj4qJRC989luu4UXk+eR/1wn7RyDwwn/vDTg3/w925Z7RH2bPNt5IrFnS1AX+MyInpVCXR0UFi+G4uDB9v0F/zxgXGmLtIT1rXkfmWDRh2XG+vXe5ZXvz8OxLWw72OcUovIcsw6kY2yS5cOPqyhm41xPNE6v96++g7pNgVxNBdDNwISN6+uh45IlLT6fiNu1esRAQR6R6d98A5sy34X57jV7NrRv/Hz1vHkwI+79NDSU5cFl+bReuBD6WqSNu9svWbo0+H3+1KmwSHBep0mToCdXdnOnToUlhvlvP3489BLlcd26snKZ8vXXsB5dDiqUTebOhQ5cec6ZNAmWcuVZWLsWBnDHzZwwAVbF7LRGjRoFlUhSz5aoLhZMmwYLBddofdxx0PHrr2HZ9tvDets8cM9D3Hth67tNfX1ZX7dwxgxY4Es/XSwmWocu4ftEn/Ia9ufIlNGjm/rSbtOnQ+eIPNdOnAg9YtyXbIxNsnxaLVoE/ZjrThs7FtYuWmSdXscJE2BjrgxqY5ZLnsmVeB00aBB06tRJ+/jCJptAYfr0Fp/75qJQg0GvBdRuvXXZd2G+Cz16QGHu3OD3Yo8e0C3u/eCMTnAdJ8yfb5w2ztRXYCzw2looFApQi4H9RectWdKi7GpxezqL/AvziLNa5vMh+HkFi9dCr15QqK8v+6x20KCW5bluXcty33bb0taAFmB9o5AZMmQItMrTK1pNCkcfDYW//Q2KQ4e673u4uqjr3x/6yq5hu0215Fq29yKs7402Kr+Pvn1hM4/66cT6R9csXOhdXsP6DvtzZMh225XegOHz8emnwZsBZZ4XLIh1X4Wdd4YCuopxJFo+XJ632X57gI03tk9v9eqWZTDffHxNmmXLlhkbGitevNbU1JgNbmimF/jWeDdASvx/gnwy35Xlu/HzAn4W9364cnJaPuikbpk2dnT40wrPF53HpR2kj5+Z5l+WR75c2rTJl/+bKYLnpRX6afH3jK+p+OMclA2Wu3fPpgtOPhlgn32g0K9fMu2HqYsC1mFSZYjXcdhPlNU3tjP2Plz0ay5Jqn90jWrMyJiwP2/RT+uMP1z7ML6vSy4B+Ne/AL76CuCLL+zTMYHPM/qHx7meqAw2sB9fk9RpqVwnlasQ1emUj7hoyLK9kpOO88qn70HHkDq6ocWqsWxMyhAt2Eks2EhzZ6o0w3D5tGAL6dGj9P+OO4LX5GlsMGlXcdsebp1+xBHWb4e82B5WVAZ77FHarjgMrVhF5MryakyvXpBr0oo24GOoLDZPK1ZEHxNiIxB0778aBRrtsOU3aUYXSFO8+hQqC7nxRgD0J8RFakSyobJEuJrMpCnu09ikoE0bgAceAJgzB2DAAKgmPJveOgZnWugvmXcRm+fZtWWn09CuXbR4TXp7WJ5qFK9xN3ogksV3i2UlhMpCunUD2HdfgNrarHNSee02rkXShO98p/n3Aw/MV5xXGbW1AAMH+jfhqxbL69ixY90n2r49wJ13ln7/yU+gYmBXdfvWyTsSfIF4DcN+mFheXboN2B5XScgiCGDZ5yneMOE3lSLCCfP6xnE6LYPNkCEAw4cD4ILnX/4SEoXivFaHeE2MPM9GdNwGfL+/OJbX5cvNLa9JilffyzoJZHEDefFKVmk/yGsb9Xl7WMI9bP3q7NLmMnboj38MqeD724ScU/niNeT22wH++c/44WLSRPaKKk+WV8v8rccOLRSv69b54fNajcjEK9ZrntphtVAp9UDPZPWQtnhNC2rDiVIhPZ0GuNr3jDNK+yb7zNVXA+ywA8Cpp8rF6xZblN9XBQ6mcw87rPmPk05K1veyUgb8tC2vLFSG1UFoYce1BC4hK1X1oiNe87TOI4T6xESpHstrXsAVg1deWfp97VrxMcceCzB1amkgOeoo8P4B3mabUnw9VpBGsHbjjaHhrrsgGCplqyjTdhuoRnTFK7kN+EHSbfnuuwE++gjg+9+vrmgDRHKk6fOaJtSGE4XEax4bP+4yhqFb8sI11wBMmWJuJe7bVy2KXC3YkuxwVnWIBgjVgi0WmgBUB5tuWvpxDbUnN+RF5LH5rFS3AeoTE4VMTr5RibM1DBGy+ebu782V5RUDWB99dCmPt9wCVYto4CPLa77I64BJC7YqniKGrkRwo4COHZu/ILcBwgKyvOY9eHM148ryihx+eOmHB4OS42tS33fWSQryeSXSgNwGKp4i9q+77lqy3GcZbYCoCEi8+gaJV7uySsr38vzzAb78shQfsBpRRRtgIZ9XP8ir5ZXcBiofrGMMps9DPq+EBSRefYPEa/puAyqwY63mLSHJ8pov8vq2hiyv1YvOzlNkeSU4cjpNr2DyOvhUquW12qEFW/kir/0HhcqqXnTabF7FK+7mhQuPr7gi65xUHGR59Q2yvPplea12ZCvL2XaK9ZBX0VRp5LUeyG3ADTvt1Pz7McdAbhgxAuAf/wA45JB8R1EQ7eaV1o5eVQaJV9/I6+DjCyRe3bHXXvJ4nrx4Jfwgr/0HRRtwA4b9u/NOgBkzAHbfHXID5lWV37yKVyIxSLwS+cXV9rBES/r0AfjNb/TEBk0Y/CGv4pUsr+7ATV1kG7vklby6DRCJQSYT3yArlj6utoclzK0bZHn1k0oRrzh5IogQsrwSHDTS+w7upoUMG5bO9bp3h9xAPq9+iA2yvPpDXsUr/yxjIHuCCCHLK8FB71h9H3zuuAPg888B9tgj2evecw/A++8D7L035AaKNpCcdSNKBJHl1U/yKl7R3Qd3uXvrLYCTT87vfRBEUtAzUQaJV98baM+eAAcckPx1MZwH/uQdsgKmA/m8+kmeBzjZLncE8aMfATzzDMCaNQAXX0zlQZB49Y48Dz5pQ5ZXP8qe/Iz9gfoPohKprQW47z6ABQsAttwy69wQHkCWV9+gwSdeWVH5pb9gi6zd/kDtn6hU8C0k/hAELdjyEBp8CB8gn9d8QlZwgiCqAIo2QOQXEvrZQT6vBEEQREaQePUNEmT2ZXXppa5rg9CBrH3+QHVBEEQVQOKVqBzxuuuuWeWk+iDLK0EQBJERJF6J/EJWaj/Knqx9/kB1QRBEFUDilcgvJF7dst9+zb8ffLD6WLK8EgRBEBlBobIIgijx3e8CLFoEsG4dwA9/qC6Vtm2bf6dJhD+Q5ZUgiCqAxKsPXHQRwB/+AHDggVnnJF+QaHJfnj/5id6xG2zQ/DuKXcIP6JkgCKIKIPHqA3vuWfohzKCBOtu96EPWrs0wI0QZ9EwQBFEFkM8rQRDmkOXVT0i8EgRRBZB4JfJL+/ZZ56B6YS2v5DbgDyReCYKoAki8EvmlSxeAH/0IoLYW4PLLs85NdUHi1U9IvBIEUQWQzyuRb047DeDUU2nQThvyeSUIgiAygiyvRP4ha1P6kM+rnxSLWeeAIAiissXr6tWr4ZJLLoEdd9wR9txzT/gDhosiCMJ/yG2AIAgiPchI44/bwE033QSjR4+GRx99FOrr6+Giiy6CXr16wQEHHJBltgiCiILEq5+Q5ZUgKpPttwfYcEOAhQtLrnJVTmbidcWKFfDcc8/B73//e9hmm22Cn/Hjx8MTTzxB4pUg8uQ2QBAEQSRvMLj7boD6eoAtt6z60s7MbWDMmDGwbt06GDZsWNNnO+ywA3zxxRfQ0NBQ9RVDELmxvBIEQRDJ07kzwFZbkQtBlpbXuXPnwoYbbght2rRp+qxbt26BH+yiRYtgo402anEOitr169ennFMibcI6prr2l0JNDRSYV9QNMZ5Lqu941PD14HkfSfVdXVB9VxcNKRkfMxOvK1euLBOuSPj3mjVrhOdMmDAhlbwRfjBq1Kiss0BI6DJ9OnRdsqTp74kjR8YuK6pvOwYy9TBr3DhYkROXDqrv6oLqm6gI8dq2bdsWIjX8u127dsJzBg0aBJ06dUolf0S2M3Xs6IYMGQKtWrWiqvCRadOg5p13mv4cOnSodVJU3/Goqatr+r128GCsDPAZqu/qguq7uli2bFkqhsbMxGvPnj1h4cKFgd9r60b/OXQlQOFax3TGLDU1NSRmqggUriRePWXTTZv9rjbbzEk9UX3HD6ET1ENOJnxU39UF1Xd1UFNTU9kLtrbaaqtAtI5kXjd++umngbUtrZsnCMKSHXYA2GcfgP79AS6+mIqRIAiCSI3MLK/t27eHn/3sZ3DllVfC9ddfD3PmzAk2KbjhhhuyyhJBECbWvnPPpfLyDYrzShBEFZBpvJsRI0YE4vW4444LfFnPPPNM+OEPf5hllgiCIAiCIAiPyVS8ovX1xhtvDH4IgiAIgiAIIgpyLiUIgiAIgiByA4lXgiCISoF8XgmCqAJIvBIEQRAEQRC5gcQrQRAEQRAEkRtIvBIEQRAEQRC5gcQrQRAEQRAEkRtIvBIEQRAEQRC5gcQrQRBEnmnXrvn3DTfMMicEQRCpQOKVIAgiz+CW2r17A+DuhIMHZ50bgiCIyt5hiyAIgojJoEEADzxAxUgQRNVAlleCIAiCIAgiN5B4JQiCIAiCIHIDiVeCIAiCIAgiN5B4JQiCIAiCIHIDiVeCIAiCIAgiN5B4JQiCIAiCIHIDiVeCIAiCIAgiN5B4JQiCIAiCIHIDiVeCIAiCIAgiN5B4JQiCIAiCIHIDiVeCIAiCIAgiN7SGHNDQ0BD8v2rVKmjVqlXW2SESZv369cH/K1asoPquAqi+qwuq7+qC6ru6WLVqVZluS4pCsVgsgufMnz8fpkyZknU2CIIgCIIgiAj69esHXbt2haoWr+vWrYPFixdD27ZtoaaGPB0IgiAIgiB8Ay2uq1evhs6dO0Pr1q2rW7wSBEEQBEEQBEJmTIIgCIIgCCI3kHglCIIgCIIgcoP34hV9Jy655BLYcccdYc8994Q//OEPWWeJiMG3334LZ511Fuy8887w3e9+F2644YagjpHp06fD8ccfD0OHDoUDDzwQ3n333bJz33//fTj44INh++23h2OPPTY4nsgPw4cPh4svvrjp76+//hp+8YtfBPV56KGHwujRo8uOf+WVV2DfffcNvj/99NNhwYIFGeSaMGHNmjVw1VVXwU477QS777473HbbbRB6plF9Vx6zZs2CU045Bb7zne/APvvsA3/605+avqP6rqzn+uCDD4aPPvqo6bO44zW2FdQAw4YNCzTeypUrK0u83nTTTcGg9uijj8Jvf/tbuOeee+Cf//xn1tkiLMBBDIUrNtInnngCbr/9dnjrrbfgjjvuCL5DgdKtWzd44YUX4Kc//SmcccYZUF9fH5yL/+P3hxxyCDz//POw0UYbwWmnndY0MBJ+8/e//x3efvvtpr8xDBqKWZyU/uUvfwk6MBwE8XPkyy+/hEsvvTRoA8888wwsWbIERowYkeEdEDpce+21waD1yCOPwK233grPPvtsUH9U35XJOeecAx06dAieYRQg2Je//vrrVN8VxOrVq+G8886D8ePHN30Wd7z+17/+FWi5q6++OtB2X3zxBdx8881mGSt6zPLly4tDhgwpfvjhh02f3XvvvcWjjz4603wRdkyYMKE4ePDg4ty5c5s+e/nll4t77rln8f333y8OHTo0qPOQ4447rnjXXXcFv99xxx1l9b5ixYrisGHDytoG4ScLFy4s7rXXXsVDDz20eNFFFwWfPffcc8V99tmn2NDQEPyN/++3337FF154Ifj7N7/5TdOxSH19fXGLLbYoTps2LaO7IHTqeeutty5+9NFHTZ89+OCDxYsvvpjquwJZtGhR0J+PHTu26bMzzjijeNVVV1F9Vwjjx48v/uQnPyn++Mc/Duo6HG/jjtdHHnlk07HIJ598Utxuu+2C43Tx2vI6ZsyYIEwWWmVCdthhh0ClJx0Al3BP9+7d4eGHHw5mayzLli0L6nTrrbcOZvFsXY8cOTL4Hb9HK11I+/btYZtttmn6nvCXG2+8MZiZDxo0qOkzrE+s30KhEPyN/+OrR1l9b7LJJtCrV6/gc8JPPv30U+jUqVPgEhSC1nV0DaL6rjzatWsX9MNodV27di1MmjQJPvvsM9hqq62oviuEjz/+GHbZZZfg7QlLnPEaN60YNWpU2ffoeoBtCDWfLl6L17lz58KGG24Ibdq0afoMhQ+asRctWpRp3ghz6urqAh+XEJyAPP7447DrrrsGdd2jR4+y4zHA8ezZs4Pfo74n/OSDDz6A//3vf8ErI5ao+pwzZw7Vd85An7bevXvDSy+9BAcccAD84Ac/gHvvvTd4zqm+Kw+Mu37FFVcEwgb9Gn/0ox/BXnvtFfixU31XBkceeWTgDoLikyXOeI0uYKjh2O8xHmyXLl2MxnOvt4dF30hWuCLh3+hATOQb9HFBp370iUHnbVFdh/UsawvUDvwFOyj0U8cBDq00LFH1iVsMUn3nC/RrnTp1Kjz99NOBtRUHMKx7HPioviuTiRMnwt577w0nnHBC4BN5zTXXwG677Ub1XeGsjOi/Vd+H28fG7d9b+z6z428m/JsfDIn8CVd01MZFW4MHDw7qmremY12H9SxrC2jNJfwEHfK33XbbMmt7iKw+o+qbtwAQ/oDWE3QBwoVaaIENF2489dRT0LdvX6rvCnyrgoYHXIiJz+2QIUOCaDL3338/bLrpplTfFUzbGOM1fhf+Had/99ptoGfPnrBw4cLA7zUEZ/NYQCRa8gvOzv/4xz8GAnb//fdvqut58+aVHYd/h68WZN+jHy3hb4SBN954I/BZx5+XX345+MHfqb4rD3wWcWAKhSvSv3//IJwS1XflgVGAcFLCGpLQDxInLFTflU3PGOM1ugdgP8F+jxoPxbDJeO61eEXHb5zNs4tycFEAzvBqarzOOqGwxuFrRYz/eNBBBzV9jj5TX331VdMrhbCu8fPwe/w7BF9LoMtB+D3hH4899lggVtEHEn8wDiT+4O9Yb59//nlT6BT8Hxd7yOobBRD+UH37C9YNuopMnjy56TNcxINiluq78kChgm4irAUN67tPnz5U3xXO9jHGa9RuqOHY71Hjodbbcsst9TNR9JzLL7+8eNBBBxW/+OKL4uuvv178zne+U/zXv/6VdbYIy1BZW221VfH2228vzpkzp+xn3bp1xQMPPLB4zjnnFMeNGxeE2MFQHDNnzgzOnT59ehA2DT/H788+++wgfEcYaonwHwx9FYa/Wrp0aXHXXXctXnPNNUE4Fvx/jz32aAq98tlnnxW32Wab4rPPPlv85ptvgrArp5xySsZ3QEQxfPjw4uGHHx7U2X//+9+gjh999FGq7wpkyZIlwTOLYe0mTZpUfPPNN4s777xz8amnnqL6rkAGM6Gy4o7Xr7zySqDlUNOhtkONh2OACd6LV4z7deGFFwYFg/FA//jHP2adJcISbMj4AIh+kClTphSPOuqo4rbbbhs05vfee6/s/P/85z/FH/7wh0E8OIwpRzE/8yteEey0fvaznwWd3GGHHVb86quvyo7HmK/f+973gmf/9NNPLy5YsCCDXBOmggbFDNbZbrvtVrz77rubBiyq78oDJ57HH398IET23XffYHym+q588epivEY9gH3EDjvsUBwxYkRx1apVRRMK+I8rUzJBEARBEARBJAk5jhIEQRAEQRC5gcQrQRAEQRAEkRtIvBIEQRAEQRC5gcQrQRAEQRAEkRtIvBIEQRAEQRC5gcQrQRAEQRAEkRtIvBIEQRAEQRC5gcQrQRCEI7755ptgm9uPPvoItthiCypXgiCIBCDxShAE4YjTTz8dpkyZAsOGDYN3332XypUgCCIBSLwSBEE4pk2bNtC9e3cqV4IgiAQg8UoQBOGAY445BmbOnAkjRoyAffbZp8ltYMaMGcHv//nPf4LP0Sp77bXXwrhx4+CQQw6BoUOHwimnnALLli1rSuvpp59uOhbTHTt2LNURQRBEIyReCYIgHHD33XfDxhtvDJdccknww/PQQw/BfffdB9dccw089thjcMYZZ8D5558PjzzyCIwcORKef/754Lh///vfcM8998Dll18OL774Iuywww5w7LHHwuLFi6meCIIgSLwSBEG4oUuXLtCqVSuora0NfnhOO+002HLLLeHggw+Grl27wkEHHQR77LFHIE532203mDRpUnDcww8/HFhi9957b+jXrx+cc8450Lt3b/jb3/5GVUUQBAEArakUCIIgkmfTTTdt+r1du3aBIGX/XrNmTfD7xIkT4eabb4bbbrut6fvVq1cHC8EIgiAIEq8EQRCpgFZZlpoasdfW+vXrA7cDtMaydOrUKdH8EQRB5AXyeSUIgvCI/v37w+zZs6Fv375NPw888EDgF0sQBEGQeCUIgnBGhw4dAt/VOIurTjjhBHj00UfhpZdegmnTpgUuBK+++ioMHDiQaoogCILcBgiCINxxxBFHwC233ALPPvusdRoHHnggzJs3D+66667g/0GDBsH9998fLN4iCIIgAArFYrFIBUEQBEEQBEHkAfJ5JQiCIAiCIHIDiVeCIAiCIAgiN5B4JQiCIAiCIHIDiVeCIAiCIAgiN5B4JQiCIAiCIHIDiVeCIAiCIAgiN5B4JQiCIAiCIHIDiVeCIAiCIAgiN5B4JQiCIAiCIHIDiVeCIAiCIAgiN5B4JQiCIAiCIHIDiVeCIAiCIAgC8sL/A74RPCPOK+O+AAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "n = 1000\n", "\n", "# parameters\n", "mean = 3\n", "sigma = 1\n", "rho = np.concatenate([np.linspace(-0.5, 0.9, 500), np.linspace(0.9, -0.5, 499)])\n", "\n", "# covariance matrix\n", "cov = np.diag(np.ones(n)*sigma**2.) + np.diag(np.ones(n-1)*rho*sigma**2., 1) + np.diag(np.ones(n-1)*rho*sigma**2., -1)\n", "\n", "# random variates\n", "np.random.seed(123456)\n", "obs_data = np.random.multivariate_normal([mean]*n, cov, check_valid='ignore')\n", "\n", "plt.figure(figsize=(8, 4))\n", "plt.plot(obs_data, c='r', alpha=0.7, lw=2)\n", "plt.xlim([0, 1000])\n", "plt.xlabel('time')\n", "plt.ylabel('data');" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Before we create an observation model to be used by `bayesloop`, we define a pure Python function that takes a segment of data as the first argument, and NumPy arrays with parameter grids as further arguments. Here, one data segment includes two subsequent data points `x1` and `x2`, and their known standard deviations `s1` and `s2`. The likelihood function we evaluate states the probability of observing the current data point `x2`, given the previous data point `x1`, the known standard deviations `s2`, `s1` and the parameters $\\mu$ and $\\rho$:\n", "\n", "$$P(x_2~|~x_1, s_2, s_1, \\mu, \\rho) = \\frac{P(x_2, x_1~|~s_2, s_1, \\mu, \\rho)}{P(x_1~|~s_1, \\mu)}~,$$\n", "\n", "where $P(x_2, x_1~|~s_2, s_1, \\mu, \\rho)$ denotes the [bivariate normal distribution](http://mathworld.wolfram.com/BivariateNormalDistribution.html), and $P(x_1~|~s_1, \\mu)$ is the marginal, univariate normal distribution of $x_1$. The resulting distribution is expressed as a Python function below. Note that all mathematical functions use NumPy functions, as the function needs to work with arrays as input arguments for the parameters:" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": true, "execution": { "iopub.execute_input": "2026-04-27T19:26:18.306760Z", "iopub.status.busy": "2026-04-27T19:26:18.306661Z", "iopub.status.idle": "2026-04-27T19:26:18.308845Z", "shell.execute_reply": "2026-04-27T19:26:18.308491Z" } }, "outputs": [], "source": [ "def likelihood(data, mu, rho):\n", " x2, x1, s2, s1 = data\n", " \n", " exponent = -(((x1-mu)*rho/s1)**2. - (2*rho*(x1-mu)*(x2-mu))/(s1*s2) + ((x2-mu)/s2)**2.) / (2*(1-rho**2.))\n", " norm = np.sqrt(2*np.pi)*s2*np.sqrt(1-rho**2.)\n", " \n", " like = np.exp(exponent)/norm\n", " return like" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As `bayesloop` still needs to know about the parameter boundaries and discrete values of the parameters $\\mu$ and $\\rho$, we need to create an observation model from the custom likelihood function defined above. This can be done with the `NumPy` class:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": { "collapsed": true, "execution": { "iopub.execute_input": "2026-04-27T19:26:18.309961Z", "iopub.status.busy": "2026-04-27T19:26:18.309874Z", "iopub.status.idle": "2026-04-27T19:26:18.311743Z", "shell.execute_reply": "2026-04-27T19:26:18.311357Z" } }, "outputs": [], "source": [ "L = bl.om.NumPy(likelihood, 'mu', bl.cint(0, 6, 100), 'rho', bl.oint(-1, 1, 100))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Before we can load the data into a `Study` instance, we have to format data segments in the order defined by the likelihood function:\n", "```\n", "[[x1, x0, s1, s0], \n", " [x2, x1, s2, s1],\n", " [x3, x2, s3, s2],\n", " ...]\n", "```\n", "Note that in this case, the standard deviation $\\sigma = 1$ for all time steps." ] }, { "cell_type": "code", "execution_count": 8, "metadata": { "collapsed": true, "execution": { "iopub.execute_input": "2026-04-27T19:26:18.312646Z", "iopub.status.busy": "2026-04-27T19:26:18.312578Z", "iopub.status.idle": "2026-04-27T19:26:18.314604Z", "shell.execute_reply": "2026-04-27T19:26:18.314138Z" } }, "outputs": [], "source": [ "data_segments = input_data = np.array([obs_data[1:], obs_data[:-1], [sigma]*(n-1), [sigma]*(n-1)]).T" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally, we create a new `Study` instance, load the formatted data, set the custom observation model, set a suitable transition model, and fit the model parameters:" ] }, { "cell_type": "code", "execution_count": 9, "metadata": { "execution": { "iopub.execute_input": "2026-04-27T19:26:18.315688Z", "iopub.status.busy": "2026-04-27T19:26:18.315615Z", "iopub.status.idle": "2026-04-27T19:26:18.635976Z", "shell.execute_reply": "2026-04-27T19:26:18.635550Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "+ Created new study.\n", "+ Successfully imported array.\n", "+ Observation model: likelihood. Parameter(s): ('mu', 'rho')\n", "+ Transition model: Gaussian random walk. Hyper-Parameter(s): ['d_rho']\n", "+ Started new fit:\n", " + Formatted data.\n", " + Cached observation likelihoods (76.2 MiB).\n", " + Set uniform prior with parameter boundaries.\n" ] }, { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "8da742f8622e467cbfbd3e6fce2da663", "version_major": 2, "version_minor": 0 }, "text/plain": [ " 0%| | 0/999 [00:00" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "plt.figure(figsize=(8, 4))\n", "S.plot('rho', label='mean inferred')\n", "plt.plot(rho, c='r', alpha=0.7, lw=2, label='true')\n", "plt.legend()\n", "plt.ylim([-.6, 1]);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Finally, we note that the `NumPy` observation model allows to access multiple data points at once, as we can pass arbitrary data segments to it (in the example above, each data segment contained the current and the previous data point). This also means that there is no check against looking at the data points twice, and the user has to make sure that the likelihood function at time $t$ always states the probability of **only the current** data point:\n", "\n", "$$ P(\\text{data}_{t}~|~\\{\\text{data}_{t'}\\}_{t'