{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 114,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import modules\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 115,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Parameters\n",
    "N = 64\n",
    "steps = 60\n",
    "loops = 3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 116,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0.01538644 0.01538644 0.01538644 ... 0.01538644 0.01538644 0.01538644]\n",
      " [0.01538644 0.01538644 0.01538644 ... 0.01538644 0.01538644 0.01538644]\n",
      " [0.01538644 0.01538644 0.01538644 ... 0.01538644 0.01538644 0.01538644]\n",
      " ...\n",
      " [0.01538644 0.01538644 0.01538644 ... 0.01538644 0.01538644 0.01538644]\n",
      " [0.01538644 0.01538644 0.01538644 ... 0.01538644 0.01538644 0.01538644]\n",
      " [0.01538644 0.01538644 0.01538644 ... 0.01538644 0.01538644 0.01538644]]\n"
     ]
    }
   ],
   "source": [
    "# Initial state\n",
    "psi = np.zeros((N,N + loops - 1))\n",
    "for i in range(N):\n",
    "    for j in range(N + loops - 1):\n",
    "            psi[i][j] = 1/np.sqrt(N*(N + loops - 1))\n",
    "print(psi)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 117,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Specify the marked vertex\n",
    "a = 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 118,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create an array to store the probabilities.\n",
    "prob = np.zeros(steps+1)\n",
    "\n",
    "# Print the probability at the marked vertex\n",
    "prob[0] = np.sum(np.square(psi[a][:]));\n",
    "#print(0,\"\\t\",prob[0])\n",
    "\n",
    "for t in range(steps):\n",
    "    # Query the oracle.\n",
    "    for i in range(N + loops - 1):\n",
    "        psi[a][i] *= -1\n",
    "    \n",
    "    # Temporary variable.\n",
    "    tmp = np.zeros((N,N + loops - 1))\n",
    "\n",
    "    # Apply the Grover coin\n",
    "    for i in range(N):\n",
    "        # Calculate the average amplitude at this vertex.\n",
    "        avg = 0\n",
    "        for j in range(N + loops - 1):\n",
    "            avg += psi[i][j]\n",
    "        avg /= N + loops - 1\n",
    "    \n",
    "        # Invert each amplitude about this average.\n",
    "        for j in range(N + loops - 1):\n",
    "            tmp[i][j] = 2*avg - psi[i][j]\n",
    "    \n",
    "    # Flip-Flop Shift\n",
    "    for i in range(N):\n",
    "        for j in range(N):\n",
    "            psi[i][j] = tmp[j][i]\n",
    "        for j in range(N, N + loops - 1):\n",
    "            psi[i][j] = tmp[i][j]\n",
    "            \n",
    "    # Print the probability at the marked vertex\n",
    "    prob[t+1] =  np.sum(np.square(psi[a][:]))\n",
    "    #print(t+1,\"\\t\",prob[t+1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 119,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.0, 1.0)"
      ]
     },
     "execution_count": 119,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAAA22klEQVR4nO29eXhb533n+3kBkOAOipu4aaE227IW25El70ucNF4mVaadpM7SJh6njpM4Tnvn3tt05rb3mcn0ee5MOjNNJ17iOB43abYmzSRO4sbZLFuWV8mWbO2iKEqiAG6iAHABCAJ47x/AoWiKC0gCONvv8zx6JBwc4vzOEfh93/e3vUprjSAIgmB/PGYbIAiCIOQHEXRBEASHIIIuCILgEETQBUEQHIIIuiAIgkMQQRcEQXAI8wq6UuoppVS/UurgLO8rpdTfK6U6lVJvK6Wuyb+ZgiAIwnzkMkN/GrhzjvfvAtZn/zwAPLZ0swRBEISFMq+ga61fBIbmOGUn8C2d4VWgVinVki8DBUEQhNzw5eEz2oCzU173ZI+Fpp+olHqAzCyeysrK91x++eV5uLwgCIJ72Ldv36DWunGm9/Ih6GqGYzP2E9BaPwE8AbBt2za9d+/ePFxeEATBPSilTs/2Xj6yXHqAFVNetwPBPHyuIAiCsADyIejPAH+SzXa5DohorS9xtwiCIAiFZV6Xi1Lqe8BtQINSqgf4f4ESAK3148CzwN1AJzAG3FcoYwVBEITZmVfQtdYfned9DXw+bxYJgiAIi0IqRQVBEByCCLogCIJDEEEXBEFwCCLogiAIDkEEXRAEwSGIoAuCIDgEEXRBEASHIIIuCILgEETQBUEQHIIIuiAIgkMQQRcEQXAIIuiCIAgOQQRdEATBIYigC4IgOAQRdEEQBIcggi4IguAQRNAFQRAcggi6IAiCQxBBFwRBcAgi6IIgCA5BBF0QBMEhiKALgiA4BBF0QRAEhyCCLgiC4BBE0AVBEByCCLogCIJDEEEXBEFwCCLogiAIDkEEXRAEwSGIoAuCIDgEEXRBEASHIIIuCILgEETQBUEQHIIIuiAIgkMQQRcEQXAIOQm6UupOpdQxpVSnUupLM7wfUEr9TCl1QCl1SCl1X/5NFQRBEOZiXkFXSnmBR4C7gI3AR5VSG6ed9nngsNZ6K3Ab8N+UUqV5tlUQBEGYg1xm6NuBTq11l9Y6AXwf2DntHA1UK6UUUAUMAcm8WioIgiDMSS6C3gacnfK6J3tsKl8DrgCCwDvAF7XW6ekfpJR6QCm1Vym1d2BgYJEmC4IgCDORi6CrGY7paa8/AOwHWoGrgK8ppWou+SGtn9Bab9Nab2tsbFygqYIgCMJc5CLoPcCKKa/byczEp3If8GOdoRM4BVyeHxMFQRCEXMhF0N8A1iulOrKBznuBZ6adcwa4A0AptRy4DOjKp6GCIAjC3PjmO0FrnVRKPQQ8B3iBp7TWh5RSD2bffxz4MvC0UuodMi6av9BaDxbQbkEQBGEa8wo6gNb6WeDZaccen/LvIPB7+TVNEARBWAhSKSoIguAQRNAFQRAcggi6YHtiiRRaT8+kFQpB/3CckXGpGbQqOfnQhYXxVz85yE/2n8Pv8+D3eSn1efD7PNxxRRP/1wckmzOf/GhfD//nDw9Q6vVQX1VKfVUpDVV+Ohoq+fd3X0GJV+Ys+eLs0Bi3/+0ukmlNY7WfjvpKOhoqWdNYyR9du4LaCun2YTYi6HlmIpXmx2/2sLapik1tARLJNOPJNEdDUZ7cfYqH71iP3+c120zH8Ow7IZbX+PnQ1W2cH0kwODLOuQsxdh0b4ANXNnPdmnqzTXQMvz7cRzKt+cJ719EXjXNqcJTfHu3nB3vHGUuk+PP3bzDbRNcjgp5n9p8NM5pI8bnb1nLnppbJ47861MsD397H/jNhdojI5IXxZIpXTp7nw9va+cu7rpg8Hhmb4Kov/4rXuoZE0PPI88f6WdtYyb/7vcvedfyev9/NG91DJlklTEXWo3nmpRODKAXXr2l41/EdHfUoBS+fPG+SZc5jX/cFYhMpbln/7jYSgYoSLm+u4fVuedb5YnQ8yWtdQ7z38qZL3rt2dR1vnQkzkbqkfZNQZETQ88yezkG2tAUIVJS863igooRNrQFe6RKRyRcvHB+gxKu4fu2ls/AdHXXsO32BRFJEJh/s6RwkkUpz+2UzC3psIsWhYNQEy4SpiKDnkeH4BG+dDXPT+oYZ379+bT1vnblALJEqsmXO5IXjA2xbVUel/1LP4XVr6ohPpHnnXLj4hjmQ54/1U+X3sW113SXvXbt6GQB7xe1iOiLoeeT1U0Ok0pob180i6GvqmUhp9p2+UGTLnEd/NM7R3mFu2TBz187tHZlZ+6tdIjJLRWvN80cHuHl9A6W+SyWjqaaMVfUVvH5KnrXZiKDnkd0nBikr8XDNymUzvn9tRx1ej+KVLmlzs1RePJF5hrdsmHnwrKssZcPyKl4TkVkyR0LD9EbjM7pbDLatqmPv6QtSD2AyIuh5ZE/nINeurqOsZOa0xCq/jy3tAQmM5oEXjg/QUOXniuZL2u5Psr2jjn3dQyQlWLcknj/WD8Btl8++h8H2jmUMjSY4OTBaLLOEGRBBzxN90Tgn+ke4aRZ3i8H1a+p5uyci1XZLIJXWvHRigFs2NODxzLT/SoYdHfWMJiRYt1SeP9rP5rYATdVls55j+NbFj24uIuh5Yk9nxgUwm//c4Ia1DaTSWvJ2l8DBcxEujE1w6yz+c4MdazIi89opWREtlgujCd48c4HbL5v7Wa9pqKS+spTX5XttKiLoeeKlzkHqKkvZ2DK7CwDgPauWUeJVvCJul0XzwvEBlGLe1VBTdRlrGip5TQKji+bFEwOkNdw+Q/75VJRSbFu9jL3dEvA3ExH0PKC15qUTg9ywtn5OFwBAeamXq1cuE0FfAi8eH2BzW4D6Kv+85+5YU8fr3ZnsI2Hh/O5oP/WVpWxtr5333GtX13FmaIy+aLzwhgkzIoKeBzr7R+gfHp93xmhw/Zp6DgUjRMYmCmyZ84hmc/2nV4fOxvaOOobjSY72ih99oaTSmheOD3DrhsZ5JyqQEXRA3IkmIoKeB17K0X9ucP3aetJafLuL4eXOQVJpPWv++XR2ZPPRxe2ycPafvUB4bGJed4vBla01VJR6xe1iIiLoeWBP5yCr6itYUVeR0/lXr6zF7/NIG4BF8MLxAar8Pq5eWZvT+a215ayoK5fBcxH87mg/Xo/KeTXk83q4emWtFBiZiAj6EplIpXm1ayjn2TmA3+dl22rxoy8UrTUvHh/kxnX1C+pzvqOjntdPDZEWP/qCeP7oAO9ZueySvkRzsW1VHUd7o0Tj4k40AxH0JfJ2T5iR8WTO/nOD69fUc7R3mPMj4wWyzHmcHBjlXDiWs7vFYEdHHRfGJjjRP1Igy5xHMBzjcCias7vFYHtHHWkNb0p7C1MQQV8iu7Ptcm+YoePfXFy/NjMASGl67jx/NFOxmKsLwMDwo78ubpeceeqlU3g9ins2t8x/8hSuWlGL16PEj24SIuhLZNexTArdQrff2tIeoKLUy8snpa9LLkyk0jz9cjdXr6zNOVZhsKKunJZAGa/K4JkTgyPj/ONrp9l5VSsr6xf2rCv9Pja11kiBkUmIoC+BfacvsP9smA9uaV3wz5Z4PWzvqGP3iUHx7ebAT946x7lwjC+8d92Cf1YpxY6OOl7rGpLmUTnw5O5TjCfTfP72hT9ryLQBOHA2zHhS2kQXGxH0JfDYrk5qK0r42I6Vi/r5nVe1cvr82GTzI2FmUmnNY7tOsrGlZs6Of3OxY009gyPjHOsbzrN1zuLCaIJvvdLNB7e0sraxalGfce3qOsaTaQ6ei+TZOmE+RNAXyZFQlN8c6ee+Gzpm3GAhF/7Vllbaast5dNdJmTnOwbPvhOgaHOWh965DqfkLXGbi/RuX4/d5eHpPd36NcxhP7TnFWCLFQ4tYCRlsy2548XKnxCyKjQj6Inls10kqS7188oZVi/6MEq+HB25Zw77TF3hDgkgzkk5rHnm+k3VNVdx5ZfOiP6ehys+Ht7Xz4zfP0S+l6TMSiU3w9J5u7t7czIbl1Yv+nIYqP9etqeO7r5+RLQCLjAj6IugeHOXnbwf5xHWrFhwMnc5Htq2gvrKUR3d15sk6Z/Hbo/0c7R3mc7etzan8fC4euHktyXSab+45lSfrnMXTe7oZHk/y0O3rl/xZn7l1LaFInJ8dCObBMiFXRNAXwddfPInP6+H+mzqW/FnlpV7uu3E1u44NcFj6dr8LrTVf+90JVtSV8/tbFx54ns7K+gru3tzCd189I4Uv0xiOT/DUnlO874rlbGydu2NoLty2oZHLllfzxItd4k4sIiLoC6Q3EudH+3r4yLZ2mmpmb/i/EP74+tVU+X089sLJvHyeU9h9YpADPRE+d9s6fAuoDJ2LB29dy/B4ku+8eiYvn+cUvvXKaSKxCR6+Y/G+86kopXjgljUc6xtm1/GBvHymU/ja707wcmdh0pVF0BfIN3Z3kdbwmVvW5u0zA+UlfHzHSn7xdpDT52ULL4OvPd9Jc00Zf3BNW94+c1NbgJvXN/DUnlPEJyStDmAskeSbL53itssa2ZJDm9xc+eDWVloCZXxdJiqTJFNp/vuvjxesj5MI+gIYGk3w3dfOsHNr64KLW+bj/ps68Hk8fP3Frrx+rl15res8r58a4jO3rsHvm3mP1sXy4K1rGRge58dvnsvr59qVlzvPMzSa4E9vXpPXzy31ZdySr3YNsf9sOK+fbVf6hsdJ60zTuEIggr4Ant5zithEis/elr/ZuUFTTRl/+J52frS3R7IwgJ/sD1Lt93HvtYvL8Z+LG9bWs7ktwBMvnpSNL4Bz4RgA65cvLu98Lu7dvpLqMh9PvCizdMj0yAFoCeTHXTsdEfQF8MN9PbzviibWLyGlay4+c8saycLIci4co6OxkvLS/M7OIePf/exta+k+P8Zzh3rz/vl2IxiJUeJVNFTOvwPUQqny+/jEdav4l4O9dA+KO9EQdFNn6EqpO5VSx5RSnUqpL81yzm1Kqf1KqUNKqRfya6b5JJJpQpE4m9oCBbvG6oZK7trUwvdeO+P6zIBQOFawWQzAB65spqOhksdfkKKuUDhOc6BsyWmhs3HfDasp8Xh48iVxJ4YimdW3aTN0pZQXeAS4C9gIfFQptXHaObXAo8Dva62vBD6cf1PNxdgnsTVQmJHV4JpVy4jGk4Rdvj1dKBKnpYDP2utR/OnNa3i7J+J6/24oEivos27KBrZ/uLeHQZe3iw6GY1SX+aguy73H/ELIZYa+HejUWndprRPA94Gd0875GPBjrfUZAK2145qTTI6stYWbNQK0Zkdu43puJBqfYGQ8WdAZOlxsedw14G5XQDAcn/zeFYpP37yG8WSa773m7nTRzLMu3OCZi6C3AWenvO7JHpvKBmCZUmqXUmqfUupPZvogpdQDSqm9Sqm9AwP2yk0NRQobzDBonhT0WEGvY2VCYWPwLOxqSJ51pvFZXzRe8Ge9rqmKNQ2VHHV5c7RQJFbQSWEugj6TY22609EHvAe4B/gA8FdKqQ2X/JDWT2itt2mttzU2LmyTArMJZkWmucAuFyNY4uYZuiGwhZ41lpV4WVZR4upnPTgyTjKtC/6sITOAhsLuHTwh83tdqIAo5CboPcCKKa/bgekNGnqAX2qtR7XWg8CLwNb8mGgNeiMZ31fVIjsr5kpDlR+fR7l61njRvVXYwROgJVBOr4sF/WIanTzrQhNLpBgaTRR08MxF0N8A1iulOpRSpcC9wDPTzvkpcLNSyqeUqgB2AEfya6q5BCOF9X0ZeD2K5TVlrp41hsIxPAqaqvOfRjedlkAZQTc/6yLFhiDzrPuGx12b+3/RbWviDF1rnQQeAp4jI9L/pLU+pJR6UCn1YPacI8AvgbeB14EntdYHC2a1CRTa9zWVzNLUvSITjMRprPZTkqf+LXPRHCij18WrIWOG3laE1VBzoIxUWrs208UYPAvpcsnJf6C1fhZ4dtqxx6e9/grwlfyZZi16I3E2FzAHfSrNgTJXd14sdBrdVFpry7kwNkF8IkVZSf6LmKxOKBKnvMRLoLwwaXRTaZmSwbU8T43t7MTFoiJzXS6uZzyZYnAkUTyRCZQRDMdcW/ASCscL+qWfSnNWWNzq2zVWnovdCWohGFlFbl0RXUysEEE3ld4CV3dNpzlQzngy7criIq11wYuKpmL8nwZdLDLFiA3BRd+xW+NDoUiMhqrSvDebm4oIeg5cLNct3gwd3CkykdgEsYlUEQdPmaEX61kvqyih1Odx7bMOFjhlEUTQc2IyOl3EoCi4U2SMZWmhv/gGbp41TqTS9A+PFyU9FDJN0dycVVTo/kQggp4ThW6oMx03FxcVqyLXoLzUS21FiSsHz75oHK0LX8A1leYad2YVaa0Jhgsf7BdBz4FQOE6gvISK0sIWFRm4ubgoWGT3FmRExo3PupgFXAatteWunKhE40lGE6mCB/tF0HOgmH5GcHdxUSgcw+dRNBahqMjArSIzmUZXzBl6oIy+aJy0y4qLJttZiA/dfDJZF8XNm3VrcVFvNkfZW6De3DORKS5y37M2Y4beEihjIqU5P5oo2jWtQLFaLIig50AoUvhudNNpCZTR68Kt6IJFXg0BtNSUcX404bpNo0Ph4vQnmopb8/4vBvvF5WIq8YnCN9SZiRaXFheZMXgaWUV9LhtAzxUxB93gYlaRu2IWoUgMr0fRVC2CbirGTKLQbXOn48biIqOoqNiDp1uziorZn8hgMiXXZYNnKBynuQiuRBH0eQgWqTf3dNxYXHR+NEEimS5oafRMuDXvv5gVuQb1laWUeJXrBs9zRchBBxH0eSnW7jnTcaPITD7rorsB3Dd4muVK9BgZXC7b6KJYrkQR9HkwlobNRe4O50Y3wMXUruI+64pSH4FydxUXmZHhYtAacFeaaDqt6Y0Up+GcCPo8BMMxllWUUF5a3NaqbiwuKnbPnKm0BNyV9x8yIQfdoNllGVznRxMkUumiBKBF0OfBDD8juLO4KBiJUer1UF9ZWvRrNwfcVS0aNHGGbgyebsngupiDLjN00zGjqMjAbcVFoXCc5kAZniIWFRm0uKy4KFREkZlOc6CMRDLNBZdkcBWrShRE0OfFjNQuA7cVF4UisaJnuBi0BMoZHEkwnnRHcVEwEqeustSUXZou7lzkjhVRMTuIiqDPQSyRIjw2YYrLBdxXXJTZbMG81RBAf9Qd+12GIrGiB58NjJoOt6yIguEYfp+HZRWF3+ZPBH0Oit3KdTotLiouSqc1fdHiV4kaTKYuuiSdLhQ2JzYE795b1A2EshtbFGObPxH0OTAz6yJzXffkRw+OjJNMa9Nm6C0uq2AMRmKmPeuGKj9eF2VwBYu4GhJBn4Nib2wxHTcVF5nRB30qzS7auWhkPMlwPGnaasjrUSyv9rviWUNxV0Mi6HNgZAKYFagzgihu2LJrMuvCJL9uld9HdZnPFRWMZma4GLTUlrtiojKRStM3XLzYkAj6HAQjcepNygSAi8VFbtiyy+wZeuba7sj7N551sfZtnQm39KCf3OavSM9aBH0Oek1MowN3FReFipgJMBstgXJX+NAtMUOvcUdxUbFbLIigz4FZVaJTcUtxUShavEyA2XDTDF0pWF7k/kRTaQ6UEZtIEY0lTbOhGBR7mz8R9DkIhs3L1TVwS3FRqEjtReeiOVDG4Mg4iWTaVDsKTSgco6naT4nXvF//yY0uos52JwaL3K1VBH0WRseTRONJU10u4J7iIiushloCZWjt/J2LrPCsm12Six6KxKgp4jZ/IuizYHzRir1F13TcUFyUTKXpixanvehcGCLn9BVRMfOiZ2OyuMjh7sRgOF7U4LMI+iyYXSVq4Ibiov7hcdLa3AwXcEe1qNba1CpRg8ZqPx6F4zO4QkXe9FwEfRbMrhI1MHxvTk7xssrg6YZCrkhsgthEyvRnXeL10OiC4qJgOFbUAi4R9FkwloLLA35T7bg4Q3fuF/9iape5IlNdVkK13+dokTkXLl4r1/lweppo/3CcC2MTtImgm09vNEZDlR+/z5yiIgM3FBeZtZfoTDi94GUy68LkGbphg5MHz394uRul4O7NLUW7pgj6LATD5m1sMZXJ4iIHB4+CkRiVpV5qyoqTCTAXTt+56GcHglSWelm/vNpsUxw9eI6OJ/n2K6f5wMZmOhoqi3ZdEfRZKHYwYy6aHTyT0Vrz1pkwK+srTS0qMnDyrPHs0Bi/eCfEx3asLFoa3Vy0BMqyjcKcl8H1/TfOEo0neeDWNUW9bk6CrpS6Uyl1TCnVqZT60hznXauUSiml/k3+TCw+E6m0ZWbo4OziotdPDbH/bJiPbV9htilAxu0zMDLORMp5xUVP7u7Co+D+m4orMrPh1I0uJlJpnnrpFNtX13HNymVFvfa8gq6U8gKPAHcBG4GPKqU2znLefwGey7eRxeaZ/UFGxpPcelmj2aYAzi4uevyFk9RXlvLhbVYRdGcWFw2NJvjB3rN86Ko204vlDJwa8P/F2yHOhWM8cEvxB85cZujbgU6tdZfWOgF8H9g5w3lfAP4Z6M+jfUUnndY89sJJLm+u5vbLmsw2B4ANy6sZT6Z5peu82abklSOhKM8fG+C+G1eb1tFyOk5NXfyHl7uJT6T5TJFdAHPRXGM8a+fELLTWfP3FLtY1VfHey4uvH7kIehtwdsrrnuyxSZRSbcC/Bh6f64OUUg8opfYqpfYODAws1Nai8KvDfXT2j/DZ29ZawqcL8MGtrTRU+Xn0+ZNmm5JXHn/hJJWlXv74utVmmzJJ+7KMG+Bo77DJluSPsUSSb73SzfuuWM66JvODoQbLa8pQylnl/y91DnIkFOWBm9fg8RRfP3IR9Jmsmr72/zvgL7TWc26ZrrV+Qmu9TWu9rbHRGu6MqWiteWxXJyvrKriniKlG81FW4uVPb+7gpc5B9p8Nm21OXjg7NMbPDgT5+HWrCJjYMnc6axur2NhSw1MvnSKVdoaL65/eOMuFsQketNDsHKDU56Ghyu+oytyvv9BFU7WfnVe3mnL9XAS9B5jq4GwHgtPO2QZ8XynVDfwb4FGl1IfyYWAxefnkeQ70RPjMrWvwmdiJbiY+ft0qAuUlPPJ8p9mm5IVv7O7C61H82xs7zDblXSileOi96+gaHOXZd0Jmm7Nkkqk039h9im2rlrFtdZ3Z5lzCNStree5QnyMyXQ6ei/BS5yD33dhhWv1KLqr1BrBeKdWhlCoF7gWemXqC1rpDa71aa70a+BHwOa31T/JtbKF5dFcnjdV+/vCadrNNuYQqv49P3bCaXx/u45jN3QGDI+P84I2z/MHV7ZYJ0E3lziubWddUxdd+10na5rP0X7yTCdB95ta1ZpsyI5+7bR2R2ATffvW02aYsmSde7KKy1MvHdqw0zYZ5BV1rnQQeIpO9cgT4J631IaXUg0qpBwttYLE4cDbMns7zfPqmDssE6KbzqRtWU1Hq5dFd9p6lP72nm0QqXfQc3VzxeBSfv30tx/qG+c2RPrPNWTRaax5/oYv1TVXcYUKALhe2rqjl1g2NPLn7FGMJ+252cWE0wS/eCfHR7SsJlJvnQszJr6C1flZrvUFrvVZr/TfZY49rrS8JgmqtP6W1/lG+DS00j+7qpKbMx8evW2W2KbOyrLKUT1y3ip8dCHL6/KjZ5iyKkfFMgO7OK5tZ21hltjmz8sEtraysq+CR5zttmy66p/N8JkB3izkBulx5+I51DI0m+M6rZ8w2ZdGcHhojldZct6beVDus5Sg2ic7+YZ471Mcnb1htiQq6ufj0TR34vB4ef6HLbFMWxfdeO0M0nuRBi7oADHxeD5+7bS0HeiLsPjFotjmL4s0zF4BMlpSVec+qOm5YW8/XX+wiPjFnXoVlMVIvzW4wJ4IOPLari7ISD5+6YbXZpsxLU00ZH9nWzj/v67FlrvSvDveytT3A1hW1ZpsyL39wTTstgTK+9jt7urhCkTgNVaWWdSFO5eE71jM4Ms73XrfnLN0q7bZdL+ixRIqf7j/HR7atoL7K3Fa5ufKZW9aS0ppv7LbfLD0YjrPGwq6WqZT6PDx461pe7x7iNRsWdfVGYpYMOs/EdWvq2b66jq+/0MV40n6z9N5InFKfh2Ump+C6XtCDkRjJtObqlbVmm5IzK+oq2HlVK9997Qwj4/YJJKXSmr6odXrk5MIfXbuChio/X7NhumgoEqe5xvyWxLny8B3r6Y3G+eHeHrNNWTDBSOZ7bXYxousF3Uq9uBfC721cTmwiRfegfYKj50fGSaa1rQTdKOrafcJ+RV2hiPn7tC6EG9fVc/XKWh7bdZJE0l7N0Xot0p3V9YJu7NVp9mbQC8UYgOxUNh20iJ9xody7PZNXvKfTPsHRsUSSSGzCNi4XyBR1PXzHes6FY/zvt+w1Sw9FzN+nFUTQJwOLZm81t1CMaLqdNmOwSibAQgmUlxAoL7HVs74YpLPXs75tQyMbllfxk7emF6Nbl3TWlWiFwdP1gh6KxGioKjV9q7mF0lDpp8SrJrcUswNBm7q3ILvxhY2eda9NV0NKKS5vrpnc+9QOnB9NMJGyhivR9YKe2cjCXl96yFQz2m27tN5oHL8FMgEWQ1ttua36dhsNr6wgMgulpTbzvbZL2wXjd9BoB2wmrhf0kI1Su6bTEii31awxGI5ZIhNgMRgiYxcmXYkWEJmF0hooZyKlGRwdN9uUnDDcW6215k8MRdAjcVptK+hlk0FdO9BrkcDRYmgJlBMemyCWsEeOdCgap77SHkVF0zFWFXaZrBiDpxUmhq4W9MwGtUlaLDCyLoaWQDl90biNlqb2ykGfipH+Z5cBNBS278rTmOnaZUUUisQp9Xqoqyg12xR3C3rIxn5GyIjMREozOGL9pWkqremNxm2X4WIwmSZqk1mjvQfPzLO2S8C/NxJjecBviQZo7hZ0m2YCGNgpF31wZJxUWk/u9G43jDoFu8zQe6P2dW8tqyjB7/PYZiejYCROi0Uqcl0u6PaeoU/6Gm0gMsYvp13jFUadgh1m6LFEivCYvYqKpqKUorW23BYTFcjGhiyy8nS1oAfDcZSyZyYA2GtpaqXA0WLw+7w0VPltMXgaNtqp7H86rbX2CPhrremNWKOoCFwu6JmiIj+lPns+BmNpageRMXK47dZiYSoZkbH+4GnMbO3UmGs6dknJHRpNkEilabHIpNCeSpYn7JyyCBeXpnYQmd5IDL/PQ60Ni4oMWgJltvDr2rXsfyqtgTL6huNMpKzdpGty8LTIRMX1gm7XwJFBc03ZZLaOlQlG4rTWltuyqMggM2uMWX5LOqNnjlXcAIuhpbYcraEvau3JitUGT9cKutba1rm6BpkKRmt/6SHjQ7dCafRSaK0tYzSRIhq3dg/6YCROnU2Ligwu5qJb+7tttYZzrhX04fEko4mUrQNHkPFJ9w+Pk7T60jQcs8yXfrFcTBO19orIEYNndqJldRdXKBLH51E0VFqjW6trBd2uG1tMp6W2jFRaM2Dh4qJUWtM3PG6ZZeliMQZ/qwfr7LaxxUy02CSDqzcSZ3lNmSWKisDFgh60eQ66wWTBi4W/+APDmaIi2w+eNikusnPDOYMqv4/qMp/lV0NBi+xUZOBaQZ+codu0j4uBHTa6sHsBl0FTtR+PsvYM3SgqsvvgCdmWxRZ+1oClctDBxYLeG4nhUbC82hq+r8Vihx4jdm+xYODzelheY+2Cl96otbIulkKLxfv9a62z7i3rfK9dK+jBSJym6jJ8Xns/gpoyHxWlXkuLjNVSu5aC1XcuMlJYrTRrXCwtteWWDoqGxyYYT6YtFYC2t5otASf4GSFTXGQHkSkrsXdRkUFLbbmlZ40hB1TkGrQGyrhg4R70VpyouFjQ7Z8JYNBqA5FpCdi7qMigLds0yqrFRSEHFBUZWL0vuhWftSsFPVNUZP8qUYOMr9HCM3SLZQIshZZAGePJNEOjCbNNmZFQJM6yihJbFxUZWL09tBVjQ64U9EhsgthEykEiU87AyDiJpDWLi0IWywRYClYXGTtv8zedyV2iLOpH743E8XoUjRZKrHCloFtxZF0KrbVllu17kUyl6R8ed4RPF6wvMkEb71Q0nebJalHrfa8hoyPLq/14LVJUBK4VdGv1X1gqVp41DmR3KpJnXRx6HRLsB+v3oO+NWu9Zu1LQg2HrRaeXQquFi4usmAmwFOorSyn1eiyZJhqfSHFhbMJSedFLxco96K0Yh3OloIciMbweRVO1M0Sm2cLl/07pmWPg8SiaA2XWfNaTG1s443sNRt6/9QZPo6jIahOVnARdKXWnUuqYUqpTKfWlGd7/uFLq7eyfl5VSW/Nvav6wou9rKVi574VTyv6nYlWRceazzhQXWS1NNBpLEptI2c/lopTyAo8AdwEbgY8qpTZOO+0UcKvWegvwZeCJfBuaT0Jh52RdGLQGrLmpbigSp7zES6Dc/kVFBlbdwNgp/Ymm0lZbbske9KGoMXha61nnMkPfDnRqrbu01gng+8DOqSdorV/WWl/IvnwVaM+vmfklFIk56ksPxkYX1ps19maXpU4oKjJoCZTRG42TSltr1mj0cXGUy8Wi8aGLW89Z61nnIuhtwNkpr3uyx2bjfuBfZnpDKfWAUmqvUmrvwMBA7lbmkcmGOhb7j1gqVt1UNxix/8YW02mpLc/0oB+2Vg/6UCRGbUUJ5aX2LyoysGrzuZBFEytyEfSZplYzTk2UUreTEfS/mOl9rfUTWuttWuttjY2NuVuZRy5kG+pYbam0VFoDZZwfTRCfsFbfi8zuOc571mC9vuhWzLpYKkYG1zmLxSyMbq1NFioqgtwEvQdYMeV1OxCcfpJSagvwJLBTa30+P+blH6MgxGoj61IxXEi9FvLtJlNp+qLO6ZljYNlZowWzLpZKU3UZXo+ypMvFit1ac7HmDWC9UqpDKVUK3As8M/UEpdRK4MfAH2utj+ffzPwxmRftMB+6FWeNAyPjpLX1/IxLxap5/71R5wm616NorrFeN9HeqDUTK3zznaC1TiqlHgKeA7zAU1rrQ0qpB7PvPw78NVAPPJoNfiW11tsKZ/biMXbpdpoP3fhyWemLb+RqO6Xs3yBQXkJ5iddSuejxiRRDownHCTpkVtNWmqhAZmK4vqnKbDMuYV5BB9BaPws8O+3Y41P+/Wng0/k1rTAEs7t011dZy/e1VKy4I70V24vmA6WU5bKKeiezLpw1eEJmNX3gbNhsMybJdGuNcfP6BrNNuQRrOYCKQCgcY3lNmWOKigzKS70sqyixVH60ITJOm6FD5p6sVJJ+cWMLZw2ekHFx9UbipC2SJjo8nmQ0kbLk99p9gu6gjS2m02Kx4qJgOFNUVFOe00LQVlitWtSpqyHIDJ6JVJrzFulB/9P9mZwQK6bjulLQnbgshWwjIwuJTG80k4PupKIig5Za6/SgH0+m+NWhPsCZgm7EBcx2caXSmr/5xWH+6icHuXFdPXdcvtxUe2bCVYKeTmt6HVhUZGClGXrPhTEOnI04MkgH0GaRHvR90TgffeJVfnmoly/esZ6KUuethozukWZOVqLxCe7/hzf4xu5TfPL6VTx933ZLFnA5739/Dn5zpI9EKu1YkWmpLSMSm2AskTT1F3vXsX7+7Af7SaU0n711nWl2FJKpfdFX1FWYYsPe7iE++503GR1P8sjHruGeLS2m2FFoLgq6OYNn9+Aon/7WXroHR/mbf72Jj+9YZYodueAKQU8k03zluaN8Y/cpLm+u5p4trWabVBCMIM0z+4Pcu31l0a+fSmu++pvj/M/nO7lseTWPfeI9dDRUFt2OYmBmLrrWmu+8dob/+LNDtNaW84/37+Cy5uqi21EsllWU4Pd5THnWwXCMnY/sQSn49v07uH5tfdFtWAiOF/Qz58f4wvfe5EBPhE9ct5L/556NjthAdybuuKKJ96xaxpd+/A6vdp3nP31oEzVlxelyODgyzhe//xZ7Os/z4fe08+UPbXLsc4aLM/RfvB3ifVcsp9JfvF+lb796mr/+6SFuv6yRv/ujqwlUOKeT5UwopWitNSer6KUTg0RiE/zsoZvY3B4o+vUXiqN96D9/O8g9f7+brsFRHvv4NfznD212tMhUl5Xwgweu4/94/wZ+9naIu/5uN6+fGirKte9/+g32dl/gv/7hFr7y4a2Ofs4AlX4fX7xjPb8+0sedX32R17qK1+3it0f6Wd9UxTc/ea3jxdygrbacfd0XOHguUtTrHg5FqSj1cmVrTVGvu1gcK+id/cM89N23WLe8imcfvpm7NjvTvzgdn9fDw3es54cPXo/Pq7j3iVf42+eOFTSHNxqf4EBPhM/fvo6PXLti/h9wCH/+/g3802euR6G49xuv8uWfHy5Kc7TDoShb2mvxOKyWYi4evmM9Ka3Z+cge/ssvjxatCd2RUJTLmqtt86wdK+j7Tmfas/+Pj1xlWtDKTK5ZuYxfPHwzH7q6ja8938nuzsGCXetIMApgiyVpvrl2dR3/8sWb+cSOVXzzpVPc8/e76ewfKdj1+ofjDAyPs9EmM8Z8sb2jjt/8+a38wdVtPLbrJHd/dTdvdBd29am15nAoyhUt9nnWjhX0d85FqC7zsarefWJuUOX38Z92bkIp2H8mXLDrHMoKul2Wpfmm0u/jyx/axD/ev4P+6DiPPN9ZsGsdCQ0D7nzWgYoSvvLhrXz7/u0kUmk+/PgrfPU3Jwp2vXPhGMPxJBtF0M3n4LkoV7bWOLKoZSFU+X2saajknXPhgl3jUDBKY7XfMZtuL5ab1jewY00db/eEC3aNQ8GMD9lOs8Z8c/P6Rp77s1u44/ImHnuhk4lUYYq7jMHTTs/akYKeTKU5EoqyqdV9LoCZ2NwW4J0CBpMOBSOunDHOxJb2WroGRxmOTxTk8w8Ho7QvK3fUHq2LodLvY+fVbcQn0hzvGy7INY6EoigFl9soJdSRgn5yYJTxZJpNbSLoAJvba+mLjtNfgKrG+ESKE/0jIuhZNrcH0DqzQiwEh0NRW7kACslV7bUAHDhbmMnK4WCU1fWVRU1JXSqOFHRjNrqpTb74kJmhAwWZpR/vGyaV1rIayrIl+6wL4XYZSyQ5NTjquoDobKyoK2dZRUnBWuse6Y1yRYt9ZufgUEE/eC5CRamXjgbrNaA3g0wsoTCCfjEgKoIOUF/lp622nLcL8KyPhIbRGpmhZ1FKsaW9lgMFGDxHxpOcPj/GFc32etaOFPRDwQgbW2oc1/N8sVT6faxtrCpIUcbBbDbRijpndrBcDFtXBAoyQz8cyg6e4kqcZGt7gON9w4wlknn93KPZZ22ngCg4UNDTac2hYFT859PY3Bbg7Z7CzNA3tkg20VQ2t9VydijGhTz37z4cjBIoL3Fst9DFsHVFLekCxCyOZAXdbu4txwn6qfOjjCVSEqSbxqa2AP3D+Q2MptKao71RcbdMY2u2wCrfbhcjICqD50W2ZAOj+V4RHQ4NEygvsV1nVscJuuFWcGPV4lxsac9/YLRrYIT4RFqCz9MwXCLv5FFkkqk0R0NR280YC01jdSZmsT/PgdEjoUxA1G6DpyMF3e/zsK5RAqJTyczs8ivoEhCdmUB5CR0NlXl1cZ0azKTiSkD0UrauCOQ1MGqsPO3mPwdHCnqUy1tq8Hkdd2tLwgiMvpNHkTEGz7WNzux5vhS2tOc3ZnExIGo/kSk0W9szMYuhPMUsus+PEp+w5+DpKNXTWnMwGGGTLEtnZEueK0YPBaNc3lwtg+cMbG4L0BuN5y1mcTgYpdTrYa2sPC/B8KPna5Z+xKYZLuAwQT8zNMZwPDlZSCO8GyMwmo99MLXWmfRQcbfMyNYVtQB5m6UfCkbZ0FxFiQyel7C5PYBS5K3A6Egois+jWL/cfoOno74dRuqSpCzOjBEozofbpedCjGg8KQHRWbiytQaPyk+mi9HG1Y4ugGJQ5fexvqkqb4Pn4WCUdU1V+H3226TFWYIejFDitefIWgw2tmREJh9uFwmIzk1FqY/1TdV5Safri44zNJoQQZ+DLe21HDgbRuulb+RyJDRsS3cLOE3Qz0XYsLzaliNrMchnxeihYASvR9mqE12x2dwe4J2eyJJF5nAo8/8lFaKzs3VFLedHE/RcWNpG0kOjCXqjcdv1cDFwjKBrrTl4LiJNouZhc1sgL26AQ8EoaxsrHb936FLY2h7g/GiCc+GlicyhrCtRBs/ZuSpPgVE7B0TBQYIejMS5MDbBJikompPN7QEG8hAYzfRAl2c9F5uzIrPUmMXhUJRV9RVUl7m7B/pcXNZcTanXs2Q/ugi6RTDcCJKyODeTrXSX8MXPDAjj0l5hHq5oqabEqziQB0EX//nclPo8bGytWXLF6OFQlKZqPw1V/vwYVmQcI+iHzmV8unYdWYvFxtalB0aNbdBkhj43fp+Xy5qrl7T933B8gtPnx0TQc+CqFbUcPBchlV58zMLOAVFwkKAfDEZZ11glPt15qCj1sa6paomCbs9OdGawpb2Wt5cQGD3am90UWtJD52VLe4CxRIrO/pFF/Xwimaaz396Cbp+9lebhnXMRblnfaLYZtmBTW4DdJwYX/fOHg1FW1Mm+lrmwpS3Ad187Q/f5MToaFtYiQWvNoezAu7FFVkPzYRRzHTgb5rJ5AsjhsQRvnrlAIqlJa00qremLxplIaVtPVGwn6OfCMfZ2D73rWCyRYmB4XIpccmRzW4Afv3mO7sFRGqr9aK3RZHrJDwyPcy4cIxiOEwzHCEXixJMpyE4wNZqXTgxy47oGU+/BLkxt7zpd0LXWnBwY4Y3uC7zRPcTe7gv0RuIZgdEaY1JfV1nK8hp7+nSLSUd9JdVlPvb3hPnItSsueX8skeQ3R/p5Zv85Xjg+wETq0lWTR8HV2YHBjthO0PefCfPF7++f8b1tq+qKa4xNMUTmtr/dNed5Xo9iebWfslIvisyWXwpoCZSz86q2QpvpCNYvr8Lv8/CV547xnVfPkM4Onqm05vT5US6MTQDQUFXKtlV13L25Ba8HPEqhlMKjMr5hu7VxNQOPR7GlPcCLxwf4n789QVqTed5ac3pojF8f7mMskaK5poz7buzgjsubqPT78HkVXqXweBQ1ZSU0Vtt38FS5+PaUUncCXwW8wJNa6/9v2vsq+/7dwBjwKa31m3N95rZt2/TevXsXbPDIeHLGlLuKUi8tAdkGLRe01vzgjbNEYhMoBQqV+VspGqpKaastp7W2nKZqvzTeygP/49fHebXrfFakmfy7uaaMazvquHZ1HavrK0S088CTu7v4z7848q5jSkFteQl3bmph51WtbF9dh8fG21MqpfZprbfN+N58gq6U8gLHgfcDPcAbwEe11oennHM38AUygr4D+KrWesdcn7tYQRcEQZiLRDKNUuBVFycqTmIuQc9l+rUd6NRad2mtE8D3gZ3TztkJfEtneBWoVUq1LMlqQRCERVDq81Di9eDxKMeJ+Xzk4kNvA85Oed1DZhY+3zltQGjqSUqpB4AHsi9HlFLHFmTtRRqAxadpWAu5F2vilHtxyn2A3IvBqtneyEXQZxripvtpcjkHrfUTwBM5XHNug5TaO9uSw27IvVgTp9yLU+4D5F5yIReXSw8wNQeoHQgu4hxBEAShgOQi6G8A65VSHUqpUuBe4Jlp5zwD/InKcB0Q0VqHpn+QIAiCUDjmdblorZNKqYeA58ikLT6ltT6klHow+/7jwLNkMlw6yaQt3lc4k4E8uG0shNyLNXHKvTjlPkDuZV5yykMXBEEQrI9UjQiCIDgEEXRBEASHYDtBV0rdqZQ6ppTqVEp9yWx7FoJS6imlVL9S6uCUY3VKqV8rpU5k/15mpo25oJRaoZR6Xil1RCl1SCn1xexxO95LmVLqdaXUgey9/MfscdvdC2Qqu5VSbymlfp59bdf76FZKvaOU2q+U2ps9Ztd7qVVK/UgpdTT7O3N9oe7FVoKebUPwCHAXsBH4qFJqo7lWLYingTunHfsS8Fut9Xrgt9nXVicJ/Dut9RXAdcDns/8PdryXceC9WuutwFXAndlMLTveC8AXganNTOx6HwC3a62vmpKvbdd7+SrwS6315cBWMv8/hbkXne1GZoc/wPXAc1Ne/yXwl2bbtcB7WA0cnPL6GNCS/XcLcMxsGxdxTz8l0+vH1vcCVABvkqmEtt29kKn/+C3wXuDn2WO2u4+srd1Aw7RjtrsXoAY4RTYBpdD3YqsZOrO3GLAzy3U2Zz/7d5PJ9iwIpdRq4GrgNWx6L1k3xX6gH/i11tqu9/J3wP8NpKccs+N9QKbS/FdKqX3ZliFgz3tZAwwA/yvrCntSKVVJge7FboKeU4sBoTgopaqAfwb+TGsdNduexaK1TmmtryIzw92ulNpkskkLRin1r4B+rfU+s23JEzdqra8h4179vFLqFrMNWiQ+4BrgMa311cAoBXQV2U3QndhioM/oTJn9u99ke3JCKVVCRsy/o7X+cfawLe/FQGsdBnaRiXPY7V5uBH5fKdVNpiPqe5VS/4j97gMArXUw+3c/8L/JdH214730AD3ZVR/Aj8gIfEHuxW6CnksbArvxDPDJ7L8/ScYfbWmyG5p8Eziitf7vU96y4700KqVqs/8uB94HHMVm96K1/kutdbvWejWZ34vfaa0/gc3uA0ApVamUqjb+DfwecBAb3ovWuhc4q5S6LHvoDuAwhboXs4MGiwgy3E1mw42TwH8w254F2v49Mi2FJ8iM3PcD9WQCWSeyf9eZbWcO93ETGVfX28D+7J+7bXovW4C3svdyEPjr7HHb3cuUe7qNi0FR290HGb/zgeyfQ8bvuR3vJWv3VcDe7HfsJ8CyQt2LlP4LgiA4BLu5XARBEIRZEEEXBEFwCCLogiAIDkEEXRAEwSGIoAuCIDgEEXRBEASHIIIuCILgEP5/mnaEG8o+ed4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(prob)\n",
    "plt.ylim([0,1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}