{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fit models for strong enhancer vs. silencer\n",
    "This notebook fits various models to classify strong enhancers vs. silencers with 5-fold cross-validation. Fit 6-mer SVMs or logistic regression models on CRX occupancy, TF occupancies, and information content (entropy). Process the White 2013 data and use as a test set. Then fit logistic regression models using random TFs (100 times) to generate a background distribution. Also fit an information content logistic regression model for strong enhancer vs. inactive.\n",
    "\n",
    "Note that fitting the background distribution takes approximately 50 minutes on a single Intel Core i7 CPU with 3.9 GHz and 16 GB of RAM."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "import sys\n",
    "import pickle\n",
    "\n",
    "import matplotlib as mpl\n",
    "import matplotlib.patches as mpatches\n",
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from scipy import stats\n",
    "from sklearn.feature_selection import RFE, RFECV\n",
    "from sklearn.linear_model import LogisticRegression\n",
    "from sklearn.model_selection import StratifiedKFold\n",
    "from IPython.display import display\n",
    "\n",
    "sys.path.insert(0, \"utils\")\n",
    "from utils import fasta_seq_parse_manip, gkmsvm, modeling, plot_utils, predicted_occupancy, sequence_annotation_processing\n",
    "\n",
    "data_dir = os.path.join(\"Data\")\n",
    "seq_bins_dir = os.path.join(data_dir, \"ActivityBins\")\n",
    "models_dir = os.path.join(\"Models\")\n",
    "figures_dir = os.path.join(\"Figures\")\n",
    "all_seqs = fasta_seq_parse_manip.read_fasta(os.path.join(data_dir, \"library1And2.fasta\"))\n",
    "# Drop scrambled sequences\n",
    "all_seqs = all_seqs[~(all_seqs.index.str.contains(\"scr\"))]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "plot_utils.set_manuscript_params()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Read in data and prepare for modeling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "activity_df = pd.read_csv(os.path.join(data_dir, \"wildtypeMutantPolylinkerActivityAnnotated.txt\"), sep=\"\\t\", index_col=0)\n",
    "# Drop sequences where the WT group is NaN from analysis\n",
    "activity_df = activity_df[activity_df[\"group_name_WT\"].notna()]\n",
    "activity_df[\"group_name_WT\"] = sequence_annotation_processing.to_categorical(activity_df[\"group_name_WT\"])\n",
    "activity_df[\"group_name_MUT\"] = sequence_annotation_processing.to_categorical(activity_df[\"group_name_MUT\"])\n",
    "\n",
    "# Read in occupancies\n",
    "occupancy_df = pd.read_csv(os.path.join(data_dir, \"predictedOccupancies.txt\"), sep=\"\\t\", index_col=0)\n",
    "# Only keep WT with an annotated class\n",
    "wt_occupancy_df = occupancy_df[occupancy_df.index.str.contains(\"WT$\")].copy()\n",
    "wt_occupancy_df = sequence_annotation_processing.remove_mutations_from_seq_id(wt_occupancy_df)\n",
    "wt_occupancy_df = wt_occupancy_df.loc[activity_df.index]\n",
    "\n",
    "# Same idea for information contents\n",
    "entropy_df = pd.read_csv(os.path.join(data_dir, \"entropyDiversityOccupancy.txt\"), sep=\"\\t\", index_col=0)\n",
    "wt_entropy_df = entropy_df[entropy_df.index.str.contains(\"WT$\")].copy()\n",
    "wt_entropy_df = sequence_annotation_processing.remove_mutations_from_seq_id(wt_entropy_df)\n",
    "wt_entropy_df = wt_entropy_df.loc[activity_df.index]\n",
    "\n",
    "# Mask to pull out the silencers and strong enhancers\n",
    "silencer_modeling_mask = activity_df[\"group_name_WT\"].str.contains(\"Strong|Silencer\")\n",
    "# Mask to pull out the inactive seqs and the strong enhancers\n",
    "inactive_modeling_mask = activity_df[\"group_name_WT\"].str.contains(\"Strong|Inactive\")\n",
    "\n",
    "# Within the data to model, mask indicating which sequences are strong enhancers\n",
    "labels_with_silencer = activity_df.loc[silencer_modeling_mask, \"group_name_WT\"].str.contains(\"Strong\")\n",
    "labels_with_inactive = activity_df.loc[inactive_modeling_mask, \"group_name_WT\"].str.contains(\"Strong\")\n",
    "\n",
    "# Sequences file for the SVM\n",
    "positives_fasta = os.path.join(seq_bins_dir, \"strongEnhancer.fasta\")\n",
    "negatives_fasta = os.path.join(seq_bins_dir, \"silencer.fasta\")\n",
    "\n",
    "# Read in PWMs for the test set\n",
    "pwms = predicted_occupancy.read_pwm_files(os.path.join(\"Data\", \"Downloaded\", \"Pwm\", \"photoreceptorAndEnrichedMotifs.meme\"))\n",
    "ewms = pwms.apply(predicted_occupancy.ewm_from_letter_prob).apply(predicted_occupancy.ewm_to_dict)\n",
    "mu = 9"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fit k-mer SVM to the data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAASQAAAEkCAYAAACG+UzsAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzsnXd4lFX2xz8nhfQeSQghIBCkSBNYithdt+ii7tpdkGXF3tZeVlGXteBPXUTXVVFZe1k7qyKiKK4Vd0VAKRpqgExIm7TJZGbO7493Mk6GJCTwJpkk9/M88zDvfe9775kh+ebec+89R1QVg8FgCAciOtsAg8FgaMAIksFgCBuMIBkMhrDBCJLBYAgbjCAZDIawwQiSwWAIG4wgGQyGsMEIksFgCBuMIBkMhrAhqrMN6GgyMzN1wIABnW2GwdCj+Prrr3er6gF7q9fjBGnAgAGsXLmys80wGHoUIrKlNfXMlM1gMIQNRpAMBkPYYATJYDCEDUaQDAZD2GAEyWAwhA1hJ0gicomIrBSROhFZtJe6fxKRXSLiFJEnRCSmg8w0GAztQNgJErADmAs80VIlEfkFcD1wDNAfGAjc1u7WGQyGdiPs9iGp6qsAIjIeyG2h6jnA46q61l//L8CzWCJlMBhawOVy8cMPP9CWENYejweHw8GuXbsCryOPPJITTzzRNrvCTpDawAjgjaDrVUCWiGSoaklwRRE5DzgPIC8vr+MsNBg6EbfbzSeffMLKlSvZtWsX5eXlqCo7duxgxYoV1NbW7ncfdXV1RpD8JAIVQdcN75OARoKkqo8CjwKMHz/eZDUwdCl8Ph9ut5uqqipUlcrKysAIZevWrZSUlODz+QL1i4uLWb16NatWrWpRdLKzs4mOjm61HREREcTHx7NlyxYOOuggpk6dym9+85v9+myhdGVBqgKSg64b3ld2gi0Gg60UFRXx6quvsmTJElavXs3WrVvxeDxtbufAAw9k+PDhHHDAASQlJQEQHx/PUUcdxdixY4mJaf06kMPhYNq0adTU1FBfX88999zTpudbQ1cWpLXAaOAl//VooCh0umYwdBV8Ph9FRUXMmzePRx55pNHoRkRISkpCRIiJiSE9PZ2UlBQyMjKIj48nIuKn9an4+HgGDBjAoEGD6N+/P/n5+aSkpBATE0NkZOQ+2VZUVMRJJ53EunXrGD58OMuWLbNdjCAMBUlEorDsigQiRSQW8Khq6J+Hp4BFIvIs1srcn4FFHWmrwbA/eL1eqqqqcDqdFBYW8txzz/HCCy9QXFwMwMiRI5kyZQrDhg2jX79+JCcnExkZSUREBKmpqaSlpeHxeKirq2vknBYRkpOTSUhIIDk5uZFY7QtFRUUcffTRfPfddwwfPpwPP/yQ3r1771ebzRF2goQlLHOCrn8P3CYiTwDfAcNVdauqvisi84APgTjglZDnDIawwev1UllZSU1NDRUVFZSVlbFz505efPFF1qxZw6ZNm3C5XAAMHjyYiy++mIMOOojo6GhycnLIyckhJSUFEelQu0PF6IMPPmg3MYIwFCRVvRW4tZnbiSF17wPua2eTDIZ9pr6+nsLCQgoKCnC5XFRVVVFVVcWWLVv4+9//zu7duwN1R4wYwTnnnMPIkSOJiYlh5MiRpKen7/cIZ38oKyujpKSk3UdGDYSdIBkMXRlVxev1EhERQXl5OZ9//jlvvPEG7777Ljt27NjDMT1q1Chmz57N0KFDiYqKora2lry8PPLz89u0AtZeDB06lOXLl5Oent7uYgRGkAwGW3A6nezatYvCwsLASOiNN95g8eLFVFT8tDslMTGRjIwMMjIyOOaYYzj++OOpqanB5/ORmJjIIYccQnJycgs9tT8Oh4OlS5dy9tlnA5YodRRGkAyG/aC8vJwVK1awZMkSCgoK2Lp1K7t372b37t14vV7Ack7PnDmTSZMmNVqZahCuiRMnkpqa2qlTswYcDgdHHXUU3333HUBAlDoKI0gGwz7g8/l46aWXuP/++/n6668D4hPM5MmTmTlzJocccgg1NTU4nc5G91WViRMnkp6e3lFmt0iwGA0fPpyf//znHW6DESSDoRW4XK7AtGzx4sU89dRTrF+/HoDIyEgOP/xwJkyYQH5+Pn369CEjI4OoqCiqq6txOBykpqYyevRoIiIiEBFEhNjY2MBmxc4mVIw6woHdFEaQDIYWKC0t5ccff6SgoIAPP/yQd955h61btwKQmprKGWecwW9/+1vi4+NxuVx4vV68Xi9Op5OYmBiysrLIzs4mPT29w5fsW0u4iBEYQTIYmsTn8/Hiiy+yaNEiNm3aREFBQWBalp2dzRlnnMFRRx1FdHQ0Ho8HVaVfv37ExMQQFxdHYmIicXFxYStCwZx99tlhIUZgBMlgCOD1evnuu+94/fXXee6551i3bl3gXkREBFOmTGHatGmMHz8ej8fDkCFDSExMJD4+noSEhE60fP9YsGABl156Kc8++2ynihEYQTIYACvWzyOPPMLdd9/Ntm3bAGtKduqppzJx4kTy8/OJjY2lrKyMyMhIJk6cSGJi4l5aDV9cLhexsbGAtay/dOnSTrbIwgiSocficrl49913KS0tZcmSJbz0knVOOysrizPPPJOjjz46sDmxpqaGmpoaDjzwQAYMGNAuB0s7CofDwdFHH83555/PpZde2tnmNMIIkqFH8v3333Pqqaeydu3aQFlERATnn38+p556KvX19QwcOJDExESioqKIiIggISGhSwsRNHZgP/LII8yePTswUgoHjCAZehwvvvgis2fPprKyksTERFJSUkhJSeHyyy8nPz8fn8/H5MmTw2ZJ3i5CV9M++OCDsBIjMIJk6EG4XC6uueYaHnzwQQAOPfRQbrrpJuLj4wFrWhYfH8/BBx/cpZ3UTRFOS/stYQTJ0O1xOp0sWLCARx99lK1btxIZGckFF1zA8ccfT2xsbCCsx9ChQ+ndu3eXWKpvC11FjMAIkqEb4/P5ePLJJ5kzZw6FhYUA9OnThxtuuIHBgwczYMAABg8eHBan6tuTiooKysvLw16MwAiSoZvy448/MmPGDD799FPAyjYzY8YMJkyYQO/evRk2bBgpKSmdbGXHkJ+fz/Lly0lJSQlrMQIjSIZuQkVFBU6nkw0bNvDKK6+waNEiamtrSUlJYfbs2Rx77LFEREQwfPhwcnJyut20LBSHw8G7777LjBkzAEuUugJGkAxdGlXl6aef5qWXXqKgoID169cHUgIde+yxnH/++SQkJJCTk8PgwYOJi4vrZIvbn9AQIg2i1BUwgmTosqxcuZLLLruMzz77LFAWFRXFEUccwXHHHcfo0aPp27dvYD9RTyDUgf3LX/6ys01qE0aQDF0Kl8vFsmXLuO+++/jggw8ASE5O5qyzzmLkyJHk5+dTX19PdnZ24KxZT6ErraY1hxEkQ5dgx44dXHnllbzxxhuB7ByxsbFMmzaNU089lcTERFQVVWX06NE9wk8UTHcQIzCCZAhzfD4fixYt4qqrrqK8vBywHLSTJ0/m17/+NUOGDCEjI4O4uDji4+OJi4sLi1CwHc306dO7vBiBESRDGOHz+SgrK6O6upqdO3eyZs0aHnjgAb799lsAxo4dy+WXX86wYcPo3bs3WVlZPcJJ3RoefPBBLrnkEp5++ukuK0ZgBMkQJtTV1TFnzhw++ugjCgoKcDgcgXvp6enMmjWLiy66iNzc3G6/kbG11NbWBgQ5Pz+fJUuWdLJF+48RJEOH4/F4qK+vR1Xx+XwUFxfz+9//PrCJESAmJobBgwczYcIEzjjjDA477LDAmTPDTyFE/vCHP3DVVVd1tjm2YQTJ0GE4HA7WrFmD2+3G5XLx/PPPs3nzZjZt2kRxcTEpKSlcccUVDB8+nIyMDOrq6sjIyGDMmDFdPuyHnQQ7sJ988kkuuuiibjN1NYJk6BAKCwtZtWpVIP/Y9ddfz6pVqwL3BwwYwJ133klKSgqRkZGkpqZy4IEHkpaW1qNWy/ZGUyFEuosYgREkQwewdetWVq9ejdfr5a233uLll1/mxx9/JCsri8svv5zevXvTp08f4uLiGDNmDCkpKT1ypWxvdJel/ZYwgmRoV7Zv384333zDv//9bxYuXBg41tGvXz/uuece0tLSAiFAhg4daqZmzdATxAiMIBnaCa/Xy+bNm/n666+59957+fLLLxERjjzySA499FCOPfZYampqmDRpUqfnsu8KVFZW4nQ6u7UYQRgKkoikA48DxwG7gRtU9bkm6sUA84GTgWjgP8AFqlrYgeYagvB4PDidTmpra9m0aRNVVVW8/PLLfPnll6SmpjJ37lwmTZoEQHl5Obm5uUaMWsmgQYNYvnw5SUlJ3VaMIAwFCXgIcANZwBjg3yKySlXXhtS7HJgMjAIqgEeBBcBvO9BWg5+ysjK+/fZbampqiIiIID4+nqKiIl555RUiIyN5+OGHAyEwvF4vdXV1DBo0qJOtDm8cDgeLFy9m1qxZAD3i+worQRKRBOB3wMGqWgV8IiJvAtOB60OqHwgsUdUi/7MvAvd1pL0Giy1btrB27VqSk5MDf71ramq44447UFV+//vfBw69lpeXIyIMGzas28WttpPQECINotTdCStBAoYAHlXdEFS2CjiiibqPA/NFJAcoB84G3ml/Ew3BFBUVsWbNGg444ABqa2v58MMP+eijj1i2bBm1tbXk5ORw7rnnUlJSgogwfPhwsrOz6dWrV2ebHraEOrBPOOGEzjapwwg3QUoEnCFlFUBT+Wg2AtuAQsALrAYuaapRETkPOA+sUKaG/cfj8VBcXMz//vc/MjIyWLp0KX/5y1+oq6sL1Bk1ahTXXnstNTU1ZGdnc9BBB5lVtL3QU1bTmiPcBKkKCPVyJgOVTdR9CIgBMoBq4FqsEdLE0Iqq+iiWj4nx48erjfb2SBwOB9988w1er5fU1FQ++ugjbrnlFnw+H6NGjWLKlCn8/Oc/p3///vh8PkpKSsjPzzditBd6uhhB+AnSBiBKRPJVdaO/bDQQ6tAGy+F9k6qWAojIAuB2EclU1d0dY27PQ1VZt24dCQkJxMbG8p///IebbroJn8/H7NmzOf/88xvVdzqd9OvXr1vtJm4vzjnnnB4tRgBhtR1WVauBV7GEJUFEDgVOBJ5uovpXwAwRSRGRaOAiYIcRo/alITxIbGws3333Hddddx1er5fp06dz3nnn7VHf7XabaXIrefDBB/nlL3/ZY8UIwm+EBJawPAE4gBLgQlVdKyKHAe+oakNM0quBB7B8Sb2ANVh7kgztyKZNm4iLi+OHH37giiuuwOVycfzxx3PZZZchIvh8PpxOZ+A0f+/evbtdSmo7CQ4hMmjQIN55p2evy4SdIPmnYCc1Ub4Cy+ndcF2CtbJm6CBqamrYsWMHb7/9NgsXLsTj8TBx4kT+/Oc/IyI4nU7cbje5ublkZWXRq1cvEzKkBRpCiEyfPp3rrruus80JC9pNkESkt6o69l7T0BXw+XysX7+e++67j+XLlwPw29/+liuuuILo6Gjcbjf19fVMnTrV7C9qBcEO7KeffprLLrvM+NmwWZBEJAVrt/SpgA9IEJHfAONVdY6dfRk6Dp/Px9q1a/nrX//K8uXLSUhIYN68eUycaC1oqiqlpaWMGzfOiFEr6O4hRPYHu53aDwMurA2Obn/ZF8CZNvdj6ACcTic//vgjn376KQ888EDgGMjdd98dECOAkpIS+vfvT3Z2dida2zUwS/stY/eU7VggV1XdIqIAquoQkSyb+zG0I+Xl5axbt46ysjIiIyN57733WLhwIQC33HJL4IAsWKtuaWlpDB06tLPM7TIYMdo7do+QnEB6cIGI9AOKbO7H0E4UFhby6aefUldXR+/evfnyyy+59957Abjuuus4/vjjA3UrKiqIiYlhzJgxREWF3fpI2FFdXU1VVZURoxaw+6foCeBlEbkRiBCRCcCdwCM292Owmbq6OjZv3swPP/xAZmYmUVFRrF69mttvvx2AK664glNPPTVQv7y8nOjoaMaPH2/OpbWSAw88MOCDM2LUNHYL0p1AHdbB11jgOSwxut/mfgw2snnzZtavX4+I0Lt3byIiIti8eTNXX3019fX1nHLKKZx99tkUFxcDlhM7OTmZ8ePHm+Mge8HhcPDGG28we/ZswBIlQ/PYLUgZqnovcG9woYhkYgVbM4QZZWVlfPfdd2RmZhIZGcnatWuZP38+//3vfwEYN24cV199NcXFxeTm5jJ48GA8Hg9xcXFmmrYXQkOINIiSoXns/okqYM/DsWCdUUtvotzQibjdbr755ptApo+CggIuueQSKisriYmJ4dhjj+VPf/oTFRUV9OnThxEjRpjg+60k1IF94okndrZJXQK7BWmPfDUikoi1J8kQZqxfv576+nqSk5NxOBxcdtllVFZWcuSRR3LrrbeSmJiI1+ulvLyc4cOHGzFqJWY1bd+xRZBEZBOgQJyIFITczgResaMfg30UFRWxbds20tLSeOaZZ1i4cCFVVVWMGjWKuXPnEhsbC1graf379zeO61ZixGj/sGuEdC7W6OhNIHiirEBRE/GwDZ1IbW0t3377Lampqdxwww189NFHAEyaNKmRGKkqHo/HnNZvA7NmzTJitB/YIkiqugxARLJVNTTioyGMaDgGEhERwdKlS/noo49ITExk7ty5HHrooYEssapKSUkJubm55lhDG3jwwQcBeOKJJ4wY7QO2+pBU1SkiBwOHYU3VJOje7Xb2Zdg3Nm3aRHFxMREREYENj1dddRVTp04N1CkpKUFVyczM7BGZLvaXmpqaQFSDAQMGsHjx4k62qOtiq5dSRP4IfAn8GrgJmICVLWSEnf0Y9o3S0lLWrVtHRkYG9957L5WVlUyZMqVREPnKykoSExM5+uijGTdunAkfshccDgc/+9nP+Otf/9rZpnQL7F42uR74tar+Bqj1/3saVsxrQyfi9XpZtWoVqamp/PDDDyxdupSYmBhuvPHGwDStvr4el8vF6NGjiY6O7mSLw58GB/batWt5/vnnqamp6WyTujx2C1KWqi73v/eJSATwb5oIuGboWHbu3InL5SI2NpbHHnsMsOIZNZzQd7lclJaWMnLkSBNCpBU0FULEjCb3H7v3IW0Xkf6qugUrtOzxWDu0623ux9AGPB4PGzZsIDU1lXXr1rF8+XJiYmI455xzAGsqFxUVxcSJE8nIyOhka8Mfs7TfftgtSPcCBwNbgLnAy0A0cKXN/RjawM6dO6mrqyMhIYH777eOFZ566qlkZmZSW1tLr169mDJlipmmtQIjRu2L3atsjwe9XywiaUCMqlbY2Y+h9bhcLtatW0daWhoLFizg66+/Ji0tjRkzZqCqVFRUMHnyZCNGrcTlclFbW2vEqJ1o19ORqurypzO6U1VvaM++DE3TcIp/yZIlPPvss0RFRTFv3jzS09MpLS2lf//+pKebY4atJS8vj+XLlxMbG2vEqB2wzaktIueIyP0icpGIRPnzpd0DbAYOsasfQ+txOBxs376d8vJy7rzzTgCuvfZaxo4dG9iFPXjw4E62MvwpKiri4YcfDlzn5eUZMWon7DrLNg+YDnyKFT97EjAZ+Bo4TFW/saMfQ+uor69n27ZtrFu3jsTERK655hrq6ur41a9+xW9/+1vA2m+Uk5MTOCZiaJqGVEUNIUQuvPDCTraoe2PXlO0M4HBV3Sgiw7BSX5+pqi/a1L5hL3g8Hv773/9SXl6O1+slIiKCjIwMHnnkEb7//nv69OnTKPeXy+UyZ9T2QqgD+3e/+11nm9TtsUuQUlV1I4Cqfi8iNUaMOpYffviBsrIyMjIyAhsd16xZwz//+U9EhL/85S8kJlp5Nt1uN7GxsaSmpnamyWGNWU3rHOwSJPEH8284u+YJuUZVt9rUlyGEsrIyCgoKOOCAAwJi5Ha7uf322/H5fEyfPp0xY8YE6judToYOHRqoa2iMEaPOwy5BSsByXgf/hG8Jeq9ApE19GbCC8peUlFBWVkZhYSFJSUlERETg8XhYvXo1L7/8MgUFBeTl5XH++ecHnqupqSEqKsrkUGuBc88914hRJ2GXIJlNLB1EbW0tBQUFbNu2DYCYmBiSkpKIjo7m/fffZ8GCBRQWFgIQERHBnDlzAo7ruro6ampqmDRpkgnO3wINIUQWLlxoxKiDsSsekteOdgzNo6ps2bKFdevWERUVRUZGRiCkrKpy9dVXBwKt9e3bl8MPP5xf/OIXHHzwwYDlxHY6nYwfP56UlJRO+xzhSnV1deAMX15eHm+++WYnW9QzMWkjuggFBQWsW7cukDMtmLfffjsQaO3SSy/lxBNPbFSnurqa2tpaJk6caDZBNkGDz+jUU0/l1ltv7WxzejRhF7VdRNJF5DURqRaRLSJyVgt1DxGRj0WkSkSKROTyjrS1o9ixYwfr1q2jd+/ee4hRVVUVDzzwAGAFWvvd737XqI7b7cblcjFlyhQjRk0Q7MB++eWXqa42kXI6k3AcIT0EuIEsYAzwbxFZFRqX25/r7V3gT8C/gF5Abgfb2q6oKoWFhXz77bdkZmY2mfXjscceo6SkhJEjRzZKc91AWVkZhxxyCElJSR1hcpeiqdU0E3qlc7FdkEQkCitSZF9V/ZeIxAGoam0rnk0AfgccrKpVwCci8ibWLvDrQ6pfCSxR1Wf913XA9zZ9jE7H5/OxZs0atm/fTkZGxh4jo7q6Ou677z5eeeUVRISrr746IFgejwev1xvYjZ2VldUZHyGsMUv74YndIWxHAOuAp4FF/uJjgCda2cQQwKOqG4LKVtF0CNxJQKmIfCoiDhF5S0S6zdbjsrIytm/f3uQ07auvvmLGjBm88sorREdHc9NNNzFihPUVVVdXU15eDkBmZibDhg0z+41CMGIUvtg9QnoYmKuqi0SkzF+2HPhHK59PBEKzllQATc03crEO7f4cWA3MA54HDg2tKCLnAecBXea4xPbt24mLi2skJnV1dcyZM4f3338fgNzcXO666y6GDh0KWCtpdXV1TJkyxUzRWsDtdlNXV2fEKAyxW5BGAv/0v1cAVa0SkdbG9qxiz1TcyUBlE3VrgddU9SsAEbkN2C0iKaHxl1T1UeBRgPHjx2srbek03G43u3btauSErq+v57rrruOTTz4hPj6emTNnctZZZwX2GHm9XioqKowYtYLc3FyWL19Or169jBiFGXavsm0BxgYXiMh44MdWPr8BiBKR/KCy0ViHdUP5Fr/o+Ql7oWktDWmIgvcZ3XzzzXzyySekpKTwxBNPMGvWrEYn9cvLy8nPzzfn05rB4XDw4IMPomr9mOTm5hoxCkPsHiHdgrUq9negl4hcA1wMtCpmg6pWi8irwO0ici7WKtuJwJQmqj8JvCIiD2AJ1s3AJ90hOuXmzZsDB2EBPv/8c95//30SEhJ46KGH9ohh5PF4EBH69+/f0aZ2CYJ9RgCXXHJJJ1tkaA5bR0iq+iYwDegH/Ac4CDhdVd9pQzMXAXGAA8sndKGqrhWRw0SkKqivD4AbsbKaOIDBQLN7lroKVVVVlJeXN8oW+/jjVmTgWbNmBfxFwZSVlZGfn0+vXr06zM6uQqgD+7TTTutskwwtYOsISUTS/D6dr/a1DVUtpYm0Saq6AsvpHVz2MJYjvduwZcuWRsLy3//+l2+++Ybk5GROOeWUPerX1tYSExNDbm632oJlC2Y1rethtw+pUETeFJHTG/YfGVpPTU0NW7duDZw1U1UWLlwIwJlnnrnHpj2Px4PT6WTs2LF7bA3o6Rgx6prYLUgHAu9j7Z7eJSJPi8ivRMSEHmkFW7duJSoqChHB6XRy9dVX8+WXX5KQkMDpp5/eqK6qUlJSwvDhw40juwnOO+88I0ZdELvTIBUBDwAPiMhALJ/O/wGZWEdBDM1QUlLCpk2byMzMpKamhpkzZ7J161YSExOZO3cuyck/7YZwu92UlpZy4IEHGkd2MzSEEHn00UeNGHUh2nOcn+J/JQHmxGIzuFwu1qxZg8PhIDk5mYiICBYvXszWrVsZMGAAf/vb3xr5h+rq6gJhRMyRkMZUVVWRkJCAiJCbm8vrr7/e2SYZ2ojdR0eGiMgcEVkPvAPEAmeo6kA7++ku1NXV8fXXX1NRUUFWVhZxcXGoKi+99BIAF1xwQSMxUlXKysoYO3asEaMQioqKmDhxIjfffHNgr5Gh62G3D+krLD/SZViHay9V1U9t7qNb0JAlxOVyNfIBffHFF2zevJnevXtz5JFHNnqmtLSUvLw8I0YhBKcqeu2110wIkS6M3VO2LFV12dxmt6SiooLy8vI9/BsvvmglawmNa+TxeIiIiOCggw7qUDvDnaZW04I3lRq6FvstSCJypqo+7788rbmT5ar61P721Z2orq7eI77Rli1b+OSTT4iOjubkk09udM/pdJKXl2c2PwZhlva7H3aMkGZi7agGmN1MHQWMIAVRUVGxR6D9f/zjH6gqJ5xwwh7RHT0ej8kUEoQRo+7JfguSqv4i6P1h+9teT6G8vLyRIK1bt46lS5fSq1cvzj333EZ1PR4PMTExjZb+ezoejwePx2PEqJth99GRr1R1QhPln6vqJDv76sp4vV6qq6s54IADAmUPPfQQAKeffvoeTmun08mAAQNMoLUgcnJy+PDDD4mKijJi1I2we5Vtz5OfFkNs7qdLU1v7UzRfr9fLfffdx2effUZCQgLnnHPOHvU9Ho/5pcOaps2fPz+wrJ+Tk2O+l26GLSMkEWkIUdsr6H0DA+hGsa7toEGQ3G43119/PR9//DFRUVHceOONexwDqaioICEhocdP10JDiFx+ebdMMNPjsWvKVtjMewW+Bl60qZ9uQWVlJZGRkfz973/n448/Jjk5mXvuuYdx48YF6qgqxcXFpKenM2rUqB49XQt1YJ955pmdbZKhnbArc+3NEPAV/duONrsz5eXlrFmzhmeeeYbIyEjmz5/PyJEjG9WpqKggJyeHgw8+mMjInns22aym9Szs2Id0qKr+x39ZKSKHN1VPVT/e3766C9u3b+fOO+8E4Nxzz91DjFQVt9vNoEGDjBgZMepR2DFCepyfnNnPNlNHga6R7qOdqa+v54MPPqC4uJjhw4fzhz/8YY86FRUV5Obm9vgdxxd8xSjSAAAgAElEQVReeKERox6GHfuQhga977e/7XV3Kisr+fzzzwE46aST9gis1jA6OvDAAzvDvLCiIYTIww8/bMSoh9CuYQZF5DDAaw7YWng8HlauXMk333wDwGGH7bmPtKysjLy8vB47OqqsrCQxMRERoU+fPrzyyiudbZKhA7E7/MhyEZnqf3818CpWZpDr7Oynq7JhwwY+//zzQJLC4I2RYAkWsEdWkZ6Cw+Fg0qRJ3HDDDSaESA/F7o2RI4HP/e/PB44EJmJlEunRFBcXs3nzZv73v/8BcPjhe/r+S0tLGTZs2B5n3HoCwQ7st956i6qqqr0/ZOh22D1liwB8/vC1Uaq6FkBE0lt+rHvjcrlYtWoVycnJrFixAthTkJxOJ8nJyeTk5HSGiZ1KU6tpJvtuz8RuQfoU+BuQA7wG4BenEpv76TKoKmvXrkVEWLt2LSUlJWRlZZGf/1Ny3pqaGnw+H2PHjt0jJEl3xyztG4Kx+6d/JuAC1gNz/GXDgQU299NlKCoqYteuXcTFxXH33XcDcMIJJwR2XtfX11NdXc2ECROIj4/vTFM7HCNGhlDszjpSDFwbUrYYWGxnP10Fn8/H+vXrSU1N5bHHHmPTpk3079+/0d6jsrIyRo0a1SPPqvl8Pnw+nxEjQwC7w49EATcA04G+WOfangbuUtV6O/vqChQXF1NTU0NFRQVPPfUUIsKcOXOIjY0FrEO2SUlJ9OnTp5Mt7Ryys7P58MMPiYiIMGJkAOyfst0NHA9cAfzM/++vgLts7ifsUVU2btxIUlISzz77LD6fj1NOOYVRo0YF6jidToYNG9aj/EYOh4N77703sKyfnZ1txMgQwG6n9mnAWFXd7b9eKyJfAd8AV9ncV1hTWlpKZWUl0dHRvPfee4gIZ599duD+7t27ycrKIiMjoxOt7FhCQ4hcdVWP+pEwtAK7/zRHAr6QMl879BP2FBYWEhsby2uvvYbb7Wbq1Knk5ubi8/koKioiOzub0aNHd7aZHUaoA3v69OmdbZIhDLFbKP4FvCkix4hIvogci7X8/y+b+wlr6uvr2blzJ7GxsYGjD6eddhpgpcweOHAgo0aN2uMcW3fFrKYZWovdgnQN8DFWBIC1wGPAf/zlrUJE0kXkNRGpFpEtInLWXur3EpHvRWT7/hhuJ2VlZfh8Pj7//HOKiorIy8tj4sSJ1NbWEhcXx+DBg3tMwDUjRoa2YPeyfx1wo/+1rzwEuIEsYAzwbxFZ1bDruwmuAYqBsNnau337duLj4/nqq68A+MUvfoGI4HQ6mTRpUo8ZGQFcfPHFRowMrcaWEZJ/evaxiJSKyPsisk+xj0QkAfgdcLOqVqnqJ8CbWNsImqp/IPB74M59td1u3G43DoeDhIQEVq1aBcDYsWNxOp3k5ubukW+tu/Pggw9yyimnGDEytAq7pmwPYu05mgnsxjo+si8MATyquiGobBUwopn6C7BGY7XN3O9wSkpKUFVqa2tZv349kZGRjBw5ErfbTd++fTvbvA7B6XQGlvWzsrJ4+eWXjRgZWoVdgjQOmKWqb2Jlr524j+0kAs6QsgqamI6JyMlApKq+trdGReQ8EVkpIiuLi4v30bTWUVBQQFJSEqtXr8br9XLQQQcRGxuLiJCSktKufYcDDoeDyZMnc9VVV5kQIoY2Y5cg9VLVWgBVrQTi9rGdKiD0DEUyUBlc4J/azQMua02jqvqoqo5X1fGhMYjspLKyEqfTSVxcXCAI25gxYwJJIbu776ioqCjgwF6yZAmVlZV7f8hgCMKu35AYEbkl6Dou5BpVvb0V7WwAokQkX1U3+stGY63YBZOPle9thX+1qheQIiK7gEmqurntH2H/2bFjR0B0GgRp7Nix1NbWMmRI986VWVRUxNFHH93Igd0Tz+cZ9g+7BOklLJFo4F8h160au6tqtYi8CtwuIudirbKdCEwJqboGCI7fPQXLj3UI1opbh+PxeNiyZQupqal4PB5Wr14NwOjRo/F4PKSlpXWGWR1CU2JkfEaGfcGuvGx2bru9CHgCcGDFUbpQVdf643O/o6qJquoBdjU8ICKlgE9VdzXZYgdQVlaG1+slMjKStWvX4nK5yMvLC4QUaThQ290IFaMPPvjAiJFhnwk7p4aqlgInNVG+Asvp3dQzy4Hc9rWsZXbv3k2vXr0AWLZsGQCHHHIINTU1DBw4sDNNa1caNniakZHBDsJOkLoqu3fvJj4+nqqqqsBxkZNPPhmv19utw7H27t2bDz/8MPDeYNgfetyh1/agIepjr169ePXVV6murmb8+PGMGGFtn+pu0zWHw8G8efMCy/q9e/c2YmSwBTNCsoHq6mrA2qX93HPPATBjxozA/e4kSKEhRK699tq9PGEwtB7bR0gicpSIPCIir/uvDxGRI+zuJ5yorKxERHjvvffYvXs3gwcPZvLkyXi9XqKiogK+pa5O6EHZmTNndrZJhm6G3YkiL8I66b8NOMpf7Ab+amc/4cbu3buJjY3l1VdfBeDMM89ERKivr+82e3GCNz0aB7ahvbB7hHQVcKyqzuWnQG3fA8Ns7idsUFVKSkrYuXMn3377LQkJCRx33HEA1NXVdQtBMvuMDB2F3YKUBGzxv2/YDBmFNUrqltTW1lJfX8+bb74JWKFG4uKskzPdZYR02WWXGTEydAh2O7U/Aa7GCvbfwMXARzb3EzZUV1dTX1/P22+/DVhL/Q2oakCcujIPPvggAAsWLDBiZGhX7BakS4HFIjIbSBKRtVijo1/b3E/YUFxczFdffUVFRQUHHXQQw4Y1np121RW2iooKkpOTEREOOOAAXnzxxc42ydADsHXKpqqFWOfJzgFmAOcD41V1p539hAuqys6dO1m71jr72+A7argnIsTExHSWefuMw+FgypQpXH755SaEiKFDsX0fklo/wf/xv7o1VVVVuN1uvv/+ewBGjhwZuFdfX098fHyXy7kWus/I6XT2iDhOhvDA7sy1m2jmZL+qdrsDXaWlpXg8HjZu3EhERARDhw4N3HO73aSmpnaidW2nqYD8RowMHYndI6RzQ677YPmVnre5n7Bgx44d7Ny5E6/Xy6BBgwIn+8Fydufn57fwdHhhsoMYwgG7s44sCy0TkWXA2+x7nO2wpK6ujoqKCjZutOLIDR8+PHCvvLycrKws+vTp01nmtQkjRoZwoSMcHLVAt5uulZeXAwR8LQ0Had1uNz6fjxEjRnSZ3GsRERFERUUZMTJ0Onb7kG4JKYoHjgfes7OfcGDTpk3Ex8fvIUjl5eWMHj26Sy33Z2ZmsmzZMnw+nxEjQ6di9wgpP+SVipX4sVslcq+oqKCsrAxVZfPmzURHRzN48GC8Xi8RERG0ZyIBu3A4HNx1112BZf3MzEwjRoZOx7YRkohEAkuBl1TVZVe74cjmzZuJiYkJLPcPGTKE6OhoysvL6devH9HR0Z1sYcuELu1ff/31nWyRwWBh2whJVb3Agu4uRjU1NRQWFpKcnBzILNLg0O4KySBDHdizZs3qbJMMhgB2T9n+LSLd9pgIWCffIyMjUVUWL14MwNSpU3G5XCQkJIT1YVqzmmYId+zehxQBvCoin2DFRApsklTVbvGnuLi4mISEBFauXElhYSHZ2dlMmjSJsrIyhg0bFrYra0aMDF0BuwVpI3CPzW2GDT6fj/LyctLT03ntNSuD97Rp0wIjpnAeHV1xxRVGjAxhjy2CJCJnqurzqnqzHe2FKzU1Nfh8PioqKli+fDkiwrRp0wL3g3dqhxsLFiwA4G9/+5sRI0PYYpcP6RGb2glrqqqqAHj77bepr69n8uTJZGdn4/F4iImJCbvY2eXl5fh8VuDOjIwMnnvuOSNGhrDGLkEKT8eJzZSWltKrVy+WLFkCEBgd1dbWkpGR0Zmm7YHD4eDQQw/l4osvDoiSwRDu2OVDihSRo2hBmFT1A5v66jSKi4spLS1l7dq1xMXFMXXqVMA61xZOgtRUCJGuFnnA0DOxS5BisLKNNCdIShc/z+Z2u6mpqWH58uUAHHHEEYHjIapKQkJCJ1r3E02tphkxMnQV7BKk6u4Y7yiY6upqRISlS5cCjaNDQng4tM3SvqGr07XCGXYilZWVbN++nQ0bNpCYmMikSZMAwsahbcTI0B0wTu1WUlZWxmeffQbAUUcdFRCgcHFoR0VFERMTY8TI0KWxZcqmqkl2tBPOOJ3OwNm1ww8/PFAeLoKUnp7O+++/j8fjMWJk6LKE3ZRNRNJF5DURqRaRLSJyVjP1rhGRNSJSKSKbROSa9rKpYTNkQ3aRsWPHAuByuYiNjSUrK6u9um4Rh8PB3LlzA8v66enpRowMXRrbs47YwENYudyygDFYB3ZXqerakHqClWrpW2AQ8J6IbFPVF+w2qK6ujvXr11NXV8fAgQMDq1bl5eVMmDCBqKiO/xpDl/b//Oc/d7gNBoPdhNUISUQSgN8BN6tqlap+ArxJEwHeVHWeqv5XVT2quh54Azi0PexyuVysWbMG+Gl05HQ6ycrK6pRgbKEO7PPOO6/DbTAY2oOwEiRgCOBR1Q1BZauAES09JNYR+8OA0FFUw/3zRGSliKwsLi5us1Eul6vJ6dqgQYM6/HS/WU0zdGfCTZASAWdIWQWwN6f5rVif5cmmbqrqo6o6XlXH78uIpmF3NliCpKpERkaSlNSxvnwjRobuTrj5kKqA0BgeyUBlcw+IyCVYvqTDVLWuPYz6+uuvqampoW/fvmRlZVFdXU1GRkaH+46uvPJKI0aGbk24CdIGIEpE8lV1o79sNM1PxWYB1wOHq+r29jLqiy++AH6artXU1DBo0KD26q5ZHnjgAQDuu+8+I0aGbklYTdlUtRp4FbhdRBJE5FDgRODp0LoicjZwB/BzVS1oL5vq6+tZt24dAKNGjWqws8NSTJeVlTVa1n/mmWeMGBm6LWElSH4uAuIAB1YK7gtVda2IHCYiVUH15gIZwFciUuV//cNuY1wuFzt27ACgf//+gaMiHXGY1uFwMHXqVM4//3wTQsTQIwi3KRuqWgqc1ET5Ciynd8P1gR1hj8vlYufOnQD069ePqqoqcnJy2n11LXSfUUVFBWlpae3ap8HQ2YTjCCms2LlzJxUVFcTExJCZmYnb7W73vUdNraYZMTL0BIwg7YWNGy3fet++fYmIsL6u9gzmb5b2DT0ZI0h74ccffwSs6ZrL5SI5ObndQo0UFRUZMTL0aIwg7YUtW7YAliBVV1fTp0+fduurV69exMfHGzEy9FjCzqkdbmzduhWA3NxcvF4v6enp7dZXWloa7733HvX19UaMDD0SM0LaCw2C1OBDsvu4iMPh4Lbbbgss66elpRkxMvRYzAhpLxQWFgJwwAEHkJmZSWRkpG1thy7tz5kzx7a2DYauiBkhtUBVVRUlJSVERUWRlJREdna2bW2HrqZdeOGFtrVtMHRVjCC1QMOSf05ODhEREbZlFjFL+wZD0xhBaoEGQcrNzQUgJiZmv9s0YmQwNI8RpBYI3oME9gjSNddcY8TIYGgG49RugR9++AGwpmy9evWyxaE9f/58AO655x4jRgZDCEaQWqDhlH9mZuZ++Y/KyspISUkhIiKC1NRU/vnPf9plosHQrTBTthYoKioCrLNr+7r/qCGEyLnnnmtCiBgMe8GMkFrA4XAAliDtS/wjE0LEYGgbZoTUDKpKQ4aS5OTkNk/ZTAgRg6HtGEFqhvLyctxuN/Hx8cTGxrZphc0s7RsM+4YRpGbYtWsXQOAwbWtDjhgxMhj2HSNIzRAqSK0dIcXExJCUlGTEyGDYB4xTuxkaBCktLY2YmJhAtMi9kZKSwpIlS6irqzNiZDC0ETNCaoaGJf+0tDQSExNbrOtwOLj55pvxer2AJUpGjAyGtmMEqRkaRkipqaktClKDz2ju3LnceuutHWSdwU5eeOEFhg0bRkJCAoMGDWLFihWdbVKA+++/n+zsbJKTk5k1axZ1de2SnDlsMILUDK0RpFAH9qWXXtqRJhpsYOnSpVx33XU8+eSTVFZW8vHHHzNw4MB279fj8ey1zpIlS7jrrrtYtmwZW7ZsoaCgoNvHzDKC1AzBPqTY2Ng97pvVtO7BnDlzuOWWW5g0aRIRERH07duXvn37ArB8+XJyc3OZN28evXv3pk+fPrz++uu8/fbbDBkyhPT0dO64445AWz6fj7vuuotBgwaRkZHBaaedRmlpKQCbN29GRHj88cfJy8vj6KOP3qtt//znP/njH//IiBEjSEtL4+abb2bRokXN1n/qqafo378/GRkZ/OUvf2HAgAG8//77AHz55ZdMnjyZ1NRU+vTpwyWXXILb7Q48KyI88MADDBw4kMzMTK655prAyYIffviBI444gpSUFDIzMzn99NPb/D23GlXtUa9x48Zpaxg1apQCumDBAq2srGx0r6ioSIcPH66ADh8+XIuKilrVpiG88Hg8Gh0drXfeeacOGjRI+/btqxdffLHW1NSoquqHH36okZGRetttt6nb7dZHH31UMzMz9cwzz1Sn06lr1qzR2NhYLSgoUFXVv/3tbzpx4kTdtm2bulwuPe+88/SMM85QVdVNmzYpoNOnT9eqqqpAHy0xatQofeGFFwLXxcXFCuju3bv3qLt27VpNSEjQFStWaF1dnV511VUaFRWlS5cuVVXVlStX6meffab19fW6adMmHTp0qN5///2B5wE98sgjtaSkRLds2aL5+fn62GOPqarqGWecoXPnzlWv16u1tbW6YsWKNn/XwEptxe9npwtER79aK0hZWVkK6LPPPqter7fRvZkzZxox2keADnm1hsLCQgV03LhxumPHDi0uLtYpU6bojTfeqKqWIMXGxqrH41FVVafTqYB+/vnngTYOOeQQfe2111RVdejQofr+++8H7u3YsUOjoqICIgDojz/+2OrvauDAgfrOO+8Ert1utwK6adOmPeredtttAfFTVa2urtbo6OiAIIVy//3360knnRS4Bhr19dBDD+nRRx+tqqrTp0/X2bNn67Zt21pteyitFSQzZWsCr9cbODaSm5u7x5L//PnzmTlzppmmdXHi4uIAuPTSS+nTpw+ZmZlceeWVvP3224E6GRkZgbAzDfWzsrIatVFVVQVYKbNOPvlkUlNTSU1NZdiwYURGRgZWbOGn2FqtITExEafTGbhueN/UQe8dO3Y0ajs+Pp6MjIzA9YYNGzjhhBMCDvIbb7yR3bt3N2oj+Pn+/fsHol3MmzcPVeVnP/sZI0aM4Iknnmj1Z2grRpCaYPfu3fh8PpKTkwOCU1paGljWT05O5sknnzRitA+05q+kHa/WkJaWRm5uLiISKAt+31b69evHO++8Q3l5eeDlcrkCPqm2tj9ixAhWrVoVuF61ahVZWVmNhKaBPn36sH379sB1bW0tJSUlgesLL7yQoUOHsnHjRpxOJ3fcccce39O2bdsC77du3UpOTg4A2dnZPPbYY+zYsYNHHnmEiy66KBArzG6MIDVB8ApbWloaDoeDww47jHPOOScgSobuwR/+8AcWLFiAw+GgrKyM+++/nxNOOGGf2rrgggu46aabAslFi4uLeeONN1p8ZsCAAc06qmfMmMHjjz/Od999R3l5OXPnzmXmzJlN1j3llFN46623+PTTT3G73dx6662NBKeyspLk5GQSExNZt24dDz/88B5t3HPPPZSVlbFt2zbmz58fcF6//PLLAbFLS0tDRFq9UbitGEFqguAVturq6sBq2v/+9z8qKio62TqDndx8881MmDCBIUOGMGzYMMaOHctNN920T21dfvnlTJs2jeOOO46kpCQmTZrEF1980Wx9t9tNSUkJkyZNavL+L3/5S6699lqOOuoo8vLy6N+/P7fddluTdUeMGMGCBQs444wz6NOnD4mJifTu3Ttw5On//u//eO6550hKSmL27NlNrpSdeOKJjBs3jjFjxnD88cfzxz/+EYCvvvqKiRMnkpiYyLRp05g/f377bY3oqCF0G4ba6cBrQDWwBTirmXoC3A2U+F93A7K39lvj1F60aJECOnXqVB02bJhxYBvahRUrVjRyRNtJZWWlRkZGBlYA9wagGzdubBdb/O23yqkdjmfZHgLcQBYwBvi3iKxS1bUh9c4DTgJGY62sLAU2Af/YXwMaRkjffvstTqfT7DMytAtTp05l6tSptrX31ltvccwxx6CqXH311YwcOZIBAwbY1n5HEFZTNhFJAH4H3KyqVar6CfAmML2J6ucA96rqdlUtBO4FZtphx6ZNmwCMGBm6FG+88QY5OTnk5OSwceNGXnjhhf1y0ncG4TZCGgJ4VHVDUNkq4Igm6o7w3wuuN8IOIxpWJ7Kzs40YGboMCxcuZOHChfv0rLZyZbK9CTdBSgScIWUVQFMR9hP994LrJYqIaMi3KyLnYU3xyMvL26sRU6dOpaqqiosvvtiIkcHQgUi4KCOAiIwF/qOq8UFlVwFHqupvQupWAD9X1S/91+OA5araYnqQ8ePH68qVK+033mAwNIuIfK2q4/dWL6x8SMAGIEpE8oPKRgOhDm38ZaNbUc9gMHQRwkqQVLUaeBW4XUQSRORQ4ETg6SaqPwVcKSJ9RSQHuApY1GHGGgwG2wkrQfJzERAHOIDngQtVda2IHCYiVUH1HgHeAlYDa4B/+8sMBkMXJdyc2qhqKdb+otDyFViO7IZrBa71vwwGQzcgHEdIBoOhh2IEyWAwhA1GkAwGQ9hgBMlgMIQNYbUxsiMQkWKsKAJ7IxPYvddanUu422js2z/C3T5ovY39VfWAvVXqcYLUWkRkZWt2lnYm4W6jsW//CHf7wH4bzZTNYDCEDUaQDAZD2GAEqXke7WwDWkG422js2z/C3T6w2UbjQzIYDGGDGSEZDIawwQiSwWAIG3q0IIlIuoi8JiLVIrJFRM5qpp6IyN0iUuJ/3S0dEKy4DfZdIyJrRKRSRDaJyDXtbVtb7Auq30tEvheR7S3V6wz7ROQQEflYRKpEpEhELg8nG0UkRkT+4betVETeEpG+TdW10bZLRGSliNSJyKK91P2TiOwSEaeIPCEiMfvSZ48WJBpnODkbeFhEmorLHZzhZBTwG+D8MLJPgBlAGvBL4BIROSOM7GvgGqC4A+xqoFX2iUgm8C5W+JoMYDDwXjjZCFwOTMb6+csByoAF7WzbDmAu0GLubBH5BXA9cAzQHxgINJ1Abm+0JldSd3wBCVg/CEOCyp4G7mqi7qfAeUHXfwQ+Dxf7mnj2AWBBONkHHAh8D/wK2B5m/793AE935M/fPtj4MDAv6Pp4YH0H2TkXWNTC/eeAO4KujwF27UtfPXmE1FyGk6b+OrVbhpMWaIt9AfxTycNo/3C+bbVvAXAjUNvOdjXQFvsmAaUi8qmIOPzTob1ng+hYGx8HDhWRHBGJxxpNvdMBNraGpn4/skQko60N9WRBsiXDSTvZ1tBna+0L5las/9cn28GmYFptn4icDESq6mvtbFMwbfn+crHy/F0O5GElHH2+Xa2zaIuNG4FtQKH/mWHA7e1qXetp6vcD9v6zugc9WZCqgOSQsmSgshV1k4Eq9Y9P24m22AdYTkgsX9LxqlrXjrZBK+3zJ/+cB1zWzvaE0pbvrxZ4TVW/UlUXlv9jioikhJGNDwExWD6uBKzY8+EyQmrq9wNa+Fltjp4sSOGe4aQt9iEis/A7FlW1I1axWmtfPjAAWCEiu7B+kfr4V2QGhIF9AN9ipWNvoKN2C7fFxjFYfpxS/x+bBcDP/A75zqap348iVS1pc0sd7cgLpxfwAtbQPAE4FGuoOaKJehdgOWT7Yq1wrAUuCCP7zgZ2AcPC7fvDitueHfT6LdbqTTbWNC4cvr+jsVatxgDRwP3AinD5Dv31ngReAVL8Nt4IFLazbVFALHAnlrM9Fohqot4v/T9/w4FU4ANasfjSZJ8d+QMcbi8gHXgdqAa2Amf5yw/DmpI11BOsaUep/zUP/7GbMLFvE1CPNXRueP0jXOwLeeZIOmCVra32ARdi+WfKsLLZ9AsnG7Gmas9iZeMpBz4BftbOtt2KNVoMft2K5WerAvKC6l4JFGH5t54EYvalT3OWzWAwhA092YdkMBjCDCNIBoMhbDCCZDAYwgYjSAaDIWwwgmQwGMIGI0gGgyFsMIJkQESeEZFbO9uOvSEi60XksBbuvyciZ3ekTQZ7MYLUjRCRzSJS6w8y1vDK6SRbnhERt9+GUr9YDNmfNlX1IFVd4W9/bmjQMFU9TlWf3Z8+QhGRKBFRfwC1KhHZLiL3iEirfndE5FgR2WynTd0ZI0jdj9+oamLQa0cn2nKHqiYC/bB2uLcY6CvMGeH/LEcD07GiAxhsxghSD0BEIkTkX/4DreUislxEhjVTt7eIvO2vVyoiHwfdy/WHWy32h8q9uDX9q2o11nmtg/3txIrIAyKyU0QKReQ+EenViv63i8iRInICcC1wtn/U8rX//iciMlNE4vyhVIcGPZvtHz1m+K+nicgqfz+fiMjBrfwsG7AC9o0JavtcsULzVorIjyJyrr88BesYSl7QiLW3///jRn/d3SLygoiktab/7o4RpJ7DYqyT99nAGqzDkk1xDVAAHOCv+2ewRM3fxldYh4x/DlwjIsfsrWMRSQLOAv7nL7oFGI8VjnUs1qHSG1rqPxhVXYx1nvBZ/yhwXMj9WqzzYWcGFZ8OLFPVEhGZADwGnIt1RuwJ4I0GUdzLZxnmt/eHoOIirAiOycBsYIGIjFLVCqxwx1uDRqwO4E/++odjxWKqwory2eMxgtT9eN3/V79cRF4HUFWfqi5S1Uq14v3cCozzxyoKpR4rokGeqrpVtWGEMhlIVtU7/OU/YEUxbCl29/UiUo4VZiMGmOUvPxu4VVWL/b+gt2NNg1rqv608R2NBOstfBlaM9L+rFf/Iq6oNU8kJLbT3rYhUA98BS7HibwOgqm+paoFafAAswzoc2xwXADeqaqH+FH/p1Nb6pbozPf4L6IacpKqp/tdJAKDuicQAAAcpSURBVCISKSLzRKRARJz89Ne9qVg6dwFbgGX+KUVDBpP+WFOPBrErx5o2Zbdgy11+O/qo6kmquslfnuPvo4EtWKOulvpvK+8DqSIyTkQGYYXGeCPos1wX8ln6BNnQFKOwIiCehSXOATEXkRNE5Av/FLMcOI6mv9sG8oC3gvpe7S/v3faP2b0wgtQzmAH8Gsshm4KVVQOssCqNUFWnqv5JVQdgZVq5TkSOwAqfujFI7FJVNUlVf7MP9uzAEoUG8rBCf7TU/x6mttSBqnqAl7FGSWcBb/p9Wfg/y20hnyVeVV/aS5s+VX0eWAncBCAiccC/sGIGZalqKlbGkobvtik7twM/D+k/VlV3tdR/T8AIUs8gCagDSoB44K/NVRSR34jIIBERrGBhXsAHfAa4ReQqv1M6UkRGisi45tpqgeeBW0QkU0QOAG4GntlL/6EUAQP89ZrjOSzfUfB0DSz/0cUiMkEsEv39NjWFbYq7gAv8tscAvbDSO3n9Dvdgv1oRkOn3ozXwD+AO8ScS8Du6p7Wy726NEaSewZNYo5IdWNEuP22h7kFYEf+qgP8A81V1hX/E8WvgZ8BmYDeWHyU0JnRruA0rM8UarPCxX2CNMJrtv4k2XsQSglIR+bKZfj4FPFgO8kCeNVX9HCsg28NYAdk2AL9vrfGq+j8sgb5aVcuxnNSvYW1tOAXL+d9Qdw1WpMfN/ilab+A+rDxwy0Sk0m9nS/6rHoMJ0GYwGMIGM0IyGAxhgxEkg8EQNhhBMhgMYYMRJIPBEDYYQTIYDGGDESSDwRA2tChIIhIjIo+LyBb/SeZvRORXIXWOEZF1IlIjIh+KSP+ge/8nIhv9z64TkRkhzz4qVtAtn4jM3JuxIjJGRP6/vbMJrauI4vjvH9I2aYM1RqkfoEWCWFsDdSGirbZSRMmiSEHEggpClBLoInVhECq40i5ExRoEKQVdWGnBD5Ao2Cxs6SYFY4QWlTSSgtaFDU386mJcnHkwublfecrjLc4PBvLumZlz3pu8eTNz5z93MvqalJQqrndG//Oqcf6MpBclTcfYZlKJgqRUnd1IQdJIQV1bJI1H5fayfRSSNkn6Osb2o6THS+K6T9JXUYbwm6SPJd2U2NdIGpP0a8zzmaRcyUNV+1X5cpxWUzVC6sS22T+ESQ5eBo4pPpNd9lzxE9hO2+uwLfUfJeUXMbXzeuz8mDcl3Z/YvwX2AWerApUpsT/BdvT2AkdZqtBexFTbdbVPwiQVvdijgIclPQkQQkjV2T3A3dhu4eMFdV0FjgHP5cTdGeP+HPuMhoAPVHxYWS/wHrARk1dcwTY2NtiPaakGME3Y79hz3vMobb8avhyntTTxeN0pYE/8ewg4ndjWAX8CdxaU/RQYybn+DfBshd9HML2Tkms/A49m8u0CLjTxvt4C3i6wHQRO1qij3z7SJde2YLuO07i/BF6tGdc9wJXk9bvA68nrQeB8M+1X5cuTp1anFa0hSdoA3IHJDwA2Y6OcRue2CPwUr2fLdmPb47/P2mqyGZgKIaRToqk8Xysl6qG258UWbU9jI7L/CxEPK4s+LkvaVpD3wUxc7wMPSLpZ0lrsKI8vkroOSzqc63R5+1X5cpyW0lk3o6RVwIfA0RDCuXi5BxMVpsxjYs4sY1jnNd5EnA1f8zV9rZRXsOlr3nRlG7ABU3Q3w3ngEnaY2RvATmwKdbKRIZhCfBmSBrDDzHYnl3/ApmEXMeHpd8BwUte+grry2q/Kl+O0lLoHlXdgJwz+Q/LPj01FsuLKa7C1iLT8IWxE8ERmhFPmM11QvrWur4K6RpO6xjK2YWwENBhC+Dun+DPA8RDCQp24s4QQrmLHaAwCvwAj2HrTXEXM/djIZ39YKi59B1OY92FT5BMkI6SCuorar8qX47SWqjkdNr04gv2id2dsQ8Cp5PU64A+SNSRM2T0N9JX4qLuGNMfStZhZ/sMaEnaC4Rxwe4G9GxuFPVyzvmVrSAX5TgPPl9hvwxT1L+TYpoHdyetrsTN3rl9p+1X58uSp1ak6g021zgA9ObYb4hd2D9AFvAacSewvYVOMGwvqXh3LncLOIu4COkryzmJ3mdZgv/SzwOpo74jlH4vXuxq2gvr2YiOWTSV5nopfVhXlifkU/d0VO4cuYE1iH4jX1gIHgJnUnqnrFmwd7kCB/Qh2t289sAoYBS422X6lvjx5anUqN9qvZwD+wqZMjbQ3ybMLOIfdXZsANia2gB0MlpYdTewTMU+adpTEsxWYjL7OAlsT246cuiZK6prBbtensY1l8oxT424Ydts86/tCYj+E3Z5fwKZG/ZnyC8D2+PfBWD6NayHJ24etBV0CLmOjy3sT+1jjfVS1X5UvT55anfw8JMdx2gaXjjiO0zZ4h+Q4TtvgHZLjOG2Dd0iO47QN3iE5jtM2eIfkOE7b4B2S4zhtg3dIjuO0Dd4hOY7TNvwLgkkiQ9+AwxEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 288x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAASQAAAEkCAYAAACG+UzsAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xl8VNX9//HXJwnZV7JBAiEiRFmURUQQEONCF0HUfqt+tVJcsK3FfcG6VRGtoi2/ytcFccGKiksBl1apxYpUEAlUVBZZDUsEAtn3TOb8/riT6STMZJ1Jbsjn+XjMg8y95957Jsube8899xwxxqCUUnYQ1NkVUEqpehpISinb0EBSStmGBpJSyjY0kJRStqGBpJSyDQ0kpZRtaCAppWxDA0kpZRshnV2BjpaUlGQyMzM7uxpKdSsbNmw4YoxJbq5ctwukzMxMcnJyOrsaSnUrIpLbknJ6yaaUsg0NJKWUbWggKaVsQwNJKWUbGkhKKduwXSCJyEwRyRGRahFZ1EzZW0XkoIiUiMhLIhLWQdVUSgWA7QIJyAPmAC81VUhEfgTcDZwL9AP6Aw8FvHZKqYCxXT8kY8xSABEZBfRpougvgReNMZtd5R8GXsMKqXZ55ZVX+Oqrr5osExERwa233kpycrN9vZRSLWS7QGqFIcC7Hu83AakikmiMOepZUESuB64HyMjIaHbHf//733nrrbeaLbdx40aee+65Y5YnJCQQFxfX7PZKqYa6ciBFA8Ue7+u/jgEaBJIx5nngeYBRo0Y1O6vBtGnTGDp0KKWlpV7X19XV8cwzz7BixQreeusthg8ffkyZsWPHEhMT02CZMYaysjLKysqora0lMjKSpKSk5qqjVLfRlQOpDIj1eF//tfcUaYULLriACy64oMkyUVFRPPzwwyxcuJDFixcTEvLfb2V5eTk5OTkMHz6curo6qqqqKC0tJT8/n/LyckQEEcHpdDJ27Fji4+PbW2Wljgt2bNRuqc3AMI/3w4BDjS/XAmXWrFn06tWLnTt38v777zdYFxUVhTGGdevWkZOTw+bNm/nhhx8ICQkhJSWF5ORkkpKSiImJ4T//+Q/V1dUdUWXlw5IlSxg0aBBRUVGceOKJrF69urOr5DZv3jx69epFbGws11xzzXH/u2K7QBKREBEJB4KBYBEJFxFvZ3J/Aa4VkcEiEg/cByzqqHpGRUXxq1/9CrB+oRvPbxcXF0dycrI7fOLi4ggLa9grISIigrq6OnJycjhy5AhOp7Ojqq9cPv74Y2bNmsXLL79MaWkpn332Gf379w/4cR0OR7NlVqxYwWOPPcbKlSvJzc1l9+7d/P73vw943TqT7QIJK1gqse6W/cL19X0ikiEiZSKSAWCM+QiYC/wL2AvkAh3607rkkktISEhg165dzd6V8yUhIYG6ujrWr1/PypUrWbNmDTt37jwm4FRg/P73v+eBBx5gzJgxBAUFkZ6eTnp6OgCffvopffr0Ye7cuaSkpNC7d2+WL1/O3//+d7KysujZsyePPvqoe19Op5PHHnuME088kcTERC699FIKCgoA+P777xERXnzxRTIyMjjnnHOardsrr7zCtddey5AhQ0hISOD+++9n0aJFPsv/5S9/oV+/fiQmJvLwww+TmZnJP//5TwC+/PJLd/NA7969mTlzJjU1Ne5tRYSnnnqK/v37k5SUxJ133un+D3Lnzp1MnDiRuLg4kpKSuOyyy1r9fW4p2wWSMeZBY4w0ej1ojNlrjIk2xuz1KPsnY0yqMSbWGHO1MaZDz2dTU1P50Y9+BMA777yD0+lk7969rT7TiYqKIiUlhbi4OIwxbN++nV27dgWiyp2uvv0s0K+WqD87zc/PZ8CAAfTp04eZM2dSWVnpLnPw4EGqqqo4cOAAs2fPZsaMGSxevJgNGzawevVqHn74Yfbs2QPA/PnzWb58OatWrSIvL4+EhAR++9vfNjjmqlWr2Lp1KytWrGi2fps3b2bYsP+2SgwbNoxDhw5x9OixrRJbtmzhhhtu4LXXXuOHH36guLiYAwcOuNcHBwczb948jhw5wtq1a1m5ciXPPPNMg30sW7aMnJwcNm7cyLvvvstLL1ldAe+//34mTZpEYWEh+/fv58Ybb2zBd7dtbBdIXUlkZCQ/+clPCAoKYuXKlUybNo1LLrnE/YNsreDgYMLDw0lOTua7775j7969lJeXu0/vjTFUVVVRUFDAvn372LVrF3V1df78SN3KoUOHqK2t5Z133mH16tV89dVX/Oc//2HOnDnuMj169ODee++lR48eXH755Rw5coSbb76ZmJgYhgwZwuDBg9m0aRMAzz33HI888gh9+vQhLCyMBx98kHfeeafB5dmDDz5IVFQUERERzdavrKysQfeR+q+93f195513mDJlCuPHjyc0NJTZs2c3CObTTjuNMWPGEBISQmZmJr/61a9YtWpVg33MmjWLnj17kpGRwS233MIbb7zh/h7k5uaSl5dHeHg448ePb8m3t0268l22ThcREUFqaioTJkxg1apVbNu2DbB+OaZPn97gzltrBAUFkZSUxObNm92/VJ535urf19XVUVZWximnnEJQUNf4v8VOl6L1oXDjjTfSu3dvAG677TbmzJnDI488AkBiYiLBwcENyqempjbYR1lZGQC5ublcfPHFDX4WwcHBHDp0yP2+b9++La5fdHQ0JSUl7vf1XzfuTgKQl5fXYN+RkZEkJia632/fvp3bbruNnJwcKioqcDgcnHbaaQ324bl9v379yMvLA2Du3Lncf//9jB49moSEBG6//XauueaaFn+O1ugav8U2FRISQlRUFNdddx2nnXYaM2fOJCMjgyNHjvDFF1+0e9/1d+SSk5NJTEwkISGhQUN5amoqBw4cYNu2bS1qJFUNJSQk0KdPnwZnEi293POmb9++fPjhhxQVFblfVVVV7jap1u5/yJAh7rMvgE2bNpGamtogaOr17t2b/fv3u99XVlY2uLT7zW9+w8knn8yOHTsoKSnh0UcfPeY/h3379rm/3rt3L2lpaQD06tWLhQsXkpeXx4IFC7jhhhvYuXNniz9Ha2ggtVNCQgJpaWksWLCA6dOnM2XKFAA++OADAKqrq722Ke3fv79Bo2JzRMTrWVBKSgp79+7ls88+44cffmjjp+i+rr76aubPn8/hw4cpLCxk3rx5TJ48uU37+vWvf829995Lbq41Wmt+fj7vvvtuk9tkZmb6bKieNm0aL774Ilu2bKGoqIg5c+Ywffp0r2X/53/+h/fff581a9ZQU1PDgw8+2CBwSktLiY2NJTo6mm3btvHss88es48nnniCwsJC9u3bx5///Gd34/Xbb7/tDruEhASfv4v+oIHUTj179mwQLD/96U8REVatWsUjjzzChAkTOO+887jlllt4+umneeutt7jmmmu46KKLuOWWW9p9q19ESE5OJiIigg0bNjT4X0417/777+f0008nKyuLQYMGMWLECO6999427evmm2/mwgsvZNKkScTExDBmzBjWrVvns3xNTQ1Hjx5lzJgxXtf/+Mc/5q677iI7O5uMjAz69evHQw95f358yJAhzJ8/n8svv5zevXsTHR1NSkqKu6vJk08+yeuvv05MTAwzZszweqds6tSpnHbaaQwfPpwLLriAa6+9FoD169dzxhlnEB0dzYUXXsif//zngHWNEDtd03eEUaNGGX8O8l9aWsrnn3/e4CHbmTNntviS7ZZbbuEXv/iFX+ricDg4evQoZ5xxhtfTemUv//73v3n66afdjcf+VFZWRnx8PDt27OCEE05otryIsGPHDgYMGOD3urj2v8EYM6q5ctqo3U5RUVEEBQVhjHG3D/zv//4vX375JUOGDOGee+4hOjqar7/+mj179nDgwAGGDBlCYmIiv/vd73j66ac544wzGDhwYLvrEhISQlxcHOvXr+eUU04hLS2tXW0iKrDGjx/v1ztW77//Pueeey7GGO644w5OOeUUutqUXxpI7RQUFER8fDxVVVXuuzDjxo1j5cqVREdHuwOh/i6Op/Xr17N06VKuv/56brrpJqZOndrua/Pw8HCCg4P5+uuv2bt3L8OGDSMyMrJd+1Rdw7vvvstVV12FMYZRo0axZMmSLvcfkl6y+cHu3bvZuXNnqy+TKisrmTVrFmvWrAGscKsfASA9PZ2pU6e2qEevLyUlJTgcDkaMGAHA4cOHiY+PJzEx8ZjHWJQKpJZesmkg+UFBQQFffvllmwZrM8bwj3/8g6eeeqpBfxWwAur//u//GD16dJvrVlVVRXFxMUFBQfTo0YPa2lpEhFGjRmk7k+owGkg+BCKQqqur+eSTT0hJSaGiogKn00l0dHSr9+NwOCgvL+fw4cMsX76cN998k7i4OBYvXuy+5MvLy8PhcNC3b982n45XVVVRWVnJuHHjWtRjWKn20kbtDhQWFkZkZCTl5eWUl5cDtCmQ6hul4+LiuO2229i7dy9r165l2rRpnHvuuezfv9999y4tLY2JEycyderUVt8ZCQ8Pp6qqiq+//pp+/fpRUVFB7969NZxUp9MzJD/ZvHkzO3bsYOzYsWzZsoWIiAh69OjRrn0WFxdz4403smXLFveysLAwwsPDKS7+72CZ9ZdeKSkpnHvuuZx//vkNegf7UlRURG1tLWA9anDGGWdo25IKCL1k8yFQgVRQUEBpaSn9+vVj165d7N69m549e7Z7v8YYvvvuO1atWkV0dDSTJ08mOjqarVu38sEHH/Dhhx+6z8o8ZWRkMGnSJK677roWPVNXVFREVFQUgwYNIiIigtDQ0HbXXal6Gkg+BCqQPJWUlPD555+TkpJCaWkpVVVVJCQktPlh26bU1NRQVFSEiLB582ZWrFjB2rVr3Q98TpkyhQceeKBF7U2eZ0yhoaEkJibSv39/rw9zKtUaGkg+dEQgGWP49NNPMcZgjOGEE05g+/bt9OjRwz1+dlFRETU1NaSkpPj9+A6Hg7Vr13L33XdTXV3Nz3/+c37729+2ql3L4XBQUVFBTU2N+ynvxup/d7paXxfV8TSQfOiIQALYsWMH27dvZ8KECcTGxlJRUcHWrVs5dOgQIuIe1qKoqChgg/yvXbuWW2+9FYfDQVhYGNnZ2VxyySWMGDGixSFSVVVFSUkJgwYNIjU1FRGhuLiYw4cPk5+fj9PppGfPnqSmppKSktLudjN1fNJA8qGjAqmqqorq6uoGA2wZY8jLy6O0tJSBAwdSU1PDZ599Rnx8fEAu58AauvSll17C8zNnZWXx+OOPt3hsHofDQXFxcYPB4MLDw4mMjEREqKqqoqKiAhEhKSmJyMhIUlJStJ+TctNA8qGjAqml9u3bx9dff01KSkqDx0ZKSkooLy/3+shJW+Tl5fHuu++yfPlyjh49SmJiIk888QShoaGUlJSQkZHhPgOq53Q6W/UoS/2UTw6Hg6qqKsaOHev1Uk91PxpIPtgtkIwx7Nixgx07dpCQkEBwcDClpaWEhYURFhZGZWVlm/o0+VJRUcEdd9zBl19+ecy6sLAwkpKSCAsL4/Dhw1RWVnL66aczceJE0tLS3G1jGzduJDQ0lJiYGGJjY+nZsyfZ2dnugfLBOkMsLy9n6NChFBUVUVFR4R7F0HPERdU9aCD5YLdAqnfw4EG2bduGMYa4uDiGDh1KSUkJ69ev93vDd01NDX/4wx9YuXKle86v3NxcCgsL27XfjIwMJkyYwKmnnsq4ceNwOp2Ul5cTGhpKjx493A3lJ5xwAieddJJ7aFigwWgJ6vijgeSDXQPJG6fTyapVqwgPD29xvyCn00lJSYn7mbX6s66WKCsro7CwkIqKCvdZzKeffsqGDRsoKCigqqqK0aNHM2HCBPeZXElJCbt372bZsmUNnsVLTU3l5ptv5vzzz28QNMYYjhw5QnBwsLsNqn6M57CwMFJTUzn55JO7zBjhqmU0kHzoSoEE1sDx27ZtIykpyb3M4XAc0wjudDopLCykrq6Ofv36ER8fT2lpKbt37w5I14LGHA4HOTk5bNq0iU8//ZQdO3YAkJ6eTnZ2NiLCkSNHGDduHD/60Y9wOp3U1dXhdDrp0aMHwcHBOBwOCgoKyMrK4sQTTwx4nVXH0UDyoasFUnV1NWvWrMHhcBAbG8vRo0fdZzyJiYk4nU5KS0upqakhMzOTzMxM9zNpTqeTjRs3Ulxc7NeuBfV323ydedXV1fHee+/x/PPPk5+ff8z6M888k9/97ndeG+zr6urIz89n9OjRbRo9QdmTBpIPXS2QwAqlXbt2sW/fPk4++WR69+7Nzp07yc3NJSQkhLS0NDIzM4mKijpm26qqKj7//HPCwsLa9fBsbW0tFRUVVFdXExISgjGGuro696VVSEgIsbGxDS616urq2LRpE2vXriUiIoKQkBD3lNVhYWFMnz6ds846C2MM6enp7h7h1dXVFBUVkZ6eTmZmJtHR0e7wM8ZQXV1NZWWlOxgjIiLcl3/KnjSQfOiKgVSvrq6uwVlJeXk5ERERzba3FBcX88UXXxAdHU14eHiLj1dWVkZFRQVgPXxbPwVT/TN6hYWFVFVVERISwtGjR90TDMTExPg8zpEjR/jjH//Ixx9/3GB5UFAQgwcP5qKLLmLq1KnuetfU1CAiREZG4nQ6qa2tPWZyTGMMYWFhDBo0iF69emkw2ZAGkg9dOZDao6CggHXr1hEeHk5MTIzXP1rPZ9mMMSQnJ5ORkUFsbGyLgqyqqor8/Hz27t3rnl3VGENSUtIxoblhwwZeeOEFioqKMMbw/fffu+eWGz9+PPfff7+7Y6UxBofD4Z5+x1sA19TUUFhYSFxcnPtMrf7Rnfj4eOLi4nx+bhV4Gkg+dNdAAuuMY+fOnRw+fJjQ0FDi4uLcs+Hm5+eTkpJCnz59AGs8p/b0f6qrq6Ouro7vv/+eHTt2uB+fcTqdXgOuoqKCTz75hD/+8Y+UlpYSFxfHLbfcwuTJk1sVIpWVlTgcjgbdCOrnxktLS+Pkk0/WIVY6gQaSD905kOqVlZWxfft293N1xhgGDBjAiSeeGJDb7Xv37mXfvn3069ePHj168N1331FeXk5ycvIxYXPo0CFmz57tns9s5MiR3HHHHWRlZbW7HgUFBQQFBZGVlUVqaqoOsdKBNJB80ED6r5KSEkSE8PDwDn0o1ul0smXLFvbt2+e1S4Ixhg8//JA//elPFBUVERQUxMCBA6msrKSiosL93FxsbCypqalkZWUxatQozj777GYDtaamhuLiYkSE6OhoIiIiiI+PJzY21t1YHxsb2+K+W6plumwgiUhP4EVgEnAE+J0x5nUv5cKAPwMXAz2Az4FfG2MONLV/DSR7cDqdbNq0iby8PEJDQwkNDSUqKqrBGVNJSQkLFy7krbfeOqYh25uBAwdy1VVXkZWVRUZGRpNnQE6nE4fDQW1tLdXV1e79OxwO+vfvz+DBg9v/IZVbVw6kN7Cm+L4WGA78DTjTGLO5Ubm7gCuxgqsYeB6INsZc0tT+NZDso77PUWVlJQUFBRw+fNg9O0p4eLi7rad+qJPIyEgiIiKIjo7GGENxcTH79+9ny5Yt/PWvf23QUzwoKIjevXvTv39/TjvtNMaNG9eiGVyNMRw+fJhRo0Z1SIfS7qJLBpKIRAGFwFBjzHbXsleBA8aYuxuVfRYoNcbc5Xp/AfAnY8xJTR1DA8m+Kioq3EMBHzx4EIfD0eJhgKurq1m2bBnr1q1jz5495OXl4XQ63etFhEsvvZSZM2c22x+rurqa8vJyxo8frxMf+ElXDaQRwOfGmEiPZXcAE40xUxqVHYV1yfZzoAh4AThsjLmlqWNoIHUNNTU1bN68mYMHD3rtNtCS7Q8cOMCWLVtYt24dK1asoK6ujvT0dB566CGGDx/e5Pb1c9mNGjXKr6MtdFddNZAmAG8bY3p5LJsBXGmMObtR2ThgAXAZUAd8A5xrjCnwst/rgesBMjIyTsvNzQ3YZ1D+Y4xh+/bt7Nq1i+Tk5HbdAdy2bRuzZ89m+/btBAUFceWVVzJx4kROOukkn2dBZWVl1NTUMHbsWA2lduqqgeTtDOl24GwvZ0iLgSistqZy4C5gsjHmjKaOoWdIXYsxhl27drmDxFNQUJB73KiW3CWsra1lwYIFvPLKK3j+3sfHx5Oens7IkSMZM2YMo0ePdjeuFxUVkZKSwimnnOLfD9bNdNVAqm9DGmKM2eFa9hcgz0sb0rfAvcaYd13v413bJhtjjvg6hgZS11RSUgJYIVR/d6y8vJzCwkJKSkrcXQHi4uKa7V/01Vdf8d5777Ft2zZ2797t7iFeb/z48dx3330kJSXhdDo5evQo2dnZ2qGyHbpkIAGIyBLAANdh3WX7O97vsr0MxALXABXAncBvjTFNzpCogXR8qqmp4dChQ+zcuZOqqiqCgoJa9Oye0+mkoKCAXbt28eWXX7J06VJ3T/HzzjuP7OxsTjzxRIYMGdLiMcjVsbpyIPUEXgLOB44CdxtjXne1L31ojIl2lUsEnnKVCwW+BW4zxhw7NqsHDaTjmzHGPdDcvn37GnT+jI6ObvYxlMY9xQFmzJjBFVdcwcSJE/VZuDbqsoEUaBpI3UtlZSUlJSXk5eWRl5dHXFxcs7fyjTFs3bqVf/3rX7z88ssEBwfz+OOPc9lll5GWlkZQUBC1tbXuUTnrtxERQkNDtZe3Fy0NpMDMvaOUTURERBAREUFqair9+vXjm2++aXYuPBFh8ODBDB48GIfDwauvvsqTTz5Jenq6+xm4+tEM6p8F9BQWFtZgrjodjrfl9AxJdSvV1dXk5ORQWVnZoimaamtrufbaa9myZQsJCQlce+211NXVcejQIYYNG8a4ceOOaaeqnwaqoqKCmJgYsrKyiI+P79YP8+olmw8aSKqmpoaNGzdSWFjoflSlfhA4bxN2Hjx4kAceeICNGzcesy48PJyePXsSGhqKw+HA4XDQp08fsrKymDRpEv3796esrAyA2NhY9zRQ3W2+Og0kHzSQFFhtPqWlpRQWFrrHUMrLy8MY4x6dsqqqyj3IW3BwMJ9//jkfffQRqamppKamsmbNGrZs2dLkcUaNGsXll1/O+PHj3WdONTU1jBs3jtjY2A76tJ1PA8kHDSTlS3V1Nd9//z3l5eUEBweTkJBAQkKC+yHgHTt2EBMTQ2Sku98uRUVF7h7dISEhiAi5ubls2LCBpUuXUl5eDkBCQgLnnHMO48eP59RTT0VEOPPMM7tN3yYNJB80kFRbFRcXs3HjRqqrqwGrD1NCQoLPtqHS0lKWL1/O+++/z+7du93LBwwYwH333UdGRgZpaWmkpKQc95dwGkg+aCCp9qgfPyksLIyCggK++eYbamtrG4zhXd8VICwszD2Tynfffcfq1av54IMPOHDgAFFRUdx4443uWVdOOukk+vfvf9z2c9JA8kEDSflTdXW1e6ZgwD0LTGVlJQcPHuSHH34gKCiI+Ph4QkJCKCsrY/bs2XzyyScAJCcnc+mll3LWWWeRlpZGcHAwwcHB9O7d2z3mOVizvnTl/k0aSD5oIKmOVFVVRV5eHjt37sTpdBITE0NYWBgffvghr7zyCrt27QIgKiqKIUOGEB0dzYgRIzjnnHMICgpyT8IQFBREXFwcUVFRxMbG0rdv3y7Vv0kDyQcNJNUZamtryc/PZ8+ePZSUlBAUFERUVBQbN25k8eLFrF+/vkH5AQMGMHnyZMCahDMyMpI+ffowYMAAKioqGDRoEP379++Mj9ImGkg+aCCpzlY/SsHu3bupqKggMTGR3NxcDhw4QH5+PosWLeLAAe9Dw59wwglceeWVjB07ljFjxrR4RM3OpoHkgwaSsguHw8GePXvcl22ePv74Y/Ly8tzPzZWVlbFu3ToKCqzxB2+77TZ++tOfMn78+C7RdUCfZVPK5kJCQhg4cCCZmZnU1NRQU1NDdXU1+/bt47zzzjtm3jqHw8E777zDk08+ybx584iMjCQpKYlhw4Z14qfwLw0kpTpZjx496NGjB1FRUYB1523r1q3k5uYSHBzsfoDXGMMFF1xAeXk5zz77LI899hjV1dX06tWL1NTUTv4U/qGBpJTNBAcHM2TIENLT03E6ne5hTZxOJ5s3b2by5MmUlpayePFinnjiCUpLS5k3b16zg9F1BRpIStmQiHjtvX3GGWewe/duLr/8chISEpg/fz7PPfccPXv2ZM6cOV2+Y6UGklJdSFBQEAMGDCA0NBQRQUR46qmn+MMf/kBmZiYzZszo7Cq2iwaSUl1QRkYGcXFxDBw4kJKSEhYtWsQNN9zAwYMHmTVrVpcde6nrdPVUSjVQH0gLFizg6quvxuFw8MADDzBp0iReeOEF9u7dS11dXWdXs1U0kJTq4kJDQ3nxxReZN28eYWFhrFq1ihkzZnD66aezcOHCLhVK2jFSqePIRx99xAsvvMDatWvJy8sDYMqUKYwYMYK0tDQuvvhiUlJSOrxe2lPbBw0kdTwzxrBz50527tzJkiVLeO211xqcIQ0dOpQ5c+YwZcqUDn04VwPJBw0k1R3UP8y7YsUKPvnkEwoKClizZg1FRUWICFdffTWPPvooKSkpHdJVQAPJBw0k1Z0YYyguLmbv3r1s376d119/neXLl2OMYdiwYTz66KNkZ2c3O1dde2kg+aCBpLqrqqoqCgoKePPNN3nooYcoLi4mIiKC6dOnc91115GRkUFCQkJABoJraSDpXTaluonw8HDS0tK46aabePvttznzzDOprKzk2Wef5ayzzmLatGmsXr2ampqaTqujdoxUqpsJDg7mnHPOoVevXrz55pssXbqUrVu38uGHH3LkyBGefPJJTj31VGJjYzt8VEoNJKW6oeDgYIYOHUpKSgrZ2dls27aNe++9l/Xr17Nw4UKuuOIKQkNDOemkk0hPT++wZ+Q0kJTqpkSE1NRUevbsSVZWFlVVVdx1110sXryYNWvWMHToULKzs8nOzmbw4MH06NEj8HXSRm2lFEBlZSWzZs1iwYIFDdqRRo4cyd13383FF1/sdarxlrBFo7aIBHm+WrhNTxFZJiLlIpIrIlc0UXakiHwmImUickhEbvZf7ZXqXiIiIrj99tt5++23WbRoEVdddRWRkZFs3LiRa665hgULFgT8MRS/B5IrJNZApAefAAAWO0lEQVSKSDlQ63o5XP+2xNNADZAKXAk8KyJDvBwnCfgIWAAkAgOAf7T/EyjVffXt25fExERSUlL4xS9+wZIlS5g4cSJlZWXcdNNNzJ07l8rKyoAd3++XbCLyDfA+8CpQ4bnOGJPbzLZRQCEw1Biz3bXsVeCAMebuRmUfBfoaY65qTf30kk2pplVVVVFcXEx1dTWFhYXk5+ezcOFC3nrrLSIjI3nmmWe48MILWzX9d2cO8t8PuNe0LemyAEd9GLlsAiZ6KTsG+EZE1mCdHa0DfmuM2duG4yqlXMLDw93D4WZkZGCM4ZRTTuHgwYN89tlnPPbYY6Snp5Odne33TpSBaENaBkxq47bRQEmjZcVAjJeyfYBfAjcDGcAe4A1vOxWR60UkR0Ry8vPz21g1pbqn+rtxS5YsITk5mW3btvHqq6+yf/9+vx8rEIEUDiwTkX+IyF88Xy3YtgyIbbQsFij1UrYSWGaMWW+MqQIeAs4UkbjGBY0xzxtjRhljRiUnJ7fy4yilAHr37s3cuXMBWLZsGRs2bKCqqsqvxwhEIG0BHgc+B3Y1ejVnOxAiIgM9lg0DNnsp+zXgeVnYvfovKNUJLrvsMoYNG0ZpaSlLly7l8OHDft2/7fohicgSrHC5DhgO/B040xizuVG5c4C/AtlYgTUXGGWMmdDU/rVRW6n2Wbx4sbtLwJo1a1o0UWWn9kMSkbNF5CURWeH6N7sVm98ARACHsdqEfmOM2SwiE0SkrL6QMeYT4B7gb66yAwCffZaUUv5xySWXMHz4cCoqKnjxxRf9uu9A9EO6DngLOAgsBX4A3hCRFs3PYowpMMZcZIyJMsZkGGNedy1fbYyJblT2WWNMujEmwRgzxRizz88fRynVSGRkJDfffDNRUVHEx8f7dd+BuO1/F3C+MWZT/QIReRPr8mphAI6nlOpgF198MaGhoZx55pl+3W8gLtkSsRq2PX0H9AzAsZRSnSAuLo7+/fv7f8fGGL++gHeB+UCk630U8BTwvr+P1cb6GV+vBQsWmHoLFizwWc76tv3XyJEjfZabMWOGu1xOTk6T+8zJyXGXnTFjhs9yI0eObHB8/Uz6mTrjMxUUFLT4MwE5pgV/n4E4Q/o11q36YhE5BBS53v8qAMdSSnWS1jw60lIBu+0vIn2B3kCeMcb/XTrbSG/7K9XxOvRZNhER12kZHsOMHHC93MuMMU5/HE8pdXzy1122Yv77yIcDjuk1La5l/p/OQCl13PBXIHmOV3SCn/aplOpm/BJIxqNDomk05pGIRABOY0y1P46llDp+BaKn9pMiMtr19QVAAVAoIlP8fSyl1PElELf9rwS+dX39APAL4ELg0QAcSyl1HAnEoyORxpgKEUkE+htj/gogIv0CcCyl1HEkEIG0XUSuxHr6/mNwD8gfuJHBlVLHhUAE0g3An7FmDrnWtexH6IwgSqlm+D2QjDHrgTMbLXsNeM3fx1JKHV/81VP7LGPMZ66vz/FVzliDqimllFf+OkN6Bhjq+trXEHIGCMB4BUqp44W/OkYO9fhae2orpdokEB0jh7ue9Pdc1ldEmh8JXCnVrQWiY+RioEejZaFYU2srpZRPgQikDGPMbs8FxphdQGYAjqWUOo4EIpD2i8hIzwWu93kBOJZS6jgSiI6R84B3RWQu1my1JwJ3AI8E4FhKqeNIIDpGLhSRIqxe2n2BfcDtxph3/H0spdTxJRBnSBhj3gbeDsS+lVLHr0Dc9hcRmSEiK0Xka9eys0TkUn8fSyl1fAlEo/ZsrMu1hUCGa9l+YFYAjqWUOo4EIpCmA5ONMUv472D/e9DHRpRSzQhEIAUDZa6v6wMp2mOZUkp5FYhA+hD4k4iEgdWmBDwMvB+AYymljiOBCKRbsWasLQbisM6M+tHCNiQR6Skiy0SkXERyReSKZsqHishWEbHN7LhKqbbx621/19lQEvBzoCdWEO0zxhxsxW6exhptMhUYDvxNRDYZYzb7KH8nkA/EtLniSilb8OsZkms67W+w5mE7bIxZ35owEpEo4GfA/caYMmPMv4H3gKt8lD8Ba1aTP7S/9kqpzhaIS7b/AFlt3DYLcBhjtnss20TDmXE9zQfuQScQUOq4EIie2p8CH4nIIqzHRurvtGGMeamZbaOBkkbLivFyOSYiFwPBxphlInJ2UzsVkeuB6wEyMjKaKqqU6kSBCKRxWP2OJjZaboDmAqkMiG20LBYo9VzgurSbC/y0JRUyxjwPPA8watQo00xxpVQn8VsgiUgkcB9WqGwEHjXGVLdyN9uBEBEZaIzZ4Vo2DGjcoD0Qa3yl1VY7OqFAnIgcBMYYY75v04dQSnUqf7YhPQ1MAbZiNUw/2dodGGPKgaXAbBGJEpFxwFSOHW3yW6yRBIa7XtcBh1xf72vrB1BKdS5/BtKPgUnGmLuAnwCT27ifG4AI4DDwBvAbY8xmEZkgImUAxhiHMeZg/QsowLqzd9AYU9f+j6KU6gz+bEOKMsb8AGCM2ScicW3ZiTGmALjIy/LVWI3e3rb5FOjTluMppezDn4EUIiLZgPh4rxNFKqWa5M9AOkzDu2hHG73XiSKVUk3yWyAZYzL9tS+lVPcUiJ7aSinVJhpISinb0EBSStmGBpJSyjY0kJRStqGBpJSyDQ0kpZRtaCAppWxDA0kpZRsaSEop29BAUkrZhgaSUso2NJCUUrahgaSUsg0NJKWUbWggKaVsQwNJKWUbGkhKKdvQQFJK2YYGklLKNjSQlFK2oYGklLINDSSllG1oICmlbEMDSSllGxpISinb0EBSStmGBpJSyjZsF0gi0lNElolIuYjkisgVPsrdKSLfikipiOwRkTs7uq5KKf8K6ewKePE0UAOkAsOBv4nIJmPM5kblBJgGfA2cCPxDRPYZY5Z0aG2VUn5jqzMkEYkCfgbcb4wpM8b8G3gPuKpxWWPMXGPMRmOMwxjzHfAuMK5ja6yU8idbBRKQBTiMMds9lm0ChjS1kYgIMAFofBZVv/56EckRkZz8/Hy/VVYp5V92C6RooKTRsmIgppntHsT6LC97W2mMed4YM8oYMyo5ObndlVRKBYbd2pDKgNhGy2KBUl8biMhMrLakCcaY6gDWTSkVYHY7Q9oOhIjIQI9lw/B9KXYNcDdwrjFmfwfUTykVQLYKJGNMObAUmC0iUSIyDpgKvNq4rIhcCTwKnG+M2d2xNVVKBYKtAsnlBiACOAy8AfzGGLNZRCaISJlHuTlAIrBeRMpcr+c6ob5KKT+xWxsSxpgC4CIvy1djNXrXvz+hI+ullAo8O54hKaW6KQ0kpZRtaCAppWxDA0kpZRsaSEop29BAUkrZhgaSUso2NJCUUrahgaSUsg0NJKWUbWggKaVsQwNJKWUbGkhKKdvQQFJK2YYGklLKNjSQlFK2oYGklLINDSSllG1oICmlbEMDSSllGxpISinb0EBSStmGBpJSyjY0kJRStqGBpJSyDQ0kpZRtaCAppWxDA0kpZRsaSEop29BAUkrZhu0CSUR6isgyESkXkVwRucJHORGRx0XkqOv1uIhIR9dXKeU/IZ1dAS+eBmqAVGA48DcR2WSM2dyo3PXARcAwwAAfA3uA5zqwrkopP7LVGZKIRAE/A+43xpQZY/4NvAdc5aX4L4E/GmP2G2MOAH8EpndYZZVSfmerQAKyAIcxZrvHsk3AEC9lh7jWNVdOKdVF2C2QooGSRsuKgRgfZYsblYv21o4kIteLSI6I5OTn5/utskop/7JbIJUBsY2WxQKlLSgbC5QZY0zjgsaY540xo4wxo5KTk/1WWaWUf9ktkLYDISIy0GPZMKBxgzauZcNaUE4p1UXYKpCMMeXAUmC2iESJyDhgKvCql+J/AW4TkXQRSQNuBxZ1WGWVUn5nq0ByuQGIAA4DbwC/McZsFpEJIlLmUW4B8D7wDfAt8DfXMqVUF2W7fkjGmAKs/kWNl6/Gasiuf2+Au1wvpdRxwI5nSEqpbkoDSSllGxpISinb0EBSStmGeOlHeFwTkXwgtwVFk4AjAa5Oe9m9jlq/9rF7/aDldexnjGm2V3K3C6SWEpEcY8yozq5HU+xeR61f+9i9fuD/Ouolm1LKNjSQlFK2oYHk2/OdXYEWsHsdtX7tY/f6gZ/rqG1ISinb0DMkpZRtaCAppWyjWweS3Wc4aUX97hSRb0WkVET2iMidga5ba+rnUT5URLaKyH671U9ERorIZyJSJiKHRORmO9VRRMJE5DlX3QpE5H0RSQ9w3Wa6RlqtFpFFzZS9VUQOikiJiLwkImFtOWa3DiQaznByJfCsiHgbl9tzhpNTgSnAr2xUPwGmAQnAj4GZInK5jepX706gI8cQblH9RCQJ+Ahr+JpEYADwDzvVEbgZGIv1+5cGFALzA1y3PGAO8FJThUTkR8DdwLlAP6A/8FCbjmiM6ZYvIArrFyHLY9mrwGNeyq4Brvd4fy3whV3q52Xbp4D5dqofcAKwFfgJsN9mP99HgVc78vevDXV8Fpjr8f4C4LsOquccYFET618HHvV4fy5wsC3H6s5nSHaf4aQ19XNzXUpOIPDD+ba2fvOBe4DKANerXmvqNwYoEJE1InLYdTmUYbM6vgiME5E0EYnEOpv6sAPq2BLe/j5SRSSxtTvqzoEUkBlO/Kg19fP0INbP9eUA1MlTi+snIhcDwcaYZQGuk6fWfP/6YM3zdzOQgTXh6BsBrZ2lNXXcAewDDri2GQTMDmjtWs7b3wc0/7t6jO4cSAGZ4cSPWlM/wGqExGpLusAYUx3AukEL6+ea/HMucFOA69NYa75/lcAyY8x6Y0wVVvvHmSISZ6M6Pg2EYbVxRWGNPW+XMyRvfx/QxO+qL905kOw+w0lr6oeIXIOrYdEY0xF3sVpav4FAJrBaRA5i/SH1dt2RybRB/QC+xpqOvV5H9RZuTR2HY7XjFLj+s5kPjHY1yHc2b38fh4wxR1u9p45uyLPTC1iCdWoeBYzDOtUc4qXcr7EaZNOx7nBsBn5to/pdCRwEBtnt+4c1bnsvj9clWHdvemFdxtnh+3cO1l2r4UAPYB6w2i7fQ1e5l4G/AnGuOt4DHAhw3UKAcOAPWI3t4UCIl3I/dv3+DQbigU9owc0Xr8fsyF9gu72AnsByoBzYC1zhWj4B65KsvpxgXXYUuF5zcT12Y5P67QFqsU6d61/P2aV+jbY5mw64y9ba+gG/wWqfKcSazaavneqIdan2GtZsPEXAv4HRAa7bg1hni56vB7Ha2cqADI+ytwGHsNq3XgbC2nJMfZZNKWUb3bkNSSllMxpISinb0EBSStmGBpJSyjY0kJRStqGBpJSyDQ0kdVwSkbM9x10Ske9F5LzOrJNqngaS6hCuQKh0DYB2UEQWiUh0Z9dL2YsGkupIU4wx0ViPaIwAftfJ9VE2o4GkOpwx5iCwAiuY6odnfVJE9rqGaH1ORCLqy4vIVBH5yjU86i4R+bFr+dWuIXFLRWS3iHTEKJ4qgDSQVIcTkT5YI0fudC16DGuwsuFYw8emAw+4yo4G/oI1/G08cBbwvWu7w8BkrOEurgbmicjIDvkQKiA0kFRHWi4ipVgDjR0Gfu8a5O564FZjDa1RijWkbP2Y4NcCLxljPjbGOI0xB4wx2wCMMX8zxuwyllVY42BP6PBPpfxGA0l1pIuMMTFYT/yfDCQByUAksEFEikSkCGvA/WTXNn2BXd52JiI/EZEvXLNwFAE/de1TdVEaSKrDuc5mFgFPAkewRmwcYoyJd73iXI3fYJ1Nndh4H65pdv7q2keqMSYe+DvWUDGqi9JAUp3l/wHnA6cAC7Haf1IARCTdNbUOWIPbXy0i54pIkGvdyUAo1pCu+YBDRH4CTOrwT6H8SgNJdQpjTD5WY/UDwCysBu4vRKQE+Cdwkqvcl7garLFGU1wF9HO1Nd0EvIU1qNoVwHsd/DGUn+kAbUop29AzJKWUbWggKaVsQwNJKWUbGkhKKdvQQFJK2YYGklLKNpoMJNdT2C+KSK7rieqvXB3QPMucKyLbRKRCRP4lIv081j0pIjtc224TkWmNtn1eRL4TEaeITG+usiIyXEQ2uI61QUSGe6zLdh2/WES+b8G+7hSRb1112yMid3qsy3CN2+P5MiJyu499DRWRFSJyRESO6UchIoNE5BNX3XaKyMVN1GuMiHzsehwiX0TeFpHeHuvDXE/DH3KVeV9E0n3sq8mfX3PHUqqjNXeGFILVdX8i1hS+9wFviWtOdrHmFV8K3I81A2cO8KbH9uXAFNe2vwT+LCJneqzfBNwAbGyuoiISCrwLLAYSgFeAd13L64/1EtZT4S0hwDTXvn4MzBSRywGMMXuNMdH1L6zexE6sRxW8qcXqoHetl3qHuOr9Adb36HpgsYhk+dhXAvA8kAn0A0qxZgKtdzMwFjgVa1rvQqx53r1p8ufXgmMp1bHaML3u18DPXF9fD6zxWBeF9VzSyT62fQ+43cvyfwPTmznuJKypjsVj2V7gx43KnQd834bP9RQw38e63wP/asE+Bljf0gbLhmJNO+xZ738AD7ewXiOBUo/3zwJzPd5fAHzXlp9fc8fSl746+tWqNiQRScUat2aza9EQrLOc+nArx3oye4iXbSOA0z22ba0hwNfGGM9Loq+9Hau1XENgTPBWN9e6aVhnZP4iWEFVf4wiERnvo+xZjer1IjBORNJEJBK4EvjQY1/PiMgzXg967M+vuWMp1aFCWlpQRHoArwGvGNd4NEA01sONnoqBGC+7eA4rvFa0oZ71xypu4bFa60Gsy1dvlyvjgVTgnTbu+zussX/uFJF5QDbWJdS/6gsY60n1Y4jIqVjPek31WLwD6zLsAFAHfAPM9NjXDT725e3n19yxlOpQLTpDEpEg4FWgBo9ffqxLkdhGxWOx2iI8t38C64zg0kZnOE0d07NBOaOlx/Kxr3s89vVco3Uzsc6ALjDGVHvZ/JfAX40xZS2pd2PGmFrgIqxLq4PA7VjtTfub2k5EBmCd+dxsjFntsepprKfcE7EukZficYbkY1++fn7NHUupjtXcNR3W5cXLWP+jRzRadz3wucf7KKACjzYk4CHgWyCxiWO0tA1pPw3bYnJpRxsScI1rn/19rI/AOgs7p4X7O6YNyUe5NcCvmljfD2uY1l97WfctMNXjfTxggKTW/vyaO5a+9NXRr+YLWJdaXwDRXtYlu/5gfwaEA48DX3is/x3WJUYvH/sOdW33OTDD9XVQE2Vzse4yhWH9T58LhLrWB7m2/4lreXj9Oh/7uxLrjGVQE2WucP2xiq8yrnLiOt5gVziEA2Ee6091LYsE7gD2eK5vtK90rHa4O3ysfxnrbl8c0AO4BzjQxp9fk8fSl746+tX0Sut/TwNUYV0y1b+u9ChzHrAN6+7ap0CmxzoDVDfa9h6P9Z+6yni+zm6iPiOADa5jbQRGeKw728u+Pm1iX3uwbtd71u25RmVW0IK7YVi3zRsf+3uP9U9g3Z4vw7o0GtBo+zJgguvr37u296xXmUfZRKy2oMNAEdbZ5WiP9c/Vf47mfn7NHUtf+urol46HpJSyDX10RCllGxpISinb0EBSStmGBpJSyjY0kJRStqGBpJSyDQ0kpZRtaCAppWxDA0kpZRv/H2MfMTetxsnPAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 288x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Hyperparameter setup\n",
    "seed = 1210\n",
    "word_len = 6\n",
    "max_mis = 1\n",
    "nfolds = 5\n",
    "\n",
    "svm_dir = os.path.join(models_dir, \"StrongEnhancerVsSilencer\")\n",
    "if not os.path.exists(svm_dir):\n",
    "    os.makedirs(svm_dir)\n",
    "\n",
    "# Fit the SVM\n",
    "svm_prefix = os.path.join(svm_dir, f\"gkmsvm_{word_len}_{word_len}_{max_mis}\")\n",
    "fig_list, xaxis, svm_tpr, svm_prec, svm_f1, svm_scores = gkmsvm.train_with_cv(positives_fasta, negatives_fasta, svm_prefix, num_folds=nfolds, word_len=word_len, info_pos=word_len, max_mis=max_mis, seed=seed)\n",
    "# Score k-mers for motif analysis\n",
    "kmer_scores = gkmsvm.score_all_kmers(word_len, word_len, max_mis, svm_prefix, os.path.join(svm_dir, f\"all{word_len}mers\"))\n",
    "# Flip the sign for the scores to generate PWMs for the negatives too\n",
    "kmer_scores = -kmer_scores\n",
    "kmer_scores.to_csv(os.path.join(svm_dir, f\"all{word_len}mersPredNeg.out\"), sep=\"\\t\", header=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fit logistic regression model using CRX occupancy, TF occupancies, and entropy\n",
    "For the TF occupancies, do grid search for the regularization hyperparameter (`C`) since there are multiple features. For the other models, there is only one feature."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimal regularization strength (C): 1.0e-02\n"
     ]
    }
   ],
   "source": [
    "cv = StratifiedKFold(n_splits=nfolds, shuffle=True, random_state=seed)\n",
    "# CRX occupancy\n",
    "crx_clf = LogisticRegression()\n",
    "crx_clf, crx_tpr_list, crx_prec_list, crx_f1_list = modeling.train_estimate_variance(crx_clf, cv, wt_occupancy_df.loc[silencer_modeling_mask, \"CRX\"], labels_with_silencer, xaxis, positive_cutoff=0)\n",
    "\n",
    "# Information content\n",
    "entropy_clf = LogisticRegression()\n",
    "entropy_clf, entropy_tpr_list, entropy_prec_list, entropy_f1_list = modeling.train_estimate_variance(entropy_clf, cv, wt_entropy_df.loc[silencer_modeling_mask, \"entropy\"], labels_with_silencer, xaxis, positive_cutoff=0)\n",
    "\n",
    "# 8 TF occupancies\n",
    "occ_clf = LogisticRegression()\n",
    "param_grid = {\"C\": np.logspace(-4, 4, 9)}\n",
    "np.random.seed(seed)\n",
    "occ_clf, occ_tpr_list, occ_prec_list = modeling.grid_search_hyperparams(occ_clf, nfolds, param_grid, \"f1\", wt_occupancy_df[silencer_modeling_mask], labels_with_silencer, xaxis, positive_cutoff=0)\n",
    "c_opt = occ_clf.get_params()[\"C\"]\n",
    "print(f\"Optimal regularization strength (C): {c_opt:1.1e}\")\n",
    "\n",
    "# Save the 8 TF occupancy regressor to file\n",
    "pickle.dump(occ_clf, open(os.path.join(svm_dir, \"logit_model.pkl\"), \"wb\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Write cross-validated model performances to file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Helper function\n",
    "def write_folds(data, filename, out_dir=svm_dir):\n",
    "    pd.DataFrame(data).to_csv(os.path.join(out_dir, filename), sep=\"\\t\", index=False, header=False)\n",
    "\n",
    "write_folds(svm_tpr, \"svmTprCv.txt\")\n",
    "write_folds(svm_prec, \"svmPrecCv.txt\")\n",
    "write_folds(crx_tpr_list, \"crxTprCv.txt\")\n",
    "write_folds(crx_prec_list, \"crxPrecCv.txt\")\n",
    "write_folds(entropy_tpr_list, \"entropyTprCv.txt\")\n",
    "write_folds(entropy_prec_list, \"entropyPrecCv.txt\")\n",
    "write_folds(occ_tpr_list, \"occupancyTprCv.txt\")\n",
    "write_folds(occ_prec_list, \"occupancyPrecCv.txt\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fit strong enhancer vs. inactive logistic model for information content"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "inactive_dir = os.path.join(models_dir, \"StrongEnhancerVsInactive\")\n",
    "if not os.path.exists(inactive_dir):\n",
    "    os.makedirs(inactive_dir)\n",
    "\n",
    "inactive_entropy_clf = LogisticRegression()\n",
    "inactive_entropy_clf, inactive_entropy_tpr_list, inactive_entropy_prec_list, inactive_entropy_f1_list = modeling.train_estimate_variance(inactive_entropy_clf, cv, wt_entropy_df.loc[inactive_modeling_mask, \"entropy\"], labels_with_inactive, xaxis, positive_cutoff=0)\n",
    "write_folds(inactive_entropy_tpr_list, \"entropyTprCv.txt\", out_dir=inactive_dir)\n",
    "write_folds(inactive_entropy_prec_list, \"entropyPrecCv.txt\", out_dir=inactive_dir)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Use of the White 2013 data as a test set\n",
    "For now, just use a test set for the SVM and TF occupancies.\n",
    "\n",
    "### Read in the sequences (Dataset S1) and format appropriately"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "white_data_dir = os.path.join(\"Data\", \"Downloaded\", \"CrxMpraLibraries\")\n",
    "white_seqs = pd.read_csv(os.path.join(white_data_dir, \"white2013Sequences.txt\"), sep=\"\\t\", header=None, usecols=[0, 8], index_col=0, squeeze=True, names=[\"label\", \"sequence\"])\n",
    "# Only keep barcode1 sequences since barcode info isn't needed\n",
    "bc_tag = \"_barcode1\"\n",
    "white_seqs = white_seqs[white_seqs.index.str.contains(bc_tag)]\n",
    "# Trim off the barcode ID\n",
    "white_seqs = white_seqs.rename(lambda x: x[1:-len(bc_tag)])\n",
    "# Only keep the 84 bp of the sequence that corresponds to the library\n",
    "seq_len = 84\n",
    "seq_start = len(\"TAGCGTCTGTCCGTGAATTC\") + 1\n",
    "white_seqs = white_seqs.str[seq_start:seq_start+seq_len]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Read in the data (Dataset S2) and compute predicted occupancies"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Function to correct off by one error in labeling\n",
    "def correct_label(name):\n",
    "    chrom, pos, group = name.split(\"_\")\n",
    "    pos = int(pos) + 1\n",
    "    return \"_\".join([chrom, str(pos), group])\n",
    "\n",
    "white_activity_df = pd.read_csv(os.path.join(white_data_dir, \"white2013Activity.txt\"), sep=\"\\t\", index_col=0, usecols=[0, 1, 2, 3], names=[\"label\", \"class\", \"expression\", \"expression_SEM\"], header=0)\n",
    "# Correct the off by one error of the labels\n",
    "white_activity_df = white_activity_df.rename(correct_label)\n",
    "white_activity_df[\"expression_log2\"] = np.log2(white_activity_df[\"expression\"])\n",
    "\n",
    "white_measured_seqs = white_seqs[white_activity_df.index]\n",
    "white_occupancy_df = predicted_occupancy.all_seq_total_occupancy(white_measured_seqs, ewms, mu, convert_ewm=False)\n",
    "white_occupancy_df = white_occupancy_df.rename(columns=lambda x: x.split(\"_\")[0])\n",
    "white_entropy_df = white_occupancy_df.apply(predicted_occupancy.boltzmann_entropy, axis=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Define strong enhancers and inactive sequences based on scrambled sequences\n",
    "Only do this for bound sequences. Strong enhancers are above the 95th percentile of scrambled, inactive sequences are the 15% centered around the mean of scrambled, silencers are anything that is at least 2 fold below the mean of scrambled."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define cutoffs\n",
    "scrambled_mask = white_activity_df[\"class\"].str.contains(\"SCR\")\n",
    "strong_cutoff = white_activity_df.loc[scrambled_mask, \"expression_log2\"].quantile(0.95)\n",
    "white_scrambled_mean = white_activity_df.loc[scrambled_mask, \"expression_log2\"].mean()\n",
    "\n",
    "# Inactive (not used)\n",
    "# inactive_middle = stats.percentileofscore(bound_activity_df[\"expression_log2\"], white_scrambled_mean) / 100\n",
    "# inactive_percentile_width = 0.15\n",
    "\n",
    "# Pull out bound sequences\n",
    "bound_mask = white_activity_df[\"class\"].str.match(\"CBR(M|NO)$\")\n",
    "bound_activity_df = white_activity_df[bound_mask].copy()\n",
    "bound_occupancy_df = white_occupancy_df[bound_mask]\n",
    "bound_entropy_df = white_entropy_df[bound_mask]\n",
    "\n",
    "# Pull out relevant sequences\n",
    "white_strong_mask = bound_activity_df[\"expression_log2\"] > strong_cutoff\n",
    "white_silencer_mask = bound_activity_df[\"expression_log2\"] < (white_scrambled_mean - 1)\n",
    "white_modeling_mask = white_strong_mask | white_silencer_mask\n",
    "white_labels = white_strong_mask[white_modeling_mask]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Make predictions with the SVM and the occupancy logistic model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Write sequences to file for the SVM\n",
    "white_modeling_seqs = white_seqs[bound_activity_df.index][white_modeling_mask]\n",
    "white_modeling_fasta = os.path.join(svm_dir, \"white2013TestSet.fasta\")\n",
    "fasta_seq_parse_manip.write_fasta(white_modeling_seqs, white_modeling_fasta)\n",
    "\n",
    "# SVM\n",
    "svm_white_tpr, svm_white_prec, svm_white_scores, svm_white_f1 = gkmsvm.predict_and_eval(white_modeling_fasta, white_labels, svm_prefix, word_len, word_len, max_mis, xaxis)\n",
    "write_folds(svm_white_tpr, \"whiteSvmTpr.txt\")\n",
    "write_folds(svm_white_prec, \"whiteSvmPrec.txt\")\n",
    "\n",
    "# Logistic model\n",
    "occupancy_probs = occ_clf.predict_proba(bound_occupancy_df[white_modeling_mask])\n",
    "occupancy_white_tpr, occupancy_white_prec, occupancy_white_f1 = modeling.calc_tpr_precision_fbeta(white_labels, occupancy_probs[:, 1], xaxis, positive_cutoff=0.5)\n",
    "write_folds(occupancy_white_tpr, \"whiteOccTpr.txt\")\n",
    "write_folds(occupancy_white_prec, \"whiteOccPrec.txt\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Assess TF occupancy feature robustness\n",
    "Do this by fitting strong enhancer vs. silencer logistic regression models using the predicted occupancies of random TFs from the HOCOMOCO database. Do this 100 times to generate a background distribution.\n",
    "\n",
    "### Read in database and convert to EWMs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "hocomoco = predicted_occupancy.read_pwm_files(os.path.join(\"Data\", \"Downloaded\", \"Pwm\", \"photoreceptorMotifsAndHOCOMOCOv11_full_MOUSE.meme\"))\n",
    "hocomoco = hocomoco.apply(predicted_occupancy.ewm_from_letter_prob).apply(predicted_occupancy.ewm_to_dict)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Fit model with random PWMs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iteration 10\n",
      "Iteration 20\n",
      "Iteration 30\n",
      "Iteration 40\n",
      "Iteration 50\n",
      "Iteration 60\n",
      "Iteration 70\n",
      "Iteration 80\n",
      "Iteration 90\n",
      "Iteration 100\n"
     ]
    }
   ],
   "source": [
    "wt_seqs = all_seqs[all_seqs.index.str.contains(\"WT\")].copy()\n",
    "wt_seqs = sequence_annotation_processing.remove_mutations_from_seq_id(wt_seqs)\n",
    "wt_seqs = wt_seqs[activity_df.index]\n",
    "modeling_seqs = wt_seqs[silencer_modeling_mask]\n",
    "\n",
    "niter = 100\n",
    "nfeatures = len(ewms)\n",
    "# Track the cross-validated mean TPR and precision for each feature set\n",
    "random_tprs = []\n",
    "random_precs = []\n",
    "# Keep track of the features selected for each round\n",
    "random_ewms = []\n",
    "\n",
    "np.random.seed(seed)\n",
    "for i in range(niter):\n",
    "    if i % 10 == 9:\n",
    "        print(f\"Iteration {i+1}\")\n",
    "        \n",
    "    # Randomly sample PWMs\n",
    "    sample = hocomoco.sample(nfeatures)\n",
    "    random_ewms.append(sample.index.str.split(\"_\").str[0].values)\n",
    "    # Do predicted occupancy scan\n",
    "    features = predicted_occupancy.all_seq_total_occupancy(modeling_seqs, sample, mu, convert_ewm=False)\n",
    "    # Fit the model\n",
    "    clf = LogisticRegression(C=c_opt)\n",
    "    clf, tpr, prec, f1 = modeling.train_estimate_variance(clf, cv, features, labels_with_silencer, xaxis, positive_cutoff=0)\n",
    "    \n",
    "    # Store the result\n",
    "    random_tprs.append(np.mean(tpr, axis=0))\n",
    "    random_precs.append(np.mean(prec, axis=0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Save results to file\n",
    "write_folds(random_tprs, \"randomFeaturesMeanTprs.txt\")\n",
    "write_folds(random_precs, \"randomFeaturesMeanPrecs.txt\")\n",
    "write_folds(random_ewms, \"randomFeaturesList.txt\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}