{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "c4c9806e-bf22-499b-8035-ecd800be5913",
   "metadata": {},
   "source": [
    "# Quantum Noise Simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e9e19b48-5c1b-4867-92f2-5b88144a084f",
   "metadata": {},
   "outputs": [],
   "source": [
    "from importlib.metadata import version\n",
    "\n",
    "print(\"qiskit-aer:\", version(\"qiskit-aer\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "203f2040-56b1-4213-a041-c975dece7dc9",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from qiskit import QuantumCircuit\n",
    "from qiskit.quantum_info import Kraus, SuperOp\n",
    "from qiskit.visualization import plot_histogram\n",
    "from qiskit.transpiler import generate_preset_pass_manager\n",
    "from qiskit_aer import AerSimulator\n",
    " \n",
    "# Import from Qiskit Aer noise module\n",
    "from qiskit_aer.noise import (\n",
    "    NoiseModel,\n",
    "    QuantumError,\n",
    "    ReadoutError,\n",
    "    depolarizing_error,\n",
    "    pauli_error,\n",
    "    thermal_relaxation_error,\n",
    ")\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c436c20e-6b69-43e9-b06b-c34ed6a1ca75",
   "metadata": {},
   "source": [
    "# bit-flip noise and phase-flip noise simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e5dea141-4517-490f-8983-bb5dd9df6859",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Construct a 1-qubit bit-flip and phase-flip errors\n",
    "p_error = 0.05\n",
    "bit_flip = pauli_error([(\"X\", p_error), (\"I\", 1 - p_error)])\n",
    "phase_flip = pauli_error([(\"Z\", p_error), (\"I\", 1 - p_error)])\n",
    "print(bit_flip)\n",
    "print(phase_flip)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5f80e811-2027-48d3-9d6f-8e3a4c04c221",
   "metadata": {},
   "source": [
    "# combination of quantum noise"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e2546310-22f0-4deb-8ab0-78b7cde48de6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Compose two bit-flip and phase-flip errors\n",
    "bitphase_flip = bit_flip.compose(phase_flip)\n",
    "print(bitphase_flip)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "23afc525-6b45-489e-bfd8-d7717fa1d9b1",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Tensor product two bit-flip and phase-flip errors with\n",
    "# bit-flip on qubit-0, phase-flip on qubit-1\n",
    "error2 = phase_flip.tensor(bit_flip)\n",
    "print(error2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c316fb18-7fad-4453-bbce-d06fbdbedf20",
   "metadata": {},
   "source": [
    "# readout error simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f546e527-d645-4611-9c11-aee3747c761b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Measurement misassignment probabilities\n",
    "p0given1 = 0.1\n",
    "p1given0 = 0.05\n",
    " \n",
    "readout_error = ReadoutError([[1 - p1given0, p1given0], [p0given1, 1 - p0given1]])\n",
    "print(readout_error)"
   ]
  },
  {
   "attachments": {
    "f74cf98a-ca61-4db6-aa26-859270981222.png": {
     "image/png": "iVBORw0KGgoAAAANSUhEUgAAAcoAAAFkCAYAAACgintCAAAAAXNSR0IArs4c6QAAAHhlWElmTU0AKgAAAAgABAEaAAUAAAABAAAAPgEbAAUAAAABAAAARgEoAAMAAAABAAIAAIdpAAQAAAABAAAATgAAAAAAAACQAAAAAQAAAJAAAAABAAOgAQADAAAAAQABAACgAgAEAAAAAQAAAcqgAwAEAAAAAQAAAWQAAAAADDMCIAAAAAlwSFlzAAAWJQAAFiUBSVIk8AAAQABJREFUeAHt3Qm4NUV5J/A20ahjVBT3hcWogCAQQEEMmyCyiOyrAmERDW4TTcxkHmd0ssw8MyYmKosouCHKIsgmi7KjRkVkU5FFBBVJENTEJG5R5/4K3i/19XfOvefee/bz1vP0qerq6uqqf/fpf79vvfVW02RIBBKBRCARSAQSgUQgEUgEEoFEIBFIBBKBRCARSAQSgUQgEUgEEoFEIBFIBBKBRGB4CDxkeJda6Uqjuu5KjcidRCARSAQSgYlC4DejaO1vjeKiec1EIBFIBBKBRGBSEBi0ZDfo+icF52xnIpAIJAKJwOAQGKikmRLl4G5c1pwIJAKJQCIwBQj0W+LrV339qmcKblF2IRFIBBKBmUGgX5Jhv+opwD90jOBvk2N7f4yamk1JBBKBRCAR6DMCyK1+7/eV7JbT1n4RZd25XtrTqXw7r73fS71ZJhFIBBKBRGB6EeiVPIM/ei0/L2JR2byFejjYaz11uTodl2jntfejXMaJQCKQCCQC04NAm9Da+3pa59Xp+VDotdx8dawk5s5bsHWwVwLrVC7yIm5VvUqbupVrn5f7iUAikAgkApOHQJvM2vvRo8iPOPLFnfLq45HutVyUL3G/VK8rVfrgTk1wvaSjjrqsvPZ+lMs4EUgEEoFEYDYQQHCduCCIz7FI9x2RThee7yK9lo9yEauzl/R85drtqutrH8v9RCARSAQSgfFCYD4iq4/VaT2o9xdK18fn632v5UodiyWbXspHmXbsgpFXpxfKKw1tnVufH8czTgQSgUQgERhfBNrk1N7X8siLeKG8Tsfrcx3vFHops+K8mqRWZM6TWKh8HO8Wq7rTsXZeXa6d7rQvrw5RX52X6UQgEUgEEoHBIrAQAbWP1/ud0vIivx3rSTuvvd+tt1Gu2/GV8hdLKPOVj2O9xHWZTmmNrPOj0ZFXH49j88X1efOVy2OJQCKQCCQCvSOwGMKpy3ZKt/NiX9wtraX1sU778tohzmnnd9zvlUAWKhfHxZF2wXo/8sXhOq/O61Y2ykR9YqHOfyCn++9iynavJY8kAolAIpAI1AgshnDqsu107NexdL3vur/ukFeXU6bej/PldwoLHS/n9MPqNUhIXKddIPLq/Dov0guVdVyo63kgZ+HfOGfhklkiEUgEEoFEYLEI9EQ2D1YaZSOWLR3bQvuErCjr3V7XE+fGO9+xTmWUW1RYiCjjgt0qjePi2KJs7C8lVkd9XuzXsbSgXLcw37Fu52R+IpAIJAKJwOIQaBNWfXb7WOzXsXSn/cjvFrev037nO09e1F2Xl47y3Y6X8gsRZbvSej8uII6045GO/Ih9CUQ64k55cayO2/XW+9KC8p1Ct/xOZTMvEUgEEoFEYHEIdCOZdn7sd4rlzbeFyrUuI69TiPrj3W9fOvI7nTNv3lKJMhoQsYu00/bbRFjvRzpi5WPrlFdfI64VsWNCe/+B3PxNBBKBRCARGCYCbVKK/U6xvNi6EaJ3e31MGk/U583triBD5eNakY5YuUWFpRCliwntOPLk22qyi3Qd1+lu5aOuOo7r1LG0oFyGRCARSAQSgdEiECQVrYj9OpZub20ydNx7PeI4rl559gV8Emn7jgn1uVFHHHugRA+/3YiyG+FEfqdYXmxBgvYjLY70bz+Yjv36WJ0X9dXx3KkrCDHy5dVBfoZEIBFIBBKB4SLQiYTkRX4dR37EQYJ1LI0TxPFejzrqY3OHV5STjqCs8zrFUUbcrrs+1nQjypUKddmJisXdtpr0ggzFNVHW+XV5dbb3NaV9XXlC5D+wt+p+5GecCCQCiUAi0H8EgsCi5npfOvbruCZF+fW+dGzBMeqWF0EaT0RemweUU6/8uK68RYWlEGXdkHbafmw1yUnbgiDbcRyvY/XEftQZ8dyhFdeR1yl0y+9UNvMSgUQgEUgElodANyKSH5srRDriIEP7kW7Hv5o7Fu/0iNUlKIsrxPUx9QnyOqXLwV5+FkOUdQPa6dgXa7A40kF2yLEmyDqtTOxH+agn9qPOOp47bSVgHKtDe78+lulEIBFIBBKB/iAQRBS11fuRFre3IET5kY4YOUoHSYqFTu915eTjC2n1RTlpwX6ndDk430+vRBkXjFidkRa3tyC3iDuRZOQFQca+c6TVGeeL62vYFyLvgb0HfuVlSAQSgUQgERgNAkFGcXX7kRckFnn2Y5PXJkfHvNODJOeS84YoX/OAeu1HGyId8bwVOtgrUdYV1Q2ItLjbVpMdAnTNIMVOsfLy6/Oko/46XV9/rsiKEPkrMjKRCCQCiUAiMHAEgoziQrEvjg2Z1ekgSrH3e8TSIU3G+79+t0cd4vp4pOeyy3XsKyPU6QdyevhdDFG6QIR22n5708lORBh5NWHW6Tju/HpTv/06nttdcV3pCMrUob1fH8t0IpAIJAKJwNIQCAKKszvtR544SLKOpWNDjLZ493t3xzaXLCHqi/2I5UfZeOeLo3y3dJzfNe5ElHGBTifVx6Q77cuLTtZxEKD4ob/5zW9u6HSBzEsEEoFEIBFIBNoIPOQhD9nowbwgPnGnLbgpygVP1fuRXuUyD9a5Un6bKKPClQq1dtplolERh9RnP4gSOUrXZNmqNncTgUQgEUgEEoGuCOCPmhjrfRJpcI4ywUfimhTb+3OHVwmrlGkT5SpndMmoG9EpHQTpWKRrkpTOkAgkAolAIpAI9IpATYxBmPgFSYojrxMnxTVq0oy8BWOVLxRctFvo1KA6L0iyjnV2qQTdrR2ZnwgkAolAIjDdCOAN/FHzSaRr3umU7oaMsguGXgmrl8rajYsORKyD9bZg47JAIpAIJAKJQCLwIAL4I4x+SIYhSQbHhPo1uGgh4JTrScJ0gV5DmyyjMRGrJ9J1HJ0QJ1H2inaWSwQSgUQgEagRqPmj5pWabyLtvEhHXNclr+ewGKJUaVywfZE6v06r337dKWkdzpAIJAKJQCKQCPSKAN5oc0nwS807dbquu1t+XaZj2kWXGuqLttPRGfnScVw6OrvU6+Z5iUAikAgkArOHQM0fNacEz9RcE8freMmIqXi5QUOEOo7GRcPF9ZYSZYEsfxKBRCARSAR6RCCErJpLcE3wTPCO6mo+6rH67sVcYL5QXyzSyke6Uyyv3upO1en5rpvHEoFEIBFIBBKBGoGaP+p0zTeRdp600I4jr1N+OaH942ILhahMOel6vz43jnWK605Fuj4304lAIpAIJAKJwHwIBHfUcSe+ibxOdbWP2V8wuOByQ/vC6os8sWtEXHdQuQyJQCKQCCQCiUAvCNT8UfNKzTd1PZFf5y0p7WJLDRpRh2hU5Me+MpGOeDnXra+Z6UQgEUgEEoHZQCDIMXhELMR+nY5jpcCDZSK96LifhBWN7TVedGPzhEQgEUgEEoGZRaBXbolyfQOqn0TZrVHR6IhdcxjX7daezE8EEoFEIBGYPASCO4JLIh54T5ZLWNFQcR0Wm1+fm+lEIBFIBBKBRKCNwGJ5pVv5dr0L7i+XKNsXiIbV+ZEXcX0s04lAIpAIJAKJwGIQCC6JuD63U159fEnpfhNlNCIaKxZiv50uB/MnEUgEEoFEIBFYAIFuPBL5wTcLVLP4w0slymhYL1eMxsc5EfdybpZJBBKBRCARSAQgENwRceT1gk59Ti/lVyqzVKKMSlx8MSEau9jzFnONLJsIJAKJQCIwfQgslT+WzTfLJcqFbkU0MOKFyufxRCARSAQSgUSgFwSCVyLu5ZwllekHUWpkbJ0aEZ2IuFOZzEsEEoFEIBFIBHpFIPgk4vZ58mNrH1v0/kMXfcb8J/StYfNfJo8mAv1B4N///d+bW2+9tfnt3/7tZs0112we85jH9KfirCURSARGgcBAOKjfRBnALNRYxzMkAiND4Ec/+lHz/e9/v/n617/efOUrX2nWXXfdQpJJlCO7JXnhRGAhBObjjYU4Z6G65z0+KKLsdtH5OtrtnMxPBPqGwG9+85vm5z//eXPLLbc0J598cvOpT32qQZqvf/3rm80226xv18mKEoFEYGAIDJ1HhkWUNdsPvZMDu11Z8cQh8Otf/7o56aSTmp/85Cdle+hDH9r81m/1Y6h+4qDIBicCk4hA8EfNKQPvx6CJcqidGThaeYGJRwAp7rrrrg3CvOKKK5rvfve7zX333df86le/KuOUE9/B7EAiMFsIDIVjBk2U3W5ZfBV0O575icBAEHjIQx7SrLXWWqXu73znO81jH/vYkkacKVkOBPKsNBHoFwIj441REWW/gMt6EoElI2C8EnEOM/ziF79oWNrajJX+9Kc/bf7Lf/kvzWqrrdZQA//gBz8oxxD3Ix/5yObxj3982YbZxrxWIpAIrIxAEuXKeOTeDCGAKCPU6cgbRPzjH/+4+da3vlWMie66667m3nvvbZ71rGc1W221VSHKiy++uPn2t79dyPIZz3hG8wd/8AfNy172sqIWHobEC4fAYhjXGwTGWWci0G8Ekij7jWjWlwjMgwAjIuOiX/ziF5sPfvCDhRxf8pKXNDfffHPzzW9+s1l//fULed54441l+sq5557b/PCHP2x22mmn5klPetI8NffnEIn3X/7lX5r/+I//aJ761Kf2p9KsJRGYcASSKCf8Bmbzl44AtWtIT8NSwSKfzTffvPm3f/u35oQTTmh++ctfNiTL3/u932v+/M//vIyfUski0U9+8pPN3Xff3ZxyyinNi170ooERJUMmKuCPfvSjpS1UwLZ/+qd/ao466qhmgw02aH7nd35n6UDnmYnAhCOQRDnhNzCbvzwEFiJIpGWe5T//8z+XCyGVUE8ik27nP+xhDyvE9ru/+7srNdB4pK0m6DXWWKOoV7fffvvm4Q9/eCn/5S9/ubnqqquKivZrX/taGc9cqaI+7ujbmWee2Xz2s58tauDf//3fL6rfK6+8sjnttNOKhXDOMe0j4FnVxCGQRDlxtywb3C8EkF6Mw4kZ0LTD/fff3/zDP/xDc8011xRy4+oOOdqoJxFekGWkxdSkpqHw+NMOPALddNNNJVu5V7ziFc0WW2yxgiSj/CMe8YhSt2sGscaxfsUk2n/8x39s3vve9zZrr712s9122zUvf/nLy5SZO++8s/nMZz7TPOEJT2ie/exnF4Ojfl0360kEJgmBJMpJulvZ1r4iEHMnER0yst8OrFDNt/zwhz+8giCVQaxRPghSfhDa05/+9OY5z3lOR6K85557mm984xuKF5Xqeuutt5JlK8I2TsjwB1ky9tG+TkFZhG0TqEj1Rzu0MT4EOp0rj7SsLYgbsZNuBVKvcdELL7ywHOcP94UvfGE5lj+JwKwhkEQ5a3c8+7soBIwdvu1tb2te85rXrDgPCQURrchsJUL12souu8b+QqJk1YpU68Dgh2WsscsnPvGJzY477ti0VbhRHkGaD8pSFilSm4oRKDXuox71qCjaMeZs4dprry3HSJSuJzj3+c9/fjn/9ttvLz5xkygLNPkzgwgkUc7gTc8uP4BASJFIL6TLNjYkK6TRdpbunFC5ts+x7xiybAcSKtUrla6AjJ7ylKesVIy6MyTOZz7zmc1+++1X1J9RSFupTM8+++yiEkaMG220USE1KtTbbrutkCUDoD/6oz+K0zrG5nMiZAGpBrGq89GPfnQx6jGFheSZIRGYVQSSKGf1zme/izozxiXFyKEd5Nk6kV67bC/7SImEhkipOZ/2tKcVxwLODQI8//zzy1QR44J77713GR9EzALr1DvuuKM566yzmhtuuKFM4TC+yTKV1OncSy+9tHnc4x5Xxj3LSfP8mA5i+omAJMOYKD4CGCyx0KUKzpAIzCoCSZSzeuez30UqQzxCGOYMGhYkSa3K687GG29cYkQcRjWf//zni5TII88ee+zR7LPPPqVJiEsZJHnGGWc073rXu8qYIiLdeuutV5C8MU31kURNQ+klxMeC84Ig4zx5sEGoGRKBWUUgiXJW7/wM99tLn/RmrNBGWqNatJGeSFUkqUEEY4l3zlmTuoYxQeRnTJLE9qUvfal561vfWo4deeSRzYEHHrjCuEZblLnooouad7/73aV95jhuueWWhRiRnY3alToVUW644YYLdoGkTMUqmAqDFOs5k/apqBFwhkRgVhEYzNtgVtHMfo89AsiEatI4oHmDFm4WzCNkjWouI+OZF7zgBX3vC3Kmev3e975XyIdDAYRMauOZx/W5qzv00EOLh55w2B4NMU2FatW46UEHHVSsYYPUkC33dwx7TDmhiu1FXUzdylWegGB9RMQ8T9Klfe0Yhleg6GfGicC4IZBEOW53JNszUAS8/DfZZJMicVFtkiBDajKuZ87goEiB2pQEizBN+fi///f/FsJEdgha/OQnP7m0DYEh0DogU3WQRrm9CyIlCVqI+m/+5m+KoRApkyu8XgIVb0ieCJxUzUE7SZfkCx9jqcZLMyQCs4rA1BOluWhM4L0AWBp6uay55prlhRhf47N682ex34gSGdmGHZAZq1dEZ47l7rvvXohTmxjidJsCEu0kNcZqI2Ln8arDQpZ0jEh/9rOfNeZlqr8O/geuryxVq2kppoL4P7CY5VrPcS7zqISRL89ACBPpPve5z62ry/SMIOA58HHnHUrjQAPiuaGF8BzPSphKovTF7kVikjTLQF/KbrIXC+nBS8Gfn9eUmGA9Kzc8+zk6BMydZGHqRbPOOusUiXExjsfXmltH05xLnnTOOeecQmbGDj3f1LbGWknFSK39IaAMtS11r6kuJNYXv/jFJa3Offfdt8zF5IFIHdpJlasuEirr3Ayzg8C//uu/FjU+94nG1X2QUcMbuvDs0MqY2uS5oKqf9jB1ROll4aYiyb/8y79svvrVrxZDBDfWV9Dll19eXiykykMOOaR59atfXY4Pynhj2h+g7N/8CHgefZz5eCPNkSg9h4jMMVtbxdqtRhP+jZ+efPLJZXrIBRdcULzlGNdk4fq+972v4ZNV/axq60AiMH/TyiXUqdddd10hQaS5+uqrN29/+9ubd7zjHeX/QbLULquZ/MVf/EWzzTbbDMy4qW5jppsiwYdBmediFFIbzQUV/8c+9rHm2GOPLR9TtA4+1Iyx0zT4mPL+POyww4oGopfx8Em+v1NHlCRJRg+vfe1riwswvisPP/zwYkLvZfWnf/qn5c/PevD4448vqiarNnhZZEgE+o0A1RV15uc+97nmC1/4QlFhGReU52ud559ev8iNE/7Zn/1Zc/TRRxcJkqMCLyxak9NPP700nWceFq/tgEhNRznggAOaI444opSJ63rJUcH6sNReH5qI8n//7/9dXtRt0m3Xnfv9Q8AH/p577lmMtUj5vY41968FTfHU9J73vKcQonv/t3/7t822225bnlMkTpvBsYUPM8RprJ1mZJqFjakjSib2H57zy0kVRcXFf6V5Zl5Ogq97nk58XV9yySVFFcWDiYVzB2XE0c+HOOuaLAQQDvLxAWdOpJjqyrOJKI0H9hqoWY2rG8v0YUfaQHKGFTz3AsMcEmU7KEeiRdh8uFKbxX/CkIQpIMjSmKT6tRuROpZheAj4mEdGxpxZXg+bKM3ztQaqjzrE94d/+IdFS4EISbjGtz3H5vt6f1pX9e/+7u+av/7rvy7P0LQ+L1NFlNRLXhiWB/JHp6qico0Xgsfdzfdlbdzl6quvLl9Ep556ahmrTKIc3gthVq7k5eL5I9EZExQ8mzYvxcW+WNSHLG3mVRo/Mq5IXSogYWSHnGu1nbElZGhMifqWJFsfLyfP/SBUW4bRIODjxEeMaUscUzDCYoU8rOBDyvvTEIGPOHN5fXh57gTvT8MG1PGeOR97tHOvfOUry7j7QgZpw+pHv68zVURpsrWxFzfZn33bOXVB248mAE0BYExB524l+csuu6wsdeQBzYnV/X7EZrs+ZOQZ7PQcLhcZWhFSImtuRjeuxeDCy9WXf02EXnTG5UmybSfsy21Hnt8/BIwZ77bbbmV8+Prrry+ENAyi9CFlY8NBneqDy/NCou2kUiWEmI/MBoSAwu8wlf60EuXKE7X6d79HUhM1QKzK4MVgvKaWJutGeVkwoxdY+PkyYuiQIRGYFAQQpA9CczD/1//6X81//+//vRCkfFaLdfCyM56ZJFmjMn5p6m/ThsRIyALeCGzQgRaCdsL70/PjWaF5C0myfX2CBu0ELQVV8XnnnVcElHa5admfKonSPDJfQ26uQWgqgzBYaN8wD0I9NcSXOYm0Pf+sfV7uJwLjgoCv+lDjhgrXvuffCyzD5CFApe69RVXPexQ1qDHlTgZa/eyd+beMIJGeQCPBKUa34P1KM4fQzVFnUGa+JRuQThJot3omJX8qJEpjPZYCYsDDsMHDRtVFjRovkPYNCYOKyI+5aLGfcSIw7giQJj3rVKziSMvvJgmMe59mvX3xvmKdzMrZWqGmAw06IEpj3aGJIGCQGKM97evLpxK2+TgTaOSQ5jSGqSFK5vaxFJAXBsOc+b6qfRGRKiMgWWObGRKBRCARGDUCm266aVGp+9iPpdOCxAbRNqpXFq+ssgXXXWjKHKmzXqfVOziJchB3p0910uET+xkyCL6ofel0+xqKMqTKCAwgzB/LkAgkAonAqBGg0jT2TAXr3fb3f//3xaI/hIF+tw9R8vgkFmgpEOF8gUBSv0OR5KDaN187hnFsKsYoqV5ZXtGPC3TkvohCJdAJSNImQo3AnD4eksjLOBFIBBKBUSHAy5IxQ7YTpErE5Z3FSJFqdD6N2WLbTMjwDo13YKjx56uHIFLbgDjX0Nc0hqlQvSJED1QQpRvogZpPonSs/aAh3KhjGm929ikRSAQmBwHjzFSwJvMjR9Mx3vzmNzcf+tCHyhqm/exJvEO9AwXXrgWJTtfyDq3Hwr074/xO5Sc5byokSjeLCiCIMSTJ+ia2b5KyUT6OUeHGuZG3lBhpU0HQ98c16rrrdsXx+rphDu5YXbZui/J1uagn4ijbrVwcb8dRp3zXbtcX5dv96Wc5dfXa7/nK+dOqa7769McfPNoPLxqJ+n4oE8ej34770Ip8ZeoQ5eT1imOndtbtiHvTqVxcW/nFlos2Rh2dYvVGW3rtz2LKxTUDz7iW/MX2Rx1RT8RRv3o7PRf19ZR1Xv1cyOulP+rp53Oh7wjrLW95SzHqMcn/uOOOKxaqVLMsn9eamw9ejxVq62KDNrPbCKvXwGi+erSrJkYSZVv4mO/8STo2FUTpoaaW8CALHq6FJEMPdPwBneNcN74fN5q/RisvUJnEH7X+k7Wv7foRei2nXm1V13z11eVco36w45p1HP13nrLq7hSinb2WUwe8F6pPuYX649q9lFsKPr32p9/lFupP4L1QOe1aSr9H+VzMd78H2e/5rlvjCPNB3O/F/B9cn2W/OY6clrMw9Z7hSecP59zM8fq0kASoH92Cj0Pv0DDG8TzYXLdbcCxUtcq4viGvaQxTQZT+TBwLhCcSN5g01+2l7EYiUuOSEXiUsMVLOPKXEnuYzSti2h31iSPtDxIk3X4Q7c9XTp/inCA1efGii2PRd/v+BHG9ulxdJtJRpzbAKPLb9SoXeVGnspEX57X7Lb+XctHeur5IRxvnw9F9i3KuF21UR7Qt7q0/eFxPe/3563J1OuoMfKKOiKPuut9x7eX0WxvUKXTqd113tFFeXNv5QrTP/rQ+F536XTr/4E+v/e70XNT1tJ+LwDgwr8sG7kt9LpzPMUoYLJrOQbpkCfvSl760vtSS0p4Z8yJZ/3t3+u/HOzT61a5YW6I9jjGgrMcs2+UneX8qiNJN5rSX6kBwk30Zxcuv0w1yg2tzazd5ueqLuA6nBQcddFBZhT4eMn+eaI/2xksvzolYmfijKRMvvTgesT+Ofqrf1ks553pJzBfUKahf2Wh/+xxtVKaXcovtdy/4LKbf8+GjX0GM0voz35f5Yvod+PTSH9futVw+F9BaOcT/ZpDPxaj+D3qqfxwQfOUrXylGN95X3BbutddeRf063zO7MlKd9/SNowNSquD9uJAFKyKtjXcIK7UVbOcrTWbu/G/NCemTFwwvEjEv0k02JzJe0J264YvMlJAIPF8sZhHdOK9T7IHj8MD1g2jipal8/Jml5QtLLRfnRVwqa/3U7ZivnNOiPeJuZN4up86oN86v9+u8Ol899f4gykX9EbtmO9T4zGq/58MHXnFvZhWfUfXbf5CUd+KJJxaSNIWNBMmv6g477FA+rlioLjfQxvF1zQ2oQNtGMzZfMJ5ZCxtrzY2VkkqnMUwFUfqT29ae83bP0QAzZ3OP3EQSACJtB3p+D2AEqyn0y02U63W6Zlwr40QgEUgEekGAE4CzzjqrrDeKyF7ykpeUZQNprRaa59hL/VEGUVppKSRC786FfF8jUu9a0qjzePKZVqfoU0GUcbOtxYfwDHpTCXznO98pN6/TA0Xfzy+s4Ktt3XXXnde3YVwj40QgEUgEhoGAD3lOyr3P9t9//0Jkpov4oJ9P27OUtgVREjSQM0GCMWK3YGiL1o5UaVxyiy22KCvTLKSZ6FbfuOdPFVFaEsZXF5+FVCW873f68jIOhyQ9EAKVgZVErLOWIRFIBBKBcUAgjGn+4A/+oNl5553L2Hm/CTL6aYzT+2+DDTYoS3wZn7QmJsmStNi+LkMiy7zR2JEijZV2W6kprjHJ8VQ4HIgbQCp83vOeV244MmQ63Ul9QJq84447ysC1L6mXvexlRW2b6tJAMuNEIBEYJQI+9Nlc7LLLLs0ee+xRLPrbZDWI9m299dZlrV4aOe9Ii0fH3Mr6eoyK7rzzzjIdhGGRNlrrdFrDVBGlOTy+vvbcc89iSGNBZtM0EGMdEOhnP/vZMq7J6Obwww8vnvrrMplOBBKBRGBUCFBh8vfaD0OdxfSBodA222xTSBpZvutd7yo+YKMOAggjyAsvvLCxrCG7kNe97nVF/bpcy9u4xjjGbSuXByZb/WdL7c+3IdpOm3rbGzWvPPFD3/GOd+wzF/c1xMNlmoebSb1q+SzqA9KiPCbWlq2hdmVe/d/+238r5tXOmVb9el9BzsoSgUSgLwhQb15//fXNOeecU949JLJaq+V9NOx3EmK2agiVKiI0XSScClAF33jjjYU8rV1pEXDLgR144IFlDuWgJd65xck/OQe8+Wv87HXarHBdb6YUdNrmshcXOhFjXUM3kgxy7ESGiJC3cRu75dge/mBa/Ig51cKpc/FAAhNq45QcCX/ta18run1TNujQechnqeUmb7fddkX3z9infkAH0qisNBFIBBKBCgHvIR/uXNKxkdhvv/2KRavhoFEGkqQxyAsuuKChlfNuZORDxRqODix6v+WWWxbp07jmMMLcR8MBc9f52dzGU4yloiKWtln6xIZMOxFqkGgn8nxgnt7ciQ+GlfanypgnekhlYY6RhU99+SBLlmNI0lxJFloveMELiqPhJMhALeNEIBEYJgJUlaZUcJRyySWXFEtTY5P8tzKQMe1iFIEV60YbbdQgQyRIirTWJInSu9V4JJJkD2J/FsJUSpSzcOOyj4lAIjAdCFC9HnPMMWW9SQaJc8NSzcYbb1xUoMh00CrNSUExJcpJuVPZzkQgEUgE+oyAqR882nz4wx9uTj311GbfffctGzeYJLdpncTfZxgHWl1KlAOAl/EQNS8rXGOjGRKBRCARmA+BmOB/5ZVXNh/72MeK4SHyXGeddYr6k4QZdhbGMBnbOMeYofmP0zyHMXAbpUSZRBl3oY+xqSfGHIw3mK5io8vP8dA+gpxVJQJThgDy84F99dVXN1/4whfK2CDvN8YqjWUyqGGdzzJVWWOGLFJNh2OYOO2S5yiJcjSjxVP2gLe7w3CIARHLMUttcXpghXLu9Xz5TfN8ozYWuZ8IJAK9IeBDOqZcMOj58pe/XAjTVDZ+VXnKYSmLICOQOhnebLXVVpGV8QAQSKIcAKi77rprkSZPOeWU4vGfOsVE3le96lXFnNpXoSkpc19IA7h6VpkIJAKTjoDVkGx8vJIyr7vuuuZTn/pUmU7Ch7V3R2isNttss5mxPh3VfW2/qTvty2tvYz2PclRgxnWpRcxFMp/TV+F73/veImE6zhctwuQNyIOeZBmoZZwIJAJtBLxHaKXe+c53Nj64w50cdSwXd6effnqxkO208EO7rknfH6XqtRMx1ni2CTL2kyhrlLqkjR8gSx6CzOc0uZj/RCoWztqNLbz4xS8u85XSBLwLiJmdCMwgAgwCrQ3JVZyYhxzvDARp3JIhzx//8R8XJwWMfGbhg3uURJmq1wH+CY1FGkOweZg5QPj85z9fpEwGP9bMtMIJsrR8zpprrjnA1mTViUAiMO4IcGtngj+DHtooTsmR4N57712minh/eG/wsUot690yCyQ56vuWRDmkO8CEe6eddipqEgukcsZucN6irMjSMZ77ecPw8I/ajdWQYMnLJAIzj4AhGOs7MvpDkhdddFFxPuAdYHoIq3nu7Xx4U7/+/Oc/L8aBOfVseI9OEuXwsC5XQpBUrrazzz67+ehHP1q+Ho1DfPzjH2+OPvroYurN+o1rq1TJDvkG5eUSgSEhgCCtxmHc0XSQ0047rahaqVVZsnI48IpXvGLFgvJUsMjUhzQ3nBmGh0AS5fCwXuVKPHL4WjR+aZLxmWee2fzJn/xJ+ZPwys9DRy4mvQpsmZEITAUCDHWs63jssceWd4DVjUwhO/LII8tiDT6qa80SNawlA82nHJYj8qkAug+dSKLsA4hLrcKfgDoFWZpQ/PKXv7w577zzikrWigKXXnpps+OOO5bFW6lZ0mHBUpHO8xKB8UHAdA+rGzHuM+Zo35J/22+/fZk+hiw7Db/cc889xYiHrUOqXYd7P5Moh4v3KlejWmXmbdzSUjsI80tf+lL5A910003FStZYJoMf86X8SXLwfhUYMyMRGGsEqFlJkIjR/5uqlSMB3nRojvh0taKRuZPdAgt6lvRW95iVVTu6YTHs/CTKYSPe5XrIzx/AVyW1ClJkGu6PZXUBC7xyU/WSl7ykmIn74uRLNkMikAiMLwK//vWvi6HOt7/97aIp4jQAQfq/+zimMbJRsy7ksYtfVx/WNFGpXRruPU+iHC7ePV3NuKQ/DyvYG264ofm7v/u7MoZx/PHHF9I84ogjijrWSuSPetSjVhrH6OkCWSgRSAQGigAJkpEOd3OMcBjtMd7jlYvmiNORvfbaq4w39toQdSLI1Cj1ilj/yiVR9g/LvtfE+fGGG27YIMiwijWG+ba3va35yEc+UuZR+bOZiJwhEUgExgMBhCbQCH3oQx8q45H2rQDypje9acUiCYv1puOjWCBZUsEuJIGWwvnTFwSSKPsC42AqoWYxRcRGFUvS3GGHHZrzzz+/EKevVF+rxi/32Wefor5Jdexg7kXWmggshACC/N73vleI0XCJhRHMj+S20jzpbbfdtnnmM59ZFkaorVkXqjeO0yA5j9ce3nnYM2QYDgJJlMPBedlXWW211RobjxxrrbVW87nPfa4YBdx2221lzMNKJQjT/CpWc/m1uWzIs4JEoGcEWK7S+rAp4FGHExG2BmwKwvMWklxOYOlKqjSNBCEnUS4HzcWdm0S5OLxGXpqVHAs5y3aRMi+44IJCmMY/WMkizG222aZ49PDHzDGNkd+ybMCUIsBZAImRoQ7Njuke5joaMrFaEGcBpEjOQ/oRfCQz4mP9armtDMNDIIlyeFj39UrUsb5WX/SiFxW3V+9///ubyy67rHnXu95VDH6YnFuhxNQTZVPC7Cv8WdkMI0DF+otf/KL4XCVBchjgI5XDcutIHnjggc3BBx/cd6MbU8PYI/DOk//n4T6AuXrIcPHu+9X8aX3ZGtyn+jn11FOLlMlx8hOf+MTmsMMOK+OX6cmj79BnhTOKgDFCXrQ+8YlPFGt0+8YguZyjzQlnAf22TvVfZ0lrGT9TyWaNLOfwPGDukfvZ3Pbzue0XVSxt++WD23/MxbZftbZfz+3bWFt12uayV4QHLLIe3E2iXIHL5Cf4iPz+979fJMzLL7+8Offcc8sfyiRmf2BTTqxSkiERSAQWhwByMg5Ja2O4gwTJOTmPOlSsDHZIfDQ4s0Zgi0Ny6aWTKJeOXZ7ZAQFfuMYqef+gGvrGN75RFok2rsmwwBgngx9fpRkSgURgfgSMO1533XVlHNL/iepz3XXXLf8jQx8M6HjKyQUM5sdxuUeTKJeLYJ6/CgLUNFQ/sWD0VVddVSzl5O+yyy5FwvQVHKuUrFJBZiQCM4yAMUhTMO68887ywWn9WBbmLM8Z0ZEiuZzjUSfDcBBIohwOzjN7FWoj/mJPPvnk5pRTTinjHOZk7rbbbmUMk9k5gx/Wesg1QyIwqwgY6zeh//7772+uuOKK5oQTTihWrQx1LH312te+tli05nzl4T8hSZTDx3ymrkiKNJ5i9XTm7CeeeGJZpcR4pjGV3XffvXn1q19dXGt5IWRIBGYVgbvvvrtYjVsb1pCFj0yGOlzOUbP6oLSlmnX4T0gS5fAxn8krctDsi/k73/lOeQkw+GGcEOvgmfNlygkT95QuZ/IRmclOU7OSIHm8uuSSS4oxHIJcf/31mz322KPEa6yxRrFmTYIc3SOSRDk67Gf2ypw1W8XAGCajH0v/ULvyR2nshcGPKSXhX3JmgcqOTy0CyPCuu+4qnnQY6fB2ZfoFK3FrxG6++ebFSpzWxX8jw2gRSKIcLf4zfXUSJqu+s846q7n66quL8YKXAunSxik7V1lWPciQCEw6AoYhzDs27OC5R5AWSOeXlecb6lWqVjGCzDA+CCRRjs+9mOmWMPg5/fTTm9NOO61YyD7+8Y8vBj+8/LD0Y8CQLvFm+hGZ6M6TIGlSjNNzFsA5x80331xUqoYb3vCGN5TnnJvIDOOHQBLl+N2TmWyRL22rsJtYzZjBi4Rq6qlPfWqRLl/3utc1z3ve84oxw0wClJ2eaATuuOOOxsLJlr7iVJxRjtV4GOqYX2xlDsZsNCoZxg+BJMrxuycz26KwkGX9ZwyTifyVV15ZVLIMGiwmbR6m1dlTHTuzj8lEdfzWW28t3nQ8yxxxkCoZrRla8BybHmV+ZIbxRiCJcrzvz8y2joSJLK2MwOCH0wJhnXXWabbaaqti7MDbT6qqZvYRGduOs2S980FnAWGsZn7keuutVzzpkCA53EiCHNtbuErDRkmUOWlulduRGYEANRRXXWuuuWaxAmQub0qJ+WUMIay5t+eeexZ/l8pw45Vqq0Av42EjQBuCIGlDSI6M00z5MH+Y1LjzzjsXZwEsWnP607DvzmRfr62M77Qvr7391lye7bdbG+K1PezB7Xfm4tge/mBa/Ii5h/rUuTjDhCFg7b1PfvKTxeBHmtS51157Na985SvL1BL+Y439zCJhmqcKD5bEYsYjHGTbvJiN/8oPbOQzkMq5ecv7EyBImCNEm/H1M844ozEmaXoTD1SWnKNmNQ6ZYTIRGKVE2YkYaxTbBBn7SZQ1SjOU9sVOhcXgx0K1xx13XPGJ6YVPnXXkkUc222+//QoymCFoGqu33HPPPUVdLTYF4dGPfnSxpDQ31WoupiF4mVuKyRQE1papul7eU+IDxRABS9aPfOQjxYEGlarFk1/zmtcUjYgPOCSZHyXLw3qUZydRjhL9vPaSEECW1rw0DkQdy6MJS0JOoi3lxaMJ5wXWxJyVcO+99zY33HBD85WvfKW5+OKLC2laesnEdS9yhMgLktUngkTf/OY3N9ttt13zzGc+c1Zg6ls/EeT1119fvEvxMEWCNH2Jswwfa54/uCLJkOL7dvGsaOgIjJIoc4xy6Ld7Oi5IgjQuaTN9xFimyds8/DDBR6CWHyI1mYM5C0YTXsjGwvjVffe7310+JKhazUdlAGV1eupWRlEsMb3YqbGR6aCJUpti7I7k+6QnPam0x+oxkxZYrcLO82azBJb+IcYwMvM8kuYzJAL9QCCJsh8ozngdFq9FAiQnyxFdeOGFZSK3ZYluvPHGIkGxMOQajMGPr/5pDCRGHw4kHS9pEqa+8vjylre8pYxHkmwQ6pe//OVyHD7U2IMMxkYZXiEUS0cha5Iv9TCVsLVJJ0HiQvC0GIzJfGww1DEO7CODFBnTlnzEZUgE+olAEmU/0ZzhurxoEcJRRx3V7Lfffs3HPvaxYlBhHNNmvOiQQw4pZEq6nFarQ4tmk3KM5RoTI+H80R/90UpOGuSHhE2F7WU/yMBC+b3vfW+RKKl6Ecrb3/72MvGeg/z/+l//69gauTDUIZWTGHnRIYGfffbZRWNhqTjP2oEHHli0FoPEMOuebQSSKGf7/g+k9xwRHHHEEcUcn4RpHUxSpgnf22yzTXPQQQeVpb2mcUkvUg+ipB6kUqV69kKvA3L84Q9/WLKoQAfteN5HCw9LsGcBSqI87LDDiqR/0UUXFYMik+/HUapEkoj+fe97X6Otd86p9KmLX//615c+kCbTGKp+ujI9CASSKAeB6ozXybLQFBEvMSvBc3t3zTXXNJ/+9KeL+s9LW9oxKjNkMS0BUVqFQsz7C5VzbWlpGoOxQi9/wZQFBlCDCCGJUbkiRypysUBFjHCofhkeIdFxIkoSJEMdzwk1K4JkGMaKFa6eKc8X6Xyc2j2I+5h1jh6BJMrR34OpbQH1qpVHEIFxMKpZ0pbxORImwrz22mvLtBKEMomGJfXNI0UiQfNLjU3q01prrVUXKeNr+vzP//zPhbS89BkAtYNxTmOIxuOMYZI6fXwYf6PWtVHrhgq3fb59kiujIecjFUZXEdRF0kWmCGlcgn4hb8u/+eCgbkWaLIONp5qCxDAqxyHH5Y7NRjuSKGfjPo+0lyQqJMiLDwmKS7HzzjuvzCk0942hiXEzcwpJOpY3CslnpA1f5MUZyvAIY1wNgTFeYvEaAXFZ0ollMGMfalkv/7qMstSyCA6BkTxJTeZdIjXGLMYVESl3bL0QJfKBaVvFS2XpYwYhGyclmY1COoMXCdyUGcR4wQUXlGcEQSJFRmK77757+dhKNWs8TRkPE4EkymGiPePX8hJea07CQoZI05JeVnJAHKQsC0Ube0KmyIPUM0kWsuaRetELVIP6gPARASIyl5KUhPz09f/8n/9TJL0Yq1XGBo8TTjihSIL77LNP8yd/8ielTipSDup9aJhygwDnC9S8999/f6mTBBbXiXPs25AUEmeNO0y84eK6rHJ9YHgeeNXRblqIAw44oBjqcJ2YIREYJQJJlKNEf0avjTBJSVzf7bjjjsVhgRckEiBdkrIOPfTQImXOJzGNG3ykvW9+85ulWQx5QlL04r/vvvuat73tbcUZATUi8iMR1uOX1kk0HveOd7yjGKiwEuYaEKHATD2IlLEUK2JTbRYK6ndubJ3Kq1/d4mEGfeGp6AMf+EBxWEEdTeq1LqRng0o6VazDvCN5rW4IJFF2QybzB4qAFzj1IynGmoCkzJe//OUrnBWQtkhQVpt3nOqwLRENtIFLqJzKNIiS8ZLxWeNtVLKkRNLlG9/4xrJUGYmz7g+piqT5F3/xF0WSfPWrX12shoMMqVrNS1W/jwyYLPQRQa1KZQtrkhvVbR0QFalUO4YpTfJOdNNNN5V5kMYiSeIMdZDjy172smItTKKkUciQCIwDAkmU43AXZrgNVH2sXldfffUypodcGPtY+cFmrI7VJiMOKllSBulo3AKSNMZGKhSsUIH8jamR1oy/kiAtS6YPbRUnFSlDJ44AuADk1YjxUwSSl+XOqEmpItVVE22Uq2PSGEcQiJW0hizr8G//9m+FPBn5IPFB42rMkXs/5Eh97CPCxxKDJnjRJKSatb5DmR4XBJIox+VOzHg7EIcXtikjxt822mijMm/O2OWHP/zhIoHw38kABoGQQMbJ4IcFrw3ZkH6pjrWR1Ebyk54vsEzVV4G0GVNGSIGcqJ9zzjmFWJDehhtuWMivrs91lKM+JR0iSWURpWsjWKrhCEiLOhjZcv02KJLUHteGjfFb0z1I2z4eTFeheg8pcpzuZ+CUcSIAgSTKfA7GDgFS5b777lteoqRKXmWMXZJCkISxO1ayxgCRQT3ON6rOmBKCDJCUaSGk5FoiXKhdppaEEwLkhSBJgPKobTluYO2KRFmB1iHICD5UtPycrrHGGgUbdcGMNEp1q155pF9TWZD4tgNwNqBNyJCnIta7nE7wqqN9+kDN+qpXvWpsNQQ1vplOBJIo8xkYWwSo5RitUMtxgxcGP29961vLckp/+Id/WAgTsY46cNJ959ykeMYoVIjGBxcTECu1I+8zpC7kj+xIfabWILeYl9omSmSKSHnbIZUZ62U9zJmDwMWbtlHt+vBgTIS4qGOt8kKi67dEqT6ryvAKZGUZ1zJdxrqQrrfWnPUzqbff110M5lk2EegVgSTKXpHKckNHAFmQGG3GsTgtIEkiBWrK97znPSXtGJWtCfTDlC6pO43zGUNlrUriY5GKuBjdGJ8MY5yFwDNuuf/++5epEoxbSGPIkcccBHPaaacVHKSVrQOjF+OhVNMIkQq2drROwmRJqp0si8UcEBj3VR98+xX03RizFWQY7JBctW3nnXcufXEPqYLTmrVfiGc9w0AgiXIYKOc1lo0A600bScRGfUclSSXLWpQHGy9+0lYnY5llN6BDBYiSapEK0zVJSiRKEi7LTiTaK1EiVRIla1fkQkWK+JGK+lnF6r9+t0kG0ZmO8trXvrb4RHVO7YWHZM44iGRK/crSFXkiy35N4GeMxDjHPFFjkIyStGHXXXddcV+QZIZEYBIRSKKcxLs2w21mHWtqBPWmlzLp0vjbGWecUUgTEVA5Ig4vakZCg1LvqZc6lFRWq3/lI9HFXtf4JqtYWwQkaTkppIvseKppB310rvFGBMjgycdEHUidjGds/QrI3IcCSdpHC7VxTI/RB9I/NS+CHKak36/+ZT2JQCCQRBlIZDwxCHjpsiy1EoYXsekGVpdAnAxHSGbG5RgEkb5IVIOwqFRnSLr9Bo/UxyCGKtNYH4mSBSspk8VqW11KgoUBQyfTaLRrUIFBDuMj444k+WOOOaYQZbico0I2fkwNnQQ5qLuQ9Q4TgYe0LtZpX157+625PJsVeOsN8dosUWBj0RCbgRBp8SPmXgKnzsUZEoFlIeClzTrUhiR5eUGY9qkgrY/JwpJhzCQFfbniiitWzK1EmqxoTZ2hYqbSrAMJFpmSLpGoeFABSTLQ+chHPlJInLRrnigvQtTP4Vd2kNL8oPqW9Y4vAnMamgPmWmcy8M/nNv4bI5a2/fLBjWcNm4Ve6+3Xc/s2Lqg6bXPZK8JKbqo6EeOKknOJNkHGfhJljVKmR44AIiHhfPe73y2T2r3ITZcwZkgSYz277dw0CGN8kxD0xbif2DxE/SMdk9Kon4e9NBkiJt3C1MovHAdom48RHyIInFpX29rS7iTgnW0cfwRGSZSpeh3/5yNb2AMCxgON0dkY1hgXM3/QlAjqwTvvvLMYm5DGjOHFWoY9VD2SIoyAejUEGnQDOSpgNMWalbTOeIpDhFhP1JgwwmSlmyERmEYEUqKcxruafSoIUBGaN8hnrHFML3wSprmDXu7cwBnza1uRJnxz+qo5CdK4J0cKDHUuu+yyIkUibyufULEaH540lXbe28lFYJQSZRLl5D432fIeETAnEVGa/H7WWWcVlaGXPWMfpOlljyxzTG1uAGduzBdefNZSs77//e8v0z1I7FSrpqCwZh2ksVCPtzWLzRgCSZQzdsOzu8NFwPieqQzG1HjQsdYjTz+83nj5G2PjmxVhzrqVpnFI47snnXRSIUp3iqr6oIMOKgRp3qWPilnHabhPcF4NAkmU+RwkAkNAICxkY7zy0ksvLWRgWgODH15wzNE0B3CQVqND6OqiLkHNaj6k+ZrmpRqPRJim2bCuZahjXqZxySTIRUGbhfuIwCiJMkff+3gjs6rxRsBL3sR8lq+MeTgKQAbG4HiTsW4khwUkKGOYyHMQ8y/HBSUfDvfcc08Ze6RmZfhkHBchhkcdJMnCNglyXO5atmMUCOQY5ShQz2uODQKmXvAmE559eMIxvYH/WL5TebpBHIyApiVQRfMFy3rVRwLJ2nQPLu4Qo6k0pGsEaWwyQyIwDgiMUqJs/ws67ctrbzmPchyenGxDXxEgVVmhxOodpEsWsSbR77nnnmWyP2mUhDmp5GGclnMA8dlnn136ykG6MUdS9tFHH10sWTkMyJAIjBsCSZTjdkeyPTOJAM82iITvUoR5+umnl32SFUvPN77xjUXqmlR1LLWqj4Bjjz22OFpH+CRI7uZIz+ZB8gk7S+OzM/mgT2inkygn9MZls6cPAWpJhGkFD44KeKEx0d50CctF7bjjjmXjnJwP2UkILH0tA8ZQxwofCNMan+FcnqEOhwEZEoFxRiCJcpzvTrZtZhHgmPxrX/taWZ0kjF2MVfL4E0t6bbrppmOpijUX0ngrIyXedMwjZcnKXyySZLCkH8ZfMyQCk4DAKIkyrV4n4QnJNo4EAUtoWYmDBx9LXLGCRToMYG655Zbm61//ejGKsaQXt3n9WttxOZ1FkCxZWe9yrE59fO+99xYfrIx0eNSxiLK+Ub1mSAQSgYURaP9TOu3La29pzLMwtlliChFAjqecckpzzjnnNHfddVcZzzv88MObffbZp3n+859fxvhGMYZpqgcJmMrYOCQrXoTOgtdUj4MPPrhYsk6T9e4UPl7ZpXkQGKVE2YkY66a2CTL2kyhrlDI9MwiwGDWGefvttxey5OLNFBMqzJ122qlYyVJrDjuQIkmP3PTdeuutZToLFStLVlNcHv/4x0+0xe6w8czrjR8CSZTjd0+yRYnAvAjwZMMohtPwiy66qDgNNwbIMIYlqZU1jAFazWSQwRxQDsu55DOeKlATU7GScK2iYlpLWrIO8i5k3cNAYJREmWOUw7jDeY2pQ8B4pA0RkSbNQzR+aV4iqY6UKWy++eZ9X5/ROOT3vve94kmHNx3OAoxDcvS+1VZbFaLkcSgtWafuscsOjQiBJMoRAZ+XnR4EEJTpIi984QuLlxsSJkMZPmSpafu1kLG6ECLVKoK84IILylJYHCNwimAllI022qhIsa6fIRFIBPqDQPvf1GlfXnvLMcr+4J+1TBkC5mH+7Gc/K0Y9/eia+jgt5wiB2zlGRIyJWN0adzSv85BDDikedZIc+4F41jGuCKTqdVzvTLYrEVgkAsiKBIng+kFcrFktD3bMMccUt3O8BnE5d+CBBzZHHnlkmbIyDtNSFglTFk8EJgqBVL1O1O3Kxk4CAv1YacNUD3MhGemY7mEqymMe85hmt912K4Y6VKxWQDHdox/XmwRcs42JwKgQSKIcFfJ53USgCwKsV82BNA7JOMjUj6233rp5wQteUJb/QpIpRXYBL7MTgQEgkEQ5AFCzykRgsQgYg0SId955Z3PxxRcXo6D777+/uJzjUcd0E9M9HvvYxy626iyfCCQCy0QgiXKZAObpicByEAgHBtSs5513XnP88ccXwx1TTl7+8pcXQx1TTDIkAonA6BBIohwd9nnlRKC5+eabm0984hPFUOe73/1uMQB69atfXVzicVjAcCdDIpAIjBaBJMrR4p9Xn0EEzIc0vcPiyVb14N3Hkl38sZIin/WsZ42Nk/UZvD3Z5URgFQSSKFeBJDMSgcEgwAGBNS6RI0Odq6++uow5MtLhl5XDAu7vMiQCicB4IZBEOV73I1szZQiYB0mCZKhj5RFeeywEHWtD8snKWIfj8klZCHrKblF2JxFYEIEkygUhygKJwOIR4HDgpz/9acN5Or+sJ598cvPJT36y+IB96lOf2uy9997Nm970puIn1tqQGRKBRGB8EUiiHN97ky2bYAR45bFwMkMdDgOsNIIgrV251157NWuuuWZxFvDQh+ZfcIJvczZ9RhDIf+mM3Ojs5nAQsGqIpa9Ij9dee21z9913l3FITstf8pKXNOutt16RIi191Y9ApRu+ZU0pyQsAZzkAADW8SURBVJAIJAL9RyCJsv+YZo0ziMC//Mu/FIL8whe+UJbbuuGGG8oqHox0GOtY/sqSXP1Us7KW5ST9xhtvbKwg8uIXv7gYA1leK9efnMGHMLs8MASSKAcGbVY87Qgw1CHN8aaDGC+99NIGUfKyY33KbbbZphjq/P7v//5A/LEiZxIr4yBp61LGUltrr712kVyTMKf9Kcz+DQOBJMphoJzXmCoEGOpYPNmqHlbzOPbYY5vPfvazDZdzz372s5sjjjiizIk0DjlIh+V8vj7jGc8oEuuHP/zh5sorryxETXo94IADivP0xz3ucUWKHWQ7purmZmcSgQ4IdFp/si7WXocy9nM9yhqlTM8UAlb2+OIXv1jWhTzjjDOKZev6669fLFkZ6pDmqFiHIc2RapE2kiZRHnfccc1Xv/rVRr7x0KOOOqrZeeedmyc84QkzdY+ys9OHwCjXo0yinL7nKXs0IAT4ZaVetfSV1T0spMwoBzluscUWKwx1rEc57GBxZxLut771rUKU2og4LcO18cYbF+mSKnittdYadtPyeolAXxAYJVGm6rUvtzArmWYEfvSjHzWWvjL+yKPOrbfeWpa52mGHHYpHnRe96EVlbch+GuosFk/S6+Mf//iymYayxhprlNVGrrrqqtLuGEfdcsstS5tJmDk1ZbEoZ/lZRSAlylm989nvBRH48Y9/3HBUTjK75JJLmssvv7wQ0XOe85xipLPddtsVKXLBikZQgOqVOpa7PPM4TVUhcT7pSU8qBj8scY2nspAdhQQ8AkjykhOOwCglyiTKCX94svn9RYAKk09W45AI8qSTTmpIZfZXX3315rWvfW1RtZLYJiWwzKWK5R2IwY8PAKrigw46qNlpp50aBj8Wgmbww1FChkRgHBFIohzHu5JtmkkETLNgqHPMMceUaRf8tJrqgVQOPPDAYqDDJ+vDHvawicGHdKkf+sbfbPTNmKu5nSxkLe31yEc+cqBWuhMDWDZ0LBFIohzL25KNmhUETPcgZRl/vPDCC5trrrmmueuuu4qjcg7LTeSnbn3a055WJK5JlboQI9+zprSYe0nKtB4mIkWYlvjadtttS3oYFruz8nxlP/uDwCiJMo15+nMPs5YJRcA4nukUn/vc5wp5cGBusWRSFoK09NUzn/nMoUz1GDSEjHdWW221spl/iRz1m4qZoRJ/tBwnRL+NYWZIBBKBpmkPSHTal9fech5lPj0TiwAJ8r777itS43XXXddcfPHFxdiFSpWalU/WV7ziFcXQZdolK8t9kaAtIo0kv//97xfn7Zb/Il1aRJoBkI+HDInAKBEYpUTZiRhrLNoEGftJlDVKmZ4IBKgYTc7nYo4nHSt7MG6Rz8vNIYccUoxbSFuzGMLgx1xRFrJwMHbJYQEvQ8Ywp/3DYRbv+6T0OYlyUu5UtnOiETAOyYKV9xqSJGvQ5z//+c1hhx3W7LjjjsV7zbA86owjkCx7fURQRX/84x9vTj/99Ib0zcMPpwqvfOUry/zMcWx7tmn6EUiinP57nD0cIQLGHc0nvOiii8p43L333lu81Wy//faNCfi81ZhPOEmWrIOCEzGyjmXwY7mw888/v6hmGQIZs2TcRMI0vpl4DeouZL2dEBglUaYxT6c7knkTjwB16g9+8IPiao6F59VXX90gTMtRHXzwwcVgZdNNNy3qxYnvbB87MPcyKutnkrSf+9znFkvfDTfcsBj78E5ktRJTTMLgh7SZDtf7eAOyqrFEIMcox/K2ZKOWigDJh4HK7bff3tx0003NmWeeWaQjKtVNNtmkjEEy1uHubdAu3LSFhIZIjO1JGyNF4vYnRSKjsoalqTOsZEmbPjistUlljVAtGm0MM0MiMCgERilRJlEO6q5mvUNFAPkYXzNPkKu5U045pbniiiuKtSbJ8dBDD2123XXXMjViUA1DguHZByma5C9GipyTI05+Y42N2jdVA4FzITcJUpn+Iczjjz++ueyyy8rc08c+9rFljHfPPfcs0jnr2LSQHdQTNtv1JlHO9v3P3i8TAWTE7Rwr1o985CPFgTly2mCDDZo3vvGNjYWTuWmz0of8QYXbbrutOe+885pPfvKTRe3LchRJ77bbbsUFnjUjv/3tb5fLI0rjfO985ztLO7mQG/cAZyQP62984xvNaaedVraf/vSnzdOf/vRi8LPffvs1VLUZEoF+I5BE2W9Es76ZQIAUSQ1osjyCMobmpb3OOusUg5PNN9+8TGtAksNQc5JozUPk1YdPWKTIWYHVRfiG5d0H2ZB4OSonYZJ0TcFA5sMIMEN0PhhiW8p1f/KTn5TxSoTJSIpXI1az1LBw33333Yvxjw+CDIlAPxAYJVGmMU8/7mDWMVQEkA3DHEY6SNKEeRaaJBmkxNDEvEjLTQ0zIAXzDalT5/7UZSzSOChn6ibwk3CpJalab7zxxmIUYx7nLrvsMnCiRJCsfTkVgBl3dQxxlirJcs6w7rrrlkWq9ZnxD7J0H0j2d9xxRyFM9wJ5pjp2mE9iXqvfCCRR9hvRrG9gCJBYuJxjqGO6B486pDaGOcYfSTFezKZ6jCoYIzVHk7SIEEm35iAyIIpAyjSZn/WoMUtS8CCDDwuLTPuwIHlTT5NwTYtZKlFGe/WRetlHgKW7LrjggmJhTGom4SNOfefxSL+V9xGRIRGYJASSKCfpbs1oW0lDSJIUyYjkgx/8YCEZY44sL02EZ0wyDi9gRElio95E4Ii7Jkm3kMEPIhVIoSTQQQb4IS1TZKyv6fry+hkQoKW7WBYbqz3xxBOLlewJJ5xQ5mL6iDn88MPLBwJypvYdh/vVTwyyrulFIIlyeu/t1PSMUUx4iiGtCSa9I0gvZ9aj4/LSNXaHkBAmNfBmm222yn0wfQWZCCbxI9RBBqSErEmxvO1w3zeoYCyYqvUv//Ivy/0544wzih/ZY489tlx7//33L8uVIdRxuWeDwiLrnR4Ekiin515OVU9IPYxivNQZi1gOigrRWJ/xNao8VqOmJyCCcQjazNm6cUASpTFTBjx1INGZYmGlDmG77bbr6hZOfxEvVbP6OCc33olgSITGCXudC6qceY6DNmrSNtewGQM96qijivcjHw/nnntuc8455xQXeTQBFo32IUEazZAIjDMCSZTjfHdmtG233HJLeZkyDrGIsvE1pMPdHCnNmBiCHDeJxHjjt771rTLmSL1IskJudTC2atK+eZMMYKhmn/CEJ9RFSpoxzLXXXluMYoy5kjqNx5r8f8899xT16dFHH12kxFVOHpOMmAJjHU+Ss/FafedLlqR56623lv7TCiBVxJ8hERhHBJIox/GuzGCbSEwIERmwBCVJhss544+Wvdp4441Haqiz0G0xNYSBjsASFEHU0yNCQkYQSMO0EFJxlCFBkhRJpFbyYBlLClQW4crnSMGHBIveAw88cFFEqe7asYHrDSPow/rrr182HzlWJzHWrH/6RMLmd5elMiMjY8/j9hE0DJzyGuOLQBLl+N6bmWgZIx2T2JEkcmD8QS1JHUeCfM1rXtNss802Y6Nene+mIEpWnoKpEyQk5KSPjHe400MMj3nMY4qVLomwJivWr6TSv/qrvyrGN9aDfMc73lEIV52OsxxFuFzHLVYCq4kSEdXXVv8wgrFJc0b32GOPMo2EEwa4XHLJJcXhOlXt2muvvcLIqSb2YbQvr5EIdEIgibITKpk3NAQ42aaGO/XUUwsBeHmTIElLXqrm303Ky5ITdmOpgvFKY6cIUh5LXWN0VMh/+qd/2uywww6lXEhO+k2Cfv3rX1/I1JQSVqK16pZzBRI3icsqHgh3UgOJ24cCB/WkZN6MzL+EkTHoI488shBqjl9O6h2ernYnUU7X/ZyI3pAgSV/m9Jlvx7sLYxOGOiQNqkYv0nEch+wGMK88JvTrl4AgvfStf0mi1Jf/+T//ZzFCMl7HUrcOpOhY5WStufmNL3zhC4uqMox1HL/++uubO++8s0iVvN8sdg6kDw5SpYCYg6Trdgwj7bqmxDBMIhUfdNBBZdyZwQ/DLZoFHxfhdN249CR9MA0Dw7zGcBFIohwu3jN9NU61Gbsw0DHX0IvRC5shx1ZbbVU8uRinmkRJCUGyaNVHY3K87egHyRIZGFNk3IMgO0nIpFGY+IiAh48FkmMELu8YwiA6x7jlW2yoidK5oyLKaLfrkxif9axnlbFW47X6BgdGT2effXYx+PEBwejJs+GDI0MiMGwEkiiHjfiMXQ8RMtShYjXN4Yq5FT1sVIjGorbeeusyTYCRxyQSZNxO/TF2SPIhERpD9OLvNVgsOaRRZGB6BVI0ZqleHm8QMVUsabMOyplGEuOjlsCyIew6uBchUdb545AmYfqQIG0jRdoGRl0Mlxj9MIAylcb4pvmgS/lQGId+ZhsmE4Ekysm8bxPRaupHS01RG37sYx8rY5EIhXTF5dzrXve6ooqcBj+gpnMgNBas5gaGyrTXGwWD+FBADnBCdtInn3xycd2H6FjSchVXB6pd5f74j/+4SInGN03sRygRSLaxDJk8hkE2EvCg51ZGG3qJSb3IUl948/nUpz5VXO6dddZZRS1rbNf4NQ0EaVTbRy0Z99KvLDPZCCRRTvb9G+vWe3nzBMNIw9w/cwEPO+yw8hInQXrRLZZQxrHDCExfedsh6ZgfuVjyp340BYbkaP4kYxZTJYzT/fVf/3Wz7777FrUsq1fTZOqAAKm0qSsRjY20XhMldTeDKZKaYIksUizJty2h1nWPMm2c2ioslu7SXs8RC1kaCXMvGT55jhY7VjvKPuW1JxOBJMrJvG9j22qSC+MchizG1EhGvvgPOeSQMs2D1SepaBpebiQyKmUfA1/5yleKNIn4GSiR1FjvmibSS4AHUuSqz0cFzDgiYPBCIucSj4Rp/LKtUjUG6lzWwzZtaJdxnqk2pE1SPkJXvzrHNbAaJqGTGn0osILmhMIcTGpZzxZDH+4MEeegXQGOK07ZrsEjkEQ5eIxn4gpIwxgkQwwvM0tfedl7uVEVGncy6Twm108DKF7kDG6oCg844IBi6Uqa82I3lriYjwHnPOUpTykbUoQdsuOJx/QJUqMxT9KqY3UwvockbMjPWF/bdV4cr8+blLT+kS5JyGvNjf/qm2fMh5gVZIxf8vbDEpiE6SOgjdGk9DXbOZ4IJFGO532ZiFZRObLSNP+PRxqSFK8rXvRUiRwF8LhCtTeNY0le4KRjqsF+BgQb47vGPk2ZoCaldiU9Mo6itq6D+2CFEKRNuqR6nbbgwyQsZI1RXnjhhcVRAWMf47ikehIzP8D6jzB9uGRIBJaLQBLlchGcwfMRpJc1S0srYZx00klF1cpq08vJy8pYpBfWNIxBjuIW+9jgrch4JQLwwmfsQ1JiCUv6rAPpkxpyKeOjdT2TkPaB8vSnP72M4zIKM42E1G3Oqo81ztbNzeTZicQPmyTMSbiz49vGlXU4TdNpX15785lms2xDvSFe28Me3H5nLo7NJ7C0+BFzL9tT5+IME4iAcUjSi694BhbUgl7cPKogSFKWF5QXWqrAlnaDGecYg7PEmA8SgaRJdW1ckUq7DvHxAvNZIgWSN7W/eahUsccdd1zBjfRJDcufLkcWbQm8xi7Tk4HA3LvkgLmW/mxu+/nc9osqlrb98sHNYq+2X7U23jZsnBx32uayV4SVHCF3IsYVJecSbYKM/STKGqUZSHsR21hPclju6928Pi8gax1uu+22RZrxpc/iMwlyeQ8FgxsbIhBgbw4kAiAh5cT7lfE1RQZZ3jnnucg4ubFLHxs0HD4qdttttzJW7iMuw2QiMEqiTNXrZD4zQ201idFkdgYUNmNCXtgsDY2HUfkxNKk9yQy1gVN4MVgmnr3f2FDH+lAzT9dYLcK0GeNFoJ5X6ljGZayJZ0ny7h3JLNkJgZQoO6GSeQUBBiSmKiBGYz/nn39+kRYZVIRHHW7FkGaGRGCcECCBGzP33HL/Z6iAipaLPNKl59bHHWvgHEcfpzvXvS2jlCiTKLvfl5k84gVjDiBjneuuu66MQZoEz+G3cchDDz20zGkz1SNDIjAJCNCIWJ2EwwVOGahpDRW86lWvKh98ptQYQsjhgvG+m0mU431/Zqp1vrotpGtdSF/jvspJkCZ8I0kvFSrB/Aqfqcdiojvr4w85ctzA4MeSZ+b50oSYe8klns0Upgzji0AS5fjem5lpma9uqzRQUzHYMaZj3p4vb/MhqaxM+EaQ+eU9M4/FVHWUYZSpNYx8TLkxdmns3YcfDQmrbc+7cc4M44fAKIkyjXnG73kYaosQ5A033LDC00m4nPPSYPxgM5aTX9tDvS15sQEg4CPviU98Ytk4JLAx9vFheMkllxR3hNb8ZKBmlRK+djMkAhDIMcoZfA6ooszPM72DkQMjHV/Ygi9rLwoGDyRIUxEyJALTiID/AWO1L33pS80555xTjNa4DESgnMUzWOM6z3ADq9oMo0VglBJlEuVo7/1Qr24eHldnNr4xOeD2gmC8E/5KedWRzpAIzAoCCFOwMsmHPvSh8tHov4Iw3/CGNxQ/xaxjOXzIsfnRPRVJlKPDfqaujCCNQ/JeQuVkQjsvL1aVsAIDR94kyPx6nqnHIjs7hwCyNAzhP8Ey1vqprL1Zw3LFyELWR6T/SIbRIJBEORrcZ+aqrP2MQ/rjcxhgqofxFyvG86pDgjR2k26+ZuaRyI52QYAkyeDHWD1H/xaOtjoJC1mESSVr84GZH5RdQBxQ9iiJMo15BnRTx6FaX8cMFXwhkyBZ+HGsbfzRslc8lKSadRzuVLZhXBDgrYfXHqpWRmy8/NDC+P9wUI80fXTG/4d6du4FPi7Nz3YMCIH2He60L6+9/dZcnq12iC6NeG3pFH0OhFEEX8QMEkzvsICyMcibb765+Ak1xYODaBuJMj3qjOIO5TUnDQEaGGTJfyzCNLfYgty0MQx+QiOTEuZg7+woJcpOxFj3tk2QsZ9EWaM0JmnedFizsuKzsge1EQfl3HWZUM1pQI6xjMnNymZMHALGMFmHH3PMMUVL479mkWyOOEynoo41fDFIgx8fwsZTSb6zJskmUU7cX2Y8G4wYEeTVV19dLFt96R599NFlErVVE3JdvvG8b9mqyUAAQcWC2jz8cInHaYG1Q02rsgbmK17xioEOZ5jOog3cSc7a3OYkysn4n4xlKy3uawySitWcyPvvv7/M/dpll10aq8CbB2bMJQ11xvL2ZaMmFAEGP+YhG6+kkqXFoXo1vOF/x0KWh59+Spf+2+9+97vL3E8aIqrfWQqjJMo05pnQJ40lK+MCVqzUQXxX+qr1B2VosNlmm6VnkQm9t9ns8UfAOpc2H6IMevz3jGNyj2djOOd/aAm6tdZaqy/SH7sDhEzlyyAvw/AQSKIcHtbLvhLHAP4k3/nOdwpBxuK0pMXtt9++2X333YtXHX/eWRu/WDa4WUEisAQELKDNk9XGG29cLMypYkPDQ9pEmAx+GP/wnbycsUUE/JOf/KQQNGvcDMNDIIlyeFgv+Uq/+tWvykRo1qykx/e85z1ljhcy9CX7yle+shgUsGJNglwyzHliIrBkBHjt2WGHHcrC0PzFHn/88c0VV1xRDH+4iNxvv/2aww47rBjXKUtNu9j/6ve+973iFMGHMNLNMDwEkiiHh/WSr2Q+5Lnnnlu8hVDvWFDZnxJBchrw6Ec/ellfqktuWJ6YCCQCKyHAypxDdZaxpmdxifeZz3ymede73tWceOKJzRFHHNHsv//+zXOe85yVzutlx7go71pWO2F3kGF4CCRRDg/rRV2JZVusn8ejjq9Ua0VusskmzZ577rlihXZ/mFmzflsUkFk4ERgiAlSrrMttG2ywQfPmN7+5fMxeeumlxUL2E5/4RFHRMsQxXEJl22uwpibtkv97/ud7Ra0/5ZIo+4NjX2u54447ijUd6fHKK68sE5wZBOy6665FtcNIgIeddBjQV9izskSgrwjQ9BibZP3Kyw/iZHxH0jSNywLpxjcZ/FgcHbkuFMyjzDB8BJIoh495xyuSIH/4wx82t99+e7Ge++xnP1v+SIwFfH2a7oEgcxC/I3yZmQiMLQII84UvfGHZNt100zKdBGHazMekITKlhL2B/zfVaqfAaM+4JsnS+6KfU086XS/z/hOBJMr/xGIkKV+IHnxuskiP73//+8sXJ9UKjzocBhiHJEFmSAQSgclGACFuueWWZZk7DgtIlu9973vLmrCcFRxwwAFlLiZjn7Zx3mqrrVbIkYcg1u85Tjm8ZyGJcnhYd7yScUhjkP4wHJj7amSoY0LxtttuWwgyxyM6QpeZicBEImAc09gkdeuRRx7ZnHDCCc15551XnAnwrmWOJAtZzgvq/z5p0wczYz5TxJIoh3f7kyiHh/WKK1GbmOrhT3HZZZeV8UiWrZyVmwtpLCM86izWhHzFRTKRCCQCY4mA/zQC5LDAVBFao2222aZMJzGlxILqjPcMtXgnUNsqz6/sE57whGLkxwsX69oMw0EgiXI4OJerIEhfguZCWoXA+ISvQ6RoUN8fw5emL8UkyCHemLxUIjACBEiWNEgI0LxIEibjH2OXvG6Zf8mwT5rTAhKlVX9uueWWon06+OCDR9Dq2bxkEuWA7ztHyra77767GOogx8svv7y58cYbi+srX4w77bRTsXyzBl6GRCARmD0EWLyySTC/EimeeeaZhQy/+tWvFumSpawFo41dIljk6aPbguu9WMvOHqL97bFls+rQaV9ee8tltmrUuqS5nDPwLqZOMXhv4VdWbf4Ub3rTm5otttiiOAzoUkVmJwKJwIwiQP162mmnNWeddVbDKw+1K4nSaiU/+MEPmte//vXFeQGNFPKc9jCnZTtgro8/m9t+Prf9ooqlbb98cPuPudj2q9Zmbo3tN122uewVQZkVoRMxrjg4l2gTZOwnUdYodUn74vOQc2fFqwY1C6s3KhNxqF5yPmQXADM7EZhhBHjh8aHNIv70009vOCtAmPJNDXnUox5ViNTHtrHOaQ+jJMqZU72ajmH8b5BjgOZCmgd5ySWXFKfISJLaxHxIA/C+AGfhwZ72P272LxEYJALUrIz+DNv4qCZN2uehy5QyC7Ufe+yxhTTZOAxjXqVhJGGQ789BYrrUumeCKLl9suApff+dc0vV7LzzzkvytdgryAbbP/3pTxeSXG+99YoEyfsGa9Y06e4VxSyXCMwmAsjIuypWH/HOop0ydxJhckLinUbaNPeaKzyWsYMkSqpelraEAB/7pq5wpDArYaqJ0heXhY0tT8PS1FQMVmTPfvazB0qUHCNzWeU6SPmlL33pTIwhzMqfJvuZCAwCAZLi97///fKBfeGFF5Y1Zu+7776ifTKV5ClPeUrD6YD3C6JkMc/+4UlPelKJ+90mdfMWhqQZD3l/fv3rX2/e8IY3FLd8SZT9RnzI9XmI4qEzRnjGGWeULzTNePKTnzzw1vCkY9zA2KOHOkMikAgkAvMhgJRIa95X/+///b9CfKaMGK4xbEO1iiiHofL0/qTaRdJXXXVVc9JJJxXSZkTEoGgWw1RKlHfddVfzvve9r+Gx30KnvooEaothPGiu081f4yw+ZNnnRCARmB8B3rkQErsGH/nmU/6P//E/ykIIPrZrDz3z17T8o4yHDB1ZIgxp3nbbbWU9XDUP4/25/B70v4apJEqGMptttllReYKMxZiHMEIMSMf+IOJZfaAGgWXWmQhMIwLeQ94TH/jAB8o7ygc+mwb5XFuybkWSrOWHGXzkc9D+1re+tWjgTFHh/IDEO6thKonSYDfHw5a38SAahGZBRr0xDJKc1Ycp+50IJAK9I8A457rrrmtOOeWUYkvB4M+0MTYVfD9zMsDilceeYQZDRgiblzDODBDkl770pSTKYd6EYVzLF9gznvGMFZdiDRZEiTiTLFdAk4lEIBEYEQIkRg5IvJMshLDvvvsWRyQ+7ElwSJTTgbXXXnuoKs/21DXvzvD+o62z+A6dSomy/dwHMcZNzsVP2wjlfiKQCAwbASuBMNZhILPJJpsUQvSuYrjjQ58kxyc0Eh3lMnvGKcMpSth5xDt12JiN6nrT7/doVMjmdROBRCARmAcBUhrnI/w9hwbMx7ywzz77NHw/U3ledNFF89SSh4aBwEwQZTx8voJsvooyJAKJQCIwSgS8l6g1qTrbVq3I08R+E/25rjO/0rzwUQTSJKlSoI3zDo136ijaM4przgRj1Dd1Fm/yKB6svGYikAgsHQES5h577FEcl/DS8/d///cNq1hL9Q071EQZwkb9Th12e0ZxvZkgyhrYWbvBdd8znQgkApODgClu++23X/G8Y144ZwQsYc2zzDBcBGaOKMGbZDnchyyvlggkAotHgJEPrzyvec1risrzb//2b5uPfvSjxcUdskyjxMVjutQzZpIolwpWnpcIJAKJwDAReOpTn1pUsOZVmhfOOcHRRx/dHHfcccVJOlVohsEjMBPTQ9owDvrh+uY3v1nMujkTDuk1BsFdmzFRN4Mi5WzOs3Urpx4D7IspB4cw825jEvsxaK9+ZaP9cTziuj/9Kjeofs+Ho/4Y94l+6vd8qzAMot+uaZvvuXB8sfd7vn7X9cEgnwsorBzG4bnQIvedA3Qec9ynm266qbjlZA1r1ZC99tqrGP50e1es3KvcWwoCM0eUXhCDDhZXtR6lQfh4eD3g0l5e/oBeuJ2CMl7U2qlMEFe7rHpYyi2mnLILGQMESWgHx8jdgnLaYJvP49Fi++16+tyt30vBZ6F+szyM+zHsfsPHJvT7uZiv3/XzM1+50rC5n17vd6/l8rmY/z1QPxf+C3GPPKeW1+KowMZzDw86tgyDQ2AmiNJDFsELot6P/H7Glp+xXpuH2os3QhCL/GiDPKHed479yIvz23FdLspGfXXZKCcvCKE+XqejvepZqGzdn7qOdnox5ZwbfWnXE/vRn4XKBRYLlYv61B/9bp8TdUUZcbuMvDospt/Kqq9dZ+w7Hu2MvPpadXqx5Zy70L1Wp6AdC5VdbL/VF31yrlDvL7Y/zq3PLxVWP+1r1NeLYsrEdeXZ76XfneqKOiNW12LqQ5SmhjDk8THlw46bTquLPPe5zy3zLbUtPrjiOv2KXT/qjrYHvv26xrjXM1NEGTd5oQd+uTdt8803b2wZEoFEIBFYLgIkSJLjn/3Zn5VVPJAWv7BHHnlkec8Mel3Imih9PHiPJlEu967m+YlAIpAIJAJ9QYDjdN55WL4ap+Te7rDDDmu22WabIlnmerd9gXnBSmZCovQFJIRKJtRIC6KTBRKBRCARGBEC//7v/958/vOfb97znvc07B6OOuqoYgHLL+zjHve4obWKBEuqFEJNHu/UoTVixBeaKaIMrGftJke/M04EEoHJQYB1q7UgjU2ybD3ooIOK0c6jHvWooXaiJspQuc7aO3QmiNJXUP1FNNSnLC+WCCQCicAiEfjhD3/YfPrTn26uvfbasqrIW97ylrKY8rBJUrPb1thBlovs0kQXn3qidJMNhodDYftIE3mmCnain91sfCIwtQgYl7ziiivKe4rPV/MlRxWMk1IDC0jSOpqDNogcVV+7Xfc/5y50KzHh+TfccENz8803ryDHH/3oR8WC7F//9V8nvGfZ/EQgEZhWBD7+8Y833/72t5stttiiqFxH2U9SLecpAoHDgtKz5m92KiVKS9Kcd955jZXC6fkRZQRfRMcee2xz+eWXF8/8JurusssuzaBNrOP6GY8XAj/+8Y+LestE+fXXX78sojteLczWzBICJDerhBiXfNrTntZsvfXWzeqrrz5UCEiQX/7yl5uzzz67ueeee0qaKlggYPAIZPkvLvXCcbt5ndOsoZtKomQybTIuFeuTn/zkYkrtJta6de6gnvjEJzZPecpT5nVZNtQnNC82NAS++93vlo+oa665pjiZ9uW+xhprJFEO7Q7khToh8JOf/KS5+OKLCxFZk/IFL3jB0N9PHBpwyL7BBhs0fM1utNFG5d3JgMc71OYd6x1q4Wnlp924ZyqJ0srgpMQMiUAbAX/yO++8s7n++uuL1oGrwfvvv7+8EGIcpn1O7icCw0KANHf++ecX4rFwsw/+YYdHPvKRxXes62d4AIGpJMq8uYlANwQYIVxwwQVljIXGIb6GaxP4budmfiIwaASoNn3EPf/5zy8ajsc85jGDvmTW3wMCSZQ9gJRFpgcBKvgjjjiiqOUvvPDCsgrDP/7jP05PB7MnE42Asb6XvvSlZR3K9dZbb6L7Mk2NT6KcpruZfVkQAWMp4fZrtdVWW7ECC+mSVJkhERglAgx4/vzP/7wY8BhCyjAeCCRRjsd9yFaMAAHEGAZew5pXaxyUpa1pSub3ihG2sSgEfsstt5Q87eGmjIHRKMapRnA78pJzCBgf3HDDDROLMUMgiXLMbkg2Z7gIBFEO66qMNawjaNqSsSiOrtddd90yodwcta9//etlztp9993XGJ9ifr/33nsXC+2HP/zhw2pmXicRSAQqBJIoKzAyOVsIIMkwax8WYfJqYk6aeXIf/OAHC+DGpG677bbGlBXLs/EeZe6vPFa59g888MAyr27QdwgOrmcj3WRIBBKBppl6zzx5kxOBcULAvLMddtih2WmnnVY0i7qVevXkk09u/uqv/qo4xDj44IPLHDXz6j760Y8Wde2KEwaY4HGF+7QPfOADM+embICwZtUTjkBKlBN+A7P5S0cgJlCrISTLdm3GEVnF3nvvveUQSYvUZXvYwx7W8TzHSGNrrbXWKsshsbqlbjWXU7BPotx///0LWapT4CiDMQfPKK7N2GhQQd0k3c985jNFguX28Xd/93eb173udYO6ZNabCEwUAkmUE3W7srH9RqAbQcZ1uOriKeWSSy4pWVzdCQjul7/85QpjoJJZ/fD4dPjhh5cxxiq7JP/pn/6puFe0s/baazfbbbddYypAkKR8xjxhYGTFCNcbVEDs+sJ4yIcAv56mKcjPkAgkAk2TRJlPwcwiQJIKAhIjpnb4xS9+UaQtlqpItd6QSgT5NbFwZNDNcTQp0RilsOWWWzbPfe5zm3r5JPWSYo1lsoQ1+Vx9gwr6jqRf9KIXFUfc3/zmN4uqt+7PoK6d9SYCk4BAEuUk3KVs40AQiLmTSK6bZx7qUyvLGzMUavJwXreAfGryi3LIGAkGUSIn/jTr4DjV7N133924Pp+fVKHtgFDD7V5IuuonHQb58ce5UKjbCgdqY+rhuq8L1ZHHE4FpRiCJcprvbvZt2QiQtJBQOClYTIWdHBhY2cbqELHMm9VrOKCuQ6x849rmUb7qVa/qKFGaQmJc8fOf/3ypD5lusskmjXxjjs973vOKtWxdd6YTgURg8QgkUS4eszxjShBAZKQyklNIl+2uhaqV1NWPcPvttxdpkSp1nXXWKWOBMT+SmjeWMVLOHErSLAIM6U7MifvVV1/dfOpTnyrHqG+pZ0mPJ554YvPFL36xSIXPeMYzFt1kHwUwIZnOJzEvuuI8IRGYYASSKCf45mXTl4cAQuAAQCCB9YsM52sVlar5kdSbpEnkhpCoS7/1rW81Fuy98cYbi3cWcydNJRGCtKhjSZEf+tCHCsm++MUvLsZAJE+BKtZ0E+OeDIUWG2DgowEhxzUXW0eWTwSmDYFREWWa003bkzQh/UEApluQyqzUbpI/YkBSVnJHDixWuZXrd3Btq9ZTvQr277jjjjIeSF2KID/96U8X6XCfffYpjrGtmRqBQdEXvvCF4qhA29/5znc2u+22W3FEQAJUn36RlBGndQQzJAJThMDIeGPQRKljI+vcFD0g2ZU+IYBMvvrVrxapjPoSsZgW8bWvfa057bTTGmOIpmtYMLffgYs6Uy8Y65AoETMn2NSwLE1de+edd24OPfTQ5lnPetYqC/Yi1XPPPbe01cr3O+644wpvPepGwuoRnv3sZ5e5mEvtA5wyJAITgMBQOGbQRBk4153Jf2CgkvHQEaBa3H777RtEQ91JmqSClc94xhZjhv1uHB+v4bjgOc95TnPmmWc25lS6JnLkYACBMhzqpAb+h3/4h+aaa64pxM4Kt5Z6Ef7b3/72QrbqTmmy33cv6xsTBII/ak4ZeNOGRZTRkehk7GecCAwdARLcIOcldusQ/61UrIxzjCGS+qhWqXu1B0HPNy5IEkWInAGYUhJ9MO6JdJEol3ecrDMUqkOoZUmtJGief+L8ulyk52tHlMk4ERgRAkPnkUER5UJsP/SOjuiG5mUTgRUImDvJ04/pIEFkSK/X8OhHP7oY/7CMZdDDcMfYJClVvRwZID9E2bZ4NQZ76aWXFmtZRM1r0Prrr7/SXE91MW5CxlS53Pe5Zifpttc2Z7lEoI8IzMcbC3HOsprRb6IcaGOX1dM8OREYAQJUuzz0mMDPAIe3HQ7QOSNARKTIcBawUPM23XTTZtttty1jrFdccUWx1HU+lS1CM73k937v94q1K6mxDuHkAFkiVmseMvjRDtKmdjJsIvEibwZBlgIzVrr66qsPTB1dtzHTiUAfEBgIB/WDKKNh4k4h8iPuVCbzEoGpRMA4KGK67LLLioWrTpL6LNhsGgfpr1ei3GqrrQppMehhOYscjbca17zooosKfuZePv3pT18FS4THSElbLr/88kLS4dUHSZJSr7rqqtKWPfbYo5yPjJG5aSxUtRkSgTFCIPgk4nbT5MfWPrbo/X4Q5XwXjU5EPF/ZPJYITB0CSIgaEzm+4Q1vKAZEjHfss7QlAfYanIdY15pzaycYR2T4w3rX+KRA6mRJ2w6ILta6tIwWq94gP9KjRaL32muvIlnGuSRNEqstQyIwhggEr0Q8sCYulygX20DlYxtYp7LiRGBcEEAyxgtJf20DGUREGuw1BDHW7vRY7JISqXURHiINAqzrpZY1Tkmd+jd/8zdltZK4tnpt2pghERhjBII7lsI7y+rWUokyGtzLxaNTcU7EvZybZRKBiUaAWhUBDYKErBvJ1R2JkorXWKJ9Ma884ghIlIUtqdI4Z7cpKFE+40RgDBEI7ohYE6V7CfU5vZRfqcxSiXKlSjrsRKOiE7GvaJ3ucGpmJQKJQC8IsE5lgIOE995773IKwjT+ueaaa65UBVUvS9uwtl3pYO4kApOBQM0dndLyBhIe0qq10768Thsv0b/dYUO+lmmvN4vp2Qx2lG1O7XTmXDpDIpAIJAKJQCKwIAJzwwO+Bn9ebb+YS9t+2dosFPurDpsFZ4NgO8Vzh1eElUh3uUsi1BdbcYW5xGLz63MznQgkAolAIpAItBFYLK90K9+ud8H95RLlgheYK1A3Vhqrr7qUfC81ZZlEIBFIBBKBWUUguKPNKQPHo59E2W78QvsD71xeIBFIBBKBRGBqEFiIU9rH+9bx5RClRtWhbqT82K/TkZcSZY1cphOBRCARSAQWQgBvBIdE7JxOaXl1aO/XxxZML4coo/K6ke08x6JzITZHHGUzTgQSgUQgEUgEFkIguCPiml/m46GF6l3weC9EqQEROjWmfSzK1HF0rI7jvIwTgUQgEUgEEoGFEKj5I9I1z7TTneqLMnHM/oJhIaKMSrpVXh93sShXx9Ghdrxg47JAIpAIJAKJQCLwIAJtDon9mm8i7RRpoR1HXqf8ckL7ZyGibJfvtN++mP3YdEQ6OhSxOS4ZEoFEIBFIBBKBXhHAG8EhEQe/BOeIhXb8QO4Sf5dDlHXD2ulunZAfnV1ik/O0RCARSAQSgRlEoOaP4Bx50jXnxLF2vGTIFuvCzoUj1F582g2K/XYnojMpUQaKGScCiUAikAj0gkAIWcEjNb8E57Tjut76WJ2/YHoxROki85Gji9UNiXTdKR0lxS5HknWdDIlAIpAIJAKzhQD+iK3mleCaOoZMvS9dh/Z+fWyVdK9EqdKaJFepaC6j3ai6I9I6qI4kyjkQMiQCiUAikAgsCoEgyYjbHNPmoIUq75kseyHK+Uiy3bD2frsjNWEu1Ik8nggkAolAIpAIBALh7LwTr7S5p70fdbRj5RYMvRBlp0rqytsNsq8jIT1Gp0KilL+QdNrpmpmXCCQCiUAiMLsIhCQZcXBLcI79TnxU5y0JvTZRqnAhEmuXqRshHSRZp4MkI17oGkvqTJ6UCCQCiUAiMLUIBEGK8UzENUHW6eCmGhB5C4VVynQjrHa+/chrp2PMUX6kxZ3WqpSHnOu4nY7z6rqko/46nssu+fLqsNB+XTbTiUAikAgkAktDoE0qnfYjTxxEVsfSsdVkKB3q1jpd57XLRz3t+us2dEpH7+NY7Je4LVGudLC1o4IgoHY6KhfHpqGdQtThWJQVRwd1HFl2IkrnBmlKR10Rz2WV0N6P/IwTgUQgEUgEBoeAd3kdYr/9ro/9eO9H7P0v3SbAhfadE3W242iP/Ah1OvK6xoshyqjEBYKIIh0XFbc3HVA+zplLdiwTQCHCOKcmy6hDLF+IvAf2Hvitr1PnZzoRSAQSgURg8AgEH8SVghPstwkt3vuR3yZK+Z1Isn1eXKMdu2bdnjrtWE+hV6JUOQKKWOWRri8sbdOJXkKUR3xxnjRgxEGKkQ5ijHiuyEoELL8O7f36WKYTgUQgEUgE+oOA93cd6v1Ii9tbEJ78SEdck6a8IMx2OsqL2/XbFyKu03VeKdTtp1eijMqDeFygTjtuX77G1sQ3t7tSiMaJqVjFysd50kGMdVr9sc0lV6SjHfLq0C2/LpPpRCARSAQSgf4gEO/2dm3yY3Ms0hEH0dmPdDvuRJp1mbouaSHyHth7YL9TOvK6xoshyqjExYOE2um6jE4EAUa+uN2JIMOI1S3dKXZ+XFsc6TpfWqiPPZCTv4lAIpAIJAKDQiDe7VF/vS8d+3WMJ+x3imsijHQQZuy3z4vrxDW0pVvasZ7CUogyKnZxZFQ3oj4mrRNtsozyYluU6UaOQYgRz52yggTrPPkR5GdIBBKBRCARGC4C8X6vrxrvenlxPPLquE169us8aUTpnPaxup52eq74iutKLzosRCidjkdep1hebEF89iNdx3W6LtNOR30R66S0EPEDe6vuR37GiUAikAgkAsNDAFnVIfbruE1o9oMY63TkdSLHOFaXr+vVBvtCO34g94HfOFbnrUgvRaJUIYJqxyqtL6YDyFCQDlKLc9vnRxnHg0SlY5tLrlRHvS8tKJshEUgEEoFEYLQI1FygJbHfKZYXW0187bz6WJ1Wrr3vmoJjQjt+ILfH34WIZb7jcUw8XzqOdyK/TnlRvo51p75GvS8txPEH9v7zt1v+f5bIVCKQCCQCicBSEQgSap/fzo/9TrG8+bY2ESrbKU8bop5OaXmdgnO6hqVIlFGZipFQXCAIKfajXMQ6pUy9RR11XtRT56mjzrcfIfJjv47nO1aXy3QikAgkAonA0hHo9t5XY/tY7NexdKf9yO81rnsQ53RqQ11uwXSvRDJfuTgmrtMuHnl1fp0X6YXKOi7U9TyQs/BvnLNwySyRCCQCiUAisFgEguB6OS/KRuycILTIm29/vmNRV8R1ffI6hSjT6diKvF5JZKFycVwcaRep9yNfHGOXdV63slEm6hMLdf4DOd1/F1O2ey15JBFIBBKBRKBGoCeiefCEumw7Hft1LF3vqybUrdL1sUhHfuxHLL9TWOh4OWexBDJf+TjWS1yX6ZTWuDo/Ohh59fE4Nl9cnzdfuTyWCCQCiUAi0DsCPRHNg9XVZTul23mxL+6WVnV9rNO+vHaIc9r5HfcXSyALlY/j3WKN6HSsnVeXa6c77curQ9RX52U6EUgEEoFEYLAILERA7eP1fqe0vMhvx3rSzmvvd+ttlOt2fKX8xRJKL+WjTDt24cir0wvlRYPrcvX5cTzjRCARSAQSgfFFoE1O7X0tj7yIF8rrdLw+1/FOoZcyK85rk8+KA10SvZaPchGrrpf0fOXaTarrax/L/UQgEUgEEoHxQmA+cqqP1Wk9qPcXStfH5+t9r+VKHUslm17Oq8v0ko5O1WXltfejXMaJQCKQCCQCk49Am7Ta+3pY5/WS7oZKfW63MqvkW71jKaFX8pqv3HzHltKmPCcRSAQSgURg8hBok1d7P3oU+RFHvrhTXn18WemlEmX7or2S3kLl2sfb++3r5n4ikAgkAonA5CPQJrr2vh7WeXV6vt73Wm6+Osp6kPMW6PHgYgmtU/l2Xnu/x6ZksUQgEUgEEoEJRKBNau19XeqUN/Cu9kuijIYmuQUSGScCiUAikAgsBoE2Cbb3l1PXYs5dpey4EmUS7iq3KjMSgUQgEZh6BJZDjgMDZ9CENOj6BwZMVpwIJAKJQCIwMQgMlGDD5+rEoJENTQQSgUQgEUgEhonAqCS+UV13mNjmtRKBRCARSAT6i8BAJcduTU2JshsymZ8IJAKJQCKQCCQCiUAikAgkAolAIpAIJAKJQCKQCCQCiUAikAgkAolAIpAIJAKDQOD/AwNdJrecwvKaAAAAAElFTkSuQmCC"
    }
   },
   "cell_type": "markdown",
   "id": "a85f4849-44de-4736-9db6-bc95b2675b43",
   "metadata": {},
   "source": [
    "P(j|0) =  [0.95 0.05]: |0> keep to be |0> with probability 0.95, and flip to |1> with probability 0.05\n",
    "\n",
    "![readout_error.png](attachment:f74cf98a-ca61-4db6-aa26-859270981222.png)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6c8c878d-d856-472e-a1b9-3a1b7d451115",
   "metadata": {},
   "source": [
    "### Other noise types see https://quantum.cloud.ibm.com/docs/en/guides/build-noise-models"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8ee98edd-a1d4-4736-886a-f46f10e3b1b2",
   "metadata": {},
   "source": [
    "# Add noise model to quantum gates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e71e3f1b-81a7-45e7-92e9-41e632c1b561",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create an empty noise model\n",
    "noise_model = NoiseModel()\n",
    " \n",
    "# Add depolarizing error to all single qubit u1, u2, u3 gates\n",
    "error = depolarizing_error(0.05, 1)\n",
    "noise_model.add_all_qubit_quantum_error(error, [\"u1\", \"u2\", \"u3\"])  #add depolarizing noise to basic gates of \"all qubits\"\n",
    " \n",
    "# Print noise model info\n",
    "print(noise_model)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c65be445-7726-44f6-82d4-4df50a0fcaed",
   "metadata": {},
   "source": [
    "## add bit-flip (X) nosie in quantum circuits"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0d691f1b-ea7f-490a-86e0-b5af452563a1",
   "metadata": {},
   "source": [
    "### noiseless circuits"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "46e5bf02-921c-4f28-b8e8-d861ee6f853d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# System Specification\n",
    "n_qubits = 4\n",
    "circ = QuantumCircuit(n_qubits)\n",
    " \n",
    "# Test Circuit\n",
    "circ.h(0)\n",
    "for qubit in range(n_qubits - 1):\n",
    "    circ.cx(qubit, qubit + 1)\n",
    "circ.measure_all()\n",
    "print(circ)\n",
    "\n",
    "# Ideal simulator and execution\n",
    "sim_ideal = AerSimulator()\n",
    "result_ideal = sim_ideal.run(circ).result()\n",
    "\n",
    "plot_histogram(result_ideal.get_counts(0))\n",
    "plot_histogram(result_ideal.get_counts(0))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d96a66a2-f334-47cd-826c-702ab26da829",
   "metadata": {},
   "source": [
    "### Add bit-flip noise in the initialization (\"reset\"), readout (\"measure\"), single-qubit gate(\"u1\",\"u2\",\"u3\")，and two-qubit gate(\"cx\", i.e., \"CNOT\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "25d78001-97be-482c-ae33-5bf3584cf417",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Example error probabilities\n",
    "p_reset = 0.03\n",
    "p_meas = 0.1\n",
    "p_gate1 = 0.05\n",
    " \n",
    "# QuantumError objects\n",
    "error_reset = pauli_error([(\"X\", p_reset), (\"I\", 1 - p_reset)])\n",
    "error_meas = pauli_error([(\"X\", p_meas), (\"I\", 1 - p_meas)])\n",
    "error_gate1 = pauli_error([(\"X\", p_gate1), (\"I\", 1 - p_gate1)])\n",
    "error_gate2 = error_gate1.tensor(error_gate1)\n",
    " \n",
    "# Add errors to noise model\n",
    "noise_bit_flip = NoiseModel()\n",
    "noise_bit_flip.add_all_qubit_quantum_error(error_reset, \"reset\")\n",
    "noise_bit_flip.add_all_qubit_quantum_error(error_meas, \"measure\")\n",
    "noise_bit_flip.add_all_qubit_quantum_error(error_gate1, [\"u1\", \"u2\", \"u3\"])\n",
    "noise_bit_flip.add_all_qubit_quantum_error(error_gate2, [\"cx\"])\n",
    " \n",
    "print(noise_bit_flip)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f51fa53b-abd8-4f2d-b582-defcb550f319",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create noisy simulator backend\n",
    "sim_noise = AerSimulator(noise_model=noise_bit_flip)\n",
    " \n",
    "# Transpile circuit for noisy basis gates\n",
    "passmanager = generate_preset_pass_manager(\n",
    "    optimization_level=3, backend=sim_noise\n",
    ")\n",
    "circ_tnoise = passmanager.run(circ)\n",
    " \n",
    "# Run and get counts\n",
    "result_bit_flip = sim_noise.run(circ_tnoise).result()\n",
    "counts_bit_flip = result_bit_flip.get_counts(0)\n",
    " \n",
    "# Plot noisy output\n",
    "plot_histogram(counts_bit_flip)\n",
    "plot_histogram(counts_bit_flip)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "64524014-8801-48c7-9515-0d950d31fa41",
   "metadata": {},
   "source": [
    "# Simulate noise model from the true hardware"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "05c806f7-7b45-4c6e-9bfa-15ad61bf3c00",
   "metadata": {},
   "source": [
    "### Optional: refresh the noise model from IBM Quantum hardware"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "87d1cebb-07dd-4b88-b7b0-54ea8a246511",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(\"Optional hardware download skipped; using the bundled IBM Kobe snapshot below.\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f89c3849-1a91-4afe-84c4-a4ef41c26354",
   "metadata": {},
   "source": [
    "# Check the noise info"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bfe21242-42b0-4a37-b8b0-d177b612d7d7",
   "metadata": {},
   "outputs": [],
   "source": [
    "import gzip\n",
    "import json\n",
    "from pathlib import Path\n",
    "\n",
    "def decode_numpy(obj):\n",
    "    if \"__complex__\" in obj:\n",
    "        return complex(*obj[\"__complex__\"])\n",
    "    if \"__ndarray__\" in obj:\n",
    "        return np.array(obj[\"__ndarray__\"], dtype=obj[\"dtype\"])\n",
    "    return obj\n",
    "\n",
    "snapshot_path = Path(\"noise_model_ibm_kobe.json.gz\")\n",
    "if not snapshot_path.exists():\n",
    "    snapshot_path = Path(\"notebooks/lecturers/xiayang\") / snapshot_path\n",
    "\n",
    "with gzip.open(snapshot_path, \"rt\", encoding=\"utf-8\") as f:\n",
    "    snapshot = json.load(f, object_hook=decode_numpy)\n",
    "\n",
    "noise_model = NoiseModel.from_dict(snapshot[\"noise_model\"])\n",
    "coupling_map = snapshot[\"coupling_map\"]\n",
    "basis_gates = snapshot[\"basis_gates\"]\n",
    "\n",
    "noise_model_dict = noise_model.to_dict()\n",
    "len(noise_model_dict['errors']) # record the probability of adding noise model to basic gates (\"sx\",\"x\",\"cz\") of all 156 qubits. \n",
    "                                # \"Rz\" rotation has error - \"virtual Z\" gate\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4c22df83-f379-4dfe-9e99-2303f7978ac8",
   "metadata": {},
   "outputs": [],
   "source": [
    "noise_model_dict['errors'][800]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9f850d93-ce85-43cc-9a87-ab065a869207",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "77409c73-d225-4055-842c-e0d162d7f356",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1c57f16a-6bc1-4595-be4b-708442629ed5",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "db6309ae-be2d-4ec4-a82e-0db1c42fa45e",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.12.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}