{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "82df20ab-b6ea-444f-bcef-dd154163ddce",
   "metadata": {},
   "source": [
    "MPI Interior Eigensolve with Lazy, Projected Operators\n",
    "========================="
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1c7a1aea-3e53-453c-b7dc-812ce2681a2f",
   "metadata": {},
   "source": [
    "This example demonstrates some 'advanced' methods for diagonalizing large Hamiltonians.\n",
    "\n",
    "First of all, assuming we are using ``slepc4py``, we can specify the 'arch' we want to use. In this case, there is an optimized version of 'petsc' and 'slepc', compiled with float scalars, named `'arch-auto-real'`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "1154f25f-2b6c-43cd-8057-3344f54022fd",
   "metadata": {},
   "outputs": [],
   "source": [
    "# optional - comment or ignore this cell to use the default arch\n",
    "#     or specify your own arch\n",
    "import os\n",
    "os.environ['PETSC_ARCH'] = 'arch-auto-real'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8198add2-fff9-4153-a6d0-b806abe2032a",
   "metadata": {},
   "source": [
    "For real problems (like below) this generally gives a bit of a speed boost. After doing that we can import `quimb`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "3488c180-b346-4c23-beac-689e591b3f77",
   "metadata": {},
   "outputs": [],
   "source": [
    "import quimb as qu"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5b4c0a6e-017c-4317-ac00-41592d5aa6aa",
   "metadata": {},
   "source": [
    "We are not going to contsruct the Hamiltonian directly, instead leave it as a `Lazy` object, so that each MPI process can construct its own rows and avoid redudant communication and memory. To do that we need to know the size of the matrix first:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "4918a5f1-792d-41bf-b21e-bf67116695f7",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Lazy(ham_mbl, shape=(262144, 262144))>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# total hilbert space for 18 spin-1/2s\n",
    "n = 18\n",
    "d = 2**n\n",
    "shape = (d, d)\n",
    "\n",
    "# And make the lazy representation\n",
    "H_opts = {'n': n, 'dh': 3.0, 'sparse': True, 'seed': 9}\n",
    "H = qu.Lazy(qu.ham_mbl, **H_opts, shape=shape)\n",
    "H"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9d11e45c-7a47-46d4-b4c4-43375ffe82ba",
   "metadata": {},
   "source": [
    "This Hamiltonian also conserves z-spin, which we can use to make the effective problem significantly smaller. This is done by supplying a projector onto the subspace we are targeting. We also need to know its size first if we want to leave it 'unconstructed':"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "91a38a3f-3544-4905-8b89-5dfc287594e2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<Lazy(zspin_projector, shape=(262144, 48620))>"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# total Sz=0 subspace size (n choose n / 2)\n",
    "from scipy.special import comb\n",
    "\n",
    "ds = comb(n, n // 2, exact=True)\n",
    "shape = (d, ds)\n",
    "\n",
    "# And make the lazy representation\n",
    "P = qu.Lazy(qu.zspin_projector, n=n, shape=shape)\n",
    "P"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e5fa9e38-df77-404f-9e32-ca2d295b71b2",
   "metadata": {},
   "source": [
    "Now we can solve the hamiltoniain, for 5 eigenpairs centered around energy `0`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "ad9b4518-6674-4e0f-a777-3f8421925b64",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "energies: [-9.10482079e-04 -7.21481890e-04 -2.40962026e-04 -1.77843488e-04\n",
      " -5.45274570e-05]\n",
      "CPU times: user 292 ms, sys: 182 ms, total: 474 ms\n",
      "Wall time: 14.1 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "lk, vk = qu.eigh(H, P=P, k=5, sigma=0.0, backend='slepc')\n",
    "print('energies:', lk)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5b59ebf1-a43c-4b44-a93b-983713498cd8",
   "metadata": {},
   "source": [
    "`eigh` takes care of projecting `H` into the subspace ($\\tilde{H} = P^{\\dagger} H P$), and mapping the eigenvectors back the computation basis once found.\n",
    "\n",
    "Here we specified the `'slepc'` backend. In an interactive session, this will spawn the MPI workers for you (using `mpi4py`). Other options would be to run this in a script using  ``quimb-mpi-python`` which would pro-actively spawn workers from the get-go, and ``quimb-mpi-python --syncro`` which is the more traditional 'syncronised' MPI mode. These modes would be more suitable for a cluster and large problems (see `docs/examples/ex_mpi_expm_evo.py`).\n",
    "\n",
    "Now we have the 5 eigenpairs, we can compute the 'entanglement matrix' for each, with varying block size. However, seeing as we have a pool of MPI workers already, let's also reuse it to parallelize the computation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "95cdea7d-fe49-4454-81d4-653f89b3011e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[[<Future at 0x7f93354b96d8 state=running>,\n",
       "  <Future at 0x7f9302eef1d0 state=running>,\n",
       "  <Future at 0x7f9302ed8f60 state=running>,\n",
       "  <Future at 0x7f9302eef0f0 state=running>],\n",
       " [<Future at 0x7f9302eef2e8 state=pending>,\n",
       "  <Future at 0x7f9302eef240 state=pending>,\n",
       "  <Future at 0x7f9302eef3c8 state=pending>,\n",
       "  <Future at 0x7f9302eef438 state=pending>],\n",
       " [<Future at 0x7f9302eef4a8 state=pending>,\n",
       "  <Future at 0x7f9302eef518 state=pending>,\n",
       "  <Future at 0x7f9302eef588 state=pending>,\n",
       "  <Future at 0x7f9302eef5f8 state=pending>],\n",
       " [<Future at 0x7f9302eef668 state=pending>,\n",
       "  <Future at 0x7f9302eef6d8 state=pending>,\n",
       "  <Future at 0x7f9302eef748 state=pending>,\n",
       "  <Future at 0x7f9302eef7b8 state=pending>],\n",
       " [<Future at 0x7f9302eef828 state=pending>,\n",
       "  <Future at 0x7f9302eef898 state=pending>,\n",
       "  <Future at 0x7f9302eef940 state=pending>,\n",
       "  <Future at 0x7f9302eef9e8 state=pending>]]"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# get an MPI executor pool\n",
    "pool = qu.linalg.mpi_launcher.get_mpi_pool()\n",
    "\n",
    "# 'submit' the function with args to the pool\n",
    "e_k_b_ij = [[pool.submit(qu.ent_cross_matrix, vk[:, [k]], sz_blc=b)\n",
    "             for b in [1, 2, 3, 4]]\n",
    "            for k in range(5)]\n",
    "\n",
    "e_k_b_ij"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0149cd2b-2511-479b-b855-df8bf8419074",
   "metadata": {},
   "source": [
    "Once we have submitted all this work to the pool (which works in any of the modes described above), we can retrieve the results:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "c4614fb7-dde8-4f32-bad7-5e8d106d3bd9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# convert each 'Future' into its result\n",
    "e_k_b_ij = [[f.result()\n",
    "             for f in e_b_ij]\n",
    "            for e_b_ij in e_k_b_ij]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "587fbca9-159f-4bbb-af13-57a0fed823b7",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1AAAAJCCAYAAADOciYVAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvIxREBQAAIABJREFUeJzs3XucnWV97/3vb45JJiQknCERMCInrdCOKPC020ItiG7RbW2hu9XHdjfqqyioraWtPelTtT7V0j5abQTEvVXYbo9ogQgiIuU4QBRCDGAI5ACSyCHHmcms9Xv+yFQTWPfkutZ9rXvd18zn/XrxSua+13VdP9Z8s+b+zT1zLXN3AQAAAAD2rafbBQAAAABALmigAAAAACAQDRQAAAAABKKBAgAAAIBANFAAAAAAEIgGCgAAAAAC0UABAAAAQCAaKAAAAAAIRAMFAAAAAIH6ygw2s7Ml/bOkXkmXuvtHp1xs/hwfPHh+y3ONna1LGdiwvUyJSMhmz2p5fOf4Mxqf2GEVl9OW2MwO9M3x2f2tM+ujY+kLRCVGtV3jPpZFZqW43B64sNePXNz69XT1zgUtj/c8NJ6gSnTSdM5s/+CQDw4tbHmu52muAWrDCuLnXjhkq57e7O4HdaiipGKvD3qHhrxvYevcDq4nt7kKfa1tu4Eys15Jn5L0aknrJd1lZle7+wNFYwYPnq/j/+VtLc8988MDWx4/+i9ua7dEJNZz3Aktj9/+489WXEl72sns7P75OvVFf9jyXGPl6qKFiouY4gsNSujpbX282Wh5+A7/bgeLSSs2t0cu7tOt1x3Rcq4z7/+tlsdnn/VImmLRMdM5s4NDC/XS37yo5Vxzv3x7x+pEHOsfaHncdxV/A+YG/8qjnaonpXauD/oWLtSiC9/T8twL38+1a65CX2vL/AjfKZIedvc17j4u6SpJ55aYD+g0MosckVvkhswiN2QWUco0UEdIWrfHx+snj+3FzJaa2YiZjUxs2VFiOaC06MyON8gsum6fud0zs5t+1vquG1ChqMzuGuPHndB10dcHzW3kdiYr00C1+jml5/18krsvc/dhdx/umzenxHJAadGZHegls+i6feZ2z8wedEDBjzMC1YnKbP/gUEVlAYWirw965pLbmaxMA7Ve0uI9Pl4kaWO5coCOIrPIEblFbsgsckNmEaXMLnx3STrGzI6WtEHSeZJ+d6oBjZ19hZtF/PY5t7Re5BsnFU94531hlSKJ5orWv0vpPlpxJW2LzqyPjhVuFvHEe05reXzOE83C+eZdyS9Ed0TBZhHTRFRuV+9cULhZxM0v/XrL42dpitdZIF5UZnue3l64WcSDlw+3PH7k/ynerGfw2rsiSkWoqTaLmAairw8G128v3CzikQ+f2vI4G6NNH203UO4+YWYXSFqu3Vs+Xu7uK5NVBiRGZpEjcovckFnkhswiVqn3gXL3ayRdk6gWoOPILHJEbpEbMovckFnEKPM7UAAAAAAwo9BAAQAAAEAgGigAAAAACFTqd6BiDWzYXrgDSdFue6+/4qbC+a4+4YAUZXXElt99ZeG5eV9iJ7asWOvdnop229v65q2FU827MklF+9R35OLCcxOPris8h/z1PDSu2Wc90vJc0W57b3xgU+F8Xz/hoCR17UvPnOL3XGvu4A2tZ6qi3fb+9pOXFY5523f+R8vjL37nnUlqAval6Fp37d+33p1Pko76S3boywl3oAAAAAAgEA0UAAAAAASigQIAAACAQDRQAAAAABCIBgoAAAAAAtFAAQAAAECgSrcxn9Kd97U8PNVW5W9Z3Xo75mXve1PhmNnLVxSe813jhedi7f+NHxWea735NWrLveXheVe23o5+qq3Kf/Lx1tvbL3lf2q3tm5ufSjofpreptip//09avzZ/4tfOKhwzsWFjdA1sVY5WBq+9q+Xxoq3KJenGcz7R8vg79H8lqQlo11RblZ9x3/aWx2986VDhGDv5xMJzvuKBghOtr2kQhztQAAAAABCIBgoAAAAAAtFAAQAAAEAgGigAAAAACEQDBQAAAACB6rMLXxuKdtv7/rJlhWN+9YK3F56bd9ujLY83ptjRrGjnPnaUQitFu+0ddOv+hWM2nfZM9Dqbzv+lwnMHXFq8CxDwXEW77f37XdcUjnntK15XeG5i3frSNQEvfuedheeKdtv77VVPFI758vGHRtew89xTCs/N/mZxfUArRbvtLd9YvHv0WYcXz7fhz05reXzx8uJrimbRzn14Hu5AAQAAAEAgGigAAAAACEQDBQAAAACBaKAAAAAAIBANFAAAAAAEooECAAAAgEBZb2M+e3nrrR2n2qr8jR+8vvDc9ace0fK4T+yKKwyINNVW5VNtvfuRf39Dy+NL/oStypHGxIaNLY9PtVX5e2+6tvDcx5a8tHRNQDum2qr8w48Ubzv+R/9wYcvjB32G11mkYyef2PL4VFuVv21167ffkaTPHdv6eDOmKBQq1UCZ2VpJWyU1JE24+3CKooBOIrfIDZlFbsgsckNmESPFHahfd/fNCeYBqkRukRsyi9yQWeSGzCIIvwMFAAAAAIHKNlAu6TtmdreZLW31ADNbamYjZjayS2MllwOSmDK3ZBY1RGaRGzKL3HBNi2Blf4TvdHffaGYHS7rezH7s7jfv+QB3XyZpmSTNs4Vecj0ghSlzS2ZRQ2QWuSGzyA3XtAhW6g6Uu2+c/PNJSV+XdEqKooBOIrfIDZlFbsgsckNmEaPtO1BmNiSpx923Tv79NyV9MFllAXzXeMvj824r3taxaKtySer91tyWx5uv2hpXGGqrDrmNVbRVuSRd+sZlrcf8yS91qhxUrK6ZnVi3vvDcVFuVv+/hlS2Pf/zFU2S22QiuC91X18xOpWirckn67l9+vOXx8z5zWqfKQcXqkFlf8UDL4xv+rDhnRVuVS9L/9+h/tDz+riNPj6oLrZX5Eb5DJH3dzP5zni+5+3VJqgI6h9wiN2QWuSGzyA2ZRZS2Gyh3XyPpZQlrATqO3CI3ZBa5IbPIDZlFLLYxBwAAAIBANFAAAAAAEIgGCgAAAAAClX0fqFpqbH6q8JxP7Co8V7Tb3pPfPK5wzKG/t6H1XFvZuW+m6jtyceG5ZkE2N51fvAPZkj+5rfBc0W57l6y9tXDMO9/ZerepwetGCsfIebuL6axnzpzCc80dO5KuVbTb3vL1dxeOOevwk5LWgJlp57nFu1If9Jni19mi3fa+uv72wjG/teS/tDzuY7z5KgoUfJ1dvPyZwiHNKaYr2m2vaHe+qcbg+bgDBQAAAACBaKAAAAAAIBANFAAAAAAEooECAAAAgEA0UAAAAAAQiAYKAAAAAAJNy23Mfdd40vmKtiqXpFfd+kTL4ze+dChpDcjHxKProscccGnxFrrtKNqqXJLG3t16K/XBa9mqfKZKvVX51Is1Wh6eaqvyv1vTeovzD/3GmwrHTKxZG1UWpr/Z37wz6XxFW5VL0twb9mt5fOuvtreNee/+81sebzzzbFvzIR/NFQ8knW+qrcq3XLuk5fF5r/lJW2tZX+s2wycm2pqvTrgDBQAAAACBaKAAAAAAIBANFAAAAAAEooECAAAAgEA0UAAAAAAQqPa78G353VcWntv/Gz9qeTz1jlLNrVsLzxXttnfSvcXzrTi5bEXA1AavGyk+V7Db3pqPnlo45oUXp90lEIhRtNvew287rHDM0R/cWHgu9U6tmJl8rHhHvcLd9r67qHBM40MHFy920z2hZQFtK9ptb90HTisc84J/KL7ekDfLllRb3IECAAAAgEA0UAAAAAAQiAYKAAAAAALRQAEAAABAIBooAAAAAAhEAwUAAAAAgWq/jfm8L91eeK7OmyNOtVX5xA0vKDw3a2nrnrYxv/V26ZLk964MrgszhLfeqnwqU21VvnzjisJzlzx9VOsxpxxROKa5fXtwXcDEmrUtj0+1Vfn/fuT7hefOP+Gslsd/ev6JhWMO/De28kd5U21VvuZN/YXnjrmpA8UAgabaqvxrj9xSeO6Nq/9by+M9b3i2cMxUbx1UJ/u8A2Vml5vZk2Z2/x7HFprZ9Wb20OSfCzpbJhCH3CI3ZBa5IbPIDZlFKiE/wneFpLOfc+xiSd9192MkfXfyY6BOrhC5RV6uEJlFXq4QmUVerhCZRQL7bKDc/WZJTz3n8LmSPj/5989LekPiuoBSyC1yQ2aRGzKL3JBZpNLuJhKHuPvjkjT5Z+EP9ZrZUjMbMbORXRprczkgiaDcklnUCJlFbsgscsM1LaJ1fBc+d1/m7sPuPtyvwU4vB5RGZpEbMovckFnkiNziP7XbQP3UzA6TpMk/n0xXEtAx5Ba5IbPIDZlFbsgsorW7jfnVkt4q6aOTf34zWUUzQNFW5ZL0ye99oeXxM757UeGYF/9B6ZJmCnLbpqKtyiXpXfuvaXn8O/1HdqiaGYXMTsF3jReeK9qqXJLecXfrLXk/dcyWturoPWBhy+ONnz33Vy1mBDIrqXf/+a1P3HRP4Ziptip/6F9e0XrMu++IqAoFyOwk6ytoC7z4jYOKtiqXpOXHf7vl8dfOf23hmOm0jfmVkm6TdKyZrTezP9TukL3azB6S9OrJj4HaILfIDZlFbsgsckNmkco+70C5+/kFp85MXAuQDLlFbsgsckNmkRsyi1Q6vokEAAAAAEwXNFAAAAAAEIgGCgAAAAACtbsLH0pozB8qPFe0297DZy0rHHOOfrl0TcBUlp9yROG5ot32Vn302MIxR33DC88NXHdXeGFAgZ+ef2LhuaLd9o6/u/hL4qpfmSg8N0N328MUGs88m3S+ot32Gr9e/PW/93vFO/4BrfhE8etckZ43FGe9aLe91RctLhyz5E82TLFYb+vjzUbxmA7hDhQAAAAABKKBAgAAAIBANFAAAAAAEIgGCgAAAAAC0UABAAAAQCAaKAAAAAAIxDbmXeD3riw89+I/aH18qq3K137o1MJzh/9H6y0p2SoaMZrbt0ePmWqr8tGFBVuRShqIXgl4vgP/7bboMVNtVb55afHr7MF3tt7Gt7nigegagBhTbVX+xEWnFZ477FMj0Wv5rvGWx/sWLyoe9Fj0MshMc+vW6HNTbVW+7q+Kc9tfsNShl9xaOKZTuAMFAAAAAIFooAAAAAAgEA0UAAAAAASigQIAAACAQDRQAAAAABCIBgoAAAAAArGN+TRQtFW5JJ3ykdbblX/5v728cMxx776v5fHm6GhcYZjRptoqf6qtyh/5aOvtoo++uHhbaj/1ZS2P220/nGIlzFS9Byxsebzxs6cKxxRtVS5JSz77k5bHHyp+mQU6bqqtyq979M6Wx7+6bV7hmMtOf0XL4xPr1scVBvQUv5VJ0VblkrR1SaPl8UPL1tMG7kABAAAAQCAaKAAAAAAIRAMFAAAAAIFooAAAAAAgEA0UAAAAAATa5y58Zna5pNdJetLdXzJ57G8l/ZGkTZMP+wt3v6ZTRWJqU+12VrTb3ifP+F+FY971/7615fFj3nVHXGFdQmbzVrTb3utWPl045rZnHmt5/KkzBlsPGLPoujqN3FZnqt32ijRXPFB4rmi3vbesXlc45sOf/52Wxxd9pGC3SS+cqmvIbL6Kdtt709wthWMumyje8TcXZLYmmq1305OkQy+5tfhcwfHTfjheOOaWzUtaHu896/HWAwJjHnIH6gpJZ7c4/k/uftLkfwQNdXKFyCzyc4XILfJyhcgs8nKFyCwS2GcD5e43S4r/dh3QJWQWOSK3yA2ZRW7ILFIp8ztQF5jZj8zscjNbkKwioHPILHJEbpEbMovckFlEabeB+rSkJZJOkvS4pI8XPdDMlprZiJmN7NJYm8sBpZFZ5Cgot2QWNUJmkRuuDxCtrQbK3X/q7g13b0r6rKRTpnjsMncfdvfhfhX8QjfQYWQWOQrNLZlFXZBZ5IbrA7SjrQbKzA7b48M3Sro/TTlAZ5BZ5IjcIjdkFrkhs2hHyDbmV0p6laQDzWy9pL+R9CozO0m7N1ZdK+ntHawRJRz37vtaHi/aqlySjjlhQ6fKqUTyzPb0tj4+xTacaJ+f+rKWx4u2Kpekdxz2vZbH//a//I/Wa9x2U3RdncZr7fRTtFW5JB1+ZsEW5x+u4X7lBZJl1kzWP9DylO8q3p4YU5vqubvs9Fe0Pj7FVuUbPtd6E+n+bx1XXMSlXyk+1wW8zk5PRVuVS9Ifv6D19cE/nf27LY/7zT8IWnOfDZS7n9/i8GVBswNdQGaRI3KL3JBZ5IbMIpUyu/ABAAAAwIxCAwUAAAAAgWigAAAAACAQDRQAAAAABNrnJhIp2exZ6jnuhJbnmiseqLKUGaM5Otry+DHvuiN6rmf/+ytbHm9cc3v0XFlht71K2W0/bHn8qTOK33OjaLe9jae33tlr130WXxgQadFHbis+WbDb3vhZwy2P+61TzJU7d3bb64C+xYsKz02sWx89X9Fue0+/NJ+dIzE99Z71eOG5ot32Hnt9wWvwj8LW5A4UAAAAAASigQIAAACAQDRQAAAAABCIBgoAAAAAAtFAAQAAAEAgGigAAAAACGTu1W0/aWabJD06+eGBkjZXtnh9TYfn4Uh3P6jbRXQCmW1pOjwPMyWz0vT4fJU1HZ4DMjuzTJfnYKbkdrp8vsqaDs9DUGYrbaD2WthsxN1bv+HFDMLzkA8+V7vxPOSFzxfPQW74fPEc5IbP124z6XngR/gAAAAAIBANFAAAAAAE6mYDtayLa9cJz0M++FztxvOQFz5fPAe54fPFc5AbPl+7zZjnoWu/AwUAAAAAueFH+AAAAAAgEA0UAAAAAATqSgNlZmeb2Woze9jMLu5GDd1gZpeb2ZNmdv8exxaa2fVm9tDknwu6WSNaI7NkNjdklszmhsyS2dyQ2Zmb2cobKDPrlfQpSa+RdIKk883shKrr6JIrJJ39nGMXS/quux8j6buTH6NGyCyZzQ2ZJbO5IbNkNjdkdmZntht3oE6R9LC7r3H3cUlXSTq3C3VUzt1vlvTUcw6fK+nzk3//vKQ3VFoUQpDZvZHZ+iOzeyOz9Udm90Zm64/M7m1GZbYbDdQRktbt8fH6yWMz1SHu/rgkTf55cJfrwfOR2b2R2fojs3sjs/VHZvdGZuuPzO5tRmW2Gw2UtTjGXuqoMzKL3JBZ5IbMIjdkdgbrRgO1XtLiPT5eJGljF+qoi5+a2WGSNPnnk12uB89HZvdGZuuPzO6NzNYfmd0bma0/Mru3GZXZbjRQd0k6xsyONrMBSedJuroLddTF1ZLeOvn3t0r6ZhdrQWtkdm9ktv7I7N7IbP2R2b2R2fojs3ubUZk19+rvNprZOZIukdQr6XJ3//vKi+gCM7tS0qskHSjpp5L+RtI3JH1Z0gskPSbpze7+3F/MQ5eRWTKbGzJLZnNDZslsbsjszM1sVxooAAAAAMhRV95IFwAAAAByRAMFAAAAAIFooAAAAAAgEA0UAAAAAASigQIAAACAQDRQAAAAABCIBgoAAAAAAvVVuVjv3CHvO2BB1JhZm5rxC42Nx4+R5M021pqGbHAg6vE7d23ReGOHdaicruofHPLBOXGZ7RlvI0c7RuPHoG2j2q5xH5uWmR2YP9tnHTovaoyt641fqM33EPSxsTYGtbVUrVlP3Pcvdza3adxHp2Vm++YMef+8hVFj+re18Tq7fWf8GElq51mfhpmVxT8RW/2pze5+UAeq6breoSHvXxCX24EN26PXsTaed0nifV7bE3p9UGkD1XfAAh36lxdGjTnu01vjF3r4sfgxkprb44M9HfUe+cKox9/26Oc7VEn3Dc5ZoJN+PS6zc9bF58jvXhk9BnuI/AJzR/OGDhXSfbMOnaeXf+a/R43pe+/c6HVsdFf0GElq/mRt9BifmGhrrTrrmbtf1ONv33Z1hyrpvv55C/XCt743aszh32/j2uDO++LHSLK++Eul2me2jYtyG4j75qokXT/6xUejB2Wif8FCLbrgPVFjjv7AndHr9Az0R4+RpGZb36yahk1Xh64P+BE+AAAAAAhUqoEys7PNbLWZPWxmF6cqCugUMosckVvkhswiN2QWMdpuoMysV9KnJL1G0gmSzjezE1IVBqRGZpEjcovckFnkhswiVpk7UKdIetjd17j7uKSrJJ2bpiygI8gsckRukRsyi9yQWUQp00AdIWndHh+vnzwG1BWZRY7ILXJDZpEbMosoZRqoVttaPG/7DjNbamYjZjbS2MYud+iq6MxOjJFZdN0+c7tnZsefbXOrZiCdqMw2dvI6i66Lv6Zl5+YZrUwDtV7S4j0+XiRp43Mf5O7L3H3Y3Yd75w6VWA4oLTqzfYNkFl23z9zumdmB+bMrLQ5oISqzvbN5nUXXxV/TDpHbmaxMA3WXpGPM7GgzG5B0nqTp+0YVmA7ILHJEbpEbMovckFlEafuNdN19wswukLRcUq+ky92ddwNFbZFZ5IjcIjdkFrkhs4jVdgMlSe5+jaRrEtUCdByZRY7ILXJDZpEbMosYpRqoWLM2NXXcp7dGjXnjVd+PXufK97w2eowkDVx3V/SYsde8PHrMrOvvjR7jExPRY9rVeGhN1ON37/g5PfWMNzVnXdwvij747sHodY67aEH0GEna9ZKjosf0/CA+f7Xnz/td3xnL1vWq771zo8a86gsj0et8/3dOjh4jSXbiMdFj/Ier2lqrzppb474Wujc7VEn39W9r6vDvxz0f+3388eh1dr75kOgxktQ47MD4QffW/OZFG6+ZPjbWgULyNbBhu47+wJ1RYzZ85bjodY74aHu/bdP3+NPRYybWrW9rrVrr0PVBmd+BAgAAAIAZhQYKAAAAAALRQAEAAABAIBooAAAAAAhEAwUAAAAAgWigAAAAACAQDRQAAAAABKKBAgAAAIBANFAAAAAAEIgGCgAAAAAC0UABAAAAQCAaKAAAAAAI1FfpamPj0sOPRQ258j2vjV7m+A/dFz1GktZsOjF6zODye6LHeLMRPQZdsmNUfvfKqCHHXbQgepkXXb8teowkXf/IRPSYI2+L/2fvE/HroEvcZaO7ooZ8/3dOjl7mHd/8dvQYSfrIX78lesy8H7a1FHKxfad0Z9zX7Z1vPiR6mVfd8HD0GEn69I1HRY958X28zk53Zqaegf6oMUd8NP6+xfuv/FL0GEl696Vvjx6z6KMbosf0zJ0bPaa5dWv0mLrhDhQAAAAABKKBAgAAAIBANFAAAAAAEKjtBsrMFpvZ98xslZmtNLMLUxYGpEZmkSNyi9yQWeSGzCJWmU0kJiS9z93vMbP9JN1tZte7+wOJagNSI7PIEblFbsgsckNmEaXtO1Du/ri73zP5962SVkk6IlVhQGpkFjkit8gNmUVuyCxiJfkdKDM7StLJku5IMR/QaWQWOSK3yA2ZRW7ILEKUfh8oM5sr6auSLnL3LS3OL5W0VJJm2VDZ5YDSojKrORVXB7Q2VW73ymzfvC5UBzxfcGZ5nUVNcE2LUKXuQJlZv3YH7Yvu/rVWj3H3Ze4+7O7DAzarzHJAabGZ7ddgtQUCLewrt3u9zvZxMYrui8ksr7Oog+hrWnI7o5XZhc8kXSZplbt/Il1JQGeQWeSI3CI3ZBa5IbOIVeYO1OmSfl/SGWa2YvK/cxLVBXQCmUWOyC1yQ2aRGzKLKG3/DpS73yLJEtYCdBSZRY7ILXJDZpEbMotYSXbhAwAAAICZoPQufDG82VRz+/aoMQPX3RW9zppNJ0aPkaS1fxbfTy7508OjxzQ2PB49xicmosegO3a95KjoMdc/0t7n93O/ckX0mL9+xR9Gj+m998HoMc0dO6LHoDwfG1PzJ2ujxtiJx0Sv85G/fkv0GEl6+9+0/N3sKX35hpdFj2ls2hQ9Bl1ikvXFXY40DjsweplP33hU9BhJuvON8b8S85Z/+4PoMY2Vq6PHoHvcXc2xsagxfY8/Hb3Ouy99e/QYSVr2R5+MHvP//NNp8Qs1GvFjpgHuQAEAAABAIBooAAAAAAhEAwUAAAAAgWigAAAAACAQDRQAAAAABKKBAgAAAIBANFAAAAAAEIgGCgAAAAAC0UABAAAAQCAaKAAAAAAIRAMFAAAAAIFooAAAAAAgkLl7ZYvNs4X+CjszaszYa14evc7g8nuix0hS3+LDo8cc97UN0WPu/5Vm9Jg6u8O/qy3+lHW7jk5oJ7PtsL6+tsY1X/GS6DE/eUf8p+rYDz4bPabx0JroMVUhs93Te9BB0WPWfvqQ6DFH/+nW6DETj66PHiNJajbaGxeBzJbX7utsz7FLoseM/fPO6DF9H9g/eoxu/1H8mArd4F+5292Hu11HJ1T2Wmvt/bPvGRyMHvPjT8VfUxz/ifjX2sbK1dFjqhL6WssdKAAAAAAIRAMFAAAAAIFooAAAAAAgUOkGysx6zexeM/t2ioKATiOzyA2ZRY7ILXJDZhEqxR2oCyWtSjAPUBUyi9yQWeSI3CI3ZBZBSjVQZrZI0mslXZqmHKCzyCxyQ2aRI3KL3JBZxCh7B+oSSe+XVLgvt5ktNbMRMxvZpbGSywGlkVnkhswiR1PmlsyihnitRbC2Gygze52kJ9397qke5+7L3H3Y3Yf7Fb8nPZAKmUVuyCxyFJJbMos64bUWscrcgTpd0uvNbK2kqySdYWZfSFIV0BlkFrkhs8gRuUVuyCyitN1Aufufu/sidz9K0nmSbnT330tWGZAYmUVuyCxyRG6RGzKLWLwPFAAAAAAE6ksxibvfJOmmFHMBVSCzyA2ZRY7ILXJDZhEiSQPVSbOuvzd6jDcbba3V2PB49Jj7f6Vws5ZCW65dEj1m9N8PiR4jSQd/8ta2xqFaPjHR1rjeex+MHnPsBw+NHrPq4oXRY150xfzoMZLUd0/8/1Nz+/a21kL1Gps2RY85+k/nRI959bd/GD3mq391VvQYSRp6ZFv0GL93ZVtroX3tvs42Vq6OHtP3gV+KHvOTN8fnfMnt0UOQmZ65c9sb2Ii/Fj7+E1ujx2z4+/gfZjvs706MHiNJPY+sjx7TeObZttbaF36EDwAAAAAC0UABAAAAQCAaKAAAAAAIRAMFAAAAAIFooAAAAAAgEA0UAAAAAASigQIAAACAQDRQAAAAABCIBgoAAAAAAtFAAQAAAEAgGigAAAAACEQDBQAAAACB+rpdwL74xMS0W2uon138AAAgAElEQVT03w+JHnPgG9e1tVbzhydHj+m5ZUXcAI9eAok0d+yIH/TQmughL7pifvSYLUfPih4jSQtX8H0d7G3i0fXRY776V2dFj1l/TiN6jCSd8KEt0WOq+8qGrrj9R9FDltwev8yW818ZP0jSwjueiB/kbXyxj/9yg+dobt1a3WIrV0cPOezvTowe88hvzYseI0mLb1gSPab3pnvjBgTGnCsVAAAAAAhEAwUAAAAAgWigAAAAACBQqQbKzPY3s6+Y2Y/NbJWZnZqqMKATyCxyRG6RGzKL3JBZxCi7icQ/S7rO3X/LzAYkzUlQE9BJZBY5IrfIDZlFbsgsgrXdQJnZPEm/Jun/liR3H5c0nqYsID0yixyRW+SGzCI3ZBaxyvwI3wslbZL0OTO718wuNbOh5z7IzJaa2YiZjezSWInlgNLILHK0z9ySWdQMmUVuuD5AlDINVJ+kX5b0aXc/WdJ2SRc/90Huvszdh919uF+DJZYDSiOzyNE+c0tmUTNkFrnh+gBRyjRQ6yWtd/c7Jj/+inaHD6grMosckVvkhswiN2QWUdpuoNz9CUnrzOzYyUNnSnogSVVAB5BZ5IjcIjdkFrkhs4hVdhe+d0n64uRuJWskva18SUBHkVnkiNwiN2QWuSGzCFaqgXL3FZKGE9UCdByZRY7ILXJDZpEbMosYpd5IFwAAAABmkrI/woc2HPzJW6PHNH94cltrrXl7/JgFx70y6vETX7stfhFkpe+eB6PHLFzR3vdnHvzXF0WPGVw9O+rx45feHr0GuqjZiB4y9Mi26DEnfGhL9BhJGr56TfSY5R/71ajHN66Z5pk1i3u8e2fqyMzCO55oa9yP33Vo9JjDf9DGcx7/TwOZ6XlkffSYxTcsaWutn120I3pM45dPjXr8rv8Vdk3LHSgAAAAACEQDBQAAAACBaKAAAAAAIBANFAAAAAAEooECAAAAgEA0UAAAAAAQiAYKAAAAAALRQAEAAABAIBooAAAAAAhEAwUAAAAAgWigAAAAACAQDRQAAAAABOqrcjEbHFDvkS+MGtN4aE2HqslLzy0r2hq34LhXRo955sVxj2/Mil4iL2Zxj3fvTB1d1Ny+vbK1BlfPjh6za7+459yn8beOrKdHPXP3ixrT3Lq1Q9V0j9+7MnrMRJtrLf/Yr0aPefZFcSFsDEYvkQ8z2cBA1BAfG+tQMZlp8+vN4T+IH/fTl7fxwvnV+CFZ4fpAjWeejR7Te9O97a31y6dGj9l+8s6oxze/2gx63DS+jAAAAACAtGigAAAAACBQqQbKzN5jZivN7H4zu9LMpvsPcyFzZBY5IrfIDZlFbsgsYrTdQJnZEZLeLWnY3V8iqVfSeakKA1Ijs8gRuUVuyCxyQ2YRq+yP8PVJmm1mfZLmSNpYviSgo8gsckRukRsyi9yQWQRru4Fy9w2S/lHSY5Iel/Ssu38nVWFAamQWOSK3yA2ZRW7ILGKV+RG+BZLOlXS0pMMlDZnZ77V43FIzGzGzkfFG3FaCQErtZHaX2CoX3RWS271eZ320G2UCPxeb2V1kFl3G9QFilfkRvt+Q9Ii7b3L3XZK+Jum05z7I3Ze5+7C7Dw/0xr+/C5BQdGb7NZ3ffAWZ2Gdu93qd5fee0X1Rme0ns+g+rg8QpUwD9ZikV5rZHDMzSWdKWpWmLKAjyCxyRG6RGzKL3JBZRCnzO1B3SPqKpHsk3Tc517JEdQHJkVnkiNwiN2QWuSGziNVXZrC7/42kv0lUC9BxZBY5IrfIDZlFbsgsYpTdxhwAAAAAZgxz9+oWM9sk6dGC0wdK2lxZMfWU63NwpLsf1O0iOoHM7lOuzwGZndlyfB7I7MyW6/MwE3Ob6+cqtVyfh6DMVtpATcXMRtx9uNt1dBPPQV74fPEc5IbP1248D/ngc7Ubz0M++FztNt2fB36EDwAAAAAC0UABAAAAQKA6NVBsF8lzkBs+XzwHueHztRvPQz74XO3G85APPle7TevnoTa/AwUAAAAAdVenO1AAAAAAUGs0UAAAAAAQqOsNlJmdbWarzexhM7u42/V0i5mtNbP7zGyFmY10ux4UI7O7kdl8kNndyGxeyC2ZzQ2ZnTmZ7ervQJlZr6QHJb1a0npJd0k6390f6FpRXWJmayUNu3uObzo2Y5DZXyCzeSCzv0Bm80FudyOz+SCzu82UzHb7DtQpkh529zXuPi7pKknndrkmYCpkFrkhs8gRuUVuyOwM0u0G6ghJ6/b4eP3ksZnIJX3HzO42s6XdLgaFyOwvkNk8kNlfILP5ILe7kdl8kNndZkRm+7q8vrU4NlP3VT/d3Tea2cGSrjezH7v7zd0uCs9DZn+BzOaBzP4Cmc0Hud2NzOaDzO42IzLb7TtQ6yUt3uPjRZI2dqmWrnL3jZN/Pinp69p9Kxj1Q2YnkdlskNlJZDYr5FZkNjNkVjMns91uoO6SdIyZHW1mA5LOk3R1l2uqnJkNmdl+//l3Sb8p6f7uVoUCZFZkNjNkVmQ2QzM+t2Q2O2R2BmW2qz/C5+4TZnaBpOWSeiVd7u4ru1lTlxwi6etmJu3+nHzJ3a/rbklohcz+HJnNBJn9OTKbEXIricxmhcxKmkGZ7eo25gAAAACQk27/CB8AAAAAZIMGCgAAAAAC0UABAAAAQCAaKAAAAAAIRAMFAAAAAIFooAAAAAAgEA0UAAAAAASq9I10e/cb8r4DFySbb/DJZrK5JEk7x5JNZYMDyeaSJB9NV5sk2UC6+nZOPKvxxk5LNmGN9A8O+eDchcnm631mZ7K5dk/Ym2wq37Ur2VydMPnGfEns9O0a99Fpmdm++XN88OD5yeazzekyJkk94wlft3eMppurA6wn3fcodza3TdvM9g4Ned/CdK+zA1vTvr9lz7Z0X3+90Ug2VydYb9rvq29p/Gyzux+UdNKa6B0a8v4F6XI7+OR4srkkSQmz5s3E19uJdeP6oNIGqu/ABTrsg3+cbL4XX5K2qdD9DyWbqudFRyebS5IaDzyYdL6+RS9INtet67+QbK66GZy7UC8566Jk883/1n3J5pKknv3TXShPbNiYbK5O6Jk1K9lct49ek2yuuhk8eL6O++c/SDbf7M/un2wuSZqzYUeyuXzk/mRzdULPnKFkc92+49vJ5qqbvoULdcR70r3OLroxbZMy59Z0X38bzzybbK5O6J07L+l8y5+9/NGkE9ZI/4KFWnTBe5LNd8y/PpZsLklqPv1Murl2pHvdliR54m9ydOH6gB/hAwAAAIBANFAAAAAAEIgGCgAAAAAC0UABAAAAQKBSDZSZnW1mq83sYTO7OFVRQKeQWeSI3CI3ZBa5IbOI0XYDZWa9kj4l6TWSTpB0vpmdkKowIDUyixyRW+SGzCI3ZBaxytyBOkXSw+6+xt3HJV0l6dw0ZQEdQWaRI3KL3JBZ5IbMIkqZBuoISev2+Hj95DGgrsgsckRukRsyi9yQWUQp00C1epfe570zlpktNbMRMxtpbN1eYjmgtOjM7hols+i6feZ2z8xOPJv4DQ+BeFGZbWzndRZdF39NS25ntDIN1HpJi/f4eJGkjc99kLsvc/dhdx/u3S/du7IDbYjObP8sMouu22du98xs3/w5lRYHtBCV2d4hXmfRdfHXtOR2RivTQN0l6RgzO9rMBiSdJ+nqNGUBHUFmkSNyi9yQWeSGzCJKX7sD3X3CzC6QtFxSr6TL3X1lssqAxMgsckRukRsyi9yQWcRqu4GSJHe/RtI1iWoBOo7MIkfkFrkhs8gNmUWMUm+kCwAAAAAzCQ0UAAAAAASigQIAAACAQDRQAAAAABCo1CYSsQafbOrFl4wlm2/sY9uSzSVJvR9+SbK5BtY9nWwuSZI/7/3cSpl4dH2yubwxnmyuuul9Zqfmf+u+ZPOt/8KRyeaSpNGH5ieba8lfPJlsLknyiYmk8zVHR5PN5Yn/PdWJbe7V7M/un2y+U//2zmRzSdINnz412VwHjiSbqiOaCd9o072ZbK66GdjqWnRjI9l8b/7H65LNJUkfv/GcZHMd+2f3J5tLSpsxSWps2ZJ0vuls8MlxHfOvjyWbb8GX034u77z5l5LN9aJ/fDDZXJLU2PyzpPM1x9L1FqHXB9yBAgAAAIBANFAAAAAAEIgGCgAAAAAC0UABAAAAQCAaKAAAAAAIRAMFAAAAAIFooAAAAAAgEA0UAAAAAASigQIAAACAQDRQAAAAABCIBgoAAAAAAtFAAQAAAEAgGigAAAAACEQDBQAAAACBaKAAAAAAIBANFAAAAAAEooECAAAAgEA0UAAAAAAQiAYKAAAAAAL1VbrazjHp/oeSTdf74Zckm0uSnn7vtmRzHfrO8WRzSVLvvHlJ52ts2ZJ0vmmrt1c9+89PNt3oQ+nmkqRfOX11srl+euZJyeaSpIHlI0nnQ5ie8abmbNiRbL4bPn1qsrkk6Y0XfC/ZXD/40SnJ5pIk3XFf2vnc0843TfVsG9OcWx9MNt/Hbzwn2VyS9O5fX55srqtPPTPZXJLUf+OKpPOp2Ug733TWaKj59DPJprvz5l9KNpcknfPqu5LNteobJySbS5K0+Wdp5+vCay13oAAAAAAgEA0UAAAAAASigQIAAACAQDRQAAAAABCIBgoAAAAAAtFAAQAAAECgthsoM1tsZt8zs1VmttLMLkxZGJAamUWOyC1yQ2aRGzKLWGXeB2pC0vvc/R4z20/S3WZ2vbs/kKg2IDUyixyRW+SGzCI3ZBZR2r4D5e6Pu/s9k3/fKmmVpCNSFQakRmaRI3KL3JBZ5IbMIlaZO1A/Z2ZHSTpZ0h0tzi2VtFSSZmlOiuWA0oIz27tfpXUBUynK7V6ZHZhfeV1AkaDM9gxVXhdQJPj6wMjtTFZ6Ewkzmyvpq5Iucvctzz3v7svcfdjdh/ttVtnlgNJiMjvQM7v6AoEWpsrtXq+zfXyjCvUQmtkB43UW9RB1fcA17YxWqoEys37tDtoX3f1raUoCOofMIkfkFrkhs8gNmUWMMrvwmaTLJK1y90+kKwnoDDKLHJFb5IbMIjdkFrHK3IE6XdLvSzrDzFZM/ndOorqATiCzyBG5RW7ILHJDZhGl7U0k3P0WSZawFqCjyCxyRG6RGzKL3JBZxCq9iQQAAAAAzBQ0UAAAAAAQiAYKAAAAAAIleSPdUDY4oJ4XHZ1svoF1TyebS5IOfed4srm2XjaQbC5J8s8cn3S+OV9/3vvDoQXftUsTGzYmm2/JXzyZbC5J+umZJyWb69HX9yabS5KOefqlSefrWbU22Vy2bRp/72jHqHzk/mTTHTiSbCpJ0g9+dEqyuZ76wGiyuSRp4H++Iul887+zKtlctiXtv8868UZDjWeeTTbfsX+WLv+SdPWpZyab67HX9CebS5KW3NBIOh/CebOp5o4dyeZ70T8+mGwuSVr1jROSzfXgW9O+59VxO45LOp8eXJtsKhsN+1W4aXwVAQAAAABp0UABAAAAQCAaKAAAAAAIRAMFAAAAAIFooAAAAAAgEA0UAAAAAASigQIAAACAQDRQAAAAABCIBgoAAAAAAtFAAQAAAEAgGigAAAAACEQDBQAAAACBaKAAAAAAIBANFAAAAAAEooECAAAAgEA0UAAAAAAQiAYKAAAAAALRQAEAAABAIBooAAAAAAjUV+ViPjqmxgMPJpzQ080lqXfevGRz+WeOTzaXJD3zlq1J5/PeVySbq3nD7cnmmu58YiLpfAPLR5LNdczTL002lyQ9+fL9ks536NiidJOtHkg3F+LccV+yqQb+Z7rXMUnaeHbaf5/73zk/3WTbe9PNNc01t29POl//jSuSzbXkhkayuSTpJ//4yqTzzdqc+PvqH/lK2vnqJuF1aGPzz5LNJUlKON9xO45LNpckPfFrC5POd9hYwtfutWHXB9yBAgAAAIBANFAAAAAAEIgGCgAAAAAC0UABAAAAQCAaKAAAAAAIRAMFAAAAAIFKN1Bm1mtm95rZt1MUBHQamUVuyCxyRG6RGzKLUCnuQF0oaVWCeYCqkFnkhswiR+QWuSGzCFKqgTKzRZJeK+nSNOUAnUVmkRsyixyRW+SGzCJG2TtQl0h6v6Rm0QPMbKmZjZjZyC6NlVwOKI3MIjdkFjmaMrdkFjXEay2Ctd1AmdnrJD3p7ndP9Th3X+buw+4+3K/BdpcDSiOzyA2ZRY5CcktmUSe81iJWmTtQp0t6vZmtlXSVpDPM7AtJqgI6g8wiN2QWOSK3yA2ZRZS2Gyh3/3N3X+TuR0k6T9KN7v57ySoDEiOzyA2ZRY7ILXJDZhGL94ECAAAAgEB9KSZx95sk3ZRiLqAKZBa5IbPIEblFbsgsQnAHCgAAAAAC0UABAAAAQCAaKAAAAAAIRAMFAAAAAIGSbCIRygYG1LfoBcnmm3h0fbK5JKmxZUuyueZ8/Y5kc0mS974i6Xw/O7E32VwTtySbqnbMTD2zZiWbrzk6mmyu1HpWrU0636Fji5LOt+XY+cnmajyaLv91Yz096pkzlGy+5vbtyeaSJLknm2r+d1Ylm0uS9r8zXcYk6cF3HJFsrtF/6U82V91Yb496585LNl/Kr+WSpGYj7XwJzdqc9vvgOxbX9/+1bpJfH4yNJZtLUtLXWj24Nt1ckg4bm0g638bfPDjZXLuuCmuNuAMFAAAAAIFooAAAAAAgEA0UAAAAAASigQIAAACAQDRQAAAAABCIBgoAAAAAAtFAAQAAAEAgGigAAAAACEQDBQAAAACBaKAAAAAAIBANFAAAAAAEooECAAAAgEA0UAAAAAAQiAYKAAAAAALRQAEAAABAIBooAAAAAAhEAwUAAAAAgWigAAAAACAQDRQAAAAABDJ3r24xs02SHg146IGSNne4nHbVuTapO/Ud6e4HVbxmJchsJchsQtMks1K96yOzCZHZSnSrtpme2zpnQqp3fbXObKUNVCgzG3H34W7X0Uqda5PqX990Vefnvc61SfWvb7qq+/Ne5/rqXNt0Vvfnvc711bm26azuz3ud66tzbRI/wgcAAAAAwWigAAAAACBQXRuoZd0uYAp1rk2qf33TVZ2f9zrXJtW/vumq7s97neurc23TWd2f9zrXV+faprO6P+91rq/OtdXzd6AAAAAAoI7qegcKAAAAAGqnVg2UmZ1tZqvN7GEzu7jb9ezJzBab2ffMbJWZrTSzC7td03OZWa+Z3Wtm3+52LTMFmS2HzHZHXXNLZlGEzLaPzHYHmS2n7rmtTQNlZr2SPiXpNZJOkHS+mZ3Q3ar2MiHpfe5+vKRXSvrjmtUnSRdKWtXtImYKMpsEma1YzXNLZvE8ZLY0MlsxMptErXNbmwZK0imSHnb3Ne4+LukqSed2uaafc/fH3f2eyb9v1e5P6hHdreoXzGyRpNdKurTbtcwgZLYEMts1tc0tmUUBMtsmMts1ZLaEHHJbpwbqCEnr9vh4vWr2Cf1PZnaUpJMl3dHdSvZyiaT3S2p2u5AZhMyWQ2a7I4vcklnsgcy2j8x2B5ktp/a5rVMDZS2O1W6LQDObK+mrki5y9y3drkeSzOx1kp5097u7XcsMQ2bbRGa7qva5JbN4DjLbBjLbVWS2Tbnktk4N1HpJi/f4eJGkjV2qpSUz69fusH3R3b/W7Xr2cLqk15vZWu2+TXyGmX2huyXNCGS2fWS2e2qdWzKLFshse8hs95DZ9mWR29q8D5SZ9Ul6UNKZkjZIukvS77r7yq4WNsnMTNLnJT3l7hd1u54iZvYqSX/i7q/rdi3THZlNg8xWq865JbNohcyWR2arRWbTqHNua3MHyt0nJF0gabl2/0Lbl+sQtD2cLun3tbsTXjH53zndLgrdQ2aRo5rnlsziecgsckNmp7/a3IECAAAAgLqrzR0oAAAAAKg7GigAAAAACEQDBQAAAACBaKAAAAAAIBANFAAAAAAEooECAAAAgEA0UAAAAAAQiAYKAAAAAAL1VblY77wh7z9o/yqXlCT1b7LK17Rml96geMfOypcc9e0a97Hqn+QK9M0e8oH9Fla+bu949fmxLTsqX1OS1IV/KqOavpnt1uusRqv/ftzA493KbPWhndaZHRry/gVdeJ0drXxJ9T61vfpFpa68zkrSVj292d0P6s7qndU7d8j7DlhQ+br9W6t/Geh5htfa56q0geo/aH8d+bG3V7mkJOmwfxuofM3+bbsqX1OSbOSByte8fWJ55WtWZWC/hTr2Te+pfN25GxuVrzn72nsqX1OSfGKi8jXv8O9WvmZV+g/aXy/4h+pfZ3313MrXPPqDd1e+piT5rvHK15zWmV2wUIveVf3r7IJVlS+pBV+6q/pFJXmXvql7Q+N/P9qVhSvQd8ACHfrnF1a+7uE3Vf/Nqv2+eW/la0qSj41Vvmboay0/wgcAAAAAgWigAAAAACAQDRQAAAAABKKBAgAAAIBANFAAAAAAEIgGCgAAAAAC0UABAAAAQCAaKAAAAAAIRAMFAAAAAIFKNVBmdraZrTazh83s4lRFAZ1CZpEjcovckFnkhswiRtsNlJn1SvqUpNdIOkHS+WZ2QqrCgNTILHJEbpEbMovckFnEKnMH6hRJD7v7Gncfl3SVpHPTlAV0BJlFjsgtckNmkRsyiyhlGqgjJK3b4+P1k8f2YmZLzWzEzEYaW7aXWA4oLTqzEzvJLLpun7nldRY1E5fZ7WQWXRd/TbuN3M5kZRooa3HMn3fAfZm7D7v7cO+8oRLLAaVFZ7ZvNplF1+0zt7zOombiMjtEZtF18de0c8ntTFamgVovafEeHy+StLFcOUBHkVnkiNwiN2QWuSGziFKmgbpL0jFmdrSZDUg6T9LVacoCOoLMIkfkFrkhs8gNmUWUvnYHuvuEmV0gabmkXkmXu/vKZJUBiZFZ5IjcIjdkFrkhs4jVdgMlSe5+jaRrEtUCdByZRY7ILXJDZpEbMosYpd5IFwAAAABmEhooAAAAAAhEAwUAAAAAgWigAAAAACAQDRQAAAAABKKBAgAAAIBANFAAAAAAEIgGCgAAAAAC0UABAAAAQKC+Khfr32Q67N8GqlxSkvTY25qVr3nkZf2VrylJfRMT1S/q1S9Zld5x19yNjcrXXfcb1X9vY+EhL698TUk64LLbq190GmdWoz3y1XMrX/ao09ZVvqaffGzla0qSRh6ofs3qX4Yq0zsqLVhV/brN3/5Z5Wtu7u/S6+zn7uzKutNZ/1bT4TdV/7V611urz+3OnS+rfE1JmvWt+uaWO1AAAAAAEIgGCgAAAAAC0UABAAAAQCAaKAAAAAAIRAMFAAAAAIFooAAAAAAgEA0UAAAAAASigQIAAACAQDRQAAAAABCIBgoAAAAAApVqoMzscjN70szuT1UQ0ElkFrkhs8gNmUVuyCxilb0DdYWksxPUAVTlCpFZ5OUKkVnk5QqRWeTlCpFZRCjVQLn7zZKeSlQL0HFkFrkhs8gNmUVuyCxi8TtQAAAAABCo4w2UmS01sxEzG9m1a3unlwNK2yuz42QW9bdnZhvbySzqb8/MToySWeRhr+uDMXI7k3W8gXL3Ze4+7O7D/f1DnV4OKG2vzA6QWdTfnpntHSKzqL89M9s3i8wiD3tdHwyS25mMH+EDAAAAgEBltzG/UtJtko41s/Vm9odpygI6g8wiN2QWuSGzyA2ZRay+MoPd/fxUhQBVILPIDZlFbsgsckNmEYsf4QMAAACAQDRQAAAAABCIBgoAAAAAAtFAAQAAAEAgGigAAAAACEQDBQAAAACBaKAAAAAAIBANFAAAAAAEooECAAAAgEB9VS5mTVf/tl1VLilJOvKy/srX3Pq+LZWvKUlz+4YrX9Nvu63yNatiW3Zo9rX3VL7uwkNeXvmaT53UrHxNSTroZcdXvqb9+JbK16zKwOM7dPQH7658XT/52MrXfOS9VvmakrT/tadUvmbjW/9R+ZpV6X1quxZ86a7K193cX/3r7JZjKl9SknRAs9Gdhaexnmd2aL9v3lv5ujt3vqzyNdef0Z37Lceuqf7rij0cdn3AHSgAAAAACEQDBQAAAACBaKAAAAAAIBANFAAAAAAEooECAAAAgEA0UAAAAAAQiAYKAAAAAALRQAEAAABAIBooAAAAAAhEAwUAAAAAgWigAAAAACBQ2w2UmS02s++Z2SozW2lmF6YsDEiNzCJH5Ba5IbPIDZlFrL4SYyckvc/d7zGz/STdbWbXu/sDiWoDUiOzyBG5RW7ILHJDZhGl7TtQ7v64u98z+fetklZJOiJVYUBqZBY5IrfIDZlFbsgsYpW5A/VzZnaUpJMl3dHi3FJJSyVp1sD8FMsBpQVnVnMqrQuYSlFuySzqiswiN1wfIETpTSTMbK6kr0q6yN23PPe8uy9z92F3H+7vHyq7HFBaVGZtsPoCgRamyu3emZ3VnQKB5wjPLK+zqIe46wNea2eyUg2UmfVrd9C+6O5fS1MS0DlkFjkit8gNmUVuyCxilNmFzyRdJmmVu38iXUlAZ5BZ5IjcIjdkFrkhs4hV5g7U6ZJ+X9IZZrZi8r9zEtUFdAKZRY7ILXJDZpEbMosobW8i4e63SLKEtQAdRWaRI3KL3JBZ5IbMIlbpTSQAAAAAYKaggQIAAACAQDRQAAAAABCIBgoAAAAAAtFAAQAAAEAgGigAAAAACEQDBQAAAACBaKAAAAAAIBANFAAAAAAE6qt0tR07ZSMPVLqkJPVNTFS+5ty+4crXlKRH3lB9Tzy+qvIlq+OSdyE/B1x2e+VrHvSy4ytfU5Ie/a/7V77m2BO9la9ZGXf5rvHq1+3Ca/v+155S+ZqStOW12ypfs3Fzs/I1K+OSN73yZQ/43J3Vr9lsVL6mJD1x4WldWVeXfKU761bBXT42Vvmys75VfW6PXXNs5WtK0k/OW1j5mmP/GtYacQcKAAAAAALRQAEAAABAIBooAAAAAAhEAwUAAAAAgWigAAAAACAQDRQAAAAABKKBAgAAAIBANFAAAAAAEIgGCgAAAAAC0UABAAAAQCAaKAAAAAAI1HYDZcvtLlMAAARqSURBVGazzOxOM/uhma00s79LWRjQCeQWuSGzyA2ZRW7ILGL1lRg7JukMd99mZv2SbjGza9399kS1AZ1AbpEbMovckFnkhswiStsNlLu7pG2TH/ZP/ucpigI6hdwiN2QWuSGzyA2ZRaxSvwNlZr1mtkLSk5Kud/c7WjxmqZmNmNnILh8rsxyQxL5yu1dmRWbRfWQWuSGzyE30NS25ndFKNVDu3nD3kyQtknSKmb2kxWOWufuwuw/322CZ5YAk9pXbvTIrMovuI7PIDZlFbqKvacntjJZkFz53f0bSTZL+/3bun8WuKgrj8LuIFoKChSmGGPQrBMTGzj8QbLRV8BsIpvSz2KUQJKCVjZW1qEGEGJUgiBMtFAvtRNgWM0XAKc7cuXfO2Xs/DwzMDCSzNvvXLO499+Y+/j+4DLqlN5qlN5qlN5pliYt8Ct/Vqnr69Psnkrya5Pt9DQaHoFt6o1l6o1l6o1nO6yKfwneU5HZVXcnJInantfbpfsaCg9EtvdEsvdEsvdEs53KRT+H7NsmNPc4CB6dbeqNZeqNZeqNZzmsvz0ABAADMwAIFAACwkAUKAABgIQsUAADAQhYoAACAhSxQAAAAC1mgAAAAFrJAAQAALGSBAgAAWKhaa5f3x6p+T/Lzjv/8mSR/7HGcLevtrM+11q6uPcQhaHax3s6q2bP1do8X0dtZNXu23u7xIno8q27/r8d73FWPZ13U7KUuUBdRVV+11l5Ye47LMNNZRzbTPc501pHNdI8znXVkM93jTGcd2Uz3OPJZvYUPAABgIQsUAADAQj0tUB+sPcAlmumsI5vpHmc668hmuseZzjqyme5xprOObKZ7HPas3TwDBQAAsLaeXoECAABY1eYXqKq6WVU/VNWDqnp/7XkOpaquV9XnVXW/qu5V1Xtrz8RuZmk20e1IZulWs+PQLL3R7Dg2/Ra+qrqS5MckryU5TvJlkrdaa9+tOtgBVNVRkqPW2t2qeirJ10neHPGsI5up2US3o5ipW82OQbOa7Y1mx2p2669AvZjkQWvtp9baP0k+SvLGyjMdRGvtt9ba3dPv/05yP8m1dadiB9M0m+h2INN0q9lhaJbeaHYgW1+griX55ZGfjzPYBZylqp5PciPJF+tOwg6mbDbRbeem7FazXdMsvdHsQLa+QNUZv9vuew73oKqeTPJxkluttb/Wnodzm67ZRLcDmK5bzXZPs/RGswPZ+gJ1nOT6Iz8/m+TXlWY5uKp6PCehfdha+2TtedjJVM0muh3EVN1qdgiapTeaHcjWP0TisZw8cPdKkoc5eeDu7dbavVUHO4CqqiS3k/zZWru19jzsZqZmE92OYqZuNTsGzdIbzY5l069Atdb+TfJuks9y8gDanRFDO/VSkneSvFxV35x+vb72UJzPZM0muh3CZN1qdgCa1WxvNDtWs5t+BQoAAGBLNv0KFAAAwJZYoAAAABayQAEAACxkgQIAAFjIAgUAALCQBQoAAGAhCxQAAMBCFigAAICF/gPiN5ewO7PYkgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1080x720 with 20 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "fig, axes = plt.subplots(4, 5, figsize=(15, 10), squeeze=True)\n",
    "\n",
    "for k in range(5):\n",
    "    for b in [1, 2, 3, 4]:\n",
    "        e_ij = e_k_b_ij[k][b - 1]\n",
    "        axes[b - 1, k].imshow(e_ij, vmax=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "82d777a8-ddc5-4b4a-b320-dbf1e49df8f3",
   "metadata": {},
   "source": [
    "Above, each column is a different spin-z=0 eigenstate, and each row a different blocking. The diagonal of each plot shows the entanglement of each block with its whole environment, and the off-diagonal shows the entanglement with other blocks."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}