symm_expand_filter.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from cryosparc.tools import CryoSPARC\n",
    "import json\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import plotnine\n",
    "from plotnine import ggplot, aes\n",
    "\n",
    "with open('/u/rposert/instance-info.json', 'r') as f:\n",
    "    instance_info = json.load(f)\n",
    "\n",
    "cs = CryoSPARC(**instance_info)\n",
    "cs.test_connection()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Problem\n",
    "It is not rare to have a target which:\n",
    " * has some type of symmetry, and\n",
    " * at each symmetry-related copy of some position, either has or does not have a feature (say, a binding partner).\n",
    "\n",
    "Depending on the symmetry involved, such a situation can become unwieldy. For a simple C2 symmetric protein (as will be used in this demonstration), only four classes are necessary to capture the full range of possible particle states:\n",
    " 1. No binder\n",
    " 2. binder on A\n",
    " 3. binder on B\n",
    " 4. binder on A and B\n",
    "\n",
    "However, even with another common symmetry like, say, the C5 symmetry of Cys-loop receptors, the number of possible particle states is 2^5 = 32 --- and this is *assuming perfect classification*. As the number of states increases as 2^N, the odds of perfect classification go down as the symmetry number goes up.\n",
    "\n",
    "# Solution\n",
    "Rather than attempt to classify the entire protein at once, we can instead perform symmettry expansion and then classify the presence or absence of the feature of interest *in only one symmetry-related position*. Since the particles have been expanded, we're really detecting presence/absence in each position. Then, we can come back and count the number of positions which have our feature of interest.\n",
    "\n",
    "In higher symmetries, it may instead be useful to develop some method of describing the relationship of the regions of interest to each other. For instance, in a C4 symmetric target, it would be useful to split the particles into those that have\n",
    " 1. 0 features of interest\n",
    " 2. 3 features of interest\n",
    " 3. 4 features of interest\n",
    " 4. 2 features of interest on adjacent subunits\n",
    " 5. 2 features of interest on opposite subunits\n",
    "\n",
    "# Implementation\n",
    "After symmetry expanding and performing 3D classification, we can count the number of expanded subunits that contain a feature of interest.\n",
    "Then, using the source id, we can figure out how many of that feature the original particle had. We can then partition particles\n",
    "into subsets based on their composition with potentially more sensitivity than other methods which consider each symmetry-related position as\n",
    "distinct.\n",
    "\n",
    "# Example dataset\n",
    "Here we will be classifying particles from [EMPIAR-10425](https://www.ebi.ac.uk/empiar/EMPIAR-10425/). Specifically, we are trying to split particles into sets which have 0, 1, or 2 copies of the MlaB subunit, like so:\n",
    "\n",
    "![Maps showing subunits with 0, 1, and 2 MlaB subunits](attachments/sym-exp-filter_mla-count.png).\n",
    "\n",
    "We accomplish this by\n",
    " 1. Symmetry-expanding the particles in C2\n",
    " 2. Create a mask around a single MlaB binding site\n",
    " 3. Performing a 3D classification using that mask\n",
    " 4. For each particle image, count how many times it has or does not have MlaB\n",
    " 5. Export the particles in their appropriate sets\n",
    "\n",
    "![The proposed workflow](attachments/sym-exp-filter_workflow.png)\n",
    "\n",
    "# Counting the subunits\n",
    "Here, we are trying to find particles which have 0, 1, or 2 MlaB. The relevant job numbers are:\n",
    "* J12 - Non-Uniform Refinement with enforced C2 symmetry\n",
    "* J14 - Symmetry Expansion in C2 (rotate all potential MlaB binding sites into same position in space)\n",
    "* J17 - 3D Classification with a mask surrounding **only one** MlaB binding site. It's important to note that the out-of-mask position is *not informative* in this job.\n",
    "    * Class 0 here has no MlaB in mask\n",
    "    * Class 1 and 2 both have an MlaB mask"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "project_number = 'P310'\n",
    "class_job_number = 'J17'\n",
    "\n",
    "project = cs.find_project(project_number)\n",
    "class_job = project.find_job(class_job_number)\n",
    "class_zero = class_job.load_output('particles_class_0')\n",
    "class_one = class_job.load_output('particles_class_1')\n",
    "class_two = class_job.load_output('particles_class_2')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, we create pandas data frames to make life easier. We combine the two classes which had MlaB."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "df_zero = pd.DataFrame(class_zero.rows())\n",
    "df_one = pd.DataFrame(class_one.rows())\n",
    "df_two = pd.DataFrame(class_two.rows())\n",
    "df_with_mlab = pd.concat([df_one, df_two])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we count the number of times a given particle had or did not have MlaB. We use `sym_expand/src_uid` here since that is unique to the original particle, not the symmetry-expanded particles we're working with."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "no_counts = df_zero['sym_expand/src_uid'].value_counts()\n",
    "yes_counts = df_with_mlab['sym_expand/src_uid'].value_counts()\n",
    "\n",
    "# get yes and no counts in the same df\n",
    "count_by_src_uid = no_counts.to_frame().join(\n",
    "    yes_counts,\n",
    "    how = 'outer',\n",
    "    lsuffix = '_no',\n",
    "    rsuffix = '_yes'\n",
    ").fillna(0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we make sure that each original (un-expanded) particle appears exactly twice in our total counts."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABQAAAAPACAYAAABq3NR5AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAB7CAAAewgFu0HU+AABb40lEQVR4nOzdeZhnV0En7s9JpxNCFiMGSZqORkAGEsNEMIgrhIgg6BVkGCEiwggIsqiIiCiOIJARlEUWxQ1EbFxA8BpFnAGDqCP7D0kCighjkk7YIWQhdJLz++Pesr8UVdW13Kpv9an3fZ567nbOPSfpyuXy6XPuKbXWAAAAAABtOmLeHQAAAAAANo8AEAAAAAAaJgAEAAAAgIYJAAEAAACgYQJAAAAAAGiYABAAAAAAGiYABAAAAICGCQABAAAAoGECQAAAAABomAAQAAAAABomAAQAAACAhgkAAQAAAKBhAkAAAAAAaJgAEAAAAAAaJgAEAAAAgIYJAAEAAACgYQJAAAAAAGiYABAAAAAAGiYABAAAAICGCQABAAAAoGECQAAAAABomAAQAAAAABomAAQAAACAhgkAAQAAAKBhAkAAAAAAaJgAEAAAAAAaJgAEAAAAgIYJAAEAAACgYQJAAAAAAGiYABAAAAAAGiYABAAAAICGCQABAAAAoGECQAAAAABomAAQAAAAABomAAQAAACAhgkAAQAAAKBhAkAAAAAAaNiR8+4AwE63f//+vUkuHQ9P3bNnz2Xz7A9A6zx3AbaeZy/MlxGAAAAAANAwASAAAAAANEwACAAAAAANEwACAAAAQMMEgAAAAADQMAEgAAAAADRMAAgAAAAADRMAAgAAAEDDBIAAAAAA0DABIAAAAAA0TAAIAAAAAA0TAAIAAABAwwSAAAAAANAwASAAAAAANEwACAAAAAANEwACAAAAQMMEgAAAAADQMAEgAAAAADRMAAgAAAAADRMAAgAAAEDDBIAAAAAA0DABIAAAAAA0TAAIAAAAAA0TAAIAAABAwwSAAAAAANAwASAAAAAANEwACAAAAAANEwACAAAAQMMEgAAAAADQMAEgAAAAADRMAAgAAAAADRMAAgAAAEDDBIAAAAAA0DABIAAAAAA0TAAIAAAAAA0TAAIAAABAwwSAAAAAANCwI+fdAQCA8847r867D+xYl867A+ws+/btK/PuAwA7jxGAAAAAANAwASAAAAAANEwACAAAAAANEwACAAAAQMMEgAAAAADQMAEgAAAAADRMAAgAAAAADRMAAgAAAEDDBIAAAAAA0DABIAAAAAA07Mh5d2ArdF23N8n9k5yb5KwkpyS5McllSS5M8uK+7y86xD0ekOSxSe6c5Ngklyb5iyTn933/yUPUvWWSpybpkuxNck2S9yR5Wd/3b1hF/+fWNgAAAACHt+ZHAHZdd2qS/0jy4gwh4GlJbkiyK8ntkzw6yXu7rnvCCvd4SZI/S3KvJCcmOZDk65M8KclFXdfdcYW6ZyS5aCx7u7HuieO9Xt913YsO0f+5tQ0AAADA4a/5ADBD0FeS/E2SH0pyct/3x2UYSXd2krdlGAn5613X3Wdx5a7rHpXkcUlqkqclOb7v+xMyjCR8f5JbJem7rjtqibpHJ+mTfPVY9qyx7vHjvWqSJ3Zd94ilOj7PtgEAAABow04IAD+T5M5939+77/t9fd9/LEn6vr+x7/t3ZRgN9/6x7FNmK47B2jPHwxf3fX9+3/fXjfXfl+R+GabU3i7JI5do+9FJbjOWud9YJ33fX9f3/flJXjKWe1bXdbu3S9sAAAAAtKP5ALDv+8/1ff/eFa5fn+TV4+FdFl0+N8nJGUbLPW+Jupcmec14+NAlbr9wbt9YdrHnjvfek+ScbdQ2AAAAAI1oPgBcpWvH7a5F5+85bi/u+/6yZeq+adzereu6YxdOdl13XIYpxrNlvsR4z0vGw3O3UdsAAAAANEIAOLjHuH3/ovOnj9uLV6i7cK0kucPM+TuO51Zb//RF5+fZNgAAAACN2PEBYNd1Zyd5wHj4u4sunzJu969wi8uXKL94fzX1T1l0fp5tAwAAANCIHR0Adl13iwzf0TsiyTuSvGJRkePG7bVZ3uy145eou9r6xy86P8+2AQAAAGjEkfPuwLx0XXdMktcnuW2STyZ5SN/3N86hK+XQRbZ32/v37987xX1gBzt5dn///pUG7gIAhzPvzuxg3nlhInv27FlurYhl7cgAsOu6o5P8WZLvTPLZJN/d9/2/L1H06nF78xVud+zM/ueXqLtQ/6pl6i/c+/OLzs+z7bVaapVhYH3eOe8OAACbyrszeOeFjVrzgK4dNwW467qjkrw2yX0yBGP37vv+vcsUv2Lc7lnhlrPXrlhmfzX1r1h0fp5tAwAAANCIHTUCsOu63Un+NMn3JrkmyX37vn/HClUuSXLfJGesUGbhWk3ywZnzHxjPlbHMB7O0hfqXbKO21+rUDdaHne7kHPxb0LOTXDnHvsC8GBED7BTendmpvPPCHO2YAHAM//4kSZfkuiTf2/f9Pxyi2luSPDnJGV3X7e37fqk51vcet2/v+/6ahZN931/ddd07k9w1w2jD1y3Rp71JTp9pa7u0vSbrmXsOHLTo+ydX+m8KANrlf+fZqbzzwnztiCnAXdcdmWG13/sn+UKSru/7C1dR9c1JPpZhJN3PLHHfU5M8ZDx89RL1943bh4xlF3vKeO8r8uUh3DzbBgAAAKARzQeAXdftyhCQPTDJ9Uke0Pf9/1lN3b7vv5jkl8bDJ3Rd99Rx9eB0XXenJBdkWIjjw0l+Z4lbvDzJR8cyF4x10nXdMV3XPTXJ48dyT+/7/sB2aRsAAACAdpRa67z7sKm6rvvOJG8dD7+Y5DOHqHJ23/df8h2irut+M8mPjYc3JLk2yQnj8ceT3LPv+4uXaf/MDCPsThpPXZVh9d2F6dcv7fv+8UvVnXfbwNbYv3//3hz8/tmppkOwE5133nltv5AAjPbt27fmlRuhBd55Yb6aHwGYL/1nPCrJrQ7xs2vxDfq+f0yGEYRvzhCiHZ1h5N0Lkpy5XAA31n1/hsU2XjjWOXq8x5uTPPBQAdw82wYAAADg8Nf8CECA7c7fhoIRgMDOYQQgO5V3XpivnTACEAAAAAB2LAEgAAAAADRMAAgAAAAADRMAAgAAAEDDBIAAAAAA0DABIAAAAAA0TAAIAAAAAA0TAAIAAABAwwSAAAAAANAwASAAAAAANEwACAAAAAANEwACAAAAQMMEgAAAAADQMAEgAAAAADRMAAgAAAAADRMAAgAAAEDDBIAAAAAA0DABIAAAAAA0TAAIAAAAAA0TAAIAAABAwwSAAAAAANAwASAAAAAANEwACAAAAAANEwACAAAAQMMEgAAAAADQMAEgAAAAADRMAAgAAAAADRMAAgAAAEDDBIAAAAAA0DABIAAAAAA0TAAIAAAAAA0TAAIAAABAwwSAAAAAANAwASAAAAAANEwACAAAAAANEwACAAAAQMMEgAAAAADQMAEgAAAAADRMAAgAAAAADRMAAgAAAEDDBIAAAAAA0DABIAAAAAA0TAAIAAAAAA0TAAIAAABAwwSAAAAAANAwASAAAAAANEwACAAAAAANEwACAAAAQMMEgAAAAADQMAEgAAAAADRMAAgAAAAADRMAAgAAAEDDBIAAAAAA0DABIAAAAAA0TAAIAAAAAA0TAAIAAABAwwSAAAAAANAwASAAAAAANEwACAAAAAANEwACAAAAQMMEgAAAAADQMAEgAAAAADRMAAgAAAAADRMAAgAAAEDDBIAAAAAA0DABIAAAAAA0TAAIAAAAAA0TAAIAAABAwwSAAAAAANAwASAAAAAANEwACAAAAAANEwACAAAAQMMEgAAAAADQMAEgAAAAADRMAAgAAAAADRMAAgAAAEDDBIAAAAAA0DABIAAAAAA0TAAIAAAAAA0TAAIAAABAwwSAAAAAANCwI+fdAQ5/Bw4ceFKSJ827H3C4Oumkk3bddNNNSZIjjjjinQcOHLhxzl0CADbJgQMHLpt3H2AevPPCdHbv3r13rXUEgEzhhCS3nncn4HBVSsmuXbsWDk+eZ18AgE3nvZkdyTsvzJcAkClcleTyeXcCDle11l033XTTyUlyxBFHXFlK8beh7ET+DzGwU3hvZkfyzgvzJQBkw3bv3v38JM+fdz/gcLV///69SS4dD8/es2ePqUHsRHXeHQDYCuuZtgUt8M4L82UREAAAAABomAAQAAAAABomAAQAAACAhgkAAQAAAKBhAkAAAAAAaJgAEAAAAAAaJgAEAAAAgIYJAAEAAACgYQJAAAAAAGiYABAAAAAAGiYABAAAAICGCQABAAAAoGECQAAAAABomAAQAAAAABomAAQAAACAhgkAAQAAAKBhAkAAAAAAaJgAEAAAAAAaJgAEAAAAgIYJAAEAAACgYQJAAAAAAGiYABAAAAAAGiYABAAAAICGCQABAAAAoGECQAAAAABomAAQAAAAABomAAQAAACAhgkAAQAAAKBhAkAAAAAAaJgAEAAAAAAaJgAEAAAAgIYJAAEAAACgYQJAAAAAAGiYABAAAAAAGiYABAAAAICGCQABAAAAoGECQAAAAABomAAQAAAAABomAAQAAACAhgkAAQAAAKBhAkAAAAAAaJgAEAAAAAAaJgAEAAAAgIYJAAEAAACgYQJAAAAAAGiYABAAAAAAGiYABAAAAICGCQABAAAAoGECQAAAAABomAAQAAAAABomAAQAAACAhgkAAQAAAKBhAkAAAAAAaJgAEAAAAAAaJgAEAAAAgIYJAAEAAACgYQJAAAAAAGiYABAAAAAAGiYABAAAAICGCQABAAAAoGECQAAAAABomAAQAAAAABomAAQAAACAhgkAAQAAAKBhAkAAAAAAaJgAEAAAAAAaJgAEAAAAgIYJAAEAAACgYQJAAAAAAGiYABAAAAAAGiYABAAAAICGCQABAAAAoGECQAAAAABomAAQAAAAABomAAQAAACAhgkAAQAAAKBhR867A1uh67oTk9w9yV2S3HncnjxefkDf929Yoe5Hk3ztIZp4ad/3j1/hHrdM8tQkXZK9Sa5J8p4kL1up7Zn6D0jy2LHvxya5NMlfJDm/7/tPHqLuhtoGAAAA4PC2U0YA3j/JG5I8Pcn9cjD8W4urknxsmZ+rlqvUdd0ZSS5K8qQkt0tyIMmJSe6V5PVd171opUa7rntJkj8by5841v/68X4XdV13x81qGwAAAIDD304JAJPkyiRvTPLsJA9cR/2f6Pv+5GV+nrZUha7rjk7SJ/nqJO9Pclbf9yckOT7J05LUJE/suu4Ry9R/VJLHjeWeluT4sf5Z4/1ulaTvuu6oqdsGAAAAoA07JQB8dd/3p/R9f9++73+h7/s/26J2H53kNhmm3d6v7/v3JUnf99f1fX9+kpeM5Z7Vdd3u2YpjqPfM8fDFfd+f3/f9dWP992UYyXhNhpF9j5yybQAAAADasSMCwL7vb5hT0w8dt/v6vr90ievPzTASb0+ScxZdOzfDVOWa5HmLK473e82idqZqGwAAAIBG7IgAcB66rjsuydnj4ZuWKtP3/WVJLhkPz110+Z7j9uKx3FIW7nu3ruuOnbBtAAAAABqxI1YBnsiTu657TpKTknwuyT8neW2SV/R9/4Ulyt8xSRn3L17hvhcnOSPJ6YvOnz5zfaW6Gdu5Q5J3T9Q2AAAAAI0wAnD1zsiwgu51GULAeyZ5WZJ3dF33NUuUP2Vmf/8K9718ifKzx6upu7j+RtsGAAAAoBECwEN7Q5L/luSWfd/fvO/7r8jw3bynJ7k+yZlJ/mqJlXiPm9m/doX7L1w7fpn6q6m7uP5G2wYAAACgEaYAH0Lf9z+5xLkrMqye+89J/jzD6MCHJ/mtdTRRDl1k00zS9v79+/dOcR/YwU6e3d+/f6WBuwDA4cy7MzuYd16YyJ49e5ZbK2JZAsAN6Pu+77rubUm+I8n35UsDwKtn9m+e5KplbnPzcfv5ReevXnR9KcfO7M/W32jba7XUKsPA+rxz3h0AADaVd2fwzgsbteYBXaYAb9zbx+1tFp2/YmZ/zwr1F65dsej8FYuur1R3cf2Ntg0AAABAI4wAnE5ddPyB8VzJMEX4g8vUO2PcXrLo/CVJ7jtzfaW6ddH9N9r2Wp26wfqw052cg38LenaSK+fYF5gXI2KAncK7MzuVd16YIwHgxn3zuP3o7Mm+76/uuu6dSe6a5D5JXre4Ytd1e5OcPh6+ZdHltyR5cpIzuq7b2/f9UvO77z1u3973/TUTtr0m65l7Dhy06PsnV/pvCgDa5X/n2am888J8mQK8gq7rVpxT3XXd/TJ8/y9JLliiyL5x+5Cu65b6m76nZBild0W+PIR7c5KPjdd/Zom2T03ykPHw1RO3DQAAAEAjdswIwK7rTlrm0vGLrn2u7/sD4/6vd11Xk7w2ybv6vr92vNetkvyPJL84lvtAklcsce+XJ/nJJKcluaDruh/u+/6fu647JslPJHn8WO7pM20mSfq+/2LXdb+U5DeSPKHruiuSvKjv++u6rrtTkj/IsAjIh5P8zpRtAwAAANCOUuviT9e1aQzyVuOcvu8vHOu8MsmPjOdrks+N+yfOlH9fku/v+/7/LdPumRlG2C2EjFdlWH13IXx9ad/3j1+q7lj/N5P82Hh4Q5Jrk5wwHn88yT37vr94M9oGtsb+/fv35uD3z041HYKd6LzzztsZLyTAjrdv3741r9wILfDOC/NlCvDKfjPJ85L8fZLLkhyV5Jgk+5P8ZZKHJ7nrcuFfkvR9//4Mi228MMNovaMzBHFvTvLAQwVwfd8/JskDx/JXjfU/nOQFSc5cLvybom0AAAAADn87ZgQgwHblb0PBCEBg5zACkJ3KOy/MlxGAAAAAANAwASAAAAAANEwACAAAAAANEwACAAAAQMMEgAAAAADQMAEgAAAAADRMAAgAAAAADRMAAgAAAEDDBIAAAAAA0DABIAAAAAA0TAAIAAAAAA0TAAIAAABAwwSAAAAAANCwSQPAUsotxp+jlrl+Rinlz0spny6lfL6U8tZSyr2n7AMAAAAAcNBkAWAp5e5JPpHk40nutMT12yf5xyTfm+TEJMcm+fYkf1VKedhU/QAAAAAADppyBOC9kpQkH6q1vmuJ689PcvxY5qYknx33S5KXlFJOnrAvAAAAAECmDQC/I0lN8ubFF0ope5N8z3j9wiS3rLXeIslDx3PHJnnkhH0BAAAAADJtALgwgu+iJa7dN8NIvyT58VrrZ5Kk1rovyf8Zr33XhH0BAAAAADJtAPjV4/aTS1y7+7i9pNb6wUXX3jhu7zBhXwAAAACATBsAHjtuj17i2rdlmOr7t0tcu2LcnjhhXwAAAACATBsAfm7cfsliHqWU05J8zXj4f5eotzA1uE7YFwAAAAAg0waA/5IhzLvfovMPnNl/2xL1Thm3S00dBgAAAAA2YMoA8G/G7d1LKY8tpdyslPJtSX42w+i+/6/WetkS9c4at/82YV8AAAAAgEwbAL48yefH/ZckuSbJ3yU5aTz3gsUVSiklyb0zBIT/OGFfAAAAAIBMGADWWj+W5MEZgr8y85Mkr6q1vnqJavfKwdWDl1ogBAAAAADYgCOnvFmt9Y2llP+S5Lwkt8sQBv5NrfVvlqlytyRvTXJTkr+fsi8AAAAAwMQBYJLUWq9I8murLPvMJM+cug8AAAAAwGDKbwACAAAAANuMABAAAAAAGjb5FOAkKaWcmORRGVb4PT3JLZIcWWs9clG5eyY5OcknV/hOIAAAAACwTpMHgKWUH0zy8iTHL5wat3WJ4mcmeUGS60opp9Rar5q6PwAAAACwk006BbiU8rAk+5KckCH4uyLJv65Q5RVJvpjkZkm+b8q+AAAAAAATBoCllD1JfjND8HdZknvVWvcm+dnl6owj/t46Ht5zqr4AAAAAAIMpRwA+PsNIvmuTnFtrffMq670jQ2j4XyfsCwAAAACQaQPA787wnb8/rLV+aA31PjxuT5uwLwAAAABApg0Av27c/sMa631+3B6/YikAAAAAYM2mXAX42HF7zTrrfWHCvsCGnHfeeUutWg1b4dJ5d4CdY9++fWXefQAAADbflCMAPzVuT15jvTuO209M2BcAAAAAINMGgJeM27Wu5vuADN8OfNeEfQEAAAAAMm0A+MYMq/l+bynlG1ZToZTyhCS3Hw//csK+AAAAAACZNgD8nSSfzvBdwQtKKXdarmAp5chSys8meX6G0X+XJXnNhH0BAAAAADLhIiC11qtKKY9N8kdJTk3yrlLK/87BVX4XRvx9Q5LvzfCtwJLkxiT/o9Z6w1R9AQAAAAAGU64CnFrrn5ZSTkzy60mOTnKfhUvj9oUzxUuS65M8stb65in7AQAAAAAMppwCnCSptf52krsmeUOG4K8s8ZMkFyT55lrrH07dBwAAAABgMOkIwAW11vcn+YFSylck+bYkpyU5McnVGb7399Za66c2o20AAAAA4KBNCQAX1Fo/l+SvNrMNAAAAAGB5k08BBgAAAAC2DwEgAAAAADRMAAgAAAAADVvzNwBLKTduRkeS1Frrpn6TEAAAAAB2mvUEbmXyXgAAAAAAm2I9AeDfJalTdwQAAAAAmN6aA8Ba6z02oR8AAAAAwCawCAgAAAAANEwACAAAAAANEwACAAAAQMPWswjIkkopRyd5SoZVgt9Ya33nKurcNcl9ktyU5Pxa641T9QcAAAAAmDAATHL/JM9IciDJb6+yzqVJfiHJriQXJ3n9hP0BAAAAgB1vyinA9xu3b621XrGaCmO5v80wavD7JuwLAAAAAJBpA8C7JKlJ3rbGen83bs+esC8AAAAAQKYNAL9m3H54jfU+sqg+AAAAADCRKQPAo8btgTXWWyh/zIR9AQAAAAAybQD46XG7Z431Fsp/bsK+AAAAAACZNgBcmPr73Wust1D+IyuWAgAAAADWbMoA8C0ZVvO9dynlbqupUEr51iT3ybB4yJsn7AsAAAAAkGkDwN/O8D2/kuQNpZRvXqnwGBK+bix/41gfAAAAAJjQkVPdqNZ6aSnleUmeluSWSf6+lPLnSd6Q5JIkVyc5LsnpSe6f5PuT7Mow+u8FtdZ/n6ovAAAAAMBgsgBw9PQkt03ygxnCvQeMP0sp4/a1SZ46cT8AAAAAgEw7BTh18JAkP5XkExlCvuV+Ppnkp2qtP1hrrVP2AwAAAAAYTD0CMElSa31RKeW3knxPku9IcmqSE5JcleSyJH+X5I211us2o30AAAAAYLApAWCSjOHen40/AAAAAMAcTDoFGAAAAADYXgSAAAAAANAwASAAAAAANGzN3wAspbxl3K211nOXOL9eX3I/AAAAAGDj1rMIyD2S1DWcX42ygboAAAAAwDLWuwrwcoFd2UBfAAAAAICJrTkArLUu+d3A5c4DAAAAAPMjtAMAAACAhgkAAQAAAKBh6/0G4JcppTxs3H1LrfWyNdTbk+S7kqTW+qqp+gMAAAAATBgAJnllhoVBHpBk1QFgkjPHujclEQAehg4cOPCkJE+adz8AWJsDBw6s5X+vAZiAZy871UknnbTrpptuSpIcccQR7zxw4MCNc+4SHLZ27969d611pgwAN8oKwoevE5Lcet6dAGDNPLsBtp5nLztSKSW7du1aODx5nn2BnWg7BIALwV+day/YiKuSXD7vTkzMixmwE2ynZ7fnLrBTbKdnL2yZWuuum2666eQkOeKII64spRgBCFtoOwSAtxy3V8+1F6zb7t27n5/k+fPux8QE0kDz1jN1YBN57gI7wjZ79sKW2b9//94kl46HZ+/Zs8d0eNhCc10FuJRydJKFxUM+Ms++AAAAAECL1jUCsJTyI0l+ZJnLzyql/OShbpHk2CR3GLc1yf9ZT18AAAAAgOWtdwrwaUnukS+frlOSnLGG+yx8/+/KJL+2zr4AAAAAAMvY6DcAl1q5dzWr+dYM3/z7SIaRf79Wa71yg30BAAAAABZZVwBYa31GkmfMniul3JQh2HtArbWfoG8AAAAAwAZNvQjIakb/AQAAAABbZKNTgP9TrXWuKwoDAAAAAF9ustCulHJjKeWGUsorp7onAAAAALAxU47auzHDFOC3TnhPAAAAAGADpgwAPzZur5nwngAAAADABkwZAF40bm874T0BAAAAgA2YMgB8TYYpwA8upVgNGAAAAAC2gSkDwD9I8rYkZyZ5oRAQAAAAAOZvsgCw1lqTfH+SNyZ5fJK3l1IeXkq5bSnl6KnaAQAAAABW78ipblRKuXH2MMldkvzuzPVD3aLWWifrDwAAAAAwYQCYIfRb6RgAAAAA2GJTBoB/l6ROeD8AAAAAYIMmCwBrrfeY6l4AAAAAwDSmXAUYAAAAANhmBIAAAAAA0DABIAAAAAA0TAAIAAAAAA2bchXg/1RKuXmS709ytyR7k5yQZNchqtVa67mb0R8AAAAA2KkmDwBLKY9J8pwkX7GWaknq1H0BAAAAgJ1u0gCwlPILSZ6RIdA7lIXAbzVlAQAAAIB1mOwbgKWUO2QI/5LkX5PcM8kx43FNcv8kxyU5M8lTk1w5XntFkpvVWg81RRgAAAAAWKMpRwA+JsNovmuTfHet9T+SpJSDA/xqrdcmuTjJxaWU30ry50kenuTYJA+esC8AAAAAQKZdBfjuGUb6/elC+LeSWutnM4wK/HSSB5VSugn7AgAAAABk2gDwa8btPy1z/ajFJ2qtn0ny+xlGDv7whH0BAAAAADJtAHj8uP3EovPXLbq+2HvH7TdN2BcAAAAAINMGgNeM28Uj/T43br8mS1v4DuGtJuwLAAAAAJBpA8CPjtvFQd6/ZJji+23L1Puv4/aLE/YFAAAAAMi0AeD7MgR9Zy46/3fj9pxSyl1mL5RSbpPkkRkWD/nghH0BAAAAADJtAHjhuL3novOvSnLD2NZbSinPLaU8upTy3CTvSnLcWO41E/YFAAAAAMjB7+9N4S+S3Jjka0sp31pr/cckqbV+uJTynCS/mCHs++kl6r4nyW9M2BcAAAAAIBOOAKy1firJ7ZOcnuQDi679UpKfzbBQSJn5SZI/SXKvWqtvAAIAAADAxKYcAZha60dWuPa8UsqvJ/mWJCdnCAPfVWu9Yso+AAAAAAAHTRYAllKOSfIVST5ba/3CUmVqrdfn4LcCAQAAAIBNtqEpwKWUE0opzy6lfCjJ1UkuT3JNKeVDpZRnlVK+YpJeAgAAAADrsu4AsJRymwyLdzw1yW3ypd/2u02Sn0vynlLKbSfoJwAAAACwDusKAEspRyT54wxBX3JwQY/MHJckX5fkj8fyAAAAAMAWW28w1yW5S5Ka5LNJHpvk1kmOSrI3yY8n+fRY9huT3H8jnQQAAAAA1me9AeCDxu0Xk5xTa315rfWKWusNtdb9tdbfTHJOkuvHcv9tox0FAAAAANZuvQHgwui/P661/vNSBWqt70/yRxmmAt95ne0AAAAAABuw3gDw5HH7d4co97ZF5QEAAACALbTeAPD4cfupQ5Rb+A7gcetsBwAAAADYgPUGgAur/tZDlFu4vniVYAAAAABgC6w3AAQAAAAADgMCQAAAAABo2JEbrP+sUspPrnD9pIWdUspbDnGvWms9d4P9AQAAAABmbDQAPGMVZRa+A3j3FcqUHPp7ggAAAADAGm0kALSwBwAAAABsc+sNAM+ZtBcAAAAAwKZYVwBYa33r1B0BAAAAAKZnFWAAAAAAaJgAEAAAAAAaJgAEAAAAgIYJAAEAAACgYQJAAAAAAGiYABAAAAAAGiYABAAAAICGCQABAAAAoGECQAAAAABomAAQAAAAABq25gCwlPJ7489Zm9AfAAAAAGBC6xkB+PAkP5Lka2ZPllL+vZTy4VLKd03RMQAAAABg446c8F6nJalJbj7hPQEAAACADVjPCMDrx+0JU3YEAAAAAJjeegLAK8ftXafsCAAAAAAwvfVMAf6HJF+b5DGllKOSvDsHRwUmyT1LKSeupzO11letpx4AAAAAsLT1BIAvSvKDSXYledT4s6AkecI6+1KTCAABAAAAYEJrngJca31nkocm+WSGwG/hZ0HZwA8AAAAAMKF1rQJca/3jUsrrk5yd5NQkRyd5RYZRfC9J8p7JeggAAAAArNu6AsAkqbV+McP3AJMkpZRXjLtvrrX2G+0YAAAAALBx61kFeCWm8QIAAADANrLuEYCL1VqnDhMBAAAAgA2aLADczrquOzHJ3ZPcJcmdx+3J4+UH9H3/hlXc4wFJHjvWPzbJpUn+Isn5fd9/8hB1b5nkqUm6JHuTXJPhO4kv2+5tAwAAAHB427RRe6WUXaWUby2lPKmU8qullJeP2yeN57dyxOD9k7whydOT3C8Hw79V6bruJUn+LMm9kpyY5ECSr0/ypCQXdV13xxXqnpHkorHs7ca6J473en3XdS/arm0DAAAAcPjblBCulPKEJB9J8rYkz0vyU0keOW6fN57/SCnlcZvR/jKuTPLGJM9O8sDVVuq67lFJHpdhheOnJTm+7/sTkpyV5P1JbpWk77ruqCXqHp2kT/LVY9mzxrrHj/eqSZ7Ydd0jtlvbAAAAALRh0gCwlHJ0KeWNSV6Y5NYZFgVZ7ufUJL9eSvmrUsqXBVgTe3Xf96f0fX/fvu9/oe/7P1tNpTFYe+Z4+OK+78/v+/66JOn7/n0ZRhNek2F03SOXuMWjk9xmLHO/sU76vr+u7/vzk7xkLPesrut2b5e2AQAAAGjH1CMAfy/JvXMw5Pv7JL+Q5EFJvmfc/nyGEYAZy9x7rLdp+r6/YZ1Vz80wXbhmGLm4+L6XJnnNePjQJeovnNs3ll3sueO99yQ5Zxu1DQAAAEAjJgsASynfluQhGUKly5Pco9b6nbXW59RaX1drfdO4Pb/Wevck35lhMYuS5CFj/e3mnuP24r7vL1umzJvG7d26rjt24WTXdcclOXtRmS8x3vOS8fDcbdQ2AAAAAI2YcgTgw8ftdUnOqbX+3UqFa61/n+S7klw7ntqO36I7fdxevEKZhWslyR1mzt9xPLfa+qcvOj/PtgEAAABoxJQB4HdmGP33ylrrv62mQq31Q0lemSGs+o4J+zKVU8bt/hXKXL5E+cX7q6l/yqLz82wbAAAAgEZMGQAuhEj/tMZ6C+W3Ywh13Li9doUys9eOX6Luausfv+j8PNsGAAAAoBFHTnivXeP2xjXWWyi/a8VS7SqHLrK9296/f//eKe4DwNby/AbYep697GAnz+7v37/SZDVgJXv27FlurYhlTRkAfizJ1yb5xhxcnXY1vnGm/nZz9bi9+Qpljp3Z//wSdRfqX7VM/YV7f37R+Xm2vVZLrTIMwPbn+Q2w9Tx7IXnnvDsAh7k1D+iacgrw/x078IhSyletpkIp5aQMi3/UrH3q8Fa4YtzuWaHM7LUrltlfTf0rFp2fZ9sAAAAANGLKEYCvSfKQJLdI8qZSyg/UWv9jucKllK9N8rokX5UhANw3YV+mckmS+yY5Y4UyC9dqkg/OnP/AeK6MZT6YpS3Uv2Qbtb1Wp26w/nbkb2aBnWA7Pb89d4GdYjs9e2ErnZyDI//OTnLlHPsCO85kAWCt9YJSyoVJ7pFhWu8lpZQ/TvKmJP+a5JoMU1Zvn+TeSX4wyTEZgqoLa60XTNWXCb0lyZOTnNF13d6+75eaY33vcfv2vu+vWTjZ9/3VXde9M8ldk9wnQ9j5Jbqu25vk9Jm2tkvba7KeuecAzJ/nN8DW8+xlp1r0zb8r/bcAW2vKKcBJ8t8yjCYrGb4v9/AMIwPfnWEU2rvH44eP18tY/kET92Mqb87wbcKS5GcWX+y67tQMox6T5NVL1F8Y1fiQsexiTxnvfUW+PISbZ9sAAAAANGLSALDW+ukMo85emuQLGQKm5X6uS/LrSe461ttUXdedNPszc+n4Rdd2L1zo+/6LSX5pPHxC13VP7brumPF+d0pyQYZRjR9O8jtLNPvyJB8dy1ww1knXdcd0XffUJI8fyz297/sDsxXn2TYAAAAA7Si11s258bAQyH2T3C3JKUmOz7Da7BUZFvz4y60I/hZ0Xbfaf9Bz+r6/cFHd30zyY+PhDUmuTXLCePzxJPfs+/7iZdo9M8MIu4XQ8aoMox8Xpl+/tO/7xy9Vd95t72TnnXfe5vyHAbCN7Nu3b82rh20Wz11gp9hOz17YSvv379+bg9/8PdUUYNhaU08B/k+11k/VWv+g1vq4WusP1FrvNW4fN57fsvBvo/q+f0ySB2aYlntVkqMzjLx7QZIzlwvgxrrvz7DYxgvHOkeP93hzkgceKoCbZ9sAAAAAHP42bQQgHM6MRAF2gu00CsVzF9gpttOzF7aSEYAwX5s2AhAAAAAAmD8BIAAAAAA0TAAIAAAAAA0TAAIAAABAwwSAAAAAANAwASAAAAAANEwACAAAAAANEwACAAAAQMMEgAAAAADQsCOnulEp5S3j7ttqrf9zqvsCAAAAAOs3WQCY5O7j9rUT3hMAAAAA2IAppwB/ctx+fMJ7AgAAAAAbMGUA+O/j9lYT3hMAAAAA2IApA8A3JClJvnfCewIAAAAAGzBlAPgbSS5N8t2llAdPeF8AAAAAYJ0mCwBrrVcl+f4klyV5VSnl10opp011fwAAAABg7SZbBbiU8pZx93NJTk3yk0l+spSyP0MoeN0hblFrredO1R8AAAAAYMIAMMk9ktRxf2FbkuwZf1ZSZuoAAAAAABOZMgBMhiBvNecAAAAAgC0wWQBYa51yQREAAAAAYAJCOwAAAABomAAQAAAAABomAAQAAACAhk29CMh/KqWckuTcJKcnuUWS3bXWH92s9gAAAACALzd5AFhK+aokL0jy4CS7Fk4nqUl+dFHZ3xjPXVprve3UfQEAAACAnW7SKcCllK9L8t4kP5QhXCzjz3JeOpY7rZTynVP2BQAAAACYMAAspexK0ifZmyH0+8Mk907y+OXq1FovSvKB8fDeU/UFAAAAABhMOQX4YUnOyDDV9wm11pclSSnl5oeod2GSOyb55gn7AgAAAABk2inADxy3f7cQ/q3SReP29hP2BQAAAADItAHgN2YY/ff6Ndb71Lj9qgn7AgAAAABk2gBwIcC7fI31FhYJmXRBEgAAAABg2tDtmnF77Brr3XrcfnrCvgAAAAAAmTYAvHTcnrXGevcctx+arisAAAAAQDJtAHhhhum8Dy6lHLOaCqWUs5LcJ8O3A98yYV8AAAAAgEwbAL4iQ5B3qyS/VUopKxUupdwmyWvHPhxI8nsT9gUAAAAAyIQBYK31fUlemWEU4HlJ/m8p5UFJTl4oU0o5vpTyLaWU5yV5X5LbZAgNX1hrvWyqvgAAAAAAgyMnvt9jM4R6d09ydpI/Gs/XcfvZmbILIwTfmORpE/cDAAAAAMi0U4BTa/1iku9K8itJrs8Q8i38ZNHx9Un+V5Ku1nrTlP0AAAAAAAZTjwBMrfXGJD9XSnl+kh9M8p1JTktyYpKrk1yWYcGQP6q17p+6fQAAAADgoMkDwAW11k8kecn4AwAAAADMwaRTgAEAAACA7UUACAAAAAAN27QpwKWUOyZ5YJK7JtmT5Pgkn0+yP8nbk7yu1vrBzWofAAAAANiEALCUcuskv5HkfssU+cbx2jNLKX+R5HG11sun7gcAAAAAMPEU4FLKmUnemyHgK6v4+b4k7ymlfMOU/QAAAAAABpMFgKWU45P8VZKTMoR7+5M8Pcm3jOeOGbd3S/ILSS4by90yyV+WUo6bqi8AAAAAwGDKEYA/keTWSWqSP09yh1rrs2utb6+1frrWev24fUet9TlJ7pDkDWPdvWN9AAAAAGBCUwaAPzBu/z3Jg2utV69UuNZ6bZKHJPnweOqBE/YFAAAAAMi0AeBtM4z+e2Wt9frVVBjLvSLDVODbTtgXAAAAACATLwIy+tc1lv+3TegDAAAAAJBpA8D/GLcnrrHeVyyqDwAAAABMZMoA8IIMU3nvu8Z698swdfgvJuwLAAAAAJBpA8AXJ/l0kq6U8qDVVCil/LckXZJPjfUBAAAAgAlNFgDWWvdnWMn3M0n+sJTy3FLKyUuVLaWcXEr5lST7MoR/D6y1XjFVXwAAAACAwZFrrVBK+b1DFHlfknOS/HSSnyqlXJLkQ0muSXJsktslOSMHw8f3JXlEKeXhtdYfXWt/AAAAAIDlrTkATPLwDN/sW8nC9V1JvmH8mVVmypwzc14ACAAAAAATWk8AmAwB3kbLLj5/qFARAAAAAFij9QSAXzd5LwAAAACATbHmALDW+v82oyMAAAAAwPQmWwUYAAAAANh+BIAAAAAA0DABIAAAAAA0bL2rAK+olHJkkrsmOSPJVya52Wrq1VqfuRn9AQAAAICdatIAsJSyO8lTkzwxyS3WcQsBIAAAAABMaLIAsJSyK8kFSb5r4dQab1Gn6gsAAAAAMJhyBODjk9xr3P9Ckj9M8rYkVya5fsJ2AAAAAIBVmjIAfNi4/UyS76i1XjLhvQEAAACAdZhyFeDbZ5jG+xLhHwAAAABsD1MGgAfG7UUT3hMAAAAA2IApA8CPjtvjJrwnAAAAALABUwaAr8uw8u89JrwnAAAAALABUwaAv5lhxd8Hl1K+acL7AgAAAADrNFkAWGv9VJLvT3JdkjeVUs4rpZSp7g8AAAAArN2RU96s1vrOUsp/TfKGJH+Q5PmllHcl+WSSmw5dvf7olP0BAAAAgJ1u0gCwlHJckp9Pcvp46quTfM8abiEABAAAAIAJTRYAllKOSfK/k9x14dSi7aHUqfoCAAAAAAymHAH42CTfPO5/PMmLk7wtw8Ig10/YDgAAAACwSlMGgA8dt5cnObvWeuWE9wYAAAAA1mGyVYCT3DbDNN6XCv8AAAAAYHuYMgA8MG7/bcJ7AgAAAAAbMGUA+O/j9isnvCcAAAAAsAFTfgPwT5N8U5LvSfLbE96Xbe7AgQNPSvKkefcDgLU5cODAZfPuA8BO49nLTnXSSSftuummm5IkRxxxxDsPHDhw45y7BIet3bt3711rnSkDwJcleWSS7y+lfG+t9YIJ7832dkKSW8+7EwCsmWc3wNbz7GVHKqVk165dC4cnz7MvsBNNFgDWWq8ppdwvyQVJXltKeVaSl9RaPztVG2xbV2VY/bklXsyAnWA7Pbs9d4GdYjs9e2HL1Fp33XTTTScnyRFHHHFlKcUIQNhCkwWApZS3jLvXJDkqyTOS/M9Syr8k+WSSmw5xi1prPXeq/rB1du/e/fwkz593PyZW590BgM22nqkDm8hzF9gRttmzF7bM/v379ya5dDw8e8+ePabDwxaacgrwPXLw5X1huyvJHVdRt8SLPwAAAABMbsoAMBmCvNWcAwAAAAC2wJTfADxiqnsBAAAAANMQ2gEAAABAwwSAAAAAANAwASAAAAAANEwACAAAAAANm2wRkFLKL270HrXWZ07RFwAAAABgMFkAmOSXktQN3kMACAAAAAATmjIATJKygbobDQ8BAAAAgEWmDADPWUWZXUlOSvKtSX44yYlJ9iX5nQn7AQAAAACMJgsAa61vXUPxPymlPCvJa5Ocl+R9tdZfnaovAAAAAMBgbqsA11o/maRL8qkkzyml3HlefQEAAACAVs0tAEySWutVSV6RYSTi4+fZFwAAAABo0VwDwNH/N25X8w1BAAAAAGANtkMAuLD678lz7QUAAAAANGg7BIB3G7fXzLUXAAAAANCguQaApZRvTvKoDKMA/3mefQEAAACAFh051Y1KKd+5yqJHJbl1ku9K8t+T7M4QAL5yqr4AAAAAAIPJAsAkF+bg9/xWq4zbv6i1vmrCvgAAAAAAmX4KcFnjz6eT/GKS/zZxPwAAAACATDsC8BmrLHd9ks8kuTjJO2qtX5ywDwAAAADAjMkCwFrragNAAAAAAGCLzHUVYAAAAABgcwkAAQAAAKBhAkAAAAAAaJgAEAAAAAAatq5FQEop/z51R5LUWuttN+G+AAAAALBjrXcV4NOS1A22XRYdb/R+AAAAAMAi6w0Aky8P8NarTngvAAAAAGDGer8BuHuDPycm+eUk16y/6wAAAADAoaxrBGCt9cb11CulHJHkUUl+McnJC6eTfCJDIAgAAAAATGgjU4DXpJTyA0mek+TrF05lGAH4/CTPq7VevVV9AQAAAICdYtMDwFLKdyR5bpK7LpxKckOS30nyjFrrxza7DwAAAACwU21aAFhKOSPJ/0py34VT4/Z1SZ5Wa/3QZrUNAAAAAAwmDwBLKXuTPDPJD2dYZGQh+Htrkp+ttb5j6jYBAAAAgKVNFgCWUr4iydOSPCHJ0TkY/F2U5Km11r+aqi0AAAAAYHU2HACWUo5K8sQkP5fkxBwM/i7NsNrvq2qtdaPtAAAAAABrt6EAsJTyIxmm++7NweDvM0nOT/LiWuv1G+seAAAAALAR6woASyn3zbDAxxk5GPx9IcmLk5xfa/3sJL0DAAAAADZkvSMAL0hSM4R/NyZ5VZJfrLVePlXHAAAAAICN2+g3AGuSLya5V5J7lVIOUXzle9Vav3aD/QEAAAAAZkyxCvDNktx6g/coGcJEAAAAAGBCGwkANzTcDwAAAADYfOsKAGutR0zdEQAAAABgeoI8AAAAAGiYABAAAAAAGiYABAAAAICGCQABAAAAoGECQAAAAABomAAQAAAAABomAAQAAACAhgkAAQAAAKBhAkAAAAAAaJgAEAAAAAAaJgAEAAAAgIYJAAEAAACgYQJAAAAAAGiYABAAAAAAGiYABAAAAICGCQABAAAAoGECQAAAAABomAAQAAAAABomAAQAAACAhgkAAQAAAKBhAkAAAAAAaJgAEAAAAAAaJgAEAAAAgIYJAAEAAACgYQJAAAAAAGiYABAAAAAAGiYABAAAAICGCQABAAAAoGECQAAAAABomAAQAAAAABomAAQAAACAhh057w4cDrque3iSV6yi6C37vv/kMvd4QJLHJrlzkmOTXJrkL5Kcv1ydmbq3TPLUJF2SvUmuSfKeJC/r+/4Nq+j/utsGAAAA4PBmBODa3JTkYyv83LRUpa7rXpLkz5LcK8mJSQ4k+fokT0pyUdd1d1yuwa7rzkhy0Vj2dmPdE8d7vb7ruhet1OGNtA0AAADA4U8AuDaX9n1/8go/n15coeu6RyV5XJKa5GlJju/7/oQkZyV5f5JbJem7rjtqibpHJ+mTfPVY9qyx7vHjvWqSJ3Zd94ilOruRtgEAAABogwBwE43B2jPHwxf3fX9+3/fXJUnf9+9Lcr8M03lvl+SRS9zi0UluM5a531gnfd9f1/f9+UleMpZ7Vtd1uyduGwAAAIAGCAA317lJTs4wAu95iy/2fX9pkteMhw9dov7CuX1j2cWeO957T5JzJm4bAAAAgAYIADfXPcftxX3fX7ZMmTeN27t1XXfswsmu645LcvaiMl9ivOcl4+G5U7UNAAAAQDusArw2t+y67j1J/st4fHmSCzNMsX3/EuVPH7cXr3DPhWslyR2SvHs8vuN4bjX1z5hpa4q2AQAAAGiEEYBrc/MMC2hcn2R3htV0H5XkvV3XPXmJ8qeM2/0r3PPyJcov3l9N/VMWnd9I2wAAAAA0QgC4OvuT/FKSM5PcrO/7WyQ5NsM027cn2ZXkeV3Xnbeo3nHj9toV7j177fgl6q62/vGLzm+kbQAAAAAaYQrwKvR9/zdJ/mbRuQNJ/rbruu9M8tYkd0vyK13X/VHf9zdtYffKoYtsrv379++ddx8AWDvPb4Ct59nLDnby7P7+/StNVgNWsmfPnuXWeliWAHCD+r7/Ytd1P5/kzUn2JvnGHPyW3tXj9uYr3GJ28Y3Pz+xfPbN/8yRXLVN/4d6fX3R+I22v1VIrFAOw/Xl+A2w9z15I3jnvDsBhbs2DwUwBnsbbZ/ZvM7N/xbjds0Ld2WtXLLO/mvpXLDq/kbYBAAAAaIQRgNOrM/uXJLlvhlV6l7NwrSb54Mz5D4znyljmg1naQv1LFp3fSNtrdeoG6m5X/mYW2Am20/PbcxfYKbbTsxe20sk5OPLv7CRXzrEvsOMIAKfxzTP7H53Zf0uSJyc5o+u6vX3fLzVH+97j9u1931+zcLLv+6u7rntnkrsmuU+S1y2u2HXd3iSnz7Q1a91tr9V65p4DMH+e3wBbz7OXnWrRN/+u9N8CbC1TgA+h67oV51V3Xbc7yS+Ph/uTvGfm8puTfCzDKL6fWaLuqUkeMh6+eonb7xu3DxnLLvaU8d5X5MsDwI22DQAAAEADBICH9rVd172967pHdV132sLJruuO7Lru7kkuTPKt4+mnzq4A3Pf9F5P80nj4hK7rntp13TFj/TsluSDDQhwfTvI7S7T98gwjCo9NcsFYJ13XHdN13VOTPH4s9/RxVeL/NEHbAAAAADSg1FoPXWoHG0O/j8ycuj7DCrvHJzlqPHdDkp/r+/5Xl7nHbyb5sZmy1yY5YTz+eJJ79n1/8TJ1z8wwuu+k8dRVGVb2XZi+/dK+7x+/VN2Ntr2TnXfeef7DAJq3b9++Na8etlk8d4GdYjs9e2Er7d+/f28OfvP3VFOAYWsZAXhoH0vyxCR/lGGhjKuTfEWSLyR5X5IXJTlzufAvSfq+f0ySB2aYlntVkqMzjLx7wVh32QCu7/v3Z1is44VjnaPHe7w5yQNXCv822jYAAAAAhz8jAGEJRqIAO8F2GoXiuQvsFNvp2QtbyQhAmC8jAAEAAACgYQJAAAAAAGiYABAAAAAAGiYABAAAAICGCQABAAAAoGECQAAAAABomAAQAAAAABomAAQAAACAhgkAAQAAAKBhAkAAAAAAaJgAEAAAAAAaJgAEAAAAgIYJAAEAAACgYQJAAAAAAGiYABAAAAAAGiYABAAAAICGHTnvDgAAALD1zjvvvDrvPrBjXTrvDrCz7Nu3r8y7D/NmBCAAAAAANEwACAAAAAANEwACAAAAQMMEgAAAAADQMAEgAAAAADRMAAgAAAAADRMAAgAAAEDDBIAAAAAA0DABIAAAAAA0TAAIAAAAAA0TAAIAAABAwwSAAAAAANAwASAAAAAANEwACAAAAAANEwACAAAAQMMEgAAAAADQMAEgAAAAADRMAAgAAAAADRMAAgAAAEDDBIAAAAAA0DABIAAAAAA0TAAIAAAAAA0TAAIAAABAwwSAAAAAANAwASAAAAAANEwACAAAAAANEwACAAAAQMMEgAAAAADQMAEgAAAAADRMAAgAAAAADRMAAgAAAEDDBIAAAAAA0DABIAAAAAA0TAAIAAAAAA0TAAIAAABAwwSAAAAAANAwASAAAAAANEwACAAAAAANEwACAAAAQMMEgAAAAADQMAEgAAAAADRMAAgAAAAADRMAAgAAAEDDBIAAAAAA0DABIAAAAAA0TAAIAAAAAA0TAAIAAABAwwSAAAAAANAwASAAAAAANEwACAAAAAANEwACAAAAQMMEgAAAAADQMAEgAAAAADRMAAgAAAAADRMAAgAAAEDDBIAAAAAA0DABIAAAAAA0TAAIAAAAAA0TAAIAAABAwwSAAAAAANAwASAAAAAANEwACAAAAAANEwACAAAAQMMEgAAAAADQsCPn3QEOfwcOHHhSkifNux8ArM2BAwcum3cfAHYaz16Ardfas3f37t1711pHAMgUTkhy63l3AoA18+wG2HqevQBbb8c/ewWATOGqJJfPuxMT2/EPB2BH2E7Pbs9dYKfw7AXYetvp2TsXAkA2bPfu3c9P8vx592Nidd4dANhs65k6sIk8d4EdwbMXYOtts2fvXFgEBAAAAAAaJgAEAAAAgIYJAAEAAACgYQJAAAAAAGiYABAAAAAAGiYABAAAAICGCQABAAAAoGECQAAAAABomAAQAAAAABomAAQAAACAhgkAAQAAAKBhAkAAAAAAaJgAEAAAAAAaJgAEAAAAgIYJAAEAAACgYQJAAAAAAGiYABAAAAAAGiYABAAAAICGCQABAAAAoGECQAAAAABomAAQAAAAABomAAQAAACAhgkAAQAAAKBhAkAAAAAAaJgAEAAAAAAaJgAEAAAAgIYJAAEAAACgYQJAAAAAAGiYABAAAAAAGiYABAAAAICGCQABAAAAoGECQAAAAABomAAQAAAAABomAAQAAACAhgkAAQAAAKBhAkAAAAAAaJgAEAAAAAAaJgAEAAAAgIYJAAEAAACgYQJAAAAAAGiYABAAAAAAGiYABAAAAICGCQABAAAAoGECQAAAAABomAAQAAAAABomAAQAAACAhgkAAQAAAKBhAkAAAAAAaJgAEAAAAAAaJgAEAAAAgIYJAAEAAACgYQJAAAAAAGiYABAAAAAAGiYABAAAAICGCQABAAAAoGECQAAAAABomAAQAAAAABomAAQAAACAhgkAAQAAAKBhAkAAAAAAaJgAEAAAAAAaJgAEAAAAgIYJAAEAAACgYQJAAAAAAGiYABAAAAAAGiYABAAAAICGCQABAAAAoGECQAAAAABomAAQAAAAABomAAQAAACAhgkAAQAAAKBhAkAAAAAAaJgAEAAAAAAaJgAEAAAAgIYJAAEAAACgYQJAAAAAAGjYkfPuAJuv67pbJnlqki7J3iTXJHlPkpf1ff+GOXYNAAAAgE1mBGDjuq47I8lFSZ6U5HZJDiQ5Mcm9kry+67oXza93AAAAAGw2AWDDuq47Okmf5KuTvD/JWX3fn5Dk+CRPS1KTPLHrukfMr5cAAAAAbCYBYNseneQ2Gab83q/v+/clSd/31/V9f36Sl4zlntV13e459REAAACATSQAbNtDx+2+vu8vXeL6czOMAtyT5Jwt6xUAAAAAW0YA2Kiu645LcvZ4+KalyvR9f1mSS8bDc7eiXwAAAABsLQFgu+6YpIz7F69QbuHa6ZvbHQAAAADmQQDYrlNm9vevUO7yJcoDAAAA0AgBYLuOm9m/doVyC9eO38S+AAAAADAnR867A8xdOXSRle3fv3/vFB0BYGt5fgNsPc9egK3X2rN3z549l621jgCwXVfP7N88yVXLlLv5uP38BtpaaoVhALY/z2+ArefZC7D1Wnv2rnkwlwCwXVfM7O/J8gHgniXK73i/+qu/Ou8uAOwonrsAW8+zF2DnEAC26wNJaoZU+IwkH1ym3Bnj9pINtHXqBuoCyclJ3jnun53kyjn2BWAn8NwF2HqevTBHAsBG9X1/ddd170xy1yT3SfK6xWW6rtub5PTx8C3rbWs9c8+Bg/bv/5KFuq/03xTA5vLcBdh6nr0wX1YBbtu+cfuQruuWGqX3lAwjBK/IBgJAAAAAALYvAWDbXp7ko0mOTXJB13V3SpKu647puu6pSR4/lnt63/cH5tNFAAAAADaTALBhfd9/IUmX5JNJ7pTkfV3XfS7DgiDnZxj999K+7393fr0EAAAAYDMJABvX9/37Myz08cIkH05ydIYA8M1JHtj3/eOXrw0AAADA4c4iIDtA3/cfT/JT4w8AAAAAO4gRgAAAAADQMAEgAAAAADRMAAgAAAAADRMAAgAAAEDDBIAAAAAA0DABIAAAAAA0TAAIAAAAAA0TAAIAAABAwwSAAAAAANAwASAAAAAANEwACAAAAAANEwACAAAAQMMEgAAAAADQMAEgAAAAADRMAAgAAAAADRMAAgAAAEDDBIAAAAAA0DABIAAAAAA0TAAIAAAAAA0TAAIAAABAwwSAAAAAANAwASAAAAAANEwACAAAAAANEwACAAAAQMMEgAAAAADQMAEgAAAAADRMAAgAAAAADRMAAgAAAEDDBIAAAAAA0LBSa513HwAAAACATWIEIAAAAAA0TAAIAAAAAA0TAAIAAABAwwSAAAAAANAwASAAAAAANEwACAAAAAANEwACAAAAQMMEgAAAAADQMAEgAAAAADRMAAgAAAAADRMAAgAAAEDDBIAAAAAA0DABIAAAAAA0TAAIAAAAAA07ct4dADhcdV13yyRPTdIl2ZvkmiTvSfKyvu/fsElt7k1y/yTnJjkrySlJbkxyWZILk7y47/uLNqNtgHnquu7EJHdPcpckdx63J4+XH7BZz915tw0wT9vh3XMe79zQIiMAAdah67ozklyU5ElJbpfkQJITk9wryeu7rnvRJrR5apL/SPLiDC9ipyW5IcmuJLdP8ugk7+267glTtw2wDdw/yRuSPD3J/XIwgGu9bYC52A7vnvN454ZWCQAB1qjruqOT9Em+Osn7k5zV9/0JSY5P8rQkNckTu657xMRN70pSkvxNkh9KcnLf98clOTbJ2UnelmFk9693XXefidsG2A6uTPLGJM9O8sAd1DbAPMz13XOO79zQJFOAAdbu0Uluk2H6wf36vr80Sfq+vy7J+V3XnZLkCUme1XXdq/u+PzBRu59Jcue+7987e7Lv+xuTvKvrunsleWeSM5M8JclfT9QuwHbw6r7vXzl7ouu6ndA2wLzM+91zXu/c0CQjAAHW7qHjdt/Ci8giz83wN5J7kpwzVaN9339u8QvYouvXJ3n1eHiXqdoF2A76vr9hJ7YNMC/b4N1zLu/c0CoBIMAadF13XIYpD0nypqXK9H1/WZJLxsNzt6JfM64dt7u2uF0AAHaeTXn3PAzeueGwIwAEWJs7ZvgWSpJcvEK5hWunb253vsw9xu37t7hdAAB2nnuM26nfPbf7OzccdgSAAGtzysz+/hXKXb5E+U3Vdd3ZSR4wHv7uVrULAMDOs8nvntv2nRsOVwJAgLU5bmb/2mVLHbx2/Cb25T91XXeLJK/J8Fx/R5JXbEW7AADsPFvw7rkt37nhcCYABNgc5dBFptF13TFJXp/ktkk+meQh4+psAAAwqW327rll79xwuBMAAqzN1TP7N1+h3MK1z29iX9J13dFJ/izJdyb5bJLv7vv+3zezTQAAdqYtfPfcVu/c0AIBIMDaXDGzv2eFcgvXrlihzIZ0XXdUktcmuU+Sq5Lcu+/7925WewAA7Fxb/O65bd65oRUCQIC1+UCSOu6fsUK5hWuXbEYnuq7bneRPk3xvkmuS3Lfv+3dsRlsAAOxsc3j33Bbv3NASASDAGvR9f3WSd46H91mqTNd1e5OcPh6+Zeo+jC9gf5KkS3Jdku/t+/4fpm4HAADm8e65Hd65oTUCQIC12zduH9J13alLXH9Khg8SX5GJX0a6rjsyw4pr90/yhSRd3/cXTtkGAAAkc3/3nNs7N7So1FoPXQqA/9R13c0yTEs4Lck/J/nhvu//eVwR7SeSPCfDy8gj+77/3SXqvzLJjyT5f33fn7aGdncl+cMkP5jk+iT37/v+rzf0DwNwGOm67qRFpz4xbh+W5I0z5z/X9/2BRXVfmXU8e6doG+BwNMW750aevRt95wa+1JHz7gDA4abv+y90Xddl+JvGOyV5X9d1V2VYhWzhufrSTXgR+bYML2DJ8LLzyqEbyzq77/tLJ+4DwDx9Ypnzr1p0fE6SCxtqG2Ae5vruOcd3bmiSKcAA69D3/fszfHT4hUk+nOToDKuhvTnJA/u+f/wK1RdWK3vnCmWWMvvMPirJrQ7xs2uN9wdo2XqfvQA71RTvnht69m7wnRuYYQowwBYaP6L82STHJDmr7/t/nm+PANrn2Quw9Tx7YXsxAhBga901w7SF13sJAtgynr0AW8+zF7YRASDA1jpn3P7yXHsBsLN49gJsPc9e2EZMAQYAAACAhhkBCAAAAAANEwACAAAAQMMEgAAAAADQMAEgAAAAADRMAAgAAAAADRMAAgAAAEDDBIAAAAAA0DABIAAAAAA0TAAIAAAAAA0TAAIAAABAwwSAAAAAANAwASAAAAAANEwACAAAAAANEwACALBjlFKOLKVcVEqppZRfm3d/mE4p5U/HP9e/nXdfAGC7EQACABxCKeW4UsoPl1JeWUq5uJTy8VLKF0spnyulfLiU8rpSyk+VUk6Zd1/nqZRyVinll8afs+bdn2X8RJIzknw2ybOXKlBKuXAMkhZ+fnU1Ny6l/P1MndMm6/EmKaU8dqa/7y6l7F5D3RfP1H3NZvZzDZ6W5IYk9yilPGTenQGA7UQACACwjFLKrlLKU5J8JMmrkvxIktOT3DLJ7iQnJLlNkh9I8vwkl5VSXlNKue2cujxvZyX5n+PPWXPtyRJKKV+Z5BfGw5fUWj+9yqqPL6V8zSZ1a25qrb+R5E3j4Z1z8N/Nikop5yR53Hh4xcz+XNVaP5Tkj8bDZ5dSjppnfwBgOxEAAgAsoZRyiyR/neRXkpw0nv63JC/JEHj8YJL/keTnM4QoX8jwbvXgJKaWbk9PTnJikuuTvGgN9Y5O8sub0aFt4EczjIZMkqeVUu6yUuFSyvFJXpGkjKceuYYgdSv8yrj9ugz/fQIAEQACAHyZcSpkn+S7xlMfS/KgJLevtT6h1vqyWuuf1FpfUWt9Tq31PkluneQ5Sa6ZT69ZSSnluCQ/Ph7+aa31k6usuvDn+dBSyp2m79l81VovT/L48fDIJK8qpRy9QpXnJ/nacf+3a61/tZn9W6ta60VJ/nE8fFIpxf/fAYAIAAEAlvK/knzbuP8fSe5Wa31trbUuV6HW+ula688nOTvJxVvQR9bmhzKM/kuS319DvReM2yMy/F40p9b6h0leNx6enuRZS5UrpdwnySPHw48m+elN79z6LPz5fn2S755nRwBguxAAAgDMKKXcOgdHRNUkP1Rr/ehq69daPzAGgcvd/+hx8YW/LqXsL6VcX0r5dCnlvaWU55ZSvu4Q/Xv4zOILDz9E2dNmyr5ymTIfHa9/dDw+Ymzjb0spHyulfKGU8h+llFeXUv7rSn3KMDV0wSsWLaRRxzKL655QSvnpmfa+WEq5alxc5R9LKc8vpdxngu+5PXzcfjrJhWuo9wdJ/nnc/55Syj022I8kSSnlVqWUZ5RS/qmU8onxn/vKcQGSJ49TbbfSYzKMdE2GkXPfvqi/Jyb5nfHwpiQPr7V+flGZm5VSfqyUckEp5dLxd+dz46rLv15Kuf2hOjHR78Prxz4mw3c7AWDHO3LeHQAA2GYel2QhXPjrWuvfT3XjUsqdk/xZDk6hXHBUkq/MsHDGT5RSfr7WuqqVZ6dUSvmqDCPB7r7o0qkZRtA9uJTy8Frrqydq7y5JLkhy8qJLu5Mcn2GBlW9J8lMZRla+a53tnJLkm8fDf6i13rCG6jcl+bkkfzkeP7eU8s0rjQZdRX8eluSlSY5bdOlW48/dkzyllPKQWuub19vOWtRaP1lKeXSSP88wSOCVpZT/WmtdmAL94gzT3JPkhbXWt87WL6XcPckfzpRZcHSGVZfPSPLjpZSn11rPX6oPU/0+1Fo/UUq5OMmZSe5XStldaz2w4r8AAGicABAA4EvNThlcy1TRFZVSviHJW3Mw9PmXDKPL/i3JVyS5b5LvzxAGPq+UcnSt9dlTtb8KR+Zg+PdPSV6b5NIMweSDkpybZFeS3y6lvH1ccXXBW5I8IMk9kzxhPPfi8fySSik3T/KGHAx73p1h5NblGb6795VJ7pjknCRLjjxcg+/OwUUr/mmtlWutf1VKuTDJPTIETw9K8ifr6cg4anN2pOTfZvjn/liGoPUhSe6SYaXpN5ZS7l1r/dv1tLVWtdZ+HCn68CS3TfK8DKHd/ZM8dCz2gQwL3/ynUsr3ZAgOd2cYNft/MiyMc1mG3+dvSvKwDFOwn1NKyeIQcBN+H/4pQwB4fIbp/Beuog4ANEsACAAwGheK+MaZU/8w0X1LhtFRC+Hf7yd51KJRSb9VSnlAkj/OEKQ8o5Ty17XWd0/Rh1W49fjzU7XWFy669vJSykszLKJxsyRPzMGgL7XW/0jyH+M00QXvqbW+YYX27ptk77j/a7XWJy9XsJRyepJPrO4fY0l3m9l/+zrv8bMzdZ9dSnn9WkeVlVK+JsMq0skQlD221vryRWVekOT8JE/J8Hvw6lLK19dar11nv9fqJzIEuV+T5DGllLfl4HcQb0jysFrrF2b6e0qSV499/VySH6i1Lg5+/6CU8r8yrKp9pyS/PP77++BMmal/H/5vkkeN+3eLABCAHc43AAEADjo5B9+Prq+1XjbRfe+bIfhIkouSPHKp8KjW+vocXIBhV4YQaCv9wRLh34KfS3LduP89E7R1u5n931upYK31klrrRgLAM2f2P7CeG9Ra35FhVGQy9P1RKxRfzhOSHDvu/9bi8G9s56Za68/m4OjJPdnC79jVWq9K8ogMAWVJsi/DtOQkeXatdfG0259Jcotx/2FLhH8L970iyX9PcmOG3+2fWFRk6t+H2T/n5lZvBoC1EgACABz0VTP7n53wvg+c2f/VQ3yD7oVJFkZ7fV8pZfeE/TiU5y93YQyGFsKf25ZSbrbBtq6Z2b/LBu91KKfN7H96A/d5WoZRcEnyi+OI0bVY+D2oOfSKws9Zot6WGEO8lyw6/e4sWh14HNn6w+Phv9Ra+0Pc91+SvGM8vPeiy1P/PnxqZv+0Ce4HAIc1U4ABAJa27kUeljA7BfVNKzZa61WllH9M8l1JjsnwvbN1LX6xRtcmed8hylw+s39ikis30N7/ycFRZr9RSvn6JH9Ua71kA/dczsIItetmp6+uVa31Q6WU307y2Ayj4n46yTNWU7eU8tVJFlZ4/tdVrCz91gwjLo/JwQVMttLPZlj4ZeHf3eOWCK5PT3LSuH/l+K3AQ7lx3H5dKeVmM38eU/8+zAa9t1i2FADsEEYAAgAcNDtq6CsnvO8p4/bztdbVhGb/MrO/Z8J+rORTq1jZ9vqZ/Q2NAKy1fiAHR5Qdm+TpSS4upXyslPK6UspPlFL+y0bamHH0uL1qgns9I8nV4/6Tx2BvNU6Z2f/XQxUew7Z/Hw+PK6WcsPoublyt9bokn5859bElip02s3/3DIt2HOrn22fq/Gcwtwm/D7N/1sesoR4ANEkACABw0BVJbhr3jy6l7F2p8BocP26vWbHUQVfP7B+/bKlp3XToItOqtf5iku/Lly628tVJfiDDVOgPllL+vpSy0RFwC8HlhkO0WuvHcnCq9HFJfnGVVWf/HLfz78FanLjB+kfNHkz8+zD7Z33dsqUAYIcQAAIAjGqt1yR578ypb5vo1gsjqY5dsdRBs9+W+/yypQ5t1wbqbola6wW11m/PsADLgzKEa+/OwSnY35bkbaWU79pAMwsjO4+Z4NuFSfK8JB8f9x9dSrndSoVHs3+OW/17sFlmA8oX1lrLGn8+uviGE/4+zE773ch3HwGgCQJAAIAvNfuNvqlWX71i3B5fSjl5FeVvP7O/f9G12Wm4R2VlJx3i+rZRa/1YrfW1tdafrrV+U4bppQur7u5O8oIN3P6jM/sb/h5crfXqJL88Hu5O8uxVVLtiZv/2y5YalVJ2JbnNeHj1uAjLdjO7SvY3THnjCX4fZhf0+eiUfQOAw5EAEADgS70sB0O2+5RSphgF+E8z+4tXP/0SpZTjk3zreHhdvnxhjs/M7N/6EO1+6yGuT212GnHZyI1qrf+R5LwknxhPfUMp5cR13u6fZ/bvsJF+zXh5kg+P+w8qpXzTSoVrrR9P8pHx8PallNMOcf+75+C3696+3k5usv8vB1fL/o5SyqYFzuv4fbjjzP6hFrcBgOYJAAEAZtRaL0/ykvGwJNlXSvna1dYvpdyxlPKsRadfO7P/06WUI1e4xU/k4BTRvtZ6YNH1i2f2l50GOU51feyh+jux2Smhq53muqzxn3125eGV/r2tZDaAvduypdZg7NvPj4clyXNXUW3h96BkWGV3JT+3RL1tpdZ6Y5JXj4dHZ3UjITfS3lp+H2b/nP9p2VIAsEMIAAEAvtzP5eBCBF+T5J9KKQ8spSw7qq2Ucosx+Htnvnw65BtzcBTSmUl+q5Sye4l7dBlWP02SG7NEqFRrvTQHQ8BvLaU8aIn7HJ3k95N8/XL93SQfmdm/80oFSylPLKU8qJSy7DTmUsp3JLnTeHhZrfWT6+zX3+TgN+Q2uqDIrD9J8q5x/5wkZx2i/ItzcAGQHyulPGpxgTJ4Tg6Gu5cnedXGu7ppnpOD39h7dCnlV5b63V5QSjmmlPKIUsqDF52f+vfhW8bt55P84yH/KQCgcev9W1QAgGbVWg+MYdwfZwhiTs4wCutDpZS/TnJJhoUljk2yJ8l3JLlHkiUXmKi11lLKD2UYiXRckkck+ZZSyquS/HuGFUu/J8kDZqr9z1rre5bp4q/kYCi0r5TyPUnekuSGDFMfH5bhm2l/mOSH1viPvxHvT/KxJLdK8tBSyicz/DNfu1Cg1vrX4+6dM3xj8XOllDcleU+Gb8rdkGHl13OSfG8O/oX1c9bbqVrrx0op/5Dk25N8eyll9xIjK9dz31pK+dkkbx5PrTjqsdZ6aSnl8UlekWEU4G+VUn4wyeszTG29dYZprgvTiQ8keWit9dql7ldKeWUOfqfyGbXWX1r/P8361FqvGEPov8zw+/+UJD9USvnTDFOvP5/h38vXZvjnOjfJzXMw6F4w2e9DKeWrkpw+Hv7lFH/WAHC4EwACACyh1vrpUsp9kjw5yc9kWFTg67PyqLobk7wmXx5upNZ6cSnlHklelyEMuUOWDjG+mOTptdZlp5TWWv+glPKtSR6T4X3uEePPrJdkWEF1ywLAWusNpZSfT/I7GRZq+Jklii2Molz4XuBXJPnv489SDmQIt35jg937/QwB4C2S3C/JGzZ4vyRJrfUtY2C14rcdZ8q/spRSk7w0QzB27viz2CeTPKTWeuEU/dxM47+Db02yL8Pv9a2T/OQKVW5McuWic1P+PjwgB1fA/v0VygHAjiEABABYxviNs18ppbw0yQ9kCGq+KcOIpBMzjGz7RIbpvX+X5I9qrR9b4X7vLqX8lyQ/muT7M0xn/KoM00L/X5L/neRltdaPLHePmXs9dhyN+JixTyeMfXn7eI83r2KhicnVWn+3lPL/xn6dneHf1VIjIx+T5JUZRnZ9e5L/kuSWGd5Pr0ryoSQXJvndWuuHJujaazKMnLxFhhGSb5jgngt+Nsm9ssrP69Raf3/8s/vxJPdJctsMf36fTfKBJBck+c1a6+cPcaubz+x/YtlSW6DW+t5SyhkZwrfvz/ANvltlCDmvTnJphhGif5vkL2qtiwPAKX8fFkZF/luG6d8AsOOVWuuhSwEAwGGulPLLSX4hwyjLvbXWuYZmG1VKuSLD9PQrktxuuanCO0kp5fQc/Ebmj08wchQAmmAREAAAdopfS/KZJEdl+FbdYauUcmaG8C9Jnin8+08Lqyt/NMnvzrEfALCtCAABANgRaq2fTfLs8fBxpZSTVyi+3S2sEvxvGb65uOOVUm6bYRGVJHlarfWL8+wPAGwnAkAAAHaSF2WYInpMkqfNuS8bca9x+/Ra6w1z7cn2cX6Gbwa+tdb6mnl3BgC2E98ABAAAAICGGQEIAAAAAA0TAAIAAABAwwSAAAAAANAwASAAAAAANEwACAAAAAANEwACAAAAQMMEgAAAAADQMAEgAAAAADRMAAgAAAAADRMAAgAAAEDDBIAAAAAA0DABIAAAAMD/344dyAAAAAAM8re+x1cYwZgABAAAAIAxAQgAAAAAYwIQAAAAAMYEIAAAAACMCUAAAAAAGBOAAAAAADAmAAEAAABgTAACAAAAwJgABAAAAIAxAQgAAAAAYwIQAAAAAMYEIAAAAACMCUAAAAAAGBOAAAAAADAmAAEAAABgLFgzMWSCUUSEAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {
      "image/png": {
       "height": 480,
       "width": 640
      },
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<Figure Size: (640 x 480)>"
      ]
     },
     "execution_count": 48,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# sanity check - things should only end up at (2,0), (1,1), or (0,2)\n",
    "count_by_src_uid['coordinate'] = count_by_src_uid['count_no'].astype(int).astype(str) + ', ' + count_by_src_uid['count_yes'].astype(int).astype(str)\n",
    "(\n",
    "ggplot(count_by_src_uid, aes('coordinate'))\n",
    "+ plotnine.theme_minimal()\n",
    "+ plotnine.geom_bar()\n",
    "+ plotnine.labs(\n",
    "    x = 'Counts (No, Yes)',\n",
    "    y = 'Number of Particles'\n",
    ")\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We've successfully split particles into whether they seem to have zero, one, or two MlaB components. It seems that it is rare to have *no* MlaB!\n",
    "\n",
    "# Export\n",
    "\n",
    "To return these particle sets to CryoSPARC, we just need to create an external job for each subset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'J25'"
      ]
     },
     "execution_count": 62,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "double_mlab_uids = count_by_src_uid.loc[count_by_src_uid['coordinate'] == '0, 2'].reset_index()['sym_expand/src_uid']\n",
    "single_mlab_uids = count_by_src_uid.loc[count_by_src_uid['coordinate'] == '1, 1'].reset_index()['sym_expand/src_uid']\n",
    "zero_mlab_uids = count_by_src_uid.loc[count_by_src_uid['coordinate'] == '2, 0'].reset_index()['sym_expand/src_uid']\n",
    "\n",
    "unexpanded_job_number = 'J12'\n",
    "unexpanded_refinement = project.find_job(unexpanded_job_number)\n",
    "unexpanded_particles = unexpanded_refinement.load_output('particles')\n",
    "\n",
    "zero_mlab_particles = unexpanded_particles.query({\n",
    "    'uid': zero_mlab_uids\n",
    "})\n",
    "single_mlab_particles = unexpanded_particles.query({\n",
    "    'uid': single_mlab_uids\n",
    "})\n",
    "double_mlab_particles = unexpanded_particles.query({\n",
    "    'uid': double_mlab_uids\n",
    "})\n",
    "\n",
    "desired_workspace = 'W2'\n",
    "project.save_external_result(\n",
    "    desired_workspace,\n",
    "    zero_mlab_particles,\n",
    "    type = 'particle',\n",
    "    title = 'No MlaB',\n",
    "    passthrough = (unexpanded_job_number, 'particles')\n",
    ")\n",
    "project.save_external_result(\n",
    "    desired_workspace,\n",
    "    single_mlab_particles,\n",
    "    type = 'particle',\n",
    "    title = 'One MlaB',\n",
    "    passthrough = (unexpanded_job_number, 'particles')\n",
    ")\n",
    "project.save_external_result(\n",
    "    desired_workspace,\n",
    "    double_mlab_particles,\n",
    "    type = 'particle',\n",
    "    title = 'Two MlaB',\n",
    "    passthrough = (unexpanded_job_number, 'particles')\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The particles are now ready for further refinement or classification!"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "cryosparc-tools",
   "language": "python",
   "name": "cryosparc-tools"
  },
  "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.10.13"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}