{
  "cells": [
    {
      "cell_type": "markdown",
      "id": "70810e0b-08f9-48ae-beb5-8c863819000b",
      "metadata": {},
      "source": [
        "---\n",
        "title: Readout error mitigation for the Sampler primitive using M3\n",
        "description: Use the M3 readout mitigation addon with the Sampler primitive\n",
        "---\n",
        "\n",
        "{/* cspell:ignore braket, Zgrm, newcommand, probs, quasis, topten */}\n",
        "\n",
        "# Readout error mitigation for the Sampler primitive using M3\n",
        "\n",
        "*Usage estimate: under one minute on a Heron r2 processor (NOTE: This is an estimate only. Your runtime might vary.)*\n",
        "\n",
        "## Background\n",
        "\n",
        "Unlike the Estimator primitive, the Sampler primitive does not have built-in support for error mitigation.\n",
        "Several of the methods supported by the Estimator are specifically designed for expectation values, and hence are not applicable to the Sampler primitive. An exception is readout error mitigation, which is a highly effective method that is also applicable to the Sampler primitive.\n",
        "\n",
        "The [M3 Qiskit addon](https://qiskit.github.io/qiskit-addon-mthree/) implements an efficient method for readout error mitigation. This tutorial explains how to use the M3 Qiskit addon to mitigate readout error for the Sampler primitive.\n",
        "\n",
        "### What is readout error?\n",
        "\n",
        "Immediately before measurement, the state of a qubit register is\n",
        "described by a superposition of computational basis states,\n",
        "or by a density matrix.\n",
        "Measurement of the qubit register into a classical bit register then proceeds in two steps.\n",
        "First the quantum measurement proper is performed.\n",
        "This means that the state of the qubit register\n",
        "is projected onto a single basis state that is characterized\n",
        "by a string of $1$s and $0$s.\n",
        "The second step consists of reading the bitstring characterizing this basis state\n",
        "and writing it into classical computer memory.\n",
        "We call this step *readout*.\n",
        "It turns out that the second step (readout) incurs more error than the first step (projection onto basis states).\n",
        "This makes sense when you recall that readout requires detecting a microscopic\n",
        "quantum state and amplifying it to the macroscopic realm. A readout resonator is coupled to\n",
        "the (transmon) qubit, thereby experiencing a very small frequency shift. A microwave pulse\n",
        "is then bounced off of the resonator, in turn experiencing small changes in its\n",
        "characteristics.  The reflected pulse is then amplified and analyzed.  This is a delicate\n",
        "process and is subject to a host of errors.\n",
        "\n",
        "The important point is that, while both quantum measurement and readout are subject to error, the\n",
        "latter incurs the dominant error, called readout error, which is the focus in this tutorial.\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "ba74fc81-30ac-4365-9ba6-a891ba082bd6",
      "metadata": {
        "jp-MarkdownHeadingCollapsed": true
      },
      "source": [
        "### Theoretical background\n",
        "\n",
        "If the sampled bitstring (stored in classical memory) differs from the bitstring characterizing\n",
        "the projected quantum state, we say that a readout error has occurred.\n",
        "These errors are observed to be random and uncorrelated from sample to sample.\n",
        "It has proven useful to model readout error as a *noisy classical channel*.\n",
        "That is, for every pair of\n",
        "bitstrings $i$ and $j$, there is a fixed probability that a true value of $j$ will\n",
        "be incorrectly read as $i$.\n",
        "\n",
        "More precisely, for every pair of bitstrings $(i, j)$, there is a (conditional) probability ${M}_{i,j}$\n",
        "that $i$ is read, given that the true value is $j.$\n",
        "That is,\n",
        "\n",
        "$$\n",
        "    {M}_{i,j} =  \\Pr(\\text{readout value is } i | \\text{true value is } j)\n",
        "    \\text{ for } i,j \\in (0,...,2^n - 1), \\tag{1}\n",
        "$$\n",
        "\n",
        "where $n$ is the number of bits in the readout register.\n",
        "For concreteness, we assume that $i$ is a decimal integer whose binary representation is\n",
        "the bitstring that labels the computational basis states.\n",
        "We call the $2^n \\times 2^n$ matrix ${M}$ the *assignment matrix*.\n",
        "For fixed true value $j$, summing the probability over all noisy outcomes $i$ must give $1$. That is\n",
        "\n",
        "$$\n",
        "    \\sum_{i=0}^{2^n - 1} {M}_{i,j} = 1 \\text{ for all } j\n",
        "$$\n",
        "\n",
        "A matrix with no negative entries that satisfies (1) is called\n",
        "*left-stochastic*.\n",
        "A left-stochastic matrix is also called *column-stochastic* because each of its columns sums to $1$.\n",
        "We experimentally determine approximate values for each element ${M}_{i,j}$ by\n",
        "repeatedly preparing each basis state $|j \\rangle$ and then computing the frequencies\n",
        "of the occurrence of sampled bitstrings.\n",
        "\n",
        "If an experiment involves estimating a probability distribution over output bitstrings by repeated sampling,\n",
        "then we can use ${M}$ to mitigate readout error at the level of the distribution.\n",
        "The first step is to repeat a fixed circuit of interest many times,\n",
        "creating a histogram of sampled bitstrings.\n",
        "The normalized histogram is the measured probability distribution over\n",
        "the $2^n$ possible bitstrings, which we denote by ${\\tilde{p}} \\in \\mathbb{R}^{2^n}$.\n",
        "The (estimated) probability ${{\\tilde{p}}}_i$ of sampling bitstring $i$\n",
        "is equal to the sum over all true bitstrings $j$, each weighted by\n",
        "the probability that it is mistaken for $i$.\n",
        "This statement in matrix form is\n",
        "\n",
        "$$\n",
        "    {\\tilde{p}} = {M} {\\vec{p}}, \\tag{2},\n",
        "$$\n",
        "\n",
        "where ${\\vec{p}}$ is the true distribution. In words, the readout error has the effect of multiplying\n",
        "the ideal distribution over bitstrings ${\\vec{p}}$ by the assignment matrix ${M}$ to\n",
        "produce the observed distribution ${\\tilde{p}}$.\n",
        "We have measured ${\\tilde{p}}$ and ${M}$, but have no direct access to ${\\vec{p}}$. In principle, we will\n",
        "obtain the true distribution of bitstrings for our circuit\n",
        "by solving equation (2) for ${\\vec{p}}$ numerically.\n",
        "\n",
        "Before we move on, it's worth noting some important features of this naive approach.\n",
        "\n",
        "* In practice, equation (2) is not solved by inverting ${M}$. Linear algebra\n",
        "  routines in software libraries employ methods that are more stable, accurate, and efficient.\n",
        "* When estimating ${M}$, we assumed that only readout errors occurred. In particular,\n",
        "  we assume there were no state preparation and quantum measurement errors —\n",
        "  or at least that they were otherwise mitigated.\n",
        "  To the extent that this is a good assumption, ${M}$ really represents\n",
        "  only readout error. But when we *use* ${M}$ to correct a measured distribution\n",
        "  over bitstrings, we make no such assumption. In fact, we expect an interesting\n",
        "  circuit to introduce noise, for instance, gate errors. The \"true\" distribution\n",
        "  still includes effects from any errors that are not otherwise mitigated.\n",
        "\n",
        "This method, while useful in some circumstances, suffers from a few limitations.\n",
        "\n",
        "The space and time resources needed to estimate ${M}$ grow exponentially in $n$:\n",
        "\n",
        "* The estimation of ${M}$ and ${\\tilde{p}}$ is subject to statistical error due to finite sampling.\n",
        "  This noise can be made as small as desired\n",
        "  at the cost of more shots (up to the timescale of drifting hardware parameters\n",
        "  that result in systematic errors in ${M}$).\n",
        "  However, if no assumptions are made on the bitstrings observed\n",
        "  when performing mitigation, the number of shots required to estimate ${M}$ grows\n",
        "  at least exponentially in $n$.\n",
        "* ${M}$ is a $2^n \\times 2^n$ matrix.\n",
        "  When $n>10$, the amount of memory required to store ${M}$ is\n",
        "  greater than the memory available in a powerful laptop.\n",
        "\n",
        "Further limitations are:\n",
        "\n",
        "* The recovered distribution ${\\vec{p}}$ may have one\n",
        "  or more negative probabilities (while still summing to one). One solution\n",
        "  is to minimize $||{M} {\\vec{p}} - {\\tilde{p}}||^2$ subject to the constraint that\n",
        "  each entry in ${\\vec{p}}$ be non-negative. However, the runtime of such a\n",
        "  method is orders of magnitude longer than directly solving equation (2).\n",
        "* This mitigation procedure works on the level of a probability distribution\n",
        "  over bitstrings. In particular, it cannot correct an error in an individual\n",
        "  observed bitstring.\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "c98d8409-66ad-452f-b87a-d8df886d77d4",
      "metadata": {},
      "source": [
        "### Qiskit M3 addon: Scaling to longer bitstrings\n",
        "\n",
        "Solving equation (2) using standard numeric linear algebra routines is limited to bitstrings no longer than about 10 bits. M3, however, can handle much longer bitstrings. Two key properties of M3 that make this possible are:\n",
        "\n",
        "* Correlations in readout error of order three and higher among collections of bits\n",
        "  are assumed to be negligible and are ignored. In principle, at the cost of more shots,\n",
        "  one could estimate higher correlations as well.\n",
        "* Rather than constructing ${M}$ explicitly, we use a much smaller effective matrix that records\n",
        "  probabilities only for bitstrings collected when constructing ${\\tilde{p}}$.\n",
        "\n",
        "At a high level, the procedure works as follows.\n",
        "\n",
        "First, we construct building blocks from which we can construct a simplified, effective, description of ${M}$.\n",
        "Then, we repeatedly run the circuit of interest and collect bitstrings which we use to construct\n",
        "both ${\\tilde{p}}$ and, with the help of the building blocks, an effective ${M}$.\n",
        "\n",
        "More precisely,\n",
        "\n",
        "* Single-qubit assignment matrices are estimated for each qubit. To do this, we repeatedly\n",
        "  prepare the qubit register in the all-zero state $|0 ... 0 \\rangle$ and then in the all-one\n",
        "  state $|1 ... 1 \\rangle$, and record the probability for each qubit that it is read\n",
        "  incorrectly.\n",
        "* Correlations of order three and higher are assumed to be negligible and are ignored.\n",
        "\n",
        "  Instead we construct a number $n$ of $2 \\times 2$ single-qubit\n",
        "  assignment matrices, and a number $n(n-1)/2$ of $4 \\times 4$ two-qubit assignment\n",
        "  matrices. These one- and two-qubit assignment matrices are stored for later\n",
        "  use.\n",
        "* After repeatedly sampling a circuit to construct ${\\tilde{p}}$,\n",
        "  we construct an effective approximation to ${M}$ using only\n",
        "  bitstrings that are sampled when constructing ${\\tilde{p}}$. This effective matrix\n",
        "  is built using the single- and two-qubit matrices described in the previous item.\n",
        "  The linear dimension of this matrix is at most of the order of the number\n",
        "  of shots used in constructing ${\\tilde{p}}$, which is much smaller than\n",
        "  the dimension $2^n$ of the full assignment matrix ${M}$ .\n",
        "\n",
        "For technical details on M3, you can refer to [*Scalable Mitigation of Measurement Errors on Quantum Computers*](https://journals.aps.org/prxquantum/abstract/10.1103/PRXQuantum.2.040326).\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "14f04ac4-5142-4436-93e1-f712d47dde1a",
      "metadata": {},
      "source": [
        "### Application of M3 to a quantum algorithm\n",
        "\n",
        "We'll apply M3's readout mitigation to the hidden shift problem. The hidden shift problem, and closely related problems such as the [hidden subgroup problem](https://en.wikipedia.org/wiki/Hidden_subgroup_problem), were originally conceived in a fault-tolerant setting (more precisely, before fault-tolerant QPUs were proved to be possible!). But they are studied with available processors as well. An example of algorithmic exponential speedup obtained for a variant of the hidden shift problem obtained on 127-qubit IBM® QPUs can be found in [this paper](https://journals.aps.org/prx/accepted/a9074K06A8e1590147da9c69f8c4b64c28247be5a) ([arXiv version](https://arxiv.org/abs/2401.07934)).\n",
        "\n",
        "In the following, all arithmetic is Boolean.\n",
        "That is, for $a, b \\in \\mathbb{Z}_2 = \\{0, 1\\}$, addition, $a + b$ is the logical XOR function.\n",
        "Furthermore, multiplication $a \\times b$ (or $a b$) is the logical AND function. For $x, y \\in \\{0, 1\\}^n$,\n",
        "$x + y$ is defined by bitwise application of XOR.\n",
        "The dot product $\\cdot: {\\mathbb{Z}_2^n} \\rightarrow \\mathbb{Z}_2$ is defined\n",
        "by $x \\cdot y = \\sum_i x_i y_i$.\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "81cb34b4-566b-4ee1-8309-7464a5fd3bb7",
      "metadata": {},
      "source": [
        "#### Hadamard operator and Fourier transform\n",
        "\n",
        "In implementing quantum algorithms, it is very common to use the Hadamard operator as a Fourier transform.\n",
        "The computational basis states are sometimes called *classical states*. They stand in\n",
        "a one-to-one relation to the classical bitstrings.\n",
        "The $n$-qubit Hadamard operator on classical states can be viewed as a Fourier transform on the Boolean hypercube:\n",
        "\n",
        "$$\n",
        "H^{\\otimes n} =  \\frac{1}{\\sqrt{2^n}} \\sum_{x,y \\in {\\mathbb{Z}_2^n}} (-1)^{x \\cdot y} {|{y}\\rangle}{\\langle{x}|}.\n",
        "$$\n",
        "\n",
        "Consider a state ${|{s}\\rangle}$ corresponding to fixed bitstring $s$.\n",
        "Applying $H^{\\otimes n}$, and using ${\\langle {x}|{s}\\rangle} = \\delta_{x,s}$,\n",
        "we see that the Fourier transform of ${|{s}\\rangle}$ can be written as\n",
        "\n",
        "$$\n",
        "   H^{\\otimes n} {|{s}\\rangle} =  \\frac{1}{\\sqrt{2^n}} \\sum_{y \\in {\\mathbb{Z}_2^n}} (-1)^{s \\cdot y} {|{y}\\rangle}.\n",
        "$$\n",
        "\n",
        "The Hadamard is its own inverse, that is,\n",
        "$H^{\\otimes n} H^{\\otimes n} = (H H)^{\\otimes n} = I^{\\otimes n}$.\n",
        "Thus, the inverse Fourier transform is the same operator, $H^{\\otimes n}$.\n",
        "Explicitly, we have,\n",
        "\n",
        "$$\n",
        "  {|{s}\\rangle} =  H^{\\otimes n} H^{\\otimes n} {|{s}\\rangle}  =  H^{\\otimes n} \\frac{1}{\\sqrt{2^n}} \\sum_{y \\in {\\mathbb{Z}_2^n}} (-1)^{s \\cdot y} {|{y}\\rangle}.\n",
        "$$\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "d10fce93-3e08-4a27-b641-eda5e7b05ce0",
      "metadata": {},
      "source": [
        "#### The hidden shift problem\n",
        "\n",
        "We consider a simple example of a *hidden shift problem*.\n",
        "The problem is to identify a constant shift in the input to a function.\n",
        "The function we consider is the dot product. It is the simplest member\n",
        "of a large class of functions that admit a quantum speedup for the hidden shift\n",
        "problem via techniques similar to those presented below.\n",
        "\n",
        "Let $x,y \\in {\\mathbb{Z}_2^m}$ be bitstrings of length $m$.\n",
        "We define ${f}: {\\mathbb{Z}_2^m} \\times {\\mathbb{Z}_2^m} \\rightarrow \\{-1,1\\}$ by\n",
        "\n",
        "$$\n",
        "  {f}(x, y) = (-1)^{x \\cdot y}.\n",
        "$$\n",
        "\n",
        "Let $a,b \\in {\\mathbb{Z}_2^m}$ be fixed bitstrings of length $m$.\n",
        "We furthermore define $g: {\\mathbb{Z}_2^m} \\times {\\mathbb{Z}_2^m} \\rightarrow \\{-1,1\\}$ by\n",
        "\n",
        "$$\n",
        "  g(x, y) = {f}(x+a, y+b) = (-1)^{(x+a) \\cdot (y+b)},\n",
        "$$\n",
        "\n",
        "where $a$ and $b$ are (hidden) parameters.\n",
        "We are given two black boxes, one implementing $f$, and the other $g$.\n",
        "We suppose that we know that they compute the functions defined above, except that we know\n",
        "neither $a$ nor $b$. The game is to determine the hidden bitstrings (shifts)\n",
        "$a$ and $b$ by making queries to $f$ and $g$. It's clear that if we play the game classically,\n",
        "we need $O(2m)$ queries to determine $a$ and $b$. For example, we can query $g$ with all pairs of strings such that one element of the pair is all zeros, and the other element has exactly one element set to $1$.\n",
        "On each query, we learn one element of either $a$ or $b$.\n",
        "However, we will see that, if the black boxes are implemented as quantum circuits, we can\n",
        "determine $a$ and $b$ with a single query to each of $f$ and $g$.\n",
        "\n",
        "In the context of algorithmic complexity, a black box is called an *oracle*.\n",
        "In addition to being opaque, an oracle has the property that it consumes the input and\n",
        "produces the output instantly, adding nothing to the complexity budget of the algorithm\n",
        "in which it is embedded. In fact, in the case at hand, the oracles implementing $f$ and\n",
        "$g$ will be seen to be efficient.\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "d5e4c1e2-d6f5-4093-afdb-f97e9f062179",
      "metadata": {
        "jp-MarkdownHeadingCollapsed": true
      },
      "source": [
        "#### Quantum circuits for $f$ and $g$\n",
        "\n",
        "We need the following ingredients in order to implement $f$ and $g$ as quantum circuits.\n",
        "\n",
        "For single-qubit classical states ${|{x_1}\\rangle}, {|{y_1}\\rangle}$, with $x_1,y_1 \\in \\mathbb{Z}_2$,\n",
        "the controlled-$Z$ gate ${CZ}$ may be written as\n",
        "\n",
        "$$\n",
        "{CZ} {|{x_1}\\rangle}{|{y_1}\\rangle}{x_1} = (-1)^{x_1 y_1} {|{x_1}\\rangle}{x_1}{|{y_1}\\rangle}.\n",
        "$$\n",
        "\n",
        "We will operate with $m$ CZ gates, one on $(x_1, y_1)$, and one on $(x_2, y_2)$, and so on, through $(x_m, y_m)$.\n",
        "We call this operator ${CZ}_{x,y}$.\n",
        "\n",
        "$U_f = {CZ}_{x,y}$ is a quantum version of ${f} = {f}(x,y)$:\n",
        "\n",
        "$$\n",
        "%\\CZ_{x,y} {|#1\\rangle}{z} =\n",
        "U_f {|{x}\\rangle}{|{y}\\rangle} = {CZ}_{x,y} {|{x}\\rangle}{|{y}\\rangle} = (-1)^{x \\cdot y}  {|{x}\\rangle}{|{y}\\rangle}.\n",
        "$$\n",
        "\n",
        "We also need to implement a bitstring shift.\n",
        "We denote the operator on the $x$ register $X^{a_1}\\cdots X^{a_m}$ by $X_a$\n",
        "and likewise on the $y$ register $X_b =  X^{b_1}\\cdots X^{b_m}$.\n",
        "These operators apply $X$ wherever a single bit is $1$, and the identity $I$ wherever it is $0$.\n",
        "Then we have\n",
        "\n",
        "$$\n",
        " X_a X_b  {|{x}\\rangle}{|{y}\\rangle} = {|{x+a}\\rangle}{|{y+b}\\rangle}.\n",
        "$$\n",
        "\n",
        "The second black box $g$ is implemented by the unitary $U_g$, given by\n",
        "\n",
        "$$\n",
        "%U_g {|{x}\\rangle}{|{y}\\rangle} = X_aX_b \\CZ_{x,y} X_aX_b {|{x}\\rangle}{|{y}\\rangle}.\n",
        "U_g = X_aX_b {CZ}_{x,y} X_aX_b.\n",
        "$$\n",
        "\n",
        "To see this, we apply the operators from right to left to the state ${|{x}\\rangle}{|{y}\\rangle}$.\n",
        "First\n",
        "\n",
        "$$\n",
        " X_a X_b  {|{x}\\rangle}{|{y}\\rangle} = {|{x+a}\\rangle}{|{y+b}\\rangle}.\n",
        "$$\n",
        "\n",
        "Then,\n",
        "\n",
        "$$\n",
        "  {CZ}_{x,y}  {|{x+a}\\rangle}{|{y+b}\\rangle} = (-1)^{(x+a)\\cdot (y+b)} {|{x+a}\\rangle}{|{y+b}\\rangle}.\n",
        "$$\n",
        "\n",
        "Finally,\n",
        "\n",
        "$$\n",
        "  X^a X^b (-1)^{(x+a)\\cdot (y+b)} {|{x+a}\\rangle}{|{y+b}\\rangle} = (-1)^{(x+a)\\cdot (y+b)} {|{x}\\rangle}{|{y}\\rangle},\n",
        "$$\n",
        "\n",
        "which is indeed the quantum version of $f(x+a, y+b)$.\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "de0bb5f9-517e-4bea-91ab-f65e012e8ae9",
      "metadata": {
        "editable": true,
        "jp-MarkdownHeadingCollapsed": true,
        "slideshow": {
          "slide_type": ""
        },
        "tags": []
      },
      "source": [
        "#### The hidden shift algorithm\n",
        "\n",
        "Now we put the pieces together to solve the hidden shift problem.\n",
        "We begin by applying Hadamards to the registers initialized to the all-zero state.\n",
        "\n",
        "$$\n",
        "H^{\\otimes 2m} = H^{\\otimes m} \\otimes H^{\\otimes m} {{|{0}\\rangle}^{\\otimes m}}{{|{0}\\rangle}^{\\otimes m}} = \\frac{1}{\\sqrt{2^{2m}}} \\sum_{x, y \\in {\\mathbb{Z}_2^m}} (-1)^{x \\cdot y} {|{x}\\rangle}{|{y}\\rangle}.\n",
        "$$\n",
        "\n",
        "Next, we query the oracle $g$ to arrive at\n",
        "\n",
        "$$\n",
        "U_g H^{\\otimes 2m} {{|{0}\\rangle}^{\\otimes m}}{{|{0}\\rangle}^{\\otimes m}}\n",
        "= \\frac{1}{\\sqrt{2^{2m}}} \\sum_{x, y \\in {\\mathbb{Z}_2^m}} (-1)^{(x+a) \\cdot (y+b)} {|{x}\\rangle}{|{y}\\rangle}\n",
        "$$\n",
        "\n",
        "$$\n",
        "\\approx \\frac{1}{\\sqrt{2^{2m}}} \\sum_{x, y \\in {\\mathbb{Z}_2^m}} (-1)^{x \\cdot y + x \\cdot b + y \\cdot a} {|{x}\\rangle}{|{y}\\rangle}.\n",
        "$$\n",
        "\n",
        "In the last line, we omitted the constant global phase factor $(-1)^{a \\cdot b}$,\n",
        "and denote equality up to a phase by $\\approx$.\n",
        "Next, applying the oracle $f$ introduces another factor of $(-1)^{x \\cdot y}$, canceling the one already\n",
        "present. We then have:\n",
        "\n",
        "$$\n",
        "U_f U_g H^{\\otimes 2m} {{|{0}\\rangle}^{\\otimes m}}{{|{0}\\rangle}^{\\otimes m}}\n",
        "\\approx \\frac{1}{\\sqrt{2^{2m}}} \\sum_{x, y \\in {\\mathbb{Z}_2^m}} (-1)^{x \\cdot b + y \\cdot a} {|{x}\\rangle}{|{y}\\rangle}.\n",
        "$$\n",
        "\n",
        "The final step is to apply the inverse Fourier transform, $H^{\\otimes 2m} = H^{\\otimes m} \\otimes H^{\\otimes m}$,\n",
        "resulting in\n",
        "\n",
        "$$\n",
        "H^{\\otimes 2m} U_f U_g  H^{\\otimes 2m} {{|{0}\\rangle}^{\\otimes m}}{{|{0}\\rangle}^{\\otimes m}}\n",
        "\\approx {|{b}\\rangle}{|{a}\\rangle}.\n",
        "$$\n",
        "\n",
        "The circuit is finished. In the absence of noise, sampling the quantum registers will\n",
        "return the bitstrings $b, a$ with probability $1$.\n",
        "\n",
        "The Boolean inner product is an example of the so-called bent functions.\n",
        "We will not define bent functions here\n",
        "but merely note that they\n",
        "\"are maximally resistant against attacks that seek to exploit a dependence of\n",
        "the outputs on some linear subspace of the inputs.\"\n",
        "This quote is from the article [*Quantum algorithms for highly non-linear Boolean functions*](https://arxiv.org/abs/0811.3208), which\n",
        "gives efficient hidden shift algorithms for several classes of bent functions.\n",
        "The algorithm in this tutorial appears in Section 3.1 of the article.\n",
        "\n",
        "In the more general case, the circuit for finding a hidden shift $s \\in \\mathbb{Z}^n$ is\n",
        "\n",
        "$$\n",
        " H^{\\otimes n} U_{\\tilde{f}}  H^{\\otimes n} U_g  H^{\\otimes n} {|{0}\\rangle}^{\\otimes n} = {|{s}\\rangle}.\n",
        "$$\n",
        "\n",
        "In the general case, $f$ and $g$ are functions of a single variable.\n",
        "Our example of the inner product has this form if we let $f(x, y) \\to f(z)$,\n",
        "with $z$ equal to the concatenation of $x$ and $y$, and $s$ equal to the concatenation\n",
        "of $a$ and $b$.\n",
        "The general case requires exactly two oracles: One oracle for $g$ and one for $\\tilde{f}$,\n",
        "where the latter is a function known as the *dual* of the bent function $f$.\n",
        "The inner product function has the self-dual property $\\tilde{f}=f$.\n",
        "\n",
        "In our circuit for the hidden shift on the inner product we omitted the middle layer\n",
        "of Hadamards that appears in the circuit for the general case. While in the general case\n",
        "this layer is necessary, we saved a bit of depth by omitting it, at the expense of a bit\n",
        "of post-processing because the output is ${|{b}\\rangle}{|{a}\\rangle}$ instead of the desired ${|{a}\\rangle}{|{b}\\rangle}$.\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "ccc0aec5-63d3-4a43-8fad-fea56ea24ec7",
      "metadata": {},
      "source": [
        "## Requirements\n",
        "\n",
        "Before starting this tutorial, ensure that you have the following installed:\n",
        "\n",
        "* Qiskit SDK v2.1 or later, with [visualization](/docs/api/qiskit/visualization) support\n",
        "* Qiskit Runtime v0.41 or later (`pip install qiskit-ibm-runtime`)\n",
        "* M3 Qiskit addon v3.0 (`pip install mthree`)\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "c090b8dd-754f-4390-93e5-663f818e61d0",
      "metadata": {
        "jp-MarkdownHeadingCollapsed": true
      },
      "source": [
        "## Setup\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "233dfc67-d280-4827-910d-8e7170ee95ad",
      "metadata": {},
      "outputs": [],
      "source": [
        "from collections.abc import Iterator, Sequence\n",
        "from random import Random\n",
        "from qiskit.circuit import (\n",
        "    CircuitInstruction,\n",
        "    QuantumCircuit,\n",
        "    QuantumRegister,\n",
        "    Qubit,\n",
        ")\n",
        "from qiskit.circuit.library import CZGate, HGate, XGate\n",
        "from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager\n",
        "from qiskit_ibm_runtime import QiskitRuntimeService\n",
        "import timeit\n",
        "import matplotlib.pyplot as plt\n",
        "from qiskit_ibm_runtime import SamplerV2 as Sampler\n",
        "import mthree"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "5cf89647-261d-4af1-af54-75548345df6e",
      "metadata": {},
      "source": [
        "## Step 1: Map classical inputs to a quantum problem\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "940cc113-5cf3-4688-9747-93268be10088",
      "metadata": {},
      "source": [
        "First, we write the functions to implement the hidden shift problem as a `QuantumCircuit`.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "675fe72b-6a86-42ea-a65f-b29a60725ecb",
      "metadata": {
        "editable": true,
        "scrolled": true,
        "slideshow": {
          "slide_type": ""
        },
        "tags": []
      },
      "outputs": [],
      "source": [
        "def apply_hadamards(qubits: Sequence[Qubit]) -> Iterator[CircuitInstruction]:\n",
        "    \"\"\"Apply a Hadamard gate to every qubit.\"\"\"\n",
        "    for q in qubits:\n",
        "        yield CircuitInstruction(HGate(), [q], [])\n",
        "\n",
        "\n",
        "def apply_shift(\n",
        "    qubits: Sequence[Qubit], shift: int\n",
        ") -> Iterator[CircuitInstruction]:\n",
        "    \"\"\"Apply X gates where the bits of the shift are equal to 1.\"\"\"\n",
        "    for i, q in zip(range(shift.bit_length()), qubits):\n",
        "        if shift >> i & 1:\n",
        "            yield CircuitInstruction(XGate(), [q], [])\n",
        "\n",
        "\n",
        "def oracle_f(qubits: Sequence[Qubit]) -> Iterator[CircuitInstruction]:\n",
        "    \"\"\"Apply the f oracle.\"\"\"\n",
        "    for i in range(0, len(qubits) - 1, 2):\n",
        "        yield CircuitInstruction(CZGate(), [qubits[i], qubits[i + 1]])\n",
        "\n",
        "\n",
        "def oracle_g(\n",
        "    qubits: Sequence[Qubit], shift: int\n",
        ") -> Iterator[CircuitInstruction]:\n",
        "    \"\"\"Apply the g oracle.\"\"\"\n",
        "    yield from apply_shift(qubits, shift)\n",
        "    yield from oracle_f(qubits)\n",
        "    yield from apply_shift(qubits, shift)\n",
        "\n",
        "\n",
        "def determine_hidden_shift(\n",
        "    qubits: Sequence[Qubit], shift: int\n",
        ") -> Iterator[CircuitInstruction]:\n",
        "    \"\"\"Determine the hidden shift.\"\"\"\n",
        "    yield from apply_hadamards(qubits)\n",
        "    yield from oracle_g(qubits, shift)\n",
        "    # We omit this layer in exchange for post processing\n",
        "    # yield from apply_hadamards(qubits)\n",
        "    yield from oracle_f(qubits)\n",
        "    yield from apply_hadamards(qubits)\n",
        "\n",
        "\n",
        "def run_hidden_shift_circuit(n_qubits, rng):\n",
        "    hidden_shift = rng.getrandbits(n_qubits)\n",
        "\n",
        "    qubits = QuantumRegister(n_qubits, name=\"q\")\n",
        "    circuit = QuantumCircuit.from_instructions(\n",
        "        determine_hidden_shift(qubits, hidden_shift), qubits=qubits\n",
        "    )\n",
        "    circuit.measure_all()\n",
        "    # Format the hidden shift as a string.\n",
        "    hidden_shift_string = format(hidden_shift, f\"0{n_qubits}b\")\n",
        "    return (circuit, hidden_shift, hidden_shift_string)\n",
        "\n",
        "\n",
        "def display_circuit(circuit):\n",
        "    return circuit.remove_final_measurements(inplace=False).draw(\n",
        "        \"mpl\", idle_wires=False, scale=0.5, fold=-1\n",
        "    )"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "6960aa96-f612-40b3-9f26-3d3afb07262b",
      "metadata": {},
      "source": [
        "We'll start with a small example:\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 2,
      "id": "8297843e-00c3-4bb5-9d33-a7e558d1698c",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Hidden shift string 011010\n"
          ]
        },
        {
          "data": {
            "text/plain": [
              "<Image src=\"/docs/images/tutorials/readout-error-mitigation-sampler/extracted-outputs/8297843e-00c3-4bb5-9d33-a7e558d1698c-1.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "execution_count": 2,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "n_qubits = 6\n",
        "random_seed = 12345\n",
        "rng = Random(random_seed)\n",
        "circuit, hidden_shift, hidden_shift_string = run_hidden_shift_circuit(\n",
        "    n_qubits, rng\n",
        ")\n",
        "\n",
        "print(f\"Hidden shift string {hidden_shift_string}\")\n",
        "\n",
        "display_circuit(circuit)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "c136f784-e105-4244-a483-b8221577afc3",
      "metadata": {
        "editable": true,
        "slideshow": {
          "slide_type": ""
        },
        "tags": []
      },
      "source": [
        "## Step 2: Optimize circuits for quantum hardware execution\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 3,
      "id": "ee4d384a-15e2-4a3f-a52a-def3d184c4fa",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "['shift 011010', 'n_qubits 6', 'seed = 12345']"
            ]
          },
          "execution_count": 3,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "job_tags = [\n",
        "    f\"shift {hidden_shift_string}\",\n",
        "    f\"n_qubits {n_qubits}\",\n",
        "    f\"seed = {random_seed}\",\n",
        "]\n",
        "job_tags"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "f2b77d93-c34a-43a4-b436-e7a25024a94a",
      "metadata": {
        "editable": true,
        "scrolled": true,
        "slideshow": {
          "slide_type": ""
        },
        "tags": []
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Using backend ibm_kingston\n"
          ]
        },
        {
          "data": {
            "text/plain": [
              "<Image src=\"/docs/images/tutorials/readout-error-mitigation-sampler/extracted-outputs/f2b77d93-c34a-43a4-b436-e7a25024a94a-1.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "execution_count": 4,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "# Uncomment this to run the circuits on a quantum computer on IBMCloud.\n",
        "service = QiskitRuntimeService()\n",
        "backend = service.least_busy(\n",
        "    operational=True, simulator=False, min_num_qubits=100\n",
        ")\n",
        "\n",
        "# from qiskit_ibm_runtime.fake_provider import FakeMelbourneV2\n",
        "# backend = FakeMelbourneV2()\n",
        "# backend.refresh(service)\n",
        "\n",
        "print(f\"Using backend {backend.name}\")\n",
        "\n",
        "\n",
        "def get_isa_circuit(circuit, backend):\n",
        "    pass_manager = generate_preset_pass_manager(\n",
        "        optimization_level=3, backend=backend, seed_transpiler=1234\n",
        "    )\n",
        "    isa_circuit = pass_manager.run(circuit)\n",
        "    return isa_circuit\n",
        "\n",
        "\n",
        "isa_circuit = get_isa_circuit(circuit, backend)\n",
        "display_circuit(isa_circuit)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "f2f426da-3601-4442-9593-16bb19abf91e",
      "metadata": {},
      "source": [
        "## Step 3: Execute circuits using Qiskit primitives\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "6ebaea63-e352-4d7c-a7ee-57a2f0a66306",
      "metadata": {},
      "outputs": [],
      "source": [
        "# submit job for solving the hidden shift problem using the Sampler primitive\n",
        "NUM_SHOTS = 50_000\n",
        "\n",
        "\n",
        "def run_sampler(backend, isa_circuit, num_shots):\n",
        "    sampler = Sampler(mode=backend)\n",
        "    sampler.options.environment.job_tags\n",
        "    pubs = [(isa_circuit, None, NUM_SHOTS)]\n",
        "    job = sampler.run(pubs)\n",
        "    return job\n",
        "\n",
        "\n",
        "def setup_mthree_mitigation(isa_circuit, backend):\n",
        "    # retrieve the final qubit mapping so mthree knows which qubits to calibrate\n",
        "    qubit_mapping = mthree.utils.final_measurement_mapping(isa_circuit)\n",
        "\n",
        "    # submit jobs for readout error calibration\n",
        "    mit = mthree.M3Mitigation(backend)\n",
        "    mit.cals_from_system(qubit_mapping, rep_delay=None)\n",
        "\n",
        "    return mit, qubit_mapping"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 6,
      "id": "9b76b943-fd71-4287-804b-980cacc078f8",
      "metadata": {},
      "outputs": [],
      "source": [
        "job = run_sampler(backend, isa_circuit, NUM_SHOTS)\n",
        "mit, qubit_mapping = setup_mthree_mitigation(isa_circuit, backend)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "674ad6e1-c496-44ca-9598-470edf4c15dc",
      "metadata": {},
      "source": [
        "## Step 4: Post-process and return results in classical format\n",
        "\n",
        "In the theoretical discussion above, we determined that for input $ab$, we expect output $ba$.\n",
        "An additional complication is that, in order to have a simpler (pre-transpiled) circuit, we inserted the required CZ gates between\n",
        "neighboring pairs of qubits. This amounts to interleaving the bitstrings $a$ and $b$ as $a1 b1 a2 b2 \\ldots$.\n",
        "The output string $ba$ will be interleaved in a similar way: $b1 a1 b2 a2 \\ldots$. The function `unscramble` below\n",
        "transforms the output string from $b1 a1 b2 a2 \\ldots$ to $a1 b1 a2 b2 \\ldots$ so that the input and output strings can be compared directly.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 7,
      "id": "aada9b41-0600-4814-a5a3-917d7fb171f9",
      "metadata": {},
      "outputs": [],
      "source": [
        "# retrieve bitstring counts\n",
        "def get_bitstring_counts(job):\n",
        "    result = job.result()\n",
        "    pub_result = result[0]\n",
        "    counts = pub_result.data.meas.get_counts()\n",
        "    return counts, pub_result"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 8,
      "id": "3930e290-89ff-404d-946b-cdf90119cb48",
      "metadata": {},
      "outputs": [],
      "source": [
        "counts, pub_result = get_bitstring_counts(job)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "138ef822-8e66-4f01-94f6-5a634465d2ac",
      "metadata": {},
      "source": [
        "The Hamming distance between two bitstrings is the number of indices at which the bits differ.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 9,
      "id": "3488f07f-7705-4ef2-b033-ebc2cbe19d45",
      "metadata": {},
      "outputs": [],
      "source": [
        "def hamming_distance(s1, s2):\n",
        "    weight = 0\n",
        "    for c1, c2 in zip(s1, s2):\n",
        "        (c1, c2) = (int(c1), int(c2))\n",
        "        if (c1 == 1 and c2 == 1) or (c1 == 0 and c2 == 0):\n",
        "            weight += 1\n",
        "\n",
        "    return weight"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 10,
      "id": "9e2f66e7-cad3-4893-9fde-ddaec077d938",
      "metadata": {
        "scrolled": true
      },
      "outputs": [],
      "source": [
        "# Replace string of form a1b1a2b2... with b1a1b2a1...\n",
        "# That is, reverse order of successive pairs of bits.\n",
        "def unscramble(bitstring):\n",
        "    ps = [bitstring[i : i + 2][::-1] for i in range(0, len(bitstring), 2)]\n",
        "    return \"\".join(ps)\n",
        "\n",
        "\n",
        "def find_hidden_shift_bitstring(counts, hidden_shift_string):\n",
        "    # convert counts to probabilities\n",
        "    probs = {\n",
        "        unscramble(bitstring): count / NUM_SHOTS\n",
        "        for bitstring, count in counts.items()\n",
        "    }\n",
        "\n",
        "    # Retrieve the most probable bitstring.\n",
        "    most_probable = max(probs, key=lambda x: probs[x])\n",
        "\n",
        "    print(f\"Expected hidden shift string: {hidden_shift_string}\")\n",
        "    if most_probable == hidden_shift_string:\n",
        "        print(\"Most probable bitstring matches hidden shift 😊.\")\n",
        "    else:\n",
        "        print(\"Most probable bitstring didn't match hidden shift ☹️.\")\n",
        "    print(\"Top 10 bitstrings and their probabilities:\")\n",
        "    display(\n",
        "        {\n",
        "            k: (v, hamming_distance(hidden_shift_string, k))\n",
        "            for k, v in sorted(\n",
        "                probs.items(), key=lambda x: x[1], reverse=True\n",
        "            )[:10]\n",
        "        }\n",
        "    )\n",
        "\n",
        "    return probs, most_probable"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 11,
      "id": "9d3dd975-b628-48fc-8565-86b0ed88c973",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Expected hidden shift string: 011010\n",
            "Most probable bitstring matches hidden shift 😊.\n",
            "Top 10 bitstrings and their probabilities:\n"
          ]
        },
        {
          "data": {
            "text/plain": [
              "{'011010': (0.9743, 6),\n",
              " '001010': (0.00812, 5),\n",
              " '010010': (0.0063, 5),\n",
              " '011000': (0.00554, 5),\n",
              " '011011': (0.00492, 5),\n",
              " '011110': (0.00044, 5),\n",
              " '001000': (0.00012, 4),\n",
              " '010000': (8e-05, 4),\n",
              " '001011': (6e-05, 4),\n",
              " '000010': (6e-05, 4)}"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "probs, most_probable = find_hidden_shift_bitstring(\n",
        "    counts, hidden_shift_string\n",
        ")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "67a88958-c900-4bd1-8329-d530ba3427ea",
      "metadata": {},
      "source": [
        "Let's record the probability of the most probable bitstring before applying readout error mitigation with M3.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 12,
      "id": "a95daf38-009e-44ea-a9c5-8041bf7b22e2",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "0.9743"
            ]
          },
          "execution_count": 12,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "max_probability_before_M3 = probs[most_probable]\n",
        "max_probability_before_M3"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "9531e7bd-e134-4110-ba00-d5a08fe612d9",
      "metadata": {},
      "source": [
        "Now we apply the readout correction learned by M3 to the counts.\n",
        "The function `apply_corrections` returns a quasi-probability distribution. This is a list of `float` objects that sum to $1$. But some values might be negative.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 13,
      "id": "3a24d85c-a980-48cf-96ab-5237d370d7af",
      "metadata": {},
      "outputs": [],
      "source": [
        "def perform_mitigation(mit, counts, qubit_mapping):\n",
        "    # mitigate readout error\n",
        "    quasis = mit.apply_correction(counts, qubit_mapping)\n",
        "\n",
        "    # print results\n",
        "    most_probable_after_m3 = unscramble(max(quasis, key=lambda x: quasis[x]))\n",
        "\n",
        "    is_hidden_shift_identified = most_probable_after_m3 == hidden_shift_string\n",
        "    if is_hidden_shift_identified:\n",
        "        print(\"Most probable bitstring matches hidden shift 😊.\")\n",
        "    else:\n",
        "        print(\"Most probable bitstring didn't match hidden shift ☹️.\")\n",
        "    print(\"Top 10 bitstrings and their quasi-probabilities:\")\n",
        "    topten = {\n",
        "        unscramble(k): f\"{v:.2e}\"\n",
        "        for k, v in sorted(quasis.items(), key=lambda x: x[1], reverse=True)[\n",
        "            :10\n",
        "        ]\n",
        "    }\n",
        "    max_probability_after_M3 = float(topten[most_probable_after_m3])\n",
        "    display(topten)\n",
        "\n",
        "    return max_probability_after_M3, is_hidden_shift_identified"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 14,
      "id": "d304c727-f163-4842-a96a-dad8a067c8d1",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Expected hidden shift string: 011010\n",
            "Most probable bitstring matches hidden shift 😊.\n",
            "Top 10 bitstrings and their quasi-probabilities:\n"
          ]
        },
        {
          "data": {
            "text/plain": [
              "{'011010': '1.01e+00',\n",
              " '001010': '8.75e-04',\n",
              " '001000': '7.38e-05',\n",
              " '010000': '4.51e-05',\n",
              " '111000': '2.18e-05',\n",
              " '001011': '1.74e-05',\n",
              " '000010': '6.42e-06',\n",
              " '011001': '-7.18e-06',\n",
              " '011000': '-4.53e-04',\n",
              " '010010': '-1.28e-03'}"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "print(f\"Expected hidden shift string: {hidden_shift_string}\")\n",
        "max_probability_after_M3, is_hidden_shift_identified = perform_mitigation(\n",
        "    mit, counts, qubit_mapping\n",
        ")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "b62fac82-e883-4030-9793-86af8299abaa",
      "metadata": {},
      "source": [
        "#### Compare identifying the hidden shift string before and after applying M3 correction\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 15,
      "id": "a3864692-582e-4fe0-8059-f0be59cfd229",
      "metadata": {},
      "outputs": [],
      "source": [
        "def compare_before_and_after_M3(\n",
        "    max_probability_before_M3,\n",
        "    max_probability_after_M3,\n",
        "    is_hidden_shift_identified,\n",
        "):\n",
        "    is_probability_improved = (\n",
        "        max_probability_after_M3 > max_probability_before_M3\n",
        "    )\n",
        "    print(f\"Most probable probability before M3: {max_probability_before_M3}\")\n",
        "    print(f\"Most probable probability after M3: {max_probability_after_M3}\")\n",
        "    if is_hidden_shift_identified and is_probability_improved:\n",
        "        print(\"Readout error mitigation effective! 😊\")\n",
        "    else:\n",
        "        print(\"Readout error mitigation not effective. ☹️\")"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 16,
      "id": "98f92ec2-50b5-4456-b40a-3f619cc42ff1",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Most probable probability before M3: 0.9743\n",
            "Most probable probability after M3: 1.01\n",
            "Readout error mitigation effective! 😊\n"
          ]
        }
      ],
      "source": [
        "compare_before_and_after_M3(\n",
        "    max_probability_before_M3,\n",
        "    max_probability_after_M3,\n",
        "    is_hidden_shift_identified,\n",
        ")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "90c3bec4-9a3d-43b7-8160-bcb615230d01",
      "metadata": {},
      "source": [
        "### Plot how CPU time required by M3 scales with shots\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "33addc38-f738-48ed-a29d-9790f446c036",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Applying M3 correction to 5000 shots...\n",
            "\tDone in 0.003321983851492405 seconds.\n",
            "Applying M3 correction to 7500 shots...\n",
            "\tDone in 0.004425413906574249 seconds.\n",
            "Applying M3 correction to 10000 shots...\n",
            "\tDone in 0.006366567220538855 seconds.\n",
            "Applying M3 correction to 12500 shots...\n",
            "\tDone in 0.0071477219462394714 seconds.\n",
            "Applying M3 correction to 15000 shots...\n",
            "\tDone in 0.00860048783943057 seconds.\n",
            "Applying M3 correction to 17500 shots...\n",
            "\tDone in 0.010026784148067236 seconds.\n",
            "Applying M3 correction to 20000 shots...\n",
            "\tDone in 0.011459112167358398 seconds.\n",
            "Applying M3 correction to 22500 shots...\n",
            "\tDone in 0.012727141845971346 seconds.\n",
            "Applying M3 correction to 25000 shots...\n",
            "\tDone in 0.01406092382967472 seconds.\n",
            "Applying M3 correction to 27500 shots...\n",
            "\tDone in 0.01546052098274231 seconds.\n",
            "Applying M3 correction to 30000 shots...\n",
            "\tDone in 0.016769016161561012 seconds.\n",
            "Applying M3 correction to 32500 shots...\n",
            "\tDone in 0.019537431187927723 seconds.\n",
            "Applying M3 correction to 35000 shots...\n",
            "\tDone in 0.019739801064133644 seconds.\n",
            "Applying M3 correction to 37500 shots...\n",
            "\tDone in 0.021093040239065886 seconds.\n",
            "Applying M3 correction to 40000 shots...\n",
            "\tDone in 0.022840639110654593 seconds.\n",
            "Applying M3 correction to 42500 shots...\n",
            "\tDone in 0.023974396288394928 seconds.\n",
            "Applying M3 correction to 45000 shots...\n",
            "\tDone in 0.026412792038172483 seconds.\n",
            "Applying M3 correction to 47500 shots...\n",
            "\tDone in 0.026364430785179138 seconds.\n",
            "Applying M3 correction to 50000 shots...\n",
            "\tDone in 0.02820305060595274 seconds.\n"
          ]
        },
        {
          "data": {
            "text/plain": [
              "Text(0.5, 1.0, 'Time to apply M3 correction')"
            ]
          },
          "execution_count": 17,
          "metadata": {},
          "output_type": "execute_result"
        },
        {
          "data": {
            "text/plain": [
              "<Image src=\"/docs/images/tutorials/readout-error-mitigation-sampler/extracted-outputs/33addc38-f738-48ed-a29d-9790f446c036-2.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "# Collect samples for numbers of shots varying from 5000 to 25000.\n",
        "shots_range = range(5000, NUM_SHOTS + 1, 2500)\n",
        "times = []\n",
        "for shots in shots_range:\n",
        "    print(f\"Applying M3 correction to {shots} shots...\")\n",
        "    t0 = timeit.default_timer()\n",
        "    _ = mit.apply_correction(\n",
        "        pub_result.data.meas.slice_shots(range(shots)).get_counts(),\n",
        "        qubit_mapping,\n",
        "    )\n",
        "    t1 = timeit.default_timer()\n",
        "    print(f\"\\tDone in {t1 - t0} seconds.\")\n",
        "    times.append(t1 - t0)\n",
        "\n",
        "fig, ax = plt.subplots()\n",
        "ax.plot(shots_range, times, \"o--\")\n",
        "ax.set_xlabel(\"Shots\")\n",
        "ax.set_ylabel(\"Time (s)\")\n",
        "ax.set_title(\"Time to apply M3 correction\")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "bdfbeddd-e90e-4a5c-ad19-810603c4f695",
      "metadata": {},
      "source": [
        "#### Interpreting the plot\n",
        "\n",
        "The plot above shows that the time required to apply M3 correction scales linearly in the number of shots.\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "4988be5c-46ad-473d-b1b1-c8280d8f6324",
      "metadata": {},
      "source": [
        "## Scaling up\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 18,
      "id": "1bf46124-9261-472e-864f-ee2ed0209079",
      "metadata": {
        "jp-MarkdownHeadingCollapsed": true
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Hidden shift string 00000010100110101011101110010001010000110011101001101010101001111001100110000111\n"
          ]
        }
      ],
      "source": [
        "n_qubits = 80\n",
        "rng = Random(12345)\n",
        "circuit, hidden_shift, hidden_shift_string = run_hidden_shift_circuit(\n",
        "    n_qubits, rng\n",
        ")\n",
        "\n",
        "print(f\"Hidden shift string {hidden_shift_string}\")"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 19,
      "id": "53cd9086-27f8-481b-ba96-9bac9b146188",
      "metadata": {},
      "outputs": [],
      "source": [
        "isa_circuit = get_isa_circuit(circuit, backend)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 20,
      "id": "b176b2b3-de30-4a9b-905d-a06c7365376a",
      "metadata": {},
      "outputs": [],
      "source": [
        "job = run_sampler(backend, isa_circuit, NUM_SHOTS)\n",
        "mit, qubit_mapping = setup_mthree_mitigation(isa_circuit, backend)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 21,
      "id": "c4dda06e-e2c6-46b9-9e74-d90bb318419b",
      "metadata": {},
      "outputs": [],
      "source": [
        "counts, pub_result = get_bitstring_counts(job)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 22,
      "id": "3421d2f3-083f-44a9-bee4-5f05e43ebb4f",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Expected hidden shift string: 00000010100110101011101110010001010000110011101001101010101001111001100110000111\n",
            "Most probable bitstring matches hidden shift 😊.\n",
            "Top 10 bitstrings and their probabilities:\n"
          ]
        },
        {
          "data": {
            "text/plain": [
              "{'00000010100110101011101110010001010000110011101001101010101001111001100110000111': (0.50402,\n",
              "  80),\n",
              " '00000010100110101011101110010001010000110011100001101010101001111001100110000111': (0.0396,\n",
              "  79),\n",
              " '00000010100110101011101110010001010000110011101001101010101001111001100100000111': (0.0323,\n",
              "  79),\n",
              " '00000010100110101011101110010001010000110011101001101010101001101001100110000111': (0.01936,\n",
              "  79),\n",
              " '00000010100110101011101110010011010000110011101001101010101001111001100110000111': (0.01432,\n",
              "  79),\n",
              " '00000010100110101011101110010001010000110011101001101010101001011001100110000111': (0.0101,\n",
              "  79),\n",
              " '00000010100110101011101110010001010000110011101001101010101001110001100110000111': (0.00924,\n",
              "  79),\n",
              " '00000010100110101011101110010001010000010011101001101010101001111001100110000111': (0.00908,\n",
              "  79),\n",
              " '00000010100110101011100110010001010000110011101001101010101001111001100110000111': (0.00888,\n",
              "  79),\n",
              " '00000010100110101011101110010001010000110011101001100010101001111001100110000111': (0.0082,\n",
              "  79)}"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "probs, most_probable = find_hidden_shift_bitstring(\n",
        "    counts, hidden_shift_string\n",
        ")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "2ec8335f-d3ef-4951-8bbd-b4d01f899be1",
      "metadata": {},
      "source": [
        "We see that the correct hidden shift string is found. Furthermore, the nine next-most-probable bitstrings are wrong in only one position.\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "98d1fb7b-e696-4cef-a8b8-18e44b96eece",
      "metadata": {},
      "source": [
        "Record the most likely probability:\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 23,
      "id": "c5a9dbdc-ebb7-470a-ac31-10b41217eeed",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "0.50402"
            ]
          },
          "execution_count": 23,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "max_probability_before_M3 = probs[most_probable]\n",
        "max_probability_before_M3"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 24,
      "id": "b4514b88-3f2d-4e96-88ce-6d291576568d",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Expected hidden shift string: 00000010100110101011101110010001010000110011101001101010101001111001100110000111\n",
            "Most probable bitstring matches hidden shift 😊.\n",
            "Top 10 bitstrings and their quasi-probabilities:\n"
          ]
        },
        {
          "data": {
            "text/plain": [
              "{'00000010100110101011101110010001010000110011101001101010101001111001100110000111': '9.85e-01',\n",
              " '00000010100110101011101110010001010000110011100001101010101001111001100110000111': '6.84e-03',\n",
              " '00000010100110101011100110010001010000110011101001101010101001111001100110000111': '3.87e-03',\n",
              " '00000010100110101011101110010011010000110011101001101010101001111001100110000111': '3.42e-03',\n",
              " '00000010100110101011101110010001010000110011101001101010101001111001100100000111': '3.30e-03',\n",
              " '00000010100110101011101110010001010000110011101001101010101001110001100110000111': '3.28e-03',\n",
              " '00000010100010101011101110010001010000110011101001101010101001111001100110000111': '2.62e-03',\n",
              " '00000010100110101011101110010001010000110011101001101010101001101001100110000111': '2.43e-03',\n",
              " '00000010100110101011101110010000010000110011101001101010101001111001100110000111': '1.73e-03',\n",
              " '00000010100110101011101110010001010000110011101001101010101001111001000110000111': '1.63e-03'}"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "print(f\"Expected hidden shift string: {hidden_shift_string}\")\n",
        "max_probability_after_M3, is_hidden_shift_identified = perform_mitigation(\n",
        "    mit, counts, qubit_mapping\n",
        ")"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 24,
      "id": "ccd14379-6c34-4be8-93ab-728356bf81e0",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Most probable probability before M3: 0.54348\n",
            "Most probable probability after M3: 0.99\n",
            "Readout error mitigation effective! 😊\n"
          ]
        }
      ],
      "source": [
        "compare_before_and_after_M3(\n",
        "    max_probability_before_M3,\n",
        "    max_probability_after_M3,\n",
        "    is_hidden_shift_identified,\n",
        ")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "96b6648b-5c26-4a0c-92f2-2d1810fa14b8",
      "metadata": {},
      "source": [
        "The results show that readout error was the dominant source of error and the M3 mitigation was effective.\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "id": "a1b8767d",
      "source": "© IBM Corp., 2017-2026"
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3"
    },
    "toc": {
      "base_numbering": 0
    }
  },
  "nbformat": 4,
  "nbformat_minor": 4
}