diff --git a/Exercises/Session 1 - Tensor Basics (with solutions).ipynb b/Exercises/Session 1 - Tensor Basics (with solutions).ipynb new file mode 100644 index 0000000..3222ddd --- /dev/null +++ b/Exercises/Session 1 - Tensor Basics (with solutions).ipynb @@ -0,0 +1,368 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "e4f44eb7", + "metadata": {}, + "source": [ + "$\\def\\tcoreleft{\\underset{\\tiny\\mid}{\\textcolor{MidnightBlue}{⦸}}}$\n", + "$\\def\\tcorecenter{\\underset{\\tiny\\mid}{\\textcolor{RedOrange}{⦿}}}$\n", + "$\\def\\tcoreright{\\underset{\\tiny\\mid}{\\textcolor{MidnightBlue}{\\oslash}}}$\n", + "

TMQS Workshop 2024 @ Zuse Institute Berlin

\n", + "

Summer School on Tensor Methods for Quantum Simulation

\n", + "

June 3 - 5, 2024

\n", + "

$\\tcoreleft - \\tcoreleft - \\tcoreleft - \\cdots - \\tcorecenter - \\cdots - \\tcoreright - \\tcoreright$

\n", + "
" + ] + }, + { + "cell_type": "markdown", + "id": "58a30939-f570-45cd-a736-d6f21aeb2a0c", + "metadata": {}, + "source": [ + "***" + ] + }, + { + "cell_type": "markdown", + "id": "6c3441f3-b5a9-42c4-ba21-2bc682b0d8ac", + "metadata": {}, + "source": [ + "## **Session 1 - Tensor Basics**" + ] + }, + { + "cell_type": "markdown", + "id": "6f6cc702", + "metadata": {}, + "source": [ + "***" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "c8096861-c33f-49b8-9bda-62529ce31bc3", + "metadata": {}, + "source": [ + "## Exercise 1.1\n", + "\n", + "... a bit of theory to begin with.\n", + "\n", + "The following exercises are intended to demonstrate to you that one can handle tensors in a very natural way, much like with classical vectors and matrices. In the end, standard linear algebra provides the foundational understanding necessary to grasp the concepts of tensors effectively. By extending basic concepts, tensors offer a powerful framework for describing and manipulating multidimensional data in various fields such as physics, engineering, and machine learning.\n", + "\n", + "Suppose $T$ and $U$ are tensors in $\\mathbb{C}^{N}$, where $N = (n_1, \\dots , n_d)$ and $d \\in \\mathbb{N}$. Furthermore, let $G, H \\in \\mathbb{C}^{N \\times N}$. Show that\n", + "\n", + "**a)**$\\quad T^\\dagger \\cdot U = \\langle \\text{vec}(T), \\text{vec}(U) \\rangle$\n", + "\n", + "**b)**$\\quad \\lVert T \\rVert_F = \\lVert \\text{vec}(T) \\rVert_2$\n", + "\n", + "**c)**$\\quad \\text{tr}(G^\\dagger \\cdot H) = \\langle \\text{vec}(G), \\text{vec}(H) \\rangle$" + ] + }, + { + "cell_type": "markdown", + "id": "bc6abad4-3b47-447c-9e87-86b44381164b", + "metadata": {}, + "source": [ + "***" + ] + }, + { + "cell_type": "markdown", + "id": "11976a26-d5e7-4ca5-b02b-ce4709527c8c", + "metadata": {}, + "source": [ + "$\\textcolor{red}{\\textbf{SOLUTION:}}$" + ] + }, + { + "cell_type": "markdown", + "id": "742e700e-8485-4923-870b-d34994d2978b", + "metadata": {}, + "source": [ + "**a)**$\\quad$$T^\\dagger \\cdot U = \\sum_{i_1 = 1}^{n_1} \\cdot \\sum_{i_d = 1}^{n_d} \\overline{T}_{i_1, \\dots, i_d} \\cdot U_{i_1, \\dots, i_d} = \\sum_{i=1}^{n_1 \\cdot \\ldots \\cdot n_d} \\overline{\\text{vec}(T)_i} \\cdot \\text{vec}(U)_i = \\langle \\text{vec}(T), \\text{vec}(U) \\rangle$\n", + "\n", + "**b)**$\\quad$$\\lVert T \\rVert_F^2 = \\sum_{i_1 = 1}^{n_1} \\cdots \\sum_{i_d = 1}^{n_d} |T_{i_1, \\dots, i_d}|^2 = \\sum_{i=1}^{n_1 \\cdot \\ldots \\cdot n_d} |\\text{vec}(T)_i|^2 = \\lVert \\text{vec}(T) \\rVert_2^2$\n", + "\n", + "**c)**$\\quad$$\\text{tr}(G^\\dagger \\cdot H) = \\sum_{i_1=1}^{n_1} \\cdots \\sum_{i_d=1}^{n_d} \\sum_{j_1=1}^{n_1} \\cdots \\sum_{j_d=1}^{n_d} G^\\dagger_{i_1, \\dots, i_d, j_1 , \\dots, j_d} \\cdot H_{j_1, \\dots, j_d, i_1, \\dots, i_d} = \\sum_{i=1}^{n_1^2 \\cdot \\ldots \\cdot n_d^2} \\overline{\\text{vec}(G)_i} \\cdot \\text{vec}(H)_i = \\langle \\text{vec}(G), \\text{vec}(H) \\rangle$" + ] + }, + { + "cell_type": "markdown", + "id": "57638421-d883-4f7f-ab0f-aae0505cffca", + "metadata": {}, + "source": [ + "***" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "41a28d39-ecbe-4a8b-9ed0-161d36c14c73", + "metadata": {}, + "source": [ + "## Exercise 1.2 \n", + "\n", + "... now onto the basic implementation\n", + "\n", + "We turn to the practical implementation of tensors and their limitations. Tensors in full format can be easily created and multiplied/contracted in Python using the standard package NumPy:\n", + "\n", + "> import numpy as np\n", + "\n", + "In addition, we also need the Matplotlib and time package for this task:\n", + "\n", + "> import matplotlib.pyplot as plt\n", + "> \n", + "> import time\n", + "\n", + "Create random order-$d$ tensors in $\\mathbb{R}^{N \\times N}$ with $N=(2,\\dots,2)$ for different dimensions $d = 1, 2, \\dots, 10$:\n", + "\n", + "> T = np.random.random_sample(2*d*[2])\n", + "\n", + "Compute the storage (```T.size * T.itemsize```) of every instance and plot the results:\n", + "\n", + "> plt.plot(np.arange(1,11), storage)\n", + "\n", + "Rerun the expriment and reshape every instance into an array in $\\mathbb{R}^{2^d \\times 2^d}$ and compute an SVD:\n", + "\n", + "> T_mat = T.reshape([2** d, 2**d])\n", + ">\n", + "> U, S, Vt = np.linalg.svd(T)\n", + "\n", + "Measure the CPU time using:\n", + "\n", + "> start_time = time.time()\n", + ">\n", + "> ...\n", + ">\n", + "> end_time = time.time() - start_time\n", + "\n", + "Plot the results and interpret both plots in terms of storage consumption and computational complexity, respectively. Given a matrix $M \\in \\mathbb{R}^{m \\times n}$, $m \\geq n$, the overall cost for computing an SVD is $O(m n^2)$." + ] + }, + { + "cell_type": "markdown", + "id": "83e132fb-95b5-43a8-98ca-0ed7c397767d", + "metadata": {}, + "source": [ + "***" + ] + }, + { + "cell_type": "markdown", + "id": "1cb59312-1033-4dd0-886c-7a034eae661b", + "metadata": {}, + "source": [ + "$\\textcolor{red}{\\textbf{SOLUTION:}}$" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "id": "f5ca5182-ecb4-4b2b-815b-4b7c44af9920", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[]" + ] + }, + "execution_count": 12, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAioAAAHHCAYAAACRAnNyAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuNSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/xnp5ZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA9wklEQVR4nO3deXhTZf7+8Tvd0tINCl0tlEWwUnYQv4oKCooM4KAjLoNaFZ0ZRRFBZmD4uTCKKI4MIypuCDiIu6ijI4sIKK4sLhSUXS1LW6DQlW7J+f1REhpaoIGWc5K8X9eVi+ack+TTppq7z/M5z7EZhmEIAADAgoLMLgAAAOBYCCoAAMCyCCoAAMCyCCoAAMCyCCoAAMCyCCoAAMCyCCoAAMCyCCoAAMCyCCoAAMCyCCoAYFH9+vVTv379zC4DMBVBBTgNFixYoBkzZphdBixo48aNeuihh/TLL7+YXQpgSTau9QM0viFDhigrK4sPI9Ty9ttva/jw4Vq+fHmt0ZOKigpJUlhYmAmVAdYQYnYBAE5OVVWVnE4nH2J+jPcWYOoHOGVFRUUaM2aMWrduLbvdroSEBF166aVat26dpOo+g48++ki//vqrbDabbDabWrdu7X58Xl6eRo4cqcTERIWHh6tr166aN2+ex2v88ssvstls+uc//6kZM2aoXbt2stvt2rhxoyoqKvTAAw+oZ8+eio2NVWRkpC688EItX768Vq379+/XjTfeqJiYGDVt2lSZmZn64YcfZLPZNHfuXI9jf/75Z1199dWKi4tTeHi4evXqpQ8++KBePxOn06l///vf6ty5s8LDwxUfH6/LL79ca9ascR9TVVWlhx9+2P29tG7dWn//+99VXl7u8VytW7fWkCFDtGrVKvXu3Vvh4eFq27atXnnlFY/jKisrNXnyZLVv317h4eFq3ry5LrjgAi1dutR9zLF6Pm6++WaP96Tmz/uZZ55R27Zt1aRJE1122WXKzs6WYRh6+OGHlZqaqoiICP3+979Xfn5+nXUvWbJE3bp1U3h4uDp27Kh3333XfczcuXM1fPhwSdLFF1/s/v1YsWLFMev19vflhRdecP+MzznnHK1evbruNw2wKEZUgFP0l7/8RW+//bbuuusudezYUfv379eqVav0008/qUePHpo0aZIKCgq0c+dO/etf/5IkRUVFSZIOHTqkfv36aevWrbrrrrvUpk0bvfXWW7r55pt18OBB3XPPPR6vNWfOHJWVlelPf/qT7Ha74uLiVFhYqJdeeknXX3+9br/9dhUVFWn27NkaOHCgvv32W3Xr1k1SdXgYOnSovv32W91xxx1KT0/X+++/r8zMzFrf04YNG9SnTx+dccYZmjBhgiIjI/Xmm29q2LBheuedd3TllVce92cycuRIzZ07V4MGDdJtt92mqqoqff755/r666/Vq1cvSdJtt92mefPm6eqrr9a4ceP0zTffaOrUqfrpp5+0cOFCj+fbunWrrr76ao0cOVKZmZl6+eWXdfPNN6tnz57KyMiQJD300EOaOnWqbrvtNvXu3VuFhYVas2aN1q1bp0svvdT7N1bSq6++qoqKCt19993Kz8/XtGnTdM011+iSSy7RihUr9Le//U1bt27VzJkzdd999+nll1/2ePyWLVt07bXX6i9/+YsyMzM1Z84cDR8+XIsWLdKll16qiy66SKNHj9ZTTz2lv//97zr77LMlyf3v0bz9fVmwYIGKior05z//WTabTdOmTdNVV12l7du3KzQ09KR+JsBpZwA4JbGxscaoUaOOe8zgwYONtLS0WttnzJhhSDLmz5/v3lZRUWGcd955RlRUlFFYWGgYhmHs2LHDkGTExMQYeXl5Hs9RVVVllJeXe2w7cOCAkZiYaNx6663ube+8844hyZgxY4Z7m8PhMC655BJDkjFnzhz39v79+xudO3c2ysrK3NucTqdx/vnnG+3btz/u9/rpp58akozRo0fX2ud0Og3DMIzvv//ekGTcdtttHvvvu+8+Q5Lx6aefurelpaUZkozPPvvMvS0vL8+w2+3GuHHj3Nu6du1qDB48+Li19e3b1+jbt2+t7ZmZmR7vj+vnHR8fbxw8eNC9feLEiYYko2vXrkZlZaV7+/XXX2+EhYV5/Lxcdb/zzjvubQUFBUZycrLRvXt397a33nrLkGQsX778hPV6+/vSvHlzIz8/333s+++/b0gy/vvf/x735wRYCVM/wClq2rSpvvnmG+3evdvrx/7vf/9TUlKSrr/+eve20NBQjR49WsXFxVq5cqXH8X/4wx8UHx/vsS04ONjdy+B0OpWfn6+qqir16tXLPf0kSYsWLVJoaKhuv/1297agoCCNGjXK4/ny8/P16aef6pprrlFRUZH27dunffv2af/+/Ro4cKC2bNmiXbt2HfN7euedd2Sz2fTggw/W2mez2dzftySNHTvWY/+4ceMkSR999JHH9o4dO+rCCy9034+Pj9dZZ52l7du3u7c1bdpUGzZs0JYtW45Zm7eGDx+u2NhY9/1zzz1XknTDDTcoJCTEY3tFRUWtn0tKSorH6FNMTIxuuukmfffdd8rJyfG6Hm9/X6699lo1a9bMfd/1M6z5cwOsjqACnKJp06YpKytLLVu2VO/evfXQQw/V+4Pg119/Vfv27RUU5Pmfomvo/9dff/XY3qZNmzqfZ968eerSpYu7NyM+Pl4fffSRCgoKPF4rOTlZTZo08XjsmWee6XF/69atMgxD999/v+Lj4z1urvCRl5d3zO9p27ZtSklJUVxc3HG/76CgoFqvnZSUpKZNm9b6vlu1alXrOZo1a6YDBw647//jH//QwYMH1aFDB3Xu3Fnjx4/Xjz/+eMwa6uPo13WFlpYtW9a5vWY9UvXP1hXOXDp06CBJJ3UGmLe/L0fX7wotR9cJWBlBBThF11xzjbZv366ZM2cqJSVFTzzxhDIyMvTxxx83+GtFRETU2jZ//nzdfPPNateunWbPnq1FixZp6dKluuSSS+R0Or1+Dddj7rvvPi1durTO29EB42Qd/SF+LMHBwXVuN2qsrnDRRRdp27Ztevnll9WpUye99NJL6tGjh1566aUTvp7D4fDqdetTjxX4Sp3A8dBMCzSA5ORk3XnnnbrzzjuVl5enHj16aMqUKRo0aJCkY39ApqWl6ccff5TT6fT4K/nnn3927z+Rt99+W23bttW7777r8TpHT72kpaVp+fLlKi0t9RhV2bp1q8dxbdu2lVQ9pTBgwIATvv7R2rVrp8WLFys/P/+YoyppaWlyOp3asmWLR+Nobm6uDh48WK/vuy5xcXG65ZZbdMstt6i4uFgXXXSRHnroId12222SqkcU6hrtOnokoqG4Rqdqvi+bN2+WJPdZRvUNa1LD/L4AvoYRFeAUOBwOj+kVSUpISFBKSorHabaRkZG1jpOk3/3ud8rJydEbb7zh3lZVVaWZM2cqKipKffv2PWENrr+aa/6V/M033+irr77yOG7gwIGqrKzUiy++6N7mdDr1zDPP1Kq/X79+ev7557Vnz55ar7d3797j1vOHP/xBhmFo8uTJtfa5avzd734nSbVW650+fbokafDgwcd9jbrs37/f435UVJTOPPNMj/ehXbt2+vnnnz2+hx9++EFffPGF169XH7t37/Y4g6mwsFCvvPKKunXrpqSkJEnVvxuSdPDgwRM+X0P8vgC+hhEV4BQUFRUpNTVVV199tbp27aqoqCh98sknWr16tZ588kn3cT179tQbb7yhsWPH6pxzzlFUVJSGDh2qP/3pT3r++ed18803a+3atWrdurXefvttffHFF5oxY4aio6NPWMOQIUP07rvv6sorr9TgwYO1Y8cOPffcc+rYsaOKi4vdxw0bNky9e/fWuHHjtHXrVqWnp+uDDz5wr/9R8y/7Z555RhdccIE6d+6s22+/XW3btlVubq6++uor7dy5Uz/88MMx67n44ot144036qmnntKWLVt0+eWXy+l06vPPP9fFF1+su+66S127dlVmZqZeeOEFHTx4UH379tW3336refPmadiwYbr44ou9fi86duyofv36qWfPnoqLi9OaNWvcp4273HrrrZo+fboGDhyokSNHKi8vT88995wyMjJUWFjo9WueSIcOHTRy5EitXr1aiYmJevnll5Wbm6s5c+a4j+nWrZuCg4P1+OOPq6CgQHa7XZdccokSEhJqPV9D/L4APse8E44A31deXm6MHz/e6Nq1qxEdHW1ERkYaXbt2NZ599lmP44qLi40//vGPRtOmTQ1JHqfC5ubmGrfccovRokULIywszOjcubPHqcKGceR00yeeeKJWDU6n03j00UeNtLQ0w263G927dzc+/PDDWqfcGoZh7N271/jjH/9oREdHG7GxscbNN99sfPHFF4Yk4/XXX/c4dtu2bcZNN91kJCUlGaGhocYZZ5xhDBkyxHj77bdP+HOpqqoynnjiCSM9Pd0ICwsz4uPjjUGDBhlr1651H1NZWWlMnjzZaNOmjREaGmq0bNnSmDhxoscpvoZRfZpvXacdH33q7iOPPGL07t3baNq0qREREWGkp6cbU6ZMMSoqKjweN3/+fKNt27ZGWFiY0a1bN2Px4sXHPD356J/38uXLDUnGW2+95bF9zpw5hiRj9erVtepevHix0aVLF8Nutxvp6em1HmsYhvHiiy8abdu2NYKDgz1OVa7rdOpT/X2RZDz44IO1tgNWxbV+gAD33nvv6corr9SqVavUp08fs8vxG61bt1anTp304Ycfml0K4NPoUQECyKFDhzzuOxwOzZw5UzExMerRo4dJVQHAsdGjAgSQu+++W4cOHdJ5552n8vJyvfvuu/ryyy/16KOP1nnqMwCYjaACBJBLLrlETz75pD788EOVlZXpzDPP1MyZMz0aTgHASuhRAQAAlkWPCgAAsCyCCgAAsCyf7lFxOp3avXu3oqOjvVqGGgAAmMcwDBUVFSklJaXWRTaP5tNBZffu3bWuYgoAAHxDdna2UlNTj3uMTwcV13LR2dnZiomJMbkaAABQH4WFhWrZsmW9Lvvg00HFNd0TExNDUAEAwMfUp22DZloAAGBZBBUAAGBZBBUAAGBZBBUAAGBZBBUAAGBZBBUAAGBZBBUAAGBZBBUAAGBZBBUAAGBZBBUAAGBZBBUAAGBZBBUAAGBZBBUAAFDLwdIKbcktksNpmFoHQQUAANSydGOuLv3XZ7p5zrem1kFQAQAAtWTtKpAkdUiMNrUOggoAAKgla3ehJKnzGbGm1kFQAQAAHhxOQxsPB5VOBBUAAGAl2/YW61ClQ5FhwWrbItLUWggqAADAw/qd1f0pGSmxCgqymVoLQQUAAHhYf7iR1uxpH4mgAgAAjpLlDioxJldCUAEAADU4nIY2WOSMH4mgAgAAath+uJG2SViw2sZHmV0OQQUAABzh6k/pmByjYJMbaSWCCgAAqMFKjbSSyUHF4XDo/vvvV5s2bRQREaF27drp4YcflmGYewEkAAAC1YZd1ulPkaQQM1/88ccf16xZszRv3jxlZGRozZo1uuWWWxQbG6vRo0ebWRoAAAHH6TS0YXf1iErnVIKKvvzyS/3+97/X4MGDJUmtW7fWa6+9pm+/NfdKjQAABKLt+0pUUuFQRGiw2lmgkVYyeern/PPP17Jly7R582ZJ0g8//KBVq1Zp0KBBdR5fXl6uwsJCjxsAAGgYrvVTOqZYo5FWMnlEZcKECSosLFR6erqCg4PlcDg0ZcoUjRgxos7jp06dqsmTJ5/mKgEACAzuRtoU8xd6czF1ROXNN9/Uq6++qgULFmjdunWaN2+e/vnPf2revHl1Hj9x4kQVFBS4b9nZ2ae5YgAA/JfVzviRTB5RGT9+vCZMmKDrrrtOktS5c2f9+uuvmjp1qjIzM2sdb7fbZbfbT3eZAAD4PafT0EbXirQWaaSVTB5RKS0tVVCQZwnBwcFyOp0mVQQAQGDasb9ExeVVCg8N0pkWaaSVTB5RGTp0qKZMmaJWrVopIyND3333naZPn65bb73VzLIAAAg4rkbas5NjFBJsnfVgTQ0qM2fO1P33368777xTeXl5SklJ0Z///Gc98MADZpYFAEDAWb/z8PopFupPkUwOKtHR0ZoxY4ZmzJhhZhkAAAS8rN3Wa6SVuNYPAAABz+k0LLd0vgtBBQCAAPdrfqmKyqtkDwlS+wTrNNJKBBUAAALeeos20koEFQAAAl6We6E366xI60JQAQAgwFn1jB+JoAIAQEAzDMOyZ/xIBBUAAALar/tLVVRWpbCQIHVIjDa7nFoIKgAABDB3I21StEIt1kgrEVQAAAhoVp72kQgqAAAENNcZP1ZspJUIKgAABCzDMJR1eEVaRlQAAIClZOcfUsGhSoUFW7ORViKoAAAQsFyNtOnJ0QoLsWYksGZVAACg0bmCSkaKNad9JIIKAAABy+qNtBJBBQCAgGQYhntEhaACAAAsZeeB6kba0GCbOiRFmV3OMRFUAAAIQK7RlLOSomUPCTa5mmMjqAAAEIB8oT9FIqgAABCQXCMqVl3ozYWgAgBAgKlekZYRFQAAYEG7Dh7SgdLqRtqzkqy5Iq0LQQUAgADjGk3pkGjtRlqJoAIAQMBx96dYeEVaF4IKAAABZr3rismpBBUAAGAhvtRIKxFUAAAIKLsLypRfUqGQIJvSLd5IKxFUAAAIKK7RlPaJ0QoPtXYjrURQAQAgoByZ9okxuZL6IagAABBAfOGKyTURVAAACBA1G2mtvnS+C0EFAIAAkVNYpn3FFQoOsunsZKZ+AACAhazfebiRNiHKJxppJYIKAAABw9emfSSCCgAAAcPXGmklggoAAAHBMIwjS+cTVAAAgJXkFpZrX3G5gmxSRx9ppJUIKgAABAT3irQJ0YoI841GWomgAgBAQFjvg420EkEFAICA4GtL57sQVAAACADuM35SGVEBAAAWkldYpryi6kZaX1mR1oWgAgCAn3ONprSLj1KTsBCTq/EOQQUAAD/niwu9uRBUAADwc764dL4LQQUAAD/nq420EkEFAAC/lldUptzCctl8bEVaF4IKAAB+bMPh6/u0i49SpN23GmklggoAAH7NlxtpJYIKAAB+zVeXznchqAAA4MeyGFEBAABWtK+4XHsKyqobaVN8r5FWIqgAAOC3XNM+bVpEKsoHG2klggoAAH4ra6dvT/tIBBUAAPyWr5/xIxFUAADwW768dL4LQQUAAD+0v7hcuwvKJEkZPtpIKxFUAADwS1m7q1ekbdsiUtHhoSZXc/IIKgAA+CF/mPaRCCoAAPil9X5wxo9EUAEAwC/5+tL5LgQVAAD8zIGSCu06eEiSlHGG7zbSSgQVAAD8jms0pXXzJorx4UZaiaACAIDf8ZdpH4mgAgCA3/H1KybXRFABAMDPZO0mqAAAAAs6WFqh7HxXIy1BBQAAWEjWruoVadOaN1FshG830koEFQAA/Io/NdJKFggqu3bt0g033KDmzZsrIiJCnTt31po1a8wuCwAAn+RPjbSSFGLmix84cEB9+vTRxRdfrI8//ljx8fHasmWLmjVrZmZZAAD4rPUElYbz+OOPq2XLlpozZ457W5s2bUysCAAA31VQWqnf8kslSRkpvr0irYupUz8ffPCBevXqpeHDhyshIUHdu3fXiy++aGZJAAD4LNdpyS3jItS0SZjJ1TQMU4PK9u3bNWvWLLVv316LFy/WHXfcodGjR2vevHl1Hl9eXq7CwkKPGwAAqOZv0z6SyVM/TqdTvXr10qOPPipJ6t69u7KysvTcc88pMzOz1vFTp07V5MmTT3eZAAD4BH8740cyeUQlOTlZHTt29Nh29tln67fffqvz+IkTJ6qgoMB9y87OPh1lAgDgEzYwotKw+vTpo02bNnls27x5s9LS0uo83m63y263n47SAADwKYVllfplf3UjbacU/wkqpo6o3Hvvvfr666/16KOPauvWrVqwYIFeeOEFjRo1ysyyAADwOa71U1KbRahZpH800komB5VzzjlHCxcu1GuvvaZOnTrp4Ycf1owZMzRixAgzywIAwOf420JvLqZO/UjSkCFDNGTIELPLAADAp60/fI0ff2qklSywhD4AADh1/jqiQlABAMDHFZZVase+EkmMqAAAAIvZcHja54ymEYrzo0ZaiaACAIDPy3Iv9OYf1/epiaACAICP88el810IKgAA+DjXxQj9rT9FIqgAAODTisur3I20jKgAAABL2bCrQIYhpcSGq3mU/11mhqACAIAP88crJtdEUAEAwIf560JvLgQVAAB8mHtEJZWgAgAALKS4vErbXSvSphBUAACAhWzcXSjDkJJiwhUf7X+NtBJBBQAAn+XvjbQSQQUAAJ+1wc8baSWCCgAAPsu9dH6q/13jx4WgAgCADyqtqNK2vcWSmPoBAAAWs3F3oZyGlBhjV0J0uNnlNBqvgsry5cv15JNP6osvvpAkPf/882rVqpXi4+N1++2369ChQ41SJAAA8OTPV0yuKaS+B7744ou644471KZNG02aNEkPPvigpkyZohtvvFFBQUGaP3++mjdvrscee6wx6wUAAAqMM34kL0ZU/v3vf+tf//qXtmzZovfee08PPPCAnnnmGc2aNUvPPPOMXnrpJb399tuNWSsAADjMtXS+vy705lLvoLJ9+3ZdccUVkqTLL79cNptNvXv3du8/99xzlZ2d3fAVAgAAD6UVVdqaV91I29lPl853qXdQKSsrU0REhPu+3W6X3W73uF9VVdWw1QEAgFp+2lPdSBsfbVdijP820kpe9KjYbDYVFRUpPDxchmHIZrOpuLhYhYWFkuT+FwAANK71OwOjkVbyIqgYhqEOHTp43O/evbvHfZvN1rDVAQCAWrJ2Vw8O+HsjreRFUFm+fHlj1gEAAOopK0BOTZa8CCp9+/ZtzDoAAEA9lFU6tMXVSBsAQYWVaQEA8CEb9xTK4TTUIsquxBj7iR/g4+o9ohIcHFyv4xwOx0kXAwAAju/ItE9MQPSGetVMm5aWpszMTI8mWgAAcPoE0hk/khdB5dtvv9Xs2bP173//W23atNGtt96qESNGqFmzZo1ZHwAAqMG1dH5GgASVeveo9OrVS7NmzdKePXs0duxYLVy4UKmpqbruuuu0dOnSxqwRAAAo8BpppZNopg0PD9cNN9ygZcuWKSsrS3l5ebr88suVn5/fGPUBAIDDfjrcSNs8MkzJsf69Iq1Lvad+atq5c6fmzp2ruXPnqrS0VOPHj1dMTExD1wYAAGqoudBbIDTSSl4ElYqKCi1cuFCzZ8/W559/rkGDBmnGjBkaNGhQvc8IAgAAJy8rwBppJS+CSnJysqKjo5WZmalnn31WCQkJkqSSkhKP4xhZAQCgcbgaaQNh6XwXm2EYRn0ODAo60s5S13CT61o/p3MdlcLCQsXGxqqgoICABADwa2WVDnV6cLGqnIa+mHCJzmgaYXZJJ82bz2+u9QMAgA/YlFOkKqehuMgwpQRII63EtX4AAPAJNad9AqWRVuJaPwAA+ISaS+cHEoIKAAA+wD2ikhI4jbQSQQUAAMsrr3Joc26RpMA640ciqAAAYHmbcopU6TDUtEmoUpv57tk+J4OgAgCAxWXtql6RtnOANdJKJ7GEfklJiR577DEtW7ZMeXl5cjqdHvu3b9/eYMUBAIDAXOjNxeugctttt2nlypW68cYblZycHHDJDgCA0+3IGT8ElRP6+OOP9dFHH6lPnz6NUQ8AAKihosqpTTnVjbSBGFS87lFp1qyZ4uLiGqMWAABwlM25RapwOBUbEXiNtNJJBJWHH35YDzzwgEpLSxujHgAAUMP6GtM+gdhu4fXUz5NPPqlt27YpMTFRrVu3VmhoqMf+devWNVhxAAAEukBupJVOIqgMGzasEcoAAAB1yXIHlcBaOt/F66Dy4IMPNkYdAADgKBVVTv28J3AbaSUWfAMAwLK25FU30saEh6hVXBOzyzFFvUZU4uLitHnzZrVo0ULNmjU7bjNPfn5+gxUHAEAgy6rRnxKIjbRSPYPKv/71L0VHR0uSZsyY0Zj1AACAw9YH8EJvLvUKKpmZmXV+DQAAGs/6w9f4CdQzfiR6VAAAsKRKh1M/7TlyMcJARVABAMCCtuQWq6LKqejwEKU1D8xGWomgAgCAJbkbaVMCt5FWIqgAAGBJ6wN8oTcXr4PKnDlzuM4PAACNLNCXznfxOqhMmDBBSUlJGjlypL788svGqAkAgIBWRSOtm9dBZdeuXZo3b5727dunfv36KT09XY8//rhycnIaoz4AAALO1r3FKq9yKsoeotbNI80ux1ReB5WQkBBdeeWVev/995Wdna3bb79dr776qlq1aqUrrrhC77//vpxOZ2PUCgBAQFi/s3raJyMlRkFBgdtIK51iM21iYqIuuOACnXfeeQoKCtL69euVmZmpdu3aacWKFQ1UIgAAgSWLFWndTiqo5Obm6p///KcyMjLUr18/FRYW6sMPP9SOHTu0a9cuXXPNNaxgCwDASXIvnZ9KUPE6qAwdOlQtW7bU3Llzdfvtt2vXrl167bXXNGDAAElSZGSkxo0bp+zs7AYvFgAAf1flcGrjHpbOd6nXtX5qSkhI0MqVK3Xeeecd85j4+Hjt2LHjlAoDACAQbdtborLK6kbaNgHeSCt5OaJSWVmpX375RS1atDjucTabTWlpaadUGAAAgcg17dORRlpJXgaV0NBQ/fjjj41SyGOPPSabzaYxY8Y0yvMDAOALai6dj5PoUbnhhhs0e/bsBi1i9erVev7559WlS5cGfV4AAHzNkUbawF4638XrHpWqqiq9/PLL+uSTT9SzZ09FRnrOn02fPt2r5ysuLtaIESP04osv6pFHHvG2HAAA/IbDaWjjblakrcnroJKVlaUePXpIkjZv3uyx72Su7jhq1CgNHjxYAwYMIKgAAALa9r3FOlTpUJOwYLVpEWV2OZbgdVBZvnx5g73466+/rnXr1mn16tX1Or68vFzl5eXu+4WFhQ1WCwAAZnNN+2SkxCiYRlpJp7gy7c6dO7Vz586Temx2drbuuecevfrqqwoPD6/XY6ZOnarY2Fj3rWXLlif12gAAWBFXTK7N66DidDr1j3/8Q7GxsUpLS1NaWpqaNm2qhx9+2Ktr/Kxdu1Z5eXnq0aOHQkJCFBISopUrV+qpp55SSEiIHA5HrcdMnDhRBQUF7huLygEA/AlL59fm9dTPpEmTNHv2bD322GPq06ePJGnVqlV66KGHVFZWpilTptTrefr376/169d7bLvllluUnp6uv/3tbwoODq71GLvdLrvd7m3JAABYnsNpaAONtLV4HVTmzZunl156SVdccYV7W5cuXXTGGWfozjvvrHdQiY6OVqdOnTy2RUZGqnnz5rW2AwDg73bsK1ZpRXUjbdt4GmldvJ76yc/PV3p6eq3t6enpys/Pb5CiAAAINO4VaZNppK3J6xGVrl276umnn9ZTTz3lsf3pp59W165dT6mYFStWnNLjAQDwVet3ciHCungdVKZNm6bBgwfrk08+cV+Y8KuvvlJ2drb+97//NXiBAAAEgizO+KmT11M/ffv21ebNm3XllVfq4MGDOnjwoK666ipt2rRJF154YWPUCACAX3M6DW3YzRk/dfF6REWSUlJS6t00CwAAjm/H/hKVVDgUHhqkdvGRJ35AADmpoFJWVqYff/xReXl5tdZOqXk2EAAAOLGsGo20IcGntBar3/E6qCxatEg33XST9u3bV2ufzWarc6E2AABwbOt3Mu1zLF7HtrvvvlvDhw/Xnj175HQ6PW6EFAAAvMfS+cfmdVDJzc3V2LFjlZiY2Bj1AAAQUJw1V6RNJagczeugcvXVV7PeCQAADeSX/SUqLq9SeGiQzmRF2lq87lF5+umnNXz4cH3++efq3LmzQkNDPfaPHj26wYoDAMDfuaZ9zqaRtk5eB5XXXntNS5YsUXh4uFasWCGb7cgyvzabjaACAIAX3Au9pTDtU5eTunry5MmTNWHCBAUFkfwAADgVrhEVzvipm9dJo6KiQtdeey0hBQCAU+R0Gtqwi2v8HI/XaSMzM1NvvPFGY9QCAEBA+S2/VEXlVQoLCVL7RBpp6+L11I/D4dC0adO0ePFidenSpVYz7fTp0xusOAAA/FnNRtpQGmnr5HVQWb9+vbp37y5JysrK8thXs7EWAAAcX5a7PyXG5Eqsy+ugsnz58saoAwCAgEMj7Ymd9DjT1q1btXjxYh06dEiSZBhGgxUFAIC/MwzjyKnJBJVj8jqo7N+/X/3791eHDh30u9/9Tnv27JEkjRw5UuPGjWvwAgEA8Ee/5ZeqsKxKYcFB6pAYbXY5luV1ULn33nsVGhqq3377TU2aNHFvv/baa7Vo0aIGLQ4AAH/lmvZJT46mkfY4vO5RWbJkiRYvXqzU1FSP7e3bt9evv/7aYIUBAODPuGJy/Xgd4UpKSjxGUlzy8/Nlt9sbpCgAAPyda6E3GmmPz+ugcuGFF+qVV15x37fZbHI6nZo2bZouvvjiBi0OAAB/ZBgGZ/zUk9dTP9OmTVP//v21Zs0aVVRU6K9//as2bNig/Px8ffHFF41RIwAAfmXngUMqOFRJI209eD2i0qlTJ23evFkXXHCBfv/736ukpERXXXWVvvvuO7Vr164xagQAwK+4RlPOSopWWAiNtMfj9YiKJMXGxmrSpEkNXQsAAAGBRtr6O6mgUlZWph9//FF5eXlyOp0e+6644ooGKQwAAH+VRX9KvXkdVBYtWqSbbrpJ+/btq7XPZrPJ4XA0SGEAAPgjGmm94/XE2N13363hw4drz549cjqdHjdCCgAAx7fzwCEdLK1UaLBNHZKizC7H8rwOKrm5uRo7dqwSExMbox4AAPyaa9qnQ2K07CHBJldjfV4HlauvvlorVqxohFIAAPB/TPt4x+selaefflrDhw/X559/rs6dOys0NNRj/+jRoxusOAAA/E3W7uoVaTnjp368DiqvvfaalixZovDwcK1YsUI2m829z2azEVQAADgGwzA448dLXgeVSZMmafLkyZowYYKCglikBgCA+tpdUKb8kgqFBNl0VhIr0taH10mjoqJC1157LSEFAAAvrd95pJE2PJRG2vrwOm1kZmbqjTfeaIxaAADwa0z7eM/rqR+Hw6Fp06Zp8eLF6tKlS61m2unTpzdYcQAA+BP30vmpBJX68jqorF+/Xt27d5ckZWVleeyr2VgLAACOoJH25HgdVJYvX94YdQAA4Nf2FJRpf0mFgoNsSqeRtt7oiAUA4DRwTfu0T4iikdYLBBUAAE6DDUz7nBSCCgAAp4F76Xwaab1CUAEAoJEZhqH1u1g6/2QQVAAAaGS5heXaV1yu4CCbOibHmF2OTyGoAADQyGikPXkEFQAAGpl7oTemfbxGUAEAoJGx0NvJI6gAANDIGFE5eQQVAAAaUW5hmfYWlSvIJhppTwJBBQCARrR+Z/VoypkJUYoIo5HWWwQVAAAaUdZupn1OBUEFAIBGRCPtqSGoAADQiNYTVE4JQQUAgEaSV1Sm3MLDjbQpNNKeDIIKAACNxDXt0y4+Sk3CQkyuxjcRVAAAaCTrd1ZfiJBpn5NHUAEAoJGw0NupI6gAANBI3Gf8pBJUThZBBQCARrC3qFw5hWWysSLtKSGoAADQCFwLvbVtEalIO420J4ugAgBAI8jayfopDYGgAgBAI6CRtmEQVAAAaAQsnd8wCCoAADSw/cXl2l1Q3UibQVA5JQQVAAAamGvap02LSEXRSHtKCCoAADQwpn0aDkEFAIAGxhWTGw5BBQCABpa1q/oaP5zxc+oIKgAANKD8kgrtOnhIktQxhRVpTxVBBQCABpRVo5E2JjzU5Gp8H0EFAIAGtHLzXklM+zQUU4PK1KlTdc455yg6OloJCQkaNmyYNm3aZGZJAACctE9/ztXLX+yQJA3unGRyNf7B1KCycuVKjRo1Sl9//bWWLl2qyspKXXbZZSopKTGzLAAAvLZjX4nuef17GYZ003lpurxTstkl+QVTV6FZtGiRx/25c+cqISFBa9eu1UUXXWRSVQAAeKe4vEp/emWNisqq1Cutmf7f4I5ml+Q3LLVcXkFBdQNSXFxcnfvLy8tVXl7uvl9YWHha6gIA4FgMw9Bf3/5BW/KKlRhj17M39FBYCC2gDcUyP0mn06kxY8aoT58+6tSpU53HTJ06VbGxse5by5YtT3OVAAB4em7ldv1vfY5Cg216dkRPJUSHm12SX7FMUBk1apSysrL0+uuvH/OYiRMnqqCgwH3Lzs4+jRUCAODps8179cTinyVJk6/opJ5pzUyuyP9YYurnrrvu0ocffqjPPvtMqampxzzObrfLbrefxsoAAKjbb/tLdfdr38lpSNf3bqk/ntvK7JL8kqlBxTAM3X333Vq4cKFWrFihNm3amFkOAAD1UlpRpT/9Z40KDlWqW8umeuiKDLNL8lumBpVRo0ZpwYIFev/99xUdHa2cnBxJUmxsrCIiIswsDQCAOhmGoQnvrNfPOUVqEWXXczf0lD0k2Oyy/JapPSqzZs1SQUGB+vXrp+TkZPftjTfeMLMsAACOafaqHfrgh90KCbLp2RE9lBRL82xjMn3qBwAAX/Hltn2a+nF18+z9Qzqqd5u6l9NAw7HMWT8AAFjZroOHdNeC7+RwGrqqxxm66bw0s0sKCAQVAABOoKzSob/8Z63ySyrU6YwYPXplZ9lsNrPLCggEFQAAjsMwDE1amKX1uwoUFxmm527oqfBQmmdPF4IKAADH8Z+vf9U763YqyCY9fX13pTZrYnZJAYWgAgDAMXy7I1//+O9GSdLff3e2zj+zhckVBR6CCgAAdcgpKNOdr65TldPQ0K4pGnkBi5KagaACAMBRyqsc+sv8tdpXXK70pGg9/geaZ81CUAEA4CgPfbBB32cfVGxEqF64sZeahFni0ngBiaACAEANC775Ta99my2bTXrq+u5q1ZzmWTMRVAAAOGzdbwf04AdZkqTxA89S3w7xJlcEggoAAJLyisp0x/y1qnQYGtQpSXf0bWd2SRBBBQAAVVQ5NerVdcotLFf7hCg9MbwrzbMWQVABAAS8Rz7aqNW/HFC0PUTP39hTUXaaZ62CoAIACGhvrcnWK1/9KkmacV03tY2PMrki1ERQAQAErB93HtSk96qbZ+8d0EH9z040uSIcjaACAAhI+4vL9Zf/rFVFlVMDzk7U3ZecaXZJqANBBQAQcKocTo1asE67C8rUtkWkpl/bVUFBNM9aEUEFABBwpn78s77enq/IsGC9cFNPxYSHml0SjoGgAgAIKO9/v0uzV+2QJD15TVedmRBtckU4HoIKACBgbNhdoL+986Mk6a6Lz9TlnZJNrggnQlABAASEAyUV+sv8tSqrdKpvh3jde2kHs0tCPRBUAAB+z+E0NPr175Sdf0it4proqeu6K5jmWZ9AUAEA+L1/Ltmkz7fsU0RodfNsbBOaZ30FQQUA4Nc++nGPZq3YJkmadnUXpSfFmFwRvEFQAQD4rU05RRr/9g+SpD9f1FZDu6aYXBG8RVABAPilgkOV+vN/1qi0wqE+ZzbX+IFnmV0STgJBBQDgd5xOQ2Ne/06/7C/VGU0jNPP6HgoJ5iPPF/GuAQD8zoxlW7R8017ZQ4L0/I09FRcZZnZJOEkEFQCAX1myIUdPLdsiSZp6VWd1OiPW5IpwKggqAAC/sTWvWGPfrG6evfn81rqqR6rJFeFUEVQAAH6hqKy6eba4vEq928Rp0uCzzS4JDYCgAgDweU6noXFv/qBte0uUHBuuZ/7YQ6E0z/oF3kUAgM97dsVWLdmYq7DgIM26oafio+1ml4QGQlABAPi05T/n6cmlmyVJjwzrpG4tm5pbEBoUQQUA4LN+2Vei0a9/J8OQRpzbStec09LsktDACCoAAJ9UUl6lP/9nrYrKqtQzrZkeHJphdkloBAQVAIDPMQxDf33nR23KLVJCtF2zRvRQWAgfaf6IdxUA4HNe+Gy7Pvpxj0KDbZp1Qw8lxISbXRIaCUEFAOBTPt+yV48v+lmS9ODQDPVMizO5IjQmggoAwGdk55fq7te+k9OQrumVqhHntjK7JDQyggoAwCccqnDoz/9Zq4OlleqaGqt//L6TbDab2WWhkRFUAACWZxiGJr77ozbuKVSLqDDNuqGnwkODzS4LpwFBBQBgeXO++EXvfb9bwUE2Pf3HHkppGmF2SThNCCoAAEv7att+TfnfT5Kk/zf4bP1f2+YmV4TTiaACALCs3QcP6a4F6+RwGrqy+xm6+fzWZpeE04ygAgCwpLJKh/4yf632l1QoIyVGj17ZmebZAERQAQBYjmEYuv+9LP24s0DNmoTquRt6KiKM5tlARFABAFjO/G9+01trdyrIJs28vodaxjUxuySYhKACALCU1b/ka/IHGyRJEwal64L2LUyuCGYiqAAALCO3sEx3vrpOVU5DQ7ok6/YL25pdEkxGUAEAWEJ5lUN3zF+rvUXlSk+K1rSru9A8C4WYXQAAIHBVOpz6dke+lmzI0ZKNudpTUKaY8BA9f2NPNQnjIwoEFQDAaVZaUaXPNu/Tkg05WvZzngoOVbr3RdtD9PQfeyiteaSJFcJKCCoAgEZ3oKRCn/yUqyUbc/X5lr0qq3S698VFhmnA2Qm6rGOSLmjfgmv4wANBBQDQKHYeKNXSjblavCFHq385IIfTcO9LbRahyzomaWBGonqmNVNIMC2TqBtBBQDQIAzD0ObcYi3ekKMlG3OUtavQY396UrQGZiTpsoxEdUyOoVEW9UJQAQCcNIfT0He/HdCSwyMnv+4vde+z2aRz0uJ0WUaiLuuYpFbNWbQN3iOoAAC8Ul7l0Jfb9mvJhhwt3ZinfcXl7n1hIUG68MwWuiwjUf3PTlSLKLuJlcIfEFQAACdUVFap5Zv2asmGHK3YtFfF5VXufdHhIbokPUEDM5J0UYd4Rdn5aEHD4bcJAFCnvKIyfbIxT4s35OjLbftU6TjSDJsQbXdP6fxf2+YKC6EZFo2DoAIAcNuxr8S9+Nq63w7IOJJN1LZFpC7LqD5Tp2tqUwUF0QyLxkdQAYAAZhiGsnYVus/U2Zxb7LG/a2qsO5ycmRBtUpUIZAQVAAgwVa5l6zfmasmGHO0uKHPvCwmy6f/aNtdlGYm6tGOikmMjTKwUIKgAQEA4VOHQZ1v2avGGHH36c54Olh5Ztj4iNFj9zorXZRmJuuSsRMU2CTWxUsATQQUA/NTB0got+6m6Gfazo5atb9YkVAPOTtTADJath7URVADAj+w6eEhLDzfDfrMj32PZ+jOaRrhXhu3FsvXwEQQVAPAxTqehgkOVyi+t0IGSCuWXVGhTTpGWbMzV+l0FHsemJ0W7m2FZth6+iKACACYyDEPF5VU6UOIZPA6Uev6b795eqYOlFaoxUOKh5rL1l3ZMVFrzyNP7DQENjKACAA2orNJxJGQcDh/5xeXKL62sDiFHhZEDJZWqcDhP/MR1iA4PUfPIMDWLDFNSTLj6dojXgI4sWw//Yomg8swzz+iJJ55QTk6OunbtqpkzZ6p3795mlwUgwFU5nDpQWlkjeNQMGpV1jHZUqLTCcVKvFREarLjIMDWLDFWzJmHVXzcJcwcR1/2ax4TSY4IAYHpQeeONNzR27Fg999xzOvfcczVjxgwNHDhQmzZtUkJCgtnlAbAIp9NQpdOpKoehKseRrysdTlU5DVU5nKp0GKpyHv738PZKx5H7lYePq/n4Q5WOOqZbKpVfUqGCQ5UnLqwOIUG26nDRpDpUNI+0q1lk6OH7R4eO6uMiwjjrBqiLzTCMY8x0nh7nnnuuzjnnHD399NOSJKfTqZYtW+ruu+/WhAkTjvvYwsJCxcbGqqCgQDExMQ1WU2lFlfYXVzTY8/m6hvgNMXRqT3J0DYbHPuMY2+t+RM3txz7es2aPxxgnPsbjeep4XqdRXbdxuH7D0OGvJaf7fvVG5+GvXcc4D3/h2ub5XEd9XeO1VOv5j/FcqlmP4fGcNZ/L4TTcQaDmB3/tQHB0iDg6SHiGiqqjwobr8cfqyWhsNpvUNCK0RvCo8e/hkY3mUZ7BI9oeQtMqcBzefH6bOqJSUVGhtWvXauLEie5tQUFBGjBggL766qtax5eXl6u8/MjlxAsLCxulrk9+ytPo175rlOcG0HBCgmwKCbYpNChIIcE2hQQHKTSo+t+6t9sUGhx0+HFBCg22KTwk+MjIRo2RDlcIiY0I5TRewESmBpV9+/bJ4XAoMTHRY3tiYqJ+/vnnWsdPnTpVkydPbvS6gm02hYf6z/+YbPL+Lztv/xg8mb8dvf2Ls9bRtrr31XxeW32O8dh+7Bc59nPV3H7i412vb7NJQYf/tdXY5v5aR/a7HhNk0+FjDh97+DjX4448X/WGII9jq/fL4/lrPkfN16/ruY56raAjQeDoD/8Q9/bqr0MPh4WQoMPHHrU9tMZjj7ffI3wE2Ri1AAKA6T0q3pg4caLGjh3rvl9YWKiWLVs2+OsM7pKswV2SG/x5AQCAd0wNKi1atFBwcLByc3M9tufm5iopKanW8Xa7XXY7p90BABAoTJ3fCAsLU8+ePbVs2TL3NqfTqWXLlum8884zsTIAAGAFpk/9jB07VpmZmerVq5d69+6tGTNmqKSkRLfccovZpQEAAJOZHlSuvfZa7d27Vw888IBycnLUrVs3LVq0qFaDLQAACDymr6NyKhprHRUAANB4vPn89p9zcAEAgN8hqAAAAMsiqAAAAMsiqAAAAMsiqAAAAMsiqAAAAMsiqAAAAMsiqAAAAMsiqAAAAMsyfQn9U+FaVLewsNDkSgAAQH25Prfrszi+TweVoqIiSVLLli1NrgQAAHirqKhIsbGxxz3Gp6/143Q6tXv3bkVHR8tms5ldjiUVFhaqZcuWys7O5npIFsD7YS28H9bC+2E9jfWeGIahoqIipaSkKCjo+F0oPj2iEhQUpNTUVLPL8AkxMTH8h28hvB/WwvthLbwf1tMY78mJRlJcaKYFAACWRVABAACWRVDxc3a7XQ8++KDsdrvZpUC8H1bD+2EtvB/WY4X3xKebaQEAgH9jRAUAAFgWQQUAAFgWQQUAAFgWQQUAAFgWQcUPTZ06Veecc46io6OVkJCgYcOGadOmTWaXhcMee+wx2Ww2jRkzxuxSAtquXbt0ww03qHnz5oqIiFDnzp21Zs0as8sKSA6HQ/fff7/atGmjiIgItWvXTg8//HC9rgODU/fZZ59p6NChSklJkc1m03vvveex3zAMPfDAA0pOTlZERIQGDBigLVu2nLb6CCp+aOXKlRo1apS+/vprLV26VJWVlbrssstUUlJidmkBb/Xq1Xr++efVpUsXs0sJaAcOHFCfPn0UGhqqjz/+WBs3btSTTz6pZs2amV1aQHr88cc1a9YsPf300/rpp5/0+OOPa9q0aZo5c6bZpQWEkpISde3aVc8880yd+6dNm6annnpKzz33nL755htFRkZq4MCBKisrOy31cXpyANi7d68SEhK0cuVKXXTRRWaXE7CKi4vVo0cPPfvss3rkkUfUrVs3zZgxw+yyAtKECRP0xRdf6PPPPze7FEgaMmSIEhMTNXv2bPe2P/zhD4qIiND8+fNNrCzw2Gw2LVy4UMOGDZNUPZqSkpKicePG6b777pMkFRQUKDExUXPnztV1113X6DUxohIACgoKJElxcXEmVxLYRo0apcGDB2vAgAFmlxLwPvjgA/Xq1UvDhw9XQkKCunfvrhdffNHssgLW+eefr2XLlmnz5s2SpB9++EGrVq3SoEGDTK4MO3bsUE5Ojsf/t2JjY3Xuuefqq6++Oi01+PRFCXFiTqdTY8aMUZ8+fdSpUyezywlYr7/+utatW6fVq1ebXQokbd++XbNmzdLYsWP197//XatXr9bo0aMVFhamzMxMs8sLOBMmTFBhYaHS09MVHBwsh8OhKVOmaMSIEWaXFvBycnIkSYmJiR7bExMT3fsaG0HFz40aNUpZWVlatWqV2aUErOzsbN1zzz1aunSpwsPDzS4Hqg7wvXr10qOPPipJ6t69u7KysvTcc88RVEzw5ptv6tVXX9WCBQuUkZGh77//XmPGjFFKSgrvB5j68Wd33XWXPvzwQy1fvlypqalmlxOw1q5dq7y8PPXo0UMhISEKCQnRypUr9dRTTykkJEQOh8PsEgNOcnKyOnbs6LHt7LPP1m+//WZSRYFt/PjxmjBhgq677jp17txZN954o+69915NnTrV7NICXlJSkiQpNzfXY3tubq57X2MjqPghwzB01113aeHChfr000/Vpk0bs0sKaP3799f69ev1/fffu2+9evXSiBEj9P333ys4ONjsEgNOnz59ap2yv3nzZqWlpZlUUWArLS1VUJDnx1FwcLCcTqdJFcGlTZs2SkpK0rJly9zbCgsL9c033+i88847LTUw9eOHRo0apQULFuj9999XdHS0ex4xNjZWERERJlcXeKKjo2v1B0VGRqp58+b0DZnk3nvv1fnnn69HH31U11xzjb799lu98MILeuGFF8wuLSANHTpUU6ZMUatWrZSRkaHvvvtO06dP16233mp2aQGhuLhYW7dudd/fsWOHvv/+e8XFxalVq1YaM2aMHnnkEbVv315t2rTR/fffr5SUFPeZQY3OgN+RVOdtzpw5ZpeGw/r27Wvcc889ZpcR0P773/8anTp1Mux2u5Genm688MILZpcUsAoLC4177rnHaNWqlREeHm60bdvWmDRpklFeXm52aQFh+fLldX5mZGZmGoZhGE6n07j//vuNxMREw263G/379zc2bdp02upjHRUAAGBZ9KgAAADLIqgAAADLIqgAAADLIqgAAADLIqgAAADLIqgAAADLIqgAAADLIqgAsLx+/fppzJgxZpcBwAQEFQAAYFkEFQAAYFkEFQCWUlJSoptuuklRUVFKTk7Wk08+aXZJAExEUAFgKePHj9fKlSv1/vvva8mSJVqxYoXWrVtndlkATBJidgEA4FJcXKzZs2dr/vz56t+/vyRp3rx5Sk1NNbkyAGZhRAWAZWzbtk0VFRU699xz3dvi4uJ01llnmVgVADMRVAAAgGURVABYRrt27RQaGqpvvvnGve3AgQPavHmziVUBMBM9KgAsIyoqSiNHjtT48ePVvHlzJSQkaNKkSQoK4m8qIFARVABYyhNPPKHi4mINHTpU0dHRGjdunAoKCswuC4BJbIZhGGYXAQAAUBfGUwEAgGURVAAAgGURVAAAgGURVAAAgGURVAAAgGURVAAAgGURVAAAgGURVAAAgGURVAAAgGURVAAAgGURVAAAgGURVAAAgGX9f9v4snQNRnpkAAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjcAAAHHCAYAAABDUnkqAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuNSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/xnp5ZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABBuElEQVR4nO3deXhU5f3+8XuyTRKyQIBsECDs+y4Q6BdQqRSXim1RaC2oWLqAgLgUrKi01rj8UKoiiK1iUdyoQFVcEFmKLAqCElAWQYiQhLBlkpB15vz+CDNkSAJJSHJmeb+uay6Zc86c+UyGi9w+5/M8x2IYhiEAAAAfEWB2AQAAAHWJcAMAAHwK4QYAAPgUwg0AAPAphBsAAOBTCDcAAMCnEG4AAIBPIdwAAACfQrgBAAA+hXADwKcMHz5cw4cPN7sMSdLixYtlsVj0ww8/1PocFotFjzzySJ3VBPgDwg2ABvfYY49pxYoVtX79nj179Mgjj1xWaPAkq1atIsAAdYhwA6DB1UW4mTNnTqXh5pNPPtEnn3xS++JMsGrVKs2ZM6fSfQUFBXrwwQcbuCLAuwWZXQAA1KWQkBCzS6hToaGhZpcAeB1GbgAPdPToUU2cOFGJiYmyWq1KTk7WH//4RxUXF7uOOXjwoMaMGaOYmBiFh4dr0KBB+uCDD9zOs27dOlksFr399tuaM2eOWrRoocjISP3qV79STk6OioqKNH36dMXGxioiIkK33367ioqK3M5hsVg0ZcoUvf766+rUqZNCQ0PVr18/bdiwwe242267TW3atKnwWR555BFZLBa38+Xn5+vVV1+VxWKRxWLRbbfdJkk6fPiw/vSnP6lTp04KCwtT06ZNNWbMGLcRmsWLF2vMmDGSpCuvvNJ1jnXr1kmqvOfm+PHjmjhxouLi4hQaGqpevXrp1VdfdTvmhx9+kMVi0f/7f/9PixYtUrt27WS1WnXFFVfoyy+/dDv2m2++0W233aa2bdsqNDRU8fHxuuOOO3Ty5MkKn/9SbrvtNs2fP9/1s3E+yv+8yl+ycv489+3bp1tvvVXR0dFq3ry5Zs+eLcMwlJ6erhtvvFFRUVGKj4/X3LlzK7xnUVGRHn74YbVv315Wq1VJSUm6//77K3z3gLdi5AbwMMeOHdOAAQN05swZTZo0SZ07d9bRo0e1bNkynT17ViEhIcrKytLgwYN19uxZTZ06VU2bNtWrr76qn//851q2bJluuukmt3OmpqYqLCxMM2fO1IEDB/Tcc88pODhYAQEBOn36tB555BFt2bJFixcvVnJysh566CG3169fv15vvfWWpk6dKqvVqhdeeEE/+9nP9MUXX6h79+41+nxLlizRnXfeqQEDBmjSpEmSpHbt2kmSvvzyS23atEljx45Vy5Yt9cMPP2jBggUaPny49uzZo/DwcA0dOlRTp07Vs88+qwceeEBdunSRJNd/L1RQUKDhw4frwIEDmjJlipKTk/XOO+/otttu05kzZzRt2jS345cuXarc3Fz9/ve/l8Vi0ZNPPqlf/OIXOnjwoIKDgyVJq1ev1sGDB3X77bcrPj5eu3fv1qJFi7R7925t2bLFLZxcyu9//3sdO3ZMq1ev1pIlS6r9ultuuUVdunTR448/rg8++ECPPvqoYmJi9OKLL+qqq67SE088oddff1333nuvrrjiCg0dOlSS5HA49POf/1wbN27UpEmT1KVLF+3atUvPPPOM9u3bd1mXCwGPYQDwKOPHjzcCAgKML7/8ssI+h8NhGIZhTJ8+3ZBk/O9//3Pty83NNZKTk402bdoYdrvdMAzDWLt2rSHJ6N69u1FcXOw6dty4cYbFYjFGjRrldv6UlBSjdevWbtskGZKMbdu2ubYdPnzYCA0NNW666SbXtgkTJlR4rWEYxsMPP2xc+E9No0aNjAkTJlQ49uzZsxW2bd682ZBk/Pvf/3Zte+eddwxJxtq1ayscP2zYMGPYsGGu5/PmzTMkGa+99pprW3FxsZGSkmJEREQYNpvNMAzDOHTokCHJaNq0qXHq1CnXsStXrjQkGe+9995F63zjjTcMScaGDRtc21555RVDknHo0KEKx5c3efLkCj8jJ0nGww8/7Hru/HlOmjTJta20tNRo2bKlYbFYjMcff9y1/fTp00ZYWJjbz3rJkiVGQECA298dwzCMhQsXGpKMzz///KK1At6Ay1KAB3E4HFqxYoVuuOEG9e/fv8J+54jAqlWrNGDAAP3kJz9x7YuIiNCkSZP0ww8/aM+ePW6vGz9+vGvUQZIGDhwowzB0xx13uB03cOBApaenq7S01G17SkqK+vXr53reqlUr3Xjjjfr4449lt9tr/4EvEBYW5vpzSUmJTp48qfbt26tx48b66quvanXOVatWKT4+XuPGjXNtCw4O1tSpU5WXl6f169e7HX/LLbeoSZMmruf/93//J6nsMmBldRYWFurEiRMaNGiQJNW6zpq68847XX8ODAxU//79ZRiGJk6c6NreuHFjderUya32d955R126dFHnzp114sQJ1+Oqq66SJK1du7ZB6gfqE+EG8CDZ2dmy2WyXvNRz+PBhderUqcJ256WZw4cPu21v1aqV2/Po6GhJUlJSUoXtDodDOTk5bts7dOhQ4b06duyos2fPKjs7+6K11kRBQYEeeughJSUlyWq1qlmzZmrevLnOnDlToabqOnz4sDp06KCAAPd/7qr7s3IGndOnT7u2nTp1StOmTVNcXJzCwsLUvHlzJScnS1Kt66ypyr7T0NBQNWvWrML28rXv379fu3fvVvPmzd0eHTt2lFTWnwR4O3puAD8QGBhYo+2GYdT4ParqM6nJyM5dd92lV155RdOnT1dKSoqio6NlsVg0duxYORyOGtdUG9X5mdx8883atGmT7rvvPvXu3VsRERFyOBz62c9+Zmqd1and4XCoR48eevrppys99sLAC3gjwg3gQZo3b66oqCilpaVd9LjWrVtr7969FbZ/9913rv11af/+/RW27du3T+Hh4WrevLmkshGOM2fOVDjuwpERqeogtGzZMk2YMMFthk9hYWGF89akYbd169b65ptv5HA43EZvavuzOn36tNasWaM5c+a4NV5X9jOqrpp8nsvVrl07ff3117r66qsb9H2BhsRlKcCDBAQEaPTo0Xrvvfe0bdu2Cvud/wd+7bXX6osvvtDmzZtd+/Lz87Vo0SK1adNGXbt2rdO6Nm/e7NZLkp6erpUrV+qaa65xjRa0a9dOOTk5+uabb1zHZWRkaPny5RXO16hRo0qDUGBgYIVRo+eee67C6E+jRo0kqdJzXOjaa69VZmam3nrrLde20tJSPffcc4qIiNCwYcMueY4La5Qqjm7NmzevRucpryaf53LdfPPNOnr0qF566aUK+woKCpSfn1/vNQD1jZEbwMM89thj+uSTTzRs2DDXVN2MjAy988472rhxoxo3bqyZM2fqjTfe0KhRozR16lTFxMTo1Vdf1aFDh/Sf//ynQn/J5erevbtGjhzpNhVcktuqumPHjtWf//xn3XTTTZo6darOnj2rBQsWqGPHjhWabPv166dPP/1UTz/9tBITE5WcnKyBAwfq+uuv15IlSxQdHa2uXbtq8+bN+vTTT9W0aVO31/fu3VuBgYF64oknlJOTI6vVqquuukqxsbEVap80aZJefPFF3Xbbbdq+fbvatGmjZcuW6fPPP9e8efMUGRlZo59FVFSUhg4dqieffFIlJSVq0aKFPvnkEx06dKhG57nw5yFJU6dO1ciRIxUYGKixY8fW+nwX89vf/lZvv/22/vCHP2jt2rUaMmSI7Ha7vvvuO7399tv6+OOPK21mB7wJ4QbwMC1atNDWrVs1e/Zsvf7667LZbGrRooVGjRql8PBwSVJcXJw2bdqkP//5z3ruuedUWFionj176r333tN1111X5zUNGzZMKSkpmjNnjo4cOaKuXbtq8eLF6tmzp+uYpk2bavny5ZoxY4buv/9+JScnKzU1Vfv3768Qbp5++mlNmjRJDz74oAoKCjRhwgQNHDhQ//jHPxQYGKjXX39dhYWFGjJkiD799FONHDnS7fXx8fFauHChUlNTNXHiRNntdq1du7bScBMWFqZ169Zp5syZevXVV2Wz2dSpUye98sorrsUDa2rp0qW66667NH/+fBmGoWuuuUYffvihEhMTa3W+X/ziF7rrrrv05ptv6rXXXpNhGPUWbgICArRixQo988wz+ve//63ly5crPDxcbdu21bRp01yNxYA3sxi16RwE4DcsFosmT56s559/3uxSAKBa6LkBAAA+hXADAAB8CuEGAAD4FBqKAVwUbXkAvA0jNwAAwKcQbgAAgE/xu8tSDodDx44dU2RkJEuPAwDgJQzDUG5urhITEy+5UKnfhZtjx45xYzgAALxUenq6WrZsedFj/C7cOJdaT09PV1RUlMnVAACA6rDZbEpKSqrWLVP8Ltw4L0VFRUURbgAA8DLVaSmhoRgAAPgUwg0AAPAphBsAAOBTCDcAAMCnEG4AAIBPIdwAAACfQrgBAAA+hXADAAB8CuEGAAD4FMINAADwKYQbAADgUwg3AADApxBuAABAnSgssWv3sRwVlzpMrYNwAwAA6sTO9DO67tmNGjlvg6l1EG4AAECdSDuaI0nqGBdhah2EGwAAUCd2nQs3PVpEm1oH4QYAANQJ58hNN38ONwsWLFDPnj0VFRWlqKgopaSk6MMPP7zoa9555x117txZoaGh6tGjh1atWtVA1QIAgKrkFZXq4Il8SX4+ctOyZUs9/vjj2r59u7Zt26arrrpKN954o3bv3l3p8Zs2bdK4ceM0ceJE7dixQ6NHj9bo0aOVlpbWwJUDAIDy9hyzyTCkhOhQNYuwmlqLxTAMw9QKLhATE6OnnnpKEydOrLDvlltuUX5+vt5//33XtkGDBql3795auHBhtc5vs9kUHR2tnJwcRUVF1VndAAD4s5c3HtJf39+jEV3i9M8J/ev8/DX5/e0xPTd2u11vvvmm8vPzlZKSUukxmzdv1ogRI9y2jRw5Ups3b26IEgEAQBXSPKSZWJKCzC5g165dSklJUWFhoSIiIrR8+XJ17dq10mMzMzMVFxfnti0uLk6ZmZlVnr+oqEhFRUWu5zabrW4KBwAALq6ZUi3Nvypi+shNp06dtHPnTm3dulV//OMfNWHCBO3Zs6fOzp+amqro6GjXIykpqc7ODQAApLPFpfo+O0+S1D3R/JEb08NNSEiI2rdvr379+ik1NVW9evXSP/7xj0qPjY+PV1ZWltu2rKwsxcfHV3n+WbNmKScnx/VIT0+v0/oBAPB332bY5DCk2EirYqNCzS7H/HBzIYfD4XYZqbyUlBStWbPGbdvq1aur7NGRJKvV6ppq7nwAAIC6s+tHz+m3kUzuuZk1a5ZGjRqlVq1aKTc3V0uXLtW6dev08ccfS5LGjx+vFi1aKDU1VZI0bdo0DRs2THPnztV1112nN998U9u2bdOiRYvM/BgAAPi1tGNl/axmL97nZGq4OX78uMaPH6+MjAxFR0erZ8+e+vjjj/XTn/5UknTkyBEFBJwfXBo8eLCWLl2qBx98UA888IA6dOigFStWqHv37mZ9BAAA/J4nzZSSPHCdm/rGOjcAANSdwhK7uj38sewOQ1tmXa346PrpufHKdW4AAID3+TbDJrvDULOIEMVFmbsysRPhBgAA1JrzklT3FtGyWCwmV1OGcAMAAGrNuXifJ6xv40S4AQAAtZZ2tGymVHcPaSaWCDcAAKCWCkvs2peVK0nq0ZJwAwAAvNy+rFyVOgw1CQ9WYj3NkqoNwg0AAKiVXR7YTCwRbgAAQC152uJ9ToQbAABQK57YTCwRbgAAQC0Ulzq0N/NcMzHhBgAAeLt9WbkqtjsUHRaslk3CzC7HDeEGAADU2PmViaM8qplYItwAAIBaKD9TytMQbgAAQI156kwpiXADAABqqMTu0Lfnmok96Z5SToQbAABQI/uz8lRc6lBkaJBaNw03u5wKCDcAAKBG0srdCdzTmoklwg0AAKihtGPnZ0p5IsINAACoEU+eKSURbgAAQA2U2h36NqPstgueOFNKItwAAIAa+D47X4UlDkVYg9SmaSOzy6kU4QYAAFSb85JU18QoBQR4XjOxRLgBAAA14MmL9zkRbgAAQLWVv6eUpyLcAACAarE7DO0+5tnNxBLhBgAAVNPB7DwVlNgVHhKo5GYRZpdTJcINAACoFufifV0TohTooc3EEuEGAABU064fyy5JeerifU6EGwAAUC3eMFNKItwAAIBqcDgM7T7m2bddcCLcAACASzp0Ml/5xXaFBgeoXXPPXJnYiXADAAAuyXlJqktClIICPTs+eHZ1AADAI3hLv41EuAEAANWw66h39NtIhBsAAHAJDoeh3UfPTQNPJNwAAAAvd+TUWeUWlSokKEAd4jx3ZWInwg0AALioXeWaiYM9vJlYItwAAIBLcN52oXui594JvDzCDQAAuChvmiklEW4AAMBFGIahtKPecU8pJ8INAACo0o+nC5RTUKKQwAB1jIs0u5xqIdwAAIAqOZuJO8VHKiTIO2KDd1QJAABM4U2L9zkRbgAAQJXSXOHGO2ZKSSaHm9TUVF1xxRWKjIxUbGysRo8erb179170NYsXL5bFYnF7hIaGNlDFAAD4j7JmYu+aKSWZHG7Wr1+vyZMna8uWLVq9erVKSkp0zTXXKD8//6Kvi4qKUkZGhutx+PDhBqoYAAD/cfRMgU6fLVFQgEWd4r2jmViSgsx8848++sjt+eLFixUbG6vt27dr6NChVb7OYrEoPj6+vssDAMCvOaeAd4yLlDUo0ORqqs+jem5ycsqGvmJiYi56XF5enlq3bq2kpCTdeOON2r17d5XHFhUVyWazuT0AAMCleeMlKcmDwo3D4dD06dM1ZMgQde/evcrjOnXqpJdfflkrV67Ua6+9JofDocGDB+vHH3+s9PjU1FRFR0e7HklJSfX1EQAA8CmumVItvSvcWAzDMMwuQpL++Mc/6sMPP9TGjRvVsmXLar+upKREXbp00bhx4/S3v/2twv6ioiIVFRW5nttsNiUlJSknJ0dRUd7T+Q0AQEMyDEP9H/1UJ/OLtfxPg9WnVRNT67HZbIqOjq7W729Te26cpkyZovfff18bNmyoUbCRpODgYPXp00cHDhyodL/VapXVaq2LMgEA8BuZtkKdzC9WYIBFXRK8azDA1MtShmFoypQpWr58uT777DMlJyfX+Bx2u127du1SQkJCPVQIAIB/2vVj2SWpDrERCg32nmZiyeSRm8mTJ2vp0qVauXKlIiMjlZmZKUmKjo5WWFiYJGn8+PFq0aKFUlNTJUl//etfNWjQILVv315nzpzRU089pcOHD+vOO+807XMAAOBr0o55180yyzM13CxYsECSNHz4cLftr7zyim677TZJ0pEjRxQQcH6A6fTp0/rd736nzMxMNWnSRP369dOmTZvUtWvXhiobAACf560zpSQPaihuKDVpSAIAwF9d8fdPlZ1bpP/8cbD6tTa3mViq2e9vj5kKDgAAPMNxW6Gyc4sUYJG6elkzsUS4AQAAF3Cub9M+NkJhId7VTCwRbgAAwAWct13wxmZiiXADAAAu4FqZOJFwAwAAfIBrppSX3XbBiXADAABcsnOLlGkrlMVLm4klwg0AACgn7VjZqE3bZo3UyOoRd2mqMcINAABwSfvRexfvcyLcAAAAF+fIjbfOlJIINwAAoBxvnwYuEW4AAMA5p/KLdfRMgSSpW6J3NhNLhBsAAHCOcwp4crNGigwNNrma2iPcAAAASeUW7/PiS1IS4QYAAJzjWryvhfdekpIINwAA4BzXTCkvve2CE+EGAADozNlipZ8610zMZSkAAODtnFPAWzcNV3SY9zYTS4QbAAAg37kkJRFuAACAfGemlES4AQAAKj9TinADAAC8XE5BiQ6fPCvJu1cmdiLcAADg53af67dp2SRMTRqFmFzN5SPcAADg53zpkpREuAEAwO/5wp3AyyPcAADg59J8aKaURLgBAMCv5RaW6OCJfElSdx9oJpYINwAA+LU9x8ouSSVGh6pphNXkauoG4QYAAD/mS4v3ORFuAADwY742U0oi3AAA4NfSjvnWTCmJcAMAgN/KLyrV99l5kgg3AADAB3ybYZNhSHFRVjWP9I1mYolwAwCA39rlg/02EuEGAAC/5YszpSTCDQAAfmu387YLiYQbAADg5QqK7dp/PFeS1KMl4QYAAHi5PRk2OQypeaRVcVGhZpdTpwg3AAD4od3HzvXb+Mj9pMoj3AAA4Id2/eibM6Ukwg0AAH7JV2dKSYQbAAD8TmGJXfuP+97KxE6EGwAA/Mx3mbmyOww1bRSihGjfaiaWCDcAAPid8pekLBaLydXUPcINAAB+Zrcr3PjeTCnJ5HCTmpqqK664QpGRkYqNjdXo0aO1d+/eS77unXfeUefOnRUaGqoePXpo1apVDVAtAAC+wVfvKeVkarhZv369Jk+erC1btmj16tUqKSnRNddco/z8/Cpfs2nTJo0bN04TJ07Ujh07NHr0aI0ePVppaWkNWDkAAN6pqNSufVllKxP7YjOxJFkMwzDMLsIpOztbsbGxWr9+vYYOHVrpMbfccovy8/P1/vvvu7YNGjRIvXv31sKFCy/5HjabTdHR0crJyVFUlG8OxwEAUJVdP+bohuc3qnF4sHbM/qnX9NzU5Pe3R/Xc5OSUDZPFxMRUeczmzZs1YsQIt20jR47U5s2bKz2+qKhINpvN7QEAgL8qf0nKW4JNTXlMuHE4HJo+fbqGDBmi7t27V3lcZmam4uLi3LbFxcUpMzOz0uNTU1MVHR3teiQlJdVp3QAAeBNfXrzPyWPCzeTJk5WWlqY333yzTs87a9Ys5eTkuB7p6el1en4AALzJ+XtK+W64CTK7AEmaMmWK3n//fW3YsEEtW7a86LHx8fHKyspy25aVlaX4+PhKj7darbJarXVWKwAA3qq41KHvMsqaiX11ppRk8siNYRiaMmWKli9frs8++0zJycmXfE1KSorWrFnjtm316tVKSUmprzIBAPAJ+7JyVWx3KCo0SEkxYWaXU29MHbmZPHmyli5dqpUrVyoyMtLVNxMdHa2wsLIf+vjx49WiRQulpqZKkqZNm6Zhw4Zp7ty5uu666/Tmm29q27ZtWrRokWmfAwAAb+C6JOXDzcSSySM3CxYsUE5OjoYPH66EhATX46233nIdc+TIEWVkZLieDx48WEuXLtWiRYvUq1cvLVu2TCtWrLhoEzIAAPD9xfucTB25qc4SO+vWrauwbcyYMRozZkw9VAQAgO/adbRsORRfnikledBsKQAAUH9K7A59m0G4AQAAPuLA8TwVlzoUaQ1S65hws8upV4QbAAD8QNq5fptuLaIUEOC7zcQS4QYAAL/gDDe+vHifE+EGAAA/4Jop1ZJwAwAAvJzdYWjPuWbibozcAAAAb/d9dp4KSxxqFBKots0amV1OvSPcAADg43b9eK6ZODHa55uJJcINAAA+L+3Y+ZlS/oBwAwCAj0vzk9suOBFuAADwYXaHod3HypqJCTcAAMDrHTqRr7PFdoUFB6pt8wizy2kQhBsAAHyY85JU18QoBfpBM7FEuAEAwKft8rN+G4lwAwCAT3PdUyrRP2ZKSXUQbux2u3bu3KnTp0/XRT0AAKCOOMo3E/vBbRecahxupk+frn/961+SyoLNsGHD1LdvXyUlJWndunV1XR8AAKilH07mK6+oVNagALX3k2ZiqRbhZtmyZerVq5ck6b333tOhQ4f03Xff6e6779Zf/vKXOi8QAADUTtq5UZsuCVEKCvSfTpQaf9ITJ04oPj5ekrRq1SqNGTNGHTt21B133KFdu3bVeYEAAKB2/G3xPqcah5u4uDjt2bNHdrtdH330kX76059Kks6ePavAwMA6LxAAANSO855S/hZugmr6gttvv10333yzEhISZLFYNGLECEnS1q1b1blz5zovEAAA1JxhGH53TymnGoebRx55RN27d1d6errGjBkjq9UqSQoMDNTMmTPrvEAAAFBzR06dVW5hqUKCAtQxLtLschpUjcONJP3qV7+qsG3ChAmXXQwAAKgbzsX7usRHKtiPmoklFvEDAMAnpR0tmynVzc/6bSTCDQAAPslfZ0pJhBsAAHyOYRh+eU8pJ8INAAA+5sfTBcopKFFwoEUd4vxnZWKnWjUUOxwOHThwQMePH5fD4XDbN3To0DopDAAA1I7zklSn+EhZg/xvDboah5stW7bo17/+tQ4fPizDMNz2WSwW2e32OisOAADUnHN9G3+8JCXVItz84Q9/UP/+/fXBBx+4FvIDAACeY5dzplQi4aZa9u/fr2XLlql9+/b1UQ8AALgMhmH49UwpqRYNxQMHDtSBAwfqoxYAAHCZMnIKdSq/WEEBFnWK96+ViZ1qPHJz11136Z577lFmZqZ69Oih4OBgt/09e/ass+IAAEDNOKeAd4iLVGiw/zUTS7UIN7/85S8lSXfccYdrm8VikWEYNBQDAGCy85ek/OtmmeXVONwcOnSoPuoAAAB1wBluuvtpv41Ui3DTunXr+qgDAABcprKVictmShFuLuG///2vRo0apeDgYP33v/+96LE///nP66QwAABQM1m2Ip3IK1JggEVdE7gsdVGjR49WZmamYmNjNXr06CqPo+cGAADzOC9JtW8e4bfNxFI1w035WyxceLsFAADgGXbRbyOJG2cCAOAzmClVhnADAICPcN5TipEbAADg9Y7nFirLVqQAi9Q1kZEbAADg5ZyXpNo1j1B4SI1XevEppoabDRs26IYbblBiYqIsFotWrFhx0ePXrVsni8VS4ZGZmdkwBQMA4KHSWN/GpVbh5vvvv9eDDz6ocePG6fjx45KkDz/8ULt3767RefLz89WrVy/Nnz+/Rq/bu3evMjIyXI/Y2NgavR4AAF/DTKnzahxu1q9frx49emjr1q169913lZeXJ0n6+uuv9fDDD9foXKNGjdKjjz6qm266qUavi42NVXx8vOsREMDVNQCAfzs/U4pwU+NUMHPmTD366KNavXq1QkJCXNuvuuoqbdmypU6Lq0rv3r2VkJCgn/70p/r8888b5D0BAPBUJ/KKlJFTKAvNxJJqcW+pXbt2aenSpRW2x8bG6sSJE3VSVFUSEhK0cOFC9e/fX0VFRfrnP/+p4cOHa+vWrerbt2+lrykqKlJRUZHruc1mq9caAQBoaM5Rm+RmjRRh9e9mYqkW4aZx48bKyMhQcnKy2/YdO3aoRYsWdVZYZTp16qROnTq5ng8ePFjff/+9nnnmGS1ZsqTS16SmpmrOnDn1WhcAAGbikpS7Gl+WGjt2rP785z8rMzNTFotFDodDn3/+ue69916NHz++Pmq8qAEDBujAgQNV7p81a5ZycnJcj/T09AasDgCA+ueaKZVIuJFqMXLz2GOPafLkyUpKSpLdblfXrl1lt9v161//Wg8++GB91HhRO3fuVEJCQpX7rVarrFZrA1YEAEDDYqaUuxqHm5CQEL300kuaPXu20tLSlJeXpz59+qhDhw41fvO8vDy3UZdDhw5p586diomJUatWrTRr1iwdPXpU//73vyVJ8+bNU3Jysrp166bCwkL985//1GeffaZPPvmkxu8NAIAvOJ1frKNnCiRJ3fz8nlJOte46atWqlVq1anVZb75t2zZdeeWVruczZsyQJE2YMEGLFy9WRkaGjhw54tpfXFyse+65R0ePHlV4eLh69uypTz/91O0cAAD4E+f9pNo0DVdUaLDJ1XgGi2EYRk1eYBiGli1bprVr1+r48eNyOBxu+9999906LbCu2Ww2RUdHKycnR1FRJFwAgHd7Yd0BPfnRXl3fM0HP/7rymcO+oCa/v2s8cjN9+nS9+OKLuvLKKxUXFyeLxVLrQgEAwOVhplRFNQ43S5Ys0bvvvqtrr722PuoBAAA1wD2lKqrxVPDo6Gi1bdu2PmoBAAA1kHO2REdOnZXENPDyahxuHnnkEc2ZM0cFBQX1UQ8AAKim3eeaiVvFhCs6nGZipxpflrr55pv1xhtvKDY2Vm3atFFwsPsP86uvvqqz4gAAQNXOr2/DBJnyahxuJkyYoO3bt+vWW2+loRgAABOxeF/lahxuPvjgA3388cf6yU9+Uh/1AACAatp9jNsuVKbGPTdJSUmsDwMAgMlshSU6dCJfEtPAL1TjcDN37lzdf//9+uGHH+qhHAAAUB27z00Bb9E4TE0ahZhcjWep8WWpW2+9VWfPnlW7du0UHh5eoaH41KlTdVYcAAConHOmFM3EFdU43MybN68eygAAADWxi5WJq1Sr2VIAAMBczJSqWrXCjc1mczUR22y2ix5LszEAAPUrr6jU1UxMuKmoWuGmSZMmysjIUGxsrBo3blzp2jaGYchischut9d5kQAA4Lw9x2wyDCkhOlTNIqxml+NxqhVuPvvsM8XExEiS1q5dW68FAQCAi+OS1MVVK9wMGzbM9efk5GQlJSVVGL0xDEPp6el1Wx0AAKhgtzPcsHhfpWq8zk1ycrKys7MrbD916pSSk5PrpCgAAFA110yplvS5VqbG4cbZW3OhvLw8hYaG1klRAACgcmeLS/V9dp4kLktVpdpTwWfMmCFJslgsmj17tsLDw1377Ha7tm7dqt69e9d5gQAA4LxvM2xyGFJspFWxkQwqVKba4WbHjh2SykZudu3apZCQ80s9h4SEqFevXrr33nvrvkIAAOCy60cW77uUaocb5yyp22+/Xf/4xz9YzwYAABPsOndPKS5JVa3GKxS/8sor9VEHAACohvP3lCLcVKXGDcUAAMAchSV27T9e1kzMZamqEW4AAPASezJssjsMNYuwKi6KlYmrQrgBAMBLuBbvaxFV6bIsKEO4AQDAS7gW7+OS1EURbgAA8BLMlKoewg0AAF6gsMSu/Vm5kgg3l0K4AQDAC+zNzFWpw1BMoxAlRrMy8cUQbgAA8AJp5da3oZn44gg3AAB4gTTnTKlE7hBwKYQbAAC8ADOlqo9wAwCAhysudWhvJs3E1UW4AQDAw+3LylWJ3VB0WLBaNgkzuxyPR7gBAMDDlb8kRTPxpRFuAADwcM5m4m4taCauDsINAAAeLo1m4hoh3AAA4MFK7A59e66ZmHBTPYQbAAA82P6sPBWXOhQZGqRWMeFml+MVCDcAAHiw84v30UxcXYQbAAA8mGumVEsuSVUX4QYAAA/mvKdUN267UG2EGwAAPFSp3aFvM2ySaCauCcINAAAe6kB2ngpLHIqwBqlN00Zml+M1TA03GzZs0A033KDExERZLBatWLHikq9Zt26d+vbtK6vVqvbt22vx4sX1XicAAGZIO1o2atM1MUoBATQTV5ep4SY/P1+9evXS/Pnzq3X8oUOHdN111+nKK6/Uzp07NX36dN155536+OOP67lSAAAaHov31U6QmW8+atQojRo1qtrHL1y4UMnJyZo7d64kqUuXLtq4caOeeeYZjRw5sr7KBADAFLsIN7XiVT03mzdv1ogRI9y2jRw5Ups3b67yNUVFRbLZbG4PAAA8nd1haM+xst9Z3bmnVI14VbjJzMxUXFyc27a4uDjZbDYVFBRU+prU1FRFR0e7HklJSQ1RKgAAl+Vgdp4KSuwKDwlUcrMIs8vxKl4Vbmpj1qxZysnJcT3S09PNLgkAgEtyXpLqlhilQJqJa8TUnpuaio+PV1ZWltu2rKwsRUVFKSwsrNLXWK1WWa3WhigPAIA645wp1S2Rfpua8qqRm5SUFK1Zs8Zt2+rVq5WSkmJSRQAA1A9mStWeqeEmLy9PO3fu1M6dOyWVTfXeuXOnjhw5IqnsktL48eNdx//hD3/QwYMHdf/99+u7777TCy+8oLffflt33323GeUDAFAvHA5Du49xT6naMjXcbNu2TX369FGfPn0kSTNmzFCfPn300EMPSZIyMjJcQUeSkpOT9cEHH2j16tXq1auX5s6dq3/+859MAwcA+JRDJ/OVX2xXaHCA2jZjZeKaMrXnZvjw4TIMo8r9la0+PHz4cO3YsaMeqwIAwFzOS1JdE6IUFOhVHSQegZ8YAAAehn6by0O4AQDAw7imgRNuaoVwAwCAB3E4DO0+Nw2ckZvaIdwAAOBBjpw6q9yiUlmDAtQhlpWJa4NwAwCAB3FekupMM3Gt8VMDAMCDnG8m5maZtUW4AQDAg6QdY6bU5SLcAADgIQzD4J5SdYBwAwCAh0g/VaCcghKFBAaoY1yk2eV4LcINAAAewnlJqlN8pEKC+BVdW/zkAADwEM6ZUt3pt7kshBsAADwEt12oG4QbAAA8QFkzsXPkhmngl4NwAwCABzh6pkCnz5YoONCiTvE0E18Owg0AAB7AOWrTMS5S1qBAk6vxboQbAAA8gHN9m+6sb3PZCDcAAHgA10yploSby0W4AQDAZOWbiZkpdfkINwAAmCzTVqiT+cUKDLCoM83El41wAwCAybb9cFqS1CE2QqHBNBNfLsINAAAmyswp1F/f3yNJGtK+mcnV+AbCDQAAJikqtesPr21Xdm6ROsVFasZPO5pdkk8g3AAAYALDMDR7RZp2pp9RdFiwFo3vp0bWILPL8gmEGwAATPDalsN6e9uPCrBIz47ro9ZNG5ldks8g3AAA0MC+OHRKc94r67O5/2edNaxjc5Mr8i2EGwAAGtCxMwX60+vbVeowdH3PBP1+aFuzS/I5hBsAABpIYYldf3xtu07kFatzfKSe/FVPWSwWs8vyOYQbAAAagGEY+svyNH39Y44ahwfrpfH9FR5CA3F9INwAANAAXt30g/7zVVkD8fPj+iopJtzsknwW4QYAgHq2+fuT+tsH30qSHri2i37SgcX66hPhBgCAenT0TIEmL/1Kdoeh0b0TNfEnyWaX5PMINwAA1JPCErt+v2SbTuUXq1tilFJ/QQNxQyDcAABQDwzD0Kx3dyntqE0xjUL04m/7KSyEm2I2BMINAAD14OXPf9DyHUcVGGDR87/uo5ZNaCBuKIQbAADq2KYDJ/TYqrIG4r9c20WD29FA3JAINwAA1KH0U2ddDcS/6NtCtw9pY3ZJfodwAwBAHSkotuv3S7br9NkS9WgRrcdu6kEDsQkINwAA1AHDMPTn/3yjPRk2NT3XQBwaTAOxGQg3AADUgZf+d1D//fqYggIseuE3fZXYOMzskvwW4QYAgMv0v/3ZevzD7yRJD93QVQPbNjW5Iv9GuAEA4DIcOXlWU5bukMOQxvRrqd8Oam12SX6PcAMAQC2dLS7VpCXblFNQol5JjfW30d1pIPYAhBsAAGrBMAzdt+wbfZeZq2YRVr14Kw3EnoJwAwBALSxcf1AffJOh4ECLFt7aV/HRoWaXhHM8ItzMnz9fbdq0UWhoqAYOHKgvvviiymMXL14si8Xi9ggN5S8UAKDhrNt7XE9+XNZA/PAN3dS/TYzJFaE808PNW2+9pRkzZujhhx/WV199pV69emnkyJE6fvx4la+JiopSRkaG63H48OEGrBgA4M9+OJGvqW/skGFI4wYk6TcDW5ldEi5gerh5+umn9bvf/U633367unbtqoULFyo8PFwvv/xyla+xWCyKj493PeLi4hqwYgCAv8orKmsgthWWqm+rxnrk591oIPZApoab4uJibd++XSNGjHBtCwgI0IgRI7R58+YqX5eXl6fWrVsrKSlJN954o3bv3l3lsUVFRbLZbG4PAABqyjAM3fv219qXlafYSKsW3NpP1iAaiD2RqeHmxIkTstvtFUZe4uLilJmZWelrOnXqpJdfflkrV67Ua6+9JofDocGDB+vHH3+s9PjU1FRFR0e7HklJSXX+OQAAvm/+2gP6aHemggMtWnBrP8VF0e/pqUy/LFVTKSkpGj9+vHr37q1hw4bp3XffVfPmzfXiiy9WevysWbOUk5PjeqSnpzdwxQAAb/fZd1mau3qfJOlvN3ZXv9ZNTK4IFxNk5ps3a9ZMgYGBysrKctuelZWl+Pj4ap0jODhYffr00YEDByrdb7VaZbVaL7tWAIB/Opidp2lv7JRhSL8Z2EpjB9BA7OlMHbkJCQlRv379tGbNGtc2h8OhNWvWKCUlpVrnsNvt2rVrlxISEuqrTACAn8otLNGkJduVW1Sq/q2b6OEbupldEqrB1JEbSZoxY4YmTJig/v37a8CAAZo3b57y8/N1++23S5LGjx+vFi1aKDU1VZL017/+VYMGDVL79u115swZPfXUUzp8+LDuvPNOMz8GAMDHOByGZrz9tQ4cz1N8VKheuLWvQoK8rpvDL5kebm655RZlZ2froYceUmZmpnr37q2PPvrI1WR85MgRBQSc/8t0+vRp/e53v1NmZqaaNGmifv36adOmTeratatZHwEA4IOe++yAVu/JUkhggBbc2lexkTQQewuLYRiG2UU0JJvNpujoaOXk5CgqKsrscgAAHmj1niz97t/bJElP/qqnbu7PTFuz1eT3N+NrAACUc+B4nu5+a6ckaXxKa4KNFyLcAABwjq2wRJOWbFNeUakGJMdo9vW0PHgjwg0AACprIL77zZ06mJ2vhOhQvfCbvgoO5NekN+JbAwBA0rw1+7Xmu+MKCQrQi7/tp2YRrJHmrQg3AAC/91Fapp5ds1+SlHpTD/Vs2djcgnBZCDcAAL+2PytX97y9U5J0+5A2+mW/luYWhMtGuAEA+K2cgrIViPOL7RrUNkYPXNvF7JJQBwg3AAC/ZHcYmvbmDh06ka8WjcM0/9c0EPsKvkUAgF96evVerdubLeu5BuKmNBD7DMINAMDvrNqVoflrv5ckPfHLnureItrkilCXCDcAAL+yNzNX977ztSTpzp8ka3SfFiZXhLpGuAEA+I0zZ4v1u39v09liu4a0b6qZozqbXRLqAeEGAOAX7A5Dd72xQ0dOnVXLJmF6flxfBdFA7JP4VgEAfuGpj/fqf/tPKDQ4QIt+219NGoWYXRLqCeEGAODz3vv6mBauL2sgfvJXvdQ1McrkilCfCDcAAJ+255hN9y/7RpL0+2Ft9fNeiSZXhPpGuAEA+KzT+cWatGSbCkrs+r8OzXT/SBqI/QHhBgDgk0rtDt31xg79eLpArWLC9dy4PgoMsJhdFhoA4QYA4JOe+Og7bTxwQuEhgVo0vp8ah9NA7C8INwAAn7Ny51G99L9DkqT/N6aXOsfTQOxPCDcAAJ+SdjTH1UD8p+HtdG2PBJMrQkMj3AAAfMbJvCL9fsl2FZU6NLxTc91zTSezS4IJCDcAAJ9QandoytIdOnqmQG2ahusfY2kg9leEGwCAT3hs1XfafPCkGoUEatH4/ooOCza7JJiEcAMA8HrvfvWjXv68rIF47s291TEu0uSKYCbCDQDAq+36MUez3t0lSZp6VXv9rHu8yRXBbEFmFwAAQG0Ulzq07fAp3fv21yoqdejqzrGaPqKj2WXBAxBuAABe4/DJfG3Yl631+7K1+fuTyi+2S5LaNm+kZ8b2VgANxBDhBgDgwfKLSrX5+5PasL8s0Bw+edZtf7OIEA3t2Fx3j+ioqFAaiFGGcAMA8BiGYWhPhk0b9p3Qhn3Z2nb4lErshmt/UIBF/Vo30dCOzTWsY3N1TYhitAYVEG4AAKY6mVekjQdOaP2+bP1v/wll5xa57U+KCdOwjs01tENzpbRrqkhGaHAJhBsAQIMqsTu0M/2M1u/N1ob92dp1NEfG+cEZhQUHanC7phrasbmGdmyuNk3DZbEwOoPqI9wAAOpd+qmz2rA/Wxv2ZWvTgZPKLSp12985PlLDOjXXsA7N1a9NE1mDAk2qFL6AcAMAqHMFxXZtOXTSNbPpYHa+2/4m4cH6SYfm5y43NVNsVKhJlcIXEW4AAJfNMAzty8rThn1ll5q2Hjql4lKHa39ggEV9khq7GoG7t4jmvk+oN4QbAECtnDlbrI0HymY1bdh3Qpm2Qrf9LRqHaWjHZhraobkGt2/GvZ7QYAg3AIBqsTsM7Uw/4xqd+Tr9jBzlGoGtQQEa1LbpudGZZmrXPIJGYJiCcAMAqFJmTqGrb2bjgRPKKShx298hNqKsb6Zjcw1IjlFoMI3AMB/hBgDgUlhi15c/nHIFmn1ZeW77o0KD9JMOzTSsY3P9X4fmSmwcZlKlQNUINwDgxwzD0PfZ+a5LTVsOnlRhyflGYItF6tWysetSU6+WjRUUGGBixcClEW4AwA+U2B06nlukzJwCZeQUKjOnUN9n52nDvhM6eqbA7di4KKuGdii71PST9s3UpFGISVUDtUO4AQAvV1hiV2ZOYVlosZWFlyzX87L/nsgrclsFuLyQwAANSI4pm9nUsbk6xUXSCAyvRrgBAA9lGIZyi0rPB5VyASazXHg5c7bk0ieTFBxoUVxUqBKiQxUXFaoWTcI0KLmpBraNUXgIvw7gO/jbDAAmMAxDp/KLy0ZZbO7hpex5WYDJL7ZX63xhwYFKiA5VvPNxLsTER4e5wkzTRiHcQRt+wSPCzfz58/XUU08pMzNTvXr10nPPPacBAwZUefw777yj2bNn64cfflCHDh30xBNP6Nprr23AigGganaHoezcIldAybSVG3VxPrcVuq3gezFRoUFKiA5TfPT5URdnkHFujwoN4lIScI7p4eatt97SjBkztHDhQg0cOFDz5s3TyJEjtXfvXsXGxlY4ftOmTRo3bpxSU1N1/fXXa+nSpRo9erS++uorde/e3YRPAMCX2R2GSuwOlToMldodKrEbKii2n+tlqTy8ZOcVye6oosHlAs0irIqPtio+Kuz8yEu58BIfHcolI6CGLIZRVYtZwxg4cKCuuOIKPf/885Ikh8OhpKQk3XXXXZo5c2aF42+55Rbl5+fr/fffd20bNGiQevfurYULF17y/Ww2m6Kjo5WTk6OoqKg6+xzHcwv1/tcZCrBIAQEWWSRZLBZZLFKApex5gMUinXseYJFrn87ts5TbJzmPOX9s2Z/Pn8vi3F7h2Irvayn3fuX3VXVsec6/IYYM13PDte/ctguOlYxyr6t4HtefKzn3xd7X/XxGpe974XkklftMFrfnzp/B+T3l91kueG3V+5zfw4XnvvC9K+47/94Xe9/Knnsjw5ArKJTYHSq1Gyp1lAUG5/Oq9pe6tp8LGY7y2yoe69pvNy74s/v7OANLaYXtZcfW9l/IwACL4iKtijs32uIWXs4FmLioUIUEMa0aqI6a/P429X8HiouLtX37ds2aNcu1LSAgQCNGjNDmzZsrfc3mzZs1Y8YMt20jR47UihUrKj2+qKhIRUVFruc2m+3yC6/Ej6cL9Nf399TLuQF4FmtQgCugOMNKQlRZf4vz0lGzCCs3hgRMYmq4OXHihOx2u+Li4ty2x8XF6bvvvqv0NZmZmZUen5mZWenxqampmjNnTt0UfBGNw4J1fc8EGSobUXA4ykYcHMa5kY5zowwOo2zUwvlfQ+7H6tw+h+t4Sca585w71nHufyXdzqXz7+Mof/4K71t5LSp3HrvDqHqE4tzntVgsumAwpOJohKXyEYrz57v4SIflIiMdKnfspWqtbOTI+fzC/yuvbCSq6pEj16sqHVWqcJ4Ltqnca6oaAatsRMubWSxScECAggItCgoMUHDAuf8GWhQceG57QNnzoHPHBQcGKCigkv1ufz5/rqBAS5Xvcf6cVZ8/ODCg3OvPn8t5XGCAhd4WwMP5/IXcWbNmuY302Gw2JSUl1fn7tG0eoed/3bfOzwsAAGrG1HDTrFkzBQYGKisry217VlaW4uPjK31NfHx8jY63Wq2yWq11UzAAAPB4pnayhYSEqF+/flqzZo1rm8Ph0Jo1a5SSklLpa1JSUtyOl6TVq1dXeTwAAPAvpl+WmjFjhiZMmKD+/ftrwIABmjdvnvLz83X77bdLksaPH68WLVooNTVVkjRt2jQNGzZMc+fO1XXXXac333xT27Zt06JFi8z8GAAAwEOYHm5uueUWZWdn66GHHlJmZqZ69+6tjz76yNU0fOTIEQUEnB9gGjx4sJYuXaoHH3xQDzzwgDp06KAVK1awxg0AAJDkAevcNLT6WucGAADUn5r8/mb1KAAA4FMINwAAwKcQbgAAgE8h3AAAAJ9CuAEAAD6FcAMAAHwK4QYAAPgUwg0AAPAphBsAAOBTTL/9QkNzLshss9lMrgQAAFSX8/d2dW6s4HfhJjc3V5KUlJRkciUAAKCmcnNzFR0dfdFj/O7eUg6HQ8eOHVNkZKQsFovZ5Xgkm82mpKQkpaenc/8tD8D34Vn4PjwP34lnqa/vwzAM5ebmKjEx0e2G2pXxu5GbgIAAtWzZ0uwyvEJUVBT/UHgQvg/PwvfhefhOPEt9fB+XGrFxoqEYAAD4FMINAADwKYQbVGC1WvXwww/LarWaXQrE9+Fp+D48D9+JZ/GE78PvGooBAIBvY+QGAAD4FMINAADwKYQbAADgUwg3AADApxBu4JKamqorrrhCkZGRio2N1ejRo7V3716zy4Kkxx9/XBaLRdOnTze7FL929OhR3XrrrWratKnCwsLUo0cPbdu2zeyy/JLdbtfs2bOVnJyssLAwtWvXTn/729+qdd8hXL4NGzbohhtuUGJioiwWi1asWOG23zAMPfTQQ0pISFBYWJhGjBih/fv3N1h9hBu4rF+/XpMnT9aWLVu0evVqlZSU6JprrlF+fr7Zpfm1L7/8Ui+++KJ69uxpdil+7fTp0xoyZIiCg4P14Ycfas+ePZo7d66aNGlidml+6YknntCCBQv0/PPP69tvv9UTTzyhJ598Us8995zZpfmF/Px89erVS/Pnz690/5NPPqlnn31WCxcu1NatW9WoUSONHDlShYWFDVIfU8FRpezsbMXGxmr9+vUaOnSo2eX4pby8PPXt21cvvPCCHn30UfXu3Vvz5s0zuyy/NHPmTH3++ef63//+Z3YpkHT99dcrLi5O//rXv1zbfvnLXyosLEyvvfaaiZX5H4vFouXLl2v06NGSykZtEhMTdc899+jee++VJOXk5CguLk6LFy/W2LFj670mRm5QpZycHElSTEyMyZX4r8mTJ+u6667TiBEjzC7F7/33v/9V//79NWbMGMXGxqpPnz566aWXzC7Lbw0ePFhr1qzRvn37JElff/21Nm7cqFGjRplcGQ4dOqTMzEy3f7eio6M1cOBAbd68uUFq8LsbZ6J6HA6Hpk+friFDhqh79+5ml+OX3nzzTX311Vf68ssvzS4Fkg4ePKgFCxZoxowZeuCBB/Tll19q6tSpCgkJ0YQJE8wuz+/MnDlTNptNnTt3VmBgoOx2u/7+97/rN7/5jdml+b3MzExJUlxcnNv2uLg41776RrhBpSZPnqy0tDRt3LjR7FL8Unp6uqZNm6bVq1crNDTU7HKgssDfv39/PfbYY5KkPn36KC0tTQsXLiTcmODtt9/W66+/rqVLl6pbt27auXOnpk+frsTERL4PcFkKFU2ZMkXvv/++1q5dq5YtW5pdjl/avn27jh8/rr59+yooKEhBQUFav369nn32WQUFBclut5tdot9JSEhQ165d3bZ16dJFR44cMaki/3bfffdp5syZGjt2rHr06KHf/va3uvvuu5Wammp2aX4vPj5ekpSVleW2PSsry7WvvhFu4GIYhqZMmaLly5frs88+U3Jystkl+a2rr75au3bt0s6dO12P/v376ze/+Y127typwMBAs0v0O0OGDKmwNMK+ffvUunVrkyryb2fPnlVAgPuvsMDAQDkcDpMqglNycrLi4+O1Zs0a1zabzaatW7cqJSWlQWrgshRcJk+erKVLl2rlypWKjIx0XRuNjo5WWFiYydX5l8jIyAq9To0aNVLTpk3pgTLJ3XffrcGDB+uxxx7TzTffrC+++EKLFi3SokWLzC7NL91www36+9//rlatWqlbt27asWOHnn76ad1xxx1ml+YX8vLydODAAdfzQ4cOaefOnYqJiVGrVq00ffp0Pfroo+rQoYOSk5M1e/ZsJSYmumZU1TsDOEdSpY9XXnnF7NJgGMawYcOMadOmmV2GX3vvvfeM7t27G1ar1ejcubOxaNEis0vyWzabzZg2bZrRqlUrIzQ01Gjbtq3xl7/8xSgqKjK7NL+wdu3aSn9fTJgwwTAMw3A4HMbs2bONuLg4w2q1GldffbWxd+/eBquPdW4AAIBPoecGAAD4FMINAADwKYQbAADgUwg3AADApxBuAACATyHcAAAAn0K4AQAAPoVwA8BnDR8+XNOnTze7DAANjHADAAB8CuEGAAD4FMINAJ+Qn5+v8ePHKyIiQgkJCZo7d67ZJQEwCeEGgE+47777tH79eq1cuVKffPKJ1q1bp6+++srssgCYIMjsAgDgcuXl5elf//qXXnvtNV199dWSpFdffVUtW7Y0uTIAZmDkBoDX+/7771VcXKyBAwe6tsXExKhTp04mVgXALIQbAADgUwg3ALxeu3btFBwcrK1bt7q2nT59Wvv27TOxKgBmoecGgNeLiIjQxIkTdd9996lp06aKjY3VX/7yFwUE8P9vgD8i3ADwCU899ZTy8vJ0ww03KDIyUvfcc49ycnLMLguACSyGYRhmFwEAAFBXGLMFAAA+hXADAAB8CuEGAAD4FMINAADwKYQbAADgUwg3AADApxBuAACATyHcAAAAn0K4AQAAPoVwAwAAfArhBgAA+BTCDQAA8Cn/H5f9/YrcckfQAAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "import numpy as np\n", + "import time\n", + "import matplotlib.pyplot as plt\n", + "\n", + "storage = []\n", + "times = []\n", + "for d in range(1,11):\n", + " T = np.random.random_sample(2*d*[2])\n", + " storage.append(T.size*T.itemsize/1e6)\n", + " T = T.reshape([2**d, 2**d])\n", + " start_time = time.time()\n", + " U, S, Vt = np.linalg.svd(T)\n", + " end_time = time.time() - start_time\n", + " times.append(end_time)\n", + "\n", + "plt.figure()\n", + "plt.title('storage consumption')\n", + "plt.xlabel('d')\n", + "plt.ylabel('memory in MB')\n", + "plt.plot(np.arange(1,11), storage)\n", + "\n", + "plt.figure()\n", + "plt.title('computational time')\n", + "plt.xlabel('d')\n", + "plt.ylabel('time in s')\n", + "plt.plot(np.arange(1,11), times)" + ] + }, + { + "cell_type": "markdown", + "id": "e827c10b-ffbd-4687-8647-ba488b100e61", + "metadata": {}, + "source": [ + "***" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "75fa57e8-6077-4b07-ad00-06251ba6a395", + "metadata": {}, + "source": [ + "## Exercise 1.3 \n", + "\n", + "Lastly, we will create a special form of tensor, sometimes called the delta tensor. We define $\\Delta \\in \\mathbb{R}^{n \\times n \\times n}$ as\n", + "\n", + "$\\hspace{1cm}$$\\Delta_{i,j,k} = \\begin{cases} 1 & \\text{if}~i=j=k\\\\ 0& \\text{otherwise.}\\end{cases}$\n", + "\n", + "Now, we want to try various contractions with this tensor. To do this, let's choose $n=3$ and additionally define the vectors $v=(1,2,3)$ and $w=(4,5,6)$.\n", + "\n", + "Use NumPy's einsum function to compute the following contractions:\n", + "\n", + "\\begin{array}{}\n", + " & & & & & & & v & & & & & & & \\\\\n", + " & & & & & & & | & & & & & | & & \\\\\n", + "- & D & = & D & - & \\qquad\\qquad\\qquad & & D & & \\qquad\\qquad\\qquad & & & D & & \\\\\n", + " & & & & & & \\diagup & & \\diagdown & & & \\diagup & & \\diagdown & \\\\\n", + " & & & & & & & & & & v & & & & w\n", + "\\end{array}\n", + "\n", + "Examine the results and try to explain how they arise.\n", + "\n", + "*Hint:* In order to contract two or more tensors with NumPy's einsum function, assign a variable (```i, j, k, ...```) to each mode of the tensors. Modes that are to be contracted receive the same variables. Example: ```np.einsum('ijk, jl -> ikl', T_1, T_2)``` contracts an order-3 tensor $T_1$ and a matrix $T_2$ on a common mode (second of $T_1$, first of $T_2$). The result is again an order-3 tensor." + ] + }, + { + "cell_type": "markdown", + "id": "e79e857f-b2db-47b6-8744-b42bba78f683", + "metadata": {}, + "source": [ + "***" + ] + }, + { + "cell_type": "markdown", + "id": "cfcbcb79-d6db-4fb0-a2d7-14236ba0ef52", + "metadata": {}, + "source": [ + "$\\textcolor{red}{\\textbf{SOLUTION:}}$" + ] + }, + { + "cell_type": "code", + "execution_count": 22, + "id": "c9a8e8d9-bf67-4e21-965c-fe0a9a224fbe", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "-D=D-\n", + " \n", + "[[1. 0. 0.]\n", + " [0. 1. 0.]\n", + " [0. 0. 1.]]\n", + " \n", + "----\n", + " \n", + " v \n", + " | \n", + "-D-\n", + " \n", + "[[1. 0. 0.]\n", + " [0. 2. 0.]\n", + " [0. 0. 3.]]\n", + " \n", + "----\n", + " \n", + " | \n", + "v-D-w\n", + " \n", + "[ 4. 10. 18.]\n" + ] + } + ], + "source": [ + "import numpy as np\n", + "\n", + "D = np.zeros([3,3,3])\n", + "for i in range(3):\n", + " D[i,i,i] = 1\n", + "\n", + "v = np.array([1,2,3])\n", + "w = np.array([4,5,6])\n", + "\n", + "print('-D=D-')\n", + "print(' ')\n", + "print(np.einsum('ijk,ljk->il', D, D))\n", + "print(' ')\n", + "print('----')\n", + "print(' ')\n", + "\n", + "print(' v ')\n", + "print(' | ')\n", + "print('-D-')\n", + "print(' ')\n", + "print(np.einsum('ijk,k->ij', D, v))\n", + "print(' ')\n", + "print('----')\n", + "print(' ')\n", + "\n", + "print(' | ')\n", + "print('v-D-w')\n", + "print(' ')\n", + "print(np.einsum('ijk,j,k->i', D, v, w))" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "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.10" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/Exercises/Session 2 - Tensor Decomposition (with solutions).ipynb b/Exercises/Session 2 - Tensor Decomposition (with solutions).ipynb new file mode 100644 index 0000000..0157287 --- /dev/null +++ b/Exercises/Session 2 - Tensor Decomposition (with solutions).ipynb @@ -0,0 +1,379 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "e4f44eb7", + "metadata": {}, + "source": [ + "$\\def\\tcoreleft{\\underset{\\tiny\\mid}{\\textcolor{MidnightBlue}{⦸}}}$\n", + "$\\def\\tcorecenter{\\underset{\\tiny\\mid}{\\textcolor{RedOrange}{⦿}}}$\n", + "$\\def\\tcoreright{\\underset{\\tiny\\mid}{\\textcolor{MidnightBlue}{\\oslash}}}$\n", + "

TMQS Workshop 2024 @ Zuse Institute Berlin

\n", + "

Summer School on Tensor Methods for Quantum Simulation

\n", + "

June 3 - 5, 2024

\n", + "

$\\tcoreleft - \\tcoreleft - \\tcoreleft - \\cdots - \\tcorecenter - \\cdots - \\tcoreright - \\tcoreright$

\n", + "
" + ] + }, + { + "cell_type": "markdown", + "id": "58a30939-f570-45cd-a736-d6f21aeb2a0c", + "metadata": {}, + "source": [ + "***" + ] + }, + { + "cell_type": "markdown", + "id": "6c3441f3-b5a9-42c4-ba21-2bc682b0d8ac", + "metadata": {}, + "source": [ + "## **Session 2 - Tensor Decomposition**" + ] + }, + { + "cell_type": "markdown", + "id": "6f6cc702", + "metadata": {}, + "source": [ + "***" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "858a33fe-e16e-4dbb-b369-511d7c71fe71", + "metadata": {}, + "source": [ + "## Exercise 2.1\n", + "\n", + "In general, the question of whether a given tensor can be represented in a low-rank form is very challenging. The following (small!) example shall illustrate this.\n", + "\n", + "**a)**$\\quad$Consider the order-3 tensor $T \\in \\mathbb{R}^{2 \\times 2 \\times 2}$ given by\n", + "\n", + "$\\hspace{1.25cm}$$T_{:, :, 1} = \\begin{pmatrix} 1 & 0 \\\\ 0 & 1\\end{pmatrix}$ and $T_{:, :, 2} = \\begin{pmatrix} 0 & 1 \\\\ 1 & 0\\end{pmatrix}$.\n", + "\n", + "$\\hspace{0.75cm}$Find a canonical decomposition with minimum rank!\n", + "\n", + "$\\hspace{0.75cm}$*Hint:* Both matrices are circulant. That is, they share the same eigenvectors (Fourier modes!).\n", + "\n", + "**b)**$\\quad$*extra task:* Let $T \\in \\mathbb{R}^{n \\times n \\times n}$ with $T_{:,:,k} = S^{k-1}$, where $S = \\begin{pmatrix} 0 & & & 1 \\\\ 1 & 0 & & \\\\ & \\ddots & \\ddots & \\\\ 0 & & 1 & 0 \\end{pmatrix}$. What is the canonical rank of $T$?" + ] + }, + { + "cell_type": "markdown", + "id": "bc6abad4-3b47-447c-9e87-86b44381164b", + "metadata": {}, + "source": [ + "***" + ] + }, + { + "cell_type": "markdown", + "id": "11976a26-d5e7-4ca5-b02b-ce4709527c8c", + "metadata": {}, + "source": [ + "$\\textcolor{red}{\\textbf{SOLUTION:}}$" + ] + }, + { + "cell_type": "markdown", + "id": "742e700e-8485-4923-870b-d34994d2978b", + "metadata": {}, + "source": [ + "**a)**$\\quad$The (normalized) eigenvectors of $T_{:, :, 1}$ and $T_{:, :, 2}$ are given by $v_0 = \\frac{1}{\\sqrt{2}}\\begin{pmatrix}1\\\\1\\end{pmatrix}$ and $v_1 = \\frac{1}{\\sqrt{2}}\\begin{pmatrix}1\\\\-1\\end{pmatrix}$. \n", + "\n", + "$\\hspace{0.75cm}$The eigenvalues of $T_{:, :, 1}$ are $\\lambda_0 = \\lambda_1 =1$ and the eigenvalues of $T_{:, :, 2}$ are $\\mu_0 = 1$ and $\\mu_1 = -1$.\n", + "\n", + "$\\hspace{0.75cm}$That is, $T_{:, :, 1} = \\sum_{i=0}^1 \\lambda_i v_i v_i^\\dagger$ and $T_{:, :, 2} = \\sum_{i=0}^1 \\mu_i v_i v_i^\\dagger$.\n", + "\n", + "$\\hspace{0.75cm}$Since $v_i v_i^\\dagger = v_i \\otimes \\overline{v_i} = v_i \\otimes v_i$, we get\n", + "\n", + "$\\hspace{1.25cm}$$\\displaystyle T = \\begin{pmatrix} 1 & 0 \\\\ 0 & 1\\end{pmatrix} \\otimes \\begin{pmatrix} 1 \\\\ 0 \\end{pmatrix} + \\begin{pmatrix} 0 & 1 \\\\ 1 & 0\\end{pmatrix} \\otimes \\begin{pmatrix}0 \\\\ 1\\end{pmatrix}$\n", + "\n", + "$\\hspace{1.25cm}$$\\displaystyle\\phantom{T} = \\sum_{i=0}^1 \\lambda_i v_i \\otimes v_i \\otimes \\begin{pmatrix} 1 \\\\ 0 \\end{pmatrix} + \\sum_{i=0}^1 \\mu_i v_i \\otimes v_i \\otimes \\begin{pmatrix}0 \\\\ 1\\end{pmatrix}$\n", + "\n", + "$\\hspace{1.25cm}$$\\displaystyle\\phantom{T} = \\sum_{i=0}^1 v_i \\otimes v_i \\otimes \\begin{pmatrix} \\lambda_i \\\\ 0 \\end{pmatrix} + \\sum_{i=0}^1 v_i \\otimes v_i \\otimes \\begin{pmatrix}0 \\\\ \\mu_i\\end{pmatrix}$\n", + "\n", + "$\\hspace{1.25cm}$$\\displaystyle\\phantom{T} = \\sum_{i=0}^1 v_i \\otimes v_i \\otimes \\begin{pmatrix} \\lambda_i \\\\ \\mu_i \\end{pmatrix} $\n", + "\n", + "$\\hspace{1.25cm}$$\\displaystyle\\phantom{T} = \\frac{1}{2} \\begin{pmatrix} 1\\\\1 \\end{pmatrix} \\otimes \\begin{pmatrix} 1\\\\1 \\end{pmatrix} \\otimes \\begin{pmatrix} 1\\\\1 \\end{pmatrix} + \\frac{1}{2} \\begin{pmatrix} 1\\\\-1 \\end{pmatrix} \\otimes \\begin{pmatrix} 1\\\\-1 \\end{pmatrix} \\otimes \\begin{pmatrix} 1\\\\-1 \\end{pmatrix}$\n", + "\n", + "Note that there exists no rank-$1$ decomposition of $T$. This means the canonical rank of $T$ is equal to $2$.\n", + "\n", + "**b)**$\\quad$Again, all \"layers\" $T_{:,:,k}$ have the same set of (complex) eigenvectors $v_0, \\dots, v_{n-1}$ and can be written as \n", + "\n", + "$\\hspace{1.25cm}$$\\displaystyle T_{:,:,k} = \\sum_{i=0}^{n-1} \\lambda_i^{(k)} v_i \\otimes \\overline{v_i}$, \n", + "\n", + "$\\hspace{0.75cm}$where $\\lambda_i^{(k)}$ is the $i$th eigenvalue of $T_{:,:,k}$.\n", + "\n", + "$\\hspace{0.75cm}$Thus, we get \n", + "\n", + "$\\hspace{1.25cm}$$\\displaystyle T = \\sum_{i=0}^{n-1} v_i \\otimes \\overline{v_i} \\otimes \\begin{pmatrix} \\lambda_i^{(1)} \\\\ \\vdots \\\\ \\lambda_i^{(n-1)}\\end{pmatrix}$.\n" + ] + }, + { + "cell_type": "markdown", + "id": "57638421-d883-4f7f-ab0f-aae0505cffca", + "metadata": {}, + "source": [ + "***" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "c7140595-9875-4201-ad30-e186a6722ff8", + "metadata": {}, + "source": [ + "## Exercise 2.2 \n", + "\n", + "We now want to consider a very well-known model in statistical physics: the Ising model. The Ising model is a simple mathematical model used to study phase transitions in magnetic systems. It consists of a lattice arrangement of spins, each of which can point \"up\" or \"down\", and interactions between neighboring spins.\n", + "\n", + "The Hamiltonian function (i.e., the energy of a configuration) for the Ising model can be expressed as follows:\n", + "\n", + "$\\hspace{1cm}$$ \\displaystyle H(\\sigma )=- J \\sum _{\\langle ij\\rangle }\\sigma _{i}\\sigma _{j}-h \\sum _{j} \\sigma _{j}$.\n", + "\n", + "Here, $\\sigma_i$ and $\\sigma_j$ represent the spins at lattice points $i$ and $j$, respectively, where $1 \\leq i,j \\leq d$. $J$ is the interaction constant between spins, $h$ represents the external magnetic field, and the summation is over neighboring lattice points $i$ and $j$.\n", + "\n", + "Since this is a system with nearest-neighbor interactions, we can represent it using a (non-cyclic) SLIM decomposition. \n", + "\n", + "**a)**$\\quad$Find the vectors $S$, $L$, $I$, and $M$ such that the tensor train defined by\n", + "\n", + "$\\hspace{1.25cm}$$\\displaystyle \\mathbf{H} = \\begin{bmatrix} S & L & I\\end{bmatrix} \\otimes \\begin{bmatrix} I & 0& 0\\\\ M & 0& 0\\\\ S & L & I \\end{bmatrix} \\otimes \\dots \\otimes \\begin{bmatrix} I & 0& 0\\\\ M & 0& 0\\\\ S & L & I \\end{bmatrix} \\otimes \\begin{bmatrix} I\\\\ M\\\\ S\\end{bmatrix}$\n", + "\n", + "$\\hspace{0.75cm}$satisfies $\\mathbf{H}_{\\nu_1, \\dots, \\nu_d} = H(\\sigma_1, \\dots, \\sigma_d)$, where\n", + "\n", + "$\\hspace{1.25cm}$$\\nu_i = \\begin{cases} 1 & \\text{if}~\\sigma_i = +1, \\\\ 2 & \\text{if}~\\sigma_i=-1.\\end{cases}$\n", + "\n", + "$\\hspace{0.75cm}$Compare your result with Ising model representation provided by Scikit-TT (use $J=h=1$ and $d=5$):\n", + "\n", + "$\\hspace{1.25cm}$1) Import ```scikit_tt.tensor_train``` as ```tt```.\n", + "\n", + "$\\hspace{1.25cm}$2) Use tt.build_core for constructing the cores of $\\mathbf{H}$, e.g., ```cores[i] = tt.build_core([[I, 0, 0], [M, 0, 0], [S, L, I]])``` for $i=1,2,3$.\n", + "\n", + "$\\hspace{1.25cm}$3) Create the class instance by ```H = TT(cores)```.\n", + "\n", + "$\\hspace{1.25cm}$4) Import ```scikit_tt.models``` as ```mdl``` and load the model from Scikit-TT by ```H_SKTT = mdl.ising(d,J,h)```.\n", + "\n", + "$\\hspace{1.25cm}$5) Compute the error in the Frobenius norm: ```(H-H_SKTT).norm()```.\n", + "\n", + "**b)**$\\quad$*extra task:* The configuration probability, i.e., the probability that the system is in a state $\\sigma$ in equilibrium, is given by\n", + "\n", + "$\\hspace{1.25cm}$$\\pi(\\sigma) = \\frac{e^{- \\beta H(\\sigma)}}{Z}$,\n", + "\n", + "$\\hspace{0.75cm}$where $Z = \\sum_{\\sigma} e^{-\\beta H(\\sigma)}$ is a normalization constant and $\\beta = 1/(k_B \\cdot T)$ is the inverse temperature with Boltzmann constant $k_B$. \n", + "\n", + "$\\hspace{0.75cm}$Suppose $\\beta=1$, can you derive a TT decomposition for $\\pi$?" + ] + }, + { + "cell_type": "markdown", + "id": "83e132fb-95b5-43a8-98ca-0ed7c397767d", + "metadata": {}, + "source": [ + "***" + ] + }, + { + "cell_type": "markdown", + "id": "1cb59312-1033-4dd0-886c-7a034eae661b", + "metadata": {}, + "source": [ + "$\\textcolor{red}{\\textbf{SOLUTION:}}$" + ] + }, + { + "cell_type": "markdown", + "id": "81208cd9-c190-46c9-95f8-d0d23778ebce", + "metadata": {}, + "source": [ + "**a)**$\\quad$The component $S$ represents the single-cell interactions, i.e., $S = -h\\begin{pmatrix}+1 \\\\-1\\end{pmatrix}$.\n", + "\n", + "$\\hspace{0.75cm}$$L \\otimes M$ represents the two-cell (nearest-neighbor) interactions, i.e., $L = -J\\begin{pmatrix}+1 \\\\-1\\end{pmatrix}$ and $M = \\begin{pmatrix}+1 \\\\-1\\end{pmatrix}$.\n", + "\n", + "$\\hspace{0.75cm}$$I$ is simply the vector $\\begin{pmatrix}1 \\\\1\\end{pmatrix}$.\n" + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "id": "21650b1d-183d-403f-80d5-9f27c21e3c69", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n", + "Tensor train with order = 10, \n", + " row_dims = [2, 2, 2, 2, 2, 2, 2, 2, 2, 2], \n", + " col_dims = [1, 1, 1, 1, 1, 1, 1, 1, 1, 1], \n", + " ranks = [1, 3, 3, 3, 3, 3, 3, 3, 3, 3, 1] \n", + "\n", + "error: 1.9946147148160272e-13\n" + ] + } + ], + "source": [ + "import numpy as np\n", + "import scikit_tt.tensor_train as tt\n", + "from scikit_tt.tensor_train import TT\n", + "import scikit_tt.models as mdl\n", + "\n", + "J = 1\n", + "h = 1\n", + "d = 10\n", + "\n", + "S = -h *np.array([1,-1])\n", + "L = -J *np.array([1,-1])\n", + "I = np.array([1,1])\n", + "M = np.array([1,-1])\n", + "\n", + "cores = [None]*d\n", + "cores[0] = tt.build_core([[S, L, I]])\n", + "for i in range(1,d-1):\n", + " cores[i] = tt.build_core([[I, 0, 0], [M, 0, 0], [S, L, I]])\n", + "cores[-1] = tt.build_core([I, M, S])\n", + "\n", + "H = TT(cores)\n", + "print(H, '\\n')\n", + "\n", + "H_SKTT = mdl.ising(d,J,h)\n", + "\n", + "print('error:', (H-H_SKTT).norm())" + ] + }, + { + "cell_type": "markdown", + "id": "4b3280ac-c09b-4d73-a4f7-0e46ab49cbce", + "metadata": {}, + "source": [ + "**b)**$\\quad$The (unnormalized) stationary distribution can be expressed as \n", + "\n", + "$\\hspace{1.25cm}$$\\mathbf{\\Pi} = \\begin{bmatrix} \\begin{pmatrix} e^\\beta \\\\ 0 \\end{pmatrix} & \\begin{pmatrix} 0 \\\\ e^{-\\beta} \\end{pmatrix} \\end{bmatrix} \\begin{bmatrix} \\begin{pmatrix} e^{2\\beta} \\\\ 0 \\end{pmatrix} & \\begin{pmatrix} 0 \\\\ e^{-2\\beta} \\end{pmatrix} \\\\[0.4cm] \\begin{pmatrix} 1 \\\\ 0 \\end{pmatrix} & \\begin{pmatrix} 0 \\\\ 1 \\end{pmatrix} \\end{bmatrix} \\otimes \\ldots \\otimes \\begin{bmatrix} \\begin{pmatrix} e^{2\\beta} \\\\ 0 \\end{pmatrix} & \\begin{pmatrix} 0 \\\\ e^{-2\\beta} \\end{pmatrix} \\\\[0.4cm] \\begin{pmatrix} 1 \\\\ 0 \\end{pmatrix} & \\begin{pmatrix} 0 \\\\ 1 \\end{pmatrix} \\end{bmatrix} \\otimes \\begin{bmatrix} \\begin{pmatrix} e^{2\\beta} \\\\ e^{-2\\beta} \\end{pmatrix} \\\\[0.4cm] \\begin{pmatrix} 1 \\\\ 1 \\end{pmatrix}\\end{bmatrix}$." + ] + }, + { + "cell_type": "markdown", + "id": "179c2c04-c06d-4478-9eca-b662e7fa21ac", + "metadata": {}, + "source": [ + "***" + ] + }, + { + "attachments": {}, + "cell_type": "markdown", + "id": "2e657e1a-f617-4540-89c8-3a51c6cc7ad6", + "metadata": {}, + "source": [ + "## Exercise 2.3\n", + "\n", + "You've already seen that any arbitrary tensor can be converted into the TT format. Let's try this out practically now. Create a random tensor with order $d$ and mode dimensions $n$:\n", + "\n", + "> T = np.random.random_sample(d*[n])\n", + "\n", + "First, we reshape $T$ into a matrix, thereby defining the \"residual tensor\" $T_{\\text{res}}$:\n", + "\n", + "> T_res = T.reshape([n, n**(d-1)])\n", + "\n", + "Now we need to write a for loop in which, at each step, one mode is separated, and the remaining matrix $T_\\text{res}$ is reshaped.\n", + "\n", + "Assume $T_\\text{res}$ is a $r*n \\times n^e$ matrix with $r \\in \\mathbb{N}$ and $e \\leq d-1$.\n", + "\n", + "To isolate the mode, we apply a (reduced) QR decomposition. The matrix $Q$ then forms the TT core for the separated mode, and the matrix $R$ represents the new residual tensor:\n", + "\n", + "> Q, R = np.linalg.qr(T_res)\n", + "> \n", + "> s = Q.shape[1]\n", + "> \n", + "> core = Q.reshape([r, n, 1, s])\n", + "> \n", + "> T_res = R.reshape([s * n, n**(e-1)])\n", + "\n", + "Try it for various values ​​of d and n, but initially don't choose them too large. What do you notice when you look at the TT ranks?" + ] + }, + { + "cell_type": "markdown", + "id": "a80b8228-1c42-4c43-8c32-6f5a0d6173c2", + "metadata": {}, + "source": [ + "***" + ] + }, + { + "cell_type": "markdown", + "id": "fb078531-7c5e-45f6-b9f2-631d27ad8b7d", + "metadata": {}, + "source": [ + "$\\textcolor{red}{\\textbf{SOLUTION:}}$" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "id": "b5325236-1967-4cae-8d61-4219d477fd2e", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n", + "Tensor train with order = 10, \n", + " row_dims = [3, 3, 3, 3, 3, 3, 3, 3, 3, 3], \n", + " col_dims = [1, 1, 1, 1, 1, 1, 1, 1, 1, 1], \n", + " ranks = [1, 3, 9, 27, 81, 243, 81, 27, 9, 3, 1]\n" + ] + } + ], + "source": [ + "import numpy as np\n", + "\n", + "d=10\n", + "n=3\n", + "cores = [None]*d\n", + "\n", + "T = np.random.random_sample(d*[n])\n", + "T_res = T.copy()\n", + "r = 1\n", + "for i in range(d-1):\n", + " e = d-i-1\n", + " T_res = T_res.reshape([r*n, n**e])\n", + " Q, R = np.linalg.qr(T_res)\n", + " s = Q.shape[1]\n", + " cores[i] = Q.reshape([r, n, 1, s])\n", + " T_res = R.reshape([s * n, n**(e-1)])\n", + " r = s\n", + "cores[-1] = R.reshape([r, n, 1, 1])\n", + "\n", + "T = TT(cores)\n", + "\n", + "print(T)" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "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.10" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}