diff --git a/scientific_library/tvb/simulator/monitors.py b/scientific_library/tvb/simulator/monitors.py index 71b931357c..469d34fefe 100644 --- a/scientific_library/tvb/simulator/monitors.py +++ b/scientific_library/tvb/simulator/monitors.py @@ -66,9 +66,9 @@ import tvb.datatypes.equations as equations from tvb.simulator.common import numpy_add_at from tvb.simulator.backend.ref import ReferenceBackend -from tvb.basic.neotraits.api import HasTraits, TVBEnum, Attr, NArray, Float, EnumAttr, narray_describe +from tvb.basic.neotraits.api import HasTraits, Range, TVBEnum, Attr, NArray, Float, EnumAttr, narray_describe from .backend import ReferenceBackend - +import tvb.simulator.integrators as integrators_module class Monitor(HasTraits): """ @@ -394,7 +394,6 @@ def sample(self, step, state): time = (step - self.istep / 2.0) * self.dt return [time, avg_stock] - class AfferentCoupling(RawVoi): """ A monitor that records the variables_of_interest from node_coupling data from a tvb simulation @@ -1005,6 +1004,226 @@ def create_time_series(self, connectivity=None, surface=None, sample_period=self.period, title='Regions ' + self.__class__.__name__) +class BoldModels(TVBEnum): + LINEAR = "linear" + NONLINEAR = "nonlinear" + +class BoldBalloonWindkessel(Monitor): + """ + Balloon-Windkessel model for hemodynamic response. + + The haemodynamic model parameters based on constants for a 1.5 T scanner. + """ + + bold_model = EnumAttr( + default=BoldModels.NONLINEAR, + label="Select BOLD model equations", + doc="""Select the set of equations for the BOLD model.""") + + period = Float( + label="Sampling period (ms)", + default=1000.0, + doc="""Repetition time (TR).""") + + integrator = Attr( + field_type=integrators_module.Integrator, + label="Integration scheme", + default=integrators_module.HeunDeterministic(dt=10.), + required=True, + doc=""" A tvb.simulator.Integrator object which is an integration + scheme with supporting attributes such as integration step size and + noise specification for stochastic methods. It is used to compute the + time courses of the balloon model state variables.""") + + RBM = Attr( + field_type=bool, + label="Revised BOLD Model", + default=True, + required=True, + doc="""Select classical vs revised BOLD model (CBM or RBM). + Coefficients k1, k2 and k3 will be derived accordingly.""") + + include_svars = Attr( + field_type=bool, + label="Include physiological state variables.", + default=False, + required=True, + doc="""If true, the state variables of the hemodynamic model are + returned (s, f, v, q). Otherwise only the BOLD signal is returned. """) + + + tau_s = Float( + label=r":math:`\tau_s`", + default=1.54, + required=True, + doc="""Balloon model parameter. Time of signal decay (s)""") + + tau_f = Float( + label=r":math:`\tau_f`", + default=1.44, + required=True, + doc=""" Balloon model parameter. Time of flow-dependent elimination or + feedback regulation (s). The average time blood take to traverse the + venous compartment. It is the ratio of resting blood volume (V0) to + resting blood flow (F0).""") + + tau_o = Float( + label=r":math:`\tau_o`", + default=0.98, + required=True, + doc=""" + Balloon model parameter. Haemodynamic transit time (s). The average + time blood take to traverse the venous compartment. It is the ratio + of resting blood volume (V0) to resting blood flow (F0).""") + + alpha = Float( + label=r":math:`\tau_f`", + default=0.32, + required=True, + doc="""Balloon model parameter. Stiffness parameter. Grubb's exponent.""") + + TE = Float( + label=":math:`TE`", + default=0.04, + required=True, + doc="""BOLD parameter. Echo Time""") + + V0 = Float( + label=":math:`V_0`", + default=4.0, + required=True, + doc="""BOLD parameter. Resting blood volume fraction.""") + + E0 = Float( + label=":math:`E_0`", + default=0.4, + required=True, + doc="""BOLD parameter. Resting oxygen extraction fraction.""") + + epsilon = NArray( + label=":math:`\epsilon`", + default=numpy.array([0.5]), + domain=Range(lo=0.5, hi=2.0, step=0.25), + required=True, + doc=""" BOLD parameter. Ratio of intra- and extravascular signals. In principle this + parameter could be derived from empirical data and spatialized.""") + + nu_0 = Float( + label=r":math:`\nu_0`", + default=40.3, + required=True, + doc="""BOLD parameter. Frequency offset at the outer surface of magnetized vessels (Hz).""") + + r_0 = Float( + label=":math:`r_0`", + default=25., + required=True, + doc=""" BOLD parameter. Slope r0 of intravascular relaxation rate (Hz). Only used for + ``revised`` coefficients. """) + + def compute_derived_parameters(self): + """ + Compute derived parameters :math:`k_1`, :math:`k_2` and :math:`k_3`. + """ + + if not self.RBM: + """ + Classical BOLD Model Coefficients [Obata2004]_ + Page 389 in [Stephan2007]_, Eq. (3) + """ + k1 = 7. * self.E0 + k2 = 2. * self.E0 + k3 = 1. - self.epsilon + else: + """ + Revised BOLD Model Coefficients. + Generalized BOLD signal model. + Page 400 in [Stephan2007]_, Eq. (12) + """ + k1 = 4.3 * self.nu_0 * self.E0 * self.TE + k2 = self.epsilon * self.r_0 * self.E0 * self.TE + k3 = 1 - self.epsilon + + return numpy.array([k1, k2, k3]) + + def balloon_dfun(self, state_variables, neural_input, local_coupling=0.0): + s, f, v, q = state_variables + + x = neural_input[0, :] + + ds = x - (1. / self.tau_s) * s - (1. / self.tau_f) * (f - 1.) + df = s + dv = (1. / self.tau_o) * (f - v ** (1. / self.alpha)) + dq = (1. / self.tau_o) * ((f * (1. - (1. - self.E0) ** (1. / f)) / self.E0) - + (v ** (1. / self.alpha)) * (q / v)) + + return numpy.array([ds, df, dv, dq]) + + + + def _bold(self, state): + s, f, v, q = state + + # BOLD models + if self.bold_model == "nonlinear": + """ + Non-linear BOLD model equations. + Page 391. Eq. (13) top in [Stephan2007]_ + """ + y_bold = numpy.array(self.V0 * (self.k1 * (1. - q) + self.k2 * (1. - q / v) + self.k3 * (1. - v))) + y_b = y_bold[numpy.newaxis, :, :] + self.log.debug("Max value: %s" % str(y_b.max())) + + else: + """ + Linear BOLD model equations. + Page 391. Eq. (13) bottom in [Stephan2007]_ + """ + y_bold = numpy.array(self.V0 * ((self.k1 + self.k2) * (1. - q) + (self.k3 - self.k2) * (1. - v))) + y_b = y_bold[numpy.newaxis, :, :] + return y_b + + + def _config_time(self,simulator): + super()._config_time(simulator) + self.k1, self.k2, self.k3 = self.compute_derived_parameters() + + self._state = numpy.ones( (4, simulator.number_of_nodes, simulator.model.number_of_modes) ) + self._state[0,:] = 0. + + # or allow providing bold dt as multiple of model dts... + assert self.integrator.dt >= simulator.integrator.dt + self._dt_istep = ReferenceBackend.iround(self.integrator.dt/simulator.integrator.dt) + self.integrator.dt = self._dt_istep * simulator.integrator.dt / 1000. + self.integrator.configure() + + # _tavg_stock covers the time between the bold timesteps + self._tavg_stock = numpy.zeros( ( self._dt_istep, self.voi.shape[0], + simulator.number_of_nodes, + simulator.model.number_of_modes)) + + def sample(self, step, state): + # push to timeaverage + self._tavg_stock[((step % self._dt_istep) - 1), :] = state[self.voi] + # if time steps over bold dt + if step % self._dt_istep == 0: + neural_activity = numpy.mean(self._tavg_stock, axis=0) + self._state = self.integrator.scheme( + self._state, + self.balloon_dfun, + neural_activity, + local_coupling=0.0, stimulus=0.0) + # if bold period, yeet current state + if step % self.istep == 0: + time = step * self.dt + bold = self._bold(self._state) + if self.include_svars: + return [time, numpy.concatenate([bold, self._state])] + else: + return [time, bold] + + + class ProgressLogger(Monitor): """Logs progress of simulation; only for use in console scripts.""" diff --git a/tvb_documentation/demos/bold_balloon_monitor.ipynb b/tvb_documentation/demos/bold_balloon_monitor.ipynb new file mode 100644 index 0000000000..2714ffc19c --- /dev/null +++ b/tvb_documentation/demos/bold_balloon_monitor.ipynb @@ -0,0 +1,370 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "id": "4cc56589-d851-433f-8d83-9a13930b83c0", + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib inline" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "34bc3e01-2790-4f92-b3d0-87cd1ae271ac", + "metadata": {}, + "outputs": [], + "source": [ + "%load_ext autoreload" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "aafed4b5-4465-4b16-9437-32402e18f520", + "metadata": {}, + "outputs": [], + "source": [ + "%autoreload 2" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "969fbcbf-164b-4843-8c3f-799473454f3b", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import matplotlib.pylab as plt" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "07a4d85a-1a44-4d6a-a038-2e75c7e85116", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "from tvb.simulator.lab import *" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "973afb69-24c8-41ce-8c12-7f08427b5761", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "from tvb.simulator.backend.nb_mpr import NbMPRBackend" + ] + }, + { + "cell_type": "markdown", + "id": "e7110756-486f-422c-bd5a-78a148afbc8a", + "metadata": {}, + "source": [ + "We start with building a small toy network to demonstrate the monitor." + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "7ac40273-2686-43e1-a1ff-f389ad433f23", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "N = 2\n", + "conn = connectivity.Connectivity()\n", + "conn.motif_all_to_all(number_of_regions=N)\n", + "conn.centres_spherical(number_of_regions=N)\n", + "conn.speed = np.r_[np.inf]" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "217226a3-f3d4-4bce-a770-0b417b09d263", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "conn.configure()" + ] + }, + { + "cell_type": "markdown", + "id": "cafb8502-9623-4ac8-8a99-df5cc8f68a6a", + "metadata": {}, + "source": [ + "Stimulus to make one of the nodes active for a brief period of time. " + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "e9433a2a-d940-4684-bf6e-8f15492bc062", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "weighting = patterns.StimuliRegion.get_default_weights(N)\n", + "weighting[0] = 4.\n", + "\n", + "# temporal profile\n", + "eqn_t = equations.Gaussian()\n", + "eqn_t.parameters[\"midpoint\"] = 11000.0\n", + "eqn_t.parameters[\"sigma\"] = 100.0\n", + "\n", + "stimulus = patterns.StimuliRegion(temporal=eqn_t,\n", + " connectivity=conn,\n", + " weight=weighting)\n" + ] + }, + { + "cell_type": "markdown", + "id": "0c3430e4-3032-48c9-9aa6-689f609fc295", + "metadata": {}, + "source": [ + "We use the Montbrio-Pazo-Roxin model with the Numba backend to generate the time series quickly." + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "id": "1aa7d93b-bcef-4504-9b01-31eca1f8f061", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "sim = simulator.Simulator(\n", + " model = models.MontbrioPazoRoxin(),\n", + " connectivity = conn,\n", + " integrator = integrators.HeunStochastic(\n", + " dt = 0.01, \n", + " noise = noise.Additive(nsig=np.r_[0.01,0.02])\n", + " ),\n", + " initial_conditions=np.zeros( (1,2,N,1) ),\n", + " stimulus=stimulus,\n", + " monitors=[monitors.Raw()]\n", + ").configure()" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "id": "09b5fbb8-5e4f-4ff9-81a3-7aa4013f4ca9", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "backend = NbMPRBackend()" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "id": "94941692-ba67-45d0-9fb2-b4d770e69efd", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "%%capture\n", + "(raw_t, raw_d), = backend.run_sim(sim, simulation_length=30000)" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "id": "e72fe8a7-ba49-40ce-872f-2c45ebf1eb4b", + "metadata": { + "tags": [] + }, + "outputs": [ + { + "data": { + "text/plain": [ + "((3000000,), (3000000, 2, 2, 1))" + ] + }, + "execution_count": 13, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "raw_t.shape, raw_d.shape" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "id": "0b71f5ad-bb6d-40fb-a38a-ffc4b164154b", + "metadata": { + "tags": [] + }, + "outputs": [ + { + "data": { + "text/plain": [ + "[,\n", + " ]" + ] + }, + "execution_count": 14, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAfUElEQVR4nO3deXxU5b3H8c8vISG2KmqhimJFra3VuoDUYrVWW6uArVarvWqt2tZyrXq19tYWtVUUveBaV0QUFGSVTaiAsgmIyBIg7AIBQtgTIAQC2fPcP84hmSSTlVmY4/f9euWVM+c8c87vmTPznTPPmcWcc4iISOJLincBIiISGQp0EZGAUKCLiASEAl1EJCAU6CIiAdEiXhtu3bq1a9++fbw2LyKSkBYtWrTLOdcm3LK4BXr79u1JT0+P1+ZFRBKSmW2qa5mGXEREAkKBLiISEAp0EZGAUKCLiASEAl1EJCAU6CIiAaFAFxEJCAW6xMz4jK3sLyqNdxkigaVAl5hYuS2fB0dk0GPM8niXIhJYCnSJicKScgB27CuKcyUiwaVAl5i7oe/nDJizMd5liASOAl1ibkn2Xnp9tCreZYgEjgJdYkq/YSsSPQp0iQmzeFcgEnwKdImpxdl7412CSGAp0EVEAkKBLiISEAp0EZGAUKBLjOisqEi0KdBFRAJCgS4iEhAKdBGRgFCgS0zog0Ui0adAFxEJCAW6iEhANBjoZnaqmX1qZqvMbKWZPRimjZnZq2aWaWbLzKxjdMoVEZG6tGhEmzLgf51zi83sGGCRmU11zoV+/2lX4Cz/74fAm/5/ERGJkQaP0J1z251zi/3p/cBq4JQaza4HBjvPPOA4M2sb8WpFRKROTRpDN7P2QAdgfo1FpwCbQy5voXboY2bdzSzdzNJzc3ObWKqIiNSn0YFuZkcDY4C/OOf2NWdjzrn+zrlOzrlObdq0ac4qRESkDo0KdDNLwQvzoc65sWGabAVODbnczp8nIiIx0ph3uRgwAFjtnHupjmYTgDv8d7t0BvKdc9sjWKeIiDSgMe9yuRT4HbDczDL8eY8C3wJwzvUDJgHdgEzgIPD7iFcqCU0fFBWJvgYD3Tk3hwYej8775d/7IlWUiIg0nT4pKiISEAp0EZGAUKBLTJi+blEk6hToIiIBoUAXEQkIBbqISEAo0EVEAkKBLjGhU6Ii0adAFxEJCAW6iEhAKNBFRAJCgS4xoc8ViUSfAl1EJCAU6BITzsW7ApHgU6CLiASEAl1iInd/cbxLEAk8BbrERFmFxlxEok2BLiISEAp0EZGAUKCLiASEAl1iQh8sEok+BbqISEAo0EVEAkKBLiISEAp0iZui0vJ4lyASKAp0iYlw50Q79poa8zpEgkyBLnFzsERH6CKRpEAXEQkIBbqISEAo0CUmTJ8sEok6BbqISEAo0EVEAkKBLiISEAp0EZGAUKCLiASEAl1EJCAaDHQzG2hmOWa2oo7lV5hZvpll+H+PR75MERFpSGOO0N8DujTQ5jPn3IX+31OHX5Z8VezcVxTvEkQCo8FAd87NBvbEoBYJsJOOTQs7f/X2fTGuRCS4IjWGfomZLTWzyWZ2bl2NzKy7maWbWXpubm6ENi2JIDlJnxQVibZIBPpi4DTn3AXAa8CHdTV0zvV3znVyznVq06ZNBDYtic7FuwCRADnsQHfO7XPOFfjTk4AUM2t92JXJV4JzinSRSDnsQDezk8z/5iUzu9hf5+7DXa98NSjPRSKnRUMNzGw4cAXQ2sy2AE8AKQDOuX7ATcCfzawMKARucTrskkaq0D1FJGIaDHTn3K0NLH8deD1iFclXyp8Gp5PV59p4lyESCPqkqIhIQCjQRUQCQoEuIhIQCnQRkYBQoIuIBIQCXUQkIBToIiIBoUCXmHD61haRqFOgi4gEhAJdRCQgFOgiIgGhQBcRCQgFuohIQCjQRUQCQoEuIhIQCnQRkYBQoIuIBIQCXUQkIBToIiIBoUAXEQkIBbqISEAo0CWu2rZKi3cJIoGhQBcRCQgFusSEq+Pr0OuaLyJNp0AXEQkIBbqISEAo0EVEAkKBLjFR5xi6fmtUJGIU6BITdQW3ToqKRI4CXWJCwS0SfQp0iSvlvEjkKNAlJuoK7tz9xTGtQyTIFOgiIgGhQJeYcBpEF4k6BbqISEA0GOhmNtDMcsxsRR3LzcxeNbNMM1tmZh0jX6YkOh2fi0RfY47Q3wO61LO8K3CW/9cdePPwyxIRkaZqMNCdc7OBPfU0uR4Y7DzzgOPMrG2kChQRkcaJxBj6KcDmkMtb/Hm1mFl3M0s3s/Tc3NwIbFoShc6JikRfTE+KOuf6O+c6Oec6tWnTJpabFhEJvEgE+lbg1JDL7fx5IiF0iC4SbZEI9AnAHf67XToD+c657RFYr4iINEGLhhqY2XDgCqC1mW0BngBSAJxz/YBJQDcgEzgI/D5axUri0hi6SPQ1GOjOuVsbWO6A+yJWkYiINIs+KSoxoQN0kehToIuIBIQCXWJCY+gi0adAl5jQty2KRJ8CXWJCcS4SfQp0EZGAUKBLTKQk664mEm16lElMfOPrqfEuQSTwFOgiIgGhQBcRCQgFuohIQCjQRUQCQoEuIhIQCnQRkYBQoIuIBIQCXWLqiV+eE+8SRAJLgS4xddzXUhjRvXO8yxAJJAW6xJzFuwCRgFKgS8yd+c2j412CSCAp0CXmWh/dMt4liASSAl1EJCAU6CIiAaFAl5jQLxaJRJ8CXWLK9B4XkahRoIuIBIQCXeLi7JOOiXcJIoGjQBcRCQgFusSFfjRaJPL0qJK4OKftsfEuQSRwFOgSF0m654lEnB5WEid6+6JIpCnQJS6SlOciEadAl7gwBbpIxCnQJS6SlOgiEadAl7hQnItEXqMC3cy6mNkaM8s0sx5hlt9lZrlmluH/3R35UiVITEfoIhHXoqEGZpYMvAH8HNgCLDSzCc65VTWajnTO3R+FGiWANOQiEnmNOUK/GMh0zm1wzpUAI4Dro1uWBN2xRzV4LCEiTdSYQD8F2BxyeYs/r6Zfm9kyMxttZqdGpDoJrN900l1EJNIidVL0P0B759z5wFRgULhGZtbdzNLNLD03NzdCm5ZE1EJvRBeJuMYE+lYg9HCqnT+vknNut3Ou2L/4DnBRuBU55/o75zo55zq1adOmOfVKUCjPRSKuMYG+EDjLzE43s1TgFmBCaAMzaxty8TpgdeRKlCDSLxeJRF6DZ6acc2Vmdj/wCZAMDHTOrTSzp4B059wE4AEzuw4oA/YAd0WxZklAzlX/VVG9yUUk8hr1VgPn3CRgUo15j4dMPwI8EtnSJIgOBbnetigSefqkqMSF4lwk8hToEhc6QBeJPAW6xIVOiopEngJd4kN5LvEw42lY/Z94VxE1+vy1xIU+VyRxMft573/P/PjWESXBO0KvqICK8nhXUc22vYW8OXN9rbfufZXp2xZFIi94gd73h9CrdbyrqObuQek8+/GXZO85GO9SjhiK88Qy+Iss7h6UHrH1bc8vZLMeDxGX0IH+Qfpm+s7MrD5z11pwFfEpqA4HSsoACPIB+sGSMi5+ZhpzM3eFXb4upwCA2Wu95U0+QK8oh+KCxrXNng8HquooKC5jf1FpEzcooR4fv5Jpq3dGbH2X9J7Bj5/71LuwYwWUa/9EQkIH+t9HL+O5j9dEbwP7d0Dh3sNeTcIF+e71kPNlk66ydmcBOfuLefbj8NdbsHEPAB+v2A40/MGigyVlDJufXTVMNfZP0Nv7ks9lW/ayafeBuq888GoYcHXlxfN7fsJ5PadUa5Kzv4jRi7bUW0O85ReW0r7HRD5ati3epRy2j1fsYPX2fbXmvzJyEvS7FKb1bPpKV4z1/ppjx4rGt107BRYOaN52YiyhAz3qXvwuvHxe7fk9W1E8+bFGr8bhhVJdGdbUsfVtewub1L7JXuvoDV2FyDtQwqTl2+u8yqGuVevJ5oWwYWazSnhm4moeHbecWWv9b+VcMaZy2XWvf85Pnp9Z+SRRqWcrmNnHm96zHj68D6Y+TkWYm/eP76Xzt1FLyd1fXHthiMXZeQ22OSxrJtd5dLoh13tFcv+wJWzcVf0JrLS8gopwHavHta9+xqC5Wc0q83DdM2QRXV/5rNb8ORn+1z5tXdT0lY7+vffXHP0ubbBJRYVj8BdZMOxmmPjX2g3KSiBvU/0rWTUeNs1tXo3NEKhA35FfFPmVFlc/qqiY3guAlvNfh4JcWDK0evttGV6wrP+01qrCvfd6+uqdnP7IJL7cUbWdu95dwH+Whj8q+3DJVn7UZwbzNuz2ZmQMh7mvVTU4uMfbfuZ0eg2byvRlG2uvpOdxMOGBWrO37S3kQHFZ2O3eM2QR9w5dzM594W/jMyfezM3JM6kIfXIacBUMvh4yhvPjbQOrtW9oyCXvYAkAB4rL2Z4f/gns0XHL/f60gvdv9KZn9q5qkDEEPn8l7HVz9nv9KK8rFLcuhr2bubHvXLq9+pn3Sm3pSC+Ac2p/91xRaTmPjltO3oESKCuGBW9DRTkVFY6HRmbw76lra11n7ecfwvBbYPI/vJP5ocpKqg0dXvnCzGqLz3psMg+PXlZ5efmW/AbHpFdu28cTE1YCsHHXARZn59Xb/pD2PSY2ql1N2/YWkplTQHFZjTcprJnML5K+4Nbk6VyW7B8pZ38Bc16G8jLYt927HzfDjvwiLuo1lfW5jRyeq8f4pVt5fPzKqhl5myBzetXR+vj74JXz6TV2IQ+PWuoN8814pvq+/OAOeLfrYdfSWIkb6OVlXJq0nFmpf2Fv/18y48uddO49vXJxh6emMHd9+PHcWpYMgZ6tyN6czcvT1tK+x0T21THmmvTZC1UXRtwG4+/l/U/msiXPfzAdejZe+0lls5z8Qk6zHeQW1D7SG7vE+ybizNVLK8dmZq7J5X+GL6HXR6tqPUgXZ+fxq6Q5nD3mKm/Gh/fAlH960xUVMOoub3rWc/xr7U0cP/qmqisX7oW5rwMOFg9id0Exw+ZnezdBdh4/6jODm/t9Ebbfh+ooCwnAgyVlXoABR+9YwPMp/cMPL314D1dsexujggMl3oO75pNbSVkFYxdvqXy1cmh5hXNc0ntG2JqqvbJZPz1sG4CT2M0phP/+/UOvntj0BSwbVbXg7Svh5e8DkLu/mNX9bodx3b0A7tu5ql1BLmxfyodLtjJsfjY9xi4jd1JvmPQ3eOoE/jlqAeOWbOWV6etg/04Y+htvPxTmMWLSNG8d6QNg1rPe9O71kL8Fnm7DKfOfqrNPAGMWVw0Z/fL1OVVj0jXt3wHZ86r67BxXvjCTG/vO9YYs6jiCXLNjP89MrP5Lk90HVz8xeqC4jG8/OokpK3fAR3/1AgzYXVDMj/rM4K1XnqTvkJEcw0HS8O//w2/h9dTX6J0ygAdbhAyZTHuCvHmD4KWzcS99jz6Tv2RfUSm5+4s5WBL+QAOAz1/1nkSBj5ZtY/eBEobOy667fRil5RWUlHlBvHnPQVZszWdfYRk/sJAhxFfOhyE3Vh2tr/Me46MWbGDUoi0w/n6Y/RxkDK25eigrxu3f2eRXVU2VmO9Dz5oDw29jaKr/XtJtOTw3YgpwQmWTvIOl3Pb2fLL6XFv7+iN+C6UH4XfjvMsfPQTAA33HkeG+DcCu/cUc6zeftTaXbx7Tku998XD19RTsAOCtmWsYvLKUqX/9SchCR9mO1ezJXsnvmUqPliO4pl8Jn/T+s/cyLDmV1TsK+NrKBVxo7fjFrMdhFuy9oyq8BszZSHrWHsbff1nlvCQzXk7tCzWHkN+6HArzYK9/R97sPYA7JmVyXo9RfPRwN5a8fCu/Sq568D40fCGz1+dz8enHc0Nfb/6q7fsgrfZNVu4cHW0tqXvawXEXAHD1v2ezJa+w2m186P66bW8hJ9dYx/DUZ5hY/kPYfxFWEnKkv28bw6YtpucC46iUZJJyVnD8qom04Epq3f0HdiErzXvS+XP5C8AVtYutYV7a/wAw+IsrWLwpj4VZeSS5ch5MHsO67LNpe96Z8G4XAM4elszE/72aM2usoyRva7XDn7GLt3BUSjJdP74cDuTwtTOe4CKDT1bCJS1Wcpf/yLpyxSPckZpDl5JnGfjsA/yhxSfs+vxdWs/pyeMpVevLXzKWYy+5F3utY+W81l8O5Sw7h7+3GMniirMguzWcfCFUVIVb+x4TGX+fN3yQTDkVT7Uh6aRzef/8QWzbW8hV3zuRi0ZeCgd38X17mhXuDLL738q41HXcUPJU1ZBFyPuyB3+RxeTlO8jL2cxbJY8yxXqwyZ0EwJRVOyF3LcXbVpD+9ctpdVQKZRWOGZNHc3WBd+S6ODuPG/vO5buWzfMp/SELHkqDne444MZ699XLk5byZApYWRFDZi2noLiUIfOySaKCqQ905syTw/yOwtR/sW3HdrZf9DDtcj/jqqRNPL7oNvjJl3Bs29rtw/hRnxnsLihmQ+9r6fNCb6ZXdKCIlmSl1fOk6h9Q/CJ5PlPKO8Hayd78CffDub9i3Kp8bjjUdsRtWOY0Ti8aFj6TIiTxAr24AN6rfYP8u+JZuvJstXmnkMus+QvpNTufp3/5HUqSjqKkrIKrvvwIgGVr13PeB5di5SW11rc9v4gz/OkThlzF/5XdwvupI6q1cYV7MeDe5PFs2H0y8JPKsYSNuwo4vV9nvglcnHShV4/t4uMVO+gy2juK+R7wfEq1VdJvQH+8r5SvqiO/sJRWR6XArkzaHag6mVM+rRfJlQ2X1nGDwfK0u+n05mhepPrL0MFbr2VMyo9xBd8Ne70563Zx2VAv2p4s/wFdWi6EwT3h5I4sOv8JLtk3mRl0gIxhVbeJc+Rlr+R3b37O9JbV19c5aTWdk1bDi+/hdXsoYPDS97gLeIn+3Ds0nY1pt3NNCjyd8i5TJ1/NiSG3B9lVryA6l6bDysa/RTX05fNNaQt4KGUMo0bmwnlVnxz8d0pflm25uFqgf9eyuSBpQ7V1/fUD7/bOSssB4LoNT3JdS2hfNAwLeRr6eXLV2LDzX3W8+Wkm/6qx31vtWwt9vlVtXlJFKc+nvMWFSeu99QwMvf9V3eY9+g4jK+0RHi+9k6SKEti2hEEbP2a3O5Y3Zx5LVpr3SvWjlv/kvKJ3OG37ZE5LglSqXoXuLyplxpc5XH/hKZW301tHvctpSTnMavlXvlM0iBJ/r/HGD2gJHFNxOht/7X0Ja5+CqnNKye/8lBR68knLHtX6c6LthXH3UJ/TrOqdNCvS7ubytR8C8FJKX87sfzvjrl9J1++35Z0Zqwj9RfqZS75kyeLePJ/Sny6p/syXzmbWBc9x6XV/4rFxKxi7ZAvrQm53N/c1OOoErMNvK8+T7F/7GW+kvgp4+7JO25ZAqTcU+H8pA/hNco1XR73b8X5xT2449BjI9F6NJVHBkHmbuL3zafXeDs2VeIG+YnTY2amU8dvkaZWXj+Egn6c9CJOhgzuKY0cUVu6gLP8I9LX3hvB2atX4bIekdWSUe0fov31nfmW785KyeD+1T61tmj++flsLb2f26/cdfviNQjoAp68fUtkuyX+Ap1BOjyEz6RLmCDjUOZZFEamMSn2Svx24h+HPvMs9Pd+F1y/i7pB2yXNeqHMdNZUV7MZSar/c+3XyZzD4Im5K/m9eSHmLR0v/WLns9cFDucx/xuiSvLDqStsW02bL3TyfkgspUDS1XeVBvbkKjh/4o1phHs6VSRlsdlVHXMvSuvNG2XXV2vy8aAo/T5tS86oA3Fk6AkaNCLssnHuTx9O//FoeaTGcn+A9Md7cYna18dquyQs5Z+R8bvA71Ia9tYIJ4KdJi8nyj1pDdbS1XJS0rtb89rad75r36un7SWHOa9ShtYX/RGNW2m1eHcUv0C15PgDdkhdULp/W8u8AXFn8YrXrtQwJ8cohEKDXiy9xQeF8Vn7z3cp517jPK6cfaTGMJ8vuIPQTBOcnbaTN7LvISptfbRsXJG3g0qQ63kWydHj4+b4/tPi42uVDn9049Mpy0egXeGjkzyv7X8V5rwZqOHHJa7xz/M8Ymb7ZewILCXTzhyqnlJ7HWbaFda4dWRmfcuhtEP9VM6RD9b+i2sW2VnvMf2zLnrXmbUi7nVsnPMbtnf9e97oPg8Xr04udOnVy6elN/6BC6cL3SJn4YLO2Oab8MnqUdmddmneE3K/sF9zT4qNqba4tfoY2ls+sivPZmHZ7s7ZTn2KXQktr+ntuy52RbIn2/sem2+ZO4OQwD45Ieb70Nzyc8kG9bb4oP4dLklfV2yZWdrjjOckad/IynKdLf8s/U6rGdIeVXVl5APKP0j/xbMrb1dpflTqUdgUZ/HfyxLC3wTXFfcI+wUXL1PKOrHbf4oEWH1bOe6DkPl5NfaNJ6/l18ROMaflkvW1uKH6ScS2faE6Z5Ljj+KbtbXT7zQ9s59QTvtasbZnZIudcp7DLEi3Qpw9/iZ+tqX/H1KfMJdHCjqwPHknszCi/kJ8mZ8S7jCPWnSX/YFDqsw03lMOysns2557cqlnXrS/QE+5dLqUVh1eywvyrTWFeP4V5bDw6dFZU1ptwgb5+89Z4lyAicljGH/hdVNabcIG+t6Cej3yLiHyFJVygV+h7+kREwkq4QE+8b7oSEYmNhAv0m5Nnx7sEEZEjUsIF+nfsyP7KUxGReEm4QN9d+Q0rIiISKuECfYNr3JftiIh81SRcoIuISHgJF+hOb1sUEQkr4QK9winQRUTCSbhA1xG6iEh4CRjoIiISTsIFekXilSwiEhMJl44achERCU+BLiISEAkX6KZRdBGRsBIu0EVEJLxGBbqZdTGzNWaWaWa1fiHWzFqa2Uh/+Xwzax/xSn0achERCa/BQDezZOANoCtwDnCrmZ1To9kfgTzn3LeBfwNR+2FCDbiIiITXmCP0i4FM59wG51wJMAK4vkab64FB/vRo4GdmFpVD6VJaRGO1IiIJrzGBfgqwOeTyFn9e2DbOuTIgH/hGzRWZWXczSzez9Nzc3GYV3KbL35t1PRGRI8W8DtEZxIjp4a5zrj/QH6BTp07NGj3pcGkXuDQ/onWJiMRS5yittzFH6FuBU0Mut/PnhW1jZi2AVsDuSBQoIiKN05hAXwicZWanm1kqcAswoUabCcCd/vRNwAzn9GvOIiKx1OCQi3OuzMzuBz4BkoGBzrmVZvYUkO6cmwAMAN43s0xgD17oi4hIDDVqDN05NwmYVGPe4yHTRcDNkS1NRESaQp8UFREJCAW6iEhAKNBFRAJCgS4iEhAWr3cXmlkusKmZV28N7IpgOfGkvhyZgtKXoPQD1JdDTnPOtQm3IG6BfjjMLN051ynedUSC+nJkCkpfgtIPUF8aQ0MuIiIBoUAXEQmIRA30/vEuIILUlyNTUPoSlH6A+tKghBxDFxGR2hL1CF1ERGpQoIuIBETCBXpDP1h9JDCzLDNbbmYZZpbuzzvBzKaa2Tr///H+fDOzV/3+LDOzjiHrudNvv87M7qxrexGufaCZ5ZjZipB5EavdzC7yb5tM/7pR+9XvOvrS08y2+vsmw8y6hSx7xK9rjZldEzI/7H3O/0rp+f78kf7XS0ejH6ea2admtsrMVprZg/78hNsv9fQlEfdLmpktMLOlfl+erG/7ZtbSv5zpL2/f3D7WyTmXMH94X9+7HjgDSAWWAufEu64wdWYBrWvMew7o4U/3AJ71p7sBkwHD+yGT+f78E4AN/v/j/enjY1D75UBHYEU0agcW+G3Nv27XGPelJ/C3MG3P8e9PLYHT/ftZcn33OeAD4BZ/uh/w5yj1oy3Q0Z8+Bljr15tw+6WeviTifjHgaH86BZjv34Zhtw/cC/Tzp28BRja3j3X9JdoRemN+sPpIFfpD2oOAX4XMH+w884DjzKwtcA0w1Tm3xzmXB0wFukS7SOfcbLzvtI947f6yY51z85x3Tx4csq5Y9aUu1wMjnHPFzrmNQCbe/S3sfc4/gv0p3o+iQ/XbJaKcc9udc4v96f3Aarzf8U24/VJPX+pyJO8X55wr8C+m+H+unu2H7q/RwM/8epvUx/pqSrRAb8wPVh8JHDDFzBaZWXd/3onOue3+9A7gRH+6rj4dSX2NVO2n+NM158fa/f5QxMBDwxQ0vS/fAPY670fRQ+dHlf8yvQPe0WBC75cafYEE3C9mlmxmGUAO3hPk+nq2X1mzvzzfrzdiGZBogZ4oLnPOdQS6AveZ2eWhC/2joIR8v2gi1+57EzgTuBDYDrwY12qawMyOBsYAf3HO7Qtdlmj7JUxfEnK/OOfKnXMX4v3W8sXA2fGsJ9ECvTE/WB13zrmt/v8cYBzejt7pv7TF/5/jN6+rT0dSXyNV+1Z/uub8mHHO7fQfhBXA23j7Bprel914QxktasyPCjNLwQvAoc65sf7shNwv4fqSqPvlEOfcXuBT4JJ6tl9Zs7+8lV9v5DIgGicLovWH95N5G/BOHBw6SXBuvOuqUePXgWNCpufijX0/T/UTWM/509dS/QTWAn/+CcBGvJNXx/vTJ8SoD+2pfiIxYrVT++Rbtxj3pW3I9EN4Y5cA51L9xNQGvJNSdd7ngFFUP/l1b5T6YHjj2i/XmJ9w+6WeviTifmkDHOdPHwV8Bvyiru0D91H9pOgHze1jnTVF88EUpRuxG96Z8fXAY/GuJ0x9Z/g3/FJg5aEa8cbKpgPrgGkhDyQD3vD7sxzoFLKuP+CdIMkEfh+j+ofjveQtxRuz+2Mkawc6ASv867yO/2nlGPblfb/WZcCEGkHymF/XGkLe5VHXfc7f1wv8Po4CWkapH5fhDacsAzL8v26JuF/q6Usi7pfzgSV+zSuAx+vbPpDmX870l5/R3D7W9aeP/ouIBESijaGLiEgdFOgiIgGhQBcRCQgFuohIQCjQRUQCQoEuIhIQCnQRkYD4f5uveUNgGPvxAAAAAElFTkSuQmCC\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "plt.plot(raw_t, raw_d[:, 0,:,0])" + ] + }, + { + "cell_type": "markdown", + "id": "32541a1b-4d23-4fb4-a36f-18f4e9aa8837", + "metadata": {}, + "source": [ + "The following function will allow us apply the Bold monitor ex-post." + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "id": "9b8e8a83-59d0-4479-8f79-77b72ca64c95", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "def replay_monitor(raw_d, sim):\n", + " \"\"\" Assumes that the raw_d was produced in consistency with the simulator used.\"\"\"\n", + " assert len(sim.monitors) == 1, \"Only one monitor at a time can be replayed.\" \n", + " res = [out for out in [ sim.monitors[0].record(step, sim.model.observe(state)) for step, state in enumerate(raw_d)] if out is not None]\n", + " b_t, b_d = zip(*res)\n", + " b_t, b_d = np.array(b_t), np.array(b_d)\n", + " return b_t, b_d" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "id": "4ccb1778-4fe2-4042-b23b-0fe7da08db38", + "metadata": { + "tags": [] + }, + "outputs": [ + { + "name": "stderr", + "output_type": "stream", + "text": [ + "/home/izaak/local_repos/megaloceros/tvb/scientific_library/tvb/simulator/monitors.py:1147: VisibleDeprecationWarning: Creating an ndarray from ragged nested sequences (which is a list-or-tuple of lists-or-tuples-or ndarrays with different lengths or shapes) is deprecated. If you meant to do this, you must specify 'dtype=object' when creating the ndarray.\n", + " return numpy.array([k1, k2, k3])\n" + ] + } + ], + "source": [ + "sim = simulator.Simulator(\n", + " model = models.MontbrioPazoRoxin(variables_of_interest=['r']),\n", + " connectivity = conn,\n", + " integrator = integrators.HeunStochastic(\n", + " dt = 0.01, \n", + " noise = noise.Additive(nsig=np.r_[0.01,0.02])\n", + " ),\n", + " initial_conditions=np.zeros( (1,2,N,1) ),\n", + " stimulus=stimulus,\n", + " monitors=[monitors.BoldBalloonWindkessel(period=10., V0=0.02, tau_s=0.65, alpha=0.32, include_svars=True)],\n", + "\n", + ").configure()\n" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "id": "6ebcb73b-74b6-47bc-9a05-0d8421784dd5", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "b_t, b_d = replay_monitor(raw_d, sim)" + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "id": "42aefae9-b422-468a-85b2-061d794dc980", + "metadata": { + "tags": [] + }, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgQAAAI/CAYAAAAIv8H4AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAB9AklEQVR4nO3dd5xcdb3/8ddnyvae3SSbuukhhSRkSei9oyKoNBELxQI/21Uv3qteO2JDFERREVEBu0ZAeockZAOk996zvZcp398fM7s7GzaV3Z2y7+fjsY85c86ZM5/z3Zlz3nOqOecQERGRwc0T7wJEREQk/hQIRERERIFAREREFAhEREQEBQIRERFBgUBEREQAX7wLiKfi4mJXVlYW7zJEREQGxNKlS6uccyW9DRvUgaCsrIyKiop4lyEiIjIgzGzbwYZpl4GIiIgoEIgkmtv/s4ay2x6LdxkiMsgoEIgkmF++uLmr+7pfL1Y4EJEBoUAgkqBe31LDKxur4l2GiAwSCgQiCerKXy6MdwkiMoikTCAws9Fm9ryZrTazVWb2mXjXJCIikixS6bTDIPBfzrk3zCwXWGpmTzvnVse7MBERkUSXMlsInHN7nHNvRLsbgTXAyPhWJSIikhxSJhDEMrMyYA6wOM6liIiIJIWUCwRmlgP8Dfisc66hl+E3m1mFmVVUVlYOfIEiIiIJKKUCgZn5iYSBPzrn/t7bOM65+5xz5c658pKSXi/nLCIiMuikTCAwMwN+A6xxzv043vWIiIgkk5QJBMCpwIeAc8zsrejfJfEuSkREJBmkzGmHzrlXAIt3HSIiIskolbYQiIiIyDFSIBBJAs3twXiXICIpToFAJAlsq26JdwkikuIUCESSwHW/0TW2RKR/KRCIJIGa5o54lyAiKU6BQCRJhMIu3iWISApTIBBJEs+v3R/vEkQkhSkQiCSJkNMWAhHpPwoEIklCV90Skf6kQCCSJB5dvifeJYhIClMgEEkSC5btjncJIpLCFAhEREREgUAkmVRsrYl3CSKSohQIRJLITQ9WxLsEEUlRCgQiIiKiQCCSTGpbAvEuQURSlAKBiIiIKBCIiIiIAoGIiIigQCAiIiIoEIiIiAgJGAgsYnS86xARERlMEi4QOOcc8Hi86xARERlMEi4QRL1hZice7YvM7H4z229mK/ujKBERkVSVqIFgPrDQzDaZ2XIzW2Fmy4/gdQ8AF/VvaSIiIqnHF+8CDuLCY3mRc+4lMyvr41pERERSXkIGAufctnjXIJKI8jP98S5BRFJUou4y6DdmdrOZVZhZRWVlZbzLETkql88ZGe8SRCRFDbpA4Jy7zzlX7pwrLykpiXc5Ikdl5a76eJcgIilq0AUCkWSW5tNXVkT6R0otXczsYWAhMMXMdprZDfGuSaQv+b0p9ZUVkQSSkAcVHivn3DXxrkGkPykQiEh/0dJFJIkMyU6LdwkikqIUCESSiM9r8S5BRFKUAoFIgnropvl89rxJ8S5DRAYJBQKRBHXKhGI+e97keJchIoOEAoFIgrvn2hPiXYKIDAIKBCIJ7tLjS+NdgogMAgoEIkmgOCdydoGLcx0ikroUCESSgs4uEJH+pUAgIiIiCgQiIiKiQCCSFD5x5ngAppXmxbkSEUlVCgQiSeCMyZFbdRdm6dLFItI/FAhEkojTeQYi0k8UCOSYhMKOjmA43mUMGjrHQET6mwJBP3lpfSX7G9viXUa/ef8vXmPyV/4T7zKSRjjseOT17SkTonbUtNAeDMW7jAHX2BagLZB88/3ph9/kN69s6ff3qWpq5/bH1xAKO4KhMPsbUncZmIoUCPrJ9fe/zhU/f63reTjsqG5qZ9Xuep5ft5+2QAjn+m7zbzAUPuIF1f6GNu57adM7ev83t9cddhznHH9dupPm9mCP/u3BEM+v3X9Er//9om3UtwSOqrZAKMwL6w4//UO97976wy/I6lsDfO8/awmEDr+Sf2zFHm77+wruenb9Mdd1pMq//TSf+uNSAJrbgzS29Wy/7dUtlN32GHO++RS761qPevptgRCnf/95vvCX5b0OL7vtMX7yzJHPZ31rgDufXk8ofPSfx+U766hv7f3zEQ47dtS0sHBT9UFfHwyFWbmr/ojfb+bXn+KcH77Adx5b3Wu4W7mrnsWb3/5+G/c3HfH3bXt1Cx/6zeIe35udtS00HfA9OhI7alp4fu1+FizbzbceXX3IccOHaf+m9iCVje1ApN3KbnuM2/7W8zNw299W8MuXNvPKxiq+/dga5n33WepaOo667qPVHgzxzOp9/f4+v3hxE+v3NR7RuM45WjtCPLum/+vqK754F5CKOr9YO2tb+Z9/rKCmqYNAKMyzB6wEi3PS8XmMRf9zLnvqW2kLhBmam05LR4jCLD+1LQFKctO7xq9vCZCf5e96Xt3UzhX3vsaMEfnsrGtl2Y46tn7vUnbVtVKck0aa18PG/U28srGKCSU5TB+Rx5CcdD71xzeo2FbL5GG5nDVlKB3BMFVN7Zzyvee4+9o5vOv4ET3qbGgL4DXD6zH+sGgbHzt1XNewD/1mMfd8MHKt/byMSG07alr4z8o9fPfxtQC8tqmKW8+eyLq9jRRkpfHBXy8i7OBvnzyZ9927kHfPGsGdV87Cou/RafnOer76z5W8vL6S+64vxzlHMOzweyM5dktVM8Py0qltCfCXih10BMN87vzJ/ODJddz30mYAstK8/OuWUxmWn8Ed/1nL/1xyHNnpkY/9wk3VTB2ey/l3vsT1J4/l0+dOYvHmaq66bxEAnzl3Ev/vnIlUNrVTnJOOx4w3t9dSXlYEwFf/uZIFy3YzdXgu750zEogsBNbta+Rnz27kzCklXFk+GoCWjsjCvHOB6pzjq/9aycyR+fz331Zw19WzuWz2yKP69Vnb3MGcbz3Ndy+fyfzxRRRmpXHTgxVUNXXw+Iq91LcEKP/O0wRCjq3fuxSAfQ1tPB1dQNW2BPjBk+u486rZvU5/V10rQ7LTyPB7eW1jFZOH51Kck961Yvr3st0cPzKfm84Yz8d/X8H8cUP42GmRz8ZPntnAZ86dhJl1zf+u2la++/ga/u/d0ykrzuaEbz1NTXP3yqJiWw2vbqxmQkk237psBtf+ejE/vnIWJ40fwv/+YwXfuXwmNc0djC7KIj/Tz/6GNt5z96vMGJnHqt0NOAdvfe18CqIHXt7xxFp+Gf0c5Gf6OXtKCU3tQeaPG8JHTi3D7/Xww6fW84sXN/G3T57Cr17azLfeO4OS3HQ6gmF+/sJGPnHmBBZtrsbn8XDqxCEA7K5v41cvb6E0P5OPnlqGmdEWCBEKO971s1cA+MenTuHD97/OC188m4cWb+OHT63nK5ceR1aaj8vnjCQzzcufK3ZQ2djOtNI8RhdlMXZIFgbc8eRaXt5QxWPL9/Ce2SP49mOr+cOi7V3T3VTZzGWzR+D3ethd18qIgkwAlm6r4bm1+/nihVO58XdLWLipmsABu/bKbnuMtd+6iGU76pg8LJc/Lt7GJ86cwOaqZi648yXu+9BcLpg+nC1VzQzJSSMYcnzuT2/x4vrKrmlk+D3cEP0/P7JkB9973/Fdwzo/v/vq23hq1V4Almyt5fxpw3hw4VYeeX0Hj336NFo6QlQ2tnPWD1/gk2dN4HPnTeaWh97guNI8SnLSKCvOZt64ItJ9Xmqjn5HC7DR21LQwqjATM+MvFTtoC4T40MllfPexNfxu4TbuvGoWU4fnUZDlx+sxhuZmRD5bW2sYlpfB6KIsym57DIA3vno+Rdk9D9LdsK+Rr/1rFfd/5EQy07z8661dfOaRt0j3eRhXnM3avY3c9cwG1nzroq7XPL92Px99YAkVXzmP4pzu5fUNv6vguegy/88fP5kRBRks2VrDe2dHlhVmxr6GNu55fiOfPncSxTnpNLQFaGkP8buFW6lp6uBr757WtbwaCNaXv1KTTXl5uauoqOiTadU2d/Dv5bvxeoyv/nMlx/Bjp4d0n4f2Xn6BnDm5hBfXV7L6mxcy7WtPHvV0v3XZdL76r1Vdz+eOLWTpttoe45wyYQivRX9VffHCKfzgyXVHNO2Xv3Q2X/rrchYe8Avp5PFD3tYPYHxJNpsrm7ueF2T5OWncELLTffztjZ09xv3au6bxzegvnAunD+PJVe8sdXuMt/2Pnvn8GZz345cO+po0r4eOUJgPnTSW3y/a1mPYrWdP5O7nN3LZ7BH8663dXf2//77j+dIBv6LmjCng9Ekl/PTZDYessXMlDpEF1fl3vsR5xw3ljvcdz5CcdB5cuJWvxfwvZ40uYNmOul6nNaEkm00xbX2gu66ezQljClm4qZqq5na+/0T3//wPN8znut8sPmStnbbcfgnjvvw4AIVZftoCYVr7YRP7+OJsNlf1Pj+3XzGTkpx0bnzw2L/bs0bls2xnPV6PHXLLxRcumMzNZ0yI++6zP918UleQPRYTh+awcX9TH1bUU+dy650oG5LF1uoWIPLd6FyxF+ekU9XU3utrTptYzMShOTzw2lYAzp5SwvPruus4sayQJVtru/7Psf/vw31nenPFnJH8/c1dvS5fDnTg8u/HV87i839e9rbxFv/PuQzLyziqOg7FzJY658p7HaZA0DeB4KYHK3h6ADZZyeARGwgeXb6bWx96M47ViEi8xC4L3qlDBQIdQ9BHFAakP+Vm+A8/kojIO6BAIJIESmL2TYqI9AcFApEkoAsSiUh/S6lAYGYXmdk6M9toZrfFux6RvjKID/URkQGSMoHAzLzAPcDFwDTgGjObFt+qRPpGht8b7xJEJMWlTCAA5gEbnXObnXMdwCPAZXGuSaRPpPtS6asqIokolZYyI4EdMc93RvuJJD2fV3czEJH+lUqB4IiY2c1mVmFmFZWV7+xCGSIDRbc9FpH+ljIXJjKzk4GvO+cujD7/MoBz7vaDvaYvL0wUCIUJhMKkeT09Lr8brYVoLYTCDp+3O4eFww6Px3DOYWZvu95552tjdb4mVkcwjD/mV2Tn6wKhcNelfmPfI/YRItcm76w7tl7nIBh2pPk8PcbvCIbxeQwzuuaps/YD5zFW5zQOnO/YdvRaZLrBcOT902I2l/c2Dwe2UzAU7nr/2Jq8Hut6TdhFrlZ44Gtjp9v5rzB7+//hYG3Z+X842Pixz2PfOxR2vdZzsPbr1NweJCvNSyDkut63t88ZdH9uwmHXNU+ddTjH2z6HvX2G2oMh0n3eHtOMLdfMaO0IkeH39PgcdQ47sP4DdX4OO8cJhyPnV8S2TW//+wO/E4f6Hh3sc3ywaff2Her8PHUEw13fjUhbRb4rndM98Hvam9jPfNiB12O0B0OkeT29LhNi36ezjYC3fZ9iP+e9tXlze7DHZXFjX9vSESTD5+2qPxzuvGx45H8TDIVxRD63nce3HPi/i51u7DLiwP9FMOzwxVwpEHp+Xw/833b+P0JhR9i56P+v5zTbg2Ey/N63ffYiwyNt1RYIke57++e0U+znIfYqhm2BENnpvq5hB5vv2Ho6p3fg/+XA9+ntf9HXBsWVCs3MB6wHzgV2AUuAa51zqw72mr4MBCIiIonuUIEgZW5u5JwLmtmtwJOAF7j/UGFAREREuqXMFoJjYWaVwLbDjnjkioGqPpxeslN79KT26Ka26Ent0U1t0VNft8dY51xJbwMGdSDoa2ZWcbBNMYOR2qMntUc3tUVPao9uaoueBrI9Bt1ZBiIiIvJ2CgQiIiKiQNDH7ot3AQlG7dGT2qOb2qIntUc3tUVPA9YeOoZAREREtIVAREREFAhEREQEBQIRERFBgUBERERQIBAREREUCERERAQFAhEREUGBQERERFAgEBEREcAX7wLiqbi42JWVlcW7DBERkQGxdOnSqoPd/nhQB4KysjIqKiriXYaIiMiAMLNtBxumXQYiIiIyuLcQiMjBdQTDvLi+kkWbq1m7t4Ha5gA+r1GQlcaUYTnMHl3IGZOLyc3wx7tUEekDCgQi0kNrR4gHXtvKr17eTE1zB+k+D1NL8yjNzyDkHFVN7fxuYTW/enkLaV4PZ00p4SOnlnHy+CGYWbzLF5FjpEAgIl2Wbqvls396kx01rZw9pYTrTy7jtEnF+L099y4GQmHe3F7HEyv38s+3dvHU6n0cV5rHFy+czNlThioYiCQhc87Fu4a4KS8vdzqoUCTi4de385V/rqQ0P4MfvH8WJ08YckSvawuEWPDWbn7+wka2Vrcwf1wRX3v3NKaPyO/nikXkaJnZUudcea/DFAgUCER+/sJGvv/EOs6aUsJPr5lD3jEcFxAIhXn49e3c9cwG6loD3Hj6OD577mQy07z9ULGIHItDBQKdZSAyyP1h0Ta+/8Q63jt7BL+6vvyYwgCA3+vh+pPLePa/zuQDc0fxyxc3c+FPXmLx5uo+rlhE+oMCgcgg9uyafXztXys5d+pQfviBWW87VuBYFGSl8b33Hc/DN52EGVz9q0X84Mm1BELhPqhYRPqLAoHIILWztoXP/ektpo3I42fXzsHXB2Eg1skThvD4p0/nyrmjuef5Tbz/3tfYUtXcp+8hIn1HgUBkEAqEwnz64TcJO7jn2hPISuufE46y033c8f7jufeDJ7C1uoVLf/oyC5bt7pf3EpF3RoFAZBC676XNvLG9ju9eMZOxQ7L7/f0unlnKE589nWmleXz64Tf52r9W0h4M9fv7isiRUyAQGWS2VjVz17MbuHjGcN4za8SAvW9pfiYP33wSN50+jgcXbuPKXyxkZ23LgL2/iByaAoHIIOKc43//uYJ0r4evv2f6gL+/3+vhfy+dxi+um8vmymYu/ekrPL92/4DXISJvp0AgMog8sXIvr26s5ksXTWFYXkbc6rhoxnAe/fRpjCzI5KMPLOEHT64lqLMQROKqXwOBmV1kZuvMbKOZ3dbL8HQz+1N0+GIzK4sZ9uVo/3VmdmFM/61mtsLM3jKzipj+RWb2tJltiD4W9ue8iSSbQCjMHU+sZfKwHK6dPzbe5TB2SDZ//9QpXH1i5CyED/56Mfsb2uJdlsig1W+BwMy8wD3AxcA04Bozm3bAaDcAtc65icCdwB3R104DrgamAxcBP49Or9PZzrnZB1xt6TbgWefcJODZ6HMRiXr49e1srW7hyxcfh9eTGPcayPB7+d77jufHV85i+c56Lvnpy7yyoSreZYkMSv25hWAesNE5t9k51wE8Alx2wDiXAb+Ldv8VONcid0W5DHjEOdfunNsCbIxO71Bip/U74L3vfBZEUkNjW4C7ntnAyeOHcNaUkniX8zZXnDCKBbeeSmFWGh+6fzF3Pr2eUHjwXlZdJB76MxCMBHbEPN8Z7dfrOM65IFAPDDnMax3wlJktNbObY8YZ5pzbE+3eCwzri5kQSQUPvLqV6uYObrt4asLeiXDSsFz+deupXDFnFHc9u4EP/WYx+xu1C0FkoCTjQYWnOedOILIr4hYzO+PAEVzkjk29/rwws5vNrMLMKiorK/u5VJH4a+kIcv+rWzhn6lBmjS6IdzmHlJXm40dXzuL77z+epdtqufSnr/DaJu1CEBkI/RkIdgGjY56PivbrdRwz8wH5QPWhXuuc63zcD/yD7l0J+8ysNDqtUqDXc5mcc/c558qdc+UlJYm36VSkrz38+g5qWwLccvaEeJdyxK4sH82/bj2V3AwfH/z1Yu54Yi0dQZ2FINKf+jMQLAEmmdk4M0sjcpDgggPGWQB8ONr9fuC56K/7BcDV0bMQxgGTgNfNLNvMcgHMLBu4AFjZy7Q+DPyrn+ZLJGl0BMP86qXNzB9XxNyxRfEu56hMHZ7HgltP48q5o7n3hU28955XWb+vMd5liaSsfgsE0WMCbgWeBNYAf3bOrTKzb5rZe6Kj/QYYYmYbgc8TPTPAObcK+DOwGngCuMU5FyJyXMArZrYMeB14zDn3RHRa3wPON7MNwHnR5yKD2j/e3MnehjZuOXtivEs5JjnReyH86vpy9jW08a6fvcKvX95MWAccivQ5i/wgH5zKy8tdRUXF4UcUSUKhsOPcH71AToaPf996WsIeTHikqpraue1vK3hmzT7mjSvi9itmMqEkJ95liSQVM1t6wCn7XZLxoEIROQKPr9jD1uoWbjlrYtKHAYDinHR+df1cvv++41m7p4GLf/Iydz2zQTdJEukjCgQiKcg5x89f2MSEkmwunD483uX0GTPjyhNH88x/ncmFM4Zz5zPrufSnr/D6lpp4lyaS9I44EJhZiZnpsHyRJPDCukrW7Gngk2dNxJMgVyXsS0NzM/jZNXP47UdOpLUjxJW/XMitD72huyeKvAOHDAQW8XUzqwLWAevNrNLMvjYw5YnI0XLOcffzGxlZkMllswfu9sbxcPbUoTz9+TP49LmTeGbNPs750Yt8/4m1NLUH412aSNI53BaCzwGnAic654qcc4XAfOBUM/tcv1cnIkft9S01LN1Wy8fPHI/fm/p7BbPSfHz+/Mk8919ncenMUn7+wibO+sHz/PrlzbR26PgCkSN1uKXFh4BrovcTAMA5txm4Dri+PwsTkWNz9/MbKc5J48ry0YcfOYWMKMjkzqtm889bTmXK8Fy+/dgaTv/+8/zmlS20BRQMRA7ncIHA75x723VDnXOVgL9/ShKRY/XWjjpe3lDFjaePJ8PvPfwLUtDs0QX88caT+PPHT2bysBy+9ehqTrvjOX7yzHqqmtrjXZ5IwjpcIOg4xmEiEgd3P7eR/Ew/1500Nt6lxN28cUU8dNNJPHLzSRw/qoCfPLOBU773HF/66zJW726Id3kiCcd3mOGzzKy3b44BGf1Qj4gco9W7G3hmzT4+d95kctIP99UePE4aP4STxg9h4/4mfvvqFv72xk7+XLGTmSPzubJ8FO+ZNZL8LG3wFNGVCnWlQkkRt/zxDV5cX8mr/32OVnCHUNfSwT/f3MWfKnayZk8DaT4P508bxqUzSzlrSglZaQpTkroOdaXCY/7km9l259yYYy9LRPrK8p11PLZiD7eePVFh4DAKstL4yKnj+Mip41i5q54/V+zgseV7eGz5HjL8Hs6aPJSLZw7njEklFGanxbtckQHzTqJw6l3tRCQJOef47uNrGJKdxsfPHB/vcpLKjJH5zBiZz/+9ezqvb6nhPyv38MTKvTyxai9mcPzIfE6fVMIZk0uYM6ZgUJzGKYPXOwkEg3dfg0gCeXbNfhZtruGbl00nN0NbB46F12OcPGEIJ08YwtffPZ23dtbx0vpKXt5Qxc9f2Mjdz28kw+/h+FEFzB1byNwxhZwwtpAibUGQFHLIYwjM7PMHGwT8r3MuuW6wfgAdQyDJrqk9yAU/fpHsdB+Pf+Z0/YLtB/WtAV7bWMWSrbUs3V7Lql31BKO3Xx5ZkMlxpbkcV5rH1OF5TC3NZWxRFj79HyRBvZNjCHIPMeyuYy9JRPrCdx5bzZ6GNv76iVMUBvpJfqafi2eWcvHMUgBaO0Ks2FXP0m21rN7TwJo9DTy3dj/RjIDPY4wszGRMURZjh2Qxtiib0UWZDM3LYFheBiU56aT59L+SxHPIQOCc+8ZAFSIiR+ePi7fx8Os7+ORZE5g7tjDe5QwamWle5o0rYt647g2kbYEQG/c3sWZPA1uqmtlW08L26haW7dhNQ9vb76swJDuNoXkZlOSmU5DpJ/+Av7xMP3mZPjL9XrLSIo+ZadE/vxdvCt6wSuLvsMcQmNnFwJeBadFeq4A7nHOP92dhIvHknKM9GI7+hWgPxHQHwwRDjlDYEXaRv1DY4RyEwo6QczjnCIXp6obIfmqPGR4j+mh4PDHdMc/TfR5yM/wUZae9bT91OOy4/9UtfOfxNZw5uYQvXDAlHk0kMTL83q4DFA9U19LBztpW9je2sa+hnX0NbexvbGd/9HFbdTP1rQEaWgNdWxkOJ83nISsaDtJ9Hvze6J/PQ5rX8Hm6u/1eDz6vB7/XSPN6enwOzWK7I5896/EZjYzjAJzDRR5wuOgjhKMdkWHd/Z2LDuvs38tru8Z/W//IfDp6TjvSj56v72yUaIfF1P62+fFwwPzZocc/oI08nreP7zGLtKkn0s9rnd2G13PAOId4DdY9D7FtMCQnjeNK847sg/EOHTIQmNlNwMeBLwGdO9vLge+Z2Sjn3H39XJ/IEXPO0dQepK4lEPlr7aC2JUB9a4Dm9iDN7UGaoo/N7aGu7qb2IM0dQVo7ulf4HcFwvGenS3FOGlOH5zGqMJN0n4eFm6tZv6+J844bxt3XztGvxQRXkJVGQVYa8PawECscdjR1BKmPfmYb24K0BUK0dIRoDYRo7QjS2uN5pLsjGCYQChMIuehj5K+1NRTzvHtYMOS6VuThaJCNBNvIY+eKtut5TEjpXG+ZGUZk5UZXPzCsa5yew+xtr+3Zr+drzazne8YMh57T6XyPTrFhJDKPPecn7Dqf95zHcEy/g81/PFw8Yzj3Xjd3QN7rcAcVrgZOc87VHNB/CPCKc+64Q07c7CIixxp4gV875753wPB04EFgLlANXOWc2xod9mXgBiAEfNo596SZjY6OP4zI//0+59xd0fG/DtwEVEYn/z+H24qhgwqTQ1sgxO66VvY1tEd/ZbVFuyO/tmqaO6hr6aCuJdB1sFdvzCAnzUd2uo/sdC856Z3dPnLSfWT4vWT4PaT7Ir+60mO7fR7So7/G0nwe0rweLJrsvZ7Ir4rYXwDemF8IZhb9FRRZyHRuWXAusgUhHLNFoXMh3R4M09geZH9DG+v2NrJuXyO769poD4SYPDyXD84fw+VzRvZYcIr0B+fcoP6cOddLaIp+dyNbBju3FtK1tbAziIRdZIthONzzux/u5TWdOlu6MzAVZvmZOPRQh/MdnXdyUKEdGAYAnHPVh/uAmJkXuAc4H9gJLDGzBc651TGj3QDUOucmmtnVwB3AVWY2DbgamA6MAJ4xs8lAEPgv59wbZpYLLDWzp2Omeadz7oeHmSdJQKGwY2t1M+v3NrK5qpnt1S1srW5me00Lexva3pbSM/1ehudH9sFOHpYT+RWW6acwK438rMhjQZafwiw/eRl+cjIi+2EH84JN5FgM9u9M54rZMwguvXO4QNBgZrOcc8tie5rZLKDxMK+dB2yM3i4ZM3sEuAyIDQSXAV+Pdv8VuNsin77LgEecc+3AFjPbCMxzzi0E9gA45xrNbA0w8oBpSoLrCIZZvaeBN7fXsnp3A2v3NrJhfyNtge7N9MU56ZQNyeLkCUO6jtIenp/B0NwMhuWlk5PuG/QLKhGRvnS4QPBfwAIz+y2wNNqvHPgwcN1hXjsS2BHzfCcw/2DjOOeCZlYPDIn2X3TAa0fGvtDMyoA5wOKY3rea2fVEjnf4L+dc7WFqlAHQ0hFk8eYaFm2u5o3ttSzfWU97dB995/7x6+aPZcrwXKYOz2N8STbZujmPiMiAOtxph6+Y2XzgU8BHor1XAyc55/b2c20HZWY5wN+AzzrnOu/GeC/wLSLHFnwL+BHwsV5eezNwM8CYMboVQ39wzrF6TwMvrq/k5fVVLN1WS0coTJrXw/SReVx30lhOGFPICWMLKM3PjHe5IiLCEZx2GF3xf63zuZkVEzkA8HB2AaNjno+K9uttnJ1m5iNyGG71oV5rZn4iYeCPzrm/x9S5L6bGXwGPHmR+7gPug8hBhUcwH3IEnHMs21nPf1bs4fGVe9hR0wrA1OG5fOTUMk6fVMyJZUVk+L1xrlRERHpzuNMOTwK+B9QQ+dX9e6AY8JjZ9c65Jw7x8iXAJDMbR2RlfjVw7QHjLCCy+2Eh8H7gOeecM7MFwENm9mMiBxVOAl6PHl/wG2CNc+7HB9Ra6pzbE316ObDy0LMufWFzZRN/WbqTBW/tZlddK36vccqEYm45ayLnHDeUobkZ8S5RRESOwOG2ENwN/A+RX+7PARc75xaZ2VTgYeCggSB6TMCtwJNETju83zm3ysy+CVQ45xYQWbn/PnrQYA2R0EB0vD8T2T0RBG5xzoXM7DTgQ8AKM3sr+ladpxd+38xmE9llsJXI9ROkHzS1B3ls+W7+XLGTpdtq8XqMMyeX8PnzJ3PeccN0+10RkSR0uOsQvOWcmx3tXhN73QEze9M5N6f/S+w/ug7B0dle3cIDr23lzxU7aGoPMqEkmyvLR3P5nJEMzdOWABGRRPdOrkMQe7m21gOGaf/7IOCcY9HmGu5/dQvPrNmH14xLjy/l+pMjBwbq1D8RkdRwuEAwy8waiFw8KTPaTfS5fhKmMOccL6yv5KfPbuDN7XUUZadxy1kT+dDJYxmmrQEiIinncKcd6pDwQcY5x7Nr9vPT5zawfGc9Iwsy+fZ7Z/D+uaN0hoCISArT1V+ky+LN1Xz3P2tZtqOO0UWZ3PG+mVw+Z5Tu3S4iMggoEAgb9zfxvf+s5Zk1+xiel8Ed75vJFSeMwu9VEBARGSwUCAaxmuYOfvTUOh5ZsoNMv5cvXjiFj506jsw07RoQERlsFAgGoXDY8aeKHdzxxFqa2oJ8cP4YPn3uJIpz0uNdmoiIxIkCwSCzanc9X/nnSt7cXse8cUV8+70zmDys7+61LSIiyUmBYJBobAvw46fX87vXtlKYlcaPPjCLK04YqesIiIgIoECQ8pxzPLp8D996dDWVTe1cO28MX7pwqi4vLCIiPSgQpLAtVc187V8reXlDFTNG5nHf9eXMHl0Q77JERCQBKRCkoLZAiJ+/sIlfvLCJdJ+Hb7xnOtedNBavR7sHRESkdwoEKeb5tfv5vwWr2F7TwmWzR/C/lx6nWxCLiMhhKRCkiF11rXzz36t4ctU+JpRk89CN8zllYnG8yxIRkSShQJDkOoJhfvPKFn767AYAvnTRFG48bbwuNywiIkdFgSBJOed4bu1+vvv4GjZVNnPBtGF87d3TGFWYFe/SREQkCSkQJKGVu+r5zmNrWLi5mvHF2dz/kXLOmTos3mWJiEgSUyBIImv3NnDP85v497LdFGWn8c3LpnPNvDG6CZGIiLxj/bomMbOLzGydmW00s9t6GZ5uZn+KDl9sZmUxw74c7b/OzC483DTNbFx0Ghuj00zrz3kbKM45lm6r5aYHK7joJy/z3Jp93HL2BF744llcf3KZwoCIiPSJfttCYGZe4B7gfGAnsMTMFjjnVseMdgNQ65ybaGZXA3cAV5nZNOBqYDowAnjGzCZHX3Owad4B3Omce8TMfhGd9r39NX/9rb41wIK3dvHHxdtZu7eR/Ew/nz1vEh85pYyCrJTIOiIikkD6c5fBPGCjc24zgJk9AlwGxAaCy4CvR7v/CtxtkYvrXwY84pxrB7aY2cbo9Ohtmma2BjgHuDY6zu+i002aQBAOOzbsb2LxlmqeXr2PhZuqCYYd00fk8Z3LZ3DZ7JHkpGsPj4iI9I/+XMOMBHbEPN8JzD/YOM65oJnVA0Oi/Rcd8NqR0e7epjkEqHPOBXsZf0Cs39fIhn1NhJ0j7BzOEe2OPoa7u120f3NHkO3VLWytbmbV7gYa2yLljy/O5sbTx3PpzFJmjsofyNkQEZFBatD95DSzm4GbAcaMGdNn0/33st387LmNR/264pw0xg7J5l3Hj6B8bCFzxxYydkiW7kIoIiIDqj8DwS5gdMzzUdF+vY2z08x8QD5QfZjX9ta/GigwM190K0Fv7wWAc+4+4D6A8vJyd/Sz1bvrTy7j3bNG4DEwMzxmeAw8Zlj00WOGx9Pdne7zkK3dACIikgD6c220BJhkZuOIrJyvpnsff6cFwIeBhcD7geecc87MFgAPmdmPiRxUOAl4HbDephl9zfPRaTwSnea/+nHe3qYkN52S3PSBfEsREZE+02+BIHpMwK3Ak4AXuN85t8rMvglUOOcWAL8Bfh89aLCGyAqe6Hh/JnIAYhC4xTkXAuhtmtG3/G/gETP7NvBmdNoiIiJyBMy5PttqnnTMrBLY1oeTLAaq+nB6yU7t0ZPao5vaoie1Rze1RU993R5jnXMlvQ0Y1IGgr5lZhXOuPN51JAq1R09qj25qi57UHt3UFj0NZHvoMnciIiKiQCAiIiIKBH3tvngXkGDUHj2pPbqpLXpSe3RTW/Q0YO2hYwhEREREWwhEREREgUBERERQIBAREREUCERERAQFAhEREUGBQERERFAgEBERERQIREREBAUCERERQYFAREREAF+8C4in4uJiV1ZWFu8yREREBsTSpUurnHMlvQ0b1IGgrKyMioqKeJchIiIyIMxs28GGaZeBSBJoC4T4/hNrqW8JxLsUEUlRCgQiSeDvb+zi5y9s4t4XN8W7FBFJUQoEIklgV10LAO3BUJwrEZFUpUAgkgQ27m8CYHdda5wrEZFUpUAgkgTW74sEgs2VzXGuRERSlQKBSIJrC4TYWh0JAtuqWwiFXZwrEpFUpEAgkuA27m/COThzcgkdoTC7arXbQET6ngKBSIJbt7cRgItmDAdgc1VTPMsRkRSlQCCS4LbVtOAxOGtK5OJiOo5ARPqDAoFIgttZ28LwvAyG52WQm+FjS5UCgYj0PQUCkQS3s7aVUYVZmBnjS3K0y0BE+oUCgUiC21XbyqjCTAAmFGdrl4GI9AsFApEEFgiF2VPfHQjGFWezp76Nlo5gnCsTkVSjQCCSwPY3thN2UFoQCQTjS3IAdByBiPQ5BQKRBLa/oQ2AobnpQGQLASgQiEjfUyAQSWCVje0ADM3NACKBwKz73gYiIn1FgUAkgVU2RQJBSXQLQWaalzFFWazf1xjPskQkBSkQiCSw/Q3tmMGQnLSuflOG5XZdvVBEpK8oEIgksMqmdoqy0vB7u7+qU4fnsrW6hbZAKI6ViUiqUSAQSWD7G9q7dhd0mjw8l1DY6TgCEelTCgQiCayy6e2BYOrwXADtNhCRPpU0gcDMLjKzdWa20cxu62X4581stZktN7NnzWxsPOoU6UtVje2U5PQMBGOHZJPm9ejAQhHpU0kRCMzMC9wDXAxMA64xs2kHjPYmUO6cOx74K/D9ga1SpO/VNHdQlJ3Wo5/f62HC0BzWaguBiPShpAgEwDxgo3Nus3OuA3gEuCx2BOfc8865lujTRcCoAa5RpE+1BUK0BkIUHhAIILLbQLsMRKQvJUsgGAnsiHm+M9rvYG4A/tOvFYn0s/rWAAD5mf63DTuuNJe9DW1URa9TICLyTiVLIDhiZnYdUA784CDDbzazCjOrqKysHNjiRI5CXUskEBRkvT0QzBxZAMCKXfUDWZKIpLBkCQS7gNExz0dF+/VgZucB/wu8xznX608n59x9zrly51x5SUlJvxQr0hdqWzoAKMx6+y6DGSPzAFi5U4FARPpGsgSCJcAkMxtnZmnA1cCC2BHMbA7wSyJhYH8cahTpU51bCHrbZZCb4Wd8STbLtYVARPpIUgQC51wQuBV4ElgD/Nk5t8rMvmlm74mO9gMgB/iLmb1lZgsOMjmRpFDfGt1C0MtBhQDHj8xnhbYQiEgf8cW7gCPlnHscePyAfl+L6T5vwIsS6Ue1nccQ9LKFAGDGyHz++dZu9je0MTQvYyBLE5EUlBRbCEQGo7qWAGleD1lp3l6HHz+qANCBhSLSNxQIRBJUXUsH+Vl+zKzX4dNH5GEGy7XbQET6gAKBSIKqawkcdHcBQHa6j4klOSzfWTdwRYlIylIgEElQda0dvZ5yGOuEMYUs3VZLOOwGqCoRSVUKBCIJqq4lQH4vFyWKdeK4IhragmzQrZBF5B0a8EBgZh8ws9xo91fM7O9mdsJA1yGS6OpaAhQeLhCUFQKwZGvNQJQkIiksHlsIvuqcazSz04DzgN8A98ahDpGEVtvSQcFhdhmMKcqiJDedCgUCEXmH4hEIQtHHS4H7nHOPAYde6okMMm2BEO3BcK9XKYxlZswrK2LJ1toBqkxEUlU8AsEuM/slcBXwuJmlx6kOkYR1qPsYHKi8rJBdda3srmvt77JEJIXFY0V8JZFLEF/onKsDioAvxqEOkYR1qDsdHujEsiJAxxGIyDsz4IHAOdfinPu7c25D9Pke59xTA12HSCI7mkAwdXguuRk+Fm2u7u+yRCSFaVO9SAKqi+4yKMg8/C4Dn9fDKROG8NL6KpzT9QhE5NgoEIgkoLrWyBaCwuzDbyEAOGNyCbvqWtlc1dyfZYlIClMgEElAtUexhQDgjEklALy8vrLfahKR1KZAIJKA6lsCpPk8ZPiP7Cs6uiiLsiFZvLyhqp8rE5FUpUAgkoA6r1J4sDsd9ub0SSUs3FxNezB0+JFFRA6gQCCSgGpbOo54d0Gns6eW0NIR4rVNOttARI6eAoFIAqprDRzRKYexTp1YTG66jydW7O2nqkQklSkQiCSgupaOow4E6T4v5x43lCdX7yUQCvdTZSKSqhQIRBJQ5BiCo7/Fx7tnjaCuJcAL63S2gYgcHQUCkQTjnKOuNUD+UW4hgMj1CIZkp/G3pTv7oTIRSWW+eBcwWDW2BdhW3UJlUzvN7UGCIUcgFCYYdgRDYXIyfJTmZzIiP5Nh+emk+7zxLlkGSGsgREcwfNQHFQL4vR4umz2S3y/aSm1zB4XZupGoiBwZBYIBtHJXPY+v2MNLGypZvbuB8FFcZbY4Jy0SEAoymFCSw5ThuUwamsuEodkKCymm8z4GhcewhQDgfXNHcv+rW/j38t1cf3JZH1YmIqlMgWAALN9Zx4+eWs+L6yvxeoy5Ywu59ZxJTCvNoyQ3jdwMPz6P4fd68HkNn8dDQ1uAPXVt7K5vZU9dG3vqW9ld38bG/U08u2Y/wWia8HqMMUVZlOZnMKowkxPLipg1uoAJJTl4PUd+Drskjq6rFB5jIJg+Ip/jSvP429KdCgQicsQUCPrRmj0N3Pn0ep5avY+CLD//fdFUrp035oj2DZfkpjOhJKfXYR3BMFurm1m3t5H1+xrZVNnE3vo2nlq9jz9XRPYdZ/g9TBmWy5gh2QzJTmNUYSYzR+YzbUQeuRnHtqKRgVHfdafDY9/c//65o/jWo6tZtbue6SPy+6o0EUlhSRMIzOwi4C7AC/zaOfe9A4anAw8Cc4Fq4Crn3NaBrhMiK+wfPrWOX728mZw0H587bzIfO62sz1bEaT4Pk4flMnlYbo/+obBjS1UTy3bUs2p3A2v3NrB8Zx01TR00tgcBMIMJJTnMGlXAkJw0JpRkc1xpHiMLMsnL9OP36jjTeKs9ilsfH8z7TxjFD59cxwOvbuUHH5jVV6WJSApLikBgZl7gHuB8YCewxMwWOOdWx4x2A1DrnJtoZlcDdwBXDXStlY3t3Pi7JSzbWc+188fwpQunvKNfekfD6zEmDs1l4tBc3je357Ddda2s29fIip31vLm9llc2VlLbEqAj2PN89aw0L3kZfvIz/eRl+ijISmN8cTaThuUyoSSb/Ew/WWk+irLTSPMpPPSHmuZ2AIrewecmP8vPFSeM5C9Ld3LbxVMZkpPeV+WJSIpKikAAzAM2Ouc2A5jZI8BlQGwguAz4erT7r8DdZmZuAG8Qv7e+jWt/vYg9dW384roTuGhG6UC99WGNKMhkREEmZ08Z2tUvHHZsr2lh7d4G9jW0U98aoKE1EHlsizxuq27mxXWVdBxwoRszKMxKY2huOiMKMhldmMnw/EyG5KSR4feS7vNQnJNGcU46WWk+cjN8+DyG12NHdX3+wai6OXIMwTs9Q+Cjp47jode3c/+rW/jihVP7ojQRSWHJEghGAjtinu8E5h9sHOdc0MzqgSHAgNz+bX9DG9f+ahH7G9t58IZ5nFhWNBBv+454PEZZcTZlxdmHHC8YCrOtpoUtlc00dwRpag9S2djO/sZ29je0sbO2lYqtNTS0BQ//nha5ol6G39P7Y5qXDJ+H9mCYjmAYr8cIhMLkZvjxeiK7YzoPqAyFHV6Pken30hEK09IRIt3nId3nIS/DT3o0mGT4vfi9kRDSGQ8dnQdlRu4o6DUjO91HfnS3iXOOzDQv2ek+mtqDhMOOsIOwc13TzPB7CYXDOBfZvO/zeMhK8+KL2e3i9xrBkKMtGKItEKYtEKItEKI9GCYcdnSEwrQFwrQHQ4TCjt11rfylYmdXHe/ExKE5XDqzlAde3cqNp43XKYj9LBgKdwXecNhR1dxOW0eYpvYgDW0BmtuDpPu8BMJhvGY0twcxi5wqmun3kpXuw4B0v4eOYJjm9hDZ6V7SfB7aAmEa2wIEQmE6gp2nKIdJ83rxGHT+6vF5DL/Pg9/joak90HXGSthBXWsHztH1WfZYZBnQ2BYk0x/5Dmb4I+8Xdo5gyBEKO3IyfNF6gngsEup93uijx0j3eekIhmloC3TV4Ih818LO0RYI0dwR+XyHwg6fx8hK80W/8x4y07xkpvnITvOSmRb5znb+aHDOEQpHvuvOOULORbsj/cIu8hcKE9Mdmb8Mv6erBuciWz8j0/fS1B4gzeslL9NHTrqPUHSZ0v1bxbqeB0JhwuHIY0coTDAUGbe2pYPalg7aAmE8FmnjUDjcVUsw5HA4/F4PHjPMIv29Zl3z2NQepLqpA6+n53LRE23bQCjMmKIsLpg+fAA+wckTCPqMmd0M3AwwZsyYPpvuvS9uYm9DGw9+bB7lSRAGjobP62FCSc5BD3LsVN8S2bIQWemFqWpqp7q5g5aOII1tQULhyIKsc3h78O2P9S0d7A+G8XmNDJ+X9qDD5/Gws7YlujL24vVYZGFmRsg59je0kxZdsNS3RnaDNHc0dq2A2wNhAuEwnd/1zoWNAaHowiKRFOekc/3JY/tkWp8+dxKPrdjDr1/ZrK0EhxAKO3bWtlDbEqC1I0RTe5Cm9gBN7SEa2wK0tIeobelgR20r1U3ttAfDVDa2Ewh1f66aO0J4PYbfawSiK9Nk0Flvf7+H12N4zQiGHe1BXVq7U3aaFwe0BUK9nop+4fRhCgQH2AWMjnk+Ktqvt3F2mpkPyCdycGEPzrn7gPsAysvL++xb8OWLj+MDc0czbUReX00y6eRn+Y/p6nrxFAyFI7/WnaO1I0Rda6DrPgCtHZFfNrkZPrzWubuDSHgJhGjpCOGPHkdR3xogGArT3BEiHP1WO+cIhFwk3HT+AvN5SY/+CvB5DJ/X0/2rwKA0P5PMtL67rsTkYblcMrOU3766lQ+fUsbQ3Iw+m3YiC4cdzR2RX1/7GtoiW7OiW7Qa2gLUNHdQ2dhOVVMHbYEQdS2Bt+0Wi+UxyM/0MzQ3g5GFmaR5PZw8fgjpPk/Xr+GcdC+h6P/c7zWG5WWQ6feSm+EjL8NPVrqva6uXc46sNB+OyPitHSFaA0FC0V+iaV4P6X4P7YHIr1K/1xPdamSkRbeCeT2RLQmO7sVY5wXOAiFHpt/LkJzIViFH5LoWPo8Hi4bpzl/TkS1dkV/yrdGtV77oL1RP9Fdsms9DTroPB4RCjmA4TCi6has9GI5slcuMfPdDoe5f215PZ709P9Od79cWCNEWDNPaEaS5PURzR5BAKLI1wAFei9Tg8US7ozV5otP2RId7PYbXEwn7na9pDYS6toJA5Pvc0hGiPRgiJ91HIORoaA3Q1B7scYp2Z2t27m32eTx4PZGDuf1eT2RcB0U5aRRmpZHh8xJ2LlKH17qWFV6PYUBHKLIV0REZJxQNRJ2fgc7vu3OOYLRdwmEiW4B6abv+lCyBYAkwyczGEVnxXw1ce8A4C4APAwuB9wPPDeTxA2k+z6AOA8nK5/V0beLPzfAzNC/1Vpj/df5knly5l588s4HvXj4z3uX0qUAozNaqZtbva6JiWw0rd9Wzs7aVfQ1tvf7aSvN6yM/yU5gVWbmfMCaLDL+X/Cw/E0pyKI4eA5Ob7icnI7I5OSc9snk71Y598WL4o+saryeyyyw7/e2rhP7Y1XSo90s1GZ63r9CzeznG1yyydSmeZ3olxX8jekzArcCTRE47vN85t8rMvglUOOcWAL8Bfm9mG4EaIqFBZNAbX5LDdSeN5cGFW7l8zsikOL6lNzXNHazZ0xD9a2Tt3gY27Gvq+mWf7vMwa1QBp04spjQ/g7wMPwVZfobnZzA0N4NheenkZ/pTbsUu0ldsAH9EJ5zy8nJXUVER7zJE+l1Te5BL7nqZsHM8/pnTyUvwi1O1dASp2FrLa5uqWbajjs1VTexraO8aXpyTznGluRxXmsfU4ZFrckwcmkOGX5fxFjkUM1vqnCvvbVhSbCEQkXcmJ93HnVfN5gO/eI2vL1jFj6+cHe+SemgLhFi6rZZXN1bx2qZqVuyqJxSO7IufNiKfUycUR1b+pblMHZ5HSa6uqyDS1xQIRAaJznto/PTZDUwamssnz5oQt1qCoTDLdtazcFMVr26sZun2WjqiB7PNHl3AJ84cz/xxQygvKyQrTYspkYGgb5rIIPKZcyexrbqZO55Yi99r3Hj6+AF770AozOtbavjbGzt5evU+GqPXrZhWmsf1J43l1InFnDiuiJxBcKCZSCLSN09kEPF6jB99YBbBkOPbj62hPRjmU2dN6LcD7Vo6gizcVM2CZbt5bs1+GtuD5Kb7uHDGcM6eMpSTJwyhSBdMEkkICgQig4zP6+EnV8/G4zF+8OQ63tpRx7cum8Hw/Hd+ymV9S4AlW2t4fWsNi7dETgMMhR0FWX4unjmcc6YO5awpQ3Xwn0gC0lkGOstABinnHPe/upXv/WcNHjM+fEoZnzxzwlGdd17V1M6SLZGV/+ItNazd24BzkfP9Z48u4MRxhZxYVsQpE4p1MyyRBHCoswwUCBQIZJDbUdPCnc+s5x9v7iInzcclM0uZPaaAccXZFGT5KchMoyMYprE9QG1zgN31rby1o47Xt9SwcX8TAJl+L3PHFjJvXBHzxxUxa3SBtgKIJCAFgoNQIBDptn5fIz97biMvb6jsuinOweSk+ygvK2T+uCHMH1/EjBH52gIgkgR0HQIROazJw3L52TVzcM6xo6aVrdXNNLUHqW3pIN3XfU3+4fkZjCnK6nH9dxFJfgoEItKDmTFmSBZjhmTFuxQRGUDaxiciIiKD+xgCM6sEtvXhJIuBqj6cXrJTe/Sk9uimtuhJ7dFNbdFTX7fHWOdcSW8DBnUg6GtmVnGwgzUGI7VHT2qPbmqLntQe3dQWPQ1ke2iXgYiIiCgQiIiIiAJBX7sv3gUkGLVHT2qPbmqLntQe3dQWPQ1Ye+gYAhEREdEWAhEREVEgEBERERQIREREBAUCERERQYFAREREUCAQERERFAhEREQEBQIRERFBgUBERERQIBARERHAF+8C4qm4uNiVlZXFuwwREZEBsXTp0irnXElvwwZ1ICgrK6OioiLeZYiIiAwIM9t2sGGDOhCIpLLKxnYWba6msrGd/Ew/x4/KZ+LQHMws3qWJSAJSIBBJMVVN7dz++Fr+8eZOwgfczHTS0BxuPmM8V5wwCq9HwUBEuikQiKSQVbvrueGBCmpaOvjoqeN47+yRjCrMpKalg4Wbqvnj4u188a/L+cOibXz3iplMH5Ef75JFJEGYc+7wY6Wo8vJyp2MIJFVsqmziA79YSLrPw68/XN7ryt45x4Jlu/n2Y2uobwnw1Xcdx3UnjdVuBJFBwsyWOufKexum0w5FUkBbIMQn/7AUAx666aSD/vI3My6bPZKnP3cGp04cwlf/tYr//ttyAqHwwBYsIglHgUAkBfzgyXWs39fEnVfNZlxx9mHHL8hK4zcfPpFPnzORP1fs5KYHK2jpCA5ApSKSqBQIRJLcpsomfvfaVq6ZN4YzJvd6enGvPB7j8xdM4fYrZvLS+kqu/dViGtoC/VipiCQyBQKRJHfHf9aS4ffyXxdMPqbXXzNvDPdeN5dVu+u5/jev06hQIDIoKRCIJLGN+xt5avU+bjx9HMU56cc8nQunD+fua09g5a56Pny/QoHIYKRAIJLEfv3yFtJ9Hj500th3PK3OULB8Zz03/K6CtkCoDyoUkWShQCCSpOpaOvj7m7t439xRDHkHWwdiXTRjOHdeNZslW2v41B/f0NkHIoOIAoFIkvr38j10BMN8cP6YPp3uu2eN4Dvvnclza/fzhb8sI3zg5Q5FJCXpSoUiServb+xk6vBcppXm9fm0r50/hvrWAHc8sZa8DD/fvGy6Ll4kkuIUCESS0JaqZt7cXsf/XDK131bUnzxrAnWtHfzyxc3kZ/r5woVT+uV9RCQxKBCIJKEnV+0F4F3Hj+jX97ntoqk0tAa4+/mN5Gf6uemM8f36fiISPwoEIkno2TX7mD4ijxEFmf36PmbGt987k4bWIN95fA15mT6uOrFvj1kQkcSggwpFkkxtcwdLt9Vy7tShA/J+Xo9x51WzOWNyCV/++woeXb57QN5XRAaWAoFIknl+3X7CDs49btiAvWeaz8MvrjuBuWML+cwjb3XtshCR1KFAIJJkXlxfSXFOOjNH9n5Hw/6Slebj/o+cyMyR+dz60Bs8v3b/gL6/iPQvBQKRJPPG9lpOLCvE4xn40wBzM/z87mPzmDI8l4//YSmvbqwa8BpEpH8oEIgkkcrGdnbUtHLCmMK41ZCf6ef3H5vPuCHZ3PC7JSzeXB23WkSk7ygQiCSRN7fXAjBnTEFc6yjMTuMPN85nREEmH31gibYUiKQABQKRJPLmjjp8HmPGAB8/0JuS3HQeuekkRhdm8dEHlvDM6n3xLklE3gEFApEk8ub2WqaPyCPD7413KQAMzcvgkZtP4rjoMQX/emtXvEsSkWOkQCCSJIKhMMt21DMnjscP9KZz98HcsYV89k9v8ftF2+JdkogcAwUCkSSxbl8jrYFQ3I8f6E1uhp8HPzaPs6cM5av/XMntj6/RXRJFkowCgUiSeGN7HUBczzA4lAy/l/s+NJfrThrDL1/azC0PvUFbIBTvskTkCCVFIDCz+81sv5mtPMjwD5rZcjNbYWavmdmsga5RpL+9ub2W4px0RhX27/0L3gmf18O3LpvBVy49jidW7eWq+xZR2dge77JE5AgkRSAAHgAuOsTwLcCZzrmZwLeA+waiKJGB9Nb2OuaMKei32x33FTPjxtPH84vr5rJubwPv/tkrLN1WG++yROQwkiIQOOdeAmoOMfw151znEmcRMGpAChMZILXNHWyuak7I4wcO5sLpw/nrJ04hzefhql8u5LevbsE5HVcgkqiSIhAcpRuA/8S7CJG+9NaOOiBxjx84mBkj8/n3radx1pQSvvHv1dz60JvUtXTEuywR6UVKBQIzO5tIIPjvQ4xzs5lVmFlFZWXlwBUn8g68sb0Wr8c4flT8L0h0tPKz/Nz3oXK+dNEUnly1lwvufInn1+nGSCKJJmUCgZkdD/wauMw5d9CLqzvn7nPOlTvnyktKSgauQJF34M3tdUwdnktWmi/epRwTj8f41FkT+ectp1KQ5eejv13CbX9brq0FIgkkJQKBmY0B/g58yDm3Pt71iPSlUNjx1o66pDp+4GBmjMzn3//vND5+5nj+snQnZ//wBf60ZLuuWSCSAJIiEJjZw8BCYIqZ7TSzG8zsE2b2iegoXwOGAD83s7fMrCJuxYr0sY37m2hqDybd8QMHk+7z8uWLj+PR/3caE4fm8N9/W8Hl977Gwk26a6JIPA3I9kcz+71z7kNm9hnn3F1H+3rn3DWHGX4jcOMxFyiSwDpP2Uu0Sxa/U8eV5vHnj5/MP97cxfefWMc1v1rE6ZOK+eKFUzh+VEG8yxMZdAZqC8FcMxsBfMzMCs2sKPZvgGoQSUoVW2sozkmnbEhWvEvpc2bGFSeM4oUvnsVXLj2Olbvqec/dr3LDA0tYuKlapymKDKCBOkLpF8CzwHhgKRB7ZRUX7S8ivViyrYYTywoT/oJE70SG38uNp4/nqhNH89tXt/K717Zyza8WMWNkHjecNo6LZ5QmzB0eRVLVgGwhcM791Dl3HHC/c268c25czJ/CgMhB7K1vY0dNK+Vlg2NDWm6Gn0+fO4lXbzuH26+YSWtHiM/9aRknfucZ/ucfK3hze622Goj0kwE9h8k598mBfD+RZLd4S+RAu/KxqXX8wOFk+L1cM28MV5WPZtGWav5asZO/v7GThxZvZ1xxNhdMH8ZF04cza1QBHk/qbjkRGUjJeVKzyCDx4vpKCrL8zBiZfBck6gsej3HKhGJOmVDMNy6bzuMr9vDo8j385uUt/PLFzQzLS+ecqcM4deIQTh4/hCE56fEuWSRpKRCIJKhw2PHiukrOnFyCV7+Cyc3wc9WJY7jqxDHUtwR4bt0+nly5j0eX7ebh17cDMHV4LvPGFTFjZD7Hj8pnYkkOPm9SnF0tEncKBCIJaun2WqqbOzhn6tB4l5Jw8rP8XD5nFJfPGUUwFGbFrnpe21TNwk3V/G3pTh5cuA2ADL+HaaV5TB+Rz8ShOUwoyWHi0ByG5aWn9EGaIsdCgUAkQf2lYgfZaV7OnzYs3qUkNJ/Xw5wxhcwZU8gtZ08kHHZsqW5mxc56VuyqZ8XOev755i4a24Ndr8lJ9zGhJJtRRVmMKshkREEmIwsyGVmYyfC8DPIz/To2QQYdBQKRBLS7rpV/vrWb950wMmnvXxAvHo8xoSSyNeC9c0YC4JyjsrGdjfub2FTZFH1sZvXuBp5evY+OYLjHNLweozDLT1F2GkXZaQzJSacoK43sdB/ZaV6yYh4z/V78XiPN68Hn9eD3Gn6vhzSfB5/ngG6fJzKex/B6LOW3UjjnCLvux3D0DJHwAf1xnf0cjsh56R6zyJ+nu9uss5uu56nehgNJSxqJK+ccgZCjLRiirSNEWyBMayBEWyDU9dgWCNMRCtMeCNEeDNMeDNMRDNMejD4PhOkIhWgP9BzWEQoTCjvC4cjCJuQc4XBkARQKRxY+nY+x/cLhyLhGZKHd+ecx8Hk8eDyG1wNeix1m+LyRR6/Huhb4kb/uFUDPRw8+b3TFQOT9g+FITa9urMJj8KmzJsb7X5QSzIyheRkMzcvglInFPYaFw46q5nZ21bayq66V/Q3t1DR3UN3cQXVTpHvN7gZqWjpoaY98rvq2tu6VnJnFrAwjjxxkJeix7hUl0HPF2WO6kdcAuOiK1xF9PGBl3fm8c3hkpd05TrR/ODoeMSvxzvHoOa2B4LFDB4bO76dZ9/fWE/O99cSM0/nd9HQ+xnyve36nu/t5PIbXLCbg9GyD2GVMb8M727BzOXVgWDplwhC+fMlxA9KWCgTSJ9oCIaqa2qlq6qCqsZ2qpnaqmzuoa+mgsS1IY1uQhrYADW1BGtsCNLYFaWkP0haMrLSPVVr011d69C/S7SXd78Hv9US//OD3eLoWGp0r967uzoWBEdNtOByh6Jc0GI4GhehKOzZMBEORANEeCHcNC4Zc1+tCMX/BcLh7xR/qHh52Dp/X8Hk8eD1GXqaPe649gdFFqXd1wkTj8RhDczMYmptxRJeH7giGae0I0dwRpKUjSGtHJLAGQ2ECIUcgFI7+Hbo7dsHfvWKIWSGHe/tV/faVdNc4MdPoXHm7mBWNI7IytIMEi0iI6Bz+9nGMmNByYIiJGX6w13Xugul+Hd1h5oBpEZ3/0MFWop3B3rlDDu9ss1Dn83D3D4OQi+l3QP/u7+vBv9ed4wZDkRrsMFszYrd29Ax0dL22czkVOzwv09/3H/qDUCCQI1bfEmBjZRPbqpvZWt3C9ujjtupmalsCvb4mw+8hL8NPboaP3Aw/eRk+RhVkkpvhIyvNR2aah0y/l4zoX6bfS2aalwy/p6tfRnQFn+b1kO6PrvCjm161n1cGWlo0eOZnDdyCWmQgKBBIryob21m6rZaVu+pZs6eBNXsa2F3f1jXcDEbkZ1JWnMVFM0oZWZBBcU565C83neKcNIpz0nW5WRGRJKFAIITCjvX7Glm6rZY3ttVSsa2W7TUtQGTf2oSSbE4cV8RxpXlMGppDWXE2owozSfdpZS8ikioUCAapnbUtvLS+ipfWV/Lqpioa2yKnZBXnpFM+tpAPnTSWuWWFTCvN0698EZFBQIFgkHDOsXpPA48t38MTq/ayubIZgBH5GVw6s5T544uYO6aI0UWZOo1HRGQQUiBIcTtqWvjr0p38e9luNlc14/UYJ48fwgfnj+XMycVMKMlRABAREQWCVNQeDPH06n38ackOXtlYBcDJ44dw4+njuXD6MN0ARkRE3kaBIIXUNHfw4MKt/H7hNqqbOxhZkMlnz53M+8tHMbIgM97liYhIAlMgSAHbqpv59ctb+MvSHbQFwpx33FCuP7mM0yYW6zx9ERE5IgoESWxLVTN3PbOeBct24/N4eO+cEdx8xngmDs2Nd2kiIpJkFAiS0K66Vn76zAb++sZO/F7jptPHc8Np4xialxHv0kREJEkpECSRhrYAP3l6A39YFLnX+4dOGsunzp7A0FwFAREReWcUCJKAc45/vrWL7zy2lurmdq6cO5pPnzdJBwqKiEifUSBIcOv3NfLVf65k8ZYaZo3K5/6PlHP8qIJ4lyUiIilGgSBBNbUHueuZ9fz21a3kZPj47uUzufrE0TprQERE+oUCQYJxzvHo8j18+7HV7Gto5+oTR/Oli6ZSlJ0W79JERCSFKRAkkI37m/j6glW8srGKGSPzuPe6uZwwpjDeZYmIyCCgQJAAWjqC/Oy5jfz65c1k+r1867LpXDt/LF7tHhARkQGiQBBHzjmeXLWXbz26hl11rbx/7ihuu3gqxbrXgIiIDDAFgjjZVBnZPfDyhiqmDs/lL584mRPLiuJdloiIDFIKBAOsqT3I3c9t5DevbCbD5+Xr757GdSeNxef1xLs0EREZxJIiEJjZ/cC7gP3OuRm9DDfgLuASoAX4iHPujYGt8tDaAiH+sGgbP39hEzXNHbzvhMjugZJc7R4QEZH4S4pAADwA3A08eJDhFwOTon/zgXujj3FX09zBw69v58GFW9nX0M7pk4r5wgVTmDW6IN6liYiIdEmKQOCce8nMyg4xymXAg845BywyswIzK3XO7RmYCrs559hR08qSrTU8sWovL66vpCMY5rSJxdx51WxOmVA80CWJiIgcVlIEgiMwEtgR83xntN+ABYLfvLKF+1/ZQk1zB62BEACl+Rl8cP4Yrpk3hsnDdEtiERFJXKkSCI6Ymd0M3AwwZsyYPpvuiPwMTho/hMIsP+NKspk1qoBppXm61LCIiCSFVAkEu4DRMc9HRfu9jXPuPuA+gPLyctdXBVw8s5SLZ5b21eREREQGVKqc67YAuN4iTgLq43H8gIiISLKyyHF4ic3MHgbOAoqBfcD/AX4A59wvoqcd3g1cROS0w4865yqOYLqVwLY+LLUYqOrD6SU7tUdPao9uaoue1B7d1BY99XV7jHXOlfQ2ICkCQbIwswrnXHm860gUao+e1B7d1BY9qT26qS16Gsj2SJVdBiIiIvIOKBCIiIiIAkEfuy/eBSQYtUdPao9uaoue1B7d1BY9DVh76BgCERER0RYCERERUSAQERERFAhEREQEBQIRERFBgUBERERQIBAREREUCERERAQFAhERESHBAoGZ3W9m+81s5UGGTzWzhWbWbmZfOGDYRWa2zsw2mtltA1OxiIhIakioQAA8QOQWxgdTA3wa+GFsTzPzAvcAFwPTgGvMbFo/1SgiIpJyfPEuIJZz7iUzKzvE8P3AfjO79IBB84CNzrnNAGb2CHAZsPpQ71dcXOzKyg76diIiIill6dKlVc65kt6GJVQgeAdGAjtinu8E5h/uRWVlZVRUVPRbUSIiIonEzLYdbFiqBIIjZmY3AzcDjBkzJs7ViAwc5xw7alrZU99Kht/LuJJs8jL88S5LRBJEqgSCXcDomOejov3exjl3H9HbSZaXl+tWj5Ly2oMh/rhoO799bQs7alq7+nsMZo8u4Op5Y3jPrBFk+L1xrFJE4i1VAsESYJKZjSMSBK4Gro1vSSLxt7O2hU/+4Q1W7Kpn/rgibj5jAmVDsmjtCLFqdwOPr9jDl/66nLue2cBtF0/lXceXYmbxLltE4sCcS5wfyWb2MHAWUAzsA/4P8AM4535hZsOBCiAPCANNwDTnXIOZXQL8BPAC9zvnvnO49ysvL3c6hkBS1Y6aFq6+bxENbQF+8P5ZXDRj+NvGcc7xysYqbn98Lav3NHDh9GF89/KZDMlJj0PFItLfzGypc66812GJFAgGmgKBpKrm9iDvufsVKhvbeeimk5gxMv+Q44fCjvtf2cIPnlxHfpafX1w3l7ljCweoWhEZKIcKBIl2HQIR6QNf/edKtlQ184vr5h42DAB4PcZNZ4znX7eeSlaal2vuW8SfK3Yc9nUikjoUCERSzIvrK/n7m7u49ZxJnDKx+Khee1xpHv+65VTmjSviS39dzg+fXMdg3oooMpgoEIikkPZgiK8vWMX44mxuOXvCMU2jICuNBz56IlefOJq7n9/I//5zJaGwQoFIqkuVswxEBPhLxU62VDXz24+eSLrv2E8j9Hk93H7FTIqy0/j5C5uobe7gJ1fPfkfTFJHEpi0EIimiIxjm3hc2ccKYAs6a3OuVSY+KmfGli6bylUuP4z8r93Lzg0tpC4T6oFIRSUQKBCIp4p9v7mJXXSv/79xJfXotgRtPH88d75vJSxsquenBCoUCkRSlQCCSApxz/G7hVqYOz+2TrQMHuurEMXz/fcfzysYqPvbAElo6gn3+HiISXwoEIilg2c56Vu1u4IMnje23Kw1+oHw0P75yFos2V/PR3y6huV2hQCSVKBCIpICHFm8jK83Le2eP6Nf3uXzOKO68ajYV22r5yG9fp0mhQCRlKBCIJLmm9iALlu3mstkjyB2AuxdeNnskP716Dm9sr+P63yymoS3Q7+8pIv1PgUAkyT2zeh9tgTBXnDBqwN7z0uNLuefaOazYVc8Hf7WY2uaOAXtvEekfCgQiSe7R5bsZnpfB3DEDe++Bi2aU8ssPzWXdvkauvm8R+xvbBvT9RaRvKRCIJLH61gAvrq/k0uNL8XgG/rbF50wdxm8/ciLba1q4+peL2F3XOuA1iEjfUCAQSWJPrdpLIOR496z+PZjwUE6dWMzvb5hHZWM7H/jFQrZXt8StFhE5dgoEIknsiZV7GVWYyaxRh7+jYX8qLyvioZtOorkjyAd++Rob9zfGtR4ROXoKBCJJqi0Q4tVNVZx33LB+u/bA0Zg5Kp8/3XwyoTC8796FLNlaE++SROQoKBCIJKmFm6ppC4Q5e+rQeJfSZcrwXP7xqVMYkpPGB3+9mP+s2BPvkkTkCCkQiCSpZ9fuI9PvZf64oniX0sPooiz+9olTmDkyn0899Ab3v7Il3iWJyBFQIBBJQs45nl9byWmTisnwJ94tiQuz0/jjjfO5YNowvvnoar7x71UEQ+F4lyUih6BAIJKE1u9rYlddK+ck0O6CA2X4vfz8g3P52Knj+O2rW/nwb1/XBYxEEpgCgUgSem7tfgDOnpK4gQDA6zG+9u5pfP/9x7NkSy3vuecV1uxpiHdZItILBQKRJPTc2n1MK81jeH5GvEs5IleWj+ZPHz+J9kCYK37+Gv94c2e8SxKRAygQiCSZupYOlm6r5dzjEnvrwIHmjCnk0f93GjNG5vG5Py3j839+S3dLFEkgCgQiSealDVWEHZyV4LsLejM0L4OHbzqJz5w7iX++uYt3/fRllu2oi3dZIoICgUjSeWHdfgqz/MweXRDvUo6Jz+vhc+dP5pGbT6YjGObyn7/K7Y+vobUjFO/SRAY1BQKRJBIOO15cV8npk0rwxuFmRn1p3rgi/vPZM7jqxNH88qXNXHTXS7y2sSreZYkMWgoEIklk5e56qps7OHtqSbxL6RP5mX5uv+J4HrppPgZc++vF3PLHN9hRoxskiQy0hAoEZna/me03s5UHGW5m9lMz22hmy83shJhhITN7K/q3YOCqFhk4L6yrxAzOmJQagaDTKROK+c9nzuCz503i2bX7OPfHL3LHE2tpbAvEuzSRQSOhAgHwAHDRIYZfDEyK/t0M3BszrNU5Nzv6957+K1Ekfp5ft5/jR+YzJCc93qX0ucw0L589bzLPf+Es3jWzlHtf2MSp33uOu57ZQH2rgoFIf0uoQOCcewk41C3SLgMedBGLgAIzKx2Y6kTiq7a5g7d21CXl2QVHozQ/kx9fNZt/33oa88YN4c5n1nPaHc/xo6fWsa+hLd7liaSsfgkEZvZ5MxvZD5MeCeyIeb4z2g8gw8wqzGyRmb23H95bJK6eW7sf50iouxv2p5mj8vn1h8t59P+dxikThnD38xs59XvP8f8efpOl22pwzsW7RJGU4uun6eYCT5lZDfAn4C/OuX399F6dxjrndpnZeOA5M1vhnNt04EhmdjOR3Q2MGTOmn0sS6TuPr9jDiPwMZo3Kj3cpA2rGyHx++aFytlU38+DCbfy5Ygf/XrabiUNzuHzOSN47ZyQjCzLjXaZI0uuXLQTOuW8456YDtwClwItm9kwfTHoXMDrm+ahoP5xznY+bgReAOQep7T7nXLlzrrykJLUOzJLU1dAW4OUNVVw8sxSz5D7d8FiNHZLNV981jUVfPpfvXj6Tgkw/P3hyHafd8RxX37eQ3766RWcniLwD/bWFoNN+YC9QDfTFds4FwK1m9ggwH6h3zu0xs0KgxTnXbmbFwKnA9/vg/UQSwhMr99IRCnPJTB0yk53u49r5Y7h2/hi2V7fwjzd3sWDZLr7x79V849+rmTo8l3OmDuWk8UOYO7aQ7PT+XsyJpAbrj/1wZvYp4EqgBPgL8Gfn3OojeN3DwFlAMbAP+D/AD+Cc+4VFfhrdTeRMhBbgo865CjM7BfglECay1eMnzrnfHO79ysvLXUVFxdHPoMgAu+zuV2jpCPHU584YtFsIDmdLVTPPrtnH06v3UbGtllDY4fMYM0bmM398EbNHFTBjZD6jCjPVhjJomdlS51x5r8P6KRDcDvzJOfdWn0+8DykQSDJYuq2W9937Gt+8bDrXn1wW73KSQlN7kKXbalm8uZrFW2pYvrOOQCiyrCvI8jNjRD7TR+YxsSSHiUNzmDA0h7wMf5yrFul/hwoE/bItzTn35f6Yrshg45zjjv+spTgnjStOGBXvcpJGTrqPMyeXcObkyHFCbYEQ6/Y2smJXPSt31bNydz33v7KlKyQADM1NZ0JJDqMKMxlRkMnIgkxKCzIYUZDJ0Nx0ctJ92rIgKU0710QS2L0vbuL1rTXcfsVMcrQv/Jhl+L3MGl3ArJgbQgVDYXbUtrJxfxObKpvYuL+JzZVNvLShkv2N7Ry48TTN56EoK42i7DSG5EQeC7PSyErzkp3uIyvNG/3zkZ0eeUz3efB7PaRFH/1eI80b7fZ58HkMjxkeA48ZZgzK0OGc62pvF9uvqxu6n4GhNusPWsJIQgqGwrQGQrQFwrQHD3gMhGgPhmkLhGgLhmgPRLoj/cLd/aKPwXCYYMjFPEa6AyFHKOwIhrq7A+EwoXBkwdO5iDGzSHe0h8X0s+hCyesxfF4Pfo/h8xp+ryfSzxNZCfi8kYW/LzrczPBGVwRm3SsFr8e6Fm5r9jTw4vpK3nV8KVefOBrpWz6vh3HF2YwrzuZ8hvUY1hEMs6+hjd11reyqa6WysZ2a5g6qmzu6HrdVt1DX0kFLR4hguO92vXZ9pqIruwMDgzcaIiymv8fA00u46PxsWtdnDcJhCEfXvmHnCLvuFXI4uhIOO0c4TMw4ncO7V9Th2NdE+4djpnHgSr5z93TsSr6v26z7+9T5vLc27NluZobHc5Tj9zZ9z9GN39mu4TCEnMO5yDKosw3D0efzxg3h8+dP7rvGOgQFAuk34bCjpqWD6qYOqpsjC9Sa5sjz2pYOmtqCNLYHaWoL0tQe+WtsC9LUHqAtED7m903zekj3eUj3e6O/0GJWyF7D6+lecaf5vfgyfNGVtQevN7LS7hT5ZXLAwiz6a6VzgRYKd4aJSLgIhh3N7cFI8DgwiESHh7sWoJHXOte9cO5cOAzPz+ALF0zmE2dO0C+gAZbm8zC6KIvRRVlHNH5HMExLR5DmjhAt7d2P7aEwgWAkcAZCYTpCYQKx/cLhyP8+ZkXgnCMUu3I9YHj3X/QzFF2hxK6gu8YNHzCuiw0LPX9he6Iht7M/sSs1Iis8eHvgMIv5xe6x7sBsdIXmyCsjHZ39O590B++e40ReE9N9wHfgwDYLOw6Y1wPasJd2O6rxe5t+dBqhsCMQOoLxY6Z/YPDrDHqRYEf0B4P1bXI6DAUCeUfagyG2VDWzfl9kc+uu2sgvqt11reyub6Mj2PuKPT/TT26Gj5x0H7kZPopz0igrzu56np0W2QSb7veQ4Ys+RlfwGX7vAd0e0n3dj8l+W2BJPmk+D2m+NAqOLD+IJCQFAjli9a0BVuysZ9nOOlbsrGf9vka21bR0b2K3yIFZIwoymTEynwtnDKc0L4MhOekMyUljSHZ6dL+rH583oW6jISIy6CkQyEHtrW/jtU1VvLapmje217K5srlrWNmQLKYMz+XS40uZNCyXSUNzGFecTYbfG8eKRUTkWCkQSJeWjiAvb6jilQ1VvLqpqisAFGT5KR9byBVzRjJrdAHHjywgP0vnbIuIpBIFgkGuqqmdZ1ZHru72ysYq2oNhstK8zBtXxDUnjuHkCUOYVpqHR/vlRURSmgLBINQWCPHMmn38/Y1dvLi+klDYMbIgk2vmjeGCacMoLysizad9/CIig4kCwSCyqbKJ3722lX+8uYvGtiDD8zK46fTxvGfWCI4rzdWpbSIig5gCQYoLhx0vbqjkgVe38uL6StK8Hi49vpT3zx3FSeOH6BQ9EREBFAhSVijseHzFHn723AbW72uiJDedz58/mWvmjaEkNz3e5YmISIJRIEgxobDj38t287PnNrCpspmJQ3O486pZXDpzhI4LEBGRg1IgSCEvb6jkO4+tYe3eRqYMy+Wea0/g4hnDdYaAiIgclgJBCli/r5HvPLaGF9dXMrook59dM4dLZ5YqCIiIyBFTIEhijW0BfvDkOv6waBvZ6T7+95LjuP6UsaT7dLVAERE5OgoESeqZ1fv4yj9Xsq+xjetPGstnz5tMYXZavMsSEZEkpUCQZKqa2vn6glU8unwPU4blcu91JzBnTGG8yxIRkSSnQJAknHP8/Y1dfOux1bS0h/j8+ZP5xJkTdOaAiIj0CQWCJLCjpoX/+ccKXt5QxdyxhXzviplMGpYb77JERCSFKBAksFDY8bvXtvLDp9ZhwDfeM50PnTRWZw+IiEifUyBIUOv2NvLff1vOWzvqOGtKCd+5fCYjCzLjXZaIiKQoBYIE0x4Mcc/zm7j3hY3kZvi56+rZvGfWCN14SERE+pUCQQJZvLma//nHCjZVNvPe2SP42runU6RTCUVEZAAoECSAupYObn98LX+q2MHookwe+OiJnDVlaLzLEhGRQUSBII4CoTCPvL6dnzyzgbrWAB8/czyfPXcymWm60qCIiAyshAoEZnY/8C5gv3NuRi/DDbgLuARoAT7inHsjOuzDwFeio37bOfe7gan66IXDjqdW7+P7T65lc2Uz88cV8X/vns60EXnxLk1ERAaphAoEwAPA3cCDBxl+MTAp+jcfuBeYb2ZFwP8B5YADlprZAudcbb9XfBSa24M8unw39720mU2VzYwvyeZX15dz3nFDddCgiIjEVUIFAufcS2ZWdohRLgMedM45YJGZFZhZKXAW8LRzrgbAzJ4GLgIe7ueSD8k5x576NhZvqeaFdZU8tWofrYEQ00rzuOvq2Vw6sxSfV1caFBGR+EuoQHAERgI7Yp7vjPY7WP8Bs2DZbh5bvpu2QJj2YIjGtiDbq1tobA8CUJDl571zRnL5nJGcWFaoLQIiIpJQki0QvGNmdjNwM8CYMWP6bLq1zR1sq24h3ech3edlaG465WMLGVeczYnjipg6PA+vrjAoIiIJKtkCwS5gdMzzUdF+u4jsNojt/0JvE3DO3QfcB1BeXu76qrAPn1LGh08p66vJiYiIDKhk24G9ALjeIk4C6p1ze4AngQvMrNDMCoELov1ERETkCCTUFgIze5jIL/1iM9tJ5MwBP4Bz7hfA40ROOdxI5LTDj0aH1ZjZt4Al0Ul9s/MAQxERETk8ixywPziZWSWwrQ8nWQxU9eH0kp3aoye1Rze1RU9qj25qi576uj3GOudKehswqANBXzOzCudcebzrSBRqj57UHt3UFj2pPbqpLXoayPZItmMIREREpB8oEIiIiIgCQR+7L94FJBi1R09qj25qi57UHt3UFj0NWHvoGAIRERHRFgIRERFRIBAREREUCERERAQFAhEREUGBQERERFAgEBERERQIREREBAUCERERQYFAREREAF+8C4hlZvcD7wL2O+dm9DLcgLuAS4AW4CPOuTeiwz4MfCU66redc7873PsVFxe7srKyPqpeREQksS1durTqYLc/TqhAADwA3A08eJDhFwOTon/zgXuB+WZWBPwfUA44YKmZLXDO1R7qzcrKyqioqOij0kVERBKbmW072LCE2mXgnHsJqDnEKJcBD7qIRUCBmZUCFwJPO+dqoiHgaeCi/q+4W1sgRGNbgJaOIG2BEB3BMKGwQ/eKEBGRZJBoWwgOZySwI+b5zmi/g/UfMD9/fiM/fW5jr8PMwGuGx2N4zcjL9DE8L4Ph+RlMHJrDjBH5zBiZz6jCTCJ7RURERAZWsgWCd8zMbgZuBhgzZkyfTffMKSXkZvgJO0fIOcJhRygMIRfZShAKd/dvaA2yu76VTZXNPLtmP8FwZCvC6KJMzpo8lHOPG8ppE4vxeRNqA46IiKSwZAsEu4DRMc9HRfvtAs46oP8LvU3AOXcf0ftLl5eX99n2/Llji5g7tuioX9cWCLF+XyNv7ajjpfWV/O2Nnfx+0TaKc9K5bPYIrpk3holDc/qqTBERkV5Zou3jNrMy4NGDnGVwKXArkbMM5gM/dc7Nix5UuBQ4ITrqG8Bc59yhjkegvLzcJdpBhe3BEC+sq+Tvb+zkubX7CYQc50wdyo2nj+Pk8UO0S0FERI6ZmS11zpX3NiyhthCY2cNEfukXm9lOImcO+AGcc78AHicSBjYSOe3wo9FhNWb2LWBJdFLfPFwYSFTpPi8XTh/OhdOHU9XUzh8WbeP3C7dx7a8WM2dMAV+8cAqnTCiOd5kiIpJiEm4LwUBKxC0EvWkLhPjr0p3c/dxG9ja0cdrEYr544RRmjS6Id2kiIpJEDrWFQIEgCQJBp7ZAiD8s2sbPX9hETXMH7zthFP990RSG5mXEuzQREUkChwoEOow9iWT4vdx4+nhe+tLZfPKsCfx72W7O/uEL/PLFTXQEw/EuT0REkpgCQRLKSffx3xdN5anPncHJE4Zw+3/WcuFPXuKFdfvjXZqIiCQpBYIkVlacza8/fCIPfPREDPjIb5dw84MV7KhpiXdpIiKSZBQIUsBZU4byn8+ezpcumsLLG6o478cv8rNnN9AWCMW7NBERSRIKBCki3eflU2dN5Nn/OpPzjhvGj55ez4U/eYnn12o3goiIHJ4CQYoZUZDJPR88gT/cMB+vx/joA0u48XfajSAiIoemQJCiTptUzBOfOYPbLp7Ka5siuxHueka7EUREpHcKBCkszefhE2dOiOxGmDaMO59ZzwV3vsSza/bFuzQREUkwCgSDQGl+JvdcewJ/vHE+aT4PN/yughseWML2au1GEBGRCAWCQeTUicU8/unT+fLFU1m4uZrz7nyRO59er90IIiKiQDDYpPk8fPzMCTz3X2dx4fTh3PXsBs6/80WeXr2PwXwZaxGRwU6BYJAanp/Bz66Zw0M3zSfD5+WmByu46r5FLN5cHe/SREQkDhQIBrlTJhTz+GdO5xvvmc7Wqmauum8R1/16MW9sr413aSIiMoB0t8Mkutthf+u8m+K9L2yiurmDMyaX8PEzxnPKhCGYWbzLExGRd0i3Pz4IBYLeNbcHeXDhNn7zyhaqmtqZVprHzWeM59LjS/F7tVFJRCRZKRAchALBobUFQix4azf3vbyZjfubKM3P4Jp5Y7iyfDTD8zPiXZ6IiBwlBYKDUCA4MuGw48X1ldz/6hZe3lCFx+CcqUO5Zt4Yzpxcgk9bDUREksKhAoFvoIuR5OPxGGdPHcrZU4eyrbqZPy3ZwZ8rdvLMmgqKc9K4ZGYp7541grljCvF4dKyBiEgy0hYCbSE4JoFQmGfX7GPBst08u2Y/7cEwpfkZXDqzlHOOG0r52CLSfNpyICKSSLTL4CAUCPpGU3uQZ1bv49/LdvPShkoCIUdOuo/TJxVz9pShnDJxCKMKs+JdpojIoJc0gcDMLgLuArzAr51z3ztg+FjgfqAEqAGuc87tjA4LASuio253zr3ncO+nQND3mtqDvLqxihfW7ef5tZXsbWgDYER+BieOK+LEssjfxKE5eLV7QURkQCVFIDAzL7AeOB/YCSwBrnHOrY4Z5y/Ao86535nZOcBHnXMfig5rcs7lHM17KhD0L+cca/c2snhzNUu21vL61hoqG9sByPB7mDo8j+kj8pg+Ip/jSnMZX5JDfqY/zlWLiKSuZAkEJwNfd85dGH3+ZQDn3O0x46wCLnLO7bDIlXLqnXN50WEKBAnOOce26haWbqtl1e4GVu2uZ/WeBhrbgl3jDMlOY1xxduSvJJtRhVmU5mdQmp/BsLwMXQdBROQdSJazDEYCO2Ke7wTmHzDOMuAKIrsVLgdyzWyIc64ayDCzCiAIfM8598/+L1mOhplRVpxNWXE275sb6eecY0dNK2v3NrC1upktVc1srmzmxfWV/GXpzgNeD8U56YzIz6AkN53CrDSKsiN/hdlpDIk+FmWlkZvhIyfDR7rPG4c5FRFJPokUCI7EF4C7zewjwEvALqDz3r1jnXO7zGw88JyZrXDObTpwAmZ2M3AzwJgxYwamajkoM2PMkCzGDHn7QYdN7UF217Wyp76NvfWt7K5rY299G7uj3St3NVDT3EFHKHzQ6af5PORl+MhJ95Gb4Y8+dnZ7yUjzkuHzkuH3kuH3dD1m+r2k+zuHRfr7vYbX48HnMTwew+cxvB7Da4bXG32M9gNwDhyOzo1wsc/DzuGi/UJhR0cwTCAUpiMUfex8HnQEQmFCzpGd5mNqaS55GdqtIpKInOv5/Q5Hn3f26/reh3uOE44sHCKs88Ewg3Sfh9wB+s4nUiDYBYyOeT4q2q+Lc243kS0EmFkO8D7nXF102K7o42YzewGYA7wtEDjn7gPug8gug76eCek7Oek+Jg/LZfKw3IOO45yjuSNEbXMH1c0d1DZ3UNPcQWNbgKb2II1tQRo7H9sCNLUF2VbdEh0WoC0YWfkmC7/XePesEfzPJcdRnJMe73KkjznnaA+GaekI0dwepLkjGHlsjzxvC4ZoD0SCY/djiPaY5x3BMO3BMIFgmGA4srIJxTx2/TlHOPoYCtPVHQ7HLBata/2EmcV0R1ZYnd3EjBN9WVf/3sbtcThxzGveNr1epnvAy7qmTa+vi+0XedI5j2HnCLnofEefd7cV3d2d7eKIaS/XFeZjpxfuhzXKxTOGc+91c/t+wr1IpECwBJhkZuOIBIGrgWtjRzCzYqDGORcGvkzkjAPMrBBocc61R8c5Ffj+QBYv8WFm5KRHtgCMLjq2UxvD4chCuC0QojUQoi0Qoi0Qpi0Y6W4PRIYFoguHYNgRCocJhSEUDkefR/6C4e4FqlnMgsy6E78nphvA7/Xg93pI83nwe4206HO/z0Oa14PPG9nqUN8a4MV1lTy0eDsvb6ji9zfMY+rwvD5pR+k/7cEQ++rb2V3fyp76VvY3tFPT0h1ea5o7qG0JUNPcQXN7kOAxrFXSfZHPT7rPQ7rP2/VZ8no8eD3gtchWrc6tWGk+D16P4Yk+jzyCz+Pp+ty66C9YAKJbt4CuX7yR3jH9Okd10Re8rb+je0jsuG8fh17HcW/vd6hhMd24zrEcHuuefzPDa3S1QVd7eCL9PTHt5vEYnphxvdHnse0amV7M97yr23p87z0H9CcmbHWXG5mrY12uHYt+CQRm9vlDDXfO/biXfkEzuxV4kshph/c751aZ2TeBCufcAuAs4HYzc0R2GdwSfflxwC/NLEzkls7fiz07QeRQPB4jM81LZpqXwngXcxhnTxnK1fNG85H7l/DBXy3mn7ecOqALDOldY1ug6/iXzZVNbK5qZntNC7vr2qhqan/b+Gk+T+SYl+hxMCMLsyjM8pOb4SMrLRJws9K8kcd0HznpXrLSfNFdWZGgmBZd+fu9pruRSp/ol7MMzOwh4ERgQbTXu4HXgQ0Azrlv9PmbHgOdZSDJalNlE++951VGFWbxz1tO0cGTA2h/Yxsrd9WzclcDK3fVs2p3A7vqWruGeyzyq27skGxGFmRQmp/J8PwMRuRnUloQOVsmO82rlbjERTzOMhgFnOCca4wW8HXgMefcdf30fiKDyoSSHO68cjY3PljBz57dyBcunBLvklKSc47NVc28vqWm6y925T++OJsTxhZy7fwxTCjJYUJJNmOGZCmgSVLqr0AwDOiIed4R7ScifeS8acN4/9xR3PviJt41q1THE/SR+tYAL62v5Pl1+3lpfVXXJv/inHTmjyviY6eNY+bIyMW0Burob5GB0F+B4EHgdTP7R/T5e4EH+um9RAatr1x6HE+v3se3H13D72+Yp83Qx2hPfSuPLd/DU6v3sXRbLaGwoyDLzxmTSjhlwhDmjStiXHG22ldSWr8EAufcd8zsP8Dp0V4fdc692R/vJTKYFWSl8dnzJvGNf6/m2TX7OW+aNsQdqeqmdh5fuZd/L9vNkq01OAfHlebxiTPHc87UocweXaj7bcig0m+nHTrn3gDe6K/pi0jEdSeN5fcLt/HDp9ZxztSheLQSO6hw2LFoczV/fH07T63aSyDkmDQ0h8+fN5l3zRrBuOLseJcoEjeJdB0CETkGfq+Hz5w3ic888hZPrNrLJTNL411SwmloC/Cn13fwx8Xb2FrdQkGWn+tPLuP9c0cxdXiudgWIoEAgkhLedfwI7np2A3c9s4GLpg/XVoKofQ1t3P/qFh5atJ3G9iAnlhXymfMmcfGMUjL8OhNAJJYCgUgK8HqMz5wb2Urwn5V7ufT4wb2VYFddK3c/t4G/Ld1FMBzmkpmlfPyMCcwclR/v0kQSlgKBSIro3Erws+c2cPGMwbmVoLqpnXue38QfFm0D4MoTR3HT6eMZO0THBogcjgKBSIrweoxbz57I5/+8jKfX7OPC6cPjXdKAae0I8YsXN/HrlzfTGgjx/rmj+PS5kxhVqMs6ixwpBQKRFPKeWd1bCS6YNizlD5ZzzvH4ir1857HV7K5v45KZw/n8+VOYODQn3qWJJB0FApEU4vN6uOWsiXzpb8t5YV0lZ08dGu+S+s36fY18fcEqXttUzXGlefzk6jnMG1cU77JEkpYn3gWISN+6/ISRjCzI5KfPbej1drLJrj0Y4sdPreOSu15m1e4GvnXZdP5966kKAyLvkLYQiKQYv9fDJ8+awFf+uZJXN1Zz2qTieJfUZ97aUceX/rqM9fuauHzOSL76rmkUZafFuyyRlKAtBCIp6APloxiel8FPn9sQ71L6RFsgxO2Pr+GKn79KY1uQ337kRO68arbCgEgfUiAQSUHpPi8fP3M8r2+pYfHm6niX845srmzi8p+/xi9f2sxVJ47hyc+dkdLHRojEiwKBSIq6Zt4YinPS+dlzG+NdyjFbsGw37/7ZK+ypb+X+j5Rz+xUzydMth0X6hQKBSIrK8Hu5+YxxvLKxiqXbauNdzlFpC4T433+s4NMPv8nU0jwe//TpnDNVd3IU6U8KBCIp7IPzx1KY5ednSXQsQecugj8u3s7HzxzPIzefxIiCzHiXJZLyFAhEUlh2uo8bTx/PC+sqWb6zLt7lHNaBuwi+fPFx+L1aTIkMBH3TRFLc9SePJS/Dl9DHEmgXgUj8JVQgMLOLzGydmW00s9t6GT7WzJ41s+Vm9oKZjYoZ9mEz2xD9+/DAVi6SuHIz/Hz01HE8vXofy3bUxbuct9EuApHEkDCBwMy8wD3AxcA04Bozm3bAaD8EHnTOHQ98E7g9+toi4P+A+cA84P/MrHCgahdJdDeePo7inHS+/u9VhMOJc/VC7SIQSRyJ9M2bB2x0zm12znUAjwCXHTDONOC5aPfzMcMvBJ52ztU452qBp4GLBqBmkaSQm+HnSxdN4c3tdfxr2a54l6NdBCIJKJECwUhgR8zzndF+sZYBV0S7LwdyzWzIEb5WZFB7/wmjmDUqn9sfX0tTezBudWgXgUhiSqRAcCS+AJxpZm8CZwK7gNDRTMDMbjazCjOrqKys7I8aRRKSx2N8/T3TqWxq57uPr4lLDdpFIJK4EumbuAsYHfN8VLRfF+fcbufcFc65OcD/RvvVHclrY6Zxn3Ou3DlXXlJS0ofliyS+OWMKuen08Ty0eDsvrR+4QNzaEeK2vy3XLgKRBJZIgWAJMMnMxplZGnA1sCB2BDMrNrPOmr8M3B/tfhK4wMwKowcTXhDtJyIH+Pz5k5lQks2X/rqc6qb2fn+/tXsbePfdr/Cnih188qwJ2kUgkqASJhA454LArURW5GuAPzvnVpnZN83sPdHRzgLWmdl6YBjwnehra4BvEQkVS4BvRvuJyAEy/F7uunoOtS0d/L+H3yQYCvfL+zjn+P2ibbzn7lepbw3w+4/N578vmqpdBCIJypxLnFOQBlp5ebmrqKiIdxkicfGXih188a/Luf7ksXzjPdMxsz6b9q66Vm7723Je3lDFmZNL+NGVsyjOSe+z6YvIsTGzpc658t6G+Qa6GBFJDB8oH82G/U3c99JmcjN8fOGCKe84FITDjode387tj6/BAd+6bDofnD8Wj6fvwoaI9A8FApFB7MsXT6WxLcg9z2+isS3I1941Dd8xbtKv2FrDNx9dzfKd9Zw2sZjbr5jJ6KKsPq5YRPqLAoHIIGZmfOe9M8jL8PHLlzazbm8jP3j/LMYMOfIV+Yqd9dzz/EaeWLWX4XkZ/PjKWVw+Z2Sf7oIQkf6nQCAyyHk8xpcvOY7Jw3L5+oJVXPCTF/nQSWP5yKnjGHmQswEa2wI8tWoff126k4Wbq8lN9/GZcyfx8TPHk5WmxYpIMtJBhTqoUKTL7rpWfvDkOv711i7CDqaPyGPGiHyKc9MwjP2Nbazd28iq3Q2Ewo7RRZl8cP5Yrp0/hrwMf7zLF5HDONRBhQoECgQib7OjpoVHl+/hpfWVrN/XSH1rgJBzlOSkM644mxPLijh7agknjCnUrgGRJKJAcBAKBCJHzjmnlb9IkjtUINAVQkTkiCgMiKQ2BQIREREZ3LsMzKwS2NaHkywGqvpweslO7dGT2qOb2qIntUc3tUVPfd0eY51zvd7Zb1AHgr5mZhUH2zczGKk9elJ7dFNb9KT26Ka26Gkg20O7DERERESBQERERBQI+tp98S4gwag9elJ7dFNb9KT26Ka26GnA2kPHEIiIiIi2EIiIiIgCQZ8xs4vMbJ2ZbTSz2+JdT38xs61mtsLM3jKzimi/IjN72sw2RB8Lo/3NzH4abZPlZnZCzHQ+HB1/g5l9OF7zc7TM7H4z229mK2P69dn8m9ncaPtujL42oa8GdJD2+LqZ7Yp+Rt4ys0tihn05Om/rzOzCmP69fn/MbJyZLY72/5OZpQ3c3B0dMxttZs+b2WozW2Vmn4n2H3Sfj0O0xWD9bGSY2etmtizaHt+I9u91HswsPfp8Y3R4Wcy0jqqdjopzTn/v8A/wApuA8UAasAyYFu+6+mletwLFB/T7PnBbtPs24I5o9yXAfwADTgIWR/sXAZujj4XR7sJ4z9sRzv8ZwAnAyv6Yf+D16LgWfe3F8Z7nY2iPrwNf6GXcadHvRjowLvqd8R7q+wP8Gbg62v0L4JPxnudDtEUpcEK0OxdYH53nQff5OERbDNbPhgE50W4/sDj6f+x1HoBPAb+Idl8N/OlY2+lo/rSFoG/MAzY65zY75zqAR4DL4lzTQLoM+F20+3fAe2P6P+giFgEFZlYKXAg87Zyrcc7VAk8DFw1wzcfEOfcSUHNA7z6Z/+iwPOfcIhf59j8YM62EdJD2OJjLgEecc+3OuS3ARiLfnV6/P9Ffv+cAf42+PrZtE45zbo9z7o1odyOwBhjJIPx8HKItDibVPxvOOdcUfeqP/jkOPg+xn5m/AudG5/mo2ulo61Qg6BsjgR0xz3dy6A9/MnPAU2a21MxujvYb5pzbE+3eCwyLdh+sXVKtvfpq/kdGuw/sn4xujW4Gv79zEzlH3x5DgDrnXPCA/gkvuol3DpFfgoP683FAW8Ag/WyYmdfM3gL2Ewl5mzj4PHTNd3R4PZF57tdlqgKBHK3TnHMnABcDt5jZGbEDo79cBu2pK4N9/qPuBSYAs4E9wI/iWs0AM7Mc4G/AZ51zDbHDBtvno5e2GLSfDedcyDk3GxhF5Bf91PhW9HYKBH1jFzA65vmoaL+U45zbFX3cD/yDyAd7X3RzJtHH/dHRD9YuqdZefTX/u6LdB/ZPKs65fdGFXxj4FZHPCBx9e1QT2YzuO6B/wjIzP5EV4B+dc3+P9h6Un4/e2mIwfzY6OefqgOeBkzn4PHTNd3R4PpF57tdlqgJB31gCTIoeMZpG5CCQBXGuqc+ZWbaZ5XZ2AxcAK4nMa+eR0B8G/hXtXgBcHz2a+iSgPrrp9EngAjMrjG4yvCDaL1n1yfxHhzWY2UnR/YXXx0wraXSu/KIuJ/IZgUh7XB09gnocMInIQXK9fn+iv6afB94ffX1s2yac6P/sN8Aa59yPYwYNus/HwdpiEH82SsysINqdCZxP5LiKg81D7Gfm/cBz0Xk+qnY66kLfyZGT+utxFOklRI6k3QT8b7zr6ad5HE/k6NVlwKrO+SSyb+tZYAPwDFAU7W/APdE2WQGUx0zrY0QOiNkIfDTe83YUbfAwkU2dASL76W7oy/kHyoksJDcBdxO9eFii/h2kPX4fnd/l0YVSacz4/xudt3XEHCF/sO9P9DP3erSd/gKkx3ueD9EWpxHZHbAceCv6d8lg/Hwcoi0G62fjeODN6HyvBL52qHkAMqLPN0aHjz/WdjqaP12pUERERLTLQERERBQIREREBAUCERERQYFAREREUCAQERERFAhEREQEBQIRERFBgUBERESA/w81+A+3dYK7tAAAAABJRU5ErkJggg==\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "fig, axs = plt.subplots(nrows=6, sharex=True, figsize=(8,10))\n", + "axs[0].plot(raw_t,raw_d[:,0,0,0])\n", + "axs[0].set(ylabel='r')\n", + "for i, (ax, label) in enumerate( zip(axs[1:], ['BOLD', 's', 'f', 'v', 'q'])):\n", + " ax.plot(b_t[:], b_d[:,i,0,0])\n", + " ax.set(ylabel=label)" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.8.10" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}