{ "cells": [ { "cell_type": "markdown", "id": "faa097af", "metadata": {}, "source": [ "# Samples and Measurement\n", "\n", "In this tutorial, we will introduce how to use Quantax to generate samples.\n", "\n", "In the beginning, we prepare ground states of transverse-field Ising models from which we will generate samples." ] }, { "cell_type": "code", "execution_count": 1, "id": "d2b6b218", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import jax\n", "import jax.numpy as jnp\n", "import quantax as qtx\n", "import matplotlib.pyplot as plt\n", "\n", "lattice = qtx.sites.Square(4)\n", "\n", "H0 = qtx.operator.Ising(h=1.0)\n", "E0, wf0 = H0.diagonalize()\n", "state0 = qtx.state.DenseState(wf0)\n", "\n", "H1 = qtx.operator.Ising(h=3.0)\n", "E1, wf1 = H1.diagonalize()\n", "state1 = qtx.state.DenseState(wf1)\n", "\n", "H2 = qtx.operator.Ising(h=10.0)\n", "E2, wf2 = H2.diagonalize()\n", "state2 = qtx.state.DenseState(wf2)" ] }, { "cell_type": "markdown", "id": "5cd5d0ea", "metadata": {}, "source": [ "## Exact sampler\n", "\n", "{py:class}`~quantax.sampler.ExactSampler` generates samples by sampling the full wavefunction. If the provided state is not a {py:class}`~quantax.state.DenseState`, it will generate a {py:class}`~quantax.state.DenseState` internally through forward pass. The sampler is defined below." ] }, { "cell_type": "code", "execution_count": 2, "id": "e50fdfbd", "metadata": {}, "outputs": [], "source": [ "sampler0 = qtx.sampler.ExactSampler(state0, nsamples=1024)\n", "sampler1 = qtx.sampler.ExactSampler(state1, nsamples=1024)\n", "sampler2 = qtx.sampler.ExactSampler(state2, nsamples=1024)" ] }, { "cell_type": "markdown", "id": "2bc17233", "metadata": {}, "source": [ "Now we can generate samples." ] }, { "cell_type": "code", "execution_count": 3, "id": "a5da46ec", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Samples(spins=Array([[ 1, 1, 1, ..., 1, 1, 1],\n", " [ 1, 1, 1, ..., 1, 1, 1],\n", " [-1, -1, -1, ..., -1, -1, -1],\n", " ...,\n", " [ 1, 1, 1, ..., 1, 1, 1],\n", " [-1, -1, -1, ..., -1, -1, -1],\n", " [ 1, 1, 1, ..., 1, 1, 1]], dtype=int8), psi=array([-0.62183999, -0.00981954, -0.62184007, ..., -0.62183999,\n", " -0.62184007, -0.62183999]), state_internal=None, reweight_factor=Array([1., 1., 1., ..., 1., 1., 1.], dtype=float64))\n", "spins shape (1024, 16)\n", "psi shape (1024,)\n" ] } ], "source": [ "samples0 = sampler0.sweep()\n", "samples1 = sampler1.sweep()\n", "samples2 = sampler2.sweep()\n", "\n", "print(samples0)\n", "print(\"spins shape\", samples0.spins.shape)\n", "print(\"psi shape\", samples0.psi.shape)" ] }, { "cell_type": "markdown", "id": "d0edd7c5", "metadata": {}, "source": [ "The generated sample is a PyTree with 4 terms: `spins`, `psi`, `state_internal`, and `reweight_factor`. Beginners usually only need `spins` and `psi`, which correspond to the generated spin configurations $s$ and their wavefunctions $\\psi(s)$. The probability of samples is given by $|\\psi(s)|^2$\n", "\n", "Let's check how magnetization changes by looking at\n", "\n", "$$\n", "\\left< M^2 \\right> = \\left< \\left( \\frac{1}{N} \\sum_i \\sigma_i^z \\right)^2 \\right>.\n", "$$" ] }, { "cell_type": "code", "execution_count": 4, "id": "a217974d", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAj0AAAG6CAYAAADjzPf8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABgUElEQVR4nO3deVxU9f7H8dcMO7KJCCigmLu540ZWalmaZWoaZu5p+251S3+ZWd0sy/JWpuWG2qK4Zq6lZYuaJor7vqSogIiyrzPn94e3Ka5ahMAM8H4+HvP4Nd/5nnM+A/fneXO+33O+JsMwDEREREQqOLO9CxAREREpCwo9IiIiUiko9IiIiEiloNAjIiIilYJCj4iIiFQKCj0iIiJSKSj0iIiISKWg0CMiIiKVgkKPiIiIVAoKPSIiIlIpOGzo+fHHH+nZsyc1a9bEZDKxbNmyv91mw4YNtG7dGjc3N+rVq0d0dHSp1ykiIiLlg8OGnszMTFq0aMGUKVOK1P/48ePceeeddOnShbi4OJ555hlGjhzJ2rVrS7lSERERKQ9M5WHBUZPJxNKlS+ndu/dV+7z44ousXLmSPXv22Nruu+8+Ll68yJo1a8qgShEREXFkzvYuoKRs3ryZrl27Fmrr1q0bzzzzzFW3yc3NJTc31/bearWSkpJCtWrVMJlMpVWqiIiIlCDDMEhPT6dmzZqYzVcfxKowoSchIYGgoKBCbUFBQaSlpZGdnY2Hh8dl20yYMIHx48eXVYkiIiJSik6dOkVoaOhVP68woac4Ro8ezahRo2zvU1NTqVWrFqdOncLHx8eOlYmIiEhRpaWlERYWhre391/2qzChJzg4mMTExEJtiYmJ+Pj4XPEqD4Cbmxtubm6Xtfv4+Cj0iIiIlDN/NzXFYe/e+qciIyNZv359obZvv/2WyMhIO1UkIiIijsRhQ09GRgZxcXHExcUBl25Jj4uL4+TJk8CloakhQ4bY+j/yyCMcO3aMf/3rXxw4cICPP/6YmJgYnn32WXuULyIiIg7GYYe3tm3bRpcuXWzvf597M3ToUKKjozl79qwtAAHUqVOHlStX8uyzz/Kf//yH0NBQZsyYQbdu3cq8dhERcXwWi4X8/Hx7lyFF4OLigpOT0zXvp1w8p6espKWl4evrS2pqqub0iIhUUIZhkJCQwMWLF+1divwDfn5+BAcHX3HeTlHP3w57pUdERKQ0/B54AgMD8fT01HPZHJxhGGRlZZGUlARAjRo1ir0vhR4REak0LBaLLfBUq1bN3uVIEf1+F3ZSUhKBgYHFHupy2InMIiIiJe33OTyenp52rkT+qd9/Z9cyD0uhR0REKh0NaZU/JfE7U+gRERGRSkGhR0RERCoFTWQWEREpopMnT5KcnHzVzwMCAqhVq1YZViT/hEKPiIhIEZw8eZKGDRuSk5Nz1T7u7u4cPHiwVILPsGHDmDNnzmXt3bp1Y82aNSV+vP/16quvsmzZMttKCeWRQo+IiEgRJCcn/2XgAcjJySE5ObnUrvZ0796d2bNnF2q70sLZcmWa0yMiIpVeZmbmVV9/F3SKs9/icnNzIzg4uNCratWqbNiwAVdXV3766Sdb34kTJxIYGEhiYiIAa9as4cYbb8TPz49q1apx1113cfTo0UL7j4+PZ8CAAfj7+1OlShXatGnDli1biI6OZvz48ezcuROTyYTJZCI6OrrY38NedKWnlGjcV0Sk/PDy8rrqZz169GDlypXF2m94ePgVzwUlvQJU586deeaZZxg8eDA7d+7k2LFjjB07loULFxIUFARcCmCjRo2iefPmZGRk8Morr9CnTx/i4uIwm81kZGTQqVMnQkJCWL58OcHBwWzfvh2r1Ur//v3Zs2cPa9asYd26dQD4+vqW6HcoCwo9pcDe474iIlIxrVix4rKANmbMGMaMGcMbb7zBt99+y0MPPcSePXsYOnQod999t61f3759C203a9Ysqlevzr59+2jatClffPEF586d49dff8Xf3x+AevXq2fp7eXnh7OxMcHBwKX7D0qXQUwoKjfuazLiFXo+TV1UsGRfIjd8LhrXUx31FRKToMjIyrvrZtazufeLEiWJveyVdunRh6tSphdp+Dyiurq58/vnnNG/enNq1a/P+++8X6nf48GFeeeUVtmzZQnJyMlarFbj0h3rTpk2Ji4ujVatWtv1VRAo9pcijQST+tz6Es091W1tB2jlS1n9K9qHNdqxMRET+rEqVKuViv1WqVCl09eV/bdq0CYCUlBRSUlIKHb9nz57Url2b6dOnU7NmTaxWK02bNiUvLw/4Y32rikwTmUuJR4NIqvceg5N3QKF2J+9qVO89Bo8GkXaqTEREKqKjR4/y7LPPMn36dNq3b8/QoUNtV3POnz/PwYMHefnll7n11ltp3LgxFy5cKLR98+bNiYuLIyUl5Yr7d3V1xWKxlPr3KE0KPaXAYjXwv/Uh4PK1QkwmM3Dpc4u1ZCeyiYhI6QkICMDd3f0v+7i7uxMQEPCXfa5Fbm4uCQkJhV7JyclYLBYGDRpEt27dGD58OLNnz2bXrl1MmjQJgKpVq1KtWjU+/fRTjhw5wnfffceoUaMK7XvAgAEEBwfTu3dvNm7cyLFjx1i8eDGbN18amQgPD+f48ePExcWRnJxMbm5uqX3P0qLhrVKwPzmv0JDW/zKZzDj7VGd/ch5ty7AuEREpvlq1anHw4EG73pm7Zs0aatSoUaitYcOG3H///fz222+sWLECgBo1avDpp58yYMAAbr/9dlq0aMH8+fN56qmnaNq0KQ0bNuSDDz6gc+fOtv24urryzTff8Nxzz9GjRw8KCgpo0qQJU6ZMAS5NhF6yZAldunTh4sWLzJ49m2HDhpXady0NJqOk75srx9LS0vD19SU1NRUfH59i7+c/yzby/i8X/7bfsx38eLp3x2IfR0RE/pmcnByOHz9OnTp1/vaqjTiWv/rdFfX8reGtUlDVvWg/1sQTh0u5EhEREfmdQk8puKlRTSzpyRiG9ap9LNnpfDD2yRK/nVFERESuTKGnFNQJr82b/Vr/d9LylTm5e9HrifHUrl27DCsTERGpvBR6SsngLs2YNqg1NXwLjzvW8HXnpvoBYDKxxdSIDYfO2alCERGRykV3b5Wi7k1rcFuTYLYeTyEpPYdAb3fa1bn0pMtnFsTx9c4zPPpZLNMHtmT13A8ZM2bMX67/IiIiIsWn0FPKnMwmIutWu6z9vagWZOTk8/3Bcwyd+QvxcxewY8cOli9fjouLix0qFRERqdg0vGUnLk5mPh4YQbtwf6xOrgT3f511W3YxcuTIEl99V0RERBR67MrD1YkZw9rQpIYPZk9fgu77N58vXcXo0aPtXZqIiEiFo9BjZz7uLswd0Y7rAqrg7FOdoP6v884HU/nPf/5j79JEREQqFIUeBxDg5ca8ke2p6euOS7UwAqNe49l/jWH+/Pn2Lk1ERK7AYjXYfPQ8X8WdZvPR85V+LUWTycSyZcuuaR/Dhg2jd+/eJVLP1Sj0OIgQPw/mjWyPfxVX3ILrUfO+16gZpmf4iIg4mjV7znLj298xYPovPD0/jgHTf+HGt79jzZ6zpXrcc+fO8eijj1KrVi3c3NwIDg6mW7dubNy4sVSPW5Eo9DiQutW9mPtAO7zcnHGu0Yjow87kW67+VGcRESlba/ac5dHPtnM2NadQe0JqDo9+tr1Ug0/fvn3ZsWMHc+bM4dChQyxfvpzOnTtz/vz5UjtmRaPQ42Cahvgyc2gb3JzNfHcgiecX7iR2+3aOHj1q79JERCocwzDIyiso0is9J59xy/dypYGs39teXb6P9Jz8Iu3vn9ype/HiRX766SfefvttunTpQu3atWnXrh2jR4/m7rvvBuC9996jWbNmVKlShbCwMB577DEyMjJs+4iOjsbPz48VK1bQsGFDPD096devH1lZWcyZM4fw8HCqVq3KU089hcVisW0XHh7O66+/zoABA6hSpQohISG2ldev5tSpU0RFReHn54e/vz+9evUqtOySxWJh1KhR+Pn5Ua1aNf71r3+VyZ3Lek6PA2p/XTWmDmrNQ3Nj+SruDAvmfYPXoVVs2riRwMBAe5cnIlJhZOdbaPLK2hLZlwEkpOXQ7NVvitR/32vd8HQt2mnYy8sLLy8vli1bRocOHXBzc7usj9ls5oMPPqBOnTocO3aMxx57jH/96198/PHHtj5ZWVl88MEHzJ8/n/T0dO655x769OmDn58fq1at4tixY/Tt25eOHTvSv39/23bvvPMOY8aMYfz48axdu5ann36aBg0acNttt11WR35+Pt26dSMyMpKffvoJZ2dn3njjDbp3786uXbtwdXVl0qRJREdHM2vWLBo3bsykSZNYunQpt9xyS5F+HsWlKz0O6pZGQUyKaoEJcG92O+drdKBHjx6kp6fbuzQRESljzs7OREdHM2fOHPz8/OjYsSNjxoxh165dtj7PPPMMXbp0ITw8nFtuuYU33niDmJiYQvvJz89n6tSptGrViptvvpl+/frx888/M3PmTJo0acJdd91Fly5d+P777wtt17FjR1566SUaNGjAk08+Sb9+/Xj//fevWOuCBQuwWq3MmDGDZs2a0bhxY2bPns3JkyfZsGEDAJMnT2b06NHcc889NG7cmGnTpuHr61uyP7Qr0JUeB9arZQjpOQW8vGwPvjfcx6HvZtKvXz++/vprXF1d7V2eiEi55+HixL7XuhWp79bjKQyb/evf9ose3ta25NDfHfuf6Nu3L3feeSc//fQTv/zyC6tXr2bixInMmDGDYcOGsW7dOiZMmMCBAwdIS0ujoKCAnJwcsrKy8PT0BMDT05O6deva9hkUFER4eHihJZCCgoJISkoqdOzIyMjL3k+ePPmKde7cuZMjR47g7e1dqD0nJ4ejR4+SmprK2bNnad++ve0zZ2dn2rRpU+pDXLrS4+AGdajNC90aAuB/ywg2JRg88MADWK2a4Cwicq1MJhOers5Fet1Uvzo1fN0xXW1f/L6odPUi7c9kutqers7d3Z3bbruNsWPHsmnTJoYNG8a4ceM4ceIEd911F82bN2fx4sXExsba5t3k5eXZtv/fZY5MJtMV267lHJORkUFERARxcXGFXocOHeL+++8v9n5LgkJPOfBY57o8dPN1APh3e4Kl207wr3/9y85ViYhULk5mE+N6NgG4LPj8/n5czyY4mf95mCmuJk2akJmZSWxsLFarlUmTJtGhQwcaNGjAmTNnSuw4v/zyy2XvGzdufMW+rVu35vDhwwQGBlKvXr1CL19fX3x9falRowZbtmyxbVNQUEBsbGyJ1Xs1Cj3lgMlkYvQdjbivbRgmsxMBPZ9ny2+phdK7iIiUvu5NazB1UGuCfd0LtQf7ujN1UGu6N61RKsc9f/48t9xyC5999hm7du3i+PHjLFy4kIkTJ9KrVy/q1atHfn4+H374IceOHWPevHlMmzatxI6/ceNGJk6cyKFDh5gyZQoLFy7k6aefvmLfgQMHEhAQQK9evfjpp584fvw4GzZs4KmnniI+Ph6Ap59+mrfeeotly5Zx4MABHnvsMS5evFhi9V6N5vSUEyaTiX/3aUZ6TgErd58ludE97D6bQUTtvx83FhGRktO9aQ1uaxLM1uMpJKXnEOjtTrs6/qV6hcfLy4v27dvz/vvvc/ToUfLz8wkLC+PBBx9kzJgxeHh48N577/H2228zevRobr75ZiZMmMCQIUNK5PjPPfcc27ZtY/z48fj4+PDee+/RrduV50J5enry448/8uKLL3LPPfeQnp5OSEgIt956Kz4+Prb9nT17lqFDh2I2m3nggQfo06cPqampJVLv1ZgMLeltk5aWhq+vL6mpqbZfjKPJLbAwcs42fjqcjI+7M/Mf6oBbdnKhiWkiInJlOTk5HD9+nDp16uDu7v73Gwjh4eE888wzPPPMM3at469+d0U9f2t4q5xxc3bik8ERRNSuSlpOAX0mr6d1p25lMhYqIiJSnin0lEOers7MGtqWhkFe5Jrd8L5rND36DuDIkSP2Lk1ERMRhKfSUU76eLswb2Z5aVd1x9gvGqesz3N6zD4mJifYuTUREKpATJ07YfWirpCj0lGOB3u58/mAk1b1ccA2oTXbbB7ijZ289tVlEROQKFHrKuTB/T754MBIfNzNuNRtwpm5Peve9V7ezi4j8Bd3DU/6UxO9MoacCqB/kzbyRkbg7m3Cv3YLd3m3Z8MOP9i5LRMTh/P704aysLDtXIv/U77+z/32C9D+h5/RUEC3C/Jg9vD2DZ/6CZ/0OrEkJoKvVwFyGTwYVEXF0Tk5O+Pn52daW8vT0LNZyEFJ2DMMgKyuLpKQk/Pz8cHL6Z2uW/ZlCTwUSWbcaUwe14ZHPYlmy4zQ+Hi68dHtdPYtCRORPgoODAS5bVFMcm5+fn+13V1x6OOGflIeHExbF0h3xPLtgJwBOB9byat+2DB482M5ViYg4FovFQn5+vr3LkCJwcXH5yys8RT1/60pPBdSnVShp2QWMW74XS6NuPPnhDKpXr0737t3tXZqIiMNwcnK6pqESKX80kbmCGnpDOM92rQ+A3y0jGfjyB/z66692rkpERMR+FHoqsKdurc+wyFoAeN3yCD0ffZnDhw/buSoRERH7UOipwEwmE6/0bEqv5kGYzE64dX6U2wc/SUJCgr1LExERKXMKPRWc2WxiUv/WdK7nh8nZBcsNIxk26lV7lyUiIlLmFHoqAWcnM9OGdqB1TU/Mrh6cadCbgwlaqkJERCoXhZ5Kwt3FiXkP30TLMD/SciwMnrmFk+f1RFIREak8FHoqkSpuzkQPb0vDIG+S0nPpNXkdjz03WmvQiIhIpaDQU8n4eboyd0Q7gr2cuZDvxLKLobz21iR7lyUiIlLqFHoqoSAfd2IevYkq5gJcq4czbb+JT2fPsXdZIiIipUqhp5KqVc2TJU/dgquRj1vNRryyNp7lK1fZuywREZFSo9BTiTUM9mb+YzdjtubjHt6Sh2f/wqbNv9i7LBERkVKh0FPJta7tz+wHOmCyWnCr1557Jy7lzJkz9i5LRESkxCn0CJ0aBfN+VHMwrLg0vJmZ2y/qji4REalwFHoEgN5twpnQpykAM38+zkffHbFzRSIiIiXLoUPPlClTCA8Px93dnfbt27N169a/7D958mQaNmyIh4cHYWFhPPvss+Tk5JRRteXfgA51GHtXEwAmfXuIqP+bQm5urp2rEhERKRkOG3oWLFjAqFGjGDduHNu3b6dFixZ069aNpKSkK/b/4osveOmllxg3bhz79+9n5syZLFiwgDFjxpRx5eXbiBvr8NQt9QDYagnnjkdfwWKx2LkqERGRa+ewoee9997jwQcfZPjw4TRp0oRp06bh6enJrFmzrth/06ZNdOzYkfvvv5/w8HBuv/12BgwY8JdXh3Jzc0lLSyv0Enj2tgbcWssJgKMBHbnvuX9rjo+IiJR7Dhl68vLyiI2NpWvXrrY2s9lM165d2bx58xW3ueGGG4iNjbWFnGPHjrFq1Sp69Ohx1eNMmDABX19f2yssLKxkv0g5ZTKZmP5IN1pXzcdkdmKrczOeeuNDe5clIiJyTRwy9CQnJ2OxWAgKCirUHhQUREJCwhW3uf/++3nttde48cYbcXFxoW7dunTu3Pkvh7dGjx5Namqq7XXq1KkS/R7lmdlsYsHzd1PXPROTsytfpdTg31Pn2bssERGRYnPI0FMcGzZs4M033+Tjjz9m+/btLFmyhJUrV/L6669fdRs3Nzd8fHwKveQPLk5mVv5fXwK5iNnNk2n7nZkRs8LeZYmIiBSLQ4aegIAAnJycSExMLNSemJhIcHDwFbcZO3YsgwcPZuTIkTRr1ow+ffrw5ptvMmHCBKxWa1mUXSG5uzix/tUovPNTcPL04eP9TsRfyLJ3WSIiIv+YQ4YeV1dXIiIiWL9+va3NarWyfv16IiMjr7hNVlYWZnPhr+PkdGkyribhXhtvdxfWjb2H2lXdSMm2MmjGFs6l61Z2EREpXxwy9ACMGjWK6dOnM2fOHPbv38+jjz5KZmYmw4cPB2DIkCGMHj3a1r9nz55MnTqV+fPnc/z4cb799lvGjh1Lz549beFHii/IrwrzH+lIiJ8HJ85n0e+jDRw49pu9yxIRESkyZ3sXcDX9+/fn3LlzvPLKKyQkJNCyZUvWrFljm9x88uTJQld2Xn75ZUwmEy+//DKnT5+mevXq9OzZk3//+9/2+goVTg1fDz4f2Z67P9jAb6kF3PnWCn5+4z5qBFazd2kiIiJ/y2Ro7McmLS0NX19fUlNTNan5L6zbdoARn+/C5FYF9wvH2DppBD5envYuS0REKqminr8ddnhLHFfXNo14+85wjPwccqpeR6cXZ5GXX2DvskRERP6SQo8Uy31d2/F8B18MSz4XvOtw64vTdZeciIg4NIUeKbYn7+3K4PpWDKuFU6616D0u2t4liYiIXJVCj1yTNx66h27VUgDYlR/ElO8O27kiERGRK1PokWv26YvDuLvWpZXY3/nmEJ9v0a3sIiLieBR6pER88NjdPNa5LgAvL9vDrHU77VyRiIhIYQo9UmJe6NaQ+9uFYRgw/psTzPk21t4liYiI2Cj0SIkxmUyM6VYX98Q9mMzOvLL2JCu3HrB3WSIiIoBCj5QwrypV2PDWcMwJ+zA5u/L4gr1sPnDK3mWJiIgo9EjJCw4KZPXL/bAmHgIXdwZO38y++PP2LktERCo5hR4pFQ3r1SHmiVsoOHccq0sVer2/jlPnM+xdloiIVGIKPVJqOkS0YGr/6ylIOU2+ixdRU38iOSPX3mWJiEglpdAjparn7V0Y3yWAqm5wNsPK0FlbScvJt3dZIiJSCSn0SKkbMeAeFj3RiWpVXNl7Jo0Rs7eSk2+xd1kiIlLJKPRImahb3Ys5D7TD08XMr79d5J5Jq8i3aIFSEREpOwo9UmaahvjSNnsb1vxc9l00M2TKOqxWw95liYhIJaHQI2Vq5lujue7segxLAZvP5PNE9E8YhoKPiIiUPoUeKVPOzs6smvEO/oe+xjCsrDqUziuLttm7LBERqQQUeqTMeXp6sn72RFx3LQVgXmwSk9fssXNVIiJS0Sn0iF34+/uz/pNXsexYAsDkDb+x4NeTdq5KREQqMoUesZtatWqx6t1n8E/aDsDoJbtZtfusnasSEZGKSqFH7KpZs2bEznyZAe3CsBrw9Pwd/HjonL3LEhGRCkihR+zOZDLxRu9m3NmsBvkWgwdmbSb2txR7lyUiIhWMQo84BCezicENrGQfi6UAJ+7/ZCP7z6bZuywREalAFHrEYXRo24aRjSEnfi+5VjNRH//IieRMe5clIiIVhEKPOJQ3Xh1LD6/fyEs8Rnq+iX5TfiAhNcfeZYmISAWg0CMOxWQy8elH/6FV2ibyU86QnG1w75QfSMnMs3dpIiJSzin0iMNxdnZm0WezCDu6lIL0ZE6lFTBo+iYycgvsXZqIiJRjCj3ikDw8PFi18DN8YufgbipgX0ImD87ZRk6+xd6liYhIOaXQIw7L39+f2O9XEvN4J7zcnNl87DxPfLGDAovV3qWJiEg5pNAjDs3Dw4PmoX5MH9IGVycT6/Yn8vzCnVitWpldRET+GYUeKRc6XOePz56FGFYLy+LOMP7rvRiGgo+IiBSdQo+UCyaTibEP9CZl9X8AmLP5NyavO2znqkREpDxR6JFyo1evXrz31H2c/2YqAP9Zf5hZPx+3c1UiIlJeKPRIufLggw/yQu92XPxxHgCvrdjHwm2n7FyViIiUBwo9Uu6MHTuW+1pUJe3XZQD8a9Eu1u5NsG9RIiLi8BR6pNwxmUxM+egjOvulYDn8Mwbw5Bc72Hgk2d6liYiIA1PokXLJycmJLz7/nB8mPUL364PJs1h5cO42dpy8YO/SRETEQSn0SLnl7u7OdeG1+c+AltxYL4CsPAtDZ23hYEK6vUsTEREHpNAj5Z6bsxN9g8+Td/YQaTkWBs34hZPns+xdloiIOBiFHqkQGtWtQ+6375N37gTnMvIYOOMXktJy7F2WiIg4EIUeqRAaN27MisULSPvq3+RfOMupC9kMnrmFi1l59i5NREQchEKPVBiRkZF8OfsTkheNoyD9PAcTMxg2+1cycwvsXZqIiDgAhR6pUHr27MnHE18jKeYVLNnpxJ26yEPztpFbYLF3aSIiYmcKPVLhjBgxgnHPPEjSwnGYrflsPHKep7+Mo8BitXdpIiJiRwo9UiGNGTOGeZNfJ3rkDbg6mVmzN4HRS3ZjtWpldhGRykqhRyokk8lEVFQUNzcI5IMBrTCbYGFsPP9etR/DUPAREamMFHqkwru1YTXqX9gCwMyfj/PRd0fsXJGIiNiDQo9UeNnZ2SRs/oqUdZ8CMOnbQ8zZdMK+RYmISJlT6JEKz8fHh9WrV+OfvJOLG78AYNzyvSzbcdrOlYmISFlS6JFKoWbNmqxduxbn/WtJ27YcgOcW7mTdvkQ7VyYiImVFoUcqjYYNG7Jy5UpyN39Oxp7vsFgNHvtiO5uPnrd3aSIiUgYUeqRSad++PTExC7j4zUdkH9lKXoGVB+duY1f8RXuXJiIipUyhRyqdO++8k9kzZxD94I1EXleNjNwChs7aypGkdHuXJiIipchk6KElNmlpafj6+pKamoqPj4+9y5EykJFbwMDpv7AzPpVgH3cWPhJJmL+nvcsSEZF/oKjnb13pkUrNy82ZMTdWxZSWQEJaDoNnbuFceq69yxIRkVKg0COV3uxpH3Jq3otY0pI4cT6LIbO2kpqdb++yRESkhCn0SKX3n//8h07tWpDw5f9hZKWy/2waD0T/SlZegb1LExGREqTQI5Wem5sbS5cu5fragZyd/3+Ql0Xsbxd45LPt5BVoZXYRkYpCoUeEP57aHFIFzi54BSx5/HjoHM8uiMOildlFRCoEhR6R/6pRowZr167FJy+ZxEWvYzKsrNx9lv9bulsrs4uIVAAKPSJ/0qBBA1atWkWnhoG83bsRZhPM//UUb605YO/SRETkGjnbuwARR9O2bVvWrFkDgOHsxouLd/PJD8fw9XDhsc717FydiIgUl0Nf6ZkyZQrh4eG4u7vTvn17tm7d+pf9L168yOOPP06NGjVwc3Oz/dUuUlz929aik++ltbkmrjnIZ7/8ZueKRESkuBw29CxYsIBRo0Yxbtw4tm/fTosWLejWrRtJSUlX7J+Xl8dtt93GiRMnWLRoEQcPHmT69OmEhISUceVSkfz000/MHTOU1E0LABj71R6W7zxj56pERKQ4HHYZivbt29O2bVs++ugjAKxWK2FhYTz55JO89NJLl/WfNm0a77zzDgcOHMDFxaVIx8jNzSU394+n76alpREWFqZlKMTGMAxGjRrF5MmTqdbtMbxa9sDZbGL6kDZ0aRRo7/JERIRyvgxFXl4esbGxdO3a1dZmNpvp2rUrmzdvvuI2y5cvJzIykscff5ygoCCaNm3Km2++icViuepxJkyYgK+vr+0VFhZW4t9FyjeTycSkSZO47777OL92KrmHfqbAavDIZ7FsPZ5i7/JEROQfcMjQk5ycjMViISgoqFB7UFAQCQkJV9zm2LFjLFq0CIvFwqpVqxg7diyTJk3ijTfeuOpxRo8eTWpqqu116tSpEv0eUjGYzWaio6O59dZbSPjqHSwn48gtsDIi+lf2nE61d3kiIlJEDhl6isNqtRIYGMinn35KREQE/fv35//+7/+YNm3aVbdxc3PDx8en0EvkStzc3FiyZAmtWjTn9MLX4NwR0nMLGDprK0fPZdi7PBERKQKHDD0BAQE4OTmRmJhYqD0xMZHg4OArblOjRg0aNGiAk5OTra1x48YkJCSQl5dXqvVK5fD7U5uvqx3GQ42sNA3x4XxmHoNnbOHMxWx7lyciIn/DIUOPq6srERERrF+/3tZmtVpZv349kZGRV9ymY8eOHDlyBKv1j7WSDh06RI0aNXB1dS31mqVyCAoKYteuXYx54VnmDG/HddWrcCY1h0Ezt5Cckfv3OxAREbtxyNADMGrUKKZPn86cOXPYv38/jz76KJmZmQwfPhyAIUOGMHr0aFv/Rx99lJSUFJ5++mkOHTrEypUrefPNN3n88cft9RWkgvL09ASgmpcbU6Oa4GXO49i5TIbO2kpaTr6dqxMRkatx2Ccy9+/fn3PnzvHKK6+QkJBAy5YtWbNmjW1y88mTJzGb/8hsYWFhrF27lmeffZbmzZsTEhLC008/zYsvvmivryAVnMViYXDfuzh45DTXPfgRe8+kMTJ6G3NHtMPdxenvdyAiImXKYZ/TYw9Fvc9f5Heff/45gwYNwiXwOsIfeJ88w4lbGgXyyeAIXJwc9kKqiEiFUq6f0yNSXgwcOJB3332X/KRjnJw3GheTwXcHknguZicWq/6eEBFxJAo9Itfoueee47nnniP39D7OLnodswmW7zzDuOV70IVUERHHodAjUgImTpzIwIEDyTyylYurJ2MCPvvlJO9+c9DepYmIyH8p9IiUALPZzKxZs7j99tvxTN7P4+39AZjy/VE+/fGonasTERFw4Lu3RMobV1dXFi1axMWLFwkLC6NK1aO8veYAb646gI+7C/e1q2XvEkVEKjVd6REpQd7e3raFax/tXJee9d0BGLN0N6t2n7VnaSIilZ5Cj0gp+eqrr5j++J34JO3CasDT83fww6Fz9i5LRKTSUugRKSVhYWG4uLiwO/pl/DNOkG8xeGReLLG/pdi7NBGRSkmhR6SUtG7dmqVLl+Li7MSOqc8QaEkmO9/C8Nm/sv9smr3LExGpdBR6REpR165dmTNnDlgL2Db5YYKdMknLKWDwzK0cT860d3kiIpWKQo9IKRswYADvvfceRkEuWyc9QLBbAckZuQyasYWE1Bx7lyciUmko9IiUgWeffZYXXngBIzeT2ie+JryaJ6cvZjNo5hZSMvPsXZ6ISKWg0CNSRt566y0+++wzvpz9KZ+NbE+wjztHkjIYNnsr6Tn59i5PRKTCU+gRKSNms5mBAwdiNpsJrerJ3Afa4ufhzK74VB6cu42cfIu9SxQRqdAUekTsIDc3l3HPPszFr/6Np6uZX46l8MQXO8i3WO1dmohIhaXQI2IH6enpbNu2jd92/IT5p09wdTKxbn8iLy7ahdWqldlFREqDQo+IHQQEBLB27VqCgoLYu+Erqu5bjJPJxJIdp3ltxT4MQ8FHRKSklWjoSUlJwWrV5XmRorjuuutYvXo13t7ebF02i9oJP2AyQfSmE7y/7rC9yxMRqXCuOfTs27ePt956ixtuuIHq1asTGBjIkCFDWLx4MZmZeviayF9p1aoVy5Ytw8XFhe+jJ3J9zj4APlh/mJk/H7dzdSIiFUuxQs/Bgwd57rnnqF+/Ph06dODXX3/lkUceITExkVWrVlG7dm1ee+01AgICuOOOO5g6dWpJ1y1SYdxyyy3MmzcPk8nEtx+/zJBW/gC8vmIfC7edsnN1IiIVh8koxuSB2bNns2XLFnr16sWtt96Kq6vrFfudOHGCr776ihUrVvDtt99ec7GlLS0tDV9fX1JTU/Hx8bF3OVLJTJs2jcaNG3PzzTfz75X7mfHzccwm+HhgBN2bBtu7PBERh1XU83exQk9FpdAjjsIwDF5YuJNF20/j6mRm9vC2dKwXYO+yREQcUlHP3/94eCs7O5vTp09f1r53795/uisRuYq9e/fy1dj7aFfDlTyLlQfnbmPHyQv2LktEpFz7R6Fn0aJF1K9fnzvvvJPmzZuzZcsW22eDBw8u8eJEKqs333yTI4cOsv6NQbSu6UFWnoVhs3/lYEK6vUsTESm3/lHoeeONN4iNjSUuLo7Zs2czYsQIvvjiCwA9V0SkBE2fPp127dqRkpxE7H8e4fpgT1Kz8xk8cwsnz2fZuzwRkXLpH4We/Px8goKCAIiIiODHH3/kk08+4bXXXsNkMpVKgSKVUZUqVVi5ciUNGjTg5LHDnJz3EvWqe5KUnsugmVtISsuxd4kiIuXOPwo9gYGB7Nq1y/be39+fb7/9lv379xdqF5Fr9/tTm4ODg9mzfSvZq9+hVlUPTqZkMXjmVi5m5dm7RBGRcuUfhZ558+YRGBhYqM3V1ZUvv/ySH374oUQLExEIDw9nzZo1+Pj4sGn9apql/ECgtxsHE9MZNvtXMnML7F2iiEi5oVvW/0S3rIuj2rBhA2+//Tbz588nMcdM1CebuZiVT8d61Zg1rC1uzk72LlFExG5K7Zb1373xxhusWrWKxMTE4u5CRIqoc+fOrFq1Cl9fXxoEeRM9vB1VXJ3YeOQ8T325gwKL1rwTEfk7xb7SYzabbZOXg4ODad26NREREbb/GxISUqKFlgVd6ZHyYtKkSVxwDSQmIYA8i5V+EaFM7Nscs1k3FIhI5VPqT2Ru3749Z8+eZfjw4QQEBLB9+3ZiY2M5cOAAFouF6tWr07p1a1atWlXsL1HWFHqkPFi9ejU9evTAbDYzfuYyZh9ywmI1GHFjHV6+s7HupBSRSqfUh7e2bNnCa6+9xvTp01m3bh0vv/wyu3btIj09nU2bNjFu3DhCQ0OLu3sRuYru3bszfPhwrFYr/340igdbeAAw8+fjfPjdETtXJyLiuK55InNGRgavvfYan3zyCY899hhjx47F09OzpOorU7rSI+VFQUEBvXv3ZuXKlVStWpUXpq/ik1/PAzD+7usZekO4fQsUESlDpX6l53deXl5MnDiRbdu2sWfPHurVq8fcuXOvdbci8hecnZ2JiYmhQ4cOXLhwganPRjG87aUHh45bvpelO+LtXKGIiOO55tADl/7qzM3NZcCAAYSGhjJ8+HBSUlJKYtcichWenp58/fXXNGzYkFOnThEzbhiD218aUn5+4S7W7dOdlSIif+Zc3A3feustdu/eze7duzlw4ADu7u40b96cdu3a8fDDD+Pr61uSdYrIFfz+1OYbb7yRRx5+mMd7NScz32DJ9tM89sV25gxvR2TdavYuU0TEIVzTLevh4eEMHTqUAQMG0KBBg5KurcxpTo+UVxkZGXh5eQFQYLHy6Ofb+XZfIlVcnfjyoQ40D/Wzb4EiIqWo1Of03HTTTZw/f57x48cTERFBx44defLJJ5k9ezY7d+7EYrEUd9ci8g/9HngA0tNSaZD8MzfUrUZmnoWhs7ZyJCndjtWJiDiGa7576/Dhw8TGxrJ9+3bb6+LFi7i5udGsWTO2bt1aUrWWOl3pkfIuLy+Ptm3bsmvXLv5v3GvsDujCzvhUgn3cWfhIJGH+5fPOShGRv1LqDyf8K8ePH2fbtm3s2LGDN998s6R3X2oUeqQimDp1Ko899hgA7388nZXZ9TmclEF4NU9iHokk0NvdzhWKiJSsUg09CQkJVK1aFTc3tyL1P3bsGNddd90/PUyZU+iRimLs2LG88cYbmM1mZn65mBnHfYi/kE2jYG8WPByJr4eLvUsUESkxpTqnZ9GiRfj7+9OnTx9mz57NuXPnLuuzZcsWxowZw/XXX0+LFi2KcxgRKabXXnuNkSNHYrVaeXToAJ5t6UR1bzcOJKTzQPSvZOUV2LtEEZEyV+zhrSNHjrB8+XK++uorfvnlF9q2bUuPHj04fvw4K1asAODOO++kV69e3Hbbbbi7O/4ldV3pkYqkoKCAe+65h6+//ho/Pz/mff0dY9YlkZZTwM0NqjNjSBtcnUvkUV0iInZVpnN6zp8/z4oVK1i1ahXh4eH06tWLyMjIcrfwoUKPVDRZWVl07dqV06dPs3btWjI9ghg0YwvZ+RbubFaDDwa0wkkrs4tIOWfXiczllUKPVEQpKSlkZ2cTEhICwE+Hz/FA9K/kWwzuaxvGhHualbs/UERE/qzM1t4SEcfm7+9vCzwAzslHeOeeJphNMP/XU7y1+gD620dEKgOFHpFKZOHChXTq1InZrz3Nv3tfD8AnPx5j6g9H7VyZiEjpU+gRqUSCg4Mxm80sX76c72b8m//r0QiAiWsO8tkvv9m5OhGR0qXQI1KJ3HTTTXz55ZeYzWZmzJjBqXVzeKJLPQDGfrWHr+JO27lCEZHSo9AjUsn06dOHjz/+GIDXX38dz6PfMSSyNoYBz8Xs5PsDSXauUESkdCj0iFRCDz/8MOPGjQPgiScep3nBIXq1rEmB1eCRz2LZejzFzhWKiJQ8hR6RSmrcuHE89NBDGIbBN2vX8O69Lbi1USC5BVZGRP/KntOp9i5RRKREKfSIVFImk4kpU6YwZ84cPvnkE1yczEwZ2Jr2dfxJzy1g6KytHD2XYe8yRURKjEKPSCXm7OzMkCFDbA8ndDHDW3fWoWmID+cz8xg8YwunL2bbuUoRkZKh0CMiAOTk5BAVFcUdXTszuXd96lavwpnUHAbP2EJyRq69yxMRuWYKPSICXFquYuvWrRw8eJDBUX349P7mhPh5cCw5k6GztpKWk2/vEkVErolCj4gAULNmTdauXUvVqlX55ZdfeOahoUQPiyDAy5W9Z9IYGb2N7DyLvcsUESk2hR4RsWnSpAlff/017u7urFixgrf+bxRzhrfD292ZrSdSeOzzWPItVnuXKSJSLAo9IlJIx44dWbBgAWazmVmzZvHl1HeYNawt7i5mvj94judidmKxaoFSESl/FHpE5DJ3330306ZNA2DSpElUJ42pgyJwNptYvvMMr3y1Ryuzi0i542zvAkTEMT344IOkpqYSGRlJeHg44cD7/Vvy1PwdfL7lJH6eLrzQrZG9yxQRKTKToT/XbNLS0vD19SU1NRUfHx97lyPicCwWCwu2nWbM0t0AjOnRiIdurmvnqkSksivq+duhh7emTJlCeHg47u7utG/fnq1btxZpu/nz52Mymejdu3fpFihSiezcuZOmTZtyvfsFXux+6QrPm6sOMH/rSTtXJiJSNA4behYsWMCoUaMYN24c27dvp0WLFnTr1o2kpL9eAfrEiRM8//zz3HTTTWVUqUjlMHbsWA4cOMAdd9zBHeFOPNLp0hWe0Ut3s3LXWTtXJyLy9xw29Lz33ns8+OCDDB8+nCZNmjBt2jQ8PT2ZNWvWVbexWCwMHDiQ8ePHc9111/3tMXJzc0lLSyv0EpErmzt3Lk2bNuXs2bN069aNEW2qMaBdLQwDnlmwgx8OnbN3iSIif8khQ09eXh6xsbF07drV1mY2m+natSubN2++6navvfYagYGBjBgxokjHmTBhAr6+vrZXWFjYNdcuUlH5+fmxZs0aatWqxaFDh7jrrrsYfVsd7mpeg3yLwSPzYon9LcXeZYqIXJVDhp7k5GQsFgtBQUGF2oOCgkhISLjiNj///DMzZ85k+vTpRT7O6NGjSU1Ntb1OnTp1TXWLVHQhISGsWbMGf39/tm7dyoD7+vN2n+vp3LA62fkWhs3+lX1ndMVURByTQ4aefyo9PZ3Bgwczffp0AgICirydm5sbPj4+hV4i8tcaN27MihUr8PDwYNWqVbz37kSmDoygTe2qpOcUMGTWVo4nZ9q7TBGRyzhk6AkICMDJyYnExMRC7YmJiQQHB1/W/+jRo5w4cYKePXvi7OyMs7Mzc+fOZfny5Tg7O3P06NGyKl2kUoiMjCQmJoYePXrwzDPP4OHqxMxhbWlSw4fkjFwGzdjC2dRse5cpIlKIQ4YeV1dXIiIiWL9+va3NarWyfv16IiMjL+vfqFEjdu/eTVxcnO11991306VLF+Li4jRXR6QU3HXXXaxYsQJvb28AfD1cmPNAO+oEVOH0xWwGz9xKSmaenasUEfmDwz6RedSoUQwdOpQ2bdrQrl07Jk+eTGZmJsOHDwdgyJAhhISEMGHCBNzd3WnatGmh7f38/AAuaxeRkmMymQAwDIO33nqL6667jnkjenLvtM0cScpg2OytfD6yPd7uLnauVETEgUNP//79OXfuHK+88goJCQm0bNmSNWvW2CY3nzx5ErPZIS9UiVQ6S5YsYcyYMbi6urJ69WrmjWhP1Ceb2RWfyoNztxE9vB3uLk72LlNEKjktQ/EnWoZCpHgsFgv33XcfixYtwtvbmx9//BGngDoMmP4LGbkFdG0cyNRBEbg46Q8VESl5FWIZChEpH5ycnJg3bx6dO3cmPT2dO+64A6/8FGYMbYObs5l1+5P416JdWK36G0tE7EehR0RKhLu7O8uWLaN58+YkJCTQrVs36npb+Xhga5zNJpbuOM34r/eii8siYi8KPSJSYnx9fVm9ejW1a9fm8OHD3Hnnndxcz59JUS0wmWDO5t94f91he5cpIpWUQo+IlKiaNWuydu1aAgMDeeCBB3BxcaFXyxBe63XpTsoP1h9m5s/H7VyliFRGDnv3loiUXw0bNuTw4cOFJhQO7lCb1Kw83v3mEK+v2IePuzP3ttEztESk7OhKj4iUij8HnuTkZKZOncrjXerx4E11AHhx8S7W7LnyWnoiIqVBV3pEpFRlZ2dz0003ceDAAXJzcxnz9NOkZRewYNspnvpyB7OGteXG+kVfM09EpLh0pUdESpWHhwdDhw4F4Nlnn2XBggW8eU8z7mgaTJ7FykPztrHj5AU7VykilYFCj4iUuhdffJGnnnoKuLSEzPffrWfyfS25qX4AWXkWhs3+lYMJ6XauUkQqOoUeESl1JpOJ999/n6ioKPLz8+nTpw97d+3kk8ERtK7lR2p2PoNnbuHk+Sx7lyoiFZhCj4iUCbPZzNy5c7nlllvIyMjgjjvu4Oyp35g9rB2Ngr1JSs9l4MxfSEzLsXepIlJBKfSISJlxc3Nj6dKltGjRgipVqmAYBr6eLswd0Y7a1Tw5lZLNkJlbuZiVZ+9SRaQC0oKjf6IFR0XKRkLCpVvVg4ODbW2nUrLoN20TiWm5tAzz4/OR7aniphtMReTvacFREXFYwcHBhQLP5s2bCfJyZt6I9vh5uhB36iIPzdtGboHFjlWKSEWj0CMidvX5559z880388ADD1CvehXmDG9HFVcnNh45z1Nf7qDAYrV3iSJSQSj0iIhdVatWDbgUfl588UVahPkxfWgbXJ3NrN2byEtLdmO1ahReRK6dQo+I2FX37t2ZNWsWAO+++y7vvfceN9QN4KMBrXAym1gUG88bK/ej6Ycicq0UekTE7gYPHszbb78NwHPPPcfnn3/O7dcHM7FvcwBmbTzOh98dsWeJIlIBKPSIiEN44YUXeOaZZwAYNmwY33zzDX0jQhnXswkA7317iOiNx+1YoYiUdwo9IuIQTCYTkyZNYsCAARQUFLBy5UoAhneswzNd6wPw6tf7WLoj3p5likg5podgiIjDMJvNREdHc9tttzFs2DBb+9O31ic1O5/ZG0/w/MJdeLm5cFuTIPsVKiLlkh5O+Cd6OKGI4ykoKCAtLQ0/v6q8sGgXi7fH4+psJnp4W9rXqcbW4ykkpecQ6O1Ouzr+OJlN9i5ZRMpYUc/futIjIg4rKyuLAQMGEB8fz4YNG3i7bzPSc/L5Zl8iw2f/ipe7M+cz/liyooavO+N6NqF70xp2rFpEHJXm9IiIw0pMTGTz5s1s376de+65B6ulgA8GtKJhkBe5BdZCgQcgITWHRz/bzpo9Z+1UsYg4MoUeEXFYderUYdWqVVSpUoV169YxbNgwnExwMTv/iv1/H6sf//U+LHqgoYj8D4UeEXFobdq0YcmSJTg7O/Pll1/ywIv/JjEt96r9DeBsag5bj6eUXZEiUi4o9IiIw7v99tuJjo4GYMnq9UXaZuXuM2TlFZRiVSJS3mgis4iUCwMHDiQhIYH/+8+cIvX/7JeTLN1+mjua1aBv61Da1/HHrDu7RCo13bL+J7plXcTxPTtqFItyW2D28sdkuvxitWFYMfKyCQ3w40z6H3N/Qqt6cE/rUPq2DqF2tSplWbKIlLKinr81vCUi5crgQYM4v+4TwIRhWAt9dum9ieRVk/mwezUWPhLJfW3D8HZzJv5CNh+sP0yndzZw77RNzN96krScK0+IFpGKSVd6/kRXekQc3/bt24mIiMCjQST+tz6Es09122cFaedIWf8p2Yc2ExsbS+vWrQHIzrPwzb4EFsXGs/FIMr/f2OXmbKbb9cH0iwilY70APdhQpJzSwwlFpELLPrSZ04e34BZ6PU5eVbFkXCA3fi/8z9UfAA9XJ3q1DKFXyxASUnNYuuM0i7fHcyQpg+U7z7B85xmCfdzp3SqEfhEh1Av0tsM3EpHSpis9f6IrPSKO7/crPX+nb9++xMTEYDZfeRTfMAx2xaeyKDae5TvPkPqnZ/+0CPWlX0QoPVvUxM/TtcRqF5HSoSs9IlKpHT16tFDg2blzJ02bNsXJyQm4tKp7izA/WoT58fJdjflufxKLt8fz/cFz7IxPZWd8Kq+v2M+tjQPp2zqUTg2r4+KkaZAi5ZlCj4hUSA8//LDtv5OTk4mIiKB69er069ePqKgoOnbsaAtFbs5O3NGsBnc0q0FyRi5fxZ1hUWw8+8+msXpPAqv3JBDg5crdLULoFxFKk5q6EixSHml46080vCXi+Io6vPXnicw//PADvXv35uLFi7bPa9asyb333ktUVBQdOnS44jDYvjNpLN4ez1dxp0n+0zpfjWv40Lf1pTlC1b3drv1Licg1Ker5W6HnTxR6RBzfyZMnadiwITk5OVft4+7uzsGDB6lVq5atLS8vj3Xr1hETE8OyZctITU21fTZr1iyGDx9+1f3lW6z8eOgci7fHs25fEnmWS5OlncwmOjeoTr+IUG5pHIibs1MJfEMR+acUeopBoUekfDh58iTJyclX/TwgIKBQ4Plfubm5fPPNN8TExLBy5UoOHjxI9eqXbn3/4osv2LFjB1FRUbRp0waTqfBt7Bez8vh65xkWbT/NzlMXbe2+Hi7c3aImfSNCaRHqe9l2IlJ6FHqKQaFHpPLJz8/HxcXF9v7mm2/mp59+Ai6t8h4VFUVUVBStWrW6LMgcSUpn8fbTLN1+moS0P6481a1ehb4RodzTKpRgX/ey+SIilZhCTzEo9IjI0qVLWbBgAV9//TVZWVm29nr16jFw4EBeffXVy7axWA02HU1mUWw8a/cmkJN/afjLZIIb6wXQLyKU25sE4+Gq4S+R0qDQUwwKPSLyu8zMTFatWsWCBQtYuXIlOTk59OjRg5UrV9r6HD16lLp16xbaLj0nn1W7z7I49jRbT6TY2r3cnLmzWQ36RoTSNryqhr9ESpBCTzEo9IjIlWRkZLBixQqCgoLo0qULAL/99hvh4eE0adLENgTWuHHjQtv9dj6TJdsvPf05/kK2rb2Wvyf3tA6hb+tQwvw9y/S7iFRECj3FoNAjIkW1ZMkSBgwYQF7eH7eyN2vWzBaAGjRoYGu3Wg22nkhhcWw8q3afJTPPYvusfR1/+kaE0qNZDbzc9Og0keJQ6CkGhR4R+SdSU1P56quviImJ4ZtvviE//4+lLJYvX07Pnj0v2yYrr4C1ey8tfrrp6Hl+/xfYw8WJ7k0vLX4aeV01zFr8VKTIFHqKQaFHRIrrwoULLFu2jJiYGDZu3Mjp06fx9r60cOm8efM4e/Ys9957L3Xq1LFtc+Zi9qXFT2PjOZacaWuv6etOn/8Of11X3avMv4tIeaPQUwwKPSJSErKysvD0/GOuTsuWLdm5cycAbdu2pX///tx77722ZwkZhsGOUxdZHBvP1zvPkJZTYNu2VS0/+kWEclfzmvh6uCAil1PoKQaFHhEpaVarlRkzZrBgwQI2bNiA1Wq1fdahQweGDh3KI488YmvLybewbn8ii2Pj+fFwMhbrpX+iXZ3N3NYkiH6tQ7mpfgDOWvxUxEahpxgUekSkNCUmJrJkyRJiYmL44YcfMAyD++67jy+//LJQn6CgIACS0nP4aselxU8PJqbb+lT3dqN3y0tPf24UrH+rRBR6ikGhR0TKytmzZ1m8eDEtWrTgpptuAmDPnj00b96cm266iaioKPr27UtwcDCGYbD3TBqLYuNZvvMMKZl/3DHWNMSHvq1DubtFTap5afFTqZwUeopBoUdE7OmTTz4pNNRlNpvp1KkTUVFR3HPPPQQGBpJXYGXDwSQWb4/nuwNJ5Fsu/RPubDbRpVEgfVuHckujQFydNfwllYdCTzEo9IiIvZ08eZJFixYRExPDli1bbO1ms5mNGzfSoUMHW1tKZh7L406zePtpdp/+Y9X4qp4u9Gp56e6vpiE+evqzVHgKPcWg0CMijuTEiRMsXLiQmJgYjh49SkJCAq6urgBER0djsVjo06cP/v7+HEpMZ3FsPEt3nCYpPde2jwZBXvRtHUqfViEE+mjxU6mYFHqKQaFHRBxVSkoK/v7+wKVb3OvVq8exY8dwdnbmtttuIyoqit69e+Pl7cPPRy4tfvrNvkTyCi7dLWY2wc0NqtO3dSi3NQnC3UWLn0rFodBTDAo9IlIe5OXl8e677xITE2N7/g+Ai4sL3bp1Y/jw4dxzzz2kZuezctdZFm+PJ/a3C7Z+3u7O3NW8Jv0iQmhdS4ufSvmn0FMMCj0iUt4cOHDANgS2Z88eAB566CE++eQT4NJzgjIzM0nONbNkezxLtp/m9MU/Fj+tE1CFe1qFcE9EKCF+Hnb5DiLXSqGnGBR6RKQ827t3LwsXLqR79+62Cc+bNm3illtuoUePHvTv358ed97J7oQcFm2PZ/XuBLLzLy1+ajJB5HXV6Ns6lDuaBePpqsVPpfxQ6CkGhR4RqWjGjx/Pq6++anvv4eHBXXfdRVRUFJ263s4PR9NYFHuKX46l2Pp4ujpxR9Ma9IsIpX0dfy1+Kg5PoacYFHpEpKIxDIOdO3cSExPDggULOHbsmO2zKlWqsG3bNho1asSplKxLi59uj+e381m2PiF+HvRtHcI9rUMJD6hij68g8rcUeopBoUdEKjLDMNixYwcLFiwgJiaGvLw8Tp06hdl86UGG0dHR+PlVJbBJe77ek8SKnWdJz/1j8dM2tavSNyKUO5vXwMddi5+K41DoKQaFHhGpLAzD4PTp04SGhgJQUFBASEgISUlJ+Pj40KtXL3r3jYKQZizblcjPh8/x37VPcXM20+36YPpGhHJjvQCcNPwldqbQUwwKPSJSWaWlpTFu3DgWLlzI6dOnbe2+vr706dOHXvcNIcG9Fotj4zmclGH7PMjHjd6tQujXOpT6Qd72KF1Eoac4FHpEpLKzWq1s3ryZmJgYFi5cyNmzZwF44YUXmDhxIoZhEHfyAou3x7NidwIXs/Jt27YI9aVvRCg9m9ekahVXe30FqYQUeopBoUdE5A8Wi4WNGzcSExPDyJEjadmyJQBr1qxh0KBB9Ol7L/U79eFgnh8/HEqm4L/jXy5OJm5tFETfiFA6N6yOi5MWP5XSVSFCz5QpU3jnnXdISEigRYsWfPjhh7Rr1+6KfadPn87cuXNtD+eKiIjgzTffvGr/K1HoERH5e48//jgff/yx7X316tW5q98Agtrdxc5UN/adTbd9Vq2K66XFTyNCuL6mrz3KlUqg3IeeBQsWMGTIEKZNm0b79u2ZPHkyCxcu5ODBgwQGBl7Wf+DAgXTs2JEbbrgBd3d33n77bZYuXcrevXsJCQkp0jEVekRE/l5BQQE//PADCxYsYMmSJZw/f972WVBQEAu/3cQPv+WwLO4MyRl/LH7aKNibfhGh9GoZQnVvN3uULhVUuQ897du3p23btnz00UfApXHmsLAwnnzySV566aW/3d5isVC1alU++ugjhgwZcsU+ubm55Ob+8f+QaWlphIWFKfSIiBRRfn4+3333HTExMSxZsoSgoCD279+PyWSiwGLl1U9iOFzgz47EAvIsl043TmYTnf67+OmtjQO1+Klcs3IdevLy8vD09GTRokX07t3b1j506FAuXrzIV1999bf7SE9PJzAwkIULF3LXXXddsc+rr77K+PHjL2tX6BER+efy8vI4efIk9erVAyA7O5vAwEAyMjIIva4BEX0eIqN6U45c+OPZP74eLvRsUYO+rUNpGeanxU+lWIoaehxydllycjIWi4WgoKBC7UFBQSQkJBRpHy+++CI1a9aka9euV+0zevRoUlNTba9Tp05dU90iIpWZq6urLfDApRNRnz598Pb2Jv7YIb6a9DzrX+qOsWI8jay/Uc3DidTsfD775SR9Pt5E1/d+YMr3Rzibmv0XRxEpPocMPdfqrbfeYv78+SxduhR3d/er9nNzc8PHx6fQS0RESkZQUBBz584lKSmJZcuWMWDAAKpUqcLJvb+y9p3H6WXexrwR7ejdsibuLmaOnsvknbUHueGt7xg8cwvLdpwmO89i768hFYhDLqMbEBCAk5MTiYmJhdoTExMJDg7+y23fffdd3nrrLdatW0fz5s1Ls0wRESkCd3d3evXqRa9evcjOzmb16tXExMRwX/8oGtSvzk31qzP3yxjGz1xOtTY9SDJ8+elwMj8dTsbLzZkezYLp2zqUdnX8Nfwl18Qh5/TApYnM7dq148MPPwQuTWSuVasWTzzxxFUnMk+cOJF///vfrF27lg4dOvzjY+ruLRER+7jvvvtYsGABAM6+QdTqHIVHo05k8MfV+lr+ntzTOoS+rUMJ8/e0V6nigMr1RGa4dMv60KFD+eSTT2jXrh2TJ08mJiaGAwcOEBQUxJAhQwgJCWHChAkAvP3227zyyit88cUXdOzY0bYfLy8vvLy8inRMhR4REfvIyMhg5cqVxMTEsGrVKnJycgATbqFNCL2pL671Isn601BXuzr+9GsdSo/mNfByc8hBCylD5T70AHz00Ue2hxO2bNmSDz74gPbt2wPQuXNnwsPDiY6OBiA8PJzffvvtsn2MGzeOV199tUjHU+gREbG/9PR0vv76a2JiYli9ejUtWrRgw8+bWLs3gcWxp9l45BwGl4a5PFyc6N700vBXZN1qWvy0kqoQoaesKfSIiDiW1NRUTp8+TZMmTQC4ePEiNeo2wbXBjVSL6IHVq7qtbw1fd/q0CqFvRCh1qxftCr9UDAo9xaDQIyLi2I4ePcqTTz7Jt99+S0FBAa41GuDVrCs+TbtguHjY+rUM86Pffxc/9fV0sWPFUhYUeopBoUdEpHw4f/48y5YtIyYmhvXr12PBjGe9drTr/xQn872w/HfxU1dnM7c1DqJvRAg316+OsxY/rZAUeopBoUdEpPxJTk5myZIlxMTEMHfuXJy9/Vked4ZPv91FUt4fV3kCvNzo06omfSNCaRSsf+MrEoWeYlDoERGpOG67/XZ+2HkUr2a3UqVJZ5w8/1jl/fqaPvRtHUqvljWp5qXFT8s7hZ5iUOgREak4zp49y+LFi4mJieHnTZtxr9OaKk1vxbNeO0xOl64AOZtNdG4YSL+IEG5pFISrs4a/yiOFnmJQ6BERqZhOnz7N4sWLWbBgAb9s302r3g9Ss2MfdsWn2vr4ujvT+793fzUL8dXTn8sRhZ5iUOgREan4Tp06xYULF2jevDmHEtOZ8+MB5vywH2evarY+11XzIKpdbfq0CiHI5+prOIpjUOgpBoUeEZHKZ/v27Tz0yKPsTS7Aq+mteNTvgNnl0jwfEwaRdapyX4c63N4kCHcXJztXK1ei0FMMCj0iIpXXsWPHiImJYf6Srzia50uVprfgHnq97XNvd2fual6Dvq1DiahdVcNfDkShpxgUekREBODw4cMsXLiQhWt+oN/zE1m17zynL2bbPq/mauG+DnUYEFmX0Kpa/NTeFHqKQaFHRESuxGo1+OX4eR54/VOy/Btgdv3j6c/XVcnngVuup0+bOlTR4qd2odBTDAo9IiLyV/bt28fnCxYRs/kQ6dWa4F67he0zszWfPm3C6RsRSoc61TBr8dMyo9BTDAo9IiJSFIZhsGfPHmYtWMbynQnk1GiJi39N2+chfh608M3hyZ7taBwaYMdKKweFnmJQ6BERkX/KMAx27tzJgeQ8dqa5s2LnWdJzC2yfe+UkcUejqrzQ/1YCq3rbsdKKS6GnGBR6RETkWuXkW5i88DumrY3DCGqIyXzpNnejII/quWfoFxHGk1G3UcXT42/2JEWl0FMMCj0iIlJSDMPgm5+2MGXFVnameWDy+2P4y8fFYEBkXfpGhNIgSFd/rpVCTzEo9IiISGmwWq18uXYTM7/bwzGLP7hWsX0W6JyNb8oBHureml7du+Lq6mrHSssnhZ5iUOgREZHSlldg5bsDSSzeHs/3B5IosF46DRuWfApOxhHhX8AjvW7itltvwcXFxc7Vlg8KPcWg0CMiImUpOT2HyUt/ZvnuJNKcfG3tlsyLWI5voUfjakyf+Iqe/vw3FHqKQaFHRETsZd/pi/zn6y18fyyDPPMfi5w2Cvamb+tQerWsyaFd24iMjMTJSWuA/ZlCTzEo9IiIiL0VWKx8fyCRT9bGEZdspcB6qd1sgozDW3GJj6VX27oMiOpHx44dMZvN9i3YASj0FINCj4iIOJLUrHy+3nWGxdvj2XHyoq3dkpNB1v4f8UjYyT2dI+gfFUVkZGSlDUAKPcWg0CMiIo7q6LkMYrb+RszWE1zI/aM9//wpMvZ8x6xXHqbvHbfar0A7UugpBoUeERFxdBarweaj54n59TdW70kg/7/DXyYTdKwbQN+IEPas+Zy0lGSioqJo06ZNhZ8IrdBTDAo9IiJSnmTkFrBq91kWbTvF1hMXbO1Gfg6Z+38iY/c6arhk0z/qXqKiomjVqlWFDEAKPcWg0CMiIuXVyfNZLNkRz5Lt8ZxMyba1519MIHPPd2TuWU94dW8effRRRo0aZcdKS55CTzEo9IiISHlnGAa/nrjA4th4Vuw6Q2aexfZZzqk9dAiCeW+OwtvdBYvFwv79+7n++uvL9RUghZ5iUOgREZGKJDvPwtq9CSzeHs/Ph5P5/YTv7mKm+/XB1HM6z5P9u9G4UUOioqLo378/jRs3tmvNxaHQUwwKPSIiUlGdTc1m6Y7TLI6N5+i5TFu7JT2ZjD3fk7FnPQUp8TRt2tQWgBo0aGDHiotOoacYFHpERKSiMwyDuFMXWbw9nq93niU1O9/2Wd7Zg6TvWkfW/h+x5mYSFxdHixYt7Fht0Sj0FINCj4iIVCa5BRbW709icWw8Gw6dw/LfxU9NVgums7v5dPQIOjcMxNnJzJtvvonZbCYqKorrrrvOzpUXptBTDAo9IiJSWZ1Lz+WruNMsio3nQEK6rT3Ay5WezYP5cNQgzh/dBUCbNm3o378/9957L7Vr17ZXyTYKPcWg0CMiIgJ7z6SyOPY0X8Wd5nxmnq3dNSuJxF+Wk7H3e6xZqQC0b9+eJ554gkGDBl22n5MnT5KcnHzV4wQEBFCrVq1rrreo52/naz6SiIiIVCjX1/Tl+pq+jO7RiA0Hz7E4Np71BxLJ8wyk6i0jqXbLCNwvHOW3DTFs2baV2w4etG2bm5tLcnIyFouFhg0bkpOTc9XjuLu7c/DgwRIJPkWh0CMiIiJX5OJk5rYmQdzWJIgLmXmXFj+NjWdnfCpZVetRvc8Y3M0WUuv5s/PURZqH+rJmzRr69OlDy5Yt/wg8JjNuodfj5FUVS8YFcuP3gmElJyeH5ORkhR4RERFxHFWruDIkMpwhkeEcTkxn8fbTLN0RT2JaLsv3p7J8/0bqBXrhm3IWc5Wq7NixAwCPBpH43/oQzj7VbfsqSDtHyvpPyT60uUy/g+b0/Inm9IiIiBSdxWrw85FkFsfGs3ZvArkFl1Y/NZuguuU8h2J/xKdNb4BCT3w2DCtg4tyyN/n5y49o3br1NdWhOT0iIiJSqpzMJjo1qE6nBtVJy8ln1a6zLIqNZ9tvF0g0V8O3bR8Mw7hsiQuTyYxhWPG/9SHbbfJlwVxmRxIREZEKy8fdhfva1WLRozew4fnOdKrtDnDVNb1MJjPOPtXZn5x3xc9Lg0KPiIiIlKjwgCq0ruFepL4XcqylXM0fFHpERESkxFV1L1rEKGq/kqDQIyIiIiXupkY1saQn/3fS8uUMw4olPZmbGtUss5oUekRERKTE1QmvzZv9WmMyXTlqmExm3uzXmjrhZbeMhUKPiIiIlIrBXZoxbVBravgWnt9Tw9edaYNaM7hLszKtR7esi4iISKnp3rQGtzUJZuvxFJLScwj0dqddHX+czFe+q6s0KfSIiIhIqXIym4isW83eZWh4S0RERCoHhR4RERGpFBR6REREpFJQ6BEREZFKQaFHREREKgWFHhEREakUFHpERESkUlDoERERkUpBoUdEREQqBYUeERERqRQUekRERKRSUOgRERGRSkGhR0RERCoFhR4RERGpFBR6REREpFJQ6BEREZFKQaFHREREKgWHDj1TpkwhPDwcd3d32rdvz9atW/+y/8KFC2nUqBHu7u40a9aMVatWlVGlIiIi4ugcNvQsWLCAUaNGMW7cOLZv306LFi3o1q0bSUlJV+y/adMmBgwYwIgRI9ixYwe9e/emd+/e7Nmzp4wrFxEREUdkMgzDsHcRV9K+fXvatm3LRx99BIDVaiUsLIwnn3ySl1566bL+/fv3JzMzkxUrVtjaOnToQMuWLZk2bdoVj5Gbm0tubq7tfWpqKrVq1eLUqVP4+PiU8DcSERGR0pCWlkZYWBgXL17E19f3qv2cy7CmIsvLyyM2NpbRo0fb2sxmM127dmXz5s1X3Gbz5s2MGjWqUFu3bt1YtmzZVY8zYcIExo8ff1l7WFhY8QoXERERu0lPTy9/oSc5ORmLxUJQUFCh9qCgIA4cOHDFbRISEq7YPyEh4arHGT16dKGgZLVaSUlJoVq1arRr145ff/31Gr7FH35PoLqCJFfTtm3bEvvfW2VU0X9+5e37OVK99qqlLI9bmscqyX2X5rnQMAzS09OpWbPmX/ZzyNBTVtzc3HBzcyvU5ufnB4CTk1OJ/1J8fHwUeuSKSuN/b5VJRf/5lbfv50j12quWsjxuaR6rPJ0L/+oKz+8cciJzQEAATk5OJCYmFmpPTEwkODj4itsEBwf/o/5/5/HHHy/WdiLFof+9XZuK/vMrb9/Pkeq1Vy1ledzSPJYj/S5LgkNPZG7Xrh0ffvghcGnoqVatWjzxxBNXnciclZXF119/bWu74YYbaN68+VUnMpeVtLQ0fH19SU1NdZi/fkRERMqSI5wLHXZ4a9SoUQwdOpQ2bdrQrl07Jk+eTGZmJsOHDwdgyJAhhISEMGHCBACefvppOnXqxKRJk7jzzjuZP38+27Zt49NPP7Xn1wAuDaONGzfusqE0ERGRysIRzoUOe6UH4KOPPuKdd94hISGBli1b8sEHH9C+fXsAOnfuTHh4ONHR0bb+Cxcu5OWXX+bEiRPUr1+fiRMn0qNHDztVLyIiIo7EoUOPiIiISElxyInMIiIiIiVNoUdEREQqBYUeERERqRQUekRERKRSUOhxACtWrKBhw4bUr1+fGTNm2LscERGRMtenTx+qVq1Kv379Su0YunvLzgoKCmjSpAnff/89vr6+REREsGnTJqpVq2bv0kRERMrMhg0bSE9PZ86cOSxatKhUjqErPXa2detWrr/+ekJCQvDy8uKOO+7gm2++sXdZIiIiZapz5854e3uX6jEUeq7Rjz/+SM+ePalZsyYmk4lly5Zd1mfKlCmEh4fj7u5O+/bt2bp1q+2zM2fOEBISYnsfEhLC6dOny6J0ERGREnGt58KyotBzjTIzM2nRogVTpky54ucLFixg1KhRjBs3ju3bt9OiRQu6detGUlJSGVcqIiJSOsrLuVCh5xrdcccdvPHGG/Tp0+eKn7/33ns8+OCDDB8+nCZNmjBt2jQ8PT2ZNWsWADVr1ix0Zef06dPUrFmzTGoXEREpCdd6LiwrCj2lKC8vj9jYWLp27WprM5vNdO3alc2bNwPQrl079uzZw+nTp8nIyGD16tV069bNXiWLiIiUqKKcC8uKw66yXhEkJydjsVgICgoq1B4UFMSBAwcAcHZ2ZtKkSXTp0gWr1cq//vUv3bklIiIVRlHOhQBdu3Zl586dZGZmEhoaysKFC4mMjCzRWhR6HMDdd9/N3Xffbe8yRERE7GbdunWlfgwNb5WigIAAnJycSExMLNSemJhIcHCwnaoSEREpO450LlToKUWurq5ERESwfv16W5vVamX9+vUlfslORETEETnSuVDDW9coIyODI0eO2N4fP36cuLg4/P39qVWrFqNGjWLo0KG0adOGdu3aMXnyZDIzMxk+fLgdqxYRESk55eVcqGUortGGDRvo0qXLZe1Dhw4lOjoagI8++oh33nmHhIQEWrZsyQcffED79u3LuFIREZHSUV7OhQo9IiIiUiloTo+IiIhUCgo9IiIiUiko9IiIiEiloNAjIiIilYJCj4iIiFQKCj0iIiJSKSj0iIiISKWg0CMiIiKVgkKPiIiIVAoKPSIixfTpp58SFhaG2Wxm8uTJvPrqq7Rs2bLI2584cQKTyURcXNxV+2zYsAGTycTFixev2sdkMrFs2bIiH1ekslLoEamkTCbTX75effVVe5fo0NLS0njiiSd48cUXOX36NA899BDPP/98oZWkRcSxaJV1kUrq7Nmztv9esGABr7zyCgcPHrS1eXl52f7bMAwsFgvOzuXvn4zSqv3kyZPk5+dz5513UqNGDVv7n39uIuJYdKVHpJIKDg62vXx9fTGZTLb3Bw4cwNvbm9WrVxMREYGbmxs///wzR48epVevXgQFBeHl5UXbtm1Zt25dof2Gh4fz5ptv8sADD+Dt7U2tWrX49NNPbZ/n5eXxxBNPUKNGDdzd3alduzYTJkwA4P7776d///6F9pefn09AQABz584FwGq1MmHCBOrUqYOHhwctWrRg0aJFtv6/Dwf9b+07d+6kS5cueHt74+PjQ0REBNu2bbNt9/PPP3PTTTfh4eFBWFgYTz31FJmZmVf82UVHR9OsWTMArrvuOkwmEydOnLji8NaMGTNo3Lgx7u7uNGrUiI8//vgvfy+rVq2iQYMGeHh40KVLF06cOPGX/X+XnJxMnz598PT0pH79+ixfvrxI24lUKoaIVHqzZ882fH19be+///57AzCaN29ufPPNN8aRI0eM8+fPG3Fxcca0adOM3bt3G4cOHTJefvllw93d3fjtt99s29auXdvw9/c3pkyZYhw+fNiYMGGCYTabjQMHDhiGYRjvvPOOERYWZvz444/GiRMnjJ9++sn44osvDMMwjBUrVhgeHh5Genq6bX9ff/214eHhYaSlpRmGYRhvvPGG0ahRI2PNmjXG0aNHjdmzZxtubm7Ghg0b/rL266+/3hg0aJCxf/9+49ChQ0ZMTIwRFxdnGIZhHDlyxKhSpYrx/vvvG4cOHTI2btxotGrVyhg2bNgVf15ZWVnGunXrDMDYunWrcfbsWaOgoMAYN26c0aJFC1u/zz77zKhRo4axePFi49ixY8bixYsNf39/Izo62jAMwzh+/LgBGDt27DAMwzBOnjxpuLm5GaNGjTIOHDhgfPbZZ0ZQUJABGBcuXLjq7w8wQkNDjS+++MI4fPiw8dRTTxleXl7G+fPni/DbF6k8FHpE5KqhZ9myZX+77fXXX298+OGHtve1a9c2Bg0aZHtvtVqNwMBAY+rUqYZhGMaTTz5p3HLLLYbVar1sX/n5+UZAQIAxd+5cW9uAAQOM/v37G4ZhGDk5OYanp6exadOmQtuNGDHCGDBgwF/W7u3tbQsb/2vEiBHGQw89VKjtp59+Msxms5GdnX3FbXbs2GEAxvHjx21t/xt66tatawt0v3v99deNyMhIwzAuDz2jR482mjRpUqj/iy++WKTQ8/LLL9veZ2RkGICxevXqq24jUhlpeEtErqpNmzaF3mdkZPD888/TuHFj/Pz88PLyYv/+/Zw8ebJQv+bNm9v++/dhs6SkJACGDRtGXFwcDRs25KmnnuKbb76x9XV2diYqKorPP/8cgMzMTL766isGDhwIwJEjR8jKyuK2227Dy8vL9po7dy5Hjx79y9pHjRrFyJEj6dq1K2+99Vah/jt37iQ6OrrQPrt164bVauX48ePF+tllZmZy9OhRRowYUWi/b7zxxmW1/m7//v20b9++UFtkZGSRjvfnn3mVKlXw8fGx/cxF5JLyNytRRMpMlSpVCr1//vnn+fbbb3n33XepV68eHh4e9OvXj7y8vEL9XFxcCr03mUxYrVYAWrduzfHjx1m9ejXr1q0jKiqKrl272ublDBw4kE6dOpGUlMS3336Lh4cH3bt3By6FLoCVK1cSEhJS6Bhubm5/Wfurr77K/fffz8qVK1m9ejXjxo1j/vz59OnTh4yMDB5++GGeeuqpy34GtWrVKtLP6n/9Xuv06dMvCzJOTk7F2udf+aufuYhcotAjIkW2ceNGhg0bRp8+fYBLJ/aiTrT9Mx8fH/r370///v3p168f3bt3JyUlBX9/f2644QbCwsJYsGABq1ev5t5777Wd0Js0aYKbmxsnT56kU6dO//i4DRo0oEGDBjz77LMMGDCA2bNn06dPH1q3bs2+ffuoV6/eP97n1QQFBVGzZk2OHTtmu1L1dxo3bnzZBORffvmlxGoSqewUekSkyOrXr8+SJUvo2bMnJpOJsWPH/uOrCe+99x41atSgVatWmM1mFi5cSHBwMH5+frY+999/P9OmTePQoUN8//33tnZvb2+ef/55nn32WaxWKzfeeCOpqals3LgRHx8fhg4desVjZmdn88ILL9CvXz/q1KlDfHw8v/76K3379gXgxRdfpEOHDjzxxBOMHDmSKlWqsG/fPr799ls++uijf/6D+q/x48fz1FNP4evrS/fu3cnNzWXbtm1cuHCBUaNGXdb/kUceYdKkSbzwwguMHDmS2NhYoqOji318ESlMc3pEpMjee+89qlatyg033EDPnj3p1q0brVu3/kf78Pb2ZuLEibRp04a2bdty4sQJVq1ahdn8xz9HAwcOZN++fYSEhNCxY8dC27/++uuMHTuWCRMm0LhxY7p3787KlSupU6fOVY/p5OTE+fPnGTJkCA0aNCAqKoo77riD8ePHA5fmw/zwww8cOnSIm266iVatWvHKK69Qs2bNf/Td/tfIkSOZMWMGs2fPplmzZnTq1Ino6Oir1lqrVi0WL17MsmXLaNGiBdOmTePNN9+8phpE5A8mwzAMexchIiIiUtp0pUdEREQqBYUeERERqRQUekRERKRSUOgRERGRSkGhR0RERCoFhR4RERGpFBR6REREpFJQ6BEREZFKQaFHREREKgWFHhEREakUFHpERESkUvh/A9UEQmN0rxsAAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "from quantax.operator import sigma_z\n", "\n", "N = lattice.Nsites\n", "M = sum(sigma_z(i) for i in range(N)) / N\n", "M2 = M @ M\n", "\n", "h = [1.0, 3.0, 10.0]\n", "\n", "exactM2 = np.array([state0 @ M2 @ state0, state1 @ M2 @ state1, state2 @ M2 @ state2])\n", "sampledM2 = np.array([\n", " M2.expectation(state0, samples0),\n", " M2.expectation(state1, samples1),\n", " M2.expectation(state2, samples2),\n", "])\n", "\n", "plt.semilogx(h, exactM2, 'ks--', label='Exact')\n", "plt.semilogx(h, sampledM2, 'o-', label='Sampled')\n", "plt.xlabel('Transverse field h')\n", "plt.ylabel(r'$\\langle M^2 \\rangle$')\n", "plt.ylim(0, 1)\n", "plt.legend()\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "91007ced", "metadata": {}, "source": [ "One can also show snapshots from generated samples." ] }, { "cell_type": "code", "execution_count": 5, "id": "81222beb", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAs0AAAD3CAYAAADmMWljAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAPmElEQVR4nO3deWxWZb7A8V8paJGhIm63HSBuVZCoHRUJE1FxGYxMxBgwRtHBGeMCcSFq3EaBiHHBJQYj40JcR525RoVojGJExEh0iCXGSPQmxjBIb1QEVIos9tw/CO/QW9pfwZZC+XwSE3ve85w+J+nT833f9/SlrCiKIgAAgBZ16+wJAADAzk40AwBAQjQDAEBCNAMAQEI0AwBAQjQDAEBCNAMAQEI0AwBAQjQDAEBCNO8AU6ZMibKysvjuu+86eypAO7GuYddmDbOtRHMXMWfOnDj22GOjoqIiBgwYEJMnT46NGze2aWxjY2Pce++9cfDBB0dFRUUcffTR8cILL3TwjIHWTJo0KY499tjo27dv7LXXXjFo0KCYMmVK/PTTT20+xqxZs2LQoEFRUVERNTU1MWPGjA6cMbClf/zjHzFu3LioqamJsrKyOOWUU1rcd926dXHjjTdGdXV19OzZM4YOHRpz585t8/f6+uuv47zzzos+ffpEZWVljB49Or788st2OAu2JJq7gDfeeCPOOeec6NOnT8yYMSPOOeecmDZtWlx11VVtGn/rrbfGjTfeGGeccUbMmDEjBgwYEBdccEG8+OKLHTxzoCX/+te/Yvjw4TF16tR46KGHYsSIEXH33XfHmWeeGY2Njen4Rx99NC699NIYPHhwzJgxI4YNGxZXX3113HPPPTtg9sDMmTNj9uzZ0b9//9hnn31a3Xf8+PHxwAMPxIUXXhgPPfRQlJeXx1lnnRXvv/9++n1++umnGDFiRMyfPz9uueWWmDp1atTV1cXJJ58cK1asaK/TISKioMNNnjy5iIji22+/7ZDjH3nkkcUxxxxTbNiwobTt1ltvLcrKyoolS5a0OnbZsmVFjx49iokTJ5a2NTY2FsOHDy/69etXbNy4sUPmDLu6jl7XW3PfffcVEVEsXLiw1f0aGhqKfffdtxg1alST7RdeeGHRq1ev4vvvv+/IacIuoaPX8NKlS4tffvmlKIqiGDx4cHHyySdvdb8PP/ywiIhi+vTppW1r164tDj300GLYsGHp97nnnnuKiCg++uij0rYlS5YU5eXlxc033/zrToImvNK8A61atSrGjx8fffr0ib333jsuueSSaGho+FXH/Oyzz+Kzzz6Lyy67LLp3717aPmHChCiKIl566aVWx8+ePTs2bNgQEyZMKG0rKyuLK6+8MpYtWxYLFy78VfODrq4j1nVLDjrooNL3bM28efNixYoVTdZ1RMTEiRNjzZo18frrr3fI/GBX1FFruH///tGtW55ZL730UpSXl8dll11W2lZRURF/+ctfYuHChfHvf/87HT9kyJAYMmRIadvAgQPjtNNOi3/+85/bfwI00z3fhfZy3nnnxcEHHxx33XVXfPzxx/HEE0/EAQccUHq7dPXq1bFhw4b0OBUVFfGb3/wmIiLq6uoiIuL4449vsk91dXX069ev9HhL6urqolevXjFo0KAm20844YTS4yeeeGLbThB2Qx2xrjfbuHFjrFq1KtavXx+ffvpp/PWvf43evXuX1mdLWvq9cNxxx0W3bt2irq4uxo0bty2nCV1WR67htqirq4vDDz88Kisrm2zfvM4XL14c/fv33+rYxsbG+OSTT+LPf/5zs8dOOOGEeOutt+LHH3+M3r17b/O8aE4070C/+93vYtasWaWvV6xYEbNmzSotzNGjR8f8+fPT4/zpT3+Kp556KiIi6uvrIyKiqqqq2X5VVVWxfPnyVo9VX18fBx54YJSVlTUbGxHpeNjddcS63mzRokUxbNiw0tdHHHFEzJkzJ/r27dvqserr66O8vDwOOOCAJtv32GOP2Hfffa1r2EJHruG2qK+vb/EaHtH6dfj777+PdevWpeOPOOKIbZ4XzYnmHeiKK65o8vXw4cPjlVdeiR9++CEqKyvj/vvvj5UrV6bHqa6uLv3/2rVrIyJizz33bLZfRUVF/PDDD60ea+3atS2O3fL4wNZ1xLre7Mgjj4y5c+fGmjVr4oMPPoi33367TZ+esXbt2thjjz22+lhFRYV1DVvoyDXcFr/mOpw1QDaebSOad6ABAwY0+XrzX9OuXLkyKisr47jjjtvmY/bs2TMiNn1czf/3888/lx5vbXxLY7c8PrB1HbGuN6usrIzTTz89Ija92vX888/H6NGj4+OPP45jjjmmxXE9e/aM9evXb/WxtvxegN1JR67htvg11+GsAbLxbBvRvAOVl5dvdXtRFBGx6W2Wli50W+rZs2fsvffeEfGft1/q6+ub3fNUX1+f3vtYVVUV8+bNi6Iomtyisfm2j+195gy7i45Y1y0599xz46KLLooXX3yx1WiuqqqKX375Jb755psmt2isX78+VqxYYV3DFnbkGt6aqqqq+Prrr5ttb8t1uG/fvrHnnnuW9t3W8Wwbn56xEzn33HOjqqoq/e+aa64pjamtrY2ITfc+bmn58uWxbNmy0uMtqa2tjYaGhliyZEmT7R9++GGT4wPbZ3vWdUvWrVsXjY2NsXr16lb3a+n3wqJFi6KxsdG6hm3Qnmt4a2pra+OLL75odjtlW67D3bp1i6OOOqrZWt88/pBDDvFHgO3IK807ke25b2rw4MExcODAeOyxx+Lyyy8vPWOeOXNmlJWVxZgxY0r7rl69uvQHB5ufDY8ePTomTZoUjzzySDz88MMRsenZ9d/+9rf47W9/G7///e/b8xRht7M963rVqlXRq1ev6NGjR5N9nnjiiYho+qkYDQ0NsXTp0thvv/1iv/32i4iIU089Nfr27RszZ86Ms846q7TvzJkzY6+99opRo0b9qnOC3UlH39M8ZsyYuO++++Kxxx6L66+/PiI2PUF+8sknY+jQoU3eRV66dGk0NDTEwIEDm4y/6aabYtGiRaXfDZ9//nm88847pePRPkTzTmR775uaPn16nH322fGHP/whzj///Pj000/j4YcfjksvvbTJR8m98sorcckll8STTz4Z48ePj4iIfv36xbXXXhvTp0+PDRs2xJAhQ+LVV1+NBQsWxN///vcW37YC2mZ71vW7774bV199dYwZMyZqampi/fr1sWDBgnj55Zfj+OOPb/JxcR999FGMGDEiJk+eHFOmTImITW8T33HHHTFx4sQYO3ZsjBw5MhYsWBDPPfdc3HnnnemnbwD/sb3X5vfeey/ee++9iIj49ttvY82aNTFt2rSIiDjppJPipJNOioiIoUOHxtixY+Pmm2+Ob775Jg477LB4+umn46uvvmryqR4RERdffHHMnz+/dOtIxKZ/l+Hxxx+PUaNGxfXXXx89evSIBx54IA488MC47rrrtmvubJ1o7gL++Mc/xssvvxxTp06Nq666Kvbff/+45ZZb4vbbb2/T+Lvvvjv22WefePTRR+Opp56KmpqaeO655+KCCy7o4JkDW3PUUUfFiBEjYvbs2VFfXx9FUcShhx4at99+e9xwww0tfjLGliZMmBA9evSI+++/P+bMmRP9+/ePBx98cLvfQga2zTvvvBNTp05tsu22226LiIjJkyeXojki4plnnonbbrstnn322Vi5cmUcffTR8dprrzXZpyW9e/eOd999NyZNmhTTpk2LxsbGOOWUU+LBBx+M/fffv31PajdXVmz5dAUAAGjGHwICAEBCNAMAQEI0AwBAQjQDAEBCNAMAQEI0AwBAQjQDAECizf+4yRndxnbkPKBLmtv43509hVY1/m9NZ0+BFoysru3sKdCCnX1dd6Xr9ZvLF3f2FNhNdPuv/8n32QHzAACAXZpoBgCAhGgGAICEaAYAgIRoBgCAhGgGAICEaAYAgIRoBgCAhGgGAICEaAYAgIRoBgCAhGgGAICEaAYAgIRoBgCAhGgGAICEaAYAgIRoBgCAhGgGAICEaAYAgIRoBgCAhGgGAICEaAYAgIRoBgCAhGgGAICEaAYAgIRoBgCAhGgGAICEaAYAgIRoBgCAhGgGAICEaAYAgIRoBgCAhGgGAICEaAYAgIRoBgCAhGgGAICEaAYAgIRoBgCAhGgGAICEaAYAgIRoBgCAhGgGAICEaAYAgIRoBgCAhGgGAICEaAYAgIRoBgCAhGgGAICEaAYAgET3zp4A0HlGVtd29hTazZvLF3f2FNpVVzqfrvRztivoSj877Ly62rqe25jv45VmAABIiGYAAEiIZgAASIhmAABIiGYAAEiIZgAASIhmAABIiGYAAEiIZgAASIhmAABIiGYAAEiIZgAASIhmAABIiGYAAEiIZgAASIhmAABIiGYAAEiIZgAASIhmAABIiGYAAEiIZgAASIhmAABIiGYAAEiIZgAASIhmAABIiGYAAEiIZgAASIhmAABIiGYAAEiIZgAASIhmAABIiGYAAEiIZgAASIhmAABIiGYAAEiIZgAASIhmAABIiGYAAEiIZgAASIhmAABIiGYAAEiIZgAASIhmAABIiGYAAEiIZgAASIhmAABIiGYAAEiIZgAASIhmAABIiGYAAEiIZgAASHTv7AkAtIeR1bWdPYV29ebyxZ09BXZRXW0tdCVdaV13pXNpK680AwBAQjQDAEBCNAMAQEI0AwBAQjQDAEBCNAMAQEI0AwBAQjQDAEBCNAMAQEI0AwBAQjQDAEBCNAMAQEI0AwBAQjQDAEBCNAMAQEI0AwBAQjQDAEBCNAMAQEI0AwBAQjQDAEBCNAMAQEI0AwBAQjQDAEBCNAMAQEI0AwBAQjQDAEBCNAMAQEI0AwBAQjQDAEBCNAMAQEI0AwBAQjQDAEBCNAMAQEI0AwBAQjQDAEBCNAMAQEI0AwBAQjQDAEBCNAMAQEI0AwBAQjQDAEBCNAMAQEI0AwBAQjQDAEBCNAMAQEI0AwBAQjQDAEBCNAMAQEI0AwBAQjQDAECie2dPAIDmRlbXdvYU2EW9uXxxZ0+h3XS1ddCVzqcr/Zy1lVeaAQAgIZoBACAhmgEAICGaAQAgIZoBACAhmgEAICGaAQAgIZoBACAhmgEAICGaAQAgIZoBACAhmgEAICGaAQAgIZoBACAhmgEAICGaAQAgIZoBACAhmgEAICGaAQAgIZoBACAhmgEAICGaAQAgIZoBACAhmgEAICGaAQAgIZoBACAhmgEAICGaAQAgIZoBACAhmgEAICGaAQAgIZoBACAhmgEAICGaAQAgIZoBACAhmgEAICGaAQAgIZoBACAhmgEAICGaAQAgIZoBACAhmgEAICGaAQAgIZoBACAhmgEAICGaAQAgIZoBACAhmgEAICGaAQAgIZoBACAhmgEAIFFWFEXR2ZMAAICdmVeaAQAgIZoBACAhmgEAICGaAQAgIZoBACAhmgEAICGaAQAgIZoBACAhmgEAIPF/a3Sgr3JsQEQAAAAASUVORK5CYII=", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "fig, axes = plt.subplots(1, 3, figsize=(9, 3))\n", "\n", "axes[0].imshow(samples0.spins[0].reshape(4, 4))\n", "axes[0].axis('off')\n", "axes[0].set_title(\"h=0.0\")\n", "\n", "axes[1].imshow(samples1.spins[0].reshape(4, 4))\n", "axes[1].axis('off')\n", "axes[1].set_title(\"h=3.0\")\n", "\n", "axes[2].imshow(samples2.spins[0].reshape(4, 4))\n", "axes[2].axis('off')\n", "axes[2].set_title(\"h=10.0\")\n", "\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "c6d4a81d", "metadata": {}, "source": [ "## Metropolis sampler\n", "\n", "`ExactSampler` in general cannot be used in large systems, in which case we need metropolis samplers. Below we show how to create a `LocalFlip` sampler, which proposes new spin configurations by flipping spins locally." ] }, { "cell_type": "code", "execution_count": 6, "id": "5f41b709", "metadata": {}, "outputs": [], "source": [ "sampler0 = qtx.sampler.LocalFlip(state0, nsamples=1024)\n", "sampler1 = qtx.sampler.LocalFlip(state1, nsamples=1024)\n", "sampler2 = qtx.sampler.LocalFlip(state2, nsamples=1024)" ] }, { "cell_type": "markdown", "id": "c341d5fd", "metadata": {}, "source": [ "You might have noticed that it takes some time to construct the Metropolis sampler. It's because the spins stored in the sampler are thermalized during initialization. You can manually set an argument `thermal_steps=0` to initialize without thermalization, but please remember to thermalize spin configurations in the future.\n", "\n", "The magnetization can be measured similarly." ] }, { "cell_type": "code", "execution_count": 7, "id": "d7c00ab6", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAj0AAAG6CAYAAADjzPf8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABb2ElEQVR4nO3dd3hUZcLG4d9MKqkQAmkkBKUYWiIICKiAoqBIU1ZEBMSuq6ioq7AiFlYURV0FRVAIqEhC700UpUkJBER6DQQSSID0OnO+P/gczFIMIclMkue+rrnWOfOeM88Elnly2msyDMNAREREpJIz2zuAiIiISHlQ6REREZEqQaVHREREqgSVHhEREakSVHpERESkSlDpERERkSpBpUdERESqBJUeERERqRJUekRERKRKUOkRERGRKsFhS8+vv/5K9+7dCQ4OxmQyMW/evL9dZ/Xq1bRo0QI3Nzfq169PdHR0mecUERGRisFhS09WVhaRkZGMHz++WOMPHz5Mt27d6NSpE/Hx8bz44os8/vjjLF++vIyTioiISEVgqggTjppMJubOnUuvXr0uO+a1115j8eLF7Ny507bswQcf5Ny5cyxbtqwcUoqIiIgjc7Z3gNKyYcMGOnfuXGRZly5dePHFFy+7Tl5eHnl5ebbnVquVM2fOULNmTUwmU1lFFRERkVJkGAYZGRkEBwdjNl/+IFalKT1JSUkEBAQUWRYQEEB6ejo5OTlUq1btonVGjx7N22+/XV4RRUREpAwdO3aMOnXqXPb1SlN6SmLYsGEMHTrU9jwtLY2wsDCOHTuGj4+PHZOJiIhIcaWnpxMaGoq3t/cVx1Wa0hMYGEhycnKRZcnJyfj4+FxyLw+Am5sbbm5uFy338fFR6REREalg/u7UFIe9eutqtW3bllWrVhVZtnLlStq2bWunRCIiIuJIHLb0ZGZmEh8fT3x8PHD+kvT4+HgSEhKA84emBg4caBv/9NNPc+jQIf71r3+xZ88evvjiC2JjY3nppZfsEV9EREQcjMMe3tqyZQudOnWyPf/z3JtBgwYRHR3NyZMnbQUIoF69eixevJiXXnqJ//73v9SpU4evv/6aLl26lHt2ERFxfBaLhYKCAnvHkGJwcXHBycnpmrdTIe7TU17S09Px9fUlLS1N5/SIiFRShmGQlJTEuXPn7B1FrkL16tUJDAy85Hk7xf3+dtg9PSIiImXhz8JTu3ZtPDw8dF82B2cYBtnZ2Zw6dQqAoKCgEm9LpUdERKoMi8ViKzw1a9a0dxwppj+vwj516hS1a9cu8aEuhz2RWUREpLT9eQ6Ph4eHnZPI1frzz+xazsNS6RERkSpHh7QqntL4M1PpERERkSpBpUdERESqBJ3ILCIiUkwJCQmkpKRc9nV/f3/CwsLKMZFcDZUeERGRYkhISKBRo0bk5uZedoy7uzt79+4tk+LzyCOPMHXq1IuWd+nShWXLlpX6+/2vt956i3nz5tlmSqiIVHpERESKISUl5YqFByA3N5eUlJQy29vTtWtXpkyZUmTZpSbOlkvTOT0iIlLlZWVlXfbxd0WnJNstKTc3NwIDA4s8atSowerVq3F1dWXNmjW2sWPGjKF27dokJycDsGzZMm655RaqV69OzZo1uffeezl48GCR7R8/fpx+/frh5+eHp6cnN910Exs3biQ6Opq3336b7du3YzKZMJlMREdHl/hz2Iv29JSRP4/7WqwGu1PyOZtrpYa7mQh/V5zMJh33FRFxIF5eXpd97Z577mHx4sUl2m54ePglzwEq7RmgOnbsyIsvvsiAAQPYvn07hw4dYsSIEcycOZOAgADgfAEbOnQozZs3JzMzkzfffJPevXsTHx+P2WwmMzOTDh06EBISwoIFCwgMDGTr1q1YrVb69u3Lzp07WbZsGT/++CMAvr6+pfoZyoNKTxn487ivKexG/O54EmefWrbXCtNPc2bVRIyEbWV23FdERCqnRYsWXVTQhg8fzvDhwxk1ahQrV67kySefZOfOnQwaNIgePXrYxt1///1F1ps8eTK1atVi165dNG3alOnTp3P69Gk2b96Mn58fAPXr17eN9/LywtnZmcDAwDL8hGVLpacMpKSkYAq7kVq9hl/0mpN3TWr1Gs7pee+V6XFfEREpvszMzMu+di2zex85cqTE615Kp06d+PLLL4ss+7OguLq68v3339O8eXPq1q3LJ598UmTc/v37efPNN9m4cSMpKSlYrVbg/C/qTZs2JT4+nhtvvNG2vcpIpacMWKwGfnc8CVx8B0mTyYxhWPG740ksVk1wLyLiCDw9PSvEdj09PYvsfflf69evB+DMmTOcOXOmyPt3796dunXrMmnSJIKDg7FarTRt2pT8/HzgwvxWlZlOZC4Du1PycfapddlbZptMZpx9arE7Jb+ck4mISGV18OBBXnrpJSZNmkSbNm0YNGiQbW9Oamoqe/fu5Y033uCOO+4gIiKCs2fPFlm/efPmxMfHc+bMmUtu39XVFYvFUuafoyyp9JSBs7nWUh0nIiL25+/vj7u7+xXHuLu74+/vX2YZ8vLySEpKKvJISUnBYrHw8MMP06VLFwYPHsyUKVPYsWMHY8eOBaBGjRrUrFmTiRMncuDAAX766SeGDh1aZNv9+vUjMDCQXr16sW7dOg4dOsTs2bPZsGEDcP6k7MOHDxMfH09KSgp5eXll9jnLig5vlYEa7sXrksUdJyIi9hcWFsbevXvtekfmZcuWERQUVGRZo0aNeOihhzh69CiLFi0CICgoiIkTJ9KvXz/uuusuIiMjmTFjBkOGDKFp06Y0atSIzz77jI4dO9q24+rqyooVK3j55Ze55557KCwspHHjxowfPx44fyL0nDlz6NSpE+fOnWPKlCk88sgjZfZZy4LJKO3r5iqw9PR0fH19SUtLw8fHp8Tb2bwljt6Tf8fJuyYm08XFxjAMrDkZzBnYkFY3tbyWyCIichVyc3M5fPgw9erV+9u9NuJYrvRnV9zvb+1qKANOZhNnVk0ETBhG0UNYhmFgMpkwu3kwe+3v9gkoIiJSBan0lAF/f3+MhG2cnvcelozUIq9ZMlLIPb4Lk5MzPyR4s2zTbjulFBERqVp0Tk8Z+Otx34vvyBxITm4ej38bR36N6xi+/DgN6oVyfa3L3w1URERErp1KTxkJCwuznczW6hKv/xYZxaDorfyemMbAbzYx+5l2BPrq+LKIiEhZ0eEtO/Hz9mDK4FbU8/ck8VwOXUYvIPH02b9fUUREREpEpceO/L3cmPZoa5wLs0jDk87vzCI9q+Sz+YqIiMjlqfTYWaifB+/dGYw1N5Mcz2A6/vtb8gsr9h0vRUREHJFKjwN4oMstvNzaC6MgjzPuwdz1xlR0+yQREZHSpdLjIF54qBv9wnMwrBaOEECfUd/bO5KIiEilotLjQN5/vj+dPE8AEJdVg2c/m23nRCIicikWq8GGg6nMj09kw8FULNaqvXfeZDIxb968a9rGI488Qq9evUolz+Wo9DiYKW8+RUT+XgCWnHBndtxxOycSEZG/WrbzJLd88BP9Jv3GCzPi6TfpN2754CeW7TxZpu97+vRpnnnmGcLCwnBzcyMwMJAuXbqwbt26Mn3fykSlx8GYTCYWjhlC38jzs/T+a/YOftqTbOdUIiIC5wvPM99t5WRa0Sttk9Jyeea7rWVafO6//362bdvG1KlT2bdvHwsWLKBjx46kpqb+/coCqPQ4JGdnZ0b3bc19N4ZgsRo8/W0cC9Zrni4RkdJmGAbZ+YXFemTkFjBywR9c6kDWn8veWrCLjNyCYm3vai5YOXfuHGvWrOGDDz6gU6dO1K1bl9atWzNs2DB69OgBwMcff0yzZs3w9PQkNDSUZ599lszMTNs2oqOjqV69OosWLaJRo0Z4eHjQp08fsrOzmTp1KuHh4dSoUYMhQ4ZgsVy4ijg8PJx3332Xfv364enpSUhIiG3m9cs5duwYDzzwANWrV8fPz4+ePXty5MgR2+sWi4WhQ4dSvXp1atasyb/+9a9yuYBHd2R2UGaziQ/6NOdgYjLbTxUyZNYealX3pm3jcHtHExGpNHIKLDR+c3mpbMsAktJzafbWimKN3/VOFzxci/c17OXlhZeXF/PmzePmm2/Gzc3tojFms5nPPvuMevXqcejQIZ599ln+9a9/8cUXX9jGZGdn89lnnzFjxgwyMjK477776N27N9WrV2fJkiUcOnSI+++/n/bt29O3b1/beh9++CHDhw/n7bffZvny5bzwwgs0bNiQO++886IcBQUFdOnShbZt27JmzRqcnZ0ZNWoUXbt2ZceOHbi6ujJ27Fiio6OZPHkyERERjB07lrlz53L77bcX6+dRUtrT48BcnMx82qcppBwCVw/6T9zA3uOn7R1LRETKmbOzM9HR0UydOpXq1avTvn17hg8fzo4dO2xjXnzxRTp16kR4eDi33347o0aNIjY2tsh2CgoK+PLLL7nxxhu57bbb6NOnD2vXruWbb76hcePG3HvvvXTq1Imff/65yHrt27fn9ddfp2HDhjz//PP06dOHTz755JJZY2JisFqtfP311zRr1oyIiAimTJlCQkICq1evBuDTTz9l2LBh3HfffURERDBhwgR8fX1L94d2CdrT4+DqhYUw64XO3DduDdQI4d6PlrNmZC8Ca2iCUhGRa1XNxYld73Qp1thNh8/wyJTNfzsuenArWtfzK9Z7X43777+fbt26sWbNGn777TeWLl3KmDFj+Prrr3nkkUf48ccfGT16NHv27CE9PZ3CwkJyc3PJzs7Gw8MDAA8PD66//nrbNgMCAggPD8fLy6vIslOnThV577Zt2170/NNPP71kzu3bt3PgwAG8vb2LLM/NzeXgwYOkpaVx8uRJ2rRpY3vN2dmZm266qcwPcWlPTwVwU7MIvurXDEtGCgXuNej87hwycvLtHUtEpMIzmUx4uDoX63Frg1oE+bpjuty2gCBfd25tUKtY2zOZLrely3N3d+fOO+9kxIgRrF+/nkceeYSRI0dy5MgR7r33Xpo3b87s2bOJi4uznXeTn3/h+8LFxeWiz3+pZVar9aqz/SkzM5OWLVsSHx9f5LFv3z4eeuihEm+3NKj0VBB3d2zLO7fXwpKTTqZrTe4YGUN+Ycn/UoqIyNVxMpsY2b0xwEXF58/nI7s3xsl89WWmpBo3bkxWVhZxcXFYrVbGjh3LzTffTMOGDTlx4kSpvc9vv/120fOIiIhLjm3RogX79++ndu3a1K9fv8jD19cXX19fgoKC2Lhxo22dwsJC4uLiSi3v5aj0VCCD77+HZxobWPNzOGX246WYrVir+A2xRETKU9emQXz5cAsCfd2LLA/0defLh1vQtWlQmbxvamoqt99+O9999x07duzg8OHDzJw5kzFjxtCzZ0/q169PQUEBn3/+OYcOHeLbb79lwoQJpfb+69atY8yYMezbt4/x48czc+ZMXnjhhUuO7d+/P/7+/vTs2ZM1a9Zw+PBhVq9ezZAhQzh+/Py951544QXef/995s2bx549e3j22Wc5d+5cqeW9HJ3TU8EMf+ohXL9bzFe7TSz+PRl/rz94q0eTEu0mFRGRq9e1aRB3Ng5k0+EznMrIpba3O63r+ZXpHh4vLy/atGnDJ598wsGDBykoKCA0NJQnnniC4cOHU61aNT7++GM++OADhg0bxm233cbo0aMZOHBgqbz/yy+/zJYtW3j77bfx8fHh448/pkuXS58L5eHhwa+//sprr73GfffdR0ZGBiEhIdxxxx34+PjYtnfy5EkGDRqE2Wzm0UcfpXfv3qSlpZVK3ssxGZrZ0iY9PR1fX1/S0tJsfzCOan58Ii/MiAdgcEt/Rv6jzZVXEBERcnNzOXz4MPXq1cPd3f3vVxDCw8N58cUXefHFF+2a40p/dsX9/tbhrQqqZ1QII7rdAMCUuBRGz1xj50QiIiKOTaWnAhvQJpQaJ86fCDZh8zkmr9xm50QiIiKOS+f0VGCurq6s/vxlWv/zE/LqtOKdFQkE1PCi200N7B1NREQqib9OH1HRaU9PBefr68tPY56C49vByZnnZuxg0/6ynelXRESkIlLpqQRCgoNY/O/7sZzcg+HsTr8Ja9h/8py9Y4mIOCxdw1PxlMafmUpPJdEkohHfPX0rhaePYHHx5KGJ6zmVnmvvWCIiDuXPuw9nZ2fbOYlcrT//zP73DtJXQ+f0VCId2rbm49Nn+CCukNM5MGjKZmKeuhkf95L/BRERqUycnJyoXr26bW4pDw8P3efMwRmGQXZ2NqdOnaJ69eo4OV3dnGV/pfv0/EVFuk/PlRxNzeL+LzeQkplHq7rV+fbxm3G/yontREQqK8MwSEpKKpc7AEvpqV69OoGBgZcsqcX9/lbp+YvKUnoAdiam0fer9WTlW4nwzmfh6z1xdtLRTBGRP1ksFgoKCuwdQ4rBxcXlint4ivv9rcNblVTTEF9uM+1mSWE9dme4MmjcMr4bcrd244qI/D8nJ6drOlQiFY9+9a/Exo98kYhzGzGsFtadNBgavdrekUREROxGpacSM5vNLPjiHQKO/QzA3L3ZfDB3k51TiYiI2IdKTyXn6urKqonv4HFgFQBfbjzN1NW77JxKRESk/Kn0VAHe3t78PP41TPt/AWDk0oOs3nvKzqlERETKl0pPFREQEMCKD57E99w+MJl55rutbEs4a+9YIiIi5UalpwppUL8+m8e9wK0N/MkpsPBo9GYOnMqwdywREZFyodJTxbg6m5nwcEsiQ6tzNruA7mNXkHhWt2MXEZHKT6WnCvJ0c2ZYO18KUo+TY3Ln3g+XcjYr396xREREypRKTxV1841NebJRPoUZKZy1utP9oyVk5xfaO5aIiEiZUempwka+8jzdPA9jycngeI4L//h0OQUWq71jiYiIlAmVniruy/ffpGXmRqwFufxxBh6d8BNWq6ZjExGRykelp4ozm83MnDCGsIQVGFYLa47lMXzmFjQPrYiIVDYqPYKLiwvLpozFZ88CAGZsO8WEXw7ZOZWIiEjpUukRALy8vNgyczxvdIsA4INle4jdfMzOqUREREqPQ5ee8ePHEx4ejru7O23atGHTpitPlvnpp5/SqFEjqlWrRmhoKC+99BK5ubnllLbic3Nz4/Fbr+OpDtcB8Nrs7SzZftzOqUREREqHw5aemJgYhg4dysiRI9m6dSuRkZF06dKFU6cuPWfU9OnTef311xk5ciS7d+/mm2++ISYmhuHDh5dz8orvtS6N8En9AwMT//w+jt8OnrZ3JBERkWvmsKXn448/5oknnmDw4ME0btyYCRMm4OHhweTJky85fv369bRv356HHnqI8PBw7rrrLvr163fFvUN5eXmkp6cXecj5k5vH9m1J7sFNGGZnHp64jl0n0uwdS0RE5Jo4ZOnJz88nLi6Ozp0725aZzWY6d+7Mhg0bLrlOu3btiIuLs5WcQ4cOsWTJEu65557Lvs/o0aPx9fW1PUJDQ0v3g1Rgd3a+gw/va0zusT8oNLnQ5/OfOXZG01WIiEjF5ZClJyUlBYvFQkBAQJHlAQEBJCUlXXKdhx56iHfeeYdbbrkFFxcXrr/+ejp27HjFw1vDhg0jLS3N9jh2TCfu/tWAfn0Z2sqD/NNHyDZc6PnJSlIy8+wdS0REpEQcsvSUxOrVq3nvvff44osv2Lp1K3PmzGHx4sW8++67l13Hzc0NHx+fIg8p6l8vPsc/aiVTmJbMmQJnen+ykozcAnvHEhERuWoOWXr8/f1xcnIiOTm5yPLk5GQCAwMvuc6IESMYMGAAjz/+OM2aNaN379689957jB49GqtVUytci09Hv027/K0YORkcyzLx1Ldx5BVa7B1LRETkqjhk6XF1daVly5asWrXKtsxqtbJq1Sratm17yXWys7Mxm4t+HCcnJwDdXfgamUwmvp/wCRMfjsLT1Yn1B1N5KSYei6arEBGRCsQhSw/A0KFDmTRpElOnTmX37t0888wzZGVlMXjwYAAGDhzIsGHDbOO7d+/Ol19+yYwZMzh8+DArV65kxIgRdO/e3VZ+pORcXFzo0uoGvhpwEy5OJpb8nsQrP2xUoRQRkQrD2d4BLqdv376cPn2aN998k6SkJKKioli2bJnt5OaEhIQie3beeOMNTCYTb7zxBomJidSqVYvu3bvzn//8x14foVK6pYE/jzVx5svtecz9PZWAxb/z+r3N7R1LRETkb5kM/apuk56ejq+vL2lpaTqp+QqOHDlC+0eG4XLzwwC8eU8jHr2tvp1TiYhIVVXc72+HPbwljis8PJwl/32dnM2zAXhn8R4WxGu6ChERcWwqPVIikZGRzHhjAFnbl4HJxAs/bGPNfk1XISIijkulR0qsU6dOjBt8G9l71mKYzDz6zQZ2HD9n71giIiKXpNIj16TvA/9g+B0h5BzZTgFOPDJlM4dOZ9o7loiIyEVUeuSavfj8c/z3H41pEuzDmax8BnyzieT0XHvHEhERKUKlR0rFg/f3YuqjrQmv6UHiuRwenLCWtGxNVyEiIo5DpUdKjb+XG1MGtcTVmsvhM3k8NOFXcgs0XYWIiDgGlR4pVQFezvhsnYY1N5M/TuXy2DfrKbRo7jMREbE/lR4pVR4eHiyb8Q3um6OxFuSx7kg6L8+I03QVIiJidyo9Uur8/f1ZOX0C1rVfY1gtzP/9FP9ZtNPesUREpIpT6ZEyUbduXZZ+/QE5v3wDwNfrEpj4ywE7pxIRkapMpUfKTLNmzZg5ZigZa78F4L2le5mzVdNViIiIfaj0SJnq0KEDE196gDtCzj9/ddYOft5zyr6hRESkSlLpkTJ3//33Memf99D7xhAsVoNnvo8j7uhZe8cSEZEqRqVHyoXZbGJMn+a0quNBboGVARPXsS85w96xRESkClHpkXLj4mQm5MhS8hL3kG0x0ffLNSSey7F3LBERqSJUeqRcffT+e9yY+Rv5KQmczTXo+8WvnMnKt3csERGpAlR6pFw5Ozsz6/uphB6cR2H6KY6nF9L/q7Vk5RXaO5qIiFRyKj1S7jw8PFg6+we8tkzFkpPO7lM5PDblN/ILNV2FiIiUHZUesQs/Pz9WzJqGsXo81vxcfjuSxiszt2O1aroKEREpGyo9YjdhYWEs+34CdY+vwMlsYsH2E7yzaJfm6RIRkTKh0iN21bRpU9bEfsXHD0QCEL3+CF+sPmjnVCIiUhmp9IhD6BkVwpv3Ngbgw+V7+WFTgp0TiYhIZaPSIw4jstoZ0tbHADB8zg6W7UyycyIREalMVHrEYbRs2ZIhna4jY/tyDEw89/0WfjuUau9YIiJSSaj0iEN5662R3BeWT/a+DRQaJgZP/o0/TqTZO5aIiFQCKj3iUEwmE1+OH8fNxm5yE34npxD6T1zP0dQse0cTEZEKTqVHHI6zszMx07/jusSV5Ccf4lyulf6TNnAqI9fe0UREpAJT6RGHVK1aNRbNnUn1Hd/jY87j+Lk8Hpm8mfTcAntHExGRCkqlRxyWn58fcWt/ZuHLXfD3cmXXyXSenLaF3AKLvaOJiEgFpNIjDq1atWrUrelJ9ODWeLo68duhM7zww1Ysmq5CRESukkqPVAhNgn3w3PodRmEBy3ed4o15v2u6ChERuSoqPVIhmEwm/v3EP0hd9BGG1cIPm47x8cp99o4lIiIViEqPVBg9e/bkk5cHcWbFFwB8/tMBotcdtnMqERGpKFR6pEJ54oknePW+dpz79VsA3lq4i/nxiXZOJSIiFYFKj1Q4I0aM4MHIGqTHLQRgaEw8v+47bedUIiLi6FR6pMIxmUyMHzeOjj6nsRzeiMWAp7+LI/7YOXtHExERB6bSIxWSk5MT07//ntXvP8qtDfzJzrcweMomDpzKtHc0ERFxUCo9UmG5u7tzfb1wvny4JZF1fDmbXcDDkzZwMi3H3tFERMQBqfRIhefl5kz/0HQKzyaSlJHPw5N+41x2vr1jiYiIg1HpkUqhWcN65C77iMKMVA6mZDN4yiZy8jVdhYiIXKDSI5VCREQEC2dMJW3+f7DkZrLtWBrPfh9HgcVq72giIuIgVHqk0mjbti3Tv/yI1DnvYi3I5ee9p3lt1g6smqdLRERQ6ZFKpnv37nw+8iVS5r2PYbUwZ1sio5futncsERFxACo9Uuk89thj/PvRXqQu+S8Ak9Yc5qtfDto5lYiI2JuzvQOIlIXhw4fToEEDztRuxPvL9jJ66R78PF35x02h9o4mIiJ2otIjlZLJZOKBBx4A4Gx2AV/9eojXZ++ghocrnRsH2DmdiIjYgw5vSaX3cufrCcg+jMWAZ7+PY/ORM/aOJCIidqDSI5Vebm4uOb9+Q/aBTeRbDB6dsok9Sen2jiUiIuVMpUcqPR8fH5YtWYzHtunkHv+DjDwLA77eyLEz2faOJiIi5UilR6qE4OBgli9ZROFP48g/fYTTmfkM+GYjKZl59o4mIiLlRKVHqoxGjRqxeO5MMhaOpjAtmSOp56eryMwrtHc0EREpByo9UqW0adOGGVO+ImXWW1iy0/g9MZ2nvt1CXqHm6RIRqexUeqTK6datG5M+/g8jO/rj6erEugOpDI3djkXTVYiIVGq6T49USQMGDADghmYpDI7exOIdJ6np6crbPZpgMpnsnE5ERMqC9vRIlXZLA39evaU2GFambTjKZ6sO2DuSiIiUEZUeqfK2zv+a1JVfAfDJj/v47rejdk4kIiJlQaVHqrz//ve/tK6Rw7l10wEYMW8nS34/aedUIiJS2lR6pMpzc3Nj7ty51M34g4xtSzGAF2ZsY/2BFHtHExGRUqTSI8L/37V56VK89y8ha+86CiwGj0/bws7ENHtHExGRUqLSI/L/goKCWLFsGayfQu7R7WTnWxg0eROHU7LsHU1EREqBSo/IXzRs2JAlixbQLG0DNwR4kpp1frqKU+m59o4mIiLXSKVH5H+0atWKlUsW8u3jbalb04PjZ3MYOHkTaTkF9o4mIiLXwKFLz/jx4wkPD8fd3Z02bdqwadOmK44/d+4c//znPwkKCsLNze38b+1LlpRTWqlsanm78e2jbfA0F7InKYMnpm4ht0DTVYiIVFQOW3piYmIYOnQoI0eOZOvWrURGRtKlSxdOnTp1yfH5+fnceeedHDlyhFmzZrF3714mTZpESEhIOSeXyuTorjj2f/0S1rwsNh05w/M/bKPQYrV3LBERKQGTYRgOOeFQmzZtaNWqFePGjQPAarUSGhrK888/z+uvv37R+AkTJvDhhx+yZ88eXFxcivUeeXl55OXl2Z6np6cTGhpKWloaPj4+pfNBpEIzDIOhQ4fy5ayVBPR9F5OzKw/cVIcP7m+u6SpERBxEeno6vr6+f/v97ZB7evLz84mLi6Nz5862ZWazmc6dO7Nhw4ZLrrNgwQLatm3LP//5TwICAmjatCnvvfceFsvlD0eMHj0aX19f2yM0NLTUP4tUbCaTibFjx9L7lmacnv8BhtVK7JbjjFm+197RRETkKjlk6UlJScFisRAQEFBkeUBAAElJSZdc59ChQ8yaNQuLxcKSJUsYMWIEY8eOZdSoUZd9n2HDhpGWlmZ7HDt2rFQ/h1QOZrOZ6Oho2tX14szyzwH4cvVBvl5zyM7JRETkajhk6SkJq9VK7dq1mThxIi1btqRv3778+9//ZsKECZddx83NDR8fnyIPkUtxc3Njzpw5NHBK4ewvUwEYtXg3c7cdt3MyEREpLocsPf7+/jg5OZGcnFxkeXJyMoGBgZdcJygoiIYNG+Lk5GRbFhERQVJSEvn5+WWaV6oGHx8fli5div/prURVOwPAqzN38PPeS59cLyIijsUhS4+rqystW7Zk1apVtmVWq5VVq1bRtm3bS67Tvn17Dhw4gNV64cqaffv2ERQUhKura5lnlqohICCA33fsYM6Ih+kVFUyh1eDZ77ayNeGsvaOJiMjfcMjSAzB06FAmTZrE1KlT2b17N8888wxZWVkMHjwYgIEDBzJs2DDb+GeeeYYzZ87wwgsvsG/fPhYvXsx7773HP//5T3t9BKmkPDw8MJtNjOkTSfvrapBTYOHRKZvZn5xh72giInIFzvYOcDl9+/bl9OnTvPnmmyQlJREVFcWyZctsJzcnJCRgNl/obKGhoSxfvpyXXnqJ5s2bExISwgsvvMBrr71mr48glZyTyWDf5FfJa3gf54JvYODkTcx+ph3B1avZO5qIiFyCw96nxx6Ke52/yJ++//57Bj7+DAH9P8DVP4zra3ky6+l21PDUIVURkfJSoe/TI1JR9O/fnzGjRnIq9k0K009z8HQWg6M3k51faO9oIiLyP1R6RK7Ryy+/zItPDiI5dgTWnHTij53j6e+2kl+o6SpERByJSo9IKRgzZgx9u3YgedbbGAW5/LrvNK/O2o7VqqPHIiKOQqVHpBSYzWYmT55Mx6Z1sf46EScTzI8/wbuLd6HT5kREHIPDXr0lUtG4uroya9Yszp07R1yKmRdj4pmy7gj+Xm78s1N9e8cTEanyVHpESpG3tzfe3t6EhkJqVj7vLtrFh8v3UtPTlQdbh9k7nohIlabDWyJlxD91B5mbZgMwfO7vLP/j0pPliohI+VDpESkjoaGh5G2eScb25VgNeP6HbWw8lGrvWCIiVZZKj0gZadGiBXPnziXjp4lk79tAfqGVx6duYdeJdHtHExGpklR6RMpQ586dmRo9hZSFH5J7bCcZeYUMmrKJhNRse0cTEalyVHpEyli/fv0YO+Z9Ts1+l/xThzmdkceAyRs5nZFn72giIlWKSo9IOXjppZd4ZciznIp9E/fCTI6mZvPIlE1k5BbYO5qISJWh0iNSTt5//32mfvU5i/91D/5ervxxIp0np8WRW2CxdzQRkSpBpUeknJjNZvr378/1tb2JHtwaLzcnNhxK5aWYeCyarkJEpMyp9IjYQQN/d8ISloGlgKU7kxgxf6emqxARKWMqPSJ2kJGRwf61izm98CMwrEzfmMAnP+63dywRkUpNpUfEDvz9/Vm+fDne5w6QuuJLAD5btZ+p64/YN5iISCVWqqXnzJkzWK3W0tykSKV13XXXsXTpUkwH13JuzXcAvLXwDxZuP2HnZCIildM1l55du3bx/vvv065dO2rVqkXt2rUZOHAgs2fPJisrqzQyilRaN954I/PmzSN782zS4xZhGDA0Np41+0/bO5qISKVTotKzd+9eXn75ZRo0aMDNN9/M5s2befrpp0lOTmbJkiXUrVuXd955B39/f+6++26+/PLL0s4tUmncfvvtfPvtt5z7aRI5e9dRYDF46ts4th87Z+9oIiKViskowSUjU6ZMYePGjfTs2ZM77rgDV1fXS447cuQI8+fPZ9GiRaxcufKaw5a19PR0fH19SUtLw8fHx95xpIqZMGEC9RveQPShaqw9kIKfpyszn27L9bW87B1NRMShFff7u0Slp7JS6RFHkJlXyEOTfmPH8TRCqldj9jPtCPR1t3csERGHVdzv76s+vJWTk0NiYuJFy//444+r3ZSIXIKXmzOv3ewNGadIPJfDoMmbSMvWdBUiItfqqkrPrFmzaNCgAd26daN58+Zs3LjR9tqAAQNKPZxIVTX+4w84/t3rWLPOsDc5g8embiYnX9NViIhci6sqPaNGjSIuLo74+HimTJnCY489xvTp0wF0N1mRUjRp0iRa3hBO0owRGPnZbDl6ln9O30qBRbeEEBEpqasqPQUFBQQEBADQsmVLfv31V7766iveeecdTCZTmQQUqYo8PT1ZvHgx9fzcSI4dCZYCftpzitdn/65fMERESuiqSk/t2rXZsWOH7bmfnx8rV65k9+7dRZaLyLX7867NNSxnOTX3PTCszN56nPeX7rF3NBGRCumqSs+3335L7dq1iyxzdXXlhx9+4JdffinVYCIC4eHhLFu2DJfTe0lZ8ikAX/16iIm/HrRvMBGRCuiqSk+dOnUIDAy85Gvt27cvlUAiUlRkZCTz58/n1jquvNQpHID3luxhVtxx+wYTEalgnEu64qhRo2jRogUtW7a0necjImWjY8eOdOjQAZPJRJbFzMRfD/Ha7B3U8HDhjgj9/09EpDhKfHNCs9lsO3k5MDDQVoD+/N+QkJBSDVoedHNCqQisVoN73v6BPXm+uLuY+e6xNtwU7mfvWCIidlNmNyf8U6tWrQgJCeGNN97g9ddfx9/fnzlz5tCnTx/CwsIIDAzknnvuKenmReQyli9fxvJRA8k5uJncAiuPRm9mb1KGvWOJiDi8EpeejRs38s477zBp0iR+/PFH3njjDXbs2EFGRgbr169n5MiR1KlTpzSzigjQtWtXBg8ayOl575N/Yg/puYUMnLyR42ez7R1NRMShXfPcW5mZmbzzzjt89dVXPPvss4wYMQIPD4/SyleudHhLKorCwkJ69erF0lW/EDzgI5z86nCdvyczn25LTS83e8cTESlXZX54609eXl6MGTOGLVu2sHPnTurXr8+0adOudbMicgXOzs7ExsbSOqopJ3/4N0ZWKodSshgcvZnMvEJ7xxMRcUjXXHrg/G+deXl59OvXjzp16jB48GDOnDlTGpsWkcvw8PBg4cKF1A/x58T04ZCXyY7jaTz9bRz5hZquQkTkf5W49Lz//vv079+f5s2b4+HhQbt27fjiiy9o3bo1EydOxNfXtzRzisgl/HnX5kAPE/3rpOHh6sTaAykMjY3HatV0FSIif3VNl6yHh4czaNAg+vXrR8OGDUs7W7nTOT1SUWVmZuLl5cWa/ad5NHozBRaDQW3r8laPJpoXT0QqveJ+f5e49HTo0IH4+HgyMjLw9PSkefPmtGjRwvZo2rQpTk5OJf4A9qDSI5XBD+v2M2zhXsDEy3c25Pk7Gtg7kohImSrz0vOn/fv3ExcXx9atW22Pc+fO4ebmRrNmzdi0adO1bL5cqfRIRZefn0+rVq047ByG351PA/Cf3k3p36aunZOJiJSd4n5/l3gaij81aNCABg0a8OCDD9qWHT58mC1btrBt27Zr3byIXAVXV1eefvppnn32Wcye1ane7kFGzNuJn4crdzcLsnc8ERG7KlHpSUpKokaNGri5Xfp+IPXq1aNevXr84x//AODQoUNcd911JU8pIsX2zDPPcOLECUaNGoWzZ3W8Irvywox4fD1caHe9v73jiYjYTYmu3po1axZ+fn707t2bKVOmcPr06YvGbNy4keHDh9OkSRMiIyOvOaiIFN8777zD448/TuryL8g98Bv5FitPTotjZ2KavaOJiNhNic/pOXDgAAsWLGD+/Pn89ttvtGrVinvuuYfDhw+zaNEiALp160bPnj258847cXd3L9XgZUHn9EhlUlhYyH333cfCJcsI6fcfnEMa4+/lyqyn2xHu72nveCIipabcTmQGSE1NZdGiRSxZsoTw8HB69uxJ27ZtK9ylsio9UtlkZ2fTuXNnEpNTqf/UOA6eySfMz4NZT7elto/j/yIiIlIc5Vp6KguVHqmMzpw5Q05ODq4+/vSZsJ6jqdlEBPkQ89TN+Li72DueiMg1K7e5t0TEsfn5+RESEkItbze+fbQNvm4mdp9M5/GpW8gtsNg7nohIuVHpEalCNv60mP1fv4jZks+mw2cY8sM2Ci2ap0tEqgaVHpEqJDAwEGtqAidi3sRsWFixK5k35u1ER7lFpCpQ6RGpQm699VZ++OEHChJ3kTT3fUwYzNh8jI9W7LV3NBGRMqfSI1LF9O7dmy+++IKc/RtIWfo5AON/PsjktYftnExEpGyp9IhUQU899RQjR44kc8cKzv06DYB3Fu1ifnyinZOJiJQdlR6RKmrkyJE8+eSTpG2IJSz3IAAvx25n9d5Tdk4mIlI2VHpEqiiTycT48eOZOnUqP499jp5RwRRaDZ75bivbEs7aO56ISKlT6RGpwpydnRk4cCBOTmY+7BPJrQ38ySmw8Gj0Zg6cyrB3PBGRUqXSIyIAWAvzyVj6MaYzRzmbXcDAbzZx4lyOvWOJiJQalR4RAc5PVxG3cT0J3w3DKSuFE2m5DJy8ibNZ+faOJiJSKlR6RASA4OBgli9fjq+7E0envopLQSYHTmXy6NTNZOcX2jueiMg1U+kREZvGjRuzcOFCXAoyODr1VZyt+WxLOMez32+lQNNViEgFp9IjIkW0b9+emJgYLGcTOfb9cJywsHrvaf41awdWq6arEJGKS6VHRC7So0cPJkyYQP6JPaTMex8nM8zdlsh/luzWPF0iUmE52zuAiDimJ554grS0NNq2bcspj7q8FLOdb9Yext/LjWc6Xm/veCIiV02lR0Qu65VXXrH9d2pmPqMW7+aDZXuo6enKA61C7ZhMROTqOfThrfHjxxMeHo67uztt2rRh06ZNxVpvxowZmEwmevXqVbYBRaqQVj4ZOO//CYDX5+xg5a5kOycSEbk6Dlt6YmJiGDp0KCNHjmTr1q1ERkbSpUsXTp268rxAR44c4ZVXXuHWW28tp6QiVcOIESM4OOdjLPvXYjXguelb2XT4jL1jiYgUm8OWno8//pgnnniCwYMH07hxYyZMmICHhweTJ0++7DoWi4X+/fvz9ttvc9111/3te+Tl5ZGenl7kISKXNm3aNJo2bcrxuWMwn9xJXqGVx6ZuZvdJ/f9GRCoGhyw9+fn5xMXF0blzZ9sys9lM586d2bBhw2XXe+edd6hduzaPPfZYsd5n9OjR+Pr62h6hoTpHQeRyqlevzrJlywgLrcOR6W/ici6BjNxCBk7exLEz2faOJyLytxyy9KSkpGCxWAgICCiyPCAggKSkpEuus3btWr755hsmTZpU7PcZNmwYaWlptsexY8euKbdIZRcSEsKyZcuo4ePFwehXcctJ5XRGHgO+2UhKZp6944mIXJFDlp6rlZGRwYABA5g0aRL+/v7FXs/NzQ0fH58iDxG5soiICBYtWoS72crByS/hZcrjSGo2j0zZREZugb3jiYhclkOWHn9/f5ycnEhOLnp1SHJyMoGBgReNP3jwIEeOHKF79+44Ozvj7OzMtGnTWLBgAc7Ozhw8eLC8ootUCW3btiU2NpYut93MjKdvoaanKzsT03nq2zjyCi32jicickkOWXpcXV1p2bIlq1atsi2zWq2sWrWKtm3bXjT+hhtu4Pfffyc+Pt726NGjB506dSI+Pl7n6oiUgXvvvZdFixbRtG5toge3xtPVifUHU3kpJh6LpqsQEQfksDcnHDp0KIMGDeKmm26idevWfPrpp2RlZTF48GAABg4cSEhICKNHj8bd3Z2mTZsWWb969eoAFy0XkdJjMpkAaBriQ0eX/SwrqMeS35Oo4bGTUb2a2l4XEXEEDlt6+vbty+nTp3nzzTdJSkoiKiqKZcuW2U5uTkhIwGx2yB1VIlXOnDlz+OLNIfg06YDfva/y/cYE/L3ceOnOhvaOJiJiYzI0e6BNeno6vr6+pKWl6aRmkatgsVh48MEHmTVrFv4398azw/nbRrzbswkD2obbN5yIVHrF/f522D09IlJxODk58e2335KSksLq1XNx8/XHOaonby74gxqertzbPNjeEUVEHPNEZhGpeNzd3Zk3bx7NmzcncfkkzAfXYhjwUkw8a/en2DueiIhKj4iUHl9fX5YuXUrdunU5PHsMbsl/UGAxeOrbLew4fs7e8USkilPpEZFSFRwczPLly6ldy5/nW/vQvn5NsvItDJ6ymUOnM+0dT0SqMJ3I/Bc6kVmk9KSnp+Pj40NmXiH9Jv7G74lphFSvxpxn2xHg427veCJSiRT3+1t7ekSkTPz5D4+XmzNje1xPDecCEs/lMPCbTaRla7oKESl/Kj0iUqZycnLo2fV2fv/8KTzNhexNzuDxaZvJLdB0FSJSvlR6RKRMVatWjUGDBmFJP8X+SS/g7mSw+chZnpu+lUKL1d7xRKQKUekRkTL32muvMWTIEApSjnL8+3/jYoYfd59i2Jzf0WmFIlJeVHpEpMyZTCY++eQTHnjgAbKO7iB14RjMJpgZd5wPlu21dzwRqSJUekSkXJjNZqZNm8btt9/O2Z2/kvvrZAAm/HKQr9ccsnM6EakKVHpEpNy4ubkxd+5cIiMjcT+5jSda+QMwavFu5mw9bud0IlLZae4tESlXPj4+LFu2DICAgAAMt918vfYwr87aQQ0PVzrdUNvOCUWkstKeHhEpd4GBgQQGBmIymRh+TwS31nHFYjV45vs44o6etXc8EamkVHpExK5++GE6M17uTo2cRHILrDwavZl9yRn2jiUilZBKj4jYVc2aNcFqYfuXQ6hppJGWU8DAbzaReC7H3tFEpJJR6RERu+ratSuTJ0/GKMgj/rOn8HPKIyk9lwHfbORMVr6944lIJaLSIyJ2N2DAAD744AOsuZns+PxJfF2sHDqdxeApm8jKK7R3PBGpJFR6RMQhvPrqq7z44otYMlLZN/F5PF1g+/E0nv4ujvxCTVchItdOpUdEHILJZGLs2LH069eP3FNHaZnxGx6uTqzZn8IrM7djtWq6ChG5Nio9IuIwzGYz0dHRTJ48mamfvMOEh1vibDaxYPsJ3lm0S/N0icg1UekREYfi6urK4MGDMZlM3NawFh/e3wyA6PVHGP/zATunE5GKTKVHRBxWdnY20e88h8eexQB8tGIfP2xKsHMqEamoNA2FiDis5ORkNmzYwOnTp4ny9edsUBv+Pfd3ani40LVpkL3jiUgFoz09IuKw6tWrx5IlS/D09CR+2rsEZO7HasCQGfFsOJiKxWqw4WAq8+MTbc9FRC7HZOjMQJv09HR8fX1JS0vDx8fH3nFE5P+tWLGCbt26UWix0vaVrzlhro27sxkvd2dSMi/cwDDI152R3RtrL5BIFVPc72/t6RERh3fXXXcRHR0NhpUNHz9FDXMuuYXWIoUHICktl2e+28qynSftE1REHJpKj4hUCP379+ejjz4Cq4WU9KxLjjH+//HmvN91qEtELqITmUWkwnj55ZeJP5HFGpeaVxx3KrOARRt307Nt43JKJiIVgfb0iEiFctOtdxRr3Lz4JFIy88o4jYhUJNrTIyIVSg334v2u9vPRHNq8t4r29f3pERlMlyYBeLu7lHE6EXFkKj0iUqFE+LtSmH4aJ++amEwXFyDDsGLNy6ZhUHUOni3k132n+XXfaYbPNXPHDbXpERlMpxtq4+7iZIf0ImJPKj0iUqE4mU2cWTWRWr2GYxjWIsXHMKyAidSl/6VGs2B++nIKi3YksWD7CQ6cymTpziSW7kzCy82ZLk0C6REVTPvra+LspCP9IlWB7tPzF7pPj4jj27p1Ky1btqRaw7b43fEkzj61bK8Vpp/mzKqJ5OzbQFRUFNu2bQPAMAzm/7KFXVmeLP49icRzObZ1anq60q15ED0ig2kRVgOz2VTun0lErk1xv7+1p0dEKqScfRtI3L8RtzpNcPKqgSXzLHnH/wDDCsBTTz1lG5uamkqfzm2pVasW9/fpw6Od7+Ow1Z+lO5NIzcpn2oajTNtwlJDq1egeGUyPyGAigrwxmVSARCoT7en5C+3pEXF8f+7p+TtxcXG0aNECgF9++YVevXpx7tw52+vBwcHc/48HiOjYi/15vqzYlUxmXqHt9fq1vegZGUyPqGDq1vQs9c8hIqWnuN/fKj1/odIj4vgSEhJo1KgRubm5lx3j7u7O3r17CQsLsy3Lz8/nxx9/JDY2lnnz5pGWlmZ7bfLkyfR7eCA/7TnFgvgT/LT3FPmFVtvrkXV86REVwr3NgwjwcS+bDyYiJabSUwIqPSIVQ0JCAikpKZd93d/fv0jh+V95eXmsWLGC2NhYFi9ezN69e6lV6/y5QdOnT2fTth2EtLmHHWlurP/LRKYmE9xcryY9o4K5u2kQvh66BF7EEaj0lIBKj0jVU1BQgIvLhfJy2223sWbNGuD8LO89HuhPrRvvYusZJ+KOnrWNc3Ey0aFhLXpEhdA5ojYerjpFUsReVHpKQKVHRObOnUtMTAwLFy4kOzvbtrx+/fr06DeY+rf3ZX58InuSMmyvVXNx4s7GAfSIDOa2hrVwddYl8CLlSaWnBFR6RORPWVlZLFmyhJiYGBYvXkxubi733HMPixcvBmB/cgZTV//Br0dzSDhzoRz5VnPhnmaB9IgMoXU9P5x0CbxImVPpKQGVHhG5lMzMTBYtWkRAQACdOnUC4OjRo4SHhxPRuDG39xmM03VtWH88j9MZF+b7CvBx497mwfSMCqZZiK8ugRcpIyo9JaDSIyLFNWfOHPr160d+fr5tWdNmzbn1vsEYYS1Zl5BNeu6FS+DDa3rQ4/8vga9f29sekUUqLZWeElDpEZGrkZaWxvz584mNjWXFihUUFBTYXps9bwHeDVqzYPsJVu5KIrfgwiXwjYN86BEVTPfIYEKqV7NHdJFKRaWnBFR6RKSkzp49y7x584iNjWXdunUkJibi7X1+j87X0d/y27FssmrewObjWRRaL/yz2yq8Bj0ig7mnWRA1vdzsFV+kQlPpKQGVHhEpDdnZ2Xh4eNieR0VFsX37dgBatr2NZvcM4KxPfbafzOLPf4GdzCZuqe9Pj8hg7moSgLe77gEkUlwqPSWg0iMipc1qtfL1118TExPD6tWrsVovHOZqddudNL1nIKc9wvk98cIdot2czdwRUZsekcF0bFQbdxcne0QXqTBUekpApUdEylJycjJz5swhNjaWX375BcMwePDBB/nhhx84dDqThdtPMicugaNnL0yx4e3mTJemgfSIDKbd9TVxdtI9gET+l0pPCaj0iEh5OXnyJLNnzyYyMpJbb70VgJ07d9K8eXPadLmfOrf05qjhz6nMCydH+3u50q1ZED2igmkRVkOXwIv8P5WeElDpERF7+uqrr3j66adtz81mJ1p360dgm3s5XFidczkXLoGvU6Ma3SOD6REZzA2B3ipAUqWp9JSASo+I2FtCQgKzZs0iNjaWjRs32pabnV0YN+tH9uf6sPyPJLLyLbbXGtT2omdUMD0iQwir6XGpzYpUaio9JaDSIyKO5MiRI8ycOZPY2FgOHjxIUlISrq6u5BZYGPFlDDvS3DiUW40Cy4V/xiNDq9MzMph7mwdR28fdjulFyo9KTwmo9IiIozpz5gx+fn4AGIZB/fr1OXToEC6evtzU81HcG93KkRxX/rwFkNkEba+vSY/IYLo2CcLXQ5fAS+Wl0lMCKj0iUhHk5+fz0UcfERsba7v/D4Cbrz839nwctwbtOJJ54SovFycTHRrWpmdUMJ0jAqjmqkvgpXJR6SkBlR4RqWj27NljOwS2c+dOAJ588kneGP0JC7afYMH2E+xNyrCN93B14s7GAfSMCubWBrVw0SXwUgmo9JSASo+IVGR//PEHM2fOpGvXrtx8880ArF+/nrv+MYgm9zxCQVBzzuRfKDnVPVy4u2kQPaOCaR3uh9msK8CkYlLpKQGVHhGpbN5++23eeust23Pv8GZEdB1Ads0IMgovlJxAH3fubR5Ez6gQmob46BJ4qVBUekpApUdEKhvDMNi+fTuxsbHExMRw6NCh8y+YzFRv2Jp/vPoBG47lkJF74R5A9fw9bfcAql/by07JRYpPpacEVHpEpDIzDINt27YRExNDbGws+fn5HDt2jAKrwS97T/PfeevYl+VGgfXCXp4mwT70iAyme2QwwdWr2TG9yOWp9JSASo+IVBWGYZCYmEidOnUAKCwsJCQkhNPnMqjZrCOht/TmXLUgrMaFAtQ63I/uUcF0axaEn6ervaKLXESlpwRUekSkqkpPT2fkyJHMnDmTxMREAMzVfKgZeQdBbXuQ5lrLNtbZbOKWBv70iAzmriaBeLk52yu2CKDSUyIqPSJS1VmtVjZs2EBsbCwzZ87k5MmTAPzz1X8Tee+jzN+eyM7EdNt4N2cznSMC6B4ZTMdGtXB30T2ApPyp9JSASo+IyAUWi4V169YRGxvL448/TlRUFADRc5YyfMJsqkfeSbbThROdvd2d6dokkB5RwbS9ribOugeQlJPifn879N/I8ePHEx4ejru7O23atGHTpk2XHTtp0iRuvfVWatSoQY0aNejcufMVx4uIyJU5OTlx2223MW7cOFvhAdi8ahEnV37D7o8e5GT0CxT8vhQ36/krwGbGHWfAN5u4efRPvLXgD+KOnkW/W4ujcNg9PTExMQwcOJAJEybQpk0bPv30U2bOnMnevXupXbv2ReP79+9P+/btadeuHe7u7nzwwQfMnTuXP/74g5CQkGK9p/b0iIj8vcLCQn755RdiYmKYM2cOqampgAm3Oo3xv+luakbeQVruhVngQ/2q0b15MD2igrkhUP+2Sumr8Ie32rRpQ6tWrRg3bhxw/jhzaGgozz//PK+//vrfrm+xWKhRowbjxo1j4MCBlxyTl5dHXl6e7Xl6ejqhoaEqPSIixVRQUMBPP/1EbGwsc+bMISAggB07/2DdgVQWbD/B4u3Hyf/LJfCNArzpEXX+HkChfh52TC6VSYUuPfn5+Xh4eDBr1ix69eplWz5o0CDOnTvH/Pnz/3YbGRkZ1K5dm5kzZ3Lvvfdecsxbb73F22+/fdFylR4RkauXn59PQkIC9evXByAnJ4fawXWwBjbGv0VXnOo0x2q6cFbFjWHV6REZTLfmQdT2drdXbKkEKvQ5PSkpKVgsFgICAoosDwgIICkpqVjbeO211wgODqZz586XHTNs2DDS0tJsj2PHjl1TbhGRqszV1dVWeOD8F1Hv7t1wStxOwvQ3OPLffqQs+S9G0m4wDLYlnOPthbu4+b1VPPz1RmI3HyMtp8COn0Aqu0p5c4X333+fGTNmsHr1atzdL//bg5ubG25ubuWYTESk6ggICGDatGnk5uayfPlyYmJiWLBgAQm/r8TsWZ1/DH2fvMCmbEs4x9oDKaw9kMIb83bSsVEtekQFc8cNAVRz1SXwUnocsvT4+/vj5OREcnJykeXJyckEBgZecd2PPvqI999/nx9//JHmzZuXZUwRESkGd3d3evbsSc+ePcnJyWHp0qXExsbyzsOdaNiwIQmp2bz3/XKW7U4l3zeIFbuSWbErGU9XJ+5qEkiPyGBuaeCPiy6Bl2vkkOf0wPkTmVu3bs3nn38OnD+ROSwsjOeee+6yJzKPGTOG//znPyxfvpybb775qt9TV2+JiNjHgw8+SExMDC7+dfFs3IHqze/A8Kxpe72Ghwt3NwuiZ2QwrcL9MJs1C7xcUKFPZIbzl6wPGjSIr776itatW/Ppp58SGxvLnj17CAgIYODAgYSEhDB69GgAPvjgA958802mT59O+/btbdvx8vLCy6t4swSr9IiI2EdmZiaLFy8mNjaWJUuWkJubi2vwDXhG3IZPs07g5m0bG+TrbpsFvkmwDyaTClBVV+FLD8C4ceP48MMPSUpKIioqis8++4w2bdoA0LFjR8LDw4mOjgYgPDyco0ePXrSNkSNH8tZbbxXr/VR6RETsLyMjg4ULFxIbG8vSpUtpHhXFx98tYkH8CZb9kURGbqFt7HX+nrZL4K+rVbxfcKXyqRSlp7yp9IiIOJa0tDQSExNp3LgxAMkpZ2hwWy/cGrbDs8HN4ORiG9s0xIcekcF0jwwmyLeavSKLHRT3+9shT2QWEREB8PX1xdfX1/Y8M+0st4R7sXLxWFLNLng0uBnPiA5Uq3cjOxPT2ZmYzuile2gV7kfPqGDuaRpEDU9XO34CcSTa0/MX2tMjIlIxpKamMm/ePGJjY1m1ahWGqycejdrTvPujJOZf2MvjbDZxawN/ekQFc2fjQLzc9Lt+ZaTDWyWg0iMiUvGkpKQwZ84cYmNjmTZtGoZHDRZtP8GUn34nKe/C4S93FzN3RATQIzKYjo1q4easewBVFio9JaDSIyJSedx11138HLcbz8a34RnRARe/C5NPe7s7c3fTQHpEhtD2+po46RL4Ck2lpwRUekREKo+TJ08ye/ZsYmNjWbt2LS61r8OzcUc8Im7F2dvfNq6WtxvdmgXRIyqYG0Or6xL4CkilpwRUekREKqfExERmz55NTEwM69dvoF3PAdw28BWW7jzJuewL832F+XnQPTKIHpEhNAr0vsIWxZGo9JSASo+ISOV37Ngxzp49S/PmzckvtDJ/416e/WAy1eq3wex64SToGwK9bTdBDPXzsGNi+TsqPSWg0iMiUvVs3bqVp556irjtv1Pt+tZ4Nu5AtetaYvrLPYBahFWnR2Qw3ZoHU8tbE1U7GpWeElDpERGpug4dOkRsbCyxsbFs37UPj0bt8IjoSLXw5sD583zMJmhf35/ukcF0bRqIj7vLlTcq5UKlpwRUekREBGD//v3MnDmTuXPnMn3uYtYcyWL+9hNsP3bONsbVyUSnG2rTIzKEOyJq4+6iS+DtRaWnBFR6RETkSm7qeDf7833wbNwBV/+6tuUerk50bRJI96hgbqnvj4uT2Y4pqx6VnhJQ6RERkSvZtWsXsbGxzIiJ4VBqLp6NO+AZcRvOvgG2MX6ertzT7Pw9gG6qWwOz7gFU5lR6SkClR0REisMwDHbu3GkrQMeynWh276NYQ6JIzcq3jQv0caNnVAjdI4NpEuyjewCVEZWeElDpERGRq2UYBjt27ACgSdNmrD+Yyndr9rBs50nMbp62cfVqetDrxjr0iAqmnr/n5TYnJaDSUwIqPSIiUhp+/vlnHn3iKZKd/PGM6IBH/daYnC/M9t4kyJveLepwb/NgAn3d7Zi0clDpKQGVHhERKS2GYbBly5bzl8HPWUBqtTp4Nu6Ae3gUJvP5K71MJmgd7kfPqBDubhpIDU/Xv9mqXIpKTwmo9IiISFkwDIONGzcSGxvLklVrePXz6Sz94zRbjp61jTGb4LYG/vRuUYfOEQF4ujnbMXHFotJTAio9IiJSno6dyabToJfJ9o/ANeB623IXM9zVOJBeLepwW0N/3Jx1D6ArUekpAZUeEREpT4ZhsHbtWmJiYpi9ch05tZvgGdEBF79g2xgfd2fubhpEz6hg2lxXEyddAn8RlZ4SUOkRERF7sVgs/Prrr8TExjLvl63kBzXHv0UX8p0uTIJa3c1Er5Zh9IwKISq0ui6B/38qPSWg0iMiIo6gsLCQ1atXU9O/FjneISzcfoKF246TWXDhK7u2h4l/tL6OXjeG0CDA245p7U+lpwRUekRExFHNiJ3FkPe+xAhrQbX6N2N2vXCpe5iPmb5t69MjMoRQPw87prQPlZ4SUOkRERFHlp+fz6pVq/ghdjZLf0/EVLcV1a5rgcnpwmzvLevWoEdkMPc0C6KWt5sd05YflZ4SUOkREZGKIi8vjx9//JHvYuey6UQerfo8y8YjZ7B9qxsGkQFuPHzbDXRpGoiPu8sVt1eRqfSUgEqPiIhURFarFbPZTFJaLgu3JzJq2lJM/vVsr5ux0jrEgwEdIrgjIgB3l8p1CbxKTwmo9IiISEVnsViYP38+0+YsZd3xPFzqt8XVP8z2uovJQvcbw+gRGcwt9f1xdjLbMW3pUOkpAZUeERGpTLKzs1m0aDFT569kU7KBW4N2OPvWtr3u5+lK2xBXBnVqyk11/TBX0HsAqfSUgEqPiIhUVpmZmSxctBjngPrsSHNjye8nSc3Kt73uSR53N63N4Nub0TjIp0LdA0ilpwRUekREpKootFgZ9t+pTF39B+7Xt8HsduFSd19TLj1vDObR25sT7u9px5TFo9JTAio9IiJS1aSlpTF73gKiV2xhb6437vVaYnK+MNt7ZB1fukcG0z0ymAAf9ytsyX5UekpApUdERKqys2fPEjNnAdN+iuekax3MQRFYrLZr4Al2zuLBtvUZ2KkZ1T1cr7it8qTSUwIqPSIiIucVFBSQlmdlye8nmbf1OFuPpV140Woh3C2bAR0i6HdbEzxcne0XFJWeElHpERERuVhubi6fffM9MRv2c9K1Dq61L9wDyGTJp6kfvNCzHbc1rIWr84VL4BMSEkhJScFiNdidks/ZXCs13M1E+LviZDbh7+9PWFjYpd7yqqj0lIBKj4iIyJWdOnWKiTMWMGvzEU57hONSI8j2mm81F+6KqMWtdavRLKAajSNuwBR2I353PImzTy3buML005xZNREjYRt79+695uKj0lMCKj0iIiLFd/LkScb/sIgs/8asPZbD6Yw822um3HSyju+m2vWtzz//yyXwhmEFTJye9x5rfxhHixYtrilHcb+/7XsQTkRERCqsoKAgRg19AgCL1WDjoVRGTlnEnqxqOLn74FG/zSXXM5nMGIYVvzue/MuJ0mWv4t97WkREROzOyWyiXX1/Vv7nEVa/eDONLQevON5kMuPsU4vdKflXHFeaVHpERESkVF0fXpcu7VsWa+zZXGsZp7lApUdERERKXQ334lWM4o4rDSo9IiIiUuoi/F0pTD/9/yctX8wwrBSmnybCv/xucqjSIyIiIqXOyWzizKqJgOmi4vPn1VtnVk3EqRxndlfpERERkVLn7++PkbCN0/Pew5KRWuQ1S0Yqp+e9h5GwDX9//3LLpEvWRUREpNSFhYWxd+/ey9yRORCnx8eV2h2Zi0ulR0RERMpEWFiYrdS0snMW0OEtERERqSJUekRERKRKUOkRERGRKkGlR0RERKoElR4RERGpElR6REREpEpQ6REREZEqQaVHREREqgSVHhEREakSVHpERESkSlDpERERkSpBpUdERESqBJUeERERqRJUekRERKRKUOkRERGRKkGlR0RERKoElR4RERGpElR6REREpEpQ6REREZEqwaFLz/jx4wkPD8fd3Z02bdqwadOmK46fOXMmN9xwA+7u7jRr1owlS5aUU1IRERFxdA5bemJiYhg6dCgjR45k69atREZG0qVLF06dOnXJ8evXr6dfv3489thjbNu2jV69etGrVy927txZzslFRETEEZkMwzDsHeJS2rRpQ6tWrRg3bhwAVquV0NBQnn/+eV5//fWLxvft25esrCwWLVpkW3bzzTcTFRXFhAkTLvkeeXl55OXl2Z6npaURFhbGsWPH8PHxKeVPJCIiImUhPT2d0NBQzp07h6+v72XHOZdjpmLLz88nLi6OYcOG2ZaZzWY6d+7Mhg0bLrnOhg0bGDp0aJFlXbp0Yd68eZd9n9GjR/P2229ftDw0NLRkwUVERMRuMjIyKl7pSUlJwWKxEBAQUGR5QEAAe/bsueQ6SUlJlxyflJR02fcZNmxYkaJktVo5c+YMNWvWpHXr1mzevPkaPsUFfzZQ7UGSy2nVqlWp/X2riir7z6+ifT5HymuvLOX5vmX5XqW57bL8LjQMg4yMDIKDg684ziFLT3lxc3PDzc2tyLLq1asD4OTkVOp/KD4+Pio9ckll8fetKqnsP7+K9vkcKa+9spTn+5ble1Wk78Ir7eH5k0OeyOzv74+TkxPJyclFlicnJxMYGHjJdQIDA69q/N/55z//WaL1REpCf9+uTWX/+VW0z+dIee2VpTzftyzfy5H+LEuDQ5/I3Lp1az7//HPg/KGnsLAwnnvuucueyJydnc3ChQtty9q1a0fz5s0veyJzeUlPT8fX15e0tDSH+e1HRESkPDnCd6HDHt4aOnQogwYN4qabbqJ169Z8+umnZGVlMXjwYAAGDhxISEgIo0ePBuCFF16gQ4cOjB07lm7dujFjxgy2bNnCxIkT7fkxgPOH0UaOHHnRoTQREZGqwhG+Cx12Tw/AuHHj+PDDD0lKSiIqKorPPvuMNm3aANCxY0fCw8OJjo62jZ85cyZvvPEGR44coUGDBowZM4Z77rnHTulFRETEkTh06REREREpLQ55IrOIiIhIaVPpERERkSpBpUdERESqBJUeERERqRJUehzAokWLaNSoEQ0aNODrr7+2dxwREZFy17t3b2rUqEGfPn3K7D109ZadFRYW0rhxY37++Wd8fX1p2bIl69evp2bNmvaOJiIiUm5Wr15NRkYGU6dOZdasWWXyHtrTY2ebNm2iSZMmhISE4OXlxd13382KFSvsHUtERKRcdezYEW9v7zJ9D5Wea/Trr7/SvXt3goODMZlMzJs376Ix48ePJzw8HHd3d9q0acOmTZtsr504cYKQkBDb85CQEBITE8sjuoiISKm41u/C8qLSc42ysrKIjIxk/Pjxl3w9JiaGoUOHMnLkSLZu3UpkZCRdunTh1KlT5ZxURESkbFSU70KVnmt09913M2rUKHr37n3J1z/++GOeeOIJBg8eTOPGjZkwYQIeHh5MnjwZgODg4CJ7dhITEwkODi6X7CIiIqXhWr8Ly4tKTxnKz88nLi6Ozp0725aZzWY6d+7Mhg0bAGjdujU7d+4kMTGRzMxMli5dSpcuXewVWUREpFQV57uwvDjsLOuVQUpKChaLhYCAgCLLAwIC2LNnDwDOzs6MHTuWTp06YbVa+de//qUrt0REpNIoznchQOfOndm+fTtZWVnUqVOHmTNn0rZt21LNotLjAHr06EGPHj3sHUNERMRufvzxxzJ/Dx3eKkP+/v44OTmRnJxcZHlycjKBgYF2SiUiIlJ+HOm7UKWnDLm6utKyZUtWrVplW2a1Wlm1alWp77ITERFxRI70XajDW9coMzOTAwcO2J4fPnyY+Ph4/Pz8CAsLY+jQoQwaNIibbrqJ1q1b8+mnn5KVlcXgwYPtmFpERKT0VJTvQk1DcY1Wr15Np06dLlo+aNAgoqOjARg3bhwffvghSUlJREVF8dlnn9GmTZtyTioiIlI2Ksp3oUqPiIiIVAk6p0dERESqBJUeERERqRJUekRERKRKUOkRERGRKkGlR0RERKoElR4RERGpElR6REREpEpQ6REREZEqQaVHREREqgSVHhGREpo4cSKhoaGYzWY+/fRT3nrrLaKiooq9/pEjRzCZTMTHx192zOrVqzGZTJw7d+6yY0wmE/PmzSv2+4pUVSo9IlWUyWS64uOtt96yd0SHlp6eznPPPcdrr71GYmIiTz75JK+88kqRmaRFxLFolnWRKurkyZO2/46JieHNN99k7969tmVeXl62/zYMA4vFgrNzxfsno6yyJyQkUFBQQLdu3QgKCrIt/+vPTUQci/b0iFRRgYGBtoevry8mk8n2fM+ePXh7e7N06VJatmyJm5sba9eu5eDBg/Ts2ZOAgAC8vLxo1aoVP/74Y5HthoeH89577/Hoo4/i7e1NWFgYEydOtL2en5/Pc889R1BQEO7u7tStW5fRo0cD8NBDD9G3b98i2ysoKMDf359p06YBYLVaGT16NPXq1aNatWpERkYya9Ys2/g/Dwf9b/bt27fTqVMnvL298fHxoWXLlmzZssW23tq1a7n11lupVq0aoaGhDBkyhKysrEv+7KKjo2nWrBkA1113HSaTiSNHjlzy8NbXX39NREQE7u7u3HDDDXzxxRdX/HNZsmQJDRs2pFq1anTq1IkjR45ccfyfUlJS6N27Nx4eHjRo0IAFCxYUaz2RKsUQkSpvypQphq+vr+35zz//bABG8+bNjRUrVhgHDhwwUlNTjfj4eGPChAnG77//buzbt8944403DHd3d+Po0aO2devWrWv4+fkZ48ePN/bv32+MHj3aMJvNxp49ewzDMIwPP/zQCA0NNX799VfjyJEjxpo1a4zp06cbhmEYixYtMqpVq2ZkZGTYtrdw4UKjWrVqRnp6umEYhjFq1CjjhhtuMJYtW2YcPHjQmDJliuHm5masXr36itmbNGliPPzww8bu3buNffv2GbGxsUZ8fLxhGIZx4MABw9PT0/jkk0+Mffv2GevWrTNuvPFG45FHHrnkzys7O9v48ccfDcDYtGmTcfLkSaOwsNAYOXKkERkZaRv33XffGUFBQcbs2bONQ4cOGbNnzzb8/PyM6OhowzAM4/DhwwZgbNu2zTAMw0hISDDc3NyMoUOHGnv27DG+++47IyAgwACMs2fPXvbPDzDq1KljTJ8+3di/f78xZMgQw8vLy0hNTS3Gn75I1aHSIyKXLT3z5s3723WbNGlifP7557bndevWNR5++GHbc6vVatSuXdv48ssvDcMwjOeff964/fbbDavVetG2CgoKDH9/f2PatGm2Zf369TP69u1rGIZh5ObmGh4eHsb69euLrPfYY48Z/fr1u2J2b29vW9n4X4899pjx5JNPFlm2Zs0aw2w2Gzk5OZdcZ9u2bQZgHD582Lbsf0vP9ddfbyt0f3r33XeNtm3bGoZxcekZNmyY0bhx4yLjX3vttWKVnjfeeMP2PDMz0wCMpUuXXnYdkapIh7dE5LJuuummIs8zMzN55ZVXiIiIoHr16nh5ebF7924SEhKKjGvevLntv/88bHbq1CkAHnnkEeLj42nUqBFDhgxhxYoVtrHOzs488MADfP/99wBkZWUxf/58+vfvD8CBAwfIzs7mzjvvxMvLy/aYNm0aBw8evGL2oUOH8vjjj9O5c2fef//9IuO3b99OdHR0kW126dIFq9XK4cOHS/Szy8rK4uDBgzz22GNFtjtq1KiLsv5p9+7dtGnTpsiytm3bFuv9/voz9/T0xMfHx/YzF5HzKt5ZiSJSbjw9PYs8f+WVV1i5ciUfffQR9evXp1q1avTp04f8/Pwi41xcXIo8N5lMWK1WAFq0aMHhw4dZunQpP/74Iw888ACdO3e2nZfTv39/OnTowKlTp1i5ciXVqlWja9euwPnSBbB48WJCQkKKvIebm9sVs7/11ls89NBDLF68mKVLlzJy5EhmzJhB7969yczM5KmnnmLIkCEX/QzCwsKK9bP6X39mnTRp0kVFxsnJqUTbvJIr/cxF5DyVHhEptnXr1vHII4/Qu3dv4PwXe3FPtP0rHx8f+vbtS9++fenTpw9du3blzJkz+Pn50a5dO0JDQ4mJiWHp0qX84x//sH2hN27cGDc3NxISEujQocNVv2/Dhg1p2LAhL730Ev369WPKlCn07t2bFi1asGvXLurXr3/V27ycgIAAgoODOXTokG1P1d+JiIi46ATk3377rdQyiVR1Kj0iUmwNGjRgzpw5dO/eHZPJxIgRI656b8LHH39MUFAQN954I2azmZkzZxIYGEj16tVtYx566CEmTJjAvn37+Pnnn23Lvb29eeWVV3jppZewWq3ccsstpKWlsW7dOnx8fBg0aNAl3zMnJ4dXX32VPn36UK9ePY4fP87mzZu5//77AXjttde4+eabee6553j88cfx9PRk165drFy5knHjxl39D+r/vf322wwZMgRfX1+6du1KXl4eW7Zs4ezZswwdOvSi8U8//TRjx47l1Vdf5fHHHycuLo7o6OgSv7+IFKVzekSk2D7++GNq1KhBu3bt6N69O126dKFFixZXtQ1vb2/GjBnDTTfdRKtWrThy5AhLlizBbL7wz1H//v3ZtWsXISEhtG/fvsj67777LiNGjGD06NFERETQtWtXFi9eTL169S77nk5OTqSmpjJw4EAaNmzIAw88wN13383bb78NnD8f5pdffmHfvn3ceuut3Hjjjbz55psEBwdf1Wf7X48//jhff/01U6ZMoVmzZnTo0IHo6OjLZg0LC2P27NnMmzePyMhIJkyYwHvvvXdNGUTkApNhGIa9Q4iIiIiUNe3pERERkSpBpUdERESqBJUeERERqRJUekRERKRKUOkRERGRKkGlR0RERKoElR4RERGpElR6REREpEpQ6REREZEqQaVHREREqgSVHhEREakS/g8wjwN//QkQ9gAAAABJRU5ErkJggg==", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "samples0 = sampler0.sweep()\n", "samples1 = sampler1.sweep()\n", "samples2 = sampler2.sweep()\n", "\n", "exactM2 = np.array([state0 @ M2 @ state0, state1 @ M2 @ state1, state2 @ M2 @ state2])\n", "sampledM2 = np.array([\n", " M2.expectation(state0, samples0),\n", " M2.expectation(state1, samples1),\n", " M2.expectation(state2, samples2),\n", "])\n", "\n", "plt.semilogx(h, exactM2, 'ks--', label='Exact')\n", "plt.semilogx(h, sampledM2, 'o-', label='Sampled')\n", "plt.xlabel('Transverse field h')\n", "plt.ylabel(r'$\\langle M^2 \\rangle$')\n", "plt.ylim(0, 1)\n", "plt.legend()\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "f73e9b3d", "metadata": {}, "source": [ "By default, the initial thermalization takes 20 * N updates, and each sweep takes 2 * N updates. This is an empirical balance between efficiency and reliability. For better thermalization, one can specify `thermal_steps` and `sweep_steps` arguments in the initialization.\n", "\n", "One can also customize metropolis samplers by an inheritance of `qtx.sampler.Metropolis`. Here we show an example of defining a sampler that flips all spins in a lattice." ] }, { "cell_type": "code", "execution_count": 8, "id": "623f7a85", "metadata": {}, "outputs": [], "source": [ "from functools import partial\n", "from jaxtyping import Key\n", "\n", "\n", "class FlipLattice(qtx.sampler.Metropolis):\n", "\n", " @partial(jax.jit, static_argnums=0)\n", " def propose(self, key: Key, old_spins: jax.Array):\n", " new_spins = -old_spins\n", " return new_spins" ] }, { "cell_type": "markdown", "id": "07ee768a", "metadata": {}, "source": [ "The sampler we just defined is obviously not ergodic. To make it work, we want to conbine it with `LocalFlip`. This can be done as shown below." ] }, { "cell_type": "code", "execution_count": 9, "id": "0859e337", "metadata": {}, "outputs": [], "source": [ "sampler_local = qtx.sampler.LocalFlip(state1, nsamples=896, thermal_steps=0)\n", "sampler_global = FlipLattice(state1, nsamples=128, thermal_steps=0)\n", "sampler_mix = qtx.sampler.MixSampler([sampler_local, sampler_global])" ] }, { "cell_type": "markdown", "id": "b8a66c96", "metadata": {}, "source": [ "`sampler_mix` generates 896+128=1024 samples per sweep. In each proposal, there is 896/1024 chance to update all spin configurations by `LocalFlip`, and 128/1024 to update by `FlipLattice`. Here we check the magnetization to verify that the generated samples are still correct." ] }, { "cell_type": "code", "execution_count": 10, "id": "99f750f4", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Exact magnetization for h=1.0: 0.29388403111101535\n", "Sampled magnetization for h=1.0: 0.30169677734375\n" ] } ], "source": [ "samples = sampler_mix.sweep()\n", "\n", "print(\"Exact magnetization for h=1.0:\", state1 @ M2 @ state1)\n", "print(\"Sampled magnetization for h=1.0:\", M2.expectation(state1, samples1))" ] }, { "cell_type": "markdown", "id": "e7d48309", "metadata": {}, "source": [ "## Measurement\n", "\n", "Now we introduce how to use the generated samples to measure physical quantities.\n", "\n", "As we have shown, one can call `op.expectation(state, samples)` to estimate the expectation value\n", "$\\braket{\\hat O}$ of an operator. To estimate the varaince $\\braket{\\hat O^2} - \\braket{\\hat O}^2$ at the same time, one can call" ] }, { "cell_type": "code", "execution_count": 11, "id": "989ca020", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "h=1.0: = 0.9398193359375 var(M2) = 0.01366935670375824\n" ] } ], "source": [ "expectation, var = M2.expectation(state0, samples0, return_var=True)\n", "print(\"h=1.0: =\", expectation, \"var(M2) =\", var)" ] }, { "cell_type": "markdown", "id": "61eef379", "metadata": {}, "source": [ "As we know, the uncertainty of an observable $O$ in Monte Carlo methods is\n", "\n", "$$\n", "u(O) \\approx \\sqrt{\\frac{\\mathrm{Var}(O)}{N_s}},\n", "$$\n", "\n", "where $N_s$ is the number of samples. Then we can estimate the uncertainty of the observable as" ] }, { "cell_type": "code", "execution_count": 12, "id": "7b53a83d", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "h=1.0: = 0.9398193359375 ± 0.0036536257547830353\n", "Exact M2 = 0.9399276867036591\n" ] } ], "source": [ "uncertainty = jnp.sqrt(var / samples0.nsamples)\n", "print(f\"h=1.0: = {expectation} ± {uncertainty}\")\n", "print(\"Exact M2 =\", exactM2[0])" ] }, { "cell_type": "markdown", "id": "e2008d26", "metadata": {}, "source": [ "Apart from observables, some other quantities are also important in quantum physics. An important one is the overlap (fidelity) between two states. One can rewrite it as\n", "\n", "$$\n", "\\begin{aligned}\n", "\\mathcal{F}(\\psi, \\phi) &= \\frac{\\braket{\\psi|\\phi} \\braket{\\phi|\\psi}}{\\braket{\\psi|\\psi} \\braket{\\phi|\\phi}} \\\\\n", "&= \\frac{\\sum_s \\psi^*(s) \\phi(s) \\sum_{s'} \\phi^*(s') \\psi(s')}{||\\psi||^2 ||\\phi||^2} \\\\\n", "&= \\sum_s \\frac{|\\psi(s)|^2}{||\\psi||^2} \\frac{\\phi(s)}{\\psi(s)} \\sum_{s'} \\frac{|\\phi(s')|^2}{||\\phi||^2} \\frac{\\psi(s')}{\\phi(s')} \\\\\n", "&= \\left< \\frac{\\phi(s)}{\\psi(s)} \\right>_{s \\sim |\\psi|^2} \\left< \\frac{\\psi(s)}{\\phi(s)} \\right>_{s \\sim |\\phi|^2}\n", "\\end{aligned}\n", "$$\n", "\n", "Here, we show the overlap between `state0` and `state1` as an example." ] }, { "cell_type": "code", "execution_count": 13, "id": "4c8d6692", "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Estimated overlap: 0.1494855468512831\n", "Exact overlap: 0.14122972655126168\n" ] } ], "source": [ "samples0 = sampler0.sweep()\n", "samples1 = sampler1.sweep()\n", "\n", "s0 = samples0.spins\n", "psi_s0 = samples0.psi\n", "phi_s0 = state1(s0)\n", "\n", "s1 = samples1.spins\n", "psi_s1 = state0(s1)\n", "phi_s1 = samples1.psi\n", "\n", "print(\"Estimated overlap: \", jnp.mean(phi_s0 / psi_s0) * jnp.mean(psi_s1 / phi_s1))\n", "print(\"Exact overlap: \", abs(state0 @ state1) ** 2)" ] }, { "cell_type": "markdown", "id": "6f67d2b3", "metadata": {}, "source": [ "Another important quantity is the entanglement entropy. In VMC, it's only convenient to compute the Rényi-2 entropy\n", "\n", "$$\n", "S_2(A, B) = -\\log \\mathrm{Tr}(\\rho_A^2)\n", "$$\n", "when we split the whole system into two subregions A and B.\n", "\n", "For a pure state $\\ket{\\psi}$, matrix elements of $\\rho_A$ is given by\n", "$$\n", "(\\rho_A)_{s_A, s_A'} = \\frac{\\sum_{s_B} \\psi(s_A, s_B) \\psi^*(s_A', s_B)}{||\\psi||^2},\n", "$$\n", "where $\\ket{s} = \\ket{s_A} \\otimes \\ket{s_B}$ splits a basis state into basis states in two subregions.\n", "\n", "Then we can rewrite $\\mathrm{Tr}(\\rho_A^2)$ into\n", "\n", "$$\n", "\\begin{aligned}\n", "\\mathrm{Tr}(\\rho_A^2) &= \\sum_{s_A, s_A'} (\\rho_A)_{s_A, s_A'} (\\rho_A)_{s_A', s_A} \\\\\n", "&= \\frac{\\sum_{s_A, s_A', s_B, s_B'} \\psi(s_A, s_B') \\psi^*(s_A', s_B') \\psi(s_A', s_B) \\psi^*(s_A, s_B)}{||\\psi||^4} \\\\\n", "&= \\sum_{s_A, s_A', s_B, s_B'} \\frac{|\\psi(s_A, s_B) \\psi(s_A', s_B')|^2}{||\\psi||^4}\n", "\\frac{\\psi(s_A, s_B') \\psi(s_A', s_B)}{\\psi(s_A, s_B) \\psi(s_A', s_B')} \\\\\n", "&= \\left< \\frac{\\psi(s_A, s_B') \\psi(s_A', s_B)}{\\psi(s_A, s_B) \\psi(s_A', s_B')} \\right>_{s \\sim |\\psi|^2, s' \\sim |\\psi|^2}.\n", "\\end{aligned}\n", "$$\n", "\n", "Therefore, the entropy is\n", "$$\n", "S_2(A, B) = -\\log \\left< \\frac{\\psi(s_A, s_B') \\psi(s_A', s_B)}{\\psi(s_A, s_B) \\psi(s_A', s_B')} \\right>_{s \\sim |\\psi|^2, s' \\sim |\\psi|^2}.\n", "$$\n", "\n", "Here, we show an example of entropy computation, where we split the whole lattice into A and B with equivalent sizes." ] }, { "cell_type": "code", "execution_count": 17, "id": "159a9ee3", "metadata": {}, "outputs": [], "source": [ "def entropy(state: qtx.state.Variational, samples: qtx.sampler.Samples):\n", " Lx, Ly = lattice.shape[1:]\n", " Ns = samples.nsamples\n", " spins = samples.spins.reshape(Ns, Lx, Ly)\n", " sA = spins[: Ns // 2, : Lx // 2]\n", " sB = spins[: Ns // 2, Lx // 2 :]\n", " sA_ = spins[Ns // 2 :, : Lx // 2]\n", " sB_ = spins[Ns // 2 :, Lx // 2 :]\n", " \n", " swap_s1 = jnp.concatenate([sA, sB_], axis=1).reshape(Ns // 2, -1)\n", " swap_s2 = jnp.concatenate([sA_, sB], axis=1).reshape(Ns // 2, -1)\n", " psi1 = state(swap_s1)\n", " psi2 = state(swap_s2)\n", "\n", " Tr_rhoA2 = jnp.mean((psi1 * psi2) / (samples.psi[: Ns // 2] * samples.psi[Ns // 2 :]))\n", " entropy = -jnp.log(Tr_rhoA2)\n", " return entropy" ] }, { "cell_type": "code", "execution_count": 18, "id": "ff430ded", "metadata": {}, "outputs": [], "source": [ "samples0 = sampler0.sweep()\n", "samples1 = sampler1.sweep()\n", "samples2 = sampler2.sweep()\n", "\n", "S0 = entropy(state0, samples0)\n", "S1 = entropy(state1, samples1)\n", "S2 = entropy(state2, samples2)" ] }, { "cell_type": "markdown", "id": "a868b363", "metadata": {}, "source": [ "From the exact wave-funciton, one can utilize quspin to compute the exact entropy." ] }, { "cell_type": "code", "execution_count": 19, "id": "95457e40", "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjkAAAG2CAYAAACd5Zf9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABVL0lEQVR4nO3dd3QU5f7H8fembXpCCGmQkNDBFDrSxIKCIgrXgh25lqsXFIleFRXQK4LlgvwUFeEqoIgJRRClCKJIEaWGXqSGlkAoCUlI253fH8hiLgEXTLKbzed1zp7jPPPMzHd3OdmPzzwzYzIMw0BERETExbg5ugARERGRiqCQIyIiIi5JIUdERERckkKOiIiIuCSFHBEREXFJCjkiIiLikhRyRERExCUp5IiIiIhLUsgRERERl6SQIyIiIi7JaULO0qVL6dmzJ1FRUZhMJmbPnv2n2yxZsoSWLVtiNptp0KABkyZNqvA6RUREpGpwmpCTl5dHUlISH3zwgV399+7dS48ePbjuuutIS0vjmWee4dFHH+W7776r4EpFRESkKjA54wM6TSYTs2bNolevXhft88ILLzB37lw2b95sa7vnnns4deoUCxYsqIQqRURExJl5OLqAK7Vy5Uq6du1aqq1bt24888wzF92msLCQwsJC27LVauXEiRPUrFkTk8lUUaWKiIhIOTIMg9OnTxMVFYWb28VPSlXZkJORkUF4eHiptvDwcHJycjhz5gw+Pj4XbDNy5Ehee+21yipRREREKtCBAweoU6fORddX2ZBzJQYPHkxycrJtOTs7m5iYGA4cOEBgYKADKxMRERF75eTkEB0dTUBAwCX7VdmQExERQWZmZqm2zMxMAgMDyxzFATCbzZjN5gvaAwMDFXJERESqmD+bauI0V1ddrvbt27N48eJSbYsWLaJ9+/YOqkhEREScidOEnNzcXNLS0khLSwPOXiKelpZGeno6cPZU00MPPWTr/8QTT7Bnzx6ef/55tm/fzocffsi0adMYNGiQI8oXERERJ+M0IWfNmjW0aNGCFi1aAJCcnEyLFi0YOnQoAEeOHLEFHoC4uDjmzp3LokWLSEpKYtSoUfz3v/+lW7duDqlfREREnItT3iensuTk5BAUFER2drbm5IiIuDiLxUJxcbGjyxA7eHp64u7uftH19v5+V9mJxyIiIvYwDIOMjAxOnTrl6FLkMgQHBxMREfGX7mOnkCMiIi7tXMAJCwvD19dXN391coZhkJ+fz9GjRwGIjIy84n0p5IiIiMuyWCy2gFOzZk1HlyN2OncrmKNHjxIWFnbJU1eX4jQTj0VERMrbuTk4vr6+Dq5ELte57+yvzKNSyBEREZenU1RVT3l8Zwo5IiIi4pIUckRERMQlaeKxiIjIRaSnp5OVlXXR9aGhocTExFTIsR9++GEmT558QXu3bt1YsGBBhRzzj1599VVmz55texJBVaSQIyIiUob09HQaN25MQUHBRft4e3uzY8eOCgs63bt3Z+LEiaXaynrQtJRNp6tERETKkJWVdcmAA1BQUHDJkZ6/ymw2ExERUepVo0YNlixZgpeXF8uWLbP1ffvttwkLCyMzMxOABQsW0KlTJ4KDg6lZsya33noru3fvLrX/gwcPcu+99xISEoKfnx+tW7fm119/ZdKkSbz22mts2LABk8mEyWRi0qRJFfY+K4pGckREpNrJy8u76Dp3d3e8vb3Ldb9+fn5XtL+Lufbaa3nmmWd48MEH2bBhA3v27GHIkCFMnz6d8PBwWy3JyckkJiaSm5vL0KFD6d27N2lpabi5uZGbm0uXLl2oXbs2c+bMISIignXr1mG1WunTpw+bN29mwYIFfP/99wAEBQWV63uoDAo5IiJS7fj7+1903S233MLcuXOvaL+xsbFljuxc6WMiv/322wtqfemll3jppZcYPnw4ixYt4vHHH2fz5s307duX2267zdbvjjvuKLXdp59+Sq1atdi6dSvx8fFMnTqVY8eOsXr1akJCQgBo0KCBrb+/vz8eHh5ERERcUe3OQCFHRETESV133XV89NFHpdrOBRIvLy+++OILEhMTqVu3Lu+++26pfr/99htDhw7l119/JSsrC6vVCpydaxQfH09aWhotWrSw7c8VKeSIiEi1k5ube9F1V/oIAYB9+/Zd8bZl8fPzKzW68r9+/vlnAE6cOMGJEydKnRbr2bMndevWZcKECURFRWG1WomPj6eoqAg4/+gEV6aJxyIiUu34+fld9HWl83Eutd+KsHv3bgYNGsSECRNo164dffv2tY3WHD9+nB07dvDKK69www030LRpU06ePFlq+8TERNLS0jhx4kSZ+/fy8sJisVRI7ZVFIUdERMRJFRYWkpGRUeqVlZWFxWLhgQceoFu3bvTr14+JEyeyceNGRo0aBUCNGjWoWbMm48ePZ9euXfzwww8kJyeX2ve9995LREQEvXr1YsWKFezZs4eZM2eycuVK4Oz8or1795KWlkZWVhaFhYWV/v7/KoUcERGRMoSGhv7pqI63tzehoaEVVsOCBQuIjIws9erUqRNvvPEG+/fv5+OPPwYgMjKS8ePH88orr7Bhwwbc3NxISUlh7dq1xMfHM2jQIN55551S+/by8mLhwoWEhYVxyy23kJCQwJtvvmk7XXfHHXfQvXt3rrvuOmrVqsWXX35ZYe+zopiMK53y7QJycnIICgoiOzubwMBAR5cjIiLlrKCggL179xIXF3dFp6Ececfj6u5S3529v9+aeCwiInIRMTExCjFVmE5XiYiIiEtSyBERERGXpJAjIiIiLkkhR0RERFySQo6IiIi4JIUcERERcUkKOSIiIuKSFHJERETEJSnkiIiIyGUxmUzMnj37L+3j4YcfplevXuVSz8Uo5IiIiPwJi9Vg5e7jfJ12iJW7j2OxVvwTkY4dO8aTTz5JTEwMZrOZiIgIunXrxooVKyr82K5Cj3UQERG5hAWbj/DaN1s5kl1ga4sM8mZYz2Z0j4+ssOPecccdFBUVMXnyZOrVq0dmZiaLFy/m+PHjFXZMV6ORHBERkYtYsPkIT05ZVyrgAGRkF/DklHUs2HykQo576tQpli1bxltvvcV1111H3bp1adu2LYMHD+a2224DYPTo0SQkJODn50d0dDT//Oc/yc3Nte1j0qRJBAcH8+2339K4cWN8fX258847yc/PZ/LkycTGxlKjRg2efvppLBaLbbvY2Fhef/117r33Xvz8/KhduzYffPDBJes9cOAAd999N8HBwYSEhHD77bezb98+23qLxUJycjLBwcHUrFmT559/nsp4PrhCjoiIVBuGYZBfVGLX63RBMcPmbKGsn+Jzba/O2crpgmK79nc5P+r+/v74+/sze/ZsCgsLy+zj5ubGe++9x5YtW5g8eTI//PADzz//fKk++fn5vPfee6SkpLBgwQKWLFlC7969mTdvHvPmzePzzz/n448/ZsaMGaW2e+edd0hKSmL9+vW8+OKLDBw4kEWLFpVZR3FxMd26dSMgIIBly5axYsUK/P396d69O0VFRQCMGjWKSZMm8emnn7J8+XJOnDjBrFmz7P48rpTJqIwo5aTsfVS7iIhUTQUFBezdu5e4uDi8vb3JLyqh2dDvHFLL1n93w9fL/lkiM2fO5LHHHuPMmTO0bNmSLl26cM8995CYmFhm/xkzZvDEE0+QlZUFnB3J6devH7t27aJ+/foAPPHEE3z++edkZmbi7+8PQPfu3YmNjWXcuHHA2ZGcpk2bMn/+fNu+77nnHnJycpg3bx5wduLxrFmz6NWrF1OmTGH48OFs27YNk8kEQFFREcHBwcyePZubbrqJqKgoBg0axL/+9S8ASkpKiIuLo1WrVhedwPy/390f2fv7rZEcERERJ3THHXdw+PBh5syZQ/fu3VmyZAktW7Zk0qRJAHz//ffccMMN1K5dm4CAAB588EGOHz9Ofn6+bR++vr62gAMQHh5ObGysLeCcazt69GipY7dv3/6C5W3btpVZ54YNG9i1axcBAQG2EaiQkBAKCgrYvXs32dnZHDlyhHbt2tm28fDwoHXr1lf82dhLE49FRKTa8PF0Z+u/u9nVd9XeEzw8cfWf9pvUrw1t40LsOvbl8vb25sYbb+TGG29kyJAhPProowwbNoxrr72WW2+9lSeffJI33niDkJAQli9fziOPPEJRURG+vr4AeHp6ltqfyWQqs81qtV52befk5ubSqlUrvvjiiwvW1apV64r3Wx4UckREpNowmUx2nzLq3LAWkUHeZGQXlDkvxwREBHnTuWEt3N1M5VrnxTRr1ozZs2ezdu1arFYro0aNws3t7EmZadOmldtxfvnllwuWmzZtWmbfli1bkpqaSlhY2EVPHUVGRvLrr79yzTXXAGdPV61du5aWLVuWW81l0ekqERGRMri7mRjWsxlwNtD80bnlYT2bVUjAOX78ONdffz1Tpkxh48aN7N27l+nTp/P2229z++2306BBA4qLi3n//ffZs2cPn3/+uW1OTXlYsWIFb7/9Njt37uSDDz5g+vTpDBw4sMy+999/P6Ghodx+++0sW7aMvXv3smTJEp5++mkOHjwIwMCBA3nzzTeZPXs227dv55///CenTp0qt3ovRiFHRETkIrrHR/LRAy2JCCo98TUiyJuPHmhZYffJ8ff3p127drz77rtcc801xMfHM2TIEB577DHGjh1LUlISo0eP5q233iI+Pp4vvviCkSNHltvxn332WdasWUOLFi0YPnw4o0ePplu3sk/z+fr6snTpUmJiYvjb3/5G06ZNeeSRRygoKLCN7Dz77LM8+OCD9O3bl/bt2xMQEEDv3r3Lrd6L0dVVurpKRMRlXeoKncthsRqs2nuCo6cLCAvwpm1cSKWdoqpssbGxPPPMMzzzzDMOraM8rq7SnBwREZE/4e5mon39mo4uQy6TTleJiIiIS9JIjoiIiNj88XEMVZ1GckRERMQlKeSIiIjLq8bX2FRZ5fGdKeSIiIjLOnd33z8+6kCqhnPf2f/eoflyaE6OiIi4LHd3d4KDg23PZvL19bU9RFKck2EY5Ofnc/ToUYKDg3F3v/zHYZyjkCMiIi4tIiIC4IKHUIpzCw4Otn13V0ohR0REXJrJZCIyMpKwsDCKi4sdXY7YwdPT8y+N4JyjkCMiItWCu7t7ufxwStWhicciIiLikhRyRERExCUp5IiIiIhLUsgRERERl6SQIyIiIi5JIUdERERckkKOiIiIuCSFHBEREXFJCjkiIiLikhRyRERExCUp5IiIiIhLUsgRERERl6SQIyIiIi5JIUdERERcklOFnA8++IDY2Fi8vb1p164dq1atumT/MWPG0LhxY3x8fIiOjmbQoEEUFBRUUrUiIiLizJwm5KSmppKcnMywYcNYt24dSUlJdOvWjaNHj5bZf+rUqbz44osMGzaMbdu28cknn5CamspLL71UyZWLiIiIM3KakDN69Ggee+wx+vXrR7NmzRg3bhy+vr58+umnZfb/+eef6dixI/fddx+xsbHcdNNN3HvvvX86+iMiIiLVg1OEnKKiItauXUvXrl1tbW5ubnTt2pWVK1eWuU2HDh1Yu3atLdTs2bOHefPmccstt1z0OIWFheTk5JR6iYiIiGvycHQBAFlZWVgsFsLDw0u1h4eHs3379jK3ue+++8jKyqJTp04YhkFJSQlPPPHEJU9XjRw5ktdee61caxcRERHn5BQjOVdiyZIljBgxgg8//JB169bx1VdfMXfuXF5//fWLbjN48GCys7NtrwMHDlRixSIiIlKZnGIkJzQ0FHd3dzIzM0u1Z2ZmEhERUeY2Q4YM4cEHH+TRRx8FICEhgby8PB5//HFefvll3NwuzG9msxmz2Vz+b0BEREScjlOM5Hh5edGqVSsWL15sa7NarSxevJj27duXuU1+fv4FQcbd3R0AwzAqrlgRERGpEpxiJAcgOTmZvn370rp1a9q2bcuYMWPIy8ujX79+ADz00EPUrl2bkSNHAtCzZ09Gjx5NixYtaNeuHbt27WLIkCH07NnTFnZERESk+nKakNOnTx+OHTvG0KFDycjIoHnz5ixYsMA2GTk9Pb3UyM0rr7yCyWTilVde4dChQ9SqVYuePXvyxhtvOOotiIiIiBMxGdX43E5OTg5BQUFkZ2cTGBjo6HJERETEDvb+fjvFnBwRERGR8qaQIyIiIi5JIUdERERckkKOiIiIuCSFHBEREXFJCjkiIiLikhRyRERExCUp5IiIiIhLUsgRERERl6SQIyIiIi5JIUdERERckkKOiIiIuCSFHBEREXFJCjkiIiLikhRyRERExCUp5IiIiIhLUsgRERERl6SQIyIiIi5JIUdERERckkKOiIiIuCSFHBEREXFJCjkiIiLikhRyRERExCUp5IiIiIhL8nB0Aa7GYjVYtfcER08XEBbgTdu4ENzdTI4uS0REpNpRyClHCzYf4bVvtnIku8DWFhnkzbCezegeH+nAykRERKofna4qJws2H+HJKetKBRyAjOwCnpyyjgWbjzioMhERkepJIaccWKwGQ2ZtwihjnfH7a+jsTVisZfUQERGRiqCQUw6+/XUbx/KKL9nnaG4x3/66rZIqEhEREYWccrAv82S59hMREZG/TiGnHNTwtu9jtLefiIiI/HX61S0HTUO9KMk5hmFYy1xvGAYlOcdoGupVyZWJiIhUXwo55cDdzcSJxeMB00WDzqnlX+h+OSIiIpVIIaecnNm5kmOzR2A5fbxUu2EpwWQy4df0GqyGrq4SERGpLLoZYDk6s3Mlh377FXOdq3D3r4El9yTWghwiHhyFT1xLZm/Po3UrR1cpIiJSPWgkp7wZVgoPbCJ/21IKD2yi+Nh+Tn4/HoCpm0+zdr+usBIREakMCjnlIDQ0FG9v74uuz924kPxtS7Ea8PSX68nOv/Q9dUREROSv0+mqchATE8OOHTvIysq6YN2GDRt47LHHyFrwPokJHTh06gwvzNzIRw+0xGTSRGQREZGKopBTTmJiYoiJibmgvWXLlpw4cYKZM2fyaq/G9J+1hwVbMpjyazoPXl3XAZWKiIhUDybDqL6X/OTk5BAUFER2djaBgYEVdhzDMCgpKcHT05P/LtvD8Lnb8PJwY/Y/O9IsquKOKyIi4ors/f3WnJxKYDKZ8PT0BOCRTnG0DPekqMTKgC/XkVdY4uDqREREXJNCTiUbMmQIX798F35uJew5lsewOVscXZKIiIhLUsipZBEREVjP5LDvy2GYgBlrDzJr/UFHlyUiIuJyFHIqWf/+/bnrrrvI37cB66ZvAXh51mb2HMt1cGUiIiKuRSGnkplMJiZMmECDBg1Inz8ev9xD5BdZGDB1PYUlFkeXJyIi4jIUchwgKCiI6dOnY/byZMfkwXibSth6JIeR87Y7ujQRERGXoZDjIM2bN+e9997DknuCgzNGADDp5318tyXDwZWJiIi4Bt0M0IEee+wx1q1bx7XXXsvugDgmLNvL8zM2El87iNrBPo4uT0REpEpTyHEgk8nEuHHjACgqsbJq7wk2HMxm4JfrSXn8ajzcNdAmIiJypfQr6iS8PNwYemM0ZjeDNftP8u73Ox1dkoiISJWmkOMkMjMzufXaqzk0+20APlyym+W/XfjATxEREbGPQo6TCA8P5+abbyZv2zJKdizBMOCZ1DSOnS50dGkiIiJVkkKOExk7dizx8fEc/vZ9PPOzyMotJHlaGlZrtX2GqoiIyBVTyHEivr6+TJ8+HV+zB/unDsEdC8t+y2Lc0t2OLk1ERKTKUchxMk2aNGH8+PEUHz/A0fkfADBq4U7W7j/p4MpERESqFoUcJ3Tffffxj3/8g9yNC/E9thWL1eDpL9eTnV/s6NJERESqDIUcJzVmzBheeuklFr/1KHVr+nLo1Bmen7kBw9D8HBEREXso5Dgpb29v3njjDSJDa/D+vS3wdDfx3ZZMpvyy39GliYiIVAkKOVVAQu0g2vtkAvD63G1sOZzt4IpEREScn0JOFfDJJ5/w+ZC/U7x/PUUlVp6aup68whJHlyUiIuLUFHKqgAcffJA2bdqQ8fU7uBXmsCcrj6Ffb3F0WSIiIk5NIacKMJvNTJs2jUCzG4dnvAGGwcx1B/lq3UFHlyYiIuK0FHKqiNjYWCZPnkzhwS2cWv4FAK/M3szuY7kOrkxERMQ5OVXI+eCDD4iNjcXb25t27dqxatWqS/Y/deoU/fv3JzIyErPZTKNGjZg3b14lVVv5brvtNp599lmyV06j+NAW8ossPDV1PQXFFkeXJiIi4nScJuSkpqaSnJzMsGHDWLduHUlJSXTr1o2jR4+W2b+oqIgbb7yRffv2MWPGDHbs2MGECROoXbt2JVdeuUaOHEn7q9uR9c1/8HO3svVIDiPnbXN0WSIiIk7HZDjJ3eXatWtHmzZtGDt2LABWq5Xo6GieeuopXnzxxQv6jxs3jnfeeYft27fj6el5RcfMyckhKCiI7OxsAgMD/1L9lenAgQMcPHiQgpD69Ju4GoBxD7Sie3yEgysTERGpePb+fjvFSE5RURFr166la9eutjY3Nze6du3KypUry9xmzpw5tG/fnv79+xMeHk58fDwjRozAYrn4qZvCwkJycnJKvaqi6Oho2rdvz3WNw3j8mnoAPD9jAwdP5ju4MhEREefhFCEnKysLi8VCeHh4qfbw8HAyMjLK3GbPnj3MmDEDi8XCvHnzGDJkCKNGjWL48OEXPc7IkSMJCgqyvaKjo8v1fThC96hiPLIPklNQwsCUNIotVkeXJCIi4hScIuRcCavVSlhYGOPHj6dVq1b06dOHl19+mXHjxl10m8GDB5OdnW17HThwoBIrrhgvD36B/V8Og+IC1u4/ybuLdjq6JBEREafgFCEnNDQUd3d3MjMzS7VnZmYSEVH2PJPIyEgaNWqEu7u7ra1p06ZkZGRQVFRU5jZms5nAwMBSr6ruk08+oYaXlWNz3wXgo592s+y3Yw6uSkRExPGcIuR4eXnRqlUrFi9ebGuzWq0sXryY9u3bl7lNx44d2bVrF1br+dMzO3fuJDIyEi8vrwqv2VlERUUxdepUzuz8mdPr52MYMCh1A8dOFzq6NBEREYdyipADkJyczIQJE5g8eTLbtm3jySefJC8vj379+gHw0EMPMXjwYFv/J598khMnTjBw4EB27tzJ3LlzGTFiBP3793fUW3CYrl27MmzYME7+MIGSrHSycgtJnpaG1eoUF86JiIg4hIejCzinT58+HDt2jKFDh5KRkUHz5s1ZsGCBbTJyeno6bm7nM1l0dDTfffcdgwYNIjExkdq1azNw4EBeeOEFR70Fh3rllVdYvnw5P80eSVS//2PZb1mMW7qbf17bwNGliYiIOITT3CfHEarqfXIuJjMzkxYtWhDV6U6y6nXD3c3EtH9cTau6IY4uTUREpNxUqfvkSPkIDw9n+fLl/DJ1NLc3j8JiNXj6yzRO5Zc9EVtERMSVKeS4mHr16uHh4cHwXvHUrenLoVNneH7GRqrxgJ2IiFRTCjkuyhMLtfd9B9YSFm7N5PNf9ju6JBERkUqlkOOi8vLyWDbnC0788CkAw7/dxpbD2Q6uSkREpPIo5LiomjVrkpqaypkN88j/7VeKLFaemrqevMISR5cmIiJSKRRyXFiHDh148803OT5vDJbTWezJymPI15sdXZaIiEilUMhxccnJyfTsdj3H5rwDhpWv1h1i5tqDji5LRESkwinkuDiTycTEiROJ9Mjj1PKpAAz5ejO7j+U6uDIREZGKpZBTDdSoUYNp06bhvWcJjYIM8ossDJi6noJii6NLExERqTAKOdVEmzZt2LtnD5/370pNPy+2HclhxLxtji5LRESkwijkVCM+Pj6EB3oz6u4kAD5buZ8Fm484uCoREZGKoZBTDVkPbaZow1wAnp+xkYMn8x1ckYiISPlTyKmGatSowYmfJlN4eDs5BSU8/eV6ii1WR5clIiJSrhRyqqEWLVrw3ph3yZrzDtaCXNaln2L0op2OLktERKRcKeRUU4899hh397iB4wveB+CjJbtZuvOYg6sSEREpPwo51ZTJZOLjjz8mxnSc0+vnATAoNY2jpwscXJmIiEj5UMipxvz9/Zk+fTpnfv6CoqN7OZ5XRHLqBqxWw9GliYiI/GUKOdVcfHw8H74/hi5eu/HxdGP5riw++mm3o8sSERH5y0yGYVTb/23PyckhKCiI7OxsAgMDHV2Ow01bc4DnZ2zE3c1E6uNX0zo2xNEliYiIXMDe32+N5IjNXa3qcHvzKCxWg6e+XM+p/CJHlyQiInLFFHLExmQysX/GmxSfOMyR7AKen7GRajzQJyIiVZxCjpRyV+/byJrzFkZJMQu3ZvLZyv2OLklEROSKKORIKffffz8P334DJ5d8CsDwb7ey+VC2g6sSERG5fAo5coExY8YQV5JO/m+/UGw1GDB1HbmFJY4uS0RE5LL8pZBz5swZDh06dEH7li1b/spuxcF8fHyYMX06hUs/oSTnGPuO5zN09mZHlyUiInJZrjjkzJgxg4YNG9KjRw8SExP59ddfbesefPDBcilOHKdhw4b898P/I2vOOxhWK1+tP8SMtQcdXZaIiIjdrjjkDB8+nLVr15KWlsbEiRN55JFHmDp1KoCuyHERd911F8P6P8D9SUEADJm9mV1Hcx1clYiIiH08rnTD4uJiwsPDAWjVqhVLly6ld+/e7Nq1C5PJVG4FimMNHjwYi9VgX/6v/Lz7OAOmrmN2/454e7o7ujQREZFLuuKRnLCwMDZu3GhbDgkJYdGiRWzbtq1Uu1R97m4mxvRpTqDZxPaM07wxd5ujSxIREflTdoecLVu2sGvXLtvy559/bhvJOcfLy4svv/ySn376qfwqFKdw5mQme798DYDPf9nP/E1HHFyRiIjIpdkdcpKTk/nwww9ty3Xq1GHNmjXcf//9DBo0iH379tnWdezYsVyLFMerW7cu/7itM9m/zADguelpHDiR7+CqRERELs7uB3RGREQwc+ZMW4DZtm0bSUlJhIWFUVhYiMlkIi0tjaioqAotuDzpAZ2Xp7i4mGuuvY599Xpjrt2E5nWCmP5kBzzddbslERGpPOX+gM7s7Gyio6Nty5999hn16tVj//79HDx4kKSkJN58882/VrU4NU9PT6alfEnJsglYC3JJO5jNqIU7HV2WiIhImewOOXXq1OHIkfPzMBYvXsxdd92Fu7s7ZrOZwYMHs3DhwgopUpxHdHQ0n330LsfnvwfAuJ92s3TnMQdXJSIiciG7Q07Xrl0ZPXo0APv372fdunXcdNNNtvX169fnwIED5V+hOJ2bb76Zp/92DafXzQUgeVoaR08XOLgqERGR0uyek3Po0CFatGiBv78/BQUFmM1m9uzZY7snzqpVq+jZsyeZmZkVWnB50pycK1dSUsJnX3zJtJN12Z5xmo4NavLZ39vh7qZ7JImISMUq9zk5tWvXZvXq1fTu3Zubb76Zr776qtRN/3744QcaNWr016qWKsPDw4O/932Qsfe1xMfTnRW7jjPup92OLktERMTmsu54XLduXUaNGlXmuq1bt3LnnXeWS1FSdTQI8+flmxvwypwdjPpuO+3iQmgdG+LoskRERP7aU8j/6LPPPmPgwIHltTupQo6vnkfu5h+wYuLJz1dxKr/I0SWJiIiUX8iR6mvAgP4kFm2l+MQhjuVZGJSyTg9pFRERh1PIkb/M3d2dLz+fhLHivxglxfy48ziTf97n6LJERKSaU8iRchEeHs6XH77DqSWfAvD6t1vYfCjbwVWJiEh1ppAj5aZLly688Lf25O9cicUw8dikX8gtLHF0WSIiUk0p5Ei5evHFF2hRvAVr7nGOnC7hlVmbND9HREQcQiFHypWbmxtfTJzA+Ievxt3NxOy0w8xYe9DRZYmISDWkkCPlLjQ0lO6tG5F849mbQw79ejO7jp52cFUiIlLdKORIhXmiS30aBlg4U2zlkU9+pqDY4uiSRESkGlHIkQrjZgLPtVOx5J1kf3YJQ79Kc3RJIiJSjSjkSIUxmUxM+e+HuK+eAsC09RnM23TYwVWJiEh1oZAjFSokJIRp7/2b06u+AmDQ1DUcOJHv4KpERKQ6UMiRCte2bVteuT2JgkPbKDTc+fuEZRRbrI4uS0REXJxCjlSKZ55+ijbFm7EU5PLbyRKGf73B0SWJiIiLU8iRSmEymZjy8f/huS4VgMmrDvPTzmMOrkpERFyZQo5UmqCgIDbMncyDV9cFIDk1jaM5BQ6uSkREXJVCjlSqgIAAXu7RlKaRgRzPK+Ifk1ZiseqxDyIiUv4UcqTSeXu689pNMZgsRaw/nM9/5m50dEkiIuKCFHLEIZrWCcE9bQYAHy1P59c9WQ6uSEREXI1CjjhEQEAAs0a/yJltP4HJjb9PWMbJvCJHlyUiIi5EIUccJiEhgX/f1oziE4fIM7z4+8c/YBianyMiIuVDIUcc6h+PPMzV1i0YJcWsP2rh/e82ObokERFxEQo54nCf/d8IzNvnAvDuD/vYfCjbwRWJiIgrUMgRh/Pz82P224MwH9uO4ebOgKnryC0scXRZIiJSxSnkiFNo1qwZv7w3gNrBPuw7ns/LszZpfo6IiPwlThVyPvjgA2JjY/H29qZdu3asWrXKru1SUlIwmUz06tWrYguUClXDz8x79zbH3c3E12mHmbB4s6NLEhGRKsxpQk5qairJyckMGzaMdevWkZSURLdu3Th69Oglt9u3bx/PPfccnTt3rqRKpSK1qhvCtTVzARj53S62Hz7l2IJERKTKcpqQM3r0aB577DH69etHs2bNGDduHL6+vnz66acX3cZisXD//ffz2muvUa9evUqsVirS8z2bU3xgE4a7F/e8t5CCYoujSxIRkSrIKUJOUVERa9eupWvXrrY2Nzc3unbtysqVKy+63b///W/CwsJ45JFH7DpOYWEhOTk5pV7ifBo3asS/e9THkneSU/jx2IcLHF2SiIhUQU4RcrKysrBYLISHh5dqDw8PJyMjo8xtli9fzieffMKECRPsPs7IkSMJCgqyvaKjo/9S3VJxHrnvLjq67QJg2RGY/IPunyMiIpfHKULO5Tp9+jQPPvggEyZMIDQ01O7tBg8eTHZ2tu114MCBCqxS/qrJbw/GZ98yAF6dt5Pdmbp/joiI2M/D0QUAhIaG4u7uTmZmZqn2zMxMIiIiLui/e/du9u3bR8+ePW1tVqsVAA8PD3bs2EH9+vUv2M5sNmM2m8u5eqkoZrOZOcMf57rhX+MR0Yi+437kx1duw9O9SmZzERGpZE7xa+Hl5UWrVq1YvHixrc1qtbJ48WLat29/Qf8mTZqwadMm0tLSbK/bbruN6667jrS0NJ2GciENG9Tj9VvqYTZZOHjGk/98t8PRJYmISBXhFCM5AMnJyfTt25fWrVvTtm1bxowZQ15eHv369QPgoYceonbt2owcORJvb2/i4+NLbR8cHAxwQbtUff3uvp3IZhk8MWUtHy/dQ/v6Nbm2cZijyxIRESfnFCM5AH369OE///kPQ4cOpXnz5qSlpbFgwQLbZOT09HSOHDni4CrFUbrHR/BQ+7oA/GPiCtKP6co4ERG5NJNRje+dn5OTQ1BQENnZ2QQGBjq6HPkTBUUlxA+aRElAJDUtx1n11oO4u5kcXZaIiFQye3+/nWYkR+TPeHt58PJ1kViLznDcvSb/HDvb0SWJiIgTU8iRKqXfnT3o4HX20v8FBz34atkGB1ckIiLOSiFHqpwpw5/C7+gmTG5uPDtjM4ePa36OiIhcSCFHqhwPDw9mD7kPa/YRDJ9gbn9jOtV4apmIiFyEQo5USQ3jYhh2YzRGSRHHvCIYu3CLo0sSEREno5AjVdYjf+vGbTFnn1D+3tJ0Nh3UYx9EROQ8hRyp0t576g66XRVOscVgwJfrOF1Q7OiSRETESSjkSJVmMpl4+44kagf7sP94Pt1enqT5OSIiAijkiAsI8vXk2Q4hGFYLhz2jGPh/KY4uSUREnIBCjriEv13TnLbeZ59i//UBM3OXr3NwRSIi4mgKOeIyvhz6d3xz0jF5mhnwxTqOnTjl6JJERMSBFHLEZXh4uDPrhV4YZ7IxgiK5dYjm54iIVGd6QKce0OlyPp69hBErT2MyudHRtINuXdpzssBKDW83moZ64e5mIjQ0lJiYGEeXKiIiV8De32+FHIUcl5Oenk6Lh4YScPVdGIYVk+n8gGVJzjFOLB6Pkb6eHTt2KOiIiFRBegq5VFtZWVmcObILwzBKBRwA94Ca1Or1EqaYFmRlZTmoQhERqQwKOeJyLFaDkBseLXPd2dBjEHLD41is1XYQU0SkWlDIEZezLasIj8BamEymMtebTG54BNZiW1ZRJVcmIiKVSSFHXM7JAqtd/SYtXMuO33ZXcDUiIuIoCjnicmp42/fPeo9XHDd+uI6eb8xg0dZMCkssFVyZiIhUJg9HFyBS3pqGelGScwz3gJoXTDwGMAwDo+gM7kYxbt5BbDoNj322hgCzB/E1rESWZJB8383UiYpwQPUiIlJeFHLE5bi7mTixeDy1er10wSXkhmEFTGTNe5elU98nx7Mmy/fnM2/TETJyCliZARDGjHeWEpS7n1viIxh0781EhIU66u2IiMgV0n1ydJ8cl5Oenk7jxo0xxbQg5IbH8QisZVt3sfvkWK0G69JPMmLKAtKOm7Caz/97sBbkEXLmAD2b1+aVR+/C7Ole6e9JRETO080A7aCQ47rS09PJysrCYjXYllV0WXc8tloNvvllCxMWrGVLjheG9/l/G4HeHtx0VQQ9EiJpWzcIPx9zZb0lERH5nUKOHRRy5M9YrQZfLUvjk4XrOWiqxemSP5z6KswjrPgIvVvFMuCuGwn093VgpSIi1YdCjh0UcuRyWKwGa/adYO6mI3y1ei+5fwg81oJcoqxHuaNtPfrf2RVfb43wiIhUFIUcOyjkyJUqsViZsuBnPluyhd1FAZh8gmzrjIJcujQI5tGbWtK+fk083XWnBhGR8qSQYweFHCkPxSUWJn67jC+WbWVfSXCpwBPs60nzmtAuypNHbu2Ml6cuaBQR+asUcuygkCPlrbComJQf1rDzjD8LNmdwPO/8oyOMglzivLJ58JqreKj71Xh66CotEZEroZBjB4UcqUglFiu/7DnO82NTOEhN3P4wwkNhLg3MufS9PoF7b2iNh05piYjYTSHHDgo5Ully8/L5cMYiZqzaS4ZHOG4+5/+91fTzolt8BLcmRNI2LkSBR0TkTyjk2EEhRxwh+/Rp3k/5jllr93OmZiPyLedDjakwl6YBhTxyUwt6dYjH3a3sJ6mLiFRnCjl2UMgRRyv+/ZTW3I1HmL1mLwXG+YnJpqJc4gOLeax7K3q0baLAIyLyO4UcOyjkiDPJOHqM/0uZz9yNRzjpF4P7H05puRXl0iOxNvdf05Q2sSEKPCJSrSnk2EEhR5zVocNHGPPlfOZtziQnMA53nwDbuloBZtpEeNK7VQzXJ8Yq8IhItaOQYweFHKkK9qcfJGXJerJ8YvhuSwY5BSW2dR7FebQKc+OJHm3p0iwaNwUeEakGFHLsoJAjVU1RiZXvNx3gmdGfUVCzEW7e/rZ1HiV5tI3w5IkebenUOEqBR0RclkKOHRRypCrbvHU7/5e6gB9+y6Y4rEmpwBMeaObm+EhuTYykZUwNBR4RcSkKOXZQyBFXYBgG6zds4r1pC/lpz2nM9dtQYDkfarws+XSM9uUft7Shbb1aCjwiUuUp5NhBIUdcjWEYFBSX8PPuE8zdeIRv1u+nmPOXpZstZ+gc68/jN7emdWxNBR4RqZIUcuygkCOu7tfVaxk1dR4rDxZiqpOEm9nXts7beoY72sRxR7t6tIgOxmRS4BGRqkEhxw4KOVJdWCwWFv/4E+Nm/8SvGSW4RzcvFXhqB/twTVwAd13dgBYxCjwi4twUcuygkCPVUXFxMd8tWsyyXccpCGvG91szySuy2Nb7UsCNjUN5+IZEmmuER0SckEKOHRRyRKCg2ML05Vv41/tf4hXbCjcvH9s6Pwrp3iyMh66LJ7FOkAKPiDgFhRw7KOSInJefn8/sb+byybyVbMkx4xVXOvDUqeFDj4RIeiRGklBbgUdEHEchxw4KOSJlO336NDNmz2Hi/F/Zke9LULPOFFnPrw90K6JHUhT3dWxCfO1ABR4RqVQKOXZQyBH5cydPnsTT249lu47z7aYjLNhwEIvJ3bY+yL2Y25rXoU+HhlwVpcAjIhVPIccOCjkil2/SlKl88NUS9llr4lO/NW6e3rZ1wR4l3N2+Abc1r6PAIyIVRiHHDgo5Ilfu8OHDTE2dwZQf0jjkHn5B4Imt6Uv3q8Lp2bw2zSIVeESk/Cjk2EEhR6R87Nu3jy9SZ3DQGgzRLfhh+1EKS85P4qnpZeFvbeLo3aouTSMDFHhE5C9RyLGDQo5IxcgrLOGjOct5a+pCfOq1xs3TbFtXy9vK39rE0atlDE0iFHhE5PIp5NhBIUekYm3ZsoUpKdOZtmIbp4Ma4FO/NSYPL9v6erX8bJelNw5X4BER+yjk2EEhR6RyGIbB+vXrmZI6g5krd3AmtClBTdpT/IfL0iN84c629ejZvA6Nwv0VeETkohRy7KCQI1L5DMNg9erVNElozg/bjzF30xG+33IY4w+XpUf5mbizXT1uTapNo/AAB1YrIs5IIccOCjkizuHZF1/hiyUbsdZOwieuFSYPT9u62gHu3NW2Hj0SI2mowCMiKOTYRSFHxHmUlJTwww8/8MW0mczfeAiiW+IT17JU4GkU7k+PhCh6JEbQIEyBR6S6Usixg0KOiHMqKipi4cKFTEmdiW/DqzHFtGTpb8cotpz/c1U32JM72sRxS0IkDcL8HVitiFQ2hRw7KOSIVB3ZZ4p5a8o8/vvdOnziWmByPz/CE1fDi7+1juWWxEjq11LgEXF1Cjl2UMgRqVpOnz7NN998wxfTZ7F8bw7mhu3xiW2Byd3D1qdJRAC3JkZyS0Ik9RR4RFySQo4dFHJEqq5Tp04xe/Zsvpg+i18PFeDTsAP+DVrzhzNaNAj1pnerutySEElcqJ/jihWRcqWQYweFHBHXcOzYMRYvXkz32/7Gwi2ZzN10hKU7MjFMbrY+jWr5cnvLaHokRBKrwCNSpSnk2EEhR8Q1GYbB9d17siajCN/GnfCObY7J7fx9eJqE+3Fbizr0SIikbk0FHpGqRiHHDgo5Iq4tPT2dadOmMXXmHH7L98G3SSe86yaVCjzxtQPPXpaeEElMTV8HVisi9lLIsYNCjkj1sWvXLlJTU5k6cw5tej9KSWQiK/ccx2I9/yfwqsgAejavTY+ESKJDFHhEnFWVDDkffPAB77zzDhkZGSQlJfH+++/Ttm3bMvtOmDCBzz77jM2bNwPQqlUrRowYcdH+ZVHIEameLBYL7u7uHM8t5LVPZpO6chfeMQmlR3iiAumZFMUtCjwiTsfe32+3i66pZKmpqSQnJzNs2DDWrVtHUlIS3bp14+jRo2X2X7JkCffeey8//vgjK1euJDo6mptuuolDhw5VcuUiUtW4u58NMzX9zfTv3px+sacxff0SxxeM5cy+DRhWC5sP5zBy/nY6v/0jt49dzviluzl4Mt/BlYvI5XCakZx27drRpk0bxo4dC4DVaiU6OpqnnnqKF1988U+3t1gs1KhRg7Fjx/LQQw/ZdUyN5IjIOYZhsGrVKlJSUpg2Zz45gfXwa9IJ39gk/nBGi6Q6QdyaGMXNCRHUqaERHhFHqFKnq4qKivD19WXGjBn06tXL1t63b19OnTrF119//af7OH36NGFhYUyfPp1bb721zD6FhYUUFhbalnNycoiOjlbIEZFSrFYry5cvZ/Xq1Tz0+AAWbMlg7sbD/LL7OJhMtn7No88FnkhqB/s4sGKR6sXekONx0TWVKCsrC4vFQnh4eKn28PBwtm/fbtc+XnjhBaKioujatetF+4wcOZLXXnvtL9UqIq7Pzc2Na665hmuuuQaAB6+uy51JYdRtksCZmo3xa9IZc/RVpB3IJu1ANsPnbqNFdDA9fr/TcpQCj4hTcJo5OX/Fm2++SUpKCrNmzcLb2/ui/QYPHkx2drbtdeDAgUqsUkSqMh8fHw7+tpWU1//JTR5bOTW5P8cXfkRB+iYMw8r6A6cYPncbHd78gb99uIJPlu/lSPYZR5ctUq05xUhOaGgo7u7uZGZmlmrPzMwkIiLiktv+5z//4c033+T7778nMTHxkn3NZjNms/kv1ysi1ZOXlxc9evSgR48enDlzhgULFpCSksLcT97jvhfe4VRAPVbvP8G69FOsSz/F699upXXdGtyScHaEJyLo4v8TJiLlzynm5MDZicdt27bl/fffB86eE4+JiWHAgAEXnXj89ttv88Ybb/Ddd99x9dVXX/YxNfFYRMpDbm4ubm5u+Pr6kplTwOAPpzFvUwbmOk0x/eHREm1izwee8EAFHpErVaXm5AAkJyfTt29fWrduTdu2bRkzZgx5eXn069cPgIceeojatWszcuRIAN566y2GDh3K1KlTiY2NJSMjAwB/f3/8/fXkYRGpPH/8mxMe6M3dzcPIWJrC4nmjMNdvd/ZOy3WuYvW+k6zed5J/f7OVNrEh3JIQwc0KPCIVxmlGcgDGjh1ruxlg8+bNee+992jXrh0A1157LbGxsUyaNAmA2NhY9u/ff8E+hg0bxquvvmrX8TSSIyIVKSsri6+++orU1FSWrdmEd8P2+DftjFftprY+JhO0iQ2hR0IkN8dHEKbAI/KnqtQl5I6ikCMilSUjI4MZM2Zw+PBhBjw/hPmbz16Wvi79lK2PCWgbF0KPxEi6x0cQFqDAI1IWhRw7KOSIiCMdPnyYmCZJ+DbuePay9NpNbOtMQLt6Z0d4usdHUitAF02InKOQYweFHBFxtD179pCamkpqaiqb9xzGt0lH/Jp0whx1PvC4maBdXE1uSYyk+1URCjxS7Snk2EEhR0ScybZt20hNTSUlJYVdR07w+Osfccg9gg0HTtn6uJng6no1uSXh7CmtUH8FHql+FHLsoJAjIs7IMAw2btxIgwYN8PPz48CJfJ77vy9YujcXc1QjWz83E7Sv/3vguSqCmgo8Uk1UuUvIRUTkLJPJRFJSkm05OsSXDiH5/PzpaA6dLsa3cUd8m3TCHNmIFbuOs2LXcYZ+vYX2v4/wdLsqXIFHBI3kaCRHRKoMq9XKzz//TGpqKtOmTeNEoQnfxp3wv+oaPMPq2/q5u5noUP9c4IkgxM/LgVWLlD+drrKDQo6IVFUWi4WffvqJlJQUvLy8+NdrbzFv09nL0jcfzrH1Oxd4evweeGoo8IgLUMixg0KOiLiaTZs20fKabudPaUU0sK1zdzPRsUEoPRIiuKmZAo9UXQo5dlDIERFXU1BQwHfffUdKSgpz5syhyCsI3yad8GvSCa/w86e0PM4FnsRIujWLIMjX04FVi1wehRw7KOSIiCvLy8tj7ty5pKSkMG/ePCy+NXlyxHh+Kwxk25Hzp7Q83Ex0ahhKj4RIblLgkSpAIccOCjkiUl3k5OQwZ84c7r77bry8vNhzLJeB/5nE2mPgFRZn6+fpbqJTg1B6JEZxY7NwgnwUeMT5KOTYQSFHRKqzv//970yePBm34Ej8GnfCt0mnCwJP54a16JEQSVcFHnEiCjl2UMgRkeouMzOTmTNnkpKSwvLly3GvURu/Jp0IjL8Wtxp1bP083U1c07AWPRLPBp5AbwUecRyFHDso5IiInHfw4EGmT59OamoqLVq0IPnVt5m76QjzNh5h59FcWz8vdzeuaXR20nLXpuEEKPBIJVPIsYNCjohI2UpKSvDwOHtT/BUrVnBdr/vwbdwJv6ad8awZbet3NvDU4tbESG5oGqbAI5VCIccOCjkiIn/u4MGDfPLJJ6SkpLB9+3Y8Q2PwbdIJ/6bX4BFy/pSWl4cbXWyBJxx/s54cJBVDIccOCjkiIvYzDINNmzaRmppKamoqu3fvxjO0Lv98879sPOXBnmN5tr5eHm5c2+jsHB4FHilvCjl2UMgREbkyhmGwdu1avvnmG4YNG4bJZGJH5mkGvPUpO8744x4caetr9nDj2sa16JEYxQ1NwvBT4JG/SCHHDgo5IiLlxzAMmjdvzsaNG/GsFYtfk04EXHUtbkERtj5mDzeuaxxGj8RIrlfgkSukkGMHhRwRkfJlsVhYunQpKSkpzJw5k+PHj+NZKw6/Jp0IaXEjFp8QW19vz9KBx9dLgUfso5BjB4UcEZGKU1xczOLFi0lJSWHWrFn844kn6DvwJeZuPMK3Gw+TfuKMra+3pxvXNwmjR0IU1zWppcAjl6SQYweFHBGRylFQUEBBQQHBwcEAzJs3j179nsKvSWeCE6/D8Au19fXxdD8beBIjua5xGD5e7g6qWpyVvb/fisoiIlLhvL298fb2ti2HhoZyW+cWzJ2byr6lk/EKr49vk04EJ97AGd8Q5m46wtxNR84GnqZh3JoQybUKPHKZNJKjkRwREYc5ffo0c+bMISUlhe+++47i4mK8wuvT/61PWZVRwsGT509p+XqdHeG5NfFs4PH2VOCprnS6yg4KOSIizuPkyZPMmjWLn3/+mQkTJgCw6VA2/xz5X454RmLxDrb19fVy54am4fRIiOTaxrUUeKoZhRw7KOSIiDi3oqIiatWqRU5ODl4RDfFr2pmQpK6UmM//zfY7F3gSI+nSSIGnOlDIsYNCjoiI8zt8+DDTp08nJSWFX375BQCvyEb4N+1MrVbdKXDzsfX1N3twQ9MweiREco0Cj8tSyLGDQo6ISNWyb98+pk2bRmpqKuvWrePf/36dHn37M+/3y9Izcgptff3NHnRtGkaPxCg6NwxV4HEhCjl2UMgREam6du7cSVBQEOHh4QB8mZLKw8++SlD8tQQnXEeRh5+tb4DZg67Nzs7h6dwoFLOHAk9VpkvIRUTEpTVq1KjUsgmDGN8Sdi0cx7GFH2Ou3Zig+OsJir+W0/gya/0hZq0/RIDZgxubnZ3D06mhAo8r00iORnJERFyGYRisX7+elJQUUlNTSU9PB0yYazeh/5ufsDw9j8w/nNIK8D4beG5NjKRTg1p4ebg5rnixm05X2UEhR0TEdRmGwS+//EJqair79+9n1qxZWK0Ga9NPMmj05xzzrkOh2/kbFAZ6e3BjswhuTYykY4NQBR4nppBjB4UcEZHq5+TJk4SHh1NcXIK5TjNqtboJv8adKDCZbX0CvT246aoIeiRG0rG+Ao+zUcixg0KOiEj1U1JSwo8//khKSgpfffUVp06dApMb5tpNiWhzM35NO5NnOT9PJ8jHk5t+n8PTsUEonu4KPI6mkGMHhRwRkeqtqKiIhQsXkpKSwtdff01ubi5j/u892t16H/M2HWHuxsMczyu29Q/2PRd4ouhQv6YCj4Mo5NhBIUdERM45c+YM8+bNo3PnzoSFhQHw8fgJDBz+PtEdeuIe14YCw9PWP9jXk27Nzp7Saq/AU6l0CbmIiMhl8PHx4Y477ijVdvjQQawZ29mVuunsKa3oq6jbqTduMS05lQ+paw6QuuYANXw96fb7HJ729WriocDjFDSSo5EcERG5hFOnTjF79mxSUlL4/vvvsVgsYHLDNzaJJ0eM58ffTnI8r8jWP8TP62zgSYjk6nohCjwVQKer7KCQIyIil+PYsWN89dVXpKam4uXlxYIFCyixWPl17wle/ngmmZ6RnLGen7Qc4udF9/izgaddnAJPeVHIsYNCjoiIXKmSkhI8PM7O+jh8+DB16tTBwIRvbHMaXH8XRWFXccZ6PtTU/GPgqVcTdzeTo0qv8hRy7KCQIyIi5eHEiRN8+umnpKamsmbNmrONJjf867WkwfV3UxzejNzzF2kR6n828NySEEm7OAWey6WQYweFHBERKW+7du0iNTWV1NRUNm3aBMAnEyfRoMMtzN14hAVbjpB9psTWP9TfzM2/B562cSEKPHZQyLGDQo6IiFSkrVu3kpqaSnJyMkFBQQC8/Z9RjPz0Kxpefze5wfXJLzkfakL9zdyScDbwtIlV4LkYhRw7KOSIiEhl69WrF19//fXZBTd3ajbrQIPr7yYnII788wM81Aowc8vvIzytFXhKUcixg0KOiIhUNsMwWLVqFSkpKUybNo3Dhw+fXeHmQXTrrtyVPIJF246SU3A+8YQFmLklIfJs4KlbA7dqHngUcuygkCMiIo5ktVpZvnw5KSkpzJgxg86dOzNz5kyKSqys2JXFyC8WctAIIb/k/E91eKCZm+Mj6ZEYSauY6hl4FHLsoJAjIiLOoqSkhJMnT1KrVi0Adu7cSePGjcHdg+jWNxJ3zR0cM0eSX3z+Zzsi0JubE85elt6yGgUehRw7KOSIiIiz2rFjB8OHD+frr7/m9OnTZxvdPYhrfwt1O/XiqEcEecVWW/+IQG9uSYikR2IELaJdO/Ao5NhBIUdERJzdmTNnmD9/PqmpqXzzzTecOXMGgGkzvqJG0w7M3XSERVszyC202LaJDPK2zeFpER3scoFHIccOCjkiIlKV5Obm8u233zJr1iwmT56Mt7c3AK8MfZXZv+4ksl0PDhkh5P9hhCfqXOBJPBt4TKaqH3gUcuygkCMiIq4gPj6eLVu2nF1w9ySp2z2Et7mZdEswZ/4QeGoH+9juw9O8CgcehRw7KOSIiIgrOH78OF999RUpKSksWbIEq/VssHHzNNP57n+Q2ONhFm/LJK/o/Cmt2sE+9EiMpEdCJIl1gqpU4FHIsYNCjoiIuJqMjAxmzJhBSkoKK1asoF+/fnz66acUFFtYsuMoH37zC7/le5ca4alTw4ceCWcvS0+o7fyBRyHHDgo5IiLiyg4cOEBxcTH16tUDYM2aNbRp0wZPHz/a9epHcML1/JbvTcEfAk90iA+3JERya0IU8bUDnTLwKOTYQSFHRESqk8WLF5OcnMzGjRttbWa/ADrc8Sj+TbvwW55XqRGemBDfs4EnMZKropwn8Cjk2EEhR0REqqNt27aRmppKSkoKO3bssLXP/W4RRF7F3I1HWLw9s9QIT0yIr20Oz58FHovVYNXeExw9XUBYgHe5P11dIccOCjkiIlKdGYbBxo0bSU1N5ccff2TZsmV4eHgA8NLQ10jLLMa7YXu2ZruVCjx1a/ra5vA0iywdeBZsPsJr32zlSHaBrS0yyJthPZvRPT6yXOpWyLGDQo6IiMiFrFYr0dHRtoeHhoRF0vGux/Gs147NJ6Cw5Hzgia15boQniv3H8/jnF+v432BxLgJ99EDLcgk6Cjl2UMgRERG5kGEYrFixgtTUVKZNm8bRo0dt68Jrx9Ct37P4NOrIjzuOlgo87iawXCJVhPl7svKlG//yqSt7f7/d/tJRRERExOWYTCY6derE+++/z6FDh/j+++957LHHCAkJIfNQOjVO72Hcg61YO+RGxvRJom2UGQ+3SwccgKO5xXz767bKeRMo5IiIiMgleHh4cMMNNzB+/HiOHDnC3LlzefLJJwHwN3sQeHIn0wfeSOHySXbtb1/myQqstjSPSjuSiIiIVGleXl7ccsstpdr27NmDr68vx/btIKLDn++jhnflja9oJEdERESu2KOPPsrRo0cZ+uR9lOQcwzCsZfYzDCslOcdoGupVabU5Vcj54IMPiI2Nxdvbm3bt2rFq1apL9p8+fTpNmjTB29ubhIQE5s2bV0mVioiIyDl+fn5073YTJxaPB0wXBJ2zyyZOLB5frvfL+TNOE3JSU1NJTk5m2LBhrFu3jqSkJLp161ZqRvcf/fzzz9x777088sgjrF+/nl69etGrVy82b95cyZWLiIgIwJmdKzk2ewSW08dLtVtOH+fY7BGc2bmyUutxmkvI27VrR5s2bRg7dixw/hr9p556ihdffPGC/n369CEvL49vv/3W1nb11VfTvHlzxo0bZ9cxdQm5iIhI+Vi3bh2tWrU6u2Byw1znKtz9a2DJPUnhwS3w++jO2rVradmy5V86lr2/304x8bioqIi1a9cyePBgW5ubmxtdu3Zl5cqyU9/KlStJTk4u1datWzdmz5590eMUFhZSWFhoW87OzgbOflgiIiJy5XJzc88vGFYKD2y6aL+/+rt7bvs/G6dxipCTlZWFxWIhPDy8VHt4eDjbt28vc5uMjIwy+2dkZFz0OCNHjuS11167oD06OvoKqhYREZHL1aVLl3Lb1+nTpwkKCrroeqcIOZVl8ODBpUZ/rFYrJ06coGbNmk7zZFVnkpOTQ3R0NAcOHNDpPCeh78S56PtwLvo+nEtFfh+GYXD69GmioqIu2c8pQk5oaCju7u5kZmaWas/MzCQiIqLMbSIiIi6rP4DZbMZsNpdqCw4OvrKiq5HAwED9wXAy+k6ci74P56Lvw7lU1PdxqRGcc5zi6iovLy9atWrF4sWLbW1Wq5XFixfTvn37Mrdp3759qf4AixYtumh/ERERqV6cYiQHIDk5mb59+9K6dWvatm3LmDFjyMvLo1+/fgA89NBD1K5dm5EjRwIwcOBAunTpwqhRo+jRowcpKSmsWbOG8ePHO/JtiIiIiJNwmpDTp08fjh07xtChQ8nIyKB58+YsWLDANrk4PT0dN7fzA08dOnRg6tSpvPLKK7z00ks0bNiQ2bNnEx8f76i34HLMZjPDhg274BSfOI6+E+ei78O56PtwLs7wfTjNfXJEREREypNTzMkRERERKW8KOSIiIuKSFHJERETEJSnkiIiIiEtSyJELjBw5kjZt2hAQEEBYWBi9evVix44dji5Lfvfmm29iMpl45plnHF1KtXXo0CEeeOABatasiY+PDwkJCaxZs8bRZVVLFouFIUOGEBcXh4+PD/Xr1+f111//02caSflZunQpPXv2JCoqCpPJdMEzJA3DYOjQoURGRuLj40PXrl357bffKqU2hRy5wE8//UT//v355ZdfWLRoEcXFxdx0003k5eU5urRqb/Xq1Xz88cckJiY6upRq6+TJk3Ts2BFPT0/mz5/P1q1bGTVqFDVq1HB0adXSW2+9xUcffcTYsWPZtm0bb731Fm+//Tbvv/++o0urNvLy8khKSuKDDz4oc/3bb7/Ne++9x7hx4/j111/x8/OjW7duFBQUVHhtuoRc/tSxY8cICwvjp59+4pprrnF0OdVWbm4uLVu25MMPP2T48OE0b96cMWPGOLqsaufFF19kxYoVLFu2zNGlCHDrrbcSHh7OJ598Ymu744478PHxYcqUKQ6srHoymUzMmjWLXr16AWdHcaKionj22Wd57rnnAMjOziY8PJxJkyZxzz33VGg9GsmRP5WdnQ1ASEiIgyup3vr370+PHj3o2rWro0up1ubMmUPr1q256667CAsLo0WLFkyYMMHRZVVbHTp0YPHixezcuROADRs2sHz5cm6++WYHVyYAe/fuJSMjo9TfraCgINq1a8fKlSsr/PhOc8djcU5Wq5VnnnmGjh076m7SDpSSksK6detYvXq1o0up9vbs2cNHH31EcnIyL730EqtXr+bpp5/Gy8uLvn37Orq8aufFF18kJyeHJk2a4O7ujsVi4Y033uD+++93dGkCZGRkANieXnBOeHi4bV1FUsiRS+rfvz+bN29m+fLlji6l2jpw4AADBw5k0aJFeHt7O7qcas9qtdK6dWtGjBgBQIsWLdi8eTPjxo1TyHGAadOm8cUXXzB16lSuuuoq0tLSeOaZZ4iKitL3ITpdJRc3YMAAvv32W3788Ufq1Knj6HKqrbVr13L06FFatmyJh4cHHh4e/PTTT7z33nt4eHhgsVgcXWK1EhkZSbNmzUq1NW3alPT0dAdVVL3961//4sUXX+See+4hISGBBx98kEGDBtke5iyOFRERAUBmZmap9szMTNu6iqSQIxcwDIMBAwYwa9YsfvjhB+Li4hxdUrV2ww03sGnTJtLS0myv1q1bc//995OWloa7u7ujS6xWOnbseMEtFXbu3EndunUdVFH1lp+fX+rhzQDu7u5YrVYHVSR/FBcXR0REBIsXL7a15eTk8Ouvv9K+ffsKP75OV8kF+vfvz9SpU/n6668JCAiwnTcNCgrCx8fHwdVVPwEBARfMh/Lz86NmzZqaJ+UAgwYNokOHDowYMYK7776bVatWMX78eMaPH+/o0qqlnj178sYbbxATE8NVV13F+vXrGT16NH//+98dXVq1kZuby65du2zLe/fuJS0tjZCQEGJiYnjmmWcYPnw4DRs2JC4ujiFDhhAVFWW7AqtCGSL/AyjzNXHiREeXJr/r0qWLMXDgQEeXUW198803Rnx8vGE2m40mTZoY48ePd3RJ1VZOTo4xcOBAIyYmxvD29jbq1atnvPzyy0ZhYaGjS6s2fvzxxzJ/M/r27WsYhmFYrVZjyJAhRnh4uGE2m40bbrjB2LFjR6XUpvvkiIiIiEvSnBwRERFxSQo5IiIi4pIUckRERMQlKeSIiIiIS1LIEREREZekkCMiIiIuSSFHREREXJJCjoiIHcaPH090dDRubm6MGTOGV199lebNm9u9/b59+zCZTKSlpV20z5IlSzCZTJw6deqifUwmE7Nnz7b7uCLVmUKOSDVgMpku+Xr11VcdXaJTy8nJYcCAAbzwwgscOnSIxx9/nOeee67U83hExPno2VUi1cCRI0ds/52amsrQoUNLPWTS39/f9t+GYWCxWPDwqHp/Hiqq9vT0dIqLi+nRoweRkZG29j9+biLifDSSI1INRERE2F5BQUGYTCbb8vbt2wkICGD+/Pm0atUKs9nM8uXL2b17N7fffjvh4eH4+/vTpk0bvv/++1L7jY2NZcSIEfz9738nICCAmJiYUg+qLCoqYsCAAURGRuLt7U3dunUZOXIkAPfddx99+vQptb/i4mJCQ0P57LPPALBarYwcOZK4uDh8fHxISkpixowZtv7nTu/8b+0bNmzguuuuIyAggMDAQFq1asWaNWts2y1fvpzOnTvj4+NDdHQ0Tz/9NHl5eWV+dpMmTSIhIQGAevXqYTKZ2LdvX5mnq/773//StGlTvL29adKkCR9++OElv5d58+bRqFEjfHx8uO6669i3b98l+5+TlZVF79698fX1pWHDhsyZM8eu7USqnUp5QpaIOI2JEycaQUFBtuVzD9dLTEw0Fi5caOzatcs4fvy4kZaWZowbN87YtGmTsXPnTuOVV14xvL29jf3799u2rVu3rhESEmJ88MEHxm+//WaMHDnScHNzM7Zv324YhmG88847RnR0tLF06VJj3759xrJly4ypU6cahmEY3377reHj42OcPn3atr9vvvnG8PHxMXJycgzDMIzhw4cbTZo0MRYsWGDs3r3bmDhxomE2m40lS5ZcsvarrrrKeOCBB4xt27YZO3fuNKZNm2akpaUZhmEYu3btMvz8/Ix3333X2Llzp7FixQqjRYsWxsMPP1zm55Wfn298//33BmCsWrXKOHLkiFFSUmIMGzbMSEpKsvWbMmWKERkZacycOdPYs2ePMXPmTCMkJMSYNGmSYRiGsXfvXgMw1q9fbxiGYaSnpxtms9lITk42tm/fbkyZMsUIDw83AOPkyZMX/f4Ao06dOsbUqVON3377zXj66acNf39/4/jx43Z8+yLVi0KOSDVzsZAze/bsP932qquuMt5//33bct26dY0HHnjAtmy1Wo2wsDDjo48+MgzDMJ566inj+uuvN6xW6wX7Ki4uNkJDQ43PPvvM1nbvvfcaffr0MQzDMAoKCgxfX1/j559/LrXdI488Ytx7772XrD0gIMAWLv7XI488Yjz++OOl2pYtW2a4ubkZZ86cKXOb9evXG4Cxd+9eW9v/hpz69evbAtw5r7/+utG+fXvDMC4MOYMHDzaaNWtWqv8LL7xgV8h55ZVXbMu5ubkGYMyfP/+i24hUVzpdJSIAtG7dutRybm4uzz33HE2bNiU4OBh/f3+2bdtGenp6qX6JiYm2/z53Guzo0aMAPPzww6SlpdG4cWOefvppFi5caOvr4eHB3XffzRdffAFAXl4eX3/9Nffffz8Au3btIj8/nxtvvBF/f3/b67PPPmP37t2XrD05OZlHH32Url278uabb5bqv2HDBiZNmlRqn926dcNqtbJ3794r+uzy8vLYvXs3jzzySKn9Dh8+/IJaz9m2bRvt2rUr1da+fXu7jvfHz9zPz4/AwEDbZy4i51W9mYUiUiH8/PxKLT/33HMsWrSI//znPzRo0AAfHx/uvPNOioqKSvXz9PQstWwymbBarQC0bNmSvXv3Mn/+fL7//nvuvvtuunbtaptXc//999OlSxeOHj3KokWL8PHxoXv37sDZkAUwd+5cateuXeoYZrP5krW/+uqr3HfffcydO5f58+czbNgwUlJS6N27N7m5ufzjH//g6aefvuAziImJseuz+l/nap0wYcIFwcXd3f2K9nkpl/rMReQ8hRwRKdOKFSt4+OGH6d27N3D2h9zeibF/FBgYSJ8+fejTpw933nkn3bt358SJE4SEhNChQweio6NJTU1l/vz53HXXXbYf8GbNmmE2m0lPT6dLly6XfdxGjRrRqFEjBg0axL333svEiRPp3bs3LVu2ZOvWrTRo0OCy93kx4eHhREVFsWfPHttI1J9p2rTpBROGf/nll3KrSUQUckTkIho2bMhXX31Fz549MZlMDBky5LJHC0aPHk1kZCQtWrTAzc2N6dOnExERQXBwsK3Pfffdx7hx49i5cyc//vijrT0gIIDnnnuOQYMGYbVa6dSpE9nZ2axYsYLAwED69u1b5jHPnDnDv/71L+68807i4uI4ePAgq1ev5o477gDghRde4Oqrr2bAgAE8+uij+Pn5sXXrVhYtWsTYsWMv/4P63WuvvcbTTz9NUFAQ3bt3p7CwkDVr1nDy5EmSk5Mv6P/EE08watQo/vWvf/Hoo4+ydu1aJk2adMXHF5ELaU6OiJRp9OjR1KhRgw4dOtCzZ0+6detGy5YtL2sfAQEBvP3227Ru3Zo2bdqwb98+5s2bh5vb+T89999/P1u3bqV27dp07Nix1Pavv/46Q4YMYeTIkTRt2pTu3bszd+5c4uLiLnpMd3d3jh8/zkMPPUSjRo24++67ufnmm3nttdeAs/NZfvrpJ3bu3Ennzp1p0aIFQ4cOJSoq6rLe2/969NFH+e9//8vEiRNJSEigS5cuTJo06aK1xsTEMHPmTGbPnk1SUhLjxo1jxIgRf6kGESnNZBiG4egiRERERMqbRnJERETEJSnkiIiIiEtSyBERERGXpJAjIiIiLkkhR0RERFySQo6IiIi4JIUcERERcUkKOSIiIuKSFHJERETEJSnkiIiIiEtSyBERERGXpJAjIiIiLun/AeR+iwHxWX2NAAAAAElFTkSuQmCC", "text/plain": [ "
" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "basis = state0.basis\n", "\n", "S0_exact = basis.ent_entropy(state0.psi, density=False, alpha=2.0)[\"Sent_A\"]\n", "S1_exact = basis.ent_entropy(state1.psi, density=False, alpha=2.0)[\"Sent_A\"]\n", "S2_exact = basis.ent_entropy(state2.psi, density=False, alpha=2.0)[\"Sent_A\"]\n", "\n", "plt.plot(h, [S0_exact, S1_exact, S2_exact], 'ks--', label='Exact')\n", "plt.plot(h, [S0, S1, S2], 'o-', label='Sampled')\n", "plt.xlabel('Transverse field h')\n", "plt.ylabel(r'$S_2$')\n", "plt.ylim(0, 1)\n", "plt.legend()\n", "plt.show()" ] } ], "metadata": { "kernelspec": { "display_name": "quantax_env", "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.12.3" } }, "nbformat": 4, "nbformat_minor": 5 }