Standardized_PGAM_Fit_with_examples.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "bfc3d69b",
   "metadata": {},
   "source": [
    "# Standardizing PGAM fits\n",
    "--------------------------------------\n",
    "\n",
    "In this notebook we propose a pipeline for fitting PGAMs (Balzani et al., 2020, NeurIPS) with a standard input and output format. The procedure has the following advantages:\n",
    "\n",
    "1. It requires minimal coding (limited to formatting the input data in a standard format)\n",
    "2. The output format is easily exportable to other tools and programing languages (e.g., MATLAB)\n",
    "3. The pipeline is compatible with the PGAM <a href=\"https://hub.docker.com/r/edoardobalzani87/pgam\"> Docker</a> and <a href=\"https://osf.io/pcsav/\"> Singularity</a> images.\n",
    "4. The pipleine can be easily parallelized for HPC usage (singularity containers).\n",
    "\n",
    "\n",
    "## Table of contents\n",
    "* [Standard input format](#standard-input-format)\n",
    "    * [Neural & behavioral inputs](#inputs)\n",
    "    * [Configuration files](#configuration-files)\n",
    "* [Example with synthetic data](#example-with-synthetic-data)\n",
    "    * [Create and save an example configuration file](#create-and-save-conf)\n",
    "    * [Fittng the model](#fitting-model)\n",
    "        * [List the fits by experiment, condition, neuron, and model configuration](#fit-list)\n",
    "        * [Load inputs, fit and postprocess](#load-fit-save)\n",
    "* [Example with real data](#example-with-real-data)\n",
    "    * [Create and save an example configuration file](#create-and-save-conf-for-real-data)\n",
    "    * [Fittng the model](#fitting-model-for-real-data)\n",
    "        * [List the fits by experiment, condition, neuron, and model configuration](#fit-list-real-data)\n",
    "        * [Load inputs, fit and postprocess](#load-fit-save-real-data)\n",
    "* [Fit via docker](#docker)\n",
    "    * [Mounting volumes](#mount-volumes)\n",
    "    * [Setting paths](#container-path)\n",
    "    * [Fit by running the container](#fit-container)\n",
    "* [HPC & singularity ](#HPC)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e9aa0044",
   "metadata": {},
   "source": [
    "\n",
    "\n",
    "## Standard input format <a name=\"standard-input-format\"></a>\n",
    "\n",
    "In order to fully specify the PGAM model we need:\n",
    "\n",
    "1. A matrix containing the population <a href=\"#counts\" >spike counts</a> (to be used as response variable and/or covariates), and a list with the <a href=\"#neu_names\">neuron identifiers</a>\n",
    "\n",
    "2. A matrix containing the <a href=\"#variables\" >task variables</a> (additional covariates), and a list with the respective <a href=\"#variable_names\" >identifiers</a>\n",
    "\n",
    "3. A vector containing <a href=\"#trial_id\" >trial IDs</a>\n",
    "\n",
    "4. A list of <a href=\"#trial_id\" >covariates</a> to be included and the ID of the unit that we want to fit.\n",
    "\n",
    "5. A <a href=\"#configuration-files\" >configuration file</a> listing the model parameters in \"YAML\" file format\n",
    "    \n",
    "We propose the following standard input format to facilitate model specification, fitting, and saving. \n",
    "\n",
    "### **Neural & behavioral inputs [1-4].** <a name=\"inputs\"></a>\n",
    "\n",
    "Inputs [1-4] will be contained in a single \".npz\" files with keys: \n",
    "\n",
    "- **counts**: <a name=\"counts\"></a>numpy.array of dimension (T x N), where T is the total number of time points (bins), and N is the number of neurons. The code has been extensively tested with bins of size 6ms, but works for any size.\n",
    "- **variables**: <a name=\"variables\"></a>numpy.array of dimension (T x M), T as above (same bin size), M the number of task variables. For continous variables, each entry is the value taken by the associated variable. For event variables, it is \"1-hot-encoding\": zeros everywhere except for when the event happens, taking a value of '1'. \n",
    "- **trial_id**: <a name=\"trial_id\"></a>numpy.array of dimension T, trial ids associated to each time point. All time points for a given trial take the same value. \n",
    "- **variable_names**: <a name=\"variable_names\"></a>numpy.array of strings of dimension M, name of each covariate in 'variables'.\n",
    "- **neu_names**: <a name=\"neu_names\"></a>numpy.array of strings of dimension N, label uniquely identifying the neurons in 'counts'\n",
    "- **neu_info**: dict, *optional\n",
    "    neu_info[neu_names[k]]: <a name=\"neu_info\"></a>dict, info about the k-th neuron, $k=0,\\dots,N-1$. keys are the information label.\n",
    "\n",
    "\n",
    "### **Configuration files [5]**<a name=\"configuration-files\"></a>\n",
    "\n",
    "Input [5], the configuration parameters for the B-spline, will be stored with the \"YAML\" file format. Python dict objects can be readily saved in YAML format as follows,\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "```python\n",
    "    import yaml\n",
    "    with open('data.yml', 'w') as outfile:\n",
    "        yaml.dump(python_dictionary, outfile, default_flow_style=False)\n",
    "```\n",
    "\n",
    "where *python_dictionary* is any python dict object whose values are strings, numeric, list containig strings or numeric, or dict.\n",
    "\n",
    "To read a \"YAML\"\n",
    "\n",
    "```python\n",
    "    with open(\"data.yaml\", \"r\") as stream:\n",
    "        try:\n",
    "            python_dicttionary = yaml.safe_load(stream)\n",
    "        except yaml.YAMLError as exc:\n",
    "            print(exc)\n",
    "```\n",
    "\n",
    "The dictionary with the \"YAML\" should be structured as follows:\n",
    "```yaml\n",
    "    var_1: # This is an example of a 2D continous variable. \n",
    "      is_temporal_kernel: false # false for continous variable, true for events (example below)\n",
    "      kernel_direction: .nan # '.nan' when \"is_temporal_kernel: false\". When \"true\", 0 = acausal (bidirectional), '1' = causal (i.e., firing change after the event happens), '-1' = anticipatory (i.e., firing change before event happens).\n",
    "      kernel_length: .nan\n",
    "      order: 4 # Order of the B-spline. A step-function is order = 1. \n",
    "      knots_num: .nan # \".nan\" will use the knots passed (below). An integer will place knots equi-spaced.\n",
    "      knots: # list of knot lists for tuning function. Provide 1 set of knots for each dimesion of the basis. The first and last entries are repeated order-times (given recursive nature of B-spline definition).\n",
    "      - # first set of knots\n",
    "        - -5.0\n",
    "        - -5.0\n",
    "        - -5.0\n",
    "        - -5.0\n",
    "        - -3.888888888888889\n",
    "        - -2.7777777777777777\n",
    "        - -1.6666666666666665\n",
    "        - -0.5555555555555554\n",
    "        - 0.5555555555555554\n",
    "        - 1.666666666666667\n",
    "        - 2.7777777777777786\n",
    "        - 3.8888888888888893\n",
    "        - 5.0\n",
    "        - 5.0\n",
    "        - 5.0\n",
    "        - 5.0\n",
    "      - # second set of knots\n",
    "        - -5.0\n",
    "        - -5.0\n",
    "        - -5.0\n",
    "        - -5.0\n",
    "        - -3.888888888888889\n",
    "        - -2.7777777777777777\n",
    "        - -1.6666666666666665\n",
    "        - -0.5555555555555554\n",
    "        - 0.5555555555555554\n",
    "        - 1.666666666666667\n",
    "        - 2.7777777777777786\n",
    "        - 3.8888888888888893\n",
    "        - 5.0\n",
    "        - 5.0\n",
    "        - 5.0\n",
    "        - 5.0\n",
    "      penalty_type: der # either 'der' or 'diff'. 'der' is derivate-based penality, 'diff' is difference based penality (see Tutorial 1:PGAM tutorial.ipynb) \n",
    "      der: 2 # order of the derivative used for penalization ('2' means second derivative, penalized wiggliness). Default is '2'.\n",
    "      is_cyclic: # note: should be a vector, one element for each dimension.\n",
    "      - false\n",
    "      - false\n",
    "      lam: 10 # initial lambda for regularization (i.e., how strong is the smoothing penalization. The larger, the more smoothing you get. This is just for initial value, hyper-parameter that is learned in fitting).\n",
    "      samp_period: 0.006 #size of bin for spike counts and behavioral variables (in seconds).\n",
    "      join_var: ['varname_1', 'varname_2'] # list of the variables that should be combined in a multidim response function. \"varname_i\" should match an entry in \"variable_names\". \n",
    "    var_2: # This is an example of an event variable (temporal kernel). Name of \"var_2\" has to be in \"variable_names\". \n",
    "      is_temporal_kernel: true\n",
    "      kernel_direction: 0 # acausal filter (meaning impact of event on firing could happen earlier or after event)\n",
    "      kernel_length: 201 # number of time-points of the total kernel (in this example 201 samples x 6ms). Suggested to be odd number of samples.\n",
    "      order: 4\n",
    "      knots_num: 8 # will create equi-spaces knots.\n",
    "      knots: .nan # '.nan' when 'knots_num' is an integer.\n",
    "      penalty_type: der\n",
    "      der: 2\n",
    "      is_cyclic:\n",
    "      - false\n",
    "      lam: 10\n",
    "      samp_period: 0.006\n",
    "    neuron_X: # Where \"neuron X\" is an input in \"neu_names\"\n",
    "      is_temporal_kernel: true\n",
    "      kernel_direction: 1 # must be '1' for auto-correlation.\n",
    "      kernel_length: 201\n",
    "      order: 4\n",
    "      knots_num: 8\n",
    "      knots: .nan\n",
    "      penalty_type: der      \n",
    "      der: 2  \n",
    "      is_cyclic:\n",
    "      - false\n",
    "      lam: 10\n",
    "      samp_period: 0.006\n",
    "      join_var: .nan # not needed for temporal variables (assumed 1 dimensional). \n",
    "    \n",
    "```\n",
    "\n",
    "where ***var_1*** is a prototypical continous (spatial) variable, ***var_2*** is a prototypical event (temporal) variable. ***spike_history*** will be used for the spike-count auto-correlation term (e.g. the neural spike history as its own predictor, as in an auto-regressive model) **and for the neural couplings (e.g. the counts of simultaneously recoded neurons as predictors)**. See the \"PGAM tutorial.ipynb\" for a definition of continous (spatial) and event (temporal) variables, as well as for a description of the B-spline parameters.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6b5ecc40",
   "metadata": {},
   "source": [
    "# Example with synthetic data <a name=\"example-with-synthetic-data\"></a>\n",
    "\n",
    "Before diving into the pipeline implementation, we will create an example synthetic dataset below.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "290dacda",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfEAAADgCAYAAAATz1/3AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAACAJ0lEQVR4nO29ebwlRXk+/rzd59xz7519mBWGnWGXRUdwFwUEXCLuGhdMzM8Y9yRqxCRGY1Di15hIXHEDVyRqFBVUgrIrMKwybDPDwDDMMPvM3c/S/f7+6K7uqurqPn3We8699Xw+M/ecPt1V1d1V9db7vs/7FjEzLCwsLCwsLPoPznQ3wMLCwsLCwqI5WCFuYWFhYWHRp7BC3MLCwsLCok9hhbiFhYWFhUWfwgpxCwsLCwuLPoUV4hYWFhYWFn0KK8S7DCI6hojuJqJRIno/EX2ViP454/znE9HD3Wxjp0BEbycij4jGiOi4NpR3dFiWR0R/1Y42WljUAxGdQURbMn7PHNOdQji+bm5zmYeEY8wNv19vx1pvwQrx7uMjAK5n5nnMfAkzv4uZP5V2MjPfxMzHdLF9qag3eeXEH5h5LjM/GJb5CSL6nqEuJqKjsgpi5keYeS6Am1psk0Ufg4jeSEQPEtE4EW0koudPZ3vqjel+AjNvDserN91taRVE9BgRnTXd7Wg3rBDvPg4FsC7PiURU6HBb5LqIiGx/sOgrENHZAP4dwF8AmAfgBQAendZGWVh0EXbS7iKI6HcAXgTgi6GJ6mgiuoyI/i38/Qwi2kJE/0BETwH4tq79hqvJDxHRfUS0n4h+RESD0u8fIaJtRLSViP4qS6MNTWMXEdEtACYAHEFEfxFqNaNE9CgR/XV47hwA1wA4MGz7GBEdSEQOEX001IB2E9GVRLS4zc9tn1TneHhPh7WzDou+xScB/Csz/5GZfWZ+kpmfBJTx9BEi2hGOi/OJ6KVE9AgR7SGij4mCiOg0IlpLRCNEtJ2IPp9VMRF9jIh2hWPyzdJx05j+e6kNf6Gd+yUi+lU45m4joiOl348lomvDtj5MRK+XfjuAiK4K23s7gOg6Q1sHieh74RjdR0R3ENHy8LfriegzRHR7OKf8XIxhIjosHG8JhYKIVobz0IfC788iolvD8u8lojMy2nMwEf2UiHaGbfpieNwhon8iosfD5/UdIlogP0utnEi7psCqd2V4zSgRrSOiNeFv3wVwCIBfhPPIR7KeST/BCvEugplfjMD0+97QRPWI4bQVABYj0NjfmVLU6wGcC+BwACcBeDsAENG5AP4OwFkAjgLwwhzNemtYzzwAjwPYAeDlAOYj0G7+k4iezszjAM4DsDVs+1xm3grg/QDOD+s6EMBeAF/KUW9uMPNCUSeALyB4hk+2sw6L/gMFfto1AJYS0YZQWH6RiIak01YAGARwEICPA/g6gLcAeAaA5wP4OBEdEZ77BQBfYOb5CATilRnVrwCwJCz3AgCXElGa22sFgAXhue8A8CUiWiT9/iYEi5FFADYAuCi8vzkArgXwAwDLwvO+TEQnhNd9CcAUgJUA/jL8l4YLwjYcDOAAAO8CMCn9/rbw+gMB1ABcklEWwkX0DQC+yMyfI6KDAPwKwL8hmL8+BOAnRLTUcK0L4JcI5pvDwudyRfjz28N/LwJwBIC5AL6Y1RYNfxaWtRDAVeJaZn4rgM0AXhHOJZ9F/WfSF7BCvPfgA/gXZi4zc1qHuoSZtzLzHgC/AHBKePz1AL7NzOuYeQLBxFAPl4Xn15i5ysy/YuaNHOAGAL9FMNml4a8B/CMzb2HmMoBPAHitaeWegdeHK+Hon+kkInoDgD8H8BpmrjZQvsXMxHIARQCvRdBHTwFwKoB/ks6pArgo7C9XIBC8X2DmUWZeh8C1dZJ07lFEtISZx5j5j3Xq/+dwnN6AQIC9PuW8KgJrQZWZrwYwBkAW+D9l5tuZuQbg+4jH88sBPMbM3w7H510AfoJgfLkAXgPg48w8zsz3A7g8o61VBILqKGb2mPlOZh6Rfv8uM98fLtb/GcGYdFPKOh7A9QjmqUvDY28BcDUzXx1aRK4FsBbASw3Xn4ZgsfDhsO1TzCwIeW8G8HlmfpSZxwBcCOCNDcwnN4dt8AB8F8DJGefWeyZ9ASvEew87mXmqzjlPSZ8nEKxWgWBgPCH9Jn9Og3IOEZ1HRH8MzXf7EAzCJRnXHwrgfyXh+yAAD8EEmxdXhtp29E8/gYhORbCqfhUz72ygbIuZC7HI/W9m3sbMuwB8Hqrg2C2RssT527UyxPh5B4CjATwUmlZfnlH33lDgCTyOYPyZsDsU0ALymAXSx/OhAE7XFrdvRqDZLwVQgDp+H89o73cB/AbAFRS42j5LREXpd72cItLH/ZsRWMJ+LB07FMDrtLY+D4GVQMfBAB7XnonAgdp9PI7gPvPOJ/qzHMxYANR7Jn0BK8R7D61sK7cNwCrp+8GN1EdEJQQr/c8BWB4K06sBUEbbngBwniaEB4Vfsh0ITXL/i8ANcXe7yrXobzDzXgBb0NqYkctbz8xvQmC6/ncAPw5N2iYs0n47BMDWdrRDwhMAbtDG1lxm/hsAOxGYveUxfkhaQaEV4JPMfDyA5yDQ8t8mnaKXUwWwK6W4T4S//UDS1p9AoM3LbZ3DzBen3NchKcJ1K4IFgdyWGoKF1ziAYfFDWHfCXJ8BpZ/keCZ9ASvEZxauBPAXRHQcEQ0j8AE2ggEAJYQTBBGdB+Al0u/bARwgiCYhvgrgIiI6FAgELhG9suk70BAO9J8A+D4z/6hd5VrMGHwbwPuIaFnoZ/4gAn9rwyCitxDRUmb2AewLD2eFVn2SiAYoCGl7OYD/aabeDPwSwNFE9FYiKob/nklEx4XWhZ8C+AQRDRPR8Qh8vEYQ0YuI6Gmh4BtBIKTle3sLER0fzhv/CuDHGWFlVQCvAzAHwHcpiGr5HoBXENE5ROSGpLEziGiV4frbESgcFxPRnPDc54a//RDA3xLR4UQ0F8CnAfwo1NofQaBZvyzUmP8JwXyVF9sR+NnzPpO+gBXiMwjMfA0CQsrvERBk/hD+VM55/SgCotqVCAhqf46AHCJ+fwjBIHs0NJkdiIAMdBWA3xLRKIA/Aji9LTcUYBUCf+cHKWaojxFRqtZhMavwKQB3IJjgHwRwN0JiWBM4F8A6IhpD0K/fmOHaegrBGNmKwI/9rnB8tA3heHwJgDeG9TyFwEIgBNd7EZjenwJwGYIFTRpWIDB/jyB4TjcgELwC3w3LeAoBEfD9ddpWAfBqBFaLbyEwr78SwMcQKAFPAPgwDDImXBy8AgH5djMCa8obwp+/FbblRgCbEBD33hdetx/AuwF8I6xvPLw2Lz4D4J/CuetDqP9M+gLE3BZLlEUPgoKsaPcDKKX4n7rdnrcC+BqACoBnc5jwpYXyViOYwAcAvJuZL2u5kRYWswxEdD2A7zHzN6a7LRaNwwrxGQYiehUCpuwcBGxVn5nPn9ZGWVhY9CysEO9vWHP6zMNfIzBnbUTg3/mb6W2OhYWFhUWnYDVxCwsLCwuLPoXVxC0sLCwsLPoUuYQ4Bflp/0RE9xDR2vDYYgpy+q4P/y6Szr8wTIP4MBGd06nGW1hYWFhYzGbkMqcT0WMA1oQZkcSxzwLYw8wXE9FHASxi5n8I4xV/iDi13v8BODoj5hBLlizhww47rKUbsbCYDbjzzjt3MXMjCS66DjueLSzyoR3juZWtLl8J4Izw8+UIcun+Q3j8ijCP9iYi2oBAoP/BUAYA4LDDDsPatWtbaIqFxewAEWWl1uwJ2PFsYZEP7RjPeX3ijCCZx51EJHbWWs7M2wAg/LssPH4Q1Dy8W8JjFhYWFhYWFm1EXk38ucy8lYiWAbiWiLIyE5HhWMJmHy4G3gkAhxxik29ZWFhYWFg0ilyaOAf7RoOZdyDYiOI0ANuJaCUQbA6PYB9qINC85WT6q2DYGICZL2XmNcy8ZunSnnbxWVjMapiIrRYWFr2BukI8TFA/T3xGkMv3fgT5skXC/QsA/Dz8fBWC/V9LRHQ4gNUIEt5bWFj0L17EzKcw85rpboiFhUWMPJr4cgA3E9G9CITxr5j51wAuBnA2Ea0HcHb4Hcy8DsEGGg8A+DWA92Qx0y0spgs/umMz1m8fne5mWFj0BG5avxM3PrJzupvREEamqvji79bD82dv0rK6PnFmfhTAyYbjuwGcmXLNRWh+JyELi67g4z9fh7c/9zBceN5x092UXocgtjKArzHzpfoJluPS//jS7zfAZ+AFR/ePe/Pm9bvwud8+grOOX45jV8yf7uZMC2zGthmGHSNTqHr+dDejL+Azw2YdzoXnMvPTAZwH4D1E9AL9BMtx6X/4PtBvabhroQbuz+IpzwrxGYRyzcOLPnc9/vfuJ6e7KX0BnwF/Fpvh8iKF2Goxw8Dov0WtWHT4/dbwNsIK8RmESs3HeMXD3vHKdDelL+Azw8rwbGQQWy1mGHzuP2EofOF91uy2opWMbRY9BiGQrGCqDw5N6ZxMYWChYjmA/yUiIJgvfhASWy1mGPrRvSSEeL8tPtoJK8RnEKxpKT/EI7KPKhtpxFaLmYdgUdtf4Ehx6beWtw9WiM8g+JFgmr0dOi/EE5rNg9/CQgZz/9mlvEhxmeaGTCOsT3wGwbcdOjd8a7WwsFDA6D/LVOwT77OGtxFWE59BEAJpFvfn3LALHgsLFf3oE2c7jq0Qn0mw/qH8sD5xCwsVvt9/PnFLbLPm9BmFWBOfvR06L+yzsrBQEZjT+2s8eFZxsUJ8JsGGmOWHbwe/hYUC7mNzer+1u52wQnwGQWQf6z+OafdhfeIWFir6MW+CNadbn/iMAltNPDesT9zCQoVvQ8z6ElaIzyCIITibV6V5wdYnbmGhoB+FuCXzWnP6jIJvtcvcsD5xCwsVHP3XP7Bx4lYTn1GI/Lyz2baUE9YnbmGhgrn/hKFntyK1mvhMQmQinuZ29ANsxjYLCxX9mHbV7hdhNfEZBWsizo+I2Da9zbCw6Bn4/chOtxY1K8RnEmza1fywyV4sLFT0Y9pVLzSjz+ZxPOPN6Q8/NYo3f+OPmKp6092UjkP4hawmXh8Rq3UW+9IsLGT047Rhc6fPAiF+75Z9uGXDbuwYKU93UzoO6+fND/usLCxUMHPfjQeb7GUWCPHZRHywCUzyo1s+8Ue2j2LfRKXDtVhYtI6+3Ip0Fs3vaZjxQryXyV7/98B2PPfi36FSa49N14ZN5Ue3fOJv/eZt+NqNj3a0DguLdsDvQ03cKi6zQoj3rmDbtGscT+6bxESl1pbyLFkrP7q1Wcx42cNEuT3v18Kik/C5N+fJLFhz+qwQ4sHfXhRsouN5bRo5opTZ3KHzols+8UC76WgVFhZtAfdhyjYbYjYbhLjfuy9ZdECvTYLEMjXzo1vPqh9NlCYQkUtEdxPRL6e7LRadAffhgnM2cZ7SkFuI64OYiBYT0bVEtD78u0g690Ii2kBEDxPROZ1oeF60W9ttJ/w2pwy0udPzg7tkoelHE2UKPgDgweluhEXnEBDb+quz2tzpjWni+iD+KIDrmHk1gOvC7yCi4wG8EcAJAM4F8GUictvT3MbRy8Q20bZ2aeK+7dC50a0Fj+9z3+eyJ6JVAF4G4BvT3RaLzqEfXT9elBtjetsxncglxFMG8SsBXB5+vhzA+dLxK5i5zMybAGwAcFpbWtsE4i0np6sF6YiT98eNe88P7sKnfvlAU+X18oKl19Bdn3jfv4//AvARADMyNc5/X7cez/r0dQCAC751Oz5wxd0dqed7f3wc7/3BXbj4modw8TUPdaSOVuD7/ddXO2lOH52q4pz/vBEPbB1pe9ntRF5N/L+QHMTLmXkbAIR/l4XHDwLwhHTelvCYAiJ6JxGtJaK1O3fubLTdCn513zb8/qEdxt96OakHG0z9G3eMYePOsZbKm82r0rzonhDv7/dBRC8HsIOZ76xzXtvGc7fxH9c+gqdGpgAANzyyEz+/Z2tH6rlvyz7ctmkP7nx8D+56fG9H6mgFfchr6yix7an9U3h4+yge2T7a/sLbiLpCPO8gli8xHEs8Yma+lJnXMPOapUuX5izajK/duBHfvvUx42+9rJ2aiG1Vz2/af9/L99pLeOipEWzePQGgswI2tgL19ft4LoA/I6LHAFwB4MVE9D39pHaO55kKn4O+4PWoxsvcf3NHJ33itT4JX8uzAYoYxC8FMAhgfjiItxPRSmbeRkQrAQhVeAuAg6XrVwHozNI2hOdz6kvs5TjxSOhKjfN8RtVrzmoZxYm33LKZjY/8+D7sn6wGXzr4sGbCooqZLwRwIQAQ0RkAPsTMb5nONvUr/EiA9+ZuYf25FWnwtxO8E7FA6EVStIy6mjgzX8jMq5j5MASEtd+Fg/gqABeEp10A4Ofh56sAvJGISkR0OIDVAG5ve8sl+NoKUt7spFss5GYgOp6sidd8Rs1rVhOfEZpfxzFR8TBRCfpIJwVsLy8gLboPZkRaeC+SHfV5tB8QJ3tpf9n9oom3Eid+MYCziWg9gLPD72DmdQCuBPAAgF8DeA8zd3QLMd/n6GVOVjwc+8+/jogjnXzJrcIU/lbzOOo8jaLbO3MxM25ev6vvFg1yf+nkAPUMi7R+BjNfz8wvn+529CvEVp9CG+819CU7vYPcFi+cSHv9mTQkxOVBzMy7mflMZl4d/t0jnXcRMx/JzMcw8zXtbrQOufONhSkuv3rDRoVtWe8lM6eb5DuFKDxCEro1n1FrUgp3m8R335b9eMs3b8OdPUjSyYLHjJrX+QHay1Ygi+7D56Dv+dybJtp+JLZ1MvpIWER78V3JmBEZ2zxJAMsT5m2b9uTyS46Xazjpk7/FDY90l1UbaeJS2zzfb9qc3q2duQSESXqyz/Zq9/zY2tHJZxUtqmZkYNbMQ+cT/4Q+8Z4ltvVmu7LQSYtav+Rlz0Ns63nI5lFZIG4fmYrDrjIm0v2TVYxO1bBl72RH26nDaE73mzend9snbgqR6wf4shDvik+8v57PbEWlSUJpXrAwpzP3pPbE3HeKOIS+01GfeI/Pb73YlxqGHIsrP2850UbWRDpdhDDTSq/mxabeRtGtnbkEvOi5dae+dsFjs0983db9eM5nrsPe8fbs/93t92HRGqaqnRXivi/M6dyTC99+TEzkd0ETb9Iw2jXMCCEuh5jp4Vp5zOlCS+/2wIrSrmptbp7Y1l3Nr19DqDxfzpYXH3905zi27p+KEn+0ihkSJz5rMNVht5BqTu9oVU0hiGOf7lY0hk4qYFYT7yKY2chSlFeWWe94ukKBfEMnqbXgE++25tevIVSeJLn1/qIfa62emcVOn+mYrHRaiAd/az3oE5eFYD8tOjsZfRSz03v7ecwIIe4xRxqV/DI9Xwq7yngRnQxTyIJObBMrdMFOn6jU8IxPXZubcNdtt0AU596BEeT7jI/8+F6s27q/7WWntbfdRDRrTu8vdJqgKcZlzes9c7o8ZfRY0zLRSd5JrU8W4TNDiPtmkpjse8rqmNO1J62ntbmm/d03UcXu8Qo27x7PVV4eq0M7IZ5pJxYNT+6bxJVrt+Cd38mb7Tc/dN5E9NlPHmsF1pzeX+i0EBf9qur5PWeilVvTT/21kwtl0wZVvYgZIcTl0Ai5A+Yltk3Xdna6lSAS5lp8Yt5Vex6rQzvRSXO66wQp+DuhsagLPem4IeSvFfQrZ2C2YqpL5vSq5/ectqu6laaxIQ2io7nTvf5wF84IIS4SKIjP0XGfcyXcMGnx3UAspIPv1VAVrEV/hVDJV163Q5q4zUJPRiTEO8E6NfjBgfZrzjZOvL/QLU285nPPmWhN3JB+QCfnvBmTO70fIO8KJE+YQehZ/dXUtIWYaYsHz9M18ZBYkbMTdT3ELHzWHY217pC/PYLGoZD/tlyPjRPvK3TeJx78rXm9lzu9X7toJ62BsyF3es+AOUzgoW3jmTftqi9N3g9uG8FzPnMd9rQpVjgLui9e9okzc8PEiugeu+YT78IquIuaeLvvp90+dovOovPs9KAfVDy/5/qESmzrrbZlwevgGOsXdvqMyNgm4sFfdsnNWL18bnxcNrNnLNXkyXvjzjFs3T+FbfsnsXjOQMfbLf+Vc6Z7Pjeeu3e6fOIdMBd3ypTlSy4WwExya5eW1K8heLMV3YgTF+g1E628CWmPyywFYqx2os2REtXj7rAZIcQFge3JfZMYHHCV43kS5MshZqYkIJ2CTnyS48Pl9Kv5zend9okHfzuhLUckvTZPdnpbVXZ6e4WuNaf3FzrvEzd/7gWkRWz0Otq98JbRL7nTZ4Q5XezP6/lqytK8W07KZu1I++3Ci9PjrOXVedXzI3NOfnO6+rfT6CQztN1M8ahc7eHIxYs1VPvZ6W0pzqLDmKx0Pne6QK8Jhr5lp3fDJ97jD2RGCHFhTveYUfVkk3S+iVQmNNW66AdJhJa1aE7verKXDg6gTllE9PdqylRl48RnJ7qpifecOV1uTm81LRPdyZ3e2w9kRghxwUL3JcEHxJsNBOfU94nLZLJurL7iHXiEEJc1cW7YLxwJjTa2Mbu+4G8nJqROha/pbTVNrO3ziYu/vT0JWAQQPvEwurHt6OUwrl62EmShkwmnojjxHltw6eh7IS6vxDzmKNYaiLf+E59zldHF2MB4K8/gu7IA8Zthp4u/3dXEO2pObzuxTftuMCO2q0o9D4BFb0Ow0wtuZ6bFXvaJy0O4x5qWie7kTm9/2e1E/wtxabJnTgrCfHHi4nx01Seum2tUTdxvWDPsdnKRTvqjOrWIShLb5M/tXThMV/4Bi+YwVetO7nSg98zpvWwlyILNnT4DhLguAKsJc3rwOe8GKJFPvAuCUBfS8u5aMjs9rybXfU28c/V16haSk6c0ebWZqMddfh8WrUFo4p0yn+r9oJfMtH3PTu+gItFL78mEvhfivmaKloltzPlesq8I8c6t7HToPmXViuCj0WQD3db4uM2aq4xOaSqJiVT2ibfZD2/jxPsLgtjWqbGf5cqZbshx4v1kT+9k1kibsa1LkHcGAqCEmAUm9vralcxK9rppTtdMQWqImayJN2hO75Ym3sFEC516/jVfF+KN+8R3jEzhsI/+Cj+9a0vmeTMhTpyIBonodiK6l4jWEdEnp7tNncJUJMQ7IxT0ftBLZlq5Kf206OxO7vS2F91W9L0QzzSn+/mIamqIWfdMKLrwrkp11rzGY9a7HifeQXNxx0yauhCXvvs53/0TeycAAN/94+PZdYXF9NBc3QzKAF7MzCcDOAXAuUT0rOltUvsgC2s5xKwT3U/vB720MY5KbOufDtvZ3On9kXa174W4PuHKsda+4hPPKEMJMQsTrHRDiGvECdUn7seaeM5tzLrtg+3kLmadevzZyV7yrernlIJEh+PlWuZ5M0ET5wBj4ddi+K/jNzRV9bB++2inq1H6mZw7vZV3dufje/H7h3corj1TmfXqGJ2q4rFd4023Q+Dx3eMYnaoafxsv17Bp13jHkr2Ua/ne4wNbR5pauLcrq9rO0TKe2j/VkbI7jf4X4vrqVjEL5TOny5nTYj9IvvqZGfc8sS93e01tFfXraVcbz9jWOfN2t+vrFjtd/paXnU4IAonHpuoI8S6GK3YSROQS0T0AdgC4lplvM5zzTiJaS0Rrd+7c2XKd/3v3k3jZJTdjopL9jFuFvOifqqquuGawb6KC13zlVvzFt+/AjY+oz6FRc/qlNz6K13/tD021Q8brv/YHfP3GR42/XXbrY3j1l29RhXgb++vP79mKl15yE8YyFrxb9k7gpZfchOsf2dFw+e1y6X385/fjb390j3Ks4b0rpgl9L8TrmcnzMKhlbd2LNoLP9+Lu2rwX53/pFqzbuj9fg5X2qT4X2V9b8xpPPNPtDGH6IqS9ZXfJnG6YvOpVLd5b1sQUlI1c5fU6mNlj5lMArAJwGhGdaDjnUmZew8xrli5d2nKdI5NVVDwf5Wpnbc56qmOBZvvfhKTN6/0jaU7PrmPPeAUjKRp0IxiZrGEkZcG5f7KKkalax/royGQVVY8zN5fZN1ENz218wdauCJk94xXsnVB3rpwxmngasYWIFhPRtUS0Pvy7SLrmQiLaQEQPE9E5nbyBehp2Hp+JJwm/Rslk+yeDDjhaRysztk8zRytC3I/jxHvWJ97BRAudGjiZG6Cw+Rwd4pp6QrzdaVynG8y8D8D1AM7tdF3iXehExHZDLl/ZxrjJar2U8oIy9e/ZZZVrflu0QLFBlPG3cI5UiW3te+Z5XHzlWvO+53blqijXfFRq6oIxVqJaK7vTyKOJpxFbPgrgOmZeDeC68DuI6HgAbwRwAoLB/mUick0FtwNZE67qE08/T/btNrr6qrZgctF34FF84l7zudNngk+8Y+b0xMSa/K3e88vrculkHH23QERLiWhh+HkIwFkAHup0ve1OvJMGmW9S9Vo3p6sL8WyhXa+Ocs1vyyJGzkSpwwsFuNxH29ldIyGbIQjLYZKdpubQNmnL5ZofLSYEZkzu9AxiyysBXB4evxzA+eHnVwK4gpnLzLwJwAYAp7Wz0TIyzemccytSaVJulNgmBG0zg03vJNWET7yxDtpt821cXz+Z09Xvxg1Q6rzLxhdVDTSw97ASwO+J6D4AdyDwif+y05WaNgXqBNKEbrMuIs9PXwg0SmwrV71AwLbYgeTMlTp8w3Nu59jLIwhb0cTblRWxXPOixYRAnPirtwdwrv3EQ036TgBHAfgSM99GRMuZeRsAMPM2IloWnn4QgD9Kl28Jj3UEWe9OMadnvAjZtxslvc/53mKh3/hkw1K9QRnShOL5DZv2u22+7aTm36m5WxcKZnZ6dhl579dr0wQznWDm+wCc2u16u6aJpwnxJt9Zliae8InXE+KhcKv5jIEWdmXxOX08mdwW7XzkeRbGgvfQaL0sW1pbnC/K1ST/Ysb4xIF8xBYJpt6WeArtYrOmDfKCQ/A4n99WFkaNptoT2nMtZxiYDH2TDyV3usxOz512tb7VoZ2IFx/tL7tTJqwsbSjyibdJE+9kRruZDt8wJjqBmuLCkrToZoW4nKeiTohZfXN682bmqM46mnDEx1Hmr3Zq4mE7MjXx5u5TPr0T5vQ4d3pLRXccDbHTNWLLdiJaCQDhXxEfsAXAwdJlqwBsNZTVFjZr2ssruJTbJ24KMcs7iL0W4soTxDbFJ+dL/vZ8UrLbPtiOmtM75hPX6pEngpwrb7ltWe/djyawxtpokVzgdqyeFA20Wc0uTbMPyte+16kj1sSbXyV7dTRhkxWwnY88z3ts1pzeztj2cs1DxfONY7vXzel52OlpxJarAFwQnnYBgJ+Hn68C8EYiKhHR4QBWA7i9ze2OkPbii64D3+dcJuZY+DXuE6+24BOPJnnDQFL2E89ZdLd9sB01p0tFtnORYHqveh/JS2wDkJpEo5HyLJKIzLwdVoPSxm07zOlZJMo8dQjzbkuaeJ0+KI7L7W5nd+Ucc1IkxBu8T09pc+uaOABUPNky0x/jN48mnkZsuRjA2US0HsDZ4Xcw8zoAVwJ4AMCvAbyHmTu2x1+aKbfoOiGhI/ie9SKUXcwafHFCe26LJq5NAI3HiYu/3el0ug/57s178d/XrW9r2UB7Taqm9yoOeTk1Z7ltIsTQXJdavkV+dCtRTipru2liW3q/1cdlPWufMDO30v+FopC2FopyVLQhRt5cfv35tBzGkDe6XlM18ebbzMxReFnZkPCn191hdYltacQWZt4N4MyUay4CcFHLrcuBVHO6QwGhI8dKUNkARdOO60EMsJbY6YaBVPPiXcx6ldgWC7/gw6u+fCsA4H1nrm65bPn5Vz0fRbc9eYlMz9JnhgPK7cOW25aVoGKmxYl3E11jp6dIjuY18Sx2ulZHXZ9465p4XXO6oc+3VYg3Yk5vySfeeNv0+oPPHoIALJs7vWtI6xxF18ntE1dDzNJN2PsmKol82Y36rWXEzMrkQkDZxazRELOGW9Ic2r3/tlK2VGa11v5JRa1L/a3e/ch9LksTz5uL3SKJvCTDVpFWfrPVNuQTr1OHzE5vFtNvTlfrMaFZn3i7Fh6qEE8uwnpcEe9/IZ5FbFO3Is0qI/zrxxugmFaF77h8Lf7tVw8ox4TwbsonrpnT9WxPjRIr8oTTtRMyl0BGO4S6wg9oozZm2kxGt9bUY9vLfS4rt3fcrxpro4VZuHQCaZp+O5K96At7fVzWZacLM3MLvIB6bgmTybidQjyew9LPidjpjRLb2tRmOT5c/txoiO90IVeceC8j7fkG5vScceKS3yYrOcFT+6ewYKioHGslY5supOVkL1V5F7OcPTT2iTfclKaQtsr3fEbBbT6uVS9T3w2qFZieZaQt5PDfAepEXclomzWnN49u+SP18h2C4oZruDwtYZMMvcj8ceItsNPr9GlxXB5j7dyKNFeyl2qz5vQ2aeKSH9y0CU6vj9++18TrmdPzkJVkv1DWi5uo1BKJ/MUAa4ZFm9TEfQwUnKi8aDOWRjXxLnW6VCHehvrlW26nOd30LMWklff5yX0ua4MOy05vHiayZyegl18IuRdtYadrc0IjGduYuS0+8TRrWfS7IL4ppummqzPUX38MxOb0xsr22iXEU8zptT4R4n2viaeZbgU7PY82FPtCkckIn6x6SSHeSu50Tduo+YyS66AS5kxulDQ3bSFmmhxrh/bUMXN6lk/cMGi/dfMm7Bor48XHLsPC4SKOWjZP+T1LE7dx4s2DpYVtJ6H31aJDqBiON1NeK7nTlVCnNvjEs3Kn63W0U2jlsYQ2nexF6hqtEdvM5nSdeNyr6HshnpqxzSX4ftyJsvy0EcuaY594IimIz5iq+oq5BWgxTjxaJTM27BjD7Zv2oOASii41xU6PiW3d0sTF36Q5vfWyO2RON7Qtyyd+4/qd2D5SxrUPbMdxK+fjkjedqvxezthiUb4HZgZRay6G2YRoYdvlOPFiwQEqXtM8hmx2en5N3ESwagZ1zemRoiCZ09v4yD2Dpq+jaWKbNr6aRbom3h+50/vGnD5RqeFzv3k4kaQ+zXQb+8SzzwNUX2ikWWvnT4aT9VRKkvxaE4JGJrSd9fkbsG7rCFzHgesEpLxqg+acPIkV2om0lI5tEeKyJt5Oc7rJJ66lhpQnhMDF4odpGYXGEL/rTE1cMfe11OxZhzyTf3vqUd9fwWnNnJ4dJ66em3VrspumFU08r09cXiy1M9oknlvTz4lzpzeqibdnfMnPWokTn0HJXnoCt2/agy/+fgP+tGW/cjzt+RZdBx43tp+476czwoUQ132grcSJxybD+JjrAEXHCTK2NboVqa+W22nEQk893hZzulREW83phqJ0n7i6EQ1HO8qZzGv6HsQy5MfQ6wzXXoPowx33iWuafjEkZDadO132ievs9AYWu6pZt/n+H1kZU83pyd/b+cTz+cTF4ri5suuVXw912elWiLcHaVt+ppvTRdrV4HuWYJPJH2kr18lKqIknfOLNk09Mde2bqKLgEmoSOz1v0e3yiY9OVXHjI/U3pemoOV3RxNsnxE2msaRPPP6t5vth9jxztIC+aYJabnsmmdmI6cidDiBKKtT8VqQclkM5cqfnM6e34lKot593FBkja7VtfOadzJ3eLjJe3TjxHl+A948QT/GR5c7YliEHZHO60Pr01ddEqhBvXmPQhQcQdCLXcULtr0mfeIsC42d3P4kLvn17Zl7woD7zoqEdK1e5jGob/aKm96RbFOQ+VfMDF0ta3H5eTdzK8MYgkz07iSQ7PdDEm61WlFcquPUztuUwMQMd9okLIaukXW26umT5OVyC7cjY1hmfeHuUok6jb4R4WhpGU+dwCHAo8CvnMecoIWYpYV0iqcdUzVc6jFjBNmryUv05jDkDbvRdENsa3k8c7el0k1UPzNkCCpAGqFZhO8hICrGt4+x09TkrQjwU4EKYy+cB2Zq43E+sJt4YYjNwt9npjvF47vJCYVgqOFoWtGR5mbHTBrNuM4hcdil1mRZL7STGplnrZMS50xsV4u0ZXzI5Vf7cL7nT+0aIp+0uZvKjOERwneDF5gnzMYWY6eUKn7jns6IZCnN6owNN7rCezzhgbin6PjzgYrziNZxsQD6tlZVpXlJR2gBth8Dqpjld18B1clItXNzpfXCw6GSb0/32TDKzEdOliRcLQhNvrt5YE3eUOHHTbeQ1p3c7d3pb2ekG3o+OSpOauGJOb2GKqMtO7/Gx2zdCPG1QmzVxihjejYaY1fOJAypDvdnUfL4mxOeU4mi/hcMD2D9RbTpOPPjcUHMUeCkuhbT69GfVjolXHvTtnMhN70kX4opP3PMjf7juEx8eKOQ2p/f4Yr7n0D2feGfY6aWiWzf2Om+IWXvY6ebfxfFO5U5vyJzeYL3tI7bVy53e24O3b4R4WlIV0+rNcQAiLe1qxouQX1Zalp4JWYhLJpdmNXE1UQErIWoLh4rYN1lpmFihCo3WNfF6ZvHY5Km1ow0Tr5JQpZ2aeMZkGlkWpPYLF4vKTg/+DhXdRMhjWl3dihiYKYjY6d2OExfs9Cb7sKKJS4Pc1O8yyV7VLrHT/eRznrZkL42a08PHUnCoLbnTidJypzdfdjfQP0I8fKB64g/Ti3eI4JIgtgXHssak3NFMbPMteyfwpyfj0DZTDGejmxTomrh8XwuGi9gnaeL5dzFrz2o6L6FO/K4LqHZozjq5rF3ISvZiWnmLHPZV34/3mhdCfMCto4mr79giP6Yrd3qxxbSrniTE62m33dXEzWXEfT97wdF8/eHfTJ94k+b0sMyCSy36xIP655YKSltMRNdeRN9kbEsTLKb37krm9DxpV+WEBCZN/Hn//nvlfFkTr7bDJ84BA/v5q5fgs689CV+/cRP2T1SjzVZyE9vaZF6KTJlNmtPbnXa1mUQ6qeUa7kkcMiXL8bxggcUsCZbInO7WSfZi/mxRH2JN3O048Sh3etMZ24LyRNpnAbMFKL2ctvvE6whxldjWPuQJe202xEycX3Ccls3pAwUHg0XXuP1rrwvxvtPEEz5xQ+8gCv7JGduy2enBX5PJ1AQ59Wps5m9s1LNsTvcZFc/HqkXDWLlgCAuHixgtx5ut5DanS2W20u/y+vnTrBztSbsaf662cSI3x4mrwjkRYhb1vXghSRRoW3k2QAH615xORAcT0e+J6EEiWkdEH+hGvdOZOx1oPkzS830UHApzPciCoH7dMhSzbgsuhfg5msswkVjbmrGtHrEunPvktuQu24818VbN6aWCE4xnQx73Xrei9Y0QTxOuRmKbE5rTfXkr0vSyZY0yT2ygTGyrRhN8C5p4aE4fCP1xC4cDDXzvRDVxbhbaRfQQroH6PvH4uSmpStvBTpcFaRs1cdN70n37aaZ8eXHjEqFUyNbE5cfQ4/NAFmoA/p6ZjwPwLADvIaLjO13ptLHTW0z2UvMZrkMoOHk08Qwh3rY4cVGX+XfTNshtZafXUQgqXvNmfFFmOzTxUsENhbjQxKV29fjg7RshnrbRiKlzROb0DLa5DPGb3JGzXtyUgXTSCjvdZ0a15kcTiDCj7xmvhG3JW6a5/EaRlxksa67ys9Ov27RrHGPlWmNt6NBKOEsTN2kt8gJCzqvvOIQBaeVer65eN8mlgZm3MfNd4edRAA8COKiVMr/0+w34918/ZPxtZKqKz1zzYDSZ/mHjbrznB3fhB7dtBgBcdssmvPcHd+GeJ/ah5vm4+JqHsG8iGCebdo3jy9dvyNWGPeMV/PuvH0q8P1Oyl3LNw2euflDpw3dv3ovv3/a4cu3ltz6Ge5/Yh4ITzD9irtq0axxf+l2yXQr3wvPxmWsexP7Jalhn8z7xWzfuwk/v2gIg6ROfqgb3IvJemHziP7x9M+58fI+x7I07x/DVGzbWbcNYuYbPXP1gbE1kxh8f3Y2f3LlFOU9erJjGyO8f3oH3fP8uXHXv1ugYM+Pzv30Y2/ZPAQjIiKZH9LO7n8StG3YZ27dl7wS+8H/rgy1fq36oibv41X3bcPP6XdFzG3CdqOzNuyfwxd+tN1oqLrtlE+4PeVPrt4/i3d+/E49sH017PG1F3whx0cl0rcz04okIjkPwfDl8KL1s0X8rcr5iqVx98ymTOb1xdrquiXOwgxJiIW5qSxbkztWKyNB9v2mIzOm+SjjUhe6rv3wLvn3zpoba4DNHTOFGMrbVPB8X/eoB7B4rG3/PsxWpfErNsJjwQ018INw2Ng3t2u+4V0BEhwE4FcBtht/eSURriWjtzp3ZKXvvenwvblpvPuc/fvMwvnbDo7jz8b0AgNs27cGv7tuGr90YCI4vXb8Rv7xvG35x71Y8vH0UX71hI25aH0zUV/9pGz7764eVcNA03PjITnzl+o14aJs60YqFtPzu7n9yP75246O447FYsL3qy7fiH//3fuXaT1/9IP746J5QE6dozrrm/m342o2PJtogT2UPbRvF1254FH/YGNxLK7nT//zrt+HvrrwXQHJh+v3bNuNrNz6Kr1y/UblPuZ//34M7cOUdqrAV+OW923DxNQ/VfcZ3bNqDr934KO4L97rwmfH92zbjkt+tV84zbf0p44e3bcav/rQN37n1sejYuq0juOR3G/DBH90DIHhnpvH1wR/dgz//RqKrAgD+v+/cif/8v0fw+O6JwJxedPCSE5YDAK5c+0RkISgVnOgZ/XrdNnzut49g1KCQfOIXD+Dl/30zAGDnWBlX/+kp7A2VsE6jb4R4mp/WSGxzgqxtzPnixMVLStMml0qJWAAzsa1RbVEltgV+ITGBLBweUM/1OZefSvHBtmCBjk1g2YXI2aAqKUQcZsa+ySr2TWancNXh+4xSIchi14g5fePOcXz9pk24OWUFbtwARfOF68le9M+eD7gOoVTMFuKKZaTHw1TqgYjmAvgJgA8y84j+OzNfysxrmHnN0qVLM8sqFdO5BOMpwkEk/JHHm9BWxd+0rI4mCOExoqUWFgtHeZEdWWAyFpPMcXsKbrATYb3rlDp89V7axk7XfNLirxDC0VakWhvT6hTPrV7Yp7he9nd7UoRHXF62Jh6924z5b3jAbdgFIFJKuw5F5vQPnnU0jlw6B16onQNBBEq8XWv9fgBIZn63O1sP940QT/ORmYSnCDFTdzGrb05P888snqMKVSHEn9gzEV3TsCYum+vCDhP5xDVNXD8/DfIprcWJmwd22nnMnKqJBwuQxv3aHjMGQstEI8Q20Y407d30XDj6LXlOTbkvP/rrOoEmPhvSrhJREYEA/z4z/7TV8koFN/W5pS2GI+FWlYR4VQi8YDw2knhJlKcL8YIhxCyPf16eO9yQ2FbvOlMd0T21PXd63Da5TeJ4cl41v594kZGtiYvrWRpTIoWxWl5cjqlKE9FMF46DRbdpf7oTCfHgvRecINOeuM/hgbhsL8XqmhZi6zrdEa99E2KWNkBNGqoTmdNldnp62b6hbD0Zi4ypmo+NO8dw5n/cEB1rJXe66KixJp4U4l5ImMksU9b82iDE84eYqZOYSZNtlGHuc/AeZbNkHkQhfymLhqw48cRfqf8EZcbPxQ194nnjxOu9j3/+2f2o+YzPvPppmed1G0REAL4J4EFm/nw7yixlcAnSNUA/1HZjgS0+x4I9/4JaXDMyWQvcImF/idjpRgtMernyoiTwicfEtlRmuIFAGd1TzcOcMP1ya7nT1boKWjKbNKtbPU08a/Fqut7neA8CGbJr0jTfREQzaVHuav7NoSaEeLSI8RnlqhcJccFlEPc5NFCAz4GcyWsNFsK+UGe+bhf6R4h7yZcJpGjiTiAAmJOmUhNM/VUZYFqd//yz+3HeiSu09jUqpMKVIMUdWQjxeYNFEOns5vrlq5pfQ81RkDfRhqy5ppnT6wnV1LJ9hkNQzJJ5IDRwPSmQqW1xXepvkY88ZYXt+UH/KhXqJXsxfzbh4e2j6NKYbxTPBfBWAH8ionvCYx9j5qubLVBmAetI1wADYRbzF/yE6bkxTTw2pw8WJSGeqYmnv2tZcxY+8Xppk3VejHwv5ZqP4VJB2UOhGehZH5OauHmRnWoR0awf9eqVv8s7M0blyeZ0Q52y5SUqSxuXQwNuw/OdbCUp13zMD62fgQUl7ltDRdEf0vtB8rv6rDuNvjGnxwNUfWCm+V3NnR4cy5KBRrKTLIh8H886YjFu+siLomPX3P+UWkaDvUicX3SdyDxflFaDy+aVjOdnQfGJt0BtM5FdTGBJY1U2hZG1mJzhaok2hJaHous0TGwDGjOnx7yJuO6gLPPE5vsM10HITm9PnHhV4kT0Epj5ZmYmZj6JmU8J/zUtwIEgr3iaTzytn1Q9VlIfi8kXkEyuDZBMxbWjUzUMSTsIxuZ0qU0ZiwPxXmWhJtjp9bglpjpkIS52NmzLfuJCE48sDb5yXPfzZllEAFWDNiG5ZXQsyNXyJHO6URMXlpckmVhgqOg2HNseZ+b0FXN6pIlXhTm9ELUtrR+kfe+WJt57s0YKxIPRV4yml+cSBVuR5vSJp20TeOXaJ/Dar9yKmsc4ZPEwDl48nFpGsz7xATdOzzgg+XqOXj4PgDTo8hDbpHHVigtWNjVlQdZcZY1UT1sa/G38+TgkkmY0YE6vozVlx4mrfSVtxe1xyE4vBBpc2gSihpipv/34zi047KO/ihZwQZ6AvhmOLUGY043jLqOfyPvbe6EZFDBo4jmEnkyGGyqq2wDr7cgy0+saNCBr4tnmfVMd0T1VPQwWXRA17qqTobPThZ9Wzwueh2sE5DenJ0zOfrALYNInLpvTTfUZNHGtjMAnntmc1PYJ03nsE6eQNCnM6W50flo/SH4PzusZTTwtYxMRLSaia4loffh3kXTNhUS0gYgeJqJz2tHQvKugoP6AoS7nv82aHEy/eT7jwW0juOeJfah6HK3QGykjC9HKWBLcsiZ27Ip5yu95Eg60K9mLzsZMry+uS/aJmzXxBs3pzHAcKGbJPIhZzPUFq0C8D7smxA3ahO8Hvj3XpWjgp01oclV6//jPax8BAOwKQ+GqNe4am3W6USoEsbem95r1rkcma8p5JrJb8FsOdroUYTKoCPGkOT1LE480aEkzLTiOponXtwqJviZr4qWi23D/16EnexFTjL5gze8Tz2dOTywKUnziEQu86JrN6YY0qHoZg8XGk71Erg5PxIkHfSD2icftAupo4ikWu0KXiG15aknL2PRRANcx82oA14XfEf72RgAnADgXwJeJyDWW3ADSsoiZNFTXCTTxvNvrmforc6Ad1fyAeV2ss6pqXBMPzpcFt2BjA7EmPmXwCaVB9aE31BwFeX2L8ipf9kGbwnPSfNRp8PxA2y04TkMLADGBN+QTlxYjQGzRMGpe4WQUZGwL3lda1rasRVVBi4Gv+r1pTu8ExIRpWvxk9TmZSS4ziFthpwPQzOnJRXMWyzzWxGOhltsnbmKnS/dUKqiLgWagm9P1vm3KnR60J4WbUI0XGZn1atcLIZjGTpdZ4Gp92WlQC6HLrdm9yEWoYqkosdN9lZ0etD/dzZammbu9EmKWkbHplQAuD0+7HMD54edXAriCmcvMvAnABgCntdpQYZbVO4dJODuhOV1Gw+Z0n1GtBccnK16kiV/xzmfhwvOONZzfuJACVCEuf14dCnH9/Cy0K1e3ia1vrk/UpZrT5U5dqeOjTi87yIpWcBsjtlXqEdsyfOKi2TEnIFlGzeMoY1skxFMmNLkqvVrxrmelOT2cMGVtWMD0zIVLaUTKNaCw02vq3NCITxyAak53RLIXqU1RvLehPxjM6QU3ZKfXiROXi0uy0wM/bcFxWtLEdXN69Ffb2jkhmFLarEcEpMG0x4UX/pPnJvHcBouucb4xa+Jx3aWCE5GYG4G8wJLN6bFPPGlOnxE+cS1j03Jm3gYEgh7AsvC0gwA8IV22BYY0jY1keALSV8Nmczol/BFZQtxoTpdinyueH63Qn3XEAXje6iWJ8xvP2Bb8LUqrNXkSP3r5XADAkjDRTC6fuCLEG2qOgrwajWySS9XEo4m1CXM6hcS2Bp5tWhSDqW0C+kSXtZd1zffjjG11zenpmrgQ4iL9ZbXGs0gTT39upj4nmMMjUzXpPD+hFTbDTgdg9Ik3rolrPnG3SU1cihMvFdyWNHEhOOW69DalLdjTfeLNstOlxZD0mykeW0BOoJOmiZeKLhxKji/OqdBEmnhoHRIhrQlzus8ScbJP2en1MjbJpxqOJZ5iIxmeANWHIcMknEXGNvW89LJNv/k+oywJpqLk3zh8yRzMG1Sj85r1iadp4sMDBaz75Dn40EuOzl2+fEpLPvGc7HTZh5ymicfm9MbNXfGWstkLgLs378UTeyaU+tJM3KZmiEelZ2xL40oI5vxAHU08W4gHHXS8HGvixcJs8Ymnm9NNfW5+ONaSmrhKBssStjpkTXLQyE6XF6Lp/aGmEdKABtjpfrIOxZxedBSCXB7UNG5KImxSs2iK42mCSIcpo5yxHbomzmx8P+K5DRvCxBSeTQrnplRw4DjJ3Ome4dka2+n5qOjsdM9kTm9AEw/b2lOaeErGpu1EtDL8fSWAHeHxLQAOli5fBWArWkRa+IhJsxLJXmRkrcbSwo6qmolMYHiggD9ceCbOPn55dKzhEKqI2CYLcbXNc0qFhvY3blecuLiXen6mWPhpyV5kdnqTceKeH2RTKjhUdwHwtz+6B5dcF+RkruTQxJNWmvg3+bvRnC4JcSGM0oS4fMv6oxQDXGjictrdmY5YE09qc9maeCDEA8Z20pyeFoZqgiyEBgsGdrqiJXeInW6oQzenN6qJ60mXEmGT2t+0zY7azk5nhryVb1xeujld/Cbed1S2NLZLBSfadlpGzVCHCSJsMfKJuzo7PVhABlyYvOz0HtPEMzI2XQXggvDzBQB+Lh1/IxGViOhwAKsB3N5qQ2spPnGTmVmkXZWRrYmrPw64DjxW/ar6BDu3VFDMcI2ai4XAlcPKioXk62gkxEw+ox27mDWriXuGVXOj7gbmINlL0a1PbBsrexgPhWFkTk8NMfMTiyXdJ55lifB8jjO2uenCSNyDXodAbE6XNPHZIsQjn7iZc6Bj/mAoxEN2+pyBgqqJ63HiuULMZHa6E21yZNqKVNfATMRNNWNbIHx9Vk3aAqKubE08MPEWGkx2JD/Tmu+rm/BIbdHvRV8op2rikbm/DjvdwNg2auKSENfHiKhLvG9T20oF1+gTVxYKGW0V80bMTg+JbVXRrliJatwn3jtpV40ZmwBcDOBKInoHgM0AXgcAzLyOiK4E8AACZvt7mLn+tkJ1kCYMTH1NZPpSz0sfCMm4QyeRwMRkGhEv2FRGPQi5NKcUvwITsUlYFBontjXUHK1t+TSa1K1IpbojXkGdlbupbFfb0jEN5ZoXlS/aUamlaD++SLAj+fCFNlInxAyINXGHKBZGeczp2j0IU7zY8KPmcWJxMVPRKDt9/lBoTg818eEBV5lom/OJy+QoN8orYUr2ogsfU4pheVEgLEhA0Kf0/usSgclchywkSwUHrtuYJq7EXWsLCNkkrLPTk4IohZ2e05yeZKeb348glRUcSgpxibm+Z9y821mp6Bh94g1r4lKcuFggDoTtCtpvXoQA6jwhp2ftliZeV4gz880w+7kB4MyUay4CcFEL7UogzU+ZZk6nhjRx9fvwQAGer2qXpjjxQUUTb0xqik53gLRDmkkTExaFPJq1muylBU0852QYmdN9fStS2ZfVnCYuBGUeTaRSi4ko1cjklS5Y9ecs7iMmuCG1zTXPhx8uMMT7n0pZ6cuXp5rTy7UwEcZsJLblY6fHmnggxOeUhCbeeu50IBYEHuSMZiZN3E9ca4oTD1w2TlSOPo4ch8DMxtTOijm92Dg7XX6mNV9lgnusauLBLo/i3Lw+8dZyp+t1ib28iQi6wU3UMadUwI7RMpgZRCpHQLDT9fnRZLI3YbwsNPHYJy7M6aJdorw0d41cV1Xa5KWnfOK9gDSClEm4OQ5BV2oyfeJahxsKmZIVxZxu0sRjIZ4nS5SpzgOkHdJMdejJGTLLbJNPPC9BKNWcLvXxZn3izMK36GS6KgSDNRLiNdFP0vzUnBhcut8wZqebfeI1LyDdCV9qWgrKrJA/sUofr3jRwmP2CPFQEzc8tzw+8UATl3OnNxMnHgs7ebJ2HUr4WPXxYNoDW9061I/6WM1PauIOBYpGWu500aebYadnauK+uqiWi80T9SOzxZvNna7/JpLaiORcyr1EqU/jWG39+lLBBVGS2KYsFDLaGvvEY3Z6LUrF6ipKVHqcuFpXzWcQIcHL6hT6ZtZIM/Ga5mqTOT1rIOgLAeGfkQWByb8xKPmwG9fEg7+qEDeY0ym/OV2+jbbsYtaIEE/RxKspi688bXAoIJpkXVvRzPVxshfzNbImLrpIIk48YxHj+Sxp4mqstw71fai/ibInyrWorbMuTjwnO31eSbDTJZ+4kuxFJUw1GicemNODz27IpzFlU/MM/u94H3BVsIv5x/OSm36IqAsTA75c86M+LczMjfBtdCuBp80Joi3Vmm+8x7TvQDCmxCUNx4mnCEE5qU2aOX1OSC4zEQUDTTz4zIaFV722xj5xTROvxsTCoP35fOLlmg9PWsR1A30za1RTTGXmEDOTOT1LiKvfhc9NEeIGLbkka+KNmosN5vQBA7HNNZj30uAzKz6cZhEJ8TpliLnFZ10DiM9pNk7c4zhOPOtaIbx1n3iWJi4sHq70rNQ85+YVtzgmJul65nTdHylDtG+84kVRELPHJ94YO12Ec0aaeEn4xNvDThcmWSDeAVHuPrqZvr4mHqfQrfl+oh85YUZJUx3lqheV1Qw7XWmblzSny/eQFntt+q6X3VTu9GjTEU0TF+b0hBAPNfGSq1yX9InHglZAXSiobZWfyUQ5zScehPgJMRJYNcwySPe/1/z620a3E32zFWmqTzyNnZ5iMpXxmWsexIkHLkiUMVR0sYcBT3r59czpDSczEeb0uXU08UY2QOEgkxh8bo3YJnzDdbTnOJVj0OnF/tpGn3iDmjiHJKN6W5Hq5r3YfG++purFiR2CwR88K5MrwvROvVC7cZyY2DaVg9imvz+x6Jio1KI2m6ITZiKyiG2mxVdkTp/UiG0txYnHAmlAFuJEcJzsOPGpOj7xmsfK4jvhEycCg1M1cVFWM7nTdSuBak5nZdGh1l/fJy6XXT93umYxTSGGRUltNPeCXF+siSevH3AdxaLmIqn06G2V78PITvdUd4Yo25SsJlFX1YPncdeY6UAfCfE0gpTwccrHyRhiluyUP167BVuPmkoK8XCSkF+O6aXIWdV2jZWDPbBzrsAiYptkTjdp4pFmndOcXnAIFbSHnd6IT1yYn3yfjXHiaclXstpQKhCKbvYkJgZkRfO9p2niU1U/WtkXHEIZwQYoSihOCmMXCCanIGNbvIhLC2GR34HuExfpYcfLXpSRrtjFgT+diDRxw3MzCfaI2DZVQyGMz1dDzHRNvFFzeqxxuQ4lBErSJ+4nftMJZVk+cdehxMJRFq6irJY1cV/v1+o9KJprYl5NvgdViDeoiaeYo4XGK0LylPqkRDAApDS2qiuEKGl91BcKafcRaeJSnLjoW6WCY/aJ95gm3jezRpYmrpu6XQOpwDQOxis1jE1VE371oTDxQD1i2/NXL8Uv3/c8XPDsQ4O2NSA5xX0MGlI+qvfSmDndbac5va4Qj/9WQg3X0RZUaZn26raBRahOdpy4bk6vlzt9supFK/voWfmqwBWRCbvHK4nrxaTsOjKxLY2dLmtB6m9VWRMX5vTZkrEtwyduOjZXyo4oCzaZKc2Sv7deX6uFGxvF7XFVTVwz7SbY6bLANvjJPd/PZqeLOgzjpFxTzelB/29AiGtWAr1fx6Fsmjk9R5y4vOiq6xPXCch+Cjs9Mqcn5xuZnS63SX13sRVFv1e9nPh7fB/pPvFG2Onq/dSsT9yMyLeqTc6+n9RggkGiXm8KQZiq+hidqiW0pIGCA85BbAOAEw9aEIWfNbJiFqfKG7WYNLHGzOnosk+co79CEy84qhZTL/lKGpgDbbcesS2RtSuqz3zNVNWLVvZyik353TEDX7/pUXzkx/clrq95HGV9K7pBP8vDTtffR0X2iQtz+iwhtsVJcpK+SlM+gcGiG02ysYnZVwRJxfNTzZ06dKuQTI4S2R6NmrhRYCeP1dfEEaYKTdZR9RiTUexym9npHGdNK9d8Y7KZ6LthzDViTm+InS7M6Wk+8YEMn3jBSRBUAXURnzCnS/1mQnrWgMZOL7rKQj9VE5f97+HiyGriBqQle/GYE75Ex1F94sJ8JUOkuxwr1xIdzg1X4mqcePpLkQdsXogOJ7fTZIqXO1GeMmU2ZbNohp0esUy1LWBjzZgTi6V6bRBx4lntEAMyJrap33VMVf1oZS8TYvQJdfPuidR2iUQ0RAG5rZk4cdHOiXItEiqzRYgX3GCxl+WrlFEqOLEQVzRxVQPK6wbStUiZ2CYSDKmmZo3YpmVFC46lsNN9P6G5iYyS8mFZEAgCX6nohObd/Atg+ZlWPT8pxKWFiGdYRKR9F9eYPpugvwNZqKqm7hzs9FIWO901EtuyNXHJJ26IE/c5WOwH7UJYdr7c6eVasJi0mrgBaatsmZEtIOIwBUzZgMRqd3SqlvjNcQLmaFbaVRlyOElexEI8+7woxCyHABSx1eG33G3RkWYC3zE6hbFyvJOUYk4PMxzpWoxsOWlEo/B8sRVpHXa6JrSrdSbyctWLthcU/YaZlQnVZ8ZouWq6PPKJi/cyWHQxlaKV+JJGpi9gYnO6N+tCzIBg0szyVernikiQOOwqWDgKF1S56kuLz2wBI+oR18p+1WDugFHAmfzfsiYuypMncV0TD6w3ZKgjbrMIpWvOJ66ON31xKrdFXnxWtWdmGnNioVJ0qWGfuCnLHRCM21LRNSd70eLETSFqMp8h3SeuLxbj74nc6SIJU0Uzp2e4a3SfuOdz1/YSB/pIiKetsuWJUkDfT9wkxMcjIV5NaEkOBYOqXtrV6DcpnCQvxP3ooXA65BV9PQQm6NY1cVnDlnHBt+7A537zsPE84dvSSYZpn+uBObCIFJ065vRqvDr3fK5LbAt84oKJGgpYJE3fo9K2lzJkTRwIcgVkmdPTLCPCQvHkvkn8y1XrAMweTRwIzOJZvkr93FgTdxUG8YKQuS6SbAA5NPGwHnGtyNgGxOFfpvzosik6+k0S7KI8WRMXIYkCIq96YrErfY408UIT7HTNJ67zMuS2CEUGSBJhfTYkXwnve8FQsX7udG2+UnY51MK/hMabZk7X2enyXKj4xKUq82riEwZ2ujiuJHtJcQck6/JCd0r3xnLfzBqpD5CToTmBOT3+bmI+CjOKMKfLMtp1KGGONaVdlc8HgMtvfSzzHn5x71bsmwjIUqK/ukQ489hlqdeITrR3vIrP//bhTJKXz4hWgHnY7GlI82PvHC1jx+hUXJ/wiYeuB0FsS/NNNcJQ95jhOMFzz9JEKl48mbzlG7fhugd3hPeQvEYQmoZ1YhsnWbwjKUJcztgGoK45XQjmZIhZfM29T+wDkO2ymU4Q0beIaAcR3d+uMksFJ9NXqZ8b+8QDE3PVD7aQFMx12ZxeT3MVk7i4VjGnGxKxJDVxs09clCfHiXs+K32x4FCUjMpUBxCH0rWFna5ZmOQFwWQdQZwWtz1/sFh3L4SkJm6+V9kNl2ZOH9I18YQ5PficFjKXEOJSPxs3xImL44I1H5Sd4RNXtH7f+sTTEGUa0oltRnO6pom7TkKoCTOKz4Egl1dODlEi9jcrEYeo/5LfbcDWfZPGc/aMV/C+H96NH93xhFJ/seDg0retwfqLzjNeJzrDv1y1Dpf8bgNuXL8ztR3Bs0hu4NAoxLUJtmjV0zYOif+KOHE913na57ptkHziaVp10Kb4tz88ujsaUKZrxDsd1szpsk+8GG44MTqZYk73WNGwS0U3UxMXk3nCnO4xzj/lQBy7Yl50rIc18csAnNvOAgMh3oA5vRCb012Homc+T2jiVd/oMzWhrF0rC4KYOR6fr2v4shYq+8nnRZp4zE4XFiKBgktxshepmaomHiwgB5vKnZ6uicsmYSCeg9KQZIsH588bKjbsE68oIXiSgA3jxImSi5VyzceA60TjwuRSFfuJA+aQPbndpu9iISP7xMVx2VQvL8YsO71JpPrEpVVPlIlL84mLVe8nrlqHt37zNgBxaAEQrHzFZFt0ybj6zTKPyHXtGisbzxGa/5a9gZB/ZPso5gy4WDl/MGQ6m8sfCMOOhC96wHWN5wGhJp7ig20Epn1/gaCDyiY4kzldD8+RfW2N5E8XgrKQM05ch0mIi7YPDbhKbmPZJ15wHPjMGZq4utIeLCY1SgERty/uR2/fQYuG8IqTD4yO9apPnJlvBLCnnWWWCm5C806zaJQKbuSzFCZmgflh+NmkFDJVXxP3lGtl36fjBP88P6nZi7/yGJDDzuZJoVDyRio1RYgHWqdDMDLgAVkTdyOroDyemZNha6KMqaraNiUToWYVmKiY+7hAzQ8iMUQZ4n3NHyxgKnzenm8mrCY0cYP1AgjN6SJO3FCfcNHJ1+khZuLd1aT2KBaHSvwufe0ZAUE0kihDtoal5k7P0MQnp4Gd3j/JXlIeoM8cEUWKroOq5yUythUcws7RKi4Lzd3MHAX5izKHw/OFz0pHlia+eU/MZE4T4qLjPBlq6g9uG8ExK+bVTQ5zyOI5GHCdyBSdNfBYskq0oIhHAq2mTQAVz1dIXEIwMYer5kJg6lTiX6VJo2qYeNIg0q4Wwv19xQ5GOtLMeiatX7yDwTAuWNlmMLyXghtMJqNTZk1ckIXEexssZJnTJZ+4bNb0413L5kpb0c6WOHEgmHzLNQ/MjFd88Wa88ZmH4J9+ZrbWF12KrCdDA64yPhcOB8mSXvOVW6NjX7thI+56fC8OWjQEhwif+LMTUPV8vOySm/AP5x4b5WZYFF47POBCrNHFgvpn92zFz+/dis+//mQpvJVx7QPb8R/XPhLVJbO9xZaphx4wLPnEVXb6YNHBQMGJmOOv/vItePPph6b4xINzN++ZwAeuuAeXvOlUMDNe9Lnr8djuCXzjbWtw1vHLAQDXPbgd7/zunYmxp7uJ5N//8rK1xuctsH+yimd/+jqMlmv4xXufFy2YFw4PYPd4BUd+7GoAwFufdSg+df6JAAKO0XlfuCkxD8q8lg/9z714wzMPwUfOOSbMLxGM8fGKhyM+djXmlgr4w4UvjhPBaJwj+R6GBwrR4uL0T18HAHjeUUvw9uccFp3zrVs24bJbN+GTf3YCLr7moYgPJVCS3LFy35Jzp+shcl/4v/XYuHMMl7zpVKU9n/rlAyg4hBMOnJ/5bNuJ3lz6G5C2yvb8OM2q0GYpjPUU0IXyWLmmaOJAbMosugSDrMhcWZ1/6kE4ePEQAGDXWDJBCBCbbbbumwQz48FtIzh2Zf0XPVBwcLzUIeR2P7ZrHOu27o++iwWN+NwsTINFDGBVE4+vmazEvi11IknXxJ/YM4Ef3r7Z2AY/fK8R4zdlAZCmBZv87+LcQIjLVotYKyq6DuRd0XSI1X7sE08ntnm+7N6I2y+HlMn7yfewOb0uiOidRLSWiNbu3Jnu8hEQ5vRyzcf9T47g+ocDLsOaQxfh8CVzAADveuGRuOwvngkiwofPORZ/d/bR+OBZqxVN/NlHHIDnHbVEKXtkqoZ7ntiHe7fsx/1PBuNjbKqGR7aP4aGnRqN+8Lo1q/CVNz8dB8wtKcle/vXPTsTfnX00HCKs3z6mzD0bdowBAP7+7KOjYwBCkt0AvnnBGlz2F6dF2RcrWmKZvz/7GHzq/BNRDBfmd23ehwe3jaSw012864VHAAAefmoUQND/HgvDHzfuHIuu2bAjaOf7X3wUPnzOMVHbdG2/5jOOXj43l5B5av8URkML4GO7x6Pn9s7nH4EPveRo/N3ZR+PgxUN4ePtodM32kSls2TuZGBPyYnvvRBX3P7lf2ujFVebcsXINO0fLUQx5UhP3sXjOAP7rDafgjGOW4twTV+Afzg36x8mrFuChp0ajZ/6hlxyN95+5Gj4DNzyyE+MVD695+ip86vwTo4WhcNUAKoG5VIxN+VVpMVbzGX96cj/u3bIv+B4uUD5w5urod+sTNyA12YsUsyseuOtASbuqT457xiuKMAIgmdOdRMpWUxkyjlw6F7/54AsAZGniQbuf3DuJrfunMDJVw3E5hDgAnHLwwujzmGRBOONz1+Nll9wcffc5vo9mfeLyHsMmv1K55qPq+XjfD+9WCIEBQSWZnEImtOgs8/9Z+wQu/OmfIs1DRmCSitmiumZ9zxP78Jt1T2UKWx3iHQyGZn+TT9zEBpbHo5gIY3N6uiYuh/zJayph6h9wnYgpL773K5j5UmZew8xrli5dWvf8UsFVFktCcP3NGUdGk+vKBYM445iA9HnKwQvx/jNX44QDF0R9Agiyub3rhUcmyg9ykMfZz+S/QnNbOq+E8562EgAUIf681Uvw/jNXYyhk0Nck4REvAA4Oj8V+8lLBwZnHLceSuSVpkxc1Vvvo5fPw9EMWoVRwMBa6bERssYAcJ37Usnn4s5MPjAReWqy2+Pz+M1fjnBOWR22Tu7IwCc8bLOLdZxyVeGY6ZKuf/K6OWjYX733xarz/zNU4YslcpR1pC1p9UV2ueXGO+EJyzhX1ydqwrAkXHML5px6EUsHFwuEB/M0ZR+L9Z67G0w9dBHlzl5ecsAIfDIWr6GNveObBeOuzDo3Gv6qJx59LBTdOrSz1gyBHQdx+Udefn35IdK1lp2vwfXnz+uDDU/unsPaxPaj5HO4fHmttesY2fVW0e7yS0MTFQy+4yc1TxPEsDA8UMDzgYtdotiY+Wq7hzsf3AoBCasqCLMTHy0lzuhhsWWlXH905hmv+tK1uXfKEI38Wg3Oy4uGp/VP4xb1bAcSs/clKQGzThbi86NL91LvCtKY7RpILH2FVEO9Uj2P96vUb8a+/eCBViOtaCBC/g6EBN8rMJeoSp5oWa3JOe+EnVIR4Wpy4RGxT88lzVO5M0cQbhWCnC6Goh1UB6fsxK9pSwYn85TLKNQ8VKQ+5/DdOaxovoOTc6XobPWXyDuLBRZ+QNXG5HfImL/KiUMztpYIb3bOoQ9yX8ImLRV0QUx+2X059WlM/F8K8CkrKVy10sub7cB1SBJf+DATGy2r5QmjJ40Fum7hfE3S3l5Ijvph0YYpF2ECYelbcD6DG4esQi0Oh9IlwvgHXUdwUQDx3ye8t0bekHffk+uX2i/c7LC3IrSauwRRr/NUbNuLt374De8bLWDxcjPaeBsJYT80nDsSdb89YRfGJA/G1BcdRiGpxGfUfldgIxQRZ878rFOIrFwzWLRMAXnzcMrw5XOWZhLgw8clEKt0p/p0/PI4PG9KI6jDljAbiCWOq5imac1Fjc6btkwwktem9kRAPwtYuvXFjFKYn/M5R0gzt2j0TFeydqGSyZHXBL/vExUYXgEoUMi3Wio4jtcNX48QzzOk+c5RKVzGnR8lGVCHewyFmPwTwBwDHENEWInpHq2WWikGyFyEYIjKXNKGbLGJAUtCaBFLVY0yaNPGqusGIQKyJS20ME9LImy9Fu25J2iFzvOuVfG1Qn5qjXNRTKjrRPQtBL4TA6FQt2J3Lic/V70Pci/xZD5Oq+UlzulgsmBY++tNWNPHwuRU1JadUdBUBnde9FZQnm9M1IR6+Ozn1qawJpyVTKYU7KYp3Jmvbch8D4n4kv7dE3xJ5/quyJu4rlglhZpejHLo5lvtCiEcTrKTlbR8Jsodt2DGGFQsGo72nASj+TiD+fPgBga9tT6iJywNWvGyRUUlHnr2el8wdwO5xsxCXO7fwpSyR9hLPwvzBIi561dMwt1RQMqYJCH9Zlk98olLDRCWZJ16HPOH4ihAPOupU1VMSochxlEITlwW3rH3rQlVsMLI9jD3/yZ1P4qd3Pxm1wyWKVst6zPre8QomKl4qAQ1ICv7YnB744GRTN7M66GW4bqy5BBq+NBmHxLb9k1XctXmvcp3P6vPRn0nRJZXY1qOaODO/iZlXMnORmVcx8zdbLVM3p4s+NVhwowVzmjKjM4jlSVjG2FQtEnTR35osPGQhLv6qAkqPPxdxzbGf1pd8u5KGKm1TqwpxROeKexaxxWJBNzJZVcsquGYhrn0WWe3kBFF6EqNo855i8pnp855MAIvN2+p1eqhgfk3cUxZTZnN68lkDgv1vHiviuU+EC3ZXWgiJ5y3nSRf1CxS0BYpsUZHZ6bKrpibJpzhpjBXiCmp+PEgEU3l3SCCbqvpYNm8wMpkAQWeU/YvixQjCzK7xMibKHhbPiYVoTGyL8+XKyEr2IrBkbindnC4NiHVbRzCvVDAOpCzMKbmRJi4L4/WSJi6zKW98ZGd03kTFi3Yby4JiClfM6aEmXvWjFS2gPheTTzwrTnxPKMSf2h8sfHaPV7BrNPgsBGCaJr53It0UL6Cb7ycjTdyJwtdEXcICYRKkBWnSq/kiYxvCsoJQqctueQxv+NoflMUap8SJi3cwUHAUE1w/+8QbhW5OFwQqWRNPM6ermqBZExdlGs3p0n7dApEmbjCny/HnJj+taVEgJvNJzW1H0uJP3LOoI9LEyzXNNO8k7sP02aSJq8leYn9yHnP6RFn2iXtKHaa2AekJe/SxWK750aJazlEu/5541lKYcZqQFM9dtF0Ie/l56zHhWex08duUFsJYCYV6LYwycCjoO3r61m6gL2YN8fDEROr5rGi8gSYumdOJFFO1EDTL55cwVHSxZyzQxOcPFjAvjBWNzOmSGUtGnpdyQJY5XfIbVWo+ls7Lp4XLmFMqRH6qUWmA/eTOLbj0xo2YrHpRJ7x142687Vu34+4wG5gshLOQ5hOXV9gyA19ncybY6Rlx4sKcvn1kCr7P2DtRwa6xchi3zSCK350sxJkZeyeChcRTI1NIg06k00PMXMnUHcWJm7aDlSY9MWhldnrF87F5zwSqHmPnaPz+ZbP7Z3/9cLSpikxsUzXx3jSndwIyO10/Lt5Bmjk9j09cwEhsM5jThXCV6xRt1H3iYic1INgvwbQoEGWPa2672IIT1y3qkF0rqmneRTVM36qY0HVNXBNOIhxSQCQscR0HJuuFbtJOaOKSyV5um9qmFHO60SceP7eEOT1cbJUKrpL9DkBmMpXouVc0TVzz4wOSJi69N9XKEwtxOTGOvHATrpB4saC+g26gL4S4mIxlbWiPtNfzivm6OZ0UU7V4nvOHilg8ZyBipw+X3EjYi5cw4JIykGMze/1HtXReCXsmKka2shDiB4b15TWly5gnmdP3h0LszacfgtXL5+LTVz8EINgLHQC27Q/i0YXFYjIS4vmzNClCXBqocupVeTANuEliWzUlTlwIbVHeyFQ1miTHyrVAAKYQ20am4p3nthuEuHh9idW/LsQFfYDjyc41mOkKTrwBh2Dbi/NEn9y8Zzy8l1iI+378fMbKNVx1b+AqSPOJd3PgTzdKoQUjuaOYK2ni5mt1BnGaOR2QhbcghgXCwyG170bmdKlOIaAUdrrE/QiOsXFREE/+qiYemdOLqhCXfeKJssJzK9qiJ+kTF2bieOGr5k7nbE1c+65o4tV4e04ZAznN6XrER7nqReNRfp5qfUGcuL5LZJYmLnhPsSau8qHkz2ZNXO1bhXBOk9+j56t72cvt0d9BN9AXyV7EhC2HbeyTTLrL5wfmdEEiEvv1Csg5fw+YGyQqKNc8DA8UsGSui0e2j0UvW7w0gTmlAkamqrkm2JNXLQBzQFx7jha7KjTgw5fOwdb9Uy1o4kFnEgLwjGOW4aJXPQ33bdmH9dvHsHLhIH77wPZokSN8xsKc36wQl6+TBZViTg9NofKqu+oFSWAqNT/auQsA9k3GG89sHylH/nEAOP9Lt2Ci4oXmdJWZCiDKPw8E5nRRvoCYt2Tz29dvehRf+L/1AALtWeZN+FKIWdHwnsXe4USxAI7M6WGfFAl/dkiLCjkNrnxO5BMPk3kI1NsMZyZBmGF1Zr/sAzVxU4B0BrEJwuQZ+8S9yLcrP285xCwqu+hgXNqqWGy6IjK86Vui6pr9QMFJJBaRuRQC5VqQoGqwGHA1mJMM8LjtGeZ0QdiSNFeF48LBQsR1zcS2bJ94ujm94oU7+zn1dzeLy1O5CdkhZho73c9ip2uauNiprqhaNgAktOfgmPr+xe/juiYu9SeZLW818RQIk6zo2DtHy0rcrTCni4xXic4YCr65gwUcMGcAO0fL2D9Zw/CAiwMXBklaRHEFh5TBPbdUiBYH9XDa4YvhOoRbN+5O/DZV9TBYdLBq4TCAgATXKOZImvi+UBNfOBzkaz5p1UK85hmrsDTU8IUGLs4X5qBGNj1I20RgpyLEDZq4XIYXaxhyeWKRMeA62D4ypVhWNu4MtFoxUQKqVi2fW/H8KN2lDuF7/vpNj+Liax6SfOIu3n3GkXj101cBMLPTl88v4eUnrcRQaDp1Q4a6EOJRxrZwctge+ubFAkfEk8vPRwjxSo2je5+tKBUc+IxElIjCTm/RJy4gCwwRuqQLMVGkyZwuZ4sU+RBEO+QJ3UT6mijrPvH4t6h9VbEHtbzRS1LoyDHuJkKZkZ2eyJ0uNHETsU39Hu/w5ShCVblHYSUQcex15heBms9RmK+e7CW4n3ix1YgmHvnEK6omLtpNJKXnzsFOF3/l9yjY6cH9Bklg4sWC9YkbofvEZXPuQMHBouGiYk4XQlhou0KAzS0V8IxDF+GBbSN4cNsInnnYYhwUCnEh6INkL3HdwwNu7nCBeYNFnLxqAW7ZuCvx22TFw1AxXjQ0Y06fWypEHV9YIhaFQlxADAbhmxeMzLw+cdn3LLuUUzVxjc2ZMKf7PobD9yab1oUgPnblPDy1f8poFt+0a0zZo1lg74RKHhS8Bh0f/vG9GC/XcMsG9X0UXQdvf+7heO6RgbUk0MTj3wDgkMXD+OKfPx0PfurciGBXcGISj5js5w+pz//aB7bjN+uewjsuvwMbdowpzyfhE59FaVZ1iIlTT/Qz4MZxwfk0cbcu6VT2g6cJozh3ulq22SceT9SelABGXxiUCm5OTdxPCFeTaV4Oy5qvbUIiM8cVdrpMbPNjdrqZ2KYrP55Sl2yyl+9RtE20Iy8itri0cIv2h68lXRfyft712Omi7brJXM6TX9AEL5DsW+Kv/B4rNVaS78iWgZ5kp5u2ISSixUR0LRGtD/8ukn67kIg2ENHDRHROOxopVmAlSRMHgon28APmBGlWKTa9ijEt/N37Q4E3b7CAC55zGBYOFzF/sIC3PfvQ6ByxMYme7GVOqdDQqur5q5fi3if2JQhuk9VAiB+0KBDizZnT3ahzCpPygiFVoxcdVGzgISZJoYXqmep0qEQYsya+SxHi6mRTcJLs9MFQE//gj+7Bhh1BOJwQ4s8+8gDUfMYdm/Yk2uL5MIaY7R1XJ/65KUL87s378LuHduCeJ/bhaQctSPxOTnzPsU88eH4HLx6OznMdigS5mHhEu048UC33pvW78Hc/uidaOMgT49b9U/jEVeuidJm9GlLWDYiJU4l0cEhxZ6UJcZO2lAU5wUsUf6wJo8gnLlVpZKdXVQKZwk53dSHuGHziScEh6pBDGU0+cdmcPn+wkMpOFwtMEUkh4HMdn3iKJj5/sBA+t6QFQzb1B3/zC3F9y1UAEUck9sE35hPXuQgyO13+C+RjpwPB85ff42RVZe2rPvFwgddjceKXIbkN4UcBXMfMqwFcF34HER0P4I0ATgiv+TIRNRZHZYDQwoSJSYQV/cfrT8blf3kagICQImdsAwLCGxBr4vMGC5g3WMRX3vwMfOUtz8C8wWJ0zpgUkuAoQtxtaLI998QV8DnIKPZoOFmPlWsYm6phsOji0AMC4bA8Z6IXGSZz+gJNE9QnPiF0Ip94CntUIC05izxhyJYQmVE9UAh3MdPixIekqIJv3/IYgJiQJrThmzVt+XvvOB2fefXTIh+1SRMXKUvnldRnIOOHt2/G6FQNr3n6QYnfRMv/7VcPRkJXsN4PkYR4wQ0EuOsSfvdQkOP76YcsBIAoZ76M8YoXWR3W7xhVfrvs1sfwyV88AED1e842iMlO1sR1c3Da2llnENeDTKATWdz066L9xDVTfbnqw/NkTTw2p4sFa8yyTpqak+z0ZLsVTbyoCh35s+wWmD9UTBLbwmuDrU7Tc6e7jmO0XuiPe7ziwaFg3kk1p0s8peBvPnM6oGbpE4vdwdB8PiWZ02W2PRCy01OTvbhR24EkkdDk/1ZzpycXT6WC+h7VTHZ9wE5P2YbwlQAuDz9fDuB86fgVzFxm5k0ANgA4rdVGipWwIBE9uisQjocdMAcrQmH4rhceiTc8M8hnLAbku84IcioLn+nccLJ/9pEH4Lkh8UwXpgMFwgmSdjWvVMw1UQgcu2Ielswt4Rs3b8Jbv3k7mBkv/cJN+PW6pzBYdLHm0EX46luegResrp9fWsfcgUJAEPN87JuoYs6AmxAEet+JhLgwp9fTxH115S4gm+Fls7i+ck3sYuZzJMQB4DfrnoLnM27esAsHLRzCs488AAOug407x5Vwq+etXoKl80oKC3i8HCSr2bZ/CgWHUAzvfdn8dKuG4CfoRENAXfD8V0h6ezK0yIjFVnCPgaZw9PJ5GCvXMFR0cfKqhQDqk9Ge2GPeXx6Y5Zq4MKdPxlqNnqwk3SceRpJIptEsJMzpVd/gEzcQ2wqm3OmyJu5EiT/ke5Kv1zVxMpnTQwa8K5vTtTjxoO0xo3v+YLo5HUC0D7nPHC20ZU1cRqT8aMcnyjWUCi4GC67io1buMcotHrP/daSFTkYbvRRjYlshtEaMTcX+eJltD+TUxMs1hd8ka9UC9TXxWHuX36OeyW662enN1rScmbcBQPh3WXj8IABPSOdtCY+1BCEUxCT/w9ufwFDRVfzBbz79ULzw6KU4+/jlePqhgXX/6YcswmMXvwwvPi5onolMtio0b7/0aSsAACccuACnHLww8lm/98VH4aJXPS13W4kI//KK4wEE244+uW8yIjQF+1gTzj1xRVMrNWFqGi/XsHOsjMWG+0lq4lVUPT8SvHk18aK2j3faCptIFeJCE7/+4R341X3bUPN8DEmhM7vGKvjeHx/HLRt24UXHLkXRdbB6+VwAwOI5A/jz0w/BR887NjpfaAwXfOt2nPAvv8Hqf7wG37x5E045eGGkIf3tWUcrbZKz7wHBzlhHLp2baLv8Doou4ejlc3FqqGEfvTzOa3/AnAEcMHcAn3jFCQCAlz5tpTLhvfjYoH8tC10kx66YhyOXBomFzjlhOV73jFX45gVr8OpTD8LvP3QGjg83vhkqtmyk6ltkauIpQkXAlGkrC4o5PUUYiW7saP3ZlDtd8Yl7Zna6+J5PE48ZzkZzeqpP3GxOBxBxUzyOBYrPQZ4Dfe5J4yBMVALzuUj7msZOB+JFvsmcnibQ5FzmYj074AbhnMpviTjxdHb6oMjYVvFSBbLeLsUnbrDyZGvinuoT11K6dgPtDjEztdyY55OI3gngnQBwyCGHZBa6atEwPv/6k3HywQtx6qGLUHIdnHb44oRJiIjw9betSVx/8atPwl897wgcYCCTlQouHrv4ZQCAOx/fi1PDzUau//AZeGjbCE40+FLr4RUnH4jBoov/7ztr8cv74k1HWp20hca5cec4ntw7ETHdszA6VVNIaZOVYJD99XfX4viVC/CBs1Yr54uBMuA6qcleZMihXYJJWvV9vP3bd0THDwnT3QLAs45YjH+5ah0A4Mxjg92WnnnYYqzbOoLTDl+MT2sLpuNXzscbn3kwhgcKWDqvhNGpKuYOFvD25xyGDTvGMFHxcNiSOfjOX56GI5bOwZ+27MdzjlyC6x/ZgRcfuwy/e2gHzjtxJVyH8Lu/f6GSMnaw6OIbb1uDZxy6CINFF0MDLsbKNdzx2B7FGvP5158MAmHBcBG//dsXRGRIgf9+06m48/G9OHblPGzdN4W5JReVGuPIZXNQlNwzZx4X3O+P/+bZuGn9rojkeMc/nlWXqzDTYPKJ66bIej7xrPhwGQo7vdqAOT1kZQ+4sk/cwE6PhLiBnZ7DJ+5zIAxcmZ0um9OLSXP6vMFCKjsdCHfj8wJzesEloKrmTpchfk+a02sQIXy7x8wWjDzmdFG+DnmjF9I08TjPeXPs9PFKLRGKKP8FGmCnaz5xfXc3z/cTWn032enNCvHtRLSSmbcR0UoAO8LjWwAcLJ23CsBWUwHMfCmASwFgzZo1RkEvsHjOQBQOZNKo6mFowMXTVtUXxs84NOLnYW6pgDWHLW64LoETDwq0rZ/fE9/+YNG8Is2L569eioJD+O26p7Bl7yReeHTSJK8T5kanqkpYmRDodz6+L5HRDJCEeMFR/NBp8eXyoC0VAoElTMgisc7hBwzjRgTbJL77jCPxk7u2YN9EFc9fHZi4/+llx+FvzzoaC4aTvu2hARcXv+YkY90nhSZtAHhB+CxWLQoWNq885SDlLwAcYeg7Zx2/XPk+t1TAi45ZphxbOBxbPGQNXWBOqRDVv2xefa7D8EAB55ywIvreDMmx3xGz02uJYxFBtQ47XZ6QHUrffjfwiavs9LlaWKJQFhViW9EFc5D/HAjit6fkjUbcOuz0Yj52OhD4cOuy00NiW8EhzBlwDT5xSRiFbfOZo1DGKHe6Zt4e0KJ6BCbKHuYPFkO3Qoo5vaCZ0w2L/bRQypGpWsQWj8zp4UJmRDKnJ9jpfgY7PTKne8p8m+0TTx6TXTU6O12Pnw+y4E0fO71ZIX4VgAsAXBz+/bl0/AdE9HkABwJYDeD2VhvZj1gxfxBL5g7gwW0j0tHWXuyCoSKec9QS/PyerdgxWo4ElozBohtpoECoiVfigTVZ9cDM2D9ZUZKmCAg260DBUXziYqUvBqlIsCIP2gEpUceA6+CGD5+B3WMVLJozgI+/4gQ4FEwUbz79UKXOgutgwfDs9Q/PRkTmdFkT10yR6RnbktpswXFS9wWQzekVz8dkxUNpgc5OD5L56O4hQLU2VTwDO72abk7X042KiAj93Eot3CJUCBujTzxOfSoEq8jJr1sXCmHbPClXQZpPXPyur5nEvcZx4gZzejFuG2D2iaeR0OSNXqIQs4KDku+qzHWJbQ/k08Qrnq9ErTTKTtcXUUof0LLmyc94OjTxujNnyjaEFwM4m4jWAzg7/A5mXgfgSgAPAPg1gPcw8+yyE4YgIrzwaFWj2z9p3hylEZx9/PIoX7jw5+v42EuPw/L5JTzvqCUYnaphQg6JqHqYrAbs6f2TSRuXrIlv2z+Fz/46SOc6VfUUJrywAsgdesncEk4+eCGWzy/hb88+GvMGizhsyRwsGCrCdSgXCclidsAUJ55kp6dp4kmTc5bmo+doH52qJbRmWRvU2xO0RTpeTGGnG8zpOky50+P7kn3iBnO6lPpUJMuphYxzn5PCSOROj7MeImKnq/UKn3iiSZFPfKpah51eTZrTRXlZPnFxb+LRF8OFjMxcl9n2QJ3c6QbimtxOk/9byZ1u6Ftp/SC4X195psVoUdJDxDbTNoTMvJuZz2Tm1eHfPdL5FzHzkcx8DDNf09nm9zb+4dxjlO9yprFm8TyJZX1QihA/atlc3Paxs/D0QxdhrFxTiBiTVS8KT9s/WUtcK/vEAeDL12/E/slqYIKUVrYvCE3hUajYUQdgTqmAN59+KG772Fn4mzAywMLCBDGZ7ld84nnZ6elmUBP0dKUjU9WEMHIoSaQzhXnJ9Ubs9LBsPVIkKyuanoMcEFve1jenlwqOogGbFhExOz1ul9iASH9W4nfToknkpg8iQ5LtzjKni9/SQilHJmsJTTz2icfMdfl+gCD9bWrudMPulUFbTD7xLD+5/PzN/UDct2zd6KYGLtAXudP7FcvmD+KGD5+BPeMVvOrLt0YxyK3gMCn0KU0TF5gfCl058cxU1Y8mzv2TFTCzoiHHmnjcWe/bsg/lqo/B8NhJqxZETPmTD16IAdfBZ197ciu3ZTHLoLOa5WN1NXE3OdHqfl4Z+kYrExUvMRk7RAktq6T5VCelDTtEO8W2lEUtSZR+vVyP3vbovuqx06VY7ThTmsRJ0bRQncgm3A1JdnpoTk+0CJE5Xb93U9vkv3K5aYJtshqb58VzKboOCg4S9cmZIPV0xjLEltQVjYUvFh/qQid9Mai/e/mzzDEKwgP9xMJTTpTVaVgh3mEcKsWyn5yDXFcPRISl80rYOVqOEtWkQaQjldOkypp41WNMVoONYATEQJG5KHdv3hdla7r9Y2di/lAxSrjy4XOOwXOOTMZgW1hkwaSlxj7xdPNu8LvJDJrfnB7Un9TEs8zpJs1MZqcb78cgqE250wVcx5F84um50wPtWBWeepnCJ+4zR1E8wvWV9ImbiW2izLTnELRT94lL3lNhTs/IhyDuTRbiZPhdsO2BbJ+4aGPF83No4gZ2ups8lnX/IlGPyPmgM+m7ASvEu4BSwcU1H3i+ksqzFfzqfc/D+h1jdXNGHxqGdv387iejY1NVTzFh7p+sGoX4tn1xVra7Nu+NNPFl4cJh5YKhKDTPwqJRmIRYNGFrG1ToMGlQ2T5x1ZwOJLVkIspvTjfkTs+6HxlxiFnytzR2uthBT059Km9PKsSFzhGI2enZmvhACrFNlGlaUOjfxQJBJ/LJ5Zsg7kMQ9ILUu0nhK9j2QDY7XZQ5Ws7hEzcSJBszp1e8wCc+WBR9VvAPuifELSW4Szhu5fxEWEuzWDZ/MMo4l4XTD1+M569egrWP7wUQmNcDIR775l/y+Rtx2S2bou+i84mtQdccugi3bNiFx/eMtxwiZ2EhYDI1Nx4nrrLT0yAnSYnrSuZO1+uTyx80sMXdUDuU86mnXR/Xk25OlzcmkX8nIokhrpnTpWx0+qKm5rGy/0A9TdzoE9d2iWsmd3q2Jh62TdqeV343g9KCKQ87PSgzmTWtVXa63h6BctXsE++mJm5n5RkMIsLfvyQm1y2eM6D4xAFgtFzDJ37xAO55Yh+AuPOJvOSffe1JcIiwfaSM809tOfmehQUAlYAk+lo72elzBgKTs+tQZIaeM2A2kYq6snzaJm1MpBgWjHEdZp+4Wr/cJlUTT2q8sRB3FeEZhbhpzGuRO11kLhRC3NWEalau+rrm9Ax2ukBa2lX5PoWpvOgQTFqw6hNPZ6fLbdLz4Ovtz8qdrlgfUvrBnIEwLa8Xs9P1PO/dgDWnz3DIfvhFcwYUn7iMn961BaccvDCKDb80zHx3xNK5+OxrT8LIVE1JnGJh0QoKrhNpV/MGixiveA3kTk83gwoEZdbghuzxcs2L6gHMQjypiWf7RQN2uteQOV3PnS63KS1OXNQp7mPhUFHRxEU6B3VRI+dOD8qqaux0oiCBjWCPm33ibqY5mYgwIOWPUOLEte19ZUEskvOI5yb2xyi65kWDwk6vo4mL+zGnUFWfkfybaJd+LK0fzAvz1yuauNg6uYvENquJz3AQEf76hUcAABYND2Cq6kV7kQuccvBC3LxhF6aqHv5nbZD6fvn8UmSyf+UpB+Gtz1ITtFjMHhDRueHWwhuI6KPtKldMiPOHCsr3QiTEzdfl8YnPHypEGqvQYEU9QNInTZTURNPMyKpPnGGKn9av1yH83HKb0tjpok7Bspd94iKNrH6NaJvHEju9pvrEi476vM0+cUfTRM33KZvT9XIKhkVZvNGLSM4S+sTD3On6eQl2epYmri0G1fpUa4V+jIgSW7Wm9YO54XawYhtZuU7rE7doKy487zg8+umXouAQ1m0dwQ9u2xxt8AIAL3vaSjy6cxwf+fF9+P3DOwF0N1mBRe8i3Er4SwDOA3A8gDeFWw63DDHZzh8MkgjJGm5Yt/E6E4M4IcQHi5HvWPjERT1yXQJGc3pqnLjGTq/mZ6cLCD+33CbXjf3BTZnTtedR8/3QnK76nXWtMXrehnbW84nLbfN9RsXzMWdANfDqzG25nEgTFz5xl4z1RWx7n8GcPT+ZsqblZaeL4yYio3zuQMHBoNiq1uQTN6S07hTsTD1L4DiEM49bhmNXBLm/5djxFxy9FETAVfduxRnHLMX7X3wUDm0Tk96i73EagA3M/CgzVwBcgWDL4ZYhfKXzh4QQj8OJgMZyp+txw8OlQqSxbtg5hvFyLapHrkvASGxL0UBlIbF3vIJt+yfNwq3Ohkelgqu0KVMTLzh4bPc49oepSkX7b1q/Czet35lob8Eh7BgpY3SqFgnRDTuCLZwjTTwKLUP4t3Fzumjb+u2j+PGdWwAAwwPqOaYtd8WxWIhzdFy+9wHJFL959zh+fFdQR1qcuFymvCgbjBYN2XHi4ng6HyJ+96WCiyf2TmBkqia5gILfrSZu0RG84ZmH4JoPPB8nr1qAtz/nMADBVpnHrJiH77/jdLzvxUfhv95wCv7uJcekbgNpMeuQa3thInonEa0lorU7d+7MVbDYxe3o5fMwb7CAgxcPRceHii4WGjbEAYKJfcncASVk81WnBhskiU11Dl08jIMXDWPF/EHc+fhelGt+tAUsELiLZKxcOIQDF6p5F2QtWb528ZxgQ5yl80p4amQKj+2ewHJDzgY5j8NZxy1P/H7I4mEcu2JeJEyWzi3hkAOGMVBwoh0Lo/YtGMS6rSPYP1nF8vmDOGDuAAoO4bJbH8Plf3gcBYeidom2PbprHLvHK1ixYBALh4u447G9UT0A8JZnBbtHHhzuwSC2al69bG60re6K+YNKW5YYdoJcuWAQdzy2Fx/5yX0AgJefdCAA4E2nH4JSId5qeMWCQRy0cAgDBSd69+IZie2jX3j0UqxYEPx2wJyBaB5aNr+Ee7fsx0d+fJ9yDyasDPNyLJU2I1oyt4TBohP1MSDoZ8vmlRIbtBy4cCh6JgCUd7t6WXAvBy0cwooFg3hk+xj2jFeiTYxOCjlIYvvrboCYu7diSMOaNWt47dq1092MWYdgRyQnkyRi0VsgojuZObnfbufqex2Ac5j5r8LvbwVwGjO/L+2avON5suJh52gZqxYNRZttEBGYAz/zYIYmW655yjav4poBNyBAORTwqjyfsX1kCg4RVi0awr6JKiqenxC6nkYCE9g3UcFk1cPKBUPYMTKFgutEwrLm+XhyX7Bj34ELh4wa5/aRKQy4DhYMFVHzWUlBKjY9GQv3N1gZCq+pqpe493LNw1P7p0AI7sNxCHvGKxgNc4zPGywqQrzq+dgatm3VomGMlWvYN1FBqeBGyaeYA/N30XGwZe8kls0vwXUCgt9EpYb9k1UctHAIRIQdo1NwiYzbOU9WPOwYDfJKFFwHBy4YjN7fVJiVbcveSSyeM4BSwQEj0LzFuxfvUNw3M+PJfZNYMFTEvHAhNVX1ohTPrkNRu0wQ72XlgiHleYu2iOt8n1H1k64Q8V7keXH7yBSKroNFw8XoXoqug23742fsaveRB+0Yz1aIW1j0EaZBiD8bwCeY+Zzw+4UAwMyfSbvGjmcLi3xox3i25nQLC4ss3AFgNREdTkQDAN6IYMthCwuLHoCNE7ewsEgFM9eI6L0AfgPABfCtcMthCwuLHoAV4hYWFplg5qsBXD3d7bCwsEjCmtMtLCwsLCz6FFaIW1hYWFhY9Cl6gp1ORDsBPJ7j1CUAdnW4Oc2gV9sF9G7berVdQG+37RhmnjfdjchCzvHcy8+4V9vWq+0CbNuawRIAc5h5aSuF9IRPPO9NENHabobX5EWvtgvo3bb1aruA3m/bdLehHvKM515/xr3Ytl5tF2Db1gzCdh3WajnWnG5hYWFhYdGnsELcwsLCwsKiT9FvQvzS6W5ACnq1XUDvtq1X2wXYtnUDvXwfvdq2Xm0XYNvWDNrSrp4gtllYWFhYWFg0jn7TxC0sLCwsLCxC9IQQJ6JziehhItpARB81/E5EdEn4+31E9PS813ahbW8O23QfEd1KRCdLvz1GRH8ionvazSrO0a4ziGh/WPc9RPTxvNd2oW0fltp1PxF5RLQ4/K2Tz+xbRLSDiO5P+X06+1m9tk1LP2sUdix3rG3TMp7tWO5I29rbz5h5Wv8hyMe8EcARAAYA3AvgeO2clwK4BgABeBaA2/Je24W2PQfAovDzeaJt4ffHACyZpmd2BoBfNnNtp9umnf8KAL/r9DMLy34BgKcDuD/l92npZznb1vV+1qF+acdyc23r+ni2Y7ljbWtrP+sFTfw0ABuY+VFmrgC4AsArtXNeCeA7HOCPABYS0cqc13a0bcx8KzPvDb/+EcCqNtbfdLs6dG0nyn8TgB+2sf5UMPONAPZknDJd/axu26apnzUKO5Y71LYOXdvusu1YztG2dvezXhDiBwF4Qvq+JTyW55w813a6bTLegWD1J8AAfktEdxLRO6ehXc8monuJ6BoiOqHBazvdNhDRMIBzAfxEOtypZ5YH09XPGkW3+lmjsGO5s23r9ni2Y7nzaLmf9ULGNjIc0ynzaefkubYV5C6fiF6E4IU8Tzr8XGbeSkTLAFxLRA+Fq7RutOsuAIcy8xgRvRTAzwCsznltp9sm8AoAtzCzvGrt1DPLg+nqZ7nR5X7WcPMMx+xYbk/bpmM827HcQbSrn/WCJr4FwMHS91UAtuY8J8+1nW4biOgkAN8A8Epm3i2OM/PW8O8OAP+LwJTTlXYx8wgzj4WfrwZQJKIlea7tdNskvBGa+a2DzywPpquf5cI09LNGYcdyh9o2TePZjuUOoa39rBEHeif+IbAGPArgcMREgxO0c14GlaRwe95ru9C2QwBsAPAc7fgcAPOkz7cCOLeL7VqBOA/AaQA2h89v2p9ZeN4CBH6jOd14ZlIdhyGdcDIt/Sxn27rezzrUL+1Ybq5tXR/Pdix3rG1t7WdtbXgLN/xSAI8gYA3+Y3jsXQDeFX4mAF8Kf/8TgDVZ13a5bd8AsBfAPeG/teHxI8IOci+Ade1uW452vTes914E5InnZF3bzbaF398O4Artuk4/sx8C2AagimBF/o4e6mf12jYt/awD/dKO5ebaNi3j2Y7ljrStrf3MZmyzsLCwsLDoU/SCT9zCwsLCwsKiCVghbmFhYWFh0aewQtzCwsLCwqJPYYW4hYWFhYVFn8IKcQsLCwsLiz6FFeIWFhYWFhZ9CivE+xhEtJCI3h1+PpCIftymci8jok1E9K52lBeW+f+I6Cki+lC7yrSwmEmw49miGfRC7nSL5rEQwLsBfJmDdH2vbWPZH2bmtkwiAMDMHyai8XaVZ2ExA7EQdjxbNAirifc3LgZwZLiB/P+ITeiJ6O1E9DMi+kW4An8vEf0dEd1NRH8kosXheUcS0a/DHXNuIqJjTZUQ0euI6P5wB6Ubw2NuuBq/I9zc/q+l8z8Sbmx/LxFd3IXnYGExE2DHs0XDsJp4f+OjAE5k5lOI6DAAv5R+OxHAqQAGEeTp/QdmPpWI/hPA2wD8F4BLEaQCXE9EpwP4MoAXG+r5OIBzmPlJIloYHnsHgP3M/EwiKgG4hYh+C+BYAOcDOJ2ZJ8QEY2FhURd2PFs0DCvEZy5+z8yjAEaJaD+AX4TH/wTgJCKaC+A5AP6HKNqdr5RS1i0ALiOiKwH8NDz2krAcYfJbgGBrxLMAfJuZJwCA1a0JLSwsmoMdzxZGWCE+c1GWPvvSdx/Be3cA7GPmU+oVxMzvClf2LwNwDxGdgmCDgfcx82/kc4noXEzT/rwWFjMYdjxbGGF94v2NUQDzmrmQmUcAbCKi1wEABTjZdC4RHcnMtzHzxwHsQrAf728A/A0RFcNzjiaiOQB+C+AviWg4PG7NbxYW+WDHs0XDsJp4H4OZdxPRLSEB5sEmingzgK8Q0T8BKAK4AsE2eDr+HxGtRrBavy485z4Ee+beRYH9bieA85n51+HKfi0RVQBcDeBjTbTNwmJWwY5ni2ZgtyK1SICILgPwy3aGpITlfgLAGDN/rp3lWlhYpMOO55kNa063MGE/gE9Rm5NDAHgLABtbamHRXdjxPINhNXELCwsLC4s+hdXELSwsLCws+hRWiFtYWFhYWPQprBC3sLCwsLDoU1ghbmFhYWFh0aewQtzCwsLCwqJP8f8D4fYC4KyfCm8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 576x216 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import sys\n",
    "sys.path.append('src/PGAM/')\n",
    "from GAM_library import *\n",
    "import gam_data_handlers as gdh\n",
    "import matplotlib.pylab as plt\n",
    "import pandas as pd\n",
    "from post_processing import postprocess_results\n",
    "import scipy.stats as sts\n",
    "\n",
    "\n",
    "# flag if you want to save the in-silico data\n",
    "save_output = True\n",
    "\n",
    "## inputs parameters\n",
    "num_events = 50\n",
    "time_points = int(0.5 * 10 ** 4)  # 30 sec at 0.006 second resolution\n",
    "rate = 150. * 0.006  # Hz rate of the final kernel\n",
    "variance = 5.  # variance of input (modeled to be Gaussian)\n",
    "int_knots_num = 20  # num of internal knots for the spline basis\n",
    "order = 4  # spline order\n",
    "\n",
    "## trial_ids: numpy.array of dim T, trial ids associated to each time point, assume 200 trials\n",
    "trial_ids = np.repeat(np.arange(2), time_points // 2)\n",
    "\n",
    "## create event (temporal) input\n",
    "idx = np.random.choice(np.arange(time_points), num_events, replace=False)\n",
    "events = np.zeros(time_points)\n",
    "events[idx] = 1\n",
    "\n",
    "## create continous (spatial) variables (spatial_var will drive the neuron, nuisance_var will not drive the neuron)\n",
    "cov = np.array([[1, 0.], [0., 1]])\n",
    "rv = sts.multivariate_normal(mean=[0, 0], cov=variance * cov)\n",
    "samp = rv.rvs(time_points)\n",
    "spatial_var_x = samp[:, 0]\n",
    "spatial_var_y = np.sqrt(5)*np.random.normal(size=time_points)\n",
    "nuisance_var = samp[:, 1]\n",
    "\n",
    "# truncate inputs to avoid jumps in the resp function\n",
    "sele_idx = (np.abs(spatial_var_x) < 5) & (np.abs(spatial_var_y) < 5)\n",
    "spatial_var_x = spatial_var_x[sele_idx]\n",
    "spatial_var_y = spatial_var_y[sele_idx]\n",
    "\n",
    "nuisance_var = nuisance_var[sele_idx]\n",
    "while (spatial_var_x.shape[0] < time_points) or \\\n",
    "        (spatial_var_y.shape[0] < time_points):\n",
    "    tmpX = rv.rvs(time_points)\n",
    "    tmpY = np.random.normal(size=time_points)\n",
    "    sele_idx = (np.abs(tmpX[:, 0]) < 5) & (np.abs(tmpY) < 5)\n",
    "    tmpX = tmpX[sele_idx, :]\n",
    "    tmpY = tmpY[sele_idx]\n",
    "    spatial_var_x = np.hstack((spatial_var_x, tmpX[:, 0]))\n",
    "    nuisance_var = np.hstack((nuisance_var, tmpX[:, 1]))\n",
    "    spatial_var_y = np.hstack((spatial_var_y, tmpY))\n",
    "\n",
    "    spatial_var_x = spatial_var_x[:time_points]\n",
    "    spatial_var_y = spatial_var_y[:time_points]\n",
    "    nuisance_var = nuisance_var[:time_points]\n",
    "\n",
    "# create a response function (for ground-truth)\n",
    "knots = np.hstack(([-5] * 3, np.linspace(-5, 5, 8), [5] * 3))\n",
    "resp_func = lambda X1,X2: sts.multivariate_normal(mean=[0,0], cov=cov).pdf(np.array([X1,X2]).T)*10.\n",
    "\n",
    "filter_used_conv = sts.gamma.pdf(np.linspace(0, 20, 100), a=2) - sts.gamma.pdf(np.linspace(0, 20, 100), a=5)\n",
    "filter_used_conv = np.hstack((np.zeros(101), filter_used_conv)) * 2\n",
    "\n",
    "# # mean of the spike counts depending on spatial_var and events\n",
    "log_mu0 = resp_func(spatial_var_x, spatial_var_y)\n",
    "for tr in np.unique(trial_ids):\n",
    "    log_mu0[trial_ids == tr] = log_mu0[trial_ids == tr] + np.convolve(events[trial_ids == tr],\n",
    "                                                                      filter_used_conv, mode='same')\n",
    "\n",
    "# adjust mean rate\n",
    "scale = np.log(np.mean(np.exp(log_mu0)) / rate)\n",
    "log_mu0 = log_mu0 - scale\n",
    "\n",
    "# generate spikes\n",
    "spk_counts = np.random.poisson(np.exp(log_mu0))\n",
    "\n",
    "# plot the firing rate and the spike counts generated\n",
    "plt.figure(figsize=(8, 3))\n",
    "plt.subplot(121)\n",
    "plt.plot(np.arange(200) * 0.006, np.exp(log_mu0)[:200] / 0.006)\n",
    "plt.title('firing rate [Hz]')\n",
    "plt.xlabel('time[sec]')\n",
    "\n",
    "plt.subplot(122)\n",
    "plt.plot(np.arange(200) * 0.006, spk_counts[:200])\n",
    "plt.title('6ms binned spike counts')\n",
    "plt.xlabel('time[sec]')\n",
    "\n",
    "### save the inputs in the standard formats\n",
    "\n",
    "## counts: numpy.array of dim (T x N), where T is the total number of time points, N is the number of neurons\n",
    "spk_counts = spk_counts.reshape(-1, 1)  # population of a single neuron\n",
    "\n",
    "## variables: numpy.array of dim (T x M), T as above, M the number of task variables\n",
    "variables = np.zeros((spk_counts.shape[0], 4))\n",
    "variables[:, 0] = spatial_var_x\n",
    "variables[:, 1] = spatial_var_y\n",
    "variables[:, 2] = nuisance_var\n",
    "variables[:, 3] = events\n",
    "\n",
    "## variable_names: numpy.array of strings of dimension M, name of each covariate in 'variables'\n",
    "variable_names = ['spatial_x','spatial_y', 'spatial_nuis', 'events']\n",
    "\n",
    "## neu_names: numpy.array of strings of dimension N, label uniquely identifying the neurons in 'counts'\n",
    "neu_names = ['neuron_A']\n",
    "\n",
    "## neu_info: dict, *optional neu_info[neu_names[k]]: dict, info about the k-th neuron, 𝑘=0,…,𝑁-1.\n",
    "## keys are the information label.\n",
    "neu_info = {'neuron_A': {'unit_type': 'SUA', 'area': 'dlPFC'}}\n",
    "if save_output:\n",
    "    np.savez('example_data_2D.npz', counts=spk_counts, variables=variables,\n",
    "             variable_names=variable_names,\n",
    "            neu_names=neu_names, neu_info=neu_info, trial_ids=trial_ids)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f19c6b24",
   "metadata": {},
   "source": [
    "## Create and save an example configuration file <a name=\"create-and-save-conf\"></a>\n",
    "\n",
    "\n",
    "Below the code for an example configuration file. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "448320fc",
   "metadata": {},
   "outputs": [],
   "source": [
    "import yaml\n",
    "import numpy as np\n",
    "\n",
    "# Create a YAML cofiguration file and save\n",
    "order = 4\n",
    "knots_1 = np.hstack(([-5] * (order - 1), np.linspace(-5, 5, 10), [5] * (order - 1)))\n",
    "knots_2 = np.hstack(([-5] * (order - 1), np.linspace(-5, 5, 8), [5] * (order - 1)))\n",
    "\n",
    "# convert to float (instead np.float64, produce a yaml human readable yaml list, optional)\n",
    "knots_1 = [float(k) for k in knots_1]\n",
    "knots_2 = [float(k) for k in knots_2]\n",
    "\n",
    "## change order of this example to match above\n",
    "cov_dict = {\n",
    "    'spatial_1': {\n",
    "        'lam': 10,\n",
    "        'penalty_type': 'der',\n",
    "        'der': 2,\n",
    "        'knots': [knots_1,knots_1],\n",
    "        'order': order,\n",
    "        'is_temporal_kernel': False,\n",
    "        'is_cyclic': [False,False],\n",
    "        'knots_num': np.nan,\n",
    "        'kernel_length': np.nan,\n",
    "        'kernel_direction': np.nan,\n",
    "        'samp_period': 0.006,\n",
    "        'join_var': ['spatial_x', 'spatial_y']  # for >1 dim response function\n",
    "    },\n",
    "    'spatial_2': {\n",
    "        'lam': 10,\n",
    "        'penalty_type': 'der',\n",
    "        'der': 2,\n",
    "        'knots': [knots_2],  # need to change this above, in the description of the variables\n",
    "        'order': order,\n",
    "        'is_temporal_kernel': False,\n",
    "        'is_cyclic': [False],\n",
    "        'knots_num': np.nan,\n",
    "        'kernel_length': np.nan,\n",
    "        'kernel_direction': np.nan,\n",
    "        'samp_period': 0.006,\n",
    "        'join_var': ['spatial_nuis']\n",
    "    },\n",
    "    'events':\n",
    "        {\n",
    "            'lam': 10,\n",
    "            'penalty_type': 'der',\n",
    "            'der': 2,\n",
    "            'knots': np.nan,\n",
    "            'order': order,\n",
    "            'is_temporal_kernel': True,\n",
    "            'is_cyclic': [False],\n",
    "            'knots_num': 10,\n",
    "            'kernel_length': 301,\n",
    "            'kernel_direction': 1,\n",
    "            'samp_period': 0.006,\n",
    "            'join_var': ['events']\n",
    "        },\n",
    "    'neuron_A':\n",
    "        {\n",
    "            'lam': 10,\n",
    "            'penalty_type': 'der',\n",
    "            'der': 2,\n",
    "            'knots': np.nan,\n",
    "            'order': order,\n",
    "            'is_temporal_kernel': True,\n",
    "            'is_cyclic': [False],\n",
    "            'knots_num': 8,\n",
    "            'kernel_length': 201,\n",
    "            'kernel_direction': 1,\n",
    "            'samp_period': 0.006,\n",
    "            'join_var': ['neuron_A']\n",
    "        }\n",
    "}\n",
    "\n",
    "# save the yaml config\n",
    "with open('config_example_data.yml', 'w') as outfile:\n",
    "    yaml.dump(cov_dict, outfile, default_flow_style=False)\n",
    "  \n",
    "\n",
    "# # load the yaml config\n",
    "# with open(\"demo/config_pgam.yml\", \"r\") as stream:\n",
    "#     try:\n",
    "#         cov_dict = yaml.safe_load(stream)\n",
    "#     except yaml.YAMLError as exc:\n",
    "#         print(exc)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "93b34659",
   "metadata": {},
   "source": [
    "## Fitting the model <a name=\"fitting-model\"></a>\n",
    "\n",
    "### List the fits by experiment, condition, neuron, and model configuration <a name=\"fit-list\"></a>\n",
    "\n",
    "Once the configuration file has been set-up, we will need a list of fits to be performed. An individual fit is identified by the experiment type, the session, the neuron ID, and the model configuation.\n",
    "\n",
    "When testing different models (including different subset of task variables, changing knots location and density, etc.) one needs to only modify the configuration files, and update the list of fits appropriately. \n",
    "\n",
    "We will again organize the list of fits as a YAML file with the following categories:\n",
    " \n",
    "- **experiment_ID**: string, experiment ID\n",
    "\n",
    "- **session_ID**: string, session ID\n",
    "\n",
    "- **neuron_num**: int, the neuron number (from 0 to N-1, whee N is the number of units in the session)\n",
    "\n",
    "- **path_to_input**: string, the path to the input data\n",
    "\n",
    "- **path_to_config**: string, the path to the configuration file\n",
    "\n",
    "- **path_to_output**: string, the path to the output folder where results will be saved\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "f09df6c9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# create and save the YAML with the fit list (1 fit only)\n",
    "fit_dict = {\n",
    "    'experiment_ID': ['exp_2D_spatial'],\n",
    "    'session_ID': ['session_A'],\n",
    "    'neuron_num': [0],\n",
    "    'path_to_input': ['example_data_2D.npz'],\n",
    "    'path_to_config': ['config_example_data.yml'],\n",
    "    'path_to_output': ['.']\n",
    "}\n",
    "# save the yaml fit list\n",
    "with open('fit_list_example_data.yml', 'w') as outfile:\n",
    "    yaml.dump(fit_dict, outfile, default_flow_style=False)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "082f96d7",
   "metadata": {},
   "source": [
    "### Load inputs, fit, and postprocess <a name=\"load-fit-save\"></a>\n",
    "\n",
    " \n",
    "The code below loads the input data, the configurations, fits the model, post-process, and save the results. The following lines of code are also saved as separate script in 'PGAM/utils/fit_from_config.py'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "0d0d53c3",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "FIT INFO:\n",
      "EXP ID: exp_2D_spatial\n",
      "SESSION ID: session_A\n",
      "NEURON NUM: 1\n",
      "INPUT DATA PATH: example_data_2D.npz\n",
      "CONFIG PATH: config_example_data.yml\n",
      "\n",
      "\n",
      "adding events...\n",
      "adding neuron_A...\n",
      "adding spatial_1...\n",
      "adding spatial_2...\n",
      "\n",
      "fitting neuron neuron_A...\n",
      "\n",
      "tot fit time: 23.41 sec\n",
      "post-process fit results...\n",
      "hstack: 0.016701799002476037 sec\n",
      "hstack: 0.008491785993101075 sec\n"
     ]
    }
   ],
   "source": [
    "# import libs\n",
    "import numpy as np\n",
    "import sys, os\n",
    "\n",
    "sys.path.append('src/PGAM/')\n",
    "\n",
    "import GAM_library as gl\n",
    "import gam_data_handlers as gdh\n",
    "from post_processing import postprocess_results\n",
    "import yaml\n",
    "import statsmodels.api as sm\n",
    "from scipy.io import savemat\n",
    "from time import perf_counter\n",
    "\n",
    "np.random.seed(4)\n",
    "\n",
    "#################################################\n",
    "# User defined input\n",
    "#################################################\n",
    "\n",
    "# frac of the trials used for fit eval\n",
    "frac_eval = 0.2\n",
    "\n",
    "# PATH to fit list\n",
    "path_fit_list = 'fit_list_example_data.yml'\n",
    "\n",
    "# save as mat\n",
    "save_as_mat = True\n",
    "#################################################\n",
    "\n",
    "\n",
    "# load the job id (either as an input form the command line, or as a default value if not passed)\n",
    "argv = sys.argv\n",
    "if len(argv) == 2:  # assumes the script is run the command \"python fit_from_config.py fit_num\"\n",
    "    fit_num = int(sys.argv[1]) - 1  # HPC job-array indices starts from 1.\n",
    "else:\n",
    "    fit_num = 0  # set a default value\n",
    "\n",
    "# load fit info\n",
    "with open(path_fit_list, 'r') as stream:\n",
    "    fit_dict = yaml.safe_load(stream)\n",
    "\n",
    "# unpack the info and load the data\n",
    "experiment_ID = fit_dict['experiment_ID'][fit_num]\n",
    "session_ID = fit_dict['session_ID'][fit_num]\n",
    "neuron_num = fit_dict['neuron_num'][fit_num]\n",
    "path_to_input = fit_dict['path_to_input'][fit_num]\n",
    "path_to_config = fit_dict['path_to_config'][fit_num]\n",
    "path_out = fit_dict['path_to_output'][fit_num]\n",
    "\n",
    "print('FIT INFO:\\nEXP ID: %s\\nSESSION ID: %s\\nNEURON NUM: %d\\nINPUT DATA PATH: %s\\nCONFIG PATH: %s\\n\\n' % (\n",
    "    experiment_ID, session_ID, neuron_num + 1, path_to_input, path_to_config))\n",
    "\n",
    "# load & unpack data and config\n",
    "data = np.load(path_to_input, allow_pickle=True)\n",
    "counts = data['counts']\n",
    "variables = data['variables']\n",
    "variable_names = data['variable_names']\n",
    "neu_names = data['neu_names']\n",
    "trial_ids = data['trial_ids']\n",
    "if 'neu_info' in data.keys():\n",
    "    neu_info = data['neu_info'].all()\n",
    "else:\n",
    "    neu_info = {}\n",
    "\n",
    "with open(path_to_config, 'r') as stream:\n",
    "    config_dict = yaml.safe_load(stream)\n",
    "\n",
    "# create a train and eval set (approximately with the right frac of trials)\n",
    "train_trials = trial_ids % (np.round(1 / frac_eval)) != 0\n",
    "eval_trials = ~train_trials\n",
    "\n",
    "# create and populate the smooth handler object\n",
    "sm_handler = gdh.smooths_handler()\n",
    "for var in config_dict.keys():\n",
    "    # if var == 'spatial_1':\n",
    "    #     continue\n",
    "    print('adding %s...' % var)\n",
    "    # check if var is a neuron or a variable\n",
    "    var_list = []\n",
    "    for join_var in config_dict[var]['join_var']:\n",
    "        if join_var in variable_names:\n",
    "            x_var = np.squeeze(variables[:, np.array(variable_names) == join_var])\n",
    "        elif join_var in neu_names:\n",
    "            x_var = np.squeeze(counts[:, np.array(neu_names) == join_var])\n",
    "        else:\n",
    "            raise ValueError('Variable \"%s\" not found in the input data!' % var)\n",
    "        var_list.append(x_var)\n",
    "\n",
    "    knots = config_dict[var]['knots']\n",
    "\n",
    "    if np.isscalar(knots):\n",
    "        knots = None\n",
    "    else:\n",
    "        knots = np.array(knots)\n",
    "\n",
    "    lam = config_dict[var]['lam']\n",
    "    penalty_type = config_dict[var]['penalty_type']\n",
    "    der = config_dict[var]['der']\n",
    "    order = config_dict[var]['order']\n",
    "    is_temporal_kernel = config_dict[var]['is_temporal_kernel']\n",
    "    is_cyclic = config_dict[var]['is_cyclic']\n",
    "    knots_num = config_dict[var]['knots_num']\n",
    "    kernel_length = config_dict[var]['kernel_length']\n",
    "    kernel_direction = config_dict[var]['kernel_direction']\n",
    "    samp_period = config_dict[var]['samp_period']\n",
    "\n",
    "    # rename the variable as spike hist if the input is the spike counts of the neuron we are fitting\n",
    "    if var == neu_names[neuron_num]:\n",
    "        label = 'spike_hist'\n",
    "    else:\n",
    "        label = var\n",
    "\n",
    "    sm_handler.add_smooth(label, var_list, knots=knots, ord=order, is_temporal_kernel=is_temporal_kernel,\n",
    "                          trial_idx=trial_ids, is_cyclic=is_cyclic, penalty_type=penalty_type, der=der, lam=lam,\n",
    "                          knots_num=knots_num, kernel_length=kernel_length, kernel_direction=kernel_direction,\n",
    "                          time_bin=samp_period)\n",
    "\n",
    "link = sm.genmod.families.links.log()\n",
    "poissFam = sm.genmod.families.family.Poisson(link=link)\n",
    "\n",
    "spk_counts = np.squeeze(counts[:, neuron_num])\n",
    "\n",
    "# create the pgam model\n",
    "pgam = gl.general_additive_model(sm_handler,\n",
    "                                 sm_handler.smooths_var,  # list of coovarate we want to include in the model\n",
    "                                 spk_counts,  # vector of spike counts\n",
    "                                 poissFam  # poisson family with exponential link from statsmodels.api\n",
    "                                 )\n",
    "\n",
    "print('\\nfitting neuron %s...\\n' % neu_names[neuron_num])\n",
    "t0 = perf_counter()\n",
    "full, reduced = pgam.fit_full_and_reduced(sm_handler.smooths_var,\n",
    "                                          th_pval=0.001,  # pval for significance of covariate icluseioon\n",
    "                                          max_iter=10 ** 2,  # max number of iteration\n",
    "                                          use_dgcv=True,  # learn the smoothing penalties by dgcv\n",
    "                                          trial_num_vec=trial_ids,\n",
    "                                          filter_trials=train_trials,\n",
    "                                          compute_MI=False)\n",
    "print('tot fit time: %.2f sec'%(perf_counter()-t0))  \n",
    "\n",
    "print('post-process fit results...')\n",
    "res = postprocess_results(neu_names[neuron_num], spk_counts, full, reduced, train_trials,\n",
    "                        sm_handler, poissFam, trial_ids,\n",
    "                          var_zscore_par=None, info_save=neu_info, bins=15)\n",
    "\n",
    "\n",
    "# saving the file: save_name will be expID_sessionID_neuID_configName\n",
    "\n",
    "config_basename = os.path.basename(path_to_config).split('.')[0]\n",
    "save_name = '%s_%s_%s_%s'%(experiment_ID, session_ID, neu_names[neuron_num], config_basename)\n",
    "\n",
    "if save_as_mat:\n",
    "    savemat(os.path.join(path_out, save_name+'.mat'), mdict={'results':res})\n",
    "else:\n",
    "    np.savez(os.path.join(path_out, save_name+'.npz'), results=res)\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6737c04e",
   "metadata": {},
   "source": [
    "# Plot the 2D tuning function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "a2c15d88",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<mpl_toolkits.mplot3d.art3d.Poly3DCollection at 0x7fcfafbbaf98>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPgAAADuCAYAAAAUXsqNAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAACfUElEQVR4nOz9d3hcWXrei/7W2rsyCjkQJAEiMYABTCDZDM3u6QmakeZey5bla10H2bJly762ZR/bj3Uk2zoOkuU0PrLkJB1rLNuSLdmSJXlS93SabrKbzUwwI5EIBIgcKlftvdb9Y2EXQRChAALdnB68z8MHJFG1a1fVfveX309ordnEJjbx6YT8pE9gE5vYxMZhk+Cb2MSnGJsE38QmPsXYJPgmNvEpxibBN7GJTzE2Cb6JTXyKYa/w+80a2iY2sfEQG3XgTQu+iU18irFJ8E1s4lOMTYJvYhOfYmwSfBOb+BRjk+Cb2MSnGJsE38QmPsXYJPgmNvEpxibBN7GJTzE2Cb6JTXyKsUnwTWziU4xNgm9iE59ibBJ8E5v4FGOT4JvYxKcYmwTfxCY+xdgk+CY28SnGSvPgm9gAKKVIpVJYloVlWdi2jRAbNhK8ie9hiBV00TcFH9YRWmscx8FxHLLZLPM/e8uy8Pl82LaNZVmbhP/ewoZ92ZsE/5igtSabzaKUQghBLpd76nfeHw+bhP+ewibBv5vhOE6e0B5R5xN8ITYJ/z2HTYJ/N2K+Sy6EyJPSs+aFktQju2f9YZPwnzJsEvy7DUopcrlcnpTzCbhagi+ER/ienh4qKyspLi7Gtu38n03Cf9dhw76szSz6OkNrjeu6eRdcymcrkVprJiYmiEQihEKhVb+Gd8NQSiGlREqZf02P2LZt4/P58pn6TcJ/b2KT4OsIrTW5XA7XdZ+x2h4ymQwdHR34fD5yuRyZTIZoNEpZWRllZWUEAoGCX08IgdY6/1rezcS7yTiOk3+sR3jbtpFSbhL+ewSbBF8nKKUYGBggEAhQWlq6KIHGx8e5f/8+u3fvpri4OO9qx2IxpqamuHPnDo7jUFxcTFlZGaWlpfj9/lWfy2IhwXzCCyGecuk3Cf/pxSbBnxPzE2nxeBzgGbIopeju7mZmZob29nYCgQDZbBYwLnxJSQklJSU0NDSglGJmZobp6WkGBwdxXZeSkpI84X0+36rPcTHCL8zsbxL+04lNgj8HFta2pZQsTFqmUik6OjqorKykvb19ReJIKfPuemNjI67rMjMzw9TUFP39/WitKS0tpaysDKXUM69XCBYjfC6XIx6PMzw8TGNj4ybhPyXYJPga4SW15sfAXkzsYWRkhO7ubvbu3UtZWdmaXseyLMrLyykvLwdMTX16eprJyUnGx8eZmZmhqqqKsrIyiouLsSxr1a8hhMg/L5VK5Rtx5lv4+SW5TcJ/92CT4KvEwtr2/Cy5R3ClFPfu3SOdTnPs2LE1xdFLwbZtKisrqaysBKC0tBSAsbExuru7sW07b+GLi4sXzeKvhPmEhyeeSiaTyd/IPMJ7ffSbhH8xsUnwVWC52jYY9zqdTvPRRx9RW1tLa2vrhl/4tm1TVlZGVVUVANlslqmpKR4/fkxnZyd+vz/v8kej0WXPx/NGFmI5woN53z6fL1+W2yT8i4NNgheAhbXtpS7gmZkZJiYmOHz4MCUlJR/3aQLg9/upqamhpqYGgHQ6zdTUFIODg8TjcYLBYJ7wkUhkTUScT3gvJMlms08lDufX4dfiRWxifbBJ8BVQSG3bcRzu3r1LPB6nvr7+YyX3Skm2YDBIbW0ttbW1aK3zhO/r6yMejxMOh/OEn98KWyjmt87OPx+P8MlkEtd1qaio2CT8J4BNgi8DpVR+rHMpcsdiMW7evEl9fT1lZWVPNZdsNNZCxlAoRCgUYuvWrWitSSaTTE1N0dvbSzweR2vN0NAQZWVlBIPB5yZ8IpEgnU4TjUY3LfwngE2CLwIvkfbgwQPq6+uXbDcdGBjg0aNHtLW1UVRUxNDQ0JrKVp8UhBBEIhEikQjbt28nFovR29uL67p0dXWRTqcpKirKW/hgMLim15FSPmXhvRh+PuEXZuk3sT7YJPgCeBef67oMDg6yY8eOZx6Ty+W4ffs2Pp+P48eP5y/exergS2FhSe1Fgd/vp66ujrq6OpRSxONxpqamuHfvHtlsNt9lV1ZWVlB1YGHibrEavNaaTCaTT9p5k3KbajfPj02Cz8N8l3ypWu/09DS3b9+mubmZLVu2PPU7bwDk0wIpJcXFxRQXF7Njxw6UUszOzjI1NcWjR48K6rJbKjPvYTHCK6VIp9P5/9scjV07NgnO8rXt+Y95+PAhIyMjHD58mHA4/MxjVmOVM5lMftLrec99vbASGaWUlJaWUlpaumSX3XzC27a94jEXohDCe99RKBTaJPwK+J4n+Eq1bTAZ4Zs3bxKJRDh+/PiSMWIhBNda09/fz8DAAAA+n4/y8vKC6tSLvd4nicW67DzCP3z4MF9Oi0QiuK675i67hYSPx+P09fXR2tqaP49NC784vmcJXmhte2Jignv37rFr1658M8lSWIngjuNw69YtfD5fvi89m80yOTnJ4OAgsVjsqbJVOBz+WC/W1VrbhbBtm4qKCioqKgCTq+jq6iIej3P16lUsy8q/t+fpsgPy3o9n4b0WW+93m4Q3+J4k+FJSSgsfk06n6enp4ejRowVlkJeLwWdnZ7l16xaNjY3U1tbmvYZAIPBUnXp+2SqZTOZnxcvLy1c1K/4iwOfzEYlEqKiooKamhmw2y/T0NCMjI3R2duLz+Z7qsiuU8J7QBTy5Mc+fhd8k/BN8zxG8EJc8nU7T0dGBEIJjx44VfEEslkXXWjM4OMjg4GC+nLYUFpatFs6K53I5SkpKKC8vz/egv+iY7xX4/X6qq6uprq4GTB5iamqKoaEhYrEYgUAgT/iioqIlP/f5BF+IQgj/vSRv9T1D8Pku+VKJNIDR0VG6urpobW3l3r17q46J5xPccRxu376NZVlPldNWc7yFWeyZmRkmJyfp7+8nlUqRTCYRQlBSUrIuCbv1vtiXO2YgEGDLli35akQqlcon7BZ22c0PV1bTcfe9rnbzPUFwLzGz3OSTUorOzk4SicSaJ8DmE9zrcNuxYwfbtm177vcAT8+KA3R3dyOlZHx8nJ6envzgyWpd3o3Eam4ai3XZTU9P8+DBAxKJBJFIJP/e10rC5dRuhoeH897Dp4Xwn3qCe7XtS5cucerUqUW/rGQySUdHBzU1NezevfuZC2A11sJzyfv7+1d0yZ8XXp3aGx31Enaey7vawZKP24Ivh/nhyrZt29Bak0gk8pNyXmPM83bZzSf8zMwMpaWlnyp5q08twRdLpC32xQwPD9Pb28v+/fufGRKRUqKUKtj1VUoxOTmJ1prjx49j2x/vx+v3+/Mur9Y67/I+fPiQRCLxVNvpWtRc14L1umkIISgqKspbV4/cU1NT3L9/n0wmQ3FxcX4Wfi0JSe+7nu/5fLfLW30qCb5QSmmxL8B1Xe7evYvjOBw/fnzRLqzVtJ7G4/F8CezAgQMrfunrdVEsdX5CCMLhMOFwOG8BvbbTzs7OPCG8hJ3f79+Q1tmN8Aq8JNvC/MRS4pVlZWUFadktlrxbzKVfqHazcHDmRSL8p47gC++2i33Y8Xicmzdvsn37drZv377kF1Jo6+nQ0BAPHz5k165dPHr0qKAveL2s2moeG41GiUaj1NfXP9V2Ojg4mC/ZgfkM18v72Ci3fyERF4pXuq771PtzXTdv3b0uu4UopBln/iy8dy6LyVv19vbS2NhIJBJZh3e8dnxqCF5IbVspxdDQEP39/Rw4cIBoNLrsMT0XfSks9AI8jbZCz/eTxMK2U8dxGBwcZGxsjOvXrz+V0FtrUwpsnAVf6QY0v6kGzHc1PT2dD1mAPNlLS0uxLGvZ8ttSWIzw2WyWv//3/z7/5J/8E3bv3r26N7fO+FQQvJDathCCjo4ObNsuOD5ezkVPJBJ0dHSwbds26urqEELgOM4nTty1wrZtotEoruvS3Nycb0rxpJ8KrVEvxEYRfLXHtCzrqS67+eKVvb29WJZFOp1mZmZmzeKV8ITwXtb/k8Z3NcELWRMEJjsaj8fZvn07dXV1BR9/KRf98ePH9Pb2sm/fvqcSc6uJ2V90LGxKWVij9kpWXsJuKcJtZAz+PFgoXpnNZrl8+TKjo6N58crnKTkmEokVPcSPA9+1BC9ESklrTV9fH48fP6akpGTFXvKFWOiie2qpmUyGY8eOPZO4Wc00mZcEfF7F1Y9rmmxhjdorWXV3d+cVW7yW2vnv6eOKwZ8Xfr8f27bzLvViXXZeDF/IUFAymVx04vDjxnclwQuRUspms9y6dYtQKMTx48e5cePGqme15xPWq5Uvp5ZaaFJudHSUzs7OfMbVI8ZqY91PKls7v2TlCUN4Gezbt2/jOE6+pdbLiawn1sOCr4TFuuy8bTOxWIxQKLRsj4HW+mMvky6GT/4MVgHPcsRiMcrKypb8kicnJ7l79y47d+7Mu5grJcwWg/ccb4HBvn37lu0BX8mCa63p6upidnaWI0eO5LeCzpc5DgQC+RHMj3uabK1YLIPttdROT09z586d/HsqKSl5bnKuJQYvBMsd0/NgvKGghT0G80OWYDC4as9KCBEE3gMCGF7+T631zy54zKvA7wMP5v7rd7XW/3C5437XENzLTsbjcYaGhvIzyAsf09vby8TEBEeOHHmqmWMtBBdC8PDhQ5RSBbWvLkdwb6toWVkZR48ezXdL+Xy+Z2JdL/HjTZN58+IbPU22Xu70/DnxZDJJU1MTqVTqqfh2rTPw3nmutwVfDSEX6zGYH7L8xb/4F0kkEnz1q1/ltddeo6GhoZDDZoDXtNZxIYQPOCeE+KbW+sKCx72vtf5yoef6XUHw+WuCvHLGQqTTaW7evElpaSnt7e2L1klXQ/BUKsXIyAjV1dXs3bu34Nr2YheK13xRyEx5KBRi27Zt+QsnFosxOTmZb94oLS39rpsm8/l8FBUV5d97JpN5rhn4jXDRn+eYC0OW999/nzNnzjA1NcUv/MIv8O///b9f8RjaXDjxuX/65v48d4LlhSb4YlJKtm0/Q9SxsTE6OzvZs2dPvgyyEKshuDdRVlFRQXV19ap60Reev5fkW+hRFHo8r1vLc33nl3ZyuRzFxcUEg8E1WcKF+Lh60ZebgU+lUvmW2qVm4DeC4K7rrtsxXdclEonwN//m31zV84QQFnAFaAH+jdb6o0UedlIIcQMYAv6W1vr2csd8YQm+VG3bi1u9x3R1dRGLxfJreZdCIQT3jhePxzl27Bh9fX1rzlLPV29ZTOZpLURaWMvt6enJN6h4ltBzfZcrXX2cKER0cbkZeK/l1PNafD7fhsTgq5k5WAleGXG10Fq7wCEhRCnwv4QQ+7XWt+Y95CqwY86N/37g94Cdyx3zhSP4SlJKnoueTCa5efMmVVVVHD16dMUvfCWCeyIPFRUVHDlyJO8xrEUlNR6P09HRsa6joovBtm2KioqoqanJW8LJycl86cojRqG92BuBtYguLjYDP1/YMZfLMTs7SyAQWDdSrqcFj8fjzzVFqLWeFkK8C3wRuDXv/2fn/f0bQoh/K4So1FqPL3WsF4rghdS2pZSkUimuXbu2YlZ74fOWIqvn4re2tj6VvFsLwXO5HB0dHQW1wq4HPA9jviX0SldeL/bAwABa63ycu5Q4xIs0Luph4Qy84zhcvXqV6elpBgYG1m0Gfj0teDKZXDXBhRBVQG6O3CHgc8A/XfCYLcCI1loLIY4DEphY7rgvDMELqW27rpsX4H/llVdWZZUWI6vWmu7ubqanpxd18VfTuKKU4v79++RyOU6fPv2x1ECXI85ivebT09N5cQhPD628vDzfevrdME3mjWm2tLRg2/Zzz8B7WM+4PpFIrKXJpRb49bk4XAK/rbX+mhDiJwC01v8e+KPAXxJCOEAK+ON6hS/tEyd4oVJK3gTYtm3bmJ6eXrXLuZDg88tWnsLpSs9ZCul0mhs3blBdXU0oFHohGhwWYmFrpreEcH7rqWVZ67rLHDa+F329ZuDXKuu8GNbiomutO4DDi/z/v5/3918Gfnk1x/1Er8RC5rYBHj16RF9fH/v376e4uJjBwcFVv9Z8snpSyMtl3Rc+Zyl4x/Lc++Hh4VWf28ZAEwr8v5D2fdLpf4Hr/qGnfrtw62gikaCvr4+xsTEmJibynWhlZWXPfcP6uFpV1zID72E9LfjzxuDriU+M4EopHj16RGVl5ZLkdhyHO3fuADy3QoplWWQyGbq7u5mcnCxICnm51tP5TTWFyip/nAj4/gS27z0QEPT/NXLZ18m6/3bRx3p1XK9ddtu2bflOtL6+vudqp90oFNqXsNIMvHcjy2az6/a+ksnkCzFJBp8AwefXtru7u/MdXAvh6YivVybadV0GBgbYsmXLoo0wi2EpC57L5fKbTgo91kZh8RCsB5/8JtoNgJVBOHH84r/h6B9Cic+ueMyFia1cLvdd304Li+clvBvZ2NhYPg/xvDPwXmjwIuBjJfjC2vZi0HOrfYaGhtZNtNBrDCktLWXXrl0FP2+xxJN342lqanpm+eBqkEgkiMfjlJWVPdfs8WII8pfRrh/hphBZwDIlx6D1MyTV0gRfKl/zIrXTrifmb2IJh8M4jkMoFHruGfh4PE5tbe0Gn31h+FgIvlJt20Mul+PWrVsEAoFldcQ913mlO6zWmgcPHjA2NsbOnTuJxWKrOu+FFnxwcJCBgQEOHjz4XC6YJ/FUXFzMgwcP8r3Z8zPaa4Wl/ge2uAAIhAAc0EIgpEbIB0jnbZR8bcnnF/Lan7Z2WjAeXiAQWJcZ+LWUyTYKG07wlaSUvCyr17nU0tJCTU3Nssf0iLccwb2FgUVFRRw7dozp6WlmZmZWde7e63jSTN7QyUq5gKUyx14pLZ1O097enn8PXm+2dyF5FnHhbHUh8PMrCDSIOWusQGQ1BEGQJiD/GSmWJvhqsVI7bTKZ5OHDh3m390V15xe7ntY6A79aNZcCJ8kE8IvA9wNJ4M9ora+udOwNJfhKUkqeJfayt4X2a1uWheu6SxLNu1k877ioEGY54MWLF5+SZlrpOYsRPJPJcOPGDSorK9mzZ0++ggDP9mbH43EmJibys9XeRbTi9hIVwxLX0cKPIIvWFmgXkKA0CInkLmgNSwhkPC8BF7bTfvTRRwSDQYaGhrh3715+sKS8vPyFaaeFlctkhc7Ax+NxZmdnV2vBC5kk+xKmLXUncAL4d3M/l8WGELzQ2raUkitXrlBSUsKxY8cKTmp4BF/sdR8+fMjo6Oi6jIvOzMwwOjpKe3v7M5rpS2GxuH16eprbt2+ze/fufB16ObljL/Pb0NCwaIOKZ9211k8dx+/8c4Q/g6IYQRZciRAuoFDZKDIXQ2ccrPD/xLV/eFWfxVohpXyqTu0Nlsy3gt77+aTaaWH1ZbKlZuC//vWv89Zbb3HlyhW+/OUv89f+2l9bMpHsocBJsj8E/Oe5x14QQpQKIWq11svWZded4IXWtsfHx4nFYrS2tq46S74YWb3MdigUWvRmsRqCe8IMk5OTVFVVFUxueJbgAwMDPHr0iMOHD69JwmexBpXJyUkePnzI9PQ0oVAoryAa4pvmsnAc83P+TcRRxk5o8Gd+jdQiBN+IppT5WDhY4llBb3RUa52P35drp90IPG+rqjcD/zf+xt/g2rVr/MzP/AyDg4MFh1gFTJJtAwbm/Xtw7v8+PoIX0m6qlKK7u5uZmZn8F7laLLTgMzMz3Lp1a9n4vVCCZ7NZbty4QVlZGa2trfT396/q3DyCK6W4c+dOPm5fry6pYDDI1q1b2bp1KwMDA2SzWZLJJGOPr3N01xBu1sZSSQC0P4Lwz+Ud3BxoECKDdG4v6aZ/nC7zfCu4WDvtYsnHjZJrWu9hk61bt9Le3l7wcwqYJFvsi1nxbrduBJ8fUy71QaVSKTo6OqisrKS9vZ3bt28v6mqvBI+sXklteHh4RQtZCMEXCjPE4/E1yTylUinu3LnD1q1bC4rb1wohRD4J5E99FTunIAd4RiM2iy6TCKkQuSxaAVYWoXNI5yLKt2II97GikHbakpKSDbHi692LvtZBo6UmyTAWe74k8HbMTPiyWDeCe7H2Uh++p2u2d+/efAPFUrH0SvC60q5fv04gECjIQi5H8Pk3ivmx+1oGMHK5HDdu3GD//v3597mR8M7PTn0LoRJoEQTM4jwpNG4yjFWUQORAaR/SzqG1xM7+HtkFBN9oF321WKyddnR0lFQqxcWLF9e1nXY9e9HT6fSqxD0KmSQD/gD4K0KI/45Jrs2sFH/DOrvoixHCkxpOp9PP6JqtleC5XI779++za9eughsKliL4csIMq43bHz58SDKZ5OjRox9LDThPRtdBuCMIC3Czc+djI4SDzKVAMefMmffmah/ZmW/QO/LjG9p+up6W1sti27ZNPB5n//7969pOu54WfA3xfCGTZN/AlMi6MWWyP1vIgTe0TOZt/1hKani1BNdaMzAwwPj4OA0NDavqFlqMrCsJMxS6yMC7Sfj9fsrKytZ9Imsl2JnfBleBBQLlnRRg/q2TJoATGDfdsiwiwRGKrEi+aysYDOZr2uuFjZokk1KuezvtesXga7mpFThJpoH/32qPvWEE97q1Fm7/mI/VENxxHG7fvo1lWdTX169pXHT+hz88PMyDBw+WFWZYbtjEg3cTq6+vZ9u2bVy/fn3DMr1LwU7/PmAsdz4Voy3AfLY6U4wQswg0OisgkEFIh9rwFaprvj8/Ztnd3c3IyAgjIyPr4v5u9KjofDxvO+16J+9elFBn3V10x3HyXV8rTYAVSvBYLMbNmzdpaGjIZ4/X4trD091ki20nmY+VXHRPnPHAgQN5y7ea9UXrdRHI9E2ElZ77h/mhdQgh5kqrSQ1hQIDIaQg6aA129ndw+f78mKU3QlldXf2M++tZw9V0o22UQkwhRCykndbzACzLWraku9rze1HIDetM8NnZWW7cuJG3Ziu9US9ZthS01jx69IiBgYGnBk+klPm+9tVAKcWlS5eorq5mz549Bem4LUZWrTU9PT1MT08/k1fYKGWUpaBVDqnHED7zmtqaM+KavDUXbu6p7LpOBdCTPmTJ04Kc3sW50P3NZrP5NT5eN5pH+OWSSRvpoq8GS7XTesIQXruwt3jwec45k8m8UAM360rwZDK5qgmw5Sy4NwsuhHim/9vbBLkaTExMkEwmaW9vX3RpwmJYzEX3NNei0eiiYo+FuPVgRCxGRkaeihNXCyEEYfUOggxaSkRCgRMAXwbhNzdArQE3bTx4P2AB01lEzkEkH4PKgVw+3PH7/dTU1Dwj7uiJKMwfLpn/PRVqbVeD9XClF7bTeosH16Oddq2KqhuFdSV4bW0tzlxypxAsRXAv+VVfX8/27duf+f1qs9sPHjxgfHw8b3kKxcIv1gsVlhsVXcmCe1WFbDZLU1MT09PT+bbNkpISKioqllxQvxhK9NfRwkJnihDxGYSVMdbazUERiBzGks99LVpZiKyLFhqZSWPFvo5b8oPmdwVY3PndaHV1dU+tKfKs4fxBmY8rBn8e+P1+fD4fra2ty7bTFpJAfZHUXOATlmxajOCePNNyya9CY3evfTUcDtPe3s6FCwu3wBQOLym3koeyHMG9gZOqqip2795NLpejqKgo37Y5MzPDxMTEUyOkFRUVy4oHBtxOcEEkZtByzitXIOIKHbIQztznpOf+pP2gUmAJ0Cnsqd/PE3wtmL+mCMiLIA4ODjIzM4NSiqGhISoqKtbFdd0Ir2D+97VcO+2jR49QSi3bTvsiqbnAC0RwbyTTdd0Vk3OFWPD1EmbQWnPv3j1SqVRBslFLEdxrp/W65BZaooVxrzdC6okHRqNRKioqnhrKkGoKnxpDJF0IgJJRBDG8SplOBBEk5k4MY9nn+tGFcNEKZOpJHL4eMfN8EcREIkFXV1c+3Fo4K76WxpKNWlu01PtebTvtakdFAYQQdcB/BrZgvr1f0Vr/4oLHvMoqFw/CBmTRVwOP4F6pqdCRzJUs+FLCDKvtZfb6vLds2cLu3bsLen+LJeaGhobo6+vj0KFD+fNZ6VgLR0hnZ2fzlhGgrKyMyuRvgs9Fa2nq34423+jcy4tU5knbKpg4PJ3K/1MokHpqxff0PPD5fNTX11NfX//MrPhahC42iuCF3myWa6e9dOkSv/3bv00kEqG/v5/6+vpCT8EB/qbW+qoQIgpcEUJ8W2t9Z8HjVrV4EF4AC55IJPKtnYU2WSxF8JWEGQoRivDgWdxAIEBzc3Nhb4ink2xKKTo7O0mlUgUJRSx3zPlWxGvqCE6cR7kSLfzYpJ9MjylAgtAOpAFPDzK92MRCGpG6gw7tXdO5LYeFHsHC5NZCoQtP+HGpnWSwMTH489w05rfTtra2kkwmeeutt/iJn/gJvvzlL/OX//JfXvEYcy2nw3N/jwkh7mImxRYSfNX4xAjuui5dXV2k02nOnj373EsMkskkHR0dyw54FJqc8zyAw4cPc/369YLPC5646Nlslo6ODkpLSzl06NC6XpReU0dgcBg7m0J7pSptml1cgljMVRmy5AmuMiEkxoJrbSGEiyCHPfV75EJ7172stdLxlhK68Nz5xYQuNiIGX68+dC/B+Nprr/FTP/VTazqGEKIB09X23IsH4RNy0T0ybtmyhVgstuqutIUW3Gs4WWmV0UoEV0o9lQdYy5cuhCCZTNLV1fWUosy6IzOClZtCSD2XKgdheb0B5mvVGkQadNRMhspUCu2zEJYLORttKyCLFbvE6rsKVsZqbhiFCl1ks9k1lRSXw4uiiS6EKAJ+B/jret4esjmsevEgfAIW/PHjx/T09LB//35KSkoYGlpx4u0ZzB8X7erqYnZ29pmGk+Wetxi87SRbtmyhvr5+zZYskUgwPT3N0aNHN7Rc4hv9r08y4yKH1gJhGxddqoxx0V1AgEqHsEKpuc5V07NONjenG+Igk13A+jemPM/xlhK6GB4eZmRkhOnp6TXr1i3Eek6SJRIJ6urqVn7gAsxJNf0O8Bta699d+Hu9hsWD8DES3GsR9bLRzyPPY1kWuVyOy5cvU1paWtB2UVia4JOTk9y9e/eZ5YMeCrlQ9dyes9nZWZqbmze8FmrFzoOQRjFVaLQKgOeWi7mi95yTY8VT6KBJqglvWaFWIMBVAYQzil7DFtWVsJ43DE/oIpPJEA6HCYVCTExMcOvWLZRST7nzn+QkWSqVWksWXQD/Ebirtf7KEo9Z9eJB+JgIPt8lL6RFdCXMzMyQTCbZvXs3VVVVBT9vIcG11vT19TEyMrLkdpKlRBTnw3EcOjo68jXtj2MRgkz2IYSLkn4sMuicD+Gk0QFAziXbvChGY8pjgnwJzcu2qVQAm1nufPjrxG0ThxcXF6+LPtpG9aJblpVvPfVKV1NTU/lQbbWTZC+Ai34a+FPATSHE9bn/+2mgHta+eBA+hhi80Pi4EMwXZgiHw6siNzzZLQ7GLbt16xa2bS8r+LhS5t3rumtsbKS2tpa+vr417RRfFdw0Ij2C8GdRhGE8g0inQIL2FyHKzJCJFhFTBxcgEguOYYOOh7EfxyAkOdzcz7XU8XxyEEwprqKiYs392R9XL7pt21RVVeWvh/mTZKlUKt9DsNSe9PV20dewePAci0syzX/MqhcPwgZY8PmaZJ2dnSQSiRXj40IuBG9c1LZtjh8/vqauNK+ElUwmuXHjBnV1dYu2ws7HctNh3l7x+dNkH8ewiTX2vxDZLCKgYcYynWzKMXF3PIG2ihAlcbQz76pJYmJuCVr5ESILk8qMkOYCyMmrBEr/RN7yeaU4rz87EonkO+sK7Uj7pIZN5k+SeZ1oExMTDAwYzcKFwhDracFfpKUHsEEuuqe95rVkLvclex/wcnfQlYQZCoWUksnJSUZGRvJJvpWw2PCInls8ODk5ueZpMtd1mZqaWlnrfBHYY99GIxBJkIn4kzo3IJSG2Tg6bCGyCUNqMF1s3t9dPzhZtGsb3UWVwZq4iy55Qsj589WeXNL8cUuPJKWlpUuSY6Nc9NUcc34nGjwrDBEMBpFSrlt76Yu0lww2gOBjY2Pcv39/yYTVQnglr6Uu8kKEGQqB1pqZmRlyuVxBGXcPC+N2T70lEAhw9OjRZy7uQqbJ0uk0169fJxgMPlUC8nZkrTjGOnMPkOi4hRRzog5yrptNgbBBxSNIMa/SksXMhAPa1YgkGN997ve5BFZ2BKh85vU8uaSioqJ8R9rU1BTj4+N0d3cTCATybbTzS1gvyrjofCy8caVSKXp6ehgbG2NsbOy5hS4+1QT3SNTe3l6wG7dUV1ohwgyFXkDe0InWmqamplWVVeZbZM+19+bdF4OUctmJOq9Dbs+ePfn2TK8E5MWM3lTZUosJZXoYrW2kcp98g9KHyaTNnXd8FuaParugtQ8hcqAkOm0hVG5OICIIWlA+8S301n0rfiaWZT1VwvLGRxdOxXnrqtYT6+lOe0IX0WiUmpoaKisrn1vowls99aJg3ZNsO3fuXFWSaTGCezXp5YQZVlpf5MFz75uamkin02uSQVZKMTExwb1791Z07Zdz0YeHh3n48CGHDx8mFArlZabna50vnCp7xrpPXkG4McjYIHgyQaazC9I0PojnoBy0shE46EwQEcwhEjGEA1qYm6YWNiKXoTh2kdk1EDIcDhMOh5+ZihsdHc0rpaw0FVcoNqoXfTGdt7UIXWQymRdqV/wn2osOzxLcI9JKLn4hbaePHz+mt7c379739/evaT/Zo0ePmJ6eLsgzWSpm7+rqIh6P53vSl7oJLLzIFlr3vdP/kYirECpt3HKpTJLN502YzB3ItSFtlh2g/ICDSKTMEsI0aFcgZBotBQIXoTP4swOLndKqMP/8I5EI6XQav9+/7FTcavBxuv1rEbpYbSvtwMAA9fX177D8JNmaFg/CC0Tw+cIMS9WkF3veYtBaP5XBz49XSrkqLTdv+snbUFrIF7fQgs+vkR8+fHjVF+dC6x745t+AtATbebLXwiWfQMtbdG8OPMXcMsIkwp0LHTKA8IPIgGODT4MFdm7WuO3rBK01Pp9vxam41ZTiNsKCF1ImK0To4s0338Tn8616Oo2VJ8nWtHgQNqhMthp4umzXrl3LCzMU8gUuZcG9Om5ZWdkzhFqNlpuXCAsEAjQ2NhZ8US0Ws+/YsYOtW7cW9PzlIIXA744gsgrh50kjizf/PTdFBkAmCxaoXBG4c7G5BK0EwtFo22+knrQfodKmxIZDcORNKP+R5z5XeNbaLjYVNzk5uapS3Ea66KvBQqGLTCbDpUuXePz4MYcOHeKVV17hl3955bL13I3vKiw7SbamxYPwAlhwx3Ho7Oxk9+7dqxJmWMyCLxRVWIhCp8m8FUZ79+7l8ePHq6pre3Vzr/210HJcQcceegdyWbBs4EmRW4u5DLpr2laB/HpwEU+a3nOJeXw2BCo5b6uVhZCukXJCExh5B1rXj+DLEcfn8z3lAhdSivukLPhKCAQC/PiP/zi/8Ru/wdWrV/MeymqwzCTZmhYPwidM8MHBQUZHR6mrq1u16spCsj569Ij+/v5ld5QVQvD+/n6GhobyYcLo6Oiq4nYhRN51KyTUWA18D3/X7Bizhek6tQJIN4OmCKFmEVltVsg/dT5zW0XnckIq7cMSIOasupgTa9M6iBYaX7xz3SbLVjtNVkgpznXdFzYzn0ql8iKNqx04WWGSbLE3XJDV+URc9PnyTE1NTWt6Hc+CeyKGXn17rVJP80dF5+86W43Ao1KKvr4+UqkUp0+fXrf2Rw/W+CVQuTyR3WwEOZoBS4MCVRJBBkxP6jzV5KcJnshhAaisce3FnPues9Eo7NTgJ0LwhViqFJdOp7l06dKaBCqXwnptNVmroupKk2SscfEgfAIWfKEww/Dw8LLa6EtBSpn/smtqatixY0dBOueLkdVbZLjYcQrtTPPWDofDYQKBwLqTG60hMYbQObS0YByzJliCcE0MzbSF9gUR/jTYc5tNXJjPWMvBPGdOo00EXIiBnnKRKosozsJ0F5SuOGpcwCmvX8bbK8UNDw9z9OjRNQlULoXn3Q3uYY196LDCJBlrXDwIHzPBFxs8sW2bRGLhJMTKyGazDA0N0dbWVrAU8mIEn9944kkJrfSchfBq7S0tLQQCgYLjr9UQQPa/AY6Z5dY6jMjGjEQy5BtWRCYNIw7UgvDWFjkmHtdaIoRCZEH7bYR0ULIISRwmQZBDqCxqRpK88mvol/5uQV116/X+VoO1CFQuh/Vy0ROJxKrFKM6fPw8rT5KtafEgfEwuujcrPTMz89wbRr2JsomJCRoaGlalcz5/mgwKj9uXs+DewElbWxvRaJTZ2dmCLL5SCqWUiSn7HhD+T/8PbutedNth1O5WWNBtZ/f9b4RjrLdIZIxai/X0NhPhZsEHKlaC9M2YJ2o/iCzkQmhfApFJol0JERCpHFoFwXURMmuOofyEpm5xq8CuuuXwca3xKUSgcqVS3Hqc51os+JkzZ9BarzRJtqbFg/AxWHDPdV1KmGEh6ZaD67rcvn0bKSU7duxY9QXnNaF4bbCZTGbFuH2p3nI9ty54fHz8qZvWSr3o3qSdUspYlyuXCP/pPwaAvLMbfvs30SUlEAiiGpugbgdu2yGsR9eMiGJCznWt+YEsWkmENfd6nhGKORAFfKAJIMhCLIeIYLrfCCAw46ViPI0SxUhmQYG2/YTjQxw4cGDlrroCRDA+DoLPx3qU4taKF23pAWwwwaenp7l9+/aSZSswBC9kG4oXu3vSyoODg6vaojL/ta5cuUJ5eXnB+8kWEta70ViW9czAyXIxu0duL6ljvfF1An/hz6C216NLShCALi5BjIyg/X78H76Hqt6CVV2Fdfy2GRZJKkQEFJZJouVMk4o5We8EJUwDVYAy70/ksuicaVnFMxgiBxkQMgs2uMJGaIWIT4KTQdqBZ7rqJiYmCu6Z/yQIvhCFlOIcx1kXN/1FGxWFDZwH94QZlnN/oTAXfXx8nPv37z8Vu1uWle/lLhTeSpq2traCxRAXWmQvIVdbW7uo7vVSBNda5zv2pJT4fvXf4fvZn8I9egyRTkMkAraNnJlGVVXiu3IJt7kFcjkoGQJbQUJgnHGNyKTAD8Kdc9MViDmOiUzmiXqLM3cuFmgniCBuknIayAXBSYI/Y6y3DCAwwyj29d/Aaf+xp95DMBh8as56Jev+IhB8PhaW4jxhx+HhYS5fvrzkVFyheNH2ksEGENxrzfSUUlZyo1dqOX3w4AETExPP9IGvNnb3xB6LiopWpXQ6v/ttpYQcLE7wp8itNb6f+dvYX/0V3JdfRY6Pomq3mvg6k0WXlmD19eIeOox15xbuwcPI6GWTDRca3DmL7Bkb76Vy5C25yGXNppN0ESIdezIHnjZ9rF5tXKeNfJvwa/N8WxshiFwO+94fPEPwhZ/LStbd69deL6y3kIYn7BgMBjl27NiSU3GFluJetFFR2CAXvaampuDGlaWI6jgON2/eJBQKLTp3XWh9ev6gR3t7O9euXSvsTSx4HW8ufSWPZGFSzqvVCyGQs7ME/uKPIj88jzr1MvLxEKp+ByIeR0ciQA4ScdT2OuSjQdx9+5H9D6A5DRmJCChTIpvXxZYneJYn46GeBz4RR8w/1XkrhMmCSKaePN8B6csaAUdATnSu6nNazLp3dXXR3d3N4ODgqmL3pbDRa4uWmoortBSXTCZXLSO20Vh3gvt8vuduOV2oc1bo8xbCW/VbXFzM4cOHAdY0TTY+Pp6Xiip0N9n8ZJoQAtl5H/9f/jHE5CSqZSdifAxdVYmYnjRxdyaFDocQ2kFbArV1K7gOuqEc6es1y20CmLoXGEUXtFFVRSPmvy2P4ApD6jnHR8RmYc7xUDKKcGMmK49rSC0ddA60HUI4cZgegNLVSwB71r2kpITa2lr8fv+qYvel8HH2oa+lFLfaMtmP/diP8bWvfY2xsbFbWuv9C38vhHiVNewjm49PvBd94Z1wZGSEnp6eFRVcVrLg8+fA17p80HEcHjx4gNaaI0eOFGR55mfqPXJbX/99fP/sH6OjJYigHx0OQXzGxKhOzpA2lwEVMG64bYHQ6HAQUfzQ1LDnYm1s0C4Ql+ikMr8LuWh3rh9dPymbkSNPcJ01zS06AyIAYjJu/i2UGVKxfUAOkQMd9AEu9pX/jPPZn1nTZwdPrONqY/elsBExfaF96CuV4gYHB/NioIXiz/yZP8Nf+St/haNHjy73sFXvI5uPDSH4WoQHPVc6FostqeAyH8tZ8NHRUbq7u59L5imVSnH9+vV8nX01F1Y6nSYWixERAv8vfwX5zpvo+h2IsRFUcwtifBRdtwMxOY7eUgOzU+jKMkgm0UURs9vbb5nGlMCkWQFsGSujcwI5rdGWa7rRHNe0ovrF3O+fxOfaEUajDcDbOZiUEFCQEeDXZjNKBnQgZxJtWkBOI4SD78F5VleneBqLEbKQ2H0p6/6iTJItVoobHh7m1q1bnDt3jvb2dv7lv/yXK+oHnj17locPHz7H2a+MT9yCg7kQrly5QklJScGWcjEL7okhTk1N0d7evuaNF9402b59+1BKMTIyUvD7kFLS0NDA429+jS3/7ddxnSz+ogg+N4euKAcUlJeBUOiSKEgFoYAxy5Ywmua5HPj8kJhCVGcQUqNtFxUzooraNbprWmMSawKYck18LSyw5m58rgWWg05iYnSA3Nxn6/KkRz0DIgQ6CXoShJoFy0LOXkPevYVqfcZ7LPjzWOm7XMm6z89qvygEXwifz8eXv/xlfv/3f5//8//8P/OjpOuEVe8jm49PnOCzs7MkEolV7/FaaMHnJ+WOHDmy5i9tYGCAR48e5SfBpqenC4rb88m0TIb6//prWBfOoQM+ck6WbDZFKjaFZfvwWQJfwIdQlmFo2jUuuiNM/I2LcHNofIjgEOQU2g8okEkgNBeGBzC95HP3MGGBngaCGrxErmtq5GLWtKpiKxCuSaxpZdz3EPlEm5gElEZI0GkXkUsS+Ed/DlW2F6qrcT73/bgvnYEC4+a1KKAuZ93D4TDZbHZddczXWxM9Go2ueYBqEaxpH9l8fKIuurc3OxQKrTr7OF+dZT2EFRZ2t3lf+krvZX4yzfrwfXy/9C8Bjc5l0AEbv3LwlRQj3BxuNII79piYLMY3OYauriGYmEFUb0HOTKLLKxHjj1EVVVidN2F/HJQ2CbbMk8TaU0ouHuw513xaPSG4dp88bi4LJyzQs+ax2p2rqksgCSIFIjz3NC9BVzKD75030D4/4uED/L/4T1F79uH8wA/iHj81N7Wy9GfzPDHzQus+PDzM4OAgV69efca6vwiZ+fUuk611H9l8fCIWfCGZrly5UpCA4nx4La5eE8zzCCvkcjlu3LhBWVnZM91tyyXzPHLrBz0EvvrvEXduQSSMiMfQFeWI2Wl0dVU+1rZGHyMbG/GPDOPuaEQ/6iNZVoGv8zbOlm0UXf0A3dqG79y3cQ4cRqpr4BfolED6FVrNNbl4WfI5/mpnXl08BzomEVH15P9sjDvuyTpNWwi/+6TU5gcxPU/uKQtICVpgJYZR4Rrc3XuRfb3IkceQzeK/extt26ijJ8j90B9HNz67Q309k2JSSoqKiigtLWX37t2r7qpbCutpwddbUVWscR/ZfHzsBM9kMty4cYOqqqo8mQpVSF2IdDpNb2/vqmSa4ekLz8u2Nzc3U1NT88xjlyK41hp3ahLfL/1zrCsfoYuKjIstNNpnGTMYDkEmjS6OmkRaSQniUT+UlmHfu4mu3UpJx0XUwXas6xeJH2gndOMik62HCTBAeI6QMq1MfD13Hp6CS37mO0u+FKYFMK3QIRAeoR0gQb5OLtIYUtuYfYUW6JQAnzYZdY25lKQLfona3YzV24VIp3B37QXXRT7oRvr8iFgMeekCQkpyP/jDOF/+w6Yrj/UlODxtbVcTuy93DutpwXO53KryPj/yIz/Cu+++C7BbCDEI/Cxzt+Hn2Uc2Hxvmoi8GL3m1sBNstV1p3l4xpVTBGm7zz8278BZOgi33+PnQ2Sz81q8T+J+/CWhIxhFOBnw+xOwUOlqMyKRRoTAyMWsGQlJJxNQEODnE8CC6fgfy9jXUsdOIoX6cV76PUHyW3MGXEJNjhALGCxNxjbbmXHPv29LeDUcCCpy5zDiYTjUb9LgF4bnPNC2elNngiSAjmMtGSUTG3ES0A0hh4nUXsDT22CU0W1C79iKHBhGzs6jabahwBDEyjNXfBz4/1td+D/+//heoYy+R/f/8CXTg+WWS56PQmrVn3Xt6ep7qSFvMuq+nBV+toup/+2//zfvroiWjte4jm4+PxYJrrRkYGGBoaIgjR448oym9GoKnUilu3LjBtm3bSCQSq777erF7X18fY2NjK245ecaCf/P3sP/bf4JMBkpLEP0P0ZXVWJ23cZt3m2GNTAYdCSOnJtHhCNbta+iyCnMjUC5q5x5EJo376hcRyQRu0x7koz7U9BQol7LZcfTLtknABfVTrvRTPz2iO3OWXfPkG43r/CYTcnPWea6dVVhP6uHaBTLKjJ+CKZ+F5v6RMwfVVQJV2YJ15SKEQrit+5H9DxED/aiaWtyarTAzjX3xAiKVxG3eReAf/j3aHw+hf/j/i/7Cl1BHjy8brxeCQj2CTyIzv9H76NaKDSe467rcuWMEIpfqTS+U4J6Q4d69eykrK8svk1sNpJTcvn0bn89XkPX3Wk91513sn/9pZNddVPNu4yLHZtA1NVgX3kedPIucmkBVVINtIUaG0WUV2B+8g7urFZFOoWu3o0rLENkMqrwSee8muiiKPdiPdh0yW+sIWRJnXxEW/SYLDgifNnH2nIUVXjSy4L4kMuQnykRG5WNxssrsL5uV4FfmMUkLAq6x1ugnNkRr4wFkmVuOYCGCDvbbb6O2NaEamrHeexuCIdz2l5DdndiPBtFVNajWfWgnh+/cdxCpFDMHjxD96APku2+CEKi2wzhf/AHc4ychWrzq724tZFzJuluWRTQaXRdL7i15eJGwoQSfb223b9++5JtfieCeB+DJ9axVyDCTyTA7O0tDQwNNTU2FdaalU2z9n7+O79L7iHQaVVmNiM2gioqRsVnk/Tuooy8hxkdRW+sgk0Y+7EZv2WrIfew05DKG+IEgsusOum4H9rvfQtVuRyQSJGu2ka2qIepk0UJixS6YeDkg0a4yU2JxAbMmo+6VzfBJdLl6EpNnn5BfawExDdF8Ah2RnUvUCQ25uc87gVmc4DdNMqZhRpvlCH7LiEEIcNur0e4W7Le/jWpqQTU0Y7/1OjoQwD39CrLzHtbtW+iyclTTTmYDQUrvdGCl0+TaT2D1dqPKKvH9639JsKcL98wrqLoduKfPotoOocsXH96Zj/Wwtgute2dnJ6lU6rkz8/N72l8kbFgMvtiI51JYjuBKKe7cuYPWmvb29mfusoW6bbOzs9y8eZOioiK2bNlSGLnPvUXw1/4N5b3duOXl2NO9ZhijtNS4qJkMekstJOOobfWI8VHkYB9q63Zk9z2cV78PMmlEMobQYF34Du6RE8j7t3APv4QKBkmNjSJLSii7fA5V34RID8EOjJX1mVq1HgUyhqy4wjS+OAKhFHoERNjE4soXxSIGgOsPYMk0Oi3NlBga7QOypoMN/5xLn9VPhlCSIErMv7WSMDkn+xQQCGsKeeEj3JNnQGnsN1/H3d+G3rIF+8030JaF++rnEQ97UA96KQqFkf4gzsuvmbg9HscaGsDq6iR37CXEzQ7k+DjyxlWsh724+9vAsnFeOgXb6nB37UbXN8K8xOt6N7pIKfMjolVVVauK3RcimUyuacR0o7EhBPdWBhWa3V6K4N6Osi1btlBfX/8MKb2+75U+fO98Dh06RE9Pz8qNK8kE9j/5GbjbAa7CFynCvnGF6d37sC0foft3ERWVEJ9FV1SiyiuwbneAk0NX1SC0wj1yAjEyhHzYi2psRt65jvPy50zyrXkPrlL4P3gX59BxQlc+RLW0osvKgQmEdEA5xkrPgIiAtuZEFLWxsNpfjNAzppPNKkFVbYFQOXS9bz6bshLEbJqpyjZKMjeMHJMFOmPcdGGBjs29X3eeQCOgYjYi7iB8oLMCUhrhy6HadoLSyM77uCdOgrSwv/0G7tFj6JJS7De/hQbSJ88QdHLox4+xPjyPmJ3BOfsZRCqJc+os8vEQVm83ubJyfFcukjt4COv8e6jmXdivfx3hutjRKLKnE9W6z8T1rQcoE4JcTS1W6z5UzRZ0TS1Enq/uPN81f57M/Is4Cw4bRPDKykrKy8sLvtsuRnBPDWa52WuvFr4UwT0tuNnZ2Xx/+0pDKqLjCvZP/SV0KIzy+7GGHhAaH8V9+TWiySTy3i0SdY34hgbJbt2OHYoQ/s630bXb0ZVVUBRFB0JY598BQLUeAKVwTr2GHB1G9naSa9qF/95NMu2n8OeyuCdfRUuJdf8KekcCIR2TzZ7hSYw8R3C1+yDWo2volgPQdQ5dXIOcGoEeB9Vq4kztCyNnTXttcGYGpA0qhxaSqe17KZ/pMI9L+JAyh/YK6vaccMSUYxpeHEyTjbQg7SKCo+h0GVpKcBVicgp3X5s5vXfeZmZnK77aLUS+8w7asozrvWUrKhzEvnzJjMIeOIjV3Y3z8ivIZBLnzKtmpVJRFG1Z+M6/R+7ocXxvfgt36zbs8++hQ2FIJChKJlAVlfhf/wN0MIiITaNtG4JBc7Pfuh3CEVTTLnRJKe6OZvSORiguXfL7Xq/M/FrUXAqYJhOscSeZhw0huM/nK3hFEDxL8MHBQQYGBhbNuC/2vMUGU7zW1XA4/FR/+5IE1xrrP/4S1r/5BVTzblQkin3+HQgGcU9/BjExhtV1D7W3jUgygW7ZjYzHCXz4LrMNO1HhMAGlsV2N/f6bUFqOu++QGQPN5bAvnYdUkvTBdpyZafTLn8NOxBBjI5BNYd2/jTpYjnAddAhISERQoTHz33pLFWLqEUI/LTGtozUwNYJIzsCsIbcu2Yac6QJA2VH0tsMwcBFd3kDZcD/aZyOUw0ztLspGbpvtwVlMIm4U0xOPnmuOmWt5swVkE8jLF1HHToAvgOjtQTc0oSyb1LY6QraF7733UC27cFt2Yp8/B2hEURHkHNzPfwmRiOMcaUd2dyGHh3Dbj2FduYR7/CXsjms4bUcQjkPu+Cl0OIxKJlDRYqyhAdxwmHDHVSgpxZqdgmwWVVeP1dOJc+I0vm/8Hu7Bo1hXLqCjJdg+Gzk7jVu7DZlK4hx/GbVnP+7Rk3nLX2hybTnrHovF+Na3vrVqFZsCpsnWvJPMw/p27q8R85cY3Llzh4mJCY4fP74suWFpsqZSKS5dukR1dTW7d+9euTNtYgz7z/1hrF/+eVRrGzoUwfedN9Db63GPnkI+7DaxdfMuEBIdCCH7H+C/ewPn7OcIb6mleOAh2h9A3rzCzM69zO7ci7xxCSbGsHo60eVVxF56hWwqRbCoCP/1S8iOq+iKSsTgAM6pM6BmTAdZQoIl0UJCUKB9IcS00bkXU4/MOWeT5mfgidWwem+gKg+jgqX5/wtGSmDWjJLpSDUiMY2qNTvAiwOhuVo1uG7EjJqmmFdSs8B2TRigNNJ1cE9uh5kY1rn3UHsP4FRUYV+7StB1sfxBY7WjJdjf/AakUrjtJ9A7GnH3H8B65y3k+XOIbA4xPIzz6mdBg/PKZ9H+AGr3XogUYT18gJQS/7tvIhyHwNuvIxIJwndvI7I5sG20tHGPngCfn9xnvwiAc7AdXAcyachl8Z17C21Z+F//A+S929jvvYH/V/8Vwb/2pwn96Jfx/fq/wRoZXnNmvqWlhWPHjrF//35s2+batWscPnyY3/iN3yjoOGfPnl1pKCW/k0xrfQEoFUIsLpCwBD7xYRN4oq925coVKisraW1tLeguuJhrP38SbLHk3jMEv/A+/p/806Bc3MMnYHwU++Y1nOOnTdb71lV0aTm6tg5CYciksB4+QIcjuGc/b8j/oBv3pVcIJOKw/zBFSmFffJ/k1npyWiB2NJPz+Sm7+D7UbkOMJgCBe/wMIpPCPXYaMfohwq8gO1fKCgp0TTNyrAu3uglr7D66qByRnARAzI56bwgAHalATE8gO84zufs4+aAmk8Z6cBPV0gBqzlLJIrTlQ/bdRe1oxeq/jNh9EDedRAzfR4oUKiuetLNmBDgabUusxCB0S5xXX0MNDhK4dpXcgYOIigrkd95F72gANOrkaVQ4jO/bb6CjxejSYigtwz18BBGL4Z44iex7iBgaRDc2Ie/fxT11Bvv8e7hH2pEXL+Du2Y+cmsTdf4iUAFUUJbhli+moa2rBmp1GRUuRA33ImWkIBpADD3GOn8G+dA5n3yFEJoNz7Aw6EsG6eQUdLcUaHUIHQ1jhInb91q8TKK/E/cM/gvN9PwhFq281ra6u5jOf+QzBYJCf+7mfY3Z24eahNWPNO8k8bIgFX225IJ1OMzQ0RENDA42NjYUvA1hA1sHBQe7fv8/Ro0eXzNznn+M4yP/wFfx/6oton41z4CiitxPr1jWcV76AEBLr2iVU026IRI01cxWy6z6qsQV39z6s999GPOjBeeUL4DqIx4OQyWDdvIZz6lX8jc1Ek3Fc16Xi2gUyZZUkAmFSdY1kj5wwtfFIMWKkG6HjkLUgGkHt2GsmyorN2h6KTByoSswgjbZ8iJm57zhrXHZVbEQthIawCqItE7aI6THznMh2iM2Yz2DgPrpqNyKdyI+jibGHyOk4otbcGoSwTatrRpiBFz9mtNUG53AturMTOTBA7vRZpGVhv/02etceVEMjpNIQj2O/9Ra6qhr3pZeMN9TUgvXWm8hrV5GPHyN7ekwDzPQMzmtfgEwG59XPmTj6pVPo2lpUZRWxSBEilSJUuxX7QQ/CyRE49y5W30Ps3vsmdKqpRYw8Jnv6NUincE6+CiWlMDMDaHzvvm46+CJFuPsP4+w7hP3eG0R77oEQ+H73Nwn9qR8g8E9+GnnvVkHX33x4k2R+vz+/bmkdsBgRPvlW1dXAmygrLy9f9UTZfNd+sUmwxSClRAwPYv+j/wPrG7+D29aOLi5F9nUjclnU8TOIZAIx+BC1Zz9COWajZ1Exsuc+6sARtD+A753XUVU1uG1HsLruI8Yeo3buAa2NWz88hOztZLq1jaLxUZyTr+Lz+bAnRkkHAljXLzK9fQeR0UF8ReNgWejSEmR6FKJl6MkwKlgG20+gVRTlq0cFtkJ9CfgEVvc5AETctLQ6vkhePTkwOYyqO4EcvIwY7zfvu/fOXF0MxOw4bsNhJLcQQ93migmVwEQaOT6Iqt5ulijMjpoMu+0aomcV2gIr/oh0rgL27iMwPo7se4j7yqtox8F++y1UUzO4CvXSSWM5334TXbMFmU6BtFHHjkMyRW7XHmRsFtXQhJiehFgMAabysGsX1ofniLcdpvTGBdyjx5Dvfhtd34h1/y5IC11djYjHcduOQnwW5/hpRCKGmBxD57LYd2/itB1B3ruFc+AourwS2deDSCXxDQ1AKMLYzn2UxWNY926hDhzGev9N5M2rqLoGcj/2V1G79xV0La5lq0kBWPNOMg+fGMG11nR2dpJMJtm/f/+au9IymQxXr15ddBJsMUSuf0T1r/wLZPcd3DOvobNZxNQ4VNXA2AgIMadzBigXbdvmwpydwd29D5HOYF94H/fwcXRpOfb5dyEYQu0/BAhIxpBddxCT44wdPE6xbSGDIUglkBevotpPEr14DvfEGcrudpDb5keh0eEAmcptREZiKLsImfZh9T9EDtzBbT6EHOlHl9dh3T2P23IMt+QAIpBGjBkCJ1NpI9kmBGJ8EDE1gtrZijVtRCZ1sBRdXI3Ve8F8EHPCD3J2FLdhD/jKoVzA2AN0ZQOy8xyqdieqqBrr4oW5erlRiNGAf4vCHX5M1smR2bWbwMgoobt3cI+fMG26g4Oooij2G6+jqqvR27ab6bPSEuTwMNpnY127jMhk0aXFiEeDqIOHTcLt9BmsD88R299G5E4HzsHDyOEhU5Hw2ejycnRZGSKdxt2yDTk6DOkU1uw0YnwU5/hprJtXyR4/bb7PhhZ0OIL/3ddRZRXoyircfYdRoTDlH72Hlcsai+/kEK4LmQxidga1o/C57g1aerDmnWQePhEXPZfLceXKFSzL4tChQwQCgVUNm3hwXZfOzk7q6upobm5eeYnBf/+PbP35v4019hh18Bikksb1tizkvZuo+kZ0KIR1+wq6dhu6tAzZ2wnJBGpnK2JsBNlxCefVL0A4jPXet1E7mnH3HIBkEuIx5EA/biTK+KETlOYy+EaGEZNjyOuXUO0nTQPM2c8i+npQTbXYAYFs2Yc/PYsVLcNREbjxASRi8LjXnPfEXGItPbfDzRfA6rkJI2lSNW0ARP3mXq1LtyIyKURyFi2fJHB0cTXMxs3fbT+y67ZJ4gFEqpFdtxDxafPvlCmQq6omrKsfoUsrULvaUCE/+AQCgcxM42MGu62NolSKQG8vk8ePM5PNwqVL5EpKkffv4Z44iWo7aG6WAR/2228hJiexHvQiZmfR9XWIWBz3zFnEzDTO574A09NMHjuJpRXq1FlEOIxuakFXV0MwjK6sRAw/QmfT2O+/hezpMjfkeBzn1S8gtDbZ9FwO38XzCNuP78L7uM17cPYfNsq0E2P4Ln+Atmxyn/l+RHwW6+pFdFkFuq6R9Ff+IwSXT/LORzKZXPWo6I/8yI9w8uRJmJsmE0L8OSHETwghfmLuId8AejE7yX4V+MuregE+AQsei8W4efPmU+OZtm2vmuDj4+OMjIxQX1+/6JjnU1AK6ys/i/Xv/jnZ5j2Qy2BJC0pKEKPDiNlZ3L1tiHQSnc3iHj5u3PTYDLplN+RcRGzaDIi8dBaSCeT927gnXzF3/clxKClF9nSS2tFE0heg/FEfwucD5ULONaW22Wnc/YcQs5OohlqEmEWmMrhFxbiNJwh1fIgu34Icf4CzpRF77AHZUAn+WRNH4xHdNSXIdLCY8I1LOIdOIWPmd7qkBh6bv4tkFlVeh5wcAG1jPbiFam4Gfwh5+xburoNYgzcgk0XEZyAZR0fLkP13UNWNiEcjCMdBbWvBHejGJzLomjqYnEHEkqiiLHJiDDnyGPczr1KcSCB7e8ns2w+jo0zV70AKiHR1IWprsd97D71tu7HCFZWovXsRqSSq/RgiNou7cxciPst0pIigVvgrq80UnVLoSAAx2I+u2Yr19hvoxmZk1z3w+XEPHkVks6i6BtNN2P8QVbcD6/YN8zvHQR07jQ4G8b/3Jqq0HF1eQW73fmYcl8oP30WkkuROvoKuqCLzD/+VSaiuAolEYtWNLgVMk615J5mHDSuTLWZNR0ZGuHnzJm1tbU+Rcr46y0rwdoL19vZSX1+/cqecUth/+89j/7t/jjrxMpZyiJeUE3NdsuOjqLJKRDpmVvsEAuD3IRKz6OJidCg4N8ddjC4uMbu8chkIBlBbtyHisyblUVyM7O0k1tBMxh+kfLgfttQa2eOqLajmnQjloKPFRoutugTKAujaWtymQ1h3r4GbQ2RT6EqTSBPlphpibW0EwA2XYM0RPTvx2PyfZd67dfsKOrwNLQT45l2YsRl0ZE7hZtrE6jqyBR0oNf9nzT02o1A1zWbSbVsrws2hShqxuq6jqupJT4zinxlFNexHb62Hhnp0YyMiGoFEL+7ePSaGBvSuXQTGR7F3tlBSXkZ0YpxMdTX2pY+I7dxJJhzCKS9HbdtqWm0rqxDJJKq4GBGPMZPNEUQRsC2zi21sFO0PYF38EMorsc+9jd6+A6RENe3COXnWfCdODuvKBayb103uJJnEOfMaOhhGdt017bXfeRO1vQF16BgEguSGBinvvW9K/q98H1paJP7W/4UTiqxaXvtFXHoAH+O46MKOsvkodJpMKcXt27cRQtDe3s7Q0NDyz9Ma+6f+EvLiedxXvwgzU+gdzZRkM7i5HFlsVF83s9sbiI4Ooptbsf0+M2Qh59pDa2rATUMgjNq5C+IJU6OuroGcYy7SVILpugZs5RK1QO3cg3WnA3ffQcT0BIRCaCeDLi5Chy2kSiNmphCJJLqsEpGYeeIOzumeI+aSpcG5i6amAR7cQAuJf8qEYe7c6qZs6VYC1z/APXoGnCcNRnKkD5IzuK37kD1mkYHsuoHaYmJL+aADXVyB7LqF2nkURnogncTdcRAxZUo9s5FySvpuoCq3QbgI+aADHFAVdYboykYncqhgNTKZQ6eS6JYmZDyG9vuhsYHIzQ7cl18mMjxEuqKCbCaDMzWFFQkTnJnGKpkjt5MjhMZv+UArxMAgurIC+/y7uC+dRt7pwD163Eg8S4mOFGFdvwyWlZeycj5r3G1t+5CPh5Dd93D3tkEyjnviDDoQwn7ndZxgGN/W7bC93vz9zg3SX/1fyNrt+U003jYay7KMtv0y9fIXcfEgfAwE91YZRSKRJRVTCymLeTvB5velz18rtBis//sfI/p7cXa1GhmlqJkAU8EANi6WELC9jkAqTnb7DtzRR8xESwhqF7skiM8SIHxoaZpOhBAQ8oOrTPumBpXOMKscgj4fgUgU7RrSu0faIZtB1deDdNH+EFLkELExRGoWSsqQPddRPsskxpKzuC1H0b4w7p6X0XYQd/cZtD+CjlZAyMR3qWgN4bkSWfGcVKoTLiEAWFfOEdv/ElFAlW9DDs257aLY6K8DuC7aMtJWIp3A2fcy9uC3EI8fAiD7bqOqW5G9HeSCRUQSkwitUbWNiLFHuHWHzGKVdBrSaSgpQxT5jZZjNIpOhhCZLFSUIWZm0dpFHTqEePgQ3dhIcOgRumYL5DLklCJjS9TIY7I+H2Enh11U9EQco7gY69olnDOvIkeGcQ8eMWSenoRoCfZb30QXl6B2NKIDQVRxMdadDsTwIOrIcVPCPP0qKI3VcQW1Zz/2hW+TqjbiFcGRIUTfA2zlkv3Hvwh72/K+slIqT3TPmjuOg5Qy/2c+XsTFg7CBBBdCEI/HuXHjxrIbSgqBNwm2e/fup2qMy1l++ev/Fus//Rvc1gNYNy6hdrVidd9FNTRhjQ6buml8Fl1ehpSKQCiMCAcISUHWdcnFp0ii8UnwhSP4fDZC2sZ6WwKhMjhaMK0TRGuK8QfmLLCyjWipACF8oHOIbAbhWlg3P0RvrYPYFDoUwjn8OURsAt1wGOvOVdT2FqxbV9BFpZCcMcSqqIXJCZwtAWa3HCAUCcHMsLnhjJoMejBsLixt+4ncvk6ivoWc9FGKF7NL3IaDWA9voLbuRg73PZF7yhlPQY4N4DbvB18AlfNjOVncloMEuy7h7mg1aqvDk4jQVsRwP2JkCF3bAP4kiCzaDoLPQvg1QlhGY70ohHZdhBToHXXGg9q+HdnXg6qrx59KYVkWs6EgUaVwo1Gyj4eQqSRWIID/QTfOq58z47m1WxFaIzuuoXbuQX7wHVRjC6q+ARGPQSqJfbfDuOavfB6RSaFqapEjI8iuu7g796CzOaZ378dfXEzk3Nsmq992BOfzP4D7pR98+vqZI7BXcvW2wno/vetOSnPjf1Fd9A2LwcfGxrh+/Tr79+9/LnI/fvyYW7ducejQoWcaCJZqVZXf/hr2P/hbqMadWJfPoxtasK58gK7bgXXpffPFD/Whq6oQ2YwZZrAtQyqfTWikj0hZCeXTwwTCIfw914lNTZAYHyI7OwGZOGpskPRoP2U2BIRCTg8jktMIN40c7cWaHUOkY8j4tInr49O4bScRUxPobc1Yty4hsims+9cgaGJhXWGaVdS2ZkPu0irkxLDZOdh9h7J7N7FVGHfXGfSWJkQubd5wai47XtOAzKYJTcxQVP4kxxGfnCQzMYMWAuWLIEf7UY1H0f4Q1tVzqJq5clAwip6Kk31stnX4k7G589qOvHIJtb0FMdyHGB1E7TP9A+JxH2KgCznchxwZREyNQWoGnBTgIHUOnDTYElEUQqQTuA2NyM57uH4/8XiMItvCX1lBZGKMsC0J2hZ2KsH0waPMjo2SSGVwxkaxzr2LatqJGBlGHT+FamzBuvIRYrAPMTaKrmvE+eyXkJ13sS6cg6IoYnoS99Sr6JJyrBtXCLkuoRtXcfcfwvnMFyFchPMX/vqK16GUEp/PRyAQwO/3PzW45Lpufn/Zi4YNI3g6nebYsWMUF69euQOexO2PHj3i2LFji2YoF7PgovMO8tf/He7Ln0P7fLhnPw/ZNO6ZzyL6e1HHX8a6/iG6ZgtifNgk1nyWqYVHi7HuXEKXlGJdP4eurCZ4433E7oOUuimCJaU4/iDWjXM4kyOIkjKs2XGsq2+ji8sQUmDduQjBCGSSJlkVLUUOdkMuh3XpTVNXv3MRHYwgu81UlxieWz01R1QvHtc1OwBwfEFCM2Y6TMxOY310DhXYhqprBUDOWXJdZG6AcmoMHBNaABTHJgg/fkhyexvZR+axuZkYzra9iFQCXWoScSoeQw8NEBkbQG3fhey/h7vnJeTNDvCHjNji6CBq3zHE1Cjy/lV0bT1qXzuqvsnIUkVLTLun32duQOkEwskgZsbRfh/YFnJsBKepmVz/A8JFRdglJcgHPWifDzE7A6Ew6uhxioGSTJpwIkbw3m0mDxwmlkqRiURxlcZ++1tg26iGFlRDE9qysD74DnJo0MziCwvtKnQ8jn3xA3INzVg1teit2wGBfeEc2V/45admzguBlBLLsvD7/fj9fn7lV36F0tLS9Vx2sG7YsFtOfX39qkpf3my3lBLHcbh16xbBYHDZTSfPZN9nprD//B+FRBzKyiEYBBS6rMKosew9iBx5hNrbhnj0EKLFKDeHTGXR5RWIqVFU0x7EUB9qzyHE5CjOqS8hUnF0MIxv4jGBgS4SW5vwbdmBv/M65LJMNbbhHxsm0nMLdegMZFImRhcC3/lvoIvLIRFDRytR9a2IxDZ0eRUiMYP2+ZGjfehIKbLvrvkspky2PIMwsmrbWxAPTPukHDEEFa6LuHMf56XXsG69632KgGl2sW5cQO07jRh/iBw0TUShZAoxPTdG+qiTad9hSgHuX8UNFpFMOISaD8KdD9HRCrTtQ9lRfBMjuIfOYF0/h9qxCzEzAcrFOfElrN5u5LVrqKZWCBchBozUlN7eiIhNoUoqwbIRiTiytxO1sxUVCCCvXcI6chw74Ef2dJpuwkf9qG116JparAfdiJHH6C21iGyG3GvfR8nsDPLmdVJ7D2Bfv8xMyx50eRXFN65AKGRuDC2t6OISrI/OIxJxcifOkBofxX/gCD6fjfX+2+jyCnQgTOaf/Rt03Y6Cr9GF0Frzq7/6q5w/f5533313zZt0NhIvxDQZPJnt9ibB5ssqr/QcwIx7/h8/jpiYhLIK5PgoYnLMjHnevWFi4tkp1JZtps2xuAQtwL55ES00crDHuMiZLGpXG2JmEl1agXXhdeSN8+C4OFMTzLSeIOgL4L/yDtLnh7ZTlA33EunuYHb3MWIT48jbl8haPkTvfdxDL+PuPIicGIFoKdbFt7BuXUT23MG6eg6RySIf9KCqW9Cldbj7XkGHS1D+EHrGSGCLaCkAqmIbYtaUo4jPIpRCTMyith5Eh4sRM3N95zWNiGQMeeU8atuTVksdLkc1HzN/D4SJzukqW5kkk5W7KX54Hz34AI1A3r+Ku+c09sX30KWVyIEutM+PjpSgy6rRxVuxO64gRh6hDr+EmBpDdt9F1zeYJiF/AF1ehXX7BsLNIWanUS27UZPj2Dcu4544hW1J5P276MpqI3nVshtdtQXr3W8bQcede4wEVH0D1s0bWJcuoNoOE0Cg9x0kUl5J6Ufvo4RgcnsjyUAQ3f8AefcWSIvsq58n82iA4p5OfKGwabw5dhp3/2H07lbU9//gmq9XrTVf/epXeeONN/gf/+N/vJDkhheM4BMTE1y9epU9e/asuLgNnlhwrTX8P78Eo6M4x06iyytx9x1E12wFKVBlVVgdl0zHWs8ddHGZIXMkimrei3XzEjpUhPXB62Za7MK3EROjiKF+dG0TuZNfIjsxgi8xQzQdRw734Rw+i67chn3h2yBtnNM/QPGdjyh9eIf0kVfRsRnk1Cjx2Vnsi+/g1rfi1jZBOIq7+xBy6AFaSmTfffNmfD5kXzckElhXLuJmbALRatzmQ2b8EdBV281PIZCeW+8PGtXWwJb841TZ3OCJUhDLoCNzCyEUyM5b6Gg5qmE/1u2LZGpbACiNFIPlIzD1mGRjG6loBbHxGUQuQ65+F2JiBLX3GNhBxMgE4vEAzE6h9h1G3rwIoTDuweMwMwtI5L2biLHH6MoKtD8AUuImE8i7t8idPIudy2LdvIZuaDQVjrodaJ8f+61vQiiCe+osKIV40IN8NIhIJnA+83nIOcjrl8H2Yd25iXv8NBw/ReWNy0Qe9pLZ3sBM9VaSoTDq2hWKBvpwXnrFtOE7pgQpO++R/bl/9VzX63/9r/+V3//93+d3f/d316wR+HHgY210WQ6ZTIbu7m6OHDmSV9FYCXkL/u638f/DnwLbwvf2N0EKZG8XuqwMVb0VqqpxD71kVExr67Df+yYikzQEjhShdrWhSytRB15CV9ai6lpQzfuQkyMowD7/TYKTo7gnvgi5LCKdNC7y2BBu61FU4158730d3XIA58BLBG98QPjhPZyjrxKeay9NuBrfuW+hUxmcokrcljbUrkOIqTFD2H5Tp45nTDlLbm3EvnYB68Z1iGVx9540CwkBXWsstDnwXIJN+BGTSdS2nQj3SeJR9t5DVe9ESwvZex+RmEFt2ZPPnudciQ5HsTsuofYY6x7KJAlEKigZ7MYNF2Pdu0Y2XEwioxBXL6IDQRPqHDyBvP4h7v52VHkt1vnvmKGcjktGTw1QjXuwPnyPbM1W6LyLPnAY28khL1/A3X8IXNfU7oXEunAO91A77olTWO+/g7zdgd5ah95ah2rZjey4gXX1Eu6R42DZ6EiRaWR6+9uona24r32R8L3blN66jhMtxdlaR7x5F7mue9gXzpGr3Y4Wkuwv/ioUr20LDsBv/dZv8d//+3/n937v91bULPik8YlbcK019+7dI5fLsW/fvlV9YEII1GA/4t9+BefMa6bl9KVXkBfPo3c0Yr/7BiKdRN67Ca6Du++oaYY48wWjbtq4G/vie+ALYF94CyzL/BQ+5M1LZHfsZSpcBsEi9LYmrA/fQDzsJvfSl5Ddd5DjjyGZAFehQ0WIx4PIsRHUvmPGOo+NYA/24R56maK5TzpXuZXAhTex7nSQUD4yrS/h7j1mGmKA8NTcJtMicwGq2gas+x1YH30IcQd370nUXJebFtIk8ACipYipccTgSH4DiqregZwYxrp1Gbftc4jYFADywR0yk2Z0tGiwE3f3KUQyjujrQgfC6KIysMPIZAx2tmFlUoi20xRf/ZBE0z6s+x3Emvchbl7COfoK8vEo1o3LqP2HsS6dQx09CbYfVVmL/e63yBxsxx3sx9regMxmsS5+gHvqFUTOQV67jK6oRkxOoA4dRYfC2G98A4qiuK98DjEzjbx+GRGLQyiM8/JnkBNjWB+8hy4qhkwW9+ARdFEx9utfRwfDjB07RWRygqJbNwiVlOMvKSV95Diu4zJUuYXLwSIGBgZIpbydyoXjd3/3d/nqV7/KH/zBH7yQGmwLIVYQbF+zmrvrujjO8tulc7kcHR0dlJSUkMlk2Lp1a8HWW2uNk0ph/eBn8XXfI7WlFmn7sCvKsCIRRDaNjkaxbl1DNTYhb1xENzabbLnEDBVYFjpSZEYjNYjEDOJxP3rHTmRfJ9NNByjtuYkORXCPvIx16yNwsuhQFOHkUPUtiLFh5HA/qrYeVdeMfdFosTlnvh/rxgXE7CTu4TPGkjbsBK2wbn6EKilHxGMIJ8dE80GEm8Mf9lPUfR0A1diKfHAXt+0U1o0P0IGQ6ZxzcjgnP4dIjSPSCWSfkWZydx3FunvFjHpOT6J27QJ/BOv6efP7trMwO4bVf5dY0yF8kxMEZodBa1TDQWRPB8J1cI+cRTzoBmEhJobB9uHuPYp19TK6sgYxNY4uiqKkTSpahZieITL4gNS+I4RvXsE59RpibALCEayLH5A68hK5yQnCgTBydsYIMnzuS8i7d5CDfTinX0W4LvLGVdS+A8iuTtSevWhfAPu9t8xn+ernEMkE8t4ddO1WxIMe3GMnIZPBvnwBVd9getIrKklMTVLS0wVC4Lz2fci+XqzOu2YrqtZkfut/k5rbaTc2NkY2m6WiooLKykpKSkqW7Vb72te+xi/+4i/y9a9/fUWl4FViw/SWPzELnkgkuHTpEtu2baOlpWVVAyfezcP/T/8BwaFBxN42igb78aGQ927j3rpGbqAP7nbgtOxGI3Bf/j50KIradQBicXQwgtV9D+E4WHc6IBRBDvbhtp0ig810yyGiWXOHVy37sc99C1yFe/QziPHHxuJaPnRNPdqy0ZVbsa5ewD1wwljdD98C24ez/wTy/k1jXR8/QsQTqPqd6Oa9CMcIIUaHH1Leewdci9mqnczU70X0dc69WXOTVPU7891osvc+1rXruFXNqMqtxpL3G6LrmnpEKom8ewdtPenTF4O9iKlpcuFi7ESc4MgAqvUkas9R02Sz7yXzQCVAC+TjAROyFJUgcgKRSaNLyk1ZbVsjIlxGOBCmaLCX3NFThG9dYfzQS6grl0krkJc/JHn4OJlEknC0FKuvFzE2Qu4LX8Z683XE8COcz34R2XUf64P3UTtbwfIbEcVUCuvcu7htR8h9/ktY776F9dGHuPsPmhtz7XbE9BT2pQu4+w6iGpoRD3tR9+9SFE+gjp1E7TuA/Z23se7dxTn7WVCK7C/83xAIEAqFqKur48iRIxw7dozS0lKGh4f56KOPuHnzJsPDw890SL7++ut85Stf4Q/+4A/Wm9wbig2z4EqpJdtIPc30AwcO5Ovk3d3dRKPRZSfDnlrV++1v4P+F/wu9vc4kacpKEbFZdMCPHH6Eqx2YHEdLcIuKCU6PoFpasTo7cNvajdu6/4i5uOubEUMPwLLJFZXinxrFOXoG68ZHqAPHIJsx1rG2ATHyCN20BxUIYN29hnBd3P3HzUjoowdowDn9RewPXjdji4fOgBDIa++jDp7EuvYh2rJx288i7l0mUb6F4jkyq6a9yJ47ZPe1o+IzZAMWwYlh/IkZsgdO4b/xAaq8Gjlh5JrcfceQ3bdwTpzFvvI2wsnhth7HunURHYqAC2r3fkhOY/WYZF5yz1FCXTeMuKNl4+49hn3jQ7Q/iGrahezuQTXuwrp/DR2J4u7Yh33tAu7ew1h3r+G0v2wsdCSKde0j3LZ25O1rOKe+gP3+O6jGnciue8w2tJDDwu/zE715Fb2tDrdlL763X0c170RV1WB1XDM3zROnkH19yAfdqMZmdFUNTE+Cz2+sdlW1mRW/+CFiahL35MsQj5lSpCWxrlwiXV6J3t9G8MZVyOXQ1TVG5dZxkP0Pyf3U38f50T+/7DWrtSYejzM2NsbExARCCK5du4ZlWfz6r/863/jGN1YtSlIgPh0WfP4kWHt7+1NNMIVsN/HILXu7Cf6VP296jM+9awj41rfMcS6eR22vQwbDiOOnkXvaEA3NJBr2EE+lSVdsxR0bRfvDiNkYOlKGqtqGu/MQ463HUcWlpuvp3g1jMYVE3r2Oe/A0urIW4eSQnTcRSqB2HZo7OwnpDKquCXXwFL73voVqPYrbtA957UOsS+/jHjqLeNgDgNpzEPuDt1EZTXB7s+kwq9yC6DV1cMv2EXzYRVHaRUZrSDW3kRszE2TJirn4G5D93Yh0CjmbQJfVo5r2Ix/cM6/RaJpYZMclnKrG/OcY8IdQrUaYUzXvQw70G9mobBpdUms+yzvXcNtOoXYeMgMxgBwdMucYTyKmp8wKp9rtyIGHuO2vYF/6EF1ZjRzsI72jCYFFaThM8c2rJFsPMFVSgfXeWySOHEdnMlid91GtB9BbapF9fYhHg7inX0GXlmF99AFiagosG3W4HV1ShvXOm2be+/NfgskJrI7rpvb9+DFTe9ugZReh77wN2Szu6bOQzSGvGkEJ5wd/eEVyg8npRKNRmpqaOHbsGG1tbQwPD/Mv/sW/IJVK8XM/93Mv7A6ypfCxEdybBIvFYosuRFiO4POb/kUygf+f/iPcQ0ehusbEYkqbjSFd91HNe7DOvQPBMPbb30SgCV74DqGSUkoediLrduAf6memrhn5sItMeTX2pfcQnbcp77lDsKcTd9ch1K6DqO2NyO47RuUjNot80Im7/xjugRNYNy8hb18ld+qLyBsfIUeHIJXNb+60bl5Gl1ajS0x3k8hmEUrg7jxANmYy67YQ+M6/g67didp1EKG1Sc71GpLqii3YPZ0Eex8SxGesrW1q14nqesSsSZoRm0H29aAzAtVgVGGfMgrXPiLdesKowXbexLp2HrftDOQUcvQRalsLqrEV6/3XUW1zrnomgxgcxLrfgTp6BqYnUTsPYt24AqXlkE6iS8pQWxuRfQ/N+1aKTEkZjraI5LJY1y/jvvQywUgxZdPTqO07CPR2kwgVkQgXIS9/hFNchqquQRcXQyKBvHYVt/0Eas8+rOtXkZ33Ea6DOnIMtXc/1rn3sO7ewTn7GspxYWSEcCBA4PYt3JNncNtfwn79m4ihR7ivfQFVVUPup392TdfspUuXeOONN/jggw+4evUqf/JP/skXcj3RctgwF11rTXZunNHbCV5dXc2OHTsW/ZAePXpELpejoaHhmeN4tW4hBIG/8Kex/vfvog4fRd68ht67DzE4gGrZCbFZqKmGbA4iQTMf7bpmrW86CT4fZNNm3NBnoVyXpBBYA71ky6so77qBW1qBTKcQ6SRu+1lQLrLvPrq8Gtl9B+3z4x45g3XjA8hl0bVNRuvrzhXU/mPIjsuofYcBjXXtIqq6FlXXiHXtAkIpUtubyEVLiXZdRx06hXXZaKupxj1mXjwcwLpitpOo5n3Irtu4bcexblxEA2rvISgKmprxlfdwgxFEJoN0HZL7jhG+dQm1dYfpBrt1gZm63ZQ8NO557uwPYL//dQTg7joIwYhJ4AmBe+IL2OdfB8BpfxnZ0wXBMGJiBHJZnOOfwfedN3CPncG6dA7n5GeQg4PoaKlRq2ltIxebJecLExkfNQIM9Q0QMJlmeeMqRKK4B48Ya6w1ibOv4bvZQWBijJnD7QSExD86jN5ej/XheVRzC6p5J9YFkyjUpWXoklJ0cbFJ0k1OkDl4GLsoinjYi66qQd7qQB08jC4uxrp0kfS33kbvbl319Xvx4kX++l//6/zv//2/qaurW/kJz4fvXhd9dnaWy5cv09TURENDw6oWEGqtcRwnv3fZ929/EXmzA/fVzxnNr5MvI27fQrXuw7p6CcrLsT58H1BY595FxOPYF94zcjxXPjSiANcvmrjto/dwY7OUXD5HODZDyfQkqd2Hmd6xB51J4YSLEHevY105j9rWgo6WGVf6wHHsD95GV27HPfYZZF831rULuIfOIO/fQuSypt4+N8ctR4cROYXacxAAFSyi+NZVdMNeI/MEqO0NyJ575tyy4LadQdU1I3qMy440E026YRfW7etYH11AJHO4B85A62HkXCLON2zaWGPBCNbVC0xvayVYVGqeKyXWvQ5U6zF0MIxIZ5DXL+C2nUS1ncJ+/3XcI2cAEOkcurQCOdCLatpj3vP1y8Ylv/oh7rGzyJ4edEUN1vVLqGOnyWSy5EIlFD3shkjE9KX7TAOIdfFDVPMuVMse7LfeQDfvxPnM5yh65y38sVlyn/8SkUeDBK9cZLaohOT4BKm9B3CrqrHf+KZZmnDsJSOh1d+H6O/HzeZIn34F23Wxzr+HDkeMCMQ+I59lXfiQzK/91zWR++rVq/zkT/4kv/d7v/dxkHtDsaEWfGBggJ6eHg4ePLhizXBsbIypqSl27dr1VLztrWSV77yF7+/+baMAMj1hfsZmUOXlyEQCVVKMTCXR0SgimUAXFZmf4TAinUD7g4hM0owvZhLQ102moprwcD/usZexLr2PLq+EZAIdLSG55wChy+8hsxkyNfUEh/txDxyFZAKr9x66qBhdXIGORJBdt9CNrWYT5+w4qmkf1sX3cY+/DLEZrDsdaCGJHztN0eVzJvnWehDZ2427/yDCkliXzhmBiaJSxNQEzrGXERLkg9ugNWJmCrd97jxLy2F22sydHz0DQqMtsK8ZS5dpOUCg6yaZskpkKk26cRd+WxK4eRkA56XPIR/cNc0qNdvR25uxLn0HgNyrX8b31tfQ0RJ0ZTW6rAZchXXlA9TWenRJGWJyGl0UNbJVR0+SSqQASbTjihG5mJpEV26BWAyr8x7umVcRExOI7i7U0XbTjXb5EurIMbRSWHduQS6He/YziMePsW53kDxwEHtwgMT2OmQ0SskH5vNxDh8lnk4TUS52ImEUXV86hQ4EsC5fQldVg99P9q//LdQf+5FVX7cdHR38hb/wF/id3/kddu7cuernrxEbZsE3jOC5XI5r166xd+/eRVcLLcTk5CQjIyPs2bPnGXKLe3cIfumz6NJSo5NWVGQGQEIRQ95g0Mj22LZZs+O66FAQMWv2d8tH/ajW/cj7t3COnca+fJ7YnjYi8RnU9h2Awrp7HbX3oBkvjZYYUYdAgNzB4/je+yZCayZ3tlHc102ueTe+UBj78nm0P4B76jXsd74JgLvvMGJ6EvmoDy0E2ZOfxX/+LQQat+04SAt596ppkb11zQx0tLUjxgbRFVVYHWb1lLunDetuB+7+wxCOIIZ6IRxBPuzCPXIK6+oHpgXU9iEScdz2M4jYNOmiKP7717CzadzDZ7CunENLi3jzPkjGCE2PInxBpOugGnch0mnk/ZvmRpGIIe/fRh0+iXX5fZxjryBGR7B67uG2n4bJCbD8yJEhyKRx65tJaolEUNRx1eRDJichXIwYGkRMT+EePYF89Ai0RvtsCEfN9xIuQoyPmiah+h2mQ3BgEDk2inPyjJk66+1Gba/DvniBTPUWZnY0EBgdoWhiDLZuRyQTqJadiEQc69JF1M5dqOotuF/+Q7g//hPLXG2L486dO/zYj/0Yv/3bv82ePXtW/fznwHcfwcHE3oViZmaGgYEBWltbnyI3Y6P4/+pPINJptGWZ8k7Ab1b3hoKmPhsJI5JJdLQIkYihoyVmsKSkBDkyjC4tMQ0vpWWIx4MIDXrrdlP3PvEy1sX3UbV16NrtRlG0pMxYSiHQW3eAP4guCiO77iFSSbIV1cSrtlLW2WHGTYPF6KYWrCvncfcfQww/gqAfp6Ia/9WLZPcdxtIO8u4thNY47aeQQwPI4QHcwyewrnyEDphVu/aFd9DFJYixUWOh209jXTqPqmswHoZQeDck98AxrI5LxvIXG8s/vfsQJYM9uAeOIhIzWHev4R45jXXFWPfk6c9j3btBYHKUmd2H8VkWod676IpqdHkNsrcTMTNJ7tXvx37/LfD5UTuazZUQKkJevYCu2YoOh0k6EsuShO/ewn3pZchkIecie7vBsnBb9yPv30dvr0M+6EU1tYDtR965ZcQWGlsQ2TSyt9fsB/f5UVu2GCveeQ/34GFEIo6uqSEXCuJ77z2EbZFqbCanNdofIBSPEezpMmuHgyHcP/SHcf/Un131tXr//n1+9Ed/lN/8zd9k//5n9gBuNL47Y/DVZByllMTjcRKJxBNyp9ME/uyfxHrzDUjEsb/zNjhZ7HfeBDT2O98GwH7rdXAd7De/CbFZ7G9/HTE4iP3e28iu+8iuTsT0LLFQMbm6nbivfgldXo3atgN5y1hMva0e68qHxtrMxHD3Hja91gMPkT33wB9B7TtqYtm6JsrvXMfZdYBkywHk+AjWxfPMHDuLdf0ScmQIlcmRmjRZbv/ta+iSaiid2xqSzSGmZ3APnUBMT5vXL4piv/0GBIpx97YbPW/bh+y8bX5fW4d14zKyvx8xOITbdnouRIii9h1GTE2gbB/Fo4NG/N91zU1txx607UdtqcPdd4TwuW/jn5nGeeUHiA73E751mXjlVmaCxVjXLqBsm9yZL2C/+zq6cSc6FDZZ+XAU6+I51M69qNJyMvE0tnYJdd7BPXnWeDyOi+ztRIfDuAeMxrluakHevWN00Swf8qMP0OUVqH1tWJ33kV3dqN2tZgJNCOTDB8iu+zhnX4VoFB0txpUWwbfehGgR+uBhgqWlRMoriMZjBHq7mWg/zjQw+MM/wvQP/tFVl7K6u7v50R/9Uf7Lf/kvnwS5NxQbasGz2WxBH7bWmlwux9DQEGNjYziOQ1VFBS1//6cJfONruC+fRd66YXZWPehBtbQgRx6bu3181gwduI7Z7uizEbPTJovaccVcSB+dI7ajiWh/LzoQhOISxNiIsTqxGXRZGda9m0b949gZrIsm3lNHTiImHkM6hRgfQ+SyOCdeQT7oRI4No7ZsN6KI8WkTN5fXkCyroPheB5M791HW/wDd1AKWjXX9CqqmFtXYjHXxvInD97aBP4js70bt3Iv10Tkj91S7HYTA3XsA+8O3zU7xbfXIR/35LLYuqzBVA62IHTgKuSyhLbVYA91mlDSdRCRiuEeN9VYVVegtdUZMIpNGdt02E15HTiMf9SMf3Ce15yApO0DZzcskWg8RGuhF7zmISCWxrl/CPXjMeBpd96G4nEBfN86pz2D19pgE2IMH6NIyo7Ty1rdRJ04h799Fte6HdBp56SPjxgdDWB+cM49t3YuYmUGMjoA/AFqhdu2GVBI5PISzpRb/B+fItu7DKikGaaFCIeTEOGJsDN3UhDpwkORP/FXGEIyPj5NIJCgrK6OqqoqysrJl20/7+vr443/8j/Nrv/Zry2353Gh8d7roKxF8sWQamPjd+Qc/i33uPXJ+P+FkHBmN4h8fBdtGjJqmD11SbATw9+1H3r6Je/IM1kfncNuN26u2bUflciQjUQL19diJODoSwf7gHVRltbF06RRu+2nkrWtGqG/kEfJBF+7RU1iXP0BbFu5nvoj1wdvmsXsOIocGUdu2gR3Aun4JXVWDs7cN33eMRzHZ1k7p/dvIdAonECSx5yAl1z8CwG1rB58P6+qHJs6+3WGI37IL68J3zMDGTbONxN13CNn/AOfYSeTQQ2TvfXR9E7KvB/e4uRGla7YSGDMijHpbPfJRH+7xlxGPB1BNuxCZFGQzRnr48vvoxl2IyXFUXSOquATr/m1EfAb38ClwHdOdduAY7tQEOQ2OFpT03CV15BQ+V6G77+NurSP40CxctM6/b2LoWBy9dTva9mO9/w6q/SVIZ8BxjPrp5Yu4Zz8DqTTWlUuonbvR5eUm+ZZOm8dFi43STjwG01M4FRVmeuxwO/bsNLqiEixpxj79PvTW7bh//E+hjp986rpSSjE1NZVP3IbDYaqqqqisrHxqbntwcJA/9sf+GP/hP/wHTpxY1Vbe9cZ3J8FzudyS+tJLkRvA/re/jP+n/w7ugTbk/bs45eXG6qTT5LbUEBh6hPvSS8jRUVRzi8mmR6PI2RmIxxCz08jHj0jtayN0pwPn6AnsKx+hS8oAbVofW3YZ99fJIcZGELkc7vEzyMsfoo6cgFQc624Hqr4J0f/QLLJv2W3iUsDde8gMVFw+j6raAo4iXV5OqOceqnkvwrIQj/pw97RhXzjHbNsRdDJOSbdpS81+5gv4Ln+IiMdwj5/GunAe1bILVVeP/f6bqO07EIP9xtIfOYF19SOco8cR0pqTYnLRXXfINrYQuXUV98hLxsXe3oAYfWR2nVVvQQ4+xG03yTZ3Txv4/RAImazz9YuAwm1/2UxtdVxCHTllNm9e+wi1ez+iv5/MjmbcsVGCfT3MNu4mKC18/gDW42FwcqgdTYhUBvx+5O0OVPNudGk51oUPUAcOwtgoVFSBkIiebvTOXZDLQc4Fvw957y6qrs4IP6RTxntCkxwfI1C7FV8mZdx120Jv3Y7asx/1ymtQ37DiNai1JpFIMDY2xvi40YZ//PgxkUiEv/f3/h7/+l//a15++eXnuMrXBZ8ugi9Hbut//S7+n/hx1N59ZitIOIIOBxGZDCpahO7tIVVSQsn1y6Tq6gmODJvZ4Koq5GA/7kunkJcuEDvUjsikCJWUmOGLm1dRLxkLr7bVI0aGEI6D88rnjEzv8IBxwzNpQ/SrF83mS9cxC+WFQDfsRFfXIK9+iN7aYNYGHz9llu5dvYiybdxXv4Dv298AMFs1pqaQfb1GH679NPLmVUR8luSWOuxUCrVtG4HRx2bMtHkXsrsTtXMPqqkJ671vG7mpqQlzAzp8zMgIHz6Off0iuYpq7EzCWL6KSkQ6hdqyBTk5jtpahxwfRtVuN+o0xaXgZI3IZKQY68YFdFGJWSKggZlJdFkVKNdkt6Olpoml+x5uzVYygwPIxp1YrsK600G8th5rdgZqtxMaGkT6A+hoKToURiiN6LyHrqlFV1YhO26gt9WZkmUmg5ieQm2vM1JOA/3oykozVBMMmJ54lWMylaK4rBTfjgYj83TgIOxpM81Kz4FsNst//s//mV/6pV8C4Mtf/jI//dM/vVE95oXi00PwhZ1pTyXisll8P/t3kZ1mH5i8dBF1+Ihx7055Vq4FMTiIjkTI1deTc12SkSLCM5NYoRDBW9fRto1TXol/ZBj31MtYH76Pu2cflJRALgN+H9blC8Y6DzxEKJUfWxQPO83PiXHcXfsQ42PoxkbTHHPxQwCcV78Pefs6cmwEd98hUrMxgtPjWD4fpHOofW3Ii++jWg8hB/uNO19caqamarfhtu7D9/YbACQOtePG4/jiM1BRSejOTVRlFWI2BsEg7snTiLERs7yw5z5oRap2O+HhQdwTJsPunjBxuXv0JeT1j3APn0DOTpje8WzaVB3A1PhLSrHu3TAlsslRkz2Pzq1wmhgzDTlKIft70dXbcLM5YpksRVphaw2BIDochbnyZaa6FiceJxuMEFAuoYkJaGox46bTk6YCkctBMGSWRGTSUF4B8ThEo6Z8JgT4bZCSnHSZTswSOX4a3//7R2BH8/Ncgs9gbGyMH/qhH+Lnf/7nOXv2LO+++y4vv/zyJz3b/d1JcMdxnupO8zrTgGUTH3m4LqK/D9Hbg3zQi7h/DzE9ifXu26g9rVgXzuPub0Pevgk+P9mKclTOIdbSQjSXxZdOmTHFTBr36HEjhdSy21jLPXvR4TD2e2+iareZslQuh3v6VdOeevsaensD8u4tdDCE23YU2deFSJkVuCBw6nfgDj0iODaC2rod1dSC/f675r1/7vtNN10qiWpsQW/ZivXhe+iSUhA2avcec3OZmUXEZnFa95P1+WBqgmxJKaV3OnDLK5CpFCKVxD1+Eh71M1W7jah2sUrLjDvr8yFScbAkYmzYLAaYnURX1RjV2KotcyOnytTTe++hjp5G3u9A7dqPDgQRqQSk54QiA0G040IkgkpnSI+P4Y+WYBeXQjqNQJrEZOUWQBgVVEAryNp+k3eZnUUGQ1jFJQQ9TQCfDbbP6LFbXpihDbGFRqgcjuswZvkJ/tW/S6T1wPNceoticnKSP/JH/gg/+7M/yw/8wA+s+/E9uK5Le3s727Zt42tf+1ohT/nuJ3h+WGSh1V4jxKNB5LUriHt3sN56A8fnI3DhPLmt27HHxxDZDDO7WvEl46iGJkKPB7Ae9uLuP4h16wburj1Y9++ZuLe5BevdNyBk7uJidsYQSkjklQ9Rx05jffg+urgE56Uz+L79dQAm9x4kHAwTuPoh6mC7yRjv2YsYHDBTWFXVc65nGHnvDu6hIxApwjr/nvl8zn7GhA9XL6Gbdpphmbod6KIiMuEIaQHl1y6S3l5PYHgQLSRqyxbs4UHctkNYt67jHjmGvH0d9+QZs4o3FAA3i5ybPhPxKbRlQyCAGH2EajuOeNyPam41+Yd02sTjNy4YCeT6FjO3PviAyR17KCotwd91B13fjHYVwp5LUk2Oo6u2gpQmKYZA+QMIBAhJLp0iIyS5VBIpLfyhIP5gEEsCaHNHyGUhFERMj+OUlDNQXU/0b/5jomuU2l4O09PT/NAP/RB/5+/8HX7wB39w3Y8/H1/5yle4fPkys7Ozn36CO46zZLy9XpicnKTn+nUOxqaJ3LuN9Vu/gS6vwOq4jgYyTc0Ee3uYPHSUgCUJDT6Amlrk7ZuoxmbEwwcm8XakHes7b5ohkvpG5IMesy4HsDqumkX1/f1k9x8iOzVG0fBjRDaDe/yUWS38oAcNuJ/7Etbbb5hZ8dOvQDqNdfkj3EPtyJ4uVOs+I+543WTLndc+j4jHjKpJ6z6sG9fMefU/RFdVk9zRQCYWw4lGCeUyWGVlBDIps/Sw6y6qvt4065w8i06nTCiiXUQihtmv5JhSW81Ws0G1uNRUEB49QNc3I6+eRzftwd3WgHX/JuLxION7j1Ec8OPvuYtb3wzREmRfLyiFatxtViXPTJkEo8+PGB8z+RLLBoFpE46WggA3mSBj+VDTkziBACE3iy0Flt+HnBrHqazh/vHPU/Nn/sqGbAeZnZ3lj/7RP8pP/uRP8sM//MPrfvz5GBwc5Ed/9Ef5mZ/5Gb7yla98+gmezWY3lNxDQ0M8evSItra2p0ZQRccN7D/4HcT9u9ivfwMtJap+B9bDB8zu2oNG4JdgB4P4Oq6harYgpmcgEsE5dhz7O28islncYyeRFy+gTpw0ywuuXgIg134SmYwbZc/2k4j+PnR5GZSUYn1wHtW6D1VUhHXtqnH9z5xFDD1C9vaYfVotO03tfnzM1IHjcZxTLxsl1FQChMDquEbmcDuB65dxt9UhJ8f4/7d35vFVlee+/6619s7OnJARCFMYEkISMhHGgMwqAokoVWuv1uK1tsd7tB5Lj4e2p55Pj1arvXr1WP14Pw5XW70FwlgGAVEIcwJJCAkQQhKSkGRvMuzMe1jrvX+8mKtWJUAGCPv7FwT2u96V/fmt913P+zy/x/D3R+hOXBZv9OAhKOHhWHQnSoC/fO9PmQLt7bKXttMBPj4oXR0Y/gFyVddU1PMlctufloFysRIRNQZhsWA6uBthMtMwfjKBJhXz6QKMpOngdEJnG4rTKYN1mhml/CxG7GRZM3/JKktIdQNMJhlfabOjdHZijJkgHyqqChZvsNbi8PbBVHkWdAMjKIjTyx5h+H97/Kr7a/eEtrY2Vq5cyU9/+lN++MMf9vr43+Tee+/l2WefpbW1lZdffnnABd6nvVaeeOIJIiMjycrKIiYmplfHFkJQVlZGe3s7qamp3T2kuv99chKuybKCy33kEOqenZj/91sIVcXf6UQtP09XbBzuehvtExMweZnxr6tDAKa8PERgGPqkOLT9n6MIAU434tw5WmLj8Q3wx3z0EEJVcS28E9Pe3SguF4a3N0ZwGEJVUSorUIePwEhNRz16COwtKNZL0hvM7UI7dgShqnJr7XahXKxBO3USpbkZfWYG6rHDdMUl0tXRiTp1JoqvD0bHCITFjGqtw6IIjK4OnM0NNPsH4FVfiyk+GZ/8oxhxiWiH9mAkpqHl7UdPz8C8bzvG0BGybFYz4Z5zJ0p7iyyndTkx5e7HPWwkjWEjCKsoRm214546D8XtRq0sxRg1HlxuGflusEJQCEpHmyzeCQxGLT+LER2L4e+PqeCYbLw4SfZCF77+gJCW1UNH4Jt/CPwCcIdFcnreCrpmLaaoqIjAwEDCw8MJDQ39h+/zWujo6OD+++/nkUce6Rdxb926lYiICNLS0vj888/7/Ho9oU9XcJvN1u0dbbVaueOOO8jKyiIuLu66VvMvzSMsFgsTJkzo+Vjt7Wjbt2D+8+sohfmImImyR9bY8YjWVjqGDceFIPRkPgD6pEQUpwPDzw9RU435kk1mX8XEoVSeR+nskGfKAQEIkwlFNy6bBsZhDB2O6fPLpoGL7kAtPYtafl4KuuI8Ykw0QtMw5eyTZagpaaAqiIBAlLqLKBXltI8cSUBZKfqsDNTDB9Bvmy8r6SIjUdraUC5dlFl8FjPu0FD0BhtdZhP+VedxjpmA35l8mcZaV4UIGyrHdnYhVA0t/+DlXPfZKC1NuMzetDudhJQWIrx90afPRz11AtVag546CwClphLh6yeDe8NGQXOjbKzY2owxIQGhmdAOfQZBQzDGyLxzAaj1NajV5eipM1FL8hEjo3GFhFM1bCxBv/oDfn5+CCGw2+3ddkkWi4Xw8HDCw8Ov3AP+W+js7OSBBx7g3nvv5bHHHrvqz18Lzz77LB9++CEmk4muri5aWlpYsWIFH3300ZU+enNu0b9KU1MTmzdvJjs7mwsXLrBo0SLuvvtuEhMTexZRv4zT6aSwsJDIyMjrqtVVjudifut1tC0bpVNIQT5GeASYzDiiRuB0dhF0sgCA5oRkTIaOb3ODzFk/ehjh74979m2Ydm6TySgz54AhUPOOYKRMQc3LxUifBl2dqCeOy9TT+QtRz5VKoaekoZ4vw5gUjxHgj+mz3eDnL80MqqtoTk4hqPI8xuRk1JIijHHjUU8cxZg1B1qaEBGRKI5O1HOnZKJJoxV9XAxqeyu6vz96qx2Hy4XicuDbWI8YPR6voly5ire3IoaEokfHop0uQLVdpHnCZAJtFxHj4hAomHL3Sc+2mQtRy06j1lRgjBqHCB2K0lCP0tqM8PGXpasmE0pLi3SOHREtjx8bbSj2Rlkl1mhDnz4PxeWErk7cXt60d3Sgv7sdv+/Yln81OUUIQVhYGOHh4fj5+V3xge5wOHjwwQe56667+PnPfz4gLiyff/75DbFF7zeBf5Uvo4vZ2dmUlpayYMECsrKySE1N/V6xd3R0UFhYyLhx43ovMaGmGvM7f8b0f97DmDgJ7cghhLc3YtQY9NBQnNUX8KmtRXW7aY2Lx+wfgKWoAGNSPFpurlytIyPR9n0hK8WmzZBOrXnHMIZHgbcPRkQEatUFmZLZ2Ih+21yUjk7U3KMysJd7BBEYhJ6SSmdLC25vLwKEkEkdqgIWM0pnB8LPD+3g53IX0NUuq+dcLlmPHhouV+rgEFnv3tGOiJ6AWlZCV+IUtAtltIcPh6BgAspPo/j4orQ0Q2c7jXFTCHJ3YTpbiB49EbXJhjFuEjRZUetqoevy+3pzA9qZk/LYMGk6StMl1LISjPGTUBqsGGNiZJHL8QMIX3+MkWPlQ0uAVn4axd6IMyENUVVO61+/wHfU2B59RU6ns9vmuLOzk5CQEMLDwwkODv4H8TqdTh5++GHmzp3LU089NWAWS7e0wL9KR0cH27ZtY/369RQVFTF37lwyMzOZNm3a197DmpubKSkpISEhoU+CMbS2YvrwfUxv/i9E9Di0A9I2qTN1Cl4WC2pxEYaPL1pdLe0jRuIKjyDoRJ7Mj1Y1xLjxKBerZQfO6mr0KekycHUgB+HjIyPn3t4oF2vA5ZL/J20KeJkBBcPPF23/PrpGjsSntgYRHCxr/cxeGMOGIvx8UXUXws8Xbf9ujClToaUFERUlj7sUwOlEeFvkNjw4FKXNjhAG6sULqPXVuNMyMOXlYJi8aBo/CbOzC5/aKtSAIDTrRfTEdISXGdPxgwiLN8aYWPDzRz2ZixE9UYo5aSpCVWXDCMCdcTuKvRG1KBcjZSbqiUMYk5Kl28uBT0FRZSPCRhuO2ES6Gi6hPPYM5qX3XdPXpOs6jY2N2Gw27HZ793t7SEgIQghWrVpFeno6q1evvpn80wavwL9KV1dXdzO3EydOkJGRQVZWFqWlpYSGhnLHHXf0fR8olwv3Xz9C+9NLmMIjsORddkGZtwClowM19whGcipaXi5dY6JpCw0jLO+YPB5Ln4bi5YVakI8xcSLasWPoqWmI0FBMuz6VO4PxExBeXgh/f9S6WtTSs7hnZmA6mIMzKgqTqoDFAt4WabDg7wtuFyIsBBEShuJ2IHx90A7sxUibipZ3CD19hixemSYdX/TYBNRzxbLkNHoCatlp3FPnoAiB0lCPCItEy5UPsJbJ0zE62gg8X4wjZjI+JfnoCWmIwCGYDu6W9z5zobS9KjyKnjIT9dQJjMQ0hKJgOibHcU+fLx1Zi3IxEqfIHP3QCPSkqWjFJ1Dra2gZPQHLhEnor/21V76qr763r1mzhnPnzhEbG8vbb7/N0KFDe+Ua/cStIfCv4nQ62b17N7/97W+x2WzMmzePe+65hzlz5vTIIeZasdvtFBcXkxAfT9CBHMz/84/Q1o56ukSmtM5bILts5B7FmD4TLScHV2IincFDCMyRq1rbtBlYVBXT2TMYsRPRDuRgjInGiJuIWlIiTQy8vVGrq3HNnUdXXS3qkBAsvt4oLjfCxyINLVRQWltQ2uwgBLg6MaLHgaYiAgJQG22IISEygh0xDNV6EREQIPtzqxpGQBBaWYkUddAQ1Lpq9HETcTmdOH398fX1xZR/GABH+hx0ezO+ZwtpGTORgMpzuJLS0VxO1OJ8FF3HPSUDpa1VprqGD5X5KtETUOprUBobUFrt6DHxiCFhaMX5GEOjZC09YEuYSrBZwf3qXyDsu73vrwVd13niiScICAggOjqa7du3s2nTpoFOP70abj2BA2zYsIFdu3bxyiuvcODAAdatW8f+/ftJS0sjKyuLefPmXVOE9buw2WycP3+eyZMnf61HmnL4EOZX/4RadBKlqUmeWc/KAAHasSMYk5NQ8/IQCQm4IiLx2rsHRQiakpKxGAKzqqD6+6MdyEEEBMi2QmYv3P5+OC5U4u10oAUHo54rxUhLRSk/j5EwWdpTBQaidLRDWxNKc5Ncwf38UFqbZETb7YTODsSw4ahV5ejps+TqOXqcLKBxOXHPuBzgamnC3daCpa4aIzgUgkPlu3v1eZT2DinQyVMxdBfmUyfQvSx0hg6DISH4nSlEjJ4gU12Tp4GmoeVebouUPB1MZtTiy+m9Z4oQgcHok6egnDmJ3tmOSTPh/sVzuO+/sj/51WAYBk899RRDhgzhxRdfvKqAbU+pqqrioYceoq6uDlVVeeyxx3jyySd78xK3psC/nNtX36V0XWf//v2sX7+evXv3kpiYSFZWFgsXLryuTo/V1dXU1dWRlJT0nTsE5dw5TP/1OtqObaAbqHV16EnJiOBgtIICjOHDUcvOISIiMOITUAvyUWtraU5Kxr+0lK7YGEyhoXiVn8fw8UU5d1be24TxqMWn0GfNQjt0AD0jA+1gDvrMDJQL5zESEmVmmLcFxd4MKqgXymXiTtEJ9GkZqCfz0BNTUVuapbWwn/RvE77+qGelVVRrfBqmpga8IiIBgVZwVP5Op94GnR1oJ4+hJ01DLTiGkToDoZkwHZFmjI3JM/Cy1eJ/sQLH5Kl4FRdgJKfLiPvRy9v0abehGDpqUZ7cph/dj6GZcGUsREPgfGsd9KIADcNg9erVmEwmXn311T4RN0BtbS21tbWkpqbS2tpKWloaGzduZNKkSb11iVtT4FdC13UOHz7M+vXr2b17NzExMdx9990sWrSoxymPQgjOnz9PW1sbCQkJPUuwsNsxffgB2icfo9hbUCsqMKKjMYZHyfRQ3UAtLQWHA33xYpT2dqivw+Xri/fJQtpGjsTcYkcJCUUdPVqWeYaEoDY3IYYEo9ZUI4YORT34BcbU6WiH96PPlKWu7oy5qGdPYcROkvZIEZGo50rk+3ZYBOrFC9IM4lgOwsuCMT4O/PzpaLXjf7mJgz45HcVWjxgqu6RoJ6QZhXv2ItTzpag1FegT4lEvWWXfL2udNNlwOuhKm4VyoQxvWy0u/yDwD0SJGolaUYricMpdwPhJuIOGYCo8ghgZjVZRhuPDHRjps675u/4mhmHw61//mq6uLt58880+E/e3kZmZyRNPPMGiRYt6a0iPwK+EYRgcP36ctWvXsmPHDqKjo1m+fDlLliz5Woukb37m9OnTqKpKbGzs1UddhUA9kIPpL39BPXEc9dQpjLHy6EcMG4bh7YPpQA44HBjTpqGWFONInEx7exveJg2Htzfm2osY4eH4NjdC5FBURydGQCCq4Zb10y1NEBCIcuaUFPPpQsSwKNn8UNdl65/qCvRZ8+RRVfhQcDnk6j02Bu3ofoSi4ogagzkgQB65VV1AbbAigkMwosbI/19bBV4+YG+SZaf1tajlpTJ4OHMBakWpfHhMTkc9fxZjYgJupxNLgUzdbUqaik9XJ95nT+KMT8Xr5HH0iGEYkyaj+PnhfO3D6/h2v/lrFzz33HNcunSJd955p1ey3npKRUUFc+bM6c686yU8Ar8aDMPg5MmTrF27lu3btzN06FCWL1/O0qVLu9sT67pOYWEhwcHB39uQocd0daF9+inqvi8wrV2LMXoUalERmM3oszJQOjtxOrpwNzfhX1GOMW263I5PmYJacIKuESMwX7KiBwdj6mxH8fNHMZzS9DAoENw6InSIbLAYEXHZNmqkfN8eORq14pxcnScloRUXYIwZjzCZabN4o1q8CPjSMiptlkxSURRZ2XZcmlno85ag5R1CsTeip89GuVgls+T8A6V1lcUb96x5mPbtQnG70WMTUJwuhFmTzRXPFKMIwaXEKXhbL+JvvYhr8hRMJYV07chDjO7ZmfeVEELwwgsvUFlZyfvvv9+v4m5ra+O2225jzZo1rFixojeH9gj8WhFCUFJSwrp169i6dSvBwcEsWLCAnTt38sYbb/xDq6RewTBQ8k+g7d2LUl2DaUM2uqriVgTeNhv6nDmopWfQJ8ahNFyC4GDobAMUdLcL0dSAy9uCX/UFnLETsZwtRk+fjpZ3+P/7zsVMRD1fKrt5xMjuHUbUCFRbPWrZGSn+0hKcYRGYXS6MCXEIVUHLPSgz7+JTpKd6XZXsYnLmFAQNQU9IlqaQTgciJAw9JgGt8JjcLQSHyXpxb4vsjFJxTubTz5qPVnQc7E10RETha71IR2IqNDXQFJ+K/ZnnepyF9n0IIXjllVcoLi7mo48+wmTq01KKr+FyuVi6dCm33347Tz/9dG8P7xF4byCE4LPPPuOhhx4iOjoaLy8vli9fzvLly4mMjOyzxIjq6moaTxWRqKiYrPXSN7ypGaWzHaWxAREaimKtwxgzRrrFRo9FuWSjKzwCh7WWVpMJb4cDHx9vfAw3WlMDQtPQ8g5jJKagHdkvXVaFjmJvxj01g47GBswWL7y8zGgFuQiTGRE1WnaEqauRri0XqxEBQdKb/cBuMAxEbCLC7IVSV40YNRYt9xDG0OHok5Iwf3a5uUN6BiBQiwsw4iajHT2IERDIpdgEwk8ckg+Q9AyUijJat+RgcxvdWWihoaGEh4cTFBR0Vb9vIQSvv/46x44d45NPPunTo9Jvu/bDDz9MSEgIr776al9cwiPw3uKll15i4cKFpKSkUFFRwfr169mwYQOaprFs2TKysrIYPnx4r4j9y3bJdrudxMTE69pOtrW1YbVauXTpEiaTiYiICCIiIrA0XkIpLkQtO4O2YyPKxSoMezOmrk5p3VRVgTFugixsOSQj4u7pc2W0+8QRjJRpaIf3Y4yLQR8xGvM+6QzrnjwFxWxGyzuEMXwUSn09RlKabBxht0ujy9FjEZHD0Y7m4LZ4o/r4w5AQRKAfakU5rv/+JO4nVnffg67rNDQ0YLPZaGlpITg4uDsL7fuCZEII3n77bT7//HPWrVv3NWfU/iAnJ4fZs2d/rW7i+eefZ8mSJb11CY/A+xIhBDU1Nd1idzqdLFu2jMzMzO/shtqTMUtLS3E6nUyaNKlXo7ydnZ1YrVasViuKohAREUF4eDiapnHqi73EnilkyBc7UMvLpM9a1CiZxjp2PEqDFaWyXNaoJ6ZerjvPQwQGI0xeiPGxqIXHEOHDUSvPoydPkb7ueTIhxj1rPkpLkzyem5qBdjgH+4Q4fCIi8DogHyD69NkIX1+cf/4IfHy/9R4Mw6C5uRmbzUZjYyP+/v7d1sZf3XoLIXj33XfZtm0bGzZs6PtMxoHBI/D+QghBfX092dnZZGdn09LSwl133UVWVhbjx4/vkdgNw6CkpASTyURMTEyf5kQ7HA6sVit1dXW0trYSERFBdHQ0fn5+qMcOYn7rT1B3Ee1UgWw/nD4LnJ1ohcfR45PRTuajp8+UppKHZStj99zFqGWlMnEmOV32CBs7AdrsskuMELjnLEA5eRytuQkjLAJ0EOMnoFSWodjtuP79Rdw/WtWjexBCfG2H4uXlhb+/P97e3uzdu5f169ezefPm68pzuMHxCHygsNlsbNy4kezsbGw2G3feeSeZmZnfWdOu6zonT54kKCiod6LzPaCrq4v8/HzGjh2L2+3GarXicDgIDQ0lMjKSwLIzeL34G2nvfPQgAO75d6IW5KI22NBj41GamuT7eYMVxW6Xbiup6SjlZah1F+Vx2W2L0I4fQWltwZGQAhcqUUePRvXyRjt22XF2/u3gcuD8aJNsBnkNdHR0cOzYMZ555hnq6+t5+umneeCBB4iOju6tX9mNhkfgNwJNTU3dBhZVVVUsXryYu+++m4SEBFRVxeVyddeqjxgxol/m9GUJ7cSJEwkODu7+udvtpqGhgfr6etrb2wkNDWXkqXyGvPw7FHsTmCzyfH5SAorNilp+Th6Xzb8Dbf9n0mtu2mxobUWxXkSMHod29BAiIpLOmIn4Xt6OG2MnIIaEoZaWyE4lFypx/td76Mvvua77ys7O5p133uG9995j7969eHt78+CDD17XmDcwHoHfaNjt9u6a9nPnzjFr1iwOHjzIe++9R2xsbL/Mob29ncLCQuLj47836eLLEkur1UrbJRsT9u4g6v0/y/PsaRnQ2YlaU4kRNQotX5pLGmFhaIX5KE4nxvAojNFjMR3aj1BVOiKG4zVyJKZT8rxdKypEhIahJ6ehNFhxbN8P17Fz2bJlC6+//nr3sWZfsGPHDp588kl0XefRRx/lX//1X/vkOj3EI/AbmeLiYpYuXUpsbCzV1dXMmzePzMxMpk6d2meJGK2trRQVFZGYmHhVTqRfBrdaD3xB5Bt/JKj0DFpXF8bQYTLz7eDlQFlyOli80I5cbqxwMAdHWjrtLichl51u9JmzobER7fQp9LhEtFMn6dywHTFrzjXf144dO/jjH//I3//+d0JCQq55nO9D13ViYmLYtWsXI0aMID09nY8//rg3c8uvlpuzffCtQlFREZ988gnbt2/n6NGjLFiwgPfff58ZM2bwL//yL+zbt6+74UNvYLfbOXXqFElJSVdtM6yqKiEhIYxedjeWTZ/R+cNHEKqKfUgE2v4v6ExOx502HS3vGNrBA7hn3IZSWio/W1JMUJcTPXaS7NVeWYVaUoI+LQMF0OcuuC5x79mzhz/84Q9s2bKlz8QNcPToUcaPH8/YsWPx8vLi/vvvZ9OmTX12vYHEI/Be4Ac/+AFTp04FwMfHh8zMTD788EPy8vJYtmwZf/vb35gxYwb//M//zGeffYbL5brmazU1NVFSUkJSUhK+vt9+BNVTFG9v1N+/jONv2wi63KHEUlRIe2MzHcOiLl+wGeFy0zJuAkpCMtrpEtQzZ3AvWIJSUyVtnp0uaOvA+dvfX/Nc9u3bx3PPPceWLVsICwu7rvu6EjU1NV/z8xsxYgQ1NTV9es2B4oYT+Msvv4yiKN2dIG9mLBYLS5Ys4d133yU/P5/77ruPrVu3MnPmTH72s5+xc+dOHA5Hj8draGjg7NmzpKSk9OqRkTFzNl2f7sd9x1JESjpBZ0/j09BAy6y5mE4WoDU2YGnvxBAKwmwGXz9pRTU2Rmbd1dRgpKUjEiZf0/UPHDjAv/3bv7F582YiI3vXDOLb+LbX0pvI3umquKEEXlVVxa5duxg1atRAT6XXMZvNLFq0iLfeeouCggIeeeQR9uzZw+zZs3n00UfZunUrnZ2d3/l5q9VKWVkZKSkpvWpy0c2QEJwf/F/05SsQl8f3PX0Ge8wkDP8AjIAgvA7m0BE+jNaJCbJr6ukSROQIRPQ4XGv+/Zoue+TIEVavXs2mTZsYPnx4b97RdzJixAiqqqq6/15dXd1v1+5vbqgg27333stvfvMbMjMzyc3N7fOt2o2AruscOnSou6Z94sSJZGVlsXjx4m7LodraWqqrq0lOTu6XHGylMB/jhd/jt/Ny7vms2WCzoZ09jT5qNNgaaBs3Dp+K86gCHA8+BC+8eNWr4PHjx/mnf/onNm3a1DdFP9+B2+0mJiaGPXv2EBUVRXp6On/961+Jj4/vtzl8g5uzs8nVsHnzZqKiokhKShroqfQrmqaRkZFBRkYGhmGQl5fH2rVreemllxg7dizh4eHdVVT9VT1lHRbFhZ8+wTQfH8zbtqKUV6DU1aHPyID2DrTySoIKCnDefifG2dOcWpZJ25EjhIaGEhERQWBg4BXFXlBQwM9//nOys7P7VdwAJpOJN954g9tvvx1d1/nJT34ykOLuU/p1BV+4cCF1dXX/8PP//M//5Pnnn+fTTz/tzgC7VVbw78IwDNasWcP69esJDAzsrmm/6667umva+wKbzUZ5eTkpKSmYzWbUv32M5X88gdLVhZ6cArqQ1Wi+fijVNbheeBH3T3/aXUhitVppbW1lyJAhREREEBwc/A95+KdOnWLVqlWsXbu233IGbnAG9zn4yZMnWbBgQXdU+Mt3oqNHj95s9re9hsPh4Ne//jW///3v8fLyori4uLumPSQkhMzMTJYuXdqrD0Gr1UpFRUW3uL9EOVWE5bFV0NyKWlmJCAmRRhWVlXQdOfoPKamGYdDU1ITVaqW5uZnAwEAiIiIICAigoqKCH//4x3z88ceDdtW8Bga3wL+JZwX/br6sUlu3bh2bN2/G19eXzMxMli1bdl017VarlcrKyu9+z+/owOuXz2D64AP05BS0Eyfo2rwZY8GCK87Xbrdz5swZVq1ahdPp5KmnnuLxxx/vk1bBNykegV8Lv/zlL9myZQteXl6MGzeO9957r89SHwcCIQTl5eWsX7+ejRs3omkay5cvJysri2HDhvVY7FcU91fQ1q/H9PLLEBmJY+PGHs+1oqKC+++/n9WrV3P69GkaGxt58803e/z5Qc6tJfDe4tNPP2X+/PmYTCZ+9atfAfDiiy8O8Kz6BiEE1dXV3TXtbrebZcuWsXz58u+taa+vr+fChQtXFaFXqqvB5UL0sLqrqqqK++67j7fffptp06b1+J6ulZvwwe4R+PWyYcMG1q1bx1/+8peBnkqfI4Sgrq6uu6a9ra2Nu+66i8zMzK/VtH8p7pSUlD6L0F+8eJGVK1fy+uuvk5GR0SfX+CY34YPdI/DrZdmyZdx333386Ec/Guip9Ds2m40NGzaQnZ1NQ0MDd955JyaTCbfbzTPPPNNn4q6rq2PlypW88sorzJ07t0+ucSVukge7R+DfxfcdvWVmZnb/OTc3l+zs7EGbkthTGhsbWbNmDZs2bSIqKooFCxZw9913Ex8f36u2UjabjXvuuYcXXnihNxsEXDU3yYN98Ce6XCu7d+/+3n//4IMP2Lp1K3v27LnlxQ0ysaauro4zZ85gGAZbt27lxRdfpKysjEWLFpGVlUVycvJ1ib2hoYGVK1fyH//xH30m7p4+2E0m02A2irgiN/0K/n3s2LGDp59+mi+++ILw8PCBns4NTVtbG9u2bWPdunWUlJQwf/58MjMzSU9Pv6qa9ubmZlasWMGzzz7bLbSB4IMPPuCtt95iz54911111w94tujXwvjx47u9yQCmT5/OW2+9dU1j3WAOIH1KZ2cnO3fuZN26deTn5zNnzhwyMzOZMWPG976vt7S0cM899/CLX/yCe++9tx9n/HVuwge7R+ADyQ3oANJvOBwOdu/ezbp16zh27BgzZswgKyuLjIyMrx2rtbW1sXLlSh5//HEeeOCBAZxx7z7Y+wmPwAeSQ4cO8bvf/Y6dO3cC8MILLwDw7LPPDuS0+h2Xy9VtY5yTk0N6ejpZWVlMmzaNBx98kIcffpiHH354oKd5M+KxbBpIbiUHkO/DbDazePFi3n77bQoKCvjxj3/Mrl27SExMZPHixR5x34Dc9FH0/uBWcgDpKSaTiblz5zJ37lz+9Kc/9Xs7IQ89wyPwHnArOYBcC33iMOOhV/Bs0XtAeno6paWllJeX43Q6+eSTT1i+fPlAT+uWZDB59vUHHoH3gK86gMTFxfGDH/ygV2qZq6qqmDdvHnFxccTHx/Paa6/1wmwHL4PZs6+v8ETRB5Da2lpqa2tJTU2ltbWVtLQ0Nm7ceEscv10Lg9izzxNFH4wMGzaM1NRUAAICAoiLi7slo/M94Vb17LtePEG2G4SKigpOnDjRL/XSNyo98ezzcHV4tug3AG1tbdx2222sWbOGFStWDPR0bjhuAc8+TybbYMXlcrF06VJuv/12nn766YGezk3BIPTs87yDD0aEEKxatYq4uDiPuD30CZ4VfADJyclh9uzZJCYmdtdfP//88yxZsuS6x9Z1nSlTphAVFcXWrVuvezwPfYrH8GEwkpGR8a1psL3Ba6+9RlxcHC0tLX0yvoebA88WfRBSXV3N3//+dx599NGBnoqHAeZKW3QPNyGKoqwDXgACgGeEEEsHeEoeBgjPCj7IUBRlKWAVQuQN9Fw8DDwegQ8+ZgHLFUWpAD4B5iuK8tHATsnDQOHZog9iFEWZi2eLfkvjWcE9eBjEeFZwDx4GMZ4V3IOHQYxH4B48DGI8AvfgYRDjEbgHD4MYj8A9eBjEeATuwcMgxiNwDx4GMR6Be/AwiPl//hiXjnf7mBMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure()\n",
    "ax = plt.subplot(projection='3d')\n",
    "ax.plot_surface(res[2]['x_kernel'][0],res[2]['x_kernel'][1],np.exp(res[2]['y_kernel']),cmap='autumn')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d05d6915",
   "metadata": {},
   "source": [
    "# Fit via docker <a name=\"docker\"></a>\n",
    "\n",
    "\n",
    "Start docker. Make sure that you have the PGAM docker image by listing all available images with the terminal command,\n",
    "```sh\n",
    "    docker images\n",
    "```\n",
    "It should list a repository named edoardobalzani87/pgam, \n",
    "\n",
    "```\n",
    "REPOSITORY              TAG       IMAGE ID       CREATED          SIZE\n",
    "edoardobalzani87/pgam   1.0       052c7daf6ff7   19 minutes ago   4.64GB\n",
    "\n",
    "```\n",
    "\n",
    "### Running the container and mounting volumes <a name=\"mount-volumes\"></a>\n",
    "\n",
    "If you are running the PGAM image in a docker container for fitting the model you must specify which of the host folders (e.g. the folders of the operating system running the container) will be mounted as volumes by the container. This is a necessary step if the container needs to read from or write into the host file system. The syntax of the command is the following,\n",
    "\n",
    "```sh\n",
    "    docker run -v <path to local folder 1>:<path to image folder 1> \\\n",
    "               -v <path to local folder 2>:<path to image folder 2> \\\n",
    "               ...\n",
    "               -ti <docker image to be run> <command> -c \"<additional commands>\"\n",
    "```\n",
    "\n",
    "More specifically, mounting a folder means that: 1) when the container is started, the content of *\\<local folder i>* in the file system of the host will be copied in *\\<image folder i>* of the container temporary file system, 2) whenever the container writes in *\\<image folder i>*  (saving/deleting files or creating/deleting directories), the same will be performed on the host  *\\<local folder i>*.  \n",
    "\n",
    "The docker image has a few folders that one may use to mount local folders to:\n",
    "\n",
    "- **/input**: mount this container folder to the host folder containing the input data and the fit list YAML file \n",
    "- **/output**: mount this container folder to the host folder that will contain the output to be saved\n",
    "- **/config**: mount this container folder to the host folder containing the config YAML file\n",
    "- **/scripts**: mount this container folder to the host folder containing the fitting script (fit_from_config.py)\n",
    "\n",
    "### Setting paths <a name=\"container-path\"></a>\n",
    "\n",
    "Make sure to edit the fit list YAML file replacing host folders with container folders, for example:\n",
    "\n",
    "```python\n",
    "fit_dict = {\n",
    "    'experiment_ID': ['exp_1'],\n",
    "    'session_ID': ['session_A'],\n",
    "    'neuron_num': [0],\n",
    "    'path_to_input': ['/input/example_data.npz'],         # '/input/', '/config/', and '/output/' are the path \n",
    "    'path_to_config': ['/config/config_example_data.yml'] # to the folder in the temp file system \n",
    "    'path_to_output': ['/output/']                        # of the container\n",
    "}                                                        \n",
    "# if saving the yaml fit list outsider the docker container\n",
    "with open('<path to local config folder>/fit_list_example_data.yml', 'w') as outfile:\n",
    "    yaml.dump(fit_dict, outfile, default_flow_style=False)\n",
    "# if saving the yaml fit list inside the docker cointainer\n",
    "#with open('/config/fit_list_example_data.yml', 'w') as outfile:\n",
    "#    yaml.dump(fit_dict, outfile, default_flow_style=False)\n",
    "```\n",
    "Similarly edit the fit script that will load the fit list YAML, pointing the \"config/\" folder of the container. In \"fit_from_config.py\" we would need to set appropriately the variable \"path_fit_list\",\n",
    "\n",
    "```python\n",
    "...\n",
    "#################################################\n",
    "# User defined input\n",
    "#################################################\n",
    "\n",
    "# frac of the trials used for fit eval\n",
    "frac_eval = 0.2 \n",
    "\n",
    "# PATH to fit list\n",
    "path_fit_list = '/input/fit_list_example_data.yml'\n",
    "...\n",
    "```\n",
    "\n",
    "### Fit by running the container <a name=\"fit-container\"></a>\n",
    "\n",
    "Enter the following command in the terminal to create and run the container,\n",
    "\n",
    "```sh\n",
    "docker run -v <path to host input>:/input \\\n",
    "           -v <path to host output>:/output \\\n",
    "           -v <path to host config>:/config \\\n",
    "           -v <path to host scripts>:/scripts \\\n",
    "           -ti edoardobalzani87/pgam:1.0 bin/bash -c \"python /scripts/fit_from_config.py 0\"\n",
    "# note the '0' at the end of the call denotes the first neuron is being fit. Can be looped to fit other neurons in the list\n",
    "```\n",
    "\n",
    "Delete unused containers by listing all the container with the command\n",
    "\n",
    "```sh\n",
    "docker ps -a\n",
    "```\n",
    "\n",
    "which will produce something similar to,\n",
    "```sh\n",
    "CONTAINER ID   IMAGE                       COMMAND                  CREATED              STATUS                     \n",
    "3f3644b47861   edoardobalzani87/pgam:1.0   \"bin/bash -c 'python…\"   About a minute ago   Exited (0)\n",
    "```\n",
    "\n",
    "Remove the container with,\n",
    "```sh\n",
    "docker rm <CONTAINER ID> \n",
    "```\n",
    "\n",
    "Fit results will be stored in \\<path to host output>."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "95a238b7",
   "metadata": {},
   "source": [
    "# HPC & singularity <a name=\"HPC\"></a>\n",
    "\n",
    "This protocol is currently tested on the NYU green HPC. Suppose that you have a script named \"run_all_fits.py\" that runs the j-th fit in the <a href=\"fit-list\">fit list YAML</a> with the following command,\n",
    "\n",
    "```sh\n",
    "python run_all_fits.py j\n",
    "```\n",
    "\n",
    "If you want to run in parallel N=100 fits on HPC, each fit with a max duration of 2 hours and memory requirement of 16GB you can follow the procedure below:\n",
    "\n",
    "\n",
    "1. Download and copy the singularity container <a href=\"https://osf.io/pcsav/\">pgam_1.0.sif</a> in one of your HPC folders.\n",
    "\n",
    "2. Copy all the configuration files, the inputs and the YAML with the fit list to your HPC folders.\n",
    "\n",
    "4. Make sure that the paths in the fit list YAML and in \"run_all_fits.py\" are appropriately set (paths should ponit to the file locations on HPC instead of the local folders).\n",
    "\n",
    "4. cd to the \"run_all_fits.py\" folder on the HPC \n",
    "\n",
    "5. Create a batch script (a text file) \"run-pgam.sh\" with following content,  \n",
    "   ```sh\n",
    "        #!/bin/bash\n",
    "\n",
    "        #SBATCH --job-name=fit_GAM\n",
    "        #SBATCH --nodes=1\n",
    "        #SBATCH --cpus-per-task=1\n",
    "        #SBATCH --mem=16GB\n",
    "        #SBATCH --time=0-02:00:00\n",
    "        #SBATCH --array=1-100\n",
    "\n",
    "        if [[ ! -z \"$SLURM_ARRAY_TASK_ID\" ]]; then\n",
    "            IID=${SLURM_ARRAY_TASK_ID}\n",
    "        fi\n",
    "\n",
    "        singularity exec $nv --bind /etc/passwd   path_to_PGAM_sif_image/pgam_1_0.sif /bin/bash -c \"python run_all_fits.py $IID\"\n",
    "   ```\n",
    "\n",
    "5. cd to the folder of the batch script and run jobs with the following command:\n",
    "```sh\n",
    "sbatch run-pgam.sh\n",
    "```\n",
    "\n",
    "\n",
    "The batch script configurations (#SBATCH --option=...) depend on specific dataset that one needs to fit and should be adjusted at need. For a complete list of options check https://slurm.schedmd.com/sbatch.html"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "32602a6d",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.13"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "autoclose": false,
   "autocomplete": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "hotkeys": {
    "equation": "Ctrl-E",
    "itemize": "Ctrl-I"
   },
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}