solution_1_basics.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "2fe29049",
   "metadata": {},
   "source": [
    "# Solutions to Tensor basics"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3c8f3846-988b-4d6b-8f4e-078c5269df10",
   "metadata": {},
   "source": [
    "Here, we cover some very basic concepts of tensor networks:\n",
    "\n",
    "- how to represent tensors numerically, and\n",
    "- what the basic operations are that we need for tensors\n",
    "\n",
    "For these (toy code) exercises will use the python [numpy](https://numpy.org/doc/stable/user/quickstart.html) librabry for simplicity. For production MPS code, you should eventually switch to a more generic library, e.g., [TeNPy](https://github.com/tenpy/tenpy), which provide a more powerfull class for representing tensors - this allows to make use of block-sparseness implied by symmetries and helps with index bookkeeping."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "82932acc-1850-4cda-a34b-2309bde3371d",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "35a02910-a23b-496f-b1af-2dd566a3a20a",
   "metadata": {},
   "source": [
    "For later use, let us also import matplotlib and set some cosmetic display options."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "b7d6aedd",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "np.set_printoptions(precision=5, suppress=True, linewidth=100, threshold=50)\n",
    "plt.rcParams['figure.dpi'] = 150"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0ec28c62-4028-42ee-80ba-6c4370111024",
   "metadata": {},
   "source": [
    "## Initializing tensors in numpy\n",
    "\n",
    "From a mathematical point of view, tensors are multi-linear function between vector spaces.\n",
    "From a numerical point of view, however, tensors are just collections of numbers arranged in a multi-dimensional grid, with element-wise addition and tensordot contractions over one or multiple legs.\n",
    "\n",
    "In this set of toycodes, we represent tensors directly as numpy arrays. \n",
    "The following examples demonstrate different ways to initialize tensors, here for the usual `2x2` Pauli matrices:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "778f138c-f9b3-4889-baca-e1000d4bf05e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1., 0.],\n",
       "       [0., 1.]])"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# np.eye = identity matrix of given dimension 2\n",
    "Id = np.eye(2)  \n",
    "Id"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "88b69153-267d-47eb-bca8-f7f0e97ad0bc",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0., 1.],\n",
       "       [1., 0.]])"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# specify each element \n",
    "X = np.array([[0, 1.], [1., 0.]])  \n",
    "X"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "3714f0a7-2b92-46df-bc61-f7ed5a91cd4c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.+0.j, -0.-1.j],\n",
       "       [ 0.+1.j,  0.+0.j]])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# initialize empty matrix of dimensions 2x2\n",
    "Y = np.zeros([2, 2], dtype=np.complex_) \n",
    "# and explicitly set non-zero entries\n",
    "Y[0, 1] = -1.j  # j = imaginary unit\n",
    "Y[1, 0] = 1.j\n",
    "Y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "3dd23ebe-5ed1-4b3e-b0e7-cfa81bd7caab",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.,  0.],\n",
       "       [ 0., -1.]])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# only specify the diagonal of a matrix\n",
    "Z = np.diag([1., -1.]) \n",
    "Z"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2cc74936-5116-4ec4-824d-f6d2f4e877cc",
   "metadata": {},
   "source": [
    "They are all simple `2x2` matrices, as we can see by checking the `shape`. Only the `Y` matrix has complex data type (`dtype`):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "630abe6c-fc4c-48c4-8382-2ddb988d4da8",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(2, 2) float64\n",
      "(2, 2) float64\n",
      "(2, 2) complex128\n",
      "(2, 2) float64\n"
     ]
    }
   ],
   "source": [
    "for M in [Id, X, Y, Z]:\n",
    "    print(M.shape, M.dtype)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1a1d2ea4-17d8-4de8-8133-175716383bee",
   "metadata": {},
   "source": [
    "If we consider the Pauli matrices as spin-1/2 operators, we already implicitly assumed the basis to be $\\lvert\\uparrow\\rangle, \\lvert\\downarrow\\rangle$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "3b336217-4c28-41ce-b351-e15906f68f0c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1., 0.])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "up = np.array([1., 0])\n",
    "up"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "a2b446cc-8f55-42e2-8bc0-c4712eb6a4f7",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0., 1.])"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "down = np.array([0., 1.])\n",
    "down"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ce9d73e7-4c8e-4894-9974-dfb9add88226",
   "metadata": {},
   "source": [
    "## Element-wise operations\n",
    "Element-wise addition of tensors and multiplication/division with scalars works as you might expect."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "e537ae4e-4fd0-4117-b4ab-96568668bc76",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.,  5.],\n",
       "       [ 5., -1.]])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Z + 5 * X"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "b0289b7d-2c78-4882-ba49-98f7f6cd8dbe",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.70711, 0.70711])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plus_x = (up + down) / np.sqrt(2)  # sqrt = square root\n",
    "plus_x"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b2697293-28ba-43e7-abdb-6e2be7cb24d2",
   "metadata": {},
   "source": [
    "Complex conjugation can be acchieved with `np.conj()` (or the `conj()` method of numpy arrays)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "1b6d2380-fea2-4425-b9f1-4eb2d7a3a5dd",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([[ 0.+0.j, -0.-1.j],\n",
       "        [ 0.+1.j,  0.+0.j]]),\n",
       " array([[ 0.-0.j, -0.+1.j],\n",
       "        [ 0.-1.j,  0.-0.j]]))"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Y, np.conj(Y)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6f9bbba5-dd3c-4f71-bf3a-1665e2cccc3f",
   "metadata": {},
   "source": [
    "### Exercise: init `minus_x`\n",
    "\n",
    "Initialize the state `minus_x`  $= \\frac{1}{\\sqrt{2}} \\big( \\lvert\\uparrow \\rangle - \\lvert\\downarrow\\rangle \\big)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "7940af87-faa4-48fb-bfd7-e676b151bd64",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.70711, -0.70711])"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "minus_x = (up - down) / np.sqrt(2)\n",
    "minus_x"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "887424c9-b3d5-4ecd-9245-cfcdc396fbcc",
   "metadata": {},
   "source": [
    "## Random vector\n",
    "We can also define a small utility function to define a random (complex) vector pointing uniformly in any direction:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "49e0f906-3578-48ae-a5b7-a531ececc671",
   "metadata": {},
   "outputs": [],
   "source": [
    "def random_vector(size, rng=np.random.default_rng(seed=1234)):\n",
    "    \"\"\"Return a normalized random vector of given shape pointing uniformly in any direction.\"\"\"\n",
    "    vector = rng.normal(size=size) + 1.j * rng.normal(size=size)\n",
    "    return vector / np.linalg.norm(vector)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "01f4dabc-1057-499d-9865-da04ec5824fe",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-0.45581+0.04337j,  0.01822+0.24548j,  0.21056+0.82791j])"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "random_vector(3)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4edce430-be9d-4c6b-8d5f-79946487b0b5",
   "metadata": {},
   "source": [
    "## Leg label convention\n",
    "\n",
    "Note that the order of tensor legs is important; and we need to keep track of it. \n",
    "Libraries like TeNPy can help with that bookkeeping, but here we keep it simple, and just write leg labels as comments in the code.\n",
    "\n",
    "- For ket (vector) indices, we just use short identifiers, e.g. `i`, `j`, `i1`, `i2`, `vL`.\n",
    "- For bra (dual) indices, we use an additional `*`, so calling `conj()` on a tensor toggles stars on all legs, and (square) operators should have       combinations of non-star and star labels.\n",
    "- When contracting (examples further down), we indicate labels that are getting contacted over with square brackets.\n",
    "- When grouping legs (examples further down), we combine labels in round brackets separated by dots, e.g. `(i.j)`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "3d961ef6-06e0-472b-9f64-040727e43039",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.70711, 0.70711])"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "plus_x # leg i\n",
    "X # legs i i*\n",
    "np.conj(plus_x) # leg i*"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "12c3ddef-da67-4a92-ade2-401cdbe32781",
   "metadata": {},
   "source": [
    "## More spins and legs, transpose and reshape to combine and split legs\n",
    "\n",
    "When we have more than one spin, the Hilbert space is a an outer product of the local Hilbert spaces. In tensor network language, this corresponds to the fact that we can have one leg per spin.\n",
    "\n",
    "For example, the singlet state $\\frac{1}{\\sqrt{2}}\\big(\\lvert \\uparrow \\rangle \\otimes \\lvert\\downarrow \\rangle - \\lvert \\downarrow \\rangle \\otimes \\lvert\\uparrow \\rangle \\big)$ on two spins could be represented by a two-leg tensor:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "71d6e1b3-4d5e-4da8-a22e-dc3d57cadf43",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.     ,  0.70711],\n",
       "       [-0.70711,  0.     ]])"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "singlet = np.zeros([2, 2])  # legs i, j\n",
    "singlet[0, 1] =  1./np.sqrt(2)\n",
    "singlet[1, 0] = -1./np.sqrt(2)\n",
    "singlet"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "38943568-c256-4fce-b05f-3487d684ab38",
   "metadata": {},
   "source": [
    "Alternatively, we can enumerate all possible combinations basis states (in lexiographical order),\n",
    "in the case of two spins $\\lvert\\uparrow\\uparrow\\rangle, \\lvert\\uparrow\\downarrow\\rangle, \\lvert\\downarrow\\uparrow\\rangle, \\lvert\\downarrow\\downarrow\\rangle$, and use a single leg.\n",
    "\n",
    "In other words, we can **combine** (or \"**group**\") the two legs to a single one.\n",
    "For numpy arrays, this corresponds to a `reshape`, taking dimensions `d1` and `d2` of the individual legs to `d1 * d2`. It is merely a re-arranging of the tensors in a different grid, the entries are still the same:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "66602ef4-bb7b-4110-a109-f9c4b611284f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([ 0.     ,  0.70711, -0.70711,  0.     ])"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "singlet_group = np.reshape(singlet, [2*2]) # legs i j -> (i.j)\n",
    "singlet_group"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1cd3ee05-027d-47bb-9339-687a4c1f7367",
   "metadata": {},
   "source": [
    "If we remember the order and dimensions of the combined legs, we can split it again by another reshape:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "2f96f59b-d094-4fad-b761-e706b6d540d3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.     ,  0.70711],\n",
       "       [-0.70711,  0.     ]])"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "singlet_split = np.reshape(singlet_group, [2, 2]) #  legs (i.j) -> i j\n",
    "singlet_split"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "955235e2-42ba-4d2a-b52d-7aa5aa2cd7a3",
   "metadata": {},
   "source": [
    "All of this generalizes to multiple legs, as long as we `transpose` legs to be combined next to each other.\n",
    "\n",
    "For example, if we initialize a random vector $\\theta_{a,b,c,d}$ and want to group legs as `a b c d -> (a.c) (b.d)`, we need to transpose first:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "5524ef3d-bd93-4b8e-a0c8-135108329411",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(2, 3, 4, 5)\n"
     ]
    }
   ],
   "source": [
    "theta = random_vector([2, 3, 4, 5])  # legs a b c d\n",
    "print(theta.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "65f18475-cd12-46a2-b295-0fe26e3f30b9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(8, 15)"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "theta_group = np.reshape(np.transpose(theta, [0, 2, 1, 3]), # a b c d -> a c b d\n",
    "                         [2*4, 3*5]) # a c b d -> (a.c) (b.d)\n",
    "theta_group.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "85009e8f-c5d6-4460-b386-9ccc878c2054",
   "metadata": {},
   "outputs": [],
   "source": [
    "theta_split = np.reshape(theta_group, [2, 4, 3, 5]) # (a.c) (b.d) -> a c b d"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c2748951-1359-4f72-a67d-b7cb52f478a4",
   "metadata": {
    "tags": []
   },
   "source": [
    "To get back to the original form, we need to transpose again:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "f4afaf91-aa36-4d31-b6a1-2c07d85a3504",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.0"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.linalg.norm(theta_split.transpose([0, 2, 1, 3])  # a c b d -> a b c d\n",
    "               -theta)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fabf5b0a-e6ff-4ed6-9fcc-7b6772230f77",
   "metadata": {},
   "source": [
    "### Exercise: reshape\n",
    "\n",
    "Initialize a random vector `psi` with 6 legs of length 2, i.e. shape `(2, 2, 2, 2, 2, 2)`. \n",
    "\n",
    "Afterwards, reshape it into a matrix with each the first and last 3 legs combined.\n",
    "\n",
    "*Tip*: In Python, you can use `[2]*5` to get a list `[2, ..., 2]` with 5 entries, \n",
    "and `2**5` $= 2^5$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "127e0ab6-71b2-4377-9afa-7c8b61ff75ad",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2, 2, 2, 2, 2, 2)"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "psi = random_vector([2]*6)  # i j k l m n \n",
    "psi.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "c8df83fa-3310-40fc-9144-db3bc9ea3fba",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(8, 8)"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "psi = np.reshape(psi, [2**3, 2**3]) # i j k l m n -> (i.j.k) (l.m.n)\n",
    "psi.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "79bd7296-654f-4188-9f0a-0c6691f353d7",
   "metadata": {},
   "source": [
    "## Contractions with tensordot\n",
    "\n",
    "The `np.tensordot` function performs contractions between two tensors along legs specified in the `axes` argument.\n",
    " \n",
    "For example for the usual matrix-vector product $\\hat{X} \\lvert \\uparrow\\rangle$, we need to contract the second leg (index 1, since we start counting at 0) of `X` with the first (and only, index 0) leg of the `up` vector:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "b1085e85-9b51-4984-b8e6-3e518c7790f6",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0., 1.])"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.tensordot(X, up, axes=(1, 0)) # i [i*], [i] -> i"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "17d06f51-a069-42e4-aff4-3afe21199865",
   "metadata": {},
   "source": [
    "Matrix multiplication looks similar, except that the second and the resulting tensor have one more leg:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "1a467af4-0950-4628-8b18-a5aeda1f285a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.+1.j, 0.+0.j],\n",
       "       [0.+0.j, 0.-1.j]])"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.tensordot(X, Y, axes=(1, 0))  # i [i*], [i] i* -> i i*"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "df541397-b01a-4eb8-be18-fd45a1d8d144",
   "metadata": {},
   "source": [
    "If we contract over all legs, we effectively perform an inner product and obtain a scalar:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "id": "5349fd63-eff5-42ee-a760-a871f62d77a8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array(0.70711)"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.tensordot(np.conj(plus_x), up, axes=(0, 0))  # [i*], [i] -> scalar"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7d663f67-bef0-448e-a1ac-2cfefcba2167",
   "metadata": {},
   "source": [
    "If we set `axes=0`, we obtain an outer product, e.g. if we want to get the operator `ZZ` $= Z_1 \\otimes Z_2$ acting on two spins:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "id": "be032d25-c88b-4268-8673-592c4b304624",
   "metadata": {},
   "outputs": [],
   "source": [
    "ZZ = np.tensordot(Z, Z, axes=0)  # i i*, j j* -> i i* j j*"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0717d10b-fb92-43ff-9155-9151354b455d",
   "metadata": {},
   "source": [
    "For printing the result, let us group the in- and outgoing legs:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "id": "64eacdba-1ab1-4bf7-8e7e-6c695bc9d2ce",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.,  0.,  0.,  0.],\n",
       "       [ 0., -1.,  0.,  0.],\n",
       "       [ 0.,  0., -1.,  0.],\n",
       "       [ 0.,  0.,  0.,  1.]])"
      ]
     },
     "execution_count": 30,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ZZ_g = np.transpose(ZZ, [0, 2, 1, 3])  # i i* j j* -> i j i* j*\n",
    "ZZ_g = np.reshape(ZZ_g, [2*2, 2*2])    # i j i* j* -> (i.j) (i*.j*)\n",
    "ZZ_g"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f0c867a0-957c-499c-9767-b0cbd21badbe",
   "metadata": {},
   "source": [
    "Finally, to contract multiple legs, we can pass two lists to `axes`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "id": "c5c4ded0-74b1-4bcd-a678-bcd593072d8e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2, 2, 3, 4)"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ZZ # i i* j j*\n",
    "Theta = random_vector([3, 2, 2, 4]) # a i j b\n",
    "# contract i* of ZZ with \n",
    "ZZ_Theta = np.tensordot(ZZ, Theta, axes=([1, 3], [1, 2])) # i [i*] j [j*], a [i] [j] b -> i j a b\n",
    "ZZ_Theta.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "153f8352-8ad1-4df6-99a9-a3001d243483",
   "metadata": {},
   "source": [
    "### Exercise: contractions\n",
    "\n",
    "- Contract the following expectation values using `tensordot` and check that you get the correct results:\n",
    "  - `<plus_x|plus_x> =?= 1.`\n",
    "  - `<plus_x|X|plus_x> =?= 1.`\n",
    "  - `<plus_x|Z|plus_x> =?= 0.`\n",
    "  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "id": "afa30bf1-dce9-4abe-9bb7-117a79988bf2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array(1.)"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.tensordot(np.conj(plus_x), plus_x, axes=[0, 0]) # [i*], [i] -> scalar"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "id": "2e933507-f357-4f44-aaf9-63633278bbf7",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array(1.)"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.tensordot(np.conj(plus_x), # i*\n",
    "             np.tensordot(X, plus_x, axes=[1, 0]),  # i [i*], [i] -> i\n",
    "             axes=[0, 0])  # [i*], [i] -> scalar"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "id": "35340495-7b0d-4981-a3a6-e9f3f6c207c5",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array(-0.)"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.tensordot(np.conj(plus_x), # i*\n",
    "             np.tensordot(Z, plus_x, axes=[1, 0]),  # i [i*], [i] -> i\n",
    "             axes=[0, 0])  # [i*], [i] -> scalar"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ab274679-bdb4-4643-8eb7-fbe43e0087eb",
   "metadata": {},
   "source": [
    "\n",
    "- Define the operator `h` $ = 4 \\vec{S_1} \\cdot \\vec{S_2} = X_1 X_2 + Y_1 Y_2 + Z_1 Z_2 $ with 4 legs.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "id": "3d450e8a-d2c0-42ed-90d6-e3ff8c8ab4ca",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.,  0.,  0.,  0.],\n",
       "       [ 0., -1.,  2.,  0.],\n",
       "       [ 0.,  2., -1.,  0.],\n",
       "       [ 0.,  0.,  0.,  1.]])"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "h = sum([np.tensordot(A, A, axes=0) for A in (X, Y, Z)])  # i i*, j j* \n",
    "h_g = np.transpose(h, [0, 2, 1, 3])  # i i* j j* -> i j i* j*\n",
    "h_g = np.reshape(h_g, [2*2, 2*2])    # i j i* j* -> (i.j) (i*.j*)\n",
    "np.real_if_close(h_g)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "45090e1d-84df-4f7e-b264-a4ad9894b049",
   "metadata": {},
   "source": [
    "- Using `tensordot`, contract and check `<singlet|h|singlet> =?= -3`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "id": "a53b92c3-6754-4805-9b05-d25af90ff4d8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array(-3.+0.j)"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "h_singlet = np.tensordot(h, singlet, axes=([1, 3], [0, 1])) # i [i*] j [j*], [i] [j] -> i j\n",
    "singlet_h_singlet = np.tensordot(np.conj(singlet), h_singlet, axes=([0, 1],[0, 1])) # [i*] [j*], [i] [j] -> scalar\n",
    "singlet_h_singlet"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4ba0a7f8-20d7-4432-aaf8-21a394ae87f8",
   "metadata": {},
   "source": [
    "## Schmidt decomposition with SVD"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "16143401-d426-4c38-88b7-96219a4b8059",
   "metadata": {},
   "source": [
    "The final tool we need are matrix operations, e.g., diagonalization, SVD or QR decomposition.\n",
    "The idea is to just group legs to view the tensor as a matrix and then perform the corresponding operation.\n",
    "\n",
    "The SVD is of particular importance because it eventually allows truncation. We can use it to calculate the entropy of the singlet:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "id": "6579de4d-83a0-4879-9223-cce6555566bc",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.     ,  0.70711],\n",
       "       [-0.70711,  0.     ]])"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# define singlet again, same as above\n",
    "singlet = np.zeros([2, 2])  # i j\n",
    "singlet[0, 1] =  1./np.sqrt(2)\n",
    "singlet[1, 0] = -1./np.sqrt(2)\n",
    "singlet"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "id": "8366fcdb-accc-4005-b4ab-2953840539ae",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.70711, 0.70711])"
      ]
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A, S, B = np.linalg.svd(singlet, full_matrices=False) # i j -> i a, a, a j\n",
    "S  # Schmidt values"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e43c8c50-8b8a-45b0-b03a-fcf5120ea0f3",
   "metadata": {},
   "source": [
    "We can also get the original tensor back by contracting `A` with `diag(S)` and `B`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "id": "ddbec137-03f3-498a-87b8-1b0924e18839",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 0.     ,  0.70711],\n",
       "       [-0.70711,  0.     ]])"
      ]
     },
     "execution_count": 39,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "contr = np.tensordot(np.diag(S), B, axes=[1, 0]) # a [a*], [a] j -> a j\n",
    "contr = np.tensordot(A, contr, axes=[1, 0]) # i [a], [a] j -> i j\n",
    "contr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "id": "1a84667a-af65-413a-b896-a232b1240dd5",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.9999999999999998"
      ]
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "np.sum(S*S)  # normalization of the state implies sum(S*S) = 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "id": "59c65227-6251-4379-9bc5-9e0b3ec660ac",
   "metadata": {},
   "outputs": [],
   "source": [
    "def entropy(sing_vals):\n",
    "    \"\"\"Return von Neumann entropy.\"\"\"\n",
    "    # first discard small singular values, since 0 * log(0) -> 0\n",
    "    # but log(0) not defined\n",
    "    S = sing_vals[sing_vals>1.e-15]\n",
    "    S2 = S*S\n",
    "    S2 /= np.sum(S2)  # ensure normalization\n",
    "    return - np.sum(S2 * np.log(S2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "id": "df8623fe-ba50-43cc-9da2-d2a95337f238",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.0"
      ]
     },
     "execution_count": 42,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "entropy(S) / np.log(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9d12f2ac-f8d3-4d93-ab5f-261836b60fea",
   "metadata": {},
   "source": [
    "## Truncation of Schmidt values\n",
    "\n",
    "Let us now look at at spin system with `L=12` sites,\n",
    "and consider a bipartition into the left 6 and the right 6 spins.\n",
    "\n",
    "We already know how to get a random state with the left and right legs combined:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "id": "b489da72-df30-45a0-baa8-b5340bba372a",
   "metadata": {},
   "outputs": [],
   "source": [
    "L = 12"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "id": "c65d6cc9-2ab4-468c-8b29-62fa7d105fe3",
   "metadata": {},
   "outputs": [],
   "source": [
    "psi_random = random_vector([2]*L)  # L legs\n",
    "psi_random = np.reshape(psi_random, [2**(L//2), 2**(L//2)])  # combined legs (i0.i1...i{L/2-1}) (i{L/2}...i{L})"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a74c2181-857f-40cb-ad45-3b873b3497b8",
   "metadata": {},
   "source": [
    "It turns out that random vectors are highly entangled with the entanglement entropy reaching nearly the maximal value $\\frac{L}{2} \\log(2)$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "id": "b6b21c51-2900-4e2e-a052-7d84241e0ab5",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "5.282231767294896"
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A_random, S_random, B_random = np.linalg.svd(psi_random, full_matrices=False)\n",
    "entropy(S_random) / np.log(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "42ef5003-bb39-4ce7-be15-b610cf6ea2fa",
   "metadata": {},
   "source": [
    "Ground states of local Hamiltonians behave very differently; in 1D the are law of the entanglement predicts a scaling constant in $L$.\n",
    "(An exception are critical, gapless points where CFT predicts $S \\propto \\log(L)$.)\n",
    "\n",
    "The following snippet uses exact diagonalization (and the Lanczos algorithm) to obtain the ground state of the transverse field ising model with $J=1, g=1.1$. (You can consider it a black bock for now - here we are interested in the properties of the state only.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "id": "16fafae8-b956-4bc2-a8c0-17d9e4c9053e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(4096,)"
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from tenpy_toycodes.tfi_exact import finite_gs_energy\n",
    "\n",
    "energy, psi_ground_state = finite_gs_energy(L=L, J=1., g=1.1, return_psi=True)\n",
    "psi_ground_state.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "id": "c08b51e3-3607-4069-966d-c69b694be19b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.44561532410538357"
      ]
     },
     "execution_count": 47,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "psi_ground_state = np.reshape(psi_ground_state, [2**(L//2), 2**(L//2)]) # combined left and right legs\n",
    "A_ground_state, S_ground_state, B_ground_state = np.linalg.svd(psi_ground_state, full_matrices=False)\n",
    "entropy(S_ground_state) / np.log(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8bcedc95-0163-4d63-b553-9954e61f09e2",
   "metadata": {},
   "source": [
    "We can see this behaviour reflected in the decay of singular values:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "id": "e76952d4-401d-459e-9721-15241b22a6c5",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7fb765146380>"
      ]
     },
     "execution_count": 48,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmUAAAG7CAYAAACCU+LlAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAABcSAAAXEgFnn9JSAABu80lEQVR4nO3deVzU1f4/8NegLLKISqDIIgoiSopIoLgkuJQLKqll5lWU7Fo3Lc3UMlGUQq+lZv1a702wW5qmqaF5c0NRUdxQDHcE0XDBBWUbQDi/P/jOXCdmgJn5MDPA6/l4zOPeOcvnnM/HT87b8zmfc2RCCAEiIiIiMiozY3eAiIiIiBiUEREREZkEBmVEREREJoBBGREREZEJYFBGREREZAIYlBERERGZAAZlRERERCaAQRkRERGRCWBQRkRERGQCGJQRERERmQAGZUREREQmgEEZERERkQlgUEZERERkApoauwNUd9q0aYPCwkK4u7sbuytERET1SnZ2NmxsbHDr1i2DtcmgrAErLCxEWXEBkHtBfYGnvAEZB0uJiIj+qqysDIWFhQZtk0FZA+bu7g7kXkD6P2zVF5h/ArCwMWyniIiI6gFfX1+Dt8lhEiIiIiITwKCMiIiIyATw8WVjJwRQVqRbXXNrQCaTtj9ERESNFIOyxq6sCIhtq1vd+Tmck0ZERCQRPr4kIiIiMgEcKWvgbjRxw1SXeAS2d8DoHi5wtLX8X6a5te6PLomIiEhSDMoauOIygT0ZhdiTUYhPEq9jbIAbFo3oAivzJnXSXm5+CTYcz0ZK5n0UlDyGrWVT9OrggJeecYOjnWWd1jdm20RERPqSCSGEsTtBdcPX1xeXbxeg7dQvVdJ7tm+FtZFBlYFZaaEkc8rkZeVYnJCOTSdvoKy86i1l3kRWbUCoT31jtq3AYJSIqGFRrFOWnp5usDYZlDVgmoIyABgf5I6lo7tKEpTJy8oRseYYUjLv11hFJSD8P/rUN2bbivoMRhmMElHDw6CMJFVdUGbeRIbk9wbC0fKxzkFZ7luZcGzVCu//kob1x67Xup4yIPw/+tQ3ZtsMRhmMMhglargYlJGkqgvKAGDO853wZh9nnYOyb/oewgs9O6LPsn1qf9w0UQaEdpa4ky/Xub6AMFrbjnaWDEYZjFapw2C0esYMRom0ZYygjBP9G7GjV+9VBmU6Op55D2VmVloFNQBQVi6w8cR1vBnqhY3Hr+tcXwhhtLZffMYVm07e0KruppPX8c5gb2VAqGt9AWG0th3tLLE4Ib1WQQ0ApGTex+KEcyoBoT71jdl2bQK6snKB9ceycTW3QKeAUFN9Y7atqK8poDt4+S4+3XNJ54Cwpvr6tq1QX0cnGYw2LgzKTEBSUhJWrFiB06dPIzs7G4sWLUJ0dHSdt1tQ8li/+qXltf6B+6ujV+/hzVAvverrSoq2jRkQMhg1fNsMRutnMKo4hrECQgajDEa1xaDMBBQUFKBLly545ZVXMHPmTIO1a2up3x+/rUUT3NMxsFMEhLoGhvoElFK0zWC09hiMMhg1VjBaX0cnGYyaRjD6Z14xKioMO8Or0a3of/LkSSxbtgyjR4+Gi4sLZDIZrKysaqwnl8uxaNEieHt7w8rKCm3btkVkZCRu3NDuLzt1hg0bhqVLl2LcuHGwtDTcvwB6dXDQq35gewedAztFPX3qG7NtfYNJfeobs219g0l96huzbX0COgB61Tdm27oGdLn5JQCgV3192wZ0C+qkqGvsthVB3fpjmv/sFUFdxJpjkJeVS1LX2G0rjvH+L2novWwvPtl1CQcv30Vqdh4OXr6Lj3+/iN7L9uL9X87Wqm5RSTnkZRVq+1FXGl1QFhMTg/fffx9btmxBTk5OrerI5XIMHDgQS5YsQUFBAUaNGgU3NzfExcWhR48eyMjIqONeS8+8iQwvPeOm1zFG93BBz/atdKqrCAj1qW/MthmMaofBKINRbevr27YxA0IGow0zGDWERheUBQcHY+HChUhISMCtW7dqVSc2NhbJyckIDg7GpUuXsGHDBqSkpGDFihXIzc1FZGSkSvn79+/jypUr1X7u3LlTF6dXa2MD9H8u72hriZcC3WDeRKZVvScDQn3qG7NtBqPaYTDKYFTb+vq2XV9HJxmM6tY2YLiAsC41uqBs3rx5WLx4McLCwtC6desay5eVleHzzz8HAHzxxRewtbVV5r3zzjvo1q0bkpKScPLkSWX6Z599ho4dO1b7mTt3rvQnV0s927fCohFdKr+YWwPzcyCfcx2TWm9FZ/maGj+TWm+FfM51wNwaTnZWGBvgqlX7TwaE+tQ3ZtsMRrXDYJTBqLb19W27vo5OMhjVrW1DB4R1pdEFZdo6dOgQ8vLy4OnpCX9//yr5Y8eOBQAkJCQo0xYuXIiysrJqP2vWrDHYOSiYN5FhfJC76qRQmQywsIGVTXN8O/VZhAd543GTZiiGVZXP4ybNEB7kjW+nPgsrm+aVdQEsGuFb6x88lYDw/+hT31htMxitPQajDEZ1qa9v2/V1dJLBqG5tGzogrCsMympw5swZAECPHj3U5ivSFeUAwMzMDE2bNq32Y2ZmmEtvZW6Gfh2fwpznOyH5vYFYOrqrxjdWrMybYOnorkh+byDmPN8J/To+BX/3FjXWtzJvgrWRQRgf5K7xR09tQChBfWO2zWC0dhiMMhjVpb6+bdfX0UkGo7q1bYyAsC5wSYwaZGdnAwBcXdX/Za5IV5TTRUFBAa5cuQIAKC0txa1bt3D69GlYWFigS5cuNdT+36rDf5WRkQFPT0/859WeWvXH0c4Sb4Z64c1Qr1rXUQR07wz2xsYT13H06j2tXkPWp76x2lYEdIsTzmHTSfX/0qpulXR96huz7UUjfHE1t7DWq+KrCwh1rW+sthUBnTa7IKgLCHWtb6y2Xwp0w+q9l7Xe9eLJgFDX+gJCr7Z7tm+Fg5fv1rquQq8ODhBC6FzX2G3rG9TpU9eYbRtjeaW6wKCsBgUFBQAAa2trtfk2NjYq5XRx4sQJhIaGKr9/8803+Oabb9CuXTtkZWXpfFxj0CWgk6q+MdpmMMpgVB0Go8YPRo0ZEDIYNXzbxggI60Kj3/tSJpPB0tIScrlcbf5rr72Gf//731iwYAFiYmKq5F++fBne3t7w9vbGxYsX67q7WjHGvl3UuOTml+gUEEpR3xhtVy5sqVtAqG99Y7Vdn/c6ra/7y+pTt77uJ6xv2/9v32V8sutSresqzHm+E4QQauvm/PsfAIDSu9e0Pq6uGJTVEJS98847WLVqFWbNmoWVK1dWyT9z5gy6d++OHj16qLyBaQoYlBHVDQajqkw1GDVmQMhg1LBt10VAyKDMCGoKyj799FPMmjULL774IjZu3Fglf8eOHQgLC0N4eDi2bNlS193VCoMyIpJSfQxG6+PopBR1GYzWrKaAkEGZEdQUlCUmJmLAgAHw9PRUTsZ/UkxMDBYuXIiFCxdi8eLFdd1drTAoIyKqVN9GJ/Wty2C0erUJCBmUGUFNQVlpaSmcnJzw8OFDnDp1qspaZX5+fkhLS8OxY8cQGBhoiC7XGoMyIqLGjcGoKm0CQgZlRlBTUAYACxYswEcffYTevXtj165dyjcuV65cidmzZ6Nv3744ePCgobpcawzKiIioMZIiIPzgb8+hokLg0c1MA/W6EQZlO3bsUHmLMiUlBTKZDEFBQcq0qKgoDB8+XPldLpcjJCQEKSkpcHZ2Rr9+/XDt2jWkpKTAwcEBR48ehZeXbssw1KU6D8qEAMqKdKtrbq3cEYCIiMjUGGNgw7QW6DCA3NxcpKSkqKQJIVTScnNzVfKtrKyQmJiIpUuXYt26ddi6dStatmyJiIgIxMTEwM3NzSB9NzllRUBsW93qzs8BLGyk7Q8REVE91uhGyhqTOo/ySwsZlBERUYNkjJEy7n1JREREZAIYlBERERGZAAZlRERERCaAQRkRERGRCWBQRkRERGQCGJQRERERmQAGZUREREQmgEEZERERkQlodCv6kwnhNk1ERERKDMrIeLhNExERkRIfXxIRERGZAAZlRERERCaAjy9Jd+bWlY8Rda2r63wyIiKiBohBGelOJuO8LiIiIonw8SURERGRCWBQRkRERGQCGJQRERERmQAGZUREREQmgEEZERERkQlgUEZERERkAhiUEREREZkABmVEREREJoBBGREREZEJ4Ir+ZDz6btNERETUgDAoI+PhNk1ERERKfHxJREREZAIYlBERERGZAAZlJmDTpk0ICgpCq1atYGVlBS8vLyxYsAClpaXG7hoREREZCOeUmYBWrVph7ty56Ny5M2xsbJCamopp06bh0aNH+Oyzz4zdPSIiIjKARheUnTx5Ert378axY8eQkpKCnJwcWFpaQi6XV1tPLpdj6dKlWL9+PbKzs9GqVSsMGTIES5Ysgaurq159GjBggMp3Dw8PJCUlYc+ePXodt0ETAigr0q2uuXXlSwZEREQmpNEFZTExMdi2bZtWdeRyOQYOHIjk5GQ4Oztj1KhRyMrKQlxcHLZv344jR47A09NTsj6eP38eO3fuxHPPPSfZMRucsiIgtq1udefn8K1PIiIyOY0uKAsODoafnx8CAwMRGBiINm3a1FgnNjYWycnJCA4Oxq5du2BrawsAWLlyJWbPno3IyEgcOHBAWf7+/fu4f/9+tcds3rw5nJycVNJsbW1RVlaG0tJSTJs2DatWrdLhDImIiKg+kgkhhLE7YUwymazax5dlZWVwcnJCXl4eTp06BX9/f5V8Pz8/pKWl4cSJEwgICAAAREdHY/HixdW2GxERgfj4eJW0K1euoLi4GCdOnMB7772H119/vcbjVMfX1xcAkJ6ervMxTFZpIUfKiIiozhjjN7TRjZRp69ChQ8jLy4Onp2eVgAwAxo4di7S0NCQkJCiDsoULF2LBggXVHtfMrOqLr15eXgCArl27wszMDJGRkZg7dy5sbBhAEBERNXQMympw5swZAECPHj3U5ivSFeWAyoBLXdClLSEEysrK9D4OERERmT4GZTXIzs4GAI1vWCrSFeV0ERMTg549e6JDhw4QQuDYsWOYN28eRo0ahRYtWtRYXzHE+lcZGRmSvoBAREREdYdBWQ0KCgoAANbW6jfAVjxaVJTTRXFxMWbMmIHr16+jadOm8PDwwKxZs/DWW2/pfEwiIiKqXxiU1UDxHoRMw7pWUrwnERsbi9jYWJ3ra5qEqGkEjYiIiEwPt1mqgZ2dHQCgsLBQbX5RUeUCpoplMoiIiIh0waCsBu7u7gCAGzduqM1XpCvKEREREemCQVkN/Pz8AACnTp1Sm69I79atm8H6RERERA0P55TVoE+fPrC3t0dGRgZSU1OrrFW2adMmAEBYWJgxuke64t6ZRERkYhiU1cDCwgLTp0/HRx99hOnTp2PXrl3KNy5XrlyJtLQ09O3bF4GBgUbuKWmFe2cSEZGJaXRB2Y4dOxATE6OSVlpail69eim/R0VFYfjw4crvCxYswJ49e5CcnIyOHTuiX79+uHbtGlJSUuDg4IC4uDiD9Z+IiIgapkYXlOXm5iIlJUUlTQihkpabm6uSb2VlhcTERCxduhTr1q3D1q1b0bJlS0RERCAmJgZubm4G6Ts9wdy6csRK17q6ProkIiKqI41+Q/KGrEFvSK4vbmhORETVMMZvKN++JCIiIjIBDMqIiIiITIBBg7LLly8bsjkiIiKiesOgQVmvXr3w+uuv4+bNmyrpJSUlOHz4sCG7QkRERGRSDD5S1qxZM/j4+OCDDz7AtWvXcPv2bWzatAkDBgwwZFeIiIiITIpBl8T48ccf8cUXX+Dx48dYunQpli1bBgCQyWSYPn26IbtCREREZFIMOlL24YcfIjY2Fg8fPkRxcTFOnjyJ/v37w9/fH/PnzzdkV4iIiIhMikGDsuLiYowcORJ2dnawtLRE9+7dsW/fPvTr1w/jxo0zZFeIiIiITIpBg7KBAwfi888/r5L+j3/8A0ePHjVkV4iIiIhMikHnlC1duhS9evXCmTNnMH36dPTs2RPm5uZYtWoVHBwcDNkVIiIiIpNi0KDMx8cHhw8fxrvvvotXXnkFT+7wtHr1akN2hRo7fffOJCIikpjBNyT39fXFzp07kZeXh1OnTuHRo0fo1q0bOnToYOiuUGMmk3H/SiIiMikGD8oUWrRowbXJiIiIiP6PUfe+vHfvnjGbJyIiIjIZBhkpKy0tRXp6Os6ePavyuX37Nh4/fmyILhARERGZNMmDsszMTKSlpakEX1euXEF5eTkAqEzub9OmjdTNExEREdVLegdlDx48wEcffYTDhw8jPT0dhYWFKvlCCFhYWMDX1xfdunWDn5+f8n8dHR31bZ6IiIioQdA7KJs8eTK2b9+uHAFzd3fHzZs3UVZWBjc3N/zwww8IDg5G06ZGe6eAiIiIyOTpPdE/MTERLVu2xNatW5GXl4esrCycPHkSQUFBuH79OsaPH4/ffvtNir4SERERNVh6B2UFBQUYNmwYRo4ciebNmwMAnn76aRw5cgQrVqzAw4cP8cILL2D8+PG4e/eu3h0mIiIiaoj0DsqSkpIwderUKukymQyzZs1CWloaBg4ciA0bNsDX1xcbNmzQt0ki4xICKC3U7fPEiy5ERERPkglhmF+J+Ph4zJ49G3l5eRg5ciS+/PJLODs7G6LpRsvX1xcAkJ6ebuSeNDClhUBsW93qzs/hTgJERPWAMX5DDbZ47OTJk3H+/HmMGTMG27Ztg6+vL+Lj4w3VPBEREZFJM+iK/k5OTti4cSO2bNmCZs2a4dVXXzVk80SmgY8/iYhIDaOsUzFq1CiEhIRg3rx5xmieyLjKivj4k4iIqjDa3pf29vb4+uuvjdU8ERERkUkx6obkRERERFSJQZkJiI6Ohkwmq/LJysoydteIiIjIQLj3kYlwdXXF8ePHVdK4NygREVHj0ehGyk6ePIlly5Zh9OjRcHFxgUwmg5WVVY315HI5Fi1aBG9vb1hZWaFt27aIjIzEjRs3JOlXkyZN0KZNG5VPkyZNJDk2ERERmb5GN1IWExODbdu2aVVHLpdj4MCBSE5OhrOzM0aNGoWsrCzExcVh+/btOHLkCDw9PfXq161bt+Dm5gYhBLp27YqoqCj07t1br2MSERFR/dHogrLg4GD4+fkhMDAQgYGBaNOmTY11YmNjkZycjODgYOzatQu2trYAgJUrV2L27NmIjIzEgQMHlOXv37+P+/fvV3vM5s2bw8nJCQDQs2dPfP/99/Dx8cHDhw/x7bffol+/fvjvf/+LwYMH63G2REREVF9Ivs3S7du3cfHiRXTq1AmtW7dWpmdmZmL+/Pk4e/Ys2rVrh0WLFiEoKEjKpnUik8lgaWkJuVyuNr+srAxOTk7Iy8vDqVOn4O/vr5Lv5+eHtLQ0nDhxAgEBAQAqJ+4vXry42nYjIiKq3dHg2WefhZWVFXbt2qXdCT2B2yzVEX23WQJ0r//+n4BMpltdc2vd6xIRNTLG+A2VfKRs2bJl+Oyzz3Du3DllUFZQUIC+ffvi1q1bEELg3LlzSEpKwpkzZ9ChQwepuyCpQ4cOIS8vD56enlUCMgAYO3Ys0tLSkJCQoAzKFi5ciAULFlR7XDOz6qfzBQUF4ddff9W949QwlRUDn3jpVpcLzxIRmTTJJ/rv378fnTt3RqdOnZRp8fHxuHnzJsaPH4+LFy9i1apVKCwsxCeffCJ185I7c+YMAKBHjx5q8xXpinJAZcDVtGnTaj81BWWpqalwc3OT6CyIiIjI1Ek+Uvbnn3+iV69eKmnbt29H06ZNsXr1ajg4OODtt9/G2rVrkZiYKHXzksvOzgZQuWSFOop0RTldzJ49G2FhYfDw8MDDhw/xzTffIDExsdYvJCiGWP8qIyND7xcQSA1z6/89htSlblmRtP3RhhC6t8/Hn0REdUryoCw/Px92dnbK70IIpKSkICAgAA4ODsr0Tp06Yfv27VI3L7mCggIAgLW1tdp8GxsblXK6yMnJwcSJE5Gbmwt7e3t07doVe/bswYABA3Q+JtUhmaz+PgbkvptERCZL8qDMxcUFmZmZyu8nTpzAw4cPERISolLu8ePHsLCwkLp5ySneg5BpGCGQ4j2J9evX61Vf0yRETSNoREREZHokD8qCg4Oxfv16bNu2DaGhofjwww8hk8kwYsQIlXLnz5+Hi4uL1M1LTjHqV1hYqDa/qKjyUZBimQyiGunz+FPal6WJiMiESB6UffDBB/jll18wevRoAJUjSaGhoSoLoWZlZeHcuXN49dVXpW5ecu7u7gCgceV+RbqiHFGN9Hn8War+HwdERFT/SR6U+fj44NChQ1i9ejVyc3MREBCAOXPmqJT5/fff4efnh/DwcKmbl5yfnx8A4NSpU2rzFendunUzWJ+IiIio4amTFf39/f2rXRh12rRpmDZtWl00Lbk+ffrA3t4eGRkZSE1NrbJW2aZNmwAAYWFhxugeERERNRB1tiH548ePsXXrVnzwwQeYNm0a1qxZo8zLyclBWloaHj9+XFfNS8bCwgLTp08HAEyfPl1lbtnKlSuRlpaGvn37IjAw0FhdJCIiogagTkbKDhw4gIkTJ+LPP/+EEAIymQxlZWWIjIwEAOzduxeTJ0/Gxo0bMWbMmLrogkY7duxATEyMSlppaanK2mpRUVEYPny48vuCBQuwZ88eJCcno2PHjujXrx+uXbuGlJQUODg4IC4uzmD9JyIiooZJ8pGys2fPYtiwYbhz5w7efvtt/Pzzz1WWjRgzZgysra2xefNmqZuvUW5uLlJSUpQf4H9rqSk+ubm5KnWsrKyQmJiIqKgoWFtbY+vWrcjKykJERARSU1Ph5aXjtjdERERE/0fykbIlS5agpKQEu3bt0rj4qbW1NTp37ozU1FSpm6/R5MmTMXnyZK3rNWvWDEuWLMGSJUuk7xQRERE1epKPlB04cAC9evWqcTV6d3d35OTouFYTERERUQMj+UjZo0eParUobElJCcrLy6Vunqhhq8/7bhIRUbUkD8qcnZ1x/vz5Gsv98ccfaNeundTNEzVs9XnfTSIiqpbkjy+fe+45pKenY8uWLRrLxMfH49q1aypvOBIRERE1ZpKPlM2fPx8//fQTxo8fjzlz5mDUqFEAKveI/OOPP7B161bExsbCwcEB77zzjtTNE1F19H38SUREdUYm/rpehQQOHTqEF198Ebdv34ZMJlPJE0LA0dERW7ZsUdkPk6Tn6+sLAEhPTzdyT4iIiOoXY/yG1snisX379sWlS5fw3XffYc+ePcjKykJ5eTlcXV0xaNAgTJs2DS1atKiLpomIiIjqpToZKSPTwJEyIiIi3RjjN7TO9r4kIiIiotpjUEZERERkAiSfU5aUlKRV+WeffVbqLhARERHVO5IHZSEhIVXeuKwOV/UnIiIiqoOgbNKkSWqDsoqKCly/fh2nTp3Co0ePMHLkSLRs2VLq5omIiIjqJcmDsvj4+Grz79+/j6lTp+LcuXM4cuSI1M0TERER1UsGn+jfqlUr/PDDD3j48CHef/99QzdPREREZJKM8valtbU1goKC8OuvvxqjeSLShRBAaaFuHy6HSERUozpZ0b82CgoK8ODBA2M1T0TaKisCYtvqVnd+DmBhI21/iIgaGKOMlCUkJCApKQne3t7GaJ6IiIjI5Eg+UhYZGakxr6CgAJcuXcLZs2chhMDs2bOlbp6IiIioXjL425cA4O7ujujoaEyaNEnq5omIiIjqJcmDssTERI15FhYWcHZ2hoeHh9TNEhEREdVrkgdl/fv3l/qQRERERA0eNyQnIiIiMgF6j5R16NBB57oymQwZGRn6doGIiIio3tM7KMvKypKgG0RERESNm95BWUVFhRT9ICIiImrUjLaiPxE1MkJU7gqgC3NrQCaTtj9ERCaGQZkJCAkJwYEDB6qky2Qy3Lp1C05OTkboFZHE9Nmm6f0/dQ/KGNARUT3BoMwE/PLLLygtLVVJGzJkCBwdHRmQEQFAWTHwiZdudbnvJhHVE3USlAkh8OOPP2Lbtm24fPky8vPzIYSoUs4Yb1+ePHkSu3fvxrFjx5CSkoKcnBxYWlpCLpdXW08ul2Pp0qVYv349srOz0apVKwwZMgRLliyBq6urXn1q1aqVyvfTp0/jzJkz2Lx5s17HJSIiovpD8qCstLQUw4cPx759+9QGYkBlMKYpr67FxMRg27ZtWtWRy+UYOHAgkpOT4ezsjFGjRiErKwtxcXHYvn07jhw5Ak9PT8n6+NVXX6Ft27YYOXKkZMck0pu5deWok651dZ1PRkTUSEi+eOyKFSuwd+9ehIWF4fLly5g4cSJkMhlKSkpw/vx5REdHw8bGBnPmzDHKm5vBwcFYuHAhEhIScOvWrVrViY2NRXJyMoKDg3Hp0iVs2LABKSkpWLFiBXJzc6tswn7//n1cuXKl2s+dO3fUtvXo0SOsW7cOr732Gpo25dNlMiEyWeVjQF0+nNNFRFQjyX/1N2zYgFatWmHdunWwsbGBmVll3Gdubo5OnTph4cKFCA0NRWhoKDp16lQloKlr8+bN06p8WVkZPv/8cwDAF198AVtbW2XeO++8g7Vr1yIpKQknT55EQEAAAOCzzz7D4sWLqz1uRESE2s3b//Of/6C4uBivvfaaVv0kIiKi+k3ykbIrV64gKCgINjaVE2sVQVl5ebmyTL9+/dCnTx98+eWXUjcvuUOHDiEvLw+enp7w9/evkj927FgAQEJCgjJt4cKFKCsrq/azZs0ate19/fXXGDlyJFxcXOrmhIiIiMgkSR6UNWnSBM2bN1d+VwRnubm5KuVcXFxw8eJFqZuX3JkzZwAAPXr0UJuvSFeUAyoD0aZNm1b7UQSrTzp48CD++OMPvPHGG3VwJkRERGTKJH986eLiguzsbOV3L6/K19iPHj2K8PBwZXpaWprKo0BTpTgXTW9YKtKfPGddffXVV/Dy8sKgQYO0qufr66s2PSMjQ9IXEIiIiKjuSD5S1qtXL6Snp6O4uBgAMGzYMADA22+/jZ07d+Ls2bOYMWMGzp8/j549e0rdvOQKCgoAANbW1mrzFSOBinK6ys3NxebNmzFt2jTIOCmaiIio0ZF8pGzMmDHYuXMndu/ejZEjR8LLywszZ87EqlWrEBYWBqByHTMbGxv885//lLp5ySmW7tAUKEm1tIejoyNKSkp0qpuenq42XdMIGhEREZkeyYOy4cOH4+bNmyppK1asQGBgILZu3YoHDx7A29sbb731Fjp27Ch185Kzs7MDABQWFqrNLyqqXHupPjyKJSIiItNlsIWwXn75Zbz88suGak4y7u7uAIAbN26ozVekK8oRkQniZuhEVA9IHpRFRUXhlVdeQefOnaU+tFH4+fkBAE6dOqU2X5HerVs3g/WJiLSkz2bo3DuTiAxE8qDso48+QmxsLPz9/TFhwgSMGzcObdvq+JehCejTpw/s7e2RkZGB1NTUKmuVbdq0CQCU8+WISAN9tmky0rZsRESGJHlQtmzZMqxbtw6nTp1Camoq5s6di5CQEEyYMAGjR49WWcOsPrCwsMD06dPx0UcfYfr06di1a5fyjcuVK1ciLS0Nffv2RWBgoJF7SmTiFNs06UII7rtJRA2eTNTRzuAXLlzADz/8gJ9++glXr16FTCaDpaUlwsLCMGHCBAwbNgzm5uZ10XS1duzYgZiYGOX3lJQUyGQyBAUFKdOioqIwfPhw5Xe5XI6QkBCkpKTA2dkZ/fr1w7Vr15CSkgIHBwccPXpUuR6bKVG8fanp7UyiRqO0kI8viUgrxvgNlXydMgUfHx98+OGHuHLlCo4cOYI333wT9vb22LRpE0aPHo3WrVvj73//e101r1Fubi5SUlKUH6ByWYsn0/66+4CVlRUSExMRFRUFa2trbN26FVlZWYiIiEBqaqpJBmRERERUv9TZSJk6FRUV2Lt3L+Li4vDTTz9BJpOp7IlJ0uJIGdH/4UgZEWmpQY2UqZOUlISff/4Zv//+uyGbJSIiIjJ5db5O2enTp/Hjjz/ip59+Qk5OjnI1/1deeQUTJkyo6+aJiPQjROVImy64xhkRaaFOgrLMzEysW7cOP/74Iy5evAghBJo2bYohQ4ZgwoQJCA8P17iXJBGRSSkrBj7Rcd4oH30SkRYkD8p69+6NlJQU5Z6QvXr1Uq5X9tRTT0ndHBEREVGDIHlQdvToUfj4+CgfT7Zv317qJoiIiIgaHMmDshMnTqBHjx5SH5aIiIioQZP87UsGZERERETaM+iSGERERESkXp0viUFEZHTcDJ2I6gEGZUTU8OmzGbqua5QREWmJjy+JiIiITACDMiIiIiITwKCMiIiIyARIPqfsnXfeQcuWLREVFSX1oYmI6h8hgLIi3epy70yiRkXyoOz//b//h1GjRkl9WCKi+qmsCIhtq1td7p1J1KhI/vjS1dUVFRUVUh+WiIiIqEGTfKTshRdewNq1a5Gfnw87OzupD09EZFj6rHFmbq37o0sianQkHymLjo6Gu7s7hg0bhtTUVKkPT0RkWIo1znT5cD4YEWlB8pGyUaNGwdLSEocPH8YzzzwDZ2dnuLu7w8rKqkpZmUyGvXv3St0FIqKGQQjdF6/lSwJE9Y7kQdn+/fuV/18IgZycHOTkqB/6l/EvDCIizcqKgU+8dKvLlwSI6h3Jg7LMzEypD0lERETU4EkelLVr107qQxIRERE1eFzRn4iIiMgESD5SBlTOJfvxxx+xbds2XL58Gfn5+RBCVCknk8mQkZFRF10gIiIiqlckD8pKS0sxfPhw7Nu3T20gBlQGY5ryiIiIiBojyR9frlixAnv37kVYWBguX76MiRMnQiaToaSkBOfPn0d0dDRsbGwwZ84crvxPRERE9H8kHynbsGEDWrVqhXXr1sHGxgZmZpVxn7m5OTp16oSFCxciNDQUoaGh6NSpEyIjI6XuAhEREVG9I/lI2ZUrVxAUFAQbm8r1cRRBWXl5ubJMv3790KdPH3z55ZdSN09ERERUL0k+UtakSRM0b95c+V0RnOXm5qJNmzbKdBcXFyQkJEjdfL1TVFSEjz76CBs2bMCNGzfg4eGBd955B3//+9+N3TUikoI+e2dy7i1RoyJ5UObi4oLs7Gzldy+vytWojx49ivDwcGV6WloabG1tpW6+3nn99deRlJSEr776Cp06dcKhQ4fw+uuvw9zcHFOmTDF294hIX4q9M3Wh6xZLRFQvSf74slevXkhPT0dxcTEAYNiwYQCAt99+Gzt37sTZs2cxY8YMnD9/Hj179pS6+WqdPHkSy5Ytw+jRo+Hi4gKZTKZ2T86/ksvlWLRoEby9vWFlZYW2bdsiMjISN27c0Ks/crkc69evR0xMDIYOHYoOHTpg0qRJmDp1KpYsWaLXsYmIlHtn6vLhKB2RwUk+UjZmzBjs3LkTu3fvxsiRI+Hl5YWZM2di1apVCAsLA1C5jpmNjQ3++c9/St18tWJiYrBt2zat6sjlcgwcOBDJyclwdnbGqFGjkJWVhbi4OGzfvh1HjhyBp6enTv0pKytDeXk5mjVrppJubW2NrKwsZGdnw93dXadjExGhrAiIbatbXe6dSWRwkgdlw4cPx82bN1XSVqxYgcDAQGzduhUPHjyAt7c33nrrLXTs2FHq5qsVHBwMPz8/BAYGIjAwUGWOmyaxsbFITk5GcHAwdu3apXzkunLlSsyePRuRkZE4cOCAsvz9+/dx//79ao/ZvHlzODk5wc7ODn379sVHH32E7t27w9PTE0eOHMF3330HAPjzzz8ZlBERETUSMtGIV3GVyWSwtLSEXC5Xm19WVgYnJyfk5eXh1KlT8Pf3V8n38/NDWloaTpw4gYCAAABAdHQ0Fi9eXG27ERERiI+PBwBcv34dr732Gnbt2gUzMzO4uLjglVdewbJly5CSkoKgoCCdz8/X1xcAkJ6ervMxiMiIhKgc7dKFuTVHyoj0YIzfUO59WY1Dhw4hLy8Pnp6eVQIyABg7diwAqLxFunDhQpSVlVX7WbNmjbK8m5sb/vvf/6KoqAhZWVnIzMyEq6srAKBDhw51fIZEZNIULwno8pHJjN17ItJSnex92VCcOXMGANCjRw+1+Yp0RTmgcl02xdps2rCysoKrq6ty39CQkBA89dRTtaqriOb/KiMjQ+f5bkRERGRYegdl+ozmmPqG5IqlPRQjV3+lSH9yCRBt7d27F0VFRfD19cXNmzexfPlypKen49ChQzofk4iIiOofvYOyrKwsCbphmgoKCgBUvg2pjmJhXEU5XeTn52PevHm4du0amjVrhtDQUBw5cgRdunSp9TE0Pe/WNIJGREREpkfvoKwhbyqueAdCpmFuhhTvSISHh6ssqktERESNE+eUVcPOzg4AUFioflXtoqLKt6K4MwERNTiKhWd1YW7NFw2IdMCgrBqKNcI0rdyvSOdaYkTU4JQVA5946VaXy2kQ6UTyoCwpKUmr8s8++6zUXZCMn58fAODUqVNq8xXp3bp1M1ifiIiIqGGSPCgLCQnROAdLnfLycqm7IJk+ffrA3t4eGRkZSE1NrbJW2aZNmwBAuX0UERERka4kD8omTZqkNiirqKjA9evXcerUKTx69AgjR45Ey5YtpW5eUhYWFpg+fTo++ugjTJ8+Hbt27VK+cbly5UqkpaWhb9++CAwMNHJPiYiIqL6TPChTbB+kyf379zF16lScO3cOR44ckbr5au3YsQMxMTEqaaWlpejVq5fye1RUFIYPH678vmDBAuzZswfJycno2LEj+vXrh2vXriElJQUODg6Ii4szWP+JiOoNfbeI4osC1AgZfKJ/q1at8MMPP8DT0xPvv/8+vv76a4O1nZubi5SUFJU0IYRKWm5urkq+lZUVEhMTsXTpUqxbtw5bt25Fy5YtERERgZiYGLi5uRmk70REWjO3rpx0rwt9l/zhvptEWjPahuSjRo3C8ePHkZOj418YVCNuSE5EOist1C+oAhiUUb3WqDYkLygowIMHD4zVPBEREZFJMUpQlpCQgKSkJHh7exujeSIiIiKTI/mcssjISI15BQUFuHTpEs6ePQshBGbPni1180RERET1ksHfvgQqV8CPjo7GpEmTpG6eiIjqO27xRI2U5EFZYmKixjwLCws4OzvDw8ND6maJiKih0HeLJ3NrLsdB9ZLkQVn//v2lPiQREVHtcTkOqqeM9vYlEREREf2PwRePJSKiekCfhWf1eXxI1IhJHpQlJSXVqpyFhQUcHBzg5eWl1QbmRERkADIZH+MRGZjkQVlISIhWQZaNjQ1efvllfPTRR3B0dJS6O0RERLXHNz/JiCQPyiZNmoS8vDz8+uuvMDMzg7+/P9zd3SGEwPXr15GamoqKigqMGDECxcXFOH36NP79739j3759yk2+iYiIjELfNz85ukh6kHyi//Lly5GWlobhw4fj0qVLOH78ODZv3oxffvkFx48fx+XLlxEWFoa0tDT85z//QVZWFiZMmICrV69i+fLlUneHiIiIqF6QfEPyyMhI7NmzB1euXIGFhYXaMiUlJejYsSMGDhyIuLg4FBYWon379nBwcMD58+el7E6jxg3JichohNB9sr8QwFIX3erquxn6u1c4UkYAjPMbKvnjy99++w0DBgzQGJABgKWlJfr06YOdO3cCqJxX5u/vj0OHDkndHSIiMgZ9XhTQdU4XUT0n+ePLR48e4e7duzWWu3fvHh49eqT83rJlS76FSURERI2W5CNlvr6+2L9/P5KTk9G7d2+1ZY4cOYLExER0795dmXbjxg04OTlJ3R0iIiLD0eexLd/ebPQkD8rmzJmDl19+GYMHD8bkyZMxduxYuLm5QSaTITs7G5s2bUJ8fDwqKiowZ84cAJWjaydPnsQLL7wgdXeIiIgMh1s8kR4kD8peeukl5OTk4L333sNXX32Fr7/+WiVfCAFzc3MsX74cL730EoDKR5nR0dEYOHCg1N0hIiIiqhfqZJulmTNnYtSoUfjuu+9w5MgR5ORUvg3j7OyM3r17Y8qUKfD09FSWb9++PebNm1cXXSEiovqGWzxRI1Vne1+2b98eH374YV0dnoiIGip9t3jSJ6iTdpUoIq1wQ3IiImpYuBwH1VN1FpQVFRXhxIkTuHnzJkpKSjSWmzRpUl11gYiIiKjeqJOgbOHChVi1ahWKijQ/1xdCQCaTMSgjIiIiQh0EZcuXL8eHH36Ipk2bIiwsDN7e3rC1tZW6GSIiIqIGRfKg7F//+heaNWuGgwcPokePHlIfnoiIqO7wzU8yIsmDsuvXr2PAgAEMyIiIqP7R981PIj1IvvdlmzZtpD4kERERUYMneVD28ssv49ixY3jw4IHUh66XkpKSMGrUKLRr1w4ymQzR0dF6lSMiIqKGSfKgLDo6Gl26dMHo0aORkZEh9eHrnYKCAnTp0gXLly+vdhSxtuWIiIioYZJ8TtmwYcNQUVGBw4cPw8fHBx4eHnB1dYVMJqtSViaTYe/evZK1ffLkSezevRvHjh1DSkoKcnJyYGlpCblcXm09uVyOpUuXYv369cjOzkarVq0wZMgQLFmyBK6urnr1adiwYRg2bBgAVLuVVG3LERERUcMkeVC2f/9+5f8vLy9HRkaGxhEzdYGaPmJiYrBt2zat6sjlcgwcOBDJyclwdnbGqFGjkJWVhbi4OGzfvh1HjhxR2aeTiIhII33f3qRGTfKgLDMzU+pD1lpwcDD8/PwQGBiIwMDAWj0GjI2NRXJyMoKDg7Fr1y7lmmorV67E7NmzERkZiQMHDijL379/H/fv36/2mM2bN4eTk5N+J0NERPUP394kPUgelLVr107qQ9aato/9ysrK8PnnnwMAvvjiC5VFbt955x2sXbsWSUlJOHnyJAICAgAAn332GRYvXlztcSMiIhAfH69d54mIiKhRk3yif31y6NAh5OXlwdPTE/7+/lXyx44dCwBISEhQpi1cuBBlZWXVftasWWOwcyAiIqKGoc42JK8Pzpw5AwAaF7pVpCvKAYCZmRnMzBp1LEtERER1QO+grEOHDpDJZNizZw/at2+PDh061LquTCYz6rIZ2dnZAKDxDUtFuqKcLgoKCnDlyhUAQGlpKW7duoXTp0/DwsICXbp00bqcOr6+vmrTMzIy+JICEVF9IYTu2zSZW1fOZ6N6Te+gLCsrC0Dl/Kwnv9cHBQUFAABra/VvvNjY2KiU08WJEycQGhqq/P7NN9/gm2++Qbt27VSuVW3LERFRA1VWBMS21a3u/By+YNAA6B2UVVRUVPvdlAkhAGhemkORr4+QkJBaHae25dRJT09Xm65pBI2IiIhMT6OeHGVnZwcAKCwsVJtfVFQ5jPzkW5lEREREdUHyoOz27dtISkrC7du3VdIzMzMxfvx4PP300xg+fDhSUlKkblpr7u7uAIAbN26ozVekK8oRERER1RXJg7Jly5YhNDQUeXl5yrSCggL07dsXGzduxLlz57Bz504MGjQIV69elbp5rfj5+QEATp06pTZfkd6tWzeD9YmIiIgaJ8mDsv3796Nz587o1KmTMi0+Ph43b97E+PHjcfHiRaxatQqFhYX45JNPpG5eK3369IG9vT0yMjKQmppaJX/Tpk0AgLCwMEN3jYiISDtCAKWFun0kmENN+pN8nbI///wTvXr1Uknbvn07mjZtitWrV8PBwQFvv/021q5di8TERKmb14qFhQWmT5+Ojz76CNOnT8euXbuUb1yuXLkSaWlp6Nu3LwIDA43aTyIiohrx7c16T/KgLD8/XzmBHqh8gzElJQUBAQFwcHBQpnfq1Anbt2+XtO0dO3YgJiZGJa20tFQlSIyKisLw4cOV3xcsWIA9e/YgOTkZHTt2RL9+/XDt2jWkpKTAwcEBcXFxkvaRiIiISB3JgzIXFxeVTclPnDiBhw8fIiQkRKXc48ePYWFhIWnbubm5VV4gUASFT5Z5kpWVFRITE7F06VKsW7cOW7duRcuWLREREYGYmBi4ublJ2kciIiKTo3j0qQsuXCsZyYOy4OBgrF+/Htu2bUNoaCg+/PBDyGQyjBgxQqXc+fPn4eLiImnbkydPxuTJk7Wu16xZMyxZsgRLliyRtD9ERET1Qlkx8ImXbnX56FMykgdlH3zwAX755ReMHj0aQOVIVWhoKHr37q0sk5WVhXPnzuHVV1+VunkiIiIyNG4RJQnJgzIfHx8cOnQIq1evRm5uLgICAjBnzhyVMr///jv8/PwQHh4udfNERERkaHzJQBKSB2UA4O/vj/j4eI3506ZNw7Rp0+qiaSIiIqJ6qU6CMqr/hBCS7P1J1NjIZDKN++kSEVWHQRkplZeX4969e8jPz0dpaamxu0NUbzVp0gTW1tZo3rw57OzsGKRR7ZhbVz7K07WurnO6yGQwKCMAlQFZdnY25HK5sbtCVO+Vl5cjPz8f+fn5aNGiBVq3bg0zM8k3UKGGRibj3KpGjkEZAQDu3bsHuVyOJk2aoHXr1rCxseGPCJEOhBAoKSlBfn4+7t+/j7y8PFhZWaFly5bG7hoRmTgGZQSgcicGAGjdujXs7e2N3Bui+s3a2hrW1tZo2rQp7ty5gwcPHjAoI9KEC9cqMSgjCCGUc8gUe38Skf6aN2+OO3fuoKSkBEIIzi0jUocL1yrx+RSpvGXJR5ZE0mnSpIny//NtZiKqCUfKiIiIGgJ93t6sz/9oaEC7CTAoIyIiagj0eXtT1zldpqAB7SbAZ1VEREREJoAjZUT1WFZWFtq3b4/+/ftj//79xu4OEdVXXLjWJHCkjIhMmkwmg4eHh2THy8rKgkwmQ0hIiGTHJKr3FI8+dfmY0Jys+o4jZWRwufkl2HA8GymZ91FQ8hi2lk3Rq4MDXnrGDY52lsbuHhERaauxvmQgMQZlZDDysnIsTkjHppM3UFau+h/hwct38emeSxgb4IZFI7rAyryJhqMQEZHJaawvGUiMjy/JIORl5YhYcwzrj12vEpAplJULrD+WjYg1xyAvKzdwD1U9+Yjr0aNHmD17Ntq3bw9zc3PMnDkTAJCXl4fPP/8czz//PNq1awdLS0s4ODhgyJAh2L17t9rjhoSEQCaTISsrC1u3bkWvXr1gY2ODVq1aYfz48bhx44baenfv3sW0adPQpk0bWFtbw9/fH99//3215/D48WN8/vnnCAgIgK2tLWxtbREUFISvvvoK5eVVr++TfduwYQMCAwNhbW0NFxcXzJ07V7nAcEZGBsaPHw8nJydYW1tjwIABSEtL0+LqVm7rNX/+fPj6+sLW1hb29vbw9vbGpEmTcOzYMQBAfHy8crHVa9euQSaTKT9PPno8ffo05s6di4CAADg6OsLS0hIdOnTAP/7xD+TkqP7LPTo6Gu3btwcAHDhwQOWYkydPVimbm5uLd999F506dVJukzR06FAkJSVpda5ERLXFoIwMYnFCOlIy79eqbErmfSxOOFfHPaqd4uJi9O/fH3FxcejevTtGjhyp3C7n6NGjeOutt3D+/Hl07NgRL7zwAjp16oRdu3bh+eefx5o1azQe98svv8SYMWMghMCQIUNga2uLn376CQMGDEBxcbFK2Xv37qFPnz749ttvYW5ujpEjR8Le3h5TpkzBxx9/rPb45eXlGDVqFN566y1cuXIFgwYNwqBBg3DhwgX84x//wIsvvoiKigq1dVevXo2//e1vMDc3x/PPP4/S0lJ8/PHHeO2113D58mX06tULx44dQ79+/eDl5YXExESEhobi9u3btbqmBQUF6NWrF5YuXYqysjI8//zzGDRoEOzt7bF+/Xr89ttvAAAvLy9EREQAqNxpIiIiQvkZMmSI8njLli3DypUrUV5ejj59+mDYsGEQQuCrr77CM888oxKYde/eHWPGjAFQuaXYk8fs27evstyFCxfg7++PFStWoLy8HMOGDUO3bt2wb98+hIaGYt26dbU6VyIibcgEl5lusHx9fQEA6enp1ZarqKjAxYsXAQCdOnWSfFX/O/ly9Fm2T+MImTrmTWRIfm+g0eaYKd5qBIDg4GD89ttvaNGihUqZzMxM3Lx5E71791ZJT01NxYABA1BRUYE///wTtra2yryQkBAcOHAANjY2+PXXXzFgwAAAQFFREQYPHozk5GR89913iIyMVNaZNm0avv32W4waNQobNmyApWXlNdm5cydGjhyJx48fV3n7csWKFXj33XfRtWtX7NmzB05OTgCAmzdvIjQ0FBcvXsQXX3yBf/zjH1X6Zmdnhx07dqBfv34AgFu3bqF79+64c+cOfHx8MHToUHz88ccwMzODEAKTJ0/G999/j4ULF2Lx4sU1Xtv4+HhMmTIFM2bMwGeffaaSd+fOHdy5cwdPP/20Mk0mk6Fdu3bIyspSe7x9+/ahc+fOcHZ2VqZVVFTgww8/xKJFizBlyhSVALmmN1bLy8vRvXt3/PHHH1i9ejVmzJihHLFLTU3F4MGDIZfLcfXqVeV11aSu/9siahBKC/VbZwyok3XKavsbKiX+DUF1buNxzY8sNSkrF9h44nod9Ug7n332WZWADADat29fJSADAH9/f7z55pt49OgREhMT1R5z1qxZyoAMqNzAevbs2QCg8nisoKAA//nPf9C0aVN89tlnyoAMAIYOHYoXX3xRY58B4NNPP1UJHJydnZWja38NiJ7smyIgA4A2bdpgwoQJyj1S//nPfyqDC5lMpuz3gQMH1B7vr+7cuQMAKuev4OTkpBKQ1caAAQNUAjKgcruwhQsXwsXFBdu2bdPqeAkJCfjjjz8wfvx4vPXWWyr7Vfr7+yMqKgqFhYX44YcftDouEVFNONGf6lxtH1v+1dGr9/BmqI6b1ErE2dkZzzzzjMb88vJy7N27F8nJybh16xbkcjkA4PLlyyr/+1fPPfdclTRvb28AlaNZCqdOnUJxcTH69OkDd3f3KnXGjx+P9evXq6RlZ2cjOzsbbdq0URv4hIWFoUWLFrh48SJyc3Ph6Oiokj948OAqdTp06ACgcjStaVPVvzY8PT2r9Ls6AQEBAID58+ejadOmGDRoEKysrGpVV5N79+7h119/xR9//IG8vDzlnLmysjLcv38f9+/fR6tWrWp1LMV8wPDwcLX5isecx48f16vPRGQChND8ooFQP8WjLjEoozpXUPLYoPWkpC4QUrhx4wbCwsJw5swZjWXy8/PVpru6ulZJUzzmLCkpUaYp5kNp6oe6dEUdTWt7KR4H5uXlIScnp0pQ5uLiUqWOjY1NjXlP9rs6AwcOxKxZs/Dpp59ixIgRsLCwQPfu3fHcc8/h1Vdf1XpNsvXr1+Pvf/87CgoKNJbJz8+vdVCmeEw6btw4jBs3TmO5u3fvatVPItLAmAvXlhUDn2j4x/9dzX+n1BUGZVTnbC11u810rSel6kZwpk6dijNnzmD06NGYN28eOnXqBDs7O5iZmeHbb7/FtGnToGnK5pOPxKqjqF/b8tq2oa5MdfV06Yc6K1euxLRp07Bt2zbs3bsXhw8fxrFjx7B8+XJs2LBB4yjVX127dg2TJ0+GEAKffvophg8fDhcXFzRr1gwA0Lt3bxw5ckTjn4M6ilG2oUOHVjtnzMfHp9bHJKJq6LOcRgNj/F89avB6tm+Fg5e1H1Xo1cGhDnojjcLCQuzevRutW7fGxo0b0aSJ6rpqV69elaSdtm0rJ69eu3ZNbX52drbGOpmZmRqPq6j317lYhtSpUyfMnTsXc+fOhVwuxxdffIF3330X06ZNq3VQ9ttvv6G0tBSzZ8/G22+/XSVflz8HxSjm66+/jpEjR2pdn4hIVwzKqM69FOiG1Xsva/325UvPuNVhr/Tz8OFDVFRUwNnZuUpA9vjxY2zZskWSdgICAmBlZYWUlBRcv34dbm6q1+Snn36qUsfd3R3u7u7Izs7Gvn37qswr27FjBx48eIBOnTpVeXRpLFZWVpg9ezZWrFiBmzdv4s6dO8pRKnNzczx+rP5R9oMHDwCgynUBKl+YULdMh4WFBQBoPOagQYPw3XffYevWrQzKiOqDBrSbAN++pDrnZGeFsQFV51BVZ2yAaW+55OTkBHt7e/zxxx84fPiwMr28vBxz587FpUuXJGnH1tYWEyZMwOPHj/H222+rzNvatWsXNm7cqLbejBkzAFS+SZmbm6tMv3XrFubMmaNSxtC2bt2Ko0ePVklPTU3F7du3YWdnp1wLDqgc+bt9+zby8vKq1FG8HPHDDz+gsPB/k3X//PNPvP7662rbf+qpp2Bubo6MjAy1i+iOHTsWPj4+iI+Pxz//+U+UlZWp5JeWluKXX37B2bNna3W+RFTHGtC+nQzK6lBSUhJGjRqFdu3aQSaTITo6Wm25TZs2ISgoCK1atYKVlRW8vLywYMEC5QrqDcGiEb7o2b52E617tm+FRSO61HGP9NO0aVPMnTtXuUbYc889h5dffhleXl74+uuv8eabb0rW1rJly+Dl5YUtW7agY8eOGD9+PAYMGIChQ4fitddeU1tn1qxZGDp0KNLS0tCxY0eMHj0aL7zwAry9vXH+/HmEh4fjjTfekKyP2ti/fz+Cg4Ph6uqKESNGYMKECQgNDUVQUBAqKioQExMDc3NzZXnFWmw9evTA3/72N0ydOlW5rMfIkSPh6+uLEydOwMvLC2PHjkVYWBi8vb3RsmVLtUuWWFhYYMiQIbh16xb8/PwwadIkTJ06FXFxcQAq/2y3bNkCNzc3vPfee2jXrh2GDBmCl156CcHBwWjdujXGjBmDjIwMw1wwImo0GJTVoYKCAnTp0gXLly9HmzZtNJZr1aoV5s6di4MHD+LChQv4+OOP8e233+Ldd981YG/rlpV5E6yNDML4IHeYN1H/LxPzJjKMD3LH2sigerH35fz587F27Vp069YNhw8fxp49e+Dn54ejR49Wu4yGtp566ikcPnwYU6dORUlJCbZu3Yp79+7hX//6F+bOnau2TpMmTfDrr79i9erV6NChA37//Xfs2rULnTp1whdffIFNmzYZbSHTyZMnY/bs2Wjbti2OHTuGzZs3IzMzE8OGDUNiYmKVuWFLly7F9OnT8fjxY2zYsAHfffcdduzYAaAywDp48CDeeOMNWFlZYfv27Th//jxmzJiB3bt3qwR3T/r3v/+NiRMn4t69e1i3bh2+++47lXXWfHx8cPr0aURHR8PJyQmHDh3Cjh07kJubi2effRZxcXEYNGhQ3V0kImqUGsyK/idPnsTu3btx7NgxpKSkICcnB5aWlsp1ozSRy+VYunQp1q9fj+zsbLRq1QpDhgzBkiVL1C5boCsPDw9MnjxZ42jZX82aNQt79uzR6xGJqazo/1e5+SXYeOI6jl69h4KSx7C1bIpeHRzw0jOm/ciSSFtc0Z/IxFWzm4Dvl5VLYqTfMdxezA1mon9MTIzWK3fL5XIMHDgQycnJcHZ2xqhRo5CVlYW4uDhs374dR44cUS6MaUjnz5/Hzp071S4w2hA42lnizVAvoy8MS0REZEoaTFAWHBwMPz8/BAYGIjAwsNrHhQqxsbFITk5GcHAwdu3apVy8c+XKlZg9ezYiIyNVHmkoVgavTvPmzWvcD08TW1tblJWVobS0FNOmTcOqVat0Og4RERHVPw0mKJs3b55W5cvKyvD5558DAL744guVTaPfeecdrF27FklJSTh58qRyW5jPPvusxg2XIyIiEB8fr13n/8/p06dRXFyMEydO4L333kPr1q1rtcEzERER1X8NJijT1qFDh5CXlwdPT0/4+/tXyR87dizS0tKQkJCgDMoWLlyIBQsWVHtcfeaMeHlVPs7r2rUrzMzMEBkZiblz5yq3sSEiIqKGq9EGZYr9Cnv06KE2X5H+5L6GZmZmBp2oK4SoskYSERERNUyNNihTbDOj6Q1LRbq6bWxqq6CgAFeuXAFQueDkrVu3cPr0aVhYWKBLl/+twxUTE4OePXuiQ4cOEELg2LFjmDdvHkaNGoUWLVrU2I7iLcu/ysjIMMqLCkRERKS9RhuUFRRUvupqbW2tNl/xyFBRThcnTpxAaGio8vs333yDb775Bu3atUNWVpYyvbi4GDNmzMD169fRtGlTeHh4YNasWXjrrbd0bpuIiIhqUN0WTRulW2+ythptUKZYnk2mYYsFKZZvCwkJqdVxYmNjERsbq3M7mtYh0zSCRkRERPjfFk1q8wy/rmCjXcnQzs4OAFT2y3tSUVERAKi8lUlERERUVxptUObu7g4AuHHjhtp8RbqiHBEREVFdarRBmZ+fHwDg1KlTavMV6d26dTNYn4iIiKjxarRBWZ8+fWBvb4+MjAykpqZWyd+0aRMAICwszNBdIyIiokao0QZlFhYWmD59OgBg+vTpKnPLVq5cibS0NPTt2xeBgYHG6iIRERE1Ig3m7csdO3YgJiZGJa20tBS9evVSfo+KisLw4cOV3xcsWIA9e/YgOTkZHTt2RL9+/XDt2jWkpKTAwcEBcXFxBus/ERERNW4NZqQsNzcXKSkpyg9QuazFk2m5ubkqdaysrJCYmIioqChYW1tj69atyMrKQkREBFJTU5XbHpEehABKC3X7SLAsCUkjPj4eMpkM0dHRxu4KEVGD1WBGyiZPnozJkydrXa9Zs2ZYsmQJlixZIn2nCCgrAmLb6lZ3fo7m9WOowdi/fz9CQ0MRERGB+Ph4SY4ZHx+PKVOmYNGiRQwkiajeaDAjZURERET1GYMyIiIiIhPAoIyoGhs3bkRgYCCaNWuG1q1bY8qUKbh9+zYmT54MmUyG/fv3q5SXyWTw8PBAaWkplixZAh8fH1haWiI8PFxZ5vr165g2bRratWsHS0tLODk5YfTo0Th+/HiV9vfv3w+ZTKbx0XxN/SgvL8fy5cvh7e0NS0tLuLm5Yd68eSgpKVF7vLS0NISFhcHe3h729vYYPHgwjhw5os0lUzp//jwmTpwIT09PWFlZwdHREd27d8fMmTNx8+ZNZf8V+8OuXbsWMplM+XnyseOOHTsQGRmJzp07o3nz5rCxsYGfnx9iY2OrnEtISAimTJkCAFi8eLHKMf/6ePTs2bOYMGECXFxcYGlpibZt22LKlCkqe9MSERlKg5lTRiS1Tz/9FLNmzUKTJk0QEhKCp556Crt378b+/furXVS4oqIC4eHhSEpKQv/+/dGtWzc4ODgAqAwCBgwYgLt378LHxwejR49GdnY2tmzZgoSEBKxbtw4vvviiZOcwYcIEbN++HUFBQejUqRMOHjyI5cuX488//8QPP/ygUjYlJQUDBgxAUVERunfvDh8fH/zxxx/o37+/1vM1T506hb59+0IulyMoKAhBQUHIz8/H1atXsXr1aoSHh8PZ2Rl9+/bFrVu38Pvvv8PT0xN9+/ZVHqN79+7K///qq6+isLAQvr6+6Nq1Kx49eoRjx47hgw8+wN69e7Fr1y40adIEADBkyBA8fvwYhw8fhp+fn8pxnnx5Z/PmzXjllVdQWlqKgIAA9O7dGxkZGYiPj0dCQgIOHDjA/WOJyLAENVhdunQRXbp0qbFceXm5OHfunDh37pwoLy+XthMlBUIsaq7bp6RA2r5oISMjQ1hYWAgrKyuRlJSkTC8uLhbDhw8XAAQAkZiYqFJPke7l5SVu3LihkldRUSG6du0qAIj3339fVFRUKPN+/vlnYWZmJuzs7MStW7eU6YmJiQKAiIiIUNvPiIiIavvRuXNnkZmZqUy/evWqaNmypQAgrly5okwvLy8XPj4+AoBYunSpyrEWLFigPN6iRYuquWpV+7V58+YqeefOnRM5OTm1PkchhNiyZYsoKFC9Hx49eiTCwsIEALF27VqVvLi4uGr7e/XqVWFtbS3s7e3FgQMHVPLWrl0rAIjAwMAazrJmdfrfFhHVqdr+hkqJjy+J1FizZg1KS0sRERGBfv36KdOtrKywevVqmJlV/5/O0qVL4eLiopK2f/9+nD17Fu3bt0dMTAxkMpkyb+zYsQgPD0d+fr6k6+N9/vnn8PDwUH5v3749/va3vwEADh48qNK3CxcuwNvbG/PmzVM5xqJFi7TeA/bOnTsAgAEDBlTJ69y5M5ydnbU6Xnh4OGxsVN/EtbOzw6pVqwAA27Zt0+p4q1evRlFREZYvX45nn31WJW/SpEkIDw/H8ePHNW7DRkRUFxiUEamRnJwMAGofJXp6esLf319jXZlMhhEjRlRJVwRB48aNUz5qe9LEiRNVyunL3NwcISEhVdK9vb0BQDmvCwAOHToEoPJ8nwwWAaBp06YYO3asVm0HBAQAqAxwjh07hoqKCq3qq3P58mWsXr0aM2bMQGRkJCZPnqxcMPry5ctaHWv37t0AgFGjRqnNVzxGVTfPj4iornBOGZEaOTk5AAA3Nze1+e7u7jh58qTaPCcnJ1haWmo85pMjV09SpCvK6cvZ2Vlt8GdrawsAKhPkFW1qGhHTdqRszpw5OHToEBISEpCQkAB7e3v07NkTYWFhmDx5Muzs7Gp9LCEE3n33XaxatQpCw4LC+fn5WvVPMZG/TZs21Za7e/euVsclItIHgzKiavx11EhBU3AAVD7i1OWYtc1/UnUjUNocR3E+2tSpTvPmzbFv3z4cPnwYCQkJ2L9/v3JC/tKlS3Hw4EF4enrW6lgbNmzAypUr4erqik8//RTBwcFwdHSEubk5SktLYWlpWe2fhzrl5eWQyWSYNGlSteU40Z+IDIlBGZEazs7OuHjxIrKzs9GxY8cq+devX9f6mG3bVu5skJmZqTb/2rVryrYVLCwsAAAFBQVq6+jSj+r6pujDX2VnZ2t9TJlMhr59+yofBebm5uLtt9/G+vXrMX/+fGzYsKFWx9myZQsA4KuvvkJYWJhK3tWrV7XuFwC4uroiIyMDn332GZo3b67TMYiIpMY5ZURq9O7dGwCwadOmKnlXr15Famqq1sdUvDCwYcMGlJeXV8lXLFHx5IsFigDt0qVLVcrfu3dPsonoisBp8+bNVUadHj9+jM2bN+vdhqOjo3LtsbNnzyrTFYHn48eP1dZ78OABAPWPkjdu3Ki2Tk3HHDRoEABg69atNXeciMhAGJQRqTFlyhSYm5sjPj5eOekfAORyOWbOnKnTxPWQkBB07doVmZmZWLhwoUrws3XrVvzyyy+wtbVVWROsffv2cHd3x9mzZ1XeMCwsLMRrr72GR48e6XaCfxEaGgpvb29cuHABn3zyiUrehx9+qHEETZOvv/5a7Yjgzp07AajOUVOM0l28eFHtsRQvJnz77bcq1+zgwYP4+OOP1dap6ZizZ89Gs2bNMGvWLCQkJFTJv3//Pr788ksUFxerrU9EVBcYlBGp4eXlhdjYWMjlcjz77LMYPHgwXn75ZXTs2BFnzpxRvl2pGJGpDZlMhh9//BEODg6IjY2Fr68vXnnlFfTt2xcvvPACzMzMsGbNmiqTzxWjS2PGjMGAAQMwcuRIeHp64o8//sDIkSMlOV8zMzPEx8ejWbNmmDt3Lnr06IFXXnkF3bp1Q2xsLKZOnarV8b7++mt06NABvr6+GDt2LF5++WX4+/tj5syZaNasGRYtWqQs6+HhgW7duuHEiRMICgrClClTMHXqVPz6668AgLfeegs2Njb48ssv8fTTT2P8+PF49tln0b9/f7z++utq2+/VqxecnJywadMmhISEIDIyElOnTlUG2B07dsQPP/yA4uJijBw5Ej4+PnjhhRcQHh4Of39/ODs7480339S48wERUZ0w6KpoZFAmsXhsRUXlIrC6fJ5YXNVY1q9fL3r06CEsLS2Fo6OjmDhxosjJyRGDBg0SAMSFCxdUygMQ7dq1q/aY165dE6+99ppwc3MT5ubm4qmnnhLh4eEiJSVFY524uDjx9NNPCwsLC9G6dWsxdepUcffu3WoXj9XUj+oWVk1NTRVDhw4VdnZ2ws7OTgwYMEAcOnSoxsVY/+rXX38VkZGRwtfXV7Ro0UJYW1sLb29v8fe//11cvny5SvnLly+L8PBw4eDgIMzMzKq0de7cOTFixAjh5OQkrK2thb+/v/j222+rPdfjx4+LwYMHC3t7eyGTyQQAERcXp1Lm0qVLYtq0aaJDhw7C0tJS2Nvbi86dO4spU6aI7du3qyzwqwsuHktUfxlj8ViZEFq+tkT1huLNsfT09GrLVVRUKB/zdOrUqcaFURu7wsJCeHh4oLi4GA8fPlS77AQRwP+2iOqz2v6GSol/QxBpcPXqVTx8+FAlraCgAK+//jru3r2rcRFYIiIiXXBJDCINNm7ciOjoaAQEBMDV1RUPHjxAamoq7t69Cw8PD8TGxhq7i0RE1IAwKCPSYODAgTh9+jSOHj2K1NRUCCHg7u6OiIgIzJs3D46OjsbuIhERNSAMyog0CAwMxE8//WTsbhARUSPBOWVEREREJoBBGREREZEJYFBGREREZAIYlBFkMpny/+uyfRARqffkHqdP/ndGRKQOgzKCTCZTbhdUWFho5N4QNRyKvUktLS0ZlBFRjfj2JQEA7OzscO/ePdy+fRsAYGNjw9XHiXQghEBJSQny8/Nx//59AEDLli2N3Csiqg8YlBEAwMHBAYWFhZDL5cjJyTF2d4gajBYtWsDe3t7Y3SCieoBBGQEAmjRpAnd3d9y7dw/5+fkoLS01dpeI6q0mTZrAxsYGdnZ2sLOz46NLIqoVBmV1LCkpCStWrMDp06eRnZ2NRYsWITo6ukq56OhoLF68uEp6ZmYmPDw86r6jqPwhcXJygpOTE4QQ4F71RNqTyWQMwohIJwzK6lhBQQG6dOmCV155BTNnzqy2rKurK44fP66SZqytfPjDQkREZFgNaib3yZMnsWzZMowePRouLi6QyWSwsrKqsZ5cLseiRYvg7e0NKysrtG3bFpGRkbhx44befRo2bBiWLl2KcePGwdLSstqyTZo0QZs2bVQ+TZo00bsPREREZPoa1EhZTEwMtm3bplUduVyOgQMHIjk5Gc7Ozhg1ahSysrIQFxeH7du348iRI/D09KyjHqu6desW3NzcIIRA165dERUVhd69exukbSIiIjKuBhWUBQcHw8/PD4GBgQgMDESbNm1qrBMbG4vk5GQEBwdj165dsLW1BQCsXLkSs2fPRmRkJA4cOKAsf//+feVr7po0b94cTk5OWvW9Z8+e+P777+Hj44OHDx/i22+/Rb9+/fDf//4XgwcP1upYREREVP/IRAOezS2TyWBpaQm5XK42v6ysDE5OTsjLy8OpU6fg7++vku/n54e0tDScOHECAQEBADRPyH9SREQE4uPjq6R7eHhg8uTJaif6q/Pss8/CysoKu3btqlX5v/L19QUApKen61SfiIiosTLGb2iDmlOmrUOHDiEvLw+enp5VAjIAGDt2LAAgISFBmbZw4UKUlZVV+1mzZo0k/QsKCkJWVpYkxyIiIiLT1qAeX2rrzJkzAIAePXqozVekK8oBgJmZmcFWuk9NTYWbm5tB2iIiIiLjatRBWXZ2NoDKpSjUUaQryumioKAAV65cAQCUlpbi1q1bOH36NCwsLNClSxdludmzZyMsLAweHh54+PAhvvnmGyQmJtbqxQXFEOtfXbhwAebm5hrziYiISL2MjAyYm5sbtM1GHZQVFBQAAKytrdXm29jYqJTTxYkTJxAaGqr8/s033+Cbb75Bu3btVB5N5uTkYOLEicjNzYW9vT26du2KPXv2YMCAATq3XVFRgbKyMp3rN1YZGRkAYLC3bhsCXjPd8Lppj9dMN7xu2isrK8Pjx48N2majDsoU7zhoWiRVincgQkJCanWc9evX69yGpkmInOivG1437fGa6YbXTXu8ZrrhddOeMZ4yNeqJ/nZ2dgCAwsJCtflFRUUAoFwmg4iIiKiuNOqgzN3dHQA0rtyvSFeUIyIiIqorjToo8/PzAwCcOnVKbb4ivVu3bgbrExERETVOjToo69OnD+zt7ZGRkYHU1NQq+Zs2bQIAhIWFGbprRERE1Mg06qDMwsIC06dPBwBMnz5dZW7ZypUrkZaWhr59+yIwMNBYXSQiIqJGokFts7Rjxw7ExMQov6ekpEAmkyEoKEiZFhUVheHDhyu/y+VyhISEICUlBc7OzujXrx+uXbuGlJQUODg44OjRo/Dy8jLoeRAREVHj06CWxMjNzUVKSopKmhBCJS03N1cl38rKComJiVi6dCnWrVuHrVu3omXLloiIiEBMTAxX1CciIiKDaFAjZURERET1VaOeU0ZERERkKhiUEREREZkABmVEREREJoBBGREREZEJYFBGREREZAIYlDVAcrkcixYtgre3N6ysrNC2bVtERkZq3OOTgJCQEMhkMo2f//73v8buolGcPHkSy5Ytw+jRo+Hi4gKZTAYrK6sa633//fcICgqCra0tWrVqhWHDhiE5OdkAPTYN2l636Ojoau+/9957z4C9N46ioiJs3boVr776Krp164bmzZvDxsYGfn5+WLJkCQoKCjTWbaz3my7XjPdapZUrV2L06NHo2LEj7O3tYWlpiXbt2iEiIgLp6eka69X1vcYlMRoYuVyOgQMHIjk5WbkYblZWFo4dOwZHR0ccOXIEnp6exu6myQkJCcGBAwcwZswY2NraVsmfPXs2unbtaoSeGVd4eDi2bdumkmZpaQm5XK6xzjvvvINVq1ahWbNmeO655yCXy7F3714IIfDzzz/jhRdeqOtuG5221y06OhqLFy9Gnz591C5WPXz4cLz44ot10ldT8e9//xuvvfYaAMDX1xddunTBo0ePkJycjPz8fPj4+ODAgQNwcnJSqdeY7zddrhnvtUpPPfUUCgsL0a1bN7i4uAAA0tPTcenSJVhYWGDr1q0YOnSoSh2D3GuCGpSoqCgBQAQHB4v8/Hxl+ooVKwQA8eyzzxqxd6arf//+AoDIzMw0dldMyrJly8TChQtFQkKCuHXrlgAgLC0tNZbfu3evACAcHBzEpUuXlOnJycnCwsJC2Nvbi/v37xui60al7XVbtGiRACDi4uIM10kTs3btWvHGG2+o3DdCCJGTkyP8/f0FADF+/HiVvMZ+v+lyzXivVTp06JAoLi6ukv7ll18KAKJt27bi8ePHynRD3WsMyhqQ0tJS0aJFCwFAnDp1qkp+t27dBABx4sQJI/TOtDEoq52agothw4YJAGLVqlVV8t566y0BQHzyySd12EPTxKBMP8nJycprWFJSokzn/aaZpmvGe61mXl5eAoBIT09XphnqXuOcsgbk0KFDyMvLg6enJ/z9/avkjx07FgCQkJBg6K5RI6AYygf+d689ifcf6crPzw8AUFJSgnv37gHg/VYTddeMaqdJkyYAAAsLCwCGvdca1N6Xjd2ZM2cAAD169FCbr0hXlKOqvvvuO9y7dw9mZmbw9vZGeHg43N3djd2teuHChQsoKSmBo6MjXF1dq+Qr7r+0tDRDd63e2LdvH06fPg25XA5XV1cMHToUAQEBxu6W0V29ehUAYG5ujlatWgHg/VYTddfsSbzX1Pv+++9x8eJFeHt7o0OHDgAMe68xKGtAsrOzAUDtTfNkuqIcVfXhhx+qfH/33XcRFRWFqKgoI/Wo/qjp/rOxsUGLFi3w4MED5Ofnw87OzpDdqxf+85//qHyPiorCmDFjEB8fr/YFlMZi9erVAIAhQ4bA0tISAO+3mqi7Zk/ivVbp448/Rnp6OgoLC3H+/Hmkp6ejbdu2WLduHczMKh8mGvJe4+PLBkTx+rO1tbXafBsbG5Vy9D/PPvss/vOf/yAjIwNFRUW4ePEiPvroIzRt2hQLFy5U/gVHmtV0/wG8BzXx8vLCJ598gvT0dBQUFOD69ev48ccf4eLigs2bN2PixInG7qLR/Pbbb/juu+9gbm6OmJgYZTrvN800XTOA99pf/f7771i7di02bdqE9PR0uLm5Yd26dSqjhga91/SelUYmY+rUqQKAWLBggdr8S5cuCQDC29vbwD2rv37//XcBQNjb24uioiJjd8foUM2E9R9++EEAEH379tVYv23btgKAyMnJqasumqTqrlt1cnJyhIODgwAgDh8+XAc9M23nzp0TLVu2FADEp59+qpLH+0296q5ZdRr7vfbgwQORlJQkBg0aJACIDz/8UJlnyHuNI2UNiGLItLCwUG1+UVERADSqoWl9Pffcc3jmmWfw8OFDHD161NjdMWk13X8A70FtOTs7Y8qUKQAq/0XfmNy4cQNDhgzBgwcP8M477+Dtt99Wyef9VlVN16w6jfleA4AWLVqgX79++O233xAQEICoqCgcP34cgGHvNQZlDYhiQrqmlfsV6Zy4rp2OHTsCAG7evGnknpi2mu6/wsJC5OXloUWLFo1ufo8+GuP9d/fuXQwePBjZ2dmYMmUKPvnkkypleL+pqs01q0ljvNf+ytzcHOPGjYMQQvk2pSHvNQZlDYjiFehTp06pzVekd+vWzWB9aggePHgAoPH8a1tXnTp1gqWlJXJzc9X+5cX7TzeN7f7Lz8/H0KFDceHCBYwePRr/+te/IJPJqpTj/fY/tb1mNWls95omTz31FAAgNzcXgGHvNQZlDUifPn1gb2+PjIwMpKamVsnftGkTACAsLMzQXau3cnNzcfDgQQCalxqhSs2aNcOAAQMA/O9eexLvP+0JIbBlyxYAaBTLFZSUlGDUqFE4ceIEnn/+eaxfv165ZtRf8X6rpM01q05ju9eqc+DAAQBQbklo0HtNrxlpZHI++OADAUD07t1bFBQUKNMV2yxVN1GxsTpy5IjYt2+fqKioUEnPzMwUffr0EQDEyJEjjdQ704IaJqzv3r1b41YklpaWonnz5uLevXuG6KpJqe665ebmirVr1wq5XK6Snp+fL6ZNmyYAiDZt2ojCwkJDdNVoHj9+LF544QUBQPTr169W59vY7zdtrxnvtUpJSUnip59+EmVlZSrppaWl4rPPPhNmZmaiWbNmIjs7W5lnqHuNG5I3MHK5HCEhIUhJSVFuSH7t2jWkpKTAwcEBR48eVbsJbWMWHx+PKVOmwNnZGd7e3mjTpg1u3LiBkydPQi6Xw9fXF/v27auyEXJjsGPHDpVX6lNSUiCTyRAUFKRMi4qKwvDhw5XfZ86cidWrV8Pa2hqDBw9GaWkpdu/ejYqKCmzcuBFjxowx6DkYgzbXLSsrC+3bt0fz5s3RuXNnuLu7Iy8vD6dOncK9e/fQokULbN++HX369DHGqRjM6tWrMXPmTADACy+8gObNm6st98knnygfLwGN+37T9prxXquk+Dv/qaeeQkBAABwcHHD37l2cPXsWN2/ehJWVFdauXYuXXnpJpZ5B7jW9wzoyOUVFRSIqKkp4enoKCwsL0bp1axEREaES9dP/nDt3TrzxxhuiR48ewtHRUTRt2lTY29uLXr16iRUrVjTqpTDi4uIEgGo/6vbQi4uLEwEBAcLa2lrY29uL559/Xhw8eNDwJ2Ak2ly3R48eiXnz5on+/fsLFxcXYWlpKaytrYWvr6+YPXu2uHHjhnFPxkAUezLW9FG3P21jvd+0vWa81ypdvXpVzJ8/X/Tp00c4OzsLc3NzYWNjI3x9fcWMGTPE5cuXNdat63uNI2VEREREJoAT/YmIiIhMAIMyIiIiIhPAoIyIiIjIBDAoIyIiIjIBDMqIiIiITACDMiIiIiITwKCMiIiIyAQwKCMiIiIyAQzKiIiIiEwAgzIiIiIiE8CgjIiIiMgEMCgjIiIiMgEMyoiIiIhMAIMyogZOJpPBw8PD2N3QWn3tty5M9VwN1a/JkydDJpNh//79dd4WkSljUEZERERkApoauwNEVLfOnz8Pc3NzY3eD6iHeO0SGxaCMqIHz8fExdheonuK9Q2RYfHxJVE+dP38eEydOhKenJ6ysrODo6Iju3btj5syZuHnzprKcunlBWVlZkMlkCAkJQXFxMd577z20a9cOlpaW8PLywj//+U8IIdS2u3HjRgQGBqJZs2Zo3bo1pkyZgtu3b2ucF7R//37IZDJMnjxZ7fG0nU+0Y8cOREZGonPnzmjevDlsbGzg5+eH2NhYlJSUaDzPR48eYfbs2Wjfvj3Mzc0xc+bMats5efIkZDIZevXqpbHM8uXLIZPJ8MEHH+jUP030vWZZWVmYNm0aPDw8YGlpCUdHR4wdOxZpaWm1al9B6ntn8+bNCAoKUt47kyZNQk5OTrV90OZcxo0bB5lMhnnz5lXJu3DhAqytrdG8eXNcvXq19hcBwG+//YbnnnsO9vb2cHBwwIgRI5CZmQkAWLp0KWQyGVJTU7U6JpFagojqnZMnT4pmzZoJmUwmevbsKV5++WUxfPhw0blzZwFAJCYmKssCEO3atVOpn5mZKQCI4OBg0bdvX9GyZUvx/PPPi+eff15YWVkJAOKDDz6o0u6qVasEANGkSRMxcOBAMW7cOOHi4iI8PDzEyJEjq7QthBCJiYkCgIiIiFB7LhEREWrrqeu3EEK0bt1a2Nraip49e4oXX3xRPP/886Jly5YCgBgwYIB4/PhxlfMMCgoS3bt3Fy1bthTh4eFi9OjRIjo6urpLLIQQwsfHRwAQV65cUZvv5+cnAIg//vhDp/5pOlddr5kQQhw8eFA0b95cABC+vr5i7NixIjg4WMhkMtGsWTOxb9++Gs9bU7+E0P3e+fzzz5X3zoABA8RLL70knJ2dhZubmwgLC1N7Ptqey/3794Wrq6swMzNTOVZpaanw9/cXAERcXFytz7+iokK8/vrrAoBwdnYWL730kvD19RUARM+ePUVRUZFwcnISI0eOrPUxiarDoIyoHlL8KG/evLlK3rlz50ROTo7ye3U/rABEv379RG5urjLv+PHjomnTpsLa2lrk5+cr0zMyMoSFhYWwsrISSUlJyvTi4mIxfPhw5fHqOijbsmWLKCgoUEl79OiR8od97dq1as8zODhYPHjwQG0fNFmyZIkAIJYsWVIl79y5cwKA6N69u879U5AqKHv48KFo06aNMDc3Fz///LNK3u7du4WFhYVwcXERJSUl1Zy15n4Jodu9k5mZKSwtLYWlpaVKnwsLC8XgwYPV3ju6nsu+ffuEmZmZcHNzU/55z507VwAQY8eOrdV5K0RHRwsAIjQ0VPlnWlZWJrp37y4AiEWLFgkA4sSJE1odl0gTBmVE9dDQoUMFgFoFGdX9sJqZmYmLFy9WqTNixIgqP5IffPCBACCmTZtWpfyVK1eEmZmZQYIyTS5fviwAiNGjRyvTngwgjh8/XutjKWRkZAgAolOnTlXyFNfj448/1rl/ClIFZYqRzPfff19tvZkzZ2oM5tWR6t6JiooSAMRrr71WpfyFCxeETCarUkefc5kzZ44AIMaPHy8SExOFmZmZaNu2rbh3716tzlsIIW7evCksLS2FlZWVyMzMVMmbN2+eACCsrKxEWFhYrY9JVBNO9CeqhwICArBz505MmjQJCxYswDPPPAMzM+2niHp4eMDb27tKuiLtyblpycnJAIAXX3yxSnlPT0/4+/vj5MmTWvdBF5cvX8Zvv/2GK1euoLCwEBUVFcp5TJcvX65S3tnZGc8884zW7XTo0AG9evXC0aNHcerUKfTo0UOZ99NPP8HMzAwvv/yy3v2Tyu7duwEA4eHhavP79u2LTz/9FMePH8fo0aP1akube+fQoUMAgJdeeqlK+U6dOsHf3x+nTp1SSdfnXD788EPs2bMH69evx44dOyCEwNq1a9GqVatan98PP/yAkpISvPHGG1Xm1bVo0QIAIJfLsXDhwlofk6gmDMqI6qE5c+bg0KFDSEhIQEJCAuzt7dGzZ0+EhYVh8uTJsLOzq9VxXF1d1abb2toCgMrEdMWEbDc3N7V13N3d6zwoE0Lg3XffxapVqzROJs/Pz1fbN11NmDABR48exY8//qgMyo4ePYqMjAyEhoaqXENd+yeVrKwsAEDPnj2rLXf37l2929Ll3tH05+Du7l4lKNPnXCwsLBAfHw8/Pz88evQIM2bMwKBBg6o9zl/VFBQCwNChQxEYGKjVcYmqw6CMqB5q3rw59u3bh8OHDyMhIQH79+/H3r17sWvXLixduhQHDx6Ep6dnjceRyWRat62pjqYgpCYVFRW1LrthwwasXLkSrq6u+PTTTxEcHAxHR0eYm5ujtLQUlpaWavthZWWlU9+Ayjf6Zs2ahZ9++gkff/wxzMzMsG7dOgCVAZsU/dOWpmtWXl4OoHI009raWmP9mgKd2tDm3lGcszZ19D2XDRs2KP9/amoqKioqtBpNVgSJvXv3rpL36NEjAMDcuXNrfTyi2mBQRlRPyWQy9O3bF3379gUA5Obm4u2338b69esxf/58lR8lKTg7O+PixYvIzs5Gx44dq+Rfv35dbT0LCwsAQEFBgdp8TfXU2bJlCwDgq6++QlhYmEqetssc1JajoyMGDx6MnTt3Yv/+/ejfvz82btwIS0tLjBkzpk76p+s1c3V1xcWLF7FgwQJ069at1u3VtbZt2+LSpUu4du2a2nsnOzu7Spo+53Lw4EEsW7YMbdu2hY+PD/bt24dly5Zh/vz5tapfUlKCu3fvws7OTjnyp3D//n18/fXXyj4SSYnrlBE1EI6OjoiOjgYAnD17VvLjK0YMNm3aVCXv6tWrGtdpcnZ2BgBcunSpSt69e/eqPLaqzoMHDwCof4S6cePGWh9HW4oRsXXr1mHv3r24ffs2hg8frpxbJHX/dL1mikd0W7durXVbhqD4h8PPP/9cJe/SpUs4ffp0lXRdz+XRo0eYOHEihBCIi4vDDz/8AAcHB0RHR+PEiRNaHauoqAilpaUqaQsWLFD+OTdp0kSr4xHVhEEZUT309ddfKxevfNLOnTsB6DeHSpMpU6bA3Nwc8fHxykn/QOVk55kzZ2p8pNa+fXu4u7vj7Nmz2LZtmzK9sLAQr732mvJRUG0oJpF/++23Ko8BDx48iI8//ljbU6q18PBw2NjYYPPmzYiLiwNQ9dGllP3T9ZpNmzYNjo6OiI2NRVxcXJVHpYWFhfj+++9x48aNWvdFClOmTIGFhQW+//57HDx4UJleXFyMt99+W+29o+u5/OMf/8C1a9cwY8YMPPfcc3B2dsa//vUvlJWVYcKECSgqKqqxv5aWlujQoQPKy8vx3XffKdM3btyIr7/+WjmSqe6/QSK9GP6FTyLSl2LR0i5duogxY8aIcePGKddOatasmUhOTlaWRTXLGvTv31/t8RXrL/11oc2PP/5YuQDooEGDxLhx44Srq6twd3dXLoVw+PDhKsdbs2aNsl5oaKgYMWKEaN26tejYsaPGRWfV9fvixYvCxsZGee4vv/yy6Nevn5DJZOLdd9+tUqem89TGK6+8olxew97eXsjl8ipltO1fdeeqyzUTQohDhw6JVq1aKY85fPhwMXr0aPHMM88o+5aamlqrc5by3lG38HDbtm2Fq6urxsVjtT2X9evXKxeaLS4uVjnWlClTNC7pos7KlSuVS38MGTJEDBo0SAAQPj4+YvXq1QKAcHV1FbNnz67V8Yhqg0EZUT3066+/isjISOHr6ytatGghrK2thbe3t/j73/8uLl++rFJWyh9WISp/+Hr06CEsLS2Fo6OjmDhxosjJyVH+aF24cEHtMePi4sTTTz8tLCwsROvWrcXUqVPF3bt3tV6n7Ny5c2LEiBHCyclJWFtbC39/f/Htt9+qrSNlULZjxw5lUBYZGamxnDb9U9CUru01U/jzzz/F7NmzhY+Pj2jWrJmwtbUV3t7eYty4cWLDhg2SLB6ry72zceNGERAQICwtLcVTTz0lXnnlFXHjxo1qz6e255KdnS1atGghLCws1Aad+fn5okOHDgKA2LZtW43nXlFRIVasWCG8vb2FhYWFsLe3F2PGjBF3794VcrlcTJw4UVhZWYlBgwbVeCyi2pIJIcGrQETUqBUWFsLDwwPFxcV4+PAh59oQEemAc8qIqNauXr2Khw8fqqQVFBTg9ddfx927dzFu3DgGZEREOuJIGRHV2rJlyxAdHY2AgAC4urriwYMHSE1Nxd27d+Hh4YGjR4+idevWxu4mEVG9xHXKiKjWBg4ciNOnT+Po0aNITU2FEALu7u6IiIjAvHnz4OjoaOwuEhHVWxwpIyIiIjIBnFNGREREZAIYlBERERGZAAZlRERERCaAQRkRERGRCWBQRkRERGQCGJQRERERmQAGZUREREQmgEEZERERkQlgUEZERERkAhiUEREREZkABmVEREREJoBBGREREZEJYFBGREREZAIYlBERERGZAAZlRERERCaAQRkRERGRCfj/2Px8jroiHzUAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 600x450 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(4,3))\n",
    "plt.plot(S_random, 'o', label='random state')\n",
    "plt.plot(S_ground_state, 's', label='ground state')\n",
    "plt.yscale('log')\n",
    "plt.ylim(1.e-15, 1.)\n",
    "plt.xlim(0, 30)\n",
    "plt.xlabel(r'singular value index $\\alpha$')\n",
    "plt.ylabel(r'singular values $\\lambda_\\alpha$')\n",
    "plt.legend(loc='lower left')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a87a7752-b8f7-4d0e-8ed7-7c2f56f770d8",
   "metadata": {},
   "source": [
    "Ultimately, this decay of singular values for ground states is the reason why DMRG and MPS are so highly successfull.\n",
    "It implies that we can truncate the state to much smaller bond dimension and still have large overlap with the initial state:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "id": "25b44bfd-6601-483b-bba3-08cfe7962404",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4.065772185590788e-08"
      ]
     },
     "execution_count": 49,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "chi = 5\n",
    "A_trunc = A_ground_state[:, :chi]\n",
    "S_trunc = S_ground_state[:chi]\n",
    "S_trunc = S_trunc / np.linalg.norm(S_trunc) # normalize to 1 again\n",
    "B_trunc = B_ground_state[:chi, :]\n",
    "psi_trunc = np.tensordot(np.diag(S_trunc), B_trunc, axes=[1, 0]) # a [a*], [a] j -> a j\n",
    "psi_trunc = np.tensordot(A_trunc, psi_trunc, axes=[1, 0]) # i [a], [a] j -> i j\n",
    "overlap = np.tensordot(np.conj(psi_trunc), psi_ground_state, axes=([0, 1], [0, 1])) # [i*] [j*], [i] [j]\n",
    "err = 1. - np.abs(overlap)\n",
    "err"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a5c56399-fc49-43c6-b05d-e7e8650beaae",
   "metadata": {},
   "source": [
    "### Exercise: truncation\n",
    "\n",
    "Play around with the \"bond dimension\" `chi` of Schmidt values kept above. What is the relation between the largest Schmidt value discared and the `err` above?\n",
    "\n",
    "Try truncating the random state. What is the (absolute value of the) overlap if you keep 32 of the original 64 Schmidt values?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "id": "a745617f-c84f-4982-991c-aa80a0771d02",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "9.037215420448774e-14"
      ]
     },
     "execution_count": 50,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "chi = 10\n",
    "A_trunc = A_ground_state[:, :chi]\n",
    "S_trunc = S_ground_state[:chi]\n",
    "S_trunc = S_trunc / np.linalg.norm(S_trunc) # normalize to 1 again\n",
    "B_trunc = B_ground_state[:chi, :]\n",
    "psi_trunc = np.tensordot(np.diag(S_trunc), B_trunc, axes=[1, 0]) # a [a*], [a] j -> a j\n",
    "psi_trunc = np.tensordot(A_trunc, psi_trunc, axes=[1, 0]) # i [a], [a] j -> i j\n",
    "overlap = np.tensordot(np.conj(psi_trunc), psi_ground_state, axes=([0, 1], [0, 1])) # [i*] [j*], [i] [j]\n",
    "err = 1. - np.abs(overlap)\n",
    "err"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "id": "fab290ec-c0a2-4063-b5d3-3aa78cda8bb2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.9439061265039328"
      ]
     },
     "execution_count": 51,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "chi = 32\n",
    "A_trunc = A_random[:, :chi]\n",
    "S_trunc = S_random[:chi]\n",
    "S_trunc = S_trunc / np.linalg.norm(S_trunc) # normalize to 1 again\n",
    "B_trunc = B_random[:chi, :]\n",
    "psi_trunc = np.tensordot(np.diag(S_trunc), B_trunc, axes=[1, 0]) # a [a*], [a] j -> a j\n",
    "psi_trunc = np.tensordot(A_trunc, psi_trunc, axes=[1, 0]) # i [a], [a] j -> i j\n",
    "overlap = np.tensordot(np.conj(psi_trunc), psi_random, axes=([0, 1], [0, 1])) # [i*] [j*], [i] [j]\n",
    "np.abs(overlap)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3fc14564-3750-4a1e-bc86-245f7e71e06d",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "344bccbb-e47a-4b63-9b69-17123b78422f",
   "metadata": {},
   "source": [
    "### Exercises (optional, if time left - for the experts, and those who want to become them)\n",
    "\n",
    "- What is the maximal entropy that one can represent with $\\chi$ Schmidt states? \n",
    "  (You can solve this analytically!)\n",
    "\n",
    "  \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "950535dd-32dc-488d-b584-1be6d4278776",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "id": "1fca3123-4735-4ae9-a296-ac92c6639768",
   "metadata": {},
   "source": [
    "- Find the half-chain entanglement entropy of the following states for differrent system sizes $L=6,8,10,12,14$ and\n",
    "  - random states as defined above\n",
    "  - ground states of the transferse field Ising model at $g = 1.2 J$, $g=0.8J$ and the critical point $g=J$.\n",
    "  \n",
    "  Plot curves of $S$ versus $L$. (Where) can you see the area law?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "70ae8fa8-a7f4-4daf-aea7-adf7e7310aab",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5e3a2b5c-166a-4d75-bcd3-d54a6e069217",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.10.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}