DMRG_TDVP.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "zel9fS3OWfAR"
      },
      "source": [
        "# DMRG TDVP(Time dependent variation principle)\n",
        "\n",
        "Edited by [Min Long](https://quantummc.xyz/members/min-long/), December.2023.\n",
        "\n",
        "This notebook contains the dynamic properties of Heisenberg chian in the formalism of matric product state(MPS).\n",
        "In the previous notebook, we evaluate the time evolution operator $e^{i\\hat{H}t}|\\psi\\rangle$, descretize it by Trottoer decomposition and perform the evolution per time step.\n",
        "\n",
        "However, there is another one-step solution which called TDVP that can evolve MPS the for large t efficiently without any trotter error(This is extremely efficient when study low temperature behavior of the system via imaginary time evolution)\n",
        "\n",
        "The idea is to construct a tangent space projector and solve the Schodinger equation in tangent space\n",
        "$$\n",
        "i\\hbar \\frac{\\partial |\\psi\\rangle }{\\partial t } = \\hat{P}\\hat{H}|\\psi\\rangle\n",
        "$$\n",
        "\n",
        "The initial state $|\\psi\\rangle$ is prepared and written in MPS form. I first calculate the ground state. Then introduce the concept of MPS manifold, then I show how to construct tangent space projector, finally I solve the Schodinger equation  \n",
        "\n",
        "First we prepare the groudstate, this process is copied from the file \"DMRG_MPS\"\n",
        "\n",
        "For the reference, reader can refer to the original paper of Jutho Haegeman [\\[ref\\]](https://arxiv.org/abs/1408.5056)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "lTBN_gKnWfAV"
      },
      "source": [
        "## Ground state calculation"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "metadata": {
        "id": "gsjvfg15-NEs",
        "outputId": "8611d4fd-d57e-41bc-b2b7-b898f3ae1b69",
        "colab": {
          "base_uri": "https://localhost:8080/"
        }
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Sweep 0 Sites 0,1    Energy -10.309016994375    States    2 Truncation   0.000000000000\n",
            "Sweep 0 Sites 1,2    Energy -10.439692620786    States    2 Truncation   0.000000000000\n",
            "Sweep 0 Sites 2,3    Energy -10.645800958670    States    2 Truncation   0.000000000000\n",
            "Sweep 0 Sites 3,4    Energy -10.810586239348    States    2 Truncation   0.000000000000\n",
            "Sweep 0 Sites 4,5    Energy -10.993669962857    States    2 Truncation   0.000000000000\n",
            "Sweep 0 Sites 5,6    Energy -11.167602338044    States    2 Truncation   0.000000000000\n",
            "Sweep 0 Sites 6,7    Energy -11.345754582728    States    2 Truncation   0.000000000000\n",
            "Sweep 0 Sites 7,8    Energy -11.521867353000    States    2 Truncation   0.000000000000\n",
            "Sweep 0 Sites 8,9    Energy -11.698938381096    States    2 Truncation   0.000000000000\n",
            "Sweep 0 Sites 9,10    Energy -11.875551828355    States    2 Truncation   0.000000000000\n",
            "Sweep 0 Sites 10,11    Energy -12.052381755193    States    2 Truncation   0.000000000000\n",
            "Sweep 0 Sites 11,12    Energy -12.229108727378    States    2 Truncation   0.000000000000\n",
            "Sweep 0 Sites 12,13    Energy -12.405884518897    States    2 Truncation   0.000000000000\n",
            "Sweep 0 Sites 13,14    Energy -12.582637122912    States    2 Truncation   0.000000000000\n",
            "Sweep 0 Sites 14,15    Energy -12.759400730100    States    2 Truncation   0.000000000000\n",
            "Sweep 0 Sites 15,16    Energy -12.936159113300    States    2 Truncation   0.000000000000\n",
            "Sweep 0 Sites 16,17    Energy -13.112919976004    States    2 Truncation   0.000000000000\n",
            "Sweep 0 Sites 17,18    Energy -13.289679661660    States    2 Truncation   0.000000000000\n",
            "Sweep 0 Sites 18,19    Energy -13.466439906026    States    2 Truncation   0.000000000000\n",
            "Sweep 0 Sites 19,20    Energy -13.643199885177    States    2 Truncation   0.000000000000\n",
            "Sweep 0 Sites 20,21    Energy -13.819959990220    States    2 Truncation   0.000000000000\n",
            "Sweep 0 Sites 21,22    Energy -13.996720035504    States    2 Truncation   0.000000000000\n",
            "Sweep 0 Sites 22,23    Energy -14.173480109155    States    2 Truncation   0.000000000000\n",
            "Sweep 0 Sites 23,24    Energy -14.350240169340    States    2 Truncation   0.000000000000\n",
            "Sweep 0 Sites 24,25    Energy -14.527000235917    States    2 Truncation   0.000000000000\n",
            "Sweep 0 Sites 25,26    Energy -14.703760299460    States    2 Truncation   0.000000000000\n",
            "Sweep 0 Sites 26,27    Energy -14.880520364443    States    2 Truncation   0.000000000000\n",
            "Sweep 0 Sites 27,28    Energy -15.057280428742    States    2 Truncation   0.000000000000\n",
            "Sweep 0 Sites 28,29    Energy -15.234040493366    States    2 Truncation   0.000000000000\n",
            "Sweep 0 Sites 29,30    Energy -15.410800557836    States    2 Truncation   0.000000000000\n",
            "Sweep 0 Sites 30,31    Energy -15.587560622380    States    2 Truncation   0.000000000000\n",
            "Sweep 0 Sites 31,32    Energy -15.764320686888    States    2 Truncation   0.000000000000\n",
            "Sweep 0 Sites 32,33    Energy -15.941080751413    States    2 Truncation   0.000000000000\n",
            "Sweep 0 Sites 33,34    Energy -16.117840815930    States    2 Truncation   0.000000000000\n",
            "Sweep 0 Sites 34,35    Energy -16.294600880451    States    2 Truncation   0.000000000000\n",
            "Sweep 0 Sites 35,36    Energy -16.471360944970    States    2 Truncation   0.000000000000\n",
            "Sweep 0 Sites 36,37    Energy -16.648121009490    States    2 Truncation   0.000000000000\n",
            "Sweep 0 Sites 37,38    Energy -16.824881074009    States    2 Truncation   0.000000000000\n",
            "Sweep 0 Sites 38,39    Energy -17.001641138529    States    2 Truncation   0.000000000000\n",
            "Sweep 1 Sites 39,40    Energy -17.266293521501    States    2 Truncation   0.000000000000\n",
            "Sweep 1 Sites 38,39    Energy -17.301814172211    States    4 Truncation   0.000000000000\n",
            "Sweep 1 Sites 37,38    Energy -17.306750465230    States    4 Truncation   0.000000000000\n",
            "Sweep 1 Sites 36,37    Energy -17.334080094395    States    4 Truncation   0.000000000000\n",
            "Sweep 1 Sites 35,36    Energy -17.341896222840    States    4 Truncation   0.000000000000\n",
            "Sweep 1 Sites 34,35    Energy -17.364006218378    States    4 Truncation   0.000000000000\n",
            "Sweep 1 Sites 33,34    Energy -17.373688405001    States    4 Truncation   0.000000000000\n",
            "Sweep 1 Sites 32,33    Energy -17.392530901625    States    4 Truncation   0.000000000000\n",
            "Sweep 1 Sites 31,32    Energy -17.403513611141    States    4 Truncation   0.000000000000\n",
            "Sweep 1 Sites 30,31    Energy -17.420359238293    States    4 Truncation   0.000000000000\n",
            "Sweep 1 Sites 29,30    Energy -17.432255166954    States    4 Truncation   0.000000000000\n",
            "Sweep 1 Sites 28,29    Energy -17.447893664001    States    4 Truncation   0.000000000000\n",
            "Sweep 1 Sites 27,28    Energy -17.460419892182    States    4 Truncation   0.000000000000\n",
            "Sweep 1 Sites 26,27    Energy -17.475328999411    States    4 Truncation   0.000000000000\n",
            "Sweep 1 Sites 25,26    Energy -17.488281105940    States    4 Truncation   0.000000000000\n",
            "Sweep 1 Sites 24,25    Energy -17.502747130697    States    4 Truncation   0.000000000000\n",
            "Sweep 1 Sites 23,24    Energy -17.515981721790    States    4 Truncation   0.000000000000\n",
            "Sweep 1 Sites 22,23    Energy -17.530176555869    States    4 Truncation   0.000000000000\n",
            "Sweep 1 Sites 21,22    Energy -17.543595855114    States    4 Truncation   0.000000000000\n",
            "Sweep 1 Sites 20,21    Energy -17.557623447949    States    4 Truncation   0.000000000000\n",
            "Sweep 1 Sites 19,20    Energy -17.571162272860    States    4 Truncation   0.000000000000\n",
            "Sweep 1 Sites 18,19    Energy -17.585086023543    States    4 Truncation   0.000000000000\n",
            "Sweep 1 Sites 17,18    Energy -17.598701692266    States    4 Truncation   0.000000000000\n",
            "Sweep 1 Sites 16,17    Energy -17.612560467573    States    4 Truncation   0.000000000000\n",
            "Sweep 1 Sites 15,16    Energy -17.626225629382    States    4 Truncation   0.000000000000\n",
            "Sweep 1 Sites 14,15    Energy -17.640042847141    States    4 Truncation   0.000000000000\n",
            "Sweep 1 Sites 13,14    Energy -17.653741326340    States    4 Truncation   0.000000000000\n",
            "Sweep 1 Sites 12,13    Energy -17.667528603075    States    4 Truncation   0.000000000000\n",
            "Sweep 1 Sites 11,12    Energy -17.681256339808    States    4 Truncation   0.000000000000\n",
            "Sweep 1 Sites 10,11    Energy -17.695007823588    States    4 Truncation   0.000000000000\n",
            "Sweep 1 Sites 9,10    Energy -17.708790453069    States    4 Truncation   0.000000000000\n",
            "Sweep 1 Sites 8,9    Energy -17.722445455229    States    4 Truncation   0.000000000000\n",
            "Sweep 1 Sites 7,8    Energy -17.736428925312    States    4 Truncation   0.000000000000\n",
            "Sweep 1 Sites 6,7    Energy -17.749711649837    States    4 Truncation   0.000000000000\n",
            "Sweep 1 Sites 5,6    Energy -17.764621497775    States    4 Truncation   0.000000000000\n",
            "Sweep 1 Sites 4,5    Energy -17.776482002484    States    4 Truncation   0.000000000000\n",
            "Sweep 1 Sites 3,4    Energy -17.796322537564    States    4 Truncation   0.000000000000\n",
            "Sweep 1 Sites 2,3    Energy -17.803801524498    States    4 Truncation   0.000000000000\n",
            "Sweep 1 Sites 1,2    Energy -17.853839396322    States    4 Truncation   0.000000000000\n",
            "Sweep 2 Sites 0,1    Energy -17.853839396322    States    2 Truncation   0.000000000000\n",
            "Sweep 2 Sites 1,2    Energy -17.853839396322    States    4 Truncation   0.000000000000\n",
            "Sweep 2 Sites 2,3    Energy -17.858309555838    States    8 Truncation   0.000000000000\n",
            "Sweep 2 Sites 3,4    Energy -17.858520176025    States    8 Truncation   0.000000000000\n",
            "Sweep 2 Sites 4,5    Energy -17.864087144428    States    8 Truncation   0.000000000000\n",
            "Sweep 2 Sites 5,6    Energy -17.864534762952    States    8 Truncation   0.000000000000\n",
            "Sweep 2 Sites 6,7    Energy -17.869775582058    States    8 Truncation   0.000000000000\n",
            "Sweep 2 Sites 7,8    Energy -17.870444067659    States    8 Truncation   0.000000000000\n",
            "Sweep 2 Sites 8,9    Energy -17.875155394193    States    8 Truncation   0.000000000000\n",
            "Sweep 2 Sites 9,10    Energy -17.876019891868    States    8 Truncation   0.000000000000\n",
            "Sweep 2 Sites 10,11    Energy -17.880237063587    States    8 Truncation   0.000000000000\n",
            "Sweep 2 Sites 11,12    Energy -17.881271937929    States    8 Truncation   0.000000000000\n",
            "Sweep 2 Sites 12,13    Energy -17.885067307994    States    8 Truncation   0.000000000000\n",
            "Sweep 2 Sites 13,14    Energy -17.886251189275    States    8 Truncation   0.000000000000\n",
            "Sweep 2 Sites 14,15    Energy -17.889694166204    States    8 Truncation   0.000000000000\n",
            "Sweep 2 Sites 15,16    Energy -17.891011030222    States    8 Truncation   0.000000000000\n",
            "Sweep 2 Sites 16,17    Energy -17.894159589101    States    8 Truncation   0.000000000000\n",
            "Sweep 2 Sites 17,18    Energy -17.895598807913    States    8 Truncation   0.000000000000\n",
            "Sweep 2 Sites 18,19    Energy -17.898497691009    States    8 Truncation   0.000000000000\n",
            "Sweep 2 Sites 19,20    Energy -17.900054685939    States    8 Truncation   0.000000000000\n",
            "Sweep 2 Sites 20,21    Energy -17.902734974166    States    8 Truncation   0.000000000000\n",
            "Sweep 2 Sites 21,22    Energy -17.904413266306    States    8 Truncation   0.000000000000\n",
            "Sweep 2 Sites 22,23    Energy -17.906891947440    States    8 Truncation   0.000000000000\n",
            "Sweep 2 Sites 23,24    Energy -17.908707494716    States    8 Truncation   0.000000000000\n",
            "Sweep 2 Sites 24,25    Energy -17.910986659633    States    8 Truncation   0.000000000000\n",
            "Sweep 2 Sites 25,26    Energy -17.912976194401    States    8 Truncation   0.000000000000\n",
            "Sweep 2 Sites 26,27    Energy -17.915041983535    States    8 Truncation   0.000000000000\n",
            "Sweep 2 Sites 27,28    Energy -17.917278850421    States    8 Truncation   0.000000000000\n",
            "Sweep 2 Sites 28,29    Energy -17.919100759473    States    8 Truncation   0.000000000000\n",
            "Sweep 2 Sites 29,30    Energy -17.921725085981    States    8 Truncation   0.000000000000\n",
            "Sweep 2 Sites 30,31    Energy -17.923256720645    States    8 Truncation   0.000000000000\n",
            "Sweep 2 Sites 31,32    Energy -17.926531243894    States    8 Truncation   0.000000000000\n",
            "Sweep 2 Sites 32,33    Energy -17.927713911453    States    8 Truncation   0.000000000000\n",
            "Sweep 2 Sites 33,34    Energy -17.932110942788    States    8 Truncation   0.000000000000\n",
            "Sweep 2 Sites 34,35    Energy -17.932882735006    States    8 Truncation   0.000000000000\n",
            "Sweep 2 Sites 35,36    Energy -17.939093633434    States    8 Truncation   0.000000000000\n",
            "Sweep 2 Sites 36,37    Energy -17.939419093715    States    8 Truncation   0.000000000000\n",
            "Sweep 2 Sites 37,38    Energy -17.947044669274    States    8 Truncation   0.000000000000\n",
            "Sweep 2 Sites 38,39    Energy -17.947044669274    States    4 Truncation   0.000000000000\n",
            "Sweep 3 Sites 39,40    Energy -17.947044669274    States    2 Truncation   0.000000000000\n",
            "Sweep 3 Sites 38,39    Energy -17.947044669274    States    4 Truncation   0.000000000000\n",
            "Sweep 3 Sites 37,38    Energy -17.947044669274    States    8 Truncation   0.000000000000\n",
            "Sweep 3 Sites 36,37    Energy -17.947253401667    States   16 Truncation   0.000000000000\n",
            "Sweep 3 Sites 35,36    Energy -17.947263693935    States   16 Truncation   0.000000000000\n",
            "Sweep 3 Sites 34,35    Energy -17.947784188082    States   16 Truncation   0.000000000000\n",
            "Sweep 3 Sites 33,34    Energy -17.947815875613    States   16 Truncation   0.000000000000\n",
            "Sweep 3 Sites 32,33    Energy -17.948568808001    States   16 Truncation   0.000000000000\n",
            "Sweep 3 Sites 31,32    Energy -17.948626868336    States   16 Truncation   0.000000000000\n",
            "Sweep 3 Sites 30,31    Energy -17.949500478419    States   16 Truncation   0.000000000000\n",
            "Sweep 3 Sites 29,30    Energy -17.949587242521    States   16 Truncation   0.000000000000\n",
            "Sweep 3 Sites 28,29    Energy -17.950492420481    States   16 Truncation   0.000000000000\n",
            "Sweep 3 Sites 27,28    Energy -17.950609935936    States   16 Truncation   0.000000000000\n",
            "Sweep 3 Sites 26,27    Energy -17.951488325298    States   16 Truncation   0.000000000000\n",
            "Sweep 3 Sites 25,26    Energy -17.951639507600    States   16 Truncation   0.000000000000\n",
            "Sweep 3 Sites 24,25    Energy -17.952456823562    States   16 Truncation   0.000000000000\n",
            "Sweep 3 Sites 23,24    Energy -17.952646077711    States   16 Truncation   0.000000000000\n",
            "Sweep 3 Sites 22,23    Energy -17.953383816001    States   16 Truncation   0.000000000000\n",
            "Sweep 3 Sites 21,22    Energy -17.953617538952    States   16 Truncation   0.000000000000\n",
            "Sweep 3 Sites 20,21    Energy -17.954266739212    States   16 Truncation   0.000000000000\n",
            "Sweep 3 Sites 19,20    Energy -17.954553886052    States   16 Truncation   0.000000000000\n",
            "Sweep 3 Sites 18,19    Energy -17.955111132773    States   16 Truncation   0.000000000000\n",
            "Sweep 3 Sites 17,18    Energy -17.955463874591    States   16 Truncation   0.000000000000\n",
            "Sweep 3 Sites 16,17    Energy -17.955928932433    States   16 Truncation   0.000000000000\n",
            "Sweep 3 Sites 15,16    Energy -17.956363250542    States   16 Truncation   0.000000000000\n",
            "Sweep 3 Sites 14,15    Energy -17.956737799146    States   16 Truncation   0.000000000000\n",
            "Sweep 3 Sites 13,14    Energy -17.957273454485    States   16 Truncation   0.000000000000\n",
            "Sweep 3 Sites 12,13    Energy -17.957560628031    States   16 Truncation   0.000000000000\n",
            "Sweep 3 Sites 11,12    Energy -17.958218845645    States   16 Truncation   0.000000000000\n",
            "Sweep 3 Sites 10,11    Energy -17.958423500257    States   16 Truncation   0.000000000000\n",
            "Sweep 3 Sites 9,10    Energy -17.959217671648    States   16 Truncation   0.000000000000\n",
            "Sweep 3 Sites 8,9    Energy -17.959347392745    States   16 Truncation   0.000000000000\n",
            "Sweep 3 Sites 7,8    Energy -17.960254403774    States   16 Truncation   0.000000000000\n",
            "Sweep 3 Sites 6,7    Energy -17.960321017789    States   16 Truncation   0.000000000000\n",
            "Sweep 3 Sites 5,6    Energy -17.961208819401    States   16 Truncation   0.000000000000\n",
            "Sweep 3 Sites 4,5    Energy -17.961229683952    States   16 Truncation   0.000000000000\n",
            "Sweep 3 Sites 3,4    Energy -17.961761029265    States   16 Truncation   0.000000000000\n",
            "Sweep 3 Sites 2,3    Energy -17.961761029265    States    8 Truncation   0.000000000000\n",
            "Sweep 3 Sites 1,2    Energy -17.961761029265    States    4 Truncation   0.000000000000\n",
            "Sweep 4 Sites 0,1    Energy -17.961761029265    States    2 Truncation   0.000000000000\n",
            "Sweep 4 Sites 1,2    Energy -17.961761029265    States    4 Truncation   0.000000000000\n",
            "Sweep 4 Sites 2,3    Energy -17.961761029265    States    8 Truncation   0.000000000000\n",
            "Sweep 4 Sites 3,4    Energy -17.961761029265    States   16 Truncation   0.000000000000\n",
            "Sweep 4 Sites 4,5    Energy -17.961765852380    States   32 Truncation   0.000000000000\n",
            "Sweep 4 Sites 5,6    Energy -17.961766370314    States   32 Truncation   0.000000000000\n",
            "Sweep 4 Sites 6,7    Energy -17.961785483666    States   32 Truncation   0.000000000000\n",
            "Sweep 4 Sites 7,8    Energy -17.961787751629    States   32 Truncation   0.000000000000\n",
            "Sweep 4 Sites 8,9    Energy -17.961826594252    States   32 Truncation   0.000000000000\n",
            "Sweep 4 Sites 9,10    Energy -17.961831998146    States   32 Truncation   0.000000000000\n",
            "Sweep 4 Sites 10,11    Energy -17.961890097781    States   32 Truncation   0.000000000000\n",
            "Sweep 4 Sites 11,12    Energy -17.961899990068    States   32 Truncation   0.000000000000\n",
            "Sweep 4 Sites 12,13    Energy -17.961972921196    States   32 Truncation   0.000000000000\n",
            "Sweep 4 Sites 13,14    Energy -17.961988607633    States   32 Truncation   0.000000000000\n",
            "Sweep 4 Sites 14,15    Energy -17.962070200406    States   32 Truncation   0.000000000000\n",
            "Sweep 4 Sites 15,16    Energy -17.962092954500    States   32 Truncation   0.000000000000\n",
            "Sweep 4 Sites 16,17    Energy -17.962176838784    States   32 Truncation   0.000000000000\n",
            "Sweep 4 Sites 17,18    Energy -17.962207900269    States   32 Truncation   0.000000000000\n",
            "Sweep 4 Sites 18,19    Energy -17.962288397447    States   32 Truncation   0.000000000000\n",
            "Sweep 4 Sites 19,20    Energy -17.962328919543    States   32 Truncation   0.000000000000\n",
            "Sweep 4 Sites 20,21    Energy -17.962401498816    States   32 Truncation   0.000000000000\n",
            "Sweep 4 Sites 21,22    Energy -17.962452389661    States   32 Truncation   0.000000000000\n",
            "Sweep 4 Sites 22,23    Energy -17.962513880647    States   32 Truncation   0.000000000000\n",
            "Sweep 4 Sites 23,24    Energy -17.962575453283    States   32 Truncation   0.000000000000\n",
            "Sweep 4 Sites 24,25    Energy -17.962624118911    States   32 Truncation   0.000000000000\n",
            "Sweep 4 Sites 25,26    Energy -17.962695432948    States   32 Truncation   0.000000000000\n",
            "Sweep 4 Sites 26,27    Energy -17.962730942051    States   32 Truncation   0.000000000000\n",
            "Sweep 4 Sites 27,28    Energy -17.962808740782    States   32 Truncation   0.000000000000\n",
            "Sweep 4 Sites 28,29    Energy -17.962832044271    States   32 Truncation   0.000000000000\n",
            "Sweep 4 Sites 29,30    Energy -17.962909430176    States   32 Truncation   0.000000000000\n",
            "Sweep 4 Sites 30,31    Energy -17.962922534556    States   32 Truncation   0.000000000000\n",
            "Sweep 4 Sites 31,32    Energy -17.962988351790    States   32 Truncation   0.000000000000\n",
            "Sweep 4 Sites 32,33    Energy -17.962994006293    States   32 Truncation   0.000000000000\n",
            "Sweep 4 Sites 33,34    Energy -17.963035646528    States   32 Truncation   0.000000000000\n",
            "Sweep 4 Sites 34,35    Energy -17.963036986246    States   32 Truncation   0.000000000000\n",
            "Sweep 4 Sites 35,36    Energy -17.963050313735    States   32 Truncation   0.000000000000\n",
            "Sweep 4 Sites 36,37    Energy -17.963050313735    States   16 Truncation   0.000000000000\n",
            "Sweep 4 Sites 37,38    Energy -17.963050313735    States    8 Truncation   0.000000000000\n",
            "Sweep 4 Sites 38,39    Energy -17.963050313735    States    4 Truncation   0.000000000000\n",
            "Sweep 5 Sites 39,40    Energy -17.963050313735    States    2 Truncation   0.000000000000\n",
            "Sweep 5 Sites 38,39    Energy -17.963050313735    States    4 Truncation   0.000000000000\n",
            "Sweep 5 Sites 37,38    Energy -17.963050313735    States    8 Truncation   0.000000000000\n",
            "Sweep 5 Sites 36,37    Energy -17.963050313735    States   16 Truncation   0.000000000000\n",
            "Sweep 5 Sites 35,36    Energy -17.963050313735    States   32 Truncation   0.000000000000\n",
            "Sweep 5 Sites 34,35    Energy -17.963050339451    States   64 Truncation   0.000000000000\n",
            "Sweep 5 Sites 33,34    Energy -17.963050364712    States   64 Truncation   0.000000000000\n",
            "Sweep 5 Sites 32,33    Energy -17.963050531355    States   64 Truncation   0.000000000000\n",
            "Sweep 5 Sites 31,32    Energy -17.963050670282    States   64 Truncation   0.000000000000\n",
            "Sweep 5 Sites 30,31    Energy -17.963051213924    States   64 Truncation   0.000000000000\n",
            "Sweep 5 Sites 29,30    Energy -17.963051595948    States   64 Truncation   0.000000000000\n",
            "Sweep 5 Sites 28,29    Energy -17.963052808328    States   64 Truncation   0.000000000000\n",
            "Sweep 5 Sites 27,28    Energy -17.963053572449    States   64 Truncation   0.000000000000\n",
            "Sweep 5 Sites 26,27    Energy -17.963055642757    States   64 Truncation   0.000000000000\n",
            "Sweep 5 Sites 25,26    Energy -17.963056924393    States   64 Truncation   0.000000000000\n",
            "Sweep 5 Sites 24,25    Energy -17.963059815491    States   64 Truncation   0.000000000000\n",
            "Sweep 5 Sites 23,24    Energy -17.963061730592    States   64 Truncation   0.000000000000\n",
            "Sweep 5 Sites 22,23    Energy -17.963065181877    States   64 Truncation   0.000000000000\n",
            "Sweep 5 Sites 21,22    Energy -17.963067794311    States   64 Truncation   0.000000000000\n",
            "Sweep 5 Sites 20,21    Energy -17.963071420366    States   64 Truncation   0.000000000000\n",
            "Sweep 5 Sites 19,20    Energy -17.963074693909    States   64 Truncation   0.000000000000\n",
            "Sweep 5 Sites 18,19    Energy -17.963078102691    States   64 Truncation   0.000000000000\n",
            "Sweep 5 Sites 17,18    Energy -17.963081839266    States   64 Truncation   0.000000000000\n",
            "Sweep 5 Sites 16,17    Energy -17.963084723224    States   64 Truncation   0.000000000000\n",
            "Sweep 5 Sites 15,16    Energy -17.963088509384    States   64 Truncation   0.000000000000\n",
            "Sweep 5 Sites 14,15    Energy -17.963090695669    States   64 Truncation   0.000000000000\n",
            "Sweep 5 Sites 13,14    Energy -17.963093927149    States   64 Truncation   0.000000000000\n",
            "Sweep 5 Sites 12,13    Energy -17.963095382839    States   64 Truncation   0.000000000000\n",
            "Sweep 5 Sites 11,12    Energy -17.963097525129    States   64 Truncation   0.000000000000\n",
            "Sweep 5 Sites 10,11    Energy -17.963098328841    States   64 Truncation   0.000000000000\n",
            "Sweep 5 Sites 9,10    Energy -17.963099337314    States   64 Truncation   0.000000000000\n",
            "Sweep 5 Sites 8,9    Energy -17.963099654844    States   64 Truncation   0.000000000000\n",
            "Sweep 5 Sites 7,8    Energy -17.963099962177    States   64 Truncation   0.000000000000\n",
            "Sweep 5 Sites 6,7    Energy -17.963100022420    States   64 Truncation   0.000000000000\n",
            "Sweep 5 Sites 5,6    Energy -17.963100068908    States   64 Truncation   0.000000000000\n",
            "Sweep 5 Sites 4,5    Energy -17.963100068908    States   32 Truncation   0.000000000000\n",
            "Sweep 5 Sites 3,4    Energy -17.963100068908    States   16 Truncation   0.000000000000\n",
            "Sweep 5 Sites 2,3    Energy -17.963100068908    States    8 Truncation   0.000000000000\n",
            "Sweep 5 Sites 1,2    Energy -17.963100068908    States    4 Truncation   0.000000000000\n",
            "Final energy expectation value -17.963100068907615\n",
            "variance =   0.000000408625\n"
          ]
        }
      ],
      "source": [
        "import numpy as np\n",
        "import scipy\n",
        "import scipy.sparse.linalg\n",
        "import scipy.sparse as sparse\n",
        "from time import perf_counter\n",
        "import numpy.linalg as npla\n",
        "import matplotlib.pyplot as plt\n",
        "import copy\n",
        "import math\n",
        "\n",
        "\n",
        "## initial E and F matrices for the left and right vacuum states\n",
        "def initial_E(W):\n",
        "    E = np.zeros((W.shape[0], 1, 1))\n",
        "    E[0] = 1\n",
        "    return E\n",
        "\n",
        "\n",
        "def initial_F(W):\n",
        "    F = np.zeros((W.shape[1], 1, 1))\n",
        "    F[-1] = 1\n",
        "    return F\n",
        "\n",
        "\n",
        "def contract_from_right(W, A, F, B):\n",
        "    # the einsum function doesn't appear to optimize the contractions properly,\n",
        "    # so we split it into individual summations in the optimal order\n",
        "    # return np.einsum(\"abst,sij,bjl,tkl->aik\",W,A,F,B, optimize=True)\n",
        "    Temp = np.einsum(\"sij,bjl->sbil\", A, F)\n",
        "    Temp = np.einsum(\"sbil,abst->tail\", Temp, W)\n",
        "    return np.einsum(\"tail,tkl->aik\", Temp, B.conjugate())\n",
        "\n",
        "\n",
        "def contract_from_left(W, A, E, B):\n",
        "    # the einsum function doesn't appear to optimize the contractions properly,\n",
        "    # so we split it into individual summations in the optimal order\n",
        "    # return np.einsum(\"abst,sij,aik,tkl->bjl\",W,A,E,B, optimize=True)\n",
        "    Temp = np.einsum(\"sij,aik->sajk\", A, E)\n",
        "    Temp = np.einsum(\"sajk,abst->tbjk\", Temp, W)\n",
        "    return np.einsum(\"tbjk,tkl->bjl\", Temp, B.conjugate())\n",
        "\n",
        "\n",
        "# construct the F-matrices for all sites except the first\n",
        "def construct_F(Alist, MPO, Blist):\n",
        "    F = [initial_F(MPO[-1])]\n",
        "\n",
        "    for i in range(len(MPO) - 1, 0, -1):\n",
        "        F.append(contract_from_right(MPO[i], Alist[i], F[-1], Blist[i]))\n",
        "    return F\n",
        "\n",
        "\n",
        "def construct_E(Alist, MPO, Blist):\n",
        "    return [initial_E(MPO[0])]\n",
        "\n",
        "\n",
        "# 2-1 coarse-graining of two site MPO into one site\n",
        "def coarse_grain_MPO(W, X):\n",
        "    return np.reshape(np.einsum(\"abst,bcuv->acsutv\", W, X),\n",
        "                      [W.shape[0], X.shape[1],\n",
        "                       W.shape[2] * X.shape[2],\n",
        "                       W.shape[3] * X.shape[3]])\n",
        "\n",
        "\n",
        "def product_W(W1, W2):\n",
        "    return np.reshape(np.einsum(\"abst,cdtu->acbdsu\", W1, W2), [W1.shape[0] * W2.shape[0],\n",
        "                                                               W1.shape[1] * W2.shape[1],\n",
        "                                                               W1.shape[2], W2.shape[3]])\n",
        "\n",
        "\n",
        "def product_MPO(M1, M2):\n",
        "    assert len(M1) == len(M2)\n",
        "    Result = []\n",
        "    for i in range(0, len(M1)):\n",
        "        Result.append(product_W(M1[i], M2[i]))\n",
        "    return Result\n",
        "\n",
        "\n",
        "# 2-1 coarse-graining of two-site MPS into one site\n",
        "def coarse_grain_MPS(A, B):\n",
        "    return np.reshape(np.einsum(\"sij,tjk->stik\", A, B),\n",
        "                      [A.shape[0] * B.shape[0], A.shape[1], B.shape[2]])\n",
        "\n",
        "\n",
        "def fine_grain_MPS(A, dims):\n",
        "    assert A.shape[0] == dims[0] * dims[1]\n",
        "    Theta = np.transpose(np.reshape(A, dims + [A.shape[1], A.shape[2]]),\n",
        "                         (0, 2, 1, 3))\n",
        "    M = np.reshape(Theta, (dims[0] * A.shape[1], dims[1] * A.shape[2]))\n",
        "    U, S, V = np.linalg.svd(M, full_matrices=0)\n",
        "    U = np.reshape(U, (dims[0], A.shape[1], -1))\n",
        "    V = np.transpose(np.reshape(V, (-1, dims[1], A.shape[2])), (1, 0, 2))\n",
        "    # assert U is left-orthogonal\n",
        "    # assert V is right-orthogonal\n",
        "    # print(np.dot(V[0],np.transpose(V[0])) + np.dot(V[1],np.transpose(V[1])))\n",
        "    return U, S, V\n",
        "\n",
        "\n",
        "def truncate_MPS(U, S, V, m):\n",
        "    m = min(len(S), m)\n",
        "    trunc = np.sum(S[m:])\n",
        "    S = S[0:m]\n",
        "    U = U[:, :, 0:m]\n",
        "    V = V[:, 0:m, :]\n",
        "    return U, S, V, trunc, m\n",
        "\n",
        "\n",
        "def Expectation(AList, MPO, BList):\n",
        "    E = [[[1]]]\n",
        "    for i in range(0, len(MPO)):\n",
        "        E = contract_from_left(MPO[i], AList[i], E, BList[i])\n",
        "    return E[0][0][0]\n",
        "\n",
        "\n",
        "class HamiltonianMultiply(sparse.linalg.LinearOperator):\n",
        "    def __init__(self, E, W, F):\n",
        "        self.E = E\n",
        "        self.W = W\n",
        "        self.F = F\n",
        "        self.dtype = np.dtype('d')\n",
        "        self.req_shape = [W.shape[2], E.shape[1], F.shape[1]]\n",
        "        self.size = self.req_shape[0] * self.req_shape[1] * self.req_shape[2]\n",
        "        self.shape = [self.size, self.size]\n",
        "\n",
        "    def _matvec(self, A):\n",
        "        # the einsum function doesn't appear to optimize the contractions properly,\n",
        "        # so we split it into individual summations in the optimal order\n",
        "        # R = np.einsum(\"aij,sik,abst,bkl->tjl\",self.E,np.reshape(A, self.req_shape),\n",
        "        #              self.W,self.F, optimize=True)\n",
        "        R = np.einsum(\"aij,sik->ajsk\", self.E, np.reshape(A, self.req_shape))\n",
        "        R = np.einsum(\"ajsk,abst->bjtk\", R, self.W)\n",
        "        R = np.einsum(\"bjtk,bkl->tjl\", R, self.F)\n",
        "        return np.reshape(R, -1)\n",
        "\n",
        "\n",
        "## optimize a single site given the MPO matrix W, and tensors E,F\n",
        "def optimize_site(A, W, E, F):\n",
        "    H = HamiltonianMultiply(E, W, F)\n",
        "    # we choose tol=1E-8 here, which is OK for small calculations.\n",
        "    # to bemore robust, we should take the tol -> 0 towards the end\n",
        "    # of the calculation.\n",
        "    E, V = sparse.linalg.eigsh(H, 1, v0=A, which='SA', tol=1E-8)\n",
        "    return (E[0], np.reshape(V[:, 0], H.req_shape))\n",
        "\n",
        "\n",
        "def optimize_two_sites(A, B, W1, W2, E, F, m, dir):\n",
        "    W = coarse_grain_MPO(W1, W2)\n",
        "    AA = coarse_grain_MPS(A, B)\n",
        "    H = HamiltonianMultiply(E, W, F)\n",
        "    E, V = sparse.linalg.eigsh(H, 1, v0=AA, which='SA')\n",
        "    AA = np.reshape(V[:, 0], H.req_shape)\n",
        "    A, S, B = fine_grain_MPS(AA, [A.shape[0], B.shape[0]])\n",
        "    A, S, B, trunc, m = truncate_MPS(A, S, B, m)\n",
        "    if (dir == 'right'):\n",
        "        B = np.einsum(\"ij,sjk->sik\", np.diag(S), B)\n",
        "    else:\n",
        "        assert dir == 'left'\n",
        "        A = np.einsum(\"sij,jk->sik\", A, np.diag(S))\n",
        "    return E[0], A, B, trunc, m\n",
        "def two_site_dmrg(MPS, MPO, m, sweeps):\n",
        "    E = construct_E(MPS, MPO, MPS)\n",
        "    F = construct_F(MPS, MPO, MPS)\n",
        "    F.pop()\n",
        "    for sweep in range(0, int(sweeps / 2)):\n",
        "        for i in range(0, len(MPS) - 2):\n",
        "            Energy, MPS[i], MPS[i + 1], trunc, states = optimize_two_sites(MPS[i], MPS[i + 1],\n",
        "                                                                           MPO[i], MPO[i + 1],\n",
        "                                                                           E[-1], F[-1], m, 'right')\n",
        "            print(\"Sweep {:} Sites {:},{:}    Energy {:16.12f}    States {:4} Truncation {:16.12f}\"\n",
        "                  .format(sweep * 2, i, i + 1, Energy, states, trunc))\n",
        "            E.append(contract_from_left(MPO[i], MPS[i], E[-1], MPS[i]))\n",
        "            F.pop();\n",
        "        for i in range(len(MPS) - 2, 0, -1):\n",
        "            Energy, MPS[i], MPS[i + 1], trunc, states = optimize_two_sites(MPS[i], MPS[i + 1],\n",
        "                                                                           MPO[i], MPO[i + 1],\n",
        "                                                                           E[-1], F[-1], m, 'left')\n",
        "            print(\"Sweep {} Sites {},{}    Energy {:16.12f}    States {:4} Truncation {:16.12f}\"\n",
        "\n",
        "                  .format(sweep * 2 + 1, i, i + 1, Energy, states, trunc))\n",
        "            F.append(contract_from_right(MPO[i + 1], MPS[i + 1], F[-1], MPS[i + 1]))\n",
        "            E.pop();\n",
        "    return MPS\n",
        "\n",
        "\n",
        "\n",
        "\n",
        "d = 2  # local bond dimension\n",
        "N = 41  # number of sites\n",
        "\n",
        "InitialA1 = np.zeros((d, 1, 1))\n",
        "InitialA1[0, 0, 0] = 1\n",
        "InitialA2 = np.zeros((d, 1, 1))\n",
        "InitialA2[1, 0, 0] = 1\n",
        "\n",
        "## initial state |01010101>\n",
        "MPS = [InitialA1, InitialA2] * int(N / 2)\n",
        "if N%2 != 0 :\n",
        "    MPS = MPS + [InitialA1]\n",
        "\n",
        "## Local operators\n",
        "I = np.identity(2)\n",
        "Z = np.zeros((2, 2))\n",
        "Sz = np.array([[0.5, 0],\n",
        "               [0, -0.5]])\n",
        "Sp = np.array([[0, 0],\n",
        "               [1, 0]])\n",
        "Sm = np.array([[0, 1],\n",
        "               [0, 0]])\n",
        "\n",
        "## Hamiltonian MPO\n",
        "W = np.array([[I, Sz, 0.5 * Sp, 0.5 * Sm, Z],\n",
        "              [Z, Z, Z, Z, Sz],\n",
        "              [Z, Z, Z, Z, Sm],\n",
        "              [Z, Z, Z, Z, Sp],\n",
        "              [Z, Z, Z, Z, I]])\n",
        "\n",
        "Wfirst = np.array([[I, Sz, 0.5 * Sp, 0.5 * Sm, Z]])\n",
        "\n",
        "Wlast = np.array([[Z], [Sz], [Sm], [Sp], [I]])\n",
        "\n",
        "# the complete MPO\n",
        "MPO = [Wfirst] + ([W] * (N - 2)) + [Wlast]\n",
        "\n",
        "HamSquared = product_MPO(MPO, MPO)\n",
        "\n",
        "D = 20\n",
        "# 6 sweeps with m=10 states\n",
        "MPS = two_site_dmrg(MPS, MPO, D, 6)\n",
        "\n",
        "Energy = Expectation(MPS, MPO, MPS)\n",
        "print(\"Final energy expectation value {}\".format(Energy))\n",
        "\n",
        "H2 = Expectation(MPS, HamSquared, MPS)\n",
        "print(\"variance = {:16.12f}\".format(H2 - Energy * Energy))\n",
        "\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "KT5EKN1X-NEt"
      },
      "source": [
        "## Measurement, transfer of singular value and benchmark"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 2,
      "metadata": {
        "id": "gDIFG0o6-NEt",
        "outputId": "1c606e3b-ee5b-4b5f-d02f-775a7631aa4a",
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 356
        }
      },
      "outputs": [
        {
          "output_type": "error",
          "ename": "KeyboardInterrupt",
          "evalue": "",
          "traceback": [
            "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
            "\u001b[0;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
            "\u001b[0;32m<ipython-input-2-74ee507341f1>\u001b[0m in \u001b[0;36m<cell line: 44>\u001b[0;34m()\u001b[0m\n\u001b[1;32m     42\u001b[0m     \u001b[0;31m# plt.savefig(\"Sz_correlation_measurement_tdvp\")\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     43\u001b[0m     \u001b[0mplt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mclose\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 44\u001b[0;31m \u001b[0mplot_Sz_correlaion_MPS\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mGs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     45\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     46\u001b[0m \u001b[0;31m# quit()\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
            "\u001b[0;32m<ipython-input-2-74ee507341f1>\u001b[0m in \u001b[0;36mplot_Sz_correlaion_MPS\u001b[0;34m(MPS_in)\u001b[0m\n\u001b[1;32m     32\u001b[0m     \u001b[0;32mfor\u001b[0m \u001b[0mi\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mN\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     33\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mi\u001b[0m \u001b[0;34m!=\u001b[0m \u001b[0mN\u001b[0m\u001b[0;34m//\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 34\u001b[0;31m             \u001b[0mSz_correlation\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmulti_site_measurement\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mMPS_in\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mGs\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mSz\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mSz\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mN\u001b[0m\u001b[0;34m//\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     35\u001b[0m         \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     36\u001b[0m             \u001b[0mSz_correlation\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmulti_site_measurement\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mMPS_in\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mGs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0mSz\u001b[0m\u001b[0;34m@\u001b[0m\u001b[0mSz\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
            "\u001b[0;32m<ipython-input-2-74ee507341f1>\u001b[0m in \u001b[0;36mmulti_site_measurement\u001b[0;34m(MPS1, MPS2, op_list, pos_list)\u001b[0m\n\u001b[1;32m     17\u001b[0m         \u001b[0mexpectation\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0meinsum\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"ik,sij->skj\"\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mexpectation\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mMPS1\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mj\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     18\u001b[0m         \u001b[0mexpectation\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0meinsum\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"skj,skl->jl\"\u001b[0m\u001b[0;34m,\u001b[0m  \u001b[0mexpectation\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mMPS2\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mj\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mconjugate\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 19\u001b[0;31m     \u001b[0mexpectation\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0meinsum\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'ik,sij,st,tkl->jl'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mexpectation\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mMPS1\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mpos_list\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mop_list\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mMPS2\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mpos_list\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mconjugate\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     20\u001b[0m     \u001b[0;32mfor\u001b[0m \u001b[0mi\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mop_list\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     21\u001b[0m         \u001b[0;32mfor\u001b[0m \u001b[0mj\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mrange\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mpos_list\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0mpos_list\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m+\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
            "\u001b[0;32m/usr/local/lib/python3.10/dist-packages/numpy/core/einsumfunc.py\u001b[0m in \u001b[0;36meinsum\u001b[0;34m(out, optimize, *operands, **kwargs)\u001b[0m\n\u001b[1;32m   1369\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mspecified_out\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1370\u001b[0m             \u001b[0mkwargs\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m'out'\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mout\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1371\u001b[0;31m         \u001b[0;32mreturn\u001b[0m \u001b[0mc_einsum\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0moperands\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   1372\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1373\u001b[0m     \u001b[0;31m# Check the kwargs to avoid a more cryptic error later, without having to\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
            "\u001b[0;31mKeyboardInterrupt\u001b[0m: "
          ]
        }
      ],
      "source": [
        "Gs = copy.deepcopy(MPS)\n",
        "\n",
        "def state_overlap(Co_List,Contro_List):\n",
        "    R = [[1]]\n",
        "    for i in range(0, len(Co_List)):\n",
        "        R = np.einsum(\"ik,sij->skj\",R,Co_List[i])\n",
        "        R = np.einsum(\"skj,skl->jl\",R,Contro_List[i].conjugate())\n",
        "    return R[0][0]\n",
        "def multi_site_measurement(MPS1,MPS2,op_list,pos_list):\n",
        "    idx = np.argsort(pos_list)\n",
        "    pos_list = np.array(pos_list)[idx]\n",
        "    op_list = np.array(op_list)[idx]\n",
        "    expectation = [[1]]\n",
        "    # print(op_list,pos_list)\n",
        "    assert len(op_list) == len(pos_list)\n",
        "    for j in range(0, pos_list[0]):\n",
        "        expectation = np.einsum(\"ik,sij->skj\", expectation, MPS1[j])\n",
        "        expectation = np.einsum(\"skj,skl->jl\",  expectation,MPS2[j].conjugate())\n",
        "    expectation = np.einsum('ik,sij,st,tkl->jl', expectation, MPS1[pos_list[0]], op_list[0], MPS2[pos_list[0]].conjugate())\n",
        "    for i in range(len(op_list)-1):\n",
        "        for j in range(pos_list[i]+1,pos_list[i+1]):\n",
        "            expectation = np.einsum(\"ik,sij->skj\", expectation, MPS1[j])\n",
        "            expectation = np.einsum(\"skj,skl->jl\", expectation,MPS2[j].conjugate())\n",
        "        expectation = np.einsum('ik,sij,st,tkl->jl', expectation, MPS1[pos_list[i+1]], op_list[i+1], MPS2[pos_list[i+1]].conjugate())\n",
        "    for j in range(pos_list[-1]+1,len(MPS1)):\n",
        "        expectation = np.einsum(\"ik,sij->skj\", expectation, MPS1[j])\n",
        "        expectation = np.einsum(\"skj,skl->jl\", expectation,MPS2[j].conjugate())\n",
        "    return expectation[0][0]\n",
        "\n",
        "def plot_Sz_correlaion_MPS(MPS_in):\n",
        "    Sz_correlation = []\n",
        "    for i in range(N):\n",
        "        if i != N//2:\n",
        "            Sz_correlation.append(multi_site_measurement(MPS_in,Gs,[Sz,Sz],[i,N//2]))\n",
        "        else:\n",
        "            Sz_correlation.append(multi_site_measurement(MPS_in, Gs, [Sz@Sz], [i]))\n",
        "    plt.plot(Sz_correlation)\n",
        "    plt.title(r\"$<\\hat{S}_z(x)\\hat{S}_z(0)>$\")\n",
        "    plt.xlabel(\"Site index\")\n",
        "    # plt.savefig(\"temp\")\n",
        "    plt.show()\n",
        "    # plt.savefig(\"Sz_correlation_measurement_tdvp\")\n",
        "    plt.close()\n",
        "plot_Sz_correlaion_MPS(Gs)\n",
        "\n",
        "# quit()\n",
        "\n",
        "def move_mix_site(start ,end , mix_MPS ,bond_dim):\n",
        "    if end > start:\n",
        "        move_dir = 'r'\n",
        "        for i in range(start,end,1):\n",
        "            M1 = mix_MPS[i]\n",
        "            M_mat = np.reshape(M1,(M1.shape[0]*M1.shape[1],M1.shape[2]))\n",
        "            A, S, B = np.linalg.svd(M_mat, full_matrices=0)\n",
        "            trunc_dim = min(bond_dim,len(S))\n",
        "            A = A[:,0:trunc_dim]\n",
        "            S = S[0:trunc_dim]\n",
        "            B = B[0:trunc_dim,:]\n",
        "            # print(npla.norm(A@np.diag(S)@B-M_mat))\n",
        "            trunc = np.sum(S[trunc_dim:])\n",
        "            A = np.reshape(A,(M1.shape[0],M1.shape[1],-1))\n",
        "            # print(\"Truncation error in moving from site {} to {} : {}\".format(i,i+1,trunc))\n",
        "            if trunc > 1e-10:\n",
        "                raise Exception(\"large truncation error\")\n",
        "            B = np.diag(S)@B\n",
        "            mix_MPS[i] =  A\n",
        "            mix_MPS[i+1] = np.einsum(\"ij,sjk->sik\",B,mix_MPS[i+1])\n",
        "    if end < start:\n",
        "        move_dir = 'l'\n",
        "        for i in range(start,end,-1):\n",
        "            M1 = mix_MPS[i]\n",
        "            M_mat = np.reshape(M1.transpose(1,0,2), (M1.shape[1], M1.shape[0] * M1.shape[2]))\n",
        "            A, S, B = np.linalg.svd(M_mat, full_matrices=0)\n",
        "            trunc_dim = min(bond_dim, len(S))\n",
        "            A = A[:, 0:trunc_dim]\n",
        "            S = S[0:trunc_dim]\n",
        "            B = B[0:trunc_dim, :]\n",
        "            # print(npla.norm(A @ np.diag(S) @ B - M_mat))\n",
        "            trunc = np.sum(S[trunc_dim:])\n",
        "            B = np.reshape(B, (-1,M1.shape[0], M1.shape[2])).transpose(1,0,2)\n",
        "            # print(\"Truncation error in moving from site {} to {} : {}\".format(i, i + 1, trunc))\n",
        "            if trunc > 1e-10:\n",
        "                raise Exception(\"large truncation error\")\n",
        "            A = A@np.diag(S)\n",
        "            mix_MPS[i] = B\n",
        "            mix_MPS[i - 1] = np.einsum(\"sjk,kl->sjl\", mix_MPS[i-1], A)\n",
        "    if \"move_dir\" not in locals():\n",
        "        raise Exception(\"Not standard input form, Nothing happend\")\n",
        "    return mix_MPS\n",
        "\n",
        "def check_othogonal(MPS):\n",
        "    infor_string = []\n",
        "    for ten in MPS:\n",
        "        matl = np.einsum(\"sij,sik->jk\",ten,ten)\n",
        "        matr = np.einsum(\"sij,slj->il\",ten,ten)\n",
        "        dim = matl.shape[0]\n",
        "        Il = np.identity(dim)\n",
        "        Ir = np.identity(matr.shape[0])\n",
        "        if npla.norm(matl-Il)<1e-10:\n",
        "            infor_string.append('A')\n",
        "        if npla.norm(matr-Ir)<1e-10:\n",
        "            infor_string.append('B')\n",
        "        if not(npla.norm(matl-Il)<1e-10 or npla.norm(matr-Ir)<1e-10):\n",
        "            infor_string.append('M')\n",
        "    return ' '.join(infor_string)\n",
        "\n",
        "\n",
        "print(\"Normalizatin check of Ground state\", state_overlap(Gs, Gs))\n",
        "\n",
        "tau = 0.2  # 0.02\n",
        "time_slice = 60  # 600\n",
        "Evolution_Sz = np.zeros((time_slice + 1, N), dtype=complex)\n",
        "Correlaion_SpSm = np.zeros((time_slice + 1, N), dtype=complex)\n",
        "\n",
        "MPSt = copy.deepcopy(Gs)\n",
        "\n",
        "\n",
        "def orthogonize(M, bond_dim):\n",
        "    M1 = copy.deepcopy(M)\n",
        "    M_mat = np.reshape(M1, (M1.shape[0] * M1.shape[1], M1.shape[2]))\n",
        "    A, S, B = np.linalg.svd(M_mat, full_matrices=0)\n",
        "    trunc_dim = min(bond_dim, len(S))\n",
        "    A = A[:, 0:trunc_dim]\n",
        "    S = S[0:trunc_dim]\n",
        "    B = B[0:trunc_dim, :]\n",
        "    S = np.array(S) / np.sqrt(sum([abs(item) ** 2 for item in S]))\n",
        "    M_mat = np.reshape(A @ np.diag(S) @ B, (M1.shape[0], M1.shape[1], M1.shape[2]))\n",
        "    return M_mat\n",
        "\n",
        "\n",
        "MPSt = move_mix_site(start=1, end=N // 2, mix_MPS=MPSt, bond_dim=D)\n",
        "MPSt[N // 2] = np.einsum(\"sij,st->tij\", MPSt[N // 2], Sp)\n",
        "MPSt[N // 2] = orthogonize(MPSt[N // 2], D)\n",
        "MPSt = move_mix_site(start=N // 2, end=0, mix_MPS=MPSt, bond_dim=D)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "LbuRCO_p-NEt"
      },
      "source": [
        "## Tangent space technique for time evolution\n",
        "\n",
        "The method of performing variation on tangent space can be traced to Dirac. Geomatrically, the TDVP corresponds to an orthogonal projection of the evolution vector of the Schrodinger equation $-iH|\\psi\\rangle$ onto the tangent space of the MPS manifold.\n",
        "\n",
        "$$\n",
        "\\frac{d|\\Psi(A)\\rangle}{dt} = -i\\hat{P}_{T_{|\\Psi(A)\\rangle}\\mathcal{M}_{MPS}}\n",
        "$$\n",
        "\n",
        "To reveal the structure of the tangent space $T_{|\\Psi(A)\\rangle}\\mathcal{M}_{MPS}$, we expand the left side equation using MPS ansatz;\n",
        "\n",
        "$$\n",
        "\\frac{d|\\Psi(A)\\rangle}{dt} = B^i\\left|\\partial_i \\psi\\right\\rangle=\\sum_{n=1}^N \\sum_{\\left\\{s_n\\right\\}=1}^d A^{s_1}(1) \\cdots B^{s_n}(n) \\cdots A^{s_N}(N)\\left|s_1 \\ldots s_n \\ldots s_N\\right\\rangle\n",
        "$$\n",
        "\n",
        "\n",
        "Where we introduce the tangent vector $B^i = \\dot{A^i}$.\n",
        "\n",
        "Now it is clear that the LHS lives in the tangent space of $|\\Psi(A)\\rangle$.\n",
        "\n",
        "\n",
        "First, to find $|\\Theta[B]\\rangle$, we need to eliminate the gauge redundancy of $B$, this redanduncy can be understood from the gauge redanduncy of MPS. We fix the redundance by constraint:\n",
        "\n",
        "$$\n",
        "\\sum_s A^s_L(n)^\\dagger B^s(n) = 0\n",
        "$$\n",
        "\n",
        "We minimize the distance between the RHS and LHS to determine the value of $B^i$. In other word, the best approximation $|\\Theta[B]\\rangle$ of a vector $|\\Xi\\rangle$ out of the manifold satisfying:\n",
        "$$\n",
        "\\min|| \\Theta[B]\\rangle - |\\Xi\\rangle ||^2\n",
        "$$\n",
        "\n",
        "We perform minimization under the constraint imposed by gauge(after some manipulation........), finally we derive:\n",
        "\n",
        "$$\n",
        "\\begin{aligned}\n",
        "|\\Theta[B]\\rangle & =\\sum_{n=1}^N \\sum_{\\left\\{s_n\\right\\}=1}^d A_L^{s_1}(1) \\cdots A_L^{s_{n-1}}(n-1) B^{s_n}(n) A_R^{s_{n+1}}(n+1) \\cdots A_R^{s_N}(N)\\left|s_1 \\ldots s_n \\ldots s_N\\right\\rangle . \\\\\n",
        "& =\\sum_{n=1}^N \\sum_{\\left\\{\\alpha, \\beta, s_n\\right\\}} B_{\\alpha, \\beta}^{s_n}(n)\\left|\\Phi_{L, \\alpha}^{[1: n-1]}\\right\\rangle\\left|s_n\\right\\rangle\\left|\\Phi_{R, \\beta}^{[n+1: N]}\\right\\rangle\n",
        "\\end{aligned}\n",
        "$$\n",
        "\n",
        "$$\n",
        "|\\Theta[B]\\rangle=P_{T_{|\\Psi[A]\\rangle \\mathcal{M}}^{\\mathrm{MPS}}}|\\Xi\\rangle=\\left(\\sum_{n=1}^N \\hat{P}_L^{[1: n-1]} \\otimes \\hat{1}_n \\otimes \\hat{P}_R^{[n+1: N]}-\\sum_{n=1}^{N-1} \\hat{P}_L^{[1: n]} \\otimes \\hat{P}_R^{[n+1: N]}\\right)|\\Xi\\rangle\n",
        "$$\n",
        "\n",
        "For the computational scheme, please refer to the appendix of [\\[ref\\]](https://arxiv.org/abs/1408.5056)."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "UimBUZcq-NEu"
      },
      "outputs": [],
      "source": [
        "\n",
        "\n",
        "\n",
        "\n",
        "print('############################')\n",
        "\n",
        "start = perf_counter()\n",
        "for i in range(N):\n",
        "    Evolution_Sz[0, i] = multi_site_measurement(MPSt, MPSt, [Sz], [i])\n",
        "    Correlaion_SpSm[0, i] = multi_site_measurement(MPSt, Gs, [Sm], [i])\n",
        "\n",
        "\n",
        "def to_matrix_Hamiltonian(E, W, F):\n",
        "    H_tensor = np.einsum('aik,abst,bjl->sijtkl', E, W, F)\n",
        "    shape1 = W.shape[2] * E.shape[1] * F.shape[1]\n",
        "    shape2 = W.shape[3] * E.shape[2] * F.shape[2]\n",
        "    return np.reshape(H_tensor, [shape1, shape2])\n",
        "\n",
        "E = construct_E(MPSt, MPO, MPSt)\n",
        "F = construct_F(MPSt, MPO, MPSt)\n",
        "\n",
        "\n",
        "for t in range(time_slice // 2):\n",
        "    print('sweep: {}'.format(t))\n",
        "    F.pop()\n",
        "    for i in range(0, N - 1, 1):\n",
        "        W = coarse_grain_MPO(MPO[i], MPO[i + 1])\n",
        "        AA = coarse_grain_MPS(MPSt[i], MPSt[i + 1]).reshape(-1)\n",
        "        H = to_matrix_Hamiltonian(E=E[-1], W=W, F=F[-1])\n",
        "        AA = scipy.linalg.expm(-1j * H * tau) @ AA\n",
        "        AA = AA.reshape(MPSt[i].shape[0] * MPSt[i + 1].shape[0], MPSt[i].shape[1], MPSt[i + 1].shape[2])\n",
        "        A, S, B = fine_grain_MPS(AA, [MPSt[i].shape[0], MPSt[i + 1].shape[0]])\n",
        "        A, S, B, trunc, m = truncate_MPS(A, S, B, D)\n",
        "        # print(trunc)\n",
        "        B = np.einsum(\"ij,sjk->sik\", np.diag(S), B)\n",
        "        MPSt[i], MPSt[i + 1] = A, B\n",
        "        E.append(contract_from_left(MPO[i], MPSt[i], E[-1], MPSt[i]))\n",
        "        if i == N - 2:\n",
        "            break\n",
        "        K = to_matrix_Hamiltonian(E=E[-1], W=MPO[i + 1], F=F[-1])\n",
        "        MPSt[i + 1] = np.reshape(scipy.linalg.expm(1j * K * tau) @ MPSt[i + 1].reshape(-1),[MPSt[i + 1].shape[0],MPSt[i + 1].shape[1],MPSt[i + 1].shape[2]])\n",
        "        F.pop()\n",
        "    for i in range(N):\n",
        "        Evolution_Sz[t * 2 + 1, i] = multi_site_measurement(MPSt, MPSt, [Sz], [i])\n",
        "        Correlaion_SpSm[t * 2 + 1, i] = multi_site_measurement(MPSt, Gs, [Sm], [i])\n",
        "    E.pop()\n",
        "    for i in range(N - 2, -1, -1):\n",
        "        W = coarse_grain_MPO(MPO[i], MPO[i + 1])\n",
        "        AA = coarse_grain_MPS(MPSt[i], MPSt[i + 1]).reshape(-1)\n",
        "        H = to_matrix_Hamiltonian(E=E[-1], W=W, F=F[-1])\n",
        "        AA = scipy.linalg.expm(-1j * H * tau) @ AA\n",
        "        AA = AA.reshape(MPSt[i].shape[0] * MPSt[i + 1].shape[0], MPSt[i].shape[1], MPSt[i + 1].shape[2])\n",
        "        A, S, B = fine_grain_MPS(AA, [MPSt[i].shape[0], MPSt[i + 1].shape[0]])\n",
        "        A, S, B, trunc, m = truncate_MPS(A, S, B, D)\n",
        "        # print(trunc)\n",
        "        A = np.einsum(\"sij,jk->sik\", A, np.diag(S))\n",
        "        MPSt[i], MPSt[i + 1] = A, B\n",
        "        F.append(contract_from_right(MPO[i + 1], MPSt[i + 1], F[-1], MPSt[i + 1]))\n",
        "        if i == 0:\n",
        "            break\n",
        "        K = to_matrix_Hamiltonian(E=E[-1], W=MPO[i], F=F[-1])\n",
        "        MPSt[i] = np.reshape(scipy.linalg.expm(1j * K * tau) @ MPSt[i].reshape(-1),[MPSt[i].shape[0],MPSt[i].shape[1],MPSt[i].shape[2]])\n",
        "        E.pop()\n",
        "    for i in range(N):\n",
        "        Evolution_Sz[2 * t + 2, i] = multi_site_measurement(MPSt, MPSt, [Sz], [i])\n",
        "        Correlaion_SpSm[2 * t + 2, i] = multi_site_measurement(MPSt, Gs, [Sm], [i])\n",
        "end = perf_counter()\n",
        "print('time used in time evolution: {:.2f}s'.format(end - start))\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "kl0Ug81pWfAc"
      },
      "source": [
        "## Evolution of spin wave"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "qOuYlSEPWfAd"
      },
      "source": [
        "Time evolution of Ground state after adding a excitation. In another word,\n",
        "$$\n",
        "|\\Psi\\rangle = S^+(\\frac{L}{2})|Gs\\rangle\n",
        "$$\n",
        "$|\\Psi(t) = e^{i\\hbar\\hat{H}}|\\Psi\\rangle$, we want to measure:\n",
        "$$\n",
        "S^z(t,x) = \\langle \\Psi(t)| S^z(x)|\\Psi(t)\\rangle\n",
        "$$\n",
        "We can see the linear dispersion of the spin wave."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "Md6603fJ-NEu"
      },
      "outputs": [],
      "source": [
        "\n",
        "# plt.plot(Evolution_Sz[-1,:])\n",
        "# plt.savefig('check')\n",
        "# quit()\n",
        "X = np.outer(np.linspace(0, N - 1, N), np.ones(time_slice + 1))\n",
        "T = np.outer(np.ones(N), np.linspace(0, tau, time_slice + 1))\n",
        "SzMatrix = np.transpose(np.array(Evolution_Sz).real)\n",
        "fig = plt.figure()\n",
        "ax = plt.axes(projection='3d')\n",
        "ax.plot_surface(X, T, SzMatrix, cmap='viridis', edgecolor='none')\n",
        "ax.set_xlabel(\"X\")\n",
        "ax.set_ylabel(\"Time\")\n",
        "ax.set_zlabel(\"Sz\")\n",
        "ax.set_title('Evolution of the magnetization')\n",
        "plt.show()\n",
        "plt.close()\n",
        "# fig.savefig(\"Evolution of the magnetization_tdvp\")\n",
        "\n",
        "# np.save('data_sz_coor.npy',Evolution_Sz)\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "Noj2bh3S-NEu"
      },
      "outputs": [],
      "source": [
        "Correlaion_SpSm.shape"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "eUbOwvxI-NEu"
      },
      "outputs": [],
      "source": [
        "plt.rcParams['xtick.direction'] = 'in'\n",
        "plt.rcParams['ytick.direction'] = 'in'\n",
        "plt.rcParams['font.family'] = 'Arial'\n",
        "plt.rcParams['figure.figsize'] = [6,4]\n",
        "plt.rcParams['axes.labelweight'] = 'bold'\n",
        "plt.rcParams[\"font.weight\"] = \"normal\"\n",
        "# plt.rcParams['xtick.fontfamily'] = 'bold'\n",
        "# plt.rcParams['ytick.fontfamily'] = 'bold'\n",
        "# plt.rcParams['font.sans-serif'] = 'bold'"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "t-Kl8MI2-NEu"
      },
      "outputs": [],
      "source": [
        "\n",
        "phase_factor = [np.exp(-Energy*1j*tau*i) for i in range(-time_slice,time_slice+1)]\n",
        "\n",
        "Correlation = np.transpose(np.array(Correlaion_SpSm))\n",
        "Correlation = np.concatenate((np.flip(Correlation[:, 1:].conjugate(), axis=1), Correlation), axis=1)\n",
        "Correlation = Correlation @ np.diag(phase_factor)\n",
        "Spec = np.absolute(np.fft.fft2(Correlation))\n",
        "Spec = Spec[:, 16:0:-1]\n",
        "Spec = Spec/np.size(Spec)\n",
        "# Spec = Spec[:, 1065:1090]\n",
        "NumOmega = np.shape(Spec)[1]\n",
        "K = np.outer(np.linspace(0, np.pi * 2, N), np.ones(NumOmega))\n",
        "W = np.outer(np.ones(N), np.linspace(0, NumOmega * np.pi / (tau * time_slice), NumOmega))\n",
        "fig, ax = plt.subplots()\n",
        "im = ax.imshow(abs(Spec).T, interpolation='Spline36', cmap=\"jet\",\n",
        "               origin='lower', extent=[0,10,0,6],\n",
        "               vmax=abs(Spec).max(), vmin=0)\n",
        "\n",
        "cbar = fig.colorbar(im, ax=ax,shrink = 0.74)\n",
        "cbar.ax.set_title(r'$S(k,\\omega)$')\n",
        "\n",
        "plt.xticks([0,5,10],['0',r\"$\\pi$\",r\"$2\\pi$\"],fontsize = 12)\n",
        "plt.xlabel(r'$\\mathbf{k \\ (1/L)}$',fontsize = 12)\n",
        "plt.yticks([0,1.5,3,4.5,6],[0,1,2,3,4],fontsize = 12)\n",
        "plt.ylabel(r'$\\mathbf{\\omega \\ (J)}$',fontsize = 12)\n",
        "ax.set_title(r'Dynamic structure factor(S = $\\frac{1}{2},D = 10$) by TDVP')\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "4pK4JLo7-QHS"
      },
      "source": [
        "For comparison, we also show the result calculated from TEBD for half integer chain here."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "nwSUVVJq-cea"
      },
      "source": [
        "![image.png](data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAfEAAAE6CAIAAACJf/xdAAAgAElEQVR4AeydCZQVxdXHa9jmAQMDDPueuKGDohEjRhQE0YkLoEHlE40kKlFjBAY1CAguoOAaIy7BjQQSIRJQEWJQMIgEE7egEhQBlUEFRaVBBBmc/rCv7z+XW1X9+r3pGWapOXXm3K6uulVdy69u367up3z351rAtYBrAdcC1aUFVHW5EHcdrgVcC7gWcC3gO6a7QeBawLWAa4Hq0wKO6dWnL92VuBZwLeBawDHdjQHXAq4FXAtUnxZwTK8+femuxLWAawHXAo7pbgy4FnAt4Fqg+rSAY3r16Ut3Ja4FXAu4FnBMd2PAtYBrAdcC1acFHNOrT1+6K3Et4FrAtYBjuhsDrgVcC7gWqD4t4JheffrSXUn1aIH777//qKOOqlOnzoQJE6rHFbmrqMgWcEyvyNZ2ZbkWSN0C8+bNe/rppwcNGuSYnrqxXAqtBRzTtSapARElJSWxX2V56Iy9kr7vV5V6XnzxxY7p5TEAqr1Ox/TSLu7Vq5dK/mVlZeXk5PzoRz/6/e9/v2fPntJE+0nq1KnTRRddFEvhTz311IUXXhiLKigpD51QTkJRUdGJJ56YnZ3dokWLHTt2iLMRD8tYzzvuuGPIkCEo66233jrvvPNatWpVt27d1q1bn3vuua+//jrOllHIgOkLFixIjl9Vu3btZs2anXDCCb/73e927dpVxspQ9m3btl1zzTX9+vVr3ry5Usq45Gzfvn348OFt2rTJzs7u1q3b448/jqLPP//82267DYdOKKcWcEwvbdhevXodddRRK4K/l1566emnn77ssstq1ap17rnn7nfj7vXXX1+7dm1pXcsg9Qr+yqDAkLU8dIpiRo0aVadOndmzZ7/00kviVPTDstRz9erVTZs23bhxIxX39ttv5+Tk9OnTZ/bs2f/85z9nzpx57LHHZmdnr1ixInp9QlJmwPSbb75ZKTV37twVK1YsW7Zszpw5w4YNq1OnzpFHHvnll1+GlBXx1Pvvv5+bm3viiSdecsklNqb369evSZMmDz744JIlSyjZn//8Z9K/cePGpk2b/u9//4tYXHkn+9e//jVo0KDWrVvXrVu3WbNmRx111MiRIzHTb7zxxkMPPfTbb79FNbZu3ZqVlXX33XcjJkbh4Ycfbtu27VdffVV2nY7ppW1onPB33XWXUmrmzJml6aq4ZLzMMl5TeegUVRo6dGjHjh1FZLqHZannmWeeefnll6PEX/7ylx07dty9ezdivvrqq/bt25922mmIKYuQAdPPOuusnJwcUIlKnz59ulLqN7/5TVkqQ3lLgj/f9z/77DMj0+lG4S9/+QvK6tevX9u2bXGne/nll/fv3x9n96Pw4IMPZmVlnXHGGY8//vgLL7wwc+bMYcOGHX300VSljz76qGHDhk888QSv4fPPP6+UKotJwbUJubi4+KCDDho/fryIz+DQMb200YwTfs+ePW3btu3Zs6fv+1dffXUikdi6dSvy3HrrrTk5OV999VWnTp3Gjx9/9dVXt2zZMpFInHLKKe+++y4l+/rrr0ePHn3ggQfWq1evUaNGJ5988htvvEGnOnXqdMMNN4wcOTIvLy8nJ+f888/fvn375MmT27Vr17hx47PPPnvLli1ICd9LSUnJ1KlTDzvssEQiccABB0yZMkVMY8ry2muv9enTp3Hjxjk5OX379n355Zd93+f+pRdeeMH3/U6dOo0YMaJPnz6NGjX61a9+JZw8jz32mFLq/fffJ53Goo06Udu9RXAlojhS+9BDDx122GH16tXr0KHDhAkTiouLKR7/9zYUvAp0yx/SqsZK6te+Z8+e++67r2vXrolEokOHDr/97W937tzJS+TN8tZbb2VlZfH5fNppp3Xs2PGbb75BFt/3n3jiienTp/MY3/fff/99VJ4LvXr1Ein5YQZM79ixIw1Ursf3/fbt27dt21ZEluXQxvRLLrkkJyeHd99f/vIXpdTy5cupuOXLl2dlZb311ltlKb3seTdu3FivXr3/+7//E6pglV977bXt2rXDISWbPHlynTp1vv76a5ErrsM77rgjNzc3Y78iquGYjqb4jnfGaXbhhRfWrVu3uLh49erVSqmHH34YeQ499NBf/OIXRMbc3NzTTz994cKFM2fOzMvL69GjByUbNGhQixYtHnnkkX/+85/Tpk1r1arVIYccQhTu1KkTsfu555675ZZblFKHHHJInz59FixYcNddd9WuXfuKK64gJRy1o0ePrl279jXXXLNo0aLJkyfXrl37pptuQpVI8DyvRYsW55577qJFi5555pkePXrk5uZu3bp11apVRwV/K1as8DyPal6nTp0RI0YsWrTopZde4gUJHPu+byzaqDOE6bw43/dvueWWrKysq6666h//+MeUKVMSicQvf/lLcTmvv/76aaed1rp16xUrVhQVFfm+H9Kqxkr6vi/qefHFF9epU2fs2LGLFi2aMmVKgwYNTjnlFKyOnTp14vUcPXp027Ztcdb3/fvvv18p9aMf/Wjq1Kn/+9//+ClR+V27dpFDT/xftWqVSEmHxcXFO3fuHDp06NixY3fu3Akj15gYkVu2bLHZ43369MnKyhKEKikpKQ79g2ZdsDG9R48exxxzDE//9ttvK6X+8Ic/UGRJSUn79u2vu+46nqbi5ZkzZyqlHnvsMWPR33zzTV5e3jXXXCPO/uxnPzvyyCMpsqSk5K677srOzv7tb38b3kEbNmxo1KjRlVdeKbRt3bpVKXXZZZch/pNPPsnKynrkkUcQk5ngmF7abjamX3PNNXtbf9OmTb7vH3fccSeeeCLl+fe//62UevHFF4mMnTt3Ru/eeOONSqktW7Z88803p5566qxZs1DMnXfeqZT6+OOPKVe7du1g1xxyyCGNGjXCfcAZZ5zRrVs3ygjUfvnll3Xr1h05ciQUFhYW9uvXD4ckrFixgt8nrl279pprrtmwYQOZq3zp6tSpU8eOHTHhURDp4SZ2SNGi6UKUiOK2bt3aoEEDPrIffvhhpdTbb78truiiiy7q1KkTRYa0akgl+bWvWrVKKTVx4kSUMmPGDKXUwoULKUbU88c//vGAAQOQmITrr78+kUiQ6d28efMhQ4bQzZBIlu7hhAkTuDlvQ49Qu2jRIqXUo48+KuJ93z/++OMbNmwo4l944QVeii7j5kxkDPG9HHTQQaeeeipP//HHHyulbrnlFkQOHDjwxz/+MQ73i/DUU08ppbp16/bMM8/wmzOqzIsvvshHAmrYqVOnSy+9lC7/9NNPz8vLw2hBGl0YOHBgXl4e2U/ibNu2bbFI0KlDDz307LPPFsnSPXRML20xASacIKZv3rx5723sww8/nJWVRcP98ssvP+CAAyhZp06dzjnnHGR55JFHlFIffvghYj7++OOlS5f+4Q9/OOmkk5RSH3zwATH9rLPOQppevXpxM4dTDJT8+9//joUEGXVh+/btLVq0aNq06WWXXfb000/zGzpxmZ06dSooKIAGFEQxnOkhRes6Q+x0XhzpXLhwIUzGzZs3K6X0J1G8NVBbvVVDKsmZTlY2x1ZxcTHd/ZBy0SzNmzcfMWIEyoXwxRdf/OUvf7n44ot/+MMfKqWysrJ+97vf4SwEXB0XYAEgWVmEW2+9VSll3HjTqlWrH/3oR0L5tm3bXgn9E24lnt1mpx900EG8c33fJ6bfeuutyD5ixIgWLVrgcL8Ie/bsOfvss2kZy87OPumkkx588EFc75QpU2DDoXp0yQ899NDSpUvbtWt3/PHH0/0iEhgFMvv4ksaT/eQnP2natCmPGTJkSKtWrXhMBrJjemmjCTDhxPnnn1+/fn2agdu3b8/Jybn55pu/+eabZs2a3XzzzZQsBIXPPvtsly5dlFKNGjU6/vjjCwoK4KEWuUQFOMWQkm4b33nnHVTPJrzzzjsXXHBBkyZNlFINGjQYNmwYmSSiFGhOeSEhRYfr5AuDKI506kaifufLW8P3fVurhlSSM33ixIlKKWGjtWrVigwxWm75srTX+TZ27FhbU1P866+/np+fn52djacgFJ+ZPz28LP3sOeecU7duXYAJCf7zn/8opa6//nrEkLBffC++748dO7ZevXqiMvvlcNWqVVOmTCkoKMjOzt57U3vyySeTA2348OFZWVlixV24cKFS6txzz61Tp84111yDe+vwmg8dOrRWrVp0i+/7/uLFi/kD5G7dujVp0oRrGDlyZFZWVkTlPCOXHdNLW0OAiU7QM9KTTz4Z6X7xi18cfvjh8+fPr1WrFixxgSpQbO3atYlE4rzzzsNOxPvuu68sTH/mmWeUUsuWLUN9NmzYsGTJEr4BA6d839+zZ8+//vWvyy+/XCk1efJkzjVKJmreuXPnCy64ABruvfde1DakaNF0IUpEcaTzz3/+szAZyU2Eaux90siZHtKqIZXk167b6bt3765duzZcvaKebdq0+fWvf436bNy4sU2bNvzJCp2aN2+eUurf//43Uvq+/80334iro8MoCzPXEy4fcMAB4kaeXrA65ZRT6tev/8knn4js5eF7ufTSS8Uz0scff5w/I/V9/4orroj3ga24rgwOP/30065du+5FNu2zHDZsmL7q3HTTTYlEgnavRi+iXbt2hx9+ONKfcMIJ3LXSunXr7t274+zeoXLdddcppbZv384j05Ud00tbTICJTkydOlUp9be//Q3pXnrpJaVUr169OOgFAsD02bNnK6X++9//IvsZZ5yhlFq3bp1uDIoKcIpB/2effVa3bt1Ro0ZBYWFhYV5enljbn3jiiebNm/OZ3KRJE6JSnz59hD+dG6RHHHEEPzt48GAwPaRooTNECS6E6r9ly5Z69erxV1FWrlzZt2/fpUuX4gJJ4K0R0qohlfR9H/XU/el/+tOflFILFiyg4kQ9e/TowTfh7dmzp0OHDkcddZSw9CdPnpxIJL744gtR+eiH+i0LYkKU0NZpelyPZLt37/7Vr36llHrwwQcRCaE8fC9kzPKnRwUFBXwvo+/7Z5555nHHHYdqVBLh0ksvhedqzJgxSimxVbx///7HHXfcs88+W6dOndGjR0ep9tdff7137g8cOJASf/7553Xr1oUH76233lJKifvRX/3qV9nZ2VGUh6RxTC9tHP7O0fLly5966qkrrriidu3a3G6l1ORLwcsUOp3B9Pfee69OnTqnnnrqokWL5s+ff/bZZ2dlZSml3nzzTT1XFKb7vn/NNdfQho3nn3/+1ltvrVu37h133FF6GYH08ccfN23a9Nhjj503b97ixYuHDRuWlZW1ZMkS3/fPPvvsNm3aLF68mNAj4DVmzJisrKxJkyYtWbKksLCwRYsWYHpI0UJniBJR3N47iXHjxmVnZ48bN+7555//4x//2LFjx737HfGgGNfFmR7eqiHtw+t50UUXUTM+99xzt99+e8OGDU866STbs+Lx48c3b96cb2555pln6tSpk5+f/8ADD/zzn/9cuHDhiBEj6tSpQzdDqHa6wq5du4YOHdq+fftGjRode+yx2AUIPWRP4JCEJUuWKKWuuuqqFStWLF++fO/Nyg033NC5c+fs7Oz7779fJC7L4cKFC5944olHH31UKXXOOec8EfzxpzX9+vVr2rTptGnTlixZQqDk73aUlJTk5eXdeOONZalDGfPS9gSu5IsvvujcuXP79u3JMKLVfeXKlTxNmzZtaO8KTe2pU6fys0Z5+/btSik8YJg0aVIikcBOtksuuWTv+4ziDax+/fodddRRRm3RIx3TS9uK77OuVatW69ate/fuPXPmTD6TKfWoUaNyc3P5TlWBKjCd9izn5+cnEom2bdueffbZS5cuzcrKuu+++zJmeklJyR133HHAAQdkZ2d36dLFNmn/85//nHLKKc2aNUskEt27d587dy5VfsmSJR07dqxXrx6tSaLmX3311aWXXtq0adOcnJzzzjvv6aef5ky3FS10higRxVGV7rvvPtqf3qpVqyFDhsCjVdo3+/pewlvVVknf93k99+zZM3HixB/+8Id169bt3Lnzddddx41uUc+VK1fqTpXXXntt8ODB7du3z87Obty4ce/evfn9HK98dPmrr7668cYbP/zww2+//Xb69Ol7N9VxYhImBg8eLBTecccdMOezs7Np6N5www06v0TGdA/5iwIokT9q3r59+1VXXdW6det69eodccQR/NsAvu+//PLLexcD2ybOdCuTWfo+ffocd9xx99xzz+LFixctWnT33XcfeOCBjRo1wq3hhg0b+P5LPOnFBqSbb765Vq1ac+bMQQWMC63v+wcffPDety4eeeSRBx54oGnTpldffXXbtm3nz59P2y4KCwuhwff9b7/9Njc3V0TyBBFlx/SIDVWarKSkJD8/X99wWprCSdWxBc444wx943x5X2jTpk25427BggVZWVl0k1feRZeH/qFDh+pbQsujoBCd06dPP+uss37wgx80aNAgOzv7oIMOuvLKK2kfGnKdcMIJ/H1g2vvIm33YsGGJRIL2MdsWWlrDDj/88Ozs7K5duz733HOff/45fbNob+m33XabMBYXL16slHrttddQjcwEx/Q02m3btm033HDD6aefXq9ePTzzTCO/S1qVW+DNN99s3Lix8R6inC5r9erViUSCPzG7+uqr9bcfy6n02NV++OGHe9+w42SMvYi4FM6ZM6d27dr4tk+42rgW2gsuuOAnP/lJeFlRzjqmR2ml79MUFxe3adOmZcuWM2bMSCObS1pdWmDy5Mm636OcLm7Hjh3du3fnb0WVU0EVpnbw4MG2ndoVVoeIBZWUlPTo0YPvdArJGMtCu3bt2rp16/L9bCElhp9yTA9vH3fWtcB+aIHdu3effvrpP//5z8Xt+X6oSk0t8q233po0aRKemZd3MyxZsgRfUChjWY7pZWxAl921QMwt8O233w4ePLh///5if2rMxTh11bQF9hvTzzrrLHy+o5q2rbss1wKZtMAll1zSq1cvvgknEy0uTyVrgc2bNw8YMCA3NzcvL2/48OH6gr1gwYKuXbs2aNCgS5cu8+fPF9V/6KGH9n7fWEQaDyMlMuYsSyR9L8kxvSxt6PJWyxb44IMP9r7xm0gkGib/aHNFtbzYGnVRvXv3HjJkyI4dO9atW5efn8/fs/N9f82aNYlEYt68ecXFxbNnz65fvz5/Qks/wFJ5mf7RRx81bdq0ffv2juk1aky7i3UtUD1aYOfOnV7on36P9d577ymlPvroI2qBWbNmiR94GTt2LP+6akFBAX4fY8eOHV27dh07dmzlZfpPf/rT8847j78WWD162l2FawHXAtW+BXbu3JmDt60sQuvWrQXWn3zyyWbNmqFx3nzzTaUU/0HBgQMH8reNCgsLsZF/6NCho0ePpu/zQEOIUNG+l4ceeqhly5afffaZjenffvvtunXrNmzYUJT8+/TTT0MXRXfStYBrAdcCqVtg69atRUVFZdzK4nned99pUWocCydocBcfTJ8xY8be3/ACiNeuXauU4l/r7du3L/925rhx4/r27ev7/owZM4499tji4uJKyvQPPvigUaNG9E6tjelFRUVa+7gI1wKuBVwLxNMCnKSAbHSBmH6jUlNYmKTUjcnw3ZcVlRJMnzt3bl5eHkohO51/1Kh///7iw3wDBw5855132rVrt379+r2fFK2MTC8pKenTp895551HF2ZjOv2kk1JXKDUyGa5WarQL1a4FuKGTsXyDUilDxsp5RjcCq0ELjFRKcZICstEFYvokpe60hEkmpq9Zs4b/zsasWbPat2/PCx0zZgy+9rX3owIFBQXjxo2bOHFigwYNcoO/hg0b7v2tmNzcXP7pQK4BcsX5Xu699176+iv91MvPf/7zTp06FRcXi1sharKAXxOUcqEat8BEpcoeJiuVMpS9lIluKFaLFhitW9BAYUSBAHWnUvdbwnc/TanZ6b7v9+zZc/Dgwdu2bVu/fn1+fj79VDoKpe9AzJ49m/a9JBIJ/Eg9pamMdjr/6mFw1d//E9fmmF4tJk+UpSgW1KYE+uQ4Vg7H9CgdWvnTxMb03yv1kCX83sL0TZs2DRo0KC8vr0WLFqNGjaLfUWrYsCG+Rfzss89269YtJycnPz8f3/EH9Csj09955x3+ay9nnHFGmzZtXnnlFezvodo7pjump0Nhx/TKT9LKU8PYmH6/Uo9Zwv0WpoPO5S1UnO9FXInNn+6Y7pjumF5jxkAF4z42pk9TaqYlTHNMF6x3TK8x8zkW38vEVP70uEqpYPq44sqjBWJj+qNKzbKE734CyuRPF6Arv8P9ZqfbLskx3TE9HTvdMb082FdddcbG9BlK/c0SZjimC7g7pjumO6bXmDFQwYtHbEyfpdTTljDLMd0xvaZO4Li8IuGPSeMqpYLp44orjxaIjelzlFpoCXMc0x3THdPTNMwFph3TywN/1VJnbEx/RqkXLOEZx3THdMf0sjE9xKUu6F+Ww2rJuJp2UbEx/R9KvWQJ/3BMd0x3TC8z041YLwvB9bw1DX/V8npjY/pipV62hMWO6Y7pjulxMF1gXYdyGWOqJeNq2kXFxvSlSr1mCUsd0x3THdNjYnoZqR2evabhr1peb2xMX67USktY7pjumO6Y7pheU8dABa8csTH9P0r9zxL+45jumF6D53O4aVypzlYwfVxx5dECsTH9daXWWMLrjumO6Y7pzk6vwWOgPNht0xkb099U6n1L+O5X6dy3ATjW3XukNWl6VypLPLwyNky4+CrUArExfbVSGy1htWM6B7rv+47pjumV0myvQuRyVbW1QGxMf7+2+ryOObxf29np+0LdMb3GMD3cLq5sZ22YcPFVqAViY3pRtvIS5lCU7ZjumF5Df6KvslE7vD5ViFyuqrYWiI3pmxqorxuaw6YGjumO6Y7pMfwwaTiRy37WhgkXX4VaIDamb2msdueaw5bGjumO6Y7pjuk1dAxU8HoQG9O9psrPMwevqWO6Y3oNnc9lt50rUkMF08cVVx4tEB/Tmyu/pTl4zR3THdMd052dXkPHQHmAO0RnfExvrfy25uC1dkx3TK+h8zkDK5t/Kj3d7GXJO7HGbEYKAWI1OBUf09srv6M5eO0d0x3THdMj2ekcykK28V0k44e2LLb4akA0dwnxMf2Hyj/IHLwfOqY7pjump2Y6x3Fcsg3ftnjHxKreAvEx/WDlH2oO3sGO6Y7pNZHpNm7a4uPiONdjK8sWX9WJ5uofH9MPU/7h5uAd5pjumO6Y7uz0mjgGKn6NiY/phyv/SHPwDndMd0yvifPZZgvb4rl9HZdsK8sWX/EMciXG2wLxMf0o5Xc3B+8ox3THdMd0Z6fXxDEQL6+jaIuP6d2V38McvO6O6Y7pNXE+22xhW3xctjnXYyvLFh+FGi5NZW6B+Jh+rPKPNwfvWMd0x3THdGen18QxUPH0j4/pxyu/lzl4xzumO6bXxPlss4Vt8dy+jku2lWWLr3gGuRLjbYH4mN5X+aeag9fXMd0x3TE9tZ0+Uam4UE56bOAOiY+XL05bxbdAfEzvp/zTzMHr55jumO6Y7pheE8dAVWb6T5Xf3xy8nzqmO6bX0PkcYhEbT+1fO73iAeRKjL0F4rPTz1T+2ebgnemY7pjumF4F7PTY+eIUVnwLxMf0gco/xxy8gY7pjumO6Y7pNXQMVDDW42P6IOX/nzl4gxzTHdNr7nw2+lhskc73UsEErH7Fxcf085R/oTl45zmmO6Y7pld2O7360a1mXlF8TD9f+UPNwTvfMd0x3THdMb3mjoGKXF3iY/ovlX+ZOXi/dEx3TK+589nmZjHG70ffS0Vyx5VVfi0QH9OHKf9Kc/CGOaY7pjumR7LTY3ztyLhmhESWH2Wc5opsgfiYfpnyh5uDd5ljumO6Y7pjes0dA1WT6Vcqf5Q5eFc6pjum19z5HGIa206V3QNj0xwSX5HccWWVXwvEZ6ePUP5vzcEb4ZjumO6YHtVOj8X9EsJu26nyo4zTXJEtEB/TRyl/jDl4oxzTHdMd079nOrfBbXjlaTKTo2gWaSqSO66s8muB+Jj+W+WPNwfvt47pjumO6d8xXQe0ACsO9ZTRY6BECLoGnqD8KOM0V2QLxMf0Mcq/yRy8MTWJ6Xv27Ln11lsPOOCARCJxxBFHzJgxY1+ef3fked53TaJGq5pLuooc5fuxLA5NHamTldUhY0x8p1IiGJPxQrmcMvF+bChXdFwtEB/Tb1T+FHPwbqxJTL/22mvr1q07efLk559/vrCwUCn15z//WWDdMb3GLGZpITUksUC5OBSw5nq4LJLRIU8QF1acnv3YAvExfZLy7zQHb5KZ6Zs3bx4wYEBubm5eXt7w4cOLi4sF+hYsWNC1a9cGDRp06dJl/vz5dPaLL74YMmRIXl5eo0aN+vTp88Ybb4hc+qHSo8opZvv27fXr17/22muhv1evXj169MAhCY7pNY/pRp5GtNMFwY2HXD/HNJd5Gi4jzX4kkSs6rhaIj+lTlH+POXhTzEzv3bv3kCFDduzYsW7duvz8/Ntuu42jb82aNYlEYt68ecXFxbNnz65fv/7GjRt93z/77LPPPPPML7/8cvfu3aNHj+7YsSPPZZQrjunFxcX//e9/N23ahHr069fvqKOOwiEJjuk1jOkcoEIGT4XAkxkhLiJ5eqEKhzyNkClNXFhxevZjC8TH9NuVP7U07Lpbebd/H4puNjD9vffeU0p99NFHRLlZs2YJOo8dO7Zfv36AYUFBwfjx433f3717986dO33f/+KLLy677LKjjz4aaWxCxTEdNSgpKfnkk09uueWWrKyshx9+GPEkJJk+MnCpjw7+j60xmNuPw72Ci9ZhSiDmPEUaXaBknN33KCUCPwu1uiqKQYLJSb88YpClgpvIFRd7C8TH9LuV/2BpmHD6dxznf57ncbI9+eSTzZo1Q8ybb76plPryyy8RM3DgwMLCQhwWFhYOGDAAh2PGjMnKymrcuPGyZcsQaRP2A9NnzpxJF3/aaVo/MXgAACAASURBVKdt375d1CzJdN4+vRzTq10LCJKm5O+EfR+ZcqYLlItDvlSAzhOD9uSHIHhITWJHjFNYwS0QH9PvVf7DpWHXA8q79/tQdNt37BJMnzFjRocOHcC6tWvXKqWKiooQ07dv3+uvvx6H48aN69u3Lw6//vrrXbt23XnnnQ0bNly3bh3ijcJ+YPp77723dOnSadOmtW7d+vDDD6c7C1QuyXRnp1fwcK/g4jjTOUaNCEbdQGEb0+9XigLInoFCXh8qiMpFNZxQRVsgPqbfr/zHzMG738D0uXPn5uXlgXJkp2/duhUx/fv3HzVqFA4LCwsHDhyIQwiHHnroXXfdhUOjsB+YjnosXrxYKTVz5kzEuL2M1c4et03+sjOdyAt2g+YQcApYx5KAWiGGLxKO6WifaibEx/Rpyp9hDt40A9PXrFmjlMLTxFmzZrVv355zb8yYMQUFBYgpKCgYN26c7/vHHXfcE088gfgf/vCH06dPx6FRqDimb968efr06Zs3b0Y9tm3bppS69dZbEeOYXjOYbgQ6IRg85fteQBaOYM50cFwIXKfN4tYrw3PdmXwZipKhJk6oii0QH9P/qPwnzMH7o4Hpvu/37Nlz8ODB27ZtW79+fX5+/oQJEzj3Vq9enUgkZs+eTfteEonEu+++6/v+iBEjDjvssA8++GDXrl3jx49v167dF198wTPqcsUxfd26dUqpSZMmoRJ//etflVILFy5EjGO6nenAWVWcS6LONoxymILpel5uU9+TdLbcr9Q0FgB36ATTdYX0FitfJJCr+jG9Og0k0ZUpD+Nj+kzlzzUHL3hcKPzpvu9v2rRp0KBBeXl5LVq0GDVq1J49e3zfb9iwIRwVzz77bLdu3XJycvLz8xcsWEBU3LVr16hRo9q0aZOXl3f66acT6DkwdbnimO77/s9//vPs7OzJkycvXrx4ypQpjRo1OvXUU0tKSni1kv509x4pH6CYhyTwU1VRBtN1jN6z77YTepjJrxFfEbgzudGF8M2BTjLFg84ZMJ1Xpjq1PIYTb9iaIMfH9FnKf9ocvFlmO51TrlzlCmX6rl27Jk6cePDBB2dnZ3fu3HncuHG7du0Sl+eYrpnqmIFcqNIzkHMZaCb4RkEw7HRKLCx0DndgPcTcDq8MrTpYD/Q1pgp1BB8/kKtQ/cte1fiYPkf5C83Bm1OTmC7wbTx0TN+X6Zh7ulD2IU4auOa4dIbooeIEl8mFAqzbEKzntRnpwlQHmkkDqqcrxDqRsjJQErtAtRJVLWMpXKcul1F5VckeH9PnKv8f5uDNdUzfF+2O6Yzp+tzjMbFMJK6Q5FjU2pSgOLzaA3TCpg7xeFB2nlcw/TGlHtO86lEUkk4APaQy5W2qo4kg2BozrXhoMwppqaq6ieNj+tPKX2wO3tOO6Y7p5i9QGueeiCz77BIK6bDsao0aUBaMdDheiKH3J13kKc1qOOI50wnoOtaJ1GT7ow5UQzrkiwR/6IqM4r6h/LCO6nHB2JhpRXJtNjkthVU0cXxMX6j8pebgLXRMd0w3MN028UR82aeWUIjDjDVDQ4gQwnTCOpiO3S/QhrywqeFAj8L0EIUpFxg9L2olhMxaTyjBYWbaeC6oChd4lmopx8f0fyj/JXPw/uGY7pgumR4+8cTZMs49oU0cpqVc5LUdAspgKDeNualOKYUeZNeZPi1wvHDfy7RgpyPM7YgKxU0Dd92UE9bFNYrDtHpBTyy0hR/q2atTTHxM/6fyXzUH75+O6Y7p+zA9fMrpZ8s45XSFekzKIvQsthiiqvB12JguPB7GvNz3AoMdAp0F06EQaMYKgQUGQE+Z13aNIr7srZdSQ3gCUZ+Uh+HaqvTZ+Jj+kvL/aw7eS47pVZ7pNEniGuspp5yeoCxF69pCYkRBISmxQRAg5gK5VoDREKZzCnP+Im90pgtzG/URleFMp8UAGfXKQAkWCVubpNV0QonIm9ahUBXlMC39IYmprJAEFX8qPqavUP5b5uCtcEyv2kwXk6Tsw1QojHiYcbkR9UdMxhlnkwVD4UIBTHWSIgsEpEEusUudxyMxsusC0hgzRsQ6Ljlic0VMVsGdm3FxyCiuC/H7XYiP6a8of7U5eK84pldhpouxS4dlGbhGhVEiMy40ivIoaYAzo6AzFIa2gCkOSUiZkSPYJqfUKRKQHhFprInxYhEZpd2ipKn4zs24RPFJZFxdWRTGmDc+pr+u/PfMwXvdMb2qMh3j1ShkNhCNqiJGVnyJVDEgjAtGAiJS4LIshzaO8/iy6Bd5cQlGgbcA5IjdF5Ks4nu2PErMTGe8ueJj+krlrzcHb6VjevVkemZbmEMmdspTmY1+XS1gJASeUpzih0bY6Va5YGVVPLRdKW8NIUdpQ56G5Lh6Vtdsi8mgRJsqxGegM/Ys8TF9lfI3mIO3yjG92jI9A6xjAmQgZDYBREGCQWkd6owrC6m5oW2To+u3aeDx0bXpKfVrT6vpRGLRKbH0rNAZfphuieHa6Gy6OssjfXxMf1f5H5uD965jepVkepRBnMFQjq7WmDKDacD1CLJEPNRxZrPKAVCdiSIGKcWTT9pvjrMiF6+JOIUsaSkUSvj+HP0UxfA6kByxGUUy3i8VbBykW5yoashhBuMz3iyxMX3r+6rkc3PY+r5jumN66f70kPkQ5VS6E4DrFExJeajDi2IE7DhMhSxS8kOkxE5zIVACZElZmXCFOAuFuoA0uiAS2yqTsklFAt47ZelZrieinFZxEXWmu1SkVYeIiWNj+udFqtgzh8+LHNOrOdPTGsrRp4ctZcTBTcm4EgEUOrThyRgv0KazT8SI9PyQUhLH+Rv/eEeUElAWY2UQSWnSUshrAllUXj9EShJQgSiCsfF572TcrVxJdDl6cdF1pjURolcgrZSxMf3TT9SuHebw6SeO6dWf6dFHc1ozxJg4rSEODTpTopCI0giWcQcFJyk3tAFE5BWqkFEAnX+fi39FgLLzq0D9qYiUCilBeH10v43wBaFK0MM3tqNKNoHXn2R0UGbdiuzpChGLKye1EUvPIFlsTP/406yvdtUyho8/zVJK6b9ztC/nyvGoQn8TI8p1VIVv7aY7lKMz3ba9N3qJ0Qc618lpYiOOiOfY4jJ4DYxympOMNMhIynEIehLEZyhFIYTp/BJI5jrTUggQ8/pEuRxKg1xCEA1oO+QXwvsos57lGqLLEcuKrhApI2oup2SxMb3oi7rennrGUPRFXcf0fVDvmK4UJoBNCH8TPeJ8gHIOERtobI89ObZA6hD8ge8Cf5y/ZOlTyseSNNexDqOY8vKrsDFdrBBYJ7DS4HL0+mBVwCUIQVw+VNmEkKbm14JuSrdbkRECqcWhTYhYUGb2R3Tl5ZEyNqZ/4CW+8OsbwwdewjG9JjA9uqlum2kUb5ztIkuUycCzQKdOGRuPRLzAmRF/3MSGywJ6dIbiI4ugORBMXnWsCpQXVwGB66Qq2ZguFNrsdA5xePZ5pN4OuMBwQW95XAXvqXR7lueFwnCDIK6ByovmcpRLKL80sTF9vdfwMz/HGNZ7DR3THdP5IOYTQJf5zCRZTxNlWiIXFAqshDMIZ3WKhQA9HOtCZ7pMFx/YwuXAvx9dIWpCAl0j2E1XIa4FZ40NIhTaDlFnEtA16KyUPctTchmqQsYM0vPRGCIjfVpCiMIKOBUb09d4jT/xc41hjdfYMb1qMT2tEcwTRxyyPIuQxczEoUiGQ2OJOMs/nShoYoMOZ5wRXimBzs1b0qCXBYama1aLq4C/KAOFvFbITuDmTOeXE451UgJ/EdfPZXEJxi5O2a3GLoYqEngaIRv165EiV/RDXVWFxcTG9Le95h/6LY3hba+5Y7pjOh/TIXNDTEtxGJLReArZBUc4YiADSeECuMYFGwShCqXApo7iKuG/dwFvCb8WqMVKE2WRENhFJXFFtstBAhKQMURADYXArwLdZOzEkEhkNAohGfloDJFDNISfCtFZ3qdiY/pKr9V6v40xrPRaOaZXIaaHD9bwsxHHa4gS4+QUkSHZcYpn4fgQZOGEDQETiCmghsMQCHK1vDg9L3emwzqm7Hq1eQylIYWoCXz0FIPiUB9eGXGBUCLyQgkXoDBc4BUmmfcL7y90YojA09vkkOxlH6ghyulUxCJiTxYb01/32r7ntzeG1722juk1hOkpnaEYwbYpYZufxnihRE/DqQE3BYdLOIYE6TjIdBkcFKdsRSAZMnIBZ0X2kMojC9cDGWeFQhwiAQl6RpFAP4SqEIHXX7/z0HswZRfrWRAj8vJDjMMQgafPQA7RXK6nYmP6K16H1X4nY3jF6+CYXlWYnsHYFVkijleRC4eYkGUUykJznVZljxGYEwoBUBLEWWyhEUpwqKcPV4iMJOjZyx4jiuCHAutRyJ7ZYMCgEkIZh6jQFnIYsaB4k8XG9Je9H7ztH2AML3s/cEyvEkwPGZ3RT0UcoDaFmc1e5BIoT8s2LzvInAa9BTjNuZyS7GKTD7o4umAbY2Ucoja1xviIZcWYLDamL/cOWOkfbAzLvQMc0ys5043DMbPIiKPTpjz6jMWPOOscN6JcuI85X3QSVcIYYXrjsBJWVa8Sb20u62TXzXbev2kND9sYK+MQtakNiY9YYizJYmP6i97Br/uHGsOL3sGO6ZWZ6SFjMYNTEcelTXOUScsnuZCNjOAQEbJOn3KKAYLx8BMFiSrRIc7yjOEysoQrpGRcFTKWt2CsGEUaO050Lj+MMk5sY6yMQ9SmNjw+YqFlTxYb05d4h/7H72oMS7xDHdMrLdPDB2IGZ6MPSpvy8OnKJzbJRhyEWOURn3xy6ukgjo4/oYcOkZ1jjtcZCSg9NrEYhRCdXL/QKSqGs+kKGejhtRKyrTf1fg93zthGV9nHp01zyvjoRZclZWxMf947fIV/pDE87x3umF45mZ5yFGaQIPpwtCkPYbqY2Pr8F4AQh1FoJQglDsM1UHEijdCAQ0pGWcSFcD38azCzlBIBiMeqE0Uh6iAEXnNeBx7PZZFdHPKUNplKsf0XzaK7ZUKGim10lX182jRHiY9eesYpY2P6Qu/opf6PjWGhd7RjumO6cYzapoFxroYD3cYFG030eIEk26GeUd+XgjRQAvjOUIoiKY14/YcoBp5SSsoraI5DOgudlFenIe5OuE6eF6tCyOXguvARAtIW8p9nCZdtPSiuRQwD41CxjSvjILRF2pSUJd5WVozxsTF9vnfMEv84Y5jvHVP1mP558LcPieM7qBzfZSzL0AzJm+7oNKoyTlQ+mfk8FywIB4d+NoRHtlO6ElEHAVBua9uYjqsD1sHNGUnzfI5SIhDZuU4sEkaFXCeWGVwmrst2OUgAPcibUuB5o8iiDrzHcWk294txUMUyMo2a04pMtxrppo+N6fO8Hov8nsYwz+tRNZi+ZcuWKVOmnHjiiYlEolbwl52dfcIJJ0yZMuWzzz6LD+m+Y7oq/Sk72+dMdazzmYwZzmd+FFLwNCkxZEvAlehWLVWJ0pAG0BMCjOJwBCP7LI3mgPus5Kd6SaeukHYBoVZigeGmOtWZNylkccm2lkkZL/SkPEQFcF0pPTBGwqZLRtuwNCpPKzKDmqSVJTamz/F6/t3vbQxzvJ6Vnelbt2696qqr6tevTyjPYn8Uk0gkhg8f/uWXX8ZC9krA9LRGYVqJ0xp/SKwXIZhuBLrxoWIIJlISJ2ICUQTnDmRKQwqBcgiC6frVkaUP/nKmP6XUUwzxVYvpooVFSxoP0aRYyAXWxVDRxxKGWVqCrieumLSqkW7i2Jj+V6/XM35fY/ir16uyM71Fixa1atUikjdt2vSkk046N/g76aSTmjZtSvG1atVq2bKlY3qqn7NIdwgivT5hxFwF+DC3Mdvh6+BQEPiI95AXZDTVKUGlYjqqhKUFAtWTEvBWhcyvN96WFNp4QXp90PUYDML9oo8iDLB0BV1VXDHp1iSt9LEx/XGvz1P+KcbwuNensjM9KyurZcuWV1999RtvvKFT+4033rj66quJ+/rZDGKcnb6v7wVDVsyZsjBdkCL6ITBHgi2jQA/YB4ESUHahk9zfABZxCpCiwxA7HV4XElLa6Vwh+cGN9eG3DrgKEsTF2tpEqLUlSxkviuNLps50MUjEEMLQykAQqmI8zKAy0bPExvQZ3sl/839qDDO8kys702+//favv/46nM5ff/31HXfcEZ4m4tlqzfTon/EyDlMxc/iMFdTjjhfdTk8JDp5AwEg/5IkhC/RwDuIUEnOdFElpuJsYwAJJue+Fu1+AdTwjxTpBeaEKApqIStfrA6ZzhiKXfkW4NFSS6xQyT5xSRlloIrQtLgeDgY8QMXiMAyxipFAV72HEOmSWLDamT/cK/uqfaQzTvYLKzvSILI4rmWO6xU6nQSzmDyYtpjF/9R+zXYAgJTiikIiDyahQFKof8lzQRpEcWEAVFio6SykpIzYvCoHOZqaT5+VM1y+E14dfVCzNKBSK0tHFaCU+EjA8xO/VZQZE5BKDMN5DlFIeQmxMn+b1n+n/zBimef0d0/dZDhzTQ5kuthzwSYvJjOmNCS9AIDDBD8HWdAWuBLIoVxwimS4gJS6BC3SWcqGegubwuvC9jMLKhk4Up9cEMUhjFJCMC6hbugJXImRROi4BnY5hwMcGx27ZWcm1xS6XvXohGsaVnbYEqAe8s6b75xrDA95ZZS9lHyameaBSpv9F6N/FF1987bXXLliwIKWeiAkqhOk0EEP6PvaRCoUhhUY/BW183mIyY3pjwgsvgcAEDtNFj0gPPVwQDBKHPCWXebKUVyGqIQ6hNi2dyEUCz6vLIjEdimqke2jUqd8uoHHQ6RgGfGxgwEQfYyEpoS12IWWhIQlSnpqoVGxMn+qd84h/vjFM9c6p7EzPysqiPYvh/3/605/u2bMnIrhDkpU/0/lAtI0DniZ22VZo9HhUic9bTGZMb0x4wXTjSzHpQseY3kYiHUYgo54Fp0IEPVeUmBCFOKXrwSkh6CkRY2ycdCOhjQu8Drx/0ekYBnxsYMBEH2MhKaEtdsFWKC/IliZlfJxMv8c7b5p/oTHc451XBZjOtqRbxVq1at11110hsI54yjE9le+Fu1/4vMVkxvTmc56zQGd6CG74k8aQZPwUZ5BRFpWxHRrzllOkrQ4iPmXpvB3CZTiLQpKJ4kRleP+i0zEM+EZGMDEl+FImgKryEGyl87JsacLjSUNsdvpd3vkP+EON4S7v/MrO9BtC/8aPH3/OOefUrl27Vq1axxxzTERwhyQrZ6bzwRGyC0Uki/cwfPBFOYv6gOl8JmN68zkvcMBhYWQKiKMLxvTGSF6KLosqiUORnvSLyBDzn2vTcxm18Sy6rCvhMcbLN0bq7YkYY3peiqgV7190On/nCMMDAybK6ApPA1XlIdiKFmXZkoXEk4bYmH67d+FU/2JjuN27sLIzPYS/ODVmzJisrKzGjRsjJmOhYpluw7oYQ/Eehoy8KKdQGcxY/gsYfG7zOc9xwDGh780AYsIFI4DCI0W5vEpCRkqbQkrAc/GLJZmfpfTh2sJXCFQpXJWtCHyXJrxV+aNdUiXK5RfFL5n3O1/gMUgwbKKMsZA00FMegrFcvSBjsvBIUhIb0yd7Q3/nDzOGyd7Q6sD0f/zjH7m5uc2aNQtHeUlJyR/+8IfDDz+8YcOGP/jBD4YPH+55nshSA5huW0jCByXOYohjumIO84nNJ7zwp3NMCAalJI5IILJHOUTpRjxRJOemsURKQInFlYrDsijkNUS1o1yjSCMuIeWhyI6ixaojrpT3PoYEBgmGDQZSZgL0lIdgrJJekDFZeCQpiY3pk7xL7vSvMIZJ3iWVnek7d+4U2DUepnwvyff9KVOm1K5de/To0c8999wDDzzQvHnzvn37lpSUcIXlyXR9cNjwakwZY2T4+As5izpgrmL28iktZruRTRlb6CFI4jASyfgpQEqHskAwlMxJfn8RliwARxfLL59kiodCqoBNoW2RQHZKwK8CqkgIOSVSRjnk2rAdk1dSbzrxihYGBhzrGDwhAyz8FDSUk2As3ViWMWVIJCmJjekTvMtv9UcYwwTv8srO9NatW999991bt27l5OWy53l3331327ZteaQuf/vtt02aNLniiitw6q9//atS6pVXXkGM75frdxnTGhzGxHFFhgy+kFMonYCOSctxFkJz8WhUUMMIGvFOpjFNupFULscTrzMYSsn0F0SNTOctwGXcoGABi6gQVUJ9uIZ0L9mWPmXz8j6iFuPtRnXTN93zFsAggRFAoyhkmIWcwggsD8FWrrEsW2JbPCmJjeljvStv9kcZw1jvSiPTN2/ePGDAgNzc3Ly8vOHDhxcXF3Pu+b6/YMGCrl27NmjQoEuXLvPnz6ezO3fuHD58eLt27Ro3bvzjH/94yZIlIpd+mHp/Ou1lrF+//qBBg6ZOnbps2bJ3gr9ly5bde++9gwYNatCgAW1z1LXzmC+//PLKK6986aWXEPnf//5XKTVr1izEOKbbN73QoMTMpLmK2QsGkYDZzgUOBQ4LnThgjS7oiTOIAdaperzyFEMAJc1UB3xtkSJDXvfX36TFh1yiKBSI5PWJ7hMPbxO9VRGjZ+Q9xXuQ9yxk3pLcbKfRgsFDY8mGv5B4I17jirSVa9RvS2yMh4bYmD7aG3GD/1tj2HvKyPTevXsPGTJkx44d69aty8/Pv+222zj31qxZk0gk5s2bV1xcPHv27Pr162/cuNH3/eHDh3fv3n3Dhg179uyZNm1agwYNPvzwQ55Rl1Mz/YYbbmjQoEHILvWsrKy9CW688UZde3jM+PHjlVJvv/02T5b0vYxUanQyjLWQzuY5MfYr3wKIPg7RwNOUh2yrpB6P0rl5rtMcExsCRwBkzggdUoBLuKCjR8TYslMyMF34T4RZDZuavqDLsQ6mE4LRGsaljjOdKmZTSPUBGaklqemozlT/8KsTTcEPbRlFPM9CMu81dCUX0OkQcBWicTImO8ZhOQn6yKcYY3G2xMZ4aIiN6dd6o673xyBct+vaa71RFIYXGez09957Tyn10UcfEetmzZrVsWNHzr2xY8f269cPMQUFBePHj/d9f9iwYQsXLkR806ZN586di0OjkJrpvu8XFRWNGDGiSZMm+u703NzcESNGbNiwwag9JHL58uXZ2dkDBw4UaZJMV+yvV01lOsaiTnPMWMxhPsO5zHEAWaeGwEr5HXJIgem4CjiIaL2ZE3wPHQgmrJP7BXm5QQp+kZBykeAK+Toh6hMF6OXRYqKb0H1c4H3NZVwChgq1SVlsdgzI8hCMOE7XFDMqQW1jY/oo79ox/vUIPSecyHj1nSh2fzz55JN8C8mbb76plOK/OTFw4MDCwkLAsLCwcMCAATgkYfHixbVq1Vq/fr2IF4eRmE55du3a9fzzz0+cOPHXwd/EiROfe+65iE9QRakvvvhibm7uYYcd9vnnn4tTSaZXJzvdOM5skRh/E5WCScXNT5qimLGYxnyeG2UBCDosDxJBJ2cxrGw87aRKUv3F5eAGQjCdNAuXOjWIoDn/4BcVgUVC14llxlYZXh+qg+3ScO0xCsaOM3Yxj8TAQNvyhjKSncaebWQa4/lwLaNs1F9JmT7CG/1b/waEUbvGjfBGU7i8aKTO9BkzZnTo0AGsW7t2rVKqqKgIMX379r3++utxOG7cuL59++LQ9/0VK1Y0a9bspptu4pFGOQ2mG/NnEPn4449nZ2cfccQRn3zyiZ49yfTRFtucd3yI54Qng2wcczgrBGPidCOFTtuhrjYlzTFjMY2NM1+PjJE1QpXAnO2QcgGjYlOHeBRJiaEKebEeELBgigqBznKdNoV8nUDb8oy6Iwi1MgpUUHn81/vUGIOBgcvhbRVis/PRaBuxIp5nyVgWOnFoU4gEKQVoiM1Ov9IbO8q/2Riu9MbqTJ87d25eXh5wR3Y633jSv3//UaNGIUFhYSH3YTz00EMNGza88847kSBEqGim33bbbVlZWb169eLXw+uXDtNT9iVPgH4VAk/DZZGMHwK4JPBTXObabDJPTzKU22xzmqI0YzGZM2OHEUYVEEm15VgH2UEiOkspqUq4RmSkphB2KGEdkWgurhCmemY1qYAmEkXg2tMS6JLFclVOZNdHshjP4ZMlxD6zabbNKT0eGmJj+uXehBH+rcZwuTdBZ/qaNWuUUps2bSLWzZo1q3379px7Y8aMKSgoQExBQcG4ceN839+zZ8+wYcOaN2/+3HPP4Wy4UKFMf/DBB5VS55577jfffGOrVoUz3TaYMA6EAOZyQaSx6QwZbSk9LWBTSj/AQqX0IBhRTod6uSIG5XKYAjoRFyowHSsBfOtAOfzyRltbJ6PAn7EmqLy4KP0QKctV0MtdaHr8wJuarouvhXgaQTaEvpk9s8FMM4LPEcj6ZKEYfXZQTLrpdT3QEBvTL/EmXeHfaQyXeJN0pvu+37Nnz8GDB2/btm39+vX5+fkTJkzgDFy9enUikZg9ezbte0kkEu+++67v+7/5zW86dOjwwQcf8MThcsUx/ZNPPqlfv36nTp2WLVu2gv19+umnvIqVm+kYl7qAcRM+QPmA41lIobDNiVCgObc0AQtM7EVKpQxIrAtQGFHQNUSPQRFgK5DKBZzlAhLAoudkp7bCf54GGbk2yDjLBZxFhaNfo54SSiIKugbEpOzoRWxdR3F0OXSBaDSY7eHeGD5ubTIfz9xGSTlZQqaM0MkPbdUwxlPG2Jg+1Js8zP+dMQz1JhuZvmnTpkGDBuXl5bVo0WLUqFH0FduGDRvOnDmTAPjss89269YtJycnPz+fvl7+2Wef1apVq27dug3ZH9JzbHK54pj+yCOPiEfDdPjYY4/xClVZpmfwUzIYoAA6zSt4D2w01+f2C0pFDDoOoI0LAAEXeAKjrCtHjEjP1QKd4QKnLWRObZuMxFwILwtneT3FJeDSjIJILA65WsgiDR3qyiN2NCWj7NBMZfHH1DTGsCsfH+aLDAAAIABJREFUZKcxiSFqBKWIROJwoIuZglxCGx3irC4Y09siKXtsTL/Qu/1if6oxXOjdbmQ6p1y5yhXH9IiXUZWZzgerbWwhHmNUBzr8BtzNQrNxYdIYx8RephSFl5UKCUi2bF/6c2Rg5kcUeN7osq4cUANxgFTQh+M4XpmXBT87qqTXNvqV8pS6nvAYnhd9/UKyr9PqcWQnnfDP0IWT44uTPWOsiyGtm+eIQUouYGpwgScQMk+WUqa8sTH9fO+uof4DxnC+d1dVYvqWLVvuuuuus84668QTTxwwYMBNN91ELztF5HWUZNWF6Sm9kDTIONC5eS5oTvN/URLHYkq/qlTEAOgD8WLCY9rbiMNZI2SogiAS2A55WYApiBNCcIHjkMMQJfqGFl4fW51FPC6ZCyINDrl+ISMN+ppzHN0XsbtfZcs89ThVD54ZbM8H2bnBLnzr4dzkwAW7bQJPDNmoH2d1wZjeFknZY2P6ed49F/rTjOE8754qw/SJEyc2bNhQ/NpRo0aNXn755SiwjpimijM9oqlOI4wDnfZTizcen0p6RTHDlyVnKWb1SqUorFLKFpBmJaM/ASIc7hwxQubwiiKL7LZDYhw2s4PF8KsghgiONSBEEMtDiCoq3VY3ER/lqnkakd12iCzoGuopvcdX2nt8VXJUGHucitDJTk72zLAO5to4LuKRHoIRxzirC8b0tkjKHhvTz/Gmnu8/YgzneFOrDNP1l0gp5vjjj4/I6yjJKjHTxaAEkUU8H3zGEUYJkB0+dCPQjTQXEF+jVJQA4lN2YgRMP27HYc4L7gA3ugAACd+OnjI8RgeNcVcMB7owdcUhsZ5jHUy33QyhzcOrajwbsR1Ew4oSSQm6hnoq3k5Hd1ODw2DnWDc6YUKGNB/YmBQw9hFDAp8mJKfULLIY09siKW9sTD/Le+Bcf7oxnOU9UGWYXr9+/XPPPffee+995plnnn/++T/96U/HHXdcVlZW/fr1o8A6YpoqxXSx/Quj1jb4EE8pscuFXmGnXXfkCoDBSLfeLyfta5hmgPh6pRCKlDIGJFifpD/xXYc7kGQEFiKRLKWALLqg56U0gjLhTEdD6ZTkywOYZdNGt0Rgq61u+lVw34ieS8QYsyMSiblhzlFOnc5709jdRWxIoMfXJI16WO4vB975F5JPaMjPLpww+NUkAWKgE0OaBEwBPsKNWBcZbb5KPRmPQTVSCpQrNqb396b9zJ9pDP29aVWG6WLToe/7v//977Oyslq2bBmR11GSVR2mA+gQ+IDmI0/IfLjrLhcOF+5pETSniY0pvUmplAGJKS8BAjfpuLUnoIAvUYSUWYAtErhOWKNADFGGEEPPSDO2rDnW8bhVaIPJT6sCFlFUjNc25EJ4MnSciAw/RInCMOcoRyem7O5NbIHHMhBOdrS5sNZphNO4FYOZH/Lxj0kBgZ/lLkpoMHIZZ42CMYsxkrLHxvQCb/qZ/l+NocCbXmWYDiIPGzbsmGOOadmyZa1atbKysoYNG4ZTZRfKh+nGAYFIfRDgFBfEoMRg5YJIQ4f61i58D1ZY6DrQXw0co2RWr0naXzSxMau3KIXgKWUMSEC5SAOMON1s54/XwJroQji86KxRG7luhKluozDtAtLXCdJP8VyVcYWgvS5QBRyL6mV8RUKP7RBrqtEw5z2OrjR2tMcGwxa20qPHqdPR49TRerMLrJPBLoa3PrApAZ8RkPW8fHKRHHEmIqOe3hZDWWJj+snejJ/6fzOGk70ZVY/p3bt3J5rXq1dv5MiRIS+FZoD4ysp0MSIxUo2CSKwPdPGFKXLshgCd05ygrE/sHUqlDEQBypsS7rhJB25sgo1TIRxEFq6TIokvsBnpoShhnWROYfKWoCBoANbJ+o6oDUDXa0W3EShICLgcIXA9ITJx3HYrxnucejBlR+9gSzuGCnocN2r8Li0c6/S6KR/tKQc5T0yyyAI0Q9CJjFNGQU9vi6HssTG9j/f4Kf5TxtDHe7zqMf3oo4/G89IGDRr8/e9/z4DdtixpMt3mhhNdaxwQPDJlejEc9fGaMoaeOBlpziH1QrCzBeY57ruLkjYXTVExt3crxUOJUiLgLOHACHdMddyhY84DOuGgF0SzcRDJdMzBA8Pta+x14fsOYVnDow21KNdmqkOhaHliurFWpFygnA55uSTrGhDDW5Kb5KuSjzqoF0K6G10purhk3zGwm63x4d1NvUwGu9HxxT8rj2GccsDrCcQk4hOQ5JTTUGQR6W2HyBUb03t5f+3rP2MMvbzvfr5NfGvXhrvyiM/knaONGze+9dZbTz/99GmnnZaVlXXkkUfGWLN0mG4cB8Z+RafaBJ7LmEYMx8lK6UMWMRj3EPC9QP11f44VMIVMNjLPcd+t09w4vX2lQgKBgDIKvsOOw0M2+NzpVj0i4nXGcRTirI45OHyWBQ/uyLjGvkZyfKO5wHQoh2YwndueMNWhh1Tx2yMQGeRFJaEcxRmBjvRcgDYBcZ3jOsqpj6i/APGU/ct7GXAXZIc3hpZwWq2p6cSCSn4YfFQHgxnDGwKmgC7oM8g40VLORJ6LJw6RkSU2pvf05vT2/24MPb05VYbpl1566auvvsrxTR9Z3H/7XqirQvoSp9CpNgEpjd9r1ocjxegDF4Mb4x4/y8DfDiVnC/AEBwLtbxFA35T0kHLbPPokTzn/wXfbnIfb3Wa/c36RDALaBGThvENezmK+GZFkvtcFRjrcJlQiyMtNdbFIQJtof7pJQsVQVdu18HgkJgFK4FcBx7k/jVBuXLOJziGdGHIKawC6GD4ZKoueo8LDzrFuNNgF3LGTnUY7Br8+L4xbX1I+KbXNVsTzaRsiI31sTO/hzevpLzKGHt68KsN0+vm6xo0bH3HEEcccc0y7du3Iq57y16X5MpBSjmyno5+iuF94YptMY8J4NiLTMaZpiHOT3Ph1QG4hwkKnOb8+uWmBgG6jech8TvcU7txhvPNpb7PcwSzBMo45o4z00EACN9WBFYI4NZcAOmc6OM4FMJ1sT6MqKOF7RlExVBV1M14R8iI9NOAuB5400b9bku7vHUn/SUnozVYGnav3Lz1NxfYYTnajB0xsHMKH0vAVME52gXXbDDJOt5CZKNKHcJyfQq7YmH6MN/84f4kxHOPNrzJMb926NbnR8SopHY4ePTolqaMn2H9Mp7UB3Q/BNhyFqc6BzmkuHu7p9/uguW6eC6CnY5vX8VV4sPpncNvObTq4ZXSfO8gFlpFgAx9nH/JCQF7iMmGdnnOK/wRr4JhznMtIJrLjEBpQMVQGAi4tykVRYuQloAua47EnVuvdyUcgFl6H92YdS659ehk9y70xfOXmWKcFjN8wcecVf3CNb8XAbKe5AKyHzyBMNAjGe2Wc5QIHt03m6WNj+tHewh/7S43haG9hlWG653ljx45t1aoVHpA2btz42muvpY9GRqd2eMr9ynTe/ZDDRyQc6zSO4WnBC0TisR55gXGzHwJ03InrHtXQCYzJn/AVhRxf6QFnKf0+k5/0c8vOaNPBFQN+gX1As4CgSICMEHhGoBlctglIKQRbeh5PWaieVDoqA0FUW1wUFgMkQ0YBdP5oxEs+wySaa31K/YJu0nswJ9m/CbZ+G/qRa84A69xmx70OHGKAO/e2l5HpmHopBRvHeTxXEhvTD/eeP9JfYQyHe89XGaYTi0tKSlavXv3SSy+tWrVq9+7d4YDO4GwlY3pKoHOmw+UCoMNjzlFOhiEgIsxzuFZhwfH78dBbckFzUKCJr0ICJSN2GPhuIzvZdALrgFp0AfjTd9QAnYLUZDxyLkeUdT0oQlQ4pFYipfGQsutA12+8tA6tkyQ19UtIxzVhSzXvvhQ2uxHrxjWb74fBY2fcOcGLhV3/YnsM2elRZhDHbloyZ7dN5gpjY/qh3pKu/n+M4VBvSRVjegaYTitLNKbzforLny500mGUEUljV38jFEAnq1ygPHy3IoBO/pZUN+ZknXHLjhOhua8QWvuKB8QTO8B3yQUiO3lj6E0W8rCvCXbg8f0wRsaVPRLwBVwI0FFozlHO9ZS9VkYNWA/wOJS85wA6LdL70pyjnPoCXcP7qzXryubJpRq9BoM9krXOV2v0rK1z6Uqp9bCmGrfHkAdGfFQgfB4Zp17ESBvKEc/1xMb0g70XD/VfN4aDvRcd0/dhfkZMT4l13q/R5fCBiLM608nfYnwEij0VYsLDPKfHZdzfEgfQwYX2vtIDznJMECBK6RCCdTJIwTIj6dKKhCouQAPIArjAeISXHOsod5dzJwnXDBlFZCxAFYz0VEAnmuckAd2cLbp6T7VnZwn6YjFOD+uwGEIejItlm1oGDpkXgk+H4nM6M5TKgOnGDTAR5ynYbRO4ntiYfoC3/GB/pTEc4C13TC9vpvNOTUsGteFgMW7Jwov+/MOKBHTY5jQTsKGN5gnf/AAjjrtZ038oSvfgMNJp2hOyCRCdfWULlACWIN3aJ4Q3lh6v0a262OAcl7UusEhw5MsG9p9EZzotA3xB5WrLqeZ8t+KmYGeLZqHXCfwnTZLWd+vkimvro87JBNSnAuvU+wYHGu9EkskDA1OdPzWFH0ZsieF+NkF2PD7FFwXgUueTiGTjs9O0JiZPbEM5xfOUE5WKjek/8F4+wH/bGH7gveyYni7TRT/RYUjXGtOnjORjEaOQBP0Ud7zgw4qLgu/ekW3Iab4++GwLHpcRzfluNjjQ42M6B/qBvrIFQgkZg2Sz04O4UlcMrHV+k479EhyO3D8OTJMQbgVTGti5tP7Rf5AduzII69w85zI30uHswh0SNEMtr7OxkrYL4fG85txI14CeYEAnmlP723rnwOR6rK++WIDTs9M503Ws43MxWLxFR6MXyAlDjnV41Y3ul5CplHJKGhOETHx9/0xsTO/gvdLJX20MHbxXHNMrOdPFKMQhN97BdP7lFv0t/3Caw0K3MX1fD2ypVyRphcGfrtvpOtO7+MoYODjIZieslxZHIIA1JxzrAuscdlwGMXUHNyUDeanR6L8wFeEBAMf57nVsLsKyyh9HQy1WCyoXFigqhqrCzOcXImQb0L1g4znrQQBd0NzYKV2SyzAWXay7WHqx+pYuwMmBUdp3wkg3Mh17WKmLQ8iOd5TI9/WUUvSpXuF+obtbTBwhwEIyIjtl5P5helvv9fb+e8bQ1nvdMb0yM52PP9gdPJJkYjptd6FHozDSBUfEy98wzznQBdOFtc64IOdqsKctEWyc4FiH7wVmIMDR1Vci0CkiO8ARCevciEtJdoAS6OTuEZAR7imYioR10ATuF8J6ONNpnbDppDrDsYOKoaoC3+KQ6syXIjwjoc5lHacDnWzzLlp3dGVLr1huacXF7hfq9xRAh8tF0JzGGHnV4VunB/W0p1YnOzqCVlaY6uR+0b/5xWeNmEqE9ZT4NibYP0xv5a1s4683hlbeSsf0ysB0jAyMGxpnQDb5B+ElFAMUn8yF44Xu+smyE35VOFtsQMfUopkm/mNaWqwwo7XOsQ58GAkCiHCs49a+dBUhKAhrnbAOaxqYE/iDiwPchAAPid5oOtZhp4c8I6U3APjKyr1euk6d6Xiyql8FXGpYhES18cQ7FdBpKRXrKx3SKeo18qfDmc6BnsKNjmEjhhMdYshxgZgOa92IdepKWlkXKmVzv2C+8HkEsgum65MRs1IXkNgoiPSx+V6afrk679uNxtD0y9WO6ZWW6Uagi0FJh5zpeDpKLOC2Id4ehO2DaYPNi3xSGWXMScxSE9zF3jjspmgf+GSBdQ6RI5mRCI4AIoSPfcxAgXX+8RA4NIB1YbnrLg5hpxMcOX+5fm5TC1OdDHba4xjyShf1hU0n1hjUSjhe+HXBNhd1NgGdPxRFX1Broy94R9ASGwL0FOY5BgmGjXFQiUgxLDFcOdZxW7Yy+JKozf1inC88kmYZxzroLIhsO0R6XRBZYmN648/W5H7ziTE0/myNY3paTBedxA/1HtUfkvD0kHlGitSNdHqfCP/5oATQ4Uw3Ol7gdeF2OjfVxUSyHfLpRxOV5q2d7MIPQyjpHLhoBU2AFYESfptvwDp2NwORMNiFzQ5TVzCd4AvHt74QhiMYWOfb1bnOdBXSHQA0CFOdA51WL5t5TvtQk13DgQ43+oGBd4W3PGTqHeFy0R3opTdPyYK+i6FRQSOEjxnbuLLF60BHL+vuF9rUiN0v5IERk0VMohCmZzB/9bmMaR7nvpeGn6zP2fGZMTT8ZL1jekUynZsD6Gx9HIDpNBzxgii+aIFxiVMpnenC5BRk5zMniowZSNM1FO7YAU2boMlmh5FITNH9MBwocKxbrXWBdf7glGOdmK4DHfQUvhe+EBLWoRmmOnZf8HeLSCEVJBwvWB5sCsk2J20w2KGK+1s40NG/MM8jAz2k8bkHLBLNdZRjqEQZV7Y0vNF0h5Vwqc9Sin8qQJ8smEc0v7BFmKakPh8xVVPO35C8cTI9UfRBfe8LY0gUfeCYXvFM5+848EFAMn6Rizaeg9r0ZS7jf/4xAPoMADEFNOF35Zj8fJ5AprMh/5ESqwJN2lRw53522hANU5FsdsCdaA6mG7FifiMJ7nWYcjDYybbNgOkC6zpQjN55gWDe/rwBdW3wp0dnOjYsUo9oT0RV8qm1aHZ4wHiD61tc+E2SdauiDeVgNK46ZGjRKaQUAs/Iu5V8L/jkvfhUgHG+UCTNLPKqA9bG+VjpmF73/aJ6n3vGUPf9Isf0CmM6jRuMHt0o4EyHMx2uFXxcVBfw8UW8Owo7UWAdW+giCmQJGhOLGUiz18h3fkue3BsjbHZutnO4E3c41vkGx0h+GLqEEKbDFha+F4CYw4WumhoEdwBw6QiBCuWOEX2FEHgC03mtjL4XdA3ptJvnieQ+dP1JdUhTi3sjuYhSn9IjjZCt5bS40s2NPopwCfqp8Bj0Kd+AZPxUgD5ZEIMtj8C68YVwAXQxf22z2Oi3ic2fnvXex7U2fWUMWe997JgeF9P1AYHRQEK6TMeuc7wgyn/sWMjiewCc6bhbp5lg+8+9tCllKMHcC0f8vn4ApZGdO2RgtgPoeEaHR6bwA4S9aAprXcBX+F44PfEoknvA+aWRDEihHfQWwylqImSkFQKHxHQ8wk2L6bDQLUAHzdNtXt0238dpTijHRnJazrHy0aVhYKAd9CYKiUEuo0AZcdfFfS/8OwFijvBDIjuwTu4XfQpzNBvnr43pYu7TYWxMV+98qj7aZQ7vfOqYXjFM5wMCvhf9Ro/GEL1PBKaTa4UegeKnLIUQ8mF0cEpYkeKQSBfxP/JiWmLuYTLT3KapboGODnfdIZOS7NKEFN8PoPrgeZrRVUJY1x0mwsQGiCHgYkMEJNYFYaRnULei4FeoyN+yu/SHLOgZhnC2pLtY7nMbRI89edvyxQmXj2GAgYGhEnFo8V5AXqNAKbGjFHY6sC7mCA4J7uKNU5qVtikJjygmMqe2noufhRwf09/+XH1YbA5vf+6YXvmZTkOQRiT/NUsRQ44X8TuZnFZiUnH7FBstxIM+2yHyQiefdZjPNMNpwkd2ERiNSvh8hc0Og30fAAE9wlSnSoIFuAreSpSGLkHUX4cyxcBEFYItvTDzqay06kaLDZqUXhNL+rj4/ha+y0j3tOC90DDDnFvltJvQWH/cavAL4bcdtrFkjBddA+8Txhs0U3a8IoBvv9BMEXOEvO30uql4O6lKMf01T73rm8NrnmN6dKZjvQ0RsGKLNFjeSaCzSMwFsgjoASm8LuRaoWGKHwQA3ylGMB0fGzHOGX37HV5wjyJg0x5XLuYhph8oKUikPc3DDhmOdYI4mA6BO9mbBy7jffwwhHXOIKATFUOFEcPvaXSs4xVQzm48BrQJSEyIFyYtaiVQCGZxnHE3Wiqg446Ho5y3nmjAJsGj1H2WxpLgdzPw5JlQrt/x8KpSk/JRgaEC7KY1wLDTn+sUMqUhtfjRV7LWxRzBIWGdTPVwfzrcL2IK425bn8ti7uMwPjv9355a5ZvDvx3T90G6H/qtXfRNiAA68zS20YDEXADT8XR0RvCCHLwrILgQ6MV0/PwF/9wr5hWfTkhPgtAWcsgzcoXGErG0iCeQxEf4ZCK86Ih90xBok0bn4AOwhPVSJJF1STwKgRGnOcnEVqxDguzcMNcfCxs3YvM3ZVATuPht5q2tYtzRT63HnlWQhU4rIt4AQHNBQLvRg1Dza7rYHoo6Y/smv6vD0IIghgQfLSGDSpziufiQNo4xlMhzCYX8kMiOLY/3KxXiT6+UTF/uqZW+OSx3TI+Z6fSYhQNd98RhhecohyyYTm50AJ2PWpuMIQ6Bp9QHN25I4XAMF2DsoFZcJ5dRLtWE+M7hzsmusalJAOv2ybeTCElH+gqha/C+DNmb2A+zD9YJTCEeGKCTo1zIOtaxaxCO7N3BR7LEf/4UkdcBQBcF8UNUjK80HOj4fG5yObQBHc2FF0SJ6djZkhM8sv7+Kai+FgqUw9xG5/Ie12U+WsLHlTiLjFwnChWgpwHGz9pk0oYtj/SYtOxMx/NVMff5YXx2+lJPveabw1LH9PiZznuR5LLY6ZzpNExBagi24UvxmBKYJJg82AmADV4pBWThXn4otAlUNLZaAu7kTFivFFFyR/D+YeAX1gnVNUnz7r6iQLTqGkDfbK1zDwx2mAjTmJMUMkcqRQqe4g6DlqIkWEvJaHNDc98FHopynw/qIAQsLby5ks9F42kucrbQjcV6paih9IeQ1Ju2vhbxYsCkHGBIgIxQiAGMIR0+8jFBIGBGCKbfo1S6/nRYZvp8t8XEx/TFnnrZN4fFjumVlOn8U4uzlOK+cgxQIfDxjUGPaUATg+YJpg1/I0O8bmc85OnJ1w9VKQWsAWTdk521TCnC2RqlaAsHe9wHTsFUPzJJ8x6+6pGUyfyEtS79wruDn4PYEujnWCdEiv/AKGc6gLsq+el5rEDs3mKf3X58+7a4V8CKYiwCHnNRMQ507HVhQE8EP1TUOvk5HVr/aOUTbQUjnR4vlz6EANCLgl8EpPUGn8fCZhIxhFL2O17pxOAxDi0RicR8mGEI0WCGoUCjnY9/MTX4IbAOlzp9QiADpqeL9fiYvshTy31zWOSYXu5M1410DAX4W7hAvhfOdA5B3Z9oHK94QKSjnGYLphD/Gga96GT8z5NR3aDBKBjTU9F465UWKnznnYhJLoUkrTiqugQf+SKsE6dAKyPWSyErsI4dhPCBcIAS1gm4/CEq0ZbWHvg9kvUsLSu5+aTUYOdMh5HOnytSWVQurwlkbBC0uKp4Kx2YqpW6BNxvHawBpUCnJtoULFqrgoX2ZaXwy3DAHx88xi42DgaKFOnxWcTw8QaFGDy0imBrLyc78ZrPCHj5SaBT2BvDt75UKaYv9NRS3xwWOqZXRqbjqwD8Ky5wIIphKvybZLPARU6jFrYSn12YS8ZVxxiJLBDEREU8vqEBPXRKbLqn5Yo+OvZqcLNPjgWv1AmTE/yyGpnqMD+J5j0Da50MdmAdz/1KgcV9xPBr0/4T4JIEbkQD6CTgfgI1ZJ6iSEzHXQIHOleOZYPXCvtksNFl3/1CAHpaTUSPlEsXHi9wf5G/5dXgF7IWBXeHc5QSn8SifkS3ksD7PWRIGPMKVfwQakknjV4Od052bq0Tu0NmChk9nOl3KgWfOLexQp6RwjizOVtEfHx2+tOeWuybw9OO6fuT6foYwgPS6Ha6znQBdNhWNDHEvOLDTgxljGaexvjIl89DMdChk5SI16loTj6llI71JDHhgemc/IJg9yTKe/pKYD3MsQBrnTYXip2FAC5hHdilbXN4rkvbMWGnCze6MNLpPR3Y6VQE+TQI5fS6L5RzBwtHOTbbaBtdCOjNg58JpT2L/D5GtA89UqYHD6XPReFy0YFO5rn+siV1JfUsHxtiGOiHPDFkjBAIOMUHG41bAXfc8/H3qKNYP0amG011VEa/HDHUkdImxMf0Jzy1wDeHJxzTqwDT8UsXNFhTWh94/sPflMOtLg1NDDtMpOgC8hqFlHr4ukXeUnpgILC+JdgcHVASpmhaWJfuBT/YmkKExf5CMttpux6wLrDLP7TAHUTJVcdqpOtM5wsGVgvOdPHmKva280eyyZUjLaAfGayIBHR65PB9tXcEb6IWBTdJ3ELHVr+U27d5jxuHBCJ5yogy8tK4FWTXX+DAc3jdS0nWj82fbtv6IiogyI6zUYT4mP64p57yzeFxx/T9w3S+wmNw07DglizMEOFS5+MVLkJ4Cbn3E5+QhjOESkGh+0VIiXU8Mk06rAXWuTUaxbe+zwZHvGXKP1QSvgubzGrs0sHj3Ih2ur6ThG/x1pnODXNOc1YcbxCbD50/byA3enujG50eigqgz1BqmlJpAb1cxxJmByaI8OOlNNUxU/jz3seUwg/dYVbiQkBqgXIcIkEUIT6mz/DU33xzmOGYvv+ZzkcDjBE+XvESKbdBaIBimFYhoNOEAdbp1Sp6rwpOmJXB1ouiYMsKwzo5GaJY6/AztA/c8aWuBjJyU240JN8Lf1+RyEs7GsH0iPtesOuG7HQ4XmCqk7M+fLtk0jyHM4p2ueiPjoXLhSz0AwP/TPPgZdHvLXQ8F12jFH4tiH4BjoCe8mUcsK9iBJopmCN47CSs9ZBNYnjLFPeywrMErPNZyV1AoDkJIln4YXxMf9RTs3xzeNQxfb8xXR89UQYrxiuchriXFNaHGKw02ipm7kUshaaKeHKwMNhu8XJAmTXBgzv6zftoThjgrHuwmR3udfhhSg123RVDpjp3j8CaxpNMMD2KS52/W0+2MLaU6Jr5TQDeaaKdnZp5bnSg07WjBejRMV/b9nG5lATrJW100YEO65VGacQOrZhkNJK5ta4zHdOEGz060Llzkj9nEmgWHOeHImX4YXxMf8hTf/bN4SHH9HJnesgizwcHZLI+uJ3OPTDkhMGQjQj0Sjgzaf5PDF7LBtbx5IA2OJK1nmp3Iz0y7emr3izRQxHZAAAgAElEQVTQs1OQnT7Sa/gsDMcuAc643ZCMdNqZg8ek5PEnU51h9zsruCQI/I0nMJ0767FUYMcL3624r78+EXzTxninYrt2AB1LmrTQ12sWunCj64/xKwbc4aUIrPOvaMADY5sjsNBpTyRMH3iZMBOjCOEQF2fjY/qDnvqjbw4POqZXIqaTpcCZjsEKrHM/jHixSHejVyofunGWCg/MY0rhBSs71uuw92u6BPZ492ADTG9fncxC7+SumO7Blm3rS0nAOj0v5UzH7hejy5t7ukF20Jx77Y3OegK60EwbFvfdrYgP5/JXivhKpl81X8yiAv2poPFhutqeGRr7sYIjwXRY65gpeO8a1g9294r9YLhS/syJW+spmS6QnfIwPqZP9dQjvjlMrZFM37BhQ25u7gsvvLAvz787KvM3vGz9GjI+gHISsL1Xf9uCv0eHF0RJIKDrW4nh4angWRexODAd9yV8dyPHOlzYgRNGx3qPwEg/2VcFLJwcRJLNLtzK+3jYYVDr8OVvewprGh4SQjD/gBftrrF90ZDr1D3pMP+TV4od+vQg4chgKyfZ5sbr7ZH0O2GXS6nHaXdylwv3odNWdNvOxYhdWcHJgHXhVYf1o08QMUewxxdMx94wTMmQaWub6SHx8TH9d576g28Ov6t5TP/ggw8OOeQQpVTFMt3ogcHQETQPwToNWbxqgU9h4I1t3EvCzqqct89AgMA65iTdOxuxHjg6gHXa9cGZfoavBgbhDF+dESCebHYY7IYHp3hgCK86OdbJVcI/D0AvBJGThLBOBjt+TJkEbKohhXyjpK5TuOmTj4XxOBTvE+GOpCC4NHGltIb1SN6XkIX+PdCxDx2v/tObohzocKNj8KCbKqFgwzo3ffQJQuY5TROYTZhxJIiJqZM9BNwhp+Jj+h2eus83hztqEtO//fbbRx99tFnwFx/T0d8hfYlTSEz3jDaUixFGNqzxJToavvxlUb4xq5Ib6YQJwXT+vAtYp8/C8L3hDOvtg53XBLuTA9IN9NUgFgYGkeAdGezA+vcuZiAPHhh8OQAfgcE7QSG7yEN2vqdUiM/IJL3zsNAPDExvrFuEcuM19gw+g9NF7HKhq+Ov/vMPudi8dpXfGoCpZLTW9dnBUa7PMj0GcOczF9M5XSE+pk/y1J2+OUwyM33z5s0DBgzIzc3Ny8sbPnx4cXGx8FIsWLCga9euDRo06NKly/z588XZq6666qKLLhKRxkNljC2nyDfeeCM7O3vkyJELFiyIiem8pwHQ8J5GFgF0mAwk6MOLYnSy4yNKGK942kM1qYQWFq+SbVrCMUobHLnBTk9NA/Algu/xAnnE9EG+GuyrC5JhcMD3gUmDHa4JwnrpJwTggcH3A/DBAHx6hQxqmOr0tidtJ9f/411QvOvP9UAVN/m1TT5koZO/pXdwCbRi6Rd4RvAsoWdA/wODZvn+0rDFBf6WZcGLu+EbpSrtc3UxeLhNAA8eJgUehPLZYZxctglI8xTTNuI0N0IgPqbf4KnJvjncYGZ67969hwwZsmPHjnXr1uXn5992220cs2vWrEkkEvPmzSsuLp49e3b9+vU3btxICbZs2TJkyBClVGVk+ueff15UVOT7/gsvvJCK6SOVGp0MY5XCMBJdxXs6Ymcji3gWKoYUfHywFPgawMlObkE+ZPmjUdS8cgrUntQmuEBcne6HWZbc40hb1wOPc5PAJu0SeJl7Byb5oIDmQ32FcEFA+UFJg52wTlu2m/Cvh/PvBxCjicuAcogAgpObJSQlnUJ6/iEX5nWhz8fblitc2tDgYvmldQ0apEnQON/dhewItoQC6MYvc1Gbw6FMPVJFbQJ6Xop5wacGLhDTCnzXJyCfg5i2ggBpHd5Q9l+V+/6B3/WemuSbw/UGpr/33ntKqY8++ogwPWvWrI4dO3Kmjx07tl+/fogpKCgYP3687/vbt2/Py8v79a9//bOf/awyMh01jsB0xf56RWZ6FKwLfsGvYhyCNPIwnuCx4TY7xiJFcqCH3zuHDEe+AIQkC9fPldhkUo4LBNZhcMEVI772RVvXd3z3q8p1greKiH1kyYLpl/jqsiBcEvD9gsBgJ3sWWG8d7BHc5ymi+H5ACNyF3W075HwnlJNOevWf759hz0VbB9+FJwsd9x8XBBdiu67egePlwKBBvjfSAXT+mij/MhchD4MH4weD2dZ3EeOjj5+QlCFl8SFks5OMV8dnE8wIMQ31jCGVDD81Wan4mH6dp270S8P1u9R13vehsEhfOZ588slmzZoBgG+++aZS6ssvv0TMwIEDCwsLcVhYWDhgwADf94uLizdt2uT7/kXBHxKECBXqe0E9IjDdaKfrfQYYQdDT6DG2wSSe3sCFgtnF3RQYvjTsxOCjQm0zQa9SLDG24mzxKBSth5bB1cHmom2Oc5Sil5JeDV403RJ8wiWAcvvgwWDPpHeC2HeZr65Mhst8dYnJpJXftKLd5XwnIrnIdbITweFtNwqUhhvmoDl++m7fF4vwXLRzcEXkQz8jWI0uCC7BdlEFwd5NMtJzyEgvCXa5rAm+nbssaDraz8dHGh85HOh81Nl60BaPno1XCC8Og4fPCH51NMxQJTqkkcaZzhsHT6eoZZA3XYHKio/p13rqer80nDiB2aDfiZ7nAXe+78+YMaNDhw6IWbt2rVKKnBYU2bdv3+uvvx4Jxo0b17dvXxxWD6aPZrY5RpLekRxGkPVkIoYPPs4s4fsD0yk7VYNkDEcYtriXpFM8C+rPNXA9qDkXKEEG/0VxxkOulhfKZVwRbmXIFUNYJyfM+uBlSGaqdw98ygMDT8vQwEK/0lcjgnBl0mCHtU5Wbdfge+Ly50wF2fVdidjGjk0sYoci3ypDrnPhN8fv3iWfiNI9Bz0XJaB3D/Zi6kDnVzQ0uNiBwYV3Dx4Xt4Y3yQt+CWRl8O1cvPRPzOK0EijnvcB7ytiVIpKnT0vmhQqZ6xHF0SEloFx82MDTQqcomcgCppOFDkc8HfJW4oscr1KILC4kPqYXeuo6vzRcs0sVet+HXxvs9Llz5+bl5YHRZKdv3boVMf379x81ahQOCwsLBw4ciMOaw3TRYTg09jH/bi0xnSwIvCuBwSSmHI0kMZSpCCoRgxgzk86KLHwo8+yodsaCfr3GohHJ00cplC4Q3hj6MswipV4OPiW4yWyqD0r6KMB0cJBb64R18TMRpX4YQXYY7HxjIqEcLxBBoHjblvZ9bXOiOf8sV5fAi8JdSZcENxx0FbiWy4LLHBTcndDT0c7BC1nfudF3B270VUFDLQp+phxfcRGwS9kL6DJ0olFAMggpNUdMAIU2dx8lIG0RZwR/uBo+DYV1xStjk43XFR/Tf+Opq31z+I3Bn75mzRqlFHlRfN+fNWtW+/btObLHjBlTUFCAmIKCgnHjxuHQMV10Mx/9GHm6dcD3XcHFiTs+roRkqBKjh+L19DwXfDgib1kOxVWHzz1b5W2mEFXszuATetOCX2mgzTCvBr/L4+3jVeemOvmdR7DRPyIgI7BeEBjC+r71BExdODH449PMmI53lPb98he9Jsq3LXa1Ax0zeUTytmNwsBP/ZOZJ/35N8oLGeTX4fs5TQaPRuCLkUZPqvWYbG0hpHFo4C6EsY8mYF5qNFYDNpOeljHouI9P5NDRaV2Jgo1Zc0OtAMfEx/QpPjfTN4QoD033f79mz5+DBg7dt27Z+/fr8/PwJEyZwZK9evTqRSMyePZv2vSQSiXfffZcnqAb+9Ci+F1vP8d4VI4BGHux0WAcYSdhpjpu+ySYvEAc0rwaK1kcwz2Kbt1xVBjJKhyCqgXhbBXgCajoeQ+1GjTYr8A6TB6bo+y+t5wT7PboGG2AKmAOaTHVOwyv39a33DrLgC+O0x5G+eFVqs9OnHMkJo79uCvMcgvEFUdqqmHS24KX/5sHWQ/6aKN/AAwvddgkFQf338aTvCLwuq5JeF/qKC//OIhoWNESMsetx1jik+Vlb5xrVRo9EEWJQ4RAJuE6KRBohpJyJKa0rFAqBly7k+Jh+sacu983hYjPTN23aNGjQoLy8vBYtWowaNWrPnj2+7zds2HDmzJnE7meffbZbt245OTn5+fkLFizgQHd2OnpXH/2YP5ODL1MLpot3l3HHJwYiDlEQDR0cIoEuII0YbfrtqkgQ5RDKuYA68EijNp4AuSDQWWD9scCZQB6YNcGTwACUTQLn+JHB00LCOnnVYaqPTtrs5F6HM5oMdmxdN3/zC1jnG9jp0Sj3wMDrIj61uC/Q9W9y0Ufhe1oeCVztK1SejPShwbpFj0aPDC78+/2LeDT6crAPfY5S9IIobvtCRiba+f/b+x4grYpj3yEoQUQh/FFeRKIpEuFCAqLmom4CKDHULa4QH8klJfeqURPNoyIBHqIuYW9EBQMJ/imK4hKlAsVbzUYhIiCKkIvGjUEvJKjJJqABUwFE5TPg6oLM4+ve86O3Z85859v9ZHfZob5aeubM9Mz09PmdPj09c9w58k6QzATt1s2YE9ZD5g+VcAnVASTdkpwDTEd4Fe+KYBtL2unc/4JyCz/PSofp/5HLL/h7f//hx3SF0R9fsnniXgLjyXzeS0BHoUkBDQA2QZnkJmYoE9+EaRqJJ4Rs0dsoOHBJ1Xl5I6XRqoo3qboRSBasjg4rgu+Z++hzDUvpg5kbKGK94WLpObTvhk3d8eRxvomWSRkWpyfgqPwwOEVAIrs+AMuFdUS/MJQD0BHx4juWS7nOcYSLjFmU5jn3mTEdgD6ewu2HJ5uMesBfxF6XrQ29LvcZ412bkRLGlBWcIJT0Et7qKjNN02S+qsJtyQ4r2tsZVUYmXUz33oYQnffmUo2qPstk6TD96py5zvp/V0dMbwjwHz+mS/eLtNO9ygTDSiqipJU+hZUehaWeMS1vJG5U5XDSrShzwD9MyCqKDvdfiY6jG1eTe+EVcjU0DFdH0Mh4smgmkoU+3Zry5AeIVK6YkWTmSye7PnkcnxhVThiEM7qBLuLgXAZ0dYRLAM1Vh7m37Eb3B6Tzh+g4eNH1ukj9UTQmTs0LkigQJlDeS6SplpuvqqNR1W0kUUASuKoItg+wSIMNbnwnquPwAo9D2dbxstP/LWf+3fp//xYx/WPEdKVDSIaVSR2PF1CmYu10KJ+6VcJJeaeFS4J/mAgwQUXIyiWUqV5F7gXEwBCs48AArJdOoBXFSWShl1tTIX7lidnuIjv2JelPvqnDYbBkyrZ5ZkDHjn8Zrcjef372uP1kI32CWBftJ48B4B1GbqwLLE1XnsiB8LNMEAp7iQAHNo2hVOGS6iraQp8VgQKSUGWQVLehOhdPYnpRppXqs0yWzk4flzPfsv7fuIjpWTFdakn4aYyS0B5FKGVSBgIOVOJl0mbHdKmUARqjLkhkYaIkJpNKeuyBkaGNAtb5wAA41icmmD6r4VkZswjiy8mQl1YwG+yDaT+nhvW0w2GwvQg7/jlsMdkg6u7452h69rSwd6jCGm8Pp1PQDtzoQ5No9Potoxy8WJMEL6pYl7AiSfsgywQ1ZZYD/MOX0KjUB0mjgCRkAUmzInGMrLwH5UczECkUlp5sLjCE0mH6mFyDU+rkaW5jIqY3BtMD08Zzz3MsFUjSCpWgT/IkdHnqqdeLxwylMoGWbUkaBbgD4VEUe1UyB53WyQBzVVf2X3LDIjM+oyFhnU5NQRhMGfmd2VSfSo6XWdbMtWa++M0llJ9FV5UtLE/F0qeuy71IOFNX7S0SR7h0Tk4c4+CckWRuy3eICuqG2zfu2FR625hAwymjXaa9aQ/tsWh0BvQNtNIArwuM9CyKFJgaqds8TXJGMHFMhPk07iqacLUi0JNAYalFsNNbA6aPzpmvW/9vdMT0lovpWb7aDi1XhKvHqkDjbqq0Woo5J2UfVIE0PhI1ZHVF46EojSyOWIcTJoF1rJfyLqRJAtMftGahNYvpt9CaBwniZ5PNDvTk/Zl8Ivk5yeeqjx3PK4PWedO/3PrfMNBFbhAdmgS34EnDaD6fuuHtVTkt895AsS4jaV00dXsRvlgkI2IDgN44I11OippfTgZmuRGXVBOydaZVASTdkijP20RknIKyq1qqnT4qZ/7V+n+jIqYXjekF1TGLMnnjXlx9Qjhj2g2JthTh6rEqUHAU2QsozpwMdyDMHAxdJrgbFazzRiTlWydY70Gek6GJqc7ul1kE3wutWWLNsuS3hMD9QTLhK8hLcxOdDyO33evjedmxDmudz+TC3iLfuig2iI4l5jeRw6eCGn2QOqC6tJC6Oov6M5GqjKZo9H70gMk/XY4U2i+KJ2VYngXf3jAvXm2UV0GHJ7qoq+DJhDsWVQBJt6RcbOcYqrR35ZaK6cNzDT7T2ODjhRHTM2E69KOg3uP+8eo9IMkbROViutwk4qqm7JWXRhXv1aLuqEBhlznaVQRKBrjxJZT0ihFXuTA25SIShq31fXm8O4ncHQNpv+g4CoCZSt7q+QSgy6yptKYq+VUSyksYZVgfTdXZsd7g40G8BR/nOPKxXL5zFqUbHUGWDOjyAaM6s4w6OZ86PJU6P456MpAGVb8fah+dkVBNm7Ckv4VX9sIKCUkWpdjeSeG5lgyZLjjXWQq4bGUfvFdlptJDiekq9qzYe9AddWA4pfOnfzlnRlj/78sR00uM6WFNggZ4MZ2/H4010vAeNtlQyemAXqpLbtPu/YMcWVjxcZOycEFaypO3mFbTwY10bEBn2o/Dpvq15MGoION3MVnoVdassGYV/VYQuC8jn8z8xAlzAx2ShVMPPR6YNExP9ot2JJuaj+XikyM5vHIqNTGXmltGTaMnq5KeLKauVlC3r6W3jaH05lF/8mKOhglAV/66gnJDAVf+KgclJZhiZkHIYkwrPoGkW7eEOeghHE1wpkvfixtSXPBdGQwLjrd0mH5xznzZ+n8XR0wvDtMDGsmXoIVShxQt/QZhfVL+UGaudAgWFpouIREer9uQGqmbRJUwZ1xFeZdQzFkO/CpdSQGOmylonTwwvWhFcSSh803kUgeSMoyuteYZa9YSslcR1i9MnDCTknOy+EwV/izcsS9pYBcSjHS205NDXfjs3F4Uo4LYSn60lFMTC5NHyyrqALqxgvIXU5lyingZT174Y0Z6HQ1wMw22kvZhYTnUFY4rQM6BqMOErK6Yq6QsmaUJt3zTc3gsio+8cWAEKDu9cZguOeN+dOVZOky/KGcutv7fRRHTWx+mK02VqqMulSQp+UvaZa7ubW8StSQrvsHCMcuoyISXOd9a8wjdqmgX5St0bAAFh/QlN/TYxP0ym1Yjl5GRzoD+DMH6MwmsL6ECs+kBAEc2zj7sYU2Dj96xSx2wLj5Bx9/ChpE+OnGjT6cQlwfJ21NFjaIDDOuM6Qup2HTq9lgaQl/EurDXZQMdk7DImHnGBIxoJUAlf05iCuC3QTFUT5O8zEdhEOCjCBQoISGbSGN7/DF9njF3uSebN4Sfwqn6TZHn58yFdGqb+/f8iOkNxejbRyrVQqqLl0ZhqeKK9trpCKVK871IxyhacfuAS6UlVENe5mqY3qSsKHlKNHEBBSVR3cucM2fRWToPUzxfNZ1NmJjqg8nOnSDcL0vIxbGK0HyjNfxjWGff+oMNfdk4/rCXxHRsK8UaKZ+NTr4Xdrz0pTtwJEWt3EDrouxGX0wOfbbQVQeqCO7ZBTSJngSjKNyl/nh0NtLZ6/IwDTkA6NKQhDAVEZgCiD3chJwUWYVp1ZxboCQ5qpXAXYNlGPddWZ73ksX3IsXrtdNZtqXD9C/mzPnW//tixPQiMF2pC8+TyoReSv2WNPRbbkrGsjsDutzDxioFmAN/8FEd8CqxqtW4pGzI5SDHGKZRVzJUgILxyjJMc/VwE7PpVN7KBqZ6V3JDl1FI+A0URjI7cWSvSJwewPS1ZL9X0hLlXPJ6M6oi5uQYpnM4I071QnB6YqozpvPuJzbSJybOHwb0tNbZrT87icAZRycW9MPx6DlaGt1gTCUNNnB+J8uK5ebKk3PC8seUhcUur6IKCNk0MktLyCZAqyaQ72I67CqJ6TJOIfBIk62gCSYg29Jh+j/lzBes//dPEdMbiemYJzl/cl6lfoPmAu57H2M6TtzFB4/gTwfGySZk05KWZUpLBx4YGGNBAl2SfVbuFxavKsBJrl6wlXl0GOFq+mzbrvwHIk6ik3gH0+cjJpB7GiullcIDo2BdAusN9DwoS3Zv1oed8Hb8HYSwfN7LK/RysLv+BGD+XGo/QmT2/LDXhd3oaYDOjxNlpF9IQ8i3eyTxpK+mYd6XfiAzBMVy84oUygxCFcOUgVtBAlUUEVAhVbIRSdVtJCUrzuSR4kMrfMfhBsSR126cQhqsyyaUqQ6plg7TP58z/a3/9/mI6Y3BdEwSE67qeDWeZx0qBX3Cqx9Uis9llPoETOfqSmnQARBKw45D0jtkb6bsDDpcckyfnXjVNx0LgOlKcMym+rWOA8SFVw4+Ycd6RbJKOZyWW3vhWxn7CMG30o78TXSgWDUdFSlOAO5BVYbTCu1EsvrnN3SjqwcJAH0WdVI+S+pP0+Vwl02JJ72gkQ7ngBQ4aKXPrrJhyrwT6s1EleNGYDheQt0yPOSi7kFmy8MpOGT0Qcq2dJj+mZw51/p/n4mYXgDTpUZ650neACjsTjlfYg5efZLfrsNZu66djj6ECfQkQLidxJ0fqJV2ycvNmyk5qFHIG0AKVhbj6l7OMpO96kspJmQr7co5kl/V7J0cw+uNVffCOjwwU8mpPZIAugdjOn96gj/4uYHaWkd/G36sgzF9JFWfSmuecKM/I5z40uGjYtJ542hvduLzJqOtSbjLguDSKGTCcpOSBB0Wu5wvcCtIyFpF0V7OWThgOFkIeQ8i9MW9B712FTqjuop8+fyQsi0dpp+VM2db/++siOlZMR2KIieJaVzCpMrJRiYXgzIpG0Hqk/y+rfKno60wgUa9hOyeS3urFMx0+aTlKFZqIFK86hInuXoac5nP7peVZEHvqPeEsFedYwoB6+yESXOGcPzJg+QEv4FcN4OxhxPh4etoSZZ3rOCUgpr6T2Azpo+i2JVyijfneJuCgH5TEhp/IS0G1BvpB+nNQK6OylGn0Sw3r0jDMpfzlcbczZe1stMuH5kT5uMdWlpmUbeha1qhJ7J70iRCu1K2c0oW93Jmzvwv6/+dGTG99Jgupxlzz3OsNIkNBEZw9XcBLXzxd6i5FrQkC4F2XUJ2L412axXMSWPlzZfcsgxHluG6XrYq0w1qpG2lOC0AUSiTkr0/WLRktH2OjOi1FBuzmAJgbiKX+jFMRzQhjljhpZGVtDy7tT6Ssge9HIymE7hm0cJsFa3KbrSGm0CkDUejVyQuF/54NM5frPekc6NwvMzL4EwH3EhJZqHlTCnxhpOyYkY6zJCvBlhlGQ7KNP1ORE9kt5GpGuL12NJherec6UGGhfu3W8T04jBdPnVBY/54Rt055gKuGqVhugvoaY4INK0I6JZLyO6l0W6tgjlprLz5ipvqfDjJdb1sVSaiX9bRSmliqvPp6nzuCsO63AG0mBB8LaH5cwnmrqLtP3MpCHIcRSXW+152E+d1FHzCR2kuoEVLbHranf/+NWP6WKo+mzzpqxryx2MDi6J8oO7I5GiXY2E27O2R+4yyONMbh+lqmpR4w0lVN0syzJCvBviE1UZdxc3LZhPfcd77EZdcAwudQc+Rw82pVkqH6V1z5lPW/+saMb2RmC4nGIDLM+qd4CyADu3hdz1vE0o105LQLUWgbwUJVbFgsiBDVUAxTBuIm88VFTdvkjGdD1jfRHEpu/MBMIa80r3ImzGYjk/h47Tk6uUKilh/zppqgvVnKIp8PoUVjiecrV8jZUxfbczSZBsnu9SW0ukrm8mPT5h+IXlRppLjpZKMdMkc7h2cATCWOsYul2OAziek48vRS7NFMbJwWG6uPNNy1AR5JRzIVNULJgOs1KU0VmkD8eYrtOVZww3ofWnGXcl1mS13Bj1E3/iqvOvvM6Z0mN45Z06z/l/niOmNwXRGW4m5gQnG7HJ51ps0pYHewI0uFYhZFfwLxXIJKF9Bwq0bzinI0C2gGBYcFxfgWi43N4cP9OCvUW9ITPV99bDOIYZqe+ckcrAsTEx1ht1qsqlXkM+knNY5y2itNX8m4j5iy5iOuVtAEM9hlGSn96JAxgnkkV9IcZMb6WnBDww20hdS04iC5+PR+WyZ+g2r/OXoHdQi7x3lrUYfh52upiYtgM+VOXJcDuEcVCxIBPhkVCEuBlhnx4i8PQM3qbrruTPoM/omm+AqfEebXC7XEHKKS9VviuyYM6dY/69jxPSGInX2kfIkeWdIThUXwIziVZdXwKE0AV0BIkilaQSgB+J/oXlZCIwlO5GFrSzj5cySTPuLKpJPGs2YvogcI/zZ0q20upjAurGGP1KBzyFdS5b4XOEeqbZmM5nqaylzFq1zjiQbP4/pOQpbZN8Lx0jw+zv7XhJ/+jlkdN9AqL2EjgF4jtgyZzh2ptPBMnJbU/1BXXz6I4Im5YdG+YCXNAnIfIguTbZKh1H+OAB6sU3Ivik6PDp1VcI6blJps3ttL9yhUmIsanQGsblACd67VCJMb5/L77fw/tpHTC8Npiv3i5xaucMo4LCDokg9UyqYJYmmXULe4QVpt3rBnII83QJpPNNGKsu73FQOMF19324HuURytHOH/DBYNR1Ly5gVyTLmM2RNb6a/7H6ZS1HqoymcMQ+4dRT8vkl8VwhrpBwXX5cPoOxH0TIT6TSuZeTV4UcFvwFUUXOzqGk+0YXPRj9mnueow7ytqZpCGPFduoyYLuWWRbayvJJqlqSsnpHOwhZlAjzTRhfIl3cckN1rfuGxjddorwnFbTFbvq+Z27ySxb3kjQmOpnX/RkwvDtPlDPE8AYt5IpW2AdBl5KJ8+CsmrAcB/QtfUq3LJO6H7ISsnpHOzhwlC3LGkFVJcPAS/IbEYse56pI+qrMAACAASURBVBsornErAfEO8pzQx0t5GbMfWdPjKeZkPi2Krk0wfTO5X6roVK/p5BkfSieY17tf4OBeSeDO7wTJ8WE9aE11HLF9kLw6GxMjvZoc65XkZOfg9+H0AOhBJlie+UHq5C7qMG9r2kCe+ipaieV3O4YVrxCQqUSXJlJVrFjzmZtzmRTMQT+zEwGeGF2xBMBd3ubybmUa92wWSw5WPx+DXCI7PWJ6Q9wOpdJ9L+5Tl6eW7yuJxdA2Vil+8sNCRxy6qx9QkWJ1kcujXS+R/W5BSS+fcCbqFkWEeQauBlqB14sDUZbSlsvVFGK4SSD7LvKf0GFbHJ0yKok4XOxzkiymqMcb6CywfuS6qYf1GnLCVAvO9EWOjuSlGUlOlQo6PWYVeXI2Jy6dVZTJRvro5IiuPM8j1DGg+WbaoQpA583rGTE9IMDwpYB4A5fCPL1XA9zSLnn5ILNxdxBuWNyzgduWTXX31gMT98FwX7TTQ+D7MV3LgOkM5Xh6y4c2QBkbyfAkCCuH1A+pJUxn/AuF9hJp90Yg38unYGaAYfhSQc7eAl6eXJJfkljybKpXkZG7OtnqyZv42WCn0PVzaDFzPHnVH6RYF+l+WZtENE4lU72M8Lre5Y0Pk+LwdEL8cyhIZhxFMc5NbH92pleTEwa2/wRq+hzyBeUBHd7z6gTN1yUWOn/MCFscoTABUXhFF8j0ssqSGeAZuJSFsyoT4MaXMt41cJ7I8tIDw3e6a4p5X9CZibzr5Y2/IGL6x4TbIbZBTJevUcB05WKTmgHHCx4D0IzAiS6usiqe3qRbS+ao+yFjEhzc8rjkJdzyGXO83BqXycKHwcUfLOUvSVUlHpJ1BJebybOxL++g7EqOcjbVEUgOCN5IISu8+WginZrLZ7/0poodkzWrjpTsTazKqNhEWh1dLCJe5LrrbOGjr98puo+6xLY5Q/lKes/AJ+0B6HhBbJyUvLUyTpZbzMsNmYHy7qUsOeDsJby3icp0K3IBCcoK091zOxRPWVcCOm9fKJXvhdf562hFR/3d1/QnRwgiC10zhQoc7+sZMD1tjvHchrUuZzetFlw3rBmuknGO0htvMq1u47yiHLoTuLVK3hy3FWCb/RLLxyt/PiiNP2fDZjvCHMm33juJO5xO3vOqZEkTgepV5C2Z3XCrZxnZ4/yBgqHEYSTtOMVJYd74SF4d5fjI4fSBvfxOUd76v5ncRKsJyrm33POirIHsEkPJwIwXvAQmLhGo2xT9dBviHO8NojLDdWEN4D0PBhnfy7hzcb9LMwLrZ1yL525R09E2Aajdxryf8tvd9FaaArutFNPl0bg8Z7CbsCAOQMHs8vmLKK+q4B0wrGpKL2UyrWLj7hnmFrgVZbym23S4Yviqy83NcQcuc5hWU6CMJp7ESrLZGdZ3HTPVR5PtrEx1DoDZSH72SgpWmUsG+HTC90lUZSK5WaaS96ZCfJqOP7vB4S4w0pfReV4Tk/MG6o303eSa35DE0uDcV4nmDChQNnfskBhfQjJMhOclfDXAOUvFcBnv1UCLrkBUTqAuVmLgEFcGmcR0TAErm3LaSECvNOahpqNtgum8DpSjRRf1d1fTW2kjmA5oBqDj+xXqExaYZuiEiyZSLfCcZz2T6is1TymlNynLS1ryzEi7nXEryiYU7RbOnqNYqaR34N5M9zbDsrb0hjGs80ZTChHjcHLeIsSRKuxV32LNNmu20PImhyHyOS38ZeoqimnBX/7MKZ8bw1Auq2+kwgtpF9K1FG9zDse6HCSviwyOhHNPdl6+FEJ/vEJQmUqYKpl9jtySipVMuoVlDpeUORlp2YSk1ZC9SVleGj2cz1Vgqqdhunqy8qRIoIDdwG9aP2862iaY/roxb6f8Xm96Kyc2prM1KmdXzROQXTna5AQD01k5JKDzAwBGutJmpXle7ZSZqjwnFc8sSfAJF0YxlwhXLHjVZcg5crAFaYXpmBFJ4HgWfD2DvOGDyXaeSHb0Eoo4rCZAf6POvFn7id0HTt6XK+r3id0HzJu15o0680eKj8Ry66TESE89QEb2VtEYYEFRyAJpsi04KeECaWwlYno5oKL3ajgTdSUhB+ulZWG3e3wVi2HAaHnzygetnBRpK0ig4AM7S4jpr9JHUXb5/r4aMb3BMyl5DE4XZ93Noq/38tRinniS+C/cnS6sw7byKgQs+jRMV/4Nr4LKTKWsnAzfFd6r4OO9ikwU8xIo1gjCyxDRRHLIARqQhzdivv1UPsN6Fa2X1uTd2SeRa3s4HQNQTj6WFQTEfzly8r5ct8Nv9rY1n7Xb+tuXBtkXvmQ3XmzXl9l1w+wa/Mrsuovt+i/ZjYPsC/3tS5+123rbmm6H3zx5X8785Uj+DBmcNJBqpPOSGjSE+yzhQ44iIAT3UppsGzFNqJLGk/NRzEvIut4CgUxZF7Q7ZJWDki6g446TznEgtbyLcWsrAi5W/hIejl9eacz/azraJgD1e2NeT/n9vumtNMDEIhMt358u7XR8kEjOkwpIkLAuJ9vVBuALK5xXvaBh0EKlnW4SJUEEbgnvJVRM6xJqyZIujWKNIFxuxQI67kmJfaAhN4Z7nB+QnLrVNTkdlzd/VlJQ+Rt1XT74+2fttiF202X2iTF2+QS76Dt2/iR79zRbcZudUW5vK7e33WZnTLMVk+zd37HzJ9hFY+zyy+wTQ+ymz9ptXT74e95UrybHC29JHUsNHfOkywMX8Q7HvUXnXQLDyUh4xduIaUIVL0NkopiXQLGC+uZWl3VBh4WAYky4PPmmg8z5Ps0C63yPs9nnAjovd5cq7uVl8tHV+P6+HDG9wWMoeQx67XQZD7cyiXReTWtZjOzsh1FrWRLZQQPQ2RkaVi+phWF9hb0vq3i1Ni1TVix4j6nCKpnWRJZ8xYqTBceuCrjYxzluMcdUN3RKV1liqi+hfaRv1vay24fYTaPto9fbB6bZijl24iI7Ybkds8JescYOW2fL1tmyNXbYCnvFcjtmkZ0wx068zc643j4w2j76Jbuxt63Je2Ceo0D1CtqI1CDcpYbeFVYme0RdR3nGEakBukmveLPMS1oZL0NkptXifBRjIlxYXVV1szz4VRXFUHYJZkEA1qWthhVRRFVhXzFviXi86WibANSLxrya8nux6a00wMQiE63ITmefOG9HZEDnr5Rx+DAjO8KHAevS75YG6ABir3opU929P1WOUtmCuKwaVdXVVZVUhVVSFS4qqVhxUo20YDKMgGgCx8JU0nN6M/koj+SP9xpozWiKY1mYx/RP7D7wObvlCrviO3b+nXbKEvvNNXbYS7b/G7bn27Zj7QFTtz//qz1g3rYd37A9X7L919myJfabd9tJ37Hzr7ArPme35N0vG8mfM5U86QPJfZ/fOCqPYud95FAMQFV4RAUFggIYuySKmiBVWPJxaVVYJVV5dbVgUlXHGNMIWT6NOZdhDixzCetem90L6KtJoxgoNhizqulom2D68xQcxd80V3+fb3orRcJ4g+KtAtP5ca22mK+j8OEN9JfnLCOsKwsd922aeilMx+2dRWWbCOgFq8vbw6UDIwpfclnJeyxt4G6+FwG5mGxCzm8VTegrFCKWmOo3UEjiM/bkfblB9oVv2Z9V2GnL7ZgX7cC9trP9m7G/M/ZJYx81din9HqVktbE7zV7b+SXb/1E7+k475Vv2Z4PsCx33v53H9Pl0vuNwnNmbo+PdcYIun+EH3ZDDDwzKlUBajhy+pMPzErgqmbh0oKKr3gW1TnFTzaUN2Z33cEPMlmt5YR0WOtBc+VvYNpdAsbbpaJtg+q+NeSnl9+umt9IApItMtDpM57g3PjZkA70pb0rAHfv9VCSMNNXTAD2LekF3S6W16t5oxN2FLnkJl3/GHC+3gg8zVywB+FNNwFTnExyr6UheOlJxIG0EnZUPfem4/+2L7fqb7bxFdsILdlDu8Mn2fwjK5xg7xdgbjJ1g7LeMvZaSc4xdnof7Ax984kU78Gf2WzfbeRfb9Z0P7DXPUDjNePKk1x/uuIMOipHf1lCAjuEHBuVKIJCjJMDJjBPkFvNyQ6ZbXuagmCRkgTAta0FKaQNXhQOcUZJZubAO77lyoMPfAvMcQFFCTF9PCsPnC6m/6yOmN3gMJY9B6U+XdhxjOj97+SgoPo6DkR0fGoZjXS6ZNg7QXahN01fOhy4yEdBadUlVzFLdW0VmqiYyJiUH0OFRe6964Q9OanBmIJhH3wzixVL+1h19zqI3HZNbnt9n1PnA3svsE7fZGSvsFW/abnlA/y9jbzH2X4wdbOzZxp5hjnQjYjBl3mLswjys/912ecJeVm5vu8w+kV8mXUsHgY2m6Jp6r8tWOoKGD3Ip1kjHiLxCSMuUwwedcYJUMVRPI1R5lfTWUmUCSVU9bbycrwoH2Mr7jusqWFeYzsdOSEDnl/hNZPkx7JYAbROAesqY51J+T0VMz4jpvEDKu8n5nBCep83JA5Nh3euBYQcc4tKUeoUVi69KXSyh1qJpyV/SKOAlZEkv7a0VzvTyKWh/eWUCTOenaZr8uUU21bFYuoEiCuryJ8CU0TbRFbbLB38fbR+dYye+YAflXS7Ljf0/xv6zqetiXiFI5nCodeRGqeuSv2S/S96YneZFO3COnTjaPtrt8JtmFe0yHU7ftKs/gV16XfhIdCUHF1YwnBJievh9MW3iVFfdZFpFznfLc064Fq6q6l5NQKYqDCZpBMor+WPHiQxu9gI6jOjNxmxsOtommL7CmKdTfiua3koDTCwy0Vp8L/fR1yblea3V9Amxzclf7+F53v2lrBzQlTRlUvkoXxDdZMnst6iqJZOqJzIpi3lpWTgL7WXCmbgtsxPFYjqXZ1hfSY/q3flY9cG06b/Kdjv85v+2SxfYa7bbXrba2B8aO8zsa29W0/4FObx5tDS2r72xlxp7u7EbzRu25yI74Zt2SU/71zymT6JD1fOnu+yjhmSsC2O0EoXCFPnOV1pMz64zGLHqqptESS/hluccb2E3U1UPq4cs7LLy5qAKpoC3qjCsK0xnHzosdJh9DBQlWL1MMB0njPIyrPxbFTG9wWMoEZnyvUhMh5G+WSw9s7WuTHWF6QwZ8JN6FSgtE4rFREBxZck0bipfVnFpVVgm3cIqRxYuSKu6KhkYctolF9Nh23IVbxPshFlKlndN/gSYfrSeWWl72r9OsIuW2zF5T/qj5D0/Iw/oaQNbafImvP1G3lQ/8MEnltsxE+yiPKavoCPa+/EXanZQQ/y1aKBzWsd4ROq1A7XS5ODNV02oZNqYvPmqrpv01uJMt7DMCVSUl2QV72C90y05FKS5CXbDsvwRCAevSwDQOTSlhJheacyvUn6VEdOLwvSVdAduIvN8K71kv0J/t1JOdRIGg3B1+NMZTaBbBXVIFZBaGzbVZUnFJC0pq7h0Wi3pcHRrcU6grnspjQnnB+7VtEsKAaVti4mAWSpbn0evZVU0p3WmN0WpE6Zfbx9YYa+we4ydn7fBt6YDOg9vszH2QmPn5GNgnrCXfcfO72W3myqKTD+HP323lfY3LCJbX/aBGcixFxxOmhy8+bItl3ZnJ5DjVlc5ja4bqCgvyea8g5XTjcKSQxYa9x1OCkkz0rHSJs2+V4z5bdPRNjE6lxrzy5Tf0qa30gATi0y0Lt8LFkhdTH8l+czNJgJ9HJG6lNABx3I2zkhnhYMuQre86iuLZdRUWcWlA0zcwionUFddUhXdpHew4cyCIMitqJ7ws2qeMUvpCb3P9KDQl2W2l91+s5233l5sXzN2Wn4ttLIQpj9s8gun+XXU/zEb7Zcm2jn5bUdV9JDoxZ+o3kQNzTMGTxfZH+4hD7PgcMLSUFddCasc2Y0wrSq6yUB1t7DKCdTFJVlFDVMmZTHUzU5wdZjqCG5GLJwKnQCgs9lXWkx/iL6czpub1N8SnP5YJIw3KH68MX3NmjUXXHDBKaec0qdPn7vvvvvIkSMNumNt8hj0+l548vBNy1fExlzGdPbA8AfGVpLQvZieXY1kSamRUlMVLYvJ6mm0LO+l0yq6djqeWOATqKsuoUoaoYbJAKcyVbIgCHJbqiecnEVhMKvz+4860+YgwvSJdk5+gbQ6H7C4r70fhhW7/AG+E/Iu9RftwEn27jymV9JDIn8ewG6yABakADokzOMqOBw1fCS9skqTM/LVOAJJVGGiKDVQdd1koF1ckrUwapeQxVC3KILNKWWny0+scNgiL7axvwWAXmNMCXbtJwC1yJhlKb8SnNKuULGo5HHF9Oeff/7kk0+eMGHCmjVr7rjjjnbt2s2aNUt1d+/evSb/b6qwwBC8HMD0GvLAANMR1Hj8MV3eUVn0VSq6lw4wUeVl03wpUFdeUny8SXWLenHKWwbe57DvRfaH6Xn0PYqa/DLp0bjDJba3rZli73zJ9rdPGfv1vOPlDmOG0V+3MnLy7pd/MXal2WI/N81WfMa+lsf0sdZ0tBQFv9KY+4SyoR4IvJaVFtPdyXLFjj6ECVXR5Ryoruq6yUBdXJK1lA7IpCyGukURAUyXRjqviErfLB/MUm2M2bt3r8KcopIJpvNbAm93Un9L8IW8orqkCh9XTL/iiisuuugi9GDatGmdO3d+//33kWOt3bUrf6K8MT8Qt5nEdN5txN+x5HnaQXcmYzp/2Z3dL+xS52VS9r2wrhelQ7Kw1EipqS6NkrJ6Go3CaURaRViRqMg9QdLrTPByk1W8tDvGJmI6Yku8/eHM2fSmddRnTiHqhOnTbMU2+1n7mLGXmQ3GTCddke90Lruji6j56Jfl5jX7mdvsjDymL6OHRD7o5RVqYrZQNpdB68J0VwfCauCdbpnpCsTNkeVdVUGOLOYyyZjjntLKdrrCdGmk46St/zbG7Nq1SwJOsXSC6fcb818pv/vbij/9gw8+6NChwz333AMhvvhi/rCbp556CjlBTOePFDOms6eshtBcYTq+IdkqMF1qeRod0HVZpSl3juTjpcGcCQb0grAuDVv1+YIsmD6LztLK29lmZP57db1tzW12xmv2M/nI9H82VdkwfamhHUkPm+22V7m97Rim53cbbSVnekHUk6NG0AXeP8JxLwFZeUUtMwNTLy/JKpgpmSkLK1oWS6NVFTcpK6IDLiGLuUwy5vAjVp68jT0ryvHCXhcA+g5j8nEvJcL0Ocbcn/Kb01Yw/dVX80fF//KXvwSCv/POO8aYBx54ADnFYDo70xnQd5BjHS51uUzawu10qeVpdEDXZRXcQjIzUFdeklW8NJgzEcApWVIem662fWXB9JmE6dUK02ts7/y5Lhfml1Cz2OmLjLED8ztO37A9j2H6KA5k3ExNSFG4tGunYyzwJslRKzogK6+oZabbGW+OrILWZaa3FmfKYml0oLrLBB1wCcm/IM+0AgrTOZARUYwISIfj5WPC9LsoUEpOLui7vJi+Z8+eMWPGdOnSpXv37rfccsuhQ4ck7llrn3zyyYEDB3bq1Klfv35PPPEErs6ZM+ess87q1KnTsGHD/vjHPyI/jTh+vpff/OY3xpinn34aXTl06JAx5q677kKOtXbnzp30Pn0DuV9+QH/LjZlDK2YPkYN1Fe0H+y2teOBEtJcp+TxdWmvM41Ty50nQyzziUGFMOeHA9Eb9Lafq/LfCmMAPJQs2hJIBIsBE1kJ/ZGagrrwkq3hpMGfiLmPwU5dk8i4S+xz6yxo/p2HFgtNxnzHr8+d5Dc+ZB3Kfzr08OTdtS+6M3AKT+4JZQMrBfjo5GEXPMyb3OZP7qflDrvv03C1n535nFuXMSP6G5EbSK1VDJctprjHeOcnNzOPiEclRKxoVj97q6pJX1DJT9SQtKaugCZmZVnG6UGlZXtGB6nxJlkcHXEIWK8gzrQBPB8/CAvrE6M/pZn+cjl1cTwiwkUzy31LwIiMDA8U6Y8zOnTsl4BRLJ76X/xS6DU1g4j+9mD58+PCrr7764MGD27dvHzBgwL333iubrqmp6dix4+OPP37o0KFHHnnklFNOefPNN621S5YsOeuss7Zt21ZbWzt58uQBAwa4cSWSj7X2+GH6c889Z4x55pln0APGdOmNsdZu376dMD3+iRKIEogSKL0Etm/fDghqBJFgejnFSsk3D9BHnzr6yxt//vOfjTF/+9vfuMXKyso+ffrI1u+4446vfvWryBk1atQPf/hDa+2ll14Kq7euru6000579tlnUcxLHD9M37ZtmzHmscceQz/Y97JgwQLkWGs/+uij7du379y5c1fyb+/evbn4L0ogSiBKoGkS2L9//65duz766CMJOMXSCab/oOG7/v8VfoXvuZi+YsWKbt26oa3f/z7/fbt3330XOWPHjp08eTKSkydPHjNmjLW2a9euv/rVr5A/ZMiQn/70p0h6ieOH6bW1te3bt5dvHLxGunHjRm/PYmaUQJRAlEBLk0BtbW2vXr3Crw+dO3eura2VPV+6dOnZZ5+NnL/85S9qtfbyyy+fMWMGCpSXl19++eXW2vbt269fvx75ZWVld955J5Je4vhhurV2xIgRQ4cOhT9o2rRpXbt2VbGM3l7GzCiBKIEogRYigdraWvXCsHfv3sStkP9///79qquPPfZY9+7dkcl2uix25ZVXTpkyBQUmT548duxYa+3pp58u10uHDBkyf/58FPMSxxXT169f365du3Hjxq1evbq8vLxdu3bSbPf2L2ZGCUQJRAm0dgnU1NQYY3bv3s0Dqays7N27txzU7bffPmrUKOSMGjWqvLzcWnvJJZfMnj2b89mfLpckUV4SxxXTrbWPPfbYF77whQ4dOpx77rlz586VXbHWFjw5QJWPySiBKIEogVYhgbKysvHjx7/33ns7duwYMGDAzJkzZbdfe+21jh07PvLIIxz30rFjxz/96U/W2sWLF5911llbtmzhuJe+ffvW1dXJii59vDHd7QFyspwcgMKRiBKIEogScCWwfv3666677rzzzuvUqdOnP/3pK6+8cvPmzW6x45+ze/fucePGde/evWfPnlOmTDl8+LC19tRTT122bBl3Zu3atYMGDercufOAAQOefPJJzjxy5MjcuXOPWsCdO3ceMWIEA3248y0I07OcHBAeTLwaJRAl0MYlMG7cuBEjRixYsGDjxo2/+MUvhg4detJJJ8llxhNePi0F0zOeHHDCz0ccYBYJ3HjjjWmBBzt27MjCIZY5USWwZ88eObR//OMfZ555JseQyPwTmG4pmJ7x5IATeCbi0LJL4NVXX/3KV77Sp0+fF5J//fr1u/DCC6urq7MziSXbiARGjBjx+c9/ngfbFqyBloLpGU8OaCNaGIdZUAJ9+/a96qqruNjhw4c7der0/e9/v2CtWKCtSWD//v1dunT5+te/zgNvC9ZAS8H0jCcHtDWNjOP1SuC9995r164dNl/84Q9/MMY8/PDD3sIxsy1L4Oqrrz7ppJPkMukJbw20FEzPeHJAW9bOOHZIYNOmTcaYVatWcc7DDz9sjNmyZQsKRCJKwFpbXp4/ekWe/NoWrIGWgunx5IB4E2aXwP335z87gBORJk6c2KFDh4Jxu9n5x5IngAQqKo4eDKmPfW0L1kBLwfR4csAJcBcdtyFcd911Z5xxBporKysbMmQIkpGIEmBAr6ioUKJoC9ZAC8L0eHKA0r+YTJPA4MGDv/a1r+Fqnz592lSwGgYeCa8EfvSjHxljeG+9KtAWrIEWhOkFTw5Q0xOTbVMCdXV1HTp0uPXWWzH86667rkOHDnzmETIj0TYlMHfuXGPMqFGjkjDX+v9ZGm3BGmhZmN42tTCOOkogSqBUEhg2bJh3P5q1to1YAxHTS6VLkU+UQJRAlEDzSyBievPPQexBlECUQJRAqSQQMb1Ukox8ogSiBKIEml8CEdObfw5iD6IEogSiBEolgYjppZJk5BMlECUQJdD8EoiY3vxzEHsQJRAlECVQKglETC+VJCOfKIEogSiB5pdAxPTmn4PYgyiBKIEogVJJIGJ6qSQZ+UQJRAlECTS/BCKmN/8cxB5ECUQJRAmUSgIR00slycgnSiBKIEqg+SUQMb3556CN9OCaa67hgzg2bNiQfch8aOoZZ5xx8OBBrvXhhx9OmzZt+PDhp512GjMcNmyYl2GfPn169ux5+PDhmTNncslrrrnGW/Lee+89+vGE7t27v/POO94CMTNKoLVIIGJ6a5mpVt/PRmD6zp07TznlFGPM3LlzMf53331XHdLkxfTq6mpjzHe/+92jFQti+sGDB3v27GmMmThxIhqKRJRAa5RAxPTWOGutss+NwPRJkyYZYz75yU9K8/nAgQM333zzQw89xN83MMZ4MX3q1KnGmKeffvqosApiurV2ypQp3NbevXtbpXxjp6MESAIR06MiHCcJKEx/6623vvjFLxpj2rVrN2/ePLcTtbW1Xbp0McZcddVV7lVr7Zo1awK+l3PPPbdHjx6HDh06WjcLpr/88svM7Z577vE2FzOjBFqFBCKmt4ppOhE6KTF9z549AwYMMMa0b9/+oYce8g7v6aefZpC9//77vQUCmP673/3u6AdLb7zxRq6YBdOPHDnStWtXY0xZWZm3uZgZJdAqJBAxvVVM04nQSWD68uXL+/fvz46Oxx9/PG1sd955J2P6r3/9a2+ZAKbfeuutxpinnnqKK2bBdGvt8OHDjTEdOnQ4cuSIt8WYGSXQ8iUQMb3lz9EJ0kNg+qmnnmqMOe2005599tnA2L73ve8xpr/66qveYgFM79u3b7du3djxcrRuRkz/xje+wS2+9dZb3hZjZpRAy5dAxPSWP0cnSA+B6Yyb/fv3D0NnozGdPePf/va3IbiI6RBFJE54CURMP+GnuKUMEJg+cuRIhvVBgwa9/fbbaf1rtO/l9ttvN8asWbMGnDNiOvteTj755Oh7gegi0eokEDG91U1Za+0wMH3Dhg233HILw/r5558v4xTl2NatW8dl7rvvPplvrX2L/j3yyCNc4JJLLuEc3pd03nnnKRFJhgAAAclJREFUfepTn6qrq0MtYPrIkSPXNPz3/vvvczGskV566aWoGIkogVYngYjprW7KWmuHJaYfOXLk+uuvZ0S+4IIL3n33XXdU77///umnn+6NZeSK7t+ZM2du3brVGHPttddKhsB0t8rrr7/OJRHLePfdd8u6kY4SaF0SiJjeuuarFfdWYrq19qOPPho/fjyD7EUXXbR//353bGzOd+zYUYG+C82cM3PmzBkzZhhjnnzyScktC6bznqMOHTrs2bNH1o10lEDrkkDE9NY1X22rt3/961/5bICf/OQnGUfev3//Ll26fPjhhxnLc7F4NkBR4oqFW7IEIqa35NmJfbN8hteZZ56JM7w+DqH8+Mc/5jO8Amu2H0e7kWeUQMklEDG95CKNDKMEogSiBJpNAhHTm030seEogSiBKIGSSyBieslFGhlGCUQJRAk0mwQipjeb6GPDUQJRAlECJZdAxPSSizQyjBKIEogSaDYJRExvNtHHhqMEogSiBEougYjpJRdpZBglECUQJdBsEoiY3myijw1HCUQJRAmUXAIR00su0sgwSiBKIEqg2SQQMb3ZRB8bjhKIEogSKLkEIqaXXKSRYZRAlECUQLNJ4P8D82IYu9qzCe8AAAAASUVORK5CYII=)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "LnymAtgC-NEv"
      },
      "outputs": [],
      "source": []
    }
  ],
  "metadata": {
    "colab": {
      "provenance": []
    },
    "kernelspec": {
      "display_name": "base",
      "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.11.4"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}