GREMLIN_TF.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "GREMLIN_TF_v2.ipynb",
      "provenance": [],
      "collapsed_sections": [],
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    },
    "accelerator": "GPU"
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/sokrypton/GREMLIN_CPP/blob/master/GREMLIN_TF.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "mu2S86VhS-8t"
      },
      "source": [
        "# GREMLIN_TF v2.1\n",
        "GREMLIN implemented in tensorflow\n",
        "\n",
        "### Change log:\n",
        "*   02Apr2019\n",
        " - fixing a few hard-coded values, to allow GREMLIN to work with any alphabet (binary, protein, rna etc)\n",
        "*   22Jan2019\n",
        " - moving [GREMLIN_TF_simple](https://colab.research.google.com/github/sokrypton/GREMLIN_CPP/blob/master/GREMLIN_TF_simple.ipynb) to a seperate notebook\n",
        "*   19Jan2019\n",
        " - in the past we found that optimizing V first, required less iterations for convergence. Since V can be computed exactly (assuming no W), we replace this first optimization step with a simple V initialization.\n",
        " - a few variables were renamed to be consistent with the c++ version\n",
        "*   16Jan2019\n",
        " - updating how indices are handled (for easier/cleaner parsing)\n",
        " - minor speed up in how we symmetrize and zero the diagional of W\n",
        "*   15Jan2019\n",
        " - LBFGS optimizer replaced with a modified version of the ADAM optimizer\n",
        " - Added option for stochastic gradient descent (via batch_size)\n",
        "  \n",
        "### Method:\n",
        "GREMLIN takes a multiple sequence alignment (MSA) and returns a Markov Random Field (MRF). The MRF consists of a one-body term (V) that encodes conservation, and a two-body term (W) that encodes co-evolution.\n",
        "\n",
        "For more details about the method see:\n",
        "[Google slides](https://docs.google.com/presentation/d/1aooxoksosSv7CWs9-ktqhUjyXR3wrgbG5a6PCr92od4/) and accompanying [Google colab](https://colab.research.google.com/drive/17RJcExuyifnd7ShTcsZGh6mBpWq0-s60)\n",
        "\n",
        "See [GREMLIN_TF_simple](https://colab.research.google.com/github/sokrypton/GREMLIN_CPP/blob/master/GREMLIN_TF_simple.ipynb) for a stripped down version of this code (with no funky gap removal, sequence weight, etc). This is intented for educational purpose,  and could also be very useful for anyone trying to modify or improve the algorithm!\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "yM3wyYU5SwYn"
      },
      "source": [
        "# ------------------------------------------------------------\n",
        "# \"THE BEERWARE LICENSE\" (Revision 42):\n",
        "# <so@g.harvard.edu> and <pkk382@g.harvard.edu> wrote this code.\n",
        "# As long as you retain this notice, you can do whatever you want\n",
        "# with this stuff. If we meet someday, and you think this stuff\n",
        "# is worth it, you can buy us a beer in return.\n",
        "# --Sergey Ovchinnikov and Peter Koo\n",
        "# ------------------------------------------------------------"
      ],
      "execution_count": 1,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "NLUvPVyxb7bo"
      },
      "source": [
        "## libraries"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "YyJpLM_tJfrY"
      },
      "source": [
        "# IMPORTANT, only tested using PYTHON 3!\n",
        "import numpy as np\n",
        "import tensorflow.compat.v1 as tf\n",
        "tf.disable_eager_execution()\n",
        "import matplotlib.pylab as plt\n",
        "from scipy import stats\n",
        "from scipy.spatial.distance import pdist,squareform\n",
        "import pandas as pd"
      ],
      "execution_count": 10,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "j0Yp7bPRmvwU"
      },
      "source": [
        "## Params"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "o3c7KURqmugY"
      },
      "source": [
        "################\n",
        "# note: if you are modifying the alphabet\n",
        "# make sure last character is \"-\" (gap)\n",
        "################\n",
        "alphabet = \"ARNDCQEGHILKMFPSTWYV-\"\n",
        "states = len(alphabet)\n",
        "a2n = {}\n",
        "for a,n in zip(alphabet,range(states)):\n",
        "  a2n[a] = n\n",
        "################\n",
        "\n",
        "def aa2num(aa):\n",
        "  '''convert aa into num'''\n",
        "  if aa in a2n: return a2n[aa]\n",
        "  else: return a2n['-']"
      ],
      "execution_count": 11,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "bX6GXKV3I2pm"
      },
      "source": [
        "## Functions for prepping the MSA (Multiple sequence alignment)"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "zA0Bne59SUIu"
      },
      "source": [
        "# from fasta\n",
        "def parse_fasta(filename,limit=-1):\n",
        "  '''function to parse fasta'''\n",
        "  header = []\n",
        "  sequence = []\n",
        "  lines = open(filename, \"r\")\n",
        "  for line in lines:\n",
        "    line = line.rstrip()\n",
        "    if line[0] == \">\":\n",
        "      if len(header) == limit:\n",
        "        break\n",
        "      header.append(line[1:])\n",
        "      sequence.append([])\n",
        "    else:\n",
        "      sequence[-1].append(line)\n",
        "  lines.close()\n",
        "  sequence = [''.join(seq) for seq in sequence]\n",
        "  return np.array(header), np.array(sequence)\n",
        "\n",
        "def filt_gaps(msa,gap_cutoff=0.5):\n",
        "  '''filters alignment to remove gappy positions'''\n",
        "  tmp = (msa == states-1).astype(np.float)\n",
        "  non_gaps = np.where(np.sum(tmp.T,-1).T/msa.shape[0] < gap_cutoff)[0]\n",
        "  return msa[:,non_gaps],non_gaps\n",
        "\n",
        "def get_eff(msa,eff_cutoff=0.8):\n",
        "  '''compute effective weight for each sequence'''\n",
        "  ncol = msa.shape[1]\n",
        "  \n",
        "  # pairwise identity\n",
        "  msa_sm = 1.0 - squareform(pdist(msa,\"hamming\"))\n",
        "\n",
        "  # weight for each sequence\n",
        "  msa_w = (msa_sm >= eff_cutoff).astype(np.float)\n",
        "  msa_w = 1/np.sum(msa_w,-1)\n",
        "  \n",
        "  return msa_w\n",
        "\n",
        "def mk_msa(seqs):\n",
        "  '''converts list of sequences to msa'''\n",
        "  \n",
        "  msa_ori = []\n",
        "  for seq in seqs:\n",
        "    msa_ori.append([aa2num(aa) for aa in seq])\n",
        "  msa_ori = np.array(msa_ori)\n",
        "  \n",
        "  # remove positions with more than > 50% gaps\n",
        "  msa, v_idx = filt_gaps(msa_ori,0.5)\n",
        "  \n",
        "  # compute effective weight for each sequence\n",
        "  msa_weights = get_eff(msa,0.8)\n",
        "\n",
        "  # compute effective number of sequences\n",
        "  ncol = msa.shape[1] # length of sequence\n",
        "  w_idx = v_idx[np.stack(np.triu_indices(ncol,1),-1)]\n",
        "  \n",
        "  return {\"msa_ori\":msa_ori,\n",
        "          \"msa\":msa,\n",
        "          \"weights\":msa_weights,\n",
        "          \"neff\":np.sum(msa_weights),\n",
        "          \"v_idx\":v_idx,\n",
        "          \"w_idx\":w_idx,\n",
        "          \"nrow\":msa.shape[0],\n",
        "          \"ncol\":ncol,\n",
        "          \"ncol_ori\":msa_ori.shape[1]}\n"
      ],
      "execution_count": 12,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "fky6gk-HlFyi"
      },
      "source": [
        "## GREMLIN"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "gx14M7Tvu-Ct"
      },
      "source": [
        "# external functions\n",
        "\n",
        "def sym_w(w):\n",
        "  '''symmetrize input matrix of shape (x,y,x,y)'''\n",
        "  x = w.shape[0]\n",
        "  w = w * np.reshape(1-np.eye(x),(x,1,x,1))\n",
        "  w = w + tf.transpose(w,[2,3,0,1])\n",
        "  return w\n",
        "\n",
        "def opt_adam(loss, name, var_list=None, lr=1.0, b1=0.9, b2=0.999, b_fix=False):\n",
        "  # adam optimizer\n",
        "  # Note: this is a modified version of adam optimizer. More specifically, we replace \"vt\"\n",
        "  # with sum(g*g) instead of (g*g). Furthmore, we find that disabling the bias correction\n",
        "  # (b_fix=False) speeds up convergence for our case.\n",
        "  \n",
        "  if var_list is None: var_list = tf.trainable_variables() \n",
        "  gradients = tf.gradients(loss,var_list)\n",
        "  if b_fix: t = tf.Variable(0.0,\"t\")\n",
        "  opt = []\n",
        "  for n,(x,g) in enumerate(zip(var_list,gradients)):\n",
        "    if g is not None:\n",
        "      ini = dict(initializer=tf.zeros_initializer,trainable=False)\n",
        "      mt = tf.get_variable(name+\"_mt_\"+str(n),shape=list(x.shape), **ini)\n",
        "      vt = tf.get_variable(name+\"_vt_\"+str(n),shape=[], **ini)\n",
        "      \n",
        "      mt_tmp = b1*mt+(1-b1)*g\n",
        "      vt_tmp = b2*vt+(1-b2)*tf.reduce_sum(tf.square(g))\n",
        "      lr_tmp = lr/(tf.sqrt(vt_tmp) + 1e-8)\n",
        "\n",
        "      if b_fix: lr_tmp = lr_tmp * tf.sqrt(1-tf.pow(b2,t))/(1-tf.pow(b1,t))\n",
        "\n",
        "      opt.append(x.assign_add(-lr_tmp * mt_tmp))\n",
        "      opt.append(vt.assign(vt_tmp))\n",
        "      opt.append(mt.assign(mt_tmp))\n",
        "        \n",
        "  if b_fix: opt.append(t.assign_add(1.0))\n",
        "  return(tf.group(opt))"
      ],
      "execution_count": 13,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "BqYqlJAXVI9N"
      },
      "source": [
        "def GREMLIN(msa, opt_type=\"adam\", opt_iter=100, opt_rate=1.0, batch_size=None):\n",
        "  \n",
        "  ##############################################################\n",
        "  # SETUP COMPUTE GRAPH\n",
        "  ##############################################################\n",
        "  # kill any existing tensorflow graph\n",
        "  tf.reset_default_graph()\n",
        "  \n",
        "  ncol = msa[\"ncol\"] # length of sequence\n",
        "\n",
        "  # msa (multiple sequence alignment) \n",
        "  MSA = tf.placeholder(tf.int32,shape=(None,ncol),name=\"msa\")\n",
        "  \n",
        "  # one-hot encode msa\n",
        "  OH_MSA = tf.one_hot(MSA,states)\n",
        "  \n",
        "  # msa weights\n",
        "  MSA_weights = tf.placeholder(tf.float32, shape=(None,), name=\"msa_weights\")\n",
        "  \n",
        "  # 1-body-term of the MRF\n",
        "  V = tf.get_variable(name=\"V\", \n",
        "                      shape=[ncol,states],\n",
        "                      initializer=tf.zeros_initializer)\n",
        "\n",
        "  # 2-body-term of the MRF\n",
        "  W = tf.get_variable(name=\"W\",\n",
        "                      shape=[ncol,states,ncol,states],\n",
        "                      initializer=tf.zeros_initializer)\n",
        "  \n",
        "  # symmetrize W\n",
        "  W = sym_w(W)\n",
        "  \n",
        "  def L2(x): return tf.reduce_sum(tf.square(x))\n",
        "  \n",
        "  ########################################\n",
        "  # V + W\n",
        "  ########################################\n",
        "  VW = V + tf.tensordot(OH_MSA,W,2)\n",
        "\n",
        "  # hamiltonian\n",
        "  H = tf.reduce_sum(tf.multiply(OH_MSA,VW),axis=(1,2))\n",
        "  # local Z (parition function)\n",
        "  Z = tf.reduce_sum(tf.reduce_logsumexp(VW,axis=2),axis=1)\n",
        "\n",
        "  # Psuedo-Log-Likelihood\n",
        "  PLL = H - Z\n",
        "\n",
        "  # Regularization\n",
        "  L2_V = 0.01 * L2(V)\n",
        "  L2_W = 0.01 * L2(W) * 0.5 * (ncol-1) * (states-1)\n",
        "\n",
        "  # loss function to minimize\n",
        "  loss = -tf.reduce_sum(PLL*MSA_weights)/tf.reduce_sum(MSA_weights)\n",
        "  loss = loss + (L2_V + L2_W)/msa[\"neff\"]\n",
        "\n",
        "  ##############################################################\n",
        "  # MINIMIZE LOSS FUNCTION\n",
        "  ##############################################################\n",
        "  if opt_type == \"adam\":  \n",
        "    opt = opt_adam(loss,\"adam\",lr=opt_rate)\n",
        "\n",
        "  # generate input/feed\n",
        "  def feed(feed_all=False):\n",
        "    if batch_size is None or feed_all:\n",
        "      return {MSA:msa[\"msa\"], MSA_weights:msa[\"weights\"]}\n",
        "    else:\n",
        "      idx = np.random.randint(0,msa[\"nrow\"],size=batch_size)\n",
        "      return {MSA:msa[\"msa\"][idx], MSA_weights:msa[\"weights\"][idx]}\n",
        "  \n",
        "  # optimize!\n",
        "  with tf.Session() as sess:\n",
        "    # initialize variables V and W\n",
        "    sess.run(tf.global_variables_initializer())\n",
        "    \n",
        "    # initialize V\n",
        "    msa_cat = tf.keras.utils.to_categorical(msa[\"msa\"],states)\n",
        "    pseudo_count = 0.01 * np.log(msa[\"neff\"])\n",
        "    V_ini = np.log(np.sum(msa_cat.T * msa[\"weights\"],-1).T + pseudo_count)\n",
        "    V_ini = V_ini - np.mean(V_ini,-1,keepdims=True)\n",
        "    sess.run(V.assign(V_ini))\n",
        "    \n",
        "    # compute loss across all data\n",
        "    get_loss = lambda: round(sess.run(loss,feed(feed_all=True)) * msa[\"neff\"],2)\n",
        "    print(\"starting\",get_loss())\n",
        "    \n",
        "    if opt_type == \"lbfgs\":\n",
        "      lbfgs = tf.contrib.opt.ScipyOptimizerInterface\n",
        "      opt = lbfgs(loss,method=\"L-BFGS-B\",options={'maxiter': opt_iter})\n",
        "      opt.minimize(sess,feed(feed_all=True))\n",
        "      \n",
        "    if opt_type == \"adam\":\n",
        "      for i in range(opt_iter):\n",
        "        sess.run(opt,feed())  \n",
        "        if (i+1) % int(opt_iter/10) == 0:\n",
        "          print(\"iter\",(i+1),get_loss())\n",
        "    \n",
        "    # save the V and W parameters of the MRF\n",
        "    V_ = sess.run(V)\n",
        "    W_ = sess.run(W)\n",
        "    \n",
        "  # only return upper-right triangle of matrix (since it's symmetric)\n",
        "  tri = np.triu_indices(ncol,1)\n",
        "  W_ = W_[tri[0],:,tri[1],:]\n",
        "  \n",
        "  mrf = {\"v\": V_,\n",
        "         \"w\": W_,\n",
        "         \"v_idx\": msa[\"v_idx\"],\n",
        "         \"w_idx\": msa[\"w_idx\"]}\n",
        "  \n",
        "  return mrf"
      ],
      "execution_count": 14,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "mppg0JLtP25z"
      },
      "source": [
        "## EXAMPLE"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "bnigmLmAlyWv"
      },
      "source": [
        "# download example fasta MSA\n",
        "!wget -q -nc https://gremlin2.bakerlab.org/db/PDB_EXP/fasta/4FAZA.fas"
      ],
      "execution_count": 15,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "osaZwTSMOicF"
      },
      "source": [
        "# ===============================================================================\n",
        "# PREP MSA\n",
        "# ===============================================================================\n",
        "# parse fasta\n",
        "names, seqs = parse_fasta(\"4FAZA.fas\")\n",
        "\n",
        "# process input sequences\n",
        "msa = mk_msa(seqs)"
      ],
      "execution_count": 16,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "CoBRuqmVVrbD",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "4c0511ab-2010-499f-cd45-394a48c9dba4"
      },
      "source": [
        "%%time\n",
        "# ===============================================================================\n",
        "# RUN GREMLIN\n",
        "# ===============================================================================\n",
        "# Note: the original GREMLIN uses the \"lbfgs\" optimizer which is EXTREMELY slow \n",
        "# in tensorflow. The modified adam optimizer is much faster, but may \n",
        "# require adjusting number of iterations (opt_iter) to converge to the same \n",
        "# solution. To switch back to the original, set opt_type=\"lbfgs\".\n",
        "# ===============================================================================\n",
        "mrf = GREMLIN(msa)"
      ],
      "execution_count": 17,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "starting 50463.14\n",
            "iter 10 29873.73\n",
            "iter 20 27907.74\n",
            "iter 30 26972.23\n",
            "iter 40 26687.46\n",
            "iter 50 26576.22\n",
            "iter 60 26540.93\n",
            "iter 70 26524.79\n",
            "iter 80 26518.27\n",
            "iter 90 26515.09\n",
            "iter 100 26513.35\n",
            "CPU times: user 1.89 s, sys: 590 ms, total: 2.47 s\n",
            "Wall time: 9.53 s\n"
          ]
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "jCrfC2Um4xww"
      },
      "source": [
        "## Explore the contact map\n",
        "### Contact prediction:\n",
        "\n",
        "For contact prediction, the W matrix is reduced from LxLx21x21 to LxL matrix (by taking the L2norm for each of the 20x20). In the code below, you can access this as mtx[\"raw\"]. Further correction (average product correction) is then performed to the mtx[\"raw\"] to remove the effects of entropy, mtx[\"apc\"]. The relative ranking of mtx[\"apc\"] is used to assess importance. When there are enough effective sequences (>1000), we find that the top 1.0L contacts are ~90% accurate! When the number of effective sequences is lower, NN can help clean noise and fill in missing contacts.\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "XxgQArVUPyPH"
      },
      "source": [
        "## Functions for extracting contacts from MRF"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "nMxp7up_P1_q"
      },
      "source": [
        "###################\n",
        "def normalize(x):\n",
        "  x = stats.boxcox(x - np.amin(x) + 1.0)[0]\n",
        "  x_mean = np.mean(x)\n",
        "  x_std = np.std(x)\n",
        "  return((x-x_mean)/x_std)\n",
        "\n",
        "def get_mtx(mrf):\n",
        "  '''get mtx given mrf'''\n",
        "  \n",
        "  # l2norm of 20x20 matrices (note: we ignore gaps)\n",
        "  raw = np.sqrt(np.sum(np.square(mrf[\"w\"][:,:-1,:-1]),(1,2)))\n",
        "  raw_sq = squareform(raw)\n",
        "\n",
        "  # apc (average product correction)\n",
        "  ap_sq = np.sum(raw_sq,0,keepdims=True)*np.sum(raw_sq,1,keepdims=True)/np.sum(raw_sq)\n",
        "  apc = squareform(raw_sq - ap_sq, checks=False)\n",
        "\n",
        "  mtx = {\"i\": mrf[\"w_idx\"][:,0],\n",
        "         \"j\": mrf[\"w_idx\"][:,1],\n",
        "         \"raw\": raw,\n",
        "         \"apc\": apc,\n",
        "         \"zscore\": normalize(apc)}\n",
        "  return mtx\n",
        "\n",
        "def plot_mtx(mtx,key=\"zscore\",vmin=1,vmax=3):\n",
        "  '''plot the mtx'''\n",
        "  plt.figure(figsize=(5,5))\n",
        "  plt.imshow(squareform(mtx[key]), cmap='Blues', interpolation='none', vmin=vmin, vmax=vmax)\n",
        "  plt.grid(False)\n",
        "  plt.show()"
      ],
      "execution_count": 18,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "RSleviAVPJ36",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 321
        },
        "outputId": "8c34f3ea-ff91-4e36-ee01-46ff05e15a3f"
      },
      "source": [
        "mtx = get_mtx(mrf)  \n",
        "plot_mtx(mtx)"
      ],
      "execution_count": 19,
      "outputs": [
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAATUAAAEwCAYAAADW0ua3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3de3xV1Zk38N9DCBDugYQYLgIKongBJHVkpGpldBCc4rS1xXdqdeStbV+r0DoqVdvPjKOjtW+rjvVSRn2lthW1aqHoUC3V8YO3GhTkVgQE5BJIouF+D8/7xznUrLU27H32ue2s/L6fDx/yLPbeZwEnT/Z51tpriaqCiMgX7YrdASKiXGJSIyKvMKkRkVeY1IjIK0xqROQVJjUi8kpWSU1ExovIShFZLSLTc9UpIqK4JO48NREpAfAhgAsBbATwLoDLVXX50c6pqKjQgQMHHfO6hw6b/WnfTmL1jz5j/w8H/Ytu33fQiHt0Ks1bf4rpYLP5r1FaEv7+ivLvR/n33nsLG1W1Muy49lm8xlkAVqvqRwAgIrMATAJw1KQ2cOAgvPFO7TEv2rT7gBGXd+mQRRcJAA41Hzbi9iXuDfpLy+qMeMKp1XntU7Fs2bbPiI/r2Sn0HPsHvwjTWjGUlcr6KMdl8/GzH4ANLeKN6TaDiFwjIrUiUtvQ2JDFyxERhcv7QIGqzlDVGlWtqawIvXMkIspKNh8/NwEY0CLun27Liv1x8721Tc4xZw4uN+IVm3Y4x5zSr7sRHz7s1g7bRajX7d53yIi7dMrmn6w4gj5u2uJ83Gzcud+IK7p1zPgaUR04ZH6E7tA+3s/jPt2P3cdXV9Y7bV8Y1seIR9w6zzlm8Z3jY/UnH77+y4VO26++MdqI1zfucY4ZWNE5J6+/eP02Ix4xsGdOrhtVNndq7wIYKiKDRaQDgMkA5uSmW0RE8cS+7VDVQyLyXQB/AFAC4HFVXZaznhERxZDVZylVfQnASznqCxFR1vhEARF5Jfbk2zhGjBqt8157669xlDloGz/d67T9cc1WI77qc4Oy7ltrsGrLLqdt6HFdM77O5U+4cwWfuqrmmOds33PQaRt00a1G3LTgnoz7AsSbO0ZtT1mpLFTVY79RwTs1IvIMkxoReYVJjYi8UtCa2ujRNdry2c+gibV9epj1lP69ykKvW/7lR5y2pue+bcSvf+g+onXuSf494bDNenYWAHry+dnI1jXsdtoGVXbJ+Dpn3zHfaXv7tnGh59mTYnM1ITaKpNc2WVMjojaJSY2IvMKkRkReYVIjIq8UdaAgyBPvrjPitjKxNoofzVvptN0+fljG1ym/5F6nrWnu9zK+zt4DzUZc1qEk42u0BrOXmIvPTDrdWTbQUb99n9NmD4JRZjhQQERtEpMaEXmFSY2IvFLUmlqUFWujTKwNsnTDdiM+bUCPqN0kj4VtorLLWukYALpaqx2/v26bc8yoQebqrn2u+KVzTP2T3zDiN1Y3OsecM6TCaaMU1tSIqE1iUiMirzCpEZFXmNSIyCtF3e/NHhQIEmdQAHAHBg5aW6wBQKm1zZq9tRcA3P/mOiN+/PKRof2Jon6Hub1c2NZtSWQXzO1ieVQPLFhjxNeNPTF2n8KE7a5uDwoECfp7TvudueeQPSgAAAtWmQMDY4fmZ1Bg9353sKNLx9a3tWNcvFMjIq8wqRGRV5jUiMgrRZ18e/iw+9rt2pk1j7gr1to1NLt+BgDfnLXYiP9r8gjnmJWbdxrxsL7dQl+bkutQs/m+aF+S+c/1KPXZXHltpfv+P39Y9is2b25yd2nrW26uMv2xtQovABwfYSXeC+9fYMSvTB2bYe+CcfItEbVJTGpE5BUmNSLyCpMaEXmlqDPy7EGBIHG3sbMLt/agAOAODPzOWuEUAHp1NLeXGwYOFBSDvZJsMVeRzdegAACUT/ypEb/20JS8vI49KBAkyqBAkFwNDMTFOzUi8gqTGhF5hUmNiLzi7VOu9sPpQRNr7RrapQG7BF391CIjjlvjo+zkqoYWZ7JtITW9eEOxu9DqJft/mIgoQ0xqROQVJjUi8kpoUhORx0WkXkSWtmjrJSKviMiq9O/l+e0mEVE0UQYKngDwcwAt9/yaDmC+qt4tItPT8c2ZvvjugO3IukRYeTQKe8Xam0vd1VTtibX2oACQu5VuKTv2ii5RJm7bK3IAyR8ooOyF/g+r6usAPrWaJwGYmf56JoBLc9wvIqJY4v7YqlLVuvTXWwBUHe1AEblGRGpFpLah0V0biogol7K+F9fUKpNHXWlSVWeoao2q1lRWcI4XEeVX3ALWVhGpVtU6EakGUB/nIrmqnwWJUguzH07P1cTamn99xYjXvPmuc8wjP/66EX9t1PHOMY07zR2nKrq1vh2nciVKDc3G+lnbFPd/fQ6AK9NfXwlgdm66Q0SUnShTOp4C8BaAYSKyUUSmALgbwIUisgrA36VjIqKiC/38p6qXH+WPxuW4L0REWfP2gfZimn/T+Ubco/OFoedMmeXOkXtscn7myO0J2MG7c8gO3tt2H3DaenbpEHAk5dvZd8w34jnXm4sy9umev9prszVf8Mbfr3CO+dmk4Ua81VrgEwB2WXNUT6zqmoPepbCSSkReYVIjIq8wqRGRV5jUiMgrHCjIkD2x1h4UAIAenUuNeO7Szc4xl5zW14g31O/KvnMRhQ0KAEDqQZHPxB0UsAvLJTEm0ZLp7duKN/HA/v+zBwWCVAWsWlzVI2ddcvBOjYi8wqRGRF5hUiMir7CmliH74fSgibV2Dc2unwV5+Xp3V+v6HeYD7fmcVGkTMWsndl+AaP2JU0N7dtEGp+2ykQMyvg61TbxTIyKvMKkRkVeY1IjIK0xqROQVbwcK8lVkt1esDRJlYCCKQg4MhClkX1rjoED5Rf9hxE0v3+Ic02StdFKeo1VOJj78lhG/+J0xObnuuobdTlv/XmVGnMTVhZPXIyKiLDCpEZFXmNSIyCtiP7icT6NH1+gb79QW7PUKIWjFWvvh9KCJtVHYNY1BlV2cYz7ZZdZpSgLmutoPo6/cvNM5Zljfbk4bJVOU94Vtc9NeI+5bXnaUIzPz+4DFGv7Bqil/8PF255jzJt9uxE1v/jT0tcpKZaGq1oQdxzs1IvIKkxoReYVJjYi8wqRGRF7hQEGGGneak3oruoVPSo2ywkXQRMcoBeCdew8a8co6dwXdm2cvMeKXp37eOcZeTeNQ82Ejfn1Vo3POBSf3Ce3f8o07jPiyBxY4xyz78QQjXrPV/Tvkcgs1ii9XKxl/3LjHiI+v6Bx6DgcKiKhNYlIjIq8wqRGRVxJXU1u1xaynDD3OraX8aN5KI759/LDsO9eGXPb4u07bs1d/7pjn2JM3AeC07zxlxJ/Oujr0tT+scyf+Vvc0dxvqVlbqHJN0761tMuIzB5cXqSfx2XXUpD2szpoaEbVJTGpE5BUmNSLyCpMaEXklcQMFbZW92gYAdLCW3GiNBfS24pQbXzTiFT+ZWKSeBCs/71YjbvqfO0PPOXjosNNW2t68Dxpz55+cY9669YLQa9srxURZJYYDBUTUJjGpEZFXQpOaiAwQkVdFZLmILBORqen2XiLyioisSv/e+ibmEJF3QmtqIlINoFpV3xORbgAWArgUwFUAPlXVu0VkOoByVb35WNeKU1PbttutNdkrufog6O+5eqv5kHvNCfF+bkx+wvw3n3VVaFkCw296yYhfuvE855goD9zvO9BsxJ06lISeYy8aAERbOCDOirB0dCs27XDarOfZcWr/7gXqTQ5raqpap6rvpb/eCWAFgH4AJgGYmT5sJlKJjoioqDKqqYnIIACjALwDoEpV69J/tAVA1VHOuUZEakWktqGxIYuuEhGFi5zURKQrgOcATFNV475UU59hAz/HquoMVa1R1ZrKisqsOktEFCZSUhORUqQS2q9V9fl089Z0ve1I3a0+P10kIoouykCBIFUz+1RVp7Vo/wmAT1oMFPRS1ZuOda0zR9fo62/++a9x0lYBSJpx975uxPO/d27oOfagAOAODCxat805ZuSgnkY8e8kmI375w0+dcx788umh/Sk/e5oRdxvxt84xH//iq6HXIYo6UNA+wrXOAXAFgCUicmSTy1sA3A3gGRGZAmA9AL4ziajoQpOaqi4AcLSFyMfltjtERNnh5z8i8kriHmi/3KoJPRUwUbT8knuNuGnu97LvXEx79h9y2jp3jPKp3mQ/4AsAQ6xVf4PqZWEr1gJuDc2un0Vhr0gMAF976E0jnnGV2xd7wvDaenfXrMF9OEmWwvGBdiJqk5jUiMgrTGpE5BUmNSLySkEHCoacOkLv+c28v8YTTq0u2GtT/j36zlqn7X//zeCCvX6UQSbb1u37jLiqR6ejHOm/ecvrjPjGJxc5xyy56+JCdcfBgQIiapOY1IjIK0xqROSVzGeJZqFHp9I2UUez65SpNQEyd6jZ3M0nygIA9oq1AHDXFSOMeNLp/TLuyx1//NBpu+3vTjLi5xfWOcc8ONc87/1///uMXzuq3t0zr4e15Rqabfxw83tz/F2t83uVd2pE5BUmNSLyCpMaEXklcQ+0b99z0Ih7dE72ruRRdruq3+HujtSne/juSH/6i7mY8MlV7i7WfcvLjNjeUQkAfvL6R0Y87Rx37thQ6+F5u4Zm18+i+u7zS43414/Oc45peulfYl2b2hbOUyOiNolJjYi8wqRGRF5hUiMirxR1oKBxp1tAH3rxD424acE9zjF7DzQbcVmHEueY963VXkfFWO01aXpNftxp+3TW1Rlf58wfvey0vXf7RbH6lCn7AXIAONhsvgf79ypzjimUoP5xgm4ycKCAiNokJjUi8gqTGhF5JXGTbyn/aj9qctoWNZg1SPvh9Jf+j7uzer68veYTp23YcebE428/s9g5pqyDuT7Do5NHOMfUrjP/7hdPf96Im579ZuR+UmGxpkZEbRKTGhF5hUmNiLzCpEZEXuFAQR40Hzb/TUvaxVv5dvnGHUY8vH/30HP2WROTAaD63BuMuOnt+0KvM+qHfzDic0a5q+X+/EunhV7HFjRIUXNCecbXSbrd+w45bV06FXShae9woICI2iQmNSLyCpMaEXmFH/LzIG4NzXbZAwuM+IVpn3eOOananJTaKeDh/m4jzImza+vd1XEH9+lixPauT+UT/q9zzg/HDTHioAe/7Rqaj/WzIKyfZa8+YHGBKHinRkReYVIjIq8wqRGRV0KTmoh0EpE/i8hiEVkmIv+Wbh8sIu+IyGoReVpEOoRdi4go30In34qIAOiiqrtEpBTAAgBTAXwfwPOqOktEHgGwWFUfPta1oky+3bLNLA4e1zPeqqMPLFhjxNeNPTHWdZJk596DTlu3suJtIbjx071GnKsVa6e+sMxpO6HCfB9M/Xy8/88P63YasT3QcvGDbzrn/Pe1uVmh5MChw0bcoT0/KGUiZ5NvNWVXOixN/1IAFwD4bbp9JoBLY/aViChnIv2oEJESEVkEoB7AKwDWANimqkeeBdkIwH2OJnXuNSJSKyK1DY0NuegzEdFRRUpqqtqsqiMB9AdwFoCTo76Aqs5Q1RpVramsqIzZTSKiaDKaIaiq20TkVQBjAPQUkfbpu7X+ADblokNxa2i2JNXQnl20wWm7bOSA0PPWbN1lxCdWdQ09J2iHropuHUPPi8OuoQWtWHv2ib1Dr2PX0O7/x1OdY+xr2w/7A9Ee+LdraLZc1c+C5KqGZtdWi1lXTaIoo5+VItIz/XUZgAsBrADwKoCvpA+7EsDsfHWSiCiqKHdq1QBmikgJUknwGVWdKyLLAcwSkTsAvA/gsTz2k4goktCkpqofABgV0P4RUvU1IqLE4EQZIvJKUZcSsCcjAoWdkGivAtAnYJWJw9Yqtu1irMARZVAgSJSBAVu+BgWisLexi8qeWBtlwGHFJnegIIpx975uxDdcNNSIJ5xaHeu6hcSBgWPjnRoReYVJjYi8wqRGRF4pak2t2A/0BtXQbHFqaG3Vt59Z7LQ9/c+fCz3Pfjg9aGKtXUM7pZ870Xbgt5814mHD3fpYVS9zhd+k19Dsmi6Q/PfkJ7sOGHHvroVdwId3akTkFSY1IvIKkxoReYVJjYi8wn28PLOuwd3+blBll4AjTZc/Ya5I3Lu7OYjy8y+dFnqNsg65eTtFWW3DHhQAgPWPXGbE1z63xDnmtnHmZNvpL64w4rsnnuKcs+9gsxF3KnW3IcyXpA0KbN9jrhDSo7M7EbjQAwM23qkRkVeY1IjIK0xqROSV0N2kcinKblJJc6jZfOi+fQl/DhyN/W8F5O/f66L/XOC0De3Xw4gf/PLpzjG1HzUZcc0J5bntmEf27D/ktO3Ya7blaqXqKHK2mxQRUWvCpEZEXmFSIyKvMKkRkVcKOvn2YLNiy7bPVpvt091dpTXOZMOgwQ4R8zpBRWxbUFGbAwPR1a5rctqibJH3Yd1OIw7axs5esdZebQNwJ9bagwJA+MDAk7XrnbYragYe8xxfde7opoegtqThdywReYVJjYi8wqRGRF7h5FsCAGy1dtaqirAqsK38sv9y2pqe/WbsPrX00rI6Iw5asTbKw+mUsr5xj9M2sKJzEXoSHSffElGbxKRGRF5hUiMirxS1pvbqynrnmC8M62PEURY93LXPffC2a6fkz6dp7XJRhwOAix9804j/+9q/jd2nMPY8tDhz0Ox5dUDw3LpieWDBGqfturHmjl33ve4e84/W7ltRamy7A773GnbuN+Ioi5RGwZoaEbVJTGpE5BUmNSLyCpMaEXmloAMFZ46u0TfefvezF5d4O+XMXrLJiCed3i+rfpEfirnrUzHZf+8d1o5PAHDvG+uM+K4JJ+ezS3nBgQIiapOY1IjIK5GTmoiUiMj7IjI3HQ8WkXdEZLWIPC0ixd3BlIgImS0SORXACgBHts/+MYB7VXWWiDwCYAqAh491AYFZRxtx6zznmMV3jg/tiF1De3/dNueYUYN6GvHBQ+4ikaXteaPqk1zU0HI1sXZz016nrW95mRH/r5kLnWN+c+XojF/L/ns/tGitc0xrrKHFFem7WkT6A5gI4NF0LAAuAPDb9CEzAVyajw4SEWUi6q3KfQBuAnDkdqc3gG2qeuQZiY0AAocgReQaEakVkdqGxoasOktEFCY0qYnIJQDqVdW9V45AVWeoao2q1lRWVMa5BBFRZFFqaucA+KKITADQCama2v0AeopI+/TdWn8Am45xDSKigsho8q2InA/gX1T1EhF5FsBzLQYKPlDVh451fiFXvp32u2VGfN+lpxbkdaMqv+g/jLjp5VtCz3lvrbs70j89ZK5wseInE7PrGOXNrxaaK4R8fbS7QsgBa0CrQ4TBrOc/2GjEXzqjv3PMH5ZvMeK/H35c6HWTphCTb28G8H0RWY1Uje2xLK5FRJQTGS06pqqvAXgt/fVHAM7KfZeIiOLjRC0i8krid5M6+475Ttuc684x4tO++4xzTP2T38isc5QX9sqoXTxYkTjKxNqkq9+x32nr071jEXoSHR9oJ6I2iUmNiLzCpEZEXmFSIyKvJH6gIK4FqxqNeOzQioK8blRNuw8YcXmXtrtyU5wJp0kTZWJtFFEGVvYdMFe6/c2iDUZ89VmDYr12HIUcNOFAARG1SUxqROQVJjUi8krrnwkJ4I3VjU5b0mpotrZcQ7O1thparlasDRJlcnKnDuZKt3YNrfyL/+mc0zTn+qz6dTRxd4TLp9b1biIiCsGkRkReYVIjIq8wqRGRVxI/ULC+cY/TNrCisxGfMyTeoED5xJ8acdOLN8S6DrUtQYMCUSYQx1mxxF6xFgCeX1pvxJXdzEGnfA0KBKnu2algrxUV79SIyCtMakTkFSY1IvJK4mtqdv0srtdWuhspv/bQlJxcmwrv8GF3IYZ27Yo3ETTKBOI4NbSgXZ+G9O5qxCdWdXWOCbNy806nbVjfbhlfJ8j2PQeNuEfn0ljX+cVba2Odxzs1IvIKkxoReYVJjYi8wqRGRF5J/EBBrpw/rDIv1w3awu/t28ZlfJ11DbudtkGVXWL1qbXbufeg09atzCw2F3NQIC57xVp7tQ3AnVhrDwoA7sDA0GmzjfjtOy52zund1ZygGzQocNCaQFwaYfDDHhQA3IGBm+aucI6555JTQq/9rTGDjXha6BkpvFMjIq8wqRGRV5jUiMgrbaamli9zrh8b67yJD79lxC9+Z0wuuuMFu35WbPsOWrWw0oBa2AcbjfhLZ/R3jomy65P9cHrQxFq7hrbqvknOMXGsrDMn5J42oEfoOVEm1kapn+US79SIyCtMakTkFSY1IvJKUXdo//ov3V15fvWN3OzKs3u/tSBfR5YPKdwDC9Y4bV87vZ8R/2rRRueY7583JOPXirvr0ye7DhixPQfNV9yhnYjaJCY1IvIKkxoReSVSoUlE1gHYCaAZwCFVrRGRXgCeBjAIwDoAX1XVpvx0k4gomkgDBemkVqOqjS3a7gHwqareLSLTAZSr6s3Huo49UECf2dy012nrW14Wel75ebcacdP/3JmzPlHKD176ixHfNeFk55goK9YmydIN2502e7Ltko/dY04/3jxmbb27EMPgPvlZiKEQAwWTAMxMfz0TwKVZXIuIKCeiJjUF8LKILBSRa9JtVapal/56C4CqoBNF5BoRqRWR2oZGd58AIqJcijp5a6yqbhKRPgBeERHjflxVVUQCP8eq6gwAM4DUx8+sektEFCJSUlPVTenf60XkBQBnAdgqItWqWici1QDqj3mRAFF2X8+VuDWrQonbF9bQsmO/B19YXuccE1RDs9k1tPod+51j+nTvmGHvcrfrk11Dy9XD6kH1s3H3vm7E8793buh1gmr7IvEWAg39+CkiXUSk25GvAVwEYCmAOQCuTB92JYDZwVcgIiqcKHdqVQBeSGfN9gB+o6rzRORdAM+IyBQA6wF8NX/dJCKKJjSpqepHAEYEtH8CIPPF+ImI8ohPFBCRV4q6SkcUW7btc9qO69kpV11qVQ41H3ba2pfk5ufSPKtAPn54dU6u29oUcvAqLnvXpzgr1uaTPSgXZRCsfrv7fd6nh/l9zlU6iKhNYlIjIq8wqRGRVxK/HGyu6mcfB9RKji9ircTekT3KbuxB9TO7vrJ66y7nmFP6dQ+99o1PLjLi8Xe1/praHmv1YwDoHLICctLqZ0HsndOj1NDsh9ODJtZG+X74wLrOGce7rx1nIrldP8sG79SIyCtMakTkFSY1IvIKkxoReSXxAwW5kq9BgebD7uTlknbhqwv07xVeTLWvHXRdu2gc0J1Iltx1cbwTE2T7noNGvPdAs3NM2ECBr+wVa6OwBwWA4IGBOH7+xkdGvGOfO6hzy7iTYl2bd2pE5BUmNSLyCpMaEXmlbRYYcujG369w2n42aXjoefZE2t8v3ewc8w+n9Q29zpg7/2TEb916Qeg5vrInlEZZudW2O6C206VTbr5N7Ae9g1Z2rY4w2dyuHUb5e9q7PkXZ8SlX9bPyiT912k46e6QRP3LF6Jy8FsA7NSLyDJMaEXmFSY2IvMKkRkReKerKt4vXb3OOGTGwZ15e+8L7Fzhtr0wdm5fXovz7ZNcBp6131w5ZX9dePQWItoJKvtiDAkC8ARCbvY0dADx51eeMOFdbSK4JWDnmxKquRnz9C0udYwZYgyY/umgoV74loraHSY2IvMKkRkReSfxuUkm3NWAXnKoYq3jm8+Fh8stNc80J3/dcckqRehI8sbb28W8ZsV0/A4C3Vn9ixJ1LS5xj7MUHzhjQjTU1Imp7mNSIyCtMakTkFSY1IvIKV+nI0q6AVR2qYtT3z5t8u9O2eM6/G3HQ6r0rN+804mF9u2X+4tSqFHNgwF6x1l5tA3AHBuxBAQAYM6S3EddtcwfcoqxYEoR3akTkFSY1IvIKkxoReYU1tSwFTSyMo+lNdxJjFKyhURxBk+4bduw34j4Bk8jtXZ+CVqy1H06fcmZ/5xi7hhZUP5v+oruqdBS8UyMirzCpEZFXmNSIyCuRkpqI9BSR34rIX0RkhYiMEZFeIvKKiKxK/16e784SEYWJOlBwP4B5qvoVEekAoDOAWwDMV9W7RWQ6gOkAbs5TP8kT9QGrmgQVpCnlF2+tddq+NWZw1tcN2p4vyv/DLeNOCj3mj2sbjdhebQNwBwaCBgXunmhOMr4/9JVTQu/URKQHgHMBPAYAqnpAVbcBmARgZvqwmQAujfiaRER5E+Xj52AADQD+n4i8LyKPikgXAFWqWpc+ZguAqqCTReQaEakVkdqGxobc9JqI6CiiJLX2AM4E8LCqjgKwG6mPmn+lqUkvgatNquoMVa1R1ZrKisps+0tEdEyhK9+KyHEA3lbVQen480gltSEAzlfVOhGpBvCaqg4LuVYDgPUAKgA0HuvYBGptfW5t/QVaX59bW3+B1tfnlv0dqKqhd0ahAwWqukVENojIMFVdCWAcgOXpX1cCuDv9++wI16oEABGpjbIsb5K0tj63tv4Cra/Pra2/QOvrc5z+Rh39vA7Ar9Mjnx8B+GekPro+IyJTkLr7+momL0xElA+RkpqqLgIQlC3H5bY7RETZKdYTBTOK9LrZaG19bm39BVpfn1tbf4HW1+eM+1vQLfKIiPKNz34SkVeY1IjIKwVPaiIyXkRWisjq9DOjiSIij4tIvYgsbdGW6If3RWSAiLwqIstFZJmITE23J7LfItJJRP4sIovT/f23dPtgEXkn/d54Oj3anhgiUpJ+qmZuOk56f9eJyBIRWSQitem2RL4njsjF4hkFTWoiUgLgQQAXAxgO4HIRGV7IPkTwBIDxVtt0pB7eHwpgPqwnKhLgEIAbVHU4gLMBXJv+d01qv/cDuEBVRwAYCWC8iJwN4McA7lXVIQCaAEwpYh+DTAXQ8snrpPcXAL6gqiNbzPVK6nviiCOLZ5wMYARS/96Z9VlVC/YLwBgAf2gR/wDADwrZh4j9HARgaYt4JYDq9NfVAFYWu48h/Z8N4MLW0G+kVnx5D8DfIDVzvH3Qe6XYvwD0T39DXQBgLgBJcn/TfVoHoMJqS+x7AkAPAGuRHsCM2+dCf/zsB2BDi3hjui3pIj28nwQiMgjAKADvIMH9Tn+UWwSgHsArANYA2KaqRxbBT9p74z4ANwE4nI57IwirLUoAAAHnSURBVNn9BVLPY78sIgtF5Jp0W2LfE8hy8YwjOFCQIU39uEjkPBgR6QrgOQDTVHVHyz9LWr9VtVlVRyJ1B3QWgJOL3KWjEpFLANSr6sJi9yVDY1X1TKTKPdeKyLkt/zBp7wlkuXjGEYVOapsADGgR90+3Jd3W9EP7SP9eX+T+OESkFKmE9mtVfT7dnPh+a2ptvleR+vjWU0SOPOWSpPfGOQC+KCLrAMxC6iPo/UhufwEAqrop/Xs9gBeQ+uGR5PfERgAbVfWddPxbpJJcRn0udFJ7F8DQ9KhRBwCTAcwpcB/imIPUQ/tAxIf3C0lSy5g+BmCFqv6sxR8lst8iUikiPdNflyFV/1uBVHL7SvqwxPRXVX+gqv01tVLNZAB/UtV/QkL7CwAi0kVEuh35GsBFAJYioe8JILV4BoANInJktZ8ji2dk1uciFAMnAPgQqRrKrcUuTgb07ykAdQAOIvWTYwpS9ZP5AFYB+COAXsXup9XnsUjdkn8AYFH614Sk9hvAGQDeT/d3KYAfpdtPAPBnAKsBPAugY7H7GtD38wHMTXp/031bnP617Mj3WlLfEy36PRJAbfq98TsA5Zn2mY9JEZFXOFBARF5hUiMirzCpEZFXmNSIyCtMakTkFSY1IvIKkxoReeX/A3TyQWBNngu9AAAAAElFTkSuQmCC\n",
            "text/plain": [
              "<Figure size 360x360 with 1 Axes>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "qWaTHLTH5rGw"
      },
      "source": [
        "## Look at top co-evolving residue pairs"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "9A2yeOJ8uPNM",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 359
        },
        "outputId": "f397715a-8fa0-40f2-d239-3c5ab1ab8c08"
      },
      "source": [
        "######################################################################################\n",
        "# WARNING - WARNING - WARNING\n",
        "######################################################################################\n",
        "# - the i,j index starts at 0 (zero)\n",
        "# - the \"first\" position = 0\n",
        "# - often in biology first position of a sequence is 1\n",
        "#   for this index use i_aa and j_aa!\n",
        "\n",
        "# adding amino acid to index\n",
        "mtx[\"i_aa\"] = np.array([alphabet[msa['msa_ori'][0][i]]+\"_\"+str(i+1) for i in mtx[\"i\"]])\n",
        "mtx[\"j_aa\"] = np.array([alphabet[msa['msa_ori'][0][j]]+\"_\"+str(j+1) for j in mtx[\"j\"]])\n",
        "\n",
        "# load mtx into pandas dataframe\n",
        "pd_mtx = pd.DataFrame(mtx,columns=[\"i\",\"j\",\"apc\",\"zscore\",\"i_aa\",\"j_aa\"])\n",
        "\n",
        "# get contacts with sequence seperation > 5\n",
        "# sort by zscore, show top 10\n",
        "top = pd_mtx.loc[pd_mtx['j'] - pd_mtx['i'] > 5].sort_values(\"zscore\",ascending=False)\n",
        "top.head(10)"
      ],
      "execution_count": 20,
      "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>i</th>\n",
              "      <th>j</th>\n",
              "      <th>apc</th>\n",
              "      <th>zscore</th>\n",
              "      <th>i_aa</th>\n",
              "      <th>j_aa</th>\n",
              "    </tr>\n",
              "  </thead>\n",
              "  <tbody>\n",
              "    <tr>\n",
              "      <th>1056</th>\n",
              "      <td>21</td>\n",
              "      <td>28</td>\n",
              "      <td>0.876234</td>\n",
              "      <td>2.976396</td>\n",
              "      <td>K_22</td>\n",
              "      <td>E_29</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>881</th>\n",
              "      <td>16</td>\n",
              "      <td>58</td>\n",
              "      <td>0.846978</td>\n",
              "      <td>2.960233</td>\n",
              "      <td>R_17</td>\n",
              "      <td>A_59</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>142</th>\n",
              "      <td>2</td>\n",
              "      <td>26</td>\n",
              "      <td>0.817712</td>\n",
              "      <td>2.942693</td>\n",
              "      <td>A_3</td>\n",
              "      <td>M_27</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>999</th>\n",
              "      <td>19</td>\n",
              "      <td>50</td>\n",
              "      <td>0.784957</td>\n",
              "      <td>2.921262</td>\n",
              "      <td>I_20</td>\n",
              "      <td>G_51</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>820</th>\n",
              "      <td>15</td>\n",
              "      <td>41</td>\n",
              "      <td>0.704959</td>\n",
              "      <td>2.859523</td>\n",
              "      <td>K_16</td>\n",
              "      <td>I_42</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>488</th>\n",
              "      <td>8</td>\n",
              "      <td>45</td>\n",
              "      <td>0.663213</td>\n",
              "      <td>2.820979</td>\n",
              "      <td>E_9</td>\n",
              "      <td>P_46</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>251</th>\n",
              "      <td>4</td>\n",
              "      <td>22</td>\n",
              "      <td>0.658282</td>\n",
              "      <td>2.816094</td>\n",
              "      <td>I_5</td>\n",
              "      <td>V_23</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>1278</th>\n",
              "      <td>27</td>\n",
              "      <td>37</td>\n",
              "      <td>0.599291</td>\n",
              "      <td>2.751493</td>\n",
              "      <td>M_28</td>\n",
              "      <td>V_38</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>934</th>\n",
              "      <td>18</td>\n",
              "      <td>26</td>\n",
              "      <td>0.583831</td>\n",
              "      <td>2.732487</td>\n",
              "      <td>V_19</td>\n",
              "      <td>M_27</td>\n",
              "    </tr>\n",
              "    <tr>\n",
              "      <th>326</th>\n",
              "      <td>5</td>\n",
              "      <td>42</td>\n",
              "      <td>0.553874</td>\n",
              "      <td>2.692869</td>\n",
              "      <td>Y_6</td>\n",
              "      <td>H_43</td>\n",
              "    </tr>\n",
              "  </tbody>\n",
              "</table>\n",
              "</div>"
            ],
            "text/plain": [
              "       i   j       apc    zscore  i_aa  j_aa\n",
              "1056  21  28  0.876234  2.976396  K_22  E_29\n",
              "881   16  58  0.846978  2.960233  R_17  A_59\n",
              "142    2  26  0.817712  2.942693   A_3  M_27\n",
              "999   19  50  0.784957  2.921262  I_20  G_51\n",
              "820   15  41  0.704959  2.859523  K_16  I_42\n",
              "488    8  45  0.663213  2.820979   E_9  P_46\n",
              "251    4  22  0.658282  2.816094   I_5  V_23\n",
              "1278  27  37  0.599291  2.751493  M_28  V_38\n",
              "934   18  26  0.583831  2.732487  V_19  M_27\n",
              "326    5  42  0.553874  2.692869   Y_6  H_43"
            ]
          },
          "execution_count": null,
          "metadata": {},
          "output_type": "execute_result"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "wqG93dC12CKx"
      },
      "source": [
        "## Explore the MRF"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "k6BsheyNx3ID",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 333
        },
        "outputId": "f5df2553-66d7-48e4-e231-cd83b237b756"
      },
      "source": [
        "def plot_v(mrf):  \n",
        "  al_a = list(alphabet)\n",
        "  v = mrf[\"v\"].T\n",
        "  mx = np.max((v.max(),np.abs(v.min())))\n",
        "  plt.figure(figsize=(v.shape[1]/4,states/4))\n",
        "  plt.imshow(-v,cmap='bwr',vmin=-mx,vmax=mx)\n",
        "  plt.xticks(np.arange(v.shape[1]))\n",
        "  plt.yticks(np.arange(0,states))\n",
        "  plt.grid(False)\n",
        "  ax = plt.gca()\n",
        "  ax.xaxis.set_major_formatter(plt.FuncFormatter(lambda x,y: mrf[\"v_idx\"][x])) \n",
        "  ax.yaxis.set_major_formatter(plt.FuncFormatter(lambda x,y: al_a[x]))\n",
        "  \n",
        "plot_v(mrf)"
      ],
      "execution_count": 21,
      "outputs": [
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1wAAAE8CAYAAADUjrC3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deZhcZZ3+//uTPaFJ0iE7CYQgAUwgMbRhGXZ+CCKKsgg4ym7cwF1ndC4VUFFHnRnBEYhsKl8lKoqAIqtsgoEOhCTsZCEJCAlko0kgJHl+f1RlbDp1Tp3Phz6BLt+v68pFL3VzP1399DnnqVN1ylJKAgAAAAB0vm5v9gAAAAAAoFGx4AIAAACAkrDgAgAAAICSsOACAAAAgJKw4AIAAACAkvQou2Bwc3MaM3KkP9jDP7QVbT39PZKa+29wZ1a1dQ91NTWFYlq3zp9ZuzbWtXGjP9O//5br2uD/db0hbW3+TK9enT+OLNH7I3qB0t69/ZlVq7Zc1/r1sa6tt47l1qzxZ/r0iXWtXOnPROdi9H4cPNifiWzforkteX8EdmOSpFdfjeW6BR5CjWQkaVCfLTfxX1kXG2T3wG46ut/s3+Tfma19NfZzRbfdkXkV7Yrc9wN6Be/8V15xRxauag5VjRgRimnZMn9m4MBYV8/AobBZrCvq5Zf9meh+M5KL7Ncl6bHHZr6QUhpS63ulL7jGjByp1unT/cFm/x/D1X/b1t8j6ZiDV7gz1/819se6336hmBYv9mfmzIl1RRYYhx++5bpWr451Rd1zjz8zalSsK3LwE70/oge5Y8f6M3/6U6xrzBh/JrIokaT994/lHnjAn9l551jXddf5M9G5uMK/WZQknXqqP/P007GuZ57xZyKP/0nS8uX+zNChsa5582K5yAMU0QcBTxw30x8aNy7U9fizsUdDIg8Ezp4dqtJh+7zkzsxZGPu5otvup57yZyIPikqxefXeMcGDlscec0dO/dNxoaqvfjUU0yWX+DNHHhnrGj7cn4kuZqILtXvv9Wd23TXWFdnfzpoV65oyxTL3ZjylEAAAAABKwoILAAAAAEriXnCZ2fvNLJnZLmUMCAAAAAAaReQM14mS7q7+FwAAAACQwbXgMrMmSftKOl3SCaWMCAAAAAAahPcM11GS/pxSekLSi2a2R60bmdlUM2s1s9Zl0UtdAQAAAEAX511wnSjpqurHVynjaYUppWkppZaUUsuQwOXdAQAAAKARFH4fLjMbJOlgSbuZWZLUXVIysy+lFH1rPAAAAABoXJ4zXMdK+kVKafuU0piU0mhJCyQF38oXAAAAABqbZ8F1oqTfd/ja1eJqhQAAAABQU+GnFKaUDqrxtfM7dzgAAAAA0Dgi78MFAAAAACjAyr7eRcv48al1+nR/cJtt/Jm5c/0ZSVevPtSdOWboXaGuJ4fHXvK2020XuzOrTvhYqGvRIn9mt/l/CHX92z1HuTPnnhuq0pIlsdz69f7MmjWxrltu8Wf69Il17bBDLPcv/+LPNF/7s1DXxa+c7M5su22oSkcOiP1Na+NGd2T6cweEqo5/31p/6Kqr6t+mlg9/OBT76RU93ZmPnvRqqEvTprkjz3/wrFDVsKt+5M7cMekzoa5+/UIxjRnjzzz/fKxrwqsz3ZkL76v5TjJ1fWLyjFBOt93mz0yZEqp6dOQh7syui24MdT049LBQLrIvmz07VKXTP+TfVt14Z99Q1157+TMDbrk61PXQ244J5SY+5e9bsmesa948f2affUJVWr48lhs26DV35qVX/PsWSbrjDn9mt91CVRozxmamlFpqfY8zXAAAAABQEhZcAAAAAFASFlwAAAAAUJLCVymUJDPbIGlONbdA0kdSSivLGBgAAAAAdHXeM1xrU0qTUkoTJC2X9KkSxgQAAAAADeGNPKXwXknB64EBAAAAQOMLLbjMrLukQyRdm/H9qWbWamaty1aseCPjAwAAAIAuy7vg6mtmsyQ9J2mYpJtr3SilNC2l1JJSahnS3PxGxwgAAAAAXVLoNVyStpdk4jVcAAAAAJAp9JTClNIaSZ+W9AUzc13pEAAAAAD+WYQvmpFSelDSbEkndt5wAAAAAKBxuM5OpZSaOnz+3s4dDgAAAAA0jrfs0wEfWjrCnZmouaGuYcMCoRSq0uzZsdxOp5zizgzQq6Gu5ube/tDb3x7q+uI+/sySJaEq7Tgqdn8sWea/P3YevirU9Y4PtflDw4eHuuYt7B7KNd92dSgXcdpp/kxra7BsfSx2d/cD3Jlx42JdWr48GAxYH7tDJk/u6c7cPzuwzZH0zkDmpptCVXr/aZ9xZwbOj3VNvOacUG7eh7/hznSPbQZCPnHE07Hgc8HCT3/an/nb30JVixb5My9udVioa98HLg3llr3vdHdm7NhQlR56oq87E7kPJemw1b9xZ9refVyoaxf/5k2S9KtHjnFn+twf6/rApAXuzIZuO4S6hq1bHMqpxyh3pF+/WNV73uPPLF0a68rzRt6HCwAAAACQgwUXAAAAAJSEBRcAAAAAlMS14DKzZGY/bPf5F83s7E4fFQAAAAA0AO8ZrlclHW1mg8sYDAAAAAA0Eu+Ca72kaZI+V8JYAAAAAKChRF7D9b+S/tXMBmTdwMymmlmrmbUuW7EiPjoAAAAA6MLcC66U0mpJP5eU+QYXKaVpKaWWlFLLkObmNzI+AAAAAOiyolcp/B9Jp0vaqhPHAgAAAAANJbTgSiktl/RrVRZdAAAAAIAa3sj7cP1QElcrBAAAAIAMPTw3Tik1tfv4eUn9On1EAAAAANAg3sgZLgAAAABADtcZri1p4p+/5w9Nnhzq6ts3EFoTqtKoUbGcrrjCn1myJFS14oPfdGde6rFTqOttA/2Ziy8OVWnvvXuHcs8+688cfnjmuybkGvLgXe7MvF2PDHWtCc7hLWnRIn/mvvtiXXvHNh/ad8Md7szvFx4Q6nrH0FAs5vzzQ7Glu/+bO7PXXqEqaZdT3JG2K2NVW29c5c488khsOzD8498I5S79kT9z/PGhqpDF3bYP5UbruVhhZA5PmRKqetvb/JnIbl2S9h0Ty919tz8TPKzSTTf5M6FjsaCmG68O5V488JhQbpdd/Jn580NVUjf/uZTu85+MdS1cGIo9tHy0OzN7dqhKH/qQP/PEE7GuPJzhAgAAAICSsOACAAAAgJKw4AIAAACAkhR+DZeZbZA0R1JPSesl/VzSf6eUNpY0NgAAAADo0jwXzVibUpokSWY2VNIvJfWXFHt1LwAAAAA0uNBTClNKSyVNlXSmmVnnDgkAAAAAGkP4NVwppfmSukva7ELFZjbVzFrNrHXZihVvZHwAAAAA0GWVctGMlNK0lFJLSqllSHNzGRUAAAAA8JYXXnCZ2VhJGyQt7bzhAAAAAEDjCC24zGyIpIsk/TillDp3SAAAAADQGDxXKexrZrP0j8vC/0LSf5UyKgAAAABoAIUXXCml7mUOBAAAAAAajZX9jMCW8eNT6/Tp/uAf/+jPTJ7sz0hSW5s/M3hwrOuGG2K5c87xZ849N9Y1dLMLT9a33XaxrkWL/JlPfjLWFc3tt58/89JLsa7Ro/2ZMWNiXT17hmLXPbGzO7PLLqEq7TT9W/7QwIGxsokTQ7EFo/zzY4fLvhbq+vsnv+nOLFkSqtKaNbHcAfd8xx86+uhQ14yV/rm454SXQ12XXrWVO3P68u+HusJ3fmTbHdkGS9Kxx/oz114b6zriiFDstcl7ujM977w11KV+/fyZ1atDVav2OiyUG/DK8/7QtGmhLp1xhj9z5ZWhqrVnfsmd2bAhVKWmh2fEgpENcVNTrOvuu/2Z006LdQ0aFIqt2DjAnWk+P3AcLEnPPuvPvOtdoSo79tiZKaWWWt8r5SqFAAAAAAAWXAAAAABQGhZcAAAAAFAS94LLzIab2VVmNs/MZprZn8xsXBmDAwAAAICuzHNZeJmZSfq9pJ+llE6ofm2ipGGSnuj84QEAAABA1+VacEk6SNJrKaWLNn0hpfRQ5w4JAAAAABqD9ymFEyTNLGMgAAAAANBoSrlohplNNbNWM2tdtmJFGRUAAAAA8JbnXXA9LGmPejdKKU1LKbWklFqGNDfHRgYAAAAAXZx3wXWbpN5mNnXTF8xsdzPbr3OHBQAAAABdn2vBlVJKkj4g6f+rXhb+YUnfkfRcGYMDAAAAgK7Me5VCpZSelfTBEsYCAAAAAA2llItmAAAAAAAkqzxLsDwt48en1unT/cEhQ/yZ2bP9GUnTlx/qzhw/8q5Q1x0bYy9327jRnzloyNxQ11+WTXBnJk0KVYU0K3blywUrYxdwefhhf+a660JV+ta3/Jl77ol1jRwZy71z0dX+0PvfH+q6+pru7kz059p7fexvOvLHeX+/A0JV75z0mj905ZWhrlVHnxrKPfCAP3PQvoGfS5Iuuqj+bTqYMeWsUNWeu6xyZ/6+ZkCoa8OGUEwzZvgzU6bEukYv9b9DzL3r6l5zq6a9uwV+MEkrxu3pzjQ/cGuo6+7eh7gzrwWn/UHzLw3lZkw43Z0ZMyZUpZUr/Zk1a2Jd73jqN/5Qt9j5hrVHHBPKLV3qz2z/2I2hrhs2HubO7LNPqEoDFseOMzfs6j/OfOKJUFVoDq9bF+saONBmppRaan2PM1wAAAAAUBIWXAAAAABQEhZcAAAAAFAS14LLzEaZ2R/M7Ekzm29mPzaz3mUNDgAAAAC6ssILLjMzSb+TdE1KaSdJO0nqK+k/SxobAAAAAHRpnjNcB0t6JaV0uSSllDZI+pykk8ysqYzBAQAAAEBX5llwjZf0uuvAppRWS1oo6W3tv25mU82s1cxal62IXcIbAAAAALq6Ui6akVKallJqSSm1DGmOvfcRAAAAAHR1ngXXI5Je966FZtZf0nBJj3fmoAAAAACgEXgWXLdK6mdmJ0mSmXWX9ENJP04prS1jcAAAAADQlRVecKWUkqQPSDrWzJ6U9KKkjSmlb5c1OAAAAADoylyv4UopLU4pva96WfgjJB1uZpPLGRoAAAAAdG09osGU0j2Stu/EsQAAAABAQwkvuMr2vIa5M/5ExdChwWDAypWx3FFHvOYPLewd6rr5Zn9m5MhQlXYe6n/bgLsfjl35smfPUEy/+IU/85GPxLo2bvRnnnsu1rVwYSz3zlGB0JVXhrradLI789Ofhqq096mxXMTf/hbLrV/vn8R7x6r02GOx3OzZ/sxOO8X+OCNTMTrvW1oGuDMvLIp1Rf+mH3nEn3nXu2JdWuqPDB4c7FoeizU3BfabQX37+jOf+lSs66FPx3Jz5/ozffrEuq66yp854ohY15bUd/XzodyiRf4j1Pm9Dgt1rQxsPwYsfTLUpedj98ezAya4M9deG6rShz7kzywPbnPylHJZeAAAAAAACy4AAAAAKA0LLgAAAAAoifs1XGa2QdKcdl+6KqX03c4bEgAAAAA0hshFM9amlCZ1+kgAAAAAoMHwlEIAAAAAKElkwdXXzGa1+3d8xxuY2VQzazWz1mUr/Jf9BgAAAIBGUMpTClNK0yRNk6SW8eNTZGAAAAAA0NXxlEIAAAAAKAkLLgAAAAAoSeQphX3NbFa7z/+cUvr3zhoQAAAAADQK94IrpdS9jIEAAAAAQKOJnOHaIoZd8T1/aPLkUNdBd57jDx18cKjrqBlfDeV0hH+Mi/vsFKr6yEf8malTQ1W649gr3Zmt9z8r1DWxaV4od955O7ozY8eGqmR/vN6d+dhph4W6Fj/XM5TTff7IHWNODlW9v0Hf8e/UU2O5pu5r/aHHYl1TpsRyu+zizwz4Xmy7eO97z3Nnzv9iqEp77eXPjBwZ69ptQuxaUyNGmDuzalWoSlsHMjv1WBAri/rBD/yZ4MTf49JPujN//ONPQl26MRabFNie9ukT6zr6aH9m0KBYl57zR866/ZhQ1QXvey2U+853/JnFi0NVmnN34I/6e5eHup7/jH8bLEkbX/FnPvvZUJV6L3zcnWnabudYWQ5ewwUAAAAAJWHBBQAAAAAlYcEFAAAAACVxLbjMbJiZ/dLM5pvZTDO718w+UNbgAAAAAKArK7zgMjOTdI2kO1NKY1NKe0g6QdKosgYHAAAAAF2Z5yqFB0tal1K6aNMXUkpPS7qg00cFAAAAAA3A85TC8ZIeKHJDM5tqZq1m1rpsxYrYyAAAAACgiwtfNMPM/tfMHjKz+zt+L6U0LaXUklJqGdLc/MZGCAAAAABdlGfB9bCk/3tn4ZTSpyQdImlIZw8KAAAAABqBZ8F1m6Q+ZvaJdl/r18njAQAAAICGUXjBlVJKkt4v6QAzW2Bm90n6maR/K2twAAAAANCVea5SqJTS31W5FDwAAAAAoI7wRTMAAAAAAPlcZ7i2pJc+6X+m4tZ/uznU9bMx33BnTtZdoa6Hjj8vlJuo19yZ4cNDVWpr82fu+PwfYl2HnOXOTPzN5aGul449NZQbGfgrsYsuDHVp9Gh3pO3VnrGqUSmU033+yOTJ9W9TyyOP+DMnv/TjWJkmhlKXzz/AnTlil1CVmq7/pT/ULfa42u9+F4rp6KP9mZsPim0XD231vw3keef5tzmStP01P/KHPvShUJfa+oRivXtv7c5E37llVCBzw2M7hLrePWhpKLfqk19xZwa03hrquuvEn7gz+z18Y6jr7p1PD+X2fdsqd2bOogGhrqFD/Zn7AvsWSdoxkLngwKtDXYufOyaU+8Uv/JltWmPzQ/0OjuUCBg2K5XoufNKdSb12ipUNHOiONK95JtaVgzNcAAAAAFASFlwAAAAAUBIWXAAAAABQEteCy8zaOnx+iplFXzABAAAAAA2NM1wAAAAAUBIWXAAAAABQEu8Fr/ua2ax2nw+SdG3HG5nZVElTJWm7ESPiowMAAACALsy74FqbUpq06RMzO0VSS8cbpZSmSZomSS3jxwff7AcAAAAAujaeUggAAAAAJWHBBQAAAAAlYcEFAAAAACVxvYYrpdTU4fMrJF3RieMBAAAAgIbBGS4AAAAAKIn3KoVbzCWX+DPnn39oqGvOHH9m4r/sF+p66D1fDeUu3v48d2bcuFCVDpp7gTuz7ISzQl09XgvFQlaujOVG33GlO7PkvZ8Idf3qV/7Mp2PTPuybjxzjzqybVf82NbumXOfOvPaxM0NdkyeHYjr/fH9mWL+XQl0rjj7dnbnsslCVjtgrlrNnn3Fn1q/fNtR18y7+7c6hYxeHui5+7DPuzMeeC+xcJM3RbqHcbndf6A/tF9uXHfylPdyZn/88VKXjPrdnKLfTTv7Mgw8eEuq64foN/tAtoSr97Gex3L67+A+sdluzJlbWw394uf1pp4Wqpt9+nDtzfLerQ12PPBKK6bC2QF9TU/3b1DDnsZ7uzKgv+Y8xJemZx0MxTRj/NnfG1scOGBevG+bOjG57NNSVhzNcAAAAAFASFlwAAAAAUJI3tOAys7bOGggAAAAANBrOcAEAAABASVhwAQAAAEBJWHABAAAAQElKWXCZ2VQzazWz1mUrVpRRAQAAAABveaUsuFJK01JKLSmlliHNzWVUAAAAAMBbHk8pBAAAAICSsOACAAAAgJK8oQVXSqmpswYCAAAAAI2GM1wAAAAAUJIeZResWtdX1y2Y4M4NHuzvOv10f0aSmrqvdWfOOKNvrOzMb4diS77uz4wbF6qSDj/cHVm0KFbVu7c/03zSSaGuHktDMend73ZHFj4aqxo61J/pvf7lUNeCZ7cK5f7yF3/mox8NVWnV/u91Z/psjHXtt18sd1C3O9yZtT0OCHXNfsCfWevfvEmSBg6M5TTYP4mj9/1ZZ/kzh+7z51DXgQf6J/Ha7XYLdT12fSimUNv224e6PvUpf2bUT78R6jr22HNCuXXr/JmePUNV0po17siyyYeFqibPD8Wk0850RxY8G9hJS9qhx2J35vGVw0Jdgwb5M3OGHxPqmjw8FJPu9Ef+vntsfvTzT0VdcUWoShP8h/eV3K6BHfX82MRf32tnd6Zt9K6hrjyc4QIAAACAkrDgAgAAAICSsOACAAAAgJKEFlxm1tbZAwEAAACARsMZLgAAAAAoCQsuAAAAAChJKQsuM5tqZq1m1rpq1bIyKgAAAADgLa+UBVdKaVpKqSWl1DJgwJAyKgAAAADgLY+nFAIAAABASVhwAQAAAEBJoguufma2pN2/z3fqqAAAAACgAfSIhFJKnBkDAAAAgDpYOAEAAABASSylVGpBy/jxqXX6dH/wj3/0ZyZP9mck3dt0qDuz9/q7Ql0zeu0Xyu0562J/6IQTQl2X/naAO3P64D+Eum7ud5Q7M2ZMqEp9+sRyTz/tzzz+eKzrssv8mdNOi3V1Cz7ccmr/q/2htrZQ14VrTnZndtstVKV9U+xvet5I/9/0U0+FqrTPPv7M1r+9PFY2aVIodvmsd7gzp37w5VBX5A/m+h3OClUNHOjPrFgRqlJzcyw3JHBR4EWLYl2HDprpziwYtEeoa4elM0I53XabPzNlSqiqbc9D3Jmmv94Y6tKSJaHYHW873Z2J7ieGDvVnovv2Rx7xZ94xP7AfkzRv0jGh3MaN/sxO82Pz4zerD3Nnjnvwq6GuJ089L5Tr1cufWb48VKWXA7uXnXeOdQ0dajNTSi21vscZLgAAAAAoCQsuAAAAACiJe8FlZm3tPj7CzJ4ws+07d1gAAAAA0PWFrlIoSWZ2iKTzJR2WUgq8ygUAAAAAGltowWVm+0v6qaQjUkrzOndIAAAAANAYIguu3pKukXRgSumxTh4PAAAAADSMyEUzXpN0j6TM64ua2VQzazWz1mXR6+ICAAAAQBcXWXBtlPRBSVPMrOaF+1NK01JKLSmlliHRNxQBAAAAgC4u9BqulNIaM3uPpLvM7PmU0qWdPC4AAAAA6PLCVylMKS03s8Ml3Wlmy1JK13biuAAAAACgy3MvuFJKTe0+Xixph04dEQAAAAA0iMhruAAAAAAABVhKqdSC3XdvSddd1+rOrVvn75o/35+RpMPe5b8Ppv/aQl3Hv29tKPejaX3dmYMPDlVpt9svcGfu3+usUNc7x63yh373u1CX9t03lpsxwx2Zt/eHQ1Vjx/oz9l8/DHXppJNCsW9dPMSd2X77UJU+8vS3/KGBA0NdF3Y/M5T7xNvv8IfWrAl1rT3w3e7MlVeGqnTiibFc04rF/lD//qGuu2YPcGf2e+ynoa62Ez/qzjTNeyjUpeHDY7nItvFd7wpV/WHuju7MUfvHrlr8+9tjF98aOdKfeeWVUJUOGPqoP7RoUajrRh0Wyh02N7CvOProUJd6+F+xsnbw6FBVt8Cpg1mzQlXas9v9sWDkd93UVP82NTw/yT8/Vq4MVWnnccE1xE03+TMHHhiqWvVKb3dmwFMzQ13W0jIzpdRS63uc4QIAAACAkrDgAgAAAICSsOACAAAAgJIUWnCZWTKzK9t93sPMlpnZ9eUNDQAAAAC6tqJnuF6WNMHMNl254VBJz5QzJAAAAABoDJ6nFP5J0nuqH58o6VedPxwAAAAAaByeBddVkk4wsz6SdpeUea1sM5tqZq1m1rp8+bI3OkYAAAAA6JIKL7hSSrMljVHl7Naf6tx2WkqpJaXUMmiQ/z17AAAAAKAReN+Z7lpJP5B0oKRtOn00AAAAANBAvAuuyyStTCnNMbMDSxgPAAAAADQM14IrpbRE0vkljQUAAAAAGkqhBVdKqanG126XdHsnjwcAAAAAGobnKoUAAAAAAAfva7jcem1Yq+1fmuvOTX94gjvTv787UnHuue7IXqd8I1S14pW+9W9Uw/77+zMLF4aqtFsgM25crGvx6gHuzOhYlfTss9Gk2443XRgLjg7/dH7nx54dvO3Yb7ozkfkrSer3MX9m+vRYVb9QLGbKlFCs7y8vDaROD3U1rQq+t/2VV/ozZ5wRqlqzJhQLadqwyh+66aZYWfQHGzrUn7nkklBVr/2/4w+tXBnq2mWX5lCurc2f+clPQlU64L+iByB+W221xaqkps2e4FRMYF71Pe20UNWPrhrmzgwaFKrSngcOjwUXLYrlAoYtf9Sd6Tdq11DXkmcslBt1443+ULfYOaIBkQPU3XcPdeXhDBcAAAAAlIQFFwAAAACUhAUXAAAAAJTE/RouM9sgaU67L70/pbSw00YEAAAAAA0ictGMtSmlSZ0+EgAAAABoMDylEAAAAABKEjnD1dfMZlU/XpBS+kDHG5jZVElTJWm7ESPewPAAAAAAoOsq5SmFKaVpkqZJUsv48SkyMAAAAADo6nhKIQAAAACUhAUXAAAAAJSEBRcAAAAAlMS94EopNZUxEAAAAABoNJzhAgAAAICSWErlXkTQrCVJre5c+u73/GWTJ/szknTPPf7MwQfHum64IZbbfnt/ZsmSUNXlY7/pzgwZEqrSkQsucGdmTDkr1NWvXyimt7/dn+k+7cJY2ejR/sx73hPr+vrXY7lJgfc9b2uLdS1e7M8MHBjrmjgxltu40Z8ZOTLWdeed/ky34ONqS5fGcqec4s/86EexrsjbjvTpE+tas8YdOfmBz4SqotPj05/2Z0ac/5VY2bHH+jPXXhuqmnPsOaHcbtd/xx+aMiXUpSb/k3+Wjd0zVDXk2ktDuS89dro7Ezn0kKQzjwpsu6+6KlY2Zow/E90uvv/9sdwLL7gjf984LFQ1Ys08d+aOJTuGug7YbkEop1Gj3JEVbT1DVc88489EDsUkaeBAm5lSaqn1Pc5wAQAAAEBJWHABAAAAQEncb3xsZhskzalmH5V0ckrJ/1wLAAAAAGhwkTNca1NKk1JKEyStk/TxTh4TAAAAADSEN/qUwrskva0zBgIAAAAAjSa84DKzHpLercrTCzt+b6qZtZpZq7TsjYwPAAAAALqsyIKrr5nNUuVa74skbXZ90pTStJRSS+XSiMHrhQMAAABAF+e+aIaqr+Hq9JEAAAAAQIPhsvAAAAAAUBIWXAAAAABQEveCK6XUVMZAAAAAAKDRcIYLAAAAAEpiKaVSC3bdtSVddlmrO7d3y2v+sttv92ckzRt7qDuz47N3hbp+/8J+odwHll7szszZ52OhrrY2f2bw4FCVRo3yZ849N9YVzc2a5c/07x/r2vnJ692ZmSOODHUNGhSKaYcHrnZnnt/3mFDXwoX+zMaNoSrtvT72N/3D+/x/0wm+LuQAABvcSURBVF8489VQ18VX9HZndt89VKWnnorlJgUuqfTcc7Gu6HYnokfgElOrVsW6xoyJ5e4KTOGm4HNW3jtypjszU3uEuvZYPyOU+/t2e7ozIx65NdT16MhD3Jl160JV4f3Lfff5M+PHx7oefNCfGTo01hXZxl15ZazrpJNiucjh6fEDb4yVjRvnjsxcvkOoao/ec0O5G5+Z4M5su22oSgMG+DPDh8e6evWymZUrtG+OM1wAAAAAUBIWXAAAAABQEhZcAAAAAFAS94LLzP7DzB42s9lmNsvM/E+SBgAAAIB/Aq6XAZvZ3pKOlDQ5pfSqmQ2W1KuUkQEAAABAF+e97tIISS+klF6VpJTSC50/JAAAAABoDN6nFN4kabSZPWFmPzGzA2rdyMymmlmrmbWuWLHsjY8SAAAAALog14IrpdQmaQ9JUyUtkzTdzE6pcbtpKaWWlFJLc/OQThkoAAAAAHQ17rdyTCltkHS7pNvNbI6kkyVd0bnDAgAAAICuz3WGy8x2NrOd2n1pkqSnO3dIAAAAANAYvGe4miRdYGYDJa2X9JQqTy8EAAAAAHTgWnCllGZK2qeksQAAAABAQ3G/8TEAAAAAoBj3RTO8tuq2VntvPdede3DuBHfmHe5ExaOP+jO9Ju4X6lo9PxTT3eM/5s6MaY51NTX5M336xLr6vrLCndm4MfaDPfdcKKZJk/yZuf4pL0laMulId6bHi7Gu9etjufu3O8ad2S5WpcGD/Zlhw2JdT/499jf985/7M2ec0TvU1SvwNvMvBN8tcYJ/EyxJWrkylou47z5/5mPdfhrq+v3gj7ozb397qEpr18Zyu+3mzyxYEOv6xrV7uDPbBTcE/fbZM5TbENg2jgg1SYsW+TNbbRXrWro0lhs/3p+J/j1HtlXjxsW6tt7an9lll1jXmjWxXKTvukWHhbp2eNmfWbIkVKU99t82lIukJuwY2zCmPn3dmba2UFUuznABAAAAQElYcAEAAABASVxPKTSzbSTdWv10uKQNqrwBsiRNSSmt68SxAQAAAECX5r1K4YuqvPeWzOxsSW0ppR+UMC4AAAAA6PJ4SiEAAAAAlIQFFwAAAACUpJQFl5lNNbNWM2tdtsJ/2W8AAAAAaASlLLhSStNSSi0ppZYhzcE3gwIAAACALo6nFAIAAABASVhwAQAAAEBJXJeFby+ldHYnjgMAAAAAGg5nuAAAAACgJJZSKrWgZfz41Dp9uj84frw/c8st/owk3XOPP3PwwbGuG24Ixa58+3nuzMCBoSodueACf2i77WJlBx7oz/zud7GusWNjucWL/ZmXXop1jR7tzzz+eKxr5cpYbtIkf6atLdYVue+jE3/ixFhu40Z/Zs2aWNezz/oz3YKPqy1dGstF7o/o38uIEf5Mnz6xrsjvbN26LdclSUOHuiNPH/GJUNX2L8x0Z+b23iPUNeHlGaHcvMF7ujM7Lrw11KV+/fyZ1atjXUuWxHKRbX50LvYIPIEqkpG028+/5M7MOfvqUJduvz2WixzrNDXFuu6+2x1J534zVGWrV4VyOuccf2bAgFjX7rv7M3/+c6jKpk2bmVJqqfU9znABAAAAQElYcAEAAABASVhwAQAAAEBJ6i64zOy/zeyz7T6/0cwuaff5D83s82UNEAAAAAC6qiJnuP4qaR9JMrNukgZLan9Fi30kBa46AQAAAACNrciC6x5Je1c/Hi9prqSXzKzZzHpL2lXSAyWNDwAAAAC6rLrX4EwpPWtm681sO1XOZt0raVtVFmGrJM1JKb3u2rdmNlXSVEnaLnLZXgAAAABoAEUvmnGPKoutTQuue9t9/teON04pTUsptaSUWoY0N3fWWAEAAACgSym64Nr0Oq7dVHlK4d9UOcPF67cAAAAAIIPnDNeRkpanlDaklJZLGqjKoosFFwAAAADUUHTBNUeVqxP+rcPXVqWUXuj0UQEAAABAA6h70QxJSiltkNS/w9dOKWNAAAAAANAoip7hAgAAAAA4FTrD9Wa4+6/mzrzW49BQ10Ff3McfeiD21mMPHndeKPfh+y52Z1a992OhLi2IxUL6969/m86yzTax3OLFnTuOt4qvfjUUu39uX3fmnS0p1KVFi/yZ666LdW1JW2+9xapuHnVqKHfo0u918khyfOUrodj9j/nvx3fO+mmoa4s64YRYrlcvd6R791hVxNq1W65LkroFHlKeOfCQUNfkyf6M3XRjqOuucaeHcvvd98NQLmLtZ/1/05HflySdFMj8at0xoa4TPzwqlNNrr/kzL78c6wqwlStCubaesSuRN33ta/7Q+eeHujR4sD/z+c/HuqZNy/wWZ7gAAAAAoCQsuAAAAACgJCy4AAAAAKAkrgWXVdxtZu9u97XjzOzPnT80AAAAAOjaXBfNSCklM/u4pN+Y2V+q+fMkHV7G4AAAAACgK3NfpTClNNfMrpP0b5K2kvTzlNK8Th8ZAAAAAHRx0cvCnyPpAUnrJLV0/KaZTZU0VZK2GzEiPDgAAAAA6MpCC66U0stmNl1SW0rp1RrfnyZpmiS1jB8ffAMeAAAAAOja3shVCjdW/wEAAAAAauCy8AAAAABQEhZcAAAAAFCS6EUzlFI6uxPHAQAAAAANhzNcAAAAAFASS6nciwiOG9eSLrig1Z278UZ/VyQjSbfc4s8cHnyr59tvj+UiY+zbN9Z15IIL/KGPfzzU9dIrPd2ZrX97eajrxfedGspts3GZO3PHI0NCXfPn+zOnnrQh1KV77w3FvvfXfWN9Af/6r/7MihWxri9/OZa76CJ/Zvs1j4a6vvbLXd2ZtrZQlcaMieVOOMGfeeCBWNfgwf7M2LGxrvPP92fO+fLLoa7Hl2wVyu088iV3ZtkrW4e6IvvA3/42VKUjj4zlTjvNn/nrX2NdF17ozwybFTtoOfWqw0K5y78w1x8aPTrUpZ/8xJ/pEXvS1V17fcmdGTkyVKV77onlPtLvandmyZ7HhLpuvdWf2X//UJXWr4/lItvu5v7BY51XXnFHNvSJbYN79LCZKaXN3i5L4gwXAAAAAJSGBRcAAAAAlKTwgsvM/mJmh3X42mfNLHAiHQAAAAAan+cM168kdXx2/gnVrwMAAAAAOvAsuH4r6T1m1kuSzGyMpJGS7ur8YQEAAABA11d4wZVSWi7pPknvrn7pBEm/TjUuc2hmU82s1cxaV63yX+ENAAAAABqB96IZ7Z9WmPl0wpTStJRSS0qpZcCA2OWxAQAAAKCr8y64/iDpEDObLKlfSmlmCWMCAAAAgIbgWnCllNok/UXSZeJiGQAAAACQK/I+XL+SNFEsuAAAAAAgVw9vIKV0jSQrYSwAAAAA0FAiZ7gAAAAAAAVYjau6d6rhw1vSySe3unPnnuvveuwxf0aSJo5Z5c7MfGpAqGuPMS+Gcvrtb/2ZJUtCVfNO+aY7s+PcP4S6tP/+7sgKNYeqml94MpTTjBn+zEsvxbpGj/ZnHn881rVyZSw3aZI/09YW61q82J8ZODDWNXFiLLdxoz9zyy2xrjFj3JF02umhKvvP74VyofvjjDNCVfcv9F8F952zfhrq0po1/sy6dVuuS5KGDvVnFi2KdR17rD9z7bWxriOOiOUCfy/z2oaFqnZc6b+G2JJhe4S6Rt14aSgX2uZH52IP9xOoYhkp9HtWt+D5hkiXpLad/b/r6C561Lb+4/p7/xZ74trIkaGYXnjBn3nqqVhX//7+zOrVsa4TTrCZKaWWWt/jDBcAAAAAlIQFFwAAAACUhAUXAAAAAJSEBRcAAAAAlIQFFwAAAACUhAUXAAAAAJSklAWXmU01s1Yza127dlkZFQAAAADwllfKgiulNC2l1JJSaunb1/8+KQAAAADQCMILLjP7lJnNqv4LvvUZAAAAADSu4Nt6Syml/5X0v504FgAAAABoKFw0AwAAAABKwoILAAAAAErCggsAAAAASmIppXILzJZJejrj24MlveD8X0YyWzrXqF3RXKN2RXN0vXm5Ru2K5hq1K5pr1K5ojq43L9eoXdFco3ZFc43aFc29Vbq2TynVvjx7SulN+yepdUtktnSuUbu6whi5P/45urrCGLk/uD/eKl1dYYyN2tUVxsj9wf3xVunqCmOMdvGUQgAAAAAoCQsuAAAAACjJm73gmraFMls616hd0VyjdkVzdL15uUbtiuYatSuaa9SuaI6uNy/XqF3RXKN2RXON2hXNveW7Sr9oBgAAAAD8s3qzz3ABAAAAQMNiwQUAAAAAZYlc2vCN/pN0uKTHJT0l6d8LZi6TtFTSXGfXaEl/kfSIpIclfaZApo+k+yQ9VM2c4+zsLulBSdcXvP1CSXMkzZLjcpOSBkr6raTHJD0qae8CmZ2rPZv+rZb02QK5z1Xvi7mSfiWpT8ExfqaaeTivp9bvV9IgSTdLerL63+YCmeOqXRsltTi6vl+9H2dL+r2kgQUy36zefpakmySN9MxbSV+QlCQNLtB1tqRn2v3ejijaJems6s/2sKT/LHh/TG/XtVDSrAKZSZL+tmkeS5pSsGuipHurfwPXSepf5G+4wPzIymXOkZxMvfmRlcucI1mZAvMjqyt3juT1Zc2RnK568yMrlzlHcjL15kfN7bWkHSTNUGU/M11Sr4K5M6uZWvd9Vub/qbJPm6vKHO9ZMHdp9WuzVdmWN9XLtPv++ZLaavyNZXVdIWlBu9/bpII5k/RtSU+osp/5dIHMXe16npV0TcGuQyQ9UM3dLeltBTIHVzNzJf1MUo+O90n1dq/bL9ebHxmZzLlRJ5c7PzIymXMjL1dvfuT05c6PjEzm3KiTy50fGZnMuVEnV3d+qMYxmOrvX2plihx/1Mrl7l9ycrnHILUy7b5Xc/+S03W28vcvNbtU//ijVlfu/iUnl3sMkpHJ3b9Ub7PZsXa9+VHzd1/vBp39T5U/hnmSxkrqpcoG5e0FcvtLmiz/gmuEpMnVj7dWZcOQ26fKRqSp+nFPVTbKezk6Py/pl/ItuDI33jm5n0k6o/pxL9X4Iy3wu3hOlTdqy7vdtqpsiPtWP/+1pFMK/P8nqLKR6yeph6RblL2R3Oz3K+k/VV2QS/p3Sd8rkNlVlUXl7cre4NXKvUvVDbGk7xXs6t/u409LuqjovFXl4PJGVd4UvONBXa2usyV90fs3Iumg6v3eu/r50KJjbPf9H0r6eoGumyS9u/rxEZJuLzjG+yUdUP34NEnf7JCp+TdcYH5k5TLnSE6m3vzIymXOkaxMgfmR1ZU7R3JymXMkb4x15kdWV+YcycnUmx81t9eqbKdOqH79IkmfKJh7h6QxqrFdzskcUf2eqfKAVNGu9vPjv9TuAcisTPXzFkm/UO0FV1bXFZKOzZkfWblTJf1cUrca86PuvlLS1ZJOKtj1hKRdq1//pKQr6mT2kbRY0rjq18+VdHrGz/e6/XK9+ZGRyZwbdXK58yMjkzk38nL15kdOX+78yMhkzo16Y8ybHxldmXMjK6fKs7jqzo9av0/V37/UyhQ5/qiVy92/5ORyj0Gy5qly9i85XWcrf/9SK1Pk+KPmGNt9f7P9S05f7jFIRiZ3/1L9+mbH2vXmR61/b8ZTCqdIeiqlND+ltE7SVZKOqhdKKd0pabm3LKX095TSA9WPX1JldbptnUxKKbVVP+1Z/ZeK9JnZKEnvkXSJd6weZjZAlYPXSyUppbQupbTS+b85RNK8lNLTBW7bQ1JfM+uhygLq2QKZXSXNSCmtSSmtl3SHpKNr3TDj93uUKhNd1f++v14mpfRoSunxvEFl5G6qjlGqPEIyqkBmdbtPt1KNOZIzb/9b0pedmVwZuU9I+m5K6dXqbZZ6+szMJH1QlYOEepkkqX/14wGqMUcycuMk3Vn9+GZJx3TIZP0N15sfNXN5cyQnU29+ZOUy50idbVPe/HBv0+rkMudIva6c+ZGVy5wjOZl68yNre32wKo9KSrXnR81cSunBlNLCjPswK/On6veSKmdhOs6PrNxq6f/ux756/fyomTGz7qo8Iv5lzxhr3bZg7hOSzk0pbazebmmBjKo/V39Vfg/XFOzKmx+1MhskrUspPVH9+mbzozqO1+2Xq/d37vyotS/Pmxt1crnzIyOTOTfycvXmR1aunoxM5two0pU1PzIydfcvNXLbqMD8yJC7f6klb99SJ5e7f8nJ1T0GyZC5f+lkdY8/8mTtX3LUnSM15O5fco613fPjzVhwbavKIw6bLFGBg4XOYGZjVHmEakaB23Y3s1mqPP3p5pRS3UzV/6gykTc6hpYk3WRmM81sasHMDpKWSbrczB40s0vMbCtHpySdoAITOaX0jKQfSFok6e+SVqWUbirw/58raT8z28bM+qnyiMNox/iGpZT+Xv34OUnDHNk34jRJNxS5oZl928wWS/pXSV8vmDlK0jMppYec4zrTzGab2WVm1lwwM06V38EMM7vDzN7p7NxP0vMppScL3Pazkr5fvT9+IOkrBTse1j8edDlOOXOkw99w4fnh+dsvkMmdHx1zReZI+4xnftQYY6E50iFXaI5k3B9150eHXKE50iFTd3503F6r8iyKle0OYmruZyLb+byMmfWU9BFJfy6aM7PLVZm/u0i6oEDmTEnXtpv7njF+uzo//tvMehfM7SjpeDNrNbMbzGynoveHKgcht3Y4MMzLnSHpT2a2pHo/fjcvo8ripYeZtVRvcqxqbz867pe3Uf35EdmX5+Zy5kfNTN7cyMnVnR85Y8ybH7UyuXOjTpeUPT9qZXLnRkbuBRWbH7WOwertXyLHbUVyWfuXmrk6+5fNMgX3L1ljzNu/1MoU2bfk3R95+5dauXr7l1qZevuXrGNt//FpqnMKrLP/qTLhL2n3+Uck/bhgdoycTylsl22SNFPS0c7cQFVeVzChwG2PlPST6scHqvhTCret/neoKk+x3L9ApkXSekl7Vj//kWqcCs3J91JlYzSswG2bJd0maYgqjypeI+nDBXtOr97vd0q6UNL/FP39qrJDbP/9FUXnhHJO6dfJ/Ycqz6E2z/xT5Q+75mv92udUOTs4Q9KA6ucLVfuUfsf7YpgqTwHtpspz5i8reB/OVWVHbaqcXV7g+dmqv7MvFOw6X9Ix1Y8/KOmWgrldVHkqwExJ35D0YkbudX/DReZHrVyROZKTyZwfebm8OdI+U3R+ZNwfRedIx1zdOZJzf2TOj4yuunOkRqbQ/KjedtP2el9Vnkmx6euja83vGrkJ7b6Wed/nZH6qnG1cTq67pJ9IOrVOZn9VXr+y6elHmU8Z69ilylM2TVJvVR6R3expOhm5tk2/4+ocvcvxc92w6fddsOt3+sc+7Utqd6yQk9lbldcE3SfpW9r89YSb7ZclDc6bH7UyHf6fNedGgdxm86NApubcyPi5RtabH1l9efMjJ5M7Nwr8bJvNj5yu3LmRk8udH9XbbHYMpjr7l1qZdt+7Xdn7lrxc3vFH7nGiauxfMn6uIscftXK5+5eMTJF9S979kXf8Uasvd/+SkcndvyjjWLve/Kg55no36Ox/1cl/Y4dJ8pWC2TEKLLhUWSTcKOnzwTF/XXVeP1O93XdUeaRsoSor3jWSrnR2nV2wa7ikhe0+30/SHx09R0m6qeBtj5N0abvPT1J1w+b82c6T9Mmiv19VXmQ8ovrxCEmPF50TCiy4JJ2iyosn+3nnn6Ttcr73fzlJu6ny6OzC6r/1qpw5HO7oKvw9VR5JPajd5/MkDSl4f/SQ9LykUQW7Vukf7+1nklYHxj9O0n01vr7Z33DB+ZH5t581R7IyBeZH7nam1hzpmHHMj3pdNe/jjPsxd47k3B/15ketrtw5UuDnqjk/Otzm66ockL2gfxx4vm6/k5P7YrvPF6rOa2vbZ1TZWV+j6utZPF3Vr+2vnAfpqplvqLJv2TQ/NqrdwsHRdWBeV/ucKi8U36Hd72xVwftjsKQXVeACS+1+Z/M6/L084vy53iXp1x2+Vmu//P/y5kdG5sp23685N/JyWfOjXlfW3MjIrag3Pwr2vW5+ZGXqzY0690fN+ZGR+WO9uVHw59psftT4HZ6tyryvu3/pmGn3+e3KOf6olVOd/UteX7v7JO9BpbMlfU0F9i8FusYU6PqiCh5/ZNwfufuXjL5CxyA5P9dm+xdlHGt75sf/ZevdoLP/Ve/E+aqcptt00YzxBbO5v+SMjKnyos7cRx07ZIaoegEKVZ47fZekI529B6rAGS5Vnne7dbuP75F0eMGOuyTt3G7yfN8xvquU8Whqjdvuqcpp137V+/Nnks4qmB1a/e92qmycMy/s0fH3q8rz0Nu/KLHWFW5qzgk5F1yqXDnzkTobg46Zndp9fJak33rnrYqf4RrR7uPPSbqq4Bg/rspz7KXKxmSxCp7hqt4ndzjuj0clHVj9+BBJMwvmNs2Rbqr8rZ7W4fY1/4brzY+sXN4cyenKnR85ucw5Um98WfMjpyt3juTkMudI3hjz5kdOV+YcycnUmx81t9eSfqPXXxThk0Vyde77rK4zVNl29824P2rl3qvqRYSqP/sPJP2g6PiqX691BiNrjCPadf2PKq+tKJL77qb7XJX92v1FxlidVz9z3B9HqrII2nSBg9MlXV0gs2l+9JZ0q6SDc/6eDtQ/znzkzo9amby5Uacrd350zFR/R5lzo8gYs+ZHzhhz50dGJnNu1Btj3vyocX/0yJsbdcaYOz+UcQymnP1LVqbd929X7Qfzsrrq7V+ycnn7l7rHlrXmcU5X5v4lJ5N7/JE3RuXvX7L68vYvWZnc/Uv1e5sda+fNj8x5We8GZfxT5bU8T6iy2v2PgplfqfL6oddUeRSj5pWIauT2VeV5m5sunVnzktodMrurclnR2aqcEs186kXO/+NAFVtwjVVl0bnpUreF7o9qdpIql76crcojZ3UvS9lusr2o6mnlgplzVFkwzVXl6ke9C+buUmVD8pCkQzy/X1WeZ3+rKpfdvEXSoAKZD1Q/flWVR0c2e1Q7I/eUKhuDTXOk49V+amWurt4fs1W5nOi23nmr2hu8Wl2/UOWypbMlXat2G786uV6qPBI5V5VL4252MJI1RlWuWvVxx+9rX1VOyz+kytMW9iiY+4wq24MnVNmBd3zKQc2/4QLzIyuXOUdyMvXmR1Yuc45kZQrMj6yu3DmSk8ucI3ljrDM/sroy50hOpt78qLm9VmXbel/1d/cbddhm5eQ+XZ0f61V50fUlBTLrVdmfbRp3x6s2bpZTZQf/1+rvbK4qZ1361+vq8P+tteDKGuNt7bquVIfLjOfkBqryiO4cVR6Bn1hkjKocdNZ88DCn6wPVnoeq+bEFMt9X5UDrcdV5ixO9/kA8d35kZDLnRp1c7vzomKk3N/K66s2PnDHmzo+MTObcqDfGvPmR0ZU5N+rkcueHMo7BlLN/ycnkHn/k5OrtX7JyefuXuseWqr1/yerK3L/kZHKPP/LGqPz9S1Zf3v4lK5O7f6neZrNj7bz5kfVv06OYAAAAAIBO9mZcpRAAAAAA/imw4AIAAACAkrDgAgAAAICSsOACAAAAgJKw4AIAAACAkrDgAgAAAICSsOACAAAAgJL8/5b8KdsanQopAAAAAElFTkSuQmCC\n",
            "text/plain": [
              "<Figure size 1098x378 with 1 Axes>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "OER2wALRvTkK",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 681
        },
        "outputId": "d4df0720-6edf-48c5-8b60-b6d13d387c31"
      },
      "source": [
        "def plot_w(mrf,i,j):\n",
        "  \n",
        "  n = int(np.where((mrf[\"w_idx\"][:,0] == i)&(mrf[\"w_idx\"][:,1] == j))[0])\n",
        "  w = mrf[\"w\"][n]\n",
        "  \n",
        "  mx = np.max((w.max(),np.abs(w.min())))\n",
        "  plt.figure(figsize=(states/4,states/4))\n",
        "  plt.imshow(-w,cmap='bwr',vmin=-mx,vmax=mx)\n",
        "  plt.xticks(np.arange(0,states))\n",
        "  plt.yticks(np.arange(0,states))\n",
        "  plt.grid(False)\n",
        "  \n",
        "  ax = plt.gca()\n",
        "  al_a = list(alphabet)\n",
        "  ax.xaxis.set_major_formatter(plt.FuncFormatter(lambda x,y: al_a[x])) \n",
        "  ax.yaxis.set_major_formatter(plt.FuncFormatter(lambda x,y: al_a[x]))\n",
        "  plt.title(f\"W for positions {i} and {j}\")\n",
        "  plt.show()\n",
        "\n",
        "for n in range(2):\n",
        "  i = int(top.iloc[n][\"i\"])\n",
        "  j = int(top.iloc[n][\"j\"])\n",
        "  plot_w(mrf,i,j)"
      ],
      "execution_count": 22,
      "outputs": [
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAATwAAAFMCAYAAABWJwDHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deZRdZZnv8e8vCZWBkJChSBgSaFBmMEK0lSWD4ISADQ4YtBuxvZe2l+gSsVtb+9poX1vbJdpXQV20CojK0NgKChf1iihDUIMCAVTmWaSSACFkrOS5f+xdcDhUpc77Zu9TVezfZ62zkrPPfs/7nF2nnnrfPTxbEYGZWROMG+kAzMy6xQnPzBrDCc/MGsMJz8wawwnPzBrDCc/MGsMJb4yQdJykByWtkvTSkY5nMJI+Junrm3n9nZJ+0s2YRitJJ0m6dqTjaBonvC6Q9E+S/m/bsjuHWLZoiLf5PHBKREyNiN/VFeuWiIh/i4j/ASBpF0khaULL69+JiNd1Kx5JEyV9Q9L9kp6SdJOkI1te75F0iaT7ylgP61ZsmyNpd0mXSuqTtELSjyXt0fK6JP1vSQ9LelLS1ZL2GcmYxwonvO74JXCQpPEAkrYHtgJe2rbsReW6g9kZuC2n84E+GmgC8CBwKDAd+GfgYkm7tKxzLfDXwKPdDm4ztgUuA/YA5gC/Bi5tef1twN8CBwMzgcXA+V2OcWyKCD9qfgA9wGrgwPL58cA5wC/alt01SNuJwCoggKeBu8vlewFXA09QJMI3tbQ5F/gqcEXZ5jWDvO/VwGcofplWUvxCzWx5/U3l+z5RrrtXy2sfAR4GngL+CBxRLj8d+Hb5/wfKmFeVj1cCJwHXtrzPQcBvgCfLfw9qi+9fgevKfn4CzC5fmwR8G1hexvcbYE6HP4tbgLcMsvwh4LBh2r4b+H0Zzz3A37W8dlj5HqcBjwF/At7d8vosiiS2stzm/9q6LYbpd2a5LWe1bP+LW17fB1g70t/zsfDwCK8LImI98CvgkHLRIcA1FKOL1mXPG91FxLqImFo+fUlE7CZpK+CHFElgO+D9wHdapz3AO4BPA9uU/QzmRIqRwvZAP/AlKKZUwAXAB4FeisT5w3IKuAdwCvCyiNgGeD1w3yDvPfC5to1iGr649UVJM4HLyz5nAV8ALpc0q+0zvLv8jD3Ah8vl76IYsc0r274XWDPEZ2ztcw6wO5kjZYpEdjQwrYzri5IOaHl9bhnXjsB7gLMkzShfOwtYS7Gt/7Z8dOoQ4NGIWF4+vxDYrZz6bkWxPa7M+0jN4oTXPb/g2SRwMEXCu6Zt2S86fK9XAFOBz0bE+oi4CvgRcELLOpdGxHURsSki1g7xPudHxK0R8TTwv4Djy+nv24HLI+KnEbGBYv/hZIoR2UaKUefekraKiPsi4u4O4251FHBnRJwfEf0RcQHwB+CYlnXOiYg7ImINcDGwoFy+gSLRvSgiNkbEjRGxcnOdlYnhO8B5EfGHjHiJiMsj4u4o/ILiD87BLatsAD4VERsi4gqKke0e5TZ9C/CJiHg6Im4FzuukT0k7USTLD7Us/hPFH7E/UiT6twGn5nympnHC655fAq8qRza9EXEncD3Fvr2ZwL4Mvf+u3Q7AgxGxqWXZ/RQjiwEPdvA+revcT7FfcXb5/vcPvFD28yCwY0TcRTHyOx14TNKFknboMO72z3B/27L2z9C6X201RZKHYn/Vj4ELJT0i6XNlQhuUpHFlm/UUo9Msko6UdEN5IOEJ4I0U22vA8ojoHyTmXp7dnzig/bMP1l8vRVL9SvkHYcAngJdRjHAnAZ8ErpI0JeNjNYoTXvcsppju/E+K/VKUo5JHymWPRMS9Hb7XI8C88hd5wHyK/WoDOimDM6+t/QZgWfn+Ow+8IEnlug+XcX83Il5VrhPAvw/y3sP1/5w+WmJ4eJB1n/vGxQjqkxGxN8Wo82iK6fnzlLF/g2Ln/1vKEWsySROB71GMdudExLYUU3110LyPYpdB+/beXH8zKJLdZRHx6baXFwAXRcRD5ej4XGAGsHcnn6XJnPC6pJyWLaGYmlzT8tK15bJOR3dQ7A9cDfyjpK3K0ymOodi3k+KvJe1djgw+BVwSERsppo9HSTqiHDmdBqwDrpe0h6TDywSwlmJKtWmQ9+4rl+86RN9XALtLeoekCZLeTvEL+6Phgpb0akn7lVPFlRSJerAYoDh4sxdwTPkzaH+viZImlU97JE0qk2S7HoqpfB/QX57e0tEpNuU2/W/gdElTJO1Nsd9tqM83jWIEe11EfHSQVX4DvE3SHEnjJP0Nxej8rk7iaTInvO76BcUO+NaDCNeUyzpOeOVBkGOAIylGZF8BTszYN3U+xRHdRymmRh8o3/+PFKdqfLl8/2MoEsZ6il/6z5bLHy1j/6dBYlxNcdDkOklPSHpF2+vLKUZmp1Ecbf1H4OiIWNZB3HOBSyiS3e8ptuvzTsuQtDPwdxQjokfLk7ZXSXpny2oD+8F2pEgya3j+yJOIeKrcPhcDj1McULmsg1gHnEIxvX2UYpufs5l1j6OYsr67JeZVkgZGhf8O3AzcRHGU+lSK0esTCfE0kiJcALSJJF1NcQrJkFdGmL3QeIRnZo3hhGdmjeEprZk1hkd4ZtYYE4ZfZcvMnDk7dtppl+R2PWueTO9s8uT0NkBs1ZPcZtATFzrpK3NAvXFjepsJuT/dnM4AxqfXKMjdHtqwPr3RVkOemzxMZ937Ya9ek9lXpindPFU554edse3vu+8+li1bNmjD2hPeTjvtwuWXL0luN+/mYU/Her799ktvA6yb+7yzEIY1sSfvN3VDf94X+omMEw5mzszqivGrMv7YAEydOvw6bdb15xVymfjosBcqPN8OOReEkP+Xo79/+HXa3HhLZlLOdOABXdyltT7jj1RP+mBk4cteNuRrntKaWWM44ZlZYyQnPEnHltVh96wjIDOzuuSM8E6guDTqhOFWNDMbTZISnqSpwKsoihsOde8FM7NRKXWE91fAlRFxB7Bc0oGDrSTpZElLJC1ZsaJvi4M0M6tCasI7gWdLEF3IENPaiDg7IhZGxMKZM3u3JD4zs8p0fIJRWZX3cGA/SQGMB0LSP4SvTzOzMSBlhPdWinsg7BwRu0TEPOBenlvT38xs1EpJeCcA329b9j18tNbMxoiOp7QR8epBln2p2nDMzOrjKy3MrDFqLx7QM66feVOWD79iu81cADykSZOGX2cQjzyS3uYvdsnqKtvtt6e3OfSgrBt0ZRUBAGBc+t/P3CIMLOvk1hdt5m/2RmHVyyg6cODuT2V19Xj/NlntuunJtROT20zP/X4MwSM8M2sMJzwzawwnPDNrjKSdDJI2AkvLdvcCf+N7YZrZWJE6wlsTEQsiYl9gBfC+GmIyM6vFlkxpF1Pcrd3MbEzISniSxgNHAJcN8foz1VL6lmeckmJmVoPUhDdZ0k3Ao8Ac4KeDrdRaLaV31qwtjdHMrBJZ+/CAnQHhfXhmNoZkTWkjYjXwAeA0SbVfrWFmVoXsgxYR8TvgFlwtxczGiKTRWURMbXt+TLXhmJnVp/7p6LhxWXcP509/Sm6yYY990/sB5maEl+uuu/La5Vz3/uTqvLvYT5+SWXQgo3hAtgULutfXGDBj3JOZLadVGsfmTJ+0LqNVtb+cvrTMzBrDCc/MGsMJz8waI/VG3CHpjJbnH5Z0euVRmZnVIHWEtw54s6TZdQRjZlan1ITXD5wNnFpDLGZmtcrZh3cW8E5J04da4TnFA3LuPWBmVoPkhBcRK4FvUVxaNtQ6zxYPmO3Zr5mNDrlHaf8DeA+wdYWxmJnVKrd4wArgYoqkZ2Y2JmzJeXhnAJ6vmtmYkV08ICL+DEypPCIzs5r4Sgsza4zuVEuZkjEQ3Gab5CYPPJDeDcBu63+f3mjPPbP6mjcvqxlTr/txeqNDDsnrrAtfi2esX5/V7NSPTkxu88UvRFZfXbV6dV67qVOHX2eELV+V/jObNbPan5lHeGbWGE54ZtYYTnhm1hgd76yRtBFYCmxFcU3tt4AvRsSmmmIzM6tUyt7pgVs0Imk74LsU9aH/pY7AzMyqlnulxWPAycApklRtSGZm9diS2zTeA4wHtmt/7TnVUvr6tiQ+M7PK1HLQ4jnVUnp76+jCzCxZdsKTtCuwEXisunDMzOqTlfAk9QJfA86MiDFw+rqZWdpR2smSbuLZ01LOB75QS1RmZjXoOOFFxPg6AzEzq1t3rhIflzFznj8/uclu5M6u8woB5Ji6dWaMr3tdtYGMFj09Wc3GRCGAHNs976SHF4yqCwHk8KVlZtYYTnhm1hhOeGbWGMkJT9JcSRdKulvSjZKukLR7HcGZmVUp6aBFed3s94HzImJRuewlwBzgjurDMzOrTupR2lcDGyLiawMLIuLmakMyM6tH6pR2X+DGOgIxM6tbLQctXC3FzEaj1IR3G3DgcCu5WoqZjUapCe8qYKKkkwcWSNpf0sHVhmVmVr2khFdWRjkOeE15WsptwGeAR+sIzsysSsnX0kbEI8DxNcRiZlYrX2lhZo1Rf7WUCFi/Pr3dihXJTZb3bJ/eD3nhbT83s/JDTme57aZMyeqqb0VeJbDe2V2shrFyZXqbadOqj6Nijz+Rd0+sGVPW5XWYWa0my6aMO7rmVFra3NtV+m5mZqOYE56ZNYYTnpk1RlLCk7STpEsl3SnpHklnSppYV3BmZlXqOOGVlVL+G/hBRLwYeDEwGfhcTbGZmVUqZYR3OLA2Is4BiIiNwKnAiZKm1hGcmVmVUhLePrRVSomIlcB9wItal7t4gJmNRrUctHDxADMbjVIS3u20VUqRNA2YC/yxyqDMzOqQkvB+BkyRdCKApPHAGcCZEbGmjuDMzKrUccJrqZTyVkl3AsuBTRHx6bqCMzOrUmp5qAcj4k3laSlvBN4g6YB6QjMzq1Z28YCIuB7YucJYzMxqVX+1lLVr4Y70Ozhq//2S28TazIoRt9+e3mbugu71BfC5jPO7zz03q6sJE/Kqpax6Oq/SR46rr5me3Oaww/L6mrp1XhWYp1alb48ZE57K6ot7Hsprt+eeee1y5FS42XbbSkPwtbRm1hhOeGbWGE54ZtYYyfvwJG0ElrYsujAiPltdSGZm9cg5aLEmIjL32JuZjRxPac2sMXIS3mRJN7U83t6+wnOqpTz+eAVhmpltuVqmtBFxNnA2wMJ99unirazMzIbmKa2ZNYYTnpk1Rs6UdrKkm1qeXxkRH60qIDOzuiQnvIjIu9DSzGyE1V48oH+ryfTNzSgE8MCD6W165iW3Abhn2kuT2+xG3rGYJ3dN7wtg+kczBtHj8vZYzJi2MatdVn+bNmV19drXpv/dndjT3eNn20zN6G/F+qy+7u7ZK6td7vc4y5Qp3etrCN6HZ2aN4YRnZo3hhGdmjZGU8CTNkfRdSfdIulHSYknH1RWcmVmVOk54kgT8APhlROwaEQcCi4Cd6grOzKxKKUdpDwfWR8TXBhZExP3AlyuPysysBilT2n2A33ayYmvxgOXL+/IiMzOrWPZBC0lnSbpZ0m/aX4uIsyNiYUQsnDWrd8siNDOrSErCuw145h60EfE+4AjAGc3MxoSUhHcVMEnS37csG/lTp83MOtRxwouIAI4FDpV0r6RfA+cBH6krODOzKiVdSxsRf6I4FcXMbMzxlRZm1hi1V0uZMAF6Z+dUZEg/n1mZlR922zWrWZbp0zKrU+y/f7WBjBaZFV26Xfmka2bOzGq228wxsD16ekY6Ao/wzKw5nPDMrDGc8MysMVKrpaxqe36SpDOrDcnMrB4e4ZlZYzjhmVljpJ6W0n6LxpnAZe0rSToZOBlg/vz5+dGZmVUodYS3JiIWDDyATwy2Umu1lN5e1xYws9HBU1ozawwnPDNrDCc8M2uM1GopU9uenwucW2E8Zma18QjPzBqj9mopuVY9reQ2U7fOrBixaVN6m8wqH9lWrkxvM21a9XGMFv396W0mjNqv+zP6lqV/7yG3IlHzeIRnZo3hhGdmjbFFCa+9mICZ2WjmEZ6ZNYYTnpk1hhOemTVGLQlP0smSlkha0tfXV0cXZmbJakl4rpZiZqORp7Rm1hhOeGbWGFuU8NqLCZiZjWYe4ZlZY4zaq6mnTtyQ0Srz44yB4gF/Xjs9uc2cab6gfKxxEYB6eYRnZo3hhGdmjeGEZ2aNkZXwXCXFzMYij/DMrDGc8MysMVw8wMwaw8UDzKwxPKU1s8ZwwjOzxshNeFMkPdTy+FClUZmZ1SDr4tOI8MjQzMYcJy4za4zaq6X090PfMiW36/3lZemdveIV6W2AdbN3TG4zkbyqFhs3pW8LgNe8Jr3N0iXrs/rKrgQzIePrlFOpBuCSS9LbLFqU11c3rVyZ1y73Zza1eyUtg/TvvjJ/z4biEZ6ZNYYTnpk1RnLCay0cIOmNku6QtHO1YZmZVS97H56kI4AvAa+PiPurC8nMrB5ZCU/SIcB/Am+MiLurDcnMrB45CW8i8APgsIj4Q8XxmJnVJuegxQbgeuA9Q63QWi1l+XJXSzGz0SEn4W0CjgdeLuljg63QWi1l1ixXSzGz0SH30rLVko4CrpH054j4RsVxmZlVLvsobUSskPQG4JeS+iIi49IIM7PuSU54ETG15f8PAn9RaURmZjXxlRZm1hi1Fw+YMAF6Z2dcAPzmN1cfzBByCwHkGD8ur6+lt+S06snqq6tyL3ofC4UAckybNtIR1KbqQgA5PMIzs8ZwwjOzxnDCM7PG6CjhSQpJ3255PkFSn6Qf1ReamVm1Oh3hPQ3sK2ly+fy1wMP1hGRmVo+UKe0VwFHl/08ALqg+HDOz+qQkvAuBRZImAfsDvxpqxdbiAX19Lh5gZqNDxwkvIm4BdqEY3V0xzLrPFA/o7XXxADMbHVJPPL4M+DxwGDCr8mjMzGqUmvC+CTwREUslHVZDPGZmtUlKeBHxEMV9LMzMxpyOEl5rhZSWZVcDV1ccj5lZbXylhZk1hhOemTWGE56ZNYYTnpk1hhOemTVGcsVjSRuBpS2Ljo2I+yqLyMysJjkl3tdExILKIzEzq5mntGbWGDkJb7Kkm8rH9wdbwdVSzGw0qmVKGxFnA2cDLFy4cORvVWRmhqe0ZtYgTnhm1hhOeGbWGMkJb7DKKWZmY4FHeGbWGE54ZtYYTnhm1hhOeGbWGFtSPGAC8HvgXRGxuurAzMyqljPCWxMRCyJiX2A98N6KYzIzq8WWTmmvAV5URSBmZnXLTniSJgBH8tzaeAOvuXiAmY062dVSgCXAA8A32leIiLMjYmFELOzt7d3SGM3MKuECoGbWGD4txcwawwnPzBrDxQPMrDE8wjOzxsg5aJFkwwb482NKbjen5/HkNuumzEhuAzBx3Ib0RhMyN93qzItSNm1KbrJ8/TZZXc3admNWO8Z18e9nf396m9yfWRc9uTL9dwVg+qR1eR329OS1y5HxHa76O+URnpk1hhOemTWGE56ZNUZywpP0cUm3SbqlvDftX9YRmJlZ1ZL24kp6JXA0cEBErJM0G+jiXk8zs3yph622B5ZFxDqAiFhWfUhmZvVIndL+BJgn6Q5JX5F06GArtVZLWbHC1VLMbHRISngRsQo4EDgZ6AMuknTSIOs9Uy1l5kxXSzGz0SH5TMyI2AhcDVwtaSnwLuDcasMyM6te0ghP0h6SXtyyaAFwf7UhmZnVI3WENxX4sqRtgX7gLorprZnZqJeU8CLiRuCgmmIxM6uVr7Qws8aovXzEVlvBnO0iud3jT6RXPpnRk94PAJu6mPenTOlaV7PI3B7+Ozhipk/L/ZmNgfP/u1lNZ6gQRjoAM7NuccIzs8ZIvZZ2FvCz8ulcYCPFCcgAL4+I9RXGZmZWqdSjtMspzr1D0unAqoj4fA1xmZlVzlNaM2sMJzwza4xaEl5rtZS+PldLMbPRoZaE11otpbfX1VLMbHTwlNbMGsMJz8waI/vSsog4vcI4zMxq5xGemTVG7cUDcs3YNvci6gyj4KJmSzRh1H51bRTzb7qZNYYTnpk1hhOemTXGsAlP0hclfbDl+Y8lfb3l+RmSPlRXgGZmVelkhHcd5X0sJI0DZgP7tLx+EHB99aGZmVWrk4R3PfDK8v/7ALcCT0maIWkisBfw25riMzOrzLDH9iPiEUn9kuZTjOYWAztSJMEngaXthT8lnUx5+8b58+dXHrSZWY5OD1pcT5HsBhLe4pbn17Wv7OIBZjYadZrwBvbj7Ucxpb2BYoTn/XdmNmakjPCOBlZExMaIWAFsS5H0nPDMbEzoNOEtpTg6e0PbsicjYlnlUZmZ1aCjCxIjYiMwrW3ZSXUEZGZWF19pYWaN4YRnZo3hhGdmjeGEZ2aN4YRnZo2RlPBUuFbSkS3L3ibpyupDMzOrVlKd7IgISe8F/kvSz8v2/wa8oY7gzMyqlHxjgIi4VdIPgY8AWwPfioi7K4/MzKxiuXdC+SRFSaj1wML2F10txcxGo6yDFhHxNHARcH5ErBvkdVdLMbNRZ0uO0m4qH2ZmY4JPSzGzxnDCM7PGyL59e0ScXmEcZma18wjPzBoje4SXZFP6sY3f3TI+uc2uuyY3AWB6z5r0RpMm5XXW35/XbllGndW5c/P6GgtWrkxvM23a8OuMtNzvx7jMsUtuuwzr1iu5zcSeqDQGj/DMrDGc8MysMTpOeJJ+Lun1bcs+KOmr1YdlZla9lBHeBcCitmWLyuVmZqNeSsK7BDhKUg+ApF2AHYBrqg/LzKx6HSe88l60vwYGauEtAi6OiOcdRpF0sqQlkpb09fVVE6mZ2RZKPWjROq0dcjrr4gFmNhqlJrxLgSMkHQBMiYgba4jJzKwWSQkvIlYBPwe+iQ9WmNkYk3Me3gXAS3DCM7MxJqfE+w+A9GtEzMxGmK+0MLPG6E7xgIwLlF+6oNqLhjcvsxBAjgmZm/yFXAggx1goBJAj9/sxBlRdCCCHR3hm1hhOeGbWGE54ZtYYTnhm1hhOeGbWGE54ZtYYtSQ8V0sxs9GoloTnailmNhplJzxJ75N0U/nYocqgzMzqsCU34j4LOKvCWMzMauWDFmbWGE54ZtYYTnhm1hga5B481XYg9QH3D/HybGBZ4lvmtOl2uxdqX7ntXqh95bZzX/W22zkiBj89JCJG7AEs6Uabbrd7ofY1FmL09mhGX7ntPKU1s8ZwwjOzxhjphHd2l9p0u90Lta/cdi/UvnLbua8Ralf7QQszs9FipEd4ZmZd44RnZo0xYglP0rGSQtKeHa6/sSxUcKukH0ratoM2IemMlucflnR6Ql+3SbpZ0mmSht1WkuZKulDS3ZJulHSFpN2HabOTpEsl3SnpHklnSpqYEOPA46PDtSnbzZH03bKvGyUtlnTcMG1WtT0/SdKZnfQ3WPuq129vI+mNku6QtPMwbULSt1ueT5DUJ+lHHfTXvv13SWhzq6T/kjRluDZlu4+X38VbyvZ/Ocz6s1rielTSwy3PewZZ/4uSPtjy/MeSvt7y/AxJHxqiL0m6VtKRLcveJunKYWL8uaTXty37oKSvbq7dFss5/6WKB3ARcA3wyQ7XX9Xy//OAj3fQZi1wLzC7fP5h4PTEvrYD/t9wcVLcnHwx8N6WZS8BDh6mza+Bd5fPxwPfAP5PSowJ23ywGHcG3p/SF3AScGZCv0mxZn62VeW/RwB3Abt10ga4CZhcPj+yfP6jOmMs//8d4EMdtHll+TObWD6fDeyQ0OfpwIeHWeetwMXl/8cBNwKLW15fDLxiM+33BX5Pcb/TqcCdw21/4GTgnLZlNwCHpG7XlMeIjPAkTQVeBbwHWJTxFouBHTtYr5/iSM6pGX0AEBGPUfxwTpGkzaz6amBDRHytpe3NEXHNZtocDqyNiHPK9TeWsZ5YbqOqHQ6sb4vx/oj4cg19dZ2kQ4D/BI6OiLs7bHYFcFT5/xOAC+qIbRDXAC/qYL3tgWURsQ4gIpZFxCMVx3I9RWIF2Ae4FXhK0oxytrEX8NuhGkfErcAPgY8AnwC+1cH2vwQ4amDEWY6Qd6DYLrUZqSntXwFXRsQdwHJJB3baUNJ4ir/il3XY5CzgnZKmp4dZiIh7KEZf221mtX0p/jKm2Ke9TUSsBO5j+F+GyW1Tqrd32N+QX9xO+wI+lfEedZsI/AA4NiL+kNDuQmCRpEnA/sCvOmzXuk2+nxKopAkUo8mlHaz+E2BeOUX/iqRDU/rqRJlA+yXNBw6iGFD8iiIJLgSWRsT6Yd7mk8A7KD7X5zrocwXF7GZgKryIYpRZ62kjI5XwTqD4olH+e0IHbSaXv2yPAnOAn3bSUZlAvgV8ICPO0WxNRCxoeVyU+gaSzir3Uf4mpS+Kv+KjzQaKkcp7UhpFxC3ALhTfwSsSmrZuk83uA20x8B1eAjxAsftiuPhWAQdSzDL6gIsknZQQZ6eup0h2Awlvccvz6zqI82mK3VTnD4xGO3ABz87wFtHh6FpbUHy46wlP0kyKqdXXJd0H/ANw/DDTRSi/YBT7nAS8L6Hb/6D4Rdg6PWKQtCuwEXhsM6vdRvHFTHF7extJ04C5wB8T36sTtwEHDDyJiPdRjJZfCHX4NwHHAy+X9LHEtpcBn6f+6Wxrknx/B6MmoNjVERFXR8S/AKcAb6khtusoktt+FFPaGyhGeAdRJMNObCofnboUOELSAcCUiOhohhQRZ7Vsx6Tp/UiM8N5K8Vdg54jYJSLmURxYOLiTxhGxmmK0dlo5NeikzQrgYhL/+gNI6gW+RrGTfnPD7auAiZJObmm7v6TNfa6fAVMknViuPx44o+xrTWqsHbgKmCTp71uWdXSkcCwovxtHUezCSPlZf5PioFQnU8yukrSHpBe3LFrA0NWHtsT1wNHAijLBrgC2pUh6nSa8JOXo9ecU278r+05HIuGdALTv8/genU1rAYiI3wG3pLShSCSzO1x3YP/MbRRHaH9CsY9iczEFcBzwGhWnpdwGfIZiCj5cm7dKuhNYDmyKiE8nxDjw+OxwDcr+jgUOlXSvpF9THPH+SAf9ddMUSQ+1PAY9JWIw5S/qG4B/lvSmDts8FBFfyg22ZlOB8yTdLukWYG+KI69VW0rx+3FD27InIyKndFOnLqA4m6ErCc+Xlo0ikg6i+MEfFxE5BxfMbDOc8MysMXxpmSVCvvMAAAAiSURBVJk1hhOemTWGE56ZNYYTnpk1hhOemTWGE56ZNcb/B+gXRLcj44XkAAAAAElFTkSuQmCC\n",
            "text/plain": [
              "<Figure size 378x378 with 1 Axes>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        },
        {
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAATwAAAFMCAYAAABWJwDHAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deZxddZnn8c83e0IBMSQQ9qggyqJpiLaiKEL3KIsKo62gDmI7k7FHdFTsduux0dbW9iVtj4LyQlQQlUVFFk0rjoIghCUoEiJKZJNFoJIYQshaqWf+OKeam2st9/erc26qcr7v16teyT33PPf31KlbT/3O9lxFBGZmTTBhWydgZtYtLnhm1hgueGbWGC54ZtYYLnhm1hgueGbWGC5444SkEyU9KGmtpL/Y1vkMRtJHJJ03zPNvkXR1N3MaqySdKukX2zqPpnHB6wJJH5b0H23Llg+x7KQhXuZzwGkR0RMRv6or19GIiH+JiP8OIGmepJA0qeX5b0XEf+lmTpJOk7RE0kZJ5w/y/AxJX5K0QtITkq7rZn5DKbfdU+UfuLWtf0gkTZV0jqTHJK2SdJWkPbdlvuOFC153XAccLmkigKTdgcnAX7Qt269cdzD7AstyBh8Yo6EeAT4JfG2I588FZgHPK/99X5fy6sQLyj9wPQN/SEr/G3gJ8HxgD+BPwBe3RYLjjQted9xKUeDml4+PAK4Bfte27J6IeKQ1sPxrvhaYCPxa0j3l8udJulbSaknLJL22JeZ8SV+WtEjSU8Ar2xMqYz8t6RZJayRdIWlWy/OvLV93dbnu81qe+6CkhyU9Kel3ko4ul58h6ZvlagOFe3U5Q3lJ+26cpMMl3VrOrG6VdHhbfv8s6YZynKslzS6fmybpm5JWlvndKmm3wTZ8RFwWEZcDKwfZBs8FXgssjIjeiNgSEbcN9jrl+m+XdFeZz72S/mfLc0dKekjS6ZIel/RHSW9veX4XSVeW2/oW4NlDjdOBZwI/jojHImIDcAlw0CherzFc8LogIjYBNwMvLxe9HLge+EXbsj+b3UXExojoKR++ICKeLWkycBVwNbAr8G7gW5IOaAl9M/ApYMdynMGcAvwtsDvQB3wBQNJzgIuA9wJzgEXAVZKmlGOcBrwwInYEXgXcP8hrD3xfM8sZyuLWJ8vi+sNyzF2AfwN+KGmXtu/h7eX3OAX4QLn8bcDOwN5l7DuB9UN8j8N5EfAA8PFyl3appNcPs/7jwPHATmVen5d0aMvzc8u89gTeAZwt6Rnlc2cDGyi29d+WXyO5TtKjki6TNK9l+VeBl0raQ9IM4C3Afwz2ArY1F7zu+TlPF4EjKAre9W3Lft7ha70Y6AE+ExGbIuJnwA+Ak1vWuSIiboiI/nIWMJgLI+LOiHgK+D/AG8vd3zcBP4yIn0TEZorjh9OBw4EtwFTgQEmTI+L+iLinw7xbHQcsj4gLI6IvIi4Cfgu8pmWdr0fE3RGxHriUp2fDmykK3X4Ds7KIWJORw17AwcATFLuGpwEXtM5mW0XEDyPinij8nOIPzhEtq2wGPhERmyNiEbAWOKDcpq8HPhYRT0XEncAFI+T2CmAe8FyK3fIftBwPXQ48CDwMrKHYHf9E2rfeTC543XMd8LJyZjMnIpYDN1Ic25tF8YvX6QHzPYAHI6K/ZdkDFDOLAQ928Dqt6zxAsds9u3z9BwaeKMd5ENgzIn5PMfM7A3hc0sWS9ugw7/bv4YG2Ze3fw6Mt/19HUeQBLgR+DFws6RFJny1nvanWUxSpT5Z/OH5Ocahh0BMrko6RdFN5omA1cCzF9hqwMiL6Bsl5DjCJP9/eQ4qI68qcVlMcs3smRWGDYrY4laLo7wBchmd4HXHB657FFLs7/wO4AaCclTxSLnskIu7r8LUeAfaW1Prz24fiL/6ATtrg7N0WvxlYUb7+vgNPSFK57sNl3t+OiJeV6wTwr4O89kjjbzVGSw4PD7Lu1i9czKA+HhEHUsw6j6fYPU91x2AvP9iKkqYC36OY7e4WETMpdvXVwTi9FIcM2rd3imgZaz5wfkSsioiNFCcsXjRwjNOG5oLXJeVu2RLg/RS7sgN+US5LuRziZorZwz9ImizpSIpdwYsT03qrpAPL40CfAL4bEVsodh+Pk3R0OXM6HdgI3CjpAElHlQVgA8UsqX+Q1+4tlz9riLEXAc+R9GZJkyS9CTiQYtd8WJJeKemQcldxDUWhHiwHyteeRnHSZ2J5wmNg1/A64A/Ah8v1XkpxgufHg7zUFIpZVS/QJ+kYhpgJtiu36WXAGSougzmQ4jjkUN/fQZLmS5ooqQc4k+IPwV3lKrcCp0jaufz5/C+KP5grOsmnyVzwuuvnFAfgW08iXF8u67jglSdBXgMcQzEj+xJwSkT8NjGfC4HzKXYdpwHvKV//d8BbKWYOK8qxXlOOOxX4TLn80TL3Dw+S4zqKkyY3lGdSX9z2/EqKmdnpFGdQ/wE4vsNf2rnAdymK3V0U2/XCIdb9R4qi/KHye1pfLqM8Pvk6il3TJ4CvMMR2jIgnKbbPpRSXgbwZuLKDXAecRrF7+yjFNv/6MOvuRnHmdQ1wL8WxvOPLfKE4ebOB4lheb5n/iQm5NJbcALSZJF0LfDMihrwzwmx74xmemTWGC56ZNYZ3ac2sMTzDM7PGmDTyKqMza9bs2HPPeclxEzJKcU4MwMR1T6YH7bhj3mB9fSOvM5j+Qa+6GN6k2n+8W1Mnl6RVEEPe5pgw+JUrHQR2b16wPucGOWDatLy4zM2fJ+OHtqkvfds/9ND9rFq1YtDvrPbfiD33nMdlly1Jjttpp/Sxcn/oO//ymvSgI4/MG2xF5qVS69alx8zOvA419xc8J27KlKyhnlyb/pu646QuV5MMdy7Lq0D775833tQpXTyktWGoOxyH9tDK6ckxxx67YMjnvEtrZo3hgmdmjZFc8CSdoKIb63PrSMjMrC45M7yTKW6NOnmkFc3MxpKkglfeyPwyiuaGQ332gpnZmJQ6w3sd8KOIuBtYKemwwVaStFDFB6csWbWqd9RJmplVIbXgnczTLYguZojd2og4NyIWRMSCWbPmjCY/M7PKdHwdXtmV9yjgEElB0V8sJP19+P40MxsHUmZ4b6D4DIR9I2JeROwN3MfWPf3NzMaslIJ3MvD9tmXfw2drzWyc6HiXNiL+7LNNI+IL1aZjZlYf32lhZo1Re/OAqetXs/9vrkgPzLg5P3baOX0cQEf92eR15LH6M8/TZN4s/6uH0s927zUjayjmTMvoHgNZN4dndYkAdrziO8kxm9/01qyxJnf0AXCDePTRkddpc/C+PSOvNIgtkzK793TT2rXJIXttTt+GU2LjkM95hmdmjeGCZ2aN4YJnZo2RdAxP0hZgaRl3H/DfImJ1HYmZmVUtdYa3PiLmR8TBwCrgXTXkZGZWi9Hs0i4G9qwqETOzumUVPEkTgaOBK4d4/j+7pfSuWTOa/MzMKpNa8KZLuh14FNgN+MlgK7V2S5mTeZ2VmVnVso7hAfsCwsfwzGwcydqljYh1wHuA0yV1+cNPzczyZJ+0iIhfAXfgbilmNk4kzc4ioqft8WuqTcfMrD71747OnAnHH58cFhMmJsco8ybv2LApIyqvCUDuzfI5nyzfs0NuI+q8G9jpyYzLsP716Y0Apk/qcmPuuXO7NtTE3AYH3TR7dnrMzJnpMcM06PCtZWbWGC54ZtYYLnhm1hipH8Qdks5sefwBSWdUnpWZWQ1SZ3gbgf8qKePoo5nZtpVa8PqAc4H31ZCLmVmtco7hnQ28RdKQHyCxVfOA3t787MzMKpRc8CJiDfANilvLhlrn6eYBc9I/fMbMrA65Z2n/HXgHsEOFuZiZ1Sq3ecAq4FKKomdmNi6M5jq8MwGfrTWzcSO7eUBEPAZkftSzmVn3+U4LM2uM+rul9PfDpvRuJJo2LTnmybVKjgHYcUJfetAwHRmG88SavBz7+7PCtluZm7+7cn5ouT/oSdtpH96+jN/NGLpzjGd4ZtYYLnhm1hgueGbWGB3v+EvaAiwFJlPcU/sN4PMR4aNLZjYupBzpHPiIRiTtCnwb2An4pzoSMzOrWu6dFo8DC4HTJOWddjQz67LRfEzjvcBEYNf257bqlrJixWjyMzOrTC0nLbbqlpLzSUVmZjXILniSngVsAR6vLh0zs/pkFTxJc4BzgLMihrms2cxsDEk5Sztd0u08fVnKhcC/1ZKVmVkNOi54ETGxzkTMzOpW/x3HEyZARiOAHDv25O5dd6/L1cTMPxs792zJiNp+b6SZOGEcHEmZkLH9c2K2Zzm1Y5ht6K1rZo3hgmdmjeGCZ2aNkVzwJM2VdLGkeyTdJmmRpOfUkZyZWZWSTlqU981+H7ggIk4ql70A2A24u/r0zMyqk3qW9pXA5og4Z2BBRPy62pTMzOqRukt7MHBbHYmYmdWtlpMWW3VL6e2tYwgzs2SpBW8ZcNhIK23VLWXOnLzMzMwqllrwfgZMlbRwYIGk50s6otq0zMyql1Twys4oJwJ/VV6Wsgz4NPBoHcmZmVUp+V7aiHgEeGMNuZiZ1cp3WphZY9TeLWXjRrjn3vTP+Xn2PptryGYIN9+cHvPSl2YNtXp1Vhg9Kx9KD5o7N2us9f1Ts+L+8If0mMwU6elJj5nYn/eeWt83OStuOuuTYzZPmp411uQJOd106G53lr6+7sT0D/3JsZ7hmVljuOCZWWO44JlZYyQVPEl7SbpC0nJJ90o6S1LeAR8zsy7ruOCVnVIuAy6PiP2B/YHpwGdrys3MrFIpM7yjgA0R8XWAiNgCvA84RVLGOTMzs+5KKXgH0dYpJSLWAPcD+7Uub20esGqVmweY2dhQy0mL1uYBs2a5eYCZjQ0pBe83tHVKkbQTMBf4XZVJmZnVIaXg/RSYIekUAEkTgTOBsyIi/ZJyM7Mu67jgtXRKeYOk5cBKoD8iPlVXcmZmVUptD/VgRLy2vCzlWODVkg6tJzUzs2plNw+IiBuBfSvMxcysVrV3S8m1ck16h4pdJj2RN9jy5ekxmd1S9lqfMRbAeeclh2z+5L/mjTV0s4lhHTBvY3rQpLy34LGvmZgcs+iHeWNNz+xEEhPSO59MXpP3Hl47ceesuJ4dIisuxwOPpP9OT5uWHrN5y9A7rr6X1swawwXPzBrDBc/MGiP5oIakLcDSlkUXR8RnqkvJzKweOUdx10fE/MozMTOrmXdpzawxcgredEm3t3y9qX0Fd0sxs7Goll3aiDgXOBfgkEMWdO9CHzOzYXiX1swawwXPzBojZ5d2uqTbWx7/KCI+VFVCZmZ1SS54EZF+E6OZ2Rigos1dfebPXxA/+cmS5Lg5s7t4rmPTpvSYKVPyxurPvDN/w4b0mBkz8sYyGyM29yk55sUvXsBtty0ZNNDH8MysMVzwzKwxXPDMrDGSCp6k3SR9W9K9km6TtFjSiXUlZ2ZWpY4LniQBlwPXRcSzIuIw4CRgr7qSMzOrUsplKUcBmyLinIEFEfEA8MXKszIzq0HKLu1BwC87WbG1ecDKlW4eYGZjQ/ZJC0lnS/q1pFvbn4uIcyNiQUQs2GWXOaPL0MysIikFbxnwn59BGxHvAo4GXNHMbFxIKXg/A6ZJ+ruWZb6U38zGjY4LXhT3oJ0AvELSfZJuAS4APlhXcmZmVUpqHhARf6S4FMXMbNzxnRZm1hg5/fCSTRjrZXXduvSY3G4pmTZO2iE5ZirbcXf9bna46abcbjpj/pcsz8aN6THDNYDaPreSmdkgXPDMrDFc8MysMVK7paxte3yqpLOqTcnMrB6e4ZlZY7jgmVljpF6W0v4RjbOAK9tXkrQQWAiw11775GdnZlah1Bne+oiYP/AFfGywldwtxczGIu/SmlljuOCZWWO44JlZY6R2S+lpe3w+cH6F+ZiZ1cYzPDNrjNq7pUyYADPGel/kmTO3dQYjmtr3VHrQlLwNv7lPWXEbNqTH7NiT19Fl84SpyTGTx0P3mO2060munonrk2MmMHTHGW9dM2sMFzwza4xRFbz2ZgJmZmOZZ3hm1hgueGbWGC54ZtYYtRQ8SQslLZG0ZMWK3jqGMDNLVkvBa+2WMnu2u6WY2djgXVozawwXPDNrjFEVvPZmAmZmY5lneGbWGF1pHjB9ypacyMpzGcqW/vSb5SdOGAc3omeaPCnve5vcxfn+5Alj+z0F8MdH099Xu8/dft9XWaZMSY8ZpgGDZ3hm1hgueGbWGC54ZtYYWQXPXVLMbDzyDM/MGsMFz8wao/bmAb29bh5gZmND7c0D5sxx8wAzGxu8S2tmjeGCZ2aNkVvwZkh6qOXr/ZVmZWZWg6x7aSPCM0MzG3dcuMysMWrvlsL69XDXXelxz3xmesykvG/nykVTk2NOPCFrKFixIivsr9+6W3LMokVZQzF5w5N5gcN0qRjStGl5Y113XXrM4YfnjZXTsQPYfcYT6UH3r8oai7lz8+Jyt3+OdevSYzZsSI/p6xvyKc/wzKwxXPDMrDGSC15r4wBJx0q6W9K+1aZlZla97GN4ko4GvgC8KiIeqC4lM7N6ZBU8SS8HvgIcGxH3VJuSmVk9cgreVOBy4MiI+G3F+ZiZ1SbnpMVm4EbgHUOtsFW3lD/9KTs5M7Mq5RS8fuCNwIskfWSwFbbqlvKMZ4wqQTOzquTeWrZO0nHA9ZIei4ivVpyXmVnlss/SRsQqSa8GrpPUGxFXVpiXmVnlkgteRPS0/P9BIOMeMDOz7vOdFmbWGIqIWgdYsGBBLLn11lrHMDMbsOCFL2TJkiUa7DnP8MysMVzwzKwxXPDMrDE6KniSQtI3Wx5PktQr6Qf1pWZmVq1OZ3hPAQdLml4+/mvg4XpSMjOrR8ou7SLguPL/JwMXVZ+OmVl9UgrexcBJkqYBzwduHmrFrZoH9PaONkczs0p0XPAi4g5gHsXsbtiPh9mqecCcOaPL0MysIqm3ll0JfA44Etil8mzMzGqUWvC+BqyOiKWSjqwhHzOz2iQVvIh4iOJzLMzMxp2OCl5rh5SWZdcC11acj5lZbXynhZk1hguemTWGC56ZNYYLnpk1hguemTVG8mdaSNoCLG1ZdEJE3F9ZRmZmNcn51LL1ETG/8kzMzGrmXVoza4ycgjdd0u3l1/cHW8HdUsxsLKpllzYizgXOheJTy3ISMzOrmndpzawxXPDMrDFc8MysMZIL3mCdU8zMxgPP8MysMXLO0qbr70+PmdDFWrxhQ3rMtGnV5zGMJ9YoOWbnnbbjE+Rj/T2V6cm16T9ngB17ttOf9aZN6TEx9LYY++8AM7OKuOCZWWOMpnnAJOAu4G0Rsa7qxMzMqpYzw1sfEfMj4mBgE/DOinMyM6vFaHdprwf2qyIRM7O6ZRc8SZOAY9i6N97Ac24eYGZjTna3FGAJ8Afgq+0rRMS5EbEgIhbMmTNntDmamVXCDUDNrDF8WYqZNYYLnpk1hpsHmFljeIZnZo1Re/OACNjYNzE5buqU9Juh1z6Vd+N1z+Tu1f0/PpqX4+47PZUe1N/dBgddvTl/XcbNPT1d3jnJaHCw47SMpghAt/qAjMZjj6e/93ebltHYY5jt7hmemTWGC56ZNYYLnpk1RnLBk/RRScsk3VF+Nu1f1pGYmVnVko50SnoJcDxwaERslDQbmFJLZmZmFUs9tbM7sCIiNgJExIrqUzIzq0fqLu3VwN6S7pb0JUmvGGwld0sxs7EoqeBFxFrgMGAh0AtcIunUQdZztxQzG3OSr1aMiC3AtcC1kpYCbwPOrzYtM7PqJc3wJB0gaf+WRfOBB6pNycysHqkzvB7gi5JmAn3A7yl2b83MxrykghcRtwGH15SLmVmtfKeFmTVG7S0W1LeZqasfSw+cPTs5pGd6+jAA3PGb9Jj5eV3uc5p8ALAp45LHPfbIGqp39eSsuDmztqQH5XZYWbMmPWbGjLyxMnNcfm96l6CenvQYgN13zdj20NUON7vNTs/xnvt3To4ZrjuTZ3hm1hgueGbWGKn30u4C/LR8OBfYQnEBMsCLImJThbmZmVUq9SztSopr75B0BrA2Ij5XQ15mZpXzLq2ZNYYLnpk1Ri0Fb6tuKatW1TGEmVmyWgreVt1SZs2qYwgzs2TepTWzxnDBM7PGyL61LCLOqDAPM7PaeYZnZo1Re/MAJk+GXXetfZhRyWwEkOPZz4rMyH0qzWM4c2bn5tjFv5+ZjRG6af/9crdjjnEwd8loVJDz+zJ16jApJL+amdk45YJnZo3hgmdmjTFiwZP0eUnvbXn8Y0nntTw+U9L760rQzKwqnczwbqD8HAtJE4DZwEEtzx8O3Fh9amZm1eqk4N0IvKT8/0HAncCTkp4haSrwPOCXNeVnZlaZES9LiYhHJPVJ2odiNrcY2JOiCD4BLG1v/ClpIeXHN+6zT/cupzAzG06nJy1upCh2AwVvccvjG9pX3qp5wJw5VeVqZjYqnRa8geN4h1Ds0t5EMcPz8TszGzdSZnjHA6siYktErAJmUhQ9FzwzGxc6LXhLKc7O3tS27ImIyPjAVDOz7uvoXtqI2ALs1Lbs1DoSMjOri++0MLPGcMEzs8ZwwTOzxnDBM7PGcMEzs8ZIKngq/ELSMS3L/kbSj6pPzcysWkkt3iMiJL0T+I6ka8r4fwFeXUdyZmZVSv5Mi4i4U9JVwAeBHYBvRMQ9lWdmZlax3A/x+ThFS6hNwIL2J90txczGoqyTFhHxFHAJcGFEbBzkeXdLMbMxZzRnafvLLzOzccGXpZhZY7jgmVlj5J60ICLOqDAPM7PaeYZnZo2RPcNL0p9xbmNCF2vxhg3pMdOm5Y2Vsy0gL8cZM/LGGg8ytuMfH5+YNdTucyMrLkvu+6Obvy9dtPYpJccMtwm3z61kZjYIFzwza4yOC56kayS9qm3ZeyV9ufq0zMyqlzLDuwg4qW3ZSeVyM7MxL6XgfRc4TtIUAEnzgD2A66tPy8yseh0XvPKzaG8BBnrhnQRcGhF/dgpL0kJJSyQt6e3trSZTM7NRSj1p0bpbO+TurJsHmNlYlFrwrgCOlnQoMCMibqshJzOzWiQVvIhYC1wDfA2frDCzcSbnOryLgBfggmdm40xOi/fLgfT7PczMtjHfaWFmjdGd5gFj/cbm3EYAOXK3xfbcCCBHxnbsahOAXGP9d6XLenZI/5kNtwm9dc2sMVzwzKwxXPDMrDFc8MysMVzwzKwxXPDMrDFqKXjulmJmY1EtBc/dUsxsLMoueJLeJen28muPKpMyM6vDaD6I+2zg7ApzMTOrlU9amFljuOCZWWO44JlZY2iQz+CpdgCpF3hgiKdnAysSXzInpttx2+tYuXHb61i5cR6r3rh9I2Lwy0MiYpt9AUu6EdPtuO11rPGQo7dHM8bKjfMurZk1hguemTXGti5453Yppttx2+tYuXHb61i5cR5rG8XVftLCzGys2NYzPDOzrnHBM7PG2GYFT9IJkkLScztcf0vZqOBOSVdJmtlBTEg6s+XxBySdkTDWMkm/lnS6pBG3laS5ki6WdI+k2yQtkvScEWL2knSFpOWS7pV0lqSpCTkOfH1opJgybjdJ3y7Huk3SYkknjhCztu3xqZLO6mS8weKrXr89RtKxku6WtO8IMSHpmy2PJ0nqlfSDDsZr3/7zEmLulPQdSR19FJ2kj5bvxTvK+L8cYf1dWvJ6VNLDLY+nDLL+5yW9t+XxjyWd1/L4TEnvH2IsSfqFpGNalv2NpB+NkOM1kl7Vtuy9kr48XNyo5Vz/UsUXcAlwPfDxDtdf2/L/C4CPdhCzAbgPmF0+/gBwRuJYuwL/b6Q8KT6cfDHwzpZlLwCOGCHmFuDt5eOJwFeB/5uSY8I2HyzHfYF3p4wFnAqclTBuUq6Z39va8t+jgd8Dz+4kBrgdmF4+PqZ8/IM6cyz//y3g/R3EvKT8mU0tH88G9kgY8wzgAyOs8wbg0vL/E4DbgMUtzy8GXjxM/MHAXcA0oAdYPtL2BxYCX29bdhPw8tTtmvK1TWZ4knqAlwHvAE7KeInFwJ4drNdHcSbnfRljABARj1P8cE6TpGFWfSWwOSLOaYn9dURcP0zMUcCGiPh6uf6WMtdTym1UtaOATW05PhARX6xhrK6T9HLgK8DxEXFPh2GLgOPK/58MXFRHboO4Htivg/V2B1ZExEaAiFgREY9UnMuNFIUV4CDgTuBJSc8o9zaeB/xyqOCIuBO4Cvgg8DHgGx1s/+8Cxw3MOMsZ8h4U26U222qX9nXAjyLibmClpMM6DZQ0keKv+JUdhpwNvEXSzulpFiLiXorZ167DrHYwxV/GFAe1x0TEGuB+Rv5lmN62S/WmDscb8o3b6VjAJzJeo25TgcuBEyLitwlxFwMnSZoGPB+4ucO41m3y/ZREJU2imE0u7WD1q4G9y130L0l6RcpYnSgLaJ+kfYDDKSYUN1MUwQXA0ojYNMLLfBx4M8X39dkOxlxFsXczsCt8EsUss9bLRrZVwTuZ4o1G+e/JHcRML3/ZHgV2A37SyUBlAfkG8J6MPMey9RExv+XrktQXkHR2eYzy1pSxKP6KjzWbKWYq70gJiog7gHkU78FFCaGt22TYY6AtBt7DS4A/UBy+GCm/tcBhFHsZvcAlkk5NyLNTN1IUu4GCt7jl8Q0d5PkUxWGqCwdmox24iKf38E6iw9m1RtF8uOsFT9Isil2r8yTdD/w98MYRdhehfINRHHMS8K6EYf+d4hdhh/SMQdKzgC3A48OstozijZniN+0xknYC5gK/S3ytTiwDDh14EBHvopgtbw99+PuBNwIvkvSRxNgrgc9R/+5sa5F8dwezJqA41BER10bEPwGnAa+vIbcbKIrbIRS7tDdRzPAOpyiGnegvvzp1BXC0pEOBGRHR0R5SRJzdsh2Tdu+3xQzvDRR/BfaNiHkRsTfFiYUjOgmOiHUUs7XTy12DTmJWAZeS+NcfQNIc4ByKg/TDTbd/BkyVtLAl9vmShvu+fgrMkHRKuf5E4MxyrPWpuXbgZ8A0SX/XsqyjM4XjQfneOI7iEEbKz/prFCelOtnF7CpJB0jav2XRfIbuPjQaNwLHA6vKArsKmElR9DoteEnK2es1FNu/K8dOt0XBOxloP+bxPTrbrQUgIn4F3JESQ1FIZne47sDxmWUUZ2ivpjhGMVxOASuYeiEAAACySURBVJwI/JWKy1KWAZ+m2AUfKeYNkpYDK4H+iPhUQo4DX58ZKaAc7wTgFZLuk3QLxRnvD3YwXjfNkPRQy9egl0QMpvxFfTXwj5Je22HMQxHxhdxka9YDXCDpN5LuAA6kOPNataUUvx83tS17IiJyWjd16iKKqxm6UvB8a9kYIulwih/8iRGRc3LBzIbhgmdmjeFby8ysMVzwzKwxXPDMrDFc8MysMVzwzKwxXPDMrDH+P11KDADyDyZSAAAAAElFTkSuQmCC\n",
            "text/plain": [
              "<Figure size 378x378 with 1 Axes>"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "xv2anXSUX9kr"
      },
      "source": [
        "## validating output\n",
        "(comparing to known output from GREMLIN_cpp version)"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "RBIvsIQ9X74s"
      },
      "source": [
        "!wget -q -nc http://files.ipd.uw.edu/krypton/4FAZA.out"
      ],
      "execution_count": 23,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "M7fK03GFX_3q",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 265
        },
        "outputId": "60366447-c199-4396-a590-d73fd6c20fb6"
      },
      "source": [
        "pd_test = pd.read_table(\"4FAZA.out\",sep=\" \")\n",
        "plt.scatter(pd_test[\"apc\"],mtx[\"apc\"])\n",
        "plt.show()"
      ],
      "execution_count": 24,
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXwAAAD4CAYAAADvsV2wAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAYNElEQVR4nO3df4xV5Z3H8feXcdTB1YIFbRmYQltkq6LYTkBjsrqtFNaN4qoVVLKaGEm6sdldumxhYWW0ddWSWpusmy11TVu1K8Xa22mgpfgrTYxYxgw/MmwRSi1ysUpbxzZ1qgN89497h14u9879cc699/z4vJIJ58dzOc/JwGeeeZ7nPMfcHRERSb4xra6AiIg0hwJfRCQlFPgiIimhwBcRSQkFvohISpzU6gqUM2HCBJ86dWqrqyEiEisvv/zyb9x9YqlzkQ38qVOn0tfX1+pqiIjEipn9qtw5demIiKSEAl9EJCUU+CIiKaHAFxFJCQW+iEhKRHaWjohI2mT6s6zZtJuDg0NMGtfBsnkzuOaiztD+fgW+iEgEZPqzrHhqJ0PDRwDIDg6x4qmdAKGFvrp0REQiYM2m3cfCfsTQ8BHWbNod2jUU+CIiEXBwcKim4/VQ4IuIRMCkcR01Ha+HAl9EJAKWzZtBR3vbccc62ttYNm9GaNfQoK2ISJOMNgtn5E/N0hERiblqZuFcc1FnqAFfTF06IiJN0IxZOJUo8EVEmqAZs3AqUeCLiDRBM2bhVKLAFxEJQaY/y6X3Pcu05Ru49L5nyfRnjzvfjFk4lWjQVkQkoGoHZKGxs3AqUeCLiAQ02oBsYaA3ehZOJaF06ZjZfDPbbWZ7zWx5ifNdZvacmfWb2Q4zuzKM64qIREEUBmSrETjwzawNeAj4G+Bc4EYzO7eo2Crgu+5+EbAI+K+g1xURiYooDMhWI4wW/mxgr7vvc/f3gCeABUVlHDgjv/0+4GAI1xURiYQoDMhWI4w+/E7gtYL9A8CcojI9wE/M7HPAacAVpf4iM1sCLAHo6uoKoWoiIsFVejFJFAZkq9GsQdsbgW+6+1fM7BLgUTM7392PFhZy97XAWoDu7m5vUt1ERMqq9sUkrR6QrUYYXTpZYErB/uT8sUK3Ad8FcPcXgVOBCSFcW0SkoaKwJEJYwgj8rcB0M5tmZieTG5TtLSqzH/gUgJl9jFzgHwrh2iIiDRWXGTjVCBz47n4YuAPYBPwfudk4A2Z2t5ldnS/2eeB2M9sO/C9wq7ury0ZEIi8uM3CqEUofvrtvBDYWHbuzYHsXcGkY1xIRaaZl82Yc14cP0ZyBUw09aSsiMoq4zMCphgJfRKSCOMzAqYZWyxQRSQkFvohISijwRURSQoEvIpISCnwRkZRQ4IuIpISmZYpIYlRa1TLtFPgikgiZ/izL1m9n+Ghu1Zbs4BDL1m8HUOjnqUtHRBKhp3fgWNiPGD7q9PQOtKhG0aPAF5FEGBwarul4GinwRURSQoEvIokwfmx7TcfTSIEvIomw+qrzaG+z4461txmrrzqvRTWKHs3SEZHIq2a6ZZKWMW4UBb6IRFamP0tP78BxA6+jTbdMyjLGjaIuHRGJpFWZnfzzum0lZ9loumV9FPgiEjmZ/iyPb9nPaC++1nTL2inwRSRy1mzaPWrYS33Uhy8ikVA4MFtN2Gu6Ze0U+CLScpn+LCue2snQ8JGqymu6ZX0U+CLScms27a467Ds13bJuCnwRabmDg0NlzxloTn1IFPgi0nKTxnWQLRH6neM6eGH5J1tQo2TSLB0Rabll82bQ0d523LGO9jaWzZvRoholk1r4ItJwlZZG0LIIzaHAF5GGKp6Bkx0cYsVTOwFOCH0FfGOFEvhmNh/4GtAGPOzu95UocwPQAziw3d1vCuPaIhI9hS36MWYc8eNn1g8NH2HNpt0K+CYLHPhm1gY8BMwFDgBbzazX3XcVlJkOrAAudfe3zOysoNcVkWgqbtEXh/2I0WbmSGOEMWg7G9jr7vvc/T3gCWBBUZnbgYfc/S0Ad38zhOuKSARVO6d+0riOJtRGCoXRpdMJvFawfwCYU1TmHAAze4Fct0+Pu/+4+C8ysyXAEoCurq4QqiYijVY8IFtqemUxzcBpjWYN2p4ETAcuByYDPzWzme4+WFjI3dcCawG6u7u1dpJIxJUakDUouRZOmxlH3TUDp4XCCPwsMKVgf3L+WKEDwEvuPgz80sxeIfcDYGsI1xeRFinVfeNwQuh3tLdx77UzFfItFkYf/lZguplNM7OTgUVAb1GZDLnWPWY2gVwXz74Qri0iLVRu4NXJPSVr+T8V9tEQuIXv7ofN7A5gE7n++UfcfcDM7gb63L03f+7TZrYLOAIsc/ffBr22iLSWlkSIl1D68N19I7Cx6NidBdsOLM1/iUiMjPaU7LJ5M05Y1lgDstGlJ21FpKxKT8lqSYR4UeCLSFmlBmWLn5LVkgjxodUyRaSscoOyeko2nhT4IlJWuadh9ZRsPCnwRaQsrVOfLOrDF5GyNCibLAp8ERmVBmWTQ106IiIpocAXEUkJdemIpECld8pKOijwRRIu059l2frtDB/NrV+ZHRxi2frtAAr9lFGXjkjC9fQOHAv7EcNHnZ7egRbVSFpFgS+ScINDwzUdl+RS4IuIpIQCXyThxo9tr+m4JJcCXyThVl91Hu1tdtyx9jZj9VXntahG0iqapSMSU9VOtdTyCDJCgS8SQ5VeTFJMyyMIqEtHJJZGezGJSDkKfJEY0otJpB7q0hGJuFJ99ZPGdZAtEe56MYmMRoEvElG5fvodDA0fPXZspK/+uk908r2Xs8d16+jFJFKJunREIijTn2Xpum3Hhf2IoeEjPPfzQ9x77Uw6x3VgQOe4Du69dqYGZmVUauGLRFBP7wAnRv2fHRwc0swbqZkCXyRCRvrrK61zo756qYcCXyQCMv1Z7vrhAG+9U3lBMwP11UtdFPgiLbYqs5PHt+zHKxcF4OaLu9SVI3UJZdDWzOab2W4z22tmy0cpd52ZuZl1h3FdkTjL9GeZdddPeKyGsF98cRdfumZmQ+slyRW4hW9mbcBDwFzgALDVzHrdfVdRudOBfwReCnpNkbgrXhqhkk6tfyMhCKNLZzaw1933AZjZE8ACYFdRuS8C9wPLQrimSCyNDMqWemiqlI72Nk23lNCEEfidwGsF+weAOYUFzOzjwBR332BmCnxJlVpDfsS4jnZ6rj5PYS+hafigrZmNAR4Abq2i7BJgCUBXV1djKybSBLV23UBuFs7N6quXBghj0DYLTCnYn5w/NuJ04HzgeTN7FbgY6C01cOvua9292927J06cGELVRFqr1KqWoxk/tp2vLpylsJeGCKOFvxWYbmbTyAX9IuCmkZPu/jYwYWTfzJ4H/sXd+0K4tkikVbt6pQZlpRkCB767HzazO4BNQBvwiLsPmNndQJ+79wa9hkhclVvVstCDC2cp6KUpzL3aGcDN1d3d7X19+iVA4mFVZiffeWk/R/P/nTrax3DvtRcAjNqHf+lHzuTx2y9pVjUlBczsZXcv+ayTVssUCSDTn+Vj//4jHtvy57AHGBo+ytJ12wCOrWpZqM2MxRd3KeylqbS0gkidRpYwLreq5VFyg7YvLP+kumwkEhT4IjWqZV69XjkoUaLAF6nBqsxOHtuyv+ryWsZYokR9+CJVyvRnawr7MWgZY4kWBb5IldZs2l112Y72MTyg6ZYSMerSEalStf3xWsJYokqBL1LCyMDswcEhJuWfgq30EFWbGTfOmaKwl8hS4IsUKDUomx0cYsVTO7nuE52s+9lrDB89/mHFMQYP3KDuG4k+9eGL5N38jRfLDsoODR/huZ8fYs1nLmRcR/ux4+PHtivsJTbUwpfUq/YF4gcHh7jmok6Fu8SWAl9SK9Of5d+e2sE7w+WelT2e5tRL3CnwJZVqfYAKNKde4k99+JI69YT99LNOU1eOxJ5a+JIa9QQ95MJ+89LLw6+QSJMp8CUV5tyzmTf+8F5Nnxk/tp3VV+kl4pIcCnxJtEpLGJeiJ2UlqRT4kkirMjt5fMt+an2fm8JekkyBL4mS6c/yr09u570jtUX92PYx/Me1F6j7RhJNgS+JkenP8vn12zlytLawV6te0kKBL4mQ6c+y9LvbqCXrDfiqljCWFFHgS+zV01+vqZaSRgp8iaVMf5ae3gEGh0Zf/6ZYR/sY7lVfvaSUAl9iJdOf5Qvf28G7h2uZaAkntxlfvv5CBb2kmgJfYqPeJ2U1KCuSo8CXWLj5Gy/ywi9+V9NnxnW003O1npQVGaHAl0ird169WvUiJ1LgS2TVM9USFPYi5YSyPLKZzTez3Wa218yWlzi/1Mx2mdkOM3vGzD4UxnUluTL9Wf5pXW1hf9rJbTy4cJbCXqSMwC18M2sDHgLmAgeArWbW6+67Cor1A93u/o6ZfRb4MrAw6LUlWepd/0bLIohUJ4wundnAXnffB2BmTwALgGOB7+7PFZTfAiwO4bqSIPXMwLn0I2fy+O2XNKhGIskTRpdOJ/Bawf6B/LFybgN+VOqEmS0xsz4z6zt06FAIVZM4UNiLNEdTB23NbDHQDVxW6ry7rwXWAnR3d9f6m73EjObVizRXGIGfBaYU7E/OHzuOmV0BrAQuc/d3Q7iuxFS9QQ8Ke5Egwgj8rcB0M5tGLugXATcVFjCzi4CvA/Pd/c0QrikxdcHqH/P7d4/U/DmtgSMSXODAd/fDZnYHsAloAx5x9wEzuxvoc/deYA3wF8B6MwPY7+5XB722xEeQVv2DWsJYJBSh9OG7+0ZgY9GxOwu2rwjjOhJP9bbqQWEvEiY9aSsN9dEVGzhcx/C7unBEwqfAl4aY+8Dz7HnzjzV9Rm+gEmksBb6Eqp5VLSG3LMI9fzdTYS/SQAp8CcXI2jf1UD+9SHMo8CWwOfds5o0/vFfz504y2Hvv3zagRiJSigJf6hakVX/GKW3suGt+yDUSkdEo8KVmQYIe9LSsSKso8KUmWhZBJL4U+FK1eh+gUl+9SDQo8KUiLYsgkgwKfBnVX67cyJ9qfIH4iFfvU6teJEoU+FJSvUsigFr1IlGlwJcTTF2+oe7PqlUvEl0KfDmm3mURQK8cFIkDBb4ECnrQdEuRuFDgp1yQ7hsFvUi8KPBTJtOfZc2m3WQHh+r+O6afdRqbl14eXqVEpCkU+CmS6c+y4qmdDA3X9/Yp0AwckThT4KfI0nXbOBrg8wp7kXhT4KdEkL56zcARSQYFfsIFeYE4aF69SJIo8BMq6BLGCnqR5FHgJ9C05Ruoc1UEQGEvklQK/ARRq15ERqPAT4ggq1qe2mb8/J4rQ66RiESNAj/mgi6LoFa9SHoo8GMsyFRLPS0rkj4K/BgK2qrXA1Qi6RRK4JvZfOBrQBvwsLvfV3T+FODbwCeA3wIL3f3VMK6dNkFa9Wec0saOu+aHWBsRiZPAgW9mbcBDwFzgALDVzHrdfVdBsduAt9z9o2a2CLgfWBj02mkRdPYNqFUvIuG08GcDe919H4CZPQEsAAoDfwHQk99+EvhPMzN3DzJdPBWChr366kVkRBiB3wm8VrB/AJhTroy7Hzazt4H3A78pLGRmS4AlAF1dXSFULd6CvFcWNANHRI43ptUVKOTua9292927J06c2OrqtEymP8vU5fWH/aUfOVNhLyInCKOFnwWmFOxPzh8rVeaAmZ0EvI/c4K0UCTIoC2rVi0h5YQT+VmC6mU0jF+yLgJuKyvQCtwAvAtcDz6r//nhz7tnMG394r+7PawljEakkcODn++TvADaRm5b5iLsPmNndQJ+79wL/AzxqZnuB35H7oSB5QVr1JxnsvVetehGpLJR5+O6+EdhYdOzOgu0/AZ8J41pJEyTs1X0jIrXQk7YtEMa8eoW9iNRKgd9kcx94nj1v/rHuzyvoRaReCvwm0rIIItJKCvwm0GJnIhIFCvwG07x6EYkKBX6DBJ1Xr6AXkbAp8BtAffUiEkUK/BAFea8sqFUvIo0VqcXT4mzq8g11h70WOxORZlALPyAtiyAicaHAr1PQp2XVoheRZlPg1yHIi0kM+KXCXkRaQIFfIy12JiJxpcCvUpA1cDTVUkSiQIFfhWnLN1DvZEu16kUkKhT4owjytOzZp5/MSyvnhlwjEZH6KfDLqLdVr0FZEYkqBX6RICtbalVLEYkyBX6BC1b/mN+/e6Tmz6n7RkTiQIGPWvUikg6pDvxMf5al67ZxtI7PqlUvInGT2sBfldnJY1v21/VZtepFJI5SF/iZ/iw9vQMMDg3X/NnpZ53G5qWXh18pEZEmSFXgr8rs5PEt+2uebqnuGxFJglQEfr2tei1fLCJJkujAz/RnueuHA7z1Tu3dN2rVi0jSJDbwM/1Zlq3fzvDR2jpwFl/cxZeumdmgWomItE6gVxya2ZlmttnM9uT/HF+izCwze9HMBsxsh5ktDHLNavX0DtQU9uPHtvPgwlkKexFJrKAt/OXAM+5+n5ktz+9/oajMO8Dfu/seM5sEvGxmm9x9MOC1S6qlv96Am9WiF5GUCBr4C4DL89vfAp6nKPDd/ZWC7YNm9iYwEQg98Gvpxhk/tp3VV52n+fQikhpBA/9sd389v/1r4OzRCpvZbOBk4Bdlzi8BlgB0dXXVXJk1m3ZXDHsz+OoNenBKRNKnYuCb2dPAB0qcWlm44+5uZmXT1sw+CDwK3OLuJVczcPe1wFqA7u7umlcnPjg4NOr59jZjzfUXKuxFJJUqBr67X1HunJm9YWYfdPfX84H+ZplyZwAbgJXuvqXu2lYwaVwH2TKh32YKexFJt0CzdIBe4Jb89i3AD4oLmNnJwPeBb7v7kwGvN6pl82bQPsZOON7eZnzlBoW9iKRb0MC/D5hrZnuAK/L7mFm3mT2cL3MD8FfArWa2Lf81K+B1S7rmok7WfOZCxnW0Hzs2fmy7WvYiIoC51/t67sbq7u72vr6+VldDRCRWzOxld+8udS5oC19ERGJCgS8ikhIKfBGRlFDgi4ikhAJfRCQlIjtLx8wOAb8K+NdMAH4TQnWiSvcXf0m/x6TfH0TvHj/k7hNLnYhs4IfBzPrKTU9KAt1f/CX9HpN+fxCve1SXjohISijwRURSIumBv7bVFWgw3V/8Jf0ek35/EKN7THQfvoiI/FnSW/giIpKnwBcRSYlEBb6ZnWlmm81sT/7P8SXKzDKzF81swMx2mNnCVtS1FmY238x2m9ne/Mvii8+fYmbr8udfMrOpza9l/aq4v6Vmtiv//XrGzD7UinrWq9L9FZS7zszczGIxxa9QNfdoZjfkv48DZvadZtcxiCr+jXaZ2XNm1p//d3plK+pZkbsn5gv4MrA8v70cuL9EmXOA6fntScDrwLhW132Ue2oj9w7gD5N7H/B24NyiMv8A/Hd+exGwrtX1Dvn+/hoYm9/+bNLuL1/udOCnwBagu9X1bsD3cDrQD4zP75/V6nqHfH9rgc/mt88FXm11vUt9JaqFDywAvpXf/hZwTXEBd3/F3ffktw+Sey1jyafSImI2sNfd97n7e8AT5O6zUOF9Pwl8ysxOfPVXNFW8P3d/zt3fye9uASY3uY5BVPP9A/gicD/wp2ZWLiTV3OPtwEPu/haAu5d8HWpEVXN/DpyR334fcLCJ9ata0gL/bHd/Pb/9a+Ds0Qqb2WxyP7F/0eiKBdAJvFawfyB/rGQZdz8MvA28vym1C66a+yt0G/CjhtYoXBXvz8w+Dkxx9w3NrFiIqvkengOcY2YvmNkWM5vftNoFV8399QCLzewAsBH4XHOqVpuKLzGPGjN7GvhAiVMrC3fc3c2s7JzT/EvXHwVucfej4dZSGsHMFgPdwGWtrktYzGwM8ABwa4ur0mgnkevWuZzcb2g/NbOZ7j7Y0lqF50bgm+7+FTO7BHjUzM6PWrbELvDd/Ypy58zsDTP7oLu/ng/0kr82mtkZwAZgpbtvaVBVw5IFphTsT84fK1XmgJmdRO5Xyt82p3qBVXN/mNkV5H6oX+bu7zapbmGodH+nA+cDz+d74T4A9JrZ1e4el3d8VvM9PAC85O7DwC/N7BVyPwC2NqeKgVRzf7cB8wHc/UUzO5XcomqR6rpKWpdOL3BLfvsW4AfFBczsZOD7wLfd/ckm1q1eW4HpZjYtX/dF5O6zUOF9Xw886/nRoxioeH9mdhHwdeDqmPX9QoX7c/e33X2Cu09196nkxijiFPZQ3b/RDLnWPWY2gVwXz75mVjKAau5vP/ApADP7GHAqcKiptaxGq0eNw/wi12/9DLAHeBo4M3+8G3g4v70YGAa2FXzNanXdK9zXlcAr5MYaVuaP3U0uGCD3j2s9sBf4GfDhVtc55Pt7Gnij4PvV2+o6h3l/RWWfJ2azdKr8Hhq5rqtdwE5gUavrHPL9nQu8QG4Gzzbg062uc6kvLa0gIpISSevSERGRMhT4IiIpocAXEUkJBb6ISEoo8EVEUkKBLyKSEgp8EZGU+H8MONNF3DlYPAAAAABJRU5ErkJggg==\n",
            "text/plain": [
              "<Figure size 432x288 with 1 Axes>"
            ]
          },
          "metadata": {
            "needs_background": "light"
          }
        }
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "-dw9TGERe3cw"
      },
      "source": [
        "## Useful input features for NN (Neural Networks)\n",
        "\n",
        "The \"apc\" values are typically used as input to the NN for contact cleaning or structure prediction. Though in recent advances (aka DeepMind/Alphafold), the entire MRF was used as the input. More specificially LxLx442. The 442 channels are the 21x21 + (raw and/or apc) value."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "nqbDMKdsYG7Q",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "bc3fee06-1c1b-4d3c-de1f-2c35669f4f8b"
      },
      "source": [
        "w_out = np.zeros((msa[\"ncol_ori\"], msa[\"ncol_ori\"], states * states + 1))\n",
        "v_out = np.zeros((msa[\"ncol_ori\"], states))\n",
        "\n",
        "mrf_ = np.reshape(mrf[\"w\"],(-1, states * states))\n",
        "mtx_ = np.expand_dims(mtx[\"apc\"],-1)\n",
        "\n",
        "w_out[(mtx[\"i\"],mtx[\"j\"])] = np.concatenate((mrf_, mtx_),-1)\n",
        "w_out += np.transpose(w_out,(1,0,2))\n",
        "v_out[mrf[\"v_idx\"]] = mrf[\"v\"]\n",
        "\n",
        "print(\"w_out\",w_out.shape)\n",
        "print(\"v_out\",v_out.shape)"
      ],
      "execution_count": 25,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "w_out (62, 62, 442)\n",
            "v_out (62, 21)\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "_4fGmBsW4q94"
      },
      "source": [
        ""
      ],
      "execution_count": null,
      "outputs": []
    }
  ]
}