tensor_manipulation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Preliminaries: Tensor Manipulation\n",
    "\n",
    "This notebook walks through the basic mathematical tools for tensor methods, using **TensorLy**.\n",
    "\n",
    "**Included is examples of the following basic operations**\n",
    "\n",
    "1. [Creating a tensor](#Creating-a-tensor)\n",
    "2. [Setting the backend](#Setting-the-backend)\n",
    "3. [Basic tensor operations](#Basic-tensor-operations)\n",
    "    1. [Tensor unfolding](#3.1-Unfolding)\n",
    "    2. [Tensor vectorization](#3.2-Tensor-Vectorization)\n",
    "    3. [Kronecker Product](#3.3-Kronecker-Product)\n",
    "    4. [Khatri-Rao Product](#3.4-Khatri-Rao-Product)\n",
    "    5. [Hadamard Product](#3.5-Hadamard-Product)\n",
    "    6. [Outer Product](#3.6-Outer-Product)\n",
    "    7. [$n$-mode Product](#3.7-$n$-mode-Product)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import tensorly as tl"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Tensors are higher order extensions of matrices that can encode multi-dimensional structure in the data\n",
    "\n",
    "In this tutorial we will show how to manipulate tensors as ndarrays, using [TensorLy](http://tensorly.github.io) with the NumPy backend to perform tensor operations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1.Creating a tensor\n",
    "\n",
    "A tensor can be represented in multiple ways. The simplest is the slice representation through multiple matrices.\n",
    "\n",
    "Let's take for this example the tensor $\\tilde{\\mathcal{X}}$ defined by its frontal slices:\n",
    "\n",
    "$$\n",
    "   \\mathbf{X}_1 = \n",
    "   \\left[\n",
    "   \\begin{matrix}\n",
    "   0  & 2  & 4  & 6\\\\\n",
    "   8  & 10 & 12 & 14\\\\\n",
    "   16 & 18 & 20 & 22\n",
    "   \\end{matrix}\n",
    "   \\right]\n",
    "$$\n",
    "\n",
    "and \n",
    "\n",
    "$$\n",
    "   \\mathbf{X}_2 =\n",
    "   \\left[\n",
    "   \\begin{matrix}\n",
    "   1  & 3  & 5  & 7\\\\\n",
    "   9  & 11 & 13 & 15\\\\\n",
    "   17 & 19 & 21 & 23\n",
    "   \\end{matrix}\n",
    "   \\right]\n",
    "$$\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "In Python, this array can be expressed as a numpy array::"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = tl.tensor(np.arange(24).reshape((3, 4, 2)), dtype=tl.float32)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[[ 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.]]], dtype=float32)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can view the frontal slices by fixing the last axis:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.,  2.,  4.,  6.],\n",
       "       [ 8., 10., 12., 14.],\n",
       "       [16., 18., 20., 22.]], dtype=float32)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X[:, :, 0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.,  3.,  5.,  7.],\n",
       "       [ 9., 11., 13., 15.],\n",
       "       [17., 19., 21., 23.]], dtype=float32)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X[:, :, 1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2.Setting the backend\n",
    "\n",
    "In TensorLy you can dynamically set the backend to use either NumPy or MXNet to represent tensors and perform the operations:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "numpy.ndarray"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "type(X)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By default, the backend is set to NumPy, here is how to change it to PyTorch:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "tl.set_backend('pytorch')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "torch.Tensor"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X = tl.tensor(np.arange(24).reshape((3, 4, 2)))\n",
    "type(X)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's change it back to NumPy for the rest of the tutorial."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "tl.set_backend('numpy')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "numpy.ndarray"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X = tl.tensor(np.arange(24).reshape((3, 4, 2)))\n",
    "type(X)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3.Basic tensor operations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.1 Unfolding\n",
    "\n",
    "Also called **matrization**, **unfolding** a tensor is done by reading the element in a given way as to obtain a matrix instead of a tensor.\n",
    "\n",
    "It is done by stacking the **fibers** of the tensor into a matrix.\n",
    "\n",
    "### Definition\n",
    "For a tensor of size $(I_1, I_2, \\cdots, I_n)$, the k-mode unfolding of this tensor will be of size $(I_k, I_1 \\times \\cdots \\times I_{k-1} \\times I_{k+1} \\cdots \\times I_n)$.\n",
    "\n",
    "\n",
    "Given a tensor $\\tilde{\\mathcal{X}} \\in \\mathbb{R}^{I_1 \\times I_2 \\times \\cdots \\times I_N}$, the mode-n unfolding of $\\tilde{\\mathcal{X}}$ is a matrix $\\mathbf{X}_{[n]} \\in \\mathbb{R}^{I_n \\times I_M}$, with $M = \\prod_{\\substack{k=1,\\\\k \\neq n}}^N I_k$ and is defined by the mapping from element $(i_1, i_2, \\cdots, i_N)$ to $(i_n, j)$, with \n",
    "\n",
    "$$\n",
    "    j = \\sum_{\\substack{k=1,\\\\k \\neq n}}^N i_k \\times \\prod_{m=k+1}^N I_m.\n",
    "$$\n",
    "\n",
    "\n",
    "\n",
    "### Convention\n",
    "\n",
    "   Traditionally, mode-1 unfolding denotes the unfolding along the first dimension.\n",
    "   However, to be consistent with the Python indexing that always starts at zero,\n",
    "   in tensorly, unfolding also starts at zero!\n",
    "\n",
    "   Therefore ``unfold(tensor, 0)`` will unfold said tensor along its first dimension!\n",
    "   \n",
    "\n",
    "### Example\n",
    "\n",
    "For instance, using the $\\tilde{\\mathcal{X}}$ previously defined:\n",
    "$$\n",
    "   \\mathbf{X}_1 = \n",
    "   \\left[\n",
    "   \\begin{matrix}\n",
    "   0  & 2  & 4  & 6\\\\\n",
    "   8  & 10 & 12 & 14\\\\\n",
    "   16 & 18 & 20 & 22\n",
    "   \\end{matrix}\n",
    "   \\right]\n",
    "$$\n",
    "\n",
    "and \n",
    "\n",
    "$$\n",
    "   \\mathbf{X}_2 =\n",
    "   \\left[\n",
    "   \\begin{matrix}\n",
    "   1  & 3  & 5  & 7\\\\\n",
    "   9  & 11 & 13 & 15\\\\\n",
    "   17 & 19 & 21 & 23\n",
    "   \\end{matrix}\n",
    "   \\right]\n",
    "$$\n",
    "\n",
    "The 0-mode unfolding of $\\tilde{X}$:\n",
    "\n",
    "$$\n",
    "   \\tilde{\\mathbf{X}}_{[0]} =\n",
    "   \\left[ \\begin{matrix}\n",
    "      0 & 1 & 2 & 3 & 4 & 5 & 6 & 7\\\\\n",
    "      8 & 9 & 10 & 11 & 12 & 13 & 14 & 15\\\\\n",
    "      16 & 17 & 18 & 19 & 20 & 21 & 22 & 23\\\\\n",
    "   \\end{matrix} \\right]\n",
    "$$\n",
    "\n",
    "The 1-mode unfolding is given by:\n",
    "\n",
    "$$\n",
    "   \\tilde{\\mathbf{X}}_{[1]} =\n",
    "   \\left[ \\begin{matrix}\n",
    "      0 & 1 & 8 & 9 & 16 & 17\\\\\n",
    "      2 & 3 & 10 & 11 & 18 & 19\\\\\n",
    "      4 & 5 & 12 & 13 & 20 & 21\\\\\n",
    "      6 & 7 & 14 & 15 & 22 & 23\\\\\n",
    "   \\end{matrix} \\right]\n",
    "$$\n",
    "\n",
    "Finally, the 2-mode unfolding is the unfolding along the last axis:\n",
    "\n",
    "$$\n",
    "   \\tilde{\\mathbf{X}}_{[2]} = \\left[ \\begin{matrix}\n",
    "      0 & 2 & 4 & 6 & 8 & 10 & 12 & 14 & 16 & 18 & 20 & 22\\\\\n",
    "      1 & 3 & 5 & 7 & 9 & 11 & 13 & 15 & 17 & 19 & 21 & 23\\\\\n",
    "   \\end{matrix} \\right]\n",
    "$$\n",
    "\n",
    "### In TensorLy\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0,  1,  2,  3,  4,  5,  6,  7],\n",
       "       [ 8,  9, 10, 11, 12, 13, 14, 15],\n",
       "       [16, 17, 18, 19, 20, 21, 22, 23]])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tl.unfold(X, mode=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0,  1,  8,  9, 16, 17],\n",
       "       [ 2,  3, 10, 11, 18, 19],\n",
       "       [ 4,  5, 12, 13, 20, 21],\n",
       "       [ 6,  7, 14, 15, 22, 23]])"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tl.unfold(X, mode=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0,  2,  4,  6,  8, 10, 12, 14, 16, 18, 20, 22],\n",
       "       [ 1,  3,  5,  7,  9, 11, 13, 15, 17, 19, 21, 23]])"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tl.unfold(X, mode=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.1 (b) Folding\n",
    "\n",
    "Folding is the inverse operation: you can **fold** an unfolded tensor back from matrix to full tensor using the ``tensorly.fold`` function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[[ 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]]])"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "unfolding = tl.unfold(X, 1)\n",
    "original_shape = X.shape\n",
    "tl.fold(unfolding, mode=1, shape=original_shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.2 Tensor Vectorization\n",
    "\n",
    "Given a tensor, $ \\mathcal{X} \\in \\mathbb{R}^{I_1 \\times I_2 \\times \\cdots \\times I_{N}}$, we can flatten it into a vector $\\text{vec}(\\mathcal{X})$ of size $\\left(I_1 \\cdot I_2 \\cdot \\ldots \\cdot I_{N}\\right)$ defined by the mapping from element $(i_1, i_2, \\cdots, i_{N})$ of $\\mathcal{X}$ to element $j$ of $\\text{vec}(\\mathcal{X})$, with $j = 1 + \\sum_{k=1}^{N} (i_k - 1) \\times \\prod_{m=k+1}^{N} I_m$.\n",
    "\n",
    "With **TensorLy**, we can use the `tensorly.base.tensor_to_vec` function for this, using our example tensor $\\mathcal{X}$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[ 0  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23]\n",
      "shape of vectorized tensor: (24,)\n"
     ]
    }
   ],
   "source": [
    "vectorized = tl.base.tensor_to_vec(X)\n",
    "\n",
    "print(vectorized)\n",
    "print('shape of vectorized tensor: {}'.format(vectorized.shape))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.3 Kronecker Product\n",
    "\n",
    "The *Kronecker Product* between two matrices $\\mathbf{A}\\in\\mathbb{R}^{I_1\\times I_2}$ and $\\mathbf{B}\\in\\mathbb{R}^{J_1\\times J_2}$ can be written in terms of block matrices as\n",
    "\n",
    "$$\n",
    "\\mathbf{A}\\otimes\\mathbf{B} = \\begin{bmatrix} a_{11} \\mathbf{B} & \\cdots & a_{1 I_2}\\mathbf{B} \\\\ \\vdots & \\ddots & \\vdots \\\\ a_{I_1 1} \\mathbf{B} & \\cdots & a_{I_1 I_2} \\mathbf{B} \\end{bmatrix}\n",
    "\\in \\mathbb{R}^{I_1 \\cdot J_1 \\times I_2 \\cdot J_2}.\n",
    "$$\n",
    "\n",
    "As a concrete example, we take two matrices\n",
    "\n",
    "$$\n",
    "\\mathbf{A} =\n",
    "\\begin{bmatrix}\n",
    "    1 & 2 \\\\ 3 & 4\n",
    "\\end{bmatrix},\n",
    "\\mathbf{B} =\n",
    "\\begin{bmatrix}\n",
    "    5 & 6 \\\\ 7 & 8\n",
    "\\end{bmatrix}.\n",
    "$$\n",
    "\n",
    "The *Kronecker Product* is then\n",
    "\n",
    "$$\n",
    "\\mathbf{A} \\otimes \\mathbf{B} = \n",
    "\\begin{bmatrix}\n",
    "    5 & 6 & 10 & 12 \\\\ 7 & 8 & 14 & 16 \\\\ 15 & 18 & 20 & 24 \\\\ 21 & 24 & 28 & 32\n",
    "\\end{bmatrix}.\n",
    "$$\n",
    "\n",
    "In **TensorLy**, we can use either the `tl.kron` method or the `tl.tenalg.kronecker` method as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "matrix A: \n",
      " [[1 2]\n",
      " [3 4]]\n",
      "matrix B: \n",
      " [[5 6]\n",
      " [7 8]]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([[ 5,  6, 10, 12],\n",
       "       [ 7,  8, 14, 16],\n",
       "       [15, 18, 20, 24],\n",
       "       [21, 24, 28, 32]])"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A = tl.tensor(np.arange(1, 5).reshape((2, 2)))\n",
    "B = tl.tensor(np.arange(5, 9).reshape((2, 2)))\n",
    "\n",
    "print('matrix A: \\n', A)\n",
    "print('matrix B: \\n', B)\n",
    "\n",
    "tl.kron(A, B)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.4 Khatri-Rao Product\n",
    "\n",
    "The *Khatri-Rao Product* between two matrices $\\mathbf{A}\\in\\mathbb{R}^{I\\times R}$ and $\\mathbf{B}\\in\\mathbb{R}^{J\\times R}$ is defined as $\\mathbf{A}\\odot\\mathbf{B}\\in\\mathbb{R}^{I\\cdot J\\times R}$ and can written in terms of the *Kronecker Product* as\n",
    "\n",
    "$$\n",
    "\\mathbf{A} \\odot \\mathbf{B} =\n",
    "\\begin{bmatrix}\n",
    "  \\mathbf{A}_{:, 1} \\otimes \\mathbf{B}_{:, 1},&\n",
    "  \\mathbf{A}_{:, 2} \\otimes \\mathbf{B}_{:, 2},&\n",
    "  \\cdots, &\n",
    "  \\mathbf{A}_{:, R} \\otimes \\mathbf{B}_{:, R}\n",
    "\\end{bmatrix}.\n",
    "$$\n",
    "\n",
    "Taking the same definitions for the matrices $\\mathbf{A}$ and $\\mathbf{B}$ above, the *Khatri-Rao product* is\n",
    "\n",
    "$$\n",
    "\\mathbf{A} \\odot \\mathbf{B} =\n",
    "\\begin{bmatrix}\n",
    "    5 & 12 \\\\\n",
    "    7 & 16 \\\\\n",
    "    15 & 24 \\\\\n",
    "    21 & 32\n",
    "\\end{bmatrix}.\n",
    "$$\n",
    "\n",
    "In **TensorLy**, we can use the `tl.kr` method to compute this, which takes a list of matrices as its input"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A khatri-rao B:\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([[ 5, 12],\n",
       "       [ 7, 16],\n",
       "       [15, 24],\n",
       "       [21, 32]])"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "print('A khatri-rao B:')\n",
    "tl.kr([A, B])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.5 Hadamard Product \n",
    "\n",
    "The Hadamard product of $\\mathbf{A}, \\mathbf{B} \\in \\mathbb{R}^{I \\times N}$ is symbolized as $\\mathbf{A} * \\mathbf{B}$ and is equal to ${a}_{(i, j)} {b}_{(i, j)}$ for the $(i, j)$ element.\n",
    "\n",
    "$$\n",
    "\\mathbf{A} * \\mathbf{B} =\n",
    "\\begin{bmatrix}\n",
    "    5 & 12 \\\\\n",
    "    21 & 32\n",
    "\\end{bmatrix}.\n",
    "$$\n",
    "\n",
    "In **TensorLy** we can compute this straight-forwardly in a backend-agnostic manner as"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A hadamard B:\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "array([[ 5, 12],\n",
       "       [21, 32]])"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "print('A hadamard B:')\n",
    "A * B"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.6 Outer Product \n",
    "\n",
    "Given a set of $N$ vectors $\\{ \\mathbf {x}^{(n)} \\in \\mathbb{R}^{I_n} \\}_{n=1}^{N}$ their outer product is defined as\n",
    "\n",
    "$$ \\tilde{\\mathcal{X}} = \\mathbf{x}^{(1)} \\circ \\mathbf{x}^{(2)} \\circ \\cdots \\circ \\mathbf{x}^{(N)}   \\doteq  \\bigcirc_{n=1}^{N} \\mathbf{x}^{(n)} \\in \\mathbb{R}^{I_1 \\times I_2 \\times \\cdots \\times I_{N}} $$\n",
    "\n",
    "and defines a rank-one $N$th-order tensor.\n",
    "\n",
    "Concretely, given three vectors $\\mathbf{x}^{(1)},\\mathbf{x}^{(2)},\\mathbf{x}^{(3)}\\in\\mathbb{R}^2$ we can compute their outer product in **TensorLy** with the `tl.tenalg.outer` method:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[[15, 18],\n",
       "        [20, 24]],\n",
       "\n",
       "       [[30, 36],\n",
       "        [40, 48]]])"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x1 = tl.tensor(np.arange(1, 3).reshape((2)))\n",
    "x2 = tl.tensor(np.arange(3, 5).reshape((2)))\n",
    "x3 = tl.tensor(np.arange(5, 7).reshape((2)))\n",
    "\n",
    "tl.tenalg.outer([x1, x2, x3])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that we can also achieve this by first computing their kronecker or khatri-rao product and then reshaping into a tensor as follows"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[[15, 18],\n",
       "        [20, 24]],\n",
       "\n",
       "       [[30, 36],\n",
       "        [40, 48]]])"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tl.tenalg.kronecker([x1, x2, x3]).reshape(2, 2, 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3.7 $n$-mode Product\n",
    "\n",
    "Also known as **tensor contraction**. This is a natural generalization of matrix-vector and matrix-matrix product. When multiplying a tensor by a matrix or a vector, we now have to specify the **mode** $n$ along which to take the product.\n",
    "### Tensor times matrix\n",
    "\n",
    "In that case we are doing an operation analogous to a matrix multiplication on the $n$-th mode. Given a tensor $\\tilde{\\mathcal{X}}$ of size $(I_1, I_2, \\cdots, I_N)$, and a matrix $\\mathbf{M}$ of size $(D, I_n)$, the $n$-mode product of $\\tilde{\\mathcal{X}}$ by $\\mathbf{M}$ is written $\\tilde{\\mathcal{X}} \\times_n \\mathbf{M}$ and is of size $(I_k, I_1 \\times \\cdots \\times I_{n-1} \\times D \\times I_{n+1} \\cdots \\times I_n)$.\n",
    "\n",
    "### Tensor times vector\n",
    "\n",
    "In that case we are contracting over the $n$-th mode by multiplying it with a vector. Given a tensor $\\tilde{\\mathcal{X}}$ of size $(I_1, I_2, \\cdots, I_N)$, and a vector $\\mathbf{v}$ of size $(I_n)$, the $n$-mode product of $\\tilde{\\mathcal{X}}$ by $\\mathbf{v}$ is written $\\tilde{\\mathcal{X}} \\times_n \\mathbf{v}$ and is of size $(I_k, I_1 \\times \\cdots \\times I_{n-1} \\times I_{n+1} \\cdots \\times I_n)$.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAboAAAGdCAYAAABkXrYLAAAKo2lDQ1BJQ0MgUHJvZmlsZQAASImVlgdUE9kax+9MeqMFkE7ovXeQEnrovdkICSWUEAJBwK4sruCKIiICyoJKVbBRZBURC7ZFQRH7giwC6rpYwILKDvAIb887+955/znf3N/5cuc/372Te84HAPkek8dLhsUASOFm8IM8nGkRkVE03DBAI5cowABjJiudRw8I8AGIFse/a/o+gObGuwZzXv/5+3+VODs2nQUAFIBwDDudlYLwGSQ6WDx+BgAoHpJXW5vBm+NihCX5SIEI181x/AKfm+OYBb4zPyckyAXh3wHAk5lMfjwApEkkT8tkxSM+ZGS1wJjL5nARdkLYgZXAZCO8GWH9lJTUOT6GsHbMv/nE/80zRujJZMYLeWEt88K7ctJ5yczs/3M7/rdSkgWL71BFgpzA9wxCRglkz+qSUr2FzI3x819kDnt+/jwnCDxDF5mV7hK1yGymq/ciC5JC6YvM5C89y8lghCwyPzVI6M9N9vMR+scyhByb7ha8yHEcd8Yi5ySEhC9yJifMb5HTk4K9l+a4CPN8QZCw5ji+u3CNKelLtbGYS+/KSAjxXKohQlgPO9bVTZjnhgrn8zKchZ685ICl+pM9hPn0zGDhsxnIH2yRE5leAUs+AcL9AQEgAvgAU2AOrIAJ8AXIijNiszLmCnZJ5WXzOfEJGTQ6cmJiaQwuy1CfZmpsYgXA3Plb+LzvH8yfK0gav5TLXg2AI/JtoA9LufBOAE6cRI5S5lJOowMAMREALnSzBPzMhRx67oYBRORcSwJZoATUgDYwQCq0BHbACbgBL+APQkAkWA1YIAGkAD5YC9aDLSAPFIDdYB8oA5XgMKgDx8Ep0AbOgYvgKrgJ7oAB8BgMgVHwCkyCaTADQRAOokBUSBZShjQgPcgUsoYcIDfIBwqCIqFoKB7iQgJoPbQNKoCKoDKoCqqHTkJnoYvQdagPeggNQxPQO+gLjILJsCSsCGvCRrA1TIe94RB4FRwPp8E5cC68Cy6Fq+FjcCt8Eb4JD8BD8Ct4CgVQJJQ0SgVlgLJGuaD8UVGoOBQftRGVjypBVaOaUB2oHtRd1BDqNeozGoumomloA7Qd2hMdimah09Ab0TvRZeg6dCv6Mvouehg9if6OoWAUMHoYWwwDE4GJx6zF5GFKMDWYFswVzABmFDONxWKlsVpYK6wnNhKbiF2H3Yk9iG3GdmH7sCPYKRwOJ4vTw9nj/HFMXAYuD3cAdwx3AdePG8V9wpPwynhTvDs+Cs/Fb8WX4Bvwnfh+/Bh+hiBG0CDYEvwJbEI2oZBwhNBBuE0YJcwQxYlaRHtiCDGRuIVYSmwiXiE+Ib4nkUiqJBtSIIlD2kwqJZ0gXSMNkz6TJci6ZBfySrKAvItcS+4iPyS/p1AomhQnShQlg7KLUk+5RHlG+SRCFTEUYYiwRTaJlIu0ivSLvBEliGqI0kVXi+aIloieFr0t+lqMIKYp5iLGFNsoVi52VmxQbEqcKm4i7i+eIr5TvEH8uvi4BE5CU8JNgi2RK3FY4pLECBVFVaO6UFnUbdQj1CvUUUmspJYkQzJRskDyuGSv5KSUhJS5VJhUllS51HmpIWmUtKY0QzpZulD6lPR96S/LFJfRl8Uu27GsaVn/so8y8jJOMrEy+TLNMgMyX2Rpsm6ySbJ7ZNtkn8qh5XTlAuXWyh2SuyL3Wl5S3k6eJZ8vf0r+kQKsoKsQpLBO4bDCLYUpRSVFD0We4gHFS4qvlaSVnJQSlYqVOpUmlKnKDsoc5WLlC8ovaVI0Oi2ZVkq7TJtUUVDxVBGoVKn0qsyoaqmGqm5VbVZ9qkZUs1aLUytW61abVFdW91Vfr96o/kiDoGGtkaCxX6NH46Omlma45nbNNs1xLRkthlaOVqPWE22KtqN2mna19j0drI61TpLOQZ07urCuhW6CbrnubT1Yz1KPo3dQr08fo2+jz9Wv1h80IBvQDTINGg2GDaUNfQy3GrYZvjFSN4oy2mPUY/Td2MI42fiI8WMTCRMvk60mHSbvTHVNWablpvfMKGbuZpvM2s3emuuZx5ofMn9gQbXwtdhu0W3xzdLKkm/ZZDlhpW4VbVVhNWgtaR1gvdP6mg3Gxtlmk805m8+2lrYZtqds/7QzsEuya7AbX661PHb5keUj9qr2TPsq+yEHmkO0w88OQ44qjkzHasfnTmpObKcapzG6Dj2Rfoz+xtnYme/c4vzRxdZlg0uXK8rVwzXftddNwi3Urcztmbuqe7x7o/ukh4XHOo8uT4ynt+cez0GGIoPFqGdMell5bfC67E32DvYu837uo+vD9+nwhX29fPf6PvHT8OP6tfkDf4b/Xv+nAVoBaQG/BGIDAwLLA18EmQStD+oJpgavCW4Ing5xDikMeRyqHSoI7Q4TDVsZVh/2Mdw1vCh8KMIoYkPEzUi5SE5kexQuKiyqJmpqhduKfStGV1qszFt5f5XWqqxV11fLrU5efX6N6BrmmtPRmOjw6Ibor0x/ZjVzKoYRUxEzyXJh7We9Yjuxi9kTsfaxRbFjcfZxRXHj8fbxe+MnEhwTShJec1w4ZZy3iZ6JlYkfk/yTapNmk8OTm1PwKdEpZ7kS3CTu5VSl1KzUPp4eL483lGabti9tku/Nr0mH0lelt2dIIo3OLYG24AfBcKZDZnnmp7Vha09niWdxs25l62bvyB7Lcc85ug69jrWue73K+i3rhzfQN1RthDbGbOzepLYpd9PoZo/NdVuIW5K2/LrVeGvR1g/bwrd15Crmbs4d+cHjh8Y8kTx+3uB2u+2VP6J/5PzYu8Nsx4Ed3/PZ+TcKjAtKCr7uZO288ZPJT6U/ze6K29VbaFl4aDd2N3f3/T2Oe+qKxItyikb2+u5tLaYV5xd/2Ldm3/US85LK/cT9gv1DpT6l7QfUD+w+8LUsoWyg3Lm8uUKhYkfFx4Psg/2HnA41VSpWFlR++Znz84Mqj6rWas3qksPYw5mHXxwJO9Jz1PpofY1cTUHNt1pu7VBdUN3leqv6+gaFhsJGuFHQOHFs5bE7x12PtzcZNFU1SzcXnAAnBCdenow+ef+U96nu09anm85onKloobbkt0Kt2a2TbQltQ+2R7X1nvc52d9h1tPxi+EvtOZVz5eelzhd2EjtzO2cv5FyY6uJ1vb4Yf3Gke03340sRl+5dDrzce8X7yrWr7lcv9dB7Llyzv3buuu31szesb7TdtLzZesviVsuvFr+29Fr2tt62ut1+x+ZOR9/yvs5+x/6Ld13vXr3HuHdzwG+g737o/QeDKweHHrAfjD9Mfvj2Ueajmcebn2Ce5D8Ve1ryTOFZ9W86vzUPWQ6dH3YdvvU8+PnjEdbIq9/Tf/86mvuC8qJkTHmsftx0/NyE+8Sdlytejr7ivZp5nfeH+B8Vb7TfnPnT6c9bkxGTo2/5b2ff7Xwv+772g/mH7qmAqWfTKdMzH/M/yX6q+2z9uedL+JexmbVfcV9Lv+l86/ju/f3JbMrsLI/JZ863Aigk4Lg4AN7VAkCJBICK9M1EkYX+eF7QQk8/T+CfeKGHnpclAEfnehEkvLoAqHYCQBNhUWQMQCLECcBmZsL4l9LjzEwXvEhtSGtSMjv7HjHA6QDwbXB2dqZtdvZbDVLsIwC6phf68jnRdwGgngb+QX8B55kBHLASDzgAAAGdaVRYdFhNTDpjb20uYWRvYmUueG1wAAAAAAA8eDp4bXBtZXRhIHhtbG5zOng9ImFkb2JlOm5zOm1ldGEvIiB4OnhtcHRrPSJYTVAgQ29yZSA1LjQuMCI+CiAgIDxyZGY6UkRGIHhtbG5zOnJkZj0iaHR0cDovL3d3dy53My5vcmcvMTk5OS8wMi8yMi1yZGYtc3ludGF4LW5zIyI+CiAgICAgIDxyZGY6RGVzY3JpcHRpb24gcmRmOmFib3V0PSIiCiAgICAgICAgICAgIHhtbG5zOmV4aWY9Imh0dHA6Ly9ucy5hZG9iZS5jb20vZXhpZi8xLjAvIj4KICAgICAgICAgPGV4aWY6UGl4ZWxYRGltZW5zaW9uPjQ0MjwvZXhpZjpQaXhlbFhEaW1lbnNpb24+CiAgICAgICAgIDxleGlmOlBpeGVsWURpbWVuc2lvbj40MTM8L2V4aWY6UGl4ZWxZRGltZW5zaW9uPgogICAgICA8L3JkZjpEZXNjcmlwdGlvbj4KICAgPC9yZGY6UkRGPgo8L3g6eG1wbWV0YT4KRvDBhAAAQABJREFUeAHsnQd8FEX7x393l7v0CiT03qSjUlQ6FlTsvb32rjT1tXex97/lFez66qtiAQUUEEGKoFQpQqiBBAgJ6f3K/mc23CWXXMjlcru5u/yGT7jd2dlnZr6zd8/OM8/MGBQRwEACJEACJEACIUrAGKL1YrVIgARIgARIQCVARccHgQRIgARIIKQJUNGFdPOyciRAAiRAAlR0fAZIgARIgARCmgAVXUg3LytHAiRAAiRARcdngARIgARIIKQJUNGFdPOyciRAAiRAAlR0fAZIgARIgARCmgAVXUg3LytHAiRAAiRARcdngARIgARIIKQJUNGFdPOyciRAAiRAAlR0fAZIgARIgARCmgAVXUg3LytHAiRAAiRARcdngARIgARIIKQJUNGFdPOyciRAAiRAAlR0fAZIgARIgARCmgAVXUg3LytHAiRAAiRARcdngARIgARIIKQJUNGFdPOyciRAAiRAAlR0fAZIgARIgARCmgAVXUg3LytHAiRAAiRARcdngARIgARIIKQJUNGFdPOyciRAAiRAAlR0fAZIgARIgARCmgAVXUg3LytHAiRAAiRARcdngARIgARIIKQJhIV07Zpx5eb/ehgLl2T6RGDa7d3Rvm2kT/eGwk1WqwPfzDmAr7/PwPr1+cjLqkBCKzNat4nAwP7xOP/sNjhzfDIMBkMoVFfzOqxak4t/UgvRo1s0RgxroXl+zIAEahKgoqtJJETOV63JwVsv7vGpNldd3KHZKrot2wrxr1vXYtOqAkTEGjF8dBLajYtAXq4Va9bmYf2yfGzcnI8J46Si8wlvs7vps2/2YebrabjqlvZUdM2u9QOjwlR0gdEOfi/F4/f1wsPTerrkmsMMag/kt+XZOGP8SrTsaEHqhlMREx0Gh0OBza64pXWdNKODDZvycdaFfyB7X4X6o/zSk33RMincRcAuGM3/NRPbUotgNFLLucDUc9ChbRT6DY9Dpw5R9aTkZRLQhoBBEUEb0ZQaiAQ8KbpALKfeZcovsGLkmcuwbU0Rbp3WGW9M7x8wyky+iNwydQMKimz4/N0TYLEE1tD6Uy9tx7I/svHFzBPRqkXVi4Hebcj8SKAuAuzR1UWG8c2KwNsf7FGVXI/B0XjukT4Bo+ScjbD+73zkZlvV3rczLlA+5fjbqqW5KC93BEqRWA4ScCMQWK+GbkXjSSAQkB3+3Hwryiv8+yMmzYBZR8pRVm6vt5o2mwO5eRXQyvhQWmbHJ5/vU8sx6fZuqjm33kLVkUDrstaRbYOii0tsyMn1nqfsUcq2Kim1NSifhiSWeeSINpafDCTgbwJUdP4mGiLypGJ7/o0daN9rAVKS56NFylxcccsaHDhU5lbDZ17Zjj7DflUdNH6YdxCnTFiGFu3notughXj0ua21FFn6wVJcKeQkd5qPdm1/QVzcXLTt8bMwF25AxqFSN9l/rc/DWZf+gbjEuUhJ+RkpXX7G5If/RnZOhVu6m6esx4ljl0CaHzdtLcC1d65Fj8GL8Oe6XLd0dZ3IHknathLV+eSs01LqSnbMeG/LKpX1dXetw+BRvyEzqxzvfbIXA0f8JpjNQ7+TFuOdj/a4fuwrhPfnjZPW4byrVmPvjhJkHyzHRdf+iXOuWKX+ra5WvxdEW8l22Crqsv9Aqcqp7/Bf8dk3+9VyS+73PLYJnfstQGLiPLRu/TM69FmA/5u5G/Klw1OQLC+4ejXiWvyktlVCQmUZX//PLpSW2vHIs1vVcixffgQVJQ7cIOrlLNu3Px5wiZy3MFONl89HzbA7rRhX3bYGia3norVoY/l51a1rIOOrh4Zwq34fj0lAEqDpks9BLQLSvV4qi+8+PYhh4xNx1x1dIH8oP/tgPzZvKcCy+SOREGdW7zuUVYadG4rxyPStWLsmH+ed1xodO0Zi3g+H8MITO2EyGfHEv3uraWXP6bwrVyP17yLccmcn9O0djxyhtNZtysOsrw9g2p3dXWX57qcDuO66dYiKM+HmuzuhY7torNmYixlvpmHRr1n4+buT0eHoFIj0jDLs2lasTgm4d9pmVJQ60H9YHGJjvHu8t+0ogkN0Vtp1iURKy4aPMTWkrLKC+9NLsWNTEa69fS327CnB6acnI61jCRbMPoypd2xCiwQLLrug3VFmDhSJHpi1zCF6tBC9KjsMRx1hqiuozGzBQLTDslVH8NQz25C9vwK9TogVY2YWVc5/v0nHzHfSMPGi1hg+KUk1zb4vlOw9d21GoRj7e2hqTzWd879FS7Nw6RV/qizPuaQNeveIRVGxFYsWZ+Gh+7dCvhBIU2WxKE+54C1DifgUnW81yGfIGdLSS/DLd4cxZmQrZ5T6KV8Ozr1kFUoKbLj4inbo0zMWO/YU4X+fZmDxr9mY880wDBmc6LqnIdxcN/GABCQB8abE0IwILF6WpZgts5U23ecrhUVWjzWf+dleNc0tU9crVqvdlWbpymwlJulH5eW3drji7rx/g5o2udM8ZcfuQlf8j78cVMIjZit9hi1SKioqZazdmKtYRNxNk9e50jkPjuSUu/LatbdIadlhrtJj8EJlX3qxM4n6+e2PGUpk3BzlzEtWKlabQ42bcPFKVW50whzlrEtX1rrHTYCHE9FDUesw7NQlit1eKdNDMo9RDS2rQ9jmxp23XM2v7/BFijDJuuRKrrJtzrtylStOHpQLft0HLVC69FuglJba3K45T6Y+8rd6r2yfwSMXKxs25zkvqZ+bt+UraTVYHjpcqrTuMl9p1Wmucji7zJU+60iZ0qnPL0pC65+UJSuyXPHyQJZfPkOlZVXlkPxlvvszStzSOk/e+XC3WraX397pjFIKCq3KcScuUlqIdl69LscVLw/+Wp+rxsvrMp0MvnBTb+R/JCAI0HTJ9x03AvJN/K33dqNlBwueFU4ZYWFVj8jI4Uk4/ZxkzBYmqJpjKRde3Bbdu8S4ZI0c3gLxrc3IzqyA7MnJ4PTv3bm7WO2ZuBKLg6REiyuvN4U5LT/Thqce7Y0O7dxd0i8Qk7XPvigFi+dmYeWfOS4RiuhA9B4ci68/GlLrHleiOg4slsqpArK3JL4TdaTyHO1LWZ2SbruxKxLiK3vGMm7CqckwiU6oND36GuJamjHnq+EY2DfeTUTfXnGiV+zOMln0Xk8cloC8gzZs3V7kSv/51+k4sLMMN97SCaNPbumKlwdykvzYES0REW5yi2/oyf++T8fOv4txx+1dMLRar03KOXFQAu64o6t6XaarGbTgVjMPnocWgapfsdCqF2vjI4H0g8IUubkIYeLH/8Gnt+IWMXbm/Lv1no3Yv78UaXtLxdhblWlKZnXCwAS3HMPEvD2zcIN3iPEfp39B/+PicMKoBKz4JQc9By2CHFv76It9bj/sopemmiYj4421fmRlBvKH9oxxKaqp8ZfF7iu/TL6tG6KjvDNXVi9s967R6umBtDLV8ab6tWMdN6asUu6g/nFu4sOF8jAKbpKZr+FfV7dH+zberWojWcbGVvISvXs1S6noF/52WD0+d0IbX4txzPtkHnLcTgbZlp7CGWOT1ej5CzJrvXxowc1TGRgXOgQa/qsQOnVnTTwQkN6N1hIF8qd2u5gYXTNERpgwYEAcFKf2OpogKqr+N3w5/+vrj4fgzns34udZh/HJu/vVP3OUAY883hv3T+qOMtH7OyKW3IqINgml5Vlmy6TKcaf9GVU9H4N4ZRvQx11x1Cx7XeeD+ycgoY1Z9GysopeYi3MntK4rqVu8r2WVQmQfMirS/18/WZe6Qp7wnv1+7kExjpeNnbuKkXW4Aum7qxjK+xzi/SVDvOzIF53WyQ0fr6wr7+rxstN84ECZ2ntNTKzq0VZPI+Nl71aWRaZ3rkKjFbfqefM49Aj4/5sWeoyaVY3ixBt+WIQBp5yShP/NHOJ13b1dJ0T2Nmb/dzjSXyrFCmF6/G1ZFv73eQaeePAf4ZwSi9PGtEJ4hBG54ke4qMSOuNjaP4RZRyq9LuWqLs4ge0KWcN8MFC3Ej+o5wonms//sx7sf7FYdLcJM9ddImnV9KauzzFp8RtTB4Itv0zHlnk0oLbTjlHFJOGlYknDmicJ/v96PtUvzXEWRCkU68dhtCgoKtZlOIPOIFi8ydiFeTnPwFGS8vC7TOZWcp3SMIwFvCPj2y+CNZKYJSgId2kWiU+8orFyRq7rra1UJuWj0Zee3w39eGYQvvzhRzeYrMR4jf6gHDYpHab4Dy/84Uit71bS2pNK0dvLQpFrXfYmQJjzpdSh7dYt/ysZ/Z6V7JUbvssp+ti9GzW07C3H77RvQsWsktm0Zj5+/ORkvPNYXd93UBV06u4/byaXNThFjsXLMc+mKbK84OBN5WzbJe9iJld6Ui5Z6zuPXo/HDhySq5mpnHvwkAV8IUNH5Qi2E7wm3CHf+6zrj0K4yPPPq9lpOJ3JcSnjd1Ro38QaJ8KCDXDS5ZmiTHCHc3QGzuXI9zrtu6aqarZ56YRsOZ5e7Jf/ltyzM++4QOveNwnlnemdidBNQx0m3ztF49aV+6hjZpLs3QvaAPDmmyPmFbwlnGdnjkD/YepRVdIDENA0DSoqEK3+1sdGaDkF1VA170krUF4eRJ7dAu9ZV43dyekJGNfOv8/7rr+qEqAQTXn9rF/buc5/PJtPINUGrm41NwsJsE9Mf8oVp1BnqK9tN13RGTAsT3n13tzr3z3mf/PxnRyHe/c9u9fqNIh0DCTSWQJXtp7GSeH/IELj9hs5YJHpNbzy3G7v3luCaSztAjsH9LebQff7VftWstWHZGPXHtyGV3rC5AGecugIThNfkqaOS0b5dhDr5++0ZcpI0cP5ZbVVx40e2xOQHuuHVZ3bhlNOX4Y6bu4jdFCLwh9iR4YP30mAWps0Z/zdImNhqmzUbUp6aaa+6uD0qhCKbOnUTrrtyHWZO2CvG69qIiczhYg6ZXWw1U4Dvvj+Iw2nlOHVsK9VTVI+ySseedu0jsXtTCR557h+cMCABcxcewilDW2CqWMmlviDHLqVS+eA/aYgX8x9PEVvlFIjJ9dJs+cei2pPqe3WLwfTn+uCeOzdh1ITluPO2rsKjNlqd4P7r71mY/20mvvluqPBurVSanTpFqWbGx57/Bxec3RbLxVw++eLy7suD6iyafLF45eV+uOPWvzFO5CG9LOU8Ojnh/Z13hNdtlg3vvDcA3TpVOgrVKYgXSMALAlR0XkBqbkmk6/g3wmnksee3Yea7e/Hjl4dcCIaOTcSj9/fyaS3IruIH8fo7OuGrLzMw96tKrzspWHpYPvFcb5cTiOwpPftwH3QVP4ZPPrMd90/Z4sp/yNgEvPF8f+GCXjWR2HWxkQcy3xtEb+YkYRKVCxXPFT3HFT9XTWGQ4geNjMerIv+eXSunUuhRVpnHraKXvXLhEcx4ba9ayy79ooTpt716XN9/7cS46NPP9MF9Uzfj2UdT1eQW4QA04fwUnH5BKyz4PquWiDuu7yx2brDggUe34NH7/nFdl+bdSf/uKpRlldn4uss74ouP0/HT/zLVP7kzxr1Te7juqevg+is7oZWY4vDvh7fgmYe3u5J1HxCN9/8zGBNP91+P3SWcB82SAHcvaJbN7n2l5aoZG8T+a3Ibn949YiDNjI0Ncj3ILdsLRQ+hQvWsHNhX9DiqOZZUly/TbhQ9SbneZjcxntSlo35v+HJNSDnn75BYqkt6m/bqHo2UVnXXX+uyyrG2dOGt2KVjFORLg1SADQlyfp5cBUYqsF5iE9QoL6ZiSPPmpn8KVBNystiZoO9xsTBXm1vpzN8pu03rcPQRq6g0ZBsjaSKW8/gOZJahbUoE+vSKaXDdnOXgJwl4IkBF54kK40iABEiABEKGAJ1RQqYpWRESIAESIAFPBKjoPFFhHAmQAAmQQMgQoKILmaZkRUiABEiABDwRoKLzRIVxJEACJEACIUOAii5kmpIVIQESIAES8ESAis4TFcaRgI8E7HY7Bg4cqLrHS/d/vf7+/e9/+1hi3kYCoU+A0wtCv41ZQ50JPPfcc3jooYdcuQ4ePBgzZswQS5x5v5KLnFtWUlKC4uJiFBUVobCwUCzXlYEdO3Zg5cqVSE1NdVuirHXr1ti6dSsSE/0/kd5VER6QQJASoKIL0oZjsQOXQEVFBc4991z88ssvrkK+8MIL8GevKycnB7NmzcJbb72FTZs2qfm89957uOWWW1x58oAESKCSABUdnwQS0IDA/v37MWzYMBw8eFCVHhkZiUWLFuHkk0/2a242mw2ff/45pk2bhk6dOmH16tWwWCr36/NrRhRGAkFMgGN0Qdx4zqKXlpbi+++/d56G7OeuXbuwatWqoKhfhw4dVHNlWFjlcrKyja677jrInpg/g5Qv5S5duhTZ2dn4/fff/SmeskggJAhQ0QV5M8rxm6uuugotW7YM8prUX/yOHTti5syZmD17dv2JAyDF2WefjcmTJ7tKIsfX7r77bkiHFX+H/v374+uvv8aHH37oNnbn73wojwSCkQBNlzq22qcbPsbMv96CQ+5qeTQYDUa3c2d89U+Z5vWJ/8EJbSo3KHVekz+YN910E0aMGIEbb7zRGR3Sn1KxT5w4Ea+88gpOOOGEgK+rdCgZP368qycqvTClY4psNy2C5HL11VcjJSVFC/G6yvT390XXwjOzwCIgvLsYdCBQYatQRs48UTn38wnK8rTflSPFRxS7w67mPHPNe0q7p8KVMe8PV8qsZa7SZBYeUl5b+bLS8ZkYZfPhTa5458EHH3ygTJgwQbFarc6oZvG5fPlyRbjwKwUFBUFR323btikJCQlyA271Ly4uTvn777+DouxNVUgtvi9NVRfm2/QEaLrU6b1j+5FtCDdF4n+XfodTOo5EUlQSZE9NhjUZq9XP03qchfCwcPVY/pcck4Lbh9yF9q16oG1MO1e8PJBjPU8++SQee+wxOMeB3BKE8MlJJ52Ezp074+233w6KWvbq1Qtvvvmma+sZoaDVcTXZO2XwTMDf3xfPuTC2uRCgotOppeds/wH3jnwQkebKXZmd2VbYK7A2XTpYGDCy0yhntOtTmrraxLVHbHisK04efPTRR2jXrh2GDh3qFt8cToxi++pbb71Vda3Py8sLiirLcVTpNOIM69atwwMPPMDxNCeQGp/+/r7UEM/TZkaAik6HBrfarfgnczNGdx5bK7f9+ftwOHcfoiLj0T9lYK3rxRXFSIxMQpixajN4YarEJ598gosvvhgmk6nWPc0hQo5LyjHK+fPnB0V1pXJ+9dVX0adPH1d5//Of/+C7775znfOgkoC/vy/kSgJUdDo8A8XWYlzU9zJEhNXenfqP/Sthd1jRM2UA4sPjPZbmvOMudIsXYz7Yvn07xo6trTirJ1SEaUwMZFWP8uuxsLyLXcLL4XDIoSf3UFBoRXlFldON+9XGn8XExGD48OHqtApZjmAIYpwOH3/8sdjZO0otrlTUt99+O3bv3q178ZWsLEC8MPkz5OZVoLS0tkep1epAXoH3efn7++LPOlJWcBKgotOh3RIiEnB+DWXlzHZZ2hL18JTOo1xjOM5r8jMxMhETe55bPQorVqxAbGwsunXr5hbvOnE4YH3mOdgefBTWMWeqx8JG5rosD2zvzUTFiPEQv0xu8d6epO4uwr/uWIfX3t2JQSN/w/6MEtetm/8pQJfjFuLS6//0qARdCRtxIE26xx9/PP7880+Ul5c3QpK+tw4ZMgTPPvusK9MsoXBuuOEGlJWVueK0PFDEXDvrVdfB9uobqEjpA/uqyvFhV55C+VkvuAzWJ552RdV3UFZuxz2PbsLLb+/ESaf9jm9mZ7huqRBK7rQLV6LngEVIP+Dds+bv74urMDxotgSq7GHNFkHTVVwdn9tf9/hcXSXbsGED2rdvj+joaI9JbF9+BUNcLMIeeRDWhx+D8ujzUO66HQbRo1CDWE3D8ZJw5CgWCuLohGaPguqILCuz46kXt+GtFwdCEb2592em4eMv9+PRe3upd8genoxfty4fpSJtdJQ2j1nPnj3VlUcOHz4MOcfOmyA8NrF3715vktZKI1c3OeussyA/GxPuvPNO/Pbbb675gHKyt1R+0rlIKnDNgngBst3zAEyPPghj926o+Ox7OF54Fabvv3Jl6di2HcoPX8PQs46XKFfKqoNX39mFs05vjfEjWyF1VxGeey0VF05sK8zqBvX9ymYTPbqDVmzdXoj2bX1n5+v3paqkPGquBLT5BWquNBtY7335acjK26+Oz/VLHuDV3dJMJ3+ok5OTIcd9agWxzqLy1bcwf/UZRHcKysIlIolQiNUWFFYOHQJ2bYXhiivc4mvJqiNizi+HMH5UMhLizPh1WRYKs21IiK9asHhA33jhDdob73+ShnCLhzLWIbeh0dIZR64reeDAAa8UnWS3ZMkSCNf+hmalppcKTi7h1VhFJ71k5fjc+vXrsW/fPlW2XAtz1KhROPXUU30qmzc3OdZvEK68LWHs3QuOPXuAjFTgjBFutzqWLVfPjSPqWKrMaQE4quyzjpRj245CPDC5B4pLbFj7Vz46dYkUz2alWNn+H759PIadshQtkjwvTSZNm4XFNiQleL7uLGBDvy+l1soeZE0HMKe8uj6zi7PQIqqlti8ddWXOeE0IUNFpgtU7oSv3ragcn2td9/hcTUnyx1q6p7dt29bjF1GOvRjGj4b4NYZ8O8dfK2C4+VYYqvX+Ks1VVpFuTE3xXp2nZ5Timss6qGm/mJWOiFgjJp7R2u3esSNaYsWfR0SHseGKTrrd33HHHUJPO/DOO+9AzDtzk+08kT1a2QPKz893Rh3zU6Z95JFHjplGr4tytwExDxJy9RSprOWfnPQv16qU17QIjn+2wXSVeLkRwfHVLPG/DcYrLlXPnf8pi5aIwwgYh7gvTqBeF2bNiuNHAu1TYFk4V43auUeMP5/TVig2A37/44jQnaWYNqWb27PZrVMU+g2OQ+eOlWOT6o1H/5PP801T1ot652LjyrHixchU/bLbcUO+L9KhZcInY9BKTNGZdfkcNznHOpHennd+cyVmXfcrhrU/6VhJeS2ICFDRNWFjOcfnRnQa7fbDUF+R5JhUXVu+GEQvJ2zy3aoI+1ffiE8HjNdc6SZS+XWJODfAONL9bd4t0TFOpt3RXb2am1+B+XMzcfLoJHRq726SysqpwClDWxxDSt2X5I/9Z5+JHqkIcpUPMSneY2K5eLFUXsE0Rle9InLFlPvvvx9PP/20Gi17d1OmTMGXX37ZoOehusxjHYddffQ5EArL8dGXQExXmE6p1nMT44TKj8KU3m8ADMJiUCuInqjx4Wkw9q3yHD3pxCQ1mVRY//1GWCeSTDj/zDZut0qnpHZtI5BYrdfvTCDb74qLOuDqSzseU8nJ9A35vkgv5TuGT0bvVlVldeZ5rM8BwvP58TNfx+A2gb/qzrHqwWvuBKjo3HnodlZuK8e69NUiPwNGdBJvyQ0I4eHhwmGuHi82MQ6nfCFc1yO7wFT97VzGL1wmTFg9YezcqQG51k76+8ojyN5fgYueqHyjr55iuXi7v/jcttWjvD6W3pRyiSzplSinEdQV5Mr9MkgewRjkj7zct04uxCzH6eS5NF/KTy2DI3UHkLoehjvFOpzVxhsdO3aKTl4aDGeKXp/T9li9IKJcLmVZPV4c5xfasHhRNoaNSES7NhFuV7duL8Kg/glqr8/twtGTCeM8KNUaCRv6fZEML+vn/oJXQ6TH084JXXDTCbd4vMbI4CVARddEbZeWv7fB43OyqPILHB8fr5rr5Ft0XT+KitweJnUTDFddLSxRVT88rvG5666F2M+lcvqBeFM3HHV5bwiOhUsPwxxhwOhT3BeUls4qe/eXoHePmIaIc6WVJkm5eHN9QZpwpXlT8vAmSF5PPPEEpDOPL0FOC5BrSUqzsb9ChGgbqdSlops6dSpuu+02f4muU47jtyXqNeNZZ7ilcY3PjRujxivCycfQqpV86KCIjV8dG/+GcfAgNzO4U8CGTfnqS8+pU5NrPZNzfj6Eyy9s50zq+szLt2LztgIM6BOHuNiqMV5XgmoHDfm+FJYXYvPhvzEgZRCiLZ4dtqqJVg9tDhvWH1qH5OhkdIrvXPNyrXM5/tfQsb9aQhihGwEqOt1Qu2e0av8f4kfahl5tBtY5f879jsozqdg6d+6sTjGQP/J1TRh37NkrbrDBMHyImxjX+Nyoyp6S7bU3YbrgPBgG9Helkz2l+pYVk0pj565itOoQjvZt3M2W83/NxBljkz2Oz3kj21WQeg4OCacaab70VvFIdnJMTK7070uQSikpqdJU58v9nu6Rm6bKHQ4uuOACyJ3JPToYebqxEXHK1m3ibhOM1dpcinONzwllJufY2abcB/OnH8CRtg/2J54BunWB/bo7Yf57Za0Xo11irE6GEwYe9exVzwDprJIt/np3d3/pkc/I3IWH0KdnHG64fR1+/3kkWidXvZAdvd314e33ZU/ubrywfDo6iZ7ZffMnYcH1yxBlrj026BIsDg4VHsTjvz0sViYag7vn3Iz7Rz+GC467qHoSt+OFu37BjV9ehPtOewZ3D5vido0ngUmAiq6J2mXp3sVqzieLdS/r6pXVVbRBgwapW7IUirdsOQnZUzA4zXliYrUrCFOg8skX6qnx+MFCDwrT3+atMN5/ryuJXDtz+vTpqkntqaeeqrNsssyxsWFCmYmxvmr+JnLc7rufDmLGa+LHskbwVnaN2+o8lRPnpZKTHqjeBrlkWqAsmyY3Zb3wwgvRtWtXdUk33TZMjZQKRXz1q/f0heeq8v0ioGN3GFq0gGN7qnA6ET1XsfKO/fX/Q9grz8P++zIou54FxHMnZr27IZfPggxmc5XZVb4MvSLm1t14dSc3s2VaegnmLcrEG9MH4EBmGe6/bwvWix7hmePrVnTefF9kfu+K3UGeHvcclu9bhpk5r6FI9O6OpejEwup4dtlTeHzM02gT2xZf/P0Jftkx75iKTnpktm3ZDW1ja/dS3aDwJGAIVPuJCpgyhWRBMosOYfHuRZBeXa/+8RKWbK/0WkvN3obZ277Hgl0/Y2+ecPn2IkgXd6nk5P5mdQX1bb1zfygLflWnGQi3PlgfeRxIShS3iB8jMd3A/qtQtgNF7+boD578oZATsGVP8d1331W9O+uSL+Ovv7ITDqeXY+mKI2oy+fY+9eFNuH9SD0RGmNxubahst5s9nEh5a9asUVdHkT2tYAvSs/TSSy9VvS2//fZbr82v/qin8RLZW7HDaapU9u+HbbJ42Rknev+RlS7+9vfeh+nKyyE8fWAYeqJqwlQ+FQ4swoFJNWfWKMipo1qhUx+xaPn3GWJsVVH/5ATyzp2ia/Xytu0owrTbuwsdaoA0eco5l71q9Ph8+b6U28txvNjKqmVUK8za8j/0az8ELaNb1Sip+2lO6RGM63oa2sa1gzR57j78D/oJL+hjheOFo8qq2zbioj6XHCsZrwUQAe5Hp1NjrBBvmL+JXpzcscBkMMFkrFQE8o3SpggHEfHDPaLDSIzpMq7eEklXdLkX2zXXXIN///vfdaaX43S2F14BMg9DTHRTXclNI06Bbcb7UBYvVU1R5scednNIkLt4yw08//nnHzz//PP1mgV/W5GNj7/Yp77JJyaYceeNXdG5g/vbvrOADZXtvM/Tp1T0csL466+/jssuu8xTkoCNk+bba6+9FnPnzsWiRYtw4okeXPk1Lr1dzJdzvP2emGIp2qqFWEt16iT1ObA99JjosRUJh5TTESYV3dGgCI9QayfhwPTZh3U6pOwX007+7/1dOHCwTJ1jec6ENjh9bKs6rQJS9D2PbcKCRVlYt3SMeIaq3rsb832R68eOfncQnjv7LVzW33uHlFViOb5LPjkNs2/8XShMel062z4kPsUPLEMQEnjppZcUYYJTxI+mJqW/9957FTE/LWBlz5kzRxGrw2hWRk0qLoSK3rIi5vIpYnqI8v3332uVjd/lWp99QSlHK8WRm6s4Dh5UlPLyRudhs9mVIeOXKLdMXd9oWdUFvLLiRaX3K+2U3NJcRYy/KcJjs/rlOo9fXP6c0ve1TkpReVGdaXghOAlUvUKFhNpuPpWQ6yNKZ4yVK1f6vdJyV+zi4mJ1PU1/C/eHbKdpddKkSXVOJvd3uf0lT+46IZ1O5Eoo559/vr/EaitHOKY4ZnwOwy2XwCA8XG2PPQWh6Bqdp1wQfMemIowdcWzzYkMykhPFZ236Amccd77q5PXs709DTk2Q4UBBBt5fOwNFFUW1RIqfbyxPW4r+7U5wG9PbkrUZ//37M0ivTGeQJs6Za99T5Tnj+BnYBKjoArt96iyd9P6TayPKycb1zqmrU4rnC3KytlzTsaFOMp6lucf6Q7ZcJ1IqebnyfzCFxYsXqyu+yCkE0tNSiyDnHsoXAX8GsYU9sDcHhnFjYP9pLgzH9YJBLCre2PDX+jzhD6Vg+BA5buyfYBU7gRSW5QkPytGYv3MuerTo5drL8Z2//g9PzJuE/1v9eq3MCsoLsO3ABjGndYzbc//Erw/h/tk346ftVaur/Lh9Np6cNxlv/flGLTmMCEwCVHSB2S5elUqO0XXo0AEffvihV+m9SSQ3BN0j1kGUis7fwR+y5Uarjz76qOqlKLfqCZYgxzwvv/xyjBs3Tp2Lp8U0Aqng5DJicpFrfwY5x9K04GPh2LQIipi2Enb3nX4Rv0SM73Y9LhrtWvvPmUh6WL5x7kws2fMr0vLScNuQO1xlvemEW3HnmEeQW5rjinMebMnaJBYgL8QpHSun3TjjHxz9OK47eQr2F+xzRuGsnhNx72nTccuJwfWi5apAczwITosrS+0kILz3FOGirginBmdUoz6FyVKzcb/GyhZb2ShXXnml8tNPPzWqjnrfnJmZqfTu3VsZOHCgcuTIEc2ynzVrlspHjgMGarDbK8tWUWFXBo5crEx+aKOuRd2Xl6a8svJFNU/hCObK+6kljyuD3uyhlFSUuOKcB0v3/qZ8u/Ub5yk/g5AAe3RB/nYjVxH5/PPPIce+/BHk6h91TUJvrPzGys7IyFBXD5GTvoMlyHaRPTm5iovcTdzfE86dHObNm6fuayd7+VqYnJ35NObz7ff3IKXzfBw6XIbZ8w/h4P4y3H5D18aIbPC9X27+L87scTZ25KSi50sp+Hzjp+qE8a/Wf4xbh0+ptdqJ9IqWpsqxncc3OC/eEDgEOGE8cNrC55LIbWPOOeccn+8PlhvlxGr5FyxBjpfJvef++usvLFiwQJOyZ4uNVB9//HHMmDFD3apo9OjRAYsnOtqE+CQznhX71e0Qq+r89N1w9Oqmn/l5a9YWtIhsieNa9sGevN2ICo/Fnxmr8OXfn2LKyIdww/E31WL3656FOLXr6eoGyLUuMiJoCHAeXdA0FQsaTASEdQfPPPMM5Ooy//3vf9XJ4f4qf6nYE05623711VfqnEfnNkUPPvig2+7l/srPn3JKS8XqPOJfVGTTv2PL3pp0QkmM9J8zjD9ZUZb/CFDR+Y8lJZGAi4BUbnJSeL9+/dSthlwXGnAgvWnl4gCyZygdTPaLFUx27typ/tX0rJTbNsmNXPv27duAHJiUBJoHASq65tHOrKWOBJYtW4YzzzxTnYuoV7Zyex85fUGr8VW96sF8SEALAk1vP9CiVpRJAk1EQJos5Xjc2LFjdS2B3OqHSk5X5MwsiAiwRxdEjcWikgAJkAAJNJwApxc0nBnvIAESIAESCCICVHRB1FgsKgmQAAmQQMMJUNE1nBnvIAESIAESCCICVHRB1FgsKgmQAAmQQMMJUNE1nBnvIAG/EJDLgvmy80RWVpa6Ua9fCkEhJNAMCFDRNYNGZhUDj4BY3FndIf3hh8UO7w0IP/zwA1q3bo1Vq1Y14C4mJYHmTYDTC5p3+7P2TURArnjy0UcfqRPLO3bs6HUp5BZKcp7e9ddfD4vF4vV9TEgCzZkAFV1zbn3WnQRIgASaAQGaLptBI7OKgUVArlspTY8NGZ+z2Wz4888/sXfv3sCqDEtDAkFAgIouCBqJRQwdAsuXL8fTTz+tmi2vuOIK1Fyc2VNNDx48iJtvvhlyl/KJEyfi22+/9ZSMcSRAAnUQoKKrAwyjScDfBMrLy/HZZ5/h5ZdfRrt27bBp0ybIntqxgty54Mknn1S3/PnXv/6FVq1aYe7cuce6hddIgARqEKCiqwGEpySgFYHc3Fyce+656g7gs2bNwoQJE+p1KJHemaeffrqqGAsLC7FlyxYMGDBAqyJSLgmEJAE6o4Rks7JSgUzg999/x6mnnqpunnriiSd6XVS52arcjmf16tU44YQTvL6PCUmguRNgj665PwGsv64E5DY+H3zwAQYPHoxBgwYhPT3dq3E6WcilS5eqc+h69+6ta5mZGQkEOwEqumBvQZY/qAjk5OTgp59+ghxvk+NvTz31lGuVk82bN+PTTz/1OG4nFeSSJUswdOhQREVFBVWdWVgSaGoCVHRN3QLMv1kRKCoqQllZGUaOHIn33nsP55xzjmvD1AcffBDXXnst5syZU4uJXC5s3bp1GDNmjDrGVysBI0iABOokwB3G60TDCyTgfwJyFRSp4N555x2MGzdOnS7gzOWJJ55Anz59kJaW5oxyfUoPzby8PFVBuiJ5QAIk4BUBKjqvMDERCfiHgMFgwNVXX63+1ZQoHUykZ2V2drZ6Sc6xMxorjS4LFy6EVJIcn6tJjeckUD8Bmi7rZ8QUJKALATlmJxdtHj9+PFJTU5GUlIRPPvkEcsK4dGCZNm0aIiMjdSkLMyGBUCLA6QWh1JqsS1AT+Pnnn1XHlDPPPBO7d+/G2LFjcdZZZ2H79u247LLL1NVRnD28oK4oC08COhOgotMZOLMjAW8JyB6edEJJTEz09hamIwES8ECAis4DFEaRAAmQAAmEDgGO0YVOW7ImJEACJEACHghQ0XmAwigSIAESIIHQIXDM6QVyNQYldOrKmpAACZAACQQoAYMol5x+o0U4pqKTGYqpPAwkQAIkQAIkoCkBo9BxGuk51KvonDXTqgBO+fwkARIgARJofgSE4VDz4LWiM0l1y0ACJEACJEACfiRgd4ghMo2VHZ1R/NhgFEUCJEACJBB4BKjoAq9NWCISIAESIAE/EqCi8yNMiiIBEiABEgg8AlR0gdcmLBEJkAAJkIAfCVDR+REmRZEACZAACQQeASq6wGsTlogESIAESMCPBKjo/AiTokiABEiABAKPABVd4LUJS0QCJEACJOBHAlR0foRJUSRAAiRAAoFHgIou8NqEJSIBEiABEvAjASo6P8KkKBIgARIggcAjQEUXeG3CEpEACZAACfiRgNeLOvsxz4AUJffecwbncV2fMl1d15zxzjQ2m03dY8kZX/PTme5Yn8e65pTnTJOSkgKjke8vkgcDCZAACUgCzUbRFRcV4elXnsYH3808uplsNcUmSUhFV30voqOKT6aSmwFWVygyee1QJU9es9jC4UiywtzfAUOY3Gip8o6qVFVHMtvKs6q46qldsbU2kDh6l/hIONQaxm1h+PzjLzF48ODK2/k/CZAACZBA81B0c36ajanPTcG+/DSXwvHY9i6NUuNqXfE1ksnTMCUMcRFxMJxcAks7W60UtXTV0RR1xcvLx7qG3HC03tUZ2aYsFAzKgl30IBlIgARIgASqCIR0jy51Ryrun34f5q+ZB7vBXo/GqILi05FQhrFKLMJ7m2DsXwiD6Zjqyacsqt+klBnQIr0Lyg9VIL1XKsJaKTCma5tn9fx5TAIkQALBQiAkFV15eTleevNFvPzZSyhWirRVcKKlLfZwxLSKgmloKYxx5SJGQ4UjFGr44RZI2JWMtPjtiDjFAWkZZSABEiABEvBMIOQU3aJfF2Ly9ElIPbxdU30jcZoUE+LNwkw5tAJhnYs9E/ZnbL4FSTs6oMiYjwODtyIikhrOn3gpiwRIIDQJhIyi27cvDQ8++wC+X/4trAYxTqWlDhC9qmhHNKJ6WICBRTCatcxMOKpUKEhM74yKdBsye+6GubUCk6YVDM2HnbUiARJongSCXtFJ9/23Zvwfnp7xFArs+doqOPGMWBzhiE6IRNiwChgTZS9OQyUnFGpYVrzoxbXG/vhdCB9pg8Y6tXl+C1hrEiCBkCYQ1Ipu+crlmPL0JGzctwHQeOqYwWFAgjkBxuOtCOtWoql+k0+cUmBGwo52KEUxDg7ahvBoDRVqSD/irBwJkEBzJxCUii4zMxMPP/sgvlz8BSrEP02VnOhVRQkzZXSXcGBwIYzh2iocxaYgPqMD7HsUZPVIE1MUaKZs7l9S1p8ESKBxBIJK0TkcDnz42Qd4+M2HkFNxpHE19+Jus8OC6NgomIWzibGV8N7U0kwppBuzotEitR0yYvfAPMoKi8a9VC8QMAkJkAAJBD2BoFF0a9auwbRnpuCP1JXa9uBkkzrESiPGBIQNssPUSwczZZEJsaltUWErw6GB22CJoYYL+m8WK0ACJBAwBAJe0eXl5eKx5x/FR3M/RJn4p6mZUjRLpD0S0R0iYTihGMZIYbfUMCh2BdEZbWDcFYacbumwdJRmSio5DZFTNAmQQDMkELCKTq4t+b9ZX+LfL9+LQ8WHtLYawqKYER0VDdOQcphaSzOltsGQHSnMlG1xKDodYaPKYNF4JRVta0PpJEACJBC4BAJS0WVlZWHyPZPw585VsIRZ0DG+Yw1FV+UQUnV0LMhVqaqOqtLbKxyIjomCvW9F5brOWc5r1VMr1dZ9rh4v14NWxOide5xTQq1PkdZ+EKgoKUdmv50wx8n7vLy3ljBGkAAJkAAJ1EcgoBRdaWkpXpz5PP7z81tQjGKgLKmy+HUZEKviq46qKizi1N0Iqq4JHSOCjK9MZS4QOwzE2hE5QkFJuJiD5wzq9ar7nOmdlys/K+XIVJWbHghFWC1B1Y4EMlJcE9WJ2NIKZWttKI4shOX0UrGiNhVcNWQ8JAESIAFNCASMopu3cC7ueX8qDii7YYjVpK4uoYYKEyLKYmGZUIGIXkIDaRzs+yyIXdEOWa0yEHZdKfCDWFGFgQRIgARIQBcCTa7oduzagYfeuR8L9/wEmMXebRp2cmSvKjw3FhF9DKIXVy72idOWsT3fiMjVrVGWW4788bsRnWxARZaGFdS2OpROAiRAAkFJQOOf+rqZlJWV4eX3X8QbC16BNVIspWWuO60/rhiFmTIqLgoRl5UjrEV1I6M/pLvLUMRSm2HbWsKyJhI5A9IQc7pR+FJSwblT4hkJkAAJ6EOgSRTdz7/Oxz0zpyDdvguGSG0rqpopS8U+cWICdkRvMT1B42Dbb0HCyo7ISshA+OVZiIngdAGNkVM8CZAACRyTgK6KbteeXXjo7fuxYPePUMx2GDTUAZVmyhih3IT5cKQwU2rcY7TnGxC9tj1KM8uRM2YHItvIHhx7ccd8+ngx4AjMmT0bGzaKtWN1CldecSW69+ihU27NK5vi4mK88urLms89dlJNTkrGbbff7jwNqE9dFJ3cCPXV91/GawteQkWEmKOmsdIxFloQFRONiEuFmbKl2Flcw6AI8eZtrYSSi0XmcbsQM8Yglu6igtMQOUVrSOD7H7/D5398pmEOVaKVMgWDBg6ioqtC4tejwsJCvPTHU4g9UZefeSQt6tp8Fd2C3xbgnhmTsc+2A4YIv7ZjLWGqmbLkqJnyOH3MlEmrOyE75hAqLjmMGG6EWqtNGBFkBOQ7moaWFjca4oXQoKX3mVtmzfPEKHSc1vtlOslqaaFz5uHrp2aqfs/e3cKb8gH8vGM2FIseZspohPcKQ9QoHcyUBQbEre+I0oxyZI1MRXgHmil9fQB5HwmQAAloTcDviq6iogKvf/AqXvn5eZRLM6XGU8akmTIyOgqRl4h94lqJLXs0DNJMGZbaEvFrk5DebTtiLzcgnGZKDYlTNAmQAAk0noBfFd2ipYuEmXIK9lZs08FMaURESZzYdVt4U/YpbzyJeiTY0i1o8WcXHIk4jIoLhJLjRqj1EONlEiABEggMAn5RdGn79gpvygcwL/V7XcyUllzhaNLzqJlS4x6jo9CA+L+7oHRvOQ6fInb67kwzZWA8uiwFCZAACXhHoFGKzmq14s2PXseLc58VZspCzc2UBulNGRWFqIuFmTJZBzPljpZosS4ZeztuRewVYlUVmim9e6qYigRIgAQCiIDPii4tLQ1jzxyFdMP+yulisqNTfcGR6sdHK6wmqRlf3/nReyOsEbCJf47eeSj+Ue4YILOVEo+Gmudq9NHrR/NQ07vl536/U5RBMSAqPw7R/RVkTtyCWHWHAefVwP105AHy5YOBBEiABEigioDPik4u4XUwYT+iOujji2zeFI4Wd4heVVtTVek1PEqfnoOYs8ODYocBh5iPVPaHGda9doSZNZ6kqCFziiYBEiABLQj4rOi0KExAydRHnzauyqJ3WrpViNgSAePxZQjvIffFYyABEiABEqhOgIquOo0gOrZlG1C+NAzGdgrCLyhTJ4WWpbvZZYOoNiwqCZAACWhHgIpOO7aaSHZUKChfbYbjiDDjjhNOOYkyG/bjNIFNoSRAAiFBgIouWJpRmim3i53K/xZmyoHliB4hCk79Fiytx3KSAAk0IQEquiaE723WthwIM6UZxmQFEecLM6WFGs5bdkxHAiRAAlR0AfwMOKzCTLlGmCkPiqXGRtuObhhLJRfATcaikQAJBCABKroAbBRZpLJUBY4NETD0E2bK4SKC+i1AW4rFIgESCHQCVHQB1kK2PNGLk2bKBCDiPGGmDKeGC7AmYnFIgASCjAAVXYA0mGITk77XhUHZZxJmSrkTg5wqQCUXIM3DYpAACQQxASq6AGi88t2AfY3wpuxTgciLxF5A1G8B0CosAgmQQKgQoKJrwpa0FQgz5e9mGKLEnLhzymDiDuVN2BrMmgRIIFQJUNE1RcuKTlvpWiOU3WGwjKyAubUsBLtxTdEUzJMESCD0CeizInPoc/S6hmW77HBsF9jFK0bUxU4l5/XtTEgCJEACJNBAAuzRNRCYr8ltuQ7kzTLAVmSA+TgHIgY6fBXF+0iABEiABBpAgIquAbB8SaoIM2XBEgesa8MQfZ4NpiQjir7xRRLvIQESIAES8IUATZe+UPPyntJUB7JfkWNvBrSY5kBkL+L2Eh2TkQAJkIDfCLBH5zeUVYKkN2X+d2IXdKsRCTfZYRa9OAYSIAESIIGmIUBF50fuqplyuTBTrgpD1EQbovrK3hyVnB8RU5TWBMQ6BYpN60wq5ctFEhSFeyhqSdteJubolujDWBFbiAVqoKLzU8uU7VZQ+K0RlgFG1UxpMFPB+QktxehIwGAzIPpItC45OqzCA9lBpyytYNvtdhg3RaAsXaffIotO+fgAjIrOB2jVb7EXKcj7QbwFFxqQcL0D5pacD1edD4+Di4BiVlCSUqxLoRXR2zCZTLrk1RwzkWwNQ8oQPcSsS/UNf7NHpwtoXTMRL6IFK+2oWGFG5AQrogfKLyyVnK5twMxIgARIwAsC7NF5AalmkvJ9Cgq+EfPhepnQYqoDBgvfSmsy4jkJkAAJBAoBKroGtIQc1C34yQhbpoL4qxVYUtiDawA+JiUBEiCBJiFARecNdmmm/NOGiqUWhJ9qReKl7MF5g41pSIAESCAQCFDR1dMKFRnC2eRrAyydw5A0WXgxRVDJ1YOMl0mABEggoAg0StE5CgBbrj6eNqZSBaV7bHCITz2CUgLkiUnfFWnCm/IyYaZsK82UNFXqwZ55kAAJkIA/Cfis6OQcjdiSeBgy/VmcumWZ7eGwrTWgIkYnZZNRAtPZdiSfL3pwRp3yrLv69V5RhHnVtsOIsjLhs81AAiRAAiTgIuCzopNzNIpb5yO8g06TBDcBMReXI7ytPqbDwsMOxA7TZ/6JqzV8PCjfB5SvFvvbxVsRERHhoxTeRgIkQAKhScBnRReaOIKrVvZiMbl3qRGGSPEScK6Y05cV+D3P4CLM0pIACYQCASq6IGxFaaYs2aDAsScMEaPEaiyt5LilVHL6jF8GITIWmQRIoBkToKILssavyADK/jAirJcYI71QbHbHTlyQtSCLSwIkoDcBKjq9ifuYn116nQozJUSLxUyU0xyo4XxEydtIgASaGQEqukBvcGGNLPlbbLWxQ5gpR4i97VrLAlPJBXqzsXwkQAKBQ4CKLnDaolZJrIeEmXKFEabuwgP0IpopawFiBAmQAAl4QYCKzgtIeidxlAtvyt9Fr81hQNSZdpii2IPTuw1kfh+8/z6WrFinS9YOewVuu/lajBw5EpmZmbihdRucqJNz0ebex+Hbf7bqUs+mykTOL73p1skwRyXoUoRIYyneeftNNa+H7r0X5lde0SXfHJHL8R98gOtvuEGX/IIlEyq6QGopaabcosD2jwkRJztgaef0pgykQjafsvy9ZTtW7Y/XpcLWkhxM3L9fzUvuun2WUHI365IzMN0m3HhDPNhsNqzaE442x1+gS03Ltn3myidctOeDrjNtD9KF+KXaZhGU0qnoAqTZrFnCTLnMBFNn4U15kdj6R6d5+AFSfRaDBEiABDQjQEWnGVrvBDsqhJlSjMMppUDk6VaExVDDeUeOqUiABEjAOwJUdN5x8n8qYZUs3abAutmE8GEOhHeUWVDJ+R80JYYKAYfDAWnWlcFoFCsCGTh2HSptq3U9+MuqNWEP8pVcA4p+EL24IoM66btSyXlIyCgSIAEXgTvuuANhYWHq3x9//OGK5wEJ1EeAPbr6CPnxuqMcyJ9vh1JiQ4SY9B0WJ99I+VbqR8QURQIkQAK1CFDR1UKiTUThWhuKfzMg/HixImWR6aiS0yavUJS6YsUKvPbaa5pUbdKkSRg1apQmsimUBEig6QlQ0WncBhWZDuR+LTZu7WxEyiQjrPnCbLlN40xDUHx6ejq+/fZbTWp2/vnnayKXQkmABAKDABWdRu2gVAC5wkxp229A0mUGmJM5HKoRaoolARIggWMSoKI7Jh7fLhZttKFokQHRYxUknafPRrG+lZR3kQAJkEDoE6Ci82MbW7OlmRIwtTYg+S4TjOF+FN7MRckxtHnz5mlCYeDAgZrIDUqhYklVpaDShV/z8peJvI5OF3DmlZWVhbVr1zpP3T73H105RkauWrUKBQUFbtflicViwbhx42rFN9cIua1XUaRNl+rHHwncFXao6PzwCCjiOcpbKHb4TjUg8VIDLG1opvQDVjcRbdq0gfxj0JaA3S40nVBAugSr2JVD5lctrFmzBmeddVa1GM+H99xzj8cLSUlJOHLkiMdrzS1SsnUUA9YsfV5c7OVUdCH7jBVvcaBwvoJo4bSXeLdQcNRxIdvWzaFiJosJhmR9prwoQqHKeXEM2hAwmcSauT2B2CFmbTKoIdX8d+C2ZeCWrAbEQDu15h41UyYqSL7TCGOkPj8OgcaB5SEBfxLo0qUL7rvvPo8iFyxYgI0bN6rXrrzySrRr165WuqioqFpxjCABKroGPgOKsLTkL7ajfLMBCZcYEd6+gQKYnARIoE4CvXv3xosvvujxuhyTcyq6O++8EyeffLLHdIwkgZoEqOhqEjnGecl2Yab8SSy+fJKClMnCm5JmymPQ8u+lrVu34ocffvCv0KPSzj33XPTr108T2RRKAiTQ9ASo6LxoA5vwQsv9RmydI8yTLW8zwBRNbF5g82uSTZs24eGHH/arTKewjh07UtE5YfCTBEKQAH+xj9Go0kxZsMyOsnUGxF9oRERnjsMdAxcvkQAJkEBAEvBZ0UnXVfsBsdVMsT4upcZ8B3IWlCMsXh9lYzsIHH7DjogTgZQpQeJN6Ucv4uLiYshtUWRwznVyfnqKc16r+elNWpkmLi4OdCSQJBhIgAT8TcBnRSf3g4IYpjJY/F0kz/Lk3lOGcGE2jNFH0ckdvlvebIQpVp/8PNfa+1hrJlD+mwn2G9znJXkvoTJlplhTcuX9DyLui/nS9OQAAB9WSURBVO8QAcW1t4KTgkHEyeA8r7yrZrrq5zK9aLuj91XeW6WRF7dpi8VJcRh42ql49hiLNvfv3x/Tp0+vzM7P/w8aNMjPEimOBEggkAj4rOik4jGlAOEd9PHIMBYbEDfagvC2+iypVbK2PCiUnKNMbN76ZzjsOYB5kCizmDvjS5A99IXvvIuYSY/gbOT7IqJB96TCgg/69ERSdi4+3LIZX5997EnCffr0gfxjaL4E5C4TF154oQqAz0LzfQ58qbnPis6XzHiPHwmITlHpVvGysTUKyoASRI0UC1pkVPWUGpLT32ITyz13TMbYDX+JXpy2oVCI/7prL6wVy8ncJDwpj4c+pm9ta0XpehDgy44elEMzDyq6IGxXW7YB1qUWGNsoCDuvGEaLuyHR2yrlZmdjyQMPoe8HH+FMaLsenlTBy9ACnw/sgNEb/8H/QfQ+vS0o05EACZBAIwhQ0TUCnt63OiqEmfIvYabMFOtpjqlAWAupPhqu5KTDyOJPPoHxpn/jLHuWDxIaVvM9CMMHPbsjorAY0zduQKuG3c7UJEACJNAoAlR0jcKn081Cn5WJzVqNm4WZsm8Zok4WEQ3Xb2pht61fj38mT8PoZUsQo3HxS4X87zp0w+/iKbs+NRXDaabUmDjFkwAJeCJAReeJSgDF2YSTSYU0U7YUC+CeWyK2/vFNwxWK5ZN+ffQxdH/zXeFsInaF1TisQhw+7t8ZwzftwNsoFX06BhIgARJoGgL8/Wka7vXm6rAKM+WacDgOCDPlSDF/MFne0nAlJ82Uy2bNQvnNU3FmfoYPEuotqluCDLEu2vtdusJeWoYnNv2N1m5Xg+tEsVtRemSfLoW2leaJ+YpiCwwR5NSdXxJb4LCizyhmamK8LnVs6kxKcvYiL221LsVQ8va78rFHRmJ6gvoFdsVpdZDjsGO42JOPwZ0AFZ07j4A4K0sVKm2jMFMeV4bIC303U+4WXo0bpt2HU36Zh0SNaya2FsPslI74JdqCa3bvxqgQMFNOnToZV15xSGNyVeK79+ihniQnJ2NOTnbVBR41mkCkUDbfzXyw0XK8FWCxjHAlfeqZZwD5x9BkBKjomgx97YxteaIX93s4jLHCTDmxBCYft/4pLS3Fwmemo8Ozr2CiDrtorkMU3u/TBYO3pgkzZZGYIRcaoUuXLpB/DMFPQM4vHX7SScFfEdbAJwJUdD5h8+9Nik1B+ToLlDQTwkaUwyymDfhippSlWr94MQ6ffzlOO7Bbc/f9w8IQ+mHHLsivqMDDW7eg9u5g/uVEaSRAAiTgCwEqOl+o+fGe8l1C2FphpuxVjoiLrJBLj/kSHGIj2DMXWzFhxgMQC9ZoGuQiY/Pi22B2i2hcIcyU40PATKkpMAonARJoUgJUdE2EX279U7EsXKzfCVjOFmbK6IY7msiiK3YFrVba8PJChxgX0z5sFmunzOjZGb13HcBb+Qc1X0lF+xoxBxIggVAnQEWncwtLxaRkGGE7HImwk8phae+7mdKa7UDKPCuuFY6B28ziT3iEtBU9LUuU54W8ZE7Vg3IM3VozrbxvS5gBByrKMS11GzpXF8RjEiABEghgAlR0OjZOyRaxQ/nXJljETO2wiWW+mykLgHEHeuFqy0BEjhUaTgSps+QeAREm4cwid5Y4Gge564N6Vvmf89jjZ/W0YlqCXLjbLrbq+flwJpan78e+vXvxx6FDmo/9VSsuD0mABEig0QSo6BqNsH4BNjF+lv9dGAyKEbFXV6BsocknJSd7gym7k/Bi+OkY3a1r/Rk3MsXGnCP4cO8edDeb8dmY8bhhzg+AUHQMJEACJBBMBKjoNGwtsUA/Cn9XYPvLgsiJFYjqa0RFlo/eJgfCMKnwZNzdaTjMPm7F421Vs0pK8MWBdOw9cgR39eyFHolJ3t7KdCRAAiQQcASo6DRqktLtDpTMNiNsoAOJ02wwmH1TcIrY12ZURk88nTwenVolaFTaSrHSTPlL9mF8t3sXzklqiUlDh6vmS00zpXASIAES0JgAFZ2fAdvyHSj4IQxKidgo9nobzK2co2ENy0j1ptydiBcsp2Nct24Nu9mH1Jtyc/BR2l50EON7rw8+ETHhwh2UgQRIgARCgAAVnZ8aURGTywqXCTPlqnBYJpQhZpBcp9A3JYeDJtxVcDImdToJFo3NlEfEKir/O3QAqZmZuLVHT/RpIVaPZmi2BIx2I6ILonWpvyK2nZI72zNoQ0CyNeyIRHmRb9akhpbKUK5PPg0tl0xPRecLtRr3lO1SUPSdCeY+ChKmWsUOA74txusoAkamd8czyaeic1dtV6eUZspF2Vn4dt9enB4ThzuGnUQzZY12bY6nDjH5v8gkHkQdgpze4vQQ1iG7ZpeFXPasPL4Upo76KCB7hrabNzemAanoGkHPXqSgcI4JtiMK4v5lhyXFtx6c4gBa7IpXzZSnduveiBJ5d+tmYab8ZP8+pIjpAy/1G4j4CM/z7ryTxlQhRUC8oxl8XLygwRxkXuIZZNCOgEVsmhDZzbcX74aWylisj0JtaLlkeio6X6gJxVSw0gHbcmGmPLUcCceLBjb6+IUVZsrbC4ZjSseTER6mbXPklJVhVlYmNh84gBvFVjoDk7VeLMwXuLyHBEiABPxLQNtfVv+WNSCkle8TvbhvhJmymxEJUypgjPDtLUaaKU9J74ZnWp2Krl21dd93iMnfvwoz5az0fRgTHonXhTelkW/SAfE8sRAkQALaE6Ci85KxvUQouLnCTJkuzJSXCzNlO9mD86EX51DQoTAOT2VPxOndKvcf87IIPiXbmpeLzzLSES/G5J7r1QdJUVE+yeFNJEACJBCsBKjo6ms5aaZcbYdtaQTMY8qQcJGwd/toprTuMuDaf4Zh2vgRiNDYTJknzJTfHsnChowMXNOhI4a2aVtfTXmdBEiABEKSABXdMZq1IkNBwddGhHUwI35SBUxRvg3q2vOAgWs74mnzeHRv3+IYOTb+kjRTLjmSjW/EyiYnmcx4/cShMB1d+7Lx0imBBEiABIKPgM+KTs7RiE6PgyHPB/OdD5zMBRYUfQaUxuiTn2mPUHL/E2tTXuJAeEcfzZSinvmry9HiiwQkJUfiVayAuiuA+E/sYVB5LNIoQjnJ4IyrTFN5XR6rR+o9R4/V1PL46D1Hr9lEm0SktEGn6Bg82aMXkmPE6tFBELYJE3CZWE+TgQRIgAS0IOCzopNzNAqTCmBWx6q0KJq7TMe2WISfVorw1r71qtyl1X9WIOaEtJkqVgfx0UzpzKF8r9j1+/o7naeafuaINSqv+WsVpnfpFhRKrljQ+CgxBe1eewH3/utfmrKhcBIggeZLwGdFJ5EZhV9DmE49LGO4QVVy4W31UXTGODH50TeHyiZ9moJhXpLsic4Vm/3svO0W3PLsdCQkajs5vkkbhJmTAAk0OYFGKbomLz0LEHQEdosS/09Mbzjv/97AhUOHBl35WWASIIHgI0BFF3xtFpQlLhWl/jgyCS3eeB6P3Xgjl34KylZkoUkgOAlQ0QVnuwVVqRcIG/Dm66/DzS++gBYtuWh0UDUeC0sCIUCAii4EGjFQq7BPFOy/g0/AGW++jgdGjAjUYrJcJEACIU6Aii7EG7gpqlcuMv3EGIeoN6bj4dtuQ5jGk+Oboo7MkwRIIHgIUNEFT1sFRUl/E2bKtVdegRtfeRkprVsHRZlZSBIggdAmQEUX2u2rW+0OiJw+6zsAo998DQ+NG6dbvsyIBEiABOojQEVXHyFePyYBq7j6GaJhevlJPDBpEsxc4eSYvHiRBEhAfwJUdPozD5kcVwgz5YqLLsT1r72Kdh06hEy9WBESIIHQIkBFF1rtqUttDotcPul5HIa/8SoemTBBlzyZCQmQAAn4SoCKzldyzfA+u1h8+jNEwj79Edx7zz0IDxdrgTKQAAmQQIAToKIL8AYKlOL9dfAAkgf0x+k//YjOXbsGSrFYDhIgARKolwAVXb2ImncCuSPCR/v2YuA1V2PGOec0bxisPQmQQFASoKILymbTvtByA9evtv2DvOMHYfLTHyMyMlL7TJkDCZAACWhAgIpOA6hOkdYcBbYdRqCXMyY4PjcezsQ8WwWumP40evTsGRyFZilJgARIoA4CVHR1gGlMtCK2siteHoYpPR9C3oTcxojS9d6CsjJ8uHsnel5xOR698EIEw952ugJiZiRAAkFJgIrOz81Wul3BsCNn4OX7XkPHDh3x1OOPA3Kn0QAOijBTfpe6HQf7HofbP/4IMTExAVxaFo0ESIAEGkaAiq5hvOpMbctzIPrPjnjtgldw5hln1Zku0C5sPZKNH4qLcMmjD+Pyfv0CrXgsj84Ehh0/HBVWud6NDsEBtG3XToeMmmcWEREROCf5MiBLn/q37hK4a9tS0TXyGVDsQNEKA+7s9CDueeO+oJlbVlxRgQ93pqKjWNnkscsvp5mykc9BqNx+2+23Q/4xBD+BhIQEfP7+F8FfET/UIKAUnbCgwZ4VieiCoTCa3CcjG8xGlH/lQLlJ1trgXdXdTIa177GV58BWWozkXhNhrLFGY4cUAwzfHs1GvHnWDkfllRrQKSoRC/ItWHDbJ65ksi4yFBcnYFHlIZTMjWiVswFDWyfjnpNPORqr74c0U/64cwf2dO+Kmz94H/Hx8foWgLmRAAmQgM4EAkbR2QtNsGT1Q4uokTAlauvKrthtKMzfgIQ2I5Hc/TzRmxGekTVDq5oRxz4XHTuPITIBsBVmos3+H5FXcghpfe9Cz7z5HtNqHbkjNwff5OXi/PvuwUXHH691dpRPAiRAAgFBoMkVnUN4KBoOtEUSxsES30ZzKKVFu2AwG9BlyFSYI5M0zc9htyJv22z0PrIaO1qORNzJ1yPKKLqkOiu6EjHm8nHqNiSfdy4evfpqGI0eFLumJCicBEiABJqOQJMpukozZQRii8YgOq6f5mNEVmGmLC/dh+Se5yMueaDmxJWCVHRwrIXDYkPGCQ8iPlpbpVpXhX7etRP/tG+Lm2bOQFJS05ShrrIxngRIgAT0INAkis5eZITlcF9hphwFU3yUpvVURJexKG8DYlsPQvvBj8KgcW/GUZqN9tiEgrxdSLUMQXQP7ZWqJ4B78/PwpfCoPHPSXThn2DBPSRhHAiRAAs2CgK6KrtJM2QaJGI9wPcyUhXtgsDjQ6cTJsES11LRBpZky0boZMcXr8E9uCmK7XYdoo2e8uaWl2J7tf5/fI0KuVZgpZ2zZhD5i0vfD118Pk0n13tG07hROAiRAAoFMwCC88Nx8E6sXVl6yH/U4DDO5ey1mZWXh9ddf9coBUuaw8Nel2F+QgJg2AzQ3U8JaBGtBBuK6nYb49oOrV0mTY6VgBzoaNuBgrg2lSafBHHNspZpy+ANcet6pfi+LbK+kxEQMHDwYycnJfpdPgSRAAiTgbwJ2hwKpI4xCxRjlfxoEnxWdt2XZvn0Hnn71U6xOLQJMEd7e5lM6aaY0Fu1FdEpfRHc5TZgpte3NOEqPoKN5CwoPpSLDMBgx7U7wqtwjE5fj2cfv9SotE5EACZBAKBPQQ9F5tq35gWqJ2N7l1bc/xudzN8EeJnzstVZywnU/MtyC2EHCXBelrdOFnJ4QX7EZieUbsSU9AVFd/4UYk8UP1CiCBEiABEjA3wQ0UXRzf16Ip9/8HtnlcTBIJadlsBbDVJGF2K6nIqJVXy1zUmUrBTvR1bwJB7JLsTNhAmJ6pGieJzMgARIgARLwnYBfFV3qjp14Rpgp/9hWKHpw8WIszveC1Xen4rDDWCzMlMnHCTPlZTqYKXPQOfwfFJdsx0ZrH8R1GA5zfYXkdRIgARIggSYn4BdFVyq8/V5/51N8MmcD7GY9zJSZwkwZhtiB1wozZQtNIcpxv7jyzUi2b8amXdGI6HI14sJoptQUOoWTAAmQgB8JNFrR/bxgMZ5841tklcWIFUc0NlPahJmyXJopxwszpfYr7UszZY+IbTiQk4+tMWMR1YMrrfvx2aMoEiABEtCFgM+K7uDBg7hr0j34c0s2ENZQb8o6ZzTUWWmHtRRh9nxEdxqHkuzd6l9V4mPIk36rIrj+Vw9qppfn4u9otHB2RftEA+LNhVib0wNxnc5GuJZ2WLWE/I8ESIAESEALAj4ruoKCAqzYlIOolD5alKuWzLJ9i5E4dDIiEzrUuqZFxPofbkKP82ciPpmr+2vBlzJDi0CF2PYpMzNTt0qFh4dzrqhGtNPT08W8tpqdAW0yk24cbdq110Z4Nak+KzopwxgWDpNZ250GnGU1mcxiEeZ4hEXGOaM0/QyzRIqte/Spm6YVoXAS0IHA1i1bMPTSE2FM1OcH8pTk0fj1p990qFnzy+K4C7shcojHvcn8DqPwTweKVouV/TUOjVJ0GpeN4kmABIKIgKm1gvC2+uyMoaFDdxAR16aoUd0tSByvjeyaUq1impYeQZ+nUo+aMA8SIAESIAES8ECAis4DFEaRAAmQAAmEDgEqutBpS9aEBEiABEjAAwEqOg9QGEUCJEACJBA6BKjoQqctWRMSIAESIAEPBKjoPEBhFAmQAAmQQOgQoKILnbZkTUiABEiABDwQoKLzAIVRJEACJEACoUOAii502pI1IQESIAES8ECAis4DFEaRAAmQAAmEDgEqOp3b0mEtg7W0QOdcmR0JkAAJNF8CXOtSx7Yv2PcHRncrwAP3TdYxV2ZFAiRAAs2bABWdDu1fUXAISSW/48V/X4aThg/TIUdmQQIkQAIk4CRAReckocGnw1aBkj3zcdt5PXDTjS/BZDJpkAtFkgAJkAAJHIsAFd2x6DTiWlHGGgxrdxhPfjgFKSkpjZDEW0mABEiABBpDgIquMfQ83GstOozY/CV45u4LMHrUSA8pGEUCJEACJKAnASo6P9FW7DYU7p6Pm87qiNtvfRFhYUTrJ7QUQwIkQAKNIsBf40bhq7y5+OAGDG65H0/NuAPt2rXzg0SKIAESIAES8BcBKrpGkLQW5yDyyCK8cPvZOO202xohibeSAAmQAAloRYCKzgeyisOOwl2/4JrxyZh01/OwWCw+SOEtJEACJEACehCgomsg5ZLMLegTm4qn374DnTt3auDdTE4CJEACJKA3ASo6L4lbS/NhzvwFT918Gs4++0Yv72IyEiABEiCBpiZARVdPCyiKAwW7FuKyEXGY9vp0REZG1nMHL5MACZAACQQSgUYpOofDBofdqkt9pMKxiwWR5aLIugSHA6VHtqNHZCo+eu1W9OjRXZdsmQkJkAAJkIB/Cfis6Fq1aoU7Lx8Og0GfDRAiIq5EaVkaDMZ0/xKoQ1rWRWdh1OgWuPDCF0UdDXWkYjQJkAAJkECgE/BZ0SUlJeHRxx4N9PqxfCRAAjoQsNvtiNobB0O2DpmJLMI6+/zTpU8BgzgXZa0ZOQcUXWoQdjBOn3x0yYWZkAAJhDQBuWB5cXIBwtroY/2wmW0hzbMpK2ftVoKEifpY6rJ/Ktelqnwt0gUzMyGB0CdgiBI9rWh9FF3o02y6GloSzAhvrU/+YbH6PC/6qG19mDEXEiABEiABEqhFgIquFhJGkAAJkAAJhBIBKrpQak3WhQRIgARIoBYBKrpaSBhBAiRAAiQQSgSo6EKpNVkXEiABEiCBWgSo6GohYQQJkAAJkEAoEaCiC6XWZF1IgARIgARqEaCiq4WEESRAAiRAAqFEgIoulFqTdSEBEiABEqhFgIquFhJGkAAJkAAJhBIBKrpQak3WhQRIgARIoBYBKrpaSBhBAiRAAiQQSgSo6EKpNVkXEiABEiCBWgS4e0EtJIwggcYRKCoqwqKFC3XbsDc+IQFjxoxRCz3jvfdwcNu2xlXAy7sLystx3W23of+AAV7eEXzJ5s+bh4qKCl0KrigKzj3vPBiNRiz45RcsnTsXFrH9kR5h4KhROP+CC/TIqknyoKJrEuzMNJQJHDp4EDfd8yqiWvXUpZqdY7Px+1FFd3DzZjyU2FKXfP8qzEB6RkZIK7opT3yA6B76KIDDm2fj7IkTVUW3Y8cO3G2JRIsosfeRDuG5DRuo6HTgzCxIIKQImKMSER6brEudTOYSXfJpjplExiUjulV3XapujtHnBUWXygRYJhyjC7AGYXFIgARIgAT8S4CKzr88KY0ESIAESCDACHCMLsAahMUhARLwL4EDBw5g2rRpqtBhw4Zh6tSp/s2A0gKeABVdwDcRC0gCJNAYAoWFhfjqq69UEXa7nYquMTCD9F6aLoO04VhsbQk88sgjGDx4sPqXmpqqbWaUTgIkoCkB9ug0xUvhwUogLS0NG4TLtQylpaXBWg2WmwRIQBBgj46PAQmQAAmQQEgToKIL6eZl5UiABEiABKjo+AyQAAmQAAmENAGO0YV087JyJKAPAZvNBttfgNXs0CVD2yBbrXyGDx8OuXRWzSA9LZ1hzpw5aNGihfPU7fPtt9/G5Zdf7hbXHE9yZhcj5zedal6oTz5UdPpwZi4BSKBcLEosXc89BXnNGfLy8pCdne08dX2axIK7iYmJrvPmfBAWFoawIUB4W32MRCK3WrhlO+Xk5NSKrx4hF2iuK031Nq9+T3M7TjovGokX61PrzK/0cfSq/bToUz/mQgJNTuDHH3/EJZdcUm85nDsD1EzYoUMH7Nu3r2Y0z5uIQEpKCuTOETWD7G1mZmaq0ZGRkUhKSqqZRD2P0mkBZY+ZM1JTAlR0muKlcBIgAb0ILF261GNW27dvR+/evdVrZ599Nr755huP6RgZugSo6EK3bVmzegi0atUKY8eO9Zhq69atrl7AkCFDEBMTUytdcrI+uxPUypgRJEACDSJARdcgXEwcSgRGjx6NxYsXe6zSNddcg88//1y9NnPmTAwcONBjOkaSAAkEPgF9Ro4DnwNLSAIkQAIkEKIEqOhCtGFZLRIgARIggUoCVHR8EkiABEiABEKaAMfoQrp5WTkSIIG4uDhce+21KgjpWMTQ/AhQ0TW/NmeNSaBZEWjTpg0+/vjjZlVnVtadABWdOw+ekYBKQE4kP+6449Tj1q1bkwoJkEAQE6CiC+LGY9G1I3DuuedC/jGQAAkEPwE6owR/G7IGJPD/7d09b1NXAAbgY4cPRQUlYgNUoQxlaSQitVV/QYeuwMDC1p2fAUuV39Gd3wAb6gIDIxmSdGg6FByF2MYXqkhVPpoq977nxnmukiixk/Oe8xzLr+1I1wQIEDhBQNGdgOMqAgQIEDj/Al66PP97aAUEDgR2dnfLi413Bz93+c3vm5vlbpcBF3js5kTUL2f7eCN0ounJ17fnWlvRzfX2WlwNgaXl5fLop2/KYJh5b7bbN787WOYvT56U7a2tg5+7/GZ1Nvi9tbUuI6qP/fP3y+Xj5HlkHj/8eK0Mh19eZLv/4EF5u7paDr/rXjdTebyy0s3APRlV0fVkI0xjfgSak0Wvr/9aZUGrszvH5tPRjsDTZ0/bGeh/jtK8BVTz6WhHwP/o2nE0CgECBAj0VEDR9XRjTIsAAQIE2hFQdO04GoUAAQIEeiqg6Hq6MaZFgAABAu0IKLp2HI1CgAABAj0VUHQ93RjTIkCAAIF2BBRdO45GIUCAAIGeCii6nm6MaREgQIBAOwKKrh1HoxAgQIBATwWcGaWnG2NaBM6bwGSnlL3L08y0r2diLmLK7ubH8v5NZuV7f2RyFF3GWQqBuRbYH4/LV6PZuRp3Mi8SLSy66+rqBnXlz8Wy/6qr0f897uLswVHicGtJKMsgMOcClxYWyujW3+XqrUzRjWOnO57zjTtiecNv98vSwyOu6OCi7d9GHYx6eMjMrfJwrksIECBAgEBEQNFFmIUQIECAQC0BRVdLXi4BAgQIRAQUXYRZCAECBAjUElB0teTlEiBAgEBEQNFFmIUQIECAQC0BRVdLXi4BAgQIRAQUXYRZCAECBAjUElB0teTlEiBAgEBEQNFFmIUQIECAQC0BRVdLXi4BAgQIRAQUXYRZCAECBAjUElB0teTlEiBAgEBEQNFFmIUQIECAQC0BRVdLXi4BAgQIRAQUXYRZCAECBAjUElB0teTlEiBAgEBEQNFFmIUQIECAQC0BRVdLXi4BAgQIRAQUXYRZCAECBAjUElB0teTlEiBAgEBEQNFFmIUQIECAQC0BRVdLXi4BAgQIRAQUXYRZCAECBAjUElB0teTlEiBAgEBE4FIkRQgBAnMvMN4oZfTXJLLO6Z1IzIUM+fB6r+xNxpG1j95MIzmD6ew4Lqm5avzP7Xboud9xTC4nQGAmcMJdSes+gzIozYejfYHkPn6e/XRQmhIazvZz2Hzp4Dj1M7pJ5oFaB0s0JAECGYFu7qSOmvvnR+fHPkQ/6i9cdnqB3D6efk5n+83/LLr5W/LZwPw1AQIECHQg0GHZnPjSZQdLMSQBAgQIEIgK+M9blFsYAQIECKQFFF1aXB4BAgQIRAUUXZRbGAECBAikBRRdWlweAQIECEQFFF2UWxgBAgQIpAUUXVpcHgECBAhEBRRdlFsYAQIECKQFFF1aXB4BAgQIRAUUXZRbGAECBAikBRRdWlweAQIECEQFFF2UWxgBAgQIpAUUXVpcHgECBAhEBRRdlFsYAQIECKQFFF1aXB4BAgQIRAUUXZRbGAECBAikBRRdWlweAQIECEQFFF2UWxgBAgQIpAUUXVpcHgECBAhEBRRdlFsYAQIECKQFFF1aXB4BAgQIRAUUXZRbGAECBAikBRRdWlweAQIECEQFFF2UWxgBAgQIpAUUXVpcHgECBAhEBT4B4Qne9nVqU/wAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<IPython.core.display.Image object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from IPython.display import Image, display\n",
    "display(Image(filename='./images/tensor_contraction.png', embed=True))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example\n",
    "\n",
    "In TensorLy, all the tensor algebra functions are located in the `tensorly.tenalg` module. For the n-mode product, you will need to use the function `mode_dot` that works transparently for multiplying a tensor by a matrix or a vector along a given mode.\n",
    "\n",
    "#### Tensor times matrix\n",
    "\n",
    "With the tensor $\\tilde{\\mathcal{X}}$ of size (3, 4, 2) we defined previously, let's define a matrix $\\mathbf{M}$ of size (5, 4) to multiply along the second mode:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(5, 4)\n"
     ]
    }
   ],
   "source": [
    "M = tl.tensor(np.arange(4*5).reshape((5, 4)))\n",
    "print(M.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Keep in mind indexing starts at zero, so the second mode is represented by `mode=1`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "res = tl.tenalg.mode_dot(X, M, mode=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As expected the result is of shape (3, 5, 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(3, 5, 2)"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Tensor times vector\n",
    "\n",
    "Similarly, we can contract along the mode 1 with a vector of size 4 (our tensor is of size (3, 4, 2).\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(4,)\n"
     ]
    }
   ],
   "source": [
    "v = tl.tensor(np.arange(4))\n",
    "print(v.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0, 1, 2, 3])"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "v"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[[ 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]]])"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "res = tl.tenalg.mode_dot(X, v, mode=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Recall that $\\tilde{\\mathcal{X}}$ is of shape (3, 4, 2). We are here contracting over its 2nd mode. \n",
    "\n",
    "Since we have multiplied by a vector, we have effectively contracted out one mode of the tensor so the result will be of matrix of size (3, 2). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(3, 2)"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 28,  34],\n",
       "       [ 76,  82],\n",
       "       [124, 130]])"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, we could form each column (fiber) of the result, by taking the dot product between the frontal slices of $\\tilde{\\mathcal{X}}$ and the vector:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 28,  76, 124])"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X[:, :, 0] @ v"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 34,  82, 130])"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X[:, :, 1] @ v"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can form the whole product by concanating these as the column of the resulting matrix:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 28,  34],\n",
       "       [ 76,  82],\n",
       "       [124, 130]])"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.stack([X[:, :, 0] @ v, X[:, :, 1] @ v]).T"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We could equivalently use the mode-0 slices:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[28 34]\n",
      "[76 82]\n",
      "[124 130]\n"
     ]
    }
   ],
   "source": [
    "for i in range(3):\n",
    "    print(X[i, ...].T @ v)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Again, we can form the whole result by stacking these as the rows of the matrix:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 28,  34],\n",
       "       [ 76,  82],\n",
       "       [124, 130]])"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.vstack([X[i, ...].T @ v for i in range(3)])"
   ]
  }
 ],
 "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.8.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}