Neural_NCPD_twitter.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook contains the Twitter data experiment results where we applied Neural NCPD to the video tensor. Due to the non-convexity of the NCPD task, results will vary run to run, so the results from this notebook will not exactly match those in the paper. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import torch\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "import sys\n",
    "sys.path.insert(1,\"./src\")\n",
    "\n",
    "from NNCPD import NNCPD, weights_H, Energy_Loss_Tensor2, Energy_Loss_Tensor, Recon_Loss, L21_Norm, outer_product, outer_product_np, PTF, random_NNCPD, Fro_Norm\n",
    "from lsqnonneg_module import LsqNonneg\n",
    "from trainNNCPD import train\n",
    "#\n",
    "import torch.nn as nn\n",
    "from torch.autograd import Variable\n",
    "from writer import Writer\n",
    "\n",
    "\n",
    "import tensorly as tl\n",
    "from tensorly import unfold as tl_unfold\n",
    "from tensorly.decomposition import parafac, non_negative_parafac"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "torch.set_default_tensor_type(torch.DoubleTensor)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load in Twitter Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "words = np.load(\"./data/words_100.npy\", allow_pickle=True)\n",
    "X = np.load(\"./data/tweets_bag_100.npy\", allow_pickle=True)\n",
    "\n",
    "\n",
    "X = X.reshape((8, 10, -1))\n",
    "\n",
    "X = torch.from_numpy(X)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Run Neural NCPD"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/jvendrow/.local/lib/python3.6/site-packages/torch/nn/_reduction.py:43: UserWarning: size_average and reduce args will be deprecated, please use reduction='sum' instead.\n",
      "  warnings.warn(warning.format(ret))\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "e73b4f92818540ee8243f9c8aa1d8072"
      }
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "r=8\n",
    "n1,n2,n3 = X.shape\n",
    "net = NNCPD([n1,4,2],[n2,4,2],[n3,4,2])\n",
    "\n",
    "loss_func = Energy_Loss_Tensor()\n",
    "\n",
    "history_unsupervised = train(net, X, loss_func, r, epoch = 15000, lr1 = 0, lr2 = 0.1, random_init=True)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_1 = history_unsupervised.get('A_X1')[-1].detach().numpy()\n",
    "X_2 = history_unsupervised.get('B_X1')[-1].detach().numpy()\n",
    "X_3 = history_unsupervised.get('C_X1')[-1].detach().numpy()\n",
    "\n",
    "A_A1 = history_unsupervised.get('A_A1')[-1].detach().numpy()\n",
    "A_S1 = history_unsupervised.get('A_S1')[-1].detach().numpy()\n",
    "B_A1 = history_unsupervised.get('B_A1')[-1].detach().numpy()\n",
    "B_S1 = history_unsupervised.get('B_S1')[-1].detach().numpy()\n",
    "C_A1 = history_unsupervised.get('C_A1')[-1].detach().numpy()\n",
    "C_S1 = history_unsupervised.get('C_S1')[-1].detach().numpy()\n",
    "\n",
    "A_A2 = history_unsupervised.get('A_A2')[-1].detach().numpy()\n",
    "A_S2 = history_unsupervised.get('A_S2')[-1].detach().numpy()\n",
    "B_A2 = history_unsupervised.get('B_A2')[-1].detach().numpy()\n",
    "B_S2 = history_unsupervised.get('B_S2')[-1].detach().numpy()\n",
    "C_A2 = history_unsupervised.get('C_A2')[-1].detach().numpy()\n",
    "C_S2 = history_unsupervised.get('C_S2')[-1].detach().numpy()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Display Top Keywords for all Layers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "RANK 8 KEYWORDS:\n",
      "------------------\n",
      "Topic 1          Topic 2     Topic 3     Topic 4  Topic 5     Topic 6        Topic 7        Topic 8    \n",
      "-------          -------     -------     -------  -------     -------        -------        -------    \n",
      "martinomalley    marcorubio  johnkasich  crooked  tedcruz     senate         berniesanders  trump      \n",
      "hillaryclinton   teammarco   kasich      hillary  cruz        florida        people         hillary    \n",
      "realdonaldtrump  vote        ohio        thank    ted         zika           bernie         donald     \n",
      "campaigning      flsen       john        great    internet    venezuela      must           president  \n",
      "maryland         click       gov         clinton  choosecruz  nicolasmaduro  change         timkaine   \n",
      "\n",
      "RANK 4 KEYWORDS:\n",
      "------------------\n",
      "Topic 1     Topic 2  Topic 3     Topic 4     \n",
      "-------     -------  -------     -------     \n",
      "tedcruz     trump    marcorubio  johnkasich  \n",
      "cruz        hillary  teammarco   kasich      \n",
      "ted         people   vote        ohio        \n",
      "internet    donald   flsen       john        \n",
      "choosecruz  vote     florida     gov         \n",
      "\n",
      "RANK 2 KEYWORDS:\n",
      "------------------\n",
      "Topic 1  Topic 2     \n",
      "-------  -------     \n",
      "trump    tedcruz     \n",
      "hillary  cruz        \n",
      "vote     ted         \n",
      "people   internet    \n",
      "donald   choosecruz  \n"
     ]
    }
   ],
   "source": [
    "keywords = np.empty((7,8), dtype=object)\n",
    "\n",
    "for i in range(keywords.shape[1]):\n",
    "    keywords[0,i] = \"Topic \" + str(i+1)\n",
    "    keywords[1,i] = \"-------\"\n",
    "    \n",
    "C = X_3\n",
    "\n",
    "for i in range(C.shape[1]):\n",
    "    col = C[:,i]\n",
    "    top = col.argsort()\n",
    "    top = top[-5:][::-1]\n",
    "\n",
    "    keywords[2:,i] = words[top] \n",
    "\n",
    "print(\"RANK 8 KEYWORDS:\")  \n",
    "print(\"------------------\")\n",
    "col_widths = [max([len(keywords[i][j]) for i in range(keywords.shape[0])])+2 for j in range(keywords.shape[1])]\n",
    "for row in keywords:\n",
    "    print(\"\".join(row[i].ljust(col_widths[i]) for i in range(len(row))))\n",
    "    \n",
    "\n",
    "C = C_A1\n",
    "\n",
    "keywords = np.empty((7,4), dtype=object)\n",
    "\n",
    "for i in range(keywords.shape[1]):\n",
    "    keywords[0,i] = \"Topic \" + str(i+1)\n",
    "    keywords[1,i] = \"-------\"\n",
    "\n",
    "for i in range(C.shape[1]):\n",
    "    col = C[:,i]\n",
    "    top = col.argsort()\n",
    "    top = top[-5:][::-1]\n",
    "\n",
    "    keywords[2:,i] = words[top]\n",
    "\n",
    "print(\"\")\n",
    "print(\"RANK 4 KEYWORDS:\")  \n",
    "print(\"------------------\")\n",
    "\n",
    "col_widths = [max([len(keywords[i][j]) for i in range(keywords.shape[0])])+2 for j in range(keywords.shape[1])]\n",
    "for row in keywords:\n",
    "    print(\"\".join(row[i].ljust(col_widths[i]) for i in range(len(row))))\n",
    "\n",
    "\n",
    "C = np.dot(C_A1, C_A2)\n",
    "\n",
    "keywords = np.empty((7,2), dtype=object)\n",
    "\n",
    "for i in range(keywords.shape[1]):\n",
    "    keywords[0,i] = \"Topic \" + str(i+1)\n",
    "    keywords[1,i] = \"-------\"\n",
    "\n",
    "for i in range(C.shape[1]):\n",
    "    col = C[:,i]\n",
    "    top = col.argsort()\n",
    "    top = top[-5:][::-1]\n",
    "\n",
    "    keywords[2:,i] = words[top]\n",
    "\n",
    "print(\"\")\n",
    "print(\"RANK 2 KEYWORDS:\")  \n",
    "print(\"------------------\")\n",
    "\n",
    "col_widths = [max([len(keywords[i][j]) for i in range(keywords.shape[0])])+2 for j in range(keywords.shape[1])]\n",
    "for row in keywords:\n",
    "    print(\"\".join(row[i].ljust(col_widths[i]) for i in range(len(row))))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Display Rank 8 Factor Matrices\n",
    "We note that the columns of the factor matrices could be permuted, topics could be in any order."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAycAAAFwCAYAAACivQTPAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABPWklEQVR4nO3dd7hcVdn+8e+dRkInBAQEiYAgCkgJHSQ0UUSKFEEEoigCllfFhliCvPoi2BARzQ8hSJEmUhVFISAgJSCoSFRKKKH3kp48vz/WmmQYppw6e8859+e65tpn9l577WcOZO/zzGqKCMzMzMzMzIo2pOgAzMzMzMzMwMmJmZmZmZmVhJMTMzMzMzMrBScnZmZmZmZWCk5OzMzMzMysFJycmJmZmZlZKQwrOgDrvhFaIkayVNFhWBto+PCiQ2gq5s0rOoQBQcOGFh1CS0PXKTqC1l6Y9uyzEbFS0XHY640YOipGDV220BgWvlWFXn/u3OL/3NLcYn8HAEsvO6voEJh3f7Hfy88fW/z9fsHC4tsm5jz4eMP7dfH/WqzbRrIUW2rnosNobkjx//haioVFR9DSsDetWnQITc2f8XjRIQwIQ5cfXXQILS1/ZvnXxLpwm//3cNEx2BuNGros26xyUKExzPxZsX/uTH+0+Jx55CMjig6Brd77j6JD4Km9RhV6/WdOWb7Q6wO89EqxvwOABw78RsP7dfGpk5mZmZmZGU5OzMzMzMysJJycmJmZmZlZKTg5MTMzMzOzUnByYmZmZmZmpeDkxMzMzMzMSsHJiZmZmZmZlYKTEzMzsw4jabqkaPHau4d1j83nT+/bqM3MWvMijGZmZp3rD8CTDY490s5AzMz6gpMTMzOzznViREwpOggzs77ibl1mZmZmZlYKhSQnknaXdI6k+yW9KmmOpMckXS3pSEnLVJWdmPu+TqypY3zeP6Xd8ZuZmXUSSbtJukLSU5LmSnpC0q8lbdjivGGSvirpPkmz8/lnS3pLu2I3s8GlrcmJpJVzMnE18BFgLvBH4LfAdGAX4HTgQUlrtjO2HF9IinZf18zMrL9IOgW4Bngf8ABwGfAEcCBwu6Tdm5x+IXA8afzKZcAc4FDgDknr9V/UZjZYtW3MiaTlgZuBdYC/AkdGxN9ryiwDHAUcB6wAPNykytuB9YGZ/RGvmZlZp5N0JPBZ4F5gv4iYVnVsb+Bi4DxJa0XECzWnrwmMAjaJiH/lc0YAvyR9wXgOsEW/fwgzG1Ta2XLyU1JicjuwU21iAhARr0TEScBmwFPNKouImRExLSI8G4mZmQ1W1zeYRniypKHAN3O5A6oTE4CIuAz4BbA8Kdmo54RKYpLPmQt8BngZ2FzStvVOknSEpKmSps5dOKtXH9DMBpe2JCeS1gYOym+PjIjZzcpHxP0R8USLOuuOOamen13J0ZLuljRT0guSLpe0Qc05E6u7c9Xe5GvKStIhkqbk+mZLekDSaZLWaBDronokfUjSX/NYm1ck/VnSds0+q5mZWQN/AM6u87oJ2BhYFbi3OsGocUPebt3g+Lm1OyLiReDK/HZ8vZMiYlJEjIuIcSOGjGr5IczMKtrVrWsPUiL0j4j4W5uuCTAZ+BBwI/BfYHNgT2C8pE0i4sFc7m7Szfyw/P7sepVJEulG/WFgHjAFeJ7UrH00cKCk90bEHQ3O/zapy9pNpHE3GwE7AdtJGh8Rf+3FZzUzs8Gn4VTCkvbPP76zC+MpV6qz78WciNQzPW9XbxWgmVl3tCs52Sxv6/7R3k/WBLYH3hkRDwBIWgK4FNgdOBb4BCxq2r5M0mH5/YQGdR5FSkyeAnaOiHtzvUOBH5Gaui+WtF5EzKlz/qeALSLiznzeEODnOY5vA7v27iObmZktMjRvZwB/alF2WovjZmZt0a7kpPKNzNNtul7FZyuJCUBEzJF0PCk52bkH9R2Tt9+oJCa53gWSvgjsRUqK9gPOq3P+tyqJST5voaRvkJKT7SUNj4h59S4s6QjgCICRLNmD0M3MbJB5NG+faPKlWzPLS1ouIl6qc2xs3s7oSWBmZo0M5EUY55OmTqxV+XZote5UJml1YC1gIWmGktfJgwQrCcn4BtVcVee8p4AXgCWAFRtdv7r/7nCW6E7oZmY2ON0OPAdsImmdHtZxcO0OScuRumtD6t5sZtZn2pWcPJO3K7fpepC+KZpfuzMiXs4/dvcv/DdX1dtoQP+DNWVrNZpZrBLTyG7GZGZmVlduiT+B1L3rMklvmPZX0ghJe0p6e4Nqvilp/aryw4FTgOWAOyPipn4I3cwGsXZ167oTOIQ0IL1dFvZTvT1epDEi+ismMzOzN4iIU/Kixp8HbpP0d9JCjHNJX6RtAixFWqCxdtzJI6Tn992SrgNeArYB1gCeJS3GaGbWp9rVcnI1KVnYUNImbbpmX6v0q10tD6yvZ62asmZmZoWKiC8AOwAXkBY4fj/wXmAMqbvxwcBf6p0KHEBqfVkL2JvUwn8usHmT6YnNzHqsLS0nEXG/pAtJa52cLmmHBrNZAYvWRZnZaq2TfjAPGC5pWG2XsIh4TNKDpBv0R0gr5C6Sm7orfXOntCFWMzMbpCJibDfL30iaVr8rZacDqtr1v/llZtbv2jkg/tOkMRlbAtdJ2rC2gKSlJH2B1Iz8pjbGVlFp8Vi/wfEf5u0J1f1z81TCJwFvAR4GLum3CM3MzMzMBqh2jTkhIp7PK6FfBGwH/F3Sv0h9XCt9X7cgDVR/irS4Ybv9ltQv98+5f+2rOfaP5+M/A7YltQDdk1enryzCuBZp1q39m7UKmZmZmZlZfW1LTgByN63tJe1B+gN/a1K/12GkGb3+BFwOnB8Rr7Uztuw4Uh/bfYAPAsPz/o8DRERIOhj4PWltki2BUcDjwOnA/0XEo7WVmpmZmZlZa21NTioi4irqrPnRoOxEYGKd/VN4fZ/Yyv7p9fbXlKl7PCJmkRZaPKbe8VwmSOucvGGtk+5er+r42K7WZWZmZmY2UA3kRRjNzMzMzKyDODkxMzMzM7NScHJiZmZmZmalUMiYEzMzMxskhghGDG9drh+NXea5Qq//yPDRhV4fYMHIKDoEnpm9dNEhwLKjCr38s08vW+j1AdZd88miQ+CBJsfccmJmZmZmZqXg5MTMzMzMzErByYmZmZmZmZWCkxMzMzMzMysFJydmZmZmZlYKTk7MzMzMzKwUnJyYmZmZmVkpeJ2TDqShQxm63ApFh9FUzJ1bdAgtLZw1u+gQWpo/4/GiQ+h4T35+m6JDaGmVH91SdAgt3XrXlkWHMOhJmg6smd9+LyK+2qTsucDB+e0NETG+f6MzM+sbbjkxMzPrPIdIGlrvgKRlgQ+2OR4zsz7h5MTMzKyzTAVWA3ZtcPxAYBRwR9siMjPrI05OzMzMOsvkvJ3Q4PgEYAFwThtiMTPrU05OzMzMOsttwH3AXpKWrz4gaT1ga+APwBP1Tpa0i6TTJN0j6TlJcyQ9LOlsSes3OGeypJA0QdJGki6W9KSkBZI+16efzswGNScnZmZmnecsYCRwUM3+CVXHG/k5cDgwH7gR+B0wFzgUmCppuybnbgvcDmwKTAGuAWZ2L3Qzs8acnJiZmXWec0hdtyZUduQB8ocCzwNXNDn3i8AqEbFZROwTEfsA6wJHAksCkySpwbkfB74PvC0iDoyI90fEpF5/GjOzzMmJmZlZh4mIJ0mtFltUdcV6D2mg/PkR0XA+94i4LCJerNkXEfEL4BZgfeAdDU6fBnwrIhb28iOYmdXldU7MzMw602Tg/aTWk6+wuBVlcqsTJa2ez307sCxQmZZ4lbxdF7i3zqmXR8SCLtR/BHAEwMhhy7Qqbma2iJMTMzOzznQF8BxpzZOTgb2Af0TEnc1OknQ88DWa/w2wbIP9D3clsNzVaxLAciNXia6cY2YG7tZlZmbWkXLXrfOBVUkD4Jeg+UB4JO0LfBOYBXwCWBtYMiIUEQJ+XSnaoIpZfRC6mVlDTk7MzMw61+S83YM0+9Z5Lcrvn7dfi4gzIuLBiKhOONbp4/jMzLrFyYmZmVmHioi7gJtJ3bsujoinW5wyOm8frT2QB9Zv0rcRmpl1T0clJ5Km50Wgxjc4vqGkJ3KZSySN6OF1JuQ6JvciXDMzs34XEdtFxJiI+HAXik/L209UPyMlrQycjceimlnBOio5aUbSlsANpJlGzgI+1GwqRTMzs0Hox8BLpJm67s8rvV8FPAAsDVxWXGhmZgMkOZG0E/AnYAXSjffwrkx12MRvSfO8H9v76MzMzMohIh4kdd26gDTo/QOk590kYGtS4mJmVpiOb76VtCdwEWmWkuMjYmJv64yIl/AN2szMSiQixnaz/CXUmXUrIh4CDmpw2gSqVp2vOqfufjOzvtbRLSeSDgZ+A4wAPlcvMZG0paSTJU2V9JSkuZIez2NStmpQb90xJ5LG5/1TJA2XdJykaZJmS3pa0rmS3tIk3jUknSLp35JmSXpZ0s35eo2mbTQzMzMzGxQ6NjmRdBRwDulboY9FxCkNin4H+DwwHLidxYtW7QvcJGn/Buc1Mxz4PfBV4P7880Lg4Fzn8nXi3RH4B/BZ0u/9GuA2YCPSGJmzexCHmZmZmdmA0anduj5HWgl3LnBQRFzapOz3gYMj4qnqnZI+QGp1+bmkqyNiZjeuvw0wFVi7Mm2jpOWA64BNgU+RkqLKtVbN11qa1Cz+q4iIfGwNUsJ0iKTrImJyN+IwMzMzMxswOrXlZK+8/XGLxISIuKY2Mcn7rwQuJs35vmM3rx+k1ppF88nncSrfy293rin/OdJg/R9ExNmVxCSf9yhplV6AzzS6oKQjcte0qXNjdjfDNTMzMzMrv05tObkB2AE4RtKdEXFRs8KSxpBWz90AWJ7Fn3uDvF0XuLob138kIv5RZ39l/vjVavbvnrcXN6jvTuBVYGNJIyPemH1ExCTSbCosN2ylqD1uZmZmZtbpOjU5mUiao/2LwHmSaJSgSPok8ENgySb1LdvN6z/SYP/LeTuyZv9aeXtHF8a9rwjM6GY8ZmZmZmYdr1OTEyLiS/kP/YYJiqTNgdOB+cCXgCuBx4CZERGSvktay6S7M2Ut7Gb5oXl7IdCqT9acbtZtZmZmZjYgdGxyAnUTlIiI6q5T+5ISj59ExPfrVLFOG8IEeDRf64SIuLdN1zQzMyvc/KWH88z2qxYbw1bTC73+wp8PL/T6AFqj+PGqL84eVXQILL38UoVef7XVni/0+gBDh3T3O/b26tQB8YtExJdIM3INA86vmRp4dN4+WnuepJWAXfs/QiBNNQzQk2mLzczMzMwGhY5PTqBugrJfPlQZoH6opKUr5SUtA5xJGhzfDieTxqN8TdKnJL2hxUrSOyV9sE3xmJmZmZmVzoBITuANCcqvc4JyFqnVZFPgQUmXSvotMB0YR0pQ2hHbo8DewCvAT4FHJF0r6TxJV0t6BPgncEA74jEzMzMzK6MBk5zAGxMU0noj40hT8L5KmuFrHHApKWF5Q3evfozteuCdwHeBp4GtSGNi3gk8SBqYf1y74jEzMzMzK5uOGhAfEWO7UOZLpJm5qn2yQfGJ+VVbx2Rgcp39U2gys1dETG9x/ElSAuIkxMzMzMysxoBqOTEzMzMzs87l5MTMzMzMzErByYmZmZmZmZWCkxMzM7MOJmmIpEckhaRnJBW/4p+ZWQ85OTEzM+tsuwJr5J/HAHsWGIuZWa84OTEzM+tsH8vbGTXvzcw6jpMTMzOzDiVpNLAXEMCBwAJgN0mrFRqYmVkPOTkxMzPrXAcDSwBTIuIm4I/AUOCweoXzuJTIPx8h6W+SZkp6TtKlkjbownmHS7pN0st5//L98cHMbHBycmJmZta5Kl24JuftWXn70WYnSfoRcDrwEnA58CywD3CbpO2anHcqMAmYA1wF3ElqtTEz6xMdtUK8mZmZJZI2ATYGXgEuybuvAJ4H3iZp+4j4S4PTjwB2jIgbc10Cvgt8FThf0roRMbvOeYcAW0fE7X33SczMFnPLiZmZWWeqtJpcFBEzASJiDnBezfF6Tq8kJvm8AL4OPEia+WvfBued5MTEzPqTW0460PC3BauePa/oMJp6bOuZRYfQWpS/J4KGlfufaMyfX3QILa3yo1uKDqGlISNHFh1CS+t+/s6iQ2jp4aIDaCNJSwAfzm/Pqjl8FvAZYH9Jn4mIV+tUcW7tjohYIOnXwHHAeBYnOdUu7WJ8R5BaZxix1ApdOcXMDHDLiZmZWSfaGxgN/Dcibq4+EBF/A+4BlgI+1OD8hxrsn563qzc43qUcMCImRcS4iBg3bORSXTnFzAxwy4mZmVknqnTZWk7STXWOr1xV7pd9ddGImNVXdZmZ1ePkxMzMrINIWgPYJb9dmcWJSD3bSFovIv5ds38sqXWl1ti8nVHnmJlZv3O3LjMzs84ygfT8vi4i1OgFXJTL1xsYf3DtDklDSQs5Akzph7jNzFpycmJmZtYh8pS/E/Lbc1oUrxw/JCce1Y6uXs8k13s8sDap1eQ3vY/WzKz73K3LzMysc4wH1gJm0TqBuAZ4BlgV2B24surY/wNukHQj8ASwKbBervdgjy0xs6K45cTMzKxzVLpoXRYRrzQrGBHzgQtqzqv4Amm64dGkmb9WBi4DtoyIG/oqWDOz7nLLiZmZWYeIiENIq7R3tfxngc/W2R/Az/KrK/Woq9c0M+sNt5yYmZmZmVkpODkxMzMzM7NScHJiZmZmZmal0K3kRNIQSR+WdLmkGZLmSHpe0lRJJ0iquxCUpImSQtKEmv0T8v6QNEvS8k2u/baqsiFpfHdib1JvSIo6+6f05XXMzMyKVrUGiplZKXU5OZG0OnA7cB6wB/AIcClwC/BW4OvAA5IO6GEsI1m8+FM9E3pYr5mZmZmZdYAuJSeSRgN/ATYjrRq7TkRsHREHRcQewCrAscCSwAWS9u1mHH8DFtAgAZE0BDgUeAG4v5t1m5mZmZlZB+hqy8lpwFjgDuB9EfFQ9cGImBcRJ5LmTRdwpqQx3YjjceBaYEtJb69zfBdgddJ87XO6Ua+ZmZmZmXWIluucSFobqHTVOjoiZjcp/hPgcGBD4NPAxG7EchbwXlLryVdrjn20qsy7G8S5EvBh4H2kVW5XJSUy9wHnAD+PiAXdiKcpSbsBnwK2BFYAniO1Kn03Iv5RVW574EZgWkSs36CuMcBjQACrR8RzfRWnmZlZkYbODZZ9uNjvFTV8RKHXX/HOoYVeH0B7vVh0CDz78lJFh8AR5/y+0OtPvH3PQq8P8PIyzf6UL15XWk72yOXujYipzQrmRZ1+ld9297d/Oanb1iGSFv0rzoPk987Xv6PJ+bsBPwbeCTwE/Ba4C9gY+CnwG0l9MghQ0inANaRE6AHSqrpPkMbM3C5p90rZiPgLcA/wdkk7Najy48ASwAVOTMzMzMxssOpKcrJZ3t7exTorCcS7JA0DiIiJeYaQyY1Oiog5wK+B1YD3VB06kDRYvuG52Z3AVhGxRkTslMfD7EgarH83sBeLW4B6TNKRpNV27wU2jIhtIuKAiNgU2IfUGnWepBWqTjs1b4+uU98Q4JP57Wm9jc/MzMzMrFN1JTlZKW+f6mKdlXJDgNHdjOesvJ1Qte+jwHxS16yGIuK+iLitzv4ngC/nt/t1M57XyS0638xvD4iIaTXXugz4BbA88JGqQ+cDzwN7SVqtptr3k8fztGqZMjMzMzMbyPpjEcYed53Kf5z/k/RH/AqS1ge2AK6JiJbJkaRhkt4j6VuSTpd0lqTJwJG5yLo9jS3bmDSW5d6I+FeDMjfk7daVHRExCziD1KpyRE35SmtK01YTSUfk9WSmzn6x3H0FzczMzMx6ouWAeODZvH1TF+usLMS4kNRa0F2Tge8DB5FaFGBxi0pDktYljf2oO+g8W7YH8VRbK2/fWW/hxhor1bw/DTgG+ISk/42I+Xmygd1Ig+kvbFZZREwCJgGMWX9Mq2ubmZmZmXWcriQnd5K6KG3VxTq3yNt7ImJ+D2I6FzgR+Bhp/MmzwJVdOO8SUmJyBXASaZaulyJiQU5c/k0vWnWyykD9GcCfWpSt7fL1iKQrSONS9s7xHpVjOrPFLGhmZmZmZgNeV5KTq4AfAOtL2rzZjFl5NqxD89uuJBRvEBFPSfo98IG86ycRMa/ZOXltlA2Bp4EP1pkyeJ2exFLHo3n7RERM6MH5p5KSk6MlXUUaT7MQOL1vwjMzMzMz61wtx5xExP2kb/kBTpM0sknxzwIbAK/Qu5mnziB1dXqOLnTpYvHA+8cbrGVycC9iqXZ7jmkTSd1OeCLietKYmh2B40lx/752UUszMzMzs8GoqwPiP0VqNdgc+J2ksdUHJQ2X9BXgh6SFBD8eEU/3NKiIuCIixuTX3V045b+kFogNJL1ukUZJHyWNX+m13IJzAql712WStqgtI2mEpD0brHQPac0VWDyD2M/6IjYzMxtYJE2XFPl1Youy51aVndKmEM3M+lyXkpOIeBbYjrReyI7A/ZJulnS+pCuBJ0njRGYBB0fERf0Ub6P4niH9kT8MuF7SdTm2fwBn5tj66lqnAD8iLfZ4m6R7JF0q6QJJfyFNAnA5iwfz1zqHtNgkwIOkxRzNzMyaed0CxdUkLQt8sM3xmJn1iy5PJRwRjwDjgEOA35MWN9yPlLRMB74DrBMRv+77MLvkf0jT9N5DGpT/PtKaK+8jz3LVVyLiC8AOwAXACqS1St4LjCGN0TkY+EuDc2cCt+S3p0fEwr6MzczMBpyppAlidm1w/EBgFIsXQTYz61hdGRC/SB7PcW5+9VpeMX5yN8/ZoMH+hcD/y6966s7UFRGN9o9vEceNwI3NytQNQlqZ9ICZRWrVMTMza2Yy6cvBCdRvbZ8ALCC1zG/erqDMzPpDfyzCaM0dB4wAzo6InqwDY2Zmg8ttpOnx95K0fPUBSeuRFv39A/BEvZMl7SLptNwN+TlJcyQ9LOnsvNhxbfk/57ErBzYKSNIPcpmTevPBzMxqOTlpA0nbSDpD0vWkGc1eIg2sNzMz64qzgJG8cYKXCVXHG/k5cDgwn9Ti/ztgLmnq/6mStqspf2reHl2vMkmj8FT4ZtZPnJy0x7qkB8MWwE3AeyPi8WJDMjOzDnIOqevWhMqOPED+UNJELFc0OfeLwCoRsVlE7BMR+5CeS0cCSwKT8jplFVcCDwPbS6rXlfog0nhLT4VvZn3OyUkbRMTkiFBELBUR20fErUXHZGZmnSMiniSNN9miqivWe0gD5c+PiLlNzr0sIl6s2RcR8QvSBC3rA++oOraAxdPc12s9qexrOBW+pCMkTZU0de7c15p+NjOzak5OzMzMOsPkvJ1Qs51MC5JWl/RJST+S9EtJkyVNBlbJRdatOeUM0sQtH5G0TFU9WwKb0WIq/IiYFBHjImLciBFLtQrPzGyRbs3WZWZmZoW5AniOtObJycBewD8i4s5mJ0k6HvgazZ/5y1a/iYjnJZ0HfJy0hEBtS4qnwjezfuGWEzMzsw6Qu26dD6xKGgC/BM0HwiNpX+CbpFaQTwBrA0vmrsYCKmuT1ZtWvzIw/qhc14rAAcBsPBW+mfUTJydmZmadY3Le7kGafeu8FuX3z9uvRcQZEfFgRMyqOr5OoxMj4u+k2b02kPRu0sQuI4ELPBW+mfUXJydmZmYdIiLuAm4mde+6OCKebnHK6Lx9tPZAHli/SYvzK60nnybN7gVwWteiNTPrPo85MTMz6yARUbsuSTPTgF2BT0j6fWVWL0krA2fT+u+Ay0iJTaUF5o6ImNq9iM3Mus4tJ2ZmZgPXj0kL/74fuF/SxZKuAh4AliYlHw1FxHxev9CiW03MrF85OTEzMxugIuJBUtetC0iD3j9AWtdkErA1KXFp5dq8fQ64sB/CNDNbxN26OtDcaQt5bOtyL2qloUOLDqGlWLCg6BBaivnziw6h4w0ZObLoEFra4JY5RYfQ0j+3X7LoEFqbV3QAfSsixnaz/CXUmXUrr+J+UIPTJlC16nwDB+ftLyNidndiMjPrLrecmJmZWV2S1iBNQTwX+GnB4ZjZIOCWEzMzM3sdSScCq5MG0y8FnBwRb5jxy8ysrzk5MTMzs1oHAm8BngC+B3yj2HDMbLBwcmJmZmav093xLmZmfcVjTszMzMzMrBScnJiZmZmZWSk4OTEzMzMzs1LwmBMzMzPrN/OWEk9uWex6Q295fI1Cr//aqm9YfqbthswZXnQILDGi+LW7Hpu7YqHXj5nF/+k9YoXi/zs045YTMzMzMzMrBScnZmZmZmZWCk5OzMzMzMysFHqdnEh6u6SfS/q3pJmSZkl6RNItkn4gade+CLQ/SJooKSRNLDoWMzMzM7PBrlejciR9CPgVMAKYAUwBXgBWAjYFtgZ2AK7tVZRmZmZmZjbg9Tg5kbQKcCYpMfk8cGpELKg6PgTYLr/MzMzMzMya6k3LyR7AksBfI+LHtQcjYiFwY36ZmZmZmZk11ZsxJyvn7dPdOUnSlpJOljRV0lOS5kp6XNIlkrZqcM6isSGS3iTpF5IekzRH0kOSTpRUdxJ1ScMlfVHSvyTNlvSkpHMkrdnFWC/I15or6RlJV0iq2xqUY4z88+GSbpP0ct6/fN6/mqSfSro/xzMzj9G5RtIRXf5FmpmZVZE0JT9vJhQdi5lZT/UmOXkkb3eWtEE3zvsOqRvYcOB24ArgOWBf4CZJ+zc5dw3gTlKrzV9JY1xWBr4CXFRbOHctuxQ4GXgrcB1wA7BzruetjS4k6Zh8jQOAJ4HLgfuB9wM3SPpEk3NPBSYBc4Cr8rVC0qr550+RWq2uAa4k/S63Ar7Q5LObmZmZmQ1ovenWdTnwOLAa8DdJfyT94X8XcEdEvNTgvO8DB0fEU9U7JX0A+A3wc0lXR8TMOud+DDgD+FREzM3nrU9Kcj4gaduIuLmq/KdIicwMYHxE3J/PGQmcCxxaL0BJ78txPg58MCJuqzq2LfA74DRJN0TEf+pUcQiwdUTcXlPv/wCrAL8AjoqIqDq2BLBlvXjMzMzMzAaDHrecRMQrwC7AVFKSszvwPdLMXM9LujnP5lV73jW1iUnefyVwMTAa2LHBZR8FPltJTPJ59wHn5Lc715T/XN5+vZKY5HNmA0cDsxpcZ2Lefrw6Mcnn3gycQGr5+WSD80+qTUyyN+XtNdWJSa53TkR4fI6ZmZmZDVq9WuckIu6LiM2BbYHvAn8mTSU8BNgGuEDS5NrzJI2RNEHS9yWdIWlyLlfpHrZug0teFxH1Eoppebta1TVWB9YCFgLn14n9aeCP9WIDtgBernc8uyFvt25w/NIG+ysJy/ck7S1pqQbl3kDSEXmcztR5zOnqaWZmNojl52vDcSiN1vvqzVhPM7Pe6NU6JxURcQtwCywa57EV8C3gPcBhuZvWxfn4J4Efkmb6amTZBvsfabD/5bytvlGunrePV7e01JheZ19lHMqywHxJjWKEtJ5LPQ832H8O6XfyYeC3wAJJ/yTNaHZB/j3WFRGTSONYWFajo1E5MzOzPlQZ6ynSc35Z0hIBXwHeAexZXGhmNhD1SXJSLU8hfIuk3UktBZsCewMXS9ocOB2YD3yJNBj8MWBmRISk7wLHkm6C9Szs63jrGJq3LwGXtSj7bL2dDVp3Kr+bgyX9H2kszLb59RngM5LOjIjDexK0mZlZP+juWE8zs17p8+SkIiIWSLqOlJxUWhj2JSUeP4mI79c5bZ0+DGFG3q4maUSD1pOxdfY9mrfzImJCH8azSET8E/gnLGpp2p3U9exjki6MiEbdyczMzNqp7lhPSecAR5HGejo5MbM+0+MxJ2rR3yl7S94+lrej8/bR2oKSVgJ27Wk8tSLiUeAh0mc8sKvXi4gZwD+AMZLG91U8jUTEwoi4ijT7GcC7+vuaZmZmXdTlsZ7VqsdJzp/5Wv9FZ2YDTm8GxB8t6SxJW9QekDQsrwOyX951Yd5WbmaHSlq6qvwywJnA8r2Ip56f5O3/Slqr6npLAKfReNzLN/L2XEnvqT0oaaikndRg0chGJB0qadM6+1dk8eD6RuNVzMzM2q07Yz0XiYhJETEuIsYNW7LLc7+YmfWqW9dwYAIwQdKTwN3A86TWkY1Y/G3KSRHxh/zzWaTpfTcFHpR0E6mb17uBuaQE5WO9iKnWqaQB6O8D7s3dzF4lDeYbCfyKOmudRMTleRHGk4A/SPoP8O987irAJqRE6ijg1m7E80HgbEkzSL+vF4EVge2BpYC/kAbKm5mZtUOrLynbMdbTzGyR3rSc/BLYh9QC8SiwIbA/KdF4FTgb2D4ivlI5ISJeAMaRZp16lbTa+jjS1LubUqe7V29ExAJgL+CrpJm5diGtoXJjvu5DTc79IbAZ6XMOJXUB+wBpFrAbgU9QZ1X6Fn4AnEJa3HEc6fe1EWnhysOBXSNiXjfrNDMza6QyVmTpBsfXbFcgZmZd0eOWk7wI42W0ntGq9rynabx44UQWL4BYfU7d/VXHJwOTGxybR1oc8ntdvV7VufcAH290vE75puNwIuIvpNYRMzOzdqhMDvP22gOSRgHj2xqNmVkLvVqE0czMzErtz3l7iKT1KjtzYnI6iyeuMTMrBScnZmZmA8uicSIRcRNwFWnxxLskXSPpKlK35veQxoKamZWGkxMzM7OBYVTe1s7duz9wIvA0sBNpjOfVedtoNi4zs0L02yKMZmZm1h557bHKlPnTq49FxGzg2PyqNZE+HutpZtYbbjkxMzPrfIcBY4BngHsKjsXMrMfccmJmZtaBJC0J/AJYm8UL+X4jIuYXF5WZWe84OTEzM+tMI4CPAK8ANwOnRMTFxYZkZtY7Tk7MzMw6UES8CDRdX8vMrNN4zImZmZmZmZWCkxMzMzMzMysFJydmZmZmZlYKHnPSgTRkCENGjWpdsEAaNbLoEFpa8NzzRYfQ8Tb+W9ERtPaPHUYUHUJLf990dtEhtDR0hfL/Hq2cll3hNXY74NZCY/jn9xYWev15y6xc6PUBxoyaU3QIPPXE8kWHwFrrP13o9ces/mKh1wd4ZeYSRYfQlFtOzMzMzMysFJycmJmZmZlZKTg5MTMzMzOzUnByYmZmZmZmpeDkxMzMzMzMSsHJiZmZmZmZlYKTEzMzMzMzKwUnJ2ZmZmZmVgpOTszMzAY4SVMkhaQJfVDXhFzX5N5HZmb2ek5OzMzMzMysFJycmJmZmZlZKTg5MTMzMzOzUhgUyYmk3SWdI+l+Sa9KmiPpMUlXSzpS0jJFx2hmZtZOkiY3G4ciaWI+PrG9kZnZYDas6AD6k6SVgYuAHfKu+4A/AnOB1YFdgN2BEySNi4iHCwnUzMzMzMwGbnIiaXngZmAd4K/AkRHx95oyywBHAccBKwBOTszMzMzMCjJgkxPgp6TE5HZgp4iYXVsgIl4BTpJ0KfBam+MzMzMzM7MqA3LMiaS1gYPy2yPrJSbVIuL+iHgin7uoD66kjSRdLOlJSQskfa62TIPrv6GfrqTxeV+r19je/wbMzMyKI+kISVMlTZ31wpyiwzGzDjJQW072ICVe/4iIv/Wwjm2BnwMzgCnAMsDMXsT0JHB2g2OrALvlnxf04hpmZmaFi4hJwCSAN71jdBQcjpl1kIGanGyWt3f0oo6PA98BvhkRC3sbUERMAybU7pe0NHBjfvvDiHi0t9cyMzPrAwOyd4WZldtAvfGslLdP96KOacC3+iIxaUTSUOBCYBPgN8CX+utaZmZmNebm7dINjq/ZrkDMzCoGastJX7g8Ivq7i9VPSVMZ3wp8pFkiJOkI4AiAkVqqn8MyM7NBYEbevr32gKRRwPi2RmNmxsBtOXkmb1fuRR39Oq2wpC8DRwIPAHt2YdD+pIgYFxHjRmhkf4ZmZmaDw5/z9hBJ61V25sTkdOAthURlZoPaQE1O7szbzXtRx6xenNv09yppf+BE4Dlg94h4pll5MzOzPrKohT4ibgKuApYF7pJ0jaSrgIeA9wBnFROimQ1mAzU5uZp0A95Q0ib9UH+P++lK2gb4Va5jr4j4Tx/HZmZmVmtU3tau6VX5suxpYCdgU9IzdFPgkbZFZ2aWDcjkJCLuJw00Bzhd0hLNyktaW9Kq3bhEj/rpSloHuBxYAjg0Im7uxjXNzMy6TZKAtfLb6dXHImJ2RBwbEW+NiBERsVpEHB4RT0bExIhQREysOWdy3j+hLR/AzAaVAZmcZJ8GHgS2BK6TtGFtAUlLSfoCqRvYm7pRd7f76UpaEfg9MAb4akRc1I3rmZmZ9dRhpGfPM8A9BcdiZtbUgJ2tKyKel7QdcBGwHfB3Sf8iTRE8F3gzsAWpFeMp4Plu1H1T7pe7B6mf7l+A+cA4Uneys4CP1px2NLAO8CrwDkmTG1T/xYh4tquxmJmZ1ZK0JPALYG1g67z7GxExv7iozMxaG7DJCUBEPAFsL2kP4CDSDfq9pM/9DPAnUjer8yOith9uK/sD3wIOJPXTfZbUT/c40ixctYbm7dKkb7EamZjrMjMz66kRwEeAV4CbgVMi4uJiQzIza21AJycVEXEVaUaSrpSdQJ2V3OuUmw0cm1+1JuZXdfk37DMzM+sPEfEioKLjMDPrroE85sTMzMzMzDqIkxMzMzMzMysFJydmZmZmZlYKg2LMiZmZmRVj9oLh3PfyKoXGsNxNcwq9/lKXFD/8Z+bbhhcdAnptaOtC/exXW25UbADnFnt5gNu2mVR0CIxucswtJ2ZmZmZmVgpOTszMzMzMrBScnJiZmZmZWSk4OTEzMzMzs1JwcmJmZmZmZqXg5MTMzMzMzErByYmZmZmZmZWC1znpRBIaUfx85c3E7GLnlO+KPzx+d9EhtLTbmzcpOoSm7t4kig6hpaErFD+vfisaVv5bccwp/79pMzPrfG45MTMz6xCShkh6RFJIekZSub+pMjPrJicnZmZmnWNXYI388xhgzwJjMTPrc05OzMzMOsfH8nZGzXszswHByYmZmVkHkDQa2AsI4EBgAbCbpNUKDczMrA85OTEzM+sMBwNLAFMi4ibgj8BQ4LDagpLG5nEp0xtVlo/XndVC0saSLpf0vKTXJN0p6WOtzjMz6y0nJ2ZmZp2h0oVrct6elbcf7cuLSNoJ+CtpPMtTwBXAy8AkSSf35bXMzGqVf/5KMzOzQU7SJsDGwCvAJXn3FcDzwNskbR8Rf+mD6ywJnAuMBL4NTIyIyMe2Af7Q22uYmTXjlhMzM7Pyq7SaXBQRMwEiYg5wXs3x3toPWBX4D3B8JTHJ17sF+FlXKpF0hKSpkqbOfWlWH4VmZoOBkxMzM7MSk7QE8OH89qyaw5X3+0taug8ut0PeXhgRC+scP78rlUTEpIgYFxHjRiw3qg/CMrPBwsmJmZlZue0NjAb+GxE3Vx+IiL8B9wBLAR/qg2u9OW8fbnC80X4zsz7hMSdmZmblVumytZykm+ocX7mq3C+7UqGkVl9ONpqNq15riplZn3FyYmZmVlKS1gB2yW9XZnEiUs82ktaLiH8Dc/O+Rl291myw//EWx8c2ub6ZWa+VoluXpOl53vTxDY5vKOmJXOYSSSPaHNfYXtYzOdczoW8iMzOzQWIC6Vl9XUSo0Qu4KJevtLI8Q0pQVpS0Up16d29wvRvz9oAGrSsH9exjmJl1TSmSk2YkbQncAKxCGvj3oYiY2/wsMzOzziZJpOQE4JwWxSvHD5E0NCLmAZWphY/PdVXq3Y40TXA9F5PWNnk7cFzNeVsCn+rWhzAz66ZSJyd5Iag/ASsAPwYOj4gFbQxhZ2B9YEYbr2lmZgYwHlgLmAX8pkXZa0itJauyuFXkm6TWk6OAeyVdLOl20hd+dacEjojXgEOAOaQE5l5J50u6DrgZOCMXndfDz2Rm1lRpkxNJewK/I/WXPT4iPl8933o7RMQDETEtfwNlZmbWTpUuWpdFxCvNCkbEfOCC6vPyuiQ7A38G1mBx0nJoRHyjSV3XAtsAV5KSnb1JXxIeDfwoF3u2m5/FzKxLSpmcSDqY9C3RCOBzETGxTpktJZ2cF3l6StJcSY/nMSlbNah3qKQjJd0i6aV8zlOS7pL0g9p+uc3GnCg5QNLvJT2d65oh6c+SPtPks62Tv4V6StIcSdMkfaULM6eYmdkgEhGH5DElH25dGiLis7n8PlX7boqIXSJimYhYKiK2iIjz8rHKeJV6dd0VEXtGxAoRsWREbBIRk0itOQB39vLjmZnVVbrZuiQdBZxGmq7wYxExuUHR75BukvcCt5OaoNcD9gX2lnRQRFxcc84vgcNITeQ3kb75GQOsDXyB1Nf2mS7EOCKX3RNYANwKPAK8CdgA2Ak4tc6pGwOn5OteT5p1ZXvgRGB1oGFSY2Zm1g6SVgZGRcTDNfu3Ak7Obye3Oy4zGxzKlpx8DtiL1Ef2oIi4tEnZ7wMHR8RT1TslfYDU6vJzSVdHxMy8f01SYvIosHmd8zZm8RSKrZxESkz+A+wVEdOq6hkKvL/Bef8DHA98u7LyrqR3kxKVoyWdFBGPdjEGMzOz/rARcK2kfwIPkZ7JawGb5OPnRESrMTBmZj1Stq5Ee+Xtj1skJkTENbUJRt5/JalVYzSwY9WhytzwdzU47+6IeLpVgPkbpaNILTsfrE5Mcj0LIuKKBqffQRo/s7Cq/I3AH0j/LXZscJ6ZmVm7TANOJz2XtiM9m9cErgMOJX3RZ2bWL8rWcnIDsANwjKQ7I+KiZoUljQH2IHWlWp7Fn2eDvF0XuDr/PA14BXi/pK8B59U2WXfRTqSxMDdHxL3dPPd3DQb1TwPeB6zW6ERJRwBHAIwcslQ3L2tmZtY1EfEYafC7mVnblS05mUjqEvVF4DxJNEpQJH0S+CGwZJP6lq38EBGvSPoYcCZpvMp3JM0A/kpKYC6IiNldiLGyau60pqXqe6TB/pfzdmSjE/NAxEkAyw1bqa2zlpmZmZmZtUPZunUREV8ijScZRkpQDqgtI2lzUpPzcOBLpMWilgaG5JlH/q9StKbuS4C3kBa1OhN4FdiPtLjjNElrdCXE7n+qRRa2LmJmZmZmNjiVLjmBugnK/jVF9iUlHj+JiO9HxL8j4rWqLlPrNKn7xYg4OyIOj4i357LXk1pEvteF8CqtH+t14yOZmZmZmVkLZevWtUhEfEkSpC5e5+cuXpWpgUfn7Rtmtsprlezajes8IOk7pMHo7+rCKdeRVsbdRtL6EXFfV69lZmY22Kww/DUOWHVqoTH875X7Fnr9UcMLvTwAC+YPLToEhqw4t+gQWLhOVzrJ9J+XXi2+XWC5IaOKDqGp4n9DTdS0oJwvab98qDLe41BJS1fKS1qG1F1r+dq6JG0i6UOS6v0X+UDethwgn2f0+jnpd/cbSevWXGdons7YzMzMzMy6obQtJxU1LSi/zj+fRVoTZVPgQUk3kbp5vZs0H/uZwMdqqloTuACYKekuUqvLCNK87WuRZvL6ZhfD+hJp4cbdgXsl/RV4jDRd8YZ5W3fVXTMzMzMzq6/ULScVNS0ovwZ2BsaRZq96lTTD1zjgUlLCUm8hw1uBY4EbSaux7w3sAswEfgBsGBFdaneOiDmk1pZDcn0bkAbWvx34O/Cp7n9KMzMzM7PBrRQtJxExtgtlvkRqsaj2yQbFJ+ZX9flPAifmV6/jygspnptfreqZQJohrNHxidTEa2ZmZmY22HREy4mZmZmZmQ18Tk7MzMzMzKwUnJyYmZmZmVkplGLMiZmZmbUm6e2k2Sp3BNYgzQz5DGnGyL8C10TEtYUFaGbWS05OzMzMOoCkDwG/Ik2DPwOYArwArESaqXJrYAfAyYmZdSwnJ2ZmZiUnaRXSGl4jgM8Dp0bEgqrjQ4Dt8svMrGM5OTEzMyu/PYAlgb9GxI9rD+bp7W/MLzOzjuUB8WZmZuW3ct4+3d0TJS0l6cuS7pD0sqRZku6VNFHS0nXKT5QUeftWSedKekrS7HzeMZL85aaZ9QsnJ2ZmZuX3SN7uLGmDrp4kaXXgduB7wJqkQfN/BFYAvgXcLGmFBqe/FZhKGnw/BbgeWAv4PnBx7kpmZtanfGMxMzMrv8uBx4Glgb9Jujq3huwiabl6J0gScBHwDuCnwJoRsVtE7AWsDZwLbAT8qME1DyUlJGtHxIci4n25/Axgb+DIPvt0ZmaZkxMzM7OSi4hXgF1ILRnDgN1JrSHXAs9LujnP5lXtvaQZvG4F/iciZlXVN4uUXDwNHNyg9WQmcHREzK4677/AN/Lbz/fFZzMzq+Y+ox0oFixgwYsvFR1GUxvdpaJDaGm31TYuOoTWhpT8+4PFkwWV1qNnrlp0CC2tts8LRYfQUsyfX3QIg15E3AdsLmkb4P3AlqQphFcAtgG2kfS+iJiQT9k9b3+TB8zX1veapKm53Oak7l7Vro2IemNczgfOANaR9OaImFFbQNIRwBEAK642onsf1MwGtZL/5WNmZmbVIuKWiDguInYBxgDbsjixOEzS/vnntfL25DzA/Q0vFicwK9W51EMNrj8HeCK/Xb1BmUkRMS4ixi29wvBufkIzG8zccmJmZtahcovILZJ2Jw1835Q0HuRiYGgudgMwvUVVD/dTiGZm3eLkxMzMrMNFxAJJ15GSk0oryKN5e3FEnNaDasfW2ylpBFDpL/mGLl1mZr3hbl1mZmYll2feauUteftY3v4+b/evU7Yr3iNpTJ39B5H+fnggIh6rc9zMrMecnJiZmZXf0ZLOkrRF7QFJwyR9Atgv77owby8D7gR2kPRzSaPrnLtKPreeJYHTJC1RVX5t4IT89pSefRQzs8bcrcvMzKz8hgMTgAmSngTuBp4HRpPWHlktlzspIv4AaTyKpL2B3wGfBD4s6R5Sd6+RwLqkNVCeBv5fnWueQ5oV7AFJNwPLkBZkHAlcCfSkq5iZWVNOTszMzMrvl6RB7bsAWwAbAisD80jduM4GzoiIm6pPiojHcmvL4cAB+bwtgedI40V+APy2wTUfJE0x/F1gJ2C5vO9M4Mf1pic2M+stJydmZmYllxdhvCy/unvubFIrR7dbOiLiQeDA7p5nZtZTHnNiZmZmZmal4OTEzMzMzMxKwcmJmZmZmZmVgsecmJmZ2SIRMRGYWHAYZjZIdUzLiaTpkqLmNVvSI5IukrRDH15rSq5/fDfPm5zPm9BXsZiZmZmZDRad2HLyB+DJ/PMKwCak1W/3l/SFiPhRYZGZmZmZmVmPdWJycmJETKm8kTSctErtUcCJki6OiMcKiu1Y4ETgiYKub2ZmZmbWsTqmW1cjETEPOAZ4BRgBvKfAWJ6IiGkR8VJRMZiZmZmZdapObDl5g4iYJek/wGbAmyr7JU0H1gTeGhHTa8+TNAXYAdixujWmpsyOwHG57uHA3cBJEXFFnbKTgcOAj0bE5JpjAj5CWqX3XcAo0uq815Bagx7t8gc2MzPrEMO1gJWGvVxoDPvu8tdCr3/N5G0KvT7ArIeXLjoEWKiiIyDu+nuh1x8+/O2FXh/g6QWvFR1CUx3fclJlubx9qg/r3Af4E7Ay8HvgHmBb4HJJX+hqJTkxORf4FbANcAdplV8BRwN3S9q8D+M2MzMzM+s4AyI5kfRO4K3AXOCPfVj1Z4GvRMRGEfHhiNgW2BOYD5wkaaMu1nMU8GFS4rRJRLwnIg4E3gacCowGLpa0RB/GbmZmZmbWUTo6OZG0gqT3AZeSPsvn+ngw/NSI+H71joi4EjgfGAp8pov1HJO334iIe6vqWgB8EXiE1P1sv0YVSDpC0lRJU+cxpxsfwczMzMysM3RicnJ9ZZ0T4Hngd6Q/7N8XEaf38bXOa7D/nLwd36oCSasDawELq85bJCLmVl2nYX0RMSkixkXEuOG4gcXMzMzMBp5OHBBfWedEwCrAu4GRwK8kbRsR9/fhtR5qsH963q7ehTrenLdPRMTsBmUerClrZmZmZjbodGJyUrvOyaqkhGVD4DxJW0VEdLGudrYcdTUmMzMzM7NBqRO7db1ORDwBHADMA7YADq46PDdvG82ft2aL6se22D+jdYSLyqzWZMD7Wt2oz8zMzMxsQOr45AQgIqYBP8tvJ0qqtAhV/th/w6TSkjYA1mhR9cEt9k/pQmyPkbptDSGtc1Ibx/Du1GdmZlYhaYikD0u6XNIMSXMkPZ8nUDlB0spFx2hm1h0DIjnJvkNaJX5t4JC87895+2VJy1YKSloDmEwat9LM5pI+X71D0u6kJGMB8NMuxvbDvD1B0qJESdJQ4CTgLcDDwCVdrM/MzAa5POHK7aRJVfYgzfx4KXALaXr9rwMPSDqgD641IU9GM7m3dZmZNTNgkpOIeAaoTPv79dx6chrwKLA58G9Jl0q6DrgPeJl0A2/mJ8D3Jd0j6XxJNwFXk8bqHBsRd3cxvJ8BvwZWBe6R9AdJvwb+A3wOeAHYPyI8R7CZmbUkaTTwF2AzUqv7OhGxdUQcFBF7kCaMORZYErhA0r6FBWtm1g0DJjnJfkha6HAt4LCIeIG0ovt5pITi/aQZtk4G3ksap9LMb4HdgOfyuZuQEpp9IuLkrgaVB+gfDBwK3AZsCXyQ9Ps/HXhXRNzR1frMzGzQO400/vEO0lT6r5tdMiLmRcSJwBdIvQTOlDSm7VGamXVTx8zWFRFju1DmVdK3RdX7HqXOWI9sfIN6avf/qWWA6bwJwIQGx4K0zskb1joxMzPrKklrkyaCATi6yTT1kHoAHE6a0fLTwMSaunYDPglsBYwhrR/2AHAF8JOImCVpOosnkDlM0mFVVZydn31mZn1ioLWcmJmZDXR7kJ7f90bE1GYF8xdjv8pv96zsV3I6cA2wD2kCmd8A95AmizkReFMufglwc/75AeDsqtdNffB5zMwW6ZiWEzMzMwPSOBNIg+G7otJt+F2ShkXEfOB/gCNJXaH3johbK4UlCdiRNB6SiPiipAmkbtI3uaXEzPqTkxMzM7POslLePtXF8pVyQ4DRkp4Hjsv7JlQnJrCoteW6XkdpZtYD7tZlZmY2sNVOmz+ONL7ksYi4pl8uKB2R11qZ+tLzC/rjEmY2QDk5MTMz6yzP5u2bmpZarLIQ40LSgPfK4PZ/92VQ1SJiUkSMi4hxy40e2l+XMbMByMmJmZlZZ7kzb7fqYvkt8vaePN4k+j4kM7O+4eTEzMyss1xFagVZX9LmzQrmwe2H5rdX5u0jebte/4RnZtZzTk7MzMw6SETcT5reF+A0SSObFP8ssAHwCmnhRkgtL88Cq+d1Trpibt56Ih0z61dOTszMzDrPp4BHgc2B30kaW31Q0nBJXwF+SOrG9fGIeBrS6vHA/+WiZ0naouZcSdpR0nJVu2fk7fp9/knMzKr4GxAzM7MOExHPStoOuJy0Jsn9km4DHgaWAbYBRgOvAZ+IiItqqvgRKdH4OHCrpKnA/fmcd5AWYnwr8FIufyvwJLBpLnsvMA+4OSLO6rcPamaDjpMTMzOzDhQRj0gaBxwEfIi0OOPmpITkQeB04KcR8WSdcwP4hKTLSYsxbgFsDDxHSlJOJSUjlfJzJL0X+A6wNbAJqffFMMDJiZn1GScnHUhDhjBk1JJFh9HUve8ZVXQILWmJV4sOoaUdpz5fdAhNTdl8xaJDaGm1ff5VdAitDemAqVYXeq2KMoqIBcC5+dWT868iDbDvStl7gD16ch0zs67ymBMzMzMzMysFJydmZmZmZlYKTk7MzMzMzKwUnJyYmZmZmVkpODkxMzMzM7NScHJiZmZmZmal4OTEzMzMzMxKweucmJmZWb959LkxfPGcjxUaw+bv/Weh11/xvrmFXh9g5qrDiw6BN98wv+gQCl+zad6/ly30+gA/XmvrokMALm14xC0nZmZmZmZWCk5OzMzMzMysFJycmJmZmZlZKTg5MTMzMzOzUihFciIpevCaXHTcZmZmZmbWd8oyW9fZdfatAuwGvAZcUuf4Tf0akZmZmZmZtVUpkpOImFC7T9J4UnLybL3jZmZmZmY2sJSiW5eZmZmZmVlHJieSJudxJxMkbSTpYklPSlog6XO1ZRrUMTEfn9hov6TVcz1PSJop6S5J+1WV3VbS7yQ9l49fL2nzOtcam+ucLmmYpK9Kuk/SbElPSTpb0lv69JdkZmYDRn5+VMZc7tGk3D9zmfHti87MrO90ZHJSZVvgdmBTYApwDTCzj+oeC9wJbA/cANwFbAJcJOlASfsA1wNjgGuBh4HxwPWS1m1S74XA8cAjwGXAHOBQ4A5J6/VR7GZmNnB9V1KnP7/NzOoqxZiTXvg48B3gmxGxsI/rPgw4BTgmIhYASDoK+BlwMrAUcHBEXJyPDQHOBz4EfAU4vE6dawKjgE0i4l/5vBHAL4GPAOcAW/Tx5zAzs4FjJrAhcDDpmWFmNqB0+jcv04Bv9UNiAjAd+HIlMckmAc8BqwPXVBITgBzD9/LbHZvUe0IlMcnnzQU+A7wMbC5p23onSTpC0lRJU+fG7J58HjMz63w/ydvj85dbZmYDSqcnJ5fXJA996fqcOCySrzU9v72mzjn/zdvVmtR7bu2OiHgRuDK/HV/vpIiYFBHjImLcCI1sUr2ZmQ1gvyF1Z34rcGRXT5I0XNKnJd0m6WVJs/LYxxMlrVhT9r153MrfmtQ3WtKc/Brd409jZlaj05OTh/ux7sca7H+10fGIqBxbosG5L+ZEpJ7pebt6V4IzM7NB69i8PU7S0q0KSxoJ/BE4FdgAuJH0hdjypG7Id0paq+qUa4HHgY0lbdSg2oOAEcCVEfF8Tz6EmVk9nZ6czOrFua0+e6uuYv3RlczMzKypiLiOlGysDBzThVO+TWqVnwasGxF7RMQBwFqklpg1gfOq6l8A/Cq/ndCgzsPydnL3ojcza67Tk5NmKl2yGn2rtGa7AqmyvKTlGhwbm7cz2hSLmZl1rmOBAI6RtFKjQpJGAUflt5+NiEXPmIiYReoa9iqwVc2Yx8l5e7Ck102eI+kdwObAk9Tv4vy6cZILZr7WrQ9mZoPbQE5OKjfgt9ceyDfr8W2NZrGDa3fkhKUyb/2UtkZjZmYdJyLuAi4ClgGOa1J0M9KXdI9HxLV16nmWOmMeI+LfwK2k1pnda06rtJqcFxHzG8S3aJzk0CWXav2BzMyygZyc/DlvD6lePyQnJqcDRS16+E1J61fFM5w0ZfFywJ0RcVNBcZmZWWf5OjAfOFJSo94Ab87bh5rU82BN2Yqz8raSjCBpKGnqe3CXLjPrBwM2Ocl/5F8FLAvcJekaSVeRbtDvYfFNt50eAW4B7pb0e0kXAA+QbvzPkhZjNDMzayki7gfOIE3C8u1WxXtwiQtJYzv3qJrRaxfSjJR3RsQ/e1CnmVlTAzY5yfYHTgSeBnYirSR/dd4+UkA8ARwAnEAaiLg3MJI0vfDm1eufmJmZdcG3SQszfkTSBnWOV7o4v7VJHZWZul435jEiXgIuI83K9eG8e0LeTu5+qGZmrZU2OYmIKRGhiBhb59iEfGxyizpmR8SxEfHWiBgREatFxOER8WRETMx1TKw5p+7+quPj8/EpDY4rItQkpvkR8b8RsV5EjIyIlSPikIiY3uyzmJmZ1YqIJ0hdg4cA361T5E7SgPc3S9q59mBuEflAfjulzvmLunbl8ZF7kyacOb9XgZuZNVDa5MTMzMy65CTgBVKS8boWkjwj18/z21MkrVo5ltc/OZ00YP7WiLi5Tt1/Bh4lDaz/Nqm132ubmFm/cXJiZmbWwfLivv+X3y5Zp8g3SK0i7wT+K+kKSReSBsLvT+rm/IaZJHPdC4Fz8tvP5u3kvojbzKweJydmZmad71TgsXoHImI2aSKYzwL/AnYE9gJeJrW6bBoRD9Y7N5tc9XPDtU3MzPrCsNZFrLfyeJKG41DMzMyaqTf+sub4bGCNJsfnkRKYU3tw7f/iZ5iZtYlbTszMzMzMrBScnJiZmZmZWSk4OTEzMzMzs1JwcmJmZmZmZqXg5MTMzMzMzErByYmZmZmZmZWCkxMzMzMzMysFRUTRMVg3SXoGeLiPqx0DPNvHdfalsscHjrGvOMa+MRhjXDMiVurD+qwP9NEzq+j/n4u+vmMox/UdQ9/F0PB+7eTEAJA0NSLGFR1HI2WPDxxjX3GMfcMx2kBS9P8rRV/fMZTj+o6hPTG4W5eZmZmZmZWCkxMzMzMzMysFJydWManoAFooe3zgGPuKY+wbjtEGkqL/Xyn6+uAYynB9cAwV/RaDx5yYmZmZmVkpuOXEzMzMzMxKwcnJICRpPUn/I+lcSdMkLZQUkvYrOjYAScMl7SzpB5KmSnpZ0lxJMyRdIml80TECSPqMpIsk3SfpOUnzJD0j6U+SPiJJRcdYS9J383/rkPTFouMBkDS5KqZ6r2lFxwggaZSkL0u6Q9KLkmZKekjSxZK2LTi28S1+h9WvtxQY5+qSTpX0b0mzJM2W9F9JP5e0VlFxWbkU/YwqyzOojM+Ydj9DyvR8KOoZUJb7ezvv38P6sjLrGEcB/1N0EE3sAFybf34SuBF4DXgHsC+wr6QTIuKbBcVX8RVgZeCfwC2kGNcEdgJ2BvaT9MGIWFhciItJ2hz4MhBA6RIn4Gbg/jr7n2h3ILUkvRX4I7AOKZ7rgfmk/957A/eQ4i/Kk8DZTY5vAawPPAA82paIakjaBLgOWB54DPhDPjQO+CRwsKTdIuKWIuKzUin6GVWWZ1CpnjEFP0MKfT4U/Awo/P7e9vt3RPg1yF7Ax4GTgAOAtYEppJvNfkXHluPbCbgE2L7OsQ+RbggB7FhwnNsBS9XZ/07SzSSAjxb9+8wxLQH8C5gB/DbH9sWi48qxTc7xTCg6lgbxLUV6KC4k/bEwtOb4isC6RcfZ4jP8K/+Ov1ZgDLfkGCYBw6v2Dwd+mY/dU/Tvyq/iX0U/o8ryDCrTM6aoZ0gZng9lfwa04/7e7vu3u3UNQhFxRkR8OSIuiogHio6nVkRcFxH7RcRf6hy7kHSzAvhIWwOrERE3RcRrdfbfC5yW3+7a3qga+jbpm5UjgZcKjqXTfJ30B9JpEfG9iFhQfTAinouI/xQTWmuStib9t1/A4n877Y5hJLB1fvutiJhXOZZ//np+u5GkJdsdn5VL0c+osjyDSvaMGczPkNI+A9pxfy/i/u3kxDrR3/J29UKjaG5+3s4pNApA0pbAMcD5EXFl0fF0EkkjgE/ktz8sMpZe+FjeXhMRjxcUwwIW/5to5jVgVj/HYtZbZXgGte0ZM5ifIR3wDGjH/b3t92+PObFO9La8LXwsQj25b+qR+e0VBccyktRX9XnKPc4IYEdJGwFLA08BNwHXRrFjdjYjNdnPiIiHJG0K7EPqB/4U8MeIuKnA+JrK32J9KL/9ZVFxRMQ8SX8GdgOOl/SpyrdvkoYDJ1RijNxXwKzECn0GtfMZU6JnSFHPh9I+A9p1fy/i/u3kxDqKpFWACfntbwoMZRFJHyUNoBxO+iZtG1Kr5Hcj4rdFxgZ8B1gPODAini04llYOrbPvX5IOjIh/tD2aZMO8nSHp+6RvD6t9Q9JlwEfqdb8ogf2BZYCngasKjuVo4BrSt5DvkzQ1798cWAH4MWmwrVlpFfEMKvgZU5ZnSFHPhzI/A9p5f2/r/dvduqxjSBoGnAssB/y5RM3L2wKHAR8G3p33fYPF3yYUQtI2wOeAy3I/6bK6G/gsaSacpYHVgD1Is5+8A/iTpDcXFNvovN2E9FD6MWm2lhWAvUiDQ/cGflZAbF1RafL/VXU/4SJExIOkP6p+T/oDa+/8ejNpQOdfio7RrJkCn0GFPGNK8gy5m2KfD2V+BrTt/t72+3d/jez3q3NelGy2riZxnpHjfARYpeh46sQ3inSzPBmYS7qprlZgLP8BXgBWrTk2mRLN1tXkM4wA/ppj/WlBMXwtXz+Ac+ocH0eawWUhsHbRv7Oa2Napin39EsSzDWmGof8AewJj8msv0kw4AXyz6Dj9Kt+rLM+oop9B7XzGlP0Z0q7nQ1mfAe2+v7f7/u2WE+sIkk4BDif949g5Ip4sOKQ3iIhZEfGviPgScCzwLuCnBYXzXVK/6C9ERCnH5rQSEXOB/8tvdy8ojFeqfv5/tQcjYipwJ2nO/x3aFVQXVb5V+2tE3FdkIJKWBy4jdUF4b0RcERHP5tflwHtJAym/IeltjWsyK0YZnkFtfsaU+hnSxudDWZ8Bbbu/F3H/dnJipSfpB6Rm3WdID4X/FhxSV0zO2w/kAWPttg/pm5zDJE2pfpFuJABH5X1nFBBfV1VW/y2qW9dDDX6uV2aVfo6lyyQNZXEf7cIGwld5P7AScGuk7gGvExH3A7eRxkGOb29oZs2V9Bk0OW/76xnTCc+QdjwfSvcMKOD+3vb7twfEW6lJOgn4AvAcsEtE/KvgkLrqBdLUe8NIfVafKiCGITT/Jmet/Fq+LdH0zIp5+2pB1/9b1c8rUn/13TF5W1SM9exGemC/CpRhvNFb8rbZ+ggv5u3oJmXM2qrEz6B2PGPK/gxpx/OhjM+Adt/f237/dsuJlZakE4EvkW7Cu0bE3wsOqTveTXpovAi0fYaTiBgbEar3Ik0LCfClvG/jdsfXDQfk7R1FXDwiZpC+EQLYufa4pBWATfPbqbXHC3R43l4UEWVImirz729W71vevG+z/LbRt5NmbVXyZ1C/PmM65BnS78+Hkj4D2n1/b/v928mJlZKk/wW+Qrrx7hoRf2t+RntJ2k7SHnn2ltpj27K4qfWXUbOarC0maeP8exxas3+YpGNIXSkAftT+6Bb5Tt5+TdK4ys48///ppJl77iQNziycpDHAB/LbMnTpgjTDy0zSN3A/krRE5UD++SfAGqQ/Av9QSIRmVYp+BvkZU6rnQ2meAQXd39t+/3a3rkEoLyJUPe3dO/L2u5K+WNkZEVu1NbBM0p7Acfnt/cBnJNUrOi0iTmxbYK+3DnAW8KKku0iDJJcB1mbx7/Nq0nSP1thY4LfA8/n3+DSp6XxD0pSRC4EvR0Rhf7BGxJW5z/kxwC2SbiV18dgixzgDOCjylCYlcAhpPYRpEXFL0cEARMTTko4mPUw/BeyT/3tD+sZtVdJK1x+LiGZdB2wQKPoZVZJnkJ8xJXk+lOwZ0Pb7exH3bycng9OywJZ19pdllpzqPovj8queG4CikpMbSHPMb0/6vW1Dmq3jSdLCXOdGxGUFxdZJ7gFOId3k30H6fQbwGOnBfFpE3FlceElEfFHSLcCnSfPdL0maTvSHwIkR8UyR8dX4aN6eWWgUNSLibEn/IK2bsD2waz40g/TQ+2GJ+vNbsYp+RpXhGeRnTImeDyV6BhRyf2/3/Vvl+bLPzMzMzMwGM485MTMzMzOzUnByYmZmZmZmpeDkxMzMzMzMSsHJiZmZmZmZlYKTEzMzMzMzKwUnJ2ZmZmZmVgpOTszMzMzMrBScnJiZmZmZWSk4OTEzMzMzs1JwcmJmZmZmZqXw/wGLzKybWTAbLwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 864x360 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axs = plt.subplots(1, 2, constrained_layout=True, figsize=(12,5))\n",
    "color = \"viridis\"\n",
    "\n",
    "plt.sca(axs[0])\n",
    "plt.yticks(range(8), ['Clinton', 'Kaine', 'O\\'Malley', 'Sanders', 'Cruz', 'Kasich', 'Rubio', 'Trump'], fontsize=22)\n",
    "plt.xticks(range(8), range(1,9), fontsize=22)\n",
    "plt.sca(axs[1])\n",
    "plt.xticks(range(8), range(1,9), fontsize=22)\n",
    "plt.yticks(range(10), ['Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov'], fontsize=22)\n",
    "\n",
    "axs[0].imshow(X_1, cmap=plt.get_cmap(color))\n",
    "axs[1].imshow(X_2, cmap=plt.get_cmap(color))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "def calculate_W(A, B, C, H_A, H_B, H_C):\n",
    "    \n",
    "    alphas = []\n",
    "    n_1, r = A.shape\n",
    "    n_2, r = B.shape\n",
    "    n_3, r = C.shape\n",
    "    weights = np.zeros((r, r, 2))\n",
    "    \n",
    "    \n",
    "    for col in range(r):\n",
    "        for it in range(r**2):\n",
    "            i = int(it % r)\n",
    "            j = int((it/r) % r)\n",
    "\n",
    "            weights[col, i, 0] += np.sum(H_A[i] * H_B[j] * H_C[col]);\n",
    "            weights[col, j, 1] += np.sum(H_A[i] * H_B[j] * H_C[col]);\n",
    "         \n",
    "\n",
    "\n",
    "    return weights"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Display Rank 4 Factor Matrices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeUAAAFwCAYAAACGm2OSAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABEnUlEQVR4nO3dd5hdVdn+8e+dRugQioBBQpGigJTQQUMXRIrSQ4mCCKiIICovrxrkFRFsiAjmh2SQIhBAqoAIhColoShoVEoIhN5LevL8/ljrJIeTc86caafM3J/rmmvP2XvttZ8J7Hlm7b2KIgIzMzNrvH6NDsDMzMwSJ2UzM7Mm4aRsZmbWJJyUzczMmoSTspmZWZNwUjYzM2sSAxodgHXcIA2ORfst0egwqvr4+u83OoR2/efvizU6hF7hPd56PSJWaHQczW75If1j2KoDGxrDP95q7H+m/jMaenkAYqm5jQ6BGU+/VPGecVJuQYv2W4ItF/1co8Oo6uZb7290CO3adeimjQ6hfTGv0RG066/zxj3X6BhawbBVB/LQras2NIY1rv5qQ6+/9KT+Db0+wLyd32p0CDyx12kV7xk/vjYzM2sSTspmZmZNwknZzMysSTgpm5mZNQknZTMzsybhpGxmZtYknJTNzMyahJOymbUkSZMlRTtfe3ey7mH5/MndG7VZdZ48xMxa3a3AyxWOTalnIGZd5aRsZq3ujIgY3+ggzLqDH1+bmZk1iYYkZUm7S7pY0lOS3pc0U9ILkm6SdLSkJYvKjs7vdkaX1DEi7x9f7/jNrPVI2lXS9ZJekTRL0kuS/ihpg3bOGyDpe5L+JWlGPv8iSR+rV+zWd9Q1KUtaMSfRm4BDgFnAX4A/AZOBnYDzgGckrVbP2HJ8ISnqfV0z61mSzgZuAXYDngauBV4CDgQekrR7ldOvAE4lvZ++FpgJHAY8LGmdnova+qK6vVOWtAxwH7AW8Dfg6Ij4e0mZJYFjgFOAZYFqq888BKwHTOuJeM2sd5B0NHAc8CSwb0RMKjq2NzAOuFTSGhFRuoTQasCiwMYR8c98ziDg96SGxcXA5j3+Q1ifUc+W8m9ICfkhYIfShAwQEe9FxJnApsAr1SqLiGkRMSki3LvSrG+7s8JwqDZJ/YEf5HL7FydkgIi4FvgdsAwpyZZzWiEh53NmAd8A3gU2k7RNpcAkHSVpgqQJr73R+HV8rfnVJSlLWhM4KH88OiKqLnUdEU9FxEvt1Fn2nXLx+EIlx0p6TNI0SW9Juk7S+iXnjC5+bF16c5eUlaRDJY3P9c2Q9LSkcyWVXSy1uB5JB0j6W36X/p6k2yVtW+1nNbOqbgUuKvN1L7ARsDLwZHFiLXFX3m5V4fglpTsi4m3ghvxxRKXAImJMRAyPiOErLNf4tYSt+dXr8fUepD8A/hERj9bpmgBtwAHA3cB/gc2APYERkjaOiGdyucdIN/Hh+fNF5SqTJNINejAwGxgPvEl6fHUscKCkz0bEwxXO/xHp0fy9pPfqGwI7ANtKGhERf+vCz2rWV1UcEiVpv/ztJ2voL7JCmX1v5wRczuS8HdpegGa1qldS3jRvyyarHrIasB3wyYh4GkDSIsA1wO7AycBXYP4jrGslHZ4/j6pQ5zGkhPwKsGNEPJnr7Q/8kvRIa5ykdSJiZpnzvwZsHhET83n9gPNzHD8Cdu7aj2xmJQrN06nAX9spO6md42Y9rl5JufAX6Kt1ul7BcYWEDBARMyWdSkrKO3aivhPz9vuFhJzrnSvp28BepD8G9gUuLXP+DwsJOZ83T9L3SUl5O0kDI2J2uQtLOgo4CmCwFu9E6GZ90vN5+1KVP7arWUbS0hHxTpljw/J2amcCMyunN08eMoc0BKJU4a/hVTpSmaShwBrAPFKPyw/JnT8KiXhEhWpuLHPeK8BbwCLAcpWuX/xuapAGdyR0s77sIeANYGNJa3WyjpGlOyQtTXotB+k1llm3qFdSfi1vV6zT9SD9ZTyndGdEvJu/XaSD9X20qN5KHdWeKSlbqlJP8UJMzrZm3Sg/eTqN9Bj7WkkLDV+SNEjSnpLWrVDNDyStV1R+IHA2sDQwMSLu7YHQrY+q1+PricChpI5W9TKvh+rt9OQiEdFTMZlZBRFxdp6M6FvAg5L+TppAZBbpD+iNgcVJE4uUvleeQvr99ZikO4B3gK2BVYHXSZOImHWberWUbyIlyQ0kbVyna3a3wnujVXKHsXLWKClrZk0gIk4APgNcTpqY6HPAZ4HlSa+VRgL3lDsV2J/U2l4D2Jv0ROsSYLMqw6zMOqUuLeWIeErSFaSxyudJ+kyF3snA/HHN09obq9wDZgMDJQ0offQdES9IeoZ0Yx5CmtFnvvxIq/DuaXwdYjXr0yJiWAfL300aHllL2cmAinb9X/4y61H17Oj1ddI71y2AO8pNAi9pcUknkB4XfaSOsRUUWrjrVTj+i7w9rfj9Ux4SdSbwMdLUoFf1WIRmZtZr1W3u64h4M89cdSWwLfB3Sf8kvcMpvNvZnNQB6xXSpBz19ifSe6fb8/uj93PsR+bjvwW2IbX4H8+ziRUmD1mD1It6v2pPAczMzCqpW1IGyI+jt5O0BymxbUV6rzOA1EP7r8B1wGUR8UE9Y8tOIb1D2gf4AjAw7z8SICJC0kjgZtLY4i1Ik9W/SFrd6icR8XxppWZmZrWoa1IuiIgbKTNmt0LZ0cDoMvvH8+F3PoX9k8vtLylT9nhETCdNEHJiueO5TJDGKS80Vrmj1ys6PqzWuszMrPfqzZOHmJmZtRQnZTMzsybhpGxmZtYkGvJO2cysr5kVc3lhzvsNjUGzq3Zv6XGzlmro5YHmT3puKZuZmTUJJ2UzM7Mm4aRsZmbWJJyUzczMmoSTspmZWZNwUjYzM2sSTspmZmZNotmHbFkZEUHMntN+wQYa8cTejQ6hXYMH1nu57k6YO7fREbRvXn0vJ2kysFr++NOI+F6VspewYJ3zuyJiRM9GZ9Y1bimbWSs7NK9nvhBJS5FWezNrGU7KZtaqJgCrADtXOH4gaWnVh+sWkVkXOSmbWatqy9tRFY6PAubSgWVWzRrNSdnMWtWDwL+AvSQtU3xA0jrAVsCtQNnOA5J2knSupMclvSFppqTnJF0kab0K57RJCkmjJG0oaZyklyXNlXR8t/501ic5KZtZKxsLDAYOKtk/quh4JecDRwBzgLuBPwOzgMOACZK2rXLuNsBDwCbAeOAWYFrHQjdbmJOymbWyi0mPqEcVduSOX4cBbwLXVzn328BKEbFpROwTEfsAawNHA4sBYyRVWlbpSOBnwMcj4sCI+FxEjOnyT2N9npOymbWsiHiZ1ErdvOiR8y6kDmCXRcSsKudeGxFvl+yLiPgdcD+wHvCJCqdPAn4YEXUeEGa9nccpm1mrawM+R2otf5cFrea29k6UNDSfuy6wFFAYXrVS3q4NPFnm1Osiot1B5JKOAo4C+OhH3Qay9jkpm1mrux54gzRm+SxgL+AfETGx2kmSTgX+h+q/B5eqsP+5WgLLj7THAGy44cCo5Rzr2/ynm5m1tPyI+jJgZVLHrkWo3sELSV8EfgBMB74CrAksFhGKCAF/LBStUMX0bgjdbCFOymbWG7Tl7R6k3tSXtlN+v7z9n4i4ICKeiYjiRLtWN8dnVhMnZTNreRHxCHAf6TH2uIh4tZ1ThuTt86UHcoexjbs3QrPatFRSljQ5D9wfUeH4BpJeymWukjSok9cZleto60K4ZlZHEbFtRCwfEQfXUHxS3n6l+PeEpBWBi3B/G2uQlkrK1UjaAriL1GtyLHBAteEQZtan/Qp4h9Tz+qk8M9eNwNPAEsC1jQvN+rJekZQl7QD8FViWdLMdUctwhSr+RBqjeHLXozOzZhMRz5AeUV9O6sz1edI9P4Y0Pec7jYvO+rKWf0QjaU/gSlKPy1MjYnRX64yId/BNadaUImJYB8tfRZle1BHxLAtPz1kwijILXURE2f1m3aWlW8qSRgJXA4OA48slZElbSDpL0gRJr0iaJenF/M55ywr1ln2nLGlE3j9e0kBJp0iaJGmGpFclXSLpY1XiXVXS2ZL+LWm6pHcl3ZevV2nohZmZ9REtm5QlHUOa91bAlyPi7ApFfwx8CxhImkC+MNHAF4F7Je1X4bxqBgI3A98DnsrfzwNG5jqXKRPv9sA/gONI/+63kFa52ZD0DvyiTsRhZma9SKs+vj6eNGvPLOCgiLimStmfASMj4pXinZI+T2plny/ppojoyAovW5MWWF+zMPRC0tLAHaRVY75G+mOgcK2V87WWID36+kNERD62KukPhUMl3RERbR2Iw8zMepFWbSnvlbe/aichExG3lCbkvP8GYBxpvOL2Hbx+kFrn88dC5vfQP80fdywpfzypE9rPI+KiQkLO5z1PmlEI4BuVLijpqPwIfsLsmNHBcM3MrBW0akv5LuAzwImSJkbEldUKS1qeNNPP+sAyLPi518/btYGbOnD9KRHxjzL7C2MfVynZv3vejqtQ30TgfWAjSYMjFs66xXPoLtVvOc+ha2bWC7VqUh5NGl/4beBSSVRKzJK+CvyCtD5qJZUmna9kSoX97+bt4JL9a+TtwzX051oOmNrBeMzMrBdo1aRMRJyUE1zFxCxpM+A80ly4JwE3AC8A0yIiJJ1OGovc0Z7PHV1DtbAc3BVAe8+eZ3awbjMz6yVaNilD2cQcEVH8iPiLpIT764j4WZkq6jXp/PP5WqdFRLm1Wc2sl5sRA/jnrGUbGsO8xTranuheQ3/yQEOvD/Ds5Rs2OoSqWrWj13wRcRKph/UA4LKSIU7VJp1fAdi55yME0pApWLAyjZmZ2UJaPilD2cS8bz5U6Hh1mKQlCuUlLQlcSOr0VQ9nkd43/4+kr0la6AmFpE9K+kKd4jEzsybUK5IyLJSY/5gT81hSK3kT4BlJ10j6EzAZGE5KzPWI7Xlgb+A94DfAFEm3SbpU0k2SpgBPAPvXIx4zM2tOvSYpw8KJmTReeDhpKNH7pB7bw4FrSIl6ocfaPRjbncAngdOBV4EtSe+8Pwk8Q+pwdkq94jEzs+bTUh29apmIPifmk0p2f7VC8dH5q7SONqCtzP7xVOmpHRGT2zn+MinxOvmamdlCelVL2czMrJU5KZuZmTUJJ2UzM7Mm4aRsZr2CpH6SpuQ1z1+TNLDRMZl1lJOymfUWOwOr5u+XB/ZsYCxmneKkbGa9xZfzdmrJZ7OW4aRsZi1P0hDSOusBHAjMBXaVVLqMqllTc1I2s95gJLAIMD4i7gX+Qlqd7fByhfN758jfHyXpUUnTJL2RZ/5bv4bzjpD0oKR38/5leuIHs77FSdnMeoPCo+q2vB2bt1+qdpKkX5KWd30HuA54HdgHeFDStlXOO4c0U+BM4EZgIqmVbtYlLTWjl5lZKUkbAxuR5pa/Ku++HngT+Lik7SLingqnHwVsHxF357pEmgr3e6TFbdaOiHJroB8KbBURD3XfT2LmlrKZtb5CK/nKiJgGEBEzgUtLjpdzXiEh5/MC+F/SfPSrkuanL+dMJ2TrCW4pt6IIYu7cRkdR1av3NH//mlVnPdfoENon/91cjaRFgIPzx7Elh8cC3wD2k/SNiHi/TBWXlO6IiLmS/kiao34EC5J7sWtqjO8oUmucFVfxr1trn+94M2tlewNDgP9GxH3FByLiUeBxYHHggArnP1th/+S8HVrheE1/0UXEmIgYHhHDlx7ipGzt8/8lZtbKCo+ml5Z0b5njKxaV+313XTQipndXXWbFnJTNrCVJWhXYKX9ckQUJuJytJa0TEf8u2T+M1JouNSxvp5Y5ZtZj/PjazFrVKNLvsDsiQpW+gCtz+XIdvkaW7pDUnzQBCcD4HojbrCInZTNrOXno0qj88eJ2iheOH5oTbrFji8cj53pPBdYktZKv7nq0ZrXz42sza0UjgDWA6bSfOG8BXgNWBnYHbig69v+AuyTdDbwEbAKsk+sd6XfHVm9uKZtZKyo8ir42It6rVjAi5gCXl5xXcAJp2NQQUk/uFYFrgS0i4q7uCtasVm4pm1nLiYhDSbNq1Vr+OOC4MvsD+G3+qqUe1XpNs85wS9nMzKxJOCmbmZk1CSdlMzOzJtGhpCypn6SDJV0naaqkmZLelDRB0mmSyg7elzQ6rzc6qmT/qML6pJKmV1uPVNLHi8qGpBEdib1KvfPXRy3ZP747r2NmzaNoDLNZU6k5KUsaCjxEmpx9D2AKaVL2+4HVSSurPC1p/07GMpgFA/bLGdXJes3MzFpCTUlZ0hDgHmBT0gw3a0XEVhFxUETsAawEnAwsBlwuqdJyZ5U8CsylQuKV1A84DHgLeKqDdZuZmbWEWlvK55Lmgn0Y2C0iPrSySkTMjogzSGP+BFwoafkOxPEicBuwhaR1yxzfibRay+XAzA7Ua2Zm1jLaHacsaU2g8Ej62IiYUaX4r4EjgA2ArwOjOxDLWOCzpNby90qOfamozKcrxLkCaV3V3Ugz8qxMSuD/Ik2zd35EdNsixJJ2Bb4GbAEsC7xBeopwekT8o6jcdsDdwKSIWK9CXcsDLwABDI2IN7orTjNrDtPmDeLR6cMaGsMyK7/b0Ov3X2v1hl4fQE8v3ugQqqqlpbxHLvdkREyoVjAPxP9D/rhnB2O5jvR4+kPz0+bOX3vn6z9c5fxdgV8BnyStkfon4BFgI+A3wNV5Xtsuk3Q2aeq+3YCnSTMAvUR6J/6QpN0LZSPiHtIqNOtK2qFClUcCiwCXOyGbmfVdtSTlTfP2oRrrLCTOT0kaABARo3Nvx7ZKJ0XETOCPwCrALkWHDiR1Aqt4bjYR2DIiVo2IHfL77u1JndAeA/ZiQYu/0yQdTZoZ6Elgg4jYOiL2j4hNgH1ITx8ulbRs0Wnn5O2xZerrB3w1fzy3q/GZmVnrqiUpr5C3r9RYZ6FcP9J8sh0xNm9HFe37EjCHdlaCiYh/RcSDZfa/BHwnf9y3g/F8SG7B/yB/3D8iJpVc61rgd8AywCFFhy4D3gT2krRKSbWfI7+vb+9JhJmZ9W49MXlIpx8R56T0BCl5LStpPWBz4JaIaPePAkkDJO0i6YeSzpM0VlIbcHQusnZnY8s2Ir2rfjIi/lmhTGES+60KO/JKMxeQWtFHlZQvtJ6rtpIlHZXHg0+Y7b5uZma9Ui0LUryetx+psc7CBCLzSK3DjmoDfgYcRGpBwoIWdEWS1ia92y3bmSpbqhPxFFsjbz9ZbsKREiuUfD4XOBH4iqT/i4g5uRPdrqROYldUqywixgBjAJbSkPaubWZmLaiWpDyR9Ch2yxrr3DxvH89LpnXUJcAZpCXWViH9UXBD1TOSq0gJ+XrgTFKv63ciYm5O2P+mC634rNABbSrw13bKlj7aniLpetJ7571zvMfkmC5sp1e7mZn1AbUk5RuBnwPrSdqsWg/o3Lv5sPyxlkS6kIh4RdLNwOfzrl9HxOxq5+SxzRsArwJfKDP0aa3OxFLG83n7UkSM6sT555CS8rGSbiS9L58HnNc94ZmZWStr951yRDxFatUBnCtpcJXixwHrA+/RtZ7EF5Ae6b5BDY+uWdCh7MUKY5FHdiGWYg/lmDaW1OFEHxF3kt6Zbw+cSor75tLJWMzMrG+qtaPX10itxM2AP0saVnxQ0kBJ3wV+QZoA48iIeLWzQUXE9RGxfP56rIZT/ktqca4v6UOTi0j6Eun9dJflFvtppMfY10ravLSMpEGS9qwwMxmkMdOwoEd4TYurm/VVkiYXLURzRjtlLykqO75OIZp1m5qSckS8DmxLGu+7PfCUpPskXSbpBuBl0nvg6cDIiLiyh+KtFN9rpOQ2ALhT0h05tn8AF+bYuutaZwO/JE1S8qCkxyVdI+lySfeQOrddx4JOaqUuJk2SAvAMaRISM6vNhyYXKiZpKeALdY7HrFvVPCQqIqYAw4FDgZtJk3LsS0rWk4Efkxaq+GP3h1mTb5KGGz1O6my2G2nM9G7kXsvdJSJOAD5Dmot7WdJY488Cy5PewY8kLeBR7txppJW1AM6LiHndGZtZLzaB1Plz5wrHDwQWZcEERmYtp5aOXvPl97WX5K8uyzN8tXXwnPUr7J8H/L/8VU7ZnteV1lSNiBHtxHE3aU7rDslrTu9MeqpwYUfPN+vD2kgNg1GUf8I0irTa3MWkV21mLacnJg+x6k4BBgEXRURnxnGb9VUPkoY67pXnxJ9P0jqkCXtuJc1DvxBJO0k6N79yekPSTEnPSbooT1RUWv72/G664jrvkn6ey5zZlR/MrMBJuQ4kbS3pAkl3knqov0PqMGZmHTOWNBd+aefNUUXHKzmftIrdHNJTrj8Ds0jDOCdI2rakfMU56wEkLYqHNVo3c1Kuj7VJvww2B+4FPhsRLzY2JLOWdDHpEfWowo7c8eswUifL66uc+21gpYjYNCL2iYh9SPfm0cBiwJiSleRuAJ4DtpNU7rXZQaQ+JR7WaN3GSbkOIqItr5K1eERsFxEPNDoms1YUES+T3idvXvTIeRdSB7DLImJWlXOvjYi3S/ZFRPyO1PlyPeATRcfmsmDIYrnWcmFfxWGNxXPWf/BWxdDM5nNSNrNW05a3o0q2bbRD0lBJX5X0S0m/l9SWF61ZKRcpXbTmAlKnzEMkLVlUzxakZW2rDmuMiDERMTwihi++7KD2wjPrWO9rM7MmcD1pZr1DJZ1FWiv9HxExsdpJkk4F/ofqv/c+tGhNRLwp6VLgSNJw0NKWs4c1WrdyS9nMWkp+RH0ZaRnVscAitDMdr6QvktZCnw58BVgTWCy/VhJQmF+h3BDJQoevY3JdywH7AzPwsEbrZk7KZtaK2vJ2D1Jv6kvbKb9f3v5PRFwQEc/kdc4LKs5lHxF/J/XWLkzjewSpB/jlHtZo3c1J2cxaTkQ8AtxHeow9roa59guL1jxfeiB3GNu4nfMLreWvk3prQ9cW3TEry++UzawlRUTpuOJqJpFm0vuKpJsLvbTzDHsX0f7vwmtJCb3Q4n44IiZ0LGKz9rmlbGZ9wa9Ik/Z8jrSgzri8pvnTwBKkpFtRRMzhwxOEuJVsPcJJ2cx6vYh4hvSI+nJSZ67Pk8YljyFNz/lODdXclrdvAFf0QJhmfnzdkiQ0sLn/033toBsaHUK7rv/xio0OoX0ebUNEDOtg+aso04s6z7pVaW31URTNElbByLz9fUTM6EhMZrVyS9nMrB2SViUNpZoF/KbB4Vgv1tzNLTOzBpJ0BjCU1ElsceCsiFioB7dZd3FSNjOr7EDgY6TlIH8KfL+x4Vhv56RsZlZBR99nm3WV3ymbmZk1CSdlMzOzJuGkbGZm1iT8TtnMrA4G9ZvDaoNeb2gMs+f2b+j153xk6YZeH2DWsnMbHUJVbimbmZk1CSdlMzOzJuGkbGZm1iS6nJQlrSvpfEn/ljRN0nRJUyTdL+nnknbujkB7gqTRkkLS6EbHYmZm1qWOXpIOAP4ADAKmAuOBt4AVgE1Iq698hgWrq5iZmVkFnU7KklYCLiQl5G8B50TE3KLj/YBt85eZmZm1oyst5T2AxYC/RcSvSg9GxDzg7vxlZmZm7ejKO+XCYrSvduQkSVtIOkvSBEmvSJol6UVJV0nassI589/9SvqIpN9JekHSTEnPSjpD0uAK5w6U9G1J/5Q0Q9LLki6WtFqNsV6erzVL0muSrpdUtvWfY4z8/RGSHpT0bt6/TN6/iqTfSHoqxzMtv4O/RdJRNf9DmtlCJI3P99uoRsdi1hldScpT8nZHSet34Lwfkx53DwQeAq4H3gC+CNwrab8q564KTCS10v9Geoe9IvBd4MrSwvkR+jXAWcDqwB3AXcCOuZ7VK11I0on5GvsDLwPXAU8BnwPukvSVKueeA4wBZgI35muFpJXz918jPaW4BbiB9G+5JXBClZ/dzMx6ua48vr4OeBFYBXhU0l9ICe8R4OGIeKfCeT8DRkbEK8U7JX0euBo4X9JNETGtzLlfBi4AvhYRs/J565GS++clbRMR9xWV/xopgU8FRkTEU/mcwcAlwGHlApS0W47zReALEfFg0bFtgD8D50q6KyL+U6aKQ4GtIuKhknq/CawE/A44JiKi6NgiwBbl4jEzs76h0y3liHgP2AmYQEruu5PWG70NeFPSfbl3dul5t5Qm5Lz/BmAcMATYvsJlnweOKyTkfN6/gIvzxx1Lyh+ft/9bSMj5nBnAscD0CtcZnbdHFifkfO59wGmklv5XK5x/ZmlCzj6St7cUJ+Rc78yI8Pt3M7M+rEvjlCPiXxGxGbANcDpwO2lIVD9ga+BySW2l50laXtIoST+TdIGktlyu8Bh87QqXvCMiyiXSSXm7StE1hgJrAPOAy8rE/irwl3KxAZsD75Y7nt2Vt1tVOH5Nhf2FRP1TSXtLWrxCuYVIOiq/h58wO2bUeppZn5d/v1R8z1xpvoKu9GUx66xuWZAiIu4H7of573G3BH4I7AIcnh9Hj8vHvwr8gtRzu5KlKuyfUmH/u3lbfIMMzdsXi1vWJSaX2Vd4z7wUMEdSpRghjccu57kK+y8m/ZscDPwJmCvpCVIP9cvzv2NZETGG9J6apfotF5XKmVm3K/RlEen33FKkoZ7fBT4B7Nm40Ky36fZVovJQqPsl7U5qGW4C7A2Mk7QZcB4wBziJ1MnpBWBaRISk04GTSf/zlzOvu+Mto7CMyjvAte2ULbvkS4XWfOHfZqSkn5DedW+Tv74BfEPShRFxRGeCNrMe09G+LGad1mNLN0bEXEl3kJJyoUX5RVLC/XVE/KzMaWt1YwhT83YVSYMqtJaHldn3fN7OjohR3RjPfBHxBPAEzH+ysDvpEfuXJV0REZUem5tZ/ZXtyyLpYuAYUl8WJ2XrFp1+p6x2nutmH8vbF/J2SN4+X1pQ0gpAt82THRHPA8+SfsYDa71eREwF/gEsL2lEd8VTSUTMi4gbSb3ZAT7V09c0sw6puS9LqeK+IO+9OadnorNepSsdvY6VNFbS5qUHJA3I43j3zbuuyNvC/8SHSVqiqPySpCk7l+lCPOX8Om//T9IaRddbBDiXyu+1v5+3l0japfSgpP6SdlCFyU4qkXSYpE3K7F+OBZ3GKr2PNrPG6Ehflg+JiDERMTwihi85pMceTFov0pX/SwYCo4BRkl4GHgPeJLWGN2TBX49nRsSt+fuxpGFKmwDPSLqX9Dj708AsUmL+chdiKnUOqWPVbsCT+XH6+6ROGoNJi2ksNFY5Iq7Lk4ecCdwq6T/Av/O5KwEbk/6AOAZ4oAPxfAG4SNJU0r/X28BywHbA4sA9pA5gZlY/7TVO6tGXxQzoWkv598A+pBbn88AGwH6kBPs+cBGwXUR8t3BCRLwFDCf1In6fNDvWcNIQok0o81i7K/ICGXsB3yP1tN6JNAb67nzdZ6uc+wtgU9LP2Z/0qPvzpF7ddwNfocwsYu34OXA2aVKS4aR/rw1JE64cAewcEbM7WKeZVVd4F7xEhePtTrlrVi+dbinnyUOupf0eyqXnvUrlSTdGs2DijuJzyu4vOt4GtFU4Nps0qclPa71e0bmPA0dWOl6mfNX37BFxD6k1bGb1U+j0uW7pAUmLAiPqGo1ZFV2aPMTMrAXcnreHSlqnsDMn5PNY0CHVrOGclM2sN5r/Hjgi7iUtDLMU8Eheke1G0uurXUh9XcyagpOymfUmi+btByX79wPOIC01uwOpD8tNeVupd7VZ3bmPvpn1CnnuhMLQx8nFx/IiNCfnr1Kj6ea+LGad5ZaymfUWhwPLA68Bjzc4FrNOcUvZzFqWpMVI65OvyYIJeL4fEZ4+y1qSk7KZtbJBwCHAe6T5p88urEhn1oqclM2sZUXE21ReVc6s5fidspmZWZNwUjYzM2sSTspmZmZNwu+UW9DSn5jNZ8e90ugwqvrz5xZaobLp9BvU3P+GAPNmzGh0CNZN5kR/3phbaU2M+lh56XfbL9SDBr7Q0MsD0G/GRxsdQlVuKZuZmTUJJ2UzM7Mm4aRsZmbWJJyUzczMmoSTspmZWZNwUjYzM2sSTspmZmZNwknZzMysSTgpm1mfIWm8pJA0qhvqGpXraut6ZGaJk7KZmVmTcFI2MzNrEk7KZmZmTaJPJGVJu0u6WNJTkt6XNFPSC5JuknS0pCUbHaOZNYaktmrvmSWNzsdH1zcy64t69SpRklYErgQ+k3f9C/gLMAsYCuwE7A6cJml4RDzXkEDNzMzoxUlZ0jLAfcBawN+AoyPi7yVllgSOAU4BlgWclM3MrGF6bVIGfkNKyA8BO0TEQgvTRsR7wJmSrgE+qHN8ZmZmH9Ir3ylLWhM4KH88ulxCLhYRT0XES/nc+e+XJG0oaZyklyXNlXR8aZkK11/oHZSkEXlfe1/Duv4vYGbNQNJRkiZImvD+W7MaHY61gN7aUt6D9AfHPyLi0U7WsQ1wPjAVGA8sCUzrQkwvAxdVOLYSsGv+fm4XrmFmTSQixgBjAD62/lLR4HCsBfTWpLxp3j7chTqOBH4M/CAi5nU1oIiYBIwq3S9pCeDu/PEXEfF8V69lZt2qVz5RtObUW/9nWyFvX+1CHZOAH3ZHQq5EUn/gCmBj4GrgpJ66lplVVHiuvESF46vVKxCz3pqUu8N1EdHTj5J/QxqS9QBwSLU/AIrfTX3w5uweDsusT5mat+uWHpC0KDCirtFYn9Zbk/JrebtiF+ro0eFRkr4DHA08DexZQ2e0MRExPCKGLz5kYE+GZtbX3J63h0pap7AzJ+TzgI81JCrrk3prUp6Yt5t1oY7pXTi36r+rpP2AM4A3gN0j4rVq5c2s281/KhUR9wI3AksBj0i6RdKNwLPALsDYxoRofVFvTco3kW66DSRt3AP1d/odlKStgT/kOvaKiP90c2xmVtmieVs6L0HhD+VXgR2ATUi/RzYBptQtOuvzemVSjoinSB2oAM6TtEi18pLWlLRyBy7RqXdQktYCrgMWAQ6LiPs6cE0z6wJJAtbIHycXH4uIGRFxckSsHhGDImKViDgiIl6OiNERoYgYXXJOW94/qi4/gPUJvTIpZ18HngG2AO6QtEFpAUmLSzqB9Lj7Ix2ou8PvoCQtB9wMLA98LyKu7MD1zKzrDifdf68Bjzc4FrOyeus4ZSLiTUnbkhak2Bb4u6R/koY6zQI+CmxOarW+ArzZgbrvze+c9iC9g7oHmAMMJz02Hwt8qeS0Y0nTfr4PfEJSW4Xqvx0Rr9cai5lVJmkx4HfAmsBWeff3I2JO46Iyq6zXJmWAPHXmdpL2IE27uRXwWdLP/RrwV9Lj5MsioqNzX+8H/BA4kPQO6nXSO6hTSL2qS/XP2yVIf7FXMjrXZWZdNwg4BHiPtEDN2RExrrEhmVXWq5NyQUTcSOpdWUvZUZSZeatMuRnAyfmr1Oj8VVx+oX1m1rMi4m1AjY7DrFa9+Z2ymZlZS3FSNjMzaxJOymZmZk2iT7xTNjNrtDdmLs5Fz27Z0BgG9W/syrAvHbRSQ68PMPC9RkdQnVvKZmZmTcJJ2czMrEk4KZuZmTUJJ2UzM7Mm4aRsZmbWJJyUzczMmoSTspmZWZPwOOUW9M6TA7j5k8s0OoyqBqza2PGQtYi58xodQrsGju/IMt8N8plGB2DWe7ilbGYtR1I/SVMkhaTXJA1sdExm3cFJ2cxa0c7Aqvn75YE9GxiLWbdxUjazVvTlvJ1a8tmspTkpm1lLkTQE2AsI4EBgLrCrpFUaGphZN3BSNrNWMxJYBBgfEfcCfwH6A4eXFpQ0LL93nlypsnw8KhzbSNJ1kt6U9IGkiZK+3N55Zp3lpGxmrabwqLotb8fm7Ze68yKSdgD+Rnpf/QpwPfAuMEbSWd15LbMCD4kys5YhaWNgI+A94Kq8+3rgTeDjkraLiHu64TqLAZcAg4EfAaMjIvKxrYFbu3oNs3LcUjazVlJoJV8ZEdMAImImcGnJ8a7aF1gZ+A9waiEh5+vdD/y2lkokHSVpgqQJc96d1k2hWW/mpGxmLUHSIsDB+ePYksOFz/tJWqIbLleYEuWKiCg3y8xltVQSEWMiYnhEDB+w1GLdEJb1dk7KZtYq9gaGAP+NiPuKD0TEo8DjwOLAAd1wrY/m7XMVjlfab9YlfqdsZq2i8Gh6aUn3ljm+YlG539dSoaT2GiaVelc3/xyt1pKclM2s6UlaFdgpf1yRBQm4nK0lrRMR/wZm5X2VHmmvVmH/i+0cH1bl+mad1hSPryVNzmP+RlQ4voGkl3KZqyQNqnNcw7pYT1uuZ1T3RGbW54wi/b66IyJU6Qu4MpcvtKpfIyXm5SStUKbe3Stc7+683b9Ca/qgzv0YZtU1RVKuRtIWwF3ASqTOHAdExKzqZ5lZbyFJpKQMcHE7xQvHD5XUPyJmA4UhUqfmugr1bksa7lTOONLY5HWBU0rO2wL4Wod+CLMaNXVSzoP3/wosC/wKOCIi6rkm4I7AeiyYX9fM6m8EsAYwHbi6nbK3kFrHK7OgFfwDUmv5GOBJSeMkPUT6Y7/s0KaI+AA4FJhJStxPSrpM0h3AfcAFuejsTv5MZmU1bVKWtCfwZ9K7oFMj4lvFYwXrISKejohJ+a9tM2uMwqPoayPivWoFI2IOcHnxeXlc8Y7A7aSVpQrJ+rCI+H6Vum4DtgZuICX5vUkNhGOBX+Zir3fwZzGrqimTsqSRpL+IBwHHR8ToMmW2kHRWHpj/iqRZkl7M75y3rFBvf0lHS7pf0jv5nFckPSLp56XvnKq9U1ayv6SbJb2a65oq6XZJ36jys62V/+J+RdJMSZMkfbeGXqBmfVJEHJrfGR/cfmmIiONy+X2K9t0bETtFxJIRsXhEbB4Rl+ZjhffR5ep6JCL2jIhlI2KxiNg4IsaQWu8AE7v445l9SNP1vpZ0DHAuacjBlyOirULRH5NujCeBh0iPmdYBvgjsLemgiBhXcs7vSZPWTwfuJf2VuzywJnAC6T3SazXEOCiX3ZO0Qs0DwBTgI8D6wA7AOWVO3Qg4O1/3TlIP0u2AM4ChQMVkbmb1JWlFYNGIeK5k/5ZAYe7rtnrHZb1bsyXl40lLss0CDoqIa6qU/RkwMiJeKd4p6fOkVvb5km4qTMUnaTVSQn4e2KzMeRuxYBhEe84kJeT/AHtFxKSievoDn6tw3jeBU4EfFWYJkvRpUoI+VtKZEfF8jTGYWc/aELhN0hPAs6TfS2sAG+fjF0dEe++4zTqk2R6Z7pW3v2onIRMRt5Qm1rz/BlIrdgiwfdGhwrjGRyqc91hEvNpegPmv52NILfkvFCfkXM/ciLi+wukPk96Pzysqfzdpcvt+JfGaWWNNAs4j3Zvbkn4/rQbcARxGmaUizbqq2VrKd5HmnD1R0sSIuLJaYUnLA3uQHhkvw4KfZ/28XRu4KX8/ibSyzOck/Q9waeljqRrtQHrXfV9EPNnBc/9cobPaJGA3oOIi7ZKOAo4CGIzn0DXraRHxAqlTl1ndNFtSHk169Ptt4FJJVErMkr4K/AKqZqilCt9ExHt5cfILSe+jfyxpKmm91JuAyyNiRg0xFmb4mVS1VHlTKux/N28HVzoxdy4ZA7CUhnhhdTOzXqjZHl8TESeR3hcPICXm/UvLSNqM9FhpIHASaYD/EkC/3IvyJ4WiJXVfBXyMNBHBhcD7pCXaxgKT8lR+7YbY8Z9qPs+Xa2ZmFTVdUoayiXm/kiJfJCXcX0fEzyLi3xHxQdGj4bWq1P12RFwUEUdExLq57J2kFvBPawiv0NpdpwM/kpmZWbua7fH1fBFxUp7Z7tvAZflRdmGI05C8Xainch5rvHMHrvO0pB+TOll9qoZT7iDN4rO1pPUi4l+1XsvM+q5FB8xmg+VeamgMU6ct3dDrz3m+8Q8LX9+k0RFU15Qt5YKSFvNlkvbNhwrvcw8rXtBc0pKkx9LLlNYlaWNJB0hatMylPp+37Xb8yj20zyf9210tae2S6/TPw7LMzMw6pGlbygUlLeY/5u/HksY0bwI8k9dWFfBp0ljCC1kwNV/BaqTp96ZJeoTUyh5EGnO4Bqln9g9qDOsk0oQju5PmxP0b8AJp2NUGeVt2hiAzM7NKmrqlXFDSYv4jaR7b4aTeyO+TemwPB64hJepyE3A8AJxMWpJtKGke252AacDPgQ0iYkKN8cwkta4PzfWtT+owti7wd7yCjJmZdUJTtJQjYlgNZU4itVCLfbVC8dH5q/j8l0nTWZ7RHXHlCUAuyV/t1TOKBUvPlTs+mpJ4zcys72mJlrKZmVlf4KRsZmbWJJyUzczMmkRTvFM2M+sISeuSRmBsD6xKGu3wGmkUxN+AWyLitoYFaNZJTspm1lIkHQD8gTSkcSowHngLWIE0+mIr0sI2TsrWcpyUzaxlSFqJNA/BIOBbwDkRMbfoeGGZxW0bE6FZ1zgpm1kr2YO0MtzfIuJXpQfzUMW785dZy3FHLzNrJSvm7asdPVHS4pK+I+lhSe9Kmi7pSUmji6frLSo/WlLk7eqSLpH0iqQZ+bwTJblhY93KSdnMWklhlbYdJa1f60mShgIPkVaCW43UGewvwLLAD4H7JC1b4fTVgQmkTmXjSavKrUGaZXBcfmRu1i38P5OZtZLrgBdJ66c/Kumm3PrdSVLZJZCUJsy/EvgE8BtgtYjYNSL2Is1hfwmwIfDLCtc8jJSI14yIAyJit1x+Kmm63qO77aezPs9J2cxaRkS8R5qzfgKpT8zupNbvbcCbku7LvbOLfZbUI/sB4JsRMb2ovumkpPoqMLJCa3kacGxEzCg677/A9/PHb3XHz2YG7ujVuvr1b3QEVcX06e0XarCp3xze6BDatcr2f2t0CE0nr2G+maStSYvRbEEaCrUssDVprfPd8pzzkBI3wNW5I1hpfR9ImpDLbUZ6rF3strxka6nLgAuAtSR9NCKmlhaQdBRwFMDiKy3esR/U+iS3lM2sJUXE/RFxSkTsBCwPbMOChHq4pP3y92vk7Vm549ZCXyxI3CuUudSzFa4/E3gpfxxaocyYiBgeEcMHLzO4gz+h9UVuKZtZy8st4Psl7U7q0LUJ6X3vOKDwWOkuYHI7VT3XQyGa1cRJ2cx6jYiYK+kOUlIutHoL66uPi4hzO1HtsHI7JQ0CVs4fF3p0bdYZfnxtZi0j96Ruz8fy9oW8vTlv9ytTtha7SFq+zP6DSL9Dn46IF8ocN+swJ2UzayXHShorafPSA5IGSPoKsG/edUXeXgtMBD4j6XxJQ8qcu1I+t5zFgHMlLVJUfk3gtPzx7M79KGYL8+NrM2slA4FRwChJLwOPAW8CQ0hjh1fJ5c6MiFshvW+WtDfwZ+CrwMGSHic91h4MrE0aw/wq8P/KXPNiUi/vpyXdByxJmkhkMHAD0JlH4mZlOSmbWSv5Pamz1k7A5sAGpKk3Z5MeV18EXBAR9xafFBEv5Nb1EcD++bwtgDdI74N/DvypwjWfIQ2VOh3YAVg677sQ+FW5YVZmneWkbGYtI08ecm3+6ui5M0it2g63bCPiGeDAjp5n1lF+p2xmZtYknJTNzMyahJOymZlZk/A7ZTOzMiJiNDC6wWFYH9MyLWVJk8vMWTtD0hRJV0r6TDdea3yuf0QHz2vL543qrljMzKzvaMWW8q3Ay/n7ZYGNSTP17CfphIiotCaqmZlZU2vFpHxGRIwvfJA0kDSjzjHAGZLGNXDKu5OBM1iwcoyZmVnNWubxdSURMRs4EXgPGATs0sBYXoqISRHxTqNiMDOz1tWKLeWFRMR0Sf8BNgU+UtgvaTKwGrB6REwuPU/SeOAzwPbFre+SMtsDp+S6B5Km9TszIq4vU7YNOBz4UkS0lRwTcAhpRqFPAYuSZhK6hdT6fx4z67WmzRnIY6+t0n7BHrT4oNkNvf67wxrfDhz0VqMjqK7x/0LdZ+m8faUb69wH+CtpGr+bgcdJC6lfJ+mEWivJCfkS4A/A1sDDpBmJBBwLPCZps26M28zMWlCvSMqSPgmsDswC/tKNVR8HfDciNoyIgyNiG2BPYA5wpqQNa6znGOBg0h8MG0fELhFxIPBx4BzSZPrjilehMTOzvqelk7KkZSXtBlxD+lmO7+ZOXhMi4mfFOyLiBuAyoD/wjRrrOTFvvx8RTxbVNRf4NjCF9Jh93zLnAiDpKEkTJE2YzcwO/AhmZtYqWjEp31kYp0xasu3PpIS2W0Sc183XurTC/ovzdkR7FUgaCqwBzCs6b76ImFV0nYr1RcSYiBgeEcMH4ga1mVlv1IodvQrjlAWsBHyatK7pHyRtExFPdeO1nq2wf3LeDq2hjo/m7Ut5lZpynikpa2ZmfVArJuXSccorkxL1BsClkraMiKixrno+Kag1JjMz66Na8fH1h0TES6RFy2eTFj0fWXR4Vt4uUeH01dqpflg7+6e2H+H8MqtU6ci1RgfqMzOzXqrlkzJAREwCfps/jpZUeAJQSHLrlp4jaX1g1XaqHtnO/vE1xPYC6fF0P9I45dI4BnakPjNbmKR+kg6WdJ2kqZJmSnozd448TdKKjY7RrBa9IilnPybN6rUmcGjed3vefkfSUoWCklYF2kjvpavZTNK3indI2p2UXOcCv6kxtl/k7WmS5v+BIKk/cCbwMeA54Koa6zOzLHemfIjUYXIP0miGa4D7SUMl/xd4WtL+3XCtUbmjaVtX6zIrp9ck5Yh4DSgMX/rf3Fo+F3ge2Az4t6RrJN0B/At4l3TTVvNr4GeSHpd0maR7gZtI7+JPjojHagzvt8AfgZWBxyXdKumPwH+A44G3gP0iwmOdzDpA0hDgHtKMe+OBtSJiq4g4KCL2IHUGPRlYDLhc0hcbFqxZDXpNUs5+QZqgYw3g8Ih4izQD16WkRPo5Uo/ps4DPkt5DV/MnYFfgjXzuxqREvk9EnFVrULnj2UjgMOBBYAvgC6R///OAT0XEw7XWZ2bznUvq4/EwaVjkh0ZMRMTsiDgDOIH0ZOxCScvXPUqzGrVM7+uIGFZDmfdJfxkX73ueMu9ysxEV6ind/9d2A0znjQJGVTgWpHHKC41VNrOOk7QmqZMnwLFVhhxCeup1BGmUxteB0SV17Qp8FdgSWJ40B8LTwPXAr/P8+pNZ0Dn0cEmHF1VxUb7/zbqkt7WUzazv2IP0O+zJiJhQrWD+o/gP+eOehf1KziMtDLMPqXPo1aR57lclLcVaWOTmKuC+/P3TwEVFX/d2w89j1jotZTOzEpvm7UM1li+8IvqUpAERMQf4JnA06bXX3hHxQKFwXkhme1KfDyLi25JGkV6J3euWsfUEJ2Uza1Ur5G2tK8MVyvUDhkh6k7QsK8Co4oQM81vXd3Q5SrMO8ONrM+srSodADie9P34hIm7pkQsWLSQz551pPXEJ62WclM2sVb2etx+pWmqBwgQi80gduQqdtv7dnUEVK15IZsDSi/XUZawXcVI2s1Y1MW+3rLH85nn7eH6f7Pnorek4KZtZq7qR1OpdT9Jm1QrmTluH5Y835O2UvF2nZ8Iz6zgnZTNrSXmZ1sLUtOdKGlyl+HHA+qSpeM/N+yaSHoEPzeOUa1FY5MadZK1HOCmbWSv7Ggum0v2zpGHFByUNlPRd0mx/ARwZEa9Cmu0L+EkuOlbS5iXnStL2kpYu2l1Y5Ga9bv9JzPBfe2bWwiLidUnbAteRxhQ/JelB0gIvSwJbA0OAD4CvRMSVJVX8kpRgjwQekDQBeCqf8wnSBCKrA+/k8g8ALwOb5LJPkqbrvS8ixvbYD2p9hpOymbW0iJgiaThwEHAAaVKRzUiJ+BnS/PK/iYiXy5wbwFckXUeaRGRzYCPSfPdPAeeQknCh/ExJnyWtSrcVaT78fqTfpU7K1mVOyi1IgwYxYOjQRodRVf+25l/w6qO7TGy/UKMNGtToCNpXbcbpOomIucAl+asz599I6jhWS9nHSVN8mnU7v1M2MzNrEk7KZmZmTcJJ2czMrEk4KZuZmTUJJ2UzM7Mm4aRsZmbWJJyUzczMmoTHKZuZ1cPbA9C1yzU0hOe3mt3Q6y8xq/0yPW2FRxs/h0K1tULdUjYzM2sSTspmZmZNwknZzMysSTgpm5mZNYmmSMqSohNfbY2O28zMrDs1S+/ri8rsWwnYlbT82lVljt/boxGZmZnVWVMk5YgYVbpP0ghSUn693HEzM7PepikeX5uZmVmLJmVJbfm98ihJG0oaJ+llSXMlHV9apkIdo/Px0ZX2Sxqa63lJ0jRJj0jat6jsNpL+LOmNfPxOSZuVudawXOdkSQMkfU/SvyTNkPSKpIskfaxb/5HMepl8/xT6lOxRpdwTucyI+kVn1j1aMikX2QZ4CNgEGA/cAkzrprqHAROB7YC7gEeAjYErJR0oaR/gTmB54DbgOWAEcKektavUewVwKjAFuBaYCRwGPCxpnW6K3ay3O11Sq//+MltIU7xT7oIjgR8DP4iIed1c9+HA2cCJETEXQNIxwG+Bs4DFgZERMS4f6wdcBhwAfBc4okydqwGLAhtHxD/zeYOA3wOHABcDm3fzz2HW20wDNgBGku4Zs16j1f/SnAT8sAcSMsBk4DuFhJyNAd4AhgK3FBIyQI7hp/nj9lXqPa2QkPN5s4BvAO8Cm0naptxJko6SNEHShFnzuuthgFlL+nXenpr/qDXrNVo9KV9XkjS70505Yc6XrzU5f7ylzDn/zdtVqtR7SemOiHgbuCF/HFHupIgYExHDI2L4oH6LVanerNe7mvTaanXg6FpPkjRQ0tclPSjpXUnTc9+OMyQtV1L2s/m99KNV6hsiaWb+GtLpn8asSKsn5ed6sO4XKux/v9LxiCgcW6TCuW/nBFzO5LwdWktwZn3cyXl7iqQl2issaTDwF+AcYH3gbtIfwsuQXjdNlLRG0Sm3AS8CG0nasEK1BwGDgBsi4s3O/BBmpVo9KU/vwrnt/eztPRLviUfmZlaDiLiDlGRXBE6s4ZQfkZ5CTQLWjog9ImJ/YA1Sy3s14NKi+ucCf8gfR1Wo8/C8betY9GaVtXpSrqbw6LnSX9Gr1SuQIstIWrrCsWF5O7VOsZi1upOBAE6UtEKlQpIWBY7JH4+LiPn3WERMJz0Cfx/YsqRPR1vejpT0oU6xkj4BbAa8TPlXWYVy8/uCzJnxQc0/mPVdvTkpF268dUsP5Jt0RF2jWWBk6Y6cqAvjLsfXNRqzFhURjwBXAksCp1Qpuinpj/MXI+K2MvW8Tpk+HRHxb+ABUmt895LTCq3kSyNiTpUY5/cFGTB48eo/kBm9OynfnreHFo//zQn5PKBRk3X8QNJ6RfEMJA29WhqYGBGe09usdv8LzAGOllTp6ddH8/bZKvU8U1K2YGzeFpIwkvqThjCCH11bN+u1STkntxuBpYBHJN0i6UbSjbkLC262epoC3A88JulmSZcDT5Nu+NdJk4iYWY0i4ingAlLnyh+1V7wTl7iC1Hdlj6Ie2juRRlhMjIgnOlGnWUW9Niln+wFnAK8CO5Bm/ropb6c0IJ4A9gdOI3Uw2RsYTBomtVnx+GUzq9mPSBOKHCJp/TLHC6+yVq9SR6Hn9Yf6dETEO6SZ9wYBB+fdo/K2reOhmlXXtEk5IsZHhCJiWJljo/KxtnbqmBERJ0fE6hExKCJWiYgjIuLliBid6xhdck7Z/UXHR+Tj4yscV0SoSkxzIuL/ImKdiBgcEStGxKERMbnaz2Jm5UXES6RXQP2A08sUmUjqyPVRSTuWHswt4M/nj+PLnD//EXbu/7E3qSPpZV0K3KyMpk3KZmYdcCbwFim5fqhFnHtYn58/ni1p5cKxPH75PFJHsAci4r4ydd8OPE/qMPYj0tMtj022HuGkbGYtL0/K85P8sdyUd98ntYI/CfxX0vWSriB18NqP9DproZERue55LJhj+7i8beuOuM1KOSmbWW9xDhVm4ouIGaQOnscB/yTNT78Xac75M4FNIuKZcudmbUXfVx2bbNYVrb5KVEvI74srvmc2s/aV619ScnwGsGqV47NJifucTlz7v/getjpwS9nMzKxJOCmbmZk1CSdlMzOzJuGkbGZm1iSclM3MzJqEk7KZmVmTcFI2MzNrEorozMIp1kiSXgOe6+ZqlyetVNWsmj0+6LsxrhYRK3Rznb1ON9y3zfD/l2PonutXvGeclA0ASRMiYnij46ik2eMDx2g9qxn+2zmGnr++H1+bmZk1CSdlMzOzJuGkbAVjGh1AO5o9PnCM1rOa4b+dY+jh6/udspmZWZNwS9nMzKxJOCn3QZLWkfRNSZdImiRpnqSQtG+jYwOQNFDSjpJ+LmmCpHclzZI0VdJVkkY0OkYASd+QdKWkf0l6Q9JsSa9J+qukQyQ13VJ/kk7P/61D0rcbHY9V1uj7tBnuw2a9x3ryPvJ6yn3TMcA3Gx1EFZ8BbsvfvwzcDXwAfAL4IvBFSadFxA8aFF/Bd4EVgSeA+0kxrgbsAOwI7CvpCxExr3EhLiBpM+A7QOC1gVtBo+/TZrgPm+4e6+n7yC3lvukJ4CzgAGAt4K7GhrOQecDVwKcjYuWI2CMiDoiIDYADgbnA9yVt39AoUyzLRsQmEfH5iDgwIrYCNgBeAfYCDm9ohJmkRYCLSHFd1+BwrDaNvk+b4T5sqnusHveRk3IfFBEXRMR3IuLKiHi60fGUiog7ImLfiLinzLErgLb88ZC6BlYiIu6NiA/K7H8SODd/3Lm+UVX0I2A94GjgnQbHYjVo9H3aDPdhE95jPX4fOSlbK3o0b4c2NIrq5uTtzIZGAUjaAjgRuCwibmh0PNZrNPo+rOs9Vq/7yEnZWtHH8/alhkZRgaTVSX9JA1zf4FgGkx63vUlz9yOw1tOw+7De91g97yN39LKWImklYFT+eHUDQ5lP0pdInWIGkloNW5P+4D09Iv7UyNiAHwPrAAdGRKMXErBeot73YRPcY3W7j5yUrWVIGgBcAiwN3N5Ej2K34cOdTeYA3wd+0ZhwEklbA8cD1+Z3gGZd1qD7sGH3WL3vIz++tlZyPmkYxPM0uJNXsYg4MiIELAZ8EvgVMBp4QNIqjYhJ0qKkjjjvAsc2Igbrtep+HzbqHmvEfeSkbC1B0tnAEaTxkjtGxMsNDmkhETE9Iv4ZEScBJwOfAn7ToHBOJ73zOyEimvLdu7WeRt+HDbjH6n4fee5rQ9J40vua/SLiqgaHsxBJPwdOAF4DRkTEPxscUrskLUdaCH0OsFhEzK7z9ScDqwILDWcB1gU+AjxDau08FRFH1i8664xG36fNdh/W4x5rxH3kd8rW1CSdSfpF8AawU6N/EXTAW6RfFgOAIaTJBuqtH+mXeCVr5K9l6hKNtawmvQ/rdY/V9T7y42trWpLOAE4i3Xw7R8TfGxxSR3ya9MvibdJf83UVEcMiQuW+SEM7AE7K+zaqd3zWOpr4Puzxe6wR95GTsjUlSf9Hmvf2bdIvgkern1FfkraVtEfuiVp6bBvg9/nj7yNibn2jM+sejbwP++o95sfXfZCkTYDfFu36RN6eXrziSURsWdfAMkl7Aqfkj08B36iwGMykiDijboF92FrAWOBtSY+QOr4sCazJgn/Pm0jDNsw6rNH3aRPch33yHnNS7puWArYos//jZfY1wpCi74fnr3LuAhqVlO8CTgO2I/27bU1aMeZl0mQKl0TEtQ2KzXqHRt+njb4P++Q95t7XZmZmTcLvlM3MzJqEk7KZmVmTcFI2MzNrEk7KZmZmTcJJ2czMrEk4KZuZmTUJJ2UzM7Mm4aRsZmbWJJyUzczMmoSTspmZWZP4/3bMHNaDINlWAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 504x360 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "w = calculate_W(A_A1, B_A1, C_A1, A_S1, B_S1, C_S1)\n",
    "\n",
    "fig, axs = plt.subplots(1, 2, constrained_layout=True, figsize=(7,5))\n",
    "color = \"viridis\"\n",
    "\n",
    "A_A1_new = A_A1 @ w[:,:,0].T \n",
    "B_A1_new = B_A1 @ w[:,:,1].T\n",
    "\n",
    "plt.sca(axs[0])\n",
    "plt.xticks(range(4), range(1,5), fontsize=22)\n",
    "plt.yticks(range(8), ['Clinton', 'Kaine', 'O\\'Malley', 'Sanders', 'Cruz', 'Kasich', 'Rubio', 'Trump'], fontsize=22)\n",
    "plt.sca(axs[1])\n",
    "plt.xticks(range(4), range(1,5), fontsize=22)\n",
    "plt.yticks(range(10), ['Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov'], fontsize=22)\n",
    "\n",
    "axs[0].imshow(A_A1_new, cmap=plt.get_cmap(color))\n",
    "axs[1].imshow(B_A1_new, cmap=plt.get_cmap(color))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Display Rank 2 Factor Matrices"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAU8AAAFwCAYAAAA44fmKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAA/1klEQVR4nO3deZhcRdn+8e9NFvbVgIABwiKLAgqEHSQssoksKvs2iiCgooK+irxiEPVFcENFMD+UQRaBALIKiEBAQIEERUWjAoZ9B9lCEpLcvz+qOmk6fXp6enq6p2eez3X1dabPqVOnOul5pk5VnSrZJoQQQu8s1O4ChBBCJ4rgGUIIDYjgGUIIDYjgGUIIDYjgGUIIDYjgGUIIDRje7gKE3hu13DCPWWVEu4tR07/+uli7i9Cj1/zyC7aXb3c5Bqphiy3uEcss19YyrL/88229/rTH3+KFl+ao2rEInh1ozCojuPemVdpdjJp2XXVsu4vQo5vfuuTRdpdhIBuxzHKMOfL4tpbh3mN+2tbrb7bL44XH4rY9hBAaEMEzhBAaEMEzhBAaEMEzhBAaEMEzhBAaEMEzhBAaEMEzhBAaEMEzhBaRNE2Se3jt3WDeY/L505pb6lAkBsmH0Ho3Ac8UHHuslQUJjYvgGULrnWZ7UrsLEfombttDCKEBbQmeknaXdIGkhyS9LmmmpCckXS/paElLlqUdn9tyxlfkMS7vn9Tq8ofQCpJ2kXSNpGclzZL0tKRfSdqgh/OGS/qKpH9ImpHPP1/Sqq0q+1DQ0uApaYUc7K4HDgFmAb8Ffg1MA3YCzgYekbRaK8uWy2dJsSJeaDtJZwI3ArsBDwNXAU8DBwD3Stq9xumXAqeQ2k+vAmYChwH3SVqn/0o9tLSszVPSMsBdwFrAH4Cjbf+lIs2SwDHAScCyQK1Zb+4F1gOm90d5Q2gXSUcDxwEPAh+zPbXs2N7AROAiSWvYfrni9NWARYGNbP89nzMS+DmpwnIBsFm/f4ghoJU1z5+QAue9wA6VgRPA9mu2Twc2AZ6tlZnt6ban2o7eydBpbisYptQtaRhwck63X3ngBLB9FfAzYBlSMKzm1FLgzOfMAj4LvApsKmnrooJJOkrSZEmT50x/o+EPOBS0JHhKWhM4ML892vaMWultP2T76R7yrNrmWT7eTcmxkv4sabqklyVdLWn9inPGl9+uV36pK9JK0qGSJuX8Zkh6WNJZkqpOslmej6T9Jf0ht/W+JukWSdvU+qxh0LkJOL/K607g/cBKwIPlAbDC7Xm7ZcHxCyt32P4vcG1+O66oYLYn2B5re+ywxRav+SGGulbdtu9BCtR/tf2nFl0ToBvYH7gD+DewKbAnME7SRrYfyen+TPryHp7fn18tM0kifTEPAt4CJgEvkW6DjgUOkLSr7fsKzv8GqUniTlK774bADsA2ksbZ/kMfPmvoHIVDlSTtm398bx3t79Vmwf9vDpTVTMvb0T0VMPSsVcFzk7ytGlT6yWrAtsB7bT8MIGlh4Epgd+BE4EiYdyt0laTD8/uugjyPIQXOZ4EdbT+Y8x0G/IB0azRR0jq2Z1Y5/9PAZran5PMWAs7J5fgG8MG+feQwCAzL2yeB3/WQdmoPx0M/alXwLP2FfK5F1ys5rhQ4AWzPlHQKKXju2EB+J+Tt10qBM+c7R9IXgb1IQftjwEVVzv96KXDm8+ZK+hopeG4raYTtt6pdWNJRwFEAq74rnm0YxErrPjxd4494LctIWtr2K1WOjcnbJxspWHi7wTxIfjZpqEel0l/rlXuTmaTRwBrAXFKP5dvkRvlSwBxXkM11Vc57FngZWBh4R9H1y9uiln/HsKJkofPdC7wIbCRprQbzOLhyh6SlSc1nkJqbQh+1KniWlsBboUXXg/SXe3blTtuv5h8X7mV+7yrLt6jD65GKtJWKRgaUyrRIL8sUBpl853Eq6fb9KkkLDCuSNFLSnpLWLcjmZEnrlaUfAZwJLA1MsX1nPxR9yGnV/d8U4FBSh02rzO2nfBseRG+7v8oUBhHbZ+aHRL4A3CPpL6SB8rNIf5g3AhYnDaCvbPd8jPT79mdJtwKvAFsBqwAvkAbLhyZoVc3zelIw20DSRi26ZrOV2olWzh1P1axRkTaEhtg+HtgOuIT0wMiHgF2BUaTmn4OB31c7FdiPVHtdA9ibdEdzIbBpjeFPoZdaUvO0/ZCkS0ljPc+WtF1BbzQwb1zo9J7GevaDt4ARkoZX3vLbfkLSI6Qv5CGkJzbmybdGpbamSS0oa+gwtsf0Mv0dpGF29aSdBqhs1zfzK/STVnYYfYbUJrg5cGu1yQ0kLS7peNJtxztbWLaSUo1xvYLj38/bU8vbm/JQpdOBVUmPlF7ebyUMIQwILRvzYvul/CTNZcA2wF8k/Z3UZlNqy9mM1JHzLGnweav9mtTOdEtuL3o9l/2T+fhPga1JNegH8tNNpUHya5B6zfetVasOIQwOLR0wmG/Dt5W0BykAbUlqxxlO6pH/HXA1cLHtdjxYexKpzWgf4CPAiLz/kwC2Lelg4AbS2MzNSZMwPEWaDer/bD9emWkIYfBpy2hr29dRZcxjQdrxwPgq+yfx9jae0v5p1fZXpKl63PabpIHwJ1Q7ntOYNM5zgbGevb1e2fEx9eYVQhgYBvMg+RBC6DcRPEMIoQERPEMIoQExw0QIoTrBnIXbuyrNE7Nfb+v1Z3lO4bGoeYYQQgMieIYQQgMieIYQQgMieIYQQgMieIYQQgMieIYQQgMieIYQQgNinGcH+teMpdl16ofaXYyaPOepdhdhwJA0jbQwIMB3bH+lRtoLmT8v7O22x/Vv6UKjouYZQmsdmud/XYCkpUizeYUOEMEzhNaZTFq19YMFxw8gTXF4X8tKFBoWwTOE1unO266C413AHHox3WFonwieIbTOPcA/gL0kLVN+QNI6pMnBbwKqrt0laSdJZ0l6QNKLkmZKelTS+eVLDVec0y3JkrokbShpoqRnJM2R9PmmfrohJoJnCK11Hmk1ywMr9neVHS9yDnAEMJu0MNxvSEvYHAZMzsvcFNkauBfYmLRA4Y3A9N4VPZSL4BlCa11AujXvKu3IHUiHkdbDuqbGuV8EVrS9ie19bO8DrA0cDSwGTJBUtGrBJ4HvAu+2fYDtD9me0OdPM4RF8AyhhWw/Q6r1bVZ2q70zqSPpYtuzapx7le3/Vuyz7Z8Bd5NWfX1PwelTga/bntvHjxCyGOcZQut1Ax8i1T6/zPxaaHdPJ0oanc9dF1gKKA17WjFv1wYerHLq1XaNySnn538UcBTA8KWX7Sn5kBbBM4TWuwZ4kTTm8wxgL+CvtqfUOknSKcBXqf17u1TB/kfrKVi+lZ8AsMi7VmnvTMgDXNy2h9Bi+db8YmAlUgfRwtTuKELSR4GTgTdJy16vCSxmW3l11l+VkhZk8WYTih7KRPAMoT2683YPUu/5RT2k3zdvv2r7XNuP5KWyS9ZqcvlCDyJ4htAGtu8H7iLdvk+0/VwPpyyXt49XHsgdTxs1t4ShJx0VPCVNywN+xxUc30DS0znN5ZJGNnidrpxHdx+KG0JNtrexPcr2QXUkn5q3R5Z/ryWtAJxP9F+0XEcFz1okbQ7cTup1PA/Yv9awjxA6zA+BV0g97Q/lJ4WuAx4GlgCual/RhqZBETwl7QD8DliW9CU7op5hGTX8mjRm7sS+ly6EvrP9COnW/BJSp9CHSd/RCaTHOl9pX+mGpo6v6kvaE7iM1GN5iu3xfc3T9ivElzE0ie0xvUx/OVV6zW3/hwUf6yzposqEI7ar7g9919E1T0kHA1cAI4HPVwuckjaXdIakyZKelTRL0lO5TXSLgnyrtnlKGpf3T5I0QtJJkqZKmiHpOUkXSlq1RnlXkXSmpH9KelPSq5LuytcrGmISQhiAOjZ4SjqG9JywgE/YPrMg6beALwAjSBMjlAYofxS4U9K+BefVMgK4AfgK8FD+eS5pBvA7K2fMyeXdHvgrcBzp3/1G0iw7G5LaaM9voBwhhDbp1Nv2z5OeypgFHGj7yhppvwscbPvZ8p2SPkyqtZ4j6XrbvZlhZivSxLZrloaYSFoauJU0a82nSUG7dK2V8rWWIN1C/dK287FVSAH9UEm32u7uRTlCCG3SqTXPvfL2hz0ETmzfWBk48/5rgYmk8XPb9/L6JtV2543Ny+2k38lvd6xI/3lSZ9b3bJ9fCpz5vMdJT4wAfLbogpKOyk0Pk996JR4WCaHdOrXmeTuwHXCCpCm2L6uVWNIo0pMc6wPLMP9zr5+3awPX9+L6j9n+a5X9pbF4K1fs3z1vJxbkNwV4HXi/pEVsz6hMUP7M8VLrvDOeOQ6hzTo1eI4njXf7InCRJIoCqKRPAd8nzXdYpGgyhSKPFex/NW8Xqdi/Rt7eV0e/0DuAJ3tZnhBCi3Vq8MT2l3IgKgygkjYFziY9O/wl4FrgCWC6bUv6NmksZ297uns7J2Jp2rBLgQVqlRVm9jLvEEIbdGzwhKoB1LbLb40/SgqMP7L93SpZtGoyhcfztU61XW2uxRAGHM2FEa+3dwTdTW+0d76TV+e+XnisUzuM5rH9JVKP+nDg4oqhR7UmU1ie4iVgm+2GvG1kWFQIYQDq+OAJVQPox/KhUgfOYZKWKKWXtCTwC1LnUSucQWoP/aqkT0taoMYv6b2SPtKi8oQQ+mhQBE9YIID+KgfQ80i1zo2BRyRdKenXwDRgLCmAtqJsjwN7A68BPwEek3SzpIskXS/pMeBvwH6tKE8Ioe8GTfCEBQMoabzlWNIQn9dJPfRjgStJAXWB2/l+LNttwHuBbwPPAVuQ2mTfCzxC6rg6qVXlCSH0TUd1GNUzwUIOoF+q2P2pguTj86syj26qLMZlexI1euZtT+vh+DOkABlBMoQON6hqniGE0CoRPEMIoQERPEMIoQERPENoE0kLSXoszxH7vKQR7S5TqF8EzxDa54PAKvnnUcCebSxL6KUIniG0zyfy9smK96EDRPAMoQ0kLUeal9bAAcAcYBdJldMZhgEqgmcI7XEwadHCSbbvBH5Lmn3r8GqJc7toafWBoyT9SdJ0SS/mJ+fWr+O8IyTdk9fOcrXlYkL9IniG0B6lW/TuvD0vbz9e6yRJPyBNs/gKcDXwArAPcI+kbWqc92PSk3YzgetIE3DHpNp90FFPGIUwGEjaCHg/aa6Dy/Pua4CXgHdL2tb27wtOPwrY3vYdOS+RHvn9CmlSnLWrrUQAHApsafve5n2SoS1qniG0XqnWeVlp4UHbM4GLKo5Xc3YpcObzDPwvaX6EVUjzJVRzegTO5oqaZwea9dpIHp20WruLUdOqjpVEqpG0MHBQfntexeHzSIsA7ivps7arzcR7YeUO23Mk/Yo0Z8I45gfhcjUXSiwr31Gk2i3Dl1q2nlOGrKh5htBae5Mm6f637bvKD9j+E/AAsDiwf8H5/ynYPy1vRxccf7SewtmeYHus7bHDF1u8nlOGrKh5htBapVvypSXdWeX4CmXpft6si9qO9aqbLIJnCC0iaRVgp/x2BeYHymq2krSO7X9W7B9Dqp1WGpO30V7SInHbHkLrdJF+5261raIXUFoFtlrH0cGVOyQNIw20B5jUD+UOVUTwDKEF8pCirvz2gh6Sl44fmgNjuWPLx3PmfE8B1iTVOq/oe2lDPeK2PYTWGAesAbxJzwHuRuB5YCVgd+DasmP/D7hd0h3A06TlZNbJ+R4cbZutEzXPEFqjdAt+le3XaiW0PRu4pOK8kuNJw5mWI/XcrwBcBWxu+/ZmFTb0LGqeIbSA7UNJT/nUm/444Lgq+w38NL/qyadwTa3QN1HzDCGEBkTwDCGEBkTwDCGEBvQqeOY1Vw6SdLWkJyXNlPSSpMmSTpVUddCvpPF5/sCuiv1dpfkGJb1Za35BSe8uS2tJ43pT9hr5zpvvsGL/pGZeJ4S+KBsDGgaIuoOnpNHAvaRJB/YAHiNNNnA3sDppZpeHJe3XYFkWYf5A32q6Gsw3hBCarq7gmZcM+D2wCekJhrVsb2n7QNt7ACsCJwKLAZdIKpoWq8ifSMsQdBVcfyHgMOBl4KFe5h1CCE1Xb83zLNKzs/cBu9l+28wutt+yfRppDJqAX0ga1YtyPAXcDGwuad0qx3cizRZzCWkm7BBCaKsex3lKWhMo3YofWzBLdcmPgCOADYDPAON7UZbzgF1Jtc+vVBz7eFmaDxSUc3nSPIm7kZ64WIkUaP9BetztHNtzelGemiTtAnwa2BxYFniRVCv/tu2/lqXbFrgDmGp7vYK8RgFPkJZFGG37xWaVM4RGeSGYvVh7V+o4d9rWbb3+87OKb3TrqXnukdM9aHtyrYR5AO8v89verkF9Nem2/G3P8+ZOpL3z9e+rcf4uwA+B95LmPPw1cD9puYOfAFfk54D7TNKZpEfodgMeJj3h8TSpzfZeSbuX0ublFB4A1pW0Q0GWnyQtBnZJBM4QOkM9wXOTvK13Cv9SgHufpOEAtsfn3sLuopPyMgS/AlYGdi47dACpM6nw3GwKsIXtVWzvkNtjtyd1Zv2ZtMxro51Z80g6mvTkx4PABra3sr2f7Y1JC3ENBy6SVD4N94/z9tgq+S0EfCq/Pauv5QshtEY9wXP5vH22zjxL6RYiPX/bG6VlCbrK9n0cmE0PM9HY/ofte6rsfxr4n/z2Y70sz9vkGvHJ+e1+tqdWXOsq4GfAMsAhZYcuJi3utVeVdbk/RG5P7qlmH0IYOPpjkHzDt8Y5ePyNFGSWlbQesBlwo+0eg7ek4ZJ2lvR1SWdLOk9SN3B0TrJ2o2XL3k9qS33Q9t8L0pQmZ9iytCPPdHMuqVZ6VEX6Um20Zq0zr9U9WdLkOdPf6G25QwhNVs/EIC/k7TvrzLM0UH4uqbbVW93Ad4EDmT87duVCWQuQtDap7bFqp0y2VAPlKbdG3r632sD6CstXvD8LOAE4UtI3bc/OnXG7kDqbLq2Vme0JpHW3WWTlVWK97RDarJ7gOYV0C7pFnXlulrcP5Km1eutC4DTSVFwrk4L3tTXPSC4nBc5rgNNJveyv5JUF1wb+SR9qxVmpI+tJ4Hc9pK28pX9M0jWkdtG9c3mPyWX6RQ+jGEIIA0w9wfM64HvAepI2rdXjnXuzD8tv6wl4C7D9rKQbgA/nXT+y/Vatc/LY0A2A54CPVBmStFYjZani8bx92nZXA+f/mBQ8j5V0Hak9dy5wdnOKF0JolR7bPG0/RKolAZwlaZEayY8D1gdeo289x+eSbmVfpI5bduZ3TD1VMJZzgXVfGnRvLtNGknodkG3fRmrT3Z60dMJywA2VDx2EEAa+ejuMPk2qdW0K/EbSmPKDkkZI+jLwfdJA70/afq7RQtm+xvao/PpzHaf8m1SDW1/S2wbRS/o4qf20z3IN+FTS7ftVkjarTCNppKQ9C56UgjTmFOaPAKhrUtvQOSRNK5vA5rQe0l5YlnZSi4oYmqCu4Gn7BWAb0njJ7YGHJN0l6WJJ1wLPkNopS+uoXFaYWT+w/TwpCA0HbpN0ay7bX4Ff5LI161pnAj8gDca/R9IDkq6UdImk35M6ya5mfmdXpQtIDwMAPEIabB8Gr2qLuAEgaSngIy0uT2iSuocq2X4MGEtaSuAG0uDzj5GC6jTgW6QJQ37V/GLW5XOkYUAPkDqtdiONOd2N3EvdLLaPB7YjPWu/LGms5q7AKFIb8cGkiVSqnTudNBMVwNm25zazbGFAmUzq9PxgwfEDgEWZ/2BJ6CC9WsMotydemF99lp846u7lOesX7J9LWlnw/xWcWrWnvWiORNvjeijHHaRn1nslz3n6QVIt/Re9PT90lG5ShaOL6ncYXaTZxC4gNYmFDhIzybfeScBI4HzbjYyDDZ3jHtKQub0qJ/qWtA7pQYqbSPMiLEDSTpLOyk1DL+bJxx+VdH5+gKQy/S257bRwXlxJ38tpTu/LBwsRPFtC0laSzpV0G2lEwiukjqcw+J1HmpuhstOyq+x4kXNIs5TNJt3l/AaYRRoOOFnSNhXpC+dQAJC0KDE8rmkieLbG2qRfgs2AO4FdbT/V3iKFFrmAiom+cwfSYaTOxWtqnPtFYEXbm9jex/Y+pO/S0aSJxydUzBR2LfAosK2kas1bB5La6GN4XBNE8GwB2915VqnFbW9r+4/tLlNoDdvPkNo7Nyu71d6Z1JF0se1ZNc69yvZ/K/bZ9s9InY7rAe8pOzaH+UPfqtU+S/sKh8e9bQ6FN2IOhVoieIbQ/7rztqti200PJI2W9ClJP5D0c0ndebKbFXOSysluziV1Rh4iacmyfDYnTS9Zc3ic7Qm2x9oeO2zxxXsq3pDWq972EEJDriE9mXaopDNIc8v+1faUWidJOgX4KrV/T9822Y3tlyRdRJpg+1AWrInG8LgmiZpnCP0s35pfTJrO8DzSqgE1HzvOiyieTKpFHgmsCSxWtgRxaTx1taF2pY6jY3Je7yBNBD6DGB7XNBE8Q2iN7rzdg9R7flEP6ffN26/aPtf2I3le2JLCuRVs/4XUO196XPkIUo//JTE8rnkieIbQArbvB+4i3b5PrGPuh9JkN49XHsgdTxv1cH6p9vkZ5k8GHsu8NFG0eYbQIrYrx2XWMpX0JNqRkm4o9crnJ9TOp+ff3atIgbdUg41lXposap4hDEw/JD1M8SHSRDwT8xywDwNLkIJjoTwReflA+Kh1NlkEzxAGINuPkG7NLyF1Cn2YNK5zAumxzlfqyObmvO1xmZfQe3Hb3oFWWO4Vjt3/+nYXo6brvjmq3UXoWbVps5vA9phepr+cKr3m+Smgorlou3j7KrPVlCYB/3ks89J8UfMMYRCStAppiNMs5k/AHZooap4hDCJ55vrRpM6mxYEzbC/QYx/6LoJnCIPLAcCqpGnuvgN8rb3FGbwieIYwiPS2vTU0Lto8QwihARE8QwihARE8QwihAdHmGUKoygvBW0u4rWXYZNQTbb3+f4YXzlUdNc8QQmhEBM8QQmhABM8QQmhAn4OnpHUlnSPpn5KmS3pT0mOS7s5rRH+wGQXtD5LG5zWsx7e7LCGEztKnDiNJ+wO/BEYCTwKTgJeB5YGNSbO/bMf82V1CCGFQaDh4SlqRtB7KSOALwI/z0qel4wsB2+RXCCEMKn2pee4BLAb8wfYPKw/mFfruyK8QQhhU+tLmuULe9rQWy9tI2lzSGZImS3pW0ixJT0m6XNIWBefMa5uU9E5JP5P0hKSZkv4j6TRJixScO0LSFyX9XdIMSc9IukDSanWW9ZJ8rVmSnpd0jaSqtelcRuefj5B0j6RX8/5l8v6VJf1E0kO5PNNzG/GNko6q+x8yDAqSJuXvR1e7yxJ6py/B87G83VHS+r0471uk2/wRwL3MX9P6o8Cdkvatce4qwBRSrfcPpDbWFYAvA5dVJs5NB1cCZwCrA7cCtwM75nxWL7qQpBPyNfYDngGuBh4iLYtwu6Qja5z7Y9KM3zOB6/K1LGml/POnSbX+G4FrSf+WWwDH1/jsIYQBpC+37VcDTwErA3+S9FtSYLqftNhU0TIB3wUOtv1s+U5JHwauAM6RdL3t6VXO/QRwLvDpsgWx1iMF4Q9L2tr2XWXpP00KtE8C42w/lM9ZBLgQOKxaASXtlsv5FPAR2/eUHdsa+A1wlqTbbf+rShaHAlvavrci388BKwI/A46x7bJjCwObVytPCGHgabjmafs1YCdgMikI706aP/Bm4CVJd+Xe+MrzbqwMnHn/tcBE0pKr2xdc9nHguFLgzOf9A7ggv92xIv3n8/Z/S4EznzMDOBZ4k+rG5+0nywNnPvcu4FRSzflTBeefXhk4s3fm7Y3lgTPnO9N2tA+H0CH6NM7T9j9sbwpsDXwbuIU0VGkhYCvgEkndledJGiWpS9J3JZ0rqTunK93+r11wyVttVwt4U/N25bJrjAbWAOYCF1cp+3PAb6uVDdgMeLXa8ez2vN2y4PiVBftLAfU7kvaWtHhBugVIOiq3E09+/aW36j0tdKD8+1DYDlo0PrkvfQOh95oyMYjtu4G7YV474xbA14GdgcPzbfjEfPxTwPdJPfVFlirY/1jB/lfztvyLMTpvnyqvqVaYVmVfqR10KWC2tMC6XOWWL9j/aMH+C0j/JgcBvwbmSPobaUTCJfnfsSrbE0jtqKy6/lLtna0hDHSlvgGRfi+XIg0Z/DLwHmDP9hVt8Gj6rEp5iNLdknYn1bQ2BvYGJkralLSW9GzgS6TOkieA6bYt6dvAiVRZSTCb2+zyVjEsb1+hh7WxgReq7SyoHZf+bQ6W9H+kttit8+uzwGcl/cL2EY0UOoQyve0bCA3otynpbM+RdCspeJZqaB8lBcYf2f5uldPWamIRnszblSWNLKh9jqmyr7RY1lu2u5pYnnls/w34G8yrqe9Oalr4hKRLbRc1F4RQj6p9A5IuAI4h9Q1E8Oyjhts81cP9bLZq3pYm5VsubxdYzU/S8qQV/5oirxj4H9JnPKDe69l+EvgrMErSuGaVp4jtubavI41eAHhff18zDHp19w1UKm9bn/PGG/1TukGiLx1Gx0o6T9JmlQckDc/jID+Wd12at6X/vMMkLVGWfknSo57L9KE81fwob78paY2y6y0MnEVxu2tpxcELJe1ceVDSMEk7qGBQfxFJh0nauMr+dzC/86movTSEevWmb+BtbE+wPdb22GGL192fOST15bZ9BNAFdEl6Bvgz8BKpdrkh8/+6nW77pvzzeaThQxsDj0i6k3Qb/wFgFimAfqIPZar0Y1IHzW7Ag7kZ4XVS4/kipElNFhjrafvqPEj+dOAmSf8C/pnPXRHYiBTojwH+2IvyfAQ4X9KTpH+v/wLvALYlrbH9e1JHUgi19FTpaUXfwJDXl5rnz4F9SDW4x4ENgH1JgfB14HxgW9tfLp1g+2VgLKnX+HXS0zpjSUN7NqbK7Xxf5IlK9gK+QupZ34k0hvSOfN3/1Dj3+8AmpM85jHSL/2FSL/4dwJFUeaqpB98DziQNvh9L+vfakPRgwRHAB23HOKRQaqtcouB4j48Wh/7XcM0zD5K/ip57pCvPe47iweXjmT9AvfycqvvLjncD3QXH3iIN3v9OvdcrO/cB4JNFx6ukr9kObPv3pNplCLWUOjvXrTwgaVFgXEtLE6qKmeRDGHhuydtDJa1T2pkD59nM74gNbRTBM4SBYV47pe07SRPKLAXcn2fcuo7UzLQzqe8gtFkEzxDaa9G8rRwXtC9wGmnKxx1IfQLX521Rb3pooVi3PYQ2yWOlS0PoppUfy5PXnJhflcbT5L6B0HtR8wyhfQ4HRgHPAw+0uSyhl6LmGUILSVqMNJ/rmsx/MOJrtme3r1ShERE8Q2itkcAhwGuk58vPLM04FjpLBM8QWsj2fymeNSx0kGjzDCGEBkTwDCGEBkTwDCGEBkSbZwd65aHFuPHDG7W7GDUtPqkD5oLctt0FGOCGmblLzmlrEaa9sVzPifrRzLnFITJqniGE0IAIniGE0IAIniGE0IAIniGE0IAIniGE0IAIniGE0IAIniGE0IAIniGE0IAIniEMYJImSbKkribk1ZXz6u57yUIEzxBCaEAEzxBCaEAEzxBCaMCQCJ6Sdpd0gaSHJL0uaaakJyRdL+loSUu2u4wh1EtSd612UEnj8/HxrS3Z0DKoZ1WStAJwGbBd3vUP4LfALGA0sBOwO3CqpLG2H21LQUMIHWfQBk9Jy5DWiFkL+ANwtO2/VKRZEjgGOAlYFojgGUKoy6ANnsBPSIHzXmCHvA7229h+DThd0pVAB0xAGUIYKAZlm6ekNYED89ujqwXOcrYfsv10Pndee5KkDSVNlPSMpDmSPl+ZpuD6C7Q5SRqX9/X0GtP3f4EQGiPpKEmTJU2e81rUJ2oZrDXPPUh/GP5q+08N5rE1cA7wJDAJWBKY3ocyPQOcX3BsRWCX/HN7p+4OQ5rtCcAEgIXHjHabizOgDdbguUne3teHPD4JfAs42fbcvhbI9lSgq3K/pCWAO/Lb79t+vK/XCkPeoLyjHGgG6z/y8nn7XB/ymAp8vRmBs4ikYcClwEbAFcCX+utaYVCZlbdLFBxfrVUFGcoGa82zGa623d+30D8hDZX6I3BIrUAt6SjgKIBFhsew1CHuybxdt/KApEWBcS0tzRA1WGuez+ftCn3Io1+HLUn6H+Bo4GFgzzo6tSbYHmt77Mhhi/Vn0cLAd0veHippndLOHDjPBlZtS6mGmMEaPKfk7aZ9yOPNPpxb899V0r7AacCLwO62n6+VPgRg3l2J7TuB64ClgPsl3SjpOuA/wM7Aee0p4tAyWIPn9aQv2waS+mOB84bbnCRtBfwy57GX7X81uWxhcFk0byvHDZX+AD8H7ABsTPrebww81rLSDWGDMnjafojUEQNwtqSFa6WXtKaklXpxiYbanCStBVwNLAwcZvuuXlwzDDGSBKyR304rP2Z7hu0Tba9ue6TtlW0fYfsZ2+Nty/b4inO68/6ulnyAQW5QBs/sM8AjwObArZI2qEwgaXFJx5Nu89/Zi7x73eYk6R3ADcAo4Cu2L+vF9cLQdDjp+/I88ECbyxIqDNredtsvSdqGNDHINsBfJP2dNARpFvAuYDNSLfBZ4KVe5H1nbmPag9Tm9HtgNjCW1FxwHvDxitOOJT0u+jrwnhqzeX/R9gv1liUMLpIWA34GrAlsmXd/zfbs9pUqVDNogydAfuRyW0l7kB7X3BLYlfS5nwd+R7qNvth2b59F2xf4OnAAqc3pBVKb00mkXvRKw/J2CVKNosj4nFcYmkYChwCvkSa2OdP2xPYWKVQzqINnie3rSL2T9aTtosqTQFXSzQBOzK9K4/OrPP0C+0KoZPu/gNpdjtCzwdzmGUII/SaCZwghNCCCZwghNGBItHmGEHpvxMjZrLJqe/su//6X9j5pOmP6yMJjUfMMIYQGRPAMIYQGRPAMIYQGRPAMIYQGRPAMIYQGRPAMIYQGRPAMIYQGxDjPDjRjpeH84+R3tLsYNa2941PtLkII/SpqniG0gKSFJD0myZKelzSi3WUKfRPBM4TW+CCwSv55FLBnG8sSmiCCZwit8Ym8fbLifehQETxD6GeSlgP2AkyaPHsOsIukldtasNAnETxD6H8Hk5Z7mZSXDf4taWWBBVYUkDQmt4tOK8osH3fBsfdLulrSS5LekDRF0id6Oi/0XgTPEPpf6Ra9O29L66pXrnPVJ5J2AP5Aak99FrgGeBWYIOmMZl4rxFClEPqVpI2A95PWJLo8776GtODguyVta/v3TbjOYsCFwCLAN4Dxtp2PbQXc1NdrhLeLmmcI/atU67zM9nQA2zOBiyqO99XHgJWAfwGnlAJnvt7dwE/ryUTSUZImS5o8+5U3m1S0wSmCZwj9RNLCwEH57XkVh0vv95W0RBMut13eXmp7bpXjF9eTie0JtsfaHjt86UWbUKzBK4JnCP1nb2A54N+27yo/YPtPwAPA4sD+TbjWu/L20YLjRftDg6LNM4T+U7olX1rSnVWOr1CW7uf1ZCippwpPUW96tdpo6IMIniH0A0mrADvltyswP1BWs5WkdWz/E5iV9xXdyq9WsP+pHo6PqXH90IABcdsuaVoegzau4PgGkp7OaS6XVLwqU/+Ua0wf8+nO+XQ1p2ShA3SRfr9uta2iF3BZTl+qpT5PCqDvkLR8lXx3L7jeHXm7X0Ht9MDGPkYoMiCCZy2SNgduB1YkNbLvb3tW7bNCaB9JIgVPgAt6SF46fqikYbbfAkpDl07JeZXy3YY0DKmaiaSxnesCJ1Wctznw6V59iNCjAR0886Df3wHLAj8EjrA9p4VF2BFYj/nPI4dQj3HAGsCbwBU9pL2RVNtcifm1ypNJtc9jgAclTZR0L6kSUXXIke03gEOBmaQA+6CkiyXdCtwFnJuTvtXgZwoVBmzwlLQn8BtS288ptr9QPnatFWw/bHtqrg2EUK/SLfhVtl+rldD2bOCS8vPyuMwdgVtIMzGVguphtr9WI6+bga2Aa0nBeG9SxeNY4Ac5WXsXYh9EBmTwlHQw6S/2SODztsdXSbO5pDPygN5nJc2S9FRuE92iIN9hko6WdLekV/I5z0q6X9L3KtuYarV5KtlP0g2Snst5PSnpFkmfrfHZ1so1gmclzZQ0VdKX6+hFDR3C9qG5TfOgnlOD7eNy+n3K9t1peyfbS9pe3PZmti/Kx0rtpdXyut/2nraXtb2Y7Y1sTyDVhgGm9PHjhWzA9bZLOgY4izS04hO2uwuSfov0hXgQuJd0u7IO8FFgb0kH2p5Ycc7PSZMxvAncSforPApYEzie1G70fB1lHJnT7kmaIeePwGPAO4H1gR2AH1c59f3Amfm6t5F6YLcFTgNGA4VBN4SeSFoBWNT2oxX7twBKz7Z3t7pcg9VAC56fJ03dNQs40PaVNdJ+FzjY9rPlOyV9mFRrPUfS9aVH4iStRgqcjwObVjnv/cwf7tGT00mB81/AXranluUzDPhQwXmfA04BvlF6CkTSB0iB9FhJp9t+vM4yhFBpQ+BmSX8D/kP6PVoD2Cgfv8B2T22woU4D7VZxr7z9YQ+BE9s3VgbAvP9aUq1wOWD7skOlcXb3F5z3Z9vP9VTA/Nf9GFLN+CPlgTPnM8f2NQWn30dqv51blv4O0qQNC1WUN4TemgqcTfoubUP6fVoNuBU4jCpT4IXGDbSa5+2kZ3RPkDTF9mW1EksaBexBulVehvmfZ/28XRu4Pv88lTSzzYckfRW4qPL2pk47kNpi77L9YC/P/U1Bp9dUYDegcHJcSUcBRwEMe8cyvbxsGApsP0HqHAotMNCC53jSLe8XgYskURRAJX0K+D6wWI38lir9YPu1PCnsL0jtpd+S9CRp/sPrgUtsz6ijjKUnOKbWTFXdYwX7X83bRYpOzI3+EwAWXn10TGgbQpsNtNt2bH+J1J45nBRA96tMI2lT0u3JCOBLpIHBSwAL5V7I/yslrcj7cmBV0gDmXwCvk6byOg+Ymh+p67GIvf9U88TzxSEMEgMueELVALpvRZKPkgLjj2x/1/Y/bb9Rdku8Vo28/2v7fNtH2F43p72NVKP8Th3FK9Ue1+nFRwohDDID7bZ9Httfyk+YfRG4ON/Cl4YeLZe3C/RM57GaH+zFdR6W9C1SZ8376jjlVtJTGltJWs/2P+q9VgidRMBCbV7yaIlVX+05UT9aaGTxA40DsuZZUlEDvVjSx/KhUnvjYeUTyUpaknQ7vkxlXpI2krS/pGozvH44b3vsQMo98ueQ/u2ukLR2xXWG5eFSIYRBbMDWPEsqaqC/yj+fRxoTujHwSJ4rUcAHSGPbfsGCyxusRnoMbrqk+0m11pGkMXBrkHriT66zWF8iDazfnfQM8R+AJ0jDoTbI26pPgIQQBocBXfMsqaiB/or03O9YUu/z66Qe+rHAlaSAWm2g+R+BE0lTd40mPfe7EzAd+B6wge3JdZZnJqm2emjOb31Sx9O6wF+IGWxCGPTU4rk2QhMsvPpor3TKZ9pdjJrWPvKv7S5Cj26edfEU22PbXY6Basm1V/RGPz20rWV4aXp711F6+PhzefOhp6reRXZEzTOEEAaaCJ4hhNCACJ4hhNCAAd/bHsJgIGld0giR7UkTHIs0/eETpEeEb8yTGYcOEcEzhH4maX/gl6ShcU8Ck4CXgeVJo0O2JE2IE8Gzg0TwDKEfSVqRNO54JPAF4Mfl63DlFQS2ya/QQSJ4htC/9iDN/PUH2z+sPJjndr2D+UsHhw4RHUYh9K/SJNw9TrRdSdLikv5H0n2SXpX0pqQHJY0vfyy5LP34vObWeEmrS7owr5U1I593gqSoMDVJBM8Q+ldpFq4dJa1fM2UZSaNJa3N9h/Ro8R+A35JWw/w6cJekZQtOXx2YTOqcmkSaNWwN0lN6E2OxweaIf8QQ+tfVpLWxlgD+JOn6XJvcSdLS1U5QmsDhMuA9wE+A1WzvYnsv0pwKF5LWK/pBtfNJS27cBqxpe3/bu+X0T5IeSz66aZ9uCIvgGUI/yuu270SqCQ4nTSbzHVLP+kuS7sq98eV2JfXA/xH4nO03y/J7kxT8ngMOLqh9TgeOLV8Zwfa/gdKa719oxmcb6qL9owPpLTH8qZHtLkZNCy29ZLuL0LMeF5lujjzn66aStiJNYrM5aYjSssBWpLlhd7PdlU/ZPW+vKF8ssCy/NyRNzuk2Jd3Ol7u5YDHDi4FzgbUkvcv2k5UJytfKWniFDvg/bKOoeYbQIrbvtn2S7Z2AUcDWzA98h5etmLBG3p6RO4AWeDE/wC5f5VL/Kbj+TODp/HZ0QZoJtsfaHjti6VrLg4WoeYbQBrlGebek3UkdQxuT2iMnAsNystuBaT1k1cgKsKEJIniG0Ea250i6lRQ8S7XI0ny0E22f1UC2Y6rtlDQSWCm/XeCWPfRO3LaH0I9yz3lPVs3bJ/L2hrytXPiwXjtLGlVl/4Gk3/mH8xrvoQ8ieIbQv46VdJ6kzSoPSBou6UjSKgQAl+btVcAUYDtJ50harsq5K+Zzq1kMOEvSwmXp1wROzW/PbOyjhHJx2x5C/xoBdAFdkp4B/gy8RFoBdkNg5ZzudNs3QWoPlbQ38BvgU8BBkh4g3c4vAqxNGgP6HPD/qlzzAlKv/sOS7gKWJA2YXwS4FmikKSBUiOAZQv/6OanTZydgM+YvEPgW6Tb9fOBc23eWn2T7iVxbPQLYL5+3OfAiqb3ye8CvC675CGkI07eBHYCl875fAD+sNvwp9F4EzxD6UR4kf1V+9fbcGaRaYq9rirYfAQ7o7XmhftHmGUIIDYjgGUIIDYjgGUIIDYg2zxAGCdvjgfFtLsaQ0TE1T0nTqjzjO0PSY5Iuk7RdE681Kec/rpfndefzuppVlhDCwNSJNc+bgGfyz8sCG5GexNhX0vG2i+Y4DCGEpunE4Hma7UmlN5JGkJ6YOAY4TdLENj56diJwGvNnrgkhDFIdc9texPZbwAnAa6QVCnduY1metj3V9ivtKkMIoTU6sea5ANtvSvoXsAnwztJ+SdNI67+sbnta5XmSJpHWy96+vDZbkWZ74KSc9wjS43Wn276mStpu4HDg47a7K44JOIT0xMj7gEVJT4rcSKpNP04IA4iBt+a2t341/aGqK5W0zNyZwwqPdXzNs0zpX/nZJua5D/A70uN0NwAPkCawvVrS8fVmkgPnhcAvSTOH30d64kTAscCfJW3axHKHEPrZoAiekt5LWjFwFgsuSdAXxwFftr2h7YNsbw3sCcwGTpe0YZ35HAMcRArsG9ne2fYBwLuBH5MmiZhYPgtOCGFg6+jgKWlZSbsBV5I+y+eb3Fk02fZ3y3fYvpa0Fsww4LN15nNC3n7N9oNlec0BvkhannY15k9NtgBJR0maLGnynDfe6MVHCCH0h04MnreVrePyEmnartWA3Wyf3eRrXVSw/4K8HddTBnn97TWAuWXnzWN7Vtl1CvMrX1tm2OKL93TZEEI/68QOo9I4TwErAh8gzVP4S0lb236oidequpAW89eVqbqIVoV35e3T5UvBVnikIm0IYYDrxOBZOc5zJVJA3QC4SNIWtl1nXq2seddbphBCB+jE2/a3sf00abLYt0iTzR5cdnhW3i5RcPpqPWQ/pof99SyiVUqzco0OoTUq0oYQBriOD54AtqcCP81vx0sq1ahLwWjdynMkrQ+s0kPWB/ewf1IdZXuCdFu+EGmcZ2U5RvQmvzA4SFpI0kGSrpb0pKSZkl7KnYKnSlqh3WUMtQ2K4Jl9i/SU0ZrAoXnfLXn7P5KWKiWUtArQTWo3rWVTSV8o35HX2T4EmAP8pM6yfT9vT5U0L5BLGgacTlo98VHg8jrzCx0sdyLeS+oo3IM02uJK4G7SkLv/Ja0/tF8TrtWVO1i7+5pXeLtBEzxtPw+UhhX9b659nkVaNGtT4J+SrsxrZP8DeJX0Za3lR8B3JT0g6WJJdwLXk9qKT7T95zqL91PgV6Q1sx+QdJOkXwH/Aj4PvAzsa3tmnfmFDpVXwvw96Ym1ScBatre0faDtPUidoCeSVsC8RNJH21bYUNOgCZ7Z90kD0dcADrf9MumJoItIAe9DpB7yM4BdSe2ktfwa2IW06NaHSDM43Q3sY/uMeguVO7AOBg4D7iEt5PUR0r//2cD7bN9Xb36ho51FajO/jzS87m0jOmy/Zfs04HjSndEvCtZgD23WMb3ttsfUkeZ10l/u8n2PU6WtMRtXkE/l/t/1WMB0Xhdpmdlqx0wa57nAWM8wNOS100u34sfWGLoG6a7nCNIoks9QMcmxpF1IyxJvAYwijXl+GLgG+FGe72Ea8ztFD5d0eFkW5+fva2jQYKt5hjCQ7UH6nXvQ9uRaCfMf21/mt3uW9is5mzShzD6kTtErSPMurEKaErE0Oc7lwF3554dJyxyXXm9b6jj0XsfUPEMYBDbJ23vrTF9qynmfpOG2ZwOfA44mNU/tbfuPpcR5AprtSW3o2P5iXtVga+DOqGk2VwTPEFpn+bytd+avUrqFgOUkvUSaHhGgqzxwwrza6q19LmWoS9y2hzBwVQ6lG0tq33zC9o39csGyCWjeeuXN/rjEoBHBM4TWeSFv31kz1XylgfJzSR1Cpc6ffzazUOXKJ6AZsfSi/XWZQSGCZwitMyVvt6gz/WZ5+0Bu74z5EQaQCJ4htM51pFrkej2tHJA7fw7Lb6/N28fydp3+KV7ojQieIbRIni6x9AjuWZIWqZH8OGB90iPHZ+V9U0i3/qPzOM96lCbHic7hJovgGUJrfZr5jwz/RtKY8oOSRkj6MulpOQOftP0czFsp9v9y0vMkbVZxriRtL6l81bTS5DjrNf2TDHHx1yiEFrL9gqRtgKtJYzIfknQPaWKYJUkLBC4HvAEcafuyiix+QAqEnwT+KGky8FA+5z2kgfKrA6Xlr/9Imjx845z2QdJjyXfZPq/fPugQEMEzhBaz/ZikscCBwP6kwfObkgLmI6T5Dn5i+5kq5xo4UtLVpMHymwHvJ82/8BBpQcFnytLPlLQradaxLUnzMyxE+t2P4NkHETw70GJLzmDsdlPbXYyaXjkrhrnUkhf/uzC/Gjn/OlIHVD1pHyA9GhqaKNo8QwihARE8QwihARE8QwihARE8QwihARE8QwihARE8QwihARE8QwihATHOM4RQ1ezpI3hhSr2z5/WPNU/6Q1uv/4LfKDwWNc8QQmhABM8QQmhABM8QQmhABM8QQmjAgAiektzAq7vd5Q4hDF0Dpbf9/Cr7VgR2IU3TdXmV43f2a4lCCKGGARE8bXdV7pM0jhQ8X6h2PIQQ2mlA3LaHEEKn6cjgKak7t3t2SdpQ0kRJz0iaI+nzlWkK8hifj48v2i9pdM7naUnTJd0v6WNlabeW9BtJL+bjt1VbFVHSmJznNEnDJX1F0j8kzZD0rKTzJa3a1H+k0Hb5/7vURl84GbGkv+U041pXutBXHRk8y2wN3AtsDEwCbgSmNynvMaTVCrcFbgfuJy1hcJmkAyTtA9wGjAJuJq1BMw64TdLaNfK9FDiFtIzsVcBM0hKz90mKJWUHr29L6vTft1BmQLR59sEnSWuznGx7bpPzPhw4EzghL5mApGOAnwJnAIsDB9uemI8tBFxMWpPmy8ARVfJcDVgU2Mj23/N5I4GfA4cAF5DWpAmDy3RgA+Bg0v9xGAQ6/S/hVODr/RA4AaYB/1MKnNkE0kJbo4EbS4ETIJfhO/nt9jXyPbUUOPN5s4DPAq8Cm0rautpJko6SNFnS5Jkvz2jk84T2+VHenpL/WIZBoNOD59UVwa2ZbsuBbZ58rWn57Y1Vzvl33q5cI98FFvyy/V/g2vx2XLWTbE+wPdb22IWXXaRG9mEAuoLUvLQ6acXLuuQ13D8j6R5Jr0p6M7eVnybpHRVpd83tpn+qkd9ykmbm13INf5oAdH7wfLQf836iYP/rRcdtl44tXHDuf3OgrGZa3o6up3Ch45yYtydJWqKnxJIWAX5LWkp4feAO0h/YZUjNQlMkrVF2ys3AU8D7JW1YkO2BwEjgWtsvNfIhwnydHjzf7MO5PX32npoC+qOpIAxStm8lBcMVgBPqOOUbpLuQqcDatvewvR+wBqkmuxpwUVn+c4Bf5rddBXkenrfdvSt9qKbTg2ctpVvuor/yq7WqIGWWkbR0wbExeftki8oSWu9EwMAJkpYvSiRpUeCY/PY42/O+E7bfJN36vw5sUdFG3p23B0t6W2ewpPcAmwLPUL3JqZRuXtv6nDeK57IMgzt4lr5w61YeyF/OcS0tzXwHV+7IAbU0DnBSS0sTWsb2/cBlwJLASTWSbkL6o/+U7Zur5PMCVdrIbf8T+COpdrt7xWmlWudFtmfXKOO8tvVhiy9e+wMNcYM5eN6St4eWj5/MgfNsoF2D0k+WtF5ZeUaQhkQtDUyxHc/sD27/C8wGjpZUdPfzrrz9T418HqlIW3Je3paCJZKGkYbCQdyyN82gDZ45CF0HLAXcL+lGSdeRvpA7M/9L1kqPAXcDf5Z0g6RLgIdJX/QXSIPlwyBm+yHgXFKn4jd6St7AJS4l9QXsUdYjvxNpBMgU239rIM9QxaANntm+wGnAc8AOpCeRrs/bx9pQHgP7AaeSGv73BhYhDV/atHz8ZxjUvkEaOH+IpPWrHC81Oa1eI49ST/vb2shtv0J6cm0kcFDe3ZW33b0vaigyYIOn7Um2ZXtMlWNd+Vh3D3nMsH2i7dVtj7S9su0jbD9je3zOY3zFOVX3lx0fl49PKjgu26pRptm2v2l7HduL2F7B9qG2p9X6LGHwsP00qalmIeDbVZJMIXUIvUvSjpUHc43yw/ntpCrnz7t1z+3pe5M6UC/uU8HD2wzY4BnCIHc68DIpCL6thpl71M/Jb8+UtFLpWB7/eTapQ+mPtu+qkvctwOOkjqdvkO5uYmxnk0XwDKEN8sMS/5ffLlYlyddItcr3Av+WdI2kS0kdRfuSmp0WGLmR857L/Gfoj8vb7maUO8wXwTOE9vkxBU+y2Z5B6tg8Dvg7ab6EvUhzIJwObGz7kWrnZt1lP9cc2xka0+mzKnWE3J5Z2A4aBqdq7fUVx2cAq9Q4/hYpwP64gWv/m/jO9auoeYYQQgMieIYQQgMieIYQQgMieIYQQgMieIYQQgMieIYQQgMieIYQQgNkNzJxS2gnSc/T/CVIRpFmdhqo+qN8q9kunJR4qGvC92wgfKf6WobC70gEzwCApMm2x7a7HEUGevnCggbC/1l/liFu20MIoQERPEMIoQERPEPJhHYXoAcDvXxhQQPh/6zfyhBtniGE0ICoeYYQQgMieA5BktaR9DlJF0qaKmmuJEv6WLvLBmlFUUk7SvpeXkP8VUmzJD0p6XJJ49pdxrCgdn+vWv29ifk8h6ZjgM+1uxA1bAeU1it/BrgDeAN4D/BR4KOSTrV9cpvKF6pr9/eqpd+bqHkOTX8DzgD2B9YCbm9vcRYwF7gC+IDtlWzvYXt/2xsABwBzgK9J2r6tpQyV2v29aun3JmqeQ5Dtc8vfSwNrwnHbtwK3Fhy7VNIHgSOAQ4DbWlm2UKzd36tWf2+i5hk60Z/ydnRbSxE6TVO/NxE8Qyd6d94+3dZShE7T1O9NBM/QUSStCHTlt1e0sSihg/TH9yaCZ+gYkoYDFwJLA7fYvrbNRQodoL++NxE8Qyc5B9gReJzU6B9CPfrlexPBM3QESWeSekqfAXa0/UybixQ6QH9+byJ4hgFP0veA44DnSb8A/25zkUIH6O/vTQTPMKBJOh04HngR2Mn239tcpNABWvG9ieAZBixJpwFfAl4GPmj7L20uUugArfreRPAMA5KkbwJfBv5L+gX4U+0zQmjt9ybm8xyCJG0M/LRs13uAJYF/Ay+VdtreosVFA0DSnsDV+e1k4MGCpFNtn9aaUoWetPt71ervTTzbPjQtBWxeZf+7q+xrh+XKfh6bX9XcDkTwHDja/b1q6fcmap4hhNCAaPMMIYQGRPAMIYQGRPAMIYQGRPAMIYQGRPAMIYQGRPAMIYQGRPAMIYQGRPAMIYQGRPAMIYQGRPAMIYQG/H+d2VNO6DSNEgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 345.6x360 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "w = calculate_W(np.dot(A_A1,A_A2), np.dot(B_A1,B_A2), np.dot(C_A1,C_A2), A_S2, B_S2, C_S2)\n",
    "\n",
    "fig, axs = plt.subplots(1, 2, constrained_layout=True, figsize=(4.8,5))\n",
    "color = \"viridis\"\n",
    "\n",
    "A_new = np.dot(A_A1,A_A2) @ w[:,:,0].T\n",
    "B_new = np.dot(B_A1,B_A2) @ w[:,:,1].T\n",
    "\n",
    "plt.sca(axs[0])\n",
    "plt.xticks(range(2), range(1,3), fontsize=22)\n",
    "plt.yticks(range(8), ['Clinton', 'Kaine', 'O\\'Malley', 'Sanders', 'Cruz', 'Kasich', 'Rubio', 'Trump'], fontsize=22)\n",
    "plt.sca(axs[1])\n",
    "plt.xticks(range(2), range(1,3), fontsize=22)\n",
    "plt.yticks(range(10), ['Feb', 'Mar', 'Apr', 'May', 'Jun', 'Jul', 'Aug', 'Sep', 'Oct', 'Nov'], fontsize=22)\n",
    "\n",
    "axs[0].imshow(A_new, cmap=plt.get_cmap(color))\n",
    "axs[1].imshow(B_new, cmap=plt.get_cmap(color))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## S Matrices\n",
    "Note that these S matrices will not look as nice as the S matrices in the paper. This is because the order of topics in NCPD is arbitrary, so to get a \"nice\" S matrix one needs to reorder the topics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW8AAADOCAYAAADxC4PrAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAOf0lEQVR4nO3de7AkZX3G8e/jclkuASQSQSnBSARTalA35SVG5Q9NRRNNaaBIajUIuDFCKiGWpkRUBCIVk6BVEqNriBKCKZUo5YXgHaVCEdgohTcERUASuUVFXJQs8Msf3auHYXZPn92Z0+eF76dqqne6Z08/Neec5/S8/fZMqgpJUlseMnYASdLSWd6S1CDLW5IaZHlLUoMsb0lq0A7LtaOdsnOtZrfl2p1WiOy049gRFvUrh9w+doRBrvn6HmNHGKTuvnvsCA8od/CD26pqn8n1y1beq9mNp6563nLtbtvde8/YCR5Qdth3/7EjLOoTn/z42BEGef6vPXfsCIPcc+utY0d4QPlMnXf9tPUOm0hSgyxvSWqQ5S1JDbK8JalBlrckNcjylqQGWd6S1CDLW5IaZHlLUoMsb0lqkOUtSQ2yvCWpQZa3JDXI8pakBg0q7ySvS/KhJNcmqSTXzTmXJGkrhr6f91uA7wNfAvaaWxpJ0iBDy/sxVXUtQJKvArvPL5IkaTGDhk02F7ckaWXwhKUkNWiun2GZZB2wDmA1u85zV5L0oDLXI++qWl9Va6pqzY7sPM9dSdKDisMmktQgy1uSGmR5S1KDLG9JatCg2SZJXgoc0N/dB9gpyUn9/eur6px5hJMkTTd0quAxwLMn1p3aL78AWN6StIwGlXdVPWfOOSRJS+CYtyQ1yPKWpAZZ3pLUIMtbkhpkeUtSgyxvSWqQ5S1JDbK8JalBlrckNcjylqQGWd6S1CDLW5IaNNcPIF7o7n1245Yjnrpcu9tmv/TOS8aO8IDy9ZP2GzvCon7rEYeOHWGQu377wLEjDHLTy/cdO8Igd/24kc/Vffl5U1d75C1JDbK8JalBlrckNcjylqQGWd6S1CDLW5IaZHlLUoMsb0lqkOUtSQ2yvCWpQZa3JDXI8pakBlnektQgy1uSGmR5S1KDBpV3kscmOSXJpUluTXJHkiuSvD7JbvMOKUm6r6FH3kcDJwDfBk4BXgN8EzgNuCTJLvOJJ0maZugn6ZwHnF5Vty9Y964k1wCvB44Bzpx1OEnSdIOOvKtqw0Rxb/aBfvn42UWSJC1me09Y7t8vb97eIJKk4ba5vJOsAt4A3A28fwuPWZdkQ5INd/9k47buSpI0YXuOvN8OPB14Y1V9c9oDqmp9Va2pqjU77OKkFEmalW0q7ySnAscD66vq9NlGkiQtZsnlneRk4CTgvcArZx1IkrS4JZV3X9xvAs4Gjq2qmkcoSdLWDS7vJG+kK+5zgKOr6t65pZIkbdWgi3SSHAe8GbgB+Azwh0kWPuTmqvr07ONJkqYZeoXlr/fLR9ENmUz6AmB5S9IyGXqF5VFVla3cnjPnnJKkBXxLWElqkOUtSQ2yvCWpQZa3JDXI8pakBlnektQgy1uSGmR5S1KDLG9JapDlLUkNsrwlqUGWtyQ1aOi7Cm7/jm7byMPXX7Zcu9tma664Z+wIg1x+6KqxIwzy2D++fOwIDxg7/WjT2BEGOeCIr4wd4QHlhi2s98hbkhpkeUtSgyxvSWqQ5S1JDbK8JalBlrckNcjylqQGWd6S1CDLW5IaZHlLUoMsb0lqkOUtSQ2yvCWpQZa3JDXI8pakBg0q7yQHJzk3yTeS3J7kziRXJTkjyX7zDilJuq+hH8awP7Af8BHgRuBu4AnAOuDIJIdW1S3ziShJmjSovKvqs8BnJ9cn+SLwQeAo4K0zTSZJ2qLtHfO+vl8+dHuDSJKGW9JnWCZZDewOrAZ+FfjrftMFM84lSdqKpR55HwvcCnwX+CSwF7C2qi6e9uAk65JsSLJhU921XUElST+31E+PPx+4iu7o+0nAC4GHbenBVbUeWA+wx0P2rm2LKEmatKTyrqob6WabAJyf5N+Ay5PsWlWnzzydJGmq7TphWVVXAl8GXjWbOJKkIWZxheUuwN4z+DqSpIGGXmG57xbWHwY8Hrh0lqEkSVs3dMz7H/rL4D9HN7d7NfAU4EjgDuDV84knSZpmaHn/K/Ay4KXAPkDRlfi7gb+pqhvmE0+SNM3Qy+M/SHcZvCRpBfAtYSWpQZa3JDXI8pakBlnektQgy1uSGmR5S1KDLG9JapDlLUkNsrwlqUGWtyQ1yPKWpAZZ3pLUoKV+huV2CKxatXy720aXH7ryMwIcc/V3xo4wyFmHHDR2hMXde8/YCQbZ4WttfM/beDbb55G3JDXI8pakBlnektQgy1uSGmR5S1KDLG9JapDlLUkNsrwlqUGWtyQ1yPKWpAZZ3pLUIMtbkhpkeUtSgyxvSWqQ5S1JDdqm8k6ya5Jrk1SSM2cdSpK0ddt65H0KsM8sg0iShltyeSd5MvDnwJtmnkaSNMiSyjvJKuA9wIXAh+eSSJK0qKV+huUJwCHAS+aQRZI00OAj7ySPBt4MnFJV1w38P+uSbEiyYVP9dBsjSpImLWXY5F3AtcAZQ/9DVa2vqjVVtWbHrF5yOEnSdIOGTZKsBZ4LPKuqNs03kiRpMYuWd5Kd6Y62LwBuSnJQv+mR/XLPft1tVfXDuaSUJN3HkGGTXejmdL8AuGbB7aJ++9r+/rFzyCdJmmLIsMlG4PAp6/cB3kk3bfAs4MoZ5pIkbcWi5d2PcZ83uT7Jgf0/v11V99suSZof35hKkhq01It0fqaf653ZRZEkDeWRtyQ1yPKWpAZZ3pLUIMtbkhpkeUtSgyxvSWqQ5S1JDbK8JalBlrckNcjylqQGWd6S1CDLW5IalKpanh0ltwLXz/jLPgy4bcZfc9ZayAjmnDVzztaDOecBVbXP5MplK+95SLKhqtaMnWNrWsgI5pw1c86WOe/PYRNJapDlLUkNar28148dYIAWMoI5Z82cs2XOCU2PeUvSg1XrR96S9KBkeUtSgyxvSWpQU+Wd5HVJPpTk2iSV5LqxM01K8tgkpyS5NMmtSe5IckWS1yfZbex8myU5OMm5Sb6R5PYkdya5KskZSfYbO9+WJNl1wff/zLHzLNRnmnb78djZFkqyd5K/TfKtJD/tf04/n+Q3x84GkOTkrTyXlWTT2Bk3S7J7khOTfKX/Xb8tySVJjkqSee57h3l+8Tl4C/B94EvAXuNG2aKjgeOAjwLnApuAw4DTgCOSPK2qfjJivs32B/YDPgLcCNwNPAFYBxyZ5NCqumXEfFtyCnC/q81WkIu5/4yDlVQ2BwAXAbsDZwFXA3sCTwQeOV6y+/gw8K0p658IvAb42PLGmS7JQ4B/B54BnA28A9gV+APgvcDjgL+cW4CqauYG/PKCf38VuG7sTFMyrgH2nLL+NKCA48fOuEj+w/ucrx07y5RsT6b7I/MXfcYzx840ka+A942dY5GMFwPfBfYbO8s2ZH93/xy/YOwsfZ6n93neNrF+J+Ba4Ifz3H9TwyZVde3YGRZTVRuq6vYpmz7QLx+/nHm2web3n3noqCkmJFkFvAe4kO7IbMVKslOS3cfOMSnJs4BnAm+tqu8l2THJrmPnGqIfcjyS7lXihSPH2WyPfvk/C1dW1f/Rvb/JxnnuvKnybtz+/fLmUVNMSLI6ycOS7J/keXRHNwAXjJlrihOAQ4Djxw6yiN8H7gTuSHJLknck2XPsUL3n98sbknwM+AmwMcnVSdaOmGuIw+nK8n1Vdc/YYXqXAT8EXpvk8CSPSnJIktOBpwAnz3PnrY15N6k/anwD3Uv+948cZ9KxdGN1m10HrK2qi8eJc39JHg28GTilqq5LcuDIkbbkMuBDdOO1e9CV5fHAs5M8o6rGPnF5cL98D3AN8Ed0L/FfDZyTZMeqeu9Y4RZxDN0QxT+NHWSzqvpBkhcC/wh8cMGmO4CXVNX589y/5b083k43PnZiVX1z5CyTzgeuojuB9STghXRva7mSvItuDPGMsYNsTVU9dWLVPye5Evgr4M/65Zh+oV/eARzWv7wnyfl0z+9bkpxdVfeOlG+qJAfTDfd8tqq+M3aeCT+mO//2UeASYG+6CQvvT/Kiqvr03PY89qD/dpwsWJEnLKfkPJXuiOHdY2cZmPeJwF3A68bO0udZC9wLPHPBugNZgScst5B/x/75vGQFZPlY/7ydNmXb2f22x42dc0q2t/bZjhw7y0SuJ9ANPb1yYv2udH8MrwNWzWv/jnnPUZKTgZPopg29ctw0w1TVlcCXgVeNnSXJznRH2xcANyU5KMlBwAH9Q/bs1+01VsbFVNUmuhNaK+HVzI398qYp277XL1faieodgJcB/0s3rXUlOQFYTTdU9jNVdSfwCbqf0wPntXPLe0764n4T3RHNsdX/SW7ELnQv/8a2C92c7hfQjdFuvl3Ub1/b3z92jHBDJFlNd7J6JZyovqxf7j9l2+Z1K21u/+8CDwf+paruGjvMhM3z4ldN2bbDxHLmLO85SPJGuuI+Bzi6VtgYIkCSfbew/jC66YyXLm+iqTbSzTKYvG1+VXBhf/+jo6RbIMkvbmHTqXS/wCvhwpLz6ca71y6cythfUft7wNVVNe3imDEd0y/PGjXFdF/vl0ctXNm/EnwR8AOmX2w0E029JWySl/Lzl8x/Snem/O/6+9dX1TmjBFsgyXHAmcANdDNMJov75prnSYyBknyE7grLz9HN7V5NN73pSLqpbs+pqitGC7gV/WyT7wB/X1UrYupgkrcBTwM+T/e9351utslhwH/SnSAc/craJOvopoN+jW7mxk7An9D9LPxOVX1qxHj3keQRdM/lf9X9TwaPrr9a9Ut0Q03nAv9B94r1FXTDJcdV1TvnFmDsQf8lniC4iO7ExbTbRWPn6zO+bysZV1LOI4CP011t91O6Ey9X0U0bfNTY+RbJfiAr7IQl3ZHWJ4H/7p/PjcAVwInA6rHzTWR9Md0rq410R+KfAn5j7FxTcp7Yf59fMXaWrWR8DN3Q6I10b4PwI+CLwIvnve+mjrwlSR3HvCWpQZa3JDXI8pakBlnektQgy1uSGmR5S1KDLG9JapDlLUkNsrwlqUH/D8rVsP2zu058AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.xticks(range(0,8), range(1,9), fontsize=18)\n",
    "plt.yticks(range(0,4), range(1,5), fontsize=18)\n",
    "\n",
    "plt.imshow(C_S1)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAB8CAYAAABwrOvwAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAALlklEQVR4nO3de4xcZR3G8e/Tbam9iKWloHIpComgjZG0alREQCwIokbxgoIWacA0hgTEgBpFMd4QTCMRK4rUyEVECWKUoChViBipgkQECkgFuUgRili0C+3PP953dVxmd8+WOfPOG55PMjm750x3np6defbdc847q4jAzMzqM6V0ADMz2zoucDOzSrnAzcwq5QI3M6uUC9zMrFIucDOzSk3t1wMNzZwV0+bM7dfDbbXpGzaXjtDIlql1/OzdsuPg78+hOzaVjtDI5nmzSkdoZPr8f5eO0Mjw5qHSERrZePvfHoqI+d229a3Ap82Zy27LTuzXw221BZc/XDpCI8Pb1/Fi/seJj5WOMKG5b1pbOkIjjxz2qtIRGtnj2FtLR2jk7se2Kx2hkesOOv0vY22rYxhnZmZP4QI3M6uUC9zMrFIucDOzSrnAzcwq5QI3M6uUC9zMrFIucDOzSrnAzcwq5QI3M6uUC9zMrFIucDOzSrnAzcwq1ajAJR0l6RuSbpA0LCkk7dduNDMzG0/Tt5P9DLAAeDDfdmotkZmZNdL0EMoxwK4RsSPw3RbzmJlZQ41G4BHx87aDmJnZ5PgkpplZpVzgZmaVarXAJR0raY2kNZsf39jmQ5mZPeO0WuARcU5ELI6IxUMz6/gjvGZmtfAhFDOzSrnAzcwq5QI3M6uUC9zMrFKNJvJIWgbskz9dnJenSFqaP/5mRFzb42xmZjaOpu+Fsg/w/lHrDur4eDXgAjcz66OmU+mXAktbTWJmZpPiY+BmZpVygZuZVcoFbmZWKRe4mVmlXOBmZpVygZuZVcoFbmZWKRe4mVmlXOBmZpVygZuZVcoFbmZWKRe4mVmlmr4b4dO2zf0b2eWz1/Xr4bbejBmlEzQy9abHS0doZIe1O5WOMKEf33dj6QiNHLL3DqUjNPL3bz9SOkIjs6gj53g8Ajczq5QL3MysUi5wM7NKucDNzCrlAjczq5QL3MysUi5wM7NKucDNzCrlAjczq5QL3MysUi5wM7NKucDNzCrlAjczq1SjApe0s6QTJF0l6R5Jw5LulXShpIVthzQzs6dqOgL/EPBlYFfgivzxDcC7gTWS9m8nnpmZjaXp+4H/Ftg3Iq7pXCnpHcD3gLOBvXqczczMxtGowCPi0jHWXyJpLbCnpO0j4qGepjMzszH14iTmcF4+2YOvZWZmDT2tApe0CFgIrImIDT1JZGZmjWx1gUuaDXwbCODkniUyM7NGtuqPGkvaBrgEeAlwakT8Yoz7HQscC/AsZm5tRjMz62LSI3BJU4GLgYOBMyPitLHuGxHnRMTiiFg8jelPI6aZmY02qQLP5X0R8FbgrIg4qY1QZmY2scYFLmkI+A5wOLAyIo5vLZWZmU2o6VT6KcAq0szLc4HlLWYyM7MGmp7EPBU4EtgA3AecKmn0fVb4UkIzs/5pWuAL8nIO8Ikx7rOKVPBmZtYHjQ6hRMTSiNAEt3UtZzUzsw5+P3Azs0q5wM3MKuUCNzOrlAvczKxSLnAzs0q5wM3MKuUCNzOrlAvczKxSLnAzs0q5wM3MKuUCNzOrlAvczKxSLnAzs0opIvrzQNJ64C89/rLbAw/1+Gu2wTl7yzl7q4acNWSEdnIuiIj53Tb0rcDbIGlNRCwunWMiztlbztlbNeSsISP0P6cPoZiZVcoFbmZWqdoL/JzSARpyzt5yzt6qIWcNGaHPOas+Bm5m9kxW+wjczOwZywVuZlapqgpc0lGSviHpBknDkkLSfqVzdZK0s6QTJF0l6Z6c815JF0paWDrfCElzJH1F0m8k/U3SJkl3S/qxpANK55uIpMvy939grg3Oeca6LSudr5OS90m6RtKjkv4p6WZJZ5fOBiDpUxPsz5B01ADknCbpg5Kul/SIpA25nz4saUbrj1/TMXBJ64AFwIPAE8BOwP4RsbpgrP8j6QvAycDtwGrgYWAhcAgwDLwxIq4uFjCTtAdwI3AdcAfwCPB84K3Ac4CPRMQZpfKNR9IRwPmk/bkxIrYvHAlIBU6arLaqy+bLI+L3/U3UnaQh4DvAEcANpOfpZuCFwOsGYX/mgdl+3TYBHwWGgF0j4t7+peoSRrocOAy4Gbgqr34D8GLgV6R+2tJagIio5ga8Htglf3wGEMB+pXONyvg24LVd1r8j572ldMacZwiY2mX984AHgMeBmaVzdsm3A7AeWAGsAx4qnakjWwCrS+dokPOUnPXDXbY95TkxSDdg/5z9igHI8sqc5RfAlI71Q6Qfiq33U1WHUCLi5xFxT+kc44mISyPimi7rLwHWAntKKj7CiYjNEfFkl/X3A78GZpDKfNB8FdgIfLx0kBpJmkUawa6OiDNHb+/2nBgwS/PyvJIhshfk5U+jY5QdEZuBK/Onrb7Wp7b5xe0phvNyYF8kkuYCrwD+AQzUD0tJbwcOBw6OiI2SSkfqZjtJxwHzSIf6ro6IOwtn6rQE2Bb4gaRtgbcAOwP3kUa1D5YMNx5Js4G3kw73/bBwHIA/5eUSSaePlHg+RHUQsAn4TZsBXOB9ImkR6Vj4mojYUDjOf0naAVhOOqH9PODNwFzgAxExPN6/7SdJ80ij7/Mj4sqJ7l/QS4GVHZ+HpG8Bywdkfy7Ky+2A24DndmzbKOm4iLig/7EaeScwC1gVEZtKh4mIm/JJ3+XATZJ+ljctIe3XIyPir22HqPLGgB4DHyPrbOCPwBbggNJ5RmVbmPfjyO0x4L2lc3XJeQFpRDuvY906BusY+OnAy4E5pIJ8A/C7vF/PLp0vZ1yZ8zwJ/Ah4Eemk9btII9sngJeVzjlG9l/l7ItKZxmV62Okk8Ajr6EtwNeAHVt/7NL/+aex06oocGAb4Iqc9ZOl84yTcyqwO/D5/AT8UulMHdkOy/vvPaPWD1SBj5H9OcDduTCfOwB5zsn78l5gxqhtx+Vt55bO2SX37vl5eVPpLB2ZppCOxT8KHA3MJ/32egTpRPtdwNw2M1R1ErM2kqYCFwMHA2dGxGmFI40pIp6MiDsj4qOkQxUnSXpd6Vz5pNtK4CcRcWHpPJMVEY8C3yddmfCKwnEglQ3AVRHxr1HbfpSXixg8S0mXEA7CycsRx5ByfSwizouI9RHxcERcBBwP7Aac0GYAF3hLcnlfRLqu+qyIOKlsokkZOZa3b9EUyXzS9emHjJ7IQZoTMC9/vqFoyvGNTDaaWTRFsjYvH+2ybWRd6xNQJkPpbPX7SId3zi8cp9PBefnLLttW5+XebQbwScwWdEyUOBxYGRHHF440Wc/Py0G4WuYx4Nwxtr0LmEZ6UT/et0STNzLy7vVfpNoaq/Nyry7bRtbd3Z8ojR0A7ApcFhHrS4fpMD0vu10qOLKu1ZOtLvAekzSFNBPv3aTiWV400BjytP61MerKCEm7kK4Thv9dy1pMRPwd6DoNXdKBwOyIKD5NXdKLgTu67M8jSZfq3QVcXyJbp4i4XdJq4PWS9o88K1jSNOBT+W4/KBRvLEfn5SAdPoE0X+JQ4GRJvx753ucB3CfyfbqNznumtqn0y4B98qeLgZeQSuaBvO6bEXFtiWwjJH0a+CSwATiLdOJltBVR+FJCSStIv5ZeSzoZuIk0lfpQ0sjiixFxSql8TeS3VpgdgzH1ewVwJOlKiZER7GLgNcA/SW+hUPS5OULSnqTymU0q6/tJs5xfClwNLIkBmdCTr1W/n/Sb2M6DkgvS+wmRfijvQXrrjJ+SrkY5kDSV/g/AqyOivd8OS5/JneRZ31X8/yVvo29LK8gYwG4DkHMf0ojmVtKknSdIkzkuI5VN8e93g//DOgbkKhTS8dBLgT+TZopuAu4Evg7sXjpfl7y7k87RrM9Z1wKnAtNLZxuVc1l+zZxROssY+eaSroi7Le/HfwO3AJ8Dnt3241c1Ajczs//xVShmZpVygZuZVcoFbmZWKRe4mVmlXOBmZpVygZuZVcoFbmZWKRe4mVmlXOBmZpVygZuZVeo/7Fs5nUpZS1wAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "plt.xticks(range(0,8), range(1,9), fontsize=19)\n",
    "plt.yticks(range(0,2), range(1,3), fontsize=19)\n",
    "\n",
    "plt.imshow(C_S2)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Reconstruction loss: 7.446640570183208\n",
      "Relative reconstruction loss: 0.8325597259468726\n",
      "\n",
      "Reconstruction loss: 7.895194077231405\n",
      "Relative reconstruction loss: 0.8827095326121577\n",
      "\n",
      "Reconstruction loss: 8.215136024200474\n",
      "Relative reconstruction loss: 0.9184801297259811\n",
      "\n"
     ]
    }
   ],
   "source": [
    "X_np = X.numpy()\n",
    "\n",
    "approx = outer_product_np(X_1, X_2, X_3)\n",
    "print(\"Reconstruction loss:\", np.linalg.norm(np.ndarray.flatten(X_np-approx), 2))\n",
    "print(\"Relative reconstruction loss:\", np.linalg.norm(np.ndarray.flatten(X_np-approx), 2)  / np.linalg.norm(np.ndarray.flatten(X_np), 2))\n",
    "print()\n",
    "\n",
    "approx = outer_product_np(A_A1 @ A_S1, B_A1 @ B_S1, C_A1 @ C_S1)\n",
    "print(\"Reconstruction loss:\", np.linalg.norm(np.ndarray.flatten(X_np-approx), 2))\n",
    "print(\"Relative reconstruction loss:\", np.linalg.norm(np.ndarray.flatten(X_np-approx), 2)  / np.linalg.norm(np.ndarray.flatten(X_np), 2))\n",
    "print()\n",
    "\n",
    "approx = outer_product_np(A_A1 @ A_A2 @ A_S2, B_A1 @ B_A2 @ B_S2, C_A1 @ C_A2 @ C_S2)\n",
    "print(\"Reconstruction loss:\", np.linalg.norm(np.ndarray.flatten(X_np-approx), 2))\n",
    "print(\"Relative reconstruction loss:\", np.linalg.norm(np.ndarray.flatten(X_np-approx), 2)  / np.linalg.norm(np.ndarray.flatten(X_np), 2))\n",
    "print()\n",
    "The "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}