Neural_NCPD_Synthetic_high_noise.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Details\n",
    "\n",
    "The notebook contains the synthetic data experiment with low noise (sigma^2 = 0.05). Due to the non-convexity of the HNCPD 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": [
    "# loading packages and functions\n",
    "import sys\n",
    "import itertools\n",
    "from utils import *\n",
    "\n",
    "sys.path.append(\"./src\")\n",
    "import torch\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "from NNCPD import NNCPD, weights_H, Energy_Loss_Tensor, Recon_Loss, L21_Norm, outer_product, outer_product_np, PTF, random_NNCPD, Fro_Norm\n",
    "from lsqnonneg_module import LsqNonneg\n",
    "from trainNNCPD import train\n",
    "#\n",
    "import torch.nn as nn\n",
    "from torch.autograd import Variable\n",
    "from writer import Writer\n",
    "\n",
    "from sklearn.decomposition import NMF\n",
    "\n",
    "import tensorly as tl\n",
    "from tensorly import unfold as tl_unfold\n",
    "from tensorly.decomposition import parafac, non_negative_parafac"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generate Data Tensor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "## set the network parameters\n",
    "torch.set_default_tensor_type(torch.DoubleTensor)\n",
    "\n",
    "X = get_synthetic_tensor(seed=14, noise=0.4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "X_original = get_synthetic_tensor(noise=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualize Original Data Tensor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA2gAAAEkCAYAAABaExIDAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABLVklEQVR4nO3deaxlyWHf97r7fu+7b3+9TXdP9zRnyKFIigtISpQMWoQlJ1YQ+A/ZAZzAgcHQSRAjNhDACBA4CBAgDiI4QCwrWmzQcBQjcmJDW0iRFiEyIilyKHL2np6e3vvty3133/OHoMQJ6vebfq/J4WnO9/PnLdY5dc6pqnOqH+dXqfl8HgAAAAAAP3zpH3YDAAAAAAB/igUaAAAAACQECzQAAAAASAgWaAAAAACQECzQAAAAACAhWKABAAAAQEJkT/I/Xl7MzC+ez/2g2vL/sT/TTRuasnlInarseK+i65llrGmKlJqZQrPrwTynC9NDfW3mssMsb6qNTVvEdeePdBvnGd2Q9ETflFnWXZsum2V0NUcdMTU53ZYU7rp9RXdMU8+cLm2e6Wna4c6l2jhsH4TJoHvKm/LOzkPvtJf3V2TZ3A0D11ey8cLU2D08U1TQhSkzD9n5a2TGcclMmOP45Jwd6CrBHC41M/NX2sxfU13PzV9pM6dM8/F6maGpU9TnSk1kkb1u17ccO2+LY87MsD7V3BV0+4fdgzAeMg/9sL28q+c81/dm+dO9mNw4cO9V1//U3JzrmnkhZ9po5ij3zenmhrkZj6e5kykz57kDunlBfVeG4J+bnX/dfZ6KdrhvKHOucVlXzIj+476HTrS8uHg+F/74C+dPUuXU/ll7SZa9OViTZWMzwlzZF3/tE7JsVJNFYdR0C5J4Wda8E1ITXTY4o3to7Q39KKdFWRT6Z0UPDSEUdvT9Gi3EZ5CLv6VnsVFDt7G4O5Rlw2W9inSDvb+oV9ZqYIagXwqlA1PJTEiD5ulWiu7jbVwxH7RmQqpsqRlJt8N+/Jt/vBjW44Wv/6tf1JUewTs5D73Tnvknn5NlM7Mwcv9AM16Lj8n8pv64zAzMPHRVr34KN/Vk4+av8l3zD3Lv68my1P1S9PeF12WVkO3r+5jr6S+jUVV39kJL1xss6vFfPNRzSudMvF79jr6PR1f0My3u6zbm+rosPTrdotXNe+ojp7uh73FlU7fRzUPZQfxcL37pH+hKj+BHeR56Jz33P/1NWeYWKr2n9DhITXW/zB3qzjKp6L5e3NH11D90r39Tfw/1VvWclzXjcWC+a9zcMFzQ49GNH/UM8m3zPWSe29DMC/0l3ZDigT6om397q/p8agHt/qHInWvng/qZVu/Hz+W+h/i/OAIAAABAQrBAAwAAAICEYIEGAAAAAAnBAg0AAAAAEuIUGYTvjD9sPSPL6jai63Ra7zH/waP5D06DCAIJIYTCbvw/TpyZNLOU+481TdLZaEEfc3h+JMuy2zqAY+rS2sR/Jzkt6gton9P/seakqMMFSrv6P7QdLOn/mrN+R9c7vqjr5Vvx6x6X9LW5/8h+VNPPrXHb/Ee9ImQjhBAKJi1zYkJhVKiKS39yYSuVbT1uMkUV1SarvOvl2rqvTE2KlVN+Mz7GXWLZYM2EB93SHWyqnvnbGDXMXGPqqfvlwj5c+I4LD5qYhK5CS9cbNE1gQc/MDTfj8/bxU3rOLu3p63aBCxMzbwczn7hwodL+yQML6rdNvzvUx+ut6fl8WItfmwtGwDvHBXO4OSpb1981+ZfLsqy/ofuYC2mzideia/aX9QVMTRBFxnzrVR/q9o9EXw/Bh3ocn9dLgYL4HkoPzTePC4Qz80Jqquu5YDTHBX6oOWBqEzZ1f82b90Dtbry/ps23F1MUAAAAACQECzQAAAAASAgWaAAAAACQECzQAAAAACAhWKABAAAAQEKwQAMAAACAhEhszP7+sCLL1gvHsmw525FltwdLsmxe0dGfxdsFWTZc1tGlw8V4DOk8ZyKgBybKva3LXJRo/r7Oh80MdJzowMTRZnrxtkwLuo3lHX3d+WN9rklFR9WmpzqitLuhb8rSy31ZdnwpnivtouidXFfXGzT1tRXMPemu6XouKjzTjreldVkfr/mGHhsjsxVAWuxykDrdbXxX6J8zMePbZhzoHSVC75IoNJHS2WOTb22kxyZS/lAfs7in63Uqev6dPBWPLl56xbTf9T+T5Jw/NluZLOhxUDzU9aZ5M/+KLUSqm3o8pkd6jnVR9Fm3zcaS2e5DRHCHEEJnQ39eqPMV93VH7pzT7zG3hUCuHz8X81AyTEsutty8X14yUfpn9TxavaXnhvazZkuiB7r/dS/Gz9e8YebYge60uY7btsd9B5oJzPR3N34y43jFoyv6ftTv6jlqltPtn5TMNgfmG3fY0Mes3T/5/jQpt7+D4d7DauuBudj6KAT+ggYAAAAAicECDQAAAAASggUaAAAAACQECzQAAAAASAgWaAAAAACQECzQAAAAACAhEhuz/3z9oSy721+UZWcbh7LsWnlLluW2dWSoi15Omajq2i0RRV/U6+L+ioucNZHMy7repKozVLObOk40d6TbORc9x8VGu1jZ/rLpiub+Fw90nG6lo+NVd39MR/RWN+PHHFVN5O/EbCHQMc/UbS9Q1udbfHUgyzrndSx5dzV+zOoD82wWT7fNgYv7R9y8rPvscMnE8bokYTFHFbf0mHOxyy4W282Hrt5wQZ+vaLYXGD07jLdjqsdO+6y+7pl5Daz+id6a4/BKfGuOEEIoHJv5d6DvyUzEL8/MHJvrmnEsxn4IIZR2db2U6VuDBd2W5g0dWd45G7/R6vcQQijt6Qmlv6z7iNwehZj9REit6XfZ/Fi/p4fLJqbefLu4uPZUX/ej7nnzjt+P15sUdScr7Z6uP/eX9bVlzHxS3NdllR3dlnk6PsZr902djPkO7JmtZMw37tBs6VPa18c8eFY/8MpD86IT3FZSakuCEELI9uLnSplvKP6CBgAAAAAJwQINAAAAABKCBRoAAAAAJAQLNAAAAABICBZoAAAAAJAQiU1xvFLclmUb+SNZ9ts775dlaRPb5AKdXIpVeqxTZyYigGiwpM+W6+rjjSuuHbrMJZq5tLbUVLdFJfTNcvra+hWdTFS/rRO/Wpd0slf7vDnmXf3gSvv6wlUSZd4kpE1K5rmVTbJdXre/+kDfk8NndGpc7b7uDCq1aLCsk46KfX3dxxf0FFLZIsbxpOrf1Qmc7csmRcyklmWfbkV/n95ryjqjxuki7lIZXa920yQrmmtLtfX5cm+olDc99hu39Pjor+j+PCnpstKBSWo042euh3/orcYLi4f6HrdNgmtxzyXK6fs1WNTXnRPJZCGE0F8xcXlC9aGe80Z13Y76LZ0EeHS1FP19zj9PJ0Lxezqp0T2jiXuvmsTY9EjXyx/oATnPumPGf3epqi7pMN92qdAmvbKt37ndDfOON2nYmUG8bFLW92rQ1GUmBDFMiuY7qqbLejM9NzRf1/N94TCeAjxa0PNoZ0Ofq3FTz1/qG8v1A6YoAAAAAEgIFmgAAAAAkBAs0AAAAAAgIVigAQAAAEBCsEADAAAAgIRggQYAAAAACZHYmP0vHLxXlhXSOhL0Fy/+piz73Y4+5t07l2XZ0XP6fIV9HSfavRKP3Czd0bHxk4qOcp1sxCNBQwghd1fHgo6buv3zjImVPdeXZem78ehit19Btq8LRw3dFftrOoa08tDcr6LbXkDXm2Xj5xtV9PEab/ZkWeuKjhF2kbn7z+oofdf+7Y/q/pU7jv9ev6v7iNp2IIQQVl7UfaRzNt4n5/pw73q9Df1cMwMXQWy2t9iL78+RO2WUfmao2+GiqFsf0vNX6S09fw2u6gj1eTc+b+Q6uj/vvV+fa+GGjpuflPT4T4/1dfdW9dzmIvjV1imTgonnNvNJf1HXc9sLqPkwhBDyZguEzEjfk8JR/Pm4bQLSE9df9ZxX2o+fyx8P75TuMzqaPLuro+FnZd3XS/f1d03vso5dz2/pcVDY1+Og/Xx8bku/bPbRMO9B9+0yNeN/VDV/czHdfVzV7Tx8Jv4Mmm/o++i2R3DfE2peCMHPJ+6eHL5H96Hydvy6swN9rsZt/R4bNvW51Lzt7hV/QQMAAACAhGCBBgAAAAAJwQINAAAAABKCBRoAAAAAJAQLNAAAAABICBZoAAAAAJAQJ4rZvzFohJ+7/nPRso8t3pb1Xuusy7Irld3o759/6g9P0rR/Q1WW/KfNO7LsF9+nI1vnZR29PDvWa9zCw3jk5qhu4kIb+lzl13XsuovnLz3Qj3m4oq879UBE6YcQJpV4PRdFXTjWbTx4j455XXzdRMDndLyqi1GeFnQ7K1vx2F8Xodo7o59NxkVwr+jrrj3QfSHfMnHg9/Qx+0vxvjAzM0HORXcv63hrF+uLuHzL9Wddb7Csb3bxfvwZue0O5hl9vOGaboiba1Ids5XGJR21nXugo9eznfhFpGa6jfU7ej6ZlMz919XCtGi2HjD3OX988ljpYV2Pbxc37coq27odrcv6uaUneh51UdWHV+NzafWhucnGqHbydszT7PeRBKVb5h3yAbEvTAih+M26LBst6L6X29P9eXRWR8dnRrqdxVvxOWpaMN9XZqscFylfPNRzW+uS/kZZellviTNs6mur3423xZ2rvKuv+9B9670mi+x2H6OqLquZLYTy7XiZu7ZhTb+Pii1zLvEdlTJTHn9BAwAAAICEYIEGAAAAAAnBAg0AAAAAEoIFGgAAAAAkBAs0AAAAAEgIFmgAAAAAkBAnitlPpeahmInHkO6PK7Lef3X2t2XZ9fHqSZrwA1PcMtHkGyauvWjKFuJxqLVXdaTppKXjPTM6iTpk+ybmWV9aGJ4yCr2wFz9orqfjVXNdnSe6/LKuNy7rf0dw8arpqS5rn9fHHDTjMarN60NZx8X2p6b6Jqdm+nn3F11cvi4rHul7qeJ7Kw8Gss7he8onPh5OZ9ww80lel+WPTP8T3cGFjM/N2EmNTJS7aWPljomHL+my2fNtWTZ8GH/v9Fb1HDuu6Pa3rsiiULlvYup39Nw2WNL18l193Z0z8bLll/RY7a6bKPBDt12GbkfKbO9gt35Y0H2ytBdvyzR/um0OXGT+SDzvOf88nQiz9+vxnXlBR+lP9c42oXBg3v3v1R9S5Tf1+HFb0eRb8d+z/dPF7B9c098Fg2XdcZuvy6Kw8+P6PV7eMXPDYvx8xQNdJ8z1tS28Yb4RRex9CCFMzRZOua5uSnqkzzdYjD/U4qFu/9hsxeK+2bJ9td2HrMJf0AAAAAAgKVigAQAAAEBCsEADAAAAgIRggQYAAAAACcECDQAAAAASggUaAAAAACTEiWL21/Ot8LfPfSFa9sX287Le/7D9M7Lsk40boqRzkqY9Nhd1WXtT3yYX9dqrxeM40/GdCv70eAsm3rOq62UHOvpzsK6zkLNHOhbUtVNFd6so0RBCmGV1G11ZduBiUnX7Bw39UIv7up2dc/G25Ns6gjdlIqD7KyZuetdsLyD6TwghFEwM7LDu4q3jDZ2WdB93z9SNm3wnfm0uXvjdrrBvwu/np9tKI9uL/965oJ9Drmv63r4+2XBJD4T0VNeblHVbsq/UZFlYiPex7FCPq3xHn6twrDu0nb8yumzxNT2RuvlLxem3z8W3AQkhhMq2O5ce481X9fu2dVVvoTNs6Osum7ktI+b07rpuo9mRxEZ+qwh+9Q7DOyv1qh7f7vsqo3ebCJOSOV9Hj7mM3knHvus6T8U7U75tYtcHeh5q3Nbz6OJruuNOKiaKXrwHQggha7ZHKhzF29I+p8fq3MyHbl4YNvUxSzt6bhst6HqTon4G+Za4z+ZZd9f0RFS/q7+1u+vxdrh3N39BAwAAAICEYIEGAAAAAAnBAg0AAAAAEoIFGgAAAAAkBAs0AAAAAEgIFmgAAAAAkBAnitk/nFTC/3744WhZOqUjQz+98Kos+/zDj0d//w/qv3uSpj02F9k6WDUR0B0dJ1pbjUcX91cWZB0XJZxv6XONa7qN5Tv6Mbt6UxN9nRnG29Jf0mv+6qaOIJ0WdL1RSV93dVPH0bbP6OvOHemo1/U/jse5TvO6jUMT6V+7p687PdL3eFI0feu+zgMeV038rShyEeLBJb+LCOsQQpjmRVnKHPBdrvuU7s+VuyYeeqSPOVqI/17c08+he0VHGud2df+q3dJtHNVlkd96pKDHSHEnPu7GZX08s1tBKIpI6RBC6JgI+LTZZmOa1/WKh7riw0/GM8YXX9N19p/TW4GsvtCXZb1z+oblem6bDX0zS5v6pdrfiF+b266ksq8LR1U9/6o4fdcPkAzzjO57qZl+gKPG6bZyGSyZ93FDfzOUHsTnvYKKcQ8hpCdmXJmY+v6ynk8GS6fbmseN495yfGzlumZe3jVjtabHar6t73HraT23ue/302w9MqqYb68H+tqKu2afhhDfHiU91veRv6ABAAAAQEKwQAMAAACAhGCBBgAAAAAJwQINAAAAABKCBRoAAAAAJAQLNAAAAABIiBPF7OfS07BeaEXLFjNdWe/6YEOW/ezqKydpwg9M/lhHXfYu6OjP3LGOlU59pRn/fdHEapq4aRcdW94yEfwVfUwX6z9b0dnd8+N4THK+bWJqS/rfA2b6NobSno41bV00UfomBjbj4u1FO130bXagj9e6rNtYv20iW01E76Ro4sxN5HRpPx75f3zedAQTR102cboyuvd0CcjvCtm2fnadZ02E71D3h9L9eP8bLOsHUbpzuv4wjScJhxBCyPZ02WysD6pi0kMIof++eHR8ZTM+P4UQQq6rDzgu6/tf2TF93cxfbruPaV5XrN+Kt3Nith2pbOprO7qq70llS28FMqrpNubNvTy+XJJlaguOXN9MDi4W35SNavHCOf88nQiDdd33ynf12BktmO2PBrpDzM7p75rMnu6zhQM9DtT2NempbqN7T7sxXr+rPxJnOT1vu7HltpyoPYw/n96Kvh9uCwE3n7sI/sqWnn/7S+Y73DyD7CDemHlKt2NcNlsSvE9vV5IRr2/3XckUBQAAAAAJwQINAAAAABKCBRoAAAAAJAQLNAAAAABICBZoAAAAAJAQLNAAAAAAICFOFLNfSo3Cc8UH0bJ/3XpO1hubDOL7g3gUfWjeOUnTHtvERGdmOiZys6YjPMfV+O82ZnRNR87WX9URqv013Y58y0S2vqXb0hnrzOyUSDx1cdOjor6P1Qc6+rZ9Pi/Lcia625kU9D0pDOP3clh3sbj6XNmefjbFPX3dw0X9vAcmVjbf0R2svxwf8uV9HWHrosfHFV2mYqyJt9bcvUl19HTt4vmz8SR6G1M9sVtzmNjinh5Xw/O6r5fe0nNNf133zcKN+MCbZXUbZ1kz9s3WFoNFPeZqd/QWCDMXo2ziqFV0dOnAxE03Tdy0ee8cX9B9y22dMi3ofjeqm+sW25yUTIS4iwKXW3qEEHJi/k3rVy3eQdW3dN8bmjnKzZXTgtnK6IHebsJF/udaemypOdH1y1nOdGhTtPc+/T1Uv6vnBjf+Bw19M1PiRs9OtHr4f+WPdUPcHKu25gjBz6NHT+v7pea28raeh3Y/oN9V5R19bUUxb7u286kEAAAAAAnBAg0AAAAAEoIFGgAAAAAkBAs0AAAAAEgIFmgAAAAAkBAnymG5t78c/s4//evRssEFndCV29VpdCr95unU87JOaqrTXFxi4aSs01Le+Nv/UJYlxl/8YTfg7X3sv/icLBtX9LOZm6SzyraJ2zJpR61Lut+5ilOR8OhSxFRSWAgh5Lq6bOfDOv5x+Xs6Gc6lQ/VMwmPtfjydaFzVddrn9b/jLL+k0456K2J6MffxUby8vxKu/eN4P5uZFNHivul/4hJnOgDKphlOdVBYyB/pdtz4T5IxD/29XZ3K+1RhT5b9+p1PyrIPLN2XZUsf6sqy33ug2/ILn/529Pd//t/8BVlnuOBSHGVRKB6ahMcV3VEyIhX27WT7J0+Gq2zpuTJ/pN/Rex8oy7KMHuIhNdSpZdO8njfqt+NtmZT1AN5/Tn+urL6g58rO2fizcYnDj+Klo5Vw+f/4bLQsv2/SLRf1Pcv04vUq98x8v2GSpJd1f8genC6GL9d1yYQuPTX++xt/Kxlz3g/Cn/vS35BlqZm+VyolO4QQ5indcd03Smasz1cyc9tEJLXW7+i+NSnp/t9b1WW5jiwKg6apZ76/SvtmvIm5eWrav/SqnhAzA/OOWI5/j7r5nL+gAQAAAEBCsEADAAAAgIRggQYAAAAACcECDQAAAAASggUaAAAAACQECzQAAAAASIgT5azOMyGMq/FYysoNHTM8+/CxLEt/qx79PWWSiYdNXTg3V+SiS/H9Ma7qstpDHcvqYt7TUxNHaxL4MwNdr35XR6Wq2NNZVsehpkwbZ3ldr/pAR8COGrozT0qny6qfi2vIdXU7Ft7Uxxss6uemouZVpP2jmmfnYbQSf/A5Ex3deY+OGS/ejc9fKX1bbJS+q5fRieCJ8cXN98iyjYqezxeK/VOdb2b2Xvjwyj1Z9mZvLfp7ZqTHY/O6HvtTESkdQgjToi7rmujo8p7uDC4WOy3eVxkTbd9b1f0/19OTZU7vchByXf3i7Kzr86XN3DwR833fzCfV+/q6D54ryLJ863TbHLyd1ETH6afH5l0x0mXZTrxsuKTbYb95hrpfzk0kftZE6Y8WzOQ20/VK24+5v8oT6PCq3uqnYPplvqPvcdpte2EejdtWon1Gd6JRI/57adf0Y/Pt1XxDv4fd/Ns1c42bR8dlfczO2XjF4oG+kW4dMcvpNpb24hOi226Bv6ABAAAAQEKwQAMAAACAhGCBBgAAAAAJwQINAAAAABKCBRoAAAAAJAQLNAAAAABIiBPF7KemIeTa8VjK0YKJsr1Rk0U5sUTMDPThyg91pubxVZ2BmTaRs/j+WHxNZ4j3V3TkrIqUfjsu3r6/ovtJoaUzZycl0U9MF6/dN9dd1dfdX9RtzJg45lxPN6bQ1mVHl+NtcbGyub4+XuecHlNLL8djZV0E+iOZpkKmFZ+6ZhdMzHvb9D8RXeyiwoOJhh+s6A49XEj+PPQL51+QZV/YfU6WvbexKcsORhVZ1p3omPS02XOlO4lvjzCu6Huc7euyYUPPC4NFXS9/bLZ+MRHQbiykR/Ex2Xpa9+PaPd1hD6/qfSGqD3S9WUFfQL6j21/a03ngaguRhRv6xd9b11v5uPlrUoy33z2XR+G+hya6q4dZRc8Ns178eNPLel5L3S/JssK+2b5Gp51b47p+5tX7+qa67ZF+VC3c1GPAbVHjlPd0/+mu62O6MeK2AlHfIS5Kf57S/aC/7L4DzXxyoK9798f0UqZx00Tmi20h3PdVpq+P11/R7VDbLdmtVnQRAAAAAOCdxAINAAAAABKCBRoAAAAAJAQLNAAAAABICBZoAAAAAJAQLNAAAAAAICFOFLM/z8/D4Fw8NrR6Q0dnds/pWMrCQTxj0tbZ0+vKzECXzTPvvpjXd5qLQnaR0lMRQfqn9fT5ZjqdO9TuuehY3Zb8cbyPd87qPt4+pxtSaOt4WBdT7WSGp7uXKva3v+ymAn2u2l19j+fikN+PeOv8kYj+3SrLekOzFUj/bPwZpXt6PpmV9bVn224ekkWJsTepyrLJTF/bK60NWbbd0dutnK21ZNm1+rYsOxzFI8bd+HZbc5hE/1AzUfRpMx5dnLY7X07MG5Uts4VDXZ8rb7bfGCydrlOq7SlCCOH4gn4XFNrxsTNY0XVGNd3vKlv62agtF067tcufmWdDGC7F7+nM3M7iQ/0eGS3E70v2po7SHy+aeajrtpPR/cEds/TAHFPv5CC3JPhRZrftaJptO9pmfxczZ5RMXP6wrs9XPNKDobQvoui7ZksrsUVICCF0N/QYT5n5ZFLU7a+/pW/K8SW3Xoj/nh6brYXO6G+lyvbJt9eZp/W44C9oAAAAAJAQLNAAAAAAICFYoAEAAABAQrBAAwAAAICEYIEGAAAAAAnBAg0AAAAAEuJEMfshFULIxuMnexsm0ryvYyT7a/F6lXt67TjWCdBhuqoz2Qu3TCY7vi9c/PukqPvBYMn8W4HuWqFw7KL0dT0X2api6otH+ly5jo5XPbxi4qZN+ycFfb9SJiK6cKwLVaRr7e5Qt6Oko4Izfd2O3mp8eplnHi9uOT0NoSBS2VtXTex/QZdVbsfbOqrr/pw5MnPUWd350vs6ZjsprnfWZNnfvfg7suzL7ffKsv2Knrj3RhVZ9rXty7Lsb1z8WvT3f9K+JOscX9T3P23mjMKR7gutiy56Wfc7tyXGYDneTjcvlA50PPcsa+Kcj0+3TUDhSJ+vcCSLpHFVn2tmvlbcnJIX86HbbuGRzELIduLnVXH5IYQwXDbvETGn5MR5QgghpPU8lDJp7SkT750x24SMq26LGn3M4eK7b5ujohmP87Tu0H0z5pzynn73T80WCDsf0m3Z+Ho8+36W08+6dVF/ay/c1JOs25qj+lBn8A8Wdfubr5sPSKFvth1x9/jwmm7H4qvx9qcnelzwFzQAAAAASAgWaAAAAACQECzQAAAAACAhWKABAAAAQEKwQAMAAACAhGCBBgAAAAAJcaKY/dQwFUo34zGY/bMm3resj5k7iMdZdi7qaMzSlomVPTQR1u++lNd3XK6rI0izZrsFFwE9quuyXE8/1LlJJi6aeOhRPT4s3DYB2Z4uq9/T5xrWdV+u3dexsu1zJip8amLxh/FxNSnrOt11PU24+58Zi7LHHIfzdAijWrws2zVx2w/1dYya8UaN67o/L7yqj9equih33ca/fPPPy7JmXu9p0BrrDOWtbj36eyGr++XvP/tbssz5VPGVU9X7fvvluu7P9dt6XM3yejzufEg/08Zb+n3VXdPHrN/RzyDXi5d1zugI62FDX3febOkxNduOZAd6wHbX3T0ZyLLO+fg1uDm2flffq+2P6HaUduLtn77weNt9hHQIExE5Pyvqe1a5bbYSECnjc5O67uY89w4cruj7Wb5nIuDXdT9Km1j/tB52P7IGS/o+5jv6Po4r5hvXbLHjtu0o75ptZrZ0f1Vx+tsfPl20/bCp65UO9cW5b5S0uSchpe/J4bX4fa7c1/fDRfAvvKkHwOE1sW3Kt832J7IEAAAAAPCOYoEGAAAAAAnBAg0AAAAAEoIFGgAAAAAkBAs0AAAAAEgIFmgAAAAAkBAnitkPqRCmpXj8ZGFXH2q45rJX4z/PczrmMj3Sh0vpdM8wWjSF+L7oL+t+UNzX/WBc1cdcfG14qra0z4vM4hBCeqKjUke1eKdMTXWfPBARqiGEUDgykctbOnt4XDVtbOho1sqOzpwdVU/+bzLlXZNhe4rI/LS5j49ing5hKlLli3v6vrSfMZHmYrsPN6+pqP8QQsgMTb61KWrkdDR5JavHwfur92XZ9cJ69PeZa8gTLt/Wc33KdL9JSY+PpVfNVjIZfS9rD3Rb5nZ7kfickhnpCxg0dfvzx7LI6m7oeahkoru7ZjsAFRVeOtBzzdxMXc3Xdb3B4jv/79CFXbNtyXl9z4q78bbOzO5BLr4+19Fl2Zt6bhubuS1/pO9ncV/X61x4931/Zft6rLpIfPcd6+S6Jt7ebcExMONYbLOz8Uf6Q9xtV3J4Vfe7tRf0VjL7z+qtZGoP9Ph393nx1Xi9/rJuf/WhPleubbau2I4f025NoYsAAAAAAO8kFmgAAAAAkBAs0AAAAAAgIVigAQAAAEBCsEADAAAAgIRggQYAAAAACXGimP3ULIRcOx5ZOa7qONHqm/o0M5GEPq7r4w2XdNncJEeXNlmP/qC5mNfems4KdlGoez+m45qXXjF7Lhgucra0Z3JP1fFMvHX1gW7juKbHRr6l21HZMh3dxIjn+vHnc3xet2Oub1XId/Tz7jfjFWcmWvxRzDMhjBvx885MP8rv6wvJHcfrqTj/EEIYLuobnW+Z+PSmvmeXynuyrJHREcTb47osG4t88mvlbVnnSZea6WfT2dDzUK5n5q8VN2foeuOyib7vmC0sxCHHC7od2YG+brV9SAh+jC++qrd3mOV0Px/X9EF7K/G2lHf18QZmznbXnW/Hy1TU/0moY7gxXrutn0N/Nd5W912THpu5pqHruesf106+FUAIIRxf0wfNdN99319Z8b4NIYTMQN+r4aKeo3omAn6e0WPExbn3VvT7v7Qfb+dgyWzTUDFbCJh+11/RWyMVWm7rLX2fZzlzT8SWJeUdfbz2OX288o6+7rS4brfty7tvxAAAAABAQrFAAwAAAICEYIEGAAAAAAnBAg0AAAAAEoIFGgAAAAAkxIlSHNOjEGp34+kmmbGOIknpQJTQPhtPRDn3JZf05hKzTIpbW8fYfPTvfk6Wta7olpQf6vPV78XPl+3pGBuXfDVPmbQmk9BVva9TBKcFXW9kUrMmhXhb/viXfknWeVI88/l4X8j09f3PtfXxtj+iUygzJoRyuKifjUpTDSGE4YJJEhrHf58V9PhtX5JFYWbSW4tibEz17XgkuU4I638UL0vN9NyQa+tx1z4XT82qbp080TOEELpremptXtf37Etf+klZNhNjLgSdmBlCCBORqPW9I32vfmdTdJQQQn9ZJ4yNqrqNZZEGFkII45JJhruky5rX48f82q/+sqzzpLj26/qdpLhUxdK2SVzUIaChfUEPWBEQGkIIYVrU/XxSiZe1L5rEyyPdfjc2hiJRcfKHssojyfZCWPmT+LFdKvE0r+9L83p8vuls6PkkO9Tj2H0zZIe6He47Kj3R5ytt6+dXOoyP1Y9/9z+SdTKmjcOGPpdLCiwc6/YfPGsS+rZ0WyriPfGVX/8V3ZAnxH+583z09ze7K7LOtapOCP7m/kVZ9gtnviXLXuqek2UPBguy7DtfvSbL0kOV3qyftSsbLOn+o+pNviqr8Bc0AAAAAEgKFmgAAAAAkBAs0AAAAAAgIVigAQAAAEBCsEADAAAAgIRggQYAAAAACXGimP15NoT+cnxNV7uvc01Tcx1LWdqPR54O63rtWDER0N01HQE9bOjLbdw+Xay0iydWEcQ908bsQN+rWdZtPaANF/V1p832CPljE4ttYsSfdDORlDptmP0iTN70cE3fx+JDHcs6WtER7+mh7kOFQ91PMoP475OyrGKjuycjXTZP6771OObpEMal+HlLByZyOm0iyA/j9VqX9H02u32EXE9f+zyj2zFY1P3h8FmzzUPHxavH2+LmmuGCHt+Zka6XMf0h23PjRxfZ/mza8qRT23pMzNYWLl58XNNl2Y4um5upftTQbcmabUnyrXhZ/lCfq33VfGNM9Lkq9+Jzc9pscfK48m3d11sX9Q0t7cevo/pQvwumRf3umZnpa+qmNj0Nhdpd/a3k2pLtx59fZ0Nv41De0f2rvKv7g3tnubLmG/qYuY55t5j79aR7vb0W/X021/dxf1SVZa1hUZb9xoOPyrJUSvcFd8yU7q5huB4vLN/Wg2Na0sdzc6XaGsltQ8Zf0AAAAAAgIVigAQAAAEBCsEADAAAAgIRggQYAAAAACcECDQAAAAASggUaAAAAACTEybLS5yGkJ/Goy2lBR26mTcyl4qInU1MT8yxixEMIYZ41kcymaOGmjl5V9yOEELK9eD0X89pv6jVzvmMi8U2sb3FP5wn31vOyzMWB/yjHyqq+lzZRzq5v5ff0zeqf1fHJqb6uN66ZiPQV3RcK+/FjDpd0H5/nTNTxXT2F9M7Fr82Ow0eQmoWQE5Ht+8/p9hT39XmLIma/sqXvS2aojzeq6nHs5q95Xtdrvq7rFY50O7vr8Wc+NfNQ3kRK91Z0v3TzkNsm5PCqfm75tunrjR/df2OcFuPXnRno+zg2EfzZrj5X94KJFz/W97h6z7zLVvT5VGJ294J5j22ZLUlM3H/ncnwemhUebx6aFkI4ejrepjNf68l6RdNnVX/OdXVbu2v6eGmz7YKb22ZZfa+HTT1WBwtmPKbi9QqtU25JYr+VdD+qPNAfia3LOkPdzfduu6In3f32QvT3v3TuJVmnZbLo37/0UJYt5Pqy7PpxPO4/hBCKZb1PyG5xVZZV3orH6ffO6rFRfqDHhtvWwuzEJP3ovt0AAAAA4AnDAg0AAAAAEoIFGgAAAAAkBAs0AAAAAEgIFmgAAAAAkBAs0AAAAAAgIU4Us5+ehFDai8eXusjTsYmcLu3F40lHdR2vmjJZlrm+i77Wx2yf07ei9kBHobfP6nqVnfjvnTM6plPd3xB8tH1mZGKxTZR+ruvOp++Xez5PumxPXJu55KxOVQ6juu6TuZZ+qMU9fcLMUJ+ve06Pt0kp3hYVvx+Cjvt+u7LqW/GxkR4+Zt9J6b658KaJm1/T13h8MV628l19o6dFE29ttt+Yp92WJKav9PVYnRZdZH78mDMzn0xK+tqyAzPH1nS93qo+4fLLei8Wt4VL1tyTJ91UjNVMS9+PWdnM5219/+d5/Uxn+vUROud1PRWlH0IIBTG3TfNmbOjXsC0LfdEnZ483D2W787D27fj8cHRFx4y793jtbnxLnOGi/uZx7+nabT2u3HfZuKKPWTiWRaEktisJIYS5OGTabDsyLus2Vrb0tXXXzDdiST8AN8e6rVN6q/Gyv3bnU7JOweyBcP1IR8NfqB3qstKBLPvu0TlZdrm6L8u+8YHflGXvKJ2yb1372udk2aQc/919D7lvLyejdrsy8yR/QQMAAACAhGCBBgAAAAAJwQINAAAAABKCBRoAAAAAJAQLNAAAAABICBZoAAAAAJAQJ4rZT83mMmK5dKDjSYtHJ4+VnhZchLUssjHVlS3dxkFTx2qmTGR27Z5uzERExDbe0nVScxNbrFNZwyyr73Gup697uODiRE3kd9dkgz7h5qLrudjoqU71tbHKKZMSPlg2Ue3mmMEcM9eO9xMZARt8LP5gRZ8sNY3XU/f3Uc1TIUxzJz92eddE8G/Eb+iorqfIwYI5mUnwdltUNG7puSHf0mXjmm5nV8TbVzf1/XDx9aO67nyVLd3GUc3EvJtbmRnotmQ7Omr7SafGz6Ss54X8volOb+j7WLth+vmSPp9/F5tBIJrp6qTNHOX6T0HcExvN/wimxVQ4vBrfg2Dxus7iHi7ol8WoEX8OhQPdz2dZczwzVu0WFW7+MnHzWfPNMBLR/aUDPQ+5b55hw3yzme+o9nm9b4TrE66Prbw4iP7+8//xd2Wde+NFWXaptCfL/vjwoiz7D8/+kSz7e4Ofk2UvH27IsnBWFz0J1PYOIYQwE9sEuSh9t62Fiu0PwX/rKfwFDQAAAAASggUaAAAAACQECzQAAAAASAgWaAAAAACQECzQAAAAACAhWKABAAAAQEKcMGY/hFw7nkM6quro0kJLZ6UePBuPiG2+YeLrizrmsrSns3h3PliSZQtvughrHefq4u1notq0oOvkj00WZ/p00fa9Zf2Yi+bZFMy9nKcLp2rLE83c/uK+LjxeNLHYR/rfSKZ5E7NvomOzPV2o4vnHFX08Fw9bODCx3rV4+x83Zj81CyEnIqK7a3qs1u/oMV44VHG7+uJLok4IIUzMNiHZvq7noqML+3o8Hp/XYzwttgmZia0KQgjh4KKO7i7t6HsyNduVFMzc1l3T7XfR/ceX9Zz+oyp3bLa9WNX3OHes+2Rvw2yTs6fr2Qhrs/VIfyXeJ12s+rimT1bcNw1Rw+0xd4pJj+dyq4r2Of1+nJguq+bGaUHfzNbT+vmkTGx87Z4uKxyZrS3EVkshhJAZmTIRXT5P62c3Ms88MzL94VB3pEnJvXNlUSibeahzJl7xH9z6tKxTL8Sj+UMI4UptV5YtFbqy7H8++AlZ9nOLL8my3Xpdlj3p3PfQuBrvr25roZG5VXPzjT5sijpmFcZf0AAAAAAgIVigAQAAAEBCsEADAAAAgIRggQYAAAAACcECDQAAAAASggUaAAAAACTEiWL2Z5kQRg1RxaTcts/p05S34nGuo5peOx5eM3HB2zrDtvpQR8fm22NZFtq6qHWxKMtKB/Go1/6SjqJOzXRM57iq601NpHHxSEfOugj1SdUcNGUe+BNO3ZN5Rj+b1jVdVtzR/TVtYpBTE3OPTdGortuiqk3KuiOUH+p+11/X9WQE92P+s9A8FcJURMTnOvra3ZwyE1PULK9vdHqozzWt63pTs0PFwpsiizqEsPORsiwrHOi2qFjs7rq+H423dMd0Wxnk22b+qujzlXdNvrrpL9UHeuuBJ11avJIGKyYSf9fcLBMrnx6ZseFeA2b+mpltQvJHYvz29PFUNH8IIQzO6oZk2vH+6q7rUcyzKfkuH1f0+F9+ScerDxfjjVLzXQghnP3DvizrrZm4f7Ndkdt2pdDS/W9ScrH48efXXzTzwp6eF0ZmPkmPzIeNue6svpVh/706g19ts5NJ63b81PIbsuyBymQPIXykfkuW5cw+Fd/qXJJl3p1T1kuGSVnPG7OiKOua7bMKZnudmn7eKRHd77Yq4S9oAAAAAJAQLNAAAAAAICFYoAEAAABAQrBAAwAAAICEYIEGAAAAAAnBAg0AAAAAEuJEMfupWQiZQTxG0kXHD5Z0juSkHC8r7eq4ynN/oGNqjy/qWFkXYX/0tI7Lz5vo7nxXt3Ncjq9/8x1d5/AZnf2bP9btGDb1PS69oqNX22dNFzDxny5O+0k3VlGpJg91biKlHRezn9GJ68Gk8IZcR7dTbSEwy5qHbS4tNdb1ciKy2CQBP5L0ZB6Kh/GDjKv635yyPT3u1HwzLel5bVYwEfwTM2eY+WSwpMd/7d7ptstQ150d6krDur7u0oGZ80puSwmz/cNcl/VW3GsqXnbln31O1nDjY1Iy22Uc6Hpua4usiWweNXS9f/8v/UH09y9uPivr/GeXvizLfuX+T8qyDzXvybLxXPeF58u63t9/7TO63upm9PfeREeZ/8zKq7Lsd7afl2XjWbz9e+XH3KJhpqPjXcz+4TX9jVI8PPl7pHNGHy8z1serPNRbC02L+pnPzXBU9yOEEAoH8fOp76S3k+vrc3U39Dyq3h0h+C2h1r5ltjPYiD+D+y+ckXV+dazLMiKSPYQQfrdq5lHzPeHmvd4ZfU9++/hj0d/T5t3vvhncVgZqu5sQQphU9EFdW37+578uy15pbUR//3Mr12Wdy/ldWfbaQD/TZrYb/f3v/2pL1uEvaAAAAACQECzQAAAAACAhWKABAAAAQEKwQAMAAACAhGCBBgAAAAAJcaIUx2k+Fdrn4+k4ua5JC4oHNoUQQijvxmNn+ku6aaOGLivv6Bib46dMPZMa6dLHdj6oj1kSYS+VTZ2Ys/yiju6bm6S9yqYu2/5xnWhUeaivrXZPp1y5JKonXaZv0omE8qZOveqvmYS6hu4LhR2TpGX+aWVa1WUqQTHfMgc0t6O4Z5Khno335dMmXv4/zZnNQ64dTwTrnCnJetOcbmt5++TptG6uGW7oc7mks3napGJ2T5fGmuvFj+nSZH06rW5j4fh0bRzVzNy2rY85WIg/n4JJXJzp6TDMs/rZjM24Gq3rRLz5nj5h1sw1Hy7fihfEg8dCCCH8w7s/Lcv+5oWvnKpeOqXvyfawJsvOLxzJssV8L/r7wbAs69wdLsmyhbyOhvtg42709+s5nQb9KDLjeag+iL8jc139zN04VkHB45p5vzT03F3d1OO4t24So83c4FJ4U1PzrluPJ3Rmh7rOLKPHxzTv5qHTtb9+V8/po4a+X73l+DNImTjD9ES3PxsfHiEEPw9Ny6c7X/5I9yEV4jox58oM9LnGNTPHXtZjsvSyfre79MdP13X662r+OPq7S2r84tF7Zdn54qEs+2brUvT37vRFWYe/oAEAAABAQrBAAwAAAICEYIEGAAAAAAnBAg0AAAAAEoIFGgAAAAAkBAs0AAAAAEiIE8Xsp2ch5DvxiMzSjo4Z7m64XOOTtOBPqWjlEHy8atXE26dN9PVgUZ9v6RV9vmk8VTaU9nR8/WBJVAoh9Fb0errQ0u1ffllHx85MdH9/VT+36kOTVfuEm4lHUL2n79VAJ0CHjN45IZSv6yE41onT9p9W3PlG9Xg/GS3ofpxzEfxG8WZ8K4bU8OTbGPybJqV02Ht/PHK3+YYeW+1zemwNxLYehZa+L+1z+tk1bukxNy2YMbeo73WuI4vs+O+txo+Z1k20cdnFfT32p0Xd/mxfH7O6ebo5atSIl6X0YwtZk66e7Z9uu4lRT78jFp7fk2Wt7y7Lsv/x3qejvw+nut/Vcnrw/+OHn5RlubR+pvu9iiw7VzmSZZ9Z1fHW/+d2PKr6qeqBrPMv33i/LHv+zENZNhT7KsxVpv0jmhRT4fCZ+BxXM3Htrj+rfts5o/tXrm1i6s3WImn9yWbntsxAn690oAfeoBkfW7X75l6ZKP2jp/Q9qd7X7SgOdNnIbGeQNnPiwlvxm9m94LaLMls7jfQ8VNox7491WRQmJd3+3LE+ZvUT8ch5N3dlzBxrt5K5qaP0XX/tndHz1z968NOyrFmI72fwrclFWecTzZuy7DvHF2RZKRO/ALeNCX9BAwAAAICEYIEGAAAAAAnBAg0AAAAAEoIFGgAAAAAkBAs0AAAAAEgIFmgAAAAAkBAnitkP8xDS43gkpItzLR7pCMzcsYieHOkI0nlGn2tsYlKPz5p44vuni40vHOjsz/5KPN63uxGP5g0hhM4ZvWYu7+h7cnRV11t8Vcd42jhqE4s9rjxeRHGiiUvrbej7kW+Z+zEz/dVF6RsuIXquh0DIizjdYVNfW0Yn14eZOddc7dLwmF0nM5qHuoixHjb1GM8O9TX2l+Ljp3Cs66z9XzoSfP9Di7JsVDMxwz19vlFd3+ypmX8r2/FBPjf/PGdS18OkYrb7OHDZ/frZdM7osoWbOrO5shm/7v33mY5pZPu6TG2bEkII9Rv6fPszvQdH5umuLDtbbkV//5nmy7LOreGqLHuhpSOg9wc6Sr+U0++4N45WZNlyXu8LsVE+jv7+V5a+Kess5PTDeaa4JcveGMSzx6eP+e/TqVkI2XhKtx0j6hsqhBCGYowv3NDPwG1tMarqsoxpR95s25Ga6TI3p6jtX/pii5MQfLT98osmd90xWzup90AIIVS3zPYiausUc678kT5X38TGlx+Ybab2TQT/hv7Ym1RlUeh8PT7G0+ZZ99fNtj0d3cbSlttSRZ8v29GN+asb35BlX2k9G/39bPFI1tkcLciyv776NVn2e8fxbUIyxOwDAAAAQPKxQAMAAACAhGCBBgAAAAAJwQINAAAAABKCBRoAAAAAJAQLNAAAAABIiJPF7Acdcd9bU5naIZT2dfRy50Ix+ruLca9s6gjSlImHbtzR7Rg0dHSpi47trevrHjTjFRdu6njYsnkiaZMqW71nYtJNvHhvRV936UDfzHnqR3dtn2vH+96kbOJQzbOZ6l0VZPRwCCH0PqBjpdMP4uMmhBDmGd3OXDv+3EzSa+id0/0gK44XQgjTQvygbjw9imkuFTrr8YFS2dFtnZntOcp78VjgaV7XOfhgU59LTwt2PObbOp44O9Bl3VU9cZRvxWPqR02dGz8pne4h9Vf1hRfMditlE/l9cE339eKR6biCikYPwcdDL7ym+8Lhe82YM7HSxW/qfOsPvP9u9PffP3yfrNOf6vt/vnQoy35q8Q1Z9uW998iy9VJblg1nuk82RGT+b+x/TNb5zs55WbZ8Ubfjxyu3or//i7SZfB/BPBXCpBQvK4n5JIS32S5DDMn2eX0v3dYcbmsRp3Zf35ujy/qFVjjW113uxr+/0mMT/17R92pm5uZ5Wpd1zuhjum/LYd1E8N+P70VT2DfbvpgtPYZ694owauhn6r5RCge6/W67pbm4hElRn6txQ5/r+LIZGw393KYlfT5X9s3207JsZxCff3+68Zqs879ufVSWDWcfkmVXS9vR3/MpvS750f3KBgAAAIAnDAs0AAAAAEgIFmgAAAAAkBAs0AAAAAAgIVigAQAAAEBCsEADAAAAgIQ4Ucx+ajIPxb14nviooQ/VXdNlla14xGRpz8TDzkxs/KqOjp7rdNVQv6Nz0vffq4/ZfF1HZM6y8cjQ1mUdhexiXnNdHU86T+uLa5/TZbmuLLJRtZnR6eJ7nwTjavzapiUXD6vvsYujzQ70PS68LjKcQwjDJd2W9FAfc7QQrzfP6TbW39DX1juj67no/seRGc9DVcwbxyaOepbT96V5Iz7+5zW3jYDZCkS0L4QQ2ud0G/vL+nzpkS7L9fXNPrwW70cLb+ko7WFDz3mlXXNtF/TcNimY98C2PmbtgS4rbsezqnc/WJN1umfM8e7p9h9fNluZ9M2Ya+o5fWxipf9g/1r098lMj8eHnbosG9RN5HdazycrxY4se3H/jCz7+Go83j6EEK6316K/X6zsyzp/5eK39fF68eOFEMLvbsa3Jdgd3Zd1HklKbxni3gfBzIuV7Xhf6a7p47mthVRs/9u1Y5rXc83CTT1vDJb0+JmV48dMT/W5suabp7+i+/NE78wRFl/T7e+abZPcfT54Nr71QO+Mbn9pR1+32yqntK3njP66LAoTE0XvtglKj/X59LnM8Sb6eC7uf1zThYV9PT7OFo5k2eagEf39l+78tKzzfPOhLLvfW5BlhXT8vTNS+xgE/oIGAAAAAInBAg0AAAAAEoIFGgAAAAAkBAs0AAAAAEgIFmgAAAAAkBAs0AAAAAAgIU4Usz/Lp0L7fDy3tXSgM0g7501c+zAej+niNke1k8d+hhBC8VAftLeqb8XyizqWdVw9efxt83UTb72oY16zHZ2FWjBL7exAFw4WdJmLJXeR/0+6wmH8uofm3zPS+pGGnImVVbH3IYRQONDnK22f7t9WVDT8qKHbMdNdMgTTDdKj+LkeN35/lkvJrTvKe2b7gbE+sdqCY57Rzy5rou2nRf183Dw0KZ3uuaq47xBCqOzEJ6LWxXg0dAgh5Lv62oYLeq4stEzc/1XdyHzHRG339f16+CkVK6/bUb1htgIw8dyTRR3Pf/6pPVn28CUdAT9d0nN6LRefVNpj/dw+e/mrsuz1/oYs+9L9Z2TZh9fuyTIXpX+hcCDLHvQXor9/feuSrHOwWJFlH27clmXf2Tkf/X02O913xJ9JzULI9UShmePSE104qotxYJqaNvHvwZ2rasbcQB902NTjJ9/W9VqX4vUmIn4/hBAat3T7K1t67LhtAgaLLp5f3+iZ2bKgflds+3JVH2+wZOao27qNUz38Q+HgdH26e/bk7838ob7HzT+/Kcs2v6f3Ahg1zDu6qNuYMluPpM1i4uMLb0V//1b6KVlHxeWHEMJOT2/v8rPLL0d/L6Z0P+YvaAAAAACQECzQAAAAACAhWKABAAAAQEKwQAMAAACAhGCBBgAAAAAJwQINAAAAABLiRDH7YR5CehqPwZzmdbxn87qJxxTHyx/ruNZ5+nQxqS66//iSXqsO6zrXNDvQsaCqLdO8zmsd103Ma1ZnQI9Lup47Zt7EYg+aul5xHi975vOfk3UyPXNtORNZXtZlhT393Fycu9oCIYQQ/vvP/lr09+OZvv+3hquybGeko1dvdZdk2WeWX5Vl/8u9j8qyv3P5C7LsXx8/F/19Z6Db+OHGHVnmImxbk3L091/5fFvWeRTzEILofvL3EIKNqu6KbTZcfL2L9Hfx/I47X65nIuyfMZHHb8R/L7ZM3L+ZR10b8219zOWXzf1K6/MNF3SEcuEwfk9m5n3ktsSY1c1cs6njxffu6Qj7rJnbZnn9LrtY2o/+fn2qY/v/xdaHZFlVxPaHEMKHVh/Isu5Ev/8+1RCdK4Twe/vPy7KRiMV+qqGj+W8fL8oyNw89u7QVP15Wx1s/inkqhKnoEhPzPi4c6WOqLYTKO/r6+kt6QFa2dL2h2qEihNBd18/8yGyXUd7Ufb0oxmrPbOfTb+qxXxnqc6VmZhybeS9jtmLJH+l49dal+P0qbZr3gCnqr+h2ZMz8NS2ae2K2+8n0dVlZXMPQbBOw+Sc6Sj/XPeX2FnM9V7rvuW8cXZZlWTFvDKb6XB9YvCvLtof6O+r+KD5/jdQeP4G/oAEAAABAYrBAAwAAAICEYIEGAAAAAAnBAg0AAAAAEoIFGgAAAAAkBAs0AAAAAEiIE8Xsp+YhZES0aXddR0Xm2zqOc7gQXyO6iHQXN13Z1NG5g2Udk7z+zZEsm5RNlLuJ91RbCBw/pW/70ks6Q3Va1Pe4tKMb0jmnY/37KzrytHr/5LHY6aE+3mhRHy9roldL2/r+9zb0McsPdb3OFd1PLufiUc+/dvBJWWctdyzLDsYVWZbP6Of2taOrsuw/v/z7suzOaEWWTUS8dSGjI4S3xzqP+UZbn+vji29Ff0+7wf0IMuN5qG7F2+vGo0mzDb31eP9r3NL9K20imcNIl+XEvBCC317EWbih5xR1zFHdxNe3dDuGpp67/+mxvpeDRROhbLZHUe+CowVdJ1swc9SKHgfVm/r9MVwwfcH8M2hpU9/LseiwtdxA1vnM0iuy7OGoKcu+vHNNlm2U9dz2G5t6u4+LFR2Zny3EO4rb7uPjq7dk2c3Osixrj+Pbo4zFXPio0tMQikfxzuli9t08VLsX73+dDT0+3Du8sq3Ple+aSHbTnWt3dGF6YrbLEZ8htft6zPVW9M06uqLHY2nfbe0ki0LGzNu7P6a32WnciV9Df00/t+5F/Q2S7pu5eVNPKKOG2V5Ef+KGtH4EobcRvycZc7zigW7HVO/gECal030bjBq63pXKriwbzuLPx8Xl3xnquaaU0c/0qcJe9Pd8St98/oIGAAAAAAnBAg0AAAAAEoIFGgAAAAAkBAs0AAAAAEgIFmgAAAAAkBAnSnFMj2ah8iCeIDVPl2S9uVkGqmSZaV6nwOSPTayXrhZGFV04aJ4uESgddHpM62L8mKvf6ck63bM6KcjpnDXJaiahZ+lVnTozquljqpSqXFufa3hGxydNJ/pcfR0UaFMjXb/LNXRa5u933xP93SWFPbO8Jcv+rcXvybKvtp+RZZmg+92N4Zos2xnp1EWVWnS2eCTrvNXV171e0g/8qfzJU4seRWoWQrYb70uDJT2OVQJtCCGsfSveHx78lE5A3fi6LAppk9R4eFUfUyW/hhBCVof3hXFZjwOVyDYp6jrDmn41ZM197K7pcVyOd4cQgm9/vqPPN2jG601W9Lw2GeuJIb+jr9slhTmnTRbdHsbHcX+q+/gv3/pJWfae5o4sc8mqI5N2OJvr5/aZhZdk2W/ufTj6+/5AJ97OzMt9JOa1EEL4ieWb0d9fy+p3wCOZz2VqYbavq+1+QPe/8mb82Va3dCcaV00nMl3WvR/zbX2+SVHf66FJESwdxBvjxke2ry+gv2zSK7d0vUlRX3jhSL+bioe6XutS/J6Ma7od5dsmFXbJJGibr/ZZTp/P1ZuPTp466uau/ppu/7Ri0ryP9EHTpp+M1vV8/9FKfPyHEMJvHXxQH1T49tEFWeZSub96FP/W60yvyzr8BQ0AAAAAEoIFGgAAAAAkBAs0AAAAAEgIFmgAAAAAkBAs0AAAAAAgIVigAQAAAEBCnChmf1ZIh86FeJz+NGcimxd02UwkjWZMAm5aRLyHEMKorqNLTZKwjdJPuVR/F2H/Svwijq7oLQnKe+aArv27Ot6zu6FjvYcNHWtaaJ08are/phtZfKCfzbSo67ko/VnBxOlWTNmOfga/s/189Pe/dkbnqudMVvD1wYYsO1s4lGXfOz4vyz5UvXOqtnyk+lb0968d67j/DzbuyrI3e6uybD3biv7u2vcopvlUaJ8vnLheeqLP21uL980LX9TZ9gfP6S0xGjfF/iEhhIU3ddm0pP/N7PCKnq7zbRNHvRQfP4UjXcfNvy4S30Xwz7K6Xs7EaY9qJ4/gL93Sc970+Y4sy96pyrLhopmjdMpzGDfMu+Wq3nLlpf34vPG3nv6yrNNt6nHxv23+uCx778KmLHtPSZet5+JjPIQQvtq+JsvWCvHtOc6Z7T6cF1tnZdnX9y9Ff+9MvnGqc/2ZeSYVBs34+7NgtgJafNW8Iw/j7//+sn5Pq++MEEJoPa3Hwbiqx9WkoMsKbfNBZLL7VYT9cEHPa2OzNZKbv2yUvtlCYNTQbcmM9PmWXorP6W6blv453Y55xsyj+/ra0lf03Da5p7ewCFe6sqj0QnxO7D6n+13xTT0PTcz3nP3WNmW5bf1t+XDclGV/UWyBtD/R74FKWl/3iz39zfb1vfg8pLY+CoG/oAEAAABAYrBAAwAAAICEYIEGAAAAAAnBAg0AAAAAEoIFGgAAAAAkBAs0AAAAAEiIE8Xsz9M6YjnXNbGgeR2rWdqJ16vf1vHWvQ0d4Vk41NGlrUs6ijOv00nDqKrXsROd1h5Ss3g0bmVLt3FmtivImXjYBz+pI79Xvqej+w+vmuhuExWcnoh4610TAWuiu3s6id7Gq46Xdb516Y6OuM21dDs/sRSPon9zuCbrfG3vaVlWzZkY5JHuQM2CjuD+vf34VgAhhLBR1NHXD4bxyNl/t/mCrPPF4/fJso/WbsmyX93+VPT3vfG/knUeRWoWQnYQ7xTFA93X+8t6/PdX4v0hM9ZzTf5Yz3luS4/0UHfo7rpuY+O2Hv9ubpjl49dm57WinofKZkuPw2t6Pqk8dPOJLApZE8GfHcTLXOz95KaJm9aXHar3TCy5eQ+MmyYO3ERff/4v/6Po77928ElZ5zsHOub5J5ZvyrKy2VfhRl/Pe72ZHh/juY6G/0TtRvT3P2pflXVeb+l2PFU9kGX/3nJ8e5TPFvQ8+ShSM903p+abJ9fX40DF9rsxMC3q+7z4uv6OOryqvxmKR+ala+R6ul5/OT43uC2amtd1v5yU9HWPzPZBY7OViXt/DBb13Ka+UcwQCPOc2X6jryvmzLdq+JOaLBpe0NeWvqXnobl4POld/X3lZLtmkjV/Mhqs6fdOYU/fr68c6u0+/utzvxX9/ZfNtkM/3/yOLPvGTH8HLhXjWxlkzQcuf0EDAAAAgIRggQYAAAAACcECDQAAAAASggUaAAAAACQECzQAAAAASAgWaAAAAACQECeK2U+P5qF2P55fPC3otd7CTR0jmRWxrIfXdASsi7DNF3RZ/Z6OGc2Y6Ov2OR3hmR7JopAbxqNxh3UTRW/ipruruh1Lr+gI0u66rle/a2Jxl3S9ubiEwaLbbkGXZQb6ualzhRDCwnd11OtEd6HQv6wf3FSc8EZ3Vdb5hTPfkmX/cvuDsuzf2fiuLHuxc06WFUxH6U/1PdkZVqO//8bkY7LO+6v3ZdnBVMfzdkREvbq/j2qeCmEqxnnnrL724oEeI7V78b45aOq2DhdNFP2OLArDRR2lXzJtdJHNk4ouzHbjxyztmznvrG5jaqrHcXlLH7PQOt32IumxLhtX4mWTitvnQBeNa7recEmXFfb0QbNt/Wyqd/W1ffaNvxr9/d/eeEnWSQfdxheOLsiybEo/m81uXZbV8zoGvVnU24R8vv2J6O9rxbas8+nV12XZrf6KLPvv7v5s9Pet0T+XdR6F23ZosKKfa3HXzBt74jmY7jwp6eOlpiePhg/hbbbZ2NZ7WIxrJt6+Eh8j5V39Ljt4Vm/j0LhlvudGZmsL840YUrqscGS+sc7E58u5+cLOHujCbN9s6WF2CUnrJobyPdMY078Gy/HC4r5u46hxum89t4VAyrwHRnXdmUdT3Sf/282/EP29P9Xvv2/3Lsuy9lh/dK4W4xeXNQ+Nv6ABAAAAQEKwQAMAAACAhGCBBgAAAAAJwQINAAAAABKCBRoAAAAAJAQLNAAAAABIiNR8bvI1////41RqN4Rw5wfXHADvAk/N53Odi/02mIcAfB8wDwH4YZPz0IkWaAAAAACAHxz+L44AAAAAkBAs0AAAAAAgIVigAQAAAEBCsEADAAAAgIRggQYAAAAACcECDQAAAAASggUaAAAAACQECzQAAAAASAgWaAAAAACQEP832NiI3sNx9rUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 864x360 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "display_tensor(X)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = torch.from_numpy(X)\n",
    "X_np = np.asarray(X)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Run Neural NCPD"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/jvendrow/.local/lib/python3.6/site-packages/torch/nn/_reduction.py:43: UserWarning: size_average and reduce args will be deprecated, please use reduction='sum' instead.\n",
      "  warnings.warn(warning.format(ret))\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "38e0bcfe59ab4eff96e0f275bfbc1c47"
      }
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "46a3ed98c9024d099a3e28fea674d97b"
      }
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "53c340384e844ab79ae41f2d669ab877"
      }
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "b5ac87997502415d851c2d67f10e5265"
      }
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "fd4c7f4c1a7745e1aa0b667e8d633ce8"
      }
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "3b13f8d4c2814a0180dd40f7473b8f47"
      }
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "1357049e25764b00860c79f56760eabd"
      }
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "97459b65cc6a4aac887dae3c1613a42a"
      }
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "52e06563089349beba89debbcf0373cb"
      }
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    },
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "65fff2ce2b9640c2ae619c3d58eb4eac"
      }
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "# unsupervised case,one layer\n",
    "r=7\n",
    "n1,n2,n3 = X.shape\n",
    "\n",
    "num_trials = 10\n",
    "histories = []\n",
    "\n",
    "for i in range(num_trials):\n",
    "    \n",
    "    torch.manual_seed(i)\n",
    "\n",
    "    net = NNCPD([n1, 4, 2],[n2, 4, 2],[n3, 4, 2])\n",
    "    loss_func = Energy_Loss_Tensor()\n",
    "    h = train(net, X, loss_func, r, epoch = 7500, lr1 = 0, lr2 = 50)\n",
    "    histories.append(h)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAWBElEQVR4nO3dfZBd9X3f8fd37917pV3JCKw1ZpBsQVCcktSx6eJCcV3ixC4wHZhOPRmUuHFcEzqJ8TQTpw2MU7tx0pm6f3hip8SEuoS0k0CI4zoag4c4Nhl3nNphMQbzYEABYkQNWvMgHvSwu9pv/7jnru6u9uFauqt7z7nv18zOnied85Wu9NFvv+cpMhNJUvmN9LsASVJvGOiSVBEGuiRVhIEuSRVhoEtSRdT7deCtW7fmjh07+nV4SSqle+655weZObHcur4F+o4dO5iamurX4SWplCLi71daZ8tFkirCQJekijDQJakiDHRJqggDXZIqwkCXpIow0CWpIkoX6DNz89w29RQ+9leSFuvbjUXH6/q79vCprzxGsz7CFW85s9/lSNLAKN0I/flXZwDYf3C2z5VI0mApXaBHtL7Pz9tykaROpQv0kSLRjXNJWqx0gd7mAF2SFitdoC+M0L3KRZIWKV2gt3vo5rkkLVa+QC++p110SVqkdIE+MtJuufS5EEkaMKUL9IXLFg10SVqkfIFeNF3mHaJL0iLlC/RYextJGkalC/RakehzRxyhS1KnNQM9Im6KiH0R8cAa250fEXMR8Z7elXesRr1V8uyR+fU8jCSVTjcj9JuBS1bbICJqwCeAv+xBTatqB/qMgS5Ji6wZ6Jn5NeD5NTb7EPDnwL5eFLWaRq0I9DkDXZI6nXAPPSLOBP4l8Jkutr06IqYiYmp6evq4jjdqy0WSltWLk6K/C/xGZq6ZsJl5Y2ZOZubkxMTEcR2sUWudFHWELkmL9eKNRZPArdG6+mQrcFlEzGXmF3qw72PYQ5ek5Z1woGfmWe3piLgZ+OJ6hTlAo1YDbLlI0lJrBnpE3AJcDGyNiL3Ax4BRgMy8YV2rW0ZxTpSZOa9Dl6ROawZ6Zu7qdmeZ+YsnVE1Xih66I3RJWqR0d4q2zcwd6XcJkjRQShjorVbLYa9ykaRFShjoLV62KEmLlTbQHaFL0mKlC/T2Y9AP20OXpEVKF+httlwkabHSBrotF0larLyBPmugS1Kn0gV6+/5Qe+iStFjpAr1t3jv/JWmR0ga6JGkxA12SKqJ0gZ62WiRpWaULdEnS8kod6OlwXZIWlDrQfSa6JB1VukB/87ZTeNPpmwHvFpWkTqUL9O2njfHeC98IeLeoJHUqXaADNOutsr1bVJKOKnmgO0KXpLaSBnoNsOUiSZ3KGeijtlwkaalyBrotF0k6hoEuSRWxZqBHxE0RsS8iHlhh/c9HxP0R8Z2I+JuI+Mnel7nY0R66LRdJautmhH4zcMkq658A/llm/kPgt4Ebe1DXqhyhS9Kx6mttkJlfi4gdq6z/m47ZbwDbelDXqhZG6Aa6JC3odQ/9A8CXVloZEVdHxFRETE1PTx/3QbzKRZKO1bNAj4ifohXov7HSNpl5Y2ZOZubkxMTEcR+rUWuVPeMIXZIWrNly6UZEvBn4LHBpZj7Xi32u5ugI3UCXpLYTHqFHxBuAzwP/OjMfPfGS1uYIXZKOteYIPSJuAS4GtkbEXuBjwChAZt4AfBR4LfD7EQEwl5mT61UwQL02Qm0k7KFLUodurnLZtcb6q4CrelZRlxq1EUfoktShlHeKQquPbg9dko4qbaA7QpekxUob6I7QJWmx0ga6I3RJWqy0gd6s17zKRZI6lDbQG3VbLpLUqbSB3jTQJWmR0gZ6o24PXZI6lTbQWz10A12S2sob6KMjzHhSVJIWlDfQa/bQJalTeQN91B66JHUqbaA3HKFL0iKlDfTmqDcWSVKn0gZ6+9b/zOx3KZI0EEob6M36CPMJc/MGuiRBiQO9Ufc1dJLUqbSB3qz7omhJ6lTaQG/Ua4AjdElqK22gHx2he6WLJEGJA90euiQtVtpAt4cuSYuVNtAbBrokLVLaQG8WJ0XtoUtSy5qBHhE3RcS+iHhghfUREZ+OiD0RcX9EnNf7Mo9lD12SFutmhH4zcMkq6y8FdhZfVwOfOfGy1mYPXZIWWzPQM/NrwPOrbHIF8D+z5RvAlog4o1cFrqTpCF2SFulFD/1M4KmO+b3FsmNExNURMRURU9PT0yd00A2jrR76wRl76JIEJ/mkaGbemJmTmTk5MTFxQvsab9YBeHVmrhelSVLp9SLQnwa2d8xvK5atq/Fma4R+wBG6JAG9CfTdwC8UV7tcAOzPzO/3YL+ratZrjNaCVw47QpckgPpaG0TELcDFwNaI2At8DBgFyMwbgDuAy4A9wAHg/etV7FJjjToHDHRJAroI9Mzctcb6BD7Ys4p+CJuadV45bMtFkqDEd4oCjDVqHPCkqCQBJQ/08WbdHrokFUod6JuadV410CUJKHmgt1ou9tAlCUoe6JtsuUjSglIH+ljTEboktZU60D0pKklHlTrQNzXqzMzNM3vEJy5KUqkDfax4QNcBby6SpHIH+qbiAV2veHORJJU70McarRH6QQNdksod6O1H6L5qy0WSyh3o7RG6L7mQpJIH+njDk6KS1FbqQB9rt1wcoUtSuQN9YYTu3aKSVO5AXxihe7eoJJU80Ed9UbQktZU60Ou1EZr1EXvokkTJAx1aD+jyKhdJqkCgjzVqjtAliQoE+njDEbokQQUCfazpCF2SoAKBPt6oe5WLJNFloEfEJRHxSETsiYhrl1n/hoi4KyLujYj7I+Ky3pe6vLFGzevQJYkuAj0iasD1wKXAucCuiDh3yWa/CdyWmW8FrgR+v9eFrmSs4XtFJQm6G6G/DdiTmY9n5gxwK3DFkm0SeE0xfQrw/3pX4urGmnUO2EOXpK4C/UzgqY75vcWyTv8JeG9E7AXuAD603I4i4uqImIqIqenp6eMo91jjjZrPQ5ckendSdBdwc2ZuAy4D/ldEHLPvzLwxMyczc3JiYqInBx5r1Dk4e4Qj89mT/UlSWXUT6E8D2zvmtxXLOn0AuA0gM/8vsAHY2osC19J+a9HBWUfpkoZbN4F+N7AzIs6KiAatk567l2zzPeCnASLiH9AK9N70VNYwtvCSC/vokobbmoGemXPANcCdwMO0rmZ5MCI+HhGXF5t9GPiliLgPuAX4xcw8KT2Q9gjdK10kDbt6Nxtl5h20TnZ2Lvtox/RDwEW9La07vldUkloqcacoOEKXpNIHum8tkqSW0ge6I3RJail9oI81HKFLElQg0MebjtAlCSoQ6AsjdK9ykTTkSh/ozfoItZHwrUWShl7pAz0ifK+oJFGBQAffKypJUJFAd4QuSVUJ9KZvLZKkagR6o+516JKGXiUCfdz3ikpSNQJ9rFm3hy5p6FUi0McbNa9ykTT0KhHoYw1H6JJUiUAfL65yOUkvSZKkgVSJQB9r1Dkynxyem+93KZLUN5UI9PGG7xWVpEoE+ljxCF2vRZc0zCoR6L61SJIqEugL7xX1ShdJQ6wagT7aCvSDjtAlDbFKBPqGItAPzRrokoZXV4EeEZdExCMRsScirl1hm5+NiIci4sGI+JPelrm6jcVVLgcNdElDrL7WBhFRA64H3gXsBe6OiN2Z+VDHNjuB64CLMvOFiHjdehW8nI22XCSpqxH624A9mfl4Zs4AtwJXLNnml4DrM/MFgMzc19syV2fLRZK6C/Qzgac65vcWyzr9KPCjEfH1iPhGRFyy3I4i4uqImIqIqenp6eOreBntlsuhWe8UlTS8enVStA7sBC4GdgH/PSK2LN0oM2/MzMnMnJyYmOjRoWFDvfXbsIcuaZh1E+hPA9s75rcVyzrtBXZn5mxmPgE8SivgT4p6bYTRWhjokoZaN4F+N7AzIs6KiAZwJbB7yTZfoDU6JyK20mrBPN67Mte2YbTmSVFJQ23NQM/MOeAa4E7gYeC2zHwwIj4eEZcXm90JPBcRDwF3Af8+M59br6KXs3G0xuE5A13S8FrzskWAzLwDuGPJso92TCfwa8VXX2xsOEKXNNwqcacowIZ6zR66pKFWnUBv1DjoZYuShlhlAn3j6AiHbLlIGmIVCvQahzwpKmmIVSfQPSkqachVJtA9KSpp2FUn0Bs1H84laahVJtA3jtZ8OJekoVapQD84e4TWPU6SNHyqE+iNGkfmk9kjBrqk4VSZQG/6CF1JQ64ygX70JRcGuqThVJ1A9zV0koZc5QLdloukYVWZQN9QtFy8W1TSsKpMoI8VI/QDBrqkIVWZQN+6uQnAvpcP9bkSSeqPygT661+zAYBn9h/ucyWS1B+VCfTxZp3NG+o8+5IjdEnDqTKBDnDmlo08+dyr/S5DkvqiUoH+k9u2cO/3XuTIvLf/Sxo+lQr0n/qx17H/4Cx/+PUnvMFI0tCp97uAXnrXuafzT37ktfzO7Q/zO7c/TLM+wqZmnS1jo1zyE6/n/RedxdZNzX6XKUnrIrp53GxEXAJ8CqgBn83M/7LCdv8K+BxwfmZOrbbPycnJnJpadZPjMntknq9+dx+PPfsyLx+a49WZOZ74wat84/Hn2dSsc/3Pncfbd27t+XEl6WSIiHsyc3LZdWsFekTUgEeBdwF7gbuBXZn50JLtNgO3Aw3gmn4F+kr27HuZD/7xveyZfoXfuvzHee8Fbzxpx5akXlkt0Lvpob8N2JOZj2fmDHArcMUy2/028AlgIK8bPOd1m/ncL1/IO3Zu5Te/8AC//mf38fKh2X6XJUk9002gnwk81TG/t1i2ICLOA7Zn5u09rK3nNm8Y5bPvO58PvfMcPv+tvVz26f/D3hcO9LssSeqJE77KJSJGgE8CH+5i26sjYioipqanp0/00MelNhJ8+N1v4rZ/eyH7D8xy1R9NMTPnu0gllV83gf40sL1jfluxrG0z8BPAX0fEk8AFwO6IOKbHk5k3ZuZkZk5OTEwcf9U9MLnjND75s2/hu8+8zO999bG+1iJJvdDNZYt3Azsj4ixaQX4l8HPtlZm5H1i4bCQi/hr49bVOig6Cnzn3dN7zj7bx3+7aw+s2N5nccRqbmnWaoyNsGK3RrI/QqI0QEf0uVZLWtGagZ+ZcRFwD3EnrssWbMvPBiPg4MJWZu9e7yPX08St+nO89f4D/+BcPLrs+ovXyjGveeQ6/cvE5J7k6SepeV9ehr4eTfdniao7MJ/ftfZFn9x/ilcNzHJ6b59DskYXvt9//fV4+PMfdH/mZfpcqacitdtlipe4UPV61keC8N5y64vrR2gif/PKjHJ47QrNeO4mVSVL3KvUsl/VyximtZ60/67PWJQ0wA70Lp441AHjhwEyfK5GklRnoXdgyNgrAiwe9s1TS4DLQu3DKxlag7zfQJQ0wA70LBrqkMjDQu/CaItBfMtAlDTADvQsbRmtsGB3hhVc9KSppcBnoXdq6qckPXvGyRUmDy0Dv0sTmJtMGuqQBZqB36fWv2cD3nvfZ6ZIGl7f+d+mic7bypQee4fz//FdsbtZJYD6T+UwyIfPo/HzCoZkjvHx4jrMnxo/upPOxOUse4Nj+9S8fmmNj8aTHpQ95nE+YmZtnbn6e8WZ91f0tZ61NcsWZjh1k64Fl3T6BcumzglZ6clCVnmfZn6cj9U6VPou2QXti6pXnb+eqf3p2z/droHfpyvO3c2BmjkeffYVDs0cYiWCkCLYIFuZHIogIXjk8x0sHZ9m0ob7oH0j7L9Z85jHLAxhv1jg4c4S5+dZ/FJ0bjUTQqI2QJIfn5oni13XzgLVuQ2a5WqEVzO1yfujnucWqs/0LwCV/vr203vHR69IXPtv12Hm/DeD/sBObm+uyXwO9S/XaCFe/40f6XYYkrcgeuiRVhIEuSRVhoEtSRRjoklQRBrokVYSBLkkVYaBLUkUY6JJUEdHNXYbrcuCIaeDvj/OXbwV+0MNy1oM1nrhBrw8Gv8ZBrw+s8Yf1xsycWG5F3wL9RETEVGZO9ruO1VjjiRv0+mDwaxz0+sAae8mWiyRVhIEuSRVR1kC/sd8FdMEaT9yg1weDX+Og1wfW2DOl7KFLko5V1hG6JGkJA12SKqJ0gR4Rl0TEIxGxJyKuPcnHviki9kXEAx3LTouIL0fEY8X3U4vlERGfLuq8PyLO6/g17yu2fywi3tfD+rZHxF0R8VBEPBgR/26QaoyIDRHxtxFxX1HfbxXLz4qIbxZ1/GlENIrlzWJ+T7F+R8e+riuWPxIR/7wX9S2ptRYR90bEFwexxoh4MiK+ExHfjoipYtlAfM7FfrdExOci4rsR8XBEXDhg9b2p+LNrf70UEb86SDUel8wszRdQA/4OOBtoAPcB557E478DOA94oGPZfwWuLaavBT5RTF8GfInWy7wuAL5ZLD8NeLz4fmoxfWqP6jsDOK+Y3gw8Cpw7KDUWx9lUTI8C3yyOextwZbH8BuCXi+lfAW4opq8E/rSYPrf47JvAWcXfiVqPP+tfA/4E+GIxP1A1Ak8CW5csG4jPudj3HwFXFdMNYMsg1bek1hrwDPDGQa2x699Lvw58nH/wFwJ3dsxfB1x3kmvYweJAfwQ4o5g+A3ikmP4DYNfS7YBdwB90LF+0XY9r/QvgXYNYIzAGfAv4x7TuwKsv/YyBO4ELi+l6sV0s/dw7t+tRbduArwDvBL5YHHPQanySYwN9ID5n4BTgCYqLLgatvmXqfTfw9UGusduvsrVczgSe6pjfWyzrp9Mz8/vF9DPA6cX0SrWelN9D8aP/W2mNggemxqKV8W1gH/BlWiPXFzNzbpljLdRRrN8PvHY96yv8LvAfgPli/rUDWGMCfxkR90TE1cWyQfmczwKmgT8s2lafjYjxAapvqSuBW4rpQa2xK2UL9IGWrf+i+34daERsAv4c+NXMfKlzXb9rzMwjmfkWWqPgtwE/1q9alhMR/wLYl5n39LuWNbw9M88DLgU+GBHv6FzZ58+5Tqs1+ZnMfCvwKq32xYJ+/z1sK86FXA782dJ1g1LjD6Nsgf40sL1jfluxrJ+ejYgzAIrv+4rlK9W6rr+HiBilFeZ/nJmfH8QaATLzReAuWu2LLRFRX+ZYC3UU608Bnlvn+i4CLo+IJ4FbabVdPjVgNZKZTxff9wH/m9Z/joPyOe8F9mbmN4v5z9EK+EGpr9OlwLcy89lifhBr7FrZAv1uYGdxxUGD1o9Ku/tc026gfWb7fbT61u3lv1CcHb8A2F/8KHcn8O6IOLU4g/7uYtkJi4gA/gfwcGZ+ctBqjIiJiNhSTG+k1d9/mFawv2eF+tp1vwf4ajFq2g1cWVxhchawE/jbE60PIDOvy8xtmbmD1t+vr2bmzw9SjRExHhGb29O0Pp8HGJDPOTOfAZ6KiDcVi34aeGhQ6ltiF0fbLe1aBq3G7vWreX8CJzAuo3X1xt8BHznJx74F+D4wS2sU8gFa/dKvAI8BfwWcVmwbwPVFnd8BJjv282+APcXX+3tY39tp/Yh4P/Dt4uuyQakReDNwb1HfA8BHi+Vn0wq7PbR+9G0WyzcU83uK9Wd37OsjRd2PAJeu0+d9MUevchmYGota7iu+Hmz/OxiUz7nY71uAqeKz/gKtK0AGpr5i3+O0fpo6pWPZQNX4w355678kVUTZWi6SpBUY6JJUEQa6JFWEgS5JFWGgS1JFGOiSVBEGuiRVxP8HuwFRCzr3qjAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot the loss curve\n",
    "histories[0].plot_scalar('loss')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "S_trues = get_synthetic_gt_topics()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Topic Modeling Loss 7 to 2 (average, best, worst): 0.0870316914580482 0.03349681369449985 0.3606750570744094\n",
      "Topic Modeling Loss 7 to 4 (average, best, worst): 0.08746390190365984 0.03403409599699758 0.29721958867942444\n",
      "Topic Modeling Loss 4 to 2 (average, best, worst): 0.08067261694603298 0.02965444532202688 0.3942359231462281\n",
      "\n",
      "Relative Recon Loss Rank 7 (average, best, worst): 0.2774291903373256 0.2774291903373149 0.27742919033733426\n",
      "Relative Recon Loss Rank 4 (average, best, worst): 0.3257536010491811 0.3136097821832106 0.43503309554226227\n",
      "Relative Recon Loss Rank 2 (average, best, worst): 0.5011643155076325 0.5011593577809417 0.5011748285664029\n"
     ]
    }
   ],
   "source": [
    "all_topic_losses, all_recon_losses = measure_all(histories, X_np, S_trues, kind=\"Neural\")\n",
    "handle_losses(all_recon_losses, all_topic_losses)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Topic Modeling Loss 7 to 2 (average, best, worst): 0.0870316914580482 0.03349681369449985 0.3606750570744094\n",
      "Topic Modeling Loss 7 to 4 (average, best, worst): 0.08746390190365984 0.03403409599699758 0.29721958867942444\n",
      "Topic Modeling Loss 4 to 2 (average, best, worst): 0.08067261694603298 0.02965444532202688 0.3942359231462281\n",
      "\n",
      "Relative Recon Loss Rank 7 (average, best, worst): 0.45432993797408994 0.454329937974077 0.4543299379741002\n",
      "Relative Recon Loss Rank 4 (average, best, worst): 0.5075737874199813 0.4929677600630872 0.6377774258872035\n",
      "Relative Recon Loss Rank 2 (average, best, worst): 0.7139180969705932 0.7133916204044575 0.7145864213797999\n"
     ]
    }
   ],
   "source": [
    "all_topic_losses, all_recon_losses = measure_all(histories, X_original, S_trues, kind=\"Neural\")\n",
    "\n",
    "handle_losses(all_recon_losses, all_topic_losses)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "def standard_HNCPD(X, ranks, seed=0):\n",
    "    \n",
    "    np.random.seed(2)\n",
    "    factors_tl = non_negative_parafac(X, ranks[0])\n",
    "    \n",
    "    X_1 =  factors_tl[0]\n",
    "    X_2 =  factors_tl[1]\n",
    "    X_3 =  factors_tl[2]\n",
    "    \n",
    "    factors = [[X_1], [X_2], [X_3]]\n",
    "    \n",
    "    for i in range(1, len(ranks)):\n",
    "        r = ranks[i]\n",
    "        for j in range(3):\n",
    "            model = NMF(n_components=r, init='random', random_state=seed, max_iter=300, solver='mu')\n",
    "            A = model.fit_transform(factors[j][2*i-2])\n",
    "            S = model.components_\n",
    "\n",
    "            factors[j].append(A)\n",
    "            factors[j].append(S)\n",
    "            \n",
    "    return factors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "all_factors = []\n",
    "for i in range(num_trials):\n",
    "    factors = standard_HNCPD(X_np, [7,4,2], seed=i)\n",
    "    all_factors.append(factors)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Topic Modeling Loss 7 to 2 (average, best, worst): 0.14537284932010094 0.02899188999856203 0.2746495389008453\n",
      "Topic Modeling Loss 7 to 4 (average, best, worst): 0.19297659518111737 0.042224271170372724 0.31005236045224255\n",
      "Topic Modeling Loss 4 to 2 (average, best, worst): 0.2040519538015505 0.043877446441610034 0.32078403987945087\n",
      "\n",
      "Relative Recon Loss Rank 7 (average, best, worst): 0.2774291903373288 0.2774291903373108 0.2774291903373438\n",
      "Relative Recon Loss Rank 4 (average, best, worst): 0.42807401228868774 0.4074058430428231 0.4600117730144047\n",
      "Relative Recon Loss Rank 2 (average, best, worst): 0.7162493143312865 0.5729056429653906 0.7980806524269276\n"
     ]
    }
   ],
   "source": [
    "all_topic_losses, all_recon_losses = measure_all(all_factors, X_np, S_trues, kind=\"Standard\")\n",
    "handle_losses(all_recon_losses, all_topic_losses)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Topic Modeling Loss 7 to 2 (average, best, worst): 0.14537284932010094 0.02899188999856203 0.2746495389008453\n",
      "Topic Modeling Loss 7 to 4 (average, best, worst): 0.19297659518111737 0.042224271170372724 0.31005236045224255\n",
      "Topic Modeling Loss 4 to 2 (average, best, worst): 0.2040519538015505 0.043877446441610034 0.32078403987945087\n",
      "\n",
      "Relative Recon Loss Rank 7 (average, best, worst): 0.45432993797409366 0.45432993797407173 0.45432993797411203\n",
      "Relative Recon Loss Rank 4 (average, best, worst): 0.611797391577996 0.5846297844192012 0.6468655065126846\n",
      "Relative Recon Loss Rank 2 (average, best, worst): 0.8920543406648231 0.7435496281173065 0.9905674873168752\n"
     ]
    }
   ],
   "source": [
    "all_topic_losses, all_recon_losses = measure_all(all_factors, X_original, S_trues, kind=\"Standard\")\n",
    "\n",
    "handle_losses(all_recon_losses, all_topic_losses)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "X1 = np.asarray(histories[0].get('A_X1')[-1])\n",
    "X2 = np.asarray(histories[0].get('B_X1')[-1])\n",
    "X3 = np.asarray(histories[0].get('C_X1')[-1])\n",
    "\n",
    "\n",
    "A_A1 = np.asarray(histories[0].get('A_A1')[-1])\n",
    "A_S1 = np.asarray(histories[0].get('A_S1')[-1])\n",
    "B_A1 = np.asarray(histories[0].get('B_A1')[-1])\n",
    "B_S1 = np.asarray(histories[0].get('B_S1')[-1])\n",
    "C_A1 = np.asarray(histories[0].get('C_A1')[-1])\n",
    "C_S1 = np.asarray(histories[0].get('C_S1')[-1])\n",
    "\n",
    "A_A2 = np.asarray(histories[0].get('A_A2')[-1])\n",
    "A_S2 = np.asarray(histories[0].get('A_S2')[-1])\n",
    "B_A2 = np.asarray(histories[0].get('B_A2')[-1])\n",
    "B_S2 = np.asarray(histories[0].get('B_S2')[-1])\n",
    "C_A2 = np.asarray(histories[0].get('C_A2')[-1])\n",
    "C_S2 = np.asarray(histories[0].get('C_S2')[-1])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### All models, Rank 7"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Relative reconstruction loss: 0.2774291903373317\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA2gAAAEkCAYAAABaExIDAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAnaUlEQVR4nO3dTagta3oX8Kpa++Pc2x06X03wI9hTByYGojgN6iCtEsRkoBPH6YlzJxkbFGcOHGjAFhE0oETaYMSRtgREQYMYUSIoRPLRdu7HOXvvtVY56KQlcT3/c9dzVq3z7nN+v+F+71v1Vq2qp+pd657/O6/rOgEAAPD2LW97AAAAAHyLCRoAAMAgTNAAAAAGYYIGAAAwCBM0AACAQZigAQAADOLmnP/4e797t37p+2+3Gsvv8l8fvlA3hpUB5rluXNe5bHv85t1nGdbmln1orIcfz0k47Cmcrtb+bj49hE7PhKUnNvXy6ZvT4/7TdHVF16xD1/ZLv/rFtz2EaZpeUxeegduP39M61K7noeMSvse9dK3sjr9BHRpDt+Y1X4faOvu7++gdqENX1H5XjRtNL8eN7V2xDp01QfvS999Ov/jz33+ZUb3Gn/3lHy3bjuFTvFvqG+LxuCvb/ufPfekzjev3al1Q4aJ48Zthghmek3OoA8cw95wvPCH83n/3jdApOPa6bfEb8PxQnJTlwnfmNE3TMb0YbbC/jguP8ev//e++wWCuW4eu7Qd/+it144W/TEl2j3WnOdyrqR6u4VqZw0M09juc7vd9/+Y9rUOHcAC7MMh9/QBZ7+sHyHwID565Mc70MpW2l85JUb++/is/U/f5DN7lOnRN3ZqXXthTHeqK+yturd//C7/Z21mqJ+k+SG37UBtuGl/CpH0dw75CvzWMY37q1Zr5MbzkpnrZ2FessYWv/7e/U7b5XxwBAAAGYYIGAAAwCBM0AACAQZigAQAADOKskJBrermv05FSSMghJE49HOrDPaQgjS0SEq/pmv/QtvsPWLv/Oj8ljCXpH6NX//D9uMG/PE7j2GJ/HZceo5TM0vJ4vX2lYI7+Nuu2a37qVXjIG4lpYOEfvqd/wJ7CUW4bj+ddHYSVzM06ut729hdDCSrdAJdqiPF5xLWkmtetUTFYKIRexXs19WuEZcVjC3kYa7rl4vjrG2hNv9V0PoN0ry7pfeL8XU3TFMcYa1S3fl1S+Mz8ggYAADAIEzQAAIBBmKABAAAMwgQNAABgECZoAAAAgzBBAwAAGMSwMfv3u33ZlmL273Z1PunczL1vxc2/Qb/O9uKurpgmfPygXh6hGzc9pXjb3eW/Y7h5WeT+bhHLnI57lBjoS49xlOMaUXfViMZyH2soDJvE1F9TiJTeRLhH5ub1Pu9D1nZjHFOqlYf6fM0p7v8hjLGzv2vWQ8t9DCHGxqfClj6+1G3Xu47S+8vaqdvt6y8dXIruDzXxHb4X5qdQo6p3y24d6ryPhn35BQ0AAGAQJmgAAACDMEEDAAAYhAkaAADAIEzQAAAABmGCBgAAMIhhY/ZTlP4hZJoe1zpKNPWbu6nMKbG1SM9spv23+yWXXwqgtxZAivyeUizuBontVbz1ulx+Z3NaQmCD/XVcfIzvcKTvmzp2K3IzebmyxHjrFEGcNpo22YywrvqNtCRGauvG4nc0xxiXR7lJGelBp99zWJKEs6WY/S2S6Lvx/OkdsTyGTu2aXvcelWpleMcNz+rYr1qyJI1xad6r6Xwdm/d4OLbRf6IafHgAAADvDxM0AACAQZigAQAADMIEDQAAYBAmaAAAAIMYNsXx4VAPLSV+xW3u620u+9YmYwpimboYAm4O971xpE8ypSTNYYp+DOmJVbLf8nQ6AfFbnZppZqktJfQ0rVXC2AZJYWs37eiKLj7GQY5rRMtTaGymj8V+DSlNNu6qmZSb0tPKhNFDN5Y36KYIpsTCzjhTcmraXkqFLJJrp2mapvu7uq1KeJumXJur/XXPcTonIYWWt69d8+JGQ1u3DoVrcz0U11/aV9heLOfpFkmNoTbEfqPo1oZUo6pz0n3nTDW28d7jFzQAAIBBmKABAAAMwgQNAABgECZoAAAAgzBBAwAAGIQJGgAAwCCGjdm/39W59ylm/3YX4oKDYzgTMVb6wvHWu1dhcykuPxz2MaQkp4jbJcWrV00hSrS1JMFrtrmJFMt6YWVM+DRNa4qOjhttLmdwLSOMYVDzoXluLlyH1rDERqx5qa0Zfb2mfp0lV7qR7N06FCObL1zbUsxzEpYCiBHc3VrTGWe3blz6HHNZ6fIK92N6ZziGzzxtM2m/B3bEa715z136vWCL53j31WuLJVAG4Bc0AACAQZigAQAADMIEDQAAYBAmaAAAAIMwQQMAABiECRoAAMAgzorZ/+VXX5j+9H/+cyfbjiFn9OlQR1k+Fm3/9o/+o3OG9sZ+6Ge/cr2dhUTQp8+HeNgUwZ0SVEPk7CFE8Kd46yrWf/nkIXS6cvx7NxZ7V1yvKa75EPJhQ6T0uq/XR1jv6w9nPoR1FZbwwVX7u3Qk9jTVkeXXXjbhGZm7McMXjtlPmdJxSYxmFHWMcj82IrOvHa2+xf46ceDd+zjVr1RPmnUv9uO9ku79dj2M0f29Taa6N8ei+MxVn096jqdznJZv2kJ8tlx4LKnmNd5x/YIGAAAwCBM0AACAQZigAQAADMIEDQAAYBAmaAAAAIMwQQMAABjEWTH777IyrnkLKYr62vHQKTo2tp0fvbqGti2OOu0vEQLP27LWK5Lkfune6ixh0YypjlH6ocbmlOpGDHeKhk82+MqyXYeO+0anFH39DkeB82zF+2OXXpbSRkNTt8Ze+r2t+a7UdtM8cN4av6ABAAAMwgQNAABgECZoAAAAgzBBAwAAGIQJGgAAwCBM0AAAAAZxVsz+Mq/T524eT7YdQzj541LHe+6WMZL+58PbHsG3xAjYkLwaY15TYmtKZQ5jWYo87cPn7sLOwjCaCdAp1rtr981GvHWSInO3iNNNquskRKdffF+UjjeXP2drcZPHZTT2dWO855ptMTE7LSFQHURnaYGNtD/Ra94/cemEcP7jcye0qQ38tvjOs8EDvqqHn6Fj3dT5qSPUqDTCOP64zEa4yedGBH+8v9PDZYN7v/uO1fncuuNv9PMLGgAAwCBM0AAAAAZhggYAADAIEzQAAIBBmKABAAAMwgQNAABgEGdl3B/XeXp1ON1lDXGojyGv/fHQiPd85rqR8inmdT6GWOxmpPGcUlmL/S0PIaI+xYw2Y7HjkXXjUDvn69px+dfeX2WDz/R9lu65tsbHkGtG73ON20ybDF8jlpfYYYsT2XQTnnH7xvou6Tym496FE5kiuLvS/V89r1Kfbmx/eDby9sWa1/3o0k8PzUt9js+zQZ7HyUg18ZrS59Y5J6kOpRrb4Bc0AACAQZigAQAADMIEDQAAYBAmaAAAAIMwQQMAABiECRoAAMAgzorZX+Z1erE7HaN+DDGj86GOuVzamfOXNcgwri+lGnfOSYo0jbHLjX1NU/6K4Zox79eOlO/GUV9zHJxt3eArs859nGp2t1bGSzY0piVcyvqVou27mvfcGtrmTgR/io3f4H6M8eJpLLtnED3OW5dq3ibvZc0am+pQq25fe6mcbk3s1JT0Prds8KF2616nX/dza+zLL2gAAACDMEEDAAAYhAkaAADAIEzQAAAABmGCBgAAMAgTNAAAgEGcFbN/WJfpm48vzt7J46GO93wKbde0Rbx1SzP2PgZ4dlNNU79LJ6Wm89+N4O/GoVbR0YvYaLZ11eU+tkg7HuUWqSLq34IYpX/sFrdqZ1f+ANRE3lCsed0alS7LtM1mv4svSbTFfXxo1ppqLGmMKUr/2jUqGWksJ4wyLQEAAHjvmaABAAAMwgQNAABgECZoAAAAgzBBAwAAGMRZKY77/3M7/ebP/YGTbSnF5hj2snt1+u8/9LWvnDGyzyiE2Pz7n/pbl9/fe+jLP/LjdWNKM0uJRimRLaXw7ML3DynR6FDtrx7/endbby+Yd/Vxz+U4XiOcr/X29M04T2FfVarl66Tz/wZ+6Ve/OP3gXztdH+bwsa6hDi0Pp/++//wZA/uswun8j391jDr0I7/0Y2Xb/W5ftj0d63vkGCIeP3f7WLa93Nf3VjmWf/Vh2addh57q456WcK13a1u1zZCcuO7OeqR/W8wyu3TS5BTqULfmJVvVof/9xekH/sbpOnQMj4OUdFrVryVceusWQdgbJB1WRql5W/jyv/yJujG8u6y7dJHUbfNjfaEcvvBB2bY81P3Wqg6FMS4vwztIuF7T/b/e1h3nh9Dvvr4Zq/O8vHwq+6R6Hsf4dH5t8wsaAADAIEzQAAAABmGCBgAAMAgTNAAAgEGYoAEAAAzCBA0AAGAQvUzeU0K8aorgr+weetHenQhbLmd9EfKFUyR+iqIO0aVxm03zxy8vvs1SirAPcdpROpedcbzH5pSM24yifu51aA1FNkXpd/slVb8lROmvKaY67ayIhv9Wx3DcaX9peY6HaumBUA+7kevdOpT6reFCvyu+G+7WoVTzNorZfw7a9WsDnffAZy8tv5Gk95puW7qPN3iPKqXnXzpdz/y52fH+Vi4AAIDBmKABAAAMwgQNAABgECZoAAAAgzBBAwAAGIQJGgAAwCDOitmfp/c0KpXPrhPx/rp+MZY1re/QjPV/DlIc9TUjc5ONzvEWdWguxjofu0sdNNsGESPxQxh96pekmP1OBH+Ktk9fS67p8+7eV839zd2I7lGMUoe2dPE6dNn9xCU9umNPH+szqG3XFOtQ7NhdZuqK99xzf4fawoXPyTN/AgAAALw7TNAAAAAGYYIGAAAwCBM0AACAQZigAQAADMIEDQAAYBBnxexvoYq33r9oxoWmlOQUOctFHF/c1o1L+HBCbHyKVI8J3Gl/we7Th1a/q2qey6t6CzHb8VppnJbjrj+WynNYquTh0Dvwp2P9nV+Ky98tdXF+bIzlg/RhH0LHbkxyutbT/pJRYqy79SSMv3rud7c3zLn6HWkVlOuNItavbh1Kz9znUNuuaT6Gl85QM9bmVRLvq3Qfp3u16rLFlZze0eO1XI9/jbVh7KVA/IIGAAAwCBM0AACAQZigAQAADMIEDQAAYBAmaAAAAIMwQQMAABjEWTH761RHrKawyhTLuhZRozcvu3HHaSC9TfLZ7T56VTcu4fuALSKUU/T1aLHMl/QuH9v023Wo8dVS6rMWyxYszYj0GEX9DJb7+PzdY9l2G07KU3Ndgg9unsq2m3DC7m/2J/++P96VfdZdPcYYurwPF0Na9iLVvRTDXdli+YrmkiRROLTquT/H5UOa49hwuY+qpnTq0zRN9VfmqXalSy+lvG/wihW924+ly0rXbPd6TvdWivyv2japQ71u5RinqXcuu++OF/7c/IIGAAAwCBM0AACAQZigAQAADMIEDQAAYBAmaAAAAIMwQQMAABjEWTH717T20pqffbz1c7fe3zY7bhCl3zQ/nY7u3sQW8dadcxLjrcfKR56P03T7yekxpXjrpU6OLyOnH7+Q4nbrpihs8k/8hx8v23YhM3sXMrM/enV/8u8v7upo+6//4D8u256DH73/S2VbjGROdegmXFydmOrXWD5+eXaf+SFc5Fu4CQ/qcC7LGttZdmCa8lIGG5kP03T3W8UxXrisp7qW3mtefXc9kBSzH5cJ6D6qx3qMXMV613vF7taM6ba+H5dXvfea5aNPz+8UljJZwvtVtdzNNE3THO7x9M623tafwVzVqOYzYk5LsaRaWfALGgAAwCBM0AAAAAZhggYAADAIEzQAAIBBmKABAAAMwgQNAABgEGdlgM5THZWaol5jWxVZGSJZY1yrmP1hpejYGCrbjODvRtVuEHxf60bYp1j8zvkaLEo/mdc6Mj/FQ6elO5YipTfWjA1i9l/c1HHBKWb/Jgz0sYheTvt67uanEHec6tAGtaZdT6pY6XTvL831adL9310KZNf4/jeNf4sxvoF5qutGuzYUDi/qtrR8yBJug7buu9l7KNah2LH57tJdriht81A8W9LSIs0yFOPy7+/qjine/kXoNzi/oAEAAAzCBA0AAGAQJmgAAACDMEEDAAAYhAkaAADAIEzQAAAABnFWzP465Rjrsl+Kvi7icTv7maYpZxqL2d/cukvrHIRY1pRGmyJnn/tXDFvEQ3cietM4BovgX+dpOt4WjWn1gXSIRSVs16EgjeN2V98IuxClfxMi+O+Kbd5tksE9hnYdai+dsMF9/FB9pumBusG9ug8PzpuQp13Fc09Tfb7S+OOxhXOSroU3EN+HLrzLtIzDelOfl2Mz7jzVvbT0yFhPircv1qGkWU+2uP3nLWpbvbPLbzPVoc44ujUq1crCc3+9BQAAeGeYoAEAAAzCBA0AAGAQJmgAAACDMEEDAAAYxFkpjrcfH6ff968/Otk2PzzVHZd6HrjenG6bH/f19lIoS7G9aZpiGtWX/+RP1Lt7cdfa3/Lxq5N/X2+b0UrhPEbhuNNY5qc65a3q99V/+rfLPh/OVfTeNP368bFsS0f9ubluvZ/ry3sXUnr++E//lZN/P7yox7GEyz+ZY3pl3RTTBTshjvXpzyms4VKuju3x184qO/+fm5fH6Xv+0ycn25aX9QexpvunqkOhrs3hvjre19f68iqc7H9Yn5vjFz4s25529bF9/psvT/59vb8v+/zo/BfLtngeg/lwfj2Zpl4d+mdf+3tln12oGR8fT9fsaZqmJVSiD5fwjGj6gb/+lZN/P9Qf27SEx2aK2EuBnil1NNah1FaMJdahUGvSOKpje/z7b1aHbj85TF/8xW8UO60HdPigrg1V7V6ewjM8fD67T+oTmpIhjx+G+vUYLpbj+cmdX/7ZHy+7rC/CuQqJecfQb/dRfY+voW4nVVrjP/naz5R97sP70BYe1vpZlmrbX/6VP3Xy7/tmxPFvPdYvUt/74vRzfZqm6dN9fb4+3df193/98+8r2+aiXs7dNMwU8Ficrsev1nXIL2gAAACDMEEDAAAYhAkaAADAIEzQAAAABmGCBgAAMAgTNAAAgEGclzM713Gic4hsTXGuVfRqitKfQ5TrGuacqV85jmmK09gUcRu3eU3daXij34u5zkJOsfcv5hQBW0vbzDH7KZf59J9j3Hxoi5Gt8frptcVrsrOvbqT/IJf/t13666ju/Z363aUI7lSjUtZ2ceCj1KfXaYxzP9VR4LtwITyt9TNiCTfy01rv7zbUxKvqftxb3OPP5NJ7DtqR4KO4SWu1NC+UVA+by4R03mNTXUi2iOB/tdZrcOzCDfnqcPo96vHYW6biVYjLr/b1rba630Pol96HLl6GLlwr/YIGAAAwCBM0AACAQZigAQAADMIEDQAAYBAmaAAAAIMwQQMAABhELyfzlPXCWa9p6hhyM1MU6sWXApimGOe6VnGuzejYOMYg9krnOSXEFv1SPGyKtr8PbSkCth2lH1w6pj7GvIaVH7qx/p0417iv7rFVbW8rYjuc66t+VZXu48d6uYlpbUYvX7o2X9szGP9yzYt6g121at4W3uVj47R9eNG4DRH8qS4cU1t6ENT7m8P+OiVqufLvI+k9Ko3lxe50PH9adiR52tXn+G6pr4V9aDs0l07oiO88V5wGAQAAcEUmaAAAAIMwQQMAABiECRoAAMAgTNAAAAAGYYIGAAAwiMvF7Mco+stubw3Jq8m6u3Le7nOf/jZi/btx090I2G6UfosI6LNsdWzz/jDtfuPjs/sthzpeed2dvo7Wz72o+4RM3RjJfFNfs/NUR+kvr07HHU/TNK0PIYL4O04fw+6Tx7JPivtfnsI4butHyhzitNebXlGfPz19ntNyH8l37T5s9dvC7cenj+02XPop5jktpTHXH2nc5uG+blvCNqcLx1Efw5vMWiyFk5Y4+UwOh2n55ien20KU+PKNtH7Nad376vidnw8brT+E5fH8MX6r4/nP43hsT71xxNoWxtjd37Q/fY39+f/yE3WXYz2Om6W+OFO/73lRXI/TNL061DXx6VBfrz//h3+ubHsOfvgf/GTdWL2jpPqU3mtSTSn6zeGSe+5TCAAAgHeGCRoAAMAgTNAAAAAGYYIGAAAwCBM0AACAQZigAQAADOL8mP1jkT8ZIltz9OTpfmuKQk37StI2H+pY6fU2REBX52OqI1vj9oL2cSfps2l8pvupzgzdhe8DntZ6IEvIeX5a6/3dzt3zXDSk0x/aUkx1jHoOca4xwr5xmXTHEY+72uYbX8ZzawmIVpT7sZnF3bxX18ZxdaV9vcOrP3DKhevJJka7KOdlWu/vzu/XXFKist6F79mvXL8u7Zr18HU671/3u3qtiWWqr4NdiNlfwg15t4QlF97lNX3eUX5BAwAAGIQJGgAAwCBM0AAAAAZhggYAADAIEzQAAIBBmKABAAAM4ryY/TVEhoc013mpY0Gr5NIqon6apmk+1G3TXX1I88vHsu3wnR/W20yWOrr08B0vzt/cqzruvy0l7abE30a/j47hHIe8+Y9SzH5K+58fyrb7ub4Wdim+txhKiqKPbc2Y/ZSKe/H04ZQgnI4tdKuO7c2HvrZioOdDfSDr7vR3VccXt/X2NkiiXj6pr+f+Rk+f8bkZwd1aruAN+l3+YuedEwvRVjtdp3lfLKUT3gvmsDRPuafm9g4v6veaVL+6ieydbe4+3qDmXVkVwf90qGve47FuuwkX7f5Y/66Stvl4rN+H0jZ5e3wqAAAAgzBBAwAAGIQJGgAAwCBM0AAAAAZhggYAADAIEzQAAIBBnBezP9dRqXN3qldEKK+3IZL5pt7ZmiKZP7wvm5aHfb3NtL8wxy0j89/h2OhdiIdN0fbp8knh3EvYZt7fhb+bCB9pDFXeIh76ipdXjGMu2jZIp99Mqgsx6j/d46lfEdv9WnHZiGJ/W+yrq7vNd7iWrkWJilHmaXvNZTvWUIDT8/a4q0fTWqIiRringzvz7xtrReaH46uuk2nqL9szd+tXZ5th+ZN2jd1CYyy3u7BcVLgJbpb6nOzCy/aLXXheBY9hSaLnLtev03/vLkERL8lGHfILGgAAwCBM0AAAAAZhggYAADAIEzQAAIBBmKABAAAMwgQNAABgEGdla877w7T7xidn72QOMarr7vQccQ2R+HFfzQjYdannqvMhbPNYH9vxxe3Jvy+PIXr104d6X83I2c75f5351emxfNfuw9b2/uDN6XP1Nty8PH1sNy/rPnNKCg4xr3NIxU3bPLyo25aUtBu2WQ8kbC5UkCpOuhWx/bs2ME/TTXFSq5jqaZrW2zDY4h6Jy31sYN53PqBs93FRU8K9v+7q454fQ3R3s57EfOJG5PcP/9RP9saRpOs2DbGZBr7+mW+c/PshbHAX4rmTx2N4/oV+NyFG/HZXj+XpcHp/u6U+yQ9P9f2bxlGl2s9f60WS/78NzNN6Xzy30hIW4V1jnYp+oc90qPe13o0fnz4/NZf7eAZ+9atfavVLSyek5+evpTLafLT8sUOjljbrYXp3Se8a6diWv/Ab9TaLcxlWwog19nHfeF/4hfqg/YIGAAAwCBM0AACAQZigAQAADMIEDQAAYBAmaAAAAIO4XMRPMz2xlJK7trCr9xcPLaUnFv3WsK8tjrpK03sj1/58ruiNUwZ/rxS+2d1X6nfh/cXLP47j0ify+mLCWDNVNUrpb2l/KeWt6hcSaKeQ4tjWPCdrrLGnj+1wl+LM0s6abRt81bkUaWFpV91Sv4ZrqxrHNOXUxZR2diyi3FKf3BbGXxWpd/cR9m3PIiHx6Q3TNAfWrkNbXJvNFMebIrE7icm1G7Sl1MubUDfmYqBlzZheV4fOfxCkQ/YLGgAAwCBM0AAAAAZhggYAADAIEzQAAIBBmKABAAAMwgQNAABgEJeL2Y8Ryo247RQBvYHl5eXjaA/fcX/y77tPQ6xsN5479dtim+9AhHrlWCSM50j50JQS0EOa+ZwS10Oe9hoipzvxvSky9xiWjCi//nlG8dbrXSiRG9wD82O4WI4pZj99SMXft4jST5rLElRR+nlfzbYgjWM9NuO0g+Px9LVQRdRPU46HnkPbIWzzeEjXSXpu1v0OxflaQ0E8FOfjdarjfocfYd+23l55uYzGSZ1vwhi7NS9p19Hz+8XaFd8Z0jmum1KtadXRtL+BnuPp3SzVjeojTUu7pPpb1bXkyqu3AAAA0GGCBgAAMAgTNAAAgEGYoAEAAAzCBA0AAGAQJmgAAACDuFzM/qUza7vR8E0pSrytyPBMEZ4pCnXe15HGa4iqnQ/1kgXrLmXApxzYdzejuIxsTSnh3VjvdIq3iMxtSMcWo8erA7jE0Kt7KH3l1KkpW9wDaRwpynkUWyzbEXTq5bLfYAmEeK9edmmLaZqm3c3pel/eV1OO0k8R/FNY+WUXlu242dXPpF2K9S9ixOMYg91SP+Oqc3LlV4ztPJflcC49lmdQK5dwX8Xnezq29jIhvX4Xj9NPl2vaV7PfbahRlbhcSbNfFc8ftxfaAAAAuCITNAAAgEGYoAEAAAzCBA0AAGAQJmgAAACDMEEDAAAYxOVi9mNm7fn5nuvNdeeO83p+FGfbtafFI0XtPgNzndhc2yBmfwmXZErhjeMv+sUI8XS9do/7TVXXdDwx5w9ofgofwhYx+91lL9I2L11v4nE3I78vnHm+PIXG5qOqfY80VXH66RLfxRUcLp8rnyL/0wkr+8WY6npX9SIz74G0NE+qX13de7UaZ6yHA62F0FjCJdWhbqR820C1rTWODd4nyhob+qTlQy5dY/2CBgAAMAgTNAAAgEGYoAEAAAzCBA0AAGAQJmgAAACDMEEDAAAYxOVi9i9s3oegyw3yNtfd9YJ6077mEDm73vTGGPt1I3MvHIs9klasbDOmNkXiH8PHti7hOllC7HKxv3aU7rWjgq9sXeoTM197+Yp0z73D92PnPO8/uPz5SONY42fT3OHh9LWXopyPIYt+DvHQT4e62BzCNtOjJY3zcDzddjjW99s+jDHF/VfH/V6sPhPqV1e81oP6/umsa7ORC78P7T+ou2xSM4I5rVMR7F419hXa2vddc5mAx/35783h9Wo6hJ3ti5qdpNPhFzQAAIBBmKABAAAMwgQNAABgECZoAAAAgzBBAwAAGIQJGgAAwCDOj9nvZGSmPp3tbRApPT8d6rZuvPLN6fnv8riv+xzryNl5X48xRel3+yXdyNbnIEXfl5pJwXFfKVY2xWl3ttn8OPO52vAaqe679JVTo27M4X7cxBZx+Zfe5hZx/2l5kdBW1eZlf/lrL13raWmLrtvb08+JFF+/W3rX6y6MP20ztd2GtqciTn8XIqzTMgG3u/D8Lv6+pGzu5yTdcxusJXDx5UWeyxIhjXEu4VUvPR/DLd7WvdxbS/CkSzJsLx33ml5VQ9n7oKij01R/Aun034S69rCc/z6dLiu/oAEAAAzCBA0AAGAQJmgAAACDMEEDAAAYhAkaAADAIEzQAAAABnFezP48T9Pd7em2FEV/H+aBRQR8ir2P0bHNyNb1rj4V3VDZuYjTX4v4/W91uq/bnkJm667e5lp9ZtPUjreulgP4I3/zK2FfdVM8yamt+RVDinP94Md+/eTf08oCKaY62R/qA0ix0ne7Ouo1xcDuq3jr0OfVU31v3IZxVOOf/0XMHn69eZ7W297yEKXd6QvieF8fe4ot7sYk715+GjYadrik7OKiX9re/V1vHFvEel+4U/ps4meaTvEGX3U+7s+/xg/H3oW3P9T72i31Nm+K+2aa8nIA1f5SHUq1cg37qupQ6vOZzHP93pDqU4zFL7bXfK+5dv3qiDXvmYvnsVmj2porxrRW7uiuthL2NYcpQbqWH0Nt6zg06to01e+P8bH+WQcFAADAtkzQAAAABmGCBgAAMAgTNAAAgEGYoAEAAAzCBA0AAGAQ58XsT1OdCZkyyFMGZiVFyneFKOo1xAVHMTL3/OOeDyFL9BlYwscW461TlGszZj/uL2yyoxvZ3ImHfu02G31SJHa332bf/qxrvQxHN4q+WDZid7hupPx0e35Jfja6S6CEfnO5hEAYRmsUr3mMbXApLMXSHVss97GudQFO21zCSdmFtrXI7k4x+/Nc37+5XzGON30IrGu5lM50SHnhjXeG5r2zS8sVdV34needrnnJBjVqE8VguksIxOUF0ktDetcLt1tadqh6/0rvXqmuHWK/039Pp8ovaAAAAIMwQQMAABiECRoAAMAgTNAAAAAGYYIGAAAwiOcZn5OSgpppR8vLp95YLpxotKakubS97jm59LnsJhM1Q0Dj8NP+LqybuJj6Ld1g0UaflMbW7dc9J29FEY13/Nxt3WeDFMfl08eLb7OU4gBTslqqC8362+5XiGlgaVfdxNgNvuqs7p8UONy9j9O9GhPNQkJaTHgs+qWPJqVJ5jp6um3eInrzd3Sfq0Wa7LTbnd9nmqbj5+/rfl0Xfue5as27tnSquimITTEpO2gGPLe2Fy+tZr9YExv1slPXpikkdqexf9ZBAQAAsC0TNAAAgEGYoAEAAAzCBA0AAGAQJmgAAACDMEEDAAAYxPkx+1WeZYw1jhmYxd9TrGyIxuxmk3cjszv94vjDiXwG8dZtKV612S9Jka1lBHTolOJV4zia20xtNynqtdpejI7txcpWR7ZpvHW6txqRuvM+fK4bxOzH+//SS2mkvPauQZb7WNOxNWP248e9wVedh+Ppgaaa0V22Y3+on7fr2qttx1AbDsfTJyxFWO8P9UnunJPU51kJF2asX10XjtnfpOYlV6xRx24d2qI0N7d589RYLiq0pbj/blvaYVVrpmmajsVJSXXoGGtUmLdU0nTg/K0BAACwBRM0AACAQZigAQAADMIEDQAAYBAmaAAAAIMwQQMAABjE+TH7VdRoN9666ncT4irT9prRq+vNFeOtkxgl2syAvmK8dTtBvXn5dJPO05E97s+PSq0isV8nxbLuQi72TYjvraJjp6mOnD08t3jrLSLuL6kbAf20r9vSMXeiqtP27u9649giH7ph95AKSnOjqTZv8FXnze70DtP9neKh877qtrTMxu3uUPeLkf+n29JHkyLL03If5b6a5+pZufYyOp39pZo3ksa70u5xgxeUrjSU9I7VGUt65UyPqg2WHojLFTV21a2x5XtPWoWstScAAAAuzgQNAABgECZoAAAAgzBBAwAAGIQJGgAAwCBM0AAAAAZxfsz+CDaIjp2f6rjg3PGyMftrisvuunSUftKMUO1G6afI1ktH1XYi5bfSjaqvo17r8ed9vdtR1WuI9p7S0grN2zh+qpe+V68dwd2sQ3Ojjh43eLLNV46Hfng6/yC6dWhfLL8xTdN0DEvopOVFbotlAqZpmp6KpTtSpH86Hzch7r9areQiy30MLtev0PGKX90/m0+hUS/XsHxFu2Y0VzmZ0+cdLI+NTs0PtfsemI7tcV/Xjaq0hRWOYmx/Z4mmuHrO2VsDAABgEyZoAAAAgzBBAwAAGIQJGgAAwCBM0AAAAAZhggYAADCIeT0jwnie51+bpul/bDcc4D3wh9Z1/WK3szoEXIA6BLxtZR06a4IGAADAdvwvjgAAAIMwQQMAABiECRoAAMAgTNAAAAAGYYIGAAAwCBM0AACAQZigAQAADMIEDQAAYBAmaAAAAIP4v81ZT4wjauKOAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 864x360 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "approx = outer_product_np(X1, X2, X3)\n",
    "print(\"Relative reconstruction loss:\", np.linalg.norm(np.ndarray.flatten(X_np-approx), 2)  / np.linalg.norm(np.ndarray.flatten(X_np), 2))\n",
    "display_tensor(approx)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Neural NCPD, Rank 4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Relative reconstruction loss: 0.3136097821832106\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA2gAAAEkCAYAAABaExIDAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAoZklEQVR4nO3dz691WVoX8LXPqXPr7SoKsIVqEjAm/pgw8I8wcaIDTTTEaAzIwJhoHDhygBpipFUICSPDDMcmRowTB6g0CcT/AgdAC2KaLqrrfd97z9kOijiQ+3xvn+fufd51qz6fYe1ea6+9z17P2euefr9rWdd1AAAA8O4d3vUAAAAA+JwFGgAAwCQs0AAAACZhgQYAADAJCzQAAIBJWKABAABM4r1r/sff99Xj+kM/cnr8YEjrv4zlmtOMMcZYUoehv9QujeO3v/kD5bHuRgTXX/UzzhUarmEgaYydsZy+fX+7kz3VZ9sunRZ6z3m0x33e0Gf3fzDenr/Tvsnf99X31o9/uKhDwSX8Paq6LcdxCW26dagex+9886vlsZegXYea7SqnP3y4vtFstp6rL6b+FjYe42f3fzDePvTr0Pd+9b31a0Ud6taG1O5W/T3VZ9e5qHvffOE1L/lC1KEvss5jvnEdTe9DVy3QfuhHTuMX/9OfevTYZa1fOl6v9cvUubjaV0v9on8O57pbzq1x/NOf+cnyWOcFYYz80tE6V3rBqS97XMK77CHUj3Cbywfxa//1t+s2S7iA5n5863vHVp/LJXyppT63lq473a9kh/u8pd/4zX/3rPYf//Bp/MIv/9mr2316eb88VtWhDw9vyjb3a10+T0s9sV5f7spjX//Zv10eu2UdCuvSfK7QLtah8HedNU3H4to+/vX/HRrdWHc+npsfwtbj6Pa5tY3H+Bu/+Uv9sYwxvvbDp/Hzv/znHj3WrQ3Vu036o06nv+f0mRzDG+q3L68e/e//+uf+VtmmVbtGv1YmnT8ife1Xf2/7gSSHcOHhnafdZ+XW7yDda+uMZeNrS+9D/i+OAAAAk7BAAwAAmIQFGgAAwCQs0AAAACZxVUjIui7jdfqX3oUYElL8I9bjof7Hdm/DvxpP//A1jSP+Q/Qg/mPUzr/xTv/G8N3nOzwpBmzc+h+ph2MpAGV5CIkrDWv4R7YxrORYD3IJAQLd893O88dQhXqkfxRftRmjDjlK/9g/1aEkjSPdmvgP5jeuQ91/nJ/adfts3ZOtAzbG2CfQJ53u7vrv2q4lXNvavbb0598dPp6rPfMzW8dS1ptUGw7h4lO9qfurb3TqL9WhNMb0jjVCuEg1llgXmsdumV0zRii/3ZpxCRPkEO5/+n5PfaaxnGd4Z3hCql/v16FcLVvXtfC+5hc0AACASVigAQAATMICDQAAYBIWaAAAAJOwQAMAAJiEBRoAAMAkrsp0XZZ1vDrcP3qsiqkeY4zLeLzNGGOci3jPU4hrTbrtYlx+ktrdMuq1OY543Z3xpwjYpBvBv0e8dYi333oca8wlD9sEpO0M0vk6H8/m9/95n9kyxjiNx7dCOC31FgkpHvpcZDa/WuralaKo78I4XoI96mHqM56usytBN1K6u93HDpa39bN3S+3ZOvt9fuaWI4exxvpQSW2Oy/U53ceQN5/6S1uS5PpVH0s1sXOv2FD3e3zruZr6S++PqV04tnkd3fpehTrkFzQAAIBJWKABAABMwgINAABgEhZoAAAAk7BAAwAAmIQFGgAAwCSuitlf1zHu1+vjvd+GNlX0dRWjPcYY5xD8ewjRsWkcMd22E/P85MFGm5R3nNqldNVecmntsMfeAs0+m3Goy/31UcfJGu7JEiJWU9z/cq7H2D3f7TxvDMtYWzH296PeguNQ1aEdtvtI9Stug9CNsG9Mn0ba9x817B3rXndLMzY6ztX3to+HHnen72JUG7nxViYtW4/x2Ze1ljH2d+H9JUbpF5eYtgiJNWqtX/GOoYZ2616K/K/6vIRXyvR+0n5X6krTuPNTR3pXau9z0jT7lhjPcbpqmXN74db7BQ0AAGASFmgAAACTsEADAACYhAUaAADAJCzQAAAAJrFZvMm5FWNTS6lFu+imIKb0tG3DAJ+Zfbdxp7cM9knJgzEJqTfIlJ7YEh6SNcZU1db3rk9THaOZNjVZitOyrGUi2Dmlna0pWe3xa0ypZOdRT/DULh079MLTcopj4+OLIWIvIU12h+TBjb/iPpfG+fZ+hxPe0OzJcM8cw7KMcSqSENP7UNUmOYZa0+lvjP4YL7FdXcCq9MpU8+I7VKoZKRkyhWg26+jmteElpKp23Xr8W9fR7virzzR81n5BAwAAmIQFGgAAwCQs0AAAACZhgQYAADAJCzQAAIBJWKABAABM4qqY/RQrm6KjU2R+FVWd4lrHWg+7inIdY4y7EeJou8mfKUl46zTRPdJJt77uFInflSJnN97KYIwxlvO2na5hK4B0rhilf+mNcel8PpNF/p7Xw/jk8pXHj4UH+vXlru6zaPd2qT+D+1CHUv1K7S7djU82TjSPyd0TJKT/P9VY9ngu454qzWPpw9l6u49bmz1m/5nPyGVdxieXV48eS7XmfvRqSqVba1LMfhpjclxP5bFPzo/X7BSJn+Lru7H3MRK/+8h2HqU93pWS7nzszJM95n63xt5w26Tel23d3wv/BgAAAPjisEADAACYhAUaAADAJCzQAAAAJmGBBgAAMAkLNAAAgElclem6rmPcF5moKbL1bchRrSL4U2x/itI+NMcRo1eTlPzZ6W+H2HiutHkE9PaRsykuP8X6bz2OVrtn3t5ljHEoJ0o9kdMWHFUBiG1GHW99F3Lq07Yjh/twuiR95J373UwgT9HX6VgcY6fP5jYUSZxzKcp5hkj5MW4bs/15p812VXdzbfcxRv2ekqLv07vNKNql96v8rlTPg0OobXGMQV2X6/PFmtesJ5cQ3R9LereOVu0membbW/p06le3/t46nn/jWPyt+QUNAABgEhZoAAAAk7BAAwAAmIQFGgAAwCQs0AAAACZhgQYAADCJq2L2x1hi3GslxUqX/YUkyxg5G7JX0zhi9GqSEjcbfabo2F3Cmrud3jI5+tYRyp3z7REBe0u3jOB+5iVfxjJer6er272+1G3O1XYfaz2JO9uHPDWO6JaPSidS+ql2e0TwV44hNjptQxGj9MP5OltbPCWM5Utpllr5R9axlDUg1afW+1DQ7S9tV5T6TFLM/v165evmE9a0xdEsU6f7zO7wrMco/aRql+pTs8a262i3Vk5WU/5/fkEDAACYhAUaAADAJCzQAAAAJmGBBgAAMAkLNAAAgElYoAEAAEziqtzTZVnHq8P9o8cuIc71Mh5vM8YY5yLm8rQ8XDO0Z7drpNt+LkW9dsbRjYdNJ0sJ6um6O2nzX7mrj4VI06V54anPruMfNJ6hG8fprimONvbZuM+TR9F+kexSh7pbiHQ061DSitO+hItOHXYLcOqz63jLv5/uMP7u1h18qcSa19yaI07H7qPe6bO7Rc3WW9uMkWviIXwI56JdGmO3xnbraOqzu71A5YZ1zS9oAAAAk7BAAwAAmIQFGgAAwCQs0AAAACZhgQYAADAJCzQAAIBJXBWzP9Ycp99xrvrbIfX3HDo99NL5exHQzf5yw96xGFXb6HP57G19rhQNf+nFk8Y+m2KE/Qs+17s43x4O4zI+PLzZtM9L8beqD5b6PKelju89jXN57HBo5t7vEFO/uZRAvMN2H2XViFHON454b0f3V0V2hw90j1jvrW09xmc+BqkOHdd6jqfacFmuf786jPpcqb/0PpTGmBzDnh5Vje3aJUo/nrDRZpa58xyd+pUi/ZN0v7p9VtsEzCLcXr+gAQAATMICDQAAYBIWaAAAAJOwQAMAAJiEBRoAAMAkLNAAAAAmcVXM/rKs47Q8nkd/Dmu9Q4icvStSNY8he/IcYmVTu3QsRqjuEX3fOddMia2dsTSj9Nt97hDrX+rG6XYjuLvna8VRzxXBvY5lvF5PV7dLbartQ1Ik/v1al89zUSefGkd7jm/8MYS07PbWInHbke6x8mQ3fi7T+dKxyebW1Ca7H5dxGJ9e3n/0WLc2VFH05XZEY5TvZE+NI0ljTI7h3ez1pah7e7xfpT73eNer2n1Z5/ch/PbT3QLli3y/Cn5BAwAAmIQFGgAAwCQs0AAAACZhgQYAADAJCzQAAIBJWKABAABM4qoM1vN6GN+6fHD1Saoo2jHqeOsPDm/KNik6NkXOpnjrkGLbjl6NfVanasZNzxJvvZ56sb5t3T8xhBjx5RIOvgTdWO9Of+/AOpYYO72lKvZ6jDHOoTAcQ7uq5o0xepHyT9l6R4kdxpj6jKnY1cHu9hW3tvU4J5urm3pBEdypNqQNQm5V18bojzHVxDRZy62YdnjnidOqO+VeSEmpLGFroc0fuz2i9PfYkmjy9yG/oAEAAEzCAg0AAGASFmgAAACTsEADAACYhAUaAADAJK6K3Putb/7A+Of/6u9efZJLOMvx7eP/vZsqk9ql1J9//89+tjyWMv1OIdDlk8vx0f/+ajmXbX7tsz8Tzra97z9+pzyWEqWqBKhf+st/sT7ZofehLg/1/VoP4QM4Pn7/xxhjnOs+W/29F46lpKBzMzEyXfd9nWQ63r+7vk2SEo2qMT4zBel3fudPjn/xMz9+dbsQ/joO90WbHRJcUx36+Z/6t+WxlP54WOrn6DtFiu6rpbjoMcavf/rny2PJJSTDHUIM2kfH163z3a+Pz7tf+e9/oWyzHrevQ+mZXkNtWML8X09FuzT3u3Mr1aF0vq5qnLdM33zmZaU61E6FLi4/TO/YX6xDoc/YrvHVmfzcT/1ieextMb/HGON+1MeOodbkhPC6z/Q+dFfclH/zq3+tbBOf9fSulOZjeld6P2VzBlVtSLUmvV+Fdrke1l/gy9v6u+zmyeLXSq9QtxsFAAAAiQUaAADAJCzQAAAAJmGBBgAAMAkLNAAAgElYoAEAAEziqvzJZdTRrCmOOsVKx/jYQvdcIXl13Iesy3PM2q4v4L6IZU1x029DFvgxnOvcXGunGNsU612e75kR6o/pxmJHjajXcKv68dCp3Q73cvpxfJfiPC/EJlW8dWrSvS3NR6Xa2mKM/Je21K7SjcvfQ2f8e9SM2GdzjsQ+L40v2+51x61AmrUh9VltPZDi/rs1asf6VdahMNQ0e6r+Ur2Ldaj5PpSOdWrvGM+olzeUovS5wh71sNsu1ZRqnFXtTW3G6L0HhjaeRgAAgElYoAEAAEzCAg0AAGASFmgAAACTsEADAACYhAUaAADAJK6K2R/r6MVRP9Xnlrox+zFSPsVq1ofeFuvfFFN9n2Lv01YAzbX2fYj1T9sLXDox+3vEHXf7PNT3a+1G5s9iwlj8rbUimzsJuM1b2Y2iTvM/zdWkio5+O+pzpbjpc2sUuV267qRsV8W4j9GfH7feiiLFQ3fGscf4D6ldo8+47UBzUsUxPs8U0fExt/9moxhj9O5H2urnPtSo+O4StiQ6xXesjd8Dw3vGLs9lqnu37HOPrZGSze9l83Pr1KhQe/2CBgAAMAkLNAAAgElYoAEAAEzCAg0AAGASFmgAAACTsEADAACYxFW5zesyRicNOUWvXqoRdJOQw5JzCanFKfo+OYY87WMRa3wM5zqGeNhDyu5eG5HMT5xvLOFmVuc79aLAo2Y89BqOLaHP5f7h8f52SLDdRWebgJcUzb+M7eOjq/66dSjtNhHavV5P5bEcK133mvqsvCkL8z5Oy/VjHCPE7KeI+m4UfbJHdP8t++u2633tjHEs7ldna4Gn7NDlGGOfOlSI8fWp1uywY0x6x8oNH//P3ZqXtuY4hJe9Y3hXytsOhQj+5fF3hnEJD197D5dm/epG31dzsrttxx7bfWzdZ+ov1ZNOHQ1N/IIGAAAwCQs0AACASVigAQAATMICDQAAYBIWaAAAAJOwQAMAAJjEVVnKy6ij6ruJoYfz4//9kiLNd4iV/eiQolfrdqcwlip69VWIgP3g8KbuMGne/3S+S1i/V5Gzy2d1f+t7vZz65aF4SMYY6yFE6R/D+c51ny23jpX9Mlt7daizBUc33rrrB4+flMdSrPRpqZ/nD5e3j/73V8t92eZPnD4tj6VtQlLcf2r3PcfX5bFL+BDORY2KNSPFTaetOZp1aBzC+VIMd9Uu1Yw96knqM133ZePtPrp1NI3xOdZRRmTHKPp0iVWT7u4JzW2HYm1LO1iknYCKsXwcat7bUde8VA9TrUnvPKnPqtaMMcbduOH7RHqeH3bYU6Kad2mM6f0qtUv1MLzPxdqctn6q7mX3Nnbqb/qom8MAAABgYxZoAAAAk7BAAwAAmIQFGgAAwCQs0AAAACZhgQYAADCJq2L219GP0y/73HiJmMaXImDfhGNvwyAvIY/zTRHZegwDeb3elceO4VwpAja1e72eymPJ/fr4o9ON0k+6sdhRiL5e2xmrk/gSRPeX87wZi19O8dQmnSs8smt4vD4K0ff3IUr/FJ7Z4+HxevMq9PfRoY6938OHze1FqljsVrTyE9Y9otxDbVvePP4sxGvbQ4q3TxMkRW1X3xMpZnu2uraM+vLT/E/vKNV2H83Y/lbNe+JYiudPT0p13R8dHt8GZIwx7sNA3i7pnSfE7BfbH42Ra2za7uNU3ZTulhipnsR2O7wrpfevchy3rVGxJnbqV6pDqb90/8v+QnfX9wYAAMAeLNAAAAAmYYEGAAAwCQs0AACASVigAQAATMICDQAAYBJXZWEua4hYbcbbl/2lCNj6UDxXaniOWbXbOm+9V8FT52uuwy+dPRBSBOlsMcmF5fL4NRSJ3n90MD14NzbTWLjKq5BhfRei9E9hal0uj0dHp5j9V4c67j85h5pxDNf2KmwvkFRbiCz3dZR2dyuQ5aG+X2uKxT6G84Uo+tY4U429dW1O111JMdXq2tW62w7dUqpDh/DSlo6lmP0PD2nboTRH6kPpfDe1xxwp3ofiVgCh/rZrVIjSX97W3x/rXWMrqU5c/hib33+/oAEAAEzCAg0AAGASFmgAAACTsEADAACYhAUaAADAJCzQAAAAJnFVzP66PBE1HtqVx4olYifh/alzLaHPU4psDXm0Kd76VMRK36Uo7SXEkwZV3PRT7kLE7TlG3Bbn68aTJt0I6BQD27tdtW50LFdbl1Af0sfQOBZ3xGieK9ahlFwcHqM066o46tymN0EOobYlp1D3Lo2/I67HZh0K8zhG6XdrVGp3Ke5lJ77+XQhbCIxqC4EXVCvXsUPdKB7bOK2a5+oe615bdSzVvPhFHR6V9M6WKkM3Lj9fw5dQqms7bOnRrveVqvY+ZeNr8wsaAADAJCzQAAAAJmGBBgAAMAkLNAAAgElYoAEAAEziqhTH0ycP4+Nv/O71Z+kkM3VT8ZopKv/gV3+81S6qkmXOISFmj/SYGyZj/dh//kZ57NVyXx779uUrrfN9eHhTHkvJlufwt4lf+PqPPfrfY3pVUwibilLKaUwQrG5Jp80T46g8/O5VZeePOX3yMD7+tUYd6rjxvPp73TqUxlnVoYeQspeuLR1LKa6ptm2cfvVP/st/KI+lOvTJ5VXrfB8ub8tj3WTLf/jT/+jR/34JIY7dtL8Q5htdwlROta15S+pxpHtSjOP5deh+fPyN/1WctJlil5JCK5ft69D6lbvy2PJZ/azHPk+P3++//5d+IjRqvutVKaFj5Lq3sb/yH/9HeSy9nxzC5EntXh1CbTv3atuPvvqtR//7dy7vl22+evzD8tj/OX9Peeyjw2flsdfrqTz2SXh//Jdf/zvlsepdqVufUs2ratT979d1yC9oAAAAk7BAAwAAmIQFGgAAwCQs0AAAACZhgQYAADAJCzQAAIBJXJczu4w6sjlFvXYilFPcbDpXt103MnvjeOibx+VvPf4XIkXVljH1e9yq9EimR6E7lka77rmqdrtu+tDdpiLFw3fO1Z3Hd3WUcNQ5X4qiTluBdB3D+TaWIqDTnyU/De2OIXv5GJ7qQ6g1l7BPRWcLi06bMUb+U22a/80a1akpW29Jskkd2noLocoOUfrTqN4pX9L5is/0zaWu52n7jWN4ruK2HeFQiqlPXhfXkPrrHjutD70+w31u1dHmV9XWr+h+QQMAAJiEBRoAAMAkLNAAAAAmYYEGAAAwCQs0AACASVigAQAATOK6mP111PHLW+dLxozebjZ5iudvrlVTrH8Vp9uN4O5G4qd2G8fsn5Zz69jdUser5vPV7VL0ddKOqq5sHFP/ZJ9b79SQ2ryLXRqWUT+33Xnc0T1XmnP3YR5052oVp/9Qz8f0oCwhgn8N0f1LON+aoq8bNTHVhbuR6lB9LMXld+vQOUVmfxmlR3y6pPnldvWmu33QS5DqUJLqYdpCpHu+xljSO0+M2Q+15hD2m0h1KI0lqWri/ejV2O4Yz+HFrHttN9XZbWz7UQAAANBhgQYAADAJCzQAAIBJWKABAABMwgINAABgEhZoAAAAk7guZj/FW28c157HsUPcfCcuv6sbKb2HW35uL0SZftu9VSklPKVsd6P006M8Qcz+s5+4dZTzJEW5xy5TLHPlEj68PeZVdwuO6tita02y8dYj33/4rGzyaodI5g8O9+WxU4jMvoTZUNWGNTyqaX5360K6XedTs8/qlqTxpz8nv4t4/nUdS7Etxnq67tVqNmuYj5tXtm4d6tbYPbYyKnx0rOtQkmLjUwT/q1CHuj46vH70vx/Wehzfd3hTHns76gL2vUvd7tWxd23pHauqiYe0A02qo+Fcl0ZJ8AsaAADAJCzQAAAAJmGBBgAAMAkLNAAAgElYoAEAAEzCAg0AAGASLzsLdkvnlHfeVEV3pyjwFAHbjYe9YbsPQ7zqq6WOST03w3vT+Y675Su/eynqdRbVGPf8VNa0XcbWOtH2Y8wVb/8SNKKvTyHvOB/rRfCnKP1jyJtPNWqWOR7H0Yy3r/psbQPyrixjrGnLnC+qVGMvjQ/wC1wr2/VkeXz7hue0uwtjOYffaqp6ebeGrQDCRL4bdbtUm8+hEKVra9evjo3P9SWsLgAAAHOyQAMAAJiEBRoAAMAkLNAAAAAmYYEGAAAwCQs0AACASVwXs7+uY1yKGMyt41AP9dpxCTH17Zjtu1OvXUcVvz9GjvvvxuXf0LfPr8pjbw/1dX9y/krrfJe1fk4OIbI1qRJidwn8TR9papcurRl9XXaX4rIb7XaN0p5kHrSFunfTc1V1fvSjxW8ZSX5J2crhUHe7j64UHX3TOhTE+doczIuK02df3e1KuvaoscV7Z5rfKYr+JejWyj1qbOozvr/cstw3Pm6/oAEAAEzCAg0AAGASFmgAAACTsEADAACYhAUaAADAJCzQAAAAJnFdzP6y1BGlW8fsh3jVNcXUp1jWrcf4Utwwnv8S1vwpEv/c/FtBbPdF/rhfeJr8Xpa0TUUQa8othXj7zSOnm3UhbnMS7mP6bLa+/6ewxcYxFIa7uH9FLUVmn8L5UruyXKaSl4bfbLemjyb1ecv6m74+dhzHcnm88/BV90SHxbx7Ke8uaZuj6p50ry29I6b6tcu9fPx8d0tdK5NUo7rt4rZD4XRVn6dwbanmpXaH2K4ef+ozzsVquqW6lr6Gm/V3wyYAAADswQINAABgEhZoAAAAk7BAAwAAmIQFGgAAwCQs0AAAACZxXcz+OupI1B2iUutxNM+V2jXjuaMqOjrEVMeY7TT+avuDp/pM7Rox3D/6/m+VTT5cHspj37q8Xx5L0bEfHO5b7ZIqRjXFTccE23CLQ8p2jL6N0bHpYyuOpXGkWNkYwV2ViuemxS+jfP7WY7gxnTnSjbZP7dIYu9L5qmtr3qs4q1L09fun1LLVZ1WjfuKn/3HvXGkY3fmYhD7f/NVvPT6O5qkO4QIeLvUFpHbHQ9jO4FC3O18ev4rU5s1DXWzSdgVLdey/1d9H352l3h6i+45y/9wxbePwh59t3me5PccpvIY27+NyDu9Ye9Tfws9//W+Wx+L3YDgW3zVSu+5rc9Fnux7usMtBvCd/4/frY8XFlTVj5J0kHs5pm6nineVX0rYDAAAATMECDQAAYBIWaAAAAJOwQAMAAJiEBRoAAMAkrktxfAm66UldKdJla3skym3c7jshjTH9OeDT1K7pEKJ9LiFmqArwaT9a3YDTZsJjatdNctpyHDecMd+9W9eNSkqTjcmQ4VgnoXbjdNcnpXaNFMd1h2+2dNntFMfgdAxJdIU9vo4OzaTGlPC4LI/fsNSmSn58ahxVItuedWgJc65MfgztUjptOtcec3VND1mnz5Rq/cKlutBNcWynIHdTHBtBpe3k6qY0ltMxvQc+fjNTcm1KeOxI5/ILGgAAwCQs0AAAACZhgQYAADAJCzQAAIBJWKABAABMwgINAABgEtuFEXdjqruxzJ3+0hjPzajXlCdaRVU3o2+XEEfbie59ql1nnD94/LRs8mqpx383GlHgY4wPDg/lsVPIlU1nKz/SFJmbOtyhXeyzMRXDjgQ5j7oxjkkC7bn05twuNt7moBthnbTnSNPDOdTmQopsjue61BM5vSSszc/t4fz4+dbwwaUx5or+eJ971qEUi791u9hmjy0xtjbLFic7aEfiJ914/ma5r64hXloaR28Ysc9Um++LWjNGXW9SlH7a0qOqa8kaLswvaAAAAJOwQAMAAJiEBRoAAMAkLNAAAAAmYYEGAAAwCQs0AACASWwXs39L3VjWW0fOljnjLyRWtjHO3zt/WB774PCmPPat8wdXn2uMMb53fV0eO4Ts1UvaHqG67BDlmtKtU7xtMxU7Z9U2+ky3o32u2R7zQ5jjlxtGTs8Sbz3LOMbIY2lsnXK4744jHArzOM6fJMyRu/eu3/olxUPHdue63XvH+sJTrP97h+vbHUObS7i0UxhjZdfHv/uuUbWrtuwZo79dxixbI72U96GG7vd7iueP233ETnvNqvPFd54d6mG3Nh9DHW3F7Ldf2h63hIv2CxoAAMAkLNAAAAAmYYEGAAAwCQs0AACASVigAQAATMICDQAAYBLbxeynGNitY1S7kczpWDeqtnu+W9ojxra4tu8/fFY2eX8JsdHH77SG8eHyUB47pMzWST6aqBuL24i+b0eIv6SE5JTTfcuo51tv9zHDuW6sG0WdnvV2vHXTQ2Prh/SRpl0mqrjpMcY4p3GE+3UOfV6qY5e6wzzGul0Zmb3jlF9CrcnP2OPt1lDU07n2qDXx+6DT5/HYGseLMMk2OmPstKVPda70PrHHViZhjKk2dGL2R9gKJNWoapqm2+sXNAAAgElYoAEAAEzCAg0AAGASFmgAAACTsEADAACYhAUaAADAJLaL2d8jpr7SjY5N7VLUa8on7mhuSbAee+vp9b1mjG0aZ+H12nukXq+nVrtjCCk9hczWMuY5SY9Wt12KnH1qPBv2uXaTjr+4Se1P2yOa/5YR/LfcWuApG4/lcrdpd2OMMc5hjsQI6CClOad46MohdJi+oR/CuVJFX5uf28P58fMdD3V/+X7UV1dNm3f29M8y7/aoNZ0+Z7kft9Yt5833iVgAgvLVLO3ukLaSSO2+qxE90i7F84dtQuqY/XS2+uLSliR1BH/dxi9oAAAAk7BAAwAAmIQFGgAAwCQs0AAAACZhgQYAADAJCzQAAIBJXJmJvm4fiXrL/sKx5e193S5kbsbo+1PR7nxunWt5qNulKP1uu05k7rcuH5RNPji8KY9961y3O6a4/KW+/4cYs1+3K2NgU0xteoxT7H0z+jadL/WZIm5vZddk/pSPm7bLCPG4rXN123XH3+2zozvGbp+N8S+hxHZ22BjjibnanP+p3funhydG9MelmP3Yroi9H2OM0zFF2IdtTg51u/tiC5djc/zvHcMHXnj2tFjqTtbmXC3bbT2Hu+N4zliqdntc2yy6kfhpJ6buO0P3Nhfna79LpHeX9PrefI9KdbTabinV0VSj3i7X71eUaqhf0AAAACZhgQYAADAJCzQAAIBJWKABAABMwgINAABgEhZoAAAAk7guZn9Zxijicccl5FxWbVK7dhR1OFeIjl2PId+zGzlbxdunMZ7qj2RN5wpx/zFKv3ufi8/tp77+k3Wbphi9usOfGNa//vuP/vd+gnh9AQ8h3vp4qNu9FyKsj+HY+VLEW4c2r+/rZzKNsYyq/cb18eF/THqmO6otMdIcSA9fNzo6bW2R4pXjJCmOdbc4iWNs1tHu+cprq5vEeOhUDlO77rEg1YaOXIfq+bSGfOs0/1O7qg6toQ6l+9F67LbYcqSqQ2n7nVhTikGlSPzwztCe43tE31djSVscvXSpDqVmaauc+D3QbJdUX407xOw3hvH5sdDnfaO2pVp5DsfuQ42qzpWmqF/QAAAAJmGBBgAAMAkLNAAAgElYoAEAAEzCAg0AAGASFmgAAACTuC5mfw/dGNitdSOgk62jam8db93RPVVKHm7GynZVcaiXZj7sIQwyRVFfwue9hHbp2KWKlQ1t0hhbqfBbTPlqC4vudh/nor9jiPPvnuuW8zHZI0p7D41Y8nXjXRjGGCNMg5v+qTNFQCfpsSu3xHjifFsf6/aXd8LZ6R1jHWPcP75lyFLVp5G3vanaddp83jDlvPfm/9rddqiSauwe74d7bIFSnqt3qpwp3+sy1q/Urvh44u1I9fDGr/zHMP+rb/BuPUxTo/P+6Bc0AACASVigAQAATMICDQAAYBIWaAAAAJOwQAMAAJjEu09xnEVKQUyO1yeM5fibHeLHbmgJgVK5YTgUwvLaCY8x3Orxg93wvT2SyVLKUDq2lte2fVJb2eee4YEpoStGLN0y6bSZgvYS0h+nGeNtT9dNSEuOh+u/k7qJhaldGkeqGyk9rVOH0jiOh/qLoLy2535myyjnQkw6DLrtSrMkNXbbzZQ02xhLuy6k18ru7e81azVMY9zjE43nS+8ojTaxrqU6VCVoly38ggYAADANCzQAAIBJWKABAABMwgINAABgEhZoAAAAk7BAAwAAmMT1MftV/GrKuUyRrWtjjdiNeY4xqSGeuNtndSzG7AczxVtXY2kOI0bHNo91o1478dYpHjrJ8daN6OiRY2AvxbEY9RoOpjG+E2lupa00unOyI50rjj/0mbb76Jyr2y7d41ki+F+Ic2Prh24q+flSfw+nWrOGE6Z26Xx1m14gdTmO5075dZRzYQnzIL3yVO06r0mfN4z7yVw9jifHkh7AW9bYSTRfC/LHlr4+vqwlNlx3qjVrFX2fbnJ450nnqqZUekT8ggYAADAJCzQAAIBJWKABAABMwgINAABgEhZoAAAAk7BAAwAAmMT1MfspRnnLNl0x5jm068Zbx8zTja97pnjrMl642V+KxA99xsji5u1/OF//d4sYyxrPdSyPVRGwY+StAFK7KgZ2DdGx6X60Hrt3lbac5kEj0px5LQ/pYLPPZh1qa2z30a1D61pfXKo1aXuRtN1HVW/Sth3HQ32TW9t9mPJz6G470t1T4suqu7tTqqVVm+57YBpj8x0x1YZOzH6qa+fY7vH/3ts8BAAAgJuyQAMAAJiEBRoAAMAkLNAAAAAmYYEGAAAwCQs0AACASVwfs1/pxqFW7TptnmrH7i7dJyp93CmydYc/MaSY+s3P1WyX4u3r4P663XLDa35n0k3bekuMl6Abb93V3QKlYU2ToNvnHnnI4ZbcP1x/ESn2PjmHi0v1MMVRX6pc6VFv3ZGi9N+G+3EMW2hUY3xnU75z4u5gb/0+1Hk3C22WUDN22dpia91nLLVLt7h5ulQ2yne6NMb02dz4fS7VjVbMfth2JG1JVJ1rDZ/aS3jEAQAAvhQs0AAAACZhgQYAADAJCzQAAIBJWKABAABMwgINAABgEst6RXzrsiy/N8b4n/sNB/gS+NPruv5gt7E6BGxAHQLetbIOXbVAAwAAYD/+L44AAACTsEADAACYhAUaAADAJCzQAAAAJmGBBgAAMAkLNAAAgElYoAEAAEzCAg0AAGASFmgAAACT+L/H9WXqVBi88AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 864x360 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "approx = outer_product_np(A_A1 @ A_S1, B_A1 @ B_S1, C_A1 @ C_S1)\n",
    "print(\"Relative reconstruction loss:\", np.linalg.norm(np.ndarray.flatten(X_np-approx), 2)  / np.linalg.norm(np.ndarray.flatten(X_np), 2))\n",
    "display_tensor(approx)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Neural NCPD, Rank 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Relative reconstruction loss: 0.5011680801673786\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA2gAAAEkCAYAAABaExIDAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAfu0lEQVR4nO3dzY4tWXoW4BWZVd0lG6Me0EJILcHQvjObKdeAGdpiwBAhRtyJLwHhESAhdaubnzaIxl0ndzCoBlwm15u1v7NW5LfzPM/wRK+fiB2xdny5u951nOc5AAAAeH9P7z0BAAAAvqNAAwAAaEKBBgAA0IQCDQAAoAkFGgAAQBMKNAAAgCa+uud//Px3f//8+qc/WTqBWcr/cdzf5rPa/XWTWvXBdz348X/59r2n8Plmn0G4t/jhfvPtr8dvX35TvprVdejKHUUeYh3acD2O0Od54fPTax1KJ/7gC/4D++x16A9+//yqsA5V14YrXTrHl7DmVceqPnKL233zi0+XjfVmu2R5nztu5OIkV9/Mqb+C33z76/HbT//z1U7vKtC+/ulPxs/+2R/fP4PwzXy7vf5wPj3d5t2F/o7UbjLWGGN89Ze/N2/3NP8Q00vHcXv9YHqJOcLzHKUXo+cwXvGFatbuH/2bn88bVVUr8qLj08vdY51P6WKF+/Ul3K/P6Ytr7QJ43NYvqLNr8hf/8V9/Vr9f//Qn42d/+if3zyeuQ68fe0rPfhjrKTxYtzCPr9M6lG71MN5sHYrfaWkdimve/Fhah8ovJBOt1qHiC8JsbYhrTdJpHZrMJc2jatc69NVPfzJ+9qf/+PUxQ7u0NszWqMWPx5t9Vtev+Ao9aXf7rz8KjUKHs3VtjDGe04tN6DOdQLo1J4/IH/7zX4WxLv4VIlnd5y1crOo7W2qXxnsOXzyp3czT2j+i/sW//1fzoZaOBAAAQJkCDQAAoAkFGgAAQBMKNAAAgCbuCgkZY5SiuNJ/nL+yzXcNQ5BDrcfY5xF6nf53thenp8XxVv/XxxtiqFKAxfnof2JI16t4LfP1mgVGrP+Pi3cEj4zx3VTTf6he6/Suf/7dweJak+ae/tvw1GU6WliHYnhQmEd5rVl9q2wIDyoHgfBhld5tVvcXv/xDuxCyEb9X07qXAlBmeQzFIJC4RlWzZuKH8wDxm13sCEDZofm6/eivtwAAAB+GAg0AAKAJBRoAAEATCjQAAIAmFGgAAABNKNAAAACauDtm/3gqxImGWM3z9nqNeDyFnNQUex/mV+xynOmcc/b13W2OYjxsisXNkbnhWCUN/Gl9zR9jhHfEpE7i4acR9Z+jGt2dnqk0z1mfxeu45Zq84TjGeCpES6cI69mh+HiHORwpbjr0GZ+5NJnFz0gea8OxpNDu+PQyPZbu2erWEFueg8kznr4j3uN5nOmyPcqu7T7GmD/n8Wu1snal/u7u7XfCe01ljm+P9/o/53sh3ENpihveeWbvBbFdp0j56lxm7apR+lWhz7zWFLYySvN/fp4fW8wvaAAAAE0o0AAAAJpQoAEAADShQAMAAGhCgQYAANCEAg0AAKCJu2P2p+mTMZc5dThplyKxYyZ+imUN8copFTSldIZmq8eKin2mucQ02tnBW8iArsbK7ohsLYwn3vq+sXZavd3HMflbVdru43iU7T5mT3Kcf+ovjHTxdh/TU3uuncD5XIyO3hGZPemz01qTzjt+3tN2jbY5eMNxzNeH6tow/e5JC0MxEn/1HN802VIpb1dS/Fx3bOkR17ZCFH1VuteL2zRFlXO4cD0co7i10A7xs7n/ufELGgAAQBMKNAAAgCYUaAAAAE0o0AAAAJpQoAEAADRxd4pjJa0xpSfOAoFimzBWSiYq582ltKPQ6zScqJriWA1yqiY8hmZTTxfX/Fcm9DyILilvO+dRChKLEYkfOE12Ml55zSikKr51bHWa7PGyIc7s6hTHR0iTrSb6LU7mzcm1G6/J5Nmqrg2V/qpnt3qOb443+feYUFnpcJe4fs2+PzasGSkNcEfi9azPNP/qPIrXJD7/z6HhvAApzaOS1Jj4BQ0AAKAJBRoAAEATCjQAAIAmFGgAAABNKNAAAACaUKABAAA0cX/MfsxmnjRJxyYl4hHGiUGcqV2KzI6TDBGesd1d//z2wUJK+FvHYmJ2JfH0tj7e+t0ilB9Ul+uV5vHZbPfx/fGu3O4j6LLdx/nc6G+PxTjq2VYBH2LNm16Tq7PT685zvj6U14ZCf6tj78fIc7xVI/+nUfSpTeqw2O7KPj/CNkCrz2HDNVm+Jjb53Bp9iwEAAHzZFGgAAABNKNAAAACaUKABAAA0oUADAABoQoEGAADQxP0x+6vNoldTtH2KeS1sA/DGcNUU2/UaJRDPrtfWaPWF2swzxGzHY1fOo0nk7P9xnjnqudbpXf/8u4PFKP24J8b8UNwSI+7BURgr3ZZhHuV465j5fX+7WUT9LjHmOUXpp3Vo0u4Ip1aOm96wDuXtPtaOldaonWv9tOfi2lDqr/jOs2P9SnM5Z/dtdfuNW/GalGP2C8/4ld/hb6l+x68+hw3XJK81hfGafG5+QQMAAGhCgQYAANCEAg0AAKAJBRoAAEATCjQAAIAmFGgAAABN3B+zX4i3TrH483TSWgTsUY6OTcdChG9oOE163RA3HZN2qzG2odm0u2rMc7AjJjnN83jpEbHaRscI/tI6lI5N+iuntYd1KPVZTRlPB2dR9NWxyjHVtT5j4vHs++N5w98eNzwH53NYhz69lPr8sJpt9zHGGGeIep+7/30oPjvhu6zy7Iwxxij3mdpNjhWj9OP8Q7vqlkrlNeqjujKa/2pN1hq/oAEAADShQAMAAGhCgQYAANCEAg0AAKAJBRoAAEATCjQAAIAm7o7ZP54K8ZkhcvOY1IjH0zyLOkXpp/mFLscZStUznXMl3zrOP/UXRorzTw3DsXRus3bPz6FRTb7+G+JQJ7H+O7YQiPOvRnenz3vWZ3msd4qjLWS9H4V9KspJwml+1Xs2zaVwG1Uf/bJiYjZ0NX22imvDbI2qPh9x/YoN75/j25OZfK9W17XCWGOMcRb3K8rt0mSgxi9oAAAATSjQAAAAmlCgAQAANKFAAwAAaEKBBgAA0IQCDQAAoIm7Y/bP2+I80VmOashXPWMu6zwm9RbaPafM2eKxo3Ktqnm61TmGY6Wk4FvYy2CHlGdczUifHDviNg0X5+zGHOcN7SaOyZYEY+y7JscxxlMhZj+2mMz1SHHTMdo+rEPpz2Jpu4x0ORfH+setGtI0iluZxHWokIp9vGxYh6przaMr7zWxeKyGZs95dW2YHgn9pbUwvitV+5w3yyZrbF5r0rmlsVKXxW2Twndd6aeOK5+rR9Hl+W9y/f2CBgAA0IQCDQAAoAkFGgAAQBMKNAAAgCYUaAAAAE0o0AAAAJq4O2Y/5iHPxHzo6+yYRY6+nvxzjwTPj6FLLCuXOc834qNn7dKxWSr7c2pUvPfiNiHzQ2m0Mx2d9VndWiQpbulRPe/KPKvbP1y9zcbx0uSL4srI6QfbymA62+LaUOmvvKFE2sooNavsezHCR5tOIE6kuo1O6DNJ46UI/ml/TZ7vTh7s+d/NL2gAAABNKNAAAACaUKABAAA0oUADAABoQoEGAADQhAINAACgibti9s9zjNttbU03i7e+pUjNGHc6n1+Mh03pnrcUK5vaVcaaH4vptimqNudzz5ulj3rW7mlDzb8jejW2e6n1yfcclejhH6jUcyX6enVc9ngjXTkduzBlOEXil3dNuTjWf6Z8X4YPLq6/XTx6rHfHCO7JwxAvdfw+vr+/o7gwlOcYO629Y1X6i2tUNZ6/vBAt1vFe/9sefT2pqn42s3ahO7+gAQAANKFAAwAAaEKBBgAA0IQCDQAAoAkFGgAAQBN3pTiOl2O8/Lcf3T3IEVMXX//ns1g6niHaJ6Uxfv0/Up9hwGqy4qzNp/Vjnc+hzytT456KHabrWO0zKaQk5WS4WtrR8fII0XDv4NPTuP3ym6Vdzp6fFFobU8SeauvQj/5q3md6HksJryk9qhpkmq7JjnVoFoz1VRrsYpVkrzHG8e3kyyC1edmQULlhHSrPc6aSKPe5U/j0NG6/+vGrh9LacAtrw+xBjv1Vg4yLc0zrV2Ud+tFfzRfZNMe41lTX7eo71my8dF9WUxBfNqRM7+hztXC9zup74GxtTmt2cT2v8AsaAABAEwo0AACAJhRoAAAATSjQAAAAmlCgAQAANKFAAwAAaOK+mP1xxIjVkll/MQt1fuiIEZhhHqnPWrPpwXRq1VjZLceSSrsYRX+xNJdZVGohfv9NKZZ1x3gVjzDHH6DybNWfx/XrUNxmY8czvlpx3esy/fJzUI1evnIdegQN16HpfVtcGyr9lc+8un4Fy9eh6jvbhne9ZDZe2n6nHA1flLcCmrt6nhX5Oi8ebHGUfuIXNAAAgCYUaAAAAE0o0AAAAJpQoAEAADShQAMAAGhCgQYAANDEfTH751gelTqLJz1vhXE2zOONZjGWNUa9Lh4szn/Hec8OpkjWahTyjnjl53m78/m6v1sc4T7vE2/bZR7/T+nZenTVhejKsbp8LukZ3hGJv2Nte3r9HHasC53iwGeWz7HHaX1fl+cnWbwOPcxaXjjv86vnLVOpWB43f7Wwxp5P4cNJn0FlvV/9/RG6e/SPDAAA4MNQoAEAADShQAMAAGhCgQYAANCEAg0AAKAJBRoAAEAT98XsH2N5SXfOsktTOumGCOgzJWemc65E9ab5V69vcY7pvGO72TlU46Zht8p2H9XU9dSuSax0Od66yfzHCOcQItnLJ7Bju4/CeGlrjtVjjTHG8XLheRcfuC3X5IeopGqv3n6nqrj9zqV9VrcI6rL+3oo35oXryZu6vNOl+adjL1fuQbOWX9AAAACaUKABAAA0oUADAABoQoEGAADQhAINAACgCQUaAABAE/fF7I+xPmFylvMeIzWrY9WaxTjXwnipv8ujr1e3+/Qyb1ONa7043vp4ea/M5u+L8daPbMFpnU+FTuKeErNx5sfiuhDmd6T9K6rbXqQuZ0tsahMegWq7OP/QrNSuGimd7Oizy3jVNfbqazJTmceSdej1f6+uDbOHK/ZX/ApcPcequJ6E9Te+IlbX7TSXQrsjbPdxPl37XHWay6W6zz9Mzy9oAAAATSjQAAAAmlCgAQAANKFAAwAAaEKBBgAA0IQCDQAAoIn7Y/ZXp5rPskur46xPXc+psoVj1UjpDc3Ksd5TzynfthrlWsgQf0shVjpG0RaVo28vtHyOn3tax6j9aSl85udsUiFuuh4Nn+690C5FTlc+hmJMdVKeY/XYlQnKaa15kC1EWoz1kUzWh/h1lp6tScPq12Oyeo5vjjeZaI70Tx3OT7zaZ1yjbveP1+U7fIyL57Jja45qn6vncuG2I35BAwAAaEKBBgAA0IQCDQAAoAkFGgAAQBMKNAAAgCYUaAAAAE3cH7NfSZGsRKWGaO+y2/xQipyOsaxpvNkOAikRtHraxfmXle6D4kTSvXDhnxhS3Hyn8VKc7upzuPqajDG+u/cWDzt77irP93cdfsZkCpavKRefW4z8r3S4I75+cYRyebwdkfhd4v53RHCz1pVrQ8i9r76zxT4Xf7Hs+H6sxuVfuqXPjrUy9JnPbfFcLvwe8AsaAABAEwo0AACAJhRoAAAATSjQAAAAmlCgAQAANKFAAwAAaOLOmP2zFFl5hFjTaR5qsXQ8Q/bqETJgzzBenH7hWOwvbAWQxjpSu3Qtq5G5s2MxEl8U8hXeJfr+Ssc5zue169BZWYdSNHxYJy9fh6YdhmPF9TdOI51bdR1afat3idJf3WYMUfSrHWN+T1fXhtlDviPafsP6lecyW2ML63Loboy8juYtJYrtCuvlju1wrvzu3xHbf+lWAEl125HF25X4BQ0AAKAJBRoAAEATCjQAAIAmFGgAAABNKNAAAACauCvF8ZtfvIw/+vP/vHYGs9STC5NSxhjj+PQyHy6lxxTGO17mkYvn87xmjnNM7VIyTmgXTT6Dv/zjvz9vUw3huYWUp0Ka3xhjhLDP8c0vX78mMUWvKKZvBjGlarF0rSrX5Lf/8uv6ZMYY3/zi0/ijP/vVZ/Xx/6msQ1dLc6k8x2Ediv2FdWg8FaMaF1/nf/tP/t7S/sYYb6Rerk9P+/HPJ89JeiCLjrDGJre0/lZTh2dSCmvh+v/1v/jcdejb8Yd//svJhFLiX+EZqb4PVRM/n5/nx27FL63Zeaf3q5RKWH13SetXtc+Jf/cn/6DWMN3P1UTP4jN+fvX6gOkZPsPtE+dfTBY+wkf68nfCwcoLzOIXov/1T+cXyy9oAAAATSjQAAAAmlCgAQAANKFAAwAAaEKBBgAA0IQCDQAAoIm7YvbHMa6LnU7jPEr09Uc2O+8HuRw7IvOXq84xtVuf0P0xfKnPMbBel/WkyzxoLSbHXzcN/ha/oAEAADShQAMAAGhCgQYAANCEAg0AAKAJBRoAAEATCjQAAIAm7ovZP8cY50Whm1eN8w7Op2L07QNE5p4bSv4jBb1uGG8awb/h8levV7Vd5RTikxg6nF3HllsczNabB3jmPoS0Jt4WfxekTOncsHRotR3PT9yZIx0szqVyDjs+tm3S+8uVa8oHfo8q67Kmhxs6PR8xEr+4RFWex9SmOo+q+hpVeD4u3MbIL2gAAABNKNAAAACaUKABAAA0oUADAABoQoEGAADQhAINAACgifti9o8xxvOkplsdhVyNXb4yrrnoCPM4n5tEwBYdt/mxcjx0+tiu/Eg3jFWOjr66XaG/6W4FC+ZwvExutGq89aTdOVvv0hzeGCtus9ElAnqH9Nm0idJv5MJz2LL1xeo46nQ9PspjU9nuIzxX8V2jut1PF4vX+t81LE/n9e52vDRU2y3O5794id0S3X/luc2uf5i8X9AAAACaUKABAAA0oUADAABoQoEGAADQhAINAACgCQUaAABAE/fF7I8xj0OO0aUFIcE6xzUX211oS8x2bHdhXPOOkj+d2pVJwRvGKqflF69J2gZhPlg4VPi8V0R6z56huM1DeO5K16Uapf+lSmvUlduj7Njuo8s6tGOOX0K8fRcfeZuN1Xa8K5Xj+SvzWNvdZ/VZ2Q4gvoOE7R1Cwy13f2W9bPLO6Rc0AACAJhRoAAAATSjQAAAAmlCgAQAANKFAAwAAaEKBBgAA0MR9MfvnWB81WumvSVz+GKM0lyPERp8p7jSNdfU1mYx3vIQ21XjSW4rFXX/esy53XOFSvPvF4iUuzH/JRzaJQ46x/zEWf2GbN9pxpysj+Mv7XoRjq3dAufrrL+2Lkb53PALvznYfje3YEqPLq3FYM+L3f5dr0mRLFb+gAQAANKFAAwAAaEKBBgAA0IQCDQAAoAkFGgAAQBMKNAAAgCbui9k/xvr46Ep/jxJh/SjzrHjwc9uQzs+VZts8xNj1+z/08znEBaftMqp/+noKDdMz9+DPY9ksRjxugVB8+NOWEs/rF5T51g/F/mLcf23+eauJ0qH5WGmXmcrz9oU+Mp/tyrXmwbcJqH8PhD6LEfDlHUSeZvsOhcHCenhWtxaK61dtLqWLUt02pXAr+wUNAACgCQUaAABAEwo0AACAJhRoAAAATSjQAAAAmlCgAQAANHFfzP455hmfKfszxbLO2lXavNVuh8J4Z4qOrUZph2M5Crl4nWdNwh0VY56DGGdejLeOYeyT61WdfznOOX00xXjrGBU+Haw4j1mbz31Ej2OM5+fX+379n7c4UyR+0fHXvw0HF8fs38LNkM7t5aXWLq0n1Ws56fOb//R1rb84Vji04U+d3/xqcmDDV9xRjL6+FZ+30k4H1XjxybHj28IcvtfBUb9vU5/3mqyF26x+x/r209r+Pkda2wp+7+e1+2PHu8bqZzw9w+V3vWKf6dw+/fpHLeYx+2yOb+cd+gUNAACgCQUaAABAEwo0AACAJhRoAAAATSjQAAAAmrgvxTG5Mj3x6qTGj6yQ1Jj7mx8qf2oxcrHWa0wRmwWLlkbK0y83LHa6PD1t8e3zg5zn8rStRxeTYSeOctRWo/V3MpdbSnctp8nOj5VTVWMy5KRhGmtx0uFbyqmpF6Y4Lm3zN1mHlkgpzUllzRtjjONl/iCfz/OHqzLPaspp2Y5nfHZJqqm2qV2aSHmNKqSRF+dxFpMyZ/yCBgAA0IQCDQAAoAkFGgAAQBMKNAAAgCYUaAAAAE0o0AAAAJq4M2b/nMeyr45rj9NIEcopA3PerhrZ+vCqkdmTaxkjTYtx00dq+FQbL91C04jYDfHWMXI2zb8a31v5uFN07Lv8iecY42nxwLMb4vJIebHdK6TtJOrbXlQf5OqAi/tLS+WGPpe3q4617REO61D1HWW16nvZlXO8fdw1r/pcVSPxL33Gr9zqpz7cpVsBxe+dwmfqFzQAAIAmFGgAAABNKNAAAACaUKABAAA0oUADAABoQoEGAADQxJ0x+0XF6Psr53G8zLPEP3QE/+rrH7Phq31uOBbMolLPFDdfVdyVIEbfp4aFc4jRsRc+vn9j1OsGfp8TvNtxK2wh0iUKvJHL4+aD6TpU7rDaMNiwvci0u2qE9eKY8O91UFkfHmFNeYQ5PoLq85GarX4Odtixjl74lVTepiV8AEfhovgFDQAAoAkFGgAAQBMKNAAAgCYUaAAAAE0o0AAAAJpQoAEAADSxLmb/yljmLzQC+hFUI2BjvHU8Vuw05tEudnV094dPSD7ma8DqeOjqFiFdthYZOYK/pFM8/2Quxy2tC7WhjurWFlWTy1wdKiZYh3PbEd1d2s5gw5Yq76LR2jDVZY5xrPUP3eq1shrXHs+suDVPdQuR2UdQf2erzSOtUWku8bugMo9wAvnzvn8efkEDAABoQoEGAADQhAINAACgCQUaAABAEwo0AACAJhRoAAAATdwZsx/irVfbEW8dnE9faHR/9fOcfAYxkjlkkMa0/BRP+rQ+53k6l+otsiGmOvUZr2WKCp/1V5zHPuf6qOdKf10isb9kk/XrDOtCeWeO8B2xZdeO2U4SG9ahs/qn2rjeh/EK5xCbvNfXd4d16GpXbi9y8bYd6RmvRPDn96Fau/g6URwvmq1DxXlUtxCIa1TaJqTwjljeoWnxPgd+QQMAAGhCgQYAANCEAg0AAKAJBRoAAEATCjQAAIAmFGgAAABN3Bmzf45xez3PMkWQVqJLq3Gn1bj8HX2Wxor5nrVj5fEKjpd0sPrZxKOlPkvj7Uj8rV7+cE0WJ72+kacbDm1LSA7bfVRjqiv3ZmqzOm76DZU16qh+QNXzvtCWeOsUs12NXq48P+nPqunWWhx7/zl9Vm6T5RHiK27V1fd7k+dnyzymfT7A1gJF5e0rUp/VtS3MpRIrX932orhE7dnmaKJ8jRffy35BAwAAaEKBBgAA0IQCDQAAoAkFGgAAQBMKNAAAgCYUaAAAAE3cF7N/HGM8P7966Hz9n990Pk1iKUPM67TNG+1K83jLc8oufX0ucaTQX2z3FNqlc0vzTyaR07c/SDn7RSlm/7kYfR2affrvX79+YEO8ddqWIEW9Vp+3Usx+uv7hmpRiwn+IY9Tv25lZhPriLTbe9O2na8dbrbgVyOpY7+OWHrra3hD52VkfFT4b74zbjgTp8hf7jFH6xe05Kv2Vtgn43I/sOMb4qrgQ83/lx2r9tiln4Z1tjDci8yftqs9VNVI+voYU51JpF9tUP9K0tdDqta24sKXvncoa5Rc0AACAJhRoAAAATSjQAAAAmlCgAQAANKFAAwAAaEKBBgAA0MR9MfvnmMeeVuNQK6pjbYhs3TJeF7Po8aQYN7/F6sv/4B/nGOPScyhF+v8Q5znGpw3bObymGpPMuypvm5Ju2hShvOFPndNY5h3raNrSI4134Zqet5kpdPi5c79yHfrAjpfaInumLVDCu1caL0XwH/F9aLZdVGiServ4GS+1S9tepJ0MNmzpEdevyvZIMbc/9Bcuiph9AACAB6ZAAwAAaEKBBgAA0IQCDQAAoAkFGgAAQBP3pTiO87FTHL9U6XpVkhp3iAFJKaIntEtJPFtikmZj1ZrFFMQmp/Y+jnmyarrXUxrrrF2lzVvtkmq7xVJCWkzhajL/PUmH1WSv6njXjVVOSCu2qya5lebRTXVNudIjzPFieU18/ZqU78vUrpiUXf3UpimI1XlUUhXf6HP5seq5LU689QsaAABAEwo0AACAJhRoAAAATSjQAAAAmlCgAQAANKFAAwAAaOK+mP3jGOPpg9Z01cjsLrHSSZpjiI6tjZWObYipftqQYb843jpGrxanX45zXfxxv0+8ddjuIza7qM2udheuNbPY6Del+V95btW45nIWdbHdlX2mr7H0saUui+1K55Zun0dKfn+EmPpHmOPFKmtieaefDetQ9RmZ7kBTXOqrcyxfr8VrTVJeDyc+aLUFAADweBRoAAAATSjQAAAAmlCgAQAANKFAAwAAaEKBBgAA0MR9MftjzGPZq7HM947z1ljVdlfHaXex+nO7pYM7curXOybnUJ1GTL2/VTstblkQP5/7hyrZ+cisfh5T/Hs1Nn7DmlGKxW8S6b/D6rjj7xpW74XUZ3UyBRdHR18afV/dVuGzx33w7/+Zi9evL1E5Uv5i03luiNKPOm2BUlA5b7+gAQAANKFAAwAAaEKBBgAA0IQCDQAAoAkFGgAAQBMKNAAAgCbuj9lfHcu+epyr5vdRVLclqKhGx75XhPJ7jjVGPUo/9lloU70NZmM1ihCmoFM8/2wuIaf6KN7Q5ej+xfH2W77hqmvshdHX5VjySkz4TiLsf7h4PdavNaXtSpLV351jw3Pwlibb7JRfea78SNMHUJiIX9AAAACaUKABAAA0oUADAABoQoEGAADQhAINAACgCQUaAABAE8d5R6zrcRy/HGP8h33TAb4A//A8z59WG1uHgAWsQ8B7m65DdxVoAAAA7OP/4ggAANCEAg0AAKAJBRoAAEATCjQAAIAmFGgAAABNKNAAAACaUKABAAA0oUADAABoQoEGAADQxP8Gv387WU9zjdoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 864x360 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "approx = outer_product_np(A_A1 @ A_A2 @ A_S2, B_A1 @ B_A2 @ B_S2, C_A1 @ C_A2 @ C_S2)\n",
    "print(\"Relative reconstruction loss:\", np.linalg.norm(np.ndarray.flatten(X_np-approx), 2)  / np.linalg.norm(np.ndarray.flatten(X_np), 2))\n",
    "display_tensor(approx)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdYAAAFwCAYAAAAFX1qpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAALIElEQVR4nO3dTYxddRnH8efM3JlpZyqhhULpCyUlRiFGUEyNEGM0hCpGjC4wcaVBF2pc4E5Z+LIwcSNLNi7cEI0SxZhgGwETFxhIUCwIWkR5b5NWQUthWmbmuHDXTJhe/P05c8fPZ3lz+pwnzTnzvaedudP1fV8AQMbU0AsAwEYirAAQJKwAECSsABAkrAAQJKwAEDQa5+ALt033l+2ZabVLHTk832x2VdVinaoz/emu6UkY1AXbpvpL94x1WY/lqcNbms2uqqruLbg8G/6InXts/Wr99TuldQeSTtZLJ/q+337262N9Bbpsz0w9dGhPbquzHNh5dbPZVVUP9vc1nc/wLt0zqvt/dVGz+Tfv/kCz2VVV3dxc0/lVVf3p081mu8fWr9Zfv1NadyDp3v6uZ1Z73T8FA0CQsAJAkLACQJCwAkCQsAJAkLACQJCwAkCQsAJAkLACQJCwAkCQsAJAkLACQJCwAkCQsAJAkLACQJCwAkCQsAJAkLACQJCwAkCQsAJAkLACQNBo6AUg6a+L2+rTf/5Ms/n9DRc1m11VNffbx5rOr6qavmRHs9ndcV9SwBMrAAQJKwAECSsABAkrAAQJKwAECSsABAkrAAQJKwAECSsABAkrAAQJKwAECSsABAkrAAQJKwAECSsABAkrAAQJKwAECSsABI2GXgCA9k71fT18+szQa6xptHvX0Cucu+dWf9kTKwAECSsABAkrAAQJKwAEjfXNS0u1UieWT7XapaYWFprNrqrqXvU+YqPrjpyp0fXPNpt/zwt3N5tdVXXj7muazq+qWjl6rNnsvl9qNhsmhdIAQJCwAkCQsAJAkLACQJCwAkCQsAJAkLACQJCwAkCQsAJAkLACQJCwAkCQsAJAkLACQJCwAkCQsAJAkLACQJCwAkCQsAJAkLACQJCwAkCQsAJAkLACQJCwAkDQaJyDl/qqf660WqWqm59vN7yqatH7iI2u2zRX0/ve3mz+rUdnms2uqupm21+j0xdtbza7O9b27wcmgdIAQJCwAkCQsAJAkLACQJCwAkCQsAJAkLACQJCwAkCQsAJAkLACQJCwAkCQsAJA0Fgfwg/AZHrh0YX6xr73D73Gmu55/pdDr3DOZneu/ronVgAIElYACBJWAAgSVgAIElYACBJWAAgSVgAIElYACBJWAAgSVgAIElYACBrrs4KffXRLfXXvda12qUMv/rrZ7Kqq/Qf+3XQ+w+sXT9fyE082m//kh9/WbHZV1e1/OdR0flXV1668vt3wpeV2s2FCeGIFgCBhBYAgYQWAIGEFgCBhBYAgYQWAIGEFgCBhBYAgYQWAIGEFgCBhBYAgYQWAIGEFgCBhBYAgYQWAIGEFgCBhBYAgYQWAIGEFgCBhBYAgYQWAIGEFgCBhBYCg0TgHd11XU5s2tdql/rXyWrPZVVXLtdJ0PsPrZkY12r6j3QlGY90yY7vpgS83nV9VdckNs81mL//mvmazYVJ4YgWAIGEFgCBhBYAgYQWAIGEFgCBhBYCgtj87AMC60E1P1/R55w29xpqOLr869Ar/M0+sABAkrAAQJKwAECSsABAkrAAQJKwAECSsABAkrAAQJKwAECSsABAkrAAQJKwAECSsABAkrAAQJKwAECSsABA01i867zfPVf/Oy1vtUj9/5W/NZldVvbT8ctP5DK9/famWjh4beo03bd9nn29+jkMvPtJs9v4DJ5rNhknhiRUAgoQVAIKEFQCChBUAgoQVAIKEFQCChBUAgoQVAIKEFQCChBUAgoQVAIKEFQCChBUAgoQVAIKEFQCChBUAgoQVAIKEFQCChBUAgoQVAIKEFQCChBUAgoQVAIJG4xx87QffXQcPfqvRKlVVLWdX3fKn7lDTEzC4aw/sr4MHfzr0Gv+3Hj7sHoOu7/tzP7jrDlbVhe3Wae5E3/cfHXoJ2tkA1+ikc4+tU13XHa+qZ4beY4PZ2/f99rNfHCusAMAb83+sABAkrAAQJKwAECSsABAkrAAQJKwAECSsABAkrAAQJKwAECSsABAkrAAQNNZvtxltWujntmxrtUtN/+NUs9lVVYt1qs70p7umJ2FQ0wsL/czWdtfo7Attr9FJ5x5bv+a3zvbn75wfeo01vfL45DzvnayXTqz2IfxjhXVuy7a64qZbc1udZesPf9dsdlXVg/19TeczvJmt22rPV9pdo5fd1vYaranptvOrqlaWm412j61f5++cry/++ENDr7GmB66aHXqFc3Zvf9eqvy1oct4aAMAEEFYACBJWAAgSVgAIElYACBJWAAgSVgAIElYACBJWAAgSVgAIElYACBJWAAgSVgAIElYACBJWAAgSVgAIElYACBJWAAgSVgAIElYACBJWAAgajXNw31WtjPUn4K31rguO10Ofv6PZ/P1HvtRsdlXVBT97rOn8qqru4gvbzX52ttlsmBSeWAEgSFgBIEhYASBIWAEgSFgBIEhYASBIWAEgSFgBIEhYASBIWAEgSFgBIEhYASBIWAEgSFgBIEhYASBIWAEgSFgBIGg09AIAtLdztFjf3P740Gus6eN7PzH0Cufu6dVf9sQKAEHCCgBBwgoAQcIKAEHCCgBBY31X8JU7j9dD37mj1S71sR9d12x2VVX3qvcRG92Rw/N1YOfVzeav3NJs9H/nnzzZ9gRVVQ3P0fdnms2GSaE0ABAkrAAQJKwAECSsABAkrAAQJKwAECSsABAkrAAQJKwAECSsABAkrAAQJKwAECSsABAkrAAQJKwAECSsABAkrAAQJKwAECSsABAkrAAQJKwAECSsABAkrAAQNBrn4L76Wu5XWu1S3ZaFZrOrqmrR+4iNrtu8qabecUWz+d/9+g+aza6quv3Oa5rOr6qauuTiZrO752ebzYZJoTQAECSsABAkrAAQJKwAECSsABAkrAAQJKwAECSsABAkrAAQJKwAECSsABA01mcFAzCZjhyerwO73jP0Gmv6/t9/MvQK5+zg3tVf98QKAEHCCgBBwgoAQcIKAEHCCgBBwgoAQcIKAEHCCgBBwgoAQcIKAEHCCgBBwgoAQWN9CP+Thxfqxl3vbbVLPXXnrmazq6oWb5tpOp/h9a8t1sofn2g2//ar9jebXVX1yd8/3XR+VdUv3tfw/fTSUrvZMCE8sQJAkLACQJCwAkCQsAJAkLACQJCwAkCQsAJAkLACQJCwAkCQsAJAkLACQJCwAkCQsAJAkLACQJCwAkCQsAJAkLACQJCwAkCQsAJAkLACQJCwAkCQsAJAkLACQNBorKO7rrqZ2UarVF2+43iz2VVVJ2aWms5neN3sTI127G53gpnxbplx3X3s/Kbzq6pe/tSeZrOX7rm/2WyYFJ5YASBIWAEgSFgBIEhYASBIWAEgqO23OAKwLnRTUzW1efPQa6zp5oe/MPQKY/j2qq96YgWAIGEFgCBhBYAgYQWAIGEFgCBhBYAgYQWAIGEFgCBhBYAgYQWAIGEFgCBhBYAgYQWAIGEFgCBhBYAgYQWAIGEFgKDROAd3c7M1tW9vq11qavqVZrOrqrqu6XjWgf7M67X03PNDr/HmfaT9KT73+CPNZn/vDy83mw2TwhMrAAQJKwAECSsABAkrAAQJKwAECSsABAkrAAQJKwAECSsABAkrAAQJKwAECSsABAkrAAQJKwAECSsABAkrAAQJKwAECSsABAkrAAQJKwAECSsABAkrAAQJKwAEdX3fn/vBXXe8qp5pt05ze/u+3z70ErSzAa7RSeceW6fcG02ser2PFVYA4I35p2AACBJWAAgSVgAIElYACBJWAAgSVgAIElYACBJWAAgSVgAI+g/eT0W9Kc/gdQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x360 with 6 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "display_topic_matrices([A_S2.T, A_S1.T, A_A2], S_trues)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Standard HNCPD, Rank 4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Relative reconstruction loss: 0.4600117730144047\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA2gAAAEkCAYAAABaExIDAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAnGklEQVR4nO3dzY8lWXoW8Ih7M6uqezx2a2YaY4zBAglvQSz4A2CB2MLOspFgYwxigcS/wAZ5a4mFhVgAGxZegPiwF1hYwniBkLGQkEBgC/mTmWm6e3qyMvNGsBjMeDz5Pln3zYjbJ6t+v2WdOXHOjRtxIt6808+Z13WdAAAA+PwdPu8JAAAA8C0KNAAAgEEo0AAAAAahQAMAABiEAg0AAGAQCjQAAIBBXJ3zP/7Kl47rD//Q9V5z+Q6/+umXex3nsG3AOpdNx49CrVp3S4esp5KmuEfJHOaY5tLpd/3p/RtMiCeb00U5yNYZxRy/efvRdHv/Wbq6olHWoXSau1/P1f8ZYx1a0um94HryaL/Ci6/fbT+RNMfD9vfjeqyuhe6J3EF7Har6deefboCH//nm9UfT7d03nsU61PVp89o7hu9hCef6EPrdF78H/PrvfPjmE9vbxrfWi49P9eHSJbs0Hy5J94G19VjJHp8tTqUxz43fvb559/F0e3r4feisAu2Hf+h6+uV//UNnT6DjR/7dj5dt8QXnainb7u/rl58v/fP3y7blqv5ClnAGD0WtcrgPi9h73Qu0burMcZryS9pcrDvf/4tfrTt1pS/8EF5oL71IdDQXzTW0zYMUaNUcf+m//cyTjnvJdehP/cJfLdvW8IQ9HOt16HSqr9mv/ItXZVt3HToWtcp8qq+Tz/5wPcf1WI91CHXRff3Rpqubui0Wi4U/9s9+o2xbr+oPMN/XL1TpXl3fe1kf8+a2PmZw+tL3FAcM69NSX3dR86WjLiKnaT6FuVTrdnf+4TlQ3ab/4Vf/QW+s/+eS61DXL93U1/MpVCMfHF6XbTdhAXhVvRhM0/TV5b0H//1v/vRPln2S5t/ipzlcYmtYR8NHK/+w/kf/5dfqPi/r83i4qV/M1mP6cOG94KZenNdXYZGtnhPpb4nh2ZKsV+Gg4ZjxnecuvOR23pXC82NKz4/Cv/+1f1S2+b84AgAADEKBBgAAMAgFGgAAwCAUaAAAAIM4KyTkkuJ/uxf+C9D0H+7H/3J0D530tPDfHybxP5gN/1HpEiYT/2PUqt8eARUx0Sh98GcQEpKE+V90hs0gk/I/3B0kxORNLM31pLt+pcu529bqs0PbLuO9i7r3T8rfOKSHUt00d//EW32GZkZIDCxoPIffFrdT/UKxhMjom7UOVrhJSRrBTZH2czg/V+FbuqGE6RpL+RtpnsVcUnjFmubRSkB9pF+3bRDt8LOtP3d6H934PPoFDQAAYBAKNAAAgEEo0AAAAAahQAMAABiEAg0AAGAQCjQAAIBBDBuzn2JS15CvOm+dRT09ks7faUuf7dCILZ/yHEOabo4K7/Q77FDzx+805eL2wujXC8bsd0caZY6teYyyjcEbyJ89rUPpoM37+JI7iHSPN8ouJ0vKsA57maR+o8Q8wxs6xj0S6qbrue53muq8+dTvRZFTH98zgjBUXGvS1hBpLmm8i+57E7fLuNgs+ltiXNh82naia3h8zCmCv8EvaAAAAINQoAEAAAxCgQYAADAIBRoAAMAgFGgAAACDUKABAAAMYtiY/ZGklPetU5JTTGdzl4AYD3s41QddUuR/Nc8UUx33Trhw3HQYb95jq4DGPJKLBtWnc7Xx8UZzvKojpdeQ5XwV+4UBu9tlNKKj0/FiW4gZXlP0dfeYoW1z3S0gdtjuY/PtKA7N7UpCvxhn3jknaY5BawuH57PbR9sHh9dl2ymcgA8O92XbzVqvba/CS8rdevPgv3fv77gepnU0rVFpB470Hlj0i9vQxPnv0K/b1rgn12X7Na+9tVB4j936eGmbrNZQmx4NAACANgUaAADAIBRoAAAAg1CgAQAADEKBBgAAMIhhUxznkAaU0tOeu0MdkJTFgK6QwpdCF0PC41zNs5uA2Ex4XMN4c4rL66QFpY8WzmPUTiZqjldJ8+8mSnWON5i01qTLa481qp0mW7Wl2yNeD71+e7RtLn+pY7Q9l/unuf52xHvjmZyuPaSkxiWcmPDon+5CfOJxql9gqrmU7xKPaabhpvUkhQ921qH0DnLxLONuenL1ucOzP757BWs31TolTaf4zYY1Ja1vPJZf0AAAAAahQAMAABiEAg0AAGAQCjQAAIBBKNAAAAAGoUADAAAYxLAx+0lK6E2R8l0xpnfjCN+QYNuOlV0PISo8xIKmfmVk63LJTOxpmrvjpVjWKrq/GwfcnEfud8Hs6BgF3phH9zN/DpZTfUOmT3EK18qyhJu8E5f/SFu1JO6wVHIJ3dj+tFQeev3m7tYj1cXXfXzEyO+i4R24/j9ZXpRtp3DSDuE6ulmPZVuK4P9kefXgv7e30Wi+D8Xr+ZLbfaTjxfs4bXvTfZ9o9Lvsq147ur+l/V627dYifkEDAAAYhAINAABgEAo0AACAQSjQAAAABqFAAwAAGIQCDQAAYBBnxex/ti7Tf3r9+sG265BBegqZp8ciK/X1V9+rJxKSLG+bMfUpQvWQIrPDcIf7hwdMY91+b4jEr9Ntp8Nd3XZ6WbcdX9fjLdd1vyrHdvlCGCzFjHZjTRvRpY+q5tKdfzd6NfRbj/WFPp/SxVz0625XUB1vqhP409zfdcNE318w7v9J45V9mvdj1zPaOgL2cgg3efWu197uoPvK0F0aumsbb+5tXpureYTpeVMCAAAYhAINAABgEAo0AACAQSjQAAAABqFAAwAAGIQCDQAAYBBnxeyv01xG5h9ClP5S5W1PU5XWHs3heGs3e3XrmOfUb4+46T2iagdJJ30WLhzl2o71HSVy9plY01qT4prD9gnr0tsaIn3nrQj75vG6UdQXXdvC9g9xa4vUr7tdRrets90H/D43a71Xzmmtr/Wbud63Jx2zNZc9Lud0W6XxmnOJx+RztV6F/ak2fh9qbSGUHkf9qQAAALAlBRoAAMAgFGgAAACDUKABAAAMQoEGAAAwCAUaAADAIM6K2b+alumDw+2DbdchKvIUkiyPVb9D3SlGmoZ+Ka45JM7Gthjneni4cY5xzb2xdomVbUTHLtch0jQdr5t2ukO87eFuacyjEZf9BHHnik6sd3OOb3u88Byz7UNTWIfS/Z/OZ7etnOel14wLrkPtey7unbDDFhWd6P4LrzXtP+P68+/n7lWIy0/31av51Byv7lfNJb5fJelSb+7sFOeS0to7nyG+V144ez1u3VGc6DSP3uXTm8djnvHWQpZQAACAQSjQAAAABqFAAwAAGIQCDQAAYBAKNAAAgEEo0AAAAAZxVsz+Ms3Tzfpw1uhpqqPJTyGz+VhFZ8ZM6bopt9XHjGnazeTlGKd/QWn+m491ChH1xbYD0zRN05Jiyetua7lPwyPCeM9eJ7qbB82H+nqew3pySNuEdK+9jbfZiEtsc9uR2HbJmP2uPWL2Lx3dX+nGeo8yx6byunsHlsJX832z3+Vi9tmQ5/t3mdPztrG2rcd6v4X4/lsesG7yCxoAAMAgFGgAAACDUKABAAAMQoEGAAAwCAUaAADAIBRoAAAAgzgrZv8wrWWM6nVI97wLOZKpX6kbybxHWnCMt364cd5lIk3dLQsuKMZzd6XI//siKlWE7TtnXeq/YaXbY17qa2XpXtDNe7XapqK7tUhnrEcNstY8C93Y+5QAHbaFyP16U2mNlYR5lNfkO3DN3az1K94pnLTrtf4iqq2WHp/LdavfKOJ2P1sPFu/xHbbEGGRLqD2s4V1v7q435QEbz/bQxS9oAAAAg1CgAQAADEKBBgAAMAgFGgAAwCAUaAAAAINQoAEAAAzirJj9++kwfXV5+WDbMQSNnkKOZNkvxF+mCOuYANuNlE+x0nGeRccdxto6gvuxtnKoq5R3nPJEmxm2KS4/CbGywvT5PXPMVn7eV8qzmX4zcbo+3tsbKQ2/X3ovSy9ZqV+37bB5pjmfq2fydZbv4ckgzwi/oAEAAAxCgQYAADAIBRoAAMAgFGgAAACDUKABAAAM4qwUx//xtT80/dg/+dsPti3Xdb/lKiTmFUkwX/zNunZMYUApmez0cADlNE3TdPPl0O9F3ZZSxK6+Uc2j7vSf/85Ph8HG9xd/7kfrxpTweN+LBFqve39jmE9hvMPDx1yP9fc239XHm0/19R/nn/q9qPsdXt+f3W8O538tzsdj89jLr371w+lH/uHfeLBtCSvaGtahKnzs/d/orUNpXUjr0HJdzzGtsWv4GtbicljCOvSXf/zflm2fLfWC+OpwV7Z93/GbrWP+wPVHZdtv3n3w4L//4s/92bLPFO7jdM/NIdlreVHHB89fCF94WveK1Nu8Dp3qtjD/dQ7X+VrPcQlzOdyGNeXFwzfqvNTzTynAaR5r1e+JKab/5bc+nP7M3/vJYszQMbUVX1E3TLadar1Hwmsx3t/6iZ8tu9yFD7CERS8lRr6a6zXqFH6zSHO5nh++bn/2Z/9c2edwW88jXevx3WUJbdf1w3G+CXO5oPDVZGFtO335i3W3amlI5ziZX9VtxbNlDbWOX9AAAAAGoUADAAAYhAINAABgEAo0AACAQSjQAAAABqFAAwAAGMRZMfvTOtXRrCktOMW5FjmXMcI6HC+m28Z59PrFpNpGZO5bLUShJjEeOvYLjen6GuTPFulzp7ZR5lHGWz9VWofax9xprmcaZBrTMSzAx3DzHMIXk6KvU780l9SvFK/L7e+r9J2GdPvLCuck3ce73ePPwTqVz5F4WjqXX3o/STH74fra+r3mUUW/F3O9LUyyxId4La0n6YMv4Ysr18SwRc0u0sXQvVdHuccPaXuUZiz+JTUuhVEeDwAAAO88BRoAAMAgFGgAAACDUKABAAAMQoEGAAAwCAUaAADAIM6L2Z+mXgpxM8K+NYc9EpS74209D77Ls9+y4IJx+W+VxhYW6VS3tsTYYT3pRl9vvaacQj733Xos2w7hJKd+KcI6zSX1K134nrPdx9ur+vxr68VgCuta2mqmHituVxTE2Te/8uo+uF3r19C4ZqQ9BILr9dQa7xRuyHK7j0tH1O8RpZ/i7TvH664Z6ZjHQRbLjb2dnwoAAOAZUqABAAAMQoEGAAAwCAUaAADAIBRoAAAAg1CgAQAADOK8mP15qku6lMSZEjeLturfp+kJyaWX7nep473D4nXyHBKgLx3D+5bb/Cvf4RpqR/eneP7GeCnt+LPlRdn2eqkfG0u4IV8druvxTvV4nx1ftvpd0kXXmmWHwd7xuPyWeZqWq/PX75QOX15HSz1OOl5Mot96+6DHFJH/aa3pxuwf0v4CoSnH7Ncn5Xouovv3iJTvSnOJ89x4K5Pu+pWu5XTMtE1I57PFrVHS5z5/KL+gAQAADEKBBgAAMAgFGgAAwCAUaAAAAINQoAEAAAxCgQYAADCI82L2m3JkfhVLmTr1xkpt3WNuvYXAWy1Gx9bxpGszcjZed8/gTxPpc8dzsnFEb3see0nbfXQjp4vrbz3Wny8lOadrb6mTnNuR2antULSlPl+5/qRs+/Twqmx7Od+Xbd939Vnd73BXtn149XHZlmKxSzusQ901vbUOHcJgRdr34xPp3ePv9Dq0hjUgTSesG+XrUNpio7kOtW28I8OHV/VaE2Pv84Je+sLhddl2m2L9w8Nl85j9JF3rabzuPbLHth6VGM/f7NdYt+d0rlKUfrdf1eXsHgAAAOxCgQYAADAIBRoAAMAgFGgAAACDUKABAAAMQoEGAAAwiIvE7I+ijLB9Qr9OjG13Hu+qOUSoOpUP2CPatxC/m52irw930/Tebz987Ji83JhPSGSe5hQ/nJKQQ0z6/ft12/IiHDOlnRff0ell3envfum/1wccytcf/Nef//RP1126MdVJN25+CdnR1f4IV2Frgde39TTu6wz+9WV9cc239RYIh/frLRfmT+ttFdZXLx/ucxNuuOp8PDKPah2K9+8bWOdpWq4bHbe+/NKll7ZdSO81YR1tv78U/X70i19tHnB8/zg1pnXhmL6A0C/cq+tVfcz5vrlGXVK6OULb6f26zKmexfNd/ZnTupGe7eV9E/r4BQ0AAGAQCjQAAIBBKNAAAAAGoUADAAAYhAINAABgEAo0AACAQZwXs79OdTRrSpwOKZ3zUsRcpijXXrp1zmTfODp2msTpf5cLxr9P0yPnPyXHXvLPFuGcpAj71La14bY5WEN89NbR0WlBSXHZ3S0GujsTxLk0j/mcde+PPWL2Q4xyXIeegXSdv+2X3Tw9EmOfOlY2XlDb71F7LOzP/Fp/Z134va3jqVtmbGbjLVz8ggYAADAIBRoAAMAgFGgAAACDUKABAAAMQoEGAAAwCAUaAADAIM6L2Z+nOiI2lHop3rrKel1DNHGK7U+R0nEe3ZjqRr93MvZ6mnIE6Q6B7ek8z6P8aSKckxRh3Y5xbxhlHt8ec5qWauVK61C6V4t/P8a1Jl1gdVP3z2Jbrxvv7Dr0XBTPwBhtH6/J5lYA7S0jGuM15x/vxZ3W+nWepvW8N6j/369Sxtunx2Na19IWR80tSbqP6mGeuZd06edj9/5/5lK90HLhd9XKu3jLAAAADEmBBgAAMAgFGgAAwCAUaAAAAINQoAEAAAzirAyil187TX/in359r7l8p2NIqrqvo4nW0G8Xh7rGreaS5v+X/tVfqcc6nVrzmNaQOtPtV/jtv1/PcZ7vw1D197aEtuPhNozXS9v5xn/8yoP/vh7q48VErBQwdt9LAlyOYcD5/boppHqVwgdI56Ry+z+f9nehFx/fTz/4b373Scf4LkVq0/Ledd0lrUPpvkq6p6azDt3V9+pf+JW/VrbNS/2dt5Mtg5joe3p4Lu//zG+UfQ4hhWsJk7xf6nN8daivhas5rNvBr/zyn3zw31PCXncdmk+daMFpWo91t2n+oG4rTldcsrupg8VHu/mpNPnHvfj4NP2Rn/vaw0OGZ2dM4Sz6pfea6h6Ypmla3n9RtiXLdX1C0/2frMeHj/nnf+yvh06tobLuK2IjSfOjn/qkNdR1iA8+LfUHuArXye19/f6VHjw/+L2fPvjvaT18dbwr2+7DopHe9a7Cy8vNqS5ljtPvlG2HYsFJ86j6TNM0vQ7zKOf/E/X34hc0AACAQSjQAAAABqFAAwAAGIQCDQAAYBAKNAAAgEEo0AAAAAZxVsz+Ok/T+vJp0bRvKkVY826JycshDjXF7Kd+Oeu5OF43ujf0i8dstlXHbO5IkO2540WKc3/Ouste/ALPP1fzXT2RGLMdIvG75rR1ShEx3o3ET/HK7bbwd9AU61/drCnCPUmXSIznj3PceC5xsU/XVvwA54/1Jta1/V3wbYfuu17qNshPD+k9I24tFLavSctvd7ykWktPYY1dwrYvaW2O8fZha6FT3Gejthb3b1qXqz7T9MhzoPHOMshlDAAAgAINAABgEAo0AACAQSjQAAAABqFAAwAAGIQCDQAAYBBnxezPyzodbu73mst3WK92qB13iOZe0zG3Hi9El8axUluMxQ5tRdbr1THEczdj7+elbjuGONpjiNNOvlENly7JZsp5imxNmqmy01z0i9NIjZ15bHFbvK0x+8kF/5wWo/RTXHvq1/zO1iJKP413CGvNIX6AuikeM7bV61D8SotjdncIif3SPJKtt/uIgzUj7auDPnUJmefy+S9+/wzh/u4fc/tDdqR3niStJ+mdM4239VyW5pp3Fd45Y8x+uP/jmh5U5yS8Vsb3ykN4MevM0S9oAAAAg1CgAQAADEKBBgAAMAgFGgAAwCAUaAAAAINQoAEAAAzirJj9aZ6n9fiM46270bcp1jQcc+3k+KbY+15q/CNR+kGM53/4n08hEj8dLsXsn5b0d4R0Unp/f6jSY9fu+U9jhfMVY7FjnHkYsPgMMYE3TGTtRPc+NVV5Waf5s5uiLXxJ+QJ8uMvxC/Xh7sNYKad3D2ldLmKs57s6i/rq42+WbfN93W891vfcfKrP13p1LNs66/avf/RB2XYM12yKeY4R0M0I/uTqGw+fyzVcW92Y/cOptw6t4Wu76DqU7rcq0v+pUezrOs135287FJ/G1TqUnoHhvppve8/AQxhvDrH4MQK+2ILn8Emxlr8FPvroi61+8zGc4/DOcAj3wRLu8eTu9PBNvoRr5PpY31yntKAEqdf9KVyvjeHizkLhHJ/iPB7uV53fafILGgAAwDAUaAAAAINQoAEAAAxCgQYAADAIBRoAAMAgFGgAAACDOC9mf5p6mZUbS1Gu0Q5zb8+lkuLTTyEX+BBq7RCnPR1TvPX5keVXRZTuNOW46RRhfQyxpsdDiBhuxluv1anc4dJP8dAx3jr9aSXF81fHi/HWvUj/cv6f/xLCNO3z57m0HnbbGlKUfloX0ilJt0iM2Q9rVFKuDZ177rGxUr+0jjbnUn3d6RzPqTXOv9FnT1tvpbPHO1k4Zvcaa80znqvmVibdbYca481hQYlR7uldKRxzTu9DcU+MWrmWhrHimhc+QH4PrI+Z+iXVsyBt+5SeH/HZ0th6xy9oAAAAg1CgAQAADEKBBgAAMAgFGgAAwCAUaAAAAINQoAEAAAzivJj9dZ3mm7udpvIHhnp1XbbNIZ90jbm/zVjWEA8b59IYbn0ZvpIltO1Raqfo/uXhyNO7+7pLiiBNsaY5erUXh5rMRZprP144jdXIxJ+mdvR1FePcPFU5gb9qa47FM5AuiD3aCmnNSEtl6hfXqHDMuRsBXXzs9u2Tlpp0+pvP1PipL7kGjLYOdePhn7vOC1H3XKWxertebD6NdP2ltSbvv5PWrzCXoFoT0xy773N5t5/eeB3d4+U5nn88v6ABAAAMQoEGAAAwCAUaAADAIBRoAAAAg1CgAQAADEKBBgAAMIjzYvbnOcbfv7Vilvi28Z7z65BTX0TbP2ke3X7lOdn++uhGr3Zj9quo2rmZyxx7NdPy2ztGlJHTW18HbK7757Sq36k7kfFdHesPF3bmiKney1J/AYdDvTYfD73tRcqmdB2k2zGlc6c5xn698cqmGPffG6uc/7aP7jfXvQCfu43flZ6DOdxXuV941wvrUBovbQmVHIu5HMLCcHWs59/dbimtoym6PzkUczmleXS3b2r08wsaAADAIBRoAAAAg1CgAQAADEKBBgAAMAgFGgAAwCAUaAAAAIM4L2Z/mi4Wq70ej3XjMczhUNec6x4xr8cQk1zF6aa0/BjBG+aRPluK578K5/kU+h0fPs/f/8VP6qFCdOx9iI49hSzn60Mdp53GS7Gs//W9Lz/47ymKOmZRh3jVJeyqEKOjr5r34VIdNMXzhnk0YoRjNPeb9L8+Tnc/8MHTDvKGrv73p2Ei4Zw115r5s5vQGI6Z1o2qX7i/1/df1W0vel/geh0eNylevFhrpmkqP8P9z3+l7hP3rwjdQtt9mOJ9c7wP/9fDa1t3R4xkjs+Wuql8xj2m8R4R16HGNH7rs/P7fMeYx8N0/+UvPO0gOzvcpS+2Fh9nYR1qPQbT/d1+37xwpH9xTr7vF95rHS4+I5tbacR7PLh/9f7DxwvzeB1eK+N7bHDa4bNVH6F7+lNpUn038yf1yfILGgAAwCAUaAAAAINQoAEAAAxCgQYAADAIBRoAAMAgzk9xfMbmZiJQTH+Mx7xwklDHxqmcKXExpfesIf4pJS6eQvpjko5ZxvSEPjFhLFwHMeAtfTVlGmOeTNnUvQz2iJR7dMx+Otm55vs6JXSfAQdZM9LF110zutdKSpMt5rKGJ1t3GikpLKauNVMcl5AQ3JIO1w1j7P6Jt/MlxDTZxhyeenrXdTrcNtaHC97j811z/eqm0Ka2at1o3N+PT2SHc5zmUoy3XPXmka7nbpppO+mwWEvjIyKkOLbnscNn21rnu0mfyy9oAAAAg1CgAQAADEKBBgAAMAgFGgAAwCAUaAAAAINQoAEAAAzivJj9dZ3m0/mxpymmvopeXa/DAUPEeIzE30Ma71C0jVQWx/N1/nf93tVd2XYIx1tCQO8xROlfHep81au5FzG8HotrMn1vzejbtL1A7HfVjR9u9Ns63vqp5mlari4z8OF42Q8417fP8zDKNgFNndvjsX7bbmTCMOZ5Wq5DnvgADgPF1K/FXOJII60nW88lpfanr22PnZ22XqTa2/bUTXkLh+ZwxUG7Wxl0518ZqVQAAAB4pynQAAAABqFAAwAAGIQCDQAAYBAKNAAAgEEo0AAAAAZxXsz+PG8eY18er05Pn6YQHVvF9rfn8RRLMZf02ZJuZO4exyw+2+2pjh1OkfhLyC49hZj9qNltLrdx6GWvzimzNVwL8ZIM/VoxsGn+8RLZOFf2TQ2Uvvyg9OXtcR8/Bxc8J2mHjfZS3423bo43V8+P7gdIt2o11pTXk3qtfKRf5+tOz/0wj/p4jTl818CNY19y7bIOfa4621JN0zSt4SLJ77ihX7jH49pwKo65x04/3fum+U7duRXTLTV33svCZ/YLGgAAwCAUaAAAAINQoAEAAAxCgQYAADAIBRoAAMAgFGgAAACDOC9mf5qm6dgIphT1ypSj9GNbjJztHTNqRNHn4+2Qq9xNty/m0o6+bWX6P9G6TofbkKO+ofnU3ROjqb3tRWir/gyXxtpj25EkjXcIbcXnbsUdTztFQLdj9quG7e+rdL7S/NdD2vImHLIVsx/awgH3WH6/deDLrUNd831v/UrbDqWY91a/eH83r/V0zK40l2K8eF9Fvfsqb0XRm0nZb4+Y/SBv99E75tbiZyvmn65Uv6ABAAAMQoEGAAAwCAUaAADAIBRoAAAAg1CgAQAADEKBBgAAMIjzY/ZPW+fjFtLMljqYMsW8XlwV9ZrK4jT/bhR1zCdNx6ybKleHOu/0EDJIYyR+iFBN4101s1fX6nOHSOl2BHe4lnO8dThmmkrxGWK6e/puwjnZLd56nqd1jxjlB1xqnN8zd9evNM/qmHusle24/G2j45f0/Gju+pIilOO13t2do9jSZo/7ao6Z+OF521yHOrtJbH7+n3oeL7gOtXW2RXrE2jxxVb94tJHOb2Muy1V3Pa+b4r2T1rZmFH21lsb7sTv/JL5HNY9Z2eMcF/3SUH5BAwAAGIQCDQAAYBAKNAAAgEEo0AAAAAahQAMAABiEAg0AAGAQ58Xsr+s0NzIyU/R9dbwYVxnm0JnfU8Thqujo9NlOoTENlmKql3TMUKM3oq+/eX9dtqXY+yVkl94v9RxP4UK5D+OtIZd5PlWdUr5q3ZRiWeNOAKnfqRff29p5oBnBv1tA8rpO88ax7JXuOCl++1Jzf5IYiX/h8RpbgcQI6DCN1C91TNf61t92nGPQjudPz53mQVufIe6B0Fibn/rFXHAd6ppb2yLts11R+W7WHat7PVxQ9300vZ90r9vuXOZqLt21srt+tRs31t2KpTGUX9AAAAAGoUADAAAYhAINAABgEAo0AACAQSjQAAAABqFAAwAAGMTZMfvT69sHm+YQDz/HCNwifPLV99R9ll6Eajs6NpWxMQK6aAvHu/uw/twp0jfFercdwzHL+N6vbT+PC1uPRcNhjyzX+hzHpN1jM/O7uk2bUfpxHns5zNPpC/V2DpsO9cl5S+S3O4ZzFu7j+e6+N97W7qu9JqZH9hYJmut2jlAv5tK8PaLmVhpJjvVvxJKnLWjCJLvR1+mctBL421HgrW5Pk9ahtBXFBf8sPt+G97Jwri8Zs3/87HXdqbvWNO+Rtmq8He6r+L2lhWjruezw2ZJ4RY6920WbX9AAAAAGoUADAAAYhAINAABgEAo0AACAQSjQAAAABqFAAwAAGEQzQ5rdhGj7mCQa4/57U0nR/VXLIeQdz6EtJdinYx7CWUn9UgpyldmcYqNTzHPsl+aRbJ2C/JwybJd1On788HYfW5tv7y4yznCuqr0mpmkK2wRE3a1AGpHZ6w5/epzDMWOkfJp+7PdwY15Pett2dNeT9ucuu3S30Gl1e5oLrkNd8xKedClufoeY/XK8tNbsEbPf1Thf8R0qHi9MIzWmNSpus5G2cKr61GPF+Yc5tt+jul931a/72dKLZWOOfkEDAAAYhAINAABgEAo0AACAQSjQAAAABqFAAwAAGMTZKY5zleDVTdup+qU0lGb6UEyqSTEwaS4hgKiTdjbfpdSl1DE0Ls2Im5S2c3p4nqcYp1UfcAn9TkuYSEwECt93mmfRFlOXeoFM7X7d260VyBjOVTwne5nnaX2ZbrztrN2EsWas1Hx/avXbXHF/T9PU/9w7rOllv2ZyYkwRq5v6SY1Bp187qbEbg7ZxWF73u2nN46lzv+A61HbfuHemaZ8UxGqo1zc7HPXSz6XGOpQad1gz4jKaGoumdIa7ydXdNXbr8xXn2E2TrNq2D2AHAABgawo0AACAQSjQAAAABqFAAwAAGIQCDQAAYBAKNAAAgEGcF7M/z9P68nqnqfwBqXQMsfExLrQ7XpLGOxRtIUm73MZgmnJsfzXWNE3tyNkUr17M5RCykI9znbN9aGaNXqVjplzmRls7Ora388Ajc0z9wlwaQ6XrpxXr/dQE53Wa5tOFYpRDFHW6V9eR/vTVWRPbseuN9XCaHtnKJJzMtB3ABbXj7TuX8Q5L/R5a62Xa5uByye9v5pLrUFd3fof0HpL6hbaqX7q/P49tXDpaa2zd1N5uYhTdNaq7Xcloa8NGRnqNAAAAeKcp0AAAAAahQAMAABiEAg0AAGAQCjQAAIBBKNAAAAAGcV7M/rpO0939g01zijtOEaRVjOqrF+f3maZpbsayriG6P0qZp9UpCadqTZGz7dj13mdbUyx2MZe707Hscj/Xn20NGar3S93vdKhP5nHtRXDP5bXQy4Cew/c2hy0X8vfWvBYa13k3gn/LLp+XuO3FHsdM61c3+r6zJi5pkWqek+4am8Yr2sLuG/1E5uauHd14+PI66TxPH+nXv857e4+U56t7bXXujWe0DnVd/H0o3XfVXO7DQ/CS18NjGvdWXuvTWOFejedk+3t8Pm17LuP7RPfx191tpfGql+aRnjvlV5p2z6qbAAAAuCQFGgAAwCAUaAAAAINQoAEAAAxCgQYAADAIBRoAAMAgzovZf+bWPaJXO1JZHKJQ01YGMZ6/GVUbz1Yxz2OIvT+EfNWlmUWdxrtKmadBOZU4xe4WCL22tfmnlRT5X47VjNLezTxN6/EyA7fvq7DWpLnPo6xRXc9g/mmpadwejx4z3iNbR73vcP7TZ+tuIVBqbiGw+TzexAXXobb29kE7XEfVrhEXnkfb1nPpPvsv/dBtxMO317zmWjPSZbIlv6ABAAAMQoEGAAAwCAUaAADAIBRoAAAAg1CgAQAADEKBBgAAMIh5PSOCfZ7n352m6df2mw7wDvjj67p+2O1sHQI2YB0CPm/lOnRWgQYAAMB+/F8cAQAABqFAAwAAGIQCDQAAYBAKNAAAgEEo0AAAAAahQAMAABiEAg0AAGAQCjQAAIBBKNAAAAAG8X8Bi4YpHQ7kbfgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 864x360 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "X1, A_A1, A_S1, A_A2, A_S2 = all_factors[1][0]\n",
    "X2, B_A1, B_S1, B_A2, B_S2 = all_factors[1][1]\n",
    "X3, C_A1, C_S1, C_A2, C_S2 = all_factors[1][2]\n",
    "\n",
    "approx = outer_product_np(A_A1 @ A_S1, B_A1 @ B_S1, C_A1 @ C_S1)\n",
    "print(\"Relative reconstruction loss:\", np.linalg.norm(np.ndarray.flatten(X_np-approx), 2)  / np.linalg.norm(np.ndarray.flatten(X_np), 2))\n",
    "display_tensor(approx)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Standard HNCPD, Rank 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Relative reconstruction loss: 0.7980806524269276\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA2gAAAEkCAYAAABaExIDAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAkMUlEQVR4nO3dza4l2ZkW4BV7n8wqYxd2223oGUj8SH0D3AEDBlwB4hYYMOUi4AIYcgeMGIGE1OoBA5BALSQkaMlCRm5jt+2qcmaevYNBtSXbnO89ub8TEbl25vMMM7x+ImLFilj7uN61rOs6AAAA+PBOH7oDAAAAfMMCDQAAYBIWaAAAAJOwQAMAAJiEBRoAAMAkLNAAAAAm8XDL//j8xbfXhx/8wV59+d22vlrqg2lngFAstvfrUGmqMx2rqgxNXT5vnkDS6eNz5QoPv7qE+pr3NB1MdXar3HrriVRfvCYHluuec+P6f/325+PtY3rAs4dvfXt99d3vd4vf1laaF3ZwenttlVvD1Vwap/Du2+dWP6I95qGi3Ktfvgv13cc81H4mW40deHNiue45N+ahx78cby9f38U81PXwVXdeD8c2Hg67fPNM4vRuj/F8cJ3dua1jh6lmbfwZ6sBP5vHmV/93vHvz5ZNFb1qgPfzgD8Yf/Yt/9vTBjf8W9/3/FLq28c0YY4zv/Y+3dZ0PdaWpvdNj0dFrfQI/+/uf1RWmj7BQ5/VVXTBNIKncUnxH/s0//UWor/7oOz2GD9PwoXJ9Xde5hDrTR+vypvi4SxPVpa4wLfjWUGcsF8ZkOu/1/HR7S+h/kvpR+dP//q9bbf3Gq+9+f/ydf/rPX1TH76sWOD/4s/Chv8MPRd/60a9CnXWl1X0dI9zbML5+8g++V/cjrSHT2ies+ar5ZIwxruFVcHp8+t//xn/4P3Whh9CRS28eGq9CJ5t1Lm8b81Dq4x4/+JzC838N5731D0WpH4U/+dG/6bX1V1599/vj7/6Tp+eh9INJR/rOSM/OD//Lm1Z714fwzVB914zePPSzv/f6/Tv2ntL7vfuN2Jnvv/O/0w/Woanmd0HSrfMa7unWuvdtCZf5sfEDQOzHxj+G/td/96/KY/4vjgAAAJOwQAMAAJiEBRoAAMAkLNAAAAAmcVNIyBij/g8bO/91XKgv/od4jfqedWoWjP+hZ1EktdVcMqf/ADRdyzX0Jf6HzuU4mCiRKV3nZihGR7o36X6v1961jO0dWd+eY6G6fd2uNhJX9wgJaQc5dOssy4RupP+AulnnLscqIVBp++TEHersBoGwrXWU42/zma97W7ufZZMMo136scM1aVUZg8MPvgHN+b6sbuMgjW8q7RVrtdcMdW/1MWUwNaoDAABgBxZoAAAAk7BAAwAAmIQFGgAAwCQs0AAAACZhgQYAADCJ22P271g3LjzGzW+s29aSltrpWGovRcA/0x/YUysWP9j8GT/6Abk2yjQj2Y+cD591ZF/2iLefJRb/U+3/i9v98O3sEkUft6LYodynaKZ5NCi3i/oI7mfn3HbZQqDgL2gAAACTsEADAACYhAUaAADAJCzQAAAAJmGBBgAAMAkLNAAAgEl8+Jj9A6M6lxTFG+Jhl1Odq7l195cUl51iWVO5dCxG5oaudGK9j5Yif+9AGq/rgRm9s/Tjt13PG1dYnUY6vXCsu6XHCOVinemntkYu/hqub3r2U1Opzngtz+GaVGMzXqt44+pjSWqvOxa2ri+9/1Kd3XJH+hD9WMKY3rg7a3q+u98Mqb3wzMWPhvit9HS5dYfvq/jOSs9/EHdyKKqMbR08ZPPUFu5BNfaaUfRHb4WTnp2yzu4rotHHVJ+/oAEAAEzCAg0AAGASFmgAAACTsEADAACYhAUaAADAJG5PcTxXqVmhzNKIbUkpKqFYOz0tJQl103aKvlRpRt+U6bbVPHZknc1knCUVbKa1rSm1s64xtBWONRMv41iO4/z21LWtx/iu1jFO74pjW3cnJZnGVK9mVNUljcu6M2sagFX8WIglWx7r6lKKYzcMMF3nUxjPZT+bib2x3EzHtizzXLl0vdI47/Slm7zbfd5eYFnHOBXjrz2fVm11p/vm9Vwuqc7Ul/CsFn05PTb72B3qpzRme3VWuucWUzub4udEuN+dl2p3vHbrTOP8eqkLliGOzYTK1lInhaLeXh0AAAB7sEADAACYhAUaAADAJCzQAAAAJmGBBgAAMAkLNAAAgEncHrN/LTImY3RpI3O2mezbjhkOYuRmKlj0Jff/fXp0Y7mjjzXKHH1PN9eM0m/bOE57CfXFuP/Uj70i+Jcx1vP2dd707+OZR6D701cqF7eUCOU6829K7U/l0vUKdcZU5lTuyJ8Y9xjP7a00DhS2K2nvq1CV697PD3Ct1hHG5sbdibHr8dNr+22H0n2NU01RZ5zL90jE796buK1KUSTO56Gpo8dztX3WaEb+7xCzn+9b2FKp0//ucqZzbs3XOgAAAAeyQAMAAJiEBRoAAMAkLNAAAAAmYYEGAAAwCQs0AACASdwUs//wy2X80b9/OhM1RVl2IkO/+NGbUGEoGGOe64Ov/+KrUC6cXIqAfrw9e/2zH9a3JEbRB9dzfd6nS13p9aEutxTlzj/+WVnmnKJ7rynXtD52fggZvSnmObV3Le5bGseXS30sSWOr6scYY5zDeae+VO2FMks679SPymPzWv3GOsZSVbF1OnFza4i1u+1CKpcarLY/GaN+DsLzUV7fMdrzb7KE846vj6rcHltzHF1ndawbwZ3mvO7cHMdkaq9RJum+HF9gGfW47exsMUY4jeaWN2n7lCg8/0tzHFXl0lyzx23d/N6EOuM8Grck6J14mkfjdiXd+fdIabeP8B27pHdjVab5jmuN18Y0CQAAwMEs0AAAACZhgQYAADAJCzQAAIBJWKABAABMwgINAABgEjfF7H/UUuR0OLZunOsdYzrbCcTNgscnF8P7qR67O4gEju79J7NZ5oxpsqEnkqL0jy5X3p/mVgAfwDpCZHuzq1vX1xWj6MOz1S3H+0uR+Pcufmt/guPn3j8HAAAAPhoWaAAAAJOwQAMAAJiEBRoAAMAkLNAAAAAmYYEGAAAwidti9pcxrg9PR12uYakXo1cL13NdKEXRx7ZSTO+pPoE19CVGxxbHUlxoPu9eBOl6Lg+N5VIfu4bRcboU7TWv8RghOzbFynajV1NXqkGU2ornFnTr7J53Ve7Ifrw0LXfZoI7f04m3PjxR/sCY4TSPtnuRrmWM524eu3efYKw0T+s+j3cRTX70HNvdAqFTZzq39MnQLZd0h0JjR4w8KMOh5hYO7WFelGvd6+fKNfgLGgAAwCQs0AAAACZhgQYAADAJCzQAAIBJWKABAABMwgINAABgErfF7K9jLJciSDJkYLaSuEO0eoqbj42lcunYNUWQ315n6v9yDcfSaYeAz2vIBU3tlVH6qVyoL55AKpe072mzXMc1DOYUYZ/KdXNlq3OL12rjfrz08q6hju7uA4dn5k+uHZPcOxan7T3ile9B9Ux+1CfNU+L8dOdz1y5bejTbS/J2BsW/N6P02x0J4pYL4Tu23HUotlUf2+ee9hosz22SZ8pf0AAAACZhgQYAADAJCzQAAIBJWKABAABMwgINAABgEhZoAAAAk7gtZn8ZYz0/nUsZo0s7uZopGjMdTOWKvn9TLkWQ1odiuY3jkLvxsP1Y1kaVp2ZjqVyM7u+2F461L3RDt/9blzuyHy+9vEtdx5G3LmrG1E+jGXsfz3uPY0faI97+wPfHLrrb2lS6749kz+tYVb1xk0fHlsf7sPU2NEE87243dthCZGspNj5uM3Wg7vVo39Md2uu0tfmz2Fx6AAAAcCALNAAAgElYoAEAAEzCAg0AAGASFmgAAACTsEADAACYxG0x++sYy6XIwUyxrJ1Y05SsHmJGUwR/2fcxcnTsNeVgNuqM0cr1oW48d4wFbcZblzGwKRI//RyQyqV70z3WaW+PuOZuZO7W5VJ9W0dpvzQleK3r2PwOdfs6RxJyX5x/W8XywT2ObW2PeOvO/DVT/P7W1yS+P7o52zsOkq2rrm553BamPtSNZF+vvfY6E3A7Lj/ZYc7o9HMJ13EN93TZYWLrjoXOee+yPcKB9y0+U81ynTL+ggYAADAJCzQAAIBJWKABAABMwgINAABgEhZoAAAAk7gtxfE0xuWzpzNMUiLNeg51Fgkmb74futZMEbq8qg9eXn+nbu6hLpfO+/zm6Y6m6/Gzf/xVeex6qRtbTnVc0KtXl/LYJdT5rc/elce+fvPqyX//w/8YLsgpxRaluKOQ2vlQX8wlJYKl9l49PfbWkCK2PNbXeFxCW6H/41qXW4s+jjHG8qa+b2V7j+m+hfEf+lGmzXXT2P7K9Vvr+MUfF+fY/cmpiHo6f/30OB+jn9CVk1q/qA8dOA/98h9+WR5Lc8YpRGY9NOehz8M89GU1D/3nv1aWub6qT/z02JuHrq/DPBTqTAlj5XOcUhxDUnFMPw51xnIP4Z0Uzns9P91eTFoOUj9KP37Z79PXz9fxiz9+fPpgTHduNHbupUx//tPPGo3luWEJr7o47xXe/aOf1/WFOMB0LDmHb6UY/hjaW4oH+cf/7Xt1fc0E7baUrJg+X14XZdIQ3yHFcQ3P1PJYN/j418PJVfWFZyr2IyW+F979SX3MX9AAAAAmYYEGAAAwCQs0AACASVigAQAATMICDQAAYBIWaAAAAJO4KWZ/HXWMaozZb6Srx5jO+tDx0aUHirsLNONo07HrDhG3LSlWeqY6KylWvtuPo8vNpppTYm75xn1ozjUp0vze56huhnIVU90ud/Q4T+2l578ZK98Rx116fzeio59t78j69hoLy6jj9Lf+Dpnpu2aSOSrNGd3vk/ht2aqxaYeOxGKbvxubGfyxzl6xw+vckL+gAQAATMICDQAAYBIWaAAAAJOwQAMAAJiEBRoAAMAkLNAAAAAmcVPM/hijzOqMqZop37NRX8oLjamZ3ZzUrfNV29cqxXNvHzkbu3Jo5uwO0glsHct8DW1VMc3Pufsb8Ola0r078rambqRo9bT9xrX+zW8N530N7aVLUpfzfHCAaph1XyEHDts0D61xn5BUaWqvqK69RVDoR5C/azbeWqj5rdr5Zn5O3Kmh2O5qjPBpmc4tvj9CR5JQZzy3dN8659apL0mfh43qAAAA2IEFGgAAwCQs0AAAACZhgQYAADAJCzQAAIBJWKABAABM4qaY/WUNcZzN1PJWfc2U6lSuXWeIJ62zUnsR1kk7Sn/zLQS6ewj06lxShP3G7cV7vYd0bknnHqQy3fPecyuAqk+nbSNwu9t97BKXHWONU7nb56EYN53aOteHUpR+mveuIbq/LBc7GewwR7Wf443FWPUUYd2MY2/HuG9d317z0Drq8bf1PW/Glre39GjXeft9jVts7PA91N1FJ/WzrHOHrZ3iO6lXZW8sHLxtVftdnF4Fnft24HvfX9AAAAAmYYEGAAAwCQs0AACASVigAQAATMICDQAAYBIWaAAAAJO4KWZ/XcZYixjlNSUhh+jlKuYy1ZeiLFMqa+xjai/VGc5trSK/Q32ncy8eejnVF+Uc66xP/OFUlzufi/bingrNKORQrrzG4wWx+EV77bZS9Hu6Jt1ynWN79GNP1bANz0G09TyUrlmjH88d23oeOj9c6m4s9UVJ89CrV3Wdl0td5+uHx9vLdcflDnNUen7WEMfeOoP0jktR7d1yaW7eeG5o17fXHLWMer7ZusnqfTtGjBGP81Ccv8KxtAVH4yf/12Guuaa2bm9qjJG/a9K5XcO1PBUZ8O8a375jtHdbimKd8d3S2Z8mdSQdC/c79OOUyj2kebv450sa42HO7sw1vWkZAACAI1mgAQAATMICDQAAYBIWaAAAAJOwQAMAAJiEBRoAAMAkborZH2PUkZAxMjTF9N4e176mnM4Y+1sXuwcpAjZFQKd06HwsXMuq3NHx1kGMxb90omN32EKAbW19qdNUE+OCmx2JdYa+1EnVZT9TFy+PdT70NUQQL51I5pGj3N+e69fU5bH6jbG7x8ax4hjq2OG0l/hu+cTntvI92CjTaWeM+GGzhG0ckiV+a6S+pDqf/ve3Ya7J3yC9b41r2Hao+z20FCeX5uU0RnbZGSK1l94f1XdUGgfN8R+Lpeuf+v+YTrz49zCPLmnbkXSsEq6Hv6ABAABMwgINAABgEhZoAAAAk7BAAwAAmIQFGgAAwCQs0AAAACZxe8x+ET+ZwiWrCNIxxliqOMuQd9qNt47xmN204EasbIpyvYY+xlTTZhxtOnbtlEsnd3AUfYwY7vSzu19B6sepG8c+Sbx1937v5c630ojS47PxT21xzk7bV4Q6T2Gsp10vUl+2z6k/VkoKn2Uox21yPnWdS3Pg5Uz3Lm2fcKQ0L+RXbnMLgTS3pQT1EL1enkP8MK4P7bIlVLMvR47x5m4S+VOjcyy+yMKxjff58Rc0AACASVigAQAATMICDQAAYBIWaAAAAJOwQAMAAJiEBRoAAMAkbo/Zn1yMJ50lLbgZFxq7H6Nje8dOzXLcoBth/alHX7dijRtjtrkzRzuavHtbO5HNaa5J2z+krUxOdRZ1nGvCT4VpHkoR3XVH6kM59v4Onrn0k2uICW+L4/z2LTja8eIfaj6smm1vn7JxfampcM32iDSvzqD7LZF37akPpvkkCs9WeQ7dKP30HHeHVrO9qlx88vfYPiRtx5DGcnp2qnNLW0zF7+LG2YUi/oIGAAAwCQs0AACASVigAQAATMICDQAAYBIWaAAAAJO4KcXx/GYd3/1fb588llJb1nN9sApEefjyMVR4ezrUGGOsD/Wx81d1e6n/KT3m9LaoM0SWnf7si/LY+VIeGmtYal/OIXXmWvf/y9d1ufPbIoXryx/XHTmFQXJNqXEhfux8rovVpXJ7RT9TQs96CX1M/X9b939c6hu+vKkf3fXdu7rOh6LcY3jewnkvVX2p3PWFcXLrMpavi+vWDJWrvP5lGpfhUEr8CmPv85/8uq4zzUPh2araWy51P/72v+ylYi3pMTjXY2V5V4/19fy6rrTq5rUez0sa6uk5DoeW9I57TBclXMtwTTaXmkpz5TXFv4VyVXuprST1o5L6915tLmP5qpiHtk5xDK+JNC7P75oJiekTq/kdcirKvf633wuN1YeSPA+FYyl9sPHa+sPie3mMMa6vwndgmJuPlt4tR8r3pr5e6Tov1RwQ5obL5+G+vQ3zeVHnT/8ypBSXRwAAADiUBRoAAMAkLNAAAAAmYYEGAAAwCQs0AACASVigAQAATOKmmP2x1JGV/Zj9pyMmc7RnL/bzeg7xmKn/IWY8dqWI08+x/aG+5rEUfRuzwmOdRbkU951i6lMfQ4TyErYsSFGpMQa96mcaB3ELgdD/dE2W2MnWsaq9PMab/SjvzQtje9dRR0t3Y/arcunxSPWl+PRmuRSnHe9DFUF8ZFtjxOjx7jVZqxuXIvHjXBPKxX4E3Tj36lp2tyvZQxxDt29lsks/dlRFr5fj8rn6qkc1PSBhq5zO1iJj5Ej5/KzefizVF6fsNDenmP3eq7N1LWMfU31HP8fJHCn7Ywkdyde5McCaYzy21eAvaAAAAJOwQAMAAJiEBRoAAMAkLNAAAAAmYYEGAAAwCQs0AACASdwWsz/GWIr4z5TEGRObGxHQMZY1RaiG6NJu1GtUnUOKgI1R1CmmvhmXv/WxcI1zlH7zKqfo62YEdNXPmDab+r/WN3yNFzk9OefUm9BeI3I9PsGhH9eq3AZRtJNE//K7dplHJxG3ouDTZEhMKW5XBHfCX9AAAAAmYYEGAAAwCQs0AACASVigAQAATMICDQAAYBIWaAAAAJO4LWZ/HXWc+KnONU3Ry1u795jneK3isRTX3jvzlFR75D2F37Zcxvjsp41tBtKYLQb76y8fN63vuXLLr+v2lvRz2ikcrJ7/S9jSIzQ1zmGuf6y3ZFhe1fdseXcpj8V4+6IvqR/rJfQ/XJO4pUe5pUTuS7K8C2PvSOm8z/W4S9vaHNmPTn3vY3kc47OfPt3u5jHv6fTCabz+1btWc9dX6RkPXQnT8lI84kvYfydumxS/h9I2OqnSVGc4Vjj/up7XTmleSHNzd0ei5rWM1+tI6RlI1+t1eO+U2w6FbqS2HtN9K7YoC8PAX9AAAAAmYYEGAAAwCQs0AACASVigAQAATMICDQAAYBIWaAAAAJO4LWa/qR2Vyu7ufVsCPk2dGOuYoLxxknCKhl+6T1aKm59EjMTvuoefEcN5x7HQiHpPsdebR9uPMc+4m6Ufv6WcNzafT+pjeV7b4ZrF/Xdury7OvWn3oFisd8Hu4Vt18y0cnlPNv2n3kDRn93YdmcbWz2IaVvfw6gMAAPgkWKABAABMwgINAABgEhZoAAAAk7BAAwAAmIQFGgAAwCRui9lfxhgh4ncG7cjWPU6rirgN1zD2Iyyn11Md1rk2l+GpL906mdSEEdZJK/I4lCnPvpuIn+LTUz9CuZjInu5fUTC21ajv2TobkfJjjGfinIs688Xq9SNpXpOOXaL0kyOvZfNaLWGMpG0JXqqah7a+Q3G+a84nyRpf/r2+VLpb/eRrkipNW6A02+v0I7Q2S6T/GKMXi3/nUfpJ995Uz2Lzkx8AAIAjWaABAABMwgINAABgEhZoAAAAk7BAAwAAmIQFGgAAwCRui9lfxxhVxG+Isu3EUvajLOtjKfG0G9m6xuz+olyIIF2u4TpeQrkQHXsNYbWxvXMod7mvWPZPSidaOT4cvcj1e4vuv1WM9j4f1489xLj8zpz3nBRTn4qF+YsbbLwVAAdy6+Cj5C9oAAAAk7BAAwAAmIQFGgAAwCQs0AAAACZhgQYAADAJCzQAAIBJ3Bazv4yxPty+pouR00VM8voQ4pNDvHXsR9gKYD2HY6FcihKvyq3n+hpeX6do+7oba7gt6ViK076+2ja/d4nXKhSstnbYqb2yXIqNT2PkWjcW+7g0fz8JdVbtrencUj/iNanKvSwa/fxmjC/+ZzEmYiZ7OFaU+/wv3ob6tt+aYPnq16HK3hy1VM9PitL/1md1fY/N5/GxOXGnOotzWL78ui5TjssxxrXXx+UcXnKpzjSG3r67vSMHx+WnMdl6Prr9T3NedaB5r3/j/HaML/78mOsd34+hC5/9pJ5PYnvhO6+cT0Z+j1TP6vWz+tnp1PdNpfWh9K2XxPYKDz+v56E1zEOxrT22tknvgjRfNsQtXJr9j3WmNUvjnq6v6vHaecelMv6CBgAAMAkLNAAAgElYoAEAAEzCAg0AAGASFmgAAACTuC3FcYw69SSmOd3cSjup8e6la5XCe1LoT7szwdaVNpMa+UR1wp5eFh75abnsMAHHpNPw/Kdy/K5uquIeuolyH4F1llM/+tnppug2xITHbmPxM3aSm3rnz1U3qXEXW6fJbsxf0AAAACZhgQYAADAJCzQAAIBJWKABAABMwgINAABgEhZoAAAAk7g9Zr8TkdmJXk1Lx6Mj+I+MBU1NpQjYZrl2X7Ze2ncjuPkkTRFjHeaF1L9lj7jm1JfiWV0u4bk6Nx/wFE98ak7qqdwkcchth75bmm3d+zXeUfmcb3xb83yyg+bcllTx6jF2vRvbn6aaGd4dY+Q+Xg7rxTdmir6vdK/XkVuPbNyWv6ABAABMwgINAABgEhZoAAAAk7BAAwAAmIQFGgAAwCQs0AAAACZxe8x+FRUZ4yVvbiVH6ae4yj3iQo9sL12r1I1wLKZ7do9tvdXBxxylv6aL5TeSjjTeS6lM5zEOD1aM0t8jtjz15WN+tu7dkRH24vI3V81DW1/pOD3tcVvTfBKOraGnVbmt6xtj5O+TPbYk6jh6u6jk6G/qju71uuM51tchAADAJCzQAAAAJmGBBgAAMAkLNAAAgElYoAEAAEzCAg0AAGASt8fs896OiuB9Tjuid50kXhV+z7J1RHH1HKTY3PggN8sduqVHd2+Og23dl259abuC0w7XsrWlzUd832azhi0stn5WU3VxO5zmPUj9D3Uu6Sf/qlyYy5f0XDW3YopbEqVY/0vjWm797H8AcTuDWWz9Lktj68C4f39BAwAAmIQFGgAAwCQs0AAAACZhgQYAADAJCzQAAIBJWKABAABM4vaY/SJ+NSaydxJnY1xrqDCV6ybfhsjZdN6tlPrU/5QkGsqlY8l6rhtcL8XJnZoXOZVLUcH3EDkds4ebTp/2byvX80ENpfHVfPZTzHOca9IzEstV/ejVt4t4nbeOLG/OGc3r3+5/p9we9607jx49hirlXPnC/i1jXM9FHVvviNGdT3Z4H6/pIyV+BxbfjtU1fLa+cCh+K/WuSexn2ZHmt+pEWxyt0zzH9aHlksqlj+NioKTtHQ68b5/2Vx4AAMBELNAAAAAmYYEGAAAwCQs0AACASVigAQAATMICDQAAYBI3x+wvVeR5Si4OB8s41Esz2jfEbS4pOjNGudfHYrpyVSzEFi+Pob4q2n6MHAuaTu1aH0snd6r6Ga5jjPtP1z8J13JN8dCpvU6cbqpvTRe5+RvJNdU5ibKPL9/+oHq2unHUE6Ual8q5d+QI6FTuSCmueUnPcbMcd2iPSP8d56HS1lUfPT/t8Fy1nlWPN79xB588W/MXNAAAgElYoAEAAEzCAg0AAGASFmgAAACTsEADAACYhAUaAADAJG6L2V/GuLy+fU23hlaWIhb/+vpcl2lGMl9fhdj4t3V760NdLrVXljnX1/DxOyGCP2w9kCLs14cUNx/O7XWda3p683SDSzi33p4Ez8Tln8M4SVH0KXO96uepPrflnPqfzjttnRCOhfMej2mviaLOJdy3bj/qTjTK/J5qt49uSndjS4xYX9x3ZIfs6K3r7NYXIv3jXb/Uz2osd+8x+7Pctz2kvjTemzMqt6nZ+vSar7L29jVhzo/fX40T32WLjbh90Ht06tb2OlIfZ3qOi3tw9NYonW/tbwqm/XXmnn/9BQ0AAGASFmgAAACTsEADAACYhAUaAADAJCzQAAAAJmGBBgAAMInbYvbHGEsV2xpTpdPBW3uQ5QjPSaJ9YzxpXSz2vntqMW6+PlT1M0XiLyEmNUbpJylKfxYp6vjefyL5UFHaRz3K6Rzaz/EOnY/XuvFsde/dubmlRHt7jknm9K6t+7/H9ejOzfd+b95D2t5mU4138Rgjb5HSFOPOG9NQNz49RrmfumM2tLf1vJ3GTrypc9jjvh2u6ssec16jznv/PAQAAPhoWKABAABMwgINAABgEhZoAAAAk7BAAwAAmMRtKY7r6AWCXcKxItkkpjF2A1ZSudTetZkWtPG1iv0IjcW0oxSCmO5bJzyxm4zTlVL2mmmTnfp2sXV7a7qh523b2sAUIVczpVEFa5Hktlx2GLPdVM89ynF/7uSZetbGpxFTYfd4BHZ45lrpfRuH0z5bJ/NKf05K3813/P7wFzQAAIBJWKABAABMwgINAABgEhZoAAAAk7BAAwAAmIQFGgAAwCRuitlfT8u4fF6s6UKS5TWkdC9FdPzpXYq9D8eKSOkxxrg+1MfWc71WXc8p4zYcqvqS6ktL5hQbn8qFrQDWcAJr2kKgiu4+hY6EY8u1k9s/xjiHwZWi9FN7VT9T/9N5pxzkVGeKgI3XOY2voly4jnHbgdSPssLbi/x/5Tt1pETdor7D462DKi5/KmGsxO0+uvHqVbmt6xtj+kjm3XSf/+6c3tGZh7bIWz/qkUzDMhRrRdvfix0i+OMc1do36SO4/tWjlR7v9DjuMC1sPs5TfbGt5nYxBX9BAwAAmIQFGgAAwCQs0AAAACZhgQYAADAJCzQAAIBJWKABAABM4qaY/WVdc/x9VS61UlS3PIaI9BQ3HyI8Y3JmjO6vD6WY+jKW+VIXiRGkxZYEY+QI2DWV67ZXlUtR1Cl2Od3TFJffjXnu9rNRX+x/aCuWi/1vnFsok7ZwiP0oK7y9yO+XPz2+sI73FOeFbqJuzMVOY6V5j6py6d5dts9Cbo/ZzvYC3Uj8dh/TONm4LzNtBdCdY4/sR+nl/VvSu3xD8T2dyjWf4zV89KQ603ZFVbnlkraMqQ9FlzCPpnJxW5XGeEnjcultEXS01hjf4bnI37HpO6oxX8YtBJpzfeOe+gsaAADAJCzQAAAAJmGBBgAAMAkLNAAAgElYoAEAAEzCAg0AAGASN8Xsj9GLGk0xl9342LK+2L9GXPOL2vv0xGj4EA+dykWteOWDNbdwuAvpvsWs+Rc2u1/V7y/1oZuzv+M1u6mtTrT9GHms73FNOvNGt750TeIeLs32jhwLR6vO7d7epwfdosPnu9Red26oyqVtjMIzEL+90rfGLHPs0Zrz0NbXK26T1W1r6+ucvsvSsfQ52ujjvX8eAgAAfDQs0AAAACZhgQYAADAJCzQAAIBJWKABAABMwgINAABgEsst8ebLsvxkjPHn+3UH+AT8rXVdf9gtbB4CNmAeAj60ch66aYEGAADAfvxfHAEAACZhgQYAADAJCzQAAIBJWKABAABMwgINAABgEhZoAAAAk7BAAwAAmIQFGgAAwCQs0AAAACbx/wB7jE2PgBGt0gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 864x360 with 3 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "approx = outer_product_np(A_A1 @ A_A2 @ A_S2, B_A1 @ B_A2 @ B_S2, C_A1 @ C_A2 @ C_S2)\n",
    "print(\"Relative reconstruction loss:\", np.linalg.norm(np.ndarray.flatten(X_np-approx), 2)  / np.linalg.norm(np.ndarray.flatten(X_np), 2))\n",
    "display_tensor(approx)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAdYAAAFwCAYAAAAFX1qpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAAK90lEQVR4nO3dTYxddR3G8d+ZuUNrCwJtB0srtioUa0EJJCWNRE2M1hgXRiQufNv50piYKFFwYdwYE2NwhQnuWGhiLG40hiawIYopiQGJRCgvpiBV0pFCatuZ6XSO+2ZCe5Pn75k7/XyWNyfP/Jvcw7dnOszt+r4vACBjaugDAMBaIqwAECSsABAkrAAQJKwAECSsABA0GufiLZum+53XzbQ6Sx15ekOz7aqq+TpVi/1C1/SLMKhJf49OOvfY6tX63kiZpHvsZJ2Y6/t+9vzXxwrrzutm6olD1+VOdZ79225ptl1Vdbh/tOk+w5v09+ikc4+tXq3vjZRJusce6Q8eXel13woGgCBhBYAgYQWAIGEFgCBhBYAgYQWAIGEFgCBhBYAgYQWAIGEFgCBhBYAgYQWAIGEFgCBhBYAgYQWAIGEFgCBhBYAgYQWAIGEFgCBhBYAgYQWAoNHQB4CkIy9urk987ivN9kc3nmq2XVU1d1/XdL+qassX55ptd29ON9uGSeGJFQCChBUAgoQVAIKEFQCChBUAgoQVAIKEFQCChBUAgoQVAIKEFQCChBUAgoQVAIKEFQCChBUAgoQVAIKEFQCChBUAgoQVAIJGQx8AgPb66utcvzz0MS6s64Y+wcXrV37ZEysABAkrAAQJKwAECSsABI31w0tHz26sr/1zX6uz1Pynb262XVXVP/bnpvusAqfna+ovzzabv+fvh5ttV1X96D23NN2vquq2b2s3ftLf1cFdAABBwgoAQcIKAEHCCgBBwgoAQcIKAEHCCgBBwgoAQcIKAEHCCgBBwgoAQcIKAEHCCgBBwgoAQcIKAEHCCgBBwgoAQcIKAEHCCgBBwgoAQcIKAEHCCgBBwgoAQaNxLl44N6qXTm5pdZaa3zTdbLuqannUNd1neN3MTE1v29ps/75X9jfbrqqaWv9m0/2qqqXtm5tt9yfG+k8KrEmeWAEgSFgBIEhYASBIWAEgSFgBIEhYASBIWAEgSFgBIEhYASBIWAEgSFgBIEhYASDIb8wGuAQ8//TG+tT2W4c+xgUdOvbk0Ee4aNPXrvy6J1YACBJWAAgSVgAIElYACBJWAAgSVgAIElYACBJWAAgSVgAIElYACBJWAAga73cFHzlbUx97pdFRqg4f+12z7aqqvU8eb7rP8PrFxVr6x9Fm+0sf7ZptV1X98uU/Nt2vqvrSje3+Pt2dWWi2DZPCEysABAkrAAQJKwAECSsABAkrAAQJKwAECSsABAkrAAQJKwAECSsABAkrAAQJKwAECSsABAkrAAQJKwAECSsABAkrAAQJKwAECSsABAkrAAQJKwAECSsABAkrAASNxrm4v2JDLd1+W6uz1PdfW262XVX16tJc031WianpZtPdB9/XbLuq6q5n39V0v6pq3fuvaLbdP3NZs22YFJ5YASBIWAEgSFgBIEhYASBIWAEgSFgBIGis/90GgMm0NLux5u7cN/QxLmjvkzcMfYQx/HjFVz2xAkCQsAJAkLACQJCwAkCQsAJAkLACQJCwAkCQsAJAkLACQJCwAkCQsAJAkLACQJCwAkCQsAJAkLACQJCwAkDQWB903o+6mt/U7rPRXzg122y7qmr+nM91X+uWZjfW8bv2Ntt/x4N/bbZdVXXZx0833a+qevngTc22F7/bbBomhidWAAgSVgAIElYACBJWAAgSVgAIElYACBJWAAgSVgAIElYACBJWAAgSVgAIElYACBJWAAgSVgAIElYACBJWAAgSVgAIElYACBJWAAgSVgAIElYACBJWAAgSVgAIGo1z8R1799TDv7671Vma65574NDQZ6CtD9+6px6+/9vtvsD97abXgu7OH7rHuOR1fd9f/MVd93BVbWl3nObm+r7/5NCHoJ018B6ddO6xVarruuNVdXToc6wxO/q+nz3/xbHCCgC8Nf/GCgBBwgoAQcIKAEHCCgBBwgoAQcIKAEHCCgBBwgoAQcIKAEHCCgBBwgoAQWN9us305Rv70earW52l1r18utl2VdV8narFfqFr+kUY1OVXz/Sbt69vtn/imZlm22uBe2z1Gr1tYz9z5aahj3FBM6+dGvoIF+1knZhb6ZfwjxXW0eara+u938qd6jy7DjzRbLuq6nD/aNN9hrd5+/q696Fbm+0/tPuaZttrgXts9Zq5clNd/4WGH6kYsvVnjw99hIv2SH9wxU8L8q1gAAgSVgAIElYACBJWAAgSVgAIElYACBJWAAgSVgAIElYACBJWAAgSVgAIElYACBJWAAgSVgAIElYACBJWAAgSVgAIElYACBJWAAgSVgAIElYACBqNc/HNV83VE5/5Rauz1P4DtzTb5tJw4ujb6zdf399sf2bP6WbbVVWX/fyNpvtVVYufPdtsu3tjutk2TApPrAAQJKwAECSsABAkrAAQJKwAECSsABAkrAAQJKwAECSsABAkrAAQJKwAECSsABAkrAAQJKwAECSsABAkrAAQJKwAEDQa+gAA/H8sTw99gkuDJ1YACBJWAAgSVgAIElYACBJWAAga66eC/3Z8tnY/cKDVWWrnB15vtl1V1R35U9N9VoH/nqmpx55qNv+vb+xrtl1Vdc1Hnmu6X1U1tWFDu/Hlvt02TAhPrAAQJKwAECSsABAkrAAQJKwAECSsABAkrAAQJKwAECSsABAkrAAQJKwAECSsABAkrAAQJKwAECSsABAkrAAQJKwAECSsABAkrAAQJKwAECSsABAkrAAQJKwAEDQa5+J1r5+td//q363OUnX8P+22q6oWz7bdZ3Dd+nU1ff2uZvsHv/eTZttVVd98cH/T/aqqhQ/tbra9/OdHm23DpPDECgBBwgoAQcIKAEHCCgBBwgoAQcIKAEHCCgBBwgoAQcIKAEHCCgBBwgoAQWP9rmAAJtPMa6dq208fH/oYF3To2FNDH+GiTV+78uueWAEgSFgBIEhYASBIWAEgSFgBIEhYASBIWAEgSFgBIEhYASBIWAEgSFgBIEhYASBorF/C3y8s1rnnX2p1lnrjy/uabVdVLf3+UNN9htfPL9S5Z55rtn9gxx3NtquqXv3tzqb7VVXv/PzTzba7s2eabcOk8MQKAEHCCgBBwgoAQcIKAEHCCgBBwgoAQcIKAEHCCgBBwgoAQcIKAEHCCgBBwgoAQcIKAEHCCgBBwgoAQcIKAEHCCgBBwgoAQcIKAEHCCgBBwgoAQcIKAEHCCgBBo3Eu7mZGNZrd2uosNb3YN9uuqurazrNadF2z6eld7222XVV1am5D0/2qqqkbdjbb7l68rNk2TApPrAAQJKwAECSsABAkrAAQJKwAEDTWTwUDMKE2rK9uz01Dn+KCbr9n39BHGMN3VnzVEysABAkrAAQJKwAECSsABAkrAAQJKwAECSsABAkrAAQJKwAECSsABAkrAAQJKwAECSsABAkrAAQJKwAECSsABAkrAASNxrt6VMubr2pzkqo6s6Vt55fH+9MygRZ2bKgjP7it2f7uu59vtl1VteurLzTdr6r6w7Gnmm3v3X+i2TZMCk+sABAkrAAQJKwAECSsABAkrAAQJKwAECSsABAkrAAQJKwAECSsABAkrAAQJKwAECSsABAkrAAQJKwAECSsABAkrAAQJKwAECSsABAkrAAQJKwAECSsABAkrAAQ1PV9f/EXd93xqjra7jjN7ej7fnboQ9DOGniPTjr32Crl3mhixff7WGEFAN6abwUDQJCwAkCQsAJAkLACQJCwAkCQsAJAkLACQJCwAkCQsAJA0P8A2sJDwpUU6P0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 576x360 with 6 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "display_topic_matrices([A_S2.T, A_S1.T, A_A2], S_trues)"
   ]
  }
 ],
 "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": 4
}