6_affine_gaps_local.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Local Alignments with Affine Gap Penalties\n",
    "\n",
    "Just like we taxed opening gaps more than extended gaps while doing global alignments, we can do the same for local alignments."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy as sp\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns; sns.set()\n",
    "import pandas as pd\n",
    "import sequence_alignment as seq"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "seq_2 = \"KKASKASPKKPHPHHSPRNGGGGLGKK\"\n",
    "seq_1 = \"KKSGGLLLGK\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Specify an opening ($g_{open}$) and an extension gap penalty ($g_{ext}$)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "g_open = -6; g_ext = -2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Initialization\n",
    "The initialized scoring matrix is identical to the linear gap penalty case: $M_{(i,j)} = 0 \\ \\forall\\ i,j$ "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "#initialize the scoring matrix\n",
    "m = len(seq_1); n = len(seq_2)\n",
    "M = np.zeros([m+1,n+1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The helper matrices $U$ and $L$ need to be initialized as follows: <br>\n",
    "$$U_{(0,0)} = L_{(0,0)} = 0$$ <br>\n",
    "$$U_{(0,j)} = L_{(0,j)} = 2g_{open}\\ \\forall \\ j>0$$ <br>\n",
    "$$U_{(i,0)} = L_{(i,0)} = 2g_{open}\\ \\forall \\ i>0$$\n",
    "\n",
    "where $g_{open}$ is the opening gap penalty."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "#initialize the helpers\n",
    "U = np.copy(M); L = np.copy(M)\n",
    "U[0,0] = 0; L[0,0] = 0\n",
    "U[0,1:] = 2*g_open; L[0,1:] = 2*g_open\n",
    "U[1:,0] = 2*g_open; L[1:,0] = 2*g_open"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Fill Step\n",
    "\n",
    "The recursion relation for local alignments with affine penalties is: <br>\n",
    "$$U_{(i,j)} = \\max \\{ M_{(i-1,j)} + g_{open},\\ U_{(i-1,j)}+g_{ext},\\ 0 \\} $$ <br>\n",
    "$$L_{(i,j)} = \\max \\{ M_{(i,j-1)} + g_{open},\\ L_{(i,j-1)}+g_{ext},\\ 0 \\} $$ <br>\n",
    "$$M_{(i,j)} = \\max \\{ M_{(i-1,j-1)} + S_{(i,j)},\\ L_{(i,j)},\\ U_{(i,j)},\\ 0 \\} $$\n",
    "\n",
    "where $M$ is the scoring matrix, $U$, $L$ are the helper scoring matrices, $S$ is the substitution matrix, $g_{open}$ is the opening gap penalty and $g_{ext}$ is the extension gap penalty. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "S = pd.read_csv('blosum62.csv',index_col=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "tracer = np.zeros([np.shape(M)[0],np.shape(M)[1],7])\n",
    "for i in range(1,m+1):\n",
    "    for j in range(1,n+1):\n",
    "        l_arr = np.array([M[i,j-1]+g_open,L[i,j-1]+g_ext])\n",
    "        L[i,j] = np.max(np.hstack([l_arr,0]))\n",
    "        l_where = (l_arr==np.max(l_arr))*(L[i,j]!=0)\n",
    "                        \n",
    "        u_arr = np.array([M[i-1,j]+g_open,U[i-1,j]+g_ext])\n",
    "        U[i,j] = np.max(np.hstack([u_arr,0]))\n",
    "        u_where = (u_arr==np.max(u_arr))*(U[i,j]!=0)\n",
    "                        \n",
    "        m_arr = np.array([M[i-1,j-1]+(S[seq_1[i-1]][seq_2[j-1]]),U[i,j],L[i,j]])\n",
    "        M[i,j] = np.max(np.hstack([m_arr,0]))\n",
    "        m_where = (m_arr==np.max(m_arr))*(M[i,j]!=0)\n",
    "            \n",
    "        idx = np.hstack([m_where,u_where,l_where])\n",
    "        tracer[i,j,idx] = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXEAAACnCAYAAADnl2wiAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO2deVxN6R/HP7eVhBJRMcg2RppGxiiDMiQUZReRZDSEYcYUsoyfJTGIscfQiGiQaazFZM06kT2iQYsWS6VFdZ/fH829LXfpPLcO5/K8Xy+vl3vu557zPc99nqd7zz3nfUSEEAIGg8FgqCUa77sABoPBYKgOm8QZDAZDjWGTOIPBYKgxbBJnMBgMNYZN4gwGg6HGsEmcwWAw1Bitd73BPhrDOOWixOFqmxdSLSzP8iz/YeQV8c4n8fJ06d8JE5a6QVtXG4/j/8UvXhuRl5PP8gLMC6kWlmf5jz1fnvd2OKV+w3r4cftkLBq6Ep7tpyP18XNMCBjN8gLMC6kWlmf5jz1fmSon8cTERGzYsAHz58/HwoULsWHDBty8eZPzBhRh7WCJhCuJSH6YBgCI3HgC37h1Z3kB5oVUC8uz/Meer4zSSTw0NBQzZ84EAHTs2BEdOnQAAMybNw/bt2/nvBF5NGrWEBnPMqWPM55loU59PejVrc3yAssLqRaWZ/mPPV8ZpcfEQ0JCEBERgdq1K65s/PjxcHV1haenJ6eNyENDQyR3ubhEzPICywupFpZn+Y89L/N6ZU9qaWmhuLhYZnlBQQG0tbU5bUAR6U8y0aCJofRxQ7MGyH6Ri4K8QpYXWF5ItbA8y3/s+cooncS9vb3h4uICf39/BAUFISgoCP7+/hg2bBi8vb05bUAR107cQPuubWDWugkAwMnbAbGHrrC8APNCqoXlWf5jz1dG6eEUZ2dndOnSBbGxsUhPTwchBJ07d8bUqVPRuHFjzhuRx6uMbKz03IB54T9AW0cLKYnPETjuV5YXYF5ItbA8y3/s+cqI3rVPXGgn0LOLfVie5VleHfKKeOeTOIPBYDBqDuZOYTAYDDVG8O6Uns4rOOVPR85SKc/n4RQh1FKdvDq2/cec72f+I6f80Ucr30k9H1ue7/GiCLVyp0z2tIPd1+2QnVMAAHia/AI/B0bWWJ5v34HQ6qHJq3vbf2x5+0GdMPRbOxBCUJhfhE2LIvDg5jO1qV/d83yPl/KojTsFACzam+HnwEh4Td8Jr+k7q9xJmvy78B0IqR7avDq3/ceWN2vZCF6zneDvsRU+TqsRtj4a/hvHqU396p4H+B0vlXlvkzitL0BbSxOtzY0x0vVLbFs7DotmD4Jxo7o1lufbdyC0emjy6t72H1u+6G0x1swOx8uMHABAws1nMGxYF1rammpRv7rn+R4vlXlvkzitL8DISB9x8U+wJeQMJkzbiTv3UrB0rqvC9dPm+fYdCK0emry6t/3Hlk9Pfokrf9+VPv52rjMunbyD4qIStahf3fN8j5fKKJ3EU1JSlP6rDrS+gLTnr+H78348TX4JAAg7eAWmJgZo0rh+jeT59h0IrR6avLq3/ceWl6BbWwdzfnWHafOGWOO3T2FOaPWre57v8SJTn7InJ02ahL59+8Ld3R1jxoyp8M/d3Z3TBhRB6wswb9EIDvafVVgmggglxfI/XdDm+fYdCK0emry6t/3HlgeARqYGWPWHD8QlYvi6bcSb/34wU4f61T3P93ipjNJJfM+ePWjZsiUCAwNx6tSpCv9OnjzJaQOKoPUFEDHBtG+/kf51culvhcSkDGRk5dZInm/fgdDqocmre9t/bHn9+rURuGcyzh+/iYDpoXhbKCuxE3L96p7ne7xURukphvr6+li8eDHCw8NhbW3NaYVcofUFPH6SiaDNJ7Fsnis0NTSQkZmDRSsV/4JLm+fbdyC0emjy6t72H1veabQtGpkawNbBArYOFtLls8dsRs6rPMHXr+55vsdLZQTvTlHHC07YxT7c8kLbX3XPs4t93m/+fV3sw9wpDAaDocYwdwqDwWCoMR+cO6Xt0tWc8glzZgAAevdcyikffXoOAG5fWSVfV2lrp61FaF8PhdD2wMd7uIC2fVh71mye7/GiiA/SnRLg5ICEjCxsv3QNGiIRZvfuie7mzaEp0sC2S9cQFhevcBu2X7eF7xxnDOr/S5X10/opVPEj0NQjFHdKVe1fE/vL3CA1605h7fl+3CnVmaskfHDulJ1uQ9CvfVvp45FfdEQLQwMM2BKCITt2w6PLF7A0kX9XIjMzQ0z6rhc0RPJP7q+QpfRTcK1f1XqE4E5pZdSAU/tXd3+ZG6Rm3SmsPd+PO6U6c1V5qpzEo6Oj8fvvv+PJkycVlu/du7fKlSuDLx/Bgfg7OHo3Qfq4T7vW2B9/GyWEILugEIfv3MdAi/Yyr9PV1YKf/0BsWs/t/HdaPwWtH4G2HiG4U0Zbf86p/au7v8wNUrPuFNae78edoupcVRmlk/jKlSuxa9cuJCUlYeTIkTh06JD0ubCwsCpXrgy+fASHbt2t8NikXl2kZudIH6dl56JJPX2Z133/Yz8cjozDo0fpnOqn9VPQ+hFo6xGCO2XRib85tb88aPaXuUFq1p3C2vP9uFNUnasqo3QSP336NIKDgzFv3jzs3r0bQUFBOHr0KACgumcm8u0jkCCC7HZKxBVrH+jSCSUlYhw7UvXxp8pw9VPQ1K9KPUJyp5RHXvtXRtX2Z26QmnGn0OaFtr9Cy/M5V8mtT9mThBCI/js+2aJFC2zevBlLlizBpUuXpMtVhW8fgYTU7BwY69eRPm5cVx9pORU/DTo4WqJdO1NsCp6ApctHQEdXC5uCJ8DISPlfQRo/BU39qtQjJHdKeeS1f2VU2V/mBqk5dwptXmj7K7Q8n3OVPJRO4o6OjnB3d0d8fOknpDZt2iAoKAjff/+9zDFyWvj2EUg4+SARQz63gKZIhLq6uhjwWVtEJzyskPHx3oGJ47fC22sb5vjuxdvCYnh7bUOWknXT+ilo6lelHiG5U8ojr/2ru7/MDVKz7hTWnu/XnSKBy1wlD6WnGPr4+MDa2hp16pT9dbC2tsaBAwewffv2KleuDL59BBJ2X7uBZgb18aeXO7Q1NRAWdxNXniRXq3ZAuZ+iJuvnipDcKeWR1/5+3/RQaR8lMDdIzbpTWHu+X3eKBFXnqg/OnSKEC07YxT7KYRf78JNnF/u83/z7utiHuVMYDAZDjWHuFAaDwVBjBO9O4fsrjfnaqi9nB4BH034AAFjOqPorUPzq0q8/jlbzOa372PVFAACbUdxqid3DvZby9fD99Zm27T+by63+O0tUa0/HBhO55V9sBcB/3xTa4TLa9mSHX2o2X1MqWrVyp/DlLwCAwN59kZCVheC4qwCAujq62Dt0BHyjj+Nm+nO5r+nVsRUmO9pATAiy8wqxcG8UnmW9VliP84gucBreBYQQpD59gTWL/sTrl28U5ls1a4iZHr2gr6eLErEYy4OjcP+x4othaOrh25VB64kZ3dUKbl0/R2FxMRLTX2Bx5Cm8zld8Ch1tW0r44dfxSLqXjP2/nlCa47tvSuDqiuG7Htr2ZK6V9zdXVUZt3Cl8+QtaGTbALtdhGNCmnXSZXfOWiBjhBnPDBgrXrautiWWj+2HGb5EYvjIUMbcT4TfYTmG+dXsTDB3XDTPGbYX30PVIfvIC46b0Urx+HS2smT0EuyKvYNzs3/HbgYv4ecqAGqmHb1cGQOeJ6dKyKSb06AzP7fsx+NdQnElIwkKXPgrztG0JAM3aNkFAxA/o7lL1HareRd8EuLti+K6Htj2Za+X9zFWKqHIST0pKwvPnpZ9Ew8PDsXjxYhw5coTzBhQhFH+Bu6UV/rh7C4cflDk9xll9gR+jjiH9jeLzOjVEGoAI0K+lCwDQ09HBWwWXKQPAw7up8BwYhLzcQmjraKGhcV1kv1b8l/kryxZIfv4asdcfAwDOXkvE3LWK31iaevh2ZdB6YjqYNUZs4hM8/++S/OjbD2D/aUtoa8rvnrRtCQDOE+wRtfs8zkYotygC/PdNgM4Vw3c9tO3JXCvvZ65ShNLDKTt27MDvv/8OsViMrl27IjU1FX369MH+/fvx+PFjTJkyhfOGKqPMRyDvawdtvry/4GnyS4x0/RJL57rC6/uQCrmFp08BAGybfiJdNv7QgSrrz39bhMXhJ/H79BF49aYAmhoijF2rXApWUiyGjf2n+H7+IBQVlSBko+Jz7ZuZGCLr9RvM+dYBrZs3Qu6bQvy6+0yN1JOe/BLp/10SDFTPlVGdtpcQ/ywNY2ysYGpQFymvcuBq3QE6Wlqor1cbmTnyv9LTtCUAbPDdAwCw6lG1UIjvvgnQuWLeRT007cl3//nY8rTjpTJKP4nv378fR44cwa5du3Ds2DFs3rwZo0ePxsaNG3H8+HFOG1C4YYH6C7jSxsQIkxy6wiUgBL0XbsXWqMtY5eFc5eti/76HEfbLsWvT31iyYaxCfYGWpgZsrVoi4mQ8POeGIvx4HFb5Doa2lvxPO6rUw5crg7btryUlY/2pi1g72hn7JrtBTAhe5eWjqIrLlLm2JS18901aV8y7cq3QtidzrQhjrlI6iYvFYujo6MDMzAyenp7Q1dWVPldSUrU3QxlC9RdwxbZdC1xPSpH+cBh27gZamxjBoE4tuXmTZg3Qwars0/6JiH9gbGIA/Xry85kvc/FvygvcSSz9Snb2WiI0NEQwNZb/xtLWw6crg7bt9XS0cfXxMwxdvxvDN+xG1K3SS41f58uvibYtaeG7b9K6YviuR5X2ZK4V4cxVSidxBwcHjBkzBiUlJZg6dSoA4N69e3Bzc0O/fv04bUARQvUXcOXus3RYtzJDA309AKVnhiRnZePVG/mduUHDuvBbPgz1DErz9v0t8e/DdOQoOPYYe/0xTBrVQ7uWxgAAq0/NQAiQmiH/bBOaevh2ZdC2vXE9fezwGoY6ujoAAO9eX+FwvHzvOEDflrTw3TdpXTF810Pbnsy1Iqy5Sukx8enTp+PKlSvQ1Cz7Cq+jo4OpU6eiZ8+enDagCKH6C7hy+eFT7Dh1Ddt9hqKoWIzXeQWYvv1Phfnbcf8iLPgMAoPHo6REjKyMHPw8Y7fC/IvXefD95RBmefZGLV1tFBWVYPbqPxX+WElTD9+uDNq2T8p8ieAzVxDmPRIaIhH++TcFiyNPKczTtiUtfPdNodVD257MtSKsuarK88S//PLLCo/Nzc1hbm7OeQPKuHw0DpePxvGWj4q5g6iYO5yyP0XLHuPvsSNY6Wv2nr+BvedvcK7ncPgVHA5X/Be5MtfvJcNrHvfJiWs9YRtOImwDtzsGSeCz7QFg98Ub2H2Rv7aU8IvPb5xyfPdNCc/TXsO538r3Xg9Ne76L/vOx5WnHS3mYO4XBYDDUGOZOYTAYDDXmg3On0Pod+HCtSDwrtC4Uvl0rtHm+2552f60ncnOtXNv6n2tFYO4UvtuT73poXTe0bh+h9WehjRdFfNDuFFq/Q025VhRB60Lh07VCWwvfbU+zrwNs2mN0n07Sx/q1ddHYQB/9fgrGixzZH9bKIxR3Ct/tyXc9tK4bWs8Qbf/kuz8LbbyU54N1p9D6HQCeXSuULhQ+XSu0tfDd9rT7ejj2LtwWhcJtUSjGLtmDrNdvsHzP30oncCG5U/huT77roXXd0HqGqL1BPPdnoY2XylBN4gEBATRxpfDtI6D1O/DtWqF1ofDpWqGthe+2V8WFImGcY2e8yMnHgTM3leaE5E7huz35rofWdUPrGaLtn3z3Z6GNl8ooPJwye/ZsmWWnTp3C69elX4GWLVvGaQOK4NtHQOt34Nu1QutCAfhzrdDWwnfb0+6rBAP9WhjTxxqjF4dWmRWSO4Xv9uS7HlrXDa1niLZ/8t2fhTheyqPwk7iBgQFiYmLw6aefokuXLujSpQv09PSk/68u78oHwdXvwLdrhdaFIoEP1wptLXy3vQRad4drj444fSMRKZnZVa6bBqG5NSQIxVVC67qh9frQ9k+++7NQx4t0e4qe8PX1xapVq3DkyBGYmprC1dUV9evXh6urK1xdXTlvQBF8+wgAOr8D364VWhcKn64V2lr4bntVXSgOndvhz/O3lWZUQWhuDUBYrhJa1w2t14e2f/Ldn4U2Xiqj9Ji4jY0NNm/ejN27d2P58uXVll6Vh28fAa3fgW/XCq0LhU/XCm0tfLe9Ki6Uunq6aGZsgPjEVKXrVgWhuTWE5iqhdd3QeoZo+yff/Vlo46UyVZ5iaGBggKCgIISHh+P+fcVvFC18+who/Q58u1ZoXSh8ulZoa+G77VVxoTQzNkDm6zcoruIrqioIza0hNFcJreuG1jNE2z/57s9CGy+V4Xye+LBhwzBsGLeT2bnCp49AFb8D364VWhcKn64V2lr4bnvafb2T9Bwuc7l5UMojFHcK3+3Jd/20rhtazxBt/+SzP/O9flXe3/IwdwqDwWCoMcydwmAwGGrMB+dOoc3T+ju4uFYknhVa1wTfrpXePZdyykefnqNSnrrtKd0mtHlaL47Q+qbQ8rT98+vonzjlz/UOBAB0OLSQU/72oNJc82Bu7++/XqXvL99uFr7HiyI+aHcKbZ7WVULjWgHofRN8u1YAwPbrtvCd44xB/bl1WK552raXwNVtQpunfa+E1jeFlufaN+d8NhyPc9Ow50npxTGRPeYjs7DsvP7d/56ukF/yhQseZKdjR+KFCsvXfDkCGQU5WHLzSIXlK3v0Q8LLTGy5eQW6mlpYbNsblo2aQEMkwvX0VPhfiK5W/armAbqxpep4AT5gdwptntbfAXBzrUig9k3w7FoBADMzQ0z6rhc0ON5gmGuetu0BOreJKnma90pofVNoea59M6jTt+jV2FL6uJleI+QW5WP8pTXSf1FpZT/+bbcdh76mHWTW49m6G6yNPpFZvqf/cDiZlzmMplp1haaGBhwP7EDfAzugq6WFKZ9/pXL9quYBurGlyngpj9JJPD6+7G7csbGxCAgIwMqVK3HjBvdfmRUhND8Frb+Dq2tFAq1vgm/Xiq6uFvz8B2LTem6/itPkadseoHOb0OZp3yuh9U2h5bn2zSMpV3Hqedkc0rF+c5RAjLWdJmHHVzPg0bI3NFA2yR18EofjKRUv3urSsAW+Nm6NfUlXZda/L+EW/npUdtrzpbRnWBcXCwJATAhuZ6XDTL+eyvWrmqcdW6qMl/IoncQXLFgAAAgNDcXSpUvRpEkTNGzYEPPnz8euXbs4b0QeyvwC7yMPlPk7dh3/ARbWLXDi0D8Ks+VdKxOm7cSdeylYOlfxlazxz9LQ1bwZTA1KJ4/yvgl5lPc1bF8yGmvnDIWmhvIvTjT1f/9jPxyOjMOjR8q/EqqSV6XtN/juwcl9FznVQpunfa+E1jeFlufaN4+nVex/miINXMl6gB/iguFzbSO6GLXFkGbdpM9HPouvkG9Uqy78LPrhp2v7USLnJLqDDyueDnw2OQmPs0u1GWb69TChgzUOP5a9toV2bNHmaceWKuOlPJwOp+zbtw8hISHw8PCAh4cHQkNDqz2JC9VPwdXfQetaofVN8OlaGejSCSUlYhw7Ei9nDbLQ5lVte76gfa+E1jeFlle1b0amXEZQwp8oIiXILS7A3idn0MPYQm5WS6SBldZDsfzWMWQW0l01bWHUGOFOo7Dzzj849fRRteunydOOFaD640XpJF5cXAyxWAwjIyPo6elJl+vo6ECjik+FVSE0PwWtv4PWtULrm+DTteLgaIl27UyxKXgCli4fAR1dLWwKngAjI32566bNq+KO4BPa90pofVNoedq+KaFvk05opd9E+lgEEUrE8t+DDgamMKtjgJ8s+mK/nTdGtOgMR7MO+NlqoNJtOJt/itB+w7D8yhmsv3GpRuqnydOOFaD640XpTGxoaIiePXvi4cOH0kMrsbGxGDlyJBwdHTltQBFC81PQ+jtoXSu0vgk+XSs+3jswcfxWeHttwxzfvXhbWAxvr23IUlA7bZ627fmG9r0SWt8UWp62b0ow12+MCeYO0IAIOhpaGNLMFiefy/997cbLZ+h9YjWGxGzCkJhN2Jt0FceSb2PBdcWX6/dv0RYLbXphzLFwHEq8W2P10+RpxwpQ/fGi9BTDkJBSl/ajR4+QnV16WpCOjg6mTZsGOzs7zhuRh9D8FLT+DlrXCq1vgm/XCp/Qtj3f0L5XQuubQsvT9k0J2x9FY+anLtjZdSa0NDTx9/N4RKZchu9nQ5W+jis/fdkDIoiwvHvZB8xrz5OrXb+q+8uV6o4XTueJm5ubS/9vbc3tlC4uCMlPAdD7O2hcKwC9b4Jv1woAPE97Ded+ym/Eqkqetu0lcHWb0OZp3yuh9U2h5bn2zaV3yrzYheIiLLuj/MKVuXERcpdvuB8jd/mPZ45K/28XLt9dNPazL2SW8e1mAejGlqrjBWDuFAaDwVBrmDuFwWAw1BjmTqH1d3DwlUhcJbRuE8sZ3Fwr8atLXSt8uxr6mf/IKX/00UqV8kLrC+pev9DGCq07yHwtN/XDo2mlrpK2S7mtP2FO6fqHXJjMKb/fdoNK66cdv7RuFkUwd4ocuPo4aF0ltPleHVthsqMNxIQgO68QC/dGSW9xpQi+fA32gzph6Ld2IISgML8ImxZF4MHNZwrXTZsXWl9Q9/qFNlZovEGBvfsiISsLwXGlV2nW1dHF3qEj4Bt9HDfTn1fIBjg5ICEjC9svXYOGSITZvXuiu3lzaIo0sO3SNYTFlZ2vTQhB8ra7qGVWBw37NQcAZJ16hpenU0CKSlCreT2YebZXef3loR27qrhZJDB3SjlofBy0rhLavK62JpaN7ocZv0Vi+MpQxNxOhN9gO6U18eVrMGvZCF6zneDvsRU+TqsRtj4a/hvHKa6DMi+0vqDu9QttrHD1BrUybIBdrsMwoE2ZD8WueUtEjHCDuWEDmfxOtyHo176t9PHILzqihaEBBmwJwZAdu+HR5QtYmjQGACQmJiIpMA6vr5T9EXh9NR0vop+ixawv0HpxV5CiEmSdeKLS+stDO3ZVcbOUp8pJ/OzZs9LTCyMiIrBo0SLs37+f8wYUITQfBEDn46B1ldDmNUQagAjQr6ULANDT0VF6ShOfvoait8VYMzscLzNyAAAJN5/BsGFdaGnLv8KNNi+0vqDu9QttrHD1BrlbWuGPu7dw+EHZ9RPjrL7Aj1HHkP5G9jzrA/F3cPRugvRxn3atsT/+NkoIQXZBIQ7fuY+BFqWfrENDQ2HQ3QT1vyybdF9dSIVR30+gpa8NkYYIpmM/hYFtE5XWXx7asUvrZqmM0sMpS5Yswd27d7F69WqsWbMGN2/exDfffIOoqCjcvXsX/v7+nDdUGWW+AHlf4/jOA6U+DgCw6iH7xshD4ir5fv4gFBWVIGTj9hrL578twuLwk/h9+gi8elMATQ0Rxq7dqzDPp68hPfkl0v+7ZB0Avp3rjEsn76BYQcekzQutL6h7/UIbK/HP0jDGxgqmBnWR8iqngjcoM6fscOLC06XXTdg2LbvyePyhAwrXe+jWXdi0aCZ9bFKvLlKzc6SP07Jz0c64IQBg/vz5uHFhMt7cKXtf3z7PQ3FOEZJ+iUPRq7eo07Y+mgxvo9L6y0M7dsu7WVo3b4TcN4X4dfcZhfnKKJ3Ez58/j8jISGhqaiImJgb79u2Djo4ORowYAScnJ84bkYfQfBCqEvv3PcT+fQ+Og62xZMNYeDoH1Ui+jYkRJjl0hUtACJ5lvYZbdyus8nDGsJWyzpryvobGTZRf+ixBlfbRra2DH1aMQCMTA/h7bK1yG1zzQu0L6lq/0MZKeW+QmAAHrt1S6g1SFRFk96NErPgMalJM8OZ2Fj6Z9jlE2hpIDr6D5/sTAfvqrZ9m7AJlbpYp/9uHO4lp6G7dCqt8B8N16lZObaT0cEqtWrWQlZUFADAyMkJeXumdl/Pz86GlVb3fRIXmg6CF1rVCm7dt1wLXk1KkP4aEnbuB1iZGMKhTfRcKQN8+jUwNsOoPH4hLxPB124g3OfKdL6rkhdgX1Ll+oY0VWm+QqqRm58BYv470ceO6+kjLUXy5u5ahLup1MoZmbS1oaGnAwKYJ8hIV//jIdf00YxdQ3UUjQekk7uPjg6FDh2L58uUwNzeHu7s7li5diuHDh2P8+PGcNqAIofkgaKF1rdDm7z5Lh3UrMzTQL8336tgKyVnZePVGtuPz7WvQr18bgXsm4/zxmwiYHoq3hcUK16tKXmh9Qd3rF9pYofUGqcrJB4kY8rkFNEUi1NXVxYDP2iI64aHCfP3Oxnh95TnEb0tACEH2Pxmo3ULWP067fpqxC6juopGg9ON0r1690KZNG0RHR+Pff/+FlZUV6tSpg4CAAFhaWip7aZUIzQdBC62rhDZ/+eFT7Dh1Ddt9hqKoWIzXeQWYvl2x/IcWmvZxGm2LRqYGsHWwgK1DmTp09pjNyHmVV+280PqCutcvtLFC6w1Sld3XbqCZQX386eUObU0NhMXdxJUnsu4UCQ16NUVJbhESf74MIgZqN68L05FtFOa5rp927FbXzVLlMZFmzZpV+1O3IoTmg5DA1cdB6yqhze89fwN7z9PdRYkPX0PYhpMI28DtrBdV8jS1vIu8utf/LvISuI4VGm/QT9HHZZb12CHfi+L3V9n56SWEYGn0abk5CU29ypTEIg0RjF3MYexirjBPu34JtGNXFTeLBOZOYTAYDDWGuVMYDAZDjfno3Sl8+COknhVK1wStu4PWvUBbj9DanuXVK289kZtL5NpW1VwitOvvcGghp/ztQaW55sErOOX/9ZoFgLlTBOGDUMUfwdUdQZundXcAqvkXuNYjtLZnefXJD7Bpj9F9Okkf69fWRWMDffT7KRgvcmR/KAbo+jLN+pd84YIH2enYkXihwvI1X45ARkEOltw8UmH5yh79kPAyE1tuXoGuphYW2/aGZaMm0BCJcD09Ff4XouXWxNwpapCncUfQ5mndHQC9f4GmHqG1PcurV/5w7F24LQqF26JQjF2yB1mv32D5nr8VTuC0fZnr+rfbjkNf0w4yr/ds3Q3WRp/ILN/TfziczMs8LlOtukJTQwOOB3ag74Ed0NXSwpTPv5KtX0julMWLF+P1a27nKtIiNB8EbZ7GHUGbp/BgJ8cAACAASURBVHV3APT+BZp6hNb2LK9e+fKMc+yMFzn5OHDmpsJMdVwiytZ/8EkcjqfcrrCsS8MW+Nq4NfYlXZXJ70u4hb8elZ3PfintGdbFxYIAEBOC21npMNOXPa/8XbtTlE7iERERGD58OE6c4HaogAYad4cQ8xt89+DkvouKd7Aa+fTkl7jyd9mNXqtydwAV/Qvbl4zG2jlDoamh+O2lqUdobc/y6pWXYKBfC2P6WOOXvTFKc7R9mev6I59V1MY2qlUXfhb98NO1/SiRc5LewYcVb+d3NjkJj7NL3Stm+vUwoYM1Dj+WvWipvDsleuFEjOr+OVb/dVZh3arurwSlyaZNm2L9+vUICQnBsGHDcOTIERQU1MylskLzQbwr1woNurV1MOdXd5g2b4g1fvuUZiX+hYiT8fCcG4rw43FY5TsY2lqKP71zRWhtz/LqlZfg2qMjTt9IREpmttKcqn2Z6/oBQEukgZXWQ7H81jFkFiq+slkeFkaNEe40Cjvv/INTTx/JPF/endJ74VZsjbqMVR7Oimup5thVOomLRCK0bt0au3btwowZM3D8+HF88803GD16NH74QfkvplUhNB8E3/4IWmhdJdX1LyhDaG3P8uqVl+DQuR3+PH9baQZQvS9zXT8AdDAwhVkdA/xk0Rf77bwxokVnOJp1wM9WA5W+ztn8U4T2G4blV85g/Y1LcjOCcqeUvw7I1tYWQUFBiImJwezZs9G7d29OG1CE0HwQfPsjaKB1dwDV9y8oQ2htz/LqlQeAunq6aGZsgPjEVKU5QLW+TLN+ALjx8hl6n1iNITGbMCRmE/YmXcWx5NtYcF3x5fH9W7TFQpteGHMsHIcS7yrMCcqdMnq07C/O2trasLCwgIWFhZxXcEdoPgi+/RE00Lo7gOr7F5QhtLZnefXKA0AzYwNkvn6DYg6HJ1XpyzTrV5WfvuwBEURY3t1Ruuza8/fvTnnnl90L7YIEdrGP4nqE1vYsr155drGPfGr6Yh/mTmEwGAw1hrlTGAwGQ41h7hQe8kKqpXze0Wo+p/yx64veST0s/2HnezpzOxxxOrL0cETvnks55aNPzwFAf/jRfC23wxePppUevmi7lNvhkYQ5M1TKM3cKy1c77zyiC5yGdwEhBKlPX2DNoj/x+uUbuVmh1c7y6pWf7GkHu6/bIfu/U2WfJr/Az4HKr0q0/botfOc4Y1D/qic7GtdQYO++SMjKQnBc6VWadXV0sXfoCPhGH8fN9OcVsgFODkjIyML2S9egIRJhdu+e6G7eHJoiDWy7dA1hcfHVyktg7hSWp863bm+CoeO6Yca4rfAeuh7JT15g3JRealE7y6tXHgAs2pvh58BIeE3fCa/pO6ucwM3MDDHpu17QEMm/sKhClqNrqJVhA+xyHYYBbcp8KHbNWyJihBvMDRvI5He6DUG/9m2lj0d+0REtDA0wYEsIhuzYDY8uX8DSpLHKeQm8ulMAIDY2FnFxpXf82L59O7y9vfHrr7/i7du3nDciD6H5HT62/MO7qfAcGIS83EJo62ihoXFdZCu436fQamd59cpra2mitbkxRrp+iW1rx2HR7EEwblRXYV5XVwt+/gOxaT23uytxdQ25W1rhj7u3cPhB2aXy46y+wI9Rx5D+RvaqzQPxd3D0boL0cZ92rbE//jZKCEF2QSEO37mPgRbtVc5L4NWdEhgYiNWrV+N///sfpk2bhhs3bmDUqFHIyMjAokWLOG9EHkLzO3xseQAoKRbDxv5T7Dr+AyysW+DEoX/UonaWV6+8kZE+4uKfYEvIGUyYthN37qVg6VxXuVkA+P7HfjgcGYdHj7gdTuDqGlp4+hQi7lW8SGf8oQOIS5N/gdChWxWzJvXqIjU7R/o4LTsXTerpq5yXUF13itJj4mfPnsWhQ4fw9u1b9OzZE+fOnYO2tjZ69OiBQYMGcd6IPITmd/jY8hJi/76H2L/vwXGwNZZsGAtP56B3XgvLf9j5tOev4fvzfunjsINXMHakDZo0ro+05xWvShzo0gklJWIcOxKPxk3olBG6tXXww4oRaGRiAH+PrVSv5YIIsvtdIlZ8hjbXvMSdMuV/+3AnMQ3drVthle9guE7diqLiqi/4qfKy+5ycHLx8+RIFBQXIzS39ylFQUICioqIqV64MofkdPra8SbMG6GBV5lA+EfEPjE0MoF9P1u8gtNpZXr3y5i0awcH+swrLRBChRM4E5eBoiXbtTLEpeAKWLh8BHV0tbAqeACMj2U+w5aF1DalCanYOjPXrSB83rquPtBzF8iyueV7dKRMnToSDgwMGDRqEWbNmwdPTE4GBgXBzc8OQIUM4bUARQvM7fGz5Bg3rwm/5MNQzKPU72Pe3xL8P05Ej57i40GpnefXKEzHBtG+/QZPGpZOSS38rJCZlICNLdkLz8d6BieO3wttrG+b47sXbwmJ4e21DlpysBFVcQ6pw8kEihnxuAU2RCHV1dTHgs7aITnhY7Tyv7pRBgwahb9++KCkpQZ06dfDll1/i3Llz+PHHH9GtWzdOG1CE0PwOH1v+dty/CAs+g8Dg8SgpESMrIwc/z9itFrWzvHrlHz/JRNDmk1g2zxWaGhrIyMzBopXcf7irClVcQ6qw+9oNNDOojz+93KGtqYGwuJu48kTWnUKbZ+4UAeaFVEv5PLvYh+XfZZ5d7KM8z9wpDAaDwWDuFAaDwVBnmDvlIzqcwod2F2Dq2o8lT/v1nzbP9+EO2nom/zOGU35Dp10A6Ov/bC63+u8smaH0eeZOYXnOeQD44dfxSLqXjP2/Kr95ttBqZ/mazdO6PmjyfLpNaOohhODupseo06w2mjuZIH71Q+Q/Lzt1MT/9LQzb1wX2qF6/hNFdreDW9XMUFhcjMf0FFkeewut8breGZO4UlueUb9a2CQIifkB3F2uFGaHWzvI1m6d1fXDN8+02oaknMTERcYvv4/nFl9JlljNa46sAC3wVYIH2E1tAu44m2nmWXWtBW7+ELi2bYkKPzvDcvh+Dfw3FmYQkLHTpozBfmSon8ejoaPj5+WHChAnw9vbG0qVLpS6V6iA0vwPLK887T7BH1O7zOBtxTWFGqLWzfM3maV0fXPN8u01o6gkNDYWJXUM07moo83pxsRi3Nz5Gm7GfoJaRrsr1S+hg1hixiU/wPLs0E337Aew/bQltTW6fsZWmNm/ejP3798PS0hIikQhWVlYwNTWFv78/9u3bx2kDihCa34Hllec3+O7ByX0X5T4n9NpZvmbztK4Prnm+3SY09cyfPx8m3RvK3W7K35nQNdSB8ZcVJ3ja+iXEP0tDV/NmMDUolYK5WneAjpYW6usp9h6VR+kkfuTIEWzYsAFubm5Yv349Lly4AA8PD4SFheG3337jtAGFGxaY34HlledpEFrtLF+zeYnrI+JkPDznhiL8eBxW+Q6GtpZmjeRVhdZVomo9T46koaWrSbXrlXAtKRnrT13E2tHO2DfZDWJC8Covn5M3BahiEi8sLER+fumPGwUFBXj16hUAQE9PDxoUli15CM3vwPLK8zQIrXaWr9k8reujum4QrrwLV0nO4zcgYsCgvWKVLi16Otq4+vgZhq7fjeEbdiPqVuml+a/zuflflM7EgwcPxqhRo7BixQqMHTsWgwcPRnJyMoYMGQInJ6dqFS40vwPLK8/TILTaWb5m87Suj+q6QbjyLlwlL+/moEGHuhBxuFkFV4zr6WOH1zDU0dUBAHj3+gqH4+9X8aoylJ5i+O2336Jjx464c+cO/Pz8YGNjgzdv3mD58uVo166dspdWidD8DiyvPE+D0Gpn+ZrN07o+qusG4cq7cJXkpRWiVkPdKnM0JGW+RPCZKwjzHgkNkQj//JuCxZGnOL+euVPYxT4K8+xiH5Yvn2cX+8jnfV/sw9wpDAaDocYwdwqDwWCoMYJ3p9DqLGnzfB5OEUIt1cnT1k/7dVVo+6vueVo1K9/10Kplafsb7eEIyxnc8vGrS/O0h0fEaW2rSJai0aT0AqUhFyZzyu+33aD0ebVyp0z2tIPd1+2Q/d+tl54mv8DPgYqvFKPN8+2bEFo9NHna2mndGkLa1w8hbz+oE4Z+awdCCArzi7BpUQQe3Hz23uqRYPt1W/jOccag/sonSNr+Ruse6dWxFSY72kBMCLLzCrFwbxSeZck/O4WLC+XQoUPYthkQiYBausDcacBnbYFfNgOnLwIaIqB5U+DnH4GGpScBgRCC5G13UcusDhr2aw4AyDr1DC9Pp4AUlaBW83ow85S94rQyauNOAQCL9mb4OTASXtN3wmv6TqVvKm2eb9+E0OqhzdPUTuvWENq+qnverGUjeM12gr/HVvg4rUbY+mj4bxz33uqR1mVmiEnf9YIGh9PzaPobrXtEV1sTy0b3w4zfIjF8ZShibifCb7CdTI6rC+XRo0dYsWIFtqwADm4DvMcC0+YB+48AdxKAA1uBP3cAn5gBy9eXviYxMRFJgXF4faVMiPX6ajpeRD9Fi1lfoPXiriBFJcg68aTKtnpvkzitr0FbSxOtzY0x0vVLbFs7DotmD4JxI8Un3NPm+fZNCK0emjxt7bRuDSHt64eQL3pbjDWzw/Eyo/QS9ISbz2DYsC60tOVfkch3PQCgq6sFP/+B2LT+pNIcQN/faN0jGiINQATo1yo9VVBPR0fu6YVcXSg6OjpYvHgxjI1KH1u0AzJflH7y/vE7QEenbHnKf3N2aGgoDLqboP6XZYKuVxdSYdT3E2jpa0OkIYLp2E9hYNukyvZ6b5M4ra/ByEgfcfFPsCXkDCZM24k791KwdK6rwvXT5vn2TQitHpo8be20bg0h7euHkE9Pfokrf5c5PL6d64xLJ++gWMF50HzXAwDf/9gPhyPj8OiR4kNqEmj7G617JP9tERaHn8Tv00cgeuFEjOr+OVb/dVYmx9WF0rRpU9jZ2QEACCn9tG3fDehiBXT47zD56xxgQwjgWBrD/PnzYWhb8dL9t8/zUJxThKRf4vBg3iWkH3oETT1thfstocpJ/OzZs5g7dy4mTJiAiRMnYu7cuTh+/HiVK65yw5S+hrTnr+H78348TS5VQ4YdvAJTEwPpHbSrm+fbNyG0emjytLXTuimEtK8fQl6Cbm0dzPnVHabNG2KNn2JhHd/1DHTphJISMY4dke/2rgxtf6N1j7QxMcIkh65wCQhB74VbsTXqMlZ5OHOqTRl5+cCMBcC/ycD/ZpUtf5IMuE8FrDsCbor/FoEUE7y5nYVmkzui1YIvUfKmGM/3J1a5XaWTeFBQEHbs2IEuXbrAy8sLnp6e6NKlC/744w8sX76c887Jg9bXYN6iERzsP6uwTAQRShS8UbR5vn0TQquHJk9bO62bQkj7+iHkAaCRqQFW/eEDcYkYvm4b8SZHsYeD73ocHC3Rrp0pNgVPwNLlI6Cjq4VNwRNgZCRrGATo+xute8S2XQtcT0qR/pAZdu4GWpsYwaBOLbl5LqSkpMBtCqChCexcA9T77+jPpX+AUZMBF0dg4Q+lP3wqQstQF/U6GUOzthY0tDRgYNMEeYlVqwCqtBhu3boVgwYNgo2NDWxsbDBo0CBs2rQJMTExNPsoA62vgYgJpn37jfSvsUt/KyQmZSAjS76nlzbPt29CaPXQ5Glrp3VTCGlfP4S8fv3aCNwzGeeP30TA9FC8LSxWmH0X9fh478DE8Vvh7bUNc3z34m1hMby9tiGrhsYurXvk7rN0WLcyQwN9PQClZ6okZ2Xj1RtuwqnKvHr1CmPGjEGfHsCqBaVnpwBA3C1g6jwgYA7gObLq9dTvbIzXV55D/LYEhBBk/5OB2i3qVfk6pacY6urqIi0tDaamphWWp6SkQEdytF5FaH0Nj59kImjzSSyb5wpNDQ1kZOZg0UrFP5bR5vn2TQitHpo8be20bgoh7euHkHcabYtGpgawdbCArYOFdPnsMZuR8yrvnddDC21/o3WPXH74FDtOXcN2n6EoKhbjdV4Bpm//U+V69+zZg9TUVESfBaLLHVo3rF96jHzVltJ/AGDWBNiwTf56GvRqipLcIiT+fBlEDNRuXhemI9tUuX2lk7ifnx9Gjx6NFi1aoFGjRgCAjIwMJCUlYdmyZdz2UAmXj8bh8lHudwmKirmDqJg7vOVp6+G7fr7rocnT1n79XjK85u3mpRaWV07YhpMI21D1WSDvqp7yPE97Ded+K6vM0fa33RdvYPfFG5zze8/fwN7z3PI/Rcv+BthjR7D0/9999x2+++47zhf7lKepV9lhI5GGCMYu5jB2MadaR5XulMLCQsTHxyM9PR2EEDRu3Biff/55tT+JMxgMBqP6KJ3EU1JSlL648mEWBoPBYLxblE7izs7OSEpKgrGxMSrHRCIRTp6k+8rGYDAYjJpF6SSem5sLNzc3LFiwANbW1u+yLgaDwWBwQOkphvr6+li8eDEiIiLeVT0MBoPBoIDdFILBYDDUGHZTCAaDwVBj2CTOYDAYagybxBkMBkONYZM4g8FgqDFsEmcwGAw1hk3iDAaDocZoLly4cOG73GBiYiLCwsIQGRmJ06dP486dO6hVqxYaN25c9Ys5EB0djfPnz8PQ0BD165f5q/fu3QsLCwuZfFJSEgoLC6Gvr4/w8HAcPHgQubm5aNOmansYAAQEBODrr7+W+1x8fLx0v2JjYxEaGoqLFy+iTp06aNJE/m2Xzp49C0NDQ+jq6iIiIgL79u1DZmYmPvvsM5ns4sWL8fnnn6NWLe4e5NjYWGRmZsLExATbt2/H5s2bkZSUBCsrK2hqyt60ITo6Glu2bMGePXtw9OhRxMfHQ19fHyYmJnLWzmAw3jXv9Dzx0NBQ7Nu3D3379q1gRTxx4gQGDhwIT0/Paq1/5cqVuHXrFlq1aoWjR4/C19cXgwYNAgC4urri4MGDFfI7duzA77//DrFYjK5duyI1NRV9+vTBqVOn0KlTJ0yZMqVCfvbs2TLbPHXqFHr16gUAMmZHyTZDQ0MRFhaGIUOGAAAOHjyIYcOGYcyYMRXyS5Yswd27d7F69WqEhobi5s2b+Oabb3DmzBk0bdoU/v7+FfKdO3eGkZERfvjhBzg4OFTZPoGBgbh69SqKi4vRtGlTiEQiDB48GKdOnUJJSQkWL15cIb9582Zcv34d3bt3x6lTp9C5c2fUqlUL4eHhGDduHIYPH17lNhkfDmfPnsWxY8eQlpYGDQ0NGBsbo0ePHujbt2+1111cXIzQ0FCkpqaid+/e6Ny5s/S5devWYerUqTKvuXDhAurWrYv27dtj3bp1uH//PqytreHp6Sn3AwkN0dHR6N27NwAgPDwcZ86cgZaWFvr06YP+/fvLrf+PP/5Anz59ULduXWzZsgU3b95Ehw4dMGnSJOjq6lbIe3t7Y+7cuWjWrFm16gQAkHeIg4MDycvLk1mel5dH+vbtK/c1ycnJSv+Vx8nJiRQVFRFCCHn8+DGxt7cnR44cIYQQMmjQIJl1Ozk5kYKCAvLs2TNiZWVFCgoKCCGEFBYWEmdnZ5l8QEAA6dq1K9mxYwc5cOAAOXDgALGzs5P+vzIuLi6EEEIGDhxIXrx4IV2ek5Mjd3/79etHiouLpfUWFhYSQggpLi4mjo6OMvlBgwaRBw8ekNGjR5OhQ4eSw4cPk/z8fDmtWLa/JSUlJD8/n3Tp0oW8ffuWEEKIWCyWu78DBw4kYrGYEEJIQUEBcXd3J4QQkp2dLbceQgg5c+YMmTNnDvH09CReXl5kzpw55NixYwproqGoqIjs2LGDLFu2jFy5cqXCc2vXrpX7mvPnz5P4+HhSVFREVq1aRSZNmkS2bNkibefqEBUVJf3/vn37iI+PD/n+++/J4cOHFda/Z88ekpmZSQoLC8m6devIt99+S4KCgqR9rzyTJk0iT548qXadNcGaNWuIp6cniYiIIBcuXCAXLlwgERERxMvLiwQEBMjkL1++rPRfZWbPnk18fX3Jtm3biL29Pdm4caP0Ock4Kk9gYCAZOXIkcXJyIp6enmTWrFkkJiaG+Pn5kUWLFnHapxkzZih8TrLNtWvXkrFjx5KoqChy4sQJMmHCBLJq1SqZ/MyZM8nMmTPJixcvyKJFi4ifnx+JiYkhCxcuJDNnzpTJd+3alfTp04ds27ZNOg5VRalPvKbR0tJCcbHsXUYKCgqgrS3/hqCTJk3iLOEihED03/2PWrRogc2bN2P8+PFo0KCBdHl5xGIxdHR0YGZmBk9Pzwp/LUtKZG9g4Ovrix49emDNmjWYOXMmvvrqK+zcuROurvJvnFdcXAyxWAwjIyPo6elJl+vo6EBDzo2Da9WqhaysLBgbG8PIyAh5eXnQ0dFBfn4+tLRk3yqRSITWrVtj165duHDhAvbu3YslS5agRYsWaNKkCX755ZcKeUIIcnJykJeXh4KCAuTm5sLQ0BAFBQUoKiqSWX9hYSHy8/Ohp6eHgoICvHr1CgCgp6cnt/6goCDEx8dj4MCBMDYuvatPeno6/vjjD1y/fh2+vr5y24kr8+fPh1gsRtu2bfHTTz9h+PDh8Pb2BlD6jajyp7UVK1bgn3/+QW5urrRNR40ahWPHjmHp0qWYN29elducOXMmVq1aJfe59evXo3fv3li3bh2uXr0Kd3d3EEKwd+9e3L9/HzNmzKiQl+x/3759sXz5cuTl5cHNzQ0xMTGYM2eOzPt148YNTJgwASNHjoS7u7vCMSKhuLgYERERqFWrFvr27Ytly5bhypUrsLCwgK+vLwwMDCrkxWIxQkJCcPLkSWRkZEBbWxuffPIJ+vfvjwEDBlTIHjlyBEePHpV5352cnODk5CTz3q5fvx7Xr1+HpaWl3HEbEhJSYdmtW7fw55+lN2ZwcXGBh4cHatWqBQ8PD5nXA8Dp06cRGRmJV69eoU+fPrh8+TI0NDTQo0cPuLi4yOTd3d1l5oBbt25h7NixACBTj4SoqCiEh4dL5wY7Ozs4OTnJvLcJCQmIjCy9ccWVK1dw6NAhiEQi9OzZU+4n98aNGyM4OBiBgYFwcHDAqFGjMGDAAJiZmcmtQxnvdBL39vaGi4sLbGxsKhxOuXjxokyjSNizZw9nCZejoyPc3d3h5+cHS0tLtGnTBkFBQfDx8cHbt29l8g4ODhgzZgxCQkKkE8C9e/fg7++Pfv36yd2GjY0N2rdvjwULFiAmJkbuZC/B0NAQPXv2hEgkwoIFCxAQEIDY2FisWLECjo6OMnkfHx8MHToUAwYMgLm5Odzd3WFjY4Nz587By8tLJl++c9va2sLW1hZFRUW4f/8+nj59KpOfOHEiHBwcQAjBrFmz4OnpCRsbG8TGxkoP9ZRn8ODBGDVqFL7++mucO3cOgwcPRnJyMqZMmQInJyeZPO1Ar8rJU3kwsoGufKD7+/sjLy8Pb9++xa5du2BpaYnVq1fj5MmTmD9/PtauXVshHxAQgKKiInh5eeH48eP49NNPYWJigpCQECQlJVU4nEh7l6+tW7di7NixGDduHL755hu59ZaHEIK8vDzo6emhQYMG2Lp1K0aNGgUjIyO5H8AA4O3btzA0NISvr6+0z71580buB0VHR0ds2bIF33//PczMzEAIwbx58+Dj4yN33Xl5ecjMzISpqSny8vKk721BQYHcD1R6enp48OAB2rRpg08++QSpqakwNTXF8+fP5baPSCRCw4YNERgYiKSkJOzbtw+enp4oLCxEkyZNEBYWVmWbSanW53gVSEtLIwcPHiSbN28mmzZtIgcPHiRpaWlKX3Pjxg3i7+/Paf0XLlwgDx8+rLAsJSWFLF68WG6+8le7xMREEhMTw2lb+/btI+PHj68yl5iYSOLi4gghhFy9epX8/fffCrNPnjwh27dvJwsWLCD+/v5kzZo15MaNGwq3T0t+fj7Jzc0lhBBy7949EhwcTM6dO6cwf+HCBRIcHEwuXLhACCEkNzeX3Lt3T27W2dlZ5hAXIaX7NHDgQJnlvr6+xMrKivj5+cn9VxknJyfy5s0b6eO0tDRib29P/vzzT7lfuQcMGCA9TFG+rXJyckj//v1l8rt27SI9evQgBw4cIJcuXSIXL14kffr0IZcuXSKXLl2SyTs4OJCMjAwyadKkCofLsrOz5a5/+PDhJCEhgRBCyJQpU6RtlZaWJvdwX/l9evz4MVm+fDlxcHAgPXv2JCNGjJDbPoSUHn7r1q1bhefktX/5Q2glJSVk5MiRhJDSw4mVD/edP3+e2NnZEQ8PDzJr1iwya9Ys4uHhQezs7EhsbKzMugkh5NGjR2TFihVyn6tMWFgYcXBwkPYzQgh5+PAhsbOzIx07dpTJh4aGEgcHhwqHxa5du0bs7e3J3r175W7jwYMHxM3NjRw8eJAQIv8wjQQ/Pz/i5ORErK2tiY+PDyGEkOPHjxM7Ozuye/dumfw///xDunXrRqZOnUp++ukn0q1bNzJ58mTSs2dPufOJvPebEEJevHihcLwr4p1P4owPF1UG+qRJk0h4eDin9bOBXoqigT5w4EDy6NEjcuPGDWJpaUmePn1KCCEkKytL7h+V/v37k8zMTEJI6R8Syb7m5uaSAQMGyOT/+usvsm7dOrJhwwYSGRlJLl++TAoLC0lYWJjcOqOiokhISAj5999/KyxXlN+yZQv59ddfK+RzcnLI5s2b5eZ/++23CutPTU0lCQkJCtdPSOkfqGXLlpGpU6fKbZPK9e/cuZM8fvyYEELI/fv3yb179xSu/+DBg2TdunVk69atZOPGjSQ8PJykpqbKzcfExFC3jyLYJM6oUQoKCsjly5fJX3/9VWGgK+L58+ckODiY8/ofP34s82k/JyeH/Pbbb3LzkolMgmSgK4NmoBNS+sO85AdIyUBXRE5ODjl+/DjZsmVLhYEuD67fCCWcPXuW9OzZk3Tv3p1ERUURR0dHMm3aNGJvby/9o1Se/fv3kx49epDp06cTOzs7EhkZSZKSkkjPnj3JH3/8USG7YsUK4uHhQRYtWkS6du1KIiIipM/J+0O3YsUKstXGJAAAAchJREFUMm7cOLJo0SJiY2Pz3vOEVPyjcv78eTJ79mxCiPxJs/z+cq2Hpn0CAwOp61cEm8QZNQbNmURCzFfm3Llz0oGuDvtbmYyMDHL06FGZw4vlefToETl69Kj002ZhYSF5+fKlTE6VM79o85JvTHysn3bSfxf7S5NXxjv9YZPxYUNzJpEQ85XvKduyZUv4+PhIl1f+Ua8m6hGJRNKzqmq6fgCwtLSUPifvR0ldXd0KGQl5eXkV8oTyzC9V8hL4WP/p06dx8OBBaGlpwd3dHZ6entDR0UG/fv3k/ij+LvaXJq8UqimfwVBCTk4OcXZ2JlevXlXLvJOTE7GwsCC9evUi9vb2xN7eXvr/Xr16qWX95fejOvl169aRUaNGVTgWf/XqVdK1a1fSqVMnmXULLT9gwIAKv40kJCSQbt26kYsXL8r9JC60+pXBJnFGjUJzJpHQ8rSTJt/10Ob5/iNBe+aXkPKqTJpCql8Z7PZsDEY54uPjER4ejv/973/vuxSVoK1f3feXhtjYWBgbG6NVq1bSZampqdi+fTvmzp37HiurHmwSZzAYDDWGqWgZDAZDjWGTOIPBYKgxbBJnMBgMNYZN4gwGg6HG/B812g/uvF2n5gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sns.heatmap(M,linecolor='white',linewidth=1,cmap=\"viridis\",square=True,annot=True,cbar=False);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's visualize the `tracer` sub-matrices:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAsMAAAILCAYAAAANN5GbAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOzdeXyU5bn/8W/IIpvIUUkoIgZZpEVAiwtRQ0qRBMhmwApBghBwqRsuhYDVQhELUitG9ODyIp6i9EhKITQURJADEgxQ2sa0lR+GQARMIMhhCQLZuH9/+CIHhMxMZssz83zefznzZK65npm5cl/ePLkmxBhjBAAAANhQi+ZOAAAAAGguNMMAAACwLZphAAAA2BbNMAAAAGyLZhgAAAC2RTMMAAAA26IZBmA5/+///T9t3bq1udPwiW3btumGG27Qt99+26x5rFu3ThUVFY0eX7BggUaMGOHHjACgedAMA7CcRx99VLt3727uNILW119/rccff1xVVVWN/kxmZqYWLVrkx6wAoHmENXcCAAD/cuW7ltq0aaM2bdr4IRsAaF7sDAOwlIyMDH399dd68cUXlZGRoQMHDuiGG27Qf/7nf+q2227Tgw8+KEn685//rOTkZN1444368Y9/rJ///Of65ptvGuJs27ZNo0aNUr9+/XT33Xdr2bJlDcf27t2riRMnql+/fho0aJDmzZunmpqahsfdfvvtevnll9W/f3/96le/uijHkydP6tlnn9Xtt9+um266SRMnTlRZWZkkafny5br99tsv+Plp06bpySefvOC+5cuX66677tLNN9+srKwsnTx5stHXZNu2bbrzzjv1l7/8RQMHDtTNN9+sWbNm6auvvtLYsWPVt29f3XvvvdqzZ0/DYxy9PoMHD5YkJScna8GCBVq+fLnS0tL0i1/8Qj/+8Y+1cOHCiy6TcPf1lKQ333xTcXFx6tOnj+655x5t2rSp0XMFAH+jGQZgKQsWLFDHjh31zDPPaMGCBQ33f/rpp8rNzdXUqVP197//Xc8995wmTpyotWvX6s0339TOnTv11ltvSZL27NmjiRMnqn///srLy9OTTz6pmTNnqrCwUNXV1Zo4caKio6O1YsUKzZs3T5s3b9bs2bMbnuvYsWPat2+fVqxYoczMzItyzM7O1p49e/T73/9eK1asUGhoqJ577rkmnWdubq7eeOMNvfvuu/rHP/6hl19+2eHPHzt2TH/84x+1aNEizZo1S0uWLNGECRM0fvx45ebmqqamRq+++qokOX19/vjHP0qS3n///Ybz++KLL9SuXTutWLFCqampFzy3J6/n+vXrtWjRIs2ZM0dr1qxRXFycnnrqKYfNPwD4E5dJALCU9u3bKzQ0VG3atFH79u0bmqZx48YpOjpakvTvf/9bL774ou655x5J0jXXXKPBgwerpKRE0nfNXo8ePTR16lRJUteuXXXixAmdPXtWq1atUnh4uF544QVJ0vXXX69f//rXuv/++xt+XpIeeughdenS5ZI5HjhwQK1atdI111yjyy+/XC+++KK+/vrrJp3nrFmzdNNNN0n6v53j5557Tq1atbrkz9fV1WnKlCnq0aOHevTooTlz5mjw4MG6++67JUkpKSkNTe5ll13m8PW58sorG17r8y+FePTRR3X11Vdf9NyevJ4HDhxQeHi4OnXqpM6dO+uxxx7TrbfeqrAwlh8A1sBvIwABoXPnzg3/3bt3b7Vp00ZvvvmmSktLVVpaqpKSEvXv31+SVFpaqj59+lzw+LFjx0qSXn75Ze3fv18333xzwzFjjM6ePdtwqYMkXXvttY3mkpmZqZ///OeKiYnRrbfeqp/+9KdKS0tz+VxCQ0MvyK93796qra1VWVmZ/vu//1v5+fkNx959991L5tSyZcuLbtfW1jbEc/T6XErr1q0v2QhLnr2eycnJ+vDDD5WQkKAf/ehH+slPfqKRI0eqZcuWDl8jAPAXmmEAAeH85umzzz7Tww8/rOHDh+uWW27RuHHjtGrVKu3atUuSFB4e3ugfidXV1emmm27SnDlzLjoWFRWlzz//XNJ3u6uNufXWW7Vx40Z9+umn+vTTT7VgwQJ9+OGH+tOf/qSQkJBLPuf3hYaGNvz3uVzDw8M1efJkTZw48ZI5nf8YSWrR4tJXujl7fS4lIiKi0WOevJ4tW7bUX/7yF23btk2bNm1Sfn6+3n//fX3wwQfq1atXo88JAP5CMwwg4CxevFjDhg274DrbN998s6Fhi46OVmFh4QWPef7559W2bVt169ZNq1evVseOHRsa3qKiIr377ruaO3euS8//X//1X+ratauGDx+u4cOHq6ysTAkJCdq1a5fCw8N1+vRpGWMaGuMDBw4oMjKy4fH19fX68ssvdcMNN0iS/vGPf+iyyy7Ttddeq8suu0xXXXWV+y+OnL8+l2rYHfHk9dy2bZv27t2r8ePH64477tCUKVM0ZMgQffrppzTDACyBP6ADYDlt2rRRaWmpjhw5csnjkZGRKi4u1r///W/t3btX8+fP16efftowwSA9PV0lJSV67bXXVFZWpj//+c9auXKl4uLilJKSohYtWigrK0tffvml/v73v2v69Omqra3V5Zdf7lJ+hw4d0osvvqi//vWv2r9/v5YvX662bduqa9euuvHGG1VdXa13331X+/fv1zvvvKMvvvjigseHhIRo+vTpKi4uVmFhoV5++WWNGzfO4W50Uzh7fVq3bi1J2rlzp8NZw+d48noaY/S73/1Oq1at0tdff62PP/5Yhw8fvuiyCwBoLjTDACxn3Lhx+vOf/3zB5QLne/LJJ9WlSxeNHTtW6enp+vLLL5WVlaXdu3erurpanTt31ltvvaVNmzYpOTlZb775pl566SXFxMSodevWysnJ0YkTJ/Szn/1Mjz76qG666Sa98sorLuf31FNPKTY2Vk8//bSGDx+uwsJCvf3222rXrp2io6M1ffp0LV68WCkpKSotLdUDDzxwweNbtWqlkSNH6qGHHtITTzyhu++++6LRa55w9vr8x3/8h372s5/p+eef1+uvv+40niev509+8hNlZWUpOztbCQkJmj9/vl544QXFxMR47XwBwBMhxpXp6wAAAEAQYmcYAAAAtkUzDAAAANuiGQYAAIBt0QwDAADAtmiGAQAAYFs0wwAAALAtmmEAAADYFs0wAAAAbItmGAAAALZFMwwAAADbohkGAACAbdEMAwAAwLZohgEAAGBbNMMAAACwLZphAAAA2BbNMAAAAGyLZhgAAAC2RTMMAAAA26IZBgAAgG3RDAMAAMC2aIYBAABgWzTDAAAAsC2aYQAAANgWzTAAAABsi2YYAAAAtkUzDAAAANuiGQYAAIBt0QwDAADAtmiGAQAAYFs0wwAAALAtmmEAAADYFs0wAAAAbItmGAAAALZFMwwAAADbohkGAACAbdEMAwAAwLZohgEAAGBbNMMAAACwLZphAAAA2BbNMAAAAGyLZhgAAAC2RTMMAAAA26IZBgAAgG3RDAMAAMC2aIYBAABgWzTDAAAAsC2aYQAAANgWzTAAAABsK6y5EwAAfKeqqkoLFy7Up59+qrq6Og0YMECPPvqoIiMjmzs1AJdAzQYHdoaD2MaNG3XDDTfohz/8oX70ox+psLCwuVMC0IiqqiplZGSoTZs2WrRokf7whz+oR48eysjIUEVFRXOnB+B7qNngEWKMMc2dBHxj5cqVmjlzpnr06KFFixbp8ssvb+6UADRi3rx5atOmjR577LEL7l++fLk2bdqk7OzsZsoMwKVQs8GDneEgl5CQoJtuukkTJ05UVVWVw5+dPHmyJk6c2HB7/vz5uuOOO3ydIgBJBQUFuvfeeyV99686GzZskCQlJyfrr3/9a8PP3XfffTp16lTD7eTkZP8mCkCS6zX7+OOPa+/evc2SI1xDMxzEunTpov79++u5555zqSGuqKhQixbffSSOHTumoqIide/e3V/pArZWX1+v8PBwSdK6dev08ccfS5JCQ0MbfsYYo9OnT6t169aSpPLycnXs2NH/yQJwqWYlad++ferSpYvf84PraIaD2M0336yf/exnkuS0Ia6pqVFYWJjatGmj6upq5eTkKDY2Vj169PB32oAt3XrrrVqzZo0k6aWXXtLcuXMlfbfj1Lt3b0nSgQMHdM011zQ8pqSkhBoFmokrNVtdXa2wsLCLGmRYC82wjThqiPfs2aPo6Gh17dpVRUVF2rt3ry6//HIWWsBPHnvsMS1evFjLly9XXV2dzp49qw0bNug3v/mNpk6dKkn68ssv1bNnz4bHfP82AP9xpWb37Nmjrl27NnOmcIZm2Gaee+45JSQkqLa29oL7z+0wde/eXS+++KLGjx+vkpISLpMA/KRDhw76/e9/r4KCAsXFxSk2Nla5ubl66623Gv6n9PsL6xdffNGwAwXAv1yp2S+//JJ1NAAwZ9iGzv8juXNKSkp066236rrrrlNcXJz69++vBQsWsOsE+FHHjh316quvNno8PDxcx48flyQdOnRI5eXl/OsN0Iyc1ezu3bvVr18/P2YEdzBaDZKkRx55RDNmzNAPfvCDhvtGjBih5cuXN2NWAM53/PhxPfPMM5KkkJAQTZkyRTfccEMzZwWgMY888oi++eYbtW3bVpKUnZ2tK664opmzwvfRDAMAAMC2uGYYAAAAtkUzDAAAANuiGQYAAIBt0QwDAADAtvw+Wq24pPKi+/r2iLTMMavlw3lYKx9fn4eV8H5Y/5jV8rHTeVgN71XwnofV8gnE83DGaTNcWlqqtWvX6uDBg2rRooUiIyMVGxurPn36uPQEAPyDWgUCB/UKWIfDyySWLFnSMNOyT58+Dd909MILLygnJ8f32QFwCbUKBA7qFbAWhzvDixcvVl5enlq1anXB/RMmTFBaWpoyMzN9mhwA11CrQOCgXgFrcbgzHBYWprq6uovuP3PmjMLDw32WFICmoVaBwEG9Atbi8Bvo8vPz9dprrykmJkYdOnSQJB0+fFhbt27V008/rcTERL8lCqBx1CoQOKhXwFqcfh3zoUOHVFhYqMrKShljFBUVpZiYGEVFRfkrRwAuoFaBwEG9AtbhtBn2NkejL/4SfsNFxxJrdzk95u6ojUvF9EZcK5yjN3L19/kH0nl48rly9BpYSbC8H5yH9Y9ZLR9GqwXXexXo52G1fALxPJzhSzcAAABgWzTDAAAAsC2aYQAAANgWzTAAAABsi2YYAAAAtuX3aRIAAACAVbAzDAAAANsKc/YD69evV0VFheLi4tSlS5eG+5cuXapRo0Y1+Ql9MYPX3fm8nsyD9fec4WCZz9vUczz/uJXOw4pzhv1Rq5I1X/NAmnfJeQTHeXjKn2trc/wus9J75ewz5+/vOLBDfVjtPJxxuDP8yiuv6IMPPlBZWZlGjx6tlStXNhz78MMPXXoCAL5HrQKBg3oFrMXhzvCmTZu0YsUKhYWFKSMjQ5mZmYqIiNCwYcPEpcaAdVCrQOCgXgFrcdgMG2MUEhIiSYqOjtbbb7+tCRMm6Morr2y4H0Dzo1aBwEG9Atbi8DKJoUOHKiMjQ8XFxZKkHj16KDs7W0899ZT27dvnlwQBOEetAoGDegWsxelotcLCQkVGRqpbt24N91VUVCgnJ0e//OUvfZ4gANdQq0DgoF4B62DOMAAAAGzL6Wg1bwukkWT+HrfizXNszlwZ4+P+58pKrDYaJ1hG/HAewXEeVmO198rq6443ntNfa5Inj7XSsebMxxm+dAMAAAC2RTMMAAAA26IZBgAAgG3RDAMAAMC2mCYBAAAA22JnGAAAALbVpNFqc+fO1bRp0zx6Ql+MzGiOkWRWH7vm63yaY/yNL94Pd3P1ZGSfv8Y1eVqvgT6OyNcjfqxwjp481krHrJZPc4xB9EW9WvW9ssPvAKvlaoVjzZmPM402w9OnT7/ovg0bNuj48eOSpDlz5rj0BAB8j3oFAgf1ClhLo81w+/btlZeXp0ceeUTt2rWTJG3dulW33Xab35ID4BrqFQgc1CtgLY1eM5yVlaVXX31Vq1evVqdOnZSWlqYrrrhCaWlpSktL82eOAJygXoHAQb0C1uLwmuGYmBj98Ic/1IwZM7Rx40bV19f7Ky8ATUS9AoGDegWsw+XRan/84x+1Zs0a5eTk+DonAB6iXoHAQb0CzYs5wwAAALCtJo1W8wYrjcZpjtFq/hq75uvnbI5RLFZ6Pzz5XPlrtJqnGElm/XP0Rq5WOGa1fPw9BtEbeK+C9zyslo+v1/PmqFe+dAMAAAC2RTMMAAAA26IZBgAAgG3RDAMAAMC2mCYBAAAA22JnGAAAALbldLTa5s2b1a9fP7Vr1055eXkqLi5W7969NXLkSLee0ArjPbwxksxfI7l8PXLJ6mPXvBHXCmPwXM3VE/6oVcle74fVj1ktH1+PHrTSMU8F6trqq7jUTnDk6uu+xN1jzjjcGX7ppZf09ttvq7q6Wq+99pry8/PVvXt3rVu3TrNnz3bpCQD4HrUKBA7qFbAWhzvDW7ZsUX5+vkJDQ7Vx40bl5uYqIiJCo0aNUlJSkr9yBOAEtQoEDuoVsBaHO8MtW7bUkSNHJElXXXWVTp06JUk6ffq0wsL8/uV1ABpBrQKBg3oFrMVh1T3++OO69957lZiYqOuvv14ZGRmKiYlRQUGBJk2a5K8cAThBrQKBg3oFrMXpaLX9+/dr/fr1+uqrr1RfX6+rr75agwYNUt++ff2VIwAXUKtA4KBeAetgzjAAAABsy+8XJ1lhvIevR6tZYVyZr+JadRSLlcZ8efK5spJgeT+oD2sd86Q+rHQeVsNnLnjPw2r5BOJ5OMOXbgAAAMC2aIYBAABgWzTDAAAAsC2aYQAAANgW0yQAAABgW+wMAwAAwLZsPVrNV6M/rDAayFdxA/H9YLSae/hcWf+Y1fJpzpGWVhvV5G985oL3PKyWTyCehzNOm+HCwkK1bNlSN998s3JycrR9+3bdeOONeuihhxQREeHSkwDwPWoVCBzUK2AdDpvhefPmaceOHaqrq1Pnzp0VEhKi9PR0bdiwQbNmzdLs2bP9lScAB6hVIHBQr4C1OGyGN2/erJUrV6qmpkZxcXEqKChQeHi4Bg4cqNTUVH/lCMAJahUIHNQrYC0O/4DOGKOqqiodPXpUZ86c0cmTJyVJZ86cUW1trV8SBOActQoEDuoVsBaHO8MPPvig4uPjZYzRlClTlJmZqZiYGBUWFmrkyJH+yhGAE9QqEDioV8BanM4ZPnPmjOrr69WmTRvt2rVLBQUF6tWrl+68805/5QjABdQqEDioV8A6+NINAAAA2BZzhi+huebaBtJ5NMf74a+ZwK68H0193PmPDZTZpXb5XFkhH84jcH63WrFWpcB6r6xwzGr5BEuuzVnnnszw5xvoAAAAYFs0wwAAALAtmmEAAADYFs0wAAAAbItpEgAAALAtdoYBAABgW5YareaLEVhNfT5vxLXCOXojVyuMlvNVXF++H84+V56Mf/GnQPpcBUt9cB7++13vaUyr4TMXvOdhtXyseh6MVgMAAADc4HRnePPmzfroo4908OBBtWjRQpGRkRo4cKASEhL8kR+AJqBegcBBvQLW4LAZzs7OVnFxsVJSUhQZ+d02dGVlpZYtW6aioiJlZWX5JUkAzlGvQOCgXgHrcNgMr169WmvWrFGLFhdeTZGUlKSkpCSKFbAQ6hUIHNQrYB0OR6ulpKTorbfeUqdOnS64f//+/Xr88ce1cuVKnycIwDXUKxA4qFfAOhzuDE+bNk3333+/oqOj1aFDB0nS4cOHVVZWpjlz5vglQQCuoV6BwEG9Atbh9Es3qqurVVxcrMrKShljFBUVpX79+ikiIsJfOQJwEfUKBA7qFbAGh81weXm5wwd//593ADQf6hUIHNQrYB0Om+Hk5GSVlZUpMjJS3/+xkJAQffLJJz5PEIBrqFcgcFCvgHU4bIZPnjypMWPGaMaMGerfv78/8wLQRNQrEDioV8A6HH4DXdu2bTV79mzl5eX5Kx8AbqJegcBBvQLW4fQP6AAAAIBg5XBnGAAAAAhmNMMAAACwLZphAAAA2FazNMP5+fkaPny4hgwZoiVLlnglZnZ2toYPH67ExES99957Xom5YcMGjRgxQkOHDtXs2bO9ElOS3nnnHSUkJCg5OVkLFy70KNbJkyeVlJSkAwcOSJKWLl2qpKQkJScna/r06aqpqfFK3OnTpys+Pl6pqalKTU3VunXrPI5ZUFCglJQUJSUlaerUqU3O9Y033lBiYqISExM1b968hvtra2v1wAMPaNu2bU3O0VFcSVqyZIkyMjI8jrlp06aG1zI1NVUDBgzQww8/7Fa+vuSLWpWoV2/Wqy9q9VJxqVd71qsvalXyTb1avVYvFdeKa6vkm3r1Ra1eKq5b9Wr87ODBg2bQoEHm6NGj5ttvvzXJycmmpKTEo5jbtm0zo0ePNrW1teb06dNm0KBBprS01KOY+/btM3fddZepqKgwNTU1Jj093WzcuNGjmMYYs2XLFpOUlGSqqqpMXV2defjhh83atWvdilVUVGSSkpJM7969zf79+82ePXvMkCFDTFVVlTl79qyZOnWqee+99zyOa4wxSUlJ5tChQ27l2VjMgQMHmt27dxtjjHniiSdMbm6uy/G2bNliRo0aZaqrq01NTY0ZN26c+fjjj01paakZNWqU6dOnj9m6dWuT82wsrjHGlJSUmNjYWDN27FivxTTGmMrKSjN48GCzd+/eJufrS76oVWOoV2/Wqy9qtbG41Ot37FSvvqhVY3xTr1av1UvFNcZ6a6sxvqlXX9Sqs7jGuF6vft8Z/uyzzzRgwAC1b99erVu3VkJCgj766COPYt52221avHixwsLCdOTIEdXX16t169YexVy3bp2GDx+ujh07Kjw8XPPnz1e/fv08iilJX3zxhe666y61bdtWoaGhio2N1fr1692KlZubqxkzZigyMlKSFBERoZkzZ6pt27YKCQlRz549nX7LkStxT506pfLycr3wwgtKTk7W66+/rrNnz3oUU5Lq6+t18uRJ1dfXq7q6WpdddpnL8Tp06KBp06YpIiJC4eHh6tatm8rLy7Vs2TJNmjTJ7feqsbg1NTX61a9+pcmTJ3st5jnz5s3T6NGjFR0d7VbOvuKLWpWoV2/Wqy9q9VJxJer1HDvVqy9qVfJNvVq9Vi8V14prq+SbevVFrTqKe46r9er3ZriyslIdOnRouB0ZGalDhw55HDc8PFyvv/66EhMTFRMTo6ioKI/iffXVV6qvr9fEiROVkpKiP/zhD7riiis8zrN3794qKCjQsWPHVF1drQ0bNuibb75xK9ZLL72kW265peH2NddcozvuuEOS9L//+79asmSJBg8e7HHcI0eOaMCAAfrNb36j3Nxc7dixQ8uWLfMopiTNnDlTGRkZio2N1dGjRzV06FCX4/Xo0UM33XSTJKmsrEyrV69WXFycpk6dqrvvvrtJubkS93e/+51Gjhypzp07ey3mudvbt2/XuHHj3M7ZV3xVqxL1KnmnXn1Rq5eKK1Gv527brV69XauSb+rV6rV6qbhWXFsl39SrL2rVUdxzt12tV783w+YSY41DQkK8EvvJJ59UYWGhKioqlJub61Gs+vp6FRYW6re//a1yc3P1z3/+UytWrPA4x5iYGI0YMUIZGRmaNGmS+vfvr/DwcI/jnu/QoUN64IEHNHLkSN1+++0ex7v22mv15ptv6qqrrlKrVq2UkZGhTZs2eRTz8OHDeuWVV7Rq1SoVFBSoX79+mjNnTpPjlJSUKDMzU1lZWV7dqTk/7tdff62KigqNHDnSazHP5bp06VKNGTNGERERXsjau3xZqxL1eo4369UXtSpRr3avV2/WquSbeg20WpWsvbZKvqlXX9Tq9+O6U69+b4ajoqIu+L+1ysrKC7b33VFaWqqdO3dKklq1aqX4+Hjt2rXLo5hXX321YmJidOWVV6ply5YaPHiwiouLPYopfXeh+5AhQ5Sfn6/3339frVq10rXXXutx3HNKS0uVnp6utLQ0PfbYY16JuWvXLq1du7bhtjFGYWFhHsXcsWOHevbsqS5duqhFixa67777tH379ibF+Nvf/qbx48fr2WefVVpamkf5OIq7atUqlZSUKDU1Vc8//7z+9a9/6amnnvJKrp988omGDx/utdy9yRe1KlGv5/N2vfqiViXq9Ry71asvalXyTb0GWq1K1l1bJd/Uqy9q1VGuTarXJl+t7KFzF/kfOXLEnDp1yqSkpJjPP//co5gbN240I0aMMNXV1aa6utpMmDDBrFq1yqOYRUVFJiEhwRw/frzhYvymXoR+KTt37jQpKSmmtrbWnDhxwgwbNszs2LHDo5iDBg0y+/fvN1VVVSYuLs7k5eV5nOf5cXfu3GkGDhxojh07ZmpqakxmZqbJz8/3KObu3btNXFycOXz4sDHGmIULF5qsrCyX45SXl5vbb7/dfPbZZ5c8PnbsWLf+IMdZ3K1btzb5Iv/GYh45csTExsY2OUd/8UWtGkO9+qJefVGr58elXu1Zr76oVWN8U6+BUqvnx7Xi2mqMb+rVF7XqKG5T69XzLYMmioqK0tNPP61x48aptrZW9957r/r27etRzLi4OH3++ee65557FBoaqvj4eCUmJnoUs1+/fpo0aZLGjBmj2tpa3XnnnV7Zyu/Vq5fi4+OVkpKi+vp6jR8/Xv379/c4riQtW7ZM33zzjXJycpSTkyNJ+ulPf+r2henn9OrVSw899JDS09NVV1en+Ph4JSUleRSzW7dumjx5ssaNG6fQ0FBdd911mjVrlsuPX7RokaqrqzV37tyG+0aPHq309HSP8vJF3MZi9u7dWx07dvQoX1/yRa1K1Os5vqhXX9SqRL3atV59UauSb+o10GpVsubaKvm3rnz1O6Cp9RpizCUuNAIAAABsgG+gAwAAgG3RDAMAAMC2aIYBAABgWzTDAAAAsC2aYQAAANgWzTAAAABsi2YYAAAAtkUzDAAAANuiGQYAAIBt0QwDAADAtmiGAQAAYFs0wwAAALAtmmEAAADYFs0wAAAAbItmGAAAALZFMwwAAADbohkGAACAbdEMAwAAwLZohgEAAGBbNMMAAACwLZphAAAA2BbNMAAAAGyLZsVwGuUAACAASURBVBgAAAC2RTMMAAAA26IZBgAAgG3RDAMAAMC2aIYBAABgWzTDAAAAsC2aYQAAANgWzTAAAABsi2YYAAAAtkUzDAAAANuiGQYAAIBt0QwDAADAtmiGAQAAYFs0wwAAALAtmmEAAADYFs0wAAAAbItmGAAAALZFMwwAAADbohkGAACAbdEMAwAAwLZohgEAAGBbNMMAAACwLZphAAAA2BbNMAAAAGyLZhgAAAC2RTMMAAAA26IZBgAAgG3RDAMAAMC2aIYBAABgWzTDAAAAsC2aYQAAANgWzTAAAABsi2YYAAAAtkUzbBOvvvqqXn311eZOA4ALqFcA8J+w5k4Avjd37lwtXrxYklRbW6usrKxmzghAY6hXAPAvmuEgt379eu3YsUNjx46VJG3fvl3r16/X3Xff3cyZAfg+6hUA/C/EGGOaOwn41qlTp7Ro0SJJ0sSJE9W6detmzghAY6hXAPAvrhm2gfMXU0cL6+TJkzVx4sSG2/Pnz9cdd9zh09wAXMjVer3vvvt06tSphtvJyck+zQsAghXNMBpUVFSoRYvvPhLHjh1TUVGRunfv3sxZAfg+Y4xOnz7d0CyXl5erY8eOzZwVAAQmrhmGJKmmpkZhYWFq06aNqqurlZOTo9jYWFVUVDR3agC+58CBA7rmmmsabpeUlKhHjx7NmBEABC52hiFJ2rNnj6Kjo9W1a1cVFRVp7969uvzyy1lgAQv68ssv1bNnz0ZvAwBcRzMMSf+3s9S9e3e9+OKLGj9+vEpKSrhMArCgPXv2qGvXrg23v/jiC/Xu3bsZMwKAwMVlEjbxxBNPODxeUlKiW2+9Vdddd53i4uLUv39/LViwgN0moBk4q9fw8HAdP35cknTo0CGVl5fzrzgA4CaaYUj67p9Z09PT9YMf/EBTpkyRJJ04cULt2rVr5swAfF9aWpqeeeYZbd68WSEhIZo1a1ZzpwQAAYs5wwAAALAtrhkGAACAbdEMAwAAwLZohgEAAGBbNMMAAACwLb9Pkyguqbzovr49Ii1zzGr5cB7WysfX52ElvB/WP2a1fOx0HgCCh9NmuLS0VGvXrtXBgwfVokULRUZGKjY2Vn369PFHfgBcRK0CANB0Di+TWLJkiZ555hlJUp8+fRq+4eiFF15QTk6O77MD4BJqFQAA9zjcGV68eLHy8vLUqlWrC+6fMGGC0tLSlJmZ6dPkALiGWgUAwD0Od4bDwsJUV1d30f1nzpxReHi4z5IC0DTUKgAA7nG4M/zII4/onnvuUUxMjDp06CBJOnz4sLZu3aqnn37aLwkCcI5aBQDAPU6/jvnQoUMqLCxUZWWljDGKiopSTEyMoqKi/JUjABdQqwAANJ3TZhgAAAAIVswZvgQr5cN5WCsf5gxb63Xlc2WtfOx0HgCCB99ABwAAANuiGQYAAIBt0QwDAADAtmiGAQAAYFtMkwAAAIBtsTMMAAAA23I6Wm39+vWqqKhQXFycunTp0nD/0qVLNWrUqCY/oZVG4wTLiB/Ow/rHXH2sJ/xRq5K1Xlc+V9bKx07nASB4ONwZfuWVV/TBBx+orKxMo0eP1sqVKxuOffjhhz5PDoBrqFUAANzjcGd406ZNWrFihcLCwpSRkaHMzExFRERo2LBh4lJjwDqoVQAA3OOwGTbGKCQkRJIUHR2tt99+WxMmTNCVV17ZcD+A5ketAgDgHoeXSQwdOlQZGRkqLi6WJPXo0UPZ2dl66qmntG/fPr8kCMA5ahUAAPc4Ha1WWFioyMhIdevWreG+iooK5eTk6Je//KXPEwTgGmoVAICmY84wAAAAbMvpaDVvc3eMzV/Cb7joWGLtLo9iWm3EjzfP0de5MnLJ/fNw9D5biV3eDyvkw3lYKx9GqwH2wpduAAAAwLZohgEAAGBbNMMAAACwLZphAAAA2BbTJAAAAGBb7AwDAADAtgJmtJq7I8maOq7MG89phbFr3nhOKxyzWj6+Pg8rCaT6aI7PlRXO0VdxA/H98PfvXQDBo0k7w3PnzvVVHgC8jHoFAMC5RneGp0+fftF9GzZs0PHjxyVJc+bM8V1WAJqEegUAwD2NNsPt27dXXl6eHnnkEbVr106StHXrVt12221+Sw6Aa6hXAADc0+hlEllZWXr11Ve1evVqderUSWlpabriiiuUlpamtLQ0f+YIwAnqFQAA9zgdrXbs2DHNmDFDnTp1UkFBgfLz8/2VG4Amol4BAGgap39A1759e2VnZ+v6669Xhw4d/JETADdRrwAANI3fv3TDSqNxAmnET1PPUbLWKCtvjFzy5gg9q46OshJGkln/HL2RqxWOWS2fQKtVAJ7hSzcAAABgWzTDAAAAsC2aYQAAANgWzTAAAABsy+9/QAcAAABYBTvDAAAAsK1Gv47ZV6w0GscuI8mskI+vz8PdsWtWG4FlJc0xWs1q74fVj1ktHzuNtAQQPJzuDG/evFknTpyQJOXl5WnWrFn605/+5PPEADQNtQoAQNM5bIZfeuklvf3226qurtZrr72m/Px8de/eXevWrdPs2bP9lSMAJ6hVAADc4/AyiS1btig/P1+hoaHauHGjcnNzFRERoVGjRikpKclfOQJwgloFAMA9DneGW7ZsqSNHjkiSrrrqKp06dUqSdPr0aYWF+f1yYwCNoFYBAHCPw9FqGzZs0MyZM5WYmKi6ujpt3bpVMTExKigo0KRJkzRixAh/5gqgEdQqAADucTpneP/+/Vq/fr2++uor1dfX6+qrr9agQYPUt29ff+UIwAXUKgAATef3L92w0mgcZyN+/DXmyU4jl6x2Hr4YyebJyD4rCaT3w2qfq2Cpj+YYaWmF0XuMVgPshS/dAAAAgG3RDAMAAMC2aIYBAABgWzTDAAAAsC2aYQAAANiW36dJAAAAAFbBzjAAAABsy+n3tBYWFqply5a6+eablZOTo+3bt+vGG2/UQw89pIiIiCY/oZXmXQbL3E7Ow1ozb5trzrA/alWy1ucjED9Xdj8PX80Z9vcMYgDBw2EzPG/ePO3YsUN1dXXq3LmzQkJClJ6erg0bNmjWrFmaPXu2v/IE4AC1CgCAexw2w5s3b9bKlStVU1OjuLg4FRQUKDw8XAMHDlRqaqq/cgTgBLUKAIB7HF4zbIxRVVWVjh49qjNnzujkyZOSpDNnzqi2ttYvCQJwjloFAMA9DneGH3zwQcXHx8sYoylTpigzM1MxMTEqLCzUyJEj/ZUjACeoVQAA3ON0tNqZM2dUX1+vNm3aaNeuXSooKFCvXr105513+itHAC6gVgEAaDrmDAMAAMC2nI5W8zYrjfgJllFFnIc1xkM192g1bwuk9yNYPld2OA9PRqv5YhQio9UA8KUbAAAAsC2aYQAAANgWzTAAAABsi2YYAAAAtsU0CQAAANgWO8MAAACwLUarXYKV8uE8/D8GzRdjnFx9Tivhc2X9Y1bLxxuj1XxR5744DwDBg51hAAAA2JbTneHNmzfro48+0sGDB9WiRQtFRkZq4MCBSkhI8Ed+AJqAegUAoGkcNsPZ2dkqLi5WSkqKIiO/+6ehyspKLVu2TEVFRcrKyvJLkgCco14BAGg6h83w6tWrtWbNGrVoceHVFElJSUpKSmJxBSyEegUAoOkcjlZLSUnRW2+9pU6dOl1w//79+/X4449r5cqVPk8QgGuoVwAAms7hzvC0adN0//33Kzo6Wh06dJAkHT58WGVlZZozZ45fEgTgGuoVAICmc/qlG9XV1SouLlZlZaWMMYqKilK/fv0UERHhrxwBuIh6BQCgaRw2w+Xl5Q4f/P1/jgXQfKhXAACazmEznJycrLKyMkVGRur7PxYSEqJPPvnE5wkCcA31CgBA0zlshk+ePKkxY8ZoxowZ6t+/vz/zAtBE1CsAAE3n8Bvo2rZtq9mzZysvL89f+QBwE/UKAEDTOf0DOgAAACBYOdwZBgAAAIIZzTAAAABsi2YYAAAAttUszXB+fr6GDx+uIUOGaMmSJV6JmZ2dreHDhysxMVHvvfeeV2Ju2LBBI0aM0NChQzV79myvxJSkd955RwkJCUpOTtbChQs9inXy5EklJSXpwIEDkqSlS5cqKSlJycnJmj59umpqarwSd/r06YqPj1dqaqpSU1O1bt06j2MWFBQoJSVFSUlJmjp1apNzfeONN5SYmKjExETNmzev4f7a2lo98MAD2rZtW5NzdBRXkpYsWaKMjAyPY27atKnhtUxNTdWAAQP08MMPu5WvL/miViXq1Zv16otavVRc6tX69QrATcbPDh48aAYNGmSOHj1qvv32W5OcnGxKSko8irlt2zYzevRoU1tba06fPm0GDRpkSktLPYq5b98+c9ddd5mKigpTU1Nj0tPTzcaNGz2KaYwxW7ZsMUlJSaaqqsrU1dWZhx9+2Kxdu9atWEVFRSYpKcn07t3b7N+/3+zZs8cMGTLEVFVVmbNnz5qpU6ea9957z+O4xhiTlJRkDh065FaejcUcOHCg2b17tzHGmCeeeMLk5ua6HG/Lli1m1KhRprq62tTU1Jhx48aZjz/+2JSWlppRo0aZPn36mK1btzY5z8biGmNMSUmJiY2NNWPHjvVaTGOMqaysNIMHDzZ79+5tcr6+5ItaNYZ69Wa9+qJWG4tLvX7HqvUKwH1+3xn+7LPPNGDAALVv316tW7dWQkKCPvroI49i3nbbbVq8eLHCwsJ05MgR1dfXq3Xr1h7FXLdunYYPH66OHTsqPDxc8+fPV79+/TyKKUlffPGF7rrrLrVt21ahoaGKjY3V+vXr3YqVm5urGTNmKDIyUpIUERGhmTNnqm3btgoJCVHPnj2dfiuZK3FPnTql8vJyvfDCC0pOTtbrr7+us2fPehRTkurr63Xy5EnV19erurpal112mcvxOnTooGnTpikiIkLh4eHq1q2bysvLtWzZMk2aNMnt96qxuDU1NfrVr36lyZMney3mOfPmzdPo0aMVHR3tVs6+4otalahXb9arL2r1UnEl6vUcq9YrAPf5vRmurKxUhw4dGm5HRkbq0KFDHscNDw/X66+/rsTERMXExCgqKsqjeF999ZXq6+s1ceJEpaSk6A9/+IOuuOIKj/Ps3bu3CgoKdOzYMVVXV2vDhg365ptv3Ir10ksv6ZZbbmm4fc011+iOO+6QJP3v//6vlixZosGDB3sc98iRIxowYIB+85vfKDc3Vzt27NCyZcs8iilJM2fOVEZGhmJjY3X06FENHTrU5Xg9evTQTTfdJEkqKyvT6tWrFRcXp6lTp+ruu+9uUm6uxP3d736nkSNHqnPnzl6Lee729u3bNW7cOLdz9hVf1apEvUreqVdf1Oql4krU67nbVq1XAO7zezNsLjHWOCQkxCuxn3zySRUWFqqiokK5ubkexaqvr1dhYaF++9vfKjc3V//85z+1YsUKj3OMiYnRiBEjlJGRoUmTJql///4KDw/3OO75Dh06pAceeEAjR47U7bff7nG8a6+9Vm+++aauuuoqtWrVShkZGdq0aZNHMQ8fPqxXXnlFq1atUkFBgfr166c5c+Y0OU5JSYkyMzOVlZXl1Z2a8+N+/fXXqqio0MiRI70W81yuS5cu1ZgxYxQREeGFrL3Ll7UqUa/neLNefVGrEvUaCPUKwH1+b4ajoqIu2FmprKy84J/i3FFaWqqdO3dKklq1aqX4+Hjt2rXLo5hXX321YmJidOWVV6ply5YaPHiwiouLPYopffdHKUOGDFF+fr7ef/99tWrVStdee63Hcc8pLS1Venq60tLS9Nhjj3kl5q5du7R27dqG28YYhYWFeRRzx44d6tmzp7p06aIWLVrovvvu0/bt25sU429/+5vGjx+vZ599VmlpaR7l4yjuqlWrVFJSotTUVD3//PP617/+paeeesoruX7yyScaPny413L3Jl/UqkS9ns/b9eqLWpWo13OsXK8APODvi5TP/VHOkSNHzKlTp0xKSor5/PPPPYq5ceNGM2LECFNdXW2qq6vNhAkTzKpVqzyKWVRUZBISEszx48cb/nCmKX8w0pidO3ealJQUU1tba06cOGGGDRtmduzY4VHMQYMGmf3795uqqioTFxdn8vLyPM7z/Lg7d+40AwcONMeOHTM1NTUmMzPT5OfnexRz9+7dJi4uzhw+fNgYY8zChQtNVlaWy3HKy8vN7bffbj777LNLHh87dqxbf5DjLO7WrVub/Ac5jcU8cuSIiY2NbXKO/uKLWjWGevVFvfqiVs+PS71av14BuM/zLYMmioqK0tNPP61x48aptrZW9957r/r27etRzLi4OH3++ee65557FBoaqvj4eCUmJnoUs1+/fpo0aZLGjBmj2tpa3XnnnR7/s5sk9erVS/Hx8UpJSVF9fb3Gjx+v/v37exxXkpYtW6ZvvvlGOTk5ysnJkST99Kc/deuPSM7Xq1cvPfTQQ0pPT1ddXZ3i4+OVlJTkUcxu3bpp8uTJGjdunEJDQ3Xddddp1qxZLj9+0aJFqq6u1ty5cxvuGz16tNLT0z3KyxdxG4vZu3dvdezY0aN8fckXtSpRr+f4ol59UasS9RoI9QrAfSHGXOLCQAAAAMAG+AY6AAAA2BbNMAAAAGyLZhgAAAC2RTMMAAAA26IZBgAAgG3RDAMAAMC2aIYBAABgWzTDAAAAsC2aYQAAANgWzTAAAABsi2YYAAAAtkUzDAAAANuiGQYAAIBt0QwDAADAtmiGAQAAYFs0wwAAALAtmmEAAADYFs0wAAAAbItmGAAAALZFMwwAAADbohkGAACAbdEMAwAAwLZohgEAAGBbNMMAAACwLZphAAAA2BbNMAAAAGyLZhgAAAC2RTMMAAAA26IZBgAAgG3RDAMAAMC2aIYBAABgWzTDAAAAsC2aYQAAANgWzTAAAABsi2YYAAAAtkUzDAAAANuiGQYAAIBt0QwDAADAtmiGAQAAYFs0wwAAALAtmmEAAADYFs0wAAAAbItmGAAAALZFMwwAAADbohkGAACAbdEMAwAAwLZohgEAAGBbNMMAAACwLZphAAAA2BbNMAAAAGyLZhgAAAC2RTMMAAAA26IZBgAAgG3RDAMAAMC2aIYBAABgWzTDAAAAsC2aYQAAANgWzXCQM8Y0dwoAAAQl1tjgQDMcxOrr6/Xss89q7dq1zZ0KAABBhTU2eNAMB6n6+nr94he/0LfffqtBgwa59JjJkydr4sSJDbfnz5+vO+64w1cpAgAQkNxZYx9//HHt3bvXx5nBHWHNnQB8Y82aNVq9erUkqU+fPpf8mTlz5mjEiBENtysqKnTFFVdIko4dO6aioiJ1797d98kCABBA3Flj9+3bpy5duvglPzQNzXCQGjZsmD755BOdOnVKCxYsUEREhMOfr6mpUVhYmNq0aaPq6mrl5OQoNjZWFRUVfsoYAIDA0NQ1trq6WmFhYQoNDfVThmgKLpMIUqGhoXrllVfUpk0b/c///I/Tn9+zZ4+io6PVtWtXFRUVae/evbr88svVo0cPP2QLAEDgcGeN7dq1qx8ygztCDH8KGdSMMQoJCXH6c/n5+frmm28UGRmphQsX6te//rXWrFmjoUOH6pZbbvFDpgAABBZX19iVK1eqvLxcP//5z/2QFZqKneEg50qRSlJJSYm6d++uPn36KC4uTv3799fu3bvVs2dPH2cIAEBgcnWN3b17N//SamE0w5Akffnll+revbu6dOmiKVOmSJJOnDihdu3aNXNmAAAEtpKSEr311lsaP368xo8fr+PHjzd3SjgPl0kAAADAttgZBgAAgG3RDAMAAMC2aIYBAABgWzTDAAAAsC2aYQAAANiW069jLi0t1dq1a3Xw4EG1aNFCkZGRio2NbfS7uJ0pLqm86L6+PSItc8xq+XAe1srH1+cBwB4CdW31VVyr/k62wjGr5ROI5+GMw53hJUuW6JlnnpEk9enTR71795YkvfDCC8rJyXHpCQAAwP9hbQWsxeHO8OLFi5WXl6dWrVpdcP+ECROUlpamzMxMnyYHAECwYW0FrMXhznBYWJjq6uouuv/MmTMKDw/3WVIAAAQr1lbAWhzuDD/yyCO65557FBMTow4dOkiSDh8+rK1bt+rpp5/2S4IAAAQT1lbAWpx+HfOhQ4dUWFioyspKGWMUFRWlmJgYRUVF+StHAACCCmsrYB1Om2EAAAAgWDkdreZtVhjvYdXRH5wH5wEA7rDa7zKrH7NaPsGSq1XPwxm+dAMAAAC2RTMMAAAA26IZBgAAgG3RDAMAAMC2mCYBAAAA22JnGAAAALbldLTa+vXrVVFRobi4OHXp0qXh/qVLl2rUqFFNfkJfjMz4S/gNFx1LrN3lVkxv5GOFc2zOXINlFAuj1QD4SiCsrc5+l3lz7bXq72QrnKOv4lr1NbfcaLVXXnlFH3zwgcrKyjR69GitXLmy4diHH37o0hMAAID/w9oKWIvDneFNmzZpxYoVCgsLU0ZGhjIzMxUREaFhw4aJS40BAGg61lbAWhw2w8YYhYSESJKio6P19ttva8KECbryyisb7gcAAK5jbQWsxeFlEkOHDlVGRoaKi4slST169FB2draeeuop7du3zy8JAgAQTFhbAWtxOlqtsLBQkZGR6tatW8N9FRUVysnJ0S9/+UufJwgAQLBhbQWsgznDAAAAsC2no9W8zQ4jyfw9iuVSz+fr57TDKBZGqwEIFFb7XeZoHQyWcWX+Og9v5GqFY82ZjzN86QYAAABsi2YYAAAAtkUzDAAAANuiGQYAAIBtMU0CAAAAtsXOMAAAAGwr6EerMZLMGmNsGOPDaDUAvmO132VW/33tyWOtug5a/Vhz5uNMk3aG586d25QfBwAALmB9BZpPozvD06dPv+i+DRs26Pjx45KkOXPm+C4rAACCFOsrYC2NNsPt27dXXl6eHnnkEbVr106StHXrVt12221+Sw4AgGDD+gpYS6OXSWRlZenVV1/V6tWr1alTJ6WlpemKK65QWlqa0tLS/JkjAABBg/UVsBano9WOHTumGTNmqFOnTiooKFB+fr6/cgMAIGixvgLW4PQP6Nq3b6/s7Gxdf/316tChgz9yAgAg6LG+Atbg9y/dsMJ4D0aSBc4xb8R1d4yPN99/V3MFAHcE0u9kd39fB9J52DVXb/RXvliXneFLNwAAAGBbNMMAAACwLZphAAAA2BbNMAAAAGyLZhgAAAC25fdpEgAAAIBVsDMMAAAA2wpz9gObN29Wv3791K5dO+Xl5am4uFi9e/fWyJEj3XpCK8y688Y8P6vPIG6O57TqefhinmFT339XcwVgD4G6tvoqbiCuLeQaOOu5Mw53hl966SW9/fbbqq6u1muvvab8/Hx1795d69at0+zZs116AgAA8H9YWwFrcbgzvGXLFuXn5ys0NFQbN25Ubm6uIiIiNGrUKCUlJfkrRwAAggZrK2AtDneGW7ZsqSNHjkiSrrrqKp06dUqSdPr0aYWFOb3CAgAAfA9rK2AtDqvu8ccf17333qvExERdf/31ysjIUExMjAoKCjRp0iR/5QgAQNBgbQWsxelotf3792v9+vX66quvVF9fr6uvvlqDBg1S3759/ZUjAABBhbUVsA7mDAMAAMC2/H5xkhXGezDCJHCOeSOuu2PwvDl2zdXnBAB3BNLvZCscs1o+wZJrUI5WAwAAAIIZzTAAAABsi2YYAAAAtkUzDAAAANtimgQAAABsi51hAAAA2JbT0WqFhYVq2bKlbr75ZuXk5Gj79u268cYb9dBDDykiIqLJT2iF8R6MMAmcY1bLx93xLs6Ouzr+BUBwCNS11VdxWVuCI1ernoczDpvhefPmaceOHaqrq1Pnzp0VEhKi9PR0bdiwQbNmzdLs2bNdehIAAPAd1lbAWhw2w5s3b9bKlStVU1OjuLg4FRQUKDw8XAMHDlRqaqq/cgQAIGiwtgLW4vCaYWOMqqqqdPToUZ05c0YnT56UJJ05c0a1tbV+SRAAgGDC2gpYi8Od4QcffFDx8fEyxmjKlCnKzMxUTEyMCgsLNXLkSH/lCABA0GBtBazF6Wi1M2fOqL6+Xm3atNGuXbtUUFCgXr166c477/RXjgAABBXWVsA6mDMMAAAA23I6Ws3brDDew6qjPzgP/59HU0eknTvW1JiuPhYA3BEsv5M5j8DO1Rvn4c2119WxpXzpBgAAAGyLZhgAAAC2RTMMAAAA26IZBgAAgG0xTQIAAAC2xc4wAAAAbIvRapdgpXw4j+bLx92xa85Gq3ky/gUALsUOv5Pteh5WyycQz8MZdoYBAABgWzTDAAAAsC2nl0ls3rxZH330kQ4ePKgWLVooMjJSAwcOVEJCgj/yAwAgKLG+AtbgsBnOzs5WcXGxUlJSFBn53XUXlZWVWrZsmYqKipSVleWXJAEACCasr4B1OBytlpCQoDVr1qhFiwuvpqivr1dSUpLWrFnj8wQBAAg2rK+AdTi8Zviyyy7TwYMHL7q/vLxcERERPksKAIBgxvoKWIfDyySmTZum+++/X9HR0erQoYMk6fDhwyorK9OcOXP8kiAAAMGG9RWwDqffQFddXa3i4mJVVlbKGKOoqCj169eP/3MFAMADrK+ANThshsvLyx0+uFOnTl5PCACAYMf6CliHw2Y4OTlZZWVlioyM1Pd/LCQkRJ988onPEwQAINiwvgLW4bAZPnnypMaMGaMZM2aof//+/swLAICgxfoKWIfDaRJt27bV7NmzlZeX5698AAAIeqyvgHU4/QM6AAAAIFg53BkGAAAAghnNMAAAAGyrWZrh/Px8DR8+XEOGDNGSJUu8EjM7O1vDhw9XYmKi3nvvPa/E3LBhg0aMGKGhQ4dq9uzZXokpSe+8844SEhKUnJyshQsXehTr5MmTSkpK0oEDByRJS5cuVVJSkpKTkzV9+nTV1NR4Je706dMVHx+v1NRUpaamat26dR7HLCgoUEpKipKSkjR16tQm5/rGG28oMTFRiYmJmjdvXsP9tbW1euCBB7Rt27Ym5+goriQtiuC34QAAB0ZJREFUWbJEGRkZHsfctGlTw2uZmpqqAQMG6OGHH3YrXwCQAmdtlXyzvrK2emdtlXyzvvpibb1UXLfWV+NnBw8eNIMGDTJHjx413377rUlOTjYlJSUexdy2bZsZPXq0qa2tNadPnzaDBg0ypaWlHsXct2+fueuuu0xFRYWpqakx6enpZuPGjR7FNMaYLVu2mKSkJFNVVWXq6urMww8/bNauXetWrKKiIpOUlGR69+5t9u/fb/bs2WOGDBliqqqqzNmzZ83UqVPNe++953FcY4xJSkoyhw4dcivPxmIOHDjQ7N692xhjzBNPPGFyc3NdjrdlyxYzatQoU11dbWpqasy4cePMxx9/bEpLS82oUaNMnz59zNatW5ucZ2NxjTGmpKTExMbGmrFjx3otpjHGVFZWmsGDB5u9e/c2OV8AMCZw1lZjfLO+srZ6Z201xjfrqy/WVmdxjXF9ffX7zvBnn32mAQMGqH379v+/vXsLafKPwwD+2NxqIgQdnFBWIIkgtKDIVq1R1iytlS7CSVtmo5sgiyC76CBBBzoQCNFVXhQFjUGBI4oOFEiGLMiKhphQRBMPKzsgvK633/9C5r+ibNvvfQPZ87l7Jz58YRfPdy+/7UVOTg7Ky8tx+/ZtqcwlS5bg8uXLyM7ORiwWg6qqyMnJkcq8e/cuKioqkJ+fD6PRiPPnz8NqtUplAsCrV6+wYsUK5ObmwmAwwG634969e2llBQIBHD16FHl5eQAAk8mEpqYm5ObmIisrC0VFRX/9YfdkcoeHhxGNRnH48GFs3LgRzc3N+P79u1QmAKiqiq9fv0JVVSiKgsmTJyedN3PmTBw8eBAmkwlGoxGFhYWIRqMIBoPw+/1pv1d/yh0ZGcGRI0fQ0NCgWWbC6dOnUVNTg3nz5qU1MxHRROlWQJ9+Zbdq062APv2qR7eOl5uQbL/+82W4v79/7DnsAJCXl4e+vj7pXKPRiObmZlRWVsJms8FisUjlvX37FqqqYufOnXC5XLh27RqmTp0qPWdJSQna2towNDQERVHw4MEDDA4OppV1/PhxLF68eOx61qxZWLZsGQDgw4cPuHr1KsrKyqRzY7EYli5dihMnTiAQCCAcDiMYDEplAkBTUxO8Xi/sdjs+fvyIdevWJZ03f/58LFy4EADw5s0b3Lp1Cw6HAwcOHMCaNWtSmi2Z3HPnzsHtdmP27NmaZSauOzo64PP50p6ZiGiidCugT7+yW/8n062APv2qR7eOl5u4TrZf//kyLH7zS25ZWVmaZO/Zswft7e3o7e1FIBCQylJVFe3t7Thz5gwCgQBevHiBGzduSM9os9lQXV0Nr9cLv9+PRYsWwWg0Suf+qK+vD9u3b4fb7UZpaal0XkFBAS5cuIDp06fDbDbD6/Xi0aNHUpkDAwM4e/YsQqEQ2traYLVacfLkyZRzuru7UV9fj8bGRk3vrP6Y+/79e/T29sLtdmuWmZj1+vXrqK2thclk0mBqIspUE6VbAX36ld06SqtuBfTpVz269dfcdPr1ny/DFovlp09r/f39P93eT0dPTw8ikQgAwGw2w+l0oqurSypzxowZsNlsmDZtGqZMmYKysjI8f/5cKhMYPei+du1atLa24sqVKzCbzSgoKJDOTejp6YHH40FVVRV2796tSWZXVxfu3Lkzdi2EQHZ2tlRmOBxGUVER5syZg0mTJmHr1q3o6OhIKePp06eoq6vD/v37UVVVJTXPeLmhUAjd3d3YtGkTDh06hJcvX2Lv3r2azHr//n1UVFRoNjsRZaaJ0q2APv3Kbh2lRbcC+vSrHt063qwp9WvKp5UlJQ75x2IxMTw8LFwul+js7JTKfPjwoaiurhaKoghFUcSOHTtEKBSSynz27JkoLy8Xnz59GjuMn+oh9N+JRCLC5XKJeDwuPn/+LNavXy/C4bBU5qpVq8S7d+/Ely9fhMPhEDdv3pSe88fcSCQiVq5cKYaGhsTIyIior68Xra2tUpmvX78WDodDDAwMCCGEuHjxomhsbEw6JxqNitLSUvH48ePf/n3btm1pfYHub7lPnjxJ+ZD/nzJjsZiw2+0pz0hE9KuJ0q1C6NOv7FZtulUIffpVj24dLzfVfpX7CJIGi8WCffv2wefzIR6PY8uWLViwYIFUpsPhQGdnJzZv3gyDwQCn04nKykqpTKvVCr/fj9raWsTjcSxfvlyTW/nFxcVwOp1wuVxQVRV1dXWaPZc+GAxicHAQLS0taGlpAQCsXr067YPpCcXFxdi1axc8Hg++ffsGp9OJDRs2SGUWFhaioaEBPp8PBoMBc+fOxbFjx5L+/0uXLkFRFJw6dWrstZqaGng8Hqm59Mj9U2ZJSQny8/Ol5iUiAiZOtwL69Cu7dZRstwL/tgf16uxU+5WPYyYiIiKijMUn0BERERFRxuIyTEREREQZi8swEREREWUsLsNERERElLG4DBMRERFRxuIyTEREREQZi8swEREREWUsLsNERERElLH+A7p9Gu5DxDsQAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x576 with 7 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(3,3,figsize = (10,8))\n",
    "\n",
    "for i,m in enumerate([[0,0,r\"$\\nwarrow_{M}$\"],[0,1,r\"$\\odot_{U}$\"],[0,2,r\"$\\odot_{L}$\"],\n",
    "                      [1,0,r\"$\\uparrow_{M}$\"],[1,1,r\"$\\uparrow_{U}$\"],\n",
    "                      [2,0,r\"$\\leftarrow_{M}$\"],[2,2,r\"$\\leftarrow_{L}$\"]]):\n",
    "    sns.heatmap(tracer[:,:,i],cmap=\"coolwarm\",vmin=-0.5,vmax=0.6,linecolor='white',\n",
    "                linewidth=1,square=True,cbar=False,ax=ax[m[0]][m[1]])\n",
    "    ax[m[0]][m[1]].set_title(m[2])\n",
    "\n",
    "fig.delaxes(ax[1][2]); fig.delaxes(ax[2][1])\n",
    "fig.suptitle(\"tracer sub-matrices\")\n",
    "fig.tight_layout(rect=[0, 0.03, 1, 0.95]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Traceback\n",
    "The traceback function now permits movement between matrices."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "WARNING: 100 iterations exceeded; <10 alignments found: to continue searching, increase max_i.\n",
      "1. \n",
      " ----------------KKSGGLLLGK- \n",
      " KKASKASPKKPHPHHSPRNGGGGLGKK \n",
      "\n"
     ]
    }
   ],
   "source": [
    "alignment = seq.traceback_iterator(tracer,seq_1,seq_2,mat=M,local=True,affine=True)\n",
    "for i,e in enumerate(alignment):\n",
    "    print(str(i+1)+\".\",\"\\n\",e[0],\"\\n\",e[1],\"\\n\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Example Alignments\n",
    "\n",
    "Aligning more sequences.\n",
    "\n",
    "### 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "seq_2 = \"KKASKASPKKPHPHHSPRNGGGGLGKK\"\n",
    "seq_1 = \"KKSGGLLLGK\"\n",
    "\n",
    "S = pd.read_csv('blosum62.csv',index_col=0)\n",
    "g=-6\n",
    "g_open = -6; g_ext = -6"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "WARNING: 100 iterations exceeded; <10 alignments found: to continue searching, increase max_i.\n",
      "Max. score = 22.0\n",
      "1. \n",
      " ----------------KKSGGLLLGK- \n",
      " KKASKASPKKPHPHHSPRNGGGGLGKK \n",
      "\n"
     ]
    }
   ],
   "source": [
    "_,_,_,_,_ = seq.affine_local_dp(seq_1,seq_2,S,g_open,g_ext,disp_alignments = True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For equal opening and extension penalties, the affine case reduces to the linear case:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "WARNING: 225 iterations exceeded; <15 alignments found: to continue searching, increase max_i.\n",
      "Max. score = 22.0\n",
      "1. \n",
      " ----------------KKSGGLLLGK- \n",
      " KKASKASPKKPHPHHSPRNGGGGLGKK \n",
      "\n"
     ]
    }
   ],
   "source": [
    "_,_,_ = seq.local_dp(seq_1,seq_2,S,g,disp_alignments = True,N=15)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "seq_1 = \"ATCCCCGGGGGGGCT\"\n",
    "seq_2 = \"AT\"\n",
    "\n",
    "S = pd.read_csv('dna_transition_transversion.csv',index_col=0)\n",
    "g_open = -5; g_ext = -2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "WARNING: 100 iterations exceeded; <10 alignments found: to continue searching, increase max_i.\n",
      "Max. score = 4.0\n",
      "1. \n",
      " ATCCCCGGGGGGGCT \n",
      " AT------------- \n",
      "\n"
     ]
    }
   ],
   "source": [
    "_,_,_,_,_ = seq.affine_local_dp(seq_1,seq_2,S,g_open,g_ext,disp_alignments = True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "seq_1 = \"KRRLHMNGDPRS\"\n",
    "seq_2 = \"PQRS\"\n",
    "\n",
    "S = pd.read_csv('PAM250.csv',index_col=0)\n",
    "g_open = -2; g_ext = -1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "WARNING: 100 iterations exceeded; <10 alignments found: to continue searching, increase max_i.\n",
      "Max. score = 12.0\n",
      "1. \n",
      " KRRLHMNGDP-RS \n",
      " ---------PQRS \n",
      "\n"
     ]
    }
   ],
   "source": [
    "_,_,_,_,_ = seq.affine_local_dp(seq_1,seq_2,S,g_open,g_ext,disp_alignments = True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 4."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "seq_2 = \"AAATTGCGCAAAATCTACTACGCGCGCGA\"\n",
    "seq_1 = \"AAATTGCGCAAAATCTTACGCGCCCTACGCGA\"\n",
    "\n",
    "S = pd.read_csv('dna_transition_transversion.csv',index_col=0)\n",
    "g_open = -10; g_ext=-1 #compare with example 2. in notebook 5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "WARNING: 100 iterations exceeded; <10 alignments found: to continue searching, increase max_i.\n",
      "Max. score = 35.0\n",
      "1. \n",
      " AAATTGCGCAAAATCT--TACGCGCCCTACGCGA \n",
      " AAATTGCGCAAAATCTACTACGCGCGCGA----- \n",
      "\n"
     ]
    }
   ],
   "source": [
    "_,_,_,_,_ = seq.affine_local_dp(seq_1,seq_2,S,g_open,g_ext,disp_alignments = True)"
   ]
  }
 ],
 "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}