Negativity_noninteracting-2D.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from scipy.linalg import block_diag\n",
    "import random\n",
    "from math import pi, sqrt\n",
    "import sys\n",
    "import os.path\n",
    "import matplotlib.pyplot as plt\n",
    "import time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def hopp_vec(Mx,My,vec,px,py):\n",
    "    # generates hoping matrix on square lattice\n",
    "    # vec= hopping vector\n",
    "    # py=1 PBC along y otherwise OBC\n",
    "    # px=1 PBC along x otherwise OBC\n",
    "    \n",
    "    hopin=np.zeros(shape=(Mx*My,Mx*My),dtype=complex)\n",
    "    for i_y in range(My):\n",
    "        for i_x in range(Mx):\n",
    "            i_px=i_x+vec[0]\n",
    "            i_py=i_y+vec[1]\n",
    "            if ((i_px>=Mx or i_px<0) and i_py< My and i_py>=0):\n",
    "                if i_px>=Mx :\n",
    "                    hx=px\n",
    "                else:\n",
    "                    hx=px.conjugate()\n",
    "                i_px=i_px%Mx\n",
    "                hopin[i_px+i_py*Mx,i_x+i_y*Mx]=hx\n",
    "            elif ((i_py>=My or i_py<0) and i_px< Mx and i_px>=0):\n",
    "                if i_py>=My :\n",
    "                    hy=py\n",
    "                else:\n",
    "                    hy=py.conjugate()\n",
    "                    \n",
    "                i_py=i_py%My\n",
    "                hopin[i_px+i_py*Mx,i_x+i_y*Mx]=hy\n",
    "                    \n",
    "            elif (i_px>=Mx or i_px<0) and (i_py>=My or i_py<0):\n",
    "                if i_px>=Mx:\n",
    "                    hx=px\n",
    "                else:\n",
    "                    hx=px.conjugate()\n",
    "    \n",
    "                if i_py>=My:\n",
    "                    hy=py\n",
    "                else:\n",
    "                    hy=py.conjugate()\n",
    "                i_px=i_px%Mx\n",
    "                i_py=i_py%My\n",
    "                hopin[i_px+i_py*Mx,i_x+i_y*Mx]=hx*hy\n",
    "                \n",
    "            elif i_px>=0 and i_px<Mx and i_py>=0 and i_py<My :\n",
    "                hopin[i_px+i_py*Mx,i_x+i_y*Mx]=1\n",
    "                    \n",
    "    return hopin\n",
    "\n",
    "\n",
    "# generates single particle states\n",
    "def wf_gen(Lx,Ly,px,py,N_pt,t):\n",
    "    H1_hop=hopp_vec(Lx,Ly,[0,1],px,py)+hopp_vec(Lx,Ly,[1,0],px,py)\n",
    "    H1=-t*(H1_hop+H1_hop.T)/2\n",
    "    energies, evecs= np.linalg.eigh(H1)\n",
    "    return evecs[:,:N_pt]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def exact_renyi_calc(r,GA,epsilon=1e-9):\n",
    "    chi0, _ =np.linalg.eigh(GA)\n",
    "    chi1=chi0[np.nonzero(np.abs(chi0)>epsilon)]\n",
    "    chi2=chi1[np.nonzero(np.abs(chi1-1)>epsilon)]\n",
    "#     return -np.sum((1-chi2)*np.log(1-chi2)+chi2*np.log(chi2))\n",
    "    return np.sum(np.log((1-chi2)**r+chi2**r))/(1-r)\n",
    "\n",
    "\n",
    "def exact_log_neg_calc(Na,Nb,GA,epsilon=1e-9):\n",
    "    \"\"\"computes log negativity of 1d free fermions\"\"\"\n",
    "    Nt=Na+Nb\n",
    "\n",
    "    Gm=np.eye(Nt)-2*GA\n",
    "    Sp=np.array([[0,1],[0,0]])\n",
    "    Gm_1=block_diag(-Gm[:Na,:Na],Gm[Na:,Na:])\\\n",
    "        +np.kron(Sp,-1j*Gm[:Na,Na:])+np.kron(Sp.T,-1j*Gm[Na:,:Na])\n",
    "    Gm_2=block_diag(-Gm[:Na,:Na],Gm[Na:,Na:])\\\n",
    "        +np.kron(Sp,1j*Gm[:Na,Na:])+np.kron(Sp.T,1j*Gm[Na:,:Na])    \n",
    "    Gx=np.dot(np.eye(Nt)-Gm_2,np.dot(np.linalg.inv(np.eye(Nt)+np.dot(Gm_1,Gm_2)),np.eye(Nt)-Gm_1))/2\n",
    "    Gxh=(Gx+np.transpose(Gx).conj())/2\n",
    "\n",
    "    chi0 =np.linalg.eigvalsh(Gxh)\n",
    "    i0=np.nonzero(np.abs(chi0)>epsilon)[0]\n",
    "    chi1=chi0[i0]\n",
    "    i1=np.nonzero(np.abs(chi1-1)>epsilon)[0]\n",
    "    nu=chi1[i1]\n",
    "    eA=np.sum(np.log(nu**0.5+(1-nu)**0.5))\n",
    "    \n",
    "    chi,VR=np.linalg.eigh(GA)\n",
    "    sA=np.sum(np.log(chi**2+(1-chi)**2))/2    \n",
    "\n",
    "    return eA+sA\n",
    "\n",
    "def exact_Renyi_neg_calc2(r,Na,Nb,GA,epsilon=1e-9):\n",
    "    \"\"\"computes log negativity of 1d free fermions only works for even r\"\"\"\n",
    "    Nt=Na+Nb\n",
    "\n",
    "    Gm=np.eye(Nt)-2*GA\n",
    "    Sp=np.array([[0,1],[0,0]])\n",
    "    Gm_1=block_diag(-Gm[:Na,:Na],Gm[Na:,Na:])\\\n",
    "        +np.kron(Sp,-1j*Gm[:Na,Na:])+np.kron(Sp.T,-1j*Gm[Na:,:Na])\n",
    "    Gm_2=block_diag(-Gm[:Na,:Na],Gm[Na:,Na:])\\\n",
    "        +np.kron(Sp,1j*Gm[:Na,Na:])+np.kron(Sp.T,1j*Gm[Na:,:Na])    \n",
    "    Gx=np.dot(np.eye(Nt)-Gm_2,np.dot(np.linalg.inv(np.eye(Nt)+np.dot(Gm_1,Gm_2)),np.eye(Nt)-Gm_1))/2\n",
    "    Gxh=(Gx+np.transpose(Gx).conj())/2\n",
    "\n",
    "    chi0 =np.linalg.eigvalsh(Gxh)\n",
    "    i0=np.nonzero(np.abs(chi0)>epsilon)[0]\n",
    "    chi1=chi0[i0]\n",
    "    i1=np.nonzero(np.abs(chi1-1)>epsilon)[0]\n",
    "    nu=chi1[i1]\n",
    "    eA=np.sum(np.log(nu**(r/2)+(1-nu)**(r/2)))\n",
    "    \n",
    "    chi,VR=np.linalg.eigh(GA)\n",
    "    sA=r*np.sum(np.log(chi**2+(1-chi)**2))/2    \n",
    "\n",
    "    return eA+sA\n",
    "\n",
    "def exact_Renyi_neg_calc(r,N,Na,Nb,V1):\n",
    "    '''computes moments of partial transpose tr(\\rho^T \\rho^T\\dag ... )\n",
    "    for 1d free fermions ''' \n",
    "    k_sw=np.arange((-(r-1)/2),(r-1)/2+0.1,1)\n",
    "    Zk=0\n",
    "    if (r % 2) == 0:\n",
    "        delta=pi\n",
    "    else:\n",
    "        delta= pi*(r-1)/r\n",
    "#     delta=pi\n",
    "\n",
    "    Nflip=Na+Nb\n",
    "    N1=int((N-Nflip)/2)\n",
    "\n",
    "    for i_k in range(len(k_sw)):\n",
    "#         op_k=block_diag(np.eye(N1),np.exp(1j*2*pi/r*k_sw[i_k])*np.eye(Na),np.exp(1j*delta-1j*2*pi/r*k_sw[i_k])*np.eye(Nb),np.eye(N1))\n",
    "        op_k=block_diag(np.exp(1j*2*pi/r*k_sw[i_k])*np.eye(Na),np.exp(1j*delta-1j*2*pi/r*k_sw[i_k])*np.eye(Nb),np.eye(N-Nflip))\n",
    "        Zk +=np.real(np.log( np.linalg.det( np.dot(np.dot(np.matrix(V1).H,op_k),V1) )))\n",
    "    return Zk\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "Lx=20\n",
    "Ly=40\n",
    "N=Lx*Ly\n",
    "Npt=int(N/2)\n",
    "t = 1  # hopping amplitude\n",
    "px=-1 # BC along x\n",
    "py=-1 # BC along y\n",
    "\n",
    "V1 = wf_gen(Lx,Ly,px,py,Npt,t)\n",
    "\n",
    "Lysub_list=np.arange(2,Ly/2,1)\n",
    "\n",
    "\n",
    "def ent_plotter(r):\n",
    "    '''r=2 # Renyi/ Negativity index'''\n",
    "\n",
    "    Rr_ex=np.zeros(len(Lysub_list)) # Renyi entropy\n",
    "    RNr_ex=np.zeros(len(Lysub_list)) # Renyi negativity\n",
    "    RNr2_ex=np.zeros(len(Lysub_list)) # Renyi negativity 2\n",
    "    LN_ex=np.zeros(len(Lysub_list)) # Log negativity\n",
    "\n",
    "    Gmat=np.dot(V1,np.matrix(V1).H)\n",
    "\n",
    "    for i_L in range(len(Lysub_list)):\n",
    "        Lysub=Lysub_list[i_L]\n",
    "        inds_A = np.arange(0,int(Lx*Lysub))\n",
    "        inds_AB = np.arange(0,int(2*Lx*Lysub))\n",
    "        Nsub =  int(Lx*Lysub)\n",
    "        Rr_ex[i_L] = exact_renyi_calc(r,Gmat[np.ix_(inds_AB,inds_AB)])\n",
    "        RNr_ex[i_L] = exact_Renyi_neg_calc(r,N,Nsub,Nsub,V1)\n",
    "        RNr2_ex[i_L] = exact_Renyi_neg_calc2(r,Nsub,Nsub,Gmat[np.ix_(inds_AB,inds_AB)])\n",
    "    plt.figure()\n",
    "#     plt.plot(Lsub_list,-Rr_ex,'o',label='(-1) Renyi Entropy')\n",
    "    plt.plot(Lysub_list,(RNr_ex),'v',label='Renyi Neg')\n",
    "    if r % 2 ==0:\n",
    "        plt.plot(Lysub_list,(RNr2_ex),'x',label='Renyi Neg 2')\n",
    "    plt.title('r=%d Renyi index' % (r))\n",
    "    plt.xlabel('$\\ell$')\n",
    "    plt.legend()\n",
    "    plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAEYCAYAAAC5q4bCAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAeEUlEQVR4nO3dfXRU9b3v8fc3IRqtCKigAiJpiVaDQkP0Gm9BbyVCLRHQk9KKKGJX2lUr0NoWvXh77Fpyq8VzCmju8XgOiFYX+HCbqwWtDfZ6sEc8NSBaUvRWi7WhKSBHoYhWIN/7xwxxDJOnecie+eXzWmtWJvvxmz17PvnNb+/Z29wdEREJU0HUBYiISPYo5EVEAqaQFxEJmEJeRCRgCnkRkYAp5EVEAqaQF8kQM/vvZvavHYwbb2avp7FsN7NRqVcnfZVCXnKemR1tZsvN7I9m9lcze9nMvtjJ9LPN7JCZ7TOzvWb2iplNyXad7v4/3f1rHYx73t3PzHYNIu0p5CUSZtavB5P3A/4EXAQMAP4H8KiZjexkng3ufhwwEPhfwGozG5hSsSJ5TCEvvcbM3jKzBWb2KvB+d4Pe3d9399vc/S13b3X3NcA2YFw35m0Ffgp8CiiN13G0md1lZm+b2Q4zu9fMjomPu9jMms3sJjPbaWYtZnZdfNx58enb6jazK81sc/z5bWb2UAd/+8Vm1txuW3zXzF41sz1m9oiZFSeM/1583X82szntltVZ/U+Z2T8kTPuIma3oajtJuBTy0tu+CnwJGOjuB81sjZm918FjTbIFmNnJwBlAU1crM7NC4DrgAPDH+OA74/OPBUYBw4AfJMx2CrFPDMOA64E6Mxvk7i8Bu4GqhGmvJvZPJBVfBiYDJcC5wOx4zZOB78bXUwpMbDdfZ/XPAWaZ2RfMbCZwHjAvxfokBO6uhx698gDeAuakuYwiYB3wz51MMxs4CLxHLNw/AL4cH2fA+8BnEqavBLbFn18cn75fwvidwAXx5wuAh+PPTwD2A6fGf78NeKiDmi4Gmttti6sTfv8xcG/8+QrgjoRxZwBOLNA7rT/++xXEurfeAT4f9euuR7QPteSlt/0p1RnNrIBYq/kj4FtdTP6iuw8EBgFPAuPjwwcDxwIbD39iAH4RH37Ybnc/mPD7fuC4+POHgGozO45YS/x5d29J8U/6SwfrGMont9MfE553p/41QCHwurv/OsXaJBAKeeltn7jsqZk9HT8LJtnj6YTpDFgOnAxc6e4HurUy933AN4l1YXyOWOv2A6DM3QfGHwM8dpC2O8vbDmwApgOzSL2rpjMtwGkJv49IeN6d+hcBW4FTzeyrWahP8ohCXiLl7l909+M6eCSeJvlPwFlAtbt/0MN17Ab+FfiBxw7E/gvwEzMbAmBmw8xsUg8W+SDwfeAcoL4ntXTTo8BsMzvbzI4F/v7wiK7qN7MJxI5BXBN/3G1mw7JQo+QJhbzkPDM7Hfg6sQONf0lo6c/swWKWAJeZ2bnE+tXfAF40s73E+vh7cg57PXA6UO/u7/dgvm5x96eJ1fureJ2/ajdJ0vrN7Hhi/4C+5e7b4101y4H745+EpA8yd900RKSnzOxN4Ovuvi7qWkQ6o5a8SA+Z2ZXEji20b2GL5JyefOtQpM8zs+eAs4FZ8f5xkZym7hoRkYCpu0ZEJGA51V1z0kkn+ciRI6MuQ0Qkr2zcuPEddx+cbFxOhfzIkSNpbGyMugwRkbxiZn/saJy6a0REAqaQFxEJmEJeRCRgOdUnn8yBAwdobm7mww8/jLqUPqm4uJjhw4dTVFQUdSkikoKcD/nm5mb69+/PyJEj0eU3epe7s3v3bpqbmykpKYm6HBFJQc6H/IcffqiAj4iZceKJJ7Jr166oSxEJ0oo759GwZxgbWsvahlUWNFE1YDtzFizNyDryok9eAR8dbXuR7PGh5dQVLaOyIHYny8qCJuqKluFDyzO2jpxvyYuIhKp66gzmL97FPUVLeOjQRK4uXMe3W+ezeNqMjK0jL1ryUSssLGTs2LGMHj2a6upq3nvvvYwuv7Gxkblz5x4x/LnnnsPM+PnPf942bMqUKTz33HMZXb+IRGPI8cWMKJ/EqtYq5vWrZ1VrFaeNm8yQ/sUZW0daLXkzqyF28+KzgPPdvTE+vAq4AziK2P04v+fuWb8s62VLn+d3LXuPGH72qcfz1LzxSebonmOOOYbNmzcDcO2111JXV8fChQtTXl57FRUVVFRUJB03fPhwFi1aRHV1dcbWJyK546bSHfjmBpYenM6swgasdCYwOmPLT7clv4XYneHXtxv+DrHbtJ0DXEt27oN5hPIRAykq/GQfclGhUX76oIyto7Kyku3btwPw5ptvMnnyZMaNG8f48eN57bXXAJg9ezZz587lwgsv5NOf/jSPP/44ALNmzeKJJ55oW9bMmTN58sknee6555gyZUrS9Y0ZM4YBAwbQ0NBwxLiNGzdy0UUXMW7cOCZNmkRLS+x+0i+99BLnnnsulZWVfO9732P06MztMCKSQdvWM+ipWupHLWLJoRrqRy1i0Npa2NY+UlOXVsi7+1Z3fz3J8Jfd/c/xX5uAYjM7Op11dcfcS0opaHegsNCMuZeMysjyDx06xLPPPsvll18OQG1tLXfffTcbN27krrvu4pvf/GbbtC0tLfz6179mzZo13HzzzQB87Wtf4/777wdgz549vPDCC1x22WVdrvfWW2/l9ttv/8SwAwcOcOONN/L444+zceNG5syZ0/bp4rrrruPee+9lw4YNFBYWZuRvF5Es2L4JalZSPXUG5408geppM6BmZWx4hvTGgdcrgZfd/W/JRppZLVALMGLEiGSTdNuQ44upGTecRxr/xIFDTlGh8XcVp6Xdv/XBBx8wduxY3nrrLcaNG0dVVRX79u3jhRdeoKampm26v/3t4z9x2rRpFBQUcPbZZ7Njxw4ALrroIm644QZ27tzJz372M6688kr69ev6JRg/PtbV9Pzzz7cNe/3119myZQtVVVVA7B/Qqaeeynvvvcdf//pXLrzwQgCuuuoq1qxZk9bfLyJZ8vn5AAwBHv16ZWxY/wlQMiFjq+gyYcxsHXBKklEL3f2JJMMT5y0D7gQu7Wgad78PuA+goqIi7TuYzL2klMc2NgOesVb84T75PXv2MGXKFOrq6pg9ezYDBw5s66tv7+ijP/7gknhjllmzZvHwww+zevVqVqxY0e0aFi5cyKJFi9r+Kbg7ZWVlbNiw4RPTvfvuuz34y0QkdF1217j7RHcfneTRVcAPJ3ZX+2vc/c1MFdyVw615MzLSik80YMAAli1bxl133cUxxxxDSUkJjz32GBAL3VdeeaXLZcyePZslS5YAUFZW1vnECS699FLefffdtnWceeaZ7Nq1qy3kDxw4QFNTE4MGDaJ///68+OKLAKxevbonf6KIBCYrp1Ca2UBgLXCLu/97NtbRmbmXlHLeyBMy1hef6HOf+xxjxoxh9erVPPzwwyxfvpwxY8ZQVlb2iYOqHTn55JM566yzuO6663q87oULF9Lc3AzAUUcdxeOPP86CBQsYM2YMY8eO5YUXXgBg+fLl1NbWUllZibszYMCAHq9LRMKQ1j1ezWw6cDcwGHgP2Ozuk8zsVuAW4PcJk1/q7js7W15FRYW3v2nI1q1bOeuss1KuMdfs37+fc845h02bNmUtfPft28dxxx0HwB133EFLSwtLl6b+FenQXgOR0JjZRndPeh52Wgde3b2eWJdM++G3A7cfOUfftm7dOubMmcN3vvOdrLau165dy49+9CMOHjzI6aefzsqVK7O2LhHJbbqsQS+aOHEib7/9dtbXM2PGDGbMyNzXokUkf+myBiIiAVPIi4gETCEvIhIwhbyISMAU8t0Q4qWGN2/eTGVlJWVlZZx77rk88sgjaS9TRHJPWCH/6yVHXr1t2/rY8DQcvqzBli1bOOGEE6irq0tree1VVFSwbNmypOMOX2o404499lgefPBBmpqa+MUvfsH8+fMz/s9LRKIXVsgPK4fHZn8c9NvWx34flrlbaYVyqeEzzjiD0tJSAIYOHcqQIUN0L1eRAIUV8iUTYpfpfGw2/GpR7GfNyoxd0S3USw3/5je/4aOPPuIzn/lM9zaEiOSN8L4MVTIBKq6H9T+GCd/PSMCHfKnhlpYWZs2axQMPPEBBQVj/80UkxJDfth4al8cCvnE5lIxPO+hDvdTw3r17+dKXvsTtt9/OBRdc0O35RCR/hNV0O9wHX7MSvrDw466bDN1KK6RLDX/00UdMnz6da6655hOfRkQkLGGFfPxWWm0t98N99Bm8lVYolxp+9NFHWb9+PStXrmTs2LGMHTu2w08lIpK/0rrUcKbpUsOZoUsNi/QtnV1qOKyWfI5bt24dn/3sZ7nxxhuzfqnhw1/eev7557n11luzti4RyW3hHXjNYbrUsIj0trxoyedSl1Jfo20vkt9yPuSLi4vZvXu3wiYC7s7u3bspLs7czdBFpHflfHfN8OHDaW5u1lfuI1JcXMzw4cOjLkNEUpTzIV9UVERJSUnUZYiI5KWc764REZHUKeRFRAKmkBcRCZhCXkQkYAp5EZGAKeRFRAKmkBcRCZhCXkQkYAp5EZGAKeRFRAKmkBcRCZhCXkQkYAp5EZGAKeRFRAKWVsibWY2ZNZlZq5kdcRNZMxthZvvM7LvprEdERFKTbkt+C3AFsL6D8T8Bnk5zHSIikqK0bhri7lsBzOyIcWY2DfgD8H466xARkdRlpU/ezD4FLAB+2I1pa82s0cwadYs/EZHM6jLkzWydmW1J8pjayWw/BH7i7vu6Wr673+fuFe5eMXjw4J7ULiIiXeiyu8bdJ6aw3P8C/J2Z/RgYCLSa2Yfufk8KyxIRkRRl5Ube7j7+8HMzuw3Yp4AXEel96Z5COd3MmoFKYK2ZPZOZskREJBPSPbumHqjvYprb0lmHiIikTt94FREJmEJeRCRgCnkRkYAp5EVEAqaQFxEJmEJeRCRgCnkRkYAp5EVEAqaQFxEJmEJeRCRgCnkRkYAp5EVEAqaQFxEJmEJeRCRgCnkRkYAp5EVEApaV2//1phV3zqNhzzA2tJa1DassaKJqwHbmLFgaYWUiEqp8yp28b8n70HLqipZRWdAExDZ0XdEyfGh5xJWJSKjyKXfyviVfPXUG8xfv4p6iJTx0aCJXF67j263zWTxtRtSliUig8il38r4lP+T4YkaUT2JVaxXz+tWzqrWK08ZNZkj/4qhLE5FA5VPu5H1LHuCm0h345gaWHpzOrMIGrHQmMDrqskQkYPmSO3nfkmfbegY9VUv9qEUsOVRD/ahFDFpbC9vWR12ZiIQqj3In/1vy2zdBzUqqTzyfZ/a/TPW0S+CdktjwkglRVyciIcqj3DF3j7qGNhUVFd7Y2Bh1GSIiecXMNrp7RbJx+d9dIyIiHVLIi4gETCEvIhIwhbyISMAU8iIiAVPIi4gETCEvIhIwhbyISMAU8iIiAVPIi4gELK2QN7MaM2sys1Yzq2g37lwz2xAf/1szy71rcIqIBC7dC5RtAa4A/jlxoJn1Ax4CZrn7K2Z2InAgzXWJiEgPpRXy7r4VwMzaj7oUeNXdX4lPtzud9YiISGqy1Sd/BuBm9oyZbTKz73c0oZnVmlmjmTXu2rUrS+WIiPRNXbbkzWwdcEqSUQvd/YlOlvt54DxgP/Bs/FKYz7af0N3vA+6D2KWGu1u4iIh0rcuQd/eJKSy3Gfg3d38HwMyeAsqBI0JeRESyJ1vdNc8A55rZsfGDsBcBv8vSukREpAPpnkI53cyagUpgrZk9A+Du7wL/CLwEbAY2ufvaNGsVEZEeSvfsmnqgvoNxDxE7jVJERCKib7yKiARMIS8iEjCFvIhIwBTyIiIBU8iLiARMIS8iEjCFvIhIwBTyIiIBU8iLiARMIS8iEjCFvIhIwBTyIiIBU8iLiARMIS8iEjCFvIhIwBTyIiIBS+umIaFYcec8GvYMY0NrWduwyoImqgZsZ86CpRFWJiKZ1Bff62rJAz60nLqiZVQWNAGxF72uaBk+tDziykQkk/rie10teaB66gzmL97FPUVLeOjQRK4uXMe3W+ezeNqMqEsTkQzqi+91teSBIccXM6J8Eqtaq5jXr55VrVWcNm4yQ/oXR12aiGRQX3yvqyUfd1PpDnxzA0sPTmdWYQNWOhMYHXVZIpJhfe29rpY8wLb1DHqqlvpRi1hyqIb6UYsYtLYWtq2PujIRyaQ++F5XSx5g+yaoWUn1iefzzP6XqZ52CbxTEhteMiHq6kQkU/rge93cPeoa2lRUVHhjY2PUZYiI5BUz2+juFcnGqbtGRCRgCnkRkYAp5EVEAqaQFxEJmEJeRCRgCnkRkYAp5EVEAqaQFxEJmEJeRCRgCnkRkYClFfJmVmNmTWbWamYVCcOLzOwBM/utmW01s1vSL1VERHoq3Zb8FuAKoP0l3GqAo939HGAc8HUzG5nmukREpIfSugqlu28FMLMjRgGfMrN+wDHAR8DedNYlIiI9l60++ceB94EW4G3gLnf/z2QTmlmtmTWaWeOuXbuyVI6ISN/UZUvezNYBpyQZtdDdn+hgtvOBQ8BQYBDwvJmtc/c/tJ/Q3e8D7oPYpYa7W7iIiHSty5B394kpLPcq4BfufgDYaWb/DlQAR4S8iIhkT7a6a94GvmAxnwIuAF7L0rpERKQD6Z5COd3MmoFKYK2ZPRMfVQccR+zsm5eA+9391bQqFRGRHkv37Jp6oD7J8H3ETqMUEZEI6RuvIiIBU8iLiARMIS8iEjCFvIhIwNI68CqftOLOeTTsGcaG1rK2YZUFTVQN2M6cBUsjrEwkv+m9lTq15DPIh5ZTV7SMyoImILYT1hUtw4eWR1yZSH7Teyt1aslnUPXUGcxfvIt7ipbw0KGJXF24jm+3zmfxtBlRlyaS1/TeSp1a8hk05PhiRpRPYlVrFfP61bOqtYrTxk1mSP/iqEsTyWt6b6VOLfkMu6l0B765gaUHpzOrsAErnQmMjroskbyn91Zq1JLPpG3rGfRULfWjFrHkUA31oxYxaG0tbGt/TxUR6RG9t1Kmlnwmbd8ENSupPvF8ntn/MtXTLoF3SmLDSyZEXZ1I/tJ7K2XmnjuXcK+oqPDGxsaoyxARyStmttHdK5KNU3eNiEjAFPIiIgFTyIuIBEwhLyISMIW8iEjAFPIiIgFTyIuIBEwhLyISMIW8iEjAFPIiIgFTyIuIBEwhLyISMIW8iEjAFPIiIgFTyIuIBEwhLyISMIW8iEjAFPIiIgFTyIuIBEwhLyISMIW8iEjAFPIiIgFLK+TNbLGZvWZmr5pZvZkNTBh3i5m9YWavm9mktCsVEZEe65fm/A3ALe5+0MzuBG4BFpjZ2cBXgDJgKLDOzM5w90Nprq9PWHHnPBr2DGNDa1nbsMqCJqoGbGfOgqURVibSM9qXo5dWS97df+nuB+O/vggMjz+fCqx297+5+zbgDeD8dNbVl/jQcuqKllFZ0ATE3hR1RcvwoeURVybSM9qXo5duSz7RHOCR+PNhxEL/sOb4sCOYWS1QCzBixIgMlpO/qqfOYP7iXdxTtISHDk3k6sJ1fLt1PounzYi6NJEe0b4cvS5b8ma2zsy2JHlMTZhmIXAQePjwoCSL8mTLd/f73L3C3SsGDx6cyt8QnCHHFzOifBKrWquY16+eVa1VnDZuMkP6F0ddmkiPaF+OXpcteXef2Nl4M7sWmAJc4u6Hg7wZOC1hsuHAn1Mtsi+6qXQHvrmBpQenM6uwASudCYyOuiyRHtO+HK10z66ZDCwALnf3/QmjngS+YmZHm1kJUAr8Jp119Snb1jPoqVrqRy1iyaEa6kctYtDaWti2PurKRHpG+3Lk0u2Tvwc4GmgwM4AX3f0b7t5kZo8CvyPWjXODzqzpge2boGYl1SeezzP7X6Z62iXwTklseMmEbi9GZzZIqjK272RoX5bU2cc9LNGrqKjwxsbGqMsIxvKfPsD0NxZyw4G5bGgtazuzoX7UIq6fdW3U5UkO076TX8xso7tXJBuXybNrJMfozAZJlfadcOiyBgHTmQ2SKu074VBLPnA6s0FSpX0nDGrJh0xnNkiqtO8EQy35kOnMBkmV9p1g6OwaEZE8p7NrJC063z5/6LWS9tQnL13SlQTzh14raU8teelSqOdMh9jqDfW1ktSpJS9dCvWc6RBbvaG+VpI6teSlW0I8ZzrUVm+Ir5WkTiEvXYufM7181CKW/O4kjjvzv3H92lo4dmVeXzCtrdW7qYl5/X7GPYeuiLTVm5Htk6HXSsKh7hrp2uFzpqfO4LyRJ1A9bQbUrIwN74Fc7B65qXQHVxXEWr1XFTRwU+lfIqslI9snQ6+VhEPnyUuv2bn3Q767+B5+UvBx98h3Wuez+Pvf6lHrOWOfCLath8dms/zUv+f2353ErWe/w/UtP4yFYgSfUDK1faTv6ew8ebXkpddk6qBgxj4R5NgnFB00lWxQn7z0qkwcFMzYAdPPzwdgCPDo1ytjw/pP6HHfdSYP4OqgqWSaWvLSezJ00atca/FmrB5dFEyyQC156T0ZvOhVrrV4M1KPLgomWaADr5J/MnTANNh6pM/RBcokLLnW4s21ekQSqCUvIpLndAqliEgfpZAXEQmYQl5EJGAKeRGRgCnkRUQCllNn15jZLuCPWVzFScA7WVx+Nqjm3qGasy/f6oX8qfl0dx+cbEROhXy2mVljR6cZ5SrV3DtUc/blW72QnzW3p+4aEZGAKeRFRALW10L+vqgLSIFq7h2qOfvyrV7Iz5o/oU/1yYuI9DV9rSUvItKnKORFRAIWXMib2Wlm9n/NbKuZNZnZvCTTXGxme8xsc/zxgyhqbVfTW2b223g9R1yK02KWmdkbZvaqmfXwhqaZZWZnJmy/zWa218zmt5sm8u1sZivMbKeZbUkYdoKZNZjZ7+M/B3Uw72Qzez2+zW+OsN7FZvZa/HWvN7OBHczb6T7UyzXfZmbbE177yzqYt9e3cSc1P5JQ71tmtrmDeSPZzilz96AewKlAefx5f+D/AWe3m+ZiYE3Utbar6S3gpE7GXwY8DRhwAfAfUdecUFsh8BdiX8jIqe0MTADKgS0Jw34M3Bx/fjNwZwd/05vAp4GjgFfa70e9WO+lQL/48zuT1dudfaiXa74N+G439pte38Yd1dxu/D8AP8il7ZzqI7iWvLu3uPum+PO/AluBYdFWlRFTgQc95kVgoJmdGnVRcZcAb7p7Nr+tnBJ3Xw/8Z7vBU4EH4s8fAKYlmfV84A13/4O7fwSsjs+XVcnqdfdfuvvB+K8vAsOzXUdPdLCNuyOSbQyd12xmBnwZWNUbtWRbcCGfyMxGAp8D/iPJ6Eoze8XMnjazst6tLCkHfmlmG82sNsn4YcCfEn5vJnf+eX2Fjt8QubadAU529xaINQqAIUmmydXtPYfYJ7pkutqHetu34l1MKzroEsvVbTwe2OHuv+9gfK5t504FG/Jmdhzwv4H57r633ehNxLoWxgB3A/+nl8tL5r+6eznwReAGM2t/3zhLMk/k57+a2VHA5cBjSUbn4nburpzb3ma2EDgIPNzBJF3tQ73pn4DPAGOBFmLdH+3l3DaO+yqdt+JzaTt3KciQN7MiYgH/sLv/rP14d9/r7vviz58CiszspF4us31Nf47/3AnUE/som6gZOC3h9+HAn3unuk59Edjk7jvaj8jF7Ry343BXV/znziTT5NT2NrNrgSnATI93DLfXjX2o17j7Dnc/5O6twL90UEtObWMAM+sHXAE80tE0ubSduyO4kI/3py0Htrr7P3YwzSnx6TCz84lth929V+UR9XzKzPoffk7sQNuWdpM9CVwTP8vmAmDP4S6HiHXY6sm17ZzgSeDa+PNrgSeSTPMSUGpmJfFPK1+Jz9frzGwysAC43N33dzBNd/ahXtPueNH0DmrJmW2cYCLwmrs3JxuZa9u5W6I+8pvpB/B5Yh/5XgU2xx+XAd8AvhGf5ltAE7Gj+S8CF0Zc86fjtbwSr2thfHhizQbUETsb4bdARQ5s62OJhfaAhGE5tZ2J/QNqAQ4QazleD5wIPAv8Pv7zhPi0Q4GnEua9jNjZWW8efk0iqvcNYn3Xh/fne9vX29E+FGHNP43vp68SC+5Tc2Ubd1RzfPjKw/tvwrQ5sZ1TfeiyBiIiAQuuu0ZERD6mkBcRCZhCXkQkYAp5EZGAKeRFRAKmkBcRCZhCXkQkYAp5kW6If9Pxnvi3jUXyhkJepHu+ARxN7BvVInlDIS/SPZOJff1+c8R1iPSIQl6kC2ZWTOwuRuXAv0VcjkiPKORFulZKLORfc/cDURcj0hP9oi5AJA8MBs6gl25NJ5JJasmLdG0osZvQFHRwGzuRnKWQF+lE/E5B5cApwL3AoWgrEukZXU9eRCRgasmLiARMIS8iEjCFvIhIwBTyIiIBU8iLiARMIS8iEjCFvIhIwP4/I0WKEtiPzfoAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAEYCAYAAAC5q4bCAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAa1UlEQVR4nO3dfZRU9Z3n8feHpmOb8KQRItpgo7ZGUEFoPcEN6iw+hUFMxhhNHEZ0ctCN2mZ3TXSCm/HMhF3jkDXGccZjRrMzE2ZRMSaG6JwJmSFDjjKRNoD4LAlqa0cQH5CoLJLv/lG3sWiqH+vW063P65w6Vt2Hut+6lJ/+1e/+7r2KCMzMLJuGVboAMzMrHYe8mVmGOeTNzDLMIW9mlmEOeTOzDHPIm5llmEPeLCWSvibp73qZN0vSM0W8d0g6cujVWb1yyFtNkPR9SV2Stkt6VtIX+1h2gaTdknYky6+XNLfUNUbE/4yIgnVFxOqIOLrUNZj15JC3ipA0fJCr/C+gJSJGAfOAb0ia0cfyj0TECGAM8DfAMkljhlKrWS1zyFvZSNos6VpJG4DfDSboI+KJiNjZ/TJ5HDGA9X4P/CPwEaA1qWM/SUskvSjpVUm3S9o/mXeapE5J/13SluTXwyXJvBOT5ffULek8SeuS5zdI+n4vn/00SZ099sU1kjZIekvS3ZKa8uZ/Jdn2K5Iu7fFefdX/oKRv5S17t6S7+ttPll0OeSu3zwN/CIyJiPclrZD0Zi+PFfkrSvobSe8ATwNdwIP9bUxSA3AJsAt4IZn8TeAoYBpwJHAo8PW81Q4GRifT/xS4TdIBEfEosA04I2/ZPyb3R2QoPgecDUwCjgcWJDWfDVyTbKcVOL3Hen3VfykwX9J/lnQRcCJw9RDrsyyICD/8KMsD2AxcWuR7NACfBK4HGntZZgHwPvAmuXB/F/hcMk/A74Aj8pafCfwmeX5asvzwvPlbgE8kz68FlibPDwTeAcYnr28Avt9LTacBnT32xR/nvb4JuD15fhdwY968o8j9cjmyv/qT138EvAS8Bnyy0v/uflT24Za8ldtLxawcEbsj4hdAM/Bf+lh0TUSMAQ4AHgBmJdPHAh8GOrp/MQD/nEzvti0i3s97/Q4wInn+feAcSSPItcRXR0TXED/Ob3vZxiHsvZ9eyHs+kPpXkPtj+Eyyr6yOOeSt3Pa67Kmkh5JRMIUeD/XxPsMZWJ/8DuBL5LowTiDXun0XmBIRY5LH6MgdpO2/+IiXgUeAzwDzGXpXTV+6gAl5ryfmPR9I/YuBp4Dxkj5fgvqshjjkraIi4lMRMaKXx6cAJI2TdKGkEZIaJJ1Frm//Xwe4jW3A3wFfj9yB2O8CN0sal7z/ocl7DtQ/AF8FjgPuH8R6A3UPsEDSZEkfBv68e0Z/9Us6hdwxiD9JHrdKOrQENVqNcMhbLQhyXTOdwBvAEuDLEfGjQbzHt4E5ko4n16/+PLBG0nZgJTCYMez3A4cB90fE7wax3oBExEPk6v3XpM6ef8wK1i9pFLk/QFdGxMtJV82dwPckKe06rTYowjcNMRssSZuAyyJiZaVrMeuLW/JmgyTpPHK/LgbUXWRWSYM969CsrklaBUwG5if942ZVzd01ZmYZ5u4aM7MMq6rumoMOOihaWloqXYaZWU3p6Oh4LSLGFppXVSHf0tLC2rVrK12GmVlNkfRCb/PcXWNmlmEOeTOzDHPIm5llWFX1yZtZNuzatYvOzk7ee++9SpeSKU1NTTQ3N9PY2DjgdYoKeUl/CZwL/J7cNbcXRMQrklrIXQWv+8bFayLi8mK2ZWa1o7Ozk5EjR9LS0oIvm5OOiGDbtm10dnYyadKkAa9XbEv+ryLifwBIaid3d5ruMN8UEdOKfP9+zbllNU92bd9n+uTxo3jw6lkF1jCzUnvvvfcc8CmTxEc/+lG2bt06qPWK6pOPiPx0/Qg9rhVeDtMnjqGxYe8vUmODmH7YAeUuxczyOODTN5R9WvSBV0mLJb0EXMTe98mcJOlXkn4uqdcmtaSFktZKWjvYv1AA7bNbGdbjgzdItM8+ctDvZWaWNf2GvKSVkjYWeJwLEBGLImICsBS4MlmtC5gYEScA/w34p+Ra1/uIiDsioi0i2saOLXjCVp/GjWri/BnNe1rzjQ3is20TGDeyqZ81zSzLGhoamDZtGsceeyznnHMOb775Zqrvv3btWtrb2/eZvmrVKiTx4x//eM+0uXPnsmrVqlS3P1D99slHRM87xffmn4CfAH8eETuBncn6Hcm1t48CSnI6a/vsVu7t6ATCrXizGlOq42r7778/69atA+Diiy/mtttuY9GiRUN+v57a2tpoa2srOK+5uZnFixdzzjnnpLa9oSqqu0ZSa97LecDTyfSxkhqS54cDrcCvi9lWX7pb8xJuxZvVmHIcV5s5cyYvv/wyAJs2beLss89mxowZzJo1i6effhqABQsW0N7ezsknn8zhhx/O8uXLAZg/fz4/+tEHNyG76KKLeOCBB1i1ahVz584tuL2pU6cyevRofvrTn+4zr6Ojg1NPPZUZM2Zw1lln0dWVuw/8o48+yvHHH8/MmTP5yle+wrHHHpvKZy+2T/7GpOtmA3AmcHUy/RRgg6T1wHLg8oh4vcht9al9disnthzoVrxZjSn1cbXdu3fzs5/9jHnz5gGwcOFCbr31Vjo6OliyZAlf+tKX9izb1dXFL37xC1asWMF1110HwBe/+EW+973vAfDWW2/x8MMPM2fOnH63e/311/ONb3xjr2m7du3iqquuYvny5XR0dHDppZfu+XVxySWXcPvtt/PII4/Q0NCQymeHIodQRsR5vUy/D7ivmPcerHGjmrjnspnl3KSZpaD7l/jda19i1+5I7bjau+++y7Rp09i8eTMzZszgjDPOYMeOHTz88MOcf/75e5bbuXPnnuef/vSnGTZsGJMnT+bVV18F4NRTT+WKK65gy5Yt/OAHP+C8885j+PD+o3PWrFxX0+rVq/dMe+aZZ9i4cSNnnHEGkPsDNH78eN58803efvttTj75ZAC+8IUvsGLFiqI+fzef8WpmFVeK42rdffJvvfUWc+fO5bbbbmPBggWMGTNmT199T/vtt9+e5/k3VJo/fz5Lly5l2bJl3HXXXQOuYdGiRSxevHjPH4WIYMqUKTzyyCN7LffGG28M4pMNjq9dY2YVV8rjaqNHj+Y73/kOS5YsYf/992fSpEnce++9QC50169f3+97LFiwgG9/+9sATJkyZcDbPvPMM3njjTf2bOPoo49m69ate0J+165dPPHEExxwwAGMHDmSNWvWALBs2bLBfMQ+OeTNrCqU8rjaCSecwNSpU1m2bBlLly7lzjvvZOrUqUyZMmWvg6q9+djHPsYxxxzDJZdcMuhtL1q0iM7OTgA+9KEPsXz5cq699lqmTp3KtGnTePjhhwG48847WbhwITNnziQiGD169KC3VUhV3eO1ra0tfNMQs9r31FNPccwxx1S6jNS88847HHfccTz22GOphW9PO3bsYMSIEQDceOONdHV1ccstt+yzXKF9K6kjIgqO53RL3sysDytXruTjH/84V111VckCHuAnP/nJnpO3Vq9ezfXXX5/K+/rAq5lZH04//XRefPHFkm/nggsu4IILLkj9fd2SN7OSqKau4KwYyj51yJtZ6pqamti2bZuDPkXd15NvahrcyCN315hZ6pqbm+ns7Bz0tc+tb913hhoMh7yZpa6xsXFQdy+y0nF3jZlZhjnkzcwyzCFvZpZhDnkzswxzyJuZZZhD3swswxzyZmYZ5pA3M8swh7yZWYY55M3MMswhb2aWYQ55M7MMc8ibmWWYQ97MLMMc8mZmGeaQNzPLMIe8mVmGOeTNzDLMIW9mlmEOeTOzDHPIm5llmEPezCzDHPJmZhnmkDczyzCHvJlZhjnkzcwyzCFvZpZhDnkzswxLJeQlXSMpJB2UN+3PJD0v6RlJZ6WxHTMzG5zhxb6BpAnAGcCLedMmAxcCU4BDgJWSjoqI3cVuz8zMBi6NlvzNwFeByJt2LrAsInZGxG+A54GTUtiWmZkNQlEteUnzgJcjYr2k/FmHAmvyXncm0wq9x0JgIcDEiROLKWfI5tyymie7tu8zffL4UTx49awKVGRmlo5+Q17SSuDgArMWAV8Dziy0WoFpUWAaEXEHcAdAW1tbwWVKbfrEMTy35W127f5g840NYvphB1SiHDOz1PQb8hFxeqHpko4DJgHdrfhm4DFJJ5FruU/IW7wZeKXoakukfXYr93Z0kv93qEGiffaRlSvKzCwFQ+6Tj4jHI2JcRLRERAu5YJ8eEb8FHgAulLSfpElAK/DLVCougXGjmjh/RjONDbkfII0N4rNtExg3sqnClZmZFack4+Qj4gngHuBJ4J+BK6p9ZE377FaGJccV3Io3s6xILeSTFv1rea8XR8QREXF0RDyU1nZKpbs1L+FWvJllRtHj5LOkfXYrz27Z4Va8mWWGQz7PuFFN3HPZzEqXYWaWGl+7xswswxzyZmYZ5pA3M8swh7yZWYY55M3MMswhb2aWYQ55M7MMc8ibmWWYQ97MLMMc8mZmGeaQNzPLMIe8mVmGOeTNzDLMIW9mlmEOeTOzDHPIm5llmEPezCzDHPJmZhnmkDczyzCHvJlZhjnkzcwyzCFvZpZhDnkzswxzyJuZZZhD3swswxzyZmYZ5pA3M8swh7yZWYYNr3QBWTLnltU82bV9n+mTx4/iwatnVaAiM6t3bsmnaPrEMTQ2aK9pjQ1i+mEHVKgiM6t3DvkUtc9uZZj2DvkGifbZR1aoIjOrdw75FI0b1cT5M5r3tOYbG8Rn2yYwbmRThSszs3rlkE9ZfmverXgzqzSHfMq6W/MSbsWbWcV5dE0JtM9u5dktO9yKN7OKc8iXwLhRTdxz2cxKl2Fm5u4aM7MsSyXkJV0jKSQdlLxukfSupHXJ4/Y0tmNmZoNTdHeNpAnAGcCLPWZtiohpxb6/mZkNXRot+ZuBrwKRwnuZmVmKigp5SfOAlyNifYHZkyT9StLPJfV64RZJCyWtlbR269atxZRjZmY99NtdI2klcHCBWYuArwFnFpjXBUyMiG2SZgA/lDQlIva5eldE3AHcAdDW1uZfA2ZmKeo35CPi9ELTJR0HTALWK3eGZzPwmKSTIuK3wM5k/Q5Jm4CjgLVpFW5mZv0b8oHXiHgcGNf9WtJmoC0iXpM0Fng9InZLOhxoBX5dbLFmZjY4pToZ6hTgLyS9D+wGLo+I10u0LTMz60VqIR8RLXnP7wPuS+u9zcxsaHzGq5lZhjnkzcwyzCFvZpZhDnkzswxzyJuZZZhD3swswxzyZmYZ5pA3M8swh7yZWYY55M3MMswhb2aWYQ55M7MMc8ibmWWYQ97MLMMc8mZmGeaQNzPLMIe8mVmGler2f1aEObes5smu7ftMnzx+FA9ePasCFZlZrXJLvgpNnziGxgbtNa2xQUw/7IAKVWRmtcohX4XaZ7cyTHuHfINE++wjK1SRmdUqd9dUoXGjmjh/RjN3r32JXbuDxgbx2bYJjBvZVOnSzCxF5eiadUu+SuW35t2KN8umcnTNOuSrVHdrXsKteLOMKkfXrEO+irXPbuXElgPdijfLqO7GXHdrvhRdsw75KjZuVBP3XDbTrXizDCt116xD3sysgkrdNevRNWZmFdY+u5Vnt+woSdesQ97MrMK6u2ZLwd01ZmYZ5pZ8hvkaOGbmlnyG+Ro4ZuaQzzBfA8fMHPIZVo4TLcysujnkM87XwDGrbw75jPM1cMzqm0fX1IFSnmhhZtXNIV8HSnmihZlVN3fXmJllWFEhL+kGSS9LWpc85uTN+zNJz0t6RtJZxZdqZmaDlUZ3zc0RsSR/gqTJwIXAFOAQYKWkoyJidwrbMzOzASpVd825wLKI2BkRvwGeB04q0bbMzKwXaYT8lZI2SLpLUvf58ocCL+Ut05lM24ekhZLWSlq7devWFMoxM7Nu/Ya8pJWSNhZ4nAv8LXAEMA3oAr7VvVqBt4pC7x8Rd0REW0S0jR07dmifwszMCuq3Tz4iTh/IG0n6LrAiedkJTMib3Qy8MujqzMysKMWOrhmf9/IzwMbk+QPAhZL2kzQJaAV+Wcy2zMxs8IodXXOTpGnkumI2A5cBRMQTku4BngTeB67wyJra5evSm9WuokI+Iub3MW8xsLiY97fqMH3iGJ7b8ja7dn9wWMXXpTerDT7j1frl69Kb1S6HvPXL16U3q10OeRsQX5ferDb5KpQ2IN2t+aW/fNGteKt7tTQYwS15G7D22a2c2HKgW/FW96ZPHLOn+7JbtQ5GcMjbgHVfl96teKt3tTQYwSFvZjZItTQYwSFvZjYEtTIYwSFvZjYE3a15iaptxYNH15iZDVn77Fae3bKjalvx4JA3Mxuy7sEI1czdNWZmGeaWvJVNLZ1AYpYVbslb2dTSCSRmWeGQt7KppRNIzLLCIW9lU0snkJhlhUPeyqpWTiAxywqHvJVVrZxAYpYVHl1jZVcLJ5CYZYVD3squFk4gsWyqx2G87q4xs7pRj8N4HfJmVjfqcRivQ97M6kY9DuN1yJtZXam3Ybw+8Go1px4Pnll66u2m9G7JW82px4Nnlq56uim9Q95qTj0ePLN01dNN6R3yVnPq8eCZ2VA55K0m1dvBM7OhcshbTfI1cMwGxqNrrGb5Gjhm/XPIW83yNXDqh4fNDp27a8ys6nnY7NA55M2s6nnY7NA55M2s6nnY7NA55M2sJnjY7NA45M2sJnjY7NB4dI2Z1QwPmx28olrykm6Q9LKkdcljTjK9RdK7edNvT6dcM6tn9XTNmbSk0ZK/OSKWFJi+KSKmpfD+ZmY2RO6uMSuST9SxapbGgdcrJW2QdJek/DMTJkn6laSfS/I33TLLJ+pYNVNE9L2AtBI4uMCsRcAa4DUggL8ExkfEpZL2A0ZExDZJM4AfAlMiYp/mjqSFwEKAiRMnznjhhReK+Dhm5bdl+3vMuunf2Pn+7/dMaxo+jH+/9g8G1XdcTb8IqqkW65+kjohoKzSv3+6aiDh9gBv5LrAiWWcnsDN53iFpE3AUsLbA+98B3AHQ1tbW918csxSlFWTdQ/vuXvsSu3bHkE/UmT5xDM9teZtduz/432CwvwjS+kxp1GLVodjRNePzXn4G2JhMHyupIXl+ONAK/LqYbZmlLc1uljRO1Enj1P20PpMvI5AdxfbJ3yTpcUkbgD8A/msy/RRgg6T1wHLg8oh4vchtmaUqzSBL40SdNE7dT+sz+TIC2VFUyEfE/Ig4LiKOj4h5EdGVTL8vIqZExNSImB4RP06nXLP0pB1kadwcuthfBGl+Jl9GIBt8WQOra2kGWRon6qTxiyCtz+TLCGSDx8lbXesOsqW/fLFqgqzYU/fT/EzF1uJROpXnkLe6V23XQ0njjldpfaZia/Eoncrrd5x8ObW1tcXatfuMsjSzGpXWOQTWt77GybtP3sxKxqN0Ks8hb2Yl5VE6leWQN7OS8iidyvKBVzMruWo7uF1PHPJmVnJpjBiyoXF3jZlZhjnkzcwyzCFvZpZhDnkzswxzyJuZZVhVXdZA0laglPf/O4jc7QpriWsuD9dcerVWL9ROzYdFxNhCM6oq5EtN0treru9QrVxzebjm0qu1eqE2a+7J3TVmZhnmkDczy7B6C/k7Kl3AELjm8nDNpVdr9UJt1ryXuuqTNzOrN/XWkjczqysOeTOzDMtcyEuaIOnfJD0l6QlJVxdY5jRJb0lalzy+Xolae9S0WdLjST373ANROd+R9LykDZKmV6LOvHqOztt/6yRtl/TlHstUfD9LukvSFkkb86YdKOmnkp5L/lvwhqOSzpb0TLLPr6tgvX8l6enk3/1+SWN6WbfP71CZa75B0st5//Zzelm37Pu4j5rvzqt3s6R1vaxbkf08ZBGRqQcwHpiePB8JPAtM7rHMacCKStfao6bNwEF9zJ8DPAQI+ATwH5WuOa+2BuC35E7IqKr9DJwCTAc25k27CbgueX4d8M1ePtMm4HDgQ8D6nt+jMtZ7JjA8ef7NQvUO5DtU5ppvAK4ZwPem7Pu4t5p7zP8W8PVq2s9DfWSuJR8RXRHxWPL8beAp4NDKVpWKc4F/iJw1wBhJ4ytdVGI2sCkiSnm28pBExL8Dr/eYfC7w98nzvwc+XWDVk4DnI+LXEfH/gGXJeiVVqN6I+JeIeD95uQZoLnUdg9HLPh6Iiuxj6LtmSQI+B/zfctRSapkL+XySWoATgP8oMHumpPWSHpI0pbyVFRTAv0jqkLSwwPxDgZfyXndSPX+8LqT3/yGqbT8DfCwiuiDXKADGFVimWvf3peR+0RXS33eo3K5Mupju6qVLrFr38Szg1Yh4rpf51baf+5TZkJc0ArgP+HJEbO8x+zFyXQtTgVuBH5a5vEL+U0RMBz4FXCHplB7zVWCdio9/lfQhYB5wb4HZ1bifB6rq9rekRcD7wNJeFunvO1ROfwscAUwDush1f/RUdfs48Xn6bsVX037uVyZDXlIjuYBfGhE/6Dk/IrZHxI7k+YNAo6SDylxmz5peSf67Bbif3E/ZfJ3AhLzXzcAr5amuT58CHouIV3vOqMb9nHi1u6sr+e+WAstU1f6WdDEwF7goko7hngbwHSqbiHg1InZHxO+B7/ZSS1XtYwBJw4E/Au7ubZlq2s8DkbmQT/rT7gSeioj/3csyByfLIekkcvthW/mq3Keej0ga2f2c3IG2jT0WewD4k2SUzSeAt7q7HCqs11ZPte3nPA8AFyfPLwZ+VGCZR4FWSZOSXysXJuuVnaSzgWuBeRHxTi/LDOQ7VDY9jhd9ppdaqmYf5zkdeDoiOgvNrLb9PCCVPvKb9gP4JLmffBuAdcljDnA5cHmyzJXAE+SO5q8BTq5wzYcntaxP6lqUTM+vWcBt5EYjPA60VcG+/jC50B6dN62q9jO5P0BdwC5yLcc/BT4K/Ax4LvnvgcmyhwAP5q07h9zorE3d/yYVqvd5cn3X3d/n23vW29t3qII1/2PyPd1ALrjHV8s+7q3mZPr/6f7+5i1bFft5qA9f1sDMLMMy111jZmYfcMibmWWYQ97MLMMc8mZmGeaQNzPLMIe8mVmGOeTNzDLMIW82AMmZjn+dnG1sVjMc8mYDczmwH7kzqs1qhkPebGDOJnf6/boK12E2KA55s35IaiJ3F6PpwM8rXI7ZoDjkzfrXSi7kn46IXZUuxmwwhle6ALMaMBY4ijLdms4sTW7Jm/XvEHI3oRnWy23szKqWQ96sD8mdgqYDBwO3A7srW5HZ4Ph68mZmGeaWvJlZhjnkzcwyzCFvZpZhDnkzswxzyJuZZZhD3swswxzyZmYZ9v8Blkkdi8O+JMcAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAEYCAYAAAC5q4bCAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAex0lEQVR4nO3df5iVZb3v8feHAUW3iPiDFBAZZTDFZILRI56gUkgyCa1DWEoitdHyF2Zlbbw6Xu7Y5a8Tkuy8bEPYzkQ0SVNrb6hj2CWkM4ZtSD1JWA4SIinmNhFmvueP9TAthvm91qwfz/q8rmtd88z9/PrOM2u+c6/7ee77VkRgZmbp1KfYAZiZWe9xkjczSzEneTOzFHOSNzNLMSd5M7MUc5I3M0sxJ3mzPJH0T5L+rZ11EyQ9n8OxQ9LInkdnlcpJ3sqKpBpJb0v6QQfbzJLUJOlNSW9IekbSOb0dW0T8S0R8tp11j0fE8b0dg1lrTvJWFJL69nDXRcBTXdhuTUQcBBwC/CuwTNIhPTynWdlykreCkfSipGsl/Rb47+4meknnA68DP+/qPhHRDPw78A9ATXKc/SXdIulPkrZKukPSAcm6D0hqlHSNpFckbZF0cbLulGT7lrglfVzSumT5+vY+Yew5bqtr8UVJv5W0Q9K9kvpnrf9Scu6XJc1udayO4n9U0q1Z294raUlXr5elj5O8FdongY8Ah0TEbkkPS3q9ndfDe3aSdDBwA3BNd04mqQq4GNgF/DEpvhEYBdQCI4GhwNeydjsSGJiUfwZYJGlQRDwFbAcmZ217IZl/Ij3xCWAKUA2cDMxKYp4CfDE5Tw0wqdV+HcU/G5gp6QxJFwCnAFf1MD5LgZ5+ZDbrqYUR8dKebyKiq23l/wwsjoiXJHVl+9MkvU6mBr8buDAiXlFm538ETo6IvwBI+hfgh8BXk313ATdExG7gUUlvAscDa4G7yCT2n0o6FDgL+HwXf4bWFkbEy0kMPyGTtCGT/L8XEeuTddeT+edIZ/FHxJ8lXZrEeQBwbkT8tYfxWQo4yVuhvdT5JnuTVEumNvvebuy2NiLeJ+kgYDEwAVgOHAEcCDRk/bMQUJW17/Ykwe/xFnBQsvwD4NnkuJ8AHo+ILd37iVr8udU5hiTLQ4CGrHV/zFruSvwPA7cDz0fEr3oYm6WEk7wV2l7Dnkr6KZkE3JbHI+LDwAeAEcCfksR2EFAl6cSIGNvhySLelPR5YGPSNv0M8DdgdERs7nbwEZslrQHOA2YC3+nuMbpgC3B01vfDs5ZfpfP45wPPAtWSPhkR9/RCjFYm3CZvRRURH46Ig9p5fTjZ7E7gODLNGbXAHcAjZJpKunKO7cC/AV9LbsR+F/iWpMEAkoZK6tKxEt8Hvgy8B1jRjf26ajkwS9KJkg4E/veeFZ3FL2kimXsQn05e35Y0tBditDLhJG8lLyLeiog/73kBbwJvR8S2bhxmAXC2pJOBa4EXgLWS3gBWkWlz76oVwDHAioj4727s1yUR8VMy8f4iifMXrTZpM/7k5vT3gcsjYnPSVLMY+J66eCPD0keeNMSs+yRtBC6JiFXFjsWsI67Jm3WTpI+TubfQuoZtVnJ849WsGyQ9BpwIzEzax81KmptrzMxSzM01ZmYpVlLNNYcffniMGDGi2GGYmZWVhoaGVyPiiLbWlVSSHzFiBPX19cUOw8ysrEj6Y3vr3FxjZpZiTvJmZinmJG9mlmIl1SZvZqVv165dNDY28vbbbxc7lIrTv39/hg0bRr9+/bq8j5O8mXVLY2MjAwYMYMSIEXhInMKJCLZv305jYyPV1dVd3q/sk/ySG69i5Y6hrGke3VI2vs8GJg/czOxrbytiZGbp9PbbbzvBF4EkDjvsMLZt6864fClok48hY1nUbyHj+2wAMgl+Ub+FxJAOhxk3sxw4wRdHT6572dfkp06bwdybt3F7vwX8oGkSF1at4urmudx87oxih2ZmVnRlX5MffHB/ho89i3uaJ3NV3xXc0zyZo8dNYfCA/p3vbGZlqaqqitraWk466SSmTp3K66+/ntfj19fXc+WVV+5T/thjjyGJn/zkJy1l55xzDo899lhez59PZV+TB7imZiuxbiW37T6PmVUrUc0FwEnFDsus4p192+P8bssb+5SfeNTBPHpVe7M+du6AAw5g3bp1AFx00UUsWrSIefPm9fh4rdXV1VFXV9fmumHDhjF//nymTp2at/P1prKvybNpNYMencOKkfNZ0DSdFSPnM+iRObBpdbEjM6t4Y4cfQr+qvduR+1WJsccMyts5xo8fz+bNmeluN27cyJQpUxg3bhwTJkzgueeeA2DWrFlceeWVnH766Rx77LHcf//9AMycOZMHH3yw5VgXXHABDz30EI899hjnnHNOm+cbM2YMAwcOZOXKlfusa2ho4P3vfz/jxo3jrLPOYsuWzBzvTz31FCeffDLjx4/nS1/6EiedVLhKaPkn+c1Pw/SlTJ02g1NGHMrUc2fA9KWZcjMrqivPrKFPq5uFVRJXnjkyL8dvamri5z//OR/96EcBmDNnDt/+9rdpaGjglltu4fOf/3zLtlu2bOFXv/oVDz/8MF/5ylcA+OxnP8v3vvc9AHbs2METTzzB2Wef3el5r7vuOr7+9a/vVbZr1y6uuOIK7r//fhoaGpg9e3bLp4uLL76YO+64gzVr1lBVVZWXn72ryr+55n1zARgMLL9kfKZswESonli0kMwsY/DB/Zk+bhj31r/ErqagX5X4X3VH53zP7G9/+xu1tbW8+OKLjBs3jsmTJ/Pmm2/yxBNPMH369Jbtdu7c2bJ87rnn0qdPH0488US2bt0KwPvf/34uu+wyXnnlFR544AE+/vGP07dv52lxwoRMU9Pjjz/eUvb888+zfv16Jk+eDGT+AR111FG8/vrr/PWvf+X0008H4FOf+hQPP/xwTj9/d5R/kjezknblmTXc19AIRN5q8Xva5Hfs2ME555zDokWLmDVrFoccckhLW31r+++/f8ty9mRJM2fO5O6772bZsmUsWbKkyzHMmzeP+fPnt/xTiAhGjx7NmjVr9trutdde68ZPln/l31xjZiVtT21eIi+1+GwDBw5k4cKF3HLLLRxwwAFUV1dz3333AZmk+8wzz3R6jFmzZrFgwQIARo8e3fHGWT70oQ/x2muvtZzj+OOPZ9u2bS1JfteuXWzYsIFBgwYxYMAA1q5dC8CyZcu68yPmzEnezHrdlWfWcMqIQ/PWFp/tve99L2PGjGHZsmXcfffdLF68mDFjxjB69Oi9bqq2513vehcnnHACF198cbfPPW/ePBobGwHYb7/9uP/++7n22msZM2YMtbW1PPHEEwAsXryYOXPmMH78eCKCgQMHdvtcPVVSc7zW1dWFJw0xK23PPvssJ5xwQrHDyJu33nqL97znPTz99NO9lnzffPNNDjroIAC++c1vsmXLFm67rWfDrrR1/SU1RESbz3y6Jm9mFWvVqlW8+93v5oorrujV2vUjjzzS0nnr8ccf57rrruu1c7XmG69mVrEmTZrEn/70p14/z4wZM5gxozhDrbgmb2aWYk7yZmYp5iRvZpZiTvJmZinmJG9mZSeNQw2vW7eO8ePHM3r0aE4++WTuvffenI8JTvJm1pt+tWDfEWE3rc6U52DPsAbr16/n0EMPZdGiRTkdr7W6ujoWLlzY5ro9Qw3n24EHHsj3v/99NmzYwM9+9jPmzp2bl39eOSV5SddL2ixpXfI6O2vdVyW9IOl5SWflHKmZlZ+hY+G+WX9P9JtWZ74fmr/pOdMy1PCoUaOoqakBYMiQIQwePLjb87m2JR81+W9FRG3yehRA0onA+cBoYArwr5IKO76mmRVf9cTM0N/3zYJfzM98nb40b6PEpnWo4SeffJJ33nmH4447rmsXogO91RlqGrAsInYCmyS9AJwKrOl4NzNLneqJUPcZWH0TTPxyXhJ8moca3rJlCzNnzuSuu+6iT5/c6+H5SPKXS/o0UA9cExGvAUOBtVnbNCZl+5A0B5gDMHz48DyEY2YlZdNqqF+cSfD1i6F6Qs6JPq1DDb/xxht85CMf4etf/zqnnXZal/frSKf/JiStkrS+jdc04DvAcUAtsAW4dc9ubRyqzZHQIuLOiKiLiLojjjiiZz+FmZWmPW3w05fCGfP+3nSTp+k50zTU8DvvvMN5553Hpz/96b0+jeSq0yQfEZMi4qQ2Xg9GxNaIaIqIZuC7ZJpkIFNzPzrrMMOAl/MWtZmVh2R6zpaa+542+jxOz5mWoYaXL1/O6tWrWbp0KbW1tdTW1rb7qaQ7chpqWNJREbElWb4a+B8Rcb6k0cAPyST9IcDPgZqIaOroeB5q2Kz0eajh7ivmUMO5tsnfJKmWTFPMi8AlABGxQdJy4HfAbuCyzhK8mVmhrVq1itmzZ/OFL3yh14ca/sY3vsHu3bs55phjWLp0aa+dqzVPGmJm3ZK2mny58aQhZtbrSqlyWEl6ct2d5M2sW/r378/27dud6AssIti+fTv9+3dvInTPDGVm3TJs2DAaGxvz0uXeuqd///4MGzasW/s4yZtZt/Tr14/q6upih2Fd5OYaM7MUc5I3M0sxJ3kzsxRzkjczSzEneTOzFHOSNzNLMSd5M7MUc5I3M0sxJ3kzsxRzkjczSzEneTOzFHOSNzNLMSd5M7MUc5I3M0sxJ3kzsxRzkjczSzEneTOzFPPMUMCSG69i5Y6hrGke3VI2vs8GJg/czOxrbytiZGZmuXFNHoghY1nUbyHj+2wAMgl+Ub+FxJCxRY7MzCw3rskDU6fNYO7N27i93wJ+0DSJC6tWcXXzXG4+d0axQzMzy4lr8sDgg/szfOxZ3NM8mav6ruCe5skcPW4Kgwf0L3ZoZmY5cU0+cU3NVmLdSm7bfR4zq1aimguAk4odlplZTlyTB9i0mkGPzmHFyPksaJrOipHzGfTIHNi0utiRmZnlxDV5gM1Pw/SlTD3sVP7jrd8w9dwz4dXqTHn1xGJHZ2bWY4qInu8sXQ/8I7AtKfqniHhU0gjgWeD5pHxtRFza2fHq6uqivr6+x/GYmVUiSQ0RUdfWunzU5L8VEbe0Ub4xImrzcHwzM+sht8mbmaVYPpL85ZJ+K2mJpEFZ5dWSfiPpl5ImtLezpDmS6iXVb9u2rb3NzMysBzptk5e0CjiyjVXzgLXAq0AA/wwcFRGzJe0PHBQR2yWNA34MjI6INzo6l9vkzcy6L6c2+YiY1MWTfBd4ONlnJ7AzWW6QtBEYBTiDm5kVUE7NNZKOyvr2PGB9Un6EpKpk+VigBvhDLucyM7Puy/Xpmpsk1ZJprnkRuCQpnwjcIGk30ARcGhF/yfFcZmbWTTkl+YiY2U75j4Af5XJsMzPLnR+hNDNLMSd5M7MUc5I3M0sxJ3kzsxRzkjczSzEneTOzFHOSNzNLMSd5M7MUc5I3M0sxJ3kzsxRzkjczSzEneTOzFHOSNzNLMSd5M7MUc5I3M0sxJ3kzsxRzkjczSzEneTOzFHOSNzNLsVwn8rYsS268ipU7hrKmeXRL2fg+G5g8cDOzr72tiJGZWaVyTT6PYshYFvVbyPg+G4BMgl/UbyExZGyRIzOzSuWafB5NnTaDuTdv4/Z+C/hB0yQurFrF1c1zufncGcUOzcwqlGvyeTT44P4MH3sW9zRP5qq+K7ineTJHj5vC4AH9ix2amVUo1+Tz7JqarcS6ldy2+zxmVq1ENRcAJxU7LDOrUK7J59Om1Qx6dA4rRs5nQdN0Voycz6BH5sCm1cWOzMwqlGvy+bT5aZi+lKmHncp/vPUbpp57JrxanSmvnljs6MysAikiih1Di7q6uqivry92GGZmZUVSQ0TUtbXOzTVmZimWc5KXdIWk5yVtkHRTVvlXJb2QrDsr1/OYmVn35dQmL+mDwDTg5IjYKWlwUn4icD4wGhgCrJI0KiKacg3YzCwtCtFLPtea/OeAb0bEToCIeCUpnwYsi4idEbEJeAE4NcdzmZmlSiF6yeea5EcBEyT9WtIvJZ2SlA8FXsrarjEp24ekOZLqJdVv27Ytx3DMzMrH1GkzmNs8l9v7LeTqvvdlvjbPZWoee8l32lwjaRVwZBur5iX7DwJOA04Blks6FlAb27f5GE9E3AncCZmna7oWtplZ+WvpJf/0Bq7q+wC3N30s773kO03yETGpvXWSPgc8EJnnMJ+U1AwcTqbmfnTWpsOAl3OM1cwsdXq7l3yuzTU/Bs4AkDQK2A94FXgIOF/S/pKqgRrgyRzPZWaWLgXoJZ9rj9clwBJJ64F3gIuSWv0GScuB3wG7gcv8ZI2ZWSsF6CXvHq9mZmXOPV7NzCqUk7yZWYo5yZuZpZiTvJlZijnJm5mlmJO8mVmKOcmbmaWYk7yZWYo5yZuZpZiTvJlZijnJm5mlmJO8mVmK5ToKpfWCQsz7aGaVwTX5ElSIeR/NrDK4Jl+Cpk6bwdybt3F7vwX8oGkSF1at4urmudycx3kfzawyuCZfglrmfWyezFV9V3BP8+S8z/toZpXBNfkS1dvzPppZZXBNvhQVYN5HM6sMrsmXogLM+2hmlcFzvJqZlTnP8WpmVqGc5M3MUsxJ3swsxZzkzcxSzEnezCzFnOTNzFLMSd7MLMVyTvKSrpD0vKQNkm5KykZI+pukdcnrjtxDNTOz7sqpx6ukDwLTgJMjYqekwVmrN0ZEbS7HNzOz3ORak/8c8M2I2AkQEa/kHpKZmeVLrmPXjAImSJoPvA18MSKeStZVS/oN8AZwXUQ83tYBJM0B5gAMHz48x3Asm2eYMusd5fS31WlNXtIqSevbeE0j809iEHAa8CVguSQBW4DhEfFe4AvADyUd3NbxI+LOiKiLiLojjjgibz+YeYYps95STn9bndbkI2JSe+skfQ54IDKjnD0pqRk4PCK2AXuacBokbSRT6/foYwXkGabMekc5/W3l2ib/Y+AMAEmjgP2AVyUdIakqKT8WqAH+kOO5rJs8w5RZ7yinv61c2+SXAEskrQfeAS6KiJA0EbhB0m6gCbg0Iv6S47msBzzDlFnvKJe/rZySfES8A1zYRvmPgB/lcmzLg2SGqcUj57Pgd4dz0PEf5DOPzIEDl3ryEbNclNHflmeGSjPPMGXWO8rob8szQ5mZlTnPDGVmVqGc5M3MUsxJ3swsxZzkzcxSzEnezCzFnOTNzFLMSd7MLMWc5M3MUsw9Xq1T5TR2tpntzTV561Q5jZ1tZntzTd46VU5jZ5vZ3lyTt06V09jZZrY31+StS8pl7Gwz25uTvHWujMbONutIJT5E4OYa69yesbOnzeCUEYcy9dwZMH1pptysjFTiQwQeT97MKsYrb7zNF2++nW/1+ftDBF9onsvNX768rO8xeTx5MzMq8yECt8mbWUWptIcIXJM3s8qRPESwYuR8FjRNZ8XI+Qx6ZA5sWl3syHqNa/JmVjnKaALufPGNVzOzMucbr2ZmFcrNNVYwldgRxazYXJO3gqnEjihmxeaavBWMR7M0KzzX5K1gKrEjilmx5ZTkJd0raV3yelHSuqx1X5X0gqTnJZ2Vc6SWCtfUbOVTfTIdUT7VZyXX1Py52CGZpVpOzTUR0fI5W9KtwI5k+UTgfGA0MARYJWlURDTlcj4rcx7N0qzg8tImL0nAJ4AzkqJpwLKI2AlskvQCcCqwJh/nszJVgR1RLD/8ZFbP5atNfgKwNSJ+n3w/FHgpa31jUrYPSXMk1Uuq37ZtW57CsZL0vrlQPZHBB/dn+SXjM23x1RMz5WYd8JNZPddpTV7SKuDINlbNi4gHk+VPAvdk79bG9m12rY2IO4E7IdPjtbN4zKzy+Mmsnus0yUfEpI7WS+oLfAwYl1XcCByd9f0w4OWeBGhm1vJk1tMbuKrvA9ze9DE/mdVF+WiumQQ8FxGNWWUPAedL2l9SNVADPJmHc5lZhfKTWT2Tjxuv57N3Uw0RsUHScuB3wG7gMj9ZY2Y95iezeiznJB8Rs9opnw/Mz/X4ZmZ+MqvnPNSwlR0/Tme2Nw81bKnix+nMus4DlFnZ8eN0Zl3nmryVHQ90ZtZ1rslbWbqmZiuxLvM43cyqlajmAuCkYodlVnJck7fykzxOt2LkfBY0TWfFyPkMemQObFpd7MjMSo5r8lZ+/DidWZf5EUozszLX0SOUrsmbWa9xn4bic5I3s33kKznHkLEsemsel+26kjXNo1v6NKwY4s7wheIbr2a2j3x1OJs6bQZzm+dye7+FXN33vszX5rlMdZ+GgnGSN7N95Cs5u09D8bm5xixF8tXMks/x292nobic5M1ylK/Emo/j5LMNPC/J2UMEF52ba8xylK/263wcJ29t4PnqcLanT8O0GZwy4tBMHNOXZsqtIFyTN8tRvgZMy8dx8tbMkq8OZ8kk7YOB5ZeMz5QNmOhafAE5yZvlKF+JNV/HyUszi5NzajjJm+VBvm4u5nwct4FbK07yZrnKV2LNx3E8ro+14iRvlqt8JdZ8HMfNLNaKBygzMytznuPVzKxCubnGKpZHSLRK4Jq8Vax8dWIyK2WuyVvFylcnJrNS5pq8VSyPkGiVwDV5q2geIdHSzjV5q1z5GoTLrITlVJOXdC9wfPLtIcDrEVEraQTwLPB8sm5tRFyay7nM8s69Q60C5K0zlKRbgR0RcUOS5B+OiG597nVnKDOz7uuoM1Re2uQlCfgEcEY+jmdmZvmRrzb5CcDWiPh9Vlm1pN9I+qWkCe3tKGmOpHpJ9du2bctTOGZmBl2oyUtaBRzZxqp5EfFgsvxJ4J6sdVuA4RGxXdI44MeSRkfEG60PEhF3AndCprmmuz+AmZm1r9MkHxGTOlovqS/wMWBc1j47gZ3JcoOkjcAowA3uZmYFlI/mmknAcxHRuKdA0hGSqpLlY4Ea4A95OJeZmXVDPm68ns/eTTUAE4EbJO0GmoBLI+IvnR2ooaHhVUl/zENM7TkceLUXj98bHHNhOObeV27xQvnEfEx7K0pqPPneJqm+vceMSpVjLgzH3PvKLV4oz5hbc49XM7MUc5I3M0uxSkvydxY7gB5wzIXhmHtfucUL5RnzXiqqTd7MrNJUWk3ezKyiOMmbmaVY6pK8pKMl/V9Jz0raIOmqNrb5gKQdktYlr68VI9ZWMb0o6b+SePbpGayMhZJekPRbSUWdiFTS8VnXb52kNyTNbbVN0a+zpCWSXpG0PqvsUEkrJf0++TqonX2nSHo+ueZfKWK8N0t6Lvm9r5B0SDv7dvgeKnDM10vanPW7P7udfQt+jTuI+d6seF+UtK6dfYtynXssIlL1Ao4CxibLA4D/B5zYapsPkBkKuejxZsX0InB4B+vPBn4KCDgN+HWxY86KrQr4M3BMqV1nMh3zxgLrs8puAr6SLH8FuLGdn2kjcCywH/BM6/dRAeP9ENA3Wb6xrXi78h4qcMzXA1/swvum4Ne4vZhbrb8V+FopXeeevlJXk4+ILRHxdLL8VzKTlwwtblR5MQ34fmSsBQ6RdFSxg0qcCWyMiN7srdwjEbEaaN3behpwV7J8F3BuG7ueCrwQEX+IiHeAZcl+vaqteCPiPyNid/LtWmBYb8fRHe1c464oyjWGjmPOGjq9dU/+spS6JJ8tmbzkvcCv21g9XtIzkn4qaXRhI2tTAP8pqUHSnDbWDwVeyvq+kdL559XW0BZ7lNp1BnhXRGyBTKUAGNzGNqV6vWeT+UTXls7eQ4V2edLEtKSdJrFSvcZtDZ2erdSuc4dSm+QlHQT8CJgb+w5x/DSZpoUxwLeBHxc4vLb8z4gYC3wYuExS6/nn1MY+RX/+VdJ+wEeB+9pYXYrXuatK7npLmgfsBu5uZ5PO3kOF9B3gOKCWzNDjt7axTcld40TrodNbK6Xr3KlUJnlJ/cgk+Lsj4oHW6yPijYh4M1l+FOgn6fACh9k6ppeTr68AK8h8lM3WCByd9f0w4OXCRNehDwNPR8TW1itK8Tontu5p6kq+vtLGNiV1vSVdBJwDXBBJw3BrXXgPFUxEbI2IpohoBr7bTiwldY1hr6HT721vm1K6zl2RuiSftKctBp6NiP/TzjZHJtsh6VQy12F74aLcJ55/kDRgzzKZG23rW232EPDp5Cmb08jMp7ulwKG2pd1aT6ld5ywPARclyxcBD7axzVNAjaTq5NPK+cl+BSdpCnAt8NGIeKudbbryHiqYVveLzmsnlpK5xln2GTo9W6ld5y4p9p3ffL+A95H5yPdbYF3yOhu4lMyQxwCXAxvI3M1fC5xe5JiPTWJ5JolrXlKeHbOARWSeRvgvoK4ErvWBZJL2wKyykrrOZP4BbQF2kak5fgY4DPg58Pvk66HJtkOAR7P2PZvM01kb9/xOihTvC2Tarve8n+9oHW9776Eixvzvyfv0t2QS91Glco3bizkpX7rn/Zu1bUlc556+PKyBmVmKpa65xszM/s5J3swsxZzkzcxSzEnezCzFnOTNzFLMSd7MLMWc5M3MUsxJ3qwLkp6Otye9jc3KhpO8WddcCuxPpke1Wdlwkjfrmilkut+vK3IcZt3iJG/WCUn9ycxiNBb4ZZHDMesWJ3mzztWQSfLPRcSuYgdj1h19ix2AWRk4AhhFgaamM8sn1+TNOjeEzCQ0fdqZxs6sZDnJm3UgmSloLHAkcAfQVNyIzLrH48mbmaWYa/JmZinmJG9mlmJO8mZmKeYkb2aWYk7yZmYp5iRvZpZiTvJmZin2/wF+CWYRwX4QjgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "ent_plotter(r=2)\n",
    "ent_plotter(r=3)\n",
    "ent_plotter(r=4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "No handles with labels found to put in legend.\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAELCAYAAAAx94awAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAAAQZElEQVR4nO3df4xldXnH8fcHXLrqEkV3hS1ghxpsqxaBTCnNGkNraZCkVVtLNIrUH0EbbDQhTY0xVVua+odirDbVJWyFRIwm2koNtlJiS2wqdpassGa1WgW7umEHFHaJWGF5+sc9K7uz8+POzpx75873/Upu9t5zzp3zcDj5zHfOfe73pKqQJLXlhHEXIEkaPcNfkhpk+EtSgwx/SWqQ4S9JDTL8JalBTxp3AcPYvHlzTU1NjbsMSZooO3fuvL+qtsy3biLCf2pqipmZmXGXIUkTJcm9C63zso8kNcjwl6QGGf6S1KCJuOYvSa179NFH2bt3Lz/5yU+OWbdx40bOOOMMNmzYMPTPM/wlaQLs3buXk08+mampKZL8bHlV8cADD7B3717OOuusoX+e4S9Ja8z0Nbdy/8M/PWrZ9t87jdNPegabjwh+gCQ885nPZHZ2dln78Jq/JK0xc4MfIIRDC0zBnzm/EIZh+EtSgwx/SWqQ4S9JE6AoFrrz4vHckdHwl6QJcO+Dj/LYjw8cE/SHu302bty4rJ9nt48krTGbN510zIe+H77jR1x94gk85fEfH7P94T7/5TD8JWmNmXnXxb3vw8s+ktQgw1+SGmT4S1KDDH9JapDhL0kNMvwlqUGGvyQ1yPCXpAYZ/pLUIMNfkhpk+EtSgwx/SWqQ4S9JDTL8JalBhr8kNai38E9yZpIvJdmT5OtJ3tYtf0+S7yfZ1T0u7asGSdL8+ryZy2PA1VV1Z5KTgZ1Jbu3WfbCq3t/jviVJi+gt/KtqH7Cve34wyR7g9L72J0ka3kiu+SeZAs4D7ugWvTXJXUl2JDllgfdcmWQmyczs7OwoypSkZvQe/kk2AZ8B3l5VB4C/A54DnMvgL4MPzPe+qtpeVdNVNb1ly5a+y5SkpvQa/kk2MAj+T1TVZwGq6r6qOlRVjwPXARf0WYMk6Vi9XfNPEuB6YE9VXXvE8q3d5wEArwB291WDJK1F09fcyv0P//SY5Zs3ncTMuy4eSQ19dvtsAy4H7k6yq1v2TuDVSc4FCrgHeHOPNUjSmjNf8C+2vA99dvt8Gcg8q27pa5+SpOH4DV9JapDhL0kNMvwlqUGGvySN2OZNJy1reR/67PaRJM1jVO2ci3HkL0kNMvwlqUGGvyQ1yPCXpAYZ/pLUIMNfkhpk+EtSgwx/SWqQ4S9JDTL8JalBTu8gSatsLdypaymO/CVpla2FO3UtxfCXpAYZ/pLUIMNfkhpk+EtSgwx/SVpla+FOXUux1VOSVtlaaedcjCN/SWqQ4S9JDTL8JalBhr8kNcjwl6QGGf6S1CBbPSVpmSZh1s6lOPKXpGWahFk7l2L4S1KDDH9JapDhL0kNMvwlqUGGvyQt0yTM2rmU3lo9k5wJ3AicBjwObK+qDyV5BvApYAq4B7isqn7UVx2StNompZ1zMX2O/B8Drq6qXwEuBK5K8jzgHcBtVXU2cFv3WpI0Qr2Ff1Xtq6o7u+cHgT3A6cDLgBu6zW4AXt5XDZKk+Y3kmn+SKeA84A7g1KraB4NfEMCzRlGDJOkJvYd/kk3AZ4C3V9WBZbzvyiQzSWZmZ2f7K1CSGtRr+CfZwCD4P1FVn+0W35dka7d+K7B/vvdW1faqmq6q6S1btvRZpiQ1p7fwTxLgemBPVV17xKqbgSu651cAn+urBknS/Pqc1XMbcDlwd5Jd3bJ3Au8DPp3kjcD3gD/ssQZJ0jx6C/+q+jKQBVa/pK/9StJKrYcpm5fiN3wlaY71MGXzUgx/SWqQ4S9JDTL8JalBhr8kNcjwl6Q51sOUzUtZstUzyUGg5ix+CJhhMGvnd/ooTJJW27AtnOulnXMxw/T5Xwv8ALiJQd/+qxjM0f9NYAdwUV/FSdJqaqGFc1jDXPa5pKo+VlUHq+pAVW0HLq2qTwGn9FyfJKkHw4T/40kuS3JC97jsiHVzLwdJkibAMOH/GgZz9OwH7uuevzbJk4G39libJKknS17z7z7Q/d0FVn95dcuRJI3CkiP/JM9NcluS3d3rc5K8q//SJGl1tdDCOaxhun2uA/4U+BhAVd2V5Cbgmj4Lk6TV1kIL57CGCf+nVNVXB/dm+ZnHeqpHkpathSmYV9swH/jen+Q5dJ09SV4J7Ou1KklaBvv3l2+Ykf9VwHbgl5N8H/gu8Npeq5Ik9WrYbp/fTvJU4ISqOth/WZKkPi0Y/klet8ByAKrqxp5qkiT1bLGR/6/NsywMev5PBwx/SZpQC4Z/Vf3J4ecZDPdfA/wZ8BXgr/ovTZKGs3nTSQt2+2h+i17zT/Ik4I+Aq4E7gFdW1TdHUJckDc12zuVb7Jr/VcDbgNsYzOx578iqktQ8e/f7tdjI/8MMJnN7EfBPR3zJK0BV1Tk91yapYfbu92ux8D9rZFVIkkZqsQ98vcwjSeuUN3CXpAYZ/pLUIMNf0prk3Pv9WnJunyR3c+y9eh8CZoBrquqBPgqTtP4sp33Tds5+DTOr5xeAQ8BN3etXdf8eAD7Owrd4lKSj2L65dgwT/tuqatsRr+9O8h9VtS2JUztL0gQa5pr/piS/fvhFkguATd1L7+glSRNomJH/m4AdSTYx+HbvAeCN3fz+f91ncZKkfgxzM5f/An41ydOAVNWDR6z+dF+FSVrbnHtnsi152SfJ05Jcy2CCt39N8oHuF4Gkhh3Ph7e2b64dw1z22QHsBi7rXl8O/D3w+30VJWl98i+CtWOYD3yfU1XvrqrvdI/3Ar+41JuS7EiyP8nuI5a9J8n3k+zqHpeupHhJ0vEZJvwfSfKiwy+SbAMeGeJ9HwcumWf5B6vq3O5xy3BlSpJW0zCXfd4C3HjEdf4fAVcs9aaquj3J1ApqkyT1ZMmRf1V9rapeCJwDnFNV5wG/tYJ9vjXJXd1loVMW2ijJlUlmkszMzs6uYHeS+uCHt5MtVXOn7RniTcn3qurZQ2w3BXy+ql7QvT4VuJ/BXEF/CWytqjcs9XOmp6drZmZm2XVKUsuS7Kyq6fnWHe+snll6k2NV1X1VdaiqHgeuAy44zv1LklZgmGv+81n+nwtAkq1Vta97+QoGLaSSxsgva7VpwfBPcpD5Qz7Ak5f6wUk+CVwEbE6yF3g3cFGSc7ufew/w5mVXLGlVOdNmmxa7h+/JK/nBVfXqeRZfv5KfKUlaHd7JS5IaZPhLUoMMf0lqkOEvNc4va7XpeFs9Ja0RK23VtJ2zTY78pQlnq6aOh+EvSQ0y/CWpQYa/JDXI8JekBhn+0oSzVVPHw1ZPacLZqqnjYfhLY+JUyhonL/tIY2J/vsbJ8JekBhn+ktQgw1+SGmT4S1KDDH9pTOzP1zjZ6imNie2cGidH/pLUIEf+0hD8QpbWG0f+0hD8QpbWG8Nfkhpk+EtSgwx/SWqQ4S9JDTL8pSH4hSytN7Z6SkOwnVPrjSN/SWqQ4S9JDTL8JalBhr8kNcjwl6QG2e2jieeka9LyOfLXxHPSNWn5egv/JDuS7E+y+4hlz0hya5Jvdf+e0tf+JUkL63Pk/3HgkjnL3gHcVlVnA7d1ryVJI9Zb+FfV7cAP5yx+GXBD9/wG4OV97V+StLBRX/M/tar2AXT/PmvE+5cksYY/8E1yZZKZJDOzs7PjLkdrmJOuScs36lbP+5Jsrap9SbYC+xfasKq2A9sBpqena1QFavLYzikt36hH/jcDV3TPrwA+N+L9S5Lot9Xzk8B/Ar+UZG+SNwLvAy5O8i3g4u61JGnEervsU1WvXmDVS/rap8bLb9pKk2PNfuCryeM3baXJYfhLUoMMf0lqkOEvSQ0y/CWpQYa/Vo3ftJUmhzdz0aqxnVOaHI78JalBhr8kNcjwl6QGGf6S1CDDX5IaZLfPhHHyNEmrwZH/hHHyNEmrwfCXpAYZ/pLUIMNfkhpk+EtSgwz/CePkaZJWg62eE8Z2TkmrwZG/JDXI8JekBhn+ktQgw1+SGmT4S1KDDH9JapDhL0kNMvwlqUGGvyQ1yPCXpAaty+kdvNuVJC1uXY78vduVJC1uXYa/JGlxhr8kNcjwl6QGGf6S1KB1Gf7e7UqSFjeWVs8k9wAHgUPAY1U1vZo/33ZOSVrcOPv8f7Oq7h/j/iWpWevyso8kaXHjCv8CvphkZ5Ir59sgyZVJZpLMzM7Ojrg8SVrfxhX+26rqfOClwFVJXjx3g6raXlXTVTW9ZcuW0VcoSetYqmq8BSTvAR6uqvcvss0scO/IihqPzYCfgQx4LI7m8XiCx+JoSx2PX6iqeUfPI//AN8lTgROq6mD3/HeAv1jsPQsVv54kmVntrqdJ5bE4msfjCR6Lo63keIyj2+dU4B+SHN7/TVX1z2OoQ5KaNfLwr6rvAC8c9X4lSU+w1XPt2D7uAtYQj8XRPB5P8Fgc7biPx9g/8JUkjZ4jf0lqkOEvSQ0y/EcoyY4k+5PsXmB9kvxNkm8nuSvJ+aOucVSGOBYXJXkoya7u8eejrnGUkpyZ5EtJ9iT5epK3zbNNE+fHkMeimfMjycYkX03yte54vHeebZZ/blSVjxE9gBcD5wO7F1h/KfAFIMCFwB3jrnmMx+Ii4PPjrnOEx2MrcH73/GTgv4HntXh+DHksmjk/uv/fm7rnG4A7gAtXem448h+hqrod+OEim7wMuLEGvgI8PcnW0VQ3WkMci6ZU1b6qurN7fhDYA5w+Z7Mmzo8hj0Uzuv/fD3cvN3SPuZ06yz43DP+15XTgf494vZeGT3rgN7o/db+Q5PnjLmZUkkwB5zEY4R2pufNjkWMBDZ0fSU5MsgvYD9xaVSs+Nwz/tSXzLGu1F/dOBvOSvBD4MPCP4y1nNJJsAj4DvL2qDsxdPc9b1u35scSxaOr8qKpDVXUucAZwQZIXzNlk2eeG4b+27AXOPOL1GcAPxlTLWFXVgcN/6lbVLcCGJJvHXFavkmxgEHafqKrPzrNJM+fHUseixfMDoKoeBP4NuGTOqmWfG4b/2nIz8Lruk/sLgYeqat+4ixqHJKelmwAqyQUMztUHxltVf7r/1uuBPVV17QKbNXF+DHMsWjo/kmxJ8vTu+ZOB3wa+MWezZZ8b47yNY3OSfJJBl8LmJHuBdzP48Iaq+ihwC4NP7b8N/Bh4/Xgq7d8Qx+KVwB8neQx4BHhVdW0N69Q24HLg7u7aLsA7gWdDc+fHMMeipfNjK3BDkhMZ/JL7dFV9Pslb4PjPDad3kKQGedlHkhpk+EtSgwx/SWqQ4S9JDTL8JalBhr8kNcjwl5aQ5OEhtnlqko90X7CR1jzDX1odbwF+DnjRuAuRhmH4S0NIMtXdXOS67oYaX+y+an/YJQzmnd81ngql5TH8peGdDfxtVT0feBD4AxjcaQk4kcHNaf59bNVJy2D4S8P7blXt6p7vBKa652czCP9vVNWjY6hLWjbDXxre/x3x/BBPTIy4BXgu8CGAJK9P8tJuhsUdcy4PSWuCs3pKK/fzDOaePyHJKcDtwBsY3EnpU1X1yDiLk+bjyF9agSRPYnCt/zTgo8ChqvofBrcePLeq/mWc9UkLcUpnqQdJbgbeVFX7x12LNB9H/tIqSvK0JB8BbjD4tZY58pekBjnyl6QGGf6S1CDDX5IaZPhLUoMMf0lqkOEvSQ0y/CWpQYa/JDXI8JekBv0/6Pz02H1UgxQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#plots log neg\n",
    "\n",
    "LN_ex=np.zeros(len(Lysub_list)) # Log negativity\n",
    "\n",
    "Gmat=np.dot(V1,np.matrix(V1).H)\n",
    "\n",
    "for i_L in range(len(Lysub_list)):\n",
    "    Lysub=Lysub_list[i_L]\n",
    "    inds_A = np.arange(0,int(Lx*Lysub))\n",
    "    inds_AB = np.arange(0,int(2*Lx*Lysub))\n",
    "    Nsub =  int(Lx*Lysub)\n",
    "    LN_ex[i_L]= exact_log_neg_calc(Nsub,Nsub,Gmat[np.ix_(inds_AB,inds_AB)])\n",
    "    \n",
    "plt.figure()\n",
    "plt.plot(np.log(Lysub_list),LN_ex,'s')\n",
    "# xplot=np.log(np.tan(pi*Lsub_list/N))\n",
    "# yplot=(xplot-xplot[0])/4+ 0.47\n",
    "# plt.plot(np.log(Lsub_list),yplot,label='CFT')\n",
    "# plt.legend()\n",
    "# plt.ylim([0,2])\n",
    "plt.xlabel('$\\ln \\ell_y$')\n",
    "plt.ylabel('Log Neg')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}