1 - Matrix Encoding.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Matrix Encoding and Permanent Estimation\n",
    "\n",
    "_This notebook was developed by Ana Filipa Carvalho._\n",
    "\n",
    "This notebook is part of a series of notebooks with the code to reproduce results from this [REFERENCE] paper. This first contains the initial steps namely the encoding function *to_unitary()* and permanent estimation. In order to do that, the notebook is organised as follows:\n",
    "- Random graph generation (to test our encoding function and later to be used to test algorithms for applications).\n",
    "- *to_unitary()* function\n",
    "- State preparation and post-selection (for permanent estimation).\n",
    "- *perm_estimation()* function to estimate the permanent and obtain some statistics around it (estimation with increase of edge connectivity, ...).\n",
    "\n",
    "To run these notebooks, you will need to install **Perceval**, a software tool developed by Quandela that we use to simulate the photonic device."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from numpy import linalg\n",
    "from scipy.linalg import sqrtm\n",
    "\n",
    "import networkx as nx\n",
    "\n",
    "import perceval as pcvl\n",
    "import perceval.components as comp\n",
    "from perceval.algorithm import Sampler\n",
    "\n",
    "import quandelibc as qc\n",
    "import matplotlib.pyplot as plt\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def random_graph_generation(n_nodes, p1, n_densest=0, p2=0):\n",
    "    '''Generates a random graph with:\n",
    "           .n_nodes as number of nodes and edge connectivity determined by probability p1\n",
    "           .n_densest - a subgraph of size n_densest with edge conections given by p2\n",
    "    '''\n",
    "    G1 = nx.Graph()\n",
    "    G1.add_nodes_from([k for k in range(n_nodes-n_densest)])\n",
    "\n",
    "    # adding edges with probability p1 for the given nodes\n",
    "    for i in range(n_nodes-n_densest):\n",
    "        for j in range(i + 1, n_nodes-n_densest):\n",
    "            r = np.random.random()\n",
    "            if r < p1:\n",
    "                G1.add_edge(i, j)\n",
    "\n",
    "    # adding edges with probability p2 for selected nodes with probability p2\n",
    "    G2 = nx.Graph()\n",
    "    G2.add_nodes_from([k for k in range(n_nodes-n_densest, n_nodes)])\n",
    "\n",
    "    # adding edges with probability p1 for the given nodes\n",
    "    for i in range(n_nodes-n_densest, n_nodes):\n",
    "        for j in range(i+1, n_nodes):\n",
    "            r = np.random.random()\n",
    "            if r < p2:\n",
    "                G2.add_edge(i, j)\n",
    "\n",
    "    G = nx.compose(G1, G2)\n",
    "    if n_densest != 0:\n",
    "        #n_connect=round(min(n_densest, n_nodes-n_densest)/2)+1\n",
    "        n_connect = round(n_densest/2)+1\n",
    "        for i in range(n_connect):\n",
    "            node1 = np.random.randint(n_nodes-n_densest)\n",
    "            node2 = np.random.randint(n_nodes-n_densest, n_nodes)\n",
    "            G.add_edge(node1, node2)\n",
    "\n",
    "    # so there is not isolated parts\n",
    "    if list(nx.isolates(G)) != [] or len(list(G.subgraph(c) for c in nx.connected_components(G))) > 1:\n",
    "        G1.clear()\n",
    "        G2.clear()\n",
    "        G.clear()\n",
    "        G = random_graph_generation(n_nodes, p1, n_densest, p2)\n",
    "\n",
    "    return G\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We construct *AdjToUnit()* now. Considering the singular value decomposition of $A$, from the theory part in the paper we have that if $$A=U\\Sigma V^T \\text{ and } s=\\sigma_{max}(A)$$\n",
    "then we can embed $A$ into a doubled-size unitary $U_A$ with the form:\n",
    "$$U_A:=\\begin{pmatrix}\n",
    "A_s & \\sqrt{I_{n \\times n} -A_s(A_s)^{T}} \\\\\n",
    "\\sqrt{I_{n \\times n} -(A_s)^{T}A_s}  & -(A_s)^{T}\\\\\n",
    "\\end{pmatrix}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def to_unitary(A):\n",
    "    ''' Input: graph A either as:\n",
    "                                 an adjacency matrix of size mxm\n",
    "                                 a networkX graph with m nodes\n",
    "        Output: unitary with size 2mx2m\n",
    "    '''\n",
    "\n",
    "    if type(A) == type(nx.Graph()):\n",
    "        A = nx.to_numpy_array(A)\n",
    "    P1, D, V = linalg.svd(A)\n",
    "\n",
    "    c = np.max(D)\n",
    "    # if it is not complex, then np.sqrt will output nan in complex values\n",
    "    An = np.matrix(A/c, dtype=complex)\n",
    "    P = An\n",
    "    m = len(An)\n",
    "    Q = sqrtm(np.identity(m)-np.dot(An, An.conj().T))\n",
    "    R = sqrtm(np.identity(m)-np.dot(An.conj().T, An))\n",
    "    S = -An.conj().T\n",
    "    Ubmat = np.bmat([[P, Q], [R, S]])\n",
    "    Ubmat = Ubmat.real #to eliminate those \n",
    "    return (np.copy(Ubmat), c)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With *to_unitary()* defined, now we can estimate the permanent of a graph by generating samples. Following the schematic on figure below, we create functions *input_state()* and *post_select()* to preparate the initial state and select the appropriate ones at the output. Finally, everything is combined in *perm_estimation()* function. \n",
    "\n",
    "![BS encoding scheme](BSdevice_graphs.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1, 1, 1, 1, 0, 0, 0, 0])"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def input_state(m):\n",
    "    '''input state for selection of our m modes\n",
    "        returns |1,1,1,...,0,0,0> m ones and m zeros'''\n",
    "    return np.append(np.ones(m), np.zeros(m)).astype(int)\n",
    "\n",
    "input_state(4)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Post selection of samples with photons only on first half modes\n",
    "def post_select(samples):\n",
    "    ''''post select on states that have all modes from m to 2*m as vacuum\n",
    "        can't have collision of first half'''\n",
    "    a = []\n",
    "    m = int(len(samples[0])/2)\n",
    "    for state in samples:\n",
    "        state = list(state)\n",
    "        if all(ele == state[m-1] for ele in state[:m]) and state[m-1] == 1:\n",
    "            # do not need to check if there is vacuum in the second half for several reasons!\n",
    "            a.append(state)\n",
    "    return a"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Computing the permanent estimation by generating samples.\n",
    "\n",
    "$$permE(A)=s^n PermE(A_s)=s^n\\sqrt{\\frac{ \\text{Accepted samples} }{\\text{Total number of samples}}}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def perm_estimation(G, nb_samples, Ns_min=0):\n",
    "    '''Function to estimate the permanent given a graph G, a number of samples Ns and the minimum number of samples Ns_min\n",
    "    Recal that Ns is the number of samples to generate but, due to post-selection, we add Ns_min so that we don't do estimations with accepted samples to be zero\n",
    "    '''\n",
    "    if Ns_min == 0:\n",
    "        Ns_min = Ns\n",
    "\n",
    "    #Sampling_Backend = pcvl.BackendFactory().get_backend(\"SLOS\")\n",
    "    m = G.number_of_nodes()\n",
    "    in_state = input_state(m)\n",
    "\n",
    "    U, c = to_unitary(G)\n",
    "    U_matrix_pcvl = pcvl.Matrix(U)\n",
    "    unitary_component = comp.Unitary(U_matrix_pcvl)\n",
    "    proc = pcvl.Processor(\"SLOS\", unitary_component)\n",
    "    proc.with_input(pcvl.BasicState( in_state))\n",
    "    #simulator = Sampling_Backend(U)\n",
    "\n",
    "    samples = []\n",
    "    i = 0\n",
    "    sampler = Sampler(proc)\n",
    "    while len(samples) < Ns_min:\n",
    "        samples = sampler.samples(nb_samples)['results']\n",
    "        samples = post_select(samples)\n",
    "        i = i+1\n",
    "    print(\"Total number of samples: \", Ns*i)\n",
    "    print(\"Number of samples post:\", len(samples))\n",
    "    perm = (c**m)*np.sqrt(len(samples)/(Ns*i))\n",
    "    return perm\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def perm_estimation(G, nb_samples, Ns_min=0):\n",
    "    if Ns_min == 0:\n",
    "        Ns_min = nb_samples\n",
    "\n",
    "    if type(G) == type(nx.Graph()):\n",
    "        m = G.number_of_nodes()\n",
    "    else:\n",
    "        m = len(G)\n",
    "    in_state = input_state(m)\n",
    "\n",
    "    U, c = to_unitary(G)\n",
    "    U_matrix_pcvl = pcvl.Matrix(U)\n",
    "    unitary_component = comp.Unitary(U_matrix_pcvl)\n",
    "    proc = pcvl.Processor(\"CliffordClifford2017\", unitary_component)\n",
    "    proc.with_input(pcvl.BasicState(in_state))\n",
    "\n",
    "    samples_accepted = []\n",
    "    i = 0\n",
    "    sampler = Sampler(proc)\n",
    "    while len(samples_accepted) < Ns_min:\n",
    "        samples_accepted.append(list(sampler.samples(nb_samples)['results']))\n",
    "        samples_accepted = post_select(samples_accepted)\n",
    "        i = i+1\n",
    "    print(\"Total number of samples: \", nb_samples*i)\n",
    "    print(\"Number of samples post:\", len(samples_accepted))\n",
    "    perm = (c**m)*np.sqrt(len(samples_accepted)/(nb_samples*i))\n",
    "    return perm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def perm_estimation(G, nb_samples, Ns_min=0):\n",
    "    if Ns_min == 0:\n",
    "        Ns_min = nb_samples\n",
    "\n",
    "    if type(G) == type(nx.Graph()):\n",
    "        m = G.number_of_nodes()\n",
    "    else:\n",
    "        m = len(G)\n",
    "    in_state = input_state(m)\n",
    "\n",
    "    U, c = to_unitary(G)\n",
    "    U_matrix_pcvl = pcvl.Matrix(U)\n",
    "    unitary_component = comp.Unitary(U_matrix_pcvl)\n",
    "    proc = pcvl.Processor(\"CliffordClifford2017\", unitary_component)\n",
    "    proc.with_input(pcvl.BasicState(in_state))\n",
    "\n",
    "    samples_accepted = []\n",
    "    i = 0\n",
    "    sampler = Sampler(proc)\n",
    "    while len(samples_accepted) < Ns_min:\n",
    "        samples = list(sampler.samples(nb_samples)['results'])\n",
    "        samples_accepted += post_select(samples)\n",
    "        i = i+1\n",
    "    print(\"Total number of samples: \", nb_samples*i)\n",
    "    print(\"Number of samples post:\", len(samples_accepted))\n",
    "    perm = (c**m)*np.sqrt(len(samples_accepted)/(nb_samples*i))\n",
    "    return perm"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plotting estimation and exact value for permanent over the edge connectivity (probability parameter in *random_graph_generation()*). The plots can vary a lot because that are a lot of variables that can change. Increasing the number of graphs to consider for the mean, increasing the number of samples for each estimation, all of these would help have a plot but since the generation is probabilistic, it should always vary a little."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.7  0.78 0.86 0.94 1.  ]\n",
      "[(150+0j)]\n",
      "[(1146.8666666666666+0j)]\n",
      "[(150+0j), (289.3333333333333+0j)]\n",
      "[(1146.8666666666666+0j), (3439.5333333333333+0j)]\n",
      "[(150+0j), (289.3333333333333+0j), (695.3333333333334+0j)]\n",
      "[(1146.8666666666666+0j), (3439.5333333333333+0j), (4152.133333333333+0j)]\n",
      "[(150+0j), (289.3333333333333+0j), (695.3333333333334+0j), (1211.4666666666667+0j)]\n",
      "[(1146.8666666666666+0j), (3439.5333333333333+0j), (4152.133333333333+0j), (9513.333333333334+0j)]\n",
      "[(150+0j), (289.3333333333333+0j), (695.3333333333334+0j), (1211.4666666666667+0j), (1854+0j)]\n",
      "[(1146.8666666666666+0j), (3439.5333333333333+0j), (4152.133333333333+0j), (9513.333333333334+0j), (14833+0j)]\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/filipagoncalvescarvalho/opt/anaconda3/lib/python3.9/site-packages/matplotlib/cbook/__init__.py:1298: ComplexWarning: Casting complex values to real discards the imaginary part\n",
      "  return np.asarray(x, float)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAEGCAYAAACUzrmNAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAqu0lEQVR4nO3dd3xV9f3H8deHPQMCYYNhy5IVEbVarHVrHXXgwlVR66jVVqu2Um2pW6vS4vppRQVFceGoonWL2oS99wgghJWEEcj4/P64B3vFhFxI7j03yfv5eOSRm+89597PlwO8c873e7/H3B0REZG9qRF2ASIikvwUFiIiUiaFhYiIlElhISIiZVJYiIhImWqFXUC8tGjRwtPS0sIuQ0SkUsnMzNzg7ql7tlfZsEhLSyMjIyPsMkREKhUzW1FSuy5DiYhImRQWIiJSJoWFiIiUSWEhIiJlUliIiEiZFBYiIlImhYWIiJRJYSEiUkUszd7K/e/Pp7CouMJfW2EhIlIF5OYX8KuxGYz/dhUbtu6q8NdXWIiIVHJFxc5vxk9j5cbtjLlgIK2b1Kvw96iyy32IiFQX970/n48XZPPX0/twaOfmcXkPnVmIiFRib0xbzROfLuWCQzty4ZAD4/Y+CgsRkUpqZtYWbpk4k8GdmjHy1N5xfS+FhYhIJbQ+N58RYzNp0aguYy4YSJ1a8f3vXGMWIiKVzM7CIq56IZOcHQVMvPpwmjeqG/f3VFiIiFQi7s4fX5/N1JVbGHPBQHq1TUnI++oylIhIJfLsl8t5JTOL64/pxol92yTsfeMWFmb2jJmtN7PZUW0vm9n04Gu5mU0P2tPMbEfUc49H7TPIzGaZ2WIze9TMLF41i4gksy8WbWDUu/M4vncrbjimW0LfO56Xof4FjAbG7m5w93N3PzazB4GcqO2XuHv/El5nDDAC+Bp4FzgBeK/iyxURSV7LN2zjmnFT6ZraiIfO6U+NGon9vTluZxbu/hmwqaTngrODc4Dxe3sNM2sDpLj7FHd3IsFzegWXKiKS1PLyC7hibAZm8NTwdBrWTfxwc1hjFkcC69x9UVRbJzObZmafmtmRQVs7ICtqm6ygrURmNsLMMswsIzs7u+KrFhFJsOJi57cvT2fphm388/yBdGzeIJQ6wgqL8/jhWcVaoKO7DwBuBMaZWQpQ0nmWl/ai7v6ku6e7e3pqamqFFiwiEoaHJi/kw3nrueOUXhzetUVodST8XMbMagFnAoN2t7n7TmBn8DjTzJYA3YmcSbSP2r09sCZx1YqIhOftmWsY/fFihh3SgeGHxW8pj1iEcWbxc2C+u39/ecnMUs2sZvC4M9ANWOrua4E8MxsSjHMMB94MoWYRkYSavTqH370yg/QDD+Cu0/oQ9kTQeE6dHQ9MAXqYWZaZXR48NYwfD2wfBcw0sxnAq8BV7r57cPxq4GlgMbAEzYQSkSouO28nI8Zm0KxBHcZcOCjuS3nEIm6Xodz9vFLaLymhbSIwsZTtM4A+FVqciEiS2lVYzNUvZLJp+y5evepwUhvHfymPWGi5DxGRJOHujHxrNhkrNvPYeQPo065J2CV9L/xzGxERAeD5r1cw/ttV/HpoF07t1zbscn5AYSEikgS+WrKBOyfN5ZiDWvK743qEXc6PKCxEREK2atN2rnlxKp1aNOTvwxK/lEcsFBYiIiHatrOQK8ZmUOzw9PB0GterHXZJJVJYiIiEpLjYuXHCdBauy2P0+QNIa9Ew7JJKpbAQEQnJIx8t4v0567j95F4c2S25lyhSWIiIhOC9WWt55KNFnDWoPZcdkRZ2OWVSWIiIJNi8tbncOGEGAzo2ZdQZ4S/lEQuFhYhIAm3atosrxmbQpH5tnrhwEHVr1Qy7pJjoE9wiIglSUFTMr1/MZH3eTl658jBaptQLu6SY6cxCRCRB7po0l6+XbuK+Xx5Mvw5Nwy5nnygsREQS4MVvVvD81yu48qjOnD6g1Bt+Ji2FhYhInH2zdCMj35zD0B6p3HzCQWGXs18UFiIicZS1eTtXvziVjs0b8MiwAdRMwqU8YqGwEBGJk+27CrlibCYFRcU8NTydJvWTcymPWGg2lIhIHLg7v39lJvO/y+WZSw6hS2qjsEsqF51ZiIjEwej/LOadWWu59cSDOLpHy7DLKTeFhYhIBftgznc8OHkhZwxoxxVHdg67nAqhsBARqUALvsvjty9Pp1/7Jtx9Zt9KsZRHLOIWFmb2jJmtN7PZUW1/NrPVZjY9+Dop6rlbzWyxmS0ws+Oj2geZ2azguUetqvzJi0iVszlYyqNh3Vo8cVE69WpXjqU8YhHPM4t/ASeU0P6wu/cPvt4FMLNewDCgd7DPP81s95/yGGAE0C34Kuk1RURCVVhUzLXjp/JdTj6PXzSI1k0qz1IesYhbWLj7Z8CmGDc/DXjJ3Xe6+zJgMTDYzNoAKe4+xd0dGAucHpeCRUTK4a/vzOPLxRv525l9GdjxgLDLqXBhjFlca2Yzg8tUu/9E2wGrorbJCtraBY/3bC+RmY0wswwzy8jOzq7oukVESjThv6v411fLufwnnThrUPuwy4mLRIfFGKAL0B9YCzwYtJc0DuF7aS+Ruz/p7ununp6amtx3nRKRqiFzxSZuf2MWR3Zrwa0nVs6lPGKR0LBw93XuXuTuxcBTwODgqSygQ9Sm7YE1QXv7EtpFREK3ZssOrnx+Ku2a1mf0eQOpVbPqTjBNaM+CMYjdzgB2z5R6CxhmZnXNrBORgexv3X0tkGdmQ4JZUMOBNxNZs4hISfILirjy+UzyC4p4+uJ0mjSovEt5xCJuy32Y2XhgKNDCzLKAkcBQM+tP5FLScuBKAHefY2YTgLlAIXCNuxcFL3U1kZlV9YH3gi8RkdC4Oze/OpPZa3J4eng6XVs2DrukuItbWLj7eSU0/99eth8FjCqhPQPoU4GliYiUy5hPl/DWjDX8/vgeHNOzVdjlJETVvcAmIhIHH81bx/3vL+DUfm359dAuYZeTMAoLEZEYLV6fx29emk7vtinc98uDq8xSHrFQWIiIxCBnewG/ei6DerVr8uRF6dSvU3WW8oiFwkJEpAy7l/JYvWUHj184kLZN64ddUsLp5kciImW45735fL5oA/f+si/pac3CLicUOrMQEdmLiZlZPP3FMi45PI1zD+kYdjmhUViIiJRi2srN3Pr6LA7v0pzbT+4ZdjmhUliIiJRgXW4+Vz6fSeuUevzj/IHUrsJLecSievdeRKQE+QVFjHg+k207C3lqeDoHNKwTdkmh0wC3iEgUd+e212YxY9UWnrxoED1aV/2lPGKhMwsRkShPf76M16at5sZju3Nc79Zhl5M0FBYiIoFPF2Zz93vzOKlva677Wdewy0kqCgsREWBp9lauHTeVHq1TeODsftVqKY9YKCxEpNrLzS/gV2MzqF2zBk8NH0SDOhrO3ZPCQkSqtaJi5/rx01i5cTtjLhhI+wMahF1SUlJ8iki1dt/78/lkQTajzujDoZ2bh11O0tKZhYhUW29MW80Tny7lwiEdueDQA8MuJ6kpLESkWpqZtYVbJs7k0E7NGHlq77DLSXoKCxGpdtbn5jNibCYtGtXlnxdoKY9YaMxCRKqVnYVFXPVCJjk7Cph49eE0b1Q37JIqhbjFqZk9Y2brzWx2VNv9ZjbfzGaa2etm1jRoTzOzHWY2Pfh6PGqfQWY2y8wWm9mjpsnPIrKf3J0/vj6bqSu38NA5/ejVNiXskiqNeJ57/Qs4YY+2yUAfdz8YWAjcGvXcEnfvH3xdFdU+BhgBdAu+9nxNEZGYPPvlcl7JzOL6Y7pxYt82YZdTqcQtLNz9M2DTHm0fuHth8OPXQPu9vYaZtQFS3H2KuzswFjg9DuWKSBX3xaINjHp3Hsf3bsUNx3QLu5xKJ8xRncuA96J+7mRm08zsUzM7MmhrB2RFbZMVtJXIzEaYWYaZZWRnZ1d8xSJSKS3fsI1rxk2la2ojHjqnPzVq6Gr2vgolLMzsdqAQeDFoWgt0dPcBwI3AODNLAUo6ol7a67r7k+6e7u7pqampFV22iFRCefkFXDE2AzN4ang6DetqXs/+SPifmpldDJwCHBNcWsLddwI7g8eZZrYE6E7kTCL6UlV7YE1iKxaRyqq42Pnty9NZumEbz182mI7NtZTH/kromYWZnQDcAvzC3bdHtaeaWc3gcWciA9lL3X0tkGdmQ4JZUMOBNxNZs4hUXg9NXsiH89Zzxym9OLxri7DLqdRiCgsz+9FE5JLa9nh+PDAF6GFmWWZ2OTAaaAxM3mOK7FHATDObAbwKXOXuuwfHrwaeBhYDS/jhOIeIyI8UFTuPfbSI0R8vZtghHRh+mJbyKC8LrgTtfSOzqe4+sKy2ZJKenu4ZGRlhlyEiCbY+N58bXp7OV0s2clr/ttx/Vj/q1NIntGNlZpnunr5n+17HLMysNZHZR/XNbAD/G3BOAXTxT0SSyqcLs7nx5els31XEfWcdzNmD2usmRhWkrAHu44FLiAwsPxTVngfcFqeaRET2SUFRMQ98sIAnPl1Kj1aN+ccFA+jasnHYZVUpew0Ld38OeM7MfunuExNUk4hIzFZt2s71L01j2sotnH9oR+44pRf1atcMu6wqJ9aps2+b2flAWvQ+7n5XPIoSEYnFv2ev5eZXZ+IO/zh/ICcfrCU84iXWsHgTyAEyCT4PISISlvyCIka9M4/nv15Bv/ZNeOy8gfoMRZzFGhbt3V0L+IlI6JZkb+XacdOYtzaXEUd15nfH9dBspwSINSy+MrO+7j4rrtWIiOzFq5lZ3PHmbOrVrsmzlxzC0Qe1DLukaiPWsPgJcImZLSNyGcoAD5YaFxGJq207C/nTG7N5bdpqhnRuxt/PHUDrJvXCLqtaiTUsToxrFSIipZizJofrxk1j+cZt3PDzblz3s27U1KqxCRdTWLj7CjP7CdDN3Z81s1SgUXxLE5HqzN0ZO2UFo96ZxwENazPuiiEM6dw87LKqrZjCwsxGAulAD+BZoDbwAnBE/EoTkeoqZ3sBN0+cwftz1nF0j1QeOLuf7pUdslgvQ50BDACmArj7GjPTxyNFpMJlrtjE9eOnsz4vnz+e3JPLjuikmxUlgVjDYpe7u5k5gJk1jGNNIlINFRc7j3+2hAc/WEjbpvV49arD6dehadhlSSDWsJhgZk8ATc3sCiK3RH0qfmWJSHWSnbeTGydM5/NFGzj54DbcfWZfUurVDrssiRLrAPcDZnYskEtk3OIOd58c18pEpFr4YtEGbnh5Onn5BfztjL6cN7iDVopNQjHfVjUIBwWEiFSIwqJiHv5wIf/8ZAldUhvxwq8Gc1DrlLDLklLEOhvqTOBeoCWRD+Tt/lCejqyI7LPVW3bwm/HTyFixmXPTOzDyF71oUCfm310lBLEenfuAU919XjyLEZGq74M53/H7V2dSWFTMI8P6c1r/dmGXJDGINSzWKShEpDx2FhZx97vz+ddXy+nTLoXR5w0krYUmVlYWsYZFhpm9DLxB1BLl7v5aPIoSkapl2YZtXDtuKnPW5HLpEWn84cSDqFtLNyiqTGJd1zcF2A4cB5wafJ2ytx3M7BkzW29ms6PampnZZDNbFHw/IOq5W81ssZktMLPjo9oHmdms4LlHTdMkRCqVN6at5pRHP2f1lh08NTydkaf2VlBUQrFOnb10P177X8BoYGxU2x+Aj9z9HjP7Q/DzLWbWCxgG9AbaAh+aWXd3LwLGACOAr4F3gROA9/ajHhFJoO27Chn55hxeyczikLQDeGTYANo2rR92WbKfYp0NVQ+4nMh/5t+vC+zul5W2j7t/ZmZpezSfBgwNHj8HfALcErS/5O47gWVmthgYbGbLgRR3nxLUMRY4HYWFSFKb/10u146bxpLsrVz3s6785phu1KqpGxRVZrEeveeB1sDxwKdAeyBvP96vlbuvBQi+775zSTtgVdR2WUFbu+Dxnu0lMrMRZpZhZhnZ2dn7UZ6IlIe78+I3Kzht9Jfk7CjghcsP5abjeigoqoBYj2BXd/8TsM3dnwNOBvpWYB0ljUP4XtpL5O5Punu6u6enpqZWWHEiUracHQVcO24at78+m8GdmvHu9UdyRNcWYZclFSTW2VAFwfctZtYH+A5I24/3W2dmbdx9rZm1AdYH7VlAh6jt2gNrgvb2JbSLSBKZvmoL146bytqcfG454SCuPKqzVoqtYmI9s3gymLn0J+AtYC6RD+rtq7eAi4PHFwNvRrUPM7O6ZtYJ6AZ8G1yqyjOzIcEsqOFR+4hIyIqLnSc/W8JZY77CHSZceRhXD+2ioKiCYp0N9XTw8FOgcyz7mNl4IoPZLcwsCxgJ3ENkBdvLgZXA2cHrzzGzCURCqBC4JpgJBXA1kZlV9YkMbGtwWyQJbNy6k5temcEnC7I5oXdr7v3lwTRpoJViqypzL3UI4H8bmTUl8lt9GlEB4+7Xx6uw8kpPT/eMjIywyxCpkqYs2cgNL09j8/YC/nRyTy4ccqBWiq0izCzT3dP3bI91zOJdIp9zmAUUV2RhIlJ5FBU7j3y0iMf+s4hOzRvyzCWH0Lttk7DLkgSINSzqufuNca1ERJLadzn5XP/SNL5dtokzB7bjL6f1oWFdrRRbXcR6pJ8P7pD3Nj9cG2pTXKoSkaTyn/nruGnCDHYWFvPg2f345aD2Ze8kVUrM9+AG7gdu53+fc3BiHOwWkcppV2Ex9/17Pk9/sYyebVIYff4AuqQ2CrssCUGsYXEjkQ/mbYhnMSKSPFZu3M5146cyIyuH4YcdyG0n9aRebS0AWF3FGhZziKw6KyLVwKQZa7jttVmYweMXDuSEPm3CLklCFmtYFAHTzexjfjhmkbRTZ0Vk3+3YVcRdb89h/LerGNixKY8MG0CHZg3CLkuSQKxh8UbwJSJV1KJ1eVwzbioL123l6qFduPHY7tTWAoASKDMszKwmcJG7/zwB9YhIgrk7EzJWMfKtOTSqW4uxlw3mqO5aiFN+qMywcPciM9tuZk3cPScRRYlIYuTlF3D767N5a8YajujanIfP7U/LxvXK3lGqnVgvQ+UDs8xsMrBtd6PGLEQqr1lZOVw7fipZm3fw++N7cNVPu1BTCwBKKWINi3eCLxGp5NydZ75czj3vzSO1UV1eGjGEQ9KahV2WJLlYV519zszqAx3dfUGcaxKRONm8bRe/f3UGH85bz897tuKBsw+maYM6YZcllUCs9+A+FXgAqAN0MrP+wF3u/os41iYiFejbZZu4fvw0Nm3bxchTe3HJ4WlaKVZiFutlqD8Dg4FPANx9enCTIhFJckXFzj8+XszfP1xIx2YNeO3Xh9OnnVaKlX0Ta1gUunvOHr+FlH0jDBEJ1frcfG54eTpfLdnIaf3bMuqMvjTSSrGyH2L9WzPbzM4HappZN+B64Kv4lSUi5fXJgvXcNGEG23cVcd9ZB3P2oPa67CT7LdaPZ14H9Cay1Mc4IAe4IU41iUg5FBQVc/d787jk2f/SolFdJl13BOekd1BQSLns9czCzOoBVwFdidwl7zB3L0xEYSKy71Zt2s5146cxfdUWzj+0I3ec0ksrxUqFKOsy1HNAAfA5cCLQE51RiCSl92at5eaJM8HhH+cP5OSDtVKsVJyywqKXu/cFMLP/A74t7xuaWQ/g5aimzsAdQFPgCiA7aL/N3d8N9rkVuJzI6rfXu/v75a1DpKrILyjir+/M5YWvV9KvfRMeO28gHZtrpVipWGWFRcHuB+5eWBHXPIMP9fWH7xcpXA28DlwKPOzuD0Rvb2a9gGFExkzaAh+aWXd3Lyp3MSKV3OL1W7l23FTmf5fHiKM687vjelCnllaKlYpXVlj0M7Pc4LEB9YOfDXB3Tynn+x8DLHH3FXsJotOAl9x9J7DMzBYT+czHlHK+t0il9mpmFn96Yzb169Tk2UsO4eiDWoZdklRhew0Ld4/3yNgwYHzUz9ea2XAgA7jJ3TcD7YCvo7bJCtp+xMxGACMAOnbsGJeCRcK2dWchd7wxm9emrWZI52b8/dwBtG6ilWIlvkI7XzWzOsAvgFeCpjFAFyKXqNYCD+7etITdS/xAoLs/6e7p7p6emqr1+KXqmbMmh1Mf+4I3pq/mhp9348VfDVFQSEKE+VHOE4Gp7r4OYPd3ADN7Cng7+DEL6BC1X3tgTaKKFEkG7s7YKSsY9c48DmhYm3FXDGFI5+ZhlyXVSJhhcR5Rl6DMrI27rw1+PAOYHTx+CxhnZg8RGeDuRgXMyhKpLHK2F3DzxBm8P2cdR/dI5YGz+9G8Ud2wy5JqJpSwMLMGwLHAlVHN9wWr2TqwfPdz7j7HzCYAc4FC4BrNhJLqInPFJq4fP531efn88eSeXHZEJ2roBkUSglDCwt23A833aLtoL9uPAkbFuy6RZFFc7Dz+2RIe/GAhbZvW49WrDqdfh6ZhlyXVmJafFEky01dt4c5Jc5i2cgsnH9yGu8/sS0q92mGXJdWcwkIkSazPzefefy9g4tQsWjSqy4Nn9+PMge20AKAkBYWFSMh2FhbxzBfLGf2fRewqKubKn3bm2qO70lhnE5JEFBYiIXF3Ppq3nr++M5flG7fz856tuP3knnRq0TDs0kR+RGEhEoLF6/O4c9JcPl+0ga4tG/HcZYP5aXd9kFSSl8JCJIFydhTw9w8XMnbKChrUqckdp/TiosMOpHZNLf4nyU1hIZIARcXOy/9dxQMfLGDz9l2cN7gjNx3bXR+uk0pDYSESZ98s3cidk+Yyd20ug9OaMfIXvejdtknYZYnsE4WFSJxkbd7O3e/N552Za2nbpB6jzx/AyX3baCqsVEoKC5EKtmNXEY9/uoTHP12CGdzw825ceVQX6tfRvbCl8lJYiFQQd+ftmWu5+915rMnJ55SD23DrST1p17R+2KWJlJvCQqQCzF6dw12T5vLt8k30apPCw+f251AtIS5ViMJCpBw2bt3JAx8s5KX/ruSABnX42xl9OfeQDtTUyrBSxSgsRPZDQVExY6es4O8fLmTHriIuPbwTvzmmG00aaIkOqZoUFiL76NOF2dw1aQ5LsrdxZLcWjDy1F11bNg67LJG4UliIxGjZhm2MemcuH85bT1rzBjw9PJ1jerbUVFipFhQWImXIyy9g9MeLeeaLZdSpWYM/nHgQlx6RRt1amgor1YfCQqQUxcXOxKlZ3Pf+ArLzdnLWoPbcfEIPWjauF3ZpIgmnsBApwdSVm7nzrTnMyMphQMemPD08Xbc1lWpNYSESZV1uPve+N5/Xpq2mZeO6PHROP07v344amgor1VwoYWFmy4E8oAgodPd0M2sGvAykAcuBc9x9c7D9rcDlwfbXu/v7IZQtVVh+QRH/98Uy/vHxYgqLnF8P7cI1R3elYV39PiUC4Z5ZHO3uG6J+/gPwkbvfY2Z/CH6+xcx6AcOA3kBb4EMz6+7uRYkvWaoad+eDuesY9c48Vm7aznG9WvHHk3vRsXmDsEsTSSrJ9GvTacDQ4PFzwCfALUH7S+6+E1hmZouBwcCUEGqUKmThujzunDSHLxdvpHurRrxw+aH8pFuLsMsSSUphhYUDH5iZA0+4+5NAK3dfC+Dua82sZbBtO+DrqH2zgrYfMbMRwAiAjh07xqt2qeS2bN/Fw5MX8sI3K2lUtxZ3/qI3FxzakVq6W51IqcIKiyPcfU0QCJPNbP5eti1pZNFL2jAInScB0tPTS9xGqq/ComLG/3cVD32wgJwdBZx/aEduPLYHzRrWCbs0kaQXSli4+5rg+3oze53IZaV1ZtYmOKtoA6wPNs8COkTt3h5Yk9CCpdKbsmQjd06aw/zv8hjSuRkjT+1NzzYpYZclUmkk/LzbzBqaWePdj4HjgNnAW8DFwWYXA28Gj98ChplZXTPrBHQDvk1s1VJZrdq0nV+/mMl5T31NXn4hYy4YyPgrhigoRPZRGGcWrYDXg/V0agHj3P3fZvZfYIKZXQ6sBM4GcPc5ZjYBmAsUAtdoJpSUZfuuQsZ8soQnPltKTTNuOrY7VxzVmXq1tUSHyP5IeFi4+1KgXwntG4FjStlnFDAqzqVJFeDuvDVjDXe/O5/vcvM5rX9b/nDiQbRporvViZRHMk2dFSmXWVk53DlpDhkrNtOnXQqjzx9AelqzsMsSqRIUFlLpbdi6k/v/vYAJmato3rAO9/6yL2cN0t3qRCqSwkIqrV2FxTz31XIe/WgROwqK+NVPOnHdMd1Iqae71YlUNIWFVEofL1jPX96ey9LsbQztkcqfTulFl9RGYZclUmUpLKRSWZq9lb+8PZePF2TTuUVDnr3kEI4+qGXZO4pIuSgspFLIzS/gsY8W8eyXy6lfuya3n9STiw9Po04tLdEhkggKC0lqxcXOK5mruP/9BWzctotzBnXgd8f3ILVx3bBLE6lWFBaStDKWb+LOSXOZtTqHQQcewLOXDKZv+yZhlyVSLSksJOmszdnBPe/N583pa2idUo9HhvXnF/3aEnzqX0RCoLCQpJFfUMRTny3ln58socid637WlauHdqFBHf01FQmb/hVK6Nydf8/+jlHvziNr8w5O7NOa207qSYdmuludSLJQWEio5q3N5a5Jc5mydCMHtW7MuCsO5fAuuludSLJRWEgoNm/bxUOTF/LiNytIqV+bv5zeh/MO6aC71YkkKYWFJFRhUTEvfrOShyYvZOvOQi4aciC/PbY7TRvobnUiyUxhIQnz5eIN3DlpDgvXbeWIrs2545Te9GjdOOyyRCQGCguJu5UbtzPq3bm8P2cdHZrV54mLBnFcr1aaCitSiSgsJG627Szkn58s5qnPl1GrhvH743tw+U866W51IpWQwkIqXHGx8+aM1dzz3nzW5e7kzAHtuPmEg2jdpF7YpYnIflJYSLmty81nVlYOs1bnMHt1DjNX55Cdt5N+7ZvwzwsGMejAA8IuUUTKSWEhMXN31uXuZNbq/wXDrCAYAGoYdEltxJFdW/DTHqmcenBbauhudSJVQsLDwsw6AGOB1kAx8KS7P2JmfwauALKDTW9z93eDfW4FLgeKgOvd/f1E113duDvfBWcMu0Nh1upcNmz9XzB0bdmII7u1oG+7JvRt14RebVO0NIdIFRXGv+xC4CZ3n2pmjYFMM5scPPewuz8QvbGZ9QKGAb2BtsCHZtbd3YsSWnUV5u6szcn/wdnC7NU5bNi6C4gEQ7eWjflp91T6tkuhb/sm9GyjYBCpThL+r93d1wJrg8d5ZjYPaLeXXU4DXnL3ncAyM1sMDAamxL3YKsjdWb1lxw/OFuaszmHjtkgw1KxhdGvZiKE9WtK3XRP6tGtCrzYp1K+jGUwi1VmovxqaWRowAPgGOAK41syGAxlEzj42EwmSr6N2y6KUcDGzEcAIgI4dO8av8ErC3cnaHB0MOcxZk8umqGDo3qoxx/SMBEPvIBg0tVVE9hRaWJhZI2AicIO755rZGOAvgAffHwQuA0oaIfWSXtPdnwSeBEhPTy9xm6pqdzBEDz7PXp3D5u0FANQKguHYnq3o0z4yxnBQ68YKBhGJSShhYWa1iQTFi+7+GoC7r4t6/ing7eDHLKBD1O7tgTUJKjUpuTurNu0RDGty2BIVDD1aN+b43q3pEww+91AwiEg5hDEbyoD/A+a5+0NR7W2C8QyAM4DZweO3gHFm9hCRAe5uwLcJLDlU7s7KTdv3OGPIJWdHJBhq14wEw4l9fhgMdWspGESk4oRxZnEEcBEwy8ymB223AeeZWX8il5iWA1cCuPscM5sAzCUyk+qaqjoTyt1ZsXH7j2Yl5eYXAlCnZg16tG7MSX3bfD9dtXvrRgoGEYk7c6+al/bT09M9IyMj7DJKVVzsrNgUFQxZkUtJeVHBcFCbxt+fLfRt14TurRpTp5bu9yAi8WNmme6evme7JsonQHGxs3zjth+cMcxZnUveziAYatWgZ+vG/KJf2++nqyoYRCSZKCwqWHGxs2zjtu/PFnZPV90aHQxtUjhtwA+DobbuECciSUxhUQ5Fxc6yDVsjg89Zucxek8PcqGCoGwTDGQPafR8M3Vo1UjCISKWjsIhRUbGzNHvrD2YlzV2Ty7ZdkbH2erVr0KtNCr8c2I7ewRhDt5aNdE9pEakSFBYlKCp2lmRv/cGy23PX5rI9CIb6tWvSq20KZ6d3+H4AuktqQwWDiFRZCos9XP6v//LVko3sKPhfMPRum8I56R0is5LaN6FLaiNqaultEalGFBZ7SGvRkI7NG3w/XbWzgkFERGGxpz+d0ivsEkREko4usouISJkUFiIiUiaFhYiIlElhISIiZVJYiIhImRQWIiJSJoWFiIiUSWEhIiJlqrI3PzKzbGDFfu7eAthQgeWEqar0par0A9SXZFVV+lLefhzo7ql7NlbZsCgPM8so6U5RlVFV6UtV6QeoL8mqqvQlXv3QZSgRESmTwkJERMqksCjZk2EXUIGqSl+qSj9AfUlWVaUvcemHxixERKRMOrMQEZEyKSxERKRM1SoszOwEM1tgZovN7A8lPP97M5sefM02syIzaxbLvolWzr4sN7NZwXMZia/+R7WW1ZcmZjbJzGaY2RwzuzTWfROtnH1JmuMSQz8OMLPXzWymmX1rZn1i3TfRytmXZDomz5jZejObXcrzZmaPBv2caWYDo54r/zFx92rxBdQElgCdgTrADKDXXrY/FfjP/uybzH0Jfl4OtAj7mMTaF+A24N7gcSqwKdi20h2X0vqSTMclxn7cD4wMHh8EfLQ/fzeTuS/JdEyCWo4CBgKzS3n+JOA9wIAhwDcVeUyq05nFYGCxuy91913AS8Bpe9n+PGD8fu4bb+XpS7KJpS8ONDYzAxoR+Q+2MMZ9E6k8fUkmsfSjF/ARgLvPB9LMrFWM+yZSefqSVNz9MyJ/X0pzGjDWI74GmppZGyromFSnsGgHrIr6OSto+xEzawCcAEzc130TpDx9gch/WB+YWaaZjYhblbGJpS+jgZ7AGmAW8Bt3L45x30QqT18geY5LLP2YAZwJYGaDgQOB9jHum0jl6QskzzGJRWl9rZBjUqtcpVUuVkJbafOGTwW+dPfdKb4v+yZCefoCcIS7rzGzlsBkM5sf/NYShlj6cjwwHfgZ0IVIzZ/HuG8i7Xdf3D2X5DkusfTjHuARM5tOJPSmETlDqozHpLS+QPIck1iU1tcKOSbV6cwiC+gQ9XN7Ir/dlWQYP7xssy/7JkJ5+oK7rwm+rwdeJ3KaGpZY+nIp8Fpwer0YWEbk2nJlPC6l9SWZjkuZ/XD3XHe/1N37A8OJjL8si2XfBCtPX5LpmMSitL5WzDEJe9AmUV9EzqKWAp343yBP7xK2a0LkumDDfd23kvSlIdA46vFXwAnJ3BdgDPDn4HErYDWRlTUr3XHZS1+S5rjE2I+m/G9g/goi18or5b+VvfQlaY5JVK1plD7AfTI/HOD+tiKPSWidDukP+iRgIZGZAbcHbVcBV0VtcwnwUiz7Vsa+EJkRMSP4mlMZ+gK0BT4gcolgNnBhZT0upfUl2Y5LDP04DFgEzAdeAw6oxMekxL4k4TEZD6wFCoicLVy+Rz8M+EfQz1lAekUeEy33ISIiZapOYxYiIrKfFBYiIlImhYWIiJRJYSEiImVSWIiISJkUFlLlBCvsTo/6Kmml0aFm9nYY9e3JzD4xs/R92P4SMxtdynNfBd/Tdq9OambpZvZo8HiomR1eEXVL9VKdlvuQ6mOHRz6NmzTMrJa7x33BQHf/URC4ewawe3ntocBWIh8wE4mZziyk2gjW9J9vZl8QLBwXtKea2WQzm2pmT5jZCjNrETx3YXCPg+nBczVLeN3lZnZvsN23ZtY1aP+XmT1kZh8D95pZfzP7OrjXwOtmdkDUy1xoZl9Z5N4jg4P9Bwdt04LvPaK272Bm/w7uUTAyqpatJdQ31MzeNrM0Ih/i+m3QnyPNbJmZ1Q62Swn6Urscf8xSRSkspCqqv8dlqHPNrB7wFJGFFY8EWkdtP5LI/T4GEln/pyOAmfUEziWymFx/oAi4oJT3zHX3wURWlf17VHt34OfufhMwFrjF3Q8m8gnbkVHbNQzOCn4NPBO0zQeOcvcBwB3A36K2HxzU0h84O5bLWO6+HHgceNjd+7v758AnRJaJgMg6YhPdvaCs15LqR5ehpCr60WUoM+sPLHP3RcHPLwC7l5z+CXAGgLv/28w2B+3HAIOA/0ZuP0F9YH0p7zk+6vvDUe2vuHuRmTUBmrr7p0H7c8Are+7v7p8Fv+E3BRoDz5lZNyKrhEb/xj/Z3TcGfXkt6MP+3MntaeBm4A0iixxesR+vIdWAwkKqk9LWtilpCefd7c+5+637+NrRj7fFUhg/rs2BvwAfu/sZwSWkT8rYfp+5+5fBYPhPgZruXuItO0V0GUqqi/lAJzPrEvx8XtRzXwDnAJjZccDusYSPgLOCexlgZs3M7MBSXv/cqO9T9nzS3XOAzWZ2ZNB0EfBp1CbnBu/xEyAn2L4JkVVpIbIoZLRjg3rqA6cDX5ZS157yiJyxRBtL5Mzm2RhfQ6ohhYVURXuOWdzj7vlELju9Ewxwr4ja/k7gODObCpxIZGXPPHefC/yRyJ3SZgKTgTalvGddM/sG+A3w21K2uRi4P3it/sBdUc9tDqa9Pk5kNVGA+4C7zexLIvdRjvYF8DyRGylNDGY8xWIScMbuAe6g7UUiAZmst96VJKBVZ6XaM7O6QJG7F5rZYcCYfZl6a2bLiSwHvSFOJcaVmZ0FnObuF4VdiyQvjVmIRGY/TTCzGsAuqtEgr5k9RuRs6qSwa5HkpjMLEREpk8YsRESkTAoLEREpk8JCRETKpLAQEZEyKSxERKRM/w/4wYOJEP53sgAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "probs = np.append(np.arange(0.7, 1, 0.08), 1)\n",
    "print(probs)\n",
    "exactList1 = []\n",
    "exactList2 = []\n",
    "for i in probs:\n",
    "    exactValues1 = [qc.permanent_cx(nx.to_numpy_array(\n",
    "        random_graph_generation(7, i))) for _ in range(15)]\n",
    "    exactList1.append(np.mean(exactValues1))\n",
    "    print(exactList1)\n",
    "    exactValues2 = [qc.permanent_cx(nx.to_numpy_array(\n",
    "        random_graph_generation(8, i))) for _ in range(15)]\n",
    "    exactList2.append(np.mean(exactValues2))\n",
    "    print(exactList2)\n",
    "\n",
    "plt.plot(probs, exactList1, label=\"Size 7\")\n",
    "\n",
    "# plt.legend()\n",
    "plt.xlabel('Edge probability')\n",
    "plt.ylabel('Permanent')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/filipagoncalvescarvalho/opt/anaconda3/lib/python3.9/site-packages/matplotlib/cbook/__init__.py:1298: ComplexWarning: Casting complex values to real discards the imaginary part\n",
      "  return np.asarray(x, float)\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaIAAAEOCAYAAADL1RP6AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA3x0lEQVR4nO3deXxU9bn48c+TEJKQhDWEsO9CWAQkLijiBiruWxf1uq/X2lZra3+99tdr+7u3t9UuqK0V1Fal0treLiIKKFYNAi5BEBL2HQIhC0tIyJ7n98c5A4dhIJmQyZlJnvfrNa/MnPPMd55vDpOHc873nK+oKsYYY4xf4vxOwBhjTPtmhcgYY4yvrBAZY4zxlRUiY4wxvrJCZIwxxlcd/E4gFqWnp+ugQYP8TsMYY2LG8uXLS1S1Z6h1VoiaYdCgQeTm5vqdhjHGxAwR2X6idXZozhhjjK+sEBljjPGVFSJjjDG+skJkjDHGV1aIjDHG+MoKkTHGGF9ZITLGGOMrK0TGGGMatXvTAVYu2kEkpg6yQmSMMeakKg5Us3BWHnk5BdRW17d4+1aIjDHGnFB9XQMLZq2mprqe6Q+OpWNSy9+QxwqRMcaYE1r8l40Ubinjktuz6NEnNSKfYYXIGGNMSGuW7CY/p4AJlw5g2MSMiH2OFSJjjDHH2butjI/+tJ5+I7txzrVDIvpZVoiMMcYc43BZDQtmrialcyKX3TuGuPjIlgorRMYYY45oqG/g3ZfyqCyvZfqDY0lKTYj4Z1ohMsYYc8TSf2ymYMMBLrp1BD0HpLXKZ1ohMsYYA8CGzwv5ctFOxl7UjxHn9G61z23VQiQiU0RkrogUiIiKyJ0niZ3lxnw3aHmiiDwnIiUiUuG21y8oppuIzBaRg+5jtoh0DYoZICJvuW2UiMizItKxJftrjDGxomRXOR+8to7ew7pw3k3DWvWzW3uPKBXIA74NVJ4oSERuAs4EdodYPQO4EbgZOB/oDMwTkXhPzBzgDGA6cLn7fLan/XjgbSDNbeNm4Cbgl83rljHGxK6qilrmv7CKxE4duOy+McRHeHBCsJa/RPYkVPUd4B0AEXklVIyIDASeAaYC84PWdQHuAe5S1ffcZbcB2934hSKShVN8JqvqUjfmAWCxiIxQ1fXApcBoYKCq7nRjHgdeEpEnVLWsRTtujDFRqqFBee/3+ZTvr+b6x84gpUtiq+cQVeeIRKQD8Cfgv1R1bYiQiUAC8G5ggVtI1gLnuosmAeXAUs/7lgAVQTFrA0XItRBIdD8jVG73i0iuiOQWFxeH2zVjjIlKn8/byo78fZz/tdPIHNLFlxyiqhABPwZKVfV3J1ifCdQDJUHL97rrAjHF6rlFrPu8KChmb1AbJW7bmYSgqrNUNVtVs3v27NnE7hhjTPTasrKY3He2kXVeb0af38e3PFr10NzJiMgFwJ3A+Oa8HfDemzzUfcqbEnOy5cYY02bsL6xg0StryBiYxpSvn4aI+JZLNO0RXQT0BvaISJ2I1AEDgZ+LyC43phCIB9KD3pvB0T2cQiBDPL9V93nPoJjgPZ90t+3gPSVjjGlTairrmP/CajokxHH5A2PpkBDf+JsiKJoK0fPA6Th7RIHHbuDXwCVuzHKgFpgWeJM7dDuLo+eEluGMzpvkaXsSkBIUkxU07HsaUO1+hjHGtEmqyvuvruVAUSWX3TuGtO5JfqfUuofmRCQVCAxQjwMGiMh4YJ+q7sA5j+ONrwUK3ZFuqOpBEXkZeFpEioBS4FfAKmCRG7NWRBYAM0XkPpxDcjOBeYF2cAY75AOvichjQA/gaeBFGzFnjGnLvli4nS0riznvpmH0HdHN73SA1t8jygZWuI9knMEJK4CfhNHGo8DfgTdwRsOVA1erqnfawFuBL3EKzkL3+W2BlW7slcBht4033DaPuXjWGGPakh35pXzy5haGn9mLcZf09zudIyQS84+3ddnZ2Zqbm+t3GsYY02QHiyv56/98Tmq3JG58fCIJia17XkhElqtqdqh10XSOyBhjTATU1tQzf+ZqAKY/OKbVi1BjrBAZY0wbpqp8MHsdpQXlTLtnNF16dvI7peNYITLGmDZs1b92sfHzvZx99RAGju7hdzohWSEyxpg2qmDDfpb8bRODx6Uz8fKBfqdzQlaIjDGmDSrfX8XCF/PompHM1DtHIXH+3TmhMVaIjDGmjamvbWD+zDzqahuY/uBYOiZHzd3cQrJCZIwxbUzOn9dTtK2MqXeMoltmit/pNMoKkTHGtCH5iwtYs2QPEy8fyJAJsTFTgBUiY4xpIwq3HiTnjQ0MGNWds64Z4nc6TWaFyBhj2oDDZTUsmJlHatdEpt0zmrgoHpwQzAqRMcbEuPr6BhbMWk11RS3THxxLUkqC3ymFxQqRMcbEuKX/u4k9mw5y0W0jSe+X5nc6YbNCZIwxMWz9p4Ws+mAX4y7uz2lnBc/3GRusEBljTIwq3nGID/64jj7DuzLpxqF+p9NsVoiMMSYGVZXXMn/mapJTE7jsvjHEx8fun/PYzdwYY9qphgbl3d/nU3GwmsvvH0unzh39TumUWCEyxpgY8+mbW9i5Zh8X3DyCXoM7+53OKbNCZIwxMWTzF0V8sXA7o87vw6jz+vidToto1UIkIlNEZK6IFIiIisidnnUJIvJzEVklIhUiskdE5ojIgKA2EkXkOREpcePmiki/oJhuIjJbRA66j9ki0jUoZoCIvOW2USIiz4pIbO/fGmPatH27K3j/1bX0GtyZKV89ze90Wkxr7xGlAnnAt4HKoHWdgDOA/3Z/Xgv0BxaIiPfWsTOAG4GbgfOBzsA8EfHOfTvHbWM6cLn7fHZgpRv7NpDmtnEzcBPwyxboozHGtLjqyjrmz1xNh8R4Lr9/LPEJbeeAVqveG1xV3wHeARCRV4LWHQSmeZeJyANAPpAFrBaRLsA9wF2q+p4bcxuwHZgKLBSRLJziM1lVl3raWSwiI1R1PXApMBoYqKo73ZjHgZdE5AlVLYtE/40xpjm0QVn0hzWUFVdy7aMTSO2W6HdKLSraS2rgLNx+9+dEIAF4NxDgFpK1wLnuoklAObDU084SoCIoZm2gCLkWAonuZxxHRO4XkVwRyS0uLm52h4wxJly587exbVUJ531lGH2Gd/U7nRYXtYXIPV/zS+AtVd3lLs4E6oGSoPC97rpATLGqamCl+7woKGZvUBslbtshL01W1Vmqmq2q2T17xsat1Y0xsW/b6hI+m7eVEWdnMvbCfo2/IQZF5bR97jmhPwJdgWua8hZAPa+1mTEnW26MMa3qQNFh3vv9GtL7pXLBrSMQiZ07aocj6vaI3CL0J+B04BJVLfWsLgTigfSgt2VwdA+nEMgQzxZzn/cMigne80l32w7eUzLGmFZXW13P/BdWI3Ew/YGxJHSMb/xNMSqqCpGIJABv4BShi1S1MChkOVCLZ1CDO3Q7i6PnhJbhjM6b5HnfJCAlKCYraNj3NKDa/QxjjPGNqvKv2WvZv6eCy+4ZQ+f0ZL9TiqhWPTQnIqnAMPdlHDBARMYD+4DdwF+BM4GrARWRwF7LQVWtVNWDIvIy8LSIFAGlwK+AVcAiAFVdKyILgJkich/OIbmZwDx3xBw4gx3ygddE5DGgB/A08KKNmDPG+G3lop1syi1i0vVD6T+qu9/pRFxr7xFlAyvcRzLwY/f5T4B+ONcO9cHZK9njeXzN08ajwN9x9pyW4IyQu1pV6z0xtwJf4hSche7z2wIr3dgrgcNuG2+4bX63JTtrjDHh2rVuH8v+vomhE3oy4dIBjb+hDRDP4DLTRNnZ2Zqbm+t3GsaYNubQvir+8tPPSU7ryE3fn0jHpKgcT9YsIrJcVbNDrYuqc0TGGNNe1dU4gxMa6hq44sGxbaoINcYKkTHG+ExV+ehP6ynecYipd42ia69OfqfUqqwQGWOMz/JzCli3rJDsKwcxeFz7u2DeCpExxvhoz+aDLP7LRgaO7cFZVw72Ox1fWCEyxhifVBysZsHM1aR2T2LaXaOQuLZ554TGWCEyxhgf1Nc1sGBmHjVVdVzx4FgSOyX4nZJvrBAZY4wPlvx1I4VbDnLx7Vn06Jvqdzq+skJkjDGtbN2yPaz+qIDx0wYwPLuX3+n4zgqRMca0oqLtZXz4+nr6jezGpOuG+J1OVLBCZIwxraTyUA3zZ64muXMCl947mrh4+xMMVoiMMaZVNNQ3sPClfCrLapn+wFiSUzv6nVLUsEJkjDGtYNk/t1Cwfj8X3DKCjIGd/U4nqlghMsaYCNuYu5eV7+1gzAV9yTq3t9/pRB0rRMYYE0GlBeX8a/Y6eg/twuSvDPc7nahkhcgYYyKk+nAt819YTcekeC67fwzxHexPbij2WzHGmAjQBuW936/h0L4qLr9/LCldEv1OKWqFVYhEZIuI9AixvKuIbGm5tIwxJrZ99vZWtueVMvkrw+k9tIvf6US1cPeIBgHxIZYnAn0be7OITBGRuSJSICIqIncGrRcReVJEdotIpYh8KCKjg2ISReQ5ESkRkQq3vX5BMd1EZLaIHHQfs0Wka1DMABF5y22jRESeFREbT2mMOWVbV5WQ+/Y2Rk7KZMwFjf5pbPeaNAWgiNzgeXmliBz0vI4HLgG2NaGpVCAPeM19BHsceAy4E1gP/Ah4T0RGqOohN2YGcC1wM1AK/AqYJyITVbXejZkDDACmAwq8BMwGrnb7Ew+87b7/fKAH8CogwDeb0A9jjAnpwN7DLPp9Pj0HpHHBLSMQaZ931A5HU+ei/V/3pwIvB62rxSlCjzXWiKq+A7wDICKveNeJs7UeAX6mqn9zl90BFAG3ADNFpAtwD3CXqr7nxtwGbAemAgtFJAu4HJisqkvdmAeAxW5BWw9cCowGBqrqTjfmceAlEXlCVcua9msxxpijaqrqeOeF1cR1iGP6g2PpkBDqAJIJ1qRDc6oap6pxwA4gI/DafSSq6ghVnXeKuQwGMoF3PZ9bCeQA57qLJgIJQTE7gbWemElAObDU0/YSoCIoZm2gCLkW4hxinHiK/TDGtEOqyr9eW8uBwgouu3c0ad2T/E4pZoR1jkhVB6tqSYRyyXR/7g1avtezLhOoB4JzCI4pVlUNrHSfFwXFBH9Oidt2JiGIyP0ikisiucXFxU3qkDGm/Vjx7g42f1HMpBuG0W9kd7/TiSlNPTR3hIicjXNOKIOgQqaq32qBnDTotYRYdlxaQTGh4psSc8LlqjoLmAWQnZ3dWD7GmHZk55p9fPLPzQzLzmD81P5+pxNzwipEIvJd4ClgE7Cbpv1hb6pC92cm4D1klsHRvZdCnMER6UBxUEyOJyZDRCSwV+Sef+oZ1M55QZ+f7rYdvKdkjDEnVFZSycKX8+jWO4WLb8uywQnNEO7w7W8D31LV01T1QlW9yPO4+BRz2YpTIKYFFohIEs6otsD5nuU4gyO8Mf2ALE/MMpzReZM8bU8CUoJisoKGfU8Dqt3PMMaYRtXV1DN/5mpQmP7gWBISbXBCc4R7aK4z7qi35hCRVGCY+zIOGCAi44F9qrpDRGYAT4jIOmAD8EOcgQdzAFT1oIi8DDwtIkUcHb69CljkxqwVkQU4o+zuwzkkNxOY546YA2ewQz7wmog8hjN8+2ngRRsxZ4xpClXlw9fXU7KrnCsfOp2uGZ38TilmhbtH9CecodHNlQ2scB/JwI/d5z9x1z+FU1h+C+QCvYFLPdcQATwK/B14A2c0XDlwtecaIoBbgS9xCs5C9/ltgZVu7JXAYbeNN9w2v3sKfTPGtCOrP9zF+k8LOeuqwQwam+53OjFNPIPLGg8WeQLnWp93cfZCar3rVfVXLZlctMrOztbc3Fy/0zDG+GTnun3Me/ZLBozpwRUPjkXi7LxQY0Rkuapmh1oX7qG5e3H2QM7l6DU5AYqzN2OMMW3W9rxS5s9cTdfMTky9a5QVoRYQViFS1cGRSsQYY6Ld5hVFvPtSPj36pnL1t8aRmBz2FTAmhGb/FkWkF86Fow0tmI8xxkSl9Z8W8v6ra+k1qDNXfdOKUEsKdxqIBBF5SkQOAQU4d+NGRH4uIg9FID9jjPFdXk4Bi15ZQ5/hXW1PKALCHTX3nzh3sP43nGtuAj7DuWO2Mca0KSsX7eCjOesZOKYHVz18Oh2TrAi1tHB/ozcDd6vqRyLiPSSXB5zWcmkZY4y/VJXcd7bx2VtbGXpGBtPuHmVTfUdIuIWoD86UC6Hasf8mGGPaBFVl2T82s+LdHYw8J5OLbhtJXLwVoUgJt3jkA1M4fhK8r2K3xjHGtAHaoOS8sYG8jwoYc0FfpnztNBuiHWHhFqIfA38Ukf44Nwj9ioiMxJm47sqWTs4YY1pTQ30DH8xex7pPCpkwbQCTbhhqNzFtBeHOR/QWzt7PpUADzuCF4Ti32FnU8ukZY0zrqK9r4N2X17Duk0LOunqwFaFWFPZ5HVVdiHP/NmOMaRPqautZOCuPbatLOffGYUyYNsDvlNqVU7mgtSvHT4y371QTMsaY1lRTVcc7v1tNwYb9XHDLCMZM6et3Su1OuBPjDQReAC4CEryrcO41Z5NxGGNiRnVlHfOe+5K9Ww8y9Y4sRpzT2++U2qVw94j+AHQF7ub4GVqNMSZmVJbX8NazX1JaUM5l941h6BkZfqfUboVbiM4CzlHVvEgkY4wxraHiYDVzn1nJweJKpj841uYT8lm4hWgrkBiJRIwxpjUc2lfFm79eQUVZDVc9PI5+I7r5nVK7F+6lwt8G/kdEhjUaaYwxUebA3sP8/RfLqSyv5dpvj7ciFCXC3SN6E2ePaL2IVAN13pWq2rmlEjPGmJZUurucuTNW0tCgXPfoBHoOSPM7JeMKtxA9HJEsXCISDzyJc3fv3sAe4HXgSVWtc2ME50La+4FuwKfAN1Q139NOIvALnJu0JgPvAw+p6i5PTDfgWeAad9Fc4JuqeiByPTTG+KF4xyHmPrOSuA7C9d85g+59UvxOyXiEO0Prq5FKxPV94BvAHcBq4HTgVZwpJ/6fG/M48BjOtBPrgR8B74nICFU95MbMAK7FKUSlOFOYzxORiapa78bMAQYA03FG/70EzMaZ5sIY00bs2XyQec+tJLFTAtc+Op4uPTv5nZIJcioXtGYCHb3LVHXHKeZzLvCWeyshgG0iMhc42/1MAR4Bfqaqf3OX3QEU4dzvbqaIdAHuAe5S1ffcmNtw7ho+FVgoIlnA5cBkVV3qxjwALHYL2vpT7IcxJgrsXLePd55fRUrXRK59ZAJp3ZP8TsmEEO4MrV1E5FURqcSZoXVr0ONUfQxc5N5IFREZBVwMvOOuHwxkAu8G3qCqlUAOThEDmIhzsa03Ziew1hMzCSgHlno+ewlQ4YkxxsSwbatLePs3q+icnsz1j51hRSiKhTtq7hfAOOA6oApnL+R7wC7gay2Qz89xDo+tEZFanGknXlXV5931me7PvUHv2+tZlwnUAyWNxBSr6pELct3nRZ6YY4jI/SKSKyK5xcXFYXfMGNN6Ni0vYv7vVtOjbwrXf+cMUrrYVSfRLNxDc9OBm1V1sYjUA8tV9Q0R2QM8APzvKebzNeB2nAKXD4wHnhGRrar6sicu+I4OEmJZsOCYUPEnbEdVZwGzALKzs+2OEsZEqXXL9vCv19aSObQLV35jHInJNmdntAt3j6grR2doPQj0cJ8vo2UOaT0N/EJV/6yqq1V1Ns5Agx+46wvdn8F7LRkc3UsqxLnnXfCl0sExGeK5x7v7vCfH720ZY2LE6g938f6ra+k7ohtXf3O8FaEYEW4h2gwMcZ+vBb7u/gG/AWiJO293wjms5lXP0Ty34hSRaYGVIpIEnM/R8z3LgdqgmH5AlidmGZCKc64oYBKQwrHnjYwxMeKLd7eT8+cNDDo9nSu/cToJiXYP5lgR7n8XXsEZUv0h8DNgHs61RXE4d104VW8B/0dEtuIcmpsAfAd4DZzzOCIyA3hCRNYBG4Af4gw8mOPGHBSRl4GnRaSIo8O3VwGL3Ji1IrIAZ5TdfTiH5GYC82zEnDGxRVX5bN5Wct/exrDsDKbeNYr4+HD/j238FO51RL/2PP+XO7otG9ioqqtbIJ9v4lwv9DzOobQ9wIvATzwxT+FcpPpbjl7QeqnnGiKAR3Hu+vAGRy9ovd1zDRHArTgXtAZG180lwhfsGmNalqqy9G+bWLloJyPP7c1F/zaSuDibVTXWiGfgmGmi7Oxszc3N9TsNY9o1bVA++vMG8nMKGHtRP87/ynDEilDUEpHlqpodal3YZ/JEZALOxHgZHD9D6+PNytAYY8LQUN/Av15bx/pPCznjsoGcc90QPGOPTIwJd4bWx3HODW3HGV3W2HBoY4xpUfV1Dbz3cj6bVxRz9jVDyL5ikN8pmVMU7h7Ro8C/q+rMSCRjjDEnU1dTz/yZeezIL2XyV4Yz7pL+fqdkWkC4hSgO58S/Mca0qpqqOt55fhUFGw9w4a0jGH1+X79TMi0k3DGOvwPuikQixhhzItWHa5n7zEp2bzrI1DtHWRFqY8LdI/ox8I6IrMSZpqHWu1JV726hvIwxBoDKQzXMfXYl+/ZUcPl9YxgyoaffKZkWFm4h+m/gUuALnGt4bICCMSZiKg5U8+aMFZSVVnHlv5/OgNE9Gn+TiTnhFqKHgFtU9Y1IJGOMMQFlJZW8OWMFlYdqufqb4+h7Wje/UzIREm4hqgRWRCIRY4wJOLD3MG/OWEFtdT3XPDKezMFd/E7JRFC4gxV+DTwiduWYMSZCSgvK+fsvv6C+roHrvjPBilA7EO4e0fnAFOBKEVnD8YMVrmmpxIwx7U/R9jLmPruSDh3iuPY7Z9AtM8XvlEwrCLcQlQB/j0Qixpj2bfemA8z7zZckpSRw7SMT6NIz2e+UTCtpciESkQ7A58A/VXV35FIyxrQ3O9fs453frSK1exLXPjKe1G5JfqdkWlGTzxGpah3ODKoJkUvHGNPebF1Vwrznv6RLRieuf+wMK0LtULiDFT4BJkYiEWNM+7Mxdy8LXlhNet9UrvvOBDp17uh3SsYH4Z4jehH4hYgMwJmSu8K7UlW/aKnEjDFt25olu/ngj+voM6wrVz50Oh2Tw56VxrQR4W75Oe7PX4VYp4BNEm+MadSqD3ax+I0N9B/VnekPjiWho/3paM/CLUSDI5KFMabdWL5gG5/8cwuDx6Vz2b1jiE8I9wyBaWvC+hegqttP9miJhESkt4i8KiLFIlIlImtE5ALPehGRJ0Vkt4hUisiHIjI6qI1EEXlOREpEpEJE5opIv6CYbiIyW0QOuo/ZItK1JfpgjDmeqvLJm5v55J9bGH5mLy6734qQcYT9r0BEpovIPLdA9HeX3Ssil5xqMm4hWAIIcCWQBXwTKPKEPQ485i4/0133noikeWJmADcCN+NchNsZmCci3v3/OcAZwHTgcvf57FPtgzHmeKrKkr9uYvn87Yw6rzdT7xpFfLwVIeMId6rwW4EXgJeASzg6lDsep0Cc6qR5jwN7VPV2z7Ktns8X4BHgZ6r6N3fZHTjF6BZgpoh0Ae4B7lLV99yY23CmN58KLBSRLJziM1lVl7oxDwCLRWSEqq4/xX4YY1wNDcpHc9az5uPdnH5xPyZ/ZTh2lzDjFe5/SR4H7lPVR4E6z/JPgPEtkM91wKci8oaIFInIShF52HNvu8FAJvBu4A2qWgnkAOe6iybiFEhvzE5grSdmElAOLPV89hKcUYDnYoxpEQ31DSz6wxrWfLybidMHWhEyIYVbiIYDy0IsL8c5/HWqhuBMNbEFuAx4BvgZ8A13fab7c2/Q+/Z61mUC9Ti3IzpZTLGqHplPyX1e5Ik5hojcLyK5IpJbXFwcZreMaX/qaxtY+GI+Gz/fyznXDeGca4daETIhhVuIdgOnhVg+Bdh86ukQB3yhqj9Q1RWq+gfgWY4WooDgCfkkxLJgwTGh4k/YjqrOUtVsVc3u2dNmiDTmZGpr6nnnd6vYsrKYyV8dzsTLB/mdkoli4RaiWcCzInKe+7q/e47mKeB3LZDPHmBN0LK1wAD3eaH7M3ivJYOje0mFOOes0huJyfBOZ+E+78nxe1vGmDDUVNUx77kv2bF2HxfdNpJxF/f3OyUT5cIdvv0U8DfgPSAF+ABn8MILqvrbFshnCTAiaNlpOAMNwBm4UAhMC6wUkSSckXGB8z3Lcaan8Mb0wxmBF4hZBqTinCsKmITTJ+95I2NMGKoqanlzxkr2bD7IpXePZtR5ffxOycSAJo2aE5FOODc8vQ5nIMBbwC/d1WtUtbyF8vk1sFREngDeACYA3wL+A5zzOCIyA3hCRNYBG4Af4pyjmuPGHBSRl4GnRaQIKMW5E8QqYJEbs1ZEFuCMsrsP55DcTGCejZgzpnkOl9Uw95mV7N9bweX3j2HIeDuEbZqmqcO3fwzcCbyOM134LUCcqn6lJZNR1c9F5Drgp8D/BXa4P5/3hD0FJAO/BboBnwKXquohT0xgVN8bbuz7wO2qWu+JuRXn/FNgdN1c4OGW7I8x7UX5/mrenLGC8n1VXPXQOPqP6u53SiaGiGfg2ImDRDYDT6jqn93XZ+EcRksK+uPeLmRnZ2tubq7faRgTFcpKKnlzxgoqy2u56hvj6DO8q98pmSgkIstVNTvUuqaeI+oPLA68UNXPcPY47ACwMe3Y/sIK/v6LL6g+XMe1j0ywImSapamH5uKBmqBldWG83xjTxpTsKmfuMysAuO47Z5DeL9XnjEysamohEeCPIlLtWZYEvCgihwMLVPWalkzOGBOd9m4t463nVpKQGM+1j0yga69OfqdkYlhTC9GrIZb9sSUTMcbEht0b9zPvN6tITkvg2kcm0Dk92e+UTIxrUiFS1bsinYgxJvrtyC9l/gurSeuRxDXfnkBqt0S/UzJtgJ3jMcY0yZaVxSx8KY/uvVO45lvjSU7r6HdKpo2wQmSMadSGzwpZ9MpaMgamcdXD40hKSWj8TcY0kRUiY8xJrfl4Nx+8vo6+w7tyxUOn0zHJ/myYlmX/oowxJ/Tl+zv5+K8bGTC6B9MfGEOHjvGNv8mYMFkhMsaElPvONj6du4UhE3py6d2jiU+wqb1NZFghMsYcQ1X55M0tfLFgO6ed3YtLbs8iLt6KkIkcK0TGmCO0QVn8142s/mAXo8/vwwU3j0DibFZVE1lWiIwxADQ0KB/+cR1rl+5h3CX9Oe+mYTa1t2kVVoiMMdTXN/D+H9awMbeI7CsHcdZVg60ImVZjhciYdq6utp6FL+azbVUJk64fyhmXDfQ7JdPOWCEyph2rra5n/gur2Ll2P1O+fhpjL+znd0qmHbJCZEw7VVNZx7zffknh5oNcfHsWWef29jsl005ZITKmHaoqr+Wt51ZSsrOcafeMZnh2L79TMu1YVF8cICL/ISIqIr/xLBMReVJEdotIpYh8KCKjg96XKCLPiUiJiFSIyFwR6RcU001EZovIQfcxW0S6tlLXjPFFTVUdxTsO8Y9ffUFpQQXTHxxrRcj4Lmr3iETkHOA+YFXQqseBx4A7gfXAj4D3RGSEqh5yY2YA1wI3A6XAr4B5IjJRVevdmDnAAGA6oMBLwGzg6gh1yZiIq6up59C+KspKqzhUUklZaRVlJVUcKnWeV5XXAtChYxxXPnw6/Ud29zljY6K0EIlIF+B14B6cQhNYLsAjwM9U9W/usjuAIuAWYKb73nuAu1T1PTfmNmA7MBVYKCJZwOXAZFVd6sY8ACx2C9r6VumoMWGqr2vg0L4qDpVUUeYWl0OlVZSVVHKotIrDZTXHxMd1ENK6J9E5PZkhA9Lo3COJzj2S6TW4s01oZ6JGVBYiYBbwv6r6LxH5kWf5YCATeDewQFUrRSQHOBeYCUwEEoJidorIWjdmITAJKAeWetpeAlS4MVaIjC8a6hso31/tFphKd2/GKTqHSqsoP1Dt7L+74uKE1O6JpPVIZuCYHnROTyKtRzKdezg/U7p0tDsjmKgXdYVIRO4DhgG3hVid6f7cG7R8L9DXE1MPlISIyfTEFKvqka+0qqqIFHligvO6H7gfYMCAAU3qizHBGhqUigPVxxSXwN5MWYlTaLThaKURgZRuiXTukUy/Ed1I6+Hs3QR+pnTpaPeBMzEvqgqRiIwAfgqcr6o1JwnVoNcSYtlxzQfFhIo/YTuqOgtnT43s7OzGPsu0U9qgHD5Uc8x5mSPnakqrKN9XRUP9sf98Urp0pHN6Mr2HdTm20PRIJrVbIvEdrNCYti2qChHOIbN0IM9ze5F4YIqIPAgERsdlAjs978vg6F5SofuedKA4KCbHE5MhIhLYK3LPP/Xk+L0tY45QVarKaykr8ezReIrNoX1V1Nc2HPOe5LQEOqcnkzEwjWFnZLjFxi003RPpkGBz/Jj2LdoK0T+B3KBlfwA24uwpbcApItOAzwFEJAk4H/ieG78cqHVj5rgx/YAsjp4TWgak4hS+wLJJQArHnjcy7YyqUn247sghsyNFZl/VkaJTV11/zHuSUhJI65FEjz4pDDo93T0/4xSatPQkEmwyOWNOKqoKkaoeAA54l4lIBbBPVfPc1zOAJ0RkHU5h+iHOwIM5bhsHReRl4Gn3nE9g+PYqYJEbs1ZEFuCMsrsP55DcTGCejZhr+2oq69xhzZVB52qcw2k1VccWmo5J8aSlJ9OlZzL9R3YnLT3JGX2Wnkxa9yQ6JkfV18iYmBOL36CngGTgt0A34FPgUs81RACPAnXAG27s+8DtnmuIAG4FnuXo6Lq5wMORTd20htrq+uOKy5FhzqWVVFfUHRPfITHeHdacRJ/hXd2RZ+4eTY8kklISfOqJMe2DeAaOmSbKzs7W3NzgI4imtdTV1nOo9OihssAw58DzykO1x8THJ8Qde7jsmJFnTqGxKQ+MiSwRWa6q2aHWxeIekWmHairrWP9pIfkf76Z0V/kx6+LiAxdtJpE+rucxgwHSeiTRqXNHKzTGRDErRCaqFe84RF5OARs+30tddT09B6Rx5lWD6ZKeRFq6c+Fmpy6JxNlFm8bELCtEJurU1tSzKbeIvJwCiraV0SEhjuFn9mL0lL70GtTZ7/SMMS3MCpGJGvsLK8jP2c26T/ZQfbiObpmdmPzV4Yw8J5PETjZgwJi2ygqR8VV9XQNbVhaTn1NAwYYDxMULQyf0ZPSUvvQZ3tXO7RjTDlghMr4oK6lkzce7WbN0D5VlNaT1SOKc64aQdW4fOnXu6Hd6xphWZIXItJqGBmVHXil5iwvYnleKAAPHpjNmSl/6j+puAw6MaaesEJmIqzhYzdole8j/uIDyfdV06tyR7OmDGDW5D2ndk/xOzxjjMytEJiJUlYL1+8nL2c3WlcU0NCj9RnZj8k3DGTQunXibusAY47JCZFpUVUUt65btIX/xbg7sPUxipw6MvbgfY87vS9denfxOzxgThawQmVOmquzdVkb+RwVsXF5EfW0DmUM6c8mdWQw7I4MOdvdpY8xJWCEyzVZTVceGz/aSv7iAkp3lJCTGM3JSb8ZM6UN6vzS/0zPGxAgrRCZspQXl5H1UwPrPCqmtqqdH31QuuGUEp53Vi45J9k/KGBMe+6thmqSutp7NXxST91EBhVsOEt8hjmHZGYyZ0pdegzvbhafGmGazQmRO6sDew+QvLmDtsj1UV9TRJSOZ824axshzepOUarfdMcacOitE5jj19Q1sW1VC3kcF7Fq3n7g4YfD4dEZP6Uu/07ohduGpMaYFWSEyRxzaV8WaJbtZ8/FuDh+sIbVbImdfM5is8/qQ0iXR7/SMMW2UFaJ2ThuUHWv3kZ9TwLZVJSgwcHQPRt/al4GjuxNnF54aYyIsqgqRiPwAuAEYAVQDnwA/UNU8T4wA/wncD3QDPgW+oar5nphE4BfAzUAy8D7wkKru8sR0A54FrnEXzQW+qaoHItW/aHK4rMa98LSAspIqktMSmHDZQEZP7kPn9GS/0zPGtCNRVYiAC4Hngc8BAX4CLBKRUaq6z415HHgMuBNYD/wIeE9ERqjqITdmBnAtTiEqBX4FzBORiapa78bMAQYA0wEFXgJmA1dHsH++UlX2bDpAXs5uNn9RREO90md4V865bihDxvckvoPt/RhjWp+oqt85nJCIpAIHgetU9S13b2g38BtV/W83JhkoAr6rqjNFpAtQDNylqq+7Mf2B7cB0VV0oIlnAGmCyqi5xYyYDi4GRqrr+ZHllZ2drbm5uJLocEdWVdaz/ZA95ObvZv6eCjskdGHlOJqOn9KV77xS/0zPGtAMislxVs0Oti7Y9omBpQByw3309GMgE3g0EqGqliOQA5wIzgYlAQlDMThFZ68YsBCYB5cBSz2ctASrcmJMWolhRtL2MvJwCNn6+l7qaBjIGpnHRbSMZfmYvEuy2O8aYKBHthegZYCWwzH2d6f7cGxS3F+jriakHSkLEZHpiitWzO6iqKiJFnphjiMj9OOelGDBgQLj9aDW11fVszN1Lfk4BRdsP0aFjHKed2YvRU/qSMbCz3+kZY8xxorYQicivgMk4h8/qg1YHH0+UEMuOazIoJlT8CdtR1VnALHAOzTXyWa1u3+4K8hYXsP6TQmoq6+jeJ4Xzv3YaI87JJDE5ajezMcZEZyESkV8DXwcuUtUtnlWF7s9MYKdneQZH95IKgXggHedckTcmxxOTISIS2Ctyzz/15Pi9rahVX9vA5pVF5OfsZvfGA8R1EIZOcG6703tYF7vtjjEmJkRdIRKRZ3CK0IWqui5o9VacIjINZ2QdIpIEnA98z41ZDtS6MXPcmH5AFkfPCS0DUnHOFQWWTQJSOPa8UVQqK6l0bruzdA+Vh2rpnJ7EpOuHknVub5LTOvqdnjHGhCWqCpGI/Ba4DbgO2C8igfM15apa7p7HmQE8ISLrgA3AD3EGHswBUNWDIvIy8LR7zicwfHsVsMiNWSsiC4CZInIfziG5mcC8xkbM+aWhvoHteaXk5RSwY80+BBh0ejpjpvSlf1Z3u+2OMSZmRVUhAh5yf74ftPzHwJPu86dwLlL9LUcvaL3Ucw0RwKNAHfAGRy9ovT3oXNOtOBe0BkbXzQUebpFetKCKA9VHbrtTvr+alC4dOfOKQYya3IfUbkl+p2eMMacsqq8jilaRvo5IG5Rd6/eTl1PA1i9L0Aalf1Y3xkzpx8DTexBvt90xxsSYWL6OqF2pKq9lrXvbnYNFlSSlJDDukv6MPr8PXTM6+Z2eMcZEhBUin6kqhVvKyM8pYNPyIurrGug9tAtnXjmYoWf0pEOCXXhqjGnbrBD5pKayjg2fFZKXs5vSgnISkuLJOq83Y6b0pUffVL/TM8aYVmOFqJUV7zxEfk4BGz7bS211Pen9U7nw1hEMP7MXHZNscxhj2h/7y9dKaqrqmPvMSvZuLSM+IY7h2RmMmdKPjEFpduGpMaZds0LUSjomdaBLz2SGZ/dixDmZJKUk+J2SMcZEBStErWja3aP9TsEYY6KOXZBijDHGV1aIjDHG+MoKkTHGGF9ZITLGGOMrK0TGGGN8ZYXIGGOMr6wQGWOM8ZUVImOMMb6y+YiaQUSKge3NfHs6UNKC6fiprfSlrfQDrC/RqK30A06tLwNVtWeoFVaIWpmI5J5ocqhY01b60lb6AdaXaNRW+gGR64sdmjPGGOMrK0TGGGN8ZYWo9c3yO4EW1Fb60lb6AdaXaNRW+gER6oudIzLGGOMr2yMyxhjjKytExhhjfGWFyBhjjK+sEJ0CEXlIRLaKSJWILBeR808S+6SI6AkeGZ64C9y2qkRki4g8GIt9EZELT7B+ZDT1xY2/TESWicghESkRkTdF5LSgmFbfLi3djxjbJl8VkZUiclhEtovI90LERP13xY0/aV/82C4iMkVE5opIgftZdzbhPWNF5CMRqXTf9yMRkaCY5m0TVbVHMx7A14Ba4D4gC3gOKAcGnCA+FcgMenwIfOCJGQxUuG1luW3XAjfGYF8uBBQYFRQXH2V9GQxUAU8Bw4DxwLvAJj+3S4T6ESvbZDpQBzwEDAGuBHYDD8fgd6UpfWn17QJcAfwUuAk4DNzZSHxnoBD4CzAGuBE4BDzWEtskYhusrT+AT4EXg5ZtBP6nie/vD9QDt3iW/RzYGBT3ErAsBvsS+HKlR/N2cb+I9d4vPXCRN3c/tkuE+hEr22QO8I+gZd8EdnJ0pG9MfFea2Bdftosnn/ImFKJ/B8qAZM+yHwIFLbFN7NBcM4hIR2Aizv84vd4Fzm1iM/cAB4C/eZZNCtHmQiBbRBLCz7RxEexLQK6I7BGR90XkomYn2gTN7Esuzv/a7hWReBFJA+4APlfVwD21WnW7RLAfR2KjfJsk4uzdeVUC/YCB7utY+a40pS8BrbZdmmESsFhVKz3LFgJ9gEGemGZtEytEzZMOxAN7g5bvxdmlPikRiQPuBl5T1WrPqswTtNnB/cxIiFRf9uD8L+pG4AZgPfC+iExpiaRPIOy+qOo2YBrwY6AaOAiMBa7yhLX2dolUP2Jim+D88bpORC4VkTj3PNdj7rre7s9Y+a40pS9+bJdwnej3HVh3sphGt0mHU82unQu+GlhCLAtlOs7hrJea2Gao5S2tRfuiqutxvlABy0RkEPBdIKf5aTZJk/siIpnAy8BrwJ+ANOAnwF9E5GJVbThJm6GWt6QW7UesbBPgRWAo8CaQgHNI6BngSZzDjydrM9TyltaiffF5u4SjKb/vZm0T2yNqnhKcf0TB/wvK4Pj/EYRyP7BUVfODlheeoM06oLQZeTZFpPoSyqfA8PDSC0tz+vINoEJVH1fVFaqaA/wbcAFHD7e09naJVD9Cibptoo7v4wyKGei+9zN39Tb3Z0x8V5rYl1AivV3CdaLfNxzte7O3iRWiZlDVGmA5zqEQr2nA0pO9V0T64IyceTHE6mXA1BBt5qpqbfOyPbkI9iWU8TiHISKimX3pxLH/y8bzOvD9aNXtEsF+hDKe6NsmgffWq2qB28bNOCe9i9zVMfVdaaQvoYwngtulGZYB54tIkmfZNJwRgNs8Mc3bJn6M0mgLD5xhnDXAvThDFZ/BGX0y0F3/P8D7Id73Q5zj951CrAsMf5zhtnmv+xmtMSS1pfvyCHAdzv/qRrttKHBDNPUFuBhoAP7TzfUMYAGwA0jxa7tEqB+xsk3Scc6ZZOH8QX4G5wT/WbH2XWliX1p9u+DsoY13H4eBH7nPB5ygH11w9nj+jDN8+wacw4yhhm+HvU0itsHawwPn2oBtOCeHlwNTPOteAbYFxQuwFXj+JG1eAHzhtrkVeDAW+wI8Dmxyv3T7gMXAFVHal6+7v/NyoBh4Cxjl93Zp6X7EyjbB+eO9zO1HBbAIODtEm1H/XWlKX/zYLhwdMh78eOUk/77G4pyzqsLZW/tP3KHbp7pN7O7bxhhjfGXniIwxxvjKCpExxhhfWSEyxhjjKytExhhjfGWFyBhjjK+sEBljjPGVFSJjmkBEypsyeVi08ky+dko3BBWRD0XkN+HENPbaGCtEpl0RkVdOMBvmJ37n1obcAPygqetFZJuIfDfiWZmoZXffNu3RIuC2oGU1fiTSEkQkQSN0f7XmUNV9p7LetD+2R2Tao2pVLQx6HPnjKCLD3MNHVSKyXkSuCm5ARM4WkS/cmBUicoW7Z3WhJ2aUiLwtIodEpEhE/uRO1xCSiAxy27hFRD52214nIpd6YgKH2K4Qkc9EpAa4TEQSRWSGiOx13/eJiEwO8THniMhKN2a5iEz0tN3DzXGXiFSKSL6I3BWijQ4i8oyI7HcfT7vzUgXaOemhN+96EfkQ567UT3v2TlNEpExEbgp63zQRqRWRXidq28QmK0TGeLh/UP+B892YhDPp35M4M20GYlKBecA6nBk7HweeDmqnN859ufKAs3DuSpwKzPX+0T6Bp4BncW5C+R7wpoj0DYr5Oc5NZ0fiTBnwFM4NOe8GJgCrgQVuHl6/AL4PZANbgLdFpJO7LgnnPmFX4dx88xlgpohcEtTGrRz9/TyAMxXII4306URuAHbhzJ3UG+itqhU4cyrdHRR7NzBPVZsyPYmJJa1xk0B72CNaHjg3c6zDuQml9/Fzd/2lONMnDPC8ZzLODSHvdF8/gHNzymRPzC1uzIXu658QdMdyoJsbc9YJchvkrn/CsywO2AD8l/v6QjfmRk9MCs6hxds9y+KBzSHed6snJhVnivd7T/L7+jPwkuf1h24+4ln2Q2BXUMxvwni9Dfhu0Odmu9upr+d3Vwlc5fe/IXu0/MP2iEx7lMPRW+AHHoE9miygQFV3eOI/xZliIWAkkKeqlUExXhOBKe5ou3IRKQd2uuuGNpLfssATdWaI/RQYFRST63k+FGf2zyWe99W77QS/z9t2Oc6e0ygAEYkXkSdEZJWIlLo53wAMCGrjE3Wrg6fNviLSuZF+NZmq5rq53eEuugXYD8xvqc8w0cMGK5j26LCqbjrBOjnB8uCYxm5bHwe8jTPdc7CWOLRUEZQPJ8gpnNvrfxd4DPg2ThEoB37K0Zk4W9tLOIf8fopzWO4Vt8CaNsb2iIw51hqc/9339yw7i2O/K2uBsSKSHBTj9QXOeZbtqrop6HGokRzOCTwREXHbXnuS+E04h+aODE4QkXicczhrTtJ2Cs4kZ4G2JwNvqepsVV2Jc2jvtBCfd7abl7fN3apadvJunVANzqHEYH/E2RYP40z094dmtm+inBUi0x4likhm0KOnu24RziCE10RkvIhMAn6Nc74i4HWc80gvuiPjpgL/4a4L7IH8FmdWyzfcEXZDRGSqiMwSkbRG8vt3EblJREbgzHY5EPjdiYLVObn/O+Bn7mi6LPd1L+D5oPAfuqPPRgO/xykCc9x1G4BLRGSyiIwEfoMz62awPsAMERnhjmz7Hs7vqLm24UxD3Vc8F9yq6kHgr8AvgRxV3XgKn2GimBUi0x5NxZlh0vtYAUfOyVyP8934FHgN+C+cGSdxY8qBq3H2eFbgnF960l1d5cbsBs7DObe0AMjHKU7V3rZO4P8A3wG+BC4HrlfVXY285/vAX3D2GlYCpwOXq+qeEG3/EmePbTjOyf/AYb7/Aj7DOQ+Tg3P47/UQn/U6zh7Mp8CLwMucWiH6EdAfZw+sOGjdy0BH96dpo2yGVmNagIhcizPsO0NVS5rZxiCc6ZXPdE/Wt3si8jVgJtBHVQ/7nY+JDBusYEwziMgdONfh7MQ5zzID5/xKs4qQOZZ7bdMgnEOeL1oRatvs0JwxzdMLmA2sxznkNh/4N18zalsexzk0uQ/4fz7nYiLMDs0ZY4zxle0RGWOM8ZUVImOMMb6yQmSMMcZXVoiMMcb4ygqRMcYYX/1/d15GstTl6ogAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(probs, exactList2, label=\"Size 8\", color='tab:purple')\n",
    "plt.yticks(fontsize = 14)\n",
    "plt.xticks(fontsize = 14)\n",
    "\n",
    "# plt.legend()\n",
    "plt.xlabel('Edge probability', fontsize = 14)\n",
    "plt.ylabel('Permanent', fontsize = 14)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The two plots above reproduce figures $4$ and $5$ of the paper. We see the mean value of the permanent of $15$ randomly generated graphs of $8$ (or $7$, respectively) vertices, plotted in function of edge probability. The edge probability represents the probability that any two vertices $i$ and $j$ of the randomly generated graph are connected by an edge.\n",
    "\n",
    "For the last code in this notebook, you can use it to reproduce figure 6 from the paper where we plot the mean value of estimation and calculation of the permanent for random graphs of sizes 6."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.7  0.78 0.86 0.94 1.  ]\n",
      "i in probs: 0.7\n",
      "Total number of samples:  6127000\n",
      "Number of samples post: 500\n",
      "Total number of samples:  29852000\n",
      "Number of samples post: 500\n",
      "Total number of samples:  1273000\n",
      "Number of samples post: 500\n",
      "Total number of samples:  2224000\n",
      "Number of samples post: 500\n",
      "i in probs: 0.7799999999999999\n",
      "Total number of samples:  1671000\n",
      "Number of samples post: 500\n",
      "Total number of samples:  23356500\n",
      "Number of samples post: 500\n",
      "Total number of samples:  3599500\n",
      "Number of samples post: 500\n",
      "Total number of samples:  2005500\n",
      "Number of samples post: 500\n",
      "i in probs: 0.8599999999999999\n",
      "Total number of samples:  1688000\n",
      "Number of samples post: 500\n",
      "Total number of samples:  1585500\n",
      "Number of samples post: 502\n",
      "Total number of samples:  1985000\n",
      "Number of samples post: 500\n",
      "Total number of samples:  1750500\n",
      "Number of samples post: 500\n",
      "i in probs: 0.9399999999999998\n",
      "Total number of samples:  2090000\n",
      "Number of samples post: 500\n",
      "Total number of samples:  2010000\n",
      "Number of samples post: 500\n",
      "Total number of samples:  1720500\n",
      "Number of samples post: 500\n",
      "Total number of samples:  1765500\n",
      "Number of samples post: 500\n",
      "i in probs: 1.0\n",
      "Total number of samples:  1817000\n",
      "Number of samples post: 500\n",
      "Total number of samples:  1706500\n",
      "Number of samples post: 500\n",
      "Total number of samples:  1779500\n",
      "Number of samples post: 500\n",
      "Total number of samples:  1820500\n",
      "Number of samples post: 500\n",
      "exact_list [(16.75+0j), (102.5+0j), (192.5+0j), (216.5+0j), (265+0j)]\n",
      "estimation_list [16.823075838408243, 103.8843088917408, 198.73325236239438, 216.68009650589966, 261.87748694919446]\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/var/folders/8v/r85qsjj96qn5kpsdfhrz1g2h0000gp/T/ipykernel_22833/2655738150.py:19: UserWarning: color is redundantly defined by the 'color' keyword argument and the fmt string \"ob\" (-> color='b'). The keyword argument will take precedence.\n",
      "  plt.plot(probs, exact_list, \"ob\", label=\"Exact Value\", color='C0')\n",
      "/Users/filipagoncalvescarvalho/opt/anaconda3/lib/python3.9/site-packages/matplotlib/cbook/__init__.py:1298: ComplexWarning: Casting complex values to real discards the imaginary part\n",
      "  return np.asarray(x, float)\n",
      "/var/folders/8v/r85qsjj96qn5kpsdfhrz1g2h0000gp/T/ipykernel_22833/2655738150.py:20: UserWarning: color is redundantly defined by the 'color' keyword argument and the fmt string \"ob\" (-> color='b'). The keyword argument will take precedence.\n",
      "  plt.plot(probs, estimation_list, \"ob\",\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAEGCAYAAACKB4k+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAhl0lEQVR4nO3de3hU5bn38e9NiCRCAEGUY42lGAyCwQYs0oLKhqBFwFI5lLpFq9ZqqXVXWrButVbf0q3VvrtYW6y+0G4U4gFwd/egHMStIhA0yhmloOWkAUw4CJiE+/1jVpYTSMiEZDI5/D7XNdfMPGutWfeTBfllHeZZ5u6IiIgANEt0ASIiUn8oFEREJKRQEBGRkEJBRERCCgUREQk1T3QBNXHmmWd6enp6ossQEWlQVq9evcfdO1Q0rUGHQnp6Onl5eYkuQ0SkQTGzDyqbpsNHIiISUiiIiEhIoSAiIqEGfU6hIsXFxWzfvp0jR44kuhSpZSkpKXTt2pXk5ORElyLSaDW6UNi+fTtpaWmkp6djZokuR2qJu7N37162b9/Oueeem+hyRBqtRhcKR44cUSA0QmZG+/btKSgoSHQpIgm14O0dPPT3TewsPEzntqlMyclgdN8utfb5jS4UAAVCI6XtKk3dgrd3MO2FNRwuLgVgR+Fhpr2wBqDWgkEnmkVEGoiH/r4pDIQyh4tLeejvm2ptHQqFOEhKSiIrKyt8TJ8+vdY+Oz8/n7/85S8ntB86dIj27dtTVFRUrn306NHk5uZW+nmtWrWqtdpEJL52Fh6m59Ekbi5qwZ2FKdxc1IKeR5PYWXi41tbRKA8fVUc8js+lpqaSn59fOwUeJz8/n7y8PK688spy7S1btmTYsGEsWLCA6667DoCioiJee+01nn766bjUIiJ1a0DzVPoXOslEDqW2cWP44WTatTyt1tbRpPcUyo7P7Sg8jPP58bkFb++o9XUVFRWRkZHBpk2R3bwJEybwxBNPAPC9732P7OxsevXqxb333hsus2rVKi655BIuvPBC+vfvT1FREffccw/z5s0jKyuLefPmlVvHhAkTmDt3bvh+/vz5DB8+nGPHjjFkyBAuuugievfuzcKFC0+o75VXXmHEiBHh++9///vMmjULgNWrVzN48GC+/OUvk5OTw65du2rt5yIisRt0pHkYCGWSMQYdqb2/75v0nsLJjs/VZG/h8OHDZGVlhe+nTZvGuHHjmDFjBpMmTeL222/nk08+4aabbgLgwQcfpF27dpSWljJkyBDeffddevbsybhx45g3bx79+vVj//79nH766dx///3k5eUxY8aME9Y7fPhwbrzxRvbu3Uv79u2ZO3cukydPJiUlhfnz59O6dWv27NnDV77yFUaOHBnTidvi4mImT57MwoUL6dChA/PmzeOnP/0pTz311Cn/fETk1JQeLKlW+6lo0qFQ2XG4mh6fq+zw0dChQ3n22We57bbbeOedd8L23NxcZs6cSUlJCbt27WL9+vWYGZ06daJfv34AtG7dusr1nnbaaYwcOZLnnnuOMWPGkJ+fz7Bhw3B37rrrLl599VWaNWvGjh07+Oijj+jYsWOVn7lp0ybWrl3L0KFDASgtLaVTp04x/iREpDa1ateCg/uOVtheW5p0KHRum8qOCgKgc9vUuKzv2LFjbNiwgdTUVPbt20fXrl3ZunUrDz/8MKtWreKMM85g0qRJHDlyBHc/pUswJ0yYwAMPPIC7M2rUKJKTk5k1axYFBQWsXr2a5ORk0tPTT/jGd/PmzTl27Fj4vmy6u9OrVy+WL19es86LSI0NGNWdpXM2UvLZ5/9Xm5/WjAGjutfaOpr0OYUpORmkJieVa0tNTmJKTkZc1vfoo49y/vnn88wzz3DDDTdQXFzM/v37admyJW3atOGjjz7ir3/9KwA9e/Zk586drFq1CoADBw5QUlJCWloaBw4cqHQdl112Ge+99x6PPfYYEyZMACLnM8466yySk5NZunQpH3xw4qi555xzDuvXr+fo0aMUFRWxePFiADIyMigoKAhDobi4mHXr1tXqz0VEYnPexR25bGLPcM+gVbsWXDaxJ+ddXPVef6ya9J5C2XmD2r766PhzCsOHD+eGG27gD3/4AytXriQtLY1BgwbxwAMP8LOf/Yy+ffvSq1cvvvjFLzJw4EAgciho3rx5TJ48mcOHD5OamsqiRYu47LLLmD59OllZWeG5imjNmjVjzJgxPPvsswwaNAiAiRMnctVVV5GdnU1WVhY9e/Y8oeZu3boxduxY+vTpQ48ePejbt29Yx3PPPccPfvADioqKKCkp4Yc//CG9evWq0c9IRE7NeRd3rNUQOJ65e9w+PN6ys7P9+JvsbNiwgfPPPz9BFUm8afuK1JyZrXb37Iqmxe3wkZl1M7OlZrbBzNaZ2e1B+31mtsPM8oPHlVHLTDOz981sk5nlxKs2ERGpWDwPH5UAP3L3t8wsDVhtZi8H0x5194ejZzazTGA80AvoDCwys/Pcvfw1oyIiEjdx21Nw913u/lbw+gCwATjZwfpRwFx3P+ruW4H3gf7xqk9ERE5UJ1cfmVk60BdYETR938zeNbOnzOyMoK0L8M+oxbZTQYiY2c1mlmdmeRpGWUSkdsU9FMysFfA88EN33w88DnQHsoBdwK/KZq1g8RPOgrv7THfPdvfsDh06xKdoEZEmKq6hYGbJRAJhjru/AODuH7l7qbsfA57g80NE24FuUYt3BXbGsz4RESkvnlcfGfAksMHdH4lqjx4j4WpgbfD6RWC8mbUws3OBHsDKeNUXT9UZOnvBggWsX78+fH/PPfewaNGiGtdQWFjIb3/72/D9zp07+eY3v1njz73vvvuYNm1aubb8/PyTXiZ633338fDDD1c6XUTqj3hefTQQuBZYY2b5QdtdwAQzyyJyaGgb8F0Ad19nZrnAeiJXLt1WF1cebV6xm+ULt3Bw31FatWvBgFHda/zFkOoMnb1gwQJGjBhBZmYmAPfff3+N1l2mLBRuvfVWADp37sxzzz1X48+dMGECV1xxBb/4xS/Ctrlz5/Ktb32rxp8tIokXz6uPXnN3c/c+7p4VPP7i7te6e++gfaS774pa5kF37+7uGe7+13jVVmbzit0snbMxHGDq4L6jLJ2zkc0rdsdlfVOnTiUzM5M+ffpw55138sYbb/Diiy8yZcoUsrKy2LJlC5MmTQp/eaenp3PXXXcxYMAAsrOzeeutt8jJyaF79+787ne/i9R88GCFw2JPnTqVLVu2kJWVxZQpU9i2bRsXXHABEBnX6Prrr6d379707duXpUuXAjBr1iy+8Y1vMHz4cHr06MGPf/zjE/qQkZFB27ZtWbFiRdiWm5vL+PHjeeKJJ+jXrx8XXnghY8aM4dNPPz1h+UsvvZSyLxzu2bOH9PR0IDLQ3pQpU+jXrx99+vTh97//fS391EWkOpr0MBfLF24pN7AUQMlnx1i+cEuN9hYqGjp76NChzJ8/n40bN2JmFBYW0rZtW0aOHMmIESMqPbTTrVs3li9fzh133MGkSZN4/fXXOXLkCL169eKWW26pdFjs6dOns3bt2nCPZdu2beFnPvbYYwCsWbOGjRs3MmzYMDZv3gxEDgW9/fbbtGjRgoyMDCZPnky3bt3K1VR234aLL76YN998k/bt29OjRw/atWsXDgd+99138+STTzJ58uSYfmZPPvkkbdq0YdWqVRw9epSBAwcybNgwzj333JiWF5Ha0aRDoaIhaE/WHquKDh+VlJSQkpLCjTfeyNe//vVyN7Q5mZEjRwLQu3dvDh48SFpaGmlpaaSkpFBYWEjLli0rHBb7ZF577bXwl3XPnj0555xzwlAYMmQIbdq0ASAzM5MPPvjghFAYP348l1xyCb/61a+YO3duOPDe2rVrufvuuyksLOTgwYPk5MT+pfSXXnqJd999N9xLKioq4r333lMoiNSxJh0KdTE2eZnmzZuzcuVKFi9ezNy5c5kxYwZLliypcrkWLSK1NGvWLHxd9r6kpIQ5c+ZUOSz28U423lX0OpKSkigpOfHmHd26dSM9PZ1ly5bx/PPPhyOoTpo0iQULFnDhhRcya9YsXnnllROWjR6iO7pOd+c3v/lNtYJERGpfkx46e8Co7jQ/rfyPoLbHJi9z8OBBioqKuPLKK/n1r38d7klUNRR2VSobFvtknzto0CDmzJkDwObNm/nwww/JyKjecOETJkzgjjvuoHv37nTt2hWIDO/dqVMniouLw88/Xnp6OqtXrwYod+I7JyeHxx9/nOLi4rCuQ4cOVasmEam5Jh0K8RqbvOycQtlj6tSpHDhwgBEjRtCnTx8GDx7Mo48+CkQOxTz00EP07duXLVu2VHtdEydOJC8vj+zsbObMmRMOi92+fXsGDhzIBRdcwJQpU8otc+utt1JaWkrv3r0ZN24cs2bNKreHEItrrrmGdevWMX78+LDt5z//ORdffDFDhw6tcHhugDvvvJPHH3+cSy65hD179oTtN954I5mZmVx00UVccMEFfPe7361wL0VE4ktDZ0uDou0rUnMJGTpbREQaHoWCiIiEGmUoNORDYlI5bVeR+Gt0oZCSksLevXv1C6SRcXf27t1LSkpKoksRadQa3fcUunbtyvbt29G9FhqflJSU8PJXEYmPRhcKycnJ+hasiMgpanSHj0RE5NQpFEREJKRQEBGRkEJBRERCCgUREQkpFEREJNToLkkVkaZjTu56dizbzemlzqdJRpfBHZk4NjPRZTVo2lMQkQZpTu56CpbsomUpGEbLUihYsos5uesTXVqDplAQkQZpx7LdJGPl2pIxdizbnaCKGgeFgog0SKeXVjy+WWXtEhuFgog0SJ8mWbXaJTYKBRFpkLoM7kgx5fcKinG6DK7Z7XSbOoWCiDRIE8dm0uHyThxKAsc5lAQdLu+kq49qSJekikiDNXFsJigEapX2FEREJKRQEBGRkEJBRERCCgUREQkpFEREJBS3UDCzbma21Mw2mNk6M7s9aG9nZi+b2XvB8xlRy0wzs/fNbJOZ5cSrNhERqVg8L0ktAX7k7m+ZWRqw2sxeBiYBi919uplNBaYCPzGzTGA80AvoDCwys/PcvTSONYokzOYVu1m+cAsH9x2lVbsWDBjVnfMu1hevJLHitqfg7rvc/a3g9QFgA9AFGAXMDmabDYwOXo8C5rr7UXffCrwP9I9XfSKJtHnFbpbO2cjBfUcBOLjvKEvnbGTzCg3mJolVJ+cUzCwd6AusAM52910QCQ7grGC2LsA/oxbbHrQd/1k3m1memeUVFBTEtW6ReFm+cAslnx0r11by2TGWL9ySoIpEIuIeCmbWCnge+KG77z/ZrBW0nTDcobvPdPdsd8/u0KFDbZUpUqfK9hBibRepK3ENBTNLJhIIc9z9haD5IzPrFEzvBHwctG8HukUt3hXYGc/6RBKlVbsW1WoXqSvxvPrIgCeBDe7+SNSkF4HrgtfXAQuj2sebWQszOxfoAayMV30iiTRgVHesefmdY2tuDBjVPUEViUTE8+qjgcC1wBozyw/a7gKmA7lm9h3gQ+AaAHdfZ2a5wHoiVy7dpiuPpLFaf1opf0st5isHm9Hajf3mvJl6jM6nlXJeoouTJs3cG+5dirKzsz0vLy/RZYhU28DpS9hRePiE9i5tU3l96uUJqEiaEjNb7e7ZFU3TN5pFEmBnBYFwsnaRuqJQEEmAzm1Tq9UuUlcUCiIJMCUng9TkpHJtqclJTMnJSFBFIhG685pIAozuG/le5kN/38TOwsN0bpvKlJyMsF0kURQKIgkyum8XhYDUOzp8JCIiIYWCiIiEFAoiIhJSKIiISEihICIiIYWCiIiEFAoiIhJSKIiISEihICIiIYWCiIiEFAoiIhKKKRTM7IQbx1bUJiIiDVusewrLY2wTEZEG7KSjpJpZR6ALkGpmfYGyO423Bk6Pc20iIlLHqho6OweYBHQFHolqPwDcFaeaREQkQU4aCu4+G5htZmPc/fk6qklERBIk1pvs/NnMvgWkRy/j7vfHoygREUmMWENhIVAErAaOxq8cERFJpFhDoau7D49rJSIiknCxXpL6hpn1jmslIiKScLHuKXwVmGRmW4kcPjLA3b1P3CoTEZE6F2soXBHXKkREpF6I6fCRu38AdAMuD15/GuuyIiLScMQ69tG9wE+AaUFTMvBf8SpKREQSI9a/9q8GRgKHANx9J5AWr6JERCQxYg2Fz9zdAQcws5bxK0lERBIl1lDINbPfA23N7CZgEfDEyRYws6fM7GMzWxvVdp+Z7TCz/OBxZdS0aWb2vpltMrOcU+mMiIjUTExXH7n7w2Y2FNgPZAD3uPvLVSw2C5gB/PG49kfd/eHoBjPLBMYDvYDOwCIzO8/dS2OpT0REakesl6QShEBVQRA9/6tmlh7j7KOAue5+FNhqZu8D/dE9G0RE6lSsVx99w8zeM7MiM9tvZgfMbP8prvP7ZvZucHjpjKCtC/DPqHm2B20V1XKzmeWZWV5BQcEpliAiIhWJ9ZzCfwAj3b2Nu7d29zR3b30K63sc6A5kAbuAXwXtVsG8XtEHuPtMd8929+wOHTqcQgkiIlKZWEPhI3ffUNOVuftH7l7q7seInKjuH0zaTuTLcWW6Ajtruj4REameWM8p5JnZPGABUUNnu/sL1VmZmXVy913B26uBsiuTXgSeNrNHiJxo7gGsrM5ni4hIzcUaCq2JDG0xLKrNgUpDwcyeAS4FzjSz7cC9wKVmlhUsuw34LoC7rzOzXGA9UALcpiuPRETqnkW+k9YwZWdne15eXqLLEBFpUMxstbtnVzQtpj0FM0sBvkPkewQpZe3ufkOtVCgiIvVCrCea/wR0BHKAZUROBB+IV1EiIpIYsYbCl9z934FD7j4b+DqgO7GJiDQysYZCcfBcaGYXAG2A9LhUJCIiCRPr1Uczg28f/zuRy0dbAffErSoREUmIWAfE+0PwchnwxfiVIyIiiRTr1UdtgX8lcsgoXMbdfxCXqkREJCFiPXz0F+BNYA1wLH7liIhIIsUaCinu/m9xrURERBIu5u8pmNlNZtbJzNqVPeJamYiI1LlY9xQ+Ax4CfsrnQ1o7OuksItKoxBoK/0bkC2x74lmMiIgkVqyHj9YRGSVVREQasVj3FEqBfDNbSvn7KeiSVBGRRiTWUFgQPEREpBGrMhTMLAm41t3/pQ7qERGRBKrynEJwB7RPzaxNHdQjIiIJFOvhoyPAGjN7GThU1qhzCiIijUusofA/wUNERBqxWEdJnW1mqcAX3H1TnGsSEZEEiel7CmZ2FZAP/C14n2VmL8axLhERSYBYv7x2H9AfKARw93zg3LhUJCIiCRNrKJS4e9FxbV7hnCIi0mDFeqJ5rZl9C0gysx7AD4A34leWiIgkQqx7CpOBXkSGuHgaKAJ+GKeaREQkQU66p2BmKcAtwJeI3HVtgLuX1EVhIiJS96raU5gNZBMJhCuAh+NekYiIJExV5xQy3b03gJk9CayMf0kiIpIoVYVCcdkLdy8xsziXIw3F5hW7Wb5wCwf3HaVVuxYMGNWd8y7umOiyRKSGqgqFC81sf/DagNTgvQHu7q3jWp3US5tX7GbRnzbgJZGrkg/uO8qiP20AUDCINHAnPafg7knu3jp4pLl786jXJw0EM3vKzD42s7VRbe3M7GUzey94PiNq2jQze9/MNplZTs27JvGy5NnNYSCU8RJnybObE1SRiNSWWC9JPRWzgOHHtU0FFrt7D2Bx8B4zywTGE7nsdTjw2+A+DlIPlRys+AK0ytpFpOGIWyi4+6vAvuOaRxG5oongeXRU+1x3P+ruW4H3iQyrIfXQfjtWrXYRaTjiuadQkbPdfRdA8HxW0N4F+GfUfNuDNqmH1rZvRvFxo5wU46xtX9f/nESkttWX/8UVXdZU4dhKZnazmeWZWV5BQUGcy5KKjB2TwZK0UorsGI5TZMdYklbK2DEZiS5NRGoo1rGPastHZtbJ3XeZWSfg46B9O9Atar6uwM6KPsDdZwIzAbKzszUoXwKM7hvZiXvo75vYWXiYzm1TmZJzftguIg1XXYfCi8B1wPTgeWFU+9Nm9gjQGeiBvihXr43u20UhINIIxS0UzOwZ4FLgTDPbDtxLJAxyzew7wIfANQDuvs7McoH1QAlwm7uXxqs2ERGpWNxCwd0nVDJpSCXzPwg8GK96RESkavXlRLOIiNQDCgUREQkpFEREJKRQEBGRkEJBRERCCgUREQkpFEREJKRQEBGRkEJBRERCCgUREQkpFEREJKRQEBGRkEJBRERCCgUREQkpFEREJKRQEBGRkEJBRERCCgUREQkpFEREJKRQEBGRkEJBRERCCgUREQkpFEREJKRQEBGRkEJBRERCCgUREQkpFEREJKRQEBGRkEJBRERCCgUREQkpFEREJNQ8ESs1s23AAaAUKHH3bDNrB8wD0oFtwFh3/yQR9YmINFWJ3FO4zN2z3D07eD8VWOzuPYDFwXsREalD9enw0ShgdvB6NjA6caWIiDRNiQoFB14ys9VmdnPQdra77wIIns+qaEEzu9nM8swsr6CgoI7KFRFpGhJyTgEY6O47zews4GUz2xjrgu4+E5gJkJ2d7fEqUESkKUrInoK77wyePwbmA/2Bj8ysE0Dw/HEiahMRacrqPBTMrKWZpZW9BoYBa4EXgeuC2a4DFtZ1bSIiTV0iDh+dDcw3s7L1P+3ufzOzVUCumX0H+BC4JgG1iYg0aXUeCu7+D+DCCtr3AkPquh4REflcfbokVUREEkyhICIiIYWCiIiEFAoiIhJSKIiISEihICIiIYWCiIiEFAoiIhJSKIiISEihICIiIYWCiIiEFAoiIhJSKIiISEihICIiIYWCiIiEFAoiIhJSKIiISEihICIioUTcoznh5uSuZ8ey3Zxe6nyaZHQZ3JGJYzMTXZaISMI1uT2FObnrKViyi5alYBgtS6FgyS7m5K5PdGkiIgnX5EJhx7LdJGPl2pIxdizbnaCKRETqjyYXCqeXerXaRUSakiYXCp8mWbXaRUSakiYXCl0Gd6SY8nsFxThdBndMUEUiIvVHkwuFiWMz6XB5Jw4lgeMcSoIOl3fS1UciIjTRS1Injs0EhYCIyAma3J6CiIhUTqEgIiIhhYKIiIQUCiIiElIoiIhIyNwb7jd5zawA+KAGH3EmsKeWykmkxtIPUF/qo8bSD1Bfypzj7h0qmtCgQ6GmzCzP3bMTXUdNNZZ+gPpSHzWWfoD6EgsdPhIRkZBCQUREQk09FGYmuoBa0lj6AepLfdRY+gHqS5Wa9DkFEREpr6nvKYiISBSFgoiIhBplKJjZcDPbZGbvm9nUCqZPMbP84LHWzErNrF0sy9a1GvZlm5mtCabl1X315eqsqh9tzOy/zewdM1tnZtfHumxdq2Ff6s02Ceqpqi9nmNl8M3vXzFaa2QWxLlvXatiXerNdzOwpM/vYzNZWMt3M7D+Dfr5rZhdFTav5NnH3RvUAkoAtwBeB04B3gMyTzH8VsORUlq3PfQnebwPObAjbBLgL+GXwugOwL5i3wW2TyvpSn7ZJNfryEHBv8LonsPhU/m3W577Uw+0yCLgIWFvJ9CuBvwIGfAVYUZvbpDHuKfQH3nf3f7j7Z8BcYNRJ5p8APHOKy8ZbTfpSn8TSDwfSzMyAVkR+kZbEuGxdqklf6ptY+pIJLAZw941AupmdHeOydakmfalX3P1VIv9mKjMK+KNHvAm0NbNO1NI2aYyh0AX4Z9T77UHbCczsdGA48Hx1l60jNekLRH45vWRmq83s5rhVWbVY+jEDOB/YCawBbnf3YzEuW5dq0heoP9sEYuvLO8A3AMysP3AO0DXGZetSTfoC9Wu7VKWyvtbKNmmMd16zCtoqu+72KuB1dy9L5eosWxdq0heAge6+08zOAl42s43BXyF1LZZ+5AD5wOVAdyL1/m+My9alU+6Lu++n/mwTiK0v04H/a2b5RALubSJ7PQ1xu1TWF6hf26UqlfW1VrZJY9xT2A50i3rflchfbBUZT/nDLdVZti7UpC+4+87g+WNgPpHdy0SIpR/XAy8Eu8TvA1uJHPdtiNuksr7Up20CMfTF3fe7+/XungX8K5FzJFtjWbaO1aQv9W27VKWyvtbONkn0SZXafhDZ+/kHcC6fn2zpVcF8bYgct2tZ3WUbSF9aAmlRr98AhtfXfgCPA/cFr88GdhAZBbLBbZOT9KXebJNq9KUtn58kv4nIsewG+X/lJH2pV9slqCOdyk80f53yJ5pX1uY2SVin4/wDvRLYTORM/E+DtluAW6LmmQTMjWXZhtgXIlcgvBM81iW6L1X1A+gMvERkt34t8O2Guk0q60t92yYx9mUA8B6wEXgBOKMBb5cK+1LftguRPf5dQDGRv/6/c1w/DHgs6OcaILs2t4mGuRARkVBjPKcgIiKnSKEgIiIhhYKIiIQUCiIiElIoiIhISKEgDVYwImx+1KOikTEvNbM/J6K+45nZK2YW843WzWySmc2oZNobwXN62WiaZpZtZv8ZvL7UzC6pjbqlaWmMw1xI03HYI99OrTfMrLm7x33wO3c/4Re+u+cBZcM+XwocJPJFLJGYaU9BGp1gTPmNZvYawQBoQXsHM3vZzN4ys9+b2QdmdmYw7dvBGPv5wbSkCj53m5n9MphvpZl9KWifZWaPmNlS4JdmlmVmbwZj3c83szOiPubbZvaGRe590T9Yvn/Q9nbwnBE1fzcz+1swRv69UbUcrKC+S83sz2aWTuTLTncE/fmamW01s+RgvtZBX5Jr8GOWRkqhIA1Z6nGHj8aZWQrwBJEBAr8GdIya/14i95u4iMj4Nl8AMLPzgXFEBkXLAkqBiZWsc7+79ycyEuqvo9rPA/7F3X8E/BH4ibv3IfKN03uj5msZ/JV/K/BU0LYRGOTufYF7gP8TNX//oJYs4JpYDj+5+zbgd8Cj7p7l7v8LvEJkeASIjJP1vLsXV/VZ0vTo8JE0ZCccPjKzLGCru78XvP8voGwo5K8CVwO4+9/M7JOgfQjwZWBV5BYIpAIfV7LOZ6KeH41qf9bdS82sDdDW3ZcF7bOBZ49f3t1fDf5ibwukAbPNrAeRUS2j/4J/2d33Bn15IejDqdwZ7A/Aj4EFRAbsu+kUPkOaAIWCNEaVjd1S0dDCZe2z3X1aNT87+vWhWArjxNoc+Dmw1N2vDg79vFLF/NXm7q8HJ6UHA0nuXuGtHkV0+Egam43AuWbWPXg/IWraa8BYADMbBpQd618MfDMYSx8za2dm51Ty+eOinpcfP9Hdi4BPzOxrQdO1wLKoWcYF6/gqUBTM34bISKoQGdww2tCgnlRgNPB6JXUd7wCRPZBofySyp/L/YvwMaYIUCtKQHX9OYbq7HyFyuOh/ghPNH0TN/zNgmJm9BVxBZCTKA+6+HribyJ233gVeBjpVss4WZrYCuB24o5J5rgMeCj4rC7g/atonweWkvyMy+iXAfwC/MLPXidxnN9prwJ+I3LTn+eAKo1j8N3B12YnmoG0OkSCsj7dslXpCo6RKk2FmLYBSdy8xswHA49W5pNXMthEZpnhPnEqMKzP7JjDK3a9NdC1Sf+mcgjQlXwByzawZ8BlN6GSrmf2GyN7RlYmuReo37SmIiEhI5xRERCSkUBARkZBCQUREQgoFEREJKRRERCT0/wH0geA9oT5XCwAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "Ns = 500\n",
    "m = 6\n",
    "probs = np.append(np.arange(0.7, 1, 0.08), 1)\n",
    "print(probs)\n",
    "estimation_list = []; exact_list = []\n",
    "for i in probs:\n",
    "    print(\"i in probs:\", i)\n",
    "    estimation_values = []; exact_values = []\n",
    "    for _ in range(4):\n",
    "        random_graph = random_graph_generation(m, i)\n",
    "        estimation_values.append( perm_estimation(random_graph, Ns) )\n",
    "        exact_values.append( qc.permanent_cx(nx.to_numpy_array(random_graph)) )\n",
    "\n",
    "    estimation_list.append(np.mean(estimation_values))\n",
    "    exact_list.append(np.mean(exact_values))\n",
    "\n",
    "print(\"exact_list\", exact_list)\n",
    "print(\"estimation_list\", estimation_list)\n",
    "plt.plot(probs, exact_list, \"ob\", label=\"Exact Value\", color='C0')\n",
    "plt.plot(probs, estimation_list, \"ob\",\n",
    "         label=\"Estimation Value\", color='tab:purple')\n",
    "\n",
    "plt.legend()\n",
    "plt.xlabel('Edge probability')\n",
    "plt.ylabel('Permanent')\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "765bd6793ea5362367183dba09c599ce34c84eeeb7f3172847d560f9d0858670"
  },
  "kernelspec": {
   "display_name": "Python 3.9.12 ('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.9.12"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}