Hadamard Walk 1D.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "033af634",
   "metadata": {},
   "source": [
    "<b>Hadamard walk. </b>\n",
    "\n",
    "Reference: \"1D and 2D quantum walks with light transverse momentum\", Laurea Magistrale in Fisica, Chiara Esposito. (page 6)\n",
    "\n",
    "#### Coin operator\n",
    "$C \\ = \\ \\frac{1}{\\sqrt{2}}\\begin{pmatrix} 1 && 1 \\\\ -1 && 1 \\end{pmatrix} $\n",
    " \n",
    "#### Translation operator \n",
    " $T \\ =  \\ \\displaystyle\\sum_{j} \\ \\lvert j \\ + \\ 1\\rangle\\langle j\\lvert  \\ \\otimes \\ \\lvert R\\rangle\\langle R\\lvert + \\ \\displaystyle\\sum_{j} \\ \\lvert j \\ - \\ 1\\rangle\\langle j\\lvert \\otimes \\ \\lvert L\\rangle\\langle L\\lvert  \\$\n",
    " \n",
    "#### Protocol \n",
    "$U \\ = \\ T (C \\otimes I_p) $ "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "6b0e70a3",
   "metadata": {},
   "outputs": [],
   "source": [
    "from numpy import *\n",
    "import matplotlib as mpl\n",
    "from matplotlib.pyplot import *\n",
    "import pandas as pd\n",
    "from scipy.sparse import csr_matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "df75285d",
   "metadata": {},
   "outputs": [],
   "source": [
    "N = 100  # number of random steps\n",
    "P = 2*N + 1   # number of positions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "e8982d3d",
   "metadata": {},
   "outputs": [],
   "source": [
    "ket0 = array([1, 0])  # |0>\n",
    "ket1 = array([0, 1])  # |1>\n",
    "\n",
    "C00 = outer(ket0, ket0)  # |0><0| \n",
    "C01 = outer(ket0, ket1)  # |0><1| \n",
    "C10 = outer(ket1, ket0)  # |1><0| \n",
    "C11 = outer(ket1, ket1)  # |1><1| \n",
    "\n",
    "I = C00 + C11  #identity"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "09b179f8",
   "metadata": {},
   "source": [
    "#### Coin\n",
    "\n",
    "$C =\\frac{1}{\\sqrt2} [ |0⟩_c⟨0|+|0⟩_c⟨1|-|1⟩_c⟨0|+|1⟩_c⟨1| ]$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "6d9f5773",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Coin\n",
    "C = (C00 + C01 - C10 + C11)/sqrt(2.) #coin operator for 1 qubit"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6969fb56",
   "metadata": {},
   "source": [
    "#### Translation operator \n",
    " $T \\ =  \\ \\displaystyle\\sum_{j} \\ \\lvert j \\ + \\ 1\\rangle\\langle j\\lvert  \\ \\otimes \\ \\lvert R\\rangle\\langle R\\lvert + \\ \\displaystyle\\sum_{j} \\ \\lvert j \\ - \\ 1\\rangle\\langle j\\lvert \\otimes \\ \\lvert L\\rangle\\langle L\\lvert  \\$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "db0b53cd",
   "metadata": {},
   "outputs": [],
   "source": [
    "LPlus_ = eye(P)\n",
    "LPlus_ = LPlus_.tolist()\n",
    "LPluss = [row[1:] + row[:1] for row in LPlus_]\n",
    "LPlus = kron(LPluss, C00)\n",
    "\n",
    "LMinus_ = eye(P)\n",
    "LMinus_ = LMinus_.tolist()\n",
    "LMinuss = [row[-1:] + row[:-1] for row in LMinus_]\n",
    "LMinus = kron(LMinuss, C11)\n",
    "\n",
    "Tx = LPlus + LMinus"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "75c429fa",
   "metadata": {},
   "source": [
    "#### Protocol \n",
    "$U \\ = \\ T (I_p\\otimes C) $ "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "70a1510c",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Protocol\n",
    "U = Tx.dot(kron(eye(P), C))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5ea5f5db",
   "metadata": {},
   "source": [
    "### Initial state\n",
    "$|ψ⟩_0=|coin⟩_0 \\otimes|posn⟩_0= \\frac{1}{\\sqrt2}(|0⟩c+i|1⟩c)⊗|0⟩p$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "4dfa560e",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Initial state\n",
    "initialPosition = zeros(P)\n",
    "initialPosition[N] = 1     # array indexing starts from 0, so index N is the central posn\n",
    "psi0 = kron(initialPosition,(ket0+ket1*1j)/sqrt(2.))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7fb98843",
   "metadata": {},
   "source": [
    "### Final State\n",
    "$|ψ⟩_N=U^{N}|ψ⟩_0$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "40916b37",
   "metadata": {},
   "outputs": [],
   "source": [
    "#State after N random quantum steps.\n",
    "psiN = linalg.matrix_power(U, N).dot(psi0) "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8f5cf958",
   "metadata": {},
   "source": [
    "### Measurement operator and position on the map\n",
    "$M_k=I_{c} \\otimes |k⟩_p⟨k|$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "078a2982",
   "metadata": {},
   "outputs": [],
   "source": [
    "prob = empty(P)\n",
    "PositionMap=zeros([1, P])\n",
    "for k in range(P):\n",
    "    posn = zeros(P)\n",
    "    posn[k] = 1     \n",
    "    M_hat_k = kron( outer(posn,posn), eye(2))\n",
    "    proj = M_hat_k.dot(psiN)\n",
    "    prob[k] = proj.dot(proj.conjugate()).real\n",
    "    PositionMap[0][k]= prob[k]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "9380f285",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAD4CAYAAADhNOGaAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABK30lEQVR4nO29e5wdV3Xn+1tV59VPSS21ZFmy9bDlh4BgHMU2AZzwSiyTRJnhznzsMJjxJBG+154Jn2EmmJDMTT4zk0xyE+7gXMePgAETiIEQiAIGQwDbYPyQ/JIty7JaL6v16G5JrX6fV9W6f+zau3bVqTp9JPfROa1e38+nP+d01T51dtXZe6+9HnttYmYIgiAICxen1RUQBEEQWosIAkEQhAWOCAJBEIQFjggCQRCEBY4IAkEQhAVOptUVSGLZsmW8du3aVldDEARh3vDss8+eYOb+s/lsWwqCtWvXYseOHa2uhiAIwryBiA6d7WfFNCQIgrDAEUEgCIKwwBFBIAiCsMARQSAIgrDAEUEgCIKwwBFBIAiCsMARQSAIgrDAEUFwDvjZvhPYNzLZ6moIwoJmpuzhG88OQlLv1yKC4Bxw5zdewn2P7Wt1NQRhQfOjV4fx8a+/iIMnp1tdlbZDBME5oOL5qHoyCxGEVlLxfABANXgVQkQQnAN8ZviijgpCS9F90JeuWIMIgnOA5wOiEAhCa/ECCeCJJKhBBME5gEUjEISWo7ug9MVaRBCcA3xmiVQQhBajBYB0xVpEEJwDfAZ88U8JQkvxRSNIRQTBOUCcxYLQekJnsfTFOCIIzgHMEqkgCK2GJWooFREE5wDPF41AEFqNjhaSvliLCIJzgJiGBKH1GB+BqAQ1iCA4B4hpSBBajywoS0cEwTlAwkcFofXoLih9sZaGBAER3UBEe4hogIjuTDhPRHRXcH4nEV0dHL+ciF6w/saJ6GNzfA9tj5iGBKH1iEaQTma2AkTkArgbwPsBDALYTkTbmPkVq9hmABuCv2sB3APgWmbeA+Aq6zpHAHxzLm9gPiDrCASh9cg6gnQa0QiuATDAzPuZuQzgIQBbYmW2AHiQFU8BWExEK2Nl3gtgHzMfesO1nkdoNdSTxicILcWXvphKI4JgFYDD1v+DwbEzLXMTgL9P+xIi2kpEO4hox8jISAPVmh/okDWxSwpCa/GlL6bSiCCghGPxJ1m3DBHlAPwGgK+nfQkz38/Mm5h5U39/fwPVmh+E6mhr6yEIC50wfLS19WhHGhEEgwAusv5fDeDoGZbZDOA5Zh46m0rOZ2RZuyC0B9IX02lEEGwHsIGI1gUz+5sAbIuV2QbgliB66DoAY8x8zDp/M+qYhc5nWDQCQWgLJMVEOrNGDTFzlYjuAPAIABfAA8y8i4huC87fC+BhADcCGAAwDeBW/Xki6oSKOPro3Fe//QlT30rrE4RWogWA9MVaZhUEAMDMD0MN9vaxe633DOD2lM9OA1j6Buo4rxF1VBDaA4kaSkdWFjcZ7ZiS/bIFobVI4EY6IgiajJiGBKE9kL6YjgiCJiOmIUFoD3xJQ52KCIImI+qoILQHso4gHREETYZFIxCEtkC083REEDSZMGSttfUQhIUOGx9BiyvShoggaDImZE1sQ4LQUnQXlPDRWkQQNBnZJ1UQ2gNPTEOpiCBoMiymIUFoCyTFRDoiCJqMOKgEoT3Q0UKyjqAWEQRNRgSBILQHpi+KSlCDCIImI+sIBKE9kL6YjgiCJsMyCxGEtkDW9KQjgqDJyIbZgtAeiJk2HREETSYMH21xRQRhgeOJaSgVEQRNRmYhgtAeSF9MRwRBk5F1BILQHkiKiXQaEgREdAMR7SGiASK6M+E8EdFdwfmdRHS1dW4xEf0DEb1KRLuJ6O1zeQPtjsxCBKE90OsIJHCjllkFARG5AO4GsBnARgA3E9HGWLHNADYEf1sB3GOd+wyA7zHzFQDeCmD3HNR73iCCQBDaA9mqMp1GNIJrAAww835mLgN4CMCWWJktAB5kxVMAFhPRSiLqBXA9gM8BADOXmfn03FW//ZEc6ILQHsg6gnQaEQSrABy2/h8MjjVSZj2AEQCfJ6LnieizRNSV9CVEtJWIdhDRjpGRkYZvoN0RjUAQ2oPQRyB9MU4jgoASjsWfZFqZDICrAdzDzG8DMAWgxscAAMx8PzNvYuZN/f39DVRrfiDb4wlCeyDZR9NpRBAMArjI+n81gKMNlhkEMMjMTwfH/wFKMCwYRB0VhPZA+mI6jQiC7QA2ENE6IsoBuAnAtliZbQBuCaKHrgMwxszHmPk4gMNEdHlQ7r0AXpmrys8HbDVUVFJBaB2SYiKdzGwFmLlKRHcAeASAC+ABZt5FRLcF5+8F8DCAGwEMAJgGcKt1if8I4MuBENkfO3feY88+fAbcJCOaIAhNR7KPpjOrIAAAZn4YarC3j91rvWcAt6d89gUAm86+ivMbe/bh+QzXEUkgCK3ArCMQOVCDrCxuMrYgEJVUEFqHRPClI4KgyfgRH0ELKyIICxxfUkykIoKgydgLyWQmIgitQ1LCpyOCoMmIaUgQ2gMxDaUjgqDJxKOGBEFoDbr/eZLupQYRBE3GXjsgYWuC0DokxUQ6IgiaTFQjkAYoCK1CTEPpiCBoMl7ER9DCigjCAkfWEaQjgqDJSIoJQWgPRCNIRwRBk/FFIxCEtkDWEaQjgqDJyDoCQWgPZB1BOiIImoysIxCE9sBsVSmqeQ0iCJqMPfbLdpWC0Dp0X5T5WC0iCJqMaASC0B6IszgdEQRN4NRUGTf878dx6ORULHxUGqAgtArP2jaWmfGhzz6FH+4eanGt2gMRBE3g0MkpvHp8Aq8NTUqKCUFoE9g4i4GKx3hi4CReOjLW2kq1CSIImoDtlJJ1BILQHoThoyy7lcVoSBAQ0Q1EtIeIBojozoTzRER3Bed3EtHV1rmDRPQSEb1ARDvmsvLtimdWMHKkoUmbE4TWYU/QtJnIk8kZgAa2qiQiF8DdAN4PYBDAdiLaxsz2JvSbAWwI/q4FcE/wqnk3M5+Ys1q3OdUgPKjqc2Twl7A1QWgdvmUaqgb/VKVPAmhMI7gGwAAz72fmMoCHAGyJldkC4EFWPAVgMRGtnOO6zhtMThOfJWpIENoEtqKGtKYupiFFI4JgFYDD1v+DwbFGyzCA7xPRs0S0Ne1LiGgrEe0goh0jIyMNVKt98SwVVLaqFIT2QI/5zHYfbWGF2ohGBAElHIsPafXKvIOZr4YyH91ORNcnfQkz38/Mm5h5U39/fwPVal98y/4oaagFoT2ww0d9673QmCAYBHCR9f9qAEcbLcPM+nUYwDehTE3nNZ6ldoppSBDaA3tBma21C40Jgu0ANhDROiLKAbgJwLZYmW0Abgmih64DMMbMx4ioi4h6AICIugD8CoCX57D+bYlpZMzRFBPS5gShZZh1BD4kaijGrFFDzFwlojsAPALABfAAM+8iotuC8/cCeBjAjQAGAEwDuDX4+AoA3yQi/V1fYebvzfldtBkRjcAXjUAQ2gE/4ixWxzxP+iTQgCAAAGZ+GGqwt4/da71nALcnfG4/gLe+wTrOOzwrNC3iIxCVQBBahi0IdIi3aAQKWVncBPyUqCGRA4LQOux1BLKyOIoIgiYQiU6QFBOC0Bbo/sfMJmxUNAKFCIImYBxRvmxVKQjtQjhBs/uodEpABEFTiDilZB2BILQFui/aJlvpkwoRBE2g6qf5CKTRCUIr4Fg/NAEdEjUEQARBU/CtqCEWjUAQWo6tmbOVdE76pEIEQRNIXUcgeU0EoSXENXPZyD6KCIImoLVNz1rKDsjsQxBaRVwQhCuLW1Wj9kIEQROwU9xKiglBaD22Ns4MSUMdQwRBE0hPQy2NThBagd0PbU1dTEMKEQRNwE5oJesIBKH1pJuGpFMCIgiaQjQNtXVcGp0gtIRozi+IsziGCIIm4EXCR8U0JAitJt4P9foBEQQKEQRNICndrX1cEIRzix8L2pCVxVEaSkMtnBleZGVxeFzWEQhCa6j1Eaj3ohEoRBA0AXtjbEa0AQqCcO6JbhAFiRqKIaahJmBvjM0MZBwCAIgcEITWoMf7jEOyeX0CDQkCIrqBiPYQ0QAR3ZlwnojoruD8TiK6OnbeJaLniejbc1XxdiaedM4NBIFEDQlCa9ADvhsIAjugQ2hAEBCRC+BuAJsBbARwMxFtjBXbDGBD8LcVwD2x878HYPcbru08wY/5CLKueswy+xCE1qD7XtZ14PscCfEWGtMIrgEwwMz7mbkM4CEAW2JltgB4kBVPAVhMRCsBgIhWA/gAgM/OYb3bmvjK4oyrNAJpc4LQGvQcLOMS2PYRyOQMQGOCYBWAw9b/g8GxRsv8bwC/D6BuzAwRbSWiHUS0Y2RkpIFqtS/2Nni+z5aPQBqdILQCrRFkYqYhieRTNCIIKOFYfERLLENEvwZgmJmfne1LmPl+Zt7EzJv6+/sbqFb7Yie0sn0EooYKQmvQA7/yEcjK4jiNCIJBABdZ/68GcLTBMu8A8BtEdBDKpPQeIvq7s67tPMFWO30GMo72EbSyVoKwcAmjhhyVdE5yDUVoRBBsB7CBiNYRUQ7ATQC2xcpsA3BLED10HYAxZj7GzJ9k5tXMvDb43I+Y+d/N5Q20I/aCMrY1Aml0gtAS2IoaYlsQyOwMQAMLypi5SkR3AHgEgAvgAWbeRUS3BefvBfAwgBsBDACYBnBr86rc/nhWjLJvrSOYTRCMTVeQzzooZN2m11EQzheYGcMTJazoLaSWia4jENNQnIZWFjPzw1CDvX3sXus9A7h9lms8CuDRM67hPESrm1WPkXUbjxr6rc8+hXdcugx/cOOVza6iIJw3/GzfSdzywDN44hPvwQWLkoWBcRa7yllclfDRCLKyuAn4NRpBY+sIRiZKGJkoNb1+gnA+MTJRguczRqfLqWXCqCEnskOZ+AgUIgiaQCTpnB9qBLO1Oc9nWekoCGdItQF7v72OAAAqkoY6ggiCJmDsj4EtstHw0Yrno+pJYLMgnAm6z1Tq9B07fBQAqsECAgngUIggaALRHcrYchbP/jnRCAThzGhEI7AXlAGiEcQRQdAEoknnAIcaSzpXtXKgCILQGI0kkNOnXCMIfHNcVvyLIGgK9u5Heh2BQ7M3uKrPddVbQRBq0X1Gbz+ZhO57OgGkXVYmXyIImoI9Q/EZgSCguvZIvchFGqUgnBlhf0ufRNnrCICoP0Eih0QQNAXdHnW6WyItCNI/I/nRBeHsOBMfgRuEclcsjUASz4kgaAp2riFmhkMAUf0IBd2YJWpIEM4Mbeap1DENxZ3FtvYgGoEIgqYQ37zeCTSCeu2tkVmNIAi1eMGgXlcjCMZ9100wDUmfE0HQDOLhow4BDtVvcJ4npiFBOBuqDfkIksNHAREEgAiCpmCnuPUZykfg1HcWV/zZIx8EQaglNKs2YhrSUUOiEdiIIGgCdmZDNhqBMg1te/EofueLO2o+00jkgyAItVTrLA67//F9+KNvvRymmEjQCGR1sQiCphD1EYTrCHxmPP/6KB5/rXYrTvERCMLZoX0ESWbVZw6M4ol9JyLZRwHxEcRpKA21cGbYm9dHw0cZ7IVmIJswX4o0SkE4Eyp1tOmq76Pqcc06AltoiCAQjaAphGmoVdZDhwgUrCOo+j6YaxufaASCcHaYQIuESVTVY1Q930o6p9cRhEJDTEMiCJpCNWYa0lFDvs9mxh9PJSELygTh7Kg3iap4PiqBrw6oTUNtf34h05AgIKIbiGgPEQ0Q0Z0J54mI7grO7ySiq4PjBSJ6hoheJKJdRPQnc30D7Yi96YVeR+AGUUPaBBRvfCZfijiLBeGM0H0m0eTqqz5XYxqyNQIRBLMLAiJyAdwNYDOAjQBuJqKNsWKbAWwI/rYCuCc4XgLwHmZ+K4CrANwQbG5/XqN9BHodAQVRQz5b9swUjcATH4EgnBHVOn2n6mkfQXxlseUjENNQQxrBNQAGmHk/M5cBPARgS6zMFgAPsuIpAIuJaGXw/2RQJhv8nfdPXY/xVZ/BDLhEJsVEmlO4KqYhQTgr6i3GrHiMiu9bUUNqyCtXJWrIphFBsArAYev/weBYQ2WIyCWiFwAMA/gBMz+d9CVEtJWIdhDRjpGR2vDK+YTtfKr6fiTFRNU02qhGkHZcEIT6VP10s6qOGtJdMr5DGSBJ54DGBAElHIuL0NQyzOwx81UAVgO4hojenPQlzHw/M29i5k39/f0NVKt9sWcYFY/hODDrCELTUFwjSI+FFgQhnXradNXjyIZPoY9ATEM2jQiCQQAXWf+vBnD0TMsw82kAjwK44UwrOd+ICgI/koY6bX9V/RlmcV4JwplQz7+mHci6v2mNoBxZUCYqQSOCYDuADUS0johyAG4CsC1WZhuAW4LooesAjDHzMSLqJ6LFAEBEHQDeB+DVuat+e2ILgqpnpaH22TIBxTQCqxEnRT8IgpBMJSUSDwj7VSnwCSRqBNLdZl9ZzMxVIroDwCMAXAAPMPMuIrotOH8vgIcB3AhgAMA0gFuDj68E8MUg8sgB8DVm/vbc30Z74TEj6xIqHhsfgQ4fjc9QNLLSURDOjnp5unRQhnYOa2dx1fdNH5X+1mCKCWZ+GGqwt4/da71nALcnfG4ngLe9wTrOO3yfkcs4qHie8hFYKSaqKasgbfVU/ASC0Dj1FpRp4aBNQXbSuc6ci4rnycpiyMripqA0gvDR2ikm0haORVY6yloCQWiYejuU1ZiGrH6p+6hoBCII5hxmFaoWFQTqj5nN7CXeaCN+BfERCELDzJZiArBMQ04Y4GgEgWgEIgjmGt0Yc27Y4BynNmqoNnxUfASCcDbUS0Otj2lB4FqCQPdRidITQTDn6IaXy4SPlqytKk3SuZoFZZaPICYkjpyewRV/9F28NjTRrGoLQlvDzHjXX/wIX9txuOZc6Hfzaz6jJ1WlqgcgqhHoPio+OREEc452PNmCwLG2qjQLx+poBPGGeWR0BsWKj9dPTjer2oLQ1lQ8xuFTMzh0cqrmXNqCMtv8Go8aAsI+KhqBCII5R89Aan0EsRQTKQvK1Pu4Izk55FQQFgqVOhs3eSk+AtvXFo8aAsRHYCOCYI7Rbc8WBC6FW1WaBp2Shlq9j57TjbgsgkBYoMSdvknn4hOlZI0gQRCIRiCCYK7Rs4uca/sIdPhoGDVUXyOICYJqeicQhIWA6QMJk6FUjcAqm+wsDkxDohGIIJhrjGkoY0UNGY0AqQvK6vkI6qnFgrAQ0AKgkjAZSvMR2P+HpiFrHUHQR0XRFkEw59TzEfi+lWIiJQ21ei8+AkGwSdviFbBDspP7DRAuKHOTfASybkcEwVyTZBpyHJVrqOKHedHPJMVEpRqNhRaEhUY9rThtQZndx7QgyLq1piGZX4kgmHN0KFp8HQERRQbyM0k6VxJnsbDA0X2nlDAZ8lJNQ7U+Aodq1xFI1JAIgjknXFkcW0dAQDlY1ALUt2fWRD9UxTQkLGzKKeZRjgRg1FtH4IEoxVksUUMiCOYaPbuoDR+lyIw+bs+M5kdPcxaLIBAWJmmToXo5uuw+VfZ80w81Ej4aIoJgjvETooZ0iomoaegMfAR1YqgPnJjC914+9sYrLghtwMtHxvCTvbV7lqc5i+uZVCsx05DWzDW6j0r4qAiCOaeaGDVENT6CmjTUsV3NbMrVdEfZl548hI9/7cU3XnFBaAPueXQf/u9tu2qOlz1lVo1PhqIm1XRncbnqG1+dJms2qRFBIIJgjvESnMUOKfNQRBDENQIvXcUt612WEkxDMxUPRYkmEs4TZioeSpXa9lzWkXN1+s1sC8r0ToEa4ywWQdCYICCiG4hoDxENENGdCeeJiO4Kzu8koquD4xcR0Y+JaDcR7SKi35vrG2g3/JTwUceJRjzUzF7qqbh1FtOUqz48n2t8DoIwHylX/cTIoDQ/WbWeSdX6v1T1zb4gGnEWh8wqCIL9hu8GsBnARgA3E9HGWLHNADYEf1sB3BMcrwL4ODNfCeA6ALcnfPa8IilqSKeYqNZzbAV7qAJJ0Q/p4aM6va6ElgrnA6WqZ9q0Tbog0KZYSnAWR4WEE3MW5yTpnKERjeAaAAPMvJ+ZywAeArAlVmYLgAdZ8RSAxUS0kpmPMfNzAMDMEwB2A1g1h/VvO7RGkM3URg3ZJGkEhYxr3tuU64SPSh4i4XyiXPUT23JaO9d9pZBxI2YioLaPOQ7B7oaubExjaEQQrAJg7wYxiNrBfNYyRLQWaiP7p5O+hIi2EtEOItoxMlIbNTBf0GN10laVNjVJ5zxGPpu85N1kH63WNthSnYU2gjDfKAWmIeYU82hCvwGAfNapu6AMCFO9aNzAZyAaQWOCgBKOxZ9c3TJE1A3gGwA+xszjSV/CzPcz8yZm3tTf399AtdqTxK0qEzSCpEabDzSC+EymXp4V0QiE84m0CDkTMBFr5zpENJ9x664jAGr7oRukfhGramOCYBDARdb/qwEcbbQMEWWhhMCXmfkfz76q84OkqCGKzUSAZFtnIZscxVCpM9hre2qSXVUQ5huhhhttz2m5hnRfKSRoBPE+RqSCNjQOEVwiSTqHxgTBdgAbiGgdEeUA3ARgW6zMNgC3BNFD1wEYY+ZjpIJ2PwdgNzN/ek5r3qYkrSyOL2QBkpLOsdEIanwEdVYW63NiGhLOB0opkx4zGUpZkZ/PuGCO2vvj/ShuGsq4ohFoMrMVYOYqEd0B4BEALoAHmHkXEd0WnL8XwMMAbgQwAGAawK3Bx98B4MMAXiKiF4Jjf8DMD8/pXbQRSUnnHKdWI6hZUOb5xkdwJmmodcy1CALhfCDUcJP9ZJ6vNqTX6wF0P9J9p+L7yDvBhCrWX+KmoXCfEPERzCoIACAYuB+OHbvXes8Abk/43E+R7D84b0nejyCukiaruLNFDcUX06hj4iMQzh/SfF62JlDxfLhOtK/ovmObVXUf05tCxYM2Qh+BCAJZWTxHvDQ4hvsf35e8H0GwjkDTmcskrCNgS1VNyzVU6wcQjUA4X2Dm1Ci4StUe4MNzto8AiO/0p8p15tR8V4WPJkcN/XjPML7x7OBc3s68QgTBHPGtF47gz7+3J0w6V8dHUMi6iXlRMmYDm+TZUFKuIdEIhPOFpM3mw3O1ewvYx41/zavVCApZdS7eD/WGUb7P+LsnD+FvHh2YozuZf4ggmCNmKh48P5zRJOUa0nTknNo01D4j4zrIOFS7MKaaHj5aqkjUkHB+YLfheHtOy9xbqxHU5vPqyKlzTmw/AtdR/bLqs8rZlZDjaKEggmCOKAYD8mSpCiC6JR7FTEMdWbfGD+D5PjIOIeNQahhc3aihhEb8reeP4MXDp8/ibgSheYwXK/jMv+yt3YmvmjzrB6Jt335vfATZWh9B1VcZR7WZNtFZHGgExUpyaouFggiCOUILgumyFgS1m9drOlJMQ65DyLhOQvbRZLup77O5TlKuof/xnd148MlDZ3tLgtAUHtszgv/3X17D7mPRtaX24J8WNRR/H4aPOpH/AaU5ZB3H9EWVhjq8ZsZxjI9gpuJjpiyCQHiDaLVysqQaUz5T30eQZBrKuk6iszgt15DdIbSJyGamXDUCShDaBT3gxttmqY4gSPMR6MVgWiOIOIs9HxmXkAm089qVxQgWlDFKQTr3eGqLhYIIgjlCN+7pUq1G4DpKBdV05JJMQ0ojyDqUmn00rkXY5qC4RsCs7J4zIgiENkO3yXjbjGoE9XwEtaahvDENRc8pc6vqi66TYhoK+opnadgLDREEc0QxaLhT2jQUSzFhq6TKNFQ748k4BNdN8hGo//ViGk3Js5xrlVptwWcsaHVXaE+MICjHNYLw/1ofQXL4aNw0FC+XdR3jr6OYZu4GPjkv8BEAYT9eaIggmCN0o9amoWydpHMdWTcxxUTGVbOXpBQT+uOViDkoXSPQ9RGNQGg30trmbD4C3QfsLLx1ncVe2KcA5auzgzacQEPwgqghACgu0ImTCII5QjfcqcA0lHEc03B1citNPpuQKdFnZJwgfNQ6x8woV310Bg09zWkW7zhps64kfrh7CAPDE7OWE4R6zJQ9fOnJg7Pm90/XCOpHDSX1gdBHULugrOL7qk9ZPgIgDCHVC8qqPhsf30INIRVBMEfoRq0FgWsN/vbSdtch5FxKiBryzZJ3+5xu2F15tTrS3q4yohHEBcEZaASf+MZL+NvHD8x+k4JQh0f3DOOP/mkXdh9PzDRvOCuNoOon9oGKlXQOiO1K5jGyLhl/nWMmZupV++5sgbRQNWgRBHNE3EfgOOHMw15HkNEhookLylSEQzRfiirXrTuBvfoyohFEG3CaQy6JyVIFk0G9BeFsmQgmQVOl+m0urW3W9xH4Vh9oPMWEXqQJhGYhsjQDl8I+C9RGMi0URBDMEWHUkHrVs3sguo4g6ypVtRJTn+0Ih4h6G9hDO/OBWhzRCOyVmNGOU2zQ5ukFarGOdhKEs0W3oalZJhVp9vhSnaihisdhH0haUJaSYiLj1NcIXIdMn7XrttAQQTAH+FZqCb2y2PYLqLA1VTbjErJOVCNg1ql11ezFPqcjg3TirEZ9BNMNmob0AripBeokE+aOqdhkKA0tAKbPwEdQrvphH6ja5p9oGupoignfTLwAWBOz6LqCyZJoBCII5gC7AYfO4nDtAFH4XjuvfA7VWpO62gmdVxptCkpSixvxEVR9rpuQbtp0StEIhDdGOKmo35bSJin1fARR01CDUUN+NGqIqNZpnHHJ9FlABIHwBrAbz3RlNtNQqKrqQV03Zjc4F/ERBB2iK0EQaI2gI+umRg3p91XPTxzsdSeYbRYnCLOhfQNJZsZy1Tf9JN1HELbnpP0IkvqAFzMNxTOYqhQTYT+0X91AI5i26iFRQ8JZYzdoZr2AjGIqqDqvZih6dyWOvOo01LZpSDf6rlyCjyCwo/YUMjV7FRQjjdvDfY/vxwfu+mlN3fXsbLZZnCDMRj0z458+vBu3fO4ZAJb/KkUj6Clkan0EVd/0geiCMvU+m1F9KlEjiIWPau1cB3TYWSXER1AHIrqBiPYQ0QAR3ZlwnojoruD8TiK62jr3ABENE9HLc1nxdiLeeFwrQgiILijLOg4ybnRLSp12Wq8jsE1DepaU6COIdJxk05B+f/DEFA6dnKrJpSIagTBXTNUxMx44MYWDJ6cAzL6yuLuQqVkgWfHY9AG7rascXaH5J+4jUFFDYfZR+zXjOJH1PYCYhlIhIhfA3QA2A9gI4GYi2hgrthnAhuBvK4B7rHNfAHDDXFS2XYk3Hj3jCMNHwxQTGZeMqqrVWL0RTf3w0VrVt2QEQbZGlZ6OmYYmS1X4XOugszWC2RJuFQMTEwC8enwc33r+SN3ywvzl+ddH8b2Xj5n/JxuIKpuuEz46Waqaa8ykOIvLVR9Zl5Sp0zLRMDPKnm/1gahpyHXCPlWbfZSsFBPqeBg1hEgOMEA0gnpcA2CAmfczcxnAQwC2xMpsAfAgK54CsJiIVgIAMz8O4NRcVrrdiAsCs5AseLp2sis164/OXjyfI+eiudfTF5TV0wjs0LzpsoeJouqE+lWjTUI+z77d5W/e/QTu+uFeAMCDTx7Cp775Ut3ywvzlsz85gP/+7d0AgGcPncJVf/J9DI5O1/1MPY1goljBdFlNJOplH825DvIZJ7ZHcZhYLuNQTf/IBumkgdr9CNTkKkw6B0SdxrUagfgI0lgF4LD1/2Bw7EzL1IWIthLRDiLaMTIyciYfbTnxxqNNQm6CjyBr2Sz17CXuI0jSCDrztaYhPXD3JmgEMzEfgV7sM1mqRMrZJqGpUhUzZQ9D48XE+zxwYgoHTqrBYGymgqmyV7MwTjg/GJupYHxGtZX9I1Oo+ozB0ZmachPFCk5MlgDU9xFMFkNtoV720XzWRS7jRDQC3Qd0oEV0hzIfbsQ0FM815CBbEzaqzrtBkkcbMQ2lQwnH4jaERsrUhZnvZ+ZNzLypv7//TD7acvQMx0QJxV+tLfIyVjZEEzUUvOpwNrsx6wE+yVFWMhFFburKYl2/iaLq1OMpGgGgNIe/eXQA/+ruJ2ruUe3g5JvBQb/GNQzh/GC8WMFEqQrPZ9Nm9G9u86cPv4pbP78dQDipSIoa0u3k1HTZtO+klcVKI3BRSgiY0JlEyzEfQcaxJleWj6Di+8g6oUZgRwup/6MageuQCII6DAK4yPp/NYCjZ1HmvEWnl1jSmQUQTWoFJKSYiM1ejEYQRBRVE9JIaNNQPGool3FQSAi3mymH/89UPDMjiw/ctp12qlzF4OgMjo4VazbH0Z8bm4kKlPFi7eCQBjPjB68MzZqUTJhbdhw8hZPBrL1R9KA/WayGwj9B6A+OTuPIaaUpTKVoBL7PJoXJyERYD7uNAlojcAKNwEo3YQmCXMat2aFMB1no/yPnrAVliSkmLB/Bks7sgk3b3ogg2A5gAxGtI6IcgJsAbIuV2QbgliB66DoAY8x8LH6h8xXdeBZ35gDUZjmMp5iIawS2j8CNp5gw4aMJuYaqPvIZBznXSVhHUEVHsMhmphL6CCZrBEH4/1TJw+npMoDa2d+40Si0JqA1g8Y1gmcOnMLvPrgDzxw8r11GbQUz40OffRqf++mZJRW0Bf1YTAu0GZtR55k51AjKtVqnjkMYnlBmx46si5lYuTQfgZ7k5FxHJWyMaQT2gB41qwZJ55xYiokE351DKuiiOIuf7HxlVkHAzFUAdwB4BMBuAF9j5l1EdBsR3RYUexjAfgADAP4WwP+lP09Efw/gSQCXE9EgEf32HN9Dy9GNJ9QIELymrSMI91eteL4Z7HX0g52G2ggCk2soml8on3GQzyYIgrKHvi4lmCaKVaOGT8Rm8HaEx3S5ajr9WKzTh4NBNfo6i0bw0uAYvrr9dQDAcDAb1LPCofGirGhuAsyMAydUqOZ4sYpS1TfPfMfBU/inF+pHezGzpQVUzG8cbxP6mOczJkpVoxHEQ5HtiKPhcVWPvq4cZipeJFKtFNEIap3FuYyDbMaJRQ35QSResEjT91H11JaT1ZQ01E5EI1DX6ci6KGRd0QjqwcwPM/NlzHwJM//P4Ni9zHxv8J6Z+fbg/FuYeYf12ZuZeSUzZ5l5NTN/rjm30jqKMY3ArdEIYlFDQcM8cGIKb/pvj2Dn4BgA5Qxz46ahmpXFVsep+MhnXORcF57PEcftTMXDki4lmGx1PB4GWKMRpAgCe2CwX+OCJc6XnjqIP/nnVwDAaBv69YP3/AyfCaKQhLnjyf0n8e6/fBQDw5PmWY9Oq9/psz85gD99eHfdzxcrvtFKx2eqdYX+6eC6w+NF6Ml4fHGibY7Uk4G+rhx8rl0XY3wECVtTKm06FlEU+Ah0XytXffzi//oRvrr9cLgxjZu8jkDvRwCoFBUdWafG17ZQkJXFc4B2MBmNoGYlYxjDbKeY2DM0gbLnY/cxlb9db50XcRbHwkfjSedygUYQPzdT8bEkEEy2INAq/+5j42BmTJU9Y1+dLof24NM1piH1uXLVx9h0xQio8RnlUPzW80dq/AoAcGpKhQ0WK54ZjEan1SzyyOkZHD5VPyRROHP0Mx0cnTbPXAuEU9NljE5VEteMPP7aCIbGi5EB39YI4mZA32dz7uhpZfLJWNk8D56YwnS5GhEEI5YgAIBiOerzymdcZRpKSEmddQk514nsUObFfAQnJ8sYnijh1eMT4VaVseANs57AiQoC0QiEN8RMRQ2m3flAECRoBHbUkG60OkxTv2aCVcdJuYb07kwVz8fzr49icHQapYpnfARA1JFcLHvoymVQyDrGLguoGfxLg2PY/JmfYPvBUUyXqujvyQNQ4aN6hpemEQDAYSuefLxYwdMHTuJjX30BTwycqHk2p6ZKwWvZXHt0uozR6TKYVccFgGNjM3j20Oy+g4HhCfzglSEAwPGxIj54z88wHDy//+9He/GzfaoOY9MVHAzMI+3I6FQZrwehuKemyvj9f3jRrPL+T3//PL738nEAauHej/cMz3q9w6em8cLh0wCAk1Nlc91RoxGEx8qeX6MZej7jt7+4HQ/89EDktx63wkjjGsFEMbT9Hw0cxv09eUyV1eTg1/76p3jgpwciWqNui1oQ2JFD5WowsckoU2e56uORXcdDZ3GCaUj7CCjoY/r6p6bKVkRRctSQ2o9ACwIHHVlX9iwWzp5ixVeqZS5ocInrCHSKibBhhoJADZY6aqji+/jZwAn8r+++ahp9PhuGzt32d8/iM/+yt0YjsNXp6UoVHTkXHVk3ahoqVvF6MGN8/dQ0psqeEQQjk2Fo39h0GRXPx+efOIBy1Y8MAnY8+fhMxdh9tdr/5acP4dlDowBUh9SvoWmoEjkOAH/9owH89heVRZGZ8eyhU2bW+peP7MGjwWB432P78fGvvQAAeOHwKJ49NIqdg2NgZnzmh3vxj88p+/enf7AHH/rs0wDUIPWbdz9hBMa2F49i75DamnNkooRXrR219o1MmvenpsqRvDhj0+EzODY2E3mvZ5LDE0XsOT5hPnPfY/vMb3j7V57Dd3aqGIq/eGQPPvJ5lXvniYET+NqOQbx4+DRKVQ/bXjyKx15T93v3j/fhE/+w01z7jq88ZwbWPccnzO/yV9/fg9u//Jyq92T4bMeMRhAI4eB5j05VMDZdwZ99dzdKVQ+npsqoeIzhiVJMI6jWhAxve/EoDp+ajkwWjo6pZ9vfk4fP6plPllRbs4WObotaW7VNk9rnlQucxY/sOo6PfulZvHxEmU7z2lns+fj2zqP428f3o+r7JvjCdcj0pZNTpWAfcDvpXK1pSPfVjpxoBMIbZKbiKdUyyIAYX1nspKSY0I1Wdw5t62QGvvXCEdz3+D4T3qntozMVD8MTJRwfLwY+glAjmChW8Vff3xMsDAuEU0wQTBSrZgHQickSpstVLA1mZ3pWByiN4On9p/An//wKfrJ3JNLpByMaQfR6zIz/8e3d+PJThwBEBYE9Oz05GRUEx8eKOD1dQbHiYfvBUXzwnifx7KFRMDPue3wftr2gopHVQFVFseKZ+xqZLGFspoKKx6Yur5+axtGxGVQ8H8+9PooXDp82AuO/fP1FPPCEiqL59A9ew79/QMXB7zh4Cu/9q8fMwPNrd/0Ef/PoPgBqoP039/0MAPDsoVG8/c9+hD3HJ8DM+PW/fgJ3/3hAlXvkNfyHL6jrPbpnGH/23Vex/eAplKoevrPzmBFog6PTOHxqGr4f1nlksmSei7634fEiTkyqQe3JfSfx7Z3HsHNwDL7P+OA9P8N9j6n6HRsrYniiCN9n80xPWs/8dODU1f+fnCrh0deGcd9j+/Hi4bHIb2ibgMZnKsYsODajTIK/99Dz+OLPDuL0TNmUMxpBd948f3W9coppSGnPMxUPX3jiAAaGJyIaAXN4Da05ZTOBj6Dq46vbD+OBJw6YFBOA6j9hIIJ6zTp2GmpEXh0n7KuFjBIEsrJYOGtKFU+plsGiL7OOwN6PgELTUDamEejG61q7KR0bK4JZzTaJYHZaGhpXx4fHS5ZGoL73p3tH8Nc/GsBP9o6gWPFUJETONaaCJZ1ZTJasgXuihOmSh658Bp05NzLLPT1dwXFLY7EHh7hGMGJdb7JUNcKq4vlmEFGCIPQR6MFqdLoM3+dwUJ8omQHgcGDj1jNVIBp5NJLyHlCDqjY92RqLHsz0QDE4Oo3j40WUqz72B6akgyenMFWq4uhYEfuGlYawd3gS+0am4Ptsju0fmTTP80CQUG3wtBJAVc+PCPqRWP2Hx0uo+ozTM5X69zFRgs/q+YX3UcTodBmTpSoOn5ox91vx1PWMaWgyfOaezzgyOmMcuqPTZdP+hieKoTCKaQT6ewA1idDPdSh4lhrddpYFguDQSS0ISiZkeUln1tStr0uVG52q4I//+RU89MzhQCNQK4sBmPUJR4Nrq3UEyjQ0PK6eVbnqm4E+KgjUvWUzdaKGLJOt0ggcWVAmnD0zwaCrB+QwPC20TdqmIa0R6JmSvWBGN0w9wzpyegZZ1wGREgRHRsOOr51rWiM4FAygx8eKmKl46AxMQ5qVizoiKQFOTJYwVa6iK5dBZy6DY6dDX8LYTMXYW4cnlANR101rBPmMg/FiBScmyuZ6evAbGi+a2ScQNw2VcSp477OarYbfVYqYzOw6AIgMqCOT4fcmDaD6NRQiRVO/+PVOTJaM6Uh9b3RAGRpXs/ITU6WIb8dcTw+q44EAmiqb7xgaL9YKssmwHrbwCn+bqGYwYj2X4fHocwaAEUvo2BrBaes32HciNHvZAnJoPPq92gSUzzimHboOBWZAXYeiMTcBMG1HmxkPnVKC8cRECRPFCoiAFb0FU177CLQAPT5eVOsIMo7ZiF5POHSb14EWZY8xPFFENdCm9ECfsSKKdN/KWP0tcR1BcDCfUX1FBIFw1hQD05AedM3WeFZOk3AdQZh0Lo6OGgLCKIwjozNmoM9nHDNLOjVVxnQpWJIf+AgOm5n0DDyf0ZFz0ZkLBcGFiwuYKFUxYgZudY3OvIuuvGuu7TqkBIE1UIzPVLBykerIuoOuXtKB8ZlqZBCxB3Q9IOn66kik09MVY8cGgJOTpcjAZ8/qhqzZvOezcT7XaAST2jasfBvmepPFSJ3CwTk+0EcHdXN8Ihz41GvJHBuaiAqP9OuVzLMcmSii4vnm2dh1j9/TTDnMETUyWYoIkyHrniLlrOd+aqpkNAIARpMB1EzfFpBamNufWbWkw/zWKxeptnNsLPxerRG4Dpm2owXBYcs0NF6sojuXQW8ha75fCwLt0B8aLwYTGyfUCIIJh762NoNOlaqmjkdPF2s0cJus1d/i63qIwv5mfASxtQ0LBREEc4DWCArZaLxyxDRkoobCvChxMpYg0NEUR08XTcfIumQGOEB1kHxW2VSBUB3XnUuHxOlrL+vO1/gIbI1AO5svXFzA6ZmKZboqYrxYxeolHQBsQdCpNALrenogG5upmEFDnxubUVrFeLGCkcnw3L6RSRMpNTJpzXytgfb0dAVHT4emjejsOfxez2fsH5ky19MmBP3eOBMD/4geUIbGo99rawKlqmdMGvYAPzReDIXCeBEzZc+YwiJaxUQJI0G5k1NlHLeei12/kYnwnsqeH3Fcj0yUIgJtxKqD7QManijWOOh1O9xvRVGdnCpHtAn9vT6rdlTIOljWnTe/9UVLOtU1gjoNjxeNILhwccG0HaMRBG2x7Pk4enoGPYUMegoZ8/061Fq3VW2es9uz/m7d5rVpyH5+OmIPgAkTtbH7WyTraCwRXSGjTLs+R9fqLBREEMwBxYpqwFojqN0k2zYNhT6COBnXgRs7V/bCqIj45+wl+UDoXNOrSjssLaU76Ii2aUjbjDvzUc1hTV+XMgNYjreJmQr6ewooZB1MlqrIuQ76e/IYn4kKAj3zBWCiZzIO4cCJKTArLYIZOHgidDi/GpQD4iaQ0KQCAK8cG4+US9IIAJh1GfqcPRvXg6nPwO5j4fcOTxQxZGki+jPlqo8BayY9NF6KagdBuVLVx97hCatcMWJG0ffBHKufJdDse0q639C/EGos02XPmFcA1Qb0JOJkELKrBbjWCDIOYXSqHDOzhd+7/8QkegtZ9Bayxj9wUZ+6xt7gGlNlD0dPzyCfcbCiJzT5aB+BdvACyufSXcigOxAE+Yxjsunqtjo0VlI+AjfUCOKr5U3SuVjGW5NmOmGCZfc3e6vK+GStI+eafrQQ9yQQQTAHGMdsNho1lJ5iImywFy4KO5GtEdjohpwkQNQMSn2v7jhaIHTmXOPAVjOyLIoV3wwo2pygNIJQc1i5qIDT05WImWNspoLeQqje93ao92MzlUgE0FHL4fxqMJCtWdppZqPrlnUBUFrAqsVqcNlTIwhCAWQPjLuO2gNjMdGkosqNheUmo+YgW1Dp6CBAm2/C+9WOckClydDUaAR2Oft61sBtC6P4fRwfK5oZ74nJ6ID8ytFkgWabhuL3q5/lqsUdyrk7UTLPXP8GFy/txClLIxgaL0Y0zf0jU+jtyKK3I5zBrw40gr1DoVDcOzSJRR1ZLOoITT5Lu4O0JlbI6MGT0+gpZI1G0JlzzboY3VbLJkzaNe05TtZNnkTpPqNNQJE+ZfU3J0EjiKwsDvpASQSBcCa8fGQMR07PGB+BFgQmbDQhxUTWDbeqBIBLV/SY925MEOi32kegZ0r27D3nhlEWGr0WwPZb9OSz6LYymOpBWF9PJ7Vb3JnF4s4sTs+ENmRt1lGDQyAICmqgmCp7qPqMNUs74TPw2tCECc/TM/31/d1moNaD0rGxIi5d3g0gHLw6c64a4CeiA61+Dq8EA153PoP9I1MoV330deUwVfbw+slpLA7MDXomXcgq5/rotDJJnZgsRcwKO60BXoVf2iak5AH++FgogGzNAYgKjMHRaWOiUc7i0JatBUEh62Dv8AQ8n7E4iKg5PlY09nOtOXTnM3j91BQmS1U4pDUCS7AcGTfPTz/LDSvUsz0+XsSavk4QKWHSmXNx4aIOHD41g6myByIlWE5OlrBmqRrs40IfQKhVWOaq14YnsLgzi0XBc+/IuhHzj25j5apvJiK6nB507VX0QJBYzmrPkbaeiQoCe18BIBQIdp+yo/QigiAWQWSHf48XK/jxq7Mv4jufEEFwljAz/v3nt+N/fueVGh9BvGFGtqq0ohgAYEMwGAK1/oNL+tU53TG0QLjigrCh2zZVAGawB0IHGBCahjRXrgyv0ZXPoDNIatcbzPCKFbWyc20wwFd9DswF6ho9HdnIQKHr9MrRcTPYDwyrGePywG4MAOuDcwBwQW8BPfmMMW1cubIXe4cnUfaUoFLmlkkjMPQAesUFPWaw37iyF4ASOlde0GvqoMr1mnIblnej6jP2DE2Y+mmNYPWSDuw+Ng7PZzWTLlVx4MSUcY5rQbByUQG7j6tyKxcVMF6s4uDJKTPgacGyekmHGZzXLevCVNnDwZPT5rfWM/jLL+g1db3ygl4wA68NTZp7euXYOIiAy1Z0m3u/dHk3psoeDpwIv/flo2NwCLhsRY95lna7WtqdN79VX1cOfV05M6CvW9aFsZkKjpyeibQrW+ire1JCohS0CUD5bWyNoCsfTijU7xlerzufMW2zYJlh9DlNvD1fGTwLABEzqEPA2qXqd9QDve53l/aH9551bB8BzKtrFqGpY3b491eePoxbv7A9ommd74ggOEuOBwt9dg6OBSuLbR9B0DATVhbHo4YigiBIQ61504WqExjTUEZdY9WSTtP54jOot6xaZN53ZEPbf2+NIAg7WEQj6MhiUbDqEwDebF2vtyNjaQSZyEBxRTAIj05XcMUFPSZnkh54NOuWhffb151DX3cOzKoOa5d2GSfjz61W37tvZBKXX9AL1yEcGyuiO5/BxUs7TXjgxuAZTZaqWLusEx1ZF6PTFXTnM1iztNM4rN90obregRNT5jN7hyeQdQlXXNCD14KVxm9epc69NjRh7n33sXFkXcLlF/SYWbp+zq8NTZjf6bWhCeQyDi5b0WP8BfoaB05MmQH+2FgRSzqzWLW4YJzV+veYLFVx6fJuZF3CRFEt9lu5qMPcx5svDJ+LfkaHTk6jryuPlYsKJuXDhuXhILykM2ucs/r30DNxfR+npytY399tJik9ltB3CEYo2vcEAIs6cljcoX7fzlwmEq68YUWPmbXb1+vIuiAiU9Zus/H2vNFqp3aermXdeVwQ1ClM36JeL1xcCE2dKZvXh+Zbx9RJT+SeOXgSQFTDO98RQXCWvBzM+AZHZzBRrKCQs30Eqoy9ib1urFk3phGsCAdG1yET+bCkM4tVgToedxb3d+fNLNv2EQDAWy9abN5HnMX5UDUHwoEbiGoEcZvvWyKdPjxnawdAdPa3ordgokfUIGQJgv5QI1hqCYn+njyW94aagx5smIELevNY1h2W67c0DPt7+7vDc/09+Ygmogd4QJkslgbZL/u781jRWzAD48+tVs+v4jHW9CmBW/EYy3sKWLmoYCJK9CBc8Rhrl3Whp5BB1Wes6M1jRW/eRDe9xfreCxd3RO/XcrJq4aTPLQ0WXC3rjt7vm6znsr6/ywx4y2PP5VKrXS3uzJnMuEs6c5Hfw/59+7vD7436g7LG7AYoTVUPmqpNhLZ/x6FInfTCsR7LWazbpJ6Bv2X1IjNbj2sEWsgSRRdcLu8Nf9+4j2Bpd9iu7JX8tl8gTAOjvsc27eogAtuXc74jguAssR2NPqsl6vGoIdeKGgpNQ45JkAWoGbKeNek01IAaTPVAoRt/zuoEutPHZ1BXXWRpBBFncTaigl9hDaBRH0EOi1MEgY4kAaLaARAVLMt7ClgeLB7q68obB6LrEFb2FkzH7evKmfQW/d15k54ACAdac73gWcTL2ZqNLSTiAsOexS7vyZv6Le8tRBY62fe7oreAFYFwUgOPVS4QGPp6+hp2XePfu8IavJb3FCL129jofVgCY4VV9/6e8LlkXcKavk5TTg3+6rda2qW0MI3WlABgWU8ey3rUOdtZ3FvIoiuXMe3UbpvKp6Q+ozPkdgZtqb8nFOA9+Qx6gqSMuk3q/rJ6SYcRQDnXTTQN5fSiykArttuYWVBm2lXetKusE/rkIqaheNSQJQh06PHLYhoSZmPX0fHIjNO2x4fpp611BMY0FDZa1yEs7ghnzK5lz7RntGYdQSbUCPTgoKIswp/x56wBKs1HQASs6eu0nM9h1FBcI7AHivjgoIWC6xAu7us0gsoe8Pq6wvtb3JGF45CZXfbFNAJ7wLMHZHv2198bDuKuQ1i/rNt0cHswtGfcRFGBsby3YK6nZ/CaN8e+Vw+0K3rCQZconKmqa4QCY4X1Gdeh2ABfSBzgO7Iu1i4LB+7+7nAAtQVGzlVmJ1O/nvA+7Ov1dSkNQD+XxZZWtqQrh77gfS7jRDTSZd05E/4ZF/qOQ0ajtJ+Z3V50G9KvtjbTXUcjWNFbwAWLgvZsrSxe0pnFhYEfJBefDPXYGkHUR2BrmhmXEjavt9YRmKghJ2LWWt6Tx+5j45E9Ps5nRBCcJbuOjuHtlyw1oWqFTKjShiqoHbscmobUq4MlnTk4DoWN1gntmWrGEw4A9qstJPJBWmsidX7looIZaG3TkK2a93XmkHEdM2h25Vwzm1tkmQG68xks6sya+sXNBVqw9HWp+1hqDV6hIMibz+voEj1wLO0KTQe2aUPVIWeub5uNbI1gWXcOuYxjBrZ+a0ZrC8ulXfnANJYx3xUKgnBmuaxbzZz1M7NnvvbgtzSYcerfe0VvwcTSL+8JhcKy7hwWdWRNOVursAfuZT05dOYy6AoGxjSNoL8nj8WdWSPAbUFlP6O+rjxch0IBbEX22MJ3eU9etYWgkS7rzoeCIOIPCgVC+PsWaq6ttUpbEOjr2eGjcY3ggt4CLgjuI2etLO7vUW3HdchMgsy9W88lNLs61j2G2lE8DbVDVpoJO2rIEgT/x8+vRrHiRxbhnc80JAiI6AYi2kNEA0R0Z8J5IqK7gvM7iejqRj87Hzk5WcKxsSLefOEiY7PtCOyj+YwTSUNtLyoDwtlLxiWjvvZ15YwNVHdKe9AI1xGE2kKoESiVOZ9xsLw3DyIyncp2Fvfkw0Fcd85lwTU688kagR5Y9KAZDR8NBwpzve7QjLLCmIZCQbKkM/balTVZKO17CmfXBfOaNIDWvHYX0N+tTT5hueWxV9ukYguF5T0FEFHk+0PTUMGqTz4oF/6/PKGuK3oLkXLLrcF6ua29dNfeT+S9dZ6IYvdhC8jwmavnrF9Dv4AtCFb0FuA4FAokWxDYGkHsVd1vgkYQ+Jn0pMIWxt22aSgbEwSLwmebtyZU/T1KoC3tytX6yaznnIn5APq6cmZSkrZVpUtxjSCcNBWyDn79rRcCwIKJHJpVEBCRC+BuAJsBbARwMxFtjBXbDGBD8LcVwD1n8Nl5ATPjn144gi13P4G/+N4eAMCbVvWaKI6Cpe5mjEbgmMZHMY0g4zimQy7tzlmfCXwE1mCgZ0NZSyOwzQX69QJrMFJ1clCwfAT5jNIewlmzeu3IRtcR9BSyILIG0OB6vQXLWdyRRU8+AyIYM0ZozsgnagR6YNI2ZVsjiAqCcIA11+utLRcXQMt6cpZQsAb43nCgj18vPsDb92ubeZTmEK3fBb21WoCtOcQFkK3J9cfuyX7t6wpNNPZn4tfVju74c9HPdGlXHjnXQWfOrYkaitdPmyn1b2ibAfVvvqgjC4dUOKp+los6ssanZGsEOddBb0cmvJ5lmjQaQc6FG6Q+SdQIrGcQX1S5oid8zqGzmNAVmEOjpiHtIwj7YjzE244a2riyFxuWK4f4t54/in9735P4y0f2nNcJ6TKzF8E1AAaYeT8AENFDALYAeMUqswXAg6yyNT1FRIuJaCWAtQ18tobXhibw/k8/dqb30lSKVQ+HT82gryuHF4OdoN504SKzkYW2axYyrpltZFyKZDoEQo0g55Jx2mn1V31Gq/0FFLIuFnVkjfDIZVR20r7OnJn96YRz+ayLFYvCAaoQaAp21BARoaeQiQyghay6ph015DqEnnwmnNn2qHK5jGNmhT0FZTfuzmcsU00euYyDRR3ZiEZQCDSTMHIlaxYVLbV8BD15taOaEQg9BXRkXXRbdenvzqO3kEEu40QGy66ci85cJjKoLurIIuc6EYGgzWC2UNDmEXsm35NX+ZfswXdZVz4oF16vp5BBR86NOIuXdufhkPIJ6Ost7swin3Ej97G0S5Wz72NJZxZZ17GEXc5ojMst81NflzKLLbcEwdJuVc7WNBd3ZkFE5tn3dYWOY1vQLx0rRrSDnkIY968H8N5CFku7ldCwfQS9MR9BVy4T1IXMvdk+goKlEWghpNuu2qoyNJHpV93PQrNYwfIVhVq23acAvYAzGrzhOrWp4m3T0JtXLULGdXDlyl489toIlnRm8cyBU/jy04eQz7io+oycq8xVyiwbRgG2it9917o39PlGBMEqAIet/wcBXNtAmVUNfhYAQERbobQJ9F64PuLEahc+ev0luOkXLsJ3XjqG4fESFnVk8Y5Ll2Hr9evx9kuWAgD+869cZha6bLnqQuND2LiyFx+9fj2uXd8HAPjY+y4zKzlv+oWLzSKYN6/qxdbr1+NdG5YBAD65+QpcEqw1+NdvW401fZ1wHMKmtUuw9fr1uGad+t6Pv/8yU+63rr3YhCNeddFibL1+Pa4L6vdff/UKXBY825uvudg4Ua++eAl+913rcN16Ve4Tm6/A5YFj8reuvRhvCaJ4fn6NKndt8L2fuOEKcw1djohw7fq+SLlP3nilcZzefO3FJipIl7tu/VIQEf7gxiuNw/ZD163B29YsARHhFy9Zqq63vg9EhD/8wJXGofzvrrsYP79mCQCE5dap633qA1eakNoPX7cGv7BWPf93blhmyjmOut5VF6tr3PL2Nbgu+J3ecakqd826PjiOut7Vptxa/OIlyyLX27R2CVyH8Icf2IhNa8PrvTP4Pd91mSp39RpV7lMf2Ihrgjrd8va1uH5Dvyq3oV+Vu3gJMq6DT914pfltPvKLa/HuK1S5X7qsH1uvX4+rLlqMbEK591yxHABwfVDu51arQe4PP3ClabO3vmOtSfD2y5ctx0evX483XxiW0/d4y9vX4L1Xquu95wpV7k0XLjLfq+/xlrevwdGxFQCAd1+uyl25steUu/4yVfcPv30N3h+Ue9+VK/DRX1qPy1Z0w3UIn9x8BX75cvVdv/PO9SZJ4TsvXYbbfukSvOlCdb1Pbr4C7w7u8cPXrTV5m9535Qrc9kuXYP2yLmSCcrruH75uLU4GWWzftaEf/+cvX4JL+lW5//qrl+P9G1Wd/tN7N+CVo+O49R1r8eLhMXx9x2FkXILrOKh6Piqe3zYJ6uyQ4LOBZku5SkT/BsCvMvPvBP9/GMA1zPwfrTLfAfBnzPzT4P8fAvh9AOtn+2wSmzZt4h07dpz9XQmCICwwiOhZZt50Np9tRCMYBHCR9f9qAEcbLJNr4LOCIAhCC2kkamg7gA1EtI6IcgBuArAtVmYbgFuC6KHrAIwx87EGPysIgiC0kFk1AmauEtEdAB4B4AJ4gJl3EdFtwfl7ATwM4EYAAwCmAdxa77NNuRNBEAThrJjVR9AKxEcgCIJwZrwRH4GsLBYEQVjgiCAQBEFY4IggEARBWOCIIBAEQVjgtKWzmIgmAOxpdT1mYRmAE62uRANIPeeO+VBHYH7Ucz7UEZhf9exi5v6z+XAjC8pawZ6z9X6fK4hoR7vXEZB6ziXzoY7A/KjnfKgjMO/qufZsPy+mIUEQhAWOCAJBEIQFTrsKgvtbXYEGmA91BKSec8l8qCMwP+o5H+oILJB6tqWzWBAEQTh3tKtGIAiCIJwjRBAIgiAscNpKELTrRvdEdBER/ZiIdhPRLiL6veD4HxPRESJ6Ifi7scX1PEhELwV12REc6yOiHxDR3uB1SYvreLn1vF4gonEi+lg7PEsieoCIhonoZetY6vMjok8GbXUPEf1qC+v4/xDRq0S0k4i+SUSLg+NriWjGeqb3nos61qln6m/cimdZp55ftep4kIheCI635HnWGX/mrm0yc1v8QaWp3ge1q1kOwIsANra6XkHdVgK4OnjfA+A1ABsB/DGA/9Lq+ln1PAhgWezYXwC4M3h/J4A/b3U9Y7/5cQBr2uFZArgewNUAXp7t+QW//4sA8gDWBW3XbVEdfwVAJnj/51Yd19rl2uBZJv7GrXqWafWMnf8rAP+tlc+zzvgzZ22znTSCaxBsdM/MZQB6o/uWw8zHmPm54P0EgN1Q+zHPB7YA+GLw/osAfrN1VanhvQD2MfOhVlcEAJj5cQCnYofTnt8WAA8xc4mZD0DtxXFNK+rIzN9n5mrw71NQOwG2lJRnmUZLniVQv55ERAD+LYC/Pxd1SaPO+DNnbbOdBEHSRvdtN9gS0VoAbwPwdHDojkAlf6DVZhcADOD7RPQsEW0Njq1gtVscgtflLatdLTch2sna6Vlq0p5fu7bX/wDgu9b/64joeSJ6jIje1apKWST9xu36LN8FYIiZ91rHWvo8Y+PPnLXNdhIElHCsrWJbiagbwDcAfIyZxwHcA+ASAFcBOAalRraSdzDz1QA2A7idiK5vcX1SIbV16W8A+HpwqN2e5Wy0XXslok8BqAL4cnDoGICLmfltAP4zgK8QUW+r6of037jtnmXAzYhOVFr6PBPGn9SiCcfqPs92EgSDaOON7okoC/UjfJmZ/xEAmHmImT1m9gH8Lc6ROpsGMx8NXocBfDOozxARrQSA4HW4dTWMsBnAc8w8BLTfs7RIe35t1V6J6CMAfg3AhzgwFAemgZPB+2ehbMWXtaqOdX7jtnqWAEBEGQD/GsBX9bFWPs+k8Qdz2DbbSRC07Ub3ga3wcwB2M/OnreMrrWL/CsDL8c+eK4ioi4h69HsoB+LLUM/wI0GxjwD4p9bUsIbIbKudnmWMtOe3DcBNRJQnonUANgB4pgX1AxHdAOATAH6Dmaet4/1E5Abv1wd13N+KOgZ1SPuN2+ZZWrwPwKvMPKgPtOp5po0/mMu2ea494LN4x2+E8ojvA/CpVtfHqtc7oVSrnQBeCP5uBPAlAC8Fx7cBWNnCOq6HihR4EcAu/fwALAXwQwB7g9e+NnienQBOAlhkHWv5s4QSTMcAVKBmVb9d7/kB+FTQVvcA2NzCOg5A2YR127w3KPvBoC28COA5AL/e4meZ+hu34lmm1TM4/gUAt8XKtuR51hl/5qxtSooJQRCEBU47mYYEQRCEFiCCQBAEYYEjgkAQBGGBI4JAEARhgSOCQBAEYYEjgkAQBGGBI4JAEARhgfP/A1luR7Dshv8LAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#Plot\n",
    "fig = figure()\n",
    "ax = fig.add_subplot(111)\n",
    "\n",
    "plot(arange(P), prob)\n",
    "xlim(0, P)\n",
    "\n",
    "show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c1a8729b",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}