Tensor1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Intro\n",
    "\n",
    "* Here we compare an PARAFAC decomposition when all the factors are unknown with when one factor is known beforehand.\n",
    "* We then compare the use of two Unfolding fuctions, one from tensorly and one adapted from another library"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Using numpy backend.\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "import tensorly as tl\n",
    "from scipy import linalg\n",
    "import tensorly.decomposition as tldecomp\n",
    "import tensorly.metrics as tlmetrics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def tensorEye(I1,I2,I3):\n",
    "    out=np.zeros((I1,I2,I3)) \n",
    "    for i1 in range(0,I1):\n",
    "        for i2 in range (0,I2):\n",
    "            for i3 in range(0,I3):\n",
    "                if (i1==i2==i3): out[i1,i2,i3]=1\n",
    "    return out"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def from_to_without(frm, to, without, step=1, skip=1, reverse=False, separate=False):\n",
    "    if reverse:\n",
    "        frm, to = (to - 1), (frm - 1)\n",
    "        step *= -1\n",
    "        skip *= -1\n",
    "    a = list(range(frm, without, step))\n",
    "    b = list(range(without + skip, to, step))\n",
    "    if separate:\n",
    "        return a, b\n",
    "    else:\n",
    "        return a + b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def unfold(tensor, mode):\n",
    "    sz = np.array(tensor.shape)\n",
    "    N = len(sz)\n",
    "    order = ([mode], from_to_without(N - 1, -1, mode, step=-1, skip=-1))\n",
    "    newsz = (sz[order[0]][0], np.prod(sz[order[1]]))\n",
    "    arr = np.transpose(tensor,axes=(order[0] + order[1]))\n",
    "    arr = arr.reshape(newsz)\n",
    "    return arr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "I1=3\n",
    "I2=3\n",
    "I3=3\n",
    "Rank=3\n",
    "Ts=tl.random.cp_tensor(shape=(I1,I2,I3),rank=Rank,full=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[[ 1.  0.  0.]\n",
      "  [ 0.  0.  0.]\n",
      "  [ 0.  0.  0.]]\n",
      "\n",
      " [[ 0.  0.  0.]\n",
      "  [ 0.  1.  0.]\n",
      "  [ 0.  0.  0.]]\n",
      "\n",
      " [[ 0.  0.  0.]\n",
      "  [ 0.  0.  0.]\n",
      "  [ 0.  0.  1.]]]\n"
     ]
    }
   ],
   "source": [
    "eye=tensorEye(Rank,Rank,Rank)\n",
    "print(eye)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## ALS with all factors Unknown"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "itermax=1000\n",
    "epsilon1=[]\n",
    "A2=np.random.rand(I2,Rank)\n",
    "A3=np.random.rand(I3,Rank)\n",
    "for ii in range(itermax):\n",
    "    krtransposed=linalg.pinv(tl.transpose(tl.tenalg.khatri_rao([A3,A2])))\n",
    "    A1=np.dot(unfold(tensor=Ts,mode=0),krtransposed)\n",
    "    krtransposed=linalg.pinv(tl.transpose(tl.tenalg.khatri_rao([A3,A1])))\n",
    "    A2=np.dot(unfold(tensor=Ts,mode=1),krtransposed)\n",
    "    krtransposed=linalg.pinv(tl.transpose(tl.tenalg.khatri_rao([A2,A1])))\n",
    "    A3=np.dot(unfold(tensor=Ts,mode=2),krtransposed)\n",
    "    error=tlmetrics.MSE(y_pred=tl.tenalg.multi_mode_dot(eye,matrix_or_vec_list=[A1, A2, A3]),y_true=Ts)\n",
    "    epsilon1.append(error)\n",
    "Xhat1=tl.tenalg.multi_mode_dot(eye,matrix_or_vec_list=[A1, A2, A3])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[[ 0.32274291,  0.37458317,  0.42271877],\n",
       "        [ 0.38404914,  0.45381235,  0.44819282],\n",
       "        [ 0.26943822,  0.33289278,  0.27810147]],\n",
       "\n",
       "       [[ 0.1935513 ,  0.16647689,  0.22589232],\n",
       "        [ 0.2106454 ,  0.18710368,  0.22454888],\n",
       "        [ 0.12795291,  0.12008356,  0.11856905]],\n",
       "\n",
       "       [[ 0.13253052,  0.15161259,  0.12736588],\n",
       "        [ 0.20590848,  0.23242478,  0.14000027],\n",
       "        [ 0.18144013,  0.20551398,  0.09138585]]])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Xhat1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[[ 0.32209341,  0.37541617,  0.42216749],\n",
       "        [ 0.38393694,  0.45382583,  0.44820815],\n",
       "        [ 0.27053073,  0.3317248 ,  0.27873834]],\n",
       "\n",
       "       [[ 0.19532735,  0.16447497,  0.22632321],\n",
       "        [ 0.2106079 ,  0.18722507,  0.22460432],\n",
       "        [ 0.12513622,  0.12310625,  0.11807863]],\n",
       "\n",
       "       [[ 0.13211993,  0.15192751,  0.12762588],\n",
       "        [ 0.20579037,  0.23259446,  0.13993744],\n",
       "        [ 0.18177979,  0.20517857,  0.0910915 ]]])"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Ts"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## ALS with one known factor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "itermax=1000\n",
    "epsilon2=[]\n",
    "A2=np.random.rand(I2,Rank)\n",
    "A3=np.random.rand(I3,Rank)\n",
    "for ii in range(itermax):\n",
    "    # krtransposto=linalg.pinv(tl.transpose(tl.tenalg.khatri_rao([A3,A2])))\n",
    "    # A1=np.dot(unfold(tensor=Ts,mode=0),krtransposto)\n",
    "    krtransposed=linalg.pinv(tl.transpose(tl.tenalg.khatri_rao([A3,A1])))\n",
    "    A2=np.dot(unfold(tensor=Ts,mode=1),krtransposed)\n",
    "    krtransposed=linalg.pinv(tl.transpose(tl.tenalg.khatri_rao([A2,A1])))\n",
    "    A3=np.dot(unfold(tensor=Ts,mode=2),krtransposed)\n",
    "    error=tlmetrics.MSE(y_pred=tl.tenalg.multi_mode_dot(eye,matrix_or_vec_list=[A1, A2, A3]),y_true=Ts)\n",
    "    epsilon2.append(error)\n",
    "Xhat2=tl.tenalg.multi_mode_dot(eye,matrix_or_vec_list=[A1, A2, A3])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[[ 0.3227448 ,  0.37458555,  0.42272614],\n",
       "        [ 0.38405194,  0.45381592,  0.44819393],\n",
       "        [ 0.26943403,  0.33288708,  0.27808707]],\n",
       "\n",
       "       [[ 0.19355341,  0.16647762,  0.22589438],\n",
       "        [ 0.21064507,  0.18710322,  0.22454573],\n",
       "        [ 0.12794715,  0.12007915,  0.11855877]],\n",
       "\n",
       "       [[ 0.1325217 ,  0.15160373,  0.1273701 ],\n",
       "        [ 0.20590868,  0.23242533,  0.14000729],\n",
       "        [ 0.18144559,  0.20551926,  0.09138943]]])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Xhat2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[[ 0.32209341,  0.37541617,  0.42216749],\n",
       "        [ 0.38393694,  0.45382583,  0.44820815],\n",
       "        [ 0.27053073,  0.3317248 ,  0.27873834]],\n",
       "\n",
       "       [[ 0.19532735,  0.16447497,  0.22632321],\n",
       "        [ 0.2106079 ,  0.18722507,  0.22460432],\n",
       "        [ 0.12513622,  0.12310625,  0.11807863]],\n",
       "\n",
       "       [[ 0.13211993,  0.15192751,  0.12762588],\n",
       "        [ 0.20579037,  0.23259446,  0.13993744],\n",
       "        [ 0.18177979,  0.20517857,  0.0910915 ]]])"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Ts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmAAAAFuCAYAAAAmiaQ9AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xd4FOXaBvB7tibZ3WRTaEFCESJFgxRBinQELAhIIqABkUNTUAE5lCNNEARU5MjhO9gQEUEitiMWEBBFBakiSFGkSQkB0naTbJv5/kiykCVLKJmS5P5d17nYnZndeZLX5Nx55t13BEmSJBARERGRYnRqF0BERERU0TCAERERESmMAYyIiIhIYQxgRERERApjACMiIiJSGAMYERERkcIMahdwPdLSshU5j9VqhsPhUuRcdG04JtrEcdEejok2cVy0Se5xqVTJFnQfO2DFMBj0apdAATgm2sRx0R6OiTZxXLRJzXFhACMiIiJSmGKXIHNycjBu3DhkZmYiNDQU8+fPR1RUVJFjVq9ejVWrVsFgMGDkyJHo2LGjUuURERERKUaxDtjq1avRqFEjfPDBB7j//vuxePHiIvvT0tKwfPlyrFq1Cm+//TZeffVVuN1upcojIiIiUoxiHbDHH38cPp8PAHD69GnExMQU2b937140adIEJpMJJpMJcXFxOHjwIBISEpQqkYiIiEgRsgSwlJQULFu2rMi22bNnIyEhAQMHDsThw4exdOnSIvsdDgdstkufFrBYLHA4HEWOsVrNikyY0+t1sNvDZD8PXTuOiTZxXLSHY6JNHBdtUnNcZAlgiYmJSExMLHbfe++9hyNHjmD48OH49ttv/dutViucTqf/udPpLBLIACj2EV67PQwZGTmKnIuuDcdEmzgu2sMx0SaOizbJPS6aWIZiyZIl+PTTTwEAYWFh0OuLdrISEhKwc+dOuFwuZGdn48iRI4iPj1eqPCIiIiLFKBbAHn74Yfzvf/9DcnIyxo0bh9mzZwMAli5dig0bNqBSpUpITk7GgAEDMGjQIIwZMwZms1mp8oiIiOg67d+/D6NGDSt2X8+e3RSupqglS/6DIUOSsWvXjmt+zdmzZ7Fly/cyVnWJYpPwY2Ji8Pbbb1+xffDgwf7HSUlJSEpKUqokIiIiukErVizDN998iZCQULVLKdbGjeuxdOkKhIVZrvk1u3Ztx/Hjx9C2bTsZK8tXpm5FRERERNpQvfotePHF+Zg5c+pVj1uy5D9wOBwYO/af6N+/D+64ozFOnDiOqKgozJo1D5IkYc6cGTh16hR8Ph/69XsUcXE18eab/4d5817D+vVf4/33l2HZspX49dc9+PrrtYiJicGZM6eRnp6O1NQzGD16LFq2bOU/59KlbyIt7RzGj38WL7/8byxc+DLOnUtFZmYm7r67NYYOHYmTJ0/g2WdnIzfXhZCQEEybNgvvv/8u8vLycMcdCahcuQoWLJgPvV4Pk8mEf/7zeUiSiAkTxiA8PAKtWrXBo48OuuHvHwMYERFRGbZ2fyo+33e2VN+z5+1VcX+jKlc9pkOHzjhz5vRVj1m06DXodALGjZsAADh9+hQWLvw/VKlSFSNHPoEDB37HoUO/IyLCjilTZiInx4knnngM//3vUpw9ewYulwvbtv0MQRBw8eIF/PjjZrRv3xH79/8Go9GEV175N7Zv34qVK1cUCWCDBw/F2rWf49VXF+HixQto1OgOTJw4BS6XC3363IehQ0fiP/95Df/4x1DcfnszbNiwHn/++Qcee+zxgg5YewwZkoyJE59HvXq34YcfvsOiRa/iqaeexcWLF/D22+/DaDTe1PeYtyIiIiKiUnfx4gUcOfIHcnJy/dsiIuyoUqUqAKBy5Spwu104duwYGjduCgAIC7OgVq3aOHXqb7Ro0Qq7d+/EuXOpuPfe7tix4xfs2bMbzZu3AADEx99W8D5V4XYHXyUhPDwcBw7sx4wZz+P111+F2+0BAJw4cRyNG98JAOjcuStatLi7yOvOn09DvXr552jcuCmOHv0LAFCtWuxNhy+AHbArLN5yFO0bVEGjaK7XQkRE2nd/oyoldqvUEBUVjVdfXYTRo4dj69afcPfdrSEIwhXH1apVC3v37kb79h2Rk+PEkSNHEBsbi3btOuCNNxajXr3b0KJFK8yfPxs1atSAwZAfXYp5q2J9+eUXsFpt+Oc//4W//z6Jzz//BJIkoWbN2ti37zc0aHAn1q37CllZmbBYrJAkEQAQE1MJf/75B+rWrYc9e3ahRo24gvOWTu+KASxAyp7TEHU6NGpdU+1SiIiIyjRBEDBp0lSMHTsKb7zxbrHH9OzZB3PnzsLIkUPgcrnwxBNDERkZBbs9EidPHsejjw5E3br1cPbsGQwYMPC6a2jW7C5Mnz4Ze/fuQUhICG65pQbOn0/DU089gwULXoLL5UFISAimTp2Js2fP4L333kF8fH1MmPAvLFiQP0dNr9dj4sQpN/ndKEqQJEkq1XeUUVpatuzn6LjoR/RpcgtGt2EA0xIuYqhNHBft4ZhoE8dFmyrEQqxlxXJMRYu0NWqXQUREROUYL0EGqIXTcLuOql0GERERlWPsgAVwCwboRY/aZRAREVE5xgAWwAsDDBIDGBEREcmHASyAG0Z2wIiIiEhWDGABPDBAzw4YERERyYgBLIAbRgYwIiKiGzBq1DAcP34Mb7+9BJ9++tENv8+uXTswbdqkUqxMexjAAnjZASMiIiKZcRmKAG4YYOAcMCIioqCcTgdeemkWHI5sZGZm4MEHe6N3775Xfc3bby9BdHQ0evXqi+PHj2H+/NlYtOgNDBrUD3fe2RRHjvwJAHjppVf9r8nLy8PkyePRvft9iImphBUr3oPRaMCZM6fRqVNXDBo0BGfOnMZLL82E1+uFIAh45pnnsHv3Dni9PgwYkIx5816EyWTGs88+h3fffQuxsbfg888/Rr16t+HkyaPIzMzCzJlzUbVqNVm/Z4EYwALkzwFzq10GERHRNTEf/AghB1aV6nvmNegHV/3ggervv/9Gly73on37Tjh/Pg2jRg0rMYAF43Q60aVLN4wZ80/MmPE8tm79EVFR0cjNzcWECWOQmNgPbdu2x65dO5CaegbvvrsSHo8HvXp1x6BBQ/Cf/7yGvn0fwT33dMAffxzCSy/NxIsvzsecOS9gwIBknDx5Anl5eQCAX37ZivnzX8Pnn3+MBg0aYdq0KZg7dz7Wr/8GycmP31D9N4oBLIAHRuglp9plEBERaVZ0dDRWr/4AmzdvQliYBV6v97peH3gXxPj42wAAlStXgdud3wTZvXsXbr21LtzuS1el6tSpC4PBAIPBALM5BABw7NgxNG7cFABQr95tOHcuFVWrVoXLlYfff9+HmjVrIzX1DA4c2A+r1QqLxVrknFWqVMGFCxdu4LtwcxjAAvBTkEREVJa46ve9ardKDitXLsfttyegd+++2LVrB37+eUuJrzGZzP6gc/jwwYC9whXHt27dBs888xyefHIoEhIa5x915WGoVasW9u7djbZt2+OPPw4hKioaANCqVRssXvxvJCUNQGrqWSxYMB89e/a6dMbi3kxBnIQfwAMjF2IlIiK6ijZt2iElZSVGjhyC1as/gF6v93euguncuSt+/vlHjB49HIcPH7qm80RFRWPIkGGYPXtG0GOeeupZfPTRh3jqqaF4+eWXMGnSFABA+/ad8Ntvv6JZs+Zo2fJuHDp0AG3bdrjmr1FughTYB9SwtLRs2c/x03+S0d6wD57hu2U/F107ue9YTzeG46I9HBNt4rhok9zjUqmSLeg+dsACeAR2wIiIiEheDGAB8ueAXd9kQiIiIqLrwQAWIH8OGJehICIiIvkwgAVgB4yIiIjkxgAWwCMYoIMIiAxhREREJA8GsACewqXRfJyIT0RERPJgAAvggREAIPhcKldCRERE5RUDWACx8FsiieoWQkREROUWA1gAnz+A+dQthIiIiMotBrAAEvQAAIEBjIiIiGTCABbAJxR8S0RegiQiIiJ5MIAFEHkJkoiIiGTGABaAAYyIiIjkxgAWoDCACfwUJBEREcmEASyAvwMmsgNGRERE8mAACyAKvARJRERE8mIACyAWLEPBhViJiIhILgxgAS7NAWMHjIiIiOTBABaAc8CIiIhIbgxgAbgMBREREcmNASzApUn4nANGRERE8mAAC8A5YERERCQ3BrAAEoSCB+yAERERkTwMSp0oJycH48aNQ2ZmJkJDQzF//nxERUUVOWbWrFnYtWsXLBYLAGDx4sWw2WxKlQgAkISCZSh4M24iIiKSiWIBbPXq1WjUqBFGjRqFjz/+GIsXL8bzzz9f5Jj9+/fjrbfeuiKYKcnHSfhEREQkM8UC2OOPPw6fLz/UnD59GjExMUX2i6KI48ePY+rUqTh//jz69u2Lvn37KlXepToEzgEjIiIieckSwFJSUrBs2bIi22bPno2EhAQMHDgQhw8fxtKlS4vsz8nJwWOPPYbBgwfD5/Nh4MCBuP3221G/fn05Sgzq0jIUvARJRERE8hAkSZKUPumRI0cwfPhwfPvtt/5tPp8Pubm5sFqtAIB58+YhPj4evXr18h+Tm+uGwaCXtbZRr76LJblj4e37PqTb7pP1XHTt9HodfD6GYq3huGgPx0SbOC7aJPe4GI3BM4tilyCXLFmCKlWqoFevXggLC4NeX7SoY8eOYcyYMfjkk08giiJ27dqF3r17FznG4XDJXqenYBycjhy4M3JkPx9dG7s9DBkcD83huGgPx0SbOC7aJPe4VKoU/IOEigWwhx9+GBMmTMCaNWvg8/kwe/ZsAMDSpUsRFxeHzp0748EHH0RSUhKMRiMeeugh1KtXT6ny/CQuxEpEREQyUyyAxcTE4O23375i++DBg/2Phw4diqFDhypVUrF8XIiViIiIZMaFWANIKFwHjAGMiIiI5MEAFoD3giQiIiK5MYAF4EKsREREJDcGsAAS54ARERGRzBjAAojgvSCJiIhIXgxgAS4tQ8EOGBEREcmDASyAjwGMiIiIZMYAdoXCOWC8BElERETyYAALwGUoiIiISG4MYAHEwm8JF2IlIiIimTCABRC5DhgRERHJjAEsgMg5YERERCQzBrAAEjtgREREJDMGsACcA0ZERERyYwALIOh0+SGMlyCJiIhIJgxgAQTkd8F4L0giIiKSCwNYMURBxzlgREREJBsGsACCUDAPjDfjJiIiIpkwgBUjfw4YO2BEREQkDwawYoiCnnPAiIiISDYMYMXwwgCIXrXLICIionKKASyAIAjwCCYIPpfapRAREVE5xQBWDI9gBLwMYERERCQPBrAAAsAOGBEREcmKASyAIBQGsDy1SyEiIqJyigGsGB6YeAmSiIiIZMMAFiD/EqSRlyCJiIhINgxgxXALJgjsgBEREZFMGMCukL8MBdgBIyIiIpkwgBXDA16CJCIiIvkwgAXwfwqSlyCJiIhIJgxgAQrXAeMlSCIiIpILA1gx3PwUJBEREcmIASyAIOSvAyb4XIAkqV0OERERlUMMYMVwC8b8B+yCERERkQwYwAIIKFgJH+BlSCIiIpKFQe0CtChdHwkAMB1dB290w/zrkoIOEPQF/wqQAp5D0EHyP9cBEACdHhIKnut0QOHjwtcQERFRhcQAFkgQsCOkFcSQSIRvGCPbaSQIl4UxXUGA012xDdBB0l16fLVji329oIdUGCChA3T6/ACJS/svBczCEFk0cErFBdDC94KuyPGSP5Dqr71WXAqoUjFhVRJ0EKyhMOV4itQsCVd+Df4QXPj99YfgawnR+uJrZYgmIqJSxgBWDBfMSE/6GvqLhyH43ABEQPRBkCT/Y0DMn6QviRAkHyBdeg7JB8H/vHCfmP/6y54DEgT/e4mAKPofC4XHXPa/ott8l85fWFNhff6apMuOLdjmcwXUF1i/79K5IF32tebvFy6vv9i6Crah9D/AEFHq73hzioRonR6SYAB0BkBnhKQv+FdnAAq2S3pj/r8Fx1x6bIBU8Nz/mmKPK9hnCIGkN0PSmy89NpiBgm2SIeSyxwXbDWZAZ2JoJCLSCAawAAIACYBoqw7RVl3tcsouf/grCJWieCksFhvgigbGywMiJB9sVjOys3IuC5uXAqJwWUAsEiADA+/Vgm1hjYUBVvTlv39A4BWKPL8slIre/LpFD+DzAqI3/7HohSB5C7Z5IIgF/3rzCp77Lh1XsO/Say/f5ymdYdGbIRktkExWSEYrJJMVov+5BZLJdmm/yQYxNApSSBTE0GiIIVGQQuwFHUYiIroZDGAB2CAoJYWXIaHPf17wT3F9sWvqldnD4DPnlE5tZVFhmPR5IPhc+R8Q8bkgeC9/nFdkOwqOE7x5Bc/z8h97ciB4HBDcTggeB3R5FyFkHb/03OMMXoagg2S2QwyNgWirBl10LYSZqsBni4VoqwFvZF1IodEKfmOIiMomBrBicPkv0hxBuHQp0xgqwwXey0hifkhzZeWHs7yL0OVegC73AoS89PzHOWnQOc5Ad2g/LDnni7xcDI2GNyoevqjb4KnSBJ5qd0G01eBfN0REl2EAC8D/i6AKT9AVXIK0QrTFXvVQuz0MGWkXoHeegS7rBAwX/4D+4iEYLh6G+WAKQn97FwDgs1SBO64D3LXuhbtGO8AYqsAXQkSkXQxgRHRzjKHw2evAZ68DT1yHS9tFH/QXD8F4dgeMp7bCfOQrhB74EKLRAld8b+Q1egzeSrerVjYRkZoUn0175MgRNGvWDC7XlYucrl69Gn369EFSUhI2bdqkdGn5eJmEqHTo9PDFNETe7QOR3W0xLjzxKzJ6roL71vsQcjAFkau7I+KzfjCc3al2pUREilO0A+ZwODB37lyYTKYr9qWlpWH58uVYs2YNXC4XBgwYgDZt2hR7rNwkTgIjKn16Izw12sJToy0cbaYh5MCHCNv9f4hc8xDy6vaEo+10SJbKaldJRKQIxTpgkiRhypQpGDt2LEJDr5z/sXfvXjRp0gQmkwk2mw1xcXE4ePCgUuX5FS5DQUTykULsyG0yHBeSf4KzxTiY//oaUSs7wvzHZ2qXRkSkCFk6YCkpKVi2bFmRbbGxsbjvvvtQv379Yl/jcDhgs9n8zy0WCxwOR5FjrFYzDAZ96Rd8GaNRD9Ejwm4Pk/U8dH30eh3HRINuflzCgK7/grdpEvRfjEb4uqfgu7ALYpdZgMFcanVWJPxZ0SaOizapOS6yBLDExEQkJiYW2da1a1esWbMGa9asQVpaGp544gmsWLHCv99qtcLpvLT+kNPpLBLIAMDhkP/m2F6PD6IEZGRU4DWnNMhuD+OYaFCpjYu+OvDAh7BsfQlhO5fAd2Yfsu57B5JZa/c/0D7+rGgTx0Wb5B6XSpVsQfcpNgds/fr1/sedOnXCO++8U2R/QkICXnvtNbhcLrjdbhw5cgTx8fFKlecn8BokkTr0RjjbTIG30h2wbRgD+8d9kPng+xCt1dSujIio1Km+DMXSpUsRFxeHzp07Izk5GQMGDIAkSRgzZgzMZnUuQUhMYESqccX3ghgag/Cv/oGIzx5BRq+PODmfiModQSpDH/lLS8uW/RyjP/oNeaKIN5May34uunZs32uTnONiOLMD9s8HwBdeAxm9UiCFRslynvKGPyvaxHHRJjUvQfKuusUoO5GUqPzyVmuOzPuXQp95DBGf94fglv8PMCIipTCABRI4BYxIKzy3tEFWjzdhuHAQ4d+MBESv2iUREZUKBrAAXAefSFvcNTvB0X42TCe+g/WHqWxRE1G5oPokfE3i73ciTclr9Cj0mUcRtvu/8NlvRW7jIWqXRER0U9gBC8BbQRJpk7PVZLhqd4PlxxdgPL1N7XKIiG4KA1gxuAwFkQYJOmR3XgBfeBxs34yEzpmqdkVERDeMASyAwFlgRJolmcOR1eNN6NzZ+ZPyfR61SyIiuiEMYMXgHF8i7fJF10d2x/kwnvkFlp/nqF0OEdENYQALwDlgRNrniu+F3NsHIezXN2A8vkntcoiIrhsDWDHYACPSPkeb5+GNro/wDWMgOM+pXQ4R0XVhACsGL0ESlQGGUGTduxiCOxvhG8YAkqh2RURE14wBLACvQBKVHb6oeDjaTofp5GaE7nlT7XKIiK5Z0AB27lzxLf1ff/1VtmK0ogzdn5yowstr9Fj++mBbX4Lh3F61yyEiuiZBA9hzzz3nfzx+/Hj/41deeUXeilQmcBY+UdkiCMju9DLEsBjY1o8CPLlqV0REVKKgAezyLtDZs2eL3V5elf+vkKh8kUIikd1pAQwZf8H684tql0NEVKLrngNW3jtE5furIyq/PDXaIidhCEJ/exfGk9+rXQ4R0VUFDWCXB63yHrquwBYYUZnkbDUR3sh6sG0YCyEvQ+1yiIiCMgTbsWvXLrRt2xYAkJGR4X+cmZmpTGUqEQTeC5KozDKEIrvLQtjX9IT1++eRfe8itSsiIipW0AC2b98+JesgIioV3soJyGn+LCy/vAx37W5w1XtQ7ZKIiK4Q9BKk2+3GsmXLIEkSUlNT8fTTT+O5555DWlqakvWpogJ8zoCoXMtpNgqeKk1g3TwROufZkl9ARKSwoAFs5syZOH36NERRxPTp01G/fn1069YN06dPV7A85VW4+W5E5ZHOgOwuCyH4XLBtHMe/qohIc4IGsNOnT2PSpEnwer3YuXMnhg4diq5du+LixYtK1qcK/qomKvt89jpwtJ4C04nNCNm/XO1yiIiKKPFTkLt27cIdd9wBo9EIAHC5XMpUphL2v4jKj7zbB8Id1x7WH2dCn/GX2uUQEfkFDWBhYWH48MMP8eabb+L++++HKIpISUlBtWrVlKxPFbxaQVROCAKyO74MSW+C7dtnANGrdkVERACuEsCmT5+OEydOoHPnzujduze2bduGjRs3Ytq0aUrWpzhOASMqX0RrNTjaz4YxdTfCdnJZCiLShqDLUERFRRW5B2SrVq3QqlUrRYpSG9cBIypfXPUeQt7RdQjb8RrcNTvCW7mx2iURUQUXNIB16tTpik8ESpIEQRCwYcMG2QtTiwBegiQqjxztXoTx9DbYvn0G6YlfAcZQtUsiogrsqgFs3759aN26NXr27InY2Fgl61IRr0ESlUdSiB3ZnRfA/nl/WH6eDWe7mWqXREQVWNAA9vzzz0MURWzZsgWLFy9GZmYmunTpgh49esBkMilZIxFRqfDUuAc5CUMQtvdtuGt3hadGO7VLIqIKKugkfADQ6XRo164d5s2bh7lz5+LHH39E69atlapNFZyET1S+Xbph9xgIeelql0NEFdRVA5goivjhhx8wceJEDBo0CDVq1MBHH32kVG2qkTgJjKj8MoQiu+u/ocu9AOv3z6tdDRFVUEEvQc6YMQPbt29HixYtkJSUhKZNmypZl2rYACMq/7yV7kDOXeNg2TYX7lpd4YrvpXZJRFTBBA1gK1euhN1ux7p167Bu3boi+7Zs2SJ7YWpi/4uo/MtpOhKm4xtg/f5f8FRrAdFWUT5oRERaEDSAHTx4UMk6NEMQuAwFUYWgMyCry2uIWnUvbBvHIrPnB4Bw1VkZRESlhr9tiKjCEiNqwXHPdJj+3oLQve+oXQ4RVSAMYMVgA4yo4shr0B+uWl1h+XkO9BcPq10OEVUQ1xzAMjMz5axDMwJX/yeick4QkN1xHiSTFbb1TwM+t9oVEVEFUGIA++WXX/DAAw+gX79+WLhwIVJSUpSoS1VchoKoYpHCKiG7wzwYz+9D2PYFapdDRBVAiQFs4cKFeP/99xETE4MRI0Zg5cqVStSlGva/iComd51uyG3QD2G7/gPDme1ql0NE5VyJAUyn08Fut0MQBJjNZlgsFiXqUhX7X0QVk7PtdIi2GghfPxqCq2JMuyAidZQYwOLi4vDKK68gIyMDb7zxRrm/KbcggAmMqIKSTFZkdX0dOudZWDdN4Jo0RCSbEgPYjBkzEBsbi2bNmiEsLAyzZs1Soi7V8BIkUcXmrdoUzpbjEXLkC4QcKN9TLohIPSUGsNzcXERGRqJx48awWq1Yu3atEnWpSmILjKhCy20yEu5b7oH1h6lcmoKIZBF0JfxCTz31FKpXr46YmBgAFWCZhvL+9RFRyQQdsru8hshVXRG+7kmk9/0CMISoXRURlSMlBjBJkjBnzpxSO+GRI0eQlJSEn376CWazuci+WbNmYdeuXf6J/osXL4bNZiu1c18rTvsgItFSBdmdFyBi7SBYf5oFR7vyPf2CiJQVNIC53fmLEdaoUQO7d+9Go0aN/PtMJtMNnczhcGDu3LlBX79//3689dZbiIqKuqH3Lw3sfxFRIXetzshp/A+E/foW3DXawV37XrVLIqJyImgA6969OwRBgCRJ2Lp1q3+7IAjYsGHDdZ9IkiRMmTIFY8eOxZNPPnnFflEUcfz4cUydOhXnz59H37590bdv3+s+T2lgA4yICjlbTYLx1FbYNoxF+iPfQLRVV7skIioHggawjRs3AgD27t2LhIQE//Zt27aV+KYpKSlYtmxZkW2xsbG47777UL9+/WJfk5OTg8ceewyDBw+Gz+fDwIEDcfvttxc53mo1w2DQl3j+m2E2539L7PYwWc9D10ev13FMNKhijEsY0PcdCO90QuS3I+FLXgsYzCW/TCUVY0zKHo6LNqk5LoIU5L47O3bswJEjR7B06VIMHjwYQH6XasWKFfjiiy+u+0Rdu3ZF1apVAQB79uxBQkICVqxY4d/v8/mQm5sLq9UKAJg3bx7i4+PRq1cv/zFpadnXfd7r9cLXh7D9ZCb+N7SF7Oeia2e3hyEjI0ftMihARRoX05G1iPh6OHLvGARHuxfVLieoijQmZQnHRZvkHpdKlYLPYw/aAQsPD0daWhrcbjfS0tIA5F9+HD9+/A0VsX79ev/jTp064Z133imy/9ixYxgzZgw++eQTiKKIXbt2oXfv3jd0rpshCFyGgoiu5L71fuTcORxhe5bAU7U5XPHK/34iovIjaACLj49HfHw8EhMTUaVKFdkKWLp0KeLi4tC5c2c8+OCDSEpKgtFoxEMPPYR69erJdt5gBE7DJ6IgnHdPhDF1N2yb/glvdEP4om9TuyQiKqOCXoLUIiUuQc765jC2Hk/HF8Nayn4uunZs32tTRRxu8T27AAAgAElEQVQXnfMsIj/sAdEcjozEtZBMVrVLKqIijklZwHHRJjUvQQZdCT87W/6wo0lsgBHRVYiWqsjqthj6zGOwbXgWkES1SyKiMihoABsxYgQAYNq0aYoVoxVlpiVIRKrwVG8FZ+vnYf7ra4T98ora5RBRGRR0DlhISAgefvhhHD9+HIcOHQKQv5aXIAhYtWqVYgUqjQ0wIroWuY3/Af3Fg7DsWAhfZD244nuV/CIiogJBA9ibb76Jc+fOYerUqZg+fTrK0FSxm1aBvlQiulGCAEf7OdBnHINt4zj4ImrCW6WJ2lURURkR9BKkTqdD1apVsXjxYmzatAlvvfUWvv32W1SqVEnJ+hTHZSiI6JrpTcjq8QbEsMoI/3IIdI7TaldERGVE0ABWaOrUqThx4gTatGmDU6dO4fnnn1eiLtVwGQoiuh5SaDQy718KweNE+NrBENwV9ANMRHRdSgxgx48fx6RJk9ClSxdMnjwZJ06cUKIuVfESJBFdD190fWR1XwLDxUMI/3II4HOpXRIRaVyJAczlciE3NxcAkJeXB5/PJ3tRahLYACOiG+CJ64DsTq/AdOon2NY/A4jl+3clEd2coJPwCw0cONC/Kv2ff/6Jp59+Wom6VMUGGBHdCNdtD8ORcx7Wn2ZC2hINxz2z+FcdERWrxADWs2dPtGvXDidPnsQtt9yCyMhIJeoiIiqTcpsMhy7nHML2LIFktMF59wSGMCK6QokBDADsdjvsdrvctWhGRVpyg4hKn7P1vyC4sxG2axEknR45LZ5jCCOiIq4pgFUk/BVJRDdN0MHR4SVAEmHZsRAQdMhpMU7tqohIQ0oMYOnp6XA4HLDZbBWqC0ZEdFMEHRwd5wGSBMv2BYAk5YcwdsKICFcJYHv37sULL7wAURQRFhYGp9MJSZIwdepUNG3aVMkaFSUIApehIKLSIejg6DQfAGDZ8Rp0rgw47nkBEEr8ADoRlXNBA9icOXPw+uuvo1q1av5tp0+fxjPPPIOUlBRFilMD/zYlolJVEMKkEDvC9iyBkJeO7M4LAL1J7cqISEVBA5jX6y0SvgCgWrVqECpA+5wNMCIqVYIOzjZTIIbGwPrzi9DlpSOr+xuQTFa1KyMilQQNYO3bt8fjjz+ONm3awGazwel0YsuWLWjXrp2S9SmuAuRLIlJJbtOREEOjYds0HvY1DyHzvncgRtRUuywiUkHQADZq1Cj8/vvv2LlzJ9LT02G1WvHcc8+hUaNGStanCi5DQURycTVIgmithvBvRiAy5X5kdV8Czy1t1C6LiBR21U9BNmzYEA0bNlSqFiKiCsFT4x6k9/0CEV8OQcTnA+BoOw15dwxmC56oAgkawLZs2RL0RW3btpWlGK1g/4uI5CbaayOj72ewrX8ath+mwnTqZ2R3zJ+sT0TlX9AAtnr1auzbtw8tW7a8Yl95DmBchoKIlCKZbMi6722E7nkDlq0vIfLDe5HVdRG8sS3ULo2IZBY0gC1YsADJyckYOnQo6tSpo2RNquIFACJSlKBDbpMR8MTejfB1T8H+aV/kNBuNnOZPA3qz2tURkUyCrgao1+sxd+5cuN1uJevRBIkXIYlIYd4qdyL9ka/huu1hWHYsROSHPWA4u1PtsohIJlddjrlGjRqoX7++UrVoAufAEpFaJJMN2Z0XIPOB9yB4HLCv6QXLlukQ3Nlql0ZEpazEe0Hee++98Pl8l15gMKBatWoYP358+V2Sgg0wIlKRu2YnpPffAMvPcxD669sIOfwZHK0mwlU/kbcxIionSvxJvvvuuzFz5kx89dVXmD17Nu644w4MHz4cs2bNUqI+IqIKSTLZ4Gg/GxmJX8AXURPhG8fBnvIAjKd+Vrs0IioFJQawo0ePonXr1jCZTGjZsiXS0tLQqlUr6HTl968wNsCISCu8lRsjo88nyOr6OnQ552D/NBERn/Xj/DCiMq7ES5AmkwkrV65EkyZNsHv3bphMJuzbt6/IZcnyRACXoSAijREEuOJ7w1WnO0L3vY+wXf9B5JqH4IrriJy7noW3ajO1KySi61RiG+vll1/GsWPH8Morr+DkyZOYN28eLly4gBdffFGJ+hTHSfhEpFmGUOTeORQXkn+Co9VkGM/tQeSah2Bf8xBMf34BiF61KySia1RiBywyMhKtWrVCTEwMateujcjISLRv316J2lTDZSiISNOMYcht+iRybx+EkIOrEfbrW4j4ZgR8thrAXUMg1HwIUlgltaskoqsosQP2yiuv4OOPP4bRaMSnn36Kl156SYm6VMMGGBGVGSYL8hIG4+Kj3yOzx5vwWWOh3zgd0cvuQviXQ2A69i27YkQaVWIHbPv27Vi1ahUAYNCgQUhKSpK9KLVxDhgRlSk6Pdx1esBdpwfs3pPwbHsXIYc+gvnoNxBDY+Cq0wOuug/AE9sS0JX4a5+IFFDiT6LX64UoitDpdJAkCUI5nyRVzr88IirvYm6Ds80UOO+eANPxDTAf/gwhhz5C6P7lEEOj4arTA+5aXeGu3gowhqldLVGFVWIAu++++9C/f380btwYe/fuxX333adEXapiA4yIyjy9yd8Vy/bkwHR8I8xH1iLk0McI3f8+JL0ZntiWcMd1gLtGO/ii4rnIK5GCSgxgTzzxBNq2bYu//voLffv2RXx8vBJ1qYgtMCIqZ4xhcNd9AO66DyDbmwfjme0wHd8E04nvYP3xBQCAaLbDU+0ueKq1gCe2BbyVEgC9UeXCicqvoAHslVdeueJy4++//w4AGDt2rLxVqUziJDAiKq8MIfDUuAeeGvfAianQZZ+C8e8fYTyzDcbTv8B8bD0AQNKb4Y1uAG+lO+CtfEf+v1G3AXqTuvUTlRNBA1idOnWUrEMzOAeMiCoS0VYdrgZJcDXI/4CV4DwH45lfYDy7C4bzv8H8x6cI3b8cACDpjPBF3gpvZD34Cv7njaoLn70OoDer+WUQlTlBA1jv3r2VrEMzmL+IqCKTLJX9lyvzN4jQZR6HMW0fDGl7ob94GMZze2H+8wsIBTNmJUEHX3hNiBFx8Nni4AuPgy+8BsSImvDZakAKsav4FRFpEz+PXAxegSQiKiDoINprw2WvDVe9By9t9+ZCn/4XDOl/QJ/+Bwzpf0KXdRLm1F+hc2UUeQvRHAGfrQZEayxEazWIlqrwWavlP7ZWg89SlZ/IpAqHASwAL0ESEV0DQyh8lRrBV6nRFbsEVxZ0WSehzz4BfeYJ6LNPFDw/CeOZX64IaEB+SBMtVfMDWUFIyw9nl4KaZArnL2kqNxjAisEGGBHRjZPM4UHDGQDAkwu98wx0zrPQOc5A5ziT/9xxFjrnWZjO/w5dTpr/Eqf/fQ1hAZ2zKx9LIZEMaVQmMIAF4I8tEZHMjKHw2evkT94PxueGznkOOucZ6AtC2uWPjX9vgdmZCkESi7xM0pv9lzh9ETUhhteCL6IWfPZa8IXXhGQOl/mLI7o2DGDF4DIUREQq05sght8CMfwWBL2bpeiFLietIJydzQ9nzoKOmuM0TMc3QZ9zruhLQqLgi6gJn/1WeKNugze6PnzR9SFaqrJzRopSLIBJkoR27dqhVq1aAIA777wT48aNK3LMokWL8N1338FgMGDy5MlISEhQqrxLhMCmNxERaZLO4L8EGZTbCX3Wcegzj0GfWfjvMRj//gEhhz7yHyaaI+CNrg9vdEN4qzaFp2oziLYaDGUkG8UC2IkTJ9CoUSP897//LXb//v378csvvyAlJQVnzpzB6NGjsWbNGqXK8+OPGhFROWKywBfTEL6YhlfsEvLSYbhwEPoLB2G4cBCGiwcRemAVhN+WAgDE0ErwVG0KT9Xm+bdrimnA2zVRqVEsgO3fvx+pqalITk5GSEgIJk2aVGSx1507d6Jt27YQBAGxsbHw+Xy4ePEioqKilCrRj1cgiYjKPykkEp7qreCp3urSRtEL/YVDMKbuzF+M9uxOmI9+A/z8IsTQGLhrtIM7rj3ctbpyPhndFFkCWEpKCpYtW1Zk29SpUzFs2DD06NEDO3bswPjx44t0uBwOB+z2S4v1WSwWZGdnKx7A2AEjIqrAdAb/Jzjzbh+Yv8mZCuPJH2A6uRmmk5sRcvhjSDoT3HHt4ap7P9y1u0MyWVUunMoaWQJYYmIiEhMTi2zLzc2FXq8HADRv3hypqamQJMl/v0mr1Qqn0+k/3ul0wmazFXkPq9UMg0EvR8l+ISH5N5+127kooJbo9TqOiQZxXLSHYyIDe22gem3g7oHwSSLE07sgHPgMpgOfwvztekhGC6RGD0NsMghSbJNi34Ljok1qjotilyAXLVoEu92OoUOH4uDBg4iNjS1ys++mTZti/vz5GDJkCM6ePQtRFK/ofjkcLtnrdLk8AICMjBzZz0XXzm4P45hoEMdFezgmCrA0BJo3BJpNgOHsLoQcWImQfath2PMePFWaIqf503DX7FxkAj/HRZvkHpdKlWxB9ykWwIYNG4bx48dj8+bN0Ov1mDNnDgBg3rx56N69OxISEtC8eXM88sgjEEURU6dOVaq0Yl3enSMiIrqCoIO3WnM4qjWHs800mA99hLA9byBi7ePwxDRCzl1j4K7djZ+kpGIJUhla9CotLVv2c7z503G88fNx/DL2HgYwDeFfj9rEcdEejonKfB6Y//gUYTtfhyHjL7irt4aj7XTY6jbnuGiQmh0wfp42iDKTSomISDv0RrjqJyK9/0Zkt3sRhgsHELm6O3QbpgHePLWrIw1hAAtU0PQqO31BIiLSHJ0BeXcMwsVHf0Beg37Qb30dkat7wHDuV7UrI41gAAvAi45ERFRapBA7HB3nwdsvBYI7G/Y1vRCy/321yyINYAALgg0wIiIqLdKtnZHebz08t7SG7buJsG4aD/jk/2Q/aRcDWADOuyciIjlIIZHIvP89OJuNRujvKxHxxSAIbofaZZFKGMCC4SQwIiIqbTo9cu6egKzOC2A89TMiPusHIS9d7apIBQxgAQTOAiMiIpm56iciq/sbMFw4APsnfRnCKiAGsCDY/yIiIjm563RD5v3vQp9xFBFfDATczpJfROUGA1gAgctQEBGRQjw17kFWt8UwnPsVEV8P48T8CoQBjIiISEXuOt2R3XE+TCc3w/bdRHYAKgjF7gVZ1vA/fyIiUoqrwSNwZv8Ny/YF8MY0Qm7jf6hdEsmMHbAAnIJPRERqyLlrDFy1u8Hy40wYT25RuxySGQNYEGXoHuVERFQeCDpkd1kIX2RdhH8zAjrHabUrIhkxgAUQuBIrERGpRDJZkdXjTQg+N2zrnwZEn9olkUwYwIiIiDTEZ6+D7PYvwnR6K8J2LVK7HJIJA1iAwv4XL0ASEZFaXLf1RV69Xgj75VUYzmxXuxySAQNYAF6BJCIi1QkCHB3mQLRVR/j60RDc2WpXRKWMASwIzsEnIiI1SSYbsrq+Dp3jNCw/vqB2OVTKGMCIiIg0ylu1GXKbjEDo7ythOr5R7XKoFDGABSFxFhgREWmAs8U4eKNug3XTeAh5GWqXQ6WEASwAl6EgIiJN0ZuR3XkBdDnnYf1xhtrVUClhAAuCc8CIiEgrvJUTkNNsFEIOpsB0dJ3a5VApYAALwP4XERFpUU7zZ+CNbgjrdxMh5KWrXQ7dJAYwIiKiskBvQlaX16DLuwjr98+rXQ3dJAawAIVTwHgJkoiItMYX0xA5zZ9FyB+fwfTnF2qXQzeBAYyIiKgMyWn6FDyVG8O2eRJ0zlS1y6EbxAAWBJehICIiTdIbkd1lIQRPDqyb/slLNmUUA1gALkNBRERa54usC2eryTAf34CQAyvVLoduAANYgML4JfIPCiIi0rDchMFwV28Dy5YZ0GUeV7scuk4MYAF0hQmMAYyIiLRM0CG786uAoINtw1hA9KldEV0HBrAAhZcgRSYwIiLSONFWHY57XoDpzDaE7nlD7XLoOjCABSjsgPESJBERlQWu2/rCVacHLNvmwZC6R+1y6BoxgAXwd8CYwIiIqCwQBGR3nA8xrDLC1z0JwZWldkV0DRjAAuj9HTAGMCIiKhukEDuy7v0PdNmnYP1uApemKAMYwAIUdsD4ny4REZUl3mrN4Ww5HiF//g8hv3+gdjlUAgawADp2wIiIqIzKbfok3LfcA+sPU6G/cFDtcugqGMAC6Ao7YMxfRERU1gg6ZHVZCMlkQ/jXwyG4s9WuiIJgAAsg8FOQRERUhkmWysjqthj6zGOwffssIIlql0TFYAALoAM/BUlERGWbp3prONtMgfnoNwjbuUjtcqgYDGAB/B0wTsMnIqIyLDdhCPLieyNs23yYjm1QuxwKwAAWQF8wC58NMCIiKtMEAdkd5sEb0xC29aOhyziqdkV0GQawAP5lKDgLn4iIyjpjKLJ6vAkIOkR8OQSCK1PtiqgAA1iAwm8IO2BERFQeiOFxyOq+BPrMvxD+9XDA51a7JAID2BXYASMiovLGc0sbZHecD9PfW2D9bhLXWtIAg1InkiQJ7dq1Q61atQAAd955J8aNG1fkmBEjRiAjIwNGoxFmsxlvvfWWUuX58WbcRERUHrnqJ8KZeQyWHQshRtRETvOn1S6pQlMsgJ04cQKNGjXCf//736ses3btWn8XSg2FC7FyJXwiIipvclo8B33WCVi2zYPPGgtX/b5ql1RhKXYJcv/+/UhNTUVycjKGDh2Kv/76q8j+8+fPIysrCyNGjED//v2xadMmpUorQlfwHWEHjIiIyh1BQHanl+Gu3ga2jeNg+utrtSuqsGTpgKWkpGDZsmVFtk2dOhXDhg1Djx49sGPHDowfPx5r1qzx7/d4PHjiiScwcOBAZGZmon///khISEB0dLT/GKvVDINBL0fJfjZrDgDAYjHDbg+T9Vx07fR6HcdDgzgu2sMx0SZtjUsY0P8DSB/0Qfi6J+F75ENItdurXZQq1BwXWQJYYmIiEhMTi2zLzc2FXp8fnpo3b47U1FRIkuS/3BgTE4N+/frBYDAgOjoaDRo0wNGjR4sEMIfDJUe5ReQ488+RlZ2HjIwc2c9H18ZuD+N4aBDHRXs4JtqkvXHRQ+jxLuyf9IV+9aPIeGglvFWbqV2U4uQel0qVbEH3KXYJctGiRf6u2MGDBxEbG1tkrtdPP/2EZ599FgDgdDrxxx9/oE6dOkqV58c5YEREVBFIIZHI7PkBxLBKiPhfMgypu9UuqUJRbBL+sGHDMH78eGzevBl6vR5z5swBAMybNw/du3dH+/btsWXLFiQlJUGn02Hs2LGIiopSqjy/wkzI/EVEROWdaKmCjIc+hP2zJER8PgCZD75fITthahCkMrTgVVpatuzn2P13JoZ9+CsW9b0DLWtGyn4+ujbaa98TwHHRIo6JNml9XHTZpxHxWRJ0OeeR+eByeKvdpXZJiqgQlyDLCp2/A1ZmcikREdFNEW2xyOyVAtFSGfbPH4Xx1E9ql1TuMYAFuDQHTOVCiIiIFCRaqyGzVwp8tlsQ8fljMB1Zq3ZJ5RoDWAAd54AREVEFJVqqIKPPGngrJyD86xEI2fee2iWVWwxgAQR+CpKIiCowKSQSGT1Xwl2rM2ybJyNs28vsSsiAASwA7wVJREQVnjEUWT3eQm79R2DZ8RqsG58DfG61qypXFFuGoqxgB4yIiAiAzgBHp5chWqvBsuM16DOPIavHm5BClV8iqjxiBywAPwVJRERUQBCQ0/I5ZHVdBOO5PYj86AHoLxxSu6pygQEsAD8FSUREVJQrvhcyeqVA8OTCvuYhmI5tULukMo8BLABvRURERHQlb9WmSE9cC19ETYR/ORhhO14HJFHtssosBrAAvBURERFR8URbLDL6fAJX3Qdh2TYX4V8OgeDKVLusMokBLIC/AwYmMCIioisYw5DddREcbWfAdGITIlffB/3539WuqsxhAAvgX4aCXVUiIqLiCQJyGw9BRq8UwJuHyDU9YT60Ru2qyhQGsACCfx0wdsCIiIiuxlvtLqQnfQVP5TsR/u0zsH43CfDmqV1WmcAAFkBfkMCYv4iIiEomWSoj86FVyGkyAqH7lyPyowehT/9T7bI0jwEsABdiJSIiuk46A5ytn0fmA+9B50xF5OoeMB9YzW7GVTCABfDPAVO3DCIiojLHXbMT0vutg6dKE4RvHAvbt09DcDvULkuTGMACCP5LkEztRERE10u0VEVmz5VwthwP8x+fwb66Owxpv6ldluYwgAUo7ID52AIjIiK6MTo9cpo/g8xeKRB8Ltg/6onQX9/iwq2XYQALoAM7YERERKXBE9sS6Y+sg7tmJ1i3TEfE/x6DznlW7bI0gQEsgK7gO8KMTkREdPOkkEhk9XgL2R1egvHMdkSu7ALTkbVql6U6BrAAOs4BIyIiKl2CgLxGjyH9kW/gi6iJiK+Hw7ZhLAR3ttqVqYYBLMClhVjVrYOIiKi88dnrIKPPp3A2fwbmQx8h8sNuMJzZrnZZqmAAC8AOGBERkYz0RuS0HI+M3vm3LrJ/8jDCts0HfB6VC1MWA1iAggYYO2BEREQy8la7C+mPfAPXbX1h2bEQ9o97QZ/xl9plKYYBLICOK+ETEREpQjLZkN35VWR2+y/0mccQ+eG9CNm7tEIsV8EAFkCnYwAjIiJSkrvuA0jv9y08sXfD9sMURHzWD7qsv9UuS1YMYAF0nIRPRESkONFaDZkPLEd2h7kwnPsVkau6IOT3leX2fpIMYAEK54BxEj4REZHCBAF5jR5Fer9v4a18B2ybxiN87aByuXgrA1gAQRAgCOyAERERqUUMr4HMhz6Eo+0MmE79hMiVnWE+/Em56oYxgBVDJwjsgBEREalJ0CG38RCkP7IOPvutCF8/GuFfD4OQc17tykoFA1gxdALgY/4iIiJSXf7irZ/A0WoyTMc2IGpVZ5iOfKl2WTeNAawYAjtgRERE2qHTI7fpk0hP+hI+SzVEfD0MtvWjIeSlq13ZDWMAK4ZeEDgHjIiISGN80fWR0fd/cN41BuY//4eoDzqV2W4YA1gxdALXASMiItIkvRE5LcYhve9a+CxV8rth34wsc3PDGMCKkX8JUu0qiIiIKBhfpUb53bCWE2D+6xtErewI8+FPy8wnJRnAisEOGBERURmgNyKn+WikP/I1fBG1EL5+FMK/+keZWDeMAawYep0ALyeBERERlQm+qHhk9PkUjtZTYDrxXf66YQdWa7obxgBWjMgwEzJyPWqXQURERNdKp0duk+FI77ce3qj6CN84FhFfPAZd9im1KysWA1gxoq0mXHS61S6DiIiIrpPPXgeZvVOQfc9MGE9vR+TKzgjZ9z4giWqXVgQDWDFiLGZcyGEHjIiIqEwSdMhLGIyL/b+Ft8qdsG2eiIjP+kGXeVztyvwYwIoRYzXhYg47YERERGWZGB6HzJ4rkd1xHgxpvyFqVWeE7nkDEH1qlwaD2gVoUbTVDIfLh01/nIfJkJ9Rhcv2C0LxryP5WC1OOJwutcugABwX7eGYaBPHRWVh3RHS9k4k7HsRVX98AZ59H+Nix1dgt9+lWkmCVIbuuZOWlq3IeQ6l5yH5nV9QZr4xREREdA0kPKDbiunGZTge0hAJz32JjIwc2c5WqZIt6D7FApjP58OcOXOwb98+uN1ujB49Gh07dixyzOrVq7Fq1SoYDAaMHDnyiv1KBTC7PQwHT1zE+YKJ+Jd/hxjK1GGzmpHt4F+PWsNx0R6OiTZxXLRF73WiVoQJ1WvWUC2AKXYJ8rPPPoPX68WqVauQmpqKr776qsj+tLQ0LF++HGvWrIHL5cKAAQPQpk0bmEwmpUosomp4CKqGh6hybrqS3R4m6w8J3RiOi/ZwTLSJ46I14WoXoFwA27JlC+Lj4zFs2DBIkoQpU6YU2b937140adIEJpMJJpMJcXFxOHjwIBISEpQqkYiIiEgRsgSwlJQULFu2rMi2yMhImM1mLFmyBNu3b8ekSZOwYsUK/36HwwGb7VKrzmKxwOFwyFEeERERkapkCWCJiYlITEwssm3MmDHo0KEDBEFAixYtcOzYsSL7rVYrnE6n/7nT6SwSyPKPMcNg0MtRchF6vQ52e5js56FrxzHRJo6L9nBMtInjok1qjotilyCbNWuGzZs3o1u3bjh48CCqVatWZH9CQgJee+01uFwuuN1uHDlyBPHx8UWOcSg0gZHX6rWHY6JNHBft4ZhoE8dFm+QeF01Mwk9KSsK0adOQlJQESZIwY8YMAMDSpUsRFxeHzp07Izk5GQMGDIAkSRgzZgzMZrNS5REREREphuuAFYN/qWgPx0SbOC7awzHRJo6LNqnZAeOtiIiIiIgUxgBGREREpDAGMCIiIiKFMYARERERKYwBjIiIiEhhDGBERERECitTy1AQERERlQfsgBEREREpjAGMiIiISGEMYEREREQKU+xekFoniiKmT5+OQ4cOwWQyYdasWahZs6baZVUYHo8HkydPxqlTp+B2uzFy5EjUrVsXEydOhCAIqFevHqZNmwadTodFixbhu+++g8FgwOTJk5GQkKB2+eXehQsX0KdPH7zzzjswGAwcFw1YsmQJNm7cCI/Hg/79+6NFixYcFxV5PB5MnDgRp06dgk6nw8yZM/mzorJff/0VL7/8MpYvX47jx49f81gEO7bUSSRJkiR988030oQJEyRJkqTdu3dLI0aMULmiiuWjjz6SZs2aJUmSJF28eFFq3769NHz4cGnr1q2SJEnSlClTpHXr1kn79u2TkpOTJVEUpVOnTkl9+vRRs+wKwe12S08++aR07733Sn/++SfHRQO2bt0qDR8+XPL5fJLD4ZD+/e9/c1xUtn79eunpp5+WJEmStmzZIo0aNYpjoqI33nhDeuCBB6TExERJkqTrGovijpUDL0EW2LlzJ+655x4AwJ133ol9+/apXFHF0r17dzzzzDP+53q9Hvv370eLFi0AAO3atcNPP/2EnTt3om3bthAEAbGxsfD5fLh48aJaZVcIc+fORb9+/VC5cmUA4LhowJYtWxAfH4+nnhr1CV8AAAW1SURBVHoKI0aMQIcOHTguKqtduzZ8Ph9EUYTD4YDBYOCYqCguLg6vv/66//n1jEVxx8qBAayAw+GA1Wr1P9fr9fB6vSpWVLFYLBZYrVY4HA48/fTTePbZZyFJEgRB8O/Pzs6+YpwKt5M8Pv74Y0RFRfn/OAHAcdGA9PR07Nu3DwsXLsSMGTPw3HPPcVxUFhYWhlOnTqFHjx6YMmUKkpOTOSYq6tatGwyGS7OsrmcsijtWDpwDVsBqtcLpdPqfi6JYZPBIfmfOnMFTTz2FAQMG4MEHH8T8+fP9+5xOJ8LDw68YJ6fTCZvNpka5FcKaNWsgCAJ+/vlnHDhwABMmTCjy1zrHRR12ux116tSByWRCnTp1YDabcfbsWf9+jovy3n33XbRt2xbjxo3DmTNnMGjQIHg8Hv9+jom6Lp/DVdJYFHesLDXJ8q5lUNOmTfH9998DAPbs2YP4+HiVK6pYzp8/jyeeeALjx49H3759AQANGzbEtm3bAADff/89mjdvjqZNm2LLli0QRRGnT5+GKIqIiopSs/RybcWKFXj//fexfPlyNGjQAHPnzkW7du04Lipr1qwZfvjhB0iShNTUVOTm5qJVq1YcFxWFh4f7g1RERAS8Xi9/h2nI9YxFccfKgSvhFyj8FOThw4chSRJmz56NW2+9Ve2yKoxZs2bhq6++Qp06dfzb/vWvf2HWrFnweDyoU6cOZs2aBb1ej9dffx3ff/89RFHEpEmTZPvhoKKSk5Mxffp06HQ6TJkyheOisnnz5mHbtm2QJAljxozBLbfcwnFRkdPpxOTJk5GWlgaPx4OBAwfi9ttv55io6O+//8bYsWOxevVqHD169JrHItixpY0BjIiIiEhhvARJREREpDAGMCIiIiKFMYAREf1/e/cPklwbxnH8axYtBklCCC3lUigR/aGhqTJRoiGwKEFoiOYciiyDlhYpGoSWQji15GJDkBA01JhuIVIQQUUQDVFRoEi90xNPQ+/w8OT7Gr/PfHNxXWf6cc597ltEpMQUwERERERKTAFMREREpMQUwESkLB0fH5NIJEgkEp8OvPwTZ2dnpNNpAEKhEIVC4W+0KCLyJR1DISJlra+vj1QqRXV19R/XiMVi2Gw2xsfH/2JnIiJf0107IlKWkskk0WiUl5cXQqEQ6+vrrK6ukk6neX9/Z2JiAp/PRzAYxGq18vT0RCwWIxKJ8Pz8zMPDAyMjI/T397O7u0tVVRVOp5Pp6WlSqRT39/csLCxQLBYxmUxEIhGam5vxeDy0t7dzeXlJXV0dsViMq6srwuEwlZWVmM1motEo9fX1//UjEpH/MQUwESlbfr+f/f191tbWODo64ubmhp2dHfL5PKOjo/T09AAwNDTEwMAA2WyWwcFBPB4Pd3d3BINBAoEAw8PD2Gw2WltbP2pHo1GCwSBut5tcLsf8/DzJZJLr62sMw8ButzM2Nsbp6SnZbBan08nc3ByZTIbHx0cFMBH5VwpgIvIjnJ+fk81mCQaDABSLRW5vbwFobGwEwGazYRgGBwcHWCwWisXil/UuLi7o6uoCoKWl5eOya6vVit1uB8But5PP5/H7/WxsbDA5OUlNTQ2hUOjb5hSRn0Gb8EWkrJlMJt7e3mhqaqK7u5vt7W0Mw8Dn89HQ0PCxBiAej9PW1sbKygper5dfW2B/1fidw+Egk8kAkMvlsNlsn2r97vDwkI6ODgzDwOv1srm5+W3zisjPoDdgIlLWOjs7mZqaYmtri5OTEwKBAK+vr7jdbiwWy6e1vb29LC0tsbe3R21tLWazmUKhgMvlIhqN4nA4PtbOzs6yuLhIPB6nWCyyvLz8ZQ8ul4uZmRlisRgVFRWEw+Fvm1dEfgb9BSkiIiJSYvoEKSIiIlJiCmAiIiIiJaYAJiIiIlJiCmAiIiIiJaYAJiIiIlJiCmAiIiIiJaYAJiIiIlJiCmAiIiIiJfYPAD2oKI/d42kAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x25b082bd390>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "sns.set_style(style='darkgrid')\n",
    "plt.figure(figsize=(10,6))\n",
    "plt.plot(np.log10(epsilon2))\n",
    "plt.plot(np.log10(epsilon1))\n",
    "plt.ylabel('log10 of the MSE')\n",
    "plt.xlabel('iterations')\n",
    "plt.legend(['1 known factor','all unknown'])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## With Tensorly Unfold"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "itermax=1000\n",
    "epsilon1=[]\n",
    "A2=np.random.rand(I2,Rank)\n",
    "A3=np.random.rand(I3,Rank)\n",
    "for ii in range(itermax):\n",
    "    krtransposed=linalg.pinv(tl.transpose(tl.tenalg.khatri_rao([A3,A2])))\n",
    "    A1=np.dot(tl.unfold(tensor=Ts,mode=0),krtransposed)\n",
    "    krtransposed=linalg.pinv(tl.transpose(tl.tenalg.khatri_rao([A3,A1])))\n",
    "    A2=np.dot(tl.unfold(tensor=Ts,mode=1),krtransposed)\n",
    "    krtransposed=linalg.pinv(tl.transpose(tl.tenalg.khatri_rao([A2,A1])))\n",
    "    A3=np.dot(tl.unfold(tensor=Ts,mode=2),krtransposed)\n",
    "    error=tlmetrics.MSE(y_pred=tl.tenalg.multi_mode_dot(eye,matrix_or_vec_list=[A1, A2, A3]),y_true=Ts)\n",
    "    epsilon1.append(error)\n",
    "Xhat3=tl.tenalg.multi_mode_dot(eye,matrix_or_vec_list=[A1, A2, A3])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[[ 0.29249063,  0.31474527,  0.31713509],\n",
       "        [ 0.34854682,  0.37506659,  0.37791441],\n",
       "        [ 0.24483378,  0.26346237,  0.2654628 ]],\n",
       "\n",
       "       [[ 0.13816617,  0.14867878,  0.14980767],\n",
       "        [ 0.16464589,  0.17717325,  0.1785185 ],\n",
       "        [ 0.11565412,  0.12445386,  0.12539881]],\n",
       "\n",
       "       [[ 0.13472973,  0.14498087,  0.14608169],\n",
       "        [ 0.16055085,  0.17276663,  0.17407842],\n",
       "        [ 0.11277759,  0.12135846,  0.12227992]]])"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Xhat3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[[ 0.32209341,  0.37541617,  0.42216749],\n",
       "        [ 0.38393694,  0.45382583,  0.44820815],\n",
       "        [ 0.27053073,  0.3317248 ,  0.27873834]],\n",
       "\n",
       "       [[ 0.19532735,  0.16447497,  0.22632321],\n",
       "        [ 0.2106079 ,  0.18722507,  0.22460432],\n",
       "        [ 0.12513622,  0.12310625,  0.11807863]],\n",
       "\n",
       "       [[ 0.13211993,  0.15192751,  0.12762588],\n",
       "        [ 0.20579037,  0.23259446,  0.13993744],\n",
       "        [ 0.18177979,  0.20517857,  0.0910915 ]]])"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Ts"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}