{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Simulating Clifford randomized benchmarking using a generic noise model\n", "\n", "This tutorial demonstrates shows how to simulate Clifford RB sequences using arbitrary $n$-qubit process matrices. In this example $n=2$." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "from __future__ import print_function #python 2 & 3 compatibility\n", "\n", "import pygsti\n", "from pygsti.extras import rb" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Get some CRB circuits\n", "\n", "First, we follow the [Clifford RB](../CliffordRB.ipynb) tutorial to generate a set of sequences. If you want to perform Direct RB instead, just replace this cell with the contents of the [Direct RB](../DirectRB.ipynb) tutorial up until the point where it creates `circuitlist`:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "- Sampling 10 circuits at CRB length 0 (1 of 6 lengths)\n", " - Number of circuits sampled = 1,2,3,4,5,6,7,8,9,10,\n", "- Sampling 10 circuits at CRB length 1 (2 of 6 lengths)\n", " - Number of circuits sampled = 1,2,3,4,5,6,7,8,9,10,\n", "- Sampling 10 circuits at CRB length 2 (3 of 6 lengths)\n", " - Number of circuits sampled = 1,2,3,4,5,6,7,8,9,10,\n", "- Sampling 10 circuits at CRB length 4 (4 of 6 lengths)\n", " - Number of circuits sampled = 1,2,3,4,5,6,7,8,9,10,\n", "- Sampling 10 circuits at CRB length 8 (5 of 6 lengths)\n", " - Number of circuits sampled = 1,2,3,4,5,6,7,8,9,10,\n", "- Sampling 10 circuits at CRB length 16 (6 of 6 lengths)\n", " - Number of circuits sampled = 1,2,3,4,5,6,7,8,9,10,\n" ] } ], "source": [ "#Specify the device to be benchmarked - in this case 2 qubits\n", "nQubits = 2\n", "qubit_labels = [0,1] \n", "gate_names = ['Gx', 'Gy','Gcphase'] \n", "availability = {'Gcphase':[(0,1)]}\n", "pspec = pygsti.obj.ProcessorSpec(nQubits, gate_names, availability=availability, \n", " qubit_labels=qubit_labels)\n", "\n", "#Specify RB parameters (k = number of repetitions at each length)\n", "lengths = [0,1,2,4,8,16]\n", "k = 10\n", "subsetQs = [0,1]\n", "randomizeout = False # ==> all circuits have the *same* ideal outcome (the all-zeros bitstring)\n", "\n", "#Generate clifford RB circuits\n", "exp_dict = rb.sample.clifford_rb_experiment(pspec, lengths, k, subsetQs=subsetQs, randomizeout=randomizeout)\n", "\n", "#Collect all the circuits into one list:\n", "circuitlist = [exp_dict['circuits'][m,i] for m in lengths for i in range(k)]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Create a model to simulate these circuits\n", "Now we need to create a model that can simulate circuits like this. Two things to note:\n", "\n", "1. RB circuits use our \"multi-qubit\" gate naming, so you have gates like `Gx:0` and `Gcphase:0:1`.\n", "2. RB circuits do gates in parallel (this only matters for >1 qubits), so you have layers like `[Gy:0Gy:1]`\n", "\n", "In this example, we'll make a model with $n$-qubit process matrices, so this will be practically limited to small $n$. We construct a model based on our standard 2-qubit X, Y, and CPHASE model, since this \n", "has all the appropriate gates. To get a model with the multi-qubit labels, we'll morph our \"standard\"\n", "module into a \"standard multi-qubit\" module, which has the same stuff in it but uses the multi-qubit\n", "naming conventions, to match point 1. above. If you can't start with a standard model, then you'll need to create an `ExplicitOpModel` object of the appropriate dimension (see the [explicit models tutorial](../../objects/ExplicitModel.ipynb)) and assign to it gates with are, for instance `('Gx',0)` rather than just `'Gx'`.\n", "\n", "Here's how we do the morphing:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "from pygsti.construction import std2Q_XYCPHASE\n", "smq2Q_XYCPHASE = pygsti.construction.stdmodule_to_smqmodule(std2Q_XYCPHASE)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now the `smq2Q_XYCPHASE` module acts just like the usual `std2Q_XYCPHASE` one but with multi-qubit conventions. We'll depolarize the target model and set one of the process matrices to a custom value as a demonstration. Here is where you can set any 2-qubit process matrices you want to any of the gates:" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "odict_keys([Label[Gx:1], Label[Gy:1], Label[Gx:0], Label[Gy:0], Label[Gcphase:0:1]])" ] }, "execution_count": 4, "metadata": {}, "output_type": "execute_result" } ], "source": [ "myModel = smq2Q_XYCPHASE.target_model().depolarize(op_noise=0.01, spam_noise=0.01)\n", "myModel[('Gx',0)] = np.kron( \n", " np.array([[1, 0, 0, 0],\n", " [0, 0.85, 0, 0],\n", " [0, 0, 0, -0.85],\n", " [0, 0, 0.85, 0]], 'd'),\n", " np.array([[1, 0, 0, 0],\n", " [0, 0.95, 0, 0],\n", " [0, 0, 0.95, 0],\n", " [0, 0, 0, 0.95]], 'd'))\n", "#print(myModel[('Gx',0)])\n", "myModel.operations.keys() #voila! you have gates like \"Gx:0\" rather than \"Gxi\"" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Since, `ExplicitOpModel` don't know how to automatically simulate multiple gates in parallel (you'd need to add an operation for each layer explicitly), we'll just *serialize* the circuits so they don't contain any parallel gates. This addresses point 2) above. Then we can simulate our circuits using our `ExplicitOpModel`, creating a `DataSet`." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Qubit 0 ---|Gy|-| |-| |-| |-|C1|-| |-|Gy|---\n", "Qubit 1 ---| |-|Gx|-|Gy|-|Gy|-|C0|-|Gy|-| |---\n", "\n", "Qubit 0 ---|Gy|-| |-| |-|C1|-|Gy|---\n", "Qubit 1 ---|Gx|-|Gy|-|Gy|-|C0|-|Gy|---\n", "\n" ] } ], "source": [ "serial_circuits = [c.serialize() for c in circuitlist]\n", "ds = pygsti.construction.generate_fake_data(myModel, serial_circuits, 100, seed=1234)\n", "\n", "#See how the DataSet contains serialized circuits (just printing the first several layers for clarity)\n", "print(ds.keys()[10][0:7])\n", "print(circuitlist[10][0:5])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Running RB on the simulated `DataSet`\n", "To run an RB analysis, we need to collect the \"success\" outcome counts for each circuit. We just build up parallel lists as show below, using the `'idealout'` key (even though it's always the same in this case) of `exp_dict` to determine what outcome we were supposed to get. These parallel lists can be used to initialize a `RBSummaryDataSet` object." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "lengthslist = []\n", "counts = []\n", "scounts = []\n", "cdepths = []\n", "c2Qgatecnts = []\n", "\n", "for m in lengths:\n", " for i in range(k):\n", " c = exp_dict['circuits'][m,i]\n", " ideal_outcome = ''.join(map(str,exp_dict['idealout'][m,i]))\n", " serial_c = c.serialize()\n", " datarow = ds[serial_c]\n", " lengthslist.append(m)\n", " counts.append(datarow.total)\n", " scounts.append(datarow[ideal_outcome])\n", " cdepths.append(c.depth())\n", " c2Qgatecnts.append(c.twoQgate_count())\n", "\n", "#Then we can create a RBSummaryDataset object, which contains all the info needed to perform \n", "# our standard RB analyses:\n", "data = rb.results.RBSummaryDataset(nQubits, lengthslist, success_counts=scounts, total_counts=counts,\n", " circuit_depths=cdepths, circuit_twoQgate_counts=c2Qgatecnts)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "One we have a `RBSummaryDataSet`, we can just follow the final steps of the [RB analysis tutorial](../RBAnalysis.ipynb) to perform a standard RB analysis and plot the results." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "#Do a standard RB analysis\n", "rbresults = rb.analysis.std_practice_analysis(data)" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfUAAAFECAYAAAA+3qfBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJzs3XecVOXZ+P/PNbM725elSVvKKr1LFSwBKzYsGEVT1OhjEuMTH02MJYkak3zVX4z18bEkRowm2DVECUQFKxhExQLSQVkQhaUuW6Zdvz9mdp3dnZ0507ZxvV+v85qZM/e5z70scJ27i6pijDHGmPbP1doFMMYYY0x6WFA3xhhjOggL6sYYY0wHYUHdGGOM6SAsqBtjjDEdhAV1Y4wxpoOwoG6MMcZ0EBbUjTHGmA7CgroxxhjTQWS1dgES1a1bNx0wYEBrF8MYY4xpEe+///5OVe3uJG27C+oDBgxg+fLlrV0MY4wxpkWIyOdO02a0+V1EZojIGhFZLyLXRfm+v4i8JiIfi8jrIlKayfIYY4wxHVnGgrqIuIH7gZOB4cD5IjK8UbI7gL+q6mjgFuDWTJXHGGOM6egyWVOfBKxX1Y2q6gWeBM5olGY4sCj8fnGU740xxhjjUCb71PsAWyI+lwOTG6X5CDgbuAc4CygSka6qWpHBchlj2jGfz0d5eTk1NTWtXRRj0io3N5fS0lKys7OTzqO1B8r9HPhfEbkIeBPYCgQaJxKRy4DLAPr169eS5TPGtDHl5eUUFRUxYMAARKS1i2NMWqgqFRUVlJeXU1ZWlnQ+mWx+3wr0jfhcGj5XT1W3qerZqno48MvwuT2NM1LVh1V1gqpO6N7d0ah+Y0wHVVNTQ9euXS2gmw5FROjatWvKLVCZDOrvAYNEpExEPMBsYF5kAhHpJiJ1Zbge+EsGy2OM6SAsoJuOKB1/rzMW1FXVD1wBLAQ+A55W1ZUicouIzAwnmwasEZG1QA/g95kqjzHGGNPRZXSeuqrOV9XBqnqYqv4+fO5GVZ0Xfv+sqg4Kp7lUVWszWR5jjEmXF198ERFh9erVrV0UE+Giiy7i2WefdZz+9ddf57TTTov63SmnnMKePaEe4cLCQgC2bdvGOeecA8CKFSuYP39+ffp58+Zx2223JVv0tLC1340xHdqLH27lyNsWUXbdyxx52yJe/HBr/IscmDt3LkcddRRz585NS35+vz8t+RwMWurPav78+ZSUlDQ417t37/qHhsZBfebMmVx3XZN11lqUBXUHtlZuZVXFqibHAd+B1i6aMSaGFz/cyvXPf8LWPdUosHVPNdc//0nKgb2yspK3336bRx55hCeffLL+/OzZs3n55ZfrP9fVGgOBANdccw0TJ05k9OjRPPTQQ0Colnj00Uczc+ZMhg8Prc115plnMn78eEaMGMHDDz9cn9cjjzzC4MGDmTRpEv/1X//FFVdcAcCOHTuYNWsWEydOZOLEibzzzjtNyrty5UomTZrE2LFjGT16NOvWrWPz5s2MHDmyPs0dd9zBzTffDMD69es5/vjjGTNmDOPGjWPDhg0A3H777YwaNYoxY8bUB68NGzYwY8YMxo8fz9FHH13fcvHMM88wcuRIxowZwzHHHNNsORorLCzkqquuYsSIERx33HHs2LEDgGnTpvE///M/TJgwgXvuuYfNmzdz7LHHMnr0aI477ji++OKL+jxeffVVJkyYwODBg3nppZcA2Lx5M0cffTTjxo1j3LhxLFmypD79vn37OPXUUxkyZAg/+tGPCAaDQGhZ8p07dzYoX92fm9fr5cYbb+Spp55i7NixPPXUU8yZMyfu7+WNN95g7NixjB07lsMPP5z9+/c3+TNIiaq2q2P8+PHa0j7Z8Yku/mJxk2Nn1c4WL4sxB7tVq1Y5Tjv11te0/7UvNTmm3vpaSmV44okn9Ac/+IGqqk6ZMkWXL1+uqqrPP/+8fv/731dV1draWi0tLdWqqip96KGH9Le//a2qqtbU1Oj48eN148aNunjxYs3Pz9eNGzfW511RUaGqqlVVVTpixAjduXOnbt26Vfv3768VFRXq9Xr1qKOO0p/85Ceqqnr++efrW2+9paqqn3/+uQ4dOrRJea+44gp94okn6stVVVWlmzZt0hEjRtSn+cMf/qA33XSTqqpOmjRJn3/+eVVVra6u1gMHDuj8+fN1ypQpeuDAgQblPPbYY3Xt2rWqqvruu+/q9OnTVVV15MiRWl5erqqqu3fvbrYcjQH1aX7zm9/U/5zf+ta39Mc//nF9utNOO03nzJmjqqqPPPKInnHGGaqqeuGFF+pJJ52kgUBA165dq3369Kn/Gaqrq1VVde3atVoXSxYvXqw5OTm6YcMG9fv9evzxx+szzzyjqqr9+/fXHTt2qKpqQUGBqmqDP7dHH320vnyNPzf3eznttNP07bffVlXV/fv3q8/na/DzR/v7DSxXhzGyteeptwsBbTJ1HoCgBlu4JMaYRGzbU53Qeafmzp3LlVdeCYRq53PnzmX8+PGcfPLJXHnlldTW1rJgwQKOOeYY8vLy+Pe//83HH39c32y7d+9e1q1bh8fjYdKkSQ3mJd9777288MILAGzZsoV169axfft2vvWtb9GlSxcAvv3tb7N27VogVCtdtWpV/fX79u2jsrKyvg8YYMqUKfz+97+nvLycs88+m0GDBjX7s+3fv5+tW7dy1llnAaEFUeruc/HFF5Ofnw9Aly5dqKysZMmSJXz729+uv762NjQ06sgjj+Siiy7i3HPP5eyzz3ZcDpfLxXnnnQfAd7/73fprgfrzAEuXLuX5558H4Hvf+x6/+MUv6r8799xzcblcDBo0iEMPPZTVq1dTVlbGFVdcwYoVK3C73fV/fgCTJk3i0EMPBeD888/n7bffru83T1Zzv5cjjzySq6++mu985zucffbZlJamd8sTC+oONBe8mwv2ydjn3Ue1L/QfjUtcdM+3+fjGpKp3SR5bowTw3iV5See5a9cuFi1axCeffIKIEAgEEBH+8Ic/kJuby7Rp01i4cCFPPfUUs2fPBkItovfddx8nnXRSg7xef/11CgoKGnx+9dVXWbp0Kfn5+UybNi3uvOVgMMi7775bH3yjueCCC5g8eTIvv/wyp5xyCg899BCDBw+ub2YGkpofHQwGKSkpYcWKFU2+e/DBB/nPf/7Dyy+/zPjx43n//fejluPYY4+NeY/IaV6Rf1ZOr6n7fNddd9GjRw8++ugjgsFggz+vaOlT1dzv5brrruPUU09l/vz5HHnkkSxcuJChQ4emfL861qfugD8YfVBGOoP6Vwe+4rNdn/HZrs9Yu3tt/AuMMXFdc9IQ8rLdDc7lZbu55qQhSef57LPP8r3vfY/PP/+czZs3s2XLFsrKynjrrbeAUG3y0Ucf5a233mLGjBkAnHTSSTzwwAP4fD4A1q5dy4EDTcfk7N27l86dO5Ofn8/q1at59913AZg4cSJvvPEGu3fvxu/389xzz9Vfc+KJJ3LffffVf44WYDdu3Mihhx7KT3/6U8444ww+/vhjevTowddff01FRQW1tbX1fc9FRUWUlpby4osvAqGad1VVFSeccAKPPvooVVVVQOjhpri4mLKyMp555hkg9PDy0UcfAaG+9smTJ3PLLbfQvXt3tmzZErUcjQWDwfoWjb///e8cddRRUX8PU6dOrR/P8Le//Y2jjz66/rtnnnmGYDDIhg0b2LhxI0OGDGHv3r306tULl8vF448/TiDwzf/fy5YtY9OmTQSDQZ566qlm79lYUVFRs33izf1eNmzYwKhRo7j22muZOHFi2mdPWFB3oLmaeuRTbjoJtrCGMelw5uF9uPXsUfQpyUOAPiV53Hr2KM48vE/Sec6dO7e+abrOrFmz6kfBn3jiibzxxhscf/zxeDweAC699FKGDx/OuHHjGDlyJD/84Q+jjuCeMWMGfr+fYcOGcd1113HEEUcA0KdPH2644QYmTZrEkUceyYABA+jUqRMQaq5fvnw5o0ePZvjw4Tz44INN8n366acZOXIkY8eO5dNPP+X73/8+2dnZ3HjjjUyaNIkTTjihQW3x8ccf595772X06NFMnTqV7du3M2PGDGbOnMmECRMYO3Ysd9xxBxAKqI888ghjxoxhxIgR/OMf/wDgmmuuYdSoUYwcOZKpU6cyZsyYqOVorKCggGXLljFy5EgWLVrEjTfeGPX3cN999/Hoo48yevRoHn/8ce6555767/r168ekSZM4+eSTefDBB8nNzeXyyy/nscceY8yYMaxevbpBrX/ixIlcccUVDBs2jLKysia/3+ZMnz6dVatW1Q+Ui9Tc7+Xuu+9m5MiRjB49muzsbE4++WRH93JKQn3w7ceECRN0+fLlLXrPJduW4A14m5wfUDyAAZ0GpOUe63evp7yyHACP28PU3lPTkq8xHc1nn33GsGHDWrsYLa6un9zv93PWWWfxgx/8wHHwaU8KCwuprKxs7WK0mmh/v0XkfVWd4OR6q6k70GxNPY0D5WzZS2NMLDfffDNjx45l5MiRlJWVceaZZ7Z2kUwbZAPlHGiJgXKRTe7W/G6MaayuubujO5hr6elgNfU4ghpskaAeGcet1m6MMSYZFtTjiBW4raZujDGmLbGgHkesfvN0jn63oG6MMSZVFtTjCARbqKYe0eRuze/GGGOSYUE9jpg1dVsm1piDktvtrt+UY+zYsWzevJnly5fz05/+FAitDBe5YUhLqK2t5bzzzmPgwIFMnjyZzZs3N0mzZcsWpk+fzvDhwxkxYkSDud0rVqzgiCOOYOzYsUyYMIFly5YlXZbbbruNv/3tbzHTbNq0icmTJzNw4EDOO+88vN6m04YBbr31VgYOHMiQIUNYuHBh/fk9e/ZwzjnnMHToUIYNG8bSpUuTLm9HYkE9DutTN8Y0lpeXx4oVK+qPAQMGMGHCBO69914guaCe6naijzzyCJ07d2b9+vVcddVVXHvttU3SZGVl8cc//pFVq1bx7rvvcv/999evT/6LX/yCm266iRUrVnDLLbc0WEs9mptvvpk5c+ZE/W7hwoWceOKJMa+/9tprueqqq1i/fj2dO3fmkUceaZJm1apVPPnkk6xcuZIFCxZw+eWX168Ed+WVVzJjxgxWr17NRx99dFCuXRCNBfU4WiWoW/O7Me3O66+/zmmnncbmzZt58MEHueuuuxg7dmz98rHRzJkzh5kzZ3Lsscdy3HHHpXT/f/zjH1x44YUAnHPOObz22ms0XlysV69ejBs3DggtcTps2DC2bg1tQysi7Nu3DwgtV9u7d++kyrFv3z68Xi/duze/f4WqsmjRovpNUy688ML6ZWkb/0yzZ88mJyeHsrIyBg4cyLJly9i7dy9vvvkml1xyCQAej6fJvucHK5unHkeLNb9HTmmzmroxzvzrOtj+SXrz7DkKTr4tZpLq6mrGjh0LQFlZWf2uahDag/tHP/oRhYWF/PznP497uw8++ICPP/64fge2SEcffXTUtcXvuOMOjj/++Abntm7dSt++fYFQjbxTp05UVFTQrVu3qPfdvHkzH374IZMnTwZCy5eedNJJ/PznPycYDCbdffDqq6/GfUCpqKigpKSErKxQCCotLa1/uGj8M9UtlRuZLi8vj+7du3PxxRfz0UcfMX78eO655x7HG750ZBbU44i1jG46l9i1QG5M+1HX/J4OJ5xwQtSADsSs5aeisrKSWbNmcffdd1NcXAzAAw88wF133cWsWbN4+umnueSSS3j11VcbXPfJJ5/wve99D4Dt27fj8Xi4++67AXjttdfo2rUrCxYs4OKLL85Iuev4/X4++OAD7rvvPiZPnsyVV17Jbbfdxm9/+9uM3rc9sKAeR0vV1G30uzFJiFOjbg9i1S4Tqan36dOHLVu2UFpait/vZ+/evXTt2rXJtT6fj1mzZtXv513nscceqx849+1vf5tLL720ybWjRo2qf5i5+eabGTBgABdddFGDNMuWLeOBBx4gEAgwfvx4AGbOnMktt9xSn6Zr167s2bMHv99PVlYW5eXl9OnTdJOdup+pTl260tJSSktL61sZzjnnHG67rf3/XUgHC+pxKDFq6jG+S5QNlDOm4ygqKqrvnwZ44YUXWLZsGbfeemtC+SRSU585cyaPPfYYU6ZM4dlnn+XYY49tUkFQVS655BKGDRvG1Vdf3eC73r1788YbbzBt2jQWLVrEoEGDEiorwMqVKxk6dChud2i72+ZaM0SE6dOn8+yzzzJ79mwee+wxzjjjjKg/0wUXXMDVV1/Ntm3bWLduHZMmTcLtdtO3b1/WrFnDkCFDeO211xg+fHjC5e2IbKBcHPFq6ulqgrdAbkzHcfrpp/PCCy/UD5TbsGFDfTN3plxyySVUVFQwcOBA7rzzzvqa67Zt2zjllFMAeOedd3j88cdZtGhR/XS8+fPnA/CnP/2Jn/3sZ4wZM4YbbriBhx9+OOEy/Otf/6rfQz6e22+/nTvvvJOBAwdSUVFRP+ht3rx59dutjhgxgnPPPZfhw4czY8YM7r///voHhvvuu4/vfOc7jB49mhUrVnDDDTckXN6OyLZejWPL/i1s2LOh2e+P7nM0bpc75ft8Wfkla3avAaBzbmfGdB+Tcp7GdETtcevV7373u9x1110xR4R3BCeccAJ//etf6dWrV2sXpd1KdetVa36PI95DT5AgblIP6g361K3WbkyH8sQTT7R2EVrEK6+80tpFOOhltPldRGaIyBoRWS8i10X5vp+ILBaRD0XkYxE5JZPlSUa8wXCZaOmwoG6MMSYZGQvqIuIG7gdOBoYD54tI45EMvwKeVtXDgdnA/2WqPMkKEjuop2sEvC0+Y4xz7a3b0Bgn0vH3OpM19UnAelXdqKpe4Emg8fBGBepGj3QCtmWwPEmJ2/yerqBugdwYR3Jzc6moqLDAbjoUVaWiooLc3NyU8slkn3ofYEvE53JgcqM0NwP/FpH/BgqA42lj4i0Fm5GaujW/G9Os0tJSysvL2bFjR2sXxZi0ys3NpbS0NKU8Wnug3PnAHFX9o4hMAR4XkZGqDSOliFwGXAbQr1+/Fi2gk4Fy6dCgpm4x3ZhmZWdnU1ZW1trFMKZNymTz+1agb8Tn0vC5SJcATwOo6lIgF2iyULGqPqyqE1R1QktPCbGBcsYYY9qLTAb194BBIlImIh5CA+HmNUrzBXAcgIgMIxTU21SbWksNlItkQd0YY0wyMhbUVdUPXAEsBD4jNMp9pYjcIiIzw8l+BvyXiHwEzAUu0jY2+iVecdK1/aqNfjfGGJOqjPapq+p8YH6jczdGvF8FHJnJMqQqXtBO2zKxtviMMcaYFNna73HEa17PRE3dGGOMSYYF9ThaaqCcBXVjjDGpsqAeR4vV1G0/dWOMMSmyoB5HSy0+E8lq7cYYY5JhQT2OeEHbgroxxpi2woJ6HC0V1K353RhjTKosqMfRYkHd1n43xhiTIgvqcbTKlDaL6cYYY5JgQT0GJ7Vw23rVGGNMW2FBPQYntfC07dJmze/GGGNSZEE9Bkc19aCNfjfGGNM2WFCPwUlQt8VnjDHGtBUW1GNwEtQV20/dGGNM22BBPQYntXB/0J+We1mfujHGmFRZUI/BScBOV/N7pEzU/o0xxnR8FtRjcBLU01VTj3w4yMSDgjHGmI7PgnoMvqAvbhp/0J+W7VcbBPWgBXVjjDGJs6Aeg9NaeDpq65FT46ymbowxJhkW1GNwGqyd1OjjiQzkmdj5zRhjTMdnQT0GxzV1Tb2mbkHdGGNMqiyox+C0Bp6W5ne15ndjjDGpsaAeg+Pm90B6m98tqBtjjEmGBfUYHNfU09D83mCgnI1+N8YYkwQL6jG05Oj3yAcD61M3xhiTjIwGdRGZISJrRGS9iFwX5fu7RGRF+FgrInsyWZ5EteiUNutTN8YYk6KsTGUsIm7gfuAEoBx4T0TmqeqqujSqelVE+v8GDs9UeRIV1KDj4JqOoG596sYYY1KVyZr6JGC9qm5UVS/wJHBGjPTnA3MzWJ6EJNKvnY556pE1dWt+N8YYk4xMBvU+wJaIz+Xhc02ISH+gDFiUwfIkJJFAnZaaerDhPPV0LD1rjDHm4NJWBsrNBp5Vjd7uLCKXichyEVm+Y8eOFilQIkE93SvKRftsjDHGxJPJoL4V6BvxuTR8LprZxGh6V9WHVXWCqk7o3r17GovYvBp/TUbSNqdxk7s1wRtjjElUJoP6e8AgESkTEQ+hwD2vcSIRGQp0BpZmsCwJqw3UOk4b0EDKtfXGNfN0belqjDHm4JGxoK6qfuAKYCHwGfC0qq4UkVtEZGZE0tnAk9rGOpGr/dUJpa/1O38IiKbxwDyrqRtjjElUxqa0AajqfGB+o3M3Nvp8cybLkKxEauoANYEaCilM+n7Wp26MMSZVbWWgXJuTcFBPsV/d+tSNMcakyoJ6MxJtfk81qFtN3RhjTKosqEfhC/gSriknWrOPFNRgk/vZQDljjDGJsqAeRU0g8Vp3MtfUXxullp+OaXLGGGMOLhbUo0im1p1KEI72QJDKQ4IxxpiDkwX1KBLtT4fQqnLJ7oNuNXVjjDHpYEE9imTnnCfbrx4tgCfzYGGMMebgZkE9imSbvpMNxNGuqw3U2rQ2Y4wxCbGgHkWyTd/J1tSbexiwJnhjjDGJsKAeRbI19XRfZ4PljDHGJMKCeiP+oD/pOeLJ1Kx9QV+z97OaujHGmERYUG+kyl+V9LUHfAcSvibWoDwL6sYYYxIRN6iLyOkictAE/ypf8kG92l9NopvNxRpcVx2wEfDGGGOccxKszwPWicj/F977vENLJagHNZhwP3is2rjV1I0xxiQiblBX1e8ChwMbgDkislRELhORooyXrhWk0vwOiT8UxKqNW1A3xhiTCEfN6qq6D3gWeBLoBZwFfCAi/53BsrWKZPrFIyUa1GMF7liD6IwxxpjGnPSpnyEiLwCvA9nAJFU9GRgD/CyzxWtZQQ2mXDtOtKYf735WWzfGGONUloM0ZwN3qeqbkSdVtUpELslMsVpHtb8aJbGBbo0lEtRVNW4ffE2ghkIKUyqTMcaYg4OT5vftjQO6iNwOoKqvZaRUrSSVQXJ1Emm+9wa9cZeCtTXgjTHGOOUkqJ8Q5dzJ6S5IW5BqfzqEFq/xBryO0jppWrfmd2OMMU412/wuIj8GLgcOE5GPI74qAt7JdMFaQ6oj3+vz8VXhcXvipnNSC7egbowxxqlYfep/B/4F3ApcF3F+v6ruymipWkk6auoQejgooSRuOkc1dVv/3RhjjEOxgrqq6mYR+UnjL0SkS0cL7Kqatv5rpw8HTgK29akbY4xxKl5N/TTgfUABifhOgUMzWK4WVxOoSdv+5U6b8Z0E7KAGqQ3UkuPOSbVYxhhjOrhmg7qqnhZ+LUs2cxGZAdwDuIE/q+ptUdKcC9xM6EHhI1W9INn7pSLeyPcla308t9RHRaXStVCYNSWbqYOzo6Z1WlNPJJ0FdWOMMfHEGig3LtaFqvpBrO9FxA3cT2j0fDnwnojMU9VVEWkGAdcDR6rqbhE5JJHCp1Os2vWStT7mLPbiDS/uVlGpzFkcGuEeLbB7A178QT9ZruYbQmr8NY5Xi6v0VtIlt4ujtMYYYw5esZrf/xjjOwWOjZP3JGC9qm4EEJEngTOAVRFp/gu4X1V3A6jq13FLnCGxBq09t9SH1w+Hyzpmud/kV/4f4PULzy31NVtbrw3Uxgzqlb5Kx2VLJK0xxpiDV6zm9+kp5t0H2BLxuRyY3CjNYAAReYdQE/3NqrogxfsmJVb/dkVlaJW5Qa5yvpv1GkuCI5gfPKL+fHP5FWQXNPt9pTeBoJ5AWmOMMQevWM3vx6rqIhE5O9r3qvp8mu4/CJgGlAJvisgoVd3TqCyXAZcB9OvXLw23bSpWTb1roVBRqTwb+BY/cC/g2qwnedU7nqLC5ueixxsEl0jtu9pfTSAYwO1yO77GGGPMwSfWinLfCr+eHuU4zUHeW4G+EZ9Lw+cilQPzVNWnqpuAtYSCfAOq+rCqTlDVCd27d3dw68TFml42a0o2niwI4uL/+S+gv+trLsx+lVlToje9A9T6a2PeL5GgrigH/OmZQ2+MMabjitX8flP49eIk834PGCQiZYSC+Wyg8cj2F4HzgUdFpBuh5viNSd4vabWB2pjT2er6zZ9b6uPNyjEsldFcnfMC7/efgZ/ogT1WTd0X9CW8Ulylt5JiT3FC1xhjjDm4ONl6tauI3CsiH4jI+yJyj4h0jXedqvqBK4CFwGfA06q6UkRuEZGZ4WQLgQoRWQUsBq5R1Yrkf5zkOAmwUwdn88cL85nzkwKyzryEXP8B+q94pvk8Y9T8D3gTr3XbYDljjDHxONl69UngTWBW+PN3gKeA4+NdqKrzgfmNzt0Y8V6Bq8NHq0m01nygywC2Dz6OPqteYuuwk6kp7tkkTayaejIBOl1L2BpjjOm4nOzS1ktVf6uqm8LH74AemS5YS0pmffVN4y5AXW4OXf7XqN8HNdjsbm3JBPX93v2EnoGMMcaY6JwE9X+LyGwRcYWPcwk1m3cYyayv7i3oypZRZ3HIpnfotH1lQvkmM0UtqEFbB94YY0xMzQZ1EdkvIvsILRDzd8AbPp4kPL2so0h2e9MvRp9NTUE3Bi39EwQDjvINajDpLV6tX90YY0wszQZ1VS1S1eLwq0tVs8KHS1U71DDsZGvAwawcNkz+AYW7NtF7zb+bfB+tWb/KV5X0xjEW1I0xxsTiZKAcItKZ0Pzx3LpzqvpmpgrVkmL1fTuxY8BUdvcaRdn7f+PrsiPx537zvBOtpp5KYLaV5YwxxsTiZErbpYRGvy8EfhN+vTmzxWo5Nf4alBQGoImw/ohLyfIeoOyDvzf4KloLQCqB2WrqxhhjYnEyUO5KYCLweXg9+MOBPbEvaT+SGfne2IEuA9g67GR6r15IQcWmmHmnEpi9AW9KrQrGGGM6NidBvUZVawBEJEdVVwNDMluslpOuILl53AX4cgoZvPQhCPeZR1sqNtXattXWjTHGNMdJUC8XkRJCS7q+IiL/AD7PbLFaTrKD1hrz5xSyceKFdPrqM3qufQ0IrdkeObc8kT3Um5PManTGGGMODnEHyqnqWeG3N4vIYqAT0Crbo2ZCuoI6wPbhMmCiAAAgAElEQVRBx9Fz3SIOWzaHin4T8eWVoCiCAOmpZe/37U85D2OMMR2Tk5o6IjJORH4KjAbKVbXDdOymNEiuMRHWHvlj3P4aDlv2aCj/iJp6Okav23KxxhhjmuNk9PuNwGNAV6AboR3VfpXpgrWYNK+8WlXSly2jzqLn+tcp2fZxg4eGdATkKl8VgSgL3RhjjDFOaurfASaq6k3h7ViPAL6X2WK1nLTW1MM+H/ttqot6MnjJA2jEXPVkV5KLpGhaRuwbY4zpeJwE9W1ELDoD5BDaH71DSGefen2eWTmsPfLH5O/dhuutu4BQM3y61m63NeCNMcZE0+xAORG5j1Dj9F5gpYi8Ev58ArCsZYqXeZmoqQPs7jOW7YdNo8c7d8OIM/F2H5y2B4hqXzXkxU6zZteahO43pMsQXOJoiIUxxpg2Ktbo9+Xh1/eBFyLOv56x0rSCTAV1gPVHXEKP7Svhxcup/t7zacs3XjO+L+DjywNfJpRnv+J+FGQXpFIsY4wxrazZoK6qj9W9FxEPMDj8cY2q+jJdsBaTwS3K/bnFeGfcSs5zl+Ja+r8wZHpa8o3X/J7MgLwqX5UFdWOMaeecjH6fBqwD7gf+D1grIsdkuFwtJhN96pECQ0+F4WdStPR+8nd/kZY84w2US2ZAnk2VM8aY9s9JJ+ofgRNV9VuqegxwEnBXZovVcjLZ/F6f/yl3EMjOZ+hb9yJpmI5W66+N+TBS5Us8qKdjZL4xxpjW5SSoZ6vqmroPqroWyM5ckVpWxoO6KhR25/Mjf0LxjnX0/eSF+BfFyxONuq1rnQP+xGvdVlM3xpj2z0lQXy4ifxaRaeHjT3wziK7di1zxLROChGrU2wZMYceAKQz4YC6FOzemnG+smnUyNfVqf3XG/yyMMcZklpOg/mNgFfDT8LEqfK5DaImaujfgJUCQNUdeji+3iGFv3InLn9pKu80NlvMH/dQGmu4OF09Qg7aojTHGtHMxg7qIuIG/qOqdqnp2+LhLVROPGm1U3WYrmVQXgP25xaw++qcU7NlC2fK/ppRnc83vqfSNJ1PDN8YY03bEDOqqGgD6h6e0dUged2Z/tNys3Aa16t2l4ygffip9V/6TzltXJJ1vczX1VAKzBXVjjGnfnDS/bwTeEZFfi8jVdYeTzEVkhoisEZH1InJdlO8vEpEdIrIifFya6A+QqrysOEuzpcAlLjwuT5MAvHHihRzoVMrQN+8lqza5rVSbC+qpDHhLx9awxhhjWo+ToL4BeCmctijiiCncdH8/cDIwHDhfRIZHSfqUqo4NH392XPI0yXXnxk+UQt4i0iQAB7Ny+GzaVWTX7GXIW/dBEgPUavw1UQe2pRLUbQS8Mca0b7GWiQVAVX8DICLFoY/qtGo5CVivqhvD1z8JnEFooF2bkZuVwaAezjtarbqy20A2Tvw+A//zF/qsepmtI05LKO+63doatzSkEpir/FWoKiKZH2dgjDEm/ZysKDdBRD4BPgY+EZGPRGS8g7z7AFsiPpeHzzU2S0Q+FpFnRaSvo1KnUY47J2MbmcQK6gDlI2ays+9EDlv2KIU7NyScf+N8fUFfUiPf6wQ1aDvAGWNMO+Ykmv0FuFxVB6jqAOAnwKNpuv8/gQGqOhp4BXgsWiIRuUxElovI8h07dqTp1vV5k+POSWuedXLdufiCPvxBf3M3Z80xP8WX24nhi/+A25vYQLXGATgdA91sZTljjGm/nAT1gKq+VfdBVd8GmolSDWwFImvepTTah11VKyKmx/0ZiNoCoKoPq+oEVZ3QvXt3B7dOTKaa4BuPfI/Gl1vMquk/I2//Vwxe8mBC/evVvoZ5p6NP3PrVjTGm/XIS1N8QkYfCq8l9S0T+D3hdRMaJyLgY170HDBKRsvCUuNnAvMgEItIr4uNM4LNEf4B0yNRgudys3CaBN5q9PUew+fDz6bHhDXqvXuA4/8aLxaQjINsIeGOMab/iDpQDxoRfb2p0/nBCG5ceG+0iVfWLyBXAQqBuEZuVInILsFxV5wE/FZGZhGr+u4CLEv8RUpepmnqeO49d1bscpZ2bfwbHyEomvfMnfv+f3gyZOpypg2Mvsd+4FSAtNXWv1dSNMaa9cjL6PelNwFV1PjC/0bkbI95fD1yfbP7pkomaulvcZLuzHS29umStjzmv+3nSfzn/9PyK2/x3c/bi3wPdYwb2uvXa60arp6OWXe2vJhAM4Ha5U87LGGNMy8rMsO92JhM19bo8fQFf3LTPLfXh9cM+Cvmx73/ozH7ulPt4cUnsB4KgBgloaCvXan918wPyEqCoNcEbY0w7ZUGdzKwqV5enkw1jKiq/SbNKB3CD7xKmuFdxae1cx/er9KYvEFtQN8aY9smCOqH139M9Vz2RJv2uhQ0Xe3k+eAyP+U/gsqyX6bFusaM89vuSW242mnQ+IBhjjGk5Thaf+baIFIXf/0pEno8z6r1dSne/ek6W87nvs6Zk42k0uuF2vsfmTiMZ8vb/Uvz1mrh5pDMQ7/em7wHBGGNMy3Ey+v3XqvqMiBwFHA/8AXgAmJzRkrWwsk5lKa3G1liX3C6O09YNhntuqY+KSqVroTBrSj7l/a6lx7xrGPnq/+P9mXdQW9h0jn5d8346g3qVv4qgBjO20p4xxpjMcBLUA+HXU4GHVfVlEfldBsvUKrrnp39RG3DWpw6hwN54pLufbD458VeMm/cLRr7yez487TaC2U1bFGoDtXiD3rSUF0ID8Kp8VRR6CtOWpzHGmMxzUhXbKiIPAecB80Ukx+F1Jg2qSvqyavrPKdz9OcPeuAuCgSZpMtEHboPljDGm/XESnM8ltIDMSaq6B+gCXJPRUpkGdvUdz/rJP6D75+8y8D9/abKUbCbWa0/HOvLGGGNalpPm917Ay6paKyLTgNHAXzNaKtPE1hGnk7v/a/qunEdN4SGUjzqj/rtMBGDb2MUYY9ofJzX154CAiAwEHia0ScvfM1oqE9WGyRezY8AUBi77C903vRM6qc1v7ZoK24LVGGPaHydBPaiqfuBs4D5VvYZQ7d04oAnsuhaXuPjsW1ex95ChDHvjLjptXwlkplZdtwStMcaY9sNJUPeJyPnA94GXwudi7zRiMiaYlcMnJ/ySmsJDGPXK7wlsW4E3kL6R7/X30aDV1o0xpp1xEtQvBqYAv1fVTSJSBjye2WKZWPy5xXw04zf4s/PwzD2PvL3bMnIfC+rGGNO+xA3qqroKuBb4IPx5k6renumCmdhqC7vz0YzfoBpkzIKbyDmwM+33sMFyxhjTvjhZJvZ0YAWwIPx5rIjMy3TBTHzVJaWUn3EvWbWVjP7XTWRX72027ZK1Pn72WBUX3X+Anz1WxZK18XePs2ltxhjTvjhpfr8ZmATsAVDVFcChGSyTScD+bofyyQm/Irfya8YsuImsmn1N0ixZ62POYm/9bnAVlcqcxd64gd0XjB/4jTHGtB2OBsqpauMqYDAThTGJCwQD7O01gk+Pv4H8veWhwF7bcDW4uv3aI3n9ofMx89amq9cZY4xpu5wE9ZUicgHgFpFBInIfsCTD5TIO1QXe3aWH8+nxN1Cw+wvGLLixQWCP3K89UnPn6wTVnt2MMaY9cRLU/xsYAdQSWnRmL/A/mSxUR+J0Q5d05L+rdByfHnc9Bbs+Z/SCm+sDe+P92us0d75OIMo688YYY9ouJ6Pfq1T1l6o6MXz8SlVrWqJwJr7GTeS7+k1g5bG/oHDXJsb860aya/ZF3a/dkxXaxz0Wq6kbY0z74mT0+ysiUhLxubOILMxssYxT0QJvRf/JfHr89eTv2cLYl3/JtNL9XDTdU18z71ooXDTd02Sr1yZ529AJY4xpV5xs6NItvDsbAKq6W0QOyWCZOhQhdhN3qjQYvXl/V98JfHLirxn1yu8Z+/INyMm/ZeqFie0ZHwxaUDfGmPbE0drvItKv7oOI9IcMdxR3ILlZuRnL2yUuXK7mf4V7eo/moxk346ney+EvXU/envKE8rfR78YY0744Ceq/BN4WkcdF5AngTeD6zBar48hx52Q073j93vt6DGPFKb/DFfAy7qXrKP5qteP8Mz3IzxhjTHo5GSi3ABgHPAU8CYxXVUd96iIyQ0TWiMh6EbkuRrpZIqIiMsFpwduLXHfmauq5WbmOHhoqux3GB6ffji+nkDH/+jVdP/+Po/w9bk+qRTTGGNOCnAyUO4vQAjQvqepLgF9EznRwnRu4HzgZGA6cLyLDo6QrAq4EnEWadiYnK3M19Vy3s6AOUFPciw9Pu50DXfoz8rXb6LV6QdxrMtnKYIwxJv2cNL/fFLmiXHjQ3E0OrpsErFfVjarqJVTLPyNKut8CtwMdcppcJmvqOe6chAKvL68TK07+HbtKxzHknQcY8P7fIMae6RbUjTGmfXES1KOlcTJqvg+wJeJzefhcPREZB/RV1ZdjZSQil4nIchFZvmPHDge3bjsyGRidNr9HCmbn8unxN7Bt8AkMWPE0Q9+6FwlEXy7WgroxxrQvToL6chG5U0QOCx93Au+nemMRcQF3Aj+Ll1ZVH1bVCao6oXv3xKZltTa3y022K/Z88GTluHOS6vdWl5u1R/2ETePOp+e6RYyd/0s8VbubpLM+dWOMaV+cLhPrJTRQ7ilCy8X+xMF1W4G+EZ9Lw+fqFAEjgddFZDNwBDCvIw6Wy1SNN8edk/yUORE+P3x2ePW5zYyb9zOKdqxrkCSTXQfGGGPSz8no9wOqel1dTVlVr1fVAw7yfg8YJCJlIuIBZgP1+7Cr6l5V7aaqA1R1APAuMFNVlyf5s7RZmZqrnmxNPdKOsiP58LTbABdjX76BQza8Uf+d1dSNMaZ9ids3LiKLibLYjKoeG+s6VfWLyBXAQsAN/EVVV4rILcByVZ0X6/qOJBM19WxXNm6XmxxJPe/Krofy/hl/ZMSi2xn++p0UVmxm44TvZnThHGOMMennZMDbzyPe5wKzAH8zaRtQ1fnA/Ebnbmwm7TQnebZHmWjGrntQcImL/Kx8qvxVKeXny+vERzN+w8B3/0y/T56naNdmcs5/Ggrb1xgGY4w5mMUN6qraeFDcOyKyLEPl6ZAyMVc9shbdKadTykEdQN3ZrDvyx1R2PZRB7/4ZHjwKznkEBhyVct7GGGMyz8niM10ijm4ichLQqQXK1mFkovk9svZf7ClOa95fDj2Jr87/G3gK4LHT4c0/gG3uYowxbZ6T5vf3CfWpC6Fm903AJZksVEeTib7pyNp/p5z0P2PllU6CH74B//wfWPQ72PwOnP0na443xpg2zEnze1lLFKQj87g8uMQVd/OVRETW1POy8shyZeEPOhrqEJcgFHmKwOWGWX+GsqPhX9eGmuPPegAOizlG0hhjTCtptvldRCaKSM+Iz98XkX+IyL0i0qVlitcxiEjap4dF1tRFhE6e9NXWCz2FuF3uusxh/EVw6WuQWwyPnwXzrwFv6n34xhhj0itWn/pDhBadQUSOAW4D/grsBR7OfNE6lnSPgG+cX3FO+vrVoz4g9BwJP3wTjrgclj0MDx0N5R1uSQFjjGnXYgV1t6ruCr8/D3hYVZ9T1V8DAzNftI4lnYPlXOJqsvRsOmvqzT4gZOfBjFvhwn+CvxYeOSHU3+73pu3exhhjkherT90tIlmq6geOAy5zeJ2JIp3T2nLcOYhIg3NFnqK09dvHHXhXdgz8+B1YcH1oZPzaBTDzf6H32JTvXafSW0m1v9pR2i65Xb7pLjDGmINYrOA8F3hDRHYC1cBbACIykFATvElAOpvfo9X63S43xZ5i9tTuSSnvguwCZ60KuZ3gzP+DIafAS1fBn6aHmuanXQ85hSmVQVX5tOJTavzOduM9rNNh9C3uGz+hMcZ0cM02v6vq7wntoDYHOEq1fuNtF6FNXkwC0jmtrbm8uuZ2TTnvLrkJjoEcdhpc8R6MuxCW/i/83xGwdmFKZfiq6ivHAR1gy/4tBIKBlO5pjDEdQczFZ1T1XVV9IXIDF1Vdq6ofZL5oHUteVl7G8+qc2znlvBMO6gB5JXD63XDxgtCCNX8/F56+EPZvTzgrVeXzfZ8ndI036OXLA18mfC8TnaoSCAYcHemcpmmMSZ31jbeQXHdu2vq8C7ILop4v9BTicXvwBpIbuOYSV2oL2fSfAj98C965J9TXvmERHHMNTP4RZDmb0vd11deO+9IjfbH/C3oX9sYlTnYTNrFs3LuRLfu3OEorCJN7TbbNf4xpI+x/wBYiImmrredn5Tf7XVI17bCSnJLUg2KWB751DVy+FPpNgVd+HWqSX7MAtMlmfw2oKl/s/yKp23oDXrYfSLxlwDTkC/rYVrnNcXpFKd9fnsESGWMSYUG9BeVnNx+MnRIkZq0olaCeyrVNdD0MvvM0fOdZEBfMPQ+emAU71jR7SUVNBQd8B5r9Pp4v9n+BxnlwMLFtq9xGQBMbn7DtwDZ8QV+GSmSMSYQF9RYUq4btVG5WbszadOec5PvV0xrU6ww6IVRrP+nW0GI1/zcFXv45VH7dJGmifemN1fhr+Krqq5TyOJgFggG27t+a8HVBDSZ1nTEm/Syot6B0BPXm+tPrZLuzQ+u2Jyg3KzctLQlRubNhyuXw0w9CS84u/wvcMwZe+y1Uh6bg7a7ZzX7v/pRvtWX/FqutJ2l71Xa8weTGY2yt3GozEIxpAyyot6B0BE0nDwbJjIJPpYbvWEE3OO3O0BS4ISfDW3eEgvs797Bl97q03OKA7wAVNRVpyetgoqqOB8dF4wv62F5lYxqMaW0W1FtQOmrqTh4MkgnQJTklyRQnOV0Pg3P+ElpLvnQCvHIjQx4/jz6rXsblr005+y/2JTfY7mBWUVOR0NoA0VgTvDGtz4J6C3K73CmvAe9kBH2xpzjhUezpmOOesF5j4LvPsens+6kp6sGgpQ8z+ekfUvrJi7h8yQeYfd597K21RQ8TkciI9+ZU+avYXbM7DaUxxiTLgnoLS7UJ3sn1dUvGOlWQXZD2rWGdCmqQ8i79+fDUW1lxyu+oKunLwGWPcsTT/0W/Fc/gTnKL14pqa4J3qspXxa6aXfETOpCOhwNjTPJs8ZkWlp+Vz26Sq814XJ4mu7M1p3NuZ8frwLdo03sj+737Q1OoRNjTaxR7eo2i+KvP6L/iGQ59/wn6ffIC24acxNYRp1Jb0M1xvrtrrcboVDpX49tZvZPaQG1adyU0xjhnQb2FpVJTT+TaRAJ1qzS9h0V78NjXYxifnHQjhTvX0+/j5+n76YuUfvoPdhx6FFtGnkFlt8Pi5lvprcQX9Dl+CDpYBTWY1kV7FOXLyi8Z0GlA2vI0xjhnQb2FpTJYLpEV6Yo9xbjFHXchEUFSWxo2RbFaEyq7DWTVsb8gd/9X9Fn5Er3WvkKPDW+wp+dIykecTkW/iWgzW64qyt7avXTLc167PxjtqNqR9oVjth3YRv/i/k22BzbGZJ4F9RbWUjV1EaEkpyTu9K5CT2Gr1WaDGnQ0oK2mqAcbjriEzeNm02vNq5Su/CcjX7uVmoKufDnkJL4cfDzegqY71O2p2WNBPY5Ut+qNxhvwUu2vzty6B8aYZmU0qIvIDOAewA38WVVva/T9j4CfAAGgErhMVVdlskytLced46gGHU2itfyhXYbGrYVluVrvua42UJvQBjcBTwHlo85g64jT6PrFe/Re/S/KPvg7Az58kp39J7Nt6Ax29x4dWpYWktoY5mCzz7svY/laUDem5WXsf3QRcQP3AycA5cB7IjKvUdD+u6o+GE4/E7gTmJGpMrUV+dn5Sa2eluh/ktnubLLdbbdPOdkHCnW52TngCHYOOILcfV/Se/VCeq19le6bl1Jd1JPtg6bz1cDpZOX3SHOJO5ZAMECVL7nZBfHs9+6nZ0HPjORtjGleJqtpk4D1qroRQESeBM4A6oO6qkZWEwqAg2J9z/ysxIO6S1zkujvW9pZZkvpfv5riXmycdBGbx11At81L6LXuNQZ88CRlH8ylqnQCTLgUhs8M7fNuGqj0VaIZ+ieXjiV/jTGJy2RQ7wNErjtZDkxunEhEfgJcDXiAY6NlJCKXAZcB9OvXL+0FbWnJNEvmZeV1uIFHIpJ0V0RjwSwPXw+cxtcDp5Gz/2t6rn+d0g1vwos/gvk/hyGnwMiz4bBjIcumW0Hmmt4h9MAQ1KDtb29MC2v1f3Gqer+qHgZcC/yqmTQPq+oEVZ3QvXv3li1gBsTblCWajto/mYk+/dqiQ/j88HPZ/oP5cPGCUDBf92+YOxv+MAhevBzWvQqBg3u70FS2uY0nqMGMNe0bY5qXyZr6VqBvxOfS8LnmPAk8kMHytBnJTGtLx7rxbVG2K5vaQPz13pes9fHcUh8VlUrXQmHWlGymDo49XiDLnQ39p4SOU/4IG1+HlS/AZ/+EFX+DvM4w7HQYcRb0PwqyWmdVvdYiZLblx2rpxrS8TAb194BBIlJGKJjPBi6ITCAig1S1bnuuU4H0bNXVxuVm5dI5t3NCW4S25lzyTCryFFHpq4yZZslaH3MWe/H6Q58rKpU5i0NbhMYK7A22oM3ywOATQ4f/blj/Gqx8Hj59Hj74K+QUw8DjYPDJoT3g8zOwt3wbk+mZD605s8KYg1XG/tWpql9ErgAWEprS9hdVXSkitwDLVXUecIWIHA/4gN3AhZkqT1viEhdjuo9p7WK0CcWe4rjLlD631Fcf0Ot4/aHzzQV1t7ib7+bIyoGhp4QOXzVsWARr/gVrF4Zq8uKGfkfA4BmhLWK7DUrmR2vzLKgb0/Fk9F+dqs4H5jc6d2PE+yszeX/T9hXnxN94pqIyeotGc+chVEt3NLAwOw+Gnho6gkHY9kE4wC+AV34dOjqXwWHT4dDpUHYM5LXeWvnp5Jboq/E1lkzXh0tc1vxuTCuwR2nTqvKz8slyZeEP+ptN07VQogbwroXNB20nDwtNuFyh/d1LJ8Bxv4Y9X8CaBbD+VfjoKVj+l9DCNr3HfRPkSye22754JzXpZLs+0jFd0RiTOPuXZ1qViFDsKY659eesKdkNAguAJyt0vjmdPGkYg1DSDyZfFjr8Xti6HDYsho2L4a0/wpt/AE8h9J0cGozXbyr0GQ/Z7WM9ASc7qSXT9QG06UWPjOnILKibVlfkKYoZ1OuCRyJNwInsJ+9Ilgf6Tw0dx/4SqvfA5rdCI+o/XwKLfhdK5/aEAnu/KaG0fSdBbtsc5FjsKcYlrphL9SbT9QHQOaf1dv4z5mBmQd20Oifz9qcOjt+PW8fj8mS+pphXEpoON+z00OeqXfDFu/DFklCQX3IvvH0nINB9aCjQl46HPhPgkOHgbv1/em6Xm2JPccxNXZLp+gAoye0Y4w6MaW9a/38Wc9ArzC5Ma36tslBPfpdvRtQDeA9A+XuhQL/1fVgzH1Y8EfouKw96jw0F+j7joOdo6HIoNLONbCZ1zu0cM6gn0/UhhHYINMa0PAvqptXlZuXGbQZORLofEpLiKYBDp4UOAFXYvQnK3w8F+a3LYdmfoG7hnex86DECeo4KH6NDNXpPZh9Q4gXfpLo+coptOpsxrcT+5ZlW5xIXeVl5aVu2NJlleDNOJFQb73IojP526JzfCztWw/ZPvjk+eS40yh5CI+27DgwF+UOGQbchoab8LmWQpu6FYk9x3PX3E+n6AOtPN6Y1WVA3bUJBdkHagnq7WSc/ywO9RoeOOqqhqXSRgX7Le/Dpc9+kcWVD18Og2+BQkO8+JHR0HRiad58AEaFLbhd2VO9I0w8VatI3pj1RVb6q+iqhVT7rFHmKKPS0gdbBMAvqpk3Iy0osGMXSrtfJF4HO/UPHsNO+OV9bCTvXho4dq2HHWvhqJax+Ceq7LQSK+4Rq8nWtAvVHWbPbz3bL65a2oO5xe9I/88CYDPuq6itW71qd1LV5WXlM7DmxzSy2ZEHdtAnpCuoucXXMOdI5haFBdX3GNTzvq4FdG0KBfuc62LUJdm2E1S9D1c6GaQt7hgJ85/7QqS+U9IVOfelS3BN3wE8gDSPyu+V163BbBJuOzRf0sWHPhqSvr/ZXU76/nH7FbWNbcAvqpk1wshCKE+ms8bcL2bmhAXY9RjT9rmbvN0F+18bw+w2w6U3Y/2V9DT8bOBrw5pVQU9idmsJDqC3oHn7fHW9+F2rzu+DLK0HjjNDvntf+t0Y2B5fNezfjC6a2DfPmfZs5JP8QcrNaf+EpC+qmTUjXP4Zcd+v/o2ozcjuFps71Htv0u4AP9m2DvVtgzxb27/iM/TtWkVu5g8Jdm+j6xXu4A94Gl6i48OZ2qg/y37x2xpvfhUDhIZQUD4Ksgna7dK45uFR6K9lWuS3lfIIaZMOeDYzoFuXhuoVZUDdtQq47F0FQEh+o0iCfNvCk3C64s7/puwdyAl4+2Lb0mz9/VbJr9pJbuQNP9W48VbvIqdqF50DoNefATop3riO7ei8S7XeW0wkKuoWO/G6N3neHgq6h1/xukN/VHgJMq9iwd0PK/+fU2VG9g721e1t9m2wL6qZNEBHys/NTHgHfbka+tzEet4ciTxH7vPtCJ0Tw5ZXgi7MjnQT9eKr34KnaxUB3IZ18NXCgAg7sCPXpH9gBuzeH5uUf2AnNTZ3LLgit0pfXGXJLQu/rXuvfdw6/7xyRtlOrLNpj2r9qfzW7a3anNc/tB7ZbUDemTrGnOOWgnpaNXA5SJbkl3wR1h9SVRW1BN2oLupHf+8jY8+eDQajZEwruVTtDrwd2QFVFaC39mj2h1+rdoTEAde/91bELkV0AOUXhozD06mn8uRByiptP4ykI5dMGlu81LeOrA1+lPc+vq75mYMlA3K34oGl/g02bUewp5ssDXyZ9vUtcbXPhmXaiJKeEL/giqWsLsgvizzpwuULL6eZ3AQY7z9xf+02Ajwz8de9r90PtPvBWht9XwoHN4N3/zWenA6Fc2aHV/bLzQqv51b3PbvTeE3k+yvdZHsjKhawccOeE3zc+lxOawmhaxVdV6Q/qAQ2ws3onPQp6pD1vp0JW7scAABEISURBVCyomzYjqT3QI6/3FNt0qhR08nRKernejK4il5UDRT1CRzJUQw8G9UF/f8P3tfvBVwXeqtCrrxp8B8Kv1aF1/H3VoQeJ+u/Dr/6a1H62uuAeGeijPgB4QjsAurJDrQn177PBFf7s+H34uubeu7LChzv0Kq6I9+FXV9uYk52svbV7qY7XApSk7VXbLagbA6FFY+ItWRpLqg8FBzu3y02Rp4i9tXsTvrZN78omEpr6l50bGqyXTsHAN8G//kGgKrQEsL8GAuFXf234iHauufM13zxM+L2h1oaAFwL+6O9bWn2AjxP8G6Rzhz9Hpos8H04nLgeHND3X4Noo34cP9e6jX+2e0HsERNDwNUr4NXy9igAuvhxygqOWlT01zW+Q1BIsqJs2Q0QozilOevCKrWSWupKckoSDuiCtPjio1bjc4T76QqAV5+irhh4wogZ+HwT9Md6H09W9D/rCeflDaxkE/aHPGj4XDJ/TQBrShQ9/bcPPGoxxaJRzgTjfRxxASfhIxJdDT3T2q0AJarDVVpizoG7alEPyDkkqqGe5suiS2yUDJTq4JPNglJ+dT7arA67i156IhJvlsxJe//+go8raXav5snIrqCLhBwHRIKBI+KFAVMOfg+HvnAtowIK6MQDd8ruxbs+6hPt1u+d1bzNrL7dnyWxMYYMTTbsiQlAEDW8PHGuW+pK1vohth6vibjtcJxAMtNqDrv0vaNqUbFd2UjXuQ/IPyUBpDj457pyE/zNqE/vXG5OAQDD+uJ0la33MWeylojIU9isqlTmLvSxZG38mRbLjgtLBgrppc3rkJzZy1OP2UJLThgdqtTOJ1rytpm7aGyeLVD231IfX3/Cc1x86H4tLXK26XLUFddPmdMntglucL97QI7+HTWVLo0Rr3m1pL2ljnOiWF38WRF0N3en5Ol1yu3TcxWdEZAZwD+AG/qyqtzX6/mrgUsAP7AB+oKqfZ7JMpu1zu9wc0fsIVJ2tyZzIA4CJL5Egne3KTtsOe8a0lCJPETnuHGoDtc2m6VooUQN418LYFQgnDwyZlLGauoi4gfuBk4HhwPkiMrxRsg+BCao6GngW+P8yVR7TvmS7svG4PY6O1nwq7ogSqalbLd20V/GC76wp2XgaVXs9WaHzzRGErrld01G8pGWy+X0SsF5VN6qqF3gSOCMygaouVtWq8Md3gdIMlscY40BBdoHjmQRFnqIMl8aYzOieF3tdgamDs7louqe+Zt61ULhouifm6PdOOZ3iL5ecYZlsfu8DbIn4XA5MjpH+EuBf0b4QkcuAywD69euXrvIZY6IQEQqyC9jv3R83bVG2BXXTPpXkltA5t3PMdTGmDnY2ha1OWaeydBQtJW1ioJyIfBeYAPwh2veq+rCqTlDVCd27t+KqTcYcJJzWwK353bRng0oGpW19i54FPdvEyoqZDOpbgb4Rn0vD5xoQkeOBXwIzVbX5UQvGmBbjpAae5coiL8tWLzPtV352PqWFqff6usXNoZ0OTUOJUpfJ5vf3gEEiUkYomM8GLohMICKHAw8BM1T16wyWxRiTgE45neIOJMrPij/X15i2rl9xP7ZXbcebwqY4AzoNwOP2pLFUyctYUFdVv4hcASwkNKXtL6q6UkRuAZar6jxCze2FwDPhecZfqOrMTJXJGONMfnY+I7uNbO1iGJNxWa4sDut0GJ/t+iyp6wuyC+hT2CfNpUpeRuepq+p8YH6jczdGvD8+k/c3xhhj4ulR0IPtVduT2kxqcOfBbWrfibZTEmOMMaaVJBOcexf2bhOD4yJZUDfGGHPQy8vKo6zY+ZQ0j9vTZgbHRbKgbowxxgClRaWONyga3HkwWa62t3u5BXVjjDGG0MJLh5UcFjedk9khrcWCujHGGBPWOaczuVmxt07tVdCrhUqTOAvqxhhjTJiI0Lugd7PfZ7my4q4b35osqBtjjDERehb0RIi+xWrP/J5temdIC+rGGGNMBI/b02yfea/Cttv0DhbUjTHGmCaKc4qbnMtyZTkeHd9aLKgbY4wxjUTbrCjHndMKJUmMBXVjjDGmkWgBvD3sSmhB3RhjjPn/27v3WDnKMo7j3x8HSisUaGm5CaFCiFqoraUQwIagJFgqNwOGIkoRE4KUCEYkJSAhVaKIGIIS5JoAIhKUS4GirYVQJRR6oVcKpWBRK1AKWGhopbSPf8x7ynY7e86e7u7Z3Tm/T7I5szPvvPO+Z3b32Xlndp4y/Tu2/Vmbj9TNzMza0E4dO21zL3gHdTMzszZVHtRbKRtbJa3fQjMzsybYAQd1MzOzYii7/4yDupmZWZtqhyBervXyxpmZmbWAYbsNY+PmjVue79Zv2xvStBoHdTMzsxz77LJPs5vQY+03tmBmZma5HNTNzMwKwkHdzMysIBzUzczMCsJB3czMrCAaGtQljZP0sqQVkibnLD9W0nxJH0s6o5FtMTMzK7qGBXVJHcBNwInAcOAsScPLiv0TOBf4faPaYWZm1lc08nfqRwIrIuI1AEl/AE4FXuwsEBEr07LNDWyHmZlZn9DI4fdPA/8qef7vNK/HJJ0vaa6kuW+//XZdGmdmZlY0bXGhXETcGhFjImLM0KFDm90cMzOzltTIoL4KOKDk+f5pnpmZmTVAI8+pzwEOkfQZsmA+AfhmrZXOmzdvjaTXa62nBkOANU3cfm9zf4urL/UV3N8iK3pfD6y2oCKiYa2QNB64AegA7oyIayRNAeZGxFRJRwAPAYOADcCbEXFowxpUB5LmRsSYZrejt7i/xdWX+grub5H1pb52p6FZ2iJiGjCtbN5VJdNzyIblzczMrEZtcaGcmZmZdc9BvedubXYDepn7W1x9qa/g/hZZX+prlxp6Tt3MzMx6j4/UzczMCsJB3czMrCAc1CuoIsPczpLuT8ufkzSs91tZH5IOkPSUpBclLZV0cU6Z4yStlbQgPa7Kq6tdSFopaXHqy9yc5ZJ0Y9q/iySNbkY7ayXpsyX7bIGk9yVdUlamrfetpDslrZa0pGTeYEkzJL2S/g6qsO7EVOYVSRN7r9Xbr0J/r5P0UnqtPiRpjwrrdvm6bzUV+nq1pFUlr9fxFdbt8jO8sCLCj7IH2e/qXwUOAvoBC4HhZWUuBH6bpicA9ze73TX0d19gdJoeCCzP6e9xwGPNbmsd+7wSGNLF8vHAE4CAo4Dnmt3mOvS5A3gTOLBI+xY4FhgNLCmZ9wtgcpqeDFybs95g4LX0d1CaHtTs/mxnf08AdkzT1+b1Ny3r8nXfao8Kfb0auLSb9br9DC/qw0fq+bZkmIuIj4DODHOlTgXuStN/BI6XpF5sY91ExBsRMT9NfwAsYzuT7xTIqcDdkZkN7CFp32Y3qkbHA69GRDPvyFh3ETELeLdsdun78y7gtJxVvwrMiIh3I+I9YAYwrmENrZO8/kbE9Ij4OD2dTUHu/1Fh31ajms/wQnJQz1dNhrktZdKbaS2wZ6+0roHSaYQvAs/lLD5a0kJJT0hq6Tv/VSGA6ZLmSTo/Z3ndsgy2kAnAfRWWFWnfAuwdEW+k6TeBvXPKFHEfA5xHNsqUp7vXfbu4KJ1quLPCqZWi7ttuOajbFpJ2Bf4EXBIR75ctnk82bDsS+DXwcG+3r87GRsRo4ERgkqRjm92gRpLUDzgFeCBncdH27VYiG4/tE7/dlXQF8DFwb4UiRXjd3wwcDIwC3gCub25zWouDer5qMsxtKSNpR2B34J1eaV0DSNqJLKDfGxEPli+PiPcjYl2angbsJGlILzezbiJiVfq7miz/wJFlRYqWZfBEYH5EvFW+oGj7Nnmr83RJ+rs6p0yh9rGkc4GTgLPTF5ltVPG6b3kR8VZEbIqIzcBt5PehUPu2JxzU823JMJeOcCYAU8vKTAU6r5Y9A3iy0hup1aVrAe4AlkXEryqU2afzmgFJR5K9dtryS4ykXSQN7Jwmu8hoSVmxqcA56Sr4o4C1JcO57egsKgy9F2nflih9f04EHskp8xfgBEmD0hDuCWle25E0DrgMOCUiPqxQpprXfcsru7bl6+T3oZrP8GJq9pV6rfogu/p5OdkVlFekeVPI3jQA/cmGMlcAzwMHNbvNNfR1LNnw5CJgQXqMBy4ALkhlLgKWkl1FOhs4ptntrqG/B6V+LEx96ty/pf0VcFPa/4uBMc1udw393YUsSO9eMq8w+5bsy8obwEayc6ffJbu+ZSbwCvBXYHAqOwa4vWTd89J7eAXwnWb3pYb+riA7h9z5/u38Zc5+wLQ0nfu6b+VHhb7ek96Ti8gC9b7lfU3Pt/kM7wsP3ybWzMysIDz8bmZmVhAO6mZmZgXhoG5mZlYQDupmZmYF4aBuZmZWEA7qZmZmBeGgbmZmVhAO6mY1kLQp5XReIunRzjzWJfMXSpov6Zhmt7U3SFpX5/r2kHRhyfNhpbm1u1l3gKSnJXX0YHv9JM1Kt342azsO6ma1WR8RoyLiMLIUkZPK5o8ELgd+1rQWtrc9gAu7LZXvPODBiNhU7QqRpemcCZy5nds0ayoHdbP6eZb89I67Ae+Vz0z34n48Hc0vkXRmybJvSXo+He3f0nm0KekKScsl/V3SfZIuLT96TfOu7qqutM4ySbdJWippuqQBqfw5Ka3lQkn3dNemSrZjuz+W9HJp34CfAwenOq5LVXfkrZ/jbEru+S7pAUm/SfW/LmmspHvS//OOkvUeTuuatR0HdbM6SAHueD5JGjEgBaKXgNuBn+SsNg74T0SMTEf6f051fZ7sSPFLETEK2AScLelwssQUo8jua31EFe3KrSstPgS4KSIOBf4LnK4sl/qVwFfSKMPFVdRTj+0eAZwOjCTLKDcmlZ0MvJpGPX5Uaf2c7fcjy8ewsmT2COC1iBgL3EKWxOgyYDjwNUk7p3JLqOJ/a9aKfN7IrDYDJC0gO0JfBsxI89enYIako4G7JR0WWydbWAxcL+la4LGI+FuafzxwODAnJU8bQJY6dDDwUKQsXJKqyTpVqa5ZwD8iYkEqNw8YBgwCHoiINQAR8W439dRru0OARyJiA7BB0qNd1J23frkhZAEfAEn9yYbyb0izArgjUuY9SZuAj1KfN0n6SNLAiPigi3aYtRwHdbParI+IUZI+RZa2cxJwY2mBiHhWWX7yoZQEwohYLmk02VH3TyXNjIgpZBni7oqIy0vrkXRJhTZ8zNajbv1LV6tQ1zDgfyWzNpEF3kpy6+lp+e3Ybp5q1l/P1v+HQ8nyyW9Oz0cCN6c27U82YlL6hWtnYEMP22XWdB5+N6uDdPT8feCH5VdOS/oc0EFZjnJJ+wEfRsTvgOuA0WnRTOAMSXulcoMlHUh2lHtauqp7IHByKv8WsJekPdMQ8kklm6lUVyVPAt+QtGdn+e2sp6flnwFOltRf0q4lffgAGNjFerki4j2yc++dgX0EWcrRTl8gS90JWYDvnCb1fU1EbOzpds2azUfqZnUSES9IWgScxSfD8pAdtU7MuQp7BHCdpM1k+aK/l+p5UdKVwHRJO6RlkyJitqT7yYLTamBOKr9R0hTgeWAV8FJJm3LrAt6s0Ielkq4Bnk5D0i8A53ZRz+sV6unpduek0wmLyL6kLAbWRsQ7kp5RdiHgE2Q57qs1HRhLlk99RPr/dA7FD0iBH7YO8ABfBh7vwXbMWobzqZu1KWVXuK+LiF82uy31IGnXiFiXTmXMAs6PiPk11Dca+EFEfLuH6z0ITI6I5du7bbNm8ZG6mbWKWyUNJzsXflctAR0gIuZLekpSR7W/VU9XzT/sgG7tykfqZmZmBeEL5czMzArCQd3MzKwgHNTNzMwKwkHdzMysIBzUzczMCsJB3czMrCAc1M3MzAri/4J3DahiSsnEAAAAAElFTkSuQmCC\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "#Plot some stuff, etc - see the RBAnalysis.ipynb for more on this.\n", "%matplotlib inline\n", "rbresults.plot()" ] } ], "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": 1 }