From 6fcd836a55c47c3eea2c8f7a489a13317682b07d Mon Sep 17 00:00:00 2001 From: Georg Wolfgang Winkler <winklerg@itp.phys.ethz.ch> Date: Tue, 31 May 2016 11:41:32 +0200 Subject: [PATCH] added ex13 solution --- exercises/ex13_solution/heisenberg_dmrg.ipynb | 331 ++++++++++++ ...rg_imaginary_time_evolution_solution.ipynb | 507 ++++++++++++++++++ 2 files changed, 838 insertions(+) create mode 100644 exercises/ex13_solution/heisenberg_dmrg.ipynb create mode 100644 exercises/ex13_solution/heisenberg_imaginary_time_evolution_solution.ipynb diff --git a/exercises/ex13_solution/heisenberg_dmrg.ipynb b/exercises/ex13_solution/heisenberg_dmrg.ipynb new file mode 100644 index 0000000..f0eb125 --- /dev/null +++ b/exercises/ex13_solution/heisenberg_dmrg.ipynb @@ -0,0 +1,331 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "import numpy as np\n", + "from numpy import transpose as tr, conjugate as co\n", + "from scipy.linalg import expm, svd\n", + "from scipy.sparse.linalg import eigsh, LinearOperator\n", + "import math\n", + "import matplotlib.pyplot as plt\n", + "%matplotlib inline\n", + "plt.rcParams['figure.figsize'] = 16, 9" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Useful functions" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "def dot(A,B):\n", + " \"\"\" Does the dot product like np.dot, but preserves the shapes also for singleton dimenstions \"\"\"\n", + " s1 = A.shape\n", + " s2 = B.shape\n", + " return np.dot(A,B).reshape((s1[0],s2[1]))\n", + "\n", + "def left_canonize(M1,M2,return_S = False):\n", + " \"\"\" Left normalizes M1 into A matrix, M2 looses its canonization\"\"\"\n", + " s, da, db = M1.shape\n", + " U, S, Vh = svd(M1.reshape((s*da,db)))\n", + " U = U.reshape((s,da,U.shape[1]))[:,:,:db]\n", + " M2 = np.tensordot(dot(np.diag(S[:db]),Vh[:db,:]),M2,axes=((1),(1)))\n", + " if return_S:\n", + " return U, M2, S\n", + " else:\n", + " return U, M2\n", + "\n", + "def right_canonize(M1,M2,return_S = False):\n", + " \"\"\" Right normalizes M2 into B matrix, M1 looses its canonization \"\"\"\n", + " s, da, db = M2.shape\n", + " U, S, Vh = svd(M2.transpose((1,0,2)).reshape((da,s*db)))\n", + " Vh = Vh.reshape((Vh.shape[0],s,db)).transpose((1,0,2))[:,:da,:]\n", + " M1 = np.tensordot(M1,dot(U[:,:da],np.diag(S[:da])),axes=((2),(0)))\n", + " if return_S:\n", + " return M1, Vh, S\n", + " else:\n", + " return M1, Vh\n", + " \n", + "def add_site_to_R(Rexp,M,W):\n", + " Rexp = np.tensordot(M.conj(),Rexp,axes=((2),(0)))\n", + " Rexp = np.tensordot(W,Rexp,axes=((0,3),(0,2)))\n", + " Rexp = np.tensordot(M,Rexp,axes=((0,2),(0,3)))\n", + " Rexp = np.transpose(Rexp,(2,1,0))\n", + " return Rexp\n", + "\n", + "def add_site_to_L(Lexp,M,W):\n", + " Lexp = np.tensordot(M.conj(),Lexp,axes=((1),(0)))\n", + " Lexp = np.tensordot(W,Lexp,axes=((0,2),(0,2)))\n", + " Lexp = np.tensordot(M,Lexp,axes=((0,1),(0,3)))\n", + " Lexp = np.transpose(Lexp,(2,1,0))\n", + " return Lexp\n", + "\n", + "def Hpsi(L,W,R,M):\n", + " \"\"\"\n", + " Calculates the \"matrix-vector product\" at arbitrary site from \n", + " L ... L-expression, belonging to the \"matrix\"\n", + " W ... W-array at site l, belonging to the \"matrix\"\n", + " R ... R-expression, belonging to the \"matrix\"\n", + " M ... M-array, the \"vector\"\n", + " \"\"\"\n", + " if len(M.shape)==1:\n", + " flatten = True\n", + " M = np.reshape(M,(s,L.shape[0],R.shape[0]))\n", + " elif len(M.shape)==3:\n", + " flatten = False\n", + " else:\n", + " raise ValueError('Unknown format for M')\n", + " hpsi = np.tensordot(L,M,axes=((2),(1)))\n", + " hpsi = np.tensordot(W,hpsi,axes=((1,2),(2,1)))\n", + " hpsi = np.tensordot(hpsi,R,axes=((1,3),(1,2)))\n", + " if flatten:\n", + " return hpsi.flatten()\n", + " else:\n", + " return hpsi" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Definitions and construction of MPO" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "J = 1.\n", + "L = 40 # Length of chain\n", + "s = 2 # Local dimension of Hilbert space\n", + "chi = 60 # The maximum matrix dimension, from which on the matrices are truncated\n", + "nmax = 10000 # Maximum number of iterations\n", + "\n", + "## Generate the MPO for the Heisenberg model\n", + "W = []\n", + "## some operators\n", + "Z = np.zeros((2,2,1,1))\n", + "\n", + "Id = np.zeros((2,2,1,1))\n", + "Id[0,0,0,0] = 1.\n", + "Id[1,1,0,0] = 1.\n", + "\n", + "Sz = np.zeros((2,2,1,1))\n", + "Sz[0,0,0,0] = +0.5\n", + "Sz[1,1,0,0] = -0.5\n", + "\n", + "Sp = np.zeros((2,2,1,1))\n", + "Sp[0,1,0,0] = 1.\n", + "\n", + "Sm = np.zeros((2,2,1,1))\n", + "Sm[1,0,0,0] = 1.\n", + "\n", + "W.append(np.concatenate((Z,J/2*Sm,J/2*Sp,J*Sz,Id),axis=3))\n", + "for i in range(1,L-1):\n", + " Wi = np.concatenate((\n", + " np.concatenate((Id,Z,Z,Z,Z),axis=3),\n", + " np.concatenate((Sp,Z,Z,Z,Z),axis=3),\n", + " np.concatenate((Sm,Z,Z,Z,Z),axis=3),\n", + " np.concatenate((Sz,Z,Z,Z,Z),axis=3),\n", + " np.concatenate((Z,J/2*Sm,J/2*Sp,J*Sz,Id),axis=3),\n", + " ),axis=2)\n", + " W.append(Wi)\n", + "W.append(np.concatenate((Id,Sp,Sm,Sz,Z),axis=2))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Random guess as groundstate" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "## random starting configuration\n", + "M = [np.random.standard_normal((s,min(chi,s**i,s**(L-i)),min(chi,s**(i+1),s**(L-1-i)))) for i in range(L)]\n", + "\n", + "## need a right canonized MPS\n", + "for i in range(L-1,0,-1):\n", + " M[i-1], M[i] = right_canonize(M[i-1],M[i])\n", + "## normalization, throw away left part\n", + "_, M[0] = right_canonize(np.ones((1,1,1)),M[0])" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### DMRG sweeps" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + ".......................................\n", + ".......................................\n", + "0 dE = [ 17.54144864]\n", + ".......................................\n", + ".......................................\n", + "1 dE = [ 1.52946031e-05]\n", + ".......................................\n", + ".......................................\n", + "2 dE = [ 9.36516908e-06]\n", + ".......................................\n", + ".......................................\n", + "3 dE = [ 2.37690756e-10]\n", + "Converged after 3 sweeps!\n", + "Ground-state energy: [-17.5414733]\n" + ] + } + ], + "source": [ + "## get initial R expression\n", + "Rexp = [np.ones((1,1,1))]\n", + "for i in range(1,L):\n", + " Rexp.append(add_site_to_R(Rexp[i-1],M[-i],W[-i]))\n", + " \n", + "enold = 0\n", + " \n", + "for n in range(nmax):\n", + " ## right sweep\n", + " Lexp = [np.ones((1,1,1))]\n", + " for i in range(L-1):\n", + " print('.', end=\"\")\n", + " dD2 = s*Lexp[i].shape[0]*Rexp[L-1-i].shape[0]\n", + " Hop = LinearOperator((dD2,dD2),matvec= lambda M: Hpsi(Lexp[i],W[i],Rexp[L-1-i],M))\n", + " en,V = eigsh(Hop,k=1,v0=M[i].flatten(),tol=1e-3 if n < 2 else 0)\n", + " M[i] = V.reshape((s,Lexp[i].shape[0],Rexp[L-1-i].shape[0]))\n", + " M[i], M[i+1] = left_canonize(M[i],M[i+1])\n", + " Lexp.append(add_site_to_L(Lexp[i],M[i],W[i]))\n", + " print('')\n", + "\n", + " ## left sweep\n", + " Rexp = [np.ones((1,1,1))]\n", + " for i in range(L-1,0,-1):\n", + " print('.', end=\"\")\n", + " dD2 = s*Lexp[i].shape[0]*Rexp[L-1-i].shape[0]\n", + " Hop = LinearOperator((dD2,dD2),matvec = lambda M: Hpsi(Lexp[i],W[i],Rexp[L-1-i],M))\n", + " en,V = eigsh(Hop,k=1,v0=M[i].flatten(),tol=1e-3 if n < 2 else 0)\n", + " M[i] = V.reshape((s,Lexp[i].shape[0],Rexp[L-1-i].shape[0]))\n", + " M[i-1], M[i] = right_canonize(M[i-1],M[i])\n", + " Rexp.append(add_site_to_R(Rexp[L-1-i],M[i],W[i]))\n", + " print('')\n", + " \n", + " # === check convergence ===\n", + " print(n,'dE = ',abs(en-enold))\n", + " if abs(en-enold) < 1e-8:\n", + " print(\"Converged after \"+str(n)+\" sweeps!\")\n", + " print(\"Ground-state energy: \" + str(en))\n", + " break\n", + " enold = en\n", + " # === check convergence ===\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Measure Entanglement Entropy" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAA7MAAAImCAYAAACBwt0rAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzsvXlsXVl23vsdkhooaqZGUhJ5Rc3ULIqiFDupsi24H9xw\nOS+J0Y2W7cJ7sJMYVpB2ynD8+jFsRrBjB4KBQI6Tf/wswKVu+8XpBLaBtFvPjYLdqFJpnseirjhq\nnkUNpMjz/li1zaurO5xhn7P3vvp+gFCly3uPFu8995z97fWttTzf90EIIYQQQgghhLhElekACCGE\nEEIIIYSQsFDMEkIIIYQQQghxDopZQgghhBBCCCHOQTFLCCGEEEIIIcQ5KGYJIYQQQgghhDgHxSwh\nhBBCCCGEEOcwImY9z/uK53mXPc+76nnebxb4+TLP8/4/z/POeJ73Q8/zGkzESQghhBBCCCHETry0\n58x6nlcF4CqAnwQwBOAYgK/5vn855zn/L4C/8H3/Y8/z3gPwf/i+/4upBkoIIYQQQgghxFpMZGbb\nAVzzfb/X9/1RAH8K4IO856wD8EMA8H3/kwI/J4QQQgghhBDyDmNCzDYC6M/5+8CXj+VyGsA/AQDP\n8/53ANM9z5uTTniEEEIIIYQQQmzHhJj1CjyW73X+DQDveZ53AsCPAxgE8DrpwAghhBBCCCGEuEGN\ngX9zAMCynL8vgdTO/j2+79/ERGa2DsA/8X3/af6BPM9Lt+CXEEIIIYQQQkiq+L5fKCFqJDN7DMAK\nz/OaPM+bDOBrAP4i9wme59V7nqcC/i0A/0+xg/m+zz8V9qerq8t4DPzDz5R/+Lm+i3/4mVbmH36u\nlfeHn2ll/uHnWvhPKVIXs77vjwH4NQA/AHABwJ/6vn/J87xuz/O++uXT3gNwxfO8ywAWAPjttOMk\nhBBCCCGEEGIvJmzG8H3/+wBW5z3WlfP//x3Af087LkIIIYQQQgghbmDCZkxISd577z3TIRDN8DOt\nTPi5Vh78TCsTfq6VBz/TyoSfa3i8cj5km/E8z3c5fkIIIYQQQgghxfE8D75FDaAIIYQQQgghhJBY\nUMwSQgghhBBCCHEOillCCCGEEEIIIc5BMUsIIYQQQgghxDkoZgkhhBBCCCGEOAfFLCGEEEIIIYQQ\n56CYJYQQQgghhBDiHBSzhBBCCCGEEEKcg2KWEEIIIYQQQohzUMwSQgghhBBCCHEOillCCCGEEEII\nIc5BMUsIIYQQQgghxDkoZgkhhBBCCCGEOAfFLCGEEEIIIYQQ56CYJYQQQgghhBDiHBSzhBBCCCGE\nEEKcg2KWEEIIIYQQQohzUMwSQgghhBBCCHEOillCCCGEEEIIIc5BMUsIIYQQQgghxDkoZgkhhBBC\nCCGEOAfFLCGEEEIIIYQQ56CYJYQQQgghhBDiHBSzhBBCCCGEEEKco8Z0AIQQQkglkM32orPzIAYH\nx9HYWIV9+z5EJtNkOqxAuBw7IYSQdxfP933TMUTG8zzf5fgJIYRUBtlsL3bvPoCenm4AdQCG0dLS\nhcOH91ovCl2OnRBCSOXjeR583/cK/Yw2Y0IIIdaQzfZiz55uvP9+F/bs6UY222s6pEB0dh7MEYMA\nUIeenm50dh40GFUwXI4dcPecIYQQEh/ajAkhhFhBoQzhkSNuZAgHB8cxIQYVdRgaGjcRTihcjt3l\nc4YQQkh8mJklhBBiBS5nCBsbqwAM5z06jIYG+2+zLsfu8jlDCCEkPvbfqQghhLwTuJwh3LfvQ7S0\ndGFCFErd6b59HxqLKSgux+7yOUMIISQ+tBkTQgixgokMYa44cSNDmMk04bvf3Yt/8A/2I5MZx+ho\nlTNW10ymCQcP7sVP/dR+LF06DqAKP/iBG7G7fM4QQgiJD7sZE0IIsQLXu+r+xV8Af/iHwO/9HvDz\nPw9cuWI6ouD82Z8B3/kO8K1vAb/yK8Dp06YjCobr5wwhhJDylOpmTDFLCCEVhsszQ7PZXmzYcBBT\npoxj2bIqfO977sT+rW8BNTXAv/t3wJw5QDYL1NebjioY//pfA4sWAd/8JjB3LnDnDlCX7961lJ6e\nXqxdexDTp49j1aoqfPe77pwzLn9XCSEkLUqJWdqMCSGkgnC9u+u0aU2YPLkLBw4A/+2/AZmM6YiC\nc+QI8NFHQHU1sH27/P1nfsZ0VMH47DPgP/5HYMoUYONG4Phx4B/9I9NRBePFiyY0NXXhX/0r4MIF\nd84Z17+rhBBiAywqIYSQCsL17q7HjokQ3LVLBJYr5puxMYm9vV3+vnOnxO8CL18C588DbW3y944O\nEeKu8Nln8n63twNHj5qOJjiuf1cJIcQGKGYJIaSCcL2769GjIkqamoCqKuDGDdMRBePiRWDx4glb\n8c6d7gjCkyeBNWsmbMUuitldu4BNm6RO+cUL0xEFw/XvKiGE2ADFLCGEVBAuzwwFRMxu3w54ngjC\nTz81HVEwPv9cRKCio0MytWNj5mIKispsKpSYdSUr/umnEv/UqcDatcCpU6YjCobr31VCCLEBXjEJ\nIaSCcHlmqO9PiFnALavukSPAjh0Tf6+vl4ZKFy6Yiyko+WJ22TL5LPr7zcUUlAcPgKEhYP16+btL\nVmOXv6uEEGILFLOEEFJBZDJNOHx4L2pr92Phwi6sW7ffmYYyPT3A9Oli1wXcE7O5mVnAjfh9X2LM\njd3z3LEaHzkimx/V1fJ3l8RsJtOEv/zLvaiu3o/587vQ3u7Od5UQQmyBYpYQQiqMadOaMHVqF/74\nj7uxYEGXM4tjVS+r2LYNuHQJGM53YlrG48dS27thw5uPd3TYL2b7+4HRUWD58jcfd0XMfvqp1Msq\nXBKzAHDvXhO2bu1CZ2c3Nm9257tKCCG2QDFLCCEVxokTwNatIkiOHwdevzYdUTDyxezUqTIm5tgx\nczEF4dgxeb8nTXrzcRcys8pi7OVN73NFzOZbpFevlhm59++biykMR45I/Nu2yfeWEEJIOChmCSGk\nANlsL/bs6cb773dhz55uZLO9pkMKzPHjMmZlzhxg6VIZu+ICuaNtFC4Iws8/f7NeVrF+PXDzpt3C\nKl8MKtragDNngJGR9GMKihqHlGuRrq6W2G3fAFGo93/zZumI/eqV6YiC4/I1khBSOVDMEkJIHtls\nL3bvPoBDhz7CJ59049Chj7B79wFnFmsnTkimB3BDDAJidT1zZiJuhQvxF6qXBURYbd9ud4azmJid\nPh1YsUI+E1s5fx5oaADmzn3z8fZ2N8SsqlfeuROYNg1oaXFn48n1ayQhpHKgmCWEkDw6Ow+ip6cb\nEzMg69DT043OzoMGowpOrph1oW4TkEV8czMwY8abj+/aJfHbOibG94uLWcBuMf7ypbzvbW2Ff267\n1Ti/XlbhSt1s75e6b9ky+a9LVmPXr5GEkMqBYpYQQvIYHBzHxCJNUYehoXET4YTi1i3g+XMgk5G/\n2yymcsmvl1UsWSK1s198kX5MQbh+XeJrbCz8c5vf/xMngDVrgLr8U/1LbBezxbLKSszaugGiyK9X\nbmtzR8y6fI0khFQWFLOEEJJHY2MVJmY/KobR0GD/JVNlZdUCed064O5d+WMzxcQsYLcg/Pzz4llZ\nQH527JjUd9qGaj5UDFfFbGOjWLz7+tKPKQz58W/bJvXuLuDyNZIQUlnwqkMIIXns2/chWlq6MLFY\nG0ZLSxf27fvQWExBybUYA0BVlYjEzz83F1MQjh6V+tJCKKuxjRw5Urj5k6K+Xubm2lgLWUwMKlav\nluZVNm6EqA2adeve/pnnyblku9U4fzNh0yYZReVCEyiXr5GEkMqCYpYQQvLIZJpw+PBeTJ++HwsX\ndmHNmv04fHivEzMgT5x4uwbS5swmADx9Knbd/Dmtip07pT7SRkrVyyp27rQvw5nbfKgYNm+EfPaZ\nbCJUFVnF2F43++IFcOHCmxtP06ZJ0y0bNz7yyWSa8L3v7UV19X7Mn9+Fjg53rpGEkMqCYpYQQgow\nd24TfL8Lf/In3aiv73JmkXb8uHsdgU+elHmykycX/vmWLVIz+/RpunGVQwmSrVtLP8/G97+/X+YP\nq9rqYthqNS4nxG0XsydOAK2tQG3tm4+7ZDW+d68JHR1d+I3f6EZHhzvXSEJIZUExSwghBTh7VjKF\nO3YAp07J6BjbuXVLOtQ2N7/5+I4dskB+/dpIWGUpVS8LiMjdvNk+cXLqFLB2rWTUSmGjmM1vPlQM\nV8VsW5tskth6zn/2WeGMvksdjZXFfvNm4PRp09EQQt5VKGYJIaQAp0/LIm3mTGD5chG3tnPihGQJ\n8wXKnDnSFMdW+2I5MQtI3axtVuNyzZ8Ura2y0XD/fvIxBaWYmMpnxw77GliNjsrmTKla5TlzZAbt\npUvpxRWGYmLcJTGrzv9Nm+R6aXv3aEJIZUIxSwghBVBiFpBFs43ZqXyOHy8+M9TG7KDi2LHyYtbG\n+Ms1f1JUV0tDIpvOoXKZTUV9PbBgAXD5cvIxBeXsWXEfzJ5d+nm2Wo1L1Su70gRKzVfesUPOj2nT\n7O8eTQipTChmCSGkALlitqPDziY4+eR3Ms7FxiZEAHD7NvDkiTS+KYWKf9yiMZZBmj8pOjrsEeMv\nX0qWvtjGRz62WY2DCnFbxawSfU0FSkxVE6hz59KNKSy9vbJJs3Sp/J1WY0KIKShmCSGJkc32Ys+e\nbrz/fhf27OlGNttrOqRAjI5KdkR113UlM1tOzNoipnI5dkyyluVqNxcvBmbNAq5eTSeucty8CTx7\nVl6EK2x6/0+cANasAerqgj3fRjG7a1f557W3y/llG+XqlV2wGqusrPodXBOzrt6bCCFvU2M6AEJI\nZZLN9mL37gPo6ekGUAdgGEeOdDkxvuHyZcmaqMY+69ZN1DzW15uNrRg3bxZu/qRYu1ayoPfuAfPm\npRpaSYLUyyrUiJ41a5KNKQiqXrCcCFd0dEzUnlZXJxtbOYJmNhUdHcB//a/JxROWTz8FOjvLP2/T\nJuDKFek6nd812CTl6pXb2qRk4J//8/RiCkt+vfjmzcB3vmMunjC4fG8ihLwNM7OEkETo7DyYs1gA\ngDr09HSjs/OgwaiCkWsxBkR8tLXZaVlUqKxsMXFVXS2i0aYMGxBezNqS3QxaL6uor5fssg1NuMKK\n2Y0bZQ7wkyfJxRSUW7eAx4+BVavKP3fqVNmIOnUq+bjCUO79dykzq3ApM+vyvYkQ8jYUs4SQRBgc\nHMfEYkFRh6Ehi4oei5AvZgH762ZLWYwVNolBQJrIHD0qNuMg2BR/mHpZhQ3xl2o+VIxJk2TWrw3z\nT1VWsyrg6sW2ulk1m7hUvfKmTeIOefkyvbjC8OqVNOHK/R1aWsT18eiRubiC4vK9iRDyNhSzhJBE\naGysAjCc9+gwGhrsv+wUErO2182W6mSssEFM5dLTA0yfDixaFOz5mzYBN26YXzC/fi2bB0Ezygob\n3v/+frE6ZzLhXmdL3WzQelnF9u12idmTJyVbXMr2XFsLrFxpbxOoM2ekVnz69InHqqokg3/mjLm4\nguLyvYkQ8jb85hJCEmHfvg/R0tKFiUXDMFpaurBv34fGYgqC7xcXs0eP2tVNN5cgmdmODhG9tswM\nDWMxBiRDuG2b+Qz5+fPSxbXcaJh8bBCz5ZoPFcMWMfvpp+GyyrZlZoNmxW22Ghebr+yK1djVexMh\npDAUs4SQRMhkmnD48F4sXLgf06Z1YfPm/U402BgYACZPBhYufPPxRYuAmTOBa9fMxFWKmzfF+ldo\n1Ecuc+YAjY121G0C4cUsYIcgLLaYL0drq9R83runP6aghLUYK5SY9X39MQVlZETEUphzZvVq4M4d\nad5mA5UgZovVi7siZjOZJnz/+3tRUyP3pp/8STfuTYSQwlDMEkISI5Npwvh4F7q7u9HS0uXEYqFQ\nVlZha93siRNiMQ6SbbNp3qmrYjZs8yeFasJl8hyKKmaXLJHM+I0b2kMKzOnTUps5Y0bw16jmbTaM\n6FH1ykE2Qtra7BWzrmdmAWBkpAnLlnXhF3+xGx984Ma9iRBSGIpZQkhi3Lkj9YUffGDHYjIIpcSs\nrXWzx4+XtxgrbBCDgMzyPXs2eNyKnTvlMzBp947S/Elh8v1/+VKy8mHfc4Vpq/Gnn4arl1XYYjXu\n65Pzttj4rFw2brSzCdTdu/Kn0His9etlFNLISPpxheXYMTkv1q+3x6lCCIkGxSwhJDHOn5fFwooV\nwNOnMufUdlzNzLomZs+fl0V9mCwbACxYAMyfD1y8mEhYZXn4UKzora3RXm/y/T9xQkRIXX4j14CY\nFrNRs8q2iNkjR4LXK9vaBEq5KQp1k66tlcZily6lH1dYVBd1illC3IdilhCSGErMep49Vr9ylBKz\nW7ZItuT583RjKoeyGQdh3TrZVDBZtwlEsxgrTArCo0flva6pifb6jg75HphowhVVDCpcF7Mm632B\n8PHbaDUuZ7F3xWp87JiI2dZWuU+ZPjcIIdGhmCWEJIYSs4AsHGwXs48fi9BbsaLwz6dOlcXPyZPp\nxlWKoSGx9S1bFuz5NtRtAvHE7K5d5sTs559Hq5dVzJ0LNDSYyQbFFbPbtkncJqyvAwMyo7XYd7MU\njY2y+dDXpz+uMIR9/7dts2O2by7lmp+5IGZHRuQ83roVmDcPmDZNzi9CiJtQzBJCEiNXzLa32y9m\nz54FNmwQwVcM2+pmlcU4zKgVG5pAxc3Mfvqp3niCEqdeVmHi/VfNh+KI2WnTpDuwCbESdaQQIK8x\nbTWOUq9sW0fj8XF5D13PzJ49K5siym5PqzEhbkMxSwhJBN+XBYKqLdy+3Q6rXylOnSpuMVbYVjcb\nxmKsMF03+/QpcP26bBxEYf16YHAQePBAb1zl8P34mVnAzPvf1yfW5kwm3nFMWY3jCnHTYvbECWDt\nWtkQCMrGjdJQyZYmUFeuAPX1UrNejE2bRMzafJ1X9bIKillC3IZilhCSCP39wPTpsvgBxFo5ZYrZ\n0R7lKFUvq7AtMxumk7HCZN0mIDbtjRtl1EsUampkMZr253DtmjSsWrw43nFMiNkwzYdK4aqYVZtp\npogSf20tsGqVZBJtIMhIqgULRLD39qYTUxRUvayCYpYQt6GYJYQkQq7FWGF73WwQMdvSIrV7g4Pp\nxFSOMJ2MFSbrNoF4FmOFibrZcvWCQWltBW7dSrcJV1wxqDAhZl++FEGXK0DC0tYmmyivX+uLKwxq\nMyEsNlmNg57/tluNKWYJqSwoZgkhiVBIzNpcNzsyIp2Ky1lfPU+yEzZYjYeGZF5r0OZPuah5rSbQ\nIWZN1M0GyUwFQTXhSvP91yVmV64EnjwRMZ4WJ0/GGykEAHPmyAaOibExceqVbRKzQc9/m8Xs06dA\nNvvmdX7dOjkvTDlVCCHxoJglxHKy2V7s2dON99/vwp493chmLfZv5VAsM2vDvMdCXL4MNDUFq2mz\npW72+HHJOEWxjpoebxNXzJqwSuto/qRI8/1XzYfC1lYXwsRmji4hbqputr9fztPm5vCvbWuzo6Px\n8LDY7Ms5VwC7xezJkyJkJ0+eeGzGDGDhQqnjdwFX1wSEJAXFLCEWk832YvfuAzh06CN88kk3Dh36\nCLt3H3Di5lVIzLa1SZMlG3fAg1iMFabnbSqiWIwVpjoa374t2ZEoI1Zyqa+X2tW07IHPn8uGx5Yt\neo6XppiN0nyoFGmf/66L2c8+k/csyqbTxo3A1avmm0AdPy4icMqU8s+1WczmW4wVrliNXV4TEJIU\nFLOEWExn50H09HQDUP66OvT0dKOz86DBqMozNiYL/3Xr3nx8zhzZAb982UxcpQgjZrdvN1t/p4gj\nZlXd5v37emMqh1pMxm1EBEjdbFpW4xMnZME7daqe43V0iEBIY2NHlxhUpClmfV8+41274h/LpJiN\n+v5PnWpHE6gw9eItLXJdefgw2Zii4LqYdXVNQEiSUMwSYjGDg+OYuGkp6jA0NG4inMD09ACLFkk3\n43xsrZsNI2ZnzwaWLjW7+PH9CZtxFKqrzXQE1mExVqSZ3dTV/EmRZhMulRnURXu7nHtpbOb09cl8\n0ygW3Xw2bZLxMs+fxz9WGKI2f1LYYDUOUy9eVSUZ5TNnko0pCseOFb7+uCJmXV0TEJIkFLOEWExj\nYxWA4bxHh9HQYPdXt5DFWGFj3azvhxOzgPkmUENDktVbujT6MUzUzboqZnU1f8oljfjjNB8qxpw5\nQGMjcOGCvmMWQ8WuI5M/dao4Ek6din+soLx8CZw7F69e2XQTKN8PXy9uo9X47l3JGK9a9fbPXBGz\nrq4JCEkSnv2EWMy+fR+ipaULEzevYbS0dGHfvg+NxRSE8+eLdwW2cTxPf7/Ugi1cGPw1putmlcU4\nziI/7Y7Gvi9iNs6IlVzWrQPu3JFFatLobP6kSEPM9vXJpkcmo/e4aTVB++wzPRZjRdrOkJMn49cr\nmxazAwOShQ+THbdRzKqZ3FUFVr6rV0sDqFev0o8rDK6uCQhJEopZQiwmk2nC4cN7MW/eftTXd2H1\n6v04fHgvMpkm06GVpFRmdssW4OJFuxYNYbOygPnM7IkT8bvT7tiRbkfgnh7pHLpokZ7jVVfL75C0\nIBwYkNFNugVhGmJWZ2Yzl7Q2cz79VG9WOe26WR0Wb9UE6sULPTGFRVnsw5xDNorZYvWygGTtm5vl\nfbaZTKYJ3/veXlRX78fs2V147z031gSEJAnFLCGW09zchLGxLhw40I2FC7ucuGmVErPTpsmsStMN\nTXKJImbXr5eM7qNHycRUDpVliEN9vQjLNOyigF6LsSINQRhlMR+Edeuku/O9e3qPm4tui7EiDTH7\n4oVsfMU9z3MxIWbjvv9Tp0rm0NQ1M4rFfv16qU8eGUkmpigUq5dVuGI1fviwCe3tMpbngw/cWBMQ\nkiQUs4RYzr17Ys/86Z8W0TVueZ+HV69kKP3q1cWfY1vdbBQxW1MDbN1qxjLt+/E6GeeSZt2pToux\nIo34k6iXBdJpwhW3+VAx1q8XC3OSmznHj0uNa22tvmOuXj1RO5k0OuuVTVqNozQ/q60Fli+XzQgb\nCFLi4IqYPXlS7j2rV8uGASHvOhSzhFjOlSty05o7V/709JiOqDRXrsgiJncofT621c1GEbOAubrZ\noSHZ1IjT/EmRtpjVnZnt6JBF/uio3uPmkkS9rCLJ9//lS1mcx7WjF6KmRgRWkt9j3fWygNRLJh23\nQtWa6rCnt7WZEbOjo9IwK8omlE1W4/5++W+paybFLCFuQjFLiOVcvjyR5dy6VW5kNlPKYqywScw+\neiRNhFasCP9aU3WzymKsw/aalpgdHZVRHToto4CMSVq2LDkLZpzFfBCSfP9PnIjffKgUSW/m6K6X\nVaRlNdZZr7xtm5nxPOfOSS3pzJnhX2uTmA0y39pFMWvjzHZC0oZilhDLuXIFWLNG/r9SxOz69cCN\nG8DTp6mEVJKzZ6XzcnV1+Neqxbzv64+rFLosxoDYOG/eTN52ef68ZKhmzNB/7CQF4dmzEneUxXwQ\nOjpkoZ3EzNak6mUVSYrZJEYKKdIUs7oy+hs2ANeupd8EKo7F3jYxW84V0tIi18Lh/Mk3FjE8LGU8\n69YBS5YAjx8DT56YjooQs1DMEmI5ymYMVI6YnTQJ2LTJ7LgJRVSLMQA0NEht2PXremMqh45OxgpV\nt5l0hjkJi7Fi167kxGyUesEwzJ0rM1uTyAglLWZ37EhuMyebleuEDit9PkrMJr0JpfP9nzpVNjXT\nbgIV5/zftEmur2lv9hUiSL1+TY3ca22p8y3E2bMiZCdPFsv8qlX2d2AmJGkoZgmxnEJi1obFQTGC\niFnAHqtxHDELpF836/t6OhnnkobVOEkxu3OnWFKTIKnmT7kk8f4nmdlUNDQAdXXAF1/oP7ayGOvu\nIA3I5kFNDdDbq//YipcvxaKrs17ZhNU4zvm/YIGcH0m+z0EYH5cNwCClArZbjZXFWMG6WUIoZgmx\nmpERWQioes6FC2WHvq/PbFzFePZMRo0sX17+uZUiZtOumx0cFKGyZIm+Y7ouZlevFrvdrVv6j51k\n8yfFzp36N0T6+mR+cHOz3uPm09GRzPmfpBD3vOStxqdOSSa1rk7fMdPuaPzggdhuW1ujH8MGq/HV\nq+KAmDev/HNdE7Nr1rBulhCKWUIs5vp1ES1Tpkw8ZrPV+OJFubkGqT+1QcyOjMhCYMOG6MdIOzOr\nLMY6M1aqbnNsTN8xc3n6VM7lOO9zKaqqZFNBtyC/f182Z9au1XvcfJLYTNDZfKgUSZ3/SWeVkxaz\nScSfdkfjo0dFQEfpJ6CwQcwGqZdVuCZmmZklhGKWEKvJtRgrbBazQS3GALByJfDwocx8NMXly0BT\nU7xur1u3AhcuiK0wDXRbjAGgvh5YtEh+jyQ4eVLq5yZNSub4QDJ1s6rOLs5iPgjr1olovndP3zGT\nFoOKJMTss2dy7ctdtOvGRTGbdhMoHfXiNojZMPOtbRazr17J9yJ3U5BilhCKWUKsppLFbFWVZBpM\nZmfjWowBEcJr1oitMA10djLOJckMc5IWY0USdbNp1MsCIpbb2/W+/2mJ2S1bxJHx/Lm+Yx47Jpsf\nuY4U3bS1yXc2iS7SgN5OxoopU+Rac+aM3uMWQ8f5b4OYVWN5grBsmThJHjxINqYonD8vJUe1tROP\nrVolGxzj4+biIsQ0FLOEWMzlyxNjeRSVImYB81ZjHWIWSK9u1veTE7NJ1s2GyYxEpb1dxMnIiL5j\nplEvq9D5/r94Id9Fnc2HilFbKzWVOq9JaQjx2bOlgdWlS/qP3d8v52GQ3gFhSatu1vflextXzLa0\niF3/4UM9cYVlZEQacQXN8nuenM9JuVTikG8xBoDp08VZY2sfDULSgGKWEIsplJldskSyCTdvmomp\nFO+qmE2rbnZgQP6rs/mTImkxm3RmduZMWTjrygKNj8sGRRqZWUDv+3/ypFiX49jnw6C7CVRaWeWk\nrMZJ1iu3taXT0fjaNRFKixfHO05VFbBxY3rZ5HyizLfesMFOq3EhMQvQakyIETHred5XPM+77Hne\nVc/zfrPAz5d6nvdDz/NOep532vO8/81EnISYppCY9Tw7s7P378tA9zBCK615j4Xwffcysyorm8Qi\nef16YGipwCtDAAAgAElEQVRIv73u9m2x7amO3Emis272yhXJeCxYoOd45dixQzZ2dNhe0xKDCp2b\nOWqk0K5deo5XiqTE7JEjyb3/aWVmdc5XNmk1juIKsbVulmKWkMKkLmY9z6sC8AcAfhpAK4Cve56X\nZ6TE/w3gz3zf3wrg6wD+MN0oSaWRzfZiz55uvP9+F/bs6UY2a3jwXQDu3ZOF7cKFb//MRjF74YIs\nAsIIrcZGqRc0YZHq75catELvb1hWrpTRMLdvxz9WKZKyGAPyOWzfrj/DrOrVku6qC+itm9W5mA/C\n3LnyfdCxiHZZzF67JuNsGhr0HK8USWdmk2DDBpnrm3QTKJ314ibFbJh6WYWNYnZ0VGLatOntn61e\n7c54HhfXYsR+TGRm2wFc832/1/f9UQB/CuCDvOeMA5j55f/PBjCYYnykwshme7F79wEcOvQRPvmk\nG4cOfYTduw9YfxFVWdlCIsBGMXvuXDiLMSC/mymrsa6sLCBWuvb25LOzaixPUnR0JNMROGmLsUKn\nVTet5k+56IhfZTbTFLPLl0s3b2WDj0OasW/aJNdZnc2rXr0Czp5N7ns6ZYqMikratlspmdkwY3kU\nSsyacAwV4/JlYOnSwnbpNWvcyMy6uhYj9mNCzDYC6M/5+8CXj+XSDeAXPM/rB/BXAPamFBupQDo7\nD6KnpxuAml5fh56ebnR2HjQYVXkKWYwVNorZsPWyiu3bkx2RUQydYhZIvm7W95MZy5PLzp36f4c0\nxezKlSJMdIiqNJs/KXSI2b4+qfdtbtYSUiA8T5/VPi2LMQBMnSrNfnR2Ij95Uq7bdXXlnxuVbduS\nrZt98UI6VG/Zoud469fL/Uxnc7YgDA9LFnvjxnCvW7BAnCq3biUTVxSKWYwBd2zGrq7FiP2YELOF\nzGb5+19fB/DHvu8vBfAzAD4udrBvf/vbf//nk08+0RclqRgGB8cxcfFU1GFoyO5e9qXE7PLlYmvV\nOZcyLlHFbHu7+5lZIPm62YEBEQ2N+Vt/GunoEPE5NqbneKojatKdjBWep0cQPnsmdled50cQdMSe\nZPOhUuhqAvXpp+lmlXVbjdPILCddN3vypIj83BEwcaitlXvWxYt6jheUkyflnjR5cvjX2mY1LiVm\nly6VbtFPn6YbU1hcXYsRM3zyySdvaLxS1KQT0hsMAFiW8/clAIbynvN/Qmpq4fv+Ec/zpnqeN8/3\n/beW7uV+QUIaG6sADOPNi+gwGhrsbuZ9+TLwi79Y+GeeJ7vmp04Bu3enG1chfD+6mG1rkxv1+LjY\nddPi9GngP/wHfcdTonxsTHb1daMsxkmKlHnzpIb44kWpzYtLT4/Y4hYtin+soChB+M/+WfRjnDgh\nFtQoi+A4rFsH3LkD3L0LzJ8f7RhpW4wVHR3Av//38Y7x5Alw/XrhusCkaG8Hvv99fcc7cgT42Z/V\nd7xCtLUB//k/J3f8JCz2ymqc5gZRlHpZhRKzNtxfAblHfpBfkPclVVXiSrl6NVnnTlxcXYsRM7z3\n3nt47733/v7v3d3dRZ9r4gw6BmCF53lNnudNBvA1AH+R95xeAD8FAJ7nrQUwpZCQJSQI+/Z9iJaW\nLshFFACG0dLShX37PjQWUxBKZWYBu6zGQ0Oy8I+yAK+vFxGVpk3q0SMRDDo77NbXi2hLYm4lkLzF\nWKGz7jRNi7FCR/wmLMaAbIK0t8ezen/2mZnYt2+X69HoaPRjHD0qm3RpbiK4mJldv17sszprfXNJ\novnZli3p181GqZdV2JSZHR+X966U7dsFq7GrazFiP6mLWd/3xwD8GoAfALgA4E9937/keV6353lf\n/fJpHwH4Zc/zTgM4BOCX0o6TVA6ZTBMOH96Lpqb9qKnpwubN+3H48F5kMk2mQyvK6Chw40ZpsWWT\nmI2alVWkXTd79qxkHnVnUJOsm02yk3EurovZ7dvl8335MvoxTDR/UsSpW37xQr6LSTYJK8asWUBT\nkzSCi0qa9bKK1aulXENHycbAgDSAWr48/rFKkXQTqCQzs2kSp8TBJjF77ZpsFM+ZU/w5LojZTKYJ\nf/3Xe1FTI2uxr37V/rUYcQMjuX3f97/v+/5q3/dX+r7/u18+1uX7/l99+f+XfN//Md/3N/u+v9X3\n/b8xESepHDKZJixe3IVf+qVurF/fZf3F8/p1GU0xdWrx51SSmE27bjYpu1tSdbO+n3wnY4VOQZ5m\nvaxi+nRZ2EX9bvi+ucwsEG8z4cQJsSpPm6Y3pqDErZtNu14WEIumroZKKiueRr1yW1sydbNDQ9I4\nSfdc6E2b5LqbVofg+/fFfVPK3VSK1lYptxi3oJyzVL2swpXxPFVVTVi0qAs/8RPd+NVftX8tRtyA\nRnXyznD9utScnD1rOpLyXLki7fZLsWqVdFt8/DidmEqhIzNbCWI2qcxsf78sutOYvbl+PTA4CDx4\nEO84o6OSOTJRwxVHEKqZx8uWlX5eUnR0iLB6/Tr8a03VyyrinP/j4yKETcSvy2qc5vufVEfjzz+X\nTTndgnz+fNlo6k1pCsvx4yIAo7pvZs2STGha8ZYiiJh1ZTyPckUtXy5rMkJ0QDFL3gmePZNOfz/x\nE9IkIU5dVxqUq5cF5Ca9caO5+X25xBWzW7fKMdIa3ZCUmN24Echm9XeVVBbjNDI+NTWS9YmbYT53\nDshkCs9FTJo4YlbVC6bdDVgxZ450rI5icTxyxF0xe+UKMHu2NCBLG1fFbBKZ2STqZRVpWo3j1Msq\nbLEaBxGzq1aJHdmGTHIpzp2T+2QmQzFL9EExS94JslmZu1hXJxkX23cwg4hZwA6r8fi42LFaW6Mf\no65OdmrTyJqPjIgdK474LsakSbJg051lTstirNBRN6tjMRmVXbvEshrF0miyXlYR5f33ffOZ2XXr\ngJs3xeIZlk8/Tb9eVqHEbBwL7KtXcv1Ky1a/fr10C9fdBCrJ8z9NMaujxMEGMev7co8vN/N3xgzZ\nDNIxYztJmJklSUAxS94JstmJphwbN8ZrUpIGly+XtxkDdojZbFa6Ec+cGe84adXNXr4sGxtJ1RUm\nUTebVidjhQ4xa6L5kyKTkRFJyjIcBpP1sooo739fn2wsNTcnElIgqqtl0yVKltOkEG9sFEdCHEvp\nqVOyAVmXP0YzIaZMkc0DneLw9WvZOEvqe5uWmPX9eGN5FDaI2Rs35JwK4lhwoW5WZWaXL5e1AyE6\noJgl7wTXr8sCF5ALqe11sy5lZuNajBVp1c0mPetQd92sav6UppjdsUMEydhY9GOYFLOeF00QjoxI\nna+JbsC5RIldiUFT9mhF1CZQJsWs58W3GpsYiaTbanzhggj7Ul1z45CWmB0cFGHeFLO3kA1iNojF\nWGF73eyLFyLOV6+esBmn1RCMVDYUs+Sd4Pr1iczshg12i9kHD8SytmhR+eeuWyc3h+Hhsk9NDIrZ\nN1GZWV036f5+yRql0fxJMX8+sGBB9Jm5T5/Kd27DBr1xhWHXrvCC8MwZYOVKaVRjkrVrgTt3pBtr\nUExbjBVRNnMePZLM8saNycQUBB1iNu33X3dH4yTrZQG5Bz94ADx8mNy/AUyUOMTd2Fm7VupQTfbY\nCCNmbR/Pc/GiXF8nT5YNk+rqaCUJhORDMUveCXLFrO02Y5WVDXIjnjRJBK1Jca5LzG7YIJ/Ts2fx\nj1WKpMXs0qXSeVhXF8w0mz/lEsdqfPKkjOKYNElvTGHYuVPqMMNgg8UYkEVee3s4UWiLmFWbOWEa\n0Rw5IsKspia5uMrhopjVnZlNul68qkruv0nNx1XoGglWWwssWQJ88UX8Y0WlksSsshgrWDdLdEEx\nS94Jcmtmm5tlZzjp3eGoBK2XVZi2GusSs5Mni6BN8nfxfRGzmzYl9294nt662bTrZRVxxKxJi7Gi\nrU1sky9eBH+NDc2fFGHe/xcv5Hc1cZ7ks3ChNKK5ejX4a2wQ4m1tcu2JMhJpYAB4+RJoadEfVylU\nEyhdzpykM7NAOlZjHfWyCpNWY1ViEkbM2lwze/bs22KWdbNEBxSzpOLxfblgqprZqio7amGKEbRe\nVmFSzI6MyK51GPFdiqStxv39wNSpyY//0Fk3m3YnY4XrYra2Vjpsh5nFaUtmFgj3/p84IZbIpJqa\nhSXs+W+DmJ09W7JwFy+Gf60aiZS2e2LyZDnHdWQ6Hz8WN0nSpQFJi9nxcfnOV4KYHRqS/zY2Bnv+\nsmVi2zVZdlSKc+fePL84nofogmKWVDy3bkkNXG4dnM1NoMKK2S1bzInZa9fkBlpbq+d4SYvZpC3G\nCl2ZWRPNnxTr10vG6cGD8K+1QcwC4epm796VhWCY716SdHTIojxIptAGMZhLmCZQY2NyvtgQf1Sr\nscn3X5fV+Ngx2RhN2uqdtJj94gvZmFiwQM/xTIpZZTEOuklSXQ2sWBHOFZEmhTKzFLNEBxSzpOLJ\ntRgrbG4CFdZmvGGDCOBXr5KLqRi6LMaK7dvj1a2VIy0x29Ym2ZKRkXjH6etLv/mToqZGfo+wovzW\nLWkAlbblshBh6mY//1w2IaosuSvOmSOZwiALaRvFbNDM7MWLIjzmzUs2piBE3Uwz0clYsW1bOPdB\nMdKy2Le2itiKe20shq56WYUNYjYMtlqNb9+WRlq59zLajIkuLLltE5Icuc2fFLY2gXr9Wi7uK1YE\nf01trTzfxA1Xt5hdvVqyY/fu6TtmLmmJ2enT5TOJa/8zlZVV7NwZ3i6t6tVMj4gBJqy6QTpL21Qv\nqwhiNfZ9+8Ts5s0iWILYHW2KPUpm9tUr+Z7rFFBh0JWZTaNeFpD71fLl0ezcQdBZLwtI992+vnC1\n97qIImZtHc+jmj/l3heYmSW6oJglFU/ujFnFhg0ixMJ03EyDbBZYvDi8bddU3axuMVtVpS/TUIi0\nxCygp27WVL2sIkrdrBqLYQPLlkmGOcjuv031soog739fnwja5uZUQgrElClyjQ3yPf7sM7GD28Cm\nTSIEnj8P/ppTp4BVq8yNc1q/Xu5xceokfT/dzZwkrca6xeykSSJoTWQ7o2ZmbRSz+RZjQK7Pg4Nm\nRx+RyoBillQ8hTKzc+YAs2bpG5+ii7D1sopKEbNAcnWzjx5JXWRa9lcddbOmOhkrOjokUxVm08eW\nellAsgBBrMZjY3LOuZiZVZlNGzLhuQStm/30U3sys1Onig321Kngr1HNn0yhmkDFEYfZrGxALFmi\nL65SJCVmR0clS677mmnCanznjpRr5G/El8NWMZvf/AmQc3fRImnMSEgcKGZJxVOoZhawswlU2HpZ\nhQkx+/y5NAgKY4kOQlJ1s2fPys20ulr/sQsRNzNrsvmTYv58+RPUEuj7+mvW4hJEEF6+LB2u6+vT\niSko69bJBszdu8WfY5NNN5cg5//9+1Jj3dqaTkxBCGs1tuH9j2s1Tttin5SYPX8eaGoCZs7Ue1wT\nYvbUqXDNnxRKzNrmOiuUmQVYN0v0QDFLKp5CmVnAziZQUTOzmzbJzTZNu86lS2KvmzRJ73Hb2yVL\nFqTOMQxpWowB2ZS4dy96/W9fn+xcm2j+lEtHR3CrcU8PMGOG7LbbQhAxa6PFGBDbfXt7aVFog5gq\nhDpvSn2PjxyRjY+0NpiCEEXMmj532triidm06mUVmzbJ9Vj3NV63xVhhQsxGsRgDIuRnzRL7ri28\nfi3rhUKbVqybJTqgmCUVzcuXYtcpZJ+ysQlUVDE7YwawdGm6dT1JWIwB+T18X7K+OklbzFZVycIq\nqtXYtMVYEaZu1iaLsWLrVvlePXtW/Dk2Nn9SlHr/X7wALlwwW1ddjKYmyQ6VshDaVC+rCCNmBwfl\nM9DtTglL3D4DaZ//8+dLjbHuMh+KWcE2q/EXX8imbKG6cs6aJTqgmCUVTW+viKNCO/822oyjilkg\nfatxUmLW85Kpm01bzALx6mZNW4wVYToa2yhmp06V73qp88nWzCxQWsyeOCFWZF1znnXieeWtxjbV\nyypWrw7uqFBZWdP1yq2tYtWM0gTq1Su5lqd9rUnCapxU87nmZrHEP3mi/9jFqCQxW8xiDDAzS/RA\nMUsqmmL1soBc8Ht7zbTcL8TDh1KHGtVWWiliFtBfNzsyIlnrpOItRpy6WdOdjBUbNkh27eHD8s+1\nUcwCpQXhkydynSi22DLNjh2SdXv9+u2f2WoxVpRqAvX6tfxetm0iVFXJ9y7IZprp5k+KyZPl2hZF\nHKpuzHV1+uMqhW4x+/y5jIPatEnfMRVVVbJpdOGC/mMX4uFDcZStXBnt9bbNmi3U/EnBmlmiA4pZ\nUtEUGsujmDRJbuJJzbsLi8rKRt3lryQxq+pmdXH5suyuT5um75hB2LEjfDdgQGzWttiMa2pkcV8u\nw5xUJ1EdlBKzx48DW7bor/3WxZw5UiZRqCTChnrNUpTazDl/HmhslN/PNoJajW3aTIhqNTZlsdct\nZk+dEsE5ZYq+Y+aSptX49GkR5VFryW2bNVsqM0ubMdEBxSwJTDbbiz17uvH++13Ys6cb2axlc20K\nUKz5k8KmJlBxLMaALMhPn06ni+GjR7J73NSUzPG3b5fMpK7fxYTFGJDasPr68AuL3l5ZlC1enExc\nYQlSN3vunCxMTM3bLIWKv1DDGZstxopC77/v2yWmCtHWJt+9kZG3f2ZjvawiiDNkZEQ2b2zp3B21\no3HazZ8UusVsUvWyijTFbByLMWCnzbhYZnbBAnHHpWnhjoqLa+B3BYpZEohsthe7dx/AoUMf4ZNP\nunHo0EfYvfuA9V/mUjZjwK4mUFHH8ijmzBHxdO2avpiKceGC1GlVJXQFmTdPfh9dv4spMQtI1iOs\n1dgWi7EiSEdjWy3GgGQ2p00rfD7Z3PxJUUjM9vaKoG1uNhJSIGbMkLnOZ868/TMb62UVQTqqnzol\njZ9s2byJ2tHY1Pm/fDnw4EGw8oUgJFUvq3BJzDY1iU05Sg21bp48kViKzXf3PDesxq6ugd8VKGZJ\nIDo7D6KnpxuAKqypQ09PNzo7DxqMqjzlMrM2NYGKm5kF0rMaJ2kxVuismzUpZkvVDRbDluZPCvU7\nlMqU2yxmgeLZTVcys/kbIiora7r5UDmKWY1tzio3Noq9vlS3Xdvib20Fbtwo3bU7n9u3RUzGve9E\noapK7r+FNjqikPR8a5fEbHW1bLSksbFdjvPn5dwsZZl2oQmUq2vgdwWKWRKIwcFxTHyJFXUYGrJs\nMncOvl+6ZhaoLJsxUFliVlfdrO9P1CCZIEpm1pZ6WcWCBZL1v3Sp+HOSzozEZdeut8XsjRvSPKfQ\n6C6bWLcOuHtX/ihsaT5UjkKbOXfuSHfYtWvNxFQOzytfN2ubmJ00SURDGOvu55/L75mUw6YcuqzG\nDx8Ct24lez4tXix9Ae7cSe7fAGQzorc3/u9ii9W4lMVY4ULdrItr4HcJilkSiMbGKgD5npVhNDTY\newo9eCCLklINRhoagLEx2aE2yevXcjGP2r1QUUliVtd4nv5+Gc+ycGH8Y0Vh82bZIQ9q+fJ9+zKz\nQOm62adP5fwtt2gxyc6dYm3NxYWsLCBio739zU0R28RUMQplZj/7TDZ5TImoIJQTszZuJoS1Gpuq\nl1XoErOqiVvUhklB8Dy55yXd0fjMGfl34jaks0XMnjtXvlO8C5lZF9fA7xL8FEgg9u37EC0tXZj4\nMg+jpaUL+/Z9aCymcqh62VI2PM+zIzt744aIrbjddrdsETFbqtYrLr4vN6ikxezWrfK5jI7GO45J\nizEgjZw2bAjeabS3V8S3Lc2fFKXE7MmTkvm2tSMwIOdAT8+bjUZcqJdV5L7/L17IotqmuupirFnz\ndlbZBSFeSswODcnm1IoV6cZUjrBNoEyf/7rEbFqukDSsxnEtxgpbxvMEycy6UDPr4hr4XYJilgQi\nk2nC4cN7sWLFflRVdWHz5v04fHgvMpmE2tlqoFy9rMKGJlA6LMaACOK6OhHHSXHnjgjaRYuS+zcA\naR7T3Bz/szEtZoFwdbO2WYwVpZpA2V4vC4ideMuWNwWKK5lZ4E0xe+KEWI9ra83GFASVVc49/10Q\ns21tIiyKzfft6LCvXjnMeJ6xMRGBJsVsa6vMhi3U7ToMSdfLKlwSszaM51Eb30HErO2Z2UymCT/4\nwV7U1Mga+Od+zv418LsExSwJTCbThDVruvCVr3Rj69Yu67/E5eplFTY0gdIlZoHkrcbKYpzGQk5H\n3awNYjZM3axtnYwVGzeKZbtQ91EXxCzwZt3sy5dyLtu4cVCIHTtEqLx+7YYYzCXXajw6Ktcn2zPi\ns2dLLXWhOeS2vv+treLsCNIE6tIlqYWfNy/5uIpRWytCJu6s96TH8ihcErOrV8tGQZIurXL094vb\nbP780s9rbpYN+DTGCsZhypQm1Nd3YfPmbvzWb9m/Bn6XoJgloRgcBH7yJ+2wr5QjaGbWBptx3LE8\nuaQlZtNAR92sDWJWLeaDLCxsrJcFpLvrtm2FrZeuiNncutnTp+U7F9fanxZz5gBLl0qmw1YxVYxc\nZ8KZM7J4nTXLaEiBKGY1tvX9nzRJrs1BrLum62UVca3GQ0PAq1fBNq7j0toq97+kBOLLl9JfQcf9\nddYscWkNDsY/VlSCWIwBuQbPmSOfpc2odVpjo9n3lbwNxSwJhRKzly6Z3fELQrkZs4rWVrlIFbKT\npYWLmdk0iCtmHz2SWr1iM+7SorlZzq+BgdLP8317bcZA4brZW7ekAZTp9zgIasTN+LhbFmOFev9t\nFVPF2LFDROHYmFuxFxKzIyMiyG3dvAlqNTZdL6uIK2ZVVjYNp1B9vQjE/v5kjn/uHLBqlfRM0IHp\nJlBBmj8pXKibpZi1F4pZEpiREbEYKiFz757ZeMoRNDM7fbp0NTY5k023mD1xIrnNhjTF7MaN4ToB\n56N2hpPschkEzwtWN3vjhljvkq5HjkohMauar9hWP1iIRYvEPnrlij2L+TC0tPTi29/uxv37XfjW\nt7qRzZYYhGoR8+aJ1fDyZffF7OnT0vhp+nQzMZUjaBOoSsnMplUvq0jSaqzLYqwwXTcbNDMLuFE3\nSzFrLxSzJDBDQ7IYrK6WL7TNVmOVBVu2LNjzTTaBevRIapwaG/UcTx0nCcuO70sX1dZW/ccuxJQp\nsng4dSra622wGCuC1M3aWi+rUII8t7Yp7cVkXFTdrGuZ2Wy2F//lvxzA3bsfYXS0G9/5zkfYvfuA\nM4JWWe0/+0w+AxfYtEnqDp8/n3jMdjEeZDzP06fS2dvU7O1cNm2S63TUzde06mUVLolZlzKzLsya\nVWJ2yZLyLiuSLhSzJDCDgxNCyXYx298vnX2nTAn2fJNNoFRWVldmy/OSsxr39UmX4blz9R+7GHGs\nxjaJ2SCZWZstxoA0jKmvlzIDhSv1soqdO4H/8T9kQR93rnOadHYexMBAN4C6Lx+pQ09PNzo7DxqM\nKjgdHcD//J8yGmnVKtPRBGPKFNm4y91MU52MbWXdOmkC9fRp8eccPy4icvLk9OIqxvz5kuXujbAn\no8oyKGYLY3I8z6tXsmGydm2w5zMzS+JAMUsCky9mTbd9L0XQelmFySZQOi3GiqTEbJoWY0WliNnt\n22VRXGpurq3Nn3LJtRr7fvqZkbgsW9aLv/qrbgBd+IVfcMeqOzg4jgkhq6jD0JDlLUC/ZOlSN9/3\nfKux7ZnZIE2gbLPYR7Ua9/SIEE6zLCMpMTs6Kq4nndlyk+u0S5dkDRY0oWB7zeyzZ8CDB+L2o5i1\nD4pZEhiXMrNBx/IoTNqMKWZLs3174Y6i5RgZkXM07XiLMXNm6bm5vu+OmFV26Z4eydTbWuObTzbb\ni29+8wCAj/DgQTcOHXLHqtvYWAUgv3h8GA0N9t/Gs9le/Pqvu/m+54rZoSGp37c9o1/OamxLvawi\nqpg1UeKwbp3cV8bG9B734kW5P9Tl71fFoLkZuH0bePFC3zGDEsZiDNhvM756Vb73VVUTYtb2Jqjv\nEvbfBYk1DAxIrQBg1r4ShKDNnxTLlwN37ogFLm0oZkuzZo18Ng8ehHvd5ctyM7dp9EqputlsVmK1\nXRjmZmZdsxh3dh7E9etuWnX37fsQLS1dmBC0w2hp6cK+fR8aiykoLr/vuWJW1Vnb3uysVEdj36+c\nzKwJV4jKBPf06D2ubosxIP1Nli8309wyTPMnQJpwPnhgRngH4fLliXXazJlyDXj82GxMZAKKWRKY\n3Mzs8uUibl++NBtTMcLajKurJ2bIpY3OGbOKTEZqpu7c0XtcE2K2ulpu8kHGTeRik8VYUapu1oWs\nLCC77b290rjMNTHrslU3k2nC4cN78Y1v7Mf773fhG9/Yj8OH9yKTaTIdWllcft9Xr5bO/ffu2W8x\nVpTqaNzXJ/9tsui0iSNmTVx/krAaJyFmAXOJh7CZ2epqOSdv3EgspFjkrtM8j1Zj26CYJYHJFbOT\nJolg+uILszEVI2xmFjDTBGpsTHZ4ddvWVBOoqF2AC/H6tWSR163Td8ygRKmbtVHMlsrM2t7JWFFT\nI4vlzz93T8y6bNUFRNB+/HEXfvjDbnz8cZcTQhZw+32vqpLv5bFj9jd/UqxbJ6K1UBMolZW1Kbu8\nfLlk5R4+DP6a16/lGm9iA9AlMWuqbvbs2XBiFrC7CVR+0oFi1i7sv5MQa8gVs4DddbNha2YBM02g\nenulm6POOhmFbqtxTw+weHEysZYjSt2sjWK2tRW4ebOwZdqVzCwg2am//VvgzBl3Ygbctuq6jOvv\ne3s78KMfyTXFhc2bSZPkflYo22lbvSwgGwYbN8r1JCgXLkjZ06xZycVVDN1idmxMfvck7lcmxvPc\nuyfjrJYuDfc6m+tmKWbthmKWBML3pfmFC2L26VO5kC5cGO51JppAJWExVugWsyYsxor29nCZWd+X\nhZwNcxRzqa4W8ZcvzF1p/qTIZHrx+78vnWn/xb9wpzOty1Zdl3H9fW9qmjjf/+W/dON8L1Y3a1u9\nrCKs1dhkF3XdYvbqVanDnT1b3zEVJmzG587JZkrY7L+tmdmxMak7zh0ptmQJxaxN1JgOgLjBvXuS\nkR3eBVoAACAASURBVKutnXhszRrg8GFzMRUjm5UdvrAX0g0b5CLs++lZsJJo/qTYuhXo7NR3vPPn\nwzV00ElTk4wuyHcHFKOvD5g6NfyGRhqoutmvfGXisWxWvl82xptPNtuL3/3dA3j5Uhr6HDo0jCNH\nupwRJ8qqS9LF1ffd1fO9rQ34m79587GREckA2jhKa/Nm4O/+LvjzTdXLAnLPzmZllmrQ0TOlSMpi\nDEisV6+mu66JYjEGRMz+6Ef644lLb6/MWM91pTU2mpuAQd6GmVkSiEIiwtbMbJR6WQCYN0+6yfb3\n64+pGEmK2ZUrpQFUmDqkUpjMzHpeuLpZGy3GikJ1sy5lZTs7D+LGDTc70xISls7Og+jtde98L9QE\n6swZuTfOmGEmplK4lJmdMkU2zHXZd0+eBLZs0XOsfObMkSTEzZvJHL8QKjMbFlszs4UcdLQZ2wXF\nLAlEITGrajFsm7UVpV5WkXYTqCTFbHW12GyjdIkshEkxC4Srm7VdzB49+ub35vhxd8Ssy51pCQmL\nq+d7oSZQNtbLKlpbJYP46lX55754IQLD5DVep9U4ycwskH7dbNTMrKqZtW1NSTFrPxSzJBCFxOzs\n2TJzbWjITEzFiJqZBdJvApVkzSygr2725UtpmZ9bM5I2YepmbRazixfL9yZ39p8rnYwBtzvTEhIW\nV8/3mhoRFLkd7W2tlwUke9jSAly8WP65p0/LfXPq1OTjKoYuMTs+Lp9RUplZIN262bExac4VZeN7\n1izJet+9qz+uOFy5QjFrO3ZfjYk1DAxIwXs+NlqNw86YzSXNJlBPnsiueZAa0KjoErNXrshCY/Lk\n+MeKyvbtksEMsmtrs5gF3pw361rzJ9c70xISBpfP93yrsc2ZWSC41dhkvaxCl5jNZsX2vWBB/GMV\nI83xPNevy+8yc2a01y9fLu+JTRRKOixcKCVcQZwEJHkoZkkgijXeMTWQuxRxMrNp2oyvXJG61qoE\nv4W6xKxpizEwcYMsN9v40SPZ2W1pSSeuKOTWzV6/LpnaJBczOnG9My0hYXD5fM/taHzvHnD7NrB2\nrdmYShFGzJpuYqVLzCZtMQbStRlHtRgrbKybvXz57XKw6moRtGnWIpPisJsxCUQxMWtbZnZ8XOyw\nzc3RXr9mjVxIdXUpLEUh64pu1q6VTnzPnolgiooNYhaYqJtdubL4c86ckZtpdXV6cYWlowP4znfk\n/12yGCtc7UxLSBRcPd/b2oDf+z35/6NH5fpp83Vx82bgL/+y/POOHgV+4zeSj6cULS3ArVvx762V\nJmajNn9S2DZr9sEDqdFevPjtnymrcdT1JtEHM7MkEK6I2Vu3pO6iLr9fR0CmTJGb1KVLeuMqRKHd\nPt1MmiQiNMww+kLYImaD1M3abjEGpD7q4kW5SbpkMSaEuMPatVIi9OSJ3fWyik2b5F5VqpTk0SNZ\nj6xbl15chaiulvVPkBrfUqQhZjMZ6W3y4kWy/w6gJzNrk81YJR0KjTVi3aw9UMySQLhSMxvHYqxI\nqwlUkp2Mc9FhNbZFzAYZz+OCmK2tle6dJ0+61cmYEOIONTVyPzt1yv56WQCYP1+ynDduFH/OiROy\nGVhjga8wrtXY99MRszU1si4qV6Kjg7Nn42VmbbMZl2rSuWSJrI2JeShmSVmeP5dutnPnvv2zZcuA\n+/ffbP9vkjhjeRRpNYFyRcw+fSrzauO+rzrYtk127kdHiz/HBTELSJbks8/ks6GYJYQkgaqbPXrU\n/swsUL5u1oZ6WUVcMTswID0zGhr0xVSMNKzGz55JBrhUGVA5bLMZlxKzzMzaA8UsKcvgoFxsC9ks\nqqpkXMvVq+nHVQgdmdk0mkCNjcloljRG3cQVsxcvil3NhlqrmTOBpUul9X8hRkbk5mNDFrkcy5f3\n4nd+pxvPn3fh13+9G9lsr+mQCCEVxrJlvfjt3+7Gixdd+Df/xv7rTDkxq2p/bSCumFVZ2UJrK92k\n0azzwgVZK8TJmi9bJk2VSm1YpwnFrBtQzJKyFKuXVdhkNY4zlkexYUPymdm+PmDevHiNI4Kyfr0I\n55cvo73eFouxolTd7KVLsrM7bVq6MYUlm+3Ff/pPB/Dw4UcYGenGoUMfYffuA9YvNAkh7pDN9uLA\nAbnOvHrlxnXmXcrMpmExVqSRmY3b/AmQPh8NDbJGsgGKWTegmCVlKVYvq7BJzOrIzC5dKtbqJAd3\np2UxBmSw/KpV0QW6bWK2VN2sKxbjzs6D6OvrBqA6ldWhp6cbnZ0HDUZFCKkkOjsPYmDAretMKTF7\n6xYwPGzP2LWlS8Vae/9+tNenKWbTmDUbt/mTwpa62ZERmQZR7HxbsoRi1hYoZklZgmRm02r7Xg4d\nNbOel3x2Nk0xC8SzGlPM6mdwcBwTC0xFHYaGxk2EQwipQFy8zixfDjx8KCNR8lFZ2TRsuUHwPLk3\nFit7KYeJzGypTtFxidv8SWFL3ez162J7LjamUWVmk3xPSTAoZklZyonZNGoxgvDiheyQloo1KEk3\ngSplXUmCShKzmzbJTfn587d/5oqYbWysAjCc9+gwGhp4SSaE6MHF60xVldx/C42Ts6leVhHVanzr\nltzD0ppROncuMHmy/LtJ4PuyZqqkzGy58Ym1tTIG8t699GIihbH3ilahZLO92LOnG++/34U9e+xv\nxgCUF7OrVklN5thYejEVordXdtF0NCpKugmUK5nZe/dkk0DHBoEupk6VGYP5VjTfl8c2bTITVxj2\n7fsQLS1dmFhoDqOlpQv79n1oLCZCSGXh6nWmmNXYpnpZRVQxe+pUes2fFEm66IaGpPHTwoXxj2XL\nrNkgSQdX6mZd1B5hsGBS17tDNtuL3bsPoKdH1bAM48iRLhw+vBeZTJPp8IpSTszW1QELFoiYjFuv\nGgcd9bKKDRuAP/ojPccqRNpidtMmsUKNjkqDhaBcuCA3a1tsXQplNd61a+Kxvj7ZKdVxM02aTKYJ\nhw/vRWfnfgwNjaOhoQr79tl9HSCEuIWr15nNm4G/+7s3H/N9ueb/8R+biakY69cDf/7n4V+XpsVY\noazG772n/9g6mj8pbMrM/viPl36OErM2O8Jc1R5hoJhNkc7OgzknEzDRjGE/Pv64y2RoJSnXAAqY\naAJlWszqmoWq6mDGxvSPpHnyBHj0SJpHpEVdndiZLl4Ml7m0zWKs2L4d+OEP33zMFYuxIpNpsvp7\nTwhxHxevM5s3AwcOvPlYNiublYsXm4mpGCoz6/vhNn1PngT+6T9NLq5CJFkSpqv5E2BPzezly8Av\n/3Lp5yxZImtkm3FVe4SBNuMUcbEZw9gYcOdO+RuIDR2NdYzlUcycKdnmnh49x8vl6lUZKl6V8rcv\nitXYZjF79Oibj7kmZgkhhLxNa6vcJ1+9mnjMxnpZQNYJNTUyGzUMJjOzSaCr+RMgYwtHR2XT3xS+\nXzk2Yxe1R1goZlPExWYMt29L44By1lQbxKxOmzGQXBOotC3Giihi9tw5O8XsunWyeMi92VHMEkKI\n+9TWyjiUixcnHrOxXlYRtm72wQPpR7FyZXIxFSLJmlldzZ8AyXCbrpu9fVvWvfX1pZ/ngph1UXuE\npXJ+EwdwsRlDuXpZRaWK2SSaQLkiZn1fbtCtrcnFFJXqamDLFuD48YnHKGYJIaQyyG8CdewY0N5u\nLp5ShBWzp07J75e2OyuTkTXdy5d6jzs6Kpn0dev0HdN03WzQiRMuiFkXtUdYKGZTRDVjaG/fD6AL\nP/VT+60vwA5SLwuYF7O+r7dmFkhu1mzaY3kUmzfLuIOgXacHB6Vz8Pz5ycYVldx5s48eAXfvFh9u\nTgghxB1yxezr1yIA29rMxlSMsGLWhMUYkExjczPwxRd6j3vlCtDUJBl1XZiumw26TluyxH4xq7TH\nwoX7MWlSF77xDfu1R1goZlMmk2nCP/7HXQC6sWdPl/UnU9DM7KJFUt9SaNB5Gty/Lxfq2bP1HbPS\nMrOzZ8vndPVqsOfbWi+ryK2bPXNGPi/dzboIIYSkT66YvXRJ+nbovL/rxBUxCyRTN6uz+ZPCtM04\nTGbW9gZQgGiPujrRHn/yJ/Zrj7BQzBrg3j1g2jTgxg3TkZQnqJj1vGTrMcqh22IMACtWyOy0Z8/0\nHXN8XGbyrlql75hhCGM1tl3MtrdPZGZPnxbbMSGEEPdRTiI1ksfWellASnEuXpT7exBMitkk1mk6\nmz8pTNuMr1wJJmbr64EXL4Dnz5OPKQ5jYyK6q6uBp09NR6MfilkD3L8vC28bhkKXI6iYBZJt+14O\n3RZjQDoUrl0rI3p00d8vDbVmzNB3zDBUkpjNZKT25+ZN1ssSQkglMW+e3Cdv3LC7XhYAZs2S+3qQ\nBMWTJ7IOWLs28bAKkkRmVmfzJ4UNNuMgDjrPAxoa7LcaDw7Kd2rRItEglQbFrAHu3ZPaDxcys0Fr\nZgGzdbNJZGYB/VbjoBfIpKgkMet58j06doxilhBCKg1lNbY9MwsEtxqfOSNZzJqa5GMqRBJJhyRs\nxs3NQF9f8B4fOnn+HLh1S2IIggt1szduyO9TXy8apNKgmDXA/fvAtm2Vl5k1KWZ1zpjNRXcTKFP1\nsootW6SRRjk71NiY1Cnp7E6YBNu3Az/6kbyvNgtvQggh4di8GThyRCy8tm9WBhWzJi3GwITN2Pf1\nHO/hQ+DxY2kApZPaWhFeQ0N6jxuEa9ekmWTQDQcX6mazWcl219czM0s0cf++XJhv3ZKW5rbi++6I\nWVcys6bF7Pz5Yt0qt5GSzcpzZ85MJ66oLFvWiz/4g254Xhd++Ze7kc32mg6JEEKIBhYv7sWBA92o\nqurCr/yK3dd3V8Rsfb2ItDt39BxPzaJPYsyQqbrZsBMnXBjPozKz8+ZRzBJN3LsnvvXFi6V2wlae\nPBErZ1BB09IC9PYCIyPJxlWIJGpmAcnMnj2rbxfT1FieXIJYjW23GANANtuL3/mdA3jx4iM8f96N\nQ4c+wu7dB6xe8BBCCClPNtuL3/s9ub4PD9t/fXdFzAJ6rcZJNH9SmKqbrUQxm5uZpc2YxGZsTCwZ\nc+bILonNdbODg8HrZQFgyhRg2TKgpye5mAoxOipWlGXL9B974UIZ+aPL6mI6MwtUjpjt7DyIGze6\nAdR9+Ugdenq60dl50GBUhBBC4tLZeRB9fe5c39euFXtqKbfd8+cy49X0vVVnE6gkmj8pTI3nCStm\nXaqZZWaWaOHhQ8l01tTILonNdbMDA8EtxgoTVuP+fslyT56czPF1WY2fPZM5vEmI7jBUipgdHBzH\nxEJHUYehoYDzEQghhFiJa9f32lpg6VIRtMU4d07WSFOmpBdXIXSO50mi+ZOCNmN9MDNLtHL/vpxM\ngBuZWRfEbFL1sgpdTaCuXgVWrkymtiQMSsyWsk67IGYbG6sADOc9OoyGBl7WCCHEZVy8vpezGttg\nMQb0ZWbHx5NdK5gQs+PjslYL46CzvQHU6Kj06Fm6lA2giCbu35c0P2B/ZjaqmNU9w6wcSdXLKnRl\nZk2P5VE0NIigLnbxHRkRq7jp2t5y7Nv3IVpaujCx4BlGS0sX9u370FhMhBBC4uPi9d0lMasj6XDj\nhpTMzZkT/1iFMFEzOzAgv8+MGcFfs3ixNNQyMUYoCP390qdn0iTajIkm7t2r7MxsEjPMypHUWB6F\nLjFrQ70sIE29SlmNr16VNvtTp6YbV1gymSYcPrwX3/jGfrz/fhe+8Y39OHx4LzIZzTMCCCGEpIqL\n13dXxOzy5SLaXr2Kd5wkLcaAiMTHj6XWOC2iJB0mTxYBfPt2MjHF5caNiYRPpdqMDY1tfndxKTM7\nMAB85SvhXqNsxr4voikNrl8Hfu7nkju+auwwMhKvLvfKFeCrX9UXVxyUmP3gg7d/5oLFWJHJNOHj\nj7tMh0EIIUQzrl3fS4nZkRGZ3Z6k+AvK5MmyYd3TE2+W/LlzyXUyBsRB1tws6+TW1uT+nVyiTpxQ\nTaAaGvTHFJdsVt5HgJlZoonczGxDg/w97u5YUkTJzNbXi5UhzR2qpGtma2vlQhDXPn3lij3W3VKZ\nWZfELCGEEGIDK1eKpfPFi7d/duGCJDDq8ntaGUJH3WzSmVkg/brZqGLW5rrZdyEzSzGbMrkNoKqr\npSC7186xaZHELJB+E6ika2aB+E2gVFOBVav0xRQHillCCCFEH5MmiaC9dOntn9liMVboKAlLOjML\npF83G0fM2trRODczO22a/DdN63YaUMymTK7NGLC3bnZkRMYILVgQ/rVpitnHjyWzPX9+sv9O3LrZ\ngQFg1iwZy2QDzc1yMbt16+2fUcwSQggh4SlmNbZRzMbJzD5/LomYpPuApD1rthLFbG5mFqhMqzHF\nbMrk2owBe+tmh4ak+1l1dfjXpilmVfOnpOtz44pZW5o/KVQTqFOn3nx8eFguyCtWmImLEEIIcRVX\nxGzcyRMXL8qaZtIkfTEVIk2b8ePHwJMn0RyJqmbWRnIzs0BlWo0pZlPGlcxsVIsxkK6YTcNiDMS3\nGUfd7UuSQlbjS5fkBlXD1nCEEEJIKAqJ2devZTN882YzMRVC2YxLzZsvRRoWYyBdm7FKOlRFUEa2\n1sy+eiXCNXc9z8wsiY0rmVmXxGySzZ8UTU2ya/fwYbTX25aZBQqLWVqMCSGEkGgUErNXrkjDz1mz\nzMRUiHnzxKF1926016fR/AmYWCNHFd1hiJN0sNVm3NsrWeNclyUzsyQ2uQ2ggMrMzDY3SzfjNArM\nk54xq6iqkptU1OwsxSwhhBBS2TQ3Aw8eyOa3wjaLMSBCNo7VOK3M7MyZ0rTozp3k/y0dYjYN0R2G\n/HpZQDQIM7MkMr4vFzlXMrNLlkR7bXW11Fxeu6Y3pkKklZkF4tXN2jSWR7FypezO5WabKWYJIYSQ\naFRVyezWCxcmHjt1yj4xC0RvAuX7wJkz6c3MTatuNs46beZM2SB48kRvTHHJr5cFaDMmMXn8WGaW\nTp488diiRXLy29Yme2AgemYW0NP2PQhp1cwC0cXs8LBYeZYt0x9THKqqpIYntwkUxSwhhBASnXyr\nsY2ZWSD6Ou32bRG0ixfrj6kQadXNXr4c3UHneXZajYtlZmkzJpHJb/4EyBegqck+q3EcmzGQTt3s\n+LjUA+TvOiVF1CZQV69KpjpKZ+ikybUaP3woGy62iW5CCCHEFXLF7Pi4bBhv2WI2pkJEzcwqi3HS\nUyQUaYznef1aBPPKldGPYWMTKGZmiXby62UVmQzFbBSGhoC5cyeGQCfNhg1ygxofD/c6G+tlFbli\n9sIFoLU1Wic/QgghhLwpZq9fB2bPfjuRYQNRa2bTav6kSMNmnM1Kprm2NvoxXMrMUsySyOR3MlY0\nN9tVN+v7IhRtF7Np1ssCckOaMyf8xoONY3kUuWKWFmNCCCEkHrli1laLMQC0tAB9fcDISLjXpdX8\nSZGGmNWxTrNRzBbKzNJmTGJRyGYM2JeZvXcPqKuLt0O1erXYa8NmMcOQZr2sIkrdrM2Z2bVrgf5+\n4OlTillCCCEkLosXi231zh27xezkycDSpUBPT7jXpZ2ZTaNmVoeYXbLELjH7/Ln05Fm06M3HaTMm\nsXAlMxvXYgwAM2ZIFrO/X09MhUhrLE8ulSZma2pEwJ45QzFLCCGExMXzJrKzNotZIHzd7OvXIvxa\nW5OLKZ+lS6XpVNgMchh0ZWZtqpm9cUN6oOSXjjEzS2LhSmZWh5gFkrcap20zBsI3gRoflwy1rWIW\nkBvtiRMUs4QQQogO1Fx628Vs2LrZa9dkfVhXl1xM+dTUSNaztze5f6MSbcaF6mUBSTa9epXs5kDa\nUMymSLEGUJWYmQUqU8yGzcwODsqFY9as5GKKy9atwP/6X/L/CxeajYUQQghxnfXr5b46aVJ6I2yi\nEHY8T9oWY0WSdbO+D1y6VHlitlC9LCDOgUprAmVEzHqe9xXP8y57nnfV87zfLPDz3/c875TneSc9\nz7vied4DE3HqppjNeN482SWxZdjywIDsgsUlDTGbds3sqlVinQ46F9hmi7Fi4cJe/PVfd+P16y78\nwi90I5tNcPuTEEIIqXDmzpX76shIF/bssfe+GtZmnHbzJ0WSdbPKcjt/frzjLFwoIw5tyXgWy8wC\nlWc1Tl3Mep5XBeAPAPw0gFYAX/c87439EN/3f933/S2+728FcADA99KOMwmK2Yw9T3ZPbLEa68rM\nRp1hFoTnz+Wi0dCQzPGLMWmSCNqLF4M933Yxm8324pvfPADgIzx+3I1Dhz7C7t0HrL3xEkIIITaT\nzfbit35L7qsPHth9Xw1rMzaZmU3KwajWaXHn5lZXi6C9eVNPXHEplpkFKq8JlInMbDuAa77v9/q+\nPwrgTwF8UOL5Xwfw3VQiS5himVlAdk9ssRq7YDO+cUO+pCZmooaxGl+5Yu9YHgDo7DyI69e7AagC\nmDr09HSjs/OgwagIIYQQN+nsPIgbN9y4r86fL709gmbpKtFmrHN8ok1NoLJZZmaTpBFAbo/bgS8f\newvP85YBaAbww+TDSp5imVmgMjOzjY1inX78OP6x8jFhMVaEaQJ1+bLdmdnBwXFM3HAVdRgaSnCm\nEiGEEFKhuHRf9bzgdbOPH4sASrtXCZCszVi3mLWlblYlfQpRaTWzNQb+zUKJfL/Ic78G4M993y/2\nc3z729/++/9/77338N5778WJLTF8v3gDKMC+zKyOmtmqqgmrcXt7/OPlYmIsj2LjRuD73w/2XNtt\nxo2NVQCG8eaNdxgNDewNRwghhITFtfuqWqf92I+Vft758zKSx4QjLkmb8eXLwI//uJ5j2TJr9skT\n4OXL4nXALtiMP/nkE3zyySeBnmtCzA4AWJbz9yUAhoo892sAfrXUwXLFrM0MD8sFoLa28M+bm4G/\n/dtUQyrI8+fAixfA3Ll6jqesxrrFrIlOxooNG8Rq4/ulayyeP5eh6cV2xmxg374PceRIF3p6lCVq\nGC0tXdi3b6/hyAghhBD3cO2+GrRu1lTzJ0DWpOPj0itlzhy9x67EzKzKyhZbo9bX2xFnKfITlN3d\n3UWfa0LMHgOwwvO8JgA3IYL16/lP8jxvNYDZvu8fSTm+RChlMQbsycwODkpTpbiF8Iqk6mavXwf+\n4T/Uf9wgLF4sQvb2bWDRouLPu3oVaGmRpgC2ksk04fDhvejs3I+hoXE0NFRh3769yGSaTIdGCCGE\nOIdr99XVq4GDB8s/z1S9LCBrUlU3u22bvuO+fCk1rrqSI42NMlvYNKXqZQHRI2HGTNpO6mLW9/0x\nz/N+DcAPIDW7f+T7/iXP87oBHPN9/6++fOrXIM2hKoJSzZ+AiZrZctm+pNFVL6tYswb4bgLtu0zW\nzHreRBOoUmLWdouxIpNpwscfd5kOgxBCCKkIXLqvBp08cfYs8PM/n3w8xVB1szrF7BdfyHEnTdJz\nPNsys8VgAygN+L7/fd/3V/u+v9L3/d/98rGuHCEL3/e7fd//v8ody+b5XbmUy8zOmSMi6eHD9GIq\nRBJiVndm1vfL7zolTZAmUK6IWUIIIYS8m6xYAfT2AqOjxZ/j+2ZtxkAydbM6LcaAPWK23BrZpQZQ\n2Wwv9uwpbjEGDIlZndg8vyuXUs2fFDZ0NB4Y0NP8SbFypeyklbpIhuXuXWDqVGDWLH3HDEuQ8Ty2\nj+UhhBBCyLvNlCmy7ivVLbivD5g+vfw6NkmSGM+ThJgdGhLxb5JymVkXGkABImR37z6AQ4c+Kvk8\n58WszfO7cilnMwbsqJvVnZmdOlWOp/P3Mtn8SREkM2v7WB5CCCGEkHLjeUxnZQE3xOy0afLHtIU3\nSGbWdIxB6Ow8mNNIrTgVIGYBW+d35VLOZgzYkZnVLWaB4PUYQTFZL6tobZWL4OvXhX/u+9IAimKW\nEEIIITZTbp1msvmTIolZs0mUg5m2Gvt++czs7NnAs2fF17C2UHhm89tUiJi1d36X4l3NzAL662ZN\nzphV1NWJLefq1cI/HxqS3bnZs9ONixBCCCEkDOXG89ggZpuagP5+YGxMz/F8PxkHnWkxq3rvlBph\nVFUl69MHD9KJKSoTM5tLY7cCDISa3/Wh6UBK4kpmVnfNLKBfzNpgMwZKW41ZL0sIIYQQFyiXmbXB\nZjx1KrBggaxTdTA0JIkJ3XNrlywxK2bLzZhVuNAEat++D9HS0oVygtZ5MfuzP7sfhw/bO79LEaQB\nlOnM7NgYcOeOzFHVSaWK2VJNoFgvSwghhBAXKFUz++qVrLts2KDXWTeru15WYTozG3TahwtNoNTM\n5kxmf8nnOS9m/+2/7bJeyALBbMa5s2ZNcPs2MHeuvnlbijVrgEuX9P1epsfyKDZuLJ2ZpZglhBBC\niO0sXChTJwqJm0uXZHzPlCnpx5WPzqRPkmJWV/Y4CuXqZRWuNIHKZJqwenXpmc3Oi9k7d0xHEIwg\nNuMZM4DaWhk9Y4Ik6mUBYP58+a+OL83ICHDzJrB0afxjxWXDhuKZWdqMCSGEEOICnle8bvbsWfMW\nYwUzs+WppMysopzWc17M3r5tOoJgBMnMAmbrZgcH9dfLAhMXSR1W474+uVDozh5HYfly+VwfP377\nZ7QZE0IIIcQVilmNbWj+pHBBzNpSM1sOVzKzQHmt57yYdSEz+/KltL+ePr38c03WzQ4MJJOZBfSJ\nWRvG8iiqqmREz/nzbz7+4gVw61awiwkhhBBCiGmKNYGyofmTQud4nnc9M+tCAyhAShSZmbUA1fyp\nXGcxwHxm1nYxa8NYnlwKNYG6dk1irKkxExMhhBBCSBiKiVnbMrM6Ej5Pn8rafNmy+MfKp74eeP5c\n/qRNkBmzCldsxo8eyajLUlDMpkBQizFgNjObpJgt1SkvDLZ0MlYUagLFellCCCGEuEShmtm7d8Vd\nmEQJWhQWLRIh+uxZvONcvQqsXCkOO914HtDQYCY7e/eujDCaObP8c12xGd++LSOZSuG8mHXBZhyk\n+ZPCdGY2qQtWuYHcQbFNzBZqAsV6WUIIIYS4xIoVkkwZHZ14TFmMgzgL08Dz9CR9kk46mKqb6fzP\njgAAIABJREFUDZqVBdzJzN65I922S+G8mGVmVh9J1swuXy7Hf/ky3nFsqpkF5CJ/7tybY4c4locQ\nQgghLjF1qmQUc9egNlmMFTrqZpOql1WYqpsNM7qSmVmLcCUzG1TMNjVJx97x8WRjysf3k7UZT5ok\nX7Avvoh3HNtqZuvrpbFXX9/EYxSzhBBCCHGNfBedTc2fFDrqZitVzIbJzLrSAOqdyMw+efKmJcJG\nwtiMp00DZs2Sbrhp8uSJ2DeC+Oyj8v+3d+/BcZ3nfcd/D0ASIsEbQIAgAYgURFmS5fgSpVbjJOow\ndVWrmUzUZOLYTtRY7SRNZmw5baKpcxkMxTK9uHGTSZQ2Eydp6EZylMS1Y9ft2FZTM6nS2lQsKVYs\niVRIECYIEoIA8AbwAgFv/3j3CMvFAtizey7v2f1+ZjgEFnt5pcMlz7PP8/5OoyFQMzM+FbrWDway\nUh4C5RzFLAAAKJ7KfJMQO7NJXJ4ni2J2bCy9519JnM5sd7cPV8q6eRbXxEQLFLM9PeF3Z+OMGUv5\n7JtNsysbabSYjbqyoezdiJSHQJ0960d1urvzXRMAAEAc5YnGCwvSiy9K3/Zt+a6pUqPF7MKCv+rE\n7bcnt6ZKRejMrlvnJwvPn09zRY1riTHjvr7wi9k4nVkpn32zY2Ppp9UlUcyGtF82Uh4CRVcWAAAU\nUXkxe+KEP8fesiXfNVVq9Bx5dNQXR52dya2pUl4BUHHPk4sQAtUSY8Z9feGHQMXZMyvRmV1JaEnG\nkfIxYy7LAwAAiqh8z2yII8bSUjFbHrwZR9ojxlI+ndnFRZ/fsndv7Y8pwr7ZlujM7twZfmc27phx\nHp3ZLIrZ6BO/ev8CCrWYvfNOf7yuXuWyPAAAoJh27fLnMtPTYYY/SX40dsuW+rNlsihmd+/2RdjC\nQrqvU+7cOZ97E6fjXIREYzqzgYg7Ztysndnt2/1fQvV+WhVqMbthg78+20svMWYMAACKyWyp8RBq\nZ1ZqbN9sFsXshg1SV1e29Umc/bKRIowZt0QA1M6d4RezRenMpr1nVlqelBdHqHtmpaUQKIpZAABQ\nVFEx+8IL4RazjZwnZzVBl/W+2XrOkUPvzM7N+SvWrLVvu/DFbOgBUPPz/mBs21b7Y/bs8YFMWY4n\njI2l35mVll/DrFYLC34vQNxPnbLy1rdKR49K4+PhFtwAAACrufNO6etf91dnuO22vFdTXSOd2ayy\nTbLeN1tPZzb0PbPRiPFaVzFpimI25M7s9LQfNWiL8X+6o0Pq7c32TZDFmLFUfwjUmTN+HOKmm5Jf\nUxLe9jbps5/1hez69XmvBgAAIL477pA+8xnpzW+W2tvzXk119Raz09PSlSt+T2vasi5m6+nMhj5m\nXEv4k9QExWzoAVBxR4wjWe6bvX5dmpmp7Q9Mo+otZk+eDLvjuW3bqMbHD2py8oAefPCgRkZG814S\nAABALJ2d/nzmzJlwz2eGhuorZqOu7FqdviQMDPipx6zU25kNecy4lvAnqQmK2dA7s3HDnyJZ7ps9\ne9Yn2GXxCVy9xezISJjhT5I0MjKqD37wMUmPaGrqoJ544hHdd99jQf4DAAAAUM3IyKg+/GF/PnPu\nXLjnM7feWt85chbhT5Ei7JktQme2JYrZ3l5pctJfXylERejMZrVfVvL7gaempEuX4j0u1CRjSRoe\nPqwTJw5KivLQO3XixEENDx/OcVUAAAC1Gx4+rJGR8M9nBgf9uf/Vq/Eel2Uxm+WY8cKCP5ffsyfe\n40LvzLbMmPGGDT7lamYm75VUV4TObFb7ZSW/d/j226Xjx+M9LuRi9syZRS39xR/p1Ph4oJ+wAAAA\nVCjK+Ux7u3TzzdJozIZxsxaz9ebKFCUAai2FL2alsEeNp6bC78xmWcxK9Y0ah7xndmCgTdJsxa2z\n6u9vircXAABoAUU6n6ln32zWxezYmORc+q9Vz35ZaamYzWKN9WiZzqwUdghUvWPGzdqZleorZkPe\nM3vo0EPat++Alv4BmNW+fQd06NBDua0JAAAgjiKdz8TdN3v9uu/k7tuX3prKbd3qf794Mf3Xqme/\nrOSvntLREX/rX1Zq7cyuS38p6Qu9M3vXXfEfNzgonTvnr1Ob9qVexsaku+9O9zXKRbHvtZqdlS5c\n8CFVIRoa2qunnnpYw8Mf1/j4ovr723To0MMaGtqb99IAAABqUqTzmbiX5zlxwo8md3Skt6ZyZksh\nUNu2pfta9XZmpaUQqKj4DkmtAVBNUczu3BluMVtvZ3b9en8drNOn0+9Iht6ZjT5xinOt3qwNDe3V\n448fyHsZAAAAdSvK+cytt0pf/Wrt989yxDgS7Zutp6kVx8iIdO+99T02CoEKcStfS40Z9/WFO2Zc\nbwCUlN2+2ayL2dtvl/72b336Wi3qHZ8AAABA84m7HS+6xmyWsrrWbCOd2VBDoF5/3U9l1tIQbJpi\nNtTObL0BUFI2+2adk8bHsy1mOzv9Jy21ptCFnGQMAACAbEVjxrWGF+XZmU1bI02fUK81OzkpdXf7\n5Oq1NEUx24wBUFI2ndmpKV9cbtyY7utUijNqTDELAACASFeX3342PV3b/fMoZqM9s2man/cZOzff\nXN/jQ73WbK3hT1KTFLOhdmYXFnyLvKurvsdn0ZkdG8u2KxuhmAUAAEC9ag2Bcq55O7OnT/uA1HrD\nYkPtzNYa/iQ1STEbagDUzIxPB1tXZ8xWFp3ZrPfLRuIUs+yZBQAAQLlamz4TE77Yq3dSsl5ZFLOn\nTjV2jhxqZ7bW8CepxjRjM1sn6b2S3lW6qVPSgqQ5Sd+Q9Cnn3NXYK01IqAFQjYQ/Sdl0ZvMsZv/w\nD9e+n3MUswAAALhRrZ3ZPLqyUjYBUCMj9Yc/SeEGQMUZM16zmDWzd0q6V9JTzrll5YeZ7ZP0z83s\nr51zfx5zrYnYvNn/fvny0tchaCT8SZL6+/2nJdeupXddrDNn/Ex/1mrtzE5M+D29W7akvyYAAAAU\nw623Ss8+u/b98ipm+/r8lOb169KGDem8RqOd2ZDHjGvtzNYyZnzVOferzrkXqv3QOXfCOfcbkk6b\nWUqHam0hhkA1Ev4k+QSvm2+uPfW3Hnntmd21S7p6de030MmTdGUBAABwo6Gh2juzd9yR/noqtbf7\ngvbs2fReI4nObIhjxokGQK1UxFa530nn3PXaXjZ5IYZANTpmLKW/bzavMWMz/ynZsWOr329khPAn\nAAAA3OjWW2vbjpdXZ1ZKf98sndk6AqDM7D2l3/+FmX3QzB6I+xxpCDEEqtHOrJT+vtm8ilmptmKW\nJGMAAABU2rvXTxi+/vrq9zt2LN9iNs19s83amU01zdg596XSl0+XfuU2WlwuxBAoOrOrq2XfLMUs\nAAAAKm3Y4LetnT698n3m5vx1WBsp+BqRZmf22jVfiDZyHr9pk/99bi6ZNSUllevMmtk9ZvZHZvY1\nM/uypI2l/bJ/Uuc6ExXqmHHIndm5Of8r66jySC3FLEnGAAAAqGat8+RXXpH27av/MpmNGhxMr5gd\nHfXP397e2POENmrsnC9me3tru3+czuxu59z7nHN/V9JDku41s3vjLzEdzRgAJaXbmY26smbpPP9a\n6MwCAACgXmtdnifP/bJSup3ZRvfLRkIbNT5/3neMb7qptvvHKWZPlPbJ9jnnxp1z/1ZSVz2LTEOo\nndlGx4zT7MzmOWIs+U/KRkd9ZHk11675Y5rHpYMAAAAQtlYuZhvdLxsJ7VqzccKfpBjFrHPubyR9\nRtK/MrPPmdmfS3qnmb3TzNbHXmnCmjUAqq9PungxnVn2vIvZjg5pzx7pxInqPx8d9Zcmyms0BAAA\nAOEqQjGbVgBUUp3Znp6wOrNxwp+kmAFQzrlvOed+zjn3gKQflPSXkv6RpE+aWY3N4HQ0awBUW5tP\na0tj1PjMmfy7nquNGrNfFgAAACtZa4IxhGJ2fNzvA01as3Zm44Q/SWsUs2bWYWZVe4vOuWnn3Bed\nc//aOfejkmrcppuO0DqzzknT01J3d+PPlda+2bGxfDuz0urFLPtlAQAAsJLVOrOLi9Lx49Idd2S7\npnKbNvlfaRSLSXZmQypmEx0zds5dk/QuM/uAmW2sdh8z225mPylpb5yFJq27W7p8eeX9l1m7cEHa\nuNHHhjcqrX2zeY8ZS/4vGIpZAAAAxLVzp9+Kd/Hi8p+dPi1t3y5t2ZL9usqltW82yc5sSGPGcTuz\na+5GdM59wcx2SfqXZrZT0k2S1kt6XdIVSacl/a5z7kJdK05IW5v/ZGFyMv8CTUpmxDiSVmc2hGL2\nzjulT3yi+s9OnpTuuSfb9QAAAKAYzHzjY2REevvbb/xZ3iPGkWjfbOX6GhEV8Lt2Nf5cO3ZIR482\n/jxJmZiQ3vGO2u+/ZjFrZg9Ier6UXhy0aNQ47wJNSib8KTI0JD3zTDLPVS6UPbPHjvmx7MpLBLFn\nFgAAAKuJJhgri8Vjx8IoZtO41uypUz5EtS1W+lF1rRAAdYuk95vZd0uSmf2ImX3IzHIdK64mpBCo\n0DuzCwv+D8vu3ck+b1w7dkjr1y/f7+ycTzlmzBgAAAArWWnfbEid2TSK2aQaPk0dAFUy7Zz7mHPu\nL83sQ5J+SdKEpIfM7D31LTMdIYVATU0l25lNes/sxITfZ7w+94sqVQ+BmpnxndquYK5kDAAAgNC0\nYjGb1H5ZqckDoErK83jfJ+nXnHOfds4dlJRjPthyIXVmkxwz7umRrl2rvrm9XiHsl41UK2ajEePK\n0WMAAAAgslLTJ6RiNulrzSa5Fa/oAVC1FLNPm9knzOxxSbdL+mzZz67EW166+vrC6swmNWZslvyo\ncQj7ZSPVilmSjAEAALCWap3ZCxd8EyiExk1ae2aT6sxu2eKbZiFcEWZuTpqfj5dAvWYx65z7uqSP\nSPoNSbc75y6Y2XeY2QckBdU3a9YxY8n/gU1y1DiEa8xGKGYBAABQj6EhX9wtLi7dduyYv/xjEgFJ\njUprzDipzqxZOPtmo/CnOJOZNR1i59xV59xR51w06PoNSSclzcVeZYpCGzNOqjMrLb1RkxL6mDHF\nLAAAANayaZO/nuzZs0u3hTJiLPlCcW5OupLgPGuSnVkpnFHjuCPGUo3FbCXn3Lxz7mvOucfreXxa\n6MzWLqRi9pZb/HGbK/tohMvyAAAAoBaV+2Zfftl3ZkNgJvX3J9edvXhRunpV6u1N5vmkcEKg4oY/\nSXUWs6EKrTObZDHbzJ3Z9nZp3z7p+PGl2+jMAgAAoBaV+2ZD6sxKyYZARV3ZJENSW64zG6reXmly\n8saZ+bwkGQAlpdOZDSUASrpx1Pj116XTp6W9wV3JGAAAAKEJvZhNMgQqjenFkPbMtnRndsMGn341\nPZ3vOpxLfsw46sw61/hzORdWAJTk/8I5dsx/PTbmP5Xp6Mh3TQAAAAhfeTH7+uu+4HvTm/JdU7kk\nQ6CS3i8rhTVm3NKdWSmMUePZWT86u3Fjcs+5fbv/fWam8ee6eNGPJmzd2vhzJaW8M8t+WQAAANSq\nfM/syIi0e3ey5+GNSrKYTaszy5hxIEIIgUq6Kyv54jOpfbMh7ZeNlBez7JcFAABArco7s6GNGEvp\n7JlNUkid2ZYeM5bC6MwmHf4USWrfbGj7ZSWfOHf8uN/vTDELAACAWvX3+2LsypUwi9ki7JmlMxuI\nvr4wOrNJhj9FkurMhrZfVvJ7nbu6fPDTyAjFLAAAAGrT3u6DQ0dHwyxmkxozdi6dziwBUAFp1jFj\nKdnObGjFrLQ0anzyJHtmAQAAULuhIX8OGWIxu3u3r08WFhp7nig7p6ur8TWVC2HMeH5eunAhfg3V\ndMVsM48ZN/OeWenGYpbOLAAAAGp1663SiRPSSy+FV8xu2OAL0EZrlDSuMSuFMWYc1U/t7fEe13TF\nbCid2TTGjJt5z6zk983+1V/5NOi48/IAAABoXbfeKh096kdxe3vzXs1yg4ONh0CldcWP7duly5f9\nZY3yUs+IsdSExWwzd2ZvuSWZa82GuGdW8p+ifelL6XziBAAAgOY1NOTPI++8M8zzyCT2zaaxX1aS\n2tp8QTs9nfxz16qe8CepSYvZZu3Mbt0q3XSTNDnZ2POEOma8adOoJicPamLigB588KBGRkbzXhIA\nAAAKoKPDn0eeOhXmeWQSxWxanVkp/xAoOrMl0Zhxo93LRqQVACU1vm/2+nW/ebyePyxpGhkZ1Y//\n+GOSHtHU1EE98cQjuu++x4L7iwgAAABhGRkZ1Uc+4s8jz50L8zwy5M6slH8I1MQEnVlJ0ubNfrRg\ndja/NaQ1Ziw1vm/27Fn/ByXu5uq0DQ8f1smTByV1lm7p1IkTBzU8fDjHVQEAACB0w8OHNTIS9nnk\nwEC4e2al/EOgGDMuk3cIVFpjxlLjndlQw5/OnFnU0l9AkU6Njy/msRwAAAAURBHOIwcHG+vMpnWN\n2QhjxgHJOwQq5M5sqOFPAwNtkirb6bPq72/KP6IAAABISBHOIxsdM56c9Nk5W7cmt6ZyPT10ZoOR\nZwjU1as+1nrz5nSeP4nObIjF7KFDD2nfvgNa+otoVvv2HdChQw/ltiYAAACErwjnkY0Ws2l2ZaXi\ndmbXJb+U/OU5ZhyFP6UVCd5oZzbUYnZoaK+eeuphDQ9/XOPji+rvb9OhQw9raGhv3ksDAABAwIpw\nHrl1qx8Vvnixvu5qmvtlJd+ZPX48vedfS70BUE1ZzOY5ZpzmiLHki9lvfUtaXPTXhIrrzBnp7rsT\nX1Yihob26vHHD+S9DAAAABRM6OeRZn7f7NiYdNdd8R+fRWc2rzFj5/wYNXtmS/LuzKYV/iRJmzb5\nT3POnavv8aHumQUAAACaWSOjxml3ZvMcM56Z8TVOR0f8xzZlMdvMnVmpsX2zoY4ZAwAAAM2skWI2\n7c5sngFQ9YY/SU1czDZrZ1aqf9+sc9L4OMUsAAAAkDU6s9XVG/4kNWkxm/eYcaid2akpqbNT2rgx\n8SUBAAAAWEW0ZzauxUWfmbM3xTyr7m7p/Hn/WlmrN/xJatJittnHjOvtzLJfFgAAAMhHvZ3Zc+d8\nZk5nZ/Jriqxb5y8tev58eq+xEsaMK3R1SZcuSdevZ//aWYwZ19uZZb8sAAAAkI96i9m098tG8ho1\nZsy4Qlub1NvrI56zFnJnlmIWAAAAyEe9xWza+2UjeYVA0ZmtIq8QqCw6s3v3+pHhhYV4j6OYBQAA\nAPLR1ydNT8efHqUzu7KmLWbzCoHKIgCqo8MXzHE/2Tlzxm88BwAAAJCt9nZf0J49G+9xWXZm8ypm\n6cxWyCsEKosxY6m+fbMEQAEAAAD5qWfUOMvOLGPGgcijMzs/L83NSdu2pf9a9eybZcwYAAAAyE89\nxWxWnVnGjAOSR2d2etonKbdl8H+1ns4sxSwAAACQn7jF7MKCn67csye9NUXyGDOenfX/jVu21Pf4\npi5ms+7MvvZa+uFPkbid2bk56cqVbEagAQAAACw3OOiL01qdOePri5tuSm9NkTzGjF991Xdlzep7\nfNMWs3mMGWcR/hSJ25k9c0bq76//DwoAAACAxsTtzGa1X1bKpzPbSPiTlFMxa2b3m9nLZnbczD66\nwn1+xMy+aWYvmNnjcV8jjzHjrMKfpPidWUaMAQAAgHzFLWaz2i8r5deZbaSYXZfcUmpjZm2SflPS\nuyWNS3rGzD7nnHu57D63SfqopHc55y6aWezh3bw6s1mNGd98s3TunA+dWr9+7ftTzAIAAAD5Crkz\nm0cAVCPhT1I+ndl7JL3inBt1zs1LelLSAxX3+UlJ/8k5d1GSnHOxPyPYudN/srC42PB6a5blmPH6\n9dKuXdLp07Xdn2vMAgAAAPmKilnnart/1p3Zqana15aERjuzeRSzA5LKS7Cx0m3lbpd0h5k9bWb/\n18zeE/dF1q+Xtm71CcNZyTIASoq3b5ZrzAIAAAD52rTJ/6q1A5plZ7ajw/+6dCmb15Ma78xmPmYs\nqVoEUWX9v07SbZL+nqQ9kv6Pmb0l6tSWe/TRR9/4ev/+/dq/f/8b30ejxlkVmFNT0l13ZfNaUrx9\ns2fOSPfem+pyAAAAAKwh6s7WUqNk2ZmVlkKgtm7N5vUmJqR3vevG244cOaIjR47U9Pg8itkx+QI1\nMii/d7byPv/PObco6ZSZHZP0Jklfr3yy8mK2UhQC9Za3NLrk2mQZACXF68yyZxYAAADIX1TMvv3t\nq99vft5n5Nx8czbrkpZCoLIqoKuNGVc2KA8ePLji4/MYM35G0m1mttfMNkh6v6TPV9znTyX9fUkq\nhT+9SdLJuC+UdQhUlgFQUvzOLHtmAQAAgHwNDtYWAnX6tM/IqSXsNSlZh0AVLgDKObcg6cOSvizp\nm5KedM69ZGYHzez7S/f5kqQpM/umpD+T9Ihzbibua2V9eZ4sA6Ck2juzCwv+D8ru3akvCQAAAMAq\nBgZ8ns1aTp3KdsRYyv5as4W7NI8kOee+KOmOitsOVHz/c5J+rpHX6evLtjObdQBUrZ3ZiQmpuzvb\nT3UAAAAALDcwIB09uvb9RkayC3+KZHmt2fl56cIFX6fUK48x48xkOWa8sOAPRldXNq8n+TfCa69J\n166tfj/2ywIAAABhqPVas3l0ZrMcM56c9K/X3l7/czR1MZvlmPHMjE/9auRgxNXe7mfuR0dXvx/F\nLAAAABCGWovZPDqzPT3ZdWYbHTGWmryYzbIzm3X4U6SWfbOEPwEAAABhGBwMd89slp3ZRsOfpCYv\nZrPszGYd/hSpZd/s2BidWQAAACAEO3ZIc3PSlSur3y+vzmxWxSyd2TVEAVDOpf9aWYc/RWrtzFLM\nAgAAAPkzk/r7Vx81vnbN1xdZn8NnGQBFZ3YNnZ3+D8vsbPqvFXJnlmIWAAAACMda+2ZHR/04cpZ5\nPFK2Y8Z0ZmuQ1eV5Xnstn2KWPbMAAABAsQwOrl7M5rFfVso2AGpigmJ2TVmFQOUVALVWZ9Y59swC\nAAAAIRkYWD0EKo/9spK0aZP/fW4u/ddizLgGWYVA5TVmvGuXdPHiyn/gLl70o9Zbt2a7LgAAAADV\nrTVmnFdnVsouBIox4xpk1ZnNKwCqrU3as2flUWP2ywIAAABhWauYzaszK2UXAkVntgbN3pmVVt83\ny35ZAAAAICyh7pmVsgmBck6anKSYXVOzB0BJq++bZb8sAAAAEJZQ98xK2YRAzcz4K890dDT2PE1f\nzDZ7AJS0dmeWYhYAAAAIx+7dvkZZWFj+s7k5n3uza1f265Ky6cwmMWIstUAxm8WYsXPS9LTU3Z3u\n66xktc4sxSwAAAAQlg0bpK6u6nXKqVM+E6ctp0otiwCoJMKfpBYoZrPozF64IG3c6P9Q5mFoiGIW\nAAAAKJKVQqBGRvLbLytlEwBFZ7ZGWXRm8xwxlnxnlgAoAAAAoDgGB6vvmz11Kr/9slI2Y8Z0ZmvU\n1SVdvixdv57ea+QZ/iRJvb3S1at+tr4SAVAAAABAeELtzGYRADUxQTFbk7Y2f0DS7M7m3Zk1q96d\nvX7dJ4Ul0cIHAAAAkJyVitlW6MwyZhxD2qPGeV5jNlJt3+zZs/6/vb09nzUBAAAAqC7kzixjxgFJ\nOwQq7zFjqXpnlv2yAAAAQJgGB8PtzBIAFZAsOrN5jhlL1Tuz7JcFAAAAwjQwsDwA6uJFn4XT25vP\nmiRpyxbp2rV0M4fozMbQ15duZzaEMeOVOrMUswAAAEB4qo0ZR11ZszxW5Jmlv2+WAKgYshgzDrEz\nSzELAAAAhGnrVsm5G69Ikvd+2Uiao8azs9LCgrR5c+PP1RLFbCsEQEWdWeeWbmPPLAAAABAms+X7\nZvPeLxtJMwQqGjFOovvcEsVsKwRAdXX532dmlm5jzywAAAAQrsp9s63QmU0q/ElqkWK2FQKgql1r\nljFjAAAAIFyV+2ZD6cymuWc2qfAnqYWK2bQ6s86FMWYs3bhv1jlpfJxiFgAAAAhVZTEbSmc2zTHj\npMKfpBYpZnt7fZt8cTH5556dldrbpY0bk3/uuMo7s1NT0qZNYawLAAAAwHLlxaxzYXVmGTMOxPr1\nPi1sejr55w6lKyvd2Jkl/AkAAAAIW3kAVJR9E2Xh5Ikx48CkFQIVQvhTpLwzS/gTAAAAELbyAKgQ\nrjEb6emhMxuUtEKgQgh/ilR2ZilmAQAAgHCVjxmHsl9WojMbnLRCoEIaMy6/1izFLAAAABC2vj6/\nFXJ+Ppz9slL6AVB0ZmNKc8w4lM7s1q3STTdJk5PsmQUAAABC197uC9qzZ8PrzKY5ZkxnNqY0x4xD\n6cxKS91Z9swCAAAA4Yv2zYbUmd2+Xbp0SXr99WSfd35eungxufqpZYrZVgiAkpb2zTJmDAAAAIQv\n2jcbUme2rc2nKid9NZjJSV87tSVUhbZMMdsKAVDSUmeWYhYAAAAIX1TMhtSZldIJgUoy/ElqsWK2\n2QOgJP9pzosvSleuhLUuAAAAAMsNDkrPPeezb7ZuzXs1S9IIgUoy/ElqoWK2FQKgJP9pztNPS/39\nYVyjCgAAAMDKBgb8+XtIXVkpnRCoJMOfpBYqZqMxY+eSfd4QO7MnTzJiDAAAABTBwIA/fw9lv2yE\nMeOAdHb6TuXly8k+b2gBUM6NSjqoY8cO6MEHD2pkZDTvJQEAAABYweKiP38/ejSs8/eennQ6s4wZ\n1ynpEKirV31c9ebNyT1nI0ZGRvXAA49JekSvvnpQTzzxiO6777Fg3hAAAAAAloyMjOonfsKfv58+\nHdb5O53ZwCQdAhWNGIeyN3V4+LBOnDgoqbN0S6dOnDio4eHDOa4KAAAAQDXDw4c1MhLm+TsBUIFJ\nOgQqtPCnM2cWtfRGiHRqfHwxj+UAAAAAWEXI5+8EQAUm6THj0MKfBgbaJM1W3Dqr/v6WOswAAABA\nIYR8/s6YcWDS6MyGVMweOvSQ9u07oKU3xKz27TugQ4ceym1NAAAAAKoL+fw96QCoxUVpclLq7U3u\nOdcl91Th6+uTjh1L7vmmpsIaMx4a2qunnnpYw8Mf1/j4ovr723To0MMaGtqb99IAAABw1We8AAAU\nU0lEQVQAVAj5/D3pzuz58/4KMx0dyT1nyxWzf/EXyT1faGPGkn9DPP74gbyXAQAAAKAGoZ6/d3f7\nAnRxUWpLYJ436fAniTHjhoQWAAUAAAAASVi3zl+C9Pz5ZJ4v6fAnqcWK2WYPgAIAAACApCQ5apx0\n+JPUYsVsswdAAQAAAEBSkgyBYsy4QV1d0uysdP16Ms8XWgAUAAAAACSFzmxA2tp8FHRSo8aMGQMA\nAABoVj09yRWzdGYTkOSoMQFQAAAAAJrVjh3JjhnTmW1QUiFQ8/PS3Jy0bVvjzwUAAAAAoWHMODBJ\ndWanpvweXLPGnwsAAAAAQkMAVGCS6swS/gQAAACgmdGZDUxfX3KdWcKfAAAAADSrpIrZ2VlpYUHa\nvLnx5yrXcsVsUmPGhD8BAAAAaGZJjRlH4U9Jb9FsuWI2yTFjOrMAAAAAmlVSndk0RoylFixmkwyA\nopgFAAAA0KyiYta5xp4njfAnqQWL2aQ6s4wZAwAAAGhmHR3+16VLjT0PndmE9Pb6QnRxsbHnoTML\nAAAAoNklMWpMZzYh69dLW7c2fkBee41iFgAAAEBzSyIEKgqASlrLFbNSMqPGXGcWAAAAQLNLojPL\nmHGCkgiBYswYAAAAQLPr6WHMOChJdGYJgAIAAADQ7HbsaHzMmM5sgvr6GuvMLixIFy5IXV3JrQkA\nAAAAQkMAVGAaHTOemfEhUu3tya0JAAAAAELTaADU/Lx08WI6WzRbsphtdMyY8CcAAAAAraDRzuzk\npK+d2lKoPFuymG20M0v4EwAAAIBW0GgAVFojxlKLFrONdmYJfwIAAADQChoNgEor/Elq0WKWziwA\nAAAArK3RMWM6swmL0oydq+/xr71GMQsAAACg+TUaADUxQWc2UZ2dPon48uX6Hk8AFAAAAIBWsGmT\n/31urr7HM2acgkZGjRkzBgAAANAqGgmBYsw4BY2EQBEABQAAAKBVNBICRWc2BXRmAQAAAGBtjYRA\n0ZlNQRQCVQ8CoAAAAAC0ikZCoAiASkEjY8YEQAEAAABoFfV2ZhcXpclJOrOJq3fM2Dlpelrq7k5+\nTQAAAAAQmnoDoGZmpM2bpQ0bkl+T1MLFbL2d2QsXfDx1WgcEAAAAAEJSbwBUmuFPUgsXs/V2Zgl/\nAgAAANBK6h0zTjP8SWrhYrbeACjCnwAAAAC0knoDoNIMf5JavJitZ8yY8CcAAAAAraTezixjxinZ\nvl2anZWuXYv3OMaMAQAAALQSxowD09Ym9fb6qOg4XnuNziwAAACA1lHvmDGd2RTVEwJFZxYAAABA\nK9myxU+0xp1qpTObonpCoAiAAgAAANBKzOobNSYAKkX1hEARAAUAAACg1dRTzDJmnCLGjAEAAABg\nbfV2ZhkzTkk9nVkCoAAAAAC0mrghULOzknPS5s3prSmXYtbM7jezl83suJl9tMrPP2hmr5rZs6Vf\n/yyNddCZBQAAAIC1xe3MRl1Zs/TWtC69p67OzNok/aakd0sal/SMmX3OOfdyxV2fdM59JM21xA2A\nco5iFgAAAEDr6emJX8ymuV9Wyqcze4+kV5xzo865eUlPSnqgyv1SrOG9uGPGs7NSe7u0cWN6awIA\nAACA0OzYEW/MOO3wJymfYnZA0umy78dKt1X6ITN73sz+2MwG01hI3DFjurIAAAAAWlG9Y8ZpynzM\nWNU7rq7i+89L+pRzbt7MfkrSJ+XHkpd59NFH3/h6//792r9/f80L6e31B2RxUWqroawn/AkAAABA\nK4obAFVvZ/bIkSM6cuRITffNo5gdk7Sn7PtB+b2zb3DOzZR9+zuSPrbSk5UXs3GtXy9t3eoL2t7e\nte9PZxYAAABAK6qnM3vbbfFfp7JBefDgwRXvm8eY8TOSbjOzvWa2QdL75TuxbzCzXWXfPiDpxbQW\nEycEimIWAAAAQCsKMQAq886sc27BzD4s6cvyxfTvOedeMrODkp5xzn1B0kfM7AckzUualvRQWuuJ\nEwLFmDEAAACAVhRiAFQeY8Zyzn1R0h0Vtx0o+/oXJf1iFmuJEwJFZxYAAABAK9q+Xbp0SXr9dWld\nDVVkFgFQeYwZB4XOLAAAAACsrq1N6uqSpqdru3+zXponKHRmAQAAAGBttYZAzc9LFy9K3d3prqfl\ni1kCoAAAAABgbbWGQL36qr9vLZc/bQTFLGPGAAAAALCmWkOgshgxlihmGTMGAAAAgBrUOmacRfiT\nRDFLZxYAAAAAatDTQ2c2KFFn1rnV73f1qo+h7uzMZl0AAAAAEBI6s4Hp7JTa2/01k1YTjRibZbMu\nAAAAAAhJrQFQExN0ZjNTy6gxI8YAAAAAWhkBUAGqJQSK8CcAAAAArYwx4wDV2pmlmAUAAADQqgiA\nClCtnVnGjAEAAAC0KjqzAerrY8wYAAAAAFbT3S3NzEiLiyvfZ3FRmpykmM0MAVAAAAAAsLp166Qt\nW6Tz51e+z8yMv8+GDemvh2JWBEABAAAAQC3WGjXOasRYopiVRAAUAAAAANRirRCorMKfJIpZSQRA\nAQAAAEAt6MwGhgAoAAAAAFhbLcUsndkMdXVJc3PStWsr34cAKAAAAACtrpYxYzqzGTKTentX3jc7\nP++L3W3bsl0XAAAAAISEzmyAVguBmpry3VuzbNcEAAAAACEhACpAq4VAEf4EAAAAAARABWm1ECjC\nnwAAAACAMeMgrTZmTPgTAAAAABAAFaS1xozpzAIAAABodat1Zi9flpyTNm/OZi0UsyVrdWYpZgEA\nAAC0uqiYdW75z6Lwp6yCcylmSwiAAgAAAIDVdXT4X5cuLf9ZluFPEsXsGwiAAgAAAIC1rTRqnGX4\nk0Qx+wYCoAAAAABgbSuFQGUZ/iRRzL6hp8d/urCwsPxndGYBAAAAwKMzG5j166Vt26Tp6eU/o5gF\nAAAAAC9qBFaKAqCyQjFbZqUQKMaMAQAAAMDbsaP6mDEBUDmqFgK1sCBduCB1deWzJgAAAAAICWPG\nAdq5c3kI1MyMHz9ub89nTQAAAAAQEgKgAlStM8t+WQAAAABYQmc2QNUuz0MxCwAAAABLqhWz169L\nly5J3d3ZrYNitky1ACjCnwAAAABgSbUx48lJf3tbhhUmxWwZxowBAAAAYHXVOrNZjxhLFLM3qBYA\nRWcWAAAAAJZU68xmHf4kUczegM4sAAAAAKxu0yb/+9zc0m10ZnMWdWadW7qNYhYAAAAAblQ5avzq\nqxSzuers9NeTvXRp6TbGjAEAAADgRpWjxhMTjBnnrnLUmM4sAAAAANyosjPLmHEAKkOg6MwCAAAA\nwI0qO7MEQAWAziwAAAAArI7ObID6+pY6s85J09NSd3e+awIAAACAkBAAFaCdO5c6sxcu+NjpDRvy\nXRMAAAAAhKR8zHhxUZqclHp7s10DxWyF8jFjRowBAAAAYLnyzuz0tLRlS/ZNQIrZCuUBUIQ/AQAA\nAMBy5Z3ZPMKfJIrZZejMAgAAAMDqyjuzeYQ/SRSzy5QHQFHMAgAAAMBy5cVsHuFPEsXsMuUBUIwZ\nAwAAAMBy5WPGExOMGQehq0uam5OuXqUzCwAAAADVbNkiXbvmfzFmHAgzHyk9OUlnFgAAAACqMVsa\nNSYAKiBRCBSdWQAAAACoLipm6cwGJAqBopgFAAAAgOrKO7MUs4GIQqAYMwYAAACA6qIQKAKgAsKY\nMQAAAACsbseOpWKWzmwgyjuzFLMAAAAAsFxPjzQ66r/u7Mz+9Slmq+jrk0ZGpHXrpI0b814NAAAA\nAIRnxw7ppZd8/WSW/etTzFbR1ye9+CJdWQAAAABYyY4dvm7KY8RYopitaudO6cQJwp8AAAAAYCU9\nPb5uyiP8SaKYraqvT1pcpDMLAAAAACvZscPXTXRmAxJ1ZClmAQAAAKC6qF6iMxuQ9ev9gWHMGAAA\nAACqi+qlvDqz6/J52bCNjIxqfv6wvvCFRZ0/36ZDhx7S0NDevJcFAAAAAMGYmRmVdFi//duLOno0\n+7rJnHOZvVjSzMwlvf6RkVHdd99jOnHioKROSbPat++AnnrqYQpaAAAAAFB2dZOZyTlX9cI/jBlX\nGB4+XHZAJKlTJ04c1PDw4RxXBQAAAADhCKFuopitcObMopYOSKRT4+OLeSwHAAAAAIITQt1EMVth\nYKBN0mzFrbPq7+d/FQAAAABIYdRNVGgVDh16SPv2HdDSgfGz34cOPZTbmgAAAAAgJCHUTQRAVTEy\nMqrh4cMaH19Ufz9pxgAAAABQKYu6abUAKIpZAAAAAECQSDMGAAAAADQVilkAAAAAQOFQzAIAAAAA\nCodiFgAAAABQOBSzAAAAAIDCoZgFAAAAABQOxSwAAAAAoHAoZgEAAAAAhUMxCwAAAAAoHIpZAAAA\nAEDhUMwCAAAAAAqHYhYAAAAAUDgUswAAAACAwqGYBQAAAAAUDsUsAAAAAKBwKGYBAAAAAIVDMQsA\nAAAAKByKWQAAAABA4eRSzJrZ/Wb2spkdN7OPrnK/HzazRTO7O8v1AQAAAADClnkxa2Ztkn5T0nsk\nvUXSB8zszir32yzpYUlfzXaFyNuRI0fyXgISxjFtThzX5sMxbU4c1+bDMW1OHNf48ujM3iPpFefc\nqHNuXtKTkh6ocr9Dkj4m6VqWi0P+eCM3H45pc+K4Nh+OaXPiuDYfjmlz4rjGl0cxOyDpdNn3Y6Xb\n3mBm75A06Jz7n1kuDAAAAABQDOtyeE2rcpt744dmJunXJH1wjccAAAAAAFqUOefWvleSL2j2nZIe\ndc7dX/r+5yU559zHSt9vlfS3ki7LF7G7JE1J+gHn3LMVz5Xt4gEAAAAAmXLOVW1u5lHMtks6Jund\nks5KOirpA865l1a4/1ck/axz7rnsVgkAAAAACFnme2adcwuSPizpy5K+KelJ59xLZnbQzL6/2kPE\nmDEAAAAAoEzmnVkAAAAAABqVR5pxIszsfjN72cyOm9lH814PGmdmp8zsr83sOTM7mvd6UB8z+z0z\nmzCzb5Td1mVmXzazY2b2JTPblucaEd8Kx/WAmY2Z2bOlX/fnuUbEY2aDZva/zexFM3vBzD5Sup33\na0FVOaYPl27nvVpgZtZhZl8rnR+9YGYHSrffYmZfLb1X/9DM8gh2RR1WOaa/b2YnS7c/a2Zvy3ut\noStkZ9bM2iQdl993Oy7pGUnvd869nOvC0BAzOynpO5xzM3mvBfUzs++RD3D7r865t5Vu+5ikKefc\nfyh9+NTlnPv5PNeJeFY4rgckXXLO/Wqui0NdzGyXpF3OuefNbLOkr8tf9/2fivdrIa1yTN8n3quF\nZmabnHNzpeyZv5T0M5J+VtKnnXN/Yma/Jel559xv57pQ1GyFY/rTkv67c+4z+a6uOIramb1H0ivO\nuVHn3LykJ+X/skaxmYr7ZxIlzrmnJVV+IPGApE+Wvv6kpH+c6aLQsBWOq0SmQWE55845554vfX1Z\n0kuSBsX7tbBWOKYDpR/zXi0w59xc6csO+UtrOknfK+m/lW7/pKQfzGFpqFOVY7pY+p73agxFLRwG\nJJ0u+35MS39Zo7icpC+Z2TNm9pN5LwaJ2umcm5D8yZak3pzXg+R8yMyeN7PfZRy1uMzsFknvkPRV\nSX28X4uv7Jh+rXQT79UCM7M2M3tO0jlJT0k6Iem8cy4qgMYk9ee1PsRXeUydc8+UfvTLpffqfzSz\n9TkusRCKWsxW+8SiePPSqPRdzrm/I+n75P/R/Z68FwRgVf9Z0j7n3Dvk/zFmhLGASuOon5b0M6Vu\nHv+eFlyVY8p7teCcc4vOuW+Xn564R9Kbq90t21WhEZXH1MzukvTzzrk3S3qnpB2SyAVaQ1GL2TFJ\ne8q+H5TfO4sCK3UA5JyblPRZ+b+s0RwmzKxPemNP16s5rwcJcM5NuqXghd+R/8cXBVIKjPm0pD9w\nzn2udDPv1wKrdkx5rzYP59xFSX8u6TslbS/lyEicCxdW2TG9v2wqZl7S74tz4TUVtZh9RtJtZrbX\nzDZIer+kz+e8JjTAzDaVPkmWmXVK+oeS/ibfVaEBphsnKD4v6aHS1x+U9LnKB6AQbjiupUIn8kPi\nPVtE/0XSi865Xy+7jfdrsS07prxXi83MeqLRcDPbKOkfSHpR0lckvbd0N96rBbLCMX05eq+amcnn\nFfBeXUMh04wlf2keSb8uX5D/nnPu3+e8JDTAzIbku7FOfhP8ExzTYjKzT0naLz8eMyHpgKQ/lfQn\nkm6W9C1J73XOnc9rjYhvheP6vfJ78hYlnZL0U9GnygifmX23pL+Q9IL8371O0i9KOirpj8X7tXBW\nOaY/Kt6rhWVmb5UPeGor/foj59y/KZ07PSmpS9Jzkh4sdfQQuFWO6Z9J6pH/4Ph5ST9dFhSFKgpb\nzAIAAAAAWldRx4wBAAAAAC2MYhYAAAAAUDgUswAAAACAwqGYBQAAAAAUDsUsAAAAAKBwKGYBAAAA\nAIVDMQsAQIDM7BNmdmfp619o5PEAADQjrjMLAEDgzOySc25L3usAACAkdGYBAMiRmW0ysy+Y2XNm\n9g0ze2/p9q+Y2d1m9u8kbTSzZ83sD0o/+zEz+1rptt8yM6vyvF8xs7sz/s8BACAzFLMAAOTrfkln\nnHPf7px7m6Qvlv/QOfcLkuacc3c75/5JaXT4fZK+yzl3t6RFST+W+aoBAMjZurwXAABAi3tB0q+U\nOrD/wzn3dJX7lHde3y3pbknPlDqyN0maSH+ZAACEhWIWAIAcOedeMbPvkPR9kn7ZzP6Xc+6XV3mI\nSfqkc+6XslkhAABhYswYAIAcmdluSVecc5+S9CvyXddK182svfT1n0n6YTPrLT2+y8z2ZLNaAADC\nQWcWAIB8vVV+zHhR0nVJP126vfxyA5+Q9IKZfb20b3ZY0pfNrK30mA9J+lbF83K5AgBAU+PSPAAA\nAACAwmHMGAAAAABQOBSzAAAAAIDCoZgFAAAAABQOxSwAAAAAoHAoZgEAAAAAhUMxCwAAAAAoHIpZ\nAAAAAEDhUMwCAAAAAArn/wNQ0bvjXp8SRgAAAABJRU5ErkJggg==\n", + "text/plain": [ + "<matplotlib.figure.Figure at 0x105cff7b8>" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "ent = np.zeros((L-1))\n", + "## measurment assumes right canonized MPS!!!\n", + "for j in range(L-1):\n", + " M[j], M[j+1], S = left_canonize(M[j],M[j+1], return_S = True)\n", + " ent[j] = -np.sum(S**2*np.log(S**2))\n", + "\n", + "\n", + "plt.figure()\n", + "plt.plot(ent,'-o')\n", + "plt.xlabel('site i')\n", + "plt.ylabel(r'$S(\\rho_i)$')\n", + "plt.xlim([0,L-2])\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "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.5.1" + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} diff --git a/exercises/ex13_solution/heisenberg_imaginary_time_evolution_solution.ipynb b/exercises/ex13_solution/heisenberg_imaginary_time_evolution_solution.ipynb new file mode 100644 index 0000000..5b3a583 --- /dev/null +++ b/exercises/ex13_solution/heisenberg_imaginary_time_evolution_solution.ipynb @@ -0,0 +1,507 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "import numpy as np\n", + "from numpy import transpose as tr, conjugate as co\n", + "from scipy.linalg import expm, svd\n", + "import math\n", + "import matplotlib.pyplot as plt\n", + "%matplotlib inline\n", + "plt.rcParams['figure.figsize'] = 16, 9" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Some useful functions" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": { + "collapsed": true + }, + "outputs": [], + "source": [ + "def dot(A,B):\n", + " \"\"\" Does the dot product like np.dot, but preserves the shapes also for singleton dimenstions \"\"\"\n", + " s1 = A.shape\n", + " s2 = B.shape\n", + " return np.dot(A,B).reshape((s1[0],s2[1]))\n", + "\n", + "def truncated_svd(thetas,chi):\n", + " \"\"\" Does an svd on two-site matrix thetas, and truncates to chi or the last nonzero singular value \"\"\"\n", + " U, S, Vh = svd(thetas,full_matrices = False)\n", + " # trunkieren\n", + " ind = np.where(np.isclose(np.cumsum(S[::-1])[::-1],0))[0]\n", + " if len(ind)>0:\n", + " chi_tr = min(chi,max(1,ind[0]))\n", + " else:\n", + " chi_tr = chi\n", + " S=S[:chi_tr]\n", + " w=1.-np.sum(S**2)\n", + " S/=math.sqrt(1-w)\n", + " U=U[:,:chi_tr]\n", + " Vh=Vh[:chi_tr,:]\n", + " return U, S, Vh\n", + "\n", + "def left_canonize(thetas,s,chi, return_S = False):\n", + " \"\"\" Splits up a two-site matrix thetas into two one-site matrices such that the left one is canonized \"\"\"\n", + " da, dg = thetas.shape[2], thetas.shape[3]\n", + " thetas = thetas.transpose((2,0,3,1)).reshape((da*s,dg*s)) # combine indizes\n", + " U, S, Vh = truncated_svd(thetas,chi)\n", + "\n", + " db = len(S)\n", + " U = U.reshape((da,s,db)).transpose((1,0,2))\n", + " Vh = Vh.reshape((db,dg,s)).transpose((2,0,1)) \n", + " for s2 in range(s):\n", + " Vh[s2] = dot(np.diag(S),Vh[s2])\n", + " if return_S:\n", + " return U, Vh, S\n", + " else:\n", + " return U, Vh\n", + "\n", + "def right_canonize(thetas,s,chi, return_S = False):\n", + " \"\"\" Splits up a two-site matrix thetas into two one-site matrices such that the right one is canonized \"\"\"\n", + " da, dg = thetas.shape[2], thetas.shape[3]\n", + " thetas = thetas.transpose((2,0,3,1)).reshape((da*s,dg*s)) # combine indizes\n", + " U, S, Vh = truncated_svd(thetas,chi)\n", + "\n", + " db = len(S)\n", + " U = U.reshape((da,s,db)).transpose((1,0,2))\n", + " Vh = Vh.reshape((db,dg,s)).transpose((2,0,1)) \n", + " for s1 in range(s):\n", + " U[s1] = dot(U[s1],np.diag(S))\n", + " if return_S:\n", + " return U, Vh, S\n", + " else:\n", + " return U, Vh\n", + "\n", + "\n", + "def apply_two_site_H(theta,H,s):\n", + " \"\"\" Applies a two-site operator on the two-site matrix theta \"\"\"\n", + " thetas = np.zeros_like(theta)\n", + " for s1p in range(s):\n", + " for s2p in range(s):\n", + " for s1 in range(s):\n", + " for s2 in range(s):\n", + " thetas[s1p,s2p] += H[s1p*s+s2p,s1*s+s2] * theta[s1,s2]\n", + "\n", + " return thetas\n", + "\n", + "def combine_two_matrices(M1,M2,s):\n", + " \"\"\" Combines the two neighbouring one-site matrices M1 and M2 into a two-site matrix theta \"\"\"\n", + " da, db, dg=M1.shape[1], M1.shape[2], M2.shape[2]\n", + " assert M1.shape[2] == M2.shape[1]\n", + "\n", + " theta = np.zeros((s,s,da,dg),dtype=complex)\n", + " for s1 in range(s):\n", + " for s2 in range(s):\n", + " theta[s1,s2]=dot(M1[s1],M2[s2])\n", + "\n", + " return theta" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Definition of system and initializiation" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "J = 1.\n", + "L = 20 # Length of chain\n", + "s = 2 # Local dimension of Hilbert space\n", + "dt = 0.05 # The imaginary timestep (should be 0.01 or smaller for good accuracy)\n", + "chi = 60 # The maximum matrix dimension, from which on the matrices are truncated\n", + "nmax = 10000 # Maximum number of iterations\n", + "nskip = 10 # Check energy convergence after nskip steps\n", + "\n", + "# Two-site Hamiltonian\n", + "H = np.array([[J/4,0,0,0],\n", + " [0,-J/4,J/2,0],\n", + " [0,J/2,-J/4,0],\n", + " [0,0,0,J/4]])\n", + "\n", + "# Imaginary time evolution operator\n", + "exp_beta_H=expm(-H*dt);\n", + "\n", + "# antiferromagnetic starting configuration\n", + "M = []\n", + "for i in range(L):\n", + " ar = np.zeros((2,1,1),dtype=complex)\n", + " ar[0,0,0] = i%2\n", + " ar[1,0,0] = (i+1)%2\n", + " M.append(ar)\n", + "\n", + "even = np.array(range(0,L-1,2))\n", + "odd = np.array(range(L+L%2-2,0,-2))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Imaginary time evolution" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "0 dE = 5.2022859299\n", + "10 dE = 2.41237932615\n", + "20 dE = 0.616048162871\n", + "30 dE = 0.19945192398\n", + "40 dE = 0.0856110936378\n", + "50 dE = 0.0446638884127\n", + "60 dE = 0.0266683370958\n", + "70 dE = 0.0176509913134\n", + "80 dE = 0.0127047223848\n", + "90 dE = 0.00977408061919\n", + "100 dE = 0.00789488605981\n", + "110 dE = 0.00658572945765\n", + "120 dE = 0.00559930606822\n", + "130 dE = 0.00480731763688\n", + "140 dE = 0.00414298945949\n", + "150 dE = 0.0035712463933\n", + "160 dE = 0.00307305023898\n", + "170 dE = 0.00263722905494\n", + "180 dE = 0.00225631290084\n", + "190 dE = 0.00192449985237\n", + "200 dE = 0.00163673118822\n", + "210 dE = 0.00138832203542\n", + "220 dE = 0.0011748521653\n", + "230 dE = 0.000992164333084\n", + "240 dE = 0.000836395474847\n", + "250 dE = 0.000704007231025\n", + "260 dE = 0.000591803752473\n", + "270 dE = 0.000496934572579\n", + "280 dE = 0.000416884722746\n", + "290 dE = 0.000349455574773\n", + "300 dE = 0.000292739848431\n", + "310 dE = 0.000245093656\n", + "320 dE = 0.000205107728403\n", + "330 dE = 0.000171579325436\n", + "340 dE = 0.000143485782885\n", + "350 dE = 0.000119960236445\n", + "360 dE = 0.000100269773123\n", + "370 dE = 8.37960556659e-05\n", + "380 dE = 7.00183361975e-05\n", + "390 dE = 5.84987014349e-05\n", + "400 dE = 4.88693502501e-05\n", + "410 dE = 4.08216810968e-05\n", + "420 dE = 3.40969704631e-05\n", + "430 dE = 2.84784768727e-05\n", + "440 dE = 2.37847443216e-05\n", + "450 dE = 1.98639157318e-05\n", + "460 dE = 1.65889561448e-05\n", + "470 dE = 1.38536333161e-05\n", + "480 dE = 1.15691357934e-05\n", + "490 dE = 9.661232097e-06\n", + "500 dE = 8.06788848529e-06\n", + "510 dE = 6.73727353018e-06\n", + "520 dE = 5.62608618715e-06\n", + "530 dE = 4.69815473814e-06\n", + "540 dE = 3.92326507992e-06\n", + "550 dE = 3.27618108287e-06\n", + "560 dE = 2.73582595511e-06\n", + "570 dE = 2.28459734686e-06\n", + "580 dE = 1.90779520892e-06\n", + "590 dE = 1.59314354775e-06\n", + "600 dE = 1.33039094941e-06\n", + "610 dE = 1.11097674527e-06\n", + "620 dE = 9.27752324742e-07\n", + "630 dE = 7.747481785e-07\n", + "640 dE = 6.46979531282e-07\n", + "650 dE = 5.40283849659e-07\n", + "660 dE = 4.51185265149e-07\n", + "670 dE = 3.76781240874e-07\n", + "680 dE = 3.1464805339e-07\n", + "690 dE = 2.62761810532e-07\n", + "700 dE = 2.19432413573e-07\n", + "710 dE = 1.83248596741e-07\n", + "720 dE = 1.53031841421e-07\n", + "730 dE = 1.27798031713e-07\n", + "740 dE = 1.06725389415e-07\n", + "750 dE = 8.91276528137e-08\n", + "760 dE = 7.44317638635e-08\n", + "770 dE = 6.21591595973e-08\n", + "780 dE = 5.19102449914e-08\n", + "790 dE = 4.33512816755e-08\n", + "800 dE = 3.62035965651e-08\n", + "810 dE = 3.02344727032e-08\n", + "820 dE = 2.52495695463e-08\n", + "830 dE = 2.10865902517e-08\n", + "840 dE = 1.76100023452e-08\n", + "850 dE = 1.47066288037e-08\n", + "860 dE = 1.22819674431e-08\n", + "870 dE = 1.02570609783e-08\n", + "880 dE = 8.56601900523e-09\n", + "Converged after 880 steps!\n", + "Ground-state energy: -8.68010558621\n" + ] + } + ], + "source": [ + "enold = 0. # for checking energy convergence\n", + "\n", + "for n in range(nmax):\n", + " # ++++ Trotter scheme: first even bonds, then odd bonds\n", + " # === apply exp_beta_H on all even bonds ===\n", + " for j in even:\n", + " # go from left to right\n", + " theta = combine_two_matrices(M[j],M[j+1],s)\n", + " thetas = apply_two_site_H(theta,exp_beta_H,s) \n", + " M[j], M[j+1] = left_canonize(thetas,s,chi)\n", + "\n", + " # advance left canonization by a further step\n", + " if j < L-2:\n", + " theta = combine_two_matrices(M[j+1],M[j+2],s) \n", + " M[j+1], M[j+2] = left_canonize(theta,s,chi)\n", + " # === apply exp_beta_H on all even bonds ===\n", + "\n", + " # === renormalize the state on the last site ===\n", + " da = M[L-1].shape[1]\n", + " theta = M[L-1].reshape((s*da,1))\n", + " U, _, _ = truncated_svd(theta,chi) # throw away norm\n", + " M[L-1] = U.reshape((s,da,1))\n", + " # === renormalize the state on the last site ===\n", + "\n", + " if L%2 == 0:\n", + " # the right-most matrix has to be right canonized in this case\n", + " theta = combine_two_matrices(M[L-2],M[L-1],s)\n", + " M[L-2], M[L-1] = right_canonize(theta,s,chi)\n", + "\n", + " # === apply exp_beta_H on all odd bonds ===\n", + " for j in odd:\n", + " # go from right to left\n", + " theta = combine_two_matrices(M[j-1],M[j],s)\n", + " thetas = apply_two_site_H(theta,exp_beta_H,s)\n", + " M[j-1], M[j] = right_canonize(thetas,s,chi)\n", + "\n", + " # advance right canonization by a further step\n", + " if j>1:\n", + " theta = combine_two_matrices(M[j-2],M[j-1],s)\n", + " M[j-2], M[j-1] = right_canonize(theta,s,chi)\n", + " # === apply exp_beta_H on all odd bonds ===\n", + "\n", + " # === renormalize the state on the first site ===\n", + " db = M[0].shape[2]\n", + " theta = M[0].reshape((1,s*db))\n", + " _, _, Vh = truncated_svd(theta,chi) # throw away norm\n", + " M[0] = Vh.reshape((s,1,db))\n", + " # === renormalize the state on the first site ===\n", + "\n", + " # ++++ Calculate energy and check convergence ++++\n", + " if n%nskip == 0:\n", + " # Measure Energy via sum of two-site Operators H\n", + " en = 0.\n", + " for j in range(L-1): \n", + " theta = combine_two_matrices(M[j],M[j+1],s)\n", + " thetas = apply_two_site_H(theta,H,s) \n", + " for s1 in range(s):\n", + " for s2 in range(s):\n", + " en += np.trace(dot(thetas[s1,s2],co(tr(theta[s1,s2])))).real\n", + " M[j], M[j+1] = left_canonize(theta,s,chi) # Proceed with left canonization\n", + "\n", + " # make right canonized MPS for next tDMRG step\n", + " for j in range(L-1,0,-1):\n", + " theta = combine_two_matrices(M[j-1],M[j],s)\n", + " M[j-1], M[j] = right_canonize(theta,s,chi)\n", + "\n", + " print(n,'dE = ',abs(en-enold))\n", + " if abs(en-enold) < 1e-8:\n", + " print(\"Converged after \"+str(n)+\" steps!\")\n", + " print(\"Ground-state energy: \" + str(en))\n", + " break\n", + " enold = en\n", + " # === check convergence ===\n", + "\n", + "if n == nmax:\n", + " print(\"Maximum iterations reached! Convergence criterium not satisfied!\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Implement the Heisenberg MPO below\n", + "\n", + "You can build it as a list of four-dimensional arrays with indexing $\\hat W^{[i]} = W_{\\sigma'_i\\, \\sigma_i\\, b_{l-1}\\, b_l}$. It can be represented as \"matrices\" with indices $(b_{l-1}, b_l)$ which have single-site operators as entries having $(\\sigma'_l,\\sigma_l)$ indices. For the Heiseberg model \n", + " \\begin{equation}\n", + " H = \\sum_{i=1}^{L-1} \\left( \\frac{J}{2} S_i^+ S_{i+1}^- + \\frac{J}{2} S_i^- S_{i+1}^+ + J S_i^z S_{i+1}^z\\right) ,\n", + " \\end{equation}\n", + " the matrices are\n", + " \\begin{equation}\n", + " \\hat W^{[i]} = \\begin{bmatrix}\n", + " \\hat I & 0 & 0 & 0 & 0 \\\\\n", + " S^+_i & 0 & 0 & 0 & 0 \\\\\n", + " S^-_i & 0 & 0 & 0 & 0 \\\\\n", + " S^z_i & 0 & 0 & 0 & 0 \\\\\n", + " 0 & (J/2) S^-_i & (J/2) S^+_i & J S^z_i & \\hat I\n", + " \\end{bmatrix},\n", + " \\end{equation}\n", + " and for the first and last site\n", + " \\begin{equation}\n", + " \\hat W^{[1]} = \\begin{bmatrix}\n", + " 0 & (J/2)S^-_0 & (J/2) S^+_0 & J S^z_0 & \\hat I\n", + " \\end{bmatrix}, \\ \\ \\hat W^{[L]} = \\begin{bmatrix}\n", + " \\hat I \\\\ S^+_L \\\\ S^-_L \\\\ S^z_L \\\\ 0\n", + " \\end{bmatrix}.\n", + " \\end{equation}" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": { + "collapsed": false + }, + "outputs": [], + "source": [ + "## Generate MPO\n", + "W = []\n", + "\n", + "## some operators\n", + "Z = np.zeros((2,2,1,1))\n", + "\n", + "Id = np.zeros((2,2,1,1))\n", + "Id[0,0,0,0] = 1.\n", + "Id[1,1,0,0] = 1.\n", + "\n", + "Sz = np.zeros((2,2,1,1))\n", + "Sz[0,0,0,0] = +0.5\n", + "Sz[1,1,0,0] = -0.5\n", + "\n", + "Sp = np.zeros((2,2,1,1))\n", + "Sp[0,1,0,0] = 1.\n", + "\n", + "Sm = np.zeros((2,2,1,1))\n", + "Sm[1,0,0,0] = 1.\n", + "\n", + "W.append(np.concatenate((Z,J/2*Sm,J/2*Sp,J*Sz,Id),axis=3))\n", + "for i in range(1,L-1):\n", + " Wi = np.concatenate((\n", + " np.concatenate((Id,Z,Z,Z,Z),axis=3),\n", + " np.concatenate((Sp,Z,Z,Z,Z),axis=3),\n", + " np.concatenate((Sm,Z,Z,Z,Z),axis=3),\n", + " np.concatenate((Sz,Z,Z,Z,Z),axis=3),\n", + " np.concatenate((Z,J/2*Sm,J/2*Sp,J*Sz,Id),axis=3),\n", + " ),axis=2)\n", + " W.append(Wi)\n", + "W.append(np.concatenate((Id,Sp,Sm,Sz,Z),axis=2))" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "collapsed": true + }, + "source": [ + "### Build iteratively the L- or R-expression and evaluate the energy of the MPS using above MPO\n", + "\n", + "You can use a three-dimensional array for $R$ with the indexing $R_{a'_l\\, b_l\\, a_l}$." + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": { + "collapsed": false + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Energy from R-expression: (-8.680105586208278+0j)\n", + "Energy from L-expression: (-8.680105586208281+0j)\n" + ] + } + ], + "source": [ + "def add_site_to_R(Rexp,M,W):\n", + " Rexp = np.tensordot(M.conj(),Rexp,axes=((2),(0)))\n", + " Rexp = np.tensordot(W,Rexp,axes=((0,3),(0,2)))\n", + " Rexp = np.tensordot(M,Rexp,axes=((0,2),(0,3)))\n", + " Rexp = np.transpose(Rexp,(2,1,0))\n", + " return Rexp\n", + "\n", + "def add_site_to_L(Lexp,M,W):\n", + " Lexp = np.tensordot(M.conj(),Lexp,axes=((1),(0)))\n", + " Lexp = np.tensordot(W,Lexp,axes=((0,2),(0,2)))\n", + " Lexp = np.tensordot(M,Lexp,axes=((0,1),(0,3)))\n", + " Lexp = np.transpose(Lexp,(2,1,0))\n", + " return Lexp\n", + "\n", + "\n", + "## Calculate energy via iteratively constructing the R-expression\n", + "Rexp = np.ones((1,1,1))\n", + "for i in range(1,L+1):\n", + " Rexp = add_site_to_R(Rexp,M[-i],W[-i])\n", + "print('Energy from R-expression: ',np.squeeze(Rexp)) # should be total energy\n", + "\n", + "## Calculate energy via iteratively constructing the L-expression\n", + "Lexp = np.ones((1,1,1))\n", + "for i in range(L):\n", + " Lexp = add_site_to_L(Lexp,M[i],W[i])\n", + "print('Energy from L-expression: ',np.squeeze(Lexp)) # should be total energy" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "collapsed": true + }, + "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.5.1" + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} -- GitLab