diff --git a/docs/tutorials/_toc.json b/docs/tutorials/_toc.json index 9fbb6f4eae1..4e1bc219575 100644 --- a/docs/tutorials/_toc.json +++ b/docs/tutorials/_toc.json @@ -178,6 +178,10 @@ { "title": "Circuit cutting for depth reduction", "url": "/docs/tutorials/depth-reduction-with-circuit-cutting" + }, + { + "title": "Readout error mitigation for the Sampler primitive using M3", + "url": "/docs/tutorials/readout-error-mitigation-sampler" } ] }, diff --git a/docs/tutorials/index.mdx b/docs/tutorials/index.mdx index ba9110988ca..871a57a65d8 100644 --- a/docs/tutorials/index.mdx +++ b/docs/tutorials/index.mdx @@ -135,6 +135,8 @@ Addons enable advanced circuit manipulation, such as cutting, backpropagating ob * [Circuit cutting for depth reduction](/docs/tutorials/depth-reduction-with-circuit-cutting) +* [Readout error mitigation for the Sampler primitive using M3](/docs/tutorials/readout-error-mitigation-sampler) + diff --git a/docs/tutorials/readout-error-mitigation-sampler.ipynb b/docs/tutorials/readout-error-mitigation-sampler.ipynb new file mode 100644 index 00000000000..21acee86a1b --- /dev/null +++ b/docs/tutorials/readout-error-mitigation-sampler.ipynb @@ -0,0 +1,1287 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "source": [ + "{/* cspell:ignore braket, Zgrm, newcommand, probs, quasis, topten */}\n", + "# Readout error mitigation for the Sampler primitive using M3\n", + "$$\n", + "\\newcommand{\\Mm}{M}\n", + "\\newcommand{\\pideal}{\\vec{p}}\n", + "\\newcommand{\\pnoisy}{\\tilde{p}}\n", + "$$\n", + "\n", + "*Estimated QPU usage: 1 minute (tested on IBM `ibm_kingston`)*\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.\n", + "However, readout error mitigation is a highly effective method that is also applicable to the Sampler primitive.\n", + "This tutorial explains how to use the [M3 Qiskit addon](https://qiskit.github.io/qiskit-addon-mthree/) to mitigate readout error for the Sampler primitive.\n", + "\n", + "The M3 Qiskit addon implements an efficient method for readout error mitigation.\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 readout error is more important than error due to to incorrectly collapsing the state.\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 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.\n", + "In the following we be concerned only with readout error.\n", + "\n", + "### 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 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 $\\Mm_{i,j}$\n", + "that $i$ is read, given that the true value is $j$.\n", + "That is,\n", + "$$\n", + " \\Mm_{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", + "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 $\\Mm$ the _assignment matrix_.\n", + "For fixed true value $j$, summing the probability over all noisy outcomes $i$ must give one. That is\n", + "$$\n", + " \\sum_{i=0}^{2^n - 1} \\Mm_{i,j} = 1 \\text{ for all } j\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 $\\Mm_{i,j}$ by\n", + "repeatedly preparing each basis state $|j \\rangle$ and then computing the frequencies\n", + "of 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 $\\Mm$ 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 $\\pnoisy \\in \\mathbb{R}^{2^n}$.\n", + "The (estimated) probability ${\\pnoisy}_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", + " \\pnoisy = \\Mm \\pideal, \\tag{2}\n", + "$$\n", + "where $\\pideal$ is the true distribution. In words, the readout error has the effect of multiplying\n", + "the ideal distribution over bitstrings $\\pideal$ by the assignment matrix $\\Mm$ to\n", + "produce the observed distribution $\\pnoisy$.\n", + "We have measured $\\pnoisy$ and $\\Mm$, but have no direct access to $\\pideal$. In principle, we will\n", + "obtain the true distribution of bitstrings for our circuit\n", + "by solving equation (2) for $\\pideal$ 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 $\\Mm$. Linear algebra\n", + " routines in software libraries employ methods that are more stable, accurate, and efficient.\n", + "- When estimating $\\Mm$, 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, $\\Mm$ really represents\n", + " only readout error. But when we _use_ $\\Mm$ 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 $\\Mm$ grow exponentially in $n$:\n", + "- The estimation of $\\Mm$ and $\\pnoisy$ 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 $\\Mm$.)\n", + " However, if no assumptions are made on the bitstrings observed\n", + " when performing mitigation, the number of shots required to estimate $\\Mm$ grows\n", + " at least exponentially in $n$.\n", + "- $\\Mm$ is a $2^n \\times 2^n$ matrix.\n", + " When $n>10$, the amount of memory required to store $\\Mm$ is\n", + " greater than the memory available in a powerful laptop.\n", + "\n", + "Further limitations are:\n", + "\n", + "- The recovered distribution $\\pideal$ may have one\n", + " or more negative probabilities (while still summing to one). One solution\n", + " is to minimize $||\\Mm \\pideal - \\pnoisy||^2$ subject to the constraint that\n", + " each entry in $\\pideal$ 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", + "### Qiskit M3 addon: Scaling to longer bitstrings\n", + "\n", + "The M3 method is an alternative to direct solution of equation (2)\n", + "for mitigating readout error in the experimental distribution over bitstrings.\n", + "While the direct method is limited to bitstrings no longer than about $10$ bits,\n", + "M3 can handle much longer bitstrings.\n", + "The two important features of M3 are\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 $\\Mm$ explicitly, we use a much smaller effective matrix that records\n", + " probabilities only for bitstrings collected when constructing $\\pnoisy$.\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 $\\Mm$.\n", + "Then, we repeatedly run the circuit of interest and collect bitstrings which we use to construct\n", + "both $\\pnoisy$ and, with the help of the building blocks, an effective $\\Mm$.\n", + "\n", + "More precisely,\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 $\\pnoisy$,\n", + " we construct an effective approximation to $\\Mm$ using only\n", + " bitstrings that are sampled when constructing $\\pnoisy$. This effective matrix\n", + " is built using the one- 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 $\\pnoisy$, which is much smaller than\n", + " the dimension $2^n$ of the full assignment matrix $\\Mm$ .\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", + "### Application of M3 to a quantum algorithm\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 QPU's can be found in [this paper](https://journals.aps.org/prx/accepted/a9074K06A8e1590147da9c69f8c4b64c28247be5a) ([arXiv version](https://arxiv.org/abs/2401.07934)).\n", + "\n", + "$$\n", + "% Define CZ as a math operator\n", + "%\\DeclareMathOperator{\\CZ}{CZ}\n", + "\\newcommand{\\CZ}{CZ}\n", + "\\newcommand{\\bra}[1]{\\langle #1|}\n", + "\\newcommand{\\ket}[1]{|#1\\rangle}\n", + "\\newcommand{\\braket}[2]{\\langle #1|#2\\rangle}\n", + "\\newcommand{\\ketbra}[2]{|#1\\rangle\\langle#2|}\n", + "\\newcommand{\\Zgr}{\\mathbb{Z}_2^n}\n", + "\\newcommand{\\Zgrm}{\\mathbb{Z}_2^m}\n", + "$$\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: \\Zgr \\rightarrow \\mathbb{Z}_2$ is defined\n", + "by $x \\cdot y = \\sum_i x_i y_i$.\n", + "\n", + "#### 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", + "H^{\\otimes n} = \\frac{1}{\\sqrt{2^n}} \\sum_{x,y \\in \\Zgr} (-1)^{x \\cdot y} \\ket{y}\\bra{x}.\n", + "$$\n", + "Consider a state $\\ket{s}$ corresponding to fixed bitstring $s$.\n", + "Applying $H^{\\otimes n}$, and using $\\braket{x}{s} = \\delta_{x,s}$,\n", + "we see that the Fourier transform of $\\ket{s}$ can be written\n", + "$$\n", + " H^{\\otimes n} \\ket{s} = \\frac{1}{\\sqrt{2^n}} \\sum_{y \\in \\Zgr} (-1)^{s \\cdot y} \\ket{y}.\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", + " \\ket{s} = H^{\\otimes n} H^{\\otimes n} \\ket{s} = H^{\\otimes n} \\frac{1}{\\sqrt{2^n}} \\sum_{y \\in \\Zgr} (-1)^{s \\cdot y} \\ket{y}.\n", + "$$\n", + "\n", + "#### The hidden shift problem\n", + "\n", + "$$\n", + "\\newcommand{\\ffunc}{f}\n", + "$$\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 \\Zgrm$ be bitstrings of length $m$.\n", + "We define $\\ffunc: \\Zgrm \\times \\Zgrm \\rightarrow \\{-1,1\\}$ by\n", + "$$\n", + " \\ffunc(x, y) = (-1)^{x \\cdot y}.\n", + "$$\n", + " Let $a,b \\in \\Zgrm$ be fixed bitstrings of length $m$.\n", + " We furthermore define $g: \\Zgrm \\times \\Zgrm \\rightarrow \\{-1,1\\}$ by\n", + "$$\n", + " g(x, y) = \\ffunc(x+a, y+b) = (-1)^{(x+a) \\cdot (y+b)},\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$\n", + " with all pairs of strings that combined have 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", + "#### Quantum circuits for $f$ and $g$\n", + "\n", + "We will need the following ingredients in order to implement $f$ and $g$ as quantum circuits.\n", + "\n", + "For single qubit classical states $\\ket{x_1}, \\ket{y_1}$, with $x_1,y_1 \\in \\mathbb{Z}_2$,\n", + "the controlled $Z$ gate $\\CZ$ may be written\n", + "\\begin{equation}\n", + "\\CZ \\ket{x_1}\\ket{y_1} = (-1)^{x_1 y_1} \\ket{x_1}\\ket{y_1}.\n", + "\\end{equation}\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 $\\ffunc = \\ffunc(x,y)$:\n", + "$$\n", + "%\\CZ_{x,y} \\ket{z} =\n", + "U_f \\ket{x}\\ket{y} = \\CZ_{x,y} \\ket{x}\\ket{y} = (-1)^{x \\cdot y} \\ket{x}\\ket{y}.\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", + " X_a X_b \\ket{x}\\ket{y} = \\ket{x+a}\\ket{y+b}.\n", + "$$\n", + "\n", + "The second black box $g$ is implemented by the unitary $U_g$, given by\n", + "$$\n", + "%U_g \\ket{x}\\ket{y} = X_aX_b \\CZ_{x,y} X_aX_b \\ket{x}\\ket{y}.\n", + "U_g = X_aX_b \\CZ_{x,y} X_aX_b.\n", + "$$\n", + "To see this, we apply the operators from right to left to the state $\\ket{x}\\ket{y}$.\n", + "First\n", + "\\begin{equation}\n", + " X_a X_b \\ket{x}\\ket{y} = \\ket{x+a}\\ket{y+b}.\n", + "\\end{equation}\n", + "Then,\n", + "\\begin{equation}\n", + " \\CZ_{x,y} \\ket{x+a}\\ket{y+b} = (-1)^{(x+a)\\cdot (y+b)} \\ket{x+a}\\ket{y+b}\n", + "\\end{equation}\n", + "Finally,\n", + "\\begin{equation}\n", + " X^a X^b (-1)^{(x+a)\\cdot (y+b)} \\ket{x+a}\\ket{y+b} = (-1)^{(x+a)\\cdot (y+b)} \\ket{x}\\ket{y},\n", + "\\end{equation}\n", + "which is indeed the quantum version of $f(x+a, y+b)$.\n", + "\n", + "#### The hidden shift algorithm\n", + "\n", + "$$\n", + "\\newcommand{\\kzero}{\\ket{0}^{\\otimes m}}\n", + "$$\n", + "\n", + "Now we put the pieces constructed above together to solve the hidden shift problem.\n", + "We begin by applying Hadamards to the registers initialized to the all-zero state.\n", + "$$\n", + "H^{\\otimes 2m} = H^{\\otimes m} \\otimes H^{\\otimes m} \\kzero\\kzero = \\frac{1}{\\sqrt{2^{2m}}} \\sum_{x, y \\in \\Zgrm} (-1)^{x \\cdot y} \\ket{x}\\ket{y}.\n", + "$$\n", + " Next, we query the oracle $g$ to arrive at\n", + " $$\n", + " U_g H^{\\otimes 2m} \\kzero\\kzero\n", + " = \\frac{1}{\\sqrt{2^{2m}}} \\sum_{x, y \\in \\Zgrm} (-1)^{(x+a) \\cdot (y+b)} \\ket{x}\\ket{y}\n", + " $$\n", + " $$\n", + " \\approx \\frac{1}{\\sqrt{2^{2m}}} \\sum_{x, y \\in \\Zgrm} (-1)^{x \\cdot y + x \\cdot b + y \\cdot a} \\ket{x}\\ket{y}.\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", + " U_f U_g H^{\\otimes 2m} \\kzero\\kzero\n", + " \\approx \\frac{1}{\\sqrt{2^{2m}}} \\sum_{x, y \\in \\Zgrm} (-1)^{x \\cdot b + y \\cdot a} \\ket{x}\\ket{y}.\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", + " H^{\\otimes 2m} U_f U_g H^{\\otimes 2m} \\kzero\\kzero\n", + " \\approx \\ket{b}\\ket{a}.\n", + " $$\n", + " The circuit is finished. In the absence of noise, sampling the quantum registers will\n", + " return the bitstrings $b, a$ with probability one.\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", + " H^{\\otimes n} U_{\\tilde{f}} H^{\\otimes n} U_g H^{\\otimes n} \\ket{0}^{\\otimes n} = \\ket{s}.\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 $\\ket{b}\\ket{a}$ instead of the desired $\\ket{a}\\ket{b}$.\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Requirements\n", + "\n", + "Before starting this tutorial, ensure that you have the following installed:\n", + "\n", + "- Qiskit SDK 2.1 or later with visualization support (`pip install 'qiskit[visualization]'`)\n", + "- Qiskit Runtime 0.41 or later (`pip install qiskit-ibm-runtime`)\n", + "- M3 Qiskit addon 3.0 (`pip install mthree`)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "jp-MarkdownHeadingCollapsed": true + }, + "source": [ + "## Setup\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Step 1: Map classical inputs to a quantum problem\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "First, we write the functions to implement the hidden shift problem as a `QuantumCircuit`." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "metadata": { + "editable": true, + "scrolled": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "outputs": [], + "source": [ + "from collections.abc import Iterator, Sequence\n", + "from random import Random\n", + "\n", + "from qiskit.circuit import CircuitInstruction, QuantumCircuit, QuantumRegister, Qubit\n", + "from qiskit.circuit.library import CZGate, HGate, XGate\n", + "\n", + "\n", + "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(qubits: Sequence[Qubit], shift: int) -> 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(qubits: Sequence[Qubit], shift: int) -> 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", + "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", + "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", + "metadata": {}, + "source": [ + "We'll start with a fairly small example:" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Hidden shift string 011010\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAPwAAADiCAYAAABqd6qMAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjYsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvq6yFwwAAAAlwSFlzAAAPYQAAD2EBqD+naQAAFoVJREFUeJzt3Xt0VOW5x/HvnpmEEKhJEGhjIXI5xAgKDSg3bYsgiEdUDl6Q0hZEFEtTD9ZFq/XYJaWo7alilVoQjuCFemxB0dKIF5ADQogSIwQwBcEQQO7EmHtmZu/zx6YKJCRDZrLfved9PmtlrUxm5t3vOzy/2XtmD/MYlmVZCCG04FM9ASGEcyTwQmhEAi+ERiTwQmhEAi+ERiTwQmhEAi+ERiTwQmhEAi+ERiTwQmhEAi+ERiTwQmhEAi+ERiTwQmhEAi+ERiTwQmhEAi+ERiTwQmhEAi+ERiTwQmhEAi+ERiTwQmhEAi+ERgKqJxArVtgEy4ztoIYPwy/PiSIyXqjBuAi8FTYJTZ4G5V/GduCU8wgsWSChF83ySg3GRyVbZuwfaLDHjPUztohPHqnB+Ai8ECIiEnghNCKBF0IjEnghNKJN4Euqqxidt/a0v128OlfZfIR+3FCDygK/e/durrzySjIzM8nOzmbz5s2qpqKFzcfgjvXw/Vy4aTX8vRQsS/WsnFNWB49sgRFv2j+PbYUv6lTPynnKAn/33XczadIkdu7cye9//3smTpyIpVMFOmjjEbhnE2wpg6oQ7K2yC37BP1XPzBlVIfjROnijFMqD9s/re2HSeqgJqZ6ds5R88Obo0aNs2rSJ3Fz7cGbkyJFYlkVBQQGXXXZZq233o/Iyrt74XquN70aWBY8XQf0Zp3LrTHjhU5jYE76RoGp2zlhZau/NQ6fsT4IWHK+D3P1wUzfn5qK6BpUEvrS0lPT0dBISvq60bt26UVpaetbAV1dXU1xc3PiA4TCXRrDd/ilprBoy7KvLkbx+KiwsBL8/gtHdKWgZ7K3KbvQ6nxXm7x/uJiup0vF5OemtYz2oNVMb/L02DG/tKqP7ic+i34iiGszKyiI5OTniaXrmo7XFxcUMGDCg0esChkH1mFtaZbuDBg0i5OWXGoZB9l+r8CW2bXBVdV09d/zgJmr3blMyNad0vfvPdBp5B0bg9EMZKxRk1bKlLHz2Z1FvQ1UNFhQU0L9//4jHUhL4jIwMDh48SDAY/GovX1JSQkZGxlnvk5WVRUFBQeNXhsMw+w+tMtf8/HxP7+EBXiqrZGNVG0KnvWVj0SUZFi5/HsNQODkHlNa35bEjfs58uZ4Q8DPv9hF0mXaWujoXimowKyvrnMYyLEXvlI0YMYLbbruNO++8k3feeYfp06ezc+dOjBZUnxUKERo3sVXmGXh1KUbAMwdCjaoNwYx82FpmX643ITURnv8ufLud6tk5Y1mJ/V6GYdjr9xvwQF8Ye2FsxvdKDSp7l37+/PksXryYzMxMZs6cydKlS1sUdtG8pADMvwKe+y5c18X+27299Qk7wM3dYNU1cENX+/IvL4ld2L1E2a6rV69ebNy4UdXmtZSVAv3Ph9dKwafhc2tKIvTrAMv2QltvH7S1mDaftBNCSOCF0IoEXgiNSOCF0IgEXgiNxEfgDR+knBf7cVPOs8cWojkeqcG4ODlh+H0Elixw/VcEi/jllRqMi8Bz8gGPlwMW4U1eqEF3z04IEVMSeCE0IoEXQiMSeCE0IoEXQiMSeCE0IoEXQiNxcx5edW/ucNjEjPHmfT7we+SDP7qvHxfUYCTiIvCqe3OHwybpw1/maFltTDffKS2Jg2smuL7odV8/LqjBSLn/kYyE4t7cpknMix3sMWO912wNuq8f1NdgpOIj8EKIiEjghdCIBF4IjUjghdCINoFX3Zs7s1sKeS9dj99vfz/07JwB3DOxj2PbV0339eOCGkRl4GfPnk1mZiY+n48VK1aomoZjdpaUs2rDfmb8sA+9e6Zy1eXpzHt5h2PbP1QNT2yDZ0+2iD5Q7dimwQXrD5p2u+jndtmXC49DyCtnAGJI2Xn4kSNHMnHiRKZMmaJqCo57dNEWNr54PbeO6sG02RswTWe6fO0qhzvet4s+eHKTi3ZC13ZwTRdHpgAK1x8yIScPtn9hd4wF+Ps+2F8NTw2GgDbHuQoDP3jwYMe3qbo3d33QZNPWIwzp15mPi487tt3fFUFNGE6NV8iCOVvhqnRIdKhXpqr1rzl4etg52R9+axmsOwTDL3BsKspr0DOftHNzf/hgKLI91aW90uib2YHCT44zfnQPXlm1p9n7fFRYSEKg5X2hQpbBxyca7w8fDod57YNP6dWmqsXj4/L1Ayw73p3acFqDv9eG4a/bT5B6qCSq8UH6w8ecu/vD++HSBU3ewjDgmQeHMn3ORvYfrmLtc9eRu34fFVXBJu83eNAgINzkbZrk89N/eR2Gr2HBVNfUMGXyj6ne+UHLxwd3rx/odu+LdPjeDzB8px+7W2aYVSvfYMFTt0c1PtIfPvbc3B8+GLIYPPXjJm8zffzF5G09QtEuu2fz4y8U8cg9l/GzR/OavN+m/Pyo93DzjlWwrfY8zDPeo01LbsP//eXPUTeWdPv6t9d+g2eOwZlPLYl+g0fGDyVrkvSHd8ywYcOYMWMGY8eObfEYqntzB4MmiQMWt8r26wtuJyEhuneVjtTApPVQGbRfy7fx2XvceYPhO+dHP0e3r9+yYM4WeHO/3RseINEHY7rC/X3txyJaqmsw4rFiMkoLPPzwwyxatIijR4+ybds2cnJy2LRpE126OPi2sSY6t4XXhsPbn8P2Mrsv/Jiu0KGN6pk5wzDgwX5wQwa8+zkYwMhvQ5/U2ITdS5QG/uGHH1a1ee0kBeyCvyFD9UzUMAzo28H+0ZlGZyCFEBJ4ITQigRdCIxJ4ITQigRdCI/EReMW9uX0++wsXY61TWhI+D/wL6b5+UF+DkVL+wZtYUf0Vwbp/TbPu68cFNRjRcPESeCFE87zz9CmEiJoEXgiNSOCF0IgEXgiNSOCF0IgEXgiNSOCF0IgEXgiNeOY77ZrjhU85ifjmhRqMi8BbYZPQ5Gmx78+dch6BJQsk9KJZXqnB+Khky4z9Aw32mLF+xhbxySM1GB+BF0JERAIvhEYk8EJoRAIvhEaUBL6srIwxY8aQmZlJv379GDVqFJ9++mmrbrOkuorReWtP+1skjfziRdCEtw7AUzvgb59BRdMt3eLSP8thfjEsKLZbaDvNDTWo5LScYRjMmDGDq6++GoCnnnqKqVOnsnbt2mbvK87d8VqY8j6cqLNbTSX54ekd8MxQuKRhU9W4Y1nw+HZ4tcRuk20Az38Kt3aHGX1Uz85ZSvbwqampX4UdYOjQoZSUxKBlr2jU74rgcI0ddk62Sa4Ow8wPwdTg+442H4PXSuy+cqYFYcv+/W+fwUfOtal3BVd88ObJJ5/kxhtvbPI2segP/1F5GVdvfO+c5hZJf3g3C1uw9mA2Jg2bqJXXhVmRv4tuidVK5uaUJccvpM7sAGc8BnWmxZLC49ChNPqNKKpBz/WHnzVrFnv27OHZZ59t8nax6A/fPyWNVUOGfXU5ktdPkfWHdy8jkED/5fWNXldTVc2kB6ZRuX2d4/NyUo/7l5E25KZGrjFY+da7PP2HCVFvQ1UNeqo//G9/+1tyc3N5++23m32WcnN/eLd79Egln9W3a7CHS26XzLvPP0mC4d0ntEjkV6fxYlmYeuv0f8dEI8z9Nw7ksgn69IdXFvhZs2Z9FfaUlJRmb5+cnHzWZzIrFCLUCnMEyM7OjllvblVmfwFTN0DItN+0AggY8Mt+fgZlZKueXqu71IQPNsCuL+33LwASDLg4zc+UoT0IxOCdLK/UoJKvqd6+fTuXXHIJPXv2pH379gAEAgE2b97covGsUIjQuIkxnqUt8OpSzwceYF8lvLgbNhyGw7UwtRfcfbHqWTmnLgxvlMJLu+FANVzfFR7oC4kxOnjzSg0qqeQ+ffogX4fvrK7t4Vf94M198FAhXNhe9Yyc1cYPt3SH9gF7/QM7xi7sXiKftBNCIxJ4ITQigRdCIxJ4ITQigRdCIxJ4ITQSH4E3fJByXuzHTTnPHluI5nikBr3/iRLA8PsILFng+q8IFvHLKzUYF4Hn5AMeLwcswpu8UIPunp0QIqYk8EJoRAIvhEYk8EJoRAIvhEYk8EJoRAIvhEbi5jy86t7c4bCJGePN+3zg98gHf3RfPy6owUjEReBV9+YOh03Sh7/M0bLamG6+U1oSB9dMcH3R675+XFCDkXL/IxkJxb25TZOYFzvYY8Z6r9kadF8/qK/BSMVH4IUQEZHAC6ERCbwQGpHAC6ERZYEfP348ffv2JTs7m4EDB7J69epW3Z7q3tyZ3VLIe+l6/H673dPsnAHcM9HZXsUlFVBUZv/udFsAN6y/KgQ7vrB/rws7umlwQQ2i8rTcggULSE1NhZPdMUeMGMGxY8fw+eLzoGNnSTmrNuxnxg/78Ob7+7nq8nS+d/s/HNl2XRh+uRnyj37dXe6pT2BAJ/hmW0emoHT9AP8ohTlbv778WBEkBWB0F8em4ArK0vWvsAOUl5ermoajHl20hQnX9mTxb75HzqN5mA41Z396B3x4FIKm3Rcd4Hgd/Dzfkc1/RdX6d5XbYa8/Zf0hC37zMeypcGQKrqH0gzf33nsvr7/+OuXl5SxfvrzJvbub+8MHQ5EVbn3QZNPWIwzp15mPi49HdJ+PCgtJCDTs7R4p04Lln3+HoHX6Y2sBu780WZlfzAUJ0Z1Dd/P6AZaWdSVkdmzQPTdsmsz/4Ci3pR2Ianx7MOkP36y5c+cyd+5cVq1axS9+8Qs2bNhAYmJio7d1d394P1y6oNmxLu2VRt/MDhR+cpzxo3vwyqo9zd5n8KBBQMtfcBoJbei/rPFA11VXMuHOn1JZtLbR6yPn3vUD9HzwdVIH3tDg72F8LH/nff77dzdHNT7SH/7cjB49mpycHIqKis4aajf3hw+GLAZP/bjJ2xgGPPPgUKbP2cj+w1Wsfe46ctfvo6Iq2OT9NuXnR7WHsyz4r0O1HA0nNbiuTXJ7cp97kna+6ALl5vUDvFvRidfKwwRp2B/+J9dezvBbpT98q6qpqeHQoUN0794dgLy8PI4fP06PHj3Oeh8394cPBk2g6YKfPv5i8rYeoWiX/Tb54y8U8cg9l/GzR/OavF//7GwSEqJ7q+X+Q/DA5q9fvwK08cGt3X18t0+/qMbGA+vPDMJ7a6Cs3n7tDhAwILWNn59ckUG7QEZU4+OCGoyUkv7wJ06c4LrrrqOiooJAIEC7du2YPXs2w4cPb9F4qntzB4MmiQMWt8r26wtuj7rgATYegXk74LNKOL8N/Pjf4JZu9p43Wl5Y/9Fa+83L9w7al0ekQ05v6NjwwKdFVNdgxGPFZJRz1KFDB/Lymn5mF7E1tLP9o6tOSfCbyF/qxq34POkthGiUBF4IjUjghdCIBF4IjUjghdCIBF4IjcRH4BX35vb57C9cjLVOaUl44T8P6r5+UF+DkVLywZvWoPorgnX/mmbd148LajCi4eIl8EKI5nnn6VMIETUJvBAakcALoREJvBAakcALoREJvBAakcALoRFXfKddLKj+0IPuHzzRff24oAYjEReBV92bW/f+6LqvHxfUYKTc/0hGQvrDK6X7+kF9DUYqPgIvhIiIBF4IjUjghdCIBF4IjWgTeNW9ud3QH702DCWVUNF0d6dW4Yb1WxZ8Xg0Hq+3fnaa6BnFD4BcvXoxhGKxYsUL1VFrVqf3Re/dM5arL05n38g5Htm1ZsPCfcPUqmLgWRq6C/yqAmtbqjdQIlesHKDoB49bAuNXwH6vh5jWwvcyxzbuG0vPwJSUlLFy4kMGDB6uchmMeXbSFjS9ez62jejBt9gbH+qP/ZTc8vwtqTzm7s+YgVIfgiUGOTAEUrv9QDUzPg5pTemburYK7N8KrI+yuNLpQFnjTNJk6dSpPP/009913X7O3l/7wLWNZsPBgX2rN0/+p603YcNjknQ+2c34gumN8N68f4LXydILhbzY4oA2GTf6Ud4gbUg5FNT5If/hmPfHEE1xxxRVnbQ99JukP3zJN9Yevr6rgpqn/SUXRuRVgQ+5dP030hw/h48W385j1mPSHb1Xbtm1j+fLlrFu3LuL7SH/4lrEsuO9gkEozocF1ie2+wauL/kiHGOzh3bp+gBXl6bxdYRI6Yw+fgMmPRw3h+lukP3yrWr9+PSUlJfTq1QuAQ4cOcdddd7F//35ycnIavY/0h2+5u3bDnz6BulNewwcMuPKbPq4eGMmBaNPcvv5v19j94UNnHCgE/D5+OuQCOiZdENX4uKAGI+WKb60dNmwYM2bMYOzYsS26v+re3G7vj25ZsHiX/RM0IWRBvzT40xBIikEduX39ANvK4OFCKK0C04JObeDxQdA7NSbTVF6DkVJ+Wk60PsOAKZnw7mjIudj+283dYhN2r7gkDZYNh5/3ti/f0zt2YfcSV/yTr127NoJbiWi18cP5bVTPQq2URNUzUEv28EJoRAIvhEYk8EJoRAIvhEYk8EJoJD4CL/3hldJ9/aC+BiPlitNy0TL8PgJLFij7imC/38fBNRO0/Zpm3dePC2owUnEReE4+4CoPWPx+X3MfuY9ruq8fF9RgJNw9OyFETEnghdCIBF4IjUjghdCIBF4IjUjghdCIBF4IjcTNeXgv9OYW8c0LNRgXgfdKb24Rv7xSg/FRyR7pzS3imEdqMD4CL4SIiAReCI1I4IXQiAReCI1oE3g39OYWenNDDSo7LTds2DD27t1LSkoKAOPGjePXv/61qunEvfcPw7xPYE+FfTnvCIzuYjep0MGRGvjjDljzuX15WQlc3gk6atQqGtXn4efOndvi9lIicmsPwq8K7BbR//LO59AhCWb0UTkzZ1QE4Ufr4It6CJ9srFZUZv/tb1dB+4Z9NuNWXHzwJlIt6c3tdZYFc7efHnaAoAWvfAaTe0FqnHdjeaMUKoNfhx3APPlE8I99ML6Hc3NRXYNKA3///ffz0EMPcdFFFzFnzhwuuuiis962urqa4uLixq8Mh4mkB2pLenMXFhY22y7azYKWwYHq7Eav81thcjfvJiup0vF5OWn1sR7UmQ0bydWGYfWnZfT64rPoN6KoBrOyskhOTo54msoC/8ILL5CRkYFlWSxevJhRo0axZ88e/GdZWHFxMQMGDGj0uoBhUD3mllaZ56BBgwipb7DbcobBd/63An9SuwZXVdfVc8f4G6ndf5Yn0jjR9a6n6XjNNHyB04/dzVA9b/51CQv/5+dRb0NVDRYUFJy1jXpjlAU+IyMDAMMwmDJlCjNnzmTv3r306NH48VVWVhYFBQWNDxYOw+w/tMo88/PzPb2HB3ilrJJ1VW0JnnJSxodFRjuTRa8vVTo3JxwIJjHnsL9B//bEQID5d15D+vSz1NW5UFSDWVlZ5zSWkv7wtbW1VFZW0rFjRwByc3OZPHkyBw4cICHh3N9B8UpvblXqw/BAAWw8Agk+uz96elu7P3zntqpn54zcfTBnC/h9gGW/hn+oH1zTJTbje6UGlVTyl19+ybXXXkt9fT0+n4+0tDRWrlzZorCL5iX64fGBsLcSdn0J32oLfVL1OSUH8O9dYVg6fHgMDODyjtDW28/jLaJkyZ07dz774bloNRe2t390lRyA739L9SzU0uaTdkIICbwQWpHAC6ERCbwQGpHAC6GR+Ai8R3pzizjmkRpU8sGb1uCFrwgW8c0LNRg3gRdCNE92X0JoRAIvhEYk8EJoRAIvhEYk8EJoRAIvhEYk8EJoRAIvhEYk8EJoRAIvhEYk8EJoRAIvhEYk8EJoRAIvhEYk8EJoRAIvhEb+H57gia9l6jP9AAAAAElFTkSuQmCC", + "text/plain": [ + "
" + ] + }, + "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_qubits, rng)\n", + "\n", + "print(f\"Hidden shift string {hidden_shift_string}\")\n", + "\n", + "display_circuit(circuit)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "editable": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "source": [ + "## Step 2: Optimize circuits for quantum hardware execution" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "['shift 011010', 'n_qubits 6', 'seed = 12345']" + ] + }, + "execution_count": 3, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "job_tags = [f'shift {hidden_shift_string}', f'n_qubits {n_qubits}', f'seed = {random_seed}']\n", + "job_tags" + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "metadata": { + "editable": true, + "scrolled": true, + "slideshow": { + "slide_type": "" + }, + "tags": [] + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Using backend ibm_kingston\n" + ] + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAIQAAACMCAYAAAC9FwHKAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjYsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvq6yFwwAAAAlwSFlzAAAPYQAAD2EBqD+naQAADzNJREFUeJzt3XtQlGXfB/DvrhxWWEAgX4FUEBcQ6eEkoryGiqDY7PREE0rGW8uTMqkcDCtmGpq0x3nIySmaMhixHoEOlJaZQ4ILJoUKgsKqKQuiMKSJIBpHUWSv949wh0t3l10OK8jvM8OM9+G697rZ715773JfPwWMMQZC+gkfdQfI2EKBIBwKBOFQIAiHAkE4FAjCoUAQDgWCcCgQhDNqgRAIBDhw4IDe+8fExCAiImJYj9nQ0ACBQACFQjHkY2zduhW+vr7D6sd4ZnAgmpqasGnTJkgkEohEIkybNg2LFi1CRkYGuru7R6eXI2jp0qUQCAQQCAQQiUSYO3cu0tPTH3W3hqS4uBj+/v4wNzeHRCJBVlbWsI9pUCAuX74MPz8/yOVypKamoqqqCqWlpUhOTkZeXh6KioqG3SFjiI2NxbVr13DhwgWsXr0acXFxyM3NfdTdMkh9fT2kUilCQkKgUCjw+uuvY926dTh8+PDwDswMEB4ezqZPn846Ozs1blepVOp/A2A//vijevns2bMsJCSEiUQiZmdnx2JjY1lHR4d6u0wmY8899xzbunUre+KJJ5iVlRV77bXX2J07d9T75Ofns0WLFjEbGxtmZ2fHpFIpq6urU2+vr69nAFhVVZXWc1iyZAnbtGkTt87NzY29+OKLjDHGtmzZwnx8fFhOTg5zdnZm1tbWLCoqirW3t+vdjzt37rC4uDjm4ODAzM3N2cyZM1lqaqp6+61bt9jatWvV5xkSEsIUCsUgv31ecnIy8/Ly4tZFRUWx8PBwg47zIL1HiNbWVsjlcsTFxcHS0lLjPgKBQOP6rq4uhIeHw9bWFhUVFdi3bx+KiooQHx/P7XfkyBFUV1ejuLgYubm52L9/P9577z3uOJs3b8apU6dw5MgRCIVCPP/881CpVPq/AjSYPHky7t69q16+dOkSDhw4gLy8POTl5eHXX3/F9u3b9e7HJ598goMHD2Lv3r2oqanB119/DRcXF3X7VatWobm5Gfn5+Th9+jT8/f0RGhqKmzdvAgOuhYqLi7X2ubS0FGFhYdy68PBwlJaWDut3ofcIUVZWxgCw/fv3c+vt7e2ZpaUls7S0ZMnJyer1A0eIzMxMZmtry40sP//8MxMKhaypqYmx/hHCzs6OdXV1qffJyMhgYrGY9fX1aexTS0sLA8DOnTvH2BBGiHv37rEvv/ySAWA7d+5krH+EsLCw4EaEt956iy1YsEDrMR/sR0JCAlu2bBk3Yt5XUlLCrK2tWU9PD7d+9uzZbNeuXYwxxq5cucI8PDzYyZMntT6mm5sbN+qw/t8pANbd3a213WCG/SmjvLwcCoUCXl5euHPnjsZ9qqur4ePjw40sixYtgkqlQk1NjXqdj48PLCws1MtBQUHo7OzEH3/8AQC4ePEi1qxZA1dXV1hbW6tfdY2NjQb1OT09HWKxGJMnT0ZsbCySkpKwYcMG9XYXFxdYWVmplx0dHdHc3KxeHqwfMTExUCgU8PDwQGJiIuRyubrtmTNn0NnZCXt7e4jFYvVPfX09Ll26BAB48sknoVQqERgYaNB5jQQTfXeUSCQQCATcEwgArq6uQP+wO9qeffZZODs7Y/fu3XBycoJKpcJTTz3FDff6iI6ORkpKCiZPngxHR0cIhfzrwtTUlFsWCATc29Jg/fD390d9fT3y8/NRVFSE1atXIywsDN9//z06Ozvh6Oio8e1gypQpep+Dg4MDrl+/zq27fv06rK2th/Vc6B0Ie3t7LF++HDt37kRCQoLW6whNPD09kZWVha6uLnW748ePQygUwsPDQ73fmTNncPv2bfUJlZWVQSwWY8aMGWhtbUVNTQ12796N4OBgAMCxY8cMOVc1GxsbSCSSIbXVtx/W1taIiopCVFQUIiMjsXLlSty8eRP+/v5oamqCiYkJd11hqKCgIBw6dIhbV1hYiKCgoCEfE4Z+7ExPT8e9e/cQEBCA7777DtXV1aipqcFXX30FpVKJSZMmaWwXHR0NkUgEmUyG33//HUePHkVCQgJefvllTJs2Tb3f3bt3sXbtWly4cAGHDh3Cli1bEB8fD6FQCFtbW9jb2yMzMxN1dXX45ZdfsHnz5mGd/FDo04+PPvoIubm5UCqVqK2txb59++Dg4IApU6YgLCwMQUFBiIiIgFwuR0NDA06cOIGUlBScOnUKAHD16lXMmTMH5eXlWvuxfv16XL58GcnJyVAqlUhPT8fevXuRlJQ0vBM09KLjzz//ZPHx8WzWrFnM1NSUicViFhgYyHbs2MFdEA71Y+e7777L7O3tmVgsZrGxsdzFV2FhIfP09GTm5ubM29ubFRcXc48z1I+dA93/2DlQWloac3Z21rsfmZmZzNfXl1laWjJra2sWGhrKKisr1e3b29tZQkICc3JyYqampmzGjBksOjqaNTY2cudx9OhRnc/F0aNHma+vLzMzM2Ourq5sz549OvfXh4DRTbZkAPrjFuFQIAiHAkE4FAjCoUAQDgWCcCgQhEOBIBwKBOFQIAiHAkE4FAjCoUAQjt43yJDh6+tTYZj3Az9EKAQmTRq51zUFwkj6+lRwXJaLlls9I3rcqbYiXPtlzYiFgt4yjESlwoiHAfj7mCM56lAgCIcCQTgUCMKhQBAOBUKHqKgoeHt7w8/PD4GBgThy5MioPp67iw1Kv3oWkyb9PUd2W/w8JEZ7jepjPog+duqwa9cu9WyqqqoqhIaG4saNGw/N9BoptQ1tKDh+Ba//nxfyj11ByHxHLP7Xz6PyWNqM6RHir7/+gkwmg5eXF/z9/ZGYmIi1a9cCAJKSkhAcHKxzYsqxY8fwzTffPLQ+Ozsbt27dGvTxB06ta2trG/J5GOL9z89gzTOzseffixH/filUKuPOkhjTI0RkZCSkUimys7PR0tICiUSC7du3o7KyEj09PSgpKcGGDRtQVVUFPz+/h9rX1dVBoVDgpZde4ta3tbUhLCwMcrkc9vb2OvuQlJSEn376CW1tbfjhhx90jg7d3d1QKpUat/Xe0++JvdurQtnZZgT5/A8Uyla92lRWVcHURHMphjlz5nATqAczZgPx22+/obm5WT0CTJ06FU5OTpg/fz7KysqwfPlyAEBYWBhKS0s1BgIA9u3bp54iN1BtbS1eeOEFnTUYACAtLQ1paWkoKChAcnIyjh8/DjMzM437KpVKzJs3T8uRJgH/2DXIWQP/cLOFt7sdqqpbEbXSFd8VXB60zcIFCwD0adx2v/6EvsZsICoqKrBgwQL1cmtrK65evQpvb2/I5XK4u7sD/RN3H5yRPtCKFSuQkpLCrautrYVMJsObb76pd39WrlyJ+Ph4nDt3TuuTPmfOHJw+fVrjtt57DAvX6S6GJhAA6Sn/i43/OYEr17tQ/F8pDpX8gY6uXp3tyk6e1DlCGGLMBmLq1KnYv38/+vr60NfXhw0bNsDDwwNmZmawsbFBe3s7AKC9vV3nNHpNM73feecd5OTk4JlnntHa7vbt22hqasKsWbOA/ootra2t6vIHmlhYWGh9Nfb2qgDoDsTGKE+Unm3GuYt/X998mHMOqYkBSHhfd1UYfz8/mJqOzOXgmJ3b2dPTg8jISCiVSjg7O6tHhIyMDFRWVuKLL77AZ599ho0bN2LdunUGDYv6uHnzJqRSKTo6OmBiYgJLS0ts27YNy5YtG9LxentVMJu3Z0T7eN/d0/8asUCM2RFCJBIhLy9PvSyTybBkyRKgvyBHTk4OgoODMW/evBEPAwDY2dkNv17TODRmA/GgiooKvPHGG+rljz/++JH253E1bgJx4cKFR92FCWFMfzFFjI8CQTgUCMKhQBAOBcJIhMK/b4gdaVNtRRjJP76O2S+mHkfj4TZ8CgTh0FsG4VAgCIcCQTgUCMKhQBAOBYJwKBCEM27+/P04GA9fTFEgjITqQxAO1Ycg4xIFgnAoEIRDgSCcCRGIbdu2wd3dHUKhEAcOHOC2Xbp0CU8//TTc3d3h5+fHzQPV1W40jIX6EBMiEMuXL0dBQQEWL1780Lb169dDJpOhtrYWH3zwAaKjo3H/FhFd7UbDwPoQc2dPQch8R+zMNe70A6MEQludh4aGBjg4OGDp0qVYsWKF1vbDrfOwcOFCjXMyW1paUFZWhpiYGKA/AIwx9YRdbe1G04SoD6GtzgP6Z1VnZWXpbD8SdR40aWxshKOjI/d/fbu4uKCxsREBAQEGH4/qQ+hBV50HACgqKkJwcDBWrVqFxMRErccZbp0HY6D6EHrQVeeBMYaamhqYm5sjIiICoaGh8PLSfBE1UnUeBpo5cyauXbuG3t5e9SjR0NCAmTNnDul4VB9CD7rqPACAubk5AEAqleL8+fNaAzHUOg+D9S0wMBBZWVmIjY1FYWEhGGM6XuW6UX0IPeiq89DR0QErKysAwCuvvIJNmzYN+cnQZevWrfj888/R0tICKysriEQilJWVYfr06bh48SJkMhlu3LgBCwsLZGZmIjAwcNB2hhov9SGMfhv+/ToPr776KuRyOd5++22IRCIsWbIEqampxuyKUVEgtJg7dy6+/fZbeHt7G/NhH7nxEgij3w9BdR7GtgnxTSXRHwWCcCgQhEOBIBwKhJFQfQjykPFwGz4FgnDoLYNwKBCEQ4EgHAoE4VAgCIcCQTgUCMKhcgBGNB6+mKJAGAnVhyAcqg9BxiUKBOFQIAiHAkE4EyIQuuo8DLV2xONqQgRCV52HodaOGA0TpmCItvoQhYWFWLhwIYKCgvDhhx9qbT9a9SF0bRusdsRomDAFQyIjI+Hr64vz58/j8OHDyM7ORkBAADw9PVFSUoITJ07g4MGD6Ozs1Ni+rq4O5eXlD62/Xx+itVW/OgqG0FU7YjQ99gVDdNWHGDhp1sTEBEIdd4uOh/oQVDBED7rqQ9xXUFAAiUSis+OjUR9Cl6HUjqCCIXoYrD5EY2MjduzYgYMHD+o8zmjUhxis34bWjngcCoY80voQ3d3dkEql2L1790NP9kjSVedhqLUjDKXP7O+4Fz3h7CRG8kcVAIBX/inBfK+pgxYMGdflAAbWh/j000+xfft2uLm5AQBycnKGXM5nrBsv5QCoPoSRjJdAUH0IwpkQ31QS/VEgCIcCQTgUCMKhQBgJ1YcgDxkPt+FTIAiH3jIIhwJBOBQIwqFAEA4FgnAoEIRDgSAcCgThUCAIhwJBOBQIwqFAEA4FgnAoEIRDgSAcCgTh/D9i+jj1q3LNcwAAAABJRU5ErkJggg==", + "text/plain": [ + "
" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager\n", + "from qiskit_ibm_runtime import QiskitRuntimeService\n", + "\n", + "# Uncomment this to run the circuits on a quantum computer on IBMCloud.\n", + "service = QiskitRuntimeService()\n", + "backend = service.backend(\"ibm_kingston\")\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", + "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", + "isa_circuit = get_isa_circuit(circuit, backend)\n", + "display_circuit(isa_circuit)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Step 3: Execute circuits using Qiskit Primitives" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "metadata": {}, + "outputs": [], + "source": [ + "# submit job for solving the hidden shift problem using the Sampler primitive\n", + "NUM_SHOTS = 50_000\n", + "\n", + "from qiskit_ibm_runtime import SamplerV2 as Sampler\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", + "import mthree\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, + "metadata": {}, + "outputs": [], + "source": [ + "job = run_sampler(backend, isa_circuit, NUM_SHOTS)\n", + "mit, qubit_mapping = setup_mthree_mitigation(isa_circuit, backend)" + ] + }, + { + "cell_type": "markdown", + "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 may be compared directly.\n" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "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, + "metadata": {}, + "outputs": [], + "source": [ + "counts, pub_result = get_bitstring_counts(job)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The Hamming distance between two bit strings is the number of indices at which the bits differ." + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "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, + "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", + "def find_hidden_shift_bitstring(counts, hidden_shift_string):\n", + " # convert counts to probabilities\n", + " probs = {unscramble(bitstring): count / NUM_SHOTS for bitstring, count in counts.items()}\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({k: (v, hamming_distance(hidden_shift_string, k)) for k, v in sorted(probs.items(), key=lambda x: x[1], reverse=True)[:10]})\n", + "\n", + " return probs, most_probable" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "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(counts, hidden_shift_string)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Let's record the probability of the most probable bitstring before applying readout error mitigation with M3." + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "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", + "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`s that sum to one. But some values might be negative." + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [], + "source": [ + "def perform_mitigation(mit, counts, qubit_mapping):\n", + "\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 = {unscramble(k): f\"{v:.2e}\" for k, v in sorted(quasis.items(), key=lambda x: x[1], reverse=True)[:10]}\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, + "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(mit, counts, qubit_mapping)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Comparing identifying the hidden shift string before and after applying M3 correction" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "metadata": {}, + "outputs": [], + "source": [ + "def compare_before_and_after_M3(max_probability_before_M3, max_probability_after_M3, is_hidden_shift_identified):\n", + "\n", + " is_probability_improved = max_probability_after_M3 > max_probability_before_M3\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, + "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(max_probability_before_M3, max_probability_after_M3, is_hidden_shift_identified)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Plot how CPU time required by M3 scales with shots" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "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": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkgAAAHHCAYAAABEEKc/AAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjYsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvq6yFwwAAAAlwSFlzAAAPYQAAD2EBqD+naQAAX3dJREFUeJzt3Xl8TOf+B/DPzCSTyTYTWScICWKJIIJElKJCUkq1Smyl6ocqirS3aFWk7S3V5WpLo3pvua0qV6uKklZtbQlBYomQorGULIjMRPbMPL8/yNRkIYkkM5N83q9XXsw5z5zzPTlhPjnneZ4jEUIIEBEREZGB1NQFEBEREZkbBiQiIiKiMhiQiIiIiMpgQCIiIiIqgwGJiIiIqAwGJCIiIqIyGJCIiIiIymBAIiIiIiqDAYmIiIioDAYkonr23HPPwdvb29RlNGj79u2DRCLBvn37TF0KVZFEIsHixYtNXQaRAQMSUS2QSCRV+jLHD+zk5GQsXrwYFy9eNHUpZmvx4sWQSCSQSqW4cuVKufVarRa2traQSCSYOXOmYXl+fj4mT54Mf39/qFQqODg4oEuXLvjoo49QXFxcz0dhejt27GAIIothZeoCiBqCr776yuj1l19+iV27dpVb3qFDB3z++efQ6/X1XGHlkpOTER0djX79+vHK1gPY2Njgm2++wauvvmq0fPPmzRW2z8/Px+nTpzF48GB4e3tDKpXi4MGDmDt3Lg4fPoz169fXU+XmYceOHVi5cmWFISk/Px9WVvxIIvPBn0aiWjB+/Hij14cOHcKuXbvKLSfLNnjw4AoD0vr16zFkyBB89913RsudnZ1x6NAho2UvvPACVCoVVqxYgQ8//BBqtbpeaq9ISUkJ9Ho95HJ5uXW5ubmwt7evt1oUCkW97YuoKniLjaiele2DdPHiRUgkErz//vtYuXIlWrVqBTs7OwwaNAhXrlyBEAJvvfUWmjdvDltbWzz55JPIysoqt92dO3eiT58+sLe3h6OjI4YMGYLTp0/ft5a1a9di5MiRAID+/ftXeCvw008/RceOHWFjY4OmTZtixowZyM7OfuBxXrp0CS+++CLatWsHW1tbuLi4YOTIkeVu5a1duxYSiQS//vorpk2bBhcXFyiVSkyYMAG3bt0yauvt7Y0nnngCP//8MwICAqBQKODn51fpFZxSUVFRsLa2xvXr18utmzp1KpycnFBQUPDAYxo7diyOHz+Os2fPGpalp6djz549GDt27APff+9xAKjS9zE7Oxtz586Ft7c3bGxs0Lx5c0yYMAE3btwwtMnMzMTkyZPh4eEBhUKBLl264L///a/Rdu79OVu+fDlat24NGxsbwy1WiUSC5ORkjB07Fk2aNEHv3r0N7123bh26desGW1tbODs7Y/To0RXeajx8+DAGDx6MJk2awN7eHp07d8ZHH30E3P25X7lyJVDmlnSpivogJSYm4vHHH4dSqYSDgwMGDBhQLnCW/vwcOHAAkZGRcHNzg729PZ566qkKzzdRVfEKEpGZ+Prrr1FUVIRZs2YhKysLy5Ytw6hRo/DYY49h3759mDdvHs6fP49PPvkEr7zyCr744gvDe7/66itMnDgRYWFhePfdd5GXl4eYmBj07t0biYmJld46e/TRR/HSSy/h448/xmuvvYYOHToAd28F4m7fm+joaISGhmL69OlISUlBTEwMjhw5ggMHDsDa2rrS4zly5AgOHjyI0aNHo3nz5rh48SJiYmLQr18/JCcnw87Ozqj9zJkz4eTkhMWLFxv2c+nSJUOH61Lnzp1DREQEXnjhBUycOBFr1qzByJEjERsbi4EDB1ZYy7PPPos333wTGzduNOojVFRUhG+//RYjRoyo0hWMRx99FM2bN8f69evx5ptvAgA2btwIBwcHDBkypNL3FRUVQavVIj8/H0ePHsX777+Pli1bok2bNvfd3+3bt9GnTx+cOXMGzz//PAIDA3Hjxg1s3boVf/31F1xdXZGfn49+/frh/PnzmDlzJnx8fLBp0yY899xzyM7OxuzZs422uWbNGhQUFGDq1KmwsbGBs7OzYd3IkSPh6+uLd955B0IIAMA///lPvPHGGxg1ahT+7//+D9evX8cnn3yCRx99FImJiXBycgIA7Nq1C0888QQ8PT0xe/ZsqNVqnDlzBtu3b8fs2bMxbdo0XLt2rcJbzxU5ffo0+vTpA6VSiVdffRXW1tb47LPP0K9fP+zfvx/BwcFG7WfNmoUmTZogKioKFy9exPLlyzFz5kxs3LjxgfsiqpAgolo3Y8YMUdk/r4kTJ4qWLVsaXqempgoAws3NTWRnZxuWL1iwQAAQXbp0EcXFxYblY8aMEXK5XBQUFAghhMjJyRFOTk5iypQpRvtJT08XKpWq3PKyNm3aJACIvXv3Gi3PzMwUcrlcDBo0SOh0OsPyFStWCADiiy++uO928/Lyyi2Li4sTAMSXX35pWLZmzRoBQHTr1k0UFRUZli9btkwAED/88INhWcuWLQUA8d133xmWaTQa4enpKbp27WpYtnfv3nLHFBISIoKDg43q2bx5c4XHXlZUVJQAIK5fvy5eeeUV0aZNG8O6Hj16iEmTJglxJ1GIGTNmlHv/N998IwAYvrp37y5Onjx5330KIcSiRYsEALF58+Zy6/R6vRBCiOXLlwsAYt26dYZ1RUVFIiQkRDg4OAitVivEPT9nSqVSZGZmVnh8Y8aMMVp+8eJFIZPJxD//+U+j5adOnRJWVlaG5SUlJcLHx0e0bNlS3Lp1q8I6xQP+XQAQUVFRhtfDhw8XcrlcXLhwwbDs2rVrwtHRUTz66KOGZaU/P6GhoUb7mjt3rpDJZEb/poiqg7fYiMzEyJEjoVKpDK9Lf0MeP368UefV4OBgFBUV4erVq8Dd39yzs7MxZswY3Lhxw/Alk8kQHByMvXv31qieX375BUVFRZgzZw6k0r//q5gyZQqUSiV+/PHH+77f1tbW8Pfi4mLcvHkTbdq0gZOTExISEsq1nzp1qtEVqenTp8PKygo7duwwate0aVM89dRThtelt+MSExORnp5eaT0TJkzA4cOHceHCBcOyr7/+Gl5eXujbt+99j+VeY8eOxfnz53HkyBHDnw+6vda/f3/s2rULmzZtwgsvvABra2vk5uY+cF/fffcdunTpYnS8pUqvqu3YsQNqtRpjxowxrLO2tsZLL72E27dvY//+/UbvGzFiBNzc3Crc3wsvvGD0evPmzdDr9Rg1apTRz5ZarYavr6/hZysxMRGpqamYM2eO4YpS2TqrQ6fT4eeff8bw4cPRqlUrw3JPT0+MHTsWv//+O7RardF7pk6darSvPn36QKfT4dKlS9XePxF4i43IfLRo0cLodWlY8vLyqnB5af+cc+fOAQAee+yxCrerVCprVE/pB0u7du2MlsvlcrRq1eqBHzz5+flYsmQJ1qxZg6tXrxpu2QCARqMp197X19fotYODAzw9Pcv1WWrTpk25D922bdsCd/vZVNbpOSIiAnPmzMHXX3+NRYsWQaPRYPv27Zg7d261PsS7du2K9u3bY/369XBycoJara70e1/Kw8MDHh4eAIBnnnkG77zzDgYOHIhz587dt5P2hQsXMGLEiPtu+9KlS/D19TUKsbjnNmnZ8+Tj41PptsquO3fuHIQQ5c5NqdJAWxo6/f3971trVV2/fh15eXnlfvZw97j0ej2uXLmCjh07GpaX/ffTpEkT4J5/J0TVxYBEZCZkMlm1lpcGjtIpA7766qsKP2xNNXR61qxZWLNmDebMmYOQkBCoVCpIJBKMHj3aJNMcNGnSBE888YQhIH377bcoLCys0UjDsWPHIiYmBo6OjoiIiCgXTh7kmWeeweuvv44ffvgB06ZNq/b+H8a9V/YetE6v10MikWDnzp0V/hw6ODjUSY018aB/J0TVxYBEZOFat24NAHB3d0doaGi131/Z1ZOWLVsCAFJSUoxucxQVFSE1NfWB+/r2228xceJEfPDBB4ZlBQUFlY7cOnfuHPr37294ffv2baSlpWHw4MFG7c6fPw8hhFHdf/zxB3DP6LDKTJgwAU8++SSOHDmCr7/+Gl27djW6ClFVY8eOxaJFi5CWllalDsdl5efnA5VcSbtX69atkZSUdN82LVu2xMmTJ6HX642CWulIu9LzWBOtW7eGEAI+Pj6Gq3SVtQOApKSk+/5cVPVKnZubG+zs7JCSklJu3dmzZyGVSstdWSWqbeyDRGThwsLCoFQq8c4771Q4O/ODhjqXznVTNriEhoZCLpfj448/Nvot/D//+Q80Gs19R23h7m/0ZX97/+STT6DT6Spsv3r1aqP6Y2JiUFJSgscff9yo3bVr1/D9998bXmu1Wnz55ZcICAh44JxCjz/+OFxdXfHuu+9i//79NZ6nqnXr1li+fDmWLFmCoKCgStvduHGjwisY//73vwEA3bt3v+9+RowYgRMnThgdb6nS7Q4ePBjp6elGo7VKSkrwySefwMHBoVr9q8p6+umnIZPJEB0dXe44hBC4efMmACAwMBA+Pj5Yvnx5uZ+je99X2c9aWTKZDIMGDcIPP/xgdIs1IyMD69evR+/evWt865ioqngFicjCKZVKxMTE4Nlnn0VgYCBGjx4NNzc3XL58GT/++CMeeeQRrFixotL3BwQEQCaT4d1334VGo4GNjQ0ee+wxuLu7Y8GCBYiOjkZ4eDiGDRuGlJQUfPrpp+jRo8cDw8UTTzyBr776CiqVCn5+foiLi8Mvv/wCFxeXCtsXFRVhwIABGDVqlGE/vXv3xrBhw4zatW3bFpMnT8aRI0fg4eGBL774AhkZGVizZs0Dv1fW1tYYPXo0VqxYAZlMZtSxubrKDp+vyLp167Bq1SpDZ+OcnBz89NNP2LVrF4YOHfrAvkv/+Mc/8O2332LkyJF4/vnn0a1bN2RlZWHr1q1YtWoVunTpgqlTp+Kzzz7Dc889h2PHjsHb2xvffvstDhw4gOXLl8PR0bHGx9i6dWu8/fbbWLBgAS5evIjhw4fD0dERqamp+P777zF16lS88sorkEqliImJwdChQxEQEIBJkybB09MTZ8+exenTp/HTTz8BALp16wYAeOmllxAWFgaZTIbRo0dXuO+3334bu3btQu/evfHiiy/CysoKn332GQoLC7Fs2bIaHxNRlZl6GB1RQ1STYf7vvfeeUbvSoeqbNm0yWl46rPnIkSPl2oeFhQmVSiUUCoVo3bq1eO6558TRo0cfWO/nn38uWrVqJWQyWblh7ytWrBDt27cX1tbWwsPDQ0yfPr3cUO6K3Lp1S0yaNEm4uroKBwcHERYWJs6ePStatmwpJk6cWO549u/fL6ZOnSqaNGkiHBwcxLhx48TNmzeNttmyZUsxZMgQ8dNPP4nOnTsLGxsb0b59+3Lfo4qG+ZeKj48XAMSgQYMeeAyl7h3mfz9lh/kfOXJEjBw5UrRo0ULY2NgIe3t7ERgYKD788EOjqRvu5+bNm2LmzJmiWbNmQi6Xi+bNm4uJEyeKGzduGNpkZGQYvtdyuVx06tRJrFmzxmg7lf2cVeX4vvvuO9G7d29hb28v7O3tRfv27cWMGTNESkqKUbvff/9dDBw4UDg6Ogp7e3vRuXNn8cknnxjWl5SUiFmzZgk3NzchkUiM/o2UHeYvhBAJCQkiLCxMODg4CDs7O9G/f39x8OBBozb3+/dQlSkciCojEezBRkQmtHbtWkyaNAlHjhx54C0nb29v+Pv7Y/v27TXe34kTJxAQEIAvv/wSzz77bI23Q0QNG/sgEVGj8vnnn8PBwQFPP/20qUshIjPGPkhE1Chs27YNycnJWL16NWbOnFmvD2IlIsvDgEREjcKsWbOQkZGBwYMHIzo62tTlEJGZYx8kIiIiojLYB4mIiIioDAYkIiIiojLYB6mG9Ho9rl27BkdHxxo9rZqIiIjqnxACOTk5aNq06X2fo8iAVEPXrl3js4CIiIgs1JUrV9C8efNK1zMg1VDp9P1XrlzhM4GIiIgshFarhZeX1wMfw8OAVEOlt9WUSiUDEhERkYV5UPcYdtImIiIiKoMBiYiIiKgMBiQiIiKiMhiQiIiIiMpgQCIiIiIqgwGJiIiIqAwGJCIiIqIyGJCIiIiIymBAIiIiIiqDM2kTERGR2dDpBeJTs5CZUwB3RwWCfJwhk9b/Q+EZkIiIiMgsxCalIXpbMtI0BYZlnioFoob6Idzfs15r4S02IiIiMrnYpDRMX5dgFI4AIF1TgOnrEhCblFav9TAgERERkUnp9ALR25IhKlhXuix6WzJ0+opa1A0GJCIiIjKp+NSscleO7iUApGkKEJ+aVW81MSARERGRSWXmVB6OatKuNjAgERERkUm5OypqtV1tYEAiIiIikwrycYanqvLwI7k7mi3Ix7neamJAIiIiIpPR6wVkUgmihvpBcjcM3av0ddRQv3qdD4kBiYiIiEwiO68II1YdxC/JGQj390TM+ECoy1xJUqsUiBkfWO/zIHGiSCIiIqp3t3KLMO7fh5GcpsUbPySht68rwv09MdBPzZm0iYiIqPHJuhuOzqRp4eogx3+fD4LCWgYAkEklCGntYuoSGZCIiIio/ty8XYhx/z6Ms+k5cHWwwTdTguHr4WjqssphQCIiIqJ6ceN2IcZ9fhgpGTlwc7TBN1N6oo27g6nLqhADEhEREdWLr+IuISUjB+6ONvhmak+0djPPcAQGJCIiIqovLw3wRV5RCcYEtUArMw5HYEAiIiKiupSVWwSlwgpWMilkUgleH+Jn6pKqhPMgERERUZ3I1BbgmVUHEfm/EyjR6U1dTrXwChIRERHVunRNAcZ+fgh/3shFQZEON3OL4KGsv2epPSwGJCIiIqpVaZp8jFl9CBdv5qGZky02TO1pUeEIDEhERERUm65l52PM54dw6WYemjexxTdTesLL2c7UZVUbAxIRERHViqvZd64cXc7Kg5fznXDUvInlhSMwIBEREVFtuXQjF+naArRwtsM3U3uimZOtqUuqMQYkIiIiqhW92rhizXM94ONqj6YWHI7AgERERNQw6PQC8alZyMwpgLujAkE+zpBJJXW+3ytZeSjS6Q2zYj/SxrXO91kfGJCIiIgsXGxSGqK3JSNNU2BY5qlSIGqoH8L9PWttP2VDmFqpwPj/HEaxTo+N00Lg42pfa/syNQYkIiIiCxablIbp6xIgyixP1xRg+roExIwPrJWQVFEIk0oAvQBaudnDXi576H2YE86kTUREZKF0eoHobcnlwhEAw7LobcnQ6StqUXWlIezecATcCUcAMO3RVnC3sHmOHoRXkIiIiCxUfGpWudByLwEgTVOAx97fhyb2cshlUsitpHi8kxrjglsCADT5xVi68yxsrKSwlkkgt5JCLpPB2koCuUyK9h6OlYawUst/OYdnunnVS5+n+sKAREREZKEycyoPR/e6lJWHS1l5htd+TZWGv2vzi/FN/OVK3xvawf2+IQx3Q1h8ahZCWrtUqR5LwIBERERkodwdq3Zba154O7Rxd0RRiR5FOh3auDka1jkqrBA5sC2KSvQo1ulRWKJHkU6P4rt/2lWxb1FVw5qlYEAiIiKyUEE+znB3tEFmTmGF6yUA1CoFpj7autLbX052crw0wLfSfcRduIlv4q88sJaqhjVLwU7aREREFkoIAaWtdYXrSuNQ1FC/h+obFOTjDE+VApVtQXJ3SoEgH+ca78McMSARERFZqPd//gPnM29DYSWFm4PcaJ1apaiVIf4yqQRRQ/2Ae0JXqdoKYeaIt9iIiIgs0O4zGVi1/wIA4INRAQj3V9fZTNrh/p6IGR9Ybh4kdR1MRmkuGJCIiIgszLXsfET+7wQA4Lle3hjS+U5AqctRZOH+nhjoV3chzNwwIBEREVkYN0cbjA7ywuE/s/Da4A71tl+ZVNKghvLfDwMSERGRhbGWSbHg8Q4oLNFBbsXuxHWB31UiIiILkXRVg6ISveG1jVXDev6ZOWFAIiIisgB/Xr+N0asPIWJ1HLJyi0xdToPHgERERGTmCop1ePHrBNwuLIFcJoVSwR4ydY0BiYiIyMwt+iEJZ9Nz4Opgg0/GdIWVjB/fdY3fYSIiIjO26egV/O/oX5BKgI/HBMBd2bAe6WGuGJCIiIjM1Nl0Ld74IQkAMDe0LXq1djV1SY2GWQSklStXwtvbGwqFAsHBwYiPj79v+02bNqF9+/ZQKBTo1KkTduzYYVhXXFyMefPmoVOnTrC3t0fTpk0xYcIEXLt2zWgb3t7ekEgkRl9Lly6ts2MkIiKqDiEE5n17EgXFevRt64YZ/duYuqRGxeQBaePGjYiMjERUVBQSEhLQpUsXhIWFITMzs8L2Bw8exJgxYzB58mQkJiZi+PDhGD58OJKS7iTsvLw8JCQk4I033kBCQgI2b96MlJQUDBs2rNy23nzzTaSlpRm+Zs2aVefHS0REVBUSiQTLR3fFY+3d8a+IAEgb6IzV5koihBCmLCA4OBg9evTAihUrAAB6vR5eXl6YNWsW5s+fX659REQEcnNzsX37dsOynj17IiAgAKtWrapwH0eOHEFQUBAuXbqEFi1aAHevIM2ZMwdz5sypUd1arRYqlQoajQZKpbJG2yAiIqL6VdXPb5NeQSoqKsKxY8cQGhr6d0FSKUJDQxEXF1fhe+Li4ozaA0BYWFil7QFAo9FAIpHAycnJaPnSpUvh4uKCrl274r333kNJSUml2ygsLIRWqzX6IiIiqm0nrmTj93M3TF1Go2fSiRRu3LgBnU4HDw8Po+UeHh44e/Zshe9JT0+vsH16enqF7QsKCjBv3jyMGTPGKCm+9NJLCAwMhLOzMw4ePIgFCxYgLS0NH374YYXbWbJkCaKjo2twlERERFWTnVeEF79OwDVNPlaODcTgTp6mLqnRatAzTRUXF2PUqFEQQiAmJsZoXWRkpOHvnTt3hlwux7Rp07BkyRLY2NiU29aCBQuM3qPVauHl5VXHR0BERI2FXi/w8v9O4Gp2Plo426G3L0esmZJJA5KrqytkMhkyMjKMlmdkZECtVlf4HrVaXaX2peHo0qVL2LNnzwP7CQUHB6OkpAQXL15Eu3btyq23sbGpMDgRERHVhs9+/RO7z2ZCbiXFp+MCoVRYm7qkRs2kfZDkcjm6deuG3bt3G5bp9Xrs3r0bISEhFb4nJCTEqD0A7Nq1y6h9aTg6d+4cfvnlF7i4uDywluPHj0MqlcLd3f2hjomIiKi6Dv95E+//nAIAWDy0I/ybqUxdUqNn8ltskZGRmDhxIrp3746goCAsX74cubm5mDRpEgBgwoQJaNasGZYsWQIAmD17Nvr27YsPPvgAQ4YMwYYNG3D06FGsXr0auBuOnnnmGSQkJGD79u3Q6XSG/knOzs6Qy+WIi4vD4cOH0b9/fzg6OiIuLg5z587F+PHj0aRJExN+N4iIqLG5nlOIWd8kQqcXeKprM4wJYvcNc2DygBQREYHr169j0aJFSE9PR0BAAGJjYw0dsS9fvgyp9O8LXb169cL69euxcOFCvPbaa/D19cWWLVvg7+8PALh69Sq2bt0KAAgICDDa1969e9GvXz/Y2Nhgw4YNWLx4MQoLC+Hj44O5c+ca9TEiIiKqD98n/oXMnEL4ujvgn0/5QyLhfEfmwOTzIFkqzoNERES1QQiBDUeuoHvLJvD1cDR1OQ1eVT+/TX4FiYiIqDGTSCQYE9TC1GVQGSZ/1AgREVFjk6bJR+T/jiM7r8jUpVAleAWJiIiojun0AvGpWcjMKYCLvRwf7voDCZezkVtYgs+e7W7q8qgCDEhERER1KDYpDdHbkpGmKTBarrCW4vXBfiari+6Pt9iIiIjqSGxSGqavSygXjgCgoFiP5DSNSeqiB2NAIiIiqgM6vUD0tmRUNlRcAiB6WzJ0eg4mN0cMSERERHUgPjWrwitHpQSANE0B4lOz6rUuqhoGJCIiojqQmVN5OKpJO6pfDEhERES1TK8XuHD9dpXaujsq6rweqj6OYiMiIqpFp69psHBLEhIvZ0OpsEJOQUmF/ZAkANQqBYJ8nE1QJT0IAxIREVEtyCkoxoe7/sB/D16EXgD2chme6NwU38RfhuRun6NSpU9bixrqB5mUz14zRwxIRERED0EIge0n0/DW9mRk5hQCAJ7o7ImFQ/ygVinwaFvXcvMgqVUKRA31Q7i/pwkrp/thQCIiInoIS3aexepf/wQAeLvY4c0n/fFoWzfD+nB/Twz0Uxtm0nZ3vHNbjVeOzBsDEhER0UN4MqAp1h26hGmPtsa0vq2gsJaVayOTShDS2sUk9VHNMCARERFVw96UTFzIvI3/69MKANCxqQpx8wdAZWdt6tKoFjEgERERVcG17Hy8tT0ZO5PSYSWV4NG2bmjr4QgADEcNEAMSERHRfRTr9FhzIBXLfzmHvCIdZFIJJj3ijaZOtqYujeoQAxIRETVqOr2otAP1kYtZWPh9ElIycgAA3Vo2wdvD/dHBU2niqqmuMSAREVGjFZuUVm4IvufdIfi92rhi0pojuF1YgiZ21lgwuAOeCWwOKUefNQoSIQQfI1wDWq0WKpUKGo0GSiV/kyAisjSxSWmYvi6h3CzXpfEnZnwgrmYX4HxmDl4Na48m9nITVEm1raqf37yCREREjY5OLxC9LbnCR4CIuyEpelsyfp/3GOcraqT4sFoiImp04lOzjG6rlSUApGkKEJ+aVa91kflgQCIiokYnM6fycFSTdtTwMCAREVGj4+6oqNV21PAwIBERUaPTqZkKclnlH4GSu6PZgnyc67UuMh8MSERE1OjYyWUI8FJVuK60S3bUUD920G7EGJCIiKjRkUol+HpKT0QN9YOnyvg2mlqlQMz4QIT7e5qsPjI9zoNUQ5wHiYjIsvx1Kw9fxV3Cq+Htja4M3W8mbWp4OA8SERHRXReu38az/z6Ma5oCyK2keHlQO8M6mVSCkNYuJq2PzA8DEhERNWinr2kw4T/xuJlbhDbuDhgX3NLUJZEFYEAiIqIG69ilLDy35ghyCkrQsakSXz4fBBcHG1OXRRaAAYmIiBqk38/dwJQvjyK/WIce3k3wn+d6QKmwNnVZZCEYkIiIqMHR5BVj+rpjyC/WoY+vK1Y/2x22cpmpyyILwmH+RETU4KjsrPFhRACe6OyJf09kOKLq4xUkIiJqMLQFxYbbaAP9PDDQz8PUJZGF4hUkIiJqED7ddx6DPvwVV7LyTF0KNQAMSEREZNGEEHg39iyWxaYgXVuAn06nm7okagB4i42IiCyWXi8QtfU0vjp0CQCw4PH2+L8+rUxdFjUADEhERGSRSnR6vPrtSWxOvAqJBHh7uD8ngaRaw4BEREQWp7BEh1nrE/FzcgZkUgk+HNUFTwY0M3VZ1IAwIBERkcUpLNHjmiYfcispPh0biFCOVqNaxoBERERmS6cXiE/NQmZOAdwdFQjycYZMKoFSYY3/TgrCheu5CPJxNnWZ1AAxIBERkVmKTUpD9LZkpGkKDMtUttZ4d0QnhPt7wsXBhs9VozrDYf5ERGR2YpPSMH1dglE4AgBNfjFeWJeA2KQ0k9VGjQMDEhERmRWdXiB6WzLEfdpEb0uGTn+/FkQPhwGJiIjMSnxqVrkrR2WlaQoQn5pVbzVR48OAREREZiUz5/7hqLrtiGqCAYmIiMyKu6OiVtsR1QQDEhERmZUgH2d4qioPPxIAnioFh/dTnWJAIiIik7tdWIKFW07hr1t5kEkliBrqB8ndMHSv0tdRQ/0gk5ZdS1R7GJCIiMikkq5qMPST37Hu0GVE/u8EhBAI9/dEzPhAqMtcSVKrFIgZH4hwf0+T1UuNAyeKJCIikxBC4IsDF7F05xkU6wSaqhR4ZVA7SCR3rgyF+3tioJ+6wpm0ieoaAxIREdW7m7cL8cqmE9ibch0AMMjPA8ue6QwnO7lRO5lUgpDWLiaqkhozBiQiIqpX5zJyMO7fh5GZUwi5lRRvDOmA8T1bGq4cEZkDBiQiIqpXXs52cLaXQ2lrjRVju6K9WmnqkojKYUAiIqI6l6bJh7ujAjKpBAprGf7zXA8428lhK5eZujSiCpnFKLaVK1fC29sbCoUCwcHBiI+Pv2/7TZs2oX379lAoFOjUqRN27NhhWFdcXIx58+ahU6dOsLe3R9OmTTFhwgRcu3bNaBtZWVkYN24clEolnJycMHnyZNy+fbvOjpGIqLH68WQaBv3rV6zce96wrJmTLcMRmTWTB6SNGzciMjISUVFRSEhIQJcuXRAWFobMzMwK2x88eBBjxozB5MmTkZiYiOHDh2P48OFISkoCAOTl5SEhIQFvvPEGEhISsHnzZqSkpGDYsGFG2xk3bhxOnz6NXbt2Yfv27fj1118xderUejlmIqLGIL9Ih/nfncSM9QnIKSjBb+euo0SnN3VZRFUiEUKY9HHIwcHB6NGjB1asWAEA0Ov18PLywqxZszB//vxy7SMiIpCbm4vt27cblvXs2RMBAQFYtWpVhfs4cuQIgoKCcOnSJbRo0QJnzpyBn58fjhw5gu7duwMAYmNjMXjwYPz1119o2rTpA+vWarVQqVTQaDRQKnn/nIjoXmfStJj1TSLOZ96GRAK82K815oS2hbXM5L+XUyNX1c9vk/6kFhUV4dixYwgNDf27IKkUoaGhiIuLq/A9cXFxRu0BICwsrNL2AKDRaCCRSODk5GTYhpOTkyEcAUBoaCikUikOHz5cC0dGRNTw6fQCcRdu4ofjVxF34SZ0egEhBL6Mu4gnVx7A+czbcHe0wbrJwfhHWHuGI7IoJu2kfePGDeh0Onh4eBgt9/DwwNmzZyt8T3p6eoXt09PTK2xfUFCAefPmYcyYMYakmJ6eDnd3d6N2VlZWcHZ2rnQ7hYWFKCwsNLzWarVVPEoiooYnNikN0duSkaYpMCzzVCkwo38b/PPHMygq0aN/Oze8P7ILXBxsTForUU006FFsxcXFGDVqFIQQiImJeahtLVmyBNHR0bVWGxGRpYpNSsP0dQko2z8jXVOAN7YkYVxwC3i72mNybx/ObUQWy6TXO11dXSGTyZCRkWG0PCMjA2q1usL3qNXqKrUvDUeXLl3Crl27jO4zqtXqcp3AS0pKkJWVVel+FyxYAI1GY/i6cuVKtY+XiKi+VHT7q7a2G70tuVw4AmBYtvtsJiY9wnBEls2kAUkul6Nbt27YvXu3YZler8fu3bsREhJS4XtCQkKM2gPArl27jNqXhqNz587hl19+gYuLS7ltZGdn49ixY4Zle/bsgV6vR3BwcIX7tbGxgVKpNPoiIjJHsUlp6P3uHoz5/BBmbziOMZ8fQu939yA2Ke2htnv8SjY+3JVidFutLAEgTVOA+NSsh9oXkamZ/BZbZGQkJk6ciO7duyMoKAjLly9Hbm4uJk2aBACYMGECmjVrhiVLlgAAZs+ejb59++KDDz7AkCFDsGHDBhw9ehSrV68G7oajZ555BgkJCdi+fTt0Op2hX5GzszPkcjk6dOiA8PBwTJkyBatWrUJxcTFmzpyJ0aNHV2kEGxGRubrf7a/p6xIQMz4Q4f6ehuW5hSVI0xQgXVOANE3+nT+1BcjQFKCwRI91//f3L41LdpzB4SoGn8ycykMUkSUweUCKiIjA9evXsWjRIqSnpyMgIACxsbGGjtiXL1+GVPr3ha5evXph/fr1WLhwIV577TX4+vpiy5Yt8Pf3BwBcvXoVW7duBQAEBAQY7Wvv3r3o168fAODrr7/GzJkzMWDAAEilUowYMQIff/xxPR45EVHtqsrtr5f/dwID/dSQSe/c/pr83yM49GfFoUciAYp1esPoswAvJ+QUlCA57cGDVNwdFQ9xJESmZ/J5kCwV50EiInMTd+Emxnx+6IHtvno+CH3augEA5m48jl/OZECtVECtUsBTpYBaZXv3TwV6t3E1Gp6v0wv0fncP0jUFFQYxCQC1SoHf5z1mCGFE5qSqn98mv4JERES1o6q3tTK0f7db9kznas1PJJNKEDXUD9PXJUByz5Up3A1HABA11I/hiCweZ+0iImogqnpbq1kTO8PfazJ5Y7i/J2LGB0KtMt6fWqUo18eJyFLxChIRUQMR5OMMT5Xigbe/gnycH3pf4f6eGOinRnxqFjJzCuDueGe7vHJEDQUDEhFRA1Hft79kUglCWrtUoSWR5eEtNiKiBuL7xL/gbG/D219EtYBXkIiIGoDUG7lYsPkUCkv02DQtBL/Pe4y3v4geAgMSEZGF0+kF/rHpBAqK9XikjQsCWzSBlLe/iB4Kb7EREVm4NQdScfTSLdjLZXh3RGdIeaWI6KExIBERWbA/r9/Gez+lAABeH+KH5vcM4SeimmNAIiKyUDq9wCubTqCwRI8+vq4YE+Rl6pKIGgwGJCIiC/Xz6XQkXM6Gg40Vlo7oDImEt9aIags7aRMRWahwfzU+GNkFUinQzMnW1OUQNSgMSEREFkoikWBEt+amLoOoQeItNiIiC7P3bCY0ecWmLoOoQWNAIiKyIOcycjBt3TEMWr4faZp8U5dD1GAxIBERWYgSnR6vbDqBohI9/DyVUCsVVXgXEdUEAxIRkYVY/dufOPGXBo4KKyx5mqPWiOoSAxIRkQVISc/B8l3nAACLh3Ys9zBaIqpdDEhERGauuPTWmk6PAe3d8XRgM1OXRNTgMSAREZm5NQdSceqqBipba7zzdCfeWiOqB5wHiYjIzI3q7oXka1r0becGD3bMJqoXDEhERGbOyU6O5aO7mroMokaFt9iIiMxUSnoOhBCmLoOoUWJAIiIyQ6evaTDk498w5ctjKCjWmbocokanRrfYiouLkZ6ejry8PLi5ucHZ2bn2KyMiaqSKSvR4ZdNJlOgFrKQS2Fjxd1mi+lblf3U5OTmIiYlB3759oVQq4e3tjQ4dOsDNzQ0tW7bElClTcOTIkbqtloioEVi59zzOpGnRxM4abw3356g1IhOoUkD68MMP4e3tjTVr1iA0NBRbtmzB8ePH8ccffyAuLg5RUVEoKSnBoEGDEB4ejnPnztV95UREDVDSVQ1W7j0PAHhruD/cHG1MXRJRo1SlW2xHjhzBr7/+io4dO1a4PigoCM8//zxWrVqFNWvW4LfffoOvr29t10pE1KDdubV2AiV6gcGd1Hiic1NTl0TUaEkEh0jUiFarhUqlgkajgVKpNHU5RNQAfPBzCj7Zcx7O9nL8PPdRuDrw6hFRbavq5/dD9/zTarXYsmULzpw587CbIiJq1Hq1doWXsy3eetKf4YjIxKodkEaNGoUVK1YAAPLz89G9e3eMGjUKnTt3xnfffVcXNRIRNQohrV2wa25fDOnsaepSiBq9agekX3/9FX369AEAfP/99xBCIDs7Gx9//DHefvvtuqiRiKhByykoNvxdYS0zaS1EdEe1A5JGozHMexQbG4sRI0bAzs4OQ4YM4eg1IqIq0OkF4i7cxA/Hr2LdoUsIeWc3voy7yFmzicxItSeK9PLyQlxcHJydnREbG4sNGzYAAG7dugWFgg9RJCK6n9ikNERvS0aapsBo+dbj1zAhxNtkdRGRsWpfQZozZw7GjRuH5s2bo2nTpujXrx9w99Zbp06d6qJGIqIGITYpDdPXJZQLRwBw7NItxCalmaQuIiqvRsP8jx07hsuXL2PgwIFwcHAAAPz4449wcnLCI488Uhd1mh0O8yei6tDpBXq/u6fCcAQAEgBqlQK/z3sMMilnziaqK1X9/K7Rs9i6deuGbt26GS0bMmRITTZFRNQoxKdmVRqOAEAASNMUID41CyGtXeq1NiIqr0q32JYuXYr8/PwqbfDw4cP48ccfH7YuIqIGJTOn8nBUk3ZEVLeqFJCSk5PRokULvPjii9i5cyeuX79uWFdSUoKTJ0/i008/Ra9evRAREQFHR8e6rJmIyOK4O1ZtEEtV2xFR3arSLbYvv/wSJ06cwIoVKzB27FhotVrIZDLY2NggLy8PANC1a1f83//9H5577jmOZiMiusfxK9mwt5HBU6VAuqYAFXX8LO2DFOTjbIIKiaisanfS1uv1OHnyJC5duoT8/Hy4uroiICAArq6udVelGWInbSJ6kBKdHp/uu4CPdp9DSxc7zB7gizkbjgN3+xyVKu2SHTM+EOH+nEWbqC7VWSdtqVSKgIAABAQEPGyNREQN1pWsPMzdeBxHL90CAPh5KtG/vTtixgeWmwdJrVIgaqgfwxGRGanRKDYiIqqYEALfJVzF4q2ncbuwBI42VnhzeEcMD2gGiUSCcH9PDPRTIz41C5k5BXB3vHNbjUP7icwLAxIRUS3JL9LhlU0n8OOpOxM+9vBugg9HBcDL2c6onUwq4VB+IjPHgEREVEtsrKTQ5BfDSirB3IFt8ULf1rwyRGShGJCIiB5CYYkOej1gK5dBKpXg/ZFdkKEtQBcvJ1OXRkQPodrPYit1/vx5/PTTT4YJJPkUaiJqbP7IyMGTKw7gze3JhmVqlYLhiKgBqHZAunnzJkJDQ9G2bVsMHjwYaWl37rVPnjwZL7/8cl3USERkVvR6gTUHUvHEJ7/jbHoOfj6djpu3C01dFhHVomoHpLlz58LKygqXL1+Gnd3fHQ8jIiIQGxtb2/UREZmVTG0Bnlt7BNHbklFUoke/dm7YOacPXBxsTF0aEdWiavdB+vnnn/HTTz+hefPmRst9fX1x6dKl2qyNiMisxCalY8Hmk7iVVwwbKykWDumA8T1bQiJhR2yihqbaASk3N9foylGprKws2NjwNygismw6vahwjqKcgmK89v0p3MorRsemSnw0OgBt3PncSaKGqtoBqU+fPvjyyy/x1ltvAQAkEgn0ej2WLVuG/v3710WNRET1IjYprdws1573zHK99OlOSLicjciBbSG3qvEYFyKyANV+FltSUhIGDBiAwMBA7NmzB8OGDcPp06eRlZWFAwcOoHXr1nVXrRnhs9iIGpbYpDRMX5dQ6YNk+Zw0ooahqp/f1f4VyN/fH3/88Qd69+6NJ598Erm5uXj66aeRmJjYaMIRETUsOr1A9LbkCsNRqehtydDpOZ0JUWNRo4kiVSoVXn/99dqvhojIBOJTs4xuq5UlAKRpChCfmsVHhBA1EjUKSAUFBTh58iQyMzOh1+uN1g0bNqy2aiMiqheZOZWHo5q0IyLLV+2AFBsbiwkTJuDGjRvl1kkkEuh0utqqjYioXqhsravUzt1RUee1EJF5qHYfpFmzZmHkyJFIS0uDXq83+qpJOFq5ciW8vb2hUCgQHByM+Pj4+7bftGkT2rdvD4VCgU6dOmHHjh1G6zdv3oxBgwbBxcUFEokEx48fL7eNfv36QSKRGH298MIL1a6diBqGPr5usJPLKl0vuTuaLcjHuV7rIiLTqXZAysjIQGRkJDw8PB565xs3bkRkZCSioqKQkJCALl26ICwsDJmZmRW2P3jwIMaMGYPJkycjMTERw4cPx/Dhw5GUlGRok5ubi969e+Pdd9+9776nTJmCtLQ0w9eyZcse+niIyHIUluiQnVcEAJBJJVj6dCfgbhi6V+nrqKF+kEk5ISRRY1HtYf7PP/88HnnkEUyePPmhdx4cHIwePXpgxYoVAAC9Xg8vLy/MmjUL8+fPL9c+IiICubm52L59u2FZz549ERAQgFWrVhm1vXjxInx8fJCYmIiAgACjdf369UNAQACWL19e49o5zJ/Icl3JysOM9QlwsLHCV5ODDcHnQfMgEZHlq+rnd7X7IK1YsQIjR47Eb7/9hk6dOsHa2vje/UsvvVSl7RQVFeHYsWNYsGCBYZlUKkVoaCji4uIqfE9cXBwiIyONloWFhWHLli3VPQx8/fXXWLduHdRqNYYOHYo33nijwhnCSxUWFqKw8O+HUWq12mrvk4hM76fT6Xhl0wnkFJTAyc4aqTdy0cbdAQAQ7u+JgX7qCmfSJqLGpdoB6ZtvvsHPP/8MhUKBffv2GT2DSCKRVDkg3bhxAzqdrtytOg8PD5w9e7bC96Snp1fYPj09vVrHMHbsWLRs2RJNmzbFyZMnMW/ePKSkpGDz5s2VvmfJkiWIjo6u1n6IyHwUlejxbuxZ/Of3VABAYAsnfDI2EM2cbI3ayaQSDuUnouoHpNdffx3R0dGYP38+pFLLnGp/6tSphr936tQJnp6eGDBgAC5cuFDpZJcLFiwwunql1Wrh5eVVL/US0cO5mp2PGV8n4PiVbADAlD4+eDW8Paxllvl/GBHVvWoHpKKiIkRERDx0OHJ1dYVMJkNGRobR8oyMDKjV6grfo1arq9W+qoKDgwEA58+frzQg2djY8GG8RBZICIGXvknE8SvZUCqs8P7ILhjU8eH+zyCihq/aKWfixInYuHHjQ+9YLpejW7du2L17t2GZXq/H7t27ERISUuF7QkJCjNoDwK5duyptX1WlUwF4erITJlFDI5FI8PZwfwT5OOPHl/owHBFRlVT7CpJOp8OyZcvw008/oXPnzuU6aX/44YdV3lZkZCQmTpyI7t27IygoCMuXL0dubi4mTZoEAJgwYQKaNWuGJUuWAABmz56Nvn374oMPPsCQIUOwYcMGHD16FKtXrzZsMysrC5cvX8a1a9cAACkpKcDdq09qtRoXLlzA+vXrMXjwYLi4uODkyZOYO3cuHn30UXTu3Lm63w4iMkNpmnwkXMrGkM53funp4KnExqk9jfpMEhHdT7UD0qlTp9C1a1cAMJp/CHd/U6uOiIgIXL9+HYsWLUJ6ejoCAgIQGxtr6Ih9+fJlo1t5vXr1wvr167Fw4UK89tpr8PX1xZYtW+Dv729os3XrVkPAAoDRo0cDAKKiorB48WLI5XL88ssvhjDm5eWFESNGYOHChdX9VhCRGdqXkom5G48jp6AEnk4KBLZoAtTg/yciatyqPQ8S3cF5kIjMS4lOj3/98gdW7r0AAOjYVImVYwPh7Wpv6tKIyIzU2TxIRETmJkNbgJe+ScTh1CwAwPieLbBwiB8U1pU/PoSI6H6qFJCefvpprF27FkqlEk8//fR9295vLiEiooeh04tykzjGXbiJ2RsScTO3CPZyGZaM6IxhXZqaulQisnBVCkgqlcpw/16lUtV1TURE5VT2GJBHWrvgZm4R2qsd8em4QLRyczBpnUTUMFS5D9Kbb76JV1555b6P42hM2AeJqP7EJqVh+roElP3PqrTbdUSQFxYP7chbakT0QFX9/K7yPEjR0dG4fft2bdVHRFQlOr1A9LbkcuEIgGHZ/pTrnBWbiGpVlf9H4WA3IjKF+NQso9tqZQkAaZoCxN/toE1EVBuq9SsX5xEhovqWmVN5OKpJOyKiqqjWMP+2bds+MCRlZfG3OCKqPecycqrUzt1RUee1EFHjUa2AFB0dzVFsRFRvvvg9FSvuTvxYGQkAterOkH8iotpSrYA0evRouLu71101RET3GNqlKT779QK6t3TGjlNpwD0ds3HPKLaooX6QSdkFgIhqT5X7ILH/ERHVNZ1eYPeZDMNrN0cb7H65H1aOC0TM+ECoVca30dQqBWLGByLc39ME1RJRQ1blK0gcxUZEdSlNk4/IjScQ9+dNfDKmK4benQ3bwebOf1Ph/p4Y6KcuN5M2rxwRUV2ockDS6/V1WwkRNVqxSemYv/kksvOKYSeXQV/JL2QyqQQhrV3qvT4ianz4sFoiMpn8Ih3e+jEZ6w9fBgB0bq7CR6O7wsfV3tSlEVEjx4BERCaRfE2LWd8k4ML1XEgkwLRHWyNyYFvIrTgjNhGZHgMSEZlEhrYAF67nwt3RBv+KCMAjbVxNXRIRkQEDEhHVG51eGDpV92/vjmXPdEZoBw8428tNXRoRkRFeyyaierE3JRMDP9yPv27lGZaN6u7FcEREZokBiYjqVEGxDtHbTmPSmiP480YuVuw5b+qSiIgeiLfYiKhW6PSi3BxFf16/jZc2HMeZNC0A4Lle3pj/eHtTl0pE9EAMSET00GKT0hC9LRlpmgLDMqXCCvnFOhTrBFzs5XhvZGc81t7DpHUSEVUVAxIRPZTYpDRMX5eAslM7agtKAAAdPB3x3+eD4O6oqPD9RETmiH2QiKjGdHqB6G3J5cLRvbLziuFib1OPVRERPTwGJCKqsfjULKPbahVJ0xQgPjWr3moiIqoNDEhEVGOZOfcPR9VtR0RkLhiQiKhGinV6/PrH9Sq1Zf8jIrI07KRNRNV2LiMHkf87gVNXNfdtJwGgVt0Z8k9EZEl4BYmIquWH41cx5JPfceqqBipba0zu7QPJ3TB0r9LXUUP9DI8XISKyFLyCRETV4uepBAD0a+eGd0d0hodSgR7eTcrNg6RWKRA11A/h/p4mrJaIqGYkQoj7jdClSmi1WqhUKmg0GiiVSlOXQ1RnhBA4dVWDzs2dDMtS0nPQ1sMBEsnfV4YqmkmbV46IyNxU9fObV5CIqFIZ2gLM++4kfv3jOr6d3guBLZoAANqpHcu1lUklCGntYoIqiYhqH/sgEVGFtp64hkH/+hX7Uq7DSibFn9dzTV0SEVG94RUkIjKSlVuEN35Iwo8n0wAAnZqp8OGoLvD1KH/ViIiooWJAIiKDvWcz8Y9vT+LG7UJYSSWY+VgbzOjfBtYyXmwmosaFAYmIDK5p8nHjdiF83R3w4agAdGquMnVJREQmwYBE1MjlFpbA3ubOfwVjg1pAAgmeDmwGhbXM1KUREZkMr5sTNQI6vUDchZv44fhVxF24CZ1eoKBYh+htpxH+0a/IKSgGAEgkEowNbsFwRESNHq8gETVwsUlp5SZxdLGXw0omQYa2EACwKzkDTwc2N2GVRETmhQGJqAGLTUrD9HUJKDsb7M3cIgCAUmGFj8Z0Rf927iapj4jIXPEWG1EDpdMLRG9LLheO7mUrl+FRX7d6rIqIyDIwIBE1UPGpWUa31SqSoS1EfGpWvdVERGQpGJCIGqjMnPuHo+q2IyJqTBiQiBqoEn3VnkPt7qio81qIiCwNO2kTNUA7TqXhje9P3beNBIBapUCQj3O91UVEZCkYkIgakGKdHkt3nsV/fk8FAPi6O+Bc5m1IAKPO2pK7f0YN9YNMKqlwW0REjRlvsRE1EOmaAoxZfcgQjqb1bYWds/tg1fhAqFXGt9HUKgVixgci3N/TRNUSEZk3XkEiagAO/XkTM9cn4MbtIjjaWOH9UV0Q1lENAAj398RAPzXiU7OQmVMAd8c7t9V45YiIqHIMSEQNQGGJHjdzi9Be7YhV47vB29XeaL1MKkFIaxeT1UdEZGkYkIgslBACEsmdq0B927ph9bPd0buNK2zlfI4aEdHDYh8kIguUfE2LJ1cewOWbeYZlA/08GI6IiGoJAxKRhdl09Aqe+vQATv6lwVs/Jpu6HCKiBom32IgsREGxDtHbTuOb+CsAgH7t3PDeM51NXRYRUYPEgERkAa5k5WH618eQdFULiQSYG9oWM/u3gZQj0YiI6gQDEpGZS76mxZjPD0GTX4wmdtb4eExX9PF1M3VZREQNGgMSkZlr7W4Pb1d7SACsHBeIZk62pi6JiKjBY0AiMjGdXpSbxFGTXwylwgpWMilsrGT4z8TuUCqsIbfiuAoiovpg8v9tV65cCW9vbygUCgQHByM+Pv6+7Tdt2oT27dtDoVCgU6dO2LFjh9H6zZs3Y9CgQXBxcYFEIsHx48fLbaOgoAAzZsyAi4sLHBwcMGLECGRkZNT6sRE9SGxSGnq/uwdjPj+E2RuOY8znhxD0z1/w2Pv78OGuPwztXB1sGI6IiOqRSf/H3bhxIyIjIxEVFYWEhAR06dIFYWFhyMzMrLD9wYMHMWbMGEyePBmJiYkYPnw4hg8fjqSkJEOb3Nxc9O7dG++++26l+507dy62bduGTZs2Yf/+/bh27RqefvrpOjlGosrEJqVh+roEpGkKjJbfzC1Cdn4xNidcRX6RzmT1ERE1ZhIhhKhCuzoRHByMHj16YMWKFQAAvV4PLy8vzJo1C/Pnzy/XPiIiArm5udi+fbthWc+ePREQEIBVq1YZtb148SJ8fHyQmJiIgIAAw3KNRgM3NzesX78ezzzzDADg7Nmz6NChA+Li4tCzZ88q1a7VaqFSqaDRaKBUKmv8PaDGSacX6P3unnLh6F5qpQIH5j/GZ6YREdWiqn5+m+wKUlFREY4dO4bQ0NC/i5FKERoairi4uArfExcXZ9QeAMLCwiptX5Fjx46huLjYaDvt27dHixYtqrUdoocRn5p133AEAOnaAsSnZtVbTURE9DeTddK+ceMGdDodPDw8jJZ7eHjg7NmzFb4nPT29wvbp6elV3m96ejrkcjmcnJyqtZ3CwkIUFhYaXmu12irvk6iszJz7h6PqtiMiotrFXp9VtGTJEqhUKsOXl5eXqUsiC+buqKjVdkREVLtMFpBcXV0hk8nKjR7LyMiAWq2u8D1qtbpa7SvbRlFREbKzs6u1nQULFkCj0Ri+rly5UuV9Et3ranY+bK1l8FQpUFnvIgkAT9WdIf9ERFT/TBaQ5HI5unXrht27dxuW6fV67N69GyEhIRW+JyQkxKg9AOzatavS9hXp1q0brK2tjbaTkpKCy5cv33c7NjY2UCqVRl9E1bX/j+t44uPfMG3dUcwd2Ba4G4buVfo6aqgfO2gTEZmISSeKjIyMxMSJE9G9e3cEBQVh+fLlyM3NxaRJkwAAEyZMQLNmzbBkyRIAwOzZs9G3b1988MEHGDJkCDZs2ICjR49i9erVhm1mZWXh8uXLuHbtGnA3/ODulSO1Wg2VSoXJkycjMjISzs7OUCqVmDVrFkJCQqo8go2ounR6gY93n8PHe85BCKBTMxUeaeOKmPGBiN6WbNRhW61SIGqoH8L9PU1aMxFRY2bSgBQREYHr169j0aJFSE9PR0BAAGJjYw0dsS9fvgyp9O+LXL169cL69euxcOFCvPbaa/D19cWWLVvg7+9vaLN161ZDwAKA0aNHAwCioqKwePFiAMC//vUvSKVSjBgxAoWFhQgLC8Onn35aj0dOjUlWbhFmb0jEb+duAADGBrfAoif8oLCWoZmTLQb6qcvNpM0rR0REpmXSeZAsGedBoqpIvHwLM75OwDVNARTWUrzzVCc8Hdjc1GURETVaVf385rPYiOrQmgMXcU1TAB9Xe8SMD0R7NcM0EZElYEAiqkP/fMofrg42mDvQF44Ka1OXQ0REVcR5kIhq0fnMHCzdeRald64dFdZYNNSP4YiIyMLwChJRLfnh+FUs2HwKeUU6eDnbYlxwS1OXRERENcSARPSQCkt0eHv7GXx16BIA4JE2LgjrWPXJS4mIyPwwIBE9hL9u5WHG+kScuHJnZvZZj7XBnNC2HKZPRGThGJCIauj3czcw85sEZOcVQ2VrjeURAejf3t3UZRERUS1gQCKqIbmVFDkFJejcXIWVYwPh5Wxn6pKIiKiWMCAR3YdOL4xmue7WsgnkVncGfwb5OOO/k4LQw6cJbKxkpi6ViIhqEQMSUSVik9LKPSdNJpXg9cEd8HxvHwBAb19XE1ZIRER1hfMgEVUgNikN09clGIUj3L2i9Ob2ZMQmpZmsNiIiqnsMSERl6PQC0duSUdlDCiUAorclQ6fnYwyJiBoqBiSiMuJTs8pdObqXAJCmKUB8ala91kVERPWHAYmojMycysNRTdoREZHlYUAiKsPdUVGr7YiIyPIwIBHdlZVbhKISPYJ8nOGpUqCyubAlADxVCgT5ONdzhUREVF8YkIgAJF6+hSEf/4Z3dpyBTCpB1FA/4G4Yulfp66ihfnycCBFRA8aARI2aEAJfH76EiM8OIU1TgF/PXUdOQTHC/T0RMz4QapXxbTS1SoGY8YEI9/c0Wc1ERFT3OFEkNVoFxTos3JKEb4/9BQAI76jGeyM7w1Fhfee1vycG+qmNZtIO8nHmlSMiokaAAYkapStZeZj+9TEkXdVCKgFeDW+PaY+2gkRiHH5kUglCWruYrE4iIjINBiRqdIp1eoz99yFcycqHs70cK8Z0Ra82fGQIERH9jX2QqNGxlknx+mA/BHg5Yfus3gxHRERUDq8gUaOgyS/G5Zt56NRcBQAI91djkJ8HpOxPREREFeAVJGrwzqRpMWzF75i4Jh7XsvMNyxmOiIioMgxI1KBtSbyKpz49gEs382BrLYMmv9jUJRERkQXgLTZqkIpK9HhnxxmsPXgRANDH1xUfj+6KJvZyU5dGREQWgAGJGpwMbQFe/DoBxy7dAgDMeqwN5oS25fxFRERUZQxIZJF0elHpBI4x+y7g2KVbcLSxwocRARjo52HqcomIyMIwIJHFiU1KQ/S2ZKRpCgzLPFUKRA31Q7i/J+aFt8etvCLMCW0LH1d7k9ZKRESWSSKEEKYuwhJptVqoVCpoNBoolUpTl9NoxCalYfq6BFT0QysB+Jw0IiK6r6p+fnMUG1kMnV4geltyheEIAASA6G3J0OmZ+YmI6OEwIJHFiE/NMrqtVpE0TQHiU7PqrSYiImqYGJDIYmTm3D8cVbcdERFRZRiQyGK4OypqtR0REVFlGJDIbF3Lzsf7P6XgTJoWABDk4wwXh8onepTcHc0W5ONcj1USEVFDxGH+ZFaEEIj78ya+PHgJu85kQKcXuHG7EEtHdIZMKsE/h/tj+rqEO23veV/pFJBRQ/04ISQRET00BiQyC7cLS/B94lV8efAizmXeNizv2coZAzr8PdFjuL8nYsYHlpsHSX3PPEhEREQPiwGJTE4IgSdX/I4L13MBAHZyGZ4ObIZne3qjndqxXPtwf08M9FNXOpM2ERHRw2JAolp3v8eAlK7f/0cm+rZ1h0wqgUQiwdAuTbH1+DU8G9ISI7o1h1Jhfd99yKQShLR2qYejISKixogzadcQZ9Ku2P0eAxLk44KNR65g3aFLuJqdj9XPdsOgjmoAQEGxDnKZFFJeBSIiojpU1c9vXkGiWlPZY0DSNAV4YV0CrKQSlNyd5drJzhrZecWGNgprWT1XS0REVDkGJKoVD3oMCACU6AU6NnXExF4+GNalKUMRERGZLQYkqhVVeQwIACwc4oeQ1q71UhMREVFNcaJIqhVVfwxIYZ3XQkRE9LAYkKhW8DEgRETUkDAgUa3o4d0EdvLK+xTxMSBERGRJGJDooQkhsHTnWeQV6Spcz8eAEBGRpWFAoof2r1/O4d+/pwIAnu3ZEp4q49toapUCMeMD+RgQIiKyGBzFRg9l1f4L+Hj3OQDA4qF+eO4RHywe1pGPASEiIovGgEQ1lnojF+/9lAIAeDW8HZ57xAfgY0CIiKgBYECiGvNxtcfKsYE4m67Fi/3amLocIiKiWsOARNVWrNPDWnan+1q4vxrh/mpTl0RERFSr2EmbquWX5AyEL/8VV7LyTF0KERFRnWFAoir77dx1vPh1Ai5cz8V/D140dTlERER1hgGJqiQ+NQtTvjyKIp0e4R3VmP94e1OXREREVGcYkOiBjl/JxvNrj6CgWI9+7dzw8ZiusJLxR4eIiBoufsrRfZ1J02LiF/G4XViCnq2csWp8N8it+GNDREQNm1l80q1cuRLe3t5QKBQIDg5GfHz8fdtv2rQJ7du3h0KhQKdOnbBjxw6j9UIILFq0CJ6enrC1tUVoaCjOnTtn1Mbb2xsSicToa+nSpXVyfJZKCIHFW09Dk1+Mri2c8O+JPaCwrvx5a0RERA2FyQPSxo0bERkZiaioKCQkJKBLly4ICwtDZmZmhe0PHjyIMWPGYPLkyUhMTMTw4cMxfPhwJCUlGdosW7YMH3/8MVatWoXDhw/D3t4eYWFhKCgoMNrWm2++ibS0NMPXrFmz6vx4LYlEIsHKcYF4umszrJ0UBAcbzgpBRESNg0QIIUxZQHBwMHr06IEVK1YAAPR6Pby8vDBr1izMnz+/XPuIiAjk5uZi+/bthmU9e/ZEQEAAVq1aBSEEmjZtipdffhmvvPIKAECj0cDDwwNr167F6NGjgbtXkObMmYM5c+bUqG6tVguVSgWNRgOlUlnDozdPRSV63kYjIqIGqaqf3yb9FCwqKsKxY8cQGhr6d0FSKUJDQxEXF1fhe+Li4ozaA0BYWJihfWpqKtLT043aqFQqBAcHl9vm0qVL4eLigq5du+K9995DSUlJpbUWFhZCq9UafTVE13MKMeTj3/C/I1dMXQoREZHJmDQg3bhxAzqdDh4eHkbLPTw8kJ6eXuF70tPT79u+9M8HbfOll17Chg0bsHfvXkybNg3vvPMOXn311UprXbJkCVQqleHLy8urBkds3rLzivDsfw7jXOZtfLT7HPKKKg+MREREDVmj7VQSGRlp+Hvnzp0hl8sxbdo0LFmyBDY2NuXaL1iwwOg9Wq22QYWknIJiTPwiHmfTc+DmaIN1/xcMO3mj/fEgIqJGzqRXkFxdXSGTyZCRkWG0PCMjA2p1xc/3UqvV921f+md1tom7faFKSkpw8WLFM0Tb2NhAqVQafTUUeUUlmLz2KE78pUETO2t8/X/B8HG1N3VZREREJmPSgCSXy9GtWzfs3r3bsEyv12P37t0ICQmp8D0hISFG7QFg165dhvY+Pj5Qq9VGbbRaLQ4fPlzpNgHg+PHjkEqlcHd3r4UjM186vUDchZv44fhVxF24ibyiEkz76hjiL2bBUWGFryYHo62Ho6nLJCIiMimT30OJjIzExIkT0b17dwQFBWH58uXIzc3FpEmTAAATJkxAs2bNsGTJEgDA7Nmz0bdvX3zwwQcYMmQINmzYgKNHj2L16tXA3aHpc+bMwdtvvw1fX1/4+PjgjTfeQNOmTTF8+HDgbkfvw4cPo3///nB0dERcXBzmzp2L8ePHo0mTJib8btSt2KQ0RG9LRprm7+kOVLbW0OQXw04uw9pJPeDfTGXSGomIiMyByQNSREQErl+/jkWLFiE9PR0BAQGIjY01dLK+fPkypNK/L3T16tUL69evx8KFC/Haa6/B19cXW7Zsgb+/v6HNq6++itzcXEydOhXZ2dno3bs3YmNjoVAogLu3yzZs2IDFixejsLAQPj4+mDt3rlEfo4YmNikN09cloOycDtr8YgDA1D6t0K2ls0lqIyIiMjcmnwfJUlnSPEg6vUDvd/cYXTm6lwSAWqXA7/Meg0wqqff6iIiI6otFzINE9SM+NavScAQAAkCapgDxqVn1WhcREZG5YkBqBDJzKg9HNWlHRETU0DEgNQJXb+VVqZ27o6LOayEiIrIEJu+kTXVHW1CMuRuOY/fZih/8W6q0D1KQDztpExERgVeQGjYHuRWy8opgJZVgoJ8HJHfD0L1KX0cN9WMHbSIiort4BakBEULg5+QMPNLGFQ42VpBKJVg2ojMkEgnauDtUOA+SWqVA1FA/hPt7mrR2IiIic8KA1EBcuH4bi7eexm/nbuCFvq0x//H2AADfe2bFDvf3xEA/NeJTs5CZUwB3xzu31XjliIiIyBgDkoW7XViCT/acwxe/p6JYJyC3ksJOLqu0vUwqQUhrl3qtkYiIyNIwIFkoIQS2nUzDP39MRoa2EAAQ2sEdbzzhh5YufNAsERHRw2BAslAr9pzHB7v+AAC0dLFD1FA/PNbew9RlERERNQgMSBbqme7NsfbgRUx6xBv/16cVFNaV31YjIiKi6mFAMiM6vaiwA7VeL/Bdwl84dVWDN5+881BeT5UtDsx/jMGIiIioDjAgmYmKhuB7qhSY1MsbsafTkXA5GwAwpJMnglvd6WTNcERERFQ3GJDMQGxSGqavS4AoszxNU4B3dp4FANjLZXhpgC+6tmhikhqJiIgaEwYkE9PpBaK3JZcLR/eytZbi57l90ayJbT1WRkRE1HjxUSMmFp+aZXRbrSL5xXpczqraA2eJiIjo4TEgmVhmzv3DUXXbERER0cNjQDIxd0dFrbYjIiKih8eAZGJBPs7wVClQ2dPQJHdHswX5ONdzZURERI0XA5KJyaQSRA31A+6GoXuVvo4a6scHyhIREdUjBiQzEO7viZjxgVCrjG+jqVUKxIwPRLi/p8lqIyIiaow4zN9MhPt7YqCfusKZtImIiKh+MSCZEZlUgpDWLqYug4iIqNHjLTYiIiKiMhiQiIiIiMpgQCIiIiIqgwGJiIiIqAwGJCIiIqIyGJCIiIiIymBAIiIiIiqDAYmIiIioDAYkIiIiojI4k3YNCSEAAFqt1tSlEBERURWVfm6Xfo5XhgGphnJycgAAXl5epi6FiIiIqiknJwcqlarS9RLxoAhFFdLr9bh27RocHR0hkfCBsmVptVp4eXnhypUrUCqVpi6HeE7MDs+HeeH5MC91eT6EEMjJyUHTpk0hlVbe04hXkGpIKpWiefPmpi7D7CmVSv5nY2Z4TswLz4d54fkwL3V1Pu535agUO2kTERERlcGARERERFQGAxLVCRsbG0RFRcHGxsbUpdBdPCfmhefDvPB8mBdzOB/spE1ERERUBq8gEREREZXBgERERERUBgMSERERURkMSERERERlMCCRwa+//oqhQ4eiadOmkEgk2LJli9F6IQQWLVoET09P2NraIjQ0FOfOnTNqk5WVhXHjxkGpVMLJyQmTJ0/G7du3jdqcPHkSffr0gUKhgJeXF5YtW1aulk2bNqF9+/ZQKBTo1KkTduzYUUdHbb6WLFmCHj16wNHREe7u7hg+fDhSUlKM2hQUFGDGjBlwcXGBg4MDRowYgYyMDKM2ly9fxpAhQ2BnZwd3d3f84x//QElJiVGbffv2ITAwEDY2NmjTpg3Wrl1brp6VK1fC29sbCoUCwcHBiI+Pr6MjN08xMTHo3LmzYeK6kJAQ7Ny507Ce58K0li5dColEgjlz5hiW8ZzUn8WLF0MikRh9tW/f3rDeIs+FILprx44d4vXXXxebN28WAMT3339vtH7p0qVCpVKJLVu2iBMnTohhw4YJHx8fkZ+fb2gTHh4uunTpIg4dOiR+++030aZNGzFmzBjDeo1GIzw8PMS4ceNEUlKS+Oabb4Stra347LPPDG0OHDggZDKZWLZsmUhOThYLFy4U1tbW4tSpU/X0nTAPYWFhYs2aNSIpKUkcP35cDB48WLRo0ULcvn3b0OaFF14QXl5eYvfu3eLo0aOiZ8+eolevXob1JSUlwt/fX4SGhorExESxY8cO4erqKhYsWGBo8+effwo7OzsRGRkpkpOTxSeffCJkMpmIjY01tNmwYYOQy+Xiiy++EKdPnxZTpkwRTk5OIiMjox6/I6a1detW8eOPP4o//vhDpKSkiNdee01YW1uLpKQkIXguTCo+Pl54e3uLzp07i9mzZxuW85zUn6ioKNGxY0eRlpZm+Lp+/bphvSWeCwYkqlDZgKTX64VarRbvvfeeYVl2drawsbER33zzjRBCiOTkZAFAHDlyxNBm586dQiKRiKtXrwohhPj0009FkyZNRGFhoaHNvHnzRLt27QyvR40aJYYMGWJUT3BwsJg2bVodHa1lyMzMFADE/v37hbj7/be2thabNm0ytDlz5owAIOLi4oS4G3qlUqlIT083tImJiRFKpdJwDl599VXRsWNHo31FRESIsLAww+ugoCAxY8YMw2udTieaNm0qlixZUodHbP6aNGki/v3vf/NcmFBOTo7w9fUVu3btEn379jUEJJ6T+hUVFSW6dOlS4TpLPRe8xUZVkpqaivT0dISGhhqWqVQqBAcHIy4uDgAQFxcHJycndO/e3dAmNDQUUqkUhw8fNrR59NFHIZfLDW3CwsKQkpKCW7duGdrcu5/SNqX7aaw0Gg0AwNnZGQBw7NgxFBcXG32v2rdvjxYtWhidk06dOsHDw8PQJiwsDFqtFqdPnza0ud/3u6ioCMeOHTNqI5VKERoa2mjPiU6nw4YNG5Cbm4uQkBCeCxOaMWMGhgwZUu77xnNS/86dO4emTZuiVatWGDduHC5fvgxY8LlgQKIqSU9PBwCjH97S16Xr0tPT4e7ubrTeysoKzs7ORm0q2sa9+6isTen6xkiv12POnDl45JFH4O/vD9z9Psnlcjg5ORm1LXtOavr91mq1yM/Px40bN6DT6XhOAJw6dQoODg6wsbHBCy+8gO+//x5+fn48FyayYcMGJCQkYMmSJeXW8ZzUr+DgYKxduxaxsbGIiYlBamoq+vTpg5ycHIs9F1bVfgcR1bsZM2YgKSkJv//+u6lLadTatWuH48ePQ6PR4Ntvv8XEiROxf/9+U5fVKF25cgWzZ8/Grl27oFAoTF1Oo/f4448b/t65c2cEBwejZcuW+N///gdbW1uT1lZTvIJEVaJWqwGg3KiDjIwMwzq1Wo3MzEyj9SUlJcjKyjJqU9E27t1HZW1K1zc2M2fOxPbt27F37140b97csFytVqOoqAjZ2dlG7cuek5p+v5VKJWxtbeHq6gqZTMZzAkAul6NNmzbo1q0blixZgi5duuCjjz7iuTCBY8eOITMzE4GBgbCysoKVlRX279+Pjz/+GFZWVvDw8OA5MSEnJye0bdsW58+ft9h/HwxIVCU+Pj5Qq9XYvXu3YZlWq8Xhw4cREhICAAgJCUF2djaOHTtmaLNnzx7o9XoEBwcb2vz6668oLi42tNm1axfatWuHJk2aGNrcu5/SNqX7aSyEEJg5cya+//577NmzBz4+Pkbru3XrBmtra6PvVUpKCi5fvmx0Tk6dOmUUXHft2gWlUgk/Pz9Dm/t9v+VyObp162bURq/XY/fu3Y3unJSl1+tRWFjIc2ECAwYMwKlTp3D8+HHDV/fu3TFu3DjD33lOTOf27du4cOECPD09LfffR7W7dVODlZOTIxITE0ViYqIAID788EORmJgoLl26JMTdYf5OTk7ihx9+ECdPnhRPPvlkhcP8u3btKg4fPix+//134evrazTMPzs7W3h4eIhnn31WJCUliQ0bNgg7O7tyw/ytrKzE+++/L86cOSOioqIa5TD/6dOnC5VKJfbt22c0dDYvL8/Q5oUXXhAtWrQQe/bsEUePHhUhISEiJCTEsL506OygQYPE8ePHRWxsrHBzc6tw6Ow//vEPcebMGbFy5coKh87a2NiItWvXiuTkZDF16lTh5ORkNOKkoZs/f77Yv3+/SE1NFSdPnhTz588XEolE/Pzzz0LwXJiFe0exCZ6TevXyyy+Lffv2idTUVHHgwAERGhoqXF1dRWZmphAWei4YkMhg7969AkC5r4kTJwpxd6j/G2+8ITw8PISNjY0YMGCASElJMdrGzZs3xZgxY4SDg4NQKpVi0qRJIicnx6jNiRMnRO/evYWNjY1o1qyZWLp0abla/ve//4m2bdsKuVwuOnbsKH788cc6PnrzU9G5ACDWrFljaJOfny9efPFF0aRJE2FnZyeeeuopkZaWZrSdixcviscff1zY2toKV1dX8fLLL4vi4mKjNnv37hUBAQFCLpeLVq1aGe2j1CeffCJatGgh5HK5CAoKEocOHarDozc/zz//vGjZsqWQy+XCzc1NDBgwwBCOBM+FWSgbkHhO6k9ERITw9PQUcrlcNGvWTERERIjz588b1lviuZCIO/8RExEREdFd7INEREREVAYDEhEREVEZDEhEREREZTAgEREREZXBgERERERUBgMSERERURkMSERERERlMCARUaMikUiwZcsWU5dBRGaOAYmIGpTr169j+vTpaNGiBWxsbKBWqxEWFoYDBw7U2j769euHOXPm1Nr2iMj8WJm6ACKi2jRixAgUFRXhv//9L1q1aoWMjAzs3r0bN2/eNHVpRGRBeAWJiBqM7Oxs/Pbbb3j33XfRv39/tGzZEkFBQViwYAGGDRtmaHfjxg089dRTsLOzg6+vL7Zu3Wq0nf379yMoKAg2Njbw9PTE/PnzUVJSAgB47rnnsH//fnz00UeQSCSQSCS4ePEibt26hXHjxsHNzQ22trbw9fXFmjVr6v17QES1gwGJiBoMBwcHODg4YMuWLSgsLKy0XXR0NEaNGoWTJ09i8ODBGDduHLKysgAAV69exeDBg9GjRw+cOHECMTEx+M9//oO3334bAPDRRx8hJCQEU6ZMQVpaGtLS0uDl5YU33ngDycnJ2LlzJ86cOYOYmBi4urrW27ETUe3iw2qJqEH57rvvMGXKFOTn5yMwMBB9+/bF6NGj0blzZ+BuJ+2FCxfirbfeAgDk5ubCwcEBO3fuRHh4OF5//XV89913OHPmDCQSCQDg008/xbx586DRaCCVStGvXz8EBARg+fLlhv0OGzYMrq6u+OKLL0x05ERUm3gFiYgalBEjRuDatWvYunUrwsPDsW/fPgQGBmLt2rWGNqVhCQDs7e2hVCqRmZkJADhz5gxCQkIM4QgAHnnkEdy+fRt//fVXpfudPn06NmzYgICAALz66qs4ePBgnR0jEdU9BiQianAUCgUGDhyIN954AwcPHsRzzz2HqKgow3pra2uj9hKJBHq9/qH2+fjjj+PSpUuYO3curl27hgEDBuCVV155qG0SkekwIBFRg+fn54fc3Nwqte3QoQPi4uJwb++DAwcOwNHREc2bNwcAyOVy6HS6cu91c3PDxIkTsW7dOixfvhyrV6+uxaMgovrEgEREDcbNmzfx2GOPYd26dTh58iRSU1OxadMmLFu2DE8++WSVtvHiiy/iypUrmDVrFs6ePYsffvgBUVFRiIyMhFR6579Mb29vHD58GBcvXsSNGzeg1+uxaNEi/PDDDzh//jxOnz6N7du3o0OHDnV8xERUVzgPEhE1GA4ODggODsa//vUvXLhwAcXFxfDy8sKUKVPw2muvVWkbzZo1w44dO/CPf/wDXbp0gbOzMyZPnoyFCxca2rzyyiuYOHEi/Pz8kJ+fj9TUVMjlcixYsAAXL16Era0t+vTpgw0bNtTh0RJRXeIoNiIiIqIyeIuNiIiIqAwGJCIiIqIyGJCIiIiIymBAIiIiIiqDAYmIiIioDAYkIiIiojIYkIiIiIjKYEAiIiIiKoMBiYiIiKgMBiQiIiKiMhiQiIiIiMpgQCIiIiIq4/8Bvtmz7GP92iwAAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "import timeit\n", + "\n", + "import matplotlib.pyplot as plt\n", + "\n", + "# 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(), 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", + "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." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Scaling up!" + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "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_qubits, rng)\n", + "\n", + "print(f\"Hidden shift string {hidden_shift_string}\")" + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "metadata": {}, + "outputs": [], + "source": [ + "isa_circuit = get_isa_circuit(circuit, backend)" + ] + }, + { + "cell_type": "code", + "execution_count": 20, + "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, + "metadata": {}, + "outputs": [], + "source": [ + "counts, pub_result = get_bitstring_counts(job)" + ] + }, + { + "cell_type": "code", + "execution_count": 22, + "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(counts, hidden_shift_string)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We see that the correct hidden shift string is found. Furthermore, the nine next-most-probably bitstrings are wrong in only one position." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Record the most probable probability" + ] + }, + { + "cell_type": "code", + "execution_count": 23, + "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, + "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(mit, counts, qubit_mapping)" + ] + }, + { + "cell_type": "code", + "execution_count": 24, + "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(max_probability_before_M3, max_probability_after_M3, is_hidden_shift_identified)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The results show that readout error was the dominant source of error and the M3 mitigation was effective." + ] + } + ], + "metadata": { + "description": "Use the M3 readout mitigation addon with the Sampler primitive", + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.13.3" + }, + "title": "Readout error mitigation for the Sampler primitive using M3", + "toc": { + "base_numbering": 0 + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/qiskit_bot.yaml b/qiskit_bot.yaml index 7218a9cb80e..1c035202747 100644 --- a/qiskit_bot.yaml +++ b/qiskit_bot.yaml @@ -683,6 +683,8 @@ notifications: - "@annaliese-estes" "docs/tutorials/quantum-phase-estimation-qctrl": - "@dsierrasosa" + "docs/tutorials/readout-error-mitigation-sampler": + - "@jlapeyre" "learning/courses/quantum-computing-in-practice/introduction": - "@livlanes" "learning/courses/quantum-computing-in-practice/running-quantum-circuits":