QBayesianRejectionSampling.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "source": [
    "# Import libraries"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from qiskit import QuantumRegister\n",
    "from qiskit import QuantumCircuit, Aer, transpile\n",
    "from qiskit.circuit.library import GroverOperator"
   ]
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Implement quantum circuit"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "outputs": [
    {
     "data": {
      "text/plain": "<Figure size 2528.36x284.278 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAB64AAADuCAYAAACJbykBAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAABhFklEQVR4nO3dd3gU1eL/8c8mhAAJLYCEovReDYjIRUFsiMAVrxQVkSJKu1YEFUEviqLY0CsCFwSxIIpeLiIoooA0pTcBpUoz1EBIIW3P7w9/7NdIEnY3uzNb3q/n2UfJlPOZnbPnzM7ZmXEYY4wAAAAAAAAAAAAAALBJhN0BAAAAAAAAAAAAAADhjYFrAAAAAAAAAAAAAICtGLgGAAAAAAAAAAAAANiKgWsAAAAAAAAAAAAAgK0YuAYAAAAAAAAAAAAA2IqBawAAAAAAAAAAAACArRi4BgAAAAAAAAAAAADYioFrAAAAAAAAAAAAAICtGLgGAAAAAAAAAAAAANiKgWsAAAAAAAAAAAAAgK0YuAYAAAAAAAAAAAAA2IqBawAAAAAAAAAAAACArRi4BgAAAAAAAAAAAADYioFrAAAAAAAAAAAAAICtGLgGAAAAAAAAAAAAANiKgWsAAAAAAAAAAAAAgK0YuAYAAAAAAAAAAAAA2IqBawAAAAAAAAAAAACArRi4BgAAAAAAAAAAAADYioFrAAAAAAAAAAAAAICtGLgGAAAAAAAAAAAAANiKgWsAAAAAAAAAAAAAgK0YuAYAAAAAAAAAAAAA2IqBawAAAAAAAAAAAACArRi4BgAAAAAAAAAAAADYioFrAAAAAAAAAAAAAICtGLgGAAAAAAAAAAAAANiKgWsAAAAAAAAAAAAAgK0YuAYAAAAAAAAAAAAA2IqBawAAAAAAAAAAAACArRi4BgAAAAAAAAAAAADYioFrAAAAAAAAAAAAAICtGLgGAAAAAAAAAAAAANiKgWsAAAAAAAAAAAAAgK0YuAYAAAAAAAAAAAAA2IqBawAAAAAAAAAAAACArRi4BgAAAAAAAAAAAADYqojdAcKBMUZpaWl2x/BIiRIl5HA4fLKucN9+UAdAHQCAcEc/AAAId8HWF9IPAoBvBVs/INEXALAHA9cWSEtLU2xsrN0xPJKSkqKYmBifrCvctx/UAVAHACDc0Q8AAMJdsPWF9IMA4FvB1g9I9AUA7MGtwgEAAAAAAAAAAAAAtuKKa4sdO3YsYH+llJqaqooVK/q1jHDfflAHQB0AgHBHPwAACHeB2hfSDwKANQK1H5DoCwDYj4Fri8XExARsp2SFcN9+UAdAHQCAcEc/AAAId/SFABDe6AcAIH/cKhwAAAAAAAAAAAAAYCsGrgEAAAAAAAAAAAAAtmLgGgAAAAAAAAAAAABgKwauAQAAAAAAAAAAAAC2YuA6QC1btkwOhyPXKzY2VgkJCXrjjTeUnZ1td0QAAAAAAAAAAAAA8IkidgdAwe666y516tRJxhglJiZq1qxZeuyxx7Rz505NnTrV7ngAAAAAAAAAAAAAUGgMXAe4hIQE9e7d2/XvIUOGqH79+po2bZrGjRunChUq2JgOAAAAAAAAAAAAAAqPW4UHmZiYGLVu3VrGGO3du9fuOAAAAAAAAAAAAABQaCE/cJ2SkqJKlSrJ4XCodOnSMsbYHanQLgxYx8XF2ZwEAAAAAAAAAAAAAAov5G8VPn78eCUmJkqSkpOTdeDAAdWoUcPmVO5LS0vTyZMnXc+4njx5sjZt2qRWrVqpbt26dscDAAAAAAAAAAAAgEIL6YHrw4cP6/XXX1dERIRiY2OVnJysrVu3BtXA9bPPPqtnn30219/uuOMOvfPOOzYlgt2Sk5P1wQcfaOXKlUpNTVXJkiV14403qmfPnipRooTd8WCBgwcPatq0adq5c6cyMzNVvnx5de/eXTfffLMiIkL+RhpA2MvJydHChQv1xRdf6NSpUypWrJiaNGmiAQMGqHLlynbHgwVSUlL00Ucfafny5UpJSVFsbKzat2+vu+++W7GxsXbHAwC/S0pK0syZM/XTTz8pLS1NpUuX1i233KLu3bsrOjra7niwwN69ezV9+nT9+uuvys7O1mWXXaZevXrp+uuvl8PhsDseAD/LysrS/PnzNX/+fCUlJal48eJq3ry5+vfvr4oVK9odDxY4e/asPvjgA61atSrX+dFevXqpePHidscD/O78+fP67LPP9M033+js2bOKiYnRNddco/vuu09lypSxOx4scOLECc2YMUMbNmxQenq6ypQpo86dO+v2229X0aJF7Y5XOCaE9e7d20gy9957r7ntttuMJDN27FjLc6SkpBhJRpJJSUlxa5mlS5caSeaBBx4w3377rVm4cKF5+eWXTVxcnGnVqpVJSkpyzduzZ0/TvXv3XMufOnXKxMfHmw8//NCvOa1c71NPPWUkmenTp180zel0mnbt2pmiRYuabdu22ZrTXzIyMsyjjz5qYmJiXDn//Cpbtqx57rnnTE5Ojt1RL0Id8I2jR4+a22+/3URERORZB2rVqmXmzJljd8w8UQcA35g5c6a54oor8mwDIiMjTY8ePcyJEyfsjgk/ycrKMk8++aQpVapUnnWgVKlS5sknnzRZWVl2R70I/QAAX0hNTTWDBw82xYsXz7MdLF++vHn55ZeN0+m0Oyr85MCBA+bWW281DocjzzpQv359M3/+fLtj5ikY+kL6QQSDd99911SuXDnPNiAqKsr06dPHnDlzxu6Y8JPz58+bhx9+uMDzo//61784P8p3opCVk5NjXnzxRVOuXLk8PwPFixc3Q4cONenp6XZHhZ8kJyeb/v37m6JFi+ZZB+Lj481bb70V1N+JQnbgev369cbhcJhixYqZ3377zYwYMcJIMnfeeaflWQozcD1hwoRcf1+1apVxOBymZ8+err+dOnXKVK5c2Xz88ceuv/Xq1cvjbQ30geuMjAzTuHFjU7p0aXPo0KFc015//XUjybz00ku25/SH8+fPmxtvvDHPhuivr549e5rs7Gy7I+dCHSi8AwcO5DtY9dfXxIkT7Y57EeoAUHgvvPCCW21AnTp1zNGjR+2OCx/LzMw0nTt3dqsOdOnSxWRmZtodORf6AQCFlZKSYtq0aeNWOzhgwICgPlGDvO3atcvEx8dfcv87HA4zbdo0u+NeJBj6QvpBBLonn3zSrX6gSZMm5uTJk3bHhY+lp6ebDh06uFUH7rrrroAbvA6GfsCXOeF7OTk5pm/fvm59Bq699lqTmppqd2T4WFJSkrnyyivdqgOPPPJI0H4nCtmB6+uuu85IMiNHjjTGGPP+++8bSaZu3bqWZ/HlwLUxxvTp08dIMqtWrXL9bdGiRSYuLs4cOXLEfPbZZyY+Pt7jA7RAH7g2xpgNGzaYIkWKmJtvvtn1t127dpnixYubq6++ulADtoHcKQ8YMMCtxujCa9SoUXZHzoU6UDgZGRmmYcOGHtWBhQsX2h07F+oAUDiffvqpR21Ay5YtA+5LOgrnoYce8qgOPPzww3ZHzoV+AEBhde/e3aN2cPz48XZHhg+lpqaaGjVquL3/IyMjzfLly+2OnUsw9IX0gwhk06dP96gfaN++fdCesEfe3B2wu/AaPXq03ZFzCYZ+wNc54Vvjxo3z6DNw11132R0ZPtaxY0eP6sCkSZPsjuyVkBy4/uKLL4wkU65cOdetYdavX28kmYiICJOWlmZpHl8PXO/evdtERkaaG264IdffBw0aZK6//npToUIFs2DBAkty2rHe0aNHG0lmypQpJjs727Rq1coUK1bM7Nq1K6By+srBgwfzvTV0fq/Y2Fhz9uxZu6O7UAcKZ86cOR7tf0mmbdu2dsfOhToAeM/pdJqmTZt63A4sWrTI7ujwkWPHjuV7C6j8XtHR0QF123j6AQCFsXPnTo/7wfLly3OLxBDyn//8x+M60KlTJ7tj5xIMfSH9IAJVTk6OqVWrlsftwMqVK+2ODh85cOCAx+dHS5YsaZKTk+2O7hIM/YA/csI3UlNTTdmyZT1uB3fv3m13dPjIunXrPN7/l19+eUA+Tu5SIhRisrKyNHLkSEnS6NGjVbp0aUlSw4YNFRERIafTqe3bt9sZsdBq166tXr166bvvvtOKFStcf3/11Ve1Z88e3XrrrbrttttsTOhfo0ePVrNmzTR8+HD985//1Nq1azVu3DjVq1fP7mh+MXXqVDmdTo+WSUlJ0YcffuinRPYLtzowadIkj5dZuXKltm7d6oc0gSHc6gDC25o1a7z6PHvTdiAwvffee8rMzPRomYyMDM2YMcNPiexHPwCEl8mTJ3u8zMmTJzV37lw/pIHVjDF65513PF5u0aJF2rdvnx8SBQb6QoSTxYsXa+/evR4vx3ei0OHN+dFz587po48+8lMi+9EPhJc5c+YoKSnJ4+W8OY5GYPKmTzt06JC++uorP6Txr5AbuH7nnXe0e/du1axZU4MHD3b9vXjx4qpRo4YkhcRgzqhRoxQREaExY8a4/hYTE6OaNWuqSZMmNibzv6ioKL3//vs6f/683n33XbVt21aPPPKI3bH85uuvv/ZquUWLFvk4SeAIpzqQnp6u5cuXe7Wst3UnGIRTHQC8bc8XLVokY4yP08AOHAtcjH4ACC+F6QsR/E6cOKHNmzd7vJwxRosXL/Z9oABBX4hwwvEwqAMXox8IL3wGEE51oIjdAXzp9OnTGjt2rCTppZdeUtGiRXNNb9Sokfbu3astW7Z4XUbLli2VmJjo0TKe/hpMktq3b1/gyeYGDRooJyfH4/W6q06dOoqI8M3vGrzZ/kspXbq0oqOjlZWVpU6dOvks6wW+3P7C8rS+XbB48WJVrVrVx2m8Qx3wXmE+588//7zeeustH6bxHnUA8N6ZM2e8Wi47O1tVqlShHoeA48ePe7XcypUrORYoBPoBIHD8/vvvXi33+eefe/0jUASO7Oxsr5cdOXKkXnjhBR+m8V6w9YX0gwgk3lxleGG5KlWqyOFw+DgRrObt+dGvv/6a70SFQF8QOE6ePOnVcrt27QqYzwAKx9t2cNasWVqwYIGP01xafHy81q9f793C9t6p3LceeeQRI8m0atUqz+lPPfWUkWTatWvndRlVqlTx+D7yf375+7kQ7dq1y/O52O748/Mr/PXyxfY7nU7Tvn17U7RoUdOgQQNTokQJs2fPnkKv14rt50Ud4EUd4MWLF69wf9EP8OLFixevcH8Fal9IP8iLFy9e1rwCtR+gL+DFi5evXlWqVPG6HQqZK6737NnjeubR2rVrC/wl3bZt27wuJz4+3uNlnE6n178Qt0ulSpV8esW1L7f/7bff1rJlyzRu3Dj9/e9/V0JCgvr3769ly5b57BeUvtz+wkpKSlJaWprHy5UsWVKlSpXyQyLPUQe8Z4zR8ePHvbrKoFy5cipWrJgfUnmOOgB4Lz09XadPn/Z4uaioKF122WV+SASrnTlzRqmpqR4vFxsbq9KlS/shkefoBwAUxunTp5Wenu7xcqVKlVLJkiX9kAhWMsbo2LFjXt2NqkKFChfdjc8uwdYX0g8ikKSmpnp1J6ro6GiVL1/e94FgOc6PXozvROHl3LlzSk5O9ni54sWLKy4uzg+JYLWTJ08qIyPD4+XKlCmjmJgYPyQqmDdjqS6F/glOgOjWrZuRZEqWLGkqVqyY56tChQqu0f5Dhw5Zlu3Pv1Ly9xXXheGvnL5c76+//mpKlChhrrrqKpOdnW2MMebFF180kszEiRMDJqcvbdq0yeNfs0RGRlpaxy+FOlA4EydO9LgO1K1b1+Tk5Ngd3YU6AHgvKyvLqzu+zJgxw+7o8JFdu3Z59evWX3/91e7oLvQDAApj+fLlHreB0dHR5uTJk3ZHh4+88MILHteBhIQE43Q67Y7uEgx9If0gAlVaWpqJi4vzuB34/PPP7Y4OH1m/fr1X50ePHDlid3SXYOgHfJ0TvnPs2DFTtGhRjz8Hq1evtjs6fGTBggUe7/9SpUqZc+fO2R3dYyExcH3hS2yRIkXMjh078p0vKyvL9eH+6quvLMsXLI19oA9c5+TkmDZt2pjo6Ohc+zk7O9u0bNmy0LdDCeT91K5dO48apF69etkdORfqQOEkJSWZMmXKeFQH3nnnHbtj50IdAArnpZde8qgNiI+PN2lpaXbHhg916tTJozrQuXNnuyPnQj8AoDCcTqdJSEjwqB0cOHCg3bHhQ4mJiSYmJsajOjBr1iy7Y+cSDH0h/SAC2dNPP+1RG1CjRg2TlZVld2z4UNu2bT2qA3fffbfdkXMJhn7Alznhe/369fPoM9CqVauA+hEfCic7O9vUq1fPozrw+OOP2x3bK0F/nwdjjB5//HFJ0qBBg9SgQYN85y1SpIjq1q0rSdq6dasl+eA7r732mlavXq2xY8fm2s+RkZGaOXOmsrOz1b9/fxljbEzpH7Nnz1a1atXcmrdJkyaaPHmynxPZI1zrQJkyZTRv3jxFR0e7Nf+9996rwYMH+zmVPcK1DgBPPPGEunXr5ta8MTExmj9/vooXL+7nVLDSzJkzVadOHbfmrVu3rmbMmOHnRPagHwDCk8Ph0Oeff+727eZatWqlN954w8+pYKWKFSvq008/VZEi7j3xbvDgwerdu7efU9mDvhDh6rnnntPNN9/s1rxlypTR/Pnz3W4zEBw++eQTXXHFFW7N26xZM7377rt+TmQP+oHw9dZbb6lly5ZuzVu5cmV99tlnPrttPOwXGRmpefPmqVy5cm7N3759e40bN87PqfzExkFzn5g1a5aRZMqWLevWbcB69uxpJJm77rrLgnR/CJZfKQXyFdc7duww0dHRpnXr1q7bn/xVYW+HEuj76fDhw+aaa64p8Bc0HTt2NKdPn7Y76kWoA76xcuXKAm8XHBkZaR577LF83x87UQeAwsvMzDSDBg0yERER+bYD1atXN+vXr7c7Kvzk2LFjl7wLS/v27c2xY8fsjnoR+gEAvrBv3z5z5ZVXFtgO3n777SY5OdnuqPCTb7/91lSsWDHf/R8VFWWeeeaZgHps0gXB0BfSDyLQpaenm/vuu884HI5824F69eqZ7du32x0VfnLo0CFz9dVXF3gscOutt5qkpCS7o14kGPoBX+WE/5w9e9Z07dq1wM9AixYtzIEDB+yOCj/ZtWuXadiwYYF14K677grqOzE6jAnen96kp6erbt26Onz4sN544w098sgjl1zm+eef15gxY9SoUSNt377d/yElpaamKjY2VpKUkpJiy4PQ3eGvnOG+/b5kjNGPP/6oSZMmafbs2crJyVFkZKTuv/9+DR48WM2aNbM7Yp6C4b2VgiNnVlaWvvzyS02ZMkXfffedcnJyVKRIET3zzDO6//77VaVKFbsj5ikY3lspeHIivP3222+aOnWqvvjiC/36669yOp2Kjo7W3LlzdeuttyoyMtLuiPCzdevW6d1339WsWbNcxwL33XefBg8e7Pavr60WLO1rsOQEwpkxRitWrNCkSZM0d+5cVzs4ePBgDR48WA0bNrQ7IvwsMzNTX3zxhaZOnaoffvjB9Z3oX//6lwYMGKCKFSvaHTFPwdDHBENGQJL27t2rKVOmaP78+dq9e7ecTqeKFSumefPm6aabblJERNDfZBQFMMZozZo1mjRpkj755BPOj/pYsOQMd9u3b9fkyZP19ddfa//+/XI6nSpevLi++eYbtW3bliutQ5zT6dT333+vSZMmaf78+a528OGHH9agQYPcvmNfoArqgetgESyNPQPXwZHzgqpVq+rIkSOqUqWKDh8+bHecAgXLexssOS+gDvhesOQELgimdgC+F0z7P1ja12DJCeAPwdQOwj+CqQ4EQx8TDBmBvwqmdgC+F0z7P1ja2GDJif8TTJ8D+F4o7n9+fgYAAAAAAAAAAAAAsBUD10Fs7ty5Gjx4sCTp5ptvVtOmTdW8eXNde+212rRpk83pAAAAAAAAAAAAAMA9RewOAO/997//VZ8+fSRJn376qcqUKeP6e9++fbVlyxYb0wEAAAAAAAAAAACAe7jiOoCdOXNGVatWVbly5dS8eXM1aNBA0dHRuv/++5WVlaVVq1apQ4cOkuQatJaks2fPyuFw2JQaAAAAAAAAAAAAADzDFdcBrEyZMrr77rtVsmRJjR49Wt98841efPFFTZs2TYsXL1abNm0UFRXlmr9Pnz5aunSpJGnhwoV2xQYAAAAAAAAAAAAAj3DFdYDbvHmzrrzySknShg0bXP8/b948devWLde8s2bN0qFDh/TCCy9o5MiRlmcFAAAAAAAAAAAAAG8wcB3g8hq4Nsbom2++0a233prnMvfdd5+WLl2qU6dOWRkVAAAAAAAAAAAAALzCwHUAO3LkiBwOh6pUqSJJ2rp1q5o0aaK1a9eqQYMGio2NlfTHs7CPHj3qWm7evHkqV66c4uLibMkNAAAAAAAAAAAAAJ7gGdcBbNOmTa6rraU/nnk9adIklS9fXrfffrvr72fPnlX37t2Vnp6uiIgIVahQQQsWLJDD4bAhNQAAAAAAAAAAAAB4hoHrANa5c2d17tzZ9e9169ZJkho1aqSlS5e6/l6tWjWtXbvW8nwAAAAAAAAAAAAA4AsMXAehn3/+2e4IAAAAAAAAAAAAAOAzDFxbLDU11e4I+bIiW7hvPwL7fQ7kbKEkkN/nQM4GAKEikNvaQM4GAAgdgdrfBGouAAg1gdzeBnI2AOGBgWuLVaxY0e4Itgr37Qd1ANQBAAh39AMAgHBHXwgA4Y1+AADyF2F3AAAAAAAAAAAAAABAeOOKawuUKFFCKSkpdsfwSIkSJXy6rnDeflAHQB0AgHBHPwAACHfB1hfSDwKAbwVbPyDRFwCwBwPXFnA4HIqJibE7hm3CfftBHQB1AADCHf0AACDc0RcCQHijHwAA93CrcAAAAAAAAAAAAACArRi4BgAAAAAAAAAAAADYioFrAAAAAAAAAAAAAICtGLgGAAAAAAAAAAAAANiKgWsAAAAAAAAAAAAAgK0YuAYAAAAAAAAAAAAA2IqBawAAAAAAAAAAAACArRi4BgAAAAAAAAAAAADYioFrAAAAAAAAAAAAAICtGLgGAAAAAAAAAAAAANiKgWsAAAAAAAAAAAAAgK0YuAYAAAAAAAAAAAAA2IqBawAAAAAAAAAAAACArRi4BgAAAAAAAAAAAADYioFrAAAAAAAAAAAAAICtGLgGAAAAAAAAAAAAANiKgWsAAAAAAAAAAAAAgK2K2B0ACAfGGKWlpfl0nU6n0/Xf1NRUn667RIkScjgcPl1nOPPH/peoA8HEX3XAn6gDAOA79AO8BwAQ7ugHAAAA4A4GrgELpKWlKTY21i/r/v33332+7pSUFMXExPh0neHMn/tfog4EA3/XAX+gDgCA79AP8B4AQLijHwAAAIA7uFU4AAAAAAAAAAAAAMBWXHENWOzYsWMB+Yvd1NRUVaxY0e4YIS9Q979EHbAKdQAAwhv9AO8BAIQ7+gEAAADkh4FrwGIxMTEB+wUN/sf+B3UAAMIb/QDvAQCEO/oBAAAA5IdbhQMAAAAAAAAAAAAAbMXANQAAAAAAAAAAAADAVgxcAwAAAAAAAAAAAABsxcA1AAAAAAAAAAAAAMBWDFwDAAAAAAAAAAAAAGzFwDUQoJYtWyaHw5HrFRsbq4SEBL3xxhvKzs62OyL8jDoAAAAAAAAAAADCRRG7AwAo2F133aVOnTrJGKPExETNmjVLjz32mHbu3KmpU6faHQ8WoA4AAAAAAAAAAIBQx8A1EOASEhLUu3dv17+HDBmi+vXra9q0aRo3bpwqVKhgYzpYgToAAAAAAAAAAABCXcjdKnzJkiUX3VrX4XAoIiJCZcuWVdu2bfXee+/JGGN3VMArMTExat26tYwx2rt3r91xYAPqAAAAAAAAAAAACDUhd8X1pk2bJEkVKlRQ3bp1XX9PTk7W3r17tWrVKq1atUp79+7VuHHj7IoJFMqFwcq4uDibk8Au1AEAAAAAAAAAABBKQu6K682bN0uS+vfvr5UrV7peW7du1dGjR9WpUydJ0htvvKGMjAwbkwLuSUtL08mTJ3XixAlt27ZNQ4cO1aZNm9SqVatcP85A6KIOAPirrKwsrVixQunp6ZKk7OxsmxPBamfOnNH58+clSefPn9eZM2fsDQQAFjt+/LirHczIyFBqaqrNiWC13bt3u+pAZmamcnJybE4EwEoZGRlavny56zsRbUD4OXXqVK7vROfOnbM5EWCt1NRUff/997SDYez333/P9Z3owv8Hu5AbuL5wxXWzZs0umla6dGk99thjkqT09HSdOHHC0myAN5599llVqFBBl112mZo2bapJkybpjjvu0P/+9z+7o8Ei1AFc8PTTT8vhcOi99967aJoxRu3bt1d0dLS2b99uQzpYITExUWPGjNEVV1yh6667TqdPn5YkHTt2TJ07d9bXX39tc0L427Zt2zRw4EBVrlxZp06dkvTHCZsqVarogQce4PMfBsK9Lwj37Yf0008/6Z577lHVqlVd7eDJkydVpUoVPfLIIzxKJ8QZY/TZZ5/p+uuvV926dV114MSJE6pdu7ZefvllfswV4ugHcPjwYT311FOqWrWq2rdv7/pOlJiYqDvuuENLly61OSH8bcOGDbrvvvtUpUqVi74TDR06VL/88ovNCQH/2rNnjx566CFVrlxZN9xwQ652sHfv3lq7dq3NCeFvK1asUPfu3XX55Zfn+k5UtWpVjRgxQgcPHrQ5YSGZEJKWlmYiIyONJPPzzz/nOc/8+fONJBMdHW2ys7MtTohwlZKSYiQZSSYlJcWtZZYuXWokmQceeMB8++23ZuHChebll182cXFxplWrViYpKck1b8+ePU337t1zLX/q1CkTHx9vPvzwQ79lhHu8fW+pA6HDV+9tRkaGady4sSldurQ5dOhQrmmvv/66kWReeukl23PCPzZt2mQqVark2kf5vR577DHjdDrtjgs/+Pjjj01UVFSB+79o0aJm9uzZdkfFX/iyffVnX+DPfoC+EL4wceJE43A4CmwHY2NjzTfffGN3VPhBZmamueeeey55LFSnTh2zd+9eu+PiL+gH4As//vijKV++/CXbgWeffZbvRCHqP//5j+v8f36v4sWLm/nz59sdFfCLRYsWmZiYmAI/AxEREeadd96xOyr8ZNy4cZfsB8uWLWtWrFhhd1SvhdTA9Y8//mgkmWLFiuU7KH3nnXcaSeaee+6xOB3CWWEGridMmJDr76tWrTIOh8P07NnT9bdTp06ZypUrm48//tj1t169epk777zTrxnhnsIOXFMHgp8v39sNGzaYIkWKmJtvvtn1t127dpnixYubq6++ulA/yqIOBK69e/eaChUqXPLA9MJr1KhRdkeGj3355ZcmIiLCrf0fERFhFixYYHdk/Imv21d/9QXBMHBtDH1huJo+fbrb/WCxYsXMmjVr7I4MH3I6naZfv35u14FatWqZ48eP2x0bf0I/gML6+eefTZkyZdxuB15++WW7I8PHZs+e7fb+j4qKMkuXLrU7MuBTq1atMtHR0W5/Dt5//327I8PHJk6c6Pb+L1mypNmyZYvdkb0SUrcKv/B868aNGysyMtL197Nnz2rt2rXq0aOH5s6dq/r16+uVV16xKSVQOG3atNG9996rOXPmaPXq1ZKkuLg4TZ8+XcOGDdPRo0c1d+5cLVu2TJMnT7Y5LfyBOhDeEhIS9NRTT2nx4sWaOnWqcnJy1KdPHxlj9P777+fq/xA6nn76aY8ecfLiiy9yq9QQkp2drSFDhsjpdLo1v9Pp1JAhQ3i+VQgL974g3Lc/HJ07d06PPPKI2/OfP39eDz30kP8CwXKrV6/WjBkz3J5/7969Gj9+vB8TwU70A+HpiSee8OhRAKNGjVJiYqL/AsFS58+f17Bhw9yePysrS0OHDpUxxo+pAOsYYzR06FBlZGS4vcxDDz2k1NRUP6aClU6ePKkRI0a4Pf+5c+dcj04ONiE1cH3h+dbr16+Xw+FwvcqUKaOrr75aS5Ys0Ysvvqgff/xRlStXtjkt4L3Ro0crMjJSY8aMcf2tY8eO6tGjh3r37q0hQ4Zo2rRpKleunI0p4U/UgfA2evRoNWvWTMOHD9c///lPrV27VuPGjVO9evXsjgY/SExM1Oeff+7RMsYYTZkyxU+JYLWvvvpKhw4d8miZgwcPauHChX5KhEAQ7n1BuG9/uPnoo4907tw5j5ZZt26d1q9f76dEsNqkSZM8XmbGjBlKT0/3QxoEAvqB8LJv3z4tWrTIo2Wys7M1bdo0PyWC1T777DPXc1zdtWPHDv3www9+SgRY66effnJduOmus2fP6pNPPvFPIFhuxowZHv1wQZK+++477dq1y0+J/CekBq4vfHDr1aunv/3tb65Xw4YNVaxYMSUlJWnWrFk6cuSIvUGBQqpdu7Z69eql7777TitWrHD9/dVXX9WePXt066236rbbbrMxIfyNOhDeoqKi9P777+v8+fN699131bZtW4+uQkJwmT17trKzsz1e7v333/dDGthh5syZXi1HHQht4d4XhPv2hxtv2zPawdCQkpKizz77zOPlkpKS9OWXX/ohEQIB/UB4+eCDD7y6ctbb42gEHm/7dOoAQoW3dZnPQOjwdl/OmjXLt0EsUMTuAL7idDq1bds2SX/88uCaa67JNf306dPq06ePvvrqK/3jH//Qzz//rIgIz8ftW7ZsyW1m4DF3b+3piVGjRmn27NkaM2aMli5dKkmKiYlRzZo11aRJk0Ktu06dOl59PpA3f+x/iToQTPxRB0qXLq3o6GhlZWWpU6dOPt9f1IHA4cnt8P7s+PHjqlKlihwOh28DwXLHjx/3arn58+eratWqPk4Db/jrWMCffYGv+wH6QhTG77//7tVy06ZN03//+18fp4HVsrOzlZWV5dWyAwcODNpbJIYa+gEURlJSklfL7d27l+9EIeLYsWNeLTd79mx9++23Pk4DWO/kyZNeLbdmzRrOC4SIo0ePerXcm2++acvgdXx8vNd3wHKYEHnQw86dO9WwYUM5HA4lJycrNjb2onl++eUX1a9fX5K0detWrwZ2qlatyhXbKJSUlBTFxMT4bf3t27dX586dNXz4cI+WS01NzfNzA9/y9/6XqAOBzhd1wBijDh06aPXq1apVq5Z+++03bd26VbVq1SrUeqkDAOB/vjoW8EdfYFU/QF8IAOGNfgAAACC0ValSRYcPH/Zq2ZC54vrC861r1KiR7wFm9erVXf9/7Ngxrwau4+PjvcqH8OZ0Or2+SsAOlSpV4lfFPhRs+1+iDviar+vA22+/rWXLlmncuHH6+9//roSEBPXv31/Lli3z2S/JqQOBIyUlRWfPnvV4uYiICFWqVMkPiWC1U6dO6fz58x4vV6xYMZUrV84PieApfxwL+Lsv8HU/QF+Iwjhx4oQyMzM9Xi4mJkZlypTxfSBYqjDtR9myZVWiRAkfJ4I36AdQGMnJyTp37pzHy0VGRnIuN0ScPHnS42e7SlKJEiVUtmxZPyQCrHXmzBmlpqZ6vFzRokVVoUIFPySC1Y4dO+bVowRjY2NVunRpPyQqWKH6XxMinnjiCSPJdOvWLd959u7dayQZSWbr1q0WpkO4S0lJcdW9lJQUv5bVrl07M2HCBI+XszJjuLH6vaUOBB5fvre//vqrKVGihLnqqqtMdna2McaYF1980UgyEydODJic8J3ExEQTFRXl2jfuvkaOHGl3dPjI/PnzPd7/ksyCBQvsjo7/z9ftq7/6An/2A/SFKIwpU6Z41Q5u2LDB7ujwkd69e3u8/+Pi4kxaWprd0fH/0Q+gMPbv328cDofH7cC4cePsjg4f+fDDD706Fvjhhx/sjg74xE8//eTVZ+C9996zOzp85NVXX/WqDvzyyy92R/dYyNwq/KabbtKSJUv07LPP6rnnnstznqeeekrjx49XlSpVdOjQIZ5vAsv8+VZTVtwq2hvBkDFYBct7Gyw5g5Gv3lun06lrr71WGzZs0KZNm9SgQQNJUk5Ojlq3bq0dO3YU6vZ41IHAdffdd2v27Nluz+9wOLR3717VqFHDj6lglZycHNWsWVMHDx50e5nq1atrz549ioyM9GMyuMuX7as/+wJ/9gP0hSiMlJQUValSRcnJyW4vc/XVV+vHH3/0YypYac2aNWrTpo1HywwfPlwTJkzwUyJ4in4AhdWlSxctWLDA7fmjoqJ06NAhVaxY0Y+pYJWMjAxVrVrVo+f8Nm7cWFu3bmUMACHBGKOrrrpKGzZscHuZsmXL6vDhw9x9JkScOnVKVatW9eiOfDfddJMWL17sx1T+ETL3u9m8ebMkqVmzZhdNS05O1tNPP61XXnlFkjRhwgQ6LABA0Hnttde0evVqjR071nWCRvrj9mczZ85Udna2+vfvrxD5TRr+5MUXX/TohMuYMWMYtA4hkZGRmjJlituD0JGRkZo8eTKD1iEq3PuCcN/+cBUbG6u3337b7flLlCjh0fwIfK1bt9bAgQPdnr9OnTp68skn/ZgIdqEfCF+vvvqq4uLi3J5//PjxDFqHkOjoaE2ePNntc/rR0dF69913GQNAyHA4HHrnnXdUvHhxt5f597//zaB1CClXrpxee+01t+cvXbq0Xn/9dT8m8iMbr/b2mUOHDrkue2/WrJn529/+5nrVqlXLFClSxEgyxYoVM++8847dcRGGguFWU8GQMVgFy3sbLDmDkS/e2x07dpjo6GjTunVr1+3w/qqwt8ejDgS2rVu3mipVqlzyFkAjR440TqfT7rjwg08//dRER0cXuP+jo6PNZ599ZndU/IWv2ld/9wWBfqtw+kJMmjTJREREFNgOlipVynz33Xd2R4UfZGVlmX79+l3yWKhBgwbmwIEDdsfFX9APwBfWrVtnLrvssku2Ay+88ALfiULUzJkzXef683vFxMSYhQsX2h0V8Itvv/3WlCxZssDPQGRkpJk6dardUeEnEyZMuGQ/WK5cObNmzRq7o3otJG4VvmDBAnXp0uWiv0dERKhUqVKqXbu2brjhBj344INcfQRbBMOtpoIhY7AKlvc2WHIGo2B5b4MlZzg7ceKEpkyZoilTpujw4cO5pnXr1k3Dhg1Thw4dbEoHK+zcuVNvv/22PvjgA6WkpLj+Hhsbqz59+mjYsGG5rj5CYAiW9jUYbhXub8GSM5xt2LBBb7/9tj755BNlZGS4/h4XF6cBAwZoyJAhql69un0B4VfGGM2fP1+TJk266LaHtWvX1pAhQzRgwACVKlXKpoTIT7C0r8GSM5z9/vvvevfddzV16lQdO3Ys17SePXtq2LBhatu2rU3pYIUtW7bo3//+tz766COlp6e7/l66dGn169dPQ4cOVe3atW1MCPjX/v379c477+i9995TUlJSrmn9+vXTsGHDlJCQYFM6WGHNmjV6++23NXfuXGVlZbn+XqFCBQ0cOFBDhgxRlSpVbExYOCExcA0EumD44hMMGYNVsLy3wZIzGAXLexssOSFlZ2dr/fr16ty5s06dOqX4+Hj9/vvvdseChZKTk1WrVi2dPHlS5cuX1759+1SyZEm7YyEfwdK+MnAdPDnxxzPe6tev72oHDx486NGtExH89u/fr6uuukqnTp1ShQoVlJiYqIiIkHkiXsgJlvY1WHJCyszM1Pr169W1a1e+E4WppKQk1a1b13UscODAAT6zCCvp6elat26dunXrptOnT6tSpUo6evSo3bFgoePHj6tRo0audvDw4cOKjo62O1ahcUQPBJiHHnpI1atXl8PhcD27/a9mzJih5s2bu17ly5fXHXfcIUn65ptvck2rXLkyv7AKMu7UAafTqeHDh6tx48aqX7++BgwYoMzMTNf0l19+WQ0bNlTz5s3VunVrrV271qL0AKxQpEgRtW7dWsWKFZMknmUchkqVKuX6MhIdHc2gNYCwU65cuVztIIPW4adGjRquY6GiRYsyaA2EmaJFi6pNmzZ8JwpjZcuWzXUswKA1wk3x4sV13XXXuY6DORYKP5dddlmudjAUBq0lBq6BgHPnnXdq5cqVqlatWr7z9OvXT5s3b3a94uPjdc8990iSbrnlllzTEhISXNMQHNypA9OnT9fGjRu1ceNG7dy5UxEREZo4caIkafPmzZo0aZLWrl2rzZs3a9iwYRo2bJhV8QEAAAAAAAAAADzGwDUQYK677jpVrVrV7fl/+uknHT9+XF27dr1o2tGjR/Xdd9/p3nvv9WVE+Jk7dWDLli268cYbVbRoUTkcDt1666364IMPJEkOh0NZWVlKTU2VJJ05c8ajOgUAAAAAAAAAAGC1InYHAFA406dP17333quoqKiLps2cOVOdOnXSZZddZkMy+FOLFi00ZcoUDRs2TMWLF9enn36qAwcOSJKaNWumRx99VDVq1FBcXJyio6P1ww8/2BsYAAAAAAAAAACgAFxxDQSx1NRUffLJJxowYMBF04wxeu+99/KchuDXt29fdezYUe3atVO7du1Ut25dFSnyx2+R9u/fry+++EJ79uzR4cOH9eijj6pnz542JwYAAAAAAAAAAMgfA9dAEPvss8/UqFEjNWzY8KJpy5cv1/nz53XLLbfYkAz+5nA49Nxzz2nTpk1avXq1GjZsqEaNGkmSPv/8czVp0kSVK1eW9Mcz0VetWqXMzEw7IwMAAAAAAAAAAOSLgWsgiE2fPj3fK6qnT5+uvn37KjIy0uJUsML58+eVlJQkSTp58qTGjx+vESNGSJJq1qypVatWKSUlRZK0YMEC1a1bV0WLFrUtLwAAAAAAAAAAQEF4xjUQYB588EF99dVXSkxM1C233KKSJUtqz549uv/++9W1a1d17dpVkvTLL79o8+bNWrhw4UXrOHv2rL744gtt27bN6vjwAXfqwNmzZ9W+fXtFRETI6XTq4YcfVpcuXSRJ3bp107p169SyZUtFR0crJiZGH3/8sc1bBQAAAAAAAAAAkD8GroEAM2XKlDz/Pm3atFz/rlevns6dO5fnvKVLl1ZqaqrPs8Ea7tSBihUraufOnXnO53A49NJLL+mll17ySz4AAAAAAAAAAABf41bhAAAAAAAAAAAAAABbMXANAAAAAAAAAAAAALAVtwoHLBaot/AO1FyhJpDf50DOFkoC+X0O5GwAECoCua21KhvvAQCEt0BuawM5GwAAQDhg4BqwWMWKFe2OABux/0EdAIDwRj/AewAA4Y5+AAAAAPnhVuEAAAAAAAAAAAAAAFtxxTVggRIlSiglJcXuGG4rUaKE3RFCSrDtf4k64GvUAQAIb/QDvAcAEO7oBwAAAOAOBq4BCzgcDsXExNgdAzZh/4M6AADhjX6A9wAAwh39AAAAANzBrcIBAAAAAAAAAAAAALZi4BoAAAAAAAAAAAAAYCsGrgEAAAAAAAAAAAAAtmLgGgAAAAAAAAAAAABgKwauAQAAAAAAAAAAAAC2YuAaAAAAAAAAAAAAAGArBq4BAAAAAAAAAAAAALZi4BoAAAAAAAAAAAAAYCsGrgEAAAAAAAAAAAAAtmLgGgAAAAAAAAAAAABgKwauAQAAAAAAAAAAAAC2YuAaAAAAAAAAAAAAAGArBq4BAAAAAAAAAAAAALZi4BoAAAAAAAAAAAAAYCsGrgEAAAAAAAAAAAAAtmLgGgAAAAAAAAAAAABgKwauAQAAAAAAAAAAAAC2KmJ3AAAAENqMMUpLS7M7hkdKlCghh8Nhd4yQEIz7/wLqAQBfCca2kDbQd4Jx/19APQB8IxjbAT7/vhOM+/8C6gF8JRg/B9R/36IOwF0MXAMAAL9KS0tTbGys3TE8kpKSopiYGLtjhIRg3P8XUA8A+EowtoW0gb4TjPv/AuoB4BvB2A7w+fedYNz/F1AP4CvB+Dmg/vsWdQDu4lbhAAAAAAAAAAAAAABbccU1AACwzLFjxwL2l4qpqamqWLGi3TFCWiDv/wuoBwD8LZDbQtpA/wvk/X8B9QDwr0BuB/j8+18g7/8LqAfwt0D+HFD/rUEdQEEYuAYAAJaJiYkJ2ANT+B/7HwBoC8Md+x8A7UB4Y/8DfA5AHUDBuFU4AAAAAAAAAAAAAMBWDFwDAAAAAAAAAAAAAGzFrcIBAACCxLFjx7RhwwZt2LBBe/fu1enTpyVJZ86c0XvvvacWLVqoYcOGioqKsjkp/MEYo3379rnqwIX9f/r0aY0cOVItWrRQixYtVLNmTTkcDpvTAoDvOZ1O/frrr9qwYYM2btyYqx185plnXO3g5ZdfTjsYojIyMrR161Zt3LhRW7ZscdWBpKQkjR8/3lUH4uLibE4KwF+OHj3qOh7ev39/ru9EM2fOVIsWLdSgQQMVKcJp71DkdDq1Z8+ePI8Fnn76aVc/UK1aNY4FEJKMMTp48KCrHTx8+LDrc3D27Fl99NFHatGiherWrauICK5bDUU5OTnatWuXNmzYoE2bNuVqB8eMGaOWLVuqRYsWqly5ctC2gw5jjLE7BAAACF2pqamKjY2VJKWkpATsM2wCNWdGRobmzp2rSZMmafXq1Zecv2zZsurXr58GDRqkOnXqWJCwYIH6vuYnEPOeOnVKM2bM0Lvvvqt9+/Zdcv5atWpp0KBB6tevn8qVK2dBQgCXEohtS14CNefRo0c1bdo0TZkyRUePHr3k/I0bN9aQIUPUu3dvlSxZ0oKEBQvU9zU/gZh369atmjRpkj788EOlpqZecv6bbrpJQ4YMUefOnRm8QkAIxM9VXgI1Z3p6uubMmaNJkyZp3bp1l5y/fPnyuv/++/Xggw+qevXq/g94CYH6vuYnEPMeP35c06dP1+TJk3Xw4MFLzl+/fn0NHjxYffr0UZkyZfwfEJcUiPUqL4GaMzk5WR988IEmTZqkHTt2XHL+KlWq6MEHH9TAgQMVHx9vQcJLC9T39q8CNefBgwc1depU/ec//9Hx48cvOf+VV16poUOH6q677lKJEiUsSOg7DFwDAAC/CtQDvr8KtJzGGH388cd67LHH3DogzUuPHj309ttv67LLLvNxOvcF2vt6KYGUNyMjQ88//7xeffVVZWRkeLx8sWLFNHz4cI0ePVpFixb1Q0IA7gqktqUggZYzJSVFTz75pKZMmaLs7GyPly9ZsqSee+45Pfzww4qMjPRDQvcE2vt6KYGU98CBAxo0aJC++eYbr5avVq2aJk2apE6dOvk4GeCZQPpcFSTQchpjNG3aND355JOuK8o84XA4dO+99+qNN96w9U4Mgfa+Xkog5U1LS9Po0aP19ttvKysry+PlS5QooVGjRmnEiBH8kMlmgVSvChJoOXNycvT6669r7NixSklJ8Xj5IkWKaPDgwXrppZds35ZAe2/zE2g5z5w5o+HDh2vGjBlyOp0eL1+mTBmNGzdOgwYNCpqr8IMjJQAACDvLli2Tw+HI9YqNjVVCQoLeeOMNr05gB4vExETdfvvt6t27t9eD1pL06aefqlGjRvr00099mM464VwHNmzYoJYtW2rcuHFeDVpL0vnz5/XCCy+oZcuW2rhxo48TArBCOLeDy5YtU9OmTfXOO+94vZ3nzp3T448/ruuuu06//vqrjxNaI1zrgDFGkydPVpMmTbwetJak3377Tbfddpv69eunM2fO+C4gYJFwbQMk6dChQ+rYsaMeeOABrwatpT/aklmzZqlRo0b68ssvfZzQGuFcB1avXq3mzZvr9ddf92rQWvpj4HvUqFG65pprtH37dh8nhBXC+TOwa9cutW3bViNGjPBq0FqSsrOz9fbbb6tp06Zavny5jxNaI5zrwKJFi9S4cWNNnz7dq0Fr6Y+B76FDh+rGG2/UgQMHfBvQT/iZEQAACGh33XWXOnXqJGOMEhMTNWvWLD322GPauXOnpk6danc8n9u1a5duuukmHT582CfrO3nypHr27Knt27frX//6V1A+3ybc6sC8efPUs2dPZWZm+mR927ZtU5s2bfTpp5+qa9euPlknAGuFWzs4ffp0PfDAA16fnPmr1atX66qrrtJXX32ltm3b+mSdVgunOpCTk6NBgwZp2rRpPlvnzJkztXbtWi1evFhVqlTx2XoBq4RTGyBJW7Zs0c0331yoH/H+WWJiorp27aqXX35ZI0aM8Mk6rRZudeDjjz9Wnz59lJOT45P1rV+/Xq1bt9a8efN04403+mSdsFa4fQaWL1+uLl266Ny5cz5Z3759+9ShQwdNnz5dffv29ck6rRZudWDixIl65JFHfLa+pUuXqmXLlvr666/VsmVLn63XLwwAAIAfpaSkGElGkklJSXF7uaVLlxpJZsKECRetr2rVqsbhcJjjx4/bntOXdu/ebeLj4105fP165plnLN+mwryvVteBC+u2sx7873//M5GRkX7Z/0WKFDELFiywfJsA0Bd6Yvr06X7rB0uUKGFWr15t+TbRF7ovJyfH9OvXz291oE6dOiYxMdHSbQKMoR/wxLZt20xcXJzf2oGXX37Z8m2iH/DMxx9/bBwOh1/2f3R0tPn+++8t3ybQDnpixYoVpnjx4n5rB2fOnGnLdlEH3Pfmm2/6bf+XKlXKbNy40Zbtche3CgcAAEElJiZGrVu3ljFGe/futTuOz6SmpqpTp05KTEz0WxkvvPCCPvjgA7+t3yqhWge2b9+uHj16+Oyqgr/Kzs7WnXfeqR07dvhl/QCsE6rt4A8//KCBAwf6bf1paWm67bbbdPToUb+VYZVQrQOvvPKKZsyY4bf17969W7fffntI31IS4SFU24AzZ87o1ltv9frW4O4YOXKk/vvf//pt/VYJ1Tqwbt069enTR8YYv6w/IyNDXbt21f79+/2yflgnVD8Dhw4dUpcuXZSenu63Mvr3769Vq1b5bf1WCdU6sHDhQp9eaf1XycnJ6tSpk06dOuW3MgorJAeuf/vtN40fP17XXnutLr/8ckVHR6t06dJq0aKFnnjiCe3atcvuiAAAoBAuHJDGxcXZnMR3nn76ae3evdujZdatW6dDhw5p3bp1bi/z0EMPhcQJ+1CrA9nZ2erXr59Hz7P2Zv+fP39e/fr144Q9EAJCrR1MTU1Vv379PLo9uDftYFJSkh588EG/nRC3UqjVge3bt2vMmDEeLeNNHfjxxx/1+uuvexoPCDih1gZI0qOPPurxI5O8aQcGDRqkkydPehov4IRaHTh//rz69u3r0XcVb/Z/SkqKBgwY4LNHksA+ofYZMMZo4MCBOnPmjEfLefo5cDqd6tevn9LS0rxIGVhCrQ4kJSV5/ENeb9rBxMRE/fOf//Q0nmVCauA6OztbY8eOVf369fXUU09p9erVkqSmTZsqLi5OGzdu1KuvvqrGjRtrwoQJNqcFAADuSEtL08mTJ3XixAlt27ZNQ4cO1aZNm9SqVSvVrVvX7ng+sWLFCr311lseLxcfH6+qVasqPj7e7WXOnDmjBx980OOy7BQOdeDVV1/V+vXrPVrGm/0vSWvXruWEPRBkwqEdfPrpp7Vv3z6PlvG2HVywYIE+/PBDj5axW6jXgZycHPXt21dZWVkeLedtHRgzZox27tzp0TKAnUK9DZCkRYsWaebMmR4v5007cPz4cT300EMel2WncKgDzz//vMd3h/K2H1i6dKmmTJni0TKwVzh8BmbOnKlvvvnG4+W8+Rzs3r1bo0eP9rgsO4VDHXj88cc9vtjE23Zw9uzZmjdvnkfLWMamW5T7XFZWlrnzzjtdz6oYPXq0+f3333PNc+DAAfPoo4+aqKgoM3LkSJuSAgAQXgr7DJu8XnfcccdF/bxdOX3hpptu8uq5NIcOHTLGGHPo0CGPl12/fr0l2+aL57lZVQcKm9dbqamppmzZspbu/7i4OJOWlmbJ9gGgL7yUxMREExUVZWk7WLt2bZOTk2PJ9tEXXtr8+fMtPxbq27evJdsGGEM/4I7WrVtb3g788ssvlmwb/cClJSUlmRIlSli6/6tUqWKysrIs2T7QDl5Kdna2qVGjhqXtYHR0tDlx4oRl20gdKNj+/fuNw+GwtB1s3ry5cTqdlm2ju0Lmiuthw4Zp7ty5Klu2rJYsWaKxY8de9AuDatWq6fXXX9eSJUvUoUMHm5ICAABPPPDAA/r222+1cOFCvfzyy4qLi9Phw4dVrFgx1zy9evVSjx49ci13+vRpVapUSR999JHVkT2ye/duffvtt5aX++6771peprdCvQ7MmTNHSUlJlpZ5+vRpffrpp5aWCcB7od4OTp8+3eMrbQtrz549WrJkiaVlFkao14FJkyZZXuYnn3zi1+foAr4U6m3Axo0b9eOPP1pe7uTJky0v01uhXgdmzZpl+W2Ljxw5oi+//NLSMuG9UP8MfPPNN5Y/ez0jI0MzZsywtMzCCPU6MGXKFMsfZ7R582b99NNPlpbpFrtHzn1h4cKFrl8IzJs3z+44AADgTwr7i8oJEybk+vuqVauMw+EwPXv2dP3t1KlTpnLlyubjjz92/a1Xr17mzjvv9HvOwho+fLhXv6hVIX9VWbx4cXP27Fm/b58vri6wqg4UNq+3WrVqZfn+l2Rat25tyfYBoC8siNPpNNWqVbOlHezWrZsl20hfWLB9+/bZciwkybz++ut+3z7AGPqBS3nggQdsaQfKlClj0tPT/b599AOX1qBBA1v6gZtuusmS7QPt4KV06dLFlnawZs2all1xSx3IX2Zmpilfvrwt7WCfPn0s2UZPBP0V18YYjRgxQpLUu3dv/f3vf7c5EQAA8Kc2bdro3nvv1Zw5c7R69WpJUlxcnKZPn65hw4bp6NGjmjt3rpYtWxYUv6BfunSpLeWmp6dr7dq1tpRdWKFUB86dO6d169bZUvbatWuVmppqS9kACieU2sH9+/frt99+s6XsZcuWWX5Vg6+EUh1Yvny5bWUvW7bMtrKBwgilNkCy7zvRmTNntGXLFlvKLqxQqgPHjh3Tzp07bSl75cqVys7OtqVsFE4ofQacTqdtxyT79u3ToUOHbCm7sEKpDvz88886efKkLWUH4vFw0A9cL168WNu3b5ckPfPMMzanAQAAVhg9erQiIyM1ZswY1986duyoHj16qHfv3hoyZIimTZumcuXK2Zjy0jIyMrR161bbyt+wYYNtZRdWqNSBzZs32zZo4nQ6g/ZEHYDQaQft7IuSkpJ04MAB28ovLOpAcJcNFFaotAFnz57V7t27bSs/mNuBUKkDdu6D9PR02wbNUXih8hnYs2ePzp07Z1v5tIP2s3MfHDx4UCdOnLCt/LwE/cD1vHnzJElXXXWV6tWrZ28YAABgidq1a6tXr1767rvvtGLFCtffX331Ve3Zs0e33nqrbrvtNhsTumfHjh2WP9Pzz4L5y0mo1IFNmzbZWn4w1wEg3IVKO7hx40Zbyw/mdjBU6oCdfeGRI0d07Ngx28oHCiNU2oDNmzfbWj79gP3sPhawu3x4L1Q+A5wX8F6o1AG72yG7y/+rInYHKKwLl7F36NDBkvJatmypxMRES8oCACAUOJ1Ov6x31KhRmj17tsaMGeO6tVxMTIxq1qypJk2aFGrdderUUUSE/3/fd/78+XynrVu3TvHx8QUuf2F6fHx8gbd2SkxM1FVXXXXR3//3v/+patWqbqb1jr/2v+TfOiBZUw+Sk5PznXapOuDu/pfyrwOjRo3Syy+/7GZaAN6iL8xfUlJSvtOsaAcHDhyoRx55xL2wXqIvLFhBA8dW1IGmTZsqKirKzbSAd+gH8peenp7vNCu+E3388cf65ptv3EzrHfqBgp05cybfaVb0Aw899JBGjRrlXlh4jXYwfykpKflOs6IdfOONNzRz5kz3whYCdSB/p0+fzneaFe3gXXfdpRIlSriZ1j3x8fFav369V8sG/cD1wYMHJUkNGjSwpLzExEQdOXLEkrIAAAhn7du3L/AWyg0aNFBOTo5fyv7999/9sl5PxMfHuz2oXKRIEa8GoDMzMwP6uMbOOiDZXw/crQPe7n/pj2ds23lLMgAFoy/0fzt45syZAk+Y242+0P914Pjx414tB1iBfsD/34nOnz/Pd6IC2F0PrOgHkpOTC/xBMexFO+j/djAtLU1paWkeL2cV6oD/28GkpKQCf1BstaAeuM7MzHR9oEqXLm1JmZf6dQsAAMjN6XQGxIGeJypVqmTZFdenTp3Kc5o7d3iJj49XkSJFlJ2dXeD8+U2Ljo5W+fLl3QvrpWDc/xdYUQ+Sk5PzHTi+VB1wd/8XtK6SJUuqVKlS7oUF4LVgbAut6guTkpLyPVFmRTtYpkwZxcTEuBfWS8G4/y+woh4cO3ZM2dnZeU6zog5cdtllXHENvwvGdsCqfiA9PT3fK82s+E5UrFgxvz/7NBj3/wVW1IMzZ84oNTU1z2lW9AOlS5dWbGyse2HhtWD8HFjVDqakpOjs2bN5TrOiHYyJiVGZMmXcyloY1IH8nT59Ot87kFjRDpYtW9YvV1x7K6gHrosWLaqiRYsqMzOzwNsp+JK3l7YDABCuUlNTg+5L4O7du/1+EvtCOXXr1s1zWl637vmrQ4cOqWrVqkpMTNTll1/ucfn33HOPpk+f7vFyngjG/X+BFfVg+vTpuv/++/Ocdqk6UNj9L0lvvfWW+vbt69WyANwXjG2hVX3huHHj9Mwzz+Q5zYp28JNPPtEtt9zi1bLuCsb9f4EV9eCWW27R4sWL85zm7zrgcDi0Z88elSxZ0uNlAU8EYztgVT+wYcMGtWzZMs9pVnwnevDBB/Xmm296vJwngnH/X2BFPXjzzTf16KOP5jnNimOB6dOn6x//+IdXy8J9wfg5sKodXLBggbp06ZLnNCvawdGjR2vkyJEeL+cp6kD+Ro4cqVdeeSXPaVa0gwsWLFCbNm28WtYfgnrgWpIaN26sjRs36scff1Tv3r3tjgMAAGy2bNkyuyO4rVatWipVqpRttyVr0aKFLeX6WzDVAbv3gd3lA/AP2sHgKd9fgq0O5Ddw7W/16tVj0BohKZjagMaNGysqKkpZWVm2lE8/YD+794Hd5cM/+AwET/n+Qh1wT0REhJo3b25b+Xnx/zXuftajRw9J0gcffKDDhw8XOO/hw4ddD2cHAACwW0REhBISEmwrP1S/nASTRo0aKTo62payixcvrgYNGthSNgBcYGdfdMUVV/j9kRm4NDvrAMdCgP2io6PVpEkT28qnHbDflVdeKYfDYUvZcXFxqlatmi1lAxdUqlRJlSpVsq18O89L4Q929kUNGjTw+W3CCyvoB66HDh2q2rVrKzk5WR06dNDy5ctzTTfGaNu2bXr00UfVoEED265oAgAAyEt+t4Pyt0qVKvHlJABERUWpY8eOtpTdsWNHFSkS9DdgAhDkKlSooKuvvtqWsjt37mxLucitQ4cOKl68uC1lUweAwGDXd6JatWqpfv36tpSN/xMbG6vrr7/elrI7d+5s26A58Gd2tYNt2rRRXFycLWXj/9SsWdO2CwsC8Xg46AeuY2NjtWjRIjVo0EC7d+9W+/btVbFiRbVq1UrNmzdX+fLl1bRpU7355ps6f/68W88EAAAAsErfvn1VrFgxy8sdOHCgoqKiLC8XFxsyZEhYlQsAfzV48GBbyqUdDAxly5bV3XffbXm5FStW1B133GF5uQAuNnDgQEVGRlpe7uDBgxUREfSnx0OCXX3y0KFDbSkX+Cu7jof5DAQGh8NhSzvocDg0aNAgy8u9lJDomWvXrq2NGzfq3Xff1Y033ihjjDZt2qTdu3crLi5Od9xxhyZNmqQjR46ocuXKdscFAABwiYuLU69evSwtMzIyUgMHDrS0TOTvxhtvVO3atS0ts27duurQoYOlZQJAfnr06GH5lR7XXnutGjVqZGmZyJ8dJ+oGDBigokWLWl4ugItVqVJFt99+u6VlFitWTH379rW0TOSva9eulp+3b9GiBRe5IWA0b95cbdq0sbTMChUq6B//+IelZSJ/9957r2JiYiwts1OnTqpevbqlZbojJAaupT8ONgYNGqRvv/1Wx48fV1ZWllJTU7V79259/vnnGjx4sC677DK7YwIAAFzkueeeU8mSJS0r7/HHH1fVqlUtKw8Fi4iI0GuvvWZpma+99hpXlwAIGMWLF9eLL75oWXkRERF69dVXLSsPl5aQkKB77rnHsvIqV66s4cOHW1YegEsbN26cpXeieuaZZ1SuXDnLykPBoqKi9Morr1ha5quvvsptwhFQJkyYYGmdHD9+vKKjoy0rDwUrXbq0nnvuOcvKK1q0qMaPH29ZeZ7gbBUAALDd7t271aZNG9WtW1dXXXWVfv7554vmWbNmjZo3b67mzZurUaNGevDBB5WRkSFJ+v7779WqVSs1bNhQjRo10ogRI+R0OiVJ27Zt03XXXaf69eurcePG6t+/v9LT0y3dvkupVq2aZSfQGzRooH/961+WlAX3de3aVb1797akrD59+gTkM4yAcOdOX+h0OvXYY4+pYcOGatq0qa6//nrt2bNHkrR//361aNFCzZs3V+PGjdW9e3clJSVJkg4cOKDIyEhXP9q8eXPt3bvX0u27lAceeEA33nijJWU98cQTatWqlSVlwX1vvfWWKlasaElZU6dOVdmyZS0pC3CXu/3A8OHD1bhxY9WvX18DBgxQZmama/qECRPUuHFjNWzYUN26ddOZM2dc0+68805VrlxZDocj198DRb169fTCCy9YUlaLFi00cuRIS8qC++6++279/e9/t6SsYcOGqX379paUBe889NBDql69uhwOhzZv3pznPAWdCwqG49+/atOmjR577DFLyurYsaP69etnSVlw36OPPqrWrVtbUtazzz6rxo0bW1KWpxi4BgAAtnvwwQf1wAMP6Ndff9XIkSPzvGVbs2bNtG7dOm3evFnbtm3T8ePHNWnSJEl/PBvxk08+0Y4dO7RhwwatXr1as2bNkvTHXVn+/e9/a9euXdqyZYtSU1P18ssvW7l5bhk4cKC6dOni0TKJiYk6fPiwEhMT3Zq/WLFiev/99215pran3DlxJ0k333yzmjZtqubNm+vaa6/Vpk2bXNPOnz+v22+/XXXr1lWzZs100003uQZ4Cppml4kTJ6patWpuz+/p/pekGjVq6M033/QiHQB/c6cvnD9/vlatWqUtW7Zo69atuuGGG/T0009L+uMK0pUrV2rz5s3avn27KleunOsX+yVLltTmzZtdr1q1alm0Ze5xOByaNm2aypcv7/Yy3rSDV155paVXMhSGu33hwoULlZCQ4PrRwvvvv++aVlA/Wb16ddWrV891MnfOnDl+36aCxMXF6b333vPojiDe1IEHHnhAt912mzcRAb9ypx+YPn26Nm7cqI0bN2rnzp2KiIjQxIkTJUnffvutZsyYoTVr1mjHjh1q0aKFRo0a5Vp20KBB+Q7+BIpHHnnE48FET9uB2NhYzZw5U0WKFPEiobXc7Qe+/vprtWzZUk2bNlXr1q21ZcsWt6e7MzhoFYfDocmTJ6tSpUpuL+NNP9CgQQO99NJL3kSEhe68806tXLmywO/IBZ0LkgL/+Dcvzz//vJo2berRMp5+DipUqKCpU6cG5R0H3GmzLnW+J9COgf8sMjJSM2fOVKlSpdxexpt2sE2bNhoxYoQ3Ea1hAAAA/CglJcVIMpJMSkrKRdOPHTtmSpYsabKysowxxjidTlOxYkWze/fufNeZnp5ubrnlFvPGG2/kOX3o0KHm2WefzXPahAkTzH333edxTiukpqaaa6+91pXDl6+oqCizYMECy7fJ2/f1+uuvNzNmzDDGGPPZZ5+Zli1b5jlfUlKS6/+/+OIL07RpU9e/09PTzVdffWWcTqcxxpi3337btGvXrsBpdteD3bt3m0qVKvmlDlSuXNns2bPH8m0C4Lu+cN68eaZZs2YmOTnZOJ1O88QTT5hHH330ovVlZ2ebAQMGmIcfftgYY8z+/ftN6dKlC53TChs2bDClS5f2SztYr149k5iYaPk2+bMvdDqdpmzZsmbLli3GmD/2dXR0tElOTjbGFNxPVqtWzWzatMlneX1l+vTpftn/kky3bt1cnzPASr7qB4YOHWrGjRvn+vfnn39umjRpYoz543vOwIEDXdM2bNhgSpYseVFZknK1DZ7ktMKZM2dMy5Yt/dIGFCtWzHz//feWb5M/+4HTp0+buLg4s337dmOMMT/88INp1KiR29OXL19uDh06lKtPsLsebNu2zZQrV84vdaB69erm0KFDlm8TvK9X+R2v5OXP54LcPf71VU5fOnr0qKldu7ZfPgNlypRx+/30NV+8t3m1WX9V0LkgYy5dpwKhDqxYscLExMT4pQ40a9bMnD592pbtchdXXAMAAFsdOnRIlSpVcv3i3eFw6IorrtDBgwcvmvfAgQNq1qyZypcvr9KlS2vIkCEXzZOYmKi5c+fmeSvk1NRUTZs2zbLbj3mqRIkSWrRokTp27OjT9cbExOjLL78MmquLjh8/rvXr17tunf2Pf/xDhw4dyvOK6DJlyrj+/+zZs7l+MVysWDF16tTJ9bfWrVvrwIEDl5xmp9q1a2vlypU+/yV47dq1tWLFiqD4hTkQjtztC7t06aL27dsrPj5elSpV0nfffaexY8e6pmdmZqp58+YqX768du/enevREKmpqbrqqquUkJCgsWPHKicnx5qN81BCQoKWLVvm81tGX3nllVq+fLllt6IuLE/6wj/f9jc5OVnlypVzPa+woH4yUPXv318ffvihoqKifLree++9V3PmzAmKqywRftztB1q0aKH58+crOTlZWVlZ+vTTT13HsC1atNCSJUuUmJgoY4w++ugjnTt3TqdPn7Z6cwqldOnSWrJkidq1a+fT9ZYpU0aLFy/W9ddf79P1+ou7/cDevXtVrlw5NWrUSJJ07bXX6uDBg9q4caNb06+77jpVrVrVqs1yS+PGjfXDDz/4PFejRo20cuXKgNte+EZe54KC5fj3rypVqqQVK1aoWbNmPl/v8uXL1bx5c5+u10rutFmBer7HE23bttWSJUsUFxfn0/W2adNGS5cuDfhH5jBwDQAAgkb16tW1ZcsWJSYmKiMjQ1988UWu6cnJyerSpYtGjBihli1b5pqWmZmpnj176uabb1a3bt2sjO2RmJgYffXVV5o4caKKFy9e6PV16NBB27dv1y233OKDdNbw5McM0h/PbL788ss1evRoffDBB/mud+LEifn+aKGgaVarWbOmNm3apEGDBvlkfUOGDNGmTZtUs2ZNn6wPgH3Wr1+v7du368iRIzp69KhuuOGGXG1F0aJFtXnzZh07dkz169fXlClTJP1xkurIkSNat26dlixZohUrVui1116zazMuqXnz5tq+fbt69epV6HVFRkbq6aef1po1a4Jm0Fpyvy90OByaM2eO7rjjDlWrVk1t27bV+++/r6JFi7rmKaif7NOnj5o0aaIBAwboxIkT/t8wN91zzz3auHGjWrRoUeh1lS1bVh988IHef/99nw+GA1br27evOnbsqHbt2qldu3aqW7euq524/vrrNXz4cHXu3FmtW7dWhQoVJCkof6xxYfB6/Pjxudozb3Xq1Enbt2/Xtdde64N01nC3H6hTp45OnTql1atXS/rjsSLnzp1zDdJcanqgatiwobZu3ar77ruv0OtyOBx6/PHHtW7dOlWpUsUH6RBo8joXFGzHv38VHx+vH3/8UU8++aRHj1HJz913361t27Z5fBvyUJDX+Z5APQb+s9atW2v79u3q2rVrodcVFRWl559/XsuWLQv4QWtJ4lbhAADAr/xxq3BjjJk9e7bp3Lmz69/JycnmmmuuMc8///xF82ZmZprbb7/d3H///a5bBXma0w67d+82d9xxh4mMjPT41j81atQwU6ZMMTk5ObZuQ17va+vWrU25cuXyfB08eNCsX7/e1K1bN9d6rrrqKvPdd98VWNbMmTPNrbfemue0cePGmdatW5vU1NQCpwVaPViyZIlp3bq1V7d/uuaaay75ngGwhr9uEbt9+3ZTuXLlPMtcs2aNady4cZ7TPv7441x9qLs57fDFF1+YJk2aeNUOdujQwaxbt87uTfBrX5iVlWXatWtnli9fbowxZu3atSY+Pt6cOHHiohx/7Sd/++03Y8wfx0kjRoxwTQukepCZmWleeeUVrx6jUbRoUdOnTx9z9OhRW7cBMMa/34natm2b57Q1a9aYqlWrXvR3KbBvFf5XP//8s+ncubNxOBwetwN16tQx77//fr7fAa3i7+9E33//vbnuuutMQkKCGTp0qGnYsKH53//+5/Z0Y0xA3Sr8rxYsWGASEhK8OhZo166dWbVqld2bAOO/W4UXdC7oz/I7/vVVTn/66aefzPXXX+/1baHnzZtn9yYYYy793l6qXfwzd28hn9e5oPyOgd3NaTWn02k++ugjU79+fY/3v8PhMB07dnQ9VihYMHANAAD8yp0Dvnbt2uV6dleLFi0ummf37t0mMzPTGGNMRkaG6dGjh3n66aeNMcacO3fOtGnTxvzrX/+6aLmsrCxzxx13mP79+xd4wiLQDkz/7NChQ2b06NGmYcOGJiIiIt8D0ri4OPP3v//dfPXVV7YPWF/gzfvq7Yk7Y4wpVqyYOXnyZK6/TZgwwbRo0SLPE3R/nRao9WDDhg3m/vvvN1dccUWBX0qqVatmBg4caDZu3Gh3ZAB/4qu+8LXXXjM33HCDycjIMMYYM378eHPzzTcbY4w5cOCA64RMTk6Oeeyxx8zdd99tjPmjXb3Qh54/f97ceeedZvTo0V7ltIPT6TQrV64099xzzyUHMGvXrm0eeughs3PnTrtju/izL1y3bp2pU6dOrr+1bNnSLF68OM/15tVPGvPHsxRjY2O9zutvmZmZ5tNPPzW33HKLKVWqVL77v0iRIqZ58+Zm/Pjx5vjx43bHBlx81Q+kp6e7nkt54sQJ06xZMzN//nzX9As/1EhNTTU33XSTeeutty5aR7ANXF+wb98+8+STT5r69esXOIhdoUIFc+edd5pvv/02LL8TnT9/3pQpUybf+fKbHsgD18b8sf0//vijue+++0zVqlULPBaoUaOGGTJkiNm2bZvdsfEn/hi4LuhckLvHv77KaYUdO3aYf/7zn6ZWrVoFfgYqV65sevfubVatWmX7D3f+zJfvrTsD1wWdC7rgz8fA/sjpS06n03z//femR48epmLFigXWgbp165rhw4ebPXv22B3bK8F3rxgAABBypkyZor59++rFF19UqVKlNGPGDEnS/fffr65du6pr1676/vvv9dZbbykyMlLZ2dm64YYbNHr0aEl/3PZn7dq1Sk1Ndd0+vHv37ho1apTmzJmjL774Qk2bNtWVV14pSfrb3/6md955x56N9ULVqlU1duxYjR07Vqmpqdq8ebP27t2r8+fPKyoqSmXKlNGVV16patWqBcWzKy/lsssuU0JCgj788EP17dtXn3/+uapWraratWvnmu/MmTNKS0tT5cqVJUnz5s1TuXLlcj0D6PXXX9fs2bO1ZMmSXM/5vNS0QJOQkKD//Oc/kqQTJ05o48aNrlvmR0dHKz4+Xi1atFD58uVtTgrAW+70hUOHDtXOnTvVrFkzRUVFKT4+XpMnT5Ykbd26VaNGjZIkOZ1OJSQk6K233pIkrVy5UmPGjHH1oR06dHDNGwwcDof+9re/6W9/+5sk6ffff9fGjRt14sQJZWZmKjo6WpdffrkSEhICvj13l7t94eWXX67ff/9dO3fuVIMGDbRnzx7t3btX9erVK7CfTE1NVVZWluv9mj17tus4KRBFRUWpe/fu6t69u5xOp/bu3astW7YoOTlZOTk5Kl68uOrWraumTZuqWLFidscFvOJOP3D27Fm1b99eERERcjqdevjhh9WlSxfXOm6++WY5nU5lZmbq3nvv1bBhw1zTbrvtNm3ZskXSH8/6rVOnjpYtW2bpNhZGjRo19NJLL+mll17SuXPntGnTJu3fv18ZGRmKiopSXFycEhISVLVq1bD6TiT90S9WqlRJkvT888+rQ4cOuea71PRg4HA4dPXVV+vqq6+WJB07dkwbNmzQ8ePHXccCVapUUUJCgs+fCQt7PPjgg/rqq6+UmJioW265RSVLltSePXtytYkFnQsK9uPfvDRo0EBvvfWW3nrrLSUlJWnjxo06fPiwMjIyVLRoUVWoUEEtWrRQfHy83VFtl9/5nmA7Bv4zh8Oh66+/Xtdff72MMTpy5Ig2btyoU6dOKSsrS9HR0apevbquvPJKlSpVyu64heIwxhi7QwAAgNCVmpqq2NhYSVJKSopiYmJsTpS3YMkZbLx9X3/55Rf17dtXp06dcp24a9KkiaT/O3nXrFkzde/eXenp6YqIiFCFChX06quvqnnz5pKkw4cP6/LLL1fNmjVVsmRJSVJ0dLR++umnfKd9//331AMAPhcsfUyw5Aw2/uwLu3btqtmzZ+vFF190DWQ99dRTuvvuu/Xbb7/l20/u27dP//jHP5STkyNjjGrWrKmJEyeqevXq1APAD4LlcxUsOYONv/uBgQMHasWKFcrOztY111yjt99+O9dATUHT/zw4WK5cOZUsWVJbtmyhHsDngqV9CZacwcgX721ebdaePXsk/V+7mJCQkO+5oIKOgX2ZE4XDwDUAAPCrYDngC5acwSbY3tdgywsgOARL2xIsOYNNsL2vwZYXCAbB8rkKlpzBJtje12DLi+AQLPUqWHIGo2B5b4MlZyjjVuEAAMAyqampdkfIVyBnCxXB8B4HQ0YAwS2Q25lAzhYqguE9DoaMQDAL5M9YIGcLFcHwHgdDRgS3QK5jgZwtlATy+xzI2cIFA9cAAMAyFStWtDsCbMT+BwDawnDH/gdAOxDe2P8AnwNQB1CwCLsDAAAAAAAAAAAAAADCG8+4BgAAfmWMUVpamt0xPFKiRAk5HA67Y4SEYNz/F1APAPhKMLaFtIG+E4z7/wLqAeAbwdgO8Pn3nWDc/xdQD+Arwfg5oP77FnUA7mLgGgAAAAAAAAAAAABgK24VDgAAAAAAAAAAAACwFQPXAAAAAAAAAAAAAABbMXANAAAAAAAAAAAAALAVA9cAAAAAAAAAAAAAAFsxcA0AAAAAAAAAAAAAsBUD1wAAAAAAAAAAAAAAWzFwDQAAAAAAAAAAAACwFQPXAAAAAAAAAAAAAABbMXANAAAAAAAAAAAAALAVA9cAAAAAAAAAAAAAAFsxcA0AAAAAAAAAAAAAsBUD1wAAAAAAAAAAAAAAWzFwDQAAAAAAAAAAAACwFQPXAAAAAAAAAAAAAABbMXANAAAAAAAAAAAAALAVA9cAAAAAAAAAAAAAAFsxcA0AAAAAAAAAAAAAsBUD1wAAAAAAAAAAAAAAWzFwDQAAAAAAAAAAAACwFQPXAAAAAAAAAAAAAABbMXANAAAAAAAAAAAAALAVA9cAAAAAAAAAAAAAAFsxcA0AAAAAAAAAAAAAsBUD1wAAAAAAAAAAAAAAWzFwDQAAAAAAAAAAAACw1f8DXLuY8kkKUiwAAAAASUVORK5CYII="
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "theta_A = 2 * np.arcsin(np.sqrt(0.25))\n",
    "theta_B_nA = 2 * np.arcsin(np.sqrt(0.6))\n",
    "theta_B_A = 2 * np.arcsin(np.sqrt(0.7))\n",
    "theta_C_nBnA = 2 * np.arcsin(np.sqrt(0.1))\n",
    "theta_C_nBA = 2 * np.arcsin(np.sqrt(0.55))\n",
    "theta_C_BnA = 2 * np.arcsin(np.sqrt(0.7))\n",
    "theta_C_BA = 2 * np.arcsin(np.sqrt(0.9))\n",
    "\n",
    "qrA = QuantumRegister(1, name='A')\n",
    "qrB = QuantumRegister(1, name='B')\n",
    "qrC = QuantumRegister(1, name='C')\n",
    "\n",
    "# Define a 3-qubit quantum circuit\n",
    "qcA = QuantumCircuit(qrA, qrB, qrC, name=\"Bayes net big\")\n",
    "\n",
    "# P(A)\n",
    "qcA.ry(theta_A, 0)\n",
    "\n",
    "# P(B|-A)\n",
    "qcA.x(0)\n",
    "qcA.cry(theta_B_nA, qrA, qrB)\n",
    "qcA.x(0)\n",
    "\n",
    "# P(B|A)\n",
    "qcA.cry(theta_B_A, qrA, qrB)\n",
    "\n",
    "# P(C|-B,-A)\n",
    "qcA.x(0)\n",
    "qcA.x(1)\n",
    "qcA.mcry(theta_C_nBnA, [qrA[0], qrB[0]], qrC[0])\n",
    "qcA.x(0)\n",
    "qcA.x(1)\n",
    "\n",
    "# P(C|-B,A)\n",
    "qcA.x(1)\n",
    "qcA.mcry(theta_C_nBA, [qrA[0], qrB[0]], qrC[0])\n",
    "qcA.x(1)\n",
    "\n",
    "# P(C|B,-A)\n",
    "qcA.x(0)\n",
    "qcA.mcry(theta_C_BnA, [qrA[0], qrB[0]], qrC[0])\n",
    "qcA.x(0)\n",
    "\n",
    "# P(C|B,A)\n",
    "qcA.mcry(theta_C_BA, [qrA[0], qrB[0]], qrC[0])\n",
    "\n",
    "\n",
    "qcA.draw('mpl', style='bw', plot_barriers=False, justify='none', fold=-1)"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Run circuit"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "outputs": [],
   "source": [
    "# Get the Aer simulator backend\n",
    "simulator_backend = Aer.get_backend('aer_simulator')\n",
    "def run_circuit(circuit, shots=10_000):\n",
    "    \"\"\"\n",
    "    Run the provided quantum circuit on the Aer simulator backend.\n",
    "\n",
    "    Parameters:\n",
    "    - circuit: The quantum circuit to be executed.\n",
    "    - shots (default=10,000): The number of times the circuit is executed.\n",
    "\n",
    "    Returns:\n",
    "    - counts: A dictionary with the counts of each quantum state result.\n",
    "    \"\"\"\n",
    "    # Transpile the circuit for the given backend\n",
    "    transpiled_circuit = transpile(circuit, simulator_backend)\n",
    "\n",
    "    # Run the transpiled circuit on the simulator\n",
    "    job = simulator_backend.run(transpiled_circuit, shots=shots)\n",
    "    result = job.result()\n",
    "\n",
    "    # Get the counts of quantum state results\n",
    "    counts = result.get_counts(transpiled_circuit)\n",
    "\n",
    "    return counts"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Rejection Sampling"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "outputs": [],
   "source": [
    "def qrejection_sampling(circ, N, query, evidence, repetitions=5):\n",
    "    \"\"\"\n",
    "    Perform rejection sampling on a quantum circuit.\n",
    "\n",
    "    Parameters:\n",
    "    - circ: The quantum circuit to be executed.\n",
    "    - N: The number of shots for each execution.\n",
    "    - evidence: The knowledge provided for the sampling.\n",
    "    - repetitions (default=5): The number of times the circuit is executed and sampled.\n",
    "    \"\"\"\n",
    "\n",
    "    for _ in range(repetitions):\n",
    "        # Get samples\n",
    "        counts = run_circuit(circ, shots=N)\n",
    "        # Count valid samples and valid samples corresponding to query\n",
    "        nval_samples = 0\n",
    "        nval_qsamples = 0\n",
    "        for c_key, c_val in counts.items():\n",
    "            # Valid samples\n",
    "            valid = True\n",
    "            for e_key, e_val in evidence.items():\n",
    "                if c_key[e_key]!=e_val:\n",
    "                    valid=False\n",
    "            if valid:\n",
    "                nval_samples+=c_val\n",
    "                # Query samples\n",
    "                qsample = True\n",
    "                for q_key, q_val in query.items():\n",
    "                    if c_key[q_key]!=q_val:\n",
    "                        qsample=False\n",
    "                if qsample:\n",
    "                    nval_qsamples+=c_val\n",
    "\n",
    "        # Compute and print the ratio\n",
    "        print(nval_qsamples/nval_samples)\n",
    "        return 0\n",
    "\n"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Run rejection sampling"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.3231552162849873\n",
      "0.3351324828263003\n",
      "0.3456998313659359\n",
      "0.33620159803318994\n",
      "0.3323004279179578\n",
      "0.33134954932707744\n"
     ]
    }
   ],
   "source": [
    "qcCopy = QuantumCircuit(qrA, qrB, qrC, name='Bayes net')\n",
    "qcCopy.append(qcA, [qrA, qrB, qrC])\n",
    "qcCopy.measure_all()\n",
    "\n",
    "query = {1: '1'}\n",
    "evidence = {0: '0', 2: '0'}\n",
    "shots = [1_000, 5_000, 10_000, 20_000, 50_000, 100_000]\n",
    "for N in shots:\n",
    "    qrejection_sampling(qcCopy,N=N, query=query, evidence=evidence)"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Implement $S_e$"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "outputs": [
    {
     "data": {
      "text/plain": "<Figure size 1458.05x284.278 with 1 Axes>",
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABHYAAADuCAYAAABcS4C9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8pXeV/AAAACXBIWXMAAA9hAAAPYQGoP6dpAAAyZklEQVR4nO3deXQUVd7G8acTICFhFyRIQBABIaxhHwVBQZHVVxTQyCI4yiZRljCgDIijDgoOEPYdRTbBURYZRTGAoJhFECGgogjBNBAISFay9PvHHHJkSEK6092V6v5+zvEPqbp1f1V9u9L9dNUti81mswkAAAAAAACm42N0AQAAAAAAAHAMwQ4AAAAAAIBJEewAAAAAAACYFMEOAAAAAACASRHsAAAAAAAAmBTBDgAAAAAAgEkR7AAAAAAAAJgUwQ4AAAAAAIBJEewAAAAAAACYFMEOAAAAAACASRHsAAAAAAAAmBTBDgAAAAAAgEkR7AAAAAAAAJgUwQ4AAAAAAIBJEewAAAAAAACYFMEOAAAAAACASRHsAAAAAAAAmBTBDgAAAAAAgEkR7AAAAAAAAJgUwQ4AAAAAAIBJEewAAAAAAACYFMEOAAAAAACASRHsAAAAAAAAmBTBDgAAAAAAgEkR7AAAAAAAAJgUwQ4AAAAAAIBJEewAAAAAAACYFMEOAAAAAACASRHsAAAAAAAAmBTBDgAAAAAAgEkR7AAAAAAAAJgUwQ4AAAAAAIBJEewAAAAAAACYFMEOAAAAAACASRHsAAAAAAAAmBTBDgAAAAAAgEkR7AAAAAAAAJhUKaMLQP5sNpvS0tKMLqPIAgICZLFYjC7DY5jt9ZecPwa8/Rh4+/6DMeDt+w8AAFBUBDslVFpamsqVK2d0GUWWkpKiwMBAo8vwGGZ7/SXnjwFvPwbevv9gDHj7/gMAABQVt2IBAAAAAACYFFfsmMC5c+dK5C+Aqampql69utFleLyS+vpL7hsD3n4MvH3/wRjw9v0HAAAoDMGOCQQGBpbYD7RwPV5/joG37z8YA96+/wAAAIXhViwAAAAAAACT4oodeLy0tDQdPnxYJ06cUFpamnx8fFShQgU1bdpUjRo1UqlSvA0AT2az2XTmzBnFxcXp/PnzysrKkp+fn2rVqqVWrVqpatWqRpcIAAAAOIxvtPBI58+f14oVK7R+/XodO3ZMOTk5+a5XtmxZtWnTRsOGDVP//v1VtmxZN1cKwBVsNpuioqK0ePFi7d69W0lJSQWuW7t2bfXs2VMjR45U06ZN3VglAAAAUHzcigWPkpiYqCFDhqhWrVqaMmWKjhw5UmCoI0np6enau3evhg4dquDgYE2bNk0ZGRlurBiAM9lsNm3atEmNGzfWAw88oE2bNhUa6kjS6dOntWjRIjVr1kydOnXS119/7aZqAQAAgOIj2IFHsNlsev/99xUSEqJ3331X165ds3sbly5d0owZM9SqVSt9++23LqgSgCudO3dOjz/+uAYMGKDjx487tI19+/bpvvvu08SJE5Wenu7kCgEAAADnI9iB6WVnZ2vYsGF6+umnlZycXOztHTt2TB06dNDixYudUB0Ad4iNjVXTpk314YcfFntbubm5mjVrltq1a6fExEQnVAcAAAC4DsEOTC0nJ0dhYWFavXq1U7ebm5urkSNHau7cuU7dLgDni46OVpcuXXThwgWnbvfIkSPq2LEj4Q4AAABKNIIdmNpLL72kTZs2uWz7L774oku3D6B4Tp06pe7du+vq1asu2f7Jkyf18MMPc1sWAAAASiyCHZjWrl27FBkZaVeb6OhonTlzRtHR0UVuM2LECH6xB0qg3NxcDR8+XJcuXSpyG0fOAUeOHNHUqVMdKREAAABwOYIdmNIff/yhZ5991u52QUFBCg4OVlBQUJHbJCcna8SIEbLZbHb3B8B1lixZot27d9vVxpFzgCS98847OnDggF1tAAAAAHfwimAnKSlJERERuvvuu+Xv769atWopPDxcqampGj58uCwWi+bPn290mbBDZGSkTp8+7bb+tm7dqn379rmtPwCFS09P1yuvvOK2/mw2myIiItzWHwDnstls+vrrrzVq1Cg9+uij6tevnyZMmKBjx44ZXRoAN7ly5Yrmz5+vJ598Un369FFYWJhWrVrF7dbwCB4f7Bw6dEhNmzbV22+/LavVqsaNGysrK0vz5s3TgAEDFB8fL0lq0aKFsYW6SFRUlCwWi2bNmlXgOhaLRb169XJjVcWTnZ2tJUuWuL3fRYsWub3P4vLE199eU6ZMkcVi0cqVK29aZrPZ1LlzZ/n5+emHH34woDrX89QxsHHjRrtuwXKG/fv36/vvv3drn87gqWPAHt5+HvB2e/bsUWhoqP7yl79o0aJF+vjjj/Xhhx9q9uzZCgkJUZcuXQh4AA+WmZmpl156STVr1tQLL7ygDRs2aNu2bVq3bp2GDRummjVrasaMGcrNzTW6VMBhHh3sJCUlqXfv3rJarRo/frwSExMVFxcnq9WqmTNnaseOHYqOjpbFYlGzZs2MLhdFtGPHDp05c8bt/W7ZskVWq9Xt/aJ4pk+friZNmmjcuHFKSEi4YdmcOXO0Z88evfrqq2rSpIlBFcIRCxcu9Kp+UTycB7zXRx99pG7duunQoUMFrhMVFaV7771XMTEx7isMgFtkZGSoR48emjNnjlJTU/NdJzk5WdOmTVNYWJhycnLcXCHgHB4d7IwdO1YJCQkaM2aMZs2apfLly+cti4iIUPPmzZWdna06deqoQoUKBlYKe2zfvt2QfrOysrRr1y5D+objypQpozVr1uTdenndiRMn9PLLL6tdu3aaOHGigRXCXufPn7dr8mNn2rFjhyH9ong4D3inI0eO6Mknn1RWVtYt1718+bJ69uypCxcuuKEyAO4yatSoIs/Ht2HDBk2fPt21BQEu4rHBTnx8vDZu3KiqVavqzTffzHedVq1aSZKaN29+w7//+uuv6tOnj8qXL6/KlStr8ODBunjxostrRtHExsZ6Zd9wXGhoqCZPnqzPPvtMS5cuVU5OjgYPHiybzaY1a9bI19fX6BJhByPfhwkJCTp//rxh/cNxnAe8z+zZs5WRkVHk9c+fP6/ly5e7sCIA7nT69GmtWbPGrjZz587V1atXXVQR4DqljC7AVdavX6/c3FyFhYWpXLly+a5TtmxZSTcGO1evXlWXLl1UpUoVrV+/Xunp6YqIiFCvXr20f/9++fiYMwtLS0tTUlKS0WUUW2ZmpqFzIJg12PGU1784pk6dqq1bt2rChAk6dOiQvv32W82ePVsNGzY0ujS38KQxEBcXZ2j/sbGxeuSRRwytwRGeNAYc5e3nAW9y8eJFbdiwwe52ixcvVkREBEEf4AGWLl1q97w5V69e1dq1azVy5EgXVQW4hjlTiiK4fsldly5dClzn+n32fw52li5dqrNnz+qjjz5Sr1699MQTT2jdunX65ptvtHXrVtcW7ULTpk1TtWrV8v3PTKxWa5EuqXYVdz6Jy5k85fUvjtKlS2vNmjXKyMjQokWLdN999+nFF180uiy38aQxYPT70Oj+HeVJY8BR3n4e8CZ79uxRZmam3e1Onz6d92ANAOa2c+dOt7YDjOSxV+z89ttvkqQ777wz3+XZ2dnav3+/pBuDne3bt+u+++5T7dq18/6tQ4cOuuuuu7Rt2zY9+uijdtfSunVruyfddfas7M8995yeeOKJfJd169at2NuvX7++W65mulWoEx0draCgoAKXX18WFBRU6ATMVqtVbdq0uenfExISFBwcXMRqHWe2119y/hhwxZMJKlasKD8/P2VlZalHjx5OH7POPAaMgYIlJycXuMxZ5wCp4PPApEmT9NprrxWxWsd5+xhw1dNJXHkecNd7ALdW0CSpRdGlSxf5+fk5sRoARnD0oSeffvqpWz7vA/8rKCjI4Yn8PTbYuf4HPT09Pd/lGzduVFJSksqXL6+6devm/fuxY8fy/eAbEhLi8KMwrVarzp4961BbZ6lfv766du3qsu0nJia6bNv2CAoKKtKJuFSpUg6dsHNzcw1/LR3h6tdfKjljoCA2m03PPPOMrl27pkaNGukf//iH+vfvr3r16jmtj5J8DLxlDLj6HCBJV65c0ZUrVxxqayRvGQOFcfV5oKTvP4rG229ZBLzdtWvXTPl5H97NY4OdoKAgJScnKy4uTh06dLhhWWJiYt7TL5o1ayaLxZK3LDk5WZUqVbppe1WqVNGJEyccrsVeubm5pvqAWKNGDbf8Snmr43KrZD4oKEilSpVSdnZ2oesWtKx06dK6/fbbi1ZsMZjt9ZecPwacfQwiIyMVFRWl119/XX379lVoaKiGDRumqKioG84BxeHMY8AYKNiVK1eUkpKS7zJnnQMK21aVKlXy5mhzJW8fA67Yf1efB9z1HsCtZWdn69y5c3a3s1gsCgoK4nUEPEBycrLS0tLsble+fHmemAxDOJIbXOexwU7Xrl0VHx+vmTNnqlu3bmrQoIGk/16mP2jQoLxfY1q0aOHyWhy5nCo1NbXASZ9Lop9++kmBgYFu6athw4b68ccf812W320Tf3bmzBkFBwfLarWqVq1advc9ZMgQLVu2zO529jLb6y85fww48xj89NNPmjx5stq0aaNJkybJ19dX06dP15QpUxQZGamxY8c6rR9nHQPGQME2b95c4C1Frj4HSP89p//5Sk9X8fYx4Oz9d8d5wJ1/C3FrjzzyiP7zn//Y1Wbs2LGaM2eOawoC4FZxcXF5T0EuKl9fX8XHx6tmzZouqgpwDY/9OSIiIkK33Xabzpw5o5CQEDVt2lT169dX27Ztddddd+mBBx6QdPOjzitXrqzLly/ftL1Lly6pSpUq7igdt2DvCdpT+oZjcnNzNXToUOXk5NzwSOOIiAi1bt1akydP1smTJw2uEvYw8n1YuXJl1alTx7D+4RjOA97J3omxS5cuzZNwAA8SGhqqjh072tWmf//+hDowJY8NdoKDg7Vv3z717NlT/v7+OnXqlKpUqaIlS5Zox44deVd8/G+w06hRo3zn0jl27JgaNWrkltpRuMKedOZqnTt3NqxvOGb27Nk6cOCAZsyYccN72NfXV6tXr1Z2draGDRsmm81mYJWwR506dQwLV7p06eK0W/fgPpwHvNPDDz+sadOmFWldi8WilStXqmHDhi6uCoA7bdiw4YaH4hSmWbNmWrRokYsrAlzDY4Md6b8hzfbt23X16lVdvXpVBw8e1HPPPafU1FSdOnVKPj4+atKkyQ1tevXqpa+++irvUeiSdPDgQZ08eVK9e/d29y4gH08++aTKly/v9n47d+6se+65x+39wnHx8fGaOnWq2rdvr/Hjx9+0PCQkRNOnT9fevXsVGRlpQIVwhMVi0fPPP29I3yNGjDCkXziO84B3mzZtmt555x35+/sXuE7lypX1wQcf6Omnn3ZjZQDc4Y477tCBAwfUvn37Qtfr3r27oqKiVLFiRTdVBjiXxeaFP08dPHhQ7du3V8OGDXX8+PEblv3xxx9q2rSpqlatqldffVUZGRmKiIhQtWrV9PXXX7ttMr0/zy2QkpJSIu/ZN7LGF154QfPnz7e73fX5NRISEuyeX2PTpk0FzuvhbGZ4/SXX1untx8Db9/9WLly4oODgYF27ds2udsU5B9SvX1/Hjx/n78D/4D1gjjq9XXJystasWaP33ntPhw4dUm5ursqUKaPFixdr4MCBbpkQHYBxbDabvvnmGy1YsEAbNmxQTk6OfH199eyzz2rkyJE33cUBmI1HX7FTkCNHjki6+TYsSapQoYJ2796tGjVqaODAgXr22Wf1l7/8Rdu3b+cJCSXI+PHj3frhuVmzZnr00Ufd1h+AwlWrVk1jxoxxa59///vf+TsAmFTlypX14osvKjY2VjVq1JD03/PIM888Q6gDeAGLxaIOHTpo7dq1eU8eCgoK0uLFiwl14BG88hNqYcGOJNWrV0/bt29XSkqKLl++rLVr16patWruLBG3UKdOHb311ltu6atUqVJavXq1Spcu7Zb+ABTNa6+9pvr167ulr969eyssLMwtfQEAAAD2INiBaY0YMUJdu3a1q43ValVCQoKsVmuR27z88stq2bKlveUBcLGAgACtWrUq7wlHReHIOeD6xPtMmgwAAICSqJTRBRhh9+7dRpcAJ/Dx8dGmTZvUuXNnff/990Vq06ZNG7v6GDRokP7+9787Uh4AN7j33nu1cuVKDRkypEjr23sOCAwM1CeffJJ36wYAAABQ0njlFTvwHJUrV9YXX3yhVq1aOX3bzzzzjFauXMmcGkAJN3jwYK1Zs8auK3eKolKlStq1a5fatWvn1O0CAAAAzsQ3Vphe1apVtWfPHqdNpBoYGKgFCxZo+fLlKlXKKy9qA0xn8ODBioqK0t133+2U7XXs2FExMTHq0KGDU7YHAAAAuArBDjxCYGCgIiMj9eWXX6phw4YOb+ehhx7SkSNHNGrUKK7UAUzmvvvu0+HDh/XSSy+pTJkyDm2jYsWKmjt3rqKiolSvXj0nVwgAAAA4H99c4VE6d+6sY8eO6dNPP1Xfvn2LFM6UK1dOI0eO1Pfff69PP/1UdevWdUOlAFwhICBA77zzjs6cOaM333xTd955Z5HatWzZUsuWLdPZs2c1duxYgl0AAACYBveZwOP4+PjooYce0kMPPaSrV6/q0KFDiomJ0fHjx/Xee+8pPT1dgYGBmjNnjlq1aqWQkBCHf90HUDLdfvvt+tvf/qZJkybpl19+UWxsrOLi4hQZGam0tDQFBAQoIiJCrVq1UqtWrZgcGQAAAKZFsAOPVr58eXXs2FEdO3aUJO3YsUNnz55VpUqV9OyzzxpcHQBXs1gsqlevnurVq6f+/ftr7dq1SktLU+XKlTVt2jSjywMAAACKjWvNAQAAAAAATIpgBwAAAAAAwKS4FcsEUlNTjS4hXyW1Lk9Tko+zu2rz9mPg7fuPkn2ceQ+U3NoAAIB3INgxgerVqxtdAgzE688x8Pb9B2PA2/cfAACgMNyKBQAAAAAAYFJcsVNCBQQEKCUlxegyiiwgIMDoEjyK2V5/yfljwNuPgbfvPxgD3r7/AAAARUWwU0JZLBYFBgYaXQYMwuvPMfD2/QdjwNv3HwAAoKi4FQsAAAAAAMCkCHYAAAAAAABMimAHAAAAAADApAh2AAAAAAAATIpgBwAAAAAAwKQIdgAAAAAAAEyKYAcAAAAAAMCkCHYAAAAAAABMimAHAAAAAADApAh2AAAAAAAATIpgBwAAAAAAwKQIdgAAAAAAAEyKYAcAAAAAAMCkCHYAAAAAAABMimAHAAAAAADApAh2AAAAAAAATIpgBwAAAAAAwKQIdgAAAAAAAEyKYAcAAAAAAMCkCHYAAAAAAABMimAHAAAAAADApAh2AAAAAAAATIpgBwAAAAAAwKQIdgAAAAAAAEyKYAcAAAAAAMCkCHYAAAAAAABMimAHAAAAAADApAh2AAAAAAAATIpgBwAAAAAAwKQIdgAAAAAAAEyKYAcAAAAAAMCkCHYAAAAAAABMyiuCnaSkJEVEROjuu++Wv7+/atWqpfDwcKWmpmr48OGyWCyaP3++0WXChXJzc/Xtt98qPT1dkpSVlWVwRQDc7fz588rIyJAkZWZmKjU11eCKALiTzWZTTEwMnwUAL3bx4sW8zwIZGRm6evWqwRUBzuHxwc6hQ4fUtGlTvf3227JarWrcuLGysrI0b948DRgwQPHx8ZKkFi1aGFsoXOKPP/7Q7Nmz1bBhQ7Vr106XLl2S9N8veB07dtT69euVm5trcJUAXOngwYMKCwtTcHCwLl68KOm/gX/NmjUVHh6ukydPGlwhAFdKTU3VvHnz1LhxY7Vp0+aGzwLt27fXu+++q+zsbIOrBOBKcXFxGjp0qGrWrJn3WeDixYu64447NHr0aJ04ccLgCoHisdhsNpvRRbhKUlKSWrZsqYSEBI0fP17Tpk1T+fLlJUlvvfWWJk2apFKlSiknJ0eXL19WhQoVDK4YznTmzBk98sgjOnr0aKHr9evXT2vXrpW/v7+bKgPgLpGRkQoPD1dhf+rKlSunLVu26KGHHnJjZQDcwWq1qmfPnoqLiyt0vZ49e2rTpk0KCAhwU2UA3GX58uUaMWKEcnJyClynbNmy2rBhg/r06ePGygDn8ehg56mnntL69es1ZswYRUZG3rS8RYsWOnz4sOrWratffvnFgArhKsnJybr33nvzrsi6lQEDBmjdunXy8fH4i9gAr7Fq1SoNGzasSOv6+/vryy+/VPv27V1cFQB3SUlJUadOnfTdd98Vaf0+ffroww8/lK+vr4srA+AuGzdu1MCBA4u0bunSpfXpp5+qS5cuLq4KcD6P/RYbHx+vjRs3qmrVqnrzzTfzXadVq1aSpObNm+f9W0JCgsaMGaO2bdvKz89PFovFLfXCud55550ihzrSf0/6X3zxhQsrAuBOKSkpCg8PL/L6GRkZeuGFF1xYEQB3W7BgQZFDHUnaunWrtm3b5sKKALhTRkaGxowZU+T1s7KyNHr06EKv8gVKKo8Ndq7PnRIWFqZy5crlu07ZsmUl3Rjs/Pzzz9qyZYuCgoLUpk0bt9QK57p27ZqWLVtmd7uFCxe6oBoARnj//fftnhAxJiZG0dHRLqoIgDvl5ORo8eLFdrfjswDgOTZv3qykpCS72sTHx2vPnj0uqghwHY8Ndnbv3i1JhV5Kl5CQIOnGYKdTp05KTEzU1q1b1bVrV9cWCZf47LPPdO7cObvbbd26NW9CRQDmtnr1aofarVmzxrmFADDE3r17derUKbvb7dq1K+/zIQBzc/RvuqOfIQAjlTK6AFf57bffJEl33nlnvsuzs7O1f/9+STcGO66YY6V169ayWq1O3y7yl5KS4lC73NxcNWrUSKVLl3ZyRQDcLTEx0aF2K1as0EcffeTcYgC4XVpamsNtQ0NDVaZMGSdWA8AIjvzQK0kbNmzQ559/7uRqgFsLCgpSTEyMQ209NthJTU2VJKWnp+e7fOPGjUpKSlL58uVVt25dl9ZitVp19uxZl/YB5zh//rzRJQAwUEZGBudrwMtduHDB6BIAGCgzM5PPAjAdjw12goKClJycrLi4OHXo0OGGZYmJiZo4caIkqVmzZi6fIDkoKMil28eNMjIydPHiRYfa1qhRgydjAR7gwoULunbtmt3tAgMDValSJecXBMCtMjMz7Z5b47qgoCCejAV4gKSkJGVmZtrdLiAgQJUrV3ZBRUDhipMbeGyw07VrV8XHx2vmzJnq1q2bGjRoIEmKjo7WoEGD8v7Yt2jRwuW1OHo5FRyTlZWl2rVr233722OPPaYtW7a4qCoA7rRs2TI999xzdrfbs2dP3hMTAZhXbm6u7r77bv366692tXv44Yf1n//8x0VVAXCndevWKSwszO52O3fuVKdOnVxQEeA6HntpQkREhG677TadOXNGISEhatq0qerXr6+2bdvqrrvu0gMPPCDpxvl14BlKly7t0Be6UaNGuaAaAEZ46qmnVKFCBbvatG3bllAH8BA+Pj4aOXKk3e34LAB4jn79+qlatWp2tQkJCVHHjh1dVBHgOh4b7AQHB2vfvn3q2bOn/P39derUKVWpUkVLlizRjh079OOPP0oi2PFU48aNU5MmTYq8flhYWF7YB8D8AgMDNX/+/CKvHxAQYNf6AEq+UaNGqXXr1kVe/7HHHlPPnj1dWBEAd/Lz89OiRYuKPO1GmTJl7FofKEk8NtiRpEaNGmn79u26evWqrl69qoMHD+q5555TamqqTp06JR8fH7u+/MM8KlasqE8//VTNmjW75boDBw7UihUrOIkDHmbQoEFatGjRLefNqlChgrZt26Y2bdq4qTIA7hAYGKhPPvmkSO/tRx99VGvXrmVuHcDD9OvXT6tWrVKpUoXPQBIYGKh///vfXK0D0/LoYKcgR48elc1mU/369RUQEHDT8s2bN2vz5s06duzYDf/PXDnmcscdd2j//v2aN2+e7rnnnpuWd+nSRR988IHef/99+fn5GVAhAFcbMWKEvv32Ww0ZMuSm93mVKlU0YcIEHT58mCv2AA9VrVo17dmzR4sWLcr3x7yOHTtq3bp12rx5s8qWLWtAhQBcbciQIYqNjdWzzz570/u8YsWKCg8P16FDh9SjRw+DKgSKz2Kz2WxGF+Fuy5cv11//+lf1799fGzduvGl5QVduDBkyRKtXr3ZxdXAFm82m7777Tg899JAuXryo6tWr2z25MgBzu3jxou655x4lJSWpatWqOn36NF/kAC9is9n0/fff68EHH9TFixd1++2369y5c0aXBcCNkpOT1aBBg7zPAqdOnVJgYKDRZQHF5rFPxSrMkSNHJBU8v44XZl0ez2KxKDQ0VP7+/pJ0y8sxAXie2267Le+qHT8/P0IdwMtYLBY1b94877NA6dKlDa4IgLtVrlz5hs8ChDrwFF55K9atgh0AAAAAAAAz8MrLFnbv3m10CQAAAAAAAMXmlVfsAAAAAAAAeAKCHQAAAAAAAJMi2AEAAAAAADApgh0AAAAAAACTItgBAAAAAAAwKYIdAAAAAAAAkyLYAQAAAAAAMCmCHQAAAAAAAJMi2AEAAAAAADApgh0AAAAAAACTItgBAAAAAAAwKYIdAAAAAAAAkyLYAQAAAAAAMCmCHQAAAAAAAJMi2AEAAAAAADApgh0AAAAAAACTItgBAAAAAAAwKYIdAAAAAAAAkyLYAQAAAAAAMCmCHQAAAAAAAJMi2AEAAAAAADApgh0AAAAAAACTItgBAAAAAAAwKYIdAAAAAAAAkyLYAQAAAAAAMCmCHQAAAAAAAJMi2AEAAAAAADApgh0AAAAAAACTItgBAAAAAAAwKYIdAAAAAAAAkypldAEAAAC4mc1mU1pamtFlFFlAQIAsFovRZXgMs73+kvPHgLcfA2/ffzAGvH3/7UGwAwAAUAKlpaWpXLlyRpdRZCkpKQoMDDS6DI9httdfcv4Y8PZj4O37D8aAt++/PbgVCwAAAAAAwKS4YgcAAKCEO3fuXIn8FTw1NVXVq1c3ugyPV1Jff8l9Y8Dbj4G37z8YA96+/7dCsAMAAFDCBQYGltgPtHA9Xn+OgbfvPxgD3r7/t8KtWAAAAAAAACZFsAMAAAAAAGBS3IoFAAAAj3XlyhXFxcUpNjZWx48f16VLlyRJycnJmj9/vlq1aqXmzZsrICDA4EoBuEJubq5+/vlnxcbGKi4uLu8ccOnSJU2ZMkWhoaFq1aqV6tSpw6PaYVoEOwAAAPAoubm52rVrlxYuXKjt27crNzf3pnXS0tL0wgsvSJL8/Pw0YMAAjRo1Sm3btuXLHeABzp8/rxUrVmjx4sU6ffr0TcvT09P15ptv5v1/w4YNNWrUKA0ePFiVKlVyY6VA8XErFgAAADzGvn37FBISou7du2vr1q35hjr/KzMzU++++67at2+vjh07Kj4+3g2VAnCF9PR0TZgwQcHBwZoyZUq+oU5+Tpw4ofDwcNWsWVNvvPGGsrOzXVwp4DwEOwAAADC9tLQ0hYeH6/7779fx48cd3s7+/fvVsmVLvfXWW8rJyXFihQBc7cCBA2rRooVmz56trKwsh7aRlpaml19+We3bt9cPP/zg5AoB1yDYAQAAgKldunRJDzzwgObNmyebzVbs7WVmZmrSpEnq37+/MjMznVAhAFdbv3697r//fv34449O2V5sbKzatWunzz//3CnbA1yJYAcAAACmdeXKFXXr1k0HDx50+rY//PBD9e/fn1sygBJuw4YNCgsLc/p7NS0tTb169dLu3budul3A2Qh2AAAAYEo2m01hYWGKi4tzWR9bt25VRESEy7YPoHhiYmI0aNAgp1ytl5/MzEz17dtXv/zyi0u2DzgDwQ4AAABMafXq1dqxY4ddbaKjo3XmzBlFR0cXuc2cOXO0b98+e8sD4GKZmZkaOnSoXVfqOHIOSElJ0fDhw4s0GTtgBIIdAAAAmM7vv/+ul156ye52QUFBCg4OVlBQUJHb2Gw2DRs2TOnp6Xb3B8B1ZsyYoaNHj9rVxpFzgCRFRUVp8eLFdrUB3MUrgp2kpCRFRETo7rvvlr+/v2rVqqXw8HClpqZq+PDhslgsmj9/vtFlwkXOnTun119/XSEhIfr9998lSRcuXNDKlSuVlpZmcHUAXM1ms+nzzz/XY489psTEREmS1WrV888/r8OHDxtcHeB6UVFRslgsmjVrVoHrWCwW9erVy41VFd/cuXN15coVt/X3888/a/369W7rz1k89fW3x5QpU2SxWLRy5cqbltlsNnXu3Fl+fn4e+wQkTx0Dly9f1pw5c9za5+uvv+7w07aM5KljwB6efh7w+GDn0KFDatq0qd5++21ZrVY1btxYWVlZmjdvngYMGKD4+HhJUosWLYwtFE5ns9n06quvqlatWnrllVd07NixvHtvr127puHDhys4OFhbt241uFIArnLy5Ek1b95c3bp107///e+8S6hzcnK0dOlStWjRQr169XLrl0MAxZeRkaEVK1a4vd8FCxa4bB4PuM706dPVpEkTjRs3TgkJCTcsmzNnjvbs2aNXX31VTZo0MahCOOLdd991+4+0v//+u7Zt2+bWPuEcnn4e8OhgJykpSb1795bVatX48eOVmJiouLg4Wa1WzZw5Uzt27FB0dLQsFouaNWtmdLlwIpvNppdeeknTp08vNFVPTk7W//3f/+mDDz5wY3UA3OGXX37RvffeqyNHjhS63o4dO/Tggw8qJSXFTZUBKK4tW7bo4sWLbu83Li5OMTExbu8XxVOmTBmtWbMm72r9606cOKGXX35Z7dq108SJEw2sEI5YsmSJIf0uWrTIkH5RPJ5+HvDoYGfs2LFKSEjQmDFjNGvWLJUvXz5vWUREhJo3b67s7GzVqVNHFSpUMLBSONvmzZs1d+7cIq2bm5urQYMG6dSpU64tCoDb2Gw29evXT+fOnSvS+rGxsQoPD3dxVQCc5csvvzSs76ioKMP6huNCQ0M1efJkffbZZ1q6dKlycnI0ePBg2Ww2rVmzRr6+vkaXCDucP39ex44dM6Tvr776yumPVYd7ePJ5oJTRBbhKfHy8Nm7cqKpVq+rNN9/Md51WrVrp8OHDat68ed6/bd68WevXr1dMTIwuXLig2rVrq1+/fpo8ebLKlSvnrvJRTEUNda7LzMzU0qVL9cYbb7ioIgDutHfvXh06dMiuNmvXrtXMmTNVtWpV1xQFGCwtLU1JSUlGl+EUsbGxXtl3cXjS6++oqVOnauvWrZowYYIOHTqkb7/9VrNnz1bDhg2NLs0tPGkMGPk+zMjI0LFjx0x5x4cnjQFHeep5wGOv2Fm/fr1yc3MVFhZWYCBTtmxZSboh2Jk1a5Z8fX31xhtvaOfOnRo5cqQWLVqk7t2783g7k/j++++1f/9+u9stX75cmZmZLqgIgLstXLjQ7jbXrl3Ld0I9wFNMmzZN1apVy/c/M8nMzDR0csu4uDjD+i4OT3n9i6N06dJas2aNMjIytGjRIt1333168cUXjS7LbTxpDBj9PjS6f0d50hhwlKeeBzz2ip3du3dLkrp06VLgOtcnTfpzsLNt27YbBvb999+vatWqKSwsTF999ZU6depkdy2tW7eW1Wq1ux0ck5qa6lC7Cxcu6M4771SpUh77tgC8hqPn3OnTp2vevHlOrgZwjLN/UHruuef0xBNP5LusW7duxd5+/fr15ePj+t8Mc3JyCr0NIjo6utDHGF9fFhQUpDNnzhS4ntVqVZs2bW76959//lnBwcF2VOwYs73+kvPHgCt+VK1YsaL8/PyUlZWlHj16OH3MOvMYMAYKdvny5QKXOescIBV8HggPD9crr7xStGKLwdvHgKsurHDleaA4+x8UFOTwPG4e+w32t99+kyTdeeed+S7Pzs7Ou6rjz8FOfmll69atJUlnz551qBar1epwW7hXUefjAOCZ0tPTOV/DY9WvX19du3Z12fYTExNdtm17BAUFFSl4KVWqlEMBjc1mM+V5wtWvv1RyxkBBbDabnnnmGV27dk2NGjXSP/7xD/Xv31/16tVzWh8l+Rh4yxhw9TlAkv744w/98ccfDrU1kreMgcK4+jxg1P57bLBz/aqN9PT0fJdv3LhRSUlJKl++vOrWrVvotq5P0NeoUSOHaiksMYbzpaWlKTk52aG21atX54odwAOcO3fOoYkNy5YtqypVqrigIsB+ubm5Jf4D8p/VqFHDbVfsFHZV3q2u2AsKClKpUqWUnZ3t0HZ8fHxUo0aNohVbDGZ7/SXnjwFnH4PIyEhFRUXp9ddfV9++fRUaGqphw4YpKipKFovFKX048xgwBgp2+fLlAq/Sd9Y5oLBtVahQ4YYH87iKt48BV+y/q88Dxdn/4uQGHvsNNigoSMnJyYqLi1OHDh1uWJaYmJj3KLNmzZoV+gKePXtWU6dOVffu3dWiRQuHauGxmO6VmJio2rVr2/2lrnHjxvrhhx+c9ocdgHHGjRunf/3rX3a3W7ZsmcLCwlxQEWC/1NRUUz244aefflJgYKDL+8nJyVGlSpWUkpKS7/L8bpv4szNnzig4OFhWq1W1atWyu//mzZu7ZX4Ns73+kvPHgDOPwU8//aTJkyerTZs2mjRpknx9fTV9+nRNmTJFkZGRGjt2rNP6cdYxYAwUbO7cuQXOi+Lqc4AkrVixQo8//rhDbe3h7WPA2fvvjvOAu94D/8tjJ0++fonZzJkz9eOPP+b9e3R0tLp06ZI3G3hhYU1KSor69u2rMmXKMKGmidSoUUOPPfaY3e1GjRpFqAN4iBEjRtjdpmrVqm75kAageHx9fdWyZUvD+m/VqpVhfcMxubm5Gjp0qHJycm54pHFERIRat26tyZMn6+TJkwZXCXsY/T40un/Yz9PPAx4b7EREROi2227TmTNnFBISoqZNm6p+/fpq27at7rrrLj3wwAOSbpxf58/S09PVu3dv/frrr/rss8/ccsktnGfixIkqU6ZMkdevXbu2Bg0a5MKKALhTgwYNNHDgQLvaREREyM/Pz0UVAXAmI79U8YXOfGbPnq0DBw5oxowZN0yt4Ovrq9WrVys7O1vDhg2TzWYzsErYo0WLFob9IFu5cmXVqVPHkL7hOE8/D3hssBMcHKx9+/apZ8+e8vf316lTp1SlShUtWbJEO3bsyLuKJ79gJysrS48//rhiYmK0c+dONW7c2N3lo5hat26ttWvXFmm+nNtvv107d+5UhQoV3FAZAHdZsWJFkZ9k+Pzzz2vChAkurgiAs/Tu3duQfn18fPTII48Y0jccEx8fr6lTp6p9+/YaP378TctDQkI0ffp07d27V5GRkQZUCEeUK1cu74d6d+vVqxdX+ZuMN5wHLDazRlLFkJKSogoVKshisejq1asKCAjIW5abm6uBAwdq69at+uSTTww7YcA5vvjiC02YMEGHDh26aZnFYlH37t21YMGCW06gDcCcMjIyNGnSJK1YsSLfSRarV6+uiRMnaty4cXxIQ4nz57kFUlJSDLln/1aMqtFms6lRo0Y6ceKE3W2vz6+RkJBg9/waffr00ccff2x3n44ww+svubZObz8G3r7/t7JlyxaHbqEuzjlAkr755hu1a9fO7naO8PYx4O37bw+PnTy5MEePHpXNZlODBg1uCHUkafTo0frggw/0t7/9TQEBAfrmm2/yltWrVy/fx6Gj5HrwwQcVFxenb775RuvWrZPValWpUqVUr149DRs2THfddZfRJQJwIX9/f82dO1evvfaa3nvvPR08eFCpqamqWLGiunfvrkcffdSu2zYBlAwWi0WjRo1SeHi4W/sdOXKkW/sDULA+ffrojjvu0O+//+62PkNDQ9W2bVu39QcUlVcGO0eOHJGU/21YO3fulCT985//1D//+c8blq1atUpDhw51eX1wLovFog4dOtz0dDQA3qNChQoaPXq0Ro8ebXQpAJzkr3/9qxYuXOjQVTuOePjhh/Xwww+7pS8At1a6dGm99dZbevrpp93W5+zZs7nCFyWSx86xU5jCgp1Tp07JZrPl+x+hDgAAQMlQtmxZrV69Wj4+rv84W758eS1btowvdEAJ89RTT6lv375u6Wv06NHq3LmzW/oC7EWwAwAAAFNq3769pkyZYlcbq9WqhIQEWa3WIrdZsGCBQ3NxAHAti8WixYsX2/UEY0fOAffcc89Nd3MAJYlX3oq1e/duo0sAAACAE8yYMUNWq1XLly8v0vpt2rSxa/szZ87UoEGDHCkNgBsEBQVp165d6ty5s5KSkm65vr3ngDp16mjXrl15k+MCJZFXXrEDAAAAz2CxWLRkyRK9+OKLTt2uj4+P5s2bp4iICKduF4DzhYSEaO/evU6/sq5x48bat2+fgoODnbpdwNkIdgAAAGBqPj4++te//qWPP/5YQUFBxd7ePffcowMHDuiFF15wQnUA3KFRo0Y6fPiwU+ZFtVgsGjdunGJiYgh1YAoEOwAAAPAIffr00dGjRzVy5EgFBATY3f62227T3//+d3333Xdq166dCyoE4EqVK1fWqlWrtGPHDoWGhjq0jU6dOmnfvn2aPXu2ypYt6+QKAdcg2AEAAIDHqFKlihYuXKizZ89q7ty5atu2rcqUKVPg+oGBgercubPWrFmjhIQEvfrqq/L393djxQCcrUePHoqJidHBgwc1ZMiQW151U7duXY0cOVJHjhzRnj17dO+997qpUsA5vHLyZAAAAHi2SpUqaezYsRo7dqyuXbumo0eP6tixY0pLS5OPj4/KlSunZs2aqUGDBvL19TW6XABOZrFY1LZtW7Vt21aSdO7cOcXGxur8+fO6du2a/Pz8VLNmTYWGhqpKlSoGVwsUD8EOAAAAPFqZMmXUsmVLtWzZ0uhSABikevXq6tGjh9FlAC7BrVgAAAAAAAAmxRU7AAAAJVxqaqrRJeSrpNblaUrycXZXbd5+DLx9/1GyjzPvAeNrI9gBAAAo4apXr250CTAQrz/HwNv3H4wBb9//W+FWLAAAAAAAAJOy2Gw2m9FFAAAA4EY2m01paWlGl1FkAQEBslgsRpfhMcz2+kvOHwPefgy8ff/BGPD2/bcHwQ4AAAAAAIBJcSsWAAAAAACASRHsAAAAAAAAmBTBDgAAAAAAgEkR7AAAAAAAAJgUwQ4AAAAAAIBJEewAAAAAAACYFMEOAAAAAACASRHsAAAAAAAAmBTBDgAAAAAAgEkR7AAAAAAAAJgUwQ4AAAAAAIBJEewAAAAAAACYFMEOAAAAAACASRHsAAAAAAAAmBTBDgAAAAAAgEkR7AAAAAAAAJgUwQ4AAAAAAIBJEewAAAAAAACYFMEOAAAAAACASRHsAAAAAAAAmBTBDgAAAAAAgEkR7AAAAAAAAJgUwQ4AAAAAAIBJEewAAAAAAACYFMEOAAAAAACASf0/2QVToXJvCegAAAAASUVORK5CYII="
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "def S_e(observation):\n",
    "    opSe = QuantumCircuit(3)\n",
    "    for obs in observation:\n",
    "        crtls = list({0,1,2} - {obs})\n",
    "        # multi control z gate\n",
    "        opSe.h(obs)\n",
    "        opSe.mcx(crtls, obs)\n",
    "        opSe.h(obs)\n",
    "        # x gate\n",
    "        opSe.x(obs)\n",
    "        # multi control z gate\n",
    "        opSe.h(obs)\n",
    "        opSe.mcx(crtls, obs)\n",
    "        opSe.h(obs)\n",
    "        # c gate\n",
    "        opSe.x(obs)\n",
    "    return opSe\n",
    "\n",
    "obs = {0,2}\n",
    "opSe = S_e(obs)\n",
    "opSe.draw('mpl', style='bw', plot_barriers=False, justify='none', fold=-1)"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Implement Grover operator"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "outputs": [],
   "source": [
    "opG = GroverOperator(opSe, qcA)"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Build circuit"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "outputs": [],
   "source": [
    "\n",
    "qcircG = QuantumCircuit(qrA, qrB, qrC, name='Grover Bayes net')\n",
    "qcircG.append(qcA, [qrA, qrB, qrC])\n",
    "qcircG.append(opG, [qrA, qrB, qrC])\n",
    "\n",
    "qcircG.measure_all()"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [
    "# Rejection Sampling"
   ],
   "metadata": {
    "collapsed": false
   }
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.3409090909090909\n",
      "0.36637931034482757\n",
      "0.317351598173516\n",
      "0.3168433451118963\n",
      "0.3361462728551336\n",
      "0.3258426966292135\n"
     ]
    }
   ],
   "source": [
    "query = {1: '1'}\n",
    "evidence = {0: '0', 2: '0'}\n",
    "shots = [1_000, 5_000, 10_000, 20_000, 50_000, 100_000]\n",
    "for N in shots:\n",
    "    qrejection_sampling(qcircG, N, query=query,evidence=evidence)"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "source": [],
   "metadata": {
    "collapsed": false
   }
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}