{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Simulating morphological Golgi cell with ``braincell``\n",
    "\n",
    "[![Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/chaobrain/brainx/blob/main/docs/getting_started/braincell_morphological_golgi_cell.ipynb)\n",
    "[![Open in Kaggle](https://kaggle.com/static/images/open-in-kaggle.svg)](https://kaggle.com/kernels/welcome?src=https://github.com/chaobrain/brainx/blob/main/docs/getting_started/braincell_morphological_golgi_cell.ipynb)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 91,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-09-16T07:08:20.871692Z",
     "start_time": "2025-09-16T07:08:17.627721Z"
    }
   },
   "outputs": [],
   "source": [
    "import os.path\n",
    "\n",
    "import braincell\n",
    "import brainstate\n",
    "import braintools\n",
    "import brainunit as u\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 92,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-09-16T07:08:20.880852Z",
     "start_time": "2025-09-16T07:08:20.876706Z"
    }
   },
   "outputs": [],
   "source": [
    "brainstate.environ.set(dt=0.01 * u.ms, precision=64)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Golgi cell model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 93,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-09-16T07:08:20.904487Z",
     "start_time": "2025-09-16T07:08:20.897872Z"
    }
   },
   "outputs": [],
   "source": [
    "class Golgi(braincell.MultiCompartment):\n",
    "    def __init__(self, popsize, morphology, el, gl, gh1, gh2, ek, gkv11, gkv34, gkv43, ena, gnarsg, V_init=-65):\n",
    "        super().__init__(\n",
    "            popsize=popsize,\n",
    "            morphology=morphology,\n",
    "            V_th=20. * u.mV,\n",
    "            V_initializer=braintools.init.Constant(V_init * u.mV),\n",
    "            spk_fun=braintools.surrogate.ReluGrad(),\n",
    "            solver='staggered'\n",
    "        )\n",
    "\n",
    "        self.IL = braincell.channel.IL(self.varshape, E=el * u.mV, g_max=gl * u.mS / (u.cm ** 2))\n",
    "        self.Ih1 = braincell.channel.Ih1_Ma2020(self.varshape, E=-20. * u.mV, g_max=gh1 * u.mS / (u.cm ** 2))\n",
    "        self.Ih2 = braincell.channel.Ih2_Ma2020(self.varshape, E=-20. * u.mV, g_max=gh2 * u.mS / (u.cm ** 2))\n",
    "\n",
    "        self.k = braincell.ion.PotassiumFixed(self.varshape, E=ek * u.mV)\n",
    "        self.k.add(IKv11=braincell.channel.IKv11_Ak2007(self.varshape, g_max=gkv11 * u.mS / (u.cm ** 2)))\n",
    "        self.k.add(IKv34=braincell.channel.IKv34_Ma2020(self.varshape, g_max=gkv34 * u.mS / (u.cm ** 2)))\n",
    "        self.k.add(IKv43=braincell.channel.IKv43_Ma2020(self.varshape, g_max=gkv43 * u.mS / (u.cm ** 2)))\n",
    "\n",
    "        self.na = braincell.ion.SodiumFixed(self.varshape, E=ena * u.mV)\n",
    "        self.na.add(INa_Rsg=braincell.channel.INa_Rsg(self.varshape, g_max=gnarsg * u.mS / (u.cm ** 2)))\n",
    "\n",
    "    def step_run(self, t, inp):\n",
    "        with brainstate.environ.context(t=t):\n",
    "            self.update(inp)\n",
    "            return self.V.value"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Utility functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-09-16T07:08:20.914787Z",
     "start_time": "2025-09-16T07:08:20.911500Z"
    }
   },
   "outputs": [],
   "source": [
    "def step_input(num, dur, amp, dt):\n",
    "    brainstate.environ.set(dt=dt * u.ms)\n",
    "    value = u.math.zeros((len(dur), num))\n",
    "    for i in range(len(value)):\n",
    "        value = value.at[i, 0].set(amp[i])\n",
    "    return braintools.input.section(values=value, durations=dur * u.ms) * u.nA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-09-16T07:08:20.931535Z",
     "start_time": "2025-09-16T07:08:20.927313Z"
    }
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "\n",
    "def plot_voltage_comparison(\n",
    "    t_vec,\n",
    "    v_vecs,\n",
    "    indices=None,\n",
    "    title=\"Voltage Comparison\",\n",
    "    xlabel=\"Time (ms)\",\n",
    "    ylabel=\"Voltage (mV)\",\n",
    "    figsize=(7, 4),\n",
    "    legend=True,\n",
    "    grid=True\n",
    "):\n",
    "    if indices is None:\n",
    "        indices = range(len(v_vecs))\n",
    "    plt.figure(figsize=figsize)\n",
    "    for idx, i in enumerate(indices):\n",
    "        plt.plot(t_vec, v_vecs[i], linewidth=2)\n",
    "    if legend:\n",
    "        plt.legend(frameon=False, fontsize=11)\n",
    "    plt.xlabel(xlabel, fontsize=13)\n",
    "    plt.ylabel(ylabel, fontsize=13)\n",
    "    plt.title(title, fontsize=15, weight='bold')\n",
    "    if grid:\n",
    "        plt.grid(alpha=0.4)\n",
    "    plt.tight_layout()\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-09-16T07:08:20.947338Z",
     "start_time": "2025-09-16T07:08:20.939424Z"
    }
   },
   "outputs": [],
   "source": [
    "# index for ion channel\n",
    "def is_basal(idx):\n",
    "    return (\n",
    "        0 <= idx <= 3\n",
    "        or 16 <= idx <= 17\n",
    "        or 33 <= idx <= 41\n",
    "        or idx == 84\n",
    "        or 105 <= idx <= 150\n",
    "    )\n",
    "\n",
    "\n",
    "def is_apical(idx):\n",
    "    return (\n",
    "        4 <= idx <= 15\n",
    "        or 18 <= idx <= 32\n",
    "        or 42 <= idx <= 83\n",
    "        or 85 <= idx <= 104\n",
    "    )\n",
    "\n",
    "\n",
    "def seg_ion_params(morphology):\n",
    "    # segment index for each type\n",
    "    index_soma = []\n",
    "    index_axon = []\n",
    "    index_dend_basal = []\n",
    "    index_dend_apical = []\n",
    "\n",
    "    for i, seg in enumerate(morphology.segments):\n",
    "        name = str(seg.section_name)\n",
    "        if name.startswith(\"soma\"):\n",
    "            index_soma.append(i)\n",
    "        elif name.startswith(\"axon\"):\n",
    "            index_axon.append(i)\n",
    "        elif name.startswith(\"dend_\"):\n",
    "            idx = int(name.split(\"_\")[-1])\n",
    "            if is_basal(idx):\n",
    "                index_dend_basal.append(i)\n",
    "            if is_apical(idx):\n",
    "                index_dend_apical.append(i)\n",
    "\n",
    "    n_compartments = len(morphology.segments)\n",
    "\n",
    "    # conductvalues\n",
    "    conductvalues = 1e3 * np.array(\n",
    "        [\n",
    "            0.00499506303209, 0.01016375552607, 0.00247172479141, 0.00128859564935,\n",
    "            3.690771983E-05, 0.0080938853146, 0.01226052748146, 0.01650689958385,\n",
    "            0.00139885617712, 0.14927733727426, 0.00549507510519, 0.14910988921938,\n",
    "            0.00406420380423, 0.01764345789036, 0.10177335775222, 0.0087689418803,\n",
    "            3.407734319E-05, 0.0003371456442, 0.00030643090764, 0.17233663543619,\n",
    "            0.00024381226198, 0.10008178886943, 0.00595046001148, 0.0115, 0.0091\n",
    "        ]\n",
    "    )\n",
    "\n",
    "    # IL\n",
    "    gl = np.ones(n_compartments)\n",
    "    gl[index_soma] = 0.03\n",
    "    gl[index_axon] = 0.001\n",
    "    gl[index_axon[0:5]] = 0.03\n",
    "    gl[index_dend_basal] = 0.03\n",
    "    gl[index_dend_apical] = 0.03\n",
    "\n",
    "    # IKv11_Ak2007\n",
    "    gkv11 = np.zeros(n_compartments)\n",
    "    gkv11[index_soma] = conductvalues[10]\n",
    "\n",
    "    # IKv34_Ma2020\n",
    "    gkv34 = np.zeros(n_compartments)\n",
    "    gkv34[index_soma] = conductvalues[11]\n",
    "    gkv34[index_axon[5:]] = 9.1\n",
    "\n",
    "    # IKv43_Ma2020\n",
    "    gkv43 = np.zeros(n_compartments)\n",
    "    gkv43[index_soma] = conductvalues[12]\n",
    "\n",
    "    # ICaGrc_Ma2020\n",
    "    gcagrc = np.zeros(n_compartments)\n",
    "    gcagrc[index_soma] = conductvalues[15]\n",
    "    gcagrc[index_dend_basal] = conductvalues[8]\n",
    "    gcagrc[index_axon[0:5]] = conductvalues[22]\n",
    "\n",
    "    # ICav23_Ma2020\n",
    "    gcav23 = np.zeros(n_compartments)\n",
    "    gcav23[index_dend_apical] = conductvalues[3]\n",
    "\n",
    "    # ICav31_Ma2020\n",
    "    gcav31 = np.zeros(n_compartments)\n",
    "    gcav31[index_soma] = conductvalues[16]\n",
    "    gcav31[index_dend_apical] = conductvalues[4]\n",
    "\n",
    "    # INa_Rsg\n",
    "    gnarsg = np.zeros(n_compartments)\n",
    "    gnarsg[index_soma] = conductvalues[9]\n",
    "    gnarsg[index_dend_apical] = conductvalues[0]\n",
    "    gnarsg[index_dend_basal] = conductvalues[5]\n",
    "    gnarsg[index_axon[0:5]] = conductvalues[19]\n",
    "    gnarsg[index_axon[5:]] = 11.5\n",
    "\n",
    "    # Ih1_Ma2020\n",
    "    gh1 = np.zeros(n_compartments)\n",
    "    gh1[index_axon[0:5]] = conductvalues[17]\n",
    "\n",
    "    # Ih2_Ma2020\n",
    "    gh2 = np.zeros(n_compartments)\n",
    "    gh2[index_axon[0:5]] = conductvalues[18]\n",
    "\n",
    "    # IKca3_1_Ma2020\n",
    "    gkca31 = np.zeros(n_compartments)\n",
    "    gkca31[index_soma] = conductvalues[14]\n",
    "\n",
    "    return gl, gh1, gh2, gkv11, gkv34, gkv43, gnarsg, gcagrc, gcav23, gcav31, gkca31"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Morphology"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "metadata": {},
   "outputs": [],
   "source": [
    "import kagglehub\n",
    "\n",
    "# Download latest version\n",
    "path = kagglehub.dataset_download(\"oujago/golgi-cell-morphology-example1\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-09-16T07:08:21.599369Z",
     "start_time": "2025-09-16T07:08:20.953345Z"
    }
   },
   "outputs": [],
   "source": [
    "Golgi_mor = braincell.Morphology.from_asc(os.path.join(path, 'golgi.asc'))\n",
    "Golgi_mor.set_passive_params()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 99,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-09-16T07:08:22.207854Z",
     "start_time": "2025-09-16T07:08:21.607376Z"
    }
   },
   "outputs": [],
   "source": [
    "Golgi_mor.visualize()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 100,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-09-16T07:08:22.372817Z",
     "start_time": "2025-09-16T07:08:22.233390Z"
    }
   },
   "outputs": [],
   "source": [
    "gl, gh1, gh2, gkv11, gkv34, gkv43, gnarsg, gcagrc, gcav23, gcav31, gkca31 = seg_ion_params(Golgi_mor)\n",
    "nseg = len(Golgi_mor.segments)\n",
    "El = -55\n",
    "Ek = -80\n",
    "Ena = 60\n",
    "V_init = -65 * np.ones(nseg)\n",
    "cell = Golgi(\n",
    "    popsize=1,  # number of cells in the population,\n",
    "    morphology=Golgi_mor,\n",
    "    el=El,\n",
    "    gl=gl,\n",
    "    gh1=gh1,\n",
    "    gh2=gh2,\n",
    "    ek=Ek,\n",
    "    gkv11=gkv11,\n",
    "    gkv34=gkv34,\n",
    "    gkv43=gkv43,\n",
    "    ena=Ena,\n",
    "    gnarsg=gnarsg,\n",
    "    V_init=V_init,\n",
    ")\n",
    "\n",
    "_ = brainstate.nn.init_all_states(cell)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 101,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-09-16T07:08:26.909178Z",
     "start_time": "2025-09-16T07:08:22.387826Z"
    }
   },
   "outputs": [],
   "source": [
    "I = step_input(num=nseg, dur=[100, 0, 0], amp=[0, 0, 0], dt=brainstate.environ.get_dt() / u.ms)\n",
    "times = u.math.arange(I.shape[0]) * brainstate.environ.get_dt()\n",
    "vs = brainstate.transform.for_loop(cell.step_run, times, I)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 102,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2025-09-16T07:08:36.965217Z",
     "start_time": "2025-09-16T07:08:36.741760Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAysAAAGbCAYAAADEAg8AAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjYsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvq6yFwwAAAAlwSFlzAAAPYQAAD2EBqD+naQAAZP9JREFUeJzt3QecXGW9//Hv9O0lbZOQTgsQIBAEgjQBCQgogoCISImiXKp4Ubj8xcJVEAQLIIpXinTx2kBaLl0IhCIdQgtpm93U7X3m/F/Pc2ZmZze7m93N7p7ZOZ/36zWemTOzu2fjsGd+51eegOM4jgAAAAAgywS9PgAAAAAA6AnBCgAAAICsRLACAAAAICsRrAAAAADISgQrAAAAALISwQoAAACArESwAgAAACArhb0+gGyTSCRUWVmp4uJiBQIBrw8HAAAAyDlmqcf6+npNnjxZwWDv+ROClW5MoDJ16lSvDwMAAADIeStXrtSUKVN6fZ5gpRuTUUn9w5WUlHh9OAAAAEDOqaurswmC1Gfv3hCsdJMq/TKBCsEKAAAAMHy21HZBgz0AAACArESwAgAAACArEawAAAAAyEr0rAAAAABDKB6Pq729XX4WiUQUCoW2+vsQrAAAAABDtHZIVVWVampqvD6UrFBWVqaJEydu1dqFBCsAAADAEEgFKhMmTFBBQYFvFxh3HEdNTU1au3atfTxp0qRBfy+CFQAAAGAISr9SgcrYsWPld/n5+XZrAhbzbzLYkjAa7AEAAICtlOpRMRkVqMu/xdb07xCsAAAAAEPEr6Vfw/VvQbACAAAAICsRrAAAAADISgQrAAAAALISwQoAAEA2euS/pHu+IrU3e30kyGHHHHOMjjjiiB6fe/bZZ23fyRtvvKEVK1boqKOOsk3zZrrXxRdfrI6OjmE/PkYXAwAAZJualdILN7r3339U2uVYr48IOWrhwoU6/vjjtWrVKk2ZMqXLc7feeqv22msv7bLLLpo7d65d4PH555/XmjVr9LWvfc2uUv/Tn/50WI+PzAoAAEC22fhR5/11S708Emzt4ohtHZ7cHMfp1zEeffTRGj9+vG677bYu+xsaGnT//ffbYOaxxx7TO++8ozvvvNMGLUceeaSuuOIK3XjjjWpra9NwIrMCAACQbRrXd96vWe7lkWArNLfHtfPlj3rys9/58QIVRLf8UT8cDtssiQlWLrvssvS4YROomIUuTz75ZF177bXaddddVVFRkf66BQsW6Oyzz9bbb7+tPfbYY9h+DzIrAAAA2aZhbcb9ai+PBD5w5pln6qOPPtLTTz/dpQTMlIeVlpaqqqqqS6BipB6b54YTmRUAAIBs01rXeb9pg5dHgq2QHwnZDIdXP7u/Zs+erf3220+33HKLDj74YH344Ye2uf7HP/6xvEawAgAAkG0yJ4A1bfTySLAVTElVf0qxssHChQt13nnn2T4Uk1XZdtttddBBB9nnTGP9kiVLury+uro6/dxwogwMAAAg23S0dN4nWMEIOPHEExUMBnX33Xfrj3/8oy0NS/WvzJ8/X2+++abWru0sT1y0aJFKSkq08847D+txEawA2WDN69LPZkr//I7XRwIAyLbMSlu91DG8E5eAoqIinXTSSbr00kvtaOLTTz89/dzhhx9ug5JTTz1Vr7/+uh599FH9v//3/3TOOecoFosN63ERrADZ4IXfSs0bpZf+R2qu8fpoAADZlFkxWjg3YGRKwTZt2mQnfU2ePDm9PxQK6cEHH7Rbk2X56le/aieIjURPy+googNyXfVbnfer3pBmHujl0QAAvNZ91fq2BkkTvDoa+MT8+fN7XZ9l+vTpeuihh0b8mMisANmgfk3n/U2feHkkAIBszKy0mmAF8B+CFcBr8Y5ui3+t9PJoAADZGKy0NXp1JICnCFYArzWZQCUj5VpLsAIAvpdI9FAGBvgPwQrgteZNva9aDADwJ6dbsNJa79WRAJ4atcHKVVddZWc/X3jhhel9LS0tdoTa2LFj7fi1448/Pr1gDZC12pq6PmalYgBA92CFMjD41KgMVl566SX97ne/02677dZl/7e//W098MADuv/++/X000+rsrJSxx13nGfHCfRL99Q+wQoAIBWshPPcLWVg8KlRF6w0NDTolFNO0e9//3uVl5en99fW1uoPf/iDrrvuOh1yyCGaN2+ebr31Vj3//PN64YUXPD1moE/tycxKwVh3S7ACAEj1MsaK3S3TwOBToy5YMWVeRx11lA477LAu+1955RW1t7d32T979mxNmzZNixcv7vX7tba2qq6urssNGFGp1H7p1M7gpXtpGADAn5mVWIm7JbMCnxpVi0Lee++9evXVV20ZWHdVVVWKRqMqKyvrsr+iosI+15srr7xSP/rRj4bleIEBZVaKJ0rVESnR7mZXogVeHxkAwOtgJY9gBf42ajIrK1eu1AUXXKC77rpLeXnJ+s0hcOmll9oSstTN/BzAk8xKpEDKTwbbLbWeHhIAIEuClWiRuyXjjmFyzDHH6IgjjujxuWeffdYOtHrjjTd0/vnn2zaLWCymuXPnaqSMmmDFlHmtXbtWe+65p8LhsL2ZJvpf//rX9r7JoLS1tammpqbL15lpYBMnTuz1+5p/8JKSki43wJPMismkpGuTGVEJAL7mOJ0XsoyOZk8PB7lr4cKFWrRokVatWrXZc6b/e6+99koPtTrzzDN10kknjejxjZpg5dBDD9Wbb76p1157LX0z/3im2T51PxKJ6PHHH09/zdKlS7VixQrNnz/f02MH+pSIu9tguLM2uZXeKQDwtVSwEi10t+3dVrQHhsjRRx+t8ePH67bbbttsqJWZsGuCGcMkCEzv+KxZszSSRk3PSnFxsebMmdNlX2FhoV1TJbXf/GNedNFFGjNmjM2QnHfeeTZQ2XfffT06amAAqf5AqLM2mcwKAPhbugysoGsWHqMv6PTq/7tIgRQIbPFlpkLpa1/7mg1WLrvsMlv2ZZhAJR6P6+STTx6Bg+3j+JRDfvGLXygYDNrFIM2UrwULFug3v/mN14cF9DNYCXZmVuhZAQB/S50b0mVgZFZGJROo/HSyNz/7vyo7M3NbYMq7rrnmGtticfDBB6dLwMxn6tLSUnlpVAcrTz31VJfHpvH+xhtvtDdgVAcrlIEBgM91LwOjZwXDxyz3sd9+++mWW26xwcqHH35om+t//OMfy2ujOlgBcqpnxQYrNNgDADIzKwQro5rJjJkMh1c/ewBMO4VpoTAX/U1WZdttt9VBBx0krxGsANlyQgqGpEhyRGULmRUA8LXuPSuUgY1Opv+jn6VYXjvxxBPtMiF33323/vjHP+rss89O9694iWAFyJoysACZFQBAl3PDL59ZrQvNHTIrGGZFRUV2LLFZg7Curk6nn356l+dNaZiZEGYWW29ubrbTeI2dd97ZLsw+XAhWAK/RswIA6GV08bJaRzKfAwlWMAJMKdgf/vAHfe5zn9PkyV0HA3z961+3Dfgpe+yxh90uW7ZMM2bMGLZjIlgBsipYIbMCAJASTsIuhtesWOeikCaAyYKyHOSu+fPny0mt8bOFwVYjZdQsCgnkerDy2qo6PfJBMqPS1ujtMQEAPJVIDl9pUl7nuSLe7u1BAR4gswJkSbDy9Icb9UJivY4w6X6CFQDwNSfhnhuanGRmJZVdCQ9fbwCQjcisAFkSrDgKdJ6UWKkYAHwtVYrTprDiqY9r9K3AhwhWAK8lU/0JJ9iZ7iezAgD+lsysmM6VVtthT7ACfyJYAbIks5JQoLORkswKAPiakz43BNXiRNydrLUCHyJYAbyWcUJKl4GZE1JqZXsAgI9LhE2TPReyRpNEMisGDcm/BQ32gNeSdckms9KYKgNLlYLlJdddAQD4smfFloGZzIqZWNxOZiWbmYURg8GgKisrNX78ePs4G1aA9+r929bWpnXr1tl/k61ZNJJgBfCak+xZUUCtitht0FxLM1fQCFYAQH4vEW5J9ayYaWDIWuZD+cyZM7VmzRobsEAqKCjQtGnT7L/NYBGsAFlUBmYunTU7MRUGWmiyBwAfC2RMikwHKzTYZz2TQTAfzjs6OhSP+7ucOxQKKRwOb3V2iWAFyKpgxa1NLlQLtckA4GsZwYqTClYoAxsNzIfzSCRib9h6NNgDWZTqNxqd1PhighUA8K2MfsbOMjCCFfgPwQrgMSe1zkoyWEmPL25r8PKwAAAeCigzWGF0MfyLYAXwWCIdrARVXhBhRCUAIKNnJWNRSIIV+BDBCpAlM8hNXfL44ljnWiuUgQGA/F4GZjbpRSHpWYEPEawAHkvEOzMrZfnRjFXsmQYGAH7VWQZGZgX+RrACZEmwEgwFVRgLdS4MSWYFAHwrkJwG5q7BRbAC/yJYAbKkZyUUDKsgFrbrrFisswIA/pUxDcwsGGwRrMCHCFYAjyWSTZShUFBF0XBGgz3BCgD4PbNi7rHOCvyMYAXIktHFwVBIBbFQZ7BCGRgA+FZQPWVWWMEe/kOwAnjMSfashEJhFUZNGViyZ4XRxQDg6xIwIxYNZywK2erdMQEeIVgBsqUMLBi0mZXm1EmpnStoAOBLyfOCkR+NZIwu5rwA/yFYATzmJNdZCYVCNrPSeQWN2mQA8HuwUhCLZEwDI7MC/yFYAbKmZyWswlg4o5GSMjAA8HsZWF7EXMSiZwX+RbACeM1JBivBoAqjmWVgZFYAwJe6lYGRWYGfEawA2XIFLWCmgYU7T0pkVgBAfg9W8mzGnZ4V+BfBCpAlJ6V0ZiW1KCQ9KwDgT5k9K2RW4HMEK0CWnJScYMjtWUnVJnMFDQB8KqNnpcvgFc4L8B+CFcBrqQb7QFAFtmcltYI9JyUA8HtmJRYx5cGpBnsyK/AfghXAY4HUSSkYVH4klL6C5hCsAID8HqxEIplTIjkvwH8IVoAsSfebnpU827OSke7PGF8JAPDj6OJQZ2bFTI+Mt3t3XIAHCFaALLmCFgiElBc2J6VoZ8Yl3ubxwQEAPC0DC0c6e1YMhq/AZwhWAI8FkuusBIIhRUIBtQXzOp8k5Q8Avg1WEk7ANtinMysGa3DBZwhWgGzJrASDCgQCCkeiijsB9zmuoAGAf4MVBZQXMR/VAmpNrbXCeQE+Q7ACeCyQrE02wUqqPrlzIhgLQwKAn4OVWDhk76fH2hOswGcIVoCsWRQylA5W0vXJpPsBwL/rbymomM2sKGNhSM4L8BeCFcBjAaVGF7vBSub4YnpWAMCHesqspMrAuIgFnyFYAbIlsxJw/3PMj4Y6Z+qzWjEA+DpYCQUDioaCZFbgWwQrQJYsChkIJXtWwqZnhcwKAPhWOlgJyoxbiYWD9KzAtwhWAI8FUotCBpI9KyazQrACAP6VHLziKKBAQLZvhcwK/IpgBciSdVZSDfb5kWBnGRjBCgD4TyKeLgMzTBkYPSvwq1ETrFx55ZX61Kc+peLiYk2YMEHHHnusli5d2uU1LS0tOuecczR27FgVFRXp+OOPV3V1tWfHDAxkdLEyRhe3pEYX07MCAL7uWXEzK6GMzArnBfjLqAlWnn76aRuIvPDCC1q0aJHa29t1+OGHq7GxMf2ab3/723rggQd0//3329dXVlbquOOO8/S4gf5OA+vMrJhgJXUFjZMSAPg1WInbnhW3wb6zZ6XV22MDRlhYo8QjjzzS5fFtt91mMyyvvPKKDjzwQNXW1uoPf/iD7r77bh1yyCH2Nbfeeqt22mknG+Dsu+++Hh050M9pYKGMRSGd1KKQBCsA4Od1VtS9Z4XzAnxm1GRWujPBiTFmzBi7NUGLybYcdthh6dfMnj1b06ZN0+LFi3v9Pq2traqrq+tyA0ZSMJlZCSQb7O3oYk5KAOBfmWVg3XtWyKzAZ0ZlsJJIJHThhRfq05/+tObMmWP3VVVVKRqNqqysrMtrKyoq7HN99cKUlpamb1OnTh324wd6Gl0cDIU2H11MbTIA+DpYMbpOA+O8AH8ZlcGK6V156623dO+9927197r00kttliZ1W7ly5ZAcIzDg0cXJBvv8qKlNJrMCAL6Vuc5KoPuikGRW4C+jpmcl5dxzz9WDDz6oZ555RlOmTEnvnzhxotra2lRTU9Mlu2KmgZnnehOLxewNyKYG+03pnhVGVAKAb9dZcZJlYJmLQnIRCz4zajIrjuPYQOWvf/2rnnjiCc2cObPL8/PmzVMkEtHjjz+e3mdGG69YsULz58/34IiB/gmmVrDvMro4dVJq8vLQAADZMLo4HFJrav0tMivwmfBoKv0yk77+/ve/27VWUn0ops8kPz/fbhcuXKiLLrrINt2XlJTovPPOs4EKk8AwGsrAnNQK9mYaWHqdFTIrAOA7ycWC3Qb7QNfMCj0r8JlRE6zcdNNNdnvwwQd32W/GE59++un2/i9+8Qtb928WgzRTvhYsWKDf/OY3nhwvMPBpYMmeFZNZYQV7APCvLj0rJrNCzwr8KzyaysC2JC8vTzfeeKO9AaNC5vs6FaxEM6aBEawAgK+DFcNkVmpTo4s5L8BnRk3PCpCTEm6qPzNYyWNEJQD4W7d1VmzPSvq8QHkw/IVgBciCE5KV0WDPCvYA4GOZ66wEuk0DI1iBzxCsANkSrKhzdHHnNDBOSgDg13ODY3pWFOjas8J5AT5DsAJkY2YlNQ2M0cUA4NsS4c7RxUG1pHpWyKzAZwhWgCwYT2mZM1I6s+JeQXM4KQGA/8Tb7KYtOQfJBiv0rMCnCFaAbMmsBMKdi0KmRxc3dZ0YBgDIfcnxxG2KpFewp8EefkWwAmRNsBLY7ApawDwfb/fq6AAAXkj+3W9zwgoETM8KvYzwL4IVwEuZWZNkz0owGFAinNe5n74VAPCXeCqzYoKVZGYllXEnswKfIVgBsmSdldQK9kY4ElPccTMtnJgAwJ89K+0mWDFlYKGM0cWm15GMO3yEYAXwUqLDbkxgkhms5EfDnc2UZFYAwF86Ug32EXcaWOZiwfZ5LmLBPwhWgCyYBtahUKplpYfxxZyUAMCvmRUlMyutqcyKfYLzAvyDYAXIhsxKckHIFDsRLJ1ZYRV7APCVtga7aXLMRauAYhFzjgh0BiwdnBfgHwQrQBb0rHSYVYozUit5EbMAWKqZkpMSAPjqvPDMNfZunQrcBvuQ+3Gtc3yx24AP+AHBCpAVqxQHbRNlilkYspnMCgD4T93q9N0qZ4w9N5ieFSOdWeG8AB8hWAGyoAzM9Kxk6tqzwkkJAHyjaUP67v3xg2zWPZVZ6cy4k1mBfxCsAFnRs2LKwLpmVjpXsSdYAQC/ZdwrAxPSZV+pzEqzQ88K/IdgBciCaWA2WMkoBDMnps7MCqOLAcBviwU7yXOCLQMLudn3zgZ7MivwD4IVIAuuoMWdrqOL6VkBAL9yumxS66wYTImEHxGsAFnRs9K1wd72rNiRlWRWAMDfmZXOnpXWdBkYmRX4B8EKkC3TwLr3rHAFDQB8qFuwEpCCwYDCwUDneYGeFfgIwQqQNdPAuq6zwjQwAPBzZqWrWDhjFXsyK/ARghUgKxrsu/asuGVgqcwKZWAA4NfMSko0HCTjDl8iWAG8FG+3m3YTrGTsZp0VAPCpbpmV1IWsWDiU0bPS4tHBASOPYAXwUjJrYq6WmYW/eu5ZIbMCAP6xeYN9KrOSWneFYAV+QrACeCmZNTGTvzIzK/nRzGlgZFYAwLfTwNKZlYxgpZ1gBf5BsAJkTWZF3RrsqU0GAP/puQzM7VmhDAz+Q7ACeCl5dcwEJps12Kd7Vho9OjgAgFeZlcwJkanMSktq8ArBCnyEYAXwUmt9RhlY5uhi1lkBAH/qq2eFzAr8h2AF8PLq2dNX2bvrVNrlIpppsGcFewDwoT6mgXVexCJYgX8QrAAeZ1WMtxIzexhdTGYFAPyqezGYzawwuhg+RLACeKVpffruosS8zUYXs84KAPhRz9PAuiwKSbACHyFYAbyScFevbwgUZlQlZ0wDSzZSOqYMLN1wCQDIaZv9vQ9kjC4mswL/IVgBsmyWfmeDvZtZCTgJKd7mzTECAEZY13NDl2lg9KzAhwhWAK+YIKSHiS/pk1IgeVIyaLIHAH9IJlY2XxQylNGzQnkw/INgBfBMzxNf3PsBhcMxtTshdwd9KwDgz3ODeupZafXm0AAPEKwA2VIG1u1pVrEHAB+fG5xUZiWzZ4VzAvyHYAXIlrrkbtFKfkbfCmVgAOAXPV/IioZMZiVVBkZmBf5BsAJkSc9Kd3atleREMK6iAYA/F4VMiUWCakmdE0zPClMi4RMEK0DWlIF1DVrchSHJrACAv/SyzkooY3SxudiV6PDsCIGRRLACZGGDvZEfNWVgZFYAwM+ZldSFrFgk1NmzYnBegE8QrACen5CCvTfYO2RWAMBf+pFZMehbgU8QrACe96yoy8SXzAZ7poEBgM/00rNiRheby1pt6fHFnBfgDwQrgGd6X8E+lfInWAEAv+n53GBGFxttqQWDyazAJwhWgGxZpbjb04wuBgAf6jZ8pWtmRay1At8hWAGyNLPi9qxwUgIAfw9fSS0KGbLbtvRaKy3eHB4wwnIyWLnxxhs1Y8YM5eXlaZ999tGSJUu8PiSgH+us9NSzQmYFAHy9gn23zEp6SiTBCnwi54KV++67TxdddJF+8IMf6NVXX9Xuu++uBQsWaO3atV4fGtD3eMrNMiv0rACA//R8bkj1rKQngrUTrMAfci5Yue666/SNb3xDZ5xxhnbeeWf99re/VUFBgW655RavDw3oxj0VJXodXRxSC6OLAcBfelkwOBWsNDuUgcFfcipYaWtr0yuvvKLDDjssvS8YDNrHixcv7vFrWltbVVdX1+UGeJNZ6WkFezIrAOAvPTfYp3pWWghW4DM5FaysX79e8XhcFRUVXfabx1VVVT1+zZVXXqnS0tL0berUqSN0tPC9bj0rPU0D6+xZIVgBAD+XCKd6VtKDVwhW4BM5FawMxqWXXqra2tr0beXKlV4fEnwjfSrqxzQwysAAwE+6X8hKlYE1pTIr9KzAJ8LKIePGjVMoFFJ1dXWX/ebxxIkTe/yaWCxmb4BXV88S6YkvgT7WWSGzAgC+OjekHnfLrKQb7MmswCdyKrMSjUY1b948Pf744+l9iUTCPp4/f76nxwZsLpnqD/RnGhiZFQDwh+6jiwNdgxXKwOAzOZVZMczY4tNOO0177bWX9t57b/3yl79UY2OjnQ4GZGXPSipaUQ/BSnoaGJkVAPDlNLDkKSIcDCgYILMC/8m5YOWkk07SunXrdPnll9um+rlz5+qRRx7ZrOkeyNYTUkp+lGlgAODffsauzMRIk11pSfcyEqzAH3IuWDHOPfdcewOyW6ouueeeFdtgzwr2AODzdVa6ji9ubU9lVriIBX/IqZ4VIKcyK3ZRSDIrAOAv3c8NnSeHLpmVjlZvDg8YYQQrgMfBSkrPDfaxztrkRHwEDw4AkBXrrGQ8ZcYXp3tWuIgFnyBYATzTPdW/+Qr2jcrr3EEpGAD4QO9Zd5tZSfUyklmBTxCsAJ7P0u9tGph7BS2emhbW1jiSRwcAyDK2ZyW1KCQ9K/AJghXAM333rERDQQUDATWlsisEKwDgwwb7bj0rZFbgMwQrgNfrrCQfds+vmKZKUwrWlOpbIVgBAB/o2rOSeXKIZQYr9KzAJwhWgCydBpaaCNbokFkBAN/o49zQpcGezAp8gmAF8Ez3npXNoxU3s0KwAgD+0TXfHuhWHtw5upjMCvyBYAXI4syKabJPTwRraxjRwwMAZMHo4oyTQyw5eMUiswKfIFgBvO5ZSU776mkmmM2sOPSsAIB/dL2Qpe6ZFXpW4DMEK0CWlIFlXj3r0rNCGRgA+HgaWLfRxWRW4DMEK0AWrlKckh81mRXKwADAP7qXganr6OLMnpXkeQTIZQQrgGe23LNSGA2TWQEAHy8YnLnOSpdpYKaUON7uzTECI4hgBfB8nZXNT0gpBSazwjorAOAffQxfMZmV1lTPitHRMuKHB2R1sPKlL31JjzzyiBzSjsDQXz3rIbNSEMtcZ4UyMADIfb0XB3fpWTEIVuADAwpWNm3apKOOOkrTpk3T5Zdfro8//nj4jgzwib5i/8JomHVWAMDHF7K6j7M3QUx7INW3QrCC3DegYOXxxx+3AcrChQt15513avvtt9chhxyiu+++W62tTKUAhnqdlQJ6VgDA3yXCga7j7I22QLI8uJ1gBblvwD0r06dP1w9/+EMbtCxatEiTJ0/WN77xDU2aNEnnnHOOXnnlleE5UiDnT0ibRyuFscxpYAQrAJD7uk+KDHQZZ2+QWYGfbFWDvcmqmAxLVVWVrrzySt17773aZ599hu7ogJzWfeJLb5mVVIM9PSsA4Oese17UDVbSTfYEK/CB8NZ+g2XLlum2226zt9raWh122GFDc2SA39ZZ6Wl0cWZmpb1pBA8OAJBtF7JSmZXOhSEJVpD7BpVZaWlpsRkVk1kxfSt//OMfbR+LCVzMtDAA/dH16llPWMEeAPxZItxTvj0VrLSkMiv0rMAHBpRZWbJkiW655Rbdd999NmD54he/aIOTQw89tMd6ewADSPX3cGIqjGVOA6MMDAD8VwaW0bMSda8xtzpkVuAfAwpW9t13X+2+++664oordMopp6i8vHz4jgzwydWzhNPXNLDMdVYa3ZMYFwYAIPfLwFLnhm7rrBhNBCvwkQEFKy+//LL23HNPe99kVkymZe3atUokUilL1+c///mhPUrAFxNf1EtmJdlgn+iQ4m1SOPkYAOCr0cX5yQb7Fif58a292YMDBLI4WEkFKo8++qi++tWvasOGDZu9xqQr4/H40B0h4JeFv3rJrKTLwFLZFYIVAPDR8JXNRxc3JqJu13EHa9wh9w2qwf7cc8/ViSeeqDVr1tisSuaNQAXor645lR57VqJhJRRUs5NspqRvBQB8O3wlFaw0pzIrHWRWkPsGFaxUV1froosuUkVFxdAfEeC3npW+VrCPJeuT02utMBEMAPxwHavHYCVdBsYK9vCPQQUrX/rSl/TUU08N/dEAvp4GtrloKKhQMMAq9gDgw56V7hexYuFg1wtY7ZwTkPsGtSjkDTfcoBNOOEHPPvusdt11V0UiyakUSeeff/5QHR+Qw3qvS1bGPjsRLMH4YgDw27kh0MM5IS+SWRrMYsHIfYMKVu655x499thjysvLsxmWzA9Z5j7BCjA0mRWjMBpWUwtlYADgr+ErwR4vYpm+lebW5AUspoHBBwYVrFx22WX60Y9+pEsuuUTB4KAqyQB0mwbW2/Ippm+loTnffdBaP2KHBwDwsgys54tYJlhpaqUMDP4xqEijra1NJ510EoEKMMzTwIzCaFgNIlgBAL9PAzPyoiE1izIw+Megoo3TTjtN991339AfDeDjWfq91YGZnpUGh2AFAPxWItxTxt2WgaWngRGsIPcNqgzMrKVy9dVX28Uhd9ttt80a7K+77rqhOj4gh3XWJffFrGJfrwL3QWvdSBwYAMAzGcFKL2utdE4DI1hB7htUsPLmm29qjz32sPffeuutLs/11AwGoO+6ZKO3/3TMXH3KwADAf6OLe8q455nMSnrtLYIV5L5BBStPPvnk0B8J4Df9ngYWUj1lYADguxLhns4LJljZkC4Do8EeuY8OecAz3YKVXlIrBdGMMrAWysAAwM8lwibb3lkGxuhi5D6CFSDbMytmdDGZFQDw3+jiHhvsg5SBwVcIVoBsqEvua52VLqOLyawAgG+mgfXSYN9lGljy9UCuIlgBPJM8ISXPM72ts1KcF6ZnBQB8p+fRxXmZZWDmPNLRMtIHBowoghUgW8rAesms2GCF0cUA4Kuse6KXS1g2s5IOVigFQ+4jWAE8b6Lse9x3USxCzwoA+HEaWA9Xscw0MNN83x5IrmLPRDDkOIIVIMt7VkxmJd2zEm+TOlpH7BABACOta9a9p8yK0RrIc3eQWUGOI1gBvJLsVemcBtbziakoZoKV5EnJYHwxAOR8ZqWvMjCjNcAq9vAHghVgFGRWHAUzSsEIVgDAD+eG3hrsDYIV+MWoCFY++eQTLVy4UDNnzlR+fr623XZb/eAHP1BbW1uX173xxhs64IADlJeXp6lTp+rqq6/27JiBgfas9Na5UpwXsdv69Phi+lYAIHelRhGbYKXn0cVGcyrjThkYclxYo8B7772nRCKh3/3ud9puu+301ltv6Rvf+IYaGxv185//3L6mrq5Ohx9+uA477DD99re/1ZtvvqkzzzxTZWVlOuuss7z+FYA+J770tYK9KQMzbGbFvIRgBQByvwzM6TmzkgpWmjLXWgFy2KgIVo444gh7S5k1a5aWLl2qm266KR2s3HXXXTbTcssttygajWqXXXbRa6+9puuuu45gBVkerPSd4IyGg4qFgxmZFcrAAMAXZWA9PJ2fLANrclLTwAhWkNtGRRlYT2prazVmzJj048WLF+vAAw+0gUrKggULbFCzadOmXr9Pa2urzcpk3oCRPCGl9DXA2E4EY3wxAPhA36OLC2NusNKYClbaGF2M3DYqg5UPP/xQ119/vb75zW+m91VVVamioqLL61KPzXO9ufLKK1VaWpq+mV4XYGQnvgT7bLBPlYLRswIA/lowONjDeaEw6hbF1CdSmZXmkTw6YMR5Gqxccskl9qpBXzfTr5Jp9erVtiTshBNOsH0rW+vSSy+1WZrUbeXKlVv9PYGh7FlJNdnXO8lV7FtqR+b4AAAeD18J9FoG1hCnDAz+4GnPyne+8x2dfvrpfb7G9KekVFZW6jOf+Yz2228/3XzzzV1eN3HiRFVXV3fZl3psnutNLBazN8Dr0cV9cddaIbMCAP4pEe65wb4wmVlpUvKzC2VgyHGeBivjx4+3t/4wGRUTqMybN0+33nqrgsGuSaH58+frsssuU3t7uyIRd9TrokWLtOOOO6q8vHxYjh8YqsxKXyVgRlHmKvY02AOALxYM7unUkBcJ2nNGOlihDAw5blT0rJhA5eCDD9a0adPs9K9169bZPpTMXpSvfOUrtrnerMfy9ttv67777tOvfvUrXXTRRZ4eO7AlvZ2QujfY1zmF7gPKwAAgh2WsYN/DycGUDBdEQmpmdDF8YlSMLjYZEtNUb25Tpkzp8pyTbEQzzfGPPfaYzjnnHJt9GTdunC6//HLGFmOUZFb6DleKY2HVpoKV5pqRODoAgKclwlKwl3NDQSys5ibKwOAPoyJYMX0tW+ptMXbbbTc9++yzI3JMwHDP0u/eYL9aqcwKwQoA+GEaWG/nhoJoSE2pYIXMCnLcqCgDA3J7dHH/elZq0pmV3tcNAgCMdhnBSm+ZlWg4owyMnhXkNoIVwPPMSnCLuRUzDaxWRe4DysAAwBeZFfWRWWkWi0LCHwhWgCyoS95SHZhpsK/NbLBPnswAADlcIhzoowxMee4DysCQ4whWAK8b7B2TWelHsJLqWXHirLUCADnL2WKDfWGXMjCCFeQ2ghXAMwPoWYlF1Kqo2uSuIUTfCgDkej+ju57KlsvACFaQ2whWgFGwgn1Jvju4L923wkQwAMj5EuFep4HFTBkYmRX4A8EKkBWji/sOWMry3StomxIF7g6a7AHAt9PACm0ZWEbPCn2MyGEEK4BXkueW/pSBlea75V/pvhUyKwDg23VW8jPLwAzGFyOHEawA2bCC/RZemhcJKhoOsoo9APgqs9LzKwpNZiVVBmZQCoYcRrACZMV4yr7DFfO8ya6kMys02ANA7mdWejk3mMyKacBvC7DWCnIfwQowCjIrhg1WHBrsAcAfwUrvDfaFsZDdtgZSfSuUgSF3EawAnkmdkIIKBrccrpTZYIUyMADIbf0ZXexOiGxJTwQjs4LcRbACZMF4yn7EKl3LwMisAIAPzg2BXtdZMdJ9K6y1ghxGsAJ4XgYW7PWElKm0gMwKAPipZ6U3qcxK51orlIEhdxGsAJ6vUrzlBnuDBnsA8IPONVMCW8isNCWSDfaUgSGHEawAWTANrL9lYDWpBnuCFQDI7QtZTrD3BvtkZqXRoQwMuY9gBciK0cXqV4N9jQhWAMAvPSu9Ntgnp4HVpzMrBCvIXQQrQBaUgfW3Z2WDU+I+aK2TOlqH+wgBACOus2dlyw32BCvIfQQrQBaMp+xXsJIfUZ0KFE/9Z9u0YbgPEADgaT9jzy/JC4fsc82UgcEHCFaALE71ZyrNj9o1WWqUzK4QrABAbpcI9/ISszZXfiSkJqUWhaTBHrmLYAXIigb7/mVWjI0qdnc0rh/e4wMAeDgNLNDnlSwzvjg9upjMCnIYwQqQBRNf+jMNrKzADVbWx5PBCpkVAMjhdVZsuNKrwlhIjU4ys9JGZgW5i2AFyOJVinvKrGxIZVYIVgAgZ7n9jL0/36UMrK1hxI4LGGkEK0AWrGDfn56VSCiowmhImxzKwADAH/2MvZ8cimJhMivwBYIVwCvpYKV/PStGeWFUG9MN9gQrAJC7ZWC9N9gbhbHMnhUyK8hdBCtAFs/S725MYVQbUpkVysAAIAdteXRxKrPSoHz3AZkV5DCCFSAbrp71L1axwUpnGRjBCgDkatbdtNf3lVsxDfZN6TIwMivIXQQrwCgqA7OZFcrAAMAf08ACfZeBNaYb7MmsIHcRrACer1IcVLCf/yWOKYhqo8OikACQu/qXdTdlYE3pFewJVpC7CFaAUTK62BhTZIKVVM/KRimRKhcAAOTWhawtlYFlZFY6WqR4x0gdITCiCFaALFjBvq/xlJnGmp6V1DorTlxqqRnOIwQAeHhuCPWx0Io7DSwZrBjtZFeQmwhWgCxYZ6U/K9gb5QVRtSushkCRu6Nx3TAeIAAgW8vAimNhtSmsDoXcHZSCIUcRrACecQZcBja2KGq361Xq7mioHr7DAwB4OilyS5kVMzGsJcD4YuQ2ghVgFGVWxhS6zZTViTJ3Rz3BCgDkbBlYoO/RxUa6FIzxxchRBCvAKOpZMdPAjKpEKrNSNXzHBwDwQObo4kCf08CMzlXsyawgNxGsAFkw8aW/mZWS/LDCwYDWOsnMCmVgAJBb3FNDsgxMWygDkxrSC0MSrCA3EawAWVCX3N+eFXOVrbww2hmsUAYGADm7YHBfPSupzEpdIpVZoQwMuYlgBcjyE1JPpWCdmRXKwAAgd6eBbanBXmoks4IcR7ACZEGDvSnt6q8xJrOicvcBmRUAyN1pYH0EKwWR7g32BCvITQQrgNeFyZJCwf7/p2hWsV/n0GAPALmpf6OLg8GACqMhNTqUgSG3EawAXmdWnGCfTZTdjcvsWWmpldpbhukAAQCeTYp0+l4UcrNV7FsJVpCbCFaALOhZCQ8gszK+OKY6Fao94I4xZiIYAORiGZj6LANLNdk3UgaGHEewAmTB6OKBNNhPKDYnpoBqgmPcHQQrAJBD+n9usJkVGuyR4whWgCxYFHIgDfbjS9z65HVKjS+mbwUA/LhgsFnFvjOzQhkYchPBCpAFJ6SBZFbGF8W6rmJPsAIAuTkNbAuf0opiETWlG+zJrCA3jbpgpbW1VXPnzrVXG1577bUuz73xxhs64IADlJeXp6lTp+rqq6/27DiBgZSBhUMDKANLZlaWtSfLwGpXDs/xAQC8nRS5xZ6VzMwKwQpy06gLVr773e9q8uTJm+2vq6vT4YcfrunTp+uVV17RNddcox/+8Ie6+eabPTlOYMsGvoK9MbYwJpOIqXTGujtqVw3XAQIAPLuQFexHGVjGNDCCFeSoURWsPPzww3rsscf085//fLPn7rrrLrW1temWW27RLrvsoi9/+cs6//zzdd1113lyrMDApoH1P1gxJWNji2Ja7YxzdxCsAEBuTgML9mMaWLrBnp4V5KZRE6xUV1frG9/4hu644w4VFBRs9vzixYt14IEHKhpNjnOVtGDBAi1dulSbNm3qs6zMZGUyb8DI9qyYdVYG9p/ihOIYmRUAyPmsu7aYWaEMDLluVAQrjuPo9NNP17e+9S3ttddePb6mqqpKFRUVXfalHpvnenPllVeqtLQ0fTO9LsBIXz0bSM9Kaq2VylRmpX6NFG8fjiMEAHjYzxjZQoe9G6zkuw9audiK3ORpsHLJJZfYesy+bu+9956uv/561dfX69JLLx3yYzDfs7a2Nn1buZJmZYx0GZjJrAQGnFlZrxJ12IUhHamucpgOEgDg2Vj7LQQrpsG+zklWm3S0SB2tI3GEwIgKy0Pf+c53bMakL7NmzdITTzxhy7xiseR4viSTZTnllFN0++23a+LEibZULFPqsXmuN+Z7dv++QDb3rKQyK3ZhyGiFxrWudEvByqcP04ECAEZ+GlhA0S1k3U1mpV759jwSNF/XUicVjR+RowR8EayMHz/e3rbk17/+tf77v/87/biystL2o9x3333aZ5997L758+frsssuU3t7uyKRiN23aNEi7bjjjiovLx/G3wLY+mBlINPAOlexl9YHx2ucTLBCRhAAcq1EuD9lYKbvsSlQoCKnUWqpJVhBzhkVPSvTpk3TnDlz0rcddtjB7t922201ZcoUe/8rX/mKba5fuHCh3n77bRvI/OpXv9JFF13k8dEDW26iDA+iDMxYnWqyryFYAYDc0Dm6eMtlYO4153olS8FMsALkmFERrPSHaY43Y42XLVumefPm2RKzyy+/XGeddZbXhwZseQX7ATbYV5S6mZWPO5JN9ps+GfrjAwB4eG7QlsvAom6wUucUujtaaob98ABflYEN1owZM+yEsO522203Pfvss54cE7A1E1+iW7h61t02Ze70lzdbxrv/FW/8aFgOEQDgVRnYlqeBpTIrdU6+aXEhs4KclDOZFWA0Z1Zi4YH9pziuKGZLxz6KJ8d1b/hwOI4QADDiOi9kbakMrDAWstuaBGVgyF0EK8AomKXfnRl1PLE0T584yUl3jes4SQFAjl3IimyhDKwoL5lZoWcFOYxgBciCFeyjA8ysGJNL8+1iYC2xZN/KBkrBAGDUyyhz39KFrFg4ZAOa9ForBCvIQQQrgNeji52BZ1aMyWVuk/3GvGnuDoIVAMgB/e9ZSY0vrlOqwZ5gBbmHYAXIghPSoDIrySb7NeFt3B30rQDA6JdMrLglwlueFFkYDZNZQU4jWAGyYFHIwQQrk5LBysfOZHfHuveG9vgAAB73rGz53GAmgtGzglxGsAJkQ7AyiDKwbZJlYK+1uwujqvqtoT0+AEDWl4GV5JvMCuusIHcRrAAeBysabGal1M2sLG6Y1Nmz0towlEcIAPB0UuSWy8BK8yOqcYrcB82bhvvogBFHsAJ4PO1lMKOLjW3Kk2VgzQVKFJn1Vhxp7TtDepgAAO8uZPUrs5IX0XqVdI6xB3IMwQqQBcHKYMrAzAlqTGHU3m8s39ndWfXG0B0jAMC7MjBny6OLjZL8iDY4JZ09Kx1tw32AwIgiWAE8vXI2+GlgxoyxblNldf727o41BCsAkAsXs8y5IdzPMrBaFSoudzV7Na0f7iMERhTBCuCJrc+sGDPGuk2VH4S3c3esenloDg8A4BFnQOcGk1kxiws3hErdHZSCIccQrACeZ1YGt4K9MWOcG6wsie/g7jA9KzRYAoBvRhebzIpRGyxzdxCsIMcQrAAeByuDXWfFmJ4sA3urNiaNNdkVR1rx4pAdJgBgZDmDKAMzNqab7CkDQ24hWAGyoGelP+MpezIzmVn5ZEOTNG2+u3PF80NzjACArF9nJRWsrE812ZNZQY4hWAGyYRrYoDMrbrCyrr5VLdvs6+5c9szQHCMAYMQFupSBBfq1KKRRHS92dzRUD+8BAiOMYAXIhsxKcHD/KZoramOT44s/KtnHzuVX5b+lujVDdqgAgJFnLmkNJLOyvL3c3VFXOdyHBowoghXA42AlHAoqGBxcGZgxe5J7Ne3t2jxpyl7uzvcf2fpjBAB4lnW3PSv9ODekgpWVibHujpqVw3d8gAcIVgBPdJ6QwqHkbPxB2rHCrVN+t6pO2vFId+d7/9y6wwMAeBqsmOb6QGDLwUp+JGSDmkonGazUrhrOIwRGHMEK4PEJKbSVwUoqs7K0ql7a6fPuzo+ekOqrtu4YAQAjy4mn7wZDbi/KlpiAxmRXVjvj3R31a6R4+3AdITDiCFYAj8vAIlsZrOw0MZlZWVMnx4wvnrK3e8J7/d6tPkwAwAhKZAQrwf4FK0Z5YVQbVKxE0PQwOvStIKcQrABeSHS4GyegcLj/J6SebF9RJFPWvKmpXWvrW6U9T3WfeOVWKe7+HADA6LqQFepHc33KuKKoXWC4KX+iu4NSMOQQghXAw6tnHQoqNsixxSl5kZC2HV9k77+5qlaac7xUMFba9In01p+H5HABACNbBhbqZxmYMa4oZre1sUnuDvP3H8gRBCuAl5kVBQe9xkqmPaaV2e2rKzZJ0UJp/rnuE0//TGpv2ervDwAY2TKwwABKhFPBSlVkqrtj/dKhPzbAIwQrgIdXzzoU6tcc/S3Zc1p5Z7Bi7P0NqahC2vix9MzVW/39AQAjXQbW/8zK+GI3WFkRmOLuWPf+0B8b4BGCFcDDq2dDlVmZN90NVl5fWauOeEKKFUuf+7n75L9+Ib33kIbFq3+U7jpBWvbM8Hx/APDtGlwDKQNzFwf+IDHZ3UFmBTmEYAXwsAzM9KxEhyCzYnpWSvLCam6P69019e7OnT8v7Xmae/K7/zTp5VulROeJcKut/0D6x/nSB49J950qNSezOgCArbuQZYevDKTB3s2svNWe0bNCCTByBMEK4OEJKW7KwIYgsxIMBvSpGWPs/X99uL7ziaOudddeibdJD14o/WZf6emrpXVDcNXNBCmpxS1baqRX79j67wkAfpYsEY4rOKAS4VSw8kFjgZRX6l6k2vDBsB0mMJIIVgAPMyvxIcqsGAft6C4I9vT7azt3hiLSCbdLn/2xFC12SwOe/Il0497SDXtLr909+GzLx0+52zHbuts37tvq3wEAfC2jRHhAwUqyZ2VDY7ucibu6Oyv/3fcX1ayQNi7bioMFRgbBCuBxg300HBiSb3nQDm6w8vInm9TQmrG+SjAoffoC6dtvSl+4Udrus1Iw4gYufztbuuckqblmYD+so0365Dn3/tG/kMziZdVvSes/HJLfBQD83LOSUECRUP/PDeOLYgoEpLZ4Qs0T9nB3rn6l9y+ofke64VPura/XAVmAYAXwuC55qDIr08cWasbYAnUkHP3rg4xSsJT8cmmPr0pf/bP03Y+kw34khfPccq7bjx5Yz8nql6X2Rnc9lxkHSDMPcve/89ch+V0AwJcGWQZmBrVUFOfZ+9XFu7g7V/URhLxym9TRIiXapWev28qDBoYXwQrgaYP90IwuTjlkdoXdPvTmmr5faGqa979QOvNRqXC8VPWmdMdxUmuyOb+/JWAmSDGZm12OdR+/84+tOn4A8LVkWa4pAwubv60DsE15vt1+FN3J3bH2Ham1oecXr1rSed9csGqpG+wRA8OOYAXwuC45P9r/hb+25Atz3bGVi96pVmNmKVhvJs+VvvYPKX+MVPmqdM/J/Zsg89GT7nbWwe52x6OkQFCqeoMaaADYyjIw2884wBLhKclg5ePWYqlsupul+eRfPZfxVr/t3jclvGYAy7Knh+DggeFBsAJ4nFkpiPZ/lv6W7Dal1JaCmRHGj71T1b8vqthZ+ur/StEi6ZNnpf9dKMX7CHSaNrplYMZ2h7nbwrHS9E+79999YGt/DQDwdRmY27MywMxKmRusrNrULG13qLvzo8c3f6HJuJgAJa9M2uvM5Oue2NojB4YNwQrgcWalcAgzK4FAQMfusY29f9cLK/r/hdvsKZ18jxSKSe89KD1wfu9TwkzJgLn6N34nqdT9WdbOX3C3BCsAsNXnhsGWga22wUryQtIHiyQnOWI+Zc1rnZn1bVNBDcEKshfBCuCF9ia7aVJMBbGhy6wYX9l7mp0i8/LyTXp95QCmfM08UPrSLW4512t3SYu+v/lJzvj3ne421aeSMvvozlrousqt+RUAQH5vsDdN8wMxpbzAbldsbHL/npuLT5uWuT2JmSqTwcqkudKM/d3pkGYRyQ0fDdEvAQwtghXAy2DFialgCDMrxoSSPB2zm9u7cvOzHw/si3c6Wvr8De79xTdID/2n1NHa+fyyZ9xSsUBImntK168tmSRN3ce9/94/t+6XAACf96wM9Nyw3YQiu/1kQ6Paw4XSjke6T7x+b++ZlVhR599tsivIUgQrgBfa3Aktjcob8mDF+PoBs+zM/X++sUZvrBrgGip7nCIdcZV7/6X/kX57gPTyre4J788L3f3zTpfKpm7+tTsd427f+fvW/goA4OOx9gMPViaX5tmy4va4o+UbmqTdv9y5YG97s3vfDFBJNddP2t3dbneIuyVYQZYiWAG80NaYLgMrHMIG+5SdJ5foi3PdfpL//ue7SiR6KOfqy75nS1/5kzslzCwe+eCF0l+/KTWulSrmSJ/9Uc9ft9Pn3e3y59zVkQEA/WfWPpHUqsiAJ0WansVUduXDtfVu30rZNKlpvfTK7Z1/m01zffFkqXymuy/Vt2Iy52ZSGJBlCFYALzRtsJtap2hYMivGRYfvoLxIUEuWbdRdLy4f+DfYYYF0/r+lQ77vTvqavKe0/7elMx+RYsU9f035dHftFVPK8PItW/07AICvJEfHtyiqgsjAzw3bTXD/Nn9Q3SCFIu7fbOPpn0kNazuz3mZamEm/GxN3kwrGuRn/zPVXgCxBsAJ4oXa13VQ6Y4e8wT6z2fKSI2bb+z996D29u2YQi37ll0kH/qd0xkPSWU9Kh/2w90AlZe+z3K25ktfbgmQAgCHNrBjbV7iZlaXVyQV+9zhVmrir1LxR+sPh0mt3u/t3P7nzi8zUsW2TpWAf/t9W/wrAUCNYAbxg5tybRkhn4pCOLu7ua/Nn6IDtx9l1V8687SVV1fZjwcettcMRUvkM9+T4wk3D//MAIMeGr7Q4UeUPokR4l8kldvvGqlp3h8muHPc/UsFYdzJYol3afoE0fb+uX5hal+XDHtZlATxGsAKMtNpV6WDl9cS2KsobnsyKEQwGdMPJe2rW+EKtqW3Rib9brOUb3H6ZYRMKu6VjxnO/cksPAABbZhbdlbRJRYMqA9t9apmt7jLjizc0JCc5TpgtnfWUdMB/Sp/9sXTCbZ0lYCmpzErVG/zNRtYhWAFG2rPX2c2LidlapzKNKYwO648rLYjo9jP21rQxBfYE9vkbntMjb/VzdfvB2uU4d4Z/W730z+8M788CgFxhmuElbXRKBtXPWJIX0Xbj3VKw1zLX2TKN9od+X/r0BVLUXY+li6IJndPBlj402KMHhgXBCjBSzAKLz/1aevkP9uEvO45XcV5YsfDwlYGlTB1ToD+fPV+7TylVbXO7vnXnKzr/nn8PX1mYqYH+/PVSMCy9+w/prb8Mz88BgFySXJjR9DMOpmfFmDu1zG5fXbFpYF+4c3Kh3zf/PKifCwwXghVgJKz/QLr7RHdVeEmLxn1NixO7aM7k0hE7hAnFebr/W/vpmwe6a7D84/VKHXLtU7rpqY/U0u7O9h9Sk3aT9r/Ivf/AhdLGZUP/MwAgly5orXrZ3n3PmaaCQY61/9SMMXb73Ifu1Ml+m3O8u/3kX1LdmkH9bEB+D1b++c9/ap999lF+fr7Ky8t17LHJqwBJK1as0FFHHaWCggJNmDBBF198sTo6Ojw7XvicOfGseEG6/3Tpxn2kDx5TIhDWDXln6RurFtiA4ayDZo3oIUXDQV36uZ30wLn7a89pZWpqi+tnj7ynQ699Wn9/bfXA12PZkoO+K03ZW2qtdf8dOpI11ACArpY+LNWuUJMT0wuJnVQYG1xm5cAdxtvt66tqtKlxAOummNHzdjV7R3o9OTUMyAKjJlj53//9X5166qk644wz9Prrr+u5557TV77ylfTz8XjcBiptbW16/vnndfvtt+u2227T5Zdf7ulxw4fWvic9eaUboNyyQHr7r5IT1zOap8NartLPaw7WmMKYbXz/zI4TPDnEOduU6s/f2k/XnrC7JpbkaXVNsy649zUd+5vn9MLHA7wa1xczieaEW93FJde8Jj38PTeIAwC4zN/E9/7pLrwr6Y/xw9WimMYVxQb17SaW5mn2xGL7bZ/5YN3AvnjeGe52ye9ZIBJZI+A42f/JwWRHZsyYoR/96EdauHBhj695+OGHdfTRR6uyslIVFRV2329/+1t973vf07p16xSN9tzE3Nraam8pdXV1mjp1qmpra1VS4o4ABLa4Gv0nz0kfP+nOqF//fvqpFkX0945P6/b44XrHmaFJpXn66r7Tdco+01RWMLyN9f3V3BbXLc8t02+e/FCNbW452GE7VeiSI2enV0Pu62v/9tpqLVvfaF971K6TVNjTujEfLJLuOsG9YnfYj6T9LxyuXwcARoeWWmnpI24f48oX7a66ir31qeX/ofyCQr12+eGD/tZXPfyefvv0RzpyzkTd9NV5/f9CE6D8clepoUo65tfSvNMGfQzAlpjP3KWlpVv8zD18M1OH0KuvvqrVq1crGAxqjz32UFVVlebOnatrrrlGc+bMsa9ZvHixdt1113SgYixYsEBnn3223n77bft1PbnyyittEAT0W3ONtPplaeVLckxt78oXFTCz65PanJCeSeymf8b31f8l5ilcUKrDd56oi+dMtGuehEPZldA0TZznfGY7nbjXVP3q8fd1z5KV+r93q/Xk0rU6ee+pOvcz29srdd2ZAOXsO1/Re1XJxcckXf3IUl121GwdO3cbBTJHY27/WWnBT6RH/0v6vx9IeaXSXskreADgl9XpV78irXheWv682xsSd7MXiVBML044Qeet+ZzM5dPPbjduq37U53efbIOVx99dq9qmdjsVsl/CUWm/86THLpOe/Ik057gtLwQMDLNRkVm59957dfLJJ2vatGm67rrrbJbl2muv1WOPPab3339fY8aM0VlnnaXly5fr0UcfTX9dU1OTCgsL9dBDD+nII4/s8XuTWUGf4u3SuqVS5b/lrFyijhUvKrzhfQVMhiDDKmecnonvqn8ldtXz2k2zpkzWftuO06e3G6dPzSjPugClLx+ubbBX5UzAYgQD0iGzJ2jBLhPtlJm8SEhPvb9OP3v4PTW0dmhcUVSf23WSnn5/nZZvaEqfKK84do5K87udIB+9TFp8g3v/oO+5t+DwT0MDgBFjPlaZtUqq35Kq33Zva992zyXJ4CSlMjJNf2/bW7e0fkbrVJ5e2PHW0z+lCSWbXyQaiCN++Yy9mGT+Fp+67/T+f6HpLTRlzGYRSRO4HP7fW3UcwKjOrFxyySX62c9+1udr3n33XSUSCXv/sssu0/HHu9Mqbr31Vk2ZMkX333+/vvlNt85zMGKxmL0BNiVf9ZZU9abaK19X++rXFdv0vkLJrInJE6Q+en+SqNCrzvZ6NbG9Xg7urpLJO2ru9HKdMGusrp45RkU9lUKNEqac639O28v2rly36H0tWbZR//fuWnvrbu8ZY3T9V/ZQRUme2joS9krerx7/wE4ae+mTjfrZ8bulmz0tc9KL5EvPXCM9/TO3fO4LN0hjZo7sLwkAW8N8LjGlUhs/dicd2u3H9gO+s3GZAq11PX7ZOqdMLyZ21EuJ2Xo+sYs+aNnGnl3MGPvDZ43Vl+ZNsReHhuIC1wl7TdUVD76jW59bpq/sPU0hc+WpP8Ix6YgrpXu+LC2+UZp9jDTNNN4D3vD0E9V3vvMdnX766X2+ZtasWVqzxh2ht/POO6f3mwDDPGcmgBkTJ07UkiVLunxtdXV1+jkgrXGD7Stx1i1Vy5p31Fq1VOGNH6ioeXX6JZGMwKTOKdA7znT9O7GdXnN20KYxczVz+gzNnVamr0wp0w8rikZV5qS/9p01Vn/65nybafnrv1dp8Ucb9EF1g1o64po1rkgnfWqqvjZ/evp3N5PGzj90e+2//Th9+77XbJbla7cs0cl7T9MlR8x2yxBMadgh/08aM0v6539Ky/8l3bi3tM+3pP2/LRW4IzcBwNNApHGtVLdaqqtM3tz7Tu0qJWorFahfo2Ci5wZ0ExLEnYA+cSbqXWea3ktMs6OIzc1k4c3aWjtMLNbcicX68qQS7TNzjHaaVNL/YKKfzN/oX/3f+/p4XaMee7tKR+46qf9fvOOR0u4nS6/f4zb+n/2cFC0c0uMDcqoMzKSJzCjiG2+8Md1g397ebjMrV1xxhS0BSzXYm8DGvNa4+eab7fjitWvX9jt70t+UFEZBbXDNcnvFq6X6fTVXvqvEuvdVUPeh8tszVvXtxpxI3k1Mt8HJqth26pgwR+OnbK8dJ5Zo9qRim3UYiUUcR7umtg5bJnb74uX2cXlBRBd9dgcbuKQDO7P42T8vkj5+yn0cKZT2/Jq079nuCE0AGCqJuNS0wS3Palzn3pL3nca16qhbq3j9WgUa1ynSvE5BZ8vLHnQ4Qa1yxmu5U9Hl9olTocpAhcaWlWj6mEK7KO/0sQWaPqbABikzxhYOeWDSm2sfW6rrn/hQM8YW6NFvHziw85fpz7xpPzdQ2/VE6bib3QtOwBDp72fuURGsGBdeeKH+/Oc/65ZbbtH06dNtc/0DDzyg9957z665YkYXm6b7yZMn6+qrr7ZN+GbU8de//nX99Kc/7ffPIVgZRZo22pR7x/qPVb/mA7Wu/UiBTcuU37BCJe19j2s0QclHicn60NlGG/Kmq618O0UmzdHEiZNsQGLGPo4d5NhIdDLZmB/84y29X91gH+9YUaxvf3YHHb5zhYLmZG3+/HzwmPT4j936bisgzTpImnuKtNMxbtkYAGT2VJi//82bkrfO+4nGDWpr2Kh4wwYlmjYq0LxRkeb1irZt2qzXsC8mM1KtclU5Y7TGGaMqZ2xyO8auLt+aX6FgySSNLyuyZbBmBLyZ9miGkZjAZHJZviJZkHGvb2nXIdc+rXX1rTr/kO100eE7DuwbmFLd24+x4/f1uZ9Le39juA4VPlSXa8GKyaRceumluuOOO9Tc3GwXh/zlL3+pXXbZJf0a02Bvpn899dRTtrH+tNNO01VXXaVwuP/VbgQrWcK8LRvXS7Ur1bzuE9Wv/UStG5ZLNasUaVitkuaVKki4H4B7U+/ka4UzwV7pWh2aqvriWXLG7qDCbWZr2sTxmjW+0F7hMg3jGD4d8YTuXrLC9r/UNLn9PyYYPPvgbXXEnInulT7z//dHT0jP/7oz05LKtmx/mLTjUdIOh0v5bgMqgFEs3iGZng7TJ2i3denHieZatTXVqL2xRvEm83iTAs01CrbWKNy6SdG2WkUSLYP6sQknoI0q1nqnVOudEm2Q2SZvKlFjZIzieeOUKJqgUPFElRcX2AEiYwuj9uKVCURMUDKhJDaqMuwPvF6p8+75t02K/OG0vXTI7M6pqf1i+lbMJMdgRPry3e7fYmAI5FywMlIIVkZIe7Ot/21dv1y11cvUvG65EjXLFaqrVH5zpUrbqhVV5zjg3lQ7ZW4wogrV5k9RW/F0BcfOVN6EbTW+YhtNtSn4fBXn9XNsI4ZNTVOb/ufZZbrt+U/sFLFUedgX95iio3abpD2mlrnZFtOs+vq90mt329Wc04JhaZu9pBn7SzMPcFdaJusCjFwZlVlTqq2h27ZRidYGtTXVqb25Th3N9Yo318oxt9Y6BVrqFGyrV6itXuGOekU7GhQdZLDRPfNRoyLVOEWbb51CNYdL1RYtVTyvXPH8cXIKxytUNF4lBXk2+BhjA5CoXXjRbM3j0RSADNSlf3lT9yxZoWgoqMuO2kk7Ty7RK8s3aWNjm8107zWjj35B8zHxL2dJb/5JCsXchX5nHzWSh48cRbAySAQrQ3DFrKFarZtWq27tSjVtWKn2mko59WsUbqxSXss6FbevV1Gic22Ovq6CrVWZVjvjtC44XvWxiWotmKRE6VSFx85QyaTtNWn8WE0b41796rKuB7KWmflvAhZz4qyq6/zQYv4//MyOE2xj/6dmjNHU8pgCZtV7s7Lzew9J697t+o1CUWnS7tLkPdzbpLnS+B0ZhQx/j1o3F4LsrUnqaLFbp61J7ebW3KSOlkZ1tDYp3tqoeEu9nNYGOW0NdmsCj2C7uTUp1NGkcEejIvEmRRPNijqdI/6HSpMTU73yVe8UqF4FNhteZ7cFagwUqC1cpPZIsVojZeqIlSmRV26zq+GiscorKlNZQcwurmtGpJsBHmVmmx9RSX4kK0qwsklrR1wX3POaHnm7qsfnj99ziv772Dl23a1e31v3ny6996D72Iw0PvhSmu6xVQhWBolgpRcdbbYJsWVTlWrXV6ppwyq11ayWU2eCkGrlNa9VUft6lSY2KdjPumBzojKBSHVgnGqjFWoqmKxE8TYKlk9VwbgZKq2YpoljSm0dcI+romNUiyccPfP+Ov3l36v11HtrVZ/MtqRMKI5ptyml2nFisR1wsGvBRk2teUXhlc9Jy56V6is3/6bhfGncdtK4HZK37aWx20tl09yFKAloMRKTpOKtbqBgVgO391vV0dasjrYWtbe62462ViXamxVvb1G8vVWJthYlzLajRU576uubFTCBh9l2tChob80KxVvszZRDhRNm22qDibDiw/7rmabyRuXZW5OTp4bk1jxuVp5ag/lqCRWqLVystlCROqJFcmIlCphbXqlCBaWKFJYpWlCqooJ8O+a9KC+skrxI+r4Z45vLWQ4v/+b+cfEnuvOF5Ta7ves2ZTY4+ecblUo4shPJbj51nh0I0PM3aJce+7704k3u49Kp0kHflXb7sruYJDBABCuD5KtgxZxA66pUt6FSjRvW2ECkva5Kifq1CjWvU7Rlg/LbNqioY5NKnC1nQlLanZDNiJgFrurC49QQHa/2ggoliiYqVDpZsTHbqGjcVI0bO16TyvPtSQr+ZtZoMeuymIUlX/5ko95cXav2+OZ/mkyVmKkZn1Ker7lFmzQ3+LFmtr2vioZ3VVrzjkIdjb3/kGixVLqNVDpFKtlGKp4kFY6TCsZ23szj/DGceAfBnErM2SThOPaDT6LLY0eJhCMnEVfC3OId9n7qsd3G41I8roQTl2O2puwoYV6XcPfZx6mvSchJdKQf263J6po1kZJbx3ywircrEe96337gMl9rAgnz+kSHAnZfuwLmvr0l7zvuNmjvd9g1l8yUqJC9tSmcaFfEbJ12e4uqbUQChv5kpZsVVYuialZMLU7qvtnG1B6IqS0YU6vJXoTy1R4qUDxSqHi4UE6kQIoWKRArVCBWrFBekcJ5xYrkFytaUKxYLF+FeREVREMqiIaT25C9oBQLB8lwj0JmTa1z735V6xvaVFYQ0Q0n72lH0PfKZLof/l5nmW7JFOlTZ7pDUYpZKgL9R7Dix2DFfDhoqVVTzVrVb6hUw8YqtdaakYzVdkRj2ExEaV2vgvaNKunYpGL18cGulytqG1SijSpTbXisGqPj1Zo/QYniSQqVTLJBSPG4aSqfMEnjSwpUkhfmxIVBaWmP6/WVNXp3TZ2WVtfbVZjfr6pXY1vvHwSDSmhaoFrbBiq1Q6hKs8NrNCtQqWlOpUqdnhdo6008EFFHuMB+iDPbjvQ2X4lgVIlgWE4wokQgrI5Aams+poYVD4SVUCDjw7rk2DlE5kN7537zh9c+l/xALyWSOxL2FjBbsy/zsb2Z75a5v/Oxu3XSrzU/1f0+7tY8n3p9ahtU3N4326BjHrtf596P28epfaHk8/ZrzOPkc2abepy5z2RZ0/cD/jvVmP+/TfjSqojaFHG3Ttjebw+kblH3/RaI2lvcvr8iigdj6gjlKW5uJmMYzpMTNoFEvoKRfAWiBQpG8xWKFigUK1QoVqBInrkV2YAiLxqyw0PywmYbVCwSUn4kpEgowN9lbGZNbbO+dccren1Vrb0odOmRO+nrB8zs/b3S1iS9fIs7FKXBXdNOgZC0wwJpzvHS9odLeaPsMxRGHMFKLgQriYRtVKzbsEaNm6rUtKlarclZ8GZhw2DzBkVbNyq/faMKO2pU4tQpqi3Phu+eBVmvUtUEylQfLldzdIza8sYpXjBBweIJChdXKK98kgrHTlL5mAqNLaYkC94wf6rMlb+Vm5q0cmOTVm1q1qpNTXYkp9m/obFV6+vb1Ny+eUCTp1ZNDmzQpMAGu90msF7jVasxgTqNCdSrXPXpbciHH6qzhbkgYkOdQCqMSoU8qcep51NhUOdjszWBo/kaEzC6WxNIhmxwmTD3TZBpgsmACTbdmwk+zVbJANRMPHJCYQXsNqJAyDxntu5NYbON2kbjQCTPrvYdiuYpFIkpEs1XOJrn3mL5ikSiioZDNuMQTd1CQbvGBgEDsvEi0WV/fUv/++oq+3j/7cbpimPnaOa4wr7XNHvrz9Krd0grX+jcb/4bmfUZafbnpFkHS+UzRuA3wGhDsDJKg5Xl776stQ//TLPqX1apUyNzjXagGp2YNqlEdaFSNdkAZKw68scpUThewaIKRcsmqqB8korGTlb52AkqL4yN2AJVwEgsSLmhoU2NbR1qaOmwvTCNre59U6fd3BZXezyhtrhjt/Z+h3mcUEdHXJH2etsLkJcwjcVNiiWa3ZvToli8SSF12JIftxQobq6R29KfSMDsN4877FIx9iN2IOBmNwKB5GN3dWuzTbX/mv/07H99ZjBAIJi8BdL3nUBIgfT9oALB1POp/QEFMr429Vr3dRlbewspEEx9bfLrgqHka8w2lN4XDAa7PA70dAuEFAylXu/uM1uFzDZsn7M/N/04bCe+BUxwYK7C2uNOHrs5BgCeMh8Jb3/+E/304ffs30UTXH9576n6j4O3s6Ob+7RuqfTGfdI7/5A2fND1ubLpbtBipjluM08aM4sewq2QSDj2nGVure3u1p7H2uNqb2u1AzU6WprV3t6shBmo0ZbsjetolZO+tWnb3T+t7bbf2bPfg2BllAYrH7/xnGb95XNd9pkJKTWBEtUFy9QULlNrrFztsbFyCsYqWDRO4aIJipVXqKi8QsVjJmpMWamtJQYAABioT9Y36vJ/vG2HoBgmK3jCvCn66r7TbSN+n8zHynXvSe8+4K6fteolt08sk1kzywQtk/eUJuwkjZ8tjd3WZipzdb2xupYO1Ta3q7apTfUN9Wqqr1FzQ63iLXVKJCfzyU7la1CwvUGh9kaFzS3eqGi8KT1MI5xoU1StpvtMeWpTLOBu89SumNoGVB2wZLf/1t7HnSevEKyM0mClpa1Db9x1qYIzD1Tp1NkqHTNB5SXFjGEEAAAjxnw8XPzRBv3i/97XS59sSu/ffUqpDt9log6ZPUE7VhS762P1pbVeWr7YXfB31RJpzRvulLzuTJZ1zMzkBMep7jAUe5vqDkQxQ1CivUwqG24dbeporlN9XY0a6mvU1FCrloZatZr1hZrM+kIm4KizAUfABhuNdvS3CTLMLd9pUmGgRUVqVqFaFA4MvGpmMFrN2A/bD2f64mKKB01fZVQdtucyotZ9L9Dsg06UVwhWRmmwAgAAkFVBy8cbdNcLK/To21XqMINAkswgnd2nlmn3KWXabkKRZowrtD0uZr2bXplJfNVvSatfkcxaWuved0vIWmu3fDBm2ETBGHdqY0G5nVxnhk/ITLEzPWRmoWDz2JSWWqbsNRVMBdTR0a7W1ma1t7aoPTlCvCNdImXWJWp2MxodnWXA+U6zW947DFoC+WoNmal8ZqCLO5HPTOZLRIrkmDVsYsUKxooUyCtWOK9IkViB2w+X3Nr70QIF7O8dS/7+MfffyWyzvNSOYGWQCFYAAAA2t7a+RYveqdbj767V8x+tV0t7zxkCs1aOWStrfHFME4rzNK4oZveZW2o9ndS463DQ3GQXjS6o+1D59csVaahUqG6Vwg2r7f1wU7WCZtS4x1qdiF1TqCWYr9ZggdqCZlKkG2DY4CJabAOLUH5Jctx3qfIKS5VfXKrC4jKFzYS0WJEbZJmbz3v16ghWBodgBQAAoG9mOMnSqnr9e2WN3qms1cfrGvXJhkZV1/VQ4rXVHFs+VW6nNjbYbZkaVBBotf0a+aZnI2D6NtxeDjvYxOZS3I+4qa2ZGmjGiQfCMQXDUQUj7iQ/M8EvYqb45RXaICNSUKJYQYnyi0pVUFSqopIyFReXqqQwX2HK8kf8Mzdd2AAAABgQ00s7Z5tSe8tkpi5W1bbYLIwZLW9vDa1dpjLWJ7dmglU84ajdLNUQd2yJmbmZtafCwYCdVGoyLyYBEQ4WKRQcr/ZgQHWRkDqSC5Ka9YRSi5Sm7pusjSlFM2VqJXYbUUm+e78o6k4lxOhBsAIAAIAhYcq8TP+KuQFDgVwWAAAAgKxEsAIAAAAgKxGsAAAAAMhKBCsAAAAAshLBCgAAAICsRLACAAAAICsRrAAAAADISgQrAAAAALISwQoAAACArESwAgAAACArEawAAAAAyEoEKwAAAACyEsEKAAAAgKwU9voAso3jOHZbV1fn9aEAAAAAOSn1WTv12bs3BCvd1NfX2+3UqVO9PhQAAAAg5z97l5aW9vp8wNlSOOMziURClZWVKi4uViAQ8CzSNMHSypUrVVJS4skxILvwnkB3vCfQE94X6I73BLL1PWFCEBOoTJ48WcFg750pZFa6Mf9YU6ZMUTYwbyD+sCAT7wl0x3sCPeF9ge54TyAb3xN9ZVRSaLAHAAAAkJUIVgAAAABkJYKVLBSLxfSDH/zAbgGD9wS64z2BnvC+QHe8JzDa3xM02AMAAADISmRWAAAAAGQlghUAAAAAWYlgBQAAAEBWIlgBAAAAkJUIVrLMjTfeqBkzZigvL0/77LOPlixZ4vUhYYRceeWV+tSnPqXi4mJNmDBBxx57rJYuXdrlNS0tLTrnnHM0duxYFRUV6fjjj1d1dbVnx4yRddVVVykQCOjCCy9M7+M94U+rV6/WV7/6Vfv/e35+vnbddVe9/PLL6efN7JzLL79ckyZNss8fdthh+uCDDzw9ZgyfeDyu73//+5o5c6b9/3vbbbfVFVdcYd8HKbwnct8zzzyjY445xq4Ib84Vf/vb37o835/3wMaNG3XKKafYxSLLysq0cOFCNTQ0yEsEK1nkvvvu00UXXWTHyb366qvafffdtWDBAq1du9brQ8MIePrpp+2HzhdeeEGLFi1Se3u7Dj/8cDU2NqZf8+1vf1sPPPCA7r//fvv6yspKHXfccZ4eN0bGSy+9pN/97nfabbfduuznPeE/mzZt0qc//WlFIhE9/PDDeuedd3TttdeqvLw8/Zqrr75av/71r/Xb3/5WL774ogoLC+35xAS3yD0/+9nPdNNNN+mGG27Qu+++ax+b98D111+ffg3vidzX2NhoPzuaC9896c97wAQqb7/9tv0c8uCDD9oA6KyzzpKnzOhiZIe9997bOeecc9KP4/G4M3nyZOfKK6/09LjgjbVr15pLYs7TTz9tH9fU1DiRSMS5//77069599137WsWL17s4ZFiuNXX1zvbb7+9s2jRIueggw5yLrjgAruf94Q/fe9733P233//Xp9PJBLOxIkTnWuuuSa9z7xXYrGYc88994zQUWIkHXXUUc6ZZ57ZZd9xxx3nnHLKKfY+7wn/keT89a9/TT/uz3vgnXfesV/30ksvpV/z8MMPO4FAwFm9erXjFTIrWaKtrU2vvPKKTcmlBINB+3jx4sWeHhu8UVtba7djxoyxW/P+MNmWzPfI7NmzNW3aNN4jOc5k3I466qgu/98bvCf86R//+If22msvnXDCCbZkdI899tDvf//79PPLli1TVVVVl/dFaWmpLS3mfZGb9ttvPz3++ON6//337ePXX39d//rXv3TkkUfax7wnsKwf7wGzNaVf5u9Linm9+TxqMjFeCXv2k9HF+vXrbc1pRUVFl/3m8XvvvefZccEbiUTC9iWYUo85c+bYfeaPTDQatX9Iur9HzHPITffee68tCzVlYN3xnvCnjz/+2Jb8mLLh//qv/7LvjfPPP9++F0477bT0//c9nU94X+SmSy65RHV1dfZiRSgUsp8nfvKTn9iSHoP3BKr68R4wW3MBJFM4HLYXTb18nxCsAFl6Jf2tt96yV8bgXytXrtQFF1xga4fN0A0gdTHDXPn86U9/ah+bzIr5e2Hq0E2wAv/505/+pLvuukt33323dtllF7322mv2gpdptOY9gdGOMrAsMW7cOHs1pPsUH/N44sSJnh0XRt65555rm9qefPJJTZkyJb3fvA9MuWBNTU2X1/MeyV2mzMsM2Nhzzz3t1S1zM030pkHS3DdXxHhP+I+Z5LPzzjt32bfTTjtpxYoV9n7q/3vOJ/5x8cUX2+zKl7/8ZTsZ7tRTT7XDN8yUSYP3BCb24z1gtt2HOnV0dNgJYV6+TwhWsoRJ38+bN8/WnGZePTOP58+f7+mxYWSYfjgTqPz1r3/VE088YUdQZjLvDzP9J/M9YkYbmw8ovEdy06GHHqo333zTXiVN3cwVdVPakbrPe8J/THlo97Hmpldh+vTp9r7522E+WGS+L0yJkKk5532Rm5qammxfQSZzAdR8jjB4T2BmP94DZmsufpkLZSnm84h5H5neFs941tqPzdx77712KsNtt91mJzKcddZZTllZmVNVVeX1oWEEnH322U5paanz1FNPOWvWrEnfmpqa0q/51re+5UybNs154oknnJdfftmZP3++vcE/MqeBGbwn/GfJkiVOOBx2fvKTnzgffPCBc9dddzkFBQXOnXfemX7NVVddZc8ff//735033njD+cIXvuDMnDnTaW5u9vTYMTxOO+00Z5tttnEefPBBZ9myZc5f/vIXZ9y4cc53v/vd9Gt4T/hjcuS///1vezMf8a+77jp7f/ny5f1+DxxxxBHOHnvs4bz44ovOv/71LzuJ8uSTT/bwt3IcgpUsc/3119sPHtFo1I4yfuGFF7w+JIwQ84elp9utt96afo35g/If//EfTnl5uf1w8sUvftEGNPBvsMJ7wp8eeOABZ86cOfYC1+zZs52bb765y/NmTOn3v/99p6Kiwr7m0EMPdZYuXerZ8WJ41dXV2b8L5vNDXl6eM2vWLOeyyy5zWltb06/hPZH7nnzyyR4/R5hgtr/vgQ0bNtjgpKioyCkpKXHOOOMMGwR5KWD+x7u8DgAAAAD0jJ4VAAAAAFmJYAUAAABAViJYAQAAAJCVCFYAAAAAZCWCFQAAAABZiWAFAAAAQFYiWAEAAACQlQhWAAAAAGQlghUAAAbo4IMP1oUXXuj1YQBAziNYAQAAAJCVCFYAAAAAZCWCFQDAkJZHnXfeebZEqry8XBUVFfr973+vxsZGnXHGGSouLtZ2222nhx9+2L5+06ZNOuWUUzR+/Hjl5+dr++2316233pr+fitXrtSJJ56osrIyjRkzRl/4whf0ySefdPmZt9xyi3bZZRfFYjFNmjRJ5557bvq5FStW2K8pKipSSUmJ/V7V1dXp53/4wx9q7ty5uuOOOzRjxgyVlpbqy1/+surr69OvMcf+ta99zX4P8/2vvfbaYf5XBACkEKwAAIbU7bffrnHjxmnJkiU2cDn77LN1wgknaL/99tOrr76qww8/XKeeeqqampr0/e9/X++8844NXt59913ddNNN9muN9vZ2LViwwAY4zz77rJ577jkbMBxxxBFqa2uzrzGvP+ecc3TWWWfpzTff1D/+8Q8bDBmJRMIGKhs3btTTTz+tRYsW6eOPP9ZJJ53U5Xg/+ugj/e1vf9ODDz5ob+a1V111Vfr5iy++2O77+9//rscee0xPPfWU/T0AACPAAQBgiBx00EHO/vvvn37c0dHhFBYWOqeeemp635o1axxz+lm8eLFzzDHHOGeccUaP3+uOO+5wdtxxRyeRSKT3tba2Ovn5+c6jjz5qH0+ePNm57LLLevz6xx57zAmFQs6KFSvS+95++237s5csWWIf/+AHP3AKCgqcurq69GsuvvhiZ5999rH36+vrnWg06vzpT39KP79hwwZ7DBdccMGg/o0AAP1HZgUAMKR222239P1QKKSxY8dq1113Te8zpWHG2rVrbdbl3nvvtaVY3/3ud/X888+nX/f666/rww8/tJkVk1ExN1MK1tLSYrMh5usrKyt16KGH9ngcJlMzdepUe0vZeeedbUmZeS7FlH+Zn5FiSr3M9zbMzzFZnH322Sf9vDmGHXfccQj+pQAAWxLe4isAABiASCTS5XEgEOiyzzxOlWkdeeSRWr58uR566CFbpmUCD1PW9fOf/1wNDQ2aN2+e7rrrrs1+hulxCQaDw3a85tgAAN4jswIA8JQJPE477TTdeeed+uUvf6mbb77Z7t9zzz31wQcfaMKECbYPJfNmGuFNNsRkRR5//PEev+9OO+1kG/TNLcX0x9TU1NgMS39su+22Nph58cUX0/vMUID3339/q39vAMCWEawAADxz+eWX28Z1U+719ttv2wZ3E2QYZkqYabY3TfKmwX7ZsmW2uf3888/XqlWr0tO8zHSuX//61zawMY3v119/vX3usMMOs+Vn5vuY/abh30z1Ouigg7TXXnv16/hM6dnChQttk/0TTzyht956S6effvqQZXUAAH3jry0AwDPRaFSXXnqp7XM58MADbY+L6WExCgoK9Mwzz2jatGk67rjjbBBjAgfTs2LGEBsmI2OyMb/5zW/s+OKjjz7aBi2pci4TCJkRyuZ7m+Bl1qxZuu+++wZ0jNdcc40OOOAAHXPMMfZ77L///rY8DQAw/AKmy34Efg4AAAAADAiZFQAAAABZiWAFAAAAQFYiWAEAAACQlQhWAAAAAGQlghUAAAAAWYlgBQAAAEBWIlgBAAAAkJUIVgAAAABkJYIVAAAAAFmJYAUAAABAViJYAQAAAKBs9P8B+VVROflCbGoAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 800x400 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, gs = braintools.visualize.get_figure(1, 1, 4, 8)\n",
    "fig.add_subplot(gs[0, 0])\n",
    "plt.plot(times, vs[:, 0, 0], label=\"V0\")\n",
    "plt.plot(times, vs[:, 0, 1], label=\"V1\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Ecosystem-py",
   "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.11.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}