DMRG_functionalities.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Linear Algebra\n",
    "import numpy as np\n",
    "from scipy.linalg import svd\n",
    "from numpy.linalg import qr\n",
    "import scipy.sparse.linalg.eigen.arpack as arp\n",
    "\n",
    "# Python Plotting\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "# Own modules\n",
    "from DMRG import * # here is the main part of the DMRG code\n",
    "\n",
    "# we'll need them a lot\n",
    "sxx = np.array([[0., 1.], [1., 0.]])\n",
    "syy = np.array([[0., -1j], [1j, 0.]])\n",
    "szz = np.array([[1., 0.], [0., -1.]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# DMRG on the transverse nearest neighbor Ising model\n",
    "\n",
    "I implement the density matrix renormalization group algorithm (DMRG) for matrix product states (MPS).  \n",
    "Exemplarly I use the transverse nearest neighbor Ising model  \n",
    "$$\n",
    "H = - J \\sum_{i=0}^{L-2} \\sigma^x_i \\sigma^x_{i+1} - g \\sum_{i=0}^{L-1} \\sigma^z_i\n",
    "$$\n",
    "We compare for small number of sites $L<20$ to exact diagonalization of the Hamiltonian"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## MPS\n",
    "First we need to define an MPS class and initialize a random MPS\n",
    "This function is defined in <DMRG.py> and looks something like this  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Help on class random_MPS in module DMRG:\n",
      "\n",
      "class random_MPS(builtins.object)\n",
      " |  random_MPS(L, d, chi_max)\n",
      " |  \n",
      " |  Initializes a random, finite dimensional, unnormalized MPS\n",
      " |  \n",
      " |  \n",
      " |  Parameters\n",
      " |  ------------\n",
      " |  L: Number of Sites\n",
      " |  d: local Hilbert space dimension (Schollwöck: sigma_j = 1,..,d)\n",
      " |  chi: local bond dimension (Schollwöch: DxD matrices)\n",
      " |  \n",
      " |  attributes:\n",
      " |  Ms: list of L ndim=3 tensors M with indices Ms[j]: sigma_j, vL, vR\n",
      " |  Ss: list of L ndim=1 lists with the Schmidt values\n",
      " |  L: number of sites\n",
      " |  chi: list of bond dimensions, -chi(i-1)--Ms[i]--chi(i)--Ms[i+1]--..\n",
      " |  normstat = None, 'left-norm', 'right-norm'\n",
      " |  \n",
      " |  Methods defined here:\n",
      " |  \n",
      " |  __init__(self, L, d, chi_max)\n",
      " |      Initialize self.  See help(type(self)) for accurate signature.\n",
      " |  \n",
      " |  get_Mi(self, i)\n",
      " |  \n",
      " |  left_normalize(self, test=False)\n",
      " |      returns a left_normalized version of the random_MPS() object using QR\n",
      " |  \n",
      " |  right_normalize(self, test=False)\n",
      " |      returns a right_normalized version of the random_MPS() object\n",
      " |      \n",
      " |      new version using QR\n",
      " |  \n",
      " |  self_norm(self)\n",
      " |      calculate the norm of the MPS\n",
      " |  \n",
      " |  ----------------------------------------------------------------------\n",
      " |  Data descriptors defined here:\n",
      " |  \n",
      " |  __dict__\n",
      " |      dictionary for instance variables (if defined)\n",
      " |  \n",
      " |  __weakref__\n",
      " |      list of weak references to the object (if defined)\n",
      "\n"
     ]
    }
   ],
   "source": [
    "help(random_MPS)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The MPS basically consist of the list Ms that containts the ndim=3 arrays at each site. The legs are ordered $a_{i-1}, s_i, a_i$ (virtual left, physical, virtual right)  \n",
    "Sites are labeled from $0$ to $L-1$  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[(1, 3, 3),\n",
       " (3, 3, 5),\n",
       " (5, 3, 5),\n",
       " (5, 3, 5),\n",
       " (5, 3, 5),\n",
       " (5, 3, 5),\n",
       " (5, 3, 5),\n",
       " (5, 3, 5),\n",
       " (5, 3, 3),\n",
       " (3, 3, 1)]"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "randmps = random_MPS(L=10,d=3,chi_max=5)\n",
    "[randmps.Ms[i].shape for i in range(10)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[(1, 3, 3),\n",
       " (3, 3, 9),\n",
       " (9, 3, 27),\n",
       " (27, 3, 81),\n",
       " (81, 3, 243),\n",
       " (243, 3, 81),\n",
       " (81, 3, 27),\n",
       " (27, 3, 9),\n",
       " (9, 3, 3),\n",
       " (3, 3, 1)]"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "randmps = random_MPS(L=10,d=3,chi_max=500)\n",
    "[randmps.Ms[i].shape for i in range(10)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## MPO\n",
    "The MPO basically consists of the list Ws that contains the ndim=4 arrays at each site.  \n",
    "The indices are set as $b_{i-1}, b_i, \\sigma_i, \\sigma'_i$  \n",
    "For the transverse Ising model the matrices read\n",
    "$$\n",
    "W_{b_{i-1} b_i} = \\begin{pmatrix}\n",
    "\\mathbb{I} & \\sigma^x & - g \\sigma^z \\\\\n",
    "0 & 0 & - J \\sigma^x \\\\\n",
    "0 & 0 & \\mathbb{I} \\\\\n",
    "\\end{pmatrix}\n",
    "$$\n",
    "The left row vector is the top row and the right column vector is the right column."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Help on class Ising_MPO in module DMRG:\n",
      "\n",
      "class Ising_MPO(builtins.object)\n",
      " |  Ising_MPO(L, d, h, J)\n",
      " |  \n",
      " |  Transverse Ising field Hamiltonian\n",
      " |  H = -J \\sum_i \\sigma_i^x \\sigma_{i+1}^x - h \\sum_i \\sigma_i^z\n",
      " |  \n",
      " |  Methods defined here:\n",
      " |  \n",
      " |  __init__(self, L, d, h, J)\n",
      " |      Initialize self.  See help(type(self)) for accurate signature.\n",
      " |  \n",
      " |  init_W(self)\n",
      " |  \n",
      " |  ----------------------------------------------------------------------\n",
      " |  Data descriptors defined here:\n",
      " |  \n",
      " |  __dict__\n",
      " |      dictionary for instance variables (if defined)\n",
      " |  \n",
      " |  __weakref__\n",
      " |      list of weak references to the object (if defined)\n",
      "\n"
     ]
    }
   ],
   "source": [
    "help(Ising_MPO)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[(1, 3, 2, 2),\n",
       " (3, 3, 2, 2),\n",
       " (3, 3, 2, 2),\n",
       " (3, 3, 2, 2),\n",
       " (3, 3, 2, 2),\n",
       " (3, 3, 2, 2),\n",
       " (3, 3, 2, 2),\n",
       " (3, 3, 2, 2),\n",
       " (3, 3, 2, 2),\n",
       " (3, 1, 2, 2)]"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "isingmpo = Ising_MPO(L=10,d=2,h=0.5,J=1.0)\n",
    "[isingmpo.Ws[i].shape for i in range(10)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## DMRG\n",
    "Now I create an abstract class that acts as my DMRG engine and takes care of all the quantities needed for the algorithm, i.e. the MPS, the LP and the RP"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Help on class DMRG in module DMRG:\n",
      "\n",
      "class DMRG(builtins.object)\n",
      " |  DMRG(MPS, MPO, eps=1e-10, chi_max=None, test=False)\n",
      " |  \n",
      " |  Abstract DMRG engine class\n",
      " |  \n",
      " |  attributes:\n",
      " |  MPS*: the current MPS\n",
      " |  MPO: the model MPO\n",
      " |  RPs*: the current RPs\n",
      " |  LPs*: the current LPs\n",
      " |  (*): get regular updates\n",
      " |  \n",
      " |  LPs[i] is LP left site i, i.e. LPs[0] is dummy and LPs[-1] = LPs[L-1] is everything but the last\n",
      " |  RPs[i] is RP right of site i, i.e. RPs[-1]=RPs[L-1] is dummy and RPs[0] is everything but the first\n",
      " |  \n",
      " |  Parameters:\n",
      " |  MPS: object, matrix product state\n",
      " |  MPO: object, matrix product operator\n",
      " |  eps: float, epsilon determining the threshold of which singular values to keep, standard set to 10^-10\n",
      " |  chi_max: int, maximal bond dimension of system\n",
      " |  test: for debugging, prints things like Ms.shape, norm, singular values etc. during sweep\n",
      " |  \n",
      " |  Methods defined here:\n",
      " |  \n",
      " |  __init__(self, MPS, MPO, eps=1e-10, chi_max=None, test=False)\n",
      " |      Initialize self.  See help(type(self)) for accurate signature.\n",
      " |  \n",
      " |  left_to_right(self, eps=None, chi_max=None)\n",
      " |      Runs through the MPS from left to right\n",
      " |  \n",
      " |  next_LP(self, i)\n",
      " |      short version: LP[i] |-> LP[i+1] # Remember: LP[i] inlcudes everythings EXCEPT i -> need to add the contributions at site i\n",
      " |      \n",
      " |      extend LP from site i to the next site i+1 \n",
      " |      .------a_j-1----M[j]--a_j\n",
      " |       |               |\n",
      " |      LP[j]--b_j-1----W[j]--b_j\n",
      " |       |               |\n",
      " |      .-----a'_j-1---M*[j]--a'_j\n",
      " |      \n",
      " |      next_LP(i): takes LPs[i] and contracts with the (right-) following M[i] ,W[i] and M*[i]\n",
      " |      \n",
      " |      result: np.array[ndim=3] # a_i b_i  a'_i\n",
      " |  \n",
      " |  next_RP(self, i)\n",
      " |      short version: RP[i] |-> RP[i-1]\n",
      " |      \n",
      " |      extend LP from site i to the next site i+1 \n",
      " |      --a_j-1----M*[j]--a_j----.\n",
      " |                 |            |\n",
      " |      --b_j-1----W[j]--b_j---RP[j]\n",
      " |                 |            |\n",
      " |      --a'_j-1---M[j]--a'_j--.\n",
      " |      \n",
      " |      next_LP(i): takes RPs[i] and contracts with the (left-)following M[i] ,W[i] and M*[i]\n",
      " |      \n",
      " |      result: np.array[ndim=3] # a_i b_i  a'_i\n",
      " |  \n",
      " |  right_to_left(self, eps=None, chi_max=None)\n",
      " |      Runs through the MPS from right to left\n",
      " |  \n",
      " |  site_update(self, i)\n",
      " |      Solves the Eigenvalue problem H_eff v = E v\n",
      " |      \n",
      " |      Returns updated matrix M\n",
      " |  \n",
      " |  ----------------------------------------------------------------------\n",
      " |  Data descriptors defined here:\n",
      " |  \n",
      " |  __dict__\n",
      " |      dictionary for instance variables (if defined)\n",
      " |  \n",
      " |  __weakref__\n",
      " |      list of weak references to the object (if defined)\n",
      "\n"
     ]
    }
   ],
   "source": [
    "help(DMRG)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "init_mps = random_MPS(L=10,d=2,chi_max=1000)\n",
    "mpo = Ising_MPO(L=10,d=2,h=0.7,J=0.3)\n",
    "dmrg = DMRG(init_mps,mpo,eps=-1.,chi_max=1000,test=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A full sweep is performed by going left-to-right and then right-to-left again  \n",
    "The engine initially right-normalizes the state  \n",
    "(also note that the dmrg engine is manipulating the mps instance itself!)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "None\n"
     ]
    }
   ],
   "source": [
    "init_mps = random_MPS(L=10,d=2,chi_max=1000)\n",
    "print(init_mps.normstat)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "dmrg = DMRG(init_mps,mpo,eps=-1.,chi_max=1000,test=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "right-norm\n"
     ]
    }
   ],
   "source": [
    "print(init_mps.normstat)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After a left-to-right (l2r) sweep all matrices become left-normalized, the engine fixes the status after l2r"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "left-norm\n"
     ]
    }
   ],
   "source": [
    "dmrg.left_to_right()\n",
    "print(init_mps.normstat)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "(actually before doing a l2r or r2l sweep the engine checks the *normstat* of the MPS, pretty clever huh)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "ename": "AssertionError",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mAssertionError\u001b[0m                            Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-13-d8302ed417c1>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0mdmrg\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mleft_to_right\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m~/Dropbox/000 PhD/Machine-Tensor/Korb_Eman_Ang/DMRG_Binz_Code/DMRG.py\u001b[0m in \u001b[0;36mleft_to_right\u001b[0;34m(self, eps, chi_max)\u001b[0m\n\u001b[1;32m    312\u001b[0m         \u001b[0;31m# Pseudo code\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    313\u001b[0m         \u001b[0;31m# start from right-normalized state BBBB..\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 314\u001b[0;31m         \u001b[0;32massert\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mMPS\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mnormstat\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;34m'right-norm'\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    315\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mtest\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    316\u001b[0m             \u001b[0mprint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'the initial norm is'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mMPS\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mself_norm\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mAssertionError\u001b[0m: "
     ]
    }
   ],
   "source": [
    "dmrg.left_to_right()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "dmrg.right_to_left()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Some important subtleties:  \n",
    "LPs and RPs have length L while LPs[0] = RPs[-1] = ones(1,1,1) is the dummy. Therefore neither contains the full contraction at any time and LP[i] and RP[i] are EXCLUDING site i!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[(1, 1, 1), (2, 3, 2), (4, 3, 4), (8, 3, 8), (16, 3, 16), (32, 3, 32), (16, 3, 16), (8, 3, 8), (4, 3, 4), (2, 3, 2)]\n",
      "[(2, 3, 2), (4, 3, 4), (8, 3, 8), (16, 3, 16), (32, 3, 32), (16, 3, 16), (8, 3, 8), (4, 3, 4), (2, 3, 2), (1, 1, 1)]\n"
     ]
    }
   ],
   "source": [
    "print([dmrg.LPs[i].shape for i in range(10)])\n",
    "print([dmrg.RPs[i].shape for i in range(10)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Eigenenergies\n",
    "At each step in l2r or r2l the engine calculates a new, better matrix and then right- or left-normalizes it.  \n",
    "The *better* matrix is obtained from solving the Eigenvalue problem of the effective Hamiltonian at site i  \n",
    "$$\n",
    "\\sum_{\\sigma'_l,a'_{l-1},a'_l} \\left(H_\\text{eff}\\right)^{\\sigma_l \\sigma'_l}_{a_{l-1} a'_{l-1} a_l a'_l} M^{\\sigma'_l}_{a'_{l-1},a'_l} = \\lambda M^{\\sigma_l}_{a_{l-1},a_l}\n",
    "$$\n",
    "where\n",
    "$$\n",
    "\\left(H_\\text{eff}\\right)^{\\sigma_l \\sigma'_l}_{a_{l-1} a'_{l-1} a_l a'_l} = \\sum_{b_{l-1},b_l} L^{a_{l-1} a'_{l-1}}_{b_{l-1}} W^{\\sigma_l \\sigma'_l}_{b_{l-1} b_l} R^{a_l a'_l}_{b_l}\n",
    "$$\n",
    "\n",
    "which graphically is represented by  \n",
    "![visualization of eigenproblem](img/eigenproblem.png)  \n",
    "(see Schollwöck chap 6)  \n",
    "We are solving for the circled matrix M  \n",
    "The Eigenvalue $\\lambda$ on the other hand corresponds to the curent ground state energy and is saved (appended to) into ```dmrg.e```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "step_es = np.array(dmrg.e)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f92f50b9ef0>"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmMAAAGDCAYAAABnZBdiAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJzt3XmclXXd//H3ZzaGbXCAYRl2FFHGEGVQtDS31Cz3BcvU6lYz28u07n6VaXeapd7ZIplb3Vmi5VaZkiaoJAooKpvsyLA7INsAw8z5/P4418AB5gwDzDnfs7yej8d5zDnnus51va8v15z58L2u63uZuwsAAABhFIQOAAAAkM8oxgAAAAKiGAMAAAiIYgwAACAgijEAAICAKMYAAAACohgDcpiZnWBm74bOkU5mNtDM3MyKQmfJFmbW38w2mVlh6CxAPqIYA3KAmS02s9N2f9/dX3b3oSEyZRszu8nM/rgP859kZjUHsL6TzCwWFUGJj+P2d5n7y93fc/dO7t6Y7nUDkPifI4CUMrMid28InSNDLXf3vqlcgZmZJHP3WCrXA2D/0TMG5LDE3hsz+46Z/WW36b8ws7uj513M7H4zW2Fmy8zsx8kOW5lZezP7vZmtM7PZZnZDYi9R1FN3o5m9LWmzmRWZ2eFmNsHMPjCzmWZ2TsL8E8zsqoTXnzWzVxJeu5lda2bzonX+OioyZGaFZvZzM3vfzBZK+sRe2uTGaPs2mtm7ZnaqmZ0p6b8ljYl6p96K5v1ctH0bzWyhmX0her+jpH9Kqkzo0ao0s4KonReYWa2ZPWpmXVvzb9VMzglmdouZTYrWP97MuidMH21m/4na8y0zO2m3z/6PmU2SVCdpsJkNMrOXomU9H7XhH6P5dzm029K+YGaHmNlEM1sftfm4/dk+ADtRjAH548+SzjKzMilexEi6RNKfoum/l9Qg6RBJR0k6XdJVzSxHkn4oaaCkwZI+JukzzczzKcULo4MkmaS/SRovqYekr0h62Mz25RDqJyWNknRklPuM6P2ro2lHSaqWdFGyBUTr+7KkUe7eOVrGYnd/VtJPJI2LDtcdGX1kdbTsMkmfk3SXmR3t7pslfVzxnq1O0WO5pK9KOk/SRyVVSlon6df7sI27+3S03h6SSiRdH21HH0n/kPRjSV2j9/9qZhUJn71c0jWSOktaovi/8+uSukm6KZqeTEv7wi2K/zuWS+or6ZcHsH0ARDEG5A13XyLpDcWLBUk6RVKdu082s56KFxdfd/fN7r5a0l2SLk2yuEsk/cTd17l7jaS7m5nnbndf6u5bJI2W1EnSbe5e7+7/lvR3xQu21rrN3T9w9/ckvShpREKW/43WtVbSrS0so1FSO0nDzKzY3Re7+4JkM7v7P9x9gcdNVLwIOaGF5X9B0vfcvcbdtyle9FxkyS8mqIx6thIfHROmP+juc6M2fDRhmz8j6Rl3f8bdY+7+L0lTJZ2V8NmH3H1mdIi4t+KF7A+i9n9F0tPNBWrFvrBd0gBJle6+NVoWgANAMQbklz9pZwH0ae3sFRsgqVjSiqaiQNJvFe+RaU6lpKUJr5c2M0/ie5WSlu523tISSX32IfvKhOd1ihd3zWVZkmwB7j5f0tcVL5JWm9kjZlaZbH4z+7iZTTaztVGbnCWpe7L5FW/HJxLacLbiBWDPJPMvd/eDdntsTpiebJsHSLo4sYiT9BHFi64mu7f/WnevSzJ9921oaV+4QfGeztejw82fT7IcAK3ECfxAfnlM0h1m1lfS+ZKartxbKmmbpO6tPNl+heKHqGZFr/s1M48nPF8uqZ+ZFSQUZP0lzY2eb5bUIWH+Xq3IkJglcf39W5rZ3f8k6U/R4drfSvqp4ofsEvPKzNpJ+qukKyQ95e7bzexJxQsR7T5/ZKmkz7v7pH3Ivz+WSvo/d7+6hXkS862Q1NXMOiQUZM39mzUtO+m+4O4rFT80LDP7iKTnzeylqNAFsB/oGQNyR7GZlSY89vjPlruvkTRB0oOSFrn77Oj9FYofgrvDzMqiE9EPNrOPJlnXo5K+a2bl0flLX95LttcUL7huMLPi6GTzsyU9Ek2fLukCM+tgZodI+q992O5HJX3VzPqaWbmk7ySb0cyGmtkpUaG1VdIWxXuuJGmVpIFm1vS9WKL4Ic01khrM7OOKnzulhPm7mVmXhPfGSvofMxsQra/CzM7dh21prT9KOtvMzrD4BQylFr9Yo9krM6ND1FMl3WRmJRYfPuPsJPO2uC+Y2cUJ61mneNHHkBjAAaAYA3LHM4oXF02Pm5LM9ydJp2nnIcomVyhegMxS/I/sX7TrYa9EN0uqkbRI0vPRvNuSBXP3eknnKH4u0vuSfiPpCnefE81yl6R6xQuc30t6ONmymvE7Sc9Jekvxc+Ieb2HedpJuizKsVPzQ239H0x6Lftaa2RvuvlHxE/IfVbw9Pq2E86yi7H+WtDA6nFcp6RfRPOPNbKOkyZKObSFP4tWYTY8L97bB7r5U0rlR9jWK92Z9Wy1/p1+meE9oreIn/o9T8n+zlvaFUZJeM7NN0bZ+zd0X7S0zgOTMvbmedgBoPTP7oqRL3T1ZTxoyTDQkxRx3/2HoLEC+o2cMwD4zs95m9uHoENZQSd+S9EToXEjOzEZFhxsLLD6u2rmSngydCwAn8APYPyWKn/w+SNIHip/79ZugibA3vRQ/hNtN8UPMX3T3N8NGAiBxmBIAACAoDlMCAAAERDEGAAAQUFadM9a9e3cfOHBg6BgAAAB7NW3atPfdvWJv82VVMTZw4EBNnTo1dAwAAIC9MrOkt2dLxGFKAACAgCjGAAAAAqIYAwAACCirzhkDAADpt337dtXU1Gjr1q2ho2Sk0tJS9e3bV8XFxfv1eYoxAADQopqaGnXu3FkDBw6UmYWOk1HcXbW1taqpqdGgQYP2axkcpgQAAC3aunWrunXrRiHWDDNTt27dDqjXkGIMAADsFYVYcgfaNhRjAAAgr02YMEH/+c9/mp320EMPqaKiQiNGjNjxmDVrVpuun3PGAABAXpswYYI6deqk448/vtnpY8aM0a9+9auUrZ+eMQAA0OamLVmnX784X9OWrGuT5f3xj3/UMcccoxEjRugLX/iCGhsbtWTJEg0ZMkTvv/++YrGYTjjhBI0fP16SdN5552nkyJGqqqrSvffeu2M5zz77rI4++mgdeeSROvXUU7V48WKNHTtWd911l0aMGKGXX365TfLuC3rGEkxbsk6TF9Zq9OBuGjmgPHQcAAAy0pjfvrrHe58c3luXHzdQW+obdeE9kzRn5UbFXCow6bBenfW5Dw/SxdX9tHZzvb74x2m7fHbcF45rcX2zZ8/WuHHjNGnSJBUXF+u6667Tww8/rCuuuEI33nijrr32Wh177LEaNmyYTj/9dEnSAw88oK5du2rLli0aNWqULrzwQsViMV199dV66aWXNGjQIK1du1Zdu3bVtddeq06dOun6669vdv3jxo3TK6+8suP1q6++qvbt2+9rsyVFMRaZtmSdLrtvsuobYiopKtDDV42mIAMAYD9s2NqgmMefxzz++kC88MILmjZtmkaNGiVJ2rJli3r06CFJuuqqq/TYY49p7Nixmj59+o7P3H333XriiSckSUuXLtW8efO0Zs0anXjiiTuGoOjatWur1p/qw5QUY5HJC2tV3xBTzKX6hpgmL6ylGAMAoBkt9WS1LynULy49SpfdN1nbG2IqLirQLy49asff1K4dS/baE7Y7d9eVV16pW2+9dY9pdXV1qqmpkSRt2rRJnTt31oQJE/T888/r1VdfVYcOHXTSSSdp69atcveMvCqUc8Yiowd3U0lhvDliLlVVlgVOBABAdho5oFwPXzVa3zx9aJscaTr11FP1l7/8RatXr5YkrV27VkuWLJEk3Xjjjbrssst088036+qrr5YkrV+/XuXl5erQoYPmzJmjyZMnS5KOO+44TZw4UYsWLdqxHEnq3LmzNm7ceEAZD0SQYszMbjGzt81supmNN7PKEDkSjRxQroevHq0rjhugwgJp3JSlcvfQsQAAyEojB5TrSycf0iZHmYYNG6Yf//jHOv300zV8+HB97GMf04oVKzRx4kRNmTJlR0FWUlKiBx98UGeeeaYaGho0fPhwff/739fo0aMlSRUVFbr33nt1wQUX6Mgjj9SYMWMkSWeffbaeeOKJpCfwjxs3bpehLZINg7G/LETBYWZl7r4hev5VScPc/dq9fa66utqnTp2a8nxjJy7Qbf+co9svGq5LqvulfH0AAGSy2bNn6/DDDw8dI6M110ZmNs3dq/f22SA9Y02FWKSjpIzqgrr6hMEaPbirfvnvedreGAsdBwAA5LBgJ/Cb2f9IukLSekkntzDfNZKukaT+/funJVthgel/xxylggKpuJDT6gAAQOqkrNIws+fNbEYzj3Mlyd2/5+79JD0s6cvJluPu97p7tbtXV1RUpCruHnp1KVWPzqVqjLmmLF6btvUCAID8krJizN1Pc/cjmnk8tdusf5J0YapyHKjfvrRAY377apuNIAwAQDbiorbkDrRtQl1NOSTh5TmS5oTI0RqXjx6gPuXt9Y1x07Vp24ENWgcAQDYqLS1VbW0tBVkz3F21tbUqLS3d72WEOmfsNjMbKikmaYmkvV5JGUrn0mLddckIXfLbV3XT0zP184uPDB0JAIC06tu3r2pqarRmzZrQUTJSaWmp+vbtu9+fD1KMuXvGHpZsTvXArvryyYfo7n/P18lDe+gTw3uHjgQAQNoUFxfvuIUQ2h63Q2qlr5w6RItq69Sry/53QwIAAOyOYqyVigsL9MtPHRU6BgAAyDEMorWPGhpj+tHfZup3Ly0MHQUAAOQAirF9VFhgWv7BFt3+3BzNXL4+dBwAAJDlKMb2kZnptguGq7xDib72yHRt3d4YOhIAAMhiFGP7obxjie645EjNX71Jtz4zO3QcAACQxSjG9tMJQyr0+Q8P0qNTa7Ry/dbQcQAAQJbiasoDcMOZQ3XZ6P4MdwEAAPYbPWMHoLS4UAdXdJK7a9L897lNBAAA2GcUY21gwrtrdNl9r+lPr78XOgoAAMgyFGNt4KOHVuiEId11y99naf7qTaHjAACALEIx1gYKCkx3XHyk2hcX6uvj3lR9Qyx0JAAAkCUoxtpIj7JS3XbhcM1YtkF3PT83dBwAAJAlKMba0BlVvfSFEwfryL5dQkcBAABZgqEt2th3zzp8x3N3l5kFTAMAADIdPWMpcv8ri3T9Y28z3AUAAGgRxViK1G1r0F/fqNGT05eFjgIAADIYxViKXHfyIaoeUK4fPDlTS9fWhY4DAAAyFMVYihQWmO4aM0KS9M1Hp6sxxuFKAACwJ4qxFOrXtYNuPq9KU5es0+uL1oaOAwAAMhBXU6bYeSP66IjKLhrSs3PoKAAAIAPRM5ZiZrajEHttYa02b2sInAgAAGQSirE0Wbq2Tp++7zXd8vdZoaMAAIAMQjGWJv26dtA1Jw7WI1OW6tkZK0PHAQAAGYJiLI2+cdqhOqJPmb77+NtatWFr6DgAACADUIylUUlRgf53zFHasr1R1/xhqn7173matmRd6FgAACAgirE0O6RHJ33u+IGauXyD7vzXXF1232QKMgAA8hjFWACdSosVc1fMpe0NMU1eWBs6EgAACIRiLIDRg7uppKhApvjQF6MHdwsdCQAABEIxFsDIAeV6+KrRGlZZJplUeVBp6EgAACAQirFARg4o19jPjFSBTHe/MC90HAAAEAjFWED9unbQp4/tr0en1mjhmk2h4wAAgAAoxgL78imHqF1Rge7419zQUQAAQADcKDyw7p3a6XufOFw9O3PeGAAA+YhiLANcduyA0BEAAEAgHKbMEHX1Dfr5c+/q9UVrQ0cBAABpRDGWIQrM9JdpNbrtn7Pl7qHjAACANKEYyxClxYX62mlD9MZ7H+iF2atDxwEAAGlCMZZBLh7ZV4O6d9TPnntXjTF6xwAAyAcUYxmkqLBA3zr9UL27aqOefmtZ6DgAACANuJoyw5x1RG997sPrdFivstBRAABAGlCMZZiCAtMPz64KHQMAAKQJhykzVM26Ov2/J99RXX1D6CgAACCFKMYy1KoN2/THye/pwUmLQ0cBAAApRDGWoUYOKNdph/fQ2IkL9EFdfeg4AAAgRSjGMtj1ZwzVpm0NumfigtBRAABAilCMZbDDepXpvBF99NCkxVq5fmvoOAAAIAW4mjLDfeO0Q9WxXaGKCy10FAAAkAIUYxmuf7cO+vF5HwodAwAApAiHKbPEtCVr9esX54eOAQAA2hjFWJZ4YfZq/ey5dzVz+frQUQAAQBuiGMsSX/jowerSvlg/f+7d0FEAAEAbohjLEl3aF+vajx6sF99do9cXrQ0dBwAAtBGKsSzy2eMHqkfndrr92Tly99BxAABAG+BqyizSvqRQ3/n4YVq1YZsaY64ihrsAACDrUYxlmQuO7hs6AgAAaEMcpsxCsZjrqenL9MLsVaGjAACAA0TPWJYaO3Gh6uobdOKhFSoupKYGACBb8Vc8CxUUmL59xqFaUluncVOWho4DAAAOQNBizMyuNzM3s+4hc2Sjk4f20KiB5br7hXnaUt8YOg4AANhPwYoxM+sn6WOS3guVIZuZmW448zCt3rhND/1nceg4AABgP4XsGbtL0g2SGDBrP40a2FVXHDdAgys6ho4CAAD2U5AT+M3sHEnL3P0ts5bHyjKzayRdI0n9+/dPQ7rscvO5R4SOAAAADkDKesbM7Hkzm9HM41xJ35P0g9Ysx93vdfdqd6+uqKhIVdysVlffoN9MmK/VG7eGjgIAAPZRynrG3P205t43sw9JGiSpqVesr6Q3zOwYd1+Zqjy5bPWGbbpz/FytXL+VnjIAALJM2s8Zc/d33L2Huw9094GSaiQdTSG2/wZ276gxo/rpz6+/p/dq60LHAQAA+4BxxnLEV08dosIC013Pzw0dBQAA7IPgxVjUQ/Z+6BzZrmdZqa48fqCenL5Mc1ZuCB0HAAC0ErdDyiFf/OjBmr9qU+gYAABgH1CM5ZCDOpTo/s+OCh0DAADsg+CHKdH2Vq7fqntfWiB3xtMFACDTUYzloPGzVuonz8zRxLlrQkcBAAB7QTGWgy4d1V/9urbXz557V7EYvWMAAGQyirEcVFJUoG+cdqhmLt+gf7yzInQcAADQAoqxHHXuiD4a2rOz7vzXXG1vjIWOAwAAkqAYy1GFBaYbzhyqo/uXq66+MXQcAACQBENb5LBTD++pUw/vqWlL1mnywlqNHtxNIweUh44FAAASUIzluGlL1ulTv5ushsaYSooK9PBVoynIAADIIBymzHGTF9aqviGmmEvbG2KavLA2dCQAAJCAYizHjR7cTcUFJkkqLCzQ6MHdAicCAACJKMZy3MgB5br/s6NUWGA6eWgFhygBAMgwFGN54MRDK3T28N56dUGttm7nykoAADIJxVieuKS6n8xMc1dtDB0FAAAk4GrKPDF6cDe99t+nqrS4MHQUAACQgJ6xPFFQYCotLpS7c6gSAIAMQjGWR7Zub9Spd07U2IkLQkcBAAARirE8UlpcqD4HtddjU2sUi3noOAAAQBRjeeeikX217IMtepXBXwEAyAgUY3nmjKpeKist0qNTl4aOAgAARDGWd0qLC3XuiD56dsZKrd+yPXQcAADyHkNb5KErjx+oDx/STe0Z5gIAgOAoxvLQIT066ZAenULHAAAA4jBl3vqgrl6/eH6e5jEiPwAAQVGM5amYS796cZ4emcKJ/AAAhEQxlqe6dizRx4b11BNvLlN9Qyx0HAAA8hbFWB67uLqf1m6u17/nrAodBQCAvEUxlsdOHFKhnmXt9OjUmtBRAADIWxRjeaywwDRmVH+1KypQI7dHAgAgCIa2yHPfOG2IzCx0DAAA8hY9Y3muqRBburZO7vSOAQCQbhRj0ItzVuuE21/U1CXrQkcBACDvUIxBxwzqqo4lhXqUMccAAEg7ijGoY7sifXJ4pf7xzgpt3tYQOg4AAHmFYgySpEtG9VVdfaP+8c6K0FEAAMgrFGOQJB3dv1yDKzrq8TcYcwwAgHRiaAtIil9V+YsxR6lPefvQUQAAyCsUY9jhQ327hI4AAEDe4TAldvHKvPd1zR+mqqGRm4cDAJAOFGPYxeb6Bo2ftUovzVsTOgoAAHmBYgy7OOWwHureqUSPTuFEfgAA0oFiDLsoLizQ+Uf10fOzV6l207bQcQAAyHkUY9jDxdX91BBzPTl9eegoAADkPIox7OHQnp11SXVf9SorDR0FAICcx9AWaNbtFx0ZOgIAAHmBnjEktWHrdk1ZvDZ0DAAAchrFGJL60dOz9PmHpmjr9sbQUQAAyFkUY0jqwpF9tHFrg56buTJ0FAAAchbFGJIaPaib+nVtr0enLg0dBQCAnEUxhqQKCkwXj+ynSfNrtXRtXeg4AADkJIoxtOjCkX1lJk2cy+2RAABIhb0ObWFmQyTdKmmYpB0DT7n74BTmQoboc1B7vXzDyepb3iF0FAAAclJresYelHSPpAZJJ0v6g6T/S2UoZJamQszdAycBACD3tKYYa+/uL0gyd1/i7jdJOiW1sZBpvv/kDN3417dDxwAAIOe0phjbamYFkuaZ2ZfN7HxJPVKcCxnGTHpy+nKtr9seOgoAADmlNcXY1yV1kPRVSSMlXS7pylSGQua5pLqf6htievptbh4OAEBb2usJ/O4+JXq6SdLnUhsHmaqqskyH9y7TY1OX6vLRA0LHAQAgZyQtxszsb5KSnrHt7uekJBEykpnpkuq++tHfZmn2ig06vHdZ6EgAAOSElnrGfh79vEBSL0l/jF5/StLiA1mpmd0k6WpJTYNX/be7P3Mgy0TqnTeij1Zv3KbyDiWhowAAkDOSFmPuPlGSzOwWdz8xYdLfzOylNlj3Xe7+873PhkxR3rFEN555WOgYAADklNacwF9hZjsGeDWzQZIqUhcJmSwWc7347mpNW7IudBQAAHJCa4qxb0iaYGYTzGyCpBclfa0N1v1lM3vbzB4ws/JkM5nZNWY21cymrlnDLXlCc0nf/es7+tW/54WOAgBATthrMebuz0oaongB9jVJQ919/N4+Z2bPm9mMZh7nKj6i/8GSRkhaIemOFtZ/r7tXu3t1RQUdcqEVFpguHNlHE+eu0cr1W0PHAQAg67Xm3pSlkq6T9BHFO0ZeNrOx7t7iX2J3P601Aczsd5L+3pp5kRkuHtlPv35xgR5/s0bXnXRI6DgAAGS11hym/IOkKkm/lPQrxW8YfkD3pjSz3gkvz5c040CWh/Qa2L2jjhnUVY9NreF+lQAAHKC99owpfljyyITXL5rZWwe43tvNbITiPW2LJX3hAJeHNLukup/uHP+uVm7Yqt5d2oeOAwBA1mpNMfammY1298mSZGbHSpp0ICt198sP5PMI79wRlTr/qD4qLLDQUQAAyGotjcD/juI9V8WSrjCz96JJ/SXNSkM2ZLDiwvgR7u2NMcXc1a6oMHAiAACyU0s9Y59MWwpkpRXrt+jsX07St884VGNG9Q8dBwCArJT0BH53X9L0kLRBUhdJ3RIeyHO9ykrVpX2RHptaEzoKAABZqzVDW9wi6bOSFmjnjcNd0impi4VsEL95eD/d+s85WrBmkw6u6BQ6EgAAWac1Q1tcIulgdz/J3U+OHhRikCSdf3T8JH56xwAA2D+tKcZmSDoo1UGQnXp0LtXJQ3vo8Tdq1NAYCx0HAICs05qhLW5VfHiLGZK2Nb3p7uekLBWyyldOOUSbtzWowBjmAgCAfdWaYuz3kn4q6R1JdH1gD0f2o+MUAID91Zpi7H13vzvlSZDVln2wRfe/vEjXnXywundqFzoOAABZozXnjE0zs1vN7DgzO7rpkfJkyCqbtzXogUmL9OSby0JHAQAgq7SmZ+yo6OfohPcY2gK7OLRnZw3p0UljJy7QUf0O0siBXUNHAgAgK+y1GHP3k9MRBNlt2pJ1WvT+ZjXEXJ++7zX96erRGjmgPHQsAAAyXmt6xmRmn5BUJam06T13vzlVoZB9Ji+sVczjYwLXN8Y0eWEtxRgAAK2w13PGzGyspDGSviLJJF0saUCKcyHLjB7cTSVF8d2pqMA0ejB3zAIAoDVa0zN2vLsPN7O33f1HZnaHpMdTHQzZZeSAcj181WhNXlir0YO70SsGAEArtaYY2xL9rDOzSkm1kgalLhKy1cgB5RRhAADso9YMbfF3MztI0s8kvSFpsaRHUhkK2WvS/Pd10s9e1Hu1daGjAACQFVpzNeUt0dO/mtnfJZW6+/rUxkK26lxapMW1dZqxfL36d+sQOg4AABkvaTFmZhe0ME3uznlj2MOhPTursMA0c/l6nfWh3qHjAACQ8VrqGTu7hWkuTuJHM0qLCzWkRyfNXL4hdBQAALJC0mLM3T+XbJqZFZjZce7+ampiIZtVVXbRxLlrQscAACArtGrQV0kys16Szoweh0qaLIliDHs4aWiFigtN2xoa1a6oMHQcAAAyWkvnjBVK+rCkj0s6WdI6SeMl3eTuc9ITD9no7CMrdfaRlaFjAACQFVrqGXtd0iRJz0q62d23tDAvsItYzLWpvkFlpcWhowAAkNFaOmdsZDqDILd88pevqF/X9vrt5dWhowAAkNFaM+grsM8GVXTkikoAAFqBYgwpUVVZppp1W/RBXX3oKAAAZDSKMaREVWUXSdIsescAAGgRxRhSoqqyTJI4VAkAwF5QjCElundqpxvOHKpjBnUNHQUAgIzW6kFfgX113UmHhI4AAEDGo2cMKVNX36DXFtZq6/bG0FEAAMhYFGNImVfmva8x907WrBWcNwYAQDIUY0iZqj7xKypnLlsfOAkAAJmLYgwpU9mlVOUdirmiEgCAFlCMIWXMTFWVXTRjOT1jAAAkQzGGlKqqLNPclZu0vTEWOgoAABmJoS2QUpce01+nV/VSgVnoKAAAZCSKMaTUoO4dNah7x9AxAADIWBymRMqNn7lS/5q1KnQMAAAyEj1jSLnfvrRQBSZ9bFjP0FEAAMg49Iwh5Y6oLNOs5RsUi3noKAAAZByKMaRcVWUXba5v1OLazaGjAACQcSjGkHLDKsskicFfAQBoBsUYUu7Qnp1VXGh6d+XG0FEAAMg4nMCPlCub5dKfAAAWP0lEQVQpKtArN56iHp3bhY4CAEDGoRhDWvQsKw0dAQCAjMRhSqTF7BUb9K1H39KqDVtDRwEAIKNQjCEt6uob9Nc3avR2DTcNBwAgEcUY0uKwXmUyk2YupxgDACARxRjSomO7Ig3u3pHhLQAA2A3FGNKmqrKLZi6jZwwAgEQUY0ib4X27qEO7Im2pbwwdBQCAjMHQFkibq04YrKtOGBw6BgAAGYWeMQAAgIAoxpBW337sLX3/yRmhYwAAkDE4TIm0Wr9lu6YtWRc6BgAAGYOeMaTVEX26aOH7m7VpW0PoKAAAZASKMaRVVWWZpPjtkQAAQMBizMy+YmbvmtlMM7s9VA6k1xF9ukgS440BABAJcs6YmZ0s6VxJw919m5n1CJED6dejczudUdVT3Tq1Cx0FAICMEOoE/i9Kus3dt0mSu68OlANpZmb67eXVoWMAAJAxQh2mPFTSCWb2mplNNLNRyWY0s2vMbKqZTV2zZk0aIyKVtm5vVGPMQ8cAACC4lBVjZva8mc1o5nGu4j1y5ZJGS/q2pEfNzJpbjrvf6+7V7l5dUVGRqrhIownvrlbVD5/jJH4AAJTCw5TuflqyaWb2RUmPu7tLet3MYpK6S6LrKw8M7NZRjTHXjGXrd5zQDwBAvgp1mPJJSadIkpkdKqlE0vuBsiDN+nftoE7tijRzOT1jAACEOoH/AUkPmNkMSfWSrox6yZAHCgpMwyrLNGM5w1sAABCkGHP3ekmfCbFuZIaqyjI98vpSNcZchQXNni4IAEBe4N6UCOLjR/RWZZf22t4YU2FBYeg4AAAEQzGGII4Z1FXHDOoaOgYAAMFxb0oEs3L9Vs1fvTF0DAAAgqIYQzDX/N9U/eCpmaFjAAAQFMUYgqmq7KKZyzeIC2kBAPmMYgzBVFWWaf2W7apZtyV0FAAAgqEYQzBVlWWSxOCvAIC8RjGGYA7vXabCAtMsBn8FAOQxhrZAMKXFhRr7mZE6rFfn0FEAAAiGYgxBfWxYz9ARAAAIisOUCGr1hq368+vv6YO6+tBRAAAIgmIMQS18f7O++/g7enPpB6GjAAAQBMUYghoWXVE5iysqAQB5imIMQZWVFmtAtw6ayRWVAIA8RTGG4KoqyzRjGT1jAID8RDGG4Koqu+i9tXXauHV76CgAAKQdQ1sguMuO7a/Lju2vzqXFoaMAAJB2FGMI7qAOJaEjAAAQDIcpkREenLRID7yyKHQMAADSjmIMGeGluWs0bsrS0DEAAEg7ijFkhCP6dNH8NZu0dXtj6CgAAKQVxRgyQlVlmRpjrjkrN4aOAgBAWlGMISNUVXaRJAZ/BQDkHYoxZIS+5e3Vu0up1m9hrDEAQH5haAtkBDPTf75ziswsdBQAANKKnjFkDAoxAEA+ohhDxnhr6Qc679eTNG8VJ/EDAPIHxRgyRoeSQk1f+oHeruEkfgBA/qAYQ8YYXNFJpcUFmrl8Q+goAACkDcUYMkZhgenw3mUMbwEAyCsUY8goVZVlmrV8g2IxDx0FAIC0YGgLZJTjBnfXmo3btLm+QZ1Li0PHAQAg5SjGkFE+Mby3PjG8d+gYAACkDYcpkZHqG2KhIwAAkBYUY8g4Vz7wuq76w9TQMQAASAuKMWScHp3baeay9XLnJH4AQO6jGEPGqaosU+3meq3asC10FAAAUo5iDBmnqk8XSWK8MQBAXqAYQ8Y5vHeZzMRI/ACAvEAxhozTqV2RvnLKEB3dvzx0FAAAUo5xxpCRvvmxQ0NHAAAgLegZQ0ZqaIxp7qqN2lLfGDoKAAApRTGGjPTqwlqdftdLevO9daGjAACQUhRjyEhVlfErKmdwRSUAIMdRjCEjde1YosoupVxRCQDIeRRjyFjDKrtQjAEAch7FGDJWVWWZFqzZpLr6htBRAABIGYa2QMY6Z0SljuzXRYUFFjoKAAApQzGGjHVwRScdXNEpdAwAAFKKw5TIaFMWr9WEd1eHjgEAQMrQM4aM9ovn5+mDLfU6aWiP0FEAAEgJesaQ0ar6lGnuyk2qb4iFjgIAQEpQjCGjVVV2UX1jTPNWbwwdBQCAlKAYQ0arqiyTJMYbAwDkLIoxZLRB3TqqY0mhZlGMAQByFCfwI6MVFJie/spH1Le8fegoAACkBMUYMh5jjQEAchmHKZHxltRu1i1/n6Wla+tCRwEAoM1RjCHjbd7WqPtfWaQ33lsXOgoAAG0uSDFmZuPMbHr0WGxm00PkQHYY0rOTSgoLOIkfAJCTgpwz5u5jmp6b2R2S1ofIgexQXFigob06M7wFACAnBT1MaWYm6RJJfw6ZA5mvqrJMM5avl7uHjgIAQJsKfc7YCZJWufu8ZDOY2TVmNtXMpq5ZsyaN0ZBJqirLVFRgWle3PXQUAADalKWqp8HMnpfUq5lJ33P3p6J57pE0393vaM0yq6urferUqW2YEtmiMeYqMCnemQoAQOYzs2nuXr23+VJ2zpi7n9bSdDMrknSBpJGpyoDcUVhAEQYAyE0hD1OeJmmOu9cEzIAscts/5+iHT80IHQMAgDYVshi7VJy4j32wYv0WjZ+1KnQMAADaVLBizN0/6+5jQ60f2eeIyi5asX6rajdtCx0FAIA2E/pqSqDVqirLJInxxgAAOYViDFljGMUYACAHUYwhaxzUoUQnD61Qp9IgN44AACAl+KuGrPLg544JHQEAgDZFzxiyjrsrFuO2SACA3EAxhqwybck6Df/ReE17b13oKAAAtAmKMWSVvuXttXFrg2YsWx86CgAAbYJiDFmlR+d26t6phCsqAQA5g2IMWcXMVFXZhWIMAJAzKMaQdaoqyzRv1UZta2gMHQUAgAPG0BbIOicf1kNm0raGmNoVFYaOAwDAAaEYQ9YZNbCrRg3sGjoGAABtgsOUyEqT5r+vnzwzW9OWpH6Ii2lL1unXL87PqXXl4jalc125uE3pXFcublM615WL25TL62oNesaQdaYtWafL739NMZd+9/JCDenRSR3bFens4ZX6/EcGaUt9oz593+Q9PndJdT996pj+Wru5Xv/1+yl7TL/yuIE676g+Wv7BFn3pT29IkjZva9C81ZvkLhUXmh655jgd1KFY1z/21h6f//pph+qjh1ZoxrL1+v5TM/aY/t2PH65jBnXVlMVr9ZNnZu8x/dOj+uv7T89QfUNMLu3YriY/u2i4DunRWeNnrtQ9Exfs8flffuoo9S3voKemL9ND/1m8x/T7rqhWt07tdPuzs3XPxIVyl8x2rufhq45Vh5IiPTRpkZ56a/ken3/8i8fLzHTPhAUaP2vlLtNKiwr152tGS5Lu/NdcvTxvzR7tV1pcoIevGq3nZq7UlMVrd/l87y6l+s1lIyVJ339yhmYs33XokkHdO+rOS0ZIkq5/7C0tWLNpl+nDepfpgqP76rL7Jmvr9tgu2yVJ1QPK9b1PDJMk/ddDU7S2rn6Xz59wSHd98/ShkqRP/26ytmzf9XzE0w7vqS+dfIgk6fS7Ju7Ypqb1XDqqf5vve4ntJ5faFRfo9ouG68FJi/f4/IHuezefc4TqG2P61L2TVd+4Z/u11b43bsp7uv+VRXu035Nf+nCb73tN7Td/dXxfKSkq0BlVvfTe2rpdPt8W+97/nP8hTVuyTheP/Y9iu/1eteW+d/5vJu3yO2Umffb4gfrh2VVtvu81td+i9zerMeYqLizQgG4ddvlOktpm3/tQ3y66/5WF+vE/Zu/xvdSW+94jU5bu2K6FazYr5q6Sovj30sgB5Xt8Np0oxpB1Ji+slUcD8LtL9Q0x9SwrUrvieEevmdSp3Z67dklhNF1JphfFpxeY7Zi+bnP9jnU1xlyTF9bqrA/1bvbzxQUmSSossGanF0bTi5JMn7F8g+obYmq6uUDTdjUpsPjni4sKmv180/SSwuanWzR9SW1ds+23sx0Km/18k3bNrD/x3L3S4oJm2297Q0yTF9aqtHjP5Xco2fm6fcme0zsmTO/QzPQOJYWavLBW9Q2xZrertHhnvo7tilTfGNs1f8L0Tu2KdvxbJW5zk/qG2B7tl4p9T9qz/d6uWZ+Sfa+wwDR5Xq0aYs23X1vteyVFBc223852aLt9T4q3Xyyh/VZv3JaSfU+Kfy/Fmvm9ast9r1O7ol32CXdp2botktp+35Pi7dfQ6HJJ2xtje3xXSG2z70nS3FUbm/1east9L/F7qTEWbVf0vRS6GDNv2vosUF1d7VOnTg0dA4FNW7JOl903WdsbYipO8f9qcnFdubhN6VxXLm5TOteVi9uUznXl4jbl8rrMbJq7V+91PooxZKNpS9Zp8sJajR7cLeX/o8nFdeXiNqVzXbm4TelcVy5uUzrXlYvblKvrohgDAAAIqLXFGFdTAgAABEQxBgAAEBDFGAAAQEAUYwAAAAFRjAEAAAREMQYAABAQxRgAAEBAFGMAAAABUYwBAAAERDEGAAAQEMUYAABAQFl1b0ozWyNpSYpX013S+yleR7agLXaiLXaiLeJoh51oi51oi51oC2mAu1fsbaasKsbSwcymtuamnvmAttiJttiJtoijHXaiLXaiLXaiLVqPw5QAAAABUYwBAAAERDG2p3tDB8ggtMVOtMVOtEUc7bATbbETbbETbdFKnDMGAAAQED1jAAAAAeVtMWZmZ5rZu2Y238y+08z0dmY2Lpr+mpkNTH/K1DOzfmb2opnNNrOZZva1ZuY5yczWm9n06PGDEFnTwcwWm9k70XZObWa6mdnd0X7xtpkdHSJnKpnZ0IR/6+lmtsHMvr7bPDm7T5jZA2a22sxmJLzX1cz+ZWbzop/lST57ZTTPPDO7Mn2pUyNJW/zMzOZE+/8TZnZQks+2+LuUbZK0xU1mtizh9+CsJJ9t8e9NNknSDuMS2mCxmU1P8tmc2ifalLvn3UNSoaQFkgZLKpH0lqRhu81znaSx0fNLJY0LnTtFbdFb0tHR886S5jbTFidJ+nvorGlqj8WSurcw/SxJ/5RkkkZLei105hS3R6GklYqPlZMX+4SkEyUdLWlGwnu3S/pO9Pw7kn7azOe6SloY/SyPnpeH3p4UtMXpkoqi5z9tri2iaS3+LmXbI0lb3CTp+r18bq9/b7Lp0Vw77Db9Dkk/yId9oi0f+dozdoyk+e6+0N3rJT0i6dzd5jlX0u+j53+RdKqZWRozpoW7r3D3N6LnGyXNltQnbKqMdq6kP3jcZEkHmVnv0KFS6FRJC9w91YMtZwx3f0nS2t3eTvw++L2k85r56BmS/uXua919naR/STozZUHToLm2cPfx7t4QvZwsqW/agwWQZL9ojdb8vckaLbVD9DfyEkl/TmuoHJCvxVgfSUsTXtdozwJkxzzRF896Sd3Ski6Q6FDsUZJea2bycWb2lpn908yq0hosvVzSeDObZmbXNDO9NftOLrlUyb9Y82WfkKSe7r5Civ8HRlKPZubJt31Dkj6veE9xc/b2u5Qrvhwdsn0gyeHrfNovTpC0yt3nJZmeL/vEPsvXYqy5Hq7dLyttzTw5w8w6SfqrpK+7+4bdJr+h+GGqIyX9UtKT6c6XRh9296MlfVzSl8zsxN2m581+YWYlks6R9Fgzk/Npn2itvNk3JMnMviepQdLDSWbZ2+9SLrhH0sGSRkhaofghut3l037xKbXcK5YP+8R+yddirEZSv4TXfSUtTzaPmRVJ6qL966LOeGZWrHgh9rC7P777dHff4O6boufPSCo2s+5pjpkW7r48+rla0hOKH2JI1Jp9J1d8XNIb7r5q9wn5tE9EVjUdjo5+rm5mnrzZN6KLEz4p6TKPTgbaXSt+l7Keu69y90Z3j0n6nZrfxrzYL6K/kxdIGpdsnnzYJ/ZXvhZjUyQNMbNB0f/+L5X09G7zPC2p6WqoiyT9O9mXTjaLjvHfL2m2u9+ZZJ5eTefLmdkxiu83telLmR5m1tHMOjc9V/xE5Rm7zfa0pCuiqypHS1rfdPgqByX9X26+7BMJEr8PrpT0VDPzPCfpdDMrjw5XnR69l1PM7ExJN0o6x93rkszTmt+lrLfb+aLnq/ltbM3fm1xwmqQ57l7T3MR82Sf2W+grCEI9FL8qbq7iV7l8L3rvZsW/YCSpVPHDM/MlvS5pcOjMKWqHjyjeZf62pOnR4yxJ10q6Nprny5JmKn4V0GRJx4fOnaK2GBxt41vR9jbtF4ltYZJ+He0370iqDp07RW3RQfHiqkvCe3mxTyhegK6QtF3xXo3/Uvx80RckzYt+do3mrZZ0X8JnPx99Z8yX9LnQ25Kitpiv+DlQTd8XTVedV0p6Jnre7O9SNj+StMX/Rd8DbyteYPXevS2i13v8vcnWR3PtEL3/UNP3Q8K8Ob1PtOWDEfgBAAACytfDlAAAABmBYgwAACAgijEAAICAKMYAAAACohgDAAAIiGIMQEqZ2a1mdpKZnWdm30kyz7VmdsU+LneCmVVHz58xs4NamLfSzP6yt+XsKzP7rJn9qpXz7shgZiPM7Kz9WSeA3EMxBiDVjlX8fqcflfRyczO4+1h3/8P+rsDdz3L3D1qYvtzdL9rf5beF3TKMUHzsKQCgGAOQGmb2MzN7W9IoSa9KukrSPWb2g2bmvcnMro+eTzCzn5rZ62Y218xOiN5vb2aPRDdlHiepfcLnF5tZ9+hz1+223G+Z2UAzm9GK5WxKeH6RmT0UPT/bzF4zszfN7Hkz67mXbf+omU2PHm+aWeemDNEo7DdLGhNNHxONTv6AmU2J5j83Wk5V1A7To7xD9vGfAUAWKAodAEBucvdvm9ljki6X9E1JE9z9w638eJG7HxMdyvuh4rda+aKkOncfbmbDFb9Z+e4ekfS/kn4Tvb5E0pna9T+erVnO7l6RNNrd3cyuknSDpG+1MP/1kr7k7pPMrJOkrU0T3L0+Kkir3f3LkmRmP1H8lmufjw63vm5mzyt+14NfuPvDURFX2IqsALIMxRiAVDpK8VvmHCZp1j58rumG9dMkDYyenyjpbkly97ejXrdduPubZtbDzColVUha5+7vmdnAhNn2upxm9JU0LroXYYmkRXuZf5KkO83sYUmPu3tNdCvPZE6XdE5T76Dit2Prr3iP4vfMrG+0nHmtyAogy1CMAWhzZjZC8XvV9ZX0vuL3ujQzmy7pOHffspdFbIt+NmrX76nW3L/tL5IuktRL8Z6y5iRbTuL7pQnPfynpTnd/2sxOknRTSwHc/TYz+4fi54VNNrPTlNA71gyTdKG7v7vb+7PN7DVJn5D0nJld5e7/bmndALIP54wBaHPuPt3dRyh+c+Rhkv4t6Qx3H9GKQiyZlyRdJklmdoSk4Unme0TSpYoXZM1dQdnSclaZ2eFmViDp/IT3u0haFj2/cm9Bzexgd3/H3X8qaariPYOJNkrqnPD6OUlfsaj7zMyOin4OlrTQ3e9W/EbUybYZQBajGAOQEmbWdJgwJukwd9+Xw5TNuUdSp+iw4g2SXm9uJnefqXihs8zdV+zjcr4j6e+KF4+Jn71J0mNm9rLiPX178/XoZP23JG2R9M/dpr8oaVjTCfySbpFULOnt6EKDW6L5xkiaEfUoHiZpv684BZC5zL01vf4AAABIBXrGAAAAAqIYAwAACIhiDAAAICCKMQAAgIAoxgAAAAKiGAMAAAiIYgwAACAgijEAAICA/j8j6BeyYd883AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axes = plt.subplots(nrows=1, ncols=1,figsize=(10,6))\n",
    "\n",
    "ax=axes\n",
    "ax.plot(step_es,'.--',label='exact E')\n",
    "ax.set(xlabel='# individual sites', ylabel='\\\\lambda',\n",
    "       title='Live ground state Energies')\n",
    "ax.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exact Diagonalization\n",
    "This part is copy-paste from TeNPy toycodes (https://github.com/tenpy/tenpy) and provides the exact diagonalization of the same transverse Ising Hamiltonian"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "\"\"\"Provides exact ground state energies for the transverse field ising model for comparison.\n",
    "\n",
    "The Hamiltonian reads\n",
    ".. math ::\n",
    "    H = - J \\\\sum_{i} \\\\sigma^x_i \\\\sigma^x_{i+1} - g \\\\sum_{i} \\\\sigma^z_i\n",
    "\"\"\"\n",
    "import numpy as np\n",
    "import scipy.sparse as sparse\n",
    "import scipy.sparse.linalg.eigen.arpack as arp\n",
    "import warnings\n",
    "import scipy.integrate\n",
    "\n",
    "\n",
    "def exact_E(L, J, g):\n",
    "    \"\"\"For comparison: obtain ground state energy from exact diagonalization.\n",
    "\n",
    "    Exponentially expensive in L, only works for small enough `L` <~ 20.\n",
    "    \"\"\"\n",
    "    if L >= 20:\n",
    "        warnings.warn(\"Large L: Exact diagonalization might take a long time!\")\n",
    "    # get single site operaors\n",
    "    sx = sparse.csr_matrix(np.array([[0., 1.], [1., 0.]]))\n",
    "    sz = sparse.csr_matrix(np.array([[1., 0.], [0., -1.]]))\n",
    "    id = sparse.csr_matrix(np.eye(2))\n",
    "    sx_list = []  # sx_list[i] = kron([id, id, ..., id, sx, id, .... id])\n",
    "    sz_list = []\n",
    "    for i_site in range(L):\n",
    "        x_ops = [id] * L\n",
    "        z_ops = [id] * L\n",
    "        x_ops[i_site] = sx\n",
    "        z_ops[i_site] = sz\n",
    "        X = x_ops[0]\n",
    "        Z = z_ops[0]\n",
    "        for j in range(1, L):\n",
    "            X = sparse.kron(X, x_ops[j], 'csr')\n",
    "            Z = sparse.kron(Z, z_ops[j], 'csr')\n",
    "        sx_list.append(X)\n",
    "        sz_list.append(Z)\n",
    "    H_xx = sparse.csr_matrix((2**L, 2**L))\n",
    "    H_z = sparse.csr_matrix((2**L, 2**L))\n",
    "    for i in range(L - 1):\n",
    "        H_xx = H_xx + sx_list[i] * sx_list[(i + 1) % L]\n",
    "    for i in range(L):\n",
    "        H_z = H_z + sz_list[i]\n",
    "    H = -J * H_xx - g * H_z\n",
    "    E, V = arp.eigsh(H, k=1, which='SA', return_eigenvectors=True, ncv=20)\n",
    "    return E[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plotting\n",
    "Now we compare the results and see that they match perfectly  \n",
    "We write a brute force function that calculates the GS energy via DMRG within some confidence (to know when it is *close enough* to the GS energy. Here in this example one full sweep is already enough usually (put test=True to see how many sweeps it needs each time)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "def dmrg_E(L,J,g,conf=1e-4,test=False):\n",
    "    # Initialize random MPS and MPO\n",
    "    init_mps = random_MPS(L=L,d=2,chi_max=100)\n",
    "    mpo = Ising_MPO(L=L,d=2,h=g,J=J)\n",
    "    dmrg = DMRG(init_mps,mpo,eps=-1.,chi_max=100,test=False)\n",
    "    \n",
    "    diff = 5.*conf\n",
    "    counter = 0\n",
    "    while diff > conf:\n",
    "        if test:\n",
    "            print(counter)\n",
    "        counter += 1\n",
    "        ### note: we could also simply use the current dmrg.e value to speed things up. This is just to give justification to the expect_mpo(mps,mpo) function\n",
    "        e_in = expect_mpo(dmrg.MPS,mpo)\n",
    "        dmrg.left_to_right()\n",
    "        dmrg.right_to_left()\n",
    "        e_end = expect_mpo(dmrg.MPS,mpo)\n",
    "        diff = abs(e_in - e_end)\n",
    "        if counter > 50:\n",
    "            break\n",
    "    return dmrg.e[-1]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "gs = np.linspace(0.,2.,20)\n",
    "dmrg_Es = np.array([dmrg_E(L=10,g=g,J=1.) for g in gs])\n",
    "exact_Es = np.array([exact_E(L=10,g=g,J=1.) for g in gs])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f92e7bd9d68>"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmoAAAGDCAYAAACbcTyoAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJzs3Xd8VfX9x/HX595sElYg7A0ie0UE3HXWhYK4ENA6wNFabYs/Wm211bbWiq0DEaXaOnEGF24BGY6gjDBkg2EmzAAJubn3+/vjXmLAAAGSnJvk/Xw87sNz7z3n3Pe5avLJ93yHOecQERERkejj8zqAiIiIiJROhZqIiIhIlFKhJiIiIhKlVKiJiIiIRCkVaiIiIiJRSoWaiIiISJRSoSYiVZKZrTazs7zOUZWY2UIzO93rHCJSdirURKRUZnalmX1lZrvNbHNk+xYzM6+zHS0zO93Mso/wGGdm7Y/hM13kO9xV4jH6aM93LJxzXZxzU734bBE5OirUROQnzOw3wL+Bh4DGQCNgFHASEHeQY/yVFrDq6eGcSy7x+Ed5f4CZxZT3OUXEeyrURGQ/ZlYH+DNwi3Pudedcngv7zjk31Dm3N7Lfc2b2pJm9b2a7gTPMrI6Z/c/McsxsjZndbWa+yP73mtkLJT6ndaS1KSbyfKqZ/cXMZppZnpl9ZGYNSuw/LHLOLWb2h8Ncw/lmtihynnVm9lszqwVMAZqWaNlqamZ9zWy2mW03sw1m9riZxUXOMz1yynmR/a+IvH6hmc2NHDPLzLof5Xd9r5m9GvnO8iK3JtNLvN/UzN6IfJ+rzOxXBxz7upm9YGY7gWvNLNHM/mtm28xssZmNLtmCWPJ2sZn5zOz/zGxF5Dt91czqR95LiJx3S+QavzGzRkdzjSJybFSoiciB+gPxwOQy7Hs18ACQAswAHgPqAG2B04DhwHVH8NlXR/ZPI9xy91sAM+sMPAkMA5oCqUDzQ5xnIjDSOZcCdAU+c87tBn4OrC/RsrUeCAJ3AA0IX/uZwC0AzrlTI+fb1yI2ycx6A/8BRkZyPAW8bWbxR3CdJV0MvALUBd4GHo9csw94B5gHNIvk+rWZnVvi2IHA65FjXwT+BLQm/P2fDVxziM/9FXAJ4X9PTYFtwBOR90YQ/vfYInKNo4D8o7w+ETkGKtRE5EANgFznXNG+FyKtRtvNLN/MTi2x72Tn3EznXAgIAFcAYyKtcKuBhwkXV2X1rHNuqXMuH3gV6Bl5/TLgXefc9EiL3j1A6BDnCQCdzay2c26bc+7bg+3onJvjnPvSOVcUyfwU4eLlYG4EnnLOfeWcCzrn/gvsBfod4phvI9/fvkfJYmuGc+5951wQeB7oEXn9BKChc+7PzrlC59xK4GngyhLHznbOZTjnQpHv7HLgr5FrzgYePUSmkcAfnHPZke/0XuCySAtngHCB1j5yjXOcczsPcS4RqSDq0yAiB9oCNDCzmH3FmnNuAEDkNlrJP/B+KLHdgHAr2JoSr60h3BpUVhtLbO8BkiPbTUt+lnNut5ltOcR5BgN3A383s/nA/znnZpe2o5kdB4wF0oEkwj8X5xzi3K2AEWb2yxKvxUUyHkxv59zyg7x34DUnRIqlVoRv024v8b4f+KLE85LfPxzwPZXyfkmtgLfMrGTBGyTcH/F5wq1pr5hZXeAFwkVd4BDnE5EKoBY1ETnQbMItRAPLsK8rsZ1LuCWmVYnXWgLrItu7CRdC+zQ+gkwbCBcOAJhZEuEWn9JDOfeNc24g4VuoGYRb5w7Mu8+TwBKgg3OuNvB74FAjW38AHnDO1S3xSHLOvXwE11MWPwCrDvicFOfc+SX2OfB6NrD/LeEWHNwPwM8POH+Cc26dcy7gnLvPOdcZGABcSPg2tohUMhVqIrIf59x24D5gnJldZmbJkY7nPYFahzguSLggesDMUsysFXAn4dYYgLnAqWbWMjJgYcwRxHoduNDMTo509P8zB/n5ZWZxZjbUzOpEWoB2Em4pAtgEpEY+f5+UyD67zOx44OYDTrmJcJ+vfZ4GRpnZiRZWy8wuMLOUI7iesvga2Glmd0UGCfjNrKuZnXCIY14FxphZPTNrBtx2iH3HE/531QrAzBqa2cDI9hlm1s3CI3l3Ei7Agwc/lYhUFBVqIvITkekj7gRGA5sJFytPAXcBsw5x6C8Jt5ytJDy44CXCHe9xzn0MTALmE761+O4R5FkI3Bo53wbCHd8PNR/aMGB1ZDTkKCKd6p1zS4CXgZWRvmJNCQ9YuBrII1yETTrgXPcC/43sf7lzLpNwP7XHIzmWA9ce5hL2jRrd9/hXGa45CFxEuJ/eKsItls8Q7uR/MH8m/L2sAj4hXODuPci+/yY8eOEjM8sDvgROjLzXOHLsTmAxMI0fC24RqUTmXGl3AkREpKozs5uBK51zhxocISJRTC1qIiLVhJk1MbOTIreqOwK/Ad7yOpeIHD2N+hQRqT7iCN+ibgNsJzw/2zhPE4nIMdGtTxEREZEopVufIiIiIlFKhZqIiIhIlKoWfdQaNGjgWrdu7XUMERERkcOaM2dOrnOuYVn2rRaFWuvWrcnMzPQ6hoiIiMhhmdmaw+8VplufIiIiIlFKhZqIiIhIlFKhJiIiIhKlqkUfNREREakcgUCA7OxsCgoKvI4S9RISEmjevDmxsbFHfQ4VaiIiIlJm2dnZpKSk0Lp1a8zM6zhRyznHli1byM7Opk2bNkd9Ht36FBERkTIrKCggNTVVRdphmBmpqanH3PKoQk1ERESOiIq0simP70mFmoiIiEjE1KlTmTVrVqnvPffcczRs2JCePXsWPxYtWlShedRHTURERCRi6tSpJCcnM2DAgFLfv+KKK3j88ccrLY9a1Mpi1XQY1w/yNu2/LSIiIoc1Z802nvh8OXPWbCuX873wwgv07duXnj17MnLkSILBIGvWrKFDhw7k5uYSCoU45ZRT+OijjwC45JJL6NOnD126dGHChAnF5/nggw/o3bs3PXr04Mwzz2T16tWMHz+eRx55hJ49e/LFF1+US95joRa1w1k1HV66nFBRIesnDqVJXhb+UACmPQgXji3/z5oyGoZNhtzvf9xOaVS+nyMiIlJOrnhq9k9eu7B7E4b1b01+YZDBT85kycY8Qg58Bsc3TuG6k9owJL0FW3cXcvMLc/Y7dtLI/of8vMWLFzNp0iRmzpxJbGwst9xyCy+++CLDhw/nrrvuYtSoUZx44ol07tyZc845B4D//Oc/1K9fn/z8fE444QQGDx5MKBTixhtvZPr06bRp04atW7dSv359Ro0aRXJyMr/97W9L/fxJkyYxY8aM4uezZ88mMTHxSL+2MlOhdjhTRhMqKsTngqRum4ffCgEoynqL7H5/oXGdBBJi/RQFQzggxmdH13mwMgtCERGRSrKzoIiQC2+HXPj5sfj000+ZM2cOJ5xwAgD5+fmkpaUBcMMNN/Daa68xfvx45s6dW3zMo48+yltvvQXADz/8wLJly8jJyeHUU08tnjqjfv36Zfr8yr71qULtcIZlsH7iNaRum0dipEjLd3HcuXMYU/45lVdH9qdvm/q8PW89d746D59BXIyP+Bg/8TE+nr/+RDo2TuH9BRt4+ouVxJd4Ly7Gx58u6kLDlHh2T/4NCYFC/OxfELpFGdiFY3HOaZSNiIhEnUO1gCXG+fn3lb0Y+syXBIpCxMb4+PeVvejTqh4A9WvFHbYF7UDOOUaMGMHf/va3n7y3Z88esrOzAdi1axcpKSlMnTqVTz75hNmzZ5OUlMTpp59OQUFBlfm9qkLtcHKXhlu3IoUTQKzP8X8dN3FO9x60bVgLgM5Na/Pbc45jb1Eo/AgE2VsUonZi+CuO9ftIjo9hbyDE9vwAewNBCotChFz4z4w3Oj9G+xl30suW7VcQhs5+iFrAPz/6npe//oEGyXE0TImnQXL4MebnxxPj97EyZxf5gSANU+JJrRWP33eY//h0m1VERCpBn1b1ePGGfny5cgv92qYWF2lH68wzz2TgwIHccccdpKWlsXXrVvLy8mjVqhV33XUXQ4cOpVWrVtx44428++677Nixg3r16pGUlMSSJUv48ssvAejfvz+33norq1at2u/WZ0pKCjt37iyPSy8XKtQOZ8ro8C1IIOCLx0+ImFCAVps+plWvp4p3O75xbY5vXPugpzm7cyPO7nzwImh4h0KCX63EH/yxIIzzOXzrZkCvwfRoXpftewLk7tpLTt5evlu7nV17i7j7gk4APP7Zct78bh0AZpBaK47WqbV4/ebwqJWM79axOa+AhinxtNv1LV2n3QRB3WYVEZGK16dVvWMu0Pbp3Lkz999/P+eccw6hUIjY2FieeOIJVq9ezTfffMPMmTPx+/288cYbPPvss1x99dWMHz+e7t2707FjR/r16wdAw4YNmTBhAoMGDSIUCpGWlsbHH3/MRRddxGWXXcbkyZN57LHHOOWUU/b7/AP7qI0bN+6gI0TLg7lIi05Vlp6e7jIzMyvm5HmbwgXMogy4YGy4JWpRBgx5DtqcWn6fM64f5CwDV1RcEPpCAUhKhdErD3v49xvzWJmzK1zI7SokJ28vPoMHLu0GwLCJX/HFslwAPowbTTtbT4yFyHdxxS14e2Lq8uEFM2nbIJm2DWuRknD0a5OJiEj1tHjxYjp16uR1jCqjtO/LzOY459LLcrwKtWhRwQWhc468vUXk5u1lx+Zsar07ihZ7FhYXaQXE8ZvALbwX7Ft8zPndGjNuaB8AXsv8Idwa1zCZpnUTD39rVUREqiUVakfmWAs13fqMFimNwrcd99167HJJud6GNDNqJ8RSOyEWdm0mWLhkv353cT7Ho/12cseJp7J8825W5u6ice0EAIqCIX7/1gICwXBRHxfjo01qLa4+sSUjBrTGOcf87B0Hb4VTfzgREZGjokKtJiql350vFIDFk2l/0SO0T0vZb3e/z/hyzJmsyNnNypxdrMzdzYrNu4iLCc+XvDlvLwOfmAlAWko8bRvWol3DZAb1bkafUBbupctxmnZERETkiKlQq4mGTS6+zRp74G3WUpgZqcnxpCbH07fNT+eZSUmI4alhfVgZKeRW5Ozi3fkbSG9djz6zR+NKmYcumJWB74KHq8TQaBEREa+oUKuJyvk2a1JcDOd2abzfa845nAPaZ7BmwtU03rnggHnormHYyi0MaNeA1bm7Wbc9n+7N62gAg4iISAkq1KRCmBlmQO5SWu5Z9JP+cHe020DT5nUBePO7dTz66TLMoENaMr1a1KNXy7pc2rsZ8TF+j65ARETEe1qUXSrWAf3hQr5Y/C7AcVs+JTk+/HfC9Se34X+/6MuvzzyOZnUT+XDRRv7y7iJifOH/PJ+duYp/fLCEjxdtIidvr2eXIiIi0cHv99OzZ0+6dOlCjx49GDt2LKFQCICpU6diZkycOLF4/++++w4z45///CcA1157LW3atKFnz5706NGDTz/9tHjfoqIifv/739OhQwd69uxJz549eeCBB0rN0bp1a7p161a8369+9atyv1a1qEnFKkN/uDqJsZx6XENOPa4hEL5tumnn3uIpQOZn7+CdeespiiwW16J+Imd3aswfL+oMQCjk8PlMo0tFRKJRBfxsTkxMLF7Lc/PmzVx99dXs2LGD++67D4Bu3boxadIkrr/+egBeeeUVevTosd85HnroIS677DI+//xzbrrpJpYtWwbA3XffzcaNG1mwYAEJCQnk5eXx8MMPHzTL559/ToMGDY76Wg5HhZpUrKPoD2dmNK6TUPz8kSt68rdB3chat4Pv1m7nux+2EYz85QRwxsNTOSVmMX/Mu4+YUIC1E66ixZ5FGl0qIuK1VdPhpcuhKABv3ADrvglvl+PP5rS0NCZMmMAJJ5zAvffeC0DLli3ZuXMnmzZtIi0tjQ8++IDzzz+/1OP79+/PunXhlX327NnD008/zerVq0lICP8eSklJKT6vF1SoSZWQEOsnvXV90lvvP+o0EAxxXtfGDJ1zO75gIT4L0Wjngv0Wtd977kMkxKqvm4hIpZsyOlyYuSLI/gaK8sOvL8oo1z+i27ZtSygUYvPmzcWvXXbZZbz22mv06tWL3r17Ex8fX+qxH3zwAZdccgkAy5cvp2XLlqSkpJS6b2nOOOMM/P7w75gRI0Zwxx13HMOV/JQKNanSYv0+xvy8E5z8ET9MHEqDbfOLR5cGfPFsHHA/Z933Ef3apnLacQ05vWND2jSopWlBREQqw7AMeOPG/Yu0mMTwCjzl7MCVli6//HKuuOIKlixZwlVXXcWsWbP2e/93v/sdo0ePZvPmzcULtR/o2Wef5d///jdbtmxh1qxZtGjR4if7VPStTw0mkOohdylN835cEgvAT4g6G2cz9MRW/LBtD39+dxE/e3gapz00le835nkYVkSkhshdGrndmf/jay4YviVajlauXInf7yctLa34tcaNGxMbG8vHH3/MmWee+ZNjHnroIZYvX87999/PiBEjAGjfvj1r164lLy/8O+K6665j7ty51KlTh2AwWK6Zy0qFmlQPpYwu9YUC1F75Hn+8qDOf/eZ0pv/uDP4ysAsdG6fQvF4iAOOnrWDYxK945ouVLN+86yd/kYmIyDHYd+sTwi1p/jgIFoZvfZaTnJwcRo0axW233faTuyV//vOfefDBB4tvTR7I5/Nx++23EwqF+PDDD0lKSuL666/ntttuo6CgAIBgMEhhYWGpx1cG3fqU6qEMo0tbpiYxrH9rhvVvXfxaYqyf9dvzuf+9xdz/3mKa10vkvC6NufvCzpV/DSIi1U2Jn82UYSWcssrPz6dnz54EAgFiYmIYNmwYd95550/2GzBgwGHPZWbcfffd/OMf/+Dcc8/lgQce4J577qFr166kpKSQmJjIiBEjaNq0aanHl+yj1r17d/73v/8d07X9JF91aEFIT093mZmZXseQKuyHrXuYujSHad9vJi7Gx7ihfQD40+QsWtRP4vSODWnXMPnHv9Y0FYiI1FCLFy+mU6dOXseoMkr7vsxsjnMuvSzHq0VNBGhRP4lh/VoxrF+r4tf2FgX5cuVW/jt7Dfe/t5hmdRM5vWNDrm26lg6fXE9IC82LiEgFU6EmchDxMX4+vONUsrftYdrSHKZ+n0PGd+v47Yo/EiplofnyHm4uIiLiyWACMxtiZgvNLGRm6Qe8N8bMlpvZ92Z2rhf5REpqXi+JoSe24unh6Xz3x3NI+MXbrK/Tm3wXVzzKtIA4vuz0B3btLfI4rYiIVCdejfrMAgYB+43PNbPOwJVAF+A8YJyZaaZSiRpxMT4Sd6ygSV7WflOB+Aix7Kv3OeH+T9i627vRQSIilaE69G+vDOXxPXly69M5txgobdLRgcArzrm9wCozWw70BWZXbkKRQzhgKhA/IeJCAa6s9S27B3Sgfq04AMa8OZ9Yv49LejWjV4u6mmRXRKqFhIQEtmzZQmpqqn6uHYJzji1bthQvRXW0oq2PWjOg5PTA2ZHXRKLHQaYCiR3yHKPatAPC/4PuLQrx5rfr+N/sNbRKTWJgz2YM6tWM1g1qeXwBIiJHr3nz5mRnZ5OTk+N1lKiXkJBA8+bNj+kcFTY9h5l9AjQu5a0/OOcmR/aZCvzWOZcZef4EMNs590Lk+UTgfefcG6Wc/ybgJoCWLVv2WbNmTYVch8ixyCsI8EHWRjLmrmPWii386mcduOPs4ygsCrGzIECD5NLXnhMRkeorKqbncM6ddRSHZQMlF9JqDqw/yPknABMgPI/aUXyWSIVLSYhlSHoLhqS3YOOOAmL94dsEn3+/mVte/JZTOzTgkl7NOLtzI5Lioq2BW0REvBZtS0i9DVxpZvFm1gboAHztcSaRctG4TgKpkRa0To1rM/LUtizdtIvbX5lL+v2fcMekueQVBDxOKSIi0cSr6TkuNbNsoD/wnpl9COCcWwi8CiwCPgBudc55swqqSAVqmZrE6POO54vRZzDppn4M7NmUZZvzSI4Pt6q9v2AD87O3h0cMrZoO4/pB3qb9t0VEpNrTElIiUcI5h5kRCjn6//1TNu3cy6B6K3hw7wP4XYD1dXr/uApCnxGaXFdEpIo6kj5q0XbrU6TG2jfM3eczPvr1afxtUDd+HXgGC5ZYBSGYD64ovAqCiIhUeyrURKJQnaRYrurbkpa3f0R27V77rYIQsHimHTeG/EL1ChARqe5UqIlEs9yltNizaL9VEMyFWJv5Af3//ikPfbiETTsLPAwoIiIVSYWaSDQ7YBWEkC+WGMKrIJzYpj7jpq7g5Ac/456MLI+DiohIRVChJhLNhk0ODxxISiV28AR8vYeHt6/8L08NS2fqb09n6ImtqJ0YHi3qnGPGslxCoao/SEhERDTqU6Ramb1iC1c9/SVtGtTiupNac1mf5ppIV0QkymjUp0gNdULrejx2VS9qJ8byx8kL6f+3z/j7lCWaSFdEpIpSoSZSjcT4fVzUoykZtwzgjZv7c1L7VDK+W0dcTPh/9e17Cg9zBhERiSa6JyJSDZkZfVrVp0+r+uwpLCI+xk9RMMQFj86gWd1EfnFyG87u3Ai/z7yOKiIih6AWNZFqbl8ftaBzXHdSa9Ztz2fUC3M4459TeW7mKnbvLfI4oYiIHIwKNZEaIj7Gzw2ntGXa705n3NDeNEiO4953FjF7xRYArSsqIhKFNOpTpAab98N2ujWrg89nvPbaS1y86NfEUqR1RUVEKpBGfYpImfRoURdfpJ/aGav+iT8U0LqiIiJRRIWaiADQ4Ob3WV+3937riu61OLhArWkiIl5RoSYiYblLabZr4X7risYQ7re2Oa+Arbs1tYeISGVToSYiYaWsK+p3AViUwd/fX8Jp//icxz9bxp5CjRIVEaksKtREJOwg64oy5DluPr0d/dql8s+PlnLaQ1N58as1BIIhrxOLiFR7GvUpImWWuXorf5+yhMw12/jlz9rzm3M6eh1JRKTKOZJRn1qZQETKLL11fV4b1Z9PFm+mR/M6AGSt28GuvUX0a5vqcToRkepHhZqIHBEz4+zOjYqfj5u6nPcXbOSMjg0Zfd7xdGpS28N0IiLVi/qoicgxGXt5T/7v58czZ802zn/0C+58dS7Z2/Z4HUtEpFpQoSYixyQh1s+o09oxffQZ3HRKW96dv4F35m3wOpaISLWgW58iUi7qJsUx5vxOjBjQmnpJcQB8kLWBFTm7+cVJbUiM83ucUESk6lGLmoiUq6Z1E4uLstkrtvDQh99z2kOf8/LXaynSlB4iIkdEhZqIVJj7BnbltVH9aVE/iTFvLuCcf01n+tIcWDUdxvWDvE37b4uIyH5UqIlIhTqhdX1eH9WfCcP64DMjZu0MeOlyQjlLyZ44lOALl0POMpj2oNdRRUSijgo1EalwZsY5XRrzwe2n0P/7BwkVFeJzQVK3zcMfzAdXBIsyvI4pIhJ1VKiJSKWJ8fuw4Rmsr9ObfBdXvAD8Xosj+POHPU4nIhJ9VKiJSOXKXUqTvKziIg2AUIgp70xi4fod3uUSEYlCKtREpHJNGY0/FAAg4Isn5Isl3oo4qXAmhUUaFSoiUpIKNRGpXMMmQ58RkJRK7OAJ+HoPh6RUkq95gV4t6wHwyMdL+XiRRoGKiJhzzusMxyw9Pd1lZmZ6HUNEykFBIMglT8xkycY8ft61Mfde3IVGtRO8jiUiUm7MbI5zLr0s+6pFTUSiSkKsn3d+eTKjz+vIZ0s2c9bD03j+yzWEQlX/j0oRkSOlQk1Eok6s38ctp7fnoztOpUeLuvz1vcVs3FngdSwRkUqntT5FJGq1Sq3F89f3ZUXOLprWTcQ5x2uZ2VzcsykJsVo7VESqP7WoiUhUMzPap6UA8O3a7Yx+Yz7n/Ws6M5fnepxMRKTiqVATkSqjT6t6vHjDiQAMfeYr7nx1Llt3Fx7mKBGRqkuFmohUKSe1b8AHvz6VW89ox9tz13P1019SHUavi4iURn3URKTKSYj187tzj+fiHs3YsmsvZkYgGGLD9gJapiZ5HU9EpNyoRU1EqqyOjVMY0L4BAP+ZsYqzH5nGE58vJxDUCgciUj14UqiZ2RAzW2hmITNLL/H62WY2x8wWRP75My/yiUjVc0mvZvzs+DQe+vB7Lnx0BnPWbPM6kojIMfOqRS0LGARMP+D1XOAi51w3YATwfGUHE5GqqVHtBJ68pg9PD09nZ0GAy8bPYsL0FbBqOozrB3mb9t8WEakCPOmj5pxbDOFh9we8/l2JpwuBBDOLd87trcR4IlKFnd25Ef3bpfLwR99zRvwSeOl6QkWFrJ94NU3yFoYXhJ/2IFw41uuoIiKHFc191AYD36lIE5EjlRwfw58u6kKHOX8hVFSIzwVJ3TYffzAfXBEsyvA6oohImVRYoWZmn5hZVimPgWU4tgvwIDDyEPvcZGaZZpaZk5NTntFFpLoYlsH6Or3Jd3EkWni+tUKLhwvUmiYiVUOF3fp0zp11NMeZWXPgLWC4c27FIc4/AZgAkJ6erkmUROSncpfSJC8Lv/04Ka4LBfnm8wy6HXeRlqESkagXVbc+zawu8B4wxjk30+s8IlLFTRkd7pMGBHzxhHyxxFsR7XI+YfNO9aoQkejn1fQcl5pZNtAfeM/MPoy8dRvQHrjHzOZGHmleZBSRamDYZOgzApJSiR08AV/v4ZCUiu/y/9IyNQnnHB9kbSAUUqO8iEQnqw5Lr6Snp7vMzEyvY4hIFTN9aQ7D//M1/dum8vDlPWhaN9HrSCJSA5jZHOdc+uH3jLJbnyIilemUDg34x+DuzMveznn/ms7b89Z7HUlEZD8q1ESkxjIzLj+hBVNuP4V2acn86uXvuO+dhV7HEhEppkXZRaTGa5Vai9dG9ueJz1fQrXltr+OIiBRToSYiAsT4fdx+Vofi5098vpzde4v49VnHERejmw8i4g399BEROYBzjvXb8xk3dQWDnpzJ8s27vI4kIjWUCjURkQOYGQ9c2o3x1/Rh3bZ8LnzsC57/cg3VYZS8iFQtKtRERA7ivK6N+fDXp9K3TSr3vr1QLWsiUunUR01E5BDSaifw3LUnsGDdDjo0SgHg+415dGyc4nEyEakJ1KImInIYPp/Ro0VdAGatyOXcf03n928tYE9hkcfJRKS6U6EmInIE+rSqx8jT2vLy12u54NH4hD9UAAAgAElEQVQZzPthu9eRRKQaU6EmInIE4mP8jPl5J1684UQKAkEGPzmLCdNXeB1LRKopFWoiIkdhQLsGfHD7qfy8WxOS42PDL66aDuP6Qd6m/bdFRI6SFmUXETkG+36G2uovKHphCBYMsKFuL5rkLcQfCkCfEXDhWI9Tikg00aLsIiKVxMwwM9yU0RAM4CdI6rb5+IP54IpgUYbXEUWkClOhJiJSDmxYBhvq9iLfxZFohQAELB4uUGuaiBw9FWoiIuUhdylN8xYWF2kAoVCQXUs+8zCUiFR1KtRERMrDlNHhPmlAwBdP0BdLvBWRvOJdAC0/JSJHRYWaiEh5GDY5PHAgKZXYwRPw9x4OSakw5Dnm/bCdS56YycocLUElIkdGoz5FRCrYrOW53PrStxQWhfj74O5c1KOp15FExEMa9SkiEkUGtG/Ae786hY6NU/jly99xT0YWe4uCXscSkSpAhZqISCVoWjeRSSP7c+MpbXj+yzW8/NVaryOJSBUQ43UAEZGaItbv4w8XdOZnxzfihNb1ANixJ0CdpFiPk4lItFKLmohIJevfLpUYv4+cvL2c/cg07n93EYFgyOtYIhKFVKiJiHikdmIMP+/amGdmrOKKp2azfnu+15FEJMqoUBMR8Uh8jJ/7Bnbl8at7sXTTLi549As+/36z17FEJIqoUBMR8diF3Zvy9m0n0ah2Ai9+uUaT44pIMQ0mEBGJAm0bJpNx60nsLQphZqzfnk+Mz0irneB1NBHxkFrURESiREKsnzqJ4RGgv3l1Huc/OoNZy3M9TiUiXlKhJiIShe4b2IW6SbFcM/ErHv10GaGQboeK1EQq1EREotBxjVKYfOtJXNyjKWM/XsqIZ79m2+5Cr2OJSCVToSYiEqVqxcfwyBU9+dugbmzZVUhcjH5ki9Q0+r9eRCSKmRlX9W3JO788mVrxMeQXBvlsyuu4cf0gbxOsmg77tkWk2tGoTxGRKsDvMwC++PhNTv76FoIWZP3TV9Ns90L8oQBMexAuHOtxShEpb2pRExGpQs5e8zDxFiSGIA13zMcfzAdXBIsyvI4mIhVAhZqISBViwzLYULc3+S6ORAsPLii0eLhArWki1ZEKNRGRqiR3KU3ysoqLNAC/hcJ91USk2lGhJiJSlUwZHe6TBgR88YR8sfhDAdyiDJ74fDm5u/Z6HFBEypMKNRGRqmTYZOgzApJSiR08AV/v4ZCUSvZZ43j002Vc9NgM5mdv9zqliJQTqw6L/6anp7vMzEyvY4iIeCpr3Q5GPj+HnF17eeCSrgxJb+F1JBEphZnNcc6ll2VftaiJiFQTXZvV4Z1fnswJrevxu9fn89CHS7yOJCLHyJNCzcyGmNlCMwuZ2U8qSjNraWa7zOy3XuQTEamq6teK47/X9WXkqW0Z0K6B13FE5Bh51aKWBQwCDjZM6RFgSuXFERGpPmL8Psac34mT2ocLtf/OWs23a7d5nEpEjoYnKxM45xZDeGmUA5nZJcBKYHclxxIRqXYKAkGenbmK9dsL+PPALlzZt6XXkUTkCERVHzUzqwXcBdzndRYRkeogIdZPxq0n0a9dKv/35gLGvLmAvUVBr2OJSBlVWKFmZp+YWVYpj4GHOOw+4BHn3K4ynP8mM8s0s8ycnJzyCy4iUs3UTYrj2WtP4JbT2/Hy12sZ9szXBENVf8S/SE1QYbc+nXNnHcVhJwKXmdk/gLpAyMwKnHOPl3L+CcAECE/PcUxhRUSqOb/PGH3e8XRtVoetuwuLF3kXkejmSR+1g3HOnbJv28zuBXaVVqSJiMjROb9bk+LtjxdtYuPOAq45sWWpfYZFxHteTc9xqZllA/2B98zsQy9yiIjUZO/MW889GVnc9cZ8CgLqtyYSjbwa9fkW8NZh9rm3ctKIiNRM/7qiJ61Tk3j0s+V8vzGPJ6/pQ9O6iV7HEpESomrUp4iIVB6fz7jznI48NawPK3J2c/HjM9i8s8DrWCJSQlT1URMRkcp3bpfGtLs1mXfnr6dhSrzXcUSkBLWoiYgI7dOS+fVZx2FmLN2Ux5g3F6jfmkgUUKEmIiL7+WrVVl75Zi2XjZ9FzvyPYVw/yNsEq6b/uC0ilUK3PkVEZD/D+rWiaZ0EXnzlBZLf/DshC7J+4lCa5GXhDwVg2oNw4VivY4rUCGpRExGRnzizUyPGp04iliJ8Lkjqtnn4g/ngimBRhtfxRGoMFWoiIlKquOveZn2d3uS7OBKtEICALx4uUGuaSGVRoSYiIqXLXUqzXQuLizQAQkH2LP3cu0wiNYwKNRERKd2U0eE+aYRb0oIWSyxF7J33JlnrdngcTqRmUKEmIiKlGzYZ+oyApFRiB0/A32c4RQn1uSf2NwwZP5v3F2zwOqFItWfOOa8zHLP09HSXmZnpdQwRkRohJ28vI5/P5Nu123nosu4MSW/hdSSRKsXM5jjn0suyr1rURETkiDRMieflm/ox8tS2nHF8mtdxRKo1FWoiInLE4mP8jDm/Ew2S4wkEQ/xpchbrtud7HUuk2lGhJiIix2RFzi7e/HYdAx+fwZw1W72OI1KtqFATEZFjcnzj2rx16wBqxcdw1YSveC3zB68jiVQbKtREROSYtU9LYfKtJ3FCm3r87vX5PPH5cq8jiVQLKtRERKRc1E2K47nr+nLDyW047biGXscRqRZUqImISLmJ9fu4+8LOdG1WB4Anp65gVe5uj1OJVF0q1EREpEJszitgwvQVXPLETGYsy/U6jkiVpEJNREQqRFpKAm/fdjKNaycw4tmveW7mKqrDJOsilemQhZqZjS6xPeSA9/5aUaFERKR6aFE/iTduGcAZHdO4951F/PndRV5HEqlSDteidmWJ7TEHvHdeOWcREZFqKDk+hgnD+nDL6e1Ib1Xf6zgiVUrMYd63g2yX9lxERKRUPp8x+rzji5+/PW89xzVK5vjGtT1MJRL9Dtei5g6yXdpzERGRwyoIBPnHB0sYPG4WHy3c6HUckah2uEKth5ntNLM8oHtke9/zbpWQT0REqpmEWD9v3DyA9mnJ3PT8HJ74fLkGGYgcxCELNeec3zlX2zmX4pyLiWzvex5bWSFFRKR6aVQ7gUkj+3Nxj6Y89OH3/HrSXEIrpsG4fpC3CVZN/3FbpAY7XB81ERGRCpEQ6+ffV/akY+MU0nK/xvfKnYSKClk/cShN8rLwhwIw7UG4cKzXUUU8o3nURETEM2bGrWe0Z0jOo4SKCvG5IKnb5uEP5oMrgkUZXkcU8ZQKNRER8d6wDNbX6UW+iyPRCgEI+OLhArWmSc2mQk1ERLyXu5QmeQuLizQAC4Vwq6Z5GErEeyrURETEe1NGh/ukEW5JCxBDDAF2f/cGhUUhj8OJeEeFmoiIeG/YZOgzApJSiR08gZj04eyJqcuERn8kxqf51aXmsuowd016errLzMz0OoaIiJSzUMjh8xnrt+ezpzBI+7RkryOJHDMzm+OcSy/LvmpRExGRqOWLtKbd9cZ8Lh03kxnLcj1OJFK5VKiJiEjU+9ugbjStk8iIZ7/mxa/WeB1HpNKoUBMRkajXvF4Sr9/cn1M6NOAPb2Xxl3cXEQxV/a47IoejQk1ERKqElIRYnhmezrUDWjNzeS4FgaDXkUQqnJaQEhGRKiPG7+Pei7uwe28RteJj2FNYxI78AE3qJHodTaRCqEVNRESqnFrx4XaGezIWcvHjM5n3w3aPE4lUDBVqIiJSZY08rS3xMT6umDCb9xds8DqOSLnzpFAzsyFmttDMQmaWfsB73c1sduT9BWaW4EVGERGJfsc1SiHj1pPo0rQOt7z4LU98vpzqMD+oyD5etahlAYOA6SVfNLMY4AVglHOuC3A6EKj0dCIiUmU0SI7nxRtOZGDPpjzzxUpydu31OpJIufFkMIFzbjGA2U+WBTkHmO+cmxfZb0slRxMRkSooIdbPv67oSfa2fNJSEnDOkbe3iNoJsV5HEzkm0dZH7TjAmdmHZvatmY0+2I5mdpOZZZpZZk5OTiVGFBGRaGRmtKifBMDjny3nosdmsHzzLo9TiRybCivUzOwTM8sq5THwEIfFACcDQyP/vNTMzixtR+fcBOdcunMuvWHDhhVwBSIiUlWd1KEBu/cWMWjcTGYu17JTUnVVWKHmnDvLOde1lMfkQxyWDUxzzuU65/YA7wO9KyqjiIhUT71b1uOtW06iSZ1ERvzna17+eq3XkUSOSrTd+vwQ6G5mSZGBBacBizzOJCIiVVCL+uFlp05q34B7MrJYnbvb60giR8yTwQRmdinwGNAQeM/M5jrnznXObTOzscA3gAPed86950VGERGp+lISYpk4Ip25P2yndYNaAARDDr/vJ4PZRKKSJy1qzrm3nHPNnXPxzrlGzrlzS7z3gnOuS+Q26UEHE4iIiJRFjN9Heuv6AHyQtYG7/zWOwGMnQt4mWDUdxvULb4tEIa31KSIiNUbjrd/wxx33YlbEiqeupHX+YvyhAEx7EC4c63U8kZ+Itj5qIiIiFaZn1l+J9wWJIUTTvCz8wXxwRbAow+toIqVSoSYiIjXHsAzW1+lNvosj0QoBCFg8XKDWNIlOKtRERKTmyF1Kk7ys4iINwG+hcF81kSikQk1ERGqOKaPDfdKAgC+ekC8WXyhAUdZb3PbSt+zeW+RxQJH9qVATEZGaY9hk6DMCklKJHTwBX+/hkJTKrN7/ZErWRoaMn83GHQVepxQpZs45rzMcs/T0dJeZmel1DBERqcKmfr+ZW1/8lpSEWP5z7Ql0blrb60hSTZnZHOdceln2VYuaiIgIcHrHNF4bNQCAIeNnsXD9Do8TiahQExERKda5aW0ybj2JIekt6JCW4nUcERVqIiIiJTWuk8C9F3chLsbHll17GTd1OaFQ1e8mJFWTCjUREZGDyJi7nn988D23vfwtBYGg13GkBtISUiIiIgfxi5Na45zjgfcXs2HHlzw9PJ0GyfFex5IaRC1qIiIiB2Fm3HBKW54c2pvFG3Zy6biZrMrd7XUsqUFUqImIiBzGeV2b8MpN/UlLSaB2gm5GSeVRoSYiIlIGPVvU5fVR/UlNjqewKMTU7zd7HUlqABVqIiIiZWRmAPx31mquffYb/vXJUqrDxPESvdR+KyIicoRGDGjNko15/OuTZazduoe/D+pOXIzaPqT8qVATERE5QnExPv45pDutUpMY+/FS1m/P56lr0qmTFOt1NKlmVP6LiIgcBTPjV2d24F9X9GTppl1s2JnvdSSphtSiJiIicgwu6dWMszo3Ijk+/Cv1h617aFE/yeNUUl2oRU1EROQY7SvSXv3mB84aO40pCzZ4nEiqCxVqIiIi5eTMTml0aVqbW176lgnTV2hEqBwzFWoiIiLlJDU5npdu7Mf5XZvw1/eXMPGF/+Ge6Ad5m2DVdBgX2RYpI/VRExERKUcJsX4eu6oXJ8csYuCi3+J8QdZNHEqTvCz8oQBMexAuHOt1TKki1KImIiJSznw+46otj5PgC+JzQVK3zcMfzAdXBIsyvI4nVYgKNRERkYowLIP1dXqT7+JItEIAAhYPF6g1TcpOhZqIiEhFyF1Kk7ys4iINIBQKsmHuRx6GkqpGhZqIiEhFmDI63CcNCPjiCfpiibci4pa+zRtzsj0OJ1WFCjUREZGKMGwy9BkBSanEDp6Av/dwQompPNnwbn7z2jy+XLnF64RSBVh1mOMlPT3dZWZmeh1DRETksAqLQrz5bTZXnNACM/M6jnjAzOY459LLsq9a1ERERCpRXIyPK/u2xMxYnbub21/5jl17i7yOJVFKhZqIiIhHFqzbwbvzN3D5+Nls2lngdRyJQirUREREPHJRj6Y8MyKd1Vt2M2jcLJZuyvM6kkQZFWoiIiIeOqNjGq+O7E9hMMTgJ2eRtW6H15EkiqhQExER8VjXZnV465YBnN25Ee0aJnsdR6KICjUREZEo0LxeEmMv70linJ+8ggAvf72W6jAzgxwbFWoiIiJR5qWv1jLmzQX8/q0sioIhr+OIh2K8DiAiIiL7u/GUtmzPD/Dk1BVs3JHP41f3pla8fmXXRGpRExERiTI+n3HXecfzwKVdmbY0hysmzGazpu+okVSoiYiIRKmhJ7bimRHp7N4bpCik/mo1kSeFmpkNMbOFZhYys/QSr8ea2X/NbIGZLTazMV7kExERiRY/O74RH99xKk3rJhIMOZZs3Ol1JKlEXrWoZQGDgOkHvD4EiHfOdQP6ACPNrHXlRhMREYkuMf7wr+unpq/g4sdm8va89R4nksriSaHmnFvsnPu+tLeAWmYWAyQChYD+dBAREQGu7tuSni3q8quXv2P8tBWavqMGiLY+aq8Du4ENwFrgn865raXtaGY3mVmmmWXm5ORUZkYRERFP1E2K43/X9+XC7k34+5Ql3DNZ03dUdxVWqJnZJ2aWVcpj4CEO6wsEgaZAG+A3Zta2tB2dcxOcc+nOufSGDRtWwBWIiIhEn4RYP49e2YuRp7bl1cxslm3e5XUkqUAVNimLc+6sozjsauAD51wA2GxmM4F0YGW5hhMREanCfD5jzPmduPrElrRKrQVA4bKpxH38fzBsMuR+D1NGh7dTGnmcVo5FtN36XAv8zMJqAf2AJR5nEhERiUr7irTpH75B6MXLCeUsJXviUIIvXA45y2Dagx4nlGPl1fQcl5pZNtAfeM/MPoy89QSQTHhU6DfAs865+V5kFBERqSpOWPJ3Ygjgc0FSt83DH8wHVwSLMryOJsfIk/UonHNvAW+V8vouwlN0iIiISBkl/uIdVk+4mkY7F5BohQAEfPHEXjDW42RyrKLt1qeIiIgcqdyltNizqLhIA/C5EKw6cLpSqWpUqImIiFR1U0bjDwWAcEtakcXidwHd+qwGVKiJiIhUdcMmQ58RkJRK7OAJxPQZDkmpLDvtMUa/Po+9RUGvE8pRsuowq3F6errLzMz0OoaIiEhUeW7mKu59ZxED2qUyflgfaifEeh1JADOb45xLP/yealETERGptq49qQ1jL+/B16u2cvn42WzcUeB1JDlCKtRERESqsUG9m/PcdX3J3pbPoHEzWb893+tIcgRUqImIiFRzJ3dowKSR/TitYxqNaid4HUeOgAo1ERGRGqBL0zr8bVA3/D5j444CPly40etIUgYq1ERERGqYf3+6lFEvzGHijFVeR5HD8GRlAhEREfHOny7qwrbdAf7y7iLWb8/nD+d3wuczr2NJKdSiJiIiUsMkxPp5Ymhvrh3QmokzVvHLV77TXGtRSi1qIiIiNZDfZ/zpos40rZvAG3PWURAIER/j9zqWHECFmoiISA1lZtx0ajtGDGhNfIyfgkCQ7XsCNK6jkaHRQrc+RUREarh9LWn3ZGRxyRMzWbJxp8eJZB8VaiIiIgLAL05uA8CQJ2cza0Wux2kEVKiJiIhIRKcmtXnzlgE0qZvAtf/5hrfnrfc6Uo2nQk1ERESKNa2byGsjB9CzZV3+NDmLHfkBryPVaBpMICIiIvupkxTL/37Rl7Vb91AnMRbnHM6hudY8oBY1ERER+YmEWD/HNUoB4LHPlnPrS99SENBca5VNhZqIiIgcUlKcnylZG7nmma/YvqfQ6zg1igo1EREROaQbTmnL41f3Yn72DgY/OYuc+R/DuH6QtwlWTf9xW8qdCjURERE5rAu7N+V/1/eldd4ckt8cSihnKdkThxJ84XLIWQbTHvQ6YrWkQk1ERETKpF/bVMbVn0Q8RfhckNRt8/AH88EVwaIMr+NVSyrUREREpMzir3ub9XV7k+/iSLRwf7WALx4uGOtxsupJhZqIiIiUXe5SmuRlFRdpAOZC4b5qUu5UqImIiEjZTRmNPxSeBDfgiydADDEuwO65bxAKOY/DVT8q1ERERKTshk2GPiMgKZXYwROI6TOc3TF1uWHPrfx60lyKgiGvE1Yr5lzVr37T09NdZmam1zFERERqJOcc46etZNPOAv50UWfMtILBoZjZHOdceln21RJSIiIickzMjJtPb4dzDjNj+eY8aifEklY7wetoVZ5ufYqIiEi5MDOCIcfI5+cw6MlZrMjZ5XWkKk+FmoiIiJQbv8945Iqe5BcGuezJWXy7dpvXkao0FWoiIiJSrro3r8sbNw+gdmIsVz/9JZ8u1vJSR0uFmoiIiJS71g1q8fqoAXRIS+HZmaupDoMXvaDBBCIiIlIhGqbE88pN/QhGBhnkFwZJiPVpVOgRUIuaiIiIVJha8THUToilIBDkmolfcXdGFkFNjFtmKtRERESkwsXH+Ojbpj4vfrWWm1+YQ0Eg6HWkKkGFmoiIiFQ4M+Ou847n3os68/HiTVzzzFds31N4+ANrOBVqIiIiUmmuPakNj1/Vm/nZO7hj0lyv40Q9DSYQERGRSnVB9yakJseRlhLvdZSopxY1ERERqXT92qbStmEyzjn+8u4ivly5xetIUcmTQs3MHjKzJWY238zeMrO6Jd4bY2bLzex7MzvXi3wiIiJSOXbmFzFtaQ7DJ37N+ws2eB0n6njVovYx0NU51x1YCowBMLPOwJVAF+A8YJyZ+T3KKCIiIhWsTlIsr4/qT7fmdbj1pW/576zVXkeKKp4Uas65j5xzRZGnXwLNI9sDgVecc3udc6uA5UBfLzKKiIhI5aibFMeLN5zIWZ0a8ae3F/LvT5Z5HSlqREMftV8AUyLbzYAfSryXHXntJ8zsJjPLNLPMnJycCo4oIiIiFSkh1s+TQ3szrF8rureo43WcqFFhhZqZfWJmWaU8BpbY5w9AEfDivpdKOVWp0xc75yY459Kdc+kNGzYs/wsQERGRShXj9/GXS7pyRsc0AL6dOpnQ4/0gbxOsmg7jIts1SIVNz+GcO+tQ75vZCOBC4Ez340qt2UCLErs1B9ZXTEIRERGJVrkLPqbT5zcQsiA/PH0VzXcvwh8KwLQH4cKxXserNF6N+jwPuAu42Dm3p8RbbwNXmlm8mbUBOgBfe5FRREREvNPgi3uI9wWJIUjajgX4g/ngimBRhtfRKpVXfdQeB1KAj81srpmNB3DOLQReBRYBHwC3Oue0GJiIiEhNMyyD9XV6k+/iSLTwUlMBi4cLak5rGni0MoFzrv0h3nsAeKAS44iIiEi0yV1Kk7ws/PbjeqBGKNxXrcslHgarXNEw6lNERERkf1NGh/ukAQFfPCFfLDEuAIsyatRi7irUREREJPoMmwx9RkBSKrGDJ+DrPRySUll39pOc9tBUnpu5yuuElcJ+HHBZdaWnp7vMzEyvY4iIiEgFKwgE+dXL3/HRok3ccno7fnduR8xKm90repnZHOdceln2VYuaiIiIVBkJsX6evKYPV/VtybipKxj9+nyKgiGvY1UYTwYTiIiIiBwtv8/466VdSUuJ59+fLuO4RinceGpbr2NVCBVqIiIiUuWYGXecfRxdmtbmtI7Vd4Ui3foUERGRKuucLo2Jj/GzdXchN/0vk/Xb872OVK5UqImIiEiVt3brHmav2MKgcbNYuinP6zjlRoWaiIiIVHk9W9Tl1VH9CTnHZU/O4pvVW72OVC5UqImIiEi10KlJbd64eQANkuO55pmvmLEs1+tIx0yFmoiIiFQbLeon8frNAzinS2M6Nk7xOs4xU6EmIiIi1Ur9WnE8dlUvGqbEEwiGmDx3HVV1gn8VaiIiIlJtvfltNre/Mpd7JmcRDFW9Yk3zqImIiEi1dXl6C1bm7uapaSvZsquQR67oSUKs3+tYZaZCTURERKotM2PMzzvRMDme+99bzNbdXzNheDp1EmO9jlYmKtREROT/27v3GDvKOozj32fb0rWXYNsFrVxLJGpJVMrGIKByMeEWqAY1JQjUSBDEW0w0ihES/EfjNcYQrEIoolwsAsWAQgRspCmyJRTaIqWUqzVSKLdCU2j78495F6eH3e0s7Z55z5nnk5x0zsw7u++z73k3v51Lx6zrnfOxg9hr6kQuXryKf7+w2YWamZmZWU7mfngfjn3/3kztLYq0ja++zvTJe9Tcq5H5ZgIzMzNrjMEi7eplT3LcT+/mgadfrLlHI3OhZmZmZo1z5Hv7mNI7ntMXLOPuR56tuzvDcqFmZmZmjTOrbzI3nH8Es/omc87CAZbcfgNceji88l94fMn/l2vmQs3MzMwaae+pvVz3pcOZ/56n6L/nPLZvWMMzl5/Btqs/Bxsehb//qO4uulAzMzOz5praO4ELdSW9PdvoiW3MeGEF47ZthtgKq2+qu3su1MzMzKzZes66ifV7zmFz7ME79DoAb/RMhJN/VnPPXKiZmZlZ0z23hpmvrHyzSAMYx/biWrWauVAzMzOzZrvt24zb/gZQHEnb3jOBnu1v+NSnmZmZWe3OvBkOOxsmzWDCaQvomXMWTJoBn72y7p6hiM57knyr/v7+GBgYqLsbZmZmZjslaXlE9Fdp6yNqZmZmZplyoWZmZmaWKRdqZmZmZplyoWZmZmaWKRdqZmZmZplyoWZmZmaWKRdqZmZmZplyoWZmZmaWKRdqZmZmZplyoWZmZmaWqa54hJSkDcCTbfhWfcBzbfg+OWpydmh2fmdvribnb3J2aHb+dmQ/ICL2qtKwKwq1dpE0UPXZXN2mydmh2fmdvZnZodn5m5wdmp0/t+w+9WlmZmaWKRdqZmZmZplyoTY6C+ruQI2anB2and/Zm6vJ+ZucHZqdP6vsvkbNzMzMLFM+omZmZmaWKRdqgKQTJD0iaa2k7wyxfaKk69L2eyUdWNr23bT+EUnHt7Pfu0uF/N+UtFrSg5L+JumA0rZtkh5Ir8Xt7fmuq5B9vqQNpYznlLadLenR9Dq7vT3fdRWy/7yUe42kF0vbOn3cr5D0rKSVw2yXpF+mn82DkuaUtnX0uEOl/Gek3A9KWirpQ6VtT0h6KI39QPt6vXtUyH60pJdKn++LSttGnDO5q5D9W6XcK9M8n562dfS4A0jaT9Jdkh6WtErS14dok9/cj4hGv4BxwGPAQcAewApgdkubLwOXpeV5wHVpeXZqPxGYlX8d02cAAAWHSURBVL7OuLozjUH+Y4BJafn8wfzp/aa6M4xx9vnAr4bYdzqwLv07LS1PqzvT7sze0v6rwBXdMO6p/x8H5gArh9l+EnAbIOBw4N5uGPdR5D9iMBdw4mD+9P4JoK/uDGOY/Wjgz0OsH9WcyfG1s+wtbU8B7uyWcU8ZZgJz0vJUYM0Qv/Ozm/s+ogYfAdZGxLqIeB24Fpjb0mYusDAtLwKOk6S0/tqI2BIRjwNr09frJDvNHxF3RcRr6e0yYN8293GsVBn74RwP3BERGyPiBeAO4IQx6udYGG3204Fr2tKzNoiIJcDGEZrMBa6KwjLgnZJm0vnjDuw8f0QsTfmgu+Z8lbEfzq78vsjCKLN31ZwHiIj/RMT9afkV4GFgn5Zm2c19F2rFID1dev8Mbx24N9tExFbgJWBGxX1zN9oMX6T4a2NQr6QBScskfWosOjiGqmY/LR0CXyRpv1Hum6vK/U+numcBd5ZWd/K4VzHcz6fTx/3taJ3zAdwuabmkc2vq01j7qKQVkm6TdEha15ixlzSJogi5obS6q8ZdxSVMhwL3tmzKbu6Pb8c3yZyGWNd6K+xwbarsm7vKGSR9HugHPlFavX9ErJd0EHCnpIci4rEx6OdYqJL9FuCaiNgi6TyKI6vHVtw3Z6Pp/zxgUURsK63r5HGvopvnfGWSjqEo1I4qrT4yjf3ewB2S/pWO1HSL+yke77NJ0knATcDBNGvsTwHuiYjy0beuGXdJUyiK0G9ExMutm4fYpda57yNqRVW8X+n9vsD64dpIGg/sSXH4uMq+uauUQdInge8Bp0bElsH1EbE+/bsOuJviL5ROsdPsEfF8Ke9vgMOq7pu50fR/Hi2nQDp83KsY7ufT6eNemaQPAr8F5kbE84PrS2P/LHAjnXe5x4gi4uWI2JSWbwUmSOqjQWPPyHO+o8dd0gSKIu33EfGnIZrkN/fbeSFfji+Ko4rrKE7tDF4gekhLmwvY8WaC69PyIex4M8E6Ou9mgir5D6W4iPbglvXTgIlpuQ94lA66uLZi9pml5U8Dy9LydODx9DOYlpan151pd2ZP7d5HcRGxumXcSzkOZPgLyk9mxwuK/9kN4z6K/PtTXHN7RMv6ycDU0vJS4IS6s+zm7O8e/LxTFCNPpc9BpTmT+2uk7Gn74IGIyV047gKuAn4xQpvs5n7jT31GxFZJXwH+SnFXzxURsUrSJcBARCwGLgd+J2ktxQd4Xtp3laTrgdXAVuCC2PH0UPYq5v8xMAX4Y3EPBU9FxKnAB4BfS9pOcXT2hxGxupYgb0PF7F+TdCrF+G6kuAuUiNgo6QfAfenLXRI7nibIWsXsUFxQfG2k31RJR487gKRrKO7u65P0DHAxMAEgIi4DbqW4+2st8BrwhbSto8d9UIX8F1Fch3tpmvNbo3hI9buAG9O68cAfIuIvbQ+wCypk/wxwvqStwGZgXvr8DzlnaojwtlXIDsUfpLdHxKulXTt+3JMjgTOBhyQ9kNZdSPGHSbZz308mMDMzM8uUr1EzMzMzy5QLNTMzM7NMuVAzMzMzy5QLNTMzM7NMuVAzMzMzy5QLNTMzM7NMuVAzMzMzy1Tj/8NbM7OhSPo+cAbFg5ifA5ZHxE/q7ZWZNY0LNTOzFpL6gdMoHp82nuJB3ctr7ZSZNZILNTOztzoKuDkiNgNIuqXm/phZQ/kaNTOzt1LdHTAzAxdqZmZD+QdwiqReSVOAk+vukJk1k099mpm1iIj7JC0GVgBPAgPAS/X2ysyaSBFRdx/MzLIjaUpEbJI0CVgCnBsR99fdLzNrFh9RMzMb2gJJs4FeYKGLNDOrg4+omZmZmWXKNxOYmZmZZcqFmpmZmVmmXKiZmZmZZcqFmpmZmVmmXKiZmZmZZcqFmpmZmVmm/gdPvuClVpbSswAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axes = plt.subplots(nrows=1, ncols=1,figsize=(10,6))\n",
    "\n",
    "ax=axes\n",
    "ax.plot(gs,exact_Es,'.--',label='exact E')\n",
    "ax.plot(gs,dmrg_Es,'X',label='DMRG E')\n",
    "ax.set(xlabel='g', ylabel='E',\n",
    "       title='Ground state Energies')\n",
    "ax.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We see that the results match perfectly so I assume the algorithm is working (in principle at least)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}