{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# About this Notebook\n", "\n", "**Bayesian Temporal Tensor Factorization** (or **BTTF** for short) is a type of Bayesian tensor decomposition that achieves state-of-the-art results on challenging the missing data imputation problem. In the following, we will discuss:\n", "\n", "- What the BTTF is?\n", "\n", "- How to implement BTTF mainly using Python Numpy with high efficiency?\n", "\n", "- How to make imputations with real-world spatiotemporal datasets?\n", "\n", "If you want to understand what is BTMF and its modeling tricks in detail, our paper is for you:\n", "\n", "> Xinyu Chen, Lijun Sun (2019). **Bayesian temporal factorization for multidimensional time series prediction**.\n", "\n", "## Quick Run\n", "\n", "This notebook is publicly available for any usage at our data imputation project. Please click [**transdim**](https://github.com/xinychen/transdim)." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "from numpy.random import multivariate_normal as mvnrnd\n", "from scipy.stats import wishart\n", "from scipy.stats import invwishart\n", "from numpy.linalg import inv as inv" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Part 1: Matrix Computation Concepts\n", "\n", "## 1) Kronecker product\n", "\n", "- **Definition**:\n", "\n", "Given two matrices $A\\in\\mathbb{R}^{m_1\\times n_1}$ and $B\\in\\mathbb{R}^{m_2\\times n_2}$, then, the **Kronecker product** between these two matrices is defined as\n", "\n", "$$A\\otimes B=\\left[ \\begin{array}{cccc} a_{11}B & a_{12}B & \\cdots & a_{1m_2}B \\\\ a_{21}B & a_{22}B & \\cdots & a_{2m_2}B \\\\ \\vdots & \\vdots & \\ddots & \\vdots \\\\ a_{m_11}B & a_{m_12}B & \\cdots & a_{m_1m_2}B \\\\ \\end{array} \\right]$$\n", "where the symbol $\\otimes$ denotes Kronecker product, and the size of resulted $A\\otimes B$ is $(m_1m_2)\\times (n_1n_2)$ (i.e., $m_1\\times m_2$ columns and $n_1\\times n_2$ rows).\n", "\n", "- **Example**:\n", "\n", "If $A=\\left[ \\begin{array}{cc} 1 & 2 \\\\ 3 & 4 \\\\ \\end{array} \\right]$ and $B=\\left[ \\begin{array}{ccc} 5 & 6 & 7\\\\ 8 & 9 & 10 \\\\ \\end{array} \\right]$, then, we have\n", "\n", "$$A\\otimes B=\\left[ \\begin{array}{cc} 1\\times \\left[ \\begin{array}{ccc} 5 & 6 & 7\\\\ 8 & 9 & 10\\\\ \\end{array} \\right] & 2\\times \\left[ \\begin{array}{ccc} 5 & 6 & 7\\\\ 8 & 9 & 10\\\\ \\end{array} \\right] \\\\ 3\\times \\left[ \\begin{array}{ccc} 5 & 6 & 7\\\\ 8 & 9 & 10\\\\ \\end{array} \\right] & 4\\times \\left[ \\begin{array}{ccc} 5 & 6 & 7\\\\ 8 & 9 & 10\\\\ \\end{array} \\right] \\\\ \\end{array} \\right]$$\n", "\n", "$$=\\left[ \\begin{array}{cccccc} 5 & 6 & 7 & 10 & 12 & 14 \\\\ 8 & 9 & 10 & 16 & 18 & 20 \\\\ 15 & 18 & 21 & 20 & 24 & 28 \\\\ 24 & 27 & 30 & 32 & 36 & 40 \\\\ \\end{array} \\right]\\in\\mathbb{R}^{4\\times 6}.$$\n", "\n", "## 2) Khatri-Rao product (kr_prod)\n", "\n", "- **Definition**:\n", "\n", "Given two matrices $A=\\left( \\boldsymbol{a}_1,\\boldsymbol{a}_2,...,\\boldsymbol{a}_r \\right)\\in\\mathbb{R}^{m\\times r}$ and $B=\\left( \\boldsymbol{b}_1,\\boldsymbol{b}_2,...,\\boldsymbol{b}_r \\right)\\in\\mathbb{R}^{n\\times r}$ with same number of columns, then, the **Khatri-Rao product** (or **column-wise Kronecker product**) between $A$ and $B$ is given as follows,\n", "\n", "$$A\\odot B=\\left( \\boldsymbol{a}_1\\otimes \\boldsymbol{b}_1,\\boldsymbol{a}_2\\otimes \\boldsymbol{b}_2,...,\\boldsymbol{a}_r\\otimes \\boldsymbol{b}_r \\right)\\in\\mathbb{R}^{(mn)\\times r},$$\n", "where the symbol $\\odot$ denotes Khatri-Rao product, and $\\otimes$ denotes Kronecker product.\n", "\n", "- **Example**:\n", "\n", "If $A=\\left[ \\begin{array}{cc} 1 & 2 \\\\ 3 & 4 \\\\ \\end{array} \\right]=\\left( \\boldsymbol{a}_1,\\boldsymbol{a}_2 \\right)$ and $B=\\left[ \\begin{array}{cc} 5 & 6 \\\\ 7 & 8 \\\\ 9 & 10 \\\\ \\end{array} \\right]=\\left( \\boldsymbol{b}_1,\\boldsymbol{b}_2 \\right)$, then, we have\n", "\n", "$$A\\odot B=\\left( \\boldsymbol{a}_1\\otimes \\boldsymbol{b}_1,\\boldsymbol{a}_2\\otimes \\boldsymbol{b}_2 \\right)$$\n", "\n", "$$=\\left[ \\begin{array}{cc} \\left[ \\begin{array}{c} 1 \\\\ 3 \\\\ \\end{array} \\right]\\otimes \\left[ \\begin{array}{c} 5 \\\\ 7 \\\\ 9 \\\\ \\end{array} \\right] & \\left[ \\begin{array}{c} 2 \\\\ 4 \\\\ \\end{array} \\right]\\otimes \\left[ \\begin{array}{c} 6 \\\\ 8 \\\\ 10 \\\\ \\end{array} \\right] \\\\ \\end{array} \\right]$$\n", "\n", "$$=\\left[ \\begin{array}{cc} 5 & 12 \\\\ 7 & 16 \\\\ 9 & 20 \\\\ 15 & 24 \\\\ 21 & 32 \\\\ 27 & 40 \\\\ \\end{array} \\right]\\in\\mathbb{R}^{6\\times 2}.$$" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "def kr_prod(a, b):\n", " return np.einsum('ir, jr -> ijr', a, b).reshape(a.shape[0] * b.shape[0], -1)" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[ 5 12]\n", " [ 7 16]\n", " [ 9 20]\n", " [15 24]\n", " [21 32]\n", " [27 40]]\n" ] } ], "source": [ "A = np.array([[1, 2], [3, 4]])\n", "B = np.array([[5, 6], [7, 8], [9, 10]])\n", "print(kr_prod(A, B))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 3) Computing Covariance Matrix (cov_mat)\n", "\n", "For any matrix $X\\in\\mathbb{R}^{m\\times n}$, cov_mat can return a $n\\times n$ covariance matrix for special use in the following." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "def cov_mat(mat):\n", " dim1, dim2 = mat.shape\n", " new_mat = np.zeros((dim2, dim2))\n", " mat_bar = np.mean(mat, axis = 0)\n", " for i in range(dim1):\n", " new_mat += np.einsum('i, j -> ij', mat[i, :] - mat_bar, mat[i, :] - mat_bar)\n", " return new_mat" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 4) CP decomposition (cp_combine)\n", "\n", "- **Definition**:\n", "\n", "The CP decomposition factorizes a tensor into a sum of outer products of vectors. For example, for a third-order tensor $\\mathcal{Y}\\in\\mathbb{R}^{m\\times n\\times f}$, the CP decomposition can be written as\n", "\n", "$$\\hat{\\mathcal{Y}}=\\sum_{s=1}^{r}\\boldsymbol{u}_{s}\\circ\\boldsymbol{v}_{s}\\circ\\boldsymbol{x}_{s},$$\n", "or element-wise,\n", "\n", "$$\\hat{y}_{ijt}=\\sum_{s=1}^{r}u_{is}v_{js}x_{ts},\\forall (i,j,t),$$\n", "where vectors $\\boldsymbol{u}_{s}\\in\\mathbb{R}^{m},\\boldsymbol{v}_{s}\\in\\mathbb{R}^{n},\\boldsymbol{x}_{s}\\in\\mathbb{R}^{f}$ are columns of factor matrices $U\\in\\mathbb{R}^{m\\times r},V\\in\\mathbb{R}^{n\\times r},X\\in\\mathbb{R}^{f\\times r}$, respectively. The symbol $\\circ$ denotes vector outer product.\n", "\n", "- **Example**:\n", "\n", "Given matrices $U=\\left[ \\begin{array}{cc} 1 & 2 \\\\ 3 & 4 \\\\ \\end{array} \\right]\\in\\mathbb{R}^{2\\times 2}$, $V=\\left[ \\begin{array}{cc} 1 & 2 \\\\ 3 & 4 \\\\ 5 & 6 \\\\ \\end{array} \\right]\\in\\mathbb{R}^{3\\times 2}$ and $X=\\left[ \\begin{array}{cc} 1 & 5 \\\\ 2 & 6 \\\\ 3 & 7 \\\\ 4 & 8 \\\\ \\end{array} \\right]\\in\\mathbb{R}^{4\\times 2}$, then if $\\hat{\\mathcal{Y}}=\\sum_{s=1}^{r}\\boldsymbol{u}_{s}\\circ\\boldsymbol{v}_{s}\\circ\\boldsymbol{x}_{s}$, then, we have\n", "\n", "$$\\hat{Y}_1=\\hat{\\mathcal{Y}}(:,:,1)=\\left[ \\begin{array}{ccc} 31 & 42 & 65 \\\\ 63 & 86 & 135 \\\\ \\end{array} \\right],$$\n", "$$\\hat{Y}_2=\\hat{\\mathcal{Y}}(:,:,2)=\\left[ \\begin{array}{ccc} 38 & 52 & 82 \\\\ 78 & 108 & 174 \\\\ \\end{array} \\right],$$\n", "$$\\hat{Y}_3=\\hat{\\mathcal{Y}}(:,:,3)=\\left[ \\begin{array}{ccc} 45 & 62 & 99 \\\\ 93 & 130 & 213 \\\\ \\end{array} \\right],$$\n", "$$\\hat{Y}_4=\\hat{\\mathcal{Y}}(:,:,4)=\\left[ \\begin{array}{ccc} 52 & 72 & 116 \\\\ 108 & 152 & 252 \\\\ \\end{array} \\right].$$" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "def cp_combine(U, V, X):\n", " return np.einsum('is, js, ts -> ijt', U, V, X)" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[[[ 31 38 45 52]\n", " [ 42 52 62 72]\n", " [ 65 82 99 116]]\n", "\n", " [[ 63 78 93 108]\n", " [ 86 108 130 152]\n", " [135 174 213 252]]]\n", "\n", "tensor size:\n", "(2, 3, 4)\n" ] } ], "source": [ "U = np.array([[1, 2], [3, 4]])\n", "V = np.array([[1, 3], [2, 4], [5, 6]])\n", "X = np.array([[1, 5], [2, 6], [3, 7], [4, 8]])\n", "print(cp_combine(U, V, X))\n", "print()\n", "print('tensor size:')\n", "print(cp_combine(U, V, X).shape)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 5) Tensor Unfolding (ten2mat) and Matrix Folding (mat2ten)\n", "\n", "Using numpy reshape to perform 3rd rank tensor unfold operation. [[**link**](https://stackoverflow.com/questions/49970141/using-numpy-reshape-to-perform-3rd-rank-tensor-unfold-operation)]" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "def ten2mat(tensor, mode):\n", " return np.reshape(np.moveaxis(tensor, mode, 0), (tensor.shape[mode], -1), order = 'F')" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "tensor size:\n", "(3, 2, 4)\n", "original tensor:\n", "[[[ 1 2 3 4]\n", " [ 3 4 5 6]]\n", "\n", " [[ 5 6 7 8]\n", " [ 7 8 9 10]]\n", "\n", " [[ 9 10 11 12]\n", " [11 12 13 14]]]\n", "\n", "(1) mode-1 tensor unfolding:\n", "[[ 1 3 2 4 3 5 4 6]\n", " [ 5 7 6 8 7 9 8 10]\n", " [ 9 11 10 12 11 13 12 14]]\n", "\n", "(2) mode-2 tensor unfolding:\n", "[[ 1 5 9 2 6 10 3 7 11 4 8 12]\n", " [ 3 7 11 4 8 12 5 9 13 6 10 14]]\n", "\n", "(3) mode-3 tensor unfolding:\n", "[[ 1 5 9 3 7 11]\n", " [ 2 6 10 4 8 12]\n", " [ 3 7 11 5 9 13]\n", " [ 4 8 12 6 10 14]]\n" ] } ], "source": [ "X = np.array([[[1, 2, 3, 4], [3, 4, 5, 6]], [[5, 6, 7, 8], [7, 8, 9, 10]], [[9, 10, 11, 12], [11, 12, 13, 14]]])\n", "print('tensor size:')\n", "print(X.shape)\n", "print('original tensor:')\n", "print(X)\n", "print()\n", "print('(1) mode-1 tensor unfolding:')\n", "print(ten2mat(X, 0))\n", "print()\n", "print('(2) mode-2 tensor unfolding:')\n", "print(ten2mat(X, 1))\n", "print()\n", "print('(3) mode-3 tensor unfolding:')\n", "print(ten2mat(X, 2))" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "def mat2ten(mat, tensor_size, mode):\n", " index = list()\n", " index.append(mode)\n", " for i in range(tensor_size.shape[0]):\n", " if i != mode:\n", " index.append(i)\n", " return np.moveaxis(np.reshape(mat, list(tensor_size[index]), order = 'F'), 0, mode)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 6) Generating Matrix Normal Distributed Random Matrix" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "def mnrnd(M, U, V):\n", " \"\"\"\n", " Generate matrix normal distributed random matrix.\n", " M is a m-by-n matrix, U is a m-by-m matrix, and V is a n-by-n matrix.\n", " \"\"\"\n", " dim1, dim2 = M.shape\n", " X0 = np.random.rand(dim1, dim2)\n", " P = np.linalg.cholesky(U)\n", " Q = np.linalg.cholesky(V)\n", " return M + np.matmul(np.matmul(P, X0), Q.T)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Part 2: Bayesian Temporal Tensor Factorization (BTTF)" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "def BTTF(dense_tensor, sparse_tensor, init, rank, time_lags, maxiter1, maxiter2):\n", " \"\"\"Bayesian Temporal Tensor Factorization, BTTF.\"\"\"\n", " U = init[\"U\"]\n", " V = init[\"V\"]\n", " X = init[\"X\"]\n", " \n", " d = time_lags.shape[0]\n", " dim1, dim2, dim3 = sparse_tensor.shape\n", " dim = np.array([dim1, dim2, dim3])\n", " pos = np.where((dense_tensor != 0) & (sparse_tensor == 0))\n", " position = np.where(sparse_tensor != 0)\n", " binary_tensor = np.zeros((dim1, dim2, dim3))\n", " binary_tensor[position] = 1\n", "\n", " beta0 = 1\n", " nu0 = rank\n", " mu0 = np.zeros((rank))\n", " W0 = np.eye(rank)\n", " tau = 1\n", " alpha = 1e-6\n", " beta = 1e-6\n", " S0 = np.eye(rank)\n", " Psi0 = np.eye(rank * d)\n", " M0 = np.zeros((rank * d, rank))\n", " \n", " mat_hat = np.zeros((dim1, dim2, dim3 + 1))\n", " U_plus = np.zeros((dim1, rank))\n", " V_plus = np.zeros((dim2, rank))\n", " X_plus = np.zeros((dim3, rank))\n", " X_new = np.zeros((dim3 + 1, rank))\n", " X_new_plus = np.zeros((dim3 + 1, rank))\n", " A_plus = np.zeros((rank, rank, d))\n", " tensor_hat_plus = np.zeros((dim1, dim2, dim3 + 1))\n", " for iters in range(maxiter1):\n", " for order in range(2):\n", " if order == 0:\n", " mat = U.copy()\n", " elif order == 1:\n", " mat = V.copy()\n", " mat_bar = np.mean(mat, axis = 0)\n", " var_mu_hyper = (dim[order] * mat_bar + beta0 * mu0)/(dim[order] + beta0)\n", " var_W_hyper = inv(inv(W0) + cov_mat(mat) + dim[order] * beta0/(dim[order] + beta0)\n", " * np.outer(mat_bar - mu0, mat_bar - mu0))\n", " var_Lambda_hyper = wishart(df = dim[order] + nu0, scale = var_W_hyper, seed = None).rvs()\n", " var_mu_hyper = mvnrnd(var_mu_hyper, inv((dim[order] + beta0) * var_Lambda_hyper))\n", "\n", " if order == 0:\n", " var1 = kr_prod(X, V).T\n", " elif order == 1:\n", " var1 = kr_prod(X, U).T\n", " var2 = kr_prod(var1, var1)\n", " var3 = (tau * np.matmul(var2, ten2mat(binary_tensor, order).T).reshape([rank, rank, dim[order]])\n", " + np.dstack([var_Lambda_hyper] * dim[order]))\n", " var4 = (tau * np.matmul(var1, ten2mat(sparse_tensor, order).T) \n", " + np.dstack([np.matmul(var_Lambda_hyper, var_mu_hyper)] * dim[order])[0, :, :])\n", " for i in range(dim[order]):\n", " inv_var_Lambda = inv(var3[ :, :, i])\n", " vec = mvnrnd(np.matmul(inv_var_Lambda, var4[:, i]), inv_var_Lambda)\n", " if order == 0:\n", " U[i, :] = vec.copy()\n", " elif order == 1:\n", " V[i, :] = vec.copy()\n", "\n", " Z_mat = X[np.max(time_lags) : dim3, :]\n", " Q_mat = np.zeros((dim3 - np.max(time_lags), rank * d))\n", " for t in range(np.max(time_lags), dim3):\n", " Q_mat[t - np.max(time_lags), :] = X[t - time_lags, :].reshape([rank * d])\n", " var_Psi = inv(inv(Psi0) + np.matmul(Q_mat.T, Q_mat))\n", " var_M = np.matmul(var_Psi, np.matmul(inv(Psi0), M0) + np.matmul(Q_mat.T, Z_mat))\n", " var_S = (S0 + np.matmul(Z_mat.T, Z_mat) + np.matmul(np.matmul(M0.T, inv(Psi0)), M0) \n", " - np.matmul(np.matmul(var_M.T, inv(var_Psi)), var_M))\n", " Sigma = invwishart(df = nu0 + dim3 - np.max(time_lags), scale = var_S, seed = None).rvs()\n", " A = mat2ten(mnrnd(var_M, var_Psi, Sigma).T, np.array([rank, rank, d]), 0)\n", " if iters + 1 > maxiter1 - maxiter2:\n", " A_plus += A\n", "\n", " Lambda_x = inv(Sigma)\n", " var1 = kr_prod(V, U).T\n", " var2 = kr_prod(var1, var1)\n", " var3 = (tau * np.matmul(var2, ten2mat(binary_tensor, 2).T).reshape([rank, rank, dim3]) \n", " + np.dstack([Lambda_x] * dim3))\n", " var4 = tau * np.matmul(var1, ten2mat(sparse_tensor, 2).T)\n", " for t in range(dim3):\n", " Mt = np.zeros((rank, rank))\n", " Nt = np.zeros(rank)\n", " if t < np.max(time_lags):\n", " Qt = np.zeros(rank)\n", " else:\n", " Qt = np.matmul(Lambda_x, np.matmul(ten2mat(A, 0), X[t - time_lags, :].reshape([rank * d])))\n", " if t < dim3 - np.min(time_lags):\n", " if t >= np.max(time_lags) and t < dim3 - np.max(time_lags):\n", " index = list(range(0, d))\n", " else:\n", " index = list(np.where((t + time_lags >= np.max(time_lags)) & (t + time_lags < dim3)))[0]\n", " for k in index:\n", " Ak = A[:, :, k]\n", " Mt += np.matmul(np.matmul(Ak.T, Lambda_x), Ak)\n", " A0 = A.copy()\n", " A0[:, :, k] = 0\n", " var5 = (X[t + time_lags[k], :] \n", " - np.matmul(ten2mat(A0, 0), X[t + time_lags[k] - time_lags, :].reshape([rank * d])))\n", " Nt += np.matmul(np.matmul(Ak.T, Lambda_x), var5)\n", " var_mu = var4[:, t] + Nt + Qt\n", " if t < np.max(time_lags):\n", " inv_var_Lambda = inv(var3[:, :, t] + Mt - Lambda_x + np.eye(rank))\n", " else:\n", " inv_var_Lambda = inv(var3[:, :, t] + Mt)\n", " X[t, :] = mvnrnd(np.matmul(inv_var_Lambda, var_mu), inv_var_Lambda)\n", "\n", " if iters + 1 > maxiter1 - maxiter2:\n", " U_plus += U\n", " V_plus += V\n", " X_plus += X\n", "\n", " tensor_hat = cp_combine(U, V, X)\n", " if iters + 1 > maxiter1 - maxiter2:\n", " X_new[0 : dim3, :] = X.copy()\n", " X_new[dim3, :] = np.matmul(ten2mat(A, 0), X_new[dim3 - time_lags, :].reshape([rank * d]))\n", " X_new_plus += X_new\n", " tensor_hat_plus += cp_combine(U, V, X_new)\n", " \n", " tau = np.random.gamma(alpha + 0.5 * sparse_tensor[position].shape[0], \n", " 1/(beta + 0.5 * np.sum((sparse_tensor - tensor_hat)[position] ** 2)))\n", " rmse = np.sqrt(np.sum((dense_tensor[pos] - tensor_hat[pos]) ** 2)/dense_tensor[pos].shape[0])\n", " if (iters + 1) % 200 == 0 and iters < maxiter1 - maxiter2:\n", " print('Iter: {}'.format(iters + 1))\n", " print('RMSE: {:.6}'.format(rmse))\n", " print()\n", "\n", " U = U_plus/maxiter2\n", " V = V_plus/maxiter2\n", " X = X_plus/maxiter2\n", " X_new = X_new_plus/maxiter2\n", " A = A_plus/maxiter2\n", " tensor_hat = tensor_hat_plus/maxiter2\n", " if maxiter1 >= 100:\n", " final_mape = np.sum(np.abs(dense_tensor[pos] \n", " - tensor_hat[pos])/dense_tensor[pos])/dense_tensor[pos].shape[0]\n", " final_rmse = np.sqrt(np.sum((dense_tensor[pos] - tensor_hat[pos]) ** 2)/dense_tensor[pos].shape[0])\n", " print('Imputation MAPE: {:.6}'.format(final_mape))\n", " print('Imputation RMSE: {:.6}'.format(final_rmse))\n", " print()\n", " \n", " return tensor_hat, U, V, X_new, A" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [], "source": [ "def OnlineBTTF(sparse_mat, init, time_lags, maxiter1, maxiter2):\n", " \"\"\"Online Bayesain Temporal Tensor Factorization\"\"\"\n", " U = init[\"U\"]\n", " V = init[\"V\"]\n", " X = init[\"X\"]\n", " A = init[\"A\"]\n", " \n", " d = time_lags.shape[0]\n", " dim1, dim2 = sparse_mat.shape\n", " dim3 = 1\n", " sparse_tensor = sparse_mat.reshape([dim1, dim2, dim3])\n", " t = X.shape[0]\n", " rank = X.shape[1]\n", " position = np.where(sparse_mat != 0)\n", " binary_tensor = np.zeros((dim1, dim2, dim3))\n", " binary_tensor[position] = 1\n", "\n", " tau = 1\n", " alpha = 1e-6\n", " beta = 1e-6\n", " nu0 = rank\n", " W0 = np.eye(rank)\n", " var_mu0 = np.matmul(ten2mat(A, 0), X[t - 1 - time_lags, :].reshape([rank * d]))\n", "\n", " X_new = np.zeros((t + 1, rank))\n", " X_new_plus = np.zeros((t + 1, rank))\n", " tensor_hat_plus = np.zeros((U.shape[0], V.shape[0], t + 1))\n", " for iters in range(maxiter1):\n", " vec0 = X[t - 1, :] - var_mu0\n", " Lambda_x = wishart(df = nu0 + 1, scale = inv(inv(W0) + np.outer(vec0, vec0)), seed = None).rvs()\n", " \n", " var1 = kr_prod(V, U).T\n", " var2 = kr_prod(var1, var1)\n", " var_mu = tau * np.matmul(var1, ten2mat(sparse_tensor, 2).reshape([dim1 * dim2])) + np.matmul(Lambda_x, var_mu0)\n", " inv_var_Lambda = inv(tau * np.matmul(var2, ten2mat(binary_tensor, 2).reshape([dim1 * dim2])).reshape([rank, rank]) + Lambda_x)\n", " X[t - 1, :] = mvnrnd(np.matmul(inv_var_Lambda, var_mu), inv_var_Lambda)\n", " \n", " mat_hat = np.einsum('ir, jr, r -> ij', U, V, X[t - 1, :])\n", " \n", " tau = np.random.gamma(alpha + 0.5 * sparse_mat[position].shape[0], \n", " 1/(beta + 0.5 * np.sum((sparse_mat - mat_hat)[position] ** 2)))\n", " \n", " if iters + 1 > maxiter1 - maxiter2:\n", " X_new[0 : t, :] = X.copy()\n", " X_new[t, :] = np.matmul(ten2mat(A, 0), X_new[t - time_lags, :].reshape([rank * d]))\n", " X_new_plus += X_new\n", " tensor_hat_plus += cp_combine(U, V, X_new)\n", "\n", " X_new = X_new_plus/maxiter2\n", " tensor_hat = tensor_hat_plus/maxiter2\n", "\n", " return tensor_hat, X_new" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "def st_prediction(dense_tensor, sparse_tensor, pred_time_steps, rank, time_lags, maxiter):\n", " start_time = dense_tensor.shape[2] - pred_time_steps\n", " dense_tensor0 = dense_tensor[:, :, 0 : start_time]\n", " sparse_tensor0 = sparse_tensor[:, :, 0 : start_time]\n", " dim1 = sparse_tensor0.shape[0]\n", " dim2 = sparse_tensor0.shape[1]\n", " dim3 = sparse_tensor0.shape[2]\n", " tensor_hat = np.zeros((dim1, dim2, pred_time_steps))\n", " \n", " for t in range(pred_time_steps):\n", " if t == 0:\n", " init = {\"U\": 0.1 * np.random.rand(dim1, rank), \"V\": 0.1 * np.random.rand(dim2, rank),\n", " \"X\": 0.1 * np.random.rand(dim3, rank)}\n", " tensor, U, V, X, A = BTTF(dense_tensor0, sparse_tensor0, init, rank, time_lags, maxiter[0], maxiter[1])\n", " X0 = X.copy()\n", " else:\n", " sparse_tensor1 = sparse_tensor[:, :, 0 : start_time + t]\n", " init = {\"U\": U, \"V\": V, \"X\": X0, \"A\": A}\n", " tensor, X = OnlineBTTF(sparse_tensor1[:, :, -1], init, time_lags, maxiter[2], maxiter[3])\n", " X0 = X.copy()\n", " tensor_hat[:, :, t] = tensor[:, :, -1]\n", " if (t + 1) % 40 == 0:\n", " print('Time step: {}'.format(t + 1))\n", "\n", " small_dense_tensor = dense_tensor[:, :, start_time : dense_tensor.shape[2]]\n", " pos = np.where(small_dense_tensor != 0)\n", " final_mape = np.sum(np.abs(small_dense_tensor[pos] - \n", " tensor_hat[pos])/small_dense_tensor[pos])/small_dense_tensor[pos].shape[0]\n", " final_rmse = np.sqrt(np.sum((small_dense_tensor[pos] - \n", " tensor_hat[pos]) ** 2)/small_dense_tensor[pos].shape[0])\n", " print('Final MAPE: {:.6}'.format(final_mape))\n", " print('Final RMSE: {:.6}'.format(final_rmse))\n", " print()\n", " return tensor_hat" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**How to transform a data set into something we can use for missing data imputation?**\n" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [], "source": [ "import scipy.io\n", "\n", "tensor = scipy.io.loadmat('../NYC-data-set/tensor.mat')\n", "dense_tensor = tensor['tensor']\n", "rm_tensor = scipy.io.loadmat('../NYC-data-set/rm_tensor.mat')\n", "rm_tensor = rm_tensor['rm_tensor']\n", "nm_tensor = scipy.io.loadmat('../NYC-data-set/nm_tensor.mat')\n", "nm_tensor = nm_tensor['nm_tensor']\n", "\n", "missing_rate = 0.0\n", "\n", "# =============================================================================\n", "### Random missing (RM) scenario\n", "### Set the RM scenario by:\n", "binary_tensor = np.round(rm_tensor + 0.5 - missing_rate)\n", "# =============================================================================\n", "\n", "sparse_tensor = np.multiply(dense_tensor, binary_tensor)" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "name": "stderr", "output_type": "stream", "text": [ "/Users/xinyuchen/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:124: RuntimeWarning: invalid value encountered in double_scalars\n", "/Users/xinyuchen/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:138: RuntimeWarning: invalid value encountered in double_scalars\n", "/Users/xinyuchen/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:139: RuntimeWarning: invalid value encountered in double_scalars\n" ] }, { "name": "stdout", "output_type": "stream", "text": [ "Imputation MAPE: nan\n", "Imputation RMSE: nan\n", "\n", "Time step: 40\n", "Time step: 80\n", "Time step: 120\n", "Time step: 160\n", "Final MAPE: 0.587582\n", "Final RMSE: 5.30109\n", "\n", "Running time: 1432 seconds\n" ] } ], "source": [ "import time\n", "start = time.time()\n", "pred_time_steps = 24 * 7\n", "rank = 30\n", "time_lags = np.array([1, 2, 24])\n", "maxiter = np.array([200, 100, 200, 100])\n", "small_sparse_tensor = sparse_tensor[:, :, dense_tensor.shape[2] - pred_time_steps : dense_tensor.shape[2]]\n", "small_dense_tensor = dense_tensor[:, :, dense_tensor.shape[2] - pred_time_steps : dense_tensor.shape[2]]\n", "tensor_hat = st_prediction(dense_tensor, sparse_tensor, pred_time_steps, rank, time_lags, maxiter)\n", "end = time.time()\n", "print('Running time: %d seconds'%(end - start))" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [], "source": [ "import scipy.io\n", "\n", "tensor = scipy.io.loadmat('../NYC-data-set/tensor.mat')\n", "dense_tensor = tensor['tensor']\n", "rm_tensor = scipy.io.loadmat('../NYC-data-set/rm_tensor.mat')\n", "rm_tensor = rm_tensor['rm_tensor']\n", "nm_tensor = scipy.io.loadmat('../NYC-data-set/nm_tensor.mat')\n", "nm_tensor = nm_tensor['nm_tensor']\n", "\n", "missing_rate = 0.1\n", "\n", "# =============================================================================\n", "### Random missing (RM) scenario\n", "### Set the RM scenario by:\n", "binary_tensor = np.round(rm_tensor + 0.5 - missing_rate)\n", "# =============================================================================\n", "\n", "sparse_tensor = np.multiply(dense_tensor, binary_tensor)" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Imputation MAPE: 0.533769\n", "Imputation RMSE: 4.811\n", "\n", "Time step: 40\n", "Time step: 80\n", "Time step: 120\n", "Time step: 160\n", "Final MAPE: 0.586588\n", "Final RMSE: 5.29504\n", "\n", "Running time: 1595 seconds\n" ] } ], "source": [ "import time\n", "start = time.time()\n", "pred_time_steps = 24 * 7\n", "rank = 30\n", "time_lags = np.array([1, 2, 24])\n", "maxiter = np.array([200, 100, 200, 100])\n", "small_sparse_tensor = sparse_tensor[:, :, dense_tensor.shape[2] - pred_time_steps : dense_tensor.shape[2]]\n", "small_dense_tensor = dense_tensor[:, :, dense_tensor.shape[2] - pred_time_steps : dense_tensor.shape[2]]\n", "tensor_hat = st_prediction(dense_tensor, sparse_tensor, pred_time_steps, rank, time_lags, maxiter)\n", "end = time.time()\n", "print('Running time: %d seconds'%(end - start))" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [], "source": [ "import scipy.io\n", "\n", "tensor = scipy.io.loadmat('../NYC-data-set/tensor.mat')\n", "dense_tensor = tensor['tensor']\n", "rm_tensor = scipy.io.loadmat('../NYC-data-set/rm_tensor.mat')\n", "rm_tensor = rm_tensor['rm_tensor']\n", "nm_tensor = scipy.io.loadmat('../NYC-data-set/nm_tensor.mat')\n", "nm_tensor = nm_tensor['nm_tensor']\n", "\n", "missing_rate = 0.3\n", "\n", "# =============================================================================\n", "### Random missing (RM) scenario\n", "### Set the RM scenario by:\n", "binary_tensor = np.round(rm_tensor + 0.5 - missing_rate)\n", "# =============================================================================\n", "\n", "sparse_tensor = np.multiply(dense_tensor, binary_tensor)" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Imputation MAPE: 0.531925\n", "Imputation RMSE: 4.92126\n", "\n", "Time step: 40\n", "Time step: 80\n", "Time step: 120\n", "Time step: 160\n", "Final MAPE: 0.587577\n", "Final RMSE: 5.39982\n", "\n", "Running time: 1087 seconds\n" ] } ], "source": [ "import time\n", "start = time.time()\n", "pred_time_steps = 24 * 7\n", "rank = 30\n", "time_lags = np.array([1, 2, 24])\n", "maxiter = np.array([200, 100, 200, 100])\n", "small_sparse_tensor = sparse_tensor[:, :, dense_tensor.shape[2] - pred_time_steps : dense_tensor.shape[2]]\n", "small_dense_tensor = dense_tensor[:, :, dense_tensor.shape[2] - pred_time_steps : dense_tensor.shape[2]]\n", "tensor_hat = st_prediction(dense_tensor, sparse_tensor, pred_time_steps, rank, time_lags, maxiter)\n", "end = time.time()\n", "print('Running time: %d seconds'%(end - start))" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [], "source": [ "import scipy.io\n", "\n", "tensor = scipy.io.loadmat('../NYC-data-set/tensor.mat')\n", "dense_tensor = tensor['tensor']\n", "rm_tensor = scipy.io.loadmat('../NYC-data-set/rm_tensor.mat')\n", "rm_tensor = rm_tensor['rm_tensor']\n", "nm_tensor = scipy.io.loadmat('../NYC-data-set/nm_tensor.mat')\n", "nm_tensor = nm_tensor['nm_tensor']\n", "\n", "missing_rate = 0.1\n", "\n", "# =============================================================================\n", "### Non-random missing (NM) scenario\n", "### Set the NM scenario by:\n", "binary_tensor = np.zeros(dense_tensor.shape)\n", "for i1 in range(dense_tensor.shape[0]):\n", " for i2 in range(dense_tensor.shape[1]):\n", " for i3 in range(61):\n", " binary_tensor[i1, i2, i3 * 24 : (i3 + 1) * 24] = np.round(nm_tensor[i1, i2, i3] \n", " + 0.5 - missing_rate)\n", "# =============================================================================\n", "\n", "sparse_tensor = np.multiply(dense_tensor, binary_tensor)" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Imputation MAPE: 0.537513\n", "Imputation RMSE: 4.8318\n", "\n", "Time step: 40\n", "Time step: 80\n", "Time step: 120\n", "Time step: 160\n", "Final MAPE: 0.58577\n", "Final RMSE: 5.26301\n", "\n", "Running time: 1084 seconds\n" ] } ], "source": [ "import time\n", "start = time.time()\n", "pred_time_steps = 24 * 7\n", "rank = 30\n", "time_lags = np.array([1, 2, 24])\n", "maxiter = np.array([200, 100, 200, 100])\n", "small_sparse_tensor = sparse_tensor[:, :, dense_tensor.shape[2] - pred_time_steps : dense_tensor.shape[2]]\n", "small_dense_tensor = dense_tensor[:, :, dense_tensor.shape[2] - pred_time_steps : dense_tensor.shape[2]]\n", "tensor_hat = st_prediction(dense_tensor, sparse_tensor, pred_time_steps, rank, time_lags, maxiter)\n", "end = time.time()\n", "print('Running time: %d seconds'%(end - start))" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [], "source": [ "import scipy.io\n", "\n", "tensor = scipy.io.loadmat('../NYC-data-set/tensor.mat')\n", "dense_tensor = tensor['tensor']\n", "rm_tensor = scipy.io.loadmat('../NYC-data-set/rm_tensor.mat')\n", "rm_tensor = rm_tensor['rm_tensor']\n", "nm_tensor = scipy.io.loadmat('../NYC-data-set/nm_tensor.mat')\n", "nm_tensor = nm_tensor['nm_tensor']\n", "\n", "missing_rate = 0.3\n", "\n", "# =============================================================================\n", "### Non-random missing (NM) scenario\n", "### Set the NM scenario by:\n", "binary_tensor = np.zeros(dense_tensor.shape)\n", "for i1 in range(dense_tensor.shape[0]):\n", " for i2 in range(dense_tensor.shape[1]):\n", " for i3 in range(61):\n", " binary_tensor[i1, i2, i3 * 24 : (i3 + 1) * 24] = np.round(nm_tensor[i1, i2, i3] \n", " + 0.5 - missing_rate)\n", "# =============================================================================\n", "\n", "sparse_tensor = np.multiply(dense_tensor, binary_tensor)" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Imputation MAPE: 0.526796\n", "Imputation RMSE: 4.89139\n", "\n", "Time step: 40\n", "Time step: 80\n", "Time step: 120\n", "Time step: 160\n", "Final MAPE: 0.57862\n", "Final RMSE: 5.31899\n", "\n", "Running time: 1084 seconds\n" ] } ], "source": [ "import time\n", "start = time.time()\n", "pred_time_steps = 24 * 7\n", "rank = 30\n", "time_lags = np.array([1, 2, 24])\n", "maxiter = np.array([200, 100, 200, 100])\n", "small_sparse_tensor = sparse_tensor[:, :, dense_tensor.shape[2] - pred_time_steps : dense_tensor.shape[2]]\n", "small_dense_tensor = dense_tensor[:, :, dense_tensor.shape[2] - pred_time_steps : dense_tensor.shape[2]]\n", "tensor_hat = st_prediction(dense_tensor, sparse_tensor, pred_time_steps, rank, time_lags, maxiter)\n", "end = time.time()\n", "print('Running time: %d seconds'%(end - start))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "**Experiment results** of short-term traffic prediction with missing values using Bayesian temporal tensor factorization (BTTF):\n", "\n", "| scenario |back_steps|rank|time_lags| maxiter | mape | rmse |\n", "|:----------|-----:|-----:|---------:|---------:|-----------:|----------:|\n", "|**Original data**| - | 30 | (1,2,24) | (1000,500,500,500) | **0.5876** | **5.30**|\n", "|**10%, RM**| - | 30 | (1,2,24) | (1000,500,500,500) | **0.5866** | **5.30**|\n", "|**30%, RM**| - | 30 | (1,2,24) | (1000,500,500,500) | **0.5876** | **5.40**|\n", "|**10%, NM**| - | 30 | (1,2,24) | (1000,500,100,100) | **0.5858** | **5.26**|\n", "|**30%, NM**| - | 30 | (1,2,24) | (1000,500,100,100) | **0.5786** | **5.32**|\n" ] }, { "cell_type": "code", "execution_count": 70, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAT8AAAECCAYAAACIQtDYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAMTQAADE0B0s6tTgAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJztnXlYVdX6x79nEFCRABGDSwJOoIIDoSim8RPTHEovWWmpJVomal0Tr3od0hJJrOsszqillpJTeTFRIQsjNEWtnBBnM0BmOXg47P37gzgK7LPYZ3MmOO/nec7zxF57rfXudY5va3j395XxPM+DIAjCypCb2wCCIAhzQM6PIAirhJwfQRBWCTk/giCsEnJ+BEFYJeT8CIKwSsj5EQRhlSjNbQBBEA0LlUqFuXPnYubMmXB1dcXOnTuRkpKCJk2aAABCQ0Px4osvIicnB6tWrUJBQQHc3Nzw/vvvo3Hjxiazk5wfQRAG4+rVq9iwYQPu3bunvZaRkYHIyEh4e3tXuXfTpk3o378/+vTpg/j4eOzZswdjx441ma207CUIwmAkJiYiPDwczs7OAACe53H9+nXs3r0bkZGRiIuLQ1lZGTQaDS5evIhevXoBAEJCQpCammpSWy3G+alUKkyfPh1ZWVnmNoUgCIlERESgQ4cO2r+LiorQvn17jBkzBkuWLMHDhw+xd+9eFBUVoXHjxlAqKxafTk5OyMvLM6mtFrHsFZoqEwRhehbIZKLuG372LLp27VrrfQ4ODpg9e7b276FDhyI2NhYvvPACZNX6kstNOxezCOdXOVVevXo1877yLz8QvK4YvQIozRGuZOcClPxZ83oTN0ClY5bZ2BUofaCjvea1lDHs0LfMzkW6HcV3hMvsPYDCTOEyh9b691frM0sdR4llQt9pbd+nqX4H5hgPPRHrfsQ4PgC4f/8+Ll++jOeffx4AwHEc5HI5HBwcUFJSgvLycigUCuTl5cHJyUlve+uCRSx7q0+VCYIwDzKRH7EolUps374dOTk54HkeCQkJ6NGjB5RKJXx9fZGSkgIASE5OFu1QDYVFzPwIgrAMDD0bcnFxQXh4OKKjo6HRaODj44OXXnoJADBhwgSsWbMG+/btg4uLCz74QHhlZyzI+RFEAyU9PV3v2ZShnN+aNWu0/927d2/07t27xj0tWrTAggULDNSj/pDzI4gGipRlpMIIdlgq5PwIgtCiz35efUdGMvYEQVTyuchQl+kNwG3Ur5kf40ifv31UsEj2TH+gIKNmwVNtJYcP8PdPCvf1dDA7bEIo5AaoCLsxdDiOvn1VlkkKdTFx+IaUEBPGM/P5VwSLZI7t6xDWZMKxYv1G9MQiwj9MRP1yfgRBGBVrWvaS8yMIQos1OQRrelaCIGqBZn4EQVgl1rTnZxGnvTt37sSpU6cAVAgdDh061MwWEYR1sk7kae975ncbdcbsM78zZ87g8uXL+Oyzz6DRaPDhhx8iICAA7u7uNW9mnXgxXtbnLu+qcVnuM0ryCRqf94dgkcypI/sEsOiWcFmzVrpPKYtu6qjjKf0EUPJpr4W8rG9g8QXm6b3B7TDCM7N+V3pidodgQsz+rAEBAejcuTMUCgVyc3PBcRxsbW3NbRZBWCW052dilEolvvrqK3z33Xfo1auXVgWWIAjTYk17fhbzrCNHjsTmzZuRm5uLY8eOmdscgqj3pKen611HLvLTEDD7zO/27dvgeR6tWrWCra0tunfvjps3dexzEQQhGhI2YGN2J3737l1s3LgRGo0GZWVlSEtLI2FTgjAT1jTzs6hQF7lcjl69emHEiBHmNokgrJJdIkNdRpnfbdQZi3B+ojFkrgi75ijfNkWwiuKt1YZ9oZ1lRy02GiVfCFMsQcrL+lLHwxj1LDgcx9A5XrRlhsvh8bVI5/d6PXIbujD7nh9BEJaDNe35kfMjCEJLQ9nPEwM5P4IgtFCQM0EQVgnN/AiCsEqsac/PIk57ly9fjuvXr8PGxgYA8Oqrr6JHjx5mtoogrI9DCnFzvyHlnJEtMT4WMfPLzMzE4sWLYW9vz77R0GEHxXeE69h7gDu7XrBI3m0iO7/HnymCRTK33uxwhYf3al5v6i5ZyYZ/cEHYjub+4P/6RbisZRD4O8eFyzz66baRlR8j55xwey5dmDYy1WwYZUL5OGSO7SXnQmHayBh/nXljhK7/XcbMQyNV3UdPZCJDXRoCZnd+xcXFKCwsxIoVK5CXl4egoCCMGDHCqr4EgrAU5HLr+XdndueXn58PPz8/TJgwAY0bN0ZMTAySkpLQr18/c5tGEFaHXOSyl4VKpcLcuXMxc+ZMuLq64qeffsL+/fsBAC1btsSkSZNgb2+P9PR0rFq1Sqvi5O3tjYiIiDr3LxazOz8PDw9ERkZq/x44cCBOnDhBzo8g6kh6erre4gZ1XXFdvXoVGzZswL17FVskubm5+PLLLxETEwMHBwfs2rULe/bswbhx45CRkYF//vOfZlNuN/vJ9rVr13D69Gnt3xzHQS43u1kEUe+Rouoik8tEfXSRmJiI8PBw7WxOJpPh3XffhYODA4CK2V1OTsXe97Vr13DmzBlERkYiJiYGDx7o2Ls0Emb3MhzHIS4uDiUlJdBoNEhMTKSTXoIwEzKZTNRHFxEREVVUmZycnBAQEAAAePToEfbt24fAwEAAQNOmTTF48GB89tlnCAgIwMqVK437cNWwiFCXQ4cO4ejRoygvL0dQUBDefPNNc5tEEFZJkmNTUfc5JacwZ5aTJ0/GRx99BFdXVwBAUVERli5dCnd3d7z33nuCdcaNG4c1a9agSZMm+hsuAbPv+QHAkCFDMGTIkNpvZCWhYSQV4vMv17zu6CM5dCYj2EewqO3Jy7UkmmGFK+hQI5GabIgRxlMvkvIYWMFHMEwHAJq61/8ERoxn0xexp736LKmzs7MRFRWFwMBA7cRGo9Hg4MGDCAsLAwDwPG/yLS+zL3sJgrAc6rrsrU5ZWRmioqLwwgsvYPTo0dq6SqUSP/74I9LS0gAAycnJaN++Pezs7IzyXEJYxMyPIAjLgHWYIYUffvgB9+/fR3JyMpKTkwE8Dmn54IMPsHHjRuzatQuOjo6YPHmyQfuuDXJ+BEFoMUScHwCsWbMGANC/f3/0799f8B4vLy9ERUUZpD8pkPMjCEKLNb1ZRc6PIAgthl72WjJmD3U5cuQIEhMTtX/n5OQgICAAU6dONaNVBGGdpHqIE0Poece0AcnGwOzO70nu3buHqKgofPLJJ9oI8SoYOvyBofTBnV4rWCQPjGCGn3BHPhGuN2AeWz1EINyiItRC9zNzmQeF+2r9MvPZhJRPgAr1E6GwoIoyH91KJSx1GZaqS+5vwmXOfuCzTgmXuXZn1xNQpZF59JP82+Gzzwr31aIbUwFH6LuRt36ZGYLECrlhKfEYUtXllGcLUfd1v5mtd9uWhkUtezdv3ozXXntN2PERBGF0rGnZazHO748//kBBQQH69u1rblMIwmqxpgMPiwlyPnLkCIYOHWpVg08QxiQ9PV3vOjK5uE9DwCIeQ6PR4MKFCwgKCjK3KQTRYJCi6iJXyEV9GgIWsey9desW3Nzc0LhxY3ObQhBWjdyKVl4Wcdp78uRJpKWl4V//+pe5TSEIq+ZCx3+Ius//j7tGtsT4WITzE41UFRChsI9aEtcYQzHl/ojegkVPx6cwVEAkhvcw1WCkhgzpUJ4xiqqLodVUGM8sVTnH0HYYYzz05LdOHqLu8/tdx7+BeoRFLHsJgrAMKNSFIAirpKEcZoiBnB9BEFqsKdSMnB9BEFpkVpQ8zCKc3/79+5GUlIRGjRohODhYK21NEIRpsaY9P7Of9l64cAFbt27FJ598Ajs7O8TExOD//u//KOCZIMzA1Z7tRN3XLvWqkS0xPmaf+V2/fh1du3bVZmzq1q0bTp06Jez8mGEChg51kRiSIDEMhv8zpcZlmVtvdtKmrNOCRTLXQHZSG1YZMwHTrZrXm7ViJ2ZijQcrxIShSsOsp0t5RnI4C+N3wBpHXWPFskPyd8awUU+sadlr9if19vbGuXPnUFxcDLVajdOnTyM/P9/cZhGEVWLoBEaWjNlnfv7+/ggJCcGCBQtgb28Pf39/XL1a/6fUBGFu0tPT9X+/14r2/Mzu/FQqFXr06IGhQ4cCqDj8qEx0TBCEdKQJGyiMYIllYvZlb3Z2NmJiYqDRaFBcXIykpCQEBweb2yyCsEpkcpmoT0PA7Ke9ALBv3z6cOHECHMdhyJAhGDBggLlNIgir5Eb/zqLu8zp63siWGB+LcH6iYZ14sU4HhU4jm3myT/lY7Uk9eWO1qeNFeFa+DYOefrPs+NuW+i1swKgj9fRb34gAYz0zy349ufmCuKWyZ6L+QqmWhtn3/AiCsBzkSrPvhJkMcn4EQTymoWjUi4CcH0EQWgxxmKFSqTB37lzMnDkTrq6uOH/+PLZv3w61Wo1evXph5MiRkMlkuHHjBtatWweVSgVfX1+88847UCpN55Ksx80TBFErdT3tvXr1KubPn4979yr2IdVqNWJjYxEZGYlly5YhMzMTv/76KwBg1apVePvtt7FixQoAQGJiovEf8AnI+REEoUWmUIj66CIxMRHh4eHa3NsZGRl4+umn8fTTT0OhUKBPnz5ITU1FdnY2Hj16BF9fXwBASEgIUlNTTfKMlZh82Vt9SgwA5eXlWLx4McLCwtCpUyfdlVmnV40ZgdHNPIWvs2S+We2x7GC9T8lqU4ctMsf20voysB3MNpl1LKSMVYf5fUoYD6l2SC2TcKqri7oueyMiIqr8nZubCycnJ+3fjo6OyM/PR15eXpXrTk5OJn+t1aTO7+rVq9iwYYN2SgwAd+/exfr163Ht2rXaG5Aa2iEhnwIzxITxcjo7NEV3uAJ3em2Ny/LACLaNt48K9/VMf/1Df4CK/0lIEA3g//pF2I6WQeDz/hAuc+oI/ub3wmWeA5mCCKw2udSVNa7Le74vPfSHVcayUUB0QuYayP59sNrLOSdcz6UL+zvTE7Hv7Yp9dY7n+RptymQycBxX5brQfcbGpMve6lNiADh27BheeukltGsnTkqHIAjjIZPLRX3EvjrXvHnzKjO6/Px8ODk5oXnz5sjLy6tx3ZSY1PlFRESgQ4cOVa6NHTsW3bt3N6UZBEHoQKaQi/qIpW3btrh79y7u3bsHjuPw448/olu3bmjRogVsbGxw8eJFAEBycjK6detmrMcShEJdCKKBIkXVxdDv7drY2GDy5MlYtmwZ1Go1unXrhp49ewIApk6divXr10OlUsHb2xuDBg0yaN+1Qc6PIBooUlRdDCVmumbNGu1/+/v7Y+nSpTXu8fLyQnR0tEH6kwI5P4IgtDQUoVIx1C/nJzW0Q2fYge46zBATxikaOzRFd7iCPDBCuIBl4zP9dfclJfQHYJ8QPtVW2I6WuvOtyJw66i7zHKi7L3sPSW3Ke74vXCA19IdVxrLRNVD4Ouv3wWrPpYvuehJOdXX2o7QePb/6perCCgVghRDoVHVhhD+wwkEKM4XLHFoL524A/s51wSrTYSMrB4aeoRbA3/8oWfbrm3OjltwTfO5vwnY4+wmHzgAVTpYxVnz+ZeE2HX2kjSMjnIU5jqwwGKExljK+wN95Rhi/VVabepI3nvE/pCdw2iwcplSfqF8zP4IgjEpDESoVAzk/giAeQ3t+BEFYIzIFOT+CIKwRmvkRBGGNWJHvM/1pb3VVl//9739aHa+AgACMHj3aqmKNCMKSKJw6VNR9Dqu+M7Ilxsesqi537tzB999/jyVLlsDGxgbz58/H+fPn0aWLcEwTM+xA3xATXeElANDME9y1/YJF8jbDmeEbkhMO6QrRYISlSLWDzz4rXNaiGzvER9c4MkItuCtfCxbJ278O7uIXwmUdxoDL2Ctc1jYM/N1kwTLZP0IElW5kz/SvJayG8Ts4UfPNBACQ953B/s3pCnVh2SE1HIel+KIn1nTaa1ZVFw8PD3z++eews7PDw4cPoVKp0LRpU1OaRBDEk8hFfszI999/j+nTp2P8+PHIycnB0qVLUVJSonc7Zld1USqVOHLkCKZMmQJHR0d4eXmZ0iSCIJ5AJpOJ+piLffv24ejRo3jllVfAcRyaNGmC0tJSbN68We+2LELGfsCAAdiyZQucnJywe/duc5tDEA2C9HQJuXVlMnEfM3Hs2DHMnDkTwcHBAIAmTZpg2rRpkp7VrM4vKysLV65U7E0pFAoEBwfj1i0d+ygEQeiFJFUXhUzUx1yUlpbC0dGxyjVbW1tJs1GzOr+ioiKsWrUKJSUl4DgOJ0+erLEsJgjChMhEfsxEhw4d8OWXX4LjOO21ffv2wcfHR++2zCJsMHnyZHz00UdwdXVFQkICjhw5AoVCgY4dO2Ls2LEmzd1JEMRjSmb9U9R9TT7dZ2RLhMnNzcWSJUvw559/Qq1Wo1mzZnjqqacwa9YsuLgwlHsEqF+qLiwVDUZyIJ0JjBghDvyNQ4JFMq8h7KRCLBUTlo1Cz9bY1SiJd1hhMNy5jYJl8i7v6B5HVkInZlgNQ7mFNY6s59aRZIn9+2B8n1JVXYSUbpq6M0OomGEwUn/7elIyW6TzizaP8wMAjuOQmZmJ7OxsODs7o23btlAw0mnqQvSy98yZM/j0008xc+ZM5OfnY+vWrdBoNHp3SBCE5VLXpOWm4OHDh1AqlWjZsiUaNWqEmzdvIjNTRzwsA1Hry+PHjyM+Ph4DBgzA/v37IZfL8fvvv2P79u0IDw/Xu1OCICwUCw9yPnjwIHbu3AmhBevXXwsH1OtClPM7ePAgZs2ahVatWuHAgQNwcHDArFmzMHv2bHJ+BNGAsPQ3Sw8cOIB///vf6NatW53jDUU5v4KCAnh4VFWFdXJyomUvQTQ0LHzmp1Qq0aVLF4MEWova82vTpg0OHap6AHD8+HF4e3vr3aFKpcL06dORlVWxgbtz505MnjwZM2bMwIwZM3D48GG92yQIwjBYeIwzXn75ZWzZsgVZWVkoLi6u8tEXUae9d+7cwaJFi2BnZ4e//voL3t7eePDgAebMmYNWrVqJ7qxS2ODOnTtYsWIFXF1d8fHHH2PMmDGSHClBEIZF/cmrou6zmbfHyJYIc+zYMWzatKlKnF8lRtnz8/DwwPLly3HmzBnk5OTA2dkZAQEBaNKkiV6dVQobrF69GgDA8zyuX7+O3bt3Izs7G506dcLo0aPRqFEj4QZYR/qMJDqSwkgkhpjweX8IFsmcOjLtFwoJkbXoxu6LoW7CDD9hKc+wQiokqLowy4yhcJJ1qsZlmWv3WkJFJH7XLFUdgd9BxW+A0ZfU74z129cXC1/2fvXVVxg/fjz8/f0hr2OOYdHRxEqlEj4+PmjfviL1XklJCUpKSvQKLIyIqJqesaioCO3bt8eYMWPQsmVLxMbGYu/evXj99ddFt0kQhAGx8BMPjuMQGhpqkD0/Uc7vp59+wqZNm6BSqWqU6TvVfBIHBwfMnj1b+/fQoUMRGxtLzo8gzERdfMqRI0e0wsQAkJOTg4CAAPj6+iI+Ph4ODg4AKkSLR40aJamPgQMH4ptvvsGwYcN0rxBFIsr57dmzB2FhYejbt6+kSGpd3L9/H5cvX8bzzz8PoMKr13UqSxBEBenp6XqLG9RFtGDAgAEYMGAAAODevXuIiorCm2++ia+//hrjx49Hjx49JLddyS+//II7d+4gPj4eNjY2VWaA27Zt06stUc4vPz8fQ4cONbhjUiqV2L59Ozp16oTmzZsjISHBIANEEIQ0VRdDLXs3b96M1157Dc7Ozrh27Rry8vKwe/dueHl5Ydy4cZJFi8ePH28Q+wCRzs/f3x/p6ekICAgwWMcA4OLigvDwcERHR0Oj0cDHxwcvvfSSQfsgCEIPDOD7/vjjDxQUFKBv377gOA7NmzfHK6+8gnbt2mHXrl2Ii4vDlClTJLXdsWPHuhv4N6JCXVauXImff/4ZrVu31q7bK5k5c6bBjCEIwryUf/6GqPsuhP5b58xy+fLl6Nq1K0JCQmqUPXz4EFOmTEFcXJwk+yZPnqzzsKMyikQsomZ+bm5uCAsL06tho8AKE2Ad9+tSI5Go9MGsxwrRYIW6CIRNyJz9mOEPzGRDLPsZSZH0Dvto7Co5LEWyYg1rTITKmrjV8vuQaIeUpFlSf3MsO1jjry8KcVtbuhyfRqPBhQsXMHHiRABAYWEhUlJSMGjQIAAV+/p1OTcYM2ZMlb+Liopw7Ngx9OnTR++2RDm/V199HPhYUFAAe3t7gx58EARhIdQxzu/WrVtwc3ND48aNAQB2dnaIj4+Hj48PWrduXed9/Z49e9a49uyzzyIqKgpDhgzRqy1Rzq+srAzbtm3DDz/8ALVaDaVSid69e2PChAmwsbHRq0OCICwYWd0ONe/fv18l9tfGxgbTpk3DunXroFar4e7uLnm/TxdNmzZFbm6u3vVEOb9du3bhxo0bmDNnDlxdXXH//n3s2LEDO3fuxNtvv613pwRBWCh1PO0NDg7WJheqxM/PDzExMXVqt5LvvquaLF2j0SAtLQ1t27bVuy1Rzi81NRVRUVFwcnICADg7O2P69OmYOXMmOT+CaEhY+HbWr7/+WuVvuVyO1q1bY/jw4Xq3Jcr5PXr0CPb29lWu2dvbC75cXBsqlQpz587FzJkzcevWrSpviOTn58PNzQ0ff/yxcGU7xqt0rPcYdcl5s2S+WX2x6rE2mRn1ZM5+wgVN3HS318yT0RfDfofW0uo1dhW+znpmXXWk9gWwx0RXGfP3IdGOZgxRD11lUn9zLDukHGzowsLf7f3oo48M1pYo5+fr64svvvhCm1xIo9Hgiy++0DtjUqWqy717FSdvgYGBCAwMBFBxKjRnzhx2EKPUE1hdp72ME8DybcL7Eoq3VjPt4M6uFyySd5vIPu29+X2NyzLPgcw63PFo4b76zQb/4IJgmay5P+79s5dgmfu+n8F9PV24zdc/Fz5dbubJFCgojxX+PhWTNoP76b/CfT33Ibgrwq9Nytu/rjPJTpNP9+Fi15r/Q+iQflPyb4f7fauwHZ3eBn/vhGCZzL0vyteOq3FdERGHgkmDBes8Ffs/cBeEE2/L/ccL/j6Av38jhhQ2sNB3e6svd4UYOnSoXm2Kcn5vvfUWPvnkEyQnJ8PJyQl5eXlo0aIFZs2apVdn1VVdnmTnzp14/vnn4enJmM0QBGFcLPT10urLXSEM6vxiY2Px4osvwtvbG8uWLcOlS5dQUFAAFxcXtGnTBkuXLq0iTFAb1VVdKsnKysLZs2excuVKvYwnCMLAWOienyGXu5UwnV9ycjJ+/vlnTJ48GUFBQfDzq7ovdenSJYMYkZiYiNDQUNja2hqkPYIgpAkbWOqy90lOnjyJpKQk5OTkwNHRESEhIVpxFH1gznFtbW0RHh6OlStXYu/evZKNrY20tDQ899xzRmufIKwRycIGFqxjf/ToUWzZsgW+vr4ICwtDx44d8cUXX+D774X3RFkwZ34ymQwhISFwcXHBsmXLcPfuXUyaNAlKpWgN1FopLCxEaWkp3N0lbM4SBGFYLHTPr5JDhw5h1qxZVeL6unbtilWrVmHgwIF6tcUUNnjrrbe0GllZWVmIiYlB48aNMWPGDDg4OFQp14fJkyfjo48+gqurKzIyMhAXF4eoqCi92yEIwrBwO/8l6j75G8uNbIkwb731FuLi4qrI63Ech/DwcGzdulWvtkRP4VxdXbFo0SKsXLkSs2fPrpOay5o1a7T/3bZtW/GOj/VCu76hDHbNme2xQi30zcUB/J2Pg1Xv/smadZ4OZodhsEIjGPklyuOED54U49aCOxMr3GbAJN3CBgyhBO7iF8LtdRgDLm2VcFmPqeAyhLdZ5G3DwF2NFy5rNwLczytqXu/1QS0CC/oJTgB/5+lg5FDhMg/WtKP1y4jXsWQcwfOC+UeAihwkzN+VIYUN6vh6m7Hx8PDAyZMnq2yTpaSk4B//+IfebTGdX/VJoZ2dHf79739j165dmD9/PtRqtd4dEgRhwVh4kPOoUaMQHR2N5ORktGjRAllZWbhy5YpeUSeVMJ3ff/7zH50GeHp6VtHrJwiiAWChp72bN2/GwIEDte8Jp6SkID8/H506dcLEiRPh6sp4A0YHTOfn6+urs0zoBWaCIOo5Fhrnl52djRkzZsDX1xcDBw7EiBEjTJe6kiAIK8BCZ36zZs1Cbm4ujh8/ji+//BJbt25F//790b9/fzg6Okpqk5wfQRCPseBQF2dnZ4wYMQKvvPIKzp07h+PHj2Pq1KkICAjAwIED9c7vISqHhyF5UtXF1dUVycnJOHDgAORyOfz8/DB27FhSiSYIM8EdmCPqPvkwywhNKywsxNatW5GSkqJ3DnGTzvyqq7rcu3cPu3btQnR0NJydnbFp0yYkJCTofkFZsqqLQL1a8iJw5zcJFsk7T2DawQpXYNoopMyhK//I33X4rNM6+goEn/eHcJlTR3CXdwmWyX1GMcNIdIa6sEKGWKooAuE9QEWIDzPsI/+KcJlje8H+5J3elvzbYYWzMMNgcs7VvO7SBdzvwol75J3GMVViWMo5zN+VvtSTiUdWVhaSk5Pxww8/QKlUYvTo0Xq3YVLnV13V5ebNm/Dx8YGzszOAikzuBw4c0FudgSAIA2HBoS5qtRqpqalISkrC5cuXERAQgIkTJ6Jz586S2jOp86uu6uLp6Ynt27cjJycHzs7OSE1NRX5+vilNIgjiSSz0wGPdunVITU2FnZ0d+vXrh6lTp2onTVIx64GHu7s73njjDcTExMDGxga9evXCtWvXzGkSQTQYJKm6WOiBR3Z2NiZNmoTu3bvXOcSlErM6P7VajbZt22qTm5w8eVJSsCJBEDWRrOpigcybN8/gbZrVzavVaixcuBAlJSUoKytDQkICevUSllgnCMIEyBXiPg0Ak4e6AFVVXZKTk3Hw4EGUlZWhT58+eO2110xtDkEQf8MdWyzqPnmo8Kuv9QmzOD/JSFTm0KnqwgpZYYQ4MEMjBEIcgIowB2Z/ulRdGM/M3zgk3JfXEPD5l4XLHH2YaipM+3WF4zDCMLhr+4X7ajOcGVbD//WLsB0tg5hlQmow8rZh7FAX1hiz1FQYzy303ci8hrDDY1jhPYx6TLUjPdGmUCppAAAYoklEQVSVFKs68n76CwlYGvSGB0EQj7HQPT9jQM6PIIjH1FHPb/ny5bh+/TpsbGwAAK+++ipcXV2xbt06qFQq+Pr64p133jGoGrxUzG8BQRCWQx0nfpmZmVi8eDHs7e2116ZPn4533nkHvr6+iI2NRWJiIgYNGlRHQ+uOZQb1EARhHuqQwKi4uBiFhYVYsWIFIiMjsWfPHmRnZ+PRo0daebyQkBCkpqaa8ol0QjM/giAeU4c9v/z8fPj5+WHChAlo3LgxYmJioFAo4OTkpL3HycnJYt7iMulp73fffYfjx49DJpOhTZs2ePfdd7Vr/zVr1qBTp04ICQkxlTkEQVSD++m/ou47b9+v1iDqtLQ0JCQkoLy8HB9//DEA4M8//8SSJUuwfLl5EiA9iclmfhkZGUhKSsLixYtha2uL1atX4/DhwwgODsamTZtw/vx5dOrUid0IM9SFpdqhv6oLS3FE35CVx/UYyYgElFbkPqPYyjN7ZwkWycM+BYpuCtdr5gnu0g7her5vMsNnUHRLoL1W7ARGAol8gIpkPty5jcJlXd5hKpww2xRIPFWRdEpaqAtrHFlhQbrUZaQmnWKFLjGfTV9ETvyEHN+1a9eQl5eHwMBAABVZ1QAgLy9Pe09+fn6VmaA5MdmeX9OmTTF+/HjY2dlBJpPB09MTOTk5OHHiBJ599ll6s4MgLIE67PlxHIe4uDiUlJRAo9EgMTERoaGhsLGxwcWLFwEAycnJ6NatmymfSCcmm/m5ubnBza0i6LKgoADff/89Jk2aBD8/PwDApUuXTGUKQRC6qMOeX7t27TB48GDMmTMH5eXlCAoKwnPPPQcPDw+sX78eKpUK3t7eFnHSC5jhwCMrKwuffvopQkNDtY6PIAjDI0nVpY5xfkOGDMGQIUOqXPPy8kJ0tLg3R0yJSZ3fjRs3EB0djeHDh1uM9yeIhkpDUnUxBiZzfoWFhYiKisKECRMQFBRkqm4JgtAH6/F9pgt12bVrFw4dOqTd9wMqZOtHjRoFgEJdCMIS4E6tFnWfvPsUI1tifOqXqgtLuYUVriAhORBTuYWhKsJsk6W+oVN5RpqNzPAHofEApI1JXWw0pEpPLTaWrx0nWEUREQcU3xFuz97DsPbX6fuUWKYn3Om1ou6TB0bUfpOFQ294EATxGNrzIwjCKrEe30fOjyCIJ6CZH0EQVkkd4/zqEyZ1fkLCBmlpafjmm2/A8zzatGmDiRMnWoTQIUFYJVY08zPZaW9GRgZiY2MRFRWlFTbw9vbGwYMHERMTA0dHRyxbtgz+/v7o37+/KUwiCKIa3PlNou6Td55gZEuMj8mmWE8KGwDQChusXbsWSqUSpaWlKCoqQtOmTXU3wgo/eXBBsEjW3F9YdcShtfRQl/wrwn05tme3qW+ISV3CQSQr4Jgw1MVUZbWE/rBUbgwaaqRLYUhbZoSQIUInJlvgu7m5oWPHjgAeCxsEBgZCqVTi9OnTiIiIQGFhIbp06WIqkwiCqI5cLu7TADD5U2RlZWHhwoVVhA0CAwOxZcsWBAQEYONGYY03giBMQB0kreobJnV+N27cwLx58/DCCy8gLCwMhYWFuHDh8XK1T58+uHVLQDCTIAi9SU9Pl1BLJvJT/zGZ86sUNggPD9cqumg0GqxcuRIPHlTsZ6SkpKBDhw6mMokgGjSSVV2sZOZndmEDT09PfPPNN5DL5XjmmWcwYcIENGnSxBQmEQRRDe7il6Luk3cYbWRLjE/DETYw8CmlUA4GoCIPA/Mk+Ob3gkUyz4FsGwsyal5/qi27r7w/hPty6ijc3t9tln/5gWCRYvQK9nPrEohgCAOwxoPL2CvcV9swZj3+zxThMrfe4O8m17z+jxDpJ/vZZ4X7atGNeeovJH4haxkE7tp+wTryNsOZghnM79qAOTx0nXxXR+77pt5tWxoUTUwQxGPoDQ+CIKyThrGfJwZyfgRBPEamMLcFJoOcH0EQj2kgJ7liIOdHEMRjrMj5mfS0V0jVJSkpCfHx8XBwcABQNa8HQRCmRdcJfHXkbcOMbInxMdnMLyMjA0lJSVi8eLFW1eXw4cO4ffs2xo8fjx49etTeiKFzeDDqMIUSWKEd908K13s6mNmfUAiEvM1wto33Tgj35d5XujADI8+IUH8y977ssSq6KdxeM89a6ul406dZK/C5vwnXc/YTHH/Z08HM0B/mGDP6Kt8mnMRH8dZq8DcO1azjNYQpbiEUpgNUhOoItadtkzFWeiOnPT+Do0vV5dq1a8jLy8Pu3bvh5eWFcePGsZVdCIIwItaz7DWrqsuzzz6L5s2bY8SIEVi6dCmcnZ0RFxdnKpMIgqiOTC7u0wAw+YFHVlYWPv30U4SGhsLf3x/+/v7asmHDhmHKlPqfD5QgLIH09HS93++V1fHAoz7t65vU+d24cQPR0dEYPnw4Bg0ahMLCQqSkpGiFDjiOg0JhPXsOBGFMJAkb1GHZa5B9fRNiMudXqeoyYcIEBAUFAQDs7OwQHx8PHx8ftG7dGgkJCRY3QARhVdQhyLm+7eubXdXF398f27dvh1qthru7O6ZMmUKqLgRhJvhbR0TdJ2s1gFleUFCA//znP3jvvffw3Xff4ZVXXkG7du2wa9cu5ObmWsT2Vj1TdWHlP9AzN0It4THcla8Fi+TtX2fXO71WuF5gBDP8pDx2fI3Likmb2cozFzYL9+U/nh0qknVKsEjm2h18/mXhMkcfwTKZow8zZIWl3MKdXS9c1m0iu80TS4Xr9Z2B+yN617j+dHwKM4SHGerCUFrhLu8StsNnlM6QG5YiDfe78GGfvNM4cMejhcv6zTaoqgt/+6io+849cNG5rK7c13/uuecQFlY1HvDhw4eYMmWKRRxsNoxjG4IgDIQ4JWddjk9IrT0hIUFbbkn7+uT8CIJ4jFwh7iOAkFp75b5+ZmZFBkVL2tend3sJgngC6ae9hw4dgkqlQnx8POLj4wFU7OtPmzYN69atq7KvbwmQ8yMI4jF1iPMbNWqUzvi9mJgYye0aC3J+BEE8gfXshJlV1WXAgAFYv/7xiV9xcTEAIDY21lQmEQTxBLqEOaojezrYyJYYH7Oquly6dAlLl1aELajVasyZMwdvvslIjKJvOAvATmDEUthgJoxh9MVQMWHW06U8w1Rn0R2Wwgx1YSTl4Y58IlgmHzBPOFyktlARhhoJM1SEMf5M+wXalPuMYqv+SE0Sxfj9oDCz5nWH1vjd/xnBKp0u3GaqBTGVW1jPpjfWM/Mzu6pLJd9++y28vb0lvpJDEIRBsCIxU7OqugQGBgIAVCoVDh8+jJEjR5rKHIIgBBEX59cQMPkcNysrCwsXLkRoaCj8/PwAAD/++CO6du0KZ2dnU5tDEA2W9PR0/StZkaSVSZ+ievR3JadOncJzzz1nSlMIosEjaQtJJhP3aQCYVdUFAHieR0ZGBnx9fU1lCkEQOmkYszoxmF3VZfDgwYiMjMTGjRtNYQZBEAz4HHFLZZlL/T+YrGeqLoxQEX1VOxq7Musw1VlYqi4XvxCu12EMO6TizvEal2Ue/diqLpkHhftq/TJ7PFhJhRgqJkJhHzKnjuz2cs4Jt+fShZlUiKk8w0jcJGS/rGUQezxYSZsEvheg4rthKf9waatqXu8xlZk8ipX8iqX4wnw2PeEfnBd1n6x5Z73btjToDQ+CIJ6gYezniYGcH0EQj2kghxliIOdHEMQTkPMjCMIaqUMOj/oGOT+CIB5jRctes6q6vPvuu7hw4QJ27NgBAGjVqhXeffdd7fu/BEGYFl2n0dWRObY3siXGx2TOLyMjA7GxsYiKitKqunh7e2P//v2YP38+WrVqhf379yM3Nxfh4eHCjeir3MIqq60OU52FoS4jVdVFSLWjWSt2HUaoiGQbWQohOpRKmHawlEpYZcxQHYZiTe5vNS7LnP2kJ79ihcGwwlZ0hQVJ/c2xxsqAqi58/lVR98kc2+ndtqVhsnDuJ1VdZDIZPD098euvv6JFixZo1aoVACAwMBCnT582lUkEQVRDJpOL+jQEzKrqMmjQIOTk5ODGjRsAgJMnTyIvL89UJhEEUR16t9d4VOb0DA0NRY8ePWBra4sNGzaA53mEhoZCqaQzGIIwBOnp6RLEDRqGYxODST3NjRs3EB0djeHDh2PQoEHgOA7NmzfH4sWLAQBXr15Fy5YtTWkSQTRYpKm6NIwlrRhM9qRCOT0BYNGiRcjJyQHP8/juu+/Qq1cvU5lEEER1rEjPz+yqLr6+vtixYwfUajX8/f0xbtw4WvoShLnQdapcHXsP49phAuqXqgsr/EEoDAMAHFoLhiTIHNtLT4jECjtghtzo2V9t4SysBDqs8AdmiIyEkCGJyYGkJ6TSM2yllpAV1m9HcniVLiUhyeFaup+ZqZyjLw/viruv6T/0b9vCoCkWQRBPQAceBEFYI3Xcz/vpp5/wzTffoLy8HIMHD8aLL75oIMMMDzk/giAeUwfnl5ubi507d2LJkiVo1KgR5s2bh44dO2pfYrA0GsaxDUEQBkJ66srz58/Dz88PzZo1g52dHYKCgpCammoasyVQvw48CIIwLroOVqpj51Lj0r59+/Do0SNt/u1jx44hIyMDEydONKSFBqPezvwk5SSVWM+UfUmt11D7klqPbJSInYvgJ/3SHWz7+pD2I2RD9XkUz/OQWfCrcDTzIwjCICQnJ+PSpUt47733AADx8fEAgBEjRpjTLJ3U25kfQRCWRefOnXHhwgUUFBSgtLQUqamp0l6xMxE08yMIwmD89NNP2Lt3L8rLy9GvXz8MGzbM3CbphJwfQRBWSb2M81OpVJg7dy5mzpwJV1dxarU7d+7EqVMVibBDQ0MxdOhQUfWWL1+O69evw8bGBgDw6quvokePHjrvP3LkCBITE7V/5+TkICAgAFOnTq21r/379yMpKQmNGjVCcHAwwsLCdN4rNAbl5eVYvHgxwsLC0KlTJ1H1/ve//2ntDQgIwOjRo2tsUlevs3PnTqSkpKBJkyYAKsZTKJj1yXq3bt3C118/TvKdn58PNzc3fPzxx7XamJycjAMHDkAul8PPzw9jx46FQlE10Y5QioTKd8TXrFmDTp06ISQkpNY6aWlp+Oabb8DzPNq0aYOJEyfWeNdcqF5SUhLi4+Ph4OCgHctRo0Yx6w0YMADr16/XlhcXFwMAYmNjmX1R6gcDwdczrly5wkdGRvIjR47k//rrL1F1fv31V37+/Pm8RqPhS0tL+YiICP7u3bui6k6dOpUvKiqSZOvdu3f5iIgI/sGDB7Xee/78ef7DDz/kHz58yJeXl/PR0dF8amqq4L1CY3Dnzh1+3rx5/BtvvMH/9ttvourdvn2bf//993mVSsWXl5fzc+bM4dPT02vta+HChXxmZibzeVjfU0FBAT9lyhT+xo0btda7e/cu/+6772rHcOPGjfy3335bpc7Vq1f5Dz/8kFepVDzHcfzKlSv5b7/9ln/w4AG/ZMkS/s033+STkpJE1XnnnXf4vLw8nud5/r///S+fmJgoqt7atWv5X375Red46KpXyaNHj/jIyEj+7NmztdYZP348f/PmTZ7neX7fvn385s2bdfZL6KbeHXgkJiYiPDwczs7OousEBARg3rx5UCgUKCwsBMdxsLW1rbVecXExCgsLsWLFCkRGRmLPnj01jvNZbN68Ga+99pooW69fv46uXbuiSZMmkMvl6Natm3amWh2hMTh27BheeukltGunO7dC9XoeHh74/PPPYWdnh4cPH0KlUqFp06bMOjzP4/r169i9ezciIyMRFxeHsrIyUTZWsnPnTjz//PPw9PSstd7Nmzfh4+Oj/TsgIKDGuAilSMjJycGJEyfw7LPPCsqk6aqzdu1aODo6orS0FEVFRTXGQ1e9a9eu4ejRo4iMjMTq1avx8OFDUfUq+fbbb+Ht7V3lgIBSPxiXeuf8IiIi0KFDB73rKZVKfPXVV5g2bRr8/PxEOaT8/Hz4+flh8uTJiIqKwqVLl5CUlCSqvz/++AMFBQXo27evqPu9vb1x7tw5FBcXQ61W4/Tp08jPzxe8V2gMxo4di+7duzP7EKqnVCpx5MgRTJkyBY6OjvDy8mLWKSoqQvv27TFmzBgsWbIEDx8+xN69e0X1BVQoeZ89exYvvfSSKBs9PT1x9epV5OTkgOM4pKam1hgXoRQJgYGBGD58OEJDQwX70VVHqVTi9OnTiIiIQGFhIbp06VJrvWeffRbNmzfHiBEjsHTpUjg7OyMuLk5Uf0DFMv/w4cPa4GBWHUr9YDjqnfOrCyNHjsTmzZuRm5uLY8eO1Xq/h4cHIiMj4ejoCFtbWwwcOBBnzpwR1deRI0cwdOhQ0UGe/v7+CAkJwYIFC7B48WL4+vqaTNdwwIAB2LJlC5ycnLB7927mvQ4ODpg9ezbc3d2hUCgwdOhQ0WMCVMzsQkNDRc28AcDd3R1vvPEGYmJiMH/+fHh6euocl6ysLCxcuBChoaHw8/MT1b5QncDAQGzZsgUBAQHYuHFjrfX8/f0xe/ZstG/fHjKZDMOGDcOvv/4qur8ff/wRXbt21fk/5Cfr9OjRA1OmTMGGDRswe/ZsODs7k/6lRKzC+d2+fRu3blWkO7S1tUX37t1x86aONIFPcO3atSpLCo7jIJfXPmQajQYXLlxAUFCQaBtVKhV69OiBzz77DAsWLIBCoRB9mCOVrKwsXLlSoXWoUCgQHBysHSdd3L9/Hz/88IP2b7FjUklaWhqee+450fer1Wq0bdsWMTExWLRoEZycnATH5caNG5g3bx5eeOEF5kERq05hYSEuXLigLe/Tp4/geAjVS0hI0JZzHFfjQIZl46lTp3SOSfU6T6Z+iI6OhqenJ6V+kIhVOL+7d+9i48aN0Gg0KCsrQ1pamqilM8dxiIuLQ0lJCTQaDRITE5knvZXcunULbm5uaNy4sWgbs7OzERMTA41Gg+LiYiQlJSE4OFh0fSkUFRVh1apVKCkpAcdxOHnyZK3jolQqsX37dm3qgYSEBFFjAlSkMigtLYW7u7toG9VqNRYuXIiSkhKUlZUhISGhxh6erhQJtdlSvY5Go8HKlSvx4EGFqGhKSkqN8RCqZ2dnh/j4eGRmVoiiCo2JLht5nkdGRgZ8fX1F2QhQ6gdDYRXz5Z49eyIzMxMzZsyAXC5Hr169RDmWdu3aYfDgwZgzZw7Ky8sRFBQkatZy//59uLjUfPGbRatWrdC7d2/MmDEDHMdhyJAhgv8gDEmbNm20z6dQKNCxY0cMGTKEWcfFxQXh4eGIjo6GRqOBj4+Pzv276mRlZek9Lvb29hg5ciTmzp2LsrIy9OnTp8Y+6qFDh6BSqRAfH699pUoo1ERMnXHjxmHx4sWQy+V45plnMGHCBFH1pk2bhnXr1kGtVsPd3R1TpkwRVW/w4MFQKpWC2wC66kycOBGffvqpNvWDJQcSWzIU5EwQhFViFctegiCI6pDzIwjCKiHnRxCEVULOjyAIq4ScH0EQVgk5P4IgrBJyfgRBWCVWEeRM1M7FixexePHiGtfVajWef/55REREmMEqgjAeFORM6CQ5ORmbN2/GggUL0KZNG3ObQxAGhWZ+hCC//fYbNmzYgKlTp6JNmzb466+/sG3bNly6dAm2trbo2bMnRo0aBRsbG+zevRv379+HWq3G+fPn4ejoiFdffRV9+vQBANy5cwdxcXHIzMyEg4MDXn75ZZ1SUwRhKmjPj6hBVlYWli1bhrCwMPTq1QsajQaLFi2Cs7MzYmNjERUVhStXrmDr1q3aOikpKQgJCUFcXBxCQ0OxadMmqNVqlJaW4pNPPkHnzp2xceNGTJs2DfHx8TolnwjCVJDzI6pQWlqKJUuWoHPnztp8q5cuXUJ+fj7eeust2NrawtnZGaNHj8YPP/wAjuMAAK1bt0ZgYCAUCgX69u0LlUqFwsJCnDlzBnZ2dhg2bBiUSiW8vLwwcODAKnlOCMIc0LKX0MLzPFauXAk7OztMmjRJez0/Px+Ojo5o1KiR9pqrqyvKyspQWFgIAHjqqae0ZZVadhzHITs7G1lZWXj77ber9GNsrUKCqA1yfoSWr776CtevX0d0dLQ2Wx1QIWOVn5+PsrIyrQP866+/oFQqYW9vz2zT2dkZXl5eiI6O1l4rKChAeXm5cR6CIERCy14CQEUuiMOHD2PmzJlwdHSsUta2bVu0aNEC27Ztw6NHj5Cbm4sdO3YgODi4Vgn1gIAA5OTk4OjRoygvL8eDBw+waNEiHDhwwJiPQxC1QjM/AkBFzpFHjx5h3rx5Nco6dOiAWbNmYevWrZg0aRJkMhl69+6NN954o9Z2mzZtirlz52Lbtm3YsWMHlEolgoKCMHr0aGM8BkGIhuL8CIKwSmjZSxCEVULOjyAIq4ScH0EQVgk5P4IgrBJyfgRBWCXk/AiCsErI+REEYZWQ8yMIwioh50cQhFXy/45heO+C3g63AAAAAElFTkSuQmCC\n", "text/plain": [ "