{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "aa06ba5d",
   "metadata": {},
   "source": [
    "# Spike-Train Distances and Synchrony\n",
    "\n",
    "This tutorial covers the pairwise spike-train measures of `braintools.analysis.spike`. Each one answers a slightly different question, and several have properties that are easy to misread, so every section checks a function against a case with a known answer.\n",
    "\n",
    "| Function | Input | Measures |\n",
    "|---|---|---|\n",
    "| `victor_purpura_distance(t1, t2, cost_factor)` | spike times | edit distance: shifting a spike by $\\Delta t$ costs $q\\lvert\\Delta t\\rvert$, adding or removing one costs 1 |\n",
    "| `van_rossum_distance(t1, t2, tau)` | spike times | distance between exponentially filtered trains |\n",
    "| `spike_train_synchrony(spikes, window_size, dt)` | spike matrix | fraction of spikes with a partner in another train within the window |\n",
    "| `correlation_index(spikes, window_size, dt)` | spike matrix | mean Pearson correlation of binned counts |\n",
    "| `spike_time_tiling_coefficient(spikes, dt, tau)` | spike matrix | STTC, a coincidence measure that is 0 for independent trains |\n",
    "| `phase_locking_value(spikes, reference_freq, dt)` | spike matrix | locking of each neuron's spikes to an external rhythm |\n",
    "| `burst_synchrony_index(spikes, burst_threshold, max_isi, dt)` | spike matrix | co-occurrence of detected bursts |\n",
    "\n",
    "The first two take arrays of spike times; the others take a `(T, N)` spike matrix plus `dt`. Pass time arguments as `brainunit` quantities and units are converted; the distance functions take plain numbers in one consistent unit."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "42e35f42",
   "metadata": {},
   "source": [
    "## 1. Setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "e85e3e05",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-10-15T08:09:47.245814Z",
     "iopub.status.busy": "2026-10-15T08:09:47.245814Z",
     "iopub.status.idle": "2026-10-15T08:09:48.806684Z",
     "shell.execute_reply": "2026-10-15T08:09:48.806684Z"
    }
   },
   "outputs": [],
   "source": [
    "import brainunit as u\n",
    "import jax.numpy as jnp\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "import braintools.analysis as ba\n",
    "\n",
    "rng = np.random.default_rng(0)\n",
    "dt = 1 * u.ms"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e941de0b",
   "metadata": {},
   "source": [
    "## 2. Victor-Purpura Distance\n",
    "\n",
    "The distance is the cheapest way to turn one train into the other. The cost factor $q$ (units of 1/time) sets the time scale: two spikes closer than $2/q$ are cheaper to shift than to delete and re-insert. For $q \\to 0$ the distance is the difference in spike counts; for large $q$ it counts the spikes that do not match exactly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "c1b552e5",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-10-15T08:09:48.808688Z",
     "iopub.status.busy": "2026-10-15T08:09:48.808688Z",
     "iopub.status.idle": "2026-10-15T08:09:48.812686Z",
     "shell.execute_reply": "2026-10-15T08:09:48.812686Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "identical: 0.0\n",
      "q =  0.01/ms: shifted by 2 ms -> 0.10  (min(5 * q * 2, 10) = 0.10)\n",
      "q =  0.10/ms: shifted by 2 ms -> 1.00  (min(5 * q * 2, 10) = 1.00)\n",
      "q =  1.00/ms: shifted by 2 ms -> 10.00  (min(5 * q * 2, 10) = 10.00)\n",
      "q = 10.00/ms: shifted by 2 ms -> 10.00  (min(5 * q * 2, 10) = 10.00)\n",
      "q -> 0, a vs c (5 vs 3 spikes): 2.0 (count difference 2)\n"
     ]
    }
   ],
   "source": [
    "a = np.array([10.0, 30.0, 50.0, 70.0, 90.0])  # ms\n",
    "b = a + 2.0  # every spike shifted by 2 ms\n",
    "c = np.array([10.0, 50.0, 90.0])\n",
    "\n",
    "print(\"identical:\", ba.spike.victor_purpura_distance(a, a, cost_factor=1.0))\n",
    "for q in [0.01, 0.1, 1.0, 10.0]:\n",
    "    d = ba.spike.victor_purpura_distance(a, b, cost_factor=q)\n",
    "    print(f\"q = {q:5.2f}/ms: shifted by 2 ms -> {float(d):.2f}  (min(5 * q * 2, 10) = {min(5 * q * 2, 10):.2f})\")\n",
    "print(\"q -> 0, a vs c (5 vs 3 spikes):\", float(ba.spike.victor_purpura_distance(a, c, cost_factor=1e-6)), \"(count difference 2)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b5723833",
   "metadata": {},
   "source": [
    "## 3. van Rossum Distance\n",
    "\n",
    "Each train is filtered with the causal kernel $K(t) = \\frac1\\tau e^{-t/\\tau}$ and the distance is $D = \\sqrt{\\int (f_1 - f_2)^2\\,dt}$. It is evaluated in closed form. For two single spikes a time $\\Delta$ apart,\n",
    "\n",
    "$$D^2 = \\frac{1}{\\tau}\\left(1 - e^{-\\Delta/\\tau}\\right),$$\n",
    "\n",
    "so small shifts cost little when $\\tau \\gg \\Delta$, and $D^2 \\to 1/\\tau$ (two unrelated spikes) when $\\Delta \\gg \\tau$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "2f8d4a5d",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-10-15T08:09:48.814137Z",
     "iopub.status.busy": "2026-10-15T08:09:48.812686Z",
     "iopub.status.idle": "2026-10-15T08:09:49.063788Z",
     "shell.execute_reply": "2026-10-15T08:09:49.063788Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  tau  Delta        D  closed form\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "    1    0.5   0.6273       0.6273\n",
      "    1    5.0   0.9966       0.9966\n",
      "    1   50.0   1.0000       1.0000\n",
      "   10    0.5   0.0698       0.0698\n",
      "   10    5.0   0.1984       0.1984\n",
      "   10   50.0   0.3152       0.3152\n"
     ]
    }
   ],
   "source": [
    "print(f\"{'tau':>5s} {'Delta':>6s} {'D':>8s} {'closed form':>12s}\")\n",
    "for tau in [1.0, 10.0]:\n",
    "    for delta in [0.5, 5.0, 50.0]:\n",
    "        D = float(ba.spike.van_rossum_distance(np.array([100.0]), np.array([100.0 + delta]), tau=tau, t_max=1000.0))\n",
    "        print(f\"{tau:5.0f} {delta:6.1f} {D:8.4f} {np.sqrt((1 - np.exp(-delta / tau)) / tau):12.4f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3f3390a6",
   "metadata": {},
   "source": [
    "The default integration end is the last spike plus $5\\tau$, which leaves out less than $e^{-10}$ of each kernel's squared tail; pass `t_max` to integrate further. Here is the distance as a function of the shift between the two 5-spike trains:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "2ca772db",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-10-15T08:09:49.065293Z",
     "iopub.status.busy": "2026-10-15T08:09:49.065293Z",
     "iopub.status.idle": "2026-10-15T08:09:49.626640Z",
     "shell.execute_reply": "2026-10-15T08:09:49.626640Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABEEAAAE1CAYAAAAI6v6kAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjksIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvJkbTWQAAAAlwSFlzAAAPYQAAD2EBqD+naQAAs1FJREFUeJzs3Qd8VGXWBvAnM+m9V0LvvfemIoiIooiiIhasn7r2wtrL6tpxbdhRKYpSRFC6gnSl917TE9J7ZvL9znszQxISSCDJnfL83bv3zmSSeYfAnTvnPe85LqWlpaUgIiIiIiIiInJwBr0HQERERERERETUEBgEISIiIiIiIiKnwCAIERERERERETkFBkGIiIiIiIiIyCkwCEJEREREREREToFBECIiIiIiIiJyCgyCEBEREREREZFTYBCEiIiIiIiIiJwCgyBERERERERE5BQYBCGyQ9OmTYOLiwuOHTum91CIiIiIiIjsBoMgRDbi6quvhre3N7Kzs6t9zC233AJ3d3ekpaVd0HPMnDkTU6ZMQUP5888/VbDGsrm5uaF58+aYOHEijhw50mDjICIiIsdz++23V7jO8PDwQOvWrfHCCy+goKBA7+ERkY1iEITIRkiAIz8/H/Pmzavy63l5efjll19wxRVX4JFHHlGPbdKkiU0HQSz+9a9/4fvvv8fnn3+OUaNG4ccff0SvXr0QHx/f4GMhIiIixyGBD7nGkO29995D06ZN8eqrr2LSpEl6D42IbBSDIEQ2lAni5+enAhVVkQBIbm6uCpYYjUZ4enqqWQ+9lZSUoKio6JyPGTRoECZMmIA77rgDH374Id555x2cPn0a33777UU/v9lsbtDZnpq8XiIiImoYrq6u6hpDtgceeABLlixB3759MWvWLCQlJek9PCKyQQyCENkILy8vXHfddVixYgWSk5PP+roERyRIIsGS6mqC/P777xgyZIh6nL+/v8q2sARVhg4dikWLFuH48ePWtFGZLbGQ55RZk4iICBVg6dKly1lBCnk++T4JYkhGSYsWLdQMzJ49e2r1Wi+99FK1P3r0qDWdtfxYLF566aWzAj1y+8EHH8SMGTPQoUMH9fyLFy+uMLb3339fZcnIn6n8eezatavCz5A/C9kqqzyOc71eCYRIum2PHj0QEBAAHx8fFez5448/avVnQUREZA9+/vln9Z64atWqs7722Wefqa9Z3m937Nih3lNlCaxcU0RGRuLOO+88azmv5X3+0KFD6vGBgYHqPVUmTSQD9kLIzxs4cCBKS0vPWnr7ySefWK8doqOjVdAkIyOjwmMOHjyIsWPHqjHL2Bs1aoTx48cjMzPT+phly5ap55Dx+vr6ok2bNvj3v/9t/Xp112mWZcKyt5DrkY4dO6o/M7lmkaXRLVu2VH/eQv68+/Tpo65p5HmWL19+QX8uRHSGq94DIKIzJMtDAg+zZ89WH/QtJGtCZjZuuukm9SZYFXnDlQsMeXOfPHmyemPeunWrChDcfPPNePbZZ9Ub+KlTp1SQQMgbt5ClNfImLBch8rzNmjXDTz/9pC5I5OLg4YcfrvBc33zzjcq+uOeee9SFRHBwcK1e5+HDh9U+JCQEF2LlypXWP6PQ0NAKgYvvvvtO1VWRCxsZ4wcffKCCLjt37lQBngtR1evNysrCl19+qX4nd999t3rOr776CiNGjMCmTZvQtWvXC3ouIiIiWyTLWeW6Qd5/5cN6ebLMVa4/5MO8JUggAQgJZkgwYffu3WpJrOw3bNhw1gTHDTfcoK493njjDWzZskW9v4aHh+PNN9+8oLFagg9BQUEVAi4vv/wyhg0bhvvvvx/79+/Hp59+ir///htr165VdctkgkPexwsLC/HQQw+pscfFxWHhwoXqekgCNPIarrrqKnTu3BmvvPKKui6Q6yf5GRcqPT1d/UwJtowbN06NS45lwkeWQN93333qWu7tt9/G9ddfj5MnT6oJLyK6QKVEZDNKSkpKo6KiSvv161fh/qlTp5bKP9clS5ao29988426ffToUXU7IyOj1M/Pr7RPnz6l+fn5Fb7XbDZbj0eNGlXapEmTs553ypQp6udNnz7del9RUZEah6+vb2lWVpa6T55PHufv71+anJx83tfzxx9/qMd//fXXpSkpKaXx8fGlixYtKm3atGmpi4tL6d9//60ed9ttt1U5rhdffFF9f3ly22AwlO7evbvC/ZaxeXl5lZ46dcp6/8aNG9X9jz76qPW+IUOGqK2yyuM41+uV31VhYWGF+9LT00sjIiJK77zzzvP+2RAREdmbm266qTQ8PFy9B1okJCSo9+VXXnnFel9eXt5Z3ztr1iz1nrp69eqz3ucrv29ee+21pSEhIecdj7xv+/j4qGsM2Q4dOlT6zjvvqGuMjh07Wq+B5D3c3d29dPjw4aUmk8n6/R999JH1OkVs3bpV3f7pp5+qfc73339fPUaerzqVr9MqXxfJ3kKuR+S+mTNnWu/bt2+f9Xpnw4YN1vvlOlDul59PRBeOy2GIbIjU+pDI//r16yukUMqSFsliuOyyy6r8PplxkUyEZ555RqVulleTuiG//fabmu2QrAYLmRGRgqY5OTlnpb5KmmhYWFiNX5dkqMjjJfVUZpKktolkvPTs2RMXQmag2rdvX+XXxowZg5iYGOvt3r17qzRSeY0XqqrXK78r6dRjqUsi2TpSL0Rek8xiEREROZobb7xRLZ8tv5xDlm3I+6B8zaJ81qpkUqampqo6HaKq90jJdChPlpfK0hnJujwfuaaQ92jZZBnJE088gQEDBqhaapZrIFlCIlkeklVhMJz5+COZnLJ8WJYLC8n0EJJ9W91yHMm0FfLz5XXXBcmwkes/C1n2Is/Trl07dQ1jYTlmhz2ii8MgCJENLokRlloesnzlr7/+Um+O8sH7XMtLLGmotSV1Qlq1alXhwkDIm6/l6+VJymp5iYmJFTZZXlOe1M6QQI0sY5E1r9IV5tZbb72gsVb1/OXJ66hM2uVVXpdbF88ngRxJh5XAkyztkQswuZAqv26YiIjIUUiHOgkUyPIXCzmWJaDyXmshEwOylFYmcCQgIu+PlvfSqt4jGzduXOG2ZRmLLBM5H3kPlmsM2WT5qly7SKCmfCDGch0jwYXyZDJD6pZYvi5jfOyxx9RyHFluK0tjPv744wpjlmCPBFnuuusu9frk+kyWCF1MQETqjlSetJI/59jY2LPuq+mfCxFVj0EQIhsjhTbbtm2rqpoL2csqEEtwxBZUrksSFRVVYSt/cSQ6deqk1uBecskl6lgqudckW8VkMtXo+WurLp5v+vTpqmaKFEuVWiBSe0UuwKT+SF3NDBEREdkSqX8hGZfz5s1T2Y9SL0NqYZTPArHU+Pjiiy9UhsfcuXOxdOlS9T4pqnqPrG6SR1sFe27yvXKNIZu8L0uBeZmQuffeey/oNb777rtqwkYKncqkjmTFSr0TmZSyXBOsXr1aZZfIhI48Vl7/5Zdfbr2OqO11RnWv/2L+XIioegyCENkgCXhIhXV5Y5WMEMlukE4v1ZEP4qJyF5TKqntTlk4qUg298oXJvn37rF8/F8sMjGWTmZPakBmfytXZq8pAqQl5HZUdOHCgQvHUung+Sf+V2SO5uJOLIHnNcgHWkO16iYiIGpp84JflLRJskCLq8oG8fBBEshTka7JEVwqRXnvttSpAIO+ZDUEmYx599FH8+uuvqghr+esYKYZaniyRkU51la9zZMLmueeeU8EOycaVYM/UqVOtX5fMWVmi/N5776mOcf/5z39UtqulQ5wlk6XytcaFXNcQUd1jEITIBlmyPmQZybZt286bBTJ8+HBVJVyqqlf+EF5+tkDauFaVhnrllVeqWZPyGRwyw/Phhx+qdaqVq8BXZpmBsWxyAVIbEsSRcUnQxyIhIUHNNNXW/Pnz1cWKhXRq2bhxI0aOHFnh+STAk5KSYr1v+/bttarsbpmdKf/nK88j9VyIiIgclbzPS5c0uWaQTWpvlV82WtX7o5BW8w1FOrtIq9n//ve/1jHL0pf//e9/FcYlmZxy/SH1yoTUIJHrn8oBEQl6SMcYy1Kfyiwd4SyPsUxOSRClfBaIdMghIv2xRS6RDZKLif79+6uiW+J8QRAp6iVtb2V9qmSMSBs1mYWQD/ZS2EtqV1iW2sgFi6x3lcdJgGP06NGq9etnn32m0kg3b96ssiYk00GCAnLRUt9t2GQ97dNPP61miyTtVMYs7eFkfXFti4xKUbSBAweq9ndyMSLjl3odTz31VIVCrTJ7I9kbkyZNUmuHZYZH0l1rUoRNSCs7yQKRMcvFk8wkyc+Qgq1STJaIiMgRSeH06667Dj/88IMqSvrOO++cdU0yePBgvPXWWyguLlbFymU5jLxPNhR535f2vJ988gn27t2r6oRMnjxZZaZIXZOrr75aZYXI1+V6aMKECer7JJvjwQcfVG1q5RpEAiLff/+9CuxIkXQhbXEluCHv/ZJBItcQ8nOkrodcfwi5npBCsPKcEjSRoJH8eVUOsBCRTi6iswwR1aOPP/5YtUHr3bt3jVuvLViwoLR///6qTay0dZXvlZZ0Fjk5OaU333xzaWBgoPr+8u1gk5KSSu+4447S0NBQ1UauU6dOZ7Vgs7SMffvtt2v0Giyt4M7Vas5i6dKlqp2dPHebNm1Uu97qWuQ+8MADZ31/+bG9++67pbGxsaUeHh6lgwYNKt2+fftZj5ef37x5c/V8Xbt2VW3nqmuRW9XrlbZ7r7/+unq8PE+3bt1KFy5cWG27XyIiIkexbNky9f4orWhPnjx51telVb20uZXrjYCAgNJx48aVxsfHq++R93YLy/t85Xaz1V3nVNcityqHDx8uNRqN6jHlW+K2bdu21M3NTbW0v//++1V7e4sjR46odr0tWrQo9fT0LA0ODi695JJLSpcvX259zIoVK0qvueaa0ujoaHUNIXtpHXzgwIGznn/YsGHqGkGe69///rf1z61yi9wOHTqcNX65lhg1atRZ91d3HURENeci/6dXAIaIqK5I9xfJoHn77bdVezwiIiIiIqLKWBOEiIiIiIiIiJwCgyBERERERERE5BQYBCEiIiIiIiIip8CaIERERERERETkFJgJQkREREREREROwVXvAVRmNpsRHx8PPz8/uLi46D0cIiIiqgVJMM3OzkZ0dDQMBvuaa+E1CBERkeNfg9hcEEQuPmJjY/UeBhEREV2EkydPolGjRrAnvAYhIiJy/GsQmwuCyOyLZeD+/v56D4eIiIhqISsrSwUSLO/n9oTXIERERI5/DWJzQRBL+qlcfPAChIiIyD7Z43ISXoMQERE5/jWIfS3WJSIiIiIiIiK6QAyCEBEREREREZFTYBCEiIiIiIiIiJwCgyBERERERERE5BRqHQRZvXo1Ro8erXrvSsGR+fPnn9Wb94UXXkBUVBS8vLwwbNgwHDx4sC7HTERERERERERU/0GQ3NxcdOnSBR9//HGVX3/rrbfwv//9D1OnTsXGjRvh4+ODESNGoKCgoPajIyIiIiIiIiKqI7VukTty5Ei1VUWyQKZMmYLnnnsO11xzjbrvu+++Q0REhMoYGT9+/MWP2MkcSj+EpceXwlxq1nsoRHUj7zSQcRyI6AAY3fUeDVGd8HXzxQ1tboC3m7feQyEiHZjNpSg2m1FiKkWxyYxiUylKym6bzHJs2Zutt+V7TJattBTmUu3nmEu1++R2qeV+tS9FqfVY+5rldqmMoey2HFu+pvbq9pn7RdkO8tUzxxXvL/84VPm9VavwPeUeVf7+M1+vpap+iC2pRWtwlxr8CJdyj6p4fxX3VfHccpf8DMuXZHfm2OWs77d8XdtrXy9/v6HsWPZyIHuD5f5Ke8v9shkNLjDI5gIY5T6DC1wN5fYuLnA1asdGg0Ht5bZ8n7vRAFejAW5GF7gZDOp7iBo8CHIuR48eRWJioloCYxEQEIA+ffpg/fr1VQZBCgsL1WaRlZVVl0Oye8+tfQ6703brPQyiupe6Se8RENUpU6kJkzpN0nsYRFSFgmITsvKLkVVQgpzCEmQXFCNbjgtKkFtUgrwiE/KLTOpY9nJbvqegxKz2hSVmFJbti2Qzle3LjiVoQUT1TwIjEhCR4Ii7qxEerrI3qNuebgZ4uBnhKZur3JZjA7zdXeHtboSPR9lebnsY4efpBj9PV/h7ulqPvdyMVQaUyLHUaRBEAiBCMj/Kk9uWr1X2xhtv4OWXX67LYTiM3OJc7D29Vx2Paz0OroY6/XURNazTR4DDK2WaS5tOkNkc2TfqBUR1qdXsCZEtOZF9Amvj1uKvuL8YBCFqIBJ0SMkuVFtydkG540Kk5RYiPbcY6XlFyMgrRkZ+EQqKGz6j1jKb7WowWD+4qRnv8jPgMvNtmSm37lE2a15xNt1gKJuBL5uJN1SYpS+btS/7mjaDX5ZJUH5mv/JMv2Ww5TIGKmcZWB5lvV3h7frcmQpn31+LTIlqHqr31UJpHSevnC9bpvJznnlMpcycarJ4LA/XsoHO3FkhM6hStpC6t+zxctlmyRjSYn0VM5AsAUCVoWQ+k5GkZTdpGU/yEC276UymkzUDyowKGVImyaYquy0ZVZVZvlf7N12CuibBlEBvNwT7uKt9kLc7gnzcEebrgXB/D4T7eSLcTzuW+yRLheyP7p+qJ0+ejMcee6xCJkhsbKyuY7IVu1N3q2UwUT5ReKHfC3oPh+jCbfwc2PiW9o7abjRw5bvA708Be+YDaUuArHzg2s8A/2i9R0pUa6eyT2Hk3JHYnrwdOUU58HX31XtIRA5BsjUOp+TicHIOTqbn4VR6PuJky8hHQmZ+lR+QzkWCBr4eZ2Z8ZZPbltlhy2yxbNoMsmUzwNNVO3YvN+sse5mFdrOk6suxpPKXpfVzNpno4khAxhIMkcBIcYlZBUrKZ2NZs7PUsWRwadlbBeWO88qyvfIKz2R7ZauMsDNZYbKXYI38XAmmylaTc0pUgBcaBcnmXbb3QrNQH7QM90WgN5d9O0UQJDIyUu2TkpJUdxgLud21a9cqv8fDw0NtdLbtKdvVvnNYZ72HQnRhZCpg5avAX+9qt3vdBYx8CzAYgXHTgK3TtWDI0dXApwOAaz4C2o7Se9REtdLIrxGa+DfB8azj2Ji4EZc1vkzvIRHZFfnwciApGztOZWJvQhYOp+SoLSnr3B9CJGtCZmLD/LQtvGwf4qPN3MoMrmUmV/YS8GBggsh+yL9XFVQ0Al4w1nvARQIlGfnFSM8tUplkp3O1bLK03KKybLMCLUCSVYjUnEIVkJGgrGwbj54+62eG+nqgZbgPWoX7oVWELzpEB6BDtL8KqJIDBUGaNWumAiErVqywBj0ks0O6xNx///11+VROwRIE6RLWRe+hENWeqQT49WFg23Tt9qXPAYOeqJhP2/1WoHFfYM4kIGE78MPNWqBk+GuAm5euwyeqjf7R/VUQZF3cOgZBiM7jRFoeNhxNw45TGdipAh/Zava1KhLYaBHmiyYh2ixrTNmMa0ygFyL8PVUghIioLgIukhUmm5xfzkeW+aTkFKoMtVNlmWqyP3k6H0dSchCfWaACJbJtOHImQCJZYm0i/dAlNhBdGgWge+MglTXCAK2NB0FycnJw6NChCsVQt23bhuDgYDRu3BiPPPIIXnvtNbRq1UoFRZ5//nlER0djzJgxdT12hybRyB0pO9QxgyBkd4rygJ9uBw4uAVwMwOgPgO4Tq35saCtg0jJgxSvA+o+Av78Ejq0Frv9K6yBDZAcGRA/ArH2zsDZ+rTp/82KG6Az5ELDucBrWHUrF2sOp6kNCZQFebujcKADto/3VrGmLMB80D/NV9xMR2Rqp2yOBWNl6NAk66+tSgFmW8h2SLSUH+xMl2y0DqTlF2B2fpbaZG7XHRvh7YEDLUAxqFar2UneEbCwI8s8//+CSSy6x3rbU87jtttswbdo0PPXUU8jNzcU999yDjIwMDBw4EIsXL4anJ3+ZtS20l16YDneDO9oFt9N7OEQ1l5sGzLoROPU34OqpLXtpU3VbbStXD2DEf4AWlwLz7gNS9gKfX6JlhPS+m0VTyeb1iuylilfH5cSp87csjyFy9myPhTvj8dvOBOyKq9j5T2ZCu8YGonuTIBX46BwTiNhgLwYPichhyPI7le0RG2i9TyZJJENk+8kMbD+VofZbT2SopX9zt8SpTbSN9MNl7cJxTdcYtI7w0/FVOC6XUktZYBshy2ekrW5mZib8/f3hrBYcXoBn1zyLrmFd8f2V3+s9HKKayTgBfH8dkHYQ8AwEbp4NNO5Tu5+RkwL88n/AwaXa7dYjgWs+BnxC6mXIRHVl0pJJ2JS4CZN7T8bN7W6Gs7Ln93F7HrstiM/Ix6IdCVi4Ix7bT2VW+Fq7KH8MaBGiZjl7NQtWHxCIiJydFG7dfDwdfx1MxZpDKWcFjSUgcnXXaIzuHI3YYG/dxulo7+N8B7JRlqUwLIpKdiNxFzB9LJCTCPg3Am6dC4S1qf3P8Q3TgicbPwOWPQ8c+B34tD9w7VSgxZksNCJbrAsiQZB18eucOghCzkXWxa86kIKv1hzFmkOp1vulVEe/FiG4qnM0Lm8foQoEEhFRRVIkVYLDsgFtVTHWvw6m4NftCVh1IBn7ErOxb/F+vLV4P3o2CcKt/ZpgVKcotua9SAyC2CgWRSW7cmwNMOtmoDATCGsHTJgDBMRc+M+TlOi+9wFNBwI/3wmk7ge+vxYY8C/gkucAV7YcI9szIGYApmyZogIhxaZiuBlZy4Acl7ScnLMlDt+sPYojKbnWU3evpsEY3TkKV3SMUp1aiIio5oJ93NUyGNky8oqweFciFmyPx/ojafjneLraJCBy58BmuLFXLLPqLhCXw9igvOI89JvVD+ZSM5ZfvxwRPhF6D4moent+AebcDZgKgcb9gJtmAV5nF4i6qCKrS58F/vlaux3dDRj7FRDSou6eg6gOyDn70tmXIq0gDV8N/wq9o3rDGdnz+7g9j72hSOvIL/46ghkbTyAzv1jd5+fhivG9YzGxX1OmaxMR1YOkrAL8+PdJfLf+mCquKvw8XXFLnya4Y0BTVaCVUOP3cebR2KBdqbvUxXSkTyQDIGTbpJPL7Nu0AEibUcCt8+o2ACLcvYGr3gdunK797PitwNRBwNYZUmGqbp+L6CIYXAxqSYyQLjFEjqTYZMa0tUcx9J0/8cmfh1UApHGwN14c3R7r/30Znh3VngEQIqJ6IkGOf13WCmuevhRvXNcJzcN8kF1QgqmrDmPwW39gyvIDqr4I1QyDIDaIS2HI5knwYeVrwKLH5QbQ43bghu8At/P3Vb9g7UYD960FmgwEinO14qlzJgH5GfX3nES11D9GC4JIXRAiRyE1P0Z+8Bde+nWPCn5Iob6pE3rgjyeG4o4BzZiOTUTUgDVEburdGMsfHYIvJ/ZU7XkLS8yYsvwghr23Ckt3J6ouNHRufNey4SBI51AWRSUbZCoBFj0KbPlOuz10MjDk6YZpYyt1Rm5bAKx5D/jjDWDXHODk38DYL2vfhYaoHvSL6qf2+07vQ2p+KkK9pNAZkX06kpKD/yzaixX7ktXtIG83PD68Dcb3imVRPiIiHRkMLhjWPkK10l20M0Gdq0+l5+Oe7zdjcOswlaXXIsxX72HaLL6D2RiJ3Fk6w3QJZyYI2RipzzH7Vi0A4mLQlqkMfaZhAiAWBiMw+Elg0lIgsAmQeQL4ZiTw55uAmWmApK8QrxC0C26njtfHr9d7OEQXfC3yw6YTKvtDAiCuBhfcOaAZ/nziEkzo24QBECIiG+Hi4qK6cK14fAgevKQl3I0GrD6QgiumrMb7yw7AZGZWSFX4LmZjTmafRHphOtwMbtYLaSKbkHca+H4MsP83wOgB3PA90PNO/cbTqCdw3xqg0w1AqQn483Vg2lVAxkn9xkRU1iVGsC4I2aPcwhI8+uM2PDN3p0qxHtQqFIsfGYwXRrdHgDc7HhER2SJvd1c8MaINlj46GJe2DUexqRQfrDiI27/ZpNruUkUMgtjoUpj2Ie3hbmQbULIRmaeAr68ATm4EPAOAifOBdlfpPSrA0x8Y+wVw7eeAuy9wYh0wdQCwe77eIyMnZimOKpkgUuSayF7sS8zC6I/WYP62eBgNLnhmZFt8e0dvtAxnSjURkT1oGuqDr2/vhfdv7AJPNwP+OpiK0R+uwfaTrKFXHoMgNoZFUcnmJO0BvrwcSN0P+McAdy4Bmmgf8mxGlxuB+/4CYnoABZnAT7cBCx4CinL1Hhk5oa5hXeHt6o3TBadVbRAie1j+Mvvvkxjz8VocSclFpL8nfrinL+4b0kKtOyciIvtybbdGmP/AADQN8UZcRj7GTV2PGRuPs2hqGQZBbLUoahiLopINOL4O+OYKIDseCG2j1eEIt9FlWsHNtQDNIOlY46LVLflsCJCg/ZsiaihuRjf0juqtjtklhmyd2VyKF37Zjafm7EBBsVkV1Fv0r4Ho1TRY76EREdFFaBvpjwUPDcTw9hEoMpnx7LxdeOInOdezhh6DIDYkrzgPB9IPqGNmgpDu9i4EvhujZVbE9gHuXAwENIJNM7oBl72gdZDxiwbSDgJfXAas+0iu9PUeHTmRAdFldUHiWBfkQrzxxhvo1asX/Pz8EB4ejjFjxmD//v3n/b6ffvoJbdu2haenJzp16oTffvutQcZrzwGQf8/bie83HFf1rZ8c0QbTbu+FEF8PvYdGRER1wN/TDZ/d2kMtb5TEvjlbTuHu7/5x+kAIgyA2ZHfabrV+PMI7ApE+kXoPh5zZP19rXWBMhUCbK4GJvwDedjQr2GwwcP9aoO1VgLkYWPosMON6IDtJ75GRkwVBtiVvQ24xl2XV1qpVq/DAAw9gw4YNWLZsGYqLizF8+HDk5lb/Z7lu3TrcdNNNmDRpErZu3aoCJ7Lt2rWrQcduTwGQyXN34oe/T6oL4/dv6IoHLmnJ5S9ERA7YQUaWN34/qQ+83IyqTsj/zdiCohLnnSBkEMSGsB4I6U7WCf7xBrDwUUAKOnafqHWBcfOC3ZGgzY3TgVHvAa6ewOEVWtHUg8v0Hhk5gVj/WMT6xaKktASbEjbpPRy7s3jxYtx+++3o0KEDunTpgmnTpuHEiRPYvHlztd/zwQcf4IorrsCTTz6Jdu3a4dVXX0X37t3x0UcfVfs9hYWFyMrKqrA5SwDkmbk78OM/ZQGQG7tiTLcYvYdFRET1aEDLUHx1e094uBqwcl8yHpq1BcUm5wyEMAhiQ7YnMwhCOjKbgIWPAKv+q90e/BQw+n+A0RV2S/K7e00C7lkFhHcAclO0jJDFk4GSQr1HR07SJYatci9eZmam2gcHV5+Rtn79egwbNqzCfSNGjFD3n2vZTUBAgHWLjY2FozOZS1X9j9n/nLIGQK7pygAIEZEz6N8iFF9M7Al3VwOW7E7CIz9uQ4kTBkIYBLERUqmXRVFJN8X5wOyJwOZpWlHRUe8Clz6rBREcQXhb4O6VQJ/7tNsbPtFqhaScv8YA0cUuiWFx1ItjNpvxyCOPYMCAAejYsWO1j0tMTERERESF++S23F+dyZMnqwCLZTt58iQcPgDy8w78vFkLgEwZ340BECIiJzO4dRg+m9ADbkYXLNqRgCd+2q7eH5wJgyA24mT2SaQXpsPN4Ib2Ie31Hg45k/x04PvrgH0LAaMHcMN3QK+74HDcPIGRbwI3zwa8Q4CknVr3GAn8sF0Y1QPpEOPq4qrO7yezHPvDdX2S2iBS1+OHH36o85/t4eEBf3//Cpsje3fpflUUz2hwwQfju+HqLtF6D4mIiHRwSdtwfHxzd7gaXDB/WzyembPDqdrnMghiIyxZIO1C2sHd6K73cMhZZMYB31wJnFgHeAQAt84F2l8Nh9Z6BHD/OqD5JUBJPvDrw1oWTN5pvUdGDsbHzQddw7uqYy6JuTAPPvggFi5ciD/++AONGp27O1VkZCSSkioWP5bbcj8Bf+xPxid/HlbHb1/fGaMZACEicmrDO0Tifzd1U4HxnzafUp3CnAWDIDaCRVGpwclSkK+GA8l7AN9I4M7fgaYD4RT8IoEJc4HLXwUMbsDeBcDUgcCxNXqPjBzMgJiyVrkMgtSKzEZJAGTevHlYuXIlmjVrdt7v6devH1asWFHhPuksI/c7u/iMfDz64zZ1fGvfJriuu423OyciogZxZacoPHtlO3X82qK92BPvHAXCGQSxETtSdqg9gyDUIE5s1AIgWaeAkFbAXcuAiA5wKgYDMOBf2msPbgFkxQHTrgJWvgaYivUeHTlYcVTpEFPMv1e1WgIzffp0zJw5E35+fqquh2z5+fnWx0ycOFHV9LB4+OGHVVeZd999F/v27cNLL72Ef/75RwVTnJlU/n9w5hZk5BWjU0wAnrtKu9glIiISdwxoisvahquWuQ/O2oK8ohI4OgZBbEBecR4OpB9QxwyCUL3b/zvw3TVAQQbQqBcwaSkQ2BhOK7obcO9qoOsEmX8GVr8NfDMSSD+m98jIAbQNbotgz2DkleRhW4o2E0/n9+mnn6pCpUOHDkVUVJR1+/HHH62PkZa5CQkJ1tv9+/dXQZPPP/9ctdX9+eefMX/+/HMWU3UGby/Zjy0nMuDn6arWf3u4GvUeEhER2RAXFxe8Pa4LIvw9cCQlFy8t2A1HxyCIDdidthumUhPCvcMR6cO1y1SPtnwH/HCLVguj1Qhg4i+Ad/UtJ52Ghy8w5mPg+q+12iin/gamDgJ2/KT3yMjOGVwM6BetLcdgl5jaLYeparv99tutj/nzzz8xbZp0tDpj3Lhx2L9/PwoLC1Ux1SuvvBLObNmeJHy++og6fvv6Lmgc4q33kIiIyAYF+7hjyo3dVGNIaaH+y7Y4ODIGQWwA64FQvSsty3BY8BBQatKyHsbPANx99B6Zbek4FrjvLyC2D1CYBcy9C5h3H1CYrffIyAFa5a6NY10QajgnT+fh8dla9tGdA5rhio6cZCEiour1axGChy5pqY6fnbcLx9Ny4agYBLEBDIJQvTKbgN+e0GpdiIGPAdd8BBjd9B6ZbQpqAtz+GzDkGcDFAGyfBXw2GIjbrPfIyE5ZMkH2nt6L1PxUvYdDTqBE6oDM2oqsghJ0jQ3EMyPb6j0kIiKyA/+6rBV6NQ1CTmEJ/jVrq6oT4ogYBNGZpPeyKCrVm+IC4Kfbgb+/lBV/wMi3gGEvyuI/vUdm24yuwCWTtWBIQCxw+ohWSHbN+4DZMd8MqP6EeoWq2iBiffx6vYdDTkBSmbefzIC/pys+urkb3F15uUdEROfnajTgg/HdEODlhu2nMvH+cq1upaPhu6LOTmWfwumC03AzuKF9SHu9h0OOJD8DmD5Wa/9qdNfqXfS5V+9R2Zcm/bTlMe3HAOYSYPlLwPdjgKwzxRiJatMlhnVBqL7J7N17y/ar40eGtUajINYBISKimosO9MKbYzur4y//OqKWVzoaBkF0ZukW0C64HdzlgypRXZAP6d9cCRxfA7j7ARPmAB2v03tU9skrCBg3Dbj6I8DNGzi6Cvi0P7DvN71HRnZYF0SCIOZSZhNR/Zn652Gk5hShaYg3JvRtovdwiIjIDl3RMRKDWoWi2FSKd5ZqgXVHwiCIzixLYTqHadE2oouWelBbupG8G/CNAO74DWg2WO9R2TdZPtT9Vq2VbmRnIP808MNNwKLHgeJ8vUdHdqBbeDd4uXqpzL/9px3vYoJsQ3xGPr74S+sG88zIdlwGQ0REF+zpK7SlvL9si8euuEw4Er472kpR1HDWA6E6cPJvLQCSeQIIbgFMWgpEMcBWZ0JbAXctB/o9qN2WWiufXwIk7dF7ZGTj3Ixu6B3ZWx2vjWeXGKof7yzZj8ISM3o3C8aIDhF6D4eIiOxYx5gAXN0lWh2/uXgfHAmDIDrKK87DgXSt2EzXsK56D4fs3YElwLejtSyF6O5aACSoqd6jcjyuHsCI/2hLjHzCgZS9wOdDgU1faK2IiarBuiBUn3aeysTcrXHq+LlR7eDCAthERHSRnhjeBm5GF/x1MBVrDjpOh7s6D4KYTCY8//zzaNasGby8vNCiRQu8+uqrqgsKVbQ7bTdMpSaEe4cj0idS7+GQPds6HZh1E1CSD7QcBtz2K+ATqveoHJv8Od+/Dmg1HDAVam2I5XeQm6b3yMhGDYjR6oJsTd6qguBEdUWusV5bpGWkXdstBp0bBeo9JCIicgCNQ7xxSx+tvtR/F++F2ewYn+nrPAjy5ptv4tNPP8VHH32EvXv3qttvvfUWPvzww7p+KsdZCsPWuHShJLj417vALw8ApSagy03ATT8AHr56j8w5+IYBN88Grviv1oHnwO9a0dTDf+g9MrJBjf0aI8Y3BiXmEmxK3KT3cMiBLNuThI1HT8PD1YAnR7TRezhERORAHrq0JXw9XLErLgsLdzpGh8Q6D4KsW7cO11xzDUaNGoWmTZvi+uuvx/Dhw7FpU9UXfIWFhcjKyqqwOVtRVAZB6IKYzcDiZ4AVr2i3BzwCjPkUMLrpPTLnIinnfe8H7l4JhLYBchKB768Flj4PlBTpPTqyIbI8wdIlZm0c64JQ3SgqMeON37W12ncNaqZaGxIREdWVEF8P3Du4ubX2lLzv2Ls6D4L0798fK1aswIEDWq2L7du3Y82aNRg5cmSVj3/jjTcQEBBg3WJjY+EsqavMBKELVlIIzLkT2DhVuy2ZCJe/rH0gJ31EdgLu+RPocYf8CwfW/Q/4ejiQdljvkZEN6R/DuiBUt2ZuPI6jqbkI9XXH/UNb6j0cIiJyQJMGNUOYnwdOnM5T7zv2rs6DIM888wzGjx+Ptm3bws3NDd26dcMjjzyCW265pcrHT548GZmZmdbt5MmTcAanck6pVomuBle0C2mn93DInhRkAtPHArvnAQY3YOxXWiYC6c/dGxg9BbhxOuAZCMRvBaYOArbNZNFUUvpE9oGriytOZJ/AyWzneL+j+mMyl+Lz1VpL3EeGtVbpykRERHXN290VjwxrpY7/t/IQsguKYc/qPAgye/ZszJgxAzNnzsSWLVvw7bff4p133lH7qnh4eMDf37/C5gwsWSDtg9vDw+ih93DIXmQnAt+MAo79Bbj7ARN+Bjpdr/eoqLJ2o7WiqU0GAsW5wPz7gTmTtAAWOTVfd190DtPaVq+LYzYIXZxVB5IRn1mAQG83XN+jkd7DISIiB3ZDz1g0D/XB6dwizNx4AvaszoMgTz75pDUbpFOnTrj11lvx6KOPqmUvdHY9EMvFMNF5pR4CvrocSNqptWa9YxHQfKjeo6LqBMQAty0ALn0ecDECu+YAUwcCJ1kQ09lZusSsjWddELo4lovQ67s3gqebUe/hEBGRA3MzGnBPWW2QH/8+adfdX+s8CJKXlweDoeKPNRqNMEsRR7JiPRCqlVObtfoSGSeA4ObApKVAFP/u2DyDERj8BHDnEiCwifb7+/oKYNVbgNmk9+hIJ5biqNIhpths3+mkpJ/4jHys3Jesjm/q01jv4RARkRO4qks0vN2NOJKai7+PpcNe1XkQZPTo0fjPf/6DRYsW4dixY5g3bx7ee+89XHvttXX9VHYrvyQfB05rhWMZBKHzOrgM+PYqIC8NiO4G3LkUCG6m96ioNmJ7AfetATrdoLUy/uM/wLSrgAzWhHBGUgcqyCMIucW52J6sBcSJauuHv0/CXAr0ax6CFmFsi05ERPXP18MVoztHW7NB7FWdB0E+/PBD1Rb3//7v/9CuXTs88cQTuPfee/Hqq6/W9VPZrd2pu1FSWoJwr3BE+kTqPRyyZdtmAbPGA8V5QItLgdsWAr5heo+KLoSnPzD2C+DazwB3X+DEOmDqAGDPL3qPjBqYwcWAvtF91TG7xNCFKDGZ8ePf2lKYm5kFQkREDeiGXlo310U745FlpwVS6zwI4ufnhylTpuD48ePIz8/H4cOH8dprr8Hd3b2un8r+l8KEd4ELW5pSVWSN3ZopwPz7AHMJ0PlG4KYfAQ/O9tm9LuOB+/4CYnpohVJnTwQWPAQU5eo9MtJhSQzrgtCFWLEvGUlZhQjxcceIDpxMISKihtO9cSBahfuioNiMBdviYY/qPAhCNS+KyqUwVCWpn7Pk38DyF7Xb/R8CxkwFXBlIdBhS10XqhAx8FIALsOU74LMhQAKXRjiL/tH91X5v2l7VLp3oQgqijusZC3dXXsoREVHDcXFxwY1l2SCz/7HPJTF852xgUkWXRVGpWiVFwNy7gQ2faLeHv6ZtlYoNkwMwugHDXgIm/gL4RQFpB4EvhwHrP9YCYeTQwrzD0DqoNUpRivXx6/UeDtmRk6fzsPpgijq+qbd2EUpERNSQru0WAzejC3acysTu+EzYG36yamBxOXFIK0iDq8FVFccjsirMBmaOA3b9DBhcgeu+0LJAyLE1HwLctxZoMwowFWlZQPL3IEfr+kCOvySGdUGoNmZtOqFWTA5qFYomIT56D4eIiJxQiK8HhrfXlmPOtsMCqQyCNDBLFki74HbwMHroPRyyFfKBd9oo4MifgJsPcPNsoPMNeo+KGopPCDB+BjDqPcDVEzi0HPi0P3Bwud4jo3rUP6a/NQgiWYJE51NUYramHt/CgqhERGQDBVLnbY1DQbEJ9oRBkAbGpTB0lrTDwFeXa/UgvEOB2xcCLS/Te1TU0KRIcq9JwD1/AuEdgNwUYMZYYPFkoKRQ79FRPege3h1erl5IzU/FgXStbTrRuSzbk4TUnCKE+XngsnYReg+HiIic2MCWoYgJ9EJWQQmW7E6EPWEQpIGxKCpVELcF+Go4kH4MCGoKTFoKxHTXe1Skp/B2wN0rgd73arelPswXlwEp+/UeGdUxd6M7ekb0VMfsEkM1MWPjcbUf3ysWbkZewhERkX6MBhdc36OROv5hk30tieE7aAMqKCnA/tPaBxkGQQiHVgDTrgLyUoGoLsCkZUBIC71HRbbAzRO48i2tLbJ3CJC0U+ses3ma1j6ZHMaAGNYFoZo5mpqLdYfTVNKYpSo/ERGRnsb1bKTel9YfScPxtFzYCwZBGtDutN0oKS1BmFcYIn20QjLkpHbMBmbeABTnAs2HArcvAnzD9R4V2Zo2VwD3r9P+jpTkA78+DMyeCOSxpaqjtcrdkrQFecV5eg+HbNjC7fFqP7hVGBoFees9HCIiIsj70aBWYXbXLpdBEJ3qgUh/ZXJS6z7S2uCaS4CO1wM3/wR4+Ok9KrJVfpHAhHnA5a9oXYP2LgCmDgSOcfmEI2jq3xTRPtEoNhfjn6R/9B4O2bBVB7S2uMM7sBYIERHZjht7atmJP/1zCiUmM+wBgyANaHsyi6I6NbMZWPIssPRZ7Xbf/9Pa4Lq66z0ysnUGAzDgYW3JVHALICsO+PYqYOVrgKlE79HRRZCAePkuMURVycwvxtaTGdZMECIiIlsxrH04/D1dkZxdiB1xmbAHDII0EGl/aM0ECWcQxOmUFAHz7gXWf6TdvvxVYMTr2odbopqSorn3rga6TgBKzcDqt4FvRmqFdcluDYjW6oKsjWN2D1Vt3aFUmMylaB7mg9hgLoUhIiLb4eFqRP8Woep47cFU2AN+Amsg8bnxSCtIg6vBFe2C2+k9HGpIhTnArBuBnbO15QzXfgYM+JfWEpWotjx8gTEfA2O/Ajz8gVObgKmDgJ0/6z0yukB9ovrA6GLEsaxjiM/R6j4Qlbf6oLYUZkhrZoEQEZHtGdhKC4L8dYhBEKpiKUzboLbwdPXUezjUUHJStGULh1cCbj5at48u4/UeFTmCTtcD960BYvsAhVnAnEnAvPuBwmy9R0a15Ofuh85hndUxW+VSVZmkq/YzCEJERLZrYEstCLL1RDpyC21/qTaDIA2ES2Gc0OmjwNfDgfitWpvT238FWg3Te1TkSIKaALf/Bgx5GnAxANtnAp8NBuI26z0yusAuMeviWBeEKjqUnIP4zAK4uxrQp1mI3sMhIiI6S5MQbzQK8kKxqRSbjtp+F0MGQXToDENOIH4b8NVw4PQRILCJVtAypofeoyJHZHQFLvm31mbZv5H2d07+7q2ZohXjJbuqC7IxYSNKpHMUUaWuMH2aBcPL3aj3cIiIiKos9G7JBlljB0tiGARpAAUlBdh/er86ZhDECRz+A5g2CshNBiI7aQGQkBZ6j4ocXZP+wP1rgPZjtPbLy18Evh8DZCXoPTKqgfYh7RHgEYDs4mzsTN2p93BswurVqzF69GhER0eri6v58+ef8/F//vmnelzlLTExEY4QBOFSGCIisoe6IGvsoDgqgyANYHfabpSUliDMKwxRPlF6D4fqkxSnnDEOKMoBmg3Wlir4Reg9KnIWXkHAuGnA1R8Cbt7A0VXAp/2B/b/rPTI6D6PBiH5R/dQxu8RocnNz0aVLF3z88ce1+r79+/cjISHBuoWHh8Ne5ReZsLEsrZhBECIismX9W4Sqvg/7k7KRnFUAW8YgSAPYkbJD7aXwncxKkYPa8KlWnNJcDHS4FrjlZ8DTX+9RkbORc0z3iVor3cjOQP5pYNZ4YNETQHG+3qOjmtQFiWddEDFy5Ei89tpruPbaa2v1fRL0iIyMtG4GO25FvvFoGopKzIgO8ETLcF+9h0NERFStYB93dIjWPvusPWzb2SD2e2VgR1gPxMGVlgLLXgAWP6Pd7n0vMPZrwNVD75GRMwttBdy1HOj3oHb77y+ALy4FkvboPTI6TxBkV+oupBek6z0cu9W1a1dERUXh8ssvx9q1586qKSwsRFZWVoXNJpfCtAnjJAoREdm8gS21rMW/bHxJDIMgDdDajkEQB2YqBubfD6z9QLt92YvAyDcBO555JAcigbgR/wEmzAF8woHkPcDnQ4FNX2jBO7IpET4RaBnYEqUoxYaEDXoPx+5I4GPq1KmYM2eO2mJjYzF06FBs2bKl2u954403EBAQYN3ke2zJ6rIgyOBWXApDRES2b1BZXZC1h1LV52BbxU9q9Sw+Nx6p+alwdXFVhe/IgRTmaMsMts8CXIzANZ8Agx7TliMQ2ZKWw4D71wEtLwdMhcBvTwCzbgJy0/QeGVXTJYZ1QWqvTZs2uPfee9GjRw/0798fX3/9tdq///771X7P5MmTkZmZad1OnjwJW3HydB4Op+TCaHBB/7KK+0RERLasR5MgeLgakJRVqFq82yoGQerZ9mQtC6RtcFt4unrqPRyqK7mpwLejgUPLtQKUN/0AdLtF71ERVc83DLjlJ+CK/wJGd+DA71rR1CN/6j0yKqd/jLYkZn38epueQbEXvXv3xqFDh6r9uoeHB/z9/StstmL1QS0LpFtsIAK83PQeDhER0Xl5uhnRu1mwzS+JYRCknu1IPVMUlRxE+jHg6xFA/BbAKxi47Veg9XC9R0V0fpKl1Pd+4O6VQGhrICcR+G4MsOxFoKRI79GRzKBE9ICn0RPJ+ck4mHFQ7+HYvW3btqllMvbIshSGXWGIiMieDGh5ZkmMrWIQpIEyQVgPxEEk7gS+Gg6kHQICYoFJS4FGPfUeFVHtRHYC7lkF9LhDKhcBa6cAX8vf68N6j8zpeRg90COyhzpeF+fcXWJycnJUEEM2cfToUXV84sQJ61KWiRMnWh8/ZcoU/PLLLyrzY9euXXjkkUewcuVKPPDAA7A3xSYz1h5KsxZFJSIishcDy4IgG46kqfczW8QgSD0qKCnAvtP71HGXcAZB7N7R1cA3VwI5SUB4B2DSMq0DB5E9cvcGRk8Bbvge8AwE4rcCnw0Gts1k0VRbqQsS79x1Qf755x9069ZNbeKxxx5Txy+88IK6nZCQYA2IiKKiIjz++OPo1KkThgwZgu3bt2P58uW47LLLYG+2HE9HTmGJajfYMTpA7+EQERHVWPsof/X+lVtkwraTGbBFrnoPwJHtSduDktIShHqFItonWu/h0MXYPQ+Yew9gKgKaDATGzwC8AvUeFdHFa381ENND+/t9fI3W7ejQCuCq9wBPfvjSMwiyJWkL8kvy4eXqBWcknV3OVRdl2rRpFW4/9dRTanMElnogUmXfYGCxbSIish8GKejdIgQLdySouiC9mmo1QmwJM0HqkaU1bufQznBhxxD7tfFz4Kc7tABIu6u1dqMMgJAjCYgBblsAXPq81ulo18/A1IHAyU16j8wpNQtohkifSBSZi/BP4j96D4d0sIr1QIiIyEFa5doiBkEaIAjCpTB2SmYgV7wC/P6kVjeh113AuGmAG7v8kAMyGIHBTwB3LgECmwAZJ4CvrwBWvQ2YTXqPzqlI0NySDbIu3rnrgjijzLxi7IrLUscDyy4iiYiI7MnAVloQX5bDZBUUw9ZwOUw9kRReaxCERVHtj6kE+PVhYNt07falzwGDntC6azT0UEwmFBfb3smDKjIajXB1dbX/rK/YXsB9fwGLHgd2/gT88Rpw5A/gus+BgEZ6j85pDIgZgDkH5zh9XRBb0xDn431x6YjxMyLczxP+bkBBQUG9Pp+zc3NzU+dvIiKqOzGBXmgW6oOjqbnYcDgNwztEwpYwCFJPEnITkJqfClcXV3QI6aD3cKg2ivKAn24HDi4BXAzA6A+A7mc6EDR0d4RTp06dc1082Q5vb2/VjtPd3R12TWqBXPcF0HKYFgw5vhb4dABw9YdaDRGqd32i+sDoYsTRzKNIyElAlK99tnl1JA11PjYUluClS8Lh5WZQHXGofkngulGjRvD19dV7KEREDtcl5mhqrloS4xRBkLi4ODz99NP4/fffkZeXh5YtW+Kbb75Bz57O00rUkgXSJrgNPF25fMJu5J0GZt4AnPobkN+bLH9pM1K3GUe54JYP1mFhYfafYeDA5EORdKZISUlRH1patWoFg8HOVxvK37cu44FGvYA5dwHxW4DZtwLdbwOueANw99F7hA7N390fnUI7YVvKNpUNcn3r6/UeklNryPNxYmY+3POLVWX9MD9eP9T3uVvO2/K7lfM2M0KIiOpOvxYh+H7DcWw+kQ5bU+dBkPT0dAwYMACXXHKJCoLIxcLBgwcRFBQEZ8KlMHZIaiBMHwukHtBaht48G2jcR7fhSMq1XKDJvyEvL+fsDmFP5HckadXHjx9XARFPTwf58BLSQqsT8ufrwJopwJZvgRPrgbFfAVGd9R6dQ+sf018FQaQuCIMg+mrI83FJVglcXF3g5+MNT087zyqzA/I7PXbsmPodMwhCRFR32kb6qf2h5ByYzaU21e2szoMgb775JmJjY1Xmh0WzZs2qfXxhYaHaLLKytGJg9m57cllnmDB+SLALSbu1AEh2AuDfCLh1LhDWBraAGSD2w+6zP6rj6g4Mewlofgkw714tUPjlZcCwl4G+9+tSK8cZSHHUT7Z9gg3xG1BiLoGrgStYHf18LIGWwmKtELGnGz+QNwS+xxIR1Y/Gwd5wdzWgoNiMU+n5aBziDVtR51fsCxYsUMtexo0bh/DwcHTr1g1ffPFFtY9/4403EBAQYN0kgGLvCkoKsO/0PnXMTBA7cGwN8PVILQAS1g6YtNRmAiBENqX5EOC+tUCbUVrL6CWTgRnjgJxkvUfmkKSeVIBHALKLs7ErdZfew6EGUGwyw1Raqj6Yy4UjERGRvXI1GtA8VFs+fTA5G7akzt9hjxw5gk8//VStrVyyZAnuv/9+/Otf/8K3335b5eMnT56MzMxM63by5EnYu72n96KktAQhniGI8Y3Rezh0Lnt+Ab6/DijMBBr3A+78HQjg70xPr732Glq0aKG2Z599ttrHJScn44orrlDnmo4dO2L16tXWrw0dOlRloHXt2lVt77//foXv/eeffzBypD61XuyeTwgwfgYw6l2tbs6hZcCn/YGDy/UemcMxGozoG9VXHbNLjHOQ2TLh4WqA4TwZCi+99JJNdI6RgrEjRoxAaGgoAgMDz/lYs9mMhx56SJ3fpV7cRx991GDjJCKihtc6QlsScyApB7akznNr5Q1OMkFef/11dVsyQXbt2oWpU6fitttuO+vxHh4eanMklqUwkgXCNEsbtukL4LcnJQEZaHsVMPZLwM12a29ImnR+WZp0ffFyM+r6d1YCGbNmzcKOHTtUu1mpL9S/f3+MGjXqrMc+88wz6Nu3LxYvXoy///4b1157rSpKKnU5hAQ+xowZU+XzzJs3r9qvUQ3I35FedwGN+wNzJgHJe4AZY4G+DwDDXgRcHeucrveSmCXHlmBd3Do80PUBvYdD9XwuTs8tQkGxCQFe2nnsXF5++WU88sgjutcfknOuFMMPDg5WAehzmT59Ovbs2YMDBw6oiS+5RpQach06sIseEZEjahWudd46mJTt2EEQaQ/Zvn37Cve1a9cOc+bMgdMVRQ3nUhibJO0NV74G/PWOdrvHHdqstsG211/LRXf7F5bU63PseWUEvN1dz7vkTS545cJXMjG+/vprlVnRtGnTi37+H3/8Ebfeeit8fLTUuTvvvFMFRaoKgsyePRuHDh1Sx7169UJ0dDRWrVqFYcOGnfd55DUsW7ZMFcOTTBGZmVy0aBGys7Mxbdo0/Pzzz/jjjz9QUlKCH374QWWaSIHn22+/Xc16SrD3mmuuUVkrTi2iPXD3SmDZC8Cmz4ENHwPHVgNjvwbCWus9OofQL7qf2u9K24XMwky1PIYc/1y86slzBxPuu+8+tR80aJAq5rl06VIsX74cH3zwgSrMLOcoOT+NHj1aPU6CExIwsQR/r7/+elx11VXqnHaxZCLr0ksvVefTmpzj7777bjVmCZrceOON6hxf+VzKczMRkWNoVZYJcjA5x7GXw8jM7f79+yvcJxH/Jk2awFlmiCxBkM6hLIpqc0wlwIKHzgRAhk4Grnrf5gMgtkKWoNxxxx0qqCnZGpLOnJaWVuVj5TxgWY5SeZOfUZUTJ05UOFdIYEXuq0yeUyr5R0ZGVvtYyRTp1KmTusiWZXoWcsHs7+9v/V6ZjezRowe2bNmivkfSuq+++mps27ZNZa/JbKuQtG350LB9+3bs3LkTjz322AX9GTocyZ668m3gph8Ar2AgcSfw+RBg87dawJEuSqRPJFoGtoS51Iz1Cev1Hg41EFkOcy6SXSv++usvda6SGmxy7tqwYQO2bt2KX375RQUbyheerykJrFR37pZWwRejpud4wXMzEZH9axXhW6FDjMNmgjz66KMqfV2Ww9xwww3YtGkTPv/8c7U5g8TcRKTkp8DVxRUdQpneaVOK8oCf7wQO/A64GLTsj553wl7IUhXJ1Kjv5zgXucDu3LmzNdtr0qRJaqauKm3atFEXq3r4/vvvVZFlCUp+/PHH6gJZUrCrWgojqeSW27KUz9fXV6Vni969e2PGjBnqePDgwXjyySfVbOOQIUNqlHHiVNqMBO5fB8y/DzjyJ/Drv4DDK4DRHwBeztUiva71j+6PQxmH1JKYK5peofdwnF59nYvNpaXYG5+NUpQiyLv2rXFlOeAtt9yCU6dOqeWEp0+fVve1bdu2Vj9HAiu2gOdmIiL710Q6xBgNKosyLiMfscHejpkJImnp8iFD0hslTfHVV1/FlClT1BuzM7BkgbQObg0vV9utL+F08k4D312jBUCMHsAN39tVAERIrQ5ZqlKfW23rgZzr8ReSCdK4cWMcP368Qkq03FdZSEiIushPTEys8rGWLlMyvgcffFBlglgyVubPn18hCFK+JpGkaJdfXy+3Je1ajB07FmvXrlXBHcvMI1XiHwVMmAdc/gog7Vyl8PCnA4Hj6/Qemd3XBbEUR5XAHjnmudjoYoCHmwE+Hq4X1Blm/PjxuOuuu1QdNglAS9DAUjhVzpflszjOVVC1PjNBanqOFzw3ExE5SIeYMG2Z+wEbqgtS55kgQt6AnPVNyFoPhK1xbUfmKa0DTOp+wDMAuOlHoIm2zp5qp1+/fiqAsW/fPjW7KPVAZP15XWWCSGvtBx54QGWXyEW7/HzpgFDdYyUlXL4uhVHj4uLULKBcGEvAIyIiQj1Olu7IsQROEhIS1GyhdJSpLVlGIx0NJk6cqGYhJeONqmAwAAMeBpoO0oqmnj4CTBsFDHoCGPI0YKyXtx2H1j2iOzyMHkjOS8bhjMNoGdRS7yFRPZCCqMLTtWYFqv38/NSSEUtHlvT0dNUVy1KAVG5byNLFjRs3qoCBZIesWbNG1QVp6EwQOW9/8cUXai9jlxohCxcuvKifyXMzEZHt1wXZl5it6oJc1k67Ptcbr0brGIMgNiZ5rxYAyY4H/KKBCXO0Yo50QcLCwlRgQjqxuLu7q8KoElyoK1K8T2p4SC0PIceWgKoUX33hhRfw22+/qdtvvvmmKqIqAQ0Zi1z0S7HW3NxcVUhV1sIbDAbVtlEKoQpZJy9ryi+EFOST55DnkuJ7ljX5VI2Y7sC9q4Hfnwa2zQBWv6Utk5EuTEHOUSOqrni6eqJnRE+VCSIbgyCOqaCkLAhynmWJFo8//jguv/xyeHt7q8KoUhRVAhsSFJFCpeUzLJ566inruVU6sfTp06dOxy7LJFNSUpCVlYVGjRqpZSuyLFFIBomct6V4tZyzJWgt520J9Ej9Dsv5/kLx3ExEZNtal3WIsaVMEJdSG8utlTfQgIAANUMgxQvtSaGpEH1n9kWJuQS/X/c7Gvk10ntIzu34emDWjUBBJhDaBrh1LhBgP78TSVeWGTuZ2dO7BeK5SJChrrrD1DcJ2kjXAFlf7sy/swa382dg4aNAYRbg4a8VI+5U9Sw0Ve273d/h7X/eVvVBPrv8M9gye34fr27sDfFv+2hqLrILihET6IUQX7aZbig8bxMR1a/FuxJw3/Qt6BQTgF8fGmgT1yB1XhPEme1N26sCIMGewYjxjdF7OM5t70Lg+zFaACS2D3DnYrsKgFD9WLx4cb0FQOgcJOBx3xqgUW8tECLLZOb/H1BoOzMCtm5AjFYXZHPSZhSUVF/PgRxgOUwNM0GIiIjsqU3uIRvqEMMgSD0thaltgUmqQ/98Dcy+FZAPCm2uBG6dD3gH6z0qh5WammoXWSCkM1kCc8fvwOCntO5MskTms8FA3Ba9R2YXmgc0R4R3hMo4lEAIOZYSkxnFJrM69nTjpRkRETluhxhbwHfaOsR6IDqTlV1/vKGl3Zeage4TtS4w7rbRionI6UlR1EufBW5bCPg30oqmfnU5sGYKYNY+AFLVJLBuyQaRuiDkWApKtL//cpFolOLCREREDtgh5mCybWQB8522jkhple3JDILoxmwCFj4CrPqvdltmm0f/j50oiGxR0wHA/WuA9tcA5hJg+Yva8rWsBL1HZtOkHohYF8eWw46GS2GIiMgZlsQcSMqBLWAQpI4k5SUhOT8Zri6u6BDaQe/hOJfifGD2RGDzNJkvBUa9q802c0kSke3yCgLGfasFK928gaOrgE/7A/t/13tkNqtvVF8YXAw4nHkYibmJeg+H6iUIwssyIiJyPK3KOsQcZBDEsWxL2ab2rYNbw8vVS+/hOI/8dOD7a4F9CwGjB3DDd0Cvu/QeFRHVhAQqe9wG3LMKiOwE5J8GZo0HFj2hBTepggCPAHQM7aiO18ZxSYwjKSi21ANhJggRETme1hFlQRAuh3EslqUwnUM76z0U55EZB3w9EjixHvAI0Frgtr9a71ERUW2FtQbuWgH0e1C7/fcXwBeXAkl79B6ZzRkQzbogjriclsthiIjIkbUM97NmgthChxgGQerIjpQdat8lnPVAGkTyPuCr4UDKXsAvCrjzd6Bp/fadpobx2muvoUWLFmp79tlnq33cHXfcgdatW6NLly4YMGAA/v777wpfT05ORocOXJpmN1w9gBH/ASbMAXzCgeQ9wOdDgU1faEWPqUJdkA0JG1RLdrJ/0hXGXFqqit+6u9b8suyll15CQYH+7ZKPHTsGo9GIrl27WrfDhw9X+Viz2YyHHnpInd9btmyJjz76qMHHS0REDa9piDfcjC420yGGVSPrgLQs3HNam7FkUdQGcGIjMPMGoCADCG2tfWgKbAyHJx8Ei/Pq9zmkNoOOtVRWr16NWbNmYceOHXB1dVXBjf79+2PUqFFnPfbaa6/FF198oR63cOFCjBs3Tl2MW/zyyy+4+mpmBtmdlsOA+9cB8+8HDi0DfnsCOLQCuOZjwCcEzk6Ww/i5+yG7KBu7Uneha3hXvYfkfOr4XFyQXwyX4jyVBWIodq3xufjll1/GI488Ak9PT+jNz88P27Zpy4LPZfr06dizZw8OHDiAzMxMdOvWDZdccgkD1kREztAhJtQX+5Oy1ZKY2GB9u3cyCFIH9qbtVTNywZ7BaOTbSO/hOLZ9vwE/3wGUFACNegE3zwa8g+EU5KL79ej6fY5/xwPuWgur6ixYsABPP/003NzccMUVV+Drr7/GP//8g6ZNm1700//444+49dZb4eOjjeHOO+9UQZGqgiDlAxx9+/ZFXFwcSkpKVFBEzJ8/Hy+++KI6lhlWyTCRsSclJWHKlCnYu3cv5syZoy7EJZgydOhQpKSk4JZbbkFCQoL6nh49euCbb7656NdFteQbBtzyE7BxKrDsBeDA71rR1Os+A5oPhTNzNbiqAqnLji/Duvh1DII4wLnYH0CnWp6L77vvPrUfNGiQysJYunQpli9fjg8++ABFRUUq40LOeaNHj1aPk/ObBEzGjBmjbl9//fW46qqrcPvtt6MhyTn+7rvvVmMODg7GjTfeqM7xMtbyJKAtGSWSNbJo0SJkZ2dj2rRp+Pnnn/HHH3+oc/0PP/yAjh074uDBg+p15ORIirUZ11xzzVk/j4iI9NcqoiwIkpSDS9tG6DoWLoepA9tTyuqBhHVWH5yonmz+FvjxFi0A0voKYOIC5wmA2AhZYiLLUCR4INkaks6clpZW5WP3799fIT26/CY/oyonTpxAkyZNrLclsCL3nY9c+F955ZXWAIhcMO/btw+9evWyPsbX1xcbN27EV199hQkTJiAqKkoFb15//XU8+eST1lnKZs2aYefOner1vfvuu7X+M6I6IufSvvcDd6/UMr5yEoHvxgDLXgRMxXBmrAtCU6dOVfu//vpLZWCEh4djxIgR2LBhA7Zu3aoy4STYUFhYWOufLYGV6s7dJpNWu6Sy3Nxcdb7t3r07XnnllWofV5tzvASoJRC9ZcsWPPPMM+r1SfBbXu9tt92mMmGELKmRgM727dvVufuxxx6r9WsmIqL619qG2uQyE6QOgyBcClOPqcer3wH+KJvZ6ToBGP0BYHSyv76SHi2zg/X9HOcgF9idO3dG+/bt1e1JkyapmbqqtGnTpkbp0RdLAhezZ89WS2ksfv/9d5WlUj4oKTOOomfPnuqCffz48ep279691UyiJaPk/fffx+OPP47Bgwern0E6k64x0j1myWStDfbaKcDR1cDYL4GQFnBGA2K0IIgsh8kszFRdY8h+z8VyMVhYYlLrpf083c48Ry0dPXpUZbKdOnVKBYRPnz6t7mvbtm2tfo4EVmpDAsqSiSeBGHlOOddKAPmpp57CxZBlPpbMFTlvSyBbls5YztszZsxQx3KulkC2ZIIMGTIEw4YNu6jnJSKi+m2Te8gGOsQwE6QOMAhSj8wmrSaAJQAy6Angmo+cLwAi5AO9pEfX51bLTKZzZT5dSCZI48aNcfz48Qop0XLfuVKrZTZw2bJliIg4k1Y3b94868WzhWXdvKRhV74tqdWiX79+KnDTp08fzJ07V81sVjejSQ3I3VsLfN7wPeAZCMRvAT4bDGyb5ZRFUyN9ItE8oDnMpWZVIJXs91xsdvVGoYsnSt284entf8HnYiGB3bvuugu7du1S5zEJGlgKp0pQpPy57FwFVWubCeLh4aECIEKWuMgyxuoCKbU5x8vPtZDzdPnaJ+XP22PHjsXatWtV4N2SFUJERLanVVkmyMFk/TvEOOEnybqVmJuI5LxkGF2M6BDCwl51qrgAmHs3sHeBXHUCI98C+tyj96icmgQJJIAhS01kdlHqgcj687rKBJHipg888IDKLpGLdvn50gGhKpL98dxzz6l18OUvomU869evx3fffXdBM6kxMTG44YYbVBaIXNjL7GJAAGfabYK0wI7pDsy9Bzi+Fph/H3B4BTDqXcAzwOm6xBzJPKLqgoxoOkLv4dAFkgyQUpTCaHCBq9Gl1sVIZclIYGCgup2enq6W81ky5OS2hSxdlOWAEjCQ89yaNWtUXZC6yASRZZJBQUGqTpQsv5EAshQ8re4cLzWYZC9jl0C2FLa+GJLJJ91mJk6cqDJEpJg2ERHZniZlHWLyirQOMXoWR2UmyEXalqJ9yGsd1BreF5C+StXIzwCmj9UCIEZ3YNw3DIDYgLCwMBWYkM4s0ppWLj5DQuquY4cU75NU6k6dOqFdu3a4/PLLrbN6Ur9D6n5YSNq3zGZKETzLLKXUJ1m5ciUGDhyoLshr688//1Rr0OVnyYX022+/zQCIrQloBNz2K3Dpc4CLEdj5EzB1EHCyYotkZ1kSszZuLUqdMBvGURQUm9VeOsPUtqaYLNuTc6ScryQQIbWRJLAhAQipC1I+OCxLU6SgqJxbJ0+erLLd6ooEVOQ55T1BaoJERkZWaG8u44uP15YPSeFrCaC3atVKZdpJ/Q4Z08WQYqnyM2QM8v5hqZdCRES2xa2sQ4w4lKxvXRCXUhu7esrKylIfOmSGwN9faqbbtrf+fgvf7/ke49uMx7N9z7zp00XIStACIMm7AXc/4KaZQLPBcDbyAV9m7GRmzxZaIFYnNDS0zrrD1AXpmiBrwqub5axP9vI7cxgnNwFzJgEZJ7SAyCWTgYGPAQZtyZMjyy/Jx8BZA1FkLsL8a+ajRaDt1Eext/fxmoy9vv5tJ2TmIyW7ECG+HogJ9Kqzn0s1x/M2EVHDeXDmFizckYDJI9vi3iEtdLsGYSZIHXaGoTqQcgD46nItAOIbAdzxm1MGQOjCySygHgEQ0kFsb+C+NUCncUCpCVj5GvDt1UDmKTg6L1cv9IjoYc0GIftUaMkEceXlGBEROb5W4WfqguiJ77oXochUhL1pe9Vx17Cueg/H/kk6+9cjgMyTQHALYNJSIIrBJVuXmppqM1kg5ISkFsh1XwBjpgLuvsDxNcCnA4A9UkvIOZbESF0Qsk8lZYXhXI28HCMiIsfXOkJbDnMwSd8OMXzXvQh70vag2FyMYM9gNPJrpPdw7NuBJcC3o4H800BMDy0AEsQP1kRUA1JLoetNwL2rgejuQEEGMPtW4NeHgaI8OHJxVPFP0j8oKKm+24e9krbXo0ePRnR0tKqXMX/+/BrV9ZG6FNJZRIqBTps2DbbMUh2/ljVRiYiI7FIrSxAkOUfXmmYMgtTRUpjaFjSjcrZOB2bdBJTkAy0v14oe+oTqPSoisjchLYA7lwADH9U6Sm2eBnw+BEjYAUfUMrAlwr3DUWgqxJakLXA0ubm5qtjmxx9/XKPHS12HUaNG4ZJLLlGdqR555BHVMnbJkiWwVaayC0CDgdcQRETk+JqE+FToEKMXBkEuwo4U7cK6S1gXvYdin+Tib/U7wC8PaOv5u9wM3DQLcPfRe2REZK9c3YFhLwET5wO+kUDqAeDLy4D1n2jnHAciwXdLNsjaeMerCzJy5Ei89tprqhtVTesBSXHLd999V3WXevDBB1V9oPfff7/a75GWrlJErfzWkEzWTBAGQYiIyDk6xDQL1T7rHUzSry4IgyB1kAnCIMgFMJuB358GVr6q3R7wCDDmE8BY+7amRERnaT4UuH8d0OZKwFQELJkMzBgH5CTDkQyIZl0Qi/Xr16vOUOWNGDFC3V+dN954Q1WRt2yxsbFoKJIGbC4LzBmZCUJERE6iVYSlOKp+dUEYBLlAibmJSMpLgtHFiA4hHfQejn0pKQTm3Als+ky7fcV/gctf1tb1k1P7+++/0b9/f3h7e2PMmDHnfGxycjKuuOIKtGrVCh07dlT1A8qTtr0yk0xOzCcEGD8TuPIdwNUTOLRMK5p6aDkcRd+ovnCBCw5lHFLvS84sMTERERERFe6T25LdkZ9fdcrt5MmTVRs9y3by5MkGzwK50OUwL730kmrvqredO3di8ODBaNu2rToX33nnnRX+vDdu3KiWNbVu3RqXXnop4uLiqv1ZkvnTokULtT377LMN9AqIiKghTezbBJ/e0h2ju0RDLwyCXGQWSOug1vB289Z7OPajIBOYPhbYPQ8wuAFjvwL63q/3qOyCzBrmFefV66ZngSIRFRWFKVOmnDN93eKZZ55B3759cfDgQXzzzTe4+eabUVxcbP36vHnzzhtIIScgwdXedwN3/wGEtwdyk7Vz0JJntYCsnQv0DETH0I7qeH189RkPVDUpoOrv719ha6hzcW5xHgpM+SgyF6CgJL/W5+KXX37ZJoIgnp6e+Oijj7Bv3z5s375d1XJ588031dfMZjNuueUWdV4/cOAArrzySlWrpSoSyJ41axZ27NiBPXv2qFouixYtauBXQ0RE9a1P8xCM7BSFqAAv6MVVt2d2oKKoVEPZicD064GknYC7HzB+upayTjWSX5KPPjP71OtzbLx543mDegsWLMDTTz8NNzc3lYnx9ddfq6yLumiT26hRI7XJBfD5zJ49G4cOHVLHvXr1Uh0kVq1aZU2Hl3EuW7YMx44dQ9euXfHQQw+pC+rs7GzVMeLnn3/GH3/8gZKSEvzwww9qBlMCKrfffjtycnLUxfs111yjZibJAUS0B+5eCSx9Hvj7C2D9R8DR1cD1XwOhrWDPpC7IztSdqi7Ita1qVj/DEUVGRiIpKanCfXJbAhteXl4OdS6+77771H7QoEEwGo1YunQpli9fjg8++ABFRUXq/CXnLumuI4YOHaqCD5bAsNRKueqqq9T57mJJNp6FjEXOx7t27VK3N2/eDFdXV1WsVtx777147rnnVPBGgifl/fjjj7j11lvh46OtFZeMEgmKSLHbqurhyOuT87z8jiXIsnfvXsyZM0dl9HzxxRfqNaekpKggTEJCgvqeHj16qKA5ERE5N2aCXCDWA6ml1EPAV8O1AIhPOHDHIgZA7JAsQbnjjjvUhabM1kkLyrS0tCofu3//fhV8qGqTn3Ex5Dkl60M+9FhIEObEiRPqWIIZ8sHH8nW5KJaL3y1btqgMEqkTcPXVV6sOErfddpuaURUymykfDGQ2U1K8H3vssYsaJ9kYNy9g1DvATT8AXsFA4g7gs8HA5m/tumjqgJgB1kwQk9kEZ9WvXz+sWLGiwn0SCJX7HY0UgRV//fWXOo+Fh4er89qGDRuwdetW/PLLL7j77rtV4dfaksBKdeduk+ncf78kC+TLL79UAWQh5+QmTZpYv+7n56fOzfHx8Wd9b+XHlj+nV8XX11cttfnqq68wYcIElUkoAfnXX38dTz75pHrM9OnTVbFcOZ/Le5YUzSUiImImyAUoMhVhb9pedcwgSA3EbdYKEualAUHNgFvnAcHN9B6V3fFy9VKzg/X9HOciF9idO3dG+/bt1e1JkyapDIuqtGnTRl2c66HyUhiZcbTc7tmzp7p4tsxM9u7dGzNmzFDHsq5dLp4lE2TIkCFnFVkkB9FmpFY0df59wJE/gV//BRxeAYz+APAKgr3pFNoJfm5+yCrKwu603Q6ToSj/Di3ZXpYWuHJOCQ4ORuPGjVU9D6kv8d1331mzIySQ+dRTT6ksgpUrV6qMsbpeUlFX5+LsghIcP50LbzdXNA/zqdW5uCry5yNZD6dOnVLZF6dPn1b3Sa2O2pDAyoWQDJQbb7wRw4cPr3FHn4shz2U5p0vwZfz48dZzugTChSyZlOWVjz/+uDq/S/YiERFRvWeC/Pe//1UpiNWtAbVHe0/vRbG5GEEeQYj1a7hK8nbp4HJg2lVaACSqKzBpGQMgF0j+HUl6dH1u8hy1HVN16jMTJCQkRF3kSyFEC1n2Ih+MxPz58ysEQWTdf/l07fJp2HJblsSIsWPHYu3atSqAY8kKIQflHwVMmAcMexkwuAJ7fgE+HQgct78uK64GV/SJ6uNwXWJkVr9bt25qE5KZJccvvPCCui1LHMpnCsiMvwQ8JPtDCnHKrL9kJUiGhC2ei90NnvA0esHbzeuiz8VCggB33XWXWooiwSIJ9lpqhsj5snwWx7lqiVxIJohk5klQQrIxZEmOhZyTjx8/br0tyxElM0+WL1ZW+bHlz+lVsZzH5Rxe+bblnC5ZQPJn0adPH8ydO1ct1TlfNgsRETk+1/ru9PDZZ5+pmWNHsj35zFKYC7lQcRrbZgELHgTMJUCLS4Ebvgc8fPUeFV0EuaCUAIYUwJPZRakHIrN/emSCjBs3TqWES4cEOdfIjLBkb8gHI5lBLr9OvaZk9lC6EkycOFHNJkqnGnJgBgMw8BGg2WBgziTg9BFg2ihg8JPA4KcAo/0kS/aP6Y/lJ5arIMh9XbR6EfZOajqcq0Co1Pap6ntkOYg9MF1ke1xZWiIBhcDAQHU7PT1dBYIsy0DktoUsXZSlIxLoleyQNWvWqLogdZEJIgEHCcBIhs7nn39e4bpIliFKgETqL0n2nVwTSp2SyvVALOf0Bx54QGUXStBG3l/k/H4x5LXGxMTghhtuUFkgsmxI3h+kHTIRETmvessEkTcZScuU4lRBQfaXXlyjeiDhXApTJbmwWzNFSzWXAEinG4CbfmQAxAGEhYWpC1NJdZaZVgkaSFZGXZHsESmMKjO+0hlAjj/55BPrrLB0FrCQ7gPr1q1TwQ4p7icX/VKsVdbCS72PCyHFUjt16qRmm2VW07LunhxcTHfg3tVA11uAUjOw6k1g2pVA+plZaVs3IFqrC7IjZQeyi7L1Hg7VgLmsRa7xAidTZInH5ZdfrjI0pF6TZGBIYEPOXxIIKp9FIUuEJBAh5zdZRiSZEXVFCppKloUlc0fGI8EMYTAY1Ln54YcfVi1yFy5cWKH7l5zT5fssASw578oY27Vrp17bxWbj/fnnnyoQI2OSoPbbb7/NAAgREcGltJ56YkqxQZkVkDc7eWOTNyCp3l2ZFO0qX7grKysLsbGxanajJq3q9DDsp2FIykvCV8O/Qu+o3noPx7aYzcDSZ4EN2gdX9H8IGPaKNuNKtSLpyjKLJTN7Vc2a2YrQ0NA66w5TF2S2T7oGyDrxhmYvvzM6h50/AwsfBQqzAA9/4Kr3gU5Vz5jbmtHzRuNY1jG8P/R9DGuiXz0beR+XD5q2/D5e27HXx7/thMx8pGQXItTXA9GB+rUJdHY8b9s5+RgjBaFNRdomk2+mYsBcXLY3afdZt7LbEvAuNWl7uU9+jhxD9qXl9pVYg5YugIuhbJNjy20jYJDNtezYoB0b3LS9ZBhabhtlcy/b3Mr9bCKq72uQesn1lXaT0oVBUtTP54033rB2ZrAHibmJKgBicDGgY2hHvYdjW0oKgfn3A7vmaLeH/wfo/6DeoyIns3jxYr2HQPZMAh6NegJz7gZObdKWyRxeCYx8E/Dwg613iZEgiLTK1TMIQjVjMl/cchgimyYBBLkuLMoBCrO1fVFu2T5POy7OLdvnA8V5QHFB2V5u5wMlshUCJXJ/gbaX2ya5r6hsLxOp9tvdqwIVGHEHXD0AV0/A1b1sb7ntqXU5s+zV5l1u7w24W/a+gLuPloVtOVZ7X05MEtVHEOTkyZMq7VEKk9Ukoi5pmeXbUFoyQWyVpBqL1kGtVfEyKlOQBfw4ATi6SjuJj/kU6DxO71FRA0hNTdV7CER1K6gpcMfv2rKYv94Bts0ATqwHxn4JxPSAreof3R8z9s7A2ri1qpYGa1bZyXIYBkHIVoMYEqDIPw3kp2tbQaa25WeUHcs+S8ucU/tsoFDuz9KCHZJxoQfJwJDMChVUsGReWDIzyjbLsSV7w5rJIf8eXSrty/48tIOKmSIqo8SyL8suUZklsrfctmSnmMplqBSfPW65TzYJDtUndz8tqG/ZPP21zEfPgLJNjgO1TTqmeZXt1X0BdlUvi6g6df63ePPmzWptavfu3a33SSXu1atXq24LsvTFUsnb0rWhfOcGu6kHwta4Z+QkA9PHAok7tAjzjd9rhVCpTtTTijWqB/xdORC5yLv0WaD5UGDu3VrR1K+GA5c+D/T/l03OpPWM6Ak3gxsSchNwNOsomgc013tIDqcu/42byn6UgcEqXTnNeVs+eOemArkpZVsqkCdbWrnt9Jm9BD2q+qB+IcpnJlj3cp8P4FZ2bMloUFkOciyZD5LhUJYBYd0kK8IDMMpelpF4VFxWogIetnd+Pu8yHvn9qH2lLBdL9oslG0YyZdSxZMmU3S6yZM/klh2XZdpYsm8KZZ+jBWSE1I2S7ULLR0kgxCsY8JYtpOw4BPAJAbxDAZ8wwCdUu883XPud8zxHjh4Eueyyy7Bz584K90k3Cekk8fTTT1cIgNgjBkEqSTsMTL8OSD+mnfgm/AxEa+0M6eJY/q1I9xUvL64Xtwd5eXlqLwVayUE0HQDctwb49WFg7wJg+Yva8phrP9Pa7NoQyU7sHtEdGxM2Yl3cOgZBbPx8zOUwtsHS4cwur0/lg7RkYWQnAtkJ2j4nCchO0vZqSwZyk7WgxoWQbAr5sKuyAizZAWUZASozwJJF4F+WWRBQLsvAsvzCDv9s65sEBSTYrrIqvBvm74oESizLk6ybJYsn60ymj3XLKMsCkuyfdC1wIixfTz9as+eWYJavBEbCtaCI2iIBvwjALwrwlX2k9nVmmVADqfO/adKyrWPHirUyfHx8VAeJyvfbmyJTEfak7VHHncMcq+3vBYnbAswYp80kSPr4hLlASAu9R+UwpEWgt7c3UlJS1IdqqbJPtjuTKAEQyYKTdpV2eTFN1ZMPADd8B2z5Dlj8jLbsb+oA4JqPgTYjYWtdYiQIInVBJrSfoPdwHEZ9nI+LCwtRajKhpMiIApeyGVpqUGazWf1O5Xcrv2ObIh9a5YNn5kkgKx7IPKXts+LK9vFa0KM2SydkuUf5mXo1W182Y2+Z1Ze9ZZZflkBIVgZn8e2f/A5V9o23FoS4EJKtIkuhZImUyhYq25fPJiqfZZSbpv39lNouGSe07ZxjNGgBEf9oLTjiH6NNNgTEAgGNtE2CJwyUUB3g36Ja2Ht6L4rNxQj0CERjvzOt55zSoRXAj7dqJ7eoLsAtP1/4SZWqJOv5o6KiVNX648ftp1WnM5MASGRkpN7DoPq6gOxxG9C4HzDnTiBxJzBrPND7HuDyV7Q0bhupC/Le5vfwT+I/KDQVwkNSxckmz8cJmQVaNki2B9xdGeTWiwS0pJ1wg9fQkSCHfICUTNqMY1pLbvmQKMEOCXxknKx5gEOyMvzkg2OE9iFRzbRHlG1lM+8yyy5BDU6o0IWSZUeS0SFbTclyHJWNlKLtLdlJOZK9JFlLCWfuk+U6KqMp4dyBEvm7LgGRoCZAYONyWxPtfhknkS0EQaRPuyOwFEWVpTBOXXBux2ytC4wUepL18jdOt/muCfbK3d0drVq1sqbrku2S2WFmgDiBsNbAXSuAFa8A6z8CNn0OHFsDXP81EN5O79Gpot1hXmFIyU/BlqQt6BfdT+8hOYy6Ph8/+L+/UFBswnd39kZMEAut6/l7rbdMSymMmXVKqylk3Y5qgQ/ZZGnC+UjwIkBmxC1b9Jm9LCGQGXOZ3SeyRVL3JbiZtp2L1EeRQIk1yynhTNZTZtyZjCipUyP/pmQ7ueHsnyNFbiUQIs8nWepBZc8d3AIIbs5/K2TFTJBaYD0QAOs+BJY+px13vF7rAiOFqajeyMVZTTotEVEDkaJ8I/4DNL8EmH8fkLwH+Hyodl/PSbqmjkuAXgIfCw4vwLr4dQyC2Oj5uMRkxuHTWjAlyN8Xnp58H7Vrsmwl9RCQegBIOwikHgTSDmkBDylweS4yq61mtC2z2pL6H6sdS7BDCoMSOTqpG6OCepFAzJnmGmcFFaW+jWRLWZbXVN7U0hvJqqomY8+/kbZ0X22tgNDWQGgr7d8cs6ScCoMgteDUQRA58Sx7Xpv5FH3/Dxj+H54wiMh5tRoG3L8OmP9/wKFlwKLHgUMrgas/1Krk61gXRIIgUhfkcTyu2zioejmFZ1qH+nnyUsxuSI2DlL1Ayj4gZf+ZvaTzn6uwqMxIyyy0dSubpZYPXgxyENWMfOawBEoa9ax6iZn8W1TZVkfL7Y9ojRyk0Ksli0Rqe5UntW9CWmpBkfC2QHh7IKyt9u+UhX0dEt95aygpNwmJuYkwuBjQMdS+C7zWmrTq+uUBYOds7basf5cWkc68JIiISMha+5tnAxunap1j9i8Cpm7Rusc0H6LLkCT7wwUuOJh+EMl5yQj3Zr0mW5NdoAVBPN0McDNyMsHmSEtSCXAk7QGSdmnZXkm7zx3skKwN9SGqlbaXWWaZbZZABws5EtU/+VxiCZI0qSILUmrwSIaW2g5rmVsqe+uw1lY4cYe2Ve5sE9ZGC4pEdAAiO2mbFA4mu8azci2zQFoFtlJtCJ2GtM+aPVFrCSm916UbQpfxeo+KiMi2Zqf6/R/QdCDw851aOvx31wADHwEuebbBi7QFeQahfUh77E7brZbEjGk5pkGfn84vq6BY7f08WcBPd9IeVAody4efhLIPQRIAkbpnZ3HRlqlI/R/5YCQzxbKX2WPWRiOybaoDUm8gtnfF+00lWo2e1P1nMrySJeNrv7a8JmGbtlVexhbZEYjsDER3BaK6arVIOEFsNxgEqSGnXAqTkwLMuF77hy+Bnxu+19K/iYjobFGdgXtXAYsnA1u+Bda8DxxZBYz9ssHbh0uXGBUEiWMQxJYzQbgUpoEVyweaHUD8FiBuCxC/VQtaVsUzUJvxtcwAyyZBDw/fhh41EdUnydQKlSyulkDbURWLtcqSGskEk02CpZIZJgGTbCneGg8cXHrm8dJiWoIhEhSJ7q4t2ZGsFLJJfPetbWeYcCcJgsj6ue+v09bSyT/qm38CGvXQe1RETiO/yISiEjO83I1sn2lvlfCv/h/Q4lLg139pH7Y+Gwxc+Y6WRddAs0QDYgbgi51fYH3CepjMJhi5ptmmZOVrmSD+zASpP1IfQK5lTv0NnNwEnNqkLW+RNpxVFUuUIKbM6lr2nNUlcm7yvmkJjrS/umL2mCyPk4BIwnZtslgyR/LSgMMrtK38uUU+P8X01IIi0d0ANy9dXg5VxCBIDRSbirEnbY/zZILEb9MyQKRVlaR9TpinnQCIqF6UlpbiVHo+Nh9PV9uWE+nYm5AFc6n2dTejC7zcjPDxcFVb81AftI3yR7tIP7VvHOwNo4EX6zalwxjtgmfuPcDxtVoXGbkwGvUu4BlQ70/fOawzfNx8kFGYgb2n9zpfLSsbx0yQeqpfJh9Gjq8DTmzQgh7yoaSqlrPSfUJmatW+G+ATqseIicgeefprNUfK1x0pLgCSd2ufoSTDTDLNpIiyFGHdI9svZwolS6ZI475AbF9tz/OPLvjuWwNyAVlkLkKgRyAa+zWGQzvyJ/DDLVrv+ohOwISfmcpFVMcKik3YFZepgh1a0CMDKdnVt1EsNpWi2FSCrLIPToeSc7B0z5kCfRIgaR3ph/ZRfmgb6Y+2ZcGRAC/OMutKZpJv+xX46z3gzzeAnT9pM9JjvwJie9XrU7sZ3NAnsg9WnlyJtXFrGQSxMdllNUGYCXIR5EOHBDok6CGBxpN/a+v3yzO6a+npUgOgUS8tMCkFTJnhQUR1Sbo8xUjGR4+KdRUlKBL3D3BKtr+14sqylw0fao+TIspSU6zZIKDJQMAvQreX4UwYBKlFPRCZWXNx5DfOnT8D8+4DzMVA00HA+BkNMmNJ5OiSsgq0YIdkepxIVwEQCWyUJ9keHaID0KNJkNq6Nw5CiK878opMyCsq0faFJmTkF+FgUg72JWZhX2I29idmI7/YhO0nM9RWXkygV1lARAuOtIvyQ9MQH7iyG0XDptMOeVLrFDNnEpBxHPh6BHDJZGDgY/Xaek+WxEgQRIqj3tvl3np7Hqo9ZoJcACleKDOs0tpSthMbAVOl4LFXMNCkf9ksax8gqgvg6qHXiInImUmxZAlsyGZZoif1RCRT7eQG7Rwm2SJSl0i2zd9oj5NCy/I5rNlgbWMnmnrBd98acIqiqOs/AZZM1o47XKu1d+SFA1GtFZvMaimLJcNDAh9xGZVmJ+U9ztcD3RsHagGPJkHoFBMAT7ezPxAHeBnOyugY1CrMemwyl+JYWq56zn0J2So4sjchWz2nZVuxL9n6eA9XA1pH+FmzRSxLaoJ93Ov8z4LKkZno+9YACx8Ddv0MrHwNOPwncN1nWsZIPRVHtbyHZRdlw8+d3StsRXYhgyA1Ih8YDi0HDq0Ajq0BCrMqft0vCmgyQAt8yF4+PEi3JiIiWyMT6cHNtK3rTWfa9kpQRM5vx1YDibvOtO795yvAxaAt2ZM6Y80v0TLaXHm9Vhf47uvsQRCzGVjxErD2A+12n/uAEW/wIoKohk7nFlkzPGS//VQGCorNFR4j5TokE6N7Ey3o0aNxMGKDveoks0xqgbQI81XbVZ3P3J+ZX6yyRPYnZmFPWXBEbktGyc64TLWVF+7noQVFovzQTpbURPmheagvi7LWJcmsk04xLYcBvz0BHF8DfDoAuPrDikXX6kgjv0Zo4t8Ex7OOY1PCJlzW5LI6fw66uOUwbJFbxRIX+XdxUAIfy8/u3CL/hmSGtPlQoNkQILQVl7YQkf2SLI+2V2qbNSiyHji6WitRIO164zZr2+q3AXdf7fzX+gqg1XAunbkIDIKcR1JuEhJzE2FwMaBTaCc4FFMx8MuDwI4ftNuXvQgMfJQXFETVkKyLg8nZ2HI8Q2V6bD2RjiOpuWc9zt/TFd0aB6Fn2dKWzrGB8PVo2NOtZI/0bhasNguzuRQnTudpWSOJWmBE9sfT8pCcXYjk7BSsPpBSYYmOBFfaRZ2pMyJBkjBfD8deGlif5M9NZoAkM2TOXVr3mNm3Aj1u1wLQ7t51ng0iQZC18WsZBLEhlvo+zAQBkJsGHFwC7FsEHF4JFOed+ZqLUVvW0mqYNgsqy1vY6YiIHDooMupMq97MOC0YIufGI39oxZ73LdQ2IVkirUYAba7Q6h/x2qzG+O57HjtStda4rQJbwdutbi9OdVWYA/x0mzbTIhcZMhPZ7Ra9R0VkU7IKirFNlrSUFTCVY0sae3ktw32tS1tkkwwKgw12a5ExNQ31UdvITlHW+3MKS1SWiAqKJJzZy2vVgiXZFX5OiI+7tc6IBEckSCJ/BlUt56FqhLQA7lwC/Pk6sGYKsHkacHw9cP1XQGTdBdwHRA/ArH2zVF0Q6ULE4JWt1QRxc95lLnsWAPt/19bGl5bLnpPCpZItJZvU0mFtMiJyVgEx2ucz2SR7P3E7cGApcGCxNokidZJkW/VfICAWaDda2yR4zIDxOTEIch7bkx1wKUxuKjBjnPaPRwI7474FWg/Xe1REupIPiMfS8qwtamVpy/6kbFXHqjxvdyO6xpbV8mgchG6NAxHobd/rMyVLxRLAKf/nIfVEpL6ILKnZm5itMkiOpeYiLbcIaw+lqa38spxm0rq3LCjSrixIEhXgyQ/e1ZF1vcNe0lJb594LpO4HvrgUuPwVbWliHfy59YrspTrFxOXEqYyQpgFN62ToVFfLYZzoMiz9OLBnPrB7nnbRXp4E/tqMAtqM1LI9eM4gIqpIShVI5odsQ58GspOAQ8u0gMihlUDmSWDDJ9rmEwa0uRLoMAZoOhgwOtF7TQ3xT6QWnWEcZvbl++uA04e1Kuq3/KS1jCNyMvlFJuw4lWGt5SFFTKW+R2VSu6NHY614qQQ95EO+M3RXkcBFoyBvtV3ePqLCn5ssCbIuqUnIxt7ELGTkFavWvbIt3JFQYWlQ+QKs8ufXJtIP3u58+7GSIMj964AFDwL7fwMWP6MVghzzKeB7pgjuhZAMxu7h3bExcaNaEsMgiG3IyneSIEhWglYIeNdcbeLFQor9SSHTdldrgY/AWD1HSURkf6QeSLcJ2lacDxz+A9grGXa/AbkpwJZvtc0nHOh4HdDxeu0zH4PMioO/+16cYlMx9qTtcZxMkIQdwIzrtR7VAY2BW+dqRcWIHJxkNcRnnmlTK5kee+KzUGKumOYhRUA7xwRYAx5SyDTcz1O3cdsiL3cjOjcKVFv5P9+krMKKtUYSsnE4JUfVPth09LTaLOT9V1r1qjojZUVYpRhroyAvm1xG1CB8QoDxM4G/vwSWPqfN7nzaH7j2U21ZwEXoH9NfBUFkScwt7bjs0ZaWw/g74nKYolxg70Jg+yytla1lqYsl8CEd6CT4cZEBPiIiKuPmdabAqtR8PPYXsOcXbctNBjZO1bagpkDHsUCXm4HQlnBmLqVy9WpDsrKyEBAQgMzMTPj7++s6lp0pO3Hzbzcj0CMQq29cbd8p3VJleNbNQFE2EN4BmDAH8D9TE4DIkRSVmLE7PtO6tEX28iG9sgh/D+uyFtm3j/aHhyvXUNaVwhITDidr7XtlaZHsZXlNas7Zvwvh425UWSKqEGtZ9ojcdrq6CUl7gDmTgGQtCI9+DwKXvXDBbcv3n96P63+9Hl6uXlgzfg3cje5O8z5uq2Nv9/xi5BebsPrJS9A4xAHqjcmlpLR43DZDq/VRXK5gtKxN7zQOaH8N4Buu5yiJiJxLSZFWUHXnz1rx6fLn5sb9tCyS9mMAD189R6nL+zgzQWq4FMauAyCy/nbuPYCpCGgyEBg/A/A6M4tLZO+SswtUx5atZQGPHXGZKhBSntSs6BDtX5bhoQU9olmvol5JQEkCS7KVl5JdqAqxls8cOZiUg9wik1qWJFt5kiEiGSOqfW/ZkpomIT7qd+qQItoDd68Elj4P/P0FsP4jLZB9/dcXlL3XOqg1Qr1CkZqfiq3JW9Enqk+9DJtqpthkVgEQh1gOI51dJPAhhX1lma2FzDZ2uQnofAMQ3FzPERIROXftsdYjtE2y9KR+yPYftMYY0opXtt+e0jL0uk/UOtc5yXWxnb/7NkwQxK6Xwmz8DPj9aZmm0dJPr/sCcGN6P9mvEpNZZRXIshYt0yNDtX2tLMjbTQU6pFWtBD66xAawDoWNCPPzUNvAVqEVPhgeTc0tV2tE2ydkFuBUer7alu9Nsj7ey82I1pI1opbUlLXvjfRHgLeb46S2jnoHaHkZMP//gMQdwGeDgZFvAt1urdVFigT6pFXugsMLVF0QBkH0lVO2FEb42mMQxJL1sfkbYO+v2gSLcPfV0qy73uJUF9JERHbB3Uc7R8uWlaAtWdw6XQtgb5uubRGdgN53a9l77g6QpXgOdvju23DsOggiFykrXgHWvKfd7nUXMPIttksiu5OZV4wtJ9OtQY/tJzNUxkB5cq3dJsKvLOARiJ5Ng9E0xJtZHnbEzWhA6wg/tV1T7v6MvCIVDFHBkbL2vRIEk5l0+bsgW3mS3WMpwGpZUiNda+y2mK0UjZSiqfPu1eorLHhIm8EZ/QHgdaabz/lYgiDr4tbhsR6P1euQqWb1QCSQJ3/v7UZxAbDzJ2DDp0Dy7jP3R3UFet6hXVh7+Ok5QiIiqgn/KGDQY8DAR4ETG4Ct3wO75gBJO4Ff/wUse16bcOl5JxDSAo6IQZBqJOclIyE3AQYXAzqGdoRdkYI4vz6spaiKS58DBj3BWRmyeWZzKY6k5pQVMNU6t0i3kcr8PFzRtbHWpla2LrGBjllgkFT74b7NQ9RmYTKX4nharjU4srcsOCLZIlIAV7aV+5IrFLxtFe571pKaEN8Lq7Ghy8XKrfOBdf8DVr6qFTo7tRkY+wXQpH+NfkS/6H5wgQv2p+9Xy2JkeQzpI8ve2uPmJAN/fwX885XWcUC4eWtLXXrcrrVrJCIi++PiAjTpp23DX9M+O0qBdukmKktxZWs9Ehj4CNC4LxyJnbwD65cF0jKwJXzcfGA3ZL3XT7cDB5dqldhltlDWeBHZoJzCEjWTr7I8TqRj64kMZJa1jiyveaiPyvKwBD1ahvs6bj0IOi/53TcP81XblZ2iKny4PGAJjJQtqZHaI5I5tDs+S23lyZIcCYZYgiISJJG/WxI0sTkGg3YR0mwQMOcu4PQRYNooYPCTwOCnAOO5386DPYPRLqSd6ngmXWKubnF1gw2dqs4EsfkgSNphYM37wI4fzyx58W8E9LlHu66oRSYSERHZOO9goP9DQN8HtIzTTZ9rneoO/K5tsX2165BWI7RrEjtn4+/A+tmRssP+lsJIgbKZNwBx/wCunsC4aVoqNZENkEZUJ0/nY/OJ09ZMD5m9r9SlFp5uBnRpFKgVL20sNT0C7WfGnnQl2UCyFEq28tlFkiGyNzGrwpKa46fzVIFW2f46mGp9vKvBRQVCLMtpLEGScD8P21heFdMDuHe1Vshs+0xg1ZvAkT+1ek9BTc75rQOiB6ggyNq4tQyC6Cjbmglio9lrKQeAv97Rlr5Y2tvG9AT6/Z9WW8xoo+MmIqKLZzAArYdrW+pBLQtViqme3ADMGg+EtQUGPKzVDbHj9wMGQRylHkjGCeD764C0g4BnIHDzbKAxi9+RfgqKTdgVp7WptRQwrao1akygV1nAQwt8yAdOu1onTzbNYHBRLUhlG9Eh0np/bmEJDiRJQORMlxrZyyy91rEmG9gWX6HQrmSKtC1bTiNFWFtF+MLTTYc6S1J34dpPtaKpCx8FTm4Epg4CRr+v1WU4R12QL3Z+gfXx62EuNavlntTwbDYTJGk3sPptYPd8rZi6aDVcW07L6wkiIucT2gq4+kNg6L+BjZ8C/3wDpOwD5t+vvV/I/R2vs8uakzb2Dmwbik3F2J26236CIIm7gOljgZxELVV1whwgvK3eoyInk5hZUBbs0IIeu+MzUWyqmObhZpQ2tQGqW0vPplrXlsgAdiuihufj4aqWWMlWPltJ6onIMpoz7XuzcSQlB+l5xVh/JE1tFrIiS4quWgqwWoIkEthrkKyRTtcDjXoCc+4GTm0Cfr4TOLRCK4Lt4XvWw7uEd1HLO9ML07H39F50COlQ/2OkamuC2EwdI1laJYXUd887c1/bq4DBT7DeBxERQdUmu/wVYNDjwD9fA+s+0t475t4F/PUucOmz2vuGLWTM1hCDIFXYd3ofisxFCPAIQBP/c6cX607a1M26GSjMBMLaaQGQgBi9R0UOTtqZ7onPsgY8pJZHXEb+WY8L9fVAjyZaAVMJeHSMCdBn5pyoBiRwIQEM2S5rF1Ehq0kK9JYvwirHEhg5nJKrtkU7EqyPlxn+dmUBEUtgRLoXSeClzgU1Be74XVsWI0sYpKiZVHof+yUQ073CQ90Mbugd2Rt/nPxDdYmxhyDIxx9/jLfffhuJiYno0qULPvzwQ/Tu3bvKx06bNg133HFHhfs8PDxQUFAAW2IzmSC5qdpMnhQ9NUtgxgXoMEbL/Ii0s4LwRERU/zwDtI4yve4GNk7Vlsqk7AV+nKB1CrvsBS1L1Q4wCHKOpTCdQzvbxhrw6kiHAJkBNBUCjfsBN81ioTKqF2k5hWo5iyXTY8epDBQUl60VLzcrLssEupcrYNooqIFmxInqkQTuJIAnW/msEaknYinAaskckWCJfMjddOy02sprEuJdrhCr1qkmNshbLdm5KFIUVWZhmg8F5t4NnD4MfDUcuOx5oN9DFQqYSV0QCYKsjV+LuzvfDVv2448/4rHHHsPUqVPRp08fTJkyBSNGjMD+/fsRHh5e5ff4+/urr1vY4vnHUhPE30unTJCiPC2tec0UoLCsWHDLy4HLXwYibD8wRkREOvPw1bIFe91V1kXmEyBhGzD9Oq1w6hVv2HxrXQZB7LUoqrQvWvSEtm5X0o9k1s/NS+9RkQOQ9qMHk7OtxUsl6HE0NfesxwV4uaF7WZtaqeUhxUzrZaabyAbJh+twf0+1DWkdZr2/qMSMwyk5Zdki2nIaCZIkZxfieFqe2pbsTrI+3sfdiDZlRVjblSvGekFFM5sOAO5bA/z6L2Dvr8CyF4DDK4FrPwP8tHoo/WO0lrrbk7cjpygHvu5nL5uxFe+99x7uvvtua3aHBEMWLVqEr7/+Gs8880y1v5fIyDO1X2yRNROkoc+XpaXakpclzwLZZfVuorpoKc4SQCMiIqoNr0Dg0ueAPvdpy2Kko8zBJdq1hxTTlg52UsfMBvETy7mKoobbYBBELmL++I+Wwip63AGMetcuC9KQ7axPl+Us0qZWAh5yLK1rK5OOGdKtpbta3hKs2tZe9Aw2kYOR9rqqcGqUP67tVjGbylJ8VYIj+5OycCApR7XvlSwr2cqTLCpLtohlSU3TEJ/zt4aWFnc3fA9s+Q5Y/IzWOebT/sA1nwBtrkCsXywa+zXGiewT2JS4CZc2vhS2qKioCJs3b8bkyZOt9xkMBgwbNgzr16+v9vtycnLQpEkTmM1mdO/eHa+//jo6dKg+u6GwsFBtFllZFdsoO8xymNNHgd+e0NoeioDGWqZQx+sdotUhERHpyCdUy/7oead27SHvNWs/0LrKDHsJ6Dze5t5rGASpJCUvBfG58apqfqfQTrApphJg0aPaxa0YOhkY8rRdFaEhfUkKv2R1WLq1SODjQHK2iq2VJ7PTXWLPZHl0jw1CgLeNFPEjskPS5nlAS9lCrfeVmMw4lpZbrs6IljUixVmlra9sy/cmVWgf3TrCr2K9kUg/BPm4V3wyeU/ocZu2THLOnUDiTmDWjUDve9Wsv3SJObH/BNbFr7PZIEhqaipMJhMiIs7UZhFye9++fVV+T5s2bVSWSOfOnZGZmYl33nkH/fv3x+7du9GoUaMqv+eNN97Ayy+/DD0KozZIi9ySIm3NtkyclBQARndg4GPamm43FqUmIqI67iZzy8/AgSXAksla8VTpJLPle+Caj2xqiUydB0HkgmLu3LnqIsXLy0tdgLz55pvq4sSeskBaBrZUVfRthqzhnTMJ2P8bIG0NJftDom1E55BfZML2U1otj60ntMDH6dyisx4XG+yFnk2CtYBH40BVxNGVbWqJ6pX8G2sZ7qe20V2irfdn5hVjb2JWWa2RbOxLysb+xCxVh2fHqUy1lRfp72kNilgyR5qH+cAtrDVw1wpg+cvAho+BTZ+pYtoDBt2PHwCsiVujAqO2WDfjQvTr109tFnL90a5dO3z22Wd49dVXq/weyTSRuiPlM0FiY2MdIxPk+Hpg4SNaO0PRbDAw6n0gtGX9Pi8RETkvFxeVeYoWlwAbPgVWvQWcWAd8OgAY9qI2IWMDWSF1/g68atUqPPDAA+jVqxdKSkrw73//G8OHD8eePXvg42NDQYXzFUUN6wybkXcamDUeOLkRMHoA138NtLtK71GRjZEPM9KhxdKtRfZ7ErJUjY/K6fpdGmltarWgRxDC/Dx0GzcRVSRZV32bh6jNQv4dnzidpxVgTcjCnrIlNSdP5yMxq0Btf+5PsT7eXQVYfFVwpF3kXeh3SXe03/g0DMm70XvBk3BtFIG4nDi1LMYWu6CFhobCaDQiKelMJoyQ2zWt+eHm5oZu3brh0KFD1T5GusfIpkdh1HrLBJHsD1k2K6nIUjfMOxQY8TrQ+QZmjhIRUcNw9QAGPgJ0uBZY8CBwdLW2VEYae1zzse5ZIXUeBFm8ePFZLeukirus7R08eDBsnc0VRc04CUwfC6Tu19oS3fQj0OTMTBc5r8ISE3bFZakMDwl4yCbFFyuTWWLrspbGgegQHaACIURkP6QWSLNQH7Vd2SmqwgfqA0nZZy2pkVojEgSVDYiT6mUIxav4wOszDCjehm75efjbyxPzd/yKhwc+CFvj7u6OHj16YMWKFRgzZoy6T+p8yO0HH6zZeGU5zc6dO3HllVfCltRrJkjqQWDOXVqVftFtAjD8NXaOIyIifQQ1ASYuADZ/Ayx9HjixXssKkXa6fe7Vra5lvdcEkXW5Ijg42GaKklWn2FSM3Wm7bScIkrRHC4BIFXe/aGDCHCCivd6jIp0kZxdYu7VIwGPnqUwUmSq2qXU1uKBDtD+6lWtTGx3IrkFEjkqyCaRQsWwWZrOWFWZp2yv7/YnZOJoGTMh/AncYl6BP3q8qCLJz16dAv3sBo+3V/JFlKrfddht69uyJ3r17qxa5ubm51m4xEydORExMjFqGK1555RX07dsXLVu2REZGBt5++20cP34cd911F2yJpSaIf11mgkhhpy3fAosnA8V5WtDj6g+BdqPr7jmIiIguhGQhShmHFpcBCx4Cjq7SaoZknABG/hcOFwSRWZtHHnkEAwYMQMeOHW2mKFl19qfvR6GpEAEeAWjq31TfwRxfpy2BKcgEQttoAZDA+l2nTLZDCibKhxcJeEjx0s0n0lXae2XBPu4qu0OyPKRzS+dGgfByZ6cgImcmXZtig73VNrzDmaUjeUUlOJgk7Xu74MDBdojO+xRenh1sMgAibrzxRqSkpOCFF15AYmIiunbtqrJNLcVST5w4oTrGWKSnp6uWuvLYoKAglUmybt06tG9vO5MHxSazqu1Sp5kgsmRWLir3LdRuNxsCXDsV8D9TZ4aIiMg2skJ+ATZPA/54Heh9t25DcSmVQgL15P7778fvv/+ONWvWVFuZvapMEClKJhkk/v7+aEgz9s7Afzf9F4NiBuGTYZ9AN3sXAj/fCZgKgdg+wE0/aG0PyWFl5BVZ63jIJsVM84pMZ9cZivCzBjxk3zTE22GKGhJRw1EFUaVbiKtnndeJkPfxgIAAXd7HbX3sUpi6+6vL1PHB/4yE28UWoE7eB8y8Acg4DhjctLa3/R6yiaJzRERE1SouqJcuZTV9H6+3TBBZs7tw4UKsXr262gCIXkXJqrM9ebv+S2H++RpY9DhQagbaXKkVQXXjcgZHIqnqh1NyrMtaZDucknvW42SWUC1rUQGPQHSNDWyYlopE5PBU8JTvLQ3OUhTVy8148QGQQ8uBn+4ACrOAwCbADd8B0V3rZqBERET1Sec27a71Mbv00EMPYd68efjzzz/RrFkz2IsdqTv06wwjCTl//hdYVbYuqvtErZWdsd7LtlA9yykswfaTWpaHZXlLVllhvPKkpaXq2NI4CD2bBqFlmK9KayciIsdgKYrq73WR7+2bvgB+fxooNQGN+wM3Tgd8znQTIiIiourV+SdsaY87c+ZM/PLLL/Dz81Nrc4WkpXh52e6sU2p+qmoX6AIXdArt1LBPbjZp2R9SNVcMeRoYOpmt7OyQBAGljaUl4LH5eAb2J2ahUpdaNQvYJVZrUyvFSyXjQ+p7EBGR47IURb3grD5TCbDk38Cmz7TbXW4GRk/RWhESERGRPkGQTz/9VO2HDh1a4f5vvvkGt99+O2yVZSlMy6CW8HX3bbgnLs7X2tmpgmYuwKh3gF62VcmeqldQbMLOuEzrshZpV5uaU3TW42ICvazdWiTw0TbK7+JToYmIyK5cVHvcwmzgp9u1ZTDisheBgY9ywoSIiMgWlsPYo+0pOtQDyU8HZt2k9Us2egBjvwTaX91wz0+1lpCZr9rUqqDHiXTsic9Esani33l3owEdYvxVLQ8V9GgShAh/fde9ERGRLQVBapkJUpQHzLwROL4WcPUCrvuc1wtEREQXiAUn9AqCZMYB08cCKXsBjwDgpllA0wEN89xU41aGe+KzKtTyiM8sOOtxYX4e1uKlkuXRMSYAnm5sU0tERFUXRq1VJohU0P/hZi0A4uEP3DoPaNSz/gZJRETk4BgEkesLczH2pO1puCBIyn7g++uArFOAXxQwYQ4Q0aH+n5fOKTWnUAU6tpzIUHtpU1tYYq7wGKPBBe2i/KwFTCXTo1GQF9vUEhFRzQuj1jQIYirWlsAc+QNw8wFu+YkBECIioovEIAiAA6cPoMBUAH93fzTxb1K/T3ZiIzDrRm0pTEgr4Na5QGDj+n1OOovJXIoDSdlalkdZpsextLyzHhfo7YZusYHo2TQY3RoHokujQPh48J8NERHVXlZ+LQqjShFUqRl24HfA1VPLGG3ct/4HSURE5OD4aQ7AtpRt1ta4Bpd6LFa5/3fgpzuAknwgpidw82y2tGsgmfnF2GZpU3s8XR1L69rKWob7omdZHQ/J8mge6sMsDyIiqtuaIOcLppvNwIIHgT3zAYOb1gK3+ZCGGSQREZGDYxCkoeqBbP4WWPgIUGoGWo0Axn0DuPvU3/M5MSnOeyQ119qtRfYHk3NQuWavj7tRtabt3jhQBT26xQYhwPsC2xYSERGdR3ZhDWqCyJvVb08A22cBLkbteqHV5Q03SCIiIgfHIAiAHSk76i8IIhczq98B/nhNu911AjB6CmDkh+26kldUgu0nM63FS2WfnqddaJbXJMS7rICpVs+jTaSfqvFBRERkM91hdvwI/PMVABetC0y70Q03QCIiIifg9EGQ1PxUxOXEwQUu6BTaqW5/uNkE/P4U8PeX2u1BTwCXPgdwecVFZXmcSs+3BjykTe3ehGxV46M8d1cDujQKsBYvlcBHqK+HbuMmIiLKshRG9aomCJJ+DFj0hHZ8yb+BTtc34OiIiIicg9MHQSxLYVoEtoCvu2/d/WBpaTf3bmDvAm02Z+RbQJ976u7nO4nCEhN2xWVpAY+yLI/k7MKzHhcV4Kl1bCmr5dE+yl8FQoiIiOyiRa4UQp17D1CUDcT2BQY93vADJHIAJrMJJaUlKDYVw1RqUl0gS8wlapPb5Y8tjy1/u/xevibH5lJzxa+Xlt1X7mulKNVuy3+lFTeZxFPH0I4t3y/U1+W/ssdYjuU/YTmucL/6X6V13jUgk75lB9bbUvvO+l/ZRK3lWGolln+M5bbs1bE8BuWOy+5XW9n9RoNRfY/RxWi9rfYuxgr3WW5XdexqcD2zNxjh6nLmdoXNpeJteQxr+1FVGASpj3og+RnADzcDx9cCRnctnbXDtXX38x1YclaBCnRIwEM2CYAUmSq2qXU1uKBDjGR5BGpZHo2DEB3opduYiYiIarccporLrzXvASc3Ah7+2nWDwdjwAySqBflALgGGQlNhha3IVFRhL5sEJNR95iJ1v3UzF6mvWe6Xn2fZy/1qb9kq3y53nwQxLPsLCQ6Q45JgiJvBzbqvcGzUblfYjG5wN7hbv+ZudNf2Bnft2Hjm2MPoob4me3VsPHNs+bplX35jYEZ/DIIk13EQJCsemH49kLxbu5AZPwNoNrhufraDKTGZsS8xu0LQQ5a6VBbi426t4yFBj86NAuDpxotDIiKyz0wQ/8o1QU7+Dfz5X+141LtAUBMdRkeOFJyQoEJ+cT7yS85seSV5KCgpUMcFpoIzxyUFFW5LsKL8fXLbcn9hSSEKzWV7U6HdBBzKZwhYMwkMRvUB1pJxYMkcsN52cT07S6FSNkPljAbLbfmQa7nPclw+a0JYsynKZ2OUZVNYvl5Vhobal3v8hfz9UPsaZJyU31uyWOS/CtktZdkrlbNfzsqQKbtdeV8h0+Yc2Tjls3gqZ/JYAmBVsTzelkgQxcP1TFDEy9XrrGNPV88Kx+q20evMsauXdfM0lt12q3gfgy3Vc+ogiPyD2ZO2Rx13Ca+DIEjKAWD6dUDmScA3ApgwB4is4zojdiw9twhbT54JeEgx0/xiU4XHyL/VNhF+1gwP2UtBU/4jJiIie1ZUYkZBsfnsTJDCbGDuXUCpCeh4PdD5Bv0GSbqQD2i5xbnIK85T+5ziHO24JFfdlk2CEOUfIwENFdgoLtuX5FU4tiy1aCjygVx9iHP10D7gVTETbplBt9xnmWUvP8NumXEvPxuv9mUz9BKYKD97X3l237Kp4EZZYEO+h9eRzrcUSgU/ypY6WTKHymcMScZR+a9VyDaqJjvJkr1UOZup8nHlrCjZyv+bVI8rKkI2suv136QES7xdveHt5q0CI5Zj70p7Hzcfdaz2lW7L5uvmCx93H/XvylE4dRDkQPoBFeX2d/dHU/+mF/fDZBZn5jggPx0IaQlMmOvUMzlmcykOpeRYa3lIAdMjKblnPU4uBFUtj7KAR5fYgHNXzSciIrLjLBDh61Hu8uv3Z7SCqAGxWhYI2RX5sJRVlIXsomxtK85GTlGOOpZghmUv95XfW4IdlgBHfZEAQuVZYzVLXHafmlU2njm2zETLfRLQsHzdMhut9mVfs8xcy20JNjDQQHpTWTswqr+XtsYSRJHsKrU3aVlWlmyr8hlXKhOrXDZW5UytClldZXvLJl8XkoFjuS+tIK1OXoOH0aNCYETqacqxn5uftnf3U/fJ1+Tzdflj+Zpsch6xBU4dBLEshekU1smaenZBDiwBZt8GyJtYTA/g5tmATyic7eJOMjssxUu3nki3VsEvr3mYj2pTKwEP2VqE+cLANrVEROQk9UC83Y1wNZZdc+z5Bdg2XatSeO1ngFegvoN0UjIbLIGMzMJMtVU+ttxWx4VawMMS+JAPJ3VFZlktHzAsm5qVda04O1vVTK7sLQEOy6yvbBKcICL9WTKW5N9xfWfDyHlJBUCK8ytkjamMMbldnK8yzaz3yb5c9pkcW4K0cmw5z1mCNacLTl9UYNYSELm3y724qvlV0INTnxnrpCjq1unAgn9paawtLwdu+BZwr9+/3HqTNYDH0/Iq1PI4kJSNSl1q4eVmRNfYQHRvohUw7RYbhCAfd72GTUREZFtFUddM0fYDHwGaDtBpZI4X0MgozEB6QXqFvVy0y96yZRZkIr0wXQU35GL/YllmRdXFvZt2gW+ZBZVj62yp+5lZU8ssqko1d/NRS0KIiC42G8bHoAVR4VW3S/ZyyzZLVpscS/ZbbpG2t2S6lc+Ak4Cx5T5ZEiRLgSQzRTbJiNELgyAXGgSRokJ/vQusfFW73eVm4Or/AUbHW8pRUGzCjlNalodskuWRlnv2X9pGQV7WWh6ytYvyOzPbRURE5MTOtMctu07ITQXit2rHfe7TcWS2P/EiF9qp+akqkCEXzqfzy/YFp7X78tNUQEOOJahxoSQ44e/hj0CPQJW+HeARoPZyX4B7gNqXT+u2HEsQQz54EBE5IleDqzofynahJAAiWSXlM+ka+zeGXpw2CCJvpnE5capoTOfQzrX7ZrMJWPwMsOlz7faAR4BhL2lVPR1AQma+NeAhNT12x2ehpFKah7vRgE6NtDa1lnoe4f62scaLiIjI1mRVzgQ5/IdatY2IToBfJJyNXBBLlkZKfgpS8lLUPjkvWV2fld8kwFHbJSdybScX6xLMCPYMVvsgzyC1l02+ZrlteZwEM7h0hIiofkjpCZUd5+6LKERBb67OngXSIrCF+mXUWEkhMPceYM98bQ3vFW8Afe+HPVer35OQpRUwPaEFPRIyz77YCPPzsNbykHa1HWP84eHKWQ8iIqILygQ5tFzbt7wMjkZSpyWwkZiXiMTcRCTlJiEpT9sk0CGbfF06N9SU1LwI8QpRQY0Qz5Azx14hKqAR7BGsbluCG8zMICKi6jhtEGRHyo7aL4UpyAR+uAU49hcgLYKu+wzoOBb2JDWn0Brw2Ho8A9tPZaCwpGIbNaPBRS1lkaCHBDwk00OWurDqNxER0cXVBPGXTBCzGTi8QvtCy2GwN1JkLyEnQWXUJuSW7XMSEJ8br4IektVR0xatErgI9w5HqFeodR/mFab2skmQQ4IeUuyTiIioLjhtEKTW9UCyE4Hp1wNJOwF3P2D8dKD5UNgyk7kU+xOzywIeWuBDCppWFujtViHgIW1qvd2d9q8GERFRPRZGdQMSdwC5KYBkosb2gS3W4ZClKMezjuNUzimcyj51Zp99qkbtFl1dXBHhE4FIn0hEeEeoTYIccp/ae0eoAId0SyAiImpITvlJV/o0707dXfMgSOohYPq1QMYJwCccmPAzEHURHWXqSWZ+sSpaKpkeW05kqOPcIlOFx0gyR6tw3zMFTJsEoXmoD7M8iIiIGmA5jMoEObREu7PZEMBVv44gUpjuaOZR63Yi+wROZJ1Qe8n2OBcpBhrtG41on2hE+UYhxjdGBTyifKLUXrI3uCSFiIhskVMGQQ6kH1BFtqQIVtOApud+8KnNwMxxQF4aENwcmDAXCG4GvZnNpTiSmlsW8NCKmB5MPrvFm6+Hq9amtnEgejQNVscBXpx1ISIi0q1F7uGV2p0tL22Q55aOKYcyDuFQ+iEczDhoDXrIspVzFbGTAEesXywa+TXSNl9tLwGPi+kSQEREpCenDIJsT9aWwnQO66ze5Kt1cBkweyJQnAdEdwNu/gnwDWu4gZaTW1ii6ndYsjwk8JGRp80qldc0xFtld1gyPVpH+KkaH0RERKSfrLJMkGBjAXByo3Zni7otiip1OGS5yp7Te7AnbY+a9DmYflAVIq1OuFc4mgU0U5NCjf0ao4l/E9W2UAIebkZOmhARkeNxdep6IKHnWNKybRaw4EHAXKJdpNzwHeBRiy4yF7kW91R6vjXDQ/Z7E7JVjY/yPFwN6NIo0Br06NY4EKG+Hg0yRiIiIqp9JkjznM3atUVIy4vOLJWuK9tStmFnyk4V+NiXtg/ZxdlVPlaWqbQMbIlWQa1UZ7xm/lrgQ7JiiYiInIlTBkHO2RmmtBRY+wGw/EXtducbgas/qtc1uwXFJuyOz9QCHsczVAHTlOzCsx4XHeCJbhLwKKvl0T7KH+6u58hkISIiIpuqCdIobd0FdYUxmU3Yn74fW5O3qoxWCX5IZ5bK3A3uaBPcBu2C26m9BD0k+MFgBxERkZMGQdLy01SFcxe4oFNYp4pflJZ1S58FNnyi3e73IHD5q4ChbgMNSVkFWpvaso4tu+OyUGSq2ErOzeiC9tEBZQEPqekRhOhArzodBxERETVkJkgpwpL+qnEQRNrNro9fj3Xx67AhYQMyCjMqfN3oYkTroNZqUqdDaAcV+Gge2JwdV4iIiM7B1VmXwkgqaIVZkZJCYN59wO652u3h/wH6P3jRz1dsMmNfQrZ1aYtscRlnV1wP9XVHt8ZB6Cm1PJoEoVNMADzdWFWdiIjIEWQVlKCFSzw8cuMBowfQZECVj9t/ej8WHVmE1adW43Dm4bM6snQN74pu4d3QNawrOoZ2hLebdwO9AiIiIsfgtEGQCkthCrKAHycAR1cBMnsy5lOg87gL+vmnc4tUa1pLLY/tJzORX1yxTa3UKZWCpT2bnilg2jjYm21qiYiIHHg5zDUG7RoETQcA7meCF1K49Lcjv+HXI7+qYqYWUrxdAh39o/urrVNoJ7ganO7SjYiIqE45bRBEOsMo2UnAjOuBxB2Auy9w4/dAi0tr3KZW2tJaC5geT1dtayvz93RVWR4S8JCtS2ygal1LREREjq+oxIzCEjOGuG2vsBRGrkk+2faJWuoinV2ELGUZ0mgIRjQbgX5R/diKloiIqI451SfxYnMxdqfuPpMJknYY+P5aIOM44B0KTPhZa4V7jlmcbSczrMVLJePDUu29vBZhPtYMD1na0jLMFwa2qSUiInJKcv3giUL0MezT7mg5DMuOL8Mzq59BkblI3SXLW0a3GI0RTUcw8EFERFSPnCoIcjD9IApMBaoWSLOcdGDmjUBeKhDUFJgwFwhpUaFN7fG0PGvxUsny2J+UrZrHlOftbkTXWK1wqaVNbaB3/XWSISIiIvsiEyZ9DXvh4VIMBMRiVuoWvLHpDZSiFENjh+LJnk+isX9jvYdJRETkFOotCPLxxx/j7bffRmJiIrp06YIPP/wQvXv3hk0shfGOgeHbq4HiXCCqC3DLz8h3D8GOI2llAY8MtcRF6ntUFhvsZQ14yL5tpB9cjWxTS0REpPe1xE8//YTnn38ex44dQ6tWrfDmm2/iyiuvhN6yCoox2LADMo/yYXRTfLHpdXX/uNbj8GyfZ2E0sBA6ERGRXQdBfvzxRzz22GOYOnUq+vTpgylTpmDEiBHYv38/wsPDoXsQ5NgmFQBJCeuHL8JfwYZpB7AnPgsl5oppHu6uBtWl5czSlkCE+3nqNHoiIiLnUdtriXXr1uGmm27CG2+8gauuugozZ87EmDFjsGXLFnTs2BF6Z4IMMGzH86HB+KXguLrvga4P4N7O97IoOhERUQNzKZV1H3VMLlZ69eqFjz76SN02m82IjY3FQw89hGeeeeac35uVlYWAgABkZmbC39+/zsaUkVeE634ajBTkYmpiMlJzuuOJ4vtQXC4OFO7nYS1eKrU8OkT7w8OVszNEREQ1VVfv47W9lrjxxhuRm5uLhQsXWu/r27cvunbtqgIpDTn2ypauX4N5227HGm8vGF2MeL7v8xjbemyd/XwiIiJCjd/H6zwTpKioCJs3b8bkyZOt9xkMBgwbNgzr168/6/GFhYVqKz/w+nD60GIVABHbcgfjfdNtaBcjtTwCVcBDAh8xgV6ckSEiItJZba8lhNwvmSPlSebI/Pnzq32ehroGWXr0QxUAcTcD7w37AENih9TL8xAREdH51XkQJDU1FSaTCRERERXul9v79pVVRS9H0lZffvll1LfA1gMxakNjZHsY0f22qdgRGwBvd6eqC0tERGQXanstIaRuSFWPl/ur01DXIDe3ugsnd76Cq30GMgBCRESkM92jADLLU37mRmZhJN21rgV7BuO/9yyq859LRERE9qmhrkG697ocP/a6vM5/LhEREdlAECQ0NBRGoxFJSUkV7pfbkZGRZz3ew8NDbUREREQXci0h5P7aPF7wGoSIiMj51HlvV3d3d/To0QMrVqyw3ifFzOR2v3796vrpiIiIyMFcyLWE3F/+8WLZsmW89iAiIqL6Xw4jqaW33XYbevbsid69e6u2dlKx/Y477qiPpyMiIiIHc75riYkTJyImJkbV9RAPP/wwhgwZgnfffRejRo3CDz/8gH/++Qeff/65zq+EiIiIHD4IIm3qUlJS8MILL6iCZNKebvHixWcVLCMiIiK6kGuJEydOqI4xFv3798fMmTPx3HPP4d///jdatWqlOsN07NhRx1dBREREtsaltLS0FHbY25eIiIhsjz2/j9vz2ImIiJxdVg3fx+u8JggRERERERERkS1iEISIiIiIiIiInEK91AS5GJbVOZLKQkRERPbF8v5tY6tta4TXIERERI5/DWJzQZDs7Gy1j42N1XsoREREdBHv57Iu157wGoSIiMjxr0FsrjCq2WxGfHw8/Pz84OLiUueRIbmwOXnypMMWPHP01+jor88ZXqOjvz7B12j/HP311edrlMsKufiIjo6u0L3FHvAa5MI5+utzhtfo6K9P8DXaP0d/fc7wGrPq8fXV9BrE5jJBZLCNGjWq1+eQP2xH/AvlTK/R0V+fM7xGR399gq/R/jn666uv12hvGSAWvAa5eI7++pzhNTr66xN8jfbP0V+fM7xG/3p6fTW5BrGvKRoiIiIiIiIiogvEIAgREREREREROQWnCoJ4eHjgxRdfVHtH5eiv0dFfnzO8Rkd/fYKv0f45+utzltdoSxz9z9vRX58zvEZHf32Cr9H+Ofrrc4bX6GEDr8/mCqMSEREREREREdUHp8oEISIiIiIiIiLnxSAIERERERERETkFBkGIiIiIiIiIyCkwCEJERERERERETsHhgiAff/wxmjZtCk9PT/Tp0webNm065+N/+ukntG3bVj2+U6dO+O2332Cr3njjDfTq1Qt+fn4IDw/HmDFjsH///nN+z7Rp0+Di4lJhk9dqi1566aWzxiq/G0f5/Qn5u1n5Ncr2wAMP2O3vb/Xq1Rg9ejSio6PV+ObPn1/h61J7+YUXXkBUVBS8vLwwbNgwHDx4sM7/Levx+oqLi/H000+rv3s+Pj7qMRMnTkR8fHyd/13X83d4++23nzXeK664wiF+h6Kqf5Oyvf3223bxO6zJe0NBQYE6z4SEhMDX1xdjx45FUlLSOX/uhf7bdWaOeg3i6Ncfgtcg9vc7dPTrD2e4BnH06w/BaxDY5DWIQwVBfvzxRzz22GOq5c6WLVvQpUsXjBgxAsnJyVU+ft26dbjpppswadIkbN26Vf3SZNu1axds0apVq9RfoA0bNmDZsmXq5Dd8+HDk5uae8/v8/f2R8P/t3QtQVNUDBvDjA3yWSBjgAzITUkhTS8NHmg9KHSFpQNMh0tQkLWE0s5nMjJkeZlrRDFkZ6tjgUIFMUloolBq+EAsfQ8kg6AQZpaZikXD+853+uy6wyy64wN17v9/Mjtzdu/eee8/evZ/3nHu2rMz8KCkpEVoVFBRUq6z79u2zOa+r1R8cPny41vahHiEyMtJl6w+fPxxrOOFYs2bNGvHee++JDz74QBw8eFCdqHFc4gvRWcdya21fZWWlKt/KlSvVv2lpaeqLPywszKmf9dauQ0DosCxvSkpKg8t0lToEy+3C45NPPlGBAidpV6hDR84N8fHx4ssvv1T/acP8CMkRERENLrcpx66R6TmDGCF/ADOIa9Wh3vOHETKI3vMHMIMIbWYQqSPDhw+XixYtMk9XV1fLnj17ytdff93q/FFRUXLq1Km1nhsxYoR8+umnpSs4f/48ft5YfvfddzbnSU5Olt26dZOuYNWqVXLw4MEOz+/q9QdLliyR/fr1kzU1NS5ff4DPY3p6unka2+Xj4yPfeust83MXL16UHTp0kCkpKU47lltr+6w5dOiQmq+kpMRpn/XW3saYmBgZHh7eqOW4ch1iW8ePH9/gPFquw7rnBhxzbm5u8rPPPjPPc+rUKTVPbm6u1WU09dg1MiNlEL3lD2AGce061Hv+MEIG0Xv+AGYQ7WQQ3fQEqaqqEnl5eaqrjEnbtm3VdG5urtX34HnL+QFXmGzNrzWXLl1S/3p6ejY435UrV4S/v7/o06ePCA8PFydOnBBahW5O6C525513itmzZ4vS0lKb87p6/eEzu3XrVjF37lx1xVcP9VdXcXGxKC8vr1VP3bp1U10TbdVTU45lrR2XqE8PDw+nfda1ICcnR3VzDAwMFLGxseKPP/6wOa8r1yG6Z2ZmZqrWXXu0Wod1zw2oC7TMWNYHus36+fnZrI+mHLtGZrQMosf8Acwgrl+HRs4fes0gRskfwAzSchlENxdBKioqRHV1tfD29q71PKaxE63B842ZX0tqampEXFycGDVqlAgODrY5H74w0K0qIyNDnezwvpEjR4pz584JrcEHG/ef7ty5UyQlJakDYMyYMeLy5cu6qz/APYEXL15U9zvqof6sMdVFY+qpKceyVqCLHu7PRRdpdCF21me9taEr6pYtW8Tu3bvFm2++qboyTp48WdWT3upw8+bN6r5We900tVqH1s4N2Ofu7u71QrG986NpHkffY2RGyiB6zB/ADOL6dWjk/KHXDGKk/AHMIC2XQdo7ZSnU4nDvFe47tXf/V0hIiHqY4OQ1YMAAsWHDBpGQkCC0BF9qJoMGDVIHOFofUlNTHboi6mo2btyothlXcfVQf0aHq9xRUVFqICeckPT0WZ85c6b5bwzAhjL369dPtc5MmDBB6AkCP1pU7A3+p9U6dPTcQNRUeswfWj6mmwsziL7oNYMYKX8AM0jL0U1PEC8vL9GuXbt6I81i2sfHx+p78Hxj5teKxYsXix07dojs7GzRu3fvRr3Xzc1NDBkyRJw+fVpoHa4YBgQE2Cyrq9YfYGCxrKwsMW/ePN3WH5jqojH11JRjWSvhA/WKQaEaaoFpymdda9D1EvVkq7yuWIewd+9eNahcY49LrdShrXMD9jm6CKPVtzHnR9M8jr7HyIySQYySP4AZxLXr0Cj5w2gZRK/5A5hBWjaD6OYiCLrZDBs2THWXsuySg2nLq9iW8Lzl/IAvD1vztzZc3cUHLD09XezZs0f07du30ctAF7GCggL1c0Nah/tQi4qKbJbV1erPUnJysrq/cerUqbqtP8BnFF9WlvX0119/qVGebdVTU45lLYQP3JuJUImf/3L2Z11r0BUa9+TaKq+r1aFlyyjKjVHcXakO7Z0bsE34z4tlfSBo4f5hW/XRlGPXyPSeQYyWP4AZxLXr0Aj5w4gZRK/5A5hBWjiDSB3Ztm2bGjV206ZN8uTJk3LBggXSw8NDlpeXq9ejo6PlihUrzPPv379ftm/fXq5du1aNUouRdjF6bUFBgdSi2NhYNUp3Tk6OLCsrMz8qKyvN89TdxtWrV8tdu3bJoqIimZeXJ2fOnCk7duwoT5w4IbVm6dKlatuKi4tV3UycOFF6eXmpUYb1UH+Wo1T7+fnJF154od5rrlh/ly9flvn5+eqBr5R169apv00jk7/xxhvqOMzIyJA//fSTGvW6b9++8tq1a+ZlYBTsxMREh49lrWxfVVWVDAsLk71795bHjh2rdVz+888/NrfP3mddS9uI15YtW6ZG8EZ5s7Ky5NChQ2X//v3l33//7fJ1aHLp0iXZuXNnmZSUZHUZWq5DR84NCxcuVN87e/bskUeOHJEhISHqYSkwMFCmpaWZpx05dskYGUTv+QOYQVyvDvWeP4yQQfSeP+xtowkziGzxDKKriyCADwh2sru7u/qJpAMHDphfGzt2rPqpJUupqakyICBAzR8UFCQzMzOlVuHAsfbAT5jZ2sa4uDjz/vD29pZTpkyRR48elVo0Y8YM6evrq8raq1cvNX369Gnd1J8JAgXqrbCwsN5rrlh/2dnZVj+Xpu3Az1ytXLlSlR8npQkTJtTbdn9/fxUgHT2WtbJ9OPnYOi7xPlvbZ++zrqVtxEksNDRU9ujRQwV8bMv8+fPrhQlXrUOTDRs2yE6dOqmfYLNGy3XoyLkBoeGZZ56R3bt3V0Fr+vTpKqTUXY7lexw5dskYGUTv+QOYQVyvDvWeP4yQQfSeP4AZRGoyg7T5/0qJiIiIiIiIiHRNN2OCEBERERERERE1hBdBiIiIiIiIiMgQeBGEiIiIiIiIiAyBF0GIiIiIiIiIyBB4EYSIiIiIiIiIDIEXQYiIiIiIiIjIEHgRhIiIiIiIiIgMgRdBiIiIiIiIiMgQeBGEqAU8+eST4tFHH21wnjvuuEO888475uny8nIxadIk0aVLF+Hh4eHU8lRWVorHHntM3HrrraJNmzbi4sWLTVpO3TJjWdu3bxfNZdy4cSIuLs7m+lvCypUrxYIFC5pt+RUVFeL2228X586da7Z1EBGRcTCDOAczCJF+tG/tAhDRfw4fPqzChsn69etFWVmZOHbsmOjWrZvIyckRDz30kLhw4cJNB5LNmzeLvXv3ih9++EF4eXmp5buCtLQ04ebm1mrrRyh89913RUFBQbOtA/XxxBNPiFWrVomNGzc223qIiIhMmEHsYwYh0g/2BCHSiB49eojOnTubp4uKisSwYcNE//791VV5Z8KyBwwYIIKDg4WPj49qPXEFnp6e4pZbbmm19X/88cdi5MiRwt/fv1nXM2fOHPHpp5+KP//8s1nXQ0REBMwg9jGDEOkHL4IQOcnnn38u7rnnHtGpUydx2223iYkTJ4qrV6/Wmmft2rXC19dXvb5o0SLx77//Wu1Wib+/+OILsWXLFhUO0JUVLTDQvXt383O24L1BQUGiQ4cOallvv/12re6cmP7+++/VcjBtK6SEh4cLb29v0bVrV3H//feLrKysZttHpu66q1evVmEM3WQXLlwoqqqqbHZFtRYQ0EK1e/duNX38+HExefJkVX5sR3R0tOrq6Uh5rNm2bZuYNm1aredQpmeffVaVC3WD9Xz00UdqOQgSCEx33XWX+Prrr83vQUva7Nmz1XZi3QiZycnJ5tdRdz179hTp6emN3sdERGQ8zCA3t4+YQZhByFh4EYTICdBl9PHHHxdz584Vp06dUt1GIyIihJTSPE92drY6qeNfdAXdtGmTetjqlvrII4+IqKgotWx0f0SogMLCQvNz1uTl5an3zZw5U3WZfOWVV9Q9pKZ1oTvn/PnzRUhIiFoOpq25cuWKmDJlijqZ5+fnq/Lg5FtaWtps+wjrMr2WkpKiyoZA4og1a9aIFStWiG+++UZMmDBB3WM8fvx4MWTIEHHkyBGxc+dO8dtvv6l942h5LKFF5OTJk+K+++6r9xrqE11IDx06pMJIbGysiIyMVC02R48eFaGhoSr84D5oQH1gWQglWHdSUpJ6v6Xhw4er7sJEREQNYQZxzj5iBrmBGYR0TxLRTcvLy8NZS545c8bq6zExMdLf319ev37d/FxkZKScMWOGeRqvr1+/3jwdHh6u3meSnZ2t1nHhwoUGyzJr1iw5adKkWs89//zzcuDAgebpJUuWyLFjxzZyK6UMCgqSiYmJNsuM8qWnpzd5H3l6esqrV6+an0tKSpJdu3aV1dXVahplRtnrrn/58uXS19dXHj9+3PxaQkKCDA0NrbWOs2fPqjIUFhbaLU9d+fn5av7S0tJaz6NMo0ePNk+jjrt06SKjo6PNz5WVlan35ubmqulp06bJOXPmNLi++Ph4OW7cOIfKRkRExsUM8h9mEGYQIkexJwiREwwePFhd+Ue3Rlx9R1dEdDe0hO6F7dq1M0+jS+r58+edXhZc1R81alSt5zD9yy+/iOrqaoeXg1aYZcuWqft20b0T3Tmx7Ka2wjiyjzCP5T3JaClCOc6ePWtzuehWi2Xt27dP7WOTH3/8UbV4odymx913361eQ2uYI+WxdO3aNfVvx44d6702aNAg89+oY3RrxXJN0D0VTPWNVhp0a7333nvF8uXL1eBwdaGLqqnVhoiIyBZmEPuYQZhBiCzxIgiRE+Ck8+2336quhQMHDhSJiYkiMDBQFBcXm+epO6I47oWtqakRWoXwgftBX3vtNdUlEiPE46RqeX9sYziyj5pizJgxKlilpqbWeh7BBV1nUW7LB4LYgw8+2OjymLqKWgsp1urW8jnToG+m+sY9wiUlJSI+Pl78+uuvKghhf9ft+or7dYmIiBrCDGIfMwgzCJElXgQhchKcZNDagftHcf+qu7u7UweVwvLAXksKWk32799f6zlMBwQE1GoFsgfvwUBh06dPV8EDI7ifOXNGNOc+QsuJqbUDDhw4oFpP+vTpY3OZuG8VIQJBCYO+mQwdOlScOHFCDcqGQcEsH6afAWxMnfXr108NlIb7aJ0B4SImJkZs3bpVDUb34Ycf1nodA6rhXmIiIiJ7mEHsYwa5gRmEjI4XQYic4ODBg+oEiMGv0FUTg2n9/vvvKgw4C34SDSfMHTt2qGWjlcGapUuXqsG9EhISxM8//6wGzHr//ffrXeW3B6OFYzvQcoFgMGvWrJtqNXJkH6GF56mnnlIn+a+++kr9Tv3ixYtF27YNf1Vh8C/MjyBhGt0eI9+jJQMDj2GQN3Q/3bVrlxotHSGusXWGMmDkdnR5vVkvv/yyyMjIEKdPn1YhCXVquV50QcXgchjMjIiIqCHMIPYxg9zADELEiyBEToGr8/i5N4xkjtaOl156Sd0nii6HztKrVy91gsXo47i/Eydma9D6gG6ZuN8zODhYnexeffXVBn/Ozpp169apn1vDyR1dOh9++GG17ObcR+iSieCDrqIzZswQYWFhamR5R4wePVpkZmaq5aJbKX7eDS1JCBs4kaMlCT8hh3uLESaaUmfz5s1T+/VmuxCjtefFF19U9/GausViuSYIJ35+fqqbLRERUUOYQexjBrmBGYRIiDYYHbW1C0FEhICEn5Tbvn270Cp8XY4YMULdR4vWnebywAMPiOeee061fBEREVHzYga5gRmEjIA9QYiIHISuwLhv9vr16822joqKChEREdGsAYeIiIhcCzMIkfOwJwgRaYIrtMIQERGR/jCDEBkLL4IQERERERERkSHwdhgiIiIiIiIiMgReBCEiIiIiIiIiQ+BFECIiIiIiIiIyBF4EISIiIiIiIiJD4EUQIiIiIiIiIjIEXgQhIiIiIiIiIkPgRRAiIiIiIiIiMgReBCEiIiIiIiIiYQT/A09DPYXzjqSyAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 1100x320 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "shifts = np.linspace(0, 20, 81)\n",
    "fig, axes = plt.subplots(1, 2, figsize=(11, 3.2))\n",
    "for q in [0.05, 0.2, 1.0]:\n",
    "    axes[0].plot(shifts, [float(ba.spike.victor_purpura_distance(a, a + s, cost_factor=q)) for s in shifts], label=f\"q = {q}/ms\")\n",
    "for tau in [1.0, 5.0, 20.0]:\n",
    "    axes[1].plot(shifts, [float(ba.spike.van_rossum_distance(a, a + s, tau=tau)) for s in shifts], label=f\"tau = {tau} ms\")\n",
    "axes[0].set_title(\"Victor-Purpura\")\n",
    "axes[1].set_title(\"van Rossum\")\n",
    "for ax in axes:\n",
    "    ax.set_xlabel(\"shift of all spikes (ms)\")\n",
    "    ax.legend(fontsize=8)\n",
    "fig.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "04c9669b",
   "metadata": {},
   "source": [
    "## 4. Coincidence Measures and Firing Rate\n",
    "\n",
    "Three functions summarize synchrony from a spike matrix:\n",
    "\n",
    "- `spike_train_synchrony`: for each pair, the fraction of spikes that have a partner in the other train within `window_size / 2`, averaged over pairs. Between 0 and 1.\n",
    "- `correlation_index`: the mean Pearson correlation of spike counts in bins of `window_size`.\n",
    "- `spike_time_tiling_coefficient` (STTC, Cutts & Eglen 2014): compares the fraction of A's spikes within $\\pm\\tau$ of a B spike with the fraction of time covered by $\\pm\\tau$ windows around B's spikes. Matrix output, 1 on the diagonal. The default is $\\tau = 5$ ms.\n",
    "\n",
    "The crucial difference shows up with **independent** trains at increasing rates. A coincidence fraction grows simply because there are more spikes to coincide with, whereas STTC and the correlation index subtract what chance predicts:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "2f307d21",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-10-15T08:09:49.627643Z",
     "iopub.status.busy": "2026-10-15T08:09:49.627643Z",
     "iopub.status.idle": "2026-10-15T08:09:51.769895Z",
     "shell.execute_reply": "2026-10-15T08:09:51.769895Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  rate  synchrony  corr. index    STTC   independent trains, window 10 ms\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   2Hz      0.023        0.001   0.001\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  10Hz      0.104        0.001   0.001\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  50Hz      0.435        0.005   0.004\n"
     ]
    }
   ],
   "source": [
    "T = 20_000  # 20 s at dt = 1 ms\n",
    "print(f\"{'rate':>6s} {'synchrony':>10s} {'corr. index':>12s} {'STTC':>7s}   independent trains, window 10 ms\")\n",
    "for rate in [2, 10, 50]:\n",
    "    s = (rng.random((T, 10)) < rate * 1e-3).astype(float)\n",
    "    sync = float(ba.spike.spike_train_synchrony(s, window_size=10 * u.ms, dt=dt))\n",
    "    ci = float(ba.spike.correlation_index(s, window_size=10 * u.ms, dt=dt))\n",
    "    sttc = np.asarray(ba.spike.spike_time_tiling_coefficient(s, dt=dt, tau=5 * u.ms))\n",
    "    print(f\"{rate:4d}Hz {sync:10.3f} {ci:12.3f} {sttc[np.triu_indices(10, 1)].mean():7.3f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b3028776",
   "metadata": {},
   "source": [
    "With a real shared component the chance-corrected measures rise above zero. Ten trains share 100 events (each train joins 60 % of them, jittered by up to $\\pm 2$ ms) on top of 10 Hz background:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "d22dd093",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-10-15T08:09:51.771899Z",
     "iopub.status.busy": "2026-10-15T08:09:51.771899Z",
     "iopub.status.idle": "2026-10-15T08:09:52.262702Z",
     "shell.execute_reply": "2026-10-15T08:09:52.262702Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "synchrony 0.258, correlation index 0.123, STTC 0.129\n"
     ]
    }
   ],
   "source": [
    "s = (rng.random((T, 10)) < 0.01).astype(float)\n",
    "shared_events = rng.choice(np.arange(10, T - 10), 100, replace=False)\n",
    "for e in shared_events:\n",
    "    who = np.flatnonzero(rng.random(10) < 0.6)\n",
    "    s[e + rng.integers(-2, 3, who.size), who] = 1\n",
    "print(f\"synchrony {float(ba.spike.spike_train_synchrony(s, window_size=10 * u.ms, dt=dt)):.3f}, \"\n",
    "      f\"correlation index {float(ba.spike.correlation_index(s, window_size=10 * u.ms, dt=dt)):.3f}, \"\n",
    "      f\"STTC {np.asarray(ba.spike.spike_time_tiling_coefficient(s, dt=dt))[np.triu_indices(10, 1)].mean():.3f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f2131b59",
   "metadata": {},
   "source": [
    "## 5. Phase Locking to a Rhythm\n",
    "\n",
    "`phase_locking_value(spikes, reference_freq, dt)` assigns each spike the phase $\\phi_k = 2\\pi f t_k$ of a reference oscillation and returns the vector strength $\\lvert\\frac1N\\sum_k e^{i\\phi_k}\\rvert$ per neuron. If spikes occur at a preferred phase with Gaussian timing jitter $\\sigma_t$, the expected value is $e^{-(2\\pi f \\sigma_t)^2/2}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "d8abaafa",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-10-15T08:09:52.263704Z",
     "iopub.status.busy": "2026-10-15T08:09:52.263704Z",
     "iopub.status.idle": "2026-10-15T08:09:52.846915Z",
     "shell.execute_reply": "2026-10-15T08:09:52.846915Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "jitter  1.0 ms: PLV 0.992, expected 0.992\n",
      "jitter  3.0 ms: PLV 0.931, expected 0.931\n",
      "jitter  6.0 ms: PLV 0.752, expected 0.753\n",
      "jitter 10.0 ms: PLV 0.457, expected 0.454\n"
     ]
    }
   ],
   "source": [
    "f_ref = 20.0  # Hz\n",
    "T = 200_000  # 200 s at 1 ms, 4000 cycles\n",
    "jitters = [1.0, 3.0, 6.0, 10.0]  # ms\n",
    "S = np.zeros((T, len(jitters)))\n",
    "cycle_starts = np.arange(0, T, 1000 / f_ref)\n",
    "for j, sd in enumerate(jitters):\n",
    "    times = cycle_starts + 10.0 + rng.normal(0, sd, cycle_starts.size)  # preferred phase 10 ms into each cycle\n",
    "    idx = np.clip(np.round(times).astype(int), 0, T - 1)\n",
    "    S[idx, j] = 1\n",
    "\n",
    "plv = np.asarray(ba.spike.phase_locking_value(S, f_ref, dt))\n",
    "for sd, value in zip(jitters, plv):\n",
    "    print(f\"jitter {sd:4.1f} ms: PLV {value:.3f}, expected {np.exp(-(2 * np.pi * f_ref * sd * 1e-3) ** 2 / 2):.3f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a4918477",
   "metadata": {},
   "source": [
    "`reference_freq` is in Hz; a plain-number `dt` here is in **seconds** (unlike the coincidence measures), which is another reason to pass a quantity. The PLV tests locking to an *external* rhythm; it says nothing about spike-to-spike synchrony.\n",
    "\n",
    "## 6. Burst Synchrony\n",
    "\n",
    "`burst_synchrony_index` first detects bursts in each train (runs of at least `burst_threshold` spikes with inter-spike intervals up to `max_isi`), then returns the fraction of bursts that overlap in time with a burst of **any** other neuron. Like $\\kappa$, it is not chance-corrected: the more neurons and the more bursts, the more likely an overlap by accident. Compare aligned bursts with bursts at independent times, which give the chance level:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "84922630",
   "metadata": {
    "execution": {
     "iopub.execute_input": "2026-10-15T08:09:52.847918Z",
     "iopub.status.busy": "2026-10-15T08:09:52.847918Z",
     "iopub.status.idle": "2026-10-15T08:09:52.859006Z",
     "shell.execute_reply": "2026-10-15T08:09:52.859006Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "aligned bursts: index 1.000\n",
      "independent bursts: index 0.409\n",
      "chance level for independent timing: 0.41\n"
     ]
    }
   ],
   "source": [
    "def bursting(aligned, n=8, T=10_000):\n",
    "    s = np.zeros((T, n))\n",
    "    shared = rng.choice(np.arange(100, T - 100), 15, replace=False)\n",
    "    for i in range(n):\n",
    "        starts = shared if aligned else rng.choice(np.arange(100, T - 100), 15, replace=False)\n",
    "        for b in starts:\n",
    "            s[b + np.arange(0, 25, 5) + rng.integers(0, 2, 5), i] = 1  # 5 spikes, 5 ms apart\n",
    "    return s\n",
    "\n",
    "\n",
    "for aligned in [True, False]:\n",
    "    s = bursting(aligned)\n",
    "    print(f\"{'aligned' if aligned else 'independent'} bursts: index \"\n",
    "          f\"{float(ba.spike.burst_synchrony_index(s, burst_threshold=3, max_isi=10 * u.ms, dt=dt)):.3f}\")\n",
    "\n",
    "# Chance level: two ~25 ms bursts overlap when their starts are within 25 ms (a 50 ms window\n",
    "# out of 10 s); each burst has 7 x 15 = 105 bursts of other neurons it could overlap.\n",
    "print(f\"chance level for independent timing: {1 - (1 - 50 / 10_000) ** 105:.2f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6d39adf5",
   "metadata": {},
   "source": [
    "## Summary\n",
    "\n",
    "- Victor-Purpura and van Rossum distances compare two trains; their parameter ($q$ or $\\tau$) sets the temporal precision that counts as \"the same\".\n",
    "- Raw coincidence fractions (`spike_train_synchrony`) grow with firing rate even for independent trains. STTC and the correlation index stay near 0 for independent trains.\n",
    "- `phase_locking_value` measures locking to an external frequency; with Gaussian jitter it matches $e^{-(2\\pi f\\sigma_t)^2/2}$.\n",
    "- `burst_synchrony_index` compares burst timing across neurons; like $\\kappa$ it has a chance level that grows with the number of neurons and bursts.\n",
    "- Pass times as quantities: plain numbers mean milliseconds for the coincidence and burst measures but seconds for `phase_locking_value`."
   ]
  }
 ],
 "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.11.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}