{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Tagging\n",
    "\n",
    "Tagging is the mechanism by which we can annotate points or regions in recorded data.\n",
    "\n",
    "![data model](resources/tag1.png)\n",
    "\n",
    "For example, at some point a stimulus is switched on and the system's response to the stimulus is recorded.\n",
    "\n",
    "To define this region we need to specify the starting point and the extent (in this example, the duration) of the stimulus-on time.\n",
    "\n",
    "![data model](resources/tag2.png)\n",
    "\n",
    "Translated to NIX we have:\n",
    "\n",
    "1. The **DataArray** that stores the *response* data\n",
    "2. The starting point (time) of stimulus onset\n",
    "3. The stimulus duration\n",
    "\n",
    "## Tag\n",
    "\n",
    "The **Tag** entity stores this information and links it to the recorded data:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import nixio\n",
    "import numpy as np\n",
    "import scipy.signal as signal\n",
    "\n",
    "def butter_lowpass(highcut, fs, order=5):\n",
    "    \"\"\"creates a butterworth lowpass filter.\n",
    "\n",
    "    Args:\n",
    "        highcut (double): the cutoff frequency in Hz\n",
    "        fs (int): the sampling rate of the data\n",
    "        order (int, optional): the order of the low-pass filter. Defaults to 5.\n",
    "\n",
    "    Returns:\n",
    "        b, a (np.array): the filter coefficients\n",
    "    \"\"\"\n",
    "    nyq = 0.5 * fs\n",
    "    high = highcut / nyq\n",
    "    b, a = signal.butter(order, high, btype='low')\n",
    "\n",
    "    return b, a\n",
    "\n",
    "def butter_highpass(lowcut, fs, order=5):\n",
    "    \"\"\"creates a butterworth highpass filter.\n",
    "\n",
    "    Args:\n",
    "        highcut (double): the cutoff frequency in Hz\n",
    "        fs (int): the sampling rate of the data\n",
    "        order (int, optional): the order of the low-pass filter. Defaults to 5.\n",
    "\n",
    "    Returns:\n",
    "        b, a (np.array): the filter coefficients\n",
    "    \"\"\"\n",
    "    nyq = 0.5 * fs\n",
    "    low = lowcut / nyq\n",
    "    b, a = signal.butter(order, low, btype='high')\n",
    "\n",
    "    return b, a\n",
    "\n",
    "def butter_bandpass_filter(data, lowcut, highcut, fs, lporder=1, hporder=1):\n",
    "    \"\"\" Applies an butterworth bandpass filter to the data. \n",
    "    Args:\n",
    "        data: the data that should be filtered. Must be one dimensional.\n",
    "        lowcut (double): the lowpass cutoff frequency in Hz\n",
    "        highcut (double): the highpass filter cutoff frequency in Hz.\n",
    "        fs (int): the sampling frequency of the data.\n",
    "        lporder (int): the order of the lowpass filter. Defaults to 1.\n",
    "        hporder (int): the order of the highpass filter. Defaults to 1.\n",
    "    \n",
    "    Returns:\n",
    "        y (np.array): the filtered data.\n",
    "    \"\"\"\n",
    "    data = np.squeeze(data)\n",
    "    if len(data.shape) > 1:\n",
    "        raise ValueError(\"data must be 1-D!\")\n",
    "    \n",
    "    lpb, lpa = butter_lowpass(highcut, fs, order=lporder)\n",
    "    hpb, hpa = butter_highpass(lowcut, fs, order=hporder)\n",
    "    y = signal.lfilter(lpb, lpa, data)\n",
    "    y = signal.lfilter(hpb, hpa, y)\n",
    "\n",
    "    return y\n",
    "\n",
    "def create_data(duration, dt, stim_on, stim_off, stim_amplitude):\n",
    "    time = np.arange(0., duration, dt)\n",
    "    stimulus = np.zeros(time.shape)\n",
    "    stimulus[(time >= stim_on) & (time < stim_off)] = stim_amplitude\n",
    "    response = butter_bandpass_filter(stimulus, .25, 10., 1. / dt)\n",
    "\n",
    "    return time, stimulus, response\n",
    "\n",
    "\n",
    "def main():\n",
    "    interval = 0.001\n",
    "    duration = 3.5\n",
    "    stim_on = 0.5\n",
    "    stim_off = 2.5\n",
    "    stim_amplitude = 1.0\n",
    "\n",
    "    time, _, response = create_data(duration, interval, stim_on, stim_off, stim_amplitude)\n",
    "\n",
    "    nixfile = nixio.File.open(\"tagging1.nix\", nixio.FileMode.Overwrite)\n",
    "    block = nixfile.create_block(\"demo block\", \"nix.session\")\n",
    "\n",
    "    # *********  This is the interesting part **********\n",
    "    data = block.create_data_array(\"response\", \"nix.sampled\", data=response, label=\"voltage\", unit=\"mV\")\n",
    "    data.append_sampled_dimension(interval, label=\"time\", unit=\"s\")\n",
    "\n",
    "    stim_tag = block.create_tag(\"stimulus\", \"nix.stimulus_segment\", position=[stim_on])\n",
    "    stim_tag.extent = [stim_off - stim_on]\n",
    "    stim_tag.references.append(data)\n",
    "    # **************************************************\n",
    "    nixfile.close()\n",
    "\n",
    "main()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The **Tag** is created with the *position*, and is given the *extent* to the stimulus-on segment. It \"refers\" to the response data. \n",
    "\n",
    "**Notes:**\n",
    "\n",
    "1. **Tags** can tag points or segments in several **DataArrays** at once, hence they store a list of *references*. \n",
    "2. The data might be n-dimensional therefore *position* and *extent* need to be lists (one entry for each dimension).\n",
    "\n",
    "## Reading tagged data\n",
    "\n",
    "Tagging alone is nice, but we may want to work with the data dring the stimulus-on period, that is, we need to read it from file.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAugAAAHwCAYAAAD0N5r7AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAABYlAAAWJQFJUiTwAABeX0lEQVR4nO3dd3Qc1d3G8eeqF8uSZbk3uXfABWxsU0wxLZQQSEILAQKhBEgCIYU0QiBvIKEFSAiBAIFA6L2DjQ02GNu49yJ3W5Zl9b667x+7Wu0KyVZZzcxK3885OquZ2Z29Go1Gj67u/K6x1goAAACAN8S43QAAAAAA9QjoAAAAgIcQ0AEAAAAPIaADAAAAHkJABwAAADyEgA4AAAB4CAEdAAAA8BACOgAAAOAhBHQAAADAQwjoAAAAgIcQ0AEAAAAPIaADAAAAHhLndgOcZIzZIqmrpByXmwIAAICOLVtSkbV2cEtf2KkCuqSuycnJmaNHj850uyEAAADouNasWaPy8vJWvbazBfSc0aNHZy5evNjtdgAAAKADmzRpkpYsWZLTmtcyBh0AAADwEAI6AAAA4CEEdAAAAMBDCOgAAACAhxDQAQAAAA8hoAMAAAAeQkAHAAAAPISADgAAAHgIAR0AAADwEAI6AAAA4CEEdAAAAMBDCOgAAACAhxDQAQAAAA8hoAMAAAAeQkAHAAAAPISAHmWstW43AQAAAO0ozu0GoHm255fpT++s0Zx1+1RrrWYMy9J1M4dpwsBubjcNAAAAEURAjwIbc0v07UcWKL+0KrjuwzW5+mhtrm45ZZSuPm6IjDEuthAAAACRwhAXj6uqqdUNz34VFs7rWCv9+d21+uv7611oGQAAANoDAd3jnvliq1bvLpIkJcTF6PkfHq0Pf3qcjsrODD7nwdkb9ezCbW41EQAAABFEQPewal+t/jVvS3D5ppNH6KjBmRrWs4ueuuIozRzZI7jtN6+u1LLtBS60EgAAAJFEQPewj9bkamdBuSQpMzVB3zs6O7gtKT5WD144UWP7dpUk1dRa/fh/S1VaWeNGUwEAABAhBHQPe3vF7uDnFxw1QMkJsWHbUxPj9PeLJqlLov9e3y15pbrj7TWOthEAAACRRUD3qMoanz5emxtcPmN830afN7B7im4/Z2xw+b9fbNPirfnt3j4AAAC0DwK6R32xOV8lgeEqg7qnaHSftCafe84R/XTS6J7B5V+9vFLVvtp2byMAAAAij4DuUQu31PeCHz+ix0HrnBtjdNvZ45QSGAKzbm+x/jl3c7u3EQAAAJFHQPeo0IB+1ODuh3x+v4xk/fTkEcHlBz/eqL1FFe3SNgAAALQfAroHVVT7tHRHQXD5yMHdmvW670/L1qje/qEw5dU+3f3euvZoHgAAANoRAd2D1uwuUlWNfwx5dvcU9UxLatbr4mJj9OszxgSXX1qyQyt3FrZLGwEAANA+COgetGZ3cfDz8f0zWvTaGcOzdOIo/w2j1kp/eHO1rLWRbB4AAADaEQHdg9bsLgp+XjdkpSV+dcZoxcX4bypduCVf763aE7G2AQAAoH0R0D1o7Z76gD6mT9cWv35ojy665OhBweW73lunGsouAgAARAUCusdYa7U2ZIjLqIPUPz+YG08crrQk/wyjm/eV6sXFOyLSPgAAALQvArrH7CqsUHFggqL05Hj17tq8G0QbykhJ0NXHDQ0u3/fhBlVU+yLSRgAAALQfArrHbM0rDX4+pEfqQScoOpTLpmerR1qiJGlPUYWenJ/T1uYBAACgnRHQPSZnf1nw88HdU9u0r5SEON1w4vDg8sNzNqmwvLpN+wQAAED7IqB7zNb99T3og9oY0CXpu0cO0KDuKZKkwvJqPfLJpjbvEwAAAO2HgO4xOSEBPTsrpc37i4+N0U2zRgaXH/9si/YWVbR5vwAAAGgfBHSP2RoyxCUSPeiS9I3xfYLlGiuqa/XARxsisl8AAABEHgHdQ6y1YQE9u3vbe9AlKSbG6JZT63vRn/tyu7aE3IwKAAAA7yCge0h+aZXKA6UQ0xLjlJGSELF9Hzeih6YOyZQk+Wqt7vlgfcT2DQAAgMghoHvI7sL6seF9MlpX/7wpxhjdcuqo4PIby3Zp1a7CiL4HAAAA2o6A7iF7QgJ67/TkiO9/4sBuOml0r+DyX95bF/H3AAAAQNsQ0D1kd2F58PO+6ZHtQa/zs1NGqm7uo9nr9unLnPx2eR8AAAC0DgHdQ3aH9aC3T0Af2TtN5xzRL7h817trZa1tl/cCAABAyxHQPSQ0oPdthyEudX5y0gjFxfi70b/MOaA56/e123sBAACgZQjoHhI6xKW9etAlaWD3FF1w1MDg8t3vrlNtLb3oAAAAXkBA95CwHvQIV3Fp6PoThikp3v/tX727SG+t2N2u7wcAAIDmIaB7SG5RZfDznl3bN6D37Jqky6YPDi7f88F6Vftq2/U9AQAAcGgEdI8orawJTlKUGBejtMS4dn/Pq48dqq5J/vfZkleqFxfvaPf3BAAAwMER0D0ir6S+9zyrS6JMXS3EdpSeEq8fHjc0uHz/hxtUEfgjAQAAAO4goHtEXklV8POsLgmOve9l07OV1SVRkrSnqEL/WbDVsfcGAADA1xHQPSK0B717IDA7ISUhTjecOCy4/PCcjSquqHbs/QEAABCOgO4R4UNcnOtBl6TvHjlQAzL9ddcPlFXr0XlbHH1/AAAA1COge8T+sCEuzvWgS1JCXIx+ctKI4PJj8zZrf8gfDAAAAHAOAd0j3BriUufsI/ppRK8ukqTSKp8emr3J8TYAAACAgO4Zbg5xkaTYGKObZ40MLj/9+VbtLCg/yCsAAADQHgjoHpHn4hCXOieP6aUJAzMkSVW+Wt3/4XpX2gEAANCZEdA9Yn/YEBfne9AlyRijn51S34v+4uId2phb4kpbAAAAOisCukcUlteXNsxMcSegS9K0oVk6ZniWJKnWSvd8sM61tgAAAHRGBHQPsNaqoKw+oHdNjnexNQrrRX97xR6t2FHoYmsAAAA6l4gFdGNMf2PM48aYXcaYSmNMjjHmPmNMtxbuZ4Yx5rXA6yuMMduMMW8bY06NVFu9prTKp5paK0lKjo9VUnysq+05rH+GThvXO7h813trXWwNAABA5xKRgG6MGSppsaTLJC2UdK+kzZJulLTAGNO9mfu5RtI8SScGHu+V9Imk4yS9Y4y5NRLt9ZqCsvobRDNS3O09r3PTrBGKMf7P523I04JN+91tEAAAQCcRqR70hyX1lHSDtfYca+0vrLUnyB+wR0q641A7MMbES/qTpApJk6y1l1hrf2mtvUTSZEmVkm41xrhT4qQdhQ5vSXd5eEudYT3T9K2J/YPLd723VtZaF1sEAADQObQ5oAd6z2dJypH0UIPNv5NUKukSY0zqIXaVKSld0nprbdididbaNZLWS0qW1KWtbfaa0IDezcUbRBv68ckjlBDrP0W+2lagD9fkutwiAACAji8SPegzA4/vW2trQzdYa4slfSYpRdLUQ+wnV9I+SSOMMcNDNxhjRkgaLmmptbbDjbUoKPfeEBdJ6peRrIumDgwu/+W9dfLV0osOAADQniIR0OtKfjQ1q82GwOOIg+3E+sdPXBdo02JjzJPGmD8ZY56Sf3z7KknnN6dBxpjFjX1IGtWc1zsttAfdSwFdkq6bOUwpCf6bVtftLdbry3a63CIAAICOLRIBPT3w2FQtvrr1GYfakbX2BUknSCqQ9D1Jv5B0ifzDZP4t/42nHU5oDfT0ZO8McZH8s5r+YMbg4PJf3luvimqfiy0CAADo2DxVB90Yc7GkD+Wv4DJa/qExoyV9JOlBSc81Zz/W2kmNfUjyZL1AL1ZxCXXlsUOUmer/w2FnQbmeWpDjboMAAAA6sEgE9Loe8vQmttetLzjYTgLjzB+XfyjLJdbatdbacmvtWvl70RdLOt8Yc3xbG+w1YUNcPFLFJVRaUrxuOGFYcPnBjzeG/VEBAACAyIlEQK+ruNLUGPO6Gz6bGqNeZ5akeEmfNHKzaa2kuYHFSa1ppJcVlIeOQffWEJc6F04ZpOzuKZKkoooa/e3jjS63CAAAoGOKRECfHXicZYwJ258xJk3SdEllkj4/xH7q6pv3aGJ73foO13Vb6ME66A0lxMXo56fW32P71IIcbdtf5mKLAAAAOqY2B3Rr7SZJ70vKlr8KS6jbJKVK+o+1trRupTFmlDGmYUWVeYHH84wxh4VuMMYcIek8SVbSx21ts9cUV9YEP09LinOxJQd36rjemjSomySp2md113ueHNIPAAAQ1SJ1k+i18tcxf8AY82qgPOLHkn4i/9CWWxs8f03gI8hau1D+Si3Jkr40xjxnjPmzMeZ/kr6QlCTpfmvtqgi12TOKK+p70LsmebMHXZKMMfrV6aODy28u362l2wvcaxAAAEAHFJGAHuhFnyzpCUlTJN0kaaik+yVNbcHkQldIukzSAkmnBPZzsqRPJV1grf1JJNrrNcUV9T3oXTzcgy5JkwZ10+njeweX73xrjfwl7AEAABAJEUuD1trt8ofr5jzXNLHeyh/yn4hUu7zOWquSkCEuXRK9HdAl6ZZTRumD1XtV7bNamJOv91fv1Sljex/6hQAAADgkT9VB74zKq33y1fp7oBPjYpQQ5/1vSXZWqi6aMii4/Od31qraV3uQVwAAAKC5vJ8GO7iSitAbRL07/ryhG04crrRAb//mvFI9t3Cbyy0CAADoGAjoLiuqiI4KLg1lpibo2pn1kxfd9+GGsJtdAQAA0DoEdJeFhtpoCuiSdNn0bPVNT5Ik7S+t0iOfbHa5RQAAANGPgO6ykiipgd6YpPhY3XzKyODyo/M2a2dBuYstAgAAiH4EdJeFlViMggouDZ1zRD+N7dtVklRZU6u73mXyIgAAgLYgoLssWm8SrRMTY/Sbb4wJLr+2dJcWbz3gYosAAACiGwHdZUUhY9CjsQddkqYO6a7TxtXXQb/9zdWqrWXyIgAAgNYgoLssdIhL1ygbgx7qV6ePVkKs/3Raur1Ary/b5XKLAAAAohMB3WVhs4hGcUAfkJmiK44ZHFz+v3fWqqyq5iCvAAAAQGMI6C4LL7MYfWPQQ117/FBldUmUJO0pqqDsIgAAQCsQ0F0W1oMepWPQ66Qlxetnp4wILj8yd5N2UXYRAACgRQjoLiut9AU/j/aALknnTRqgMX38ZRcrqim7CAAA0FIEdJeVV9UH9OSEWBdbEhmxMUa/PbO+7OKrS3dpyTbKLgIAADQXAd1lpSE3UqYmRH8PuvT1sot/eIOyiwAAAM1FQHdZWUgPekpi9Peg1/nlaZRdBAAAaA0CustKKzteD7okDeyeostnUHYRAACgpQjoLuuoPeiSdN3M8LKLD83e6HKLAAAAvI+A7iJrbdgY9JT4jhXQ05Li9fNTRwaXH527RVvySl1sEQAAgPcR0F1UWVMrG7h3MjEuRnGxHe/b8a2J/TVhYIYkqcpXqz+8scrdBgEAAHhcx0uEUSR0/HlKByix2JiYGKM/nDVOxviXZ6/bp4/W7HW3UQAAAB5GQHdR2PjzDnSDaEPj+6fru0cODC7f9sZqVVT7DvIKAACAzouA7qKwGugd7AbRhn52ykilJ8dLkrbll+lf8za73CIAAABvIqC7qLSyc/SgS1JmaoJunjUiuPzg7I3aWVDuYosAAAC8iYDuorJO1IMuSRdOGaTRfbpKkiqqa3XHW6tdbhEAAID3ENBd1FnGoNeJjTH6w9ljg8tvr9ijzzbmudgiAAAA7yGguyisB72DVnFp6MjsTJ1zRN/g8u9eX6VqX62LLQIAAPAWArqLQsegJ3eCHvQ6vzx9dPAPko25JXpyfo67DQIAAPAQArqLOmMPuiT16pqkG04cHly+78MN2ltU4WKLAAAAvIOA7qKwKi6JnacHXZIumz5YQ3ukSpJKKmt0+5vcMAoAACAR0F3VWXvQJSkhLka3nz0uuPzm8t2at2Gfiy0CAADwBgK6i0KruCR3soAuSdOGZYXdMPqbV1cywygAAOj0COguqqiur16SHN/5Arok3XrGGKUl+Yf35Owv09/nbHK5RQAAAO4ioLsotLc4qZMG9B5pibrllJHB5b/P2aQteaUutggAAMBdBHQXlYcE9M7agy75Zxg9rH+6JKnKV6vfvrZS1lqXWwUAAOAOArqLQnvQO+MY9DqxMUZ3nDNeMca/PG9Dnt5cvtvdRgEAALiEgO6i8rAhLp37WzG+f7oumToouHz7m6tVVFHtYosAAADc0blToctCbxLtrGPQQ910ykj1SEuUJOUWV+qe99e73CIAAADnEdBdVMEY9DBdk+L1m2+MCS4/tSBHK3cWutgiAAAA5xHQXVReRRWXhs48rI9mDMuSJNVa6dZXVshXyw2jAACg8yCgu6iihh70howx+sPZY5UQ6z81l+0o1JPzc9xtFAAAgIMI6C4q7+QziTZlSI8uunbm0ODyX95fpx0HylxsEQAAgHMI6C6prbWqrKm/STQxjm9FqGuOH6phPbtIksqqfPrNq9RGBwAAnQOp0CWh4TwpPkbGGBdb4z2JcbH6v3PHB5dnr9unN6iNDgAAOgECukvCa6AzvKUxk7MzdfHUgcHlP7yxSgVlVS62CAAAoP0R0F1STonFZrnl1FHq1dVfGz2vpEp3vLXG5RYBAAC0LwK6S6iB3jxdk+L1h7PHBZdfWLxD8zfmudgiAACA9kVAd0loBZdEAvpBnTK2t04b1zu4/MtXVoT9gQMAANCRENBdUhlWA51vw6HcdtZYpSXFSZK27i/T/R9tcLlFAAAA7YNk6JLyqvoqLtRAP7SeXZP0y9NGB5f/OXezVu8qcrFFAAAA7YOA7pKwKi5xBPTm+O6RA3RUdqYkyVdr9cuXl8tXS210AADQsRDQXRI6hjqJHvRmiYkxuvPc8UqI9Z+2y3YU6vFPt7jcKgAAgMgioLuEMoutM6xnF11/wrDg8l/eX6fN+0pcbBEAAEBkEdBdEtaDzk2iLXL18UM1uk9XSf4ZWX/+0nLVMtQFAAB0ECRDl1AHvfXiY2P0l/MPU1yMkSR9mXNATy7IcbdRAAAAEUJAd0lFdX0VlyQCeouN7Zuua48fGlz+87trtXV/qYstAgAAiAwCuktC66AnxvFtaI0fnTBcI3ulSfL/wXPLiwx1AQAA0Y9k6JLKkB70RMostkpCXIzuPv8wxQaGunyxJV/PfLHV5VYBAAC0DQHdJZU1IQGdm0Rb7bD+GfrhsUOCy396Z62255e52CIAAIC2IRm6hCEukXPDicM1rGcXSVJZlU+/eHm5rGWoCwAAiE4kQ5eE9aAzxKVNkuJjdfd5hykw0kWfbdyvZxdud7dRAAAArURAd0n4GHS+DW01YWA3XXlM/VCXO99eo50F5S62CAAAoHVIhi4JG+LCGPSI+MnJIzQkK1WSVFJZo59T1QUAAEQhkqFLqnwMcYm0pPhY3X3+YTKBoS6fbszT01R1AQAAUYaA7pLQIS4JDHGJmEmDMnXVseFDXTbvK3GxRQAAAC1DMnRJ+E2ifBsi6acnjwibwOimF5apJuQ/FgAAAF5GMnRJeJlFhrhEUmJcrO75zuGKj/WPdflqW4EembvZ5VYBAAA0DwHdJfSgt6+xfdP145NGBJfv/WC9Vu4sdLFFAAAAzUMydElYmUWquLSLHx47RBMGZkiSamqtbnp+mSqqfQd/EQAAgMtIhi5hiEv7i4uN0T3fPkLJ8f7ju25vse79YL3LrQIAADg4ArpLqhji4ojBWan61emjgsv/nLdZC7fku9giAACAgyMZuiR0DDplFtvXxVMH6ZjhWZIka6WbXliqksoal1sFAADQOJKhC2p8taoJzHAZY6S4GONyizo2Y4zuOu8wpSXFSZK255frjrdWu9wqAACAxhHQXdBwFlFjCOjtrU96sm4/e1xw+dmF2/Xeqj0utggAAKBxBHQXUMHFHWcf0VdnjO8TXP7FS8u1t6jCxRYBAAB8XcTSoTGmvzHmcWPMLmNMpTEmxxhznzGmWyv2NdEY819jzI7AvvYaYz4xxnwvUu11EzXQ3WGM0R3fHKc+6UmSpANl1brp+WWqDQw3AgAA8IKIpENjzFBJiyVdJmmhpHslbZZ0o6QFxpjuLdjXjyR9KWmWpI8k/VXSK5JiJZ0eifa6jRKL7slISdA93z5CdaOKPt2Yp8c+3eJuowAAAELERWg/D0vqKekGa+3f6lYaY+6R9BNJd0i6+lA7McbMkvSApA8knWetLW6wPT5C7XUVJRbddfTQ7rr6uKH6+5xNkqS73luracO6a2zfdJdbBgAAEIEe9EDv+SxJOZIearD5d5JKJV1ijEltxu7ullQu6cKG4VySrLXVbWutN4QNcWEMuit+ctIIHdbfH8irfVY3PPuVyquYZRQAALgvEulwZuDxfWttbeiGQMj+TFKKpKkH24kxZpykwyS9LynfGDPTGHOzMeYmY8yJxpgOk2RDh7gkxHaYLyuqJMTF6L7v1M8yumlfqe54m9KLAADAfZFIhyMDj03Nob4h8DjiEPs5MvCYK2mOpI/l71H/i6QPJS01xgxrToOMMYsb+5A06pAvdkBYFRfGoLtmSI8u+v1ZY4LLT3++TR+s3utiiwAAACIT0OsG7hY2sb1ufcYh9tMz8HiFpGxJZwT2PULS05LGS3rLGJPQ2oZ6BUNcvOPbkwfo1LG9g8s/f2m5cim9CAAAXOSldFjXllhJ37XWvm2tLbLWbpD0PUmL5A/r3zrUjqy1kxr7kLS23VrfAuFVXLz0Leh8jDH6v2+NV++u/tKL+aVVuukFSi8CAAD3RCId1vWQN1UCo259wSH2U7d9j7V2QegGa62V9Fpg8agWts9zwuugM8TFbRkpCbrnO4cHSy/O25CnR+dtdrdRAACg04pEQF8XeGxqjPnwwGNTY9Qb7qegie0HAo/JzWuWd4WPQacH3QumDc3SVccOCS7f/d46Ldl24CCvAAAAaB+RSIezA4+zGlZaMcakSZouqUzS54fYz+fyl2TMbqIk47jAY9TPKlPpqw/o8QR0z7h51khNGJghSaqptbr+v1+psKxDVPYEAABRpM3p0Fq7Sf7SiNmSrmuw+TZJqZL+Y60trVtpjBlljAmrqGKtLZP0mKQkSX80pm7AgWSMGS/p+5JqJL3Y1ja7rTpkiAtlFr0jPjZGD3x3grom+efv2llQrp+/tFz+EVYAAADOiFQ6vFb+8ogPGGNeNcb8yRjzsfyziK6XdGuD568JfDT0G0lLJf1Y0gJjzF+NMU9L+kL+4H5z4A+CqFblY4iLVw3ITNFd5x0WXH531R795/OtLrYIAAB0NhFJh4HQPFnSE5KmSLpJ0lBJ90uaaq3d38z9FEk6RtKdkjIl/UjSNyR9KukUa+39kWiv20J70OPpQfecU8f10aVHDwou//HNNVq1q6kqogAAAJEVsXRord1urb3MWtvHWptgrR1krf2xtfZrd9pZa4211jSxnxJr7a3W2hHW2kRrbYa1dpa19v1ItdVtoT3oCfSge9IvTx+tMX26SvJ/v370369UUlnjcqsAAEBnQDp0QWhApwfdm5LiY/XghROUmuAvg7klr1S/fmUF49EBAEC7Ix26oKqGHvRoMKRHF9157vjg8qtLd+mFxTtcbBEAAOgMSIcuqA4d4hLb6EgfeMTZR/TTtyf3Dy7/9rWV2rC32MUWAQCAjo6A7oIqbhKNKr8/a6yG9+wiSaqortU1zyxRKePRAQBAOyEduqDaVz+OmSEu3peSEKcHL5yopHj/92pjbol+8TLj0QEAQPsgHbqAm0Sjz8jeabrzm/Xj0d9YtktPLaA+OgAAiDzSoQu4STQ6nTuxvy44amBw+Y9vrdaSbV+rIgoAANAmpEMXhN8kyrcgmvzuzDEa189fH73aZ/WjZ5Yov7TK5VYBAICOhHToAm4SjV5J8bH6+0WTlJ4cL0naVVihG5/7Sr5axqMDAIDIIB26oJqZRKPagMwU3fudw4PL8zbk6YGPNrjYIgAA0JGQDl0Q3oNOHfRodMKoXrpu5tDg8gMfb9CcdbkutggAAHQUBHQXVFFmsUP4yUkjdPSQ7pIka6Uf/2+pdhaUu9wqAAAQ7UiHLuAm0Y4hLjZGD1wwQT3TEiVJBWXVuvbpxaqo9rncMgAAEM1Ihy7gJtGOo0daoh66aKJiY/xDlZbtKNRvX1vJJEYAAKDVSIcu4CbRjuXI7Ezdevro4PLzi3bo6S+2udgiAAAQzUiHLqAHveO5bHq2zp3QL7h82+ur9GVOvostAgAA0Yp06IIqetA7HGOM7jx3fHASo5paq2ueXqI9hRUutwwAAEQb0qELQnvQuUm040iKj9U/Lp6kzNQESVJeSaWufnqxKmu4aRQAADQf6dAFjEHvuPp3S9GDF04I3jS6dHuBfvfaKm4aBQAAzUY6dJiv1qpuVvgYo2CQQ8cxbWiWfhVy0+hzX27Xfxdy0ygAAGgeArrDuEG0c7h8erbOOaJvcPn3r6/S4q3cNAoAAA6NhOgwbhDtHIwx+tO5h2lsX/9No9U+q6ufXqLdhcw0CgAADo6E6DBuEO08khP8N412S4mXJO0rrtRVTy1WeRU3jQIAgKaREB3GDaKdy4DMFD100UTFBe41WLGzUDe/uIybRgEAQJNIiA5jDHrnM21oln5/1tjg8lvLd+v+jza42CIAAOBlJESHhfagx8dSwaWzuHjqIF169KDg8n0fbtBby3e72CIAAOBVBHSHVYaOQY+LdbElcNpvvjFGM4ZlBZdvemGpVuwodLFFAADAiwjoDgsbg04PeqcSFxujhy6cqMFZqZKkiupaXfnUIuUWVbjcMgAA4CUEdIdV++pvDuQm0c4nPSVe/7p0stKS4iRJe4oqdOV/FquimsouAADAj4ToMG4SxdAeXfTQhRODs8gu216gn7+0nMouAABAEgHdceE3iXL4O6tjR/TQb84YHVx+bekuPTR7o4stAgAAXkFCdFgVAR0Bl07L1gVHDQwu/+X99Xpt6U4XWwQAALyAhOiw8ImKuEm0MzPG6LazxuroId2D6372wnJ9mZPvYqsAAIDbCOgOqwm5STQuhsPf2SXExegfF0/S0B7+yi5VPn9lly15pS63DAAAuIWE6DCGuKCh9JR4/fv7R6l7aoIkqaCsWpf9e6HyS6tcbhkAAHADCdFhoT3ozCSKOgO7p+hfl05WYqD0Zs7+Ml311CLKLwIA0AkR0B1WU1vfgx5HQEeICQO76b7vHCETOC0WbT2gn724XLW1lF8EAKAzIaA7LLQOOmPQ0dBp4/vol6eNCi6/sWyX/vrBOhdbBAAAnEZCdFhNLTOJ4uCuPGaILppSX37xodmb9PyX211sEQAAcBIJ0WE1vtAedIa44Ovqyi8eP7JHcN2vXlmhuev3udgqAADgFAK6w6pDyyxSxQVNiIuN0YMXTtToPl0l+f/zcvXTi7V8R4G7DQMAAO2OhOiwsImKuEkUB9ElMU6Pf3+y+qQnSZLKqny6/IkvtXU/NdIBAOjICOgOCx2DTg86DqVPerKeuvwopSfHS5LySqr0vccXKq+k0uWWAQCA9kJCdFg1Y9DRQsN7pemxkBrpW/eX6bJ/f6nSyhqXWwYAANoDAd1h4RMVcfjRPJOzM/XABRNU9zfdip2FuvrpxWFlOwEAQMdAQnRYaA86AR0tccrY3rr9nHHB5Xkb8vTzl5jICACAjoaE6LDwKi4McUHLXDRlkG48cXhw+ZWvdurP7611sUUAACDSCOgOq6kN7UEnoKPlfnzScF1w1IDg8iOfbNZjn25xsUUAACCSCOgOY4gL2soYo9vPHqeTRvcKrrv9zdV6afEOF1sFAAAihYTosLAhLjEcfrROXGyM/nbBBE0a1C247paXluv9VXtcbBUAAIgEEqLDanwMcUFkJCfE6vFLj9So3mmSJF+t1Y+e/UrzN+W53DIAANAWBHSHMVERIik9JV5PXXGUBnVPkSRV1dTqyicXaen2AncbBgAAWo2E6LDQutX0oCMSeqYl6ekrpqhX10RJUmmVT9//90Kt31vscssAAEBrENAdFtqDzk2iiJQBmSl6+oop6pYSL0kqKKvWJY99oe35ZS63DAAAtBQJ0WGhY9DjYuhBR+QM75WmJy47SqkJsZKkvUWVuuhfXyi3qMLllgEAgJYgoDsstIoLPeiItMMHZOhflx6phDj/ubUtv0yXPLZQBWVVLrcMAAA0FwnRYdRBR3s7emh3PXzhRMUG/kOzbm+xvvf4QhVVVLvcMgAA0BwkRIeFV3FhiAvax0ljeukv5x8WXF6+o1Dff3yhSiprXGwVAABoDgK6w8J60JmoCO3omxP6645vjgsuL9lWoMuf+FLlVT4XWwUAAA6FhOiwsIAeRw862tdFUwbp92eOCS4v3JKvK59apIpqQjoAAF5FQHdYTchNonH0oMMB358+WL86fVRw+dONebr66cWqrCGkAwDgRSREh4XfJEoPOpxx1bFDdfOsEcHlOev26Uf//SrsfAQAAN5AQHdY+E2iHH4450cnDNf1JwwLLn+weq9ufO6rsNr8AADAfSREh9GDDjf99OQR+uGxQ4LLb6/Yo5teWCZfyB+OAADAXQR0B1lrwycqYgw6HGaM0S9OG6XLpmcH1722dJd+RkgHAMAzSIgOCg1AMUaKiaEHHc4zxui33xiji6YMDK57+aud+unzSxnuAgCABxDQHRTWe874c7jIGKPbzx6nC44aEFz32tJd+snzywjpAAC4jJTooOra0PHnHHq4KybG6I5zxof1pL+xbJdufG4p1V0AAHARKdFBYTXQuUEUHhATY/THc8bpe0cPCq57a8Vu3fAsJRgBAHALAd1BoUMHmKQIXmGM0W1njdX3p2UH172zco+ue2aJqmoI6QAAOI2U6KCqkICeQA86PMQYo9+dOUZXzBgcXPf+6r269pklzDgKAIDDCOgOCh/iwqGHtxhj9OszRuuqkDrpH67Zq2ueJqQDAOAkUqKDakJuEmUMOrzIGKNfnjZKVx83NLju47W5+sGTi1RWVeNiywAA6DwI6A5ikiJEA2OMfn7qSF03sz6kz9uQp+89tlBFFdUutgwAgM6BlOig0KoY8XH0oMO7jDG6edZI3XTyiOC6RVsP6MJHP1d+aZWLLQMAoOOLWEA3xvQ3xjxujNlljKk0xuQYY+4zxnRrwz6PNcb4jDHWGPPHSLXVLaE96FRxgdcZY3T9icP1m2+MCa5bubNI33lkgfYWVbjYMgAAOraIpERjzFBJiyVdJmmhpHslbZZ0o6QFxpjurdhnmqQnJZVFoo1e4KsNnUmUHnREhytmDNafvzVeJnDKbsgt0fn/WKDt+R3mRxMAAE+JVDfuw5J6SrrBWnuOtfYX1toT5A/qIyXd0Yp93i8pXdKfItRG14XWQY+NIaAjenznyIG6/7sTFBc4b7fll+n8fyzQxtwSl1sGAEDH0+aAHug9nyUpR9JDDTb/TlKppEuMMakt2OfZ8vfG3yBpV1vb6BU1tQxxQfQ66/C++sfFk5QQ5z939xRV6DuPLNDqXUUutwwAgI4lEilxZuDxfWtt2LSD1tpiSZ9JSpE0tTk7M8b0lPSopFettU9HoH2eETrEhTKLiEYnjemlf3//SCXHx0qS9pdW6bv/XKBFOfkutwwAgI4jEgF9ZOBxfRPbNwQeRzSxvaFH5W/X1a1tkDFmcWMfkka1dp+REN6DTkBHdJo+LEtP/+AopSXFSZKKKmp00b++0Edr9rrcMgAAOoZIBPT0wGNhE9vr1mccakfGmMslnSXpWmtth/tt76tlDDo6hkmDMvXslVOV1SVBklRZU6ur/rNYLyza7nLLAACIfp4ZCG2MyZZ0n6QXrLXPt2Vf1tpJjX1IWhuBprYaY9DRkYzrl64Xr56mAZnJkvxDuH724nI98skml1sGAEB0i0RKrOshT29ie936gkPs53FJ5ZKujUCbPKkmpA46PejoCLKzUvXS1dM0uk/X4Lo/vbNWd7y1WrUhf5ACAIDmi0RAXxd4bGqM+fDAY1Nj1OtMlL9U477AxETWGGMl/Tuw/dbAulfb1FoXMQYdHVHPrkn63w+nasrgzOC6R+dt0c0vLAubPRcAADRPXAT2MTvwOMsYExNaySUw2dB0+Scb+vwQ+3lK/movDQ2XdKykpfJPhvRVWxvsltAx6FRxQUfSNSleT15+lG587iu9t8p/+8jLX+1UflmVHr5oolISInGpAQCgc2hzD7q1dpOk9yVlS7quwebbJKVK+o+1trRupTFmlDEmrKKKtfYGa+0PGn6ovgf9rcC6hrXWo0ZoD3osY9DRwSTFx+rhiybpgqMGBNfNWbdPF/3rC+WXVrnYMgAAokukUuK1knIlPWCMedUY8ydjzMeSfiL/0JZbGzx/TeCjU/ExxAUdXGyM0Z3fHK8bThgWXPfVtgKd+/BnyskrPcgrAQBAnYgE9EAv+mRJT0iaIukmSUMl3S9pqrV2fyTeJ9pVc5MoOgFjjH46a6RuO2usTOA0z9lfpnP/Pl+Ltx5wt3EAAESBiI2zsNZut9ZeZq3tY61NsNYOstb+2Fr7td/I1lpjrW1WQrXWPhF4/q8j1Va3hI1BJ6Cjg7t0Wrb+ftEkJcb5LzP5pVW68NHP9e7K3S63DAAAb2MgtIPCqrjEcujR8Z06rreevWqqMlPrJzS65pkleuzTLS63DAAA7yIlOsjnYww6Op+JA7vp5WumaXBWqiTJWun2N1fr96+vCrsvAwAA+BHQHRRexYWAjs4jOytVL10zTZMGdQuue2J+jq55erHKq3wutgwAAO8hoDuohjHo6MQyUxP0zA+m6LRxvYPr3l+9Vxc8+rnySipdbBkAAN5CQHdQWA86ExWhE0qKj9VDF07UlccMDq5bur1AZz/4mdbuKXKxZQAAeAcB3UGhY9DjmagInVRMjNGtZ4zRbWeNVd0/knYWlOtbD8/XR2v2uts4AAA8gJToIMagA/UunZatf106WakJsZKk0iqffvDUIj06d7Os5eZRAEDnRUB3UNhMogxxAXTCqF56+drp6t8tWZK/wssdb6/Rz19arqqa2kO8GgCAjomA7qDQm0TpQQf8RvZO06vXTdfkkAovzy/aoYsf+0L5pVUutgwAAHcQ0B1UQx10oFFZXRL1zJVTdO7EfsF1C7fk65yHPtOGvcUutgwAAOcR0B0UNsSFm0SBMIlxsfrr+YfrF6eNkgn8/botv0znPjxfs9fmuts4AAAcREp0UA1j0IGDMsbo6uOG6h8XT1JK4ObR4soaXf7kl3po9kZuHgUAdAoEdAf5qOICNMspY3vrxaunqW96kiT/zaN3v7dO1/13iUora1xuHQAA7YuA7qBqHzOJAs01pm9XvX79DB01ODO47u0Ve3Tuw/O1dX+piy0DAKB9EdAdFN6DzqEHDiWrS6Ke+cEUXXr0oOC6dXuLddaDn+mT9ftcbBkAAO2HlOigsDHo9KADzRIfG6Pbzh6nu847TAmx/ktWYXm1Lvv3Qv3jk02MSwcAdDgEdAcxURHQet+ePEDPX320enf1j0uvtdL/vbNW1z/7lcqqGJcOAOg4COgOYqIioG2OGJCh16+friOz6yc1enP5bp378Hxt3lfiYssAAIgcArqDqIMOtF3PtCQ984OpunjqwOC6tXv849LfWbHbxZYBABAZpEQHVfsoswhEQkJcjP54znj9+VvjlRDnv4yVVNbommeW6PY3V4dVTAIAINoQ0B3k4yZRIKK+c+RAvXzNNA3ITA6ue+zTLfruPz/XnsIKF1sGAEDrEdAdxEyiQOSN65euN390jE4a3Su4bvHWAzrjgXn6bGOeiy0DAKB1COgO8tWGTlTEoQciJT0lXv+8ZJJ+cdoo1f1zan9plS557As9+PEG1dZSihEAED1IiQ6qqWUMOtBeYmKMrj5uqP575VRldUmU5C/F+Jf31+uKJ7/UgdIql1sIAEDzENAdVONjDDrQ3qYO6a63b5ihowZnBtfNXrdPp90/T19s3u9iywAAaB4CuoN89KADjujZNUn//cEU/fDYIcF1e4oqdMGjn+v+DzeE/SwCAOA1BHQHhU5UFB/LoQfaU1xsjH55+mg9dulkdUuJl+Qf8nLvh+t14aNUeQEAeBcp0UH0oAPOO3F0L7194zFhQ16+2JKv0x+Yp4/X7nWxZQAANI6A7qAa6qADruiTnqxnr5yqG08cHqzykl9apcufWKQ/vrlaVTVMbAQA8A4CuoNCbxKNpQ464KjYGKOfnDxCz/xgqnp1TQyu/9enW3TeP+YrJ6/UxdYBAFCPgO6gmrA66AR0wA1HD+2ud248VieM6hlct3xHoc54YJ6eX7Rd1nIDKQDAXQR0B/nChrhw6AG3ZKYm6LFLJ+s33xij+MB/s0qrfLrlxeW65ukl1EwHALiKlOggxqAD3mGM0RUzBuvla6ZrSFZqcP27q/bolPvmau76fS62DgDQmRHQHVJba1X3n3Nj/LMeAnDf+P7pevOGGbp46sDgutziSn3v8YX6/eurVFHtc7F1AIDOiIDukGrGnwOelZIQpz+eM16Pf3+ysrokBNc/MT9HZz34qVbvKnKxdQCAzoaA7hBqoAPed8KoXnr3x8fqpNH1N5Cu31uicx76TP+cu4kZSAEAjiCgOyR0/Hk8N4gCnpXVJVGPfm+y7vjmOCXF+39Wq3y1uvPttfruPxdQjhEA0O5Iig7xUQMdiBrGGF00ZZDeuuEYHdY/Pbj+y5wDOu3+eXpyfo5q6U0HALQTArpDqOACRJ+hPbropWum6YYThweHppVX+/S711fpon99oe35ZS63EADQERHQHcIYdCA6xcfG6Kcnj9Cr107XyF5pwfULNu/XqffN1TNfbGVyIwBARBHQHVLtC63iwmEHos34/ul6/frpuvb4oar7G7u0yqdbX1mp7z2+ULsKyt1tIACgwyApOoQedCD6JcbF6pZTR+nla6draI/6yY3mbcjTKffO1f++3EZvOgCgzQjoDgkbg85NokBUO2JAht664RhddewQmcCPc3FljX7+0gpd8thCbdvP2HQAQOsR0B3i4yZRoENJio/Vr04frRevPlqDs+p70z/dmKdZ932iR+duVk3I0DYAAJqLgO6QmpCZRGMZgw50GJMGZertG47RlccMDo5Nr6iu1R1vr9G5f5+vNbuZhRQA0DIkRYfU+OhBBzqq5IRY3XrGGL1y7XSN6l1f6WX5jkKd+bdPdfd7a1VR7XOxhQCAaEJAd0gNN4kCHd7hAzL0xvUz9LNTRiohzn95ram1emj2Jp1+/zx9sXm/yy0EAEQDArpDQsegx3OTKNBhxcfG6LqZw/TOjcfoqOzM4PrNeaX6zj8/169eWaHCsmoXWwgA8DoCukPCx6AT0IGObmiPLnruqqm645vjlJYYF1z/3y+26cR75uiVr3ZQkhEA0CgCukPCq7hw2IHOICbG6KIpg/TBT4/TSaN7BdfnlVTpJ/9bpgse/Vwbc4tdbCEAwItIig4JvUmUHnSgc+mdnqRHvzdJ/7h4kvqkJwXXf745X6fdP093vbtW5VXcRAoA8COgO6SGOuhAp2aM0anjeuvDnx6nq44dEvxDvdpn9fCcTTr53k/08dq9LrcSAOAFBHSH+ELGoDOTKNB5pSbG6Venj9ab18/QpEHdgut3HCjX5U8s0g//s0i7CspdbCEAwG0EdIfUMAYdQIjRfbrqhR8erT9/a7wyUuKD699btVcn/vUTPTR7I7XTAaCTIik6xEcddAANxMQYfefIgfr4puP17cn9g+vLq326+711mnXvXH2wei/VXgCgkyGgO6SamUQBNCEzNUF3nXe4Xrj66LCZSLfll+nKpxbp0n9/qY25JS62EADgJAK6Q3zUQQdwCEdmZ+rN62fo9rPHKj25ftjL3PX7dOp9c3XHW6tVVMEkRwDQ0RHQHRI2Bj2Www6gcXGxMbrk6GzNufl4XTx1oOr+nq+ptXp03had8Jc5en7RdtXWMuwFADoqkqJDfJRZBNAC3VIT9MdzxuuN62foqOzM4Pq8kird8uJyffPhz7QoJ9/FFgIA2gsB3SFMVASgNcb2Tdf/fjhVf7tgQtgkR8t2FOq8fyzQNU8v1tb9pS62EAAQaQR0h1DFBUBrGWN05uF99dFNx+n6E4YpIa7+0v3Oyj066Z5PdPubq1VYxvh0AOgICOgOYSZRAG2VkhCnm2aN1Ec/PU5nHt43uL7aZ/XYp1t07N2z9dinW1RVU3uQvQAAvI6A7pDakDrGMQR0AG0wIDNFf7tggl65dlrYbKSF5dW6/c3VmnXvJ3p35W7qpwNAlCKgOyRsiIshoANouwkDu+nFq4/WwxdN1MDMlOD6nP1luvrpJfr2Iwv01bYDLrYQANAaBHSH1DAGHUA7MMbo9PF99MFPj9WvzxitrklxwW1f5hzQNx+erx/+Z5E25ha72EoAQEsQ0B1SS0AH0I4S42L1g2OGaO4tM3X59MGKj62/zry3aq9m3TtXP3thmXYWlLvYSgBAcxDQHeKzBHQA7S8jJUG/PXOMPvhJ+I2ktVZ6YfEOzbx7jm5/c7XyS6tcbCUA4GAI6A4J7UGPYQw6gHaWnZWqv10wQW9eP0PHjegRXF/lq/VXfLlrth74aINKK2tcbCUAoDEEdIcwkygAN4zrl64nLz9Kz145VRMGZgTXl1TW6J4P1uvYu2br359tUUW1z71GAgDCENAdEnqTKGUWATjt6KHd9fI10/TPSyZpeM8uwfX7S6t02xurdfzdc/SfBTmqrCGoA4DbCOgOCa2DHks+B+ACY4xmje2td398rO4+7zD1y0gObttTVKHfvLZKM++eo2e+2MpkRwDgIgK6Q3xUcQHgEbExRudPHqCPbz5OvztzjHqkJQa37Sqs0K2vrNTMv8zRcwu3qdpHUAcApxHQHcJMogC8JjEuVpdNH6x5t8zUr88YrawuCcFtOwvK9YuXV+iEv87R84u2q4agDgCOIaA7pMbHTaIAvCkpvr6G+q9OH6XM1Pqgvj2/XLe8uFwn3vOJXli0nR51AHAAAd0hoXXQKbMIwItSEuJ01bFDNe+Wmfr5qaOUkRIf3LZ1f5l+9uJy/82kn2+l6gsAtKOIBXRjTH9jzOPGmF3GmEpjTI4x5j5jTLdmvj7VGHORMea/xpi1xphSY0yxMWaRMeYmY0zCoffiXcwkCiBapCbG6Zrj/UH95lkjlJ5cH9R3FpTrN6+u1LF3zdajczdTRx0A2kFEAroxZqikxZIuk7RQ0r2SNku6UdICY0z3ZuzmGElPSzpF0kpJf5P0X0n9JP1F0mxjTFIk2uuGkBEuBHQAUSEtKV4/OmG45v18pn52ykh1C+lRzy2u1B1vr9GMP3+sBz7aoMLyahdbCgAdS6R60B+W1FPSDdbac6y1v7DWniB/UB8p6Y5m7GOPpIsl9bHWnhfYxw8ljZC0RNI0SddFqL2O89XWj9skoAOIJl2T4nXdzGH67Bcn6NdnjFavrvVVXw6UVeueD9Zr+v99rD+/u1Z5JZUuthQAOoY2B/RA7/ksSTmSHmqw+XeSSiVdYoxJPdh+rLVLrbXPWGurGqwvlvTXwOLxbW2vW8LKLDIGHUAUSkmIC95Mesc3x2lAZn0d9ZLKGv19zibN+PPH+v3rq7Q9v8zFlgJAdItED/rMwOP71tqw2/sD4fozSSmSprbhPer+dxq1gx1DCx9QZhFANEuMi9VFUwZp9k3H655vH65hITOTVlTX6on5OTru7tn60X+XaMWOQhdbCgDRKS4C+xgZeFzfxPYN8vewj5D0USvf4/LA47vNebIxZnETm0a18v3bLHwmUQI6gOgXFxujcyf21zlH9NN7q/bowdkbtWpXkSSp1kpvLt+tN5fv1tFDuuuq44bo+BE9ZLj+AcAhRSKgpwcem+omqVuf0ZqdG2N+JOlUSUslPd6afXhB2BCXWH5BAeg4YmKMThvfR6eO6615G/L0z7mb9enGvOD2BZv3a8Hm/RrZK01XHjtEZx3eVwlxVPkFgKZEIqC3G2PMuZLuk/8G0m9Za5tVJsBaO6mJ/S2WNDFiDWwBxqAD6OiMMTp2RA8dO6KHVu4s1KPzNuvN5buD1791e4t18wvLdPd7a3XZ9MG6cMpAdU2KP8ReAaDziUQXRl0PeXoT2+vWF7Rkp8aYcyQ9JylX0vHW2s2taZxX+KiDDqATGdcvXfd/d4I++dnxunz6YKUkxAa37S2q1P+9s1ZH3/mRfvfaSm3eV+JiSwHAeyIR0NcFHkc0sX144LGpMepfY4w5X9ILkvZKOs5au+4QL/E8ZhIF0Bn175ai3545Rgt+caJuOXWkeqTVl2gsrfLpyQVbdcJfP9Fl/16ouev3yYZcKwGgs4rEEJfZgcdZxpiY0Eouxpg0SdMllUn6vDk7M8ZcJOlJSTslzYz2nvM6oTOJxjEGHUAnk54Sr2uPH6YrZgzWa1/t0r8+3az1e+t7zmev26fZ6/ZpWM8u+v60bJ07sZ9SEjw9ChMA2k2be9CttZskvS8pW1+fSOg2SamS/mOtLa1baYwZZYz5WkUVY8ylkp6StE3SsR0lnEtSTS096ACQGBerbx85QO/9+Fg9fcUUnTS6p0IviRtzS/TrV1dq6p0f6c6312jHAeqpA+h8ItU9ca2k+ZIeMMacKGmNpCny10hfL+nWBs9fE3gMXpaNMTPlr9ISI3+v/GWNlOMqsNbeF6E2OyqszCJj0AF0csYYzRiepRnDs5STV6onF+TohUU7VFLpn+6iqKJG/5y7Wf+at1mzxvTWJUcP0rSh3SnTCKBTiEhAt9ZuMsZMlvQH+Usini5pt6T7Jd1mrT3QjN0MUn2P/uVNPGer/FVdog5VXACgcdlZqfrdmWP105NH6MXFO/Tk/Bzl7Pf3nNda6d1Ve/Tuqj0akpWqC6cM1PmTBig9heovADquiA3ws9Zul3RZM5/7tYRqrX1C0hORao/XhAb0GMr/AsDXpCXF67Lpg3Xp0dmasz5X//4sR/M21NdT35xXqj++tUZ3v7dOZx7eVxdPHaTD+6fTqw6gw+EOHIeEBvQ4EjoANCkmxuiEUb10wqhe2rC3WM98sU0vLd6h4sDwl8qaWr24eIdeXLxD4/p11cVTBumsI/pyUymADoOk6BBf2Bh0FxsCAFFkeK80/f6ssfri1hP1p3PHa2zfrmHbV+4s0i9eXqEpd36k37++Suv3FrvUUgCIHLobHFJLFRcAaLWUhDhdcNRAfffIAVq6vUBPf75Nby7fpcoaf2Xf4ooaPTE/R0/Mz9GEgRn67pEDdMZhfdUlkV9zAKIPVy6H+KjiAgBtZozRhIHdNGFgN/3mG6P14uIdeuaLbdqSF6zkq6+2FeirbQW67Y3V+sZhffSdIwdq4sAMxqoDiBoEdIf4fAR0AIikjJQE/eCYIbp8+mAt2Lxfz3yxVR+s3qvqwPW2rMqn5xft0POLdmhYzy76zuQB+ubEfsrqkniIPQOAuwjoDqEHHQDaR0yM0fRhWZo+LEv7Syr1ylc79b8vt2tDbv1MpRtzS3TH22v053fX6uQxvfTtIwfo2OE9uB4D8CQCukN8tfWfUwcdANpH9y6J+sExQ3TFjMH6anuB/rdwu95YvktlVT5J/lmd31m5R++s3KNeXRN1zoR++tbE/hrRK83llgNAPQK6Q0JnEo2hxwYA2pUxRhMHdtPEgd302zPH6K3lu/Xcl9u0ZFtB8Dl7iyr1yCeb9cgnmzWuX1d9c0J/nXV4X/VIYwgMAHcR0B0SXgedgA4ATklNjNO3jxygbx85QBv2Fuv5Rdv18pKd2l9aFXzOyp1FWrlzte58e42OG9FD507sp5NG91JSfKyLLQfQWRHQHRI+kygBHQDcMLxXmm49Y4xuOXWU5q7fp5eX7NQHa/aqKlCu0Vdr9fHaXH28NldpiXE647A+Ondif00e1I1rNwDHENAdEhrQGYMOAO6Kj43RiaN76cTRvVRYXq23V+zWy0t26MucA8HnFFfW6Lkvt+u5L7erf7dknXl4X515WF+N7pNGyUYA7YqA7hCquACAN6Unx+uCowbqgqMGatv+Mr3y1U69/NUObd1fFnzOjgPl+vucTfr7nE0a1rOLzjysr848vI+G9OjiYssBdFQEdIcwkygAeN/A7im68aThuuHEYVqy7YBeWrJTby7bpaKKmuBzNuaW6N4P1+veD9drXL+uOuvwvjrjsL7ql5HsYssBdCQEdIfUcJMoAEQNY4wmDcrUpEGZ+t2ZYzR3fZ7eWLZLH6zeq/JqX/B5/ptLi3Tn22s1eVA3nXVEX502rg+VYAC0CQHdAaG95xI3iQJANEmMi9XJY3rp5DG9VFZVo4/W5Or1Zbv0ybp9qgqZ5GLR1gNatPWAfv/6Kk0d0l2nje+jU8b2Us+0JBdbDyAaEdAdwPhzAOgYUhLi/DeLHt5XheXVen/VHr2+bJfmb9ofLAZQa6X5m/Zr/qb9+u1rKzV5UDedNq6PTh3XW30ZBgOgGQjoDqCCCwB0POnJ8Tp/8gCdP3mA8koq9c7KPXpj6S4tzMkPPsda6cucA/oy54D+8OZqHT4gQ6eN663TxvXWoO6pLrYegJcR0B0QFtDpQQeADierS6IumTpIl0wdpD2FFXpv1R69s3K3Fm7JV+gox2XbC7Rse4H+7521GtOnqz+sj++tYT3T3Gs8AM8hoDuAIS4A0Hn0Tk/SpdOydem0bOWVVOqD1Xv19ordWrBpf1jBgNW7i7R6d5H++sF6DemRqpNH99JJY3pp4sBu/K4AOjkCugPCSyy62BAAgKOyuiQGa6wXlFXpwzW5enflbs3dkBecvVSSNu8r1SP7NuuRuZuVmZqgmSN76uQxvXTM8CylJvKrGuhs+Kl3AENcAAAZKQk6b1J/nTepv4orqvXx2ly9u3KP5qzbF1a6Mb+0Si8t2aGXluxQQlyMpg/trpPG9NJJo3upV1cqwgCdAQHdAeEBPcbFlgAAvCAtKV5nH9FPZx/RTxXVPs3flKcPVufqozV7lVtcGXxeVU2tZq/bp9nr9unWV1bqsP7pOml0L50wqqfG9u0qQ+EBoEMioDsgfAy6iw0BAHhOUnysThjVSyeM6qXa2nFasbNQH6zeqw/X7NXaPcVhz12+o1DLdxTqng/Wq2daoo4b0UMzR/XUjOFZ6poU79JXACDSCOgOoMwiAKA5YmKMDh+QocMHZOjmU0Zqe36ZPlyzVx+s3qsvtuSH/T7JLa7UC4t36IXFOxQbYzRpUDcdP7KHjh/RU6P7pNG7DkQxAroDauvvA2IWUQBAsw3ITNFl0wfrsumDVVhWrTnrc/XRmlzN3bBPBWXVwef5aq0WbsnXwi35uuvdderdNSnQu95D04dlKY3edSCqENAdEDrEJY6ADgBohfSU+nHrvlqrpdsL9Mm6XM1et08rdhaGPXdPUYX+t2i7/rdou+ICvevHDM/SjOE9NL5fOgULAI8joDvAF9KFTg86AKCt6oa0TBrUTT+dNVL7iis1d/0+zV6Xq3kb8lRYXt+7XlNr9cWWfH2xJV9/eX+9uibFadrQLM0YnqVjhmcxoyngQQR0B/hChrgwBh0AEGk90hL1rUn99a1J/VXjq9XS7QWas26f5qzP1cqdRWHPLaqo0bur9ujdVXskSQMykzVjWA/NGJal6cO6KyMlwY0vAUAIAroDqIMOAHBKXGyMJmdnanJ2pm4+ZaRyiyv02cY8zduQp0835IWVcZSk7fnlenbhNj27cJuMkcb3Sw+E9SxNGtRNSfGxLn0lQOdFQHdArQ2dSZSADgBwTs+0JH1zQn99c0J/WWu1IbdEn27I06cb8/T55v0qq6qfJMna+lKOD8/ZpITYGB0xIENTh3bX1CGZmjiQwA44gYDugJqQHvS4WAI6AMAdxhiN6JWmEb3SdPmMwaqqqdVX2w74e9g35mnZ9gKF/MpSla9WC3PytTAnXw98JCXExWjCgAxNHdJdRw/triMGZBDYgXZAQHdA6BAXetABAF6REBejKUO6a8qQ7vrprJEqLK/Wgk379dnGPC3YvF8bc0vCnl9VUxu84fT+jzYoIS5GEwcGAvuQ7jpiYIYS4wjsQFsR0B1QaxmDDgDwvvTkeJ06rrdOHddbkrSvuFJfbNmvBZv26/PN+7VpX2nY86tqavX55nx9vjlf98kf2I/on6HJ2d10ZHamJg7qpvRkarADLUVAdwAziQIAolGPtER947C++sZhfSVJucUVgUDuD+ybGwnsdUNipE0yRhrZKy0Y2CdnZ6pfRrILXwkQXQjoDqCKCwCgI+iZlqSzDu+rsw73B/a9RRWBsO4P7VvywgO7tdLaPcVau6dYT3++TZLUNz1Jk7MzdWR2N03OztSIXmn8bgQaIKA7gIAOAOiIenVNCs5uKvmHxCzemq8vcw5oUU6+Vu4qCvsdKEm7Civ0+rJden3ZLklSWlKcJgzspgkDMjRhYIaOGJBBLXZ0egR0B/hCyywS0AEAHVSPtESdOq6PTh3XR5JUWlmjZdsL/IF9a76WbD2g0pCyjpJUXFGjuev3ae76fcF1Q3qkasKAbjpiYIYmDMjQqN5piouNcfRrAdxEQHdAbdgYdBcbAgCAg1IT4zRtWJamDcuSJNX4arV2T7EWbsnXokBP+74GEydJ0uZ9pdq8r1QvLdkhSUqOj9X4/umaMDBDEwZ008SBGerZNcnRrwVwEgHdAeFDXOgBAAB0TnGxMRrXL13j+qXr8hmDZa3VtvwyLd1eoK+2FeirbQe0aldR2PwhklRe7dPCLflauCU/uK5fRrIOH5Cu8f0ydFj/dI3rm670FCrGoGMgoDsgPKC72BAAADzEGKNB3VM1qHtqcBx7RbVPq3YVBgK7P7TvKqz42mt3FpRrZ0G53l6xJ7guu3uKxvfP0Ph+XTW+X4bG9euqtCRCO6IPAd0BPuqgAwDQLEnxsZo0KFOTBmUG1+0prNDS7QeCoX35zgJVVNd+7bU5+8uUs79MbwRuQJX849kP65eu8f39Pe1j+3ZVSgLxB97GGeoAZhIFAKD1eqcn6dT0+ptPq321WrenWCt2Fmr5jkKt2FmgtbuLvzY0Rqofz/7qUn9ojzHS0B5dNLZvV43p21Vj+6ZrdJ+uykylcgy8g4DugNCZROPoQQcAoE3iQ8ayX3CUf11FtU/r9hRr+c5CrdhRoBU7i7R+b/HXyjzWWmlDbok25JYEQ7sk9UlP0pg+/tBe9zigWwrV1+AKAroDanyUWQQAoD0lxcfq8AEZOnxAhqRBkvyhffXuIq3YUd/TvjG3RI10tGt3YYV2F1boo7W5wXVdEuM0uk9aSHBP1/BeXZQUH+vMF4VOi4DugNAe9FiGuAAA4Iik+FhNHNhNEwd2C64rrazRmt1FWrO7SKt3F2n1riKt3VOsypqvj2kvqazRlzkH9GXOgeC62BijIVmpGtE7TaN6pWlE7zSN7JWmgZn0tiNyCOgO8IX8zHOTKAAA7klNjNPk7ExNzq6/CbXGV6steaVavbtIq3b5Q/vq3UXKL6362ut9tTY4ROYt7Q6uT46P1fBeXTSyV5pG9k7TiF5pGtU7TT3SEmXonEMLEdAdwEyiAAB4V1xsjIb3StPwXmnBco/WWu0tqtTq3YXBwL56V5Fy9pc1uo/yap+WB4bShMpIiQ+G9RGB8D68ZxdlpHBTKppGQHeAL6QLnZtEAQDwPmOMeqcnqXd6kk4Y1Su4vrSyRuv3Fmv93mKt21OidXuLtG5PifJKvj4jqiQVlFV/bZIlScrqkqAhPbpoWM8uGlb32LOL+qQn0eMOAroTQu4RpcwiAABRLDUxThMGdtOEkHHtkrS/pFLr9hZr/Z5irdtbonV7irR+b4lKKmsa3U9eSZXySr4e3FMSYjU0JLD7P/dP5hTPbIedBgHdAbW1TFQEAEBH1r1LoqZ1SdS0oVnBddZa7Sqs0Lo9/l52f697sTbnlTQ60ZIklVX5tGJnoVbsDB8qExdjNKh7iob17KIhPbpocFaqBmelKrt7qrK6JNDr3sEQ0B3ATKIAAHQ+xhj1y0hWv4zksGEytbVWOwvKtXFfiTbllmhj3ce+EhWUVTe6r5paq037SrVpX6mkvWHb0hLjlF0X2LNSNSTwOLh7qtJT4tvzS0Q7IaA7wEcPOgAACIiJMRqQmaIBmSmaObJncL21VvmlVcGwXhfcN+8r1c6C8ib3V1xZ02ivuyRlpiYou3uKBmd10eAs/2N2Voqyu6cqNZEY6FV8ZxwQFtD5FxQAAGiEMUbduySqe5dETRnSPWxbaWWNNu8r1cZ9xdqyr1Rb9pdpS16JcvLKmhznLkn5pVXKL63Skm0FX9uW1SVBAzJTNLDhR/cU9UpLovKciwjoDggN6JzsAACgpVIT4zS+f7rG908PW2+t1b6SSuXk+QP7lrz64J6zv7TRCZjq+G9UrdJXjYT3hLgYDeiWHAztdUF+UPdUDchMVkoCEbI9cXQdwEyiAACgPRhj1DMtST3TknTU4MywbbW1VruLKpSTV6rNeaXKySvVlsDHjgNlqg4tM9dAVU1tyJj3r8vqkqiBmckakJniH2ffzT/Wvn+3ZPXLSFFyQmxEv87OhoDugNAe9LhYAjoAAGh/MTH1N6lOH5YVts1Xa7WnqELb9pdpW36ptuWXaVt+ubbll2l7flmjs6iGyiupVF5JZaNDZyT/2Pf6wF4f4Pt1S1b/jBR1TY6j8sxBENAdEDbEhZMRAAC4LDYkvB89tPvXthdXVAfD+rbAx9b9/uUdB8pVU9t077tUP/a9sRtXJalLYlx9YA+E9z4ZyeobmByqZ1qSEuI6b913AroDwqu4uNgQAACAZkhLitfYvuka2zf9a9t8tVa7C8u1bX+ZdhSUa+eBcu0sKNeOA2XaWVCu3QUVhwzwJZU1Wre3WOv2Fje63Rj/MJo+6Unq3TVJfTOS1Ts9KbjcJz1ZvdITlRjXMYfSENAdEFoHnR50AAAQzWJjjPp3S1H/bimNbvfVWuUWV4QE95DHQIhvaqKmOtZK+4orta+4UsvVeC+85K9E0zs9Sb27JvvDe12IDwT53ulJUXlDa/S1OAoxkygAAOgsYmOM+qQnq096siY3sr2u3vvOsN73cu0prNDuwnLtLqzQvpJK2YN3wkuqr0SzcmdRk89JS4pTr65JevzSIzWwe+N/VHgNAd0Bof/miSOgAwCATiy03vth/TMafU61r1a5xZXaEwjsuwsqtLuwQnuK/Mt7Ciu0t6hChxhJI0kqrqhRcUWJuiRFT+yNnpZGsdAyi9RBBwAAOLj42JjgTaxNqfHVal9JZTCw+x8DgT4Q4HOLKlXlq1VCbIy6pcQ7+BW0DQHdAcwkCgAAEFlxsTHBoTRNsdbqQFm19pdURlVZRwK6A3wh90HQgw4AAOAMY4wyUxOUmZrgdlNahKJ/DvDV1id0xqADAADgYAjoDgidSZcqLgAAADgYAroDaplJFAAAAM1EQHeAjzroAAAAaCZuEnXAVccN0VlH9JWv1uqIARluNwcAAAAeRkB3wMSB3dxuAgAAAKIEQ1wAAAAADyGgAwAAAB5CQAcAAAA8hIAOAAAAeAgBHQAAAPCQiAV0Y0x/Y8zjxphdxphKY0yOMeY+Y0yLSpgYYzIDr8sJ7GdXYL/9I9VWAAAAwKsiUmbRGDNU0nxJPSW9JmmtpKMk3SjpVGPMdGvt/mbsp3tgPyMkfSzpOUmjJF0m6QxjzNHW2s2RaDMAAADgRZHqQX9Y/nB+g7X2HGvtL6y1J0i6V9JISXc0cz93yh/O77HWnhjYzznyB/2egfcBAAAAOqw2B/RA7/ksSTmSHmqw+XeSSiVdYoxJPcR+uki6JPD83zfY/KCkrZJOMcYMaWubAQAAAK+KRA/6zMDj+9ba2tAN1tpiSZ9JSpE09RD7mSopWdJngdeF7qdW0nsN3g8AAADocCIxBn1k4HF9E9s3yN/DPkLSR23cjwL7OShjzOImNo061GsBAAAAN0WiBz098FjYxPa69RkO7QcAAACIWhGp4uI11tpJja0P9KxPdLg5AAAAQLNFoge9rmc7vYntdesLHNoPAAAAELUiEdDXBR6bGhs+PPDY1NjySO8HAAAAiFqRCOizA4+zjDFh+zPGpEmaLqlM0ueH2M/nksolTQ+8LnQ/MfLfaBr6fgAAAECH0+Yx6NbaTcaY9+UP0NdJ+lvI5tskpUp6xFpbWrfSGDMq8Nq1IfspMcb8R9JV8tdBvylkPz+SlC3pvTbOJJq9Zs0aTZrU6BB1AAAAICLWrFkj+fNrixlrbZsbEJisaL78s32+JmmNpCny1yxfL2matXZ/yPOtJFlrTYP9dA/sZ4SkjyUtlDRa0tmScgP72dSGdm6R1FX+SZWcVlfice1Bn4U6HK+W4Xi1DMerZTheLcPxajmOWctwvFrGreOVLanIWju4pS+MSECXJGPMAEl/kHSqpO6Sdkt6RdJt1toDDZ7baEAPbMuUfwbScyT1kbRf0juSfmut3RGRxrqgrjZ7UxVmEI7j1TIcr5bheLUMx6tlOF4txzFrGY5Xy0Tj8YpYmUVr7XZJlzXzuV8L5iHb8iXdGPgAAAAAOpVI3CQKAAAAIEII6AAAAICHENABAAAADyGgAwAAAB4SsSouAAAAANqOHnQAAADAQwjoAAAAgIcQ0AEAAAAPIaADAAAAHkJABwAAADyEgA4AAAB4CAEdAAAA8BAC+iEYY/obYx43xuwyxlQaY3KMMfcZY7q1YB9zjDH2IB9JTbxujDHmeWNMrjGmwhizzhhzmzEmOXJfYWS19XgZY44/xLGq+xjQ4HUHe+7n7fPVto0x5jxjzN+MMfOMMUWBtj7dyn21+LhH2/kVieNljOlujPmBMeYVY8xGY0y5MabQGPOpMeYKY8zXronGmOxDnF/PRe6rjJxInV+Bc6mpr33PQV43zRjztjEmP3CclxtjfmyMiW3bV9Y+InR+fb8Z1y5fg9dE6/nV4p+lQ+yvQ1/DInW8Oss1LJLnV7Rew+KceJNoZYwZKmm+pJ6SXpO0VtJRkm6UdKoxZrq1dn8LdnlbE+trGnnvKZI+lhQv6UVJ2yWdIOm3kk40xpxora1swXu3uwgdrxw1fZzGSzpX0kpr7fZGtm+V9EQj63ccsvHu+LWkwyWVyN/GUa3ZSWuOezSeX4rM8Tpf0t8l7ZY0W9I2Sb3kP6/+Jek0Y8z5tvEZ3JZJerWR9Stb0Q4nROT8CiiUdF8j60sae7Ix5mxJL0mqkPQ/SfmSzpR0r6Tp8n8fvCYSx2upmr5+HSP/z9g7TWyPtvOrLT9LYTrJNSxSx6uzXMMidn4FRN81zFrLRxMfkt6TZCVd32D9PYH1/2jmfub4D3Wz3zdW0urAe5wVsj5G/guRlfQLt49Pex2vg+z/2cB+bmhkm5U0x+1j0MKvZ6ak4ZKMpOMDX8PT7X3co/j8avPxkv8X+JmSYhqs7y3/LwAr6VsNtmUH1j/h9jFw6fzKkZTTgud3lZQrqVLS5JD1SfKHMCvpu24fn/Y6XgfZ/4KGP3NRfn61+GfpIPvq8NewSB2vznINi/D5FZXXMNe/CV79kDQ08E3Y0sgJkib/X12lklKbsa85allAPyHw3p80sm1IYFuOJOP2cWqP49XE/rPk/0u2TFJGI9ujLqA3aP/xal3gbPFxj8bzK1LH6xD7/FVgn39rsD7qfrlF8ni14pfb5YH3erKRbU2ee176iPT5Jf9//6z8PfOxDbZF/fnVyNfb6M9SE8/tlNew1h6v1uyno51jLT1e0XoNY4hL02YGHt+31taGbrDWFhtjPpM0S9JUSR81Z4fGmO9IGiypStIaSR/bxv8Fd0Lg8d2GG6y1m40x6yWNkP9CtKk57+2AiB+vBi6VlCjpKWttQRPPyTDGXC7/X9iFkhZbaz05/jyCWnPco/H8ckJ14PFrQ84C+hpjfiipu6T9khZYa5c70jL3JRpjLpY0UP6wtFzSXGutr5HnNnl+SZor/x/Z04wxiU1c/zqiqwKPjzVxzKSOdX4d6mcpFNewlh2vtuyno5xjrTleUXcNI6A3bWTgcX0T2zfIf9EYoeYHzoY3YuQaY66z1r7YivceEfjwysWnPY5XqCsDj48c5DmHS3osdIUxZpmkS6y1K1rxntGgNcc9Gs+vdmWMiZP0vcBiYxdlSTo58BH6ujmSLrXWbmu/1nlCb0n/abBuizHmMmvtJw3WN3l+WWtrjDFbJI2VPzytiXhLPSZww+LFknzyj51tSoc4v5r5sxSqU1/DWnG82rKfqD/H2nC8ou4aRhWXpqUHHgub2F63PqMZ+3pN/rFU/SUly3/z0Z8Cr/2fMebUdnxvp7Rbm40xx8n/A7PSWju/iafdI/+NGz3k/7fokfKPRTxc0sfGmH4tfd8o0ZrjHo3nV3v7P0njJL1trX2vwbYySbdLmiSpW+DjOPlvXDpe0kfGmFTnmuq4f0s6Uf5fcKnyD9d4RP5/m79jjDm8wfM5v8J9W/6v9V3b+M3tHe38OtjPUmM6+zWspcerNfvpSOdYa45XVF7DCOgOsNbea61901q701pbYa1dZ639laSb5P8e/MnlJnpd3b+H/9nUE6y1N1lr51tr86y1JdbaRdba8+W/CztL0s1ONBTRxxhzg/w/i2slXdJwu7U211r7W2vtEmttQeBjrvy9el9IGibpB4422kHW2tustR9ba/daa8ustSuttVfL/0dxsqTfu9tCz6u7fjX637+OdH4d6mcJ4SJ1vDrLNay1xytar2EE9KbV/YWU3sT2uvUFbXiPf8k/huoIY0yaw+8dae3SZmNMpqRvSSrX1/891Rz/CDwe24rXRoPWHPdoPL/ahTHmR5Lul78ixExrbX5zX2utrVH9kIWOen4dTFM/W5xfAcaYsZKmyX9z6NsteW20nV9t+FnqlNewtlx7IrWfaDrHInW8GvD0NYyA3rR1gccRTWwfHnhsagzcIVlrKyQVBxZD/73U7u/dDtqrzXU3hz5/kJtDD2Zf4DFa/n3XUq057tF4fkWcMebHkv4mfw3gmdbaJiesOIiOfn4dTFNfe5PnV2D86GD5OyY2t1/TPKM5N4ceTFScX238Wep017AIXXs6zTUsUserEZ6+hhHQmzY78Dir4YxVgd7u6fKP62p1lRBjzEj5x4IVS8oL2fRx4LHh2HQZY4bIf9Jslbd+wbXX8aq7ObTJ4S2HMDXw6KVjFUmtOe7ReH5FlDHm5/JPOLFU/gt+bit31dHPr4Np6mtv8vySv6cqRdL8jl7BxfhniL5E/ptDHzvE05vi+fMrAj9LneoaFqlrT2e5hkXw62yMt69h7VG7saN8qOWTJ4ySNKrBusGSMhvZdw/VF7z/Z4NtB5uE4QV5cBKGSB2vBtuPCbxuxSHe9zBJ8U2szwvs40K3j88hvobjdZC6y/LPljdK0tAIHPeoPL8ieLx+E3jtosZ+Nht5/kQ1qM8cWH+i/LX5raRpbh+T9jhekkarkbkL5L+5akNgn79qsK2r/D1TUTVRUaTOr5DnXBLYxxsd9fxqyc8S17CIHq9OcQ2LxPGK5muYCbwpGtHI9MNrJE2Rv27revlP6P0hz7eSZK01Ieu+L/84p0/l/ystX/46nKfLP45pkaSTbYPhG41MY7xN/h+myZI+k+S1aYwjcrwa7O8/8pcnu8Fa+7eDvO8T8lfJmSf/dM+V8v+gnir/hfxRST+0HjvZjTHnSDonsNhb0inynyPzAuvyrLU3B56bLf9EHluttdkN9tOi4x54TTSeX+eojcfLGHOppCfk79X8mxq/Sz/HWvtEyGvmyP8v8/nyjyWW/H/81dXK/Y219o+t/8raR4SO1+/lvylrrvw9ksXyTyxzhvy/rN6W9E1rbVUj7/2i/L/8n5P/uneW/NWYXpT07Y768xiyv3mSZsgfIN84yPvOUXSeXy36Wers17BIHa/Ocg2L4PH6vaL1Gub2X0he/5A0QP4SPbvln2Boq6T7JHVr5LlWDWYMlb+czxOSVsg/MUB14Bs9T9L1khIO8t5j5O8NyJM/dK6XdJukZLePS3sdr5Bt3eS/MbTRmUMbPPccSS9L2iipKPC+uyW9oQbTanvpQ/47x+1BPnJCnpvdcF1rj3u0nl+ROF7N2IdVgxlpJV0h6U35Z6MrCRyrbZL+J+kYt49LOx+v4yQ9K3/VhAL5r1/7JH0gfy3iJmdqlH9owtuSDgR+lldI+okazKTplY8I/zyODmzffqivtwOfX2E/S804Zh36Ghap49XS/UTrORbB4xW11zB60AEAAAAP4SZRAAAAwEMI6AAAAICHENABAAAADyGgAwAAAB5CQAcAAAA8hIAOAAAAeAgBHQAAAPAQAjoAAADgIQR0AAAAwEMI6AAAAICHENABAAAADyGgAwAAAB5CQAcAAAA8hIAOAAAAeAgBHQAAAPAQAjoAAADgIQR0AAAAwEP+H9QUXF44p7axAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "image/png": {
       "height": 248,
       "width": 372
      },
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import nixio\n",
    "\n",
    "%config InlineBackend.figure_formats = ['retina'] # only for users with a high resolution display\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "nixfile = nixio.File.open(\"tagging1.nix\", nixio.FileMode.ReadOnly)\n",
    "block = nixfile.blocks[0]\n",
    "\n",
    "stimulus_tag = block.tags[\"stimulus\"]\n",
    "stimulus_response_data = stimulus_tag.tagged_data(\"response\")[:]\n",
    "time = stimulus_tag.references[\"response\"].dimensions[0].axis(len(stimulus_response_data), start_position=stimulus_tag.position[0])\n",
    "\n",
    "plt.plot(time, stimulus_response_data)\n",
    "nixfile.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Tags** work in the unit of the data, respectively its dimension. It is, however, possible to provide diverging units to the **Tag**.\n",
    "\n",
    "```\n",
    "    # *********  This is the interesting part **********\n",
    "    data = block.create_data_array(\"response\", \"nix.sampled\", data=response, label=\"voltage\", unit=\"mV\")\n",
    "    data.append_sampled_dimension(interval, label=\"time\", unit=\"s\")  # dimension is defined in seconds\n",
    "\n",
    "    stim_tag = block.create_tag(\"stimulus\", \"nix.stimulus_segment\", position=[stim_on * 1000])\n",
    "    stim_tag.extent = [(stim_off - stim_on) * 1000]\n",
    "    stim_tag.units = [\"ms\"]  # the tag is specified in milli seconds\n",
    "    stim_tag.references.append(data)\n",
    "    # **************************************************\n",
    "```\n",
    "\n",
    "Retrieving the tagged data works regardless. Units are scaled and handled transparently by the library. \n",
    "\n",
    "**Note:** This works only for SI units."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Tagging in 2-D (and n-D)\n",
    "\n",
    "The same mechanism extends to n-dimensions. In this case the **Tag's** position and extent have to be lists of the length n.\n",
    "\n",
    "![2dtag](resources/2d_tag.png)\n",
    "\n",
    "With this we can create a \"Region of Interest\" (ROI) for example in image data. Image data is usually 3-D with the dimensions representing height, width and the three (sometimes four) color channels\n",
    "\n",
    "![roi](resources/single_roi.png)\n",
    "\n",
    ">  \"Lenna\" by Original full portrait: \"Playmate of the Month\". Playboy\n",
    ">  Magazine. November 1972, photographed by Dwight Hooker.This 512x512\n",
    ">  electronic/mechanical scan of a section of the full portrait:\n",
    ">  Alexander Sawchuk and two others[1] - The USC-SIPI image\n",
    ">  database. Via Wikipedia -\n",
    ">  http://en.wikipedia.org/wiki/File:Lenna.png#mediaviewer/File:Lenna.png"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "from PIL import Image as img\n",
    "\n",
    "\n",
    "def load_image():\n",
    "    image = img.open(os.path.join(\"resources\", \"lenna.png\"))\n",
    "    pix = np.array(image)\n",
    "    channels = list(image.mode)\n",
    "    return pix, channels\n",
    "\n",
    "\n",
    "img_data, channels = load_image()\n",
    "\n",
    "# create a new file overwriting any existing content\n",
    "file_name = 'single_roi.nix'\n",
    "nixfile = nixio.File.open(file_name, nixio.FileMode.Overwrite)\n",
    "block = nixfile.create_block(\"session 1\", \"nix.session\")\n",
    "\n",
    "data = block.create_data_array(\"lenna\", \"nix.image.rgb\", data=img_data)\n",
    "# add descriptors for width, height and channels\n",
    "data.append_sampled_dimension(1, label=\"height\")\n",
    "data.append_sampled_dimension(1, label=\"width\")\n",
    "data.append_set_dimension(labels=channels)\n",
    "\n",
    "# create a Tag, position and extent must be 3-D since the data is 3-D\n",
    "position = [250, 250, 0]\n",
    "extent = [30, 100, 3]\n",
    "tag = block.create_tag('Region of interest', 'nix.roi', position)\n",
    "tag.extent = extent\n",
    "tag.references.append(data)\n",
    "\n",
    "nixfile.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Tagging multiple points or segments\n",
    "\n",
    "Sometimes we do not only want to tag single points, segments or regions but sets of such. In the following situation a signal has been recorded and within this signal certain events, threshold crossings have been detected (figure below).\n",
    "\n",
    "![multiple points](resources/multiple_points.png)\n",
    "\n",
    "For storing this kind of data we need two **DataArrays**, the first stores the recorded signal, the other the events. Finally, a **MultiTag** entity is used to link both. One can use the event times stored in one of the DataArrays to tag multiple points in the other signal.\n",
    "\n",
    "``` python\n",
    "    signal_array = block.create_data_array(\"signal\", \"nix.sampled\", data=signal, label=\"voltage\", unit=\"mV\")\n",
    "    signal_array.append_sampled_dimension(sampling_interval, label=\"time\", unit=\"s\")\n",
    "\n",
    "    event_array = block.create_data_array(\"threshold crossings\", \"nix.events.threshold_crossings\", data=events, label=\"time\", unit=\"s\")\n",
    "    event_array.append_range_dimension_using_self()\n",
    "\n",
    "    mtag = block.create_multi_tag(\"event tag\", \"nix.tag.events\", event_array)\n",
    "    mtag.references.append(signal_array)\n",
    "```\n",
    "\n",
    "## Tagging multiple segments\n",
    "\n",
    "A very similar approach is taken for tagging multiple segments in which, for example, a stimulus was switched on.\n",
    "\n",
    "![multiple regions](resources/multiple_regions.png)\n",
    "\n",
    "For storing such data we again need one **DataArray** to store the recorded signal. Storing the regions is similar to the approach for the simpler Tag, i.e. positions and the extents need to be provided. Accordingly, two additional **DataArrays** are required. The first of which stores the positions and the second the extents of the tagged regions.\n",
    "\n",
    "```python\n",
    "    block = nixfile.create_block(\"multiple regions\", \"nix.session\")\n",
    "\n",
    "    data_array = block.create_data_array(\"signal\", \"nix.data.sampled\", data=signal)\n",
    "    data_array.label = \"voltage\"\n",
    "    data_array.unit = \"mV\"\n",
    "    data_array.append_sampled_dimension(sampling_interval, label=\"time\", unit=\"s\")\n",
    "\n",
    "    positions = block.create_data_array(\"stimulus onsets\", \"nix.region.onsets\", data=stim_onsets)\n",
    "    positions.append_set_dimension()\n",
    "\n",
    "    extents = block.create_data_array(\"stimulus extents\", \"nix.region.extents\", data=stim_extents)\n",
    "    extents.append_set_dimension()\n",
    "\n",
    "    mtag = block.create_multi_tag(\"stimulus segments\", \"nix.segments.stimulus\", positions=positions, extents=extents)\n",
    "    mtag.references.append(data_array)\n",
    "```\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Tagging in n-d\n",
    "\n",
    "The same principle as demonstrated above applies also to n-dimensional data. Tagging in n dimensions requires positions and extents stored in DataArrays of appropriate shapes.\n",
    "\n",
    "The following figures show the tagging of multiple regions in 2- and 3-D.\n",
    "\n",
    "![tagging 2d](resources/2d_mtag.png)\n",
    "\n",
    "Tagging multiple regions in n-D data requires the DataArrays for storing positions and extents to be two-dimensional. The first dimension represents the number of regions, the second has as many entries as the referenced data has dimensions.\n",
    "\n",
    "According to the number of dimensions of the data (here, width and height) each starting point and the extent of a tagged region is defined by two numbers. Thus, the position and extent DataArrays are two dimensional. The first dimension represents the number of tagged regions, the second the number of dimensions.\n",
    "\n",
    "This approach can be extended into n-D. The following figure illustrates the 3-D case.\n",
    "\n",
    "![tagging 3d](resources/3d_mtag.png)\n",
    "\n",
    " Again, the **DataArrays** for positions and extents are always 2-D, the first dimension represents the number of tagged regions, the second the number of dimensions.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Time for a break and some exercises"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "aee8b7b246df8f9039afb4144a1f6fd8d2ca17a180786b69acc140d282b71a49"
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}