lectures/lecture17/lecture-17.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Lecture 17:  Tensors and tensor decompositions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Recap of the previous lecture\n",
    "- Large scale eigenvalue problem\n",
    "- Krylov methods\n",
    "- Jacobi-Davidson method\n",
    "- LOBPCG\n",
    "- Intro to streaming algorithms"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Today lecture \n",
    "\n",
    "- Tensor decompositions and applications\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Tensors\n",
    "\n",
    "By **tensor** we imply a **multidimensional array**:\n",
    "\n",
    "$$\n",
    "A(i_1, \\dots, i_d), \\quad 1\\leq i_k\\leq n_k,\n",
    "$$\n",
    "\n",
    "where $d$ is called dimension, $n_k$ - mode size.\n",
    "This is standard definition in applied mathematics community. For details see [[1]](http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.153.2059&rep=rep1&type=pdf), [[2]](http://arxiv.org/pdf/1302.7121.pdf), [[3]](http://epubs.siam.org/doi/abs/10.1137/090752286).\n",
    "\n",
    "* $d=2$ (matrices) $\\Rightarrow$ classic theory (SVD, LU, QR, $\\dots$)\n",
    "\n",
    "* $d\\geq 3$ (tensors) $\\Rightarrow$ under development. Generalization of standard matrix results is not  **straightforward**.\n",
    "<img src=\"./tensor_dogs.png\">\n",
    "\n",
    "Picture is taken from [this presentation](http://users.cecs.anu.edu.au/~koniusz/tensors-cvpr17/present/anandkumar_anima_tmcv2017.pdf)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## More formal definition\n",
    "\n",
    "Tensor is a multilinear form. When you fix the basis, you get a $d$-dimensional table."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Curse of dimensionality\n",
    "\n",
    "The problem with multidimensional data is that number of parameters grows <font color='red'> exponentially </font> with $d$:\n",
    "\n",
    "\n",
    "$$\n",
    "    \\text{storage} = n^d.\n",
    "$$\n",
    "For instance, for $n=2$ and $d=500$\n",
    "$$\n",
    "    n^d = 2^{500} \\gg 10^{83}  -  \\text{ number of atoms in the Universe}\n",
    "$$\n",
    "\n",
    "Why do we care? It seems that we are living in the 3D World :)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Applications\n",
    "\n",
    "#### Quantum chemistry\n",
    "\n",
    "Stationary Schroedinger equation for system with $N_{el}$ electrons\n",
    "\n",
    "$$\n",
    "    \\hat H \\Psi = E \\Psi,\n",
    "$$\n",
    "\n",
    "where\n",
    "\n",
    "$$\n",
    "\\Psi = \\Psi(\\{{\\bf r_1},\\sigma_1\\},\\dots, \\{{\\bf r_{N_{el}}},\\sigma_{N_{el}}\\})\n",
    "$$\n",
    "\n",
    "3$N_{el}$ spatial variables and $N_{el}$ spin variables. \n",
    "<img src=\"./large_mol.jpg\" width=600>\n",
    "\n",
    "* Drug and material design\n",
    "* Predicting physical experiments"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Uncertainty quantification\n",
    "\n",
    "Example: oil reservoir modeling. Model may depend on parameters $p_1,\\dots,p_d$ (like measured experimentally porocity or temperature) with uncertainty\n",
    "\n",
    "$$\n",
    "u = u(t,{\\bf r},\\,{p_1,\\dots,p_d})\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### And many more\n",
    "\n",
    "* Signal processing\n",
    "* Recommender systems\n",
    "* Neural networks\n",
    "* Language models\n",
    "* Financial mathematics\n",
    "* ..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Working with many dimensions\n",
    "\n",
    "How to work with high-dimensional functions?\n",
    "\n",
    "- **Monte-Carlo**: class of methods based on random sampling. Convergence issues\n",
    "- **Sparse grids**: special types of grids with small number of grid points. Strong regularity conditions\n",
    "- **Best N-term approximation** : sparse expansions in certain basis.\n",
    "- Promising approach based on <font color='red'>tensor decompositions </font>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Tensor decompositions\n",
    "\n",
    "## 2D\n",
    "\n",
    "Skeleton decomposition:\n",
    "$$\n",
    "A = UV^T\n",
    "$$\n",
    "or elementwise:\n",
    "$$\n",
    "a_{ij} = \\sum_{\\alpha=1}^r u_{i\\alpha} v_{j\\alpha}\n",
    "$$\n",
    "leads us to the idea of **separation of variables.**\n",
    "\n",
    "**Properties:**\n",
    "* Not unique: $A = U V^T = UBB^{-1}V^T = \\tilde U \\tilde V^T$\n",
    "* Can be calculated in a stable way by **SVD**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Canonical decomposition\n",
    "\n",
    "The most straightforward way to generize separation of variables to many dimensions is the **canonical decomposition**: (CP/CANDECOMP/PARAFAC)\n",
    "\n",
    "$$\n",
    "a_{ijk} = \\sum_{\\alpha=1}^r u_{i\\alpha} v_{j\\alpha} w_{k\\alpha}\n",
    "$$\n",
    "\n",
    "minimal possible $r$ is called the **canonical rank**. Matrices $U$, $V$ and $W$ are called **canonical factors**. This decomposition was proposed in 1927 by Hitchcock, [link](https://onlinelibrary.wiley.com/doi/abs/10.1002/sapm192761164)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "**Properties**:\n",
    "\n",
    "* For a $d$ dimensional tensor memory is $nrd$\n",
    "* Unique under mild conditions\n",
    "* Set of tensors with rank$\\leq r$ is not closed (by contrast to matrices): <br>\n",
    "  $a_{ijk} = i+j+k$, $\\text{rank}(A) = 3$, but\n",
    "  \n",
    "  $$a^\\epsilon_{ijk} = \\frac{(1+\\epsilon i)(1+\\epsilon j)(1+\\epsilon k) - 1}{\\epsilon}\\to i+j+k=a_{ijk},\\quad \\epsilon\\to 0 $$\n",
    "  \n",
    "  and $\\text{rank}(A^{\\epsilon}) = 2$\n",
    "* No stable algorithms to compute best rank-$r$ approximation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Alternating Least Squares algorithm\n",
    "\n",
    "0. Intialize random $U,V,W$\n",
    "1. fix $V,W$, solve least squares for $U$\n",
    "2. fix $U,W$, solve least squares for $V$\n",
    "3. fix $U,V$, solve least squares for $W$\n",
    "4. go to 2."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### DNN compression ([Lebedev, et. al 2015](https://arxiv.org/pdf/1412.6553.pdf))\n",
    "\n",
    "- Convolution is represented as 4D tensor\n",
    "- Apply CP decomposition to this tensor\n",
    "- Apply one-by-one convolutions with smaller kernels given by CP decomposition\n",
    "- Fine tune such model\n",
    "\n",
    "<img src=\"./cp_compress_dnn.png\" width=600>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[[ 0  1]\n",
      "  [ 2  3]\n",
      "  [ 4  5]\n",
      "  [ 6  7]]\n",
      "\n",
      " [[ 8  9]\n",
      "  [10 11]\n",
      "  [12 13]\n",
      "  [14 15]]\n",
      "\n",
      " [[16 17]\n",
      "  [18 19]\n",
      "  [20 21]\n",
      "  [22 23]]]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Approximation error')"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEGCAYAAAB7DNKzAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAzu0lEQVR4nO3dd3yV5f3/8dcnizBC2CuDPWUTEhVEVBREhlsBrQNFW2dtax39trXL1rZuqlVEHIAiVplCHQiCSEjYKxBWEgIk7Jn9+f2RQ39pSsJJck7uMz7Px+N+JOfmnPt+HyDnk+u+rvu6RFUxxhhjKhLidABjjDG+zQqFMcaYSlmhMMYYUykrFMYYYyplhcIYY0ylwpwO4A3NmjXTdu3aOR3DGGP8Smpq6iFVbV5+f0AWinbt2pGSkuJ0DGOM8Ssisvd8++3SkzHGmEr5fKEQkQ4i8o6IzHY6izHGBCNHCoWITBWRHBHZVG7/CBFJE5F0EXkKQFV3qepEJ3IaY4xxrkUxDRhRdoeIhAKTgWuBHsA4EelR+9GMMcaU5UihUNVlwJFyuxOBdFcLogD4CBjr7jFFZJKIpIhISm5urgfTGmNMcPOlPooYILPM4ywgRkSaisibQD8RebqiF6vqW6qaoKoJzZv/z+guY4wx1eRLw2PlPPtUVQ8DD7p1AJHRwOhOnTp5NJgxxgQzX2pRZAFxZR7HAtkOZTHGGL9yKr+I5+ZtJudEnseP7UuFYjXQWUTai0gEcDswtyoHUNV5qjopOjraKwGNMcZXvbV0J++u2EP28QApFCIyE1gJdBWRLBGZqKpFwMPAYmArMEtVN1fxuKNF5K3jx497PrQxxviogyfyePu73VzXuzV94xp5/PiO9FGo6rgK9i8EFtbguPOAeQkJCfdX9xjGGONvXvpyO0UlJfxyeDevHN+XLj057kxBEWcKipyOYYwxbtt+8CSzUjK58+J2xDet55VzBFShqOmlp9/P38KoV5ezLvOYZ4MZY4yXPL9wK/XrhPHIld4b7RlQhaKmndlj+sSQV1jMTW98zytf7aCouMTDCY0xxnO+Tz/EkrRcHrqiE43rR3jtPAFVKGraorikY1O+eHwIY/q04aWvtnPzmyvZfei0h1MaY0zNlZQof/piKzGN6nL3pe28eq6AKhSeGB4bXTecl27ry+vj+7H70GlGvvId01ftRVU9mNQYY2pm7vpsNu07wc+u6UJkeKhXzxVQhcKTRvVuw+LHh5DQrjHPfraJ+95LIfdkvtOxjDGGvMJi/ro4jYvaNOT6vjFeP19AFQpP30fRKjqS9+5J5Deje7A8/RDDX17Gvzcf8MixjTGmut5fuYd9x87yzMjuhIScb/YjzwqoQuGNO7NDQoR7BrVn/iODaR0dyaQPUvnl7A2cyrdhtMaY2nfsTAGvf5PO5V2aM6hTs1o5Z0AVCm/q3DKKz34yiJ8M7cgnqZmMfOU7UveWnyndGGO86/Vv0jmVX8TTI71zc935WKGogoiwEJ4c0Y2PH7gERbnlzZX8dfE2CopsGK0xxvsyj5zh/ZV7uXlALN1aNay18wZUoaituZ4GtmvCwkcv4+YBsUxespMb31hBes5Jr57TGGNeWJxGSAg8cXXXWj1vQBWK2pw9NioynBdu7sObdwwg+1ge1726nGkrdlNSYsNojTGetz7zGPPWZ3Pf4A60io6s1XMHVKFwwoierVj0+GVc2rEpv523hbveTeaAF6b5NcYEL1XlTwu30rR+BA9c3qHWz2+FwgNaREUy9e6B/OH6nqTsOcrwl5exYMN+p2MZYwLE11tzWLX7CI8P60xUZHitn98KhYeICHdc3JYFjw6mXbP6PDRjDU98vI4TeYVORzPG+LGi4hL+vGgbHZrV5/bEeEcyWKHwsA7NGzD7wUt4fFhn5qzP5tqXv+OHXYedjmWM8VOzUrJIzznFkyO6ER7qzEd2QBUKX1nhLjw0hMeHdWH2g5cQERbCuLd/4PmFW8kvKnY0lzHGv5zOL+LFL7eT0LYxwy9q6ViOgCoUvrZmdr/4xix4dDDjE+P557JdjH19BdsOnHA6ljHGT7y1bBeHTuXz9MjuiHh/qo6KBFSh8EX1IsL44w29mHp3AodOFTDmtRVM+W6XDaM1xlQq50Qeb3+3i5G9WjGgbWNHs1ihqCVXdmvJ4scvY2jX5vxhwVYmTFnFvmNnnY5ljPFRL321g4KiEp700jrYVWGFohY1bVCHf945gBdu7s2GrGOMeHkZn6/dZ2tdGGP+y46DJ/l4dQZ3XNyWds3qOx3HCkVtExFuTYjji8eG0LVlFI9/vI5HZq7l2JkCp6MZY3zEXxZto35EGI9e1dnpKIAVCsfEN63Hxw9cwi+Gd2XRpgOMePk7lu845HQsY4zDfth1mK+25vDjKzrSxIvrYFeFzxcKEakvIu+JyNsiMsHpPJ4UGiI8dEUnPn9oEA0iw7jjnVU8N28zeYU2jNaYYFRSUjpVR+voSO4d1N7pOP/hSKEQkakikiMim8rtHyEiaSKSLiJPuXbfCMxW1fuBMbUethb0jIlm/iODufvSdry7Yg+jX1vOpn3O3gtijKl98zfuZ0PWcX52TVevr4NdFU61KKYBI8ruEJFQYDJwLdADGCciPYBYINP1tID9VTsyPJTfjrmI9+9N5EReITf8YwWTl6RTbMNojQkK+UXFvLBoG91bN+SGft5fB7sqHCkUqroMKL88XCKQrqq7VLUA+AgYC2RRWiygkrwiMklEUkQkJTc31xuxa8WQLs1Z/PgQrunRir8uTuO2f64k88gZp2MZY7zsg5V7yTp6lmdGdiO0FtbBrgpf6qOI4f+3HKC0QMQA/wJuEpE3gHkVvVhV3wKeA9ZERPhGB1B1NaoXwevj+/HSbX1IO3CSES8vY1ZKpg2jNSZAHT9TyGvfpHNZ52Zc1rm503H+hy8VivOVUFXV06p6j6r+WFWnV3YAX5vCoyZEhBv6xbLop0PoFRvNk7M38OCHqRw5bcNojQk0k79N50ReIc+M7O50lPPypUKRBcSVeRwLZFflAL4yKaAnxTSqy4z7LuaZkd1Ysi2Xa15axpJtOU7HMsZ4SOaRM0xbsYeb+sfSvXXtrYNdFb5UKFYDnUWkvYhEALcDc6tygEBqUZQVEiJMGtKROQ8PolmDCO6Ztppffb6RMwVFTkczxtTQ3/+dhgj87JouTkepkFPDY2cCK4GuIpIlIhNVtQh4GFgMbAVmqermKh434FoUZXVv3ZDPHxrEpCEdmL4qg1GvLmdd5jGnYxljqmlj1nE+X5fNxMHtaR1d1+k4FZJA7CBNSEjQlJQUp2N41cqdh/nZrHUcPJnPo1d25qErOhLm0KImxpiqU1XGv72KtIMn+fYXQ2nowBKn5YlIqqomlN8fUJ8sgd6iKOuSjk354vEhjOnThpe+2s7Nb65k96HTTscyxrjp27RcVu46zKNXdvKJIlGZgCoUgdpHUZHouuG8dFtfXh/fj92HTjPyle+YvmqvDaM1xscVFZfw/Bdbade0HuOT2jod54ICqlAEq1G927D48SEktGvMs59t4r73Usg9me90LGNMBWanZrH94Cl+OaIbEWG+/zHs+wmrIJguPZXXKjqS9+5J5Leje7A8/RDDX17GvzcfcDqWMaacMwWl62D3j2/EiJ6tnI7jloAqFMF26am8kBDh7kHtmf/IYFpHRzLpg1R+OXsDp/JtGK0xvmLKd7vJOZnPs9c5uw52VQRUoTClOreM4rOfDOKhKzrySWomI1/5jtS95afWMsbUttyT+fxz6U5GXNSKAW2bOB3HbQFVKIL50lN5EWEh/GJ4N2Y9cAmKcsubK/nr4m0UFJU4Hc2YoPXK19vJLyrhyRFdnY5SJQFVKIL90tP5JLRrwhePDeHmAbFMXrKTG99YQXrOSadjGRN00nNOMTM5k/FJ8XRo3sDpOFUSUIXCnF+DOmG8cHMf3rxjANnH8rju1eVMW7GbElvrwpha85dF26gbHuoz62BXhRWKIDKiZysWPX4Zl3Zsym/nbeGud5M5cDzP6VjGBLzk3Uf4cstBHry8A80a1HE6TpUFVKGwPooLaxEVydS7B/KH63uSsucow19expI0m43WGG9RLV0Hu1XDSCYO7uB0nGoJqEJhfRTuERHuuLgtCx4dTKuGkfxy9gYKi62T2xhvWLBxP+syj/HENV2oG+E762BXRUAVClM1HZo34MkRXck5mc/XWw86HceYgFNQVMILi9Lo1iqKm/rHXvgFPsoKRZAb2rUFbaIjmb4qw+koxgScD3/YS8aRMzx1re+tg10VViiCXGiIcNvAeL7bcYi9h232WWM85fjZQl79ZgeDOzXj8i6+tw52VVRaKEQkVES+qq0wxhm3DYwjNESYmZzpdBRjAsYb3+7k+NlCnrq2m99M1VGRSguFqhYDZ0TEL3qHbdRT9bSKjuSqbi34JCXT7tw2xgP2HTvL1BW7uaFvDD1j/OLjs1LuXHrKAzaKyDsi8uq5zdvBqsNGPVXfhIvbcvh0AYttxlljauzvi9MA+Nlw/5qqoyJhbjxngWszAeyyTs2Ia1KXGasyGN2njdNxjPFbm/Yd57N1+3hgSEdiGvnuOthVccFCoarviUgE0MW1K01VC70by9S2kBBhXGI8LyxKY2fuKTr62Vw0xvgCVeX5L7bSqG44Px7a0ek4HnPBS08iMhTYAUwG/gFsF5Eh3o1lnHDLgDjCQoSZNlTWmGpZuj2XFemHeeTKzkTX9e11sKvCnT6KvwPXqOrlqjoEGA685N1YxgnNo+ow/KJWzF6TRV5hsdNxjPErxSXK8wu3Ed+kHndc7PvrYFeFO4UiXFXTzj1Q1e1ArZVKEeng6kifXVvnDGYTkuI5dqaQLzbtdzqKMX7l0zVZpB08yZMjuvrFOthV4c67SXV9UA91bW8Dqe4cXESmikiOiGwqt3+EiKSJSLqIPFXZMVR1l6pOdOd8puYu6diU9s3qM8MuPxnjtrMFxfz932n0jWvEdb1aOx3H49wpFA8Cm4FHgceALa597pgGjCi7Q0RCKe3vuBboAYwTkR4i0ktE5pfbWrh5HuMhIsL4xHhW7znK9oO2wJEx7nhn+S4OnsjnmZH+sw52VVzozuwQIFVVX1TVG1X1BlV9SVXz3Tm4qi4Dyi/WnAiku1oKBcBHwFhV3aiqo8ptNv+1A24aEEtEaIi1Koxxw6FT+by5dBdX92hJYnv/WQe7Ki50Z3YJsF5E4j14zhig7FwRWa595yUiTUXkTaCfiDxdyfMmiUiKiKTk5uZ6Lm0QalI/gmt7teLTNVmcLbBObWMq8+rXOzhbWMxT13ZzOorXuHPDXWtgs4gkA/+ZNU5Vx1TznOdrl1W4JqeqHsaNS12q+paI7AdGR0REDKhmNuMyIaktc9ZlM29DNrcmxDkdxxiftCv3FDNWZTAuMS6g7z1yp1A85+FzZgFlP3ligWxPHFhV5wHzEhIS7vfE8YLZwHaN6dSiATNWZVihMKYCLyxKo05YCI9d1eXCT/Zj7vRRTFbVpeW3GpxzNdBZRNq77vi+HZhbg+OVzWuTAnqIiDAhKZ51mcfYnG1/n8aUl7LnCIs2H+CByzvSPMr/1sGuCq/2UYjITGAl0FVEskRkoqoWAQ8Di4GtwCxV3Vyd458nr00K6EE39oulTph1ahtT3rl1sFtE1eG+y9o7HcfrvNpHoarjKti/EFjobkh3ichoYHSnTp08feigFF0vnNF92vD52n08PbI7Deq489/FmMC3aNMB1mQc4y839aJeROD/XDjRR+E11kfheeOT4pmdmsXcddmMT/Lk4Ddj/FNBUQl/WbSNLi0bcPOA4Oi/u+ANd67+iD2UTuWxlNI+hjVezlUt1kfhef3iGtGtVRTTV+1FtcLBacYEjRmr9rLn8Bmevra7X6+DXRXuzB57PzAb+KdrVwzwuRczVZv1UXieiDDh4rZszj7BhiwrwCa4ncgr5NVv0rm0Y1OGdvXvdbCrwp0pPB4CBgEnAFR1B+CTU2tYi8I7ru/bhnoRodapbYLem9/u5MjpAp6+NjCn6qiIO4Ui3zXVBgAiEkYlN8g5yVoU3hEVGc7Yvm2Yuz6bE3m2ZpUJTtnHzvLO8t1c37cNvWKD6zPGnUKxVESeAeqKyNXAJ8A878YyvmZ8YlvOFhbz+dp9TkcxxhEvfrkdVfh5gKyDXRXuFIqngFxgI/AApcNaf+XNUMb39IqNpldMNNN/yLBObRN0tmSf4NM1Wdw9qB2xjes5HafWuTPqqURV31bVW1T1Ztf3PvlJYX0U3jUhKZ60gydZk3HU6SjG1Krnv9hKw8hwHhoanPdoBdQyTNZH4V2j+7ShQZ0wpv9gndomeCzbnst3Ow7xyJWdiK4XOOtgV0VAFQrjXfXrhHFDvxjmb9zPsTMFF36BMX6uuER5/ottxDWpy52XBNY62FVhhcJUyfikeAqKSvh0jXVqm8D32dp9bN1/gl8M70adsFCn4zjGnRvuuojI2yLybxH55txWG+GqyvoovK9764b0j29kd2qbgJdXWLoOdp/YaEYF4DrYVeFOi+ITSqfs+BXwizKbz7E+itoxPqktu3JPs2p3+VVujQkcU1fsZv/xPJ4e2Z2QIJmqoyLuFIoiVX1DVZNVNfXc5vVkxmeN6t2ahpFhTLc7tU2AOnwqnzeW7GRY9xZc3KGp03Ec506hmCciPxGR1iLS5Nzm9WTGZ0WGh3LTgFgWbdrP4VP5TscxxuNe+yadMwG+DnZVuFMo7qL0UtP3QKprS/FmKOP7JiTFU1iszE7NcjqKMR6159BpPvxhL7cNjKNTiyin4/gEd264a3+erUNthKsq68yuPZ1aRJHYvgkzkjMoKbFObRM4Xli8jYiwEB4f1tnpKD7DnVFP4SLyqIjMdm0Pi4hP3nVindm1a0JSPHsPn+H7nYedjmKMR6zJOMrCjQeYNKQDLaIinY7jM9y59PQGMAD4h2sb4NpngtyInq1oXC+c6av2Oh3FmBpTVf60YCvNo+pw/2U+edHEMe4shTpQVfuUefyNiKz3ViDjP+qEhXJLQhxTl+8m50QeLRrab2DGfy3efJCUvUf50w29qG/rw/8Xd1oUxSLS8dwDEekAFHsvkvEn4xLjKSpRZqVkOh3FmGorLC5dB7tTiwbcmhDrdByf406h+AWwRES+FZGlwDfAz7wby/iL9s3qM6hTU2YmZ1JsndrGT32UnMHuQ6d5+tpuhIXazEbluTPq6WugM/Coa+uqqku8Hcz4j/GJbdl37CzLduQ6HcWYKjuZV8jLX+0gqX0Truzmk6s8O67CQiEiV7q+3ghcB3QCOgLXufbVGhG53jXf1BwRuaY2z20u7OoeLWnWoI5NP2780j+X7uLw6QKevS641sGuispaFJe7vo4+zzbK3ROIyFQRyRGRTeX2jxCRNBFJF5GnKjuGqn6uqvcDdwO3uXtuUzsiwkK4NSGWb7YdZP/xs07HMcZtB47nMWX5Lsb0aUPv2EZOx/FZFXbtq+pvXN/+TlV3l/0zEWlfhXNMA14H3i/z+lBgMnA1kAWsFpG5QCjwfLnX36uqOa7vf+V6nfEx4xLjeWPpTj5KzuSnV3dxOo4xbnnxyzRKSuAXQbgOdlW402vz6Xn2zXb3BKq6DCg/zWgikK6qu1S1APgIGKuqG1V1VLktR0r9BfhCVdec7zwiMklEUkQkJTfXrpXXtrgm9RjSuTkfr86kqLjE6TjGXNC2Ayf4JDWLH13SlrgmwbcOdlVU1kfRTURuAqJF5MYy291ATQfMxwBlx1NmufZV5BFgGHCziDx4vieo6luqmqCqCc2bN69hPFMd45PiOXAij2+25Vz4ycY47M9fbCOqThgPXxmc62BXRWV3lXSltC+iEaX9EuecBO6v4XnP12NU4dhKVX0VePWCBxUZDYzu1Mn+4Z1wVbcWtGxYhxnJGVxzUSun4xhToRXph/g2LZdnRnajUb0Ip+P4vMr6KOYAc0TkElVd6eHzZgFxZR7HAtkePoepZWGhIdw2MJ7XvtlB5pEz1pw3PqmkRPnTwq3ENKrLjy5p53Qcv+BOH8VaEXlIRP7hGsE0VUSm1vC8q4HOItJeRCKA24G5NTymTQroA24fGIcAH622obLGN81Zv4/N2Sd4ckRXIsODdx3sqnCnUHwAtAKGA0sp/e3/pLsnEJGZwEqgq4hkichEVS0CHgYWA1uBWaq6uarhz3Mum2bcYW0a1eXKbi34eHUWhdapbXxMXmExf1u8nV4x0Yzu3cbpOH7DnULRSVX/Dzitqu9RevNdL3dPoKrjVLW1qoaraqyqvuPav1BVu6hqR1X9Y/Xi/8+5rEXhA8YnxXPoVD5fbjnodBRj/su07/ew79hZnh7ZLejXwa4KdwpFoevrMRHpCUQD7byWyPi9y7u0IKZRXWbYmtrGhxw9XcDkJelc2a0Fl3Zs5nQcv+JOoXhLRBoD/0dpP8IW4AWvpqomu/TkG0JDhNsHxrE8/RB7Dp12Oo4xQOk62Kfzi2wd7GpwZ1LAKap6VFWXqmoHVW2hqm/WRriqsktPvuO2gXGEhggzk61VYZy39/BpPvhhD7cmxNGlpa2DXVUXXJ1DRBoBP6L0ctN/nq+qj3otVTXZfRS+o0XDSK7u3pJPUrN44pou1Amz0SXGOX9dnEZYSAhP2PQy1eLOpaeFlBaJjUBqmc3nWIvCt4xPiufI6QIWbTrgdBQTxNZlHmP+hv3cP6SDrcJYTe6s9xepqk94PYkJOIM7NSO+ST1mrMpgbN/KZmgxxjvOrYPdrEEEk4bYOtjV5dZ9FCJyv4i0FpEm5zavJ6sG68z2LSEhwrjEeFbtPkJ6jtu33hjjMV9tzSF5zxEeH9aFBrYOdrW5UygKgL9SetPcuctOKd4MVV126cn33JIQS3ioMGOVraltaldRcQl//mIrHZrX5/aBcRd+gamQO4XiCUpvumunqu1dm7XhjFuaNajD8Ita8emaLPIKi52OY4LIR6sz2Zl7mqdG2DrYNeXO395m4Iy3g5jANSGpLcfPFrJgw36no5ggcfhUPi9/tYPEdk24ukdLp+P4PXcu2hUD60RkCZB/bqcNjzXuurhDEzo0r8+M5AxuGhDrdBwT4E7mFXLXu8mcyi/k16N72DrYHuBOi+Jz4I/A99jwWFMNIsL4xHhS9x5l24ETTscxASyvsJj73kth2/6TvHnHAHrG2GeBJ1ywReGaCNCYGrmpfywvLE5jxqoMfje2p9NxTAAqKi7h4RlrSd5zhFdu78fQri2cjhQwKlsKdZbr60YR2VB+q72IJhA0rh/Bdb1a89mafZwpKHI6jgkwJSXKk59u4KutB/nd2J6M6WNTiHtSZS2Kx1xfR9VGEBP4JiTF89nafcxbn81tA+OdjmMChKry+wVb+Neaffz8mi7ceXFbpyMFnApbFKp6bohKfVXdW3YD2tdOPBNIBrRtTJeWDWz6ceNRr32Tzrsr9jBxcHseusIGsniDO53Zs0Tkl1Kqroi8Bjzv7WDVYXdm+7Zzndrrs46zaZ/9G5mae3/lHl78cjs39Y/l2ZHdbYSTl7hTKJKAOEpHPa0GsoFB3gxVXTbqyffd0D+WyPAQplurwtTQnHX7+PWczVzdoyV/uamXrVjnRe6ucHcWqAtEArtV1RZDNtUSXTec0b3bMGfdPk7mFV74Bcacx5JtOfxs1nou6dCU18b1szuvvcydv93VlBaKgcBgYJyIzPZqKhPQJlzcljMFxcxZl+10FOOHkncf4cEPU+neuiFv/WgAkeG21om3uVMoJqrqr1W1UFUPqOpYYI63g5nA1Sc2mh6tGzJ9VQaq6nQc40c2Zx9n4rTVxDauy7R7BhIVGe50pKDgTqFYLyKPishs1/Yw8JG3g5nAJSJMuDierftPsC7zmNNxjJ/Yfeg0d01NJioyjA8mJtG0QR2nIwUNdwrFG8AA4B+u7dz3xlTb2L4x1I8ItaGyxi0Hjudxx5RVqMIH9yXRplFdpyMFFXcKxUBVvUtVv3Ft9wCJ3g52joh0F5E3Xa2ZH9fWeY13NagTxpi+MczbkM3xs9apbSp29HQBd76ziuNnC3nv3kQ6Nm/gdKSg406hKBaRjuceiEgHSmeUvSARmSoiOSKyqdz+ESKSJiLpIvJUZcdQ1a2q+iBwK5DgznmNf5iQFE9eYQmfrclyOorxUafyi7h72mr2HjnDlLsSbJI/h7hTKH4OLBGRb0VkKfAN8DM3jz8NGFF2h4iEApOBa4EelI6i6iEivURkfrmthes1Y4DlwNduntf4gZ4x0fSJjbZObXNe+UXFPPBBCpv2Hecf4/tzcYemTkcKWpXOHuv6UO8DdAa6AgJsU9X8yl53jqouE5F25XYnAumqust1jo+Asar6PBXMK6Wqc4G5IrIAmFFB1knAJID4eJtHyF9MSGrLk59uIGXvUQa288ml2I0DiopLeGzmOlakH+bFW/swzBYfclSlLQpVLQbGqGq+qm5Q1fXuFolKxABlF1DOcu07LxEZKiKvisg/gYWVZH0LeA5YExERUcOIpraM6tOaqDph1qlt/kNVefazTSzafIBfj+rBjf1tsSunuXPp6XsReV1ELhOR/ue2GpzzfPfZV3jdQVW/VdVHVfUBVZ1c2YFtCg//Uy8ijBv6x7Bg436Oni5wOo5xmKry/Bfb+Dglk0ev6sy9g23+UV/gTqG4FLgI+B3wd9f2txqcM4vSuaPOiaV0/qgas0kB/dP4pHgKikr41Dq1g94bS3fy1rJd3HVJW346rLPTcYzLBQuFql5xnu3KGpxzNdBZRNqLSARwOzC3Bscrm9VaFH6oW6uGDGjbmBnWqR3UZqzK4IVFaYzt24bfjL7IZoL1IRcsFCLS1NVHsEZEUkXkFRFxa/iBiMwEVgJdRSRLRCaqahHwMLAY2ArMUtXNNXkTZc5nLQo/NSEpnl2HTrNy12GnoxgHzN+QzbOfb+TKbi342y19bCZYHyMX+g1ORL4ElgEfunZNAIaq6jAvZ6u2hIQETUlJcTqGqYK8wmKS/vQ1l3Vuxuvja9IFZvzN0u253PfeavrGNeL9e5OoG2GT/DlFRFJV9X/uV3Onj6KJqv5eVXe7tj8AjTye0AOsReG/IsNDuXlALIs3H+DQqZoOrDP+InXvUR78IJVOLaKYctdAKxI+yp1CsUREbheRENd2K7DA28Gqw/oo/Nu4xHgKi5VPUqxTOxikHTjJvdNW07JhHd6/N5HoujYTrK9yp1A8QOlNbgWu7SPgCRE5KSInvBnOBJdOLRqQ1L4JM5MzKCmxTu1AlnH4DHe+s4q64aF8MDGJ5lE2E6wvc2fUU5SqhqhqmGsLce2LUtWGtRHSXXbpyf9NuLgtGUfOsDz9kNNRjJfknMjjjndWUVBcwgcTE4lrUs/pSOYC3Fo/UERuFJEXReTvInK9lzNVm1168n/DL2pJk/oRTF+11+koxguOnynkR1OTOXQqn2n3JNK5ZZTTkYwb3Bke+w/gQWAjsAl4UEQqvUPamOqqExbKLQmxfLU1h4Mn8pyOYzzoTEER9763ml25p3n7Rwn0jWvkdCTjJndaFJcDw1X1XVV9FxgJDPVqqmqyS0+BYdzAeIpLlFmrMy/8ZOMXCopKePDDNazNOMqr4/oyqFMzpyOZKnCnUKQBZadjjQM2eCdOzdilp8DQrll9BndqxszkDIqtU9vvFZcoT8xax7Ltufz5xt6M6Nna6UimitwpFE2Bra71KL4FtgAtRGSuiHhk6g1jypuQFE/28TyWbs9xOoqpAVXl13M2MX/Dfp4Z2Y1bB8Zd+EXG51S6HoXLr72ewphyhvVoSfOoOkz/IYMru9laBP7qb/9OY/qqDH48tCOThnS88AuMT3JneOzSshtQBNxa5rHPsD6KwBEeGsJtCXEsScsh+9hZp+OYapjy3S4mL9nJuMR4nhze1ek4pgbcHR7bV0ReEJE9wB8onczP51gfRWC5PTEOBT6yTm2/Myslkz8s2Mp1vVrzh+t72kywfq7CQiEiXUTk1yKyFXid0lXpxDXN+Ou1ltAErdjG9RjapTkfr86gqLjE6TjGTYs2HeCpTzdwWedmvHRbX0JtJli/V1mLYhtwFTBaVQer6mtAce3EMqbU+KS2HDyRz9fbrFPbH3yffohHZ66lT1wj/nnnACLC3LpoYXxcZf+KNwEHKJ0U8G0RuYrzL2NqjNdc0bU5rRpGMt3W1PZ56zOPcf/7KbRvVp937x5IvQh3xsoYf1BhoVDVz1T1NqAb8C3wU6CliLwhItfUUj4T5MJCQ7g9MY7vduSScfiM03FMBdJzTnL3u8k0aRDBBxMTaVQvwulIxoPcGfV0WlWnq+ooSte3Xgc85e1g1WGjngLTbQPjEGDmamtV+KKso2e4Y0oyYaEhfDgxiRYNI52OZDysShcQVfWIqv6zhmtme42NegpMraPrclX3lnySkklBkXVq+5JDp/K5851kzhQU8f69ibRtWt/pSMYLrKfJ+IXxSfEcOlXAl1sOOh3FuJzIK+RH7ySz//hZ3r1nIN1b+9SqA8aDrFAYvzCkc3NiG9e16cd9RF5hMfdNS2FHzknevGMAA9o2cTqS8SIrFMYvhIYI4xLj+X7nYXblnnI6TlArLC7hoelrWL33CC/e2pehXVs4Hcl4mRUK4zduSYglLESYmWyd2k4pKVGenL2Br7fl8PuxPRndp43TkUwt8ItCISL1RSRVREY5ncU4p0VUJNdc1JLZqVnkFdq9n7VNVfnd/C18tnYfvxjelTsubut0JFNLvFooRGSqiOSIyKZy+0eISJqIpIuIO0NtfwnM8k5K40/GJ7bl6JlCFm064HSUoPPK1zuY9v0e7hvcnp8MtZlgg4m3WxTTgBFld4hIKDAZuBboAYwTkR4i0ktE5pfbWojIMErXwLDhLoZLOzalXdN6zLA7tWvVtBW7efmrHdwyIJZnr+tuk/wFGa/eY6+qy0SkXbndiUC6qu4CEJGPgLGq+jzwP5eWROQKoD6lReWsiCxUVRtMH6RCXJ3az3+xjR0HT9K5ZZTTkQLeZ2uz+O28LVzToyXP39jLikQQcqKPIobSmWjPyXLtOy9VfVZVHwdmAG9XVCREZJKIpIhISm5urifzGh9z84BYIkJDbP6nWvD11oP8/JMNXNqxKa+O60dYqF90axoPc+Jf/Xy/jlxwYWRVnaaq8yv587eA54A1ERE2z0wga9qgDiN6tuJfa7I4W2Cd2t6yatdhfjJ9DT3bNOStHyUQGR7qdCTjECcKRRZQduHcWCDbEwe2KTyCx/ikeE7kFTF/g0f+65hyNu07zn3vpRDXpB7v3pNIgzo2E2wwc6JQrAY6i0h7EYkAbgfmeuLANilg8Ehq34SOzeszw+6p8Lhduae4a2oyDeuG88HERJrUtxZ6sPP28NiZwEqgq4hkichEVS0CHgYWU7qk6ixV3eyJ81mLIniICOOT2rI24xhbsk84HSdg7D9+ljvfSQbgg4mJtI6u63Ai4wu8WihUdZyqtlbVcFWNVdV3XPsXqmoXVe2oqn/01PmsRRFcbuofQ52wEGYk2/xPnnDkdAF3TFnFibOFvHdvIh2aN3A6kvERATWEwVoUwaVRvQiu692az9dmczq/yOk4fu1UfhF3v5tM1tGzTLkrgZ4x9jNk/r+AKhTWogg+E5LiOZVfxNz11qldXXmFxUx6P4XN2SeYPL4/SR2aOh3J+JiAKhTWogg+/eMb061VlN2pXU1FxSU8OnMt3+88zN9u6c2wHi2djmR8UEAVChN8Sju149m47zgbso45HcevqCpP/2sj/95ykN+O7sEN/WKdjmR8VEAVCrv0FJyu7xdD3fBQa1VUgaryxwVb+SQ1i8eu6szdg9o7Hcn4sIAqFHbpKTg1jAxnTJ82zF2fzYm8Qqfj+IV/fLuTKct3c/el7Xh8WGen4xgfF1CFwgSv8UnxnCkoZs7afU5H8Xkf/rCXvy5O44Z+Mfx6VA+b5M9cUEAVCrv0FLx6x0bTM6Yh01dloHrBqcOC1rz12fzfnE1c1a0FL9zcm5AQKxLmwgKqUNilp+AlIoxPbMu2AydZk3HM6Tg+6du0HH768ToGtm3C5An9CbeZYI2b7H+KCRhj+rahQZ0w69Q+j9S9R3jww1S6tIxiyt02E6ypGisUJmA0qBPG2L5tmL8hm+NnrFM7r7CYBRv288AHKYx7axWto+vy/sREGkaGOx3N+JmAmjtYREYDozt16uR0FOOQCUltmb4qg0/XZHHv4OAb8llYXMLy9EPMW5fN4s0HOF1QTPOoOky4OJ4HhnSkWYM6Tkc0fkgCseMvISFBU1JSnI5hHHL95BWcyi/iy58OCYoRPSUlSsreo8xdv4+FGw9w5HQBDSPDGNmrNWP6tCGpQ1NCrdPauEFEUlU1ofz+gGpRGAOlQ2WfnL2B5N1HAnbeIlVlc/YJ5q7PZv76bLKP5xEZHsLVPVoxpk8bhnRpRp0w64cwnmGFwgSc0b3b8Pv5W5iRnBFwhWJX7inmrs9m7vpsduWeJixEuLxLc355bTeGdW9JfVuJzniB/a8yAaduRCg39Y9lxqoMfjO6wO9XaNt//Czz1+9n7vpsNu47jkjpCn/3De7AtT1b0djP35/xfVYoTEAanxTPtO/3MDs1k0lDOjodp8qOni5g4ab9zF2XTfKeI6iW3lT4q+u6M6p3G1pFRzod0QSRgCoUNurJnNOlZRQD2zVmZnIm9w3u4Bd3IJ/OL+LLLQeZuz6bZdtzKSpROjavz+NXdWFM3za0b1bf6YgmSAVUoVDVecC8hISE+53OYpw3Pimen368npW7DjOoUzOn45xXflExS9NymbM+m6+3HiSvsIQ20ZFMHNyeMX3b0KN1w6AYuWV8W0AVCmPKurZna56bt4UZqzJ8qlAUlyg/7DrMnHX7WLTpACfyimhSP4JbBsQxpm8bBsQ39osWkAkeVihMwIoMD+Xm/rFM+34POSfzaBHl3HV9VWVt5jHmrstmwcb95J7Mp35EKMN7lg5nHdSpmc29ZHyWFQoT0MYlxTNl+W4+ScnioStqv+9q+8GTzFm3j3nr95Nx5AwRYSFc0bU5Y/vGcGW3FjbnkvELVihMQOvYvAGXdGjKzOQMfnx5x1q5pJN55Axz12czb3022w6cJERgUKdmPHJlJ4b3bGVzLRm/4/OFQkSGAr8HNgMfqeq3TuYx/mfCxfE8PGMty3bkMrRrC6+cI/dkPgs2lN4Id26a8wFtG/PcmIsY2as1zaNsjiXjv7xaKERkKjAKyFHVnmX2jwBeAUKBKar650oOo8ApIBLI8mJcE6Cu6dGKZg0imLEqw6OF4vjZQhZvPsC89dmsSD9EiUK3VlE8OaIro3u3Ia5JPY+dyxgnebtFMQ14HXj/3A4RCQUmA1dT+sG/WkTmUlo0ni/3+nuB71R1qYi0BF4EJng5swkwEWEh3Dwgjre/28WB43k1ulktr7CYr7fmMGfdPr5Ny6WguIT4JvX4ydBOjOnbhi4tozyY3Bjf4NVCoarLRKRdud2JQLqq7gIQkY+Asar6PKWtj4ocBSpsv4vIJGASQHx8fE1imwA0PjGeN5fu5OPVmTw2rHOVXntu6u6567L5d7mpu8f2jaFPbLTd62ACmhN9FDFAZpnHWUBSRU8WkRuB4UAjSlsn56WqbwFvQek0454IagJHfNN6XNa5GR+tzuChKzoSdoGhqCUlyuo9R5i7PpuFG/dz9EwhDSPDGNW7DWP72tTdJrg4USjO99NV4Qe7qv4L+JdbB7YpPEwlJiS15cEPU/k2LZdhPVr+z5+Xnbp73vps9rum7h7WvSVj+8bY1N0maDlRKLKAuDKPY4FsTxzYpvAwlbmqewtaRNVhRnLGfxWKiqbufsqm7jYGcKZQrAY6i0h7YB9wOzDeEwe2FoWpTHhoCLcNjOP1Jemk7j3Cmr3HbOpuY9zg7eGxM4GhQDMRyQJ+o6rviMjDwGJKRzpNVdXN3sxhzDm3J8YzeUk6N72xErCpu41xh62ZbYLOuyt2c+JskU3dbUw5tma2MS73DGrvdARj/EpATVcpIqNF5K3jx487HcUYYwJGQBUKVZ2nqpOio6OdjmKMMQEjoAqFtSiMMcbzAqpQWIvCGGM8L6AKhTHGGM8LqEJhl56MMcbzAqpQ2KUnY4zxvIAqFMYYYzwvIO/MFpFcYG81X94MOOTBOP7A3nNwsPcc+Gr6ftuqavPyOwOyUNSEiKSc7xb2QGbvOTjYew583nq/dunJGGNMpaxQGGOMqZQViv/1ltMBHGDvOTjYew58Xnm/1kdhjDGmUtaiMMYYUykrFMYYYyplhcJFRKaKSI6IbHI6S20RkTgRWSIiW0Vks4g85nQmbxKRSBFJFpH1rvf7nNOZaouIhIrIWhGZ73SW2iAie0Rko4isE5GgWO5SRBqJyGwR2eb6mb7EY8e2PopSIjIEOAW8r6o9nc5TG0SkNdBaVdeISBSQClyvqlscjuYVIiJAfVU9JSLhwHLgMVX9weFoXiciTwAJQENVHeV0Hm8TkT1AgqoGzc12IvIe8J2qThGRCKCeqh7zxLGtReGiqsuAI07nqE2qul9V17i+PwlsBWKcTeU9WuqU62G4awv435REJBa4DpjidBbjHSLSEBgCvAOgqgWeKhJghcK4iEg7oB+wyuEoXuW6BLMOyAG+VNWAfr8uLwNPAiUO56hNCvxbRFJFZJLTYWpBByAXeNd1iXGKiNT31MGtUBhEpAHwKfC4qp5wOo83qWqxqvYFYoFEEQnoy4wiMgrIUdVUp7PUskGq2h+4FnjIdWk5kIUB/YE3VLUfcBp4ylMHt0IR5FzX6j8Fpqvqv5zOU1tczfJvgRHOJvG6QcAY1zX7j4ArReRDZyN5n6pmu77mAJ8Bic4m8rosIKtMC3k2pYXDI6xQBDFX5+47wFZVfdHpPN4mIs1FpJHr+7rAMGCbo6G8TFWfVtVYVW0H3A58o6p3OBzLq0SkvmtwBq7LL9cAAT2aUVUPAJki0tW16yrAY4NSwjx1IH8nIjOBoUAzEckCfqOq7zibyusGAXcCG13X7QGeUdWFzkXyqtbAeyISSukvSbNUNSiGiwaZlsBnpb8HEQbMUNVFzkaqFY8A010jnnYB93jqwDY81hhjTKXs0pMxxphKWaEwxhhTKSsUxhhjKmWFwhhjTKWsUBhjjKmUFQpj3CAirUTkIxHZKSJbRGShiHQRkXYictY1S+kWEXlTRDz+cyUivxWRn3v6uMa4wwqFMRfgujHxM+BbVe2oqj2AZygdrw+w0zUtSG+gB3B9ZcfyRiExxpvsP6wxF3YFUKiqb57boarrVPW7sk9S1SLge6BT2f2uVsdWEfkHsAaIE5E3RCSl/LoYrnUUnhORNa71FLqVDyMi94vIF667y43xOisUxlxYT0rX6qiUiNSjdOqEjef5466UrnXST1X3As+qagKlrZDLRaR3mececk1o9wbwX5ebRORhYDSl64acrda7MaaKbAoPY2quo2sKFAXmqOoX53nO3nILJN3qmv46jNKpRXoAG1x/dm5yxlTgxjKvuZPSyd+uV9VCD+Y3plJWKIy5sM3AzZX8+bk+isqcPveNiLSntKUwUFWPisg0ILLMc/NdX4v575/RTUBfSqdI3+1OcGM8wS49GXNh3wB1ROT+cztEZKCIXF7N4zWktHAcF5GWlK6Z4I61wAPAXBFpU81zG1NlViiMuQAtnTnzBuBq1/DYzcBvgexqHm89pR/6m4GpwIoqvHY5pa2RBSLSrDrnN6aqbPZYY4wxlbIWhTHGmEpZoTDGGFMpKxTGGGMqZYXCGGNMpaxQGGOMqZQVCmOMMZWyQmGMMaZS/w98bAx+IXj+7AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import tensorly as tl\n",
    "import tensorly.decomposition as tldec\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "\n",
    "tensor = tl.tensor(np.arange(24).reshape((3, 4, 2)))\n",
    "# tensor = tl.random.cp_tensor(shape=(3,3,3), rank=3, full=True)\n",
    "print(tensor)\n",
    "rank_range = [1, 2, 3, 4, 5, 6]\n",
    "error_rec = []\n",
    "for r in rank_range:\n",
    "    factors = tldec.parafac(tensor.astype(\"f\"), rank=r)\n",
    "    error_rec.append(np.linalg.norm(tl.kruskal_to_tensor(factors) - tensor))\n",
    "plt.semilogy(rank_range, error_rec)\n",
    "plt.xlabel(\"CP rank\")\n",
    "plt.ylabel(\"Approximation error\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Tucker decomposition\n",
    "\n",
    "Next attempt is the decomposition proposed by ([Tucker, 1963](https://link.springer.com/content/pdf/10.1007/BF02289464.pdf)) in [Psychometrika](https://link.springer.com/journal/11336):\n",
    "\n",
    "$$\n",
    "a_{ijk} = \\sum_{\\alpha_1,\\alpha_2,\\alpha_3=1}^{r_1,r_2,r_3}g_{\\alpha_1\\alpha_2\\alpha_3} u_{i\\alpha_1} v_{j\\alpha_2} w_{k\\alpha_3}.\n",
    "$$\n",
    "\n",
    "Here we have several different ranks. Minimal possible $r_1,r_2,r_3$ are called **Tucker ranks**.\n",
    "\n",
    "**Properties**:\n",
    "\n",
    "* For a $d$ dimensional tensor memory is <font color='red'>$r^d$</font> $+ nrd$. Still **curse of dimensionality**\n",
    "* Stable SVD-based algorithm:\n",
    "    1. $U =$ principal components of the unfolding `A.reshape(n1, n2*n3)`\n",
    "    2. $V =$ principal components of the unfolding `A.transpose([1,0,2]).reshape(n2, n1*n3)`\n",
    "    3. $W =$ principal components of the unfolding `A.transpose([2,0,1]).reshape(n3, n1*n2)`\n",
    "    4. $g_{\\alpha_1\\alpha_2\\alpha_3} = \\sum_{i,j,k=1}^{n_1,n_2,n_3} a_{ijk} u_{i\\alpha_1} v_{j\\alpha_2} w_{k\\alpha_3}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Shape of core = (3, 4, 2)\n",
      "Shape of factors = (3, 3)\n",
      "Shape of factors = (4, 4)\n",
      "Shape of factors = (2, 2)\n",
      "Approximation error = 3.181718859384126e-14\n"
     ]
    }
   ],
   "source": [
    "tensor = tl.tensor(np.arange(24).reshape((3, 4, 2)))\n",
    "core, factors = tldec.tucker(tensor, ranks=[3, 4, 2])\n",
    "print(\"Shape of core = {}\".format(core.shape))\n",
    "for f in factors:\n",
    "    print(\"Shape of factors = {}\".format(f.shape))\n",
    "print(\"Approximation error = {}\".format(np.linalg.norm(tensor - tl.tucker_to_tensor(core, factors))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Application in recommender systems [(Frolov, Oseledets 2016)](https://dl.acm.org/citation.cfm?id=2959170)\n",
    "\n",
    "- Represent User-Item **matrix** as binary User-Item-Rating **tensor**\n",
    "- After Tucker decomposition of this tensor we have three factors\n",
    "- Two of them represent projections on the space of ratings and space of items\n",
    "- This approach takes into account negative feedback to items from some users"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## CP and Tucker decompositions implementations\n",
    "\n",
    "- Matlab: [Tensorlab](https://www.tensorlab.net/) and [Tensor Toolbox](https://www.sandia.gov/~tgkolda/TensorToolbox/index-2.5.html)\n",
    "- Python: [TensorLy](http://tensorly.org/stable/home.html) and [Scikit-tensor](https://github.com/mnick/scikit-tensor)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Tensor Train decomposition\n",
    "\n",
    "* Calculation of the canonical decomposition is unstable\n",
    "* Tucker decomposition suffers from the curse of dimensionality\n",
    "\n",
    "Tensor Train (**TT**) decomposition ([Oseledets, Tyrtyshnikov 2009](http://www.mat.uniroma2.it/~tvmsscho/papers/Tyrtyshnikov4.pdf) and [Oseledets, 2011](https://epubs.siam.org/doi/abs/10.1137/090752286)) is both stable and contains linear in $d$ number of parameters:\n",
    "\n",
    "$$\n",
    "a_{i_1 i_2 \\dots i_d} = \\sum_{\\alpha_1,\\dots,\\alpha_{d-1}} \n",
    "g_{i_1\\alpha_1} g_{\\alpha_1 i_2\\alpha_2}\\dots g_{\\alpha_{d-2} i_{d-1}\\alpha_{d-1}} g_{\\alpha_{d-1} i_{d}}\n",
    "$$\n",
    "\n",
    "or in the matrix form\n",
    "\n",
    "$$\n",
    "    a_{i_1 i_2 \\dots i_d} = G_1 (i_1)G_2 (i_2)\\dots G_d(i_d)\n",
    "$$\n",
    "\n",
    "* The storage is $\\mathcal{O}(dnr^2)$ \n",
    "* Stable TT-SVD algorithm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "**Example**\n",
    "$$a_{i_1\\dots i_d} = i_1 + \\dots + i_d$$\n",
    "Canonical rank is $d$. At the same time TT-ranks are $2$:\n",
    "$$\n",
    "i_1 + \\dots + i_d = \\begin{pmatrix} i_1 & 1 \\end{pmatrix} \n",
    "\\begin{pmatrix} 1 & 0 \\\\ i_2 & 1 \\end{pmatrix}\n",
    "\\dots\n",
    "\\begin{pmatrix} 1 & 0 \\\\ i_{d-1} & 1 \\end{pmatrix}\n",
    "\\begin{pmatrix} 1  \\\\  i_d \\end{pmatrix}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Implementations\n",
    "\n",
    "- Matlab: [TT Toolbox](https://github.com/oseledets/TT-Toolbox)\n",
    "- Python: [ttpy](https://github.com/oseledets/ttpy)\n",
    "- TensorFlow: [t3f](https://github.com/Bihaqo/t3f)\n",
    "- PyTorch: [tt-pytorch](https://github.com/KhrulkovV/tt-pytorch) and [tntorch](https://github.com/rballester/tntorch)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Using TT in riemannien optimization (example is from t3f examples)\n",
    "\n",
    "- Given some TT tensor $A$ with large tt-ranks\n",
    "- We would like to find a tensor $X$ (with small prescribed tt-ranks $r$) which is closest to $A$ (in the sense of Frobenius norm):\n",
    "\\begin{equation*}\n",
    "\\begin{aligned}\n",
    "& \\underset{X}{\\text{minimize}} \n",
    "& & \\frac{1}{2}\\|X - A\\|_F^2 \\\\\n",
    "& \\text{subject to} \n",
    "& & \\text{tt_rank}(X) = r\n",
    "\\end{aligned}\n",
    "\\end{equation*}\n",
    "\n",
    "- It is known that the set of TT tensors with elementwise fixed TT ranks forms a manifold. \n",
    "- Thus we can solve this problem using the so called Riemannian gradient descent. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Riemannian gradient descent\n",
    "\n",
    "- Given some functional $F$ on a manifold $\\mathcal{M}$ Riemannian gradient descent is defined as\n",
    "\n",
    "$$\\hat{x}_{k+1} = x_{k} - \\alpha P_{T_{x_k}\\mathcal{M}} \\nabla F(x_k),$$\n",
    "\n",
    "$$x_{k+1} = \\mathcal{R}(\\hat{x}_{k+1})$$\n",
    "\n",
    "with $P_{T_{x_k}\\mathcal{M}}$ being the projection onto the **tangent space** of $\\mathcal{M}$ at the point $x_k$ and $\\mathcal{R}$ being a **retraction** - an operation which projects points to the manifold, and $\\alpha$ is the learning rate.\n",
    "\n",
    "- We can implement this in `t3f` using the `t3f.riemannian` module. As a retraction it is convenient to use the rounding method (`t3f.round`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "import t3f\n",
    "import tensorflow as tf\n",
    "tf.set_random_seed(0)\n",
    "np.random.seed(0)\n",
    "sess = tf.InteractiveSession()\n",
    "\n",
    "# Initialize A randomly, with large tt-ranks\n",
    "shape = 10 * [2]\n",
    "init_A = t3f.random_tensor(shape, tt_rank=16)\n",
    "A = t3f.get_variable('A', initializer=init_A, trainable=False)\n",
    "\n",
    "# Create an X variable and compute the gradient of the functional. Note that it is simply X - A.\n",
    "\n",
    "init_X = t3f.random_tensor(shape, tt_rank=2)\n",
    "X = t3f.get_variable('X', initializer=init_X)\n",
    "\n",
    "gradF = X - A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "# Let us compute the projection of the gradient onto the tangent space at X\n",
    "riemannian_grad = t3f.riemannian.project(gradF, X)\n",
    "\n",
    "# Compute the update by subtracting the Riemannian gradient\n",
    "# and retracting back to the manifold\n",
    "alpha = 1.0\n",
    "\n",
    "train_step = t3f.assign(X, t3f.round(X - alpha * riemannian_grad, max_tt_rank=2))\n",
    "\n",
    "# let us also compute the value of the functional\n",
    "# to see if it is decreasing\n",
    "F = 0.5 * t3f.frobenius_norm_squared(X - A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "81.62205\n",
      "58.53462\n",
      "56.27\n",
      "56.08319\n",
      "51.73249\n",
      "50.77673\n",
      "50.77672\n",
      "50.77671\n",
      "50.77672\n",
      "50.77672\n"
     ]
    }
   ],
   "source": [
    "sess.run(tf.global_variables_initializer())\n",
    "\n",
    "log = []\n",
    "for i in range(100):\n",
    "    F_v, _ = sess.run([F, train_step.op])\n",
    "    if i % 10 == 0:\n",
    "        print (F_v)\n",
    "    log.append(F_v)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "It is intructive to compare the obtained result with the quasioptimum delivered by the TT-round procedure. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "52.40742\n"
     ]
    }
   ],
   "source": [
    "quasi_sol = t3f.round(A, max_tt_rank=2)\n",
    "\n",
    "val = sess.run(0.5 * t3f.frobenius_norm_squared(quasi_sol - A))\n",
    "print (val)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "We see that the value is slightly bigger than the exact minimum, but TT-round is faster and cheaper to compute, so it is often used in practice."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0xd261fa4a8>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZkAAAEKCAYAAADAVygjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzt3Xl8VOXZ//HPlT2EEEjYA0pQFlkTJKCCSute0Fp3qrVuULUu1bZP7VOtS+vz68+f2qq1ta7Ypxaq4galalXQSmtZFAREcQEhgIBBwhqyXb8/ZjINIcskzGTIzPf9eo3hnDnnPteZE+fKfd/n3Le5OyIiItGQFOsAREQkfinJiIhI1CjJiIhI1CjJiIhI1CjJiIhI1CjJiIhI1CjJiIhI1CjJiIhI1CjJiIhI1KTEOoBY69q1q/fr1y/WYYiItCuLFy/+0t27NbddwieZfv36sWjRoliHISLSrpjZ5+Fsp+YyERGJGiUZERGJGiUZERGJmoTvkxGJhMrKSkpKSigvL491KCIRlZGRQZ8+fUhNTW3V/koyIhFQUlJCdnY2/fr1w8xiHY5IRLg7paWllJSUUFBQ0Koy1FwmEgHl5eXk5eUpwUhcMTPy8vIOqIauJCMSIUowEo8O9PdaSaaV/vr+RqbNXx3rMEREDmpKMq302spN/OGtz2IdhkhIcnIyhYWFDBs2jNNPP51t27YBsGHDBs4555wYR7e/n//857z22muxDqNBHTt2BA78s/vNb37D7t27m91u3rx5TJo0qdXHOVDbtm3jd7/7XVTKVpJppUE9s9lYVk7Z7spYhyICQGZmJkuWLGH58uXk5uby4IMPAtC7d2+effbZGEe3vzvuuIMTTzyxzY5XXV3d4n0O9LMLN8nEmpLMQWhQj2wAPtq0I8aRiOzv6KOPZv369QCsWbOGYcOGAYEv2h//+McUFxczYsQI/vCHPwCBv6SPP/54zjvvPAYOHMhNN93EU089xZgxYxg+fDiffvopALNmzWLs2LEUFRVx4oknsmnTJgBuu+02LrvsMiZMmED//v25//77Q8c+4ogjmDJlCkOHDuXkk09mz549AFxyySWhL/A77riD4uJihg0bxtSpU3F3ACZMmMBPfvITxowZw8CBA/nHP/6x37nW1NRw9dVXM3ToUCZNmsQ3vvGNULn9+vXjjjvuYPz48TzzzDM88sgjFBcXM3LkSM4+++xQAli9ejVHH300xcXF3HLLLaGyw/3sJkyYwDnnnMPgwYO58MILcXfuv/9+NmzYwNe+9jW+9rWv7Rf3yy+/zODBgxk/fjzPPfdcaP2uXbu47LLLKC4upqioiBdffBGAFStWMGbMGAoLCxkxYgQff/wxAH/84x8ZMWIEI0eO5Dvf+Q4AW7Zs4eyzz6a4uJji4mLmz5/f5HW66aab+PTTTyksLOTHP/5xc79eLaJbmFtpUM//JJkxBbkxjkYOJrfPWsEHG7ZHtMwhvTtx6+lDw9q2urqa119/ncsvv3y/9x577DFycnJYuHAhe/fuZdy4cZx88skALF26lJUrV5Kbm0v//v254oorWLBgAffddx8PPPAAv/nNbxg/fjzvvPMOZsajjz7KXXfdxT333APAhx9+yNy5c9mxYweDBg3iqquuAuDjjz9m+vTpPPLII5x33nnMnDmTiy66aJ+4rrnmGn7+858D8J3vfIfZs2dz+umnA1BVVcWCBQuYM2cOt99++35NbM899xxr1qxh2bJlbN68mSOOOILLLrss9H5GRgZvv/02AKWlpUyZMgWAm2++mccee4xrr72W66+/nquuuoqLL744VANsyWf33nvvsWLFCnr37s24ceOYP38+1113Hffeey9z586la9eu+5RVXl7OlClTeOONNzj88MM5//zzQ+/deeedfP3rX+fxxx9n27ZtjBkzhhNPPJGHHnqI66+/ngsvvJCKigqqq6tZsWIFd955J/Pnz6dr165s3boVgOuvv54bbriB8ePHs3btWk455RRWrlzZ6HX61a9+xfLly1myZEmD534glGRaqVdOBtkZKXz0RWS/TERaa8+ePRQWFrJmzRqOPPJITjrppP22efXVV3n//fdDf+mXlZXx8ccfk5aWRnFxMb169QLgsMMOC32BDh8+nLlz5wKB54HOP/98Nm7cSEVFxT7PTkycOJH09HTS09Pp3r17qJZTUFBAYWEhAEceeSRr1qzZL665c+dy1113sXv3brZu3crQoUNDSeass85qct+3336bc889l6SkJHr27LlfraHuF/jy5cu5+eab2bZtGzt37uSUU04BYP78+cycORMIJLmf/OQnLfrsxowZQ58+fQBC12D8+PH7lVHrww8/pKCggAEDBgBw0UUX8fDDD4eO89JLL3H33XcDgYS0du1ajj76aO68805KSko466yzGDBgAG+88QbnnHNOKInl5gb+4H3ttdf44IMPQsfbvn07O3YEWl0au07REpdJxsz6Az8Dctw9Kj2eZsbgntl89IWay2Rf4dY4Iq22T6asrIxJkybx4IMPct111+2zjbvzwAMPhL5ca82bN4/09PTQclJSUmg5KSmJqqoqAK699lpuvPFGzjjjDObNm8dtt90W2qfu/snJyaF96q+vbS6rVV5eztVXX82iRYvo27cvt9122z7PZdTuX7fM+ufUlKysrNC/L7nkEl544QVGjhzJtGnTmDdvXui95m7VDfezayzO+ho7nrszc+ZMBg0atM/6I444grFjx/LXv/6VU045hUcffRR3b7Ccmpoa/vWvf5GZmbnfe62J9UBEtU/GzG4wsxVmttzMpptZRivLedzMNpvZ8gbeO9XMPjKzT8zsJgB3/8zd928riLCBPbL58Isdzf6Si7SlnJwc7r//fu6++24qK/e9MeWUU07h97//fWj9qlWr2LVrV9hll5WVkZ+fD8CTTz4ZkXhrE0rXrl3ZuXNnizvax48fz8yZM6mpqWHTpk37JI76duzYQa9evaisrOSpp54KrR83bhwzZswA2Gd9Xa357LKzs0M1iLoGDx7M6tWrQ31d06dP3+c4DzzwQOh75b333gPgs88+o3///lx33XWcccYZvP/++5xwwgk8/fTTlJaWAoSay04++WR++9vfhspsrhmssTgjIWpJxszygeuA0e4+DEgGLqi3TXczy6637vAGipsGnNrAMZKBB4HTgCHAZDMbEpETCMPgntnsKK/ii+0ar0oOLkVFRYwcOTL0xVnriiuuYMiQIYwaNYphw4bxve99r0V/yd52222ce+65HHvssfv1M7RW586dmTJlCsOHD+fMM8+kuLi4RfufffbZ9OnTJ3Q+Y8eOJScnp8Ftf/GLXzB27FhOOukkBg8eHFp/33338eCDD1JcXExZWVmD+7bms5s6dSqnnXbafk14GRkZPPzww0ycOJHx48dz6KGHht675ZZbqKysZMSIEQwbNix0I8Jf/vIXhg0bRmFhIR9++CEXX3wxQ4cO5Wc/+xnHH388I0eO5MYbbwTg/vvvZ9GiRYwYMYIhQ4bw0EMPNRlnXl4e48aNY9iwYRHv+Mfdo/IC8oF1QC6BZrnZwMn1tjkXeAPICC5PAeY0Ul4/YHm9dUcDr9RZ/inw0zrLzzYX55FHHumt9e/PSv3Qn8z2Nz7c1OoyJD588MEHsQ4hoe3YscPd3b/88kvv37+/b9y4McYRxZeGfr+BRR5GLohaTcbd1wN3A2uBjUCZu79ab5tngJeBGWZ2IXAZcF4LDlObyGqVAPlmlmdmDwFFZvbThnY0s9PN7OHG/moJR+g2ZvXLiMTUpEmTKCws5Nhjj+WWW26hZ8+esQ5JgqLW8W9mXYBvAgXANuAZM7vI3f9Udzt3v8vMZgC/Bw5z950tOUwD69zdS4Erm9rR3WcBs0aPHj2lBcfbR06HVHp2ylCSEYmxpvphJLai2fF/IrDa3be4eyXwHHBM/Y3M7FhgGPA8cGsLj1EC9K2z3AfY0LpwW2eQ7jATEWlUNJPMWuAoM+tggXvsTgBW1t3AzIqARwjUeC4Fcs3sly04xkJggJkVmFkagRsLXopI9GEa3DObT7bspKq6pi0PKyLSLkSzT+bfwLPAu8Cy4LEerrdZB+Bcd//U3WuA7wKf1y/LzKYD/wIGmVmJmV0ePEYVcA3wCoEE9rS7r4jSKTVoYI9sKqpqWFMa/m2gIiKJIqoPY7r7rTTRBObu8+stVxKo2dTfbnITZcwB5hxAmAekdniZD7/YweHds5vZWkQksWiAzAN0ePeOJCcZq9QvIzFUWlpKYWEhhYWF9OzZk/z8fIYPH05hYSFDhgwJravdpqKiImax1g6jD7Bx48b9hri//vrryc/Pp6bmP03Qs2fP5tZbW9plKwcDJZkDlJGaTL+8DnyoJCMxlJeXx5IlS1iyZAlXXnklN9xwA8uWLWPJkiV88MEHoXW126SlpTVaVrSHGanr3nvvDQ1YCYHhUJ5//nn69u3LW2+9FVo/ceJEXnrppXYxbL7sS0kmAgb1zNaQ/9Ku3XbbbUydOpWTTz6Ziy++mPLyci699FKGDx9OUVFRaIDMadOmcc0114T2mzRpUuj24Y4dO/Kzn/2MkSNHctRRR4UGXmxsGH2AmTNncuqp/xnMY+7cuQwbNoyrrrpqn6FWzIwJEyYwe/bsaH0EEiVKMhFweLeOrNu6m4oq3WEm7dfixYt58cUX+fOf/xwa7n7ZsmVMnz6d7373u/sMWtmQXbt2cdRRR7F06VKOO+44Hnkk0L1aO4z+woUL93lIcvXq1XTp0mWfARunT5/O5MmT+da3vsXs2bP3GXtt9OjRDc4nIwe3uByFua0VdMuixmHt1t0c3r1j8ztI3Js3bx5vvvlmaLm2Saj2ixfg+OOPZ8KECdxzzz3s3Bl4BrlXr15MnTqVWbNm8e6774a2vfHGG8nOju6NJWeccUZo1N63336ba6+9FggM5njooYeyatWqJvdPS0sL9a8ceeSR/P3vfwcaH0Z/48aNdOvWLbR/RUUFc+bM4de//jXZ2dmMHTuWV199lYkTJwLQvXt3Nmxo08fgJAKUZCKgX15gKPHVX+5SkhEgMKPjhAkT9lvfUOf1D3/4w/3WnX766aH5VKLhwQcfDCW8OXMCN2fWHRLfGxlZPCUlZZ8O+bq1m9TU1NCw8/WHkG9oOPrMzMx99n/55ZcpKytj+PDhAOzevZsOHTqEkkx5eXmDQ9fLwU3NZRFQ0DXwP+eaL/WsjLQP3//+90M3AfTu3Xu/94877rjQkPerVq1i7dq1DBo0iH79+rFkyRJqampYt24dCxYsaPZYjQ2jP3DgwH0mIZs+fTqPPvooa9asYc2aNaxevZpXX3011Nm/atWq0FTI0n4oyURA5w5pdOmQymdKMhInrr76aqqrqxk+fDjnn38+06ZNIz09nXHjxlFQUMDw4cP50Y9+xKhRo5otq7Fh9LOysjjssMP45JNP2L17N6+88kqo1lL7/vjx45k1axYQuCmg7vvSPlhj1eJEMXr0aF+0aNEBl3PW7+aTlpLEjKlHRyAqaW9WrlzJEUccEesw2p3nn3+exYsX88tfNj2a1KZNm/j2t7/N66+/3kaRSV0N/X6b2WJ3H93cvqrJREhB146s+VL38Iu0xLe+9S369evX7HZr167lnnvuiX5AEnHq+I+Qgq4dmPluObsrquiQpo9VJFxXXHFFs9u0dLZMOXioJhMhBV0Dd5WpNpO4Er3pWeLTgf5eK8lESL+uHYDAbcySeDIyMigtLVWikbji7pSWlpKRkdHqMtSuEyG1z8poyP/E1KdPH0pKStiyZUusQxGJqIyMDPr06dPq/ZVkIiQrPYUendL5bIuSTCJKTU2loKAg1mGIHHTUXBZBBV2zVJMREalDSSaCCrpmqU9GRKQOJZkIKuiaxdZdFZTtrmx+YxGRBKAkE0GhgTLVZCYiAijJRFT/brWjMe+McSQiIgcHJZkI6pvbgSSD1XogU0QEUJKJqPSUZPK7ZKrzX0QkSEkmwgIDZSrJiIiAkkzEFeR1UE1GRCRISSbCeuRksHNvFXsqqmMdiohIzCnJRFheVhoApbv2xjgSEZHYU5KJsC4dAklm666KGEciIhJ7SjIRltdRSUZEpJaSTITlZqUDSjIiItDEUP9mtgxoaAYmA9zdR0QtqnYsN0s1GRGRWk3NJzOpzaKII50yUkhJMkqVZEREGk8y7v55WwYSL8yMLllpfKUkIyLSfJ+MmR1lZgvNbKeZVZhZtZltb4vg2qu8rDTVZERECK/j/7fAZOBjIBO4AnggmkG1d7lZaeqTEREhzLvL3P0TINndq939CeBr0Q2rfeuiJCMiAjTd8V9rt5mlAUvM7C5gI5AV3bDatzwlGRERILyazHeC210D7AL6AmdHM6j2LjcrjbI9lVRW18Q6FBGRmGq2JlPnLrNy4PbohhMfascv+2p3Bd2zM2IcjYhI7IRzd9k4M/u7ma0ys89qX20RXHvVpTbJ7KqMcSQiIrEVTp/MY8ANwGJA49eHIXefkZizYxuMiEgMhZNkytz9b1GPJI7kafwyEREgvCQz18z+H/AcEJokxd3fjVpU7ZzGLxMRCQgnyYwN/hxdZ50DX498OPGhc4dUQElGRCScu8v04GULpSYnkZOZqiQjIgkvnLvLcszsXjNbFHzdY2Y5bRFce6bxy0REwnsY83FgB3Be8LUdeCKaQR0oM+tvZo+Z2bOxiiE3K42tO5VkRCSxhZNkDnP3W939s+DrdqB/czuZ2SAzW1Lntd3MftCaIM3scTPbbGbLG3jvVDP7yMw+MbObAIJxXt6aY0VKl6w0vtqtJCMiiS2cJLPHzMbXLpjZOGBPczu5+0fuXujuhcCRwG7g+brbmFl3M8uut+7wBoqbBpxaf6WZJQMPAqcBQ4DJZjak2TNqA2ouExEJL8lcBTxoZmvM7HMCQ/9f2cLjnAB82sBEaMcDL5pZBoCZTQHur7+zu78FbG2g3DHAJ8GaSwUwA/hmC2OLitzgxGXuDc1gLSKSGJpNMu6+xN1HAiOA4e5e5O5LW3icC4DpDZT9DPAyMMPMLgQuI9DvE658YF2d5RIg38zyzOwhoMjMftrQjmZ2upk9XFZW1oLDhS83K42qGmd7eVVUyhcRaQ8avYXZzC5y9z+Z2Y311gPg7veGc4DgNAFnAA1+2bv7XWY2A/g9gf6fnWHGDmANF+mlNFPbcvdZwKzRo0dPacHxwlb3gcyczNRoHEJE5KDXVE2mds6Y7AZeHVtwjNOAd919U0NvmtmxwDAC/TW3tqBcCNRc+tZZ7gNsaGEZUfGfJLO3mS1FROJXozUZd/9D8J+vufv8uu8FO//DNZkGmsqC5RQBjwATgdXAn8zsl+5+c5hlLwQGmFkBsJ5As9y3WxBb1NSOX1aq25hFJIGF0/H/QJjr9mNmHYCTCIx71pAOwLnu/qm71wDfBerfHICZTQf+BQwysxIzuxzA3asITKb2CrASeNrdV4QTW7R1yQo0kek2ZhFJZE31yRwNHAN0q9cv0wlIDqdwd98N5DXx/vx6y5UEajb1t5vcRBlzgDnhxNOWQjUZ3cYsIgmsqbHL0gj0vaSw76Qo24FzohlUPMhMSyYzNVlP/YtIQmuqT+ZN4E0zm9bA8y0ShtysNA2SKSIJLZw+mUfNrHPtgpl1MbNXohhT3MjNSmOr+mREJIGFk2S6uvu22gV3/wroHr2Q4odqMiKS6MJJMjVmdkjtgpkdSmDSMmlGXlaabmEWkYQWzsyYPwPeNrM3g8vHAVOjF1L8yNVIzCKS4MKZGfNlMxsFHEVgGJcb3P3LqEcWB3I7prG7opo9FdVkpoV117eISFwJp7kMIJ3AKMhlwBAzOy56IcWPnp0yAPhie3mMIxERiY1mazJm9n+B84EVQE1wtQNvRTGuuNC7cyYA67/aQ0HXrGa2FhGJP+H0yZwJDHJ3jfTYQvnBJLNhW7NzvImIxKVwmss+AzRWfSv06JSBGaxXkhGRBBVOTWY3sMTMXgdCtRl3vy5qUcWJtJQkumenqyYjIgkrnCTzUvAlrdC7cyYbypRkRCQxhXML85NtEUi8yu+cyYoN22MdhohITIRzd9lqGnjC3937RyWiOJPfOZNXP9iEu4emrhYRSRThNJeNrvPvDOBcIDc64cSf3p0zqaiqoXRXBV07psc6HBGRNtXs3WXuXlrntd7dfwN8vQ1iiwu9dRuziCSwcJrLRtVZTCJQs8luZHOpp3fnwFP/67/aw4g+nZvZWkQkvoTTXHZPnX9XAauB86ITTvypfSBTz8qISCJqNMmY2fXufh9wi7u/3YYxxZWczFQ6pCWzYZvGLxORxNNUn8ylwZ/3t0Ug8crMyO+cqT4ZEUlITTWXrTSzNUA3M3u/znoD3N1HRDWyOKIHMkUkUTWaZNx9spn1BF4Bzmi7kOJP786ZrNhQFuswRETaXJMd/+7+BTCyjWKJW/mdM/hyZwXlldVkpGryMhFJHOFOWiYHoPZZmY1l6vwXkcSiJNMG6k5eJiKSSMJOMmamqR1bSZOXiUiiajbJmNkxZvYBsDK4PNLMfhf1yOKIJi8TkUQVTk3m18ApQCmAuy8FjotmUPEmLSWJHtkZqsmISMIJq7nM3dfVW1UdhVjiWu/OGXpWRkQSTjhJZp2ZHQO4maWZ2Y8INp1J+Hp3ztTQMiKScMJJMlcC3wfygRKgMLgsLZDfOZP12/bgvt/8byIicSuc6Ze/BC5sg1jiWu3kZZt37KVHp4xYhyMi0ibCmU+mGzAF6Fd3e3e/LHphxZ/ifoHJROcs28il4wpiHI2ISNsIp7nsRSAHeA34a52XtMCQ3p0Y0SeHGQvWqclMRBJGOJOWdXD3n0Q9kgRwQfEh/Pfzy3hv3TZGHdIl1uGIiERdODWZ2Wb2jahHkgDOKOxNh7RkZixYG+tQRETaRKNJxsx2mNl24HoCiWaPmW2vs15aqGN6CqeP6M2spRvZUV4Z63BERKKu0STj7tnu3in4M8ndM+ssd2rLIOPJBWP6sqeympeWboh1KCIiURfO2GWvh7NOwlPYtzODe2YzY0H9QRREROJPU81lGWaWB3Q1sy5mlht89QN6t1WA8cbMuKC4L8vWlzFn2cZYhyMiElVN3V32PeAHBBLKYsCC67cDD0Y5rrh2wZhDeHHpBn4wYwmdO6RyzGFdYx2SiEhUNNUnc5+7FwA/cvf+7l4QfI1099+2YYxxJyM1mScuKebQvA5M/eNilq8vi3VIIiJR0WyfjLs/0BaBJJrOHdL44+Vj6JSRwiVPLGD1l7tiHZKISMRp+uUY6pWTyR8vH0uNw7cfeYd1W3fHOiQRkYhqquN/XPBnetuFk3gO796RP10+lt0V1Ux+5B3NnikicaWpmsz9wZ//aotAEtmQ3p340+VjKdtTybcfeYfN2zXvjIjEh6aSTKWZPQHkm9n99V9tFWCiGN4nhycvG8OWHXuZ8r+LKa/U5KMi0v41lWQmAa8A5QRuYa7/kggbdUgX7rugiPdLtvGjZ5ZqtGYRafcafU4mOFnZDDNb6e5L2zCmhHbSkB785NTB/OpvHzKgezbXnzgg1iGJiLRaOHeXlZrZ82a22cw2mdlMM+sT9cgS2PeO68/Zo/rw69dW8eaqLbEOR0Sk1cJJMk8ALxF48j8fmBVcJ1FiZvzPWcPo0Smd//3X57EOR0Sk1cJJMt3d/Ql3rwq+pgHdohxXwktPSeabhfnM+2gzX+2qiHU4IiKtEk6S2WJmF5lZcvB1EVAa7cAEvlnYm6oa568aSFNE2qlwksxlwHnAF8BG4JzgOomyIb06MbBHR154b32sQxERaZVwxi5b6+5nuHs3d+/u7me6+0HdUWBm/c3sMTN7NtaxHAgz48yifBZ9/hVrSzXkjIi0P1Edu8zMOpvZs2b2oZmtNLOjW1nO48G725Y38N6pZvaRmX1iZjcBuPtn7n75gcZ/MDhjZGDqnheXqDYjIu1PtAfIvA942d0HAyOBlXXfNLPuZpZdb93hDZQzDTi1/kozSyYwt81pwBBgspkNiUzoB4c+XTowpiCXF5as18OZItLuRC3JmFkn4DjgMQB3r3D3bfU2Ox540cwygvtM4T9jpoW4+1vA1gYOMwb4JFhzqQBmAN8MM77TzezhsrKDfy6XbxXl8+mWXSxfvz3WoYiItEizScbMegT7N/4WXB5iZuE0RfUHtgBPmNl7ZvaomWXV3cDdnwFeJjCywIX85yaDcOUD6+oslxAYay3PzB4Ciszspw3t6O6z3H1qTk5OCw4XG98Y1ovUZFOTmYi0O+HUZKYRGMOsd3B5FYFpmZuTAowCfu/uRcAu4Kb6G7n7XQTGR/s9cIa77wyj7FrWwDp391J3v9LdD3P3/9OC8g5KOR1SOW5AN+Ys26gmMxFpV8JJMl3d/WmgBsDdq4BwhgguAUrc/d/B5WcJJJ19mNmxwDDgeeDWcIKud4y+dZb7ABtaWEa78I3hvdhQVs576+q3OIqIHLzCSTK7zCwPcAAzOwpotiPD3b8A1pnZoOCqE4AP6m5jZkXAIwT6US4Fcs3sl+GHz0JggJkVmFkacAGBIXDizolDepCabMx5Xw9mikj7EU6SuZHAF/dhZjYf+CNwbZjlXws8ZWbvA4XA/9R7vwNwrrt/6u41wHeB/Z7BMbPpBCZPG2RmJbV9QsFa1TUEmvNWAk+7+4owY2tXcjJTOXZAN/62/As1mYlIu9HoUP+13P1dMzseGESgD+Qjd68Mp3B3XwKMbuL9+fWWKwnUbOpvN7mJMuYAc8KJp72bOLwXb3y4mSXrtlF0SJdYhyMi0qxmk4yZXVxv1Sgzw93/GKWYpBG1TWZ/fX+jkoyItAvhNJcV13kdC9wGnBHFmKQRajITkfYmnOayffpfzCwH+N+oRSRN+oaazESkHWnNE/+7Ac0JHCMnBZvMnvznmliHIiLSrHD6ZGYRvH2ZQFIaAjwdzaCkcTmZqVx5/GE88MYnjDu8K+eO7tv8TiIiMdJskgHurvPvKuBzdy+JUjwShh+cOJBFa77ilheXM7xPDoN7dop1SCIiDQpnPpk367zmK8HEXnKScd/kQrIzUrn6T++yc29VrEMSEWlQo0nGzHaY2fa7xOvWAAAQ2klEQVQGXjvMTMMBx1j37AwemFzEmtJdTHlyEVt3VcQ6JBGR/TSaZNw92907NfDKdne1zxwEjuqfxz3njWTx2q84/YG3WVZy8E9bICKJJey7y4ITjB1S+4pmUBK+bxX1YeaVxwBw9kP/ZNr81VRV18Q4KhGRgHDmkznDzD4GVgNvAmuAv0U5LmmB4X1ymHXteI7qn8dtsz5g4v1v84+Pt8Q6LBGRsGoyvwCOAla5ewGB0ZTnN72LtLXcrDSevLSYhy4axe7KKr7z2AIum7aQD79Q95mIxE44SabS3UuBJDNLcve5BEZUloOMmXHqsF68duPx/PS0wSxcs5XT7vsHP3x6KSVf7Y51eCKSgMJ5TmabmXUE3iIwbP9mAs/LyEEqPSWZ7x1/GOcX9+V38z5l2j/X8NLS9ZxZmM+VEw7jsG4dYx2iiCQIa26gRTPLIjA9sgEXAjnAU8HaTbs3evRoX7RoUazDiKoN2/bwhzc/ZcbCdVRU13DykB5cPr4/xf26YNbQDNYiIk0zs8Xu3uhULqHtGksyZvZb4M/u/s9IB3cwSYQkU+vLnXt5Yv5q/vTOWsr2VDK0dycuOaYfE0f0okNaOJVaEZGASCSZ6wlMZ9wL+AswPTgJWVxJpCRTa09FNc+/t54n5q/m48076ZiewsThvThndB+OPKQLSUmq3YhI0w44ydQp6FACyeYCIAOYDsxw91WRCDTWEjHJ1HJ3Fq75imcWreOvyzayu6Kanp0yOG14TyYO70XRIV1IVsIRkQZELMnUK7QIeBwY4e7JBxDfQSORk0xdu/ZW8fcPNvHXZRt5c9UWKqpqgpOkdeW4gd04un8efbpkqg9HRIDwk0w4Q/2nAqcSqMmcQOCBzNsPOEI5qGSlp3BmUT5nFuWzo7ySeR9t4a1VW3hz1RZmv78RgJ6dMhjdrwsj+3RmaH4nhvbOISczNcaRi8jBrKk+mZOAycBEYAEwA3jB3Xe1XXjRp5pM09ydVZt2smB1KQvWfMWiNVvZWFYeer9Hp3T6d+1I/25ZHJLbgd6dM+ndOYPu2RnkZqXRIS35gGo/NTXOtj2VlO7cy5ade9myI/jauZetOyvYuquCsj2VoQmPks3o3imdPl060Dc3k68P7k6vnMwD/BREpL5IdPzPBf4MzHT3rRGO76ChJNNypTv3smLDdlZs2M4nm3fy2Zc7+XTzTraX7//4VEZqEp0yUumYnkLHjBQyUpNJS04iNdn26e+pqnEqqmrYW1XDnopqduytZEd5FTvKq6iu2f93NC05idysNHKz0sjJTA2VVVldwxfby9mwbQ+V1U6SwXEDu3H+6L6cNKQHKcmtmQxWROqLSp9MPFKSiZzt5ZVs3FbOhrI9fLljL6W7AjWN7Xsq2bm3ip17qyivrKay2qmsrqGq2qmt5CQnGekpSaSlJJGRkkx2RgrZGankZKaS1zGNvI7pdM1Ko3undLp1zKBTZkqTNaSaGmd16S6ef3c9zy4u4Yvt5YwpyOW3k4vo3imjjT4RkfilJBMmJZn4V13jzHy3hFtfXEFWegoPTC7i6MPyYh2WSLsWbpJR24HEveQk47zRfXnxmnF0ykzhwkff4cUl62MdlkhCUJKRhDGwRzYvXTOeIw/tws3PL2dj2Z5YhyQS95RkJKF0TE/h7nNHUlXj/Pdzy0j05mKRaFOSkYRzaF4WPz5lEHM/2sJz76rZTCSalGQkIV1yTD9GH9qF22etYPP28uZ3EJFWUZKRhJSUZNx1zgj2VtVwz6txMQyfyEFJSUYSVv9uHTlrVB9eXLqest2VsQ5HJC4pyUhCu+ioQyivrGHmuyWxDkUkLinJSEIb2juHokM689S/P9edZiJRoCQjCe+isYfy6ZZd/OuzuJhRXOSgoiQjCW/iiF507pDKU++sjXUoInFHSUYSXkZqMuce2YdXVnyh25lFIkxJRgT49thDqapx/rJwXaxDEYkrSjIiQEHXLMYW5PLS0g2xDkUkrijJiASdNqwnH2/eyadbdsY6FJG4oSQjEnTy0J4AvLLiixhHIhI/lGREgnp3zmRknxxeWbEp1qGIxA0lGZE6ThnWk6XrtmmuGZEIUZIRqeOUYJPZq6rNiESEkoxIHYd168iA7h15ebn6ZUQiQUlGpJ5ThvZkwZqtbN1VEetQRNo9JRmRek4d1pPqGue1lWoyEzlQSjIi9Qzt3Yn8zplqMhOJACUZkXrMjEkje/HWqi18uXNvrMMRadeUZEQacFZRH6pqnFkaZkbkgCjJiDRgUM9shvbuxPPvrY91KCLtWkqsA4i1HTt2cPvtt4eWp0yZAsAjjzwSWnf88cczYcIE7rnnHnbuDIxr1atXL6ZOncqsWbN49913Q9veeOONbNiwgRkzZoTWTZo0iSOPPHKf4wwcOJDJkyczffp0Vq1aFVp/6623snjxYmbPnh1ad8EFF9C7d2/uvffe0LpRo0Zx+umn8/DDD7Nx40YAOnbsyA9/+EPmzZvHm2++qXM6wHMaA6ytzuGTzSNZ+PrsuDineLxOOqfYnFO4LNGnnB09erQvWrQo1mHIQWjzjnKO/j9v8L3j+vNfpw6OdTgiBxUzW+zuo5vbTs1lIo3onp3BsQO68sJ766mpSew/xkRaS0lGpAnfKspnQ1k576wujXUoIu2SkoxIE04e0pOO6SnMXKwbAERaQ0lGpAmZacmcNSqfF5asZ/n6sliHI9LuKMmINOPGkwbSpUMaP31uGVXVNbEOR6RdUZIRaUbnDmncfsZQlq0v44n5a2Idjki7oiQjEoZvDO/JiUd0556/f8Ta0t2xDkek3VCSEQmDmfGLM4eRkpTE9X95TzNnioRJSUYkTL1yMvnV2cNZuXE7J937Fn9653M9PyPSjIQfVkakJSaN6M2I/M7c9Nz73PzCcv7877WMH9CVwr6dGdKrE7kd08hOT8HMYh2qyEFBSUakhQ7J68BTV4zlmcUl/Pnfa5k2fw0Vde46S04yOqankJqcRGqykZJsJFngZQDB/GPQaDKKVIpqqp6lNCi/u3AUA3pkR/UYSjIirWBmnDe6L+eN7sveqmpWbtzBqk07KNtdSdmeSnaUV1JZ41RW1VBV47g7NQ41wbECPfSf/XmTqaEVsTaQTiJ9DGmf0lOSo34MJRmRA5Sekkxh384U9u0c61BEDjrq+BcRkaiJyyRjZv3N7DEzezbWsYiIJLKoJhkzW2Nmy8xsiZm1etIWM3vczDab2fIG3jvVzD4ys0/M7CYAd//M3S8/kNhFROTAtUVN5mvuXtjQ5DZm1t3MsuutO7yBMqYBpzawfzLwIHAaMASYbGZDIhK1iIgcsFg3lx0PvGhmGQBmNgW4v/5G7v4WsLWB/ccAnwRrLhXADOCb4RzYzE43s4fLyjSyrohItEQ7yTjwqpktNrOp+73p/gzwMjDDzC4ELgPOa0H5+cC6OsslQL6Z5ZnZQ0CRmf20wcDcZ7n71JycnBYcTkREWiLatzCPc/cNZtYd+LuZfRislYS4+11mNgP4PXCYu+9sQfkNPU/m7l4KXNn6sEVEJBKiWpNx9w3Bn5uB5wk0b+3DzI4FhgXfv7WFhygB+tZZ7gNsaFWwIiISceYenSd/zSwLSHL3HcF//x24w91frrNNETAdmAisBv4EfObuNzdQXj9gtrsPq7MuBVgFnACsBxYC33b3FS2IcwvweYtPMKAr8GUr923PEvG8E/GcITHPOxHPGVp+3oe6e7fmNopmc1kP4Png2EwpwJ/rJpigDsC57v4pgJl9F7ikfkFmNh2YAHQ1sxLgVnd/zN2rzOwa4BUgGXi8JQkGIJwPqTFmtqihu+biXSKedyKeMyTmeSfiOUP0zjtqScbdPwNGNrPN/HrLlcAjDWw3uYky5gBzWhmmiIhEUaxvYRYRkTimJHNgHo51ADGSiOediOcMiXneiXjOEKXzjlrHv4iIiGoyIiISNUoyrdTQwJzxxsz6mtlcM1tpZivM7Prg+lwz+7uZfRz82SXWsUaamSWb2XtmNju4XGBm/w6e81/MLC3WMUaamXU2s2fN7MPgNT863q+1md0Q/N1ebmbTzSwjHq91Q4MMN3ZtLeD+4Hfb+2Y26kCOrSTTCgk0MGcV8EN3PwI4Cvh+8DxvAl539wHA68HleHM9sLLO8v8Ffh0856+AeBzl+z7gZXcfTODO0JXE8bU2s3zgOmB08Pm7ZOAC4vNaT2P/QYYbu7anAQOCr6kERmNpNSWZ1mn1wJztibtvdPd3g//eQeBLJ5/AuT4Z3OxJ4MzYRBgdZtaHwAPCjwaXDfg6UDs/UTyecyfgOOAxAHevcPdtxPm1JvAYR2bwwe4OwEbi8Fo3MshwY9f2m8AfPeAdoLOZ9WrtsZVkWqfBgTljFEubCI64UAT8G+jh7hshkIiA7rGLLCp+A/wXUBNczgO2uXtVcDker3d/YAvwRLCZ8NHgSB1xe63dfT1wN7CWQHIpAxYT/9e6VmPXNqLfb0oyrdPgwJxtHkUbMbOOwEzgB+6+PdbxRJOZTQI2u/viuqsb2DTerncKMAr4vbsXAbuIo6axhgT7IL4JFAC9gSwCTUX1xdu1bk5Ef9+VZFonYQbmNLNUAgnmKXd/Lrh6U231Ofhzc6zii4JxwBlmtoZAM+jXCdRsOgebVCA+r3cJUOLu/w4uP0sg6cTztT4RWO3uW4KjjTwHHEP8X+tajV3biH6/Kcm0zkJgQPAulDQCnYUvxTimiAv2RTwGrHT3e+u89RLw3eC/vwu82NaxRYu7/9Td+7h7PwLX9Q13vxCYC5wT3CyuzhnA3b8A1pnZoOCqE4APiONrTaCZ7Cgz6xD8Xa8957i+1nU0dm1fAi4O3mV2FFBW26zWGnoYs5XM7BsE/sKtHZjzzhiHFHFmNh74B7CM//RP/DeBfpmngUMI/I96rrs3NHNpu2ZmE4AfufskM+tPoGaTC7wHXOTue2MZX6SZWSGBmx3SgM+ASwn8IRq319rMbgfOJ3An5XvAFQT6H+LqWtcdZBjYRGBalRdo4NoGE+5vCdyNthu41N0XtfrYSjIiIhItai4TEZGoUZIREZGoUZIREZGoUZIREZGoUZIREZGoUZIRiRAz2xn82c/Mvh3hsv+73vI/I1m+SLQoyYhEXj+gRUkmOLJ3U/ZJMu5+TAtjEokJJRmRyPsVcKyZLQnOV5JsZv/PzBYG5+f4HgQe9gzO1/NnAg+8YmYvmNni4BwnU4PrfkVgpOAlZvZUcF1trcmCZS83s2Vmdn6dsufVmR/mqeBDdiJtKqX5TUSkhW4iOFIAQDBZlLl7sZmlA/PN7NXgtmOAYe6+Orh8WfCp60xgoZnNdPebzOwady9s4FhnAYUE5n/pGtznreB7RcBQAuNOzScwLtvbkT9dkcapJiMSfScTGAtqCYEhefIITAgFsKBOggG4zsyWAu8QGKRwAE0bD0x392p33wS8CRTXKbvE3WuAJQSa8UTalGoyItFnwLXu/so+KwNjo+2qt3wicLS77zazeUBGGGU3pu54W9Xo/3eJAdVkRCJvB5BdZ/kV4KrgtAmY2cDghGD15QBfBRPMYAJTXteqrN2/nreA84P9Pt0IzG65ICJnIRIB+stGJPLeB6qCzV7TgPsINFW9G+x830LDU/q+DFxpZu8DHxFoMqv1MPC+mb0bnHqg1vPA0cBSAhNL/Ze7fxFMUiIxp1GYRUQkatRcJiIiUaMkIyIiUaMkIyIiUaMkIyIiUaMkIyIiUaMkIyIiUaMkIyIiUaMkIyIiUfP/AcVa3qdMXJ8cAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.semilogy(log, label='Riemannian gradient descent')\n",
    "plt.axhline(y=val, lw=1, ls='--', color='gray', label='TT-round(A)')\n",
    "plt.xlabel('Iteration')\n",
    "plt.ylabel('Value of the functional')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Exponential machines ([Novikov et al 2017](https://arxiv.org/pdf/1605.03795.pdf))\n",
    "\n",
    "- Linear model depends on vector \n",
    "\n",
    "$$ \\hat{y}_l(x) = \\langle w, x \\rangle + b$$ \n",
    "\n",
    "- Exponential machine depends on tensor $W$\n",
    "\n",
    "$$\\hat{y}_{expm}(x) = \\sum_{i_1 = 0}^1 \\ldots \\sum_{i_d = 0}^1 W_{i_1, \\ldots, i_d} \\prod_{k=1}^d x_k^{i_k}$$\n",
    "\n",
    "- Represent tensor $W$ in TT-format and control TT-rank \n",
    "- Learning of this model is loss minimization with constraint on the tensor rank - riemannien optimization again"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Quantized Tensor Train\n",
    "\n",
    "Consider a 1D array $a_k = f(x_k)$, $k=1,\\dots,2^d$ where $f$ is some 1D function calculated on grid points $x_k$.\n",
    "\n",
    "Let $$k = {2^{d-1} i_1 + 2^{d-2} i_2 + \\dots + 2^0 i_{d}}\\quad i_1,\\dots,i_d = 0,1 $$ \n",
    "be binary representation of $k$, then\n",
    "\n",
    "$$\n",
    "    a_k = a_{2^{d-1} i_1 + 2^{d-2} i_2 + \\dots + 2^0 i_{d}} \\equiv \\tilde a_{i_1,\\dots,i_d},\n",
    "$$\n",
    "\n",
    "where $\\tilde a$ is nothing, but a reshaped tensor $a$. TT decomposition of $\\tilde a$ is called **Quantized Tensor Train (QTT)** decomposition. \n",
    "\n",
    "Interesting fact is that the QTT decomposition has relation to wavelets, more details see [here](https://epubs.siam.org/doi/abs/10.1137/100811647).\n",
    "\n",
    "Contains <font color='red'>$\\mathcal{O}(\\log n r^2)$</font> elements!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Cross approximation method\n",
    "\n",
    "- If decomposition of a tensor is given, then there is no problem to do basic operations fast. \n",
    "\n",
    "- However, the question is if it is possible to find decomposition taking into account that typically tensors even can not be stored. \n",
    "\n",
    "- **[Cross approximation](https://ac.els-cdn.com/S0024379509003747/1-s2.0-S0024379509003747-main.pdf?_tid=0739e049-650f-4fed-8087-61563927f6ac&acdnat=1544040611_b7cd24f6dd7a48d85665d6b562129ab4)** method allows to find the decomposition using only few of its elements (more details in the blackboard)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Tensor networks\n",
    "\n",
    "- Universal way for writing tensor decompositions\n",
    "- Can be interpreted as special neural network architecture\n",
    "- Some theoretical results on DNN expressive power can be derived with tensor network framework, see an example [here](https://arxiv.org/pdf/1711.00811.pdf)\n",
    "- [tntorch](https://github.com/rballester/tntorch) combines tensor networks and PyTorch\n",
    "- More details see [here](https://arxiv.org/pdf/1609.00893.pdf) and [here](https://arxiv.org/abs/1708.09165)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Summary \n",
    "\n",
    "* Tensor decompositions - useful tool to work with multidimensional data\n",
    "* Canonical, Tucker, TT, QTT decompositions\n",
    "* Cross approximation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Next week\n",
    "\n",
    "- Exam: Monday and Tuesday"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<link href='http://fonts.googleapis.com/css?family=Fenix' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=Alegreya+Sans:100,300,400,500,700,800,900,100italic,300italic,400italic,500italic,700italic,800italic,900italic' rel='stylesheet' type='text/css'>\n",
       "<link href='http://fonts.googleapis.com/css?family=Source+Code+Pro:300,400' rel='stylesheet' type='text/css'>\n",
       "<style>\n",
       "    @font-face {\n",
       "        font-family: \"Computer Modern\";\n",
       "        src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');\n",
       "    }\n",
       "    div.cell{\n",
       "        /*width:80%;*/\n",
       "        /*margin-left:auto !important;\n",
       "        margin-right:auto;*/\n",
       "    }\n",
       "    h1 {\n",
       "        font-family: 'Alegreya Sans', sans-serif;\n",
       "    }\n",
       "    h2 {\n",
       "        font-family: 'Fenix', serif;\n",
       "    }\n",
       "    h3{\n",
       "\t\tfont-family: 'Fenix', serif;\n",
       "        margin-top:12px;\n",
       "        margin-bottom: 3px;\n",
       "       }\n",
       "\th4{\n",
       "\t\tfont-family: 'Fenix', serif;\n",
       "       }\n",
       "    h5 {\n",
       "        font-family: 'Alegreya Sans', sans-serif;\n",
       "    }\t   \n",
       "    div.text_cell_render{\n",
       "        font-family: 'Alegreya Sans',Computer Modern, \"Helvetica Neue\", Arial, Helvetica, Geneva, sans-serif;\n",
       "        line-height: 1.2;\n",
       "        font-size: 120%;\n",
       "        /*width:70%;*/\n",
       "        /*margin-left:auto;*/\n",
       "        margin-right:auto;\n",
       "    }\n",
       "    .CodeMirror{\n",
       "            font-family: \"Source Code Pro\";\n",
       "\t\t\tfont-size: 90%;\n",
       "    }\n",
       "/*    .prompt{\n",
       "        display: None;\n",
       "    }*/\n",
       "    .text_cell_render h1 {\n",
       "        font-weight: 200;\n",
       "        font-size: 50pt;\n",
       "\t\tline-height: 110%;\n",
       "        color:#CD2305;\n",
       "        margin-bottom: 0.5em;\n",
       "        margin-top: 0.5em;\n",
       "        display: block;\n",
       "    }\t\n",
       "    .text_cell_render h5 {\n",
       "        font-weight: 300;\n",
       "        font-size: 16pt;\n",
       "        color: #CD2305;\n",
       "        font-style: italic;\n",
       "        margin-bottom: .5em;\n",
       "        margin-top: 0.5em;\n",
       "        display: block;\n",
       "    }\n",
       "    \n",
       "    li {\n",
       "        line-height: 110%;\n",
       "    }\n",
       "    .warning{\n",
       "        color: rgb( 240, 20, 20 )\n",
       "        }  \n",
       "\n",
       "</style>\n",
       "\n",
       "<script>\n",
       "    MathJax.Hub.Config({\n",
       "                        TeX: {\n",
       "                           extensions: [\"AMSmath.js\"]\n",
       "                           },\n",
       "                tex2jax: {\n",
       "                    inlineMath: [ ['$','$'], [\"\\\\(\",\"\\\\)\"] ],\n",
       "                    displayMath: [ ['$$','$$'], [\"\\\\[\",\"\\\\]\"] ]\n",
       "                },\n",
       "                displayAlign: 'center', // Change this to 'center' to center equations.\n",
       "                \"HTML-CSS\": {\n",
       "                    styles: {'.MathJax_Display': {\"margin\": 4}}\n",
       "                }\n",
       "        });\n",
       "</script>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from IPython.core.display import HTML\n",
    "def css_styling():\n",
    "    styles = open(\"./styles/custom.css\", \"r\").read()\n",
    "    return HTML(styles)\n",
    "css_styling()"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "celltoolbar": "Slideshow",
  "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.7"
  },
  "nav_menu": {},
  "toc": {
   "navigate_menu": true,
   "number_sections": false,
   "sideBar": true,
   "threshold": 6,
   "toc_cell": false,
   "toc_section_display": "block",
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}