{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Drift detection and characterization\n",
"\n",
"The drift code is currently part way through a major update. The core function works (see below), but the functions are either not documented or incorrectly documented, and most of this tutorial is commented out as it currently doesn't work.\n",
"\n",
"This tutorial shows how to implement drift detection, and a limited characterization, on time-stamped data. This data can be from almost any quantum circuit based experiment on one or more qubits, as long as the data is taken using a suitable time-ordering of the experiments, and the data is recorded as a time series. For example, possible experiments include suitably time-ordered GST, RPE, Ramsey or RB experiments.\n",
"\n",
"This notebook is an introduction to these tools, and it will be either augmented with further notebooks, or updated to be more comprehensive, at a later date."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"from __future__ import print_function\n",
"\n",
"# Importing the drift module is essential\n",
"from pygsti.extras import drift\n",
"\n",
"# Importing all of pyGSTi is optional, but often useful.\n",
"import pygsti"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## A quick overview\n",
"We now give a quick overview of the drift detection and characterization methods in the `drift` module. Further details are given later in this tutorial. As we demonstrate below, the analysis can be implemented with only two steps:\n",
"\n",
"1. Import some timestamped data, into a `pyGSTi` dataset.\n",
"2. Pass this data to a single analysis function, `drift.do_drift_characterization()`.\n",
"\n",
"Here we demonstrate this with time series GST data, on the $G_i$, $G_x$, $G_y$ model, generated from a simulation (the code required to run these simulations is not currently available in `pyGSTi`). In this simulation the $G_i$ gate has low-frequency drift, the $G_x$ has high-frequency drift, and the $G_y$ gate is drift-free (where \"low\" and \"high\" frequency are with respect to the sample rate). More details on the input data format are given later."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Loading ../tutorial_files/timeseries_data.txt: 100%\n"
]
}
],
"source": [
"ds = pygsti.io.load_tddataset(\"../tutorial_files/timeseries_data.txt\")"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
" - Formatting the data...complete.\n",
" - Calculating power spectra...complete.\n",
" - Implementing statistical tests for drift detection...complete.\n",
" - Estimating the probability trajectories...complete.\n"
]
}
],
"source": [
"# This takes 5 - 10 minutes\n",
"results_gst = drift.do_drift_characterization(ds)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Thats it!**\n",
"\n",
"Everything has been calculated, and we can now look at the results. \n",
"\n",
"**Is there any detectable drift?**\n",
"\n",
"One useful result is printed below: a yes/no outcome for whether or not drift is detected. This is calculated using multiple statistical tests on the data at a specified global confidence level (which defaults to 0.95 when no user-specified value is passed to `drift.do_basic_drift_characterization`). That is, here there is a probability of at most 0.05 that this function will report drift when there is none."
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"#results_gst.any_drift_detect()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**We can plot power spectra**\n",
"\n",
"These spectra *should* be flat - up to statistical flucations, due to finite-sampling noise, around the mean noise level - if there is no drift. There are a range of power spectra that we can plot, but the most useful for an overview of the data is the \"global power spectrum\", obtained from averaging power spectra calculated from the individual data for each of the different operation sequences (again, details on exactly what this is are given later). This is plotted below. If there are peaks above the significance threshold, this power spectra provides statistically significant evidence of drift."
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"results_gst.plot_spectrum()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**We can extract the drift frequencies**\n",
"\n",
"If we have detected drift, we would probably likely like to know the frequencies of the drift. This information can be extracted from the results object as shown below. All frequencies will be in Hz if the timestamps have been provided in seconds (again, details later). Note that these are the frequencies in the drifting outcome probabilities -- they are *not* directly the frequencies of drift in, say, a Hamiltonian parameter. However, they are closely related to those frequencies."
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"collapsed": false
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"[0.4 0.011 0.016 0.001 0.008 0.005 0.003 0.2 ]\n"
]
}
],
"source": [
"print(results_gst.get_drift_frequencies())"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Is there drift for a particular sequence?**\n",
"\n",
"There are individual power spectra for all of the sequences. E.g., if we are interested in whether the $G_xG_i^{128}G_y$ sequence shows signs of drift, we can plot the power spectrum:"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAABCwAAADMCAYAAAC4GTIBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJzs3WlgVOX59/HvLFkJJGRjVyCyqSgoIApWBRVUZBW3KlpbecS2tFYqiJZNENmqgiJSrfy1ooKgFESxgtYKguJKgIAEAgQCIQmB7JOZM8+LyZxMSIDsk+X3eZPMzJlzrtnPuc51X7fF7Xa7ERERERERERGpQ6z+DkBERERERERE5ExKWIiIiIiIiIhInaOEhYiIiIiIiIjUOUpYiIiIiIiIiEido4SFiIiIiIiIiNQ5SliIiIiIiIiISJ2jhIWIiIiIiIiI1DlKWIiIiIiIiIhInaOEhYiIiIiIiIjUOUpYiIiIiIiIiEido4SFiIiIiIiIiNQ5dn8HUBX5+fnEx8cTExODzWbzdzgiIiIiIiIicgaXy8WJEye49NJLCQ4OLvf96nXCIj4+nl//+tf+DkNEREREREREzuPtt9+mV69e5V6+1hIWAwYMIDAwkKCgIAAmTJjAtddey48//siUKVMoKCigTZs2zJs3j6ioqHKtMyYmBvA86JYtW9ZY7CIiIiIiIiJSOceOHePXv/61eQxfXrVaYbFw4UI6d+5sXjYMg7/+9a/Mnj2bXr16sXjxYubPn8/s2bPLtT7vMJCWLVvStm3bGolZRERERERERKquoq0c/Np0Mz4+nqCgILMk5O677+aTTz7xZ0giIiIiIiIiUgfUaoXFhAkTcLvdXHnllfzlL38hJSWF1q1bm7dHRkZiGAaZmZlERETUZmh+k5CUwY7ENLrHRdO1faS/wxERERERERGpE2otYfH222/TqlUrHA4Hs2bNYsaMGdx00021tfk6KSEpg6eWbKbQaRBgtzLrkX5KWoiIiIhIg1FYWEhycjL5+fn+DkVEaklwcDBt27YlICCgyuuqtYRFq1atAAgMDOTee+9l3LhxjBkzhqNHj5rLZGRkYLVaG011xY7ENByFBgBOp8GOxDQlLERERESkwUhOTqZp06a0b98ei8Xi73BEpIa53W7S09NJTk6mQ4cOVV5frfSwyM3NJSsrC/A8gPXr19OtWzcuvfRS8vPz2b59OwDvvvsugwcPro2Q6oTucdFYrZ4vbrvdSve4aD9HJCIiIiJSffLz84mKilKyQqSRsFgsREVFVVtVVa1UWKSnp/PHP/4Rl8uFYRjExcUxdepUrFYrc+fOZerUqSWmNW0suraPpHe3FnyXcFzDQURERESkQVKyQqRxqc7PfK0kLNq1a8eHH35Y5m1XXHEFa9eurY0w6qSmoYFYLBYlK0RERERERER81OosIVKa4XZjGG5/hyEiIiIi0uANGDCAwMBAAgMDycvL46KLLuLhhx/miiuuOO99V69eTc+ePatlXH5N2L17NwcOHODWW2/1dygi1aZWeljI2blcblxKWIiIiIiI1IqFCxfy73//m//85z+MGDGCsWPH8tNPP533fh988AFJSUk1H2Al7d69m08++eSstzudzlqMRqR6qMLCzwy3J1lhGG6zAaeIiIiISGOWkJTBjsQ0usdF1+jQ6Ztvvpmff/6Z119/nYULF/L111/zwgsvUFBQgMvl4pFHHuG2225j1apVxMfHM3PmTF544QUmTpxIVFQU06dPJy8vj4KCAu68804efPDBUttIT0/n8ccfJz09HYCrr76ayZMns3r1atauXUtQUBCHDh0iOjqaefPm0aJFCwCWLl3Kp59+isvlokWLFjzzzDPExMTgcDh4/vnn+d///ofVaqVdu3bMnDmThQsXkp2dzbBhw+jduzdPP/00Xbp04Q9/+ANffPEF1157LTabjdzcXCZOnAjAokWLzMuLFi1i//79ZGdnk5SUxCWXXMLYsWN57rnnOHr0KDfddJN5P5HaooSFn7kMo+ivEhYiIiIi0nBt2n6I/3xz6LzL5eYXcuDoadxusFigQ+tmhAYHnPM+N/W5gAG9LqhUXJdffjmbNm0C4OKLL2b58uXYbDbS0tIYOXIk/fv3Z9SoUXz44Yc89NBD3HDDDQBkZ2ezbNkyAgMDycnJYfTo0Vx77bXExcWVWP/atWu54IILWLZsGQCnTp0yb/vuu+/48MMP6dixIy+99BKzZs1i4cKFrFmzhsOHD7NixQqsVivLly/nueeeY8GCBSxdupTDhw+zevVqAgMDycjIoHnz5owfP54vvviChQsXlth+UFAQq1atAjwJinPZuXMnq1atIjQ0lBEjRrBgwQJee+01nE4nAwcO5K677qJ9+/aVep5FKkMJCz/z9q/wVlqIiIiIiDRmOXlOvLvGbrfn8vkSFlXh9tkPz8jIYPLkyRw8eBCbzcapU6c4cOAAPXr0KHW//Px8pk2bxp49e7BYLKSmppKQkFAqYXH55ZezbNky5syZQ58+fejfv79525VXXknHjh0BGD16NLfffjsAmzZtIj4+nhEjRgDgcrkICwsD4PPPP2fSpEkEBgYCEBl57goU7zrKo3///jRt2hSALl260LVrV7PnR4cOHTh06JASFlKrlLDwM2//CpfLgACbn6MREREREakZA3qVrwoiISmDp5Zsxuk0sNutTPj1lTU6LGTHjh106tQJgGnTpjFgwABeeuklLBYLgwYNoqCgoMz7/f3vfycmJobnnnsOu93OQw89VOayPXv25IMPPmDLli2sWbOGpUuX8s4775wzJrfbzbhx47jjjjuq/PhCQ0PN/202G0ZRhTdQKt6goKASy5552eVyVTkekYpQ000/Myss1HhTRERERISu7SOZ9Ug/7rulG7Me6VejyYrPPvuMd955h4ceegiArKws2rRpg8ViYfPmzRw8eNBctkmTJmRlZZmXs7KyaNmyJXa7nb1797J9+/Yyt3H48GHCwsK47bbbePLJJ9m5c6eZNPj+++/NRp6rVq2ib9++gGc2k+XLl5vDRxwOBwkJCQDccMMN/N///R8OhwPwVIUAhIWFlYivLBdeeKG5/ezsbL744ouKPF0itU4VFn5mVlgoYSEiIiIiAniSFjWVqBg/frw5rWlcXBxLly7l8ssvB+Dxxx9n+vTpLFq0iO7du9OlSxfzfnfddRfPPfccr7/+OhMnTmTcuHE88cQTvP/++3To0IHevXuXub1vvvmGZcuWYbVaMQyD6dOnY7V6zhtfccUVzJkzh4MHD5pNNwGGDx9OZmYm9913H+CpuLjnnnvo2rUrY8eOZcGCBQwfPpyAgAAuvPBCFi5cyNVXX80///lPhg4dSp8+fXj66adLxXLTTTexfv16brnlFlq3bs0ll1xSrc+tSHWzuN31t3lCcnIyAwcOZOPGjbRt29bf4VTKU69s5ud9abw5dRDNmwX7OxwRERERkWqze/duunXr5u8w6qTVq1eX2SRTpCE487Nf2WN3DQnxM5eaboqIiIiIiIiUoiEhfmaYTTeVsBARERERaSxGjhzJyJEj/R2GSJ2mCgs/07SmIiIiIiIiIqUpYeFnLreaboqIiIiIiIicSQkLPzNc3iEhxnmWFBEREREREWk8lLDwM+9QEBVYiIiIiIiIiBRTwsLPXIanskIVFiIiIiIiIiLFlLDwMzXdFBERERGpPadOneKyyy5j5syZ/g6lWk2aNIl//etfFb7fsmXLSE9Pr/Z4nnrqKbZv314t61q0aBFz5syplnWdKTk5mauuuqpG1i1Vp4SFn3mbbarppoiIiIhIzVu3bh2XX345H330EQ6Ho1rX7XQ6q3V9teHNN9+skYTFrFmz6NWrV7WvVxoXu78DaOzMCgslLERERESkAduSnM7m5Oo/MAbo1zaKa9pGlWvZVatW8de//pVXX32VjRs3csstt5CXl8f111/Pxx9/TGRkJABz5syhSZMm/OEPf+Cnn35i/vz55OTkADB+/Hiuv/56kpOTGTVqFCNHjmTr1q3ceeedtG/fnhdeeIGCggJcLhePPPIIt912GwD79u3jySefJC8vj65du3Lo0CHGjRvHDTfcQGpqKjNnzuTo0aMUFBRw22238cgjj5SK//vvv+eZZ57BMAycTifjxo1jyJAhAOzdu5cxY8Zw7NgxevTowZw5c7BYLKSlpTF16lQOHToEwG9/+1uGDx/OK6+8QmpqKuPHjycoKIgFCxZw0UUXldjeokWLOHDgAFlZWRw+fJgLLriAF198kZCQEHJycpg5cyY7duwAYNiwYTz88MMA3H///Tz00EPccMMNvPfeeyxbtozAwEAMw+CFF14gLi6O/fv38+yzz3Ly5EkKCwt54IEHGDVq1Hlfw6VLl/Lpp5/icrlo0aIFzzzzDGFhYZV6DaVuq/WExUsvvcSiRYtYu3YtnTt35scff2TKlCkUFBTQpk0b5s2bR1RU+b5sGgJVWIiIiIiI1I6EhAQyMzPp27cvJ06cYNWqVdxyyy2EhIRw4403sm7dOsaMGYPT6WTt2rW8++67nD59mqlTp7J06VJiY2NJTU3ljjvuYN26dQBkZmbSvXt3Jk6cCHiGnCxfvhybzUZaWhojR46kf//+hIeH88QTT/DAAw8wbNgwduzYwZ133mnGNnHiRB599FF69+6Nw+HgwQcfpHv37vTr16/EY/jHP/7Bb3/7W4YMGYLb7SYrK8u87ZdffmHZsmVYLBZGjBjBli1b6NevHzNnzqRTp068/PLLpKamMnLkSC6++GLGjRvHypUrWbhwIZ07dz7r8xYfH8/7779P06ZN+e1vf8vatWu58847Wbx4MYZhsHbtWnJycrjrrrvo3Lkz1113XYn7z507l48//pjY2FgcDgculwun08mECROYN28ecXFxZGdnM2rUKHr06EFcXNxZY1mzZg2HDx9mxYoVWK1Wli9fznPPPceCBQsq/RpK3VWrCYudO3fy448/0qZNGwAMw+Cvf/0rs2fPplevXixevJj58+cze/bs2gzLr1RhISIiIiKNwTUVqIKoKe+//z7Dhg3DYrFw8803M3PmTI4fP06LFi0YMWIEs2bNYsyYMXz55Zd07NiRtm3b8t///pfk5GSzcgDAYrFw8OBBmjdvTlBQELfccot5W0ZGBpMnT+bgwYPYbDZOnTrFgQMHuOiii9i7dy+33347AN27d6dLly4A5Obm8s0335CRkWGuJycnh8TExFIJi6uuuopXXnmFQ4cO0a9fPy6//HLzthtvvJGgoCAALr74YnOZr7/+mkmTJgEQGxvLddddx7Zt286ZpPDVv39/mjVrBsBll11mVmp8/fXXTJ48GYvFQlhYGLfddhtff/11qYRF3759mTRpEjfccAPXX3897dq1Y9++fSQmJvKXv/zFXK6wsJD9+/efM2GxadMm4uPjGTFiBAAul4uwsDCASr+GUnfVWsLC4XAwY8YMFixYwJgxYwBPpi4oKMgc23T33XczcODARpWwUIWFiIiIiEjNczgcrFu3jsDAQNasWQN4DpBXr17NuHHj6NWrFzk5OezZs4cPPviAkSNHAuB2u+nSpQtvv/12qXUmJycTEhKCxWIxr5s2bRoDBgzgpZdewmKxMGjQIAoKCszbfZf1MgwDi8XC+++/T0BAwDkfx4MPPsiAAQPYsmULzzzzDP369eOxxx4DMJMVADabDZfLVYFnyGP06NE4HA6aNGnC8uXLy1yv7+Mpj5deeokdO3awdetWxowZw7Rp02jdujXNmzc3X4vycrvdjBs3jjvuuKPUbZV9DaXuqrWmmy+++CJDhw6lbdu25nUpKSm0bt3avBwZGYlhGGRmZtZWWH6nCgsRERERkZq3ceNGOnTowJdffsmmTZvYtGkT//znP/nggw/MZYYPH84bb7zBt99+y6BBgwDo2bMnBw8eZOvWreZyP//8M+6zzPKXlZVFmzZtsFgsbN68mYMHDwIQFhZGp06dzGEIO3fuZO/eveZtV155JUuXLjXXk5KSwokTJ0qt/8CBA1xwwQXcfffdjBkzxuwfcS5XX301K1asAODEiRP897//pW/fvgA0adKkxLCSlStXsmbNGjNZcb71rlq1CrfbTXZ2NuvXr+eaa64psYzT6eTw4cNcdtlljB07ln79+rF79246dOhAcHAwH374oblsYmIi2dnZ59zmgAEDWL58OadOnQI8iaiEhATz9up4DaXuqJUKix9++IH4+HgmTJhQG5urV1xKWIiIiIiI1LhVq1aZwzG8evbsiWEYfPPNN/Tp04fhw4czcOBARo4cSUhICADh4eEsXryYefPm8eyzz1JYWEi7du1YsmRJmdt5/PHHmT59OosWLSox7AM8TSAnT57M0qVL6dy5M507d6Zp06YA5tB4b4xNmjRh1qxZxMTElFj/W2+9xbZt2wgICCAwMJCnn376vI/96aefZsqUKea6J0yYQKdOnQAYM2YMkydPJjg4uMymm+fy6KOP8swzz5jrHTp0KL/61a9KLGMYBpMmTSIrKwuLxUKrVq14/PHHsdvtLFmyhGeffZbXX38dwzCIiorihRdeOOc2hw8fTmZmJvfddx/gqZ6455576Nq1q3l7VV9DqTss7lpIKy1dupQ333yTwMBAAI4dO0ZUVBT3338///73v80sY0ZGBgMHDuSHH34o13qTk5MZOHAgGzduLFG5UZ/c8eQ6ChwuJj/Ym6u7tz7/HURERERE6ondu3fTrVs3f4dRZ+Tk5BAaGorFYmHfvn3cf//9fPLJJ4SHh/s7NJFqdeZnv7LH7rVSYTF27FjGjh1rXh4wYABLlizhoosuYsWKFWzfvp1evXrx7rvvMnjw4NoIqc5wubwVFn4OREREREREatQPP/zA3LlzzaEIzzzzjJIVIudQ69Oa+rJarcydO5epU6eWmNa0MTHc3qabyliIiIiIiDRk/fv3p3///v4OQ6Te8EvCYtOmTeb/V1xxBWvXrvVHGH7ndrvVdFNERERERESkDLU2S4iU5puj0LSmIiIiIiIiIsWUsPAjw2cYiCosRERERERERIopYeFHvlUVqrAQERERERERKaaEhR8ZSliIiIiIiNSaAwcOcNdddzFo0CDuuusukpKSylzO5XIxffp0brzxRm666SZWrlxp3rZo0SKuvvpqhg0bxrBhw5g+fXqF49i2bRsjR44873KrV69m/PjxFV7/pEmT+Ne//lXu5ZOTk3nvvffKtezp06f5xz/+UeGYfB0/fpz777+/SusAGDZsGPn5+QB8//33DBkyhOHDh7N161YefvhhDh06VOVt1ITdu3ezfv36Etd16dKFnJycGtvm/fffz+eff17h+50rrgEDBrB3796qhnZOSlj4kW+SQkNCRERERERq1tSpU7n33nvZsGED9957L1OmTClzubVr13Lo0CE+/fRT3nvvPRYtWkRycrJ5+/Dhw1mzZg1r1qxh6tSptRV+jTly5EiFEhavvfZalbbXokUL3nrrrSqtA2DNmjUEBweb/w8fPpwPP/yQvn378o9//IMLLrigytuoCbt37+aTTz6p1H2dTmc1R1O3KWHhR6qwEBERERGpHenp6ezatYshQ4YAMGTIEHbt2kVGRkapZdevX8/o0aOxWq1ERkZy4403VuoA8/nnn+emm25i9OjRzJs376xVFR9++CG33347t99+O7///e9JT083b8vKyuKRRx7h1ltvZcyYMRw/fhyAPXv2cO+99zJixAhuvfVWli1bdt548vLyGD9+PLfeeitDhw7lT3/6EwAzZswgMTGRYcOGmRUdc+bMYdSoUQwdOpQHHniAI0eOmMtmZWUxbNgw7r777nNuzzAMpk2bxuDBgxk6dKi5fHJyMldddZW53IYNGxg8eDDDhw9nyZIlJc7qd+nShSVLljBq1CgGDhzIhg0bzPt5l3vttdf4+OOPefPNN82qC9+z/8ePH+ePf/yj+Ry/+uqrgCcxNXr0aIYPH87w4cP5+uuvzXUPGDCAF198kbvuuosBAwaUqFhJTEzkoYceMtf3wQcfAJCamsr48eO54447uP3221myZEmp5+TkyZMsXLiQLVu2MGzYMGbOnGne9tZbb531cS5atIhRo0bx0ksvAbB06VLuuOMORowYwSOPPMKJEycA+Oyzz7j99tsZNmwYQ4YMYdu2beZ6vvnmG+655x4GDhzI/PnzzesPHjzIAw88wO23386IESP48ssvy3w9t2/fbj7mGTNm4HbX/DGsX6Y1FQ9VWIiIiIhIY5G66QuOb9xUI+tuMXAAsQOuP+cyKSkptGjRApvNBoDNZiM2NpaUlBQiIyNLLdu6dWvzcqtWrTh27Jh5+aOPPuKrr74iJiaGP/7xj/Ts2bPU9jZt2sTnn39uVgGcbWjH3r17mT9/PqtXryY2NpYXXniBZ555hhdeeAGA7777jg8//JCOHTvy0ksvMWvWLBYuXEibNm1YtmwZgYGB5OTkMHr0aK699lri4uLO+hx89dVX5OTkmMMRTp06BcCUKVOYM2cOq1evNpd9+OGHmThxIgArV65k/vz5PP/880yZMoVRo0axZs2asz/ZRRISEti2bRvr16/HarWa2/OVlpbGlClTeO+992jfvn2ZiZewsDBWrVrFd999x5///GcGDRpU4vbf/e537Nu3j0svvZT77ruv1P0nTJjAddddx6JFiwDMJFX//v0ZMmQIFouF/fv38+CDD5Y4WM/Pz+e9994jOTnZPJgPCgri0Ucf5c9//jO33HIL4ElCAEycOJFHH32U3r1743A4ePDBB+nevTv9+vUz19m8eXPGjx/PF198wcKFC8v9OIOCgli1ahXgqSY5fPgwK1aswGq1snz5cp577jkWLFjAwoULmTFjBj179sTlcpGXl2euIyUlhbfffpucnBxuvPFG7rjjDtq3b8+ECRO48847GT16NPv27ePXv/41H3/8cYnPhcPh4LHHHmP+/PlcddVVrF+/nrfffrvUc13dlLDwo5IVFsY5lhQRERERkbrg7rvv5pFHHiEgIIDNmzfz6KOPsn79epo3b15iuW3btnHLLbcQGhoKeIaRLF68uNT6tm3bxnXXXUdsbKy5/mHDhpm3X3nllXTs2BGA0aNHc/vttwOeg+lp06axZ88eLBYLqampJCQknDNh0bVrVxITE5k+fTp9+vTh+uuvP+uyX375JcuXLyc3N7fSwxDatWuH0+nkqaee4qqrruKGG24otcxPP/3ExRdfTPv27QEYNWoUs2fPLrHMrbfeCkCPHj1ITU2loKCAoKCgcsWQk5PDDz/8wBtvvGFe5z0QP3z4MI8//jjHjx/HbreTlpbGiRMniImJKbHdtm3b0qxZM44dO4ZhGDidTjNZAZ4kRG5uLt98802Jip2cnBwSExNLJCzO5VyPc8SIEeZymzZtIj4+3rzO5XIRFhYGQN++fZk9ezY333wzv/rVr+jcubN5v8GDB2O1WmnatClxcXEcOnSI6Ohodu/ezahRowC46KKL6NatGz/++CMDBgww77t//35CQkLMyphbb731rEOqqpMSFn5UosKiFsppRERERET8JXbA9eetgqhJrVq14vjx47hcLmw2Gy6Xi9TUVFq1alXmskePHuWyyy4DSlZceA9mAfr160erVq345Zdf6NOnT+08kCJ///vfiYmJ4bnnnsNut/PQQw9RUFBwzvu0a9eOdevWsXXrVr788kuef/551q5dW2q5I0eOMHv2bN5//33atWvH999/z4QJEyocY9OmTfnoo4/Ytm0bW7ZsYf78+ebwiYrwHrR7q2OcTme5Exbn8pe//IVJkyZx4403YhgGl19+eYnn0Hcb3veMxWIpc12GYWCxWHj//fcJCAioVDznepzexBeA2+1m3Lhx3HHHHaXWMXnyZPbs2cPWrVv505/+xG9+8xvuvPPOsz6eqjjbc1Gd1MPCj3wrLAyXEhYiIiIiIjUlKiqKbt26sW7dOgDWrVtHt27dSg0HAc+Z6JUrV2IYBhkZGXz22Wdmeb63hwR4miceOXKEDh06lFpHnz592LBhA3l5eRiGwb///e8y47rqqqv473//a/YgWLFiBddcc415+/fff2/OZrJq1Sr69u0LeHpbtGzZErvdzt69e9m+fft5n4Njx45hs9m48cYbefLJJ8nIyCAzM5OwsDCys7PN5bKzswkICCAmJgbDMHj33XfN28LCwsjPzy9X1UVGRgZ5eXlce+21TJgwgaZNm3L48OESy1x++eXs2rXLnNGjMgmNc2nSpAk9e/YsMdTEWwWRlZVF27ZtAc9z63A4zru+Dh06YLfb+fjjj83rTp48SVhYGFdeeSVLly41r09JSTFfV19hYWFkZWVV9iExYMAAli9fbg6xcTgcJCQkAJ5KiC5duvDAAw8wdOhQduzYcc51hYWF0a1bN/N5T0xMJCEhgR49epRYrmPHjuTn55vvs08++YTTp09X+jGUlyos/Mh3GIhLFRYiIiIiIjVq2rRpTJo0icWLF9OsWTPmzJlj3vbwww8zfvx4unfvzrBhw/jpp5+4+eabAfj9739Pu3btAE9lw86dO7FarQQEBDB37twSVRdeAwcO5IcffmDo0KGEh4fTo0ePMns4dO7cmQkTJvDQQw8BniqIGTNmmLdfccUVzJkzh4MHDxIdHc28efMAGDduHE888QTvv/8+HTp0oHfv3ud9/Hv27GHBggWApyJg7NixtGjRgqioKDp06MCQIUPo2LEjCxcuZPDgwdx66600b96c6667zjxQjYiIMBsvhoeH8+677/Liiy8SGxvLPffcU2J7KSkp/O1vf8PpdOJyufjVr35Fjx49OHr0qLlMdHQ006ZN4+GHHyYkJITrr7+egIAAQkJCzvt4ymv+/PlMnz6dIUOGYLVaGTJkCGPHjuXJJ5/k0UcfJTw8nGuvvZaIiIjzrstut7N48WJmzJjB4sWLsVgsPPTQQwwfPpz58+cze/Zsc9hOkyZNmDVrVqn3x9VXX80///lPhg4dSp8+fXj66acr9HiGDx9OZmam2a/D7XZzzz330LVrVxYsWMDBgwex2Ww0a9aMWbNmlev5mTJlCsuWLcNutzN37txSibzAwEDCBscMAAAgAElEQVT+/ve/m9P49u7du0Sfl5picddGa88akpyczMCBA9m4caOZGatPDh07ze/neebCHT2wE2NuvdjPEYmIiIiIVJ/du3fTrVs3f4fhN9nZ2YSFhWEYBk899RSxsbE89thj/g6rzvE+T+CpdHj//fd55513/ByVVMWZn/3KHrurwsKPfHtYuDQkRERERESkQZk4cSJHjhwhPz+fSy65hIcfftjfIdVJb731Fp988gkul4vw8PASU31K46aEhR8ZaropIiIiItJgvfzyy/4OoV4YN24c48aN83cYUgep6aYflaiwMJSwEBEREREREfFqEBUWe+Y/z8kmTfwdRoXl5hdyb7Kn8U7k+mB2fBvm54hERERERKqPcdtgcpKTa2X6QxGpG9xuN46TJ9nx1BTzutScnEqtSxUWfqRRICIiIiLSoJ06TWZODvW4z7+IVIDb7SYzJwdOVc+Upw2iwqLLhMfq5SwhOxLTWL54MwA3X3Uhg+7scZ57iIiIiIjUH4WFhSQnJ3MoK8vfoYhILQkODubiGwcScMtg87rk5GT4eF2F11VrCYtHH32U5ORkrFYroaGh/O1vf6Nbt24cOHCASZMmkZmZSUREBHPmzKF9+/a1FZZfGS7fHhaGHyMREREREal+AQEBdOjQwd9hiEg9VWsJizlz5tC0aVMAPvvsMyZPnswHH3zA1KlTuffeexk2bBhr1qxhypQpvPnmm7UVll+5fErjDDXdFBERERERETGVq4eFy+Vi4sSJOByOSm/Im6wAyM7OxmKxkJ6ezq5duxgyZAgAQ4YMYdeuXWRkZFR6O/WJoVlCRERERERERMpUrgoLm83G5s2bq9zd96mnnmLz5s243W5ee+01UlJSaNGiBTabzdxObGwsKSkpREZGVmlb9YESFiIiIiIiIiJlK/csIQ888ACLFi2isLCw0hubNWsWX3zxBY899hhz586t9HoaCt++FRoSIiIiIiIiIlKs3D0s/vWvf5GWlsYbb7xBZGRkiWqLL774okIbHT58OFOmTKFly5YcP34cl8uFzWbD5XKRmppKq1atKrS++sqbr7BZLUpYiIiIiIiIiPgod8Ji3rx5ld5ITk4Op0+fNhMRmzZtIjw8nKioKLp168a6desYNmwY69ato1u3bo1iOAgUV1gE2K0aEiIiIiIiIiLio9wJiz59+lR6I3l5efzpT38iLy8Pq9VKeHg4S5YswWKxMG3aNCZNmsTixYtp1qwZc+bMqfR26htvVUWA3aoKCxEREREREREf5U5YOBwOXn75ZdatW0dmZibfffcdX331FUlJSdx3333nvG90dDQrVqwo87a4uDhWrlxZsagbCJcSFiIiIiIiIiJlKnfTzWeffZa9e/cyf/58s39Fp06deOedd2osuIbOm6Sw220aEiIiIiIiIiLio9wVFp999hmffvopoaGhWK2ePEeLFi04fvx4jQXX0JkVFjZriRlDRERERERERBq7cldYBAQE4HK5SlyXkZFBREREtQfVWBhuDQkRERERERERKUu5ExaDBw9m4sSJHD58GIDU1FRmzJjBbbfdVmPBNXQuV3HCQkNCRERERERERIqVO2Hx2GOP0bZtW4YOHcrp06cZNGgQsbGx/P73v6/J+Bq0EhUWbiUsRERERERERLzK3cMiMDCQyZMnM3nyZDIyMmjevLnZfFMqx7fCIjff6edoREREREREROqOcldYzJ49m40bN3L69GkiIyOVrKgGxRUWNlVYiIiIiIiIiPgod4VFaGgob7zxBn/5y1+48MIL6dOnD71796Z3795ERkbWZIwNlndmkAC71ay2EBEREREREZEKJCz+9Kc/AeBwOPjxxx/573//y+TJk8nNzWX37t01FmBD5p3JNMCmHhYiIiIiIiIivsqdsMjJyeH777/n22+/5ZtvvuHo0aP079+fPn361GR8DZq3wsKuWUJERERERERESih3wqJPnz60adOG+++/n1mzZhEXF1eTcTUKhuHGZrVgtVgwXIa/wxERERERERGpM8qdsPjDH/7A9u3befXVV/niiy/M/hWXXXYZAQEBNRljg2UYbqxWCzabxRwSkpCUwY7ENLrHRdO1vXqDiIiIiIiISONU7oTFuHHjAHC5XOzcuZPPPvuMsWPHYhgGP/zwQ40F2JC5fCosXIabhKQMJr+yGZfLwG63MuuRfkpaiIiIiIiISKNU7oRFZmam2b9i27ZtHDhwgEsuuUQ9LKrArLCwWjAMNzsS0yh0eoaGOJ0GOxLTlLAQERERERGRRqncCYvrrruOyy67jF69evHkk0/Ss2dPgoODazK2Bs+ssLB5Kiy6x0Wbt9nt1hKXRURERERERBqTcicsvv32WwIDA2sylkbHW2FhtXgqLLq2j8Rus9I6pgl/HN1D1RUiIiIiIiLSaJU7YREYGMi2bdv48MMPSU1NJTY2lmHDhtG3b9+ajK9B81ZY2KyeCgu3241hGMQ2D1WyQkRERERERBo1a3kXXLlyJX/+85+JiYnhpptuIjY2lscff5wVK1bUZHwNmmG4sVosWIt6WLgMN4YbCp0uf4cmIiIiIiIi4lflrrB47bXXeOONN+jatat53S233ML48eO58847z3nfkydP8sQTT3Do0CECAwO58MILmTFjBpGRkfz4449MmTKFgoIC2rRpw7x584iKiqr8I6pHXIaB1WbFZrXiMtw4Cj2JCkeh4efIRERERERERPyr3BUWmZmZxMXFlbiuY8eOnDp16rz3tVgs/O53v2PDhg2sXbuWdu3aMX/+fAzD4K9//StTpkxhw4YN9OrVi/nz51f8UdRThgG2ogoLgAKHJ2GhCgsRERERERFp7MqdsOjZsyezZ88mLy8PgNzcXObOnUvPnj3Pe9+IiAiuuuoq83KPHj04evQo8fHxBAUF0atXLwDuvvtuPvnkk4o+hnrLZRjmtKYA+UUJiwJVWIiIiIiIiEgjV+6ExYwZM9i7dy+9evXimmuuoXfv3iQkJDB9+vQKbdAwDN555x0GDBhASkoKrVu3Nm+LjIzEMAwyMzMrtM76ynAXzRJiJiycgCosRERERERERM7bwyIvL49XXnmFvXv30qdPH+bPn2/OEtKyZcsKb/CZZ54hNDSU++67j//85z+VCrqhcLmKZwmB4iEh6mEhIiIiIiIijd15ExYzZswgPj6ea6+9lk8//ZRTp07xt7/9rVIbmzNnDgcPHmTJkiVYrVZatWrF0aNHzdszMjKwWq1ERERUav31zZkVFnkFqrAQERERERERgXIMCfnf//7H66+/zhNPPME//vEPPv/880pt6O9//zvx8fG8/PLLBAYGAnDppZeSn5/P9u3bAXj33XcZPHhwpdZfH7mMkhUW3h4WDqcqLERERERERKRxO2+FRW5uLrGxsQC0atWK7OzsCm/kl19+4dVXX6V9+/bcfffdALRt25aXX36ZuXPnMnXq1BLTmjYWhuEu0XSzoKiHhaPQhdvtxmKx+DM8EREREREREb85b8LC5XKxdetW3G43AE6ns8RlgKuvvvqc6+jUqRN79uwp87YrrriCtWvXViTmBsMoqrCwnlFh4XaD0+UmwK6EhYiIiIiIiDRO501YREVFMXnyZPNyREREicsWi4WNGzfWTHQNnGdIiNVnSIjTvK3Q6SLAXu5JXEREREREREQalPMmLDZt2lQbcTRKhuHGHuDbdLO42aaj0CA02F+RiYiIiIiIiPiXTuH7kVFUYWG1el6GAp8KC4dmChEREREREZFGTAkLP3IZhqfppqVkDwvwNN4UERERERERaayUsPAjw8DTdNNWVg8LTW0qIiIiIiIijZcSFn7krbCweissClRhISIiIiIiIgJKWPiV4XZ7hoSUUWHhUIWFiIiIiIiINGJKWPiRy+XGZvGpsPDpYVFYqISFiIiIiIiINF5KWPiR4XZjtVmwWUtXWBRoSIiIiIiIiIg0YkpY+JHLcGO1+AwJ8elhUahpTUVERERERKQRU8LCjwzD7ZklxFJcYVH0Lw4NCREREREREZFGTAkLP3IZRU03rZ6XId/hIjTIDqjCQkRERERERBo3JSz8yKywKHoVChxOmoQEAJolRKpPfGIa7322l4SkDH+HIiIiIiIiUm52fwfQmJ1ZYeF0uT0Ji5N5ONR0U6pBQlIGkxdvxg2sCLAy65F+dG0f6e+wREREREREzksVFn7kqbCwYi2aJQQgNLiowkI9LPwmISmDlRsbRkXCjsQ03EX/O50GOxLT/BqPiIiIiIhIeanCwo+KKyyKExaBdisBdqt6WPhJQlIGk1/ZTKHTILABVCRc3KE4drvdSve4aD9GIyIiIiIiUn6qsPCj4h4WPgmLABuBATb1sPCTHYlpOIue+4ZQkXBhy2YAtIgMrffJFxERERERaVyUsPAjwzBKVVgE2K0E2q3qYeEn3eOizQSSzVb/KxJyC5wANAkJULJCRERERETqlVpJWMyZM4cBAwbQpUsX9u7da15/4MAB7rrrLgYNGsRdd91FUlJSbYRTJ7jdbgw3ZVZYBATYlLDwk67tI7nmstYA/G7YpfX+ID8v31nir4iIiIiISH1RKwmLgQMH8vbbb9OmTZsS10+dOpV7772XDRs2cO+99zJlypTaCKdOMAxPK0TrGQkLs8JCQ0L8JjjQBkBMRIifI6m63KJERW5BoZ8jERERERERqZhaSVj06tWLVq1albguPT2dXbt2MWTIEACGDBnCrl27yMio/zMzlIfh9iQsbGUOCbFRqFlC/MZ7kJ/TAKoSvImK3AbwWEREREREpHHxWw+LlJQUWrRogc3mOZtts9mIjY0lJSXFXyHVKperqMLCcsaQELuNwAArDs0S4je5+YUl/tZneUU9LAqdhmaeERERERGRekVNN/3ErLCwWbBZi1+GgAArgQE2CjUkxG/MCou8+p+w8K2sUJWFiIiIiIjUJ35LWLRq1Yrjx4/jcnnO+rpcLlJTU0sNHWmoXMbZKywC7FYK1HTTb3LyG84wCt/H4K22EBEREWlMEpIyWLlxLwlJjWPouUhD4reERVRUFN26dWPdunUArFu3jm7duhEZWb9nZSgvb9PNMntYBNgoVMLCb4p7WNT/CgvfJEVDSMCIiIiIVERCUgZPvbKZtz7ezVNLNitpIVLP2GtjIzNnzuTTTz8lLS2N3/zmN0RERPDRRx8xbdo0Jk2axOLFi2nWrBlz5sypjXD8KiEpgx2JaVzYshlQepaQQLuVAM0S4ldmD4u8+n+A79uHoyH05BARERGpiB2JaeZ+tdNpsCMxrd5PWy/SmNRKwuLpp5/m6aefLnV9XFwcK1eurI0Q6oSEpAyeWrKZQqeB3eYpbrFarSUrLAJsRbOEqMLCH1yGm3yH57lvcBUWGhIiIiINnPfEUPe4aB2UCgDd46KxWiwYbjd2m5XucdH+DklEKqBWEhbisSMxjUKngdsNTpcn02uzevpYeAXarUWzhKjCwh/yGlhFQm6+E6sFDLeGhIiISMOWkJTB5Fc243QZBNitzHqkn5IWQtf2kXS+IIKEgycZf1dPvSdE6hnNElKLusdFF1dWFCUprFYrVqsFb84iIMBGYIANhyos/KKhzaqRV+CkebNgz/8NIAEjIiJyNiVODBWV/kvF/bg3lXc/3dOgej14Z+eLCg/2cyQiUlFKWNSiru0jmfxgbwCuuqQlgNm/wpvACFAPC7/yDgMJDrQ1iCEhufmFRIeHFP1f/xMwIiIiZ+Nb6m+3q/S/MhKSMpiy9Gve3pDQoBpUnsp2eP7mOPwciYhUlBIWtaxnlxZYLeAuuuztX+H9G1g0S4hhuHG5lLSobd6D+uiIEHLzGkLCwklE0yCsVot6WEit2rU/nRWfaQo5Eak9nS5ojsXiOemg4SCV8/O+NIqKERpUlcrpnALP3+wCP0ciIhWlHha1zGa10CwsiIxT+UBxhYXNZgEnBNhtBNo9eSSH0yDEppxSbfL2rYiJCOHIiWwMw11iFpf6Jq/ASWiwndAge4PoySH1Q0JSBk8u3ozhdhMYoHHkIo3Nz/vSSEjK4LKLarfxZVaOA7cbHIUuOl3QvNa225C0b93M/L+hVKk4Cl3kFXiGWqvCQrzi96exa3/tf0+djRoGn50SFn4QERZE+mlPwsJ2xpCQwABPhQV4vmBDgvQS1SZvhUVM81Dcbsh3OAkNDvBzVJWXm+8kJMhOaLBdQ0Kk1uxITDPHC2sKOZHG5adfTvD0ki1YgIBaTliezPLsWxluzxn15k0bZr+Cmjywad40yPz/qQf7NIjv7tM+SYpTqrAQ4Ic9qUxZ+jUWC3WiQa/vTJJ1IZ66Rqfv/SA8LJDMrJIVFlar56UItNsIsHsTFo1nSEhCUgYrN/q/fNxbhRAd4en7kJNXvw/y8woKCQ0OIDQ4oMQUp1J/1ZXPyrlc2jHK/L+hnKGT0urDe1Fq37adxwDP0NeKDimo6nvqZFbxwWhmVsM8ME1IymDSy1/x1vrdVe4xsftAOss3JJRYR1pmvvl/WGhglWKtK3yTFKezVWEhsDU+BaDONOhVw+Bz0+l7P4gIC8bp8px9PLOHRUDRtKYAhc7GMVNIXcoq5ngrLCK8jSoLgZAKr6culHUVOl04XW5Cg+2EaEhIgxC/P43JL28Gav/MZUW0jgkz/x9/p6aQa4jM7+1Co06/F6X2RfvMwlCRhGV17Av4JilOni6gQ+sK3b1e2JGYhsuoegWbZwrYLThdBqs2/cKscZ7nOy0zz1zm6IlsOjeAoTXeYSBWq4VTOQ0zkSUVExVevG9fF06sdI+Lxmqx4HJ7hqL7O566RhUWfhDetDhjXTy9qTdhYSPQW2HhNPh6R0qDP4O1IzGNwsK6kVXMzS/EarUQWTQVaGVmCqnOsx9V4R0CoiEhDce2+GO4qdyZy9p04mTxDm9osPLiDZH5vQ0U1uH3otS+oKJhrQF2KzP/3zUA5dqPqY59AW/1KkBmdv45lqy/Lu5QPRVsOxLTcBY1dy90FT/faZl52G0WLBY4mpZT9YDrAO+QkFZRoSWGh0jj5T053DIytE4k3Lu2j+Tq7q0AuOHKdn6Pp67RnqQfRIQVjw+02UomLAIDrAQUfYg2bT/EB18k1pnxVTWle1w0VqsFl1E9WcWqVDfk5jsJDbITGmI3L1fUz/uq5+xHVXljDw22ExocwLH03FqPQaqXd6gSVH5HtTaqf1JPFr/Xqvq+qwvVSlJa97hoLBYLbrcbi0Vng2pKfXz/pxYlLAudBqdyHMx9azuFTuO8DXirY1/gZFYBdpsFp8vNydMN80x6jM/vwN8euqrS7wvPZ9hTEm/zeb7TTuURHRGCYbg5eqKBJCyKhoS0a9GUvYcy/RyN1AXe7ymL1VJnvlu9CcSCwsZRYV8RSlj4gW/Cwlth4TskJKiowmJbfNE4UHfDblzXtX0kF3eIYkdiGtd0b12lx+itbjDc7koleXLzCwkNCaBJUaPNnEpMbdq+Vc132C7PTqx3CEhIUAChwXbyCjQkpL4LKJpBqGloAFN+27fCnxVPCfDmch08VMWJopJim9XCsYzK7/D6Y9hBfTxA9Ieu7SNp3tTTQLp5WBBdLqzbZeP18XX1Z/PKqvBNWG6NT6HQWXQW/zz7MV3bR9L9oih+3JtmnmGs6Ot28nQB0REhZGYVlOhnURU796fx8740enaOrRPPv+93akQVmop2bR9JRFgQJ7MK+FXPNuZjS8v0JCzsVitH07KrHG9dcCrHgdVqoXV0GNt3HzcTrWeqj98TUjmpGZ7vqbTMvDozI6A3iXL0RMP43FUnJSz8INynA7OtqNlmcYWFzayw8P3Rrwvjq2pSdp6nRK+qYwu/3XWsStUNZoVFsLfCouIH+U1CimcVqcxB5fnsPpDOxJe/Ave5d2K9TTaLe1hoSEh9561WyM13clG7iArf39vUCWo2CZp6MpfgQBstIkM5XoUKi7KGHdTkTqQ34ekyNB3r+ZzMyif9dD6topuQkpbDkRPZtI1t6u+wymROs2sY2OtRteKWn48CJYeA1Ye4T5zMo2ObcPYfOUVwYPFuZnmqJrKLGl3nFjgr9Xk8mZVP86bBWLCYM4ZURUJSBk+9sgWX4eb9Tb/UifeOb9VaakZuiZMkFZFf4CSzqPKg0Ok2r0/LzOPiDlGEBNv53w9Hqhasn22NT+Hw8SwOHjtNs9BAIpoG4XS5yc13lthXgzMS+nar2dPjzGWU0GgYTpSoBKsbMwqdKDruO5qWc9akWmOlHhZ+UKLCougVsFk94wVtVovZw8LbmDMkyF4nfiRritvtJqVonOSBo6dxu93nucfZNW1S3B+kMkke74+YWWFRiYP8Iz6Z0WZNqr/D9tfxKbjd5+9jkFvg28MigHyHy0zmSP10LN3zOXEZ7hJ9IsrL9/1otdZcEvTEyTximofQMqqJGXNldI+LxuLTmLimk7ZlNbPz0owYJSUknQTg3kFdAdi++3gl11Pzz+u2nSk4XQZGHeiTVBFhPgdU9emkRerJXOLahNO8aRD7kovL73t1bXHO/Ri3201K0e9n8vGsc34ez+ZkVgERTYOIaBpULbOEVCaGmub7nVqVCrbkE9m43Z5EUkrROg3DTfqpfKIjQmgdHUZ2XmG97fmwafshZr3xDW99vJtvdx4jKNBGeJjnN7Csk2O+CX3fnh5e3sTnWx/7tz+ZVI/Uk7lEFTUI9m006y+5+YVk5RYSFR5Mbn5xMlE8lLDwgxI9LHwqLALsNiwWi1lhAXBJxyjyCpy0iAwttZ6GsgOdmVVAvsNFm5gmnM5xkHG68mdF8h2ecV9WC0x/+OoKJ3ly8gsJCbITFGjDarVUqsLCt5SrJsZ/ehuagWeu+dSM3DLfAyV7WHjOcuVpppB67Vh6jnkQk1LBZmgJSRms++oAdpvn+2Vgr+pt6rT7QDrv/mcPCUkZpJ7MJaZ5qCdhkZFb6SRk1/aRtG/lOWvf5+KWNZ607dg63Pzfbis+QPzplxNMevkr/tVAdlR/+uUEb2/YXaXHsedgBnablWu6t6JFZCgffXWgwuurzgOAc/0eBtiKf1Pr04F/QdHU5t7mlZWdCaI29xMchS5OZhUQGxnKBS2bsuegJ7HVMiq0xNSSZcnKLSQn30lIkI0jJ7K5oEVxxY7NVr7XLTOrgOZNg2jeLKhaKix8p2j2Z+d+39fxeHouraKaEBhgIzXj/AdaZ3sPHD6eBcDFHSI5VvR7kpldgMtwexIWMU0A+NfHFf+uqAv7p5t/KqpQcnv2lQCaNfHsf5c1tWnLqCbm/9Yy+vJsTziO01U3GsTXR7v2p7Pis5p/T5TnvZebX0h2XiGXFDWwrcwJoIpu83y8MfToHAPUzPFDfaYhIX7QLMz3LGfx2cPAovHp3goLgGG/6sjO/ensS86k98Utzeur2qvBu47KlrZVZ1mctwv1NZe1ZuXGXzhw9HSJ6YYq4peiZkqGm1LlfuWRl++kSWwAFouFJsH2SvWwOJqWQ3R4MGmn8ktUW1SXUzkOAgNsxLVpxu6kk2zYdpBN3x0u9R4oHhISQGhQcRPRhjKvennt2JfGdwnH6Xtpq3K/V6vj/X2udezcn8aOxHR6dIop9/rdbjfH0nO5okssm38+SkpaNhBb7lieemUzDqdRNAtOENnVmLzyltI6XW5WfrYXu91Kp3bNaRkVSoHD5TmIaFbxcku3223ukJ+qhbN8dp8D27tv7mK+Np98nVQnGulWh4SkDP726hbcblj9+b5KV+99l3Cc5k0D2bT9MGmZebgMN08t2Vyh9f28r3iWgrKe1/J+Dj/ecoBXVv+M202ZpdxppzwHroFVOPD3h4MppwFPyXJ0ROmTFufjSQh9hdNVe0OcvGcqY5uHcGHLZvz0SxrNmgQS1zqc7/eeOGeZc0pRv4QrurRg889HOVR0QA1wx4BO54290GmQleugebNgLBYLP/9y9gPK8r63fL+3fPs81JaEpAw2fnuIT785hFE0NKZFZCgto0Kx2y0lhg6f7f5PLv4Kl8tdagjp4eNZ2KwWrugSS3xiOtm5DvP1i4kIIa/opMcnXyexcfuhcr9/amOq+vK8fr4FpRYgtnlocYVFGcmzL39MxoJnaHbHNuGl1uvdj4L6lfisDrsPpLM9IZXe3c5dJXWu+09a/JXnO7qS30Xlec2/232c6a9thfNMVuBNDlzcMYovfzxSpQoL7/5Voatq73fvZ7lHpxg2fnuYoyeyucQnYerdVmMdkqSEhR8EB9oJCbKRV+Aym23arBZzih3v34imQYQGB2CxwL7kUyUSFt8lHK/SDrTvWL0Au5Vnyxird777Oqv44fTy7qRc092bsDhFr24tyn1/7wf40o5R/HL4JBe1i2Df4UwOHD1NB58zpuWRk19oViOEBgdUqu9DSloOHdqE46bk8JDq+qLZfSCDiztEcmlcNLuTTp61Kau3miK0aEgIFA8TaSwSkjJ4+tUtGIabtV/tL9d71XvW1+Uq3eixvK/h93tSmbb061I/mrsPpPP2hgR+KtqRXrlxb7k/P6dzHOQVOOnWIZLtCcc5WoGhFr6lrrjdNAkJIOno6XLfvzzr9w5hK3QaOJwGsUVDQsAz5royCYv0U/lk5xUSGGAj6eip847pPN/rc77bE4+cAjwVWnk+n5UCR3HH7ro8O0t5/LD3BN6Cl8omX+IT00hK8RxMLln9M4ZR/NpXZH0to4oPwm0+Q5QSkjL49/8S+eqnoyWSEECp5/CnX06weNXP5nrKiiE+MR2b1YLDaVTqfegvB4+dpmVUKMfSczl0/DQxzcufyE9IymDpBztKfC7P9tqU971ZnuWOFzWyi2keyqHjnt+/0zkOtu06htPlJvVkHmmZecTvT+fyi0qux1s11ucST8LiP9sOEWi34qa4evJcsXgPQps3DcICZOcVUuh0EeBzEsh73/ImcvYXfScE2q1klnFWHuDHvan8cjizRDFZjPsAACAASURBVCxnxleZz7/3wN9RVGkDntfxeEYel3aMxmazms/32Xy9I8V8D5z5eT90LIvWMWG0jQ0DICU9xzxwiwoP5tuiYV4V7aFS072HypsQ8T4Wu81CgM1K2xZhhBdVWJyZAP8u4Thbd3ga3Rc6jTKbHvo2cf3d0Ev99j1e2ycbfU9IrNq0l5uvas+AClZo/q/ouxxq9pjls28P4QY4z2QF3uRAXNtwAgNsZqPwytiRmIajGnqDeRtuXtIxGrvNUuqE5+4D6Ux6eTNud+nkY2OghIWfRIQFk1eQY1ZYeIeEACQVnVXJzCpgxutbiYkIITG55DRM3qEkULJ0ubxKjNWr4Aesuhr3eb84k1M9Wf4OrZvRvGkgX/5wpNxfpiV+uGxWHE6DkTd04mDKaQ4cPQW0K3ObZ1t/br7TTFg0KUfCwru+pqGBZOU6uKRDFEdPZNOjcwwFDpf5o1ddsx3k5BVy8Nhprunehcsvima5xXMWoayDqNx8J1YLBAXaCKlCE9H66Od9aew5mEFqRm6FD6bOdtZ394F0Jr+ypVTjvjPfA93jovl068FSP5oAkxZvNuM5c/3n4x233Cq6Ca2imnAsrfzNLLvHRXtOMRW9V7pcEMmm7YcoKHSVGGJ0NmU9Rt+Yu/n87310MREh5gHpv//nmZ75fI8zPjGNnfvTCQ8LIivXYVY89L20JV/+cIS0zHzzoK2sA4Jz7dCUJ9G6/8gposODadYkyKzWAkhO9XyOQ4PsTHu4crOz+KuZZ/z+NHYfyDCfp7Dg4sqzypa4r/58n/m/t7u6y6j49KbZRRVsVquF7hdFlXodvRxOg43fHmLj9sOlXr93P91Tar2+MZw4mUdKeg4DerVj0/bD/HI4s0Tpd111Ktszy8XogZ1YufEXDh3L4squ5Uvk+1ZUeZVV4m4uW46DP+/r4nKdu3Gpd6c7tnkoDp+p+bwnWDb/dIT/+2gXhhtWnPFZOJqWg8WC+ThT0nO4pGMUBQ4nB4oSB+BppDjrjW9KTfnuHQLi2zwvM8tRKtHj+Y73xOMoNFi+IYF7B3Ut8/HsP3IKq9VC30tb8dO+0hUin207yIsrfiwxkwvg+a5xGthsFq7o0oLv9xzHMNwVavrqPfD3ZcGCo9BFy6hQLBbOWYKekJTBlh1HzctnDqs5fDyL9q2bFSeW03LN5zA6IoQenWJ459M9nrgrsJ/pO11qTUx5vCMxzUziOAoN3vp4N/ff0q3Ec5rvcHLIJ+HndLlo1iTQrHA+s8Ji/ZYk83+3282pHAensgsI9xnCHZ+YRsfW4Rw8dtp8n/uqjaS097PtPM/n8Ew796fzfx/tZHfSyQrPOuR7QsJlwMdfn73i5mzPgcMn4Vje4V1nxuB7zFLWaw4lTzTY7VaahgaycuPeUvH4fk/FRIRUKWFxoU/T28o8Nq8TJ3Ox26xEhQfTKrqJWX3u9eWPRzDcDaPSszKUsPCT8LBAUtJzfCosrOaUhYnJxT/M3h+8n/elkZCUYb45j5zIxmb1fHmMuOGiCr9pLz2jzKhD6+IP3Pm+dH2/wMv6MXK73Xz101FS0nK4zOcMyo59acTv90wN5na7zYM/gMhmwfxyOJNT2Q5OZjnM0mIofUbNl+8Pl/fLrOuFzbmwZdNSZ5DNAxZn2UkDR6ELp8swh5KEhtjJyS9k94F0diT+//bOPEqq6tr/33tr6Krq7uru6m56YOgBhEYZJDYCooKgIlHEgIiJotEnxmiiK2/pT/IziUuieT9dmsTHQgz8jOuneb74CC8SFfEFUBRRZmmmZuh5Hqq6uubp1v39ceveuvfWrepqpm5wf9ZySXffYd9z9tlnn33O2ceu+Bb588T3MhAMZDgSRWlBJgIhTtpDmU4Hmw4nm4QVFRMrbKgqt+EHN43Dph1n8OQ90wBAYZh9wQjMGXowDCM79WR4rLDQ0rGzjfyr79m6uxFrNx0Gw0BxTBUDwZEa6D2ZJu2s9lu/aUoIZABQzICJjsDIwviASOw018lmorWen+qbAKAjlhm+2GZBSUEmWrvdSa9V/744PxPgganjCnD/womwuwLYtq8ZLZ1u6bSRVM+SD360nB2xEx2RZ5YcgcI8C/pcglO463A79h7vTOkgHau343+/8ZUU8GAQr7853xuFLw61obGjH4V5Zs3gw0BBWPnxisk6+/p2JypH5iLPmhGb3efR6wygw+7F6KIstHR5NGfoB9Kpgye7pQGbWrYjZ3pR2+SQ9OB8Orxbdjdg3aYaxcDO4Q6AjbWN6VcWJdT1jv3NAJikM2j/3NOE/Se6wDBxm7dy8WRs2nEarOo8+4ECXYdP96Ag14xrqkZgx/4WvPdpLZzuoCJYIfLFobaE+utx+nGs3g4mFowDAzCM0oEU2+ntsyvwxaE2nGlx4oarR55bweLCD06aOoX+a1JlAbbtbUZzp3uAO+LIZ/wYCNsjsyyGpH1oOrPhO/a3xNtXOIod+5uTLLX2gWWEGfobrx6JT79ulAZXoXAU/7OnWVqqr26HHb1eFOaakZOVgcI8M3r6/LiywganO4g9xzqlYME/vqgHkHjkuzgDnis7ia3PHZBshlhfebK/A8C3p3pwrMGu6XPUtfVj1IgsTKyw4Ytv2+BwBRRbVrd80yjIIvueUISL6yrHY+/xTun6VAONw6d7cKLRIW0VlPtpeh2D/ByzNDNcFAsyePxheP3hlCdeiMGDm6fH23Q4wqGj14v8XLOUVLPD7kVThws6lkF7jwcTK/Lx7Ipq/J939uHKSpvUlgbS97GjcmEw6BAMcbBmGgZ95PFAdmPsKOUJWTVnehO2ozW0uRDlge9fV4E/f3gMAJCTmQGTUchR5vKGFDrR2esFA8SS37MIc1HUt/Vj2oQRUjnXt/Vj+S0TkGUxYNe3bbCY9AkBc3mfBCDld5wNu2vaBz2bf+hkN57f8LW0woGH4Icma8MiYvnodYmrGtXvPlrfi092N2J3TbtmYK6uvR85WUb0e0K4/bryQZfD+DFKHdKqc57n0dDeD52OAcfxuGf+eKz/4IjmMe5icCA3KwOFuWb0ppHDIpnNd/THc+XMvSb5ypNj9XYcqO1Kmo+ru8+PwlwzWJZBtsWI4w12xbhP7r+fzUT1pQ7Dn8uRDOeJhoYGrFq1Ck6nE7m5uXj55ZdRXl4+4H2tra2YP38+lr/6f5FdmP4WguFAU6cLbm8I48fkwWjQ4UyLExEuijHFQpIp8bQMhmFiRkb4d0WpFeYMPWob+5CdKRjBLLMBo4uy4fWH4QsIHZfFpOy8fIF4p2YxGeAPRFDX5oxtewgj02xEkc0MnkdsZYIQjMiLZduWthQEwmjr8SIU5mCIDc4Lc83IsgjP9frDaO/1SEuoRZnVz7Wo8kMYDTrkZWcoljfarCb0uYNSOVSUWhO+q8+VmCeiojQHTncQLl9IMfPb7fAp9nsW2Swwmwzwx8rMaNChttGBkoJM5OeYUd/Wj0AoIg0y1TL09PmSLscsL8lBIBRBp92LieU29Dr9CRHcZN+UitZuD5zuAMpLrMiyGBEKczjV3If8HDMcroCirByuILz+MCaU5SEYiuB0ixM5WRnIzzENqB+p/ib+rNOx4GIBHrl+JHuO/Hly/RaDZerfpSoXXyCMPncwNiBW3nOyqQ/hiHL5cIZRj1CYEwJZHe6Ee5Rl7Ea/JwSGEe4bO1LYVlTX6pSi9+K9Hn9YOstbibDFKxTmkJdtgtMTVCSeFGfozBl6VI6Mb1tKVjZefxihSBR9rgCurMxHt8OHXmcg5ogqv8XnD6O+3aX4xmCIQ1uPB+NG5cKUoUcwzOF0cx9GFmYhz2rSfK9YLt19Ps1vLLJZUJgnrKDotHth7w9g7KgcnGkRViaMKbYiGIoo2kiW2YgRNrP0TXI9aelya+4rZlkGE8qEFS7iO9t7PQonochmgdGgk5LIAYIdkDvx8ucryiums2aTAY3t/SjMs8CgZ9He48H4MXnwBcJo7fZgdFE2WrrcUpmJ9ZVMD+V0ypZaA8K++LzsDPgCEWnljNzWg2Fgy86AKUOvaGPy9qUuQ622d6qpDyFZWyiyWeDyhsAwDFhWmKkVHcF+T1BRflrf4vGFpBWA6v5B1JOqchuCIQ72/kBsICR3MRiUFmbCZjUBPHCi0YHsTCMsGXq0x7YGSqOr2DuyLAa4Vcu3GYZBSX6mdI8oS4ZRh45eL0YXZUuB9cYOF7z+MCpKrOiwC6saxa2CyezVQHZM3V7U/eT5wN7vR0evFxPKbGjtdiMa5RMGasnkdboD0qoghmGQZ82Aoz+ACWU2aWJExOsPS30zAKlvUdPQ3q/os9X2SXx3a7cbXn9EGqTKZWvudEtB39hTUDkyrmN1rf1gWaHt1rX1wx8Iozg/EwwDqSz0Oga1jX3gYpMdDMOgpCATHBcFF+XR6/RjQlkeIhEedW1OlBVbodMxivqyZgoDJ4tJr1h1mGHUIRiKQt6Wm7vcyDQbkG81ob6tH2XFVuk0sijPo7bBIQVsAWFlgtsXRjCkPTkgl1euX+q2JfebrJkZKMg1IRSOorVbaKNjR+UiFObQ0uXGuFG5iPK8dCpKbnZGQt+k17FSbgZAGGTJ2w8b88vcvpBCBovJoKh7ufxaPgAgDO4b2/sFH9UbwthRuTBnJJ8fVds19WlxalvU1uNBnyuADKNeUc6iH2sxGRTt50yLE1w0KtmFk019yDDoYiu8RD+bhy3HBIOOhSlDj6YOF4psmdLqHLc3hKZOYZtxvzcER79fIZu6vOX+q6z24/YvDdS+Fssy6Hb4Jd0HhNVEeVbB1mnZLK8/jOZOt+IeOTarSdN2CXYh7kfwPA9rprDyUfimeNtV9x3yZxv0LDKMejR3ulBks8DhCsKUoUNZceqjeMVvZ1kG0SgPno8dl56hR0C2ikL0KYST8CI40+JEQa4ZvU6/pH8icp+lpcsNfzCC8WPy0NbtgdsXwpjibEV5MwwDLhqFQa+DPxiBM8l4pLHdhVCEg17HIhrlMW50bkLd6VhW0d60/ATR/hXZLKhvU/pwFpMBp5r7wPM8wpEo8nPMyMkyDuhvD0fcPV14/+lHsH37dowaNSrt+4bFCovnn38eP/rRj7B48WJs3rwZv/nNb/DOO+8MtVgXFHGpM8MIDTMQM7pC3gUrKkqt8PrDCEei0qkZPM+jyyE0WC4aRYaBBW8xwOMNKwbj6g7FFwhLHZmo/IIjyaAwz4ymjjC8/hDq28PQy2akeZ6HwxVAnzuYMKgEhE6/1+kXHFUnA5NRpzAk4jOc7qBiaSrP8wnJLBmGgU7HSoYREDKGi/8Wk+9ZszjpHbnZQuJAhmFgyRBWQwCCg2qzmsBxUXTavbBmGmExGRBRGewoDzTKgijikY/eQAQME4hFM+OdjVj+RTYLLCaD4pvUZBhZ8Lyw1L6l2wOPL3Hvq/A8P3KyOEXnr9Oxim8UByP2/oA04GrqdEtGzGjQyTqRePmGIxyiPA9fICwt5+v3BOHyhpBlMYBlGJgz9PAHI+j3CAMLtSF1uAKxTMXCICrbbIDbF1aUi+iwsyyDXmcA8gGXXH6vPwwdy6LH6VPI2uP0g+P4BPnVQRCxXMJcFB6Zbsj1w2gIagYrSgsy0dDeH0viFn+PWP5ieRsNOjjdIWRnGpBh0KHX6UeXwwuAgT8YgcmoRyAUETrNEJfiuDcepQVZaOl2w+UNKWQVO1inRxjoRqM8AqEIvH4hQ768HDp6vfAHOUlmhmEQCEZgNOgA8PEZE1mbau3xJHxjJObsiI61MXYiUZ87CB6C86p+rzlDj9zsDPByf0tGOBKFLyDUk8sTgkHPKmYAWrrcKCnIVLRpjz8ET1u8zOSDPX+S/Comox662Ja5QIiDxx9OOK4wGOak7QWivnFRPu70MAxcsWS1oXAU2Ra54xsPNgBCEElMgNzp8MVmKBnodSx0OhbeQAR5VkEvRadCqnVZPThjZWvO0MPpCUKvY5Fh1MHrD6MvZleh0mHZD4rTkhiGQabZAI+i7TGS3CX5meiwx89tryi1IhDkFMEK8Tn+YASFeRboWAadvhA67V4Y9Il74eXf4vGHoWdZRcCXjyV8FttpltmAbghbMOz9yWareLT3eBAIRoRZsKgwG6lwpmPlkJttgs2aAa8/rHA62ViwRR4AEmWxWc3ocvjQ3SfYFJcvJNnehg4XsiwGeH1h9PT5oNOx6JCddS/qIaAMnmZZDNIgRrTPfS5lexH7SfW1yQIhclsvXhflhUGr0aADx0Xh9gl9WzjCwWTUo88VgM8fhjegDBZHo3zCQNfeHwDDMMjJyoDNmgEdy0gD1GyLUfHuQGzbhtgXdDl8CIQ48DwvXZNh1MPrjygGiGLfLp9UEINw4rdaTAbpP6EfitezQa9DNMrDnBHvHwKhCExGHRyugJTwscshnIgBAIFgBP5gBFw0ipysDPR7grHgYqyPihEKR2M2EnC4A8Ix4LL6cnmFYEVxvkUxOJbnqhHrNRKJQscyMBkFd7m3348wFwXHRREKRxHleUkWQHk8YnamMd5XxQJxpgy9ok8VA5M9qrYlDkwZhoEtNsmQYYiXH8dFpQTtbl9I0X4dscG8iKjHTncIXJQHyyChLTMsE+vb478T+2HRbom/V+d3UAft3LGgaGl+Jk56BRsj+k3qNhCJ8rDH/AaGYWDNMiacKiWWxwibsLWoL2YbQ2FO0b94/SE0tIdRUWqFPxiBXiesXDYadfAHooiIq44YoS+SPx8Q+mbRD9TrWDg9QbCsoE9C/xTbwi1bcCCWk3r1pNwnk12N9h4P/MEIzEY9ItEo9Bq2wGIyKAJKiTCxPHgR9LkDsQkRSGUo+nDqYKQWch+f54U8bnFbH28zLMtgTHG2FMDnOB7mDGHiMVlyeUesbxUx6HUxPQyiO7bCQfS95DaTi/JSwmE5LMPAlm1CR8ir8Cm8sToXbU9+jgn9nlDCiuJQOO6ziDlxfIEwDHoWES6aEChLhrxteHwhePwh5GRlwJyhR6fdi+Yud+wUGu1nJfP9Q2EO2ZnGWFtT+sOhMIdQmJMCol5/OGGS8lIKWpwNQx6wsNvtOH78ON5++20AwB133IHf/va3cDgcsNnSWzL06LSKQUVphgN//M+DqDvYg/lTK9Dc5caBQ0IiNJYBbi/Kx7L54wEkJl2S71b0GFgsnTsO/7nzFHpUz+9jlBmS5RT2htHV64XFbEBxhgXfHoonYUuGpdUnzJrFlk2xDDDmikI4TqnfnIgzZtx5KPxz6HQMMvQsfEEODAC7QVha/PXRdhysTXyu1k5NBsDsqaWoLMrBu1tOgI/JNnZiMU4f74EDgpG7fmoJmo92oqggE0V5Fuw/0QVnLHIrYpe9Rxfbk60lQzsLjBudh7rYcW06lsFdc8aiod2Fgye7hZ9nTYDdFcArW+oUcovL3MVnOwC0pCg7HSsshTvZ5FTM4kh6MnM8/q12L3Yf6VDclxlbvg4A/YYe3HTNaHx7MF6mdiSnsDeMsrH5ON3sxNEzyizryXbLaj3PIZP/VLNTszyTPXPcVUUoYTKw73iXlPRsIOTPEZcwRjgeC2eVYU5ZMVZ9dEbzHq3yd+tZLLlpHN4/2JzwXI7j4QDQFvsdA+C2WeWoHJmDLbsb0NDugtGgw4yKYnx7oC6+QgfCVopfxpYlHj7Vg1/9aTda2/1o7vQo6jcV7ppeLJxVjqMHlW3EXJiJk3ZfguMkl99l6JWWqu7adAI8D7SmeBfLCEv+87JMuHVmGayZRvzjyzq0dHnQB0GXSwsz0dnlAQOgR5wNid17+8KJmHzdKLyz5TiO1GlrnR2QtrfNu2YURhVlI8tswH/+z0n0uYMoHZuPRdcUoGdfJ+rautEYW76uYxmMH5OLE419CXZRzzKoP9kfcxrjfzPoWUy9ogCHDrejOcpr2sifPFiOnCwjVv3jtOK5bkMvRo3IQqfdjmkzy7Ftf3Ms0KekqDwPB5v6Ep6tYxnccu0YbD3YnnBPOqRqsx6jTjHYsrT60NrtAc/HyxYA+hjB1i9fOgXmDD1+/8/GpG0aAIx5FpzQOIVA3GLyy8fi+/45Lorlv9oChyo5ohaK+tKzWHnXZGzY2aJILugxsPhpTFefe1PYyqfXs7h+6kjs2BdvtXJZAOCr/65FNMpDaREFfSypyEdTvR29kFK6SIh6aDaJgaH479NFfa2OBa69qgRjirNx+FQPamN9xmBw1/Ri0fWV2HSgGXZ0K2RmY98elJWb/6hd6lPdsTLkeR7/a/PplHWt1wlt18HHbRsQG6LFFr2sWDgR7287KdVT0KSHVzYokH9/n2oJ9sbtp3DgYI9kG26YNhI7D7Yhu8OPxnaXor10xoKrok90yy3jcfRgDzLyvdLWOJdO6Ee1zKbniB0PL7oKXx/sSfrNeh2Dpx6/HigvwX9srcW3pxN9jr7Yd7v1LBbcNRlfxZ6n7jP6ZXZPhGWAOxZOxOSrC6Rl5P/v4+M4Wq9s/8nkU9se0W7v/vtJ8DzgMdjxywen44UPT8PJMEn7j9tmlmH+9DHgojxWrd2FcIBBhIui+2QPdCwjDUblZalu3/IcPMkQ2o+wzayz3o7sTCO+V1GCbw6dgoMX/Kaykhw0tPen9DfdeiEXmbp9OgC0sQAQ959YBlgwswwddh++lfmihfYw7D1eZGcaMbmsGF990wEuysOtt+PRuybjy72d2npjiG8L2vW3E4jySLAjniO9WLl4Mt78vCWWtweY+r1R2H2kAzZrBkYVZkvJSpORqh3qWAaVI3NwusWZ9BqWAcquKMS3p3o0h8Q5nUFUVo3Ax181Su2TgXBc5nVTSrFh8xFpG5jEGRfs/cmPAb6q0oZnfiiMSw7WCltMcjqDqG9zSm2XiW3JS7KYAx6DHfOuGY2tB7o17aroH7Mso7ktEABcsXHCrsNtUvJyAMjvCcHtCyHij2DJ7Cr8V20/9p3oQk6mERPKbNh7vFPwSXUMRhZmoSO2xc5hYHH9lFLUHhx4LCOXs7QyH728sC03GrMR994yHu8eTN2niiTzPXNG5WLezDJFf5hdnI3mLrfQ7vUsbpw2Ctv2NUv3yMcDlwKtrRa8fxb3sQNfcmHp6OhAUVERdDohUqzT6TBixAh0dKjNxOVDbaMDOw8Jw4RX/rIf2RYjDHpWGhzI9yVVldvw0mOzMS12Lq+cSCSKriRHWqXoV1Db1Id+bwgdvV78c28T9DoW8h1qDIDxo3MVR/y1dHmkvelMTM7ZU0qlE02gul8egY7ywn8sy+CK0fElrXyUx8gRwhYYcf+n2xfCpMoCxbNGpMiMzgPYc6xTKENDvAxzZEfHRnkeX3wr7Pvr7PWieqKwJ1E9sJMTjfLQsYzwPJ0wMBI/iYtCOlseECKgmWYDls2/AgyEGYBfr9+NwyoHSBywPrZkimZ9asFFgRONfQpnRNw3LuZj2Htc6BxZBpgeO11FDFYAQrmyDDTrSovapj5s2nEGNWfO/YxxUf5Ujo6cilJBH/Ye68L7/zyVMljBQKibsSMTT4LhOF5657Z9Ldh5sFWqPwZC4srkZ00Ig6+TGoMLjuMTHASGAQrzzLhtVjnuvWUCAGHmZ/3fj4KXfffV4wsVDrxOL0jQ2OFOO1gBCPUpntAgp63HqwiOlBQkJhYU950eqetNcNgYICEZYZQXZiUc7gCmjCvAbbPKMWtyKQCh7XFRXtI1HvEEjHJbVlVuwwPfvzKl/okD6l017Zg8tgALr6vAD+aOAwAcq7Pjl2/swslmJ0IyJ4vneZiMeqgPDQmHo/DHVr+ov5Hjosgw6MBx2sEKAHjtvQPYebBN87mN7S74AhFs+vyMFKxgZDYCENtr4nOFmZCB2+H40bkw6NmU+qkmqAoStHR5pG/XWLyAt/5xFMfrE11G0fbPuLIYDIBOrWAFhGPX1PlITrc4EZYlWRTLZeGsciy9aZyUr0lNhBPsvtjPiVfJ90i/9Nhs3L9wIl56bLZCr9WyHKnrhdZIRLSZ4tJpQHvui4tCEaxIhVY/qfW8r4904P1/njqrYAUglIM406qWOcpDEawAAF+Qk64Vy/CoRl0nyJqkTcgnGt7fdhIrF0/GgpllYBgoghVacos5DwAhf5C8jxa3a9a3uRLeG43ywkqa2LX5VjMYxPP4iPImM5uRSBR7jnZq/1H2DlG/7rutCkZRttiABrLvjnBR7K5JHmgU7Z7Ux6js37L541FVbpO2/CaDgbAcXOt7RLut9Y3J+g+WAUbYLLF2Klyz93gXDp7skeRcMLMct1xbJimXVvuuKrfh1hllKWUHhP7gm6OdcPvCaO/x4s3YUcPC34QEpgN1daFIFDpWmAR44u6puLIibmO4aDx5q1jG86rH4L4FVQq7Wtso+Lit3R7FKUacqh4ZCDmhRFL1j/Jr3L4QHvvBZGn32mcHWhEMceh1BlCQa1b4GtPGF6a0f1plmCpYIX73dVNKYUjSl9S39eO/tp+OBysYwff80YIq3DarHC89Nhu3zSpTbA9LFawAgJNNTinBq8moAwPgTKssWAFBbxbMKE/ad0UiUbT1JD/dTPQpkgUrxGe4fSHcf9tEaSUkAJxqdqKj1weXN4RfvrELB092AxBWjOTnxLfgcByvyAcUDkfx+SF5iFYbsT8bWZAJHsDRejuONzik7+e4KI4mmZSRo+WXyTnT6sSGzUewcvFk3DqjDAyEFdVyW+RXbTc7l0SflxJDvsLiu8iRul7JgIqNT0wap5WUp6rchh8tqMKxBrsUFRWNVoZBL9/yq4mQRCgWtW12KpyeaJTHgpljAAgDOzED+Mq7JgMA/uPTWkXkGhCMkphRu7zEih37m6V7WZbBLdeWoXJkTkIUl4/yGDdKyLAszpbdOqNM8bPY6IwGVvrdsvnjtSPCMTgusQwBIUmYeqDMcUInL5YZA+CK0bmob3eB4+Jla4glkhMTJQGQTvmQP1EePDhS1ytN28mDBJFIvFzERHblJVZF1ghR2QAAFIxJREFUfaYDAyE6LH/Oxu2npMSlgGDMWNnqGlG+edVjMK96DHbsF85z5zhZACSFfoiISfpEx0ys4/o2J7bta0EkEv8O8Vqtwb38fQzDSPeJZV5Vlo/GDveATo1Ox+DWWDkAyrphGGFVjejEcZywJ9kg06klc8cl1Sl5QO64qs1pyS3X27Yej6RbYlnxvJCASp2F/kSDQ7Ptyt8jtqnqiUU4cKIbXFScZS7FiUY7IrHtCuJspFjnBkPiN6plNRpYxd8MehZLbxqHP/29RtpCJMJx8cHj9IlF+ODzMwlbosT3ytuN3OF96bHZkq2IxBLBMTE9EZEPUiNcPFmcWpdEvb5uSmmCXZTXvVa9WjOVCfdYBhg3OhenW5zgY20XELYYDPRc0Un70YIqHKnrlZ6R7N3ydii3mWJdy22vmKxNbGOcWB4QdEo+EyfasuKCTHT2ehWbRnQ6pc6K5cwwwlGNUnJG0e7dNRlH6nqx70TiLKR4jdaJCvLBlLxcxOtmTipR1L9aJ+X9nLo/qCq3Kd63cfsp6Rr5OyaPLRASH6va7C0yW/FVTYciWD1+dC4aOlyK8hH7hbr2foV+qstBrCstu3o+EMsnNyv1nvdkdZ+svWu9R23btBD9lSKbRZoBF8tKXoZqWwPEbYDYRx+p601q/9T9r7xvFb9N1GuFjaztltqRlm3QqdqaXL/kstWc6cW7n5yQZGIZJuF5yeRNlWDxpmtGY9veZikRpryfFJ9z97wrUtptgz7ejzEMk+DL2HJMinIQ7zuu0d9EozwK88yYPLYA2/c3a7YpkXnVo6VrBvIB5M/Xaaw+0YKJ1a8YIBMnAcpLrAm+l5Z9eemx2Xgv5q+qfVx5PyzWo/itS2PlrbY5Bo32Iq+LI3W9CatAtHwNUUYt+5cO4greTLMhQbfKS6zCMelJVlokKyvRns6rHiOU2enEVdbyPk/wZ+J987EGe0J7FAMiABR6om6XN04biZNNDoQHKAMG8RUh8j5Q3l+89FPtOpf7C4Kt52M5xRJ9FvlWWfk75W1ZXu5q2yDep9YtsY1YTHp8sLMO0di2xVS+p4hoZ8WTgOS+nV7PwmoxKtq9PJnu5cyQByxKSkrQ1dUFjuOg0+nAcRy6u7tRUlIy1KJdMESnSm4g1Q6ZGnmHKm88gNI4yDtwsUNJNuiWO9Gi8VIHTe5bUIWjsiON1M6q3PCp71UHM+QOu/za8hJrwr3qAI54jejA97mDCZ2yugwfWzIFb/53TUJEXt1hqQcIyZwN+YBLXsbybPpyZ0LrW1PVp/z/6m/UepeWLmkZTPk9Yl2JJwEMpB9qZ1/rW8RvVOtlskCW/Jpk+qwlg+gcQeMEA62ylDshqfROLG91eaivSSW3fMAk1wGtwbu87sRrU5WPPAt5MvnFb9UKjCWTVcueyIOQWm1M0t2fDtwW1GjZCrWeyN+jLh+1bUv2jVrloS5XUcdYlsFjS6ZITrGWvqifqw7yyO2hXN7qiUXIyzYl6JS8HcrLQd22tNqY+lqxzaoDcWpd0CpnrW+Uv1+fQje16lhti9SDHnX9a71TPXDUek+qa5L1k/JrfhrrF0QHUh500Araq21lMrnVdlU9kFPbMrcvBK8/LDmz6r8ls4tsbADjC4QV8iar+2TtXcsZF3VL3RdptU95Xafbh6r76FT2T32v+lqtPimVjUxmV5PJJgZ4RBshDp61yjDd0x9E26luy1o2eCC7rfb/xHrQKgcgbk+1AiHn2u7kQVj1JMTKxZMVOqSu62R9tjqYpG6fWvZFa2JPHfzS8jlT+aCp6lodHE3l9yWzf8naWrp96n0LqnBc5vNpTQQlO7JXHSTWqhutOtHyXbRsTTrtUu57qW2NXO5kZapV5+rAZLKJAvnkaqpgnRq9jpHGRPLJs2S6NHNSSUq7lI6dVds9ebufVz0mpbyXC8PilJAVK1bg7rvvlpJu/u1vf8O777474H3iKSGDzTQ6HDifR6LJnwWkPhZvoOOikj1/oKPu0pXvfEYB03lusu89F5lS3Xu+v3Uw33g+v+18HcN1LkeVnqsMF+NM9HN97/mScaja2IWWfzC27VzlPhd7ku57LwSDaf8XWjeHqs0NlnSPAz5f70lly9J9Zzo6ezHtwPl411DYygsh21ByNmV4Po/ZTPf5Z+PLDPbvF+MbL/S7ztXuAoM/SvVs+oyh9OuT3ZNOYPJ82LJzHRMN9B3pyHap2CctznbsPiwCFnV1dVi1ahVcLhesVitefvllVFZWDnjfpRywIAiCIAiCIAiCIIjvAmc7dh/yLSEAMHbsWGzcuHGoxSAIgiAIgiAIgiAIYpgwLAIWZwvHCdmwOztTZ4MmCIIgCIIgCIIgCGJoEMfs4hg+XS7pgEVPj3B6xX333TfEkhAEQRAEQRAEQRAEkYqenh6UlQ18XLLIsMhhcbYEAgEcPXoUhYWF0Ol0Qy0OQRAEQRAEQRAEQRAqOI5DT08PJk2aBJMp9ZHdci7pgAVBEARBEARBEARBEJcn7FALQBAEQRAEQRAEQRAEoYYCFgRBEARBEARBEARBDDsoYEEQBEEQBEEQBEEQxLCDAhYEQRAEQRAEQRAEQQw7KGBBEARBEARBEARBEMSwgwIWBEEQBEEQBEEQBEEMOyhgQRAEQRAEQRAEQRDEsIMCFheAhoYGLF++HAsWLMDy5cvR2NiYcA3HcXjhhRdw880345ZbbsHGjRsvvqAEcRako9+7du3CkiVLMGnSJLz88ssXX0iCOEvS0e+1a9fi9ttvx6JFi7BkyRJ8+eWXF19QgjgL0tHvTZs2YdGiRVi8eDEWLVqEd9555+ILShBnQTr6LVJfX4+pU6eSj0JcEqSj22vWrMGsWbOwePFiLF68GC+88MLFF/QCwfA8zw+1EJcbDzzwAJYuXYrFixdj8+bN2LRpU0KH/8EHH+DDDz/Ehg0b4HQ6cdddd+G9997DqFGjhkhqgkiPdPS7qakJPp8PW7duRSgUwrPPPjtE0hLE4EhHv7/88ktUV1fDbDajtrYW999/P3bt2gWTyTREUhNEeqSj3x6PB5mZmWAYBh6PB4sWLcK6detQVVU1RFITRHqko9+AMGn44x//GCNGjMCIESPIRyGGPeno9po1a+Dz+S5LfaYVFucZu92O48eP44477gAA3HHHHTh+/DgcDofiui1btmDZsmVgWRY2mw0333wztm7dOhQiE0TapKvfZWVlmDhxIvR6/VCISRBnRbr6fcMNN8BsNgMAJkyYAJ7n4XQ6L7q8BDEY0tXvrKwsMAwDAAgEAgiHw9LPBDFcSVe/AWD9+vWYO3cuysvLL7KUBDF4BqPblysUsDjPdHR0oKioCDqdDgCg0+kwYsQIdHR0JFxXWloq/VxSUoLOzs6LKitBDJZ09ZsgLkXORr8/+OADjBkzBsXFxRdLTII4Kwaj39u3b8ftt9+Om266CY888ggmTJhwscUliEGRrn7X1tZi165d+PGPfzwEUhLE4BmM7f7444+xaNEiPPzwwzh06NDFFvWCQQELgiAIgjgL9u7di9dffx2vvfbaUItCEOeV+fPn4+OPP8ann36KzZs3o76+fqhFIohzJhwO49e//jVeeOEFafBHEJcL9957L7Zv344PP/wQ//Iv/4LHH38cfX19Qy3WeYHWa59nSkpK0NXVBY7joNPpwHEcuru7UVJSknBde3s7pkyZAiBxxQVBDEfS1W+CuBQZjH4fOnQIzzzzDN544w1UVlYOgbQEMTjOxn6XlpZi8uTJ+Pzzz0nPiWFNOvrd09OD5uZmPProowAAl8sFnufh8Xjw29/+dqhEJ4iUpGu7CwsLpX/Pnj0bJSUlOH36NK699tqLLfJ5h1ZYnGfy8/MxceJEfPTRRwCAjz76CBMnToTNZlNcd9ttt2Hjxo2IRqNwOBzYtm0bFixYMBQiE0TapKvfBHEpkq5+19TU4Be/+AX+/d//HVddddVQiEoQgyZd/a6rq5P+7XA4sGfPHowfP/6iykoQgyUd/S4tLcWePXuwY8cO7NixAw8++CDuueceClYQw5p0bXdXV5f07xMnTqCtrQ0VFRUXVdYLBZ0ScgGoq6vDqlWr4HK5YLVa8fLLL6OyshIrV67Ek08+icmTJ4PjOKxevRpfffUVAGDlypVYvnz5EEtOEAOTjn7v378f//qv/wqPxwOe55GdnY2XXnoJN9xww1CLTxApSUe/ly5dira2NhQVFUn3vfLKK7TPnxj2pKPfv/vd7/DVV19Br9eD53ksW7YMK1asGGrRCWJA0tFvOZfzqQrE5UU6uv3ss8/i2LFjYFkWBoMBTz75JObMmTPUop8XKGBBEARBEARBEARBEMSwg7aEEARBEARBEARBEAQx7KCABUEQBEEQBEEQBEEQww4KWBAEQRAEQRAEQRAEMeyggAVBEARBEARBEARBEMMOClgQBEEQBEEQBEEQBDHsoIAFQRAEQRCXDX/961/x0ksvDbUYEj//+c+xc+fOoRaDIAiCIC5JKGBBEARBEJcw8+bNw5QpUzBt2jTpv66urqEWa0gIhUJYt24dHnnkEQBAa2srJkyYgEgkorhu1apV+MMf/nDW71mxYgU2btyY1rUrV67E66+/ftbvIgiCIIjvMvqhFoAgCIIgiHPjzTffxHXXXZfymkgkAr3+8u72t2/fjsrKShQVFV2Q5/M8D57nB3XPlClT4PF4cOTIEUyePPmCyEUQBEEQlyu0woIgCIIgLkPE1QUbN27E3Llz8eCDDwIAvv32W9x7772orq7GnXfeiT179kj3tLS04P7778e0adPw0EMPYfXq1Xj66acBAHv27MGNN96oeMe8efOwe/duAEA0GsX69etx8803Y8aMGXjqqafgdDoVsvz973/H3LlzMWPGDKxbt056DsdxePPNN3HzzTdj2rRpWLJkCTo6OgAAdXV1eOihh3DttddiwYIF2LJlS9Jv/uKLLzB9+vRBl1WqMlmxYgX+8Ic/4N5778XUqVPxzDPPYP/+/Vi9ejWmTZuG1atXY8OGDYoVLldddRVWrVolPePaa6+lbSEEQRAEcRZc3lMtBEEQBPEdZ9++fdiyZQtYlkVXVxd+8pOf4JVXXsENN9yAr7/+Gk8++SQ++eQT2Gw2PP3007j66qvx5z//GYcPH8ajjz6K+fPnp/Wed999F9u2bcNf/vIX2Gw2vPjii1i9ejV+//vfS9ccOHAAW7duRWNjI+6++27ceuutGDt2LN5++218/PHHWL9+PSoqKnDy5EmYTCb4fD48/PDDePLJJ7FhwwacOnUKDz30EMaPH49x48YlyHDq1KmEoMpADFQmALB582Zs2LABFRUV4Hke3d3duPPOO7Fs2TLpOStXrgQAdHR0YNmyZVi4cKH0t7Fjx+LAgQODkosgCIIgCFphQRAEQRCXPE888QSqq6tRXV2Nxx9/XPG3n//857BYLDCZTNi8eTNuvPFGzJkzByzLYvbs2Zg0aRJ27tyJ9vZ2HDlyBE899RSMRiOmT5+OefPmpS3DX//6V/ziF79AcXExjEYjfvazn+HTTz9V5I/42c9+BpPJhKqqKlRVVaG2thYAsHHjRjz11FOorKwEwzCoqqpCXl4ePv/8c4wcORJLly6FXq/HlVdeiQULFmDr1q2aMrjdbmRmZib8fubMmVL5VFdX46OPPpL+lqpMRH7wgx/giiuugF6vh8FgSFoGgUAATzzxBB544AHMmTNH+n1mZiZcLlfaZUkQBEEQhACtsCAIgiCIS5y1a9cmzWFRXFws/bu9vR1bt27FZ599Jv0uEolgxowZ6O7uhtVqhcVikf5WWloqbc0YiPb2djzxxBNg2fhcCMuysNvt0s8FBQXSv81mM3w+HwCgs7MTY8aMSXhmW1sbampqUF1dLf2O4zjceeedmjJYrVZ4vd6E33/zzTeK/B3y7RqpykSkpKRE+6NVPPfcc6ioqMCjjz6q+L3X64XVak3rGQRBEARBxKGABUEQBEFcxjAMI/27pKQEixcvxosvvphwXVtbG1wuF3w+nxS0aG9vl+43m80IBALS9RzHweFwSD8XFxfjd7/7Ha655pqEZ7e2tqaUsbi4GM3NzRg/frzi9yUlJZg+fTrefvvtNL4UmDBhAhobG9O6Vv6OZGUiIi/DZKxfvx4NDQ147733Ev5WV1eHqqqqQclFEARBEARtCSEIgiCI7wx33nknPvvsM3z55ZfgOA7BYBB79uxBZ2cnRo4ciUmTJmHNmjUIhULYv3+/YtVBRUUFgsEgPv/8c4TDYaxbtw6hUEj6+w9/+EP88Y9/RFtbGwDA4XBg27Ztacm1bNkyvP7662hsbATP86itrUVfXx/mzp2LxsZGfPDBBwiHwwiHw6ipqUFdXZ3mc+bMmYN9+/adtzJJRkFBAVpaWqSfd+7ciXfeeQdr166FyWRKuH7fvn2Dzq1BEARBEAQFLAiCIAjiO0NJSQneeOMN/OlPf8KsWbMwZ84cvPXWW4hGowCA1157DYcPH8aMGTOwdu1a3HXXXdK92dnZeP755/GrX/0KN954I8xms2K7yQMPPIB58+bh4YcfxrRp03DPPfegpqYmLbkeeughLFy4EA8//DC+973v4bnnnkMwGERWVhbeeustbNmyBTfccAOuv/56vPrqq4pAiZybbroJ9fX16OrqOm9losUDDzyATz/9FNOnT8eLL76ITz75BH19ffj+978vnRTym9/8BgBQU1MDi8WCKVOmpC0TQRAEQRACDD/YA8UJgiAIgvhOsGbNGjQ1NeHVV18dalHS5v3338eZM2fw3HPPDbUoAISkp3fffbciCSdBEARBEOlBOSwIgiAIgrhsWL58+VCLoGDNmjVDLQJBEARBXLLQlhCCIAiCIAiCIAiCIIYdtCWEIAiCIAiCIAiCIIhhB62wIAiCIAiCIAiCIAhi2EEBC4IgCIIgCIIgCIIghh0UsCAIgiAIgiAIgiAIYthBAQuCIAiCIAiCIAiCIIYdFLAgCIIgCIIgCIIgCGLYQQELgiAIgiAIgiAIgiCGHf8fBkQL3kmAX9AAAAAASUVORK5CYII=\n",
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# The circuit we are interested\n",
"opstr = pygsti.objects.Circuit(None, stringrep='Gx(Gi)^128Gy')\n",
"# We hand the circuit to the plotting function\n",
"results_gst.plot_spectrum(sequence=opstr,loc='upper right')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Box-plots for GST data**\n",
"\n",
"If the data is from GST experiments, or anything with a GST-like structure of germs and fudicials, we can create a box-plot which shows the maximum power in the spectrum for each sequence. This maximum power is a reasonable proxy for comparing how \"drifty\" the data from the different sequences appears to be. But note that the maximum power should *not* be used to directly compare the level of drift in two different datasets with different parameters, particularly if the number of timestamps is different - because this maximum power will increase with more data, for a fixed level of drift. More on this at a later date.\n",
"\n",
"In the plot below we see that the amount of drift appears to be increasing with sequence length, as would be expected with gate drift. Without performing a detailed analysis, by eye it is clear that the $G_i$ gate is the most drifty, that the $G_x$ gate has some drift, and that the data looks consistent with a drift-free $G_y$ gate."
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# This box constructs some GST objects, needed to create any sort of boxplot with GST data\n",
"from pygsti.construction import std1Q_XYI # The model used with the GST data we imported\n",
"\n",
"# This manually specifies the germ and fiducial structure for the imported data.\n",
"fiducial_strs = ['{}','Gx','Gy','GxGx','GxGxGx','GyGyGy']\n",
"germ_strs = ['Gi','Gx','Gy','GxGy','GxGyGi','GxGiGy','GxGiGi','GyGiGi','GxGxGiGy','GxGyGyGi','GxGxGyGxGyGy']\n",
"log2maxL = 9 # log2 of the maximum germ power\n",
"\n",
"# Below we use the maxlength, germ and fuducial lists to create the GST structures needed for box plots.\n",
"fiducials = [pygsti.objects.Circuit(None,stringrep=fs) for fs in fiducial_strs]\n",
"germs = [pygsti.objects.Circuit(None,stringrep=mdl) for mdl in germ_strs]\n",
"max_lengths = [2**i for i in range(0,log2maxL)]\n",
"gssList = pygsti.construction.make_lsgst_structs(std1Q_XYI.gates, fiducials, fiducials, germs, max_lengths) "
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"collapsed": false
},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"\n",
"
\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
},
{
"data": {
"text/plain": [
""
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Create a boxplot of the maximum power in the power spectra for each sequence.\n",
"w.ColorBoxPlot('driftpwr', gssList[-1], None, None, driftresults = (results_gst,None))"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# NOTHING BELOW HERE WORKS"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**Constructing estimates of the drifting probabilities**\n",
"\n",
"The analysis creates estimates of the time-dependent probability, $p(t)$, to obtain a given outcome, for each sequence. Below, we will explain how to access and plot all of these estimates. Here, we demonstrate how to plot the estimated $p(t)$ that the analysis concludes is \"the most drifty\" (the $p(t)$ where the estimated drift has the highest power). Obviously, the more this function oscillates, the more drifty the data for this sequence appears to be."
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"#results_gst.plot_most_drifty_probability(plot_data=True)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"It can also be informative to display multiple reconstructions on a single plot. Below we plot the reconstructions for the $G_i$ and $G_x$ germs, both with varying germ power and fixed fudicials."
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Pick a list of sequence labels (here circuits) or indices to plot the estimated p(t) for\n",
"gstrs = [pygsti.objects.Circuit(None,stringrep='Gx(Gi)^'+str(2**l)+'Gy') for l in range(4,10)]\n",
"# Hand this list to the plotting function\n",
"#results_gst.plot_probability_trajectory_estimates(gstrs)"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {
"collapsed": false
},
"outputs": [],
"source": [
"# Pick a list of sequence labels (here circuits) or indices to plot the estimated p(t) for\n",
"gstrs = [pygsti.objects.Circuit(None,stringrep='(Gx)^'+str(2**l)) for l in range(4,10)]\n",
"# Hand this list to the plotting function\n",
"#results_gst.plot_multi_estimated_probabilities(gstrs,loc='upper right')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"**More details please?**\n",
"\n",
"Having provided an overview of the tools in the `drift` module, we now give a more detailed introduction to:\n",
"1. The types of data that can be analyzed, and the format they must be imported in to use the `drift` module.\n",
"2. The analysis methods that are used inside the `drift.do_basic_drift_characterization()`.\n",
"3. How to access the various outputs returned by the analysis, and some details on how to interpret them. "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Input data types and formats\n",
"#### What types of data can be analyzed?\n",
"\n",
" The methods in this notebook can be used on data from any experiment that satisfies the following criteria: \n",
"\n",
"- The experiments consist of $S \\geq 1$ different circuits, or sequences, each with $M \\geq 2$ possible outcomes with $M$ the same for all $S$ circuits.\n",
"\n",
"\n",
"- Sequence $s$ is repeated, and one of the $M$ outcomes is recorded, $N$ times during each of $T$ different time-intervals, with $N$ and $T$ the same for all sequences $s$. Let $\\tau_{s,1},\\tau_{s,2},\\dots,\\tau_{s,T}$ denote times associated with these time-intervals (e.g., the start time or mid-point of each interval).\n",
"\n",
"\n",
"- The time-gap between consecutive sets of $N$ repeats of a sequence is approximately constant, and approximately independent of sequence. That is $t_{s,gap} \\approx \\tau_{s,t+1} - \\tau_{s,t}$ for all $s,t$ and some constant $t_{s,gap}$ (the time gap is constant for each sequence), and $t_{gap} \\equiv t_{s,gap} \\approx t_{s',gap}$ for all $s,s'$ (the time gap is the same for every sequence).\n",
"\n",
"#### What type of circuits can the data come from?\n",
"\n",
"Many characterization routines, including GST, RPE, restricted sets of GST sequences (e.g., Ramsey sequences) or RB, satisfy the above criteria *if* a suitable time-ordering for the repeats of all of the sequences is chosen.\n",
"\n",
"An example is a full set of GST sequences, with sequence 1 repeated $N=5$ times, then sequence 2 repeated $N=5$ times, etc, with this entire procedure looped through $T = 500$ times. As long as each circuit follows on from the following circuit with a roughly constant time gap, then -- even if the different circuits take varying times to implement -- the time-gaps between consecutive sets of $N$ repeats of each circuit are all the same.\n",
"\n",
"To obtain best performance, it is preferable to minimize $N$ and maximize $T$ for the same fixed $N \\times T$ (this will provide better detection for high frequency drift).\n",
"\n",
"\n",
"#### What input formats can be used?\n",
"\n",
"Data can be provided to the analysis function in a `pyGSTi` dataset, which can be imported using the `pygsti.io.load_tddataset()` function, demostrated above. The automated input for this currently has limited functionality, and the analysis is much slower when using this form of input. However, this method has the advantage of simplicity: minimal user specification is required when the analysis uses such a dataset.\n",
"\n",
"2. Data can be provided to the analysis function as an ordinary numpy array. The analysis is faster when using this format. The array can have a range of different dimensions, but the most flexible format is when the array is 4 dimensional. \n",
"\n",
"Here we will explain the format when the array is 4 dimensional, below we will demonstrate the method with lower dimensional arrays. When the input array is 4 dimensional, the input data format is a numpy array $A$ of dimension $(S \\times E \\times M_{\\rm marg} \\times T)$, where\n",
"1. $S$ is the number of sequences,\n",
"2. $E$ is the number of \"entities\" (often 1; and more on this below),\n",
"3. $M$ is the number of possible measurement outcomes for each circuit,\n",
"4. $T$ is the number of timesteps.\n",
"\n",
"The matrix element $A[s,e,m,t]$ should record the number of counts (a value in $0,1,\\dots,N-1$) for measurement outcome $m$, on entity $e$, for sequence $s$, with the $t^{\\rm th}$ set of repeats for that sequence.\n",
"\n",
"The aspect of this array that is perhaps least inuitive is the \"entities\" concept. This allows the storing of independent data for the same circuits, timestamps, and the same possible measurement outcomes for more than one \"entity\". For example, if multiple single-qubit experiments are implemented in parallel, and any correlations between measurement outcomes are discarded, the data from each qubit can be stored with different entity indexes. This then allows for two things:\n",
"\n",
"1. The analysis function independently looks for drift in each entity.\n",
"2. The analysis function combines the data from all entities, to see if that has evidence for drift which was not evident in the analysis of each entity independently. For example, weak evidence for drift in each entity can combine to convincing evidence for drift in the device as a whole.\n",
"\n",
"Later we will discuss further how to use this \"entities\" concept to optimize data analysis, via examples.\n",
"\n",
"As shown below, the input data is stored in the results object as an $(S \\times E \\times M \\times T)$ array, \n",
"which is the same format as can be used for the input. In our GST example there are 3121 sequences ($S = 3121$), a single entity ($E = 1$), two-outcome measurements ($M=2)$, and 500 timesteps ($T=500$)."
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# import numpy as np\n",
"# np.shape(results_gst.data)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## The analysis methods\n",
"\n",
"#### What is the aim of the methods in this tutorial?\n",
"\n",
"The analysis presented in this tutorial has two purposes: to detect and characterize drift. To explain this further, we need to define what we mean by drift.\n",
"\n",
"There is an underlying probability for circuit $s$ to output the measurement outcome $m$ at time $t$. Denote this by $p_{s,m,t}$, and assume these probabilities satisfy $\\sum_m p_{s,m,t} = 1$ for all $s$ and $t$ (i.e., one of the $M$ measurement outcomes is always observed). By definition, there is no drift if $p_{s,m,t} = p_{s,m,t'}$ for all $s,m,t,t'$. Otherwise there is drift. \n",
"\n",
"In the experiments considered herein, we are essentially taking $N$ samples from a multinomial probability distribution with probability vector $(p_{s,0,t},p_{s,1,t},\\dots)$ for each $s$ and $m$, at a set of equally-spaced times. As such, the data allows us to estimate properties of the vectors:\n",
"\n",
"$$ \\boldsymbol{p}_{s,m} = (p_{s,m,t_1},p_{s,m,t_2},\\dots,p_{s,m,t_3}).$$\n",
"\n",
"The aim of our analysis is:\n",
"\n",
"1. **Drift Detection:** Starting from the assumption that there is no drift, decide whether there is statistically significant evidence in the data that $p_{s,m,t}$ is *not* independent of $t$.\n",
"2. **Drift characterization:** Estimate $p_{s,m,t}$ for all sample times $t$.\n",
"\n",
"The drift detection and characterization works via spectral analysis, by implement a type of discrete Fourier transform on the data and looking for peaks in the power spectrum that are larger than would be expected if there was no drift.\n",
"\n",
"Specifically, consider a data vector $\\boldsymbol{x}=(x_0,x_1,\\dots,x_{T-1})$, of the number of counts observed for measurement outcome $m$ with sequence $s$ as a function of time $t \\in [0,1,\\dots,T]$. We remove the data mean from this vector, and renormalize it in a particular way, to obtain the $\\boldsymbol{z}$ vector, defined by:\n",
"\n",
"$$\n",
"z_{s,m,t} = \\frac{x_{s,m,t} - \\bar{x}_{s,m}}{\\sqrt{ N\\bar{x}_{s,m}(1- \\bar{x}_{s,m})}},\n",
"$$\n",
"\n",
"We won't fully explain the reasons for this mapping on $\\boldsymbol{x}$ here, except to note that (1) it is not essential, and (2) the motivation is that -- after this mapping -- the statisical behaviour of the random variable from which $\\boldsymbol{z}$ is drawn, when there is no drift, is largely independent of $N$ and the value of the time-independent probabilitiy $p_{s,m}$.\n",
"\n",
"We then perform a Fourier transform on the data -- specifically, we implement a Type-II orthogonal discrete cosine transform (DCT) -- to obtain Fourier modes $ \\tilde{\\boldsymbol{z}}_{s,m} = F \\boldsymbol{z}_{s,m}$. The precise definition of this transform is not important for obtaining an intuitive understanding of this method, but for completeness it is defined by:\n",
"\n",
"\\begin{align}\n",
"F_{\\omega,t} =\\sqrt{\\frac{2^{1-\\delta_{\\omega,0}}}{N}} \\cos\\left(\\frac{\\omega \\pi}{N}\\left( t + \\frac{1}{2}\\right)\\right), \n",
"\\end{align}\n",
"with $\\omega,t=0,\\dots,N-1$.\n",
"\n",
"We then square these Fourier modes, to obtain a power spectrum $\\tilde{\\boldsymbol{z}}_{s,m}^2$ for each measurement outcome and each sequence. It is these power spectra that we analyze to detect and characterize drift. We can look at the individual power spectra for a particular sequence and measurement outcome. Using the central limit theorem and some statistics (which we will not cover here), it can be show that - when there is no drift - each Fourier mode is approximately a zero-mean unit-variance $\\chi^2_1$ random variable. As such, we can calculate the probability that a spike in a power spectrum is due to random statistical flucations, and if this probability is low enough then we can be confident the peak is due to drift. I.e., the underlying $p(t)$ is not constant.\n",
"\n",
"We can also average these power spectra over measurement outcomes, and/or over sequences (and over entities, which haven't been discused above, and on which this entire analysis is independently performed). When averaging over all of these quantities, we obtain the global power spectrum presented above. The advantage of this averaging is that it supresses noise from random flucations, but does not supress any \"signal\" due to nonconstant $p(t)$ - as long as the same frequencies are present in many of the power spectra that are being averaged.\n",
"\n",
"Before explaining any more of the analysis methods, we turn to a simple example."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Example 1 : Single sequence data\n",
"We demonstrate this on single-sequence data, for a single entity, with two-outcome measurements. This is the circumstance under which the analysis function can be given data in a 1D array, and so we will use this method below. We now demonstrate the method with simulated data for a drifting and drift free probability.\n",
"\n",
"Below we use a drifting probability of\n",
"$$ p(t) = 0.5 + 0.2 \\cos(0.1 t).,$$\n",
"with integer $t$ in the range $0,1,2,\\dots,T$.\n",
"This is similar in form to the drifting probability that would be obtained with certain sorts of Ramsey experiment with drifting $\\sigma_z$ over-rotation angle."
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# Imports for generating fake data to demonstrate the methods.\n",
"from numpy.random import binomial\n",
"from numpy.random import multinomial"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# N = 5 # Counts per timestep\n",
"# T = 100 # Number of timesteps\n",
"\n",
"# # The confidence of the statistical tests. Here we set it to 0.999, which means that\n",
"# # if we detect drift we are 0.999 confident that we haven't incorrectly rejected the\n",
"# # initial hypothesis of no drift.\n",
"# confidence = 0.999\n",
"\n",
"# # A drifting probability to obtain the measurement outcome with index 1 (out of [0,1])\n",
"# def pt_drift(t): return 0.5+0.2*np.cos(0.1*t)\n",
"\n",
"# # A drift-free probability to obtain the measurement outcome with index 1 (out of [0,1])\n",
"# def pt_nodrift(t): return 0.5\n",
"\n",
"# # If we want the sequence to have a label, we define a list for this (here, a list of length 1).\n",
"# # The labels can, but need not be, pyGSTi OpString objects.\n",
"# sequences = [pygsti.objects.OpString(None,'Gx(Gi)^64Gx'),]\n",
"\n",
"# # If we want the outcomes to have labels, we define a list for this.\n",
"# outcomes = ['0','1']\n",
"\n",
"# # Let's create some fake data by sampling from these p(t) at integer times. Here we have\n",
"# # created a 1D array, but we could have instead created a 1 x 1 x 1 x T array.\n",
"# data_1seq_drift = np.array([binomial(N,pt_drift(t)) for t in range(0,T)])\n",
"# data_1seq_nodrift = np.array([binomial(N,pt_nodrift(t)) for t in range(0,T)])\n",
"\n",
"# # If we want frequencies in Hertz, we need to specify the timestep in seconds. If this isn't\n",
"# # specified, the frequencies are given in 1/timestep with timestep defaulting to 1.\n",
"# timestep = 1e-5\n",
"\n",
"# # We hand these 1D arrays to the analysis function, along with the number of counts, and other\n",
"# # optional information\n",
"# results_1seq_drift = drift.do_basic_drift_characterization(data_1seq_drift, counts=N, outcomes=outcomes,\n",
"# confidence=confidence, timestep=timestep, \n",
"# indices_to_sequences=sequences)\n",
"# results_1seq_nodrift = drift.do_basic_drift_characterization(data_1seq_nodrift, counts=N, outcomes=outcomes, \n",
"# confidence=confidence, timestep=timestep, \n",
"# indices_to_sequences=sequences) "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can now compare the global power spectrum obtained when there is drift, and when there is not drift. Except in highly unusual cases (1 in 1000) there will be no drift detected when there is indeed no drift. "
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# results_1seq_drift.plot_power_spectrum()"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# results_1seq_nodrift.plot_power_spectrum()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"We can look to see what the p-value associated with the largest peak in the power spectrum is. "
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# print(results_1seq_drift.global_pvalue)\n",
"# print(results_1seq_nodrift.global_pvalue)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The analysis function performs:\n",
"1. an analysis on a per-sequence, per-entity, per-outcome basis (properties starting with `pspepo`), an analysis on a per-sequence, per-entity, outcome-averaged basis (properties starting with `pspe`), \n",
"2. an analysis on a per-sequence, entity-averaged, outcome-averaged basis (properties starting with `ps`), \n",
"3. an analysis on a sequence-averaged, per-entity, outcome-averaged basis (properties starting with `pe`),\n",
"4. a global analysis on a sequence-averaged, entity-averaged, outcome-averaged basis (properties starting with `global`)\n",
"\n",
"That is, it inspects all the power spectra, after the relevant averaging has been performed, for evidence of drift (e.g., for case 1 there are $S \\times E \\times M \\times T$ spectra).\n",
"\n",
"In the example here there is a single sequence, a single entity, and two-outcome measurements. Hence, in this case all of the analyzes are exactly equivalent, as is clear by noting that all of the power spectra are the same:"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# # The power spectrum obtained after averaging over everthing\n",
"# print(results_1seq_drift.global_power_spectrum[:4])\n",
"# # The power spectrum obtained after averaging over everthing except sequence label\n",
"# print(results_1seq_drift.ps_power_spectrum[0,:4])\n",
"# # The power spectrum obtained after averaging over everthing except entity label\n",
"# print(results_1seq_drift.pe_power_spectrum[0,:4])\n",
"# # The power spectrum obtained after averaging over everthing except sequene and entity label\n",
"# print(results_1seq_drift.pspe_power_spectrum[0,0,:4])\n",
"# # The two power spectra obtained after averaging over nothing\n",
"# print(results_1seq_drift.pspepo_power_spectrum[0,0,0,:4])\n",
"# print(results_1seq_drift.pspepo_power_spectrum[0,0,1,:4])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"As we can see above, the power spectra for the two different measurement outcomes are identical. This is because the Fourier modes for all of the measurement outcomes must sum to zero. Hence, there are only $M-1$ independent power spectra, which here is 1.\n",
"\n",
"As already demonstrated above, the analysis function also creates an estimate for the drift probability for each measurement outcome (and for each sequence, and each entity, but here there is only one of each of these). This estimate is created in the following way:\n",
"\n",
"1. Inspect the relevant power spectrum (the power spectrum for the sequence, entity, and measurement outcome under consideration).\n",
"2. Take the modes obtained from the data, and set all modes with power below the \"significance threshold\" to zero. This significance threshold is generally adjusted for the number of sequences and entities considered, which is not relevant in this case but is discussed below.\n",
"3. Keep all modes with power above the significance threshold and invert the Fourier transform.\n",
"4. This is the estimate of $p(t)$, up to some final adjustments to guarantee that $p(t)$ is within $[0,1]$.\n",
"\n",
"Because we have created our data from a known underlying $p(t)$ we can compare our estimate of $p(t)$ with the true $p(t)$. We do this below."
]
},
{
"cell_type": "code",
"execution_count": 21,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# # Lets create an array of the true probability. This needs to be\n",
"# # of dimension S x E x M x T\n",
"# parray_1seq = np.zeros((1,1,2,T),float)\n",
"# parray_1seq[0,0,0,:] = np.array([pt_drift(t) for t in range(0,T)])\n",
"# parray_1seq[0,0,1,:] = 1 - parray_1seq[0,0,0,:]\n",
"\n",
"# # The measurement outcome index we want to look at (here the esimated p(t) \n",
"# # for one index is just 1 - the p(t) for the other index, because we are\n",
"# # looking at a two-outcome measurement).\n",
"# outcome = 1\n",
"\n",
"# # If we hand the parray to the plotting function, it will also plot\n",
"# # the true probability alongside our estimate from the data\n",
"# results_1seq_drift.plot_estimated_probability(sequence=0,outcome=outcome,parray=parray_1seq,plot_data=True)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Example 2 : Single-sequence multi-qubit data\n",
"Single sequence data with two measurement outcomes is (a) often not will be obtained in experiments, and (b) is not sufficient to demonstrate all of the methods contained within the analysis function, or to understand all of the ouput. Hence, we now consider an example with a single-sequence and 4 measurement outcomes, which could represent a two-qubit single-sequence experiment.\n",
"\n",
"Let us consider the case of 4 possible measurement outcomes:"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# We only explicitly name them for labelling purposes.\n",
"# outcomes = ['00','01','10','11']"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"This could correspond to computational basis measurements of two qubits. Now let us assume that the drifting $p(t)$ for obtaining each measurement outcome factorizes into probabilities for outcomes 0 and 1 for each qubit (e.g., because it is parallel single-qubit experiments and there is no crosstalk). Specifically, take the first qubit to have drifting probability for outcome 0 of\n",
"$$ p_{0}(t) = 0.5+0.05 \\cos(0.08t),$$\n",
"and the second qubit to have drifting probability for outcome 0 of\n",
"$$ p_{0}(t) = 0.5+0.05 \\cos(0.2t),$$\n",
"with $t \\in [0,1,\\dots, T]$.\n",
"\n",
"Let's create some data from these drifting probabilities"
]
},
{
"cell_type": "code",
"execution_count": 23,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# N = 10 # Counts per timestep\n",
"# T = 1000 # Number of timesteps\n",
"\n",
"# # The drifting probabilities for the 4 outcomes\n",
"# def pt00(t): return (0.5+0.07*np.cos(0.08*t))*(0.5+0.08*np.cos(0.2*t))\n",
"# def pt01(t): return (0.5+0.07*np.cos(0.08*t))*(0.5-0.08*np.cos(0.2*t))\n",
"# def pt10(t): return (0.5-0.07*np.cos(0.08*t))*(0.5+0.08*np.cos(0.2*t))\n",
"# def pt11(t): return (0.5-0.07*np.cos(0.08*t))*(0.5-0.08*np.cos(0.2*t))\n",
"\n",
"# # Because of the type of input (>2 measurement outcomes), we must record the\n",
"# # data in a 4D array (even though some of the dimensions are trivial)\n",
"# data_multiqubit = np.zeros((1,1,4,T),float)\n",
"\n",
"# # Generate data from these p(t)\n",
"# for t in range(0,T):\n",
"# data_multiqubit[0,0,:,t] = multinomial(N,[pt00(t),pt01(t),pt10(t),pt11(t)])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Example 2.1 : A simple analysis of the full dataset\n",
"The simplest way to analyze this data is just to hand it to the drift characterization method, as before (note that now we don't need to specify the counts number, as the input is a 4D array). We do this below."
]
},
{
"cell_type": "code",
"execution_count": 24,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# results_multiqubit_full = drift.do_basic_drift_characterization(data_multiqubit,outcomes=outcomes,confidence=0.99)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"By looking at the global power spectrum, we can clearly see that there is drift."
]
},
{
"cell_type": "code",
"execution_count": 25,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# results_multiqubit_full.plot_power_spectrum()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"As always, we can also return the frequencies, found from inspecting the global power spectrum. As we have not specified a timestep, these are integers between $1$ and $T-1$."
]
},
{
"cell_type": "code",
"execution_count": 26,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# print(results_multiqubit_full.global_drift_frequencies)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Because there are multiple measurement outcomes, we can look at the power spectra for each of these, and also we can look at the drift frequencies found from just analyzing this spectrum. One example of such a power spectrum is given below"
]
},
{
"cell_type": "code",
"execution_count": 27,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# outcome = '00'\n",
"# results_multiqubit_full.plot_power_spectrum(sequence=0,entity=0,outcome=outcome)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The drift frequencies found via an independent analysis of each of these spectra. Generally, these need not contain the same frequencies as the global analysis. In this example, because all of the frequencies appear in all of the spectra, the global analysis is more sensitive (reduced noise without reducing signal strength). "
]
},
{
"cell_type": "code",
"execution_count": 28,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# print(results_multiqubit_full.pspepo_drift_frequencies[0,0,0])\n",
"# print(results_multiqubit_full.pspepo_drift_frequencies[0,0,1])\n",
"# print(results_multiqubit_full.pspepo_drift_frequencies[0,0,2])\n",
"# print(results_multiqubit_full.pspepo_drift_frequencies[0,0,3])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"As always, we can also look at the reconstruction for a given drifting probability. An example is shown below, where again we construct the true underlying probability with which to compare it."
]
},
{
"cell_type": "code",
"execution_count": 29,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# # Creates an array of the true probability.\n",
"# parray_multiqubit_full = np.zeros((1,1,4,T),float)\n",
"# parray_multiqubit_full[0,0,0,:] = np.array([pt00(t) for t in range(0,T)])\n",
"# parray_multiqubit_full[0,0,1,:] = np.array([pt01(t) for t in range(0,T)])\n",
"# parray_multiqubit_full[0,0,2,:] = np.array([pt10(t) for t in range(0,T)])\n",
"# parray_multiqubit_full[0,0,3,:] = np.array([pt11(t) for t in range(0,T)])\n",
"\n",
"# results_multiqubit_full.plot_estimated_probability(sequence=0,outcome=1, plot_data=True,\n",
"# parray=parray_multiqubit_full)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Example 2.2 : Marginalizing over qubits\n",
"\n",
"The above analysis looks for drift in the $2^2$ measurements outcomes of a two-qubit circuit. This would be problematic for many qubits, as for $Q$ qubits there is $2^Q$ possible measurement outcomes, and the probability of any one outcome being observed will - for many circuits - converge to zero as $Q$ increases. Then, even for very large $T$, there could be, e.g., at most 2 counts for any one measurement outcome, and the power spectra obtained will be useless.\n",
"\n",
"However, it is likely in many circumstances that any drift will show up in the marginalized probabilities for obtaining 0 or 1 of each qubit. E.g., in the example above this is precisely what is drifting. We can implement an analysis like this by specifying that the analysis function marginalizes.\n",
"\n",
"To do this in an automated way, we specify ` marginalize = 'std'`. This then assumes that the data is such that the 0th indexed outcomes is the bitstring 0...0000, the 1st indexed outcome is the bitstring 0...0001, etc. It then associates the first qubit with entity 0 (where first qubit has its measurement outcome recorded as the first bit in the string), the second qubit with entity 1 and so on.\n",
"\n",
"Note that, when using this method, the input data array grows exponentially in the number of qubits. As such, pre-marginalizing will be preferable for a large number of qubits, in which case the data for different qubits should be stored with different \"entity\" indices.\n",
"\n",
"Below we demonstrate this method on the same data as analyzed above"
]
},
{
"cell_type": "code",
"execution_count": 30,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# results_multiqubit_marg = drift.do_basic_drift_characterization(data_multiqubit, outcomes=outcomes, \n",
"# marginalize = 'std', confidence=0.99)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"As we can see from the returned data shape, the data has been split into separate sets for 2 different entities, each have two possible measurement outcomes."
]
},
{
"cell_type": "code",
"execution_count": 31,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# np.shape(results_multiqubit_marg.data)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"As always, we can plot the global power spectrum (averaged over everything, including entity, i.e., qubit)."
]
},
{
"cell_type": "code",
"execution_count": 32,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# results_multiqubit_marg.plot_power_spectrum()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"However, now we can look at the power spectrum for each individual qubit. This then shows us that each qubit has a different drift frequency. Something that was not obvious from the analysis implemented without marginalization. Moreover, these spectra show higher peaks (they are approximately twice as high), because -- earlier -- when we constructed the global power spectrum we were averaging pure noise with \"signal\" at both drift frequencies."
]
},
{
"cell_type": "code",
"execution_count": 33,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# results_multiqubit_marg.plot_power_spectrum(sequence=0,entity=1)"
]
},
{
"cell_type": "code",
"execution_count": 34,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# results_multiqubit_marg.plot_power_spectrum(sequence=0,entity=0)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"With this analysis we can isolate the drift frequencies of each qubit"
]
},
{
"cell_type": "code",
"execution_count": 35,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# # Drift frequencies for the first qubit\n",
"# print(results_multiqubit_marg.pe_drift_frequencies[0])\n",
"# # Drift frequencies for the second qubit\n",
"# print(results_multiqubit_marg.pe_drift_frequencies[1])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Finally, we can again plot the estimated drift probabilities. Again, these are more informative that what was obtained without marginalization. In particular, the oscillations are noticably stronger, as they are each the sum of two of the drifting probabilities from the \"raw\" unmarginalized case."
]
},
{
"cell_type": "code",
"execution_count": 36,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# # Creates an array of the true probability.\n",
"# parray_multiqubit_marg = np.zeros((1,2,2,T),float)\n",
"# parray_multiqubit_marg[0,0,0,:] = np.array([pt00(t)+pt01(t) for t in range(0,T)])\n",
"# parray_multiqubit_marg[0,0,1,:] = np.array([pt10(t)+pt11(t) for t in range(0,T)])\n",
"# parray_multiqubit_marg[0,1,0,:] = np.array([pt00(t)+pt10(t) for t in range(0,T)])\n",
"# parray_multiqubit_marg[0,1,1,:] = np.array([pt01(t)+pt11(t) for t in range(0,T)])\n",
"\n",
"# results_multiqubit_marg.plot_estimated_probability(sequence=0,entity=0,outcome=0,parray=parray_multiqubit_marg)"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"#### Example 2.3 : Marginalizing over qubits loses information and can miss detectable drift\n",
"Marginalizing involves throwing away information about correlations between measurement outcomes. As such, it is insensitive to certain types of drift. In particular, if only the correlations between measurement outcomes are drifting. This is perhaps a rather contrived situation, but it is possible. For example, consider the drifting probabilities:\n",
"$$ p_{00}(t) = 0.25-0.05 \\cos(0.05*t) $$\n",
"$$ p_{01}(t) = 0.25+0.05 \\cos(0.05*t) $$\n",
"$$ p_{10}(t) = 0.25+0.05 \\cos(0.05*t) $$\n",
"$$ p_{11}(t) = 0.25-0.05 \\cos(0.05*t) $$\n",
"In this case, as demonstrated below, the analysis of the full data clearly demonstrates drift, yet the analysis of the marginalized data does not."
]
},
{
"cell_type": "code",
"execution_count": 37,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# N = 10 # Counts per timestep\n",
"# T = 1000 # Number of timesteps\n",
"\n",
"# outcomes = ['00','01','10','11']\n",
"\n",
"# def pt_correlated00(t): return 0.25-0.05*np.cos(0.05*t)\n",
"# def pt_correlated01(t): return 0.25+0.05*np.cos(0.05*t)\n",
"# def pt_correlated10(t): return 0.25+0.05*np.cos(0.05*t)\n",
"# def pt_correlated11(t): return 0.25-0.05*np.cos(0.05*t)\n",
"\n",
"# data_1seq_multiqubit = np.zeros((1,1,4,T),float)\n",
"# for t in range(0,T):\n",
"# pvec = [pt_correlated00(t),pt_correlated01(t),pt_correlated10(t),pt_correlated11(t)]\n",
"# data_1seq_multiqubit[0,0,:,t] = multinomial(N,pvec)\n",
" \n",
"# results_correlatedrift_marg = drift.do_basic_drift_characterization(data_1seq_multiqubit,\n",
"# outcomes=outcomes, marginalize = 'std')\n",
"\n",
"# results_correlatedrift_full = drift.do_basic_drift_characterization(data_1seq_multiqubit, \n",
"# outcomes=outcomes, marginalize = 'none')"
]
},
{
"cell_type": "code",
"execution_count": 38,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# results_correlatedrift_marg.plot_power_spectrum()"
]
},
{
"cell_type": "code",
"execution_count": 39,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"# results_correlatedrift_full.plot_power_spectrum()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"What this demonstrates is that it is important to carefully decide how to analyze data from a given experiment, with physical motivations for the analysis used."
]
}
],
"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.7.0"
}
},
"nbformat": 4,
"nbformat_minor": 2
}