2_global_dp.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Global Sequence Alignments\n",
    "We can perform optimal global sequence alignments with a dynamic programming approach."
   ]
  },
  {
   "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 = \"AAAAATGCGCC\"\n",
    "seq_1 = \"AAAATGCCA\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Initialization Step\n",
    "\n",
    "First, we need to initialize an $(m+1) \\times (n+1)$ matrix and assign gap penalties. We will call this matrix the scoring matrix $(M)$. <br>\n",
    "*Note:* We will start off with linear gap penalties and modify the approach for affine penalties in another notebook."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = len(seq_1); n = len(seq_2)\n",
    "M = np.zeros([m+1,n+1])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's specify a gap penalty $g$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "g = -6"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$M$ is initialized as: <br>\n",
    "$$ M_{(0,0)} = 0 $$ <br>\n",
    "$$ M_{(0,j)} = gj \\ \\forall j>0 $$ <br>\n",
    "$$ M_{(i,0)} = gi \\ \\forall i>0 $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "#initialize the matrix\n",
    "M[0,:] = g*np.arange(0,n+1)\n",
    "M[:,0] = g*np.arange(0,m+1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAASYAAAD7CAYAAADHEzmfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO2deVhUZfv4PwMMq4oICijuWKjkhhvupuKG+74XWq6VVm9qi5a5a5lli5bmnmGmZq/lmqaApmZqmrgjJggoyCLrzPn9wTdMZRmGGd5zrt/9ua6uiznnzOe+eWxunnPmnOfWKYqiIAiCoCJs/tcJCIIgPI4UJkEQVIcUJkEQVIcUJkEQVIcUJkEQVIcUJkEQVIcUJkEQVIddaQc0xj5lMZeN1yUAjtzwtZgToG2NKwAsuNDDot6Z9XYDMCB8kkW921p9BkD9ne9azHm+T66r+ldLLOYEiBr3HwCemr/Mot5Lb04DoME0y3rPLsv1Bg77wGLOiG9eA6Bz+/kWcwLsP/wmAN1rvW5R70/XlgLQxWaQRb37jFsL3CczJkEQVIcUJkEQVIcUJkEQVIcUJkEQVEepX/wuikMRsGwVZGXD07Vg7nQo41Jy763rsPkzG9LTwMYGRr1ipEYd832KonD0k3jcqtnj37c8ADmZRo6tSiDhSiYoCh51HGn5ogd2DqbXf0VR+Hv1XzhWccGje/XcbUaFmI2RpEUmAlDmGQ+8hvii0+mK9AX7NCDEtxUKkG7IZsG5nzifdPuRY6b7d6WaizuTj282Oc8x9Rozsm4jFEUhKiWJGUf2cjfjATY6He+06Eg7nxrY6WxYde4Emy6eMdk7MqAhw5o0REHhZuJ93t69j3sP0gEY3qQBgxr542hnx5+xcbz5331kGwwmuzv612beiK60mpn7ZYFOB1OD29K2Xk0Uo0JUQiLvhx4gMS3dZCdAu6a+zJrYjc5jV+RtGzsgkM6BT2MwKkRev8Oir/aRlW1aruMndaJ9Bz+SkzMAuBV9l7nv7XjkmImTO1PFx423ZxZ8AflxajztxcTZ/XAp64jBaOSTt77jyp9/AzBk0rN06tcUW1sbDu78nU3L95rm9K/GlI9DcHF1xmgw8tGEVVz+/RoAbfq3YNjM/ugd7IiLSmDRmE9IuZdqkldVM6Z7SfDWQlj+Pvy0EXwqwwcrS+7NzIBlb9rQbZCR2Z8ZCR5u5KuF5v/qSdFZ7Jl1mxthjw7y2e8SUYwKfZb50HtZVQxZCue2JZnszbidxo3Fp7l/4s6j8cJjyIx5gO/7LfF9rwUPIhNJPhlXpK9GGXder9+FFyM2MuDQF6yM/JXlzYY8ckzXyvUJ9mlgco4A/u6evPBMM/r/sImg79dy434SrwW0AWCEX0NquroRtO1reu3cQIh/AA0repnkre9ViZAWAQxZv4XgLzcQdS+Rqe1bARD0tC+jmjbiuc3b6LFqPY52djzfvInJOVfzKM9rvdth869i3q+5P/V8KjFk6SYGLNlAdMJ9Xu/TrhgjAT5e5XlpRDt0Ng+9jev60CXQj+dmbmTkG+twcbJnUNfGJjvr+1dh7ns7mDBuNRPGrX6iKLXvWJdOQf7FytPBUc+8dS/y3apfmNJrGd98sp83lo0AoFkHP9p2b8hLvT9iQrelNGxZm7Y9GhbtdLJn4Z63CV2yk4kBb7Bx7nfM2PgyAE8F1GLKJ2OZM3ApLzZ4jVuXbxMyb5jJ+RY5Y7p69Sp79uwhNjYWGxsbKlWqRNu2bXnmmWdMDmIqYSfA3w9q+OS+HtYH+o6FWdNy/7qZy4XfoaI3NGie+7pRIHh4Gc32XfzpPnU6lcOl4qN/WT3rO1Gmoh06Gx06oEIte5JuZpnsvXfgFuXbeqN3d3xku2IEY5YBJduIoigoBgWdvujCmmUwMOuPH0jIzC2g55Nu4+FYBr3OlmzFQK0yHoTUac3nkYdpXcn0Wy7+vHuHDqFfkaMYcbC1xdOlDNEp9wHoWqMOmy+ewaAoJGdlsuvaRfr51uNMfGyR3vOxcQR9sZYcoxF7W1s8y5bhVlIyAH3967Lm+O/cz8gEYNbPB7C3Me2Pi6PejgUju7F052EWjuqet/1K7F0+3HUkb9Z1PvoOQ1sX/YH8Bwd7O96d3IPlGw/z3pSHt5bY2thgb2+Lg70dBqMRe70dmSbOlvR6W3x9vRg0tCWVK5fn9t+JfL5iP3FxueNQrbo7g4e2ZOO6ozRtVtPkXJu0fYqYm3c5cegiAMf2nyf21j0AAoP8+eWH02Sm5/6/uve7EzzbtwlHdhc+0w0IakjM1Tv89tNpACJ+OEns9dw/mJ1GtuPnNQe5ExUPwPp3t1LOvYzJ+RZamDZt2kRoaChdu3bNK0Tx8fG888479O7dm5CQEJMDmUJsHHhXevjasyKkpulIe6CU6HQu9pYOVzdY+6GO6Gs6nF1g4DjzC1PLFysCcPvso4WpSiPnvJ9T47K5sOs+rSZWNNlbedTTAKRdSHxku1sbb5JP3OHiq0fBqFCmfgXKNSraezs9idvpD2dsb/h35ZfYSLIVA8629iwI6M9bv+/Av3xlk3P8hxzFSFB1Xxa17UqWwcCHp44C4O1SlpjUlLzjYtNSqFvB9DHIMRrp/FRt5vXoQpbBwPJfc09Vari74R4Ty1dD+lGprAunov9m8cEjJjnfGdyZreHnuHQ74ZHtZ6Ni8n4u6+TAhKAWhIafNTnX6eO6sOPAWa7834fvH06ev8lv56LYvuIFcnIMRN1OZMd+005n3d3LcPr0DVav+oVb0fcYPLQFc+YPZMK4NTg66Zn+Zm+WLPyRp542bRb6D1VqViQxPoWpCwdRs25l0pLTWb3wvwBU9C7PmfAreccmxN7Hw8u1SKfPU97ci03i1a8mUrtBdVKT0vhy+sbcfXW8uXbuJu9tfwPPGhW5fu4mX7y6zuR8C/2Ts379erZs2cKkSZMYNGgQgwYNYtKkSXzzzTeEhoaaHMRUjAXUChP/MBaIIQfOnYB2PRTeWWHk2T5Glr9jQ7bpk5likXA1k5/euk3dHq5UbVbyC2RxO69hW9Yev+VtefqDNhjSskn4Ocrk9zvZ6vmw6SCquVRg1ukfAJjTuDebrx3nSkrRp4QFsTfqCo03fsqy38PZ0G0QOnjkVOkfDMbirUW4/9JVWnz0BZ8ciWDN0P7oADsbG1rVrM4r2//LgDWbcXV0ZFr71kW6hrRugMFgZMdv5ws8xsfdlbVTBvP79dtsOWpaAenfpSEGg5EfD/35xL7gDv5UruhKr4lfEDxxJTHx93l5VAeTvLGx93lreii3onNnM6FbjuNd2Q0vL1def6MnO74/yY3r8UVYnsTOzpamHfz46ZvjvNJnOT+sC2POmrHo7W2xsXny38xowr+Znd6O5j0as3vVPiY3n8GOFT8x779vore3w1ZvR8vgAD6asIqJTd4g8U4S01aNNz3fwn8ZO3Jycp7YnpGRgV6vNzmIqXh7wtm/Hr6+kwCuZRWcnYrv2rFOx5ljuQOemgxeVaGWX+6+xq1g3UcQHwuVqxXtOr35HjdPpAFQrZkLjYdXKPDYa0dSOLYqgZYveFCrXdlCvXe2XyXldO5f8bKNPfDsVzvf45JPxeM94ils7GzAzobyrb1JPhmHR7fqTxw7xa8jHb1yZ16/xEayLep3Pm0xjKspCTwftpZMYw6ejuUIcK9OzTIejK4diKveiTJ6Bz5vOYKJxzblm8OrTVrTuXru6d6VxATW//UHJ+/kXjgNvXSO+a274OrgyO3UFCo5P5yyezqXJSYtJV8nwMvtAulUp1auN+Eem0+d4dSt3Ivz3505z3vdOuHq5Ehcahr7Iq+QlpX71+SHPy8yuU2LfJ2TugXSwT/XmZ1jxMnejtDXR6C3tcVBn/vz5FU7iE9Oo5mvD0tG9+TrgydZd+hUgXkCvDCwFW0Ccv+NcnIMODjoWbdgFHq73NO2dQtG8dri7+nQzJc9YX/xICMbgB0HzvLa850K9I4JaUdgq9xvYmJikgg7Esn+vQ8Lnk4HOQYj/g2q4lPVnQGDmlOunCMuLg7MWzSYt6bnP0kYNbUrLTrXA6BceRduXY0j8sxNIPdUbuqCQXhVdSfudhJulR7+v+rh6UpCTP7XRse8N4TAXk0BKOtehuiLf3Pxt9zZVsQPJ3n1y4l41/Lk7u17XP/zJol3cj17vv6FJQdmFzq+/6bQwjRhwgT69u1LYGAgFSvmTsfj4+M5duwY06ZNMzmIqbRuBos/gxu3cq8zffsDPFv0H8V86TtGoe+Y3Kp//x6884INNy5DjTpw6RzoABOvydJ4eIVCi9E/3AhP5bevEgia7Y2Hr2ORx3v2q11gMfo3TtXLknwijjJ1K6DkGEk5nYBTrfyn2isu/sKKi78A4Kp3IrT9i+yI/oPPIw/nHXMnI5mOex4+YtG3aiO6VK5X6LdyH/4exoe/hwHQzLMKK57tRffv15GYmU7f2vWITEwgKTODfVGXGfy0P/tvXsFFb0/v2n68ebTgb3g+/jWCj3+NAKBp1Sos69uDPqs3kpieQe/6flyOv0tSegZ7Ll6me92nCP3jHJk5Bjo/VZtzMXfydX72cwSf/RzxxPbKbuX4fvooBi/NLb4Na3izLKQX09fvJuxi0TPQL78L58vvwp/Y7uVRjk1LxjBm5gYAIq/H0aFZHX4+cgGDUaFD8zqcvxzzxPv+Yd2aX1m35lcAatSsyLJPRvLn2WhiY+/Tu28Trl2NIyE+haEDPsl7T1C3Z2jX3q/Qb+U2fLSHDR/tAcDNoywr9/4HX/8qXPnzb/yb1UJRIDb6Hsf2nWf4K1346ZvjGA0GOg9oyv5tJ/PPdfa3rJv9ba7TszxfnV9GnSa1uPz7NZ5pWxdFUYi5HseRbcd4ccloNs/bRsq9VNr0b0HkiatFjPBDCi1MvXr1onnz5kRERBAXF4eiKDRt2pSXXnoJT09Pk4OYirsbzJsBU2dBdjZUrQIL3yy517UCTHnXyKZPbMjMADs9TJplRG9fcve/ObXxHgoQ9unDqbannyMtx5t+jSU/vIY9RczGSC7NjEBno8OlnhsVezw5W3qcITWb4u3sSmfvunT2rpu3PSRsHfezi/eV+L85cedvVvxxjG97DiFHUYh7kMqL+3K/Odrw1x9UK1een/s/h97Gls0Xz3A89pZJ3pPRf/N5+G9sGDkIg9FIXEoak77LPfXcdOoMro6ObA8ZgY1Ox4XYOBb89KvZvwPkzqx06HgluA2vBOd+q/j33WSmfb2rRN51O4/z8qgObF76HNnZBi7fjGfpmgMmvffG9XhWLN/L+wsGY2OrIyE+hflzdpYoH4DEhBTmjF/L5DkDcHSyJzsrh7mT1pKdlcPxgxeo4efF8u0vY2dvy7F959n/ff6F6RHnnSTe7beYlz4dh6OLA9mZObw3YCnZmdkc+/EUHj7ufHDoPWxsbLgTFc+H4z43OV9daTcjkId45SFekId4QR7ilYd4BUHQFFKYBEFQHVKYBEFQHVKYBEFQHVKYBEFQHVKYBEFQHaV+u4AgCEJRyIxJEATVUeoLxVnyZsh/boS05E2b8PDGzf+fvVrKVWteLeVaGt5891k0kiAIggWQwiQIguqQwiQIguqQwiQIgupQXZcUsHxHE7BO9xVrdXTRkldLuWrNq6VcLe1V3YzJ0h1NwDrdV6zV0UVLXi3lqjWvlnK1hld1hSm/jibj3zK/cQDk333lx/1QkltLreHUmldLuWrNq6VcreFVXWH6d0eT96fY8OEMG4ym9zTM31lg9xV1ObXm1VKuWvNqKVdreAu9xnT79u3CdlO5cvHb/hTFPx1NXl+sUMtP4XQ4LH/HhkXrzV8K1xrdV6zV0UVLXi3lqjWvlnK1hrfQwjR+/Hhu3LhBpUqVePyROp1Ox4EDpq1jXBSW7GiSH5bsvmJNp9a8WspVa14t5WoNb6H17JtvvqFmzZosXryYgwcPPvKfpYoS5HY0mf25kdmfG3n7EyMJsXDjcu6+4nY0yY/WzeDMhdzuK1Cy7ivWdGrNq6VctebVUq7W8Ba5usDZs2fZunUr77//vvlR/oUpz8pdOgdbv3zY0WTYRCN18mnVXpxn5Q4fy/0q89/dV8qXy/9YU58NKo5Ta14t5ao1r5ZyLQ1vfpT6sifyEK82vFrKVWteLeVaGt5891k0kiAIggWQwiQIguqQwiQIguqQwiQIguqQwiQIguqQZgSCIKgOmTEJgqA6Sn09pgUXeljMNbPebkB7921owaulXLXm1VKupeHNd59FIwmCIFgAKUyCIKgOKUyCIKgOKUyCIKgOKUyCIKgO1XRJURSFo5/E41bNHv++5QHIyTRybFUCCVcyQVHwqONIyxc9sHMofj2VLhYyBlrzailXS3tVMWNKis5iz6zb3AhLfWT72e8SUYwKfZb50HtZVQxZCue2JRXbL10sZAy05tVSrtbwFlmY9u/fz4YNG7h58+Yj27/99lvzoz7GxZ/uU6dTOWq0LvPIds/6TjQY6IbORoeNrY4KtexJjc8utl+6WMgYaM2rpVyt4S20MC1dupSNGzdy48YNhg4dys6dO/P2bdmyxbyI+dDyxYrU7lD2ie1VGjnjWiW3A0FqXDYXdt2nRqsyTxxXFNLFQsZAa14t5WoNb6HXmA4fPsz27duxs7Nj1KhRhISEYG9vT/fu3Z9oTmBNEq5m8svCWOr2cKVqs+KftEoXCxkDrXm1lKs1vIW+TVEUdLrc7iU1atRg5cqVzJs3j+PHj+dtN4fTm++xc1o0O6dFc3rzvUKPvXYkhb3v3iZgVAUaDHQzK563J8TfffjaUl0sLO3UmldLuWrNq6VcreEttDB169aNUaNGcfbsWQDq1KnD8uXLmTp16hPXnIpD4+EV6LOsKn2WVaXx8AoFHncjPJXfvkogaLY3tdo9eapnKtLFQsZAa14t5WoNb6GnclOmTCEgIAAXl4enTwEBAXz//fesWbPG/KgmcmrjPRQg7NP4vG2efo60HF+xWB53N5g3A6bOerSDQ0mwhlNrXi3lqjWvlnK1hrfU12OS1QW04dVSrlrzainX0vDmu8+ikQRBECyAFCZBEFSHFCZBEFSHFCZBEFSHFCZBEFSHdEkRBEF1yIxJEATVUerrMQ0In2Qx17ZWnwHau29DC14t5ao1r5ZyLQ1vvvssGkkQBMECSGESBEF1SGESBEF1SGESBEF1qKoZwd+r/8Kxigse3avnbjMqxGyMJC0yEYAyz3jgNcTXrLWgZLF4GQOtebWUq6W9qpgxZdxO48bi09w/ceeR7UnhMWTGPMD3/Zb4vteCB5GJJJ+MK7ZfFouXMdCaV0u5WsNbZGG6ceMGd+7kFoytW7cyd+5cdu/ebX7EfLh34Bbl23rj2szzke2KEYxZBpRsI8YcI4pBQacvfi2VxeJlDLTm1VKu1vAWeiq3du1aNmzYgNFopGXLlsTExNClSxe2bdvG9evXmTx5snlRH6PyqKcBSLuQ+Mh2tzbeJJ+4w8VXj4JRoUz9CpRrVLxF4qCwhdIVs6ea1nBqzaulXLXm1VKu1vAWWpi2bdvG7t27SUhIIDg4mGPHjuHg4MCgQYMYOHCgxQpTQcTtvIZtWXv8lrdFyTJy85MzJPwchUe36sXyyGLxMgZa82opV2t4C32b0WjE3t6eKlWqEBISgoODQ94+g8FgXkTgzvarXJl1nCuzjnNn+9UCj0s+FY9bW29s7GywdbajfGtv0i4mFnh8Qchi8TIGWvNqKVdreAstTEFBQYwcORKDwcBLL70EwMWLFxk+fDjdu3c3LyLg2a82vnNa4DunBZ79ahd4nFP1siSfyL3YreQYSTmdgFMt12LHk8XiZQy05tVSrtbwFrm6wIkTJ2jWrFne62vXrhEdHU379u3NCljYs3K3vrrwyO0COanZxGyMJD0qBZ2NDpd6bngPqYPOLreeFudZucPHcr/K/PdC6eXL5X+sqc8GFcepNa+WctWaV0u5loY3P0p92RN5iFcbXi3lqjWvlnItDW+++ywaSRAEwQJIYRIEQXVIYRIEQXVIYRIEQXVIYRIEQXVIYRIEQXVIlxRBEFSHzJgEQVAdpb5QXP2d71rMdb5PrktrN5RpwaulXLXm1VKupeHNd59FIwmCIFgAKUyCIKgOKUyCIKgOKUyCIKgOVXRJCfZpQIhvKxQg3ZDNgnM/cT7p9iPHTPfvSjUXdyYf32xWDOliIWOgNa+WcrW0938+Y6pRxp3X63fhxYiNDDj0BSsjf2V5syGPHNO1cn2CfRqYHUO6WMgYaM2rpVyt4S1WYVq4cKH5kQogy2Bg1h8/kJCZCsD5pNt4OJZBr7MFoFYZD0LqtObzyMNmx5AuFjIGWvNqKVdreAs8lZs5c+YT2w4ePMj9+/cBWLBggXkRH+N2ehK305PyXr/h35VfYiPJVgw429qzIKA/b/2+A//ylc2OIV0sZAy05tVSrtbwFjhjKl++PIcOHcLPz4/mzZvTvHlznJ2d8362NE62ej5sOohqLhWYdfoHAOY07s3ma8e5klL8Jpf/RrpYyBhozaulXK3hLfBt06dP58MPP2T37t1UrlyZfv364erqSr9+/ejXr5950f6PKX4d2dZhAts6TGCKX0e8nVzZ1HYsBkXh+bC1pORk4OlYjgD36oyuHZh3XIB7NT5vOaLY8aSLhYyB1rxaytUa3kLrWWBgICtXrmTz5s0sWrSoRC2b/s2Ki78w4NAXDDj0BRuuHmNt6+fYF/MX/zn1HZnGHADuZCTTcc8HecetuPgLp+7eZOKxTcWOJ10sZAy05tVSrtbwFnm7QPny5Vm+fDlbt24lMjLS/EgFMKRmU7ydXensXZfO3nXztoeEreN+drpFYri7wbwZMHXWox0c1ObUmldLuWrNq6VcreEt9WVP5CFebXi1lKvWvFrKtTS8+e6zaCRBEAQLIIVJEATVIYVJEATVIYVJEATVIYVJEATVIYVJEATVIV1SBEFQHTJjEgRBdZT6QnHVv1piMVfUuP8A2ruhTAteLeWqNa+Wci0Nb777LBpJEATBAkhhEgRBdUhhEgRBdUhhEgRBdaiiS8qYeo0ZWbcRiqIQlZLEjCN7uZvxABudjndadKSdTw3sdDasOneCTRfPmBVDuljIGGjNq6VcLe39n8+Y/N09eeGZZvT/YRNB36/lxv0kXgtoA8AIv4bUdHUjaNvX9Nq5gRD/ABpW9Cp2DOliIWOgNa+WcrWGt9DCdPbs2byfIyIiWLhwIUuXLuXMGfNmLfnx5907dAj9ipTsLBxsbfF0KUNiZu4CcV1r1CH00jkMikJyVia7rl2kn2+9YseQLhYyBlrzailXa3gLLUyzZ88GYNOmTcyfPx8vLy88PDyYNWsWGzduNC9iPuQoRoKq+3Js2ARaePmw9dI5ALxdyhKTmpJ3XGxaCt4uZYvtL7iDg/k5W8OpNa+WctWaV0u5WsNr0jWm0NBQ1q9fj5ubGwADBw5k4MCBjBw50ryo+bA36gp7o64w9OkGbOg2iHahX2Kj0z1xnMFY/BIsXSxkDLTm1VKu1vAW+racnByMRiPu7u44Ozvnbbe3t8emBL/Jq01as7vfGHb3G8PHHXrS1LNK3r7QS+eoUqYcrg6O3E5NoZJzmbx9ns5liUlLyU9ZKNLFQsZAa14t5WoNb6HVxc3Njfbt23PlypW807qIiAiGDh1Kt27dzIsIfPh7GD22r6PH9nVs+OsPVjzbCzeH3N+gb+16RCYmkJSZwb6oywx+2h9bnY5y9g70ru3H3qjLxY4nXSxkDLTm1VKu1vCatLrAtWvXSE5OplGjRpw6dYqUlBQ6dOhgVsD8npUbWbcRo+s2IkdRiHuQyjth+4lOvY+tTsdbLTrQtkoN9Da2bL54hlXnTuS9rzjPyh0+lvtV5r87OJQvl/+xpj4bVByn1rxaylVrXi3lWhre/Cj1ZU/kIV5teLWUq9a8Wsq1NLz57rNoJEEQBAsghUkQBNUhhUkQBNUhhUkQBNUhhUkQBNUhzQgEQVAdMmMSBEF1lPp6TE/NX2Yx16U3pwHau29DC14t5ao1r5ZyLQ1vvvssGkkQBMECSGESBEF1SGESBEF1SGESBEF1SGESBEF1qKJLysiAhgxr0hAFhZuJ93l79z7uPchd93t4kwYMauSPo50df8bG8eZ/95FtMBQ7hnSxkDHQmldLuVra+z+fMdX3qkRIiwCGrN9C8JcbiLqXyNT2rQAIetqXUU0b8dzmbfRYtR5HOzueb96k2DGki4WMgda8WsrVGt4iC9ORI0dITk4GYMeOHcyZM4dt27aZH/ExzsfGEfTFWlIzs7C3tcWzbBmSHmQA0Ne/LmuO/879jEwUYNbPB9h57kKxY0gXCxkDrXm1lKs1vIUWpnnz5rFy5UoyMzP56KOP2LVrF76+vuzbt4+5c+eaFzEfcoxGOj9VmyMvvUCzaj5sO3segBrubri7OPHVkH78MG4kL7dtSXJmZrH90sVCxkBrXi3lag1vodeYwsLC2LVrF7a2thw6dIjQ0FDs7e0ZMmQIwcHB5kUsgP2XrrL/0lUGN/JnzdD+dP58DXY2NrSqWZ2JW38gKyeHRb26Mq19a+bvP1wst3SxkDHQmldLuVrDW+jbHB0duXs3t/WBu7s7Dx7klr/09HTs7My/bv5yu0B2jh3BzrEj+KBPdwJ8Kuft++7MeSq7lsXVyZG41DT2RV4hLSuLbKORH/68SOMq3sWOJ10sZAy05tVSrtbwFlqYpkyZwsCBA1m0aBG1atVi1KhRzJ8/n8GDB/P888+bFxH4+NcI+qzeRJ/Vm/jm97Ms69sDNydHAHrX9+Ny/F2S0jPYc/Ey3es+hYOdLQCdn6rNuZg7xY4nXSxkDLTm1VKu1vAWuexJdHQ0+/fvJyoqCoPBgIeHBx07dqRBgwZmBczvId5hTRowIqAhBqORuJQ03ttzkFv3k7HR6ZjUugU96j2FjU7Hhdg43vnpAGlZWUDxHuKVLhYyBmrxainX0vDmR6mvxySrC2jDq6VctebVUq6l4c13n0UjCYIgWAApTIIgqA4pTIIgqA4pTIIgqA4pTIIgqA7pkiIIguqQGZMgCKqj1NdjajDNcvcxnV0m9zFZy6ulXLXm1aFpsMQAAA+ISURBVFKupeHNd59FIwmCIFgAKUyCIKgOKUyCIKgOKUyCIKgOVTQj+IeO/rWZN6IrrWZ+BoBOB1OD29K2Xk0Uo0JUQiLvhx4gMS292G5ZLF7GQGteLeVqaa9qZkzVPMrzWu922Oh0edv6Nfennk8lhizdxIAlG4hOuM/rfdoV2y2LxcsYaM2rpVyt4S20MM2dO5f79++bbzcRR70dC0Z2Y+nOR5fMvRJ7lw93Hclr13Q++g7eboUs8FIAsli8jIHWvFrK1RreQgvTjh07GDx4MHv37jXPbiLvDO7M1vBzXLqd8Mj2s1Ex/HUrDoCyTg5MCGrB3jMF3/tQELJYvIyB1rxaytUa3kILk4+PD59++inr169n0KBB7N69m4yMDPMiFcCQ1g0wGIzs+O18wXm4u7J2ymB+v36bLUfPFDuGLBYvY6A1r5ZytYa30IvfOp0OX19fNm7cSHh4ON9++y3z5s2jRo0aeHl58cEHH5gVdFK3QDr41wIgO8eIk70doa+PQG9ri4M+9+fJq3YQn5xGM18flozuydcHT7Lu0Cmz4nl7wtm/Hr621GLxlnZqzaulXLXm1VKu1vAWWs/+/Xxvq1atWL58OYcOHWLmzJl07tzZvIjAZz9HMHjpJgYv3cSIj76h/+INDF66icmrdpCZncPgpZuIT06jYQ1vloX04q3NP5tdlEAWi7eWV0u5as2rpVyt4S10xjRixIgntun1evz9/fH39zc/qolM6haIDh2vBLfhleA2APx9N5lpX+8qlsfdDebNgKmzHl0ovSRYw6k1r5Zy1ZpXS7law1vqy57IQ7za8GopV615tZRraXjz3WfRSIIgCBZACpMgCKpDCpMgCKpDCpMgCKpDCpMgCKpDCpMgCKpDuqQIgqA6ZMYkCILqKPWF4gKHmfd8XX5EfPMaoL0byrTg1VKuWvNqKdfS8Oa7z6KRBEEQLIAUJkEQVIcUJkEQVIcUJkEQVIequqS0a+rLrInd6Dx2Rd62sQMC6Rz4NAajQuT1Oyz6ah9Z2YZiu6WLhYyB1rxaytXSXtXMmHy8yvPSiHbobB52SWlc14cugX48N3MjI99Yh4uTPYO6Ni62W7pYyBhozaulXK3hLbIwRUREcPr0aQDWrFnDhAkTWLFiBVlZWeZHfQwHezvendyD5Rsf7ZJia2ODvb0tDvZ22NnaYK+3I9OM2ZJ0sZAx0JpXS7law1voqdzixYs5efIkOTk5+Pj4oNPpGDZsGAcPHmTOnDnMnTvXvKiPMX1cF3YcOMuVqPhHtp88f5PfzkWxfcUL5OQYiLqdyI79xW9GUHAHB8XsqaY1nFrzailXrXm1lKs1vIUWpiNHjrBz506ysrJo3749R48eRa/X065dO/r06VP8aPnQv0tDDAYjPx76Ey+PR3vGBXfwp3JFV3pN/ILsHCNvT+jKy6M68OHag8WKIV0sZAy05tVSrtbwFlqYFEUhJSWFBw8ekJGRQWpqKm5ubmRkZJCdnW1eROCFga1oE1AbgJwcAw4OetYtGIXeLve0bd2CUby2+Hs6NPNlT9hfPMjIjbXjwFlee75TseNJFwsZA615tZSrNbyF1rMXXniBoKAg+vTpw3/+8x9CQkJYvHgxw4cPZ8CAAeZFBL78LpwxMzcwZuYGxr6zmZFvrGPMzA28uuh7MrNyGDNzAwmJaURej6NDszrY/t8F8Q7N63D+ckyx40kXCxkDrXm1lKs1vEWuLpCRkYHBYMDFxYXIyEiOHj2Kn58frVubF7WwZ+W8PMqxackYOj3/CQD2elteHtWBZv7VyM42cPlmPEvXHCAtPffCe3GelTt8LPerzH93cChfQLdxU58NKo5Ta14t5ao1r5ZyLQ1vfpT6sifyEK82vFrKVWteLeVaGt5891k0kiAIggWQwiQIguqQwiQIguqQwiQIguqQwiQIguqQwiQIguqQLimCIKgOmTEJgqA6Sn2huM7t51vMtf/wm4D2bijTgldLuWrNq6VcS8Ob7z6LRhIEQbAAUpgEQVAdUpgEQVAdUpgEQVAdquiSMn5SJ9p38CM5OQOAW9F3mfvejkeOmTi5M1V83Hh75lazYkgXCxkDrXm1lKulvaqYMdX3r8Lc93YwYdxqJoxb/URRat+xLp2C/M32SxcLGQOtebWUqzW8RRam/fv3M2PGDMaOHcuECROYP39+XtcUS6DX2+Lr68WgoS1ZuXoss+f0p1Klh6tLVavuzuChLdm47qjZMaSLhYyB1rxaytUa3kIL08qVK9m2bRsNGjRAp9PRqFEjKleuzNtvv01oaKh5ER/D3b0Mp0/fYPWqXxg/djV/XfibOfMHAuDopGf6m71ZsvBHHjzINDtGwR0czM/bGk6tebWUq9a8WsrVGt5CC9Pu3bv57LPPGD58OJ9++inh4eE899xzbNmyha+//tq8iI8RG3uft6aHciv6HgChW47jXdkNLy9XXn+jJzu+P8mN6/FFWApHuljIGGjNq6VcreEt9OJ3ZmYm6enpODs7k5GRQVJSEgDOzs7YlOA3GRPSjsBWdQCIiUki7Egk+/f+mbdfp4McgxH/BlXxqerOgEHNKVfOERcXB+YtGsxb04s3W5MuFjIGWvNqKVdreAutLv3792fYsGEsWbKE0aNH079/f/7++28GDBhAcHCweRGBdWt+zbvQvW7Nr0x+uQteXq4A9O7bhGtX40iIT2HogE/yjlu75lfOnY0udlEC6WJhLa+WctWaV0u5WsNb5OoCERERXLhwgXr16hEYGEhaWhq3bt3i6aefNitgfs/KdepSn6HDW2FjqyMhPoUPFv2XuLjkR44J6vYM7dr7PXK7QHGelZMuFjIGavFqKdfS8OZHqS97Ig/xasOrpVy15tVSrqXhzXefRSMJgiBYAClMgiCoDilMgiCoDilMgiCoDilMgiCoDmlGIAiC6pAZkyAIqkMKkyAIqkMKkyAIqkMKkyAIqkMKkyAIqkMKkyAIqkMKkyAIqkMKkyAIqkMKkyAIqkMKkyAIqkN1hWnXrl306NGDLl26sGnTJot5U1NTCQ4O5tatWxZzrlixgp49e9KzZ08WL15sMe/y5cvp0aMHPXv2tFjTh3+zaNEiZsyYYTHf6NGj6dmzJ3369KFPnz6cOXOmxM6DBw/Sv39/unXrxty5cy2QJWzdujUvxz59+hAQEMCcOXMs4t65c2fe/wuLFi2yiBNg1apVdO3alV69evH555+XyPX4ZyA8PJxevXoRFBTEsmXLLOYFmD59Ot9//735ySoqIjY2VunYsaOSmJiopKWlKb169VIuX75cYu8ff/yhBAcHK/Xr11eio6MtkKmihIWFKUOGDFEyMzOVrKwsZfTo0crevXtL7D1+/LgydOhQJTs7W0lPT1c6duyoXL161QIZ5xIeHq60aNFCmT59ukV8RqNRad26tZKdnW0Rn6Ioys2bN5U2bdooMTExSlZWljJs2DDl0KFDFvMriqJcunRJ6dKli3L37t0Sux48eKA0a9ZMuXv3rpKdna0MHDhQCQsLK7E3LCxMCQ4OVlJSUpScnBxl/Pjxyp49e8xyPf4ZSE9PV9q3b6/cvHlTyc7OVkJCQswa48e9sbGxyvjx45UGDRoo27ZtMytXRVEUVc2YwsPDadmyJeXLl8fZ2ZmuXbvy888/l9gbGhrK7NmzqVSpUtEHm0jFihWZMWMG9vb26PV6ateuze3bt0vsbd68OevXr8fOzo67d+9iMBhwdna2QMaQlJTEsmXLmDBhgkV8ANeuXUOn0/HCCy/Qu3dvNm7cWGLnvn376NGjB15eXuj1epYtW0bDhg0tkO1D3n33XaZNm0aFChVK7DIYDBiNRtLT08nJySEnJwcHB4cSey9cuECbNm0oU6YMtra2tG3blv3795vlevwzcPbsWapXr07VqlWxs7OjV69eZn3WHvfu2rWLTp060b17d7Py/IdC2zeVNnFxcVSsWDHvdaVKlTh79myJvfPmzSux43Hq1KmT9/ONGzfYvXs3W7ZssYhbr9fz8ccfs2bNGrp164anp6dFvLNmzWLatGnExMRYxAeQnJxMYGAg7777LhkZGYwePZqaNWvSurX5LTKioqLQ6/WMHTuW+Ph4OnbsyNSpUy2Wc3h4OBkZGSX+8PxDmTJleOWVV+jevTuOjo40b96cJk2alNhbv3595s+fz/jx43FycuLgwYMoZi4G8vhnIL/P2p07d0rsHTduHACnTp0yI8uHqGrGlN+g63S6/0EmpnP58mVCQkKYPn06NWrUsJj35ZdfJiIigpiYGIt0Pd66dSve3t4EBgZaILuHNG7cmMWLF+Ps7EyFChUYOHAghw8fLpHTYDAQERHBkiVLCA0N5dy5c2zfvt1CGcOWLVt4/vnnLea7ePEi27Zt45dffuHo0aPY2NiwevXqEnsDAwPp378/o0aNYty4cQQEBKDX6y2Qsfo/a6oqTJ6eniQkJOS9jouLs+jpl6U5deoUzz33HK+99hr9+vWziPPq1av89Vdu50AnJyeCgoKIjIwssXf37t2EhYXRp08fPv74Yw4ePMj8+SXvWHPy5EkiIiLyXiuKgp1dySbiHh4eBAYGUqFCBRwdHenUqZNFZs4AWVlZnDhxgmeffdYiPoCjR48SGBiIu7s79vb29O/fn99++63E3tTUVLp06cKuXbvYsGEDTk5OVK1a1QIZq/+zpqrC1KpVKyIiIrh37x7p6ens3buXdu3a/a/TypeYmBgmT57M0qVL6dmzp8W8t27d4u233yYrK4usrCwOHDhAQEBAib1ff/01P/74Izt37uTll1/m2Wef5c033yyxNyUlhcWLF5OZmUlqairbt2+nS5cuJXJ27NiRo0ePkpycjMFg4MiRI9SvX7/EuQJERkZSo0YNi123A/Dz8yM8PJwHDx6gKAoHDx7kmWeeKbH31q1bTJ48mZycHFJSUti6davFTj8bNmzI9evXiYqKwmAw8OOPP6rqs6aqa0yenp5MmzaN0aNHk52dzcCBA2nQoMH/Oq18Wb16NZmZmSxcuDBv29ChQxk2bFiJvO3bt+fMmTP07dsXW1tbgoKCLFr4LE3Hjh3z8jUajQwfPpzGjRuXyNmwYUPGjRvH8OHDyc7OpnXr1gwYMMAi+UZHR+Pl5WUR1z+0adOGCxcu0L9/f/R6Pc888wwvvvhiib1+fn4EBQXRu3dvDAYDzz33nEX+SAE4ODiwcOFCXnrpJTIzM2nfvj3dunWziNsSyNK6giCoDlWdygmCIIAUJkEQVIgUJkEQVIcUJkEQVIcUJkEQVIcUJkEQVIcUJkEQVIcUJkEQVMf/A9gcx3abL9yFAAAAAElFTkSuQmCC\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": [
    "As you can see, the initialized matrix has linearly decreasing scores along the first row and first column."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Fill Step\n",
    "To fill the matrix with scores, we need to iterate over the elements of the matrix and fill in values according to the following recursive relation:\n",
    "\n",
    "$$ M_{i,j} = \\max \\{ (M_{i-1,j-1} + S_{i,j}), \\ (M_{i-1,j} + g), \\ (M_{i,j-1} + g) \\} $$\n",
    "\n",
    "Where $M$ is the scoring matrix, $S$ is the substitution matrix and $g$ is the gap penalty."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Substitution Matrix\n",
    "Consider the following substition matrix for DNA sequences."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>A</th>\n",
       "      <th>G</th>\n",
       "      <th>C</th>\n",
       "      <th>T</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>A</th>\n",
       "      <td>2.0</td>\n",
       "      <td>-5.0</td>\n",
       "      <td>-7.0</td>\n",
       "      <td>-7.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>G</th>\n",
       "      <td>-5.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>-7.0</td>\n",
       "      <td>-7.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>C</th>\n",
       "      <td>-7.0</td>\n",
       "      <td>-7.0</td>\n",
       "      <td>2.0</td>\n",
       "      <td>-5.0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>T</th>\n",
       "      <td>-7.0</td>\n",
       "      <td>-7.0</td>\n",
       "      <td>-5.0</td>\n",
       "      <td>2.0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "     A    G    C    T\n",
       "A  2.0 -5.0 -7.0 -7.0\n",
       "G -5.0  2.0 -7.0 -7.0\n",
       "C -7.0 -7.0  2.0 -5.0\n",
       "T -7.0 -7.0 -5.0  2.0"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "S = pd.read_csv('dna_transition_transversion.csv',index_col=0)\n",
    "S"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "tracer = np.zeros([np.shape(M)[0],np.shape(M)[1],3])\n",
    "for i in range(1,m+1):\n",
    "    for j in range(1,n+1):\n",
    "        arr = np.array([M[i-1,j-1]+(S[seq_1[i-1]][seq_2[j-1]]),(M[i-1,j]+g),(M[i,j-1]+g)])\n",
    "        M[i,j] = np.max(arr)\n",
    "        idx = np.where(arr==np.max(arr))\n",
    "        tracer[i,j,idx] = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAASYAAAD7CAYAAADHEzmfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOydd1QVV/u2r9MoSu9FEQEVFUtERRSxIPYee4nRxBpNonmTGGOJRhNjYtdYE3uLLcYWRewBe0PsCAhI7yDtlO+PUZFIOcghH7+15lrLJeyZuWczs+eZPXv2PLdEo9FoEBEREalCSP9/V0BERETk34iBSUREpMohBiYREZEqhxiYREREqhxiYBIREalyiIFJRESkyiEGJhERkSqH/L/e4YnwBjrT6lL7HgDz7vbWmSbAbI+/AHg/aJJOdfe3/hWApkdn6VT3Vo/vAXDe+pPONCM++BoA1yVLdKYJEDZtGgANZizVqe69H6YC0Gy8bnVvrBN02/b7RWeaFw7+D4DO3t/rTBPgZLDQrro5fa5T3ePPlgHgLx2oU90A9d4Sl4k9JhERkSqHGJhERESqHGJgEhERqXKIgUlERKTK8Z8PfpdF6GUNhzepURaAQ20YOlWKYXVJhXVTI/O49ls8BS/USKTQcrwdlq4G76yn0WiI+e0+Bo7VsepWCwB1vorn2x6SE54BGjB0McFhZD2kerJyac9r3I8nWQlsffoPAFIkTPfoiaeFMwAXEx+x9P4JrbT61m7A+IZeaNCQo1Ty3dVThCTHFVlnVvOOOJuY89Hp/VrXcWTTpgxr3BiAyLQ0vg0IIDknB6lEwrft2tHW2Rm5VMqGa9fYdeeO1rrDWjVhiFcTNGiISk5n9sEAUrJzABji1ZgBzT3QV8i5F5PAzAMBFKhUWmu3b+LKvNFd8P1ceAkhkcCn/dri06g2ao2GqIRU5m8PJC0rR2tNgLYt3fj2s250Hb7yddnowa3p6FMPtVrDw7B4fllzkvwC7eo6bkonfDs2IDNDqEfUs2R+mHWgyDoTPvPHoaYFs/+3R+t6OtezZ+K896lubIBKrWblN3/wJCQagMGfdMLv/RbI5FJOH7zOjqV/a6fp4cTkFWOobloNtUrNsgnreXzjKQA+/b0Y+k1/FPpyEiKT+GnUSjJTsrTSrVI9psw0DTuWqBkzS8rM32RY2ks4vKniyQ+UeWpOz4uiQR8Luv/ijMcAS4KWxb6zXu7zbCIW3ST9anyR8sTDEaDW4DbPC7fvvdAUqEk8GqG1bm0ja9Z7jcbfwaNIec8aTXGubsXA86sYfGE1zS2c8bdrWKaei4kFMzw78EHgH3Q/spmVIUGsbd+vyDo9arnT16VsrTfxsLHhY09PBu7eTbetW4lMS2NqmzYADG3cGGdzc7pt2ULfHTsY3awZje3stNJt4GDD6LaeDFu7mz7LtxGZnMqn/q0B6NTQjeHeTfno9/30Xr4VfYWcUW2aaV3nmjZmTB3gi1RSeJPr09qD+rVsGLZgB4PnbSMqIZ1pA3zLcSSghr0Zkz5sh+QN3aYNa+LnU4+PvtjGqM82U91Qj/d7aF/XBo1q8sPsA0wctYGJoza8FZR8/RrQsUujctVT30DBgh0T2Lc2kMndf2HX8pN8tXwkAC061Kdtj6ZM6bGYCf4/0cTbjbY9m5ataajHwhMz+ePnQ0z0/Irt8/cxffunANT1dGHyyo+YN+AXxjX+gujHzxmzYKjW9S2zxxQWFsaJEyeIi4tDKpViY2ND27ZtadSofAdGGx7c0OBUF2wchZPs00PCT5PUDPxEU+TEl5fYW9kY2+nh6GkEQI0WRhjZKN5ZLyUwGrO29igsi/a4qtUzQ8/KEIlUqKuBkzF5MdrdIQAG12rJoegbxOWkFSmXSiQYyhToSeVIJBLkUjl5amWZevkqJV8HHycxJxuAkOQ4rA2qo5BKKVCrcTW1ZHzDlqy48w++DrW1rufdhAT8Nm1CqVajJ5Nha2REVHo6AJ3d3Nh95w4qjYaMvDyOPHxI3/r1uRMXV4Yq3HueQLfFmwVduQwbEyNiUjMA6PNefTZfvEF6Th4Acw8FopBpd181UMiZP6Yri/ee44ePur0ufxqbzLJ9FyhQCj2Ze5HxDGrfROvjoK8nZ9bnPVi16Syzp/Z4XS6TStDTk6OvJ0etVqOnJyc/v+zzBaBQyHCra8eAYd7YO5rzPDqFtctPkhgvHIeatawYNNybHb9fwLOVi9Z1bebrTmxkMlfP3AfgUsBd4qKSAfDu2pgzh66Tl5MPwMk/rtCxX3MuHLlVqqZn5ybEhsVz5fhNAIL/ukZceAIAfiN8+fv308RHJgKw9bu9mFgaaV3fUgPTjh07+OOPP+jSpcvrQJSYmMisWbPo3bs3Y8aM0XpH2pCWCObWhQHIzBpyXwj/DKu/u25GbD4GZjIurY4jNTIXvWoy3vvA+p31HEbWAyD7XmqRcmMPy9c/5yflkBzwDIdR9bXWXRh6FAAvy6IN7q+om/jbeXCy05fIJFKCE59wPuFhmXrR2RlEZ2e8/n1m846cin5CgVpNNbmCpW168L+gYzS21K5H8yZKtRp/V1d+6NyZfJWKZUFBANgbGxObVRiM4zIzcbeyKpeuX31X5vX3J1+pYtUpYa6Ls5U5ltXjWPdhP2yMq3M9MobFxy9opfntiE4cOB/C45ikIuV3nhb2mo2r6TO2hxf7z2v/2PnlRH8OnbxNWERikfLrIc+4djuSfevHoVSqePY8lUMnb2ulaWllzK3r4fy+5jTRz5IZONybuYsGM2nUBgwMFXw9pw+/zP+LOu72WtcTwNHFmtTEDD5fNITaDRzIzsjhtx8OA2Btb8bti49er5sUl4aVnWmZmjXq2pMSl8a0jRNxbVyLrLRsNny9XVhWx56nIc+Ye/ArbJ2tCQ95xtppW7Sub6m3nK1bt7J7924mTZrEwIEDGThwIJMmTWLXrl388ccfWu9EW0pKWSct3xDN27pKeH4jGzd/U7otcqZudzPOLIhGVaCumHAJ5ERkEP7jdSz8amLSVPuLsiTG1+1Aan42HQN+okvgz5jqGTKydmuttzeUK1jt2wdnY3OmBx0HYFHrbmx5cINHaUllbF0yAWFhtFizhhXBwWzq3x8JFHlUeoWqnLkIA++H0WbBWlYHBrN+dH8kEpBLpXjXqcW0XUcZ9OtOTA0N+KxzmzK1BrZrjEqt5lBQaInr1LAyZeP/BnHryXP2nNUugPTt2hSVSsOxwLtvLevu54G9jSl9x6yh75i1xManM/nD9lrpxsWmMfOL3UQ/E3oze3cEY+9ojp29GdNm9OLQvqtEPE0sQ+Vt5HIZzTs04PjOYD7ruYS/Nl1g3uZxKPRkSKVvnzO1uuxzJlfIadn9PY6tD+CTltP5c9VxFhydgUJPjkwhp1VPT5ZNWM/EZl+RGp/G1PXjta9v6X+MHKXy7S5obm4uCsW7PwqVhLk1RDwoPCDpSVDNCPQNyv8Yd3tXEjHXhDt3XqYKE0c9rOoaAlCzpTGX18STFV+AaQ39MrXiD4aReVO4gI3fs8K2n2uJ66ZdjiN220Psh9fDzLv0nsjEuh1pb+MOwNmEB6x5dLrY9fzsGrAw9ChKjYospYrD0bfoZNeQbeFBb607tYkP/jXdAAiIesKeJ7f5rcMAnqQnM+TkLvJUSuyqGdPCpgYuJhZ81KA5pvqGGCv02dRxAKNP7yu2Dp+3bo2fi9CTe5KSwvZbt7j+/DkAe+/e5Xs/P0wNDHiekYFN9cLura2REXGZmSUeg8mdvOno/lI3MYXdl25zI1LQPXA9lDl9/TAxMCAhM5vA0Cdk5wmPG4dvPWBSR69iNSf08qZdE0GzQKnGQE/OrpnDUchk6L/8ecrKP0lKz6Z53RosHNuDLSevsS3geon1BPhoaBvatHB9qavCQF/B70s+QCEXdH9f8gFffr+fdq3qEHD+Pjm5BQD8dfIOU8f5laj7wdh2ePvUBSD2eRr/nHtA4N8hr5dLAKVSRaMmTtR0sqT/YC+MTQypbqTP/MVDmPnF7mJ1R07rhpe/MF5pYl6N6CfxPLwVCQiPcp8vGoKdkxUJMamY25i83s7KzpSk2LRiNUfNHYx3r+YAGFsaEfUghgdXngDCo9y0DROxd7El+XkK4XefkRov6JzYdIafA+eUenzfpNTANGHCBPr27Yu3tzfW1sKjT2JiIpcuXWLq1Kla70Rb3D0l/LlBQ0KMBhtHCRePamjk/W5jS02GWtFkqNBbyUlVcuTzcJLDcrF0NSA+9AWA1uNMtv1cSw1Gr0i/Gk/sjkc4f/EehrVNylx/zaPTJQajN7mfHktnew+uJYcjl0hpZ+vOnbSoYtddevsiS29fBMBUz4AjPUaxL+wuy+/883qduBeZeO379fXvA1w96FarXqlv5ZYFBb1+XGvu6MjyHj3otW0bqbm59HF351FyMmm5uZwKC2NAw4YEhoVRTU+PnvXqMSswsETdVaeCWXUqGABPZ0d+GdKd/iu3k/Yil55N3Xkcn0x6Ti4n7z6ma6O67L0aQp5ShV8DV0Ki44vVXHs4mLWHg98qt7c0Ye/skQydvwOAxi72LJ7Yi282HiMoNLLEOr7it13/8Nuuf94qt7M2YcvyDxkzbSsAj54m4NuqDifOhqJSa2jXqg73Hj4vUXfrhnNs3XAOAGcXaxavGUXo7SjiYtPo1d+T8LAEkhIzGdp72ett/Ls3pm3H+qW+ldu25Djblgg9ZHNrY9YFTsetUQ2ehETj0dIFjUZDXFQylwLuMuzzLhzfGYxapaLTgJac2nelWM0tc/awZY6wT3NbMzaGLqVOMxce33hKo7b10Wg0xIYncGH/Jcb9/AE7F+wnMyULn/5ePLwaVsYRLqTUwNSrVy9atmxJcHAwCQkJaDQamjdvzpQpU7C1tdV6J9pibCZh2DQpv89Xo1KClT2M+LLiLw4NzeX4fuXI1Q3xKHPVyBQSfL9yQKan25eS8fvCQAMxm+6/LqtWxxSHke4V0v3l3jGmN+zJwXafotZouJwUxuawssdXRtR7D4fqJnRxqkMXpzqvy4cF7CYtL/ed63MtJoZfL19m56BBKNVqErKzmXDoEAA7bt/GycyMIyNHoieTsevOHa5ER2ulez0ihnVnr7Bl7EBUKjUJmdlM2S58t7jr0m1MDQ3YN3k4UomEe88TWHTs/Dv/DSD0rJBImNLPhyn9fACIScrgf2sPV0h3275LTB7dgW0rR5NfoCIsIpEl609ptW3E00R+XXKCeT8PRiqTkpSQwQ+zD5S9YRmkJmYy7+Pf+GT+QAyq6VGQr2T++N8pyFNy+VQozvXsWX54KnKFnEsnQzi172rZmvFpfNdvEVNWf4xBdX0K8pTMff8XCvIKuHTkOlY1LFl8di5SqZT4yESWfLxG6/pK/mszAvEjXvEjXhA/4gXxI17xI14REZH/U4iBSUREpMohBiYREZEqhxiYREREqhxiYBIREalyiIFJRESkyvGfTxcQERERKQuxxyQiIlLl+M8TxelyMuSriZBXnznrTBOghVMEAKsedNSp7mR34fOT0VdH61R3U4tNAHifnK4zzeDOCwFw3f2DzjQBwobMAKD2qsU61Q2f/AUA7nN1O8HywRxhgmXTKbrTvbVS0PQeqttjELxLOAZ+HX7UqW7gmW8A6FZPd+0L4PjDhSUuE3tMIiIiVQ4xMImIiFQ5xMAkIiJS5RADk4iISJWjyrmkgO4dTS6eknJsr5AGU98ARk5S4lJPt7Mkwi7lELAslQm7HSqko9FoeLL+CdVqVMOxh6NQptYQsT2CtJA0NCoNDj0csPMrXzrcmQ0H8jQrjp2RhelS+tdsRW/HFuhLFTzIiOGH0H0UaLRz8uhTqyFj3VuhAXJVBcy7fpKQ1DgkwFdNOtDewQ2NRkNEZiozrx0nJe+FVrp969ZnXLPmaDSQo1Qy98JpQhKEvEuftfSmp1s9VBoNdxPjmXEmgHwtXFKGt2jCkOYvnVdS0pl1OICUF0WdUFYM6kliZjbfHz+jVT3fpENjV74f0QWfrwpzXH3UuQU9WzZAJpVy7Op91h6/VG5d3+ZuzJ7YlU4frXpd1r5lHUb18UKhkBGXlMG8X4+TkaVdCpsJEzvi286dzExh/aioZObPO4REAmPHdcCrlStqtYaY6FSWLjlOenrZbjEff92Dtl0bkZkunN/o8EQWTt0FwPDJnfDt3hi1Ss2T0BhWzD5IgZa5z6tcYHrlaNJqkh2OnkZEXckkaFksvVZqnyz/TZ5HSdi1Qc78X/Mxt4Rbl6Usn6tg+c58ndU57bmSfzalQwVj3YuYFzzd8pSsJ1lUq1HtdXlcYBw58Tk0XdgUVa6KkO9CqO5cHWNX4zI1a1W35n/1++JhWpOnTwoNAdrZNGRgzdaMv7KGTGUuC5oMZ0gtH7ZFnCtTs7axBdOb+tH7xG8k5mbT3t6VX33ep+3h1Qx0aYKHhT19TvxOvlrF1006MKOpH/+7XHaOIxczc75p40vPPdtJfJFN+1q1WdOtNz5bNuDlWINeddzpsXsbeSola7v15sPG77H+5rVSNRva2zCmtSd91m4nKy+fr/zb8lnH1sw5Upi87qPWzWnu5Mjx0EelKBWPk7UZU/v6FklP69PAGf+mdRm6aAdqtYZfJ/Wnc1wKJ29qr1/Dzowpw31fG1sAuLvY8sWHHRk7exdxSRl8NrI9Ewb7sOg37XI9NWxYg/nfH+JeaEyR8m7dmlCnrh0Txm2ioEDFuPEdmDDJj59+PFKmZoP3nFg4bSf3bz4rUt6opQvtejRmct8V5OcpmbVqJL1Htmb/b9rl0Kpyj3LFOZr4fFG+xOtvolBo+HiaEvOXPgG166pJSwVlgS5qCwV5ak4uScFnTNnJ28si7lQcNr42WHpZFilPuZ6Cja8NEpkEeXU5lq0sSfxHu7zPA2p6czTmGoFxIUXKuzk0Y2fkBTKUOWjQsOjeQf6OvamVZr5axTdXjpKY+9J9JSUWKwMjFFIpj9OTWHgrkHy16uWyOByqa3ds8lQqpp8+SeKLl7oJcVhXE1xdZBIp+jIZBnI5CqkUfbmcPC16S6GxCXRZuZmsvHzB0cXYiLQXhT0ML+catHWrxe7r2psQvMJAIWfBB11ZfKBoMO/YxI3j1x+Qm68kX6nir8uhdG+hfbJAfT05333SneXbi+p29anP4bN3iUsSDCY27gti++GyE7rBS/eVOrYMGuTF+o1jmDO3HzYv0+lGRCSxfu1pCl763j18GIutbdnnTKGQ4drAgffH+LL60Gd8u2IE1vbCdoJTjAI9AwVyhQyFvpyCPO16S1AFe0y6djSxtgNrO8F0QKOBHevkNPNWI9dRyvIzv6bh0bU6Vs4VF3QZJeSpTg9NL1Ken5yPvmVhbnJ9C31eRGn3aLT4gTDXq7mFW5Fyp2pW3NerztJmo7HSN+F2agSrHh/TSjMmO52Y7MI6znivE4HPH1OgVnMzufBubKIwYEpDH3aG3dBONzODmMw3XF182hMYHkaBWk1Q9DMuRkXyz6ixFKhVPE1LZedd7YwDlGo1fvVcmd9bcF5ZcVZIUGZjVJ0ZXdvz8faDDPYsvx3ZzCGd2PdPCI+fFzV0sDUz5vLDwh5EfFoWtmbaWxd9/bE/fwbe4Ulk0ZtPTTtznjxL4qcv+mBvbUJYVBLLt57VStPS0oibNyLZuPEs0VEpDBrsxbz57zNh3Cbu3Ss8Z0ZGBoz8wIfDf5V9k7KwNeH2pTA2LfmbmPAk3v/Ilzm/jmJyvxXcuhTGjaDHbD0zHWWBiujwRI7tuaz1MSi1x/T8+fNS/1UGleVokpsDK7+XEx8j4eNp2kfu0rhzLAupVEKDThXwltKC4r4akhTjbFEe5BIZLSzr8O3tnYy+tApjhSET3LqUS8NQpmBl637UMjLnmytHiyxzMjJjt98IriVFse1x6Un+39KVy1ndtSe1TM34+vRJAAbW96CGiSlem9bi9fs6ojLS+danvdaagQ/D8P55LavOBbNxRH8UMhmLB3Tnh7/PkZiVXa76AQxq+9J95dLb7ivFOcVo4zoC0N+/CSqVmiNn33Zfkctl+Hi68NPGAEZ9s42UtBdMH9tZK924uHRmfPMH0VEpAPyx5zIODubYvWHTZO9gxtLlw7kbEsWhP8s+Z/HRqcwet5mYcCEw7//tPPZOFtjWMKfz+82xq2HO8LYLGO6zgLjoVMZ+3aMMxTf+1tIWjh8/noiICGxsbN66OCQSCYGlJJkvD7p0NCmOpARYMkuBg5OGb38pQO/dZAC4tCOD8KvCo4BUDso8Dbs+T0BVoEGZL/zca5YlRpZle0492/eMlBtCQ7FoZoHTAKdi19O31Cc/tXBMLD81Hz0LvWLXHevqj4+14GV3MfE+G8ICil0vKS+DcwmhvFAJBpInYm8yxrVkJ4/PPXzxcxTyhgfGPGbP01tsaDuQsIxkhp/ZQZ6qMNi3sqnFitZ9WX//Ehsfln6XnNqyNZ1qC0YPp8LD2HMvhI09+/IkJYWhB/e+1u3i6sahR/fJLhCewXeF3mGub/H1ndLem471hN5nWGIKO6/e5kaUcCPdfzOU73r44eFgSw0zU6Z3EZx3rYyqI5NI0JPLmHW4+DGbid29ad/opfuKSo2BQs6er4cjl8nQf/nz5DV/EpuagbVp4c3KxtSI+LSSjU/HDmiNj6dwDJRKFfr6Crb8OPK1+8qWH0fyxaIDJKVmEfYskZSXA81Hzt1l1cyS091+OLot3q2Fcxb7PI2LFx9xKqAw4EkkoFQKN/ymTZ2YObsve3ZfYu8fxRsRAIz81B+vjkL7iotKIehUKKcPvdG7kkhQFaho7d+QM4dvkZMttNvjf1xh0iztv/ooNTDt2rWLYcOGMWfOHDw9PbUWLS+6dDT5N1kZsOALPdp2VtF/pPY+9yXRargJrYa/7YCSEa9k56cJDF1mo7WW0wCnEoPRm1h4WpBwPgGLZhaoclUkXUrCZXTxLqwbwgJKDEZvcjr+Ln52jfgr+gp5aiW+Ng25n16yacCyu+dZdlcYuDTVM+BQ5zHsD7/DytCLRdZrZunIGp/3+SzoT87HPS2zHkuvBLH0iuC+YqpvwOHBI9h3P5QVV4u6nIQmJNDVpQ4HH9xDpdHQ1aUOt+KKt3lfeTaYlWdfOq84ObJkQHf6rt1OWk4uvRq58zghmZtRz+mwbOPrbSa3a4V5NcNS38qtORbMmmNvu684WJiwb8ZIBv8kuK+cDXnK+G6t2PdPCCq1ht5eDfjr8r0SdTfsC2LDvretuOysTNjx8yhGfbMNgDOXHzF5eDs2/3mZjKxc2reow/2wkh2ON2+6wOZNwltYZ2crlq0Ywd2QKOLi0undpxlPnyaSlJRJg4aOzP3+febPO8TVq6Wfs20rAti2QmhfterY8vOOCYRejyA+OpUew1oR/jCWpPgMwu49p42/B4GHbqJWqWnj35AHt4t39imOUgOTkZER8+fPZ+/evZUamN5E144mgYdlJCXAtYtSrl0s1Pjm5wKMy3ZYqhLY+dmRG5/LrRm30Cg12Ha0xbR+xQbbD0QFY6IwZFOrKUglUh5mxLDi4dGyNwSGuzXDoZoJnWvUo3ONeq/LR57ZyWeNfJEAXzbpwJdNOgAQnZ3GxIslW0O9YoRHExyMjOni4kYXl8IxseGH9rL6+mVm+rQnYPiH5KlUPEhKZPa5snvs15/FsPbCFbZ+OBCVWnBe+WTPX1r9ne/K+btPqeNgyY7/DUMhk3ImJIzDV0oOTNpy8cZTrC2M+XX2YKQSCXFJGfyw/oRW20ZEJLFqRQDzfxiIVCohKTGTBd8LzjYfftgWgI/Htefjce0BwXhzThnuLJGP41kz/xDfrRklOLrEpfPTNGGqwO61Zxj3TQ/WHZ1KQb6K8IexrJ77p9Z/63+e9kT8iFf8iBfEj3hB/IhX/IhXRETk/xRiYBIREalyiIFJRESkyiEGJhERkSqHGJhERESqHKIZgYiISJVD7DGJiIhUOf7zj3jfD5qkM639rYX8N5U13+hp9LtnNSgOlxrCTOWdT7x0qjvMTfjsY+qtITrTXNp0NwD+Z6fqTBMgoL0wH8j9wDyd6j7oPxsA5w0/61Q3YuyXANT5UXfzmB5/IxzTRl/ods5VyGJBt+WHS3Sqe2XzNADad1ukU92zx78qcZnYYxIREalyiIFJRESkyiEGJhERkSqHGJhERESqHGJgEhERqXJUmdS6Go2GmN/uY+BYHatutQBQ56t4vu0hOeEZoAFDFxMcRtZDqld2EraS0JWbyZsEXZTzy0+GHDicqTPNO6fzCTogJKRT6EvoNsEQhzoVP11xV9J4vDcWJBIU1WU0Gu9Edbt3z5z3pftQwrNj2Rd19q1lcxqOJjk/nVWPS0+f8Sa9ajbiozrer91XFtz+m7tpwtvM8fV86OPUGJlEyuGoEFbdL9s4AaCvWwPGN26BBshRFvBdUCAhSYLzSlfnunzS1As9mYyYrAymnT1GWp52riMjPJsw7D3BfeVZajrfHhfcV/TlMr7r3JFG9nZIJRJuP4/lu5OnyVNqlw9saJsmDGrdBI1GQ3RyOt/tDSAlKwepRMKM/h1o7lIDgAsPwll8+EIZam/Trpkrc8Z2pePE1W8tmzq0HTVtzZm2TLsUJRM/7kD7tvXIzBQcVZ5FpzJvoZD1Y/jgVnTxa4hMJiXg9D027/hH6zpWiR5T7vNsIhbdJP1qfJHyxMMRoNbgNs8Lt++90BSoSTwa8c770ZWbyZvEREvZuM4ATcUy/xYhKVpFwO85DJ9nxIRVJvgOMWDPgvKnf/03qnw1t1dF0uwLF9oucse2uSn3NpecHK40nKrZsKjJJHytmxa7fFDNjniYFp/MriRqG1nyZaNOjA3aSb/T61nz4AIrWg0CwNfWjS6O9Xn/9AZ6nVqDl5UzXR0blKnpYmrODK92fHB8H90PbGHlzWDW+vcFoJGVLfPa+DHx1CG67N9MeHoqX7Zoq1VdG9rZ8FFLTwZt202PjduISE1lqm9rACa19kImldLrt230/G0bBgo5E7xbaqXboIYNo9p7MnLlbvr/so3IpFQmdxV0ezWvj7O1Of1/2caAxdtp7lKDzo3raKX7ipq2Znw6uF2x6QTMV9wAACAASURBVH87tahLV+/65dLzaODAvIV/8fHkLXw8ecvroOTVwoX2besxbspWRk/4naZNnGjftl4ZaoWUGZhOnTrFtm3bePasqD3Lnj17yvUHlEZKYDRmbe0xbWFbpLxaPTOse9VGIpUgkUowcDKmIEm7u9m/0aWbyStyc+HnHw0ZN/Hd6lQScgX0+rQaxhbC6XGoIyMrVYOqoGIRVaPWoNFoUL4Q7tzKXDVSxbvlDu/t6MPJuCucT7z11rImZm40t3DnyPO3szKWRr5ayawbR0jMFdLQ3k17LrivSKR0cnDnSNRdclQF5KtVHIi8RW+nsg0E8lUqvj5/gsScl84rifFYGwrOK/3cGrDnYQjRWYIBwrLr/7D2dslpZd8kNC4B/3WF7it2xkak5gjt4GpUNL8GXUYDqDUa7sUl4mCqXVbCe9EJ9PxxM1m5+ejJZdiYGpGW/TKVs0SCoZ4CPbkMxct/2vbCQHBfmTuuG8t2nX1rmbO9BSO7N+e3v7T3v1MoZNRxtWXw+y3ZuPpD5n7bBxtrwVKsrXcdTp25T25eAfkFKv4+GYJ/x4Zaa5camH755Re2b99OREQEQ4YM4dChQ6+X7d69W+udlIXDyHqYt357MqOxhyX6doK/Wn5SDskBzzD5V/DSFl26mbxi5VIDuvXMp7ZLxVP2vomZrYy6LYV6ajQaTmzIoZ6XAtk7BpFXyA1keHxck+BZjwicEELkiUTchzu+k9aqxwc4Ff+2p5ulngmT3Pqx8P521JSvGxnzIp1zcY9f/z69UWfOxD6kQKPGvpoJcTmFDipxORnYGpZ9sUdnZXAmqjBd7MxWHTgV+YQCtZraphbIJVI2+PfleP9RzGvjT3aB9n6DSrWaTnVcuTB5LM1r1mD/HcGY4GL4MyJS0gBwMDFmVIv3+PuB9p5ySrWajh6unJo9Fk+XGvx5VdA9dPUeGTl5BM4ey5k543iWlMa5e2WnL37FNx924uDZOzyJLurqYqivYO64rszdeILsXO3/fksLI27cjmTDpvN8/Mlm7j2IZcGc/gDYWBuTmFR4vhKTMrG2KtsH8RWlBqZz586xceNGZs2axc6dO1m+fDnHjx8HinfuqCxyIjII//E6Fn41MWlqVe7tK8PN5MghBTIZdOmmI4O6YsjP1bDvxxekxKrp/alhhfUynuXwZH8cbRfXx29tI9z62XJj8VOdnUuZRMqMBh+w5slBUvIzyt6gBAxlCpa1HICTkQUzbwhGmVKKcR4px/OzoVzBar/eOJuYMf2CkI5WLpXiV8uVGRcD6H5gC4k52SxsWz6nmFOPw/BavpaVF4LZNLh/kVo2tLNh14hBbL9+izNPwsule/puGL6z17LmRDDrxvVHIoGJnVuRmpVDu+/W0en7DZhWM+CDds200nu/o+C+cvjC264uM8d05o9Tt3gak1yuOsbFpzN99n6iYgRDjT37r+Bgb4adrWmxLj5qtfbnq9TRVI1Gg+Tls6izszPr1q1j9OjRWFhYvC5/F+IPhpF5U4jaxu9ZYdvPtcR10y7HEbvtIfbD62Hmrb0tti7dTIoj4IQeeXnwybjqFBRAfr7w87wfXmBpVf4L/cy2HB5eFoJcPS8Fzbros2teFlY1ZYz60QiF/rsd70d/xBJ/TfCAy47Nxb6V+evB7lpdrLm3JYaCTBV6JqUPrI9y7oq3lQcAwUl32RLx91vr1DWuiZ2BBRPchDEccz1jpEjRkypY8rD4R/8p9dvT0b4uAKdjH7Ev4gZrvIcQlpnEqPNbyVMLLinPc9KxNij0ZrM1NCYup/iXDVM92+BfS8gZHhD5hD0P7vBbl/48SUtmyNE9r51XEl5k8TA16fVj3t5HIezsPrjEY/BZW2/86gjjZk+SUthx4zbXowX3lX13QpnX1Q9TQwPScnLpUb8u33XxY97J0xy+97BETYBPunjTvqGg+zQ+hd1Bt7kZLugevBLKrAF+mBga4NfIjR8PnkGpUpOlyueva/fwb1yHreeK9+0b1681vu+9dHVRqjHQk7N93gjB1eXlz18sO0TTuo7UsjdnaJdmmFQ3wMhQn6VT+zF16cG3NEeP9KGNl3C9Po9L50LQIwJOF+YzlyBBpVKRkJCJhXnh+bKyMiIxSfuXQ6W2xq5duzJy5EimT59O48aNqVOnDsuXL2fy5Mnk57+7xbZtP9dSg9Er0q/GE7vjEc5fvIdh7fI5B+jSzaQ4lv9aOBgdHydhwkdGrF7/7gPUHUYa0mGk0CvKyVSz/rMsmnTSo/0wgwrVs+4ge+oOEh6Tk+5mErL2GXlpBeibKYi7mk41G70ygxLAloi/iw1Gb3I/I5Lhlwq/gRvp3AVTRfVS38qtvH+WlffPAmCqMGB/x7EcjLzN6gdFraRPP3/EJ/V9+SP8BkqNmn61mnIw8u3xLYCl1/9h6XXhDZCpvgFH+o1k36NQlt8oOuZ1LPwRM7zas+pmMGl5uXR1rsudpJJdR5ZfCGb5BcElpXkNR5b16U7v37eTmpNL74buPEpMJi0nl6716jDLvwOjdx/gblx8iXqvWH0imNUnBN1mtR35eUR3BizZTlp2Lj2aufMkLpn0F7ncj0mgS9O6XA2LRi6V0r6BC3cii3eKAVh/MIj1B98e57O3MmHX/A8YMXs7AD2mrn+9rIdPA/ya1y3xrdymbRfZtE1wxqldy4oVPw8jJDSGuPh0+vRoSlh4AolJWVy89JgPh7XhyPHbqNRqunZqxN+nQorVLI5SW+TkyZPx9PSkevXCRyBPT08OHDjA77//rvVO3pX4fWGggZhN91+XVatjisNI7e2W/y9y9Wg+6YlqHgTl8yCo8AbwwQ9GVDN59xepVh7G1O5lw6W5j5HKJSiM5Hh+Wb43Z5XJEJfm2FczpZODO50cCs/x6IvbOBP3iLqmNvzR4SMUUhmnnz/iz2dl23qPqN8Uh+omdHGuQxfnwjdYw47uIfBZGPbVjdnTcwhSiYSYrAy+Ol968H3FtegY1gRdYftwwX0lPiubSfuFN1JftG+DRAI/dO/0ev3r0c+Ze7Jka6hX3AiPYf2pK/w+UdBNzMjms02C7qJD55jRrwN/fT0KlVrD5cfP+P302+N8/xXhkUmsWHOKH7/rj1QqJTEpk+9/Eh69gy+H4eJszZrlgj/eP5cec+LU24+RJfGf52MSswuI2QVAzC4AYnYBMbuAiIjI/ynEwCQiIlLlEAOTiIhIlUMMTCIiIlUOMTCJiIhUOUSXFBERkSqH2GMSERGpcvzn+ZiaHp2lM61bPb4HYPTV0TrTBNjUYhNQefON0p7X1KmumUMUACfCy04Doi1dagufGcy721tnmgCzPYTJgrqczwaFc9p02b6gsI05b/1JZ5oRH3wNgOsS3c43CpsmzDdqMEO386Pu/SDMj2o2Xre6N9aVPEdO7DGJiIhUOcTAJCIiUuUQA5OIiEiVQwxMIiIiVY4qY0YAMK9xP55kJbD1qZCyQoqE6R498bRwBuBi4iOW3j9RLk2NRsOT9U+oVqMajj2EbI0atYaI7RGkhaShUWlw6OGAnZ/2uZ7epDJMA548lbB4hR5Z2RKkUg3Tp+VTv17FZ3WEXtZweJMaZQE41IahU6UYVq9YVkyA1Mg8rv0WT8ELNRIptBxvh6Xru6drqUxjCl22sb61GzC+oRcaNOQolXx39RQhyUXTpsxq3hFnE3M+Or1f6zqObNqUYY0bAxCZlsa3AQEk5whmBN+2a0dbZ2fkUikbrl1j152yMyy8YlirJgzxEswTopLTmX0wgJRswURgiFdjBjT3QF8h515MAjMPBFCg0j4za/smrswb3QXfz4WXEBIJfNqvLT6NaqPWaIhKSGX+9kDSsnK00qsSPabaRtas9xqNv4NHkfKeNZriXN2KgedXMfjCappbOONvp33e4BcxLwj9MZTky0Uz88UFxpETn0PThU1p/H1jYv+OJTOs/A4nlWEakJsLn35pwIghBWzbkMuYkUrmLHh3F5NXZKZp2LFEzZhZUmb+JsPSXsLhTRUPdso8NafnRdGgjwXdf3HGY4AlQctKzhFUFpVlTKHrNuZiYsEMzw58EPgH3Y9sZmVIEGvb9yuyTo9a7vR10b69AnjY2PCxpycDd++m29atRKalMbVNGwCGNm6Ms7k53bZsoe+OHYxu1ozGdtrdUBs42DC6rSfD1u6mz/JtRCan8qm/YHLQqaEbw72b8tHv++m9fCv6Cjmj2miXGROgpo0ZUwf4FjE46NPag/q1bBi2YAeD520jKiGdaQN8tdYsMzBFREQQHy80kr179zJ//nyOHTum9Q60YXCtlhyKvkHA87tFKyeRYChToCeVo5DKkUvlrzMaakPcqThsfG2w9LIsUp5yPQUbXxskMgny6nIsW1mS+E9iuetdGaYBl6/JcHRQ06aVkIbUt42KBXPy3lnvFQ9uaHCqCzaOQuPx6SHh2mlNhdPqxt7KxthOD0dPIVthjRZG+Hzx7uliKsuYQtdtLF+l5Ovg44UmB8lxWBsIJgcArqaWjG/YkhV3tLcsAribkIDfpk1k5QsmB7ZGRqTmCL2Mzm5u7Lt7F5VGQ0ZeHkcePqRvfe1cTe49T6Db4pfmCXIZNiZGpL0Qjl+f9+qz+eIN0nPy0Ghg7qFA/rp1rwxFAQOFnPljurJ4b1ErraexySzbd4GCl2YJ9yLjsbfUPtljqc8cmzdvZtu2bajValq1akVsbCz+/v7s37+f8PBwPvnkE613VBoLQ48C4GVZNGnZX1E38bfz4GSnL5FJpAQnPuF8QulpSt/EZZSglx6aXqQ8PzkffcvCXoi+hT4vol6Uu95mtjLMbIVHCV2ZBjyLkmBpoWH+Ij0eh0kwNoLJ4989W+gr0hLB3LqwXmbWkPtC+GdYgVToGbH5GJjJuLQ6jtTIXPSqyXjvA+t31nMYKVj8ZN9LLVJu7FF4c3llTOEwSnurIV23sejsDKKzC/Oaz2zekVPRgslBNbmCpW168L+gYzS2LP8QgVKtxt/VlR86dyZfpWJZkJCF0t7YmNisrNfrxWVm4m6lfQ58pVqNX31X5vX3J1+pYtWpvQA4W5ljWT2OdR/2w8a4OtcjY1h8XDu/um9HdOLA+RAexxQ1OLjztLDXbFxNn7E9vNh/XvvHzlJ7TPv37+fYsWNs376dv//+m3Xr1jF8+HDWrFnDiRPlG+t5F8bX7UBqfjYdA36iS+DPmOoZMrJ26wrrFtdLKC55urbo0jRAqYKgyzL69lSyZV0eA/spmTrdgApkMgagpI6R9N29QwVdJTy/kY2bvyndFjlTt7sZZxZEoyrQodHeG1TUmOLfVLSNGcoVrPbtg7OxOdODBKOORa27seXBDR6lJZWxdckEhIXRYs0aVgQHs6m/YHJQnBecqpw93sD7YbRZsJbVgcGsHy2YHMilUrzr1GLarqMM+nUnpoYGfNa5TZlaA9s1RqVWcyio5MyUNaxM2fi/Qdx68pw9Z29rXc9Se0xqtRo9PT0cHR0ZM2YM+vqFvQxVOQbG/s3Euh1pbyOkTj2b8IA1j04Xu56fXQMWhh5FqVGRpVRxOPoWnewasi28eL+yZ/uekXJDcGywaGaB0wCnYtfTt9QnP7XwSs9PzUfPQk+ruleWacArrCw1ODtp8GggXNjtfFT88AvExEqoXevdH7vMrSHiQeH26UlQzQj0Dcpf39u7koi5Jty58zJVmDjqYVVXCMg1WxpzeU08WfEFmNYoe2yssowpKqONTW3ig3/NlyYHUU/Y8+Q2v3UYwJP0ZIac3EWeSoldNWNa2NTAxcSCjxo0x1TfEGOFPps6DmD06X3F1uHz1q3xc3lpcpCSwvZbt7j+XDAj2Hv3Lt/7+WFqYMDzjAxs3khzbWtkRFxmyWOjkzt509H9pW5iCrsv3eZGpKB74Hooc/r6YWJgQEJmNoGhT8jOE66Jw7ceMKlj8V89TOjlTbsmRQ0Ods0cjuKlwcGumcOZsvJPktKzaV63BgvH9mDLyWtsC7heYj2Lo9TA1LlzZ0aMGMHWrVuZMmUKAA8ePGDmzJl069atXDt6kzWPTpfYUN7kfnosne09uJYcjlwipZ2tO3fSokpc32mAU4nB6E0sPC1IOJ+ARTMLVLkqki4l4TJau9zXlWUa8IrWXipWrNHj/kMJ9etpuHlbikQCDvYVGwty95Tw5wYNCTEabBwlXDyqoZH3uwXRJkOtaDJU6K3kpCo58nk4yWG5WLoaEB8qPBIb2Wjn31dZxhSV0caW3r7I0ttCIn5TPQOO9BjFvrC7LH9jHCnuRSZe+359/fsAVw+61apX6lu5ZUFBrx/Xmjs6srxHD3pt20Zqbi593N15lJxMWm4up8LCGNCwIYFhYVTT06NnvXrMCgwsUXfVqWBWnRJMDjydHfllSHf6r9xO2otcejZ153F8Muk5uZy8+5iujeqy92oIeUoVfg1cCYku3kRh7eFg1h4Ofqvc3tKEvbNHMnT+DgAau9izeGIvvtl4jKDQyBLrWBKlBqbPPvuMq1evIpMV9vf19PSYMmUK7dq1K/fOyssv944xvWFPDrb7FLVGw+WkMDaHld+r/d/Y+dmRG5/LrRm30Cg12Ha0xbR++R16K8M0wNICFn2fx8/L9MjJlaBQwMJ5eehr16ErEWMzCcOmSfl9vhqVEqzsYcSXFX8pa2gux/crR65uiEeZq0amkOD7lQMyPd2+8K0sY4p3bWMj6r0nmBw41aGL0xsmBwG7Sct7d2fmazEx/Hr5MjsHDUKpVpOQnc2El0azO27fxsnMjCMjR6Ink7Hrzh2uRGtn8X49IoZ1Z6+wZexAVCo1CZnZTNkufLe469JtTA0N2Dd5OFKJhHvPE1h07HwZiqUzoZc3SCRM6efDlH4+AMQkZfC/tYe12v4/T3sifsQrfsQL4ke8IH7EK37EKyIi8n8KMTCJiIhUOcTAJCIiUuUQA5OIiEiVQwxMIiIiVQ4xMImIiFQ5RJcUERGRKofYYxIREaly/OeJ4ipjopr3yek60wQI7rwQgKm3huhUd2nT3YBuJ0JC4WRIdVxdnWlK7R4BcPWZs840AVo4RQCw6kFHnepOdhc+P6msyba6bGOv2pfr7h90pgkQNmQGALVXLdapbvjkLwBwn6vbCZYP5ogTLEVERP4PIQYmERGRKocYmERERKocYmASERGpclQJl5TKcpt4xcyGA3maFcfOyMJ0Fv1rtqK3Ywv0pQoeZMTwQ+g+CjTlT34XdyWNx3tjQSJBUV1Go/FOVLermHlAZbmZBJyHVZtAKgUTI/j+K3ByrJjmxVNSju0V0uLoG8DISUpcdODo8iZhl3IIWJbKhN0OFdKpLMccXbevPrUaMta9FRogV1XAvOsnCUmNQwJ81aQD7R3c0Gg0RGSmMvPacVLytEsL3bdufcY1a45GAzlKJXMvnCYkQci79FlLb3q61UOl0XA3MZ4ZZwLI1yIZ5PAWTRjS/KXzSko6sw4HkPKiqBPKikE9SczM5vvjZ7SqJ1SBHlNluU0A1KpuzcrmY/Gza1SkvJ1NQwbWbM2n1zYyLGgp+jIFQ2r5lFtfla/m9qpImn3hQttF7tg2N+XeZu3y45REZbmZ5ObB1wtgxfdw8Dfo0AYWrKiY5vMoCbs2yPnyhwJ+WFdAn2Eqls/VLkGctqQ9V/LPpnSo4CGoDMecymhftY0tmN7Uj9HndtPrxG+sDv2HX33eB2CgSxM8LOzpc+J3uv+9kcisFGY09dNK18XMnG/a+DLqrwP02LONVdcusaabkNLGy7EGveq403PPdrru2oKRQo8PG79XpmZDexvGtPZk6O+76b1mG5EpqXzWsWha4o9aN6f5O9z9yhWYFi5cWO4dlEVluU0ADKjpzdGYawTGhRQp7+bQjJ2RF8hQ5qBBw6J7B/k79ma59TVqwWVE+UK4syhz1UgrYEQAledmolIJeb+zXrpLvcihwsnnFAoNH09TYv7SJ6B2XTVpqaAsqJjuKwry1JxckoLPmPIn8fs3leGYUxntK1+t4psrR0nMfXk9pMRiZWCEQirlcXoSC28Fkq9WvVwWh0N17Y5NnkrF9NMnSXzxUjchDutqwnUmk0jRl8kwkMtRSKXoy+XkadFbCo1NoMvKl84rMhm2xoXOKwBezjVo61aL3de1NyF4RYmPct98881bZadPnyY9XXAc+fHHH8u9s+KoTLeJxQ+EpGTNLdyKlDtVs+K+XnWWNhuNlb4Jt1MjWPW4/JZUcgMZHh/XJHjWIxTGMjRq8J5XsblEleVmUr0azJkGQz8BMxNQq2HHqgpVFWs7sLYTcpNrNLBjnZxm3mrkOuo0nfk1DY+u1bFyrrhgZTjmVEb7islOJya7sI4z3utE4PPHFKjV3EyOeV1uojBgSkMfdobd0E43M4OYzDeuM5/2BIaHUaBWExT9jItRkfwzaiwFahVP01LZeVc74wClWo1fPVfm9xacV1acFZxXbIyqM6Nrez7efpDBno3KUHmbEntMZmZmnD17Fnd3d1q2bEnLli2pVq3a6591TWW5TRSHXCKjhWUdvr29k9GXVmGsMGSCW5dy62Q8y+HJ/jjaLq6P39pGuPWz5cbipxXq3VSWm8mjMFizFY5sgfMHYPwI+Gx2yfsrD7k5sPJ7OfExEj6epr3vX2ncOZaFVCqhQacKRGMt0LVjDuimfRnKFKxs3Y9aRuZ8c+VokWVORmbs9hvBtaQotj0uX5J/Q7mc1V17UsvUjK9PnwRgYH0PapiY4rVpLV6/ryMqI51vfdprrRn4MAzvn9ey6lwwG0f0RyGTsXhAd374+xyJWe9mAFtij+nrr7/G19eXZcuWMW3aNLy8vNiyZQv9+vUraROtqSy3ibGu/vhYC15jFxPvsyEsoNj1kvIyOJcQyguVYCR5IvYmY1y1e1Z/9Ecs8deEO1p2bC72rcxfD3bX6mLNvS0xFGSq0DN5t/cKunQzeZOLV+E9j8LB7mH9YOFqSEsHc7N3101KgCWzFDg4afj2lwL0KjDuf2lHBuFXhUcBqRyUeRp2fZ6AqkCDMl/4udcsS4wsy47SleGYU1nt63MPX/wchbzhgTGP2fP0FhvaDiQsI5nhZ3aQpyoM9q1sarGidV/W37/ExoeXS9QEmNqyNZ1qC0YPp8LD2HMvhI09+/IkJYWhB/e+1u3i6sahR/fJLhCewXeF3mGub/H1ndLem471hN5nWGIKO6/e5kaU4Lyy/2Yo3/Xww8PBlhpmpkzvIjjvWhlVRyaRoCeXMevwqVLr/IpSrx5vb2/q16/PnDlzOHv2bIUsm96kstwmNoQFlNhY3uR0/F387BrxV/QV8tRKfG0acj9du0HruoPsqTtIcJpNuptJyNpn5KUVoG+mIO5qOtVs9N45KIFu3UzepEFd2HEQklLAygICL0IN+4oFpawMWPCFHm07q+g/suJto9VwE1oNf9sBJSNeyc5PExi6zEZrrcpwzKms9rXs7nmW3RWS/5vqGXCo8xj2h99hZejFIus1s3Rkjc/7fBb0J+fjnpZZj6VXglh6RXBfMdU34PDgEey7H8qKq0VdTkITEujqUoeDD+6h0mjo6lKHW3HF27yvPBvMyrMvnVecHFkyoDt9124nLSeXXo3ceZyQzM2o53RYtvH1NpPbtcK8mmG53sqVeQWZmZmxfPly9u7dy8OH2rvgaktluU2UxoGoYEwUhmxqNQWpRMrDjBhWPDxa9ob/wsrDmNq9bLg09zFSuQSFkRzPL7WzgSqJynIzadUMxgyBUZ+BQgGmxrBqQcU0Aw/LSEqAaxelXLtYWMdvfi7AWHs36P+v6Mox500q0r6GuzXDoZoJnWvUo3ONeq/LR57ZyWeNfJEAXzbpwJdNOgAQnZ3GxItlT6EZ4dEEByNjuri40cWlcExs+KG9rL5+mZk+7QkY/iF5KhUPkhKZfa5kW6hXXH8Ww9oLV9j64UBUasF55ZM9f2n1d5bFf572RPyIV/yIF8SPeEH8iFf8iFdEROT/FGJgEhERqXKIgUlERKTKIQYmERGRKocYmERERKocYmASERGpcoguKSIiIlUOscckIiJS5fjPE8W5LlmiM62wadMEzUqaqOZ/tuQJYO9CQHthgtq8u711qjvbQ5htq8vJkK8mQp6L0N2kTYB2zsLEzafR9jrVdakhfEKx84mXTnWHuQnfo+lysu2ribaV1b7cD8zTqe6D/rMBcN7ws051I8Z+WeIyscckIiJS5RADk4iISJVDDEwiIiJVDjEwiYiIVDmqhEvKyKZNGda4MQCRaWl8GxBAck4OUomEb9u1o62zM3KplA3XrrHrjvb5gyvLbeIVX7oPJTw7ln1RZ99aNqfhaJLz01n1+EC5NF+RGpnHtd/iKXihRiKFluPtsHQ1eCetV1SWo4lGA5sXg2Mt6DywwnJFCLoo55efDDlwWDuDAG24czqfoANCSh2FvoRuEwxxqFPxS0HXjjm6bl+9ajbiozrer6+HBbf/5m6a8NJgfD0f+jg1RiaRcjgqhFX3z2ml2detAeMbt0AD5CgL+C4okJAkwXmlq3NdPmnqhZ5MRkxWBtPOHtM6ldH/9x6Th40NH3t6MnD3brpt3UpkWhpT27QBYGjjxjibm9Ntyxb67tjB6GbNaGynXe7vynKbAHCqZsOiJpPwtW5a7PJBNTviYfrueZmUeWpOz4uiQR8Luv/ijMcAS4KWFZ+4S1sqy9Ek9hks+Rquna+w1FvEREvZuM4AjVp3mknRKgJ+z2H4PCMmrDLBd4gBexa8W/rXN9GlY05ltK/aRpZ82agTY4N20u/0etY8uMCKVoMA8LV1o4tjfd4/vYFep9bgZeVMV8eyU/O4mJozw6sdHxzfR/cDW1h5M5i1/n0BaGRly7w2fkw8dYgu+zcTnp7Kly3aal3fUm8Td+7cofHLnkxwcDDnzp1DLpfj7+9PkyZNtN5JadxNSMBv0yaUarXgtGBkRNRLw4PObm7sjBvqAwAAIABJREFUvnMHlUZDRl4eRx4+pG/9+tyJiytDtfxuEyPqeGpd596OPpyMu0JiXupby5qYudHcwp0jz4MwVhhqrfkmsbeyMbbTw9HTCIAaLYwwsqlYECnN0aQi5gFn/oLWncFC++SSWpGbCz//aMi4ibn8tODdjmNxyBXQ69NqGFsI92SHOjKyUjWoCjTIKuBwo0vHnMpoX/lqJbNuHCExNwuAu2nPhetBIqWTgztHou6SoxJS6x6IvEVvp0b8HXOvdE2Viq/Pnyh0OEqMx9pQcF7p59aAPQ9DiM4SDBCWXf8HMwPt61tqYJozZw4HDx5kx44d7N69m/ffF3ocs2fPZuDAgYwYMULrHZWGUq3G39WVHzp3Jl+lYlmQkA7U3tiY2Kys1+vFZWbibmWllWZluU0Ar7vP75nXKVJuqWfCJLd+fHNnHT0cvLXW+zcZsfkYmMm4tDqO1Mhc9KrJeO8D63fWg8pzNBk2Wfj/wa2K6fyblUsN6NYzn9ouuknn/AozWxlmtsLjrEaj4cSGHOp5KSoUlOD/tXfe4VFU/fu+t6U3kpBGgNBLgNAhEGpIAoQuKL0JgooF0VdERF8URBABqYqAdJAmgigEQgsJXXoPEEjvve7u/P4YSAikbLIb3/3+rrmvi4vdyc4zZ3fPfPbMmZnnMWxiTlX0r6jsNKKyi/aHmc39OR5zlwJBi6uFDWcTHhX+LTYnHWfz8i1IIzPTCwsPwOyOPTgaISYc1bG1505yAmv9BuFubcud5ES+Phusc3t1OpT77bff2LRpE+PHj2f8+PFs3bqVLVu26LwRXQgKD6fd6tX8GBbGhiFDkAFy2audRVPBO2iqKm3iZRQyObOajmX1g30k56eXv0IZCGqIvpxFfT9b+iz0oGFfO47Pi0RToP8xTVUkmhiag/tVKBQQ0MdAAXUlkJ8rsPvbbJJjtAx4X/8RWVUk5ryIofqXuULF0vZDqWVlz+zLBwCQ8+p+pq3A8bO5UsVK3wF42Ngx8/RhAJRyOb616zErJIi+ezeSkJPFgi66J8WUOWJSq9VotVocHBywsLAoXG5iYoJcXvnpqQ87dcK3rniM/CA5mS1XrnApWkxa2HXjBl/7+mJrZkZ0ejpOlkXxPc5WVsRmlD4JWlVpE+M8euPt2AyAsMQbbHz89yuvaWhdExcze6bWF4+xq5lYI0eOiVzFD3d3lvuZXN2eSNRFcXSYl6HBpoYJjg3FHaZme2vOrY4jM64AW/fKT6YaMtGkKgk6bEJeHrz7liUFBZCfLz6eOz8bB8eK7+jHN+dw95xY5Bp1UNE6wJTtczNxrKlg3LdWqEwrN1oyVGJOVfWv95p0p6erOGoLjrnH7seXWe09nPCMRMad2kSeVtwfonPSqG5mVbies7k1sTkl72fT23TGr/azhKOIB+y8c411AUN4kJrE8D93Fu5j8dmZ3E1JLDzM23XvOtv6vlHm5/AiZX5i1apVo1u3bshkMr788ksWLFhAWFgYixYtonfv3jpv5GWWhoYWHq61rVGDZYGB9N+8mZTcXAY2bsy9pCRSc3M5Gh7OUE9PjoWHY2FiQr9GjfjiWOkm6VWVNrHx8d8ldpYXuZ0ewaizRbcCjPEIwFZlqfNZE68RjniNEA9Tc1LUHPzwEUnhuTjUMyPupni2UJ95JkMnmlQly1YVTUbHxcqY+qYVK3+u/AR1jzHm9BgjFvmcDC0/f5CJVy8Tuo/U7yynoRJzqqp/Lb99guW3TwBgqzJjT8/J7Iu4yso7xc9UBEff490mXfnt0WXUgpbBtVuyL6LkY/Mll86w5JKYZmRrasbBwWPYfe8myy6HFnvdoUf3mNWhOyv+CSM1L5feHg25llj+3PBzyvzUNm3aBMDDhw9JTxeHjyYmJrz//vt0795d542UxcWoKFadO8e2119HrdUSn5XF1P37Adh69Sq17Ow4OGYMJgoF269d43ykbmc6qipt4t/AvJqSrv+pwYW1cahztShUMrr+xw2FSeVHqf8/JJoYggt/5pOWoOVOaD53Qovy5MbOt8LCpvKfb1Uk5hiS4XXb4mphSy+3xvRya1y4fELIZo7H3qOhrRO/9XgTlVxBcPQ9fn9S/mU5o5u0FBOOPBoQ4PFCwtGfOzn2JBxXS2t29huOXCYjKjOd/5wqu/i+yL9ueyLdxCvdxAvSTbwg3cQr3cQrISHxfwqpMElISBgdUmGSkJAwOqTCJCEhYXRIhUlCQsLokMIIJCQkjA5pxCQhIWF0/Ot+TE1nLTGY1q354nUgdVYsNpgmwKNpM4Cqux7ktdB3DKq7p9MqAFbc6WkwzWmNxRsuf7mnu1WFLkxqeBqouuuNUqNrGlTXzu0pAIcflW8DoisBdcS79qvqeraq6l8t//zCoLpXAr8u9W/SiElCQsLokAqThISE0SEVJgkJCaNDKkwSEhJGh1SYJCQkjA6jSEkZ2dGL4R28EBB4mpTGnH1BJGflADC8QwuGtm2GqUrJrah4Zu8NokCjm5/QoIZNeKt1WwQBctRq/ns6mOvxYoLDB+296Ve/ERpB4EZCHLOOB5Gvo25VpE2AaPUate42ZjUscexTGwBtvobozXfJeZQOApjXtcFtTCPkJgqddUsi/GwOQUtTmLrDTS8dgMsHs7lyKAeZDOxcFfhPs8HSTv/fvKpIM3nwUMbiH03IzJIhlwvM/CifJgZIirl5TuDABi3qAnCrAyOmyzG31M+uFwyfllOVfWxui8E8yIxn00PRr0mOjJnN+tHG3gOAkIR7LLl9WCet//mIqambExO6tGHkmh0MXLaZiKQU3vfrBEAvz/qM8m7Jm+v3MGDZJkxVSsZ1bq2Tbl27anzWuSvj/thL4M7NrLh4ltV9xNOzHWq4079BY/rt3ELv7RuxUpkwvkUrnXSrIm0CIDc6i8cL/yHtQlyx5QkHHoNWoP7cDtT/ugNCgZaEPx/rpFkaqdFqzmxIAwNcWhv7oIAL+7IZtagaE1Y6YOeqIGRLZvkrlkNVpJnk5sL7n5gxengBm9fmMnGMmi/n6W/jmZEqsPUHLRO/kDN7nQIHVxkHNuj/4Ro6Laeq+lgdq+r83GECfm7Nii3v594SD0tHhp1awRunV9LW3gM/F0+dNMv9+Tl9+jReXl7Y2Njw+++/c+3aNTw9PQuDCfTlVnQ8fRb/KqakKBU42VgRlSKa0g1s1YRfQy6TlpMHwH/3H0Ol0K2W5mk0zAw+QkL2swSH+FiqW4gJDgqZHFOFAjOlEq2gxVSpJE/H0VJVpE0AJB+LxK6LKyqH4r+GFo3sMHE0RyYXf33NalmTF1X5Hb8gT8uRH5LxmWjLkR9eTeGoKC71VUz6yQGFUoY6XyAzWYuts36jOaiaNJNzFxXUcNPSuaPoZ921swY3V/191O9cFqjVEJxqiO3yCZTx3Ttahr0rICvBt15XDJ2WU1V97I3a7dkfeZnYnNRiy+UyGeYKFSZyJTKZDKVcWWjnWx5lFqZ58+Zx+/ZtlixZwtKlS7l+/Tq+vr4EBQVx+/ZtZs+erXPjy0Kt1eLbpB5zh/iRr9aw4uguADwcq+FgGctP4wfjZG3JpYgoFv91WifNqIx0ojJeSHDw6c6xR+EUaLWERj4h5GkEZ8ZNpkCr4WFqCttuXNVNtwrSJgDcxogum1m3ihcL62YOhY/zE3NICnqC27gmOmmWxPFVqTTrbYmjh/6Zcs9RKGXcD8vj8PJ0FCoZPqMsy1+pHKoizeTJUxkO9gLfLDThfrgMayuYNiW//BXLITUBqlUvapdddcjNFv+Z6/FRGDotp6r62IKbYshHB4fijp1/PP0HP5dmHOn1CQqZnLCEB5yKv6uTZpnDjzNnzrBx40aqV6/OiRMnWL16NSNHjmTlypWcOXNG54brwrHb4XSet4aVx8L4ecIQZDIxacG7QW0+2v4nr6/ahq25GR/4d66QrrlSycre/ahta8enwUcAGNakGe42tnTYsIYO63/iaXoan/t0r5huFaRNlEXO43QefXsJe9+a2LTULcLqZa4dykQul9G0l/6F42UaeJsybVt1Oo20ZNecVAStYW7BNGSaiVoDoecUDOqnZuNPeQwbrGb6TDPy9axNpd1tKtdz4FiVaTklYYg+9iJTGvYgJT+LnkHfEXBsEbYm5oyp00mndcscMZmZmZGUlISTkxMODg5kZ2djYmJCTk4OSmXlJyGn9fKmZ+NnKSkJyew4e5XLEWJKyt5LN/lykC82ZmbEZ2Rx7OYDsvLEnnPgyh3e6Vn6rQzT23eiV516ABx9FM7OW9f5pd8gHiQnM2LfrsIEh4B69dl/7zZZBeIh1/ab1/hv19KTeKsibSJuXzgZ/yQCYN3KEefB9Urdfuq5WGI238V1VCPsvHVLIn7O2a3pPLogTiDLlaDOE9j+YTyaAgF1vvi4/xcOWDnovheFbMnkwXnxO3GsqaBlX3PcPU0AaN7LjKBVGeRmCpjbVGx0U1VpJs9xdBDwqCXQrKm4Y3fz0TD/e4iKkVGnduULabXq8PhO0fppiWBhBaZmFW+vIdNyqqqPvd2wJ92dRN/wE/F3WH2v5Lw4X5emLLj5J2pBQ6Zaw4HIK/Ry8WTzo9ASX/8iZVaXadOmMXToUAIDA6lbty5jxozB29ubkJAQJk2aVK54aaw4GsaKo2EAtPGowffD+zJk+RZSs3Pp17Ix9+OSSMvJ5ciN+/Ru3pBdF66Tp9bg27Qe1yPjStVdcj6UJefFN21rasaBN0az+/ZNfrwQVux1N+Pj6V23Afvu3EIjCPSu24ArsaVPKlZF2oTz4HpldpTnpF2II2brPTxmtMK8TsVTAzqOsqHjqFfXS49Ts+39eEYsrXiErs9oK3yeZZ1G3sznwKJ0xi2zx8JWzq2TuTjWUmJeCWP/qkozeU6nDhp+XG3C7bsymjQS+OeqHJkM3Fz1G901biPj97UC8VECTjVkhPwp0Ny7ckXUkGk5VdXHVt8LLrUYvcjttBj8XZtxMekRSpmcbs6NuZb6VKe2l1mYevbsSYMGDTh69CgRERG0bNkSS0tLFixYUBgdri+XHkfx04nzbJw8DI1GS3xGFu9tEW9G3H72KrbmZuyeNgq5TMat6HgWHjpVjqLI6GZeuFlZE1C3PgF16xcuH7V/FysvnWO2T3eCRo0nT6PhTmICc06WHgv1IlWRNlEWcbvDQYCoDbcLl1k0sMVtTOMy1vr3cPc0oePrFuyclYJMIcPKXs6gz2311q2KNBMHe1j4dR6LlpqQkytDpYIFc/MwNdGvrdZ2MkZ+JGf9N1o0anB0hdGf6H/CuyrSckqiqvrY97cOMdOzH/u6vY9WEDiXGM6v4brNEf/rfkySu4DkLgCSuwBI7gKSu4CEhMT/KaTCJCEhYXRIhUlCQsLokAqThISE0SEVJgkJCaNDSkmRkJAwOqQRk4SEhNHxr/sxtZ5iuOuYLv8kXsfU+L+G0wS486Wo67F2kUF1H0/+BKi660EmXJhgMM0N7TYA8M7l0QbTBFjVegsA068MN6jukpY7AMNebwRF1xxpYxsaTFPucg+AC088DKYJ0K7WY8Cw17NB0TVthuxfUNTHSkIaMUlISBgdUmGSkJAwOqTCJCEhYXRIhUlCQsLoMIowgud096rH3AkBdP1QvGlQJoP3B3fBp3kdtILA0/gUvtlyjNTMHJ30RrXzYnjbZyEHyWl8cSCI5Ozi6/74ej8SMrL4+q/jOrdzUP2mTGnRDgHIURfwVegxrieKdiy9PRrybssOmCgURGWm89GJQ6Tm5eqsDYY1dQfRAfLBzw+wcLegRmANcZlW4PGWx6ReT0XQCLgFuuHiWzGvp+c8/TuOyCPxyE1kWLqZ02hibVRW+net2POp3N8VAzIZKksFzafUwtJFP4/uqgoNCDoFKzaAXA42VvD1f6BWDf00Q47KObRL9MkyNYMx76ipa4DghBcxVCiFofuY0YyYajrZMX1oV+QveCQP7NSMJrWdGDlvK2/M3czT+DQ+GtpVJz1PVycmdmrDiPU7GLB6MxHJKXzQs7h73pud2tK2gr2nrm01ZnXoxti/dtN370aW/xPGGr9BADR3dGZuZ1/ePrqfgD2/8igthU/a6X53flWYumdHZXPz25sknUsqtjz2WCw5cTm0XNCSFl+3IObvGDLCSza1K4vkm+lEHIih1eeN6LCgGQ6tbLm99nGFdV5Gk6/l6ooIWs+oS5eFjXFua8utXyP10qyq0IDcPPh0Hvz4NexbBz06w7wf9dOMfipj+1oln8wvYP5PBQwcqWHZfw1nhwyGC6Woij5WZmH65ptvSEtLK+slBsFMpeSbib1ZvKt41NHDmCSW7j5NgVoMCrgVEYerg25maTdj4glY/iuZefmYKBQ4W1uRml00cung4U6X+rXZcalifkn5Gg2fnjpMQs6zkIOEOKqbiyEHg+s3Zefd60Rmil7jSy+dYc3V8zprPzd1D4q+UWz5i6buKrmyQqbusUdjcerqhEMHh2LLky8l49TVCZlChtJSiUNHBxLOJOjc1udkPMrCvpkNZg6iqZFTu2okXk5Fq9bP/lXQCgiCgDpb/O7VuVrkevh9Q8mhAReDxe3og0Yj2utmPgtxyc5Bb48nlUpg0kdqqj372uo01JKaAuoC/XSf82Iohb5URR8rc7z9+++/c/r0aWbMmIG/v3/lW14On4/uxd5T17kflVhs+bWHRa6S1hamTA7swJ5TuhcStVaLb6N6fDNADDn48YQYcuBkZcms3t2ZtGUfb7RpXqG2RmamFxYegNkde3A04gEFWi11bO25k5zAWr9BuFvbcic5ka/Plu/095yqMHWvO07USrtZ/AcmPykfU4eiwyJTe1Oyn2br3Nbn2NSz4unf8eQk5GFe3ZTok4kIaoGCDDWm1Sq/dyrNFDSbVJOwL+6hslYgaMF7rn7XElVVaIClBXz5EYx4F+xsQKuFrSv0airVXaC6i1jcBQG2/qSktbcWpYEGTYYMpaiKPlbmiMnd3Z2VK1eyadMmhg0bxqFDh8jNrdh8SXkM69YCjVbL/tCbpbfD0ZZfPn6dKw+i2XlCtzST5xy7G473ojWsOBnGL6OHoFIoWDy0L/P/PklCZuVzysyVKlb6DsDDxo6Zp8X5HqVcjm/teswKCaLv3o0k5GSxoEtApbfxHH1M3UujpFHC8/ieilCtiTV1X3Pj+g8POD/rphjTY6VArtRvliD9SQ4P9sTSZXETfNc0p/5gZy4vfqjX6KaqQgPuhcPqTXBwI5zaC1NGwwdzSt9eRcjNgeVfK4mLkjHpI91GyeVRlaEUL6JPHytzxCSTyahfvz5btmwhNDSUnTt3Mm/ePDw8PHBxcWHx4so5R07t7003L7HKFqi1mJko2T57FCqFAtNnj99b/juJaVm0bejOgsmBbDxykc1Bl8rUfa+7Nz0bibrhCclsu3CVy0/FkIM9/9zkq0Bfmrk5425ny8wAca7K0coShUyGiVLBFweOlqg7vU1n/GqL9rxBEQ/Yeeca6wKG8CA1ieF/7iwMOYjPzuRuSmLhYd6ue9fZ1veNUttbFabuT3Y/IflyMgD2re2pNbRWiZqmDqbkpxRZ1uan5GNir9sIJ3xXFImXxAigak1tcPd3wq2HGCuUl1pA+K5IlFYV39vv/RZD3EXxVzcrJhfXjtUKJ7trB1Tn1sYoCjI0mNhUbmLdkKEBLxJyAVo1K5rsHjkYFqyE1DSoZld53cR4+OELFW61BD7/vgATPeb9DRlK8W/0sTK/4RcrXqdOnejUqRMFBQXcvXuXp091MxUviTUHwlhzIOyV5a4ONuyaM4YR32wFoEVdVxa/3Z/PfjlE6M2IcnWXnwhj+YlnIQe1avDD0L4MWrOF1Jxc+jdvzP34JP55Gk2Ppb8UrjOtW0eqWZiXeVZuyaUzLLkkniGzNTXj4OAx7L53k2WXixeGQ4/uMatDd1b8E0ZqXi69PRpyLTG2VN2qMHWvNbRWqR3lRezb2BN/Kh771vZocjUknk2k7oS65a4HUG9YDeoNE/fCrOgcLn99h46LmqO0UPB4XzQunRwqFfTY8HVXGr7uCkDijQyur3lCXmoBpnYqYi+kYeFkUumiBIYNDXiRpg1h6z5ITAZHezgWAu6u+hWlzHSYN8OELv4ahozRLYy1LAwZSvFv9LEyv+VRo0a9skylUtGsWTOaNWtWwhqGZWp/b5DJeG+wD+8N9gEgKjGdj9ccKHfdS0+iWHP6PJvGD0OjFUMO3t35h95tGt2kJW6WNgR4NCDAo0Hh8pF/7uTYk3BcLa3Z2W84cpmMqMx0/nPqb723qY+pe2m4+LqQG5fLlVlXENQCzj2dsW1S8YlQSzdzag9w5cIXt0AA20ZWNJpQW6+2ATg2s6ZOfyfO/vc+cqUMlZWSNp/o1qlLo6pCAzq2honDYdwHoFKBrTWsmKef5rEDChLj4WKInIshRW38bFEB1hUPy/mfoE8fK7MwDRs2zCAN1JWYpHR8PlhZ+PydZXv10ttx8Ro7LpY9Wb7i5NkKaa68cpaVV0pfZ8vtK2y5XXJkk67Mubav2PO0ghw+u7JLL80GUxoUey5TyKgzpo5ems+pGeBMzQBng2i9iEdAdTwCKp88WxKe7WV4ttc/wvxlRg0W/xmKgaM0DByl/0ipPGyclUzdqd81TM8xZB8zmuuYJCQkJJ4jFSYJCQmjQypMEhISRodUmCQkJIwOqTBJSEgYHVJhkpCQMDqklBQJCQmjQxoxSUhIGB3/ulFcl8HfG0zr9L6PAWj5nmFTUq4sF1NSGnxrWN37nz1LX9n0nUF1H4/9FADvIzMNphnmvwCA7sc+NpgmwAlf8fv3OzHdoLpB3cXvau6NAQbVndNMvFvAkIkmz9NMTj42XPIKQDcPMX3lYaSrQXXruosuH9sedDCo7sj650r9mzRikpCQMDqkwiQhIWF0SIVJQkLC6JAKk4SEhNFhVCkpXdrX5/MP+tB71PLCZRPe6ERPn0ZotQJ3w+P4fvUR8gsqdtd1jxb1+Hp0AD7/WVW47E3/dvRr3xSFXM6hC7dZ85fuLgOj23gxspWYvvIkJY3P/xLTV0yVCr7y70lzVxfkMhlXo2P46kgweery2zuoTlOmeHZAQCBHrearC0e5nlTcy+mLtj3xsKnGm8F7dH/zz5jtOYyHmbFsiyiySxlSsyMDarTDVK7iTnoU82/upkCo2Gc7s8kbPMqKZeeTIr/237t8RWJekc3qjogTHI37p0K6nzQewaOsGHY/PfHK3770nEBSfhor7lfOfSIlIo+L6+IoyNYik0P7KS441DOrlNZzqirRRBDg18VQozb4G9jsIzREyfffmbP3QMVDKErjWnA+oXtFQzqVqYw+U81xa1DxMmM0IyZ3VzveGd+tmMFYS8+a+Po04s0Zmxn3wa9YmpvwWmDrCunWqm7H9EFdkb9g6enT1AO/lg0ZsXArQ+dvom2Dmvi30u0MiaeLE2+2b8Prm3cQ+MtmHqekML2raHP7TqcOKORy+q/bTL91mzFTKZnq3b5czbo29sxq04Oxx36j78FfWX49lDXdi3toBNZuzKC6uiWjvEhty+osbzsZX5fi3ubdnDwZVrMT71/8hZGhSzBVqBhe20dn3VoWTvzQairdnb2KLa9pUZ0MdQ6Tzi8p/FeRolTLwomFXu/QtXrLEv/+es2eNLOtvC+TOk9L8NynNB1oT9/vPWg21IHQpTHlr1gGVZVoEvMEfvgULp7SW+oVoiLl/PKTGYJ+mRHFSIzUELQ+h1FzrZi6woauw83YOa9y9tXlFqawsDD++UfsWOvXr2fq1KmsWLGC/Pz8ctbUHVMTJV98GMiKDSeKLVfIZZiYKDE1UaJUyjExUZKfr7vvsZlKybyxvVm8t3j6Sk+v+vx16Q65+Wry1Rr+OHeTvu0a66R5MzYev5+K0ldcrK1IyRF/IS48jWRV6DkEQCsI3IpNwM22fFevfI2aT8P+KkpeSYqlupmYvAJQz9aBKZ7t+fHaGZ3f+3OG1vTmz6iLHIu9Xmx5H7fWbIs4Tbo6BwGBhbf28XeM7gVksHtn/oo5z4m44h7snrYeaAUtS1pPZV37jxhbxw85urtEDqjhw5HY85xKeNXTysuuPm3tG3Mw+lVLYV2JuZKFtYsJNdpYAeDezgqfGfqdXq+qRJPjf0Anf2irW2KZzuTmwqJvzXnrbcP69ytV0P99C6ztxX7r1kBBZoqApqDiI8cyx1gLFy7k4sWLqNVq3N3dkclkjBgxguDgYObOncs333xTuXfwEp+87cf+I1cJf1w82uXS9SdcvBrB7p/fQq3W8CQ6hf1HdA8jmD28F7vPXOd+dPH0FWc7a87dfVL4PC41E2c7K5111VotvRrUY15fP/I1GpaeEk3cQh4VabrZWDOuXSu++LtkH/EXicxKJzLrheSVtj05Gikmr1goVSzpHMjHoYdo4VDxQMrFd8TrcNra1y+2vJaFI7dNLFnSegKOpjZcTXnMivuHdNZddk80s2tTrbg5mEIm52LyPdbcP4ipQsW3Xm+Src5l91PdHDefH561eknXwcSGd+oP5rNrPxHo5q1zO18mPSYfMzsFZ1fGkhKRi4mFglZj9TOjq6pEk5HTxP/v6Oc7+ArLl5jRp18+deoa1ojOzlmBnbN4OCsIAofX5tCogwpFJWK3yhwxnT59mh07drBt2zbOnTvH999/T7du3fjqq6+4dq1ieWylMah3SzQagUPHbrzyt76+zXB1smXQxNUMmriGmLg0po3vrpPu612epa+cfTV9RV6CH7VWW7GqfvR+OB2WrWH56TA2vDGk2JjA08WJ7aNfZ8ulKxx/8EhnTXOlipVdB+JhXY2ZoX8BsLBTHzbeucy91MRy1q4YSpmCdg4N+PzqNiacXYG1ypyp9fVPdPkz+hzL7+2nQNCQqc5l15NT+FTXz4ZZIZMzq+lYVj/YR3J+evkrlIGghujLWdT3s6VWhYhQAAAXe0lEQVTPQg8a9rXj+LxINAX6H9NURaKJoTm4X4VCAQF9DBRQVwL5uQK7v80mOUbLgPfNK6VRbhhBRkYG2dnZ5ObmkpmZSbVq1cjNzaWgoPJv7M0Rnencrh4ABWoNZqYq1v8wFpVSTElZ/8NYPvl6D906NiDo1G1ycsVt/XHkGtPf8i1V9+2+3nRv/ix9RaPFTKVk56ejUCoUmD57PG3178SkpFPdtii6xsnWirjUzFJ1P+jijW8DUfdBYjJbL1/lUqSYvrL72k3m9vbF1tyM1JxcAps05KsAX+YeCebArdKz36Z7+eBX81nyytMH7HxwlXU9hvIgLYnhR7aTp1HjYmFNOyd36trY82bTttiammOtMmVDz6FMCN5dou7ken74VG8CQEjCbdaGB5X4usS8dE7G3yRbkwfA4Zh/mFiv9M92Qt0AOjs2BeBM4i02PCw5otzPpTXhmTE8zBTnbWQyGZoyJtTHefTG21EsXGGJN9j4+FWP9IbWNXExs2dqfTHxuJqJNXLkmMhV/HB3Z6naz7m6PZGoi+L3m5ehwaaGCY4NxR2mZntrzq2OIzOuAFv3yseQGDLRpCoJOmxCXh68+5YlBQWQny8+njs/GwfHih9yHd+cw91z4v7ZqIOK1gGmbJ+biWNNBeO+tUJlWrmwhzIL0+TJk/H390cQBD755BMmTpyIt7c3YWFhvPbaa5XaIMC67WdYt/3V+RKX6jZsXDaeiR9tAuDew3i6dmzA4RM30WgFunVswK270aXqrj4UxupDr6avuNnbsHvWGN74TkxfOXH9IVP6dGT3metotAIDOjTlj3O3StVddjqMZadF3bbuNVg6sC8D1m8hJSeXAZ6NuZeQRGpOLr0bNeALvx5M2LGXG7FxZX4GS66GsORqCAC2JmYcDBzH7vAbLHthHik2O4MOu4vOJA6t14w+tRuVeVZubXhQqcXoRYLjbuDr0pw/Is+Tp1XT1cmT22mlR3BveHi41GL0InUsXenm1II51zailCsY7N6ZoNjLpb5+4+O/SyxGL3I7PYJRZ+cWPh/jEYCtylLns3JeIxzxGuEIQE6KmoMfPiIpPBeHembE3RQDGK2cKn/cZehEk6pk2aqiyei4WBlT37Ri5c+Vz1fsMcacHmPEIp+ToeXnDzLx6mVC95H6neUsszANHDiQgIAANBoNlpaWtGvXjpCQED7++GM6d+6s14Z1YfPus0yb0IPNyyeQX6Ah/HECP/xc/pxNeZy68ZAGbg5s/XgkKoWc49fDOXC+9ML0Ihcjo1gdep4to8T0lbjMLN7ZI87jzOjeGZkM5vftVfj6S5HR/PdI6dFQAKMbtRKTV2o1IKDWC8krQTtIzTPsBOVz9j4Nw0ZlzoaO7yGXybmbHsWPd//UW3fjoyN80Ggw6zt+jFIm50T8Nf6MLv2eqH8b82pKuv6nBhfWxqHO1aJQyej6HzcUJpU/Qf3/Q6KJIbjwZz5pCVruhOZzJ7To5NjY+VZY2FTs8/3XbU+km3ilm3hBuokXpJt4pZt4JSQk/k8hFSYJCQmjQypMEhISRodUmCQkJIwOqTBJSEgYHVJhkpCQMDqklBQJCQmjQxoxSUhIGB3/ulGcv/fXBtM6EvYFAN4jFhtMEyBs+wwAms8w7AWW1xeLFxXW++EHg+qGf/SRqLtjvuE0h88Cqu5CyMZ755bzyopxZ8gcAF4Lfceguns6ibcErbjT02Ca0xoHA/DLvS4G0wSY1FB0cKiqCyGrqr0lIY2YJCQkjA6pMElISBgdUmGSkJAwOqTCJCEhYXQYRUrKW+/1omvPpmSk5wDw9EkS878o7rUz9QM/3GraM+fj8o3BXqZr2/rMebs3vd5cUbise/sGjBvYAZVKQWxiOnNX/UV6pm4WIyM6e/F6Jy8EQSAyKY2vdgWRnJmDXCZj1pAetK3rDsDpO49YfEA3S1mAMS1bMrJFCwAiUlP5PCiIpBxR9/Nu3eji4YFSLmftxYts19FBdGBtTyY37ogA5GoKmHvpCNdTYpEB//HqQXe3+giCwOOMFGZf/IvkvGyd2wuGTzPpX7M5bzbwLmzvvKt/cyNVvLt9SiMfBtZqgUIm58DT66y4fbJssRcQBIGodbcxq2GJY5/aAGjzNURvvkvOo3QQwLyuDW5jGiE3UeisWxLhZ3MIWprC1B1ueukAXD6YzZVDOchkYOeqwH+aDZZ2+o8nDJVm8jKGaq9RjJiaNq/J/Dl7eXvcWt4et/aVotTVtyk9A5qXsnbZuLvY8d6orsheSElpXNeZGeN78tmSPxj9n408jUlh6hu6JYQ0dXdiXPc2jFm+gyHfbyYiMYVpvcWUlP5tm+BRvRpDvt/M0MVbaFvXHf8WDcpRFGnm5MSkNm0YtmMHfTZtIiI1lenPPK9GtGiBR7Vq9Nm4kUFbtzKhdWtauJTv/13H2p6ZLX2ZcHIH/Q+vY+XNM6zyEQ3+htX1opm9KwMPr6fv378QkZnMrJalO1i+TFWkmdSxcuCT5r2YHLqNwcE/s/rOaX7s+DoAXZ3rE1CjCa8Fr6X/0dV0cPSgd42mOunmRmfxeOE/pF0obt6XcOAxaAXqz+1A/a87IBRoSfjzcYXa/DKp0WrObEgDA1wdGPuggAv7shm1qBoTVjpg56ogZEvpTqu6Ysg0k6pqb7kl8ujRoxw9epSEhARUKhW1atWiT58+tGrVqlIbfBmVSkH9hi4MHemNa41qREcms2bZERLiRG/nmrUdeX2UN1vXn6ZNx4p1dFMTJV+925dlW07y32l9C5f39mnCgRM3iE0Ut/HL7lBsrXXzJr4VGU+/b39FrdViolTgZGtFVJKoI5fJMDdRYaJUIJPJUCkVOmXKAdyIj8d3wwZRV6HA2cqKp2liNpt//frsuHYNjSCQnpfHwbt3GdSkCddiY8vUzNdq+Oz8nyTkPktfSY7B0cwKlVzO/bREFlw5Rr5W8+xvsYxu0EantkJRmklCXsorf3sxzcRapbvnc75WzReXD5KQK3bmG6nRYntlcnq5Nebg0xvkaEQb170RVxhQqzl/R5Vv8Jd8LBK7Lq6oHIq7Klo0ssPE0bzwR8usljV5UZXf8QvytBz5IRmfibYc+eHVz6WiuNRXMeknBxRKGep8gcxkLbbO+o3moOw0k8oEB1RFe8scMf3000/s2bOHFi1aIJPJaNmyJW5ubsyePZvffvutUht8GQdHa65cesT61cG8PfZn7tyM4r8L3wDAzFzFp18O5Ptv/iA7O6/C2p9O8uP3Y9d4EFE8faWmSzUUcjnfzRjIpgVj+HiiL9k5usdRqbVaejarx9E5k2lT153fL4iBB/sv3CI9J49jcyZz/Mu3eJKYyslbDyuk61evHmfeeot27u7suSnqulpbE5NZtMPEZmTgYlV+qktUVhonYsILn89q1Ytj0fcp0Gr5JymKmyniCMJGZcZ7nj789fS2zm1dcX8vR+MuvrL8eZrJgttb0FIxg/+o7DROxt4vfD6zuT/HY+5SIGhxtbAhNqcoiCA2Jx1nc93sId3GNKJap1fN06ybOWDqYgFAfmIOSUFPsGnnXKE2v8jxVak0622Jo4f+mXLPUShl3A/LY834RCJvFNC8l36WtSCmmTRsL7ZR3zSTlzFUe8ssTIcOHWLVqlWMHDmSlStXEhoayvjx49mxYwcbNmyo1AZfJjYmldkzdhD5JAmAXVvDcK1RDRdXOz6a1Z/9uy/w+GFCOSqvMsTPC41Gy8ETr6avKJUKfNrU5btfghj32WaSU7OZOdm/QvrBN8LpOmcNqw+H8dNbQ5DJ4G3/jqRk5tDtq5/o9fVabC3MGNutYgGdQeHhtFu9mh/DwtgwRExfKSnVRVOBO4nMFSqWdxpMbatqfHa+uH1uLSs7dviO5mLiUzbfv1Shtr6ModJMzBUqlrYfSi0re2ZfPgBQYjad1kBpjTmP03n07SXsfWti09KxUhrXDmUil8to2suy/BdXkAbepkzbVp1OIy3ZNScVoYKJPqVhiDSTkjBEe8s8lMvLyyMnJwcLCwtyc3NJTU0FwMLCArm88tNTYyd3w9tHtBWNiU7lzMk7HPu7KJBRBqjVGpp71aJmLQeGvNEBaxtzLK1M+WbxcGbP2FGi7uShnfBpI6avqNUaTE1VbPx2TGH6ysZvxzBj4V4SUzIJf5JAcpo40Xvw5A1WzC49f/ndAG+6e4qHkQ/jktkRepV/HomhCPvO3+SLob7YmJvh27w+3+47jlqjJVOTzx8Xb+HXogGbTpZsxv9hp0741n2WvpKczJYrV7gULeruunGDr319sTUzIzo9HSfLog7vbGVFbEbJsc4fNuuKbw1xXutY1H12PrzC2i7DCE9PYtTxreRpimKFOjrV5sdOg/j59ll+uVu2L3dVpZm816Q7PV3FvhAcc4/djy+z2ns44RmJjDu1iTyt2N7onDSqmxWNEp3NrYnNKfkziNsXTsY/YtyVdStHnAfXK/V9pZ6LJWbzXVxHNcLOu2K5fWe3pvPogjiBLFeCOk9g+4fxaAoE1Pni4/5fOGDloPvhTMiWTB6cF0fvjjUVtOxrjrunCQDNe5kRtCqD3EwBc5uKjW6qKs2kqtpbZmEaMmQII0aMwMfHh5CQEIYMGUJUVBTvvvsu/fr1q9QbAdi09iSb1opnVDzqVmfx6nHcvPqU2JhU+g9pw6PweBITMhgxYGnhOn59W9ClZ5Myz8qt3R3K2t2vprS6ONqwddE4xn22GYDj5+4xbVQ3fv39HOmZuXRv14Db4aXP16w8HMbKw2JKSus6NVg0ui9Df9hCalYuga0b8yA2ibTsXG5HxRPQsiEXwiNRyuV0b1qXaxGlx08vDQ1laajY3rY1arAsMJD+mzeTkpvLwMaNuZeURGpuLkfDwxnq6cmx8HAsTEzo16gRXxw7VrLmjVMsvSFmStuamLHffyJ7Hl1j+c2QYq9r7VCD1T6v8UHo75yKLf9ws6rSTJbfPsHy2yfE9qrM2NNzMvsirrLyTvFc7ODoe7zbpCu/PbqMWtAyuHZL9kWUnATpPLhemcXoOWkX4ojZeg+PGa0wr1Px1ICOo2zoOOrV9dLj1Gx7P54RS50qrOkz2gqf0eLjyJv5HFiUzrhl9ljYyrl1MhfHWkrMK2jsD1WXZlJV7S2zML311ls0b96cW7duMXPmTLy9vcnKyuK7776jUaNGlXojL/P4YQKrfjjM3EVvIFfISYxPZ/4c3U8vV4aQyw+pbm/NqjlvIJfJiE1MZ/7P5UcTAVx+FMXPR8+z/m0xJSUhPYsPNoiG9Qv3n2TW4B788ek4NFqBc/efsD741XmYkrgYFcWqc+fY9vrrqLVa4rOymLp/PwBbr16llp0dB8eMwUShYPu1a5yPLD1q6Tmj6rfGzcIGf/dG+LsXfV9jjm/jg+ZdkQGfePXgE68eAERmpfJ2SOnRUFXN8LptcbWwpZdbY3q5FUW2TwjZzPHYezS0deK3Hm+ikisIjr7H70/0C12N2x0OAkRtKJpbs2hgi9sY3eLiqxp3TxM6vm7BzlkpyBQyrOzlDPrcVm9dQ6aZVFV7/3XbE+kmXukmXpBu4gXpJl7pJl4JCYn/U0iFSUJCwuiQCpOEhITRIRUmCQkJo0MqTBISEkaHFEYgISFhdEgjJgkJCaNDKkwSEhJGh1SYJCQkjA6pMElISBgdUmGSkJAwOqTCJCEhYXRIhUlCQsLokAqThISE0SEVJgkJCaNDKkwSEhJGh9EVpgMHDtC3b1/8/PzYunWrwXQzMzPp168fkTo4P+rKihUrCAwMJDAwkIULFxpMd9myZfTt25fAwECDhT68yHfffcfMmTMNpjd27FgCAwMZOHAgAwcO5OrVq3prBgcHM2TIEHr37s0333xjgFbCrl27Cts4cOBA2rRpw9y5hjGs279/f2Ff+O677wyiCfDzzz8TEBBA//79Wb16tV5aL+8DoaGh9O/fH39/f5YsqbwpYkn71qeffsrevXo40QpGRGxsrNCjRw8hJSVFyMrKEvr37y/cv39fb90rV64I/fr1Ezw9PYWnT58aoKWCcObMGeGNN94Q8vLyhPz8fGHs2LHCkSNH9NY9d+6cMHz4cKGgoEDIyckRevToIYSHhxugxSKhoaFChw4dhE8//dQgelqtVujcubNQUFBgED1BEIQnT54IPj4+QkxMjJCfny+MGDFCOHHihMH0BUEQ7t27J/j5+QlJSUl6a2VnZwvt2rUTkpKShIKCAmHo0KHCmTNn9NY9c+aM0K9fPyEjI0NQq9XClClThMOHD1dK6+V9ICcnR+jWrZvw5MkToaCgQJg4cWKlPuOXdWNjY4UpU6YILVq0EPbs2VOptgqCIBjViCk0NJSOHTtiZ2eHhYUFAQEB/P132Qb4uvDbb7/x5Zdf4uRUcXP40qhevTozZ87ExMQElUpFvXr1iH6WcKIP7du3Z9OmTSiVSpKSktBoNFhYWBigxZCamsqSJUuYOnWqQfQAHj58iEwmY/LkyQwYMIAtW7borRkUFETfvn1xcXFBpVKxZMkSvLy8DNDaIr766iumT5+Ovb293loajQatVktOTg5qtRq1Wo2pqaneurdu3cLHxwcrKysUCgVdunTh6NGjldJ6eR+4du0atWvXpmbNmiiVSvr371+pfe1l3QMHDuDr60ufPn0q1c7n6B9WbkDi4+OpXr164XMnJyeuXdPPcB5g3rx5emu8TIMGRdHfjx8/5tChQ+zYUXKsVEVRqVT8+OOPrF+/nt69e+PsXPkQxheZM2cO06dPJyam9OSWipKeno63tzdfffUVubm5jB07ljp16tD5Wbx5ZYiIiEClUvHmm2+SkJBAjx49+PDDDw3W5tDQUHJzc/XeeZ5jZWXFBx98QJ8+fTAzM6N9+/a0bl2xPMGS8PT0ZP78+UyZMgVzc3OCg4MRKmkG8vI+UNK+FhcX9/JqFdadNGkSAJcu6ZdRaFQjppI+dFkJYY/GxP3795k4cSKffvopHh4eBtN9//33CQsLIyYmxiCpx7t27cLV1RVvb28DtK6IVq1asXDhQiwsLLC3t2fo0KGcPHlSL02NRkNYWBiLFi3it99+4/r16+zbt89ALYYdO3YwYcIEg+nduXOHPXv2cPz4cUJCQpDL5axbt05vXW9vb4YMGcKYMWOYNGkSbdq0QaUyTMqvse9rRlWYnJ2dSUxMLHweHx9v0MMvQ3Pp0iXGjx/PjBkzGDx4sEE0w8PDuX1bjBMyNzfH39+fu3fv6q176NAhzpw5w8CBA/nxxx8JDg5m/nz9U1UuXrxIWFhY4XNBEFAq9RuIOzo64u3tjb29PWZmZvj6+hpk5AyQn5/PhQsX6NnTcKknISEheHt74+DggImJCUOGDOH8+fN662ZmZuLn58eBAwfYvHkz5ubm1KxZ0wAtNv59zagKU6dOnQgLCyM5OZmcnByOHDlC165d/9fNKpGYmBjeffddvv/+ewIDAw2mGxkZyezZs8nPzyc/P59jx47Rpk0bvXU3bNjAwYMH2b9/P++//z49e/Zk1qxZeutmZGSwcOFC8vLyyMzMZN++ffj5+eml2aNHD0JCQkhPT0ej0XD69Gk8PT31bivA3bt38fDwMNi8HUDjxo0JDQ0lOzsbQRAIDg6mefPmeutGRkby7rvvolarycjIYNeuXQY7/PTy8uLRo0dERESg0Wg4ePCgUe1rRjXH5OzszPTp0xk7diwFBQUMHTqUFi1a/K+bVSLr1q0jLy+PBQsWFC4bPnw4I0aM0Eu3W7duXL16lUGDBqFQKPD39zdo4TM0PXr0KGyvVqtl5MiRtGrVSi9NLy8vJk2axMiRIykoKKBz58689tprBmnv06dPcXGpWBR4efj4+HDr1i2GDBmCSqWiefPmvPXWW3rrNm7cGH9/fwYMGIBGo2H8+PEG+ZECMDU1ZcGCBbz33nvk5eXRrVs3evfubRBtQyBZ60pISBgdRnUoJyEhIQFSYZKQkDBCpMIkISFhdEiFSUJCwuiQCpOEhITRIRUmCQkJo0MqTBISEkaHVJgkJCSMjv8HJ1rAN1AuR8cAAAAASUVORK5CYII=\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": [
    "For every element in the above matrix, as we were filling up, we recorded the traceback direction in a 3D matrix called `tracer`. Let's visualize the three channels of the `tracer` matrix: $\\nwarrow$, $\\uparrow$ and $\\leftarrow$ "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAvcAAADkCAYAAAAcogunAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAapklEQVR4nO3de3BU9f3/8ddCNkCAbxEkQSqWUjPRQgKYQl1AadBwCwiEBIECI+E6g6BUO1BqxWEAgSIUbadDZ8BOgCmGicAgacEYYQgJIs5oOkNNKcqtQoNUueaySc7vj9804iUbsrvZPedzno+/ksXznnfOntfJ209OPvFYlmUJAAAAgOO1inYDAAAAAMKD4R4AAAAwBMM9AAAAYAiGewAAAMAQDPcAAACAIRjuAQAAAEMw3AMAAACGYLgHAAAADMFwDwD4mg0bNmjDhg3RbgMAEISYaDcAALCPNWvWKDc3V5Lk9/u1ZMmSKHcEAGgOVu7R4NChQ0pKStKDDz6oH//4xyotLY12SwAiqLCwUCdOnNC0adM0bdo0HT9+XIWFhdFuC0AEWJYV7RYQJgz3aHD16lXFxcUpOTlZ7733nnw+X7RbAhBBjz/+uHJzc9WxY0d17NhR27Zt0+OPPx7ttgC0sLq6Oj333HM6cOBAtFtBGDDc42tGjBihfv36adasWbp+/Xq02wEQYXFxcd/5MQAz1dXV6fnnn9fNmzeVlpYW7XYQBjxzjwb33XefampqlJ2drdWrV2vWrFnasmWLOnbsGO3WAABAC/jrX/+qgoICSVJycvJ3/jcvv/yyMjMzI9kWQuCxeMgKjVi9erU+/PBDBnzAZV577TVJ0sKFC6PcCYCW9r+V+1u3bum1115TbGxstFtCiHgsB41atmwZj+gAAGCw1q1ba/369Wrfvr3efffdaLeDMOCxHAS0bNkybdmyRX6/P9qtAACAFtC6dWu98sor8ng80W4FYcBjOQAAAIAheCwHAAAAMATDPQAAAGAIhnsAAADAEAz3AAAAgCEY7gEAAABDMNwDAAAAhoj4PvdlpyqCPjYlMT7kGtRxTi92q2OnXm6v4yR2OHd2fR/tUMdOvZhah/wHz47n36Q6durF1DqRyj8r9wAAAIAhGO4BAAAAQzDcAwAAAIZguAcAAAAM0eQv1J4+fVoHDhzQpUuX1KpVK8XHx+uRRx5RcnJyJPoDEEXkH3Av8g84U8CV+x07dugXv/iFJCk5OVm9e/eWJP3mN7/R1q1bW747AFFD/gH3Iv+AcwVcuc/NzdWePXvUrl27r70+c+ZMTZgwQTk5OS3aHIDoIf+Ae5F/wLkCrtzHxMSotrb2W69XVVXJ6/W2WFMAoo/8A+5F/gHnCrhyP3/+fI0fP14+n09du3aVJF2+fFnHjh3T4sWLI9IggOgg/4B7kX/AuQIO92PHjtXAgQNVWlqqiooKWZaln/zkJ1q4cKESEhIi1SOAKCD/gHuRf8C5mtwtJyEhQePHj49ELwBshvwD7kX+AWdin3sAAADAEAz3AAAAgCEY7gEAAABDeCzLsqLdBAAAAIDQsXIPAAAAGKLJ3XLCrexURdDHpiTGS5L2e5NC6iHDX06dJmqE8j5JX71XJtWxUy+313ESO5w7u76Pdqhjp15MrUP+gxfO82+H77d2q2OnXkytE+5eGsPKPQAAAGAIhnsAAADAEAz3AAAAgCEY7gEAAABDMNwDAAAAhmC4BwAAAAzBcA8AAAAYIuA+95999lnAg7t37x7WZgDYB/kH3Iv8A84VcLifN2+ezpw5o/j4eFmW9bV/83g8euedd1q0OQDRQ/4B9yL/gHMFHO7/8pe/aOrUqVq+fLlSU1Mj1RMAGyD/gHuRf8C5Aj5z36FDB61cuVJ79uyJVD8AbIL8A+5F/gHnCrhyL0kpKSlKSUmJRC8AbIb8A+5F/gFnYrccAAAAwBAM9wAAAIAhGO4BAAAAQ3isb+5xBQAAAMCRWLkHAAAADNHkbjnhVnaqIuhjUxLjJUn7vUkh9ZDhL6dOEzVCeZ+kr94rk+rYqZfb6zhJuL5mrvGWqWOnXkytQ/6DZ8fzb4fv2+GqY6deTK0T7l4aw8o9AAAAYAiGewAAAMAQDPcAAACAIRjuAQAAAEMw3AMAAACGaHK4Lyws1LZt23Tu3Lmvvf7GG2+0WFMA7IH8A+5F/gFnCjjcr1+/Xtu3b9eZM2c0efJk7d27t+Hfdu7c2eLNAYge8g+4F/kHnCvgPveHDx/W7t27FRMTo+nTpysnJ0exsbEaNWqU+MO2gNnIP+Be5B9wroDDvWVZ8ng8kqSePXtq8+bNmjlzpjp37tzwOgAzkX/Avcg/4FwBH8sZOXKkpk+frrKyMklSYmKiNm3apGefffZbz+ABMAv5B9yL/APOFXDl/umnn1Zqaqrat2/f8FpqaqrefPNNbd26tcWbAxA95B9wL/IPOFfA4V6SfD7ft16755579Otf/7pFGgJgH+QfcC/yDzgT+9wDAAAAhmC4BwAAAAzBcA8AAAAYwmOxYS0AAABgBFbuAQAAAEM0uVtOuJWdqgj62JTEeEnSfm9SSD1k+Mup00SNUN4n6av3yqQ6durl9jpOEq6vORzXeLjyZrfrwcRr3KQ65D94djz/drinSeH9/m+HXkytE+5eGsPKPQAAAGAIhnsAAADAEAz3AAAAgCEY7gEAAABDMNwDAAAAhmhyt5wzZ86oXbt2SkhI0K5du1ReXq6HHnpIo0ePjkR/AKKI/APuRf4BZwo43P/5z3/Wtm3bVF9fr4cfflgXL15Uenq68vPz9emnn2rBggWR6hNAhJF/wL3IP+BcAYf7/Px8FRQU6PPPP9eYMWN07NgxtWnTRtnZ2crKyiLcgMHIP+Be5B9wroDP3NfX1ys2Nlbf//73lZOTozZt2jT8W11dXYs3ByB6yD/gXuQfcK6Aw/3w4cM1bdo01dXVaeHChZKkjz/+WFOnTtWoUaMi0iCA6CD/gHuRf8C5Aj6W88wzz+j9999X69atG16LjY3VwoULNXTo0BZvDkD0kH/Avcg/4FxN7pYzYMCAr33eq1cv9erVq8UaAmAf5B9wL/IPOBP73AMAAACGYLgHAAAADMFwDwAAABjCY1mWFe0mAAAAAISOlXsAAADAEE3ulhNuZacqgj42JTFekrTfmxRSDxn+cuo0USOU90n66r0yqY6derm9jpOE62sOxzVuh7zdXscO15Vdr3GT6pD/4Nnx/Nutjon3RpPqhLuXxrByDwAAABiC4R4AAAAwBMM9AAAAYAiGewAAAMAQzRru16xZ01J9ALA58g+4F/kHnKPR3XJ+9atffeu1oqIiXb16VZL08ssvt1xXAKKK/APuRf4BZ2t0uO/UqZP27Nmj+fPn6//+7/8kSceOHdPAgQMj1hyA6CD/gHuRf8DZGn0sZ8mSJdqwYYMKCgrUvXt3TZgwQd/73vc0YcIETZgwIZI9Aogw8g+4F/kHnC3gH7Hy+Xx68MEHtXz5ch06dEh1dXWR6gtAlJF/wL3IP+BcTf5CbadOnbRp0yb16tVLXbt2jURPAGyC/APuRf4BZwq4cn+77OxsZWdnt2QvAGyK/APuRf4BZ2GfewAAAMAQDPcAAACAIRjuAQAAAEN4LMuyot0EAAAAgNCxcg8AAAAY4o53ywmXslMVQR+bkhgfcg031NnvTQq6Roa/POQat9ex27kx8fpzEhOvcbvVMfEaN6kO+Q+eHc+/SXXsNItI9r3H2un7T2NYuQcAAAAMwXAPAAAAGILhHgAAADAEwz0AAABgCIZ7AAAAwBABd8spKytTSkqKJKm0tFSHDx9WTEyM0tPT1bdv34g0CCA6yD/gXuQfcK6AK/fLly+XJO3YsUOrV69Wt27ddPfdd+vFF1/U9u3bI9IggOgg/4B7kX/Aue5on/u8vDzl5ubqrrvukiRlZWUpKytL06ZNa9HmAEQf+Qfci/wDzhNw5b62tlb19fXq0qWL4uLiGl6PjY1Vq1Y8rg+YjPwD7kX+AecKmNC77rpLQ4cO1b/+9a+GH9GVlpZq8uTJGjlyZEQaBBAd5B9wL/IPOFfAx3Jyc3MlSZ988omuXbsm6f//X/uiRYv0s5/9rMWbAxA95B9wL/IPONcdPXPfq1evho9TU1NbrBkA9kP+Afci/4Dz8OAcAAAAYAiGewAAAMAQDPcAAACAITyWZVnRbgIAAABA6Fi5BwAAAAxxR7vlhFPZqYqgj01JjA+5hhvq7PcmBV0jw19um17s1k+4ewlXHSex0/sYrmvTbnXCcY+1Qy+392NSHfIfPDuef5Pq2KmX2+vY7R5rp+8/jWHlHgAAADAEwz0AAABgCIZ7AAAAwBAM9wAAAIAhGO4BAAAAQzQ53B85ckTXrl2TJO3Zs0crVqxQfn5+izcGIPrIP+Be5B9wpoDD/apVq7R582ZVV1frd7/7nfbt26f7779fb7/9tlauXBmpHgFEAfkH3Iv8A84VcJ/7o0ePat++fWrdurUOHTqkvLw8xcbG6sknn9SYMWMi1SOAKCD/gHuRf8C5Aq7ct23bVleuXJEkdenSRbdu3ZIkVVZWKiYm4n//CkAEkX/Avcg/4FwBE/r0008rKytLGRkZ6tWrl6ZPny6fz6fi4mLNnj07Uj0CiALyD7gX+QecK+BwP2zYMCUmJqqwsFBnz55Vv3791L59e61Zs0YpKSmR6hFAFJB/wL3IP+BcTf5srUePHpo5c2YkegFgM+QfcC/yDzgT+9wDAAAAhmC4BwAAAAzBcA8AAAAYwmNZlhXtJgAAAACEjpV7AAAAwBAR/0sUZacqgj42JTE+5BrUubMa+71JIfWS4S8PuRe79RPuXsJVx0nscD2E+/ybVMdOvdxexw73xnDVCXcvTmLi+Tepjp16aYk6drivRer7Pyv3AAAAgCEY7gEAAABDMNwDAAAAhmC4BwAAAAzBcA8AAAAYIuBwv3LlSl29ejVSvQCwEfIPuBf5B5wr4HC/Z88eTZo0SQcPHoxUPwBsgvwD7kX+AecKONzfe++9+sMf/qDc3FxlZ2eroKBAVVVVkeoNQBSRf8C9yD/gXAH/iJXH49H999+v7du3q6SkRG+88YZWrVqlnj17qlu3bnrllVci1SeACCP/gHuRf8C5Ag73lmU1fDxo0CANGjRIfr9f5eXlOn/+fIs3ByB6yD/gXuQfcK6Aw/3Pf/7zb73m9XrVp08f9enTp8WaAhB95B9wL/IPOFfAZ+6zs7Mj1QcAmyH/gHuRf8C52OceAAAAMATDPQAAAGAIhnsAAADAEB7r9l+JBwAAAOBYrNwDAAAAhgi4FWZLKDtVEfSxKYnxIdegzp3V2O9NCqmXDH95yL20RD+h1AlHjZao4yR2uB7s+j7aoY6demmJOna6x4arFycx8fybVMdOvZhaJ9wzTWNYuQcAAAAMwXAPAAAAGILhHgAAADAEwz0AAABgCIZ7AAAAwBBN7pZTWlqqtm3bqn///tq6dauOHz+uPn36aO7cuYqNjY1EjwCihPwD7kX+AWcKONyvW7dOJ06cUG1tre699155PB5NmTJFRUVFWrFihVauXBmpPgFEGPkH3Iv8A84VcLg/cuSI9u7dq5qaGg0dOlTFxcXyer169NFHNW7cuEj1CCAKyD/gXuQfcK6Az9xblqXr16/riy++UFVVlW7cuCFJqqqqkt/vj0iDAKKD/APuRf4B5wq4cj9nzhwNHz5clmXpl7/8pXJycuTz+VRaWqqJEydGqkcAUUD+Afci/4BzBRzux40bpxEjRqiurk7t27fXgAEDVFxcrOeff16DBw+OVI8AooD8A+5F/gHnanK3nLZt2zZ8nJSUpKSkpBZtCIB9kH/Avcg/4Ezscw8AAAAYguEeAAAAMATDPQAAAGAIj2VZVrSbAAAAABA6Vu4BAAAAQzS5W064lZ2qCPrYlMR4SdJ+b2i/sZ/hL6dOEzVCeZ+kr94ru9Wx07kJ1/vtJHY4d3bKW0vU4R7beB07XX/h6sVJ7PD9wK7fm+xQJ9z5t8PXZLc64e6lMazcAwAAAIZguAcAAAAMwXAPAAAAGILhHgAAADBEk79QW1hYqMLCQl2+fFler1f33XefRo0apf79+0eiPwBRRP4B9yL/gDMFXLnfvHmz8vPzlZKSIo/Ho379+ql79+564YUXlJeXF6keAUQB+Qfci/wDzhVw5b6goEB79uyRx+PRxIkTNWfOHOXm5mrixImaNGmSJk2aFKk+AUQY+Qfci/wDzhVw5b66ulqVlZWSpKqqKn355ZeSpLi4OLVqxeP6gMnIP+Be5B9wroAr95mZmZoyZYqGDBmi4uJiZWZm6t///rcWLFigMWPGRKpHAFFA/gH3Iv+AcwUc7ufOnavk5GSdPHlSS5culc/n082bN7V27VolJYX2F8wA2Bv5B9yL/APO1eRuOT6fTz6fr+Hz9u3bE2zAJcg/4F7kH3AmHpwDAAAADMFwDwAAABiC4R4AAAAwhMeyLCvaTQAAAAAIHSv3AAAAgCEY7gEAAABDMNwDAAAAhmC4BwAAAAzBcA8AAAAYguEeAAAAMATDPQAAAGAIhnsAAADAEAz3AAAAgCEY7gEAAABD2Ga437dvn0aPHq309HTt2LEj6Do3btzQmDFjdOHChaBr/P73v1dGRoYyMjK0bt26oOts2rRJo0ePVkZGhl5//fWg6/zP2rVrtXTp0qCPnzFjhjIyMjRu3DiNGzdOH330UVB1ioqKlJmZqZEjR2rlypXNPn7Xrl0NPYwbN06pqalasWJFUL3s3bu34b1au3ZtUDUk6U9/+pNGjBihsWPH6o9//GOzjv3mNVdSUqKxY8dq+PDh2rhxY0i1JGnJkiV68803m1XHacKVfyn0ewD5DyzU/EvhuweQfzOQ/6aR/28j/wFYNnDp0iUrLS3N+uKLL6ybN29aY8eOtU6dOtXsOh9++KE1ZswYq3fv3tb58+eD6uXo0aPWk08+aVVXV1s1NTXWjBkzrIMHDza7znvvvWdNnjzZ8vv9VmVlpZWWlmadPn06qJ4sy7JKSkqsn/70p9aSJUuCOr6+vt4aPHiw5ff7g+7Bsizr3Llz1pAhQ6yLFy9aNTU11pQpU6xDhw4FXe+f//ynlZ6ebl25cqXZx966dcsaMGCAdeXKFcvv91tZWVnW0aNHm13n6NGj1pgxY6zr169btbW11rx586wDBw7c0bHfvOYqKyutoUOHWufOnbP8fr+Vk5Nzx+fnm7UuXbpkzZs3z0pJSbHy8/Ob/XU5Rbjyb1mh3wPIf2Dhzr9lBX8PIP9mIP9NI//fRv4Ds8XKfUlJiR5++GF16tRJcXFxGjFihP72t781u05eXp6WL1+u+Pj4oHvp2rWrli5dqtjYWHm9Xv3oRz/SZ5991uw6AwcOVG5urmJiYnTlyhXV1dUpLi4uqJ6+/PJLbdy4UfPnzw/qeEn65JNP5PF4NGfOHD3xxBPavn17UHXefvttjR49Wt26dZPX69XGjRvVt2/foPt66aWXtHjxYnXu3LnZx9bV1am+vl6VlZWqra1VbW2t2rRp0+w6J0+e1JAhQ9ShQwe1bt1ajzzyiAoLC+/o2G9ec2VlZfrBD36gHj16KCYmRmPHjr3ja/mbtfbt26fHHntMo0aNavbX5CThyr8U+j2A/AcW7vxLwd8DyL8ZyH9g5P+7kf/AYoI6KswqKirUtWvXhs/j4+NVVlbW7DqrVq0KuZfExMSGj8+cOaOCggLt3LkzqFper1evvvqqtm7dqpEjRyohISGoOi+++KIWL16sixcvBnW8JF27dk0+n08vvfSSqqqqNGPGDP3whz/U4MGDm1Xn7Nmz8nq9mjVrli5fvqy0tDQ9++yzQfVUUlKiqqqqoC/eDh066JlnntGoUaPUtm1bDRw4UA899FCz6/Tu3VurV6/WvHnz1K5dOxUVFcmyrDs69pvX3Hddy//5z3+CqjV79mxJ0gcffHBHxztVuPIvhX4PIP+BhTP/Umj3APJvBvIfGPn/buQ/MFus3H/XifR4PFHo5CunTp1STk6OlixZop49ewZdZ9GiRSotLdXFixeVl5fX7ON37dqle+65Rz6fL+geJKl///5at26d4uLi1LlzZ2VlZenw4cPNrlNXV6fS0lL99re/VV5env7+979r9+7dQfW0c+dOzZw5M6hjJenjjz9Wfn6+3n33XRUXF6tVq1basmVLs+v4fD5lZmZq+vTpmj17tlJTU+X1eoPqyY7Xst3Z8ZyR/+8WzvxLod0DyL8Z7HjOyP93I/+B2elatsVwn5CQoM8//7zh84qKipAerQnVBx98oKeeekrPPfecJkyYEFSN06dP6x//+IckqV27dho+fLjKy8ubXaegoEBHjx7VuHHj9Oqrr6qoqEirV69udp0TJ06otLS04XPLshQT0/wf3Nx9993y+Xzq3Lmz2rZtq8ceeyyoVZaamhq9//77GjZsWLOP/Z/i4mL5fD516dJFsbGxyszM1PHjx5td58aNG0pPT9e+ffu0bds2tWvXTj169AiqJ7tdy05gt3NG/hsXrvxLod8DyL8Z7HbOyH/jyH9gdrqWbTHcDxo0SKWlpfrvf/+ryspKHTx4UI8++mhUerl48aIWLFig9evXKyMjI+g6Fy5c0AsvvKCamhrV1NTonXfeUWpqarPrvP7663rrrbe0d+9eLVq0SMOGDdOyZcuaXef69etat26dqqurdePGDe3evVvp6enNrpOWlqbi4mJdu3ZNdXV1OnLkiHr37t3sOuXl5erZs2fQzyFK0gMPPKCSkhLdunVLlmWpqKhIycnJza5z4cIFLViwQLW1tbp+/bp27doV9KNCffv21aeffqqzZ8+qrq5Ob731VtSuZacg/40zNf9S6PcA8m8G8t848t848h+YLZ65T0hI0OLFizVjxgz5/X5lZWUpJSUlKr1s2bJF1dXVWrNmTcNrkydP1pQpU5pVZ+jQofroo480fvx4tW7dWsOHDw/pZhGqtLS0hn7q6+s1depU9e/fv9l1+vbtq9mzZ2vq1Kny+/0aPHiwJk6c2Ow658+fV7du3Zp93O2GDBmikydPKjMzU16vV8nJyZo7d26z6zzwwAMaPny4nnjiCdXV1empp54K6kYsSW3atNGaNWu0cOFCVVdXa+jQoRo5cmRQtdyC/Lc8u+VfCv0eQP7NQP5bHvlvnKn591h3+psDAAAAAGzNFo/lAAAAAAgdwz0AAABgCIZ7AAAAwBAM9wAAAIAhGO4BAAAAQzDcAwAAAIZguAcAAAAMwXAPAAAAGOL/AevUYd3TbzMuAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 936x432 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, ax = plt.subplots(1,3,figsize = (13,6))\n",
    "for i,p in enumerate([r\"$\\nwarrow$\",r\"$\\uparrow$\",r\"$\\leftarrow$\"]):\n",
    "        sns.heatmap(tracer[:,:,i],linecolor='white',linewidth=1,\n",
    "                    cmap=\"coolwarm\",vmin=-0.5,vmax=0.6,square=True,cbar=False,ax=ax[i]);\n",
    "        ax[i].set_title(p)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Traceback\n",
    "When we pass `tracer` to the `traceback` function, it constructs 2x2 blocks iteratively from the 3 channels shown above and moves the block along the direction specified by the arrows."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "WARNING: 100 iterations exceeded; <10 alignments found: to continue searching, increase max_i.\n",
      "1. \n",
      " AAA-ATGC-CA \n",
      " AAAAATGCGCC \n",
      "\n",
      "2. \n",
      " -AAAATGC-CA \n",
      " AAAAATGCGCC \n",
      "\n",
      "3. \n",
      " AAAA-TGC-CA \n",
      " AAAAATGCGCC \n",
      "\n",
      "4. \n",
      " A-AAATGC-CA \n",
      " AAAAATGCGCC \n",
      "\n",
      "5. \n",
      " AA-AATGC-CA \n",
      " AAAAATGCGCC \n",
      "\n"
     ]
    }
   ],
   "source": [
    "alignment = seq.traceback_iterator(tracer,seq_1,seq_2)\n",
    "for i,e in enumerate(alignment):\n",
    "    print(str(i+1)+\".\",\"\\n\",e[0],\"\\n\",e[1],\"\\n\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are atleast 5 possible optimal alignments for the sequences `AAAAATGCGCC` and `AAAATGCCA`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Example Alignments\n",
    "\n",
    "Let's align some more sequences.\n",
    "\n",
    "### 1."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "seq_2 = \"AAATTGCGCAAAATCTACTACGCGCGCGA\"\n",
    "seq_1 = \"AAATTGCGCAAAATCTTACGCGCCCTACGCGA\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Max. score = 12.0\n",
      "1. \n",
      " AAATTGCGCAAAATCTTAC-GCGCCCTACGCGA \n",
      " AAATTGCGCAAAATCT-ACTACGC---GCGCGA \n",
      "\n",
      "2. \n",
      " AAATTGCGCAAAATCT--TACGCGCCCTACGCGA \n",
      " AAATTGCGCAAAATCTACTACGCG-C----GCGA \n",
      "\n",
      "3. \n",
      " AAATTGCGCAAAATCTTAC-GCGCCCTACGCGA \n",
      " AAATTGCGCAAAATC-TACTACG--C-GCGCGA \n",
      "\n",
      "4. \n",
      " AAATTGCGCAAAATCT--TACGCGCCCTACGCGA \n",
      " AAATTGCGCAAAATCTACTACGCGC-----GCGA \n",
      "\n",
      "5. \n",
      " AAATTGCGCAAAATCTTAC-GCGCCCTACGCGA \n",
      " AAATTGCGCAAAATCT-ACTACG-C--GCGCGA \n",
      "\n",
      "6. \n",
      " AAATTGCGCAAAATCT--TACGCGCCCTACGCGA \n",
      " AAATTGCGCAAAATCTACTACGCG-----CGCGA \n",
      "\n",
      "7. \n",
      " AAATTGCGCAAAATCTTAC-GCGCCCTACGCGA \n",
      " AAATTGCGCAAAATC-TACTACGC---GCGCGA \n",
      "\n",
      "8. \n",
      " AAATTGCGCAAAATCT--TACGCGCCCTACGCGA \n",
      " AAATTGCGCAAAATCTACTACGCG--C---GCGA \n",
      "\n",
      "9. \n",
      " AAATTGCGCAAAATCTTAC-GCGCCCTACGCGA \n",
      " AAATTGCGCAAAATC-TACTACG-C--GCGCGA \n",
      "\n",
      "10. \n",
      " AAATTGCGCAAAATCTTAC-GCGCCCTACGCGA \n",
      " AAATTGCGCAAAATCT-ACTACG--C-GCGCGA \n",
      "\n"
     ]
    }
   ],
   "source": [
    "_,_,_ = seq.global_dp(seq_1,seq_2,S,g,disp_alignments=True,max_i=200)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "seq_2 = \"ATCCCCGGGGGGGCTCTCTC\"\n",
    "seq_1 = \"AT\""
   ]
  },
  {
   "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 = -104.0\n",
      "1. \n",
      " A-------------T----- \n",
      " ATCCCCGGGGGGGCTCTCTC \n",
      "\n",
      "2. \n",
      " AT------------------ \n",
      " ATCCCCGGGGGGGCTCTCTC \n",
      "\n",
      "3. \n",
      " A-----------------T- \n",
      " ATCCCCGGGGGGGCTCTCTC \n",
      "\n",
      "4. \n",
      " A---------------T--- \n",
      " ATCCCCGGGGGGGCTCTCTC \n",
      "\n"
     ]
    }
   ],
   "source": [
    "_,_,_ = seq.global_dp(seq_1,seq_2,S,g,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
}