samples/characterization/randomized-benchmarking/randomized-benchmarking.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Randomized Benchmarking with Q# and Python"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**NB**: This notebook requires additional packages to analyze randomized benchmarking results offline. Dependencies for this sample can be installed as a new virtual environment by using `conda`:\n",
    "\n",
    "```sh\n",
    "conda env create -f environment.yml\n",
    "conda activate randomized-benchmarking\n",
    "```\n",
    "\n",
    "$\\newcommand{\\Tr}{\\operatorname{Tr}}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Abstract ##"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This sample demonstrates how to use Q# and Python together to simulate randomized benchmarking, a common protocol for characterizing quantum systems. In particular, this sample shows how to generate random sequences within your quantum program, analyze the results offline, and how to use online statistical inference to adaptively choose RB experiments.\n",
    "\n",
    "> **NOTE**: This sample is somewhat more advanced than most, and assumes familiarity with open quantum systems as well as some statistics knowledge.\n",
    "> If you would like a refresher on these topics, please check out the [**process-tomography**](../process-tomography/README.md) and [**phase-estimation**](../phase-estimation/README.md) samples."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Preamble ##"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can enable Q# support in Python by importing the `qsharp` package."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Preparing Q# environment...\n"
     ]
    }
   ],
   "source": [
    "import qsharp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll also be using numerics in this sample, so we import NumPy as well."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we import plotting support and the QuTiP library, since these will be helpful to us in manipulating the quantum objects returned by the quantum process tomography functionality that we call later."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "import qutip as qt\n",
    "qt.settings.colorblind_safe = True"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, to analyze randomized benchmarking results offline, we'll be using the [QInfer](http://qinfer.org) library, so we go on and import it here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "ERROR:duecredit:Failed to import duecredit due to No module named 'duecredit'\n",
      "C:\\Users\\cgran\\Anaconda3\\envs\\iqs-dev\\lib\\site-packages\\IPython\\parallel.py:12: ShimWarning: The `IPython.parallel` package has been deprecated since IPython 4.0. You should import from ipyparallel instead.\n",
      "  warn(\"The `IPython.parallel` package has been deprecated since IPython 4.0. \"\n",
      "C:\\Users\\cgran\\Anaconda3\\envs\\iqs-dev\\lib\\site-packages\\qinfer-1.0-py3.8.egg\\qinfer\\parallel.py:61: UserWarning: Could not import IPython parallel. Parallelization support will be disabled.\n"
     ]
    }
   ],
   "source": [
    "import qinfer as qi"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Randomized Benchmarking ##"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider the problem of learning the noise acting on an $n$-qubit register. Without any prior assumptions, density operators on that register are matrices of $2^n \\times 2^n$ complex numbers, such that noise channels acting on those density operators can be expressed as $4^n \\times 4^n$ complex numbers. Even for a two-qubit register, it takes 256 complex numbers to write down a noise process acting on that register.\n",
    "\n",
    "This quickly gets out of hand, such that using tomography to learn noise quickly becomes infeasible as systems grow in size. Moreover, process tomography rests on the critical assumption that state preparation and measurement (also known as \"_SPAM_\") are either perfect, or at least only subject to noise that itself is fully characterized. If this assumption is violated, SPAM errors can be conflated with the noise being estimated and vice versa.\n",
    "\n",
    "One way to solve this is to use a technique known as _randomized benchmarking_, which uses a total of three real numbers to parameterize noise models instead of a tomographically complete description. By choosing sequences of gates of different lengths, randomized benchmarking can also deal with state preparation and measurement errors — errors due to gates grow with sequence lengths, while state preparation and measurement errors appear as constant offsets.\n",
    "\n",
    "> **NOTE**: Though randomized benchmarking is works well for larger systems, for simplicitly this sample considers the single-qubit case."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Sequence generation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Concretely, randomized benchmarking works by randomly choosing sequences of _Clifford operations_, those operations that map Pauli operations to other Pauli operations. In Q#, we can draw random elements of the single-qubit Clifford group by using the `DrawRandomSingleQubitClifford` operation.\n",
    "\n",
    "> **💡 TIP**: The `%%qsharp` magic command tells IPython to treat that cell as Q# and to pass it to the Q# compiler, adding any new operations as Python variables.\n",
    "> For more info, see the [Python interoperability sample](../../interoperability/python/)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%qsharp\n",
    "open Microsoft.Quantum.Arrays;\n",
    "open Microsoft.Quantum.Random;\n",
    "open Microsoft.Quantum.Synthesis;\n",
    "\n",
    "operation DrawRandomSingleQubitRBSequence(nOperators : Int) : SingleQubitClifford[] {\n",
    "    return DrawMany(DrawRandomSingleQubitClifford, nOperators, ());\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[(0, 0, 0, 0),\n",
       " (0, 2, 0, 4),\n",
       " (1, 3, 1, 0),\n",
       " (2, 0, 0, 4),\n",
       " (0, 2, 0, 6),\n",
       " (0, 3, 0, 1),\n",
       " (2, 3, 0, 3),\n",
       " (2, 0, 1, 4),\n",
       " (1, 0, 1, 7),\n",
       " (2, 3, 0, 3)]"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "DrawRandomSingleQubitRBSequence.simulate(nOperators=10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Each of the tuples in the above list represents a different Clifford operation, using the notation defined by the [`Microsoft.Quantum.Synthesis.SingleQubitClifford` user-defined type](https://docs.microsoft.com/qsharp/api/qsharp/microsoft.quantum.synthesis.singlequbitclifford).\n",
    "\n",
    "> **📝 NOTE**: In randomized benchmarking, the term \"sequence length\" generally refers to the number of distinct Clifford operations performed in a given sequence. For instance, the tuple `(0, 1, 1, 0)` refers to the operation `BoundCA([X, S])`. For more details, see the documentation for the [`Apply1C` operation](https://docs.microsoft.com/qsharp/api/qsharp/microsoft.quantum.synthesis.apply1c) and the [`BoundCA` function](https://docs.microsoft.com/qsharp/api/qsharp/microsoft.quantum.canon.boundca).\n",
    "\n",
    "In practice, however, randomized benchmarking requires that each different sequence multiplies to identity (written as the tuple `(0, 0, 0, 0)`), but this isn't true in general for the sequences chosen by the `DrawRandomSingleQubitRBSequence` that we defined above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%qsharp\n",
    "\n",
    "operation CheckThatRandomSequenceMultipliesToIdentity(nOperators : Int) : Unit {\n",
    "    let seq = DrawRandomSingleQubitRBSequence(nOperators);\n",
    "\n",
    "    // Fold uses a function representing a binary operator to reduce\n",
    "    // a Q# array to a single value. In this case, we use the FlippedCall\n",
    "    // function below to implement the opposite product 𝑎 ×ₒₚ 𝑏 ≔ 𝑏 × 𝑎.\n",
    "    // This allows us to multiply a sequence \"backwards,\" since the matrix\n",
    "    // product 𝐵𝐴 represents applying the operation A and then the operation B.\n",
    "    let product = Fold(FlippedCall(Times1C, _, _), Identity1C(), seq);\n",
    "    \n",
    "    // At this point, we check if the product was the identity (0, 0, 0, 0),\n",
    "    // and if not, then we can report the actual product.\n",
    "    EqualityFact1C(Identity1C(), product, \"Random sequence did not multiply to identity.\");\n",
    "}\n",
    "    \n",
    "function FlippedCall<'TLeft, 'TRight, 'TOutput>(\n",
    "    fn : (('TLeft, 'TRight) -> 'TOutput),\n",
    "    right : 'TRight,\n",
    "    left : 'TLeft\n",
    ")\n",
    ": 'TOutput {\n",
    "    return fn(left, right);\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "\r\n",
       "                    <details>\r\n",
       "                        <summary>\r\n",
       "                            Unhandled exception of type Microsoft.Quantum.Simulation.Core.ExecutionFailException: Random sequence did not multiply to identity.\n",
       "\tExpected:\tSingleQubitClifford((2, 0, 1, 1))\n",
       "\tActual:\tSingleQubitClifford((0, 0, 0, 0))\r\n",
       "                        </summary>\r\n",
       "                        <table>\r\n",
       "                            \r\n",
       "                    <thead>\r\n",
       "                        <tr>\r\n",
       "                            <th>Source</th>\r\n",
       "                            <th>Callable</th>\r\n",
       "                        </tr>\r\n",
       "                    </thead>\r\n",
       "\r\n",
       "                    <tbody>\r\n",
       "                        \r\n",
       "                        <tr>\r\n",
       "                            <td>(notebook)</td>\r\n",
       "                            <td>CheckThatRandomSequenceMultipliesToIdentity</td>\r\n",
       "                        </tr>\r\n",
       "                    \r\n",
       "                    </tbody>\r\n",
       "                \r\n",
       "                        </table>\r\n",
       "                    </details>\r\n",
       "                "
      ],
      "text/plain": [
       "Unhandled exception. Microsoft.Quantum.Simulation.Core.ExecutionFailException: Random sequence did not multiply to identity.\n",
       "\tExpected:\tSingleQubitClifford((2, 0, 1, 1))\n",
       "\tActual:\tSingleQubitClifford((0, 0, 0, 0))\r\n",
       " ---> SNIPPET.CheckThatRandomSequenceMultipliesToIdentity on C:\\snippet_.qs:line 0\r\n"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Random sequence did not multiply to identity.\n",
      "\tExpected:\tSingleQubitClifford((2, 0, 1, 1))\n",
      "\tActual:\tSingleQubitClifford((0, 0, 0, 0))\r\n"
     ]
    }
   ],
   "source": [
    "CheckThatRandomSequenceMultipliesToIdentity.simulate(nOperators=10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can modify the definition of DrawRandomSingleQubitRBSequence to correctly generate identity sequences by drawing one less Clifford operator and then appending the inverse of the all the other operations:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%qsharp\n",
    "\n",
    "operation DrawRandomSingleQubitRBSequence(nOperators : Int) : SingleQubitClifford[] {\n",
    "    let root = DrawMany(DrawRandomSingleQubitClifford, nOperators - 1, ());\n",
    "    let product = Fold(FlippedCall(Times1C, _, _), Identity1C(), root);\n",
    "    return root + [Inverse1C(product)];\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With this new definition in place, we can verify that our random sequence generation returns identity sequences."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "()"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "CheckThatRandomSequenceMultipliesToIdentity.simulate(nOperators=10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### RB measurements"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have an operation that can randomly generate sequences that do nothing on an ideal quantum device, the core insight behind randomized benchmarking is that if we run random no-op programs on a noisy device, it's easy to see when we have an error. That is, if we prepare $\\left|{0}\\right\\rangle$ and observe anything other than `Zero` when measuring in the $Z$-basis at the end, we can conclude that we saw at least one error."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%qsharp\n",
    "open Microsoft.Quantum.Measurement;\n",
    "\n",
    "operation MeasureSingleQubitRBSequence(nOperators : Int) : Result {\n",
    "    use q = Qubit();\n",
    "    let ops = DrawRandomSingleQubitRBSequence(nOperators);\n",
    "    for op in ops {\n",
    "        // The Apply1C operation takes a single-qubit Clifford operator\n",
    "        // and applies it as an operation.\n",
    "        Apply1C(op, q);\n",
    "    }\n",
    "    return MResetZ(q);\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Indeed, when we simulate how this protocol would work on an ideal device using the full-state simulator, we see that we always get `Zero` back."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sum(MeasureSingleQubitRBSequence.simulate(nOperators=10) for _ in range(50))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In general, randomized benchmarking requires us to draw a new sequence for each measurement [arXiv:1404.5275](https://arxiv.org/abs/1404.5275); though some generalizations exist that relax this limitation ([arXiv:1802.00401](https://arxiv.org/abs/1802.00401)), they require more sophisticated statistical analysis such that we focus on randomized benchmarking as originally presented here.\n",
    "\n",
    "To that end, to collect multiple \"shots\" of RB, we can't pick a sequence then measure it many times, but must include a call to our sequence generation operation in our loop. Doing so, we can estimate the probability with which sequences of a given length correctly produce a result of `Zero` at the end. This probability is often called the \"_survival probability_.\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%qsharp\n",
    "\n",
    "operation MeasureSingleQubitRBSequenceLength(nOperators : Int, nSequences : Int) : Int {\n",
    "    mutable nUp = 0;\n",
    "    for idxSequence in 1..nSequences {\n",
    "        set nUp += MeasureSingleQubitRBSequence(nOperators) == Zero ? 1 | 0;\n",
    "    }\n",
    "    return nUp;\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1.0"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "MeasureSingleQubitRBSequenceLength.simulate(nOperators=10, nSequences=100) / 100"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As expected, the survival probability for a system without noise is 100%."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Noise models for RB"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Of course, we want to use randomized benchmarking to study how imperfect devices are affected by noise, so let's switch from using the full-state simulator to the open systems simulator provided with the Quantum Development Kit. To do so, we'll need to provide a _noise model_ that describes the noise we want to simulate.\n",
    "\n",
    "> **TIP**: For a refresher on open systems concepts, including channels and density operators, check out https://github.com/microsoft/qsharp-runtime/blob/main/documentation/examples/open-systems-concepts.ipynb."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We start by asking the `qsharp` package for an ideal noise model, representing the case of a noiseless device. Then, we can fill in each part of the noise model with channels that we want to use to simulate each different operation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "noise_model = qsharp.get_noise_model_by_name('ideal')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using the [QuTiP library](https://qutip.org/), we can randomly choose channels from the space of all possible channels with [the `rand_super_bcsz` function](https://qutip.org/docs/latest/apidoc/functions.html#qutip.random_objects.rand_super_bcsz)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "Quantum object: dims = [[[2], [2]], [[2], [2]]], shape = (4, 4), type = super, isherm = False\\begin{equation*}\\left(\\begin{array}{*{11}c}0.609 & (-0.006+0.003j) & (-0.006-0.003j) & 0.491\\\\(0.185+0.308j) & (0.032-0.170j) & (0.035-0.213j) & (0.221-0.074j)\\\\(0.185-0.308j) & (0.035+0.213j) & (0.032+0.170j) & (0.221+0.074j)\\\\0.391 & (0.006-0.003j) & (0.006+0.003j) & 0.509\\\\\\end{array}\\right)\\end{equation*}"
      ],
      "text/plain": [
       "Quantum object: dims = [[[2], [2]], [[2], [2]]], shape = (4, 4), type = super, isherm = False\n",
       "Qobj data =\n",
       "[[ 0.60897679+1.84314369e-18j -0.00557146+3.19611640e-03j\n",
       "  -0.00557146-3.19611640e-03j  0.49051966+1.08962318e-16j]\n",
       " [ 0.18490256+3.08259005e-01j  0.03182878-1.70194656e-01j\n",
       "   0.03455891-2.13070916e-01j  0.22083071-7.41498439e-02j]\n",
       " [ 0.18490256-3.08259005e-01j  0.03455891+2.13070916e-01j\n",
       "   0.03182878+1.70194656e-01j  0.22083071+7.41498439e-02j]\n",
       " [ 0.39102321+1.30104261e-18j  0.00557146-3.19611640e-03j\n",
       "   0.00557146+3.19611640e-03j  0.50948034+1.13515967e-16j]]"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "qt.rand_super_bcsz()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In general, these channels will have very poor fidelity, so to get a more realistic noise model we can mix channels drawn from the `rand_super_bcsz` function with channels describing the ideal action of each operation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.41912961921549163"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "qt.average_gate_fidelity(qt.rand_super_bcsz())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "def mix_with_bcsz(ideal, pr_success=0.99):\n",
    "    return pr_success * qt.to_super(ideal) + (1 - pr_success) * qt.rand_super_bcsz()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.9942763044478301"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "qt.average_gate_fidelity(mix_with_bcsz(qt.qeye(2)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's use this to fill out our noise model, making sure to define what happens for each Pauli and Clifford operation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "noise_model.x = mix_with_bcsz(qt.sigmax())\n",
    "noise_model.y = mix_with_bcsz(qt.sigmay())\n",
    "noise_model.z = mix_with_bcsz(qt.sigmaz())\n",
    "noise_model.h = mix_with_bcsz(qt.qip.operations.hadamard_transform())\n",
    "noise_model.s = mix_with_bcsz(qt.qip.operations.phasegate(np.pi / 2))\n",
    "noise_model.s_adj = mix_with_bcsz(qt.qip.operations.phasegate(-np.pi / 2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "qsharp.set_noise_model(noise_model)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Running our RB measurements using this noise model, we note that we no longer get that every sequence at a given length returns `Zero`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "85"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "MeasureSingleQubitRBSequenceLength.simulate_noise(nOperators=10, nSequences=100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Indeed, if we plot the number of sequence measurements that correctly return `Zero` at each length, we can confirm that as we run longer and longer sequences, the effects of the noise increase."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([  3,   3,   4,   5,   6,   8,  10,  13,  16,  19,  24,  30,  37,\n",
       "        45,  56,  69,  86, 106, 131, 162, 200])"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "seq_lens = np.logspace(np.log10(3), np.log10(200), 21).astype(int)\n",
    "seq_lens"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> **📝 NOTE**: We use exponentially spaced sequences here since longer sequences tell us more about the fidelity (that is, longer sequences depend more strongly on $p$) than shorter sequences. This allows us to both use short sequences to learn $A$ and $B$, as well as longer sequences to learn $p$. We'll see a better approach to choosing sequence lengths later when we see how to include adaptivity in our randomized benchmarking program."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "est_pr_survival = [\n",
    "    MeasureSingleQubitRBSequenceLength.simulate_noise(nOperators=seq_len, nSequences=200) / 200\n",
    "    for seq_len in seq_lens\n",
    "]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Estimated survival probability')"
      ]
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEKCAYAAAAfGVI8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAejElEQVR4nO3de5gcdZ3v8fcnAeSacElUzGVJBoQnaOQSIcYLiIoBBQ7CuoCuynqMsERkedwD6HG9u+JRVhEkRg6KIIIc8IAIAgpEjwEkETIQMJAEhBE0BDQBUSDwPX9UdazpdPdUJl3VPV2f1/PMM91V1VXfqZnpb9fv96vvTxGBmZlV16hOB2BmZp3lRGBmVnFOBGZmFedEYGZWcU4EZmYV50RgZlZxhSUCSRdIWiXpnibrJelsScsl9Uvap6hYzMysuSKvCL4LzG6x/hBgt/RrDnBegbGYmVkThSWCiPgF8GSLTY4AvheJ24DtJe1cVDxmZtbYZh089gTgkczzgXTZY/UbSppDctXANttss+8ee+xRSoBmZr1i8eLFqyNifKN1nUwEarCsYb2LiJgPzAeYMWNGLFq0qMi4zMx6jqTfNVvXyVFDA8CkzPOJwKPtPsi8BStYuGL1oGULV6xm3oIV7T6UmdmI1MlEcDXwvnT00ExgTURs0Cy0qaZPHMvcS+5cnwwWrljN3EvuZPrEse0+lJnZiFRY05CkHwAHAuMkDQCfAjYHiIh5wLXAocBy4Bng+CLimNU3jnOO25u5l9zJe/efzMW3P8w5x+3NrL5xRRzOzGzEKSwRRMSxQ6wP4KSijp81q28c791/MmfftJyTD9rVScDMLKMSdxYvXLGai29/mJMP2pWLb394gz4DM7Mq6/lEUOsTOOe4vTn14N3XNxM5GZiZJTo5fLQU/QNrOOe4vekfWAP8vc+g9rx/YA0nHNDXyRDNzDqq568ITjigj1l94waNHuofWMPoUQwaPeQhpWZWVT2fCGqyo4eW/WEtX/zJbznxwKnM6hvnIaVmVmk93zQ0b8EKpk8cy6y+cYNGD+35ijGcd8tKnvrrOg8pNbNK6/krgmyT0MIVq/nOwofYcvNRPPzkMxzwyvGcfdNy3rv/ZCcBM6usnr8iqDUJffiixTz/wotsPnoUF3zgtSx9dA1f/MlvOXLvCVx8+8PM7NvJycDMKqnnrwggSQbTJ4zlb8+/yPGzdgHgvFtW8vF37MHuL9/OQ0rNrNJ6/ooAkhFB9/3hqfU3lD3+9LMb9AnUhpT6qsDMqqbnE0H2hrJZfeN4/Olnuab/MQ57zSsGbeP7Ccysqnq+aah2Q1ntk34tAfx4SVLx2kNHzazqev6KoP5T/qy+cXzrn/dl7iV3Mn7bZR46amaV1/NXBI1k7yfw0FEzq7pKJgJXIzUz+7vKJQJXIzUzG6xyiaC+87h/YA0nHjh1fTVScAE6M6uWyiWCWjXSmukTx3LeLSsHVSH1KCIzq5KeHzU0FM9pbGZVV7krgkY8isjMqsyJAI8iMrNqq3wiyI4i2volm3HigVMHjSJyx7GZ9brKJ4LsKKJax3FtFJE7js2sChQRnY5ho8yYMSMWLVpU2P5rb/7uODazXiJpcUTMaLSu8lcE9dxxbGZV40RQxx3HZlY1TgQZLj9hZlXkRJBRX36idrNZtvyEmVmvcWexmVkFuLPYzMyaGjIRSLpC0jskVTJpzFuwYoM+At9kZma9JM+b+3nAccADkr4kaY+CY+oq0yeO3eBOY99kZma9ZMjqoxHxM+BnksYCxwI3SnoE+DZwcUQ8X3CMHeXqpGbW63I190jaCfgA8N+BO4GvA/sANxYWWRfxTWZm1svy9BFcCfwS2Bo4LCIOj4jLIuIjwLZFB9gNfJOZmfWyPFcE50fEtIj4z4h4DEDSSwCaDUWqkTRb0jJJyyWd3mD9WEk/lrRE0lJJxw/rpyiQbzIzs16XJxF8vsGyW4d6kaTRwLnAIcA04FhJ0+o2Owm4NyJeAxwIfFXSFjliKo1vMjOzXte0s1jSy4EJwFaS9gaUrhpD0kw0lP2A5RGxMt3fpcARwL2ZbQLYTpJImpmeBNZt7A9RpBMO6Ntg2ay+ce4nMLOe0WrU0NtJOognAmdllj8FfDzHvicAj2SeDwD7121zDnA18CiwHfBPEfFi/Y4kzQHmAEyePDnHoc3MLK+miSAiLgQulHRURFwxjH2rwbL6ehZvB+4CDgL6SIam/jIi1tbFMh+YD0mJiWHEYmZmTTTtI5D03vThLpJOrf/Kse8BYFLm+USST/5ZxwNXRmI58CDQ9Tes+W5jM+slrTqLt0m/b0vSbFP/NZQ7gN0kTUk7gI8haQbKehh4C4CklwG7AytzR98hvtvYzHpJodVHJR0KfA0YDVwQEV+QdAJARMyT9Argu8DOJE1JX4qIi1vts1uqj3pKSzMbSVpVH201aujsVjuNiJOHOnBEXAtcW7dsXubxo8DBQ+2nG2XvNj75oF2dBMxsxGo1amhxaVGMQPV3G8/s22mTksG8BSuYPnHsoH0sXLGa/oE1DYewmpm1y1CjhqyB7N3Gs/rGMbNvp0HPh6PW71DbR/YYZmZFatpHIOlrEXGKpB+z4bBPIuLwooNrpBv6CIr69O5+BzMryrD6CICL0u9faX9II1v9m30tMWSXDycxuN/BzDqh6fDRiFicfl9AUlvoTyQlIG5Nl1mqXcNJXeXUzDphyIlpJL0DmAesIBniOUXShyPiuqKDGynaMXlNEf0OZmZ55Kk++lXgzRFxYEQcALwZ+K9iwxp5NnXyGlc5NbNOGfKKAFiVln+oWQmsKiieEWtTh5O6yqmZdUqrG8relT5cKula4Icko4f+kaR8hKXcrGNmI1mrK4LDMo//CByQPn4c2KGwiEagVs06TgRm1u0KrTVUhG64j8DMbKQZ7n0EtRdvCXwQ2BPYsrY8Iv6lbRGamVnH5Bk1dBHwcpJJZBaQzCvwVJFBmZlZefIkgl0j4pPAX9L6Q+8AXl1sWGZmVpY8ieD59PufJb0KGAvsUlhEZmZWqjyJYL6kHYBPkswwdi9wZqFR9RhPbWlm3WzIRBAR50fEnyJiQURMjYiXRsS3ygiuV3hqSzPrZnlGDe0EfBp4PckNZb8EPhcRTxQbWu9oRy0iM7Oi5GkaupSkpMRRwNHAauCyIoPqRZtai8jMrCh5EsGOEfG5iHgw/fo8sH3BcfUcl5g2s26VJxHcLOkYSaPSr3cDPyk6sF6SrUV06sG7r28mcjIws27QaqrKp0j6BARsA7yYrhoFPB0RY0qJsM5ILDHhienNrNNalZhwrSEzswrYpFpD6Q4OB96UPr0lIq5pV3BmZtZZQ/YRSPoS8FGSG8nuBT6aLjMzsx6Q54rgUGCviHgRQNKFwJ3A6UUGZmZm5cgzaggGDxf17bBmZj0kzxXBF4E7Jd1MMoLoTcAZhUZlZmalaZkIJI0iGTY6E3gtSSI4LSL+UEJsZmZWgpaJICJelDQ3In5IUnnUzMx6TJ4+ghslfUzSJEk71r4Kj8zMzEqRp4+gNjfxSZllAUxtfzhmZla2IRNBREwpIxDLzyUrzKyd8txQtqWkUyVdKekKSadI2rKM4KyxjZ3oxjOkmVkrefoIvgfsCXwDOAeYBlyUZ+eSZktaJmm5pIY3oEk6UNJdkpZKWpA38CrLTnRz1g3L1lc2bTbHgWdIM7NW8vQR7B4Rr8k8v1nSkqFeJGk0cC7wNmAAuEPS1RFxb2ab7YFvArMj4mFJL92o6CssO9HNyQft2nKiG8+QZmat5LkiuFPSzNoTSfsDv8rxuv2A5RGxMiKeI5np7Ii6bY4DroyIhwEiYlW+sG1jJ7rxDGlm1kyeRLA/sFDSQ5IeAm4FDpB0t6T+Fq+bADySeT6QLst6JbCDpFskLZb0vkY7kjRH0iJJix5//PEcIfe24Ux04xnSzKyZPE1Ds4e5bzVYVj/5wWbAvsBbgK2AWyXdFhH3D3pRxHxgPiTzEQwznp7RP7BmUNNOremnf2BNw0/62cQxq28cM/t2GrJfwcyqI8/w0d8Nc98DwKTM84nAow22WR0RfwH+IukXwGuA+7GmGg0RndU3rumb+sYmDjOrllwT0wzTHcBukqYAvweOIekTyLoKOEfSZsAWJM1Q/1VgTJW0sYnDzKqlsEQQEeskzQWuB0YDF0TEUkknpOvnRcR9kn4K9JMUtzs/Iu4pKiYzM9uQ5yw2M6uAYc1ZLOkpNuzchaQTOCJiTJviMzOzDmqaCCJiuzIDMTOzzsjdR5De9bu+xlDtJjAzMxvZ8hSdO1zSA8CDwALgIeC6guMyM7OS5Lmz+HMkU1Xen5akfgv5SkxYD3EFU7PelScRPB8RTwCjJI2KiJuBvYoNy7qNK5ia9a48fQR/lrQt8Avg+5JWAeuKDcu6jSuYmvWuPFcERwDPAP8G/BRYARxWZFDWnVzB1LqFmyrbK08imAO8IiLWRcSFEXF22lRkFdPOCqb+R7ZN4abK9sqTCMYA10v6paSTJL2s6KCs+wyn9HUr/ke2TbGxs/RZa0Mmgoj4TETsCZwEvAJYIOlnhUdmXaVVBdPhqMo/sq98iuOmyvbJc0VQswr4A/AE4CklK+aEA/o2+Eeb1TeuYWXTvKrwj+wrn+J4sqX2yXND2YmSbgF+DowDPhQR04sOzHpfWf/InfxUXpUrn7K1u6my6vJcEfwDcEpE7BkRn8pOPm82XGX+I3f6U3kVrnzK1u6mym5WxgeZpmWoJY2JiLWSdmy0PiKebFsUG8FlqHvDvAUrmD5x7KA3xYUrVtM/sGaTmpuaqb35d+IeiE4e20a++qlm65/n1aoMdatEcE1EvFPSgyTlqLNzEEdETN2In6VtnAhsuM66YRln37Sckw/alVMP3r2UY7brn7jdyk7Etmna8WGiVSJo2jQUEe9Mv0+JiKnp99pXR5KA2XB1qmOxW5swOt1cZhun6ObFIWcok3QVcClwVUQ809ajD4OvCGxjdeun8k5zk9XI0bErgoyzgDcC90m6XNLRkrYc6kVm3aJbP5V3mjuxR4YyBlbknrNY0mjgIOBDwOxOTVXpKwKz9vAVwcjQrv6cYc1ZXLeDrUgKzf0TsA9wYe6jm1nXqW8em9m3k5vLulSjN/tZfePa+nvKc0PZZcB9JFcD5wJ9EfGRtkVgZqVzc5lltbwikDQKuBs4LiJeKCckMytaGZ8ybeRoeUUQES8C73ASMDPrXXlGDd0g6ShJGnpTMzMbafJ0Fp8KbAOsk/Q3kjuMo1OjhszMrL2GTAQRsV0ZgZiZWWcMmQgkvanR8oj4RfvDMTOzsuVpGvr3zOMtgf2AxSTDSc3MbITL0zR0WPa5pEnAlwuLyMzMSrUxU1XWDACvancgZmbWGXn6CL5BMh8BJIljL2BJgTGZmVmJ8vQRZCu8rQN+EBG/KigeMzMrWZ4+gvUF5iTtAEwqNCIzMytVnqJzt0gak85dvAT4jqSzig/NzMzKkKezeGxErAXeBXwnIvYF3ppn55JmS1omabmk01ts91pJL0g6Ol/YZmbWLnkSwWaSdgbeDVyTd8fpRDbnAocA04BjJU1rst2ZwPV5921mZu2TJxF8luRNenlE3CFpKvBAjtftl75mZUQ8RzLv8RENtvsIcAWwKmfMZmbWRnk6iy8HLs88XwkclWPfE4BHMs8HgP2zG0iaABxJcpfya5vtSNIcYA7A5MmTcxzazMzyGs4NZXk1KltdP0Hy14DThprvICLmR8SMiJgxfvz4dsVnZmbknLN4mAYYPNR0IvBo3TYzgEvTqQ7GAYdKWhcR/7fAuMzMLKPIRHAHsJukKcDvgWOA47IbRMSU2mNJ3wWucRIwMytX00Qg6dRWL4yIlvcSRMQ6SXNJOppHAxdExFJJJ6Tr5w0jXjMza7NWVwS1CWl2J+nIvTp9fhiQay6CiLgWuLZuWcMEEBEfyLNPMzNrr6aJICI+AyDpBmCfiHgqff5pMqOIzMxsZMszamgy8Fzm+XPALoVEY2ZmpcvTWXwR8GtJPyIZ/nkk8L1CozIzs9LkuaHsC5KuA96YLjo+Iu4sNiwzMytL3hvKtgbWRsTXgYF0SKiZmfWAPGWoPwWcBpyRLtocuLjIoMzMrDx5rgiOBA4H/gIQEY/y96GlZmY2wuVJBM9FRJDWCZK0TbEhmZlZmfIkgh9K+hawvaQPAT8Dzi82LDMzK0ueUUNfkfQ2YC3JXcb/ERE3Fh6ZmZmVYshEIOnMiDgNuLHBMjMzG+HyNA29rcGyQ9odiJmZdUar6qMnAv8KTJXUn1m1HfCrogMzM7NytGoaugS4DvhP4PTM8qci4slCozIzs9K0qj66BlgDHAsg6aXAlsC2kraNiIfLCdHMzIqU587iwyQ9ADwILAAeIrlSMDOzHpCns/jzwEzg/nRqybfgPgIzs56RJxE8HxFPAKMkjYqIm4G9ig3LzMzKkmc+gj9L2pZkesrvS1oFrCs2LDMzK0ueK4IjgL8C/wb8FFhBMm+xmZn1gDwlJv4CIGkM8OPCIzIzs1LlKTHxYeCzJFcFLwIiqUQ6tdjQzMysDHn6CD4G7BkRq4sOxszMypenj2AF8EzRgZiZWWfkuSI4A1go6Xbg2drCiDi5sKjMzKw0eRLBt4CbgLtJ+gjMzKyH5EkE6yLi1MIjMTOzjsjTR3CzpDmSdpa0Y+2r8MjMzKwUea4Ijku/n5FZ5uGjZmY9Is8NZVPKCMTMzDqj1QxlB0XETZLe1Wh9RFxZXFhmZlaWVlcEB5CMFmpUVygAJwIzsx7QaoayT6UPPxsRD2bXSXJzkZlZj8gzauiKBsv+T7sDMTOzzmjVR7AHsCcwtq6fYAzJ3MVDkjQb+DowGjg/Ir5Ut/49wGnp06eBEyNiSf7wzcxsU7XqI9gdeCewPYP7CZ4CPjTUjiWNBs4F3gYMAHdIujoi7s1s9iBwQET8SdIhwHxg/436CczMbJO06iO4CrhK0usi4tZh7Hs/YHlErASQdCnJJDfrE0FELMxsfxswcRjHMTOzTZCnj+BISWMkbS7p55JWS3pvjtdNAB7JPB9IlzXzQeC6RivSO5sXSVr0+OOP5zi0mZnllScRHBwRa0maiQaAVwL/nuN1arAsGm4ovZkkEZzWaH1EzI+IGRExY/z48TkObWZmeeUpMbF5+v1Q4AcR8aTU6D1+AwPApMzzicCj9RtJmg6cDxwSEU/k2bGZmbVPniuCH0v6LTAD+Lmk8cDfcrzuDmA3SVMkbQEcA1yd3UDSZJIb0/45Iu7fuNDNzKwd8tQaOl3SmcDaiHhB0jMknb5DvW6dpLnA9STDRy+IiKWSTkjXzwP+A9gJ+GZ6lbEuImYM/8cxM7ONpYiGzfZI+h8R8eX08T9GxOWZdV+MiI+XFOMgM2bMiEWLFnXi0GZmI5akxc0+aLdqGjom8/iMunWzNzkqMzPrCq0SgZo8bvTczMxGqFaJIJo8bvTczMxGqFadxa+RtJbk0/9W6WPS57lqDZmZWfdrVWJidJmBmJlZZ+S5j8DMzHqYE4GZWcU5EZiZVZwTgZlZxTkRmJlVnBOBmVnFORGYmVWcE4GZWcU5EZiZVZwTgZlZxTkRmJlVnBOBmVnFORGYmVWcE4GZWcU5EZiZVZwTgZlZxTkRmJlVnBOBmVnFORGYmVWcE4GZWcU5EZiZVZwTgZlZxTkRmJlVnBOBmVnFORGYmVWcE4GZWcU5EZiZVZwTgZlZxTkRmJlVnBOBmVnFFZoIJM2WtEzSckmnN1gvSWen6/sl7VNkPGZmtqHCEoGk0cC5wCHANOBYSdPqNjsE2C39mgOcV1Q8ZmbWWJFXBPsByyNiZUQ8B1wKHFG3zRHA9yJxG7C9pJ0LjMnMzOpsVuC+JwCPZJ4PAPvn2GYC8Fh2I0lzSK4YAJ6WtCxnDOOA1XkDLlG3xgWObTi6NS5wbMPRrXHBpsX2D81WFJkI1GBZDGMbImI+MH+jA5AWRcSMjX1d0bo1LnBsw9GtcYFjG45ujQuKi63IpqEBYFLm+UTg0WFsY2ZmBSoyEdwB7CZpiqQtgGOAq+u2uRp4Xzp6aCawJiIeq9+RmZkVp7CmoYhYJ2kucD0wGrggIpZKOiFdPw+4FjgUWA48Axzf5jA2ujmpJN0aFzi24ejWuMCxDUe3xgUFxaaIDZrkzcysQnxnsZlZxTkRmJlVXE8mgqFKW5QcyyRJN0u6T9JSSR9Nl39a0u8l3ZV+HdqB2B6SdHd6/EXpsh0l3SjpgfT7Dh2Ia/fMeblL0lpJp3TqnEm6QNIqSfdkljU9T5LOSP/2lkl6e8lx/S9Jv01LtvxI0vbp8l0k/TVz7uYVFVeL2Jr+/so6Zy1iuywT10OS7kqXl3beWrxXFP+3FhE99UXSMb0CmApsASwBpnUwnp2BfdLH2wH3k5Tc+DTwsQ6fq4eAcXXLvgycnj4+HTizC36ffyC5GaYj5wx4E7APcM9Q5yn93S4BXgJMSf8WR5cY18HAZunjMzNx7ZLdrkPnrOHvr8xz1iy2uvVfBf6j7PPW4r2i8L+1XrwiyFPaojQR8VhE/CZ9/BRwH8nd093qCODC9PGFwH/rXCgAvAVYERG/61QAEfEL4Mm6xc3O0xHApRHxbEQ8SDIibr+y4oqIGyJiXfr0NpJ7c0rX5Jw1U9o5Gyo2SQLeDfygqOM30+K9ovC/tV5MBM3KVnScpF2AvYHb00Vz00v4CzrRBENyF/cNkhanZTwAXhbpvRzp95d2IK6sYxj8T9npc1bT7Dx109/fvwDXZZ5PkXSnpAWS3tihmBr9/rrpnL0R+GNEPJBZVvp5q3uvKPxvrRcTQa6yFWWTtC1wBXBKRKwlqbTaB+xFUlvpqx0I6/URsQ9JFdiTJL2pAzE0peRGxMOBy9NF3XDOhtIVf3+SPgGsA76fLnoMmBwRewOnApdIGlNyWM1+f11xzlLHMviDR+nnrcF7RdNNGywb1nnrxUTQdWUrJG1O8ov9fkRcCRARf4yIFyLiReDbFHgp3ExEPJp+XwX8KI3hj0orwKbfV5UdV8YhwG8i4o/QHecso9l56vjfn6T3A+8E3hNpY3LafPBE+ngxSXvyK8uMq8Xvr+PnDEDSZsC7gMtqy8o+b43eKyjhb60XE0Ge0halSdsc/zdwX0SclVmeLbd9JHBP/WsLjmsbSdvVHpN0Mt5Dcq7en272fuCqMuOqM+jTWafPWZ1m5+lq4BhJL5E0hWSujV+XFZSk2cBpwOER8Uxm+Xglc4QgaWoa18qy4kqP2+z319FzlvFW4LcRMVBbUOZ5a/ZeQRl/a2X0hpf9RVK24n6S7P2JDsfyBpLLtX7grvTrUOAi4O50+dXAziXHNZVkxMESYGntPAE7AT8HHki/79ih87Y18AQwNrOsI+eMJBk9BjxP8insg63OE/CJ9G9vGXBIyXEtJ2k3rv2tzUu3PSr9PS8BfgMc1oFz1vT3V9Y5axZbuvy7wAl125Z23lq8VxT+t+YSE2ZmFdeLTUNmZrYRnAjMzCrOicDMrOKcCMzMKs6JwMys4pwIbESQ9Im0ImN/WgVy/07HtCkkfVfS0QXs9+OZx7tkK2yaNeNEYF1P0utI7pTdJyKmk9z480jrV1XWx4fexGwwJwIbCXYGVkfEswARsTrS8hiS9k2LgS2WdH3mVvx9JS2RdKuSGv33pMs/IOmc2o4lXSPpwPTxwen2v5F0eVrzpTZvw2fS5XdL2iNdvq2k76TL+iUd1Wo/zbT4GW6RdKakX0u6v1bwTNLWkn6YHvMySbdLmiHpS8BW6RVTrcbQaEnfTq+mbpC0VVt+I9ZTnAhsJLgBmJS+GX5T0gGwvi7LN4CjI2Jf4ALgC+lrvgOcHBGvy3MASeOA/wm8NZJCfItIiozVrE6Xnwd8LF32SWBNRLw6vVK5Kcd+6o/b6meAZG6B/YBTgE+ly/4V+FN6zM8B+wJExOnAXyNir4h4T7rtbsC5EbEn8GeSO2XNBtms0wGYDSUinpa0L0mJ4DcDlymZeW4R8CrgxqRMC6OBxySNBbaPiAXpLi4iKWDXykySiT5+le5rC+DWzPpaAbDFJIXJIGmiOiYT558kvXOI/dTbvdHP0OS4u6SP3wB8PT3mPZL6W+z/wYi4q8E+zNZzIrARISJeAG4BbpF0N0nxrcXA0vpP/UqmZ2xWO2Udg6+Et6y9DLgxIo5t8rpn0+8v8Pf/GzU4zlD7qSca/Aw5jpvXs5nHLwBuGrINuGnIup6SOYx3yyzaC/gdSaGt8WlnMpI2l7RnRPwZWCPpDen278m89iFgL0mjJE3i76WQbwNeL2nXdF9bSxqq3PANwNxMnDsMYz8Nf4Yhjvv/SGbRQtI04NWZdc+nzU1muTkR2EiwLXChpHvTZpBpwKcjmYr0aOBMSUtIqjXOSl9zPHCupFuBv2b29SvgQZIqmF8hqShJRDwOfAD4QXqM24A9hojr88AOku5Jj//mjd3PED9DM98kSR79JCWn+4E16br5QH+ms9hsSK4+aj1PybR/10TEqzodSzsoqY+/eUT8TVIfSWniV6ZJxWyjuY/AbOTZGrg5bQIScKKTgG0KXxGYmVWc+wjMzCrOicDMrOKcCMzMKs6JwMys4pwIzMwq7v8D9XNgRxA7GuMAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(seq_lens, est_pr_survival, 'x')\n",
    "plt.ylim(0, 1)\n",
    "plt.xlabel('Sequence length')\n",
    "plt.ylabel('Estimated survival probability')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By learning the rate at which the survival probability decays to 0.5 (all measurements are just random), we can learn how strongly noise has affected our quantum device. In the next part of this sample, we look at ways we can learn this rate."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Inferring noise strengths from RB results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Curve fitting"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Following the derivation of [arXiv:1109.6887](https://arxiv.org/abs/1109.6887), the average survival probability for a sequence of length $m$ is given by\n",
    "$$\n",
    "    \\Pr(\\text{survive} | m; A, B, p) = Ap^m + B,\n",
    "$$\n",
    "where $A$, $B$, and $p$ parameterize the strength of the noise in terms of the initial state $\\rho$, the final measurement $E$, and the average gate fidelity $F(\\Lambda)$ over the noise in each different operation in the Clifford group,\n",
    "$$\n",
    "\\begin{aligned}\n",
    "    p & = \\frac{dF(\\Lambda) - 1}{d - 1} \\\\\n",
    "    A & = \\Tr(E\\Lambda(\\rho - 𝟙 / d)) \\\\\n",
    "    B & = \\Tr(E\\Lambda(𝟙 / d)),\n",
    "\\end{aligned}\n",
    "$$\n",
    "and where $d$ is the dimension of the Hilbert space that our state is defined on ($d = 2$ for a single qubit, as in this case).\n",
    "\n",
    "Effectively, $Ap^m + B$ acts as a model we can use to learn $p$, $A$, and $B$ from the results of measuring random RB sequences.\n",
    "\n",
    "To see how this works, let's use _curve fitting_ to find $p$, $A$, and $B$\n",
    "\n",
    "> **⚠ WARNING**: Curve fitting RB data can introduce some subtle errors into your estimate — we'll see more rigorous ways of learning $p$, $A$, and $B$ later in the sample.\n",
    "\n",
    "We can start by defining a Q# operation that runs an entire nonadaptive RB experiment."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%qsharp\n",
    "open Microsoft.Quantum.Convert;\n",
    "\n",
    "operation RunRBExperiment(sequenceLengths : Int[], nSequencesPerLength : Int) : Double[] {\n",
    "    mutable results = [];\n",
    "    for length in sequenceLengths {\n",
    "        let survivalAtLength = MeasureSingleQubitRBSequenceLength(length, nSequencesPerLength);\n",
    "        set results += [IntAsDouble(survivalAtLength) / IntAsDouble(nSequencesPerLength)];\n",
    "    }\n",
    "    return results;\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "est_pr_survival = RunRBExperiment.simulate_noise(sequenceLengths=seq_lens, nSequencesPerLength=200)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Estimated survival probability')"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEKCAYAAAAfGVI8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAel0lEQVR4nO3de7xVdbnv8c9XxLyCF2hv5ZKw8vLCMhESYleYZYGlHtO9j5pdrBNh0WV72ketU3Y/2q5OmSaRBzXNNHcXqXSrpVAnvEEI3hWQdIWKaAFlaeiz/xhj0VjTeRks5phzrjm+79drvtYclznmw1iT+azx+43f81NEYGZm5bVduwMwM7P2ciIwMys5JwIzs5JzIjAzKzknAjOzknMiMDMrucISgaT5ktZJurvGdkk6T9JKSSskHVpULGZmVluRVwSXADPqbJ8J7Jc+ZgEXFhiLmZnVUFgiiIhfAU/X2eVY4LuRuBXYXdLeRcVjZmbVbd/G9x4FPJpZ7k3XPVa5o6RZJFcN7LLLLpMOPPDAlgRoZtYtli5duj4iRlbb1s5EoCrrqta7iIh5wDyAyZMnx5IlS4qMy8ys60j6Xa1t7bxrqBcYk1keDaxtUyxmZqXVzkSwAHhXevfQVGBDRLyoWcjMzIpVWNOQpO8DhwMjJPUCZwNDASJiLnAtcBSwEngGOLWoWMzMrLbCEkFEnNRgewAfKur9+8xdtIqDRw9nWs+ILesWr1rPit4NzJ7eU/Tbm5l1vK4fWXzw6OHMuWIZi1etB5IkMOeKZRw8enibIzMz6wztvGuoJVb0buC0w8cz54plnDJlLJff9ginHT6eFb0b+l0lmJmVVSmuCC5cuJrp+4/gvJtWMn3/EVy4cLWvCMzMUl2fCKb1jOC0w8fzk2VrOWzfPfjJsrWcdvh4Xw2YmaW6vmlo8ar1XLhwNRP2Gcbta/7AcRNHceHC1Ry0T3JF4E5jMyu7rr8i6OsjeOTpZ9hx6Hb84r4nOO3w8fx0+Vp3GpuZUYIrgr67hr79zkkAfOCypXz1hgcZOmQ7vv3OSW4iMrPSK8UVwfknT2Razwim9Yzg4FHD+evfXuCVo/4+tmDxqvXMXbSqzZGambVH118RZNv/F69az4rfb2DHodtx1+83bBlbMOeKZZx/8sR2hWhm1lZdnwj69A0kyzYRvfeSO9xEZGal1/VNQ30qm4hOnbbvi5qIzMzKqDSJYPb0nn59Apff9ggfOeLl3P/4pi1NRGZmZVSaRNCnr4no/JMncvqbD+D8kyf2q0VkZlY2pUsE2SYiSEYen3/yRFb0bmhzZGZm7aGkGvTg4akqzcy2nqSlETG52rbSXRFAMkdBZVOQxxKYWVmVMhFk5yiYu2gV3/n1qn7lJpwUzKxMSpkI+voF5lyxjAce38iXfn7/loqknrjGzMqmNAPKKk3rGcEpU8Zy3k0rOW7iPly4cDWb/rKZy297pF9nsplZtyvlFQH0H0uw6MH1TN9/JOfdtJJTpox1EjCzUinlFUF2LMG0nhHsttP2fOnn93PQPsO4ePEapvbs1W/wmecsMLNuVsorguxYgr6Jaz7x1gO39At84LKlLF613v0FZlYKpbwiyP51XznA7OhX7cMHLlvK+Tet5P7HN7m/wMy6XimvCLKyNYiALQXpFq96yv0FZlYKpU8ElbKdyJff9ohrEJlZ13MiyHBBOjMrIyeCDBekM7MyctE5M7MScNG5JnPROjPrJk4EA5AtWgd4vIGZDWoNxxFI+iEwH7guIl4oPqTOly1ad8qUsa5PZGaDWp4rgguBk4GHJJ0j6cCCYxoUskXrPN7AzAazhokgIn4REe8ADgXWADdKWizpVElDiw6wU3m8gZl1i1x9BJL2At4D/A9gGfANksRwY2GRdTCPNzCzbtIwEUj6EfBrYGfg6Ig4JiKuiogPA7sWHWAn8ngDM+smDccRSDoqIq6tWPeSiHi24cGlGSRXD0OAiyLinIrtw4HLgbEkHddfiYiL6x3T4wjMzLbeto4j+EKVdbfkeNMhwAXATGACcJKkCRW7fQi4NyJeBRwOfFXSDjliMjOzJql5+6ikfwRGATtJmggo3TSMpJmokcOAlRGxOj3elcCxwL2ZfQLYTZJImpmeBjZv7T/CzMwGrt44greQdBCPBr6WWb8J+ESOY48CHs0s9wJTKvY5H1gArAV2A/57tbEKkmYBswDGjh2b463NzCyvmokgIi4FLpV0fET8cADHVpV1lR0SbwHuBI4AekhuTf11RGysiGUeMA+SPoIBxGJmZjXU7COQdEr6dF9Jp1c+chy7FxiTWR5N8pd/1qnAjyKxEngY6KgBa64rZGbdrl5n8S7pz11Jmm0qH43cAewnaVzaAXwiSTNQ1iPAGwEk/QNwALA6d/Qt4LpCZtbtCi1DLeko4Oskt4/Oj4gvSpoNEBFzJe0DXALsTdKUdE5EXF7vmO24fbTvy991hcxssKp3+2i9u4bOq3fQiPhIozdOxx9cW7Fubub5WuDNjY7Tbtm6Qh854uVOAmbWVerdNbS0ZVF0uMq6QlN79nIyMLOu0eiuodLL1hWa1jOCqT179VuGpEP54NHD+yWHxavWs6J3A7On97QrdDOzXOrdNfT19OdPJS2ofLQswjbLU1eoUYey7zwys05Ws7NY0qSIWCpperXtEbGo0Mhq6NRaQ/U6lCuvKiqXzcyKNqDO4ohYmv5clN7+eSDJgLAHIuK5QiIdxOp1KHtGMzPrZHnKUL8VWAWcR1ISYqWkmUUHNtg0mqjGM5qZWafKU330q8AbIuLwiJgOvAH4v8WGNbjkmajGM5qZWafKkwjWpeUf+qwG1hUUz6DUqEPZM5qZWSer11n89vTpkcDLgB+Q9BH8M0k/wf9sSYQVOrWzuB7fXmpm7Vavs7heIqg3U1hExHubEdzWGoyJwMys3QZ619CpxYVkZmadol6JCQAk7Qi8DzgI2LFvfbuuCMzMrLnydBZfBvwjySQyi0jmFdhUZFBmZtY6eRLByyPiU8Cf0/pDbwVeWWxY5eMyFGbWLnkSwd/Sn3+U9ApgOLBvYRGVlCfAMbN2adhHAMyTtAfwKZIZxnZNn1sTuQyFmbVLw0QQERelTxcB44sNp9w8AY6ZtUOeWkN7SfqmpN9KWirp65L2akVwZeMyFGbWDnn6CK4kKSlxPHACsB64qsigyshlKMysXfIkgj0j4vMR8XD6+AKwe8FxlU6eCXDMzIqQp7P4ZkknktQaguSq4OfFhVRO1WoOTesZ4X4CMytczUQgaRNJkTkBpwOXp5u2A/4EnF14dGZmVrh6tYZ2a2UgZmbWHnmahpB0DPD6dHFhRPysuJDMzKyV8tw+eg7wUeDe9PHRdJ2ZmXWBPHcNHQUcGRHzI2I+MCNdZ23iukRm1kx5EgH0v13UxW/azHWJzKyZ8vQRfAlYJulmkjuIXg+cVWhUVpfrEplZM9VNBJK2A14ApgKvJkkEZ0TE4y2IzepwXSIza5a6TUMR8QIwJyIei4gFEXGNk0BncF0iM2uWPH0EN0r6uKQxkvbsexQemdXkukRm1kx5+gj65ib+UGZd4JLUbVOvLpGbiMxsayki2h3DVpk8eXIsWbKk3WGYmQ0qkpZGxORq2xpeEUjaEfgg8FqSK4FfA3Mj4q9NjdLMzNoiTx/Bd4GDgG8C5wMTgMvyHFzSDEkPSFop6cwa+xwu6U5J90halDdwMzNrjjx9BAdExKsyyzdLWt7oRZKGABcARwK9wB2SFkTEvZl9dge+BcyIiEckvXSrojczs22W54pgmaSpfQuSpgC/yfG6w4CVEbE6Ip4jmens2Ip9TgZ+FBGPAETEunxhm5lZs+RJBFOAxZLWSFoD3AJMl3SXpBV1XjcKeDSz3Juuy9of2EPSwnQ+5HdVO5CkWZKWSFry5JNP5gjZzMzyytM0NGOAx1aVdZW3KG0PTALeCOwE3CLp1oh4sN+LIuYB8yC5a2iA8ZiZWRUNE0FE/G6Ax+4FxmSWRwNrq+yzPiL+DPxZ0q+AVwEPYmZmLZG3+uhA3AHsJ2mcpB2AE4EFFftcA7xO0vaSdiZphrqvwJjMzKxCYYkgIjYDc4DrSb7cfxAR90iaLWl2us99wH8CK4DbgYsi4u6iYrKB8xwIZt2ryCsCIuLaiNg/Inoi4ovpurkRMTezz79HxISIeEVEfL3IeMqqGV/ingPBrHvVTASSNknaWOWxSdLGVgZp26YZX+LZORC+dsMDW4reubaR2eBXs7M4InZrZSBWnGZNZOM5EMy6U+6mIUkvlTS271FkUNZ82S/xU6aMHdCXuOdAMOtODROBpGMkPQQ8DCwC1gDXFRyXNdm2fol7DgSz7pXniuDzJFNVPhgR40gGf+UpMWEdohlf4vXmQDCzwa3hfASSlkTE5LTQ3MSIeEHS7RFxWGtC7M/zEWy9uYtWcfDo4f2agxavWs+K3g3Mnt7TxsjMrFW2aT4C4I+SdgV+BXxP0jpgczMDtGJV+7Kf1jPCnb1mBuRrGjoWeAb4V5LBX6uAo4sMyrqfB6i1l8+/ZeVJBLOAfSJic0RcGhHnRcRTRQdm3c0D1NrL59+y8jQNDQOul/Q0yZwC/xERTxQblnW7Zo1tsIHx+beshlcEEfHZiDgI+BCwD7BI0i8Kj8y6XjPGNuThZpDqWnX+rfNtTa2hdcDjwFOAp5S0bdaqAWpuBqnOAwStT54BZadJWgj8EhgBvD8iDi46MOturRyg1s46SZ16NeIBgpaV54rgZcDHIuKgiDg7O/m82UC1eoBau5pBOvVqxAMELavmgDJJwyJio6Q9q22PiKcLjawGDyizgej7Am5Hx2g739usz0AHlF0BvA1YSjLXcHYO4gDGNy1CswJlm0Gm9Yxgas9eLW0ectVW63Q1m4Yi4m3pz3ERMT792fdwErBBo93NIO6UtU6Xp7P4GkknpXMKmw06s6f3vOiv8Gk9I1pSZ8mdsratWnHDQZ7O4q8BrwPuk3S1pBMk7di0CMy6WLuvRmzwa8UNBw2rj27ZURoCHAG8H5gREcOaFsVWcGexmZVNM244qNdZnGtAmaSdgOOB2cCrgUu3KgIzMxuwom9/ztNHcBVwH8nVwAVAT0R8uKlRmJlZTUXfcFC36Jyk7YC7gJMj4vmmvrOZmTXUituf614RRMQLwFudBMysk3Rq6Y4itOKGgzx9BDdIOl6SGu9qZla8Ti3dUYRW3P6cZz6C04FdgM2S/koywjjaddeQmZnnU2iuhokgInZrRSBmZlvDpTuap2EikPT6ausj4lfND8fMLJ/KO2mm9uzlZDBAeZqG/i3zfEfgMJJCdEcUEpGZWQPtLiTYbfI0DR2dXZY0BvhyYRGZmTVQ704aJ4Ktl+eKoFIv8IpmB2Jmlle1O2am9YxwEhigPH0E3ySZfwCS200PAZYXGJOZmbVQniuCbIW3zcD3I+I3BcVjZmYtlqePYEuBOUl7AGMKjcjMzFoqT9G5hZKGpXMXLwculvS14kMzM7NWyFNiYnhEbATeDlwcEZOAN+U5uKQZkh6QtFLSmXX2e7Wk5yWdkC9sMzNrljyJYHtJewP/Avws74HTiWwuAGYCE4CTJE2osd+5wPV5j21mZs2TJxF8juRLemVE3CFpPPBQjtcdlr5mdUQ8B1wJHFtlvw8DPwTW5YzZzMyaKE9n8dXA1Znl1SSzlTUyCng0s9wLTMnuIGkUcBzJKOVX1zqQpFnALICxY8fmeGszM8sr11SVA1StbHXlBMlfB85oNN9BRMyLiMkRMXnkyJHNis/MzBjYyOK8eul/q+loYG3FPpOBK9OpDkYAR0naHBE/KTAuMzPLKDIR3AHsJ2kc8HvgRODk7A4RMa7vuaRLgJ85CZiZtVbNRCDp9HovjIi6YwkiYrOkOSQdzUOA+RFxj6TZ6fa5A4jXzMyarN4VQd+ENAeQdOQuSJePBnLNRRAR1wLXVqyrmgAi4j15jmlmZs1VMxFExGcBJN0AHBoRm9Llz5C5i8jMzAa3PHcNjQWeyyw/B+xbSDRmZtZyeTqLLwNul/Rjkts/jwO+W2hUZmbWMnkGlH1R0nXA69JVp0bEsmLDMjOzVsk7oGxnYGNEfAPoTW8JNTOzLpCnDPXZwBnAWemqocDlRQZlZmatk+eK4DjgGODPABGxlr/fWmpmZoNcnkTwXEQEaZ0gSbsUG5KZmbVSnkTwA0nfBnaX9H7gF8BFxYZlZmatkueuoa9IOhLYSDLK+NMRcWPhkZmZWUs0TASSzo2IM4Abq6wzM7NBLk/T0JFV1s1sdiBmZtYe9aqPngZ8EBgvaUVm027Ab4oOzMzMWqNe09AVwHXA/wHOzKzfFBFPFxqVmZm1TL3qoxuADcBJAJJeCuwI7Cpp14h4pDUhmplZkfKMLD5a0kPAw8AiYA3JlYKZmXWBPJ3FXwCmAg+mU0u+EfcRmJl1jTyJ4G8R8RSwnaTtIuJm4JBiwzIzs1bJMx/BHyXtSjI95fckrQM2FxuWmZm1Sp4rgmOBvwD/CvwnsIpk3mIzM+sCeUpM/BlA0jDgp4VHZGZmLZWnxMQHgM+RXBW8AIikEun4YkMzM7NWyNNH8HHgoIhYX3QwZmbWenn6CFYBzxQdiJmZtUeeK4KzgMWSbgOe7VsZER8pLCozM2uZPIng28BNwF0kfQRmZtZF8iSCzRFxeuGRmJlZW+TpI7hZ0ixJe0vas+9ReGRmZtYSea4ITk5/npVZ59tHzcy6RJ4BZeNaEYiZmbVHvRnKjoiImyS9vdr2iPhRcWGZmVmr1LsimE5yt1C1ukIBOBGYmXWBejOUnZ0+/VxEPJzdJsnNRWZmXSLPXUM/rLLuP5odiJmZtUe9PoIDgYOA4RX9BMNI5i5uSNIM4BvAEOCiiDinYvs7gDPSxT8Bp0XE8vzhm5nZtqrXR3AA8DZgd/r3E2wC3t/owJKGABcARwK9wB2SFkTEvZndHgamR8QfJM0E5gFTtupfYGZm26ReH8E1wDWSXhMRtwzg2IcBKyNiNYCkK0kmudmSCCJicWb/W4HRA3gfMzPbBnn6CI6TNEzSUEm/lLRe0ik5XjcKeDSz3Juuq+V9wHXVNqQjm5dIWvLkk0/meGszM8srTyJ4c0RsJGkm6gX2B/4tx+tUZV1U3VF6A0kiOKPa9oiYFxGTI2LyyJEjc7y1mZnllafExND051HA9yPiaanad/yL9AJjMsujgbWVO0k6GLgImBkRT+U5sJmZNU+eK4KfSrofmAz8UtJI4K85XncHsJ+kcZJ2AE4EFmR3kDSWZGDaOyPiwa0L3czMmiFPraEzJZ0LbIyI5yU9Q9Lp2+h1myXNAa4nuX10fkTcI2l2un0u8GlgL+Bb6VXG5oiYPPB/jpmZbS1FVG22R9L/iogvp8//OSKuzmz7UkR8okUx9jN58uRYsmRJO97azGzQkrS01h/a9ZqGTsw8P6ti24xtjsrMzDpCvUSgGs+rLZuZ2SBVLxFEjefVls3MbJCq11n8KkkbSf763yl9Trqcq9aQmZl1vnolJoa0MhAzM2uPPOMIzMysizkRmJmVnBOBmVnJORGYmZWcE4GZWck5EZiZlZwTgZlZyTkRmJmVnBOBmVnJORGYmZWcE4GZWck5EZiZlZwTgZlZyTkRmJmVnBOBmVnJORGYmZWcE4GZWck5EZiZlZwTgZlZyTkRmJmVnBOBmVnJORGYmZWcE4GZWck5EZiZlZwTgZlZyTkRmJmVnBOBmVnJORGYmZWcE4GZWck5EZiZlVyhiUDSDEkPSFop6cwq2yXpvHT7CkmHFhmPmZm9WGGJQNIQ4AJgJjABOEnShIrdZgL7pY9ZwIVFxWNmZtUVeUVwGLAyIlZHxHPAlcCxFfscC3w3ErcCu0vau8CYzMyswvYFHnsU8GhmuReYkmOfUcBj2Z0kzSK5YgD4k6QHcsYwAlifN+AW6tS4wLENRKfGBY5tIDo1Lti22F5Wa0ORiUBV1sUA9iEi5gHztjoAaUlETN7a1xWtU+MCxzYQnRoXOLaB6NS4oLjYimwa6gXGZJZHA2sHsI+ZmRWoyERwB7CfpHGSdgBOBBZU7LMAeFd699BUYENEPFZ5IDMzK05hTUMRsVnSHOB6YAgwPyLukTQ73T4XuBY4ClgJPAOc2uQwtro5qUU6NS5wbAPRqXGBYxuITo0LCopNES9qkjczsxLxyGIzs5JzIjAzK7muTASNSlu0OJYxkm6WdJ+keyR9NF3/GUm/l3Rn+jiqDbGtkXRX+v5L0nV7SrpR0kPpzz3aENcBmfNyp6SNkj7WrnMmab6kdZLuzqyreZ4knZV+9h6Q9JYWx/Xvku5PS7b8WNLu6fp9Jf0lc+7mFhVXndhq/v5adc7qxHZVJq41ku5M17fsvNX5rij+sxYRXfUg6ZheBYwHdgCWAxPaGM/ewKHp892AB0lKbnwG+Hibz9UaYETFui8DZ6bPzwTO7YDf5+Mkg2Hacs6A1wOHAnc3Ok/p73Y58BJgXPpZHNLCuN4MbJ8+PzcT177Z/dp0zqr+/lp5zmrFVrH9q8CnW33e6nxXFP5Z68YrgjylLVomIh6LiN+mzzcB95GMnu5UxwKXps8vBf5b+0IB4I3Aqoj4XbsCiIhfAU9XrK51no4FroyIZyPiYZI74g5rVVwRcUNEbE4XbyUZm9NyNc5ZLS07Z41ikyTgX4DvF/X+tdT5rij8s9aNiaBW2Yq2k7QvMBG4LV01J72En9+OJhiSUdw3SFqalvEA+IdIx3KkP1/ahriyTqT/f8p2n7M+tc5TJ33+3gtcl1keJ2mZpEWSXtemmKr9/jrpnL0OeCIiHsqsa/l5q/iuKPyz1o2JIFfZilaTtCvwQ+BjEbGRpNJqD3AISW2lr7YhrH+KiENJqsB+SNLr2xBDTUoGIh4DXJ2u6oRz1khHfP4kfRLYDHwvXfUYMDYiJgKnA1dIGtbisGr9/jrinKVOov8fHi0/b1W+K2ruWmXdgM5bNyaCjitbIWkoyS/2exHxI4CIeCIino+IF4DvUOClcC0RsTb9uQ74cRrDE0orwKY/17U6royZwG8j4gnojHOWUes8tf3zJ+ndwNuAd0TamJw2HzyVPl9K0p68fyvjqvP7a/s5A5C0PfB24Kq+da0+b9W+K2jBZ60bE0Ge0hYtk7Y5/j/gvoj4WmZ9ttz2ccDdla8tOK5dJO3W95ykk/FuknP17nS3dwPXtDKuCv3+Omv3OatQ6zwtAE6U9BJJ40jm2ri9VUFJmgGcARwTEc9k1o9UMkcIksanca1uVVzp+9b6/bX1nGW8Cbg/Inr7VrTyvNX6rqAVn7VW9Ia3+kFStuJBkuz9yTbH8lqSy7UVwJ3p4yjgMuCudP0CYO8WxzWe5I6D5cA9fecJ2Av4JfBQ+nPPNp23nYGngOGZdW05ZyTJ6DHgbyR/hb2v3nkCPpl+9h4AZrY4rpUk7cZ9n7W56b7Hp7/n5cBvgaPbcM5q/v5adc5qxZauvwSYXbFvy85bne+Kwj9rLjFhZlZy3dg0ZGZmW8GJwMys5JwIzMxKzonAzKzknAjMzErOicAGBUmfTCsyrkirQE5pd0zbQtIlkk4o4LifyDzfN1th06wWJwLreJJeQzJS9tCIOJhk4M+j9V9VWp9ovItZf04ENhjsDayPiGcBImJ9pOUxJE1Ki4EtlXR9Zij+JEnLJd2ipEb/3en690g6v+/Akn4m6fD0+ZvT/X8r6eq05kvfvA2fTdffJenAdP2uki5O162QdHy949RS59+wUNK5km6X9GBfwTNJO0v6QfqeV0m6TdJkSecAO6VXTH01hoZI+k56NXWDpJ2a8huxruJEYIPBDcCY9MvwW5Kmw5a6LN8EToiIScB84Ivpay4GPhIRr8nzBpJGAP8beFMkhfiWkBQZ67M+XX8h8PF03aeADRHxyvRK5aYcx6l833r/BkjmFjgM+Bhwdrrug8Af0vf8PDAJICLOBP4SEYdExDvSffcDLoiIg4A/koyUNetn+3YHYNZIRPxJ0iSSEsFvAK5SMvPcEuAVwI1JmRaGAI9JGg7sHhGL0kNcRlLArp6pJBN9/CY91g7ALZntfQXAlpIUJoOkierETJx/kPS2BsepdEC1f0ON9903ff5a4Bvpe94taUWd4z8cEXdWOYbZFk4ENihExPPAQmChpLtIim8tBe6p/KtfyfSMtWqnbKb/lfCOfS8DboyIk2q87tn05/P8/f+NqrxPo+NUElX+DTneN69nM8+fB9w0ZC/ipiHreErmMN4vs+oQ4HckhbZGpp3JSBoq6aCI+COwQdJr0/3fkXntGuAQSdtJGsPfSyHfCvyTpJenx9pZUqNywzcAczJx7jGA41T9NzR43/9PMosWkiYAr8xs+1va3GSWmxOBDQa7ApdKujdtBpkAfCaSqUhPAM6VtJykWuO09DWnAhdIugX4S+ZYvwEeJqmC+RWSipJExJPAe4Dvp+9xK3Bgg7i+AOwh6e70/d+wtcdp8G+o5VskyWMFScnpFcCGdNs8YEWms9isIVcfta6nZNq/n0XEK9odSzMoqY8/NCL+KqmHpDTx/mlSMdtq7iMwG3x2Bm5Om4AEnOYkYNvCVwRmZiXnPgIzs5JzIjAzKzknAjOzknMiMDMrOScCM7OS+y9boUHALKLf7gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(seq_lens, est_pr_survival, 'x')\n",
    "plt.ylim(0, 1)\n",
    "plt.xlabel('Sequence length')\n",
    "plt.ylabel('Estimated survival probability')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once we have these results, we can use the `curve_fit` function to roughly estimate $p$, $A$, and $B$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.optimize import curve_fit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "model = (lambda m, p, A, B: A * p**m + B)\n",
    "p_opt, p_cov = curve_fit(\n",
    "    model,\n",
    "    seq_lens,\n",
    "    est_pr_survival,\n",
    "    [0.95, 0.5, 0.5]\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.96854501, 0.52046551, 0.4814992 ])"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p_opt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x1a2209d9fa0>"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEKCAYAAAAfGVI8AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAuj0lEQVR4nO3deXxU5dn/8c/FJjsooGVfIoKgEGQVRZZqRcS94l63arFFbf3VovWxtdb2UR/1qYqKS6lUUFHQB+pSsSpoDVpAFlkECaBGEAElILIk5Pr9cU7iJCaTQ8jMJJnv+/WaV2buOXPONSfJXHPOfe7rNndHRETSV61UByAiIqmlRCAikuaUCERE0pwSgYhImlMiEBFJc0oEIiJpLmGJwMwmmdmXZrasjOfNzB4wszVmttTMjklULCIiUrZEHhE8CYyM8/wpQNfwdjXwSAJjERGRMiQsEbj728BXcRY5A/i7B94DmptZ60TFIyIipauTwm23BT6LeZwTtm0suaCZXU1w1ECjRo36du/ePSkBiojUFAsXLtzi7q1Key6VicBKaSu13oW7PwY8BtCvXz9fsGBBIuMSEalxzOyTsp5L5VVDOUD7mMftgA0pikVEJG2lMhHMAn4SXj00CMh19++dFhIRkcRK2KkhM3sGGAa0NLMc4PdAXQB3nwi8AowC1gDfApcnKhYRESlbwhKBu19QzvMO/CJR2y80cW42vdo1Y3BGy6K2rOwtLM3JZezQjERvXkTKkZeXR05ODrt37051KDVC/fr1adeuHXXr1o38mlR2FidFr3bNGPf0IiZc2IfBGS3Jyt5S9FhEUi8nJ4cmTZrQqVMnzEq7hkSicne2bt1KTk4OnTt3jvy6Gl9iYmlOLtcM68K4pxdx3+xVjHt6EdcM68LSnNxUhyYiwO7du2nRooWSQCUwM1q0aLHfR1c1PhH0ateMR+asZegRLXngzTUMPaIlj8xZS692zVIdmoiElAQqT0X2ZY1PBIMzWnLNsC7MX7yE/20xk38s+oxrhnUp1mcgIpLOanwiyMrewiNz1nJyi82ctXMav++yikfmrCUrewtZ2VuYODc71SGKSIqZGZdccknR4/z8fFq1asXo0aNTGFX5GjduXCnrqfGJoLCP4PlvjmK1t2PQhslcM7QT/1iygXFPL9IpIpFqZOLcbLKytxRrq4wvdI0aNWLZsmXs2rULgNdff522bdse0DorKj8/P+nbrPGJoLCPYOIl/eH4G+jKZyx6/VleWrqx6EoiEakeCq8CLEwGhVcBVsYXulNOOYWXX34ZgGeeeYYLLvjuCvidO3dyxRVX0L9/f/r06cPMmTMBWL9+PUOGDOGYY47hmGOOISsrC4CNGzdywgknkJmZyVFHHcU777wDFP8GP336dC677DIALrvsMm644QaGDx/O+PHjyc7OZuTIkfTt25chQ4bw0UcfAbBu3TqOPfZY+vfvz6233nrA77mIu1erW9++fX1/PDJnjb+7ZnPwID/Pv7i9my++NdMveDSraJl312z2R+as2a/1ikjlWLFixX4t/+6azd7n9tl+72sfeZ/bZ3/3/30AGjVq5EuWLPFzzjnHd+3a5b179/a33nrLTz31VHd3v/nmm/2pp55yd/evv/7au3bt6t98843v3LnTd+3a5e7uq1ev9sLPp3vuucfvuOMOd3fPz8/37du3F22n0PPPP++XXnqpu7tfeumlfuqpp3p+fr67u48YMcJXr17t7u7vvfeeDx8+3N3dTzvtNJ88ebK7u0+YMKHY+mKVtk+BBV7G52qNH0cQO2gsa/02/rnvDG6vNZFmG+aSld0VQOMKRKqRwRktuXhgBx54cw3XjTi80o7qe/Xqxfr163nmmWcYNWpUsedmz57NrFmzuOeee4DgktdPP/2UNm3aMG7cOBYvXkzt2rVZvXo1AP379+eKK64gLy+PM888k8zMzHK3f+6551K7dm2++eYbsrKyOPfcc4ue27NnDwDvvvsuM2bMAOCSSy5h/PjxlfHWa34iKFR4CPnQhdeze+Y/GLtzOuc92Yu6tWvz6CV9dYpIpJrIyt7ClPc/5boRhzPl/U8ZlNGi0v5/Tz/9dH79618zZ84ctm7dWtTu7syYMYNu3boVW/62227jsMMOY8mSJRQUFFC/fn0ATjjhBN5++21efvllLrnkEm688UZ+8pOfFLu0s+S1/o0aNQKgoKCA5s2bs3jx4lJjTMSltjW+j6DQ0pxcJlzYh2OPaEP9Yf+P3nzMgH2LObptMyUBkWoitjLADT/qxoQL+xTrMzhQV1xxBb/73e84+uiji7WffPLJPPjggwRnWGDRokUA5Obm0rp1a2rVqsVTTz3Fvn37APjkk0849NBDueqqq7jyyiv54IMPADjssMNYuXIlBQUFvPjii6XG0LRpUzp37szzzz8PBEloyZIlABx33HE8++yzAEydOrVS3jOkUSIYOzSj6AN/XrORbKQldx3yEh9t3F5pf0QikliFX+gK/5cHZ7RkwoV9Kq1SQLt27bj++uu/137rrbeSl5dHr169OOqoo4o6an/+858zefJkBg0axOrVq4u+1c+ZM4fMzEz69OnDjBkzitZ55513Mnr0aEaMGEHr1mVPyDh16lT++te/0rt3b3r27FnUOX3//ffz0EMP0b9/f3JzK686ghVmuOriQCemKfxG8Vz/VRz+3i2sHP4EF73dXFcQiaTIypUrOfLII1MdRo1S2j41s4Xu3q+05dPmiKBQ4TeKw0/6GTTvyJErH2DCBb1Ve0hE0lbadBYXKlZ6evhv4cWfMXjX2wwe+uPUBSUikkJpd0QAMaMTjz4XDu0Jb97BvI83qtyEiKSltEwERaMT133NK4ddDV+v462n7y0anagaRCKSTtIyERReaTDu6UW8nteL+QXd+GXdFxncvkGlDlkXEakO0jIRwHejE19cvJGsztfScO8W/j3l9mKzmYmIpIO0TQSxoxMn5/yADxsfR+9PJnNVn8ZKAiJp5osvvuD8888nIyODHj16MGrUqKJyEcmwefNmBg4cSJ8+fXjnnXcYNWoU27ZtY9u2bTz88MMJ335aJoKSoxOvGdaFX249iwa2lxYL7ik2wEz9BSI1m7tz1llnMWzYMLKzs1mxYgV//vOf2bRpU+R1FI4orqg33niD7t27s2jRIoYMGcIrr7xC8+bNlQgSKXZ0YuHENeePGsH7Lc/mHN7g3qdeKJq4Rv0FIjXbW2+9Rd26dRk7dmxRW2ZmJkOGDGHOnDnFJqcZN24cTz75JACdOnXi9ttv5/jjj+fuu+9mwIABRcutX7+eXr16AbBw4UKGDh1K3759Ofnkk9m4cWOx7S9evJjf/OY3vPLKK2RmZrJr1y46derEli1buOmmm8jOziYzM5Mbb7wxYfsg7cYRQPGxBMWGrPe9m7y/zObGvX/ngTeO5qNN36i/QCSZXr0Jvviwctf5g6PhlDvLfHrZsmX07du3QquuX78+//73vwGYNm0aa9eupUuXLkybNo0xY8aQl5fHtddey8yZM2nVqhXTpk3jlltuYdKkSUXryMzM5Pbbb2fBggVMmDCh2PrvvPNOli1bVmYBusqSlkcEsWJrENHwEOr+8BYG8SEN17/OxQM7KAmISJnOO++8ovtjxozhueeeA4KkcN5557Fq1SqWLVvGSSedRGZmJnfccQc5OTmpCrdMaXlEEM+8Q07nB9zP/zR9jpHv9avUErciUo4439wTpWfPnkyfPr3U5+rUqUNBQUHR47JKR0OQFM4991zOPvtszIyuXbvy4Ycf0rNnT+bNm5eY4CtJ2h8RxMrK3sIvnl3Grh/ewcG7P+P5o9+v1BK3IlL1jBgxgj179vD4448Xtc2fP5+5c+fSsWNHVqxYwZ49e8jNzeWNN94ocz0ZGRnUrl2bP/7xj0VHCt26dWPz5s1FiSAvL4/ly5dHjq1Jkybs2LGjgu8sOiWCGIX9BT2GnA09z6LDskd4YvTBKkgnUoOZGS+++CKvv/46GRkZ9OzZk9tuu402bdrQvn17xowZQ69evbjooovo0yf+TIbnnXceU6ZMYcyYMQDUq1eP6dOnM378eHr37k1mZmbRvMZRtGjRguOOO46jjjoqoZ3FaVeGOrLtG2FCf2jfHy5+ARIwK5CIqAx1IqgMdWVp2hp+eCtkvwnLi88kVFS0LobGG4hIdaVEEE//n0LrTPjnzbD7u9NDRUXrwmSg8QYiUp2VmwjMbIaZnWpm6Zc0atWG0f8L32yCN/9U1BxbtO6+2atUn0jkAFW3U9RVWUX2ZZQP90eAC4GPzexOM+u+31upztoeAwOugvmPw+cfFDUXFq174M01Gm8gcgDq16/P1q1blQwqgbuzdetW6tevv1+vK3ccgbv/C/iXmTUDLgBeN7PPgMeBKe6eV5GAq5UR/wUrZsKsa+Gqt6BOvWJF66a8/6nGG4hUULt27cjJyWHz5s2pDqVGqF+/Pu3atduv10QaUGZmLYCLgUuARcBU4HjgUmDYfm2xOqrfDEb/BZ69AN6+m6yOY4udDhqU0UKnh0QqqG7dunTu3DnVYaS1KH0ELwDvAA2B09z9dHef5u7XAo0THWCV0X0UZF4E79zHFyveLfahX9hnoPEGIlIdlTuOwMxGufsrJdoOcvc95a7cbCRwP1AbeMLd7yzxfDNgCtCB4OjkHnf/W7x1Jm0cQWl258LDg6FeQ/jZ21C3QWriEBHZTwc6juCOUtrKLZxhZrWBh4BTgB7ABWbWo8RivwBWuHtvglNM95pZvQgxpUb9ZnDGg7BlNbxZ2m4REal+yuwjMLMfAG2BBmbWBygcWtuU4DRReQYAa9x9bbi+Z4EzgBUxyzjQxMyM4DTTV0D+/r6JpMoYAf2uhHkPQfdToePgVEckInJA4nUWnwxcBrQD7otp3wH8NsK62wKfxTzOAQaWWGYCMAvYADQBznP3ghLLYGZXA1cDdOjQIcKmE+yk2yH7Dfi/a2Dsu3BQ+nSViEjNU+apIXef7O7DgcvcfXjM7XR3fyHCuksrzlOyQ+JkYDHQBsgEJphZ01Jieczd+7l7v1atWkXYdIId1BjOfAS2fQov3wC6/llEqrEyE4GZXRze7WRmN5S8RVh3DtA+5nE7gm/+sS4HXvDAGmAdUKUGrJVZV2j9YTB0PCydBoumpCg6EZEDF6+zuHDGhcYEp21K3sozH+hqZp3DDuDzCU4DxfoU+CGAmR0GdAPWRo4+CeLWFTrhRuh8ArxyI3y5MsWRiohUTELLUJvZKOAvBJePTnL3P5nZWAB3n2hmbYAngdYEp5LudPe4X69Tcflo4Yf/xQM7MOX9T4sPHNuxCSYeDw0OhqvfgnqN4q9MRCQF4l0+WmYiMLMH4q3U3a+rhNj2W6rGEdw3exUPvLmG60Yczg0/6lb8ybVz4O9nQu8L4KxHkh6biEh54iWCeFcNLUxQPNVOuXWFugwLThO9fTd0HgKZF6YsVhGR/VVmInD3yckMpKoqPC0Ur67QxLnZ9Gr/UwZ3mgcv/QpadiNrT0eW5uQydmhGit+BiEh88a4a+kv48x9mNqvkLWkRpljhPMbx6gr1ateMcc8uZX6/e6DxoeydMoY7ps4umqhGM5qJSFUWr4+gr7svNLOhpT3v7nMTGlkZUlprKI7CI4frjs7nnMVXUKt5expd8wYc1OR7RxUlH4uIJFqFag25+8Lw51yC2kJfE5SAmJeqJFCVFU5Uc9t7BczucSeNctfA9CuhYJ9mNBORKi1KGepTgWzgAYKSEGvM7JREB1bdxHYo/2lVG7IH3AYfvwav3QJoRjMRqbqiTExzLzA8HPmLmWUALwOvJjKw6qS0DuVzn4aXul9Om/cfgRYZZLU4SzOaiUiVFCURfFmYBEJrgS8TFE+1VFaH8j8+68jPCjbir45nqm1jwkWXa0YzEaly4nUWnx3ePQnoCDxHUDTuXGCVu/+/pERYQlXtLC7Tnh1seWA4zfdupM5V/4JDjwSCowhdXioiyVLRkcXxZgpzd7+iMoLbX9UuEQDk5sDjI6DOQfDTN6FxFaigKiJppUIji9398sSFlGaatYMLnoW/jYIpZ8FPZkHDQ1IdlYgIEO2qofpm9gsze9jMJhXekhFcjdL2GDh/CmxeDZNPh51bUx2RiAgQbc7ip4AfEEwiM5dgXoEdiQyqxjr8RLjwWdj6MUw+DXZuKf81IiIJFiURHO7utwI7w/pDpwJHJzasGixjBFw4Db5aC0+Ohm82AypDISKpEyUR5IU/t5nZUUAzoFPCIkoHXYYFyeDr9TB5NHzzZfwJcEREEihKInjMzA4GbiWYYWwFcFdCo0oHXYbCxdODeY+fPJXBh+arDIWIpES5icDdn3D3r919rrt3cfdD3f3RZARX43U6Hi6eAbmfB8mgVZ7KUIhI0kW5aqiFmT1oZh+Y2UIz+4uZtUhGcGmh4+AgGez4gl2Pj+S19xYVlaEo2WcgIpIIUU4NPUtQUuIc4MfAFmBaIoNKOx2P5cPhk9i3YxP/1+jP3DCwcdFpIiUDEUm0KIngEHf/o7uvC293AM0THFfaeXfv4awd+RQN9n4Nk0YyuOGG702AIyKSCFESwVtmdr6Z1QpvYwiqj0olGjs0g16DToJLZ0FBPvz1Rwz+do5qEYlIwsWbqnKHmW0HfgY8DewNb88Cv0pOeGmoTR/42VxokwkzroTXfwcF+1IdlYjUYPFmKGvi7k3Dn7XcvU54q+XuTZMZZNppfGhQj6jflfDu/TD1x/DtV6mOSkRqqCinhjCz083snvA2OtFBCVCnHoy+D057ANa9A48Ph03LUx2ViNRAUS4fvRO4nmAg2Qrg+rBNkqHvpXD5K5C3G544CVbMTHVEIlLDRDkiGAWc5O6T3H0SMDJsk2RpPwCungOH9YDnfsLCv91A1sebii2iukQiUlGRTg1R/HJRFb9Jhaat4bKXoc8l9P3kr+ybej7vr1wLqC6RiByYKHMW/xlYZGZvAQacANyc0KikdHUOgtMfhDaZHPfKeD6bNorJmfdx/4d1VZdIRCosbiIws1pAATAI6E+QCMa7+xdJiE1KYwb9f0qtQ3vQcspFnL/oErp2uILBHU9IdWQiUk3FPTXk7gXAOHff6O6z3H2mkkDVkJV/BGfsu5t1LYcz+LPH+PbBwfDJvFSHJSLVUJQ+gtfN7Ndm1t7MDim8JTwyKVNhn8DtFw2n+7XTWTHiCXJzt8HfRsJLv4Jd21IdoohUI+bu8RcwW1dKs7t7l8SEFF+/fv18wYIFqdh0lTFxbja92jUr1ifw3kef0ODfd9H782eg0aEw6m448vTgVJKIpD0zW+ju/Up9rrxEUNUoEZRjwyKYdR18sRS6jYJR/wPN2qU6KhFJsXiJIMqAsvpmdoOZvWBmM8zsl2ZWv/LDlErRpg9c9Rac9EfIfgseGgjvP6p6RSJSpih9BH8HegIPAhOAHsBTUVZuZiPNbJWZrTGzm8pYZpiZLTaz5WY2N2rgEkftOnDcdfCL96D9QHj1N/DXk+CLZamOTESqoCh9BEvcvXd5baW8rjawGjgJyAHmAxe4+4qYZZoDWcBId//UzA519y/jrVenhvaTO3w4Hf55E+zeBsf+Ao77JTRUf79IOjmgU0MEg8kGxaxsIPBuhNcNANa4+1p3LyxffUaJZS4EXnD3TwHKSwJSAWbQ61wYNx96nRdUM/3fnvDqTZCbk+roRKQKiJIIBgJZZrbezNYD84ChZvahmS2N87q2wGcxj3PCtlhHAAeb2ZxwPuSflLYiM7vazBaY2YLNmzdHCFm+p+EhcObD8PP3oMcZMP9xuL83vHgNfPlRqqMTkRSKUmJiZAXXXdp1iyXPQ9UB+gI/BBoA88zsPXdfXexF7o8Bj0FwaqiC8QjAoUfCWRNh+G9h3kOwcDIseRq6nQrH/wra9091hCKSZOUmAnf/pILrzgHaxzxuB2woZZkt7r4T2GlmbwO9CfoWJJGad4BT7oITfgP/eTS4smjVy9Dx+CAhHP5DjUEQSRNRq49WxHygq5l1NrN6wPnArBLLzASGmFkdM2tIcBpqZQJjkpIatQiODn61HE7+M3y1FqaeAxOHBJ3M+/JTHaGIJFjCEoG75wPjgNcIPtyfc/flZjbWzMaGy6wE/gksBf4DPOHuusYxFQ5qHFxRdP0SOONh2LcnmDN5Ql+Y/1cef3M5Wdlbir1EcyCI1AwaWZwGSitJkZW9haU5uYwdmlH6iwoKYNUr8O/74POF7K3fgkf3nMyAMTcy8MguRfWOVP5apHqo0OWjZrbDzLaXctthZtsTF65Utl7tmjHu6UVF3+gjTWRTqxYcORp++gZc+hL12mZyrT9Nz2mDWfDoNTw09XkmXJCpJCBSA+iIIE0UfvhfPLADU97/tGLf5Dcu4aPpf+TwLW9QxwqgWXvoPhqOPA06DIJatRMTvIgcsHhHBFEuHy1cyaFAUY2hwkFgUj0MzmjJxQM78MCba7huxOEV+iaf9W1bxm37GT8dcCObFszk+iarOGTBJHj/EWjYErqfGlQ87XwC1KmXgHchIolQbiIws9OBe4E2wJdAR4LO356JDU0qU1b2Fqa8/ynXjTicKe9/yqCMFvuVDEr2CWR1z+DEpxfx8Ln3MmjfIlj5D1g2Az6YDAc1hSNODo4UDj8R6jVK4DsTkQMVqdYQMAL4l7v3MbPhBDWDrk5GgCXp1ND++96HeAU6eiN1OOfthnVzYeUs+OgV2PUV1KkfJIMjTwuSQ4ODE/EWRaQcBzQfgZktcPd+YULo4+4FZvYfdx+QiGDLo0Sw/yp01dCB2pcPn2bBypeCo4UdG6BWneC00ZGnBSOZmxyWmG2LyPccaCL4F3Am8N9AS4LTQ/3dfXAlxxmJEkE1VFAQTJizclZw+2otYEEHc/fR0H0UHNxZI5lFEuhAE0EjYBfBpaYXAc2Aqe6+tbIDjUKJoJpzhy9XMv/VyfTInUujr8OB5A1b8nXzI1lXtyvHDBwObTKDq5KUHBIiJUeJklIHetXQ1cDz7p4DTK7UyCT9mMFhPcgbciNDnj6Rx89sQd+9C9m06n22Zc8n0+bBJ5OCZRu2gNa9oXVmkBhaZwY1kpQcDljh2JLS+o0k/UQ5Ivg9MAb4imBOgenuvikJsZVKRwQ1R6ljGzo0hk3LYeMi2LAYNi6GL1dCQVjzqMEhQXIoTAxtMqF5RyWHCqiUsSVSbRzQEYG7/wH4g5n1As4D5ppZjrufWMlxSpopc2xDu77BrVDebvhy+XeJYcNiyHowJjkc/P0jh4M7FSUHnQYpXWWMLZGaIfKAMoJO4i+ArcChiQlH0knksQ1160PbvsGtUP6e8Mhh8XcJYt5DUJAXPF+/edGRw/BaGdw4tRY3XTiSwYe30mmQ0IGOLZGaI8qpoWsIjgRaAdOBabHzDiebTg3VDJUxtuF78vfAlyuKHzlsWl6UHHbQkG8bteejnU044ojutG7fBZq2DW9tglvdBpX1FoGqezSSkP0vVdqBXjV0J/Csuy9OQGz7TYmgZkjaB2T+3iA5bFzM4v+8zVcbsjm6yTe0KtgCu7d9f/mGLcKkEJsg2kKzmMf7kSyq6gduVU1QkjgVSgRm1tTdt5vZIaU97+5fVWKMkSkRSEWU2jHavgFs3wjbc2D7Btj+OeR+Ht7fELTv+vr7K2twyHdJoVnbEokjfFyvYfxt61u3JFlFO4ufBkYDCwnmGo69LMOBLpUWoUgClfwWPiijRczjw6Hl4WW/eO+3sGMj5MYki+2ff3c/Z35QSqOkBgcXJYXBTdvyYJt6vDZ3N386shOD9+6BtU3goCZBXaaDmgQTA9VtqKufJCVUhlpqvISfBsnbFZMkSh5ZfE7e1znU3RPhANpqlUgOJW/ltNdr/F2bSoJLCQfaRzCTYPzATHf/NgHx7RclAqlOCo9GHjqvB8e2PYiFH3/C/8xayH+d2I6jWtSCPTtgz/bwZ+yttLYdkLcz2obrNoqeTOo1gtp1oXY9qFX3u/u164X368bcD9tr1Sm+jI5kEqayvsgc6Mji+wiuGrrTzP4DTANecvfdkSMQSVNLc3KZcGEfjg3/iftmtuS6Jh35d04uR3WrwNFIwb7vksLeb+InjZLtOzcXb/eCynujRYkhTrL4XmKJskyEBFWrTpCIrFY5tzjLYOUvE2U9xZapnOSYjFHgkU8NmVltgnLUVwEj3b1ppUWxH3REIFIJ3CHv2zCh7AwG5+3bG97yYn7mfddeKcuEPwtKvG5fKa8rHBNSnR1QQvnu+V15BXyxYy+fdjyXX+UMqdAFBwc8Q5mZNQBOIzgyOAbVHBKp3syCU0JVedIg9xKJpYxksS8f8OAIp8xblOfLWybqemIelxrX/m+ngRewd2MuMz7O5+JhHSr9qrMoM5RNAwYC/wQeAua4V+YxpYhIKcy+OwVEFU5YSZCVvYVxKxdx8bAOCRkFHjcRmFkt4EPgQnffV2lbFRGRSOJf/lw5yaBWvCfDb/6nKgmISFUycW42WdlbirVlZW9h4tzsFEWUOIUXHBR+6A/OaMmEC/uwNCe30rYRNxGEZpvZOWa6PkxEqobCK2kKk0Hht+Ze7ZqlOLLKN3Zoxve++Q/OaFmppUCidBbfQHCCLt/MdhOMMPZUXTUkIlL4rVilOypHlPkImiQjEBGR/aH5FCpPlKuGTiit3d3frvxwRESi0XwKlSfKqaEbY+7XBwYQFKIbkZCIRETKkYwradJJlFNDp8U+NrP2wN0Ji0hEpBzxrqRRIth/+zNVZaEc4KjKDkREJKrSrpgZnNFSSaCCovQRPEgw/wAEl5tmAksSGJOIiCRRlCOC2Apv+cAz7v5uguIREZEki9JHUFRgzswOBtonNCIREUmqckcWm9kcM2sazl28BPibmd2X+NBERCQZopSYaObu24Gzgb+5e1/gxCgrN7ORZrbKzNaY2U1xlutvZvvM7MfRwhYRkcoSJRHUMbPWwBjgpagrDieyeQg4BegBXGBmPcpY7i7gtajrFhGRyhMlEdxO8CG9xt3nm1kX4OMIrxsQvmatu+8lmPf4jFKWuxaYAXwZMWYREalEUTqLnweej3m8FjgnwrrbAp/FPM4hmOCmiJm1Bc4iGKXcv6wVmdnVwNUAHTp0iLBpERGJKsoRQUWVVra65ATJfwHGlzffgbs/5u793L1fq1atKis+ERGhYiOLo8qh+KWm7YANJZbpBzwbTnXQEhhlZvnu/n8JjEtERGIkMhHMB7qaWWfgc+B84MLYBdy9c+F9M3sSeElJQEQkucpMBGZ2Q7wXunvcsQTunm9m4wg6mmsDk9x9uZmNDZ+fWIF4RUSkksU7IiickKYbQUfurPDxaUCkuQjc/RXglRJtpSYAd78syjpFRKRylZkI3P0PAGY2GzjG3XeEj28j5ioiERGp3qJcNdQB2BvzeC/QKSHRiIhI0kXpLH4K+I+ZvUhw+edZwN8TGpWIiCRNlAFlfzKzV4EhYdPl7r4osWGJiEiyRB1Q1hDY7u73AznhJaEiIlIDRClD/XtgPHBz2FQXmJLIoEREJHmiHBGcBZwO7ARw9w18d2mpiIhUc1ESwV53d8I6QWbWKLEhiYhIMkVJBM+Z2aNAczO7CvgX8ERiwxIRkWSJctXQPWZ2ErCdYJTx79z99YRHJiIiSVFuIjCzu9x9PPB6KW0iIlLNRTk1dFIpbadUdiAiIpIa8aqPXgP8HOhiZktjnmoCvJvowEREJDninRp6GngV+G/gppj2He7+VUKjEhGRpIlXfTQXyAUuADCzQ4H6QGMza+zunyYnRBERSaQoI4tPM7OPgXXAXGA9wZGCiIjUAFE6i+8ABgGrw6klf4j6CEREaowoiSDP3bcCtcyslru/BWQmNiwREUmWKPMRbDOzxgTTU041sy+B/MSGJSIiyRLliOAMYBfwK+CfQDbBvMUiIlIDRCkxsRPAzJoC/0h4RCIiklRRSkz8DLid4KigADCCSqRdEhuaiIgkQ5Q+gl8DPd19S6KDERGR5IvSR5ANfJvoQEREJDWiHBHcDGSZ2fvAnsJGd78uYVGJiEjSREkEjwJvAh8S9BGIiEgNEiUR5Lv7DQmPREREUiJKH8FbZna1mbU2s0MKbwmPTEREkiLKEcGF4c+bY9p0+aiISA0RZUBZ52QEIiIiqRFvhrIR7v6mmZ1d2vPu/kLiwhIRkWSJd0QwlOBqodLqCjmgRCAiUgPEm6Hs9+Hd2919XexzZqbTRSIiNUSUq4ZmlNI2vbIDERGR1IjXR9Ad6Ak0K9FP0JRg7uJymdlI4H6gNvCEu99Z4vmLgPHhw2+Aa9x9SfTwRUTkQMXrI+gGjAaaU7yfYAdwVXkrNrPawEPASUAOMN/MZrn7ipjF1gFD3f1rMzsFeAwYuF/vQEREDki8PoKZwEwzO9bd51Vg3QOANe6+FsDMniWY5KYoEbh7Vszy7wHtKrAdERE5AFH6CM4ys6ZmVtfM3jCzLWZ2cYTXtQU+i3mcE7aV5Urg1dKeCEc2LzCzBZs3b46waRERiSpKIviRu28nOE2UAxwB3BjhdVZKm5e6oNlwgkQwvrTn3f0xd+/n7v1atWoVYdMiIhJVlBITdcOfo4Bn3P0rs9I+478nB2gf87gdsKHkQmbWC3gCOMXdt0ZZsYiIVJ4oRwT/MLOPgH7AG2bWCtgd4XXzga5m1tnM6gHnA7NiFzCzDgQD0y5x99X7F7qIiFSGKLWGbjKzu4Dt7r7PzL4l6PQt73X5ZjYOeI3g8tFJ7r7czMaGz08Efge0AB4OjzLy3b1fxd+OiIjsL3Mv9bQ9ZvYbd787vH+uuz8f89yf3f23SYqxmH79+vmCBQtSsWkRkWrLzBaW9UU73qmh82Pu31ziuZEHHJWIiFQJ8RKBlXG/tMciIlJNxUsEXsb90h6LiEg1Fa+zuLeZbSf49t8gvE/4OFKtIRERqfrilZioncxAREQkNaKMIxARkRpMiUBEJM0pEYiIpDklAhGRNKdEICKS5pQIRETSnBKBiEiaUyIQEUlzSgQiImlOiUBEJM0pEYiIpDklAhGRNKdEICKS5pQIRETSnBKBiEiaUyIQEUlzSgQiImlOiUBEJM0pEYiIpDklAhGRNKdEICKS5pQIRETSnBKBiEiaUyIQEUlzSgQiImlOiUBEJM0pEYiIpDklAhGRNKdEICKS5pQIRETSXEITgZmNNLNVZrbGzG4q5XkzswfC55ea2TGJjEdERL4vYYnAzGoDDwGnAD2AC8ysR4nFTgG6hrergUcSFY+IiJQukUcEA4A17r7W3fcCzwJnlFjmDODvHngPaG5mrRMYk4iIlFAngetuC3wW8zgHGBhhmbbAxtiFzOxqgiMGgG/MbFXEGFoCW6IGnERVNS5QbBVRVeMCxVYRVTUuOLDYOpb1RCITgZXS5hVYBnd/DHhsvwMwW+Du/fb3dYlWVeMCxVYRVTUuUGwVUVXjgsTFlshTQzlA+5jH7YANFVhGREQSKJGJYD7Q1cw6m1k94HxgVollZgE/Ca8eGgTkuvvGkisSEZHESdipIXfPN7NxwGtAbWCSuy83s7Hh8xOBV4BRwBrgW+DySg5jv08nJUlVjQsUW0VU1bhAsVVEVY0LEhSbuX/vlLyIiKQRjSwWEUlzSgQiImmuRiaC8kpbJDmW9mb2lpmtNLPlZnZ92H6bmX1uZovD26gUxLbezD4Mt78gbDvEzF43s4/DnwenIK5uMftlsZltN7NfpmqfmdkkM/vSzJbFtJW5n8zs5vBvb5WZnZzkuP7HzD4KS7a8aGbNw/ZOZrYrZt9NTFRccWIr8/eXrH0WJ7ZpMXGtN7PFYXvS9lucz4rE/625e426EXRMZwNdgHrAEqBHCuNpDRwT3m8CrCYouXEb8OsU76v1QMsSbXcDN4X3bwLuqgK/zy8IBsOkZJ8BJwDHAMvK20/h73YJcBDQOfxbrJ3EuH4E1Anv3xUTV6fY5VK0z0r9/SVzn5UVW4nn7wV+l+z9FuezIuF/azXxiCBKaYukcfeN7v5BeH8HsJJg9HRVdQYwObw/GTgzdaEA8EMg290/SVUA7v428FWJ5rL20xnAs+6+x93XEVwRNyBZcbn7bHfPDx++RzA2J+nK2GdlSdo+Ky82MzNgDPBMorZfljifFQn/W6uJiaCsshUpZ2adgD7A+2HTuPAQflIqTsEQjOKebWYLwzIeAId5OJYj/HloCuKKdT7F/ylTvc8KlbWfqtLf3xXAqzGPO5vZIjOba2ZDUhRTab+/qrTPhgCb3P3jmLak77cSnxUJ/1uriYkgUtmKZDOzxsAM4Jfuvp2g0moGkElQW+neFIR1nLsfQ1AF9hdmdkIKYiiTBQMRTweeD5uqwj4rT5X4+zOzW4B8YGrYtBHo4O59gBuAp82saZLDKuv3VyX2WegCin/xSPp+K+WzosxFS2mr0H6riYmgypWtMLO6BL/Yqe7+AoC7b3L3fe5eADxOAg+Fy+LuG8KfXwIvhjFssrACbPjzy2THFeMU4AN33wRVY5/FKGs/pfzvz8wuBUYDF3l4Mjk8fbA1vL+Q4HzyEcmMK87vL+X7DMDM6gBnA9MK25K930r7rCAJf2s1MRFEKW2RNOE5x78CK939vpj22HLbZwHLSr42wXE1MrMmhfcJOhmXEeyrS8PFLgVmJjOuEop9O0v1PiuhrP00CzjfzA4ys84Ec238J1lBmdlIYDxwurt/G9PeyoI5QjCzLmFca5MVV7jdsn5/Kd1nMU4EPnL3nMKGZO63sj4rSMbfWjJ6w5N9IyhbsZoge9+S4liOJzhcWwosDm+jgKeAD8P2WUDrJMfVheCKgyXA8sL9BLQA3gA+Dn8ekqL91hDYCjSLaUvJPiNIRhuBPIJvYVfG20/ALeHf3irglCTHtYbgvHHh39rEcNlzwt/zEuAD4LQU7LMyf3/J2mdlxRa2PwmMLbFs0vZbnM+KhP+tqcSEiEiaq4mnhkREZD8oEYiIpDklAhGRNKdEICKS5pQIRETSnBKBVAtmdktYkXFpWAVyYKpjOhBm9qSZ/TgB6/1tzP1OsRU2RcqiRCBVnpkdSzBS9hh370Uw8Oez+K9KW78tfxGR4pQIpDpoDWxx9z0A7r7Fw/IYZtY3LAa20MxeixmK39fMlpjZPAtq9C8L2y8zswmFKzazl8xsWHj/R+HyH5jZ82HNl8J5G/4Qtn9oZt3D9sZm9rewbamZnRNvPWWJ8x7mmNldZvYfM1tdWPDMzBqa2XPhNqeZ2ftm1s/M7gQahEdMhTWGapvZ4+HR1Gwza1ApvxGpUZQIpDqYDbQPPwwfNrOhUFSX5UHgx+7eF5gE/Cl8zd+A69z92CgbMLOWwH8BJ3pQiG8BQZGxQlvC9keAX4dttwK57n50eKTyZoT1lNxuvPcAwdwCA4BfAr8P234OfB1u849AXwB3vwnY5e6Z7n5RuGxX4CF37wlsIxgpK1JMnVQHIFIed//GzPoSlAgeDkyzYOa5BcBRwOtBmRZqAxvNrBnQ3N3nhqt4iqCAXTyDCCb6eDdcVz1gXszzhQXAFhIUJoPgFNX5MXF+bWajy1lPSd1Kew9lbLdTeP944P5wm8vMbGmc9a9z98WlrEOkiBKBVAvuvg+YA8wxsw8Jim8tBJaX/NZvwfSMZdVOyaf4kXD9wpcBr7v7BWW8bk/4cx/f/d9YKdspbz0lGaW8hwjbjWpPzP19gE4Nyffo1JBUeRbMYdw1pikT+ISg0FarsDMZM6trZj3dfRuQa2bHh8tfFPPa9UCmmdUys/Z8Vwr5PeA4Mzs8XFdDMyuv3PBsYFxMnAdXYD2lvodytvtvglm0MLMewNExz+WFp5tEIlMikOqgMTDZzFaEp0F6ALd5MBXpj4G7zGwJQbXGweFrLgceMrN5wK6Ydb0LrCOognkPQUVJ3H0zcBnwTLiN94Du5cR1B3CwmS0Ltz98f9dTznsoy8MEyWMpQcnppUBu+NxjwNKYzmKRcqn6qNR4Fkz795K7H5XqWCqDBfXx67r7bjPLIChNfESYVET2m/oIRKqfhsBb4SkgA65REpADoSMCEZE0pz4CEZE0p0QgIpLmlAhERNKcEoGISJpTIhARSXP/H1h+VCcMEYyPAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(seq_lens, est_pr_survival, 'x', label='Measured')\n",
    "plt.plot(seq_lens, model(seq_lens, *p_opt), label='Curve fit')\n",
    "plt.plot()\n",
    "plt.ylim(0, 1)\n",
    "plt.xlabel('Sequence length')\n",
    "plt.ylabel('Estimated survival probability')\n",
    "plt.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "While this gives us qualitative agreement with our RB data, there are two problems left by using curve fitting:\n",
    "\n",
    "- The curve that best minimizes residual errors may or may not be the curve that gives us the most accurate value for $p$.\n",
    "- Because all estimation is done in post-processing, we cannot use measurements online to choose ideal sequence lengths."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Thankfully, we can use Q# to solve both of these problems."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Particle filtering"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, let's see how we can improve our estimation given RB results. Essentially, curve fitting answers the question of what parameters to a model minimize the distance between the model and the observed data (in particular, where \"distance\" is formalised as \"mean squared error\"). This is a subtlely different question than that of what parameters minimize the average error with respect to the true parameters; the former is concerned with distance between datasets, while the latter is concerned with the distance between model parameters.\n",
    "\n",
    "These two questions coincide when every data point follows the same normal distribution (homoskedastic normally distributed errors), but since that does not hold in randomized benchmarking, we'll turn to Bayesian inference instead. In particular, we use a technique known as _particle filtering_, where each different hypothesis about $p$, $A$, and $B$ is represented by a weight $w_i$ associated with that hypothesis,\n",
    "$$\n",
    "    \\Pr(p, A, B) \\approx \\sum_i w_i \\delta(p - p_i) \\delta(A - A_i) \\delta(B - B_i).\n",
    "$$\n",
    "This allows us to write down everything we currently know about our parameters in terms of a list $\\{(w_i, p_i, A_i, B_i)\\}_{i=0}^{N - 1}$.\n",
    "\n",
    "Using the QInfer library to implement particle filtering, we can go on and analyze the results above."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "df = pd.DataFrame({'m': seq_lens, 'counts': 100 * np.array(est_pr_survival), 'n_shots': 100, 'interleaved': False})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([0.96547397, 0.47530413, 0.49100553]),\n",
       " array([[ 4.06717136e-05,  9.36237953e-05, -9.99758858e-05],\n",
       "        [ 9.36237953e-05,  9.28988116e-04, -6.22238620e-04],\n",
       "        [-9.99758858e-05, -6.22238620e-04,  5.96174614e-04]]))"
      ]
     },
     "execution_count": 36,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "est, cov, extra = qi.simple_est_rb(data=df.to_records(), return_all=True)\n",
    "est, cov"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using the data we get back from QInfer, we can plot our uncertianty about $p$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x1a220d88640>]"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAEICAYAAABGaK+TAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAp/UlEQVR4nO3deXxc1X338c9vFu27NJJlS7a8W94XYXZCbIwhJBhooYTQuGlS2jSlWZsH2teTZumT0KQpSdssJaSJGwghSSEmS8HGLHZIMMgbeJPlfdO+76OZOc8fc02Eke2RNDN37szv/XqJmbkz1v0dZvTV1bnnniPGGJRSSjmPy+4ClFJKjY8GuFJKOZQGuFJKOZQGuFJKOZQGuFJKOZQnnjsrKSkxVVVV8dylUko53o4dO1qNMb7zt8c1wKuqqqitrY3nLpVSyvFE5MRo27ULRSmlHEoDXCmlHEoDXCmlHEoDXCmlHEoDXCmlHEoDXCmlHEoDXCmlHCqu48CVUmHGGDbvb6K5Z4h+f4DhoOGelVMpzE6zuzTlIBrgStmg9kQH9/1ox9u2ed3CfdfNtKki5UQa4ErZ4PdH2hCBFz59PaW56dz48Fb2nO6yuyzlMBrgStngtWPtzJuUx/SSbACWVObzxulOe4tSjqMnMZWKs+FgiB0nOrh8etFb2xZXFHCqfYD2Pr+NlSmn0QBXKs7ePNPFwHCQlW8L8HwAPQpXY6IBrlScvXasHeBtAb5oSj4isOeU9oOryGmAKxVn24+2MdOXTUlO+lvbcjO8zPTl6BG4GhMNcKXiKBgy1B7vYOX04nc8t7ginz2nuzDG2FCZcqJLBriIzBWR3SO+ukXkEyJSJCKbRaTeui2MR8FKOdmBhm56hgJcMaPoHc8tqSigtXeIhq5BGypTTnTJADfG1BljlhpjlgIrgH7gaeABYIsxZjawxXqslLqI7Vb/92VVowR4ZQGgJzJV5MbahbIaOGKMOQGsAzZY2zcAt0WxLqWS0mvH2qgsymRyQeY7nqsuz8XrFnbriUwVobEG+N3AE9b9MmNMA4B1WzraPxCR+0SkVkRqW1paxl+pUg5njOG1Y+2srHpn/zdAusfNvEl5bzsCr2vsod8fiFOFymkiDnARSQNuBX42lh0YYx4xxtQYY2p8vncsqqxUyjjU1EtH/zCXj9L/fc7iinzePN1FKGT49kuHWfuNrfz14zv1xKYa1ViOwG8GdhpjmqzHTSJSDmDdNke7OKWSydO7zuASuG72hQ9kllQW0DMU4M83vM5Xn61jblkuL9W1sHH32ThWqpxiLAH+fv7QfQLwDLDeur8e2BitopRKNkOBID+rPcXq6jIm5Wdc8HVLKgoAeKmuhU/eMIdf/+01LK0s4Au/3Edb71CcqlVOEVGAi0gWsAZ4asTmh4A1IlJvPfdQ9MtTKjk8u7eRtj4/914x7aKvm1Wawwcun8q37lnOx2+Yjcft4qt/vJjeoQBf+tX+OFWrnCKi2QiNMf1A8Xnb2giPSlFKXcLjr55kWnEW184quejr3C7h/92+6G3b5pTl8tHrZ/FvW+q5bdkUrp876ngBlYL0SkylYqyusYfXjrdzz8qpuFwyru/xsXfPxJebztO7zkS5OuVkGuBKxdjj20+Q5nFxZ03luL9HusfNLF8Op9r7o1iZcjoNcKViqG8owFM7z3DLonKKJrjeZWVRJqc7BqJUmUoGGuBKxdDzB5roHQpwz+VTJ/y9KgqzaO4ZYnA4GIXKVDLQAFcqhrYfayc33cPyqROf662yKHz5/ZlOPQpXYRrgSsXQjuMdLJtWiHucJy9HqijMAtB+cPUWDXClYqSrf5i6ph5qpkVnpuVKK8C1H1ydowGuVIzsPNkBQE1VdAK8NDedNLeLUx16BK7CNMCVipHaE+24XcJSa57viXK5hCmFOhJF/YEGuFIxUnu8gwWT88hKi+iC54hUFGZyWvvAlUUDXKkY8AdC7D7VyYoo9X+fU1GYpUfg6i0a4ErFwL6zXQwFQqMunTYRFYWZtPX56RvSRR6UBrhSMbHjhHUCM8pH4JVF4ZEoOhZcgQa4UjFRe7yDyqJMSvMuPPf3eFQWhi/m0bHgCjTAlYo6Ywy1J9qpmRbd7hP4w8U82g+uQANcqag70dZPa68/auO/RyrJSSPD69IjcAVogCsVdbVv9X9H/whcRKgozNKLeRSgAa5U1NUebycvw8Ps0pyYfP9KvZhHWSJdE7NARH4uIgdF5ICIXCkiRSKyWUTqrdvo/72olAO9frydmqqica++cykVhVnahaKAyI/Avwk8a4yZBywBDgAPAFuMMbOBLdZjpVJaW+8QR1r6YtL/fU5lUSbdgwG6BoZjtg/lDJcMcBHJA64Dvg9gjPEbYzqBdcAG62UbgNtiU6JSznFu/He0L+AZ6Q8jUfQoPNVFcgQ+A2gBfiAiu0TkURHJBsqMMQ0A1u2oS2WLyH0iUisitS0tLVErXKlEVHuigzS3i0VT8mO2D51WVp0TSYB7gOXAd4wxy4A+xtBdYox5xBhTY4yp8fl84yxTKWd4/Xg7iyvyyfC6Y7aPCr2YR1kiCfDTwGljzHbr8c8JB3qTiJQDWLfNsSlRKWcY8AfZe6aLmhh2nwAUZHnJSffoEbi6dIAbYxqBUyIy19q0GtgPPAOst7atBzbGpEKlHGL3qU6Gg4aV02M7IEtEqCrJoq6xJ6b7UYkv0omK7wceF5E04CjwIcLh/1MR+TBwErgzNiUq5Qy1x9sBWDE1tkfgANfM8vHotqP0DA6Tm+GN+f5UYopoGKExZrfVj73YGHObMabDGNNmjFltjJlt3bbHulilEtnrJzqYW5ZLflbsA3V1dSmBkGFbfWvM96USl16JqVQUBEOGnSc6Yjr+e6RllQUUZHnZckBPPaUyDXClouBgYze9Q4GYjv8eyeN2cf0cHy/VNRMMmbjsUyUeDXClouD1Y+EexHgdgQOsqi6jrc/P7lOdcdunSiwa4EpFwe+OtDGlIJMpBZlx2+e7Zvtwu4QXDjbFbZ8qsWiAKzVBgWCI3x9t49rZJYjEZgKr0eRneamZVqj94ClMA1ypCXrjTBc9gwGunlUS932vri7lYGOPrpGZojTAlZqgV6yhfHYE+Kp5ZQC8cFCPwlORBrhSE/Tbw60smJxHUXZa3Pc905fNtOIsntp5WkejpCANcKUmoG8owM6THVwzO/5H3xC+rP5v3j2LXSc7eXjzIVtqUPbRAFdqAl473s5w0HCNDd0n59xZU8ldNRX8x4uHdURKitEAV2oCXqlvJc3jitsFPBfyxXULmV+exyef3KPTzKYQDXClJuC3h1u5rKowpvN/RyLD6+Y79y4nZAyf/tkeW2tR8aMBrtQ4NfcMcrCxx5bRJ6OZVpzNp9bM4bVj7ew62WF3OSoONMCVGqffH2kD4NpZibPS1J01leSme/jBK8ftLkXFgQa4UuP04sFmCrO8zJ+cZ3cpb8lJ93DXZZX85s0GGrsG7S5HxZgGuFLjMDgc5PkDzdw4fxJuV/wun4/E+iurCBrDj149bncpKsY0wJUah9/Wt9I7FODmRZPsLuUdphZnsaa6jB9vP8ngcNDuclQMRRTgInJcRN4Ukd0iUmttKxKRzSJSb93Gbx5NpWz2mzcbyM/0JswJzPN96OrpdPQPs3H3GbtLUTE0liPwdxtjlhpjaqzHDwBbjDGzgS3WY6WS3lAgyOb9TayZX4bXnZh/xF4xo4jq8jx++LsTdpeiYmgin751wAbr/gbgtglXo5QDvHK4lZ6hALcsKre7lAsSEe6qqeBAQzfHW/vsLkfFSKQBboBNIrJDRO6ztpUZYxoArNvSWBSolN2GAkGGAn/oS/71G43kZngStvvknDXzwzMVbt6vl9cnq0gD/GpjzHLgZuBjInJdpDsQkftEpFZEaltaWsZVpFJ2+uhjO7n6oRfYtK8RfyDE5v2NrJlfRponMbtPzqkozGJ+eR6b9jfaXYqKkYg+gcaYs9ZtM/A0sBJoEpFyAOt21AmJjTGPGGNqjDE1Pl/iXPCgVCSOtPTywsFmhoZD3PejHXzg0VfpHkzs7pOR1swvY8eJDtp6h+wuRcXAJQNcRLJFJPfcfeBGYC/wDLDeetl6YGOsilTKLk9sP4nHJTz3yeu4f9Usdp7sJDfdY9v0sWO1Zn4ZIQNbdMGHpOSJ4DVlwNPWWn8e4MfGmGdF5HXgpyLyYeAkcGfsylQq/gaHg/xsx2nWLpjE5IJMPn3jXG5eWM7AcJB0j72TV0VqweQ8phRksmlfE3fVVNpdjoqySwa4MeYosGSU7W3A6lgUpVQi+PUbDXQNDPOBK6a+tS2RLpuPhIiwZn4ZP3n9JAP+IJlpzvjFoyKT2GdhlLLRY9tPMMOXzZUziu0uZULWzC9jcDjEtnodRJBsNMCVGsW+s13sOtnJBy6fhtV96FgrpxeRl+Fhkw4nTDoa4EqN4sfbT5LucfFHy6fYXcqEed0u3j2vlBcONmOMLnycTDTAlTpPKGR4dm8jaxdMoiAr/ivNx8KyygLa+/y06HDCpKIBrtR59jd009bn5/q5yXPdwgxfDgBHW/Sy+mSiAa7UeV4+FD7Zd+3s5Anw6SXZABzTeVGSiga4UufZeqiF+eV5+HLT7S4laiYXZJLmcWmAJxkNcKVG6B0KsONEB9fOccaVlpFyu4Tpxdkcbem1uxQVRRrgSo3w6pE2AiHDu5Ko++Sc6SXZHNUj8KSiAa7UCFvrW8j0ullRlXwLTE33ZXOyrZ9AMGR3KSpKNMCVGmHroRaunFnsmLlOxmJGSTaBkOFUx4Ddpago0QBXynKyrZ/jbf1c55CZBsdqhu/cSBTtB08WGuBKWV625gq5dk7y9X8DTC/RseDJRgNcKcvLdS1MKchkhjVmOtkUZadRkOXVoYRJRANcKeC7Lx/h+QNN3LK43PGTV13M9JJsPQJPIhrgKqUZY/j6pjoe+t+DvG/JZP5u7Vy7S4qp6SXZegSeRDTAVUr76nN1/PsLh7n7skq+8SdL8bqT+0dipi+Hxu5B+oYCdpeioiC5P61KXYQ/EOJ7W49y65LJfOWORbhdydt1co7OiZJcIg5wEXGLyC4R+ZX1uEhENotIvXWbfFc+qKR2pKWXQMhww/yypO73HkkDPLmM5Qj848CBEY8fALYYY2YDW6zHSjnGoaYeAOZNyrW5kvjRAE8uEQW4iFQAtwCPjti8Dthg3d8A3BbVypSKsYONPXjdQlVxcg4bHE2G182Ugkyd1CpJRHoE/g3gs8DISRTKjDENANZt6Wj/UETuE5FaEaltadFFVVXiONTYw4ySHNI8qXUqSEeiJI9LfnJF5L1AszFmx3h2YIx5xBhTY4yp8fmS8wo35UwHG3uYm0LdJ+fM8IVnJdT1MZ0vkkOPq4FbReQ48BNglYg8BjSJSDmAddscsyqVirKewWHOdA6kZIDPLsulZzDAyfZ+u0tRE3TJADfGPGiMqTDGVAF3Ay8YY+4FngHWWy9bD2yMWZVKRVl9c7gPeG5Z6gX4NbPCk3VtPaRdmk43kc6/h4A1IlIPrLEeK+UIdY3hESipeAReVZxFZVHmW2t/KufyjOXFxpiXgJes+23A6uiXpFTs1TX2kJUWHpGRakSEd83x8fTOM/gDoZQ7iZtM9J1TKamusYc5Zbm4UuDqy9FcN9tHnz/IjhMddpeiJkADXKWkQ009Kdn/fc6VM4vxuES7URxOA1ylnNbeIdr6/MxJwf7vc3IzvKyYVqgnMh1OA1ylnHMnMFPpEvrRXDfHx/6Gbpp7Bu0uRY2TBrhKOecCfE4Kd6EAvMtaOm7boVabK1HjpQGuUk5dYw/F2Wn4ctPtLsVW88vzKMlJY2u9dqM4lQa4Sjl1TT0pf/QN4HIJ1872sa2+lVBIL6t3Ig1wlVJCIUN9U2rOgTKa6+f6aO/zs/t0p92lqHHQAFcppbF7kD5/kFmlOXaXkhCun1OKxyU8t6/R7lLUOGiAq5RyxJoHe6ZPAxwgP8vLlTOL2bSvSWcndCANcJVSjraE58Ge6UudRRwu5cYFkzjW2vfWBF/KOTTAVUo52tJLTron5UegjHTj/DIAntur3ShOowGuUsqRlj5m+rJTZhHjSJTlZbBsagHP7dcAdxoNcJVSjrb0MkP7v99h7YJJ7D3TzekOXeTBSTTAVcro9wc42zXIjBLt/z7f2gWTANi0r8nmStRYaICrlPHWCUwdQvgO00uymVOWo8MJHUYDXKWMo9ZK7DN0BMqo1i6YxOvH22nv89tdioqQBrhKGUeaexGBqmIN8NFcN8dHyMCuk7rIg1NcMsBFJENEXhORPSKyT0S+YG0vEpHNIlJv3RbGvlylxu9oax8VhZlkeN12l5KQzk2ve6Ch2+ZKVKQiOQIfAlYZY5YAS4GbROQK4AFgizFmNrDFeqxUwjrS3MuMEu3/vpDcDC+VRZkcaOixuxQVoUsGuAk7d4mW1/oywDpgg7V9A3BbLApUKhpCIcOx1j69hP4SqiflcaBRj8CdIqI+cBFxi8huoBnYbIzZDpQZYxoArNvSC/zb+0SkVkRqW1p03mFlj8buQQaGg3oC8xKqy/M43trHgD9odykqAhEFuDEmaIxZClQAK0VkYaQ7MMY8YoypMcbU+Hy+cZap1MToJFaRqS7PI2TCc6arxDemUSjGmE7gJeAmoElEygGs2+ZoF6dUtOgkVpGpLtcTmU4SySgUn4gUWPczgRuAg8AzwHrrZeuBjTGqUakJO6KTWEWksjCL7DS3BrhDeCJ4TTmwQUTchAP/p8aYX4nI74GfisiHgZPAnTGsU6kJOaqTWEXE5RLmledxUEeiOMIlA9wY8wawbJTtbcDqWBSlVLQdbenl8hnFdpfhCNXluWzcfRZjjP7CS3B6JaZKek3dgzqJ1RhUl+fRMxjgdMeA3aWoS9AAV0ltKBDko4/tINPr5uZFk+wuxxHmTcoD9ESmE2iAq6T2+Wf2s/NkJ1+7czGzSnUl+kjMm5SLCHpFpgNogKuk9fj2Ezzx2kk+ev1M3rt4st3lOEZ2uodpRVkc1CsyE14ko1CUchR/IMR/vnyEb26p511zfHzmxrl2l+Q41eV52oXiABrgKqnsONHBg0+9waGmXm5ZXM6Xb1+E26UjKcaqujyPZ/c10jcUIDtdYyJR6TujksarR9u453uvMikvg++vr2F1dZndJTlWdXkexoRPZNZUFdldjroADXCVFIIhwxd+uZ/y/Eye/cS15GZ47S7J0VZMK0QEfnu4VQM8gelJTJUUflZ7igMN3Tz4nnka3lFQlJ3GssoCXjioUxwlMg1w5Xg9g8P8y6Y6aqYVcsuicrvLSRqr5pXyxukumnsG7S5FXYAGuHK8b790hNZeP59733y99DuKVs0Ln0N4qU7n8U9UGuDK0U539PP9bce4Y/kUFlcU2F1OUqkuz6U8P4MXDmg3SqLSAFeO9uTrpwiEQjrWOwZEhHfPK2VbfQv+QMjuctQoNMCVYxlj2Lj7LFfPKmFyQabd5SSlVXNL6fMHee1Yu92lqFFogCvH2nmyk5Pt/axbOsXuUpLWVbOKSfO4dDRKgtIAV461cfcZ0j0u1i7QC3ZiJSvNw1Uzi3mxTgM8EWmAK0caDob41RsN3DC/TMd9x9iqeaUca+3jqLUwtEockayJWSkiL4rIARHZJyIft7YXichmEam3bgtjX65SYdvqW2jv83Obdp/E3OrqMtwu4Ye/O253Keo8kRyBB4BPG2OqgSuAj4nIfOABYIsxZjawxXqsVFz8YtdZCrK8vGuOz+5Skt6UgkzuWTmVx7efpL5J5whPJJcMcGNMgzFmp3W/BzgATAHWARusl20AbotRjUq9Te9QgE37G7llUTlpHu0FjIdP3DCbrDQ3X/7NAbtLUSOM6dMvIlWEFzjeDpQZYxogHPJAadSrU2oUm/Y1Mjgc4rZl2n0SL8U56dy/ahYv1rWw9ZBemZkoIg5wEckB/gf4hDEm4pneReQ+EakVkdqWFn3j1cQ9vesMlUWZrJiqp13iaf1VVUwtyuKffr2fQFAv7EkEEQW4iHgJh/fjxpinrM1NIlJuPV8OjDrOyBjziDGmxhhT4/Npf6WamMauQX57uJXbl07BpQs1xFW6x82DN8/jUFMvT9aesrscRWSjUAT4PnDAGPOvI556Blhv3V8PbIx+eUq93cbdZzAGbl9eYXcpKemmhZOomVbIw5vr6RsK2F1OyovkCPxq4E+BVSKy2/p6D/AQsEZE6oE11mOlYsYYw1M7z7BsagHTS7LtLicliQh/f0s1rb1DPLL1qN3lpLxLrshjjPktcKG/VVdHtxylLmx/Qzd1TT186baFdpeS0pZPDc+7/sjWo3zg8qmU5mXYXVLK0jFYyjGe2nkGr1t4ry7aYLvP3jSXQCjEw88fsruUlKYBrhwhEAyxcfdZVs0rpTA7ze5yUt604mzuvWIaT75+ikN6cY9tNMCVI2yrb6W1d4jbl+nJy0Txt6tmk+F16yX2NtIAVwkvFDI8/PwhyvLSefc8HYqaKAqz03j33FI2728iFDJ2l5OSNMBVwvv5ztO8cbqLB26eR7rHbXc5aoQbF5TR0jPErlMddpeSkjTAVULrGRzmq8/WsXxqgc48mIBWzSslze3i2b2NdpeSkjTAVUL79xcO09Y3xOdvXaArzieg3AwvV88q5tl9jRij3SjxpgGuEtbRll5+8Mox7lxRoSvOJ7C1CyZxqn2AAw06GiXeNMBVwvrac3Wke9x8Zq2uOJ/Ibphfhkvg2X3ajRJvGuAqIe0728X/7m3kz6+ZTmmuXumXyEpy0rmsqojntB887jTAVUJ6eHM9uRkePnzNdLtLURFYu2ASdU09HGvts7uUlKIBrhLOG6c7ef5AE39x7QzyM3XBYidYu3ASAM/sPmtzJalFA1wlnIc3H6Igy8uHrq6yuxQVoSkFmdxQXcq3XzrM4WY9mRkvGuAqobx+vJ0X61q477oZ5Gbo0beTfPmORWSne/jEk7vxB3TFnnjQAFcJYe+ZLj7xk128/5FXKclJZ/2VVXaXpMaoNDeDr9yxiL1nuvmGzlIYF5ecD1ypWBoKBHngf97k6V1nyE5z88Erq/jwtdPJTtePphOtXTCJuy+r5DsvH+Fdc3xcPqPY7pKSmv6UKNv0DQX4q8d2sK2+lftXzeIjetIyKfzf987n1aNtfOzHu/jFx66iojDL7pKSlnahKFt09Pm559Ht/O5IG/9y5xI+feNcDe8kkZ3u4XsfrGEoEOQjG2rpGRy2u6SkFcmixv8lIs0isnfEtiIR2Swi9dZtYWzLVMnmYz/eyYGGbr577wr+eIXO8Z1sZpfl8p0PrKC+uZf7n9hFIKgnNWMhkiPwHwI3nbftAWCLMWY2sMV6rFRE9p7p4ndH2vjMjXNYM7/M7nJUjFwzu4QvrVvIS3UtfH2zntSMhUsGuDFmK9B+3uZ1wAbr/gbgtuiWpZLZht8dJ9Pr5k8um2p3KSrG7rl8Ku9dXM5jr55gcDhodzlJZ7x94GXGmAYA67b0Qi8UkftEpFZEaltaWsa5O5Us2vv8bNxzljuWT9E+7xRx92VT6RkM8PyBJrtLSToxP4lpjHnEGFNjjKnx+XQ5rFT3xGsn8QdC/NlVVXaXouLkypnFTMrL4OmdZ+wuJemMN8CbRKQcwLptjl5JKlkNB0M89uoJrplVwuyyXLvLUXHidgnrlk3mpUMttPYO2V1OUhlvgD8DrLfurwc2Rqcclcw27WuioWuQ9Xr0nXLuWFZBMGT45R6d7CqaIhlG+ATwe2CuiJwWkQ8DDwFrRKQeWGM9VuqCugeH+fcX6qksymTVvAueMlFJau6kXBZMzuMp7UaJqkteiWmMef8Fnlod5VpUkurqH+aD/7Wdw829fPfeFbhdurZlKrp92RT+6dcHqG/q0S60KNErMVVMtff5uefRVznQ0MN3713BDTruO2XdunQybpfw1C49Co8WDXAVE0OBIE++fpLbvvUKh5t7eeSDGt6prjQ3g+vn+HjitZN09vvtLicpaICrqOrs9/Pdl49w7T+/yP/5nzfJzfDw33++kuvnar+3gs+snUv3wDBf36RXZkaDzkaoouJgYzc/fOU4v9h9hsHhEFfNLObrdy3hmlkliGiftwqrLs/jT6+Yxo9ePcH7V05l/uQ8u0tyNA1wNSGhkOE7Lx/h65vqSPO4uH3ZFD54ZRXV5fqDqUb3qTVz+eUbDXz+mX08+ZdX6C/4CdAAV+PW2e/nUz/dwwsHm3nfksl8ad0CCrLS7C5LJbj8LC9/t3YuDz71Js/sOcu6pVPsLsmxtA9cjUtj1yDv+4/fsq2+hS+uW8C/3b1Uw1tF7K6aShZNyeeLv9zPibY+u8txLA1wNWaDw0H+8rEdtPf6efIvr+SDV1bpn8FqTNwu4Zt3LyVkDH/2g9dp00vsx0UDXI2JMYbPbdzLnlOdfP2upSyfqmt5qPGZ4cvh0fWXcbZzgA9vqGXAr9PNjpUGuBqTx149wU9rT3P/qlnctHCS3eUoh1sxrZB/e/8y9pzu5P4ndhEMGbtLchQNcBWRAX+Qf91Uxxd+uZ/V80r55A1z7C5JJYm1Cybx+fct4PkDTXz1uYN2l+MoOgpFXZQ/EOK5fY185TcHONs1yK1LJvNPty/EpfOZqChaf1UV9c09/OfLR5lTmssf6TqpEdEAV29jjOHNM108v7+J7cfa2X2qk6FAiOryPL5x9zJWTi+yu0SVpP7xfQs42tLHg0+9SVVJNium6fmVSxFj4tfnVFNTY2pra+O2PxW5M50D/Oj3J/j1m2c51T6AS2DB5HwuqyriihlFrK4u01kEVcx19vtZ961XaOv1c9uyydyxvIJllQUpP8pJRHYYY2resV0DPLUN+IN89+Uj/OfWIwSChqtnlXDL4nJunF+m47qVLU609fH1TYfYtL+RweEQs0pz+OKtC7hqVondpdlGA1y9TUefn427z/DI1qOc7RrklsXlPHjzPCoKs+wuTSkAegaH+d+9jXz7xcMcb+vn/SsrefA91eRlpN5i2BrgSaytd4i6ph4ONfZwpnOAlp4hWnqHMAZKctIpyUknL9ODSwQB6pp62LSvCX8wxNLKAv7+PdXat60S1uBwkIc3H+J7245SlJ3GqnmlXDGjmMtnFDM5PyMlulc0wOPkZFs/Bxq7CQQNgVCINLeLyqIsppdkk5Xmpql7iCMtvdQ39XCgoYeDjd0ca+0jK81DQZaX/EwvRdlpFGankZ/ppXcwQEvPEK29QwwMBwkEDcOhEEPDIQaHg/T7gwwM/+ECiAyvC19uOr6cdESE1t4hWnqG6B9xkURBlpfblk7hrppKnQ1OOcYbpzv51ouHefVoO10Dw0D4szynLJd5k3JZNa+Uq2eV4HW/c3R0KGToGQqQm+5x5AiqmAS4iNwEfBNwA48aYy66NmYyBnj34DC7T3byyuFWthxs5nBz7wVfm+F1MTgceutxUXYa1eW5zCjJYSgQpGtgmI7+YTr6/HT0++noHyY3w4MvJx1fbjpZaW7cLsHjcpHudZHpdZPpdTMpP4N5k/KYU5aDLzd91COSUMhggJAxuEUc+SFWCsKf5YONPdSeaOdgYw91jT0cbOimzx+kMMvL6uoyBGjoGuRs1wDtfX66BoYxBoqz07hmdgnXzfZx3Rwfvtx0u5sTkagHuIi4gUOEFzU+DbwOvN8Ys/9C/2a8AW6MifqfScYYBodDnOkc4GR7H6c7BhgaEa6ZaW7yMr3kZXgQEQb8QQaHg7T3+WnsHqSxa5BDTT3UNfVgDHjdwuXTi1k1r5TLqopI97pwu4TB4SAn2vo51tpHe5+facVZzPLlMKv0wmGrlBqboUCQrYda+dUbZ3nhYDOZXjflBZlMzs+gJCedgiwvuRke9p/tZlt9K2194RWBllTks2peGZMLMujsH6ZzwE8gZMjwuMnwuvHlprNwSh4zfTmjHtnHy4UCfCLjwFcCh40xR60d/ARYB1wwwMfrcxv38dj2E3jdLrwuwetx4XG5SHMLIkIgFMIfCDFsdVuEQmAw4dd4XHjdLlwC57JywB+kzx8c92W7aR4X5fkZTCvO5uaF5ayYVsjSqQXkpI/+v3PB5PzxNl0pFYF0j5s188tYE8GyfaGQYX9DNy8ebOaFuma+seUQ545j3S7B7RL8gdDb/k2ax0VFQSajHW8Z6z+hUQ6GXSIg4duv3LGIy6qie65pIgE+BTg14vFp4PLzXyQi9wH3AUydOnVcO7p+ro+CLC/DQcNwMGR9GQLBEEFjSHOHQ9rrduFxS/hknUAgGA52f9BgjHnrTcpMc5Od7iYrzcPkggymFmVRWZRFVlr4f4cxhgF/kO7B4bf62jK84dfnZ3opzPLqkbNSDuVyCQun5LNwSj73r55Ne5+fvqEA+VlectPDf3GHQobBQJCznQPsO9vN3jNdnO0avPD3tAYIjIwFY3ir2xIDWWnuqLdlIl0odwJrjTEfsR7/KbDSGHP/hf5NMvaBK6VUrF2oC2UinTqngcoRjyuAsxP4fkoppcZgIgH+OjBbRKaLSBpwN/BMdMpSSil1KePuAzfGBETkb4DnCA8j/C9jzL6oVaaUUuqiJjQboTHmN8BvolSLUkqpMdAFHZRSyqE0wJVSyqE0wJVSyqE0wJVSyqHiOhuhiLQAJy7ykhKgNU7lxEOytQeSr03ansSm7QmbZozxnb8xrgF+KSJSO9rVRk6VbO2B5GuTtiexaXsuTrtQlFLKoTTAlVLKoRItwB+xu4AoS7b2QPK1SduT2LQ9F5FQfeBKKaUil2hH4EoppSKkAa6UUg4VtwAXkZtEpE5EDovIA6M8XygiT4vIGyLymogstLZnWI/3iMg+EflCvGq+mPG2Z8TzbhHZJSK/il/VFzaR9ojIcRF5U0R2i0hCrNgxwfYUiMjPReSgiBwQkSvjW/07TeDnZ671vpz76haRT8S9AaOY4Hv0SSsP9orIEyKSEd/q32mC7fm41ZZ9Y3p/wkuNxfaL8HSzR4AZQBqwB5h/3mu+BvyjdX8esMW6L0COdd8LbAeuiEfdsWjPiOc/BfwY+JWdbYlGe4DjQInd7YhiezYAH7HupwEFTm7Ped+nkfBFIY59jwgv53gMyLQe/xT4Mwe3ZyGwF8giPEPs88DsSPYbryPwtxZANsb4gXMLII80H9gCYIw5CFSJSJkJ67Ve47W+7D7zOu72AIhIBXAL8Gj8Sr6oCbUnAY27PSKSB1wHfN96zm+M6Yxb5aOL1vuzGjhijLnY1dDxMtE2eYBMEfEQDj67VwObSHuqgVeNMf3GmADwMnB7JDuNV4CPtgDylPNeswe4A0BEVgLTCC/Tdq67YTfQDGw2xmyPdcGXMKH2AN8APguESAwTbY8BNonIDgkvYm23ibRnBtAC/MDq4npURLJjX/JFTfT9Oedu4IkY1ThW426TMeYM8C/ASaAB6DLGbIp5xRc3kfdoL3CdiBSLSBbwHt6+XOUFxSvAR1vC/fyj6IeAQiuo7wd2AQEAY0zQGLOUcGNXnt+fbINxt0dE3gs0G2N2xLbEMZnQ+wNcbYxZDtwMfExErotVoRGaSHs8wHLgO8aYZUAf8I7+zDib6PuDhJc9vBX4WYxqHKuJ/AwVEj66nQ5MBrJF5N4Y1hqJcbfHGHMA+GdgM/As4aAPEIEJrcgzBpdcANkY0w18CEBEhHAf17HzXtMpIi8BNxH+rWWXibTnbuBWEXkPkAHkichjxhg7P4ATen+MMWet22YReZrwn5NbY1/2BU2kPVnA6RF/5f0c+wM8Gj8/NwM7jTFNsS01YhNp01rgmDGmxXruKeAq4LHYl31BE/0Z+j5Wt52IfNn6fpcWpw5+D3CU8G/Mcx38C857TQGQZt3/C+C/rfs+rJNIQCawDXhvvE5ORLs9573mehLjJOZE3p9sIHfE/d8BNzm1PdbjbcBc6/7nga85uT3Wtp8AH7L7sxalz9zlwD7Cv2yF8Enn+53aHutxqXU7FTgIFEa03zg28D3AIcJnav/B2vZXwF9Z968E6q3inzrXAGAx4T813iB81P05uz98E2nPed/jehIgwCf4/sywPqx7rB+qf7C7LRN9f4ClQK31mftFpD9MCdyeLKANyLe7HVFs0xes7XuBHwHpDm/PNmC/9XO0OtJ96qX0SinlUHolplJKOZQGuFJKOZQGuFJKOZQGuFJKOZQGuFJKOZQGuFJKOZQGuFJKOVS8LqVXKiGJyE8IX81XBUwC/toY82tbi1IqQnoErlLdEuCoMeZy4APAP9pcj1IR0ysxVcoSkUzCU5JWGmMGRaQI2G6MmW1zaUpFRI/AVSpbCNQbYwatx8sJz0WhlCNoH7hKZUuAqdZ6im7CEyR91t6SlIqcBrhKZUuAx4GXgDzgy8aYV2ytSKkx0D5wlbJEZCvwF8aYOrtrUWo8NMBVyhKRM4RPYCbK2qRKjYkGuFJKOZSOQlFKKYfSAFdKKYfSAFdKKYfSAFdKKYfSAFdKKYfSAFdKKYfSAFdKKYf6//BdLzi5E4DgAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "extra['updater'].plot_posterior_marginal(idx_param=0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Because the particle filter approximation captures our entire state of knowledge about $p$, we can even use that knowledge to make adaptive decisions about what sequence lengths to measure. As shown in [arXiv:1404.5275](https://arxiv.org/abs/1404.5275), the optimal sequence length to measure is approximately $m = 1 / (1 - F)$. To use that with a statistical inference engine written in Python, though, we'd need to make a round trip back and forth between our quantum device and our Python host with every single measurement — this is clearly not practical in general. Fortunately, the particle filtering algorithm is simple enough we can actually write it directly in Q#!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Adaptive sequence lengths in Q\\#"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to perform the statistical inference we need to choose sequence lenghts adaptively, we need to define our RB model in terms of Q# functions:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%qsharp\n",
    "\n",
    "function RBLikelihood(datum : Result, modelparams : Double[], seqLength : Int) : Double {\n",
    "    let p = modelparams[0];\n",
    "    let A = modelparams[1];\n",
    "    let B = modelparams[2];\n",
    "    let pr0 = A * p^IntAsDouble(seqLength) + B;\n",
    "    return datum == Zero ? pr0 | 1.0 - pr0;\n",
    "}\n",
    "\n",
    "function IsRbModelValid(modelparams : Double[]) : Bool {\n",
    "    let p = modelparams[0];\n",
    "    let A = modelparams[1];\n",
    "    let B = modelparams[2];\n",
    "    \n",
    "    return 0.0 <= p and\n",
    "            p <= 1.0 and\n",
    "            0.0 <= A and\n",
    "            A <= 1.0 and\n",
    "            0.0 <= B and\n",
    "            B <= 1.0 and\n",
    "            A + B <= 1.0 and\n",
    "            A * p + B <= 1.0;\n",
    "}\n",
    "\n",
    "function RBModel()\n",
    ": (\n",
    "    ((Result, Double[], Int) -> Double),\n",
    "    (Double[] -> Bool)\n",
    ") {\n",
    "    return (RBLikelihood, IsRbModelValid);\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we define our prior distribution over $p$, $A$, and $B$, representing what we know about our system before we start. Here, we'll assume that $p \\ge 90\\%$, as we probably think our device has small enough noise for RB to be interesting in the first place!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%qsharp\n",
    "open Microsoft.Quantum.Math;\n",
    "\n",
    "// We'll need to open the namespace containing functions and operations defined locally in this sample.\n",
    "// For details on how the particle filtering algorithm is implemented in Q#,\n",
    "// see Inference.qs and Math.qs in this folder.\n",
    "open Microsoft.Quantum.Samples;\n",
    "\n",
    "operation DrawFromPrior() : Double[] {\n",
    "    mutable model = [0.0, 0.0, 0.0];\n",
    "    let dist = StandardNormalDistribution();\n",
    "\n",
    "    repeat {\n",
    "        // Assume p is at least 90%.\n",
    "        let p = DrawRandomDouble(0.9, 1.0);\n",
    "        let A = 0.5 + Sqrt(0.1^2.0) * dist::Sample();\n",
    "        let B = 0.5 + Sqrt(0.1^2.0) * dist::Sample();\n",
    "\n",
    "        set model = [p, A, B];\n",
    "    }\n",
    "    until IsRbModelValid(model);\n",
    "\n",
    "    return model;\n",
    "}\n",
    "\n",
    "operation DrawInitialPrior(nParticles : Int) : SmcApproximation {\n",
    "    return SmcApproximation(\n",
    "        DrawMany(DrawFromPrior, nParticles, ()),\n",
    "        ConstantArray(nParticles, 1.0 / IntAsDouble(nParticles))\n",
    "    );\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once we have our RB model defined, we can use it to interactively estimate $p$, $A$, and $B$ and decide sequence lengths  on-the-fly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%qsharp\n",
    "open Microsoft.Quantum.Samples;\n",
    "\n",
    "operation RunAdaptiveRBExperiment(nParticles : Int, nShots : Int) : (Double[], Double[][]) {\n",
    "    let dim = 2.0; // We're  working with a single qubit, so 𝑑 = 2.\n",
    "\n",
    "    Message(\"starting with initial prior...\");\n",
    "    mutable dist = DrawInitialPrior(nParticles);\n",
    "    Message(\"...done.\");\n",
    "\n",
    "    let model = RBModel();\n",
    "    let options = DefaultSmcOptions();\n",
    "\n",
    "    for idxShot in 1..nShots {\n",
    "        Message($\"Starting shot {idxShot}...\");\n",
    "        let currentMean = Mean(dist);\n",
    "        let estF = 1.0 - (1.0 - currentMean[0]) * (dim - 1.0) / dim;\n",
    "        let seqLength = Ceiling(1.0 / (1.0 - estF));\n",
    "        let survival = MeasureSingleQubitRBSequence(seqLength);\n",
    "        Message($\"Survival at length {seqLength}: {survival}\");\n",
    "        set dist = Update(dist, survival, seqLength, model, options);\n",
    "    }\n",
    "\n",
    "    return (Mean(dist), Cov(dist));\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "starting with initial prior...\n",
      "...done.\n",
      "Starting shot 1...\n",
      "Survival at length 41: Zero\n",
      "Starting shot 2...\n",
      "Survival at length 45: Zero\n",
      "Starting shot 3...\n",
      "Survival at length 51: Zero\n",
      "Starting shot 4...\n",
      "Survival at length 60: Zero\n",
      "About to resample...\n",
      "Starting shot 5...\n",
      "Survival at length 70: One\n",
      "Starting shot 6...\n",
      "Survival at length 55: One\n",
      "Starting shot 7...\n",
      "Survival at length 47: Zero\n",
      "Starting shot 8...\n",
      "Survival at length 51: Zero\n",
      "Starting shot 9...\n",
      "Survival at length 57: Zero\n",
      "Starting shot 10...\n",
      "Survival at length 64: Zero\n",
      "Starting shot 11...\n",
      "Survival at length 72: Zero\n",
      "Starting shot 12...\n",
      "Survival at length 81: One\n",
      "Starting shot 13...\n",
      "Survival at length 67: One\n",
      "Starting shot 14...\n",
      "Survival at length 58: One\n",
      "Starting shot 15...\n",
      "Survival at length 51: One\n",
      "Starting shot 16...\n",
      "Survival at length 46: Zero\n",
      "Starting shot 17...\n",
      "Survival at length 50: Zero\n",
      "Starting shot 18...\n",
      "Survival at length 53: One\n",
      "Starting shot 19...\n",
      "Survival at length 48: Zero\n",
      "Starting shot 20...\n",
      "Survival at length 52: One\n",
      "Starting shot 21...\n",
      "Survival at length 47: Zero\n",
      "Starting shot 22...\n",
      "Survival at length 50: Zero\n",
      "Starting shot 23...\n",
      "Survival at length 54: Zero\n",
      "Starting shot 24...\n",
      "Survival at length 57: Zero\n",
      "Starting shot 25...\n",
      "Survival at length 61: Zero\n",
      "Starting shot 26...\n",
      "Survival at length 65: Zero\n",
      "Starting shot 27...\n",
      "Survival at length 69: Zero\n",
      "Starting shot 28...\n",
      "Survival at length 74: One\n",
      "Starting shot 29...\n",
      "Survival at length 68: Zero\n",
      "Starting shot 30...\n",
      "Survival at length 72: One\n",
      "Starting shot 31...\n",
      "Survival at length 66: Zero\n",
      "Starting shot 32...\n",
      "Survival at length 70: One\n",
      "Starting shot 33...\n",
      "Survival at length 64: Zero\n",
      "Starting shot 34...\n",
      "Survival at length 68: Zero\n",
      "Starting shot 35...\n",
      "Survival at length 72: Zero\n",
      "Starting shot 36...\n",
      "Survival at length 76: One\n",
      "Starting shot 37...\n",
      "Survival at length 70: One\n",
      "Starting shot 38...\n",
      "Survival at length 65: Zero\n",
      "Starting shot 39...\n",
      "Survival at length 69: Zero\n",
      "Starting shot 40...\n",
      "Survival at length 72: Zero\n",
      "Starting shot 41...\n",
      "Survival at length 76: One\n",
      "Starting shot 42...\n",
      "Survival at length 71: One\n",
      "Starting shot 43...\n",
      "Survival at length 66: Zero\n",
      "Starting shot 44...\n",
      "Survival at length 69: Zero\n",
      "Starting shot 45...\n",
      "Survival at length 73: Zero\n",
      "Starting shot 46...\n",
      "Survival at length 76: One\n",
      "Starting shot 47...\n",
      "Survival at length 71: One\n",
      "Starting shot 48...\n",
      "Survival at length 67: One\n",
      "Starting shot 49...\n",
      "Survival at length 64: Zero\n",
      "Starting shot 50...\n",
      "Survival at length 66: Zero\n",
      "Starting shot 51...\n",
      "Survival at length 69: Zero\n",
      "Starting shot 52...\n",
      "Survival at length 72: One\n",
      "Starting shot 53...\n",
      "Survival at length 68: Zero\n",
      "About to resample...\n",
      "Starting shot 54...\n",
      "Survival at length 70: Zero\n",
      "Starting shot 55...\n",
      "Survival at length 73: Zero\n",
      "Starting shot 56...\n",
      "Survival at length 75: Zero\n",
      "Starting shot 57...\n",
      "Survival at length 78: One\n",
      "Starting shot 58...\n",
      "Survival at length 74: Zero\n",
      "Starting shot 59...\n",
      "Survival at length 77: Zero\n",
      "Starting shot 60...\n",
      "Survival at length 79: One\n",
      "Starting shot 61...\n",
      "Survival at length 76: One\n",
      "Starting shot 62...\n",
      "Survival at length 73: One\n",
      "Starting shot 63...\n",
      "Survival at length 70: Zero\n",
      "Starting shot 64...\n",
      "Survival at length 72: One\n",
      "Starting shot 65...\n",
      "Survival at length 69: One\n",
      "Starting shot 66...\n",
      "Survival at length 66: One\n",
      "Starting shot 67...\n",
      "Survival at length 63: Zero\n",
      "Starting shot 68...\n",
      "Survival at length 65: Zero\n",
      "Starting shot 69...\n",
      "Survival at length 67: One\n",
      "Starting shot 70...\n",
      "Survival at length 65: One\n",
      "Starting shot 71...\n",
      "Survival at length 62: One\n",
      "Starting shot 72...\n",
      "Survival at length 59: Zero\n",
      "Starting shot 73...\n",
      "Survival at length 61: Zero\n",
      "Starting shot 74...\n",
      "Survival at length 63: One\n",
      "Starting shot 75...\n",
      "Survival at length 61: Zero\n",
      "Starting shot 76...\n",
      "Survival at length 63: Zero\n",
      "Starting shot 77...\n",
      "Survival at length 65: One\n",
      "Starting shot 78...\n",
      "Survival at length 62: Zero\n",
      "Starting shot 79...\n",
      "Survival at length 64: Zero\n",
      "Starting shot 80...\n",
      "Survival at length 66: One\n",
      "Starting shot 81...\n",
      "Survival at length 63: One\n",
      "Starting shot 82...\n",
      "Survival at length 61: One\n",
      "Starting shot 83...\n",
      "Survival at length 59: Zero\n",
      "Starting shot 84...\n",
      "Survival at length 60: One\n",
      "Starting shot 85...\n",
      "Survival at length 58: One\n",
      "Starting shot 86...\n",
      "Survival at length 56: One\n",
      "Starting shot 87...\n",
      "Survival at length 53: Zero\n",
      "Starting shot 88...\n",
      "Survival at length 55: One\n",
      "Starting shot 89...\n",
      "Survival at length 53: Zero\n",
      "Starting shot 90...\n",
      "Survival at length 55: Zero\n",
      "Starting shot 91...\n",
      "Survival at length 56: One\n",
      "Starting shot 92...\n",
      "Survival at length 54: One\n",
      "Starting shot 93...\n",
      "Survival at length 52: Zero\n",
      "Starting shot 94...\n",
      "Survival at length 54: Zero\n",
      "Starting shot 95...\n",
      "Survival at length 55: Zero\n",
      "Starting shot 96...\n",
      "Survival at length 57: One\n",
      "Starting shot 97...\n",
      "Survival at length 55: One\n",
      "Starting shot 98...\n",
      "Survival at length 53: One\n",
      "About to resample...\n",
      "Starting shot 99...\n",
      "Survival at length 51: One\n",
      "Starting shot 100...\n",
      "Survival at length 49: Zero\n"
     ]
    }
   ],
   "source": [
    "est_online, cov_online = RunAdaptiveRBExperiment.simulate_noise(nParticles=6000, nShots=100)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.linalg import sqrtm\n",
    "est_err = sqrtm(cov_online)[0, 0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "𝑝 = 0.9601871511798478 ± 0.02050192476404032 (online, 100 shots)\n",
      "𝑝 = 0.9654739702794455 ± 0.005356729111985975 (offline, 2100 shots)\n"
     ]
    }
   ],
   "source": [
    "print(f\"𝑝 = {est_online[0]} ± {est_err} (online, 100 shots)\")\n",
    "print(f\"𝑝 = {est[0]} ± {sqrtm(cov)[0, 0]} (offline, {100 * len(seq_lens)} shots)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "By implementing our sequence generation and statistical inference in Q#, we were able to learn our RB parameters using only 100 shots nearly as well as we did with 2100 shots when we chose sequence lengths non-adaptively."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Epilogue"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for component, version in sorted(qsharp.component_versions().items(), key=lambda x: x[0]):\n",
    "    print(f\"{component:20}{version}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "print(sys.version)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}