doc/SBRG.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "c1de376d-03ac-4eaf-a3a3-5b5c3040b627",
   "metadata": {},
   "outputs": [],
   "source": [
    "# will autoupdate any of the packages imported:\n",
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "15fb0499-affc-4fb1-be90-30f6a969c4fb",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy\n",
    "from context import *"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fe0a2582-9df2-41a6-a932-0fc4f24e1cd8",
   "metadata": {},
   "source": [
    "# Spectrum Bifurcation Renormalization Group (SBRG)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "37f1f5dd-fcc3-4d93-90ef-7dd6dcec7de4",
   "metadata": {},
   "source": [
    "Spectrum bifurcation renormalization group (SBRG) is an approximate method to diagonalize strongly disordered Pauli Hamiltonians by Clifford unitary transformations [arXiv:1508.03635](https://arxiv.org/abs/1508.03635). Given a many-body Hamiltonian $H=\\sum_a h_a P_a$ as a linear combination of a (sparse) set of Pauli operators $P_a$, the SBRG finds an approximate Clifford circuit $U$, such that the transformed Hamiltonian takes the diagonal form"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5e87769d-1e24-41f3-b531-a7437d0de6d0",
   "metadata": {},
   "source": [
    "$$H_\\text{eff}=U^\\dagger HU\\simeq \\sum_{i}\\epsilon_{i} Z_i+\\sum_{ij}\\epsilon_{ij}Z_i Z_j+\\sum_{ijk}\\epsilon_{ijk}Z_i Z_j Z_k+\\cdots$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "39bcded0-06e9-4d88-8e75-f3a99b2aa9a6",
   "metadata": {},
   "source": [
    "This is a strong disorder renormalization group method, which relies on the strong randomness of the Hamiltonian coefficients $h_a$, such that different terms in the Hamiltonian have seperating energy scales. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4f0b68dc-f0d0-4c13-8da9-240690949706",
   "metadata": {},
   "source": [
    "For example, consider a two-qubit test Hamiltonian as follows. Here $Z_0Z_1$ is a leading term and $X_0$, $X_1$ are treated as perturbations. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "91d789bd-4532-4a52-9595-527363a49eef",
   "metadata": {},
   "outputs": [],
   "source": [
    "H = - paulialg.pauli('ZZ') + 0.2 * paulialg.pauli('XI') + 0.1 * paulialg.pauli('IX')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0d137f7f-a49c-40b5-9c5f-5b5bcb4b13b1",
   "metadata": {},
   "source": [
    "Use SBRG to find the effective Hamiltonian and the Clifford unitary circuit that approximately diagonalizes the Hamiltonian."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "8ea6aef1-e20f-4907-a981-2e33b9aaaf47",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.02 IZ +1.02 ZI"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Heff, U = circuit.SBRG(H)\n",
    "Heff"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "35a9e5d2-f6bc-44b7-94a2-45742ba16bf5",
   "metadata": {},
   "source": [
    "One can check the accracy of $H_\\text{eff}$ by comparing its eigenvalues with those of $H$ by exact diagonalization."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "609b9cda-5014-4198-bfd0-d2ec91e5d781",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-1.04403065, -1.00498756,  1.00498756,  1.04403065])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "numpy.linalg.eigvalsh(H.to_qutip())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "6496db87-7764-45f4-a10c-420f0b62c06c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([-1.045, -1.005,  1.005,  1.045])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "numpy.linalg.eigvalsh(Heff.to_qutip())"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a8a1612a-cf21-4e4b-868a-f35db21ce6f7",
   "metadata": {},
   "source": [
    "One can see that the eigenvalues matches to the 2nd order of perturbations as expected."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1e976179-7e76-455b-b049-dd44640aae54",
   "metadata": {},
   "source": [
    "Consider a larger example of diagonalizing random transverse field Ising model Hamiltonian by the SBRG algorithm."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "id": "7b53ec52-f7c5-415c-9b6d-7181c45bb281",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.18 IIIIIIIIIIIIIIIX -1.06 IIIIIIIIIIIIIIZZ +1.62 IIIIIIIIIIIIIIXI +0.02 IIIIIIIIIIIIIZZI +0.43 IIIIIIIIIIIIIXII -1.70 IIIIIIIIIIIIZZII +0.56 IIIIIIIIIIIIXIII +0.50 IIIIIIIIIIIZZIII +0.21 IIIIIIIIIIIXIIII +0.92 IIIIIIIIIIZZIIII -2.35 IIIIIIIIIIXIIIII +2.32 IIIIIIIIIZZIIIII -0.11 IIIIIIIIIXIIIIII -0.87 IIIIIIIIZZIIIIII +0.50 IIIIIIIIXIIIIIII +1.42 IIIIIIIZZIIIIIII +0.39 IIIIIIIXIIIIIIII +1.03 IIIIIIZZIIIIIIII -0.11 IIIIIIXIIIIIIIII -0.45 IIIIIZZIIIIIIIII +1.99 IIIIIXIIIIIIIIII -0.59 IIIIZZIIIIIIIIII +1.21 IIIIXIIIIIIIIIII +0.81 IIIZZIIIIIIIIIII +0.44 IIIXIIIIIIIIIIII -0.39 IIZZIIIIIIIIIIII -0.95 IIXIIIIIIIIIIIII -0.85 IZZIIIIIIIIIIIII -0.88 IXIIIIIIIIIIIIII -1.80 ZIIIIIIIIIIIIIIZ +0.49 ZZIIIIIIIIIIIIII +1.15 XIIIIIIIIIIIIIII"
      ]
     },
     "execution_count": 45,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "N = 16\n",
    "H = paulialg.pauli_zero(N)\n",
    "for i in range(N):\n",
    "    H += numpy.random.randn() * paulialg.pauli({i:'Z',(i+1)%N:'Z'}, N)\n",
    "    H += numpy.random.randn() * paulialg.pauli({i:'X'}, N)\n",
    "H"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e18c83b9-d458-4a45-9a48-b1ce77f22da4",
   "metadata": {},
   "source": [
    "Many-body localization effective Hamiltonian."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "id": "52a4f8eb-ec0a-4a12-8fbf-30423387f9f5",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.87 IIIIIIIZIZIZIIII +1.56 IIIIIZIIIIIIIIII +1.04 IIIIIZIZIIIIIIII +1.97 IIIIZIIIIIIIIIII -1.49 IIIIZIZIIIIIIIII -1.85 IIIZIIIIIIIIIIII +0.49 IIIZIIIIIIIIZZII -0.52 IIIZIIIIIIIZZIII -2.18 IIZIIIIIIIIIIIII +1.41 IIZIIIIIZIIIIIII +2.13 IZIIIIIIIIIIIIII +1.08 IZIIIIIZIIZIIIII -3.67 ZIIIIIIIIIIIIIII +0.94 ZIIIIIIIIZIIIIII -0.00 ZIIZIZZZIZZZIZZZ -0.12 ZIZIIIIIZIIIIIZI -0.00 ZIZZIIZZIIZIIZZI"
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Heff, U = circuit.SBRG(H)\n",
    "Heff"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "79dabbd4-3467-45b4-b7b3-fc5eef5317bf",
   "metadata": {},
   "source": [
    "A representative approximate eigenstate (in the original basis)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "id": "ead3b4d7-5f9c-4de3-ab1e-161b669685e2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "StabilizerState(\n",
       "   +IIIIIIIIIIXIIIII\n",
       "   +IIIIIXIIIIIIIIII\n",
       "   +ZIIIIIIIIIIIIIIZ\n",
       "   +IIIIIIIIIIIIZZII\n",
       "   +IIIIIIIIIIIIIIXI\n",
       "   +IIIIIIIZZIIIIIII\n",
       "   -IIIIXIIIIIIIIIXI\n",
       "   +IIIIIIZIZIIIIIII\n",
       "   -ZIXIIIIIIIIIIIIZ\n",
       "   -IIIIIIIIIZIZIIII\n",
       "   -IXIIIXZIZIIIIIII\n",
       "   -IIIIIIZIIIIZIIII\n",
       "   +IIIIIIZIIIIIIZII\n",
       "   +IIIXIIZIIIIIZIII\n",
       "   -YIXIIIIIIIXIIIIY\n",
       "   -YXXXXXYYYYXXXXXY)"
      ]
     },
     "execution_count": 47,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rho = U.backward(stabilizer.zero_state(N))\n",
    "rho"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3eb311da-091e-497c-80f8-94d35a77eb0d",
   "metadata": {},
   "source": [
    "Half-system-size entanglement entropy (in unit of bit) of the state."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "id": "dde62d81-14cd-46c2-90ff-200a28fdea6e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2"
      ]
     },
     "execution_count": 48,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rho.entropy(numpy.arange(N//2))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d61e1c11-6bf6-414e-aca7-51638f37ed78",
   "metadata": {},
   "source": [
    "Perform the calculation of entanglement entropy over different sizes of entanglement regions for eigenstates of different random Hamiltonians."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "id": "0399b595-d5bf-4e11-b434-5680a4776d07",
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 128\n",
    "Ls = numpy.arange(N//2)\n",
    "S_lst =[]\n",
    "samples = 20\n",
    "for _ in range(samples):\n",
    "    H = paulialg.pauli_zero(N)\n",
    "    for i in range(N):\n",
    "        H += numpy.random.randn() * paulialg.pauli({i:'Z',(i+1)%N:'Z'}, N)\n",
    "        H += numpy.random.randn() * paulialg.pauli({i:'X'}, N)\n",
    "    Heff, U = circuit.SBRG(H)\n",
    "    rho = U.backward(stabilizer.zero_state(N))\n",
    "    S_lst.append(numpy.array([[rho.entropy((i + numpy.arange(L)) % N) for i in range(N)] for L in Ls]).mean(1))\n",
    "Ss = numpy.stack(S_lst).mean(0)\n",
    "Serr = numpy.stack(S_lst).std(0)/numpy.sqrt(samples)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "48dd0970-92b6-46d9-8209-16b3f23a6217",
   "metadata": {},
   "source": [
    "The entanglement entropy should scale logarithmically for random transverse field Ising model at the strong-disorder critical point."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "id": "36484b04-cfae-49bb-ab5d-9320276c9a36",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjoAAAGwCAYAAACgi8/jAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAA9hAAAPYQGoP6dpAABAD0lEQVR4nO3deXiV9Z3//9fZT3aWQCAQVhFR3EhcAJm2WqnUajtdZKadYkftDFOtRdr+Wurvqq3tNI611loFl2o7/MZWplPt2O/Q0cy3igs6soSKggUBCUtiSMienP3+/XHnnOwQ4Jz7PsvzcV3nOif3fZ9z3rlLzev6rA7DMAwBAABkIafdBQAAAKQKQQcAAGQtgg4AAMhaBB0AAJC1CDoAACBrEXQAAEDWIugAAICs5ba7ALvFYjEdPXpURUVFcjgcdpcDAABGwTAMdXR0qLy8XE7nyO02OR90jh49qoqKCrvLAAAAp+HQoUOaOnXqiOdzPugUFRVJMm9UcXGxzdUAAIDRaG9vV0VFReLv+EhyPujEu6uKi4sJOgAAZJiTDTthMDIAAMhaBB0AAJC10irovPzyy7ruuutUXl4uh8Oh3//+9yd9z6ZNm1RZWSm/369Zs2bpkUceSX2hAAAgI6RV0Onq6tKFF16ohx56aFTXHzhwQB//+Me1ZMkS1dbW6jvf+Y5uv/12/e53v0txpQAAIBOk1WDkZcuWadmyZaO+/pFHHtG0adP0wAMPSJLmzZunrVu36r777tNnPvOZFFUJAAAyRVq16Jyq119/XUuXLh1w7GMf+5i2bt2qcDg87HuCwaDa29sHPAAAQHbK6KDT0NCgsrKyAcfKysoUiUTU1NQ07Huqq6tVUlKSeLBYIAAA2Sujg440dP68YRjDHo9bs2aN2traEo9Dhw6lvEYAAGCPtBqjc6omTZqkhoaGAccaGxvldrs1fvz4Yd/j8/nk8/msKA8AANgso1t0Fi5cqJqamgHHXnjhBVVVVcnj8dhUFQAASBdpFXQ6Ozu1Y8cO7dixQ5I5fXzHjh2qq6uTZHY7rVixInH9ypUrdfDgQa1evVq7d+/Wk08+qSeeeELf+MY37CgfAACkmbTqutq6das+8pGPJH5evXq1JOnGG2/Ur371K9XX1ydCjyTNnDlTGzdu1B133KGHH35Y5eXlevDBB5laDgAAJEkOIz56N0e1t7erpKREbW1tbOoJAECGGO3f77Rq0QEAIOtFw1Is0vuI2V1N6nkLJJd9cYOgAwBAf5GQFO6RIgEpHOh7DvcMPD74dbhHigT7XsffGw1J0UhfuMk1n6iWplxk29cTdAAAmSkSHBhEIoGBAWPA6+GCSUCK9D73f4+RA60sOYSgAwBIHcMYpmVkhIAx4PwwoSXeWhK/jkCCUSDoAACkWHSY7pd4IAmeIIyM0GKSuD4kKafnvMBmBB0AyHThgBTqlIL9Hx19r/ufC3f3BpdB3Tq5OHYEOYGgAwDpJhyQelql7happ/fR3Sr1HDefA61SoKMvwETD9tYLpDGCDgBYIRLqDSwtvSHm+AhhpsXs8gGQFAQdADhd8fAyILD0ex0PLj0tUqjb5mKB3ETQAYD+hoSX1kEtMf0CTajL3loBnBRBB0D2M2JmV1FX80laXloJL0CWIegAyGxGzAwqXcekzqbe52NSV1Pfz93HzenTAHIOQQdA+jIMs6VlQHDpfXT2Bpru40yNBjAigg4Ae/W0SR0NUnt9v0eD1NlodjURYgCcAYIOgNSKRaWOD/qFmd7njgbzNWNiAKQQQQfAmQt29bXGdPQLM/GWGfYkAmATgg6Ak4tFzXEx/buW2hv6Qk2ww+4KAWBYBB0AfYKdUushqfVwv+fDZrhhrAyADETQAXJRJCi1HJSaD0jH3zefW+rMGU4AkEUIOkC262mTju2VmvZKTfvMYNNez7gZADmBoANkk0C7GWrij6b3zMHAAJCjCDpApgp0mK00/YMNoQYABiDoAJmi4wPpcK109M/SB++a69AAAE6IoAOkq2CHdOTPZrg5UmuOqwEAnBKCDpAuoiGp/h3pyA7z0fQeA4YB4AwRdAC7GIY5C+pIb4tNwy5z2jcAIGkIOoCVAu3Soa1S3Rbp8HbzZwBAyhB0gFQyDHNmVN0W83FsL91RAGAhgg6QbMFO6dA2s+Xm0DZWGwYAGxF0gGRo2me22BzaYk79ptUGANICQQc4HaEuc9p33Raz1aa72e6KAADDIOgAo9V8wGyxqdsifbBbikXtrghABjIMQzLMZ8OwuxoLBANyRiNyuOyJHAQdYCShbnPad91Wc7xNV5PdFQFIMsMwZMQMxaIx89H/dbT3dcx8bSReDz034vtjMRkx8zvMYGOGnJxS81nN/s4vNebSq235eoIO0N/x93uDzRZzXZtYxO6KgJyUCB+Dg8TgMNHvtdEbLAYEkcGfMejzjFiupY7cQ9BBbjMMqeEdad/L0sE3pM5jdlcEpC3DMBIhwYgNbLUYGi5GCiPDt3oMbiWxpdXDIblcTjlcTjmdTjldDjldTvPh7Pd6yM8OOZ297+t9PeC402E+HIOfJYfDYcMvarFl35djWpVtX0/QQW5q2ie995IZcNjxG1nEMAxFIzHFIlFFIzFF489h83lgK8fAoGGcJKjY1fqRCA39AsiQMDLk9cCg4hgmmAz+XIczB0KHHbx+yabxORJBB7mkvUF670Uz4LTU2V0NMCIjZigSjioSjigaMp8joah5LBRRJBxVNBJVrF+AiQeaWNSapQ0SwcE5fKuHY7gAcqJwMlwg6dcSApwugg6ymxGTDr4pvfN/zC0Xcm4UIOwWi8Z6Q0ok8RwdFFri56K9x6KRMw8rDodDLo9TLrdLLrf57HS75OrfyjEgnPRv5ThRC0rvz4QPZAiCDrJTT6v07vPSro10TSEpDMNQLBIbGFDirSvxIBM/1i/InEl3j8vjkjv+8Lrl8sZ/dvcLMf2DTF+wcbqcSfztgcxF0EF2adgtvfMH6cCrUjRsdzVIU4Zh9LagxMNJ/9aVqKKDgky8Jea0GwQdktvjltvbF1rcHpdcvc/m8YHnXR4XrSZAEhB0kPnCAWnvn6Rd/yU177e7GtjEiBkKhyIKByOKBMMKB83X5rGwIr0/R8JmF9HpcroccnncfYGlfyvLoCATf3a6nYQWwCYEHWSu1sPm2Js9/2NuyYCsYxiGouFoIrCYYaVfiAn2hZjIaYQXl9uZaD1JtKoMaHHpF1p6W1zoEgIyC0EHmSUWld5/w+yeOvpnu6vBGTAMQ+FARKFASKFAWOFAWKFAWKEe89lsmYmYK8mOlkPyeN3y+Dzy+Nxy+9zyxB9ej9w+94BAw3RiIPsRdJAZAu3S7j+a3VMs6pf2DMNQJBgxg0vvI9wvxJhBJjzqMS8uj6tfYIkHGE/vMU/iHONaAAxG0EF6a9onvf2cufZNNGR3NVDfQN5wvxATCoQV7gn1CzWjbIlxSF6fR548j7z+fo88T1+rjNdNdxGA00bQQfoxYtKB16Sdz0kNb9tdTc6JdykFe0IKdYcU7DEf4X6tMaOdMu3xueXN88jr98rjd8vr9/YFGb8ZZGiBAZBKBB2klwOvSVvWs3JxihmGoVBPWIGuoIJdQQW6Qwp2BXvDTXhUrTFubzzEmKGlf2uM12+2yDAGBoDdCDpID0f+LL35S6nxL3ZXklVi0ZgCXUEFOoPq6Qwo0Bk0w0136MStMg7Jl+eVN88rX743EV7MIOOVx0d3EoDMQNCBvY7tkd78lXS41u5KMpoRMxTo6gszPR0B9XQGFewOjjjg1+FwyFfglb/AJ1+BT/78eKgxu5dojQGQDQg6sEfrYTPgHHjN7koyimEYCnaH1NMRUKDTDDM9HQEFu0Ijdje53E75i/zKK/TLX+CTv9Anf4FP3jwP42MAZD2CDqzV0yptfUp694/mmjg4oWgkpu62bnW2dqurxXweaVVfp8upvEJfItTEXzPgF0AuI+jAGpGg9NYz0o7fSuEeu6tJS/EBwl2t8WDTpe6OwJCuJ4fTYQaZIp/8vYEmr8gvj58WGgAYjKCD1DJi0p7/a86k6mqyu5q0EovG1N3eYwablm51tXYrHIwMuc7jc6twbIEKxuarcEy+8or9cjoZCAwAo0HQQeocfUva/CgbbfaKRqLqaulWR0uXOo93qautZ+jMJ4eUX5ynwjH5iWDjzfPaUzAAZAGCDpIv1C298YS5ZcNo1/jPQpFwVJ29oabjeJe623uG3A6316WCMfkqHGO22BSU5DFtGwCSiKCD5KrbIr3y85zcjyoSiqjjeF+w6ekIDLnGm+dR4bgCFY0tUOG4AvnyvYyrAYAUIuggOQIdZjfV3v9rdyWWCQfCiW6ojuNdCnQGh1zjy/eawWZcb7ChGwoALEXQwZnb/6r06lqpp8XuSlIq1BNSR2+o6TzepWD30E1G/YW+RKgpHFsgr99jQ6UAgLi0Czpr167Vj3/8Y9XX1+u8887TAw88oCVLlox4/VNPPaV7771Xe/fuVUlJia655hrdd999Gj9+vIVV56hoSHr5QXNWVZaJL8zXebwr0WoT6gkPuS6vyD8g2Hh8afd/KQDIaWn1X+UNGzZo1apVWrt2rRYvXqxHH31Uy5Yt065duzRt2rQh17/66qtasWKFfvrTn+q6667TkSNHtHLlSt1yyy169tlnbfgNckigXXr+B1mzu7hhGAp0BgcEmyFTvXtnRPUPNm6Py56CAQCj4jBGs02xRS677DItWLBA69atSxybN2+ePvWpT6m6unrI9ffdd5/WrVunffv2JY79/Oc/17333qtDhw6N6jvb29tVUlKitrY2FRcXn/kvkQvajkh/vMt8zmDRSFTtTZ1qa2xX27EORUIDVxx2OBwqGJOXGDxcMDZfLjfBBgBOySeqpSkXJf1jR/v3O21adEKhkLZt26Zvf/vbA44vXbpUmzdvHvY9ixYt0p133qmNGzdq2bJlamxs1H/8x3/o2muvHfF7gsGggsG+QaPt7e3J+QVyRf3b0gs/MFt0MlCwJ2QGm8YOdTR3DdgfyuF0qHBsvgrHmoOHC8bkM9UbADJc2gSdpqYmRaNRlZWVDTheVlamhoaGYd+zaNEiPfXUU1q+fLkCgYAikYiuv/56/fznPx/xe6qrq/X9738/qbXnjL0vSpt+KkWHjlVJZ8HukI7Xt6qlvm3IlG9fvlclE4s0ZmKxCsbms+IwAGSZtAk6cYPXFDEMY8R1Rnbt2qXbb79d3/3ud/Wxj31M9fX1+uY3v6mVK1fqiSeeGPY9a9as0erVqxM/t7e3q6KiInm/QLba9mtp6/9ndxWjFglF1NLQpuajrepq6R5wrnBsvkomFmvMxCL5CnysYwMAWSxtgk5paalcLteQ1pvGxsYhrTxx1dXVWrx4sb75zW9Kki644AIVFBRoyZIl+uEPf6jJkycPeY/P55PP50v+L5CtDMOcWfXuf9tdyUnFojG1Nrbr+NFWtR/rHNAtVTS+QOPKx2jMxGK5vWnzzx4AkGJp8198r9eryspK1dTU6K//+q8Tx2tqavTJT35y2Pd0d3fL7R74K7hc5mDRNBpjndk2P5LWIccwDHU0d+n40Va1NLQpFo0lzuUV+TV+yhiNnTyG9WwA9HE4JZdXcnslp0dye3qfvX3PLq/k8vQ+vIOehzvW/9yga5xuyekyv3fAs0vKhRZlb4GtX582QUeSVq9erS9+8YuqqqrSwoUL9dhjj6murk4rV66UZHY7HTlyROvXr5ckXXfddfryl7+sdevWJbquVq1apUsvvVTl5eV2/irZ4c1/ld5+zu4qhjAMQz3tATUfbVVLfeuAaeDePI/GlY/RuMljlFfkt7FKAMNyugcFgUGBwzVcyBjumuGCR7/j7sHH+13rZPZkLkmroLN8+XI1Nzfr7rvvVn19vebPn6+NGzdq+vTpkqT6+nrV1dUlrv/Sl76kjo4OPfTQQ/r617+uMWPG6Morr9S//Mu/2PUrZI8dv5Vqn7a7igFCPWE1H23R8SOtCnT1zZxzeVwaO6lE48vHqGBsPmNugJE4nIPCg2dQC0b/oOAdPigMCSMnek//0NH7s4MB/7BWWq2jYwfW0RnGrv+SXnnI7iok9Y67+aBdzUda1N7UmTjucDo0ZmKxxpWPUfGEQmZLITs5nJI332z69xYMfO3xS26/5MkzX3vyen+Ov44fH3QOyBIZt44O0sTeP5n7VtnIMAx1tfao+UiLWupbFY30jbspHFtgjruZVCIXqxIjnQ0JKQXDB5bB1/j6h5k8u38LIOMRdNDn/dell+6XjNjJr00BI2boeH2rGvYfG7ATuDfPo/FTxmp8+Rj5CpgxBys4JG/eCAGlcGggGXyNr5CQAqQJgg5Mh2ul/6mWYtGTX5tksWhMTYdb9MGBY4mNMx1Oh8ZOKlHp1LEqHFfAuBucOqdb8hdJ/hLJVyT5i81nX+HAwBIPJgPCSh5jSYAsQdCB1FJnbutg8YrHkXBUx+qa1fh+U2KfKbfXpYkzSjVh2ng2zEQfl9cMKvGHL/66aOCxvH7nvPl2Vw0gDRB0cl2wS3r+bincY9lXhoNhfXCgScfqjifWvfHmeVQ2c4JKp45lf6lc4HCaYSRvjJQ3Vsof2+9177O/pC/EMIgWwGki6OQyIyb96V8s24U8FAjrg/3HdOzQcRkxc7Kfv9CnSbMmaNzkMXI46Z7KbA4zlOSPNYNK3pi+1/HwEg80/hLWMgFgCYJOLntzvVS3JeVfE+oJqWH/MTUdakmsWJ1fkqfJsyeqZGIR42/SnbdAKhg/MKgkwkw8xIyT8ggvANIPQSdX7XtF2rEhpV8R7A6pYV+jmo+0JgJO4dh8TT5roorGFxJw0oE3XyqcKBWVmc8F46WC0n7PpXQbAchoBJ1c1HzAnEaeIpFwVPXvfaDGg81S73KUReMLNHm2GXBgIZdHKpokFU/ufUzqCzZFE81ZSACQxQg6uSbQYQ4+jgSS/tGGYajpcIuO7mlIzKIqLi3U5LMmqnCsvZu6ZTVvYb8gM1kq6fe6YDzTpAHkNIJOLolFzbVyOhqS/tGdLV2q23VUPe1mgPIX+FRx7mQVl9JikBwOszVm/Cxp/EypZEpvmCk3p1gDAIZF0Mkl//ukdKQ2qR8Z6gnr8F/q1VLfJklyuZ2aPKdME6eNZxbV6XL7pHEzekNN72PcDNaFAYDTQNDJFUffkt56Jqkf2XT4uA7tOqpY1ByIU1oxVuVzJsnj45/VqBWU9rXSjJ/d21pTTncTACQJf5FyQTQkvfxg0j4uFo2pbtdRNR9ukSQVjM3XtHnlyi9hb58TKiqTJsyRSudIE86SSs8y150BAKQMQScXbN+QtEUBg11B7autU0+HORan/OwyTZo1ganigxVONENN/EGoAQBbEHSyXUudtOPfk/JRrR+06/23DikaicntdWnmhdNUXMp0cRVO6GulmXC2GWrySuyuCgAggk52Mwzp5Z9LsciZfUzM0JG9H+iD/cckSQVj8jXr4mny+j3JqDKzFJQObKWZMMdcKRgAkJYIOtns3eelhrfP6CPCwYj276hT5/EuSdLEGeM1de7k3JhRlT9uUPfTHHPrAwBAxiDoZKvuFul/nzijjwgFwtrz5gEFu4Jyupyafv4UjZs8Jjn1pRt/iTRxbl8rzYSzpYJxdlcFADhDBJ1s9fpjUrDztN8e6gmZIac7JI/fozmXzFBeYRbteeQvlsovkCZfYD6Pm253RQCAFCDoZKND26T3Xjrttwe7Q9rz5n6FesLy5nl09qWz5Mv3Jq8+O3jypCkXS1MuNIPN2OkSM8UAIOsRdLJNOCC98tBpvz3YFdSeNw8oFAjLl+/V2ZfOlDcvQ0NOUZk0/TJp2qVmuHHl4OBpAMhxBJ1ss+2p097LKtAV1J7/3a9wMCJfgc8MOZk0s8rhlMrOkaZdZgYcuqMAIOcRdLJJ2xFp57On9daezoD2vHlAkWBE/kIz5Hh8GRByvPnS1MrelptLWJQPADAAQSebbN9g7lB+ino6Atrz5n5FQlHlFfk155KZ6b1fVfGkvlabyedLrjSuFQBgK/5CZIuOD6T3/nTKbwv1hLR3ywFFQlHlF5shx+1Nw38W42ZIMxZJs64wN74EAGAU0vAvGk5L7b+fcmtOJBzV3q3vK9zbXTXn0llye1wpKvA0lM6RZi2WZi6Wxky1uxoAQAYi6GSDriZpT80pvSUWi2l/7UEFOoPy+NyaUzUjPUJOwXjpnGXS3I+as6YAADgDBJ1ssOM/pGh41JcbhqGDbx9RR3OXnC6nzqqaYfMUcoc09SLp3Gul6ZdLzjQIXACArEDQyXQ9rdK7/31Kb6l/r1HHj7RKDmnWxdOUX5yXktJOylckzb3aDDgl5fbUAADIagSdTPfn30mR4Kgvbzp8XPXvNUqSpp83RSUTilJV2chKz5LmXy/N/pDkztDFCAEAGYGgk8kCHdKu/xr15e1NHTr49hFJ0qTZE1RaYeGmlU63Oah4/vXSpHOt+14AQE4j6GSync9K4Z5RXdrd3qN92+skQxpXPkblcywa6Js3Vpq3zOyeYjdwAIDFCDqZKtQlvf2HUV0aDUe1b/tBxaIxFY4r0PT5U+RI9YaWRZOkBX8jzbmSPaYAALYh6GSqt5+TQp2juvTwX+oTO5HPXjBdTpczdXUVTpAu/htp7lJWLAYA2I6/RJkoHJB2/n5Ul7Yd61DToRZJ0ozzp6ZurZz8cdLFy6V510guBhgDANIDQScT7fo/UqD9pJdFwlEd3HlYkjRx+ngVjS9Mfi3+EuniG8wxOG5f8j8fAIAzQNDJNJGQ9NYzo7r08O56hYMR+fK9mjJ3UvJrmf1X0hW3SX4bpqgDADAKBJ1Mc+BVqbvlpJe1Nrar+Uhvl9UFU5M7LsdfLF1xqxl0AABIYwSdTPPeppNeEglFEuvllM0sVeHYguR9//TLpL/6mpQ/NnmfCQBAihB0MkmgQzq8/aSXHdp1VJFgRP4CX/LWy/EWSIv+0dyyAQCADEHQyST7X5FikRNe0tLQpuP1bZKS2GVVfqH0ka+bU8cBAMggBJ1Msu/E3VbhYER17/Ru8TBrggrG5J/5d8692uyqYkdxAEAGIuhkiq5mqf7tE15yaNdRRUJR+Qt9mnzWxDP/zgV/K12y4sw/BwAAmxB0MsW+lyUjNuLpjuZOtTS0SQ5p5gUVZ9Zl5XBKV3zFXBsHAIAMRtDJFO+9NOIpwzB0ZE+DJGlCxTjll+Sd/ve4vNJV35JmLjr9zwAAIE0QdDJB21Hp2J6RTzd2qKu1Rw6nQ5Nnn0GXla9QuuZ70qTzTv8zAABIIwSdTHCCQcj9W3PKZpTK4z/NncILJ0of/4E0dtrpvR8AgDRE0MkEJ+i2On60VYHOoFxup8pmneb076Iy6fofM30cAJB1krgvAFKi+YDUUjfsqVgspqN7GyVJZbMmnN7O5PnjpU9UE3IAAFmJoJPuTtCa03yoRaGekNxetyZOLz31z/aXSJ/4kVQ8+fTrAwAgjRF00t0I43Ni0Zjq95mtOZPPmiCX+xT/p/QWStf+M2NyAABZjaCTzhp2Sx0fDHuq8f1mhYMRefM8Kq0Yd2qf68mTPn63VDo7CUUCAJC+CDrpbN9Lwx6OhKNqOHBMklQ+p0xO5yn8z+j2mVPIy+adeX0AAKQ5gk66ikWlfa8Me+qDA8cUDZtbPYwrHzP6z3S6pavvlMovSE6NAACkOYJOujr6ltTTMuRwOBhW4/vNkqQpZ5fJ4XCM8gMd0lX/jzTtkiQWCQBAeiPopKsRZlvV7zumWDSmgpI8lUwsHv3nXfBpadaS5NQGAECGIOiko2hYOrB5yOFQT1hNdcclSeVzJ42+NWfCHOnSLyWxQAAAMkPaBZ21a9dq5syZ8vv9qqys1CuvDD9OJS4YDOrOO+/U9OnT5fP5NHv2bD355JMWVZsih2ulUOeQw81HjsswDBWOzVfx+MLRfZYnT7rq25KLRbABALknrf76bdiwQatWrdLatWu1ePFiPfroo1q2bJl27dqladOGX+/lhhtu0AcffKAnnnhCZ511lhobGxWJRCyuPMk+2DXkkGEYajrcKkmnNp38iq9IJeVJKgwAgMySVkHn/vvv180336xbbrlFkvTAAw/o+eef17p161RdXT3k+v/+7//Wpk2btH//fo0bZ/7xnzFjxgm/IxgMKhgMJn5ub29P3i+QLMfeG3Ko83iXQj0hOV1OjZ1UMrrPOesj0tkfTXJxAABkjrTpugqFQtq2bZuWLl064PjSpUu1efPQ8SqS9Nxzz6mqqkr33nuvpkyZorPPPlvf+MY31NPTM+L3VFdXq6SkJPGoqKhI6u+RFE17hx46Ys7AGldeIqdrFP+zFU+SltyW7MoAAMgoadOi09TUpGg0qrKysgHHy8rK1NDQMOx79u/fr1dffVV+v1/PPvusmpqa9JWvfEXHjx8fcZzOmjVrtHr16sTP7e3t6RV2OhqlwMBWpmgkqtaGNknS+CljT/4ZTrc5Lsebn4oKAQDIGGkTdOIGzyQyDGPE2UWxWEwOh0NPPfWUSkrM7pz7779fn/3sZ/Xwww8rLy9vyHt8Pp98Pl/yC0+WY0Nbc1rq2xSLGvIV+FQwZhThpeqL0sS5KSgOAIDMkjZdV6WlpXK5XENabxobG4e08sRNnjxZU6ZMSYQcSZo3b54Mw9Dhw4dTWm/KDNdtddjstiqdOvbkU8qnXCRd9LkUFAYAQOZJm6Dj9XpVWVmpmpqaAcdramq0aNGiYd+zePFiHT16VJ2dfVOx9+zZI6fTqalTp6a03pQZNBA50BlUV2u35JDGn2y7B2+h9JFvSqNeLRkAgOyWNkFHklavXq1f/OIXevLJJ7V7927dcccdqqur08qVKyWZ42tWrFiRuP7zn/+8xo8fr7//+7/Xrl279PLLL+ub3/ymbrrppmG7rTLCoBad5t5ByCWlRfL4PSd+b9UXpIJT3MkcAIAsllZjdJYvX67m5mbdfffdqq+v1/z587Vx40ZNnz5dklRfX6+6urrE9YWFhaqpqdFXv/pVVVVVafz48brhhhv0wx/+0K5f4cwMGohsxIxE0Bk/9SSDkMdUSOddl8rqAADIOA7DMAy7i7BTe3u7SkpK1NbWpuLiU9g7KhUOvCa90BfS2hrb9d62g3J7XDr/ynPkdJ6gAW7ZD6RpVRYUCQCA/Ub79zutuq5y3qDxOfFByOPKx5w45Ey7hJADAMAwCDrppKkv6ERCEbU1dkg6SbeV0y0t/IdUVwYAQEYi6KSTfmvoHD/aKsMwlF/sV37xCQZWn3edNCZDZ5gBAJBiBJ100XlMCrQlfox3W52wNcdfLFV+PtWVAQCQsQg66eLYnsTL7rYe9XQE5HA4NG7ymJHfc8kKyVeY+toAAMhQBJ100W8gcnwDzzFlxXJ7R1gBYNxM6ZxrrKgMAICMRdBJF70DkWPRmI4fbZV0km6rRf8oOV0WFAYAQOYi6KSL3oHI7U0dioaj8vjcKi4doVtqxiJpyoUWFgcAQGYi6KSDfgORO453SzK7rYbdwNPlkRZ+2crqAADIWASddNBvWnlXqxl0CsbkD3/tnKuk4klWVAUAQMYj6KSD+PicWEzd7T2SThB0zvuEVVUBAJDxCDrpoLdFp6cjICNmyOVxyZfvHXpd2blS6WyLiwMAIHMRdNJBb4tOX7dV3vDjc+azOzkAAKeCoGO3zmNST6skqau1t9uqZJhuq7yx0szFFhYGAEDmI+jYbbQDkeddY864AgAAo0bQsVtvt1UkFFGwOyRpmKDjdEnnftzqygAAyHgj7C9wcuFwWA0NDeru7taECRM0bty4ZNaVO3pbdOKtOb4Cn9yeQSseT79cKii1ujIAADLeKbXodHZ26tFHH9WHP/xhlZSUaMaMGTr33HM1YcIETZ8+XV/+8pe1ZcuWVNWanRIDkc3xOYVj8oZecx6DkAEAOB2jDjo//elPNWPGDD3++OO68sor9cwzz2jHjh36y1/+otdff1133XWXIpGIrr76al1zzTXau3fvyT801/UfiNzWOz5n8EDksdPY7gEAgNM06q6rzZs368UXX9T5558/7PlLL71UN910kx555BE98cQT2rRpk+bMmZO0QrNSb2uOYRgjD0SmNQcAgNM26qDz29/+dlTX+Xw+feUrXzntgnJK7/icQFdQ0UhMDqdDeUX+vvOePGnOlTYVBwBA5mPWlZ0Gjc8pKMmTw9lvocCzPyp5R9gKAgAAnFRSg86ePXu0eDGL2o3ascErIg/utmJfKwAAzkRSg044HNYbb7yRzI/MXl1NUk+L+XK4oFN+oTkQGQAAnDa6ruzSOz4nGomppyMgaVDQYRAyAABn7JQWDFy5cqUqKyt18cUX64ILLpDXO8wO2xid4+9Lkrp7p5V7fG55/b1bPBROkGZcblNhAABkj1MKOm+99ZaeeuopdXV1yePx6Nxzz9WCBQtUWVmpBQsWyOmkgWjUQmbA6WrrHYjcvzVn5mJz2wcAAHBGTinobN68WYZh6N1339X27dsTj2eeeUZtbW2SJIfDcZJPgSQpbHZXDTs+Z9J8OyoCACDrnPJeVw6HQ/PmzdO8efP0hS98IXF837592rZtm3bs2JHM+rJXJDjyQoGTz7OpKAAAssuog05dXZ2mTRt5FtDs2bM1e/Zs3XDDDZKkI0eOaMqUKWdeYbaKBBQOhBUORiSHuYaOJGnMVClvjK2lAQCQLUY9qOaSSy7Rl7/8Zb355psjXtPW1qbHH39c8+fP1zPPPJOUArNWJJBYKDCvyC+nq/d/CrqtAABImlG36OzevVs/+tGPdM0118jj8aiqqkrl5eXy+/1qaWnRrl279M4776iqqko//vGPtWzZslTWnfkiwcRGnoV0WwEAkBKjbtEZN26c7rvvPh09elTr1q3T2WefraampsQu5V/4whe0bds2vfbaa4Sc0QgH1MlAZAAAUuqUByP7/X59+tOf1qc//elU1JMzjFCPuuNTy0t6g07BeKl4ko1VAQCQXc5o4ZtAIJCsOnJOd1OTjJghl8clX0Hvwou05gAAkFRnFHT6b+BZU1NzxsXkkq5jzZJ6dyyPrz00maADAEAynVHQicViidff+ta3BpxbunTpmXx01utqbpU0eHwOA5EBAEimMwo6/VdBNgxjwLljx46dyUdnva7jHZL6BR1fkTRuhn0FAQCQhc4o6Bw7dkx/+MMf9P777w/Z+oGtIEYWaWlQsDskqd9CgZPOlbhnAAAk1SnPuurvjjvu0LPPPqsf/OAH2r9/vxYtWqS5c+dq7ty5On78eLJqzDpd726VJPkKvHJ7e/8nYCAyAABJd0pB55133pHP59NZZ50lSVq9evWA8/v379fbb7+tt99+e8BAZQzUtXeHpH7TyiUWCgQAIAVOKeisXr1a5513nu6///7Esf/8z//Ur3/9a02cOFF33HGHrr/+el1//fVJLzSbRNqaJEnePI95wO2TSufYWBEAANnplMbo/PnPf9ZnPvOZxM+7d+/W5z73Ob322mt6+umnddlll+no0aNJLzLbxAKdkiRXfH+riedIrjPqRQQAAMM4paDT1tamioqKxM/r16/XrFmzdPDgQR0+fFgXXXSR7rnnnqQXmW1iPV2S1G8jT7qtAABIhVMKOlOnTlV9fX3i5//5n//RDTfcIJfLJZ/PpzVr1uiFF15IepHZJhow97hyuntvPwsFAgCQEqcUdK6++urE+JyDBw+qtrZWV199deL87NmzdejQoeRWmIVi8aDjckpOl1Q2z+aKAADITqc0MOTOO+/UxRdfrFmzZikQCKiiokJXXHFF4vwHH3ygwsLCpBeZbWJBM+i4XE5p/GzJ47e5IgAAstMpBZ0pU6Zoy5YtevDBB9Xa2qrbbrttwMKAf/rTn3T22WcnvchsEw2am6E63S66rQAASKFTnuozffp0/eQnPxn23K5du/TZz372jIvKdrFQb9BxORmIDABACiV1TvP69euT+XFZKxYyt39wEXQAAEipM9rrCqfOMAxFQ0FJknPcFCmvxOaKAADIXgQdixmRkBSLSZKcU863uRoAALIbQcdi8anlkuSquMi+QgAAyAEEHYvFg47D4ZBj6oU2VwMAQHYj6FgssSqyxyUVldlcDQAA2Y2gY7FYoHefKy+LBAIAkGoEHYvFu65c/nybKwEAIPsRdCyW6Lry5dlcCQAA2S/tgs7atWs1c+ZM+f1+VVZW6pVXXhnV+1577TW53W5ddNFFqS3wDCW6rgg6AACkXFoFnQ0bNmjVqlW68847VVtbqyVLlmjZsmWqq6s74fva2tq0YsUKXXXVVRZVevrougIAwDppFXTuv/9+3Xzzzbrllls0b948PfDAA6qoqNC6detO+L5//Md/1Oc//3ktXLjwpN8RDAbV3t4+4GGlaLxFh6ADAEDKpU3QCYVC2rZtm5YuXTrg+NKlS7V58+YR3/fLX/5S+/bt01133TWq76murlZJSUniUVFRcUZ1n6pYoEeS5MwrtPR7AQDIRWkTdJqamhSNRlVWNnBtmbKyMjU0NAz7nr179+rb3/62nnrqKbndo9ufdM2aNWpra0s8Dh06dMa1n4pYj9mC5MorsvR7AQDIRUndvTwZHA7HgJ8NwxhyTJKi0ag+//nP6/vf/77OPvvsUX++z+eTz+c74zpPV7S7QxItOgAAWCFtgk5paalcLteQ1pvGxsYhrTyS1NHRoa1bt6q2tla33XabJCkWi8kwDLndbr3wwgu68sorLan9VMR6eoNOPruWAwCQamnTdeX1elVZWamampoBx2tqarRo0aIh1xcXF2vnzp3asWNH4rFy5UrNnTtXO3bs0GWXXWZV6ack1tMpSXIVFNtcCQAA2S9tWnQkafXq1friF7+oqqoqLVy4UI899pjq6uq0cuVKSeb4miNHjmj9+vVyOp2aP3/+gPdPnDhRfr9/yPF0klhHJ5+gAwBAqqVV0Fm+fLmam5t19913q76+XvPnz9fGjRs1ffp0SVJ9ff1J19RJd9Ge3pWRmV4OAEDKOQzDMOwuwk7t7e0qKSlRW1ubiotT38qy+/YPq7vuPZ31//6rSqrSf4FDAADS0Wj/fqfNGJ1ckVhHhxYdAABSjqBjsWiQoAMAgFUIOhaLhQKSJJe/wOZKAADIfgQdi8VCQUmS00eLDgAAqUbQsVAsHJIRjUqSnHkEHQAAUo2gY6FYsDvxmhYdAABSj6BjoVjvGjoOp0tOj9fmagAAyH4EHQvFW3ScXkIOAABWIOhYKBrf/sFr3+7pAADkEoKOhWIBs0XH5cuzuRIAAHIDQcdC8aDj9BN0AACwAkHHQtF40GHGFQAAliDoWCjWO0bHlceqyAAAWIGgY6F40HGy/QMAAJYg6FgoGt+5PK/Q5koAAMgNBB0LJbqu8ottrgQAgNxA0LFQYtZVfpHNlQAAkBsIOhaK9nRKkpx5tOgAAGAFgo6FYt0dklhHBwAAqxB0LBQLmC06LmZdAQBgCYKOhaI98enlLBgIAIAVCDoWSuxeTosOAACWIOhYKNbbouOiRQcAAEsQdCwUC/YuGEjQAQDAEgQdC0UTQYeuKwAArEDQsVAsGJBE1xUAAFYh6FjEiIRlRCOS6LoCAMAqBB2LRHu3f5DougIAwCoEHYvEp5bL6ZTT47W3GAAAcgRBxyLxDT1dPr/NlQAAkDsIOhZJrIrsJegAAGAVgo5FEqsi+9jQEwAAqxB0LMKqyAAAWI+gYxFWRQYAwHoEHYtEA71jdPIKba4EAIDcQdCxSGLWFUEHAADLEHQsEou36OQX21wJAAC5g6BjkfjKyKyKDACAdQg6FunruiLoAABgFYKORRJdV7ToAABgGYKORRJdVz6mlwMAYBWCjkXiLTquPIIOAABWIehYJNbdKYkWHQAArETQsUiUMToAAFiOoGMRuq4AALAeQcciMQYjAwBgOYKORRKzrlhHBwAAyxB0LBILBSTRogMAgJUIOhYwohEZ4ZAkycVgZAAALEPQsUB8fI4kORmMDACAZQg6FogGe4OO0ymH22tvMQAA5BCCjgViPb0bevry5HA4bK4GAIDcQdCxQKy3Rcfp9dtcCQAAuYWgY4FoT3xVZMbnAABgJYKOBfpWRWbGFQAAViLoWCDRdcXUcgAALEXQsUC0J74qcpHNlQAAkFvSLuisXbtWM2fOlN/vV2VlpV555ZURr33mmWd09dVXa8KECSouLtbChQv1/PPPW1jt6CRadOi6AgDAUmkVdDZs2KBVq1bpzjvvVG1trZYsWaJly5aprq5u2OtffvllXX311dq4caO2bdumj3zkI7ruuutUW1trceUnlhijQ9cVAACWchiGYdhdRNxll12mBQsWaN26dYlj8+bN06c+9SlVV1eP6jPOO+88LV++XN/97ndHdX17e7tKSkrU1tam4uLi06r7ZI7++seq//efacLHb9S0f/jnlHwHAAC5ZLR/v9OmRScUCmnbtm1aunTpgONLly7V5s2bR/UZsVhMHR0dGjdu3IjXBINBtbe3D3ikWmLnclp0AACwVNoEnaamJkWjUZWVlQ04XlZWpoaGhlF9xk9+8hN1dXXphhtuGPGa6upqlZSUJB4VFRVnVPdoxHrougIAwA5pE3TiBm+RYBjGqLZN+M1vfqPvfe972rBhgyZOnDjidWvWrFFbW1vicejQoTOu+WT6ppezYCAAAFZy211AXGlpqVwu15DWm8bGxiGtPINt2LBBN998s37729/qox/96Amv9fl88vl8Z1zvqejruiLoAABgpbRp0fF6vaqsrFRNTc2A4zU1NVq0aNGI7/vNb36jL33pS/r1r3+ta6+9NtVlnpZYIL6pJ0EHAAArpU2LjiStXr1aX/ziF1VVVaWFCxfqscceU11dnVauXCnJ7HY6cuSI1q9fL8kMOStWrNDPfvYzXX755YnWoLy8PJWUlNj2ewwWn17OOjoAAFgrrYLO8uXL1dzcrLvvvlv19fWaP3++Nm7cqOnTp0uS6uvrB6yp8+ijjyoSiejWW2/Vrbfemjh+44036le/+pXV5Y8o2tMpiVlXAABYLa3W0bGDFevo7PzyZQodO6Jz7n1OBWcvSMl3AACQSzJuHZ1slui6okUHAABLEXQsEA32SJKcDEYGAMBSBJ0UM6JRGaGgJMnFYGQAACxF0EmxWKgn8Zp1dAAAsBZBJ8Xi2z/I4ZTDY+1ChQAA5DqCTopFEwOR80a1lQUAAEgegk6KxQJm15XLl2dzJQAA5B6CTooxtRwAAPsQdFIsGt+5nBlXAABYjqCTYvHByC5/oc2VAACQewg6KRbfudyZx9RyAACsRtBJsWg86LAqMgAAliPopFh8MLKLwcgAAFiOoJNisfhgZFZFBgDAcgSdFIv2EHQAALALQSfF4i06dF0BAGA9gk6K9S0YSIsOAABWI+ikGF1XAADYh6CTYgxGBgDAPgSdFGN6OQAA9iHopFhiZWSCDgAAliPopFhiZWR/ns2VAACQewg6Kda3qSctOgAAWI2gk2IMRgYAwD4EnRQyYjHFgj2SGKMDAIAdCDopFA85kuSiRQcAAMsRdFIoPrVcDoccXr+9xQAAkIMIOinUN7U8Xw6Hw+ZqAADIPQSdFEpMLffRbQUAgB0IOimUWBU5j4HIAADYgaCTQjFadAAAsBVBJ4WivS06Tlp0AACwBUEnheItOkwtBwDAHgSdFKLrCgAAexF0UoiuKwAA7EXQSaFE1xUtOgAA2IKgk0L9FwwEAADWI+ikUDS+czldVwAA2IKgk0Kxnt4FA+m6AgDAFgSdFIrRogMAgK0IOikU7YlPL8+zuRIAAHITQSeFEntd+WnRAQDADgSdFKLrCgAAexF0Uqiv64rByAAA2IGgk0LxFh32ugIAwB4EnRSKj9FhwUAAAOxB0EkRIxbrtzIyY3QAALADQSdFYqFA4rWLwcgAANiCoJMi8W4rORxyeP32FgMAQI4i6KRIotvKly+Hw2FzNQAA5CaCTopEGYgMAIDtCDopEgv0SGJqOQAAdiLopAhTywEAsB9BJ0X6uq6YcQUAgF0IOikSH4zMhp4AANiHoJMifYsF5tlcCQAAuYugkyKsigwAgP0IOinC9HIAAOyXdkFn7dq1mjlzpvx+vyorK/XKK6+c8PpNmzapsrJSfr9fs2bN0iOPPGJRpScWn3XF9HIAAOyTVkFnw4YNWrVqle68807V1tZqyZIlWrZsmerq6oa9/sCBA/r4xz+uJUuWqLa2Vt/5znd0++2363e/+53FlQ8VX0eHrisAAOyTVkHn/vvv180336xbbrlF8+bN0wMPPKCKigqtW7du2OsfeeQRTZs2TQ888IDmzZunW265RTfddJPuu+8+iysfiq4rAADslzZBJxQKadu2bVq6dOmA40uXLtXmzZuHfc/rr78+5PqPfexj2rp1q8Lh8LDvCQaDam9vH/BIBaaXAwBgv7QJOk1NTYpGoyorKxtwvKysTA0NDcO+p6GhYdjrI5GImpqahn1PdXW1SkpKEo+Kiork/AKDORxyeHxy+pheDgCAXdx2FzDY4J2+DcM44e7fw10/3PG4NWvWaPXq1Ymf29vbUxJ2Zn/rsQH1AAAA66VN0CktLZXL5RrSetPY2Dik1SZu0qRJw17vdrs1fvz4Yd/j8/nk8/mSU/QonCikAQCA1Eqbriuv16vKykrV1NQMOF5TU6NFixYN+56FCxcOuf6FF15QVVWVPB5PymoFAACZIW2CjiStXr1av/jFL/Tkk09q9+7duuOOO1RXV6eVK1dKMrudVqxYkbh+5cqVOnjwoFavXq3du3frySef1BNPPKFvfOMbdv0KAAAgjaRN15UkLV++XM3Nzbr77rtVX1+v+fPna+PGjZo+fbokqb6+fsCaOjNnztTGjRt1xx136OGHH1Z5ebkefPBBfeYzn7HrVwAAAGnEYeT4aNn29naVlJSora1NxcXFdpcDAABGYbR/v9Oq6woAACCZCDoAACBrEXQAAEDWIugAAICsRdABAABZi6ADAACyFkEHAABkLYIOAADIWgQdAACQtdJqCwg7xBeGbm9vt7kSAAAwWvG/2yfb4CHng05HR4ckqaKiwuZKAADAqero6FBJScmI53N+r6tYLKajR4+qqKhIDocjaZ/b3t6uiooKHTp0iD20hsH9GRn35sS4PyfG/RkZ9+bEMu3+GIahjo4OlZeXy+kceSROzrfoOJ1OTZ06NWWfX1xcnBH/YOzC/RkZ9+bEuD8nxv0ZGffmxDLp/pyoJSeOwcgAACBrEXQAAEDWIuikiM/n01133SWfz2d3KWmJ+zMy7s2JcX9OjPszMu7NiWXr/cn5wcgAACB70aIDAACyFkEHAABkLYIOAADIWgQdAACQtQg6KbJ27VrNnDlTfr9flZWVeuWVV+wuyRYvv/yyrrvuOpWXl8vhcOj3v//9gPOGYeh73/ueysvLlZeXpw9/+MN655137CnWYtXV1brkkktUVFSkiRMn6lOf+pT+8pe/DLgmV+/PunXrdMEFFyQWLlu4cKH++Mc/Js7n6n0ZTnV1tRwOh1atWpU4lsv353vf+54cDseAx6RJkxLnc/nexB05ckR/93d/p/Hjxys/P18XXXSRtm3bljifbfeIoJMCGzZs0KpVq3TnnXeqtrZWS5Ys0bJly1RXV2d3aZbr6urShRdeqIceemjY8/fee6/uv/9+PfTQQ9qyZYsmTZqkq6++OrEHWTbbtGmTbr31Vr3xxhuqqalRJBLR0qVL1dXVlbgmV+/P1KlTdc8992jr1q3aunWrrrzySn3yk59M/Mc2V+/LYFu2bNFjjz2mCy64YMDxXL8/5513nurr6xOPnTt3Js7l+r1paWnR4sWL5fF49Mc//lG7du3ST37yE40ZMyZxTdbdIwNJd+mllxorV64ccOycc84xvv3tb9tUUXqQZDz77LOJn2OxmDFp0iTjnnvuSRwLBAJGSUmJ8cgjj9hQob0aGxsNScamTZsMw+D+DDZ27FjjF7/4BfelV0dHhzFnzhyjpqbG+NCHPmR87WtfMwyDfzd33XWXceGFFw57LtfvjWEYxre+9S3jiiuuGPF8Nt4jWnSSLBQKadu2bVq6dOmA40uXLtXmzZttqio9HThwQA0NDQPulc/n04c+9KGcvFdtbW2SpHHjxkni/sRFo1E9/fTT6urq0sKFC7kvvW699VZde+21+uhHPzrgOPdH2rt3r8rLyzVz5kz9zd/8jfbv3y+JeyNJzz33nKqqqvS5z31OEydO1MUXX6zHH388cT4b7xFBJ8mampoUjUZVVlY24HhZWZkaGhpsqio9xe8H98rsE1+9erWuuOIKzZ8/XxL3Z+fOnSosLJTP59PKlSv17LPP6txzz835+yJJTz/9tLZv367q6uoh53L9/lx22WVav369nn/+eT3++ONqaGjQokWL1NzcnPP3RpL279+vdevWac6cOXr++ee1cuVK3X777Vq/fr2k7Pz3k/O7l6eKw+EY8LNhGEOOwcS9km677Ta99dZbevXVV4ecy9X7M3fuXO3YsUOtra363e9+pxtvvFGbNm1KnM/V+3Lo0CF97Wtf0wsvvCC/3z/idbl6f5YtW5Z4ff7552vhwoWaPXu2/vVf/1WXX365pNy9N5IUi8VUVVWlH/3oR5Kkiy++WO+8847WrVunFStWJK7LpntEi06SlZaWyuVyDUm+jY2NQxJyrovPhMj1e/XVr35Vzz33nF588UVNnTo1cTzX74/X69VZZ52lqqoqVVdX68ILL9TPfvaznL8v27ZtU2NjoyorK+V2u+V2u7Vp0yY9+OCDcrvdiXuQq/dnsIKCAp1//vnau3dvzv/bkaTJkyfr3HPPHXBs3rx5icky2XiPCDpJ5vV6VVlZqZqamgHHa2pqtGjRIpuqSk8zZ87UpEmTBtyrUCikTZs25cS9MgxDt912m5555hn96U9/0syZMwecz/X7M5hhGAoGgzl/X6666irt3LlTO3bsSDyqqqr0hS98QTt27NCsWbNy+v4MFgwGtXv3bk2ePDnn/+1I0uLFi4csY7Fnzx5Nnz5dUpb+d8euUdDZ7OmnnzY8Ho/xxBNPGLt27TJWrVplFBQUGO+//77dpVmuo6PDqK2tNWpraw1Jxv3332/U1tYaBw8eNAzDMO655x6jpKTEeOaZZ4ydO3caf/u3f2tMnjzZaG9vt7ny1Punf/ono6SkxHjppZeM+vr6xKO7uztxTa7enzVr1hgvv/yyceDAAeOtt94yvvOd7xhOp9N44YUXDMPI3fsykv6zrgwjt+/P17/+deOll14y9u/fb7zxxhvGJz7xCaOoqCjx399cvjeGYRhvvvmm4Xa7jX/+53829u7dazz11FNGfn6+8W//9m+Ja7LtHhF0UuThhx82pk+fbni9XmPBggWJKcO55sUXXzQkDXnceOONhmGYUxnvuusuY9KkSYbP5zP+6q/+yti5c6e9RVtkuPsiyfjlL3+ZuCZX789NN92U+P/PhAkTjKuuuioRcgwjd+/LSAYHnVy+P8uXLzcmT55seDweo7y83Pj0pz9tvPPOO4nzuXxv4v7whz8Y8+fPN3w+n3HOOecYjz322IDz2XaPHIZhGPa0JQEAAKQWY3QAAEDWIugAAICsRdABAABZi6ADAACyFkEHAABkLYIOAADIWgQdAACQtQg6AAAgaxF0AABA1iLoAMhKixYt0j/8wz/YXQYAmxF0AGSdWCymt956SwsWLLC7FAA2I+gAyDrvvvuuurq6CDoACDoAss/27dvldrt1wQUX2F0KAJsRdABkne3bt+vcc8+V3++3uxQANiPoAMg627dvp9sKgCSCDoAsYxiGduzYocrKSrtLAZAGCDoAssq+ffvU1tZGiw4ASZLb7gIAIJm2b98uSXK5XHr77bcTxz0ej+bOnWtXWQBsQtABkFVqa2slSZdffvmA45dffrlef/11O0oCYCOHYRiG3UUAAACkAmN0AABA1iLoAACArEXQAQAAWYugAwAAshZBBwAAZC2CDgAAyFoEHQAAkLUIOgAAIGsRdAAAQNYi6AAAgKxF0AEAAFnr/wdMokxeHzkKkQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "plt.plot(Ls, Ss,  color='#CC4F1B')\n",
    "plt.fill_between(Ls, Ss-Serr, Ss+Serr, facecolor='#FF9848')\n",
    "plt.xlabel('$L$'); plt.ylabel('$S_E(L)$');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c9a31e94-3e3c-4aa6-883a-34af55fb4412",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}