iss_example.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Example In Situ Sequencing Data\n",
    "Here is a small example on how to run ISTDECO on an In Situ Sequencing dataset from [1]. (References are in the end of the notebook)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Set path variables\n",
    "The data and codebook can be downloaded from [2]."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-02-25T14:31:54.569656Z",
     "start_time": "2021-02-25T14:31:54.554445Z"
    }
   },
   "outputs": [],
   "source": [
    "# Where to load the images and codebook from\n",
    "data_path = '../../Data/ISTDECO_ISS_DATA' \n",
    "\n",
    "# Where to save .csv with results\n",
    "csv_save_path = '../../Results/ISS/results.csv'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-02-25T14:31:59.704592Z",
     "start_time": "2021-02-25T14:31:54.571602Z"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from PIL import Image\n",
    "from scipy.io import loadmat\n",
    "\n",
    "rounds = 5\n",
    "channels = 4 # A, C, G, T\n",
    "tile_height, tile_width = 512, 512\n",
    "ntile_y, ntile_x = 10, 18\n",
    "\n",
    "data_shape = (rounds, channels, ntile_y*tile_height, ntile_x*tile_width)\n",
    "image_data = np.zeros(data_shape, dtype='float32')\n",
    "\n",
    "channels = ['A', 'C', 'G', 'T']\n",
    "# Load images and tile\n",
    "for r in range(rounds):\n",
    "    for c in range(len(channels)):\n",
    "        filepath = f'{data_path}/Round{r+1}/{channels[c]}.tif'\n",
    "        im = np.array(Image.open(filepath), dtype='float32')\n",
    "        image_data[r,c] = im   \n",
    "\n",
    "# Load gene names and codebook\n",
    "temp = loadmat(f'{data_path}/codebook.mat')\n",
    "codebook = temp['codebook']\n",
    "target_names = np.array([temp['target_names'][i][0][0] for i in range(len(temp['target_names']))])\n",
    "codebook = codebook / codebook.sum(axis=(1,2), keepdims=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We need to set appropriate variables before we can run ISTDECO.\n",
    "\n",
    "#### Notation\n",
    "- $r$ Number of imaging rounds\n",
    "- $c$ Number of channels\n",
    "- $m$ Number of codes in codebook\n",
    "- $h$ Image height\n",
    "- $w$ Image width\n",
    "- $\\sigma_{psf}$ Standard deviation of the Gaussians describing the spatial shape of the spots\n",
    "- $Y \\in \\mathbb{R}_+^{r\\times c \\times h \\times w}$ Image data\n",
    "- $D \\in \\mathbb{R}_+^{m \\times r\\times c}$ Codebook (Combinatorial labels)\n",
    "- $b \\in \\mathbb{R}_+^{m \\times r\\times c}$ Background"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-02-25T14:31:59.709120Z",
     "start_time": "2021-02-25T14:31:59.705929Z"
    }
   },
   "outputs": [],
   "source": [
    "r = 5\n",
    "c = 4\n",
    "m = codebook.shape[0]\n",
    "h = 512\n",
    "w = 512\n",
    "sigma_psf = 1.\n",
    "D = codebook\n",
    "Y = image_data\n",
    "\n",
    "# The images are already pre-processed\n",
    "# with a baseline removal, hence we \n",
    "# put the baseline to an small-nonzero\n",
    "# number. Mainly to avoid divisions by\n",
    "# zero in the algorithm.\n",
    "b = 1e-5\n",
    "\n",
    "# Standard deviation of the point-spread-function\n",
    "psf_std = 1.\n",
    "\n",
    "# Number of iterations\n",
    "niter = 75\n",
    "\n",
    "# Suppression radius\n",
    "local_max_radius = 1\n",
    "\n",
    "# We want to use GPU\n",
    "device = 'cuda'\n",
    "\n",
    "# The images are too large to fit all on the GPU at once,\n",
    "# hence we need to tile the images.\n",
    "tiles = [(tiley,tilex) for tiley in range(ntile_y) for tilex in range(ntile_x)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Decode dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-02-25T14:37:40.606766Z",
     "start_time": "2021-02-25T14:31:59.710941Z"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 180/180 [05:23<00:00,  1.80s/it]\n"
     ]
    }
   ],
   "source": [
    "from tqdm import tqdm\n",
    "import pandas as pd\n",
    "from istdeco import ISTDeco\n",
    "\n",
    "results = pd.DataFrame()\n",
    "tau = np.percentile(image_data,99)\n",
    "\n",
    "# Loop over the tiles\n",
    "for iy, ix in tqdm(tiles):\n",
    "    \n",
    "    # Slice the image\n",
    "    sy = slice(iy*tile_height, (iy+1)*tile_height)\n",
    "    sx = slice(ix*tile_width, (ix+1)*tile_width)\n",
    "    \n",
    "    # Pick a tile\n",
    "    tile = image_data[:,:, sy, sx] \n",
    "    \n",
    "    # Create the model\n",
    "    model = ISTDeco(tile,codebook, (psf_std,psf_std), b=b).to(device)\n",
    "    \n",
    "    # Run for few iterations and pick out \n",
    "    # deconvolved images\n",
    "    X, Q, _ = model.run(niter=niter)\n",
    "    \n",
    "    # Pick out barcodes.\n",
    "    # 99th percentile might seem like a strict choice, but recall \n",
    "    # that the images are very sparse, especially X. \n",
    "    code_id, y, x = np.where(np.logical_and(X>tau,Q>0.5))\n",
    "        \n",
    "    intensity = X[code_id, y, x]\n",
    "    quality = Q[code_id, y, x]\n",
    "                                 \n",
    "    # Store in dataframe\n",
    "    df = pd.DataFrame(data={\n",
    "            'x': x.astype('float') + ix*tile_width, \n",
    "            'y': y.astype('float') + iy*tile_height,\n",
    "            'target_id': code_id + 1, # Add a one here so id match with the benchmark\n",
    "            'intensity': intensity, \n",
    "            'quality': quality,\n",
    "            'target_name': target_names[code_id]\n",
    "    })\n",
    "                                 \n",
    "    results = results.append(df)\n",
    "\n",
    "results.to_csv(csv_save_path)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Plot results\n",
    "In the codebook we added nine additional \"fake\", or \"non-targeted\" codes. We can use these codes to assess the quality of the decoded barcodes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-02-25T14:37:45.744489Z",
     "start_time": "2021-02-25T14:37:40.609102Z"
    },
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7fa96d019b38>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmQAAAE9CAYAAACleH4eAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAABjGklEQVR4nO3dd3hU1dbH8e9K7wkJSYAktFADhAChK6KigAgoRURFQQTRa702vOq1e732ix1FUVGwIYJYUFRQASH03nsntNATst8/ZsIbEcIkmcmZM7M+zzPPZE5mzvkNCSt79tlnbzHGoJRSSimlrBNgdQCllFJKKX+nDTKllFJKKYtpg0wppZRSymLaIFNKKaWUspg2yJRSSimlLKYNMqWUUkopiwVZHaA8KleubGrWrGl1DKVUBZo7d+4eY0yi1TnKS+uXUv6npPpl6wZZzZo1ycnJsTqGUqoCichGqzO4g9YvpfxPSfVLT1kqpZRSSllMG2RKKaWUUhbTBplSSimllMVsPYZMKXfIz89ny5YtHDt2zOooqpiwsDBSU1MJDg62OopSXk1rmPcpS/3SBpnye1u2bCE6OpqaNWsiIlbHUYAxhtzcXLZs2UKtWrWsjqOUV9Ma5l3KWr/0lKXye8eOHSMhIUELmRcRERISEvQTv1Iu0BrmXcpav7ymh0xEAoAngRggxxjzgcWRlB/RQuZ97PYz0RqmrGS3/y++riw/D4/2kInIeyKyS0SWnLa9i4isFJE1IjLcubknkArkA1s8mUspb5Kbm0tWVhZZWVlUqVKFlJSUU49PnDjh1mPt37+fN954o9Sve+yxx3jhhRfcmsUOtIYpdW5aw9zD06csRwNdim8QkUDgdaArkAH0F5EMoD4wwxjzT+AWD+dSymskJCSwYMECFixYwLBhw7j77rtPPQ4JCTnr6woKCkp9rLIWMz82Gq1hSpVIa5h7eLRBZoyZDuw9bXMrYI0xZp0x5gQwDscnyy3APudzTrozx8lCw8jpazl6wq27Vcpj3nnnHVq2bEnTpk3p3bs3R44cAWDgwIEMGzaM1q1bc//997N27VratGlDkyZNePjhh4mKijq1j+eff56WLVuSmZnJo48+CsDw4cNZu3YtWVlZ3HfffWd9HsDTTz9NvXr1OO+881i5cmUFvnvv4S01bMrSHfy6cpc7d6mUR2kNKz0rxpClAJuLPd4CtAb+B7wqIucD08/2YhEZCgwFqF69uksHXLB5H89+t4LfVu/hneuzCQsOLGt2pSpEr169GDJkCAAPP/wwo0aN4vbbbwccV1TNmDGDwMBALr/8cu6880769+/PW2+9der1U6ZMYfXq1cyePRtjDD169GD69Ok8++yzLFmyhAULFpT4vMjISMaNG8eCBQsoKCigefPmtGjRosL/HbxUmWtYWeoXwIifV1MpIoSO9ZPKmlmpCqU1rPS8ZlC/MeYIMNiF540Uke1A95CQEJf+dVvUiOfZ3pnc/8Uibvogh5HXtyAixGveuvIij09ayrJtB926z4xqMTzavVGpXrNkyRIefvhh9u/fz6FDh+jcufOp7/Xt25fAQMeHipkzZzJhwgQArrnmGu69917AUaSmTJlCs2bNADh06BCrV6/+WyPgbM/Ly8vjyiuvJCIiAoAePXqU/o37GVdqWFnqF0BG1Rh+XLYTY4wO3lYl0hpm3xpmxbQXW4G0Yo9TndtcZoyZZIwZGhsb6/JrrspO4/k+mcxYu4f+7/zJzoN6Ob3yXgMHDuS1115j8eLFPProo3+5fDoyMvKcrzfG8OCDD54ax7FmzRoGD/57W8HV56m/KFcNK0v9al69EvuO5LNuz2HXUyplIa1hpWdFN9EcoK6I1MJRxK4GrqmIA/fNTiM2PJi7Pl1A91d/5+0BLWhWvVJFHFrZRGk/BXpKXl4eVatWJT8/n48//piUlJQzPq9NmzZ8+eWX9OvXj3Hjxp3a3rlzZx555BGuvfZaoqKi2Lp1K8HBwURHR5OXl3fO53Xo0IGBAwfy4IMPUlBQwKRJk7j55ps9/r5tosJrWJvaCQD8vnoP6YlR53i28mdaw+xbwzw97cVYYCZQX0S2iMhgY0wBcBvwA7Ac+MwYs7SU++0uIiMPHDhQ6kyXNqrCl7e0IyQogH5vz+K939dTWGhKvR+lPOnJJ5+kdevWtG/fngYNGpz1ea+88govvfQSmZmZrFmzhqJel0svvZRrrrmGtm3b0qRJE/r06UNeXh4JCQm0b9+exo0bc9999531ec2bN6dfv340bdqUrl270rJly4p6617FEzWsLPWrZuVIaleO5MdlO0v5DpSyhtaw0hNj7NsYyc7ONjk5OWV67d7DJ7jv84VMXbGLNrXj+U+vTGpVPnc3qvI9y5cvp2HDhlbHKJMjR44QHh6OiDBu3DjGjh3L119/bXUstznTz0ZE5hpjsi2K5DalrV///X4F70xfR87DnYiLOPtUAsr/aA3zTqWtX7ZcOqk8PWRF4iNDePeGbJ7t1YSlWw/S+eXpvPDDSvKO5bsxqVKeNXfuXLKyssjMzOSNN97gxRdftDqSOoey1q8ujapQUGj4ablOf6F8h9aw/2fLSw2NMZOASdnZ2UPKsx8R4epW1bmoQRLPfLuc135Zw5g/NzK0Q21uaFuTyFBb/vMoP3L++eezcOFCq2OoUihr/cpMjSUlLpxJC7fRp0Wqh9IpVbG0hv0/W/aQuVtSTBivXN2Mibe1p1laHM99v5Lzn/uF//20mj2HjlsdTymlEBGuaFaN31bvZscBvUpcKV9jywaZO05ZnklmahzvD2rF+Fvb0TQ1lpd/WkXb/0xlyIc5fLt4O8fydaZ/pVT5lKd+XZWdRqGBcXM2eSCZUspKtmyQlWUen9JoXr0S7w9qxU//vIAb2tZkweb93PrxPFo+/RPDv1zErHW5emWmUqpMylO/aiREckG9RMbO3kT+yUIPpFNKWUUHSZWgTlIUD1+ewYOXNeSPNXuYMH8rExduY9yczSRFh3JxwyQuapBMu/QEHW+mlKoQA9rU4KYPc5i6fCddGle1Oo5Syk1s2UPmqVOWZxMYIHSol8hL/bLIebgT/7s6i+yalZi0cDtDPswh64kp9Ht7Jq/9vJq5G/fpJ1dVaiLCPffcc+rxCy+8wGOPPeaWfU+YMIFly5a5ZV8l2b9/P2+88UapX/fYY4/xwgsveCCRdypv/bqwQRIpceG8+9t6NydTquy0hpW/htmyQebpU5YliQgJomdWCm9c24J5j1zCJze15sbzapF3rIAXpqyi95szyHxsCgNG/clrP69m1rpcHXumzik0NJTx48ezZ88et++7LMWsoKCg1McpazHzN+WtX4EBwtAOtcnZuI+Za3PdnE6pstEaVn62bJB5i5CgANrVqcyDXRvy7Z3nM/fhTrx5bXP6tUxjd95xXpiyiqtHzqLxoz9wxet/8NQ3y/h+yXZ25+mVm+qvgoKCGDp0KC+//PLfvrdhwwYuuugiMjMzufjii9m0yTGge+DAgdxxxx20a9eO2rVr88UXX/zttTNmzGDixIncd999ZGVlsXbtWt555x1atmxJ06ZN6d27N0eOHDm1v2HDhtG6dWvuv/9+1q5dS5s2bWjSpAkPP/wwUVH/v2TP888/T8uWLcnMzOTRRx8FYPjw4axdu5asrCzuu+++sz4P4Omnn6ZevXqcd955rFy50n3/kH6iX8s0kqJDeXHKSuw8ubfyHVrD3MAYY9tbixYtjDfbd/i4+XHpDvPsd8tN3zdnmLoPfWtqPPCNqfHAN+aC5342d38634yZtcEs23bAFJwstDqu31q2bJnVEUxkZKQ5cOCAqVGjhtm/f795/vnnzaOPPmqMMebyyy83o0ePNsYYM2rUKNOzZ09jjDE33HCD6dOnjzl58qRZunSpSU9PP+O+b7jhBvP555+ferxnz55TXz/00ENmxIgRp57XrVs3U1BQYIwxplu3buaTTz4xxhjz5ptvmsjISGOMMT/88IMZMmSIKSwsNCdPnjTdunUz06ZNM+vXrzeNGjU6te+zPS8nJ8c0btzYHD582Bw4cMCkp6eb559//ozZz/SzAXKMF9Sf8t7KW7/GzNpgajzwjflhyfZy7UfZn9Yw76xhpa1fthyJLiLdge516tSxOkqJ4iJC6JSRTKeMZACOF5xkydaDzN24lzkb9jF91W7Gz9sKQFRoEFlpcTSvUYkWNSrRrHocMWHBVsb3T98Nhx2L3bvPKk2g67PnfFpMTAzXX389I0aMIDw8/NT2mTNnMn78eAAGDBjA/ffff+p7V1xxBQEBAWRkZLBzp2vrHC5ZsoSHH36Y/fv3c+jQITp37nzqe3379iUwMPDUcSdMmADANddcw7333gvAlClTmDJlCs2aNQPg0KFDrF69murVq//lOGd7Xl5eHldeeSUREREA9OjRw6XcvsJd9atfdhrv/b6e/3y3ggvqJxIaFOiegMretIbZtobZskFm3DRTf0ULDQqkhbPBNbSDo3dy094jzNu0j7kb9zFv435e+3k1hQZEoH5yNM1rVKJlzUpk14gntZJjvS/lu+666y6aN2/OoEGDXHp+aGjoqa+N89TVQw89xOTJkwFYsGDB314zcOBAJkyYQNOmTRk9ejS//vrrqe9FRp57PVdjDA8++CA333zzX7Zv2LDBpee98sor5zyGL3NX/QoKDODf3Rtxw3uzeXvaOu64uK6bEipVdlrDys6WDTJfISLUSIikRkIkVzZzLIVy6HgBCzfvJ2fDPuZu2sekBdv45E/H+fZqsWG0qhVP69oJtK4VT63KkdpAczcXPgV6Unx8PFdddRWjRo3ixhtvBKBdu3aMGzeOAQMG8PHHH3P++eeXuI+nn36ap59++tTj6Oho8vLyTj3Oy8ujatWq5Ofn8/HHH5OSknLG/bRp04Yvv/ySfv36MW7cuFPbO3fuzCOPPMK1115LVFQUW7duJTg4+G/HOdvzOnTowMCBA3nwwQcpKChg0qRJfyt4yjUX1Evk8syqvPbLGi5rUoU6SdFWR1JW0xp2it1qmDbIvExUaBDt61SmfZ3KAJwsNKzamcecDXv5c/1efl+Ty4QF2wBIig6lde0E2qU7btXjI7SB5gPuueceXnvttVOPX331VQYNGsTzzz9PYmIi77//fqn2d/XVVzNkyBBGjBjBF198wZNPPknr1q1JTEykdevWfylAxb3yyitcd911PP3003Tp0oWiqwIvvfRSli9fTtu2bQGIiopizJgxpKen0759exo3bkzXrl15/vnnz/i85s2b069fP5o2bUpSUhItW7Ysyz+Tcvp39wz+WLOHuz9dyJe3tCMkSK/VUtbSGlY2UtRFaEfZ2dkmJyfH6hgVyhjD2t2H+XN9Ln+u28usdbnscl61mVopnAvrJ3FRwyTa1k4gLFjHlLhi+fLlNGzY0OoYXufIkSOEhztOk48bN46xY8fy9ddfV2iGM/1sRGSuMSa7QoN4gDvr1/dLtjNszDxu7ZjO/V0auGWfyj60hp2Z1TWstPXLlj1kdhnU7wkiQp2kKOokRXFt6xrOBtohZq7NZdqqPXwxdwsfzdpIWHAA7dMrc1HDJC5qkETV2PBz71ypYubOncttt92GMYa4uDjee+89qyP5BE/Ury6Nq3J1yzTe+HUtLWvFc2H9JLftWym7slsN0x4yH3Ms/ySz1uXyy4pdTF2xiy37jgLQsGoMFzVI5KIGSWSlVSIwQE9tFtFPl95Le8hcdyz/JFe8/gc7Dh7jm9vPI7VShNv2rbyb1jDv5Bc9ZOrswoID6Vg/iY71k3ish2HNrkP87GycvTVtHa//spZKEcF0rJ/EhQ2SuKBuIrEROr2GUnYXFhzIW9e1oPurv3PzR3P5Ylg7wkN02IJSdqENMh8mItRNjqZucjQ3X5DOgSP5TFu9m19W7OLXlbv4av5WAgOEFjUq0bVxFS5rUpXkmDCrY1vCGKMXRHgZO/feW6Vm5Uj+1z+LwR/kcO/nC3m1fzMCtDfcL2gN8y5lqV/aIPMjsRHB9GhajR5Nq3Gy0LBg8z5+XrGLn5bt4vFJy3jim2W0S0+gV7NUujSuQmSof/x6hIWFkZubS0JCghY0L2GMITc3l7Aw//yAUB4XNUhmeJcG/Oe7FdROjOSeS+tbHUl5mNYw71LW+qVjyBQAa3YdYuKCrUxYsI1Ne48QERJI18ZVuSo7lVa14n36P3l+fj5btmzh2LFjVkdRxYSFhZGamkpw8F9PqesYsnMzxjD8y8V8mrOZ5/pkclV2mkeOo7yD1jDvU5b6pQ0y9RfGGHI27uPLuVv4ZtF2Dh0voHZiJNe0qk7fFmk63kxZThtkrsk/Wcig9+cwa10uH9zY6tTchkop6/hcg6zYZeNDVq9ebXUcn3XkRAHfLNrO2NmbmL9pP+HBgfRqnsKg9rWokxRldTzlp+zeIKvI+nXgaD5935rB9v3H+GxYWxpWjfHo8ZRSJfO5BlkR7SGrOEu3HWD0Hxv4euE2ThQUcnGDJIZ1TKdlzXiroyk/Y/cGWZGKql9b9x+l9xszKDSGL29pR1q8ToehlFVKql+6xoZySaNqsTzftykzh1/EXZ3qMn/zfvq+NZOr3prJ9FW79Yo4pbxUSlw4H9zYimP5J7nhvdnsPXzC6khKqTPQBpkqlYSoUO7qVI/fH7iQR7tnsHnfEa5/bzZ93prJn+tyrY6nlDqD+lWiefeGlmzdf5RBo+dw5ESB1ZGUUqfRBpkqk4iQIAa1r8Wv93XkqSsas2XfEfqNnMXg0XNYu/uQ1fGUUqdpVSueV/s3Y/GW/dwyZh75JwutjqSUKkYbZKpcQoMCua5NDabddyHDuzZg9vq9dH55Oo9PWsrBY/lWx1NKFXNpoyo8c2UTpq3azQNfLNKhBkp5EW2QKbcICw5k2AXp/HJfR65qmcboGRu4+MVpfLNomxZ9pbzI1a2q889L6jF+/lae/X6F1XGUUk7aIFNuVTkqlGeubMLX/2hPlZgwbvtkPkM+zGF33nGroymlnG6/qA7XtanO29PWMer39VbHUUqhDTLlIZmpcXx1azse7taQ31bvofMr05m2arfVsZRSONa5fbxHY7o0qsKT3yxjwvytVkdSyu95TYNMRDqKyG8i8paIdLQ6jyq/oMAAbjq/NpPvOI+k6FAGvj+bl35cxclCPYWpfI/dalhggPDK1Vm0rZ3AvZ8v5NeVu6yOpJRf82iDTETeE5FdIrLktO1dRGSliKwRkeHOzQY4BIQBWzyZS1WsOknRfHVre3o3T2XE1NXc8N5scg/pKUzl/Xy9hoUFBzLy+hbUTY7m1o/nsWDzfqsjKeW3PN1DNhroUnyDiAQCrwNdgQygv4hkAL8ZY7oCDwCPeziXqmDhIYG80Lcpz/XOZM6GvXQb8TtzN+61OpZS5zIaH69h0WHBfHBjSypHOXqxddoapazh0QaZMWY6cPpf3VbAGmPMOmPMCWAc0NMYUzQpzj4g1JO5lHWuapnG+FvbERIUQL+3Z/HhzA1WR1LqrPylhiVFh/HR4FYEijDw/dl6EY5SFrBiDFkKsLnY4y1Aioj0EpG3gY+A1872YhEZKiI5IpKze7cOErejRtVimXT7eVxQL5F/f72UZ75drlNjKDspcw3z5vpVIyGSUQNbsjvvODeOnsPh4zqbv1IVyWsG9RtjxhtjbjbG9DPG/FrC80YaY7KNMdmJiYkVmFC5U2x4MCOvz2ZAmxqMnL6OB75cRIHOHK5szJUa5u31Kystjtf6N2fptgPc9sk8/T+pVAWyokG2FUgr9jjVuc1lItJdREYeOHDArcFUxQoMEJ7o2Yg7Lq7LZzlbuO2T+RwvOGl1LKXOpVw1zNvrV6eMZJ68ojG/rNzNI18v0d5rpSqIFQ2yOUBdEaklIiHA1cDE0uzAGDPJGDM0NjbWIwFVxRER/nlJPR65PIPvl+5g8OgcPVWivF25apgd6te1rWvwjwvTGTt7M6/9vMbqOEr5BU9PezEWmAnUF5EtIjLYGFMA3Ab8ACwHPjPGLC3lfr36E6YqvcHn1eL5PpnMWLuHa9/9k/1HTlgdSSmP1DC71K97L61Pr2YpvPjjKj7P2XzuFyilykXs3B2dnZ1tcnJyrI6h3OiHpTu4/ZP51KwcwUeDW5McE2Z1JOVlRGSuMSbb6hzlZYf6daKgkEGjZ/Pnur28N7AlHep537g3peykpPrlNYP6S8MunzBV6XVuVIXRg1qydd9R+rw1g425h62OpJRb2al+hQQF8OZ1LaiTFMUtY+ayZKv3Z1bKrmzZILPDGAxVdu3qVObjIW3IO1ZAn7dmsmLHQasjKeU2dqtfMWHBjB7UiriIEAaNnsPW/UetjqSUT7Jlg0z5vqy0OD6/uS0BAv3ensW8TfusjqSU36oSG8boQS05duIkg0fP4ZBeeKOU29myQWanLn9VdnWTo/liWDsqRQRz7Tt/8ttq75pIU6mysGv9qpsczevXNmf1rkPcMXY+JwvtO/5YKW9kywaZ3br8VdmlxUfw2bC21EiIYPDoHKYu32l1JKXKxc71q0O9RB7v0YifV+zisYlLdY4ypdzIlg0y5V+SosMYN7QNDapGM2zMXH5ducvqSEr5reva1GDI+bX4aNZG3v1tvdVxlPIZtmyQ2bXLX5VdXEQIHw1uTb3kaG7+aC4z1uyxOpJSZeIL9evBrg3p1qQqz3y3nClLd1gdRymfYMsGmZ27/FXZxYYH89Hg1o7Tlx/kkLNhr9WRlCo1X6hfAQHCC32bkpkSy53jFuh0GEq5gS0bZMp/xUeG8PFNbagSG8ag0XNYuk3/EChlhfCQQN65IZtKEcHc9EEOOw4cszqSUramDTJlO4nRoYy5qTXRoUHc8N4cnTxWKYskRYcxamBL8o7lc9OHczhyQqfDUKqsbNkg84UxGKp8UuLC+XBwa04WFjJg1Gx2HdRP58oefK1+Nawaw4j+zVi27SD3fb6IQp0OQ6kysWWDzBfGYKjyq5MUxfuDWrHn0HFueH8OB47mWx1JqXPyxfp1ccNkHujSgMmLt/PK1NVWx1HKlmzZIFOqSFZaHG9d14I1u/IY8mEOx/JPWh1JKb80tENt+rZIZcTU1Xw1f4vVcZSyHW2QKdvrUC+RF/o2Zfb6vTw4frFOVqmUBUSEp69sQpva8TzwxWK9ClqpUtIGmfIJPbNSuOeSenw1fyuv/rzG6jhK+aWQoADeuq4F1eLCGDZmHjt1bKdSLrNlg8zXBsUq97jtojr0apbCSz+u4vsl262Oo9QZ+Xr9iosIYeT12Rw5UcDNH83VYQRKuciWDTJfHBSryk9EeKZXE5qmxfHPzxayfPtBqyMp9Tf+UL/qJUfz0lVZLNyyn+FfLtJhBEq5wJYNMqXOJiw4kJEDWhAdFsRNH+SwK09PmShlhS6Nq3B3p3pMWLCNj2ZttDqOUl5PG2TK5yTHhDHqhpbsPXyCW8fM40RBodWRlPJLt11Yh4saJPHEpGXM27TP6jhKeTVtkCmf1Dgllv/2ySRn4z6enrzM6jhK+aWAAOHlq7KoFhfOLWPm6gTOSpVAG2TKZ/VoWo2bzqvFBzM36rxISlkkNiKYtwe04ODRAoaNmcvxAh3kr9SZaINM+bThXRvQqlY8D321hPV7dM1LpazQsGoML/RtyrxN+3likvZYK3Um2iBTPi0oMID/XZ1FSFAAt348Ty/BV8oi3TKrcvMFtfn4z02M0UH+Sv2NLRtkvj6Pj3KvqrHhvNi3Kcu3H+SJb/TTubKWP9ev+zs34ML6iTw+aSnzdZC/Un9xzgaZiPQVkWjn1w+LyHgRae75aGfnD/P4KPe6uGEywy5I55M/N/HNom1Wx1EVyNtqmD/Xr8AA4eV+WSTHhHHrx/PIPXTc6khKeQ1XesgeMcbkich5QCdgFPCmZ2O5WWEhrP4JThyxOomy0D2X1iMrLY5/jV/M1v1HrY6jKo79a9j2RbBjsdUp3CIuIoS3rmvB3sMnuGPcfApO6rQ0SoFrDbKiQTfdgJHGmMlAiOciecCW2fBxb1j6ldVJlIWCnePJCg3cNW4+Jwt19nA/Yf8a9tUw+OUZq1O4TeOUWJ66ojF/rMnluR9WWh1HKa/gSoNsq4i8DfQDvhWRUBdf5z3SWkN8bVg41uokymI1EiJ5omcj5mzYx7u/rbM6jqoY9q9hsSlwwLembumbncZ1baozcvo6Ji/StWeVcqUoXQX8AHQ2xuwH4oH7PBnK7UQg82rY8Dsc2Gp1GmWxK5ulcElGMi9OWcXKHXlWx1GeZ/8aFlMNDvpe7Xrk8gxa1KjEvZ8v1P+Lyu+ds0FmjDkCfA0cFpHqQDCwwtPB3K5JH8DAki+sTqIsJiL8p1cTYsKDuOfzBTqGxcf5RA2LSYUjuZDvWzPdhwYF8sa1zYkOC+LG0XPYnaeD/JX/cuUqy9uBncCPwGTn7RsP53K/hHRIaQGLPrM6ifIClaNCeaxHI5ZsPahzIvk4n6hhsSmOex/sJStaezb38HHu+lTHdir/5copyzuB+saYRsaYJs5bpifCiEikiOSIyOWe2D9N+8POJT5ztZIqn25NqnJencq89OMqvfzet9m/hsWmOe73++aHhyapsTzR0zHI/7Wf11gdRylLuNIg2wyUaQZDEXlPRHaJyJLTtncRkZUiskZEhhf71gOA57qwGveGgGBYoIP7lePU5WM9Mjhy4iQv/rjK6jjKc+xfw+JrOe73rnf7rr1F3xap9GqWwitTV/HLyl1Wx1GqwrnSIFsH/CoiD4rIP4tuLu5/NNCl+AYRCQReB7oCGUB/EckQkUuAZYDn/idGxEO9zo5xZIW6hI6COknRXNemBuNmb2L1Th1U7KPsX8Oiq0FwBOSudfuuvYWI8PSVTWhQJYY7x85n816dN1L5F1caZJtwjL0IAaKL3c7JGDMd2Hva5lbAGmPMOmPMCWAc0BPoCLQBrgGGiIhnLktv0hcO7YT10zyye2U/d1xcl8iQIP77vb3GeSuX2b+GBQQ4pu7Z67sNMoDwkEBGDmgBwO1j55OvF9woPxJ0ricYYx4HEJEo5+ND5TxmCo5TCEW2AK2NMbc5jzMQ2GOMOeP/RBEZCgwFqF69eumPXq8LhMU6TlumX1T61yufEx8Zwi0XpvPc9yuZtS6XNrUTrI6k3Mibali56ld8bdi5tGyJbSQtPoJnejXhtk/m89/vVvDw5RlWR1KqQrhylWVjEZkPLAWWishcEWnkqUDGmNHGmLNeAWWMGWmMyTbGZCcmJpb+AMFhjrFkyyfBcT1FpRxubF+LqrFhPPvdCozRq7x8iTfVsHLVr6SGsG89FPj+BSiXZ1bjhrY1ePf39br2rPIbrnSpjwT+aYypYYypAdwDvFOOY24F0oo9TnVuc5mIdBeRkQcOlGmcruNqy4KjsGxi2V6vfE5YcCC3X1SXBZv389NyHVDsY7yqhpW5fsWngymEfRtK9zqbeqjb/08au2jLfqvjKOVxrjTIIo0xvxQ9MMb8CkSW45hzgLoiUktEQoCrgVK1jIwxk4wxQ2NjY8uWILWlLqWk/qZvdirpiZE8+91ynSzWt3hVDStz/UpId9z78MD+4kKCAnh7QAsSIkO5Zcw89h0+YXUkpTzKpassReQREanpvD2M46qlcxKRscBMoL6IbBGRwcaYAuA2HEuZLAc+M8aUamBEuXvIRBy9ZBt+g/2byrYP5XOCAwN4oEsD1u4+zKjffXd6AT/kVTWs7D1ktR33e/1nDdbKUaG8cW1zducd565PF+ikscqnudIguxFIBMY7b4nObedkjOlvjKlqjAk2xqQaY0Y5t39rjKlnjEk3xjxd2tDl7iEDx9WWAIt1KSX1/y7JSKZzo2Re/HEV63aXd+y38hJeVcPKXL8i4iG8EuT618SpTdPi+Hf3DKat2s2IqautjqOUx7iyluU+Y8wdxpjmztudxph9FRHOo+JrQWorWPQp6CBu5SQiPNmzMaGBATw6cakO8PcBPlXDEur4XYMM4Lo2NejVLIVXf17N7PWnz0KilG84a4NMRF5x3k8SkYmn3yos4Zmzle+UZZHm18PuFbBppnuCKZ+QFBPGPZfW47fVe/hqvu+tHegvvLWGlat+JdTxq1OWxT15RWOqx0dw+9h57DroW4usKwUl95B95Lx/AXjxDDfLuOWUJUDjXhASDfM/dk8w5TMGtK1Jdo1KPDpxqc4Ybl9eWcPKVb8S0h0LjJ847P5gXi4yNIg3rm3BwaMF3PbJfL3wRvmcszbIjDFznV9mGWOmFb8BWRWS7izc1kMWEgkZPWHZBDiu44XU/wsMEF7ulwXAHeO0+NuRt9awcveQgd/2kmVUi+E/vZowe8Nenv9hpdVxlHIrVwb133CGbQPdnKNU3NZDBtDsWjhxyDFRrFLFpMVH8MyVTZi/ab8uPm5vXlXDylW/4oumvvC/cWRFrmiWwnVtqvP29HV8v2S71XGUcpuzLp0kIv1xrMlW67TxFtH8fW03+6reFuJqwIKPIau/1WmUl+netBoz1uby5q9raVUzngsbJFkdSbnIJ2tY0dQXfjIX2dn8+/JGLN5ygPu/WEST1DhS4sKtjqRUuZW0luUMYDtQmb+Ot8gDFnkyVIUSgWYD4JenYP9miEs792uUX3m0ewbzN+3j3s8X8t2d55MUE2Z1JOUa36thoVEQXc3vG2QhQQH87+pmdH/1d4Z+mMMXw9oRHhJodSylyqWkMWQbnTNaXwv8WWzsxXIcS4VYxm1jyIo0usJxr6ct1RmEBQfy2jXNOHyigH9+tpBCnZzSFry1hpW7fiWk+/UpyyI1K0fyv/5ZLNt+kOHjF+kUNcr2XBlD9hlQfETzSeBzz8RxjVvHkAFUrgvVmsPCT9yzP+Vz6iRF81j3Rvy+Zg9vTffv3gkb8qoaVu765adzkZ3JRQ2SueeSeny9YBtj/tRVV5S9udIgCzLGnFpEzPl1iOciWSSzH+xYDDuWWJ1Eeal+LdPollmVF6esYt4me84r6qd8q4Yl1IGje+GIPYfBudutHevQoV4iT32zjGXbDlodR6kyc6VBtltEehQ9EJGewB7PRbJI5lUQGApzR1udRHkpEeGZK5tQJSaMO8bO58DRfKsjKdf4Vg0rmvpCe8kACAgQXrqqKbHhwdw2dh77j+gi5MqeXGmQDQP+JSKbRWQT8ABws2djlcztY8jAsU5cw+6w+DMoOO6+/SqfEhsezKvXNGP7gWPc+7mOJ7MJr6ph5a5fles67vfoVCxFKkeFMqJ/M7bsPcotY+bpIuTKllxZy3KtMaYN0BDIMMa0M8ZY+tHM7WPIijS7Fo4d0MH9qkTNq1fiX5c15MdlOxnz50ar46hz8LYaVu76FVcDAoK0h+w0bWon8NQVjZm5Lpe3dZynsqFzNshEJFlERgGfG2MOiUiGiAyugGwVr1ZHiEmBxV9YnUR5uRvb1+SCeok8NXk5K3fkWR1HlcDnalhgkKNR5udTX5xJ3+xUujWpyss/rmLuRh1jp+zFlVOWo4EfgGrOx6uAuzyUx1oBAY71Ldf8CHk7rU6jvJiI8OJVTYkJC+KuTxdwvOCk1ZHU2Y3G12pYUkPYtczqFF6naJxnckwYd326gEPHC6yOpJTLXGmQVTbGnLps3BhTgOOycd/UYhAUFkDOKKuTKC9XOSqU//bOZPn2gzw+Sf84ejHfq2GJDWDvejipF5acLjYimJf7ZbF131Ee+mqxzk+mbMOVBtlhEUkADICItAHcOJq+9DwyqL9IQjrU7Qw57+ngfnVOFzdMZtgF6Xzy5yYmLtxmdRx1Zl5Vw9xSv+JrgzkJ+3XurTNpWTOeuzs55icb9ft6q+Mo5RJXGmT/BCYC6SLyB/AhcIdHU52Dxwb1F2l9MxzeDSu+8cz+lU+599J6NE2L44lJSzlwRHssvJBX1TC31K/K9Rz3u1e6J5QPuu2iOnRulMx/v1/Bkq2W9iEo5RJXGmRLgQuAdjguFW8ErPBkKMvV7ugY3L9grNVJlA0EBQbw9BWNyT18gie+0VOXXsj3alhCuuN+rw7sPxsR4T+9MkmIDOUfn8zjyAkdT6a8mysNspnGmAJjzFJjzBJjTD4w09PBLBUQCE2vhrVTHQuOK3UOjVNiuf3COnw5bwuT9NSlt/G9GhYRD+HxOvXFOcRHhvDK1VlszD3Cs9/Zuw2ufN9ZG2QiUkVEWgDhItJMRJo7bx2BiIoKaJkWAwGBOe9YnUTZxB0X16VZ9Tj+9dVituw7YnUcv+fzNSwhXae+cEGb2gnc2L4WH87cyLRVu62Oo9RZldRD1hl4AUgFXgJedN7uBv7l+WgWi6sO9bvC/DFwQv+4qnMLCgzgf/2aYQzcOW4BJwoKz/0i5Um+XcPi02HvOqtT2ML9XepTLzmKez5byN7DurSS8k5nbZAZYz4wxlwIDDTGXFjs1tMYM74CM1qnza1wJBeW6ESxyjXVEyJ4plcT5m7cx1OTdTyZlXy+hiXUgYNb4fghq5N4vbDgQF7p14yDR/O5cfQcjuXbe9YT5ZtcGUP2h4iMEpHvALxhlmuPTntRXI12kFDX0UumlIt6NK3GTec5TpF8NEuXVvICXlXD3Fa/ita01HFkLsmoFsOI/lks2LxfPywpr+RKg+x9vGyWa49Pe1FEBJpdB5v/dEzCqJSLHujagPPrVubpycvYmHvY6jj+zqtqmNvqV9GVlvs2lDuTv+jSuCpDO9RmzKxNfL9ku9VxlPoLnan/XBpd6bhf9rW1OZStBAcG8FyfTIICAnh4whKdLdxavlnDKtV03O/TD4ulce+l9WmaGsv9XyzSi2+UV7HlTP0VqlINSMmGRZ+C/lFVpVA1NpwHutTnt9V7GPOnzqhuId+sYaHREJmkpyxLKSQogBH9HRff/OPjeRSc1ItvlHco60z9t3s0lbfJ6u9YyHfnEquTKJu5rk0NWteK55UfV+ks/tbx3RqWkA65eqVladVIiOTpXk1YuOUAI6autjqOUoALDTJjzDxOm+XaGLPI08G8SsaVEBgC8z6yOomyGRHhkcszOHA0n3/pQseW8OkalpCuPWRl1D2zKr2bpzLi5zX8tGyn1XGUKnFi2F5FN6AHUB+oB3R3bvMfkQnQsDssGgf5R61Oo2ymcUos/7y0HpMXb+fbxTusjuM3/KKGJdSBw7vgmP3PwFY0EeHpKxvTOCWGuz9bwPYDWtuVtUrqIevuvA0GRgHXOm/vAjd6PpqXaTHQUfR0cL8qg5s7pNOgSjTPfLtc19SrOL5fwxLqOO51xv4yCQsO5PVrmpN/spAHvlys48mUpUqaGHaQMWYQEAxkGGN6G2N641iYN7iiAnqNmudDfG2Y+4HVSZQNBQYIT/RszNb9R3njF/3jWRH8oobF13bc64z9ZVYjIZKHumUwfdVuHU+mLOXKoP40Y0zxCVt2AtXdHUREGorIWyLyhYjc4u79l5sINL0GNs2Ag7p4tCq9VrXiubJZCiOnr9O5ySqW79YwbZC5xYA2NejVLIURP6/hh6U6rEBZw5UG2VQR+UFEBorIQGAy8JMrOxeR90Rkl4gsOW17FxFZKSJrRGQ4gDFmuTFmGHAV0L50b6OCNO4FCMz70OokyqaGd21AQAA8NXm5DvCvOL5bw4LDIbY67FlVIYfzZc/0akLT1Fju/nQBq3fmWR1H+SFXrrK8DXgLaOq8jTTGuHrJ+GigS/ENIhIIvA50BTKA/iKS4fxeDxzF8lsX91+xEtKh9gWOBtlJHQekSi85Joy7OtXjx2U7+Txni9Vx/ILP17DKdbVB5gZhwYG8PSCbiJBA7v5sga53qSqcKz1kGGO+Msbc7bx95erOjTHTgb2nbW4FrDHGrDPGnADGAT2dz59ojOmKY+DtGYnIUBHJEZGc3bt3uxrFfbIHOxb0Xf9rxR9b+YQh59emVa14nvxmGXsPn7A6jl/wlhrmkfpVuR7sWQ2FOiC9vKrEhvHUFU1YsvUgt348jxMF+m+qKo5LDTI3SwE2F3u8BUgRkY4iMkJE3qaET5fGmJHGmGxjTHZiYqKns/5d3UshIgHmjq74YyufEBggPHVFYw6dKNBBxPZU5hrmkfpVuS7kH4E8HdvqDl0aV+GJno34ecUu3vxVL8BRFSfI6gBFjDG/Ar+68lwR6Q50r1OnjicjnVlwGDTu7ZgktuA4BIVWfAZle/WSoxnQpgajZ2ygU8Nkzqtb2epIqpxcrWFur19FU1/sWQWxqe7Zp5+7vm1NZq/fy2u/rOayJlWomxxtdSTlB0qaGHaq8/6/bj7mViCt2ONU5zaXGWMmGWOGxsbGujWYy9IvhoKjsH66NcdXPmF41wakxYfz5DfLdP4jD/DWGub2+pXYwHG/W8eRudNjPRoRFhTIsDFzOXhMlz1TnlfSKcuqItIO6CEizUSkefFbOY45B6grIrVEJAS4Gsc6cy4Tke4iMvLAAYtmp06/EEKidZJYVS4RIUE8dFkGK3fm8dovuvyNB3hlDXN7/YpKgpAo2LfePftTAFSOCmVE/2ZsyD3CvZ8t1KuilceV1CD7N/AIjk9/LwEvFru94MrORWQsMBOoLyJbRGSwMaYAuA34AVgOfGaMWVqa0Jb3kAWFQkYPWPIlHDl9vK9SruvcKJlezVL439TV/LbagotUfJtX1jC31y8RiK+lc5F5wIUNkniwawOmLNvJhzM3Wh1H+Tg5V6tfRB4xxjxZQXlKJTs72+Tk5Fhz8B2L4a3zoPMz0PYf1mRQPuHoiZNcNuI3AgOEH+7qQGCAWB3Jq4nIXGNMdime75U1zK3167PrYedSuH2ue/anTjHGcOPoOfyxNpdJt51H/So6nkyVXUn1y5V5yJ4UkR4i8oLzdrn7I5aO5acsAao0gdRWkPMeaFe2KofwkEDu71yfNbsO8cmf+inc3bythnmkfsXXhn0bdX5EDxARnuvTlJiwIO4cN1/XolUec84GmYj8B7gTWOa83Skiz3g6WEksP2VZpPn1kLsGNs20NoeyvS6Nq9AuPYHnflhJng4gditvq2EeqV/xtaEwHw7qZMOekBgdyrO9MlmxI4+Xf9SLJ5RnuDIPWTfgEmPMe8aY93DMWm15L5lXaNzLMZhWFxxX5SQiDO/agLxjBbw4RQu+m/l+DdM1LT2uU0Yy/bLTeO+PDczdqGOHlfu5OjFsXLGvLe6W8pJTlgAhkZDZD5aOh6P7rc2ibC8zNY4b2jrmJvt15S6r4/iauGJfW1rDPHbKEiBXJzL1pIcub0i1uDDuGLuAA0e0J1u5lysNsv8A80VktIh8AMwFnvZsrJJ5zSlLgOYD4OQJWPy51UmUD3ioWwbVYsN45afVFBbq2EQ38aoa5pH6FV0VQmNg13L37VP9TUxYMK/2b87Og8cYPn6RToWh3MqVQf1jgTbAeOBLoK0x5lNPB7ONas2galPHguNKlVNIUAB3dqrLgs37+X3NHqvj+AS/qGEijguNdiyyOonPy0qL497O9fluyQ7Gzdl87hco5SJXFxff7lw0d6IxZoenQ52L15yyLNK4t6MQ7tMr5FT5XdEshfjIEN77Qyf6dBdvqmEeq19VsxzT8RSedO9+1d8MPb827dITeGbycjbvPWJ1HOUjrFhcvNy86pQlQN3Ojvu1U63NoXxCaFAgQ86vza8rdzNrXa7VcZSbeax+JTWEgmOwXz8YelpAgPBsr0wA+r41k10Hj1mcSPkCWzbIvE5ifUhsCHNG6Zxkyi0Gta9J1dgwnv1uhY5TUa5JrO+437XC2hx+onpCBJ8MacPeIyd4cvJy/X+qyq3EBpmIBIqI/u8+FxHHbP07l8Dan61Oo3xAWHAgt15YhwWb95OzcZ/VcWzLr2pYUkNAHDP2qwrRJDWWge1qMmnhNt6erlOOqPIpsUFmjDkJrBSR6hWUxyVeN4YMIPMqiKoCs96wOonyEb2bpxAbHsx7v+tYsrLyxhrmsfoVGu2Y/kIH9leo4V0acFmTKrw4ZSVrdx+yOo6yMVdOWVYClorIVBGZWHTzdLCSeN0YMnAsOJ49CNb8BHv1D6gqv4iQIK5pXZ0flu7QgcPl41U1zKP1q0pjx8B+VWECAoTHujciLDiQf3w8j/yThVZHUjblSoPsERyzWj8BvFjspk7XbABIAMwfY3US5SOua1ODQgOf5ejl9eXgPzWsShPYtx6OHbQ6iV9Jignj+T6OpZVG/7HB6jjKplyZh2wasAEIdn49B5jn4Vz2FJsCdS5xNMh0kV/lBilx4XRulMzI6evYsOew1XFsya9qWBXHlX86jqzidW5UhU4Nk3jpx1Us3eZFw2mUbbiyuPgQ4AvgbeemFGCCBzPZW+ZVcGgHbJtvdRLlI/7dvRHGwFvTdFmcsvCrGlalieNeT1tWOBHhmSubEBUWxPAvF3NSV9pQpeTKKct/AO2BgwDGmNVAkidDnYtXDuovUvtCCAqDBXraUrlHSlw4XZtU4av5Wzl0XHtey8CraphH61d0VYhIgJ3aILNCUkwYD3dryOKtB/ho5gar4yibcaVBdtwYc6LogYgEAZY2/b1yUH+RyARo1AsWfALHvLDBqGxpUPtaHC8o5B29tL4svKqGebR+iUDlerrIuIV6NK3GBfUS+c93K1i9M8/qOMpGXGmQTRORfwHhInIJ8DkwybOxbK7FQMeC40vGW51E+YistDguzUjmw5kbOJavS+OUkn/VsMT6jjFkOlGpJUSE5/tmEhESyB3jFpB3LN/qSMomXGmQDQd2A4uBm4FvgYc9Gcr20lo5PqUu/sLqJMqHDGpfi31H8vnf1NVWR7Eb/6ph1ZrBsf2wb4PVSfxWUnQYL/XLYvXOPIZ8mKMfopRLXLnKshD4AHgSeBz4wOgaESUTgUZXwqYZcGCr1WmUj2ibnkD3ptV497d1Oi9ZKfhdDaua5bjfvsDKFH7vwvpJvNC3KbPW7eW571daHUfZgCtXWXYD1gIjgNeANSLS1dPBbC/rGscpg7mjrU6ifMi9l9aj0MAnszdZHcU2/K6GJWVAYIhe6e0FrmiWwjWtqzN6xnrmb9Il0FTJXDll+SJwoTGmozHmAuBC4GXPxvIBlWpC3Uth3gdQcOKcT1fKFTUSIulQtzLj522hQGcEd5V/1bCgEEejbNsCq5MoYHjXBiTHhHHv5ws5qOPJVAlcaZDlGWPWFHu8DrD00hGvnvaiuJY3waGdsPJbq5MoH3JdmxrsPHhce8lc51U1rELqV7Us2L5QB/Z7gZiwYP7bO5O1uw9zx9j5+PLZclU+Z22QiUgvEekF5IjItyIyUERuwHF10pwKS3gGXj3tRXF1LobwSrDoU6uTKB9yUYMk2tSOZ8TUNRzWecnOyltrWIXULx3Y71U61EtkeNcG/LpyN7eMmcfREzrIX/1dST1k3Z23MGAncAHQEcfVSuEeT+YLAgIdc5Kt/QVO6LI3yj1EhPu7NGDPoeO885vOS1YC/61hRQP7dRyZ1xh6fm2GdqjN90t38OIUHeSv/i7obN8wxgyqyCA+K6MH5IyCFZMdyyop5QbNq1fi8syqvP7LGro0rkKDKjFWR/I6fl3Dig/sb9zL6jQKCAgQ/nVZQ/KO5TPqj/VkVIuhV/NUq2MpL+LKVZa1ROQlERkvIhOLbhURzifU7ABx1WHhWKuTKB/zeI9G5J80fDVPp1YpiV/WsKAQqNYcNv5hdRJ1mn9f3oistDju+XyhXnmp/sKVQf0TgA3AqziuViq6KVcEBEDjPrBuGhzZa3Ua5UMSokK5oF4iX87bqgsZl2wC/ljDUrN1xn4vFB4SyEeDW5McHcaT3yyzOo7yIq40yI4ZY0YYY34xxkwrunk8mS+p1wXMSfjzLauTKB/TLbMqew4dZ93uQ1ZH8Wb+WcMS6kDBMV3X0gtFhQYxpENt5m3az0ezNlodR3kJVxpk/xORR0WkrYg0L7p5PJkvSWsFtTvCtP/C+ulWp1E+pG3tBAAmLdxmcRKv5p81LK21414H9nula1tX5/y6lXlkwhIeHL/Y6jjKC7jSIGsCDAGe5f+7+l/wZCifIwK9R0FgKEx9Qk8hKLdJi4+gaWosU1fssjqKN/PPGhZfCxDYqz1k3igsOJC3B7SgcUoMY2dv4s1f9efk71xpkPUFahtjLjDGXOi8XeSJMCJyhYi8IyKfisilnjiGZSIrwyWPw5Y5sGOR1WmUD7m4YTJLtx1k1U5L52v2ZhVSw7yufgWHQ3xtxwSxyitFhAQx4db2XJKRzH+/X8FjE5daHUlZyJUG2RIgrqwHEJH3RGSXiCw5bXsXEVkpImtEZDiAMWaCMWYIMAzoV9Zjeq3MfiCBsPwbq5MoH3J1yzRCgwIYo2NRzqbMNcz29avmeY4LirRX3msFBQbwxrXNubplGqNnbOCucfM5pBM++yVXGmRxwAoR+aGMl4yPBroU3yAigcDrQFcgA+gvIhnFnvKw8/u+JSIeEtJ1HJlyq6SYMFrWjGfG2lyro3irOMpew0Zj5/qV3BjyD8MhPaXtzYIDA3jyisZc07o6ExZs457PFugSS37orBPDFvNoeQ5gjJkuIjVP29wKWGOMWQcgIuOAniKyHMc4j++MMfPKc1yv1eQq+OUpOLAVYlOsTqN8xIUNknjym2XMWLOHdnUqWx3H25S5htm+flWq6bjftx6iky2NokoWHBjAM1c2IbVSOM99v5Lx87bSu4VOHOtPztlDVvwycTdeMp4CbC72eItz2+1AJ6CPiAw70wtFZKiI5IhIzu7du8sZwwIZPRz3v/n+mGJVca5umUZUaBAT9WrLv/FADbNP/Yqv7bjfq0ts2cWwDunUT47m1Z9Xay+Zn3Flpv48ETnovB0TkZMictATYZxzBbUwxgwzxpxx0i5jzEhjTLYxJjsxMdETMTyrcj2oeT7M/QD2bz7385VyQWRoEG1qO05bahH/q4qqYV5ZvyrVcFzdvVMHi9tFQIAwpENtNuQeYdY6nUzcn7jSQxZtjIkxxsTgWJC3N/BGOY+7FUgr9jjVuc0lItJdREYeOHCgnDEsIAI9XwNTCAs+tjqN8iEd6yexae8Rlmz1yOcl2/JADbNP/QoMhqpNHVd3K9vo3CiZ2PBgHpu4lGP5J62OoyqIK4P6TzEOE4DO5TzuHKCuc425EOBqwOULBYwxk4wxQ2NjY8sZwyKVakL1trDMt5fTUxWrW5OqBAiMn7/F6ihey001zF71q1YHR4NMe+RtIzosmCd6NmLlzjyueWcWhbo0ml9w5ZRlr2K3PiLyLHDM1QOIyFhgJlBfRLaIyGBjTAFwG/ADsBz4zBjjcp+6rXvIimT0hF1LYc9qq5MoH1EpMoSLGiTz47KdVkfxKuWpYT5RvzKvcvTIr/f91aJ8Sc+sFB7u1pB5m/Yzfr7LHbDKxlzpIete7NYZyAN6unoAY0x/Y0xVY0ywMSbVGDPKuf1bY0w9Y0y6Mebp0oS2fQ8ZOAb3S4CetlRu1ahaDFv3H9XTHH9V5hrmE/Ur1nml3u4VFXM85TY3tq9Fo2oxvDN9nY4N9QOujCEbVOw2xBjztDHG0kltfKKHLKYapLaCpROsTqJ8SP0q0RgD8zbtszqK1/C2Glbh9SskEtIvghXfVszxlNsEBAj9W1Vn5c48bhkzjxMFhVZHUh501nnIROTfJbzOGGOe9EAelxhjJgGTsrOzh1iVwS0a94bv7oONM6FGW6vTKB/QsX4iYcEB/LBkB+3S/Xs+Mm+tYZbUr8SGjjpjjOPCImUbV2WnsWpnHh/O3MjXC7bSp0Uqoj9Dn1RSD9nhM9wABgMPeDiXf2h2HYTHw+y3rU6ifERESBDN0irx2+o9VkfxBlrDilSqAQVHIW+H1UlUKYUEBTC8awOqx0dw3xeLGDld55TzVWdtkBljXiy6ASNxXC4+CBgH1K6gfGfkE6csAUIioO6ljlMJJw6f+/lKuaBDvUTW7TnMih3+Pf2Ft9YwS+pXknNlp93LK+6Yym0iQoL47s7zOb9uZV6csophH83l+yXbrY6l3KzEMWQiEi8iTwGLcJzebG6MecDqMWQ+Mai/SJM+cPI4zH7H6iTKR/RqnkJQgPBFjk5/4Y01zJL6FVfdca9XddtWZGgQr/TL4vLMqny/dAd3f7qQg8fyrY6l3OisDTIReR7HfDt5QBNjzGPGGB0p7G51OjmWN9nwu9VJlI9IjgnjkoxkvvbzZZS0hhUTmwqhMbBzidVJVDkkRIXyUr8sPryxFUfzT/Kv8Yv16ksfUlIP2T1ANeBhYFuxpUfyPLV0kqt85pQlOAbYJjbUteaUW2VUjWF33nGOF/j19BdeWcMsqV8BgVClCexeWXHHVB7ToV4i3TKr8s2i7Tpe1IeUNIYswBgTXnzZEect2rkEiWV86pQlQFpL2LsW9q63OonyEUkxoQDsOnjc4iTW8dYaZln9SqgDe1Y5rrRUtvdCn6akVgrnzV/XWh1FuUmplk5SHpJ+seN+2zxrcyifkRQTBsCuPP9tkKnTVG8DR/fBmqlWJ1FuEB4SSLcmVcnZuNfvL+DxFdog8wYJ6Y77PWuszaF8RvX4CAB+XqHLKCmnhj0gJAqWfGF1EuUmN51fm9jwEO4YO5+jJ/x6eIJPsGWDzKfGkIFjJu1qzWDZ15B/1Oo0ygekJ0ZxSUYy7/2+gfyTOru3N7GsfoVGQYNusPwbKDhRscdWHpEYHcpLVzVl1c5DPPOtTmlid7ZskPncGDKAljc5Fhv/apjVSZSPuK5NDY7mn+TVn7Xn1ZtYWr/qXgon8nRdSx/SoV4i17auztjZm9i894jVcVQ52LJB5pOaXQetb4FlE2DeR1anUT7ggnqJ9GqWwms/r2bDHp14WAFVmzruty+0Nodyq6EdanPSGG54fzYLN++3Oo4qI22QeZNLn4TkJjDxNvj1v1anUT7gns71CQoMYMRUnRBUAfHpEBAE+/SKbl9SIyGS53pncuBIPj1f/4OZa3OtjqTKQBtk3iQwGG76yXFaYfpzUKiDNFX5pMSF0z2zGt8v3cHh4wVWx1FWCwiA8EpwRP9g+5q+2Wn89M8LABgza6NOGGtDtmyQ+dyg/uKCwxxXQxUWwPrpVqdRPuDihkkcOXGSDbl62tIbWF6/KtWEXJ27yhdVigxhQJsaTF68nRFTdeyo3diyQeaTg/qLa9LHUTS/ux+O7rc6jbK5lLhwAOZv2m9tEAV4Qf1Kagg7l+oEsT7q8R6N6JZZlVd/Xq0LkNuMLRtkPi84HC570bEQ8IhmsG2+1YmUjTVOiaVuUhQT5m+1OoryBikt4Ohe7SXzUQEBwlM9G1OrciSPTVxmdRxVCtog81Z1O8HN06DgOMwfY3UaZWOBAcJ5dSuzaOsBjuXruES/l9bGcb95lrU5lMdUigyhX8s0dhw8xp5DulqHXWiDzJtVbQopzWHbAquTKJtrWTOeEwWFrN55yOooymqV6zlm7N+x2OokyoMaVnUs1/r9kh0WJ1Gu0gaZt4us7Fh/TqlyyEqLQwS+nLfF6ijKagEBEJvmGBKhfFaDKtHER4bw8IQlPPe9TgRsB9og83bJjWHvWse8ZMd0AVlVNtXiwrmudQ1Gz9jA5EU60NfvJWfAvg1Wp1AelBAVyi/3dOS8OpX5LGez1XGUC7RB5u1aDYEGl8Ovz8DHfaxOo2zs0e4Z1EmK4vFJS/lhqZ7G8GuhMdrz7gdiI4JpUaMSuYdPMG3VbqvjqHOwZYPM8nl8KlJYLFz9MbT5B2ydBwe3WZ1I2VRQYACv9MuiUkQIN380lz/X6eSgVvCK+pXU0HGl5Y4l1mVQFeLa1tWpnxzN4NFz2HnwmNVxVAls2SCzfB4fK7QcDAGBMPkenT9IlVnjlFjGDXVcZZezUXtIrOAV9atJX8eHvd9etC6DqhBJMWE82zuTgkLDj8t2Wh1HlcCWDTK/lJAOFz4EK7+FjX9YnUbZWKXIEDJTY3nv9/VMXa4F2i9FxEOjXrDqeygstDqN8rCmqbG0rFmJ/01dTWGhfqD3Vtogs5OWgyEgGFb/aHUSZXMv9m1KYnQogz/I4ZM/N1kdR1khKQPyj0CeDoPwdSJC3+w0ducdZ70uoea1tEFmJyGREF0FDuqM66p86iZH8/Vt7elQL5HHJi5l894jVkdSFa22YyFqZr6uwyD8QHxECACHjhVYnESdjTbI7KZ6W1j5HRzzgwsalEeFBgXyj47pnDhZqAuP+6PE+tD8Bpj1BvzyjNVplIclRDkaZNsP6MB+b6UNMrtp2g9OHILti6xOonxAaHAgACcKdByRX7r8FajdERZ/bnUS5WH1q0STGB3KU5OX6XJKXkobZHYTn+6412VPlBskRDo+NX+zaDtGT1v5n4AAiE11jCVTPi0iJIh3r89md95xHhyvfz+8kdc0yESktoiMEpEvrM7i1eJrQdUsmP+RXh2lyi0tPoK7OtXlq/lbGaOD+8vFtjUssSEc2gkHdFktX9c0LY67L6nHj8t2Ml0nivU6Hm2Qich7IrJLRJactr2LiKwUkTUiMhzAGLPOGDPYk3l8Rtt/wK5lsHS81UmUD7jz4ro0rBrDhPlbtZfsNH5RwxpcBgFBOieZn7ixfS0iQwJ1tQ4v5OkestFAl+IbRCQQeB3oCmQA/UUkw8M5fEvjPpDUCH5/xeokygeICJc1rsLcjfvoNuJ3Ji/azkmdq6jIaHy9hsXXhuwbYe5o2DTL6jTKw0KCAmhdO4EZa3WlDm/j0QaZMWY6sPe0za2ANc5PkyeAcUBPT+bwOQEB0HwA7FwMu1ZYnUb5gFs6pvN8n0yO5Z/kH5/M45KXpjFXZ/L3nxp20SMQVx2+HAL7NlqdRnlYu/QE1u85zOg/1rNm1yHtGfcSVowhSwGKLz2/BUgRkQQReQtoJiIPnu3FIjJURHJEJGf3bj8+B964NwSGwoxXrU6ifEBQYAB9s9P48Z8X8Po1zTly4iQPfbVYC/WZlbmGeW39CouB3qMcC46/2Q7mjNK5yXxYt8yq1EiI4LFJy+j00jRaPj2Vf3w8jw9nbmDljjydzd8iQVYHKGKMyQWGufC8kcBIgOzsbP/9rYlKgpY3wZ9vwnl3Q+U6VidSPiAwQOiWWZWj+Se59/OFfL9kB12bVLU6li24UsO8un6lZsOtM2HibTD5n45llfqPc6yhq3xK1dhwfr23Ixtyj/Dnulz+XL+XP9flMnnxdgAqRQTTu3kqD19u3zPxdmRFD9lWIK3Y41TnNpeJSHcRGXnggJ9PjnreXSABMPd9q5MoH9MzqxoNqkRz/5eLWLPrkNVxvE25aphX16+4NBgwATo9DqunwKLPrE6kPEREqFU5kqtbVeflfln8Mfwifrv/Qp7vk0n1+AjG/KmnriuaFQ2yOUBdEaklIiHA1cDE0uzAGDPJGDM0NjbWIwFtIyoJ6neFheMg/6jVaZQPCQ4M4J3rswkJDOCmD+Zw6Lgut1JMuWqY19cvEWh3B1TJhF//o7XFT4gIafER9M1Oo1PDZI7lF+qE0RXM09NejAVmAvVFZIuIDDbGFAC3AT8Ay4HPjDFLS7lf7/2EWdFaD4MjuTDhVh3zodwqLT6CF/o2ZUPuEX5fvcfqOJbwRA2zRf0KCIBLnoD9m+CrYTrnoZ+pWTkSgJyNp1/PojzJ01dZ9jfGVDXGBBtjUo0xo5zbvzXG1DPGpBtjni7Dfr37E2ZFqnkedHrUMSfZtOesTqN8TNv0BIIChIVb9lsdxRKeqGG2qV/pFzoaZcsmwM9PWp1GVaBLMpKJCQti3OzN536ychuvmalflUP7u6Bpf/j1GfjlP9pTptwmLDiQJqmxfDpnMzN13iL/0+52aDEQfn8J1vxkdRpVQcKCA+mZlcL3S3eQf1J7RyuKLRtktujyr0gi0H0EZF0L056F8UMg/5jVqZSPeL5PU+Iigrlu1J+8+9s6nQqjnGxVv0Sg63OQUAcm36vjyfxIk5RYThQUsn2//i2pKLZskNmmy78iBYVAz9cdEzwu/hw+7AnHbFDwlderkxTF1/9oT6eGSTw1eTl3jFvAsfyTVseyLdvVr6BQ6PYS7FsPv79sdRpVQaonRACwce9hi5P4D1s2yGz1CbMiiUCHex0TPG6eBQs/tTqR8hHRYcG8dV0L7utcn0kLt/HU5GVWR7ItW9av2hdAg8th9js6JMJP1E2KIjBAmLp8l9VR/IYtG2S2+4RZ0Rr3hqgqsDXH6iTKh4gI/7iwDkM71GbMrE1MWrjN6ki2ZNv6VacTHN0Le9dZnURVgISoUPo0T+WT2ZvYcUBPW1YEWzbI1DmIQEoL2DrX6iTKB93XuT4talRi+JeLWLdbJ431G6ktHfdaV/zGbRfVobDQ8Novq62O4he0QearUppD7ho4dtDqJMrHBAcG8Gr/ZoQEBXDrx/N0PJm/SGoIYXHww7/gt5d0jKofKJoodsysTVz0wq888+1yZq/fy0ld69IjbNkgs+UYjIoWX9txv3+TtTmUT6oWF85L/bJYsSOPxyeVal5nv2fb+hUQCAPGO2bwn/o4vNwEpj4Bh7xokXTldo92z+DJno1IqRTO+3+s56q3Z5L91I/c89lCvl+yncO6iofbiJ0vYc/OzjY5OTpO6oy25MC7FzsWB67f1eo0ykf99/sVvPnrWl7u15Qrm6VWyDFFZK4xJrtCDuZBtq5f2+Y7rrhcNhGCwqD59Y45y+LSzv1aZVt5x/KZvmoPPy3fyc8rdnHgaD4hQQG0T0+gU0YynRomkxwTZnVMr1ZS/Qqq6DCqgsQ6C+N+nWlZec49l9Rj7oZ9/Gv8EpqkxFInKdrqSKoiVGsGV30Ie1bD769AzijHrclVcP4/oXJdqxMqD4gOC6ZbZlW6ZVYl/2QhORv28dPynfy4bCe/fLWEh75aQmZqLJ0aJnNxwyQaVokhIECsjm0b2kPmq4yBp6tA9mDo8ozVaZQP23nwGJf97zcyqsXw0eDWHj+e9pB5oQNbYMZrMHc0BAbD3UshLMbqVKqCGGNYs+sQU5bt5KflO1mweT/GQEJkCG3SE2ifXpl26QnUSIhAxL8baD7XQyYi3YHuderUsTqK9xKBxAawS8f3KM9KjgmjT4tU3vtjPYeOFxAVasuyUmF8sn7FpkLXZ6FBN/jgclg7FRpdaXUqVUFEhLrJ0dRNjuYfF9ZhV94xpq3czcy1ufyxdg+TF20HICUunLbpCbSvk0C79Mp6evM0tqycxphJwKTs7OwhVmfxalUaw8rvHb1lfv6pRHnWBfUSeXv6OmauzeWSjGSr43g1n65fNdpBRAKsmKwNMj+WFB1G3+w0+manYYxh3Z7DzFizhxlrc/lp+U6+mLsFgPTESNqlV6Z9nQTa1E4gLiLE4uTWsmWDTLkouQnMHwN5OyCmqtVplA/LrhlPREgg01bt0gaZPwsIhHpdYfkkOJnvOH2p/JqIkJ4YRXpiFAPa1qSw0LBs+8FTvWdfztvCR7M2IgKNqsXQznl6s2XNeCL9rLfdv96tv6nSxHG/c4k2yJRHhQQF0C69Mr+u3I0xxu/Hifi1Bt1gwRjY8DukX2h1GuVlAgKEximxNE6JZUiH2uSfLGTh5v3MWJvLH2v2MPqPDYycvo6gAKFZ9TjaOhtozarHERoUaHV8j9IGmS9LbuS437EY6l5ibRbl8zrWT+Sn5TuZtGg7PZpWszqOskr6hRAUDiu/0waZOqfgwACya8aTXTOeOy6uy9ETJ5m7cR9/rHWc4nzt59WMmLqaqrFhTL7jfOIjffe0pi0bZD45KNYTwuMgJgV2r7Q6ifIDVzZLYcL8rdw5bj4HjuYzoE0NqyN5JZ+vX8HhUL01bJphdRJlQ+EhgZxXtzLn1a0MwIGj+UxbtZu7P13Ai1NW8vSVTSxO6Dm2nKnftovzWqFyPdijDTLleZGhQXw0uDUX1U/ikQlLePnHVdh5Wh1P8Yv6ldYadi6F43lWJ1E2FxseTI+m1RjQpgZjZ29i6TabrXBRCrZskKlSSGwAu1dBYaHVSZQfCA8J5O0BLejTIpX/TV3NI18v0XXv/FFaazCFsGWO1UmUj7i7Uz3iIkJ4fOIyn/2gpw0yX5dYD/IPw8GtVidRfiIoMIDn+2Ry8wW1GTNrE7ePncfxAl2A3K+ktgQENs+2OonyEbERwdx7aX1mb9jLN855zXyNNsh8XWIDx72OI1MVSER4sGtDHu7WkG8X72Dge3PIO5ZvdSxVUcJiHBcVbZpldRLlQ/q1TKNRtRie+XY5R0743qLm2iDzdZXrO+51HJmywE3n1+blfk2Zs2EvV4+cxe6841ZHUhUlrTVsyYFC7R1V7hEYIDzWoxHbDxzjrV/XWh3H7bRB5usiEyCiMuxeYXUS5aeubJbKOzdks273Yfq+NcMnP9mqM6jeBk7kwdqfrU6ifEjLmvH0zKrGW9PXsWXfEavjuJUtG2Qi0l1ERh444LtXW7hVciPYvsjqFMqPXVg/iWd7N2FD7hGWbz9odRxL+U39qtMJKtWEsVfDjFcdS7gp5QYPdGnAyULDhzM3Wh3FrWzZIPOLy8bdKaWFY7b+/GNWJ1F+LD0xCoDcQycsTmItv6lfEfEwdBrU6wJTHoZx18DRfVanUj6gWlw4XRpV4dM5mzl6wndOiduyQaZKKaUFFBY4ZuxXyiIJUY4ZtnMP+3eDzK+Ex0G/MdDlWVj9I7zdAbbOtTqV8gED2tbgwNF8Ji3cZnUUt9EGmT9IzXbcb82xNofya0VLnuQe0oH9fkUE2twCN/4ABhjVGWa9pacwVbm0rhVPveQoPpi5wWfmJdMGmT+IruJYQkk/mSoLhQYFEh0apD1k/iq1Bdw8zTG27PsH4LPr4ZiPj6NTHiMiXN+2Jku3HWTepv1Wx3ELbZD5i5Tm2iBTlouPCvH7MWR+LSIe+o+FS56EFZPh7Qtg2wKrUymburJZCtGhQXw0c4PVUdxCG2T+IqUF7F0HR/ZanUT5sYTIEHIP6ylLvyYC7e+AQd9CwXEYdQnMeVdPYapSiwwNoneLVCYv3u4Tcxxqg8xfpBSNI5tnbQ7l1xKiQrWHTDlUbwPDfodaHWDyPfDFjboYuSq1AW1rkH/S8OmcTVZHKTdtkPmLalmA6GlLZSlHD5k2yJRTZAJc8zlc/G9YNgFGdoQdS6xOpWwkPTGK8+tW5uM/N1FwstDqOOXiNQ0yEYkUkQ9E5B0RudbqPD4nNNqxrqU2yJSFEqJC2Hv4BIWFvnd6SmtYGQUEwPn3wA2T4PghePdimPuBnsJULhvQpgbbDxzjp+U7rY5SLh5tkInIeyKyS0SWnLa9i4isFJE1IjLcubkX8IUxZgjQw5O5/FZqC8fUF1rolEUSIkM5WWg4aJOFxrWGVaCa58Gw3xynMifdAV8NgxOHrU6lbODihsmkxIXzwQx7z9zv6R6y0UCX4htEJBB4HegKZAD9RSQDSAU2O5/mO1PvepOUFnAkF/bb+5dW2VfR5LB77DOObDRawypOVBJcNx46PgiLPoWRF8LuVVanUl4uMEC4tk11Zq7LZfVO+45D9GiDzBgzHTj9sr5WwBpjzDpjzAlgHNAT2IKjoJWYS0SGikiOiOTs3r3bE7F9V5VMx/2u5dbmUH4rLsLRINt/xB4NMnfXMK1fLggIhI7D4foJcGgn/PAvqxMpG+iXnQbAlGX2PW1pxRiyFP7/UyQ4ilgKMB7oLSJvApPO9mJjzEhjTLYxJjsxMdGzSX1NUKjjvrDA2hzKbwWKAI4J222szDVM61cp1O4IVZroaUvlkqKVQPJtPLA/yOoARYwxh4FBrjxXRLoD3evUqePZUEop5SJXa5jWL6XUmVjRQ7YVSCv2ONW5zWXGmEnGmKGxsbFuDaaUUi4oVw3T+qWUOhMrGmRzgLoiUktEQoCrgYml2YGIdBeRkQcO6DpoSqkKV64apvVLKXUmnp72YiwwE6gvIltEZLAxpgC4DfgBWA58ZoxZWpr96idMpVRF8EQN0/qllDoTj44hM8b0P8v2b4FvPXlspZQqL61hStmLnafZ9JqZ+ktDu/yVUnal9Usp9xPnFdx2ZssGmXb5K6XsSuuXUupMbNkgU0oppZTyJbZskGmXv1LKrrR+KaXOxJYNMu3yV0rZldYvpdSZ2LJBpsrJzpehKJ+gv4JKKU+wc2kRY+PKKCK7gY1AZWCPxXHORrOVnTfn02xlV958NYwxtl8Islj9coW3/0w9Qd+zf/C393zW+mXrBlkREckxxmRbneNMNFvZeXM+zVZ23p7PG/njv5m+Z//gj+/5bPSUpVJKKaWUxbRBppRSSillMV9pkI20OkAJNFvZeXM+zVZ23p7PG/njv5m+Z//gj+/5jHxiDJlSSimllJ35Sg+ZUkoppZRt2bpBJiJdRGSliKwRkeFW5ylORNJE5BcRWSYiS0XkTqsznU5EAkVkvoh8Y3WW4kQkTkS+EJEVIrJcRNpanak4Ebnb+TNdIiJjRSTMwizvicguEVlSbFu8iPwoIqud95W8KNvzzp/rIhH5SkTirMhmJ95c586kNL+T4jDC+d4WiUjzYq+5wfn81SJyQ7HtLURksfM1I8S5qrSVv/dnq/e+/L5FJExEZovIQud7fty5vZaI/OnM+amIhDi3hzofr3F+v2axfT3o3L5SRDoX237G3/2zHcP2jDG2vAGBwFqgNhACLAQyrM5VLF9VoLnz62hglTflc+b6J/AJ8I3VWU7L9QFwk/PrECDO6kzFsqUA64Fw5+PPgIEW5ukANAeWFNv2HDDc+fVw4L9elO1SIMj59X+tymaXm7fXuVL83M/4OwlcBnwHCNAG+NO5PR5Y57yv5Py6kvN7s53PFedru5Z0jAp6z2es9778vp05opxfBwN/OvN9Blzt3P4WcIvz61uBt5xfXw186vw6w/l7HQrUcv6+B5b0u3+2Y9j9ZuceslbAGmPMOmPMCWAc0NPiTKcYY7YbY+Y5v84DluP4Y+4VRCQV6Aa8a3WW4kQkFkdBHwVgjDlhjNlvaai/CwLCRSQIiAC2WRXEGDMd2Hva5p44GrU476+oyExFzpTNGDPFGFPgfDgLSK3wYPbi1XXuTEr5O9kT+NA4zALiRKQq0Bn40Riz1xizD/gR6OL8XowxZpZx/DX+8LR9WfJ7X0K999n37cx+yPkw2HkzwEXAF2fIUzznF8DFzl6+nsA4Y8xxY8x6YA2O3/sz/u47X3O2Y9ianRtkKcDmYo+34EUNnuKcXbPNcHyC8BavAPcDhRbnOF0tYDfwvjhOp74rIpFWhypijNkKvABsArYDB4wxU6xN9TfJxpjtzq93AMlWhinBjTg+6auzs02dO4ez/U6e7f2VtH3LGbaXdIwKdVq99+n3LY5hLwuAXTgaj2uB/cU+dBXPeeq9Ob9/AEig9P8WCSUcw9bs3CCzBRGJAr4E7jLGHLQ6D4CIXA7sMsbMtTrLGQThON3xpjGmGXAYRze8V3COz+iJo+FYDYgUkeusTXV2zk/TXncptYg8BBQAH1udRVWsividtOr3vqR674vv2xhz0hiThaOnuxXQoKKO7Yvs3CDbCqQVe5zq3OY1RCQYx3/Oj40x463OU0x7oIeIbMDRDXyRiIyxNtIpW4Atxpii3sQvcDTQvEUnYL0xZrcxJh8YD7SzONPpdjpPceC832Vxnr8QkYHA5cC1zj8g6uy8vs656Gy/k2d7fyVtTz3D9pKOUSHOUu99/n0DOIeV/AK0xXH6NegMOU+9N+f3Y4FcSv9vkVvCMWzNzg2yOUBd59UWITgGCU60ONMpzvPco4DlxpiXrM5TnDHmQWNMqjGmJo5/t5+NMV7Ry2OM2QFsFpH6zk0XA8ssjHS6TUAbEYlw/owvxjFexJtMBIquzroB+NrCLH8hIl1wnCrvYYw5YnUeG/DqOlcKZ/udnAhc77zqsA2OIQDbgR+AS0WkkrNX+lLgB+f3DopIG+f/v+tP25clv/cl1Huffd8ikijOq6RFJBy4BEct/AXoc4Y8xXP2wfF3xzi3X+28CrMWUBfHBQxn/N13vuZsx7A3K68oKO8Nx5Uqq3Cct37I6jynZTsPR9fxImCB83aZ1bnOkLMj3neVZRaQ4/y3m4DzKiNvuQGPAyuAJcBHQKiFWcbiGMuWj6N3cTCOMRZTgdXAT0C8F2Vbg2NcSNH/ibes/nl6+82b61wpfu5n/J3EcaXe6873thjILrafG52/L2uAQcW2Zzv/760FXuP/Jzi37Pf+bPXel983kAnMd77nJcC/ndtr42hQrQE+L6qPQJjz8Rrn92sX29dDzve1EufVoyX97p/tGHa/6Uz9SimllFIWs/MpS6WUUkopn6ANMqWUUkopi2mDTCmllFLKYtogU0oppZSymDbIlFJKKaUspg0yVSYikioiX4vIahFZJyKviUhoOfb3q4hkO7/+VkTinLdb3ZdaKaVKR0QOufCcu0Qkws3HrSYiXzi/zhKRy9y5f+V9tEGmSs05MeF4YIIxpi6OifzCgefcsX9jzGXGMfNzHGBpg6zYbNBKKXU2dwFubZAZY7YZY4omP83CMSeX8mHaIFNlcRFwzBjzPjjWMwPuxjHzdJSIDBSR14qeLCLfiEhH59dvikiOiCwVkcfPtHMR2SAilYFngXQRWSAiz4vIhyJyRbHnfSwiPU97bVURme58zRIROd+5vYuIzBORhSIy1bktXkQmiMgiEZklIpnO7Y+JyEci8gfwkXNG6i9FZI7z1t5N/45KKZsQkY7OnvwvRGSFs/6IiNyBY13bX0TkF+dzLxWRmc6a87k41rgsqm2PO7cvFpEGzu0XOGvWAhGZLyLRIlLTWcNCgCeAfs7v93OemUh0vjZARNYUPVb2pZ/+VVk0Av6yMLkx5qA41sasc47XPmSM2SsigcBUEck0xiw6y3OHA42NY/FaROQCHA2/CSISi2MNyRtOe801OJYaedp5jAhnoXoH6GCMWS8i8c7nPg7MN8ZcISIXAR/i+CQKkAGcZ4w5KiKfAC8bY34Xkeo4ljdpeI73qZTyPc1w1L9twB9Ae2PMCBH5J3ChMWaP88Pkw0AnY8xhEXkA+CeORhXAHmNMc+dwjHuBm5z3/zDG/OFsvB0rOqAx5oSI/BvHbP63ATgbctcCr+BYX3ehMWa3x9+98ijtIVMV7SoRmYdjyY1GOBo+LjHGTMOxtlki0B/40hhTcNrT5gCDROQxoIkxJg9oA0w3xqx37mev87nn4Vj6CGPMz0CCiMQ4vzfRGHPU+XUn4DURWYBj3bWYok+8Sim/MtsYs8UYU4hjeaSaZ3hOGxx17Q9nzbgBqFHs+0ULj88t9vo/gJecvW1xZ6hrp3sPx3qW4Fhq6f1SvQvllbSHTJXFMv5/YVcAnA2ZKjjWImvMXxv7Yc7n1MLxSbClMWafiIwu+l4pfAhch2Oh2UGnf9MYM11EOgDdgNEi8hKwr5THADhc7OsAoI0x5tjZnqyU8gvHi319kjP/DRXgR2NM/3Ps49TrjTHPishkHOPE/hCRzhTrJTudMWaziOx09uy3wtFbpmxOe8hUWUzFcSrwegDnqcEXgdecvUobgCzn2IY0HAUDIAZHQ+eAiCQDXc9xnDwg+rRto3EMoMUYs+z0F4hIDWCnMeYd4F2gOTAL6OBsEFLslOVvOAuZc4zbHmPMwTPkmALcXuwYWefIrZTyL8Vr1SygvYjUARCRSBGpV9KLRSTdGLPYGPNfHL38DUrYf5F3gTHA585xvMrmtEGmSs04VqS/EugjIquBXKDQGPO08yl/AOtx9KSNAOY5X7cQx6nKFcAnzueVdJxcHJ8Wl4jI885tO4HlnL2LviOwUETmA/2A/znHVgwFxovIQuBT53MfA1qIyCIcFxCcPh6tyB1AtnPw/zJgWEm5lVJ+ZyTwvYj84qw3A4Gxztoyk783sE53l7POLQLyge9O+/4vQEbRoH7ntolAFHq60meI42+rUmUnIu2AscCVxph5Hj5WBLAYaG6MOeDJYymllLcSx7yNLxtjzrc6i3IPHUOmys0YM4O/Dlr1CBHpBIzCUYS0MaaU8ksiMhy4BR075lO0h0wppZRSymI6hkwppZRSymLaIFNKKaWUspg2yJRSSimlLKYNMqWUUkopi2mDTCmllFLKYtogU0oppZSy2P8B4gIG9OBEp6IAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 720x360 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "# Name of nontargeted genes\n",
    "non_targeted_names = [ f'blank-{i+1}' for i in range(9)]\n",
    "\n",
    "# Load the raw decoding from the csv\n",
    "results = pd.read_csv(csv_save_path)\n",
    "\n",
    "# Split the dataframe into targeted and non targeted detections\n",
    "ntargeted = results[results['target_name'].isin(non_targeted_names)]\n",
    "targeted = results[~results['target_name'].isin(non_targeted_names)]\n",
    "\n",
    "# Sort data by intensity and quality\n",
    "ntargeted_quality = ntargeted.sort_values('quality',ascending=False)\n",
    "ntargeted_intensity = ntargeted.sort_values('intensity',ascending=False)\n",
    "targeted_quality = targeted.sort_values('quality',ascending=False)\n",
    "targeted_intensity = targeted.sort_values('intensity',ascending=False)\n",
    "\n",
    "# Plot curves\n",
    "plt.figure(figsize=(10,5))\n",
    "plt.subplot(121)\n",
    "plt.semilogy(targeted_quality['quality'].to_numpy(), \\\n",
    "             np.arange(0,targeted_quality.shape[0]), \\\n",
    "             label='Targeted')\n",
    "\n",
    "plt.semilogy(ntargeted_quality['quality'].to_numpy(), \\\n",
    "             np.arange(0,ntargeted_quality.shape[0]), \\\n",
    "             label='Non-targeted')\n",
    "plt.xlabel('Quality score')\n",
    "plt.ylabel('Number of detections')\n",
    "plt.legend()\n",
    "\n",
    "plt.subplot(122)\n",
    "plt.semilogy(targeted_intensity['intensity'].to_numpy(), \\\n",
    "             np.arange(0,targeted_quality.shape[0]), \\\n",
    "             label='Targeted')\n",
    "\n",
    "plt.semilogy(ntargeted_intensity['intensity'].to_numpy(), \\\n",
    "             np.arange(0,ntargeted_quality.shape[0]), \\\n",
    "             label='Non-targeted')\n",
    "plt.xlabel('Intensity')\n",
    "plt.ylabel('Number of detections')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From the above images we note that a quality threshold set to 0.6 seems to filter out many of the non-targeted detections while still keeping many the targeted detections. Let's use this threshold and check correlation with benchmark computed using pipeline in [1]."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2021-02-25T14:51:58.549251Z",
     "start_time": "2021-02-25T14:51:57.926805Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0.5, 1.0, 'r = 0.9761286938922521')"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZEAAAGDCAYAAADqPy3HAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAABBbElEQVR4nO3deXxU1fnH8c9DQgIJIKuURYwKFpfiQkBttS51l1btotBSFrW0ilutRaD91Z9agriLioqCiBtYW6rFBXH7WRdWBRRxiRoUQRRQESKBhOf3xz3RMU7CMGQymeT7fr3uK/eeu8xzZ5J5cs8591xzd0RERJLRJN0BiIhI5lISERGRpCmJiIhI0pREREQkaUoiIiKSNCURERFJmpKIiIgkTUlEMo6Z7W9mC82sNPzcv4Zt9zKzZ8zsCzMrNrNTY9b9xsw2xEylZuZm1jtmmwPN7PmwfrWZXRCz7goze83Mys3sf6u87klm9oKZfW5mH5vZnWbWMmZ9WzObbmZrzWyNmd1nZq2qHOMCM3vfzDaa2TIz2zOUm5n9xcw+MLP1ZjYtdl8zu8rMPgzrlpvZ6CrHPcrMXgnr3zOzYTHrtnXsa8zsHTP70szeNLNBMev2NLOHzexTM1tnZrPM7Psx64eYWUWV9/yIsG5nM3vAzFaGz+pFMzsoZt9OZvZIWO9mVlDdZy51zN01aUrpBGTX4rFygOXAH4Fc4PywnBPvdYG3gYuALOAoYCOwZzXHHgK8C1hYbg98AvwmvFZLYK+Y7QcDJwAPA/9b5Vi/Bo4H8oA2wOPAbTHrJwBPAq2AnYCngOti1p8FLAH2BgzYA2gb87pvArsALcLr3x2z7/eB/DDfBVgK/DwsNwW+AH4fjtsH2ADsl+CxLwN6Ev0DehDwGfDDsK4vcCbQNrzOFcCbVd7fF6p573cPn1On8FkNA9YALcL6jsA5wCGAAwXp/r3WFD67dAegqWFOQAlwSfgiLKutRAIcC3xU+UUfyj4Ajo+z7b7hCzJ22yeBK6o59rPApTHLRcA9CcR0b9UkEmebnwOvxSw/DpwTszwcmBXmmwAfAj+p5lgPAX+OWf4hsAnIi7NtF+A1YERY7hi+hPNitpkPDNjeY4f1jwB/qmZd2/Ba7cJytUmkmv3XA72rlGUridSvSdVZkkoDgJOA1u5eXnWlmS0J1T3xpgnVHHMfYImHb5RgSShPhBEll6qx7Ar8GJgaU3wwsM7MXjKzT8zsP2bWLcHXqerHRFcElW4B+plZGzNrA/yCKLEAdA3TvqFa6n0zu8zMYv9ercp8LtAj5nxGmtkGYAWQD9wP4O6rgQeAoWaWZWaHALsCLyR67JjXaE50JbO06rqYc/7Y3dfGlB0Qqu/eNrP/MbPseDtaVEWZAxRXc2ypJ5REJJXGu/uH7v5VvJXu3svdW1cznVPNMVsQVcfE+oKoqqmqt4iqo/5sZk3N7FjgcKIqpqoGAf919/djyroSVe9cAHQD3if6At4uZnZMOM7fYopfIfqSXBumCqIqrsrXheiq6wfAkUQJ+cxQ/gRwlpkVmNlORFd8xJ6Xu19J9J4cCNzDt9+zB0IsZcB/gb+4+4eJHjvGbcBiYFacc+5KlCgviil+niiB70yUNAcAf46zb6sQ82XuXvWzlnpGSURS6cNtb7LdNhC1I8RqBXxZdUN33wKcQnQ19DHwJ+BBov/OqxoE3F2l7CtghrvPd/dNRO0BPwxfrgkxs4OJrgJ+6e5vx6x6kKi9pmWI/12iarHK1wW4yt0/d/cS4HbgxFA+mSgRPEd0FfBsKP/WeXnk1XC8y0I8PYFp4XxziK7gRpjZSdtzbDO7mighnFblqhAz60BUbTjB3b9Ouu7+nru/7+5b3f014HLgl1X2bQ78B5jj7mORek9JRFKpxiGizWxplZ46sdNt1ey2FOhlZrFVLr2opkrF3Ze4++Hu3s7djyNqwJ1XJY4fAZ2J2gNiLalyDts15LWZHUDUZnCGuz9dZfX+wO3uvtHdNxD9V1+ZJN4CNlf32uFL+FJ3L3D3rkTn/lGY4skmapiH6Iv/bXefFY7zFvAoUQeBhI5tZpeF7Y919/VVzrkNUQJ5xN3HbOMtcmKqzswsF/g3UcL6/Tb2lfoi3Y0ymhrmRNSwfnQKjlvZO+sCorr6c6mmd1bYvhfQjKg65mKiKqncKttMBKbG2fcoot5H+xP1NrqeqMqrcn3TcOz7gb+H+aywbl9gNXB6NXE9C9wENA/TBOClmPVTgZlEVypdiXpMnRnWtSVKCkbUe+t1YFhY14ToC7hNWN8XWAWcH9bvQXQ1dxTf9Poqjtm/2mOH9aOAd4DvxTmnVkQJ+uZqzvkEoGOY7xmOfWnMe/kfoiQStxNGeH/ziZLP94Fm6f4916TeWZpSNKUqiYRjHwAsJKqmeQU4IGbdaODxmOWrQyLYQNRw3b3KsZoBn1N9T6izif4L/yx8ye0Ss25K+EKLnYaEdXcBW8PrVk5LY/bdLRxvLbCOqC2iR8z6VkTVTl8SVQv+jW+6Hu9JdLVSSpRAL4rZr0k41rrwmm+H9yS2h9pp4Qv8S6L/+scBTbZ17LDeidpSYs9rdFg3OKzfWGV9t7D+GqLEuhF4j6g6q2lYd3jYt7TKvodVee1vTen+PdfkX/9SioiIbDe1iYiISNKUREREJGkpTSJm1trMHgpj7Cwzs0MsGjNodhh/Z3bozVE5Zs94i8Y3WmJmB8YcZ3DY/h0zGxxT3tuisYuKw74WLw4REUmNVF+J3Ag84e49gf2AZcBI4Gl37wE8HZYh6rnRI0zDgFshGqgOuJRonJ6+wKWViSds87uY/Y5P8fmIiEiMlCWRcEPWj4FJAO6+2d0/B07mm5u67ia6GYxQPtUjc4DWZtYJOA6Y7e7r3P0zYDZwfFjXyt3neNQ7YGrMsUREpA7EHbemluwGfArcZWb7EXXJvICon/iqsM3HRAPCQTRQXOwdzitCWU3lK+KUf4dFQ10PA8jPz+/ds2fP5M9KRKSRWbhw4Rp37xBvXSqTSDbRuD3nuftcM7uRb6qugKiTt5mlvI+xu08kuqGMwsJCX7BgQapfUkSkwTCz5dWtS2WbyApghbvPDcsPESWV1aEqivDzk7D+I6JnGFTqyjfDLVRX3jVOuYiI1JGUJRF3/xj4MObJZj8B3iAaS6iyh9VgoofeEMoHhV5aBwNfhGqvWcCxMUNmH0v03IVVwHozOzj0yhoUcywREakDqazOAjgPuM/McoiGORhKlLgeNLMziYZVOC1s+xjRAHTFREMfDAVw93VmdgXRg3MALnf3dWH+HKKhJ5oTDWlR+TwGERGpA41u2BO1iYiIbB8zW+juhfHW6Y51ERFJmpKIiIgkTUlERESSpiQiIiJJS3XvLBERSaMNZeXMXLySkrUbKWiXT7/9OtMit/a++pVEREQaqPkl6xhy1zzcoXRzBXk5WVzx6BtMGdqXPgVta+U1VJ0lItIAbSgrZ8hd89hYVkHp5gogSiQbyypCeXmtvI6SiIhIAzRz8Uqquw3QHWYuWVkrr6MkIiLSAJWs3fj1FUhVpZsrKFlTWiuvoyQiItIAFbTLJy8nK+66vJwsCtrn1crrKImIiDRA/fbrTHUPDDeDfr0618rrKImIiDRALXKzmTK0L/m5WV9fkeTlZJGfmxXKa6dzrrr4iog0UH0K2jJv9NHMXLKSkjWlFLTPo1+vzrWWQEBJRESkQcvPzeb0Pt1SdnxVZ4mISNKUREREJGlKIiIikjQlERERSZqSiIiIJE1JREREkqYkIiIiSVMSERGRpCmJiIhI0pREREQkaUoiIiKSNCURERFJmpKIiIgkTUlERESSpiQiIiJJUxIREZGkKYmIiEjSlERERCRpSiIiIpI0JREREUmakoiIiCRNSURERJKmJCIiIklTEhERkaQpiYiISNKUREREJGlKIiIikjQlERERSVpKk4iZlZjZa2a2yMwWhLK2ZjbbzN4JP9uEcjOz8WZWbGZLzOzAmOMMDtu/Y2aDY8p7h+MXh30tlecjIiLfVhdXIke6+/7uXhiWRwJPu3sP4OmwDHAC0CNMw4BbIUo6wKXAQUBf4NLKxBO2+V3Mfsen/nRERKRSOqqzTgbuDvN3A6fElE/1yBygtZl1Ao4DZrv7Onf/DJgNHB/WtXL3Oe7uwNSYY4mISB1IdRJx4EkzW2hmw0JZR3dfFeY/BjqG+S7AhzH7rghlNZWviFMuIiJ1JDvFxz/U3T8ys52B2Wb2ZuxKd3cz8xTHQEhgwwC6deuW6pcTEWk0Unol4u4fhZ+fADOI2jRWh6oows9PwuYfAbvE7N41lNVU3jVOebw4Jrp7obsXdujQYUdPS0REgpQlETPLN7OWlfPAscDrwCNAZQ+rwcDDYf4RYFDopXUw8EWo9poFHGtmbUKD+rHArLBuvZkdHHplDYo5loiI1IFUVmd1BGaEXrfZwP3u/oSZzQceNLMzgeXAaWH7x4ATgWKgFBgK4O7rzOwKYH7Y7nJ3XxfmzwGmAM2Bx8MkIiJ1xKKOTY1HYWGhL1iwIN1hiIhkDDNbGHObxrfojnUREUmakoiIiCRNSURERJKmJCIiIklTEhERkaQpiYiISNKUREREJGlKIiIikjQlERERSZqSiIiIJE1JREREkqYkIiIiSVMSERGRpCmJiIhI0pREREQkaUoiIiKSNCURERFJmpKIiIgkTUlERESSpiQiIiJJUxIREZGkKYmIiEjSlERERCRpSiIiIpI0JREREUmakoiIiCRNSURERJKmJCIiIklTEhERkaQpiYiISNKUREREJGlKIiIikjQlERERSdo2k4iZ5ZtZkzC/p5n9zMyapj40ERGp7xK5EnkeaGZmXYAngd8CU1IZlIiIZIZEkoi5eynwc2CCu/8K2Ce1YYmISCZIKImY2SHAb4BHQ1lW6kISEZFMkUgSuQAYBcxw96VmtjvwbGrDEhGRTJCdwDYd3f1nlQvu/p6Z/TeFMYmISIZI5EpkVIJlIiLSyFR7JWJmJwAnAl3MbHzMqlZAeaoDExGR+q+m6qyVwALgZ8DCmPIvgT+mMigREckM1SYRd18MLDaz+919Sx3GJCIiGSKRNpG+ZjbbzN42s/fM7H0zey/RFzCzLDN71cxmhuXdzGyumRWb2XQzywnluWG5OKwviDnGqFD+lpkdF1N+fCgrNrORiZ+2iIjUhkSSyCTgOuBQoA9QGH4m6gJgWczyOOB6d+8OfAacGcrPBD4L5deH7TCzvYH+RDc4Hg9MCIkpC7gFOAHYGxgQthURkTqSSBL5wt0fd/dP3H1t5ZTIwc2sK3AScGdYNuAo4KGwyd3AKWH+5LBMWP+TsP3JwDR3L3P394FioG+Yit39PXffDEwL24qISB1J5D6RZ83sauBfQFllobu/ksC+NwAjgJZhuR3wubtX9u5aAXQJ812AD8Oxy83si7B9F2BOzDFj9/mwSvlBCcQkIiK1JJEkUvnFXBhT5kRXFNUys37AJ+6+0MyOSCq6WmJmw4BhAN26dUtnKCIiDco2k4i7H5nksX8E/MzMTgSaEd1fciPQ2syyw9VIV+CjsP1HwC7ACjPLBnYC1saUV4rdp7ryqucwEZgIUFhY6Emej4iIVLHNJGJmf4tX7u6X17Sfu48i3NkerkQudvffmNk/gF8StWEMBh4OuzwSll8O659xdzezR4D7zew6oDPQA5gHGNDDzHYjSh79gV9v63xERKT2JFKdtTFmvhnQj2/3ttpelwDTzOzvwKtEvb8IP+8xs2JgHVFSIAz6+CDwBtGd8sPdvQLAzM4FZhGNKjzZ3ZfuQFwiIrKdzH37anfMLBeY5e5HpCSiFCssLPQFCxakOwwRkYxhZgvdvTDeumSesZ5H1P4gIiKNXCJtIq8R9caCqNqoA1Bje4iIiDQOibSJ9IuZLwdWx9znISIijdg2q7PcfTnQGvgpcCrRECMiIiLbTiJmdgFwH7BzmO4zs/NSHZiIiNR/iVRnnQkc5O4bAcxsHNG9HDelMjAREan/EumdZUBFzHJFKBMRkUYukSuRu4C5ZjYjLJ/CNzcIiohII5bI2FnXmdlzRM8TARjq7q+mNCoREckIidwncjCwtHLodzNrZWYHufvclEcnIiL1WiJtIrcCG2KWN4QyERFp5BJqWPeYAbbcfSuJtaWIiEgDl0gSec/MzjezpmG6AHgv1YGJiEj9l0gS+QPwQ6JndlQ+gnZYKoMSEZHMkEjvrE8Iz/YQERGJlcxQ8CIiIoCSiIiI7AAlERERSVoio/h2NLNJZvZ4WN7bzM5MfWgiIlLfJXIlMgWYBXQOy28DF6YoHhERySCJJJH27v4gsBUgPNWwouZdRESkMUgkiWw0s3aE56yHsbS+SGlUIiKSERIZvuQi4BFgDzN7EegA/DKlUYmISEZI5GbDV8zscOD7RA+jesvdt6Q8MhERqfcSGQq+GXAO0fNEHPivmd3m7ptSHZyIiNRviVRnTQW+5Jtnqv8auAf4VaqCEhFJtw1l5cxcvJKStRspaJdPv/060yJXA5hXlcg7sq+77x2z/KyZvZGqgERE0m1+yTqG3DUPdyjdXEFeThZXPPoGU4b2pU9B23SHV68k0jvrldAjCwAzOwhYkLqQRETSZ0NZOUPumsfGsgpKN0d3M5RurmBjWUUoL09zhPVLIkmkN/CSmZWYWQnwMtDHzF4zsyUpjU5EpI7NXLySbx7D923uMHPJyroNqJ5LpDrr+JRHISJST5Ss3fj1FUhVpZsrKFlTWscR1W+JJJHzgUnurnYQEWnwCtrlk5eTFTeR5OVkUdA+Lw1R1V+JVGctA+4ws7lm9gcz2ynVQYmIpEu//TpjFn+dGfTr1Tn+ykZqm0nE3e909x8Bg4ACYImZ3W9mR6Y6OBGRutYiN5spQ/uSn5tFXk4WEF2B5OdmhXJ1842V0LthZllAzzCtARYDF5nZ791dj84VkQalT0Fb5o0+mplLVlKyppSC9nn069VZCSSORO5Yvx7oBzwDFLn7vLBqnJm9lcrgRETSJT83m9P7dEt3GPVeIml1CfBXd98YZ13fWo5HREQySCIN6wOrJhAzexrA3TUkvIhII1btlUgYeDEPaG9mbYhG8AVoBXSpg9hERKSeq6k66/dEj8HtDCzkmySyHrg5tWGJiEgmqDaJuPuNwI1mdp6731TddiIi0nglcp+IEoiIiMSVSMO6iIhIXEoiIiKStG0mEYsMNLO/heVuZrbN+0PMrJmZzTOzxWa21MwuC+W7hXG4is1supnlhPLcsFwc1hfEHGtUKH/LzI6LKT8+lBWb2cgkzl9ERHZAIlciE4BDgAFh+UvglgT2KwOOcvf9gP2B48PDrcYB17t7d+Az4Myw/ZnAZ6H8+rAdZrY30B/Yh2hY+glmlhWGYrkFOAHYGxgQthURkTqSSBI5yN2HA5sA3P0zIGdbO3lkQ1hsGiYHjgIeCuV3A6eE+ZPDMmH9T8zMQvk0dy9z9/eBYqI75fsCxe7+nrtvBqaFbUVEpI4kkkS2hP/6HcDMOgBbEzl4uGJYBHwCzAbeBT5398rnS67gmxsXuwAfAoT1XwDtYsur7FNdebw4hpnZAjNb8OmnnyYSuoiIJCCRJDIemAHsbGZjgBeAokQO7u4V7r4/0JXoyqFnknHuEHef6O6F7l7YoUOHdIQgItIgbXMARne/z8wWAj8humv9FHdftj0v4u6fm9mzRG0rrc0sO1xtdAU+Cpt9BOwCrDCzbGAnYG1MeaXYfaorFxGROlDtlYiZta2ciKqjHgDuB1aHshqZWQczax3mmwPHED0l8Vngl2GzwcDDYf6RsExY/4y7eyjvH3pv7Qb0AOYB84EeobdXDlHj+yMJn7mIiOywmq5EFhK1gxjQjagnlQGtgQ+A3bZx7E7A3aE9pQnwoLvPNLM3gGlm9nfgVWBS2H4ScI+ZFQPriJIC7r7UzB4E3gDKgeHuXgFgZucCs4AsYLK7L92OcxcRkR1k0T/7NWxgdgcww90fC8snEFVp/b4O4qt1hYWFvmDBgnSHISKSMcxsobsXxluXSMP6wZUJBMDdHwd+WFvBiYhI5krkyYYrzeyvwL1h+TfAytSFJCIimSKRK5EBQAeibr4zgJ355u51ERFpxBLp4rsOuMDMWkaLX9+FLiIijVwiAzD+wMxeBV4HlprZQjPbN/WhiYhIfZdIddbtwEXuvqu77wr8CZiY2rBERCQTJJJE8t392coFd38OyE9ZRCIikjES6Z31npn9D3BPWB4IvJe6kEREJFMkciVyBlHvrH8B/wTahzIREWnkEumd9Rlwfh3EIiIiGSaR3lmzKwdSDMttzGxWSqMSEZGMkEh1Vnt3/7xyIVyZ7JyyiEREJGMkkkS2mlm3ygUz25XwlEMREWncEumdNRp4wcz+j2go+MOAYSmNSkREMkKNScTMmhA9YfBA4OBQfKG7r0l1YCIiUv/VmETcfauZjXD3B4GZdRSTiIhkiETaRJ4ys4vNbJcqj8wVEZFGLpE2kdPDz+ExZQ7sXvvhiIhIJknkZsNtPUtdREQaqUSuRAhDv+8NNKssc/epqQpKRKShWL58OS+//DL9+/dPdygpkcgd65cCN4XpSOAq4GcpjktEJKO99dZbDB06lO7du/O73/2O9evXpzuklEikYf2XwE+Aj919KLAfUbdfERGp4tVXX+W0005jr732Yvr06ZxzzjksXbqUVq1apTu0lEikOuur0NW33MxaAZ8Au6Q4LhGRjPLiiy8yZswYHn/8cVq1asWoUaO44IIL2Hnnhj1KVCJJZEEYgPEOYCGwAXg5lUGJiGQCd+fJJ5+kqKiI559/nvbt2zNmzBiGDx/OTjs1jgqbRHpnnRNmbzOzJ4BW7r4ktWGJiNRfW7duZcaMGRQVFfHKK6/QtWtXbrzxRs466yzy8vLSHV6dSrR31s+BQ4nuD3kBUBIRaWQ2lJUzc/FKStZupKBdPv3260yL3IS+QhqMLVu28MADDzB27FjefPNNunfvzp133slvf/tbcnJy0h1eWmzzN8DMJgDdgQdC0e/N7Gh3H17DbiLSgMwvWceQu+bhDqWbK8jLyeKKR99gytC+9Clo+ANYfPXVV9x1111cddVVLF++nF69evHAAw/wq1/9iqysrHSHl1aJ/BtxFLCXuzuAmd0NLE1pVCJSb2woK2fIXfPYWFbxdVnp5mh+yF3zmDf6aPIb6BXJl19+ya233sp1113H6tWrOeSQQ7jllls48cQTMbN0h1cvJNLFtxjoFrO8SygTkUZg5uKVeDVPEHKHmUtW1m1AdWDt2rVceumldOvWjUsuuYRevXrx7LPP8uKLL3LSSScpgcSo9t8HM/sPURtIS2CZmc0LywcB8+omPJHGp761PZSs3fj1lUdVpZsrKFlTWscRpc7KlSu59tpruf3229m4cSOnnnoqo0aNok+fPukOrd6q6TfzmjqLQkSA+tn2UNAun7ycrLiJJC8ni4L2md8b6b333uOqq67irrvuoqKiggEDBjBy5Ej22WefdIdW71WbRNz9/+oyEJHGrr62PfTbrzNXPPpG3HVm0K9X5zqOqPa8/vrrXHnllUybNo2srCyGDh3KiBEj2H13DVKeqETaRESkDtTXtocWudlMGdqX/Nws8nKinkh5OVnk52aF8sxrVJ8/fz6nnnoqP/jBD/j3v//NhRdeyPvvv89tt92mBLKdMu/TF2mg6nPbQ5+CtswbfTQzl6ykZE0pBe3z6Nerc0YlEHfnueeeo6ioiKeeeorWrVvzt7/9jfPPP5927dqlO7yMlch9Is2I7hMBKHb3TakNSaRxqu9tD/m52Zzep9u2N6xn3J1HH32UoqIiXn75ZTp27Mi4ceM4++yzadmyZbrDy3jVVmeZWbaZXQWsAO4GpgIfmtlVZta0rgIUaSz67deZ6nqOZnrbQzpUVFQwffp09t9/f37605+ycuVKJkyYwPvvv8+IESOUQGpJTW0iVwNtgd3cvbe7HwjsAbRGPbdEal1DbHtIh82bNzNp0iR69uxJ//792bx5M3fffTfvvPMOZ599Ns2bN093iA2KeTUteWb2DrCnV9nAzLKAN929Rx3EV+sKCwt9wYIF6Q5DpFoby8ozuu0hXTZu3Midd97JNddcw4oVK+jduzejR4/mlFNOoUkT9SHaEWa20N0L462r6TfTqyaQUFhhZtX0IRGRHZWpbQ/p8vnnn3PLLbdwww03sGbNGn784x8zadIkjjnmGN1ZXgdqSiJvmNmgqs9SN7OBwJupDUtEpGaffPIJN9xwA7fccgvr16/nhBNOYPTo0Rx66KHpDq1RqSmJDAf+ZWZnED2MCqAQaA6cmurARETi+fDDD7nmmmu444472LRpE7/61a8YOXIkBxxwQLpDa5RqumP9I+AgMzsKqLz3/zF3f7pOIhMRifH2228zbtw47rnnHtydQYMGMWLECL7//e+nO7RGraYBGPsA7d39ceCZmPITgE/cfWF1+4qI1JZFixYxduxY/vGPf5Cbm8sf/vAHLr74Yrp1U7tRfVBTl4VxQLwBc94g6v5bIzPbxcyeNbM3zGypmV0Qytua2Wwzeyf8bBPKzczGm1mxmS0xswNjjjU4bP+OmQ2OKe9tZq+FfcabWtFEGoyXXnqJk046iQMOOIAnnniCkSNHsnz5csaPH68EUo/UlERauvvyqoWhrH0Cxy4H/uTuewMHA8PNbG9gJPB06CL8dFgGOAHoEaZhwK0QJR3gUqIh6PsCl1YmnrDN72L2Oz6BuESknnJ3nnzySY444gh+9KMfMXfuXP7+97+zfPlyioqK2HnnndMdolRRUxJpU8O6bY6/4O6r3P2VMP8lsAzoApxMdAc84ecpYf5kYKpH5gCtzawTcBww293XuftnwGzg+LCulbvPCV2Rp8YcS0QyyNatW5kxYwZ9+/bluOOOo7i4mOuvv57ly5fzl7/8hdatW6c7RKlGTUnkKTMbE1tFFKqcLiemjSQRZlYAHADMBTq6+6qw6mOgY5jvAnwYs9uKUFZT+Yo45fFef5iZLTCzBZ9++un2hC4iKVReXs4999zDD37wA37+85/z2Wefcccdd/Duu+9y4YUXkp+fn+4QZRtqSiJ/AnYHis3sn2b2T+AdYM+wLiFm1gL4J3Chu6+PXReuIFJ+46K7T3T3Qncv7NChQ6pfTkS2YdOmTdx222306NGDQYMG0aRJE+677z7efPNNzjrrLHJzc9MdoiSopi6+G4EBZrY733TxXeru7yV68DBQ4z+B+9z9X6F4tZl1cvdVoUrqk1D+EdHz2yt1DWUfAUdUKX8ulHeNs72I1FMbNmzg9ttv59prr2XVqlUcdNBBjB8/npNOOklDk2SomkbxHQgQksY6d/9PZQIxs3O3deBQDTYJWObu18WsegSo7GE1GHg4pnxQqDI7GPgiVHvNAo41szahQf1YYFZYt97MDg6vNSjmWCJSj6xbt47LLruMbt26cfHFF7PPPvvwzDPP8PLLL/PTn/5UCSSTuXvcCXgl3ny85Wr2P5SoqmoJsChMJwLtiHplvQM8BbQN2xtwC/Au8BpQGHOsM4DiMA2NKS8EXg/73EwYULKmqXfv3i4idWPlypV+8cUXe35+vgN+8skn+5w5c9IdlmwnYIFX851a07AnVs18vOV4yemFGrb7SZztnWiolXjHmgxMjlO+ANh3W7GISN16//33ueqqq5g8eTLl5eUMGDCAkSNHsu+++nNtaGocxbea+XjLIiK88cYbjB07lgceeICsrCwGDx7MJZdcwh577JHu0CRFakoiPc1sCdHVxB5hnrCsJ9mLyNcWLFhAUVERM2bMIC8vjwsuuICLLrqILl3i9rqXBqSmJLJXnUUhIhnH3Xn++ecpKiriySefpHXr1vzP//wP559/Pu3bJzKohTQENXXxXQ5gZq2JhhQBeNvdv6iDuESknnJ3HnvsMYqKinjppZfo2LEjV155JWeffTatWrVKd3hSx2oaxTcXuJ1oKJH3iaqxdjWzGcAf3H1znUQoIvVCRUUFDz30EGPHjmXx4sV069aNm2++mTPOOEPPLW/Eauqc/VegKbCLux/g7vsD3YgSz//UQWwiUg9s3ryZyZMns9dee9G/f382bdrElClTKC4uZvjw4UogjVxNbSKnAn3dvbSywN2/NLNzgDkokYg0aKWlpdx5551cffXVrFixggMOOIB//OMfnHrqqWRlZaU7PKknakoiW2MTSCV332Bm6uIr0kB98cUXTJgwgeuvv55PP/2Uww47jDvuuIPjjjsOPbJHqqrxPpEwzEi835qtKYpHRNLk008/5cYbb+Tmm2/miy++4Pjjj2f06NEcdthh6Q5N6rGakshOwELiJxFdiYg0ECtWrOCaa65h4sSJbNq0iV/84heMHj2aAw44IN2hSQaoKYkc7nGebCgiDUNxcTHjxo3j7rvvxt0ZOHAgl1xyCT179kx3aJJBakoiM4ADa1gvIhloyZIljB07lgcffJCmTZsybNgw/vznP7PrrrumOzTJQIkOwCgiGW7OnDmMGTOGmTNn0rJlS/785z9z4YUX8r3vfS/doUkGqymJdDGz8dWtdPfzUxCPiNQid+fpp5+mqKiIZ599lrZt23L55Zdz7rnn0qZNm3SHJw1ATUnkK6KGdRHJMFu3buU///kPRUVFzJs3j86dO3Pdddfxu9/9jhYtWqQ7PGlAakoia9397jqLRER2WHl5OdOnT2fs2LEsXbqU3Xffndtvv53BgwfrueWSEjUNe6KxsUQyRFlZGbfffjvf//73GThwIAD33nsvb731FsOGDVMCkZSp6Uqkv5ntVDlqr5kdSTQY43LgZg3AKJJ+GzZsYOLEiVxzzTWsWrWKvn37cv3119OvXz89t1zqRE2/ZdOBfAAz2x/4B/ABsB8wIeWRiUi11q1bx+WXX86uu+7Kn/70J/baay+eeuop5syZw89+9jMlEKkzNV2JNHf3lWF+IDDZ3a81sybAopRHJiLf8fHHH3Pddddx6623smHDBn72s58xatQoDj744HSHJo1UoveJHAWMAnD3rRqETaRulZSUcPXVVzNp0iS2bNnC6aefzsiRI+nVq1e6Q5NGrqYk8oyZPQisAtoAzwCYWSfU6C4ZZENZOTMXr6Rk7UYK2uXTb7/OtMit6Ve//li2bBlXXnkl9913H02aNGHIkCGMGDGC7t27pzs0EaDmJHIhcDrQCTjU3beE8u8Bf0lxXCK1Yn7JOobcNQ93KN1cQV5OFlc8+gZThvalT0HbdIdXrYULF1JUVMSMGTNo3rw5559/PhdddBFdu3ZNd2gi32LujWtA3sLCQl+wYEG6w5A6sKGsnIOKnmJjWcV31uXnZjFv9NHk16MrEnfnv//9L0VFRcyaNYuddtqJ8847j/PPP58OHTqkOzxpxMxsobsXxltX0zPWvyT+kO8GuLu3qqX4RFJi5uKVVPc/kjvMXLKS0/t0q9ug4sbiPP744xQVFfHiiy/SoUMHxo4dyznnnEOrVvozk/qt2iTi7i3rMhCR2laydiOlm797FQJR1VbJmu88uLNOVVRU8K9//YuioiIWLVrELrvswvjx4znzzDPJy8tLa2wiiVJncmmwCtrlk5cT/1ngeTlZFLRPzxf1li1bmDJlCvvssw+nnXYapaWlTJ48meLiYs477zwlEMkoSiLSYPXbrzPV9UY3g369OtdpPF999RU333wz3bt3Z+jQoTRv3pzp06fzxhtvMHToUHJycuo0HpHaoCQiDVaL3GymDO1Lfm7W11ckeTlZ5OdmhfK6aVRfv34948aNo6CggPPOO4+uXbvy6KOP8sorr3DaaaeRlRX/akkkE9SfrikiKdCnoC3zRh/NzCUrKVlTSkH7PPr16lwnCWTNmjXceOON3HTTTXzxxRcce+yx/OUvf+Gwww5DN+xKQ6EkIg1efm52nfbCWrFiBddeey0TJ06ktLSUn//854wePZrevXvXWQwidUVJRKSWFBcXc9VVVzFlyhS2bt3Kb37zGy655BL23nvvdIcmkjJKIiI76LXXXmPs2LFMnz6dpk2bcuaZZzJixAh22223dIcmknJKIiJJmjt3LkVFRTzyyCO0aNGCP/3pT/zxj3+kU6dO6Q5NpM4oiYhsB3fnmWeeoaioiGeeeYa2bdty2WWXce6559K2bf0di0skVZRERBKwdetWZs6cSVFREXPnzqVTp05ce+21DBs2jBYtWqQ7PJG0URIRqUF5eTkPPvggY8eO5fXXX2e33XbjtttuY/DgwTRr1izd4YmknZKISBxlZWVMnTqVcePG8e6777L33ntzzz330L9/f7Kz9WcjUkl/DSIxNm7cyMSJE7nmmmtYuXIlhYWFzJgxQ88tF6mGkogI8Nlnn3HLLbdwww03sHbtWo444gimTJnC0UcfrbvLRWqgJCKN2urVq7n++uuZMGECX375Jf369WPUqFH88Ic/THdoIhlBSUQapeXLl3PNNddw5513snnzZk477TRGjhzJfvvtl+7QRDKKkog0Km+99RZXXnkl9957L2bGoEGDuOSSS+jRo0e6QxPJSClrKTSzyWb2iZm9HlPW1sxmm9k74WebUG5mNt7Mis1siZkdGLPP4LD9O2Y2OKa8t5m9FvYZb6q4lhq8+uqrnHbaaey1115Mnz6d4cOH8+6773LnnXcqgYjsgFR2N5kCHF+lbCTwtLv3AJ4OywAnAD3CNAy4FaKkA1wKHAT0BS6tTDxhm9/F7Ff1tUR44YUXOPHEEznwwAOZNWsWo0aNoqSkhBtuuIFddtkl3eGJZLyUJRF3fx5YV6X4ZODuMH83cEpM+VSPzAFam1kn4Dhgtruvc/fPgNnA8WFdK3ef4+4OTI05ljRy7s4TTzzBj3/8Yw477DDmz5/PmDFj+OCDDxgzZgw777xzukMUaTDquk2ko7uvCvMfAx3DfBfgw5jtVoSymspXxCmPy8yGEV3h0K1b3T1XQuLbUFbOzMUrKVm7kYJ2+fTbrzMtauEhUVu3bmXGjBkUFRXxyiuv0LVrV2688UbOOussPbdcJEXS1rDu7m5mXkevNRGYCFBYWFgnrynxzS9Zx5C75uEOpZsryMvJ4opH32DK0L70KUhuAMMtW7Zw//33c+WVV/Lmm2/So0cPJk2axMCBA/XccpEUq+tbcFeHqijCz09C+UdAbAV111BWU3nXOOVSj20oK2fIXfPYWFZB6eYKIEokG8sqQnn5dh3vq6++4pZbbqF79+4MGTKEnJwcpk2bxrJlyzjjjDOUQETqQF0nkUeAyh5Wg4GHY8oHhV5aBwNfhGqvWcCxZtYmNKgfC8wK69ab2cGhV9agmGNJPTVz8Uq8mutAd5i5ZGVCx1m/fj1XXXUVu+22G+eeey5dunRh5syZLFq0iNNPP52srKxajFpEapKy6iwzewA4AmhvZiuIelldCTxoZmcCy4HTwuaPAScCxUApMBTA3deZ2RXA/LDd5e5e2Vh/DlEPsObA42GSeqxk7cavr0CqKt1cQcma0hr3X7NmDePHj+emm27i888/5+ijj2batGkcfvjhGppEJE1SlkTcfUA1q34SZ1sHhldznMnA5DjlC4B9dyRGSZ14jecF7fLJy8mKm0jycrIoaB+/8fujjz7i2muv5fbbb6e0tJRTTz2VUaNG0adPn1Sfhohsg+5Yl1pXXeP5rb/pTXUXDGbQr1fnb5W9++67XHXVVUyZMoWKigoGDBjAyJEj2WeffergLEQkERrbWmpVTY3nZ9+3kFsH9iY/N4u8nKjdIi8ni/zcLKYM7Ut+6Ob7+uuvM3DgQPbcc0+mTJnCGWecwdtvv80999yjBCJSz+hKRGrVthrPV33+FfNGH83MJSspWVNKQfs8+vXqTH5uNvPmzaOoqIiHH36Y/Px8LrroIi666CI6depUtychIglTEpFalUjjeX5uNqf3iW76dHeee+45xowZw9NPP02bNm249NJLOe+882jXrl1dhi4iSVASkbhququ8ct3bq7/k89IttM5ryp4dW25X47m7M3PmTIqKipgzZw7f+973uPrqq/n9739Py5Yt6/RcRSR5SiLyHTXdVQ4w5K55lFc4ZeVbv94nN7tJQo3nJ+zTkWnTplFUVMRrr71GQUEBEyZMYOjQoTRr1qwuTk9EapGSiHzL6vWb+M2dc9hc/k3DRuVVxZDJ83Cc0s1bv7NfWflWysr5uvH87HsXfisJGfCL9qvpvd++FBcXs9deezF16lT69+9P06ZN6+r0RKSWKYnI1+aXrGPgnXO/lUBiban4bvKoqmrj+TurPmf50gU8MbGIK0repXfv3vzzn//klFNOoUkTdQ4UyXRKIgJ80zU3toqqqs0V2x67srLxfMtXG3h39r3ceP31rFmzhsMPP5xJt0/gmGOO0d3lIg2IkogANXfNrZSTFX3515RMmjdtwivPP8Gu/S9l/fr1nHTSSYwaNYof/ehHtRmuiNQTSiINXKLP7qipa26lpllNcLzGJLJx4wYeumU0vzzlp4waNYr9999/R09BROoxJZEGbHue3VFT11yAnGxjyhnV987auqUMtlZQWDqPp15bxJ577pm6ExOResN8W3UYDUxhYaEvWLAg3WGk3Iaycg4qeoqNZd9NCvk5WYw4/vus+mLT11cnQLXb52Y34b8jjmTnVlEX3I1l5UyYOZeHZr/Iux+uwjaX8sN9d2fcuafRc4/dUntiIlLnzGyhuxfGW6crkQZq5uKVlFdT7bRxcwVjHl3G5gr/1tXJlKF9v3PlYgZThvb9OoG89NJLjBkzhscee4xWrVoxfPhwLrxwlJ5bLtJIKYk0UG+v/jKhnlZf3wNy1zzmjT467rhWeTlZPPnkkxQVFfF///d/tG/fnr///e8MHz6c1q1b18XpiEg9pSTSQH1eumW7tq98suDpfbp9Pa7V1q1b+fe//01RURELFy6kS5cu3HDDDZx11lnk5+enImwRyTBKIg1U67ztuws89smCW7ZsYdq0aYwdO5Zly5bRvXt37rzzTgYOHEhubm4qwhWRDKVbhhuoPTu2JDc78Y83LyeLLjvlcOutt7LnnnsyaNAgsrOzuf/++1m2bBlnnnmmEoiIfIeSSAN1ZM+dcRLvebdl82b+fPpRnHPOOXTs2JFHHnmExYsXM2DAALKzdcEqIvHp2yHDJHLzYOX9IYZBTCLJzW7y9Qi7Tcwo3VxBNhVs2byZDx+8lEN7duf+qXdxxBFHaGgSEUmIkkgGqUwOW7fCV1sqyG5iXPrIUsYPOICPP/+KZ976hPIKZ37JumrvKn/+z0fy5WdrGHHzdP776jK++vQDjti9FX/952QOOuigOj4jEcl0utkwQ9R082CimmUb3dcvYvbtl1FeXs6AAQMYOXIk++67by1GKiINjW42bABmLl7J1m2PxF6jTeXOnNeLGTJkCCNGjGCPPfaoneBEpNFSEskQJWs38tWW5K9CALK8nL+PvJA/HNurlqISkcZOvbMyREG7fLKb7Fhjd7Nmufz28L1rKSIRESWRemdDWTnT5n3AlY8vY9q8D9hQVg5Av/06k7XdSSRq78rLySI/N4spQ/uSH2cYeBGRZOkbpZ7YUFbOzc+8w53/fY8mZt8ZHLFPQVvuGFTIoMnzEj5mTnYWgw7elR4dW9CvV2clEBGpdeqdVQ/ML1nHkMnz2FjNszxys5sw8OBu9Ni5JW3zczj73oVs60m1zZo24Z4zD/rOc0NERLaXemfVY5XPNq8ugQCUlW9l0gsl0dDsQK/mn/HK+hZY9nfHx2oC/HT/ThSd2ktXHiKScvqWqSNV7zQ/sufOPPvmJzz++iq21DBke6zKYdsXbsilSZP4+zTPzVICEZE6o2+aOlD1MbW52U0Y+a/XyM1uUuMzP6rTvHkeQw/dnXvmLI/7ACklEBGpK/q2SbGvq6ti7jSvTBzJJBCAsgowLO4DpJRARKQu6RsnxWYuXklt913Iy8mioH0e+bnZXz9ASkQkHXSfSIq9vfrLr9syaosZ9OvVuVaPKSKSDF2JpNDzb3/K1JeX7/BxcrK+uW9E7R4iUp/omygFVq/fxAUPvMqc99ft8LFyspvw15N6surzMrV7iEi9o2+jWrShrJy/zHiNhxetrLVjnrjv9xh0yG61djwRkdqkJFJL5pesY/CkuZRu2cHx2mM0b9qEQ/ZoV2vHExGpbUoiSdpQVs7kF9/jvpc/YOPmcko3V7C1lnthNWliakAXkXpNSSQJDy/6iAumLar146oBXUQyjb6htsPq9Zv46U3P88mXW1JyfDPjrEMLNOquiGQMfUslaNAdL/H8u5/t8HGymxjZTYxNce5Wz2pi9OjYQjcQikjG0M2GCSgY+WitJJBmTZtwUq9OcRMIRGNglawp3eHXERGpKxmfRMzseDN7y8yKzWxkbR+/YOSjSe1X+QzCNnlNOWT3dlz+s31Y+NdjOGT3duTlZMXdp3I4ExGRTJHR1VlmlgXcAhwDrADmm9kj7v5GbRw/2QQC8PMDu3D5yft+p12j336dueLR+OFpOBMRyTSZfiXSFyh29/fcfTMwDTg5zTHRtInFTSAALXKzQ6+rrK+vSPQMdBHJVJn+jdUF+DBmeQVwUNWNzGwYMAygW7fUN1pPGtKnxmTQp6CthnEXkQahUXxruftEYCJEz1hP1etkG0w+oy8/7tFhm9tqGHcRaQgyPYl8BOwSs9w1lNUpA35/+O6cd1QPXU2ISKOS6d9484EeZrYbUfLoD/y6tg5ecuVJ22xc/4OSh4g0Yhn9zefu5WZ2LjALyAImu/vS2nyN6hLJ5SfvzS8O3EXJQ0QatYz/BnT3x4DHUvkaJVeelMrDi4hkrEzv4isiImmkJCIiIklTEhERkaQpiYiISNKUREREJGlKIiIikjQlERERSZqSiIiIJE1JREREkmbuKRvUtl4ys0+B5Unu3h5YU4vh1Bc6r8zTUM9N51U/7erucYcnb3RJZEeY2QJ3L0x3HLVN55V5Guq56bwyj6qzREQkaUoiIiKSNCWR7TMx3QGkiM4r8zTUc9N5ZRi1iYiISNJ0JSIiIklTEkmAmR1vZm+ZWbGZjUx3PIkwsxIze83MFpnZglDW1sxmm9k74WebUG5mNj6c3xIzOzDmOIPD9u+Y2eA0nctkM/vEzF6PKau1czGz3uG9Kg77WhrP63/N7KPwuS0ysxNj1o0KMb5lZsfFlMf9/TSz3cxsbiifbmY5dXReu5jZs2b2hpktNbMLQnlGf2Y1nFfGf2Y7xN011TARPXb3XWB3IAdYDOyd7rgSiLsEaF+l7CpgZJgfCYwL8ycCjwMGHAzMDeVtgffCzzZhvk0azuXHwIHA66k4F2Be2NbCviek8bz+F7g4zrZ7h9+9XGC38DuZVdPvJ/Ag0D/M3wacXUfn1Qk4MMy3BN4O8Wf0Z1bDeWX8Z7Yjk65Etq0vUOzu77n7ZmAacHKaY0rWycDdYf5u4JSY8qkemQO0NrNOwHHAbHdf5+6fAbOB4+s4Ztz9eWBdleJaOZewrpW7z/HoL3dqzLFSqprzqs7JwDR3L3P394Fiot/NuL+f4T/zo4CHwv6x71FKufsqd38lzH8JLAO6kOGfWQ3nVZ2M+cx2hJLItnUBPoxZXkHNvzj1hQNPmtlCMxsWyjq6+6ow/zHQMcxXd471+dxr61y6hPmq5el0bqjWmVxZ5cP2n1c74HN3L69SXqfMrAA4AJhLA/rMqpwXNKDPbHspiTRch7r7gcAJwHAz+3HsyvAfXIPomteQzgW4FdgD2B9YBVyb1mh2gJm1AP4JXOju62PXZfJnFue8GsxnlgwlkW37CNglZrlrKKvX3P2j8PMTYAbRJfTqUBVA+PlJ2Ly6c6zP515b5/JRmK9anhbuvtrdK9x9K3AH0ecG239ea4mqhbKrlNcJM2tK9EV7n7v/KxRn/GcW77waymeWLCWRbZsP9Ai9JnKA/sAjaY6pRmaWb2YtK+eBY4HXieKu7OEyGHg4zD8CDAq9ZA4GvgjVDrOAY82sTbhEPzaU1Qe1ci5h3XozOzjUSQ+KOVadq/ySDU4l+twgOq/+ZpZrZrsBPYgal+P+fob/9J8Ffhn2j32PUn0OBkwClrn7dTGrMvozq+68GsJntkPS3bKfCRNR75G3iXpU/CXd8SQQ7+5EPT4WA0srYyaqc30aeAd4Cmgbyg24JZzfa0BhzLHOIGoQLAaGpul8HiCqJthCVE98Zm2eC1BI9If/LnAz4SbcNJ3XPSHuJURfQp1itv9LiPEtYnojVff7GX4P5oXz/QeQW0fndShRVdUSYFGYTsz0z6yG88r4z2xHJt2xLiIiSVN1loiIJE1JREREkqYkIiIiSVMSERGRpCmJiIhI0pREpNaYWUUYxXSxmb1iZj+s5eMfYWYza/OYNbzWhm2sL7CY0Xfrmpk9Z2YJP7PbzIaYWedUxpRq4T3/dbrjkG9TEpHa9JW77+/u+wGjgLHpDmh7hRveUv53EXNXcl0ZAmR0EgEKACWRekZJRFKlFfBZ5YKZ/dnM5odB6i4LZQVmtszM7rDo+QxPmlnzsK67mT0Vc1WzRzhUCzN7yMzeNLP7wl3Elc9PGRuuhBaY2YFmNsvM3jWzP4RtWpjZ0+F4r5nZyTFxvGVmU4luYNslJu72ZvaymZ0U5xyzQwzLQkx5YZ+/hXN93cwmxsT4nJndYNHzXS4wsz5m9lI4x3lm1tLMmpnZXSG+V83syLBvczObFl5rBtA8JsZjQ4yvmNk/LBrbiZj1vyS6Oe++8P6cZGb/jll/TDgmZrbBzK4Pn8fTZtYhlO9hZk9YNKDnf82sZ9U3I7y/lbEvMbNfhPIBoex1MxsXs/2GmPlfmtmUMD/FomeEvGRm74X4Aa4EDgvn8Ecz2ye8b4vC6/WI8xlJqqX7bkdNDWcCKoju4n0T+ALoHcqPJXrGtBH94zKT6FkaBUA5sH/Y7kFgYJifC5wa5psBecAR4bhdw3FeJhpoEqLnp5wd5q8nunu4JdABWB3Ks4mGEAdoT3RXsIU4tgIHx5zLBqJRZucCx8Q51wKiu5d/FJYnE54pQbgTO8zfA/w0zD8HTAjzOUTPx+gTlluF+P4ETA5lPYEPwvlfFFPeK7xvheE8ngfyw7pLgL/Fifc5wp3g4ZzfBDqE5ftjYnTgN2H+b8DNYf5poEeYPwh4Js5rjANuiFluQ3T180H4HLKBZ4BTKt/jmG1/CUwJ81OI7tZuQvRMjuJQfgQwM2afm2JizQGap/tvoDFOuhKR2lRZndWT6LkjU8N/4ceG6VXgFaIvx8r/Gt9390VhfiFQYNG4X13cfQaAu29y99KwzTx3X+HRYHeLiL7MK1WOafYa0YONvnT3T4EyM2tN9OVZZGZLiIbd6MI3w5Ev9+hZFpWaEn1xjnD32dWc74fu/mKYv5doWAyAIy16Ot1rRM+H2Cdmn+nh5/eBVe4+P5zjeo+GAD80HAt3fxNYDuxJlHQry5cQJUmIHsy0N/CimS0iGm9p12riJezvRMltYHhfDiF6sBNEybQyxnuBQ8OVzQ+Bf4TXuJ3oAU1VHU00fEnl63wG9AGec/dPw/ndF85lW/7t7lvd/Q2++YyqehkYbWaXALu6+1cJHFdqWV3Xy0oj4e4vm1l7ov9ADRjr7rfHbmPRMxnKYooqiKmmqUbV7bPjrNtaZbutYbvfhHh6u/sWMysh+i8fYGOV1yknSmrHAf9XTSxVxwxyM2sGTCD6r/9DM/vfmNeI9zo7yoge3DRgO/e7C/gPsAn4h3/zDIuqnOiK4HN33z/pKKs/dqVmVdbFfn5xH33r7veb2VzgJOAxM/u9uz9TyzHKNuhKRFIi1JlnEQ1vPQs4o7Ku3sy6mNnO1e3r0VPjVpjZKWH73Mr2hh20E/BJSCBHUvN/7E40+F/P8J9uPN3M7JAw/2vgBb75MlwTzveXcfeMBuTrZGZ9AEJ7SDbwX6Jkh5ntCXQL2z4fXgMz25eoSgtgDvAjM+se1uWH/ar6kqh6Lzo595XASuCvRAmlUpOYmH8NvODRMzPeN7NfhdcwM9svzmvMBoZXLlg08u484PDQtpQFDOCbpLzazPayqCPDqdW8T9Weg5ntDrzn7uOJRrvtVd2OkjpKIlKbmodGzkVEVSKDPXrOwpNE9e4vhyqeh4j5MqjGb4HzQ9XTS8D3aiG++4DCEMMgonaBarl7BdGX3lFmdk6cTd4ieuDXMqL6/1vd/XOiZ0q8TpQ851dz7M3A6cBNZraY6Au48iqmSYhxOjDE3cuIHnzUIrzW5URXSYTquiHAA+G9epmourCqKcBt4fOpvNq7j6hKblnMdhuBvhZ1Xz4qvBZEie3MEOtS4j8i+u9Am9CAvhg40qNh20cSDXG+GFjo7pXDm48kah97iWg0421ZAlRY1BHhj8BpwOvh921fosfkSh3TKL4ijZSZ3Qy86u6TYso2uHuLGnYT+RYlEZFGyMwWEl11HBOudCrLlURkuyiJiIhI0tQmIiIiSVMSERGRpCmJiIhI0pREREQkaUoiIiKSNCURERFJ2v8DaUUOa5LjcCMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "n_genes = 170\n",
    "def count_targets(target_ids):\n",
    "    return np.array([np.sum(target_ids == i+1) for i in range(n_genes)])\n",
    "        \n",
    "# Load benchmark genes\n",
    "benchmark = pd.read_csv('benchmark/iss_benchmark.csv')\n",
    "    \n",
    "# remove barcodes near the border of the image\n",
    "fov = [100,100,9215-200,5120-200]\n",
    "results = results[ (results['x'] >= fov[0]) & (results['x'] <= fov[2])  & (results['y'] >= fov[1]) & (results['y'] <= fov[3])]\n",
    "benchmark = benchmark[ (benchmark['x'] >= fov[0]) & (benchmark['x'] <= fov[2])  & (benchmark['y'] >= fov[1]) & (benchmark['y'] <= fov[3])]\n",
    "\n",
    "\n",
    "# Filter detections such that less than n non targeted barcodes are detected\n",
    "n = 10\n",
    "threshold = ntargeted['quality'].nlargest(n).iloc[-1]\n",
    "result_counts = count_targets(results[results['quality'] > threshold]['target_id'])\n",
    "\n",
    "# Load the benchmark data and filter so that less than n non targeted barcodes are detected\n",
    "threshold = benchmark[benchmark['is_non_targeted'] == 1]['quality'].nlargest(n).iloc[-1]\n",
    "benchmark_counts = count_targets(benchmark[benchmark['quality'] > threshold]['target_id'])\n",
    "\n",
    "\n",
    "# Plot the count of the different detected barcodes\n",
    "tmp = np.vstack((result_counts, benchmark_counts)).T\n",
    "r = np.corrcoef(tmp[:, 1], tmp[:, 0])[0, 1]\n",
    "x = np.linspace(50, 25000)\n",
    "f, ax = plt.subplots(figsize=(6, 6))\n",
    "ax.scatter(tmp[:, 1], tmp[:, 0], 50, zorder=2)\n",
    "ax.plot(x, x, '-k', zorder=1)\n",
    "\n",
    "plt.xlabel('Benchmark barcode type counts')\n",
    "plt.ylabel('ISTDECO barcode type counts')\n",
    "plt.title(f'r = {r}')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# References\n",
    "[1] *Qian, Xiaoyan, et al. \"Probabilistic cell typing enables fine mapping of closely related cell types in situ.\" Nature methods 17.1 (2020): 101-106.*\n",
    "\n",
    "[2] Axel Andersson, Ferran Diego, Fred A. Hamprecht, & Carolina Wählby. (2021). In Situ Sequencing data used in \"ISTDECO: In Situ Transcriptomics Decoding by Deconvolution\" [Data set]. Zenodo. http://doi.org/10.5281/zenodo.4560687"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "istdecotest",
   "language": "python",
   "name": "istdecotest"
  },
  "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"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}