diff --git a/docs/tutorials/readout_mitigation.ipynb b/docs/tutorials/readout_mitigation.ipynb new file mode 100644 index 0000000000..d641bda3cc --- /dev/null +++ b/docs/tutorials/readout_mitigation.ipynb @@ -0,0 +1,419 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Readout Mitigation\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Readout mitigation is part of `qiskit-terra`. The readout mitigators can be initalized based on existing backend data, or via a readout mitigation experiment.\n", + "\n", + "In a readout mitigation experiment, simple circuits are generated for various combinations of \"0\" and \"1\" readout values. The results give us a matrix describing the probability to obtain a wrong measurement. This matrix is used to initialize the readout mitigation object, which is given as the result of the expriment." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "from qiskit import QuantumCircuit\n", + "from qiskit.visualization import plot_histogram\n", + "from qiskit_experiments.library import ReadoutMitigationExperiment\n", + "# For simulation\n", + "from qiskit.providers.aer import AerSimulator\n", + "from qiskit.test.mock import FakeParis\n", + "\n", + "from qiskit.result.mitigation.utils import (\n", + " expval_with_stddev,\n", + " str2diag,\n", + " counts_probability_vector\n", + ")\n", + "\n", + "backend = AerSimulator.from_backend(FakeParis())" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [], + "source": [ + "SHOTS = 1024\n", + "qubits = [0,1,2,3]\n", + "num_qubits = len(qubits)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Standard mitigation experiment\n", + "\n", + "The default mitigation experiment is *local*, meaning error probability is measured individually for each qubit. The experiment generates two circuits, one for all \"0\" and one for all \"1\" results." + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + " ░ ┌─┐ \n", + " q_0: ─░─┤M├─────────\n", + " ░ └╥┘┌─┐ \n", + " q_1: ─░──╫─┤M├──────\n", + " ░ ║ └╥┘┌─┐ \n", + " q_2: ─░──╫──╫─┤M├───\n", + " ░ ║ ║ └╥┘┌─┐\n", + " q_3: ─░──╫──╫──╫─┤M├\n", + " ░ ║ ║ ║ └╥┘\n", + "meas: 4/════╩══╩══╩══╩═\n", + " 0 1 2 3 \n", + " ┌───┐ ░ ┌─┐ \n", + " q_0: ┤ X ├─░─┤M├─────────\n", + " ├───┤ ░ └╥┘┌─┐ \n", + " q_1: ┤ X ├─░──╫─┤M├──────\n", + " ├───┤ ░ ║ └╥┘┌─┐ \n", + " q_2: ┤ X ├─░──╫──╫─┤M├───\n", + " ├───┤ ░ ║ ║ └╥┘┌─┐\n", + " q_3: ┤ X ├─░──╫──╫──╫─┤M├\n", + " └───┘ ░ ║ ║ ║ └╥┘\n", + "meas: 4/═════════╩══╩══╩══╩═\n", + " 0 1 2 3 \n" + ] + } + ], + "source": [ + "exp = ReadoutMitigationExperiment(qubits)\n", + "for c in exp.circuits():\n", + " print(c)" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": {}, + "outputs": [], + "source": [ + "result = exp.run(backend).block_for_results()\n", + "mitigator = result.analysis_results(0).value" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The resulting measurement matrix can be illustrated by comparing it to the identity." + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAASoAAAEDCAYAAACCmUnRAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAYLklEQVR4nO3dfbBdVXnH8e/vnhBAICbhRkVAwUKrEZ1EXqe1WHkzzFhgLChIIViUsZbaGUcFtKIDOC3TqY5MEYkVEN9fqJoZoamIdDpVMAEiEKwloIVEVPICIWACSZ7+cfYJh5tzz9573XNv9ln5fWb23HP33mutfTJzn+y19l7rUURgZtZkIzv7AszMyjhQmVnjOVCZWeM5UJlZ4zlQmVnjOVCZWeMNJFBJuqDp5XJtK7Vcrm2llhuGtnYGSQsk/ULSSkkX9zj+AUkPSLpX0g8lvbLr2EJJDxbbwq79h0u6r6jzKkkqvZCImPAGLGt6uVzbGoZr9L/HYNqa6g1oAQ8BrwKmAz8D5o45583Ai4rPfw18o/g8G3i4+Dmr+DyrOPZT4BhAwC3AyWXX4q6fmY3nKGBlRDwcEc8CXwdO7T4hIn4UEc8Uv94BHFB8fgvwg4hYFxHrgR8ACyTtB8yIiDuiHbVuBE4ruxAVEa6S0Zkz4qCXzdlh/+NPbGDOzBnjFxxp9dz9+PonmTPrxb3LbNs6bnV923tR7/2Pr1nDnNHR3mW2bhm/rbXrmLPv7N4HW9PGaWstc0b3HbfOcdvqd419y/Vpb8tzvcusW8+c2bPGr3TabuO0Vf8aJ+V7DbhcU9r61SOPsGbN2vKuUB8Halpsotrf9Rq2rQA2de1aFBGLACSdDiyIiHcXv58DHB0RF/aqS9K/AL+JiCskfRDYIyKuKI59DPg9cDvwjxFxQrH/T4GLIuKt/a6z91/aOA562Rzu/PyVdYq07dUniI3n6Q31ywAj84+vX2jjuqS2NKP+Hx9Anf8ctrdVoRvfs611jyWVY9bL0solSP1uOTrijX824To2E5zBXpXOvYanNkXEERNtU9JfAkcAb5poXb2462eWoRGp0lZiNXBg1+8HFPteQNIJwEeBUyJic0nZ1TzfPRy3zh2+T9kJZjZcRPsPu8pWYilwqKSDJU0HzgQWv6AtaT5wLe0g9buuQ0uAkyTNkjQLOAlYEhGPARskHVM87TsX+F7ZhdTq+pnZcBip2pvuMwoREVskXUg76LSA6yJihaTLaD+5XAz8E7A38K2iC/9IRJwSEeskXU472AFcFhGdMZb3ATcAe9J+6ndL2WU6UJllRohpAxr3i4ibgZvH7Lu06/MJfcpeB1zXY/8y4LA61+FAZZah3MZ0SgNV8RbtBQCveGnaUy4zmzqiRtdvSJQG3ohYFBFHRMQRfd+VMrPGGNBgemO462eWG+X3bpoDlVlmOq8n5MSByixD0/K6oXKgMstNezA9r0jlQGWWoV2767fHXoy8+sjajWz77udrl2md/aHaZQC23n1r/bbeMO47a33FhjVJ5VImM6dMZAbQ7P2SysX639RvK3Ei81RO0t4V5Ph6gu+ozDK0a99RmVnjCQY2haYpHKjMMuQ7KjNrNMljVGY2BEbIK1I5UJllKLc7qtKurKQLJC2TtOzxtWlri5vZ1BngCp+NUW/1hPGysZhZY3QWzquyDQt3/cwytMt1/cxs+KjiVlpPeUr3YyXdLWlLkQews//NkpZ3bZsknVYcu0HSL7uOzSu7Dt9RmWVmUFNoJLWAq4ETgVXAUkmLI+KBrtMeAc4DPthdNiJ+BMwr6pkNrAT+o+uUD0XEt6teiwOVWYYG9HrC9pTuAJI6Kd23B6qI+FVxbFufek4HbulK/V5b/UA1Tnr2flImGG+5+iO1ywCMnPi2pHIpkjMlJ0xmTm0rVcoE45SJzADMfGlaOetpgC987g882vX7KuDohHrOBD41Zt8nJV0K/BC4uCtxaU8eozLLUKviBox2Xj8qtgsGeR2S9gNeRzs3YMclwKuBI4HZwEVl9bjrZ5aZmgvnrYmII8Y5Vimle4m3A9+JiOc6O4psyQCbJV3PmPGtXnxHZZahAT31K03pXsFZwNdecG3tuyyKlO6nAfeXVeJAZZahQQSqiNgCdFK6/xz4Zielu6RTACQdKWkVcAZwraQV269BOoj2Hdl/jqn6K5LuA+4DRoEryr6Pu35mGRrU+54VUrovpd0l7FX2V7QH5MfuP67udThQmWUmx+QOnpRsliFPSjazxpOqbcPCXT+zDMkL55lZk1WdcDxMHKjMMuRAZWYNJ1rDNABVgQOVWWbc9RtpoT33qd1IbNpYu0zrvA/XLgPw3CXvqV1m+lXfSmorVVJK9w1rE9vaN6lcUlupKd2f+G39Ql5xYXxD9kSvCt9RmWUoszjlQGWWI+f1M7NG8xiVmQ2F3LLQOFCZZchvpptZow0qC02T1Fs9YU3aI3Izm1qDyuvXFPVWTxidundyzCxdboHKXT+zDOW2cJ4DlVlmxHAtildFbt/HzBhc10/SAkm/kLRS0sU9jh8r6W5JWySdPubYVknLi21x1/6DJd1Z1PmNIsNNXw5UZhmSVGkrqaMFXA2cDMwFzpI0d8xpjwDnAV/tUcXvI2JesZ3Stf9K4NMRcQiwHji/7PvU6/pt3UI8VX/ddO0zdUsYp0ww3vrTm8tP6iV1wu/TG2oXGZl/fFpbQ0AJE4xj3WPlJ/WSOHE6RVkgmNS2B1PNUcDKiHgYQNLXgVOBBzonFJlmkLSt0nW1/1GOA95Z7Poi8Angmn7lfEdllpmq3b4imPVL6b4/8GjX76vokf6qjz2KOu+QdFqxb1/giSJnYOU6PZhulhuJVvU3PvuldJ+oV0bEakmvAm4rko4+mVKR76jMMqQRVdpKrKad6bjjgGJfJRGxuvj5MHA7MB9YC8yU1LlJqlSnA5VZZsTA0mUtBQ4tntJNB84EFpeUaV+DNEvS7sXnUeBPgAciIoAfAZ0nhAuB75XV50BllpuKQaosUBXjSBcCS4CfA9+MiBWSLpN0CoCkIyWtAs4ArpW0oij+GmCZpJ/RDkz/GBGdQfiLgA9IWkl7zOoLZV/JY1RmGRrUE8eIuBm4ecy+S7s+L6XdfRtb7sfA68ap82HaTxQrc6Ayy1BmM2jKA1XxuPICgFcc8PJJvyAzmxgBI5mt81Jv9YTZU/fippklUntScpVtWLjrZ5ahIYpBlThQmWWnfB7fsHGgMsuMAGX24pEDlVlulN9ger1A1ZqWtBJCbFxfu4z2nlW7TKqR1/5xUrltN30uqVzr3B2W9Sm19e5b09p6wwlJ5ZpOs/dLKhfrf5PWXsKqC+2XsHcOd/3MrPEyi1MOVGa5aafLyitSOVCZ5abahOOh4kBlliGPUZlZo7Wn0OzsqxgsByqz3KjSonhDxSndzTI0oIXzGsMp3c0y5EnJZtZonaWIc+JAZZah3J76ZfZswMw6c/2qbKVVJaZ0lzRP0k8krZB0r6R3dB27QdIvu9K9zyu7Dt9RmWVoEDdUXSndT6SdKHSppMVdSRrg+ZTuHxxT/Bng3Ih4UNLLgbskLYmIJ4rjH4qIb1e9Fgcqs8y0x6gG0vVLTukeEf/b9fnXkn4HzAGeSLmQKQlUKSshpKy4kNoWI62ktlJWQQDYctVFtcuMLDgjqS17oZRVECBx1YWZL01qa8JUaz2qUUnLun5fFBGLis+9UrofXftypKOA6cBDXbs/KelS4IfAxRGxuV8dvqMyy06tFT4nM6U7kvYDvgQsjIjOXdclwG9oB69FtPP8XdavHg+mm+WoNVJt629CKd0lzQC+D3w0Iu7o7I+Ix6JtM3A9FXL8OVCZ5UbtMaoqW4mJpHSfDnwHuHHsoHlxl4XaF3AacH9ZfQ5UZjkaUbWtjwmmdH87cCxwXo/XEL4i6T7gPmAUuKLs63iMyiw7g5vIN4GU7l8GvjxOncfVvQ4HKrPMSHj1BDMbApktn+DVE8wypNZIpW1YuOtnlhuVD5QPGwcqswzltnqCA5VZjnxHZWaNluHKeY0NVKkp3acyfXxs2phUrnX+JbXLPPvBv0pqa/dr/i2pnL1QUkr3J35bv6Gtz9Uv04NaDlRm1mQZZqFxoDLLkbt+ZtZ4vqMysyaT/HqCmQ0D31GZWbMJjQzP9JgqSgOVpAuACwBeceCBJWeb2U4nsruj8qRkswwNaIXPxnDXzyxHmd1ROVCZ5WbI1pqqIq8RNzMD2it8VtlK60lM6V4cWyjpwWJb2LX/cEn3FXVepQp9UAcqs9yIgaTL6krpfjIwFzhL0twxp3VSun91TNnZwMdpJyw9Cvi4pM6k2muA9wCHFtuCsq/kQGWWoQENpm9P6R4RzwKdlO7bRcSvIuJeYNuYsm8BfhAR6yJiPfADYEGRKmtGRNwREQHcSDtlVl/1xqgiiOf6Zl7uSbvtXrtMqtSVEJLa2mPvKWsrdRWErTdfl1ROf3BY7TKx8cmktkbmH1+/0MZ1SW1pxmhSufbfVM22UlK6t3arX2bHlusMpk9WSvdeZfcvtlU99vflwXSzHDUkpfuguOtnlpvOwnkTz0IzkZTu45VdzQvzAFaq04HKLDuCVqva1l9ySnfa2ZVPkjSrGEQ/CVgSEY8BGyQdUzztOxf4XlllDlRmORrAHdVEUrpHxDrgctrBbilwWbEP4H3AvwIrgYeAW8q+jseozHIzwDXTU1O6F8euA3Z4mhMRy4BaT2scqMxylNmb6TVXT+gZOM2sUQSZLfNSb/WEfb16gtlQGMxTv8Zw188sNyK7OyoHKrPs5Nf1c6Ayy9EQdeuqcKAyy41TupvZUNilA5WUtBJC01dcyNnIm96WVG7rtZfXLjPtA/+c1tby22qXac07Lqmt2LA2qZxmDM8Tb+2KWWjMbMj4qZ+ZDYVduutnZkPAryeY2TDwHZWZNdqu+HqCU7qbDRtVWRRvqDilu1mOPCnZzBptV+z6mdmwye+pX17fxszaBtT1q5DSfXdJ3yiO3ynpoGL/2ZKWd23bJM0rjt1e1Nk59pKy6/AdlVmOBtD160rpfiLtRKFLJS2OiAe6TjsfWB8Rh0g6E7gSeEdEfAX4SlHP64DvRsTyrnJnF2unV+I7KrPcaGDpskpTuhe/f7H4/G3geO2YK/6somyyKbmjSprIvOXZtLamTU8ql62RtMfUKROMt1x6flJbI+e8L6lcitTJxSmTmXfqRObqd1QTTem+/ZyI2CLpSWBfYE3XOe9gxwB3vaStwE3AFRER/S7SXT+zHDUkpbuko4FnIuL+rt1nR8RqSfvQDlTnADf2q8ddP7PcCNBIta2/Kindt58jaRrwYqD79vNM4GvdBSJidfHzKeCrtLuYfTlQmWVHMFJx669KSvfFwMLi8+nAbZ1unKQR4O10jU9JmiZptPi8G/BW4H5KuOtnlqPyu6VSxZhTJ6V7C7iuk9IdWBYRi4EvAF+StBJYRzuYdRwLPBoRD3ft2x1YUgSpFnAr8Pmya3GgMstN56nfAFRI6b4JOGOcsrcDx4zZ9zRweN3rcKAyy9GuNoXGqyeYDaEBdP2axKsnmOXIqyeYWaMpv0nJDlRmOUqckdBUDlRmuVGld6SGigOVWY4yG0x3oDLL0RANlFfR2ECVugpCyqoLOa+4oD33SSoXmzbWLtO6KC2l+8ZzT69dZp+bbk1qK1XKSghJ6eO3bqlfZgfyHZWZNZzwGJWZDQE/9TOzRvNTPzMbCh6jMrPG81M/M2u2XfCpn1dPMBsyYmDrUTWFV08wy5FXTzCzZstv9YS8vo2ZFVlodnpK94Mk/b4rbfvnusocLum+osxVPRKW7sCByixHA0iX1ZXS/WRgLnCWpLljTtue0h34NO2U7h0PRcS8Yntv1/5rgPcAhxbbgrKv40Bllp2Kd1PlNzKDSun+/JVJ+wEzIuKOIq3WjcBpZReS3RhVygTjeG5zWlsJqeqnWjy1Lqmc9pk94CsZX8oE462LF5Wf1IP+aF5SuXh6Q+0yI/OOq99QawB/kvWe+k1WSneAgyXdA2wA/j4i/qs4f9WYOvcvu8jsApWZ1XqParJSuj8GvCIi1ko6HPiupNemVuZAZZajwbx6UCel+6rulO5Ft24zQETcJekh4A+L8w8oqXMHHqMyy9EABtOZQEp3SXOKwXgkvYr2oPnDEfEYsEHSMcVY1rnA98ouxHdUZrkZ0OoJE0zpfixwmaTngG3AeyOiM2D6PuAGYE/glmLry4HKLEcDmuuXmtI9Im4CbhqnzmXAYXWuw4HKLDva9RbO86Rks+FT4WXvoeJJyWa5EYMaTG8Md/3MsrMLrkdlZkPIa6abWaOJXW8w3cyGjbt+ZjYMMnvq50BF+ioIw7DqQuoqCLFxff229p6V1FaKkePfnlRu69UfTyo37cOfqd/W8tvqN/TMU/XL9OI7KjNrNCcgNbOh4DsqM2u2XXAKjZkNIQ+mm1mjdabQZMSByiw7+eX18+oJZhny6glm1nxePcHMGk35PfUbnpBqZtXt/JTuJ0q6q0jdfpek47rK3F7U2Un3/pKy6/AdlVmOBjCY3pXS/UTaiUKXSlocEQ90nbY9pbukM2mndH8HsAb484j4taTDaCeI6E40enaxdnolvqMyy03Vu6lJTOkeEfdExK+L/SuAPSUlT3J1oDLL0WAG03uldB+bfv0FKd2B7pTuHX8B3B0R3bP4ry+6fR9ThUeU7vpNQPKqC1uerd/WtOlJbaVKWQkhZcWF1LZSB4tTVkEA2PKRheUnjTHyrvfXb6g1oEHw6k/0RiV1d8EWRcSiwVwEFGncrwRO6tp9dkSslrQP7ZRa5wA39qvHgcosO9UGygtrIuKIcY4lp3QHkHQA8B3g3Ih4qFMgIlYXP5+S9FXaXcy+gcpdP7McDWaMaiIp3WcC3wcujoj/fv6yNE3SaPF5N+CtwP1lF+I7KrMs7fSU7hcChwCXSupkVj4JeBpYUgSpFnAr8Pmya3GgMsuNGNjqCRNI6X4FcMU41R5e9zocqMxylNdUP09KNstPflloPCnZLEcDmkLTFO76mWVpeIJQFQ5UZjkaorulKhyozLLkQGVmTTZk409VOFCZ5Sizp34OVDtBygTjlInMqW2lSk3pHhufSGhrZlpbmzYmlWt9tP5k5qfe+bbaZbY+uqp2mV5yWzPdgcosRw5UZtZswoPpZtZ8vqMys0ZzpmQzGwq+ozKzxssrTpVPSpZ0gaRlkpY9vmbtVFyTmU2IamzDwasnmOXIqyeYWaMNcIXPpnCgMstRZk/98vo2Zsb2dFkD6PpJWiDpF5JWSrq4x/HdJX2jOH6npIO6jl1S7P+FpLdUrbMXByqzLE18MF1SC7gaOBmYC5wlae6Y084H1kfEIcCnaScbpTjvTOC1wALgs5JaFevcgQOVWY4Gc0d1FLAyIh6OiGeBrwOnjjnnVOCLxedvA8cXKdpPBb4eEZsj4pfAyqK+KnXuoNYY1V33LF+jvWb+X49Do8CaOnXthHK5tpVaLte2Uss1pa1XJlzDC9x1z/Il2mvmaMXT9+iT0n1/4NGuY6uAo8eU335OkQfwSWDfYv8dY8ruX3wuq3MHtQJVRMzptV/Ssj5pocc1leVybSu1XK5tpZYbhraqiogFk1X3zuKun5mNZzXQnSPvgGJfz3MkTQNeDKztU7ZKnTtwoDKz8SwFDpV0sKTptAfHF485ZzGwsPh8OnBbRESx/8ziqeDBwKHATyvWuYNBvUe1qPyUnV4u17ZSy+XaVmq5YWhrShVjThcCS4AWcF1ErJB0GbAsIhYDXwC+JGklsI524KE475vAA8AW4G8iYitArzrLrkXt4Gdm1lzu+plZ4zlQNYikkHRen+PvLs75yRRcS6etw6ten9lkcaAaEpL2Bi4vfj1Mk59mZD7tsYX7J7kds1IOVMPjYuBlwM3A3sCrJrm9ecDPI2LzJLdjVsqBaghIOhD4APBd4LPF7tdPYnsq6r9nstowq8OBajj8A+1XST7M812xSQtUtN952RsHKmsIr0fVcJKOBN4JfCYiHizudp5icgPVvOLn8klsw6wy31E136eA9cBlAMVbvw/QJ1BJOqF4Qle23T5OFfOKn8sH9i3MJsB3VA0m6XTgjcDHgJA0szj0IHCkpL0i4ukeRX8MvKZCE8+Ms38+8MuIeKLeFZtNDgeqhirmQV1Z/Ho5z7+a0O0w4M6xOyPiGeB/JtD8PGDS39Uyq8qBqrneT/sVhL8D7h1z7DW0n/69nh6BaiIkvZT2axAeSLfGcKBqIEmjwEeBf4+Iq3ocX87zgWrQ5hc/HaisMTyY3kyfAPYE/rbXwWLsaBWTE6jmFT+XT0LdZkm8ekKDSArgXRFxw86+ll6afn2WL99RmVnjOVCZWeM5UJlZ43mMyswaz3dUZtZ4DlRm1ngOVGbWeA5UZtZ4DlRm1ngOVGbWeA5UZtZ4/w9mehBeG0SepgAAAABJRU5ErkJggg==\n", + "text/plain": [ + "
" + ] + }, + "execution_count": 5, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "result.figure(0)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Mitigation matrices\n", + "\n", + "The individual mitigation matrices can be read off the mitigator." + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[[ 1.01443299 -0.04123711]\n", + " [-0.01443299 1.04123711]]\n", + "\n", + "[[ 1.01139896 -0.04974093]\n", + " [-0.01139896 1.04974093]]\n", + "\n", + "[[ 1.01616162 -0.01818182]\n", + " [-0.01616162 1.01818182]]\n", + "\n", + "[[ 1.0060241 -0.02208835]\n", + " [-0.0060241 1.02208835]]\n", + "\n" + ] + } + ], + "source": [ + "for m in mitigator._mitigation_mats:\n", + " print(m)\n", + " print()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Mitigation Example" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "metadata": {}, + "outputs": [], + "source": [ + "qc = QuantumCircuit(num_qubits)\n", + "qc.h(0)\n", + "for i in range(1, num_qubits):\n", + " qc.cx(i - 1, i)\n", + "qc.measure_all()" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "metadata": {}, + "outputs": [], + "source": [ + "counts = backend.run(qc, shots=SHOTS, seed_simulator=42, method=\"density_matrix\").result().get_counts()\n", + "unmitigated_probs = {label: count / SHOTS for label, count in counts.items()}" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "metadata": {}, + "outputs": [], + "source": [ + "mitigated_quasi_probs = mitigator.quasi_probabilities(counts)\n", + "mitigated_stddev = mitigated_quasi_probs._stddev_upper_bound\n", + "mitigated_probs = (mitigated_quasi_probs.nearest_probability_distribution().binary_probabilities())" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Probabilities" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAApMAAAFLCAYAAACKkwciAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABIbUlEQVR4nO3deXxVxf3/8ddkJYGEJGwhYQmbyKIhEIsogtW2uHytrVAtLtRdsUIpv7rU2mqte6sVcafVaq0Lol9tsdraLypiEcsqi7KHnbAkQDAQsszvjzkXr7dJyD259yaB9/PxOI/knnPuzNw55+R+MmdmjrHWIiIiIiLiR1xTF0BEREREWi4FkyIiIiLim4JJEREREfFNwaSIiIiI+KZgUkRERER8UzApIiIiIr4lNHUBmlL79u1tXl5eUxdDRESakQULFuyy1nYIet0xISHhD8BA1Agjxx5rjNlbXV39XE1NzZNDhgw5FLrDMR1M5uXlMX/+/KYuhoiINCPGmA3BrxMSEv6QnZ3dr0OHDqVxcXGanFmOKdZaDh06lLh169YJ+/btGwz8KHQf/YclIiJSv4EdOnTYp0BSjkXGGJKTkyu7d+++Fxhe2z4KJkVEROoXp0BSjnXeNRBf67YYl0VEREREjiIKJkVERI5yF198cbebbrqpc13bb7311uyLLrqoeyzLVJfRo0fnTZw4MSfa+UyePDnn/PPP7+HnvY8++mi7IUOG9K1r+4gRI/pMnTq1XW37pqamFqxYsSKprvf27t17wMyZM9P8lKupHNMDcERERPy45hGGRDP9aZNY0JD9cnNzT9ixY0fixo0bP+vcuXNVYH2/fv36f/HFFylffPHF0r59+x566aWXNga2zZw5M+2qq67qUVxc/Flg3f333789UmU3xgxZunTpsoEDB1ZEKs2ARx99tN1Pf/rTvOTk5Jq4uDi6dOlScccdd2wZO3bs3kjn1RizZ89eXde28vLyRYHfR48enZebm3vo0Ucf3RpYt2bNmuXRLl+kqWVSRESkBcvNzT307LPPZgVef/rppykHDhw4ar/fBw0atL+8vHzR3r17F1122WW7rrzyyp47d+78r758lZWVTVG8Y9JRe7KJiIgcCy688MLdL7/8crvA6z/84Q/tLrrool3B+wRuHe/bty9uzJgxfXbu3JmYmppakJqaWlBUVJQYesv3sccea5eTk3NCRkbGoJtuuqlzbm7uCW+++WYawPvvv586aNCg49PS0gZ16NDhxHHjxnU7ePCgASgsLOwLcNJJJ/VPTU0tmDZtWibAyy+/3Pb444/vn5aWNqigoOD4efPmpQTy+vjjj1P69+/fr3Xr1gXnnntuz4qKigbFJvHx8UyYMGHXwYMH4z7//PPkyZMn55x11lk9zz///B5t2rQpmDp1avuioqLEM844o3fbtm0HdevWbeBDDz3UPjiNiooKc+655/Zs3bp1Qf/+/fvNnTv3cLluu+227K5duw5s3bp1Qa9evQa88MILGcHvtdaacePGdUtLSxvUo0ePAW+99dbhW9Pf+MY3+j788MNfyyvAGDNk2bJlyb/73e/av/XWW1lPPvlkdmpqasEZZ5zRG1xrc6Cuq6urD5cjIyNj0DnnnNOzuLg4HqC8vNycf/75PTIyMgalpaUNGjhwYL9NmzY1yR1nBZMiIiIt2CmnnLJ///798QsXLmxVVVXFW2+9lXXVVVeV1LZvenp6zYwZM1Z36NChsry8fFF5efmivLy8rzXhLViwoNXNN9/c7bnnnlu/ffv2JXv37o0vLi5ODGxPSEjgoYce2lRSUrJ4zpw5X8yZMyftwQcf7AAwf/78lQD/+c9/VpSXly+65pprSj/++OOUH//4x3lPPPHEhtLS0sVXXnnlzgsuuKD3gQMHzMGDB80PfvCD3hdddNHukpKSxWPGjCl99913MxryuSsrK5kyZUr71NTUmgEDBlQA/Otf/8oYM2ZM6d69exdde+21u8eMGdMzJyfn0LZt25a88sora+++++7cv/71r4eDvsD+Xt4lY8aM6V1RUWEAevfuXfHRRx+t3Ldv36Jbb71163XXXddjw4YNh+vhs88+a92rV6+Du3btWnLbbbdtvfTSS3sFAr2G+NnPfrbr/PPPLxk/fvz28vLyRbNmzVoTus+9997b8e2338744IMPVm7btm1JRkZG9dVXX90N4PHHH29XVlYWv2nTps9KS0sXP/nkkxtat25d09D8I0nBpIiISAt34YUX7v7jH//Y7s0330zv1avXgR49evzXU0oa6uWXX84888wz94waNWp/q1at7EMPPbTVGHN4+2mnnVZ+5plnfpmYmEjfvn0PXX755Ts/+uijOgeMPPHEEx0uu+yynWecccaXCQkJTJgwYXdiYqKdNWtW6/fff791VVWV+eUvf7kjOTnZXnHFFaUnnHBCeX3lW7JkSZu0tLRBHTt2zJ8xY0bWX/7ylzXt2rWrBhg0aNCXl1122Z74+Hi2b9+esGjRojZTp07dnJqaak855ZQDF1988a7nn3/+cCvugAEDyq+44orS5ORke8cddxQfOnTIvP/++60BrrzyytK8vLzK+Ph4rrnmmtLu3btXfPTRR60D783KyqoMlPuaa64pzcvLq5gxY0Zbv/Vem+eee67DXXfdtaVXr16VKSkp9r777tv6zjvvZFZWVpKYmGhLS0sTVqxYkZyQkMBpp51WnpWV1STBpAbgiIiItHBXX3317tNPP73vhg0bki+55JLdjUlr69atibm5uYdbK9PS0moyMjIOD+757LPPkn/yk590Xbp0aeuDBw/GVVdX079//zoDwM2bNye98cYb7Z599tmOgXVVVVVm8+bNScYY27Fjx8q4uK/atrp06VLvwJ38/Pz9CxYsWFnbtpycnMNB9MaNG5PS09OrMjMzDwdY3bt3P7Ro0aLU2vaPj4+nU6dOlZs2bUoEd6v/scce67Rly5YkgAMHDsTv3LnzcNxUW7m3bt1a5yhtP7Zt25Z0ySWX9DbGHJ7nND4+ns2bNyeOHz++ZNOmTUkXX3xxz7KysvgLLrigZMqUKVuSk5NjPieqWiZFRERauOOOO+5Qly5dDn3wwQdtL7vssj317RscmNSmc+fOlVu2bDl8O3f//v1mz549h4Oo6667rnufPn0Orl69eun+/fsX/fznP99SX3q5ubmVEydO3FZWVrY4sBw4cGDRddddV5Kbm1u5Y8eOxJqarxrUtmzZknykz1vPZzv8e7du3Q7t27cvobS09HCss3HjxqTOnTsfDpSDg7/q6mqKi4sTu3btWrlq1aqkyZMnd58yZcrG0tLSxWVlZYt79+59wNqvqq6WcicFB6fhlrc2nTp1qnzjjTdWBdddRUXFwh49elQmJyfbhx56aNvatWuXf/TRR1+89957bZ944ol29SYYJQomRUREjgLPPfdc0d///veV6enp9d7qzMnJqdq7d2/C7t27a+3fN3bs2NJZs2ZlvPfee60PHjxobrrpppzgIGr//v3x6enp1W3btq1ZtGhRq+AWR4B27dpVrVq16nBAeP311+98/vnnO86aNat1TU0N+/bti3vllVfalpaWxp155plfxsfH23vuuadjRUWFef755zM+++yzVCKgd+/elYMGDdr/k5/8pEt5ebmZN29eyssvv9z+sssuO9xyu3z58tTnn38+o7Kykt/85jedkpKS7De/+c0vy8rK4owxZGdnVwJMmTKl3Zo1a1KC0y8pKUkMlPvZZ5/NXLduXcro0aPDmqKoY8eOlevXr68zeL7iiit23H777V1WrVqVBLB169aEF198MQPgb3/7W9qnn36aUlVVRUZGRnVCQoJtqic1KZgUERE5CgwYMKBixIgR9fY3BCgoKDh43nnnlfTq1euEtLS0QUVFRYnB2wsLCw/ed999G8eNG9czOzs7v02bNjVZWVlVrVq1sgAPPvjgptdffz2rTZs2BVdffXX3733ve18b7HPzzTdvve666/LS0tIG/eEPf8gcMWJE+WOPPVY0ceLEbm3bth3Uq1evgYF+i61atbKvvvrq2pdeeql9VlbWoOnTp2eNGjVqT6TqZPr06es2bdqU1Llz5/wxY8b0uuWWW7Z+73vfKwts/9a3vrVn+vTpWRkZGQWvvvpqu1dffXVtcnKyHTJkyMFrr722eMSIEf06dOiQv3Tp0pSCgoL9wWmfeOKJX65evbpV+/bt8++6667cF154YW12dnZ1OOUbP378rtWrV6ekpaUN+ta3vtUrdPvtt9++45xzztnzne9857jWrVsXDB069PhPPvmkNbjuCBdeeGGvtLS0gv79+w8cNmxY2Q033NCoLg5+meD/No41hYWFdv78+U1dDBERaUaMMQustYWB10uWLCnKz8/fVd97jmZ79+6Na9euXcGyZcuWHn/88b4H9kjLt2TJkvb5+fl5oevVMikiIiJf89JLL7UtKyuL27dvX9z48eO79OnT58Bxxx2nQFJqpWBSREREvuatt97KyMnJOTE3N/fE9evXt3rllVfWBo9cFgmmqYFERETka1599dUNwIamLoe0DPo3Q0RERER8UzApIiIiIr4pmBQRERER3xRMioiIiIhvCiZFRERExDcFkyIiItIgvXv3HjBz5sy0uraPGDGiz9SpU5vk+dChcnNzT3jzzTfrLGukfOMb3+j78MMPt/fz3tGjR+dNnDgxp67tqampBStWrEgK3ffdd99tk5eXN7Cu961evTopNTW1oKqqyk+xwqapgURERMK0tu+dQ6KZfq+Vdy5oyH7GmCFLly5dNnDgwIrAusmTJ+esXbs2+a233lof6XKtWbNmeX35zJ49e3Uk8pk5c2baVVdd1aO4uPizSKQXavTo0Xl//etfsxITE21iYqIdMGDAl48//vimgoKCg9HIz6/y8vJFta0/66yz9hcVFS0LvM7NzT3h8ccfLwo8KrJPnz6H6npvNKhlUkRERI4548eP315eXr5o8+bNn7Vv377q8ssvzwvdp6amhurqsB63fUxSMCkiInKUmjlzZlqnTp1OvOOOOzplZWXld+jQ4cQpU6Ycvg09evTovEsvvbTbiBEj+qSmphYMHjz4+I0bNyZceeWVXdPT0wf16NFjwMcff5wS2D9w63jGjBnpU6dOzX777bczU1NTC/r27dsfvn7Lt6qqimuuuaZLZmZmfm5u7gn33ntvB2PMkMrKSgCmTJnSrmfPngNat25d0KVLlxN++9vftgfYt29f3JgxY/rs3LkzMTU1tSA1NbWgqKgosbq6mttuuy27a9euAzMyMgadc845PYuLi+MDZXv88cezcnJyTsjIyBh0yy23ZDe0jtLS0mouvvjiktWrV6cEPsOECRNyBw8efHxqaurgzz//PPm9995rPXDgwH5paWmDBg4c2O+9995rHZzG2rVrk0844YR+bdq0KTjzzDN7BZfr7LPP7tm+ffv8tLS0QYWFhX3nz5/fKvi9u3btSjjllFP6tG7duuCkk07qu2rVqqTANmPMkGXLliXXdVwBvve97/XYtm1b0g9/+MM+qampBbfffnunlStXJgXX9e7du+MvvPDC7h06dDixY8eOJ06cODEncAt82bJlySeddFLftLS0QZmZmfnnnntuz4bWXYCCSRERkaPY7t27E/fu3Ru/bdu2zx577LENt956a7edO3ceDnbefvvtzHvuuWfLrl27FiclJdUMGzas3+DBg8tLSkoWn3feeaWTJ0/uGprmmDFj9k2YMGH7ueeeW1peXr5o5cqVK0L3efjhhzvMmjWr7fz581csXrx4xcyZMzODt3fq1Knqb3/725qysrJFTz/99Ppf/epXXefMmZOanp5eM2PGjNUdOnSoLC8vX1ReXr4oLy+v8t577+349ttvZ3zwwQcrt23btiQjI6P66quv7gawYMGCVjfddFP3P/7xj+u3bdu2ZPfu3QnFxcVJoWWqzd69e+P+8pe/ZPXr1688sG7GjBlZzzzzTFFZWdnCtm3bVo8ePbrP+PHji0tKShZPmDChePTo0X22b99+uA5fe+21ds8+++z6rVu3LklISODaa6/tFtg2atSovatXr166Y8eOJSeeeGL5pZde+rVg7a233mr3y1/+ctuuXbsWDxw4sHzs2LE9GlLugDfffHN9586dD73yyiury8vLF919993Fofv88Ic/zEtISGDt2rXLFi1atOL9999v+/vf/749wM9//vOcM844Y++ePXsWb9my5bOJEyfuCCd/UDApIiJyVEtISLC//e1vtyYnJ9uLLrpob0pKSs1nn312uHVs1KhRe0477bTy1NRUe9555+1JTk6uufHGG3cnJCRw6aWXlq5YsSLVT75vvPFG5vXXX1/cq1evyg4dOlTffPPN24K3//CHP9w7YMCAiri4OM4999z9p5566r7333+/TV3pPffccx3uuuuuLb169apMSUmx991339Z33nkns7KykpdffjnzjDPO2Hv22WfvT0lJsQ8//PBWY4ytr3xPP/10dlpa2qBevXqd8OWXX8a/8MILRYFtF1100e7CwsKDiYmJ/PWvf03v3r17xY9//OOSxMRErrvuupKePXsenD59ekZg/zFjxuw+6aSTDqanp9fce++9W/7+979nBlr+Jk2atDszM7MmJSXFPvjgg1tXrlyZsnv37sOB6De/+c3D5X7kkUe2LF68uM2aNWsSw67wOmzatCnhww8/bPvMM89sTE9Pr8nNza268cYbi2fMmJEF7vzYuHFjclFRUWJqaqodNWrU/nDz0AAcERGRFio+Pp5Dhw6Z4HWVlZUmMTHxcCDVtm3bqsTEr2KTlJSUmrKyssONSR07dqwM3ta+ffvDQ4BTU1NrDhw4cDjwCUdxcXFit27dDqfdo0ePQ8Hbp0+fnn7PPffkFBUVtaqpqeHgwYNxAwYMOFBXetu2bUu65JJLegcHifHx8WzevDlx69atibm5uYfTT09Pr8nIyKh3KPN11123/dFHH91a27auXbseTmvr1q1JXbp0qQje3qVLl0NbtmxJrG3/Pn36HKqqqjLbtm1L6Ny5c9XEiRNz//a3v2WWlpYmBsq+ffv2hHbt2lUDBJe7bdu2Nenp6VUbN25M6t27dyURsGbNmqSqqirTuXPn/MA6a63Jzs4+BDBlypTNN998c+6wYcP6paenV994443bJ02atDucPNQyKSIi0kJlZ2cfWrNmzddu5xYVFSV169btUF3viZQjtfx17NixctOmTYcDrvXr1x8u54EDB8yPfvSjXpMmTSresWPHkrKyssUjR47ca62tM+1OnTpVvvHGG6vKysoWB5aKioqFPXr0qOzcuXPlli1bDqdfVlYWt2fPHt8NZsZ8FZ/n5OQc2rx589f6LW7ZsiUpNzf3cLC3adOmw3mvWbMmKSEhwXbu3Lnq6aefznr33Xcz3nvvvVX79u1btH79+qUAgc8ZSCvw+969e+P27duXEMnj17Nnz8qkpCRbUlJyuN7279+/KDAyv1u3blWvvPLKhh07dnz2+OOPb7jlllu619ZPsz4KJkVERFqo888/v+T+++/PWbt2bWJ1dTVvvvlm2qxZszLGjh1bEu28O3XqVLV58+akukY7X3DBBaVPPfVUp/Xr1yfu2rUr/sEHHzw8KObgwYPm0KFDcR07dqxMTEy006dPT//444/TA9tzcnKq9u7dmxB8O/iKK67Ycfvtt3cJDFDZunVrwosvvpgBMHbs2NJZs2a1/cc//tHm4MGD5v/9v/+XY639WoutX6NHj95bVFSU/NRTT2VVVlYybdq0zDVr1rT6wQ9+sDewz+uvv95uwYIFrcrKyuJ+8Ytf5Jx11lmlCQkJlJWVxSclJdmOHTtW7d+/P27SpEm5oel/8MEHh8s9efLk3Pz8/C/DbZVs37595Zo1a2oNALt371556qmn7r322mu7lpSUxFVXV7N8+fLkt99+uw3As88+m7l27dpEgHbt2lUZY4iLi6v3H4VQCiZFRERaqAceeGDrSSedtH/EiBHHZ2RkDLrtttu6PPPMM+tOOumkqM+XOG7cuBKAzMzMQf379+8Xun3y5Mk7R44cua+goGBAfn5+/1GjRu2Nj4+38fHxZGZm1tx9990bx40b16tt27aDXnrppXZnnnnm4eCsoKDg4HnnnVfSq1evE9LS0gYVFRUl3n777TvOOeecPd/5zneOa926dcHQoUOP/+STT1oDFBYWHnzggQc2Xn755T2ys7PzMzMzqzp16hSR1r3s7OzqGTNmrJk6dWqnrKysQY888kj2jBkz1nTu3PnwbfQxY8bsvvzyy3t07tw5v6KiIu6ZZ57ZBDB+/Pjdubm5FV27ds0//vjjB5x88slfhqb/3e9+d/evf/3rzllZWYOWLFmS+tJLL60Lt4w33XTT9oceeqhzWlraoF/96ledQrdPnz696NChQ6Zfv34DMzIyBo0ZM6ZX4Db9p59+2nrYsGH9UlNTC77//e/3vvvuuzf2798/rLozwU2tx5rCwkI7f/78pi6GiIg0I8aYBdbawsDrJUuWFOXn5+9qyjIdDaZPn54+adKk7lu3bl3a1GURf5YsWdI+Pz8/L3S9WiZFREQk4vbv329effXVtpWVlaxfvz7xnnvuyTnrrLP2NHW5JPIUTIqIiEjEWWvNb37zm5yMjIyCIUOG9O/Tp8/B3/3ud1uaulwSeZoaSERERCIuLS2tZtmyZZ83dTkk+tQyKSIiIiK+KZgUERGpX01NTU1EppkRaam8a6DWeaAUTIqIiNRv2c6dO9sqoJRjkbWWioqKxA0bNmQAc2rbR30mRURE6lFVVXX19u3b/7B9+/aBqBFGjj01xpi91dXVj9bU1DxZ2w4KJkVEROoxZMiQHcB3m7ocIs2V/sMSEREREd8UTIqIiIiIb7rNHQPXPFL/9mmTYlEKERERkchTy6SIiIiI+KZgUkRERER8UzApIiIiIr4pmBQRERER3xRMioiIiIhvCiZFRERExDcFkyIiIiLim4JJEREREfEt5sGkMeYGY8x6Y8xBY8wCY8xpDXzfcGNMlTFmWcj6y40xtpalVXQ+gYiIiIgExPQJOMaYi4ApwA3AHO/nO8aY/tbajfW8LxN4Afg/ILeWXcqBXsErrLUHI1XuI1nb9876dxh/hO0iIiIiLVSsWyYnA3+y1k6z1n5urZ0AbAPGH+F9fwSeB+bWsd1aa7cHLxEss4iIiIjUIWbBpDEmCRgC/DNk0z+BU+p53w1AJ+DuepJPMcZsMMZsNsbMNMYUNLrAIiIiInJEsWyZbA/EA8Uh64uB7NreYIw5AbgDuNRaW11HuiuBK4HzgbHAQeBjY0yfSBRaREREROoW0z6T4TDGJAOvAj+z1q6vaz9r7VyCbn8bY/4NLAYmABNrSfda4FqAnJwcPvjgAwB69uxJWloaS5YsAaBdu3YMGDCA2bNnA5CQkMDw4cNZuHAh+/btA6CwsJDi4tDYOHyrV69my5YtAPTt25f4+HhWrFgBQHZ2Nj169GDuXPcRU1JSGDp0KPPmzePAgQMADBs2jPXr17N9u7u7379/f6qrq1m5ciUAubm5dOnShXnz5gHQpk0bCgsLmTt3LhUVFQAMHz6cVatWsWPHDgAGDhxIRUUFq1evBqBr16506tSJ+fPnA5Cens7gwYOZM2cOVVVVAIwYMYLly5eze/duAPLz8ykrK2PdunUA5OXlkZWVxcKFCwHIzMwkPz+fDz/8EGstxhhGjhzJkiVLKC0tBWDw4MGUlJRQVFTU6OO0adMmAPr06UNycjLLlrmxXB07duS4445jzpw5ACQnJzNs2DDmz5/P/v37ARg6dCibN2/WcdJx0nE6Bo6TiITHWGtjk5G7zV0OjLXWvha0/nFgoLV2ZMj+ecB6ILhFMg4w3rpzrLWht8wD730OyLbWnl1fmQoLC23gj3ljHGkAzv1HGIAzbVKjiyAiIhFijFlgrS1s6nKItBQxu81trT0ELAC+HbLp28C/a3nLFuAEYFDQ8hSwxvu9tvdgjDHAibiBPSIiIiISRbG+zf0w8GdjzKfAx8D1QA4uSMQY8wKAtXactbYSCJ1TcgdQYa1dFrTuDuATYDWQjru1fSJHHiEuIiIiIo0U02DSWvuqMaYdcDvQGRcsnmOt3eDt0s1HshnAM7hBPHuBRcAIa+2njS+xiIiIiNQn5gNwrLVPAE/Use30I7z3TuDOkHU/BX4amdK1TNc8Uv929ckUERGRaNGzuUVERETENwWTIiIiIuKbgkkRERER8U3BpIiIiIj4pmBSRERERHxTMCkiIiIivimYFBERERHfFEyKiIiIiG8KJkVERETENwWTIiIiIuKbgkkRERER8U3BpIiIiIj4pmBSRERERHxTMCkiIiIivimYFBERERHfFEyKiIiIiG8KJkVERETENwWTIiIiIuKbgkkRERER8U3BpIiIiIj4pmBSRERERHxTMCkiIiIivimYFBERERHfEpq6AHJka/veWf8O44+wXURERCRK1DIpIiIiIr4pmBQRERER3xRMioiIiIhvCiZFRERExDcFkyIiIiLim4JJEREREfFNwaSIiIiI+BZWMGmMiTPGxAW9zjbGXG2MOTXyRRMRERGR5i7clsm3gQkAxpg2wHzgt8AHxphxES6biIiIiDRz4QaThcAs7/cLgH1AR+Aa4GcRLJeIiIiItADhBpNtgD3e798B/tdaW4kLMHtFsFwiIiIi0gKEG0xuBE41xrQGRgHveeuzgPJIFkxEREREmr+EMPd/GPgzsB/YAMz21o8AlkawXCIiIiLSAoQVTFprnzbGLAC6Au9Za2u8TWuBX0a6cCIiIiLSvIXbMom1dj5uFHfwurcjViIRERERaTHCnrTcGHODMWa5MabcGNPTW3eLMebCyBdPRERERJqzcCctnwTcDjwDmKBNW4EbI1csEREREWkJwm2ZvB64xlo7BagKWr8QGBCxUomIiIhIixBuMNkdWFbL+kogpfHFEREREZGWJNxgch0wuJb15wArGl8cEREREWlJwh3N/TvgMWNMKq7P5DBjzGXAzcCVkS6ciIiIiDRv4c4z+ZwxJgG4F0jFTWC+FZhorX01CuUTERERkWbMzzyT04Bpxpj2QJy1dkfkiyUiIiIiLUHYwWSAtXZXJAsiIiIiIi3PEYNJY8xnwEhrbakxZilg69rXWntiJAsnIiIiIs1bQ1omXwcqgn6vM5gUERERkWPLEYNJa+2vg36/M6qlEREREZEWJdzHKc4yxmTUsj7dGDMrYqUSERERkRYh3EnLTweSalnfCjit0aURERERkRalQaO5jTHBT7050RhTEvQ6HhgFbIlkwURERESk+Wvo1EDzcQNvLPDPWrYfACZEqlAiIiIi0jI0NJjsgXt84jrgG8DOoG2HgB3W2uoIl01EREREmrkGBZPW2g3er+H2sRQRERGRo1hDJi2/APibtbbS+71O1to3IlYyEREREWn2GtIyOQPIBnZ4v9fF4gbjiIiIiMgxoiGTlsfV9ruIiIiISMyDQ2PMDcaY9caYg8aYBcaYOuenNMaMNMb82xiz2xhzwBjzhTHmZ7XsN9oYs8IYU+H9/H50P4WIiIiIQMP7TDbIkfpMGmMuAqYANwBzvJ/vGGP6W2s31vKW/cCjwFKgHDgVeNoYU26tfcJLcxjwKnAH8AZwAfCaMeZUa+28hpZdRERERMLX0D6TDdGQPpOTgT9Za6d5rycYY84CxgM//68ErV0ALAhatd4Lbk8DnvDWTQLet9be472+xxjzTW/92AaWXURERER8OOJtbmttXAOXegNJY0wSMIT/nvT8n8ApDSmsMabA2/fDoNXDaknzHw1NU0RERET8i2Wfyfa4lsvikPXFuNHidTLGbDbGVOCexPOEtfapoM3ZftIUERERkcZrKfNMnga0AU4GHjDGrLfW/tlPQsaYa4FrAXJycvjggw8A6NmzJ2lpaSxZsgSAdu3aMWDAAGbPng1AQkICw4cPZ+HChezbtw+AwsJCiotD49jwrV69mi1b3KPN+/btS3x8PCtWrAAgOzubxAik36VLF+bNc11I27RpQ2FhIXPnzqWiogKA4cOHs2rVKnbs2AHAwIEDqaioYPXq1QB07dqVTp06MX/+fADS09MZPHgwc+bMoaqqCoARI0awfPlydu/eDUB+fj5lZWWsW7cOgLy8PLKysli4cCEAmZmZ5Ofn8+GHH2KtxRjDyJEjWbJkCaWlpQAMHjyYkpISioqKgMYdp02bNgHQp08fkpOTWbZsGQAdO3bkuOOOY86cOQAkJyczbNgw5s+fz/79+wEYOnQomzdvrvc49ejRg7lz5wKQkpLC0KFDmTdvHgcOHABg2LBhrF+/nu3btwPQv39/qqurWblyJQC5ubk6TjpOOk7N4DiJSHiMtbb+HYypAbKttTu83+ti67vV7d3mLgfGWmtfC1r/ODDQWjuyQQU25nbgCmttL+/1RmCqtfa3QfvcBNxore1eX1qFhYU28Me8Mdb2vbPe7fePr3/7tElNm76IiHzFGLPAWlvY1OUQaSka2mdyR9DvvvpMWmsP4QbTfDtk07eBf4dZ5uSg13MjkKaIiIiI+NCgZ3NH0MPAn40xnwIfA9cDOcBTAMaYFwCsteO81xOA9cBK7/0jgJ/x1UhucFMNzTbG3Aq8CXwf+CYwPMqfRUREROSYF3YwaYwZjJt2p7+36nPg99bahUd6r7X2VWNMO+B2oDOwDDjHWrvB26VbyFvigQeAPKAKWAvcihd8emn+2xjzQ+Bu4C5vn4s0x6SIiIhI9IUVTBpjLgFeAGYBf/dWnwx8aoy53Fr74pHS8CYbf6KObaeHvH4EeKQBac6g4fNhioiIiEiEhNsyeQ/wS2vtvcErjTE/x7UMHjGYFBEREZGjR7jzTHYAptey/jWgY+OLIyIiIiItSbjB5PvA6bWsP52vP5VGRERERI4BDZ20POAd4D5jTCHwibfuZOAC4M6Il05EREREmrWG9JmsbWDL4afIBJlKHQNrREREROTodMRg0loby+d3i4iIiEgLokBRRERERHzzM2l5JnA2boLxpOBt1tq7IlQuEREREWkBwp20/GTgbaACN03QFtyTbCqAItwTaERERETkGBHube7fAn8BcoGDwBm4Fsr5uMceioiIiMgxJNxg8kTgMWutBaqBZGttMXALmhpIRERE5JgTbjB5KOj3YqC79/t+ICciJRIRERGRFiPcATgLgZOAVcAHwN3GmE7ApcBnkS2aiIiIiDR34bZM/gLY6v1+O7ATN1l5Jv89ibmIiIiIHOXCapm01s4P+n0nboogERERETlGhT3PJIAxphfQz3u5wlq7LnJFEhEREZGWItx5JtsBfwS+C9R8tdrMBK601u6OcPlEREREpBkLt8/kH4DewGlAK28ZAfQApkW2aCIiIiLS3IV7m3sUcKa1dm7Quo+NMdcB/4pcsURERESkJQi3ZXIn8GUt68sB3eIWEREROcaEG0zeBTxijMkNrPB+fwg9l1tERETkmHPE29zGmKWADVrVAygyxmzxXgee090R16dSRERERI4RDekzOSPqpRARERGRFumIwaS19texKIiIiIiItDx+Jy0/A+iPu/293Fr7QSQLJSIiIiItQ7iTlucC/wsM4atndOcYY+YD37fWbq3zzSIiIiJy1Al3NPejQDXQ21rb1VrbFejjrXs00oUTERERkeYt3Nvc3wZOt9auD6yw1q4zxkwE/i+iJRMRERGRZi/clkn4+jRB9a0TERERkaNcuMHk/wFTjTFdAyuMMd2AR1DLpIiIiMgxJ9xgciLQGlhnjNlgjNkArPXWTYx04URERESkeQu3z+Ru4BvA6cDx3rrPrbX/imShRERERKRlaHAwaYyJB/YC+dba94D3olYqEREREWkRGnyb21pbDWwAkqJXHBERERFpScLtM/kb4H5jTPtoFEZEREREWpZw+0z+DOgBbDHGbAa+DN5orT0xUgUTERERkeYv3GByBm5OSROFsoiIiIhIC9OgYNIYkwr8FvgekIibU3KCtXZX9IomIiIiIs1dQ/tM/hq4HHgbeBn4FvBklMokIiIiIi1EQ29zXwBcZa19BcAY8xfgY2NMvDfKW0RERESOQQ1tmewKfBR4Ya39FKgCcqJRKBERERFpGRoaTMYDh0LWVRH+AB4REREROYo0NBg0wIvGmIqgda2AacaY8sAKa+13I1k4EREREWneGhpMPl/LuhcjWRARERERaXkaFExaa6+IdkFEREREpOUJ93GKIiIiIiKHKZgUEREREd8UTIqIiIiIbwomRURERMQ3BZMiIiIi4puCSRERERHxTcGkiIiIiPimYFJEREREfFMwKSIiIiK+KZgUEREREd8UTIqIiIiIbwomRURERMQ3BZMiIiIi4puCSRERERHxTcGkiIiIiPgW82DSGHODMWa9MeagMWaBMea0evbtbIx5yRjzhTGm2hjzp1r2udwYY2tZWkX1g4iIiIhIbINJY8xFwBTgXqAA+DfwjjGmWx1vSQZ2AfcD8+pJuhzoHLxYaw9GqtwiIiIiUrtYt0xOBv5krZ1mrf3cWjsB2AaMr21na22RtXaitfZPQEk96Vpr7fbgJfJFFxEREZFQMQsmjTFJwBDgnyGb/gmc0sjkU4wxG4wxm40xM40xBY1MT0REREQaIJYtk+2BeKA4ZH0xkN2IdFcCVwLnA2OBg8DHxpg+jUhTRERERBogoakL0FjW2rnA3MBrY8y/gcXABGBi6P7GmGuBawFycnL44IMPAOjZsydpaWksWbIEgHbt2jFgwABmz54NQEJCAsOHD2fhwoXs27cPgMLCQoqLQ2Pj8K1evZotW7YA0LdvX+Lj41mxYgUA2dnZJEYg/S5dujBvnut22qZNGwoLC5k7dy4VFRUADB8+nFWrVrFjxw4ABg4cSEVFBatXrwaga9eudOrUifnz5wOQnp7O4MGDmTNnDlVVVQCMGDGC5cuXs3v3bgDy8/MpKytj3bp1AOTl5ZGVlcXChQsByMzMJD8/nw8//BBrLcYYRo4cyZIlSygtLQVg8ODBlJSUUFRUBDTuOG3atAmAPn36kJyczLJlywDo2LEjxx13HHPmzAEgOTmZYcOGMX/+fPbv3w/A0KFD2bx5c73HqUePHsyd607FlJQUhg4dyrx58zhw4AAAw4YNY/369Wzf7nph9O/fn+rqalauXAlAbm6ujpOOk45TMzhOIhIeY62NTUbuNnc5MNZa+1rQ+seBgdbakUd4/0xgl7X28gbk9RyQba09u779CgsLbeCPeWOs7XtnvdvvH1//9mmTmjZ9ERH5ijFmgbW2sKnLIdJSxOw2t7X2ELAA+HbIpm/jRnVHhDHGACfiBvaIiIiISBTF+jb3w8CfjTGfAh8D1wM5wFMAxpgXAKy14wJvMMYM8n5NB2q814estSu87XcAnwCrvX0m4oLJWkeIi4iIiEjkxDSYtNa+aoxpB9yOmw9yGXCOtXaDt0tt800uCnl9HrAByPNeZwDP4Abx7PX2H2Gt/TSihRcRERGR/xLzATjW2ieAJ+rYdnot68wR0vsp8NOIFE5EREREwqJnc4uIiIiIbwomRURERMQ3BZMiIiIi4puCSRERERHxTcGkiIiIiPimYFJEREREfFMwKSIiIiK+KZgUEREREd8UTIqIiIiIbwomRURERMQ3BZMiIiIi4puCSRERERHxTcGkiIiIiPimYFJEREREfFMwKSIiIiK+KZgUEREREd8UTIqIiIiIbwomRURERMQ3BZMiIiIi4puCSRERERHxTcGkiIiIiPimYFJEREREfFMwKSIiIiK+KZgUEREREd8UTIqIiIiIbwomRURERMQ3BZMiIiIi4puCSRERERHxTcGkiIiIiPimYFJEREREfFMwKSIiIiK+KZgUEREREd8UTIqIiIiIbwomRURERMQ3BZMiIiIi4puCSRERERHxTcGkiIiIiPimYFJEREREfFMwKSIiIiK+KZgUEREREd8UTIqIiIiIbwomRURERMQ3BZMiIiIi4puCSRERERHxTcGkiIiIiPimYFJEREREfFMwKSIiIiK+KZgUEREREd8UTIqIiIiIbwlNXQBpemv73lnv9l4r698uIiIixy61TIqIiIiIb2qZlCO65pEj7zNtUrRLISIiIs2RWiZFRERExDcFkyIiIiLim4JJEREREfFNwaSIiIiI+KYBOBITmn5IRETk6KRgUo4K0Q5WFQyLiIjUTsGkSDPQ1MFwLPK4f3z92zW9lIhIy6Q+kyIiIiLiW8xbJo0xNwA3AZ2B5cAka+1H9ew/EngYGABsBR601j7VmDSl+TnSxOhqtZJot3yqK4OIiD8xDSaNMRcBU4AbgDnez3eMMf2ttRtr2b8H8HfgWeBSYDjwhDFmp7X2dT9pyrFJwaqIiEh0xLplcjLwJ2vtNO/1BGPMWcB44Oe17H89sNVaO8F7/bkxZijwM+B1n2mKSC0UcNevsfWjlk8ROVrFrM+kMSYJGAL8M2TTP4FT6njbsFr2/wdQaIxJ9JmmiIiIiERILFsm2wPxQHHI+mLgW3W8Jxv4Vy37J3jpGR9piogcdZrDiP3mnr6IRIex1sYmI2NygC3ASGvt7KD1vwIusdb2reU9q4AXrbV3Ba0bAXwI5OCCyXDTvBa41nvZF1gZgY93JO2BXUq/SfNQ+k2fh9Jv2vRjkUdLTz+gu7W2QwzyETkqxLJlchdQDXQKWd8J2F7He7bXsX+Vl54JN01r7TPAMw0udQQYY+ZbawuVftPlofSbPg+l37TpxyKPlp6+iPgTsz6T1tpDwALg2yGbvg38u463za1j//nW2kqfaYqIiIhIhMR6NPfDwJ+NMZ8CH+NGa+cATwEYY14AsNaO8/Z/CrjRGPMI8DRwKnA5MLahaYqIiIhI9MQ0mLTWvmqMaQfcjptgfBlwjrV2g7dLt5D91xtjzgF+j5vqZyswMTDHZAPTbA6ifVu9pacfizyUftPnofSbNv1Y5NHS0xcRH2I2AEdEREREjj56NreIiIiI+KZgUkRERER8UzApIiIiIr4pmIwBY4wJ+l11HmMh9W/q21ekqeg8PTL9/RRpnnRhxoC11nrPEcdaWxNYry8MJ9pfEF79pwd+j2ZeLZUxJr6ln4/GmPgop2+Cf0ZaSz9Po13/4P5+GmM6efm1MsbEeno7EamFRnNHmTEmGxgDDMY9vvETYIa1dm4E80iw1lZFKr068oi31lZHOQ+DOydratvm5wvWGNMbNy/pN4E83ET4fwPet9YWNybtkHxiUT+dgTTgAO6Z9NuttQcjnEcc7hhE/LMYY9KA1sAOIBUor+1YRyCfOPj6P25B23wfay9waW2t3dvIItaWdqzO06gfgyjW/yBgHHAOkA3MB94D/g9YZK2tjkQdiUj4FExGmTHmbaA38DlunsxTgIHAGuB+4M+R+uIOtAxEM6jx8rCR+gIyxvwYWA7Ms9YeCFof5+XT2C/PD3EB2BzcIzbPAIYDu4EpwEOR/BKKdP0EpXsDcCXu3KnEPfnpE2AW8KG1tqIRAfdzXlrTrbWlQesTgJpIfBZjzMXAFbh/qsAFS+8A/7LWrvT2aUygcQ+wEPintbYsaH087jM09jw6yyt/PpCEC2DewgV7XzYmbS/9qJ+n0TwG0a5/L60FQBkuyC4GzvWWauBPwO3W2i8VUIo0AWutligtuC+EnUAX73UroC0uoJwGrAMmNyL9U4ElwNVAUsi2BFw3BgNk4f3j4COPQmAmrnU1sZY8fKXrvX84UAPMBl4EJgL5IfskA78EOvtI/5te/WeGrM8B7gC2AE8A8Y34DFGrn6B0TvfK+gDQDzjbO39WeevvARIaeQxWARuB14DzQ/ZJAZ4F+vvMYwRQBLzglf1HuCCmAtgE/LiR9RP4DIuAj4DfASNr+QyPAt18pr8S+Cdwg3c+zsMFMcuAHzSy/LE4T6N2DKJd/0HXwE6gVS3brgI2eJ8nrTHHQosWLf6WJi/A0bwAdwKz6tiWDvwC2A8M9pn+894X2jagCngX+J+QfU711vsNNp4HDnnBxlrgD7V8UZwCvBLulx3uyUb/Bu7CBWT/Ad7HPTpzHO5231Dvi6qNj7L/DNfiluK9jgfigrZfCuwFzmjEMY5a/QS9/yXg6VrWJ+IeH7oD+KPPtO8G/gGc79XXTC/oWAU8DpwMnOQdA19f1MB04Jla1qd618hu4JeNOAYPeufN9cBjwAfAYu9z/RzXmvsNv58BmAFMq2X98bgnsmwBxjWi/LE4T6N2DKJd/14e13tpdvJeJxP0DzQwEtda+X2/daRFixb/izovR9e/gB8bY86y1r4bvMFau88Y8wDuj+BI3C2icOXhHiM5E/fH+gfAa8aYSlwL06PAD4Ec679P5XG4lq/5Xh4jgBeNMaXAX4E/A5cAA234t9fbA/Ottb/yBiidAXwbKMDdirsQ94X9f9ba/T7K/nfgFuAC4C+B8gVug1lrXzTG/ABX/7N8pA/RrZ+AQ0CWMaaVtfagMaYVUGWtrQSe8o73TcaYAdba5WGm3Rr3z8jb1toqY8zruC//YbgWp1eAXOAdG3T7MkxJuH96ADDGJONufZYDd3pdGn5ojPmztbbIR/rtgXXW2qe8tAZ75T8J+C7wfaAH8K7Pz5CF66YSKH/g9v8XxpiJuKB+kjHmXWvtDh/px+I8jeYxiHb9g/sb9wvctfSwtbbC+xyB2+gfel0FTgP+12ceIuJXU0ezR/OCu7XzAq5/5K3AEFwH/sD2DGAzcIGPtHNwrWDXea/jgUzcH/Bbca18FbjWgPN8lj8P90d8vPe6FW4Q0YW4flz/wd0iqwG+6yP9gcDZtazviBuM8Ccv7XN8lj8eeAgowbUgnQO0C9qe7ZV/THOsn6B8RgG7gAtD1id4P1vjblGP9JF2G+DkOs7d/rguFL6PgZfWJbiWr1NCj4/3MwtYDwzzmX42cHot69vibiHf1cjzaCKuH+NxIesDfc67ete43zsMgfN0dzTO02gfgxjUf6Ceb8HdyfkXrv9wTtA+Pb06alSXAy1atPhbNAAnyowxXXG3er4FlOICjO24VoJhQF9rbV+fabcDUq21m0LWJ+ACywnARGttRiPK3wkXAK8LWZ+G+wN+A3CR3zwCI9G9Fo04ggZ8GGPOA1601rZtRPlTgfG41pEUXNBVgrttOBR3a3Fw3SkcMf1o14/B3dK7F3c85+Nurb9mrd1tjGmPa/15xFqb7vdzBOdng/4oGGO+ixuY06oR5W+N6+N5Nq4V7nVca/Me77hfhLsFm9bY8nt5fm3wlncevWytbeMzvQ64rgY9cS21/8CNHi7ztl8APN+Y8nsthRNwxzIF909mRM7TkGNwFq5vYUSOQS3nS+A6ro5U/Yfk931cYNwT9zd0L65lvQAosdaOaGweIhI+BZMxYow5Edf36RTcwIxMXH/Bh621S6OU55u4P+qjI5hm6JfHm0CFtfaiSOaBGzg0A0i31n4rAml2A87D9QFsD3TCtXA8aa1d39j0vTzicddUVdC6N4lQ/Rhj/ge4GPfF2QHXR6wKF2z+wVr7u8bmEZKfAX6N66d2XSPTao0bSXw+ruW5ChcIxOOOxSvW2jsbVeDa843DDV7JstZe2Ih0euP+MRiJK/tmoBwXpB0PzLTW3hyB8h4P/A8wCPc3ojMROk+9f3Aux42A7oRrLYzaMQi6jhtd/yHpdsX1BR8AdPGWd3EzY/jpZiAijaRgMgq8VoYTcX2gyoDPcH0Dt3vb++IGORyyUToAxpg2wFRgirV2cZTyyADeAG6x1v4njPcF6uf7uC+z5bjbhJusm9ojzrrJiROAttba3WGWqw2u7+JYYA9uJO4C3DGoNMZ0sNbuDCfNI+SXZkP6gnlBTDo+6qeW9A8H8MbNNdkf6Ibrh9YKN9J6tfXRJ/NI0/94n6N16OfzkU+yddMXZeOOzQDc7eFk3KCNBdbaQz7Trnc6Jm97mrV2j6/Cfz2t/rh/Svrjbg2n4AaSvW9d/0M/aQYG3FSGrI/oeRqUbl/cXZEeuECsFT6OQch1VgKsxg3cWmqt3Ra0X6Pr3ztPsSF9vwPnld90RSQyFExGgTHmUVwguQ33hZOHu7X9FvCAtXZDjMrh6w9tSPBypGAjNdwv0VrqpztuDs43cbdr19X97gal/zzuS261l34XYB9uNOjj1toPG5O+l0c/YDKulXAN7vb5YuCj4G4HxpgUGzR/ZiPyi+jcecaY4dbaOSHrvnasA0F9I/IIrqN1uD55HwOzIxTYDbHWLghZF8l5Dbvi+uZ9A1f25cC/rbVLvHxSrL+BYYH0ayt/Ei4wrqzjbeGkH3odfy0Yi8Dxre06K8VdB89Ya9/3X/rDedR2nn6tjkwMHtogIkdQX4dKLeEvuBaLfcB38DrR425J3or7QjoAXIf/eR87AZfhbhvVt1+rwPeIjzw68N/T2xhC5k3kq877Dc6jgfVzrZefn7L3x7UGnxxUvrbANXw1N+CduFt7fo9BL1xr52zgPtzo0bm4EfnTge+E1l2kjjPe3KGB+geSfaR9PO4WZxmuD+CptRzrJFwLXKcI1tEnuHlRXwe+5SfdoPT7eJ9hGfAwUFDLZ0jEBYJJPtLv4R3Pz3F9VBfipgBajpu6Ki8G5U/yW34vjdqu4zhca3BCUD5hT1nVwOvsDkKmOYrweRr4LL7PUy1atERmafICHG0LcBuu5SXwOiFk+73AFwSNRAwz/aneH9gS3PQ/54QGFLhboD/zE2h473/cy2O793v/kO1xXh4/CP18zaB+JgFzgl6HTuZ+PW7U53F+0vfSeBL3FI60oHWdcP3RZgNfAlc18jyK2nH2jsEi3MCwj3B9ALfjJpvu5e3T0cu/a5TqqLwxdQT8Ctci9ntc3+MtuGD15kCZcVMa1eA9NCDM9J/yyp8dUt+34CbI3kXI5O7Nqfze+490HccHXcfhzhEbi+ss6uepFi1aIrPEIZH2OdDZ67CPdSOVE4ybGxDciMpy3BNT/CjEPQnl/+E66P8vsN4YM9UYExjteQ1wvfXfl+gk3Jfpk7i5BpcZY9YYY24zxmRZd2vsR7hb9uHeXop2/SwBuhtjzvTSP+Sln+Jtfw0XDIz1mT642/ILrbVlxph4457LXWyt/ZN1o0mfAq7xRpL7Fc3jnIu73fw0bkDMGbh+l+cCq40xn+Fagj63ITMFhOFIdfQkjaujvriR4Q/g6uE2XCvfpcBcY8zfvM/3ubV2s4/0B+AeU7ndGJPo3UrdaK19wFrbHfdM6OuNMXHeQJPmVn448nVczVfXcbj9bWNxncXiPBWRSGjqaPZoW4B2uIDpC9x8g//VaoT7Q3ydj7RzcH+kr/VeJ+BuBd2C66dUDSzFtYz9xGf5u+OmPrkc1wKZi5vn8HFcv8Aa3B/4EuCnzal+vPe2wn3Rb8O1jqTUss9iGvf4uJ94dXF80LokvNYZ3C3A9cA3faYfteOMa406B7ghZH0SrpXqf3ABfQ1wRXOsI68+LgZ+HrI+Czew5EZcd4MafLZ+4rpCLObrLauJfPWUmuG4Jx791xydzaT80b6Oo3qdxeo81aJFS2QWDcCJAmNMDu721Qm4KUQ+xT25YjPwY9wXSZ619ssw022N++98h7V2Xsi2VNwk4D/DDW5Jsz4Gfhhj0nGjrIts0EAVr8UhBzfx+g24L1O/eUSlfkLKeg+udfMA7vbhW7gJj6/wyt7X+h992wM3SjsD+I219tmQ7QNxfewy/OQRi+MclN5/DcIwxozCzUXYprnWUUhaifa/R0JfgJtaytdnMMYMwd3m3g7caa39a8j243HBUlYzLX8sruOoXmcheUXlPBWRyFAwGSXGTSj+P7jHA/bE3dbKBD7EPWf5lQjk8V8jfI0xf8L1JzotEunj+lKFTsfxFyDXWnt6I9KOSv14t1OrvWlLhuMer3Yy7hFv8bg5+6ZZa9/xW3YvnzTcwJJLcC1W//TSHujlu9haO64xeQTlFbHj7B1TE/rFHLLPnbgnoYzyUdzgdKJSR3WNQvZGLFdba60x5ndAoZ9zNFDfXleMB3Hnz25cv713cK2qY4D11sfcidEufy3pRvw6jvZ1FsvzVEQaT8FkBBljugC9vZdfAitw/7H3xD22rhzYZa0t8Zn+157sUcv2FFzLwJPW2v/1k0cd6RrcF0Q1bl692cB91trXw0wnqvVTR55JuFGtB3C35vb6bfH00jO40anVXj/PE3DTo5yB+yJdD7wIvGG9eUV95NEkxzko/dOBnTb853wH3h/VOgo6jwzuNufK4HS8/M8HtthGzO/ppdUK9/Sqb+NGVg/A3Rqehns6U9jTfMWy/CH5RuQ6riPtiF5nDczzdBpxnopI5CiYjBBjzHjcnHT5uKBoHe627fvADBuDDuLGmERca8Zcn++Pw32JdQBScSNMP7RBT5UwbsLxb1lr3w4z7ajWjwmZz/FIAVkkmaC5DY0xba21e6Ocn6/jHFpHsRSpOgo5j77EzfG5GXeL9U1r7cpGljP4GkjB9Yv8yFq71wssLe628K7mWP5aPkOkr+OoX2dNeZ6KiD8KJiPAu2W7BngIN3KyA64143TcLbGtuGdkr6jtlmUD0k/EzXu3wUbpaQ/eLck/At/EtZZsxn1xHsTdev6ztfYLn2lHu34ycYN23sa1eP07kEbwl51xk2hvtj6e5nKkY+Cn3OHm0ci066uj4Mmt+wHbrI9JxaNdR0c4j/rhztmfeudRvA1zhHIt18AWXOthOe627YvW2tXevmFP+B3t8tfxGSJ5HcfiOov6eSoiUWCbwSiglr4AE4B5dWwbjutrtQ5o7zP9SbhWjOdwE/RmEzIvHO7Rfefif4LjX+D+iJ/kvT4eN03Jk8B/gL8CHZpp/UwAKnADeapxrUl34Tr/B/bpipuzrmeUj8HZQGJzO85HQx3F4Dyq7xpYgBuQ4+saiEX5G/AZInEdR/scinoeWrRoifzS5AU4GhbcE21WAAO918nBX/a4qSxWABf7TH8u7nbwR94f2PW4p2YMxz27Gtz0HJ804jN8BEyuZX08Xz0y7d1mWj/TcPPpdcQ98/s+r7zVuNuH1+ImPt7fiPqJxTGIWh5HQx3F4DyK2jUQi/JH+zPE6ByKeh5atGiJ/KJJyyNjBu6W0iRjTJq1tsK6SXzjAKy1G4E9uGfXhsUY0wGoxI2MPA03f9wfcSOhZwOzjDG34FqF5tWVzhHySMBNmDzayw/jJpqOs9ZWW2tn44KALsaYfB9ZRLN+knFfwJustTustZ9Za3+Om/R7lLftTtwUJg/4KHusjkHU8jha6ojonkfRvgaiWv5of4YYnUNRz0NEoqSpo9mWvuD6VBnge7gnPpThvkSH8NVjBy/11uf5SL8z8FNgVC3bCnCTEO/GfUnlNuJznIy7pfQAtTznFndraX+4eUS7frw8kvEee0ctzwLG9UlrzGPpon4Mop1HS6+jGJ1HUbkGYlX+GHyGqJ5DscpDixYtkV80ACdCjDEZuC+EU3CTBZ/qbdqO+xL5s7X2Tp9pp+A6tx/0pvcAb4W3/R7gHGttgc/043BfaFfgno2dALwOvIp7WsaJuBam/tbak3zmkUEU6ifQKd8Y0xP40lpbXMu2XwGXW2t7+im7l1ZUj0E08zjK6iiD6JxHUb8Goln+aH+GWJxDsTpPRSTyFEw2gjGmI3AZ7vnJu3BzrO0B5uD69yTi5pN711q7qpF51ToS1rgnoiwEnrPWNvrWj/dldznuKTSDcC0lB3Gd9++zIU9kOUJaUa2foPQnAzuAKtzj3V7DzWH4pRfUXANstdbODDePkPyifgwincfRUEexvM68/DKI0DXQFOX38swg8tdx1M6hWJ+nIhJZCiYbwbinkAzAjfIswT1b9wTgONwfxNvD/eIJST8dKKvtyzlon1bARcDL1tpDkcjDa+FohZtIfCCulSDszxGD+qkt/QLcCNbNwG+ttf/0m76XR5Mcg0jlcTTUUVNcZ5G6BmJR/mh/hhidQ1HPQ0SiR8GkT95/yWW423azg9Z1A4YCV+Oe7HKhtXahzzyexk2R8Slu7r59teyTYRsx11oD88i01pbW1epUR7pRrZ960u+C6zd2DW4QyFi/9e+l2VyOQdh5HA111Iyus7CvgViVP5qfIRbnUKzOUxGJItsMOm62xAX3X/RS4OQ6ticD83G3lPykPxbX0XwPbu65p3F9rHoBKd4+gcfqDYxgHhfgbrkF8mgDvAmc0Mzq50jpJzUm/SY+BhHJ42iooya6ziJyDcSi/NH+DDE6h6KehxYtWqK7qGXSJ2+wwUzc48rGAWttyBMxjDETgKustYN8pD8NN7fag7gvhh/hvqBXAn8H/g/oC0yx1ib5/AxRyyMG9RPV9L336xg0Yflj8Rlaevmj/RliVP6o5yEi0aV5Jn2y7tmxv8C1urwAjDPGdDXGtIHDAw5G4uZ9C4tx88WtB/ZYa9dZa39nrT0BOAn3SLQfAdOBqcCf/ZQ/2nlEs35ikb6OQdOXP9qfoaWXPxafIdrlj1UeIhJlTd002tIXXMf2V3EjNHfhOpA/i5uKYx4+bo156WYCx3u/J+H1bw3afhHu1tagRpQ9FnlEpX5ikb6OQfMofzQ/Q0sv/9FwDsUyDy1atERn0W3uCDFuaotzcZMSH8T9F/2atfaLCOYRh/uiqDbGXIO7bZUaqfSjmUe06ycW9e/lo2Nw5HyiVkct/Tpr6ccgRvUfkzoSkchRMBkFxj2+rObIezYqj8lAvLX2ty0tj2jXTyzq38tHx+DI+UStjlr6ddbSj0GM6j8mdSQijaNgsoUyxiQC1VEOyqKeR0umY3BkKn/TOxo+g4g0bwomRURERMQ3jeYWEREREd8UTIqIiIiIbwomRURERMQ3BZMiIiIi4puCSRERERHxTcGkiIiIiPj2/wEDPsJe4okfRAAAAABJRU5ErkJggg==\n", + "text/plain": [ + "
" + ] + }, + "execution_count": 10, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "legend = ['Mitigated Probabilities', 'Unmitigated Probabilities']\n", + "plot_histogram([mitigated_probs, unmitigated_probs], legend=legend, sort=\"value_desc\", bar_labels=False)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Expectation value" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [], + "source": [ + "diagonal_labels = [\"ZZZZ\", \"ZIZI\", \"IZII\", \"1ZZ0\"]\n", + "ideal_expectation = []\n", + "diagonals = [str2diag(d) for d in diagonal_labels]\n", + "qubit_index = {i: i for i in range(num_qubits)}\n", + "unmitigated_probs_vector, _ = counts_probability_vector(unmitigated_probs, qubit_index=qubit_index)\n", + "unmitigated_expectation = [expval_with_stddev(d, unmitigated_probs_vector, SHOTS) for d in diagonals]\n", + "mitigated_expectation = [mitigator.expectation_value(counts, d) for d in diagonals]" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 12, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD4CAYAAAD8Zh1EAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAb30lEQVR4nO3de3iU1dnv8e9NoI0CYovRrQYMUs4JSSCiFoWkQEWrUEFtULsBC1RFbbUo2Coo+u6rfcVDPbHFqnjg5GFL04r1UKFo3RawQIuglkPEaKuIEIwETer9/jHJdAiZZBImhCx+n+ua68qsZ83z3LMy/HjyHNaYuyMiIi1fq+YuQEREkkOBLiISCAW6iEggFOgiIoFQoIuIBKJ1c234qKOO8oyMjObavIhIi/Tmm29+4u5ptS1rtkDPyMhg1apVzbV5EZEWyczei7dMh1xERAKhQBcRCYQCXUQkEM12DF2kJaioqKCkpIQ9e/Y0dylyiElNTSU9PZ02bdok/BoFukgdSkpKaN++PRkZGZhZc5cjhwh3Z/v27ZSUlNClS5eEX6dDLiJ12LNnDx07dlSYywFlZnTs2LHBfxkq0EXqoTCX5tCYz50CXUQkEEEHen5+Pvn5+c1dhgTELLmPxLZpXHzxxdHnlZWVpKWlcfbZZwNQVFTEL3/5SwAWL17M+vXro32nT5/Oyy+/3Kj3umbNGpYsWdLg1+Xn59d602B+fj49evQgJyeHnJwczjvvvEbV1VB33XUXu3fvbnC/s846i507dzZhZckXdKCLhKBt27asW7eO8vJyAF566SWOP/746PIRI0Ywbdo0YN9AnzlzJkOHDm3Udhsb6HWZN28ea9asYc2aNTz99NNJXXc8jQ30JUuWcOSRRzZhZcmnQBdpAc466yyee+45ABYsWMCYMWOiy+bOncsVV1zB66+/TlFREddeey05OTls2rSJcePGRYNzyZIl9OzZk/79+3PVVVdF9/BXrFjBqaeeSm5uLt/+9rd55513+PLLL5k+fTqLFi0iJyeHRYsW8fnnn3PJJZcwYMAAcnNz+e1vfwtAeXk5hYWF9OrVi3PPPTf6H0+iRo4cyWOPPQbAAw88wEUXXQRE9uh/8pOfkJOTQ2ZmJitWrACIW8e///1vpkyZQmZmJn379uWee+7h7rvv5sMPP6SgoICCggIALrvsMvLy8ujTpw8zZswAqLVfRkYGn3zyCQB33HEHmZmZZGZmctdddwFQXFxMr169mDhxIn369OG73/1ug9970rl7szz69+/vTW3w4ME+ePDgJt+OhGv9+vV7PYfkPhLRtm1bX7t2rY8ePdrLy8s9Ozvbly5d6t/73vfc3f2RRx7xyZMnu7v72LFj/amnnoq+tvp5eXm5p6en++bNm93dvbCwMPr60tJSr6iocHf3l156yUeNGrXPet3dr7/+en/88cfd3X3Hjh3erVs3Lysr89tvv93Hjx/v7u5r1671lJQUX7ly5T7vY/Dgwd69e3fPzs727OxsnzJliru7/+tf//KuXbv68uXLvVu3br59+/Zo/wkTJri7+5/+9Cfv06dPnXXcf//9Pnr06Oh7qV7PCSec4Nu2bYvWUd1eWVnpgwcP9rVr19bar/r5qlWrPDMz08vKyvyzzz7z3r17+1//+lffsmWLp6Sk+OrVq93d/fzzz4/WlSw1P3/u7sAqj5Orug5dpAXo27cvxcXFLFiwgLPOOqvBr3/77bc58cQTo9c0jxkzhjlz5gBQWlrK2LFj+cc//oGZUVFRUes6XnzxRYqKipg1axYQuaRz69atLF++nKuuuipaZ9++fePWMW/ePPLy8vZqO+aYY5g5cyYFBQU8++yzfPOb34wuq/5LZNCgQezatYudO3fGrePll1/m0ksvpXXrSKzFrifWk08+yZw5c6isrOSf//wn69evr7Pm1157jXPPPZe2bdsCMGrUKF599VVGjBhBly5dyMnJAaB///4UFxfHXc+BoEAXaSFGjBjBlClTWLZsGdu3b0/aem+88cZomBYXF8e9kMDdeeaZZ+jRo0fStl3t73//Ox07duTDDz/cq73mpXtmtl91bNmyhVmzZrFy5Uq+8Y1vMG7cuP26C/jrX/969OeUlJRmP+SiY+jSKLqC6MC75JJLmDFjBllZWXH7tG/fns8++2yf9h49erB58+boHuSiRYuiy0pLS6MnWefOnRt3XWeccQb33HMPkb/6YfXq1UBk73n+/PkArFu3jr/97W8Nel8rVqzg+eefZ/Xq1cyaNYstW7ZEl1XX+dprr9GhQwc6dOgQt45hw4bxwAMPUFlZCcCnn366z/vYtWsXbdu2pUOHDnz00Uc8//zz9Y7d6aefzuLFi9m9ezeff/45zz77LKeffnqD3uOBokAXaYBkH0VviPT09OihjXgKCwu57bbbyM3NZdOmTdH2ww47jPvvv5/hw4fTv39/2rdvT4cOHQC47rrruP7668nNzY2GIUBBQQHr16+PnhS98cYbqaiooG/fvvTp04cbb7wRiJxkLCsro1evXkyfPp3+/fvHre+iiy6KXrY4dOhQvvjiCyZOnMjDDz/Mcccdx+23384ll1wSDevU1FRyc3O59NJLeeihhwDi1jFhwgQ6d+5M3759yc7Ojv4nM2nSJIYPH05BQQHZ2dnk5ubSs2dPLrzwQgYOHBitLbZfrH79+jFu3DgGDBjAySefzIQJE8jNza37l9VMzBv6qUqSvLw8b+ovuKjeg1y2bFmTbudQdKiM7YYNG+jVq1dzl5EUZWVltGvXDndn8uTJdOvWjauvvrq5y4orPz+fWbNm7XPM/VBS2+fPzN5091oHpUUeQ2/oHbGJ9m+m/9tEDogHH3yQRx99lC+//JLc3Fx+/OMfN3dJkmQtMtBFpOGuvvrqg3qPvKbQ//prCjqGLiISCAW6iEggFOgiIoFQoIuIBEInRWUvuoKobnZzcr/swmfUPTDFxcWcffbZrFu3Ltp200030a5dO6ZMmbLf258+fTqDBg1i6NCh3HXXXUyaNInDDz8ciEwINn/+/EbNOLh48WK6d+9O7969G/S6du3aUVZWtk97SkrKXjdUFRYWRmeYbCo7d+5k/vz5XH755Q3q9+GHH3LVVVcdsNkkY2kPXeQQFju9bjKnj605je/+Ouyww6LT7q5Zs6bJwxwiQX3//fc3uN9xxx3XLGEOwQf6sqqHSJjy8/OZOnUqAwYMoHv37rz66qtA5Bb+73//+wwbNoyMjAzuvfde7rjjDnJzcznllFOit8VXT69b3/Sxt9xyCz169OC0005jzJgx0YmxHnzwQU466SSys7MZPXo0u3fvrnUa302bNkXvUj399NN5++23gcjcKqeeeipZWVnccMMNDXrvpaWl9OjRg3feeQeITOT14IMPApE9/auvvpo+ffowZMgQtm3bBhC3jo8++ohzzz2X7OxssrOzef3115k2bRqbNm0iJyeHa6+9lrKyMoYMGUK/fv3IysqKTttbs19xcTGZmZlAZOKw8ePHk5WVRW5uLkuXLo3+fkaNGsXw4cPp1q0b1113XUN/9bWLNw1jUz/2Z/rc5N+AnfhUpqFLfLwGVz3CHtt9ps+9iaQ+6rNly5botLHVZsyY4bfddpu7R6aYveaaa9zd/bnnnvMhQ4a4e2Tq265du/quXbv8448/9iOOOMJnz57t7u4//elP/c4773T3vafbjTd97IoVKzw7O9vLy8t9165d/q1vfSu6/U8++STa/xe/+IXffffd+6zX3f073/mOv/vuu+7u/sYbb3hBQYG7u59zzjn+6KOPurv7vffe623btq11HFq1ahWddjc7O9sXLlzo7u4vvviin3LKKb5gwQI/44wz/vN7An/iiSfc3f3mm2+OTgMcr44LLrggOiaVlZW+c+fOfca+oqLCS0tL3d1927Zt3rVrV//qq6/26Rf7fNasWdGphTds2OCdOnXy8vJyf+SRR7xLly6+c+dOLy8v986dO/vWrVv3ed+aPlckIPG+KDi2fdSoUcC+07cWFBTQvn376Lwt55xzDgBZWVkNmkDrz3/+MyNHjiQ1NZXU1NToeiAyGdcNN9zAzp07KSsr44wzztjn9WVlZbz++uucf/750bYvvvgiuu5nnnkGgB/+8IdMnTq11hqqD7nUNGzYMJ566ikmT57M2rVro+2tWrXiBz/4AQAXX3wxo0aNqrOOV155JfolGykpKXTo0IEdO3bstS135+c//znLly+nVatWfPDBB3z00UfxB47IpGJXXnklAD179uSEE07g3XffBWDIkCHR+XR69+7Ne++9R6dOnepcX30U6CIHsY4dO+4TLJ9++ml0XnP4zxSuKSkpe02uFTu1a6tWraLPW7VqtVe//TFu3DgWL15MdnY2c+fOrfXuzq+++oojjzyy1kCGxn27fey6N2zYwOGHH86OHTtIT0+Pu4366qjPvHnz2LZtG2+++SZt2rQhIyMjqVPvJuN3EvgxdJGWrV27dhx77LG88sorQCTM//CHP3DaaaclfVvxpo8dOHAgv/vd79izZw9lZWX8/ve/jy777LPPOPbYY6moqGDevHm1ruuII46gS5cuPPXUU0BkT7d6b3rgwIEsXLgQYK/XJ+rOO++kV69ezJ8/n/Hjx0e/nOOrr76KnpicP38+p512Wp11DBkyhNmzZwORr7IrLS3dZzxKS0s5+uijadOmDUuXLuW9996rc9wgMvVu9ft699132bp1a5PMJ19Ne+giDVDfZYZN4bHHHmPy5Mlcc801AMyYMYOuXbsmfTvV08ced9xx0ZN3ACeddBIjRoygb9++HHPMMWRlZUUPFdxyyy2cfPLJpKWlcfLJJ0eDrbCwkIkTJ3L33Xfz9NNPM2/ePC677DJuvfVWKioqKCwsJDs7m1//+tdceOGF/OpXv2LkyJFxaysvL49+MxDA8OHDGT9+PL/5zW9YsWIF7du3Z9CgQdx6663cfPPNtG3blhUrVnDrrbdy9NFHR+dVr6uOSZMm8dBDD5GSksLs2bM59dRTGThwIJmZmZx55plMnTqVc845h6ysLPLy8ujZsycQ+Ssqtt/kyZOjdV5++eVcdtllZGVl0bp1a+bOnbvXnnmyJTR9rpkNB34NpAC/cfdf1ljeGXgUOLKqzzR3r/Prwvdn+tz9+AutTqFcK70/NLZ7C2n63P1RPfXu7t27GTRoEHPmzKFfv37NXVZc8a5nb2mSPn2umaUA9wHDgBJgpZkVuXvsRaY3AE+6+2wz6w0sATIa9xZE5GAzadIk1q9fz549exg7duxBHeaHskQOuQwANrr7ZgAzWwiMBGID3YEjqn7uAOz9xYAiCTpUvjijpan+9p+WIoS988ZI5KTo8cD7Mc9Lqtpi3QRcbGYlRPbOr6xtRWY2ycxWmdmq6gv9RQ52iRyWFEm2xnzuknWVyxhgrrunA2cBj5vZPut29znunufueWlpaUnatEjTSU1NZfv27Qp1OaDcne3bt5Oamtqg1yVyyOUDIPZq9/Sqtlg/AoZXFfL/zSwVOAr4uEHViBxk0tPTKSkpQX9RyoGWmpoa97r6eBIJ9JVANzPrQiTIC4ELa/TZCgwB5ppZLyAV0L8AiUp4lsLihvVv6ssI27Rps9dNPCIHs3oPubh7JXAF8AKwgcjVLG+Z2UwzG1HV7WfARDNbCywAxrn+RhUROaASurGo6pryJTXapsf8vB4YmNzSRESkIXTrfyPk5+dHL68TETlY6NZ/ObiMb+4CRFou7aGLiARCgS4iEggFuohIIHQMPUZLvVZaRAS0hy4iEgwFuohIIBToIiKB0DH0xtC10iJyENIeuohIIBToIiKBUKCLiARCgS4iEggFuohIIBToIiKBUKCLiARCgS4iEggFuohIIBToIiKBUKCLiARCgS4iEggFuohIIBToIiKBUKCLiARCgS4iEggFuohIIBToIiKBUKCLiARCgS4iEggFuohIIBToIiKBSCjQzWy4mb1jZhvNbFqcPheY2Xoze8vM5ie3TBERqU/r+jqYWQpwHzAMKAFWmlmRu6+P6dMNuB4Y6O47zOzopipYRERql8ge+gBgo7tvdvcvgYXAyBp9JgL3ufsOAHf/OLlliohIfRIJ9OOB92Oel1S1xeoOdDezP5vZG2Y2PFkFiohIYuo95NKA9XQD8oF0YLmZZbn7zthOZjYJmATQuXPnJG1aREQgsT30D4BOMc/Tq9pilQBF7l7h7luAd4kE/F7cfY6757l7XlpaWmNrFhGRWiQS6CuBbmbWxcy+BhQCRTX6LCayd46ZHUXkEMzm5JUpIiL1qTfQ3b0SuAJ4AdgAPOnub5nZTDMbUdXtBWC7ma0HlgLXuvv2pipaRET2ldAxdHdfAiyp0TY95mcHrql6iIhIM9CdoiIigVCgi4gEQoEuIhIIBbqISCAU6CIigVCgi4gEQoEuIhIIBbqISCAU6CIigVCgi4gEQoEuIhIIBbqISCAU6CIigVCgi4gEQoEuIhIIBbqISCAU6CIigVCgi4gEQoEuIhIIBbqISCAU6CIigVCgi4gEQoEuIhIIBbqISCAU6CIigVCgi4gEQoEuIhIIBbqISCAU6CIigVCgi4gEQoEuIhIIBbqISCAU6CIigUgo0M1suJm9Y2YbzWxaHf1Gm5mbWV7yShQRkUTUG+hmlgLcB5wJ9AbGmFnvWvq1B34C/CXZRYqISP0S2UMfAGx0983u/iWwEBhZS79bgF8Be5JYn4iIJCiRQD8eeD/meUlVW5SZ9QM6uftzda3IzCaZ2SozW7Vt27YGFysiIvHt90lRM2sF3AH8rL6+7j7H3fPcPS8tLW1/Ny0iIjESCfQPgE4xz9Or2qq1BzKBZWZWDJwCFOnEqIjIgZVIoK8EuplZFzP7GlAIFFUvdPdSdz/K3TPcPQN4Axjh7quapGIREalVvYHu7pXAFcALwAbgSXd/y8xmmtmIpi5QREQS0zqRTu6+BFhSo216nL75+1+WiIg0lO4UFREJhAJdRCQQCnQRkUAo0EVEAqFAFxEJhAJdRCQQCnQRkUAo0EVEAqFAFxEJhAJdRCQQCnQRkUAo0EVEAqFAFxEJhAJdRCQQCnQRkUAo0EVEAqFAFxEJhAJdRCQQCnQRkUAo0EVEAqFAFxEJhAJdRCQQCnQRkUAo0EVEAqFAFxEJhAJdRCQQCnQRkUAo0EVEAqFAFxEJhAJdRCQQCnQRkUAo0EVEApFQoJvZcDN7x8w2mtm0WpZfY2brzexvZvZHMzsh+aWKiEhd6g10M0sB7gPOBHoDY8ysd41uq4E8d+8LPA38d7ILFRGRuiWyhz4A2Ojum939S2AhMDK2g7svdffdVU/fANKTW6aIiNQnkUA/Hng/5nlJVVs8PwKer22BmU0ys1Vmtmrbtm2JVykiIvVK6klRM7sYyANuq225u89x9zx3z0tLS0vmpkVEDnmtE+jzAdAp5nl6VdtezGwo8AtgsLt/kZzyREQkUYnsoa8EuplZFzP7GlAIFMV2MLNc4AFghLt/nPwyRUSkPvUGurtXAlcALwAbgCfd/S0zm2lmI6q63Qa0A54yszVmVhRndSIi0kQSOeSCuy8BltRomx7z89Ak1yUiIg2kO0VFRAKhQBcRCYQCXUQkEAp0EZFAKNBFRAKhQBcRCYQCXUQkEAp0EZFAKNBFRAKhQBcRCYQCXUQkEAp0EZFAKNBFRAKhQBcRCYQCXUQkEAp0kUNEfn4++fn5zV2GNCEFuohIIBToIiKBUKCLiARCgS4iEggFushBRicvpbEU6CIigWjd3AWIyP6xmy2xjsUN6+8zvHEFSbPRHrqISCAU6CIigVCgi8ghI/QTzgp0EZFAKNBFRAKhq1xEDhXjm7sAaWoKdJEDxBK8urDB/W9qaCUSKh1yERHZTwfLyVbtoYuIxNHSbtpKaA/dzIab2TtmttHMptWy/Otmtqhq+V/MLCPplYqIHKzGc1Cco6g30M0sBbgPOBPoDYwxs941uv0I2OHu3wLuBH6V7EJFRKRuiRxyGQBsdPfNAGa2EBgJrI/pM5L/nJp5GrjXzMzdNRmEiDQ5nXCOSCTQjwfej3leApwcr4+7V5pZKdAR+CS2k5lNAiYBdO7cuZElQ9P9N6H/fzS2TSfRsa0+t7ZsWcJrbngxgdHYRhzQq1zcfY6757l7Xlpa2oHctIhI8BLZQ/8A6BTzPL2qrbY+JWbWGugAbE9KhSKHmGWJ7z6K7CWRPfSVQDcz62JmXwMKgaIafYqAsVU/nwe8ouPnIiIHVr176FXHxK8AXgBSgIfd/S0zmwmscvci4CHgcTPbCHxKJPRFROQASujGIndfAiyp0TY95uc9wPnJLU1ERBpCt/6LiARCt/6LyCEj9BPO2kMXEQmEAl1EJBAKdBGRQCjQRUQCoUAXEQmEAl1EJBAKdBGRQCjQRUQCoUAXEQmENdekiGa2DXivWTYe31HU+FIOSRqNbdPR2Dadg3FsT3D3Wr9QotkC/WBkZqvcPa+56wiRxrbpaGybTksbWx1yEREJhAJdRCQQCvS9zWnuAgKmsW06Gtum06LGVsfQRUQCoT10EZFAKNBFRAIR7DcWmdm5wIwazX2BK4GJMW2tgT5Ab+Bi4Hsxyw4HugIdgHuA3JhlRwKHufsxSS28BapjrCcDk90908z+i/hjex6Q5+5XHIh6WyIzK3P3dmY2mfif32OAKe5+tpmNQ2O6DzN7GDgb+Ljqc5kCvFmjWzrwR2AK8Lsay04EZrv7VDPrD8wFDiPyncs/8eY+hu3uh8QDmAT8CWhVo/3/AE/Eec084NZa2lsBy4EJzf2+DsZHzFifCKyrb2yBccC9zV33wfwAyuK0Rz+/QD7we41pneM4COhXx+fyWOB9ILOWZVlVy/5X1fMVwCmAAc8DZzb3+wt2Dz2WmXUHpgPfdvevYtoHARcQ+QXXfM3FwLeAsbWs8ufANnf/TdNU3HLFjjVxDunVM7aSoLo+v1I7d19uZhm1LTMzAx4FbnP3dTWWpQLzifzF+S8zOxY4wt3fqFr+GPB9IsHebIIPdDNrQ+QX8TN33xrTfiSRP5d+6O67arwmA/glkO/ulTWWDQAmoH9E+6g51rX9w6lrbCVxdX1+pdGuBiqJHF6t6b+B19y9qOr58UBJzPKSqrZmFXygA7cAb7n7ohrt/xd43N3/HNtYdUztCeBGd99YY1m7qmU/cvdPm7DmlireWAN1j600WK2fX2kcM8sGfgqc5FXHU2KWnQkMBfo3Q2kNEnSgm1k+MJoae9NmNhY4gchJ0JpuAP7p7o/Usuwe4Lfu/sfkVtryxRvrGuoaW0lQPZ9faSAzO4zIOZ3L3P2jGsuOBh4ARrp7ecyiD4icPK2WXtXWrIINdDP7BvAIcKG7fxbTfiKRE0mn13I45RQiJ5NqO6Z+HpBN5CSIxIg31jX6xB1bSVxdn19ptFnAn9z9uVqWPQzc4+6rYxvd/Z9mtqvqc/0X4H9T+6GaAyrYQAcuBY4GZkfOdUR1IHLJ3P+r0X4lkT3Iw4GlNZaNBv6ratmKGstOrfE/96Eo3lgviPn5ZuKPrSRuKvE/v1IPM1tA5Gqgo8ysBJgNXA68bWZrYrq+BdxL5FLbTmZ2Ucyyl9z92qrXzSVy2eLzNPMJUdCt/yIiwdCdoiIigVCgi4gEQoEuIhIIBbqISCAU6CIigVCgi4gEQoEuIhKI/wG75swnAbDp3QAAAABJRU5ErkJggg==\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "mitigated_expectation_values, mitigated_stddev = zip(*mitigated_expectation)\n", + "unmitigated_expectation_values, unmitigated_stddev = zip(*unmitigated_expectation)\n", + "legend = ['Mitigated Expectation', 'Unmitigated Expectation']\n", + "fig, ax = plt.subplots()\n", + "X = np.arange(4)\n", + "ax.bar(X + 0.00, mitigated_expectation_values, yerr=mitigated_stddev, color='b', width = 0.25, label=\"Mitigated Expectation\")\n", + "ax.bar(X + 0.25, unmitigated_expectation_values, yerr=unmitigated_stddev, color='g', width = 0.25, label=\"Unmitigated Expectation\")\n", + "ax.set_xticks([0.125 + i for i in range(len(diagonals))])\n", + "ax.set_xticklabels(diagonal_labels)\n", + "ax.legend()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Correlated readout mitigation\n", + "\n", + "In correlated readout mitigation on $n$ qubits, a circuit is generated for each of the possible $2^n$ combinations of \"0\" and \"1\". This results in more accurate mitigation in the case where the readout errors are correlated and not independent, but requires a large amount of circuits and storage space, and so is infeasible for more than a few qubits." + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + " ░ ┌─┐ \n", + " q_0: ─░─┤M├───\n", + " ░ └╥┘┌─┐\n", + " q_1: ─░──╫─┤M├\n", + " ░ ║ └╥┘\n", + "meas: 2/════╩══╩═\n", + " 0 1 \n", + " ┌───┐ ░ ┌─┐ \n", + " q_0: ┤ X ├─░─┤M├───\n", + " └───┘ ░ └╥┘┌─┐\n", + " q_1: ──────░──╫─┤M├\n", + " ░ ║ └╥┘\n", + "meas: 2/═════════╩══╩═\n", + " 0 1 \n", + " ░ ┌─┐ \n", + " q_0: ──────░─┤M├───\n", + " ┌───┐ ░ └╥┘┌─┐\n", + " q_1: ┤ X ├─░──╫─┤M├\n", + " └───┘ ░ ║ └╥┘\n", + "meas: 2/═════════╩══╩═\n", + " 0 1 \n", + " ┌───┐ ░ ┌─┐ \n", + " q_0: ┤ X ├─░─┤M├───\n", + " ├───┤ ░ └╥┘┌─┐\n", + " q_1: ┤ X ├─░──╫─┤M├\n", + " └───┘ ░ ║ └╥┘\n", + "meas: 2/═════════╩══╩═\n", + " 0 1 \n" + ] + } + ], + "source": [ + "qubits = [0,3]\n", + "num_qubits = len(qubits)\n", + "exp = ReadoutMitigationExperiment(qubits, method=\"correlated\")\n", + "for c in exp.circuits():\n", + " print(c)" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "

This code is a part of Qiskit

© Copyright IBM 2017, 2022.

This code is licensed under the Apache License, Version 2.0. You may
obtain a copy of this license in the LICENSE.txt file in the root directory
of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.

Any modifications or derivative works of this code must retain this
copyright notice, and modified files need to carry a notice indicating
that they have been altered from the originals.

" + ], + "text/plain": [ + "" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "import qiskit.tools.jupyter\n", + "%qiskit_copyright" + ] + } + ], + "metadata": { + "interpreter": { + "hash": "42e50baa22dbe56fc7f4f6bf32ac20952839a6a23dcf2ef84eded7e8cac03444" + }, + "kernelspec": { + "display_name": "qiskit38", + "language": "python", + "name": "qiskit38" + }, + "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.12" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/qiskit_experiments/library/__init__.py b/qiskit_experiments/library/__init__.py index 9dd3e5ac17..08dd8b22dc 100644 --- a/qiskit_experiments/library/__init__.py +++ b/qiskit_experiments/library/__init__.py @@ -142,6 +142,7 @@ class instance to manage parameters and pulse schedules. from .randomized_benchmarking import StandardRB, InterleavedRB from .tomography import StateTomography, ProcessTomography from .quantum_volume import QuantumVolume +from .mitigation import ReadoutMitigationExperiment # Experiment Sub-modules from . import calibration diff --git a/qiskit_experiments/library/mitigation/__init__.py b/qiskit_experiments/library/mitigation/__init__.py new file mode 100644 index 0000000000..25d6a23a35 --- /dev/null +++ b/qiskit_experiments/library/mitigation/__init__.py @@ -0,0 +1,41 @@ +# This code is part of Qiskit. +# +# (C) Copyright IBM 2021. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. + +""" +=============================================================================================== +Readout Mitigation Experiments (:mod:`qiskit_experiments.library.mitigation`) +=============================================================================================== + +.. currentmodule:: qiskit_experiments.library.mitigation + +Experiment +=========== +.. autosummary:: + :toctree: ../stubs/ + :template: autosummary/experiment.rst + + ReadoutMitigationExperiment + + +Analysis +======== + +.. autosummary:: + :toctree: ../stubs/ + :template: autosummary/analysis.rst + + CorrelatedMitigationAnalysis + LocalMitigationAnalysis +""" +from .mitigation_experiment import ReadoutMitigationExperiment +from .mitigation_analysis import CorrelatedMitigationAnalysis +from .mitigation_analysis import LocalMitigationAnalysis diff --git a/qiskit_experiments/library/mitigation/mitigation_analysis.py b/qiskit_experiments/library/mitigation/mitigation_analysis.py new file mode 100644 index 0000000000..8ec1049044 --- /dev/null +++ b/qiskit_experiments/library/mitigation/mitigation_analysis.py @@ -0,0 +1,161 @@ +# This code is part of Qiskit. +# +# (C) Copyright IBM 2021. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. +""" +Measurement calibration analysis classes +""" +from typing import List, Tuple +import numpy as np +import matplotlib.pyplot as plt +from qiskit.result import CorrelatedReadoutMitigator +from qiskit.result import LocalReadoutMitigator +from qiskit.result import marginal_counts +from qiskit_experiments.framework import ExperimentData +from qiskit_experiments.framework.matplotlib import get_non_gui_ax +from qiskit_experiments.framework import BaseAnalysis, AnalysisResultData, Options + + +class CorrelatedMitigationAnalysis(BaseAnalysis): + """ + Measurement correction analysis for a full calibration + """ + + def _run_analysis( + self, experiment_data: ExperimentData, **options + ) -> Tuple[List[AnalysisResultData], List["matplotlib.figure.Figure"]]: + data = experiment_data.data() + qubits = experiment_data.metadata["physical_qubits"] + labels = [datum["metadata"]["label"] for datum in data] + matrix = self._generate_matrix(data, labels) + result_mitigator = CorrelatedReadoutMitigator(matrix, qubits=qubits) + analysis_results = [AnalysisResultData("Correlated Readout Mitigator", result_mitigator)] + ax = options.get("ax", None) + figures = [self._plot_calibration(matrix, labels, ax)] + return analysis_results, figures + + def _generate_matrix(self, data, labels) -> np.array: + list_size = len(labels) + matrix = np.zeros([list_size, list_size], dtype=float) + # matrix[i][j] is the probability of counting i for expected j + for datum in data: + expected_outcome = datum["metadata"]["label"] + j = labels.index(expected_outcome) + total_counts = sum(datum["counts"].values()) + for measured_outcome, count in datum["counts"].items(): + i = labels.index(measured_outcome) + matrix[i][j] = count / total_counts + return matrix + + def _plot_calibration(self, matrix, labels, ax=None) -> "matplotlib.figure.Figure": + """ + Plot the calibration matrix (2D color grid plot). + + Args: + matrix: calibration matrix to plot + ax (matplotlib.axes): settings for the graph + + Returns: + The generated plot of the calibration matrix + + Raises: + QiskitError: if _cal_matrices was not set. + + ImportError: if matplotlib was not installed. + + """ + + if ax is None: + ax = get_non_gui_ax() + figure = ax.get_figure() + ax.matshow(matrix, cmap=plt.cm.binary, clim=[0, 1]) + ax.set_xlabel("Prepared State") + ax.xaxis.set_label_position("top") + ax.set_ylabel("Measured State") + ax.set_xticks(np.arange(len(labels))) + ax.set_yticks(np.arange(len(labels))) + ax.set_xticklabels(labels) + ax.set_yticklabels(labels) + return figure + + +class LocalMitigationAnalysis(BaseAnalysis): + """ + Measurement correction analysis for a full calibration + """ + + @classmethod + def _default_options(cls) -> Options: + """Return default analysis options. + + Analysis Options: + plot (bool): Set ``True`` to create figure for fit result. + ax(AxesSubplot): Optional. A matplotlib axis object to draw. + """ + options = super()._default_options() + options.plot = True + options.ax = None + return options + + def _run_analysis( + self, experiment_data: ExperimentData + ) -> Tuple[List[AnalysisResultData], List["matplotlib.figure.Figure"]]: + data = experiment_data.data() + qubits = experiment_data.metadata["physical_qubits"] + matrices = self._generate_matrices(data) + result_mitigator = LocalReadoutMitigator(matrices, qubits=qubits) + analysis_results = [AnalysisResultData("Local Readout Mitigator", result_mitigator)] + if self.options.plot: + figure = assignment_matrix_visualization( + result_mitigator.assignment_matrix(), ax=self.options.ax + ) + figures = [figure] + else: + figures = None + return analysis_results, figures + + def _generate_matrices(self, data) -> List[np.array]: + num_qubits = len(data[0]["metadata"]["label"]) + counts = [None, None] + for result in data: + for i in range(2): + if result["metadata"]["label"] == str(i) * num_qubits: + counts[i] = result["counts"] + matrices = [] + for k in range(num_qubits): + matrix = np.zeros([2, 2], dtype=float) + marginalized_counts = [] + for i in range(2): + marginalized_counts.append(marginal_counts(counts[i], [k])) + # matrix[i][j] is the probability of counting i for expected j + for i in range(2): + for j in range(2): + matrix[i][j] = marginalized_counts[j][str(i)] / sum( + marginalized_counts[j].values() + ) + matrices.append(matrix) + return matrices + + +def assignment_matrix_visualization(assignment_matrix, ax=None): + """Displays a visualization of the assignment matrix compared to the identity""" + if ax is None: + ax = get_non_gui_ax() + figure = ax.get_figure() + n = len(assignment_matrix) + diff = np.abs(assignment_matrix - np.eye(n)) + im2 = ax.matshow(diff, cmap=plt.cm.Reds, vmin=0, vmax=0.2) + ax.set_yticks(np.arange(n)) + ax.set_xticks(np.arange(n)) + ax.set_yticklabels(n * [""]) + ax.set_xticklabels(n * [""]) + ax.set_xlabel(r"$|A - I|$", fontsize=16) + figure.colorbar(im2, ax=ax) + return figure diff --git a/qiskit_experiments/library/mitigation/mitigation_experiment.py b/qiskit_experiments/library/mitigation/mitigation_experiment.py new file mode 100644 index 0000000000..a5febc3372 --- /dev/null +++ b/qiskit_experiments/library/mitigation/mitigation_experiment.py @@ -0,0 +1,114 @@ +# This code is part of Qiskit. +# +# (C) Copyright IBM 2021. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. +""" +Measurement calibration experiment classes. +""" +from typing import Iterable, List + +from qiskit import QuantumCircuit +from qiskit.exceptions import QiskitError +from qiskit_experiments.framework import BaseExperiment +from .mitigation_analysis import CorrelatedMitigationAnalysis, LocalMitigationAnalysis + + +class ReadoutMitigationExperiment(BaseExperiment): + """Class for readout mitigation experiments""" + + METHOD_LOCAL = "local" + METHOD_CORRELATED = "correlated" + ALL_METHODS = [METHOD_LOCAL, METHOD_CORRELATED] + + def __init__(self, qubits: Iterable[int], method=METHOD_LOCAL): + """Initialize a mitigation experiment. + + Args: + qubits: The qubits being mitigated + method: A string denoting mitigation method + + Raises: + QiskitError: if the given mitigation method is not recoginzed + + Additional info: + The currently supported mitigation methods are: + * "local": each qubit is mitigated by itself; this is the default method, + and assumes readout errors are independent for each qubits + * "correlated": All the qubits are mitigated together; this results in an exponentially + large mitigation matrix and so is useable only for a small number of qubits, + but might be more accurate than local mitigation. + """ + super().__init__(qubits) + if method not in self.ALL_METHODS: + raise QiskitError("Method {} not recognized".format(method)) + if method == self.METHOD_LOCAL: + self.helper = LocalMitigationHelper(self.num_qubits) + if method == self.METHOD_CORRELATED: + self.helper = CorrelatedMitigationHelper(self.num_qubits) + + self.analysis = self.helper.analysis() + + def circuits(self) -> List[QuantumCircuit]: + """Returns the experiment's circuits""" + return [self._calibration_circuit(self.num_qubits, label) for label in self.helper.labels()] + + @staticmethod + def _calibration_circuit(num_qubits: int, label: str) -> QuantumCircuit: + """Return a calibration circuit. + + This is an N-qubit circuit where N is the length of the label. + The circuit consists of X-gates on qubits with label bits equal to 1, + and measurements of all qubits. + """ + circ = QuantumCircuit(num_qubits, name="meas_mit_cal_" + label) + for i, val in enumerate(reversed(label)): + if val == "1": + circ.x(i) + circ.measure_all() + circ.metadata = {"label": label} + return circ + + +class CorrelatedMitigationHelper: + """Helper class for correlated mitigation experiment data""" + + def __init__(self, num_qubits: int): + """Creates the helper class + Args: + num_qubits: The number of qubits being mitigated + """ + self.num_qubits = num_qubits + + def analysis(self): + """Returns the analysis class for the mitigation""" + return CorrelatedMitigationAnalysis() + + def labels(self) -> List[str]: + """Returns the labels dictating the generation of the mitigation circuits""" + return [bin(j)[2:].zfill(self.num_qubits) for j in range(2 ** self.num_qubits)] + + +class LocalMitigationHelper: + """Helper class for local mitigation experiment data""" + + def __init__(self, num_qubits: int): + """Creates the helper class + Args: + num_qubits: The number of qubits being mitigated + """ + self.num_qubits = num_qubits + + def analysis(self): + """Returns the analysis class for the mitigation""" + return LocalMitigationAnalysis() + + def labels(self) -> List[str]: + """Returns the labels dictating the generation of the mitigation circuits""" + return ["0" * self.num_qubits, "1" * self.num_qubits] diff --git a/releasenotes/notes/readout-mitigation-experiment-4ea5392ee955a54c.yaml b/releasenotes/notes/readout-mitigation-experiment-4ea5392ee955a54c.yaml new file mode 100644 index 0000000000..d630cd3637 --- /dev/null +++ b/releasenotes/notes/readout-mitigation-experiment-4ea5392ee955a54c.yaml @@ -0,0 +1,5 @@ +--- +features: + - | + Added a new experiment (`qiskit_experiments.library import ReadoutMitigationExperiment`) + for calibrating readout mitigators. diff --git a/test/test_mitigation.py b/test/test_mitigation.py new file mode 100644 index 0000000000..92acaa5a7b --- /dev/null +++ b/test/test_mitigation.py @@ -0,0 +1,133 @@ +# This code is part of Qiskit. +# +# (C) Copyright IBM 2022. +# +# This code is licensed under the Apache License, Version 2.0. You may +# obtain a copy of this license in the LICENSE.txt file in the root directory +# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0. +# +# Any modifications or derivative works of this code must retain this +# copyright notice, and modified files need to carry a notice indicating +# that they have been altered from the originals. + +""" +A Tester for the Readout Mitigation experiment +""" + + +import unittest +from test.base import QiskitExperimentsTestCase +import numpy as np +from qiskit.quantum_info.operators.predicates import matrix_equal +from qiskit_experiments.library import ReadoutMitigationExperiment +from qiskit_experiments.framework import ExperimentData + + +class TestMitigation(QiskitExperimentsTestCase): + """Test ReadoutMitigationExperiment""" + + def test_local_analysis(self): + """Tests local mitigator generation from experimental data""" + qubits = [0, 1, 2] + run_data = [ + { + "counts": {"000": 986, "010": 10, "100": 16, "001": 12}, + "metadata": {"label": "000"}, + "shots": 1024, + }, + { + "counts": {"111": 930, "110": 39, "011": 24, "101": 29, "010": 1, "100": 1}, + "metadata": {"label": "111"}, + "shots": 1024, + }, + ] + expected_assignment_matrices = [ + np.array([[0.98828125, 0.04003906], [0.01171875, 0.95996094]]), + np.array([[0.99023438, 0.02929688], [0.00976562, 0.97070312]]), + np.array([[0.984375, 0.02441406], [0.015625, 0.97558594]]), + ] + run_meta = {"physical_qubits": qubits} + expdata = ExperimentData() + expdata.add_data(run_data) + expdata._metadata = run_meta + exp = ReadoutMitigationExperiment(qubits) + result = exp.analysis.run(expdata) + mitigator = result.analysis_results(0).value + + self.assertEqual(len(qubits), mitigator._num_qubits) + self.assertEqual(qubits, mitigator._qubits) + self.assertTrue(matrix_equal(expected_assignment_matrices, mitigator._assignment_mats)) + + def test_correlated_analysis(self): + """Tests correlated mitigator generation from experimental data""" + qubits = [0, 2, 3] + run_data = [ + { + "counts": {"000": 989, "010": 12, "100": 7, "001": 15, "101": 1}, + "metadata": {"label": "000"}, + "shots": 1024, + }, + { + "counts": {"001": 971, "101": 15, "000": 36, "011": 2}, + "metadata": {"label": "001"}, + "shots": 1024, + }, + { + "counts": {"000": 30, "010": 965, "110": 15, "011": 11, "001": 2, "100": 1}, + "metadata": {"label": "010"}, + "shots": 1024, + }, + { + "counts": {"011": 955, "111": 15, "010": 26, "001": 27, "110": 1}, + "metadata": {"label": "011"}, + "shots": 1024, + }, + { + "counts": {"100": 983, "101": 8, "110": 13, "000": 20}, + "metadata": {"label": "100"}, + "shots": 1024, + }, + { + "counts": {"101": 947, "001": 34, "100": 32, "111": 11}, + "metadata": {"label": "101"}, + "shots": 1024, + }, + { + "counts": {"100": 26, "110": 965, "010": 21, "111": 11, "000": 1}, + "metadata": {"label": "110"}, + "shots": 1024, + }, + { + "counts": {"111": 938, "011": 23, "110": 35, "101": 27, "100": 1}, + "metadata": {"label": "111"}, + "shots": 1024, + }, + ] + + expected_assignment_matrix = np.array( + [ + [0.96582031, 0.03515625, 0.02929688, 0.0, 0.01953125, 0.0, 0.00097656, 0.0], + [0.01464844, 0.94824219, 0.00195312, 0.02636719, 0.0, 0.03320312, 0.0, 0.0], + [0.01171875, 0.0, 0.94238281, 0.02539062, 0.0, 0.0, 0.02050781, 0.0], + [0.0, 0.00195312, 0.01074219, 0.93261719, 0.0, 0.0, 0.0, 0.02246094], + [0.00683594, 0.0, 0.00097656, 0.0, 0.95996094, 0.03125, 0.02539062, 0.00097656], + [0.00097656, 0.01464844, 0.0, 0.0, 0.0078125, 0.92480469, 0.0, 0.02636719], + [0.0, 0.0, 0.01464844, 0.00097656, 0.01269531, 0.0, 0.94238281, 0.03417969], + [0.0, 0.0, 0.0, 0.01464844, 0.0, 0.01074219, 0.01074219, 0.91601562], + ] + ) + run_meta = {"physical_qubits": qubits} + expdata = ExperimentData() + expdata.add_data(run_data) + expdata._metadata = run_meta + exp = ReadoutMitigationExperiment(qubits, method="correlated") + result = exp.analysis.run(expdata) + mitigator = result.analysis_results(0).value + + self.assertEqual(len(qubits), mitigator._num_qubits) + self.assertEqual(qubits, mitigator._qubits) + self.assertTrue(matrix_equal(expected_assignment_matrix, mitigator.assignment_matrix())) + + +if __name__ == "__main__": + unittest.main()