11_toric_code.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Toric Code\n",
    "\n",
    "This notebook shows how to implement Kitaev's Toric Code model on an infinite cylinder geometry. \n",
    "The Hamiltonian of the Toric Code is:\n",
    "\n",
    "$$\\mathcal{H}=-\\sum_p B_p - \\sum_s A_s$$ \n",
    "\n",
    "where $B_p = \\prod_{i \\in p} \\sigma^z_i$ denotes the product of $\\sigma^z$ around an elementary plaquette and $A_s = \\prod_{i \\in s} \\sigma^x_i$ denotes the product of $\\sigma^x$ on the four links sticking out of a site of the lattice. The $B_p$ and $A_s$ operators commute with each other and have eigenvalues $\\pm 1$, hence the ground state will have eigenvalue $+1$ for each of them.\n",
    "\n",
    "On an infinite cylinder, the model exhibits a four-fold degeneracy of the ground state, characterized by the \n",
    "different values $\\pm 1$ of the incontractible Wilson and t'Hooft loops.\n",
    "These loop operators are defined by the products\n",
    "\n",
    "$$W=\\prod_{\\text{vert links}} \\sigma^z\\qquad \\qquad H=\\prod_{\\text{hor links}} \\sigma^x$$\n",
    "\n",
    "around the cylinder.\n",
    "\n",
    "TeNPy already implements this hamiltonian in `tenpy.models.toric_code.ToricCode`.\n",
    "This notebook demonstrates how to extend the existing model by additinoal terms in the Hamiltonian\n",
    "and convenient functions of the loop operators for simplified measurements.\n",
    "Then we run iDMRG to obtain the four different ground states and check that they are indeed orthogonal and degenerate.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy\n",
    "import matplotlib.pyplot as plt\n",
    "np.set_printoptions(precision=5, suppress=True, linewidth=120)\n",
    "plt.rcParams['figure.dpi'] = 150"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import tenpy\n",
    "import tenpy.linalg.np_conserved as npc\n",
    "from tenpy.algorithms.dmrg import TwoSiteDMRGEngine\n",
    "from tenpy.networks.mps import MPS\n",
    "from tenpy.networks.terms import TermList\n",
    "from tenpy.models.toric_code import ToricCode, DualSquare\n",
    "from tenpy.models.lattice import Square\n",
    "\n",
    "tenpy.tools.misc.setup_logging(to_stdout=\"WARNING\")  # don't show info text"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's plot the DualSquare lattice first to get the geometry clear. The unit cell consists of two sites, the blue one on the vertical links ($u=0$), and the orange one on the horizontal links ($u=1$).\n",
    "We also plot lines illustrating where the Wilson and t'Hooft loops go."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAx0AAAJ2CAYAAAA63e++AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAXEQAAFxEByibzPwAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzde3hU1d33/8/K5CQJATmpGDAYq5DImUSrULUq9a7lLrVA64F6QHvAVqvVPsVfb9Hqr09blVofq3dbD6CoyKHaW22lpQ/QAEoQPNBEUZFwMCoYRELCkNN+/pgkEDJDZkLW3ivJ+3VduSaZvfbMd11fhuQzs/bexvM8AQAAAIAtSUEXAAAAAKBrI3QAAAAAsIrQAQAAAMAqQgcAAAAAqwgdAAAAAKwidAAAAACwitABAAAAwCpCBwAAAACrCB0AAAAArCJ0AAAAALCK0AEAAADAKkIHAAAAAKuSgy7AT8aYjyX1kLQ96FoAAAAABw2SVO153vEd+aDG87yOfDynGWP2pqWl9czNzQ26FAAAAMA5mzdv1oEDByo9z8vqyMftVp90SNqem5ubV1JSEnQdAAAAgHPy8/NVWlra4auCOKYDAAAAgFWEDgAAAABWEToAAAAAWEXoAAAAAGAVoQMAAACAVYQOAAAAAFYROgAAAABY1d2u0wEAANCpeZ6n7nRxZ7SPMUbGmKDLaEboAAAAcFx9fb0qKipUWVmpmpqaoMtBJ5GamqqePXuqb9++CoVCgdZC6AAAAHBYfX29tm3bpnA4HHQp6GRqampUUVGhqqoqDR48ONDgQegAAABwWEVFhcLhsEKhkI477jhlZGQoKYnDcnFkDQ0Nqqqq0ieffKJwOKyKigoNGDAgsHoIHQAAAA6rrKyUJB133HHq1atXwNWgs0hKSmr+91JeXq7KyspAQwcxGQAAwFGe5zUfw5GRkRFwNeiMmv7d1NTUBHoCAkIHAACAow79I5ElVWiPQ//dEDoAAAAAdFmEDgAAAABWEToAAAAAWEXoAAAAAGAVoQMAAACdxoYNG2SM0UknnRR1e11dnTIyMmSM0Y033hh1zEsvvSRjjMaNGydJmjt3rowxuuOOO1qMu+OOO2SM0dy5cztyCt0SoQMAAACdxsiRI5WVlaVt27Zp27ZtrbavX79e1dXVkqSioqKoj7Fq1SpJ0oQJE+wViha4OCAAAEBnt25dfOMGD5aOO+7onmvXLqmsLL6xY8dKHXyq31AopC9+8YtaunSpioqKdPnll7fY3hQoRo4cqTfffFN79+5VVlZW1DFNoeMb3/iGzjzzTPXr169Da8VBhA4AAIDOrrAwvnEPPihdf/3RPdfzz0vf/W58Y2tqOjx0SJGwsHTpUq1atSpq6DjmmGP0wx/+UNddd51eeeUVfeUrXzmkpBq99tprkqTx48dLknr16sXV3i1jeRUAAAA6laZPKKItn1qzZo0KCgp03nnnRR2zbt06hcNhnXbaaRowYICk2Md0xFJVVaVf//rXGjVqlHr37q3MzEzl5uZq6tSpWrp0aavx27dv1/e+9z2ddNJJSktL04ABA3TJJZdoXZRPqMrKymSM0bnnnqv9+/frZz/7WfN+p5xyin79618HepG/9uKTDgAAAHQqhYWFSktLU2lpqXbv3q0+ffpIkjZt2qSdO3fq2muvVW5uro4//vhWoeNoj+eor6/XxIkTtWbNGmVnZ+vcc89VamqqduzYoRdffFEZGRktPlnZuHGjvvzlL+vTTz/V0KFDdckll2jbtm167rnn9MILL+jpp5/W1KlTWz1PTU2NJk6cqJKSEhUWFmrYsGFauXKlfvazn6myslJ33313u+oPCp90AAAAoFNJT0/XuHHj5HmeVq9e3Xx/U6A4++yzJUlnnXWWiouLVVNT02pMe0NHUVGR1qxZo69//evaunWrnn/+eS1cuFBr1qzRxx9/rB/96EfNYz3P0+WXX65PP/1Us2bNUmlpqZ555hmtXr1aixYtUkNDg2bMmKFPPvmk1fO88sorMsbo3Xff1csvv6yXX35ZRUVFSk5O1m9/+1vt27evXfUHhdABAACATifaEqtVq1bJGKOzzjpLUiR8hMPh5mM4PM/TmjVrWuyfqJ07d0qSzj33XCUddrxKr169NHbs2OafV6xYoY0bN2rIkCG66667ZIxp3jZlyhRNnjxZlZWVevzxx1s9T1JSkh555JEWB7ePGzdO//Ef/6Hq6urmOXUWhA4AAAB0Ok2hoemTi6bvTz/9dPXu3VvSwU88msY0Lcc68cQTNWTIkHY976hRo5SUlKR77rlHCxYsUGVlZcyxTYHoW9/6lkKhUKvt06dPbzHuUDk5OTr11FNb3d9030cffdSu+oNC6AAAAECnc/bZZyspKUmvvfaa9u/fr507d+r9999vDhqSNGbMGB1zzDHNf9R3xPU5Tj31VN1zzz3atWuXLr30Uh177LEaNWqUbrnlFpWUlLQYW15eLikSIKJpur9p3KGys7Oj7pOZmSlJOnDgQDtnEAxCBwAAADqdXr16afjw4aqtrdXatWubg0XTaXAlKSUlRQUFBVq9erU8z+uwiwLefPPN2rx5sx544AF99atf1datW3XfffdpxIgR+v3vf99q/KHLqqKJtr2tfTobQgcAAAA6pUOXWB1+EHmTs88+W5999plKSko69ErkgwYN0o9+9CP9z//8j3bt2qUnn3xSSUlJuvnmm7Vnzx5J0sCBAyVJW7ZsifoYW7dulSSdcMIJR12P6wgdAAAA6JQOPZh81apVOvHEE1stZWoKIc8++6zKysrUu3dvnX766R1aR3Jysq644goVFBSopqZG7777bov6nn32WdXX17fab/78+S3GdWVcpwMAAKCzKy6Ob9zgwUf/XJMnS6NGxTc2ysHTHanpj/U1a9YoHA7rm9/8ZqsxZ511lowxevDBByVFll8dzdKl5cuXq76+Xl/+8pdbnL1q69atevvtt2WMaT4e49xzz9Xw4cO1ceNG3X777br77rubn/v555/Xn//8Z2VmZuqqq65qdz2dBaEDAACgsyso8O+5+vePfDnghBNOUG5urjZv3iyp9dIqSTr22GM1bNgwlZaWSjr6TxXefPNN3XTTTerfv7/Gjh2rvn37ateuXfrXv/6lcDisH//4x83Lqowxeuqpp3Teeefpl7/8pZ577jmNGjVK27Zt0+rVq5WcnKzHHntMxx9//FHV1BkQOgDAZe+9J916a3xj77lH+sIX7NZzFLrQVAA4ZMKECc2h49CDyA919tlnN4eOWGPi9bWvfU0VFRVavny53nzzTVVUVKh///6aMGGCZs6cqcmTJ7cYP3z4cG3YsEF33323Xn75ZS1evFi9evXS5MmTNWvWLBUWFh5VPZ2F8Twv6Bp8Y4wpycvLyzv8dGYA4Kx166R4fyEVF/v7bmeCutBUAN80NDRo06ZNkqTTTjut1cXogLYk+m8oPz9fpaWlpZ7n5XdkHfzLBQAAAGAVoQMAAACAVYQOAAAAAFYROgAAAABYRegAAAAAYBWhAwAAAIBVhA4AAAAAVhE6AAAAAFjFFckBwGVZWdIFF8Q/1mFdaCoAgAQROgDAZaedJv3jH0FX0SG60FQAAAlieRUAAAAAq5wNHcaYPsaYncYYzxjzTtD1AAAAAGgfZ0OHpDmS+gVdBAAAAICj42ToMMacL+lKSX8KuhYAAAB0Ljk5OcrJyWlxX1lZmYwxre4/3FVXXSVjjObOnWutviOpr6/X7bffrtzcXKWmpsoYo6uuuiqQWjqScweSG2OOkfTfkkol3Svpu8FWBFvq6uq0Y8eOqNuys7OVnOzcP89ugb64ib64ib4AwSgrK9OQIUN0zjnnaMWKFUGX06F+97vf6a677tLAgQN1ySWXKD09XePHj5cUCVNbt26V53kBV5k4F/83nC0pV9K5kmqDLQU27dixQ0OGDIm6bcuWLW2+EwE76Iub6Iub6AtcsW5dfOMGD5aOO+7onmvXLqmsLL6xY8dKSU6uq3HX888/L0kqKirSySefHHA1Hcep0GGMGSHpJ5Ie9zzvX8aYnGArAgAAcF9hYXzjHnxQuv76o3uu55+XvhvnOpSaGkJHopo+Pe1KgUNy6JgOY0ySIsdw7JH006N8rJJoX4p8ggIAAIBO6o477mj+hHHlypUyxjR/2T72oaKiQrfeequ+8IUvKD09XX369NFFF12kv//97zH3KS0t1eWXX64TTjhBqampOvHEE/Wd73xHmzZtajGu6ViSLVu2SFKLec2dO1fGGG3durXVts7yiapLn3T8SFKhpKs9z6sIuhgAAAC4Z9SoUfrmN7+pJUuW6LjjjtNFF13UvK3p2AcbPvzwQ33pS1/SBx98oMGDB2vy5MnatWuXli1bpqVLl2rOnDm66aabWuzzz3/+U5MmTdL+/fs1ZswYnXvuuXrnnXf05JNP6rnnntNf//pXTZgwoUXtixcvVlVVla688srmxznllFN05ZVXRt3Wr1/nONmrE6HDGDNI0t2SVnqeN/doH8/zvPwYz1MiKe9oHx8AAADBmDx5skaNGqUlS5Zo6NChUc8yVRbvQScJ+P73v68PPvhA06dP16OPPqqUlBRJ0qpVq/SVr3xFt956q84//3yNGDFCklRVVaXLL79c+/fv18MPP6zvf//7zY/129/+VjfffLMuu+wyvf/++0pLS9O1116ra6+9VitWrFBVVVWreY0fPz7mts7AleVVD0lKlfSDoAsBAABA17V169YWy5MO/5o3b16rfT744AO9+OKLysrK0gMPPNAcOKRIGPj+97+v+vp6PfTQQ833L1y4UJ988okmTJjQInBI0k033aSxY8dqx44deu655+xN1iFOfNIh6WuKHMvxsDHm0PvTG28HG2NWNI31PG+fj7UBAACgi8jIyNCUKVNibl+1apU2b97c6j5J+upXv6revXu32mf69OmaM2eOioqKmu9r+v7yyy+P+jxXXHGF1q9fr6KiIn37299OeB6djSuhQ5J6SzonxrZjDtnmUs0AAADoRPr163fE5UlXXXVVq9BRXl4uSTEP2m66v2lce/fpypxYXuV5non2Janp5OebDrl/T5C1AgAAoHs6bEVOq/ujbY+1T7zbuwonQgcAAADgqoEDB0pS8+lsD9d04PoJJ5wQ9z5Np789dJ+ujNABAC6rrZV27VJSRYX6SVG/kioqIpcIrq0NstI2NU4lri/HpwIgYKmpqZKkuro6X56v6XS2L730kvbsab3oZv78+ZLUfPrbQ79/6qmnoj5m0/2H7tMWv+fdkTg+AgBc9sYbUmGhBkvaFWvMuHGR2+JiqaDAp8IS1ziVuDg+FcA5xcXxjRs8+Oifa/JkadSo+MaGQkf/fNH069dPKSkp2rx5s+rr6xWy9USNTj75ZF188cV66aWXdOONN+qRRx5pPoPVK6+8oocfflihUEgzZ85s3mfatGmaNWuWioqK9Mc//lHfPeQy7g888IDWrVun7OxsfeMb34i7joEDB+q9997Tpk2blJ8f9QoRznI6dHieVyapeyx0AwAAaCc/Q3r//pGvIKWmpuqiiy7SCy+8oJEjR2rMmDFKTU3V2WefrauvvtrKc/7hD3/QhAkT9MQTT2jlypX64he/qF27dmnFihWqr6/Xfffd13yNDilylqynnnpKkyZN0ve+9z398Y9/1Kmnnqp33nlHr7/+ujIyMvT0008rLS0t7hr+8z//UytXrtT555+v8847TxkZGerXr59+9atf2Zhyh2J5FQAAADqdRx55RNOnT1dFRYWefvppPfroo1q5cqW15zvxxBO1bt06/eQnP1FycrL+/Oc/a/369Tr//PO1dOlS3Xzzza32Of/887Vu3Tpdeuml2rFjhxYvXqyPP/64+XS5iSytkqQbbrhBP//5z5WZmaklS5bo0Ucf1YIFCzpqilY5/UkHAAAAEM2AAQP0xBNPxD0+JydHnue1OW7u3LkxT6nbt29f3Xvvvbr33nvjft78/Hw9/fTTcY8/0tXUk5OTddddd+muu+6K+/FcwScdAAAAAKwidAAAAACwitABAAAAwCpCBwAAAACrCB0AAAAArCJ0AAAAALCKU+YCgMtyc6WFC7Vz505d/8MfRh3y+wcf1IABAyJjHdY4lbjHAgC6DkIHALisTx9p6lRVl5VpcYzQcc/FF0s5Of7W1Q6NUwGQAGNM8/cNDQ1KSmKRChLT0NDQ/P2h/578xr9cAAAARxljlJqaKkmqqqoKuBp0Rk3/blJTUwMNHXzSAQAA4LCePXuqoqJCn3zyiSQpIyODTzzQpoaGBlVVVTX/u+nZs2eg9RA6AAAAHNa3b19VVVUpHA6rvLw86HLQCaWnp6tv376B1kDoAAAAcFgoFNLgwYNVUVGhyspK1dTUBF0SOonU1FT17NlTffv2VSgUCrQWQgcAAIDjQqGQBgwYoAEDBsjzPHmeF3RJcJwxJtBjOA5H6AAAAOhEXPtjEogHRyEBAAAAsIrQAQAAAMAqQgcAAAAAqwgdAAAAAKwidAAAAACwitABAAAAwCpCBwAAAACrCB0AAAAArCJ0AAAAALCKK5IDgMvKy6XHH1evzz7TbTGG9HrwQenYY6Wrr5YGDvS1vEQ0TiUujk8FAJAgQgcAuOzDD6Wf/1zHSvr/Y425777I7cSJTv+l3jiVuDg+FQBAglheBQAAAMAqQgcAAAAAqwgdAAAAAKwidAAAAACwitABAAAAwCpCBwAAAACrCB0AAAAArOI6HQhUUo/eyhw5UemDTldSag811FQrvG2jdlfXKSfo4rox+uKmXT1669mRE7V20Onal9pDmTXVOnPbRk176x/qX70n6PK6LV4vANA2QgcCEa6t170ry5U983GZUEqLbccMGaNp8zdp6rj9mj0pT+kpoYCq7H7oi4NCIYX79NMdZ03X4qFfUt1hfSkaMkb3j79cUzat1GwlKT2gMuMRCkm9esU/1nW8XgAgfsbzvKBr8I0xpiQvLy+vpKQk6FK6tXBtva58rFhrt+xuc+wZQ/po3jWF/ML2AX1xE31xE30B0FXl5+ertLS01PO8/I58XI7pgO/ufKEkrl/UkrR2y27d+UKp5Yog0RdX0Rc30RcASAyhA77aWRnW4vU7Etpn8frt2lV5wFJFkOiLq+iLm+gLACSO0AFfLVy3XbX1iS3pq633tPC17ZYqgkRfXEVf3ERfACBxhA74Kt7lCId79YOKDq4Eh6IvbqIvbqIvAJA4Qgd8te9Ana/7IT70xU30xU30BQASR+iArzLT2neW5vbuh/jQFzfRFzfRFwBIHKEDvjpjSJ927XfmyX07uBIcir64ib64ib4AQOIIHfDVtIJBSgmZhPZJCRlNGzfIUkWQ6Iur6Iub6AsAJI7QAV8N6JmuKWOzE9pnythB6t8zzVJFkOiLq+iLm+gLACSO0AHfzZ6UH/fyhDOG9NHsSXmWK4JEX1xFX9xEXwAgMYQO+C49JaR51xRq0rBj5dXXRh2TnCRdWjhY864pVHpKyOcKuyf64ib64ib6AgCJMZ6X2AWOOjNjTEleXl5eSUlJ0KVAUllZmXLzRytzxIVKHzxcSak91FBTrfC2jVr+6C81Ju+UoEvsluiLm+iLm+gLgK4mPz9fpaWlpZ7n5Xfk43L+PgSqoXqP9r66SHtfXdTi/j49fhNQRZDoi6voi5voCwC0jeVVAAAAAKwidAAAAACwitABAAAAwCqO6QAAl735pvT1ryu7rk5bYgzJHj9eSk6W/vIXaeRIX8tLRONU4uL4VAAACSJ0AIDLamqkrVuVLCkn1pgPPzw41mGNU4l7LACg62B5FQAAAACrCB0AAAAArCJ0AAAAALCK0AEAAADAKkIHAAAAAKsIHQAAAACsInQAAAAAsIrQAQAAAMAqLg4IAC47/njpllv0+eef649/+lPUId+97jr16tUrMtZhjVOJeywAoOsgdACAywYNku65R5+VlemnMULH1NtuU6+cHH/raofGqQAAuiGWVwEAAACwitABAAAAwCpCBwAAAACrCB0AAAAArCJ0AAAAALCK0AEAAADAKkIHAAAAAKsIHQAAAACsInQAAAAAsIrQAQAAAMAqQgcAAAAAqwgdAAAAAKwidAAAAACwitABAAAAwCpCBwAAAACrkoMuAABwBJ9/LhUXK/3jj3VBjCHpRUXSe+9JhYVSr16+lpeIxqnExfGpAAASROgAAJe9+640caKOl/SPWGO+853IbXGxVFDgU2GJa5xKXByfCgAgQSyvAgAAAGAVoQMAAACAVYQOAAAAAFY5FTqMMTcbY/5sjHnPGPO5MeaAMWarMWaeMSY/6PoAAAAAJM61A8lvk5Qh6S1JGxvvy5f0HUnfNsZM9jzvb0EVB7SpNiyt/e/I92f+QEpOC7YeRNAXJ6WFwrrhjEhffrf2B6qppy9O4PUCwALXQsfXJa33PC986J3GmB9IekjSI8aYwZ7n1QdSHdCW15+Uls2OfJ+aIRVeF2w9iKAvTpox5kn95sJIX/bVZOjh1+iLE3i9ALDAqeVVnuetPjxwNN7/sKT3JQ2UdJrvhQHxqA1LRXMO/lw0R6o7EFw9iKAvTkoLhTVr/MG+3DZhjlJD9CVwvF4AWOLaJx1H0vTpRk2gVaDDHDgQ+xdZWVmZf4V0kJ5vL1DfyvKDd1SWq2LZ/aoc+q3gimqHzZs3t2ubqzp7X1LLyzUwzrHl5eWqcfi1U16eKjXOZsaYJ5WddbAv2VnlmjH6ieZPO8rLy1VW5v5/97xeOp/s7GwlJ3emP3+ArsF4nhd0DW0yxnxH0jxJ70oa5nleQzsfpyQvLy+vpKSkQ+tD+yxatEjTpk0LuowOkRaS3r8hU9lZLT883LG3QbkP7FMNCwID0RX6Mk7SujjHFkh6zWItRy8ym7RQWO/fMLpF6JCkHXsHKveBNxqP7XB/Nl1NV3i9xGPLli3KyckJugzAWfn5+SotLS31PK9DT+Lk1PKqJsaYW40xc40xi4wx/1YkcJRLuqy9gQOwacaYlFa/qCUpOytJM0anBFARpK7Rl42SvhDn18YYj+GOyGxmjMltFTikpk87TlZnmU1X0xVeLwDc5erni1+RdP4hP2+XNN3zvPXx7GyMifVRRu7RFgYcLi0kzRof++wut01I06Ov13aZdwk7i67SlwOKHNDWNRxQWuh9zRqfqVjved02Ya8efb3c+b50NV3l9QLAXU5+0uF53gWe5xlJx0r6kqRNklYYY/6/YCsDWov17mAT3iUMBn1xE31xE30BYJuToaOJ53l7PM8rkvRVSesl3WWMKYhjv/xoX5I631F9cFpqG+8ONpk1Pk2pIR8KgiT64ir64ib6AsAPToeOJp7n1Up6VpKRNCngcoBmM0Yf+d3BJoN6Jeka3iX0DX1xE31xE30B4AdXj+mI5tPG2/6BVoEOM3z48Jjbli9f7v7ZReprlL3kYql6Z1zD/8/UHP2vBS9KoVTLhR2dzZs364ILLoi6bdmyZcrNdfzQKPriJvripi7alyPJzs4OugSgW+pMoeOcxluWSHUR6enpMbfl5OS4HzqK/xT3L2pJSq7+RDm7/yUVXGuxKLtyc3Ppi4Poi5voCwAc5MzyKmPMBGPMt4wxyYfdn2KM+ZGk6ZL2K7LMCghW3YGWV+2NF1f3tYu+uIm+uIm+APCRM6FDkdPZLpD0kTHmZWPMU8aYpZK2SnpAkSuRX+V53vYgiwQkSRuekCpbX2egTXs/lF5/suPrQQR9cRN9cRN9AeAjl5ZXrZT0S0WWUY2Q1E+RoFEmabGkBzzP6zqnq0fnNnq6NOqy9u1rOP2LNfTFTfTFTfQFgI+cCR2e522RxHU40DmkxD4eBQGiL26iL26iLwB85NLyKgAAAABdEKEDAAAAgFWEDgAAAABWEToAAAAAWEXoAAAAAGAVoQMAAACAVc6cMhcAEMX770uzZql/VZUWxhjSf+ZMKSND+t//WzrlFF/LS0TjVOLi+FQAAAkidACAyz77TFq8WBmSpsYa87e/RW5/+lOfimqfxqnExfGpAAASxPIqAAAAAFYROgAAAABYRegAAAAAYBWhAwAAAIBVhA4AAAAAVhE6AAAAAFhF6AAAAABgFaEDAAAAgFVcHBAAXNazp3TOOQqHw3p17dqoQ8484wylp6dHxjqscSpxjwUAdB2EDgBw2dCh0ooV+risTOcNGRJ1yJYFC5STk+NvXe3QOBUAQDfE8ioAAAAAVhE6AAAAAFhF6AAAAABgFaEDAAAAgFWEDgAAAABWEToAAAAAWEXoAAAAAGAVoQMAAACAVYQOAAAAAFYROgAAAABYRegAAAAAYBWhAwAAAIBVhA4AAAAAVhE6AAAAAFhF6AAAAABgVXLQBQAAjqCuTtq3T0l796pXjCFJe/dKe/ZImZlSsrv/rTdOJS6OTwUAkCD+SwcAl73+ulRYqMGS9sQaM3Jk5La4WCoo8KmwxDVOJS6OTwUAkCBCBwKV1KO3MkdOVPqg05WU2kMNNdUKb9uo3dV1ygm6uG6MvrhpV4/eenbkRK0ddLr2pfZQZk21zty2UdPe+of6V8eMJLCM1wsAtI3QgUCEa+t178pyZc98XCaU0mLbMUPGaNr8TZo6br9mT8pTekoooCq7H/ripnByqu48/7taPPx81R7Wl6IhY3T/+Ms0ZeMyza73lB5Qjd0RrxcAiB8HksN34dp6XflYsV58+7NWv6ib1DVIzxRv05WPFStcW+9zhd0TfXFTuN7TlVPv1DOjLmoVOJrUhlL0zKj/0JWr9tAXn/B6AYDEEDrguztfKNHaLbvjGrt2y27d+UKp5Yog0RdX3flGpdYOHi553pEHep7WflpLX3zC6wUAEkPogK92Voa1eP2OhPZZvH67dlUesFQRJPriqp2VYS3eGo78YMyRBzdupy/28XoBgMQROuCrheu2q7a+jXdsD1Nb72nha9stVQSJvrhq4brtqk2sLfTFB7xeACBxhA74Kt7lCId79YOKDq4Eh6IvbqIvbqIvAJA4Qgd8te9Ana/7IT70xU30xU30BQASR+iArzLT2neW5vbuh/jQFzfRFzfRFwBIHP8DwldnDOmjovc+TXi/M0/ua6EaNKEvbupqffrSoToAACAASURBVDn5ZGn+/PjHuqqr9QUA/EDogK+mFQzS7/75XkIHYaaEjKaNG2SxKtAXN3W1vvTtK11+edBVHL2u1hcA8APLq+CrAT3TNWVsdkL7TBk7SP17plmqCBJ9cRV9cRN9AYDEETrgu9mT8nXGkD5xjT1jSB/NnpRnuSJI9MVV9MVN9AUAEkPogO/SU0Kad02hJg07Vl59bdQxyUnSpYWDNe+aQqWnhHyusHuiL26iL26iLwCQGON5CV55qhMzxpTk5eXllZSUBF0KJJWVlSk3f7QyR1yo9MHDlZTaQw011Qpv26jlj/5SY/JOCbrEbom+uIm+uIm+AOhq8vPzVVpaWup5Xn5HPi4HkiNQDdV7tPfVRdr76qIW9/fp8ZuAKoJEX1xFX9xEXwCgbSyvAgAAAGAVoQMAAACAVYQOAAAAAFYROgAAAABYRegAAAAAYBWhAwAAAIBVhA4AAAAAVhE6AAAAAFhF6AAAAABgFaEDAAAAgFXJQRcAADiCjz6S5s9X1u7duiXGkKw//EHq00e64grphBN8LS8RjVOJi+NTAQAkiNABAC7bsUP66U/VR9I9scb86leR23PPdfov9capxMXxqQAAEsTyKgAAAABWEToAAAAAWEXoAAAAAGAVoQMAAACAVYQOAAAAAFYROgAAAABYRegAAAAAYBWhAwAAAIBVXBwQAFyWlCT16KEGz9P+/fujDjnmmGOUZExkrMMapxL3WABA10HoAACXjR0rVVVpW1mZhgwZEnXIltJS5eTk+FtXOzROBQDQDfFeEgAAAACrCB0AAAAArCJ0AAAAALCK0AEAAADAKkIHAAAAAKsIHQAAAACsInQAAAAAsIrQAQAAAMAqQgcAAAAAqwgdAAAAAKwidAAAAACwitABAAAAwCpCBwAAAACrkoMuoIkxpoekiZImSSqQlCMpJOl9SUskzfE8b19gBQLxqA1La/878v2ZP5CS04KtBxH0xU30xU30BYAFzoQOSZdJ+lPj9yWSXpaUJeksSXdKutQYc47neTsDqg9o2+tPSstmR75PzZAKrwu2HkTQFzfRFzfRFwAWuLS8qkbSw5JO9TzvdM/zpnmed5Gk0yS9LmmopPuDLBA4otqwVDTn4M9Fc6S6A8HVgwj64ib64ib6AsASZz7p8DzvCUlPRLn/I2PM9ZLWSLrEGJPqeV6N7wWiwx04EPsXWVlZmX+FdJCeby9Q38ryg3dUlqti2f2qHPqt4Ipqh82bN7drm6s6e19S3n5bA2bO1IDaWr0bY8yA8eNVm5KinQ89pNphw3ytLxFvv52imTMHSJKmnzZXvzizZV/+62tzNX/TNZKkhx7aqWHDaoMoMyG8Xjqf7OxsJSc78+cP0G0Yz/OCrqFNjcd7VDX+ONDzvI/a+TgleXl5eSUlJR1XHNpt0aJFmjZtWtBldIi0kPT+DZnKzmr54eGOvQ3KfWCfauoDKqyb6wp9GSdpXZxjCyS9ZrGWoxeZTVoorPdvGK3srPIWW3fsHajcB95QTX2aOsNsupqu8HqJx5YtW5STkxN0GYCz8vPzVVpaWup5Xn5HPq5Ly6uO5OTG21pJu4MsBIhmxpiUVr+oJSk7K0kzRqcEUBEk+uKqGWOebBU4JCk7q1wzRrf6wBs+4fUCwKbOEjpubLx92fO8NheXGmNKon1JyrVbJrqjtJA0a3zss7vcNiFNqSEfC4Ik+uKqtFBYs8bPibn9tglzlBriGAK/8XoBYJvzocMY81VJMxT5lOO/Ai4HaCXWu4NNeJcwGPTFTbE+5WjCpx3B4PUCwDanQ4cxZpik+ZKMpFs9z3sznv08z8uP9iWp8x3VB6eltvHuYJNZ43mX0E/0xU2poYYjfsrRZNb43yo11OBDRZB4vQDwh7OhwxiTrci1Oo5V5MKAvwu4JKCVGaOP/O5gk0G9knQN7xL6hr64acbovUf8lKPJoF4f6prRe32oCBKvFwD+cPLsVcaYfpKKFLk2x+OSZngdUChnr3LLO++8o2ExTu+5fPly988uUl+j7CUXK7k6vutV1vU4Tju++aIUSrVc2NHZvHmzLrjggqjbli1bptxcxw+N6mJ9SX3zTQ2cPDmuseXPP6+akSMtV9RO9TU6bsHXdEztJ3EN359yvD759gvO9qUJr5fOh1PmAkdm6+xVzr3qjDE9Jf1NkcDxZ0nXdUTggHvS09NjbsvJyXE/dBT/SYrzF7UkJVd/opzd/5IKrrVYlF25ubn0xW+7dsU9dODAgZKr/Sn+kxRn4JCkY2o/drsvceD1AgAHObW8yhiTJukvipzMfamkSz3P6yJnBkeXUneg5VV748XVfe2iL26iL26iLwB85MwnHcaYkKRnJJ2nyNKqS7jyOJy14Qmpsu216a3s/VB6/UneJbSlK/bluOOkG2/U3r179djjj0cdcs3VVysrKysy1kVdsS9dAX0B4CNnQoekH0r6RuP3n0p6yBgTbdwtnud96ltVQDSjp0ujLmvfvobTv1jTFfsyeLB0//3aXVamm2KEjsm3364sl5fxdMW+dAX0BYCPXAodxx7y/TdijpLuUCSUAMFJiX08CgJEX9xEX9xEXwD4yJljOjzPu8PzPBPHV1nQtQIAAACInzOhAwAAAEDXROgAAAAAYBWhAwAAAIBVhA4AAAAAVhE6AAAAAFhF6AAAAABgFaEDAAAAgFWEDgAAAABWEToAAAAAWEXoAAAAAGAVoQMAAACAVYQOAAAAAFYROgAAAABYlRx0AQCAI9i7V9qwQekff6xzYgxJf/VVqaxMGjNGysrys7qENE4lLo5PBQCQIEIHALhs0ybpvPN0vKQVscZcemnktrhYKijwp652aJxKXByfCgAgQSyvAgAAAGAVoQMAAACAVYQOAAAAAFYROgAAAABYRegAAAAAYBWhAwAAAIBVhA4AAAAAVhE6AAAAAFhF6AAAAABgFVckBwCXnX669O9/68MPP9TEr3wl6pC/L12qE088UTr5ZJ+LS0zjVOLi+FQAAAkidACAy445RsrPV21GhkpjDKk99VQpJ8fPqtqlcSoAgG6I5VUAAAAArCJ0AAAAALCK0AEAAADAKkIHAAAAAKsIHQAAAACsInQAAAAAsIrQAQAAAMAqQgcAAAAAq7g4IAKV1KO3MkdOVPqg05WU2kMNNdUKb9uo3dV1ygm6uG6MvriJvriJvgBA24zneUHX4BtjTEleXl5eSUlJ0KV0e+Haet3y9Fq98O+dMqGUVtuTk6Sp4wZr9qQ8paeEAqiwe6IvbqIvbqIvALqi/Px8lZaWlnqel9+Rj8vyKvguXFuvKx8r1otvfxb1F7Uk1TVIzxRv05WPFStcW+9zhd0TfXETfXETfQGAxBA64Ls7XyjR2i274xq7dstu3flCqeWKINEXV9EXN9EXAEgMoQO+2lkZ1uL1OxLaZ/H67dpVecBSRZDoi6voi5voCwAkjtABXy1ct1219YkdR1Rb72nha9stVQSJvriKvriJvgBA4ggd8FW8yxEO9+oHFR1cCQ5FX9xEX9xEXwAgcZwyF77ad6DO1/0QH/ripq7Wl82bpdtvj2/sL34h5ebarae9ulpfAMAPhA74KjOtff/k2rsf4kNf3NTV+rJ7t/T00/GN/fGP3Q0dXa0vAOAHllfBV2cM6dOu/c48uW8HV4JD0Rc30Rc30RcASByhA76aVjBIKSGT0D4pIaNp4wZZqggSfXHVtIJBSkmsLfTFB7xeACBxhA74akDPdE0Zm53QPlPGDlL/nmmWKoJEX1w1oGe6ppyUHvnBa+NsSY3b6Yt9vF4AIHGEDvhu9qT8uJcnnDGkj2ZPyrNcEST64qrZo3rqjG0bJdPGO+vG6Ix+KfTFJ7xeACAxhA74Lj0lpHnXFGrSsGPl1ddGHZOcJF1aOFjzrilUekrI5wq7J/ripvSQ0bxFs3XpG39TSoy+pNTX6tI3/qZ543vTF5/wegGAxBivrY/suxBjTEleXl5eSUlJ0KVAUllZmXLzRytzxIVKHzxcSak91FBTrfC2jVr+6C81Ju+UoEvsluiLY9atkwoLJUm7evTWwhEX6tXBw7UvtYcya6p15raNmvbWP9S/eo9UXCwVFARccGyHTKVNjk+lGa8XAF1Nfn6+SktLSz3Py+/Ix+X8fQhUQ/Ue7X11kfa+uqjF/X16/CagiiDRF1f1r96j619dpOsP6wuCxesFANrG8ioAAAAAVvFJBwC4LDNTOusshQ8c0Pr166MOGTt2rNLT0iJjHdY4lbjHAgC6DkIHALhs2DBp9Wp9XFam8UOGRB2yZfFi5eTk+FtXOzROBQDQDbG8CgAAAIBVhA4AAAAAVhE6AAAAAFhF6AAAAABgFaEDAAAAgFWEDgAAAABWEToAAAAAWEXoAAAAAGAVoQMAAACAVYQOAAAAAFYROgAAAABYRegAAAAAYBWhAwAAAIBVhA4AAAAAVhE6AAAAAFiVHHQBAIAjqK+XwmGZ6mr1iDHEVFdLVVVSeroUCvlaXiIapxIXx6cCAEgQoQMAXLZhg1RYqJMkVcUak58fuS0ulgoKfCoscY1TiYvjUwEAJIjlVQAAAACsInQAAAAAsIrQAQAAAMAqQgcAAAAAqwgdAAAAAKwidAAAAACwitABAAAAwCpCBwAAAACrCB0AAAAArOKK5ADgsiFDpMcf16effqpbbr016pB777lH/fr1i4x1WONU4h4LAOg6CB0A4LJ+/aSrrtK+sjLNixE67pgyRf1ycvytqx0apwIA6IZYXgUAAADAKkIHAAAAAKsIHQAAAACsInQAAAAAsMqpA8mNMWMlXSipUNIZkgZKOuB5XnqghQHxqg1La/878v2ZP5CS04KtBxH0xU30xU30BYAFToUOSf8l6etBFwG02+tPSstmR75PzZAKrwu2HkTQFzfRFzfRFwAWuLa86hVJv5A0SdLxAdcCJKY2LBXNOfhz0Ryp7kBw9SCCvriJvriJvgCwxKlPOjzP+/WhPxtjgioFPjhwIPYvsrKyMv8K6SA9316gvpXlB++oLFfFsvtVOfRbwRXVDps3b27XNlfRFzfRFzd1lb4cSXZ2tpKTnfrzB+gWjOd5QdcQkzHGUwce02GMKcnLy8srKSnpiIfDUVq0aJGmTZsWdBkdIi0kvX9DprKzWn54uGNvg3If2Kea+oAK6+boi5voi5u6S1+2bNminE5wMU0gKPn5+SotLS31PC+/Ix/XteVVQKc0Y0xKq1/UkpSdlaQZo1MCqAgSfXEVfXETfQFgU5cMHcaYkmhfknKDrg1dT1pImjU+9tldbpuQptSQjwVBEn1xFX1xE30BYFuXDB2An2K9O9iEdwmDQV/cRF/cRF8A2NYlQ4fnefnRviR1vqP64LTUNt4dbDJrPO8S+om+uIm+uIm+APBDlwwdgF9mjD7yu4NNBvVK0jW8S+gb+uIm+uIm+gLAD5y9CoF55513NGzYsKjbli9f7v7ZReprlL3kYiVX74xreF2P47Tjmy9KoVTLhR2dzZs364ILLoi6bdmyZcrNdfzQKPriJvripi7alyPhlLnAkdk6exWvOgQmPT12lszJyXE/dBT/SYrzF7UkJVd/opzd/5IKrrVYlF25ubn0xW8ffywtWKCs3bt1Y4who1euVJ+NG6Vvf1s63tHrqna1vsSB1wsAHEToANqj7kDLq/bGq2iONHq6lNz2+mm0Q1fsy/bt0k03qY+k+2ONueuuyO3ZZ7sZOrpiX7oC+gLARxzTAbTHhiekQ6/aG6+9H0qvP9nx9SCCvriJvriJvgDwkVOfdBhjLpb0X4fdnWqMefWQn+/yPO8lH8sCWhs9XRp1Wfv2NZz+xRr64qbGvqxfL33pnPh2+ddKaexY0RebeL0A8JFToUNSf0lnHHafOey+/v6VA8SQ0iHnNkBHoy9uauxLQ7JUXRvfLg3JkjrvscqdA68XAD5yKnR4njdX0tyAywAAAADQgTimAwAAAIBVhA4AAAAAVhE6AAAAAFhF6AAAAABglVMHkgMADmOMlJwsT1JdXV3UIcnJyTJNYx3WOJW4xwIAug5CBwC4bNw4qbZWW8vKNGTIkKhDtrz3nnJycvytqx0apwIA6IZYXgUAAADAKkIHAAAAAKsIHQAAAACsInQAAAAAsIrQAQAAAMAqQgcAAAAAqwgdAAAAAKwidAAAAACwiosDIibP81Sn6FdA7gi1XuyrhNV6tUfcjohkJctw6WYAAOA4Qgdi2la5TecsOEdGRjOnzVRWZlaHPv6n9Z/G3LawfqH61ffr0Ofriq4JXaMUpQRdBgAAwBGxvAoxPff2c9r60VaVfVSmNza9EXQ5AAAA6KQIHYhpccni5u/Xv7M+wEoAAADQmRE6EFV5ZbnWbF/T/PO7W9/V3qq9AVYEAACAzorQgaiWlC6RJ6/5Z8/z9MY7LLECAABA4ggdiGph6cJW97HECgAAAO1B6EAr5ZXlWr1tdav7WWIFAACA9uCUuWjl8KVVTZqWWH1p7JcCqAropjZulL79bQ2srdW/YwwZOHGilJIiLVggDR/ua3mJaJxKXByfCgAgQYQOtBJtaVWT9e+sJ3QAfgqHpdJSpUrKjzXmvfcOjnVY41TiHgsA6DpYXoUWPtz7YdSlVU1YYgUAAIBEETrQwpK3oy+tasJZrAAAAJAoQgdaWFS6qM0xnMUKAAAAiSB0oFlbS6uasMQKAAAAiSB0oFlbS6uasMQKAAAAiSB0oFk8S6uasMQKAAAA8eKUuZAkNXgNumbUNbpq5FWSpPkb52tF2YoWY4bmDFVBfoEkyRgjz/NkjPG5UgAAAHQ2hA5IkpJMkq4efXXzzxs+2tAqdJzQ/wSdPepsnysDAABAZ0foAACXDRggzZypvXv36sn586MOmX7FFcrKyoqMdVjjVOIeCwDoOggdAOCyk06Sfv977S4r0w9jhI6L77pLWTk5/tbVDo1TAQB0Q4QOBCqpR29ljpyo9EGnKym1hxpqqhXetlHV1fzTDFKsvuyurlNO0MV1Y/TFTfQFANrGX3YIRG2tVFSUreyZj8uEUlpsO2bIGD2zoEEjxlTrwov3KDklxoOgw4Vr63XvyvKYfZk2f5Omjtuv2ZPylJ4SCqjK7oe+uIm+AED8OGUufFdbKy2c11/vbOrX6hd1k4aGJL3xWqaenddftbU+F9hNhWvrdeVjxXrx7c9i9qWuQXqmeJuufKxY4dp6nyvsnuiLm+gLACSG0AHfLXupt7aVpUltXojQ07ayNC17qbcfZXV7d75QorVbdsc1du2W3brzhVLLFUGiL66iLwCQGEIHfLWvMklvbcho/Kmta3xEtr+1IUNV+/inatPOyrAWr9+R0D6L12/XrsoDliqCRF9cRV8AIHH8JQdfvbk+Qw0NiV1QsKHB6M3XMtoeiHZbuG67auvb+uSppdp6Twtf226pIkj0xVX0BQASR+iAr7ZtSWvXflvbuR/iE+8ykcO9+kFFB1eCQ9EXN9EXAEgcoQO+qqlJ7FOOo90P8dl3oM7X/RAf+uIm+gIAiSN0wFepqYktSTja/RCfzLT2nT27vfshPvTFTfQFABJH6ICvBg9p34GUJ7VzP8TnjCF92rXfmSf37eBKcCj64ib6AgCJI3TAVyPHVikpKbFPLZKSPI0cV2WpIkjStIJBSgkltoQtJWQ0bdwgSxVBoi+uoi8AkDhCB3yV2bNBI8Y0BYi2r9MhSSPGVCkjs8FqXd3dgJ7pmjI2O6F9powdpP49OcDfJvriJvoCAIkjdMB3F1y8R4NzDiie63QMzjmgCy/e40dZ3d7sSflxLxs5Y0gfzZ6UZ7kiSPTFVfQFABJD6IDvUlKkaVfu0tChn8qrr406JimpQaPG7dO3rtyl5BSfC+ym0lNCmndNoSYNOzZmX5KTpEsLB2veNYVKTwn5XGH3RF/cRF8AIDGcSgOBSEmRJozfoX9c/2NljrhQ6YOHKym1hxpqqhXetlE//MNYDc7rHXSZ3U56Skg/OWegHpp5cdS+LH/0lxqTd0rQZXY79MVN9AUA4kfoQKAaqvdo76uLtPfVRS3u79FjZEAVQYrdlz49fhNQRd1YZaX01ltK++gjnVm9R3p1UeTrECeUTpM++0QaMULq2TOgQtvWOJW4OD6VFni9AEDbCB0A4LJ33pHGj9cJklbHGjN1auS2uFgqKPCpsMQ1TiUujk8FAJAgjukAAAAAYBWhAwAAAIBVhA4AAAAAVhE6AAAAAFhF6AAAAABgFaEDAAAAgFWEDgAAAABWEToAAAAAWEXoAAAAAGAVVyQHAJfl50uvv67y8nJ99eKLow7560svaeDAgdKpp/pcXGIapxIXx6cCAEgQoQMAXNajhzRqlGp699abMYbU5OVJOTl+VtUujVMBAHRDLK8CAAAAYBWhAwAAAIBVhA4AAAAAVhE6AAAAAFhF6AAAAABgFaEDAAAAgFWEDgAAAABWEToAAAAAWEXoAAAAAGAVoQMAAACAVYQOAAAAAFYROgAAAABYRegAAAAAYBWhAwAAAIBVhA4AAAAAViUHXQAA4Ag++ED6xS/Ub98+PR5jSL9bbpEyM6Xbb5dOPtnX8hLROJW4OD4VAECCCB0A4LKKCmnePGVKuirWmCVLIrfXX+/0X+qNU4mL41MBACSI5VUAAAAArCJ0AAAAALCK0AEAAADAKkIHAAAAAKsIHQAAAACsInQAAAAAsMq50GGMSTfG3GmMedcYEzbGlBtjHjPGZAddG9CWUF2tRq5ZoZFrViipri7octCkNiytuj/yVXcg6GrQKC0U1q1n3a9bz7pfqSH64gxeLwAscOo6HcaYdEn/lHSWpI8k/UVSjqSrJX3NGPNFz/M2B1chcGRD31inM//vXyVJtampKh13VsAVQZL0+pPSstmR71MzpMLrgq0HkqQZY57Uby6M9GVfTYYefo2+OIHXCwALnAodkm5TJHC8Immi53n7JMkYc7Ok+yQ9Jumc4MoDYgvV1WrU6uXNP49evVzvjCpUQ7JrL7NupjYsFc05+HPRHGnMd6TktOBqSkSPHtK4cTpQU6O33nor6pARI0YoLTU1MtZhjVORJKUmhXX7uQf7MvvLc/RG0ndU25DWPBYB6OyvFwDOcuavIWNMiqQfNf54fVPgkCTP8+YYY66U9CVjzFjP89YHUiQ6VEN9Q8xtn+34zMdKOsbYTRuUWfl588+ZlZ9r0D+Xa/1po609Z1moTCkmpWMfs6ysXdtc1fPtBepbWX7wjspyVSy7X5VDvxVcUYnIyJAWLVJZWZnOO++8qEOW/+53ysnJifzgcI8apyKpsS9rD/bluGPK9cKdLfvi8FSa8XrpfLKzs5XMm0GA71x61Y2X1FvSZs/zXo+yfbGkEZImSSJ0dAGff/R5zG33nXefj5UcvbSQ9P4NmVJWy8OkTv/nUn3960tUU2/neX+un9t54Bhi/dHrqlh92f/3u3Xa135mrS9+oy9uoi9u2rJly8GQDsA3Lh1IPrLxdkOM7RsOGwc4Y8aYFGVntX45ZWclacbojv0kAvGjL26iL26iLwBscil0DG683RFj+47DxsVkjCmJ9iUptyMKBQ6VFpJmjY+93vm2CWlKDflYECTRF1fRFzfRFwC2uRQ6Mhtvq2NsrzpsHOCEWO8ONuFdwmDQFzfRFzfRFwC2uRQ6TOOt18b2Nnmelx/tSxKn20WHSm3j3cEms8bzLqGf6Iub6Iub6AsAP7gUOiobbzNibG86geK+GNsB380YfeR3B5sM6pWka3iX0Df0xU30xU30BYAfXDp71bbG21hXHs8+bBw6uZPGnaSfLP9J1G1Zx2cplOz2W2qh+jrN/Msfper4cvCvpw7QSfOvU32o415200LTOvyUueFwWBs3boy6bfjw4UpPT+/Q5+tw9TXKXnKxVL0zruH/Z2qO/teCF6VQquXCjg59cRN96Xyys2P9mQHAJpdCx5uNt2NibG+6P/rVsdDppKan6gvjvxB0Ge2W99oaZcUZOCQpq7pSZ+/ZotKxX+ywGnJCOR0eOiRp6NChHf6Yvin+U9x/QElScvUnytn9L6ngWotFdQz64ib6AgBtc2l51WpJn0vKNcZEu5ralMbbF/0rCYguqa5Oow+5+ni8Rq/6v0qqq7NQESRJdQdaXk05XkVzIvvCDvriJvoCwEfOhA7P82okPdj444PGmOZjO4wxNytyYcBVnuetC6I+4FBD3yhucfXxeGVWfq6hb/JP2JoNT0iHXk05Xns/lF5/suPrQQR9cRN9AeAjl5ZXSdLdki6QdJak94wxRZJOknSGpApJVwdYG9Bs06gCvTtiXLv29ZLiPhEbEjV6ujTqsvbta9w+hqhToy9uoi8AfORU6PA8L2yMOU/SLEmXSZos6TNJ8yT9l+d524OsD2hSn8wZXJyU4vhBu90VfXETfQHgI6dChyR5nrdf0u2NX8D/a+/+gzWr6zqAvz/AooCApKQQILaVwJq2gCSKv9KppmxYcbFyknKUGk3HGfwxmL9DLZvUEkccFPyRM426iTqZP6ppkdJaWQ1/oC5tCv5ADFHRDRZyv/1xnht4uZd2ufd7z7P3eb1m7pz7POc8z/k8c+b58T7f8/1+AQDYy01Nnw4AAGB1EjoAAICupu7yKgBuZ9eu5Ec/2r1t99032Wd6zyWtopcCwB7ykQ4wzbZuTfbff/f+tm4du9o7tYpeCgB7SOgAAAC6EjoAAICuhA4AAKAroQMAAOhK6AAAALoSOmBEN3z/huzatWvsMgAAujJPB4zoA5s/kC9+5Ys58fgTc9LxJ2XtUWuzj8kJAIBVRuiAkd2448ZsvnxzNl++OYccdIgAAgCsOkIHTBEBBABYjYQOdtuWz2/Jtqu3jV3GqvLdG7+76DoBhCTJsccmF164+9tOsVX0UgDYQ0IHu23HTTuy46YdY5cxk+YHkPXHrc8jT3xkcsTYldHd4YcnZ589dhXLYhW9FAD2kNABe5kbd9yYH5jOSQAADENJREFUS7demv3X7C90AAB7BddosKCjDz06D77Pg8cug0VUKofe49CxywAA2C1aOljQuaedm2ee/Mwc+ho/bKdFpfLwYx6ejSdszIbjN+TIg4/Mft7CAMBewC8WFnXAmgPyno3vGbuMVe38LefnsmsuW3R9pXLaMaflSeuelDOOPyNHHnzkClYHALA8hA4WtWbfNTlz3Zljl7GqfXDbB+8QOgQNAGC1ETpgCswFjTNPODNPPOGJggYAsKoIHTCi9fddn1OOPEXQAABWNaEDRnTOqeeMXQIAQHeGzAUAALoSOgAAgK6EDgAAoCuhAwAA6EroAAAAuhI6AACAroQOAACgK/N0AEyz665LNm3avW03bkzuc5++9SzBKnopAOwhoQNgml1zTfKsZ+3etqecMtW/1FfRSwFgD7m8CgAA6EroAAAAuhI6AACAroQOAACgK6EDAADoSugAAAC6EjoAAICuhA4AAKAroQMAAOjKjOQA0+whD0laG7uKZbGKXgoAe0hLBwAA0JXQAQAAdCV0AAAAXQkdAABAV0IHAADQldABAAB0JXQAAABdCR0AAEBXQgcAANCV0AEAAHQldAAAAF0JHQAAQFdCBwAA0JXQAQAAdCV0AAAAXQkdAABAV0IHAADQldABAAB0JXQAAABdCR0AAEBXQgcAANCV0AEAAHQldAAAAF0JHQAAQFdCBwAA0JXQAQAAdCV0AAAAXQkdAABAV0IHAADQldABAAB0JXQAAABdCR0AAEBXQgcAANCV0AEAAHQldAAAAF0JHQAAQFdCBwAA0JXQAQAAdCV0AAAAXQkdAABAV0IHAADQldABAAB0JXQAAABdTUXoqKqDquopVXV+VW2pqp1V1arq3LFrAwAAlma/sQuY+Nkk7xy7CAAAYPlNRUtHkh8kuSjJHyQ5Mcmrxi0HAABYLlPR0tFa257k6XO3q+r0EcsBAACW0bS0dAAAAKuU0AEAAHQ1FZdXLbeq+sIiq47bvn171q1bt6L1AADA3mD79u1JcvRyP++qDB13Yp+dO3fuuvLKK780diHcwdrJcvuoVTCf4zKdHJfp5LhMJ8dlejk20+m4JAcs95MuS+ioqk1JHriHDzurtbZlOfY/X2ttwaaMuRaQxdYzHsdmOjku08lxmU6Oy3RyXKaXYzOd7uSKoSVZrpaOY5M8YA8fc+Ay7RsAAJhiyxI6WmsnL8fzAAAAq4/RqwAAgK6EDgAAoCuhAwAA6Kpaa2PXkCSpqkuSHDG5eVSSn0rytSTfnNx3bWvtCWPUBgAA3HXTFDq+muR+d7LJ1a21Y1emGgAAYLlMTegAAABWJ306AACAroQOAACgK6EDAADoSugAAAC6EjoAAICuhA4AAKCrmQ4dVXVQVT2lqs6vqi1VtbOqWlWdO3Zts6Cq7l5Vr6iqbVV1c1V9s6ourqqjxq5tVlXVSVV1blW9r6q+MXk/3Dx2XbOuqg6sqg1VdVFVfbaqbqyqHVV1RVW9tKruMXaNs6qqzpm8X66qqu9Pvkeurqp3VNW6setjUFU/UVXfnnymfWnsemZZVW2eHIfF/n517BpnWVXdt6peP/ltdlNV3VBVW6vqz5b83LM8T0dV/UKSzyyw6oWttT9d6XpmSVXdPck/JnlYkmuTXJbk2CSnJPmvJKe21raPVuCMqqr3Jzl93t07W2t3H6MeBlX19CRvmdz8QpIrkxyS4f1zcJIvJXlUa+3b41Q4u6rq+iQHJflskm9M7l6X5OeS3JJkQ2vtwyOVx0RVvT3JWUkqyZdba8eNW9HsqqrNSR6V5G+S/HCBTV7bWvvcihZFkqSqTk3yd0numeF75vMZvmNOSHJUa22/pTz/kh68CvwgyUVJtiT5VJInJnnRqBXNjj/K8IPpk0l+ubX2w2Q4a5jktUkuzvChxMr6ZJIrMrwfPpXkW+OWw8QtSS5I8vrW2lVzd1bVEUk+lGR9kr9I8uRxyptppyfZ2lr7sRbBqnpGkjcleWtVHdNa+9Eo1ZGqemyS301yYZLfH7kcbvO81tpXxy6CQVUdmSFw3C3JGa21S+atP2XJ+5jllo75qurlSV4WLR1dVdWaJN/OkKRPbK19Zt76K5I8KMnJrbWtI5TIRFW1aOmYapMzU59IsjPJIa21W0YuiYmquirJzyRZ11q7cux6ZlFVHZChFeqWJBuSbIuWjlHdrqXj/kLH9KiqdyZ5SpJnt9be2GMfM92ng9GcliFwbJ8fOCY2TZa/sXIlwV7risnybknuNWYh3MFc64YgOJ6XJVmb5BlJbh25FphKVXVYkicl+X6St/baz6xfXsU4HjxZfnqR9Z+etx2wuJ+eLG9NcsOYhXCbqjoryQMynFn/z5HLmUlV9aAkz03yttbax6vq2HErYp6nVdW9kuzK8D55f2vtmpFrmlUPz3Di6h+S3FpVGzOcIF6Toc/ge1pr1y11J0IHYzhmsvz6Iuu/Pm87YHHPmSw/0lrbOWolM6yqnp+hA/lBSY6f/P/NJE9ure0as7ZZVFX7ZBh84XtJXjByOSzsxfNu/3lVnddaO2+Uambb3Eh712UY2OfUeev/pKqe2lp771J24vIqxjA3vOd/L7J+x7ztgAVU1a8leVqGVo6XjFzOrPuVDJ2VN2b4Av9ahsChX9o4np1hNMTnt9a+M3Yx/JiPZ+g7sDbJgRlaBF+U5H+S/HFVPedOHksfh02WZ2XoU/u0JIcnuX+S12U4mfKuSevhXbZXt3RU1aYkD9zDh53VWtvSox52W02Wi41iUIvcD0xU1fFJ3pXh/fL81toV/89D6Ki19rgkqap7Jvn5JC9NsrmqXtxae9Woxc2Yqjo6ySuTXNpae/vI5TBPa+2l8+7aluTVVXV5ko8meUVVXdhau2nlq5tZ+06W+yX5w9baxZPb1yd5blUdk+GEyguS/M5d3cleHToyzOvwgD18zIEd6mDP/GCyPGiR9XPHaKHxu2HmTSbQ/EiGs1Ova6395cglMdFa+16SyyatUJ9Mcl5Vfay19qmRS5slb0qyf4bO4+wlWmsfmwSPk5M8NMk/jVzSLJn7XbYryTsWWH9xhtDx6KXsZK8OHa21k8eugbtkrqPYYjOPHzVvO2Ciqu6d5O8z9Hl6W5LnjVsRC2mt3VpV705yUoaR+ISOlfP4DH05Lqj6sYbzuaG/j5kM25okj5+bJ4qpcFWG0HHE2IXMmK9Olt9apG/g3PqfXMpO9urQwV5r7jKQExdZP3f/Z1egFthrVNXBST6c5Lgk70tydjPZ0jS7frI8fNQqZtM9s/gEswfcbp3fQdNlrm+BILiy5qYvOKyqaoHvlbnh2Jd0XHQkZwz/kmEs6LVVtX6B9Rsny79duZJgulXV3ZJ8IMNZwI8m+W2zXE+9uR+220etYsa01mqhvwydYpNhcsC5+783Zq3cpqoOT/KIyc3FhtSng9ba55J8JUMg/8UFNnn0ZLmk4yJ0sOImMybPzXb5xqr6v74dVXVOhpET/tk10DCoqn2T/HWSx2QYzvAMM4+Pr6oeUVW/WVX7zbt/TVU9O8MIPTclefcoBcKUqaqHVtVjat51b5M5VC7J0Nfzg621xYbUp5/XTJZvmFzGmySpqpMyzHeTJG9eyg5mvlmxqi7JbdcOzvUleGZVbZj8f21r7QkrX9mq98okj0vysCRXVdVlSe6XIWF/J8lTR6xtZlXVr+eOQ6/uX1X/ervb57XWPrSCZZE8K8nc59D1Sd407zt7zvNaa9cvtIIu1mboV3N9VW3N8Nl17wyjVx2R5OYkv9da+9p4JcJUOS7De+baqtqW5FsZfnudlKHPzReSnD1eeTPtLUkem+TMJF+uqk9kmLrgYRkGZnhLa23TUnYw86EjyfoMP3Zv7+jJX5JcvbLlzIbW2s1V9ZgkL0zy5CQbknw3w6gJL/ElPZrDc8em1Zp3n+vTV95ht/v/zk6CvDy39SOgv0uTvDrDZVQPyhA4bsnQ6XJTkje01v5jtOpg+vxbkgsyfKeckGEm7B1J/j3Je5NcYKjccbTWdlXVbyXZnOTpSX4pw9QGlyd5c2vtr5a6j9IHEQAA6EmfDgAAoCuhAwAA6EroAAAAuhI6AACAroQOAACgK6EDAADoSugAAAC6EjoAAICuhA4AAKAroQMAAOhK6AAAALoSOgAAgK6EDgAAoCuhAwAA6EroAAAAuhI6AACAroQOAACgq/8FuTrw+EB9gY4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1050x750 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(7, 5))\n",
    "ax = plt.gca()\n",
    "lat = DualSquare(4, 4, None, bc='periodic')\n",
    "sq = Square(4, 4, None, bc='periodic')\n",
    "sq.plot_coupling(ax, linewidth=3.)\n",
    "ax.plot([2., 2.], [-0.5, 4.3],  'r:', linewidth=5., label=\"Wilson\")\n",
    "ax.plot([2.5, 2.5], [-0.5, 4.3], 'b:', linewidth=5., label=\"t'Hooft\")\n",
    "lat.plot_sites(ax)\n",
    "lat.plot_basis(ax, origin=-0.25*(lat.basis[0] + lat.basis[1]))\n",
    "ax.set_aspect('equal')\n",
    "ax.set_xlim(-1, 6)\n",
    "ax.set_ylim(-1)\n",
    "ax.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Extending the existing model\n",
    "\n",
    "\n",
    "We will implement additional parameters which allow to optionally add terms for the Wilson and t'Hooft loop operators $W, H$ defined above, $$ H \\rightarrow H - \\mathtt{J_{WL}} W -  \\mathtt{J_{HL}} H $$\n",
    "\n",
    "Note that here we only want to add a single loop for each of them, not one at each possible starting point. Hence, we use the method `add_local_term`, not `add_multi_coupling_term` or `add_coupling_term`. (In the case of infinite MPS, it will still be one loop per MPS unit cell, but that's okay as well.)\n",
    "\n",
    "The sign of $\\mathtt{J_{WL}},\\mathtt{J_{HL}}$ will allow us to select sectors with $\\langle \\psi | W |\\psi \\rangle = \\pm 1$ and $\\langle \\psi | H |\\psi \\rangle = \\pm 1$, respectively.\n",
    "Note that so far the constraints commute with the Hamiltonian and with each other.\n",
    "However, this is no longer the case if we add another global field  $H \\rightarrow H -\\mathtt{h} \\sum_i \\sigma^z_i$, which makes the Hamiltonian no longer exactly solvable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "class ExtendedToricCode(ToricCode):\n",
    "\n",
    "    def init_terms(self, model_params):\n",
    "        ToricCode.init_terms(self, model_params)  # add terms of the original ToricCode model\n",
    "        \n",
    "        Ly = self.lat.shape[1]\n",
    "        J_WL = model_params.get('J_WL', 0.)\n",
    "        J_HL = model_params.get('J_HL', 0.)\n",
    "        # unit-cell indices:\n",
    "        # u=0: vertical links\n",
    "        # u=1: horizontal links\n",
    "       \n",
    "        # Wilson Loop\n",
    "        x, u = 0, 0 # vertical links\n",
    "        self.add_local_term(-J_WL, [('Sigmaz', [x, y, u]) for y in range(Ly)])\n",
    "\n",
    "        # t'Hooft Loop\n",
    "        x, u = 0, 1 # horizontal links\n",
    "        self.add_local_term(-J_HL, [('Sigmax', [x, y, u]) for y in range(Ly)])\n",
    "            \n",
    "        h = model_params.get('h', 0.)\n",
    "        for u in range(2):\n",
    "            self.add_onsite(-h, u, 'Sigmaz')\n",
    "        \n",
    "    def wilson_loop_y(self, psi):\n",
    "        \"\"\"Measure wilson loop around the cylinder.\"\"\"\n",
    "        Ly = self.lat.shape[1]\n",
    "        x, u = 0, 0 # vertical links\n",
    "        W = TermList.from_lattice_locations(self.lat, [[(\"Sigmaz\",[x, y, u]) for y in range(Ly)]])\n",
    "        return psi.expectation_value_terms_sum(W)[0]\n",
    "    \n",
    "    def hooft_loop_y(self, psi):\n",
    "        \"\"\"Measure t'Hooft loop around the cylinder.\"\"\"\n",
    "        Ly = self.lat.shape[1]\n",
    "        x, u = 0, 1 # horizontal links\n",
    "        H = TermList.from_lattice_locations(self.lat, [[(\"Sigmax\",[x, y, u]) for y in range(Ly)]])\n",
    "        return psi.expectation_value_terms_sum(H)[0]\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## iDMRG\n",
    "We now run the iDMRG algorithm for the four different sectors to obtain the four degenerate ground states of the system. \n",
    "\n",
    "To reliably get the 4 different ground states, we will add the energetic constraints for the loops.\n",
    "In this case, the original Hamiltonian will completely commute with the loop operators, such that you directly get the correct ground state. However, this is a very peculiar case - in general, e.g. when you add the external field with $h \\neq 0$, this will no longer be the case.\n",
    "We hence demonstrate a more \"robust\" way of getting the four different ground states (at least if you have an idea which operator distinguishes them): \n",
    "For each sector, we run DMRG twice: first with the energy penalties for the $W,H$ loops to get an initial $|\\psi\\rangle$ in the desired sector. Then we run DMRG again for the clean toric code model, without the $W,H$ added, to make sure we have a ground state of the clean model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "dmrg_params = {\n",
    "    'mixer': True,\n",
    "    'trunc_params': {'chi_max': 64, \n",
    "                     'svd_min': 1.e-8},\n",
    "    'max_E_err': 1.e-8,\n",
    "    'max_S_err': 1.e-7,\n",
    "    'N_sweeps_check': 4,\n",
    "    'max_sweeps':24,\n",
    "}\n",
    "model_params = {\n",
    "    'Lx': 1, 'Ly': 4, # Ly is set below\n",
    "    'bc_MPS': \"infinite\",\n",
    "    'conserve': None,\n",
    "}\n",
    "\n",
    "\n",
    "def run_DMRG(Ly, J_WL, J_HL, h=0.):\n",
    "    print(\"=\"*80)\n",
    "    print(f\"Start iDMRG for Ly={Ly:d}, J_WL={J_WL:.2f}, J_HL={J_HL:.2f}, h={h:.2f}\")\n",
    "    model_params_clean = model_params.copy()\n",
    "    model_params_clean['Ly'] = Ly\n",
    "    model_params_clean['h'] = h\n",
    "    model_clean = ExtendedToricCode(model_params_clean)\n",
    "    model_params_seed = model_params_clean.copy()\n",
    "    model_params_seed['J_WL'] = J_WL\n",
    "    model_params_seed['J_HL'] = J_HL\n",
    "    model = ExtendedToricCode(model_params_seed)\n",
    "    psi = MPS.from_lat_product_state(model.lat, [[[\"up\"]]])\n",
    "\n",
    "    eng = TwoSiteDMRGEngine(psi, model, dmrg_params)\n",
    "    E0, psi = eng.run()\n",
    "    WL = model.wilson_loop_y(psi)\n",
    "    HL = model.hooft_loop_y(psi)\n",
    "    print(f\"after first DMRG run: <psi|W|psi> = {WL: .3f}\")\n",
    "    print(f\"after first DMRG run: <psi|H|psi> = {HL: .3f}\")\n",
    "    print(f\"after first DMRG run: E (including W/H loops) = {E0:.10f}\")\n",
    "\n",
    "    E0_clean = model_clean.H_MPO.expectation_value(psi)\n",
    "    print(f\"after first DMRG run: E (excluding W/H loops) = {E0_clean:.10f}\")\n",
    "\n",
    "    # switch to model without Wilson/Hooft loops\n",
    "    eng.init_env(model=model_clean)\n",
    "\n",
    "    E1, psi = eng.run()\n",
    "\n",
    "    WL = model_clean.wilson_loop_y(psi)\n",
    "    HL = model_clean.hooft_loop_y(psi)\n",
    "    print(f\"after second DMRG run: <psi|W|psi> = {WL: .3f}\")\n",
    "    print(f\"after second DMRG run: <psi|H|psi> = {HL: .3f}\")\n",
    "    print(f\"after second DMRG run: E (excluding W/H loops) = {E1:.10f}\")\n",
    "    print(\"max chi: \", max(psi.chi))\n",
    "    \n",
    "    return {'psi': psi, \n",
    "            'model': model_clean,\n",
    "            'E0': E0, 'E0_clean': E0_clean, 'E': E1, \n",
    "            'WL': WL, 'HL': HL}\n",
    "    print(\"=\"*80)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "================================================================================\n",
      "Start iDMRG for Ly=4, J_WL=-5.00, J_HL=-5.00, h=0.00\n",
      "after first DMRG run: <psi|W|psi> = -1.000\n",
      "after first DMRG run: <psi|H|psi> = -1.000\n",
      "after first DMRG run: E (including W/H loops) = -2.2500000000\n",
      "after first DMRG run: E (excluding W/H loops) = -1.0000000000\n",
      "after second DMRG run: <psi|W|psi> = -1.000\n",
      "after second DMRG run: <psi|H|psi> = -1.000\n",
      "after second DMRG run: E (excluding W/H loops) = -1.0000000000\n",
      "max chi:  8\n",
      "================================================================================\n",
      "Start iDMRG for Ly=4, J_WL=-5.00, J_HL=5.00, h=0.00\n",
      "after first DMRG run: <psi|W|psi> = -1.000\n",
      "after first DMRG run: <psi|H|psi> =  1.000\n",
      "after first DMRG run: E (including W/H loops) = -2.2500000000\n",
      "after first DMRG run: E (excluding W/H loops) = -1.0000000000\n",
      "after second DMRG run: <psi|W|psi> = -1.000\n",
      "after second DMRG run: <psi|H|psi> =  1.000\n",
      "after second DMRG run: E (excluding W/H loops) = -1.0000000000\n",
      "max chi:  8\n",
      "================================================================================\n",
      "Start iDMRG for Ly=4, J_WL=5.00, J_HL=-5.00, h=0.00\n",
      "after first DMRG run: <psi|W|psi> =  1.000\n",
      "after first DMRG run: <psi|H|psi> = -1.000\n",
      "after first DMRG run: E (including W/H loops) = -2.2500000000\n",
      "after first DMRG run: E (excluding W/H loops) = -1.0000000000\n",
      "after second DMRG run: <psi|W|psi> =  1.000\n",
      "after second DMRG run: <psi|H|psi> = -1.000\n",
      "after second DMRG run: E (excluding W/H loops) = -1.0000000000\n",
      "max chi:  8\n",
      "================================================================================\n",
      "Start iDMRG for Ly=4, J_WL=5.00, J_HL=5.00, h=0.00\n",
      "after first DMRG run: <psi|W|psi> =  1.000\n",
      "after first DMRG run: <psi|H|psi> =  1.000\n",
      "after first DMRG run: E (including W/H loops) = -2.2500000000\n",
      "after first DMRG run: E (excluding W/H loops) = -1.0000000000\n",
      "after second DMRG run: <psi|W|psi> =  1.000\n",
      "after second DMRG run: <psi|H|psi> =  1.000\n",
      "after second DMRG run: E (excluding W/H loops) = -1.0000000000\n",
      "max chi:  8\n"
     ]
    }
   ],
   "source": [
    "results_loops = {}\n",
    "for J_WL in [-5., 5.]:\n",
    "    for J_HL in [-5., 5.]:\n",
    "        results_loops[(J_WL, J_HL)] = run_DMRG(4, J_WL, J_HL)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we can see from the output, the idea worked and we indeed get the desired expectation values for the W/H loops.\n",
    "Also, we can directly see that the ground states are degenerate in energy. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "To check that we found indeed four orhogonal ground states, let's calculate mutual overlaps $\\langle \\psi_i | \\psi_j \\rangle$ for $i,j \\in range(4)$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "overlaps\n",
      "[[ 1.+0.j  0.+0.j -0.+0.j -0.+0.j]\n",
      " [ 0.+0.j  1.+0.j -0.+0.j  0.+0.j]\n",
      " [-0.+0.j -0.+0.j  1.+0.j  0.+0.j]\n",
      " [-0.+0.j  0.+0.j  0.+0.j  1.+0.j]]\n"
     ]
    }
   ],
   "source": [
    "psi_list = [res['psi'] for res in results_loops.values()]\n",
    "overlaps= [[psi_i.overlap(psi_j) for psi_j in psi_list] for psi_i in psi_list]\n",
    "print(\"overlaps\")\n",
    "print(np.array(overlaps))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The correlation length is almost vanishing, as it should be:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Correlation lengths: \n",
      "[0.02917 0.02959 0.0292  0.02962]\n"
     ]
    }
   ],
   "source": [
    "print(\"Correlation lengths: \")\n",
    "model = results_loops[(+5., +5.)]['model']\n",
    "print(np.array([psi.correlation_length() / model.lat.N_sites_per_ring for psi in psi_list]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Intermezzo**: Let's quickly check for a single case that this still works for non-zero $h = 0.1$, where the t'Hooft loop no longer commutes with the Hamiltonian:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "================================================================================\n",
      "Start iDMRG for Ly=4, J_WL=5.00, J_HL=-5.00, h=0.10\n",
      "after first DMRG run: <psi|W|psi> =  1.000\n",
      "after first DMRG run: <psi|H|psi> = -1.000\n",
      "after first DMRG run: E (including W/H loops) = -2.2516147149\n",
      "after first DMRG run: E (excluding W/H loops) = -1.0018727133\n",
      "after second DMRG run: <psi|W|psi> =  1.000\n",
      "after second DMRG run: <psi|H|psi> = -0.995\n",
      "after second DMRG run: E (excluding W/H loops) = -1.0025283892\n",
      "max chi:  64\n"
     ]
    }
   ],
   "source": [
    "res_h = run_DMRG(4, +5., -5., 0.1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The t'Hooft loop no longer is exactly 1 in the ground state, since it no longer commutes with the Hamiltonian, but it is still close to one, so it still distinguishes the degenerate states. As expected, the local perturbation of the field can not destroy the global, topological nature of the ground state (manifold)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Topological entanglement entropy\n",
    "\n",
    "Back to $h=0$. For the toric code, the half-cylinder entanglement entropy should be $(L_y-1)\\log 2$.\n",
    "Indeed, we get the correct value at bond 0:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Entanglement entropies / log(2)\n",
      "[3. 3. 3. 3.]\n"
     ]
    }
   ],
   "source": [
    "print(\"Entanglement entropies / log(2)\")\n",
    "print(np.array([psi.entanglement_entropy()[0]/np.log(2) for psi in psi_list]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Recall that bond 0 is usually a sensible cut for infinite DMRG on a common lattice.\n",
    "This is especially important if you want to extract the [topological entanglement entropy](https://arxiv.org/abs/hep-th/0510092) $\\gamma$ from a fit $S(L_y) = \\alpha L_y - \\gamma$, as corners in the cut can also contribute to $\\gamma$.\n",
    "\n",
    "Let us demonstrate the extraction of $\\gamma = \\log(2)$ for the toric code:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "================================================================================\n",
      "Start iDMRG for Ly=2, J_WL=5.00, J_HL=5.00, h=0.00\n",
      "after first DMRG run: <psi|W|psi> =  1.000\n",
      "after first DMRG run: <psi|H|psi> =  1.000\n",
      "after first DMRG run: E (including W/H loops) = -3.5000000000\n",
      "after first DMRG run: E (excluding W/H loops) = -1.0000000000\n",
      "after second DMRG run: <psi|W|psi> =  1.000\n",
      "after second DMRG run: <psi|H|psi> =  1.000\n",
      "after second DMRG run: E (excluding W/H loops) = -1.0000000000\n",
      "max chi:  4\n",
      "================================================================================\n",
      "Start iDMRG for Ly=3, J_WL=5.00, J_HL=5.00, h=0.00\n",
      "after first DMRG run: <psi|W|psi> =  1.000\n",
      "after first DMRG run: <psi|H|psi> =  1.000\n",
      "after first DMRG run: E (including W/H loops) = -2.6666666667\n",
      "after first DMRG run: E (excluding W/H loops) = -1.0000000000\n",
      "after second DMRG run: <psi|W|psi> =  1.000\n",
      "after second DMRG run: <psi|H|psi> =  1.000\n",
      "after second DMRG run: E (excluding W/H loops) = -1.0000000000\n",
      "max chi:  8\n",
      "================================================================================\n",
      "Start iDMRG for Ly=4, J_WL=5.00, J_HL=5.00, h=0.00\n",
      "after first DMRG run: <psi|W|psi> =  1.000\n",
      "after first DMRG run: <psi|H|psi> =  1.000\n",
      "after first DMRG run: E (including W/H loops) = -2.2500000000\n",
      "after first DMRG run: E (excluding W/H loops) = -1.0000000000\n",
      "after second DMRG run: <psi|W|psi> =  1.000\n",
      "after second DMRG run: <psi|H|psi> =  1.000\n",
      "after second DMRG run: E (excluding W/H loops) = -1.0000000000\n",
      "max chi:  16\n",
      "================================================================================\n",
      "Start iDMRG for Ly=5, J_WL=5.00, J_HL=5.00, h=0.00\n",
      "after first DMRG run: <psi|W|psi> =  1.000\n",
      "after first DMRG run: <psi|H|psi> =  1.000\n",
      "after first DMRG run: E (including W/H loops) = -2.0000000000\n",
      "after first DMRG run: E (excluding W/H loops) = -1.0000000000\n",
      "after second DMRG run: <psi|W|psi> =  1.000\n",
      "after second DMRG run: <psi|H|psi> =  1.000\n",
      "after second DMRG run: E (excluding W/H loops) = -1.0000000000\n",
      "max chi:  32\n",
      "================================================================================\n",
      "Start iDMRG for Ly=6, J_WL=5.00, J_HL=5.00, h=0.00\n",
      "after first DMRG run: <psi|W|psi> =  1.000\n",
      "after first DMRG run: <psi|H|psi> =  1.000\n",
      "after first DMRG run: E (including W/H loops) = -1.8333333333\n",
      "after first DMRG run: E (excluding W/H loops) = -1.0000000000\n",
      "after second DMRG run: <psi|W|psi> =  1.000\n",
      "after second DMRG run: <psi|H|psi> =  1.000\n",
      "after second DMRG run: E (excluding W/H loops) = -1.0000000000\n",
      "max chi:  64\n"
     ]
    }
   ],
   "source": [
    "model_params['order'] = 'Cstyle' # The effect doesn't appear with the \"default\" ordering for the toric code.\n",
    "# This is also a hint that you need bond 0: you want something independent of what order you choose\n",
    "# inside the MPS unit cell.\n",
    "\n",
    "results_Ly = {}\n",
    "\n",
    "for Ly in range(2, 7):\n",
    "    results_Ly[Ly] = run_DMRG(Ly, 5., 5.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA5kAAAKZCAYAAAAyIe+SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAXEQAAFxEByibzPwAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdeVxU1fsH8M9lmWERUBZBBFGRRUBAXHIHTfpa7mluuaH+svq6pmVlpaZZmrlkZZl79S33XLDMBRRLzFRcQEUUXFAMRbZhnzm/P3BGRwYEBIbl83695oU+995znzsXkWfOuedIQggQERERERERVQQDfSdAREREREREtQeLTCIiIiIiIqowLDKJiIiIiIiowrDIJCIiIiIiogrDIpOIiIiIiIgqDItMIiIiIiIiqjAsMomIiIiIiKjCsMgkIiIiIiKiCsMik4iIiIiIiCoMi0wiIiIiIiKqMCwyiYiIiIiIqMKwyCQiIiIiIqIKY6TvBOoaSZKSAJgBuKnvXIiIiIiIiHRwBpAlhHAoz8GSEKKC86GSSJKULpfLLVxdXfWdChERERERURFXr15Fbm5uhhDCsjzHsyez6t10dXX1io6O1nceRERERERERXh7eyMmJqbcIy/5TCYRERERERFVGBaZREREREREVGFYZBIREREREVGFYZFJREREREREFYZFJhEREREREVUYFplERERERERUYbiESQ0hhADXNCWi8pAkCZIk6TsNIiIiqiNYZFZjKpUKaWlpePDgAXJzc/WdDhHVYHK5HA0aNICVlRUMDDiIhYiIiCoPi8xqSgiBpKQkpKWl6TsVIqoFcnNzkZSUhJycHDg4OLBnk4iIiCoNi8xqKiMjQ1NgNmzYEJaWljA0NNRzVkRUEymVSqSnp+Pff/9FamoqzM3NYWlpqe+0iIiIqJZikVlNpaenAwCsra1hY2Oj52yIqCYzMDCAjY0NCgoKkJKSgoyMDBaZREREVGn4YE41lZWVBQCwsLDQcyZEVFuof54oFAo9Z0JERES1GYvMakgIAaVSCaBwsg4iooqg/nmiVCo5WzURERFVGhaZ1dDjv/xxcg4iqiiP/zxhkUlERESVhUUmERERERERVRgWmURERERERFRhOLssEREREVFNJgTwcNLIMjMzA/h4FlUwFpk6SJLkAGAWgN4AnAFkA4gHcEgI8Y4+cyMiIiIi0pKVBdSrV75jMzMBc/OKzYfqPA6XfYIkSR0BXAQwDUA+gN0AIgHYAHhLj6lVmrt3gU8+AV54AejYsfDrwoWF8Zpk7NixkCQJ4eHh+k6l2goPD4ckSRg7dqy+Uymz8t7fW7duYdy4cXB0dISJiQnc3d3x0UcfIScnp3ISJSIiIqrjWGQ+RpIkRwD7AMgBvCyE8BZCDBVCvCSEaAqgk14TrGDZ2cBrrwHOzsAHHwAHDgCRkYVfZ88ujE+cCNS138VrQiEmSRKaNm2q7zSqvatXryIgIADr16+HjY0N+vfvD6VSifnz56NHjx7Izc3Vd4pEREREtQ6LTG2fAagP4B0hxM4nNwoh/q76lCpHdjbw4ovA998D+fm698nPB1avBnr1KtyfqKYZN24ckpOTMWXKFJw/fx6bN2/G5cuXMXDgQBw/fhwLFy7Ud4pEREREtQ6LzIckSWoAYAiANABr9JxOpZs6FThypHT7HjkCTJtWufkQVbSTJ0/i6NGjaNiwIRYvXqyJGxkZYdWqVTA2NsbKlSuRX9ynLERERERVJPZuBlSq2rOGNYvMRzqjcJjsMQD5kiQNliRpuSRJX0uSNFmSJHs951dhkpKADRvKdsz69dXnGc3t27ejffv2MDU1hb29PUaPHo3bt28Xu39ERAQmTZoEX19fNGjQAKampvD09MS7776L1NRUrX3Hjh2L7t27AwA2btwISZI0r7lz52r2Cw0Nxbhx49CyZUtYWlrC3Nwcfn5+WLhwYZmHYMbFxWHu3Lno2LEjHBwcIJPJ4OTkhNGjRyM2NlZr3w0bNkB6OAPc9evXtfILCgoq03nv3LmDsWPHwt7eHqampggICMCmTZuK3T8mJgavvvoqGjVqBJlMhsaNG2P06NG4fPlykX0fH3KckpKCN954A40aNYJcLoePjw/WrVtX7HnKen+Ls3fvXgBA3759IZfLtbbZ29uja9euePDgAf78888yt01ERERUEQqUKqw4eAUvrojAuj/j9Z1OheHsso94P/x6F0AEgI5PbP9UkqQQIcTW0jQmSVJ0MZtcy5lfhVm7tvghssXJzy887v33Kyen0vrqq68wefJkGBoaIjAwELa2tjh48CA6dOgAPz8/nce8/fbbiIqKgo+Pj+Y5vNOnT2PRokXYu3cvIiMjUe/hjGxdunRBUlIS9u/fD1dXV3Tp0kXTjr+/v+bP48ePh0KhgLe3N1q1aoX09HT8/fffmD17Ng4dOoQ//vgDhoaGpbqmNWvWYNGiRfDy8kLbtm1hYmKCmJgY/PDDD9i1axciIiLg6+sLAGjRogXGjBmDjRs3wtzcHIMHD9a04+npWer3MSUlBR06dEBubi6CgoLw4MEDhIWFYcyYMYiPj8ecOXO09j906BD69u2L7OxsBAQEICgoCJcuXcIPP/yAnTt3Yt++fejatWuR86SmpqJjx45IS0tD+/btkZmZiaNHj2L8+PFQqVSYMGGC1v7lub/FOXv2LAAgICBA5/aAgAAcPnwYZ8+eLXOBTkRERPSs4u8pMH1zFKJuFnZ6LP79Mrq62cHDwULPmVUAIQRfQgCFz2MKFM4omwlgHABbAE0BfPFwWy4A31K2F13MK8fLy0uURKlUipiYGBETEyOUSmWJ+5ZHcLAQhQsqle0VHFzhqZRJfHy8kMvlQi6Xi7CwME1coVCI4OBg8fAeaW0TQojQ0FCRkpKiFcvJyRGvvfaaACDmzZuntS0sLEwAEGPGjCk2l507d4rMzEytWHp6uujTp48AIDZu3Fjq6zp+/LiIi4srEl+3bp0AILp3715kGwDh4uJS6nOoqa8NgAgODta6hr///lvUq1dPGBgYiDNnzmjimZmZwt7eXgAQq1at0mpv6dKlAoBwcnISOTk5Os8zaNAgrfP8+uuvAoBo0qSJVlvlvb/Fad26tQAgdu3apXP78uXLBQDx1ltvlaq92qCyf7YQEZGeZGaW75c7oPBYqlIqlUr8GJkgPD/4TbjM2qt5dfr0kDh9PeXpDVQBLy8vASBalLO24nDZR9TdTkYA3hJCrBNC3BNCJAghZgDYBkAGoFTrZIrCmWmLvABcrZz0Sy8jo2qPqyjr1q1Dbm4uRo8erdXzZGZmhpUrV2qGkT7ppZdeQoMGDbRicrkcy5cvh5GREXbt2lXmXAYMGADzJ9aUsrCwwLJlywCgTG126NABrq5FO7hDQkLQuXNnhIeHIy0trcw5lkSSJKxcuVLrGtq1a4f//ve/UKlUWLVqlSa+ZcsW3L17F127dsXrr7+u1c706dPRpk0b3Lp1Czt3FpkrC5aWlli9erXWefr3749WrVrhxo0bSEhI0MTLe3+Lk5mZqTleF3VO6v2IiIiIKtu/GTkYt+EkZu+8gOx8pSY+uI0Tfp/WFa2bNCjh6JqDw2UfUZdQKgAbdWxfB2AwgKCqSqiyWJSzB768x1WUY8eOAQCGDBlSZJuHhwdat26N06dP6zw2MTERe/bswaVLl5Ceng6VSgUAkMlkuHLlSrnyuXLlCvbt24e4uDgoFAqoVCp1L3aZ28zMzMSePXsQFRWFlJQUzWQ0d+7cgRBCsxRHRWndujU8PDyKxIcPH45FixZp3mug8JlWAHj11Vd1tjVy5EicOnUKERERGDZsmNa2tm3bwtrausgx7u7uOH/+PO7cuaNZiuVZ7q8u6ntRXHGq3k5ERERUFX6/kIT3dpzDg6xHz601MDPGpy+3Qi+fRnrMrOKxyHwk4eHXJCGErplb1NsbVkk2lSgwsHAtzLLS92Nr6slfmjRponN7kyZNdBYhS5cuxXvvvYe8vLwKyUMIgZkzZ2LZsmXFFioZZej2PXz4MIYNG4bk5ORi9ylLe6Xh4uKiM64u+B6faEf95+LW5dR1jJqTk5POY9TPwD4+SVJ5729xLB5+KqJQKHRuz8rK0sqFiIiIqDJk5OTj4z0x2Hrqlla8u4cdFg32RUMLEz1lVnk4XPaRMw+/NpB0d33YPPxa48fWjR8PGBuX7Rhj48Lj9OlpPVO6REZGYsaMGTA1NcWGDRuQkJCAnJwczXjxRo3K/qnR5s2bsXTpUjRu3Bjbtm1DYmIi8vLyIITQFE2l7SXLzMzEkCFDkJycjA8//BAxMTFavaLDhw8vU3uV6Wnvu67tZblX5bm/JVEXq7du3dK5XR0vrqglIiIielZ/x6fgxRURWgWmqbEhPhnog3Vj29XKAhNgkakhhDgPIB6AKYDndOwS9PBr6btSqikHB2Ds2LIdExIC2Ot5ERdHR0cAhUt36HLjxo0iMfVzggsWLMCYMWPg4uKiWc4iOzsbSUlJZc5D3eaqVaswaNAgODo6wvhh1X7t2rUytRUREYH79+9j0KBB+Pjjj9GyZUuYmZlpCq2ytldaxb2H6rj6vX78z/HxuqfVVh9TnoL9ceW5vyVRz0ZbXO+nOq6euZeIiIioouQWKPHZb5cwdPVx3HqQrYn7O9fHvqld8epzLhX2wXp1xCJT26KHX7+UJMlWHZQkqQ2AGQ//+m2VZ1UJVqwoHDZbGoGBhfvrm3o5ka1bi64iExsbi6ioqCLxBw8eAACcnZ2LbNu6davOHkKZTAYAKCgo0JlHSW1u2bKluPTL3FZcXFyxBZKxsXGx+ZVGVFRUkTU4AeDnn38GAHTu3FkTUy9N8tNPP+lsSx3XtYRJWZTn/pakd+/eAIA9e/YUWbv07t27iIiIgJWVldYyNURERETP6nJSBgZ8/Re+PXIV6l81DQ0kvBXsjm2vd0QzW/OSG6gFWGRq+x7AVgDtAFyWJGmPJElhAP4CUB/A90KIbfpMsKKYmgK//Qa89lrxQ2eNjQu3//47YFINevJDQkIgk8mwadMmzWQ0QGGP5NSpUzWT+TzO3d0dALB27VrNZDoAEBMTg1mzZuk8j7pH7fLlyzq3q9tcvXq1VpEaERGBzz//vEzXpG5rx44dWs9kpqamYvz48Vo5P5nj3bt3kZqaWqbzqalUKkyZMkXzXCIAnDp1Cl9//TUMDAwwceJETXzIkCGwt7dHREQEVq9erdXOl19+iZMnT8LJyQkDBw4sVy5q5bm/JWnfvj06d+6Mf//9V+teFxQU4M0330R+fj4mT56s6YUmIiKqsczMgMzM8r2KmYWdyk6lElgTcQ19Vx7DxTvpmnhzW3PseKMTpjzvBiPDulF+1Y2rLCUhhArAMAD/BXAdQA8UFpz/ABgthHhNj+lVOFNT4LvvgJs3gU8+AYKDgQ4dCr9+8klh/LvvqkeBCQDNmzfHokWLkJOTg+7du6Nnz54YNmwYWrRogQsXLqBPnz5FjgkJCYGDgwP27NkDDw8PDB06FMHBwfD390fXrl11ToDTtGlT+Pr64p9//kH79u0REhKCCRMmYPfu3QCAKVOmwNzcHN988w18fHwwfPhwdOvWDYGBgUWW+Hiatm3bIjg4GDdu3IC7uzsGDhyIgQMHolmzZrh9+zb69++v87h+/fqhoKAAAQEBGDlyJCZMmFCmArdPnz64ePEiXF1dMXToUPTq1QsdO3ZERkYG3n//fbRp00azr7m5OX766SeYmppi4sSJaNu2LUaMGIGAgABMnToV5ubm+N///qcZhlxe5bm/T7N+/XrY2NhgxYoV8PX1xbBhw+Dh4YEdO3bgueeew+zZs58pZyIiompBkgBz8/K9avGQzaqUmJqNV9ecwILQi8hTPvpgfHRHF4RO6Qo/5/p6zE4PyrvAJl/lewGI9vLyEiXhgukl27Jli2jTpo2Qy+XC1tZWjBgxQty6dUuMGTNGABBhYWFa+9+8eVOMGDFCNG7cWJiYmIiWLVuKTz/9VBQUFAgXFxdR+M9A25UrV8SAAQOEjY2NMDAwEADEnDlzNNtjYmJE3759RcOGDYWZmZlo3bq1WL16tRBCCADCxcWl1NeTlZUlZs+eLdzc3IRcLhfOzs7i9ddfF/fu3Sv2mjIzM8WkSZOEs7OzMDIyEgBEYGDgU88VFhYmAIgxY8aIxMREMXLkSGFnZyfkcrnw8/MT69evL/bYCxcuiOHDhwt7e3thbGwsGjVqJEaOHCkuXbpU4nl0Ke66hCj7/X2aGzduiLFjxwoHBwchk8mEq6ur+OCDD0RWVlaZ2qkN+LOFiIioYqlUKrHz9C3hM+d34TJrr+bVdsEBEXbprr7TKzcvLy8BIFqUs+aRRDWYtbIukSQp2svLyys6OrrYfVQqlWaopoeHBwwM2OFMRM+OP1uIiIgqTmpWHmb/egGh5+5oxV/0ccDCga3QwFymp8yenbe3N2JiYmKEEN7lOZ7rZBIREREREZXB0dhkvL3tLO6mP5pc0EJuhHn9vTGwdeNaPXNsabDIJCIiIiIiKoXsPCU+++0iNh7XXnLtuWbW+GKIH5wacCIlgEUmERERERHRU527lYrpm6NwNVmhickMDfD2fzwwvkszGBjU7d7Lx7HIJCIiIiIiKkaBUoVV4Vex4tAVFKgezWfj6WCB5cP84elgqcfsqicWmURERERERDok3FNg+pYonLnxaG1ySQJe69YcbwW7Q25kqMfsqi8WmURERERERI8RQuDnv29i/t4YZOcrNfHG9U3xxRA/dGhuo8fsqj8WmURERERERA8lZ+Ti3e3ncOjSv1rxQQFOmNPPC5YmxnrKrOZgkUlERERERARgf3QS3ttxHimKPE2sgZkxFg5shRdbNdJjZjULi0wiIiIiIqrTMnML8PGeaGz555ZWPMjDDosH+aKhpYmeMquZWGQSEREREVGddTIhBW9ticLNlGxNzMTYALN7e2Hkc00gSVyapKxYZBIRERERUZ2TV6DCsoOx+PbIVYhHK5PAz8kKy4b6o7ldPf0lV8OxyCQiIiIiojol9m4Gpv0ShZg76ZqYoYGEyT1a4L/dW8DY0ECP2dV8LDKJiIiIiKhOUKkE1v0Zj8X7LyOvQKWJN7M1x7Kh/vB3rq/H7GoPluh1jRCAQlG+1+PjCKqhsWPHQpIkhIeH6zsVqqEkSULTpk31nUaZbdiwAZIkYe7cufpOhYiIqNq6nZqNkWtPYEHoRa0Cc2SHJgid0oUFZgVikVnXZGUB9eqV75WVpe/sq0R4eDgkScLYsWP1nUqxamoxRERERKQPu6IS8Z/lR/HX1fuamJ2FHOtD2mHBgFYwk3GAZ0Xiu0lERERERLVSalYePtwVjT1nb2vFe3k7YOHLrWBtLtNTZrUbi0wiIiIiIqp1Iq4k4+2t55CUnqOJ1ZMbYV4/b7wc0JhLk1QiDpelGmf79u1o3749TE1NYW9vj9GjR+P27dvF7h8REYFJkybB19cXDRo0gKmpKTw9PfHuu+8iNTVVa9+xY8eie/fuAICNGzdCkiTN6/Hn3UJDQzFu3Di0bNkSlpaWMDc3h5+fHxYuXIjc3NwyXU9cXBzmzp2Ljh07wsHBATKZDE5OThg9ejRiY2O19lU/ewcA169f18ovKCioxPPEx8dDkiTN9T3Ox8cHkiQhJCREK65UKmFlZYWGDRtCPHwm9/HhxElJSZgwYQKcnJxgZGSE5cuXa47NysrC/Pnz4ePjA1NTU1hZWaFbt2745ZdfdObXtGlTzbWtWbMGvr6+MDU1hYODAyZOnFjkXqldu3YNQ4YMgbW1NerVq4cuXbpg//79zzTsOS8vD3PmzIGrqytMTEzQvHlzfPTRR8jJydG5//379/H222/Dzc0NJiYmsLa2Rq9evfDHH3/o3F893FmpVGLx4sVwd3eHXC6Hs7MzZs2aVez30Llz59CnTx9YWVnBysoKwcHBOH78eJmvj4iIqDbLyVdi7u5ojFr7t1aB2b6ZNX6b2hWD2jixwKxk7MmkGuWrr77C5MmTYWhoiMDAQNja2uLgwYPo0KED/Pz8dB7z9ttvIyoqCj4+PujRowdyc3Nx+vRpLFq0CHv37kVkZCTq1StcB6lLly5ISkrC/v374erqii5dumja8ff31/x5/PjxUCgU8Pb2RqtWrZCeno6///4bs2fPxqFDh/DHH3/A0NCwVNe0Zs0aLFq0CF5eXmjbti1MTEwQExODH374Abt27UJERAR8fX0BAC1atMCYMWOwceNGmJubY/DgwZp2PD09SzxPs2bN4OLiguPHjyMnJwcmJiYAgHv37iEmJgYAEBYWpnXMmTNnkJ6ejkGDBhX5YZycnIx27dqhoKAAXbp0QU5ODszMzAAAGRkZ6N69O06dOgU7Ozv06dMHCoUChw8fRkREBCIjI7UK0se98847WLFiBdq1a4devXrhr7/+wurVq3Hx4kUcOXJEK48rV66gU6dOuHfvHjw8PNC6dWskJCTgpZdewptvvlmat78IIQQGDx6MgwcP4vnnn4e/vz8OHTqE+fPn46+//sL+/fu17m1iYiK6deuGa9euoUmTJhgwYACSk5Nx8OBB7N+/H0uXLsX06dN1nuvVV1/F3r170b59e3h4eCAiIgKLFy9GYmIifvzxR619T5w4gR49eiArKwv+/v7w9PTEhQsXEBgYWK2fHyYiIqpK52+lYdrmM7iarNDEZIYGmPkfd4zv0hyGBiwuq4QQgq8qfAGI9vLyEiVRKpUiJiZGxMTECKVSWeK+ZZaZKUThPLFlf2VmVmwuZRQfHy/kcrmQy+UiLCxME1coFCI4OFgAEAC0tgkhRGhoqEhJSdGK5eTkiNdee00AEPPmzdPaFhYWJgCIMWPGFJvLzp07ReYT70d6erro06ePACA2btxY6us6fvy4iIuLKxJft26dACC6d+9eZBsA4eLiUupzqI0ePbrIe7R161YBQHh7ewsAIj4+XrPt888/FwDEypUrNTH1+wNADBw4UGRnZxc5z6RJkwQA0bNnT5GRkaGJX7x4UTRs2FAAEKGhoVrHuLi4CACiUaNG4syZM5p4cnKyaNGihQAgDh06pHXM888/LwCIyZMna/1b2bRpkybHku7jk9THODk5iatXr2ri//77r/Dx8REAxIoVK7SOUd/zUaNGiby8PE08IiJCmJmZCUNDQ3H27Fmd52nZsqXW+33t2jXRoEEDAUDre0KpVApPT08BQHz66adabX3wwQea9ubMmVPi9VXqzxYiIiI9yi9QipWHYoXre6HCZdZezes/y46ImNtp+k6vxvHy8hIAokU5ax4Ol6UaY926dcjNzcXo0aO1hoaamZlh5cqVxQ57eOmll9CgQQOtmFwux/Lly2FkZIRdu3aVOZcBAwbA3NxcK2ZhYYFly5YBQJna7NChA1xdXYvEQ0JC0LlzZ4SHhyMtLa3MOeoSGBgIAFrLvKh7Bz/88EOd2wDoHIorl8uxcuVKTY+omkKhwNq1a2FgYIBvvvlG00sMFPa2fvDBBwCAL7/8UmeO8+fP1+o1trW1xRtvvAEAOHr0qCYeFxeHQ4cOwdraGp9++ikMDB79OBs1ahS6detW7PvwNB999BGaN2+u+budnR0+//xzAMDXX3+tiV+7dg179+6FpaUlvvzySxgbG2u2denSBa+//jqUSiW++eYbnedZuXKl1izBzZo1w8iRIwEUDvNWCw8Px6VLl+Du7o5Zs2ZptTFnzhw0adKk3NdKRERU012/r8CQ745jyR+xKFAVPt4jScBr3Zrj1/92RstGlnrOsO7hcFmqMY4dOwYAGDJkSJFt6qGSp0+f1nlsYmIi9uzZg0uXLiE9PR0qVeHaSDKZDFeuXClXPleuXMG+ffsQFxcHhUIBlUqleW6xrG1mZmZiz549iIqKQkpKCvLz8wEAd+7cgRACV69eRUBAQLnyfJy6WHy8kAwPD4ePjw/69esHmUyG8PBwjB07FiqVCseOHYOtrS28vb2LtBUQEIDGjRsXiZ86dQrZ2dno0KED3NzcimwfNWoUpkyZgj///BNCiCIfDrzwwgtFjnF3dwdQ+H6o/fXXXwAKP0R4suAHgFdeeUWrKC2LYcOGFYn16tULDRo0QGxsLJKTk2FnZ6f5nnzppZdQv37RtbVGjRqFpUuXahWMasbGxjqLd13Xqj7PK6+8UuT9MjIywuDBg7F06dLSXyAREVEtIITALydvYv7eGGTlKTXxxvVNseQVP3R0tdFjdnUbi0yqMdST+xTXa9OkSROdRebSpUvx3nvvIS8vr0LyEEJg5syZWLZsmaaofFJGRkap2zt8+DCGDRuG5OTkYvcpS3slad68OZydnREZGYmcnBxkZmYiOjoakydPhqmpKdq1a6cpQM+cOYPU1FS8/PLLOnuJi7sP6vtU3Dqe9evXh5WVFdLS0pCeng4rKyut7U5OTkWOUfeGPj4hjvo8zs7OOs9T3t69Bg0awMLCQuc2FxcXPHjwALdv34adnd1Tr1Ud1zUxVaNGjXQ+t1vStZb0vU9ERFSXJGfk4t3t53Do0r9a8ZdbN8bc/t6wNDEu5kiqChwuSzWGuqAry2xgkZGRmDFjBkxNTbFhwwYkJCQgJydHM168UaNGZc5j8+bNWLp0KRo3boxt27YhMTEReXl5EEJoCoPiis8nZWZmYsiQIUhOTsaHH36ImJgYrV7R4cOHl6m90ggMDERubi4iIyNx9OhRCCE0PWpBQUG4fv06EhISShwqC6DIMNknleY+6dqnrLO9Fbd/Rb5nT2uzuBzU8We9zvJ87xMREdVWB2Luotfyo1oFZn0zY3w9IgBLh/qzwKwGWGRSjeHo6AigcOkOXW7cuFEktnPnTgDAggULMGbMGLi4uEAulwMAsrOzkZSUVOY81G2uWrUKgwYNgqOjo+ZZvGvXrpWprYiICNy/fx+DBg3Cxx9/jJYtW8LMzExTTJS1vdJ4fMisepkP9bOaT24DHj3HWVrq+xQfH69ze1paGtLS0mBubl5sj2FpqD8g0HXfAeDmzZvlavfBgwfF9hyrz6U+99OuNSEhQWv/8irP9z4REVFtk5lbgFnbzuH/Nv2D+4pHI9S6udth/7Ru6O37bP/fUsVhkaYiiXUAACAASURBVEk1hno5ka1btxbZFhsbi6ioqCLxBw8eANA9pHLr1q06e6ZkMhkAoKCgQGceJbW5ZcuW4tIvc1txcXHFPmNqbGxcbH5P82Qh6evrC2trawBAp06dIJPJNEuN2NjYoFWrVmVqv02bNjA1NcXff/+t89lU9dIcXbp0eaaeuU6dOgEA9u3bh6ysrCLbt23bVu62N2/eXCS2f/9+PHjwAG5ubmjYsCGAR9+ToaGhOtfxVF9r165dy53L4+fZvn17ke/ZgoICbN++/ZnaJyIiqu7+SUjBSysisPmfRx8imxgbYH5/b2wMaQd7y5JHWFHVYpFJNUZISAhkMhk2bdqkNZFKdnY2pk6dqpnM53HqSVTWrl2rmUwHAGJiYorM0qmm7jW6fPmyzu3qNlevXq31C39ERIRmBtLSUre1Y8cOrWcyU1NTMX78eK2cn8zx7t27Ogubp3F1dYWTkxOOHz+OCxcuFJmpt127dti6dStSU1PRrVu3MheC5ubmGDduHFQqFf773/9CoXi0TlVsbCwWLFgAAJg8eXKZc3+cm5sbunfvjpSUFLz//vta9+Knn37SDPctj48//ljTCwkUriX6zjvvAIDW+pvNmzdH7969kZGRgalTp2rdr+PHj2PVqlUwNDQs95qdat27d4e7uzsuXbqEJUuWaG1bsGBBsT2cRERENV1egQqf77+EId8dx42URx8q+zlZIXRKV4zq2JSPk1RDLDKpxmjevDkWLVqEnJwcdO/eHT179sSwYcPQokULXLhwAX369ClyTEhICBwcHLBnzx54eHhg6NChCA4Ohr+/P7p27QoXF5cixzRt2hS+vr74559/0L59e4SEhGDChAnYvXs3AGDKlCkwNzfHN998Ax8fHwwfPhzdunVDYGAgXn/99TJdU9u2bREcHIwbN27A3d0dAwcOxMCBA9GsWTPcvn0b/fv313lcv379UFBQgICAAIwcORITJkwoU4EbGBioeY70yWcug4KCkJOTo9mvPD799FO0adMGBw4cQPPmzTFkyBD07t0bfn5+SEpKwpQpU9C7d+9ytf24VatWwdraGitWrICXlxdGjBiBTp06YdSoUZrCTt0zXVpNmjSBn58fvL290a9fPwwaNAhubm44d+4cunfvjkmTJmnt/91336FZs2bYtGkT3NzcMHz4cPTs2RNdu3aFQqHA4sWL4evr+0zXaWBggA0bNsDU1BTvvPMOAgICMGLECPj6+mLhwoWYMGHCM7VPRERUHV25m4GB3/yJr8Ou4uHKJDA0kDD1eTdse6MTXO3qldwA6Q2LzLrGzAzIzCzfy8xM39lj2rRp2LJlC/z9/XHs2DEcOnQIQUFBiIyMhI1N0WmqbWxscPLkSYwYMQJ5eXnYvXs3EhMT8fHHH+Pnn38u9jzbt2/HgAEDcO3aNWzatAlr167VDF11d3fHyZMn0bdvX9y7dw+7d+9GZmYmvvvuuzL3ZAKFa2rOnj0bdnZ2+O2333Dq1CkMGzYMkZGROpfFAAqLuEmTJqGgoACbN2/G2rVrERoaWupzqgtLSZKKrCf5eNFZ3KQ/T2NhYYEjR45g3rx5sLW1xe7duxEREYG2bdvif//7H1asWFGudp/k4eGBEydOYPDgwUhKStKsT7p3717NcFpd3xclkSQJ27Ztw7Rp03D+/Hns3bsXVlZWmD17NkJDQ2FkpD0pd+PGjXHy5EnMmDEDRkZG2LFjB06dOoXnn38e+/fvx1tvvVUh19qxY0f89ddfePHFFxEXF4e9e/fCzs4OYWFh6Ny5c4Wcg4iIqDpQqQTWHYtH75XHEH07XRNvamOGba93xPRgdxgbsoypzqTKmIGRiidJUrSXl5dXdHR0sfuoVCrNUE0PDw+tReaJqHTeeOMNfPvtt/jll18wdOhQfadTLfBnCxERVXd30rIxc+tZ/Bl3Xyv+6nNNMLt3S5jJuAJjVfD29kZMTEyMEKLoYumlwLtERDVWTk4O4uPj0bJlS6349u3bsWbNGlhZWVXIsFwiIiKqfLuiEvHhrxeQnvNockPbenIsHtwKPTzt9ZgZlRWLTCKqsVJTU+Hl5QVPT0+4ubnB2NgYFy9exMWLF2FgYIBvvvkG9erxeQ0iIqLqLC0rHx/suoA9Z29rxf/jbY+FA1vBpp5cT5lRebHIJKIay8rKCjNmzMCBAwdw7NgxZGRkwNraGv3798fMmTM1S38QERFR9XTsyj3M3HoWSek5mlg9uRHm9PXC4DZOnDm2hmKRSUQ1lqmpaZElPYiIiKj6y8lXYtHvl7D+zwStePum1vhiiB+crfU/4SSVH4tMIiIiIiKqMhcS0zBtcxTi/s3UxIwNJcx4wQP/17U5DA3Ye1nTscgkIiIiIqJKp1QJfHvkKpYdiEWB6tEKFx72Flg21B9ejpZ6zI4qEotMIiIiIiKqVDfuZ2H6liicuv5AKz6hSzPM/I8HTIwN9ZQZVQYWmUREREREVCmEENjyz018vCcGijylJu5oZYIlr/ihUwtbPWZHlYVFJhERERERVbh7mbl4d/t5HLx4Vys+sHVjzO3nDStTYz1lRpWNRSYREREREVWoAzF38e72c7ivyNPErEyN8clAH/TxddRjZlQVWGQSEREREVGFyMwtwIK9Mfjl5E2teFc3W3w+2A8OViZ6yoyqEotMIiIiIiJ6Zqeup2D65rO4kZKlicmNDDC7d0uM6uACSeLSJHUFi0wiIiIiIiq3vAIVVhyKxarwq3hsZRL4Ollh6RB/tGhYT3/JkV6wyCQiIiIionKJ+zcD0zZH4UJiuiZmIAGTurfA5OfdYGxooMfsSF9YZBIRERERUZmoVAIbjyfgs98uIbdApYm72Jhh2VB/BDRpoL/kSO/40QIR6TR27FhIkoTw8PAKa/P111+HJEn47LPPKqzNynbq1Cl89tlnePnll9G4cWNIkgQTE05aQEREddedtGyMXvc35u2J0Sowh7dvgn1TurLAJPZkEtUECQkJaNasGQIDAyu06KtqZ8+eBQD4+fnpOZPSmz9/Pnbt2qXvNIiIiKqFPWdvY/bO80jPKdDEbOvJsGiQL55vaa/HzKg6YZFJRFVCCIHz588DAHx9ffWcTel17NgRfn5+aNeuHdq1awcHBwd9p0RERFTl0rLy8dHuC9gVdVsr/oKXPT59uRVs6sn1lBlVRywyiahKxMXFQaFQwMbGBo0bN9Z3OqU2a9YsfadARESkV3/G3cPMrWdxJy1HEzOXGWJOP2+80saJS5NQEXwmk2qkhIQETJw4EU2bNoVcLoednR0GDx6Mc+fOae03dOhQSJKks1C4dOkSzMzMYGlpiWvXrmm1LUkSgoKCkJ6ejqlTp8LZ2RkmJiZo2bIlli1bBpVKVaS9suT1uBs3bmDSpElwc3ODiYkJbGxs0L59eyxcuBDZ2dmYO3cumjVrBgA4cuQIJEnSvMaOHfvM59++fTvat28PU1NT2NvbY/To0bh9+3ax+5eXeqhsTerFJCIiqsty8pX4eE8MXl1zQqvAbNe0AX6f1g1D2jqzwCSdWGRSjXPs2DH4+flh9erVqFevHvr16wc3Nzfs2LEDHTp0QFhYmGbfb7/9Fk5OTliyZInWs4z5+fkYMWIEsrOz8eWXX6J58+ZFzpObm4sePXpg06ZNaN++PYKDg3H9+nW89dZbGD9+/DPlpXb06FH4+vri66+/hkqlQv/+/dGxY0fcu3cPs2fPxt27d+Hv749BgwYBAOzt7TFmzBjNq0uXLs90/q+++gqDBw/G6dOn0alTJwQFBeHgwYPo0KED7t+/r/P9nzt3rs4C92lq4vOYREREddWFxDT0XXkM6/6M18SMDSXM6uWJX17rCGdrMz1mR9WeEIKvKnwBiPby8hIlUSqVIiYmRsTExAilUlnivnVNWlqacHBwEMbGxmLr1q1a2w4cOCBkMplo3LixyM3N1cQPHz4sDAwMhLOzs3jw4IEQQoh33nlHABCDBw8uco74+HgBQAAQvr6+Ijk5WbMtLi5OODo6CgBi165dz5RXSkqKsLOzEwDEsmXLhEql0jruyJEjIjU1VSunwMDACntf4uPjhVwuF3K5XISFhWniCoVCBAcHa96Dx7cJIcScOXMEADFmzBiduRSnb9++AoBYv359mY6rbgAIuVyu7zTKhT9biIjoaQqUKvHV4SuixfuhwmXWXs2r5xfh4vytVH2nR1XEy8tLAIgW5ax52JNZAwkhkJadX+NeorDIfibr1q1DUlISZs6cicGDB2tt69mzJ958800kJiZi7969mnj37t0xY8YM3Lx5E2+++SbCw8OxZMkSODo64rvvvivxfEuWLIGtra3m766urvjwww8BAF9//fUz5fX9998jOTkZffr0wbRp04oMN+nWrRusrKwq7X1Zt24dcnNzMXr0aAQFBWniZmZmWLlyZbHDX2xtbeHh4YFGjRqVKjc1DpclIiKq3m7cz8LQ747j8/2Xka989Hvb+C7NsGdyF/g0Lt3vJUSc+KcGSs8pgN+8P/SdRpmdnfMCrEyNn6mNAwcOAAAGDBigc3uXLl2wfPlynDx5Ei+//LImvmDBAhw8eBA///wzQkNDIYTAxo0bYW1tXey5rK2tERwcXCQ+YsQIvPHGG/jrr78ghIAkSeXK6+DBgwCAiRMnluLKS1ae8x87dgwAMGTIkCL7e3h4oHXr1jh9+nSRbZMmTcKkSZPKlF9qaipu3LgBIyMjeHt7l+nYshg8eDAuXLhQpmPUw6GJiIjqKiEEtv5zC/P2REORp9TEG1mZYMkrfujcwraEo4mKYpFJNUpCQgIA4Lnnnitxv3v37mn9XSaTYcOGDfDz80N6ejomT56Mnj17ltiGi4uLzrilpSXq16+P1NRUpKenw8rKqlx53bx5E0Bh7+izKs/51ZP7NGnSROe+TZo00Vlkloe6F9Pd3R1yeeVNcZ6QkIDLly+X6ZisrKxKyoaIqBYQAijvz0kzM4CTwlR79zJz8d6O8zgQc1crPsDfEfP6+zxzBwHVTSwyHyNJUjiAwBJ2eVEI8XsVpUM6KJWFn6698sorMDMr/oFzXcXW5s2bNX8+c+YMVCoVDAzKN2L8yaG/z5JXRczKVp7zq6+hKmaFq6pJf/75559KbZ+IqM7JygLq1SvfsZmZgLl5xeZDpXb3LrBmDXDkCJCRAVhYAEFBwPjxgL194T4HY+7i3R3ncC8zT3OcpYkRPhnYCn39HPWTONUKLDJ12w4gU0c8saoT0cXSxAhn57yg7zTKzNLk2b/dnJyccPnyZXzwwQdlerYvIiICn332GRwdHeHp6YnDhw/js88+w/vvv1/sMTdu3NAZT09PR1paGszNzWFpaVnuvJydnXHp0iXExcXB09Oz1NeiS3nO7+joiNjYWFy/fh1ubm5Fthd3/eVR1iKzZcuW6NOnDz7//HNNLC8vDz4+Pnj11VcxZ86cCsuNiIioNsnOBqZOBTZsAPLztbcdOADMnQuMCimA9fMx2Hr6ptb2rm62+HywHxysTKosX6qdOPGPbjOFEGN1vM7rOzGgsOfJytS4xr0qosdMPcT1119/LfUx6enpGDVqFIQQWL9+PX788UfY2Nhg7ty5JfZ83b9/X/Pc5ON+/vlnAECnTp0011SevNTHrF69+qn7ymQyAEBBQUGJbZXl/OrlT7Zu3VpkW2xsLKKiokrd1tOUddKfzp07IzIyUiu2ZMkSFBQU6FzzlIiIiAoLzBdfBL7/vmiBqSbZPcDviNAqMOVGBpjb1wsbQ9qzwKQKwSKTapSJEyfCzs4OCxcuxPr164sMW1UoFNi0aRNu3bqlib355pu4fv06Jk+ejBdeeAGNGjXC999/j/z8fLz66qslPpP39ttva60XGR8fj/nz52vafZa8JkyYAFtbW+zZswdfffVVkWMiIiKQlpYGoHBGV2NjY1y9elUzNPZZ35eQkBDIZDJs2rQJERERmnh2djamTp0KlUql8z356quv4Onpiffee0/3m/YEpVKJ6OhoAKXvyezcuTNOnTqlKapv3ryJTz75BMuWLYOJCf/zIyIi0mXq1MLhsToZqGDV9TIcXv0Lxg0e/e7j09gSoVO6YGznZjAw4DO0VDGkilhWorZ47JnMZkKIhEo6R7SXl5eX+pduXVQqlWbyEg8Pj3I/N1hb/fnnn+jXrx9SUlLg4uICHx8fyOVy3LhxAxcvXoRCocCZM2fg7++PX375BcOHD4e3tzf++ecfrQJl3LhxWL9+PSZOnIhvv/1WE09ISECzZs3QoUMH5OXlIT4+Hj169EBeXh4OHTqErKwsjBw5Ej/88EO581ILCwtD//79kZGRAVdXV7Rp0wZZWVmIjo5GfHw84uPj0bRpUwBAv379sGfPHnh7eyMgIAAymQydO3dGSEhIuc+/fPlyTJ8+HYaGhggKCoKtrS0iIiJgYGAAf39/7N27F2FhYVpLnMydOxfz5s3DmDFjsGHDhqfer5iYGM2MsoGBxT/yvHjxYs0sr7GxsfDw8MCpU6cQEBCAwYMHQ6FQ4Lfffnvq+SpaaGio5oMFADhx4gQkSdKakfbDDz9E7969qzy3suLPFiIqM4WCz2TWEElJQJMmunswjWwyYNv7LOSN0jQxoQIy/26Bs/9zg5Mj/z8gbd7e3oiJiYkRQpRrWQA+k6nbeEmSbACoAMQC+FUIUXEPqNEz6dy5M86fP4+lS5ciNDQUhw8fhqGhIRwdHdGnTx+8/PLL8PLyws2bN/HGG29AJpPhxx9/LNID9uWXX+LIkSP47rvv8NJLL6Ffv35a2+VyOX7//Xe8//77+PXXX3Hv3j00a9YM//d//4dp06aVO6/Hde/eHVFRUVi0aBH279+PX3/9FZaWlnB1dcVrr70GBwcHzb5r1qzBzJkzceDAAfzvf/+DUqlEQUGBpsgsz/mnTZuGxo0bY9GiRTh27BgsLCzwwgsvYPHixZg9e/Yz3Sc19VBZADhS7Mer0Fp3093dHba2toiMjMS9e/ewZ8+eMi9NUlGSk5Nx4sQJrZgQQiuWnJxc1WkRERFpWbtWV4EpYBGQgPpBl2Bg/GiEUv4DM9wP9UNuojU2bQBKmKKCqFzYk/mYEmaXzQcwXwgxX8e24toqrqvS1cvLS86ezOpL3ZMZGBiI8PBwfadTZ/Xv3x9mZmY4c+YMBg4ciE8//VTfKdV4/NlCRGXGnswa44UXCif2UTOslwObl87CtJn2sm4ZUc54cNgLIr+wryk4GPij5i2/TpXsWXsy+RuGtqMARgFwBWAGwAPAbAAFAD6WJGmqHnMjqlM6d+6MzZs3Q6FQ4IMPPtB3OkRERNVaRsajP5t53kajcUe1CkylQoZ/t7VFyn5fTYH55HFEFYXDZR8jhPjoiVAsgIWSJP0DYD+AeZIkrRZCZJeiLZ1V/8MeTi9d24jokYCAAAgh8MUXX8Ccn4QTERGVyMICkOT5sA6+gHret7W2ZV2xx/3fW0GVJdd5HFFFY5FZCkKIPx4Wmm0BdAAQpueUiGq9jRs3onv37hgyZIi+UyEiIqr2WnS+h5gWZ2FkmaOJqXINkXLIG4rzTgB0zxz72Nx+RBWGRWbpXUFhkdnoaTtSzda0adMiS4BQ1VAqlbh//z62b9+OnTt34syZM/pOiYiIqFrLyVfi8/2XsS87HkaWj8VvNcD9vf4oSDMr9lhjY2D8+CpIkuocFpml1+Dh10y9ZkFUix05cgQ9e/ZE8+bNsXXrVri5uek7JSIiomor+nYapm+OQuzdR7+eCqWE1GPuSD/hCoiS170MCQHs7Ss7S6qLWGSWgiRJdgC6PvzraX3mQlSb9ejRAyqV6uk7EhER1WFKlcDqo9ew9MBl5Csfjb5ytauHrEP+uBFp9dQ2AgOBFSsqM0uqyzi77EOSJHWQJKm7JEnSE/GmAHYCMAewWwhxSw/pERERERHhZkoWhq0+jkW/X9IqMMd1bobQKV1weLsVXnutcCisLsbGwGuvAb//DjyxhDhRhWFP5iOeANYDuCNJUiyAJABOANoAMAEQDeD/9JceERERURUzMytc77K8x1KFEUJg66lbmLc7Goo8pSbuYGmCJa/4oYubbWHAGPjuO+Djj4G1a4Hw8MJlSiwsCif5GT+eQ2Sp8rHIfOQEgFUAnkPhEiOdASgARAHYCmBVaZYuISIiIqo1JAngMlJ6dz8zF+/tOI8/Yu5qxfv5OWJ+fx9YmRXttrS3B95/v/BFVNVYZD4khLgI4E1950FEREREpHb40l28s+087mXmamKWJkaYP8AH/f0b6zEzouKxyCQiIiIiqmYUuQX4ZN9F/O/EDa145xY2WPKKHxpZmeopM6KnY5FJRERERFSNnL7xAG9tjkLC/SxNTG5kgHdf9MSYjk1hYFDy0iRE+sYisxp6fIJblUoFAwNOAkxEz+7x5WGemEibiIiqgXylCisPXcFXYXFQPZo4Ft6Ollg+1B9u9hb6S46oDFhkVkOSJEEmkyEvLw8KhQJWVk9f64iI6GkUCgUAQCaTscgkIqpm4v7NxPTNUTifmKaJGUjAm0EtMOV5N8iM2OlANQeLzGrKwsIC9+/fx927hbOImZubs0eTiMpFpVJBoVBofp5YWPCTcCKi6kIIgU3Hr2PhvovILXg04qSJtRmWDvFD26bWesyOqHxYZFZTNjY2UCgUyMnJwe3bt/WdDhHVEiYmJrCxsdF3GkREBOBueg5mbj2LiCv3tOLD2jnjgz5eqCfnr+pUM/E7t5oyNDREkyZNcP/+fWRkZCAvL0/fKRFRDSaTyWBhYQEbGxsYGhrqOx0iojov9NwdvL/zPNKy8zUxG3MZPhvki2Avez1mRvTsWGRWY4aGhmjYsCEaNmwIIQSEEE8/iIjoCZIk8RlMIqJqIi07H3N3R2PnmUSteM+WDfHpy76ws5DrKTOiisMis4bgL4lERERENdtfV+9h5pazuJ2Wo4mZyQzxUR8vDG3nzN/1qNZgkUlEREREVIly8pX44o/LWHMsHo8PTAtoUh/LhvrDxcZcf8kRVQIWmURERERElSTmdjqmb47C5bsZmpiRgYTpwe6Y2K05jAy5egDVPiwyiYiIiIgqmFIl8H3ENXzxx2XkKx91X7ZoWA/Lh/rDpzHXQafai0UmEREREVEFupmShRlbzuLvhBSteEjnppjVyxMmxpzlm2o3FplERERERBVACIFtp25h3p4YZOYWaOIOliZY8oofurjZ6jE7oqrDIpOIiIiI6Bndz8zF+zvPY3/0Xa14Xz9HLOjvAyszYz1lRlT1WGQSERERET2Dw5fu4p1t53EvM1cTszQxwvwBPujv31iPmRHpB4tMIiIiIqJyyMorwCehF/HTiRta8U6uNljyih8c65vqKTMi/WKRSURERERURmduPMD0zVFIuJ+licmMDDCrlydCOjWFgYGkx+yI9ItFJhERERFRKeUrVVh5OA5fh8VBqXq0NIlXI0ssH+YPd3sLPWZHVD2wyCQiIiIiKoWryZl4a3MUzt5K08QMJOD1QFdM6+kOmZGBHrMjqj5YZBIRERERlUAIgR8jr+OTfReRk6/SxJ2tTbF0iD/aNbXWY3ZE1Q+LTCIiIiKiYtxNz8Hb287haGyyVnxoW2d82NcL9eT8dZroSfxXQURERESkw77zd/D+zvNIzcrXxGzMZfj05VZ4wdtBj5kRVW8sMomIiIiIHpOek4+5u6Kx40yiVvx5z4b4bJAv7CzkesqMqGZgkUlERERE9NDxq/cxc+tZJKZma2JmMkN81McLQ9s5Q5K4NAnR07DIJCIiIqI6LydfiS/+uIw1x+IhHq1MgoAm9bF0iD+a2prrLzmiGoZFJhERERHVaRfvpGP65ihcSsrQxIwMJEzr6YbXA11hZMilSYjKgkUmEREREdVJSpXAmohr+OKPWOQpHy1N4mpnjuVDW6OVk5UesyOquVhkEhEREVGdczMlCzO2nsXf8Sla8bGdmmJWL0+Yygz1lBlRzccik4iIiIjqDCEEdpxOxJzd0cjMLdDE7S3l+HywH7q52+kxO6LagUUmEREREdUJKYo8zN55Hr9dSNKK9/ZthE8G+KC+mUxPmRHVLiwyiYiIiKjWC7v8L97Zdg7JGbmamIWJERYM8EE/P0cuTUJUgVhkEhEREVGtlZVXgIX7LuLHyBta8U6uNljyih8c65vqKTOi2otFJhERERHVSmduPMBbW84i/p5CE5MZGWBWL0+EdGoKAwP2XhJVBhaZRERERFSr5CtV+OpwHL4Ki4NSJTRxr0aWWD7MH+72FnrMjqj2Y5FJRERERLXG1eRMvLU5CmdvpWlikgS8HuiK6T3dITMy0GN2RHUDi0wiIiIiqvGEEPgx8jo+2XcROfkqTdypgSmWDfVHu6bWesyOqG5hkUlERERENdq/6Tl4e9s5HIlN1ooPaeuED/t4wcLEWE+ZEdVNLDKJiIiIqMb67fwdvLfzPFKz8jUxa3MZFg5shV4+DnrMjKjuYpFJRERERDVOek4+5u6Oxo7TiVrxHp4N8dmgVmhoYaKnzIiIRSYRERER1SiR1+5jxpazSEzN1sRMjQ3xYR8vDG/vDEni0iRE+sQik4iIiIhqhNwCJZb+EYvVEdcgHq1MgtZN6mPZEH80tTXXX3JEpMEik4iIiIiqvYt30jF9cxQuJWVoYkYGEqY+74Y3glxhZMilSYiqCxaZRERERFRtKVUCa49dw5L9schTPlqapLmdOZYP9YevU309ZkdEurDIJCIiIqJq6daDLMzYchYn4lO04mM7NcWsXp4wlRnqKTMiKgmLTCIiIiKqVoQQ2HE6EXN3RyMjt0ATt7eU4/PBfujmbqfH7IjoaVhkEhEREVG1kaLIw+yd5/HbhSSteO9WjfDJQB/UN5PpKTMiKi0WmURERERULYRf/hdvbzuH5IxcTczCxAjz+/ugv78jlyYhqiFYZBIRERGRXmXnKbFw30X8EHldK96xuQ2WPxckzQAAIABJREFUDPFD4/qmesqMiMqDRSYRERER6U3UzVS8tTkK1+4pNDGZoQHe6eWBcZ2bwcCAvZdENQ2LTCIiIiKqcgVKFb4Ou4ovD1+BUiU08ZaNLLF8qD88HCz0mB0RPQsWmURERERUpa4lZ2L6lrM4ezNVE5MkYGI3V0wPdoPciEuTENVkLDKJiIiIqEoIIfDTiRv4JPQisvOVmrhTA1MsHeKP9s2s9ZgdEVUUFplEREREVOn+Tc/BO9vPIfxyslb8lTZO+KivFyxMjPWUGRFVNBaZRERERFSpfr9wB+/tOI8HWfmamLW5DAsHtkIvHwc9ZkZElYFFJhERERFVivScfMzbHYPtp29pxXt4NsRng1qhoYWJnjIjosrEIpOIiIiqnhBAVlb5jjUzK5wlhqq1E9fu460tZ5GYmq2JmRob4oM+LTGifRNIvIdEtRaLzBJIkmQN4BIAOwCXhRCeek6JiIiodsjKAurVK9+xmZmAuXnF5kOldvcusGYNcOQIkJEBWFgAQUHA+PGAvf3/s3ffUXaV9f7H38/0zGTSC4QQQnqfARSRGkCQjrREvVevXBS9V72QGFACSqSJlADegvoDxXavCU16E+lFEJhJJYFAIAnpdZJMn+f3xwwnOZBAIGdmT3m/1jprku85++zPgrVm5pNnn/1AdV090x9dyK+feYu4bWcSSvfuxg0TS9m3l//vpPbOkvnRpgO9kg4hSZKUtMpKOO88uO02qK1Nf+6xx2DaNDjznE2sHVbGgpUVqeeyswLnHT2Ufx8/mJzsrBbNLCkZlsydCCEcDfwL8Gvg3ITjSJIkJaayEo4/vnH1csciBaVv81zRAsLKhtR0UK8ibphYSsne3Vokp6TWwZK5AyGETsAvgXnAdVgyJUlSB3beeTsvmNnFlfQ6sYyCfdalzb/++X246PiRdMrLboGEkloTS+aOXQoMBsYDtR/9UkmSpPZrxYrGS2Q/LFI0ahk9jp1LVn5dalpXkc/Gx8bx3Ul96JTXUikltSZeGP8BIYRxwA+A38YYn046jyRJUpJuvfXDn8HMKqih16mv0evk8rSCueX1PVj+m8PZ/EYfbr21hYNKajVcydxOCCEL+H/ABuDC3XyvuTt5avDuvK8kSVJL+uBlsgX7rqbn8eXkFFenZg3VOax7dDRb5u0FNG5N8uSTMHVqy+WU1HpYMtN9HzgQODvGuDbpMJIkSUmraLpRbMipp9v4+XQ54J2056ve6cGaB0qpr+i0w+MkdTyWzCYhhL2BK4CnYoy37e77xRhH7+Q8c4FRu/v+kiRJLaG4GPL22ECvk8rI7bklNY91Wax/ejgVL+/L+6uXHzxOUsdkydzmf4A84N+SDiJJktQa1NU3UHzQIvYofYOQFVPzmlXFrLmvlNo1XXZ67PjxLRBQUqtkydzmJBo/i3lzCGn/GlfQ9HVACOHJ918bY9zcgtkkSZJa1NtrtjBpRhllNRsITbeKjBE2/X0QG54dBvU735okNxfOOaeFgkpqdSyZ6boBR+zkuU7bPed/N0mS1C7FGPnfl97livvnU1lbn5rXbezEmvtLqF7a82Pf4+yzoW/f5kwpqTWzLDWJMX74wwRACGEg8DawIMY4oiUzSZIktaRVFVX88I5ZPLFgddr8SyX9efWWUSxbmvux73HEEXDTTc2VUFJb4D6ZkiRJ4uE5K/jiDU+nFczuhbnc/E/7c+NXSnjk/lzOPbfxUtgdyc2Fc8+Fhx+GgoIdv0ZSx+BKpiRJanmFhbD5U97eoLAws1k6uIqqWi67bx63v7I0bT5+eG+uOWMcfbo0NsZOneBXv4LLLoNbb23cB7OiovEusuPHN34G00tkJYElU5IkJSEEKCpKOkWH99Lb65g8s4yl6ytTs0652Vx84kj+6XMD+MDNEIHGIjl1auNDknbEkvkxYoyL2dHmT5IkSW1UdV09Nzz2Br96ehFx284klOzdjRsmlDCod+fkwklq8yyZkiRJHciCFRWcP6OM+cs3pWbZWYH/OGoo3z1yMDnZ3rJD0u6xZEqSJHUADQ2R3zz3Ntc8vICa+obUfFCvIm6YWErJ3t0STCepPbFkSpIktXPLNlQyZWY5L7y1Nm3+tYP2YeoJI+mUl51QMkntkSVTkiSpnYox8peyZfzkL3OpqK5LzXsX53PNmeM4cnifBNNJaq8smZIkSe3Qhq01XHz3HB6YvTxtfvyYPbjytLH0KMpLKJmk9s6SKUmS1M48vXA1F9xRzspN1alZcX4OPz11NKftt9cOtyaRpEyxZEqSJLUTlTX1XP3QfH73wjtp8wP37cH0CSX0716YUDJJHYklU5IkqR2YtXQDk2aUsWj1ltQsLzuLKV8cxjmHDiI7y9VLSS3DkilJktSG1dU3cPOTi7jp8Teoa4ip+Yg9irlhYikj9+ySYDpJHZElU5IkqY1avGYLk2aW8dq7G1KzEODcwwYx+dhh5Oe4NYmklmfJlCRJamNijPzfS0u4/P55VNbWp+Z7devE9RNKOGhQzwTTSeroLJmSJEltyKqKKn5052z+9vqqtPkZ+/fn0lNG0aUgN6FkktTIkilJktRGPDJ3BRfdNZt1W2pSs26FufzstLEcP3bPBJNJ0jaWTEmSpFauoqqWy+6bx+2vLE2bHzGsN9eeOY4+XQoSSiZJH2bJlCRJasVeXryOSTPKWLq+MjUryM3i4hNH8c+fG0AIbk0iqXWxZEqSJLVCNXUN3PDXhfzyqUXEbTuTUNK/K9MnljK4d+fkwknSR7BkSpIktTILV1Zw/p/LmLd8U2qWnRX43pFD+N5RQ8jNzkownSR9NEumJElSK9HQEPnNc29zzSMLqKlrSM337VXE9Akl7Dege4LpJGnXWDIlSZJagfc2VDLl9nKeX7Q2bf7PBw1g6gkjKczz1zZJbYPfrSRJkhIUY+Te8ve45C9zqKiqS817F+dzzZnjOHJ4nwTTSdInZ8mUJElKyIatNVzylzncP2t52vy40Xtw1elj6VGUl1AySfr0LJmSJEkJeOaN1Uy5vZyVm6pTs875Ofz0lNGcvv9ebk0iqc2yZEqSJLWgypp6fv7w69z2/OK0+YEDe3D9hBL27lGYTDBJyhBLpiRJUguZvXQj5894jUWrt6RmudmBKccO55uHDSI7y9VLSW2fJVOSJKmZ1dU3cPOTi7jp8Teoa4ip+fC+xdwwsZRR/bokmE6SMsuSKUmS1IzeWbuFSTPKePXdDalZCPCtwwYx+ZhhFORmJ5hOkjLPkilJktQMYoz8+eUlXH7/PLbW1Kfme3XrxHVnlfD5wT0TTCdJzceSKUmSlGGrK6r50Z2zePz1VWnz0/fbi2mnjqZLQW5CySSp+VkyJUmSMujRuSu46K7ZrN1Sk5p1K8zlyi+N5cRxeyaYTJJahiVTkiQpAzZX13H5ffOY8Y8lafPDh/Xm2jPH0bdLQULJJKllWTIlSZJ20z8Wr2PSzDKWrKtMzQpys7j4hJH880H7EIJbk0jqOCyZkiRJn1JNXQM3/nUhv3xqEdvtTEJJ/65Mn1jK4N6dkwsnSQmxZEqSJH0KC1dWcP6fy5i3fFNqlp0V+O6RQ/j+UUPIzc5KMJ0kJceSKUmS9Ak0NER++/xifv7w69TUNaTmA3sWMn1iKfsP6J5gOklKniVTkiRpFy3fWMmU28t57s21afN/+twALj5xJIV5/molSX4nlCRJ2gX3lC3jx3+Zw6aqutSsV+d8rjlzLEeN6JtgMklqXSyZkiRJH2Hj1louuWcO95W/lzb/4ui+XHXaWHp2zk8omSS1TpZMSZKknXj2jTVMub2cFZuqUrPO+TlcevIozjygv1uTSNIOWDIlSZI+oKq2nqsfep3bnl+cNj9wYA+un1DC3j0KkwkmSW2AJVOSJGk7c5Zt5PwZZby5anNqlpsd+MGxw/nWYYPIznL1UpI+iiVTkiQJqKtv4FdPv8UNjy2kriGm5sP7FnPDxFJG9euSYDpJajssmZIkqcN7Z+0WJs8s55V31qdmIcA3D92XHxw7nILc7ATTSVLbYsmUJEkdVoyRGS8v4bL757G1pj4179e1gOsmlHDw4F4JppOktsmSKUmSOqTVFdVcdNcs/jp/Vdr8tP32Ytopo+naKTehZJLUtlkyJUlSh/PYvJX86M5ZrN1Sk5p17ZTLlaeN4aRx/RJMJkltnyVTkiR1GJur67ji/nn8+eUlafPDhvbi2jNL2KNrQULJJKn9sGRKkqQO4ZV31jFpRjnvrtuamuXnZHHxiSP52kH7EIJbk0hSJlgyJUlSu1ZT18BNjy/k5icXsd3OJIzr35XpE0oZ0qdzcuEkqR2yZEqSpHbrjZUVnD+jjLnvbUrNsgJ878ghfP/ooeRmZyWYTpLaJ0umJElqdxoaIr97YTFXP/Q61XUNqfnAnoVMn1jK/gO6JxdOkto5S6YkSWpXlm+s5ILbZ/Hsm2vS5l/93AAuPmEkRfn++iNJzcnvspIkqd24t/w9Lrl7Npuq6lKzXp3zuebMsRw1om+CySSp47BkSpKkNm/j1lp+fM8c7i1/L21+7Ki+/Oz0sfTsnJ9QMknqeCyZkiSpTXv2jTVMub2cFZuqUrOivGwuPWU0Zx3Q361JJKmFWTIlSVKbVFVbz88ffp3fPrc4bf7Zgd2ZPqGUvXsUJhNMkjo4S6Ykqf2JEbZu/XTHFhaCK1+t3pxlG5k0o4w3Vm1OzXKzA5OPGc65hw8iO8v/h5KUFEvmdkIIk4FDgbFAH6AAWAE8CVwTY5ybXDpJ0i7buhU6d/50x27eDEVFmc2jXbJyJdxyCzz1FFRUQHExjB8P55wDfZvu2VPfEPnlU4u48a8Lqa2PqWOH9e3MDRNLGd2vazLhJUkplsx0U4EiYBYwu2k2Gvg68OUQwpdijA8lFU6SpPaoshLOOw9uuw1qa9Ofe+wxmDYNzj4bply6lYvuKeMf76xPe803D92XKV8cTkFudotlliTtnCUz3anAKzHGqu2HIYR/A/4HuCWEMCDGWJ9IOkmS2pnKSjj++MbVy52prY3874tLePSGecTsbT+C9+xawPVnlXDwkF4tkFSStKssmduJMT63k/nNTZfSDgGGA/NaNJgkSe3Ueed9dMHMKqym5xdnUzhsJXG7+ZdK+/HTU8fQtVNus2eUJH0ylsxd9/4/ndYkmkKSpHZixYrGS2R3ptPglfQ8fhbZRdt+9DZU5XLV6WP458P7NX9ASdKnYsncBSGEr9O4grkQeCvhOJIktQu33vrhz2AChNw6uh81j+LSJWnzyrd7sfbBEt7tUwCHt1BISdIn1mwlM4SQS2Mx6w10BTYCq4EFMcYd/EhpPUIIF9B4w58iYGTTn98DvhpjbNjF99jZnWgHZySkJElt3I4uk83fax09Tywnt/u2LWgaarPY8OQIKl4dCASefBKmTm2plJKkTyqjJTOE0Bv4BnAicCCQv4OXVYUQXgIeAH4XY1ydyQwZ8kXg6O3+vgT4WozxlYTySJLU7lRUbPeXrAa6HfIGXQ56k5C1bVy9ogtr7i+lbm3xjo+TJLU6GSmZIYShwGXAaUBe03gN8AqwDthE42pmd2AEcETT44oQwl3AT2KMb2YiSybEGL8AEELoRuOemT8BngwhXBJjvHIX32P0juZNK5yjMpVVkqS2qripN+b0rKDXieXk77kx9VxsgE0vDmHDc0OhIWuHx0mSWqfdLpkhhP8EzgWygSeA/wWejDG+/RHHDAKOBL4KTADOCCH8Osb4/d3Nk0kxxg3AMyGEE4AXgMtDCI/GGF9OOJokSW3e4YdHXly7mG7jXycrd9unUWrXF7Lm/lJq3uu+w+PGj2+hgJKkTyUTK5nnADcD18QY39uVA2KMb9F4A51bQwh7ARcC3wRaVcl8X4yxNoQwAzgAOBmwZEqStBtWbKxiVo9yehyzJm1eUTaA9X8bSazd8a8oublwzjktkVCS9GllomQOijGu+LQHxxiXAeeFEH6WgSzN6f2fgr0TTSFJUht3X/l7XPKXOWys3HYfwPoteax9aByVi/p+5LFnnw19P/olkqSE7XbJ3J2C2Rzv04yOaPq6KNEUkiS1URu31vKTe+dwT1n6hU+d1vVl4R/H0lC5o/sFbnPEEXDTTc2ZUJKUCVkf/5KOIYRwWAhhYggh5wPz3BDC94GvAZXAjEQCSpLUhj335hqOu+nptIJZlJfNNWeM45UbD+CbX8snN3fHx+bmwrnnwsMPQ0FBCwWWJH1qmd7CJAfoCaz7uL0wQwg9gM4xxnczmWE3DAZ+C6wJIbwCrAV60Xh32T2BKuAbMcYlO38LSVKrUFgImzd/+mOVMVW19Vz7yAJufTb9foCf2ac70yeUMqBn43/vX/0KLrsMbr0VnnyycZuS4uLGm/ycc46XyEpSW5KpLUx6ATcCp9O4N2ZtCOEhGrcmmb2Tw66ncXUwo0V3NzwFXEXjZbHjaCyYNcBi4A7gF61pmxVJ0kcIAYqKkk7R4c1ZtpFJM8p4Y9W2wp+bHZh0zDC+ffhgsrNC2uv79oWpUxsfkqS2KxNbmBQBTwPDgfd/WuQBpwLHhxCmxBj/a2eH7+75M6Vpy5WLk84hSVJbV98Q+dXTi7jhsYXU1sfUfGifztwwsZQxe3VNMJ0kqbllYhVxMjACeA34LlAODALOB/4VuCmEMCDGeGEGziVJklqxJeu2MnlmGS8vXp82/9dD9uXC44ZTkJudUDJJUkvJRMk8A9gEnBBjXNk0mwt8K4RwN/BH4AchhO7AuTHGuJP3kSRJbVSMkdtfWcpP753Llpr61HzPrgVcd1YJhwzplWA6SVJLykTJHAI8vV3BTIkxPhhCOAR4mMZVza4hhK/GGOsycF5JktQKrN1czUV3zebReem/Cpxa2o/LThlD18Kd3DZWktQuZaJk1tO4krlDMcb5TUXzMRpXPe8JIZyRgfNKkqSEPT5/JT+8cxZrNtekZl0KcrjitLGcUtIvwWSSpKRkomQuBkZ91AtijEtDCIfSuKJ5XNPXtRk4tyRJSsCW6jqueGA+//dS+k5khw7pxbVnjWPPrp0SSiZJSlomSuY/gG+EEAbFGN/a2YtijGtDCEcC9wLjAT+bKUlSG/TKO+uZPLOMd9ZuTc3yc7K46PgRfP3zA8nKajU3j5ckJSArA+9xH41bkUz+uBfGGDfTuJL5F1rR9iWSJOnj1dY3cP2jCzjrl8+nFcwxe3Xh/u8fyjcO2deCKUnKyErmI8C3gJqPeyFAjLEmhHAm8D2gewbOL0mSmtmbqzYzaUYZs5dtTM2yAvz7+CH8x9FDycvJxL9bS5Lag90umTHGSuDWT3hMA/CL3T23JElqXjFGfv/CO1z14Hyq6xpS8wE9Cpk+oYTPDOyRYDpJUmuUiZVMSZLUDq3cVMWU28t55o01afMvf3ZvLjlpFJ3z/TVCkvRhu/3TIYQwIsb4emt5H0mStPsemLWcqXfPZmNlbWrWsyiPq88YxzGj+iaYTJLU2mXinyDnhBBmAD+LMc75pAeHEEqBH9G4h6a7NUuSlKCNlbVMu3cud7+2LG3+hZF9ufqMsfTqnJ9QMklSW5GJknk58APgyyGEcuBPwFNAeYyx9oMvDiHkA/sBRwJfpXGPzS3AZRnIIkmSPqXnF61hysxy3ttYlZoV5mVz6cmjmPCZvQnBO8dKkj5eJm7889MQws3AxcDXgWtp3AOzNoSwGFgPVABdgB7APk3nDcBG4CYaV0FX724WSZL0yVXV1nPdIwu45dm30+YH7NOd6RNK2KdnUULJJEltUUY+sR9jXAWcF0L4ETABOAk4BBi2g5evAJ4BHgBmxhirdvAaSZLUAua+t5FJM8pYuHJzapaTFZh0zDC+c8Rgst33UpL0CWX0tnBN25n8rulBCKE30AfoSuOq5SpXLCVJSl59Q+TXT7/F9McWUFsfU/MhfTpz48RSxuzVNcF0kqS2rFnvPd5UKC2VkiS1IkvWbWXyzDJeXrw+bX72IQP54XEjKMjNTiiZJKk9cIMrSZI6iBgjd7yylJ/eN4/N1XWp+R5dCrjurBIOHdorwXSSpPYi4yUzhHAa8FqMcXGm31uSJH06azdXM/Xu2Twyd2Xa/OSSflxx6hi6FrqLmCQpM5pjJfNOIIYQNgJlwGtNjzJgfoyxvhnOKUmSduJvr6/kwjtms2ZzdWrWpSCHy780hlNL90owmSSpPWqOknkNMBo4ARjf9Hj/jgLVIYQ5NJbOl4DHXfGUJKl5bK2p48oH5vOnv7+bNj9kSE+uPbOEft06JZRMktSeNUfJvA54HNgA3AYsAjrRWDxPBD4DHAB8EyCE8ARwYYzx1WbIIklSh/Tau+uZNKOMxWu3pmZ5OVn86LgRfOPggWS5NYkkqZk0R8mcDuwDlH5wlTKEkAtMAS4D7gNqgVOAF0IIZ8cY/7cZ8kiS1GHU1jfwn397k/9+4k3qG7ZtTTK6XxdunFjK0L7FCaaTJHUEzVEyjwee2tFlsDHGWuBnIYSlwP8D9gcuBR4Cbg0hvBpjfL0ZMkmS1O4tWr2ZSTPKmLV0Y2qWFeDfxg/mvKOHkZeTlWA6SVJH0Rw/bfKBj7wGJ8b4B2ABMLWpVH6t6bjzmyGPJEntWoyR37+wmBN/8Uxawdy7RydmfvvzXPDFERZMSVKLaY6VzFeB8SGE/jHGpR/xunnA0QAxxmebbgh0dDPkkSSp3Vq5qYoL7pjF0wtXp80nfmZvfnzyKDrnuyW2JKllNcdPnpto3Mbk9hDCKTHG1Tt53Sig23Z/X0jjHWklSdIueHD2cqbePZsNW2tTs55Fefzs9LEcO3qPBJNJkjqyjJfMGOPdIYQbabz0dX4I4RfA7THG+QAhhG40fg5zLFD+gSzuoSlJ0sfYVFXLtHvmctdry9LmR4/ow9VnjKN3cX5CySRJap6VTGKMk5suf70OmAZcGkKoBbawbfUyNj33vqHAqubII0lSe/HCorVMub2cZRsqU7PCvGx+ctIoJn52b0JwaxJJUrKa7YMaMcbfhBDuAL4NnAyUAN2BdcArwGUxxucAQgjjgJHA75orjyRJbVlVbT3XP7qAW559m7htZxL2H9CN6RNKGdirKLlwkiRtp1nvBhBj3ARc2/QghJATY6zbwUtn01gyq5ozjyRJbdH85ZuYNKOM11dUpGY5WYHzvzCU7xwxmJxs7xwrSWo9mqVkhhBygG/SuGdmT2A9UAY8DDz3wdfHGCONW5pIkqQm9Q2RW555i+sfXUhNfUNqPrh3ETdO3I+x/bsmmE6SpB3LeMkMIRQCjwEHkb5f5onA1BDCLODcGOPLmT63JEntxZJ1W/nB7eW89Pa6tPk3Dh7Ij44fQUFudkLJJEn6aM2xknkp8HkaVyx/DrwFFACjgTOBk4AnQwgTYowPNMP5JUlqs2KM3PnqMqbdO5fN1ds+YdK3Sz7XnVXCYUN7J5hOkqSP1xwl8yxgKXBsjLFyu/mrwB9CCIcCdwF/CiGUxhgXN0MGSZLanHVbaph612wenrsibX7iuD258ktj6FaYl1AySZJ2XXOUzH7AfR8omCkxxmdDCKcBzwCTgf9ohgySJLUpTyxYxYV3zGJ1RXVqVlyQwxVfGsMpJf3cmkSS1GY0R8lcDXT5qBfEGJ8LIfydxhsDSZLUYW2tqeOqB+fzxxffTZsfPLgn151VQr9unRJKJknSp9McJfN54NQQwl4xxmUf8brFNO6dKUlSh/Tau+uZPLOct9dsSc3ycrL44XEjOPvggWRluXopSWp7mqNk/gI4A7grhHBijHHNTl43Ctiyk+ckSWq3ausb+K+/vcl/PfEm9Q0xNR+1Zxdu/HIpw/oWJ5hOkqTdk/GS2XQp7DTgMmBOCOFm4E/A2zHG+hBCN+ASYAxwe6bPL0lSa7Zo9WYmzyijfOnG1CwE+M4Rg5n0hWHk5WQlmE6SpN3XHCuZxBivCCGsBa4FftL0qA8hbAa60rh/5gYay6YkSe1ejJE/vvgOVz44n6rahtS8f/dO3DCxlM8O7JFgOkmSMqdZSiZAjPHmEMI9wLeBE4CxQDdgI/AwcEmMcVFznV+SpNZi1aYqLrhjFk8tXJ02n/CZ/vz4pFEUF+QmlEySpMxrtpIJEGN8D7i06UEIIS/GWNOc55QkqTV5aPZyLrp7Nhu21qZmPYry+NnpY/ni6D0STCZJUvPY7ZIZQvgv4LWmx5yPKpEWTElSR7GpqpZp987lrlfTb7R+1Ig+XH3GWPoUFySUTJKk5pWJlcx/B96/NV5dCGEe20rna0BZjHFzBs4jSVKb8OJba/nBzHKWbahMzTrlZvPjk0bxlQP3JgS3JpEktV+ZKJknAPs1PfYHxtG4/+W/ND0fQwiL2FY6XwVe+4itTSRJapOq6+qZ/uhCfv3MW8RtO5Ow34Bu3DChlIG9ipILJ0lSC9ntkhljfJjGG/kAEELoApSyrXTuB4wAhgITaFr1DCG8B7waYzx1dzNIkpS0+cs3MWlGGa+vqEjNcrIC5x09lH8bP5icbLcmkSR1DM2xT+Ym4OmmBwAhhHwa7y77funcr+nvJ2X6/JIktaT6hsitz77FdY8spKZ+29Ykg3oXcePEUsb175ZgOkmSWl6z3l32fTHGauAfTQ8AQghZNK5wSpLUJi1dv5UfzCzn72+vS5t/4+CB/PC4EXTKy04omSRJycnE3WXfAO4B7geeiTHW78pxMcYGYN7unl+SpJYWY+SuV5cx7d65VFTXpeZ9u+Rz7ZklHD6sd4LpJElKVqZWMicDk4CNIYQHgfuAh2OMGzP0/pIktQrrttRw8d2zeWjOirT5iWP35MrTxtCtMC+hZJIktQ6ZuPHP0BDCcOBU4GRgIvAVGrczeRa4F7g/xrhod88lSVKSnlywigvumMXqiurUrLggh8tPHcOppf3cmkSSJDK0khljXABcA1wTQugJnAicAhwDHAlMDyEsoLFw3gc8H+P2N3eXJKn1qqyp56rWRF8SAAAgAElEQVQH5/OHF99Jm39+UE+um1DCXt06JZRMkqTWpznuLrsW+D3w+xBCLnAUjSucJwEXAhcA60IID9BYOB+JMW7OdA5JkjKhbMkGJs8o4601W1KzvOwsLjxuOP96yL5kZbl6KUnS9pr17rIxxlrgkabH90II42i8rPYk4GvA14GaEML3Y4z/rzmzSJL0SdTVN/DfTyziF397g/qGbRffjNyzCzdOLGX4HsUJppMkqfVqkS1M3hdjnAXMAi4PIexB4wrnyUDXlswhSdJHeWv1ZibNLKd8yYbULAT49uGDmXTMUPJz3JpEkqSdycQWJqOAI2jcvmTOrh4XY1wB/L+mhyS1DTHC1q2f7tjCwsamolYrxsif/v4uVz4wn8rabTty9e/eiekTSjlw3x4JppMkqW3IxErmV4GLaNzC5EMlM4SQBZwFjAa2Ao/HGF/OwHkzKoRQCBxL48rqZ4GBQDbwJnAnMN3Pjkpi61bo3PnTHbt5MxQVZTaPdtnKlXDLLfDUU1BRAcXFMH48nHMO9O0LqzZVceGds3hyweq04846oD8/OXkUxQW5yQSXJKmNyUTJPBKoAW774BMhhALgSRpL2/uuDCHcCXw9xliVgfNnylfZtqo6F3gY6AIcDPwU+EoI4YgY46qE8kmSPoXKSjjvPLjtNqitTX/uscdg2jQ44dzlLO49mw2V217QoyiPq04by3Fj9mjRvJIktXWZKJkDgFdijJt28Nx/AAc2/XkusAz4PHAGkAWcmYHzZ0oNcDNwQ4zxjfeHIYQ9gQeA/YAbaSyjkqQ2oLISjj++cfVyR0JeLV2+MI+yzkuhctv8yOG9+fmZ4+hTXNAyQSVJakcyUTJ7A0/v5LlvAhH4nxjj9wGa9tF8EDgthHBMjPGxDGTYbTHG39O49coH58tDCN8FngdODyHkxRhrWjygJOkTO++8nRfM/P5r6XVSOTldt7XL7JjNZaeP5KsHDiD4+VlJkj6VrAy8RzWNl5WmCSEMBYY0/fXa9+dN+2h+GwjANzJw/pZQ3vQ1H+iZZBBJ0q5ZsaLxEtkPya6n2xHz6fvVF9MKZvV73Xjvt4fxhYH7WDAlSdoNmSiZbwKjdjD/QtPX12OM727/RIyxjMbi9vkMnL8lDGr6WgusSzKIJGnX3Hrrhz+DmdtrE3t+/Tm6HvRW6ka/sSGw4ZlhrPjj56laXcStt7Z8VkmS2pNMlMyHgYEhhC99YH46jZfK7uxy2LeAtnI3hfOavj4cY6zelQNCCHN39AAGN19MSdL70i+TjRR/9i32/JfnyOtTkZrWri1ixR8OZuPzQyE2/kh88skWjSlJUruTic9k/ifw78BvQgh9gBeAw4CjaSyZd+/kuHqgLgPnb1YhhBOAc2hcxfxxwnEkSbuooqlLZhdX0uvEMgr2Sb8QZdMr+7DhyZHEuuwdHidJkj6d3S6ZMcYVIYSvAXfQeHfW7ZXHGHdyywWGAct39/zNKYQwEvgjjZ8fvSDGWP4xh6TEGEfv5D3nsuPLiyVJGdS5OFI0ahk9jp1LVv62f9Osq8hn7UPjqHq7zw6PKy5uqYSSJLVPmVjJJMZ4fwjhc8A0YDyNN8h5nsYb/HxICKE/MAb4aybO3xyaMj4MdAemxxhvSjiSJGkXrd9SQ92Bc+hVl/5vmVte34N1j4yloSpvp8eOH9/M4SRJaucyUjIBmlb5TtvFl59H4+dBW8X2JR8UQuhFY7YBwG+BKckmkiTtqqcWruaC28tZVbftI/QN1Tmse3Q0W+btRePFKTuWmwvnnNMCISVJascyVjI/oQbgOeCBhM6/UyGEYuAhYARwF/CtGGNMNpUk6eNU1tTzs4fm8/sX3kmbV73TgzUPlFJf0elj3+Pss6Fv3+ZKKElSx5BIyYwx/jCJ836cEEI+cA/wGeAR4CsxxvpkU0mSPk75kg1MmlHGW2u2pGZ52Vmcf9RwZly+LysrPn7fyyOOgJv8YIQkSbstE1uYtAshhGzg/4AjgWeA02OMNcmmkiR9lLr6Bm766xucfvPzaQVzxB7F3Pv9Q/j3owfx8EOBc89tvBR2R3Jz4dxz4eGHoaCghYJLktSOJXW5bGv0PbZ9pnQN8D8h7PBfvqfEGNe0WCpJrUthIWze/OmPVca8vWYLk2aUUbZkQ2oWApx7+CAmHzOM/JzGrUk6dYJf/QouuwxuvbVxH8yKisa7yI4f3/gZTC+RlSQpcyyZ23Tf7s8fdQOjaTSWUEkdUQhQVJR0ig4txsj/vvQuV9w/n8rabZ9o2KtbJ6ZPKOFzg3ru8Li+fWHq1MaHJElqPpbMJjHGaTQWSElSK7Wqooof3jGLJxasTpufsX9/pp0yiuKCnVwTK0mSWowlU5LUJjw8ZwUX3TWL9VtrU7PuhblcddpYjh+7Z4LJJEnS9iyZkqRWraKqlsvum8ftryxNm48f3ptrzhhHny7erUeSpNbEkilJarVeensdk2eWsXR9ZWrWKTebi08cyT99bgA7uUGbJElKkCVTktTqVNfVc8Njb/CrpxcR47Z5yd7duGFCCYN6d04unCRJ+kiWTElSq7JgRQXnzyhj/vJNqVl2VuD7Rw3he0cOISfbLZ4lSWrNLJmSpFahoSHym+fe5pqHF1BT35Ca79uriBsmllK6d7cE00mSpF1lyZQkJW7ZhkqmzCznhbfWps2/dtA+XHTCCArz/HElSVJb4U9tSVJiYozcU/YeP75nDhVVdal57+J8rjlzHEcO75NgOkmS9GlYMiVJidiwtYaL/zKHB2YtT5sfP2YPrjxtLD2K8hJKJkmSdoclU5LU4p5euJoL7ihn5abq1Kw4P4efnjqa0/bby61JJElqwyyZkqQWU1lTz9UPzed3L7yTNj9w3x5Mn1BC/+6FCSWTJEmZYsmUJLWIWUs3MGlGGYtWb0nN8rKzmPLFYZxz6CCys1y9lCSpPbBkSpKaVV19Azc/uYibHn+DuoaYmo/Yo5gbJpYycs8uCaaTJEmZZsmUJDWbxWu2MGlmGa+9uyE1CwHOPWwQk48dRn5OdoLpJElSc7BkSpIyLsbI/720hMvvn0dlbX1qvle3Tlw/oYSDBvVMMJ0kSWpOlkxJUkatqqjiR3fO5m+vr0qbn77/Xkw7ZTRdCnITSiZJklqCJVOSlDGPzF3BRXfNZt2WmtSsW2EuV502lhPG7plgMkmS1FIsmZKk3ba5uo7L7pvLzH8sTZsfMaw31545jj5dChJKJkmSWpolU5K0W15evI7JM8tYsq4yNSvIzeLiE0fxz58bQAhuTSJJUkdiyZQkfSo1dQ3c8NeF/PKpRcRtO5NQ0r8r0yeWMrh35+TCSZKkxFgyJUmf2MKVFZz/5zLmLd+UmmVnBb535BC+d9QQcrOzEkwnSZKSZMmUJO2yhobIb557m2seWUBNXUNqvm+vIqZPKGG/Ad0TTCdJkloDS6YkaZe8t6GSKbeX8/yitWnzfz5oAFNPGElhnj9SJEmSJVOS9DFijNxb/h6X/GUOFVV1qXnv4nyuOXMcRw7vk2A6SZLU2lgyJUk7tWFrDZf8ZQ73z1qeNj9u9B5cdfpYehTlJZRMkiS1VpZMSdIOPfPGaqbcXs7KTdWpWef8HKadMpoz9t/LrUkkSdIOWTIlSWmqauu5+qHXue35xWnzAwf24PoJJezdozCZYJIkqU2wZEqSUmYv3cj5M15j0eotqVludmDKscP55mGDyM5y9VKSJH00S6Ykibr6Bn751CJu/Osb1DXE1Hx432JumFjKqH5dEkwnSZLaEkumJHVw76zdwqQZZbz67obULAT41mGDmHzMMApysxNMJ0mS2hpLpiR1UDFG/vzyEi6/fx5ba+pT8726deK6s0r4/OCeCaaTJEltlSVTkjqg1RXV/OjOWTz++qq0+en77cW0U0fTpSA3oWSSJKmts2RKUgfz6NwVXHTXbNZuqUnNuhXmcuWXxnLiuD0TTCZJktoDS6YkdRCbq+u4/L55zPjHkrT54cN6c+2Z4+jbpSChZJIkqT2xZEpSB/CPxeuYNLOMJesqU7OC3CymnjCSrx20DyG4NYkkScoMS6YktWM1dQ3c+NeF/PKpRWy3Mwnj+ndl+oRShvTpnFw4SZLULlkyJamdWriygvP/XMa85ZtSs+yswHePHML3jxpCbnZWgukkSVJ7ZcmUpHamoSFy2/OLufrh16mpa0jNB/YsZPrEUvYf0D3BdJIkqb2zZEpSO7J8YyVTbi/nuTfXps3/6XMDuPjEkRTm+W1fkiQ1L3/bkKR24p6yZfz4L3PYVFWXmvXqnM81Z47lqBF9E0wmSZI6EkumJLVxG7fWcsk9c7iv/L20+RdH9+Wq08bSs3N+QskkSVJHZMmUpDbs2TfWMOX2clZsqkrNOufncOnJozjzgP5uTSJJklqcJVOS2qCq2nqufuh1bnt+cdr8wIE9uH5CCXv3KEwmmCRJ6vAsmZLUxsxZtpHzZ5Tx5qrNqVluduAHxw7nW4cNIjvL1UtJkpQcS6YktRF19Q386um3uOGxhdQ1xNR8WN/O3DhxP0b165JgOkmSpEaWTElqA95Zu4XJM8t55Z31afNvHrovU744nILc7ISSSZIkpbNkSlIrFmNkxstLuOz+eWytqU/N+3Ut4LqzSjh4SK8E00mSJH2YJVOSWqk1m6v50Z2z+ev8lWnz0/bbi2mnjKZrp9yEkkmSJO2cJVOSWqHH5q3kR3fOYu2WmtSsa6dcrjxtDCeN65dgMkmSpI9myZSkVmRzdR1X3D+PP7+8JG1+2NBeXHtmCXt0LUgomSRJ0q6xZEpSK/HKO+uYNKOcd9dtTc3yc7K4+MSRfO2gfQjBrUkkSVLrZ8mUpITV1DVw0+MLufnJRWy3Mwnj+ndl+oRShvTpnFw4SZKkT8iSKUkJemNlBefPKGPue5tSs6wA3ztyCN8/eii52VkJppMkSfrkLJmSlICGhsjvXljM1Q+9TnVdQ2o+sGch0yeWsv+A7smFkyRJ2g2WTElqYcs3VnLB7bN49s01afOvHDiAS04cSVG+35olSVLb5W8yktSC7i1/j0vuns2mqrrUrFfnPH5+xjiOHtk3wWSSJEmZYcncTgjhAOAY4EDgc0A/oDrG6J4BknZq5Uq45RZ46imoqIDiYhg/Hs45B/o29caNW2v58T1zuLf8vbRjjx3Vl5+dPpaenfNbPrgkSVIzsGSm+zFwatIhJLUNlZVw3nlw221QW5v+3GOPwbRpcPbZ8OXz1jD1nnKWb6xKPV+Ul82lp4zmrAP6uzWJJElqVyyZ6V4AyoGXmx4rko0jqbWqrITjj29cvdyZuljP7W8u4JHfv502/+zA7kyfUMrePQqbOaUkSVLLs2RuJ8b48+3/7uqCpJ0577yPLpi5fTbS6+Qy8npt3jbLDkw+ZjjnHj6I7Cy/v0iSpPbJkilJn9CKFY2XyO5QiHT53CK6HbqQkB1T49o1nfnDeaUcNrZri2SUJElKiiVTkj6hW2/98GcwAXK6bqXnSWUU9F+fNt/00r6sf3o4zwzK5rCxLRRSkiQpIZbMZhJCmLuTpwa3aBBJGffhy2QjncctoftR88jKr09N6zYVsPaBEqre7QXAk0/C1KktFlOSJCkRlkxJ+oQqKrb9Oauwmp5fnE3hsJVpr9k8tx/rHhtDrM7d4XGSJEntlSWzmcQYR+9o3rTCOaqF40jKoOLixq+dBq+k5/GzyC6qST1XX5nLukfHsPX1fjs9TpIkqT2zZErSJ3Tw4XW8kjWP4tIlafPKt3ux9sES6jcX7PC48eNbIJwkSVLCLJmS9Am88s56Hssuo7h0a2rWUJvFhidHUPHqQGDHW5Pk5sI557RMRkmSpCRZMiVpF9TWN/CLx9/gv594k4ZtO5NQvaILa+4vpW7tR18Le/bZ0LdvM4eUJElqBSyZkvQx3lxVwaQZ5cxetjE1ywpQvGQIs/5vKDRkfeTxRxwBN93U3CklSZJaB0umJO1EQ0Pk9y8s5mcPvU51XUNqvk/PQqZPKGVUn+6cXwS//e2O983MzW1cwbzpJijY8cc0JUmS2h1L5nZCCCcCP/7AOC+E8OJ2f788xvhAC8aSlIAVG6u44I5ynnljTdr8KwcO4JITR1KU3/jt81e/gssug1tvbdwHs6Ki8S6y48c3fgbTS2QlSVJHY8lM1xv43Adm4QOz3i0XR1IS7it/j0v+MoeNlduWJ3t1zuPnZ4zj6JEfbo19+8LUqY0PSZKkjs6SuZ0Y423AbQnHkJSQjVtr+cm9c7in7L20+TGj+nL16WPp2Tk/oWSSJElthyVTkoDn3lzDlNvLWb6xKjUrysvm0pNHc9Zn+hPCjrcmkSRJUjpLpqQOraq2nmsfWcCtz76dNv/MPt2ZPqGUAT0LE0omSZLUNlkyJXVYc9/byKQZZSxcuTk1y8kKTDpmGN85YjDZWa5eSpIkfVKWTEkdTn1D5NdPv8X0xxZQWx9T86F9OnPDxFLG7NU1wXSSJEltmyVTUoeyZN1WJs8s4+XF69Pm/3rIvlx43HAKcrMTSiZJktQ+WDIldQgxRm5/ZSk/vXcuW2rqU/M9uxZw3VklHDKkV4LpJEmS2g9LpqR2b+3mai66azaPzluZNj+1tB+XnTKGroW5CSWTJElqfyyZktq1x+ev5Id3zmLN5prUrEtBDlecNpZTSvolmEySJKl9smRKape2VNdxxQPz+b+X3k2bHzqkF9eeNY49u3ZKKJkkSVL7ZsmU1O688s56Js8s4521W1Oz/JwsfnT8CP7l8wPJcmsSSZKkZmPJlNRu1NY38IvH3+C/n3iThm07kzC6XxdunFjK0L7FyYWTJEnqICyZktqFN1dtZtKMMmYv25iaZQX49/FD+I+jh5KXk5VgOkmSpI7DkimpTYsx8vsX3uGqB+dTXdeQmg/oUcj0CSV8ZmCPBNNJkiR1PJZMSW3Wyk1VTLm9nGfeWJM2//Jn9+aSk0bROd9vcZIkSS3N38AktUkPzFrO1Ltns7GyNjXrWZTH1WeM45hRfRNMJkmS1LFZMiW1KRsra5l271zufm1Z2vwLI/ty9Rlj6dU5P6FkkiRJAkumpDbk+UVrmDKznPc2VqVmhXnZXHryKCZ8Zm9CcGsSSZKkpFkyJbV6VbX1XPfIAm559u20+QH7dGf6hBL26VmUUDJJkiR9kCVTUqs2972NTJpRxsKVm1OznKzApGOG8Z0jBpOd5eqlJElSa2LJlNQq1TdEfv30W0x/bAG19TE1H9KnMzdOLGXMXl0TTCdJkqSdsWRKanWWrNvKD2aW89LidWnzsw8ZyA+PG0FBbnZCySRJkvRxLJmSWo0YI3e8spSf3jePzdV1qfkeXQq47qwSDh3aK8F0kiRJ2hWWTEmtwtrN1Uy9ezaPzF2ZNj+5pB9XnDqGroW5CSWTJEnSJ2HJlJS4v72+kgvvmM2azdWpWZeCHC7/0hhOLd0rwWSSJEn6pCyZkhKztaaOKx+Yz5/+/m7a/JAhPbn2zBL6deuUUDJJkiR9WpZMSYl47d31TJpRxuK1W1OzvJwsfnTcCL5x8ECy3JpEkiSpTbJkSmpRtfUN/Off3uS/n3iT+oZtW5OM2rMLN325lKF9ixNMJ0mSpN1lyZTUYhat3sykGWXMWroxNcsK8J0jBnP+F4aRl5OVYDpJkiRlgiVTUrOLMfKHF9/hqgfnU1XbkJrv3aMT0yeU8tmBPRJMJ0mSpEyyZEpqVis3VXHBHbN4euHqtPnEz+zNj08eRed8vw1JkiS1J/52J6nZPDh7OVPvns2GrbWpWc+iPH52+liOHb1HgskkSZLUXCyZkjJuU1Ut0+6Zy12vLUubHz2iD1efMY7exfkJJZMkSVJzs2RKyqgXFq1lyu3lLNtQmZoV5mXzk5NGMfGzexOCW5NIkiS1Z5ZMSRlRVVvP9Y8u4JZn3yZu25mE/Qd0Y/qEUgb2KkounCRJklqMJVPSbpu/fBOTZpTx+oqK1CwnK3D+F4bynSMGk5Pt1iSSJEkdhSVT0qdW3xC55Zm3uP7RhdTUb9uaZHDvIm6cuB9j+3dNMJ0kSZKSYMmU9KksWbeVH9xezktvr0ubf+PggfzwuBF0ystOKJkkSZKSZMmU9InEGLnr1WVceu9cNlfXpeZ9u+Rz7ZklHD6sd4LpJEmSlDRLpqRdtm5LDRffPZuH5qxIm584bk+u/NIYuhXmJZRMkiRJrYUlU9IueWLBKi68YxarK6pTs+KCHK740hhOKenn1iSSJEkCLJmSPsbWmjquenA+f3zx3bT5wYN7ct1ZJfTr1imhZJIkSWqNLJmSduq1d9czeWY5b6/Zkprl5WTxw+NGcPbBA8nKcvVSkiRJ6SyZkj6ktr6B//rbm/zXE29S3xBT81F7duHGL5cyrG9xgukkSZLUmlkyJaVZtHozk2eUUb50Y2oWAnzniMFM+sIw8nKyEkwnSZKk1s6SKQlo3Jrkjy++w5UPzqeqtiE179+9EzdMLOWzA3skmE6SJElthSVTEqs2VXHBHbN4auHqtPmEz/TnxyeNorggN6FkkiRJamssmVIH99Ds5Uy9ezbrt9amZj2K8rjqtLEcN2aPBJNJkiSpLbJkSh3Upqpapt07l7teXZY2P2pEH64+Yyx9igsSSiZJkqS2zJIpdUAvvrWWH8wsZ9mGytSsU242Pz5pFF85cG9CcGsSSZIkfTqWTKkDqa6rZ/qjC/n1M28Rt+1Mwn4DunHDhFIG9ipKLpwkSZLaBUum1EHMX76JSTPKeH1FRWqWkxU47+ih/Nv4weRkuzWJJEmSdp8lU2rn6hsitz77Ftc9spCa+m1bkwzqXcSNE0sZ179bgukkSZLU3lgypXZs6fqt/GBmOX9/e13a/BsHD+SHx42gU152QskkSZLUXlkypXYoxshdry5j2r1zqaiuS837FOdz7VklHDGsd4LpJEmS1J5ZMqV2Zv2WGqbePZuH5qxIm584dk+u+NIYuhflJZRMkiRJHYElU2pHnlywigvvmMWqiurUrLggh8tPHcOppf3cmkSSJEnNzpIptQOVNfVc9eB8/vDiO2nzzw/qyXUTStirW6eEkkmSJKmjsWR+QAihALgI+AowAFgHPAz8JMa4NMls0o6ULdnA5BllvLVmS2qWl53FhccN518P2ZesLFcvJUmS1HIsmdtpKpiPAwcDy4F7gIHA2cBJIYTPxxgXJZdQ2qauvoH/fmIRv/jbG9Q3xNR85J5duHFiKcP3KE4wnSRJkjoqS2a6qTQWzBeAY2OMmwFCCJOB64HfAEckF09q9NbqzUyaWU75kg2pWQjw7cMHM+mYoeTnuDWJJEmSkmHJbBJCyAW+3/TX775fMAFijNNDCP8CHB5COCDG+EoiIdXhxRj509/f5coH5lNZW5+a9+/eiekTSjlw3x4JppMkSZIsmds7FOgGLIoxvraD5+8AxgEnA5ZMtbhVm6q48M5ZPLlgddr8zAP6c+nJoyguyE0omSRJkrSNJXObkqavr+7k+Vc/8LpPLcbIhg0bPv6FTbp27fqhrSc2btxIjHEnR6TLy8ujsLAwbVZTU8PWrVt3OUNxcTHZ2emXYFZUVFBfX7+TI9Ll5OTQuXPntFldXR2bN2/eyREfVlRURG5uepHasmULtbW1u3R8VlYWXbp0SZs1NDSwadOmXc7QqdP/b+/eo6us73yPf74JISHcb4pyFajchIC0glaL2NZWbUWUwhxrtV3artMzulTUGRlUrFprtZJpl0tdp6dSZzrTAaxg1VGPetR6m1ERgoBcRFCUiyB3yIWQ7/nj2dlhJ9mQwEN+e2e/X2vt9SS/J8/2s9zJw/7s5/Jrp8LCwpSx8vJyVVZWptmioS5dujQYO9Lvw0srNusXT6/QrvIDUn6B8goK1bW4QL+6dKS+e9pJqqysbNbvVKdOnZSXl5cytnv3btXU1DRp+4KCArVv3z5l7MCBA9q3b1+aLRrq0KGD2rRJ3QXt3btX1dXVTdo+Pz9fHTumXnd68OBB7dmzp8kZiouL1bZt6ryh+/fvV1VVVZO2NzN17tw5ZczdtWvXriZnKCoqUlFRUcpYRUWFKioqmvwc7CMiubyPOFRhYaHatUu9o3RlZaXKy8ub/BzsIyLsIyLsIyLsI+qwj4i05n1EU1/LtNydR7STnS3JJc1Os74ksX5RE59veZpHRc+ePT3xXE167Nixw+vr3Llzk7e/6qqrGmw/Z86cZmVYvHhxg+coKSlp8vYTJkxosP0rr7zSrAwLFixo8ByTJk1q8vb9+/dvsP26deualaG0tLTBc1x//fXNeo7GNGf7jmMv9h8/9t++ZXd5cvvS0tJmPce6desaZOjfv3+Tt580aVKD7RcsWNCsDK+88kqD55gwYUKTty8pKWmw/eLFi5uVYc6cOQ2e46qrrmry9p07d26w/Y4dO5qVYdasWQ2eY9asWc16DvYR0YN9RPS4/vrrG2zPPiLCPqIO+4gI+4gI+4gI+4g6eXl5Lmm5H2W34khmndqPx9J9LLev3s8BQZw1uLse+/HXGnziBAAAAGQCSmad2nfsfoT1TeLuIxp9ErPlkoY357mAQw3t1YmCCQAAgIxl3sTz8Vs7M5st6UZJpe4+vZH1JZKWSHrf3ccew39n+bBhw4a/9dZbTd6GaykiuXAtxeotuzXjyQ+0Zkvd/5e8PNP/nDBI15x9itrk53EtRQLXUtRhHxHJhX1EU7CPiLCPqMM+IsI+IsI+IsI+ok5j+4hhw4Zp5cqVK9IdODsSSmaCmd0gqVTSfHef2sj6iyQ9I2mhu08+hv/O8uHDhw9fvnz50YdFq1NT4/rDG+v0wAurVHWwbgc9sEd7lU4brZK+Df9BAQAAAI6HESNGaMWKFUddMjldtk5ZYnl6mvW140tbIAtyyOc7y3XTvCX6r4+3p4xfeWZ/zbhgmNq1zU+zJQAAAJB5KJl13pS0S9IgMxvjDefKnJJYPtOysdBaua0NDGwAABeWSURBVLsWLvlcdyxcrj2Vdad2nNCxUPdPGaVzh5wQMB0AAABwdCiZCe5eZWYPSZop6SEzO9/d90mSmU2XNErSG+7+bsicaB127KvSbQuX6dkPNqWMXziyl355yUh1bd82zZYAAABAZqNkprpH0rcknSVpjZm9Lqm/pHGSvpT0k4DZ0Eq8tnqrbplfpi/21F3k37GwjX4xaYQmj+nNnWMBAACQ1SiZh3D3CjObKGmGpMslXSJph6THJd3u7htC5kN2K686qPue+1CPv/1Jyvi4U7rpwakl6tO1OM2WAAAAQPagZNbj7uWS7kg8gFgs/Wynbpi7RB9vrbs1d9v8PN3ynSG6+uxTlJfH0UsAAAC0DpRM4DiqPlijh19dq9+9vEbVNXXTBQ3t1VH//HejNbRXp8NsDQAAAGQfSiZwnKzbtk83zl2iJRvqJkg2k372jYGa/u1TVdiGqUkAAADQ+lAygZi5u/79nU91zzMfqvzAweR47y7tNHtqicYN7B4wHQAAAHB8UTKBGH2xp0L/+MRSvbJqa8r4Zaf30Z0XD1fHooJAyQAAAICWQckEYvL8ss2a8eRS7dh/IDnWtbhA904eqQtGnhQwGQAAANByKJnAMdpTcUB3Pb1C8xd9ljJ+7pCeuv+yUTqhU1GgZAAAAEDLo2QCx+Cddds1fd4SfbajPDnWriBfMy8aph+O6yczpiYBAABAbqFkAkehqrpGpS+t1qOvrZXXzUyikr5dVDq1RAN7dggXDgAAAAiIkgk006rNe3TD3CX6cNPu5Fh+num68wbr2omD1SY/L2A6AAAAICxKJtBENTWux95cp/tfWKWq6prk+Ck92qt02miN7tslYDoAAAAgM1AygSb4fGe5bp5Xprc//jJl/Efj+2vGhUNV3JY/JQAAAECiZAKH5e56aslG3f7UMu2pqE6O9+xYqPunjNLEIScETAcAAABkHkomkMbO/VWauXCZnl26KWX8gtN66ZeTR6pb+7aBkgEAAACZi5IJNOJvq7fqlifKtGV3ZXKsY2Eb/WLSCE0e05upSQAAAIA0KJnAIcqrDuq+5z7U429/kjJ+xindNHtqifp0LQ6UDAAAAMgOlEwg4YPPdumGuYu1duu+5Fjb/Dzd/J1TdfXZA5Wfx9FLAAAA4Egomch51Qdr9Mira/Xbl9eousaT40N7dVTptNEadlKngOkAAACA7ELJRE5bv22fps9bovc/3ZkcM5N+ds5ATT//VBW2yQ+YDgAAAMg+lEzkJHfXf7y7QXc/s0L7qw4mx3t3aacHp5Zo/MDuAdMBAAAA2YuSiZyzdU+lbv3LUr288ouU8UtP7607Lx6hTkUFgZIBAAAA2Y+SiZzywvLNmvHkB9q+ryo51qW4QPdOHqkLR54UMBkAAADQOlAykRP2VlbrrqeXa957n6WMTzi1px6YMkondCoKlAwAAABoXSiZaPXeXb9d0+ct0Ybt5cmxooI8zbxouK4Y109mTE0CAAAAxIWSiVarqrpGpS+t1qOvrZXXzUyikj6dNXvaaA3q2SFcOAAAAKCVomSiVVq9ZY9u+I8lWrFpd3IsP8907cTBuva8wSrIzwuYDgAAAGi9KJloVWpqXHPeWq9fP79SVdU1yfFTerTX7KklGtOva8B0AAAAQOtHyUSrsXFnuW6eX6a31n6ZMn7F+H76pwuHqbgtv+4AAADA8ca7brQKTy35XLctXKY9FdXJsZ4dC3X/lFGaOOSEgMkAAACA3ELJRFbbub9Ktz+1XE+XbUwZ/+6IXrr30pHq1r5toGQAAABAbqJkImu9vmarbpm/VJt3VyTHOhS20Z0Xj9Blp/dmahIAAAAgAEomsk7FgYO677mV+uNb61PGzxjQTQ9OLVHfbsVhggEAAACgZCK7fPDZLt0wd7HWbt2XHCvIN918/hBdc85A5edx9BIAAAAIiZKJrFB9sEaPvrZW//zSGlXXeHJ8yIkdVTpttIaf3ClgOgAAAAC1KJnIeJ98uU/T55Vp0Sc7kmNm0k/PGajp3z5VRQX5AdMBAAAAOBQlExnL3TX33Q2665kV2l91MDneu0s7/eYHJTpzUPeA6QAAAAA0hpKJjLR1T6VmPLlUL334Rcr4pWN6685JI9SpqCBQMgAAAACHQ8lExnlxxRbd+pel+nJfVXKsS3GBfnnJSF006qSAyQAAAAAcCSUTGWNvZbXufnqF5r63IWX8G6f21ANTRunETkWBkgEAAABoKkomMsJ767dr+rwyfbp9f3KsqCBPMy8cpivG95cZU5MAAAAA2YCSiaCqqmv025dX65FX1+qQmUk0qk9nlU4brUE9O4QLBwAAAKDZKJkIZs2WPbph7hIt37g7OZafZ/r7iYN13XmDVZCfFzAdAAAAgKNByUSLq6lx/fGt9brv+ZWqqq5Jjg/oXqzZ00br9H5dA6YDAAAAcCwomWhRm3aV65b5S/XGR9tSxn84rp9mXjRMxW35lQQAAACyGe/o0WL+WrZRty34QLsrqpNjPToU6v4pI3Xe0BMDJgMAAAAQF0omjrtd+w/o9qeW6a9lG1PGvzPiRN07eaS6dygMlAwAAABA3CiZOK7eWLNNN88v0+bdFcmxDoVtNOv7wzVlbB+mJgEAAABaGUomjouKAwf16+dXas6b61PGzxjQTQ9OLVHfbsVhggEAAAA4riiZiN2yz3fphrlL9NEXe5NjBfmmm84fop+eM1D5eRy9BAAAAForSiZic7DG9ehra1X64mpV13hyfMiJHVU6bbSGn9wpYDoAAAAALYGSiVh8+uV+3ThviRZ9siM5ZiZdc/Ypuun8ISoqyA+YDgAAAEBLoWTimLi75r23QXc9vUL7qg4mx0/uXKTfTC3RWYN6BEwHAAAAoKVRMnHUtu2t1IwnP9CLK7akjE8e01t3XjxCndsVBEoGAAAAIBRKJo7KSyu26NYnl2rb3qrkWOd2Bfrl5NP0vVEnB0wGAAAAICRKJpplX2W17nl2hf78zoaU8XO+0kMPTClRr85FgZIBAAAAyASUTDTZok+268a5Zfp0+/7kWGGbPM28aJh+NL6/zJiaBAAAAMh1lEwcUVV1jX738ho9/OpHOmRmEo3q01mzp47W4BM6hAsHAAAAIKNQMnFYH32xRzfMXaJln+9OjuWZdO3Ewbrum19RQX5ewHQAAAAAMg0lM8HM2ku6VNIZksZJKpHUVtIMd78vZLYQampcj7+9Xvc9t1KV1TXJ8QHdizV72mid3q9ruHAAAAAAMhYls85XJP1L6BCZYPOuCt3yRJleX7MtZfzycf0088Jhal/Irw0AAACAxtEW6uyR9AdJ70h6V9JlkmYGTRTA02UbddvCZdpVfiA51qNDoe6fMlLnDT0xYDIAAAAA2YCSmeDuayVdU/u9mU0KGKfF7dp/QHf8dZmeWrIxZfz84SfqV5eOVPcOhYGSAQAAAMgmlEzozY+26eb5Zdq0qyI51r5tvmZdPEI/GNuHqUkAAAAANBklM4dVHDioB15YpT+8sS5l/GsDumr21NHq2604UDIAAAAA2YqSeZyY2fI0qwa1aJDDeP+THSkFsyDfNP3bQ/SzbwxUfh5HLwEAAAA0H5Mc5rCzBvfQFeP7SZJOPbGDFv791/XzcwdRMAEAAAActVZzJNPMnpB0WjM3u9Ld3zkeedx9RGPjiSOcw4/Hf/No/NOFw9SrU5GuOWegigryQ8cBAAAAkOVaTcmUNEDSkGZuk/MXHRa3baNrz/tK6BgAAAAAWolWUzLd/auhMwAAAABAruOaTAAAAABAbCiZAAAAAIDYUDIBAAAAALFpNddkxsHMFkg6KfFtn8Tyf5nZJYmvN7n75JZPBgAAAADZgZKZaoyk/vXG+iYekvRJy8YBAAAAgOxCyTyEuw8InQEAAAAAshnXZAIAAAAAYkPJBAAAAADEhpIJAAAAAIgNJRMAAAAAEBtKJgAAAAAgNpRMAAAAAEBsKJkAAAAAgNhQMgEAAAAAsaFkAgAAAABiQ8kEAAAAAMSGkgkAAAAAiA0lEwAAAAAQG0omAAAAACA2lEwAAAAAQGzM3UNnyClmtruwsLDjoEGDQkcBAAAAgAbWrl2rysrKPe7e6Wi2p2S2MDPbLKlY0obQWQ5R23jXBk2Bw+E1yg68TtmB1ynz8RplB16n7MDrlPky8TXqK2m/u/c6mo0pmZCZLZckdx8ROgsax2uUHXidsgOvU+bjNcoOvE7Zgdcp87XG14hrMgEAAAAAsaFkAgAAAABiQ8kEAAAAAMSGkgkAAAAAiA0lEwAAAAAQG+4uCwAAAACIDUcyAQAAAACxoWQCAAAAAGJDyQQAAAAAxIaSCQAAAACIDSUTAAAAABAbSiYAAAAAIDaUTAAAAABAbCiZAAAAAIDYUDJzlJkVmdkvzGy1mVWY2UYze8zM+oTOhoiZjTWzW83sSTP73MzczCpC50IdMys2s0vM7A9mttTMdpvZPjMrM7M7zKxD6IyImNn0xN/SGjPbZWaVZvaJmT1uZiNC50NDZtbNzL5I7PtWhs6DiJm9mnhN0j2+Gzoj6phZLzMrTbzfKzez7Wa2yMzuD50t15nZuUf4W6p93BE669Ewdw+dAS3MzIokvSzpLEmbJL0uaYCkMyRtlXSmu68NFhCSJDNbKGlSveFKdy8KkQcNmdk1kn6f+Ha5pBWSOin62+ooaaWkCe7+RZiEqGVm2yS1l7RU0ueJ4RGSTpVUJekSd38uUDw0wsz+KOlKSSZplbsPDZsIUlQyJU2Q9BdJexv5kQfd/YMWDYVGmdmZkv5TUhdF/z4tU/Rv03BJfdy9TcB4Oc/Mhkq6Nc3qfElXJL4+z91faZlU8aFk5iAzu0vS7ZLelnS+u+9NjE+X9KCkv7n7hIARIcnM/lFSsaR3E4/NomRmFDO7UtJ4SaXuvuaQ8ZMkPStpjKQ/u/vlgSIiwcy+LmmRu1fUG/+5pIclbZTUz90PhsiHVGb2TUkvSfrfkn4mSmbGOKRknuLu68OmQTpmdrKiDz8LJf3Q3RfUW3+Gu78TJByOyMwuUPQBwQZJA9y9JnCkZqNk5hgzK5D0haJPtU5398X11pdJGiXpq+6+KEBEpGFmLkpm1kh8gvyWpEpJndy9KnAkpGFmayQNljTC3VeEzpPrzKydoiPOVZIukbRalMyMQcnMDmb2L5J+JOk6d38odB40j5n9m6TLJd3n7jNC5zkaXJOZe85WVDDX1i+YCU8klt9vuUhAq1SWWBZK6h4yCI6o9uglHwRkhlmSBkn6uaQDgbMAWcfMukqaKmmXpP8TOA6ayczaq+5yqT+FzHIsOBc795Qklu+nWf9+vZ8DcHQGJpYHJG0PGQTpJU55HqLoaNnHgePkPDMbJekmSXPc/W9mNiBsIhzG1WbWXVKNor+fhe7+aeBMiHxd0QecL0k6YGZTFB1kKFB0r4B57r4lYD4c3qWK7iGw2N2Xhw5ztCiZuadfYvlZmvWf1fs5AEfn+sTyeXevDJoESWZ2i6Ib/rSXNCzx9UZJl2fjNS+tiZnlKbqR1k5J/xA4Do7stnrf/8bM7nb3u4OkwaFq75i9RdHNHc+st/5XZvYTd5/fsrHQRLU3/PnXoCmOEafL5p7aKRX2p1m/r97PAWgmM7tQ0tWKjmLeHjgOUn1H0lWSpih6I7ZBUcHkGvTwrlN0l/Nb3P3L0GGQ1t8UXes3SNHN6YZImimpWtJdZnb9YbZFy+iaWF6p6D4bV0vqKekUSbMVfcj2p8SZA8ggZtZL0jcVXcbx58BxjgklM/dYYpnujk+WZhxAE5jZMEXXUJiiN8tlR9gELcjdv+XupuhN2DckrZL0qpnNDJsst5lZX0n3SHrN3f8YOA4Ow93vcPc/ufvH7l7u7qvd/V5FN2mSpF8kbt6EcPITyzaSprv7Y+6+zd3Xu/tNiu6/0VacMZCJLlf0+r3o7ptDhzkWlMzcsyexbJ9mfXFi2djcVwAOw8z6SHpeUYGZ7e6/DRwJabj7Tnd/XdKFkhZJutvMvhY4Vi57WNGb3p+HDoKj4+7/V9J7kjormtoJ4dS+16uR9Hgj6x9LLM9tkTRojlZxqqzENZm5qPai/D5p1vep93MAmsDMekh6UdH1zHMk3Rw2EZrC3Q+Y2VxJYxXdVfvdwJFy1fcUXYv5iFnKCTW1Uzb1S0ydIUnfq53fGRlnjaSvSjopdJActz6x3JzmngC1609okTRoksSZUGMUHehZGDjOMaNk5p7aU/dOT7O+dnxpC2QBWgUz6yjpOUlDJT0p6afOJMTZZFti2TNoCnRRNP9iY9odso73Lpmr9lpAPgQIq3aKuq5mZo38e1Q7rRavU2b5UWL5pLunu3dK1uB02dzzpqJ5kwaZ2ZhG1k9JLJ9puUhA9jKzQklPKfr0/gVJ/8PdDx5+K2SY2vKyNmiKHObu1thD0Y1KJGnVIeM7Q2ZF48ysp6RzEt+mmyYNLcDdP5C0TtGHM+Ma+ZFzE0tepwxh0Skclye+zfpTZSVKZs5x9ypJDyW+fSgx4askycymK7oL2RvuziljwBGYWb6iu79NVHSb+EsTf2PIIGZ2jplNM7M29cYLzOw6RZ8el0uaGyQgkCXMbLyZTbR65zQn5jNdoOh+D39193TTpKHl/Dqx/F3icg5JkpmNVTQXrSQ92uKpkM45kvormlLr/wXOEgtOOclN90j6lqSzJK0xs9cV/WKPk/SlpJ8EzIYEM7tIDae/aGtm/3XI93e7+7MtGAuprpU0OfH1NkkP13vvVetmd9/W2Aq0iEGKrpPdZmaLFO3nekgaqejasQpJP3b3DeEiAllhqKK/pU1mtlrSZkX3chir6PrZ5ZJ+Gi4eDvF7RVNh/EDSKjN7S9H0dGcpusnW7939iYD5kKr2hj//1lrmbKZk5iB3rzCziZJmKDo0f4mkHYruQHY7b7QyRk81PM3F6o1xDVlYXQ/5enLan5LuVN11f2h5r0m6V9FpsaMUFcwqRTe/eELS79z9o2DpgOzx35IeUfTv0HBJX1c0v/YSSfMlPeLu5eHioZa715jZ30l6VdI1ks5TNH3de5IedfdWcUpma5C47Kb2crU/hcwSJ+PeFAAAAACAuHBNJgAAAAAgNpRMAAAAAEBsKJkAAAAAgNhQMgEAAAAAsaFkAgAAAABiQ8kEAAAAAMSGkgkAAAAAiA0lEwAAAAAQG0omAAAAACA2lEwAAAAAQGwomQAAAACA2FAyAQAAAACxoWQCAAAAAGJDyQQAIAuY2aNm5mZ2a+gsAAAcDiUTAIDsUJJYlgVNAQDAEZi7h84AAAAOw8xM0h5J7SX1cffPA0cCACAtjmQCAJD5BisqmF9SMAEAmY6SCQBA5qs9VXZp0BQAADQBJRMAgMzH9ZgAgKxByQQAIPNRMgEAWYOSCQBA5uN0WQBA1uDusgAAZDAz6yJph6RqSR3cvTJwJAAADosjmQAAZLbao5irKZgAgGxAyQQAILNxPSYAIKtQMgEAyGxNLplm9qGZPVBvrK2ZrTazWcclHQAA9VAyAQDIbM256c+bksbXG7tZUhtJv44zFAAA6VAyAQDIUGaWL2lE4tumnC77pqSxZtYmsX1fSTMl3ejuFccnJQAAqbi7LAAAGcrMhktanvj2tcP86D+4+ztmdqqkVZLGuvv7ZvaEpPbufsHxzgoAQK02oQMAAIC0Sg75esJhfm6TJLn7ajPbJmm8mfWQ9H1Jpx3HfAAANMCRTAAAWhEze0rSfkljJC1w9xmBIwEAcgzXZAIA0Lq8KWmapPaS7gmcBQCQgyiZAAC0Lu9LMkk3ufu+0GEAALmH02UBAGhFzOxfJfV29/NCZwEA5CZu/AMAQJZLTHXSXdJlkiYruh4TAIAgKJkAAGS/CZJekvSxpB+4+5rAeQAAOYzTZQEAAAAAseHGPwAAAACA2FAyAQAAAACxoWQCAAAAAGJDyQQAAAAAxIaSCQAAAACIDSUTAAAAABAbSiYAAAAAIDaUTAAAAABAbCiZAAAAAIDYUDIBAAAAALGhZAIAAAAAYkPJBAAAAADEhpIJAAAAAIgNJRMAAAAAEBtKJgAAAAAgNpRMAAAAAEBs/j/mzsNoBFDtmwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1050x750 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "Lys = np.array(list(results_Ly.keys()))\n",
    "S0s = np.array([res['psi'].entanglement_entropy()[0] for res in results_Ly.values()])\n",
    "S1s = np.array([res['psi'].entanglement_entropy()[2] for res in results_Ly.values()])\n",
    "\n",
    "plt.figure(figsize=(7, 5))\n",
    "ax = plt.gca()\n",
    "ax.plot(Lys, S0s / np.log(2), 'bo', label='data at bond 0')\n",
    "ax.plot(Lys, S1s / np.log(2), 'rs', label='data at wrong bond')\n",
    "Lys_all = np.arange(0., max(Lys) + 2.)\n",
    "ax.plot(Lys_all, Lys_all - 1., '-', label='expected: $L_y - 1$')\n",
    "ax.axhline(0., linestyle='--', color='k')\n",
    "ax.set_xlabel('$L_y$')\n",
    "ax.set_ylabel('$S / log(2)$')\n",
    "\n",
    "ax.legend()\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As expected, we find the scaling $S(L_y) = \\log(2) L_y - \\log(2)$ at bond 0, indicating that $\\gamma = \\log(2)$. \n",
    "Cutting the MPS at bond 1 (red points) is a \"weird\" cut of the cylinder and does not show the expected scaling."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}