Standard_HNCPD_twitter.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook contains the Twitter data experiment results where we applied Neural NCPD to the Twitter dataset. 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": 1,
   "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 weights_H, outer_product, outer_product_np\n",
    "#\n",
    "import torch.nn as nn\n",
    "from torch.autograd import Variable\n",
    "from writer import Writer\n",
    "\n",
    "import tensorly as tl\n",
    "from tensorly import unfold as tl_unfold\n",
    "from tensorly.decomposition import parafac, non_negative_parafac\n",
    "\n",
    "from sklearn.decomposition import NMF"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "torch.set_default_tensor_type(torch.DoubleTensor)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load in Twitter Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "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",
    "results = np.load(\"./data/results_twitter_ncpd.npy\", allow_pickle=True).item()\n",
    "\n",
    "# Here we load in the initial factor matrices to make sure they are the same accross comparisons\n",
    "X_1 = results['X_1']\n",
    "X_2 = results['X_2']\n",
    "X_3 = results['X_3']\n",
    "\n",
    "X = X.reshape((8, 10, -1))\n",
    "\n",
    "#X = torch.from_numpy(X)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Run Standard  HNCPD"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def standard_HNCPD(X, r):\n",
    "    \"\"\"\n",
    "    Given a 3-mode tensor X and a list of rank r, computes the Standard HCNPD factorization.\n",
    "    \"\"\"\n",
    "    np.random.seed(0)\n",
    "    #factors = non_negative_parafac(X, r[0])[1]\n",
    "    factors = [X_1, X_2, X_3] # This is to make sure all methods have the same initial factorization\n",
    "    \n",
    "    As = [[], [], []]\n",
    "    Ss = [[], [], []]\n",
    "    \n",
    "    for i in range(1, len(r)):\n",
    "        \n",
    "        for k in range(3):\n",
    "            \n",
    "            model = NMF(n_components=r[i], init='random', random_state=0)\n",
    "            if(i==1):\n",
    "                A = model.fit_transform(factors[k])\n",
    "            else:\n",
    "                A = model.fit_transform(Ss[k][-1])    \n",
    "            S = model.components_\n",
    "\n",
    "            As[k].append(A)\n",
    "            Ss[k].append(S)\n",
    "            \n",
    "    return factors, As, Ss\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "Xs, As, Ss = standard_HNCPD(X,[8,4,2])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Display Top Keywords for all Layers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "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",
      "johnkasich  hillary          tedcruz     marcorubio  \n",
      "kasich      trump            cruz        teammarco   \n",
      "ohio        people           ted         vote        \n",
      "john        realdonaldtrump  internet    flsen       \n",
      "gov         crooked          choosecruz  florida     \n",
      "\n",
      "RANK 2 KEYWORDS:\n",
      "------------------\n",
      "Topic 1   Topic 2     \n",
      "-------   -------     \n",
      "tedcruz   marcorubio  \n",
      "cruz      teammarco   \n",
      "ted       vote        \n",
      "internet  florida     \n",
      "obama     flsen       \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 = Xs[2]\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 = As[2][0]\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(As[2][0], As[2][1])\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": 7,
   "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(Xs[0], cmap=plt.get_cmap(color))\n",
    "axs[1].imshow(Xs[1], cmap=plt.get_cmap(color))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "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": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeUAAAFwCAYAAACGm2OSAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABEvUlEQVR4nO3dd7hcVdn+8e+dkAQSOgEBg1QpCkgJHTR0QaRIJ5QoGAEVEURFXjXIKy8CFkQE80NykCJdqoIIBASkJCgiEhQwBEIH6ZD6/P5Ya5JhMu20Kefcn+s61z6z99prPxOY88zaexVFBGZmZtZ8A5odgJmZmSVOymZmZi3CSdnMzKxFOCmbmZm1CCdlMzOzFuGkbGZm1iIWanYA1nmDNSQW1rBmh1HVQmsNbHYINc15sg2+kw5o/RjfnPHCKxGxbLPjaHULLTIsBi+2dFNjmDOkqZeHgc0fgrvCoq83OwSefvSdip8ZJ+U2tLCGsflCOzc7jKqWvGDxZodQ01v7L9LsEGqKoQs3O4Sabnn8R083O4Z2MHixpfnofsc1NYa3Vpvb1OvPWXJ2U68P8L2tr292CBy+1r0VPzOt/zXczMysn3BSNjMzaxFOymZmZi3CSdnMzKxFOCmbmZm1CCdlMzOzFuGkbGZm1iKclM2sLUmaKilq/OzZxbpXyedP7dmozarz5CFm1u5uAV6ocGxaIwMx6y4nZTNrd6dFxMRmB2HWE3z72szMrEU0JSlL2lXSRZKekPS2pBmSnpV0k6QjJS1WVHZcfrYzrqSOUXn/xEbHb2btR9LOkq6X9KKkmZKel/RbSevVOG8hSd+W9Jik9/P5F0r6SKNit/6joUlZ0nI5id4EHAzMBP4I/A6YCuwAnAs8JWnlRsaW4wtJzV/GxMx6lKSzgJuBXYAngWuB54EDgAck7Vrl9MuBk0nPp68FZgCHAg9KWqv3orb+qGHPlCUtCdwDrAH8BTgyIv5eUmYx4CjgJGApoNrqMw8A6wDv9ka8ZtY3SDoSOAZ4FNgnIqYUHdsTuBK4RNJqEfHfktNXBhYBNoyIf+ZzBgO/JjUsLgI27fU3Yf1GI1vKvyAl5AeA7UoTMkBEvBURpwMbAy9Wqywi3o2IKRHh3pVm/dsdFYZDdUgaCHwvl9uvOCEDRMS1wK+AJUlJtpxTCgk5nzMT+CrwJrCJpK0qBSZprKRJkibNfu+dLr9B6z8akpQlrQ4cmF8eGRHvVysfEU9ExPM16iz7TLl4fKGSoyX9TdK7kv4r6TpJ65acM674tnXph7ukrCQdImliru99SU9KOkfSShVinVePpP0l/SU/S39L0m2Stq72Xs2sqluAC8v83A1sAKwAPFqcWEvcmbdbVDh+cemOiHgduCG/HFUpsIgYHxEjI2LkQosMq/omzKBxt693I30BeCQi/tqgawJ0APsDdwH/BjYBdgdGSdowIp7K5f5G+hAfll9fWK4ySSJ9QA8CZgETgddIt6+OBg6Q9OmIeLDC+T8g3Zq/m/RcfX1gO2BrSaMi4i/deK9m/VXFIVGS9s2/fryO/iLLltn3ek7A5UzN2xG1AjSrV6OS8sZ5WzZZ9ZKVgW2Aj0fEkwCShgDXALsCJwJfhHm3sK6VdFh+PaZCnUeREvKLwPYR8WiudyDwU9ItrSslrRURM8qc/2Vg04iYnM8bAJyX4/gBsGP33rKZlRiYt9OBP9UoO6XGcbNe16ikXPgG+lKDrldwTCEhA0TEDEknk5Ly9l2o7/i8/W4hIed650j6BrAH6cvAPsAlZc7/fiEh5/PmSvouKSlvI2lQRMwqd2FJY4GxAAsztAuhm/VLz+Tt81W+bFezpKQlIuKNMsdWydvpXQnMrJy+PHnIbNIQiFKFb8MrdqYySSOA1YC5pB6XH5A7fxQS8agK1dxY5rwXgf8CQ4BlKl2/+NnUIA3pTOhm/dkDwKvAhpLW6GIdo0t3SFqC9FgO0mMssx7RqKT8ct4u16DrQfpmPLt0Z0S8mX/tbGb7cFG9lTqqPVVStlSlnuKFmBbuZExmVkW+83QK6Tb2tZIWGL4kabCk3SWtXaGa70lap6j8IOAsYAlgckTc3QuhWz/VqNvXk4FDSB2tGmVuL9Xb5clFIqK3YjKzCiLirDwZ0deB+yX9nTSByEzSF+gNgWGkiUVKnytPI/39+puk24E3gC2BlYBXSJOImPWYRrWUbyIlyfUkbdiga/a0wnOjFXOHsXJWKylrZi0gIo4DPgVcRpqY6DPAp4HhpMdKo4E/lzsV2I/U2l4N2JN0R+tiYJMqw6zMuqQhLeWIeELS5aSxyudK+lSF3snAvHHN79Yaq9wLZgGDJC1Ueus7Ip6V9BTpg3kwaUafefItrcKzp4kNiNWsX4uIVTpZ/i7S8Mh6yk4FVLTrf/OPWa9qZEevr5CeuW4G3F5uEnhJwyQdR7pd9KEGxlZQaOGuU+H4T/L2lOLnT3lI1OnAR0hTg17VaxGamVmf1bC5ryPitTxz1RXA1sDfJf2T9Ayn8GxnU1IHrBdJk3I02u9Iz51uy8+P3s6xH5GP/xLYitTifzjPJlaYPGQ1Ui/qfavdBTAzM6ukYUkZIN+O3kbSbqTEtgXpuc5CpB7afwKuAy6NiGZMFHsS6RnSXsDngEF5/xEAERGSRgN/II0t3ow0Wf1zpNWt/i8inimt1MzMrB4NTcoFEXEjZcbsVig7DhhXZv9EPvjMp7B/arn9JWXKHo+I90gThBxf7nguE6RxyguMVe7s9YqOr1JvXWZm1nf15clDzMzM2oqTspmZWYtwUjYzM2sRTXmmbGbW7wQMnNXlCQF7xNyFm3t9Zje/Hbjfos82OwQOr3Ks+f9CZmZmBjgpm5mZtQwnZTMzsxbhpGxmZtYinJTNzMxahJOymZlZi3BSNjMzaxEep9yWBAMHNjuIqu5/ZI1mh1DTx2Y/3ewQanp18+WaHUJtjzf2cpKmAivnlz+KiG9XKXsx89c5vzMiRvVudGbd45aymbWzQ/J65guQtDhptTeztuGkbGbtahKwIrBjheMHkJZWfbBhEZl1k5OymbWrjrwdU+H4GGAOnVhm1azZnJTNrF3dDzwG7CFpyeIDktYCtgBuAZ4vd7KkHSSdI+lhSa9KmiHpaUkXSlqnwjkdkkLSGEnrS7pS0guS5kg6tkffnfVLTspm1s4mAAsDB5bsH1N0vJLzSGsDzAbuAn4PzAQOBSZJ2rrKuVsBDwAbAROBm4F3Oxe62YKclM2snV1EukU9prAjd/w6FHgNuL7Kud8Alo+IjSNir4jYC1gTOBIYCoyXpArnHgGcCXw0Ig6IiM9ExPhuvxvr95yUzaxtRcQLpFbqpkW3nHcidQC7NCJmVjn32oh4vWRfRMSvgHuBdYCPVTh9CvD9iJjbzbdg9gEep2xm7a4D+Ayptfwt5reaO2qdKGlEPndtYHGgMLxq+bxdE3i0zKnXRcScOuofC4wFGLToUrWKmzkpm1nbux54lTRm+QxgD+CRiJhc7SRJJwPfofrfwcUr7K9r5pl8S3s8wNBlV4p6zrH+zbevzayt5VvUlwIrkDp2DaF6By8k7Q18D3gP+CKwOjA0IhQRAn5bKFqhivd6IHSzBTgpm1lf0JG3u5F6U19So/y+efudiDg/Ip6KiOJE2/rzxFqf5KRsZm0vIh4C7iHdxr4yIl6qccrSeftM6YHcYWzDno3QrD5tlZQlTc0D90dVOL6epOdzmaskDe7idcbkOjq6Ea6ZNVBEbB0RwyPioDqKT8nbLxb/nZC0HHAh7m9jTdJWSbkaSZsBd5J6TU4A9q82HMLM+rWfAW+Qel4/kWfmuhF4ElgUuLZ5oVl/1ieSsqTtgD8BS5E+bIfXM1yhit+Rxiie2P3ozKzVRMRTpFvUl5E6c32W9JkfT5qe843mRWf9WdvfopG0O3AFqcflyRExrrt1RsQb+ENp1pIiYpVOlr+KMr2oI+I/LDg9Z8EYyix0ERFl95v1lLZuKUsaDVwNDAaOLZeQJW0m6QxJkyS9KGmmpOfyM+fNK9Rb9pmypFF5/0RJgySdJGmKpPclvSTpYkkfqRLvSpLOkvS4pPckvSnpnny9SkMvzMysn2jbpCzpKNK8twK+EBFnVSj6Q+DrwCDSBPKFiQb2Bu6WtG+F86oZBPwB+DbwRP59LjA617lkmXi3BR4BjiH9u99MWuVmfdIz8Au7EIeZmfUh7Xr7+ljSrD0zgQMj4poqZc8ERkfEi8U7JX2W1Mo+T9JNEdGZFV62JC2wvnph6IWkJYDbSavGfJn0ZaBwrRXytRYl3fr6TUREPrYS6YvCIZJuj4iOTsRhZmZ9SLu2lPfI25/VSMhExM2lCTnvvwG4kjRecdtOXj9IrfN5YyHzc+gf5Zfbl5Q/ltQJ7ccRcWEhIefzniHNKATw1UoXlDQ234KfNCve72S4ZmbWDtq1pXwn8CngeEmTI+KKaoUlDSfN9LMusCTz3/e6ebsmcFMnrj8tIh4ps78w9nHFkv275u2VFeqbDLwNbCBp4YgFs27xHLqLD1jGc+iamfVB7ZqUx5HGF34DuEQSlRKzpC8BPyGtj1pJpUnnK5lWYf+bebtwyf7V8vbBOvpzLQNM72Q8ZmbWB7RrUiYiTsgJrmJilrQJcC5pLtwTgBuAZ4F3IyIknUoai9zZns+dXUO1sBzc5UCte88zOlm3mZn1EW2blKFsYo6IKL5FvDcp4f48Is4sU0WjJp1/Jl/rlIgotzarmfV1A2D2Is0d+bjQ8OYubjVn1sDahXrZXiM2bXYIVL7Z2r4dveaJiBNIPawXAi4tGeJUbdL5ZYEdez9CIA2Zgvkr05iZmS2g7ZMylE3M++RDhY5Xh0patFBe0mLABaROX41wBul583ckfVnSAncoJH1c0ucaFI+ZmbWgPpGUYYHE/NucmCeQWskbAU9JukbS74CpwEhSYm5EbM8AewJvAb8Apkm6VdIlkm6SNA34B7BfI+IxM7PW1GeSMiyYmEnjhUeShhK9TeqxPRK4hpSoF7it3Yux3QF8HDgVeAnYnPTM++PAU6QOZyc1Kh4zM2s9bdXRq56J6HNiPqFk95cqFB+Xf0rr6AA6yuyfSJWe2hExtcbxF0iJ18nXzMwW0KdaymZmZu3MSdnMzKxFOCmbmZm1CCdlM+sTJA2QNC2vef6ypEHNjsmss5yUzayv2BFYKf8+HNi9ibGYdYmTspn1FV/I2+klr83ahpOymbU9SUuT1lkP4ABgDrCzpNJlVM1ampOymfUFo4EhwMSIuBv4I2l1tsPKFc7PnSP/PlbSXyW9K+nVPPPfunWcd7ik+yW9mfcv2RtvzPoXJ2Uz6wsKt6o78nZC3n6+2kmSfkpa3vUN4DrgFWAv4H5JW1c572zSTIEzgBuByaRWulm3tNWMXmZmpSRtCGxAmlv+qrz7euA14KOStomIP1c4fSywbUTclesSaSrcb5MWt1kzIsqtgX4IsEVEPNBz78TMLWUza3+FVvIVEfEuQETMAC4pOV7OuYWEnM8L4H9I89GvRJqfvpzTnZCtN7il3IYkoYVa+z/dR25sdgS1zVxjhWaHUNNSv53c7BBamqQhwEH55YSSwxOArwL7SvpqRLxdpoqLS3dExBxJvyXNUT+K+cm92DV1xjeW1Bpn0GJL1XOK9XNuKZtZO9sTWBr4d0TcU3wgIv4KPAwMA/avcP5/KuyfmrcjKhx/up7gImJ8RIyMiJELLTKsnlOsn2vt5paZWXWFW9NLSLq7zPHlisr9uqcuGhHv9VRdZsWclM2sLUlaCdghv1yO+Qm4nC0lrRURj5fsX4XUmi61St5OL3PMrNf49rWZtasxpL9ht0eEKv0AV+Ty5Tp8jS7dIWkgaQISgIm9ELdZRU7KZtZ28tClMfnlRTWKF44fkhNusaOLxyPnek8GVie1kq/ufrRm9fPtazNrR6OA1YD3qJ04bwZeBlYAdgVuKDr2/4A7Jd0FPA9sBKyV6x3tZ8fWaG4pm1k7KtyKvjYi3qpWMCJmA5eVnFdwHGnY1NKkntzLAdcCm0XEnT0VrFm93FI2s7YTEYeQZtWqt/wxwDFl9gfwy/xTTz2q95pmXeGWspmZWYtwUjYzM2sRTspmZmYtolNJWdIASQdJuk7SdEkzJL0maZKkUySVHbwvaVxeb3RMyf4xhfVJJb1XbT1SSR8tKhuSRnUm9ir1zlsftWT/xJ68jpm1jqIxzGYtpe6kLGkE8ABpcvbdgGmkSdnvBVYlrazypKT9uhjLwswfsF/OmC7Wa2Zm1hbqSsqSlgb+DGxMmuFmjYjYIiIOjIjdgOWBE4GhwGWSKi13VslfgTlUSLySBgCHAv8Fnuhk3WZmZm2h3pbyOaS5YB8EdomID6ysEhGzIuI00pg/ARdIGt6JOJ4DbgU2k7R2meM7kFZruQyY0Yl6zczM2kbNccqSVgcKt6SPjoj3qxT/OXA4sB7wFWBcJ2KZAHya1Fr+dsmxzxeV+WSFOJclrau6C2lGnhVICfwx0jR750XEnE7EU5WknYEvA5sBSwGvku4inBoRjxSV2wa4C5gSEetUqGs48CwQwIiIeLWn4jSz1jBgFgx9cW5TY1jxQ6809fpTnm7+GuYvH7lFs0OAc6+qeKielvJuudyjETGpWsE8EP83+eXu9caXXUe6Pf2B+Wlz56898/UfrHL+zsDPgI+T1kj9HfAQsAHwC+DqPK9tt0k6izR13y7Ak6QZgJ4nPRN/QNKuhbIR8WfSKjRrS9quQpVHAEOAy5yQzcz6r3qS8sZ5+0CddRYS5yckLQQQEeNyb8eOSidFxAzgt8CKwE5Fhw4gdQKreG42Gdg8IlaKiO3y8+5tSZ3Q/gbswfwWf5dJOpI0M9CjwHoRsWVE7BcRGwF7ke4+XCJpqaLTzs7bo8vUNwD4Un55TnfjMzOz9lVPUl42b1+ss85CuQGk+WQ7Y0Lejina93lgNjVWgomIxyLi/jL7nwe+mV/u08l4PiC34L+XX+4XEVNKrnUt8CtgSeDgokOXAq8Be0hasaTaz5Cf19e6E2FmZn1bb0we0uVbxDkp/YOUvJaStA6wKXBzRNT8UiBpIUk7Sfq+pHMlTZDUARyZi6zZ1diyDUjPqh+NiH9WKFOYxH7eg4u80sz5pFb02JLyhdZz1VaypLF5PPikmVUf65uZWbuqZ0GKQs+AD9VZZ2ECkbmk1mFndQBnAgeSWpAwvwVdkaQ1Sc92y3amyhbvQjzFVsvbj5ebcKTEsiWvzwGOB74o6X8jYnbuRLczqZPY5dUqi4jxwHiAJQYOr3VtMzNrQ/Uk5cmkW7Gb11nnpnn7cF4yrbMuBk4jLbG2IulLwQ1Vz0iuIiXk64HTSb2u34iIOTlhP043WvFZoQPadOBPNcqW3tqeJul60nPnPXO8R+WYLqjRq93MzPqBepLyjcCPgXUkbVKtB3Tu3XxofllPIl1ARLwo6Q/AZ/Oun0fErGrn5LHN6wEvAZ8rM/Rpja7EUsYzeft8RIzpwvlnk5Ly0ZJuJD0vnwuc2zPhmZlZO6v5TDkiniC16gDOkbRwleLHAOsCb9G9nsTnk27pvkodt66Z36HsuQpjkUd3I5ZiD+SYNpTU6UQfEXeQnplvC5xMivsPpZOxmJlZ/1RvR68vk1qJmwC/l7RK8UFJgyR9C/gJaQKMIyLipa4GFRHXR8Tw/PO3Ok75N6nFua6kD0wuIunzpOfT3ZZb7KeQbmNfK2nT0jKSBkvavcLMZJDGTMP8HuF1La5u1l9Jmlq0EM1pNcpeXFR2YoNCNOsxdSXliHgF2Jo03ndb4AlJ90i6VNINwAuk58DvAaMj4opeirdSfC+TkttCwB2Sbs+xPQJckGPrqWudBfyUNEnJ/ZIelnSNpMsk/ZnUue065ndSK3URaZIUgKdIk5CYWX0+MLlQMUmLA59rcDxmParuIVERMQ0YCRwC/IE0Kcc+pGQ9FfghaaGK3/Z8mHX5Gmm40cOkzma7kMZM70LutdxTIuI44FOkubiXIo01/jQwnPQMfjRpAY9y575LWlkL4NyIaO68e2btYxKp8+eOFY4fACzC/AmMzNpOPR295snPay/OP92WZ/jq6OQ561bYPxf4f/mnnLI9ryutqRoRo2rEcRdpTutOyWtO70i6q3BBZ88368c6SA2DMZS/wzSGtNrcRaRHbWZtpzcmD7HqTgIGAxdGRFfGcZv1V/eThjrukefEn0fSWqQJe24hzUO/AEk7SDonP3J6VdIMSU9LujBPVFRa/rb8bLriOu+SfpzLnN6dN2ZW4KTcAJK2lHS+pDtIPdTfIHUYM7POmUCaC7+08+aYouOVnEdaxW426S7X74GZpGGckyRtXVK+4pz1AJIWwcMarYc5KTfGmqQ/BpsCdwOfjojnmhuSWVu6iHSLekxhR+74dSipk+X1Vc79BrB8RGwcEXtFxF6kz+aRwFBgfMlKcjcATwPbSCr32OxAUp8SD2u0HuOk3AAR0ZFXyRoWEdtExH3NjsmsHUXEC6TnyZsW3XLeidQB7NKImFnl3Gsj4vWSfRERvyJ1vlwH+FjRsTnMH7JYrrVc2FdxWGPxnPWzZrxT9b2ZgZOymbWfjrwdU7LtoAZJIyR9SdJPJf1aUkdetGb5XKR00ZrzSZ0yD5a0WFE9m5GWta06rDEixkfEyIgYOWjIsFrhmXWu97WZWQu4njSz3iGSziCtlf5IREyudpKkk4HvUP3v3gcWrYmI1yRdAhxBGg5a2nL2sEbrUW4pm1lbybeoLyUtozoBGEKN6Xgl7U1aC/094IvA6sDQ/FhJQGF+hXJDJAsdvo7KdS0D7Ae8j4c1Wg9zUjazdtSRt7uRelNfUqP8vnn7nYg4PyKeyuucF1Scyz4i/k7qrV2YxvdwUg/wyzys0Xqak7KZtZ2IeAi4h3Qb+8o65tovLFrzTOmB3GFswxrnF1rLXyH11obuLbpjVpafKZtZW4qI0nHF1UwhzaT3RUl/KPTSzjPsXUjtv4XXkhJ6ocX9YERM6lzEZrW5pWxm/cHPSJP2fIa0oM6VeU3zJ4FFSUm3ooiYzQcnCHEr2XqFk7KZ9XkR8RTpFvVlpM5cnyWNSx5Pmp7zjTqquTVvXwUu74UwzXz7uh3F3LnMfe/9ZodRVbTB/1kD7n2k2SHUpGFDmx1CbRWn6+gZEbFKJ8tfRZle1HnWrUprq4+haJawCkbn7a8jorU/gNa23FI2M6tB0kqkoVQzgV80ORzrw9qgPWNm1hySTgNGkDqJDQPOiIgFenCb9RQnZTOzyg4APkJaDvJHwHebG471dU7KZmYVdPZ5tll3+ZmymZlZi3BSNjMzaxFOymZmZi3Cz5TNzBpE0dzrr7boq029/mOzV2zq9QFmLlluIbDW4ZaymZlZi3BSNjMzaxFOymZmZi2i20lZ0tqSzpP0uKR3Jb0naZqkeyX9WNKOPRFob5A0TlJIGtfsWMzMzLrV0UvS/sBvgMHAdGAi8F9gWWAj0uorn2L+6ipmZmZWQZeTsqTlgQtICfnrwNkRMafo+ABg6/xjZmZmNXSnpbwbMBT4S0T8rPRgRMwF7so/ZmZmVkN3nikvl7cvdeYkSZtJOkPSJEkvSpop6TlJV0navMI58579SvqQpF9JelbSDEn/kXSapIUrnDtI0jck/VPS+5JekHSRpJXrjPWyfK2Zkl6WdL2ksq3/HGPk3w+XdL+kN/P+JfP+FSX9QtITOZ538zP4myWNrfsf0swWIGli/ryNaXYsZl3RnaQ8LW+3l7RuJ877Iel29yDgAeB64FVgb+BuSftWOXclYDKplf4X0jPs5YBvAVeUFs630K8BzgBWBW4H7gS2z/WsWulCko7P19gPeAG4DngC+Axwp6QvVjn3bGA8MAO4MV8rJK2Qf/8y6S7FzcANpH/LzYHjqrx3MzPr47pz+/o64DlgReCvkv5ISngPAQ9GxBsVzjsTGB0RLxbvlPRZ4GrgPEk3RcS7Zc79AnA+8OWImJnPW4eU3D8raauIuKeo/JdJCXw6MCoinsjnLAxcDBxaLkBJu+Q4nwM+FxH3Fx3bCvg9cI6kOyPiX2WqOATYIiIeKKn3a8DywK+AoyIiio4NATYrF4+ZmfUPXW4pR8RbwA7AJFJy35W03uitwGuS7sm9s0vPu7k0Ief9NwBXAksD21a47DPAMYWEnM97DLgov9y+pPyxefs/hYScz3kfOBp4r8J1xuXtEcUJOZ97D3AKqaX/pQrnn16akLMP5e3NxQk51zsjIvz83cysH+vWOOWIeCwiNgG2Ak4FbiMNiRoAbAlcJqmj9DxJwyWNkXSmpPMldeRyhdvga1a45O0RUS6RTsnbeROrShoBrAbMBS4tE/tLwB/LxQZsCrxZ7nh2Z95uUeH4NRX2FxL1jyTtKWlYhXILkDQ2P4efNIsZ9Z5m1u/lvy8VnzNXmq+gO31ZzLqqRxakiIh7gXth3nPczYHvAzsBh+Xb0Vfm418CfkLquV3J4hX2T6uw/828Lf6AjMjb54pb1iWmltlXeM68ODBbqjp5+bIV9j9dYf9FpH+Tg4DfAXMk/YPUQ/2y/O9YVkSMJz2nZnEt3eRp7c36lUJfFpH+zi1OGur5LeBjwO7NC836mh5fJSoPhbpX0q6kluFGwJ7AlZI2Ac4FZgMnkDo5PQu8GxEh6VTgRNL//OXM7el4yxiYt28A19Yo+0q5nRVa84V/m9GS/o/0rHur/PNV4KuSLoiIw7sStJn1ms72ZTHrsl5bujEi5ki6nZSUCy3KvUkJ9+cRcWaZ09bowRCm5+2KkgZXaC2vUmbfM3k7KyLG9GA880TEP4B/wLw7C7uSbrF/QdLlEVHptrmZNV7ZviySLgKOIvVlcVK2HtHlZ8qqcV83+0jePpu3S+ftM6UFJS0L9Ng82RHxDPAf0ns8oN7rRcR04BFguKRRPRVPJRExNyJuJPVmB/hEb1/TzDql7r4spT7QF2TGO70TnfUp3enodbSkCZI2LT0gaaE8jnefvOvyvC38T3yopEWLyi9GmrJzyW7EU87P8/Z/Ja1WdL0hwDlUfq793by9WNJOpQclDZS0nSpMdlKJpEMlbVRm/zLM7zRW6Xm0mTVHZ/qyfEBEjI+IkRExctCQuvt1Wj/WndvXg4AxwBhJLwB/A14jtYbXZ/63x9Mj4pb8+wTSMKWNgKck3U26nf1JYCYpMX+hGzGVOpvUsWoX4NF8O/1tUieNhUmLaSwwVjkirsuTh5wO3CLpX8Dj+dzlgQ1JXyCOAu7rRDyfAy6UNJ307/U6sAywDTAM+DOpA5iZNU6txkkj+rKYAd1rKf8a2IvU4nwGWA/Yl5Rg3wYuBLaJiG8VToiI/wIjSb2I3ybNjjWSNIRoI8rc1u6OvEDGHsC3ST2tdyCNgb4rX/c/Vc79CbAx6X0OJN3q/iypV/ddwBcpM4tYDT8GziJNSjKS9O+1PmnClcOBHSNiVifrNLPqCs+CF61wvOaUu2aN0uWWcp485Fpq91AuPe8lKk+6MY75E3cUn1N2f9HxDqCjwrFZpElNflTv9YrOfRg4otLxMuWrPmePiD+TWsNm1jiFTp9rlx6QtAgwqqHRmFXRrclDzMzawG15e4iktQo7c0I+l/kdUs2azknZzPqiec+BI+Ju0sIwiwMP5RXZbiQ9vtqJ1NfFrCU4KZtZX7JI3paOP9oXOI201Ox2pD4sN+Vtpd7VZg3Xa5OHmJk1Up47oTD0cWrxsbwIzYn5p9Q4ergvi1lXuaVsZn3FYcBw4GXg4SbHYtYlbimbWduSNJS0PvnqzJ+A57sRMbt5UZl1nZOymbWzwcDBwFuk+afPKqxIZ9aOnJTNrG1FxOtUXlXOrO34mbKZmVmLcFI2MzNrEU7KZmZmLcLPlNvV3DnNjqCqWUNb//ve3L1GNjuEmhZ/5JVmh1Dbm7WLGBDBwBnNXXDq4GXuber1fz9jg6ZeH2CJp1p70a/W/8tpZmbWTzgpm5mZtQgnZTMzsxbhpGxmZtYinJTNzMxahJOymZlZi3BSNjMzaxFOymZmZi3CSdnM+g1JEyWFpDE9UNeYXFdH9yMzS5yUzczMWoSTspmZWYtwUjYzM2sR/SIpS9pV0kWSnpD0tqQZkp6VdJOkIyUt1uwYzaw5JHVUe84saVw+Pq6xkVl/1KdXiZK0HHAF8Km86zHgj8BMYASwA7ArcIqkkRHxdFMCNTMzow8nZUlLAvcAawB/AY6MiL+XlFkMOAo4CVgKcFI2M7Om6bNJGfgFKSE/AGwXEe+XFoiIt4DTJV0DvNPg+MzMzD6gTz5TlrQ6cGB+eWS5hFwsIp6IiOfzufOeL0laX9KVkl6QNEfSsaVlKlx/gWdQkkblfbV+Vun+v4CZtQJJYyVNkjRp9gx/77fa+mpLeTfSF45HIuKvXaxjK+A8YDowEVgMeLcbMb0AXFjh2PLAzvn3Od24hpm1kIgYD4wHWHSpEdHkcKwN9NWkvHHePtiNOo4Afgh8LyLmdjegiJgCjCndL2lR4K788icR8Ux3r2VmPapP3lG01tRX/2dbNm9f6kYdU4Dv90RCrkTSQOByYEPgauCE3rqWmVU0M28XrXB85UYFYtZXW8o94bqI6O1byb8gDcm6Dzi42hcASWOBsQALM7SXwzLrV6bn7dqlByQtAoxqaDTWr/XVlvLLebtcN+ro1eFRkr4JHAk8CexeR2e08RExMiJGDmJIb4Zm1t/clreHSFqrsDMn5HOBjzQlKuuX+mpSnpy3m3Sjjve6cW7Vf1dJ+wKnAa8Cu0bEy9XKm1mPm3dXKiLuBm4EFgceknSzpBuB/wA7AROaE6L1R301Kd9E+tCtJ2nDXqi/y8+gJG0J/CbXsUdE/KuHYzOzyhbJ29LxSYUvyi8B2wEbkf6ObARMa1h01u/1yaQcEU+QOlABnCup6v1eSatLWqETl+jSMyhJawDXAUOAQyPink5c08y6QZKA1fLLqcXHIuL9iDgxIlaNiMERsWJEHB4RL0TEuIhQRIwrOacj7x/TkDdg/UKfTMrZV4CngM2A2yWtV1pA0jBJx5Fud3+oE3V3+hmUpGWAPwDDgW9HxBWduJ6Zdd9hpM/fy8DDTY7FrKw+2/s6Il6TtDVpQYqtgb9L+idpqNNM4MPApqRW64vAa52o++78zGk30jOoPwOzgZGk2+YTgM+XnHY0adrPt4GPSeqoUP03IuKVemMxs8okDQV+BawObJF3fzciZjcvKrPK+mxSBshTZ24jaTfStJtbAJ8mve+XgT+RbidfGhGdnQNvX+D7wAGkZ1CvkJ5BnUTqVV1qYN4uSvrGXsm4XJeZdd9g4GDgLdICNWdFxJXNDcmssj6dlAsi4kZS78p6yo6hzMxbZcq9D5yYf0qNyz/F5RfYZ2a9KyJeB9TsOMzq1ZefKZuZmbUVJ2UzM7MW4aRsZmbWIvrFM2Uzs2ZTwIDZzV298eJXt2zq9Yc+M7B2oV5244/PbHYIDL+88jG3lM3MzFqEk7KZmVmLcFI2MzNrEU7KZmZmLcJJ2czMrEU4KZuZmbUIJ2UzM7MW4XHKbWjOGkN47aw1mx1GVcsd9WyzQ6hp9rPTaxdqtiWXbHYEZtZAbimbWduRNEDSNEkh6WVJg5odk1lPcFI2s3a0I7BS/n04sHsTYzHrMU7KZtaOvpC300tem7U1J2UzayuSlgb2AAI4AJgD7CxpxaYGZtYDnJTNrN2MBoYAEyPibuCPwEDgsNKCklbJz52nVqosHy+7UoSkDSRdJ+k1Se9ImizpC7XOM+sqJ2UzazeFW9UdeTshbz/fkxeRtB3wF9Lz6heB64E3gfGSzujJa5kVeEiUmbUNSRsCGwBvAVfl3dcDrwEflbRNRPy5B64zFLgYWBj4ATAuIiIf2xK4pbvXMCvHLWUzayeFVvIVEfEuQETMAC4pOd5d+wArAP8CTi4k5Hy9e4Ff1lOJpLGSJkmaNGvmOz0UmvVlTspm1hYkDQEOyi8nlBwuvN5X0qI9cLlP5e3lETG3zPFL66kkIsZHxMiIGDlo8LAeCMv6OidlM2sXewJLA/+OiHuKD0TEX4GHgWHA/j1wrQ/n7dMVjlfab9YtfqZsZu2icGt6CUl3lzm+XFG5X9dToaRaDZNKvavLtZ7Nus1J2cxanqSVgB3yy+WYn4DL2VLSWhHxODAz76t0S3vlCvufq3F8lSrXN+uylrh9LWlqHvM3qsLx9SQ9n8tcJWlwg+NapZv1dOR6xvRMZGb9zhjS36vbI0KVfoArcvlCq/plUmJeRtKyZerdtcL17srb/Sq0pg/s2tswq64lknI1kjYD7gSWJ3Xm2D8iZlY/y8z6CkkiJWWAi2oULxw/RNLAiJgFFIZInZzrKtS7NWm4UzlXksYmrw2cVHLeZsCXO/UmzOrU0kk5D97/E7AU8DPg8IiY08AQtgfWYf78umbWeKOA1YD3gKtrlL2Z1Dpegfmt4O+RWstHAY9KulLSA6Qv+2WHNkXEO8AhwAxS4n5U0qWSbgfuAc7PRWd18T2ZldWySVnS7sDvSc+CTo6IrxePFWyEiHgyIqbkb9tm1hyFW9HXRsRb1QpGxGzgsuLz8rji7YHbSCtLFZL1oRHx3Sp13QpsCdxASvJ7khoIRwM/zcVe6eR7MauqJZOypNGkb8SDgWMjYlyZMptJOiMPzH9R0kxJz+VnzptXqHegpCMl3SvpjXzOi5IekvTj0mdO1Z4pK9lP0h8kvZTrmi7pNklfrfLe1sjfuF+UNEPSFEnfqqMXqFm/FBGH5GfGB9UuDRFxTC6/V9G+uyNih4hYLCKGRcSmEXFJPlZ4Hl2urociYveIWCoihkbEhhExntR6B5jczbdn9gEt1/ta0lHAOaQhB1+IiI4KRX9I+mA8CjxAus20FrA3sKekAyPiypJzfk2atP494G7St9zhwOrAcaTnSC/XEePgXHZ30go19wHTgA8B6wLbAWeXOXUD4Kx83TtIPUi3AU4DRgAVk7mZNZak5YBFIuLpkv2bA4W5rzsaHZf1ba2WlI8lLck2EzgwIq6pUvZMYHREvFi8U9JnSa3s8yTdVJiKT9LKpIT8DLBJmfM2YP4wiFpOJyXkfwF7RMSUonoGAp+pcN7XgJOBHxRmCZL0SVKCPlrS6RHxTJ0xmFnvWh+4VdI/gP+Q/i6tBmyYj18UEbWecZt1SqvdMt0jb39WIyETETeXJta8/wZSK3ZpYNuiQ4VxjQ9VOO9vEfFSrQDzt+ejSC35zxUn5FzPnIi4vsLpD5Kej88tKn8XaXL7ASXxmllzTQHOJX02tyb9fVoZuB04lDJLRZp1V6u1lO8kzTl7vKTJEXFFtcKShgO7kW4ZL8n897Nu3q4J3JR/n0JaWeYzkr4DXFJ6W6pO25Gedd8TEY928tzfV+isNgXYBai4SLukscBYgMHLLt7Jy5pZZ0XEs6ROXWYN02pJeRzp1u83gEskUSkxS/oS8BNgaJX65mWviHgrL05+Ael59A8lTSetl3oTcFlEvF9HjIUZfqZULVXetAr738zbhSudmDuXjAcY9tEVvLC6mVkf1Gq3r4mIE0jPixciJeb9SstI2oR0W2kQcAJpgP+iwIDci/L/CkVL6r4K+AhpIoILgLdJS7RNAKbkqfxqhtj5dzWP58s1M7OKWi4pQ9nEvG9Jkb1JCffnEXFmRDweEe8U3Rpeo0rdr0fEhRFxeESsncveQWoB/6iO8Aqt3bU68ZbMzMxqarXb1/NExAl5ZrtvAJfmW9mFIU5L5+0CPZXzWOMdO3GdJyX9kNTJ6hN1nHI7aRafLSWtExGP1XstM+u/QjB3UNnh0A2z7rBnm3r92wZv3NTrA+zyyKHNDoE0Cra8lmwpF5S0mC+VtE8+VHiee2jxguaSFiPdll6ytC5JG0raX9IiZS712byt2fEr99A+j/Rvd7WkNUuuMzAPyzIzM+uUlm0pF5S0mH+bf59AGtO8EfBUXltVwCdJYwkvYP7UfAUrk6bfe1fSQ6RW9mDSmMPVSD2zv1dnWCeQJhzZlTQn7l+AZ0nDrtbL2+Z+JTYzs7bT0i3lgpIW829J89iOJPVGfpvUY3skcA0pUZebgOM+4ETSkmwjSPPY7gC8C/wYWC8iJtUZzwxS6/qQXN+6pA5jawN/xyvImJlZF7RESzkiVqmjzAmkFmqxL1UoPi7/FJ//AulGfuWb+Z2IK08AcnH+qVXPGOYvPVfu+DhK4jUzs/6nLVrKZmZm/YGTspmZWYtwUjYzM2sRLfFM2cysMyStTRqBsS2wEmm0w8ukURB/AW6OiFubFqBZFzkpm1lbkbQ/8BvSkMbpwETgv8CypNEXW5AWtnFStrbjpGxmbUPS8qR5CAYDXwfOjog5RccLyyxu3ZwIzbrHSdnM2slupJXh/hIRPys9mIcq3pV/zNqOO3qZWTtZLm9f6uyJkoZJ+qakByW9Kek9SY9KGlc8XW9R+XGSIm9XlXSxpBclvZ/PO16SGzbWo5yUzaydFFZp217SuvWeJGkE8ABpJbiVSZ3B/ggsBXwfuEfSUhVOXxWYROpUNpG0qtxqpFkGr8y3zM16hP9nMrN2ch3wHGn99L9Kuim3fneQtES5E5QmzL8C+BjwC2DliNg5IvYgzWF/MbA+8NMK1zyUlIhXj4j9I2KXXH46abreI3vs3Vm/56RsZm0jIt4izVk/idQnZldS6/dW4DVJ9+Te2cU+TeqRfR/wtYh4r6i+90hJ9SVgdIXW8rvA0RHxftF5/wa+m19+vSfemxm4o1dbGvjkTJbZu+Yqk00196OrNjuEmuasvmGzQ6jplksvaHYINQ1cobHXy2uYbyJpS9JiNJuRhkItBWxJWut8lzznPKTEDXB17ghWWt87kiblcpuQbmsXuzUv2VrqUuB8YA1JH46I6aUFJI0FxgIMWWTJTr1P65/cUjazthQR90bESRGxAzAc2Ir5CfUwSfvm31fL2zNyx60FfpifuJctc6n/VLj+DOD5/HJEhTLjI2JkRIxcaMiwTr5D64/cUjaztpdbwPdK2pXUoWsj0vPeK4GBudidwNQaVbX2LSjr85yUzazPiIg5km4nJeVCq7ewvvqVEXFOF6pdpdxOSYOBws37BW5dm3WFb1+bWdvIPalr+UjePpu3f8jbfcuUrcdOkoaX2X8g6W/okxHxbJnjZp3mpGxm7eRoSRMkbVp6QNJCkr4I7JN3XZ631wKTgU9JOk/S0mXOXT6fW85Q4BxJQ4rKrw6ckl+e1bW3YrYg3742s3YyCBgDjJH0AvA34DVgadLY4RVzudMj4hZIz5sl7Qn8HvgScJCkh0m3tRcG1iSNYX4J+H9lrnkRqZf3k5LuARYjTSSyMHAD0JVb4mZlOSmbWTv5Namz1g7ApsB6pKk3Z5FuV18InB8RdxefFBHP5tb14cB++bzNgFdJz4N/DPyuwjWfIg2VOhXYDlgi77sA+Fm5YVZmXeWkbGZtI08ecm3+6ey575NatZ1u2UbEU8ABnT3PrLP8TNnMzKxFOCmbmZm1CCdlMzOzFuFnymZmZUTEOGBck8OwfqZtWsqSppaZs/Z9SdMkXSHpUz14rYm5/lGdPK8jnzemp2IxM7P+ox1byrcAL+TflwI2JM3Us6+k4yKi0pqoZmZmLa0dk/JpETGx8ELSINKMOkcBp0m6solT3p0InMb8lWPMzMzq1ja3ryuJiFnA8cBbwGBgpybG8nxETImIN5oVg5mZta92bCkvICLek/QvYGPgQ4X9kqYCKwOrRsTU0vMkTQQ+BWxb3PouKbMtcFKuexBpWr/TI+L6MmU7gMOAz0dER8kxAQeTZhT6BLAIaSahm0mt/2cwsz5rwOxgyKuzmhrDP99dsXahXjTsuWjq9QHu2+CqZocwby3Rctq+pVxkibx9sQfr3Av4E2kavz8AD5MWUr9O0nH1VpIT8sXAb4AtgQdJMxIJOBr4m6RNejBuMzNrQ30iKUv6OLAqMBP4Yw9WfQzwrYhYPyIOioitgN2B2cDpktavs56jgINIXxg2jIidIuIA4KPA2aTJ9K8sXoXGzMz6n7ZOypKWkrQLcA3pvRzbw528JkXEmcU7IuIG4FLSHYiv1lnP8Xn73Yh4tKiuOcA3gGmk2+z7lDkXAEljJU2SNGlWvN+Jt2BmZu2iHZPyHYVxyqQl235PSmi7RMS5PXytSyrsvyhvR9WqQNIIYDVgbtF580TEzKLrVKwvIsZHxMiIGDlIC9e6rJmZtaF27OhVGKcsYHngk6R1TX8jaauIeKIHr/WfCvun5u2IOur4cN4+n1epKeepkrJmZtYPtWNSLh2nvAIpUa8HXCJp84iot4tfI+8UNL/boZmZtbR2vH39ARHxPGnR8lmkRc9HFx2embeLVjh95RrVr1Jj//TaEc4rs2KVjlyrdaI+MzPro9o+KQNExBTgl/nlOEmFOwCFJLd26TmS1gVWqlH16Br7J9YR27Ok29MDSOOUS+MY1Jn6zGxBkgZIOkjSdZKmS5oh6bXcOfIUScs1O0azevSJpJz9kDSr1+rAIXnfbXn7TUmLFwpKWgnoID2XrmYTSV8v3iFpV1JynQP8os7YfpK3p0ia9wVB0kDgdOAjwNNA80e1m7WZ3JnyAVKHyd1IoxmuAe4lDZX8H+BJSfv1wLXG5I6mHd2ty6ycPpOUI+JloDB86X9ya/kc4BlgE+BxSddIuh14DHiT9KGt5ufAmZIelnSppLuBm0jP4k+MiL/VGd4vgd8CKwAPS7pF0m+BfwHHAv8F9o2IGXXWZ2aApKWBP5Nm3JsIrBERW0TEgRGxG6kz6InAUOAySXs3LVizOvSZpJz9hDRBx2rAYRHxX9IMXJeQEulnSD2mzwA+TXoOXc3vgJ2BV/O5G5IS+V4RcUa9QeWOZ6OBQ4H7gc2Az5H+/c8FPhERD9Zbn5nNcw6pj8eDpGGRHxgxERGzIuI04DjSnbELJA1veJRmdWqb3tcRsUodZd4mfTMu3vcMZZ7lZqMq1FO6/081A0znjQHGVDgWpHHKC4xVNrPOk7Q6qZMnwNFVhhxCuut1OGmUxleAcSV17Qx8CdgcGE6aA+FJ4Hrg53l+/anM7xx6mKTDiqq4MH/+zbqlr7WUzaz/2I30N+zRiJhUrWD+Uvyb/HL3wn4l55IWhtmL1Dn0atI89yuRlmItLHJzFXBP/v1J4MKin7t74P2YtU9L2cysxMZ5+0Cd5QuPiD4haaGImA18DTiS9Nhrz4i4r1A4LySzLanPBxHxDUljSI/E7nbL2HqDk7KZtatl87beleEK5QYAS0t6jbQsK8CY4oQM81rXt3c7SrNOcFI2s/6idAjkSNLz42cj4uZeuaA0FhgLMGTIEjVKm/mZspm1r1fy9kNVS81XmEBkLqkjV6HT1uM9GVSx4oVkBg8a1luXsT7ESdnM2tXkvN28zvKb5u3D+Xmy56O3luOkbGbt6kZSq3cdSZtUK5g7bR2aX96Qt9Pydq3eCc+s85yUzawt5WVaC1PTniNVXWj8GGBd0lS85+R9k0m3wEfkccr1KCxy4/441iuclM2snX2Z+VPp/l7SKsUHJQ2S9C3SbH8BHBERL0Ga7Qv4v1x0gqRNS86VpG0lFffQKixys06PvxMz/G3PzNpYRLwiaWvgOtKY4ick3U9a4GUxYEtgaeAd4IsRcUVJFT8lJdgjgPskTQKeyOd8jDSByKrAG7n8fcALwEa57KOk6XrviYgJvfZGrd9wUjazthYR0ySNBA4E9idNKrIJKRE/RZpf/hcR8UKZcwP4oqTrSJOIbApsQJrv/gngbFISLpSfIenTpFXptiDNhz+A9LfUSdm6zUm5HUUQM1p7QalnTh7Y7BBqWnzoG7ULNdkuH92q2SHU4YlmB0BEzAEuzj9dOf9GUsexeso+TJri06zH+ZmymZlZi3BSNjMzaxFOymZmZi3CSdnMzKxFOCmbmZm1CCdlMzOzFuGkbGZm1iI8TtnMrAE0czaDp71Su2AvuumujZt6/VWfaO35FVqBW8pmZmYtwknZzMysRTgpm5mZtQgnZTMzsxbREklZUnThp6PZcZuZmfWkVul9fWGZfcsDO5OWX7uqzPG7ezUiMzOzBmuJpBwRY0r3SRpFSsqvlDtuZmbW17TE7WszMzNr06QsqSM/Vx4jaX1JV0p6QdIcSceWlqlQx7h8fFyl/ZJG5Hqel/SupIck7VNUditJv5f0aj5+h6RNylxrlVznVEkLSfq2pMckvS/pRUkXSvpIj/4jmfUx+fNT6FOyW5Vy/8hlRjUuOrOe0ZZJuchWwAPARsBE4Gbg3R6qexVgMrANcCfwELAhcIWkAyTtBdwBDAduBZ4GRgF3SFqzSr2XAycD04BrgRnAocCDktbqodjN+rpTJbX73y+zBbTEM+VuOAL4IfC9iJjbw3UfBpwFHB8RcwAkHQX8EjgDGAaMjogr87EBwKXA/sC3gMPL1LkysAiwYUT8M583GPg1cDBwEbBpD78Ps77mXWA9YDTpM2PWZ7T7N80pwPd7ISEDTAW+WUjI2XjgVWAEcHMhIQPkGH6UX25bpd5TCgk5nzcT+CrwJrCJpK3KnSRprKRJkibNwvPHWr/287w9OX+pNesz2j0pX1eSNHvSHTlhzpOvNTW/vLnMOf/O2xWr1Htx6Y6IeB24Ib8cVe6kiBgfESMjYuQghlSp3qzPu5r02GpV4Mh6T5I0SNJXJN0v6U1J7+W+HadJWqak7Kfzc+m/VqlvaUkz8s/SXX43ZkXaPSk/3Yt1P1th/9uVjkdE4VilrPl6TsDlTM3bEfUEZ9bPnZi3J0latFZhSQsDfwTOBtYF7iJ9EV6S9LhpsqTVik65FXgO2EDS+hWqPRAYDNwQEa915U2YlWr3pPxeN86t9d5r3RLvjVvmZlaHiLidlGSXA46v45QfkO5CTQHWjIjdImI/YDVSy3tl4JKi+ucAv8kvx1So87C87ehc9GaVtXtSrqZw67nSt+iVGxVIkSUlLVHh2Cp5O71BsZi1uxOBAI6XtGylQpIWAY7KL4+JiHmfsYh4j3QL/G1g85I+HR15O1rSBzrFSvoYsAnwAuUfZRXKzesLMnNOd9oQ1l/05aRc+OCtXXogf0hHNTSa+UaX7siJujDucmJDozFrUxHxEHAFsBhwUpWiG5O+nD8XEbeWqecVyvTpiIjHgftIrfFdS04rtJIviYjZVWKc1xdk8MBFqr8hM/p2Ur4tbw8pHv+bE/K5QLMm6/iepHWK4hlEGnq1BDA5Ijynt1n9/geYDRwpqdLdrw/n7X+q1PNUSdmCCXlbSMJIGkgawgi+dW09rM8m5ZzcbgQWBx6SdLOkG0kfzJ2Y/2FrpGnAvcDfJP1B0mXAk6QP/CukSUTMrE4R8QRwPqlz5Q9qFe/CJS4n9V3ZraiH9g6kERaTI+IfXajTrKI+m5SzfYHTgJeA7Ugzf92Ut9OaEE8A+wGnkDqY7AksTBomtUnx+GUzq9sPSBOKHCxp3TLHC4+yVq1SR6Hn9Qf6dETEG6SZ9wYDB+XdY/K2o/OhmlXXskk5IiZGhCJilTLHxuRjHTXqeD8iToyIVSNicESsGBGHR8QLETEu1zGu5Jyy+4uOj8rHJ1Y4rohQlZhmR8T/RsRaEbFwRCwXEYdExNRq78XMyouI50mPgAYAp5YpMpnUkevDkrYvPZhbwJ/NLyeWOX/eLezc/2NPUkfSS7sVuFkZLZuUzcw64XTgv6Tk+oEWce5hfV5+eZakFQrH8vjlc0kdwe6LiHvK1H0b8Aypw9gPSHe3PDbZeoWTspm1vTwpz//ll0PLFPkuqRX8ceDfkq6XdDmpg9e+pMdZC4yMyHXPZf4c28fkbUdPxG1WyknZzPqKs6kwE19EvE/q4HkM8E/S/PR7kOacPx3YKCKeKndu1lH0e9WxyWbd0e6rRLWF/Ly44nNmM6utXP+SkuPvAytVOT6LlLjP7sK1/40/w9YAbimbmZm1CCdlMzOzFuGkbGZm1iKclM3MzFqEk7KZmVmLcFI2MzNrEU7KZmZmLUIRXVk4xZpJ0svA0z1c7XDSSlWtqtXjg/4b48oRsWwP19nn9MDnthX+/3IMPXP9ip8ZJ2UDQNKkiBjZ7DgqafX4wDFa72qF/3aOofev79vXZmZmLcJJ2czMrEU4KVvB+GYHUEOrxweO0XpXK/y3cwy9fH0/UzYzM2sRbimbmZm1CCflfkjSWpK+JuliSVMkzZUUkvZpdmwAkgZJ2l7SjyVNkvSmpJmSpku6StKoZscIIOmrkq6Q9JikVyXNkvSypD9JOlhSyy31J+nU/N86JH2j2fFYZc3+nLbC57BVP2O9+Tnyesr901HA15odRBWfAm7Nv78A3AW8A3wM2BvYW9IpEfG9JsVX8C1gOeAfwL2kGFcGtgO2B/aR9LmImNu8EOeTtAnwTSDw2sDtoNmf01b4HLbcZ6y3P0duKfdP/wDOAPYH1gDubG44C5gLXA18MiJWiIjdImL/iFgPOACYA3xX0rZNjTLFslREbBQRn42IAyJiC2A94EVgD+CwpkaYSRoCXEiK67omh2P1afbntBU+hy31GWvE58hJuR+KiPMj4psRcUVEPNnseEpFxO0RsU9E/LnMscuBjvzy4IYGViIi7o6Id8rsfxQ4J7/csbFRVfQDYB3gSOCNJsdidWj257QVPoct+Bnr9c+Rk7K1o7/m7YimRlHd7Lyd0dQoAEmbAccDl0bEDc2Ox/qMZn8OG/oZa9TnyEnZ2tFH8/b5pkZRgaRVSd+kAa5vciwLk263vUZr9yOw9tO0z2GjP2ON/By5o5e1FUnLA2Pyy6ubGMo8kj5P6hQziNRq2JL0hffUiPhdM2MDfgisBRwQEc1eSMD6iEZ/DlvgM9awz5GTsrUNSQsBFwNLALe10K3YrfhgZ5PZwHeBnzQnnETSlsCxwLX5GaBZtzXpc9i0z1ijP0e+fW3t5DzSMIhnaHInr2IRcURECBgKfBz4GTAOuE/Sis2ISdIipI44bwJHNyMG67Ma/jls1mesGZ8jJ2VrC5LOAg4njZfcPiJeaHJIC4iI9yLinxFxAnAi8AngF00K51TSM7/jIqIln71b+2n257AJn7GGf44897UhaSLpec2+EXFVk8NZgKQfA8cBLwOjIuKfTQ6pJknLkBZCnw0MjYhZDb7+VGAlYIHhLMDawIeAp0itnSci4ojGRWdd0ezPaat9DhvxGWvG58jPlK2lSTqd9IfgVWCHZv8h6IT/kv5YLAQsTZpsoNEGkP6IV7Ja/lmyIdFY22rRz2GjPmMN/Rz59rW1LEmnASeQPnw7RsTfmxxSZ3yS9MfiddK3+YaKiFUiQuV+SEM7AE7I+zZodHzWPlr4c9jrn7FmfI6clK0lSfpf0ry3r5P+EPy1+hmNJWlrSbvlnqilx7YCfp1f/joi5jQ2OrOe0czPYX/9jPn2dT8kaSPgl0W7Ppa3pxaveBIRmzc0sEzS7sBJ+eUTwFcrLAYzJSJOa1hgH7QGMAF4XdJDpI4viwGrM//f8ybSsA2zTmv257QFPof98jPmpNw/LQ5sVmb/R8vsa4ali34fmX/KuRNoVlK+EzgF2Ib077YlacWYF0iTKVwcEdc2KTbrG5r9OW3257Bffsbc+9rMzKxF+JmymZlZi3BSNjMzaxFOymZmZi3CSdnMzKxFOCmbmZm1CCdlMzOzFuGkbGZm1iKclM3MzFqEk7KZmVmLcFI2MzNrEf8fKkYQceDBgQsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 504x360 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "w = calculate_W(As[0][0], As[1][0], As[2][0], Ss[0][0], Ss[1][0], Ss[2][0])\n",
    "\n",
    "fig, axs = plt.subplots(1, 2, constrained_layout=True, figsize=(7,5))\n",
    "color = \"viridis\"\n",
    "\n",
    "A_A1_new = As[0][0] @ w[:,:,0].T \n",
    "B_A1_new = As[1][0] @ 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": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAU8AAAFwCAYAAAA44fmKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAA/mElEQVR4nO3dd7xcRf3/8dc7IaH3jiChSFFASuigoQsiRaWGEgURUBFBvop81SDqFwELIoL5IVykE0CqgAoEpBOaigQFDIFA6D2kv39/zGyybHb37t27d8vN5/l47ONkz5kzZ/Zm7+fOmZkzI9uEEELomQGtLkAIIXSiCJ4hhFCHCJ4hhFCHCJ4hhFCHCJ4hhFCHCJ4hhFCH+VpdgNBzgwcs6AUHLtrqYlQ338BWl6Bb70yZ9JrtZVtdjna1zFIDPWSVQS0twz9fa+1/z/S33mDm+++r3LEInh1owYGLsuUy+7S6GNUtsVirS9CtW8ed+lyry9DOhqwyiAdvXaWlZVjnvKNaev0J5/yy4rG4bQ8hhDpE8AwhhDpE8AwhhDpE8AwhhDpE8AwhhDpE8AwhhDpE8AwhhDpE8AyhSSSNl+RuXnvVmfeQfP74xpY6VBKD5ENovluBSRWOTWhmQUL9IniG0Hyn2h7T6kKE3onb9hBCqENLgqek3SRdJOlpSe9JmirpBUk3STpS0qJFaUfmtpyRJXkMy/vHNLv8ITSDpF0kXS/pZUnTJL0k6TJJ63dz3nySvivpSUlT8vkXSvpos8o+L2hq8JS0XA52NwEHAdOAPwN/BMYDOwLnAM9KWrWZZcvls6RYES+0nKQzgVuAXYFngGuBl4D9gQcl7Vbl9CuAk0ntp9cCU4FDgIckrd13pZ63NK3NU9ISwD3AmsB9wJG2/16SZlHgKOAkYEmg2qw3DwLrApP7orwhtIqkI4FjgCeAL9oeV3RsL2A0cImk1W2/WXL6qsCCwEa2/5XPGQz8nlRhuQjYrM8/xDygmTXP35AC54PA9qWBE8D2u7ZPAzYBXq6Wme3JtsfZjt7J0GnuqDBMqUvSQOAHOd2+xYETwPa1wO+AJUjBsJxTCoEznzMN+AbwDrCppK0rFUzSEZLGShr76usz6/6A84KmBE9JawAH5LdH2p5SLb3tp22/1E2eZds8i8e7KTla0mOSJkt6U9J1ktYrOWdk8e166Ze6JK0kHSxpTM5viqRnJJ0tqezkh8X5SNpP0n25rfddSbdJ2qbaZw39zq3AhWVedwMbAisCTxQHwBJ35u2WFY5fXLrD9lvADfntsEoFsz3K9lDbQ5dduv0ntG6lZt22704K1P+w/WiTrgnQBewH3AX8B9gU2AMYJmkj28/mdI+RvryH5vcXlstMkkhfzAOB6cAY4A3SbdDRwP6SPmP7oQrn/4jUJHE3qd13A2B7YBtJw2zf14vPGjpHxaFKkgqzXH+ihvb3ctOsv5UDZTnj83bl7goYutes4LlJ3pYNKn1kVWBb4BO2nwGQND9wDbAbcCLwFZh9K3StpEPz+xEV8jyKFDhfBnaw/UTOdyDwS9Kt0WhJa9ueWub8rwGb2X44nzcAODeX40fATr37yKEfKFT3JgJ/7SbtuG6Ohz7UrOBZ+Av5SpOuV3BMIXAC2J4q6WRS8NyhjvyOz9vvFwJnznempG8De5KC9heBS8qc/8NC4MznzZL0fVLw3FbSINvTy11Y0hHAEQALDFikjqKHDvF83r5U5Y94NUtIWtz222WODcnbifUULHxYfx4kP4M01KNU4a/1Sj3JTNLKwOrALFKP5YfkRvlCwBxWIZsby5z3MvAmMD+wdKXrF7dFDR6wYE+KHjrLg8DrwEaS1qwzj+GlOyQtTmo+g9TcFHqpWcHz1bxdrknXg/SXe0bpTtvv5H/O38P8PlKUb6UOr2dL0paqNDKgUKYFelim0M/kO49TSLfv10qaa1iRpMGS9pC0ToVsfiBp3aL0g4AzgcWBh23f3QdFn+c067b9YeBgUodNs8zqo3zrHkRvu6/KFPoR22fmh0S+BTwg6e+kgfLTSH+YNwIWJg2gL233nED6fXtM0u3A28BWwCrAa6TB8qEBmlXzvIkUzNaXtFGTrtlohXailXLHUzmrl6QNoS62jwM+DVxOemDks8BngGVIzT/Dgb+VOxXYl1R7XR3Yi3RHczGwaZXhT6GHmlLztP20pCtIYz3PkfTpCr3RwOxxoZO7G+vZB6YDgyTNV3rLb/sFSc+SvpAHkZ7YmC3fGhXamsY0oayhw9ge0sP0d5GG2dWSdjygol0/zq/QR5rZYfR1Upvg5sDt5SY3kLSwpONItx3LN7FsBYUa47oVjv8ib08pbm/KQ5VOAz5KeqT0qj4rYQihLTTt2Xbbb+Qnaa4EtgH+LulfpDabQlvOZqSOnJdJg8+b7Y+kdqbbcnvRe7nsh+fjvwW2JtWgH89PNxUGya9O6jXfp1qtOoTQPzR1MuR8G76tpN1JAWhLUjvOfKQe+b8C1wGX2n6/mWXLTiK1Ge0NfB4YlPcfDmDbkoYDN5PGZm5OmoThRdJsUP9n+/nSTEMI/U9LZpK3fSNlxjxWSDsSGFlm/xg+3MZT2D++3P6SNGWP2/6ANBD++HLHcxqTxnnONdazp9crOj6k1rxCCO2hPw+SDyGEPhPBM4QQ6hDBM4QQ6hCrZ4YQ2pZbPaVold6KqHmGEEIdIniGEEIdIniGEEIdIniGEEIdIniGEEIdIniGEEIdIniGEEIdYpxnJxogNHhwq0tR1WublVsVt800ae1JSeNJCwMC/Mz2d6ukvZg588LeaXtY35Yu1CtqniE018F5/te5SFqMNJtX6AARPENonrGkVVt3qnB8f9IUhw81rUShbhE8Q2ierrwdUeH4CGAmPZjuMLROBM8QmucB4ElgT0lLFB+QtDZpcvBbgbJrd0naUdLZkh6X9LqkqZKek3Rh8VLDJed0SbKkEZI2kDRa0iRJMyUd29BPN4+J4BlCc11AWs3ygJL9I4qOV3IucBgwg7Qw3J9IS9gcAozNy9xUsjXwILAxaYHCW4DJPSt6KBbBM4Tmuoh0az6isCN3IB1CWg/r+irnfhtYwfYmtve2vTewFnAksBAwSlKleYAOB84APmZ7f9uftT2q159mHhbBM4Qmsj2JVOvbrOhWe2dSR9KltqdVOfda22+V7LPt3wH3klZ9/XiF08cBP7Q9q5cfIWQxzjOE5usCPkuqfX6HObXQru5OlLRyPncdYDGgMOxphbxdC3iizKnX2Z5ZQ/5HAEcAfPQjER6qiZ9OCM13PfA6aczn6cCewD9sP1ztJEknA9+j+u/tYhX2P1dLwfKt/CiAoZ9cwLWcM6+K2/YQmizfml8KrEjqIJqf6h1FSPoC8APgA9Ky12sAC9lWXp31skLSCll80ICihyIRPENoja683Z3Ue35JN+n3ydvv2T7P9rN5qeyCNRtcvtCNCJ4htIDtR4B7SLfvo22/0s0pS+Xt86UHcsfTRo0tYehORwVPSePzgN9hFY6vL+mlnOYqSXXNnpEHFFtSVy+KG0JVtrexvYztA2tIXpjG5CvF32tJywEXEv0XTddRwbMaSZsDd5J6HS8A9qs27COEDvMr4G1ST/vT+UmhG4FngEWAa1tXtHlTvwiekrYH/gosSfqSHVbLsIwq/kgaM3di70sXQu/ZfpZ0a345qVPoc6Tv6CjSY51vt65086aOr+pL2gO4ktRjebLtkb3N0/bbxJcxNIjtIT1MfxVles1t/5e5H+ssGEGZCUdsl90feq+ja56ShgNXA4OBY8sFTkmbSzpd0lhJL0uaJunF3Ca6RYV8y7Z5ShqW94+RNEjSSZLGSZoi6RVJF0v6aJXyriLpTElPSfpA0juS7snXqzTEJITQhjo2eEo6ivScsIAv2z6zQtKfAN8CBpEmRigMUP4CcLekfSqcV80g4Gbgu8DT+d+zSDOA3106Y04u73bAP4BjSD/3W0iz7GxAaqO9sI5yhBBapFNv248lPZUxDTjA9jVV0p4BDLf9cvFOSZ8j1VrPlXST7Z7MMLMVaWLbNQpDTCQtDtxOmrXma6SgXbjWivlai5Buof5g2/nYKqSAfrCk22139aAcIYQW6dSa5555+6tuAie2bykNnHn/DcBo0vi57Xp4fZNqu7PH5uV20p/ltzuUpD+W1Jn1c9sXFgJnPu950hMjAN+odEFJR+Smh7HTZsbDIiG0WqfWPO8EPg0cL+lh21dWSyxpGdKTHOsBSzDnc6+Xt2sBN/Xg+hNs/6PM/sJYvJVK9u+Wt6Mr5Pcw8B6woaQFbE8pTVD8zPHi8y8fzxyH0GKdGjxHksa7fRu4RBKVAqikrwK/IM13WEmlyRQqmVBh/zt5u0DJ/tXz9qEa+oWWBib2sDwhhCbr1OCJ7RNyIKoYQCVtCpxDenb4BOAG4AVgsm1L+ilpLGdPe7p7OidiYdqwK4C5apUlpvYw7xBCC3Rs8ISyAdS2i2+Nv0AKjL+2fUaZLJo1mcLz+Vqn2C4312IIbefFGQvww1c/0dIyLPBaa0fwaUblY53aYTSb7RNIPerzAZeWDD2qNpnCslReArbRbs7beoZFhRDaUMcHTygbQL+YDxU6cA6RtEghvaRFgfNJnUfNcDqpPfR7kr4maa4av6RPSPp8k8oTQuilfhE8Ya4AelkOoBeQap0bA89KukbSH4HxwFBSAG1G2Z4H9gLeBX4DTJD0F0mXSLpJ0gTgn8C+zShPCKH3+k3whLkDKGm85VDSEJ/3SD30Q4FrSAF1rtv5PizbHcAngJ8CrwBbkNpkPwE8S+q4OqlZ5Qkh9E5HdRjVMsFCDqAnlOz+aoXkI/OrNI8uyizGZXsMVXrmbY/v5vgkUoCMIBlCh+tXNc8QQmiWCJ4hhFCHCJ4hhFCHCJ4htIikAZIm5DliX5U0qNVlCrWL4BlC6+wErJL/vQywRwvLEnoogmcIrfPlvJ1Y8j50gAieIbSApKVI89Ia2B+YCewiqXQ6w9CmIniG0BrDSYsWjrF9N/Bn0uxbh5ZLnNtFC6sPHCHpUUmTJb2en5xbr4bzDpP0QF47y+WWiwm1i+AZQmsUbtG78vaCvP1StZMk/ZI0zeLbwHXAa8DewAOStqly3lmkJ+2mAjeSJuCOSbV7oaOeMAqhP5C0EbAhaa6Dq/Lu64E3gI9J2tb23yqcfgSwne27cl4iPfL7XdKkOGuVW4kAOBjY0vaDjfsk87aoeYbQfIVa55WFhQdtTwUuKTlezjmFwJnPM/C/pPkRViHNl1DOaRE4Gytqnp1o4EBmLd3TlUOaa8nLHmp1EdqSpPmBA/PbC0oOX0BaBHAfSd+w/V6ZLC4u3WF7pqTLSHMmDGNOEC5WdaHEovIdQardstiKC9Zyyjwrap4hNNdepEm6/2P7nuIDth8FHgcWBvarcP5/K+wfn7crVzj+XC2Fsz3K9lDbQxdecnAtp8yzouYZQnMVbskXl3R3mePLFaX7faMuajvWq26wCJ4hNImkVYAd89vlmBMoy9lK0tq2nyrZP4RUOy01JG9j5dUmidv2EJpnBOl37nbbqvQCCqvAlus4Gl66Q9JA0kB7gDF9UO5QRgTPEJogDykakd9e1E3ywvGDc2AsdnTxeM6c78nAGqRa59W9L22oRdy2h9Acw4DVgQ/oPsDdArwKrAjsBtxQdOz/AXdKugt4ibSczNo53+HRttk8UfMMoTkKt+DX2n63WkLbM4DLS84rOI40nGkpUs/9csC1wOa272xUYUP3ouYZQhPYPpj0lE+t6Y8Bjimz38Bv86uWfCquqRV6J2qeIYRQhwieIYRQhwieIYRQhx4Fz7zmyoGSrpM0UdJUSW9IGivpFEllB/1KGpnnDxxRsn9EYb5BSR9Um19Q0seK0lrSsJ6UvUq+s+c7LNk/ppHXCaE3isaAhjZRc/CUtDLwIGnSgd2BCaTJBu4FViPN7PKMpH3rLMsCzBnoW86IOvMNIYSGqyl45iUD/gZsQnqCYU3bW9o+wPbuwArAicBCwOWSKk2LVcmjpGUIRlS4/gDgEOBN4Oke5h1CCA1Xa83zbNKzsw8Bu9r+0MwutqfbPpU0Bk3A+ZKW6UE5XgT+AmwuaZ0yx3ckzRZzOWkm7BBCaKlug6ekNYDCrfjRFWapLvg18A9gMeDrPSxLYW7DEWWOfakkzVwkLSvpm5JukfRfSVMkvS3pfklfK/OYW69I2kXS9ZJeljRN0kuSLpO0fkm6bXPb6ZNV8loml/cDSUs3spwh1GuWBzB55uCWvphFa19VFiqppea5e073hO2x1RLmAbx/yG97ugb1daTb8g89z5s7kfbK1682w+4uwK+AT5DmPPwj8AhpuYPfAFfn54B7TdKZpEfodgWeIT3h8RKpzfZBSbsV0ublFB4H1pG0fYUsDyctBna57dcbUcYQQt+qJXhukre1TuFfCHCflDQfgO2Rubewq9JJeRmCy4CVgJ2LDu1P6kyqeG72MLCF7VVsb5/bY7cjdWY9Rlrmtd7OrNkkHUl68uMJYH3bW9ne1/bGpIW45gMukbRk0Wln5e3RZfIbAHw1vz27t+ULITRHLcFz2bx9ucY8C+kGkJ6/7Ylyt+5fAmbQzUw0tp+0/UCZ/S8B/5PffrGH5fmQXCP+QX67r+1xJde6FvgdsARwUNGhS0mLe+1ZZl3uz5Lbk7ur2YcQ2kdfDJKv+9Y4B49/koLMkpLWBTYDbrHdbfCWNJ+knSX9UNI5ki6Q1AUcmZOsVW/Zsg1JM908YftfFdIUJmfYsrAjz3RzHqlWekRJ+kJttGqtM6/VPVbS2GkzJve03CGEBqtlYpDX8nb5GvMsDJSfRapt9VQXcAZwAHNmx67YUVQgaS1S2+O6VZL1dtW01fP2E+UG1pdYtuT92cDxwFck/dj2jNwZtwvwOnBFtcxsjyKtu83iC60U622H0GK1BM+HSbegW9SY52Z5+3ieWqunLgZOJU3FtRIpeN9Q9YzkKlLgvB44DXgSeDuvLLgW8BS9qBVnhY6sicBfu0lbeks/QdL1pHbRvXJ5j8plOr+bUQwhhDZTS/C8Efg5sK6kTav1eOfe7EPy21oC3lxsvyzpZuBzedevbU+vdk4eG7o+8ArwedszS5KsWU9Zyng+b1+yPaKO888iBc+jJd1Ias+dBZzTmOKFEJql2zZP20+TakkAZ0taoEryY4D1gHfpXc/xeaRb2dep4ZadOR1TL5YJnFBm3Zc6PZjLtJGkHgdk23eQ2nS3Iy2dsBRwc+lDByGE9ldrh9HXSLWuTYE/SRpSfFDSIEnfAX5BGlZ6uO1X6i2U7ettL5Nfj9Vwyn9INbj1JH2qpGxfIrWf9lquAZ9Cun2/VtJmpWkkDZa0R4UnpSCNOYU5IwBqmtQ2dA5J44smsDm1m7QXF6Ud06QihgaoKXjafg3YhjRecjvgaUn3SLpU0g3AJFI7ZWEdlSsrZtYHbL9KCkLzAXdIuj2X7R/A+blsjbrWmcAvSYPxH5D0uKRrJF0u6W+kTrLrmNPZVeoi0sMAAM+SBtuH/qvcIm4ASFoM+HyTyxMapOahSrYnAENJSwncTBp8/kVSUB0P/IQ0YchljS9mTb5JGgb0OKnTalfSmNNdyb3UjWL7OODTpGftlySN1fwMsAypjXg4aSKVcudOJs1EBXCO7VmNLFtoK2NJnZ47VTi+P7Agcx4sCR2kR2sY5fbEi/Or1/ITR109PGe9CvtnkVYW/H8VTi3b015pjkTbw7opx13AXdXSlC1EmvN0J1It/fyenh86ShepwjGC8ncYI0iziV1EahILHSRmkm++k4DBwIW26xkHGzrHA6Qhc3uWTvQtaW3SgxS3kuZFmIukHSWdnZuGXleafPw5SRfmB0hK09+W204rzosr6ec5zWm9+WAhgmdTSNpK0nmS7iCNSHib1PEU+r8LSHMzlHZajig6Xsm5wGGkx5PvAv4ETCMNBxwraZuS9BXnUACQtCAxPK5hIng2x1qkX4LNgLuBz9h+sbVFCk1yESUTfecOpENInYvXVzn328AKtjexvbftvUnfpSNJE4+PKpkp7AbgOWBbSeWatw4gtdHH8LgGiODZBLa78qxSC9ve1vb9rS5TaA7bk0jtnZsV3WrvTOpIutT2tCrnXmv7rZJ9tv07UqfjusDHi47NZM7Qt3K1z8K+isPjiudQmPxWzDteTQTPEPpeV96OKNl20Q1JK0v6qqRfSvq9pK482c0KOUnpZDfnkTojD5K0aFE+m5Oml6w6PM72KNtDbQ9daIn5uyvePK1Hve0hhLpcT3oy7WBJp5Pmlv2H7YernSTpZOB7VP89/dBkN7bfkHQJaYLtg5m7JhrD4xokap4h9LF8a34paTrDC0irBlR97FhpEcUfkGqRXwHWABYqWoK4MJ663FC7QsfRUTmvpUkTgU8hhsc1TATPEJqjK293J/WeX9JN+n3y9nu2z7P9bJ4XtqDi3Aq2/07qnS88rnwYqcf/8hge1zgRPENoAtuPAPeQbt9H1zD3Q2Gym+dLD+SOp426Ob9Q+/w6cyYDj2VeGijaPENoEtul4zKrGUd6Eu0rkm4u9MrnJ9QupPvf3WtJgbdQg41lXhosap4htKdfkR6m+CxpIp7ReQ7YZ4BFSMGxojwRefFA+Kh1NlgEzxDakO1nSbfml5M6hT5HGtc5ivRY59s1ZPOXvO12mZfQc3Hb3oH8wRRmPVZp/bn2oPk7YIxg1fUJ6md7SA/TX0WZXvP8FFCluWhH8OFVZsspTAL++1jmpfGi5hlCPyRpFdIQp2nMmYA7NFDUPEPoR/LM9SuTOpsWBk63PVePfei9CJ4h9C/7Ax8lTXP3M+D7rS1O/xXBM4R+pKftraF+0eYZQgh1iOAZQgh1iOAZQgh1iDbPEEJZkhmkmS0tw8wFW3r5qtXLqHmGEEIdIniGEEIdIniGEEIdeh08Ja0j6VxJT0maLOkDSRMk3ZvXiN6pEQXtC5JG5jWsR7a6LCGEztKrDiNJ+wF/AAYDE4ExwJvAssDGpNlfPs2c2V1CCKFfqDt4SlqBtB7KYOBbwFl56dPC8QHANvkVQgj9Sm9qnrsDCwH32f5V6cG8Qt9d+RVCCP1Kb9o8l8vb7tZi+RBJm0s6XdJYSS9LmibpRUlXSdqiwjmz2yYlLS/pd5JekDRV0n8lnSppgQrnDpL0bUn/kjRF0iRJF0latcayXp6vNU3Sq5Kul1S2Np3L6PzvwyQ9IOmdvH+JvH8lSb+R9HQuz+TcRnyLpCNq/kGGfkHSmPz9GNHqsoSe6U3wnJC3O0harwfn/YR0mz8IeJA5a1p/Abhb0j5Vzl0FeJhU672P1Ma6HPAd4MrSxLnp4BrgdGA14HbgTmCHnM9qlS4k6fh8jX2BScB1wNOkZRHulPSVKueeRZrxeypwY76WJa2Y//01Uq3/FuAG0s9yC+C4Kp89hNBGenPbfh3wIrAS8KikP5MC0yOkxaYqLRNwBjDc9svFOyV9DrgaOFfSTbYnlzn3y8B5wNeKFsRalxSEPydpa9v3FKX/GinQTgSG2X46n7MAcDFwSLkCSto1l/NF4PO2Hyg6tjXwJ+BsSXfa/neZLA4GtrT9YEm+3wRWAH4HHGXbRcfmBzYvV54QQvupu+Zp+11gR2AsKQjvRpo/8C/AG5Luyb3xpefdUho48/4bgNGkJVe3q3DZ54FjCoEzn/ckcFF+u0NJ+mPz9n8LgTOfMwU4GviA8kbm7eHFgTOfew9wCqnm/NUK559WGjiz5fP2luLAmfOdajvah0PoEL0a52n7SdubAlsDPwVuIw1VGgBsBVwuqav0PEnLSBoh6QxJ50nqyukKt/9rVbjk7bbLBbxxebtS0TVWBlYHZgGXlin7K8Cfy5UN2Ax4p9zx7M683bLC8Wsq7C8E1J9J2kvSwhXSzUXSEbmdeOx0ptZ6WuhA+fehYjtopfHJvekbCD3XkIlBbN8L3Auz2xm3AH4I7Awcmm/DR+fjXwV+Qeqpr2SxCvsnVNj/Tt4WfzFWztsXi2uqJcaX2VdoB10MmCHNtS5XsWUr7H+uwv6LSD+TA4E/AjMl/ZM0IuHy/HMsy/YoUjsqi2kpV0oXAnP6BkT6vVyMNGTwO8DHgT1aV7T+o+GzKuUhSvdK2o1U09oY2AsYLWlT0lrSM4ATSJ0lLwCTbVvST4ETKbOSYDar0eUtY2Devk03a2MDr5XbWaF2XPjZDJf0f6S22K3z6xvANySdb/uwegodQpGe9g2EOvTZlHS2Z0q6nRQ8CzW0L5AC469tn1HmtDUbWISJebuSpMEVap9DyuwrLJY13faIBpZnNtv/BP4Js2vqu5GaFr4s6QrblZoLQqhF2b4BSRcBR5H6BiJ49lLdbZ7q5n42+2jevpC3S+XtXKv5SVqWtOJfQ+QVA/9L+oz713o92xOBfwDLSBrWqPJUYnuW7RtJoxcAPtnX1wz9Xs19A6WK29Ynvxlt69X0psPoaEkXSNqs9ICk+fI4yC/mXVfkbeE/7xBJixSlX5T0qOcSvShPOb/O2x9LWr3oevMDZ1O53bWw4uDFknYuPShpoKTtVWFQfyWSDpG0cZn9SzOn86lSe2kItepJ38CH2B5le6jtoQstOX/jS9aP9Oa2fRAwAhghaRLwGPAGqXa5AXP+up1m+9b87wtIw4c2Bp6VdDfpNv5TwDRSAP1yL8pU6ixSB82uwBO5GeE9UuP5AqRJTeYa62n7ujxI/jTgVkn/Bp7K564AbEQK9EcB9/egPJ8HLpQ0kfTzegtYGtiWtMb230gdSSFU012lpxl9A/O83tQ8fw/sTarBPQ+sD+xDCoTvARcC29r+TuEE228CQ0m9xu+RntYZShraszFlbud7I09UsifwXVLP+o6kMaR35ev+t8q5vwA2IX3OgaRb/M+RevHvAr5CmaeauvFz4EzS4PuhpJ/XBqQHCw4DdrI9vYd5hv6n0Fa5SIXj3T5aHPpe3TXPPEj+WrrvkS497xUqDy4fyZwB6sXnlN1fdLwL6KpwbDpp8P7Par1e0bmPA4dXOl4mfdV2YNt/I9UuQ6im0Nm5TukBSQsCw5pamlBWzCQfQvu5LW8PlrR2YWcOnOcwpyM2tFAEzxDaw+x2Stt3kyaUWQx4JM+4dSOpmWlnUt9BaLEIniG0VmFx3fdL9u8DnEqa8nF7Up/ATXlbqTc9NFGs2x5Ci+Sx0oUhdOOLj+XJa07Mr1IjaXDfQOi5qHmG0DqHAssArwKPt7gsoYei5hlCE0laiDSf6xrMeTDi+7ZntK5UoR4RPENorsHAQcC7pOfLzyzMOBY6SwTPEJrI9ltUnjUsdJBo8wwhhDpE8AwhhDpE8AwhhDpEm2cnWmRBZm20YatLUdWgV95tdRG691SrC9De3p8xmIffaO2ToNVni2itqHmGEEIdIniGEEIdIniGEEIdIniGEEIdIniGEEIdIniGEEIdIniGEEIdIniGEEIdIniG0MYkjZFkSSMakNeInFdX70sWIniGEEIdIniGEEIdIniGEEId5ongKWk3SRdJelrSe5KmSnpB0k2SjpS0aKvLGEKtJHVVaweVNDIfH9ncks1b+vWsSpKWA64EPp13PQn8GZgGrAzsCOwGnCJpqO3nWlLQEELH6bfBU9ISpDVi1gTuA460/feSNIsCRwEnAUsCETxDCDXpt8ET+A0pcD4IbJ/Xwf4Q2+8Cp0m6Bni/yeULIXSwftnmKWkN4ID89shygbOY7adtv5TPnd2eJGkDSaMlTZI0U9KxpWkqXH+uNidJw/K+7l5Dev8TCKE+ko6QNFbS2Olvf9Dq4rS1/lrz3J30h+Efth+tM4+tgXOBicAYYFFgci/KNAm4sMKxFYBd8r9n9uIaIfSK7VHAKIDF1l7eLS5OW+uvwXOTvH2oF3kcDvwE+IHtWb0tkO1xwIjS/ZIWAe7Kb39h+/neXivM8/rlHWW76a8/5GXz9pVe5DEO+GEjAmclkgYCVwAbAVcDJ/TVtUK/Mi1vF6lwfNVmFWRe1l9rno1wne2+voX+DWmo1P3AQdUCtaQjgCMA5p9/8T4uVmhzE/N2ndIDkhYEhjW1NPOo/lrzfDVvl+tFHn06bEnS/wBHAs8Ae9TQqTXK9lDbQwcPWrgvixba3215e7CktQs7c+A8B2jtkpfziP4aPB/O2017kUdvuhqr/lwl7QOcCrwO7Gb71WrpQwBm35XYvhu4EVgMeETSLZJuBP4L7Axc0Joizlv6a/C8ifRlW1/SRn2Qf91tTpK2Av6Q89jT9r8bXLbQvyyYt6XjkAt/gF8Btgc2Jn3vNwYmNK1087B+GTxtP03qiAE4R9L81dJLWkPSij24RF1tTpLWBK4D5gcOsX1PD64Z5jGSBKye344vPmZ7iu0Tba9me7DtlWwfZnuS7ZG2ZXtkyTldef+IpnyAfq5fBs/s68CzwObA7ZLWL00gaWFJx5Fu85fvQd49bnOStDRwM7AM8F3bV/bgemHedCjp+/Iq8HiLyxJK9NvedttvSNqGNDHINsDfJf2LNARpGvARYDNSLfBl4I0e5H13bmPandTm9DdgBjCU1FxwAfClktOOJj0u+h7w8SqzeX/b9mu1liX0L5IWAn4HrAFsmXd/3/aM1pUqlNNvgydAfuRyW0m7kx7X3BL4DOlzvwr8lXQbfantnj7bvg/wQ2B/UpvTa6Q2p5NIveilBubtIqQaRSUjc15h3jQYOAh4lzSxzZm2R7e2SKGcfh08C2zfSOqdrCXtCMo8CVQm3RTgxPwqNTK/itPPtS+EUrbfAtTqcoTu9ec2zxBC6DMRPEMIoQ4RPEMIoQ7zRJtnCKHn5hswiyUX6M0sjL33Qhu3/kbNM4QQ6hDBM4QQ6hDBM4QQ6hDBM4QQ6hDBM4QQ6hDBM4QQ6hDBM4QQ6hDjPDvQ1KXF04cManUxqvr4T6Z1nyiEDhY1zxCaQNIASRMkWdKrktr7r1/oVgTPEJpjJ2CV/O9lgD1aWJbQABE8Q2iOL+ftxJL3oUNF8Ayhj0laCtgTMGny7JnALpJWamnBQq9E8Ayh7w0nLfcyJi8b/GfSygJzrSggaUhuFx1fKbN83BWObSjpOklvSHpf0sOSvtzdeaHnIniG0PcKt+hdeVtYV710natekbQ9cB+pPfVl4HrgHWCUpNMbea0QQ5VC6FOSNgI2JK1JdFXefT1pwcGPSdrW9t8acJ2FgIuBBYAfASNtOx/bCri1t9cIHxY1zxD6VqHWeaXtyQC2pwKXlBzvrS8CKwL/Bk4uBM58vXuB39aSiaQjJI2VNHbqm1MaVLT+KYJnCH1E0vzAgfntBSWHC+/3kbRIAy736by9wvasMscvrSUT26NsD7U9dP4lF2hAsfqvCJ4h9J29gKWA/9i+p/iA7UeBx4GFgf0acK2P5O1zFY5X2h/qFG2eIfSdwi354pLuLnN8uaJ0v68lQ0ndVXgq9aaXq42GXojgGUIfkLQKsGN+uxxzAmU5W0la2/ZTQGFSgEq38qtW2P9iN8eHVLl+qENb3LZLGp/HoA2rcHx9SS/lNFdJGtzkcg3pZT5dOZ8RjSlZ6AAjSL9ft9tWpRdwZU5fqKW+SgqgS0tatky+u1W43l15u2+F2ukB9X2MUElbBM9qJG0O3AmsQGpk3892TNkT2pYkkYInwEXdJC8cP1jSQNvTgcLQpZNzXoV8tyENQypnNGls5zrASSXnbQ58rUcfInSrrYNnHvT7V2BJ4FfAYbZnNrEIOwDrMud55BBqMQxYHfgAuLqbtLeQapsrMqdW+QNS7fMo4AlJoyU9SKpElB1yZPt94GBgKinAPiHpUkm3A/cA5+Wk0+v8TKFE2wZPSXsAfyK1/Zxs+1vFY9eawfYztsfl2kAItSrcgl9r+91qCW3PAC4vPi+Py9wBuI00E1MhqB5i+/tV8voLsBVwAykY70WqeBwN/DIne62HnyVU0JbBU9Jw0l/swcCxtkeWSbO5pNPzgN6XJU2T9GJuE92iQr4DJR0p6V5Jb+dzXpb0iKSfl7YxVWvzVLKvpJslvZLzmijpNknfqPLZ1sw1gpclTZU0TtJ3auhFDR3C9sG5TfPA7lOD7WNy+r2L9t1te0fbi9pe2PZmti/JxwrtpeXyesT2HraXtL2Q7Y1sjyLVhgEe7uXHC1nb9bZLOgo4mzS04su2uyok/QnpC/EE8CDpdmVt4AvAXpIOsD265JzfkyZj+AC4m/RXeBlgDeA4UrvRqzWUcXBOuwdphpz7gQnA8sB6wPbAWWVO3RA4M1/3DlIP7LbAqcDKQMWgG0J3JC0HLGj7uZL9WwCFZ9u7ml2u/qrdguexpKm7pgEH2L6mStozgOG2Xy7eKelzpFrruZJuKjwSJ2lVUuB8Hti0zHkbMme4R3dOIwXOfwN72h5XlM9A4LMVzvsmcDLwo8JTIJI+RQqkR0s6zfbzNZYhhFIbAH+R9E/gv6Tfo9WBjfLxi2x31wYbatRut4p75u2vugmc2L6lNADm/TeQaoVLAdsVHSqMs3ukwnmP2X6luwLmv+5HkWrGny8OnDmfmbavr3D6Q6T221lF6e8iTdowoKS8IfTUOOAc0ndpG9Lv06rA7cAhlJkCL9Sv3Wqed5Ke0T1e0sO2r6yWWNIywO6kW+UlmPN51svbtYCb8r/HkWa2+ayk7wGXlN7e1Gh7UlvsPbaf6OG5f6rQ6TUO2BWoODmupCOAIwAGLrVEDy8b5gW2XyB1DoUmaLfgOZJ0y/tt4BJJVAqgkr4K/AJYqEp+ixX+YfvdPCns+aT20p9Imkia//Am4HLbtUwjU3iCY1zVVOVNqLD/nbytOBNDbvQfBTD/kJVjQtsQWqzdbtuxfQKpPXM+UgDdtzSNpE1JtyeDgBNIA4MXAQbkXsj/KyQtyfsq4KOkAcznA++RpvK6ABiXH6nrtog9/1SzxfPFIfQTbRc8oWwA3ackyRdIgfHXts+w/ZTt94tuideskvdbti+0fZjtdXLaO0g1yp/VULxC7XHtHnykEEI/02637bPZPiE/YfZt4NJ8C18YerRU3s7VM53Hau7Ug+s8I+knpM6aT9Zwyu2kpzS2krSu7SdrvVYInWTGrAG8OaVaq1jfKz+atT2u35Y1z4KSGuilkr6YDxXaGw8pnkhW0qKk2/ElSvOStJGk/SQtWOZSn8vbbjuQco/8uaSf3dWS1iq5zsA8XCqE0I+1bc2zoKQGeln+9wWkMaEbA8/muRIFfIo0tu185l7eYFXSY3CTJT1CqrUOJo2BW53UE/+DGot1Amlg/W6kZ4jvA14gDYdaP29b/DczhNCX2rrmWVBSA72M9NzvUFLv83ukHvqhwDWkgFpuoPn9wImkqbtWJj33uyMwGfg5sL7tsTWWZyqptnpwzm89UsfTOsDfiRlsQuj32qLmaXtIDWlOINX4in21QvKR+VV8/iTSY5CnNqJceaD7xfnVXT4jmDNFWbnjIykpbwihvXVEzTOEENpNBM8QQqhDBM8QQqhDW7R5htDfSVqHNEJkO9IExyJNf/gC6RHhW/JkxqFDRPAMoY9J2g/4A2lo3ERgDPAmsCxpdMiWpAlxInh2kAieIfQhSSuQxh0PBr4FnFW8DldeQWCb/AodJIJnCH1rd9LMX/fZ/lXpwTzk7S7mLB0cOkR0GIXQtwqTcHc70XYpSQtL+h9JD0l6R9IHkp6QNLL4seSi9CPzmlsjJa0m6eK8VtaUfN7xkqLC1CARPEPoW4VZuHaQtF7VlEUkrUxam+tnpEeL7wP+TFoN84fAPZKWrHD6asBYUufUGNKsYauTntIbHYsNNkb8EEPoW9eR1sZaBHhU0k25NrmjpMXLnaA0gcOVwMeB3wCr2t7F9p6kORUuJq1X9Mty55OW3LgDWMP2frZ3zeknkh5LPrJhn24eFsEzhD6U123fkVQTnI80mczPSD3rb0i6J/fGF/sMqQf+fuCbtj8oyu8DUvB7BRheofY5GTi6eGUE2/8BCmu+f6sRn21eF+0fHWjQ22Llm9v7795N993Q6iJ0a+CKzblOnvN1U0lbkSax2Zw0RGlJYCvS3LC75jkQIAVYgKuLFwssyu99SWNzuk1Jt/PF/lJhMcNLgfOANSV9xPbE0gTFa2UtsPyiPfug85j2/g0MoR+xfa/tk2zvCCwDbM2cwHdo0YoJq+ft6bkDaK4XcwLssmUu9d8K158KvJTfrlwhzSjbQ20PHbR4ualvQ0HUPENogVyjvFfSbqSOoY1J7ZGjgYE52Z3A+G6yqmcF2NAAETxDaCHbMyXdTgqehVpkYT7a0bbPriPbIeV2ShoMFBor5rplDz0Tt+0h9KHcc96dj+btC3l7c96WLnxYq50lLVNm/wGk3/ln8hrvoRcieIbQt46WdIGkzUoPSJpP0ldIqxAAXJG31wIPA5+WdK6kpcqcu0I+t5yFgLMlzV+Ufg3glPz2zPo+SigWt+0h9K1BpFUERkiaBDwGvEFaAXYDYKWc7jTbt0JqD5W0F/An0moJB0p6nHQ7vwCwFmkM6CvA/ytzzYtIvfrPSLoHWJQ0YH4B4AagnqaAUCKCZwh96/ekTp8dgc2Ys0DgdNJt+oXAebbvLj7J9gu5tnoYsG8+b3PgdVJ75c+BP1a45rOkIUw/BbYHFs/7zgd+VW74U+i5CJ4h9KE8SP7a/OrpuVNItcQe1xRtPwvs39PzQu2izTOEEOoQwTOEEOoQwTOEEOoQbZ4h9BO2RwIjW1yMeUbH1DwljS/zjO8USRMkXSnp0w281pic/7AenteVzxvRqLKEENpTJ9Y8bwUm5X8vCWxEehJjH0nH2a40x2EIITRMJwbPU22PKbyRNIj0xMRRwKmSRrfw0bMTgVOZM3NNCKGf6pjb9kpsTweOB94lrVC4cwvL8pLtcbbfblUZQgjN0Yk1z7nY/kDSv4FNgOUL+yWNJ63/sprt8aXnSRpDWi97u+LabEma7YCTct6DSI/XnWb7+jJpu4BDgS/Z7io5JuAg0hMjnwQWJD0pcgupNv08IbSRqdPn4+kXy00X2jxLvuKWXn/A9CrHmleMPldYD+blBua5N/BX0uN0NwOPkyawvU7ScbVmkgPnxcAfSDOHP0R64kTA0cBjkjZtYLlDCH2sXwRPSZ8grRg4jbmXJOiNY4Dv2N7A9oG2twb2AGYAp0naoMZ8jgIOJAX2jWzvbHt/4GPAWaRJIkYXz4ITQmhvHR08JS0paVfgGtJnObbBnUVjbZ9RvMP2DaS1YAYC36gxn+Pz9vu2nyjKaybwbdLytKsyZ2qyuUg6QtJYSWOnT32vBx8hhNAXOjF43lG0jssbpGm7VgV2tX1Og691SYX9F+XtsO4yyOtvrw7MKjpvNtvTiq5TMb8PrS0z/yLdXTaE0Mc6scOoMM5TwArAp0jzFP5B0ta2n27gtcoupMWcdWXKLqJV4iN5+1LxUrAlni1JG0Joc50YPEvHea5ICqjrA5dI2sJ2rV10zax5t7bbMITQUJ142/4htl8iTRY7nTTZ7PCiw9PyttJ97qrdZD+km/21LKJVSLNSlQ6h1UvShhDaXMcHTwDb44Df5rcjJRVq1IVgtE7pOZLWA1bpJuvh3ewfU0PZXiDdlg8gjfMsLcegnuQX+gdJAyQdKOk6SRMlTZX0Ru4UPEXScq0uY6iuXwTP7Cekp4zWAA7O+27L2/+RtFghoaRVgC5Su2k1m0r6VvGOvM72QcBM4Dc1lu0XeXuKpNmBXNJA4DTS6onPAVfVmF/oYLkT8UFSR+HupNEW1wD3kobc/S9p/aF9G3CtEbmDtau3eYUP6zfB0/arQGFY0f/m2ufZpEWzNgWeknRNXiP7SeAd0pe1ml8DZ0h6XNKlku4GbiK1FZ9o+7Eai/db4DLSmtmPS7pV0mXAv4FjgTeBfWxPrTG/0KHySph/Iz2xNgZY0/aWtg+wvTupE/RE0gqYl0v6QssKG6rqN8Ez+wVpIPrqwKG23yQ9EXQJKeB9ltRDfjrwGVI7aTV/BHYhLbr1WdIMTvcCe9s+vdZC5Q6s4cAhwAOkhbw+T/r5nwN80vZDteYXOtrZpDbzh0jD6z40osP2dNunAseR7ozOr7AGe2ixjulttz2khjTvkf5yF+97njJtjdmwCvmU7v9rtwVM540gLTNb7phJ4zznGusZ5g157fTCrfjRVYauQbrrOYw0iuTrlExyLGkX0rLEWwDLkMY8PwNcD/w6z/cwnjmdoodKOrQoiwvz9zXUqb/VPENoZ7uTfueesD22WsL8x/YP+e0ehf1KziFNKLM3qVP0atK8C6uQpkQsTI5zFXBP/vczpGWOC68PLXUceq5jap4h9AOb5O2DNaYvNOV8UtJ8tmcA3wSOJDVP7WX7/kLiPAHNdqQ2dGx/O69qsDVwd9Q0GyuCZwjNU5jfrdaZvwrpBgBLSXqDND0iwIjiwAmza6u397qUoSZx2x5C+yodSjeU1L75gu1b+uSCRRPQzHzn/b64RL8RwTOE5nktb5evmmqOwkD5WaQOoULnz1ONLFSx4gloBi62cF9dpl+I4BlC8zyct1vUmH6zvH08t3fG/AhtJIJnCM1zI6kWuW53Kwfkzp9D8tsb8nZC3q7dN8ULPRHBM4QmydMlFh7BPVvSAlWSHwOsR3rk+Oy872HSrf/KeZxnLQqT40TncINF8Ayhub7GnEeG/yRpSPFBSYMkfYf0tJyBw22/ArNXiv2/nPQCSZuVnCtJ20lavGh3YXKcdRv+SeZx8dcohCay/ZqkbYDrSGMyn5b0AGlimEVJCwQuBbwPfMX2lSVZ/JIUCA8H7pc0Fng6n/Nx0kD51YDC8tf3kyYP3zinfYL0WPI9ti/osw86D4jgGUKT2Z4gaShwALAfafD8pqSA+SxpvoPf2J5U5lwDX5F0HWmw/GbAhqT5F54mLSg4qSj9VEmfIc06tiVpfoYBpN/9CJ69EMGzA01fGCZt3t4tLp/d5DOtLkINap1RsPHy4n8X51c9599I6oCqJe3jpEdDQwO1929gCCG0qQieIYRQhwieIYRQhwieIYRQhwieIYRQhwieIYRQhwieIYRQhxjnGUIoa/DrYpVLWxsi5r/5vpZefz5XntM0ap4hhFCHCJ4hhFCHCJ4hhFCHCJ4hhFCHtgieklzHq6vV5Q4hzLvapbf9wjL7VgB2IU3TdVWZ43f3aYlCCKGKtgietkeU7pM0jBQ8Xyt3PIQQWqktbttDCKHTdGTwlNSV2z1HSNpA0mhJkyTNlHRsaZoKeYzMx0dW2i9p5ZzPS5ImS3pE0heL0m4t6U+SXs/H7yi3KqKkITnP8ZLmk/RdSU9KmiLpZUkXSvpoQ39IoeXy/3ehjb7iZMSS/pnTDGte6UJvdWTwLLI18CCwMTAGuAWY3KC8h5BWK9wWuBN4hLSEwZWS9pe0N3AHsAzwF9IaNMOAOyStVSXfK4CTScvIXgtMJS0x+5CkWFK2//qppE7/fQtF2qLNsxcOJ63N8gPbsxqc96HAmcDxeckEJB0F/BY4HVgYGG57dD42ALiUtCbNd4DDyuS5KrAgsJHtf+XzBgO/Bw4CLiKtSRP6l8nA+sBw0v9x6Ac6/S/hOOCHfRA4AcYD/1MInNko0kJbKwO3FAInQC7Dz/Lb7arke0ohcObzpgHfAN4BNpW0dbmTJB0haayksbPer/y8bWhLv87bk/Mfy9APdHrwvK4kuDXSHTmwzZavNT6/vaXMOf/J25Wq5DvXgl+23wJuyG+HlTvJ9ijbQ20PHbDwwlWyD23oalLz0mqkFS9rktdw/7qkByS9I+mD3FZ+qqSlS9J+JrebPlolv6UkTc2vper+NAHo/OD5XB/m/UKF/e9VOm67cGz+Cue+lQNlOePzduVaChc6zol5e5KkRbpLLGkB4M+kpYTXA+4i/YFdgtQs9LCk1YtO+QvwIrChpA0qZHsAMBi4wfYb9XyIMEenB88PenFud5+9u6aAvmgqCP2U7dtJwXA54PgaTvkR6S5kHLCW7d1t7wusTqrJrgpcUpT/TOAP+e2ICnkemrddPSt9KKfTg2c1hVvuSn/lV21WQYosIWnxCseG5O3EJpUlNN+JgIHjJS1bKZGkBYGj8ttjbM/+Ttj+gHTr/x6wRUkbeVfeDpf0oc5gSR8HNgUmUb7JqZBudtv69GnRtl5Nfw6ehS/cOqUH8pdzWFNLM8fw0h05oBbGAY5pamlC09h+BLgSWBQ4qUrSTUh/9F+0/Zcy+bxGmTZy208B95Nqt7uVnFaodV5ie0aVMs5uWx80ONrWq+nPwfO2vD24ePxkDpznAK0alP4DSesWlWcQaUjU4sDDtuOZ/f7tf4EZwJGSKt39fCRv/1sln2dL0hZckLeFYImkgaShcBC37A3Tb4NnDkI3AosBj0i6RdKNpC/kzsz5kjXTBOBe4DFJN0u6HHiG9EV/jTRYPvRjtp8GziN1Kv6ou+R1XOIKUl/A7kU98juSRoA8bPufdeQZyui3wTPbBzgVeAXYnvQk0k15O6EF5TGwL3AKqeF/L2AB0vClTYvHf4Z+7UekgfMHSVqvzPFCk9NqVfIo9LR/qI3c9tukJ9cGAwfm3SPytqvnRQ2VtG3wtD3GtmwPKXNsRD7W1U0eU2yfaHs124Ntr2T7MNuTbI/MeYwsOafs/qLjw/LxMRWOy7aqlGmG7R/bXtv2AraXs32w7fHVPkvoP2y/RGqqGQD8tEySh0kdQh+RtEPpwVyj/Fx+O6bM+bNv3XN7+l6kDtRLe1Xw8CFtGzxD6OdOA94kBcEP1TBzj/q5+e2ZklYsHMvjP88hdSjdb/ueMnnfBjxP6nj6EenuJsZ2NlgEzxBaID8s8X/57UJlknyfVKv8BPAfSddLuoLUUbQPqdlprpEbOe9ZzHmG/pi87WpEucMcETxDaJ2zqPAkm+0ppI7NY4B/keZL2JM0B8JpwMa2ny13btZV9O+qYztDfTp9VqWOkNszK7aDhv6pXHt9yfEpwCpVjk8nBdiz6rj2f4jvXJ+KmmcIIdQhgmcIIdQhgmcIIdQhgmcIIdQhgmcIIdQhgmcIIdQhgmcIIdRBdj0Tt4RWkvQqjV+CZBnSzE7tqi/Kt6rtipMSz+sa8D1rh+9Ub8tQ8TsSwTMAIGms7aGtLkcl7V6+MLd2+D/ryzLEbXsIIdQhgmcIIdQhgmcoGNXqAnSj3csX5tYO/2d9VoZo8wwhhDpEzTOEEOoQwXMeJGltSd+UdLGkcZJmSbKkL7a6bJBWFJW0g6Sf5zXE35E0TdJESVdJGtbqMoa5tfp71ezvTcznOW86CvhmqwtRxaeBwnrlk4C7gPeBjwNfAL4g6RTbP2hR+UJ5rf5eNfV7EzXPedM/gdOB/YA1gTtbW5y5zAKuBj5le0Xbu9vez/b6wP7ATOD7krZraSlDqVZ/r5r6vYma5zzI9nnF76X2mnDc9u3A7RWOXSFpJ+Aw4CDgjmaWLVTW6u9Vs783UfMMnejRvF25paUInaah35sInqETfSxvX2ppKUKnaej3JoJn6CiSVgBG5LdXt7AooYP0xfcmgmfoGJLmAy4GFgdus31Di4sUOkBffW8ieIZOci6wA/A8qdE/hFr0yfcmgmfoCJLOJPWUTgJ2sD2pxUUKHaAvvzcRPEPbk/Rz4BjgVdIvwH9aXKTQAfr6exPBM7Q1SacBxwGvAzva/leLixQ6QDO+NxE8Q9uSdCpwAvAmsJPtv7e4SKEDNOt7E8EztCVJPwa+A7xF+gV4tPoZITT3exPzec6DJG0M/LZo18eBRYH/AG8UdtreoslFA0DSHsB1+e1Y4IkKScfZPrU5pQrdafX3qtnfm3i2fd60GLB5mf0fK7OvFZYq+vfQ/CrnTiCCZ/to9feqqd+bqHmGEEIdos0zhBDqEMEzhBDqEMEzhBDqEMEzhBDqEMEzhBDqEMEzhBDqEMEzhBDqEMEzhBDqEMEzhBDqEMEzhBDq8P8BUTSPjbhRfpUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 345.6x360 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "w = calculate_W(np.dot(As[0][0],As[0][1]), np.dot(As[1][0],As[1][1]), np.dot(As[2][0],As[2][1]), Ss[0][1], Ss[1][1], Ss[2][1])\n",
    "\n",
    "fig, axs = plt.subplots(1, 2, constrained_layout=True, figsize=(4.8,5))\n",
    "color = \"viridis\"\n",
    "\n",
    "A_new = np.dot(As[0][0],As[0][1]) @ w[:,:,0].T\n",
    "B_new = np.dot(As[1][0],As[1][1]) @ 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": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW8AAADOCAYAAADxC4PrAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAOfElEQVR4nO3de6xlZX3G8e8DDAwDBUSxohPBSkRbJV6mUaxVSWybaquNiqHNaAngxApNS4024g2BSmpbNBEtjvVCKTaKVeIFEW8oKaV2qoRqxds4IFVuKoiD4Ci//rHW6GGzZ846M3vPOi98P8nOmr3WnrOe7HPOc9Z+17v2TlUhSWrLbmMHkCQtneUtSQ2yvCWpQZa3JDXI8pakBu2xq3a0Z/aqleyzq3Z3r7flQW08lyuu3zx2hHuNJGNHGMQZbLN1Gz+8uaoOmly/y8p7JfvwxN2esat2t+Ma+cG77rgnjx1hkNVnXj52hHuN3VauHDvCIHfdccfYEYZp5I/hp+664Jpp6x02kaQGWd6S1CDLW5IaZHlLUoMsb0lqkOUtSQ2yvCWpQZa3JDXI8pakBlnektQgy1uSGmR5S1KDLG9JapDlLUkNGlTeSV6Z5IIkG5NUkk1zziVJ2o6h7+f9BuAHwBeBA+aWRpI0yNDyfnhVbQRI8mVg3/lFkiQtZtCwydbiliQtD56wlKQGzfUzLJOsA9YBrGTVPHclSfcpcz3yrqr1VbWmqtasYK957kqS7lMcNpGkBlnektQgy1uSGmR5S1KDBs02SfJC4JD+7kHAnkle3d+/pqrOm0c4SdJ0Q6cKHg88bWLd6f3yc4DlLUm70KDyrqqnzzmHJGkJHPOWpAZZ3pLUIMtbkhpkeUtSgyxvSWqQ5S1JDbK8JalBlrckNcjylqQGWd6S1CDLW5IaZHlLUoPm+gHEC9V+q9hy5BN21e522IpLNowdYZDVZ14+doRBfnDckWNHWNSB7/qPsSMMkkNWjx1hkBuf8cCxIwyy/6YtY0cY5qMXTF3tkbckNcjylqQGWd6S1CDLW5IaZHlLUoMsb0lqkOUtSQ2yvCWpQZa3JDXI8pakBlnektQgy1uSGmR5S1KDLG9JapDlLUkNGlTeSR6R5LQkVyS5KcltSa5M8qok+8w7pCTp7oYeeR8HnAx8CzgNeDnwNeAM4PIke88nniRpmqGfpPMB4MyqunXBunOSfAN4FXA8cPasw0mSpht05F1VGyaKe6v39ctHzy6SJGkxO3vCcuuH6t2ws0EkScPtcHkn2R14DfAz4L3beMy6JBuSbNjy0807uitJ0oSdOfJ+M3Ak8Nqq+tq0B1TV+qpaU1VrVuzppBRJmpUdKu8kpwMnAeur6szZRpIkLWbJ5Z3kVODVwLuBl8w6kCRpcUsq7764XwecC5xQVTWPUJKk7Rtc3kleS1fc5wHHVdVdc0slSdquQRfpJDkReD1wLfAp4E+SLHzIDVX1ydnHkyRNM/QKy9/slw+lGzKZ9DnA8pakXWToFZbHVlW2c3v6nHNKkhbwLWElqUGWtyQ1yPKWpAZZ3pLUIMtbkhpkeUtSgyxvSWqQ5S1JDbK8JalBlrckNcjylqQGWd6S1KCh7yq40/Kj21lxyYZdtTstEwd97FtjR1jURd+9cuwIgzzziPuPHWGQB771m2NHuE/wyFuSGmR5S1KDLG9JapDlLUkNsrwlqUGWtyQ1yPKWpAZZ3pLUIMtbkhpkeUtSgyxvSWqQ5S1JDbK8JalBlrckNcjylqQGDSrvJIcnOT/JV5PcmuT2JFcnOSvJwfMOKUm6u6EfxrAaOBj4EHAd8DPgMcA64Jgkj62qG+cTUZI0aVB5V9WngU9Prk/yeeD9wLHAG2eaTJK0TTs75n1Nv7zfzgaRJA23pM+wTLIS2BdYCfw68Lf9potmnEuStB1LPfI+AbgJ+A7wCeAAYG1VXTbtwUnWJdmQZMMW7typoJKkX1rqp8dfCFxNd/T9OODZwAO29eCqWg+sB9gvB9aORZQkTVpSeVfVdXSzTQAuTPJvwH8lWVVVZ848nSRpqp06YVlVVwFfAl46mziSpCFmcYXl3sCBM/g6kqSBhl5h+aBtrD8KeDRwxSxDSZK2b+iY9z/2l8F/hm5u90rgCcAxwG3Ay+YTT5I0zdDy/lfgRcALgYOAoivxtwN/V1XXzieeJGmaoZfHv5/uMnhJ0jLgW8JKUoMsb0lqkOUtSQ2yvCWpQZa3JDXI8pakBlnektQgy1uSGmR5S1KDLG9JapDlLUkNsrwlqUFL/QzLe73f/8otY0cY5OO/ccDYEQb5+Q03jh1hUb/34MeOHWGQzc8/bOwIg+zzge+PHWGQ3fbZZ+wIw/x4+mqPvCWpQZa3JDXI8pakBlnektQgy1uSGmR5S1KDLG9JapDlLUkNsrwlqUGWtyQ1yPKWpAZZ3pLUIMtbkhpkeUtSgyxvSWrQDpV3klVJNiapJGfPOpQkaft29Mj7NOCgWQaRJA235PJO8njgL4HXzTyNJGmQJZV3kt2BdwAXAx+cSyJJ0qKW+hmWJwOPBJ43hyySpIEGH3kneRjweuC0qto08P+sS7IhyYYt3LmDESVJk5YybHIOsBE4a+h/qKr1VbWmqtasYK8lh5MkTTdo2CTJWuB3gKdW1Zb5RpIkLWbR8k6yF93R9kXA9UkO6zc9pF/u36+7uapumUtKSdLdDBk22ZtuTvezgG8suF3ab1/b3z9hDvkkSVMMGTbZDBw9Zf1BwNvopg2+E7hqhrkkSduxaHn3Y9wfmFyf5ND+n9+qqntslyTNj29MJUkNWupFOr/Qz/XO7KJIkobyyFuSGmR5S1KDLG9JapDlLUkNsrwlqUGWtyQ1yPKWpAZZ3pLUIMtbkhpkeUtSgyxvSWqQ5S1JDUpV7ZodJTcB18z4yz4AuHnGX3PWWsgI5pw1c87WfTnnIVV10OTKXVbe85BkQ1WtGTvH9rSQEcw5a+acLXPek8MmktQgy1uSGtR6ea8fO8AALWQEc86aOWfLnBOaHvOWpPuq1o+8Jek+yfKWpAZZ3pLUoKbKO8krk1yQZGOSSrJp7EyTkjwiyWlJrkhyU5LbklyZ5FVJ9hk731ZJDk9yfpKvJrk1ye1Jrk5yVpKDx863LUlWLfj+nz12noX6TNNuPx4720JJDkzy90m+meSO/uf0s0l+e+xsAElO3c5zWUm2jJ1xqyT7Jjklyf/0v+s3J7k8ybFJMs997zHPLz4HbwB+AHwROGDcKNt0HHAi8GHgfGALcBRwBvCCJE+qqp+MmG+r1cDBwIeA64CfAY8B1gHHJHlsVd04Yr5tOQ24x9Vmy8hl3HPGwXIqm0OAS4F9gXcCXwf2B44AHjJesrv5IPDNKeuPAF4OfGTXxpkuyW7Ax4EnA+cCbwFWAX8MvBt4FPDXcwtQVc3cgF9b8O8vA5vGzjQl4xpg/ynrzwAKOGnsjIvkP7rP+Yqxs0zJ9ni6PzJ/1Wc8e+xME/kKeM/YORbJeBnwHeDgsbPsQPa398/xs8bO0uc5ss/zpon1ewIbgVvmuf+mhk2qauPYGRZTVRuq6tYpm97XLx+9K/PsgK3vP3O/UVNMSLI78A7gYrojs2UryZ5J9h07x6QkTwWeAryxqr6XZEWSVWPnGqIfcjyG7lXixSPH2Wq/fvndhSur6qd072+yeZ47b6q8G7e6X94waooJSVYmeUCS1Ul+l+7oBuCiMXNNcTLwSOCksYMs4vnA7cBtSW5M8pYk+48dqvfMfnltko8APwE2J/l6krUj5hriaLqyfE9V/XzsML0vALcAr0hydJKHJnlkkjOBJwCnznPnrY15N6k/anwN3Uv+944cZ9IJdGN1W20C1lbVZePEuackDwNeD5xWVZuSHDpypG35AnAB3XjtfnRleRLwtCRPrqqxT1we3i/fAXwD+FO6l/gvA85LsqKq3j1WuEUcTzdE8a6xg2xVVT9M8mzgn4D3L9h0G/C8qrpwnvu3vHeNN9ONj51SVV8bOcukC4Gr6U5gPQ54Nt3bWi4n59CNIZ41dpDtqaonTqz65yRXAX8D/EW/HNOv9MvbgKP6l/ckuZDu+X1DknOr6q6R8k2V5HC64Z5PV9W3x84z4cd0598+DFwOHEg3YeG9SZ5TVZ+c257HHvTfiZMFy/KE5ZScp9MdMbx97CwD8x4B3Am8cuwsfZ61wF3AUxasO5RleMJyG/lX9M/n5csgy0f65+2MKdvO7bc9auycU7K9sc92zNhZJnI9hm7o6SUT61fR/THcBOw+r/075j1HSU4FXk03begl46YZpqquAr4EvHTsLEn2ojvavgi4PslhSQ4DDukfsn+/7oCxMi6mqrbQndBaDq9mruuX10/Z9r1+udxOVO8BvAj4Pt201uXkZGAl3VDZL1TV7cDH6H5OD53Xzi3vOemL+3V0RzQnVP8nuRF70738G9vedHO6n0U3Rrv1dmm/fW1//4Qxwg2RZCXdyerlcKL6C/1y9ZRtW9ctt7n9fwj8KvAvVXXn2GEmbJ0Xv/uUbXtMLGfO8p6DJK+lK+7zgONqmY0hAiR50DbWH0U3nfGKXZtoqs10swwmb1tfFVzc3//wKOkWSHL/bWw6ne4XeDlcWHIh3Xj32oVTGfsrav8I+HpVTbs4ZkzH98t3jppiuv/tl8cuXNm/EnwO8EOmX2w0E029JWySF/LLl8x/Tnem/B/6+9dU1XmjBFsgyYnA2cC1dDNMJov7hprnSYyBknyI7grLz9DN7V5JN73pGLqpbk+vqitHC7gd/WyTbwNvraplMXUwyZuAJwGfpfve70s32+Qo4D/pThCOfmVtknV000G/QjdzY0/gz+h+Fv6gqi4ZMd7dJHkw3XP533XPk8Gj669W/SLdUNP5wL/TvWJ9Md1wyYlV9ba5BRh70H+JJwgupTtxMe126dj5+ozv2U7G5ZTzBcBH6a62u4PuxMvVdNMGHzp2vkWyH8oyO2FJd6T1CeD/+udzM3AlcAqwcux8E1mfS/fKajPdkfglwG+NnWtKzlP67/OLx86ynYwPpxsavY7ubRB+BHweeO68993UkbckqeOYtyQ1yPKWpAZZ3pLUIMtbkhpkeUtSgyxvSWqQ5S1JDbK8JalBlrckNej/AVR0r/qWuY4rAAAAAElFTkSuQmCC\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(Ss[2][0])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAB8CAYAAABwrOvwAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAALeklEQVR4nO3de4xcZR3G8e/TC5RSsbYUELFVi5FL4yWtmCgqoBSEqERRQUGrEDCNISFiAA0XNfEKppGIFUVr5CIiBDRI0ApVCBJFURQvBaSAiFKkrQja0vbnH++7Oiwzu2frnHnnleeTTM7uOdOdp2dnnj17zvvOKiIwM7P6TCodwMzMto0L3MysUi5wM7NKucDNzCrlAjczq5QL3MysUlMG9UDbafuYxo6DejgbEpo6tXSEcb1w7w2lIzRy5293Kh2hkdiyuXSEZioZQf0o6x6OiDndtg2swKexI6+YvHhQD7fttm4pneD/ypTdnlM6wriuue6a0hEaOewlB5eO0MjWdetKR2gkttbR4Cu3XHZvr20+hWJmVikXuJlZpVzgZmaVcoGbmVXKBW5mVikXuJlZpVzgZmaVcoGbmVXKBW5mVikXuJlZpVzgZmaVcoGbmVXKBW5mVqlGBS7pWElflnSbpE2SQtIB7UYzM7OxNH072Y8D84CH8m343yPUzOz/XNNTKMcBcyNiV+CbLeYxM7OGGh2BR8QP2w5iZmYT44uYZmaVcoGbmVWq1b+JKekE4ASAaUxv86HMzJ52Wj0Cj4gLImJRRCyayvZtPpSZ2dOOT6GYmVXKBW5mVikXuJlZpVzgZmaVajQKRdLxwP7500V5eZqkJfnjr0TETX3OZmZmY2g6jHB/4D2j1h3S8fEqwAVuZjZATafSLwGWtJrEzMwmxOfAzcwq5QI3M6uUC9zMrFIucDOzSrnAzcwq5QI3M6uUC9zMrFIucDOzSrnAzcwq5QI3M6uUC9zMrFIucDOzSrX6R42fREJTB/dw2+p9v76rdIRGLnzRC0pHaGTrI+tKRxjXIbu/tHSERv7+zvmlIzQy8+pfl47QyKQdppWO0Mza3pt8BG5mVikXuJlZpVzgZmaVcoGbmVXKBW5mVikXuJlZpVzgZmaVcoGbmVXKBW5mVikXuJlZpVzgZmaVcoGbmVXKBW5mVqlGBS5pD0knS1op6X5JmyQ9IOkSSQvaDmlmZk/V9Aj8A8DngLnAtfnj24CjgFslHdhOPDMz66XpG3T/FHhNRNzYuVLS24BvAecDe/c5m5mZjaFRgUfElT3WXy5pNbCXpJ0j4uG+pjMzs576cRFzU15u7sPXMjOzhv6nApe0EFgA3BoR6/uSyMzMGtnmApc0A/g6EMCpfUtkZmaNbNNfGZa0HXA5sC9wVkRc3+N+JwAnAExj+rZmNDOzLiZ8BC5pCnAZcChwbkR8rNd9I+KCiFgUEYumqpK/AG1mVokJFXgu70uBI4DzIuKUNkKZmdn4Ghe4pMnAN4AjgeURcVJrqczMbFxNp9JPAlaQZl5eCCxtMZOZmTXQ9CLmWcAxwHrgz8BZkkbfZ5mHEpqZDU7TAp+XlzOBM3rcZwWp4M3MbAAanUKJiCURoXFua1rOamZmHfx+4GZmlXKBm5lVygVuZlYpF7iZWaVc4GZmlXKBm5lVygVuZlYpF7iZWaVc4GZmlXKBm5lVygVuZlYpF7iZWaVc4GZmlVJEDOaBpLXAvX3+sjsDD/f5a7bBOfvLOfurhpw1ZIR2cs6LiDndNgyswNsg6daIWFQ6x3ics7+cs79qyFlDRhh8Tp9CMTOrlAvczKxStRf4BaUDNOSc/eWc/VVDzhoywoBzVn0O3Mzs6az2I3Azs6ctF7iZWaWqKnBJx0r6sqTbJG2SFJIOKJ2rk6Q9JJ0saaWk+3POByRdImlB6XwjJM2U9HlJt0j6q6SNku6TdI2kg0rnG4+kq/L3f2jGBuc8vW7Hl87XScm7Jd0oaYOkf0i6Q9L5pbMBSDp7nP0Zko4dgpxTJb1f0s8krZO0PvfTByXt0Prj13QOXNIaYB7wEPAE8BzgwIhYVTDWk0j6FHAqcCewCngEWAAcBmwC3hARNxQLmEnaE/gl8BPgLmAdsDtwBPBM4EMRcU6pfGORdDRwEWl/PhYROxeOBKQCJ01WW9Fl83ci4heDTdSdpMnAN4CjgdtIz9MtwAuA1w7D/swHZgd02wScDkwG5kbEA4NL1SWM9B3gjcAdwMq8+mBgH+DHpH7a2lqAiKjmBrwOeG7++BwggANK5xqV8S3Aq7usf1vO+7vSGXOeycCULuufDfwFeByYXjpnl3y7AGuBZcAa4OHSmTqyBbCqdI4GOU/LWT/YZdtTnhPDdAMOzNmvHYIsr8hZrgcmdayfTPqh2Ho/VXUKJSJ+GBH3l84xloi4MiJu7LL+cmA1sJek4kc4EbElIjZ3Wf8gcDOwA6nMh80XgMeAj5QOUiNJO5KOYFdFxLmjt3d7TgyZJXn5tZIhsufn5fej4yg7IrYA1+VPW32tT2nzi9tTbMrLoX2RSJoF7Af8HRiqH5aS3gocCRwaEY9JKh2pm2dJOhGYTTrVd0NE3F04U6fFwE7AFZJ2At4M7AH8mXRU+1DJcGORNAN4K+l039WF4wD8Ni8XS/rMSInnU1SHABuBW9oM4AIfEEkLSefCb42I9YXj/IekXYClpAvazwbeBMwC3hcRm8b6t4MkaTbp6PuiiLhuvPsX9GJgecfnIemrwNIh2Z8L8/JZwB+A3Tq2PSbpxIi4ePCxGnk7sCOwIiI2lg4TEbfni75Lgdsl/SBvWkzar8dExJ/aDlHljSE9B94j6wzgN8BW4KDSeUZlW5D348jtUeBdpXN1yXkx6Yh2dse6NQzXOfDPAC8HZpIK8mDg53m/nl86X864POfZDHwXeBHpovU7SEe2TwAvLZ2zR/Yf5+wLS2cZlevDpIvAI6+hrcAXgV1bf+zS//n/YadVUeDAdsC1OeuZpfOMkXMKMB/4ZH4CfrZ0po5sb8z7752j1g9VgffI/kzgvlyYuw1BngvyvnwA2GHUthPztgtL5+ySe35+Xt5eOktHpkmkc/EbgPcCc0i/vR5NutB+DzCrzQxVXcSsjaQpwGXAocC5EfGxwpF6iojNEXF3RJxOOlVxiqTXls6VL7otB74XEZeUzjNREbEB+DZpZMJ+heNAKhuAlRHxz1HbvpuXCxk+S0hDCIfh4uWI40i5PhwRX4uItRHxSERcCpwEPA84uc0ALvCW5PK+lDSu+ryIOKVsogkZOZf3mqIpkjmk8emHjZ7IQZoTMDt/vr5oyrGNTDaaXjRFsjovN3TZNrKu9QkoE6F0tfrdpNM7FxWO0+nQvPxRl22r8vJlbQbwRcwWdEyUOBJYHhEnFY40Ubvn5TCMlnkUuLDHtncAU0kv6scHlmjiRo68+/0XqbbFqrzcu8u2kXX3DSZKYwcBc4GrImJt6TAdts/LbkMFR9a1erHVBd5nkiaRZuIdRSqepUUD9ZCn9a+OUSMjJD2XNE4Y/juWtZiI+BvQdRq6pNcDMyKi+DR1SfsAd3XZn8eQhurdA/ysRLZOEXGnpFXA6yQdGHlWsKSpwNn5blcUitfLe/NymE6fQJovcThwqqSbR773+QDujHyfbkfnfVPbVPrjgf3zp4uAfUkl85e87isRcVOJbCMkfRQ4E1gPnEe68DLasig8lFDSMtKvpTeRLgZuJE2lPpx0ZPHpiDitVL4m8lsrzIjhmPq9DDiGNFJi5Ah2EfAq4B+kt1Ao+twcIWkvUvnMIJX1g6RZzi8GbgAWx5BM6Mlj1R8k/Sa2x7DkgvR+QqQfynuS3jrj+6TRKK8nTaX/FfDKiGjvt8PSV3IneNV3BU8e8jb6tqSCjAE8bwhy7k86ovk9adLOE6TJHFeRyqb497vB/2ENQzIKhXQ+9Ergj6SZohuBu4EvAfNL5+uSdz7pGs3anHU1cBawfelso3Ien18z55TO0iPfLNKIuD/k/fgv4HfAJ4BntP34VR2Bm5nZf3kUiplZpVzgZmaVcoGbmVXKBW5mVikXuJlZpVzgZmaVcoGbmVXKBW5mVikXuJlZpVzgZmaV+jcnxi5XGrhwegAAAABJRU5ErkJggg==\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(Ss[2][1])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Reconstruction loss: 7.446640570183208\n",
      "Relative reconstruction loss: 0.8325597259468726\n",
      "\n",
      "Reconstruction loss: 8.166827191042334\n",
      "Relative reconstruction loss: 0.9130790379832163\n",
      "\n",
      "Reconstruction loss: 8.728341867521513\n",
      "Relative reconstruction loss: 0.9758582873317787\n",
      "\n"
     ]
    }
   ],
   "source": [
    "approx = outer_product_np(Xs[0], Xs[1], Xs[2])\n",
    "print(\"Reconstruction loss:\", np.linalg.norm(np.ndarray.flatten(X-approx), 2))\n",
    "print(\"Relative reconstruction loss:\", np.linalg.norm(np.ndarray.flatten(X-approx), 2)  / np.linalg.norm(np.ndarray.flatten(X), 2))\n",
    "print()\n",
    "\n",
    "X1_approx = As[0][0] @ Ss[0][0]\n",
    "X2_approx = As[1][0] @ Ss[1][0]\n",
    "X3_approx = As[2][0] @ Ss[2][0]\n",
    "\n",
    "approx = outer_product_np(X1_approx, X2_approx, X3_approx)\n",
    "print(\"Reconstruction loss:\", np.linalg.norm(np.ndarray.flatten(X-approx), 2))\n",
    "print(\"Relative reconstruction loss:\", np.linalg.norm(np.ndarray.flatten(X-approx), 2)  / np.linalg.norm(np.ndarray.flatten(X), 2))\n",
    "print()\n",
    "\n",
    "X1_approx = As[0][0] @ As[0][1] @ Ss[0][1]\n",
    "X2_approx = As[1][0] @ As[1][1] @ Ss[1][1]\n",
    "X3_approx = As[2][0] @ As[2][1] @ Ss[2][1]\n",
    "\n",
    "approx = outer_product_np(X1_approx, X2_approx, X3_approx)\n",
    "print(\"Reconstruction loss:\", np.linalg.norm(np.ndarray.flatten(X-approx), 2))\n",
    "print(\"Relative reconstruction loss:\", np.linalg.norm(np.ndarray.flatten(X-approx), 2)  / np.linalg.norm(np.ndarray.flatten(X), 2))\n",
    "print()\n"
   ]
  }
 ],
 "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
}