jupyter_notebooks/Tutorials/18 Randomized Benchmarking Data Analysis.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Randomized Benchmarking Data Analysis\n",
    "This tutorial demonstrates how to analyze RB data using `pyGSTi`. This data analysis is for any RB method whereby *average success/survival probabilities* (ASPs) as a function of *RB sequence length* ($m$) are fit to the function\n",
    "$$ P_m = A + B p^m,$$\n",
    "where $A$, $B$ and $p$ are fit parameters (and $A$ may be fixed to a constant). So this includes [Direct RB](https://arxiv.org/abs/1807.07975) and [Clifford RB](http://journals.aps.org/prl/abstract/10.1103/PhysRevLett.106.180504), as well as [Unitarity RB](http://iopscience.iop.org/article/10.1088/1367-2630/17/11/113020/meta) (although note that in that case the desired quantity is $p$ rather than $r \\approx 1-p$). This analysis is not suitable for any types of RB that require more complicated analysis, e.g., [Dihedral RB](https://journals.aps.org/pra/abstract/10.1103/PhysRevA.92.060302), although note that with minor post-processing these tools can be used to analyzes [interleaved RB](https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.109.080505) data (full support for analyzing interleaved RB data may be added at a later date)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from __future__ import print_function #python 2 & 3 compatibility\n",
    "import pygsti\n",
    "from pygsti.extras import rb"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generating some fake data to analyze\n",
    "To show how to do use the RB analysis, we'll generate some fake [\"Direct randomized benchmarking\"](https://arxiv.org/abs/1807.07975) (DRB) data. But note that the analysis is the same for DRB and Clifford RB: fitting observed average success probability (ASP) data to the exponential decay $P_m = A + Bp^m$, so it is not of much importance for this tutorial that the data is from DRB.\n",
    "\n",
    "We'll create 5-qubit DRB data, for an (imaginary) device containing 5 qubits with ring connectivity. We'll assume that the errors consist of local depolarizing after every gate with a gate- and qubit-independent depolarizing error rate of 0.1%. I.e., each circuit layer is followed by the channel $\\mathcal{D}^{\\otimes 5}_{\\lambda}$ where\n",
    "$$ \\mathcal{D}_{\\lambda}[\\rho] = \\lambda \\rho + (1-\\lambda) \\mathbb{I} ,$$\n",
    "is a 1-qubit uniformly depolarizing channel. To obtain a per-qubit error rate of 0.1% we would set\n",
    "$\\lambda = 3 (1-0.001) / 4$.\n",
    "\n",
    "Because these simulations are on 5 qubits, they take a while to run. So we have already saved data from these simulations. To re-run the simulations set `runsims = True`, which will overwrite the old data. Otherwise, the data previously generated by this code will be imported."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "runsims = False"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, we create a `ProcessorSpec` object encoding the device to simulate: see previous tutorials for more information."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "if runsims:\n",
    "    \n",
    "    nQubits = 5 \n",
    "    qubit_labels = ['Q0','Q1','Q2','Q3','Q4'] \n",
    "    gate_names = ['Gi', 'Gxpi2', 'Gxmpi2', 'Gypi2', 'Gympi2', 'Gcphase'] \n",
    "    availability = {'Gcphase':[('Q0','Q1'), ('Q1','Q2'), ('Q2','Q3'), \n",
    "                               ('Q3','Q4'), ('Q4','Q0')]}\n",
    "    pspec = pygsti.obj.ProcessorSpec(nQubits, gate_names, availability=availability, \n",
    "                                     qubit_labels=qubit_labels)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we create the error model, and specify the RB parameters (RB lengths etc)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "if runsims: \n",
    "    \n",
    "    # The local error rate of every qubit\n",
    "    errorrate = 0.001\n",
    "    \n",
    "    # Put these into a dictionary, mapping the qubit label to the error rate (which is qubit-indep here)\n",
    "    gate_errorrate_dict = {q : errorrate for q in qubit_labels}\n",
    "    \n",
    "    # The error type : 'uniform' means locally uniform depolarization.\n",
    "    ptype = 'uniform'\n",
    "    \n",
    "    # This creates this error model in the format needed for the simulator\n",
    "    errormodel = rb.simulate.create_locally_gate_independent_pauli_error_model(pspec, gate_errorrate_dict,\n",
    "                                                                               ptype='uniform')\n",
    "    \n",
    "    # The DRB lengths\n",
    "    lengths = [0,10,20,30,50,100,200,400]\n",
    "\n",
    "    # The number of circuits per length\n",
    "    k = 30\n",
    "    \n",
    "    # The counts for each circuit. Below we use a basic stochastic-unravelling simulator,\n",
    "    # so the time taken is linear in the number of counts.\n",
    "    counts = 50"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we run the simulations using `rb.simulate.rb_with_pauli_errors()`. This samples DRB circuits as specified (here using the default sampler as we leave it unspecified) and writes the data to file. This simulator is of only tangential relevance to this tutorial, so is not explained any further here. See the docstrings for more information."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "if runsims:\n",
    "    \n",
    "    filename = 'tutorial_files/MySimulatedDRBData.txt'\n",
    "    rbdata = rb.simulate.rb_with_pauli_errors(pspec, errormodel, lengths, k, counts, \n",
    "                                              rbtype='DRB', filename=filename, verbosity=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Import RB summary data\n",
    "Currently, `pyGSTi` only has functions for analyzing RB data that has been summarized into an \"RB length\" vs. \"success counts\" format. This format is demonstrated in the file `tutorial_files/MySimulatedDRBData.txt` (note that the final two columns are optional and contain auxillary information about the RB circuits). So, if the data is stored in a `DataSet` object, this has to be pre-processed manually into the above format. Functions to implement this conversion and/or analysis functions that can handle `DataSet` objects will be added in a future `pyGSTi` release."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Importing tutorial_files/MySimulatedDRBData.txt...Complete.\n"
     ]
    }
   ],
   "source": [
    "rbdata = rb.io.import_rb_summary_data('tutorial_files/MySimulatedDRBData.txt')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Implementing the analysis\n",
    "Once we have an `RBSummaryDataset` object, implement the analysis is simple. To run a \"standard practice\" analysis, this data just needs to be passed to `rb.analysis.std_practice_analysis()`, as we do below. This function implements an ordinary unweighted least-squares analysis, and estimates \"error bars\" on the extracted parameters using a standard non-parameteric boostrap. Note that this is *not* necessarily the best way to analyze RB data: there has recently been some developments in RB statistics (in particular, see https://arxiv.org/abs/1802.00401), some of which may be significant improvements on this method. But, we expect that this analysis method performs reasonably well in most circumstances.\n",
    "\n",
    "This function also has some useful optional arguments, which it is sometimes important to adjust:\n",
    "\n",
    "- `seed` is a list [$B$,$p$] that provides an initial guess, or *seed*, for the fit of $B$ and $p$ in $P_m = A + Bp^m$. It is important to adjust these when the default values of [0.8,0.95] are far from the optimum.\n",
    "- `bootstrap_samples` is the number of samples used for the bootstrapped standard deviation of $r$. Increase this value, from the default value of 200, for a slower analysis that reports a more accurate standard deviation.\n",
    "- The analysis performs two fits: the first one to $P_m = A + Bp^m$ with $A$ fixed and $B$ and $p$ seeded as above; the second one with all of $A$, $B$ and $p$ as variables and $B$ and $p$ seeded from the fixed-$A$ fit. The optional parameter `asymptote` sets the value of $A$ for the fixed-asymptote fit, and it is also the seed for $A$ in the full fit. If not specified, it defaults to $1/2^n$ for $n$ qubits, which is the correct value for $A$ whenever the measurement errors are not biased and the gate errors are not strongly non-unital (if their is a small non-unital error in the gates, e.g., weak amplitude damping, the true $A$ is a small perturbation on $1/2^n$ so this value is probably fine).\n",
    "- `rtype` is the \"type\" of RB number. The default is `rtype = 'EI'`, giving an $n$-qubit RB $r$ defined by:\n",
    "$$ r = \\frac{(4^n - 1)(1 - p)}{4^n}. $$\n",
    "This corresponds to an $r$ that is, for Clifford or direct RB with stochastic errors, approximately the probability of error in a gate. More generally, for Clifford and direct RB this $r$ is directly related to the entanglement infidelity averaged over all gates [(modulo some \"gauge\" problems)](https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.119.130502).\n",
    "This is our preferred definition of the RB error rate, as a function of $p$, for both Clifford and direct RB (and other similar forms of RB). But note that it is not the conventional definition in the Clifford RB literature. The convention in Clifford RB corresponds to setting `rtype = 'AGI'`, giving an $n$-qubit RB $r$ defined by:\n",
    "$$ r = \\frac{(2^n - 1)(1 - p)}{2^n}. $$\n",
    "In this case, the Clifford and direct RB $r$ is directly related to the average gate infidelity (AGI) averaged over all gates [(again, modulo to the same \"gauge\" problems)](https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.119.130502)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "rbresults = rb.analysis.std_practice_analysis(rbdata)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And that's it.\n",
    "\n",
    "## Looking at the results\n",
    "\n",
    "The analysis returns an `RBResults` object, that encapsulates the analysis results and the original RB summary data.\n",
    "\n",
    "### Plotting the data and the fit\n",
    "\n",
    "Once we've done the analysis, we can plot the data and our fit to the data, using the `.plot()` method. By default, this plots the \"full\" fit to $P_m = A + B p^m$ with all of $A$, $B$ and $p$ as variables. (If you would like to plot a different fit then set the optional argument `fitkey` to the key in `rbresults.fits` dictionary that this fit corresponds to: see below for more info on `rbresults.fits`)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfUAAAFECAYAAAA+3qfBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAAIABJREFUeJzs3Xd8VFX+//HXJyEhoUjoSoISpEgPXQUVUaoKiIq49roWVld3+apbLOvuqj9dC6hrQ3BFAQuCBcGCDTsdpFcliGKUJi3l/P64kxBCyk2YyUxm3k8e95GZO2fufO7c6Cen3HPMOYeIiIhUfXHhDkBERESCQ0ldREQkSiipi4iIRAkldRERkSihpC4iIhIllNRFRESihJK6iIhIlFBSFxERiRJK6iIiIlGiWrgDKK8GDRq4Zs2ahTsMERGRSjFv3ryfnXMN/ZStckm9WbNmzJ07N9xhiIiIVAoz2+i3rJrfRUREooSSuoiISJRQUhcREYkSVa5PXURiW3Z2Nps2bWLv3r3hDkUkqJKSkkhLSyMhIaHCx1BSF5EqZdOmTdSuXZtmzZphZuEORyQonHNkZWWxadMm0tPTK3wcNb+LSJWyd+9e6tevr4QuUcXMqF+//mG3QCmpi0iVo4Qu0SgYv9dK6iIiIlFCSV1EpAKmTZuGmbFixYpwhyKFXHbZZbz66qu+y3/00UeceeaZxb42ePBgtm3bBkCtWrUA2Lx5M+eeey4ACxcuZMaMGQXl33jjDe67776Khh4USuoiEtWmLcik132zSb/tbXrdN5tpCzKDctxJkybRu3dvJk2aFJTj5eTkBOU4saCyvqsZM2aQkpJy0L4mTZoU/NFQNKkPGTKE2267rVJiK0nMJ/W9ObotRiRaTVuQye1Tl5C5bQ8OyNy2h9unLjnsxL5r1y7mzJnDuHHjmDx5csH+kSNH8vbbbxc8z6815ubmMnr0aLp3707Hjh156qmnAK+WeNJJJzFkyBDatm0LwLBhw+jatSvt2rXj6aefLjjWuHHjaNWqFT169ODqq69m1KhRAGzdupVzzjmH7t270717dz777LND4v3222/p0aMHGRkZdOzYkdWrV7Nhwwbat29fUObBBx/krrvuAmDNmjWcfvrpdOrUiS5durB27VoA7r//fjp06ECnTp0KktfatWsZOHAgXbt25aSTTipouXjllVdo3749nTp14uSTTy4xjqJq1arFzTffTLt27TjttNPYunUrAH369OGPf/wj3bp149FHH2XDhg307duXjh07ctppp/Hdd98VHOP999+nW7dutGrVirfeeguADRs2cNJJJ9GlSxe6dOnC559/XlB+x44dnHHGGbRu3Zprr72WvLw8wJuW/Oeffz4ovvzvbf/+/dxxxx1MmTKFjIwMpkyZwoQJE8q8Lh9//DEZGRlkZGTQuXNndu7cech3cFicc1Vq69q1qwum7fu2B/V4IhJay5Yt8132xHs/cMfc+tYh24n3fnBYMUycONFdccUVzjnnTjjhBDd37lznnHNTp051l1xyiXPOuX379rm0tDS3e/du99RTT7l77rnHOefc3r17XdeuXd26devchx9+6GrUqOHWrVtXcOysrCznnHO7d+927dq1cz///LPLzMx0xxxzjMvKynL79+93vXv3djfccINzzrkLLrjAffrpp8455zZu3OiOO+64Q+IdNWqUmzhxYkFcu3fvduvXr3ft2rUrKPPAAw+4O++80znnXI8ePdzUqVOdc87t2bPH/fbbb27GjBnuhBNOcL/99ttBcfbt29etWrXKOefcl19+6U499VTnnHPt27d3mzZtcs459+uvv5YYR1FAQZm777674DxPOeUUd9111xWUO/PMM92ECROcc86NGzfODR061Dnn3KWXXuoGDBjgcnNz3apVq1xqamrBOezZs8c559yqVatcfi758MMPXfXq1d3atWtdTk6OO/30090rr7zinHPumGOOcVu3bnXOOVezZk3nnDvoexs/fnxBfEWfl3RdzjzzTDdnzhznnHM7d+502dnZB51/cb/fwFznM0fG/H3q3vclItFo87Y95drv16RJk7jpppsAr3Y+adIkunbtyqBBg7jpppvYt28fM2fO5OSTTyY5OZl3332XxYsXFzTbbt++ndWrV5OYmEiPHj0Oui95zJgxvP766wB8//33rF69mi1btnDKKadQr149AM477zxWrVoFeLXSZcuWFbx/x44d7Nq1q6APGOCEE07gX//6F5s2bWL48OG0bNmyxHPbuXMnmZmZnH322YA3IUr+51x++eXUqFEDgHr16rFr1y4+//xzzjvvvIL379u3D4BevXpx2WWXMWLECIYPH+47jri4OM4//3wALrroooL3AgX7Ab744gumTp0KwMUXX8z//d//Fbw2YsQI4uLiaNmyJc2bN2fFihWkp6czatQoFi5cSHx8fMH3B9CjRw+aN28OwAUXXMCcOXMK+s0rqqTr0qtXL2655RYuvPBChg8fTlpa2mF9TlFK6krqIlGrSUoymcUk8CYpyRU+5i+//MLs2bNZsmQJZkZubi5mxgMPPEBSUhJ9+vRh1qxZTJkyhZEjRwLe/2fGjh3LgAEDDjrWRx99RM2aNQ96/v777/PFF19Qo0YN+vTpU+Z9y3l5eXz55ZcFybc4v/vd7+jZsydvv/02gwcP5qmnnqJVq1YFzcxAhe6PzsvLIyUlhYULFx7y2pNPPslXX33F22+/TdeuXZk3b16xcfTt27fUzyh8m1fh78rve/KfP/zwwzRu3JhFixaRl5d30PdVXPnDVdJ1ue222zjjjDOYMWMGvXr1YtasWRx33HGH/Xn5Yr5PPY88JXaRKDV6QGuSE+IP2pecEM/oAa0rfMxXX32Viy++mI0bN7Jhwwa+//570tPT+fTTTwGvNjl+/Hg+/fRTBg4cCMCAAQP473//S3Z2NgCrVq3it99+O+TY27dvp27dutSoUYMVK1bw5ZdfAtC9e3c+/vhjfv31V3JycnjttdcK3tO/f3/Gjh1b8Ly4BLtu3TqaN2/OjTfeyNChQ1m8eDGNGzfmp59+Iisri3379hX0PdeuXZu0tDSmTZsGeDXv3bt3069fP8aPH8/u3bsB74+bI444gvT0dF555RXA++Nl0aJFgNfX3rNnT/7xj3/QsGFDvv/++2LjKCovL6+gReOll16id+/exV6HE088sWA8w4svvshJJ51U8Norr7xCXl4ea9euZd26dbRu3Zrt27dz1FFHERcXxwsvvEBubm5B+a+//pr169eTl5fHlClTSvzMomrXrl1in3hJ12Xt2rV06NCBW2+9le7duwf97omYT+oiEr2GdU7l3uEdSE1JxoDUlGTuHd6BYZ1TK3zMSZMmFTRN5zvnnHMKRsH379+fjz/+mNNPP53ExEQArrrqKtq2bUuXLl1o3749v//974sdwT1w4EBycnJo06YNt912G8cffzwAqamp/OUvf6FHjx706tWLZs2aUadOHcBrrp87dy4dO3akbdu2PPnkk4cc9+WXX6Z9+/ZkZGSwdOlSLrnkEhISErjjjjvo0aMH/fr1O6i2+MILLzBmzBg6duzIiSeeyJYtWxg4cCBDhgyhW7duZGRk8OCDDwJeQh03bhydOnWiXbt2TJ8+HYDRo0fToUMH2rdvz4knnkinTp2KjaOomjVr8vXXX9O+fXtmz57NHXfcUex1GDt2LOPHj6djx4688MILPProowWvHX300fTo0YNBgwbx5JNPkpSUxPXXX8/zzz9Pp06dWLFixUG1/u7duzNq1CjatGlDenr6Ide3JKeeeirLli0rGChXWEnX5ZFHHqF9+/Z07NiRhIQEBg0a5Ouz/LKqVkvt1q2bmzt3btCOl7Uni7pJdYkz/X0jUhUsX76cNm3ahDuMSpffT56Tk8PZZ5/NFVdc4Tv5VCW1atVi165d4Q4jbIr7/Tazec65bn7eH/OZzOHU/C4iEe+uu+4iIyOD9u3bk56ezrBhw8IdkkSgmB8oB15iFxGJZPnN3dEulmvpwaCaunNK6iIiEhViPqnnOY1+FxGR6BDzSd3hr6aenZdNTp6/+Yb35+4/3LBERETKLeaTOvibgCY7N9t3st6bq/nkRUSk8sV8Us9zeeS5vDLLlaffXc35ItEtPj6+YFGOjIwMNmzYwNy5c7nxxhsBb2a4wguGVIZ9+/Zx/vnn06JFC3r27MmGDRuKLTdz5kxat25NixYtDlomdP369fTs2ZMWLVpw/vnns3+/V4mZMGECDRs2LDjXZ5999qDj7dixg7S0tIKFTCoiOzubLl26lFnu3nvvpUWLFrRu3ZpZs2YVW6ak8yjt+ynpuNu2bePcc8/luOOOo02bNnzxxRcVPsfKEvNJ3W+y1oA6EcmXnJzMwoULC7ZmzZrRrVs3xowZA1QsqR/ucqLjxo2jbt26rFmzhptvvplbb731kDK5ubnccMMNvPPOOyxbtoxJkyYVzE9+6623cvPNN7NmzRrq1q3LuHHjCt53/vnnF5zrVVddddAx//73vxeswlaWZs2aFbt/zpw59OrVq9T3Llu2jMmTJ/Ptt98yc+ZMrr/++oNmhctX0nmU9P2UdtybbrqJgQMHsmLFChYtWlQl5keI+aRenpq6nxq4aukisemjjz7izDPPZMOGDTz55JM8/PDDZGRkFEwfW5wJEyYwZMgQ+vbty2mnnXZYnz99+nQuvfRSAM4991w++OCDQ/5/9PXXX9OiRQuaN29OYmIiI0eOZPr06TjnmD17dsEiJpdeemnBNLGlmTdvHj/++CP9+/c/rNhnzpxZ5sxq06dPZ+TIkVSvXp309HRatGjB119/fVCZ0s6jpO+npONu376dTz75hCuvvBKAxMTEQ9ZWj0S6T935q607HMbhT/IvIkH0zm2wZUlwj3lkBxh0X6lF9uzZQ0ZGBgDp6ekFq6qBVxu99tprqVWrFn/+85/L/Lj58+ezePHighXYCjvppJOKnVv8wQcf5PTTTz9oX2ZmJk2bNgWgWrVq1KlTh6ysLBo0aFBsGYC0tDS++uorsrKySElJoVq1agX7MzMPrDn/2muv8cknn9CqVSsefvhhmjZtSl5eHn/605+YOHEi77//fpnnWZoPP/yQO++8s9QymZmZBdPmFhcjUOp5lPT9lHTc5ORkGjZsyOWXX86iRYvo2rUrjz76qO9FZcIl5pO639HvzjmcBbcWvi93H9Xjqwf1mCISevnN78HQr1+/YhM6UGotv7KcddZZXHDBBVSvXp2nnnqKSy+9lNmzZ/PEE08wePDgMpcOveGGG/jss88A2Lx5c8EfQ+eddx5//etfyczMpF69egVLukaKnJwc5s+fz9ixY+nZsyc33XQT9913H/fcc0+4QyuVkno5pon1W6P3ezw/zf4iUooyatRVQWk1v/LU1FNTU/n+++9JS0sjJyeH7du3U79+/WLL5Nu0aROpqanUr1+fbdu2kZOTQ7Vq1Qr2Awcd46qrripYt/yLL77g008/5YknnmDXrl3s37+fWrVqHTT4DuDxxx8veNysWbND/hiaOXNmwZK0jz/+OM888wwAM2bMoEmTJmXGXlhp51HS91PScdPS0khLS6Nnz56A12Rf9NwikfrUfU4+41zw54hX/7tIdCq6JOfrr7/O7bffXu7jfPrppwcNyMvfiiZ0gCFDhvD8888D3vKwffv2PWRd8O7du7N69WrWr1/P/v37mTx5MkOGDMHMOPXUUwuWPH3++ecZOnQoAD/88EPB+994442CwWIvvvgi3333HRs2bODBBx/kkksuqVDSK9yffsMNNxScY+GEnn9+kydPZt++faxfv57Vq1fTo0ePg8qUdh4lfT8lHffII4+kadOmrFy5EoAPPviAtm3blvv8KlvMJ3Xfze8a+S4iPp111lm8/vrrBQPl1q5dyxFHHBHSz7zyyivJysqiRYsWPPTQQwUJdvPmzQwePBjw+pIfe+wxBgwYQJs2bRgxYgTt2rUD4P777+ehhx6iRYsWZGVlFQwQGzNmDO3ataNTp06MGTOGCRMmBC3m3Nxc1qxZc9CyryVp164dI0aMoG3btgwcOJDHH3+c+Ph4AAYPHszmzZtLPY+Svp/Sjjt27FguvPBCOnbsyMKFC/nLX/4StHMPlZhfenXttrWkVE+hfnL9Usv9svcX4ogjJan00Y+5ebns2L+Dukl1y/zs3dm7qZEQWf1IIpGuKi69etFFF/Hwww/TsGHDcIcSUebMmcPEiROLXQM+Vh3u0qvqUy/H/efBLicisWHixInhDiEi9e7dm969e4c7jKiipF6OyWf83NFWrpnnlPxFRCSIYr5P3S/fk8+UZzR9Fev6EIkU+m9HolEwfq9jPqkH+5Y2EQmtpKQksrKylNglqjjnyMrKIikp6bCOo+Z3n33qzrlDbg8pvqD63kVCKS0tjU2bNrF169ZwhyISVElJSWVO5lOWmE/q4K/Jw+80sUrUIqGVkJBAenp6uMMQiUhK6j75rtGrOV9ERMIkppP6tAWZ/PudDWzdkUOTlGRGD2jNsM6pxZZ1OH+ruZXnFjn1CYqISBDFbFKftiCT26cuYU+2t25u5rY93D7VW+2ppMTuh99pZ0FJXUREgitmR78/MGtlQULPtyc7lwdmrSy2vCapERGRSBezSX3ztj0AtLLvqceOQ/YXVZ5ErVXaREQkHGI2qTdJSSaJfUxMvJdxiQ+SxL6C/SXxu5pbHv6StWrqIiISTDGb1EcPaI0l1ODv2ZfTydYyJuExaiYYowe0Lra839Xc8vDfpy4iIhJMMTtQLn8w3L/fMf7x2y/clfA8b6a/SfOMQcWW9z33u0a/i4hImIQ0qZvZQOBRIB541jl3X5HXjwaeB1ICZW5zzs0IZUyFDeucSuujd1A36S8wpzbNv3gMPm8DvW48pKzzporzRX3qIiISDiFrfjezeOBxYBDQFrjAzNoWKfY34GXnXGdgJPBEqOIpU797oN3Z8N7fYelrh77uc/rXPJen0e8iIhIWoexT7wGscc6tc87tByYDQ4uUccARgcd1gM0hjKd0cXEw7Ek4+kR4/VrY8NlBL4di9Lua30VEJJhCmdRTge8LPd8U2FfYXcBFZrYJmAH8IYTxFOugRVoSkmDki1C3GUy+ALYeuGfdb/N7eUa/+y0nIiLiR7hHv18ATHDOpQGDgRfM7JCYzOwaM5trZnODvTKTBf4VqFEPLnwV4qvDxHNh55aCl/zO/e63VT0vT0ldRESCJ5RJPRNoWuh5WmBfYVcCLwM4574AkoAGRQ/knHvaOdfNOdetYcOGQQ/0kCVV6x4DF74Mu7PgpRGwb1e5lmj1W041dRERCaZQJvVvgJZmlm5miXgD4d4oUuY74DQAM2uDl9QrdZFkMyt+SdUmneG8CbBlKbx8CeRm+zqe3z718vS9i4iI+BGypO6cywFGAbOA5Xij3L81s3+Y2ZBAsT8BV5vZImAScJkLQ6Y7pKaer1V/OOsRWPsBKx+7iVP/vZhe981m2oKiDQ4H+K2pl2eUvIiIiB8hvU89cM/5jCL77ij0eBnQK5QxlKXYWnoh0+w01uZdwJ+YxB3VanLXtktLXc3N78xzfsuJiIj4FbMzyuUrK6k/MGslmfvPpHa1HVxT7W1+cUcwJns4D8xaWfLa634aG5xuaRMRkeCK+aQeZ3HEHTrgvoC3apvx75zfUc92ckvCq/xCbV7c1u+wPldzxIuISLDFfFI3M+ItvsTXm6QkkxlI7LdmX00ddvGPahOwGvWBMyr8ueWZI15ERMSPcN+nHnal1dLBW80tOcFL+rnEMyr7RubTmrtzHoE1HxxS3vdscjjN/S4iIkGlpE5cqTX1YZ1TuXd4BxodUQ0DGqTU4cczJhDXqA1MuRg2zS0oO21BJmePXcLA/7fG1yh5ERGRYFLzu1nJt7QFDOucSsfm+4gjjmZ1mnk727wGz/WHF8+Fy2cyLbM2t09dwp7sXAAyt+0pdZQ8BH9BlzyXV2bLg4iIRK+YzwBxFlfmCPh8ByXh2o3h4tchLgFeOJuJ73xSkNDz7cnO5YFZKylOKCaf2b5ve1CPJyIiVYuSehmj3/MZdmi5es29xJ69m4f3/p0jyTrkfd7o+UOFYqBcrsstu5CIiEStmE/qxSbrEsoVW6M/sj1c/Dp1437jxcR/04CDa8tNUpKLPV4oJp9RP72ISGxTUvfRp15QrqRm+tQuzO/9NEfZL0xM/Dcp7AQgOSGe0QNaF/uW8qzm5pdq6iIisS3mk7rfPvU44iit2Mmnn8X8Xk+QHreF/yXeR6s6edw7vEOps86ppi4iIsEU86Pfoex71QGwsqeU7d3/XH45MocO027g3UZjoe3UUssHOwmrpi4iEttUU/dZU/fb977/2D5sHvhP2PQNTBoJ2ZU3UE41dRGR2KakTtzh96kXsatFXxj2JGyY401Qk7PvkDKhGCinmrqISGxTUo/z9xX4Tf4FOp3vrcW+5j149QrIzT7o5VDcp65pZ0VEYpuSus+voFwJPV/Xy2Dg/bDirWITe7ApqYuIxDYl9XJMq+q37/2gPwCOvxYG3AvL3wh5Ylfzu4hIbIv5pO6X34Fy+WUPcsL1hRL75SFL7FrKVUQktimp+3TYk9SccD0MvA+WvwmvXBaSxK6auohIbIv5+9T99pVXuJZe2PHXAQYzb2Xxqp+4ZvcoGqb8wugBrUucpKY8dEubiEhsi/maut/b1MBfYi+rRj+t+ln8K+8y+uR9zZiEsfy0bSe3T11S6trrfqmmLiIS25TUy1FT9/sHQGnlHpi1kmf29+fO7EsZGP8NjyWMISd7X4lLtJaHRr+LiMS2mE/qvm9pK8dAudJyf/5SrM/nDuDO7EsZED+XxxMeZeu2nf6OXQoldRGR2BbzfeoJ8Qm+ygWjlg7eUqyZhRK7A/6R8DzPJzwC2f0hofilWv1QUhcRiW0xX1P3qzyj30ur/Y8e0JrkhPiC5//LHcAdeddwfN4CePE82FfxGruSuohIbFNS96msZF1QrujkM0UM65zKvcM7cFSd6hiQmpJMl7P/iA1/GjZ+Di+cDXt+rVCMGignIhLbYr753S/fC79Q9sIvwzqn0qdNTb7b+R0dG3YM7B3hNb2/cjk8fxZcPA1qNihXjLqlTUQktqmmXg6+pon12Uxf7PHanAUXTIafV8P4wbDjh4qEKSIiMUpJ3SffidrnEq0lJv+Wp8NFr8GOTBg/EH7d6DtGTRMrIhLblNR9Ks/o9wqt6FZYs95wyXSvb338IPh5ja+3qfldRCS2KamXQ7D61POPVWq5tG5w2duQs89L7FuWlno8JXQREVFSLwc/ydrvzHOGlTpJDQBHdoDL34G4al4f+8bPSy2u5ncRkdhWZlI3s7PMyrHoeJQq1wC4INboadgKrnwXajXybndbMaPYYkroIiLiJ1mfD6w2s/9nZseFOqBI5btP3edAufIck5SmcMUsaNQWplwECyb6e5+IiMSUMpO6c+4ioDOwFphgZl+Y2TVmVjvk0UUQ36PafQ6UK0/NH4Ca9eHSNyH9ZJh+A8x5xP97RUQkJvhqVnfO7QBeBSYDRwFnA/PN7A8hjC3iBGUAXKFjlWfZVwCq14LfvQzthsP7d8K7f4M8b2pY55wGy4mIxLgyZ5Qzs6HAZUAL4H9AD+fcT2ZWA1gGjA1phBGiUm9pK021RDhnnDfb3Odj4bcsGDIGzN+MdyIiEr38TBM7HHjYOfdJ4Z3Oud1mdmVowopQPnJmudZnr2gSjouDQf8PajaED/8Fu7Ow856r2LFERCRq+Gl+31I0oZvZ/QDOuQ9CElUV5nfd9fIMqCvhAHDK/8GZD8Pqd7Hnh5KwZ3vFjyciIlWen6Ter5h9g4IdSLSIsziS4pPKLFehPvXidLsCzp8IPy4l483RvmefExGR6FNiUjez68xsCXCcmS0utK0HFldeiJHBbwKOszhqJNTwdbyg9YG3ORN36RvE7/8NxvWD774KznFFRKRKKa2m/hJwFjA98DN/6xq4zU0Oh/nvf/d1uLQeLBzyICTX9ZZu/XZa0I4tIiJVQ2lJ3TnnNgA3ADsLbZhZvdCHFt2C0vQeMG1BJr3vn82I/9Vi8G93kFWnLbxyGXz+GOg2NxGRmFHa6PeXgDOBeYDj4LHfDmgewrgiUjATsWHEBWHq/WkLMrl96hL2ZOcCsGx7AqftvoW3014g9d2/wrbvYOC9EBd/2J8lIiKRrcSk7pw7M/AzvfLCiWBBvgW83DPKleCBWSsLEnq+bdnxnP/L75lzQiv44jHYvgnOeRYSy+7rFxGRqqvEpG5mXUp7o3NufvDDiR3BGv2+edueYvdnbt8HA/4FKUfDO7fChMEwchIccdRhf2ao/Zb9G3tz9lI/uX64QxERqVJKa37/TymvOaBvkGOJKcGqqTdJSSazmMTeJCXZe9Dz915if/VKeKYvXDAJmmQc9ueGUk5eDtl52eEOQ0SkyimxU9c5d2opW0wm9GD2qQfreKMHtCY54eD+8uSEeEYPaH1gR+tB3vKtcfHw3EBYNv2wPzfUNI+9iEj5lXafet/Az+HFbX4ObmYDzWylma0xs9tKKDPCzJaZ2bdm9lLFTqNq8jPzXFmGdU7l3uEdSE1JxoDUlGTuHd6BYZ1TDy54ZHu4ejYc2QFevgQ+eSBiR8YroYuIVExpze+nALPx7k0vygFTSzuwmcUDj+PNSLcJ+MbM3nDOLStUpiVwO9DLOfermTUqZ/xVWrBq/sM6pzKscypfbP6CE5qcUHLBWo285Vvf+APM/idsXQVDxkJC2TPgVSYX+CciIuVT2uj3OwM/L6/gsXsAa5xz6wDMbDIwFG9lt3xXA487534NfNZPFfysKinYq6r5Ol5CEgx/Ghq2htn3wK/rYeRLXsKPEErqIiIVU2b7r5nVN7MxZjbfzOaZ2aNm5mdYcirwfaHnmwL7CmsFtDKzz8zsSzMbWEIM15jZXDObu3XrVh8fXTWErY/eDE7+M4x4AX78Fp4+FbYsCWosh0Nrw4uIVIyfTt3JwFbgHODcwOMpQfr8akBLoA9wAfCMmaUULeSce9o51805161hw4ZB+ujwC3ZNvdx99G2HwOXvgMuDcf2KmOwpAAAgAElEQVRhaak9KpVGNXURkYrxkwWOcs7d45xbH9j+CTT28b5MoGmh52mBfYVtAt5wzmU759YDq/CSfEwIdk29QppkwDUfeQPoXr0c3rsT8nLLelfoKaeLiJSbn6T+rpmNNLO4wDYCmOXjfd8ALc0s3cwSgZHAG0XKTMOrpWNmDfCa49f5jr4ShSIBh72mnq92Y7j0LW8Z188egRfPg92/BDW28lJNXUSk/Eq7pW2nme3AG8z2ErA/sE0GrinrwM65HGAU3h8Ay4GXnXPfmtk/zGxIoNgsIMvMlgEfAqOdc1mHc0JVSUTd914tEc58GM56FNZ/As+c6vW3h4H600VEKqa00e+1D/fgzrkZwIwi++4o9NgBtwS2iBZnccFPwuEY/V6WrpdBo7Yw5WJ4th8MewLaDTv845aD+tRFRCrGV3utmdU1sx5mdnL+FurAIk2cxQU/CUdSTb2wpj28fvbG7eCVS+H9uyu1n101dRGRivFzS9tVwCd4TeV3B37eFdqwpCKC+kfHEUfBZW95Nfc5D8FLIyq1n12JXUSk/PzU1G8CugMbnXOnAp2BbSGNKkZEbE09X7XqXh/7mY94/exPnQyb5gX3M4qh5ncRkYrxk9T3Ouf2AphZdefcCqB1Ge8RHyK9Ob9At8vhilmAwXMD4OtnInbeeBGRWOYnqW8KTAgzDXjPzKYDG0MbllREsP9IOEhqF/j9x3BsX5jxZ3jtSti3K2Qfp5q6iEj5lbagCwDOubMDD+8ysw+BOsDMkEYlFRLSpA5Qox5cMBk+e9hbEGbLEm+q2UbHBfVjnHOafEZEpAL8jn7vYmY3Ah2BTc65/aENSyqiUmaoi4uDk/4El0yHPb9697MvfiXoH6OauohI+fkZ/X4H8DxQH2gAjDezv4U6MCm/Sp12Nv1k+P2ncFQGTL0K3roZsvcG5dBK6CIiFeOnpn4h0N05d2dgOdbjgYtDG1ZsiMjJZ8rjiKPg0jfgxBth7nPw7GneGu0iIhIWfpL6ZiCp0PPqHLowi1RAxN/S5kd8AvS/B373Cuz8AZ4+BRZMPKzR8Q4tvSoiUhGlzf0+1szGANuBb81sgpmNB5ai+9SDolpcmeMUy6XSa+qFteoP186B1K4w/QaYejXs2xm+eEREYlBpWWVu4Oc84PVC+z8KWTQxJrlaclCPF/alXI9o4g2g+/Qh+OjfkDkPzn0OmnQOb1wiIjGitAVdns9/HFg6tVXg6UrnXHaoA5PyC3tSB4iLh1NGQ7Ne8NpV3qIw/e6G468Hny0JanoXEakYP6Pf+wCrgceBJ4BVsbigS1UQ1ub3oo450WuOb9kfZv0FXjofdm31/XaNgBcRKT8/A+X+A/R3zp3inDsZGAA8HNqwpCIioqZeWI16MPJFGPQArPsI/nsCrJoV7qhERKKWn6Se4Jxbmf/EObcKSAhdSFJhEZbTAa/Jvec13lKutRp7q729dQvs313iW1RLFxGpGD9Jfa6ZPWtmfQLbMxwYRCcRJOJq6oU1bgtXz4YT/+Dd0/7USZA5v8TiSuwiIuXnJ6lfBywDbgxsywL7JMJEdFIHbynX/v/0JqzJ3gPj+sEnD0Be7iFFI/5cREQiUKk3SptZPPCcc+5C4KHKCUkqKqIGypUm/WS47jN4+0/ewjCr34fhT0HdZuGOTESkSis1qTvncs3sGDNL1CIukS/O3/o8kSG5rncPe6tBXnL/b28Y8C+m2WncO3M5P+3IpknK94we0JphnVPDHa2ISJXgZ0qzdcBnZvYG8Fv+Tuecau6RpopU1A/S8Tw4+niYdh28eSP1XAa270oc9cnctofbpy4BUGIXEfHBT9VuLfBWoGztQptEmDirQjX1wlKawiVv8J9qV9ON5bxb/VbOifsEcOzJzuWBWSvLPISIiPioqTvn7gYwsyO8p04TekeoKtX8XlRcHI/tOpXp1pYHEp7iP4lPMij3K/6SfRWbtdKAiIgvfmaU62ZmS4DFwBIzW2RmXUMfmpRXlRkoV4ImKcl85xozcv/fuDv7YnrFfcu71f+Py2p/fVirvomIxAo/VbvngOudc82cc82AG4DxIY1KKqSq3wY2ekBrkhPiccQxPncQg/ffy3pSuTP7EZhyEez8MdwhiohEND9JPdc592n+E+fcHCAndCFJRVX1mvqwzqncO7wDjeskYMD+Os3ZOPQ16HcPrH4PHu8O819QrV1EpARW1opYZvYIkAxMAhxwPrAXmAjgnCt5WrAQ6Natm5s7VxPaFefnPT/TILlBuMM4bJm7MtmTvYcWdVsc2PnzanjjRvjuc0g/Bc56FOqlhy9IEZFKYmbznHPd/JT1c0tbp8DPO4vs74yX5PuWIzaRMhl2aKtDg5Zw2dswbzy8dyc8cQL0/Sv0vA7i/fwai4hEPz+j30+tjEDk8FX1PvV8JZ5HXBx0vxJaDfQmrHn3b7D0NRgyFo7sULlBiohEoCp8D5QUVWXvUy9GqX+g1EmFCybBueNh+yZ4ug988A/I3ltp8YmIRKLoyQISNTV1oOzZ8cyg/XC44WvoMAI+/Y+3Xvva2ZUSnohIJFJSjyJVffR7vnKdR416cPZ/4eJpgMELZ8OrV8DOLSGLT0QkUvmZfOY8M6sdePw3M5tqZl1CH5pIORx7Klz3OfS5HZa/CY91h6+fKXZZVxGRaOWnpv5359xOM+sNnA6MA/4b2rCkIqKlpg4V7EpISII+t8F1X0BqF5jxZ3j2NNi8IPgBiohEIF+TzwR+ngE87Zx7G0gMXUhSUdHSp37Y59Gghdccf8442LEZnukLM0bD3u3BCVBEJEL5SeqZZvYU3qQzM8ysus/3SSWLlqQOQTgXM+hwrjeQrtuVXlP82G6w8CXIywtOkCIiEcZPch4BzAIGOOe2AfWA0SGNSiokWprfzYqZfKaiklPgjAfh6tmQcrS3bvtz/dUkLyJRyU9SPwp42zm32sz6AOcBX4c0KqmQaKmpW+BfUKV2gSvfg6FPwK8b4OlTvWlnf8sK7ueIiISRn6T+GpBrZi2Ap4GmwEshjUoqJD4uPtwhBEXIWhzi4qDzhfCHeXD89bBgIoztDF89Dblao0hEqj4/ST3POZcDDAfGOudG49XeJcJUj68e7hCCJqStDkl1YOC/vVvgjsqAd0bD06fAhjmh+0wRkUrgJ6lnm9kFwCXAW4F9CaELSYSyZ5QLhkbHwSXTYcT/vJHxE86AyRdC1tpK+HARkeDzk9QvB04A/uWcW29m6cALoQ1LYllI+tRL/DCDtkNh1DfQ9++w7iN4vCfMvB12/1I5MYiIBEmZSd05twy4FZgfeL7eOXd/qAOT2BWWUfwJyXDyn+EP8yHjd/DVkzCmM3zxBOTsr/x4REQqwM80sWcBC4GZgecZZvZGqAOT2Ba22/NqN4YhY+DaOdCkM8y6HZ7oCcvfAufCE5OIiE9+mt/vAnoA2wCccwuB5iGMSWJcpTa/l6RxO7j4dbjwVYhLgCkXen3uur9dRCKYr4Fyzrmi82tqSi4JmYhI6uD1t7fs542SP+Mh2LrSW7t96u/h143hjk5E5BB+kvq3ZvY7IN7MWprZWODzEMclMSyoM8oFQ3w16H4l3Dgfev0Rlk2Dx7rBO7fBrq3hjk5EpICfpP4HoB2wD2/Sme3AH/0c3MwGmtlKM1tjZreVUu4cM3Nm1s3PcUXCIqkO9LvbG0zXaSR8/RSMyYAP74W9O8IdnYgI5kI0+MfM4oFVQD9gE/ANcEFgNH3hcrWB/JXfRjnn5pZ23G7durm5c0stIlXc9n3b2ZOzhyNrHhnuUEq3dRV8+E9YNh1q1IeT/gzdrvCWgBURCRIzm+ec81Xp9TP6/T0zSyn0vK6ZzfJx7B7AGufcOufcfmAyMLSYcvcA9wN7/QQs0c8sQvrUy9KwlTdxzdUfwpEdvZHyj3WDBS9CXm7Z769iQlUBEJHg8dP83iCwOhsAzrlfgUY+3pcKfF/o+abAvgJm1gVoGlijvURmdo2ZzTWzuVu3qg8z2hkR1qdeltQucMk0b3a6mg1g+vXw3xOj6ja47LxsVv66MtxhiEgZfM39bmZH5z8xs2OAw/4/lZnFAQ8BfyqrrHPuaedcN+dct4YNGx7uR0uEi5jR7+XVvI9Xax/xP6+mPuVCb075le9U/eTuVFMXqQr8JPW/AnPM7AUzmwh8Atzu432ZeCu65UsL7MtXG2gPfGRmG4DjgTc0WE6qtPxpZ6//0lvmde92mDQSnjkVVs2qssndBf6JSGTzM03sTKALMAWvX7yrc85Pn/o3QEszSzezRGAkUDATnXNuu3OugXOumXOuGfAlMKSsgXIS/SLulraKiK/mLfM6ai4MecybR/6lEfBMX1j1bpVL7kroIlWDn4FyZ+NNQPOWc+4tIMfMhpX1vsByraOAWcBy4GXn3Ldm9g8zG3K4gUv0qpJN7yWJT4AuF3truA8ZC7t/hpfOg2dPh9XvV5nk7pxq6iJVQZm3tJnZQudcRpF9C5xznUMaWQl0S1v02529m905u2mQ3CDcoQRfzn5Y9BJ88iBs/x7SukOf2+HYvl7TfYTam7OXtdvX0q5+u3CHIhJzgnpLWwllqpUvJJHyiaraemHVEqHrZd4ENmc+DDt+gInD4dnTYMXbkKcZmEWk4vwk9blm9pCZHRvYHgLmhTowiV1R0adelmqJ3kQ1NwaS++4smPw771a4xS9Dbk64IxSRKsjvNLH78QbKTcGbLvaGUAYlEjOqVfeS+6h5MPxZrwl+6tUwtgvMfQ6yNSeTiPhXZjO6c+43oMR520WCrcrep3444qtBx/Og/TmwaiZ8+iC8dTN8dD+cOAq6Xg7Va4UtvKhvORGJEmUmdTP7kGImm3HO9Q1JRCKxLC4OjhsMrQfB+k/g0//Au3/zBtb1vBZ6/h5q1AtLaDH3h5ZIFeRnwNufCz1OAs4B1OEnIVNl5n4PJTNofoq3bZoHcx6Cj++Dzx717n8//nqof2y4oxSRCOOn+b3ooLjPzOzrEMUjIkWldYWRL8JPK+CLx2D+/+CbcXDcGXDiH6Bpz4i+HU5EKo+fyWfqFdoamNkAoE4lxCYxqsot6FJZGh0HQx+DPy6Fk/4EGz+D5wZ4E9l8+7pGzIuIr+b3eXh96obX7L4euDKUQYlIKWo3htP+DifdAgtfgi8eh1cug5RjvGb5zheFdVCdiISPn+b39MoIRCRfzPen+5VYE3pc7d0St3IGfD4WZt4KH/3bm+Cm+1WQcnSZhxGR6FFi87uZdTezIws9v8TMppvZGDMLz/BbETlUXDy0OQuufBeufA/ST/ES/KOdYPKF3ij6w5xjXnO/i1QNpfWpP4U36QxmdjJwH/A/YDvwdOhDk5hlui+6wpr2gPNfgJsWQ6+bYOPn8PxZ3kx1c8fD/t8qfmzldJGIV1pSj3fO/RJ4fD7wtHPuNefc34EWoQ9NRCospSmcfhfcssxb+jUuHt76IzzUBmb9FX7dUK7DqZYuUjWU1qceb2bVAkuongZc4/N9IoclJmeUC5WEZG/p184XwXdfwtdPwZf/9QbXtR7k9cmn9/EmvSnBtAWZ3D9zOVu276NJyo+MHtCaYZ1TK+8cRMS30pLzJOBjM/sZ2AN8CmBmLfCa4EWkqjCDY07wtu2Z3rzy8yZ4A+zqpnsD6zpfBDUPXu522oJMbp+6hD3ZuQBkbtvD7VOXACixi0SgUtdTN7PjgaOAdwNzwGNmrYBazrn5lRPiwbSeevTLzctlT84eaiXqtqyQyt4Ly9/0Evx3n0N8ojfgrtsVcEwvMKPXfbPJ3LbnkLempiTz2W2aKVqkMpRnPfVSm9Gdc18Ws29VRQMT8cPMqB5fPdxhRL+EJG8RmY7nebPVzRsPiybB0tegQSvoejm/basPHPrH1eZiEr2IhJ+fpVdFKlWcxZEQnxDuMGJLo+Ng0P1wywoY+gQk1YFZt/NV0ij+k/BfutpKCg9/b5KSHL5YRaREGvAmIgck1vAWjOl8IWxZQubMxxiwfjrnVP+UtXlH8XJuH2bE9eFPAzLCHamIFKPUPvVIpD51kcr15jerWDxrPP33v0/3uFXkWTxxLft7A+taDQC1qoiEVHn61JXURaRM2/dtZ932dXSOqwULJnp977t+hJoNoeP50PlirwlfRIKuPEldfeoi4l+DltDvbrh5GVwwxVv29asn4Yme8Mxp3kj6Pb+GO0qRmKU+dRHx5aBWvfhq0Hqgt+3aCktehvkvwFs3wzu3QquB0GkktOgH1RLDF7RIjFFSF5HDU6shnHCDt+zr5gWweAoseRWWvwHJ9aD9cOg4EtK6eZPgiEjIKKmLSHCYQWoXb+v/T1j7ISye7PXBf/Ms1Gvu9b93HOE9FpGgU1IXEV/KtXJefAK06u9te3d4tfZFk+Gj++Cje72++PbnQtuhULtx6IIWiTFK6iJSpsNaCjfpCO/2t84XwfZNsOQVWPwyvDMaZt4KzXpDu+HQZgjUrB+8oEVikG5pE5Ey7di/g3Xb1pHRKIiTzvy0HJZOhW+nQtYasHg49lQvwR93BiSnBO+zRKqwoM39LiIChGYp3EZtoO9f4dS/wJYl3pzz306F6dfDW4nQ4nQvwbceCNVrB//zRaKQkrqIlMmww2uCL/XgBkd19LbT74LM+V5yXzrVWxq2WhK07A/thnk/leBFSqSkLiJlMrPQ1NYP/SBI6+pt/e6B77/yEvy307zBdvHVvSb6NmdBq0HqgxcpQkldRMpkVFJSLywuDo45wdsG3ucl+OVvwvK3YNVMsDhv3fc2Q7w++DqplRufSATSQDkRKdOenD2s27aOdg3ahTsUcA5+WBRI8G/Czyu9/andoM2ZXpKvf2x4YxQJIi3oIiJBtTdnL2u3r6Vd/QhI6kVtXQUrAgl+8wJvX6O2Xu291SBo0tmr9YtUUUrqIhJU+3P3s3bbWtrUbxPuUEq37XtY8bbX//7dF+DyoGajwEQ4A6H5qVC9VrijFCkX3dImIkFlZsRZFajtpjSF46/1tt2/wJr3YeU7sOxNb7ra+ERodhK0HuStBZ9ydLgjFgkq1dRFpEw5eTms3baW1vVahzuUisnN9mruq2Z5Sf6Xtd7+Ru28++BbDYTUrhAXH944RYqh5ncRCao8l8fabWtpWbdluEMJjp9XeyPoV84MNNPnQo36cGxfb9KbY/tCrUbhjlIEUPO7iARZnMWREJcQ7jCCp0FLbzvxD7DnV1jzgVeLXzvbm5se4KhOcOxpXpJv2sNbpEakDHkuL6xdVaqpi4jky8uDLYu8JL/mA+/eeJcLibWh+SnQ4jQv0dc9JtyRSoTauX8ntRODO+uhauoiIhURF+fdAtekM5z8Z9i7HdZ/4g24W/MBrHjLK1e/pZfgm/fxJsBJOiKcUUsEcYS3oqykLiJSkqQ63pS0bc7yJr35eVWgFv8+zJsAXz3prS6X2tWryaef4jXVV6se7sglTMLd+q3mdxGRisje6zXPr/8Y1n0Mm+d798VXS4ajj/dq8c1PgSM7alR9DPl176/UTaob1GOq+V1EJNQSkryk3fwUOA2vqX7DZweS/Pt3euWSUiD9JK8W3+wkaNjaW7hGopKa30VEokFSHThusLcB7Nzi9cev+9hL9Mvf9PbXaADHnAjNenv98Y3aahrbKJLn8sL6+UrqIiKhUPtI6DjC25yDX9bBxs9g4+dejX75G165pBQvyR/TC5r1UnN9FRfuLm0ldRGRUDPzVo6rfyx0ucTbt+07L7lvnOP9XDnD21/9CK9P/pgT4egToUmGBt5VIVFdUzezgcCjQDzwrHPuviKv3wJcBeQAW4ErnHMbQxmTiEhESDkaMo6GjAu85zs2H0jyGz+H1e96++MTvVvsmvb0kn1aD6jVMHxxS6nCndRDNvrdzOKBVUA/YBPwDXCBc25ZoTKnAl8553ab2XVAH+fc+aUdV6PfRSQm7PrJG13/3Zfez80LIS/be63esYEk39P72aC1+uUjROauTFJrpQb1mJEy+r0HsMY5ty4Q1GRgKFCQ1J1zHxYq/yVwUQjjERGpOmo1OnCPPHi30G1e4CX477+C1bNg0Uvea0l1vBp8fpJP7QqJNcMXewwLd009lEk9Ffi+0PNNQM9Syl8JvFPcC2Z2DXANwNFHa6lEEYlBCUlwzAneBt7gu6y18H2gJv/dV7DmPe81i/NG1ad2PbA1PA7iNYwq1DRQDjCzi4BuwCnFve6cexp4Grzm90oMTUQkMplBgxbe1jnQyLn7F9j0jbdlzoNl02D+895rCTW8vvnULgcSfZ2mumc+yHJdblg/P5RJPRNoWuh5WmDfQczsdOCvwCnOuX0hjEdEJLrVqAetBngbHLiVLnMebJrr/fzqKcjd771es1Gh2nwXb0sO7mxosSaaa+rfAC3NLB0vmY8Efle4gJl1Bp4CBjrnfgphLCIisafwrXQdR3j7cvbDj0u9BJ85HzLnwqpCPZ91m3nLzhZsGVCzQVjCr4qitqbunMsxs1HALLxb2p5zzn1rZv8A5jrn3gAeAGoBr5jXBPSdc25IqGISEYl51RIP1Mrz7d3uDcLLnA8/LPK2ZdMPvH5Empfgm2QcSPa1j6z82KuAcA+U04IuIiJyqD3bYMviA0l+80LIWgP5c5vXauzV4vOT/JHtIeWYmO+jX561nDb12wT1mJFyS5uIiFRVySmQfrK35du3E7YsPZDof1jojbjPr50m1obG7bwE37gdNO4AjdpA9VrhOYcwyCMP5xwWpj9ulNRFRMSf6rUPvq0OYP9u+GkZbFkCP37r9dcvfhn27QgUMKiXfiDJ5yf9aK3VO2+lNkNJXUREqprEGpDWzdvyOefNbZ+f5H9c6tXwl79FQfN9fq2+URtva3ic97Nmwyqd7PNcHnkujzgLzwx/SuoiIhJcZlD3GG/LX4oWYP9v8NPyA0n+x6Xw7eswb/yBMsn1Agn+OGjY5sDPKjLfvQv8CxcldRERqRyJNYuv1e/60Uv2W1d4208rYOlr3qj8fDXqe8m9YesiNfvIut0uz+WF9V51JXUREQkfM+/2uNpHwrGnHtjvHOzcAluXe0l+63LYuhKWvAr7iiT7+i0Ds+u1Cjxu6d1vH59Q6acDqKYuIiJyEDM44ihvO7bvgf3Owc4fCtXsV3q32q16FxZMPFAurhrUTfcSfP1Awm/Q0kv6NeuHLGyHU01dRETEFzM4oom3tTjt4Nf2bPMS/M+r4edVkLUafl4Da94/MDUueP32+Qm+QQuo19zb6qYf9u13zqlPXURE5PAlpxzaZw+QlwvbNnoJPiuQ8H9e491jv3DiwWVrNT6Q5OulF3rc3FvitgyqqYuIiIRSXPyBxEz/g1/buwN+Xe8tfFOwrYe1s2HhDweXrVH/4CRfUMNvBjXqM23hZv454yeydm6hSUoyowe0Zljn1Mo6S0BJXUREYlnSEQemui1q/2/w64YiCX8dbPzcm2CnUDN7drWatM5uwD/zGjEpri8fb+vE7VOXAFRqYldSFxERKU5izcBMeO0OfS17r9ekn7UWtm3k9fc+pV7eZprbZlLYCcCe7FwemLVSSV1ERCSiJSR598w3bA3ArdOOLnZ43OZteyo1rPDMYyciIhJFmqQkl2t/qCipi4iIHKbRA1qTnBB/0L7khHhGD2hdqXGo+V1EROQw5feb/3PGYrJ25mn0u4iISFU2rHMqqUdtpk29NtRKDM8a8mp+FxERCSIL49KxSuoiIiJBYoF/4aKkLiIiEiThrKWDkrqIiEjQGKbmdxERkWhgpuZ3ERGRqKCauoiISBRRTV1ERCQKmKmmLiIiEhXiiFNNXUREJBpooJyIiEiU0EA5ERGRKKGauoiISJSIt3jV1EVERKJBnIU3rSqpi4iIBEl8XHxYP19JXUREJEjiTUldREQkKtRLqhfWz1dSFxERCZLaibXD+vlK6iIiIlFCSV1ERCRKKKmLiIhECSV1ERGRKKGkLiIiEiWU1EVERKKEkrqIiEiUUFIXERGJEkrqIiIiUUJJXUREJEooqYuIiEQJc86FO4ZyMbOtwMYgH7YB8HOQjxkO0XIeoHOJRNFyHhA95xIt5wE6l9Ic45xr6KdglUvqoWBmc51z3cIdx+GKlvMAnUskipbzgOg5l2g5D9C5BIua30VERKKEkrqIiEiUUFL3PB3uAIIkWs4DdC6RKFrOA6LnXKLlPEDnEhTqUxcREYkSqqmLiIhECSV1ERGRKBHTSd3MBprZSjNbY2a3hTue8jKzDWa2xMwWmtncwL56Zvaema0O/Kwb7jiLY2bPmdlPZra00L5iYzfPmMB1WmxmXcIX+cFKOI+7zCwzcF0WmtngQq/dHjiPlWY2IDxRF8/MmprZh2a2zMy+NbObAvur1HUp5Tyq3HUxsyQz+9rMFgXO5e7A/nQz+yoQ8xQzSwzsrx54vibwerNwxp+vlPOYYGbrC12TjMD+iPzdKszM4s1sgZm9FXgeGdfEOReTGxAPrAWaA4nAIqBtuOMq5zlsABoU2ff/gNsCj28D7g93nCXEfjLQBVhaVuzAYOAdwIDjga/CHX8Z53EX8OdiyrYN/J5VB9IDv3/x4T6HQvEdBXQJPK4NrArEXKWuSynnUeWuS+C7rRV4nAB8FfiuXwZGBvY/CVwXeHw98GTg8UhgSrjPoYzzmACcW0z5iPzdKhLjLcBLwFuB5xFxTWK5pt4DWOOcW+ec2w9MBoaGOaZgGAo8H3j8PDAsjLGUyDn3CfBLkd0lxT4U+J/zfAmkmNlRlRNp6Uo4j5IMBSY75/Y559YDa/B+DyOCc+4H59z8wOOdwHIglSp2XUo5j5JE7HUJfLe7Ak8TApsD+gKvBvYXvSb51+pV4DQzs0oKt0SlnEdJIvJ3K5+ZpQFnAM8GnhsRck1iOamnAt8Xer6J0v/Dj0QOeNfM5pnZNYF9jZ1zPwQebwEahye0Cikp9qp4rUYFmg2fK9QFUmXOI9BE2BmvRtQL3xsAAAd8SURBVFVlr0uR84AqeF0CzbwLgZ+A9/BaErY553ICRQrHW3Augde3A/UrN+LiFT0P51z+NflX4Jo8bGbVA/si+poAjwD/B+QFntcnQq5JLCf1aNDbOdcFGATcYGYnF37Ree09VfKexaocO/Bf4FggA/gB+E94wykfM6sFvAb80Tm3o/BrVem6FHMeVfK6OOdynXMZQBpeC8JxYQ6pQoqeh5m1B27HO5/uQD3g1jCG6IuZnQn85JybF+5YihPLST0TaFroeVpgX5XhnMsM/PwJeB3vP/gf85upAj9/Cl+E5VZS7FXqWjnnfgz8DywPeIYDTbkRfx5mloCXCF90zk0N7K5y16W486jK1wXAObcN+BA4Aa85ulrgpcLxFpxL4PU6QFYlh1qqQucxMNBV4pxz+4DxVI1r0gsYYmYb8Lpt+wKPEiHXJJaT+jdAy8CIxUS8AQxvhDkm38ysppnVzn8M9AeW4p3DpYFilwLTwxNhhZQU+xvAJYERsccD2ws1B0ecIn1/Z+NdF/DOY2RgNGw60BL4urLjK0mgn28csNw591Chl6rUdSnpPKridTGzhmaWEnicDPTDGyPwIXBuoFjRa5J/rc4FZgdaV8KqhPNYUeiPRcPrgy58TSLudwvAOXe7cy7NOdcML2/Mds5dSKRck1COwov0DW+E5Sq8Pqq/hjuecsbeHG/E7iLg2/z48fpqPgBWA+8D9cIdawnxT8JrAs3G63+6sqTY8UbAPh64TkuAbuGOv4zzeCEQ52K8/6CPKlT+r4HzWAkMCnf8Rc6lN17T+mJgYWAbXNWuSynnUeWuC9ARWBCIeSlwR2B/c7w/PNYArwDVA/uTAs/XBF5vHu5zKOM8ZgeuyVJgIgdGyEfk71Yx59WHA6PfI+KaaJpYERGRKBHLze8iIiJRRUldREQkSiipi4iIRAkldRERkSihpC4iIhIllNRFRESihJK6iIhIlFBSFzkMZpYbWAd6qZm9WWjWrPz9i8xsvpmdGO5YK4OZ7Sq7VLmO9//bu7cQq6o4juPfXwPmhIbp5EMEDURQjjbDmA+RL+VLUEJgIWE3eghK6EIXjAzCigoLIoqw8sEMIgQzQqqpKbREcWg0L2lCFx+KCs1CcUyb/j2sNbnncM7MOWcOzMzh93maWWet/1r7vPxnrb1n/2dIuq/we7sKtetHGdsqaYuklhrmmyJpa+F1n2aTipO62dgMRERXRMwllWBdXtLeSSpa8dy4rXBym0GqR12Pu4GNETFY7YBIZZh7gaV1zmk2rpzUzRpnO+XLQ54PHCttzO/v35x38/skLS18dpuknXm3v2ZotynpCUmHJH0l6V1Jj5TuXnPbUyPFymMOSHpT0n5JPfmd3Ei6I5fC/EbS+tHWVEkd8z4p6bvitQHPA5fmGKtz6JZy48tYRqH2gaQNkl7N8Q9LWihpff4+1xbGbcpjzSYdJ3WzBsgJbhFniwK15kR0EHgLeLrMsOuBXyKiM+/0P86xriDtFK+JVKpyEFgmaT6pgEQX6V3mC6pYV9lY+ePLgNciogP4E1giqQNYCVyXTxkeqCJOI+ZdACwBOkmlhK/KfVcA3+dTj0crjS8z/xTSO7Z/KjTPA36IiIXAGlLRl8eAOcANOlvLex9VfLdmE5HvG5mNTauk3aQd+gHg09w+kJMZkq4G3pY0N4YXW9gLvCTpBVJRiC9z+yJgPtCXilfRSip3OhN4PyJO5rjVVBWsFGsr8GNE7M79vgbagQuADRFxBCAi/hglTqPmbQM+iIhTwClJH44Qu9z4Um2khA+ApKmko/yXc1MAayNX/pI0CJzO1zwo6bSk6RFxfIR1mE04TupmYzMQEV2SzgM+Id1Tf6XYISK2S2oDLqSQCCPikKRu0q77GUm9EbGKVKFqXUQ8Xowj6cEKa/iH4aduU4vDKsRqB/4uNA2SEm8lZePU2r+OecupZvwAw7+HDqA/Ui11SCcCr+c1XUw6MSn+wXUucKrGdZmNOx+/mzVA3j3fDzxc+uS0pMuBFuBoSftFwMmIeAdYDXTnj3qBmyXNzv1mSrqEtMu9KT/VPR1YnPv/BsyWNCsfId9YmKZSrEo+B26RNGuof51xau2/DVgsaaqkaYVrOA5MH2FcWRFxjHTvfSixzyOVKR5yJakMKKQEP/Qz+dqPRMSZWuc1G2/eqZs1SETskrQHuJWzx/KQdq13lnkKex6wWtK/pHrs9+Y430paCfRIOid/tjwidkh6j5Scfgf6cv8zklaRajX/DBwsrKlsLODXCtewX9KzwJZ8JL0LuGuEOIcrxKl13r58O2EP6Y+UvcBfEXFU0jalBwE/ItXYrlYPqbb6Z6Tveif8fxTfmhM/DE/wANcCm2uYx2zCcD11s0lK6Qn3ExHx4nivpREkTYuIE/lWxlbgnojoH0O8buChiLi9xnEbgRURcajeuc3Gi3fqZjZRvCFpDule+LqxJHSAiOiX9IWklmr/Vz0/Nb/JCd0mK+/UzczMmoQflDMzM2sSTupmZmZNwkndzMysSTipm5mZNQkndTMzsybhpG5mZtYknNTNzMyaxH9RIUiXkEdE/gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x113b625f8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# This requires matplotlib.\n",
    "%matplotlib inline\n",
    "rbresults.plot()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Getting at the estimates\n",
    "The various estimated paramaters (e.g., $r$) are stored in the `rbresults.fits` dictionary. Here, because we have implemented a standard practice analysis, there are two keys: `'full'` corresponding to the full fit, and `'A-fixed'` corresponding to the fit with the asymptote $A$ fixed to the value specified in the analysis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "dict_keys(['full', 'A-fixed'])"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rbresults.fits.keys()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Each element of `rbresults.fits` is a `FitResults` object, which is just a container for the results of a fit (and a few simple methods. The most important information is the estimates of $A$, $B$, $p$ and $r$, accessed as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'A': 0.05253778812382854,\n",
       " 'B': 0.7691985059092118,\n",
       " 'p': 0.9945673921330701,\n",
       " 'r': 0.005427302585809846}"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rbresults.fits['full'].estimates"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that even when a parameter was fixed to a constant, it will appear in the estimates list. To check which parameters were actually variables look at:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'A': True, 'B': True, 'p': True}"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rbresults.fits['full'].variable"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can look at the bootstrapped standard deviations:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'A': 0.029585903246459568,\n",
       " 'B': 0.02763279399963954,\n",
       " 'p': 0.0005596198798488627,\n",
       " 'r': 0.0005590733760599479}"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rbresults.fits['full'].stds"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And we can access the distribution of the bootstrapped values for any of these parameters, most importantly $r$. Below we histogram this distribution:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD8CAYAAABn919SAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvNQv5yAAADRNJREFUeJzt3W2MpeVdx/Hvz6WAoVqeJmTDls7WEhreCM1K2mAaA1aBbQATYiBq9gVmX9gmbappF/umGk0WE6W+MDUIlI3VQsUHCNtEEWiMiaEO5amAyEKXdAmwWy3W+oKG8vfF3JBhmdlzZuacGc6f7yeZzP18/a9cs7+95j7n3JOqQpI0+35iswuQJE2GgS5JTRjoktSEgS5JTRjoktSEgS5JTRjoktSEgS5JTRjoktTEcRvZ2Omnn17z8/Mb2aQkzbwHH3zwe1U1N+q4DQ30+fl5FhYWNrJJSZp5SZ4b5zhvuUhSEwa6JDVhoEtSEwa6JDVhoEtSEwa6JDVhoEtSEwa6JDVhoEtSExv6SVFplPk9+zel3YN7d25Ku9IkOUOXpCYMdElqwkCXpCYMdElqwkCXpCYMdElqwkCXpCYMdElqwkCXpCYMdElqYuxAT7IlyUNJ7h7Wtyd5IMmBJLcnOX56ZUqSRlnNDP1TwJNL1q8HbqiqDwDfB66dZGGSpNUZK9CTbAN2AjcN6wEuAu4YDtkHXDmNAiVJ4xl3hv5F4LPAa8P6acDLVfXqsH4IOHO5E5PsTrKQZOHIkSPrKlaStLKRgZ7k48DhqnpwLQ1U1Y1VtaOqdszNza3lEpKkMYzzPPQLgcuTXAacCPw08KfAyUmOG2bp24Dnp1emJGmUkTP0qrquqrZV1TxwNXBfVf0acD9w1XDYLuDOqVUpSRppPe9D/xzwmSQHWLynfvNkSpIkrcWq/gRdVX0D+Maw/CxwweRLkiSthZ8UlaQmDHRJasJAl6QmDHRJasJAl6QmDHRJasJAl6QmDHRJasJAl6QmDHRJasJAl6QmDHRJasJAl6QmDHRJasJAl6QmDHRJasJAl6QmDHRJamJVf4JO6mp+z/5Na/vg3p2b1rZ6cYYuSU0Y6JLUhIEuSU0Y6JLUhC+Kvo1t1gt1vkgnzSZn6JLUhIEuSU0Y6JLUhIEuSU0Y6JLUhIEuSU0Y6JLUhIEuSU0Y6JLUhIEuSU0Y6JLUhIEuSU0Y6JLUhIEuSU0Y6JLUxMhAT3Jikm8meSTJ40l+b9i+PckDSQ4kuT3J8dMvV5K0knFm6K8AF1XVzwLnAZck+TBwPXBDVX0A+D5w7fTKlCSNMjLQa9EPh9V3DV8FXATcMWzfB1w5lQolSWMZ6x56ki1JHgYOA/cAzwAvV9WrwyGHgDNXOHd3koUkC0eOHJlEzZKkZYwV6FX146o6D9gGXAB8cNwGqurGqtpRVTvm5ubWWKYkaZRVvculql4G7gc+Apyc5PU/Mr0NeH7CtUmSVmGcd7nMJTl5WP5J4GPAkywG+1XDYbuAO6dVpCRptONGH8JWYF+SLSz+B/C1qro7yRPAbUn+AHgIuHmKdUqSRhgZ6FX1KHD+MtufZfF+uiTpbcBPikpSEwa6JDVhoEtSEwa6JDVhoEtSEwa6JDVhoEtSEwa6JDVhoEtSEwa6JDVhoEtSEwa6JDVhoEtSEwa6JDUxzvPQ9Q4zv2f/ZpcgaQ2coUtSEwa6JDVhoEtSEwa6JDVhoEtSEwa6JDVhoEtSEwa6JDVhoEtSEwa6JDVhoEtSEwa6JDVhoEtSEz5tUdpkm/V0y4N7d25Ku5oeZ+iS1ISBLklNGOiS1ISBLklNGOiS1ISBLklNGOiS1ISBLklNGOiS1ISBLklNjAz0JO9Ncn+SJ5I8nuRTw/ZTk9yT5Onh+ynTL1eStJJxZuivAr9dVecCHwY+keRcYA9wb1WdDdw7rEuSNsnIQK+qF6rqW8Py/wJPAmcCVwD7hsP2AVdOq0hJ0miruoeeZB44H3gAOKOqXhh2vQicscI5u5MsJFk4cuTIOkqVJB3L2IGe5N3A3wKfrqofLN1XVQXUcudV1Y1VtaOqdszNza2rWEnSysYK9CTvYjHM/6qq/m7Y/FKSrcP+rcDh6ZQoSRrHOO9yCXAz8GRV/cmSXXcBu4blXcCdky9PkjSucf5i0YXAbwCPJXl42Pa7wF7ga0muBZ4DfnU6JUqSxjEy0KvqX4GssPviyZYjSVorPykqSU0Y6JLUhIEuSU0Y6JLUhIEuSU0Y6JLUhIEuSU0Y6JLUhIEuSU2M89H/d7T5Pfs3uwRJGoszdElqwkCXpCYMdElqwkCXpCYMdElqwkCXpCYMdElqwkCXpCYMdElqwkCXpCYMdElqwkCXpCYMdElqwkCXpCYMdElqwkCXpCYMdElqwkCXpCYMdElqwkCXpCYMdElqwkCXpCYMdElqwkCXpCYMdElqwkCXpCYMdElqwkCXpCYMdElqYmSgJ7klyeEk316y7dQk9yR5evh+ynTLlCSNMs4M/VbgkqO27QHuraqzgXuHdUnSJhoZ6FX1L8B/H7X5CmDfsLwPuHLCdUmSVmmt99DPqKoXhuUXgTMmVI8kaY2OW+8FqqqS1Er7k+wGdgOcddZZ621O0oTM79m/aW0f3Ltz09rubK0z9JeSbAUYvh9e6cCqurGqdlTVjrm5uTU2J0kaZa2Bfhewa1jeBdw5mXIkSWs1ztsWvwr8G3BOkkNJrgX2Ah9L8jTwi8O6JGkTjbyHXlXXrLDr4gnXIklaBz8pKklNGOiS1ISBLklNGOiS1ISBLklNGOiS1ISBLklNGOiS1ISBLklNGOiS1ISBLklNGOiS1ISBLklNGOiS1ISBLklNGOiS1ISBLklNGOiS1ISBLklNjPybopLUxfye/ZvS7sG9OzekHWfoktSEgS5JTRjoktSEgS5JTfiiqKQNt1kvTnbnDF2SmjDQJakJA12SmpiZe+jec5OkY3OGLklNGOiS1ISBLklNGOiS1ISBLklNGOiS1ISBLklNGOiS1ISBLklNGOiS1MS6Aj3JJUmeSnIgyZ5JFSVJWr01B3qSLcCfAZcC5wLXJDl3UoVJklZnPTP0C4ADVfVsVf0IuA24YjJlSZJWaz2Bfibw3SXrh4ZtkqRNMPXH5ybZDeweVn+Y5KlptzlBpwPf2+wipqBjv+zTbHhH9inXr7uN941z0HoC/XngvUvWtw3b3qSqbgRuXEc7mybJQlXt2Ow6Jq1jv+zTbLBP07WeWy7/DpydZHuS44GrgbsmU5YkabXWPEOvqleTfBL4R2ALcEtVPT6xyiRJq7Kue+hV9XXg6xOq5e1oJm8VjaFjv+zTbLBPU5Sq2uwaJEkT4Ef/JamJ9oE+6vEESU5Icvuw/4Ek80v2XTdsfyrJLx913pYkDyW5e8m2W5N8J8nDw9d5s9KnJAeTPDbUvbBk+6lJ7kny9PD9lAZ9+kKS55eM02Uz1KeTk9yR5D+SPJnkI8P2WR6nlfo0k+OU5JwlNT+c5AdJPj3sm+44VVXbLxZfrH0GeD9wPPAIcO5Rx/wW8OfD8tXA7cPyucPxJwDbh+tsWXLeZ4C/Bu5esu1W4KpZ7BNwEDh9mfb+CNgzLO8Brm/Qpy8AvzOj47QP+M1h+Xjg5AbjtFKfZnacjrr+i8D7NmKcus/Qx3k8wRUs/kAB3AFcnCTD9tuq6pWq+g5wYLgeSbYBO4GbNqAPR5tKn45h6bX2AVdOoA9H2+g+bYSJ9ynJe4CPAjcDVNWPqurlZa41M+M0ok8bYdo/excDz1TVc8tca+Lj1D3Qx3k8wRvHVNWrwP8Ap40494vAZ4HXlmnzD5M8muSGJCesuwdvNa0+FfBPSR7M4qd7X3dGVb0wLL8InDGJTqxU7zJ1veWYCfQJ4JPDON0ypdsT0+jTduAI8OUs3u67KclJwzGzOk7H6hPM5jgtdTXw1SXrUx2n7oE+cUk+DhyuqgeX2X0d8EHg54BTgc9tZG3r9PNV9SEWn575iSQfPfqAWvw9cZbeFrVSn74E/AxwHvAC8MebVN9qHQd8CPhSVZ0P/B+Lv7a/yYyN07H6NKvjBEAWP3B5OfA3y+2fxjh1D/RxHk/wxjFJjgPeA/zXMc69ELg8yUEWfz27KMlXAKrqhVr0CvBlpvOr/zT6RFW9/v0w8PdLan8pydbhWluBwxPsy1vqPbqu5Y5Zb5+q6qWq+nFVvQb8BbMzToeAQ1X1wLD9DhbDEGZ3nFbs0wyP0+suBb5VVS8t2TbVceoe6OM8nuAuYNewfBVw3/A/513A1cMr3NuBs4FvVtV1VbWtquaH691XVb8ObwwQw/21K4Fvz0KfkpyU5KeG2k8CfmlJ7UuvtQu4c9b79Po4DX6FGRmnqnoR+G6Sc4ZzLgaeWOZaMzNOx+rTrI7TkvOu4c23W46+1uTHaZKvsL4dv4DLgP9k8RXozw/bfh+4fFg+kcVfiQ4Mg/H+Jed+fjjvKeDSZa79C7z5XS73AY+x+IP3FeDds9AnFl/hf2T4evz1aw77TgPuBZ4G/hk4tUGf/nIYp0dZ/Ae2dRb6NGw/D1gYav8H4JRZHqcRfZrlcTqJxVn8e45qa6rj5CdFJamJ7rdcJOkdw0CXpCYMdElqwkCXpCYMdElqwkCXpCYMdElqwkCXpCb+Hxs39lAXHOVjAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x114491a58>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "% matplotlib inline\n",
    "\n",
    "plt.hist(rbresults.fits['full'].bootstraps['r'], 10)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Because we are doing this analysis on simulated data, we know what the DRB error rate should be according to the theory of DRB from [\"Direct randomized benchmarking for multi-qubit devices\"](https://arxiv.org/abs/1807.07975)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The theory-predicted r is: 0.004990009995000988\n",
      "The (simulated) observed r is: 0.005427302585809846 +/- 0.0005590733760599479\n"
     ]
    }
   ],
   "source": [
    "theory_predicted_r = 1 - (1-0.001)**5\n",
    "observed_r = rbresults.fits['full'].estimates['r']\n",
    "observed_r_std = rbresults.fits['full'].stds['r']\n",
    "print(\"The theory-predicted r is: {}\".format(theory_predicted_r))\n",
    "print(\"The (simulated) observed r is: {} +/- {}\".format(observed_r,observed_r_std))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}