{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Simulating Clifford randomized benchmarking using implicit models\n",
    "\n",
    "This tutorial demonstrates shows how to simulate Clifford RB sequences using $n$-qubit \"implicit\" models which build $n$-qubit process matrices from smaller building blocks.  This restricts the noise allowed in the $n$-qubit model; in this tutorial we take $n=3$ and use a `LocalNoiseModel`."
   ]
  },
  {
   "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\n",
    "import numpy as np"
   ]
  },
  {
   "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 = 3\n",
    "qubit_labels = list(range(nQubits)) \n",
    "gate_names = ['Gx', 'Gy','Gcphase'] \n",
    "availability = {'Gcphase':[(i,i+1) for i in range(nQubits-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 = qubit_labels\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. The RB circuits use pyGSTi's \"multi-qubit\" conventions, which mean:\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",
    "\"Implicit\" models in pyGSTi (see the [implicit model tutorial](../../objects/ImplicitModel.ipynb)) are designed to efficiently describe multi-qubit processors.  There are numerous ways of constructing implicit models, all of which can simulate the type of circuits described above.  Here we'll demonstrate the simplest type: a \"local noise model\" (class `LocalNoiseModel`) where the noise on a gate can only act on that gate's target qubits - so, for instance, 1-qubit gates are still given by 1-qubit operators, not $n$-qubit ones.\n",
    "\n",
    "The construction of a local noise model follows the same pattern as building the `ProcessorSpec` above (in fact, `pspec.models['target']` *is* essentially the same model we build below except it was built with the default `parmeterization=\"static\"` argument."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "myModel = pygsti.obj.LocalNoiseModel.build_from_parameterization(nQubits, gate_names,\n",
    "                                                                 availability=availability, \n",
    "                                                                 qubit_labels=qubit_labels,\n",
    "                                                                 parameterization=\"full\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Setting `parameterization=\"full\"` is important, as it lets us assign arbitrary numpy arrays to gates as we'll show below.  If you need to use other gates that aren't built into pyGSTi, you can use the `nonstd_gate_unitaries`\n",
    "argument of `build_from_parameterization` (see the docstring).\n",
    "\n",
    "The `build_from_parameterization` function creates a model with ideal (perfect) gates.  We'll now create a 1-qubit depolarization superoperator, and a corresponding 2-qubit one (just the tensor product of two 1-qubit ones) to add some simple noise.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "depol1Q = np.array([[1, 0,   0, 0],\n",
    "                    [0, 0.99, 0, 0],\n",
    "                    [0, 0, 0.99, 0],\n",
    "                    [0, 0, 0, 0.99]], 'd') # 1-qubit depolarizing operator\n",
    "depol2Q = np.kron(depol1Q,depol1Q)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As detailed in the [implicit model tutorial](../../objects/ImplicitModel.ipynb), the gate operations of a `LocalNoiseModel` are held in its `.operation_blks['gates']` dictionary.  We'll alter these by assigning new process matrices to each gate.  In this case, it will be just a depolarized version of the original gate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "myModel.operation_blks['gates'][\"Gx\"] = np.dot(depol1Q, myModel.operation_blks['gates'][\"Gx\"])\n",
    "myModel.operation_blks['gates'][\"Gy\"] = np.dot(depol1Q, myModel.operation_blks['gates'][\"Gy\"])  \n",
    "myModel.operation_blks['gates'][\"Gcphase\"] = np.dot(depol2Q, myModel.operation_blks['gates'][\"Gcphase\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's what the gates look like now:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "FullDenseOp with shape (4, 4)\n",
      " 1.00   0   0   0\n",
      "   0 0.99   0   0\n",
      "   0   0   0-0.99\n",
      "   0   0 0.99   0\n",
      "\n",
      "FullDenseOp with shape (4, 4)\n",
      " 1.00   0   0   0\n",
      "   0   0   0 0.99\n",
      "   0   0 0.99   0\n",
      "   0-0.99   0   0\n",
      "\n",
      "FullDenseOp with shape (16, 16)\n",
      " 1.00   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0\n",
      "   0   0   0   0   0   0   0   0   0   0   0   0   0 0.99   0   0\n",
      "   0   0   0   0   0   0   0   0   0   0   0   0   0   0 0.99   0\n",
      "   0   0   0 0.99   0   0   0   0   0   0   0   0   0   0   0   0\n",
      "   0   0   0   0   0   0   0 0.99   0   0   0   0   0   0   0   0\n",
      "   0   0   0   0   0   0   0   0   0   0 0.98   0   0   0   0   0\n",
      "   0   0   0   0   0   0   0   0   0-0.98   0   0   0   0   0   0\n",
      "   0   0   0   0 0.98   0   0   0   0   0   0   0   0   0   0   0\n",
      "   0   0   0   0   0   0   0   0   0   0   0 0.99   0   0   0   0\n",
      "   0   0   0   0   0   0-0.98   0   0   0   0   0   0   0   0   0\n",
      "   0   0   0   0   0 0.98   0   0   0   0   0   0   0   0   0   0\n",
      "   0   0   0   0   0   0   0   0 0.98   0   0   0   0   0   0   0\n",
      "   0   0   0   0   0   0   0   0   0   0   0   0 0.99   0   0   0\n",
      "   0 0.98   0   0   0   0   0   0   0   0   0   0   0   0   0   0\n",
      "   0   0 0.98   0   0   0   0   0   0   0   0   0   0   0   0   0\n",
      "   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0 0.98\n",
      "\n"
     ]
    }
   ],
   "source": [
    "print(myModel.operation_blks['gates'][\"Gx\"])\n",
    "print(myModel.operation_blks['gates'][\"Gy\"])\n",
    "print(myModel.operation_blks['gates'][\"Gcphase\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that our `Model` object is set to go, generating simulated data is easy:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "ds = pygsti.construction.generate_fake_data(myModel, circuitlist, 100, seed=1234)"
   ]
  },
  {
   "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": 8,
   "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",
    "        datarow = ds[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": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Do a standard RB analysis\n",
    "rbresults = rb.analysis.std_practice_analysis(data)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfUAAAFECAYAAAA+3qfBAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJzs3XmcVNWd///Xp7qqd5qlQfZVQGUXGhDcADWiMbjgmkmUbEYTxyRmkWQmmnGSmWSS75joz0SJW0zirlFHUWMUjYobKIIgsi/NZtOy9F7b+f1R3U13U9VdVV3V6/v5eNxHV9177rmf6kY/dc499xxzziEiIiKdn6e9AxAREZHUUFIXERHpIpTURUREuggldRERkS5CSV1ERKSLUFIXERHpIpTURUREuggldRERkS5CSV1ERKSL8LZ3AInq27evGzFiRHuHISIi0iZWrly53znXL56ynS6pjxgxghUrVrR3GCIiIm3CzLbHWzat3e9mNt/MPjGzTWa2OMrxYWa2zMw+MLPVZnZuOuMRERHpytKW1M0sA7gDOAcYB1xhZuOaFPt34FHn3InA5cDv0xWPiIhIV5fOlvoMYJNzbotzzg88DJzfpIwDCmpf9wR2pzEeERGRLi2d99QHAzsbvC8GZjYp8zPg72b2r0AecGa0iszsauBqgGHDhqU8UBHpPAKBAMXFxVRXV7d3KCIplZ2dzZAhQ/D5fEnX0d4D5a4A7nfO/T8zmwX82cwmOOfCDQs555YASwCKioq0ALxIN1ZcXEyPHj0YMWIEZtbe4YikhHOO0tJSiouLGTlyZNL1pLP7fRcwtMH7IbX7Gvoa8CiAc+4tIBvom8aYRKSTq66uprCwUAlduhQzo7CwsNU9UOlM6u8BY8xspJllEhkI90yTMjuAMwDM7AQiSb0kjTGJSBeghC5dUSr+XactqTvngsB1wIvAx0RGua81s1vMbEFtse8D3zCzD4GHgEXOOXWvi4iIJCGtz6k755Y658Y65451zv2idt9Nzrlnal+vc86d7Jyb7Jyb4pz7ezrjERFJlaeeegozY/369e0dijSwaNEiHn/88bjLv/rqq5x33nlRj5177rkcPHgQgPz8fAB2797NxRdfDMCqVatYunRpfflnnnmGX/7yl8mGnhKa+11EurSnPtjFyb98hZGLn+PkX77CUx80HdqTnIceeohTTjmFhx56KCX1BYPBlNTTHbTV72rp0qX06tWr0b5BgwbVf2lomtQXLFjA4sVHzbPWppTUm1FcVsy60nUtbhWBivYOVUSieOqDXfz4yTXsOliFA3YdrOLHT65pdWIvLy/njTfe4J577uHhhx+u33/55Zfz3HPP1b+vazWGQiF++MMfMn36dCZNmsRdd90FRFqJp556KgsWLGDcuMjcXBdccAHTpk1j/PjxLFmypL6ue+65h7FjxzJjxgy+8Y1vcN111wFQUlLCwoULmT59OtOnT+fNN988Kt61a9cyY8YMpkyZwqRJk9i4cSPbtm1jwoQJ9WV+85vf8LOf/QyATZs2ceaZZzJ58mSmTp3K5s2bAfjVr37FxIkTmTx5cn3y2rx5M/Pnz2fatGmceuqp9T0Xjz32GBMmTGDy5MmcdtppMeNoKj8/n+9973uMHz+eM844g5KSyDCrOXPm8N3vfpeioiJ+97vfsW3bNubNm8ekSZM444wz2LFjR30d//jHPygqKmLs2LE8++yzAGzbto1TTz2VqVOnMnXqVJYvX15f/vDhw3z+85/nuOOO45prriEcjjyANWLECPbv398ovrrfm9/v56abbuKRRx5hypQpPPLII9x///0t/l1ee+01pkyZwpQpUzjxxBMpKys76nfQKs65TrVNmzbNtZV397zrlu1Y1uJWXFbcZjGJdHfr1q2Lu+zs/37ZDb/x2aO22f/9cqti+Mtf/uK++tWvOuecmzVrlluxYoVzzrknn3zSXXnllc4552pqatyQIUNcZWWlu+uuu9x//ud/Ouecq66udtOmTXNbtmxxy5Ytc7m5uW7Lli31dZeWljrnnKusrHTjx493+/fvd7t27XLDhw93paWlzu/3u1NOOcV9+9vfds45d8UVV7jXX3/dOefc9u3b3fHHH39UvNddd537y1/+Uh9XZWWl27p1qxs/fnx9mV//+tfu5ptvds45N2PGDPfkk08655yrqqpyFRUVbunSpW7WrFmuoqKiUZzz5s1zGzZscM459/bbb7u5c+c655ybMGGCKy6O/L/xwIEDMeNoCqgv8x//8R/1n/P000931157bX258847z91///3OOefuueced/755zvnnLvqqqvc2Wef7UKhkNuwYYMbPHhw/Weoqqpyzjm3YcMGV5dLli1b5rKystzmzZtdMBh0Z555pnvsscecc84NHz7clZSUOOecy8vLc865Rr+3++67rz6+pu9j/V3OO+8898YbbzjnnCsrK3OBQKDR54/27xtY4eLMke39nHqHFXZhKgOVcZUt95enORoRScbug1UJ7Y/XQw89xHe+8x0g0jp/6KGHmDZtGueccw7f+c53qKmp4YUXXuC0004jJyeHv//976xevbq+2/bQoUNs3LiRzMxMZsyY0ei55Ntuu42//e1vAOzcuZONGzeyd+9eTj/9dPr06QPAJZdcwoYNG4BIq3TdunX15x8+fJjy8vL6e8AAs2bN4he/+AXFxcVcdNFFjBkzJuZnKysrY9euXVx44YVAZEKUuut85StfITc3F4A+ffpQXl7O8uXLueSSS+rPr6mpAeDkk09m0aJFXHrppVx00UVxx+HxeLjssssA+NKXvlR/LlC/H+Ctt97iySefBODLX/4yP/rRj+qPXXrppXg8HsaMGcOoUaNYv349I0eO5LrrrmPVqlVkZGTU//4AZsyYwahRowC44ooreOONN+rvmycr1t/l5JNP5oYbbuBf/uVfuOiiixgyZEirrtOUknoM5YFyHPENxC/zp7j7RERSYlCvHHZFSeCDeuUkXednn33GK6+8wpo1azAzQqEQZsavf/1rsrOzmTNnDi+++CKPPPIIl19+ORDpEb399ts5++yzG9X16quvkpeX1+j9P/7xD9566y1yc3OZM2dOi88th8Nh3n777frkG80Xv/hFZs6cyXPPPce5557LXXfdxdixY+u7mYGkno8Oh8P06tWLVatWHXXszjvv5J133uG5555j2rRprFy5Mmoc8+bNa/YaDR/zavi7ivecuve33nor/fv358MPPyQcDjf6fUUr31qx/i6LFy/m85//PEuXLuXkk0/mxRdf5Pjjj2/19eronnoMibS+K4OVhBtPgiciHcAPzz6OHF9Go305vgx+ePZxSdf5+OOP8+Uvf5nt27ezbds2du7cyciRI3n99deBSGvyvvvu4/XXX2f+/PkAnH322fzhD38gEAgAsGHDBioqjh6Lc+jQIXr37k1ubi7r16/n7bffBmD69Om89tprHDhwgGAwyBNPPFF/zuc+9zluv/32+vfREuyWLVsYNWoU119/Peeffz6rV6+mf//+fPrpp5SWllJTU1N/77lHjx4MGTKEp556Coi0vCsrKznrrLO47777qKyM9GB+9tlnFBQUMHLkSB577DEg8uXlww8/BCL32mfOnMktt9xCv3792LlzZ9Q4mgqHw/U9Gg8++CCnnHJK1L/D7Nmz68cz/PWvf+XUU0+tP/bYY48RDofZvHkzW7Zs4bjjjuPQoUMMHDgQj8fDn//8Z0KhUH35d999l61btxIOh3nkkUdiXrOpHj16xLwnHuvvsnnzZiZOnMiNN97I9OnTU/70hJJ6DIm0vsMurMFyIh3QBScO5r8vmsjgXjkYMLhXDv990UQuOHFw0nU+9NBD9V3TdRYuXFg/Cv5zn/scr732GmeeeSaZmZkAfP3rX2fcuHFMnTqVCRMm8M1vfjPqCO758+cTDAY54YQTWLx4MSeddBIAgwcP5ic/+QkzZszg5JNPZsSIEfTs2ROIdNevWLGCSZMmMW7cOO68886j6n300UeZMGECU6ZM4aOPPuLKK6/E5/Nx0003MWPGDM4666xGrcU///nP3HbbbUyaNInZs2ezd+9e5s+fz4IFCygqKmLKlCn85je/ASIJ9Z577mHy5MmMHz+ep59+GoAf/vCHTJw4kQkTJjB79mwmT54cNY6m8vLyePfdd5kwYQKvvPIKN910U9S/w+233859993HpEmT+POf/8zvfve7+mPDhg1jxowZnHPOOdx5551kZ2fzrW99iz/96U9MnjyZ9evXN2r1T58+neuuu44TTjiBkSNHHvX3jWXu3LmsW7eufqBcQ7H+Lr/97W+ZMGECkyZNwufzcc4558R1rXiZ62RzvRQVFbkVK1ak/Tor961MKLEf1/s4BuYPTGNEIgLw8ccfc8IJJ7R3GG2u7j55MBjkwgsv5Ktf/Wrcyaczyc/Pp7y8+45Tivbv28xWOueK4jlfLfUokml5lwV0X11E0udnP/sZU6ZMYcKECYwcOZILLrigvUOSDkgD5aKoDCR+j1wj4EUkneq6u7u67txKTwW11KMoDyT+j6o8UE5nu5UhIiJdi5J6FMk8ohZ2YSqD8T3XLiIikg5K6lEk01IHPa8uIiLtS0m9Cedc0slZ99VFRKQ9Kak3URWsSnoimWRb+CLSuWRkZNQvyjFlyhS2bdvGihUruP7664HIzHANFwxpCzU1NVx22WWMHj2amTNnsm3btqPKVFdXM2PGjPpnym+++eb6Yy+//DJTp05lypQpnHLKKWzatCnpWK655pqoC8s0tHLlSiZOnMjo0aO5/vrro45Jcs5x/fXXM3r0aCZNmsT7779ff6zh32DBggVJx9rVKKk30Zou9DJ/mQbLiXQDOTk5rFq1qn4bMWIERUVF3HbbbUBySb21y4nec8899O7dm02bNvG9732PG2+88agyWVlZvPLKK3z44YesWrWKF154oX7WumuvvZa//vWvrFq1ii9+8Yv8/Oc/b/Z6ixYt4tVXX4167O23366fOCeWa6+9lj/+8Y9s3LiRjRs38sILLxxV5vnnn68/vmTJEq699tr6Yw3/Bs8880yz1+pOlNSbaE1rO+RCVAVbt1CEiHROr776Kueddx7btm3jzjvv5NZbb2XKlCn108dGc//997NgwQLmzZvHGWec0arrP/3001x11VUAXHzxxbz88stHNTLMrH6hl0AgQCAQqJ/n3Mw4fPgwEJmudtCgQUnF8fHHHzN27FgyMjJiltmzZw+HDx/mpJNOwsy48sor66elbfqZrrzySsyMk046iYMHD7Jnz56k4uou9Jx6E60d7FYeKCfXl5uiaESkWc8vhr1rUlvngIlwzi+bLVJVVcWUKVMAGDlyZP2qahBZg/uaa64hPz+fH/zgBy1e7v3332f16tX1K7A1dOqpp0adW/w3v/kNZ555ZqN9u3btYujQoQB4vV569uxJaWkpffv2bVQuFAoxbdo0Nm3axLe//W1mzpwJwN133825555LTk4OBQUF9S34RD3//PP1c97HsmvXrkarkw0ZMoRdu45e477hZ2pYbuDAgVRXV1NUVITX62Xx4sWajKeWknoTrb0vXu4v55jcY1IUjYh0RHVdv6lw1llnRU3oQLOt/GRlZGSwatUqDh48yIUXXshHH33EhAkTuPXWW1m6dCkzZ87k17/+NTfccAN33313o3NffPHF+m79HTt28MYbb5Cfn09WVhbvvPNOfZn77rsv5XE3tX37dgYPHsyWLVuYN28eEydO5Nhjj037dTs6JfUGqoJVBMOtu6+l6WJF2lALLerOoLnlRBNpqQ8ePJidO3cyZMgQgsEghw4dorCwMGbdvXr1Yu7cubzwwgv1S5LWtdovu+yyqK3ts88+u3752EWLFrFo0SLmzJlTf7yyspKDBw8yaNAgdu7cyRe+8AUgMnDummuuaRRrcXFx/fvi4mIGDz56kZ26zxStXN3PUaNGMWfOHD744AMldXRPvZFUrLSmx9pEpOmSnH/729/48Y9/nHA9r7/+eqMBeXVb04QOsGDBAv70pz8BkeVh582bd9S64CUlJRw8eBCI3EJ46aWXOP744+nduzeHDh1iw4YNALz00ktJLZqzbNky5s6dC8DQoUPr422Y0AEGDhxY38XvnOOBBx7g/PPPj/qZHnjgAZxzvP322/Ts2ZOBAwdy4MABampqANi/fz9vvvkm48aNSzjerkgt9QZSMXlMIBygJlRDVkZWCiISkc7oC1/4AhdffDFPP/00t99+O5s3b6agoCCt1/za177Gl7/8ZUaPHk2fPn3q1xrfvXs3X//611m6dCl79uzhqquuIhQKEQ6HufTSSznvvPMA+OMf/8jChQvxeDz07t2be++9N+EYnn/+eS6++OK4yv7+979n0aJFVFVVcc4559QvQVq3ROk111zDueeey9KlSxk9ejS5ubn13foff/wx3/zmN/F4PITDYRYvXqykXktLrzawpmQNpdWlra5nQt8J9M3p23JBEUlYZ1x69Utf+hK33nor/fr1a+9Q0mrq1Km88847+Hy+9g6l02rt0qtqqTeQqvvh5f5yJXURqfeXv/ylvUNoEw0nh5H2oXvqtfwhP/6QPyV1aWY5ERFpD0rqtVK5GIsWdhFJr85221AkHqn4d53WpG5m883sEzPbZGaLoxy/1cxW1W4bzOxgOuNpTipb1zWhGgKhQMrqE5EjsrOzKS0tVWKXLsU5R2lpKdnZ2a2qJ2331M0sA7gDOAsoBt4zs2ecc+vqyjjnvteg/L8CJ6Yrnpak+lG0skAZfTKiTyghIskbMmQIxcXFlJSUtHcoIimVnZ3daKa9ZKRzoNwMYJNzbguAmT0MnA+si1H+CuDmGMfSLtX3wcv95fTJVlIXSTWfz8fIkSPbOwyRDimd3e+DgZ0N3hfX7juKmQ0HRgKvxDh+tZmtMLMV6fh2HggHUr4Qi2aWExGRttZRBspdDjzunAtFO+icW+KcK3LOFaXjOc/qYHXK66wKaLU2ERFpW+lM6ruAoQ3eD6ndF83lwENpjKVZNaGalNeZqsfjRERE4pXOpP4eMMbMRppZJpHEfdRK9mZ2PNAbeCuNsTQrHQnYH/YTduGU1ysiIhJL2pK6cy4IXAe8CHwMPOqcW2tmt5jZggZFLwcedu34fEp1KPXd76DWuoiItK20ThPrnFsKLG2y76Ym73+Wzhjika7kWxOqIdvbumcORURE4tVRBsq1q3TcUwe11EVEpG0pqZPelrqIiEhbUVJHLXUREekaun1SD4VDBMPBtNStlrqIiLSlbp/U05l4ldRFRKQtdfukns4ucnW/i4hIW+r2SV0tdRER6Sq6fVJPZ2s65NJ3v15ERKSpbp/U092aVhe8iIi0lW6f1IMuvS3pQDiQ1vpFRETqdPukblh667f01i8iIlKn2yf1dEv3lwYREZE63T6peyy9v4J01y8iIlKn22ecdCddtdRFRKStdPukrnvqIiLSVSippznpqvtdRETaSrfPOOp+FxGRrqLbJ3V1v4uISFehpJ7u7nf9ikVEpI10+4yjlrqIiHQV3T6p6zl1ERHpKrp9xklnS12D5EREpC0pqaexe1xd7yIi0pa87R1Ae4u3e3z5hgBPvBWgtNxRmG8snOVj9lhfSuoWERFJhW6f1H2e5hMzRBL6/cv8+GtXaS0td9y/LLJOenOJPZ66RUREUqXbNyUzMzJbLPPEW4H6hF7HH4zsb23dIiIiqZLWpG5m883sEzPbZGaLY5S51MzWmdlaM3swnfFEk+lpOfGWlruE9idSt4iISKqkrfvdzDKAO4CzgGLgPTN7xjm3rkGZMcCPgZOdcwfM7Jh0xROL1+PFMByxE3RhvkVN4IX5zQ+E82Wo+11ERNpOOlvqM4BNzrktzjk/8DBwfpMy3wDucM4dAHDOfZrGeKIysxaT78JZPjKbfP3J9Eb2N0ctdRERaUvpTOqDgZ0N3hfX7mtoLDDWzN40s7fNbH60iszsajNbYWYrSkpKUh5oS8l39lgfi+ZmcnreZn7kfZjCPFg0N7PF0e+6py4iIm2pvUe/e4ExwBxgCPBPM5vonDvYsJBzbgmwBKCoqKj5G9lJyMzIhObHvDF7rI+LQ7sYu/wZis4/i8rew+KrV0REpI2ks6W+Cxja4P2Q2n0NFQPPOOcCzrmtwAYiSb5NxZt89w+bAUDfHe/GVV6PtImISFtKZ1J/DxhjZiPNLBO4HHimSZmniLTSMbO+RLrjt6QxpqjivfftzyvkcL8x9N3+dnz1qqUuIiJtKG1J3TkXBK4DXgQ+Bh51zq01s1vMbEFtsReBUjNbBywDfuicK01XTLEkMkp9//CTKCjZSGZFy2FqoJyIiLSltD6n7pxb6pwb65w71jn3i9p9Nznnnql97ZxzNzjnxjnnJjrnHk5nPLEkknz3D58JtNwF7zEPGZ6MVsUlIiKSiG4/oxwk1k1e2XMIlQWDWuyCV9e7iIi0NSV1EkzAZuwfPpNeez4iw18Ru051vYuISBtTUifxBLx/+Ew84SCFO1fGrlMtdRERaWNK6kSmik1kmdTD/cbiz+5J3+3vxCyjpC4iIm1NSZ3IVLEJJWFPBvuHzaBP8UosFH3WmqyMrBRFJyIiEh8l9VqJJuH9w0/CG6ii9541KalPRESktZTUayWahA8OmkTIm01hjC54JXUREWlrSuq1Ek3CYW8mnw05MXJfPRxqdX0iIiKt1WJSN7MvmCUwiqyTSmZgW8mI2WRVHaDnvvVHHVNSFxGRthZPsr4M2Ghm/2Nmx6c7oPaSnZGd8Dmlw6YTysjkmK2vN9rv9Xg1m5yIiLS5FpO6c+5LwInAZuB+M3urdn3zHmmPrg0l01IP+XIoHVpEv63LG3XBq5UuIiLtIa5udefcYeBx4GFgIHAh8L6Z/WsaY2tT2d7EW+oAJaNOIbP6EL32rq3fp6QuIiLtIZ576ueb2d+AVwEfMMM5dw4wGfh+esNrO5meTAxL+LzSIdMIebPot/XN+n1K6iIi0h7iaalfBNxau4rar51znwI45yqBr6U1ujZkZgktwVon7MuOdMFvewur7YJXUhcRkfYQT1Lf65z7Z8MdZvYrAOfcy2mJqp0km4w/HXVqpAu+diIaJXUREWkP8ST1s6LsOyfVgXQEyYyAB/hsyFRC3uz6Lvgsr5K6iIi0vZhJ3cyuNbM1wPFmtrrBthVY3XYhtp1kF2EJe7PYP2xGbRd8UC11ERFpF95mjj0IPA/8N7C4wf4y59xnaY2qnSTbUgcoGXky/bf8k16715A59PQURiUiIhKf5rrfnXNuG/BtoKzBhpn1SX9obS/Zx9og0gUf9OXQf+ub+DyJD7gTERFprZZa6ucBKwEHjZ73csCoNMbVLlpzLzzszWT/8JPou205BKrAl5PCyERERFoWs6XunDuv9udI59yo2p91W5dL6AA5Ga1LxPtGz8Hrr4BPnk9RRCIiIvGL2VI3s6nNneicez/14bQvX4YPj3kIu3BS5x8YOJFgXj+8qx+BCRelODoREZHmNdf9/v+aOeaAeSmOpUPI8eZQEahI7mRPBlUnnEeP9/8MFfshr29qgxMREWlGzKTunJvbloF0FNkZ2ckndSA46RJYcR989ATM/GYKIxMREWlec93v85xzr5hZ1H5k59yT6Qur/bRmBDyAb8BkGDARPnxISV1ERNpUc4+01T1s/YUo23nxVG5m883sEzPbZGaLoxxfZGYlZraqdvt6gvGnXGueVYfaLwWTLofdH0DJhhRFJSIi0rLmut9vrv35lWQqNrMM4A4i08wWA++Z2TPOuXVNij7inLsumWukQ2ta6j6PD6/HCxMvgZd+CqsfhjNuSmF0IiIiscWz9Gqhmd1mZu+b2Uoz+52ZFcZR9wxgk3Nui3POT2Qt9vNbG3C6tSap108P26M/HDsPVj8K4eRG0ouIiCQqngVdHgZKgIXAxbWvH4njvMHAzgbvi2v3NbWwdk75x81saBz1plVrut9zvA2ec590ORzaCdvfjH2CiIhICsWT1Ac65/7TObe1dvs50D9F1/8/YIRzbhLwEvCnaIXM7GozW2FmK0pKSlJ06eh8GT4yLCOpcxu18o//PGT2gFUPpigyERGR5sWT1P9uZpebmad2uxR4MY7zdgENW95DavfVc86VOudqat/eDUyLVpFzbolzrsg5V9SvX784Lt06yXbBN1qdLTM3MgHN2r9B9aEURSYiIhJbc0uvlpnZYeAbROaB99duDwNXx1H3e8AYMxtpZpnA5cAzTa4xsMHbBcDHiYWfHsl2wR/1ZWDaVRCsgjWPpSAqERGR5jU393sP51xB7U+Pc85bu3mccwUtVeycCwLXEWnVfww86pxba2a3mNmC2mLXm9laM/sQuB5Y1PqP1HrJttQb3VMHGDQV+k+ElVHvKoiIiKRUc9PE1jOz3sAYoD7bOef+2dJ5zrmlwNIm+25q8PrHwI/jDbatJJvUj2rhm0Va60t/EHlufdCJKYhOREQkungeafs68E8iLe7/qP35s/SG1b6OanHHIdOTSYYnygC7iZeANxvefyAFkYmIiMQWz0C57wDTge2188GfCBxMa1TtLJl76jFb9zm9YNwFsPox8Cc/p7yIiEhL4knq1c65agAzy3LOrQeOS29Y7SuZ7vdmz5l2FfjLIiPhRURE0iSepF5sZr2Ap4CXzOxpYHt6w2pfXo8Xn8eX0DnNJvVhs6DvWA2YExGRtGoxqTvnLnTOHXTO/Qz4KXAPcEG6A2tvibbWczKauQ9vBlOvhOJ3YV/Tqe9FRERSI56WOmY21cyuByYBxbVzuXdpiSb1FstPvgIyMmHFva2ISkREJLZ4Rr/fRGT61kKgL3Cfmf17ugNrb822vKNoMann9YXxF0XWWa8+3IrIREREoounpf4vwHTn3M21y7GeBHw5vWG1v0Ra6oY1niI2lplXg79c88GLiEhaxJPUd9Ng0hkgiyZzuHdFiST1zIxMPBbHr3LwNBgyHd5doiVZRUQk5Zqb+/12M7sNOASsNbP7zew+4CO6+HPqkNiz6glNVjPjm/DZZtj8ShJRiYiIxNbcNLEran+uBBo+YP1q2qLpQBJpqSeU1MedD3/gs1IqAAAgAElEQVT/N3j3LhhzZhKRiYiIRBczqTvn6h+qrl1lbWzt20+cc4F0B9bePOYhKyOLmlBNi2UTGinvzYRpX4HXfgmlm6Hw2FZEKSIickSLC7qY2Rwio9+3AQYMNbOr4lnQpbPL9mbHldTjGiTXUNFXcK//hk9f/Tk7T/52i8X75/ZnaMHQFsuJiEj3Fs8qbf8P+Jxz7hMAMxsLPARMS2dgHUFmRmZc5RJO6j0GUDH2LArXv8CGEy8l5Gu++z5UEVJSFxGRFsUz+t1Xl9ABnHMbgMTmUO2k4h0sl3BSB3accC7eQCUDNvyjxbJVwSoqAloMRkREmhdPUl9hZneb2Zza7Y8cGUTXpaWrpV4ZqOTTPsM42H8cQ9c8jYWDLZ5TUlmS0DVERKT7iSepXwusA66v3dbV7uvy4mmp+zy+6OuoN2N/1X4Adk66iOyKEvptfbPFc0qqlNRFRKR5zd5TN7MM4F7n3L8A/9s2IXUcWd6WW+DxtuYb+rTyUwBKh06jotdQhq3+G5+OOi2y8EsMFYEKKgOV5PpyE76eiIh0D8221J1zIWB47SNt3U483eqJTFIDka738kB55I152DnxQvI/20rvXR+0eK5a6yIi0px4ut+3AG+a2U/N7Ia6Ld2BdQSZnkyM2K1nSLylXtdKr7Pv2NOoye3DsNV/i3FG7HNFREQaiiepbwaerS3bo8HW5ZlZi0k70SVamyZml+GjeMICeu9ZTY+Sjc2eWxGooNxfntD1RESk+2jxOXXn3H8AmFlB5K0rS3tUHUhLE9AkMvK93F9OZbDyqP27jzub4R88ytA1f2PdvB81W8enlZ+Sn5kf9zVFRKT7iGc99SIzWwOsBtaY2Ydm1uUnnqnTUks9kaS+r3Jf1P2hzFx2nzCfftveIudQ8wvg7avch3Mu7muKiEj3EU/3+73At5xzI5xzI4BvA/elNaoOpKWBcPEm9bALs7dib8zjxeMXEPZ4Gb7qsWbrqQnVcKDmQFzXFBGR7iWepB5yzr1e98Y59wbQ8mwpXUSqWuqlVaUEwrHXwfHn9mb3CefQf/Nr5Bza3Wxde8r3xHVNERHpXuJJ6q+Z2V21s8mdbma/B141s6lmNjXdAba35lrqiUw8s7ui+UQNsHPihZHW+ofNt9ZLq0vxh/xxXVdERLqPeBZ0mVz78+Ym+08EHDAv1olmNh/4HZAB3O2c+2WMcguBx4HpzrkONQVtYU4hJw08Keoxa2aymIaqg9UcqG65y9yf25vdx89nyLpn2T7lUqoKBkYtF3Zh9lXs0yIvIiLSSDyj3+cmU3HtbHR3AGcBxcB7ZvaMc25dk3I9gO8A7yRznXTzmCfhx9aaijVALpqdky5i0PoXGL7qUdaf9p2Y5fZW7lVSFxGRRuLpfk/WDGCTc26Lc84PPAycH6XcfwK/AqrTGEu7iqeVXqeutd5/06vkHI5977wiUKEueBERaSSdSX0wsLPB++LaffVq78kPdc49l8Y42lUoHOKw/3BC5+ycdBFhj5dhqx5ttlwiXxZERKTrS2dSb5aZeYgsEvP9OMpebWYrzGxFSUnnmv/8YM1Bwi6c0Dl1rfUBm14l92BxzHJ6tE1ERBqKZ/KZS2rve2Nm/25mT8Y56n0X0PCm75DafXV6ABOIjKTfBpwEPGNmRU0rcs4tcc4VOeeK+vXrF8elO46DNQfjKrd8Q4Dv/6mSRXdU8P0/VfJE/vmEvFmMXPHnmOd8Vv1ZqsIUEZEuIJ6W+k+dc2VmdgpwJnAP8Ic4znsPGGNmI2tXebsceKbuoHPukHOub4NJbd4GFnS00e+tVeZveVbd5RsC3L/MT2l5ZKa40nLH79/M4c1B59Nv+9sU7Fsf9Tx/yE8gFPvZdxER6V7ieaQtVPvz88AS59xzZvbzlk5yzgXN7DrgRSKPtN3rnFtrZrcAK5xzzzRfQ9fg8/haLPPEWwH8Tabz8QfhJ8VnsyznBUa99ydWff6/oq63Hus5+b0Ve9l8cHNcMQ4rGMbQHhpJLyLS2cWT1HeZ2V1EHk37lZllEee9eOfcUmBpk303xSg7J546Oxuvp+VfcV0LvaldFVlsm3k5Y5ffSeHOFZQOm35U3R6L/qco95c3O4NdQ/H0JoiISMcXT3K+lEhr+2zn3EGgD/DDtEbVhfgyWm6pF+ZHn8SmMN/Yc9xZVBYMYtR7D0A41Oh4c70AFcGKuGOsCMRfVkREOq54kvpA4Dnn3EYzmwNcAryb1qi6kHi63xfO8pHZpEGf6Y3sdx4vW4u+RN7BHQzY9GrcdSeSqKuCVVr5TUSkC4gnqT8BhMxsNLCEyIj2B9MaVReS681tsczssT4Wzc2sb7EX5huL5mYye2wkaZeMmM3hfmMYufIvZASq6s/L8eZErS8QDiQ0MU3YhakKVrVcUEREOrR47qmHawe9XQTc7py73cw+SHdgXUWf7D5kZmS2mGRnj/XVJ/GjmLFp5teZ+uyNDPvwcbYWfRmAAXkDohavDFQmHGdFoIJcX8tfQEREpOOKp6UeMLMrgCuBZ2v3tdynLEBk0ZdYyTcRh/sfz95j5zB0zVNkH95DjjeH3tm9o5ZNptVdHeyys/SKiHQb8ST1rwCzgF8457aa2Ugg9owocpSBedFXW0vUlulX4jxeRr9zX7N11oRqEq67OqSkLiLS2bWY1GtXVbsReL/2/Vbn3K/SHVhX0lyrOhH+vEK2T7mEvjveYeDej2OWS6bVncwXARER6VjimSb2C8Aq4IXa91PMrFtMHJNKg/MHt1woDjsnnE9Nz8H4/v5TiDGbXDKtbg2UExHp/OLpfv8ZkWVUDwI451YBo9IYU5dUmF3Y6nXZAVyGj+BZt8D+T+Cdu6KWqQkm3urWMq4iIp1fXAPlnHOHmuxLbNkxwcxS0lovyCwgb/xCGPM5WPZfcHDnUWWSaakHwgFCTSa3ERGRziWepL7WzL4IZJjZGDO7HVie5ri6pIF5A8mw6HO1x2tIjyGROeDP/Q3gYOkPoMHEMYFQIOGlXuvovrqISOcWT1L/V2A8UENk0plDwHfTGVRX5fV4WzUSPisji345tUvP9h4Oc38CG16AdU/Xl2nNKHbdVxcR6dziGf1e6Zz7N+fc9Nrt351zev4pSYPyB7XqXGu4UtvMa2HAJHj+RqiKrNvemufN1VIXEenc4hn9/pKZ9WrwvreZvZjesLquXF8ueb68pM7tm9O38Y4MLyy4DSo+hZf/A2hdS10T0IiIdG7xdL/3rV2dDQDn3AHgmPSF1PUVZhcmfE62Nzv6l4FBJ0Za7CvuhR1vJzXyvY4moBER6dzimvvdzIY553YAmNlwQEt6tUJhTiE7ynYkdk5zXwTm/gTW/x88dS01lyxJOq6O1lLfX7Wfw/7DzZYZmj80ruVtRUS6g3iS+r8Bb5jZa4ABpwJXpzWqLq4gswCfx0cgHH3ymGgKc5pJ6ln5cMEf4P7PU/jG7ZSc9NWk4qoMJr4QTDrtKt/FgeoDzZYpyCw4+raEiEg3Fc9AuReAqcAjwMPANOec7qm3gpnRJ7tP3OU95qFXVq/mC404BTfzWgasfZpeuz9MKq5gONihBsuV+ctSUkZEpLuIZ6DchUQmoHnWOfcsEDSzC9IfWtfWM6tn3GULMgvwWMvDHypO+z6VPQdz/D9vI8NfkVRcFYHkzku1ykAlwXCwxXJK6iIiR8QzUO7mhjPK1Q6auzl9IXUPCSX1rIK4ylUQ5uPTv0tW5WeMfvvupOLqKEk93mStpC4ickQ8ST1amXjuxUszcr25eD3x/Rp7Zsb3BaAiWEFZv7Fsn3wxAze+Qt+tiU/8VxnoGPfV403WgXBAk+aIiNSKJ6mvMLP/NbNja7f/BVamO7CuzswoyIyvBR53S722y337lEs53G8Mx7/x/5Fdti+huMoD5QmVT5eyQPwtcLXWRUQi4mkq/ivwUyID5QBeAr6dtoi6kWE9hrW4zrrXvPg88T2yVZfcXIaPdXN+QNFT3+OEZf+PVef9Fy7OXoGKQAVhF47rHn66OOcSStRl/jKOydXUCSIiLf6f3jlXASxug1i6nV7ZveiV3cKo9jhVB6vxh48sn1pdMIBPTr2O8a/8DyNX/JUtM66Kq56wC1MeKI+7FyEdKoOVCS1Ko5a6iEhEi0ndzJYRZbIZ59y8tEQkSYk2SUvJyJPZffzZDFvzJAcHTeSzIVPjqqvMX9auSf1wTfMTzjRV5i/DOdd4XnwRkW4onj7ZHzR4nQ0sBFp+1kjaVKzW6qaZX6Ng33qOf+23rLzgf6nJa3milsM1h1Oy9nuyWppFrqmQC1EZrEx6Tn0Rka4inslnVjbY3nTO3QDMiadyM5tvZp+Y2SYzO6oL38yuMbM1ZrbKzN4ws3GJfwSB2Ikw7M1i3bwf4QI19Hj4v/jmHQf4/p8qWb4h9mx27d2dnWhSh8Rb9yIiXVE8k8/0abD1NbOzgRafsTKzDOAO4BxgHHBFlKT9oHNuonNuCvA/wP8m/hGkpYFl//i0Pzf4r2UCm7nFez+l5WHuX+aPmdgrg5UJTWGbSqFwKKnH6pL5IiAi0tXE0/2+ksg9dSPS7b4V+Foc580ANjnntgCY2cPA+cC6ugLOuYb/J85DC8UkpW7EeixPvBWgNDCd49wFXO99ijVuJH8JnsUTbwWYPTb6yPoyf1lCU9mmSnmgHJfEPwMldRGR+Ea/j0yy7sHAzgbvi4GZTQuZ2beBG4BMIOrgOzO7mtpFZIYNG5ZkOF1XSwmttDySJG8NXsx4287N3gdYHx7KivLjY57TXkk92W70ykAloXCIDE9GiiMSEek8Yna/m9l0MxvQ4P2VZva0md1mZin7v71z7g7n3LHAjcC/xyizxDlX5Jwr6tevX6ou3WW0dA+8MD8yKtzh4XuBb7HT9eMPmb9lQt7+mOe01z3qw4HkrutwHWbiHBGR9tLcPfW7AD+AmZ0G/BJ4ADgExLNo9y5gaIP3Q2r3xfIwoIViktBSS33hLB+ZtX0yh8njG4Hvk0mQ+7L+B29N9ETYXt3Z5f7kE7O64EWku2suqWc45z6rfX0ZsMQ594Rz7qfA6Djqfg8YY2YjzSwTuBx4pmEBMxvT4O3ngY3xhy4QWS61pUVYZo/1sWhuZn2L/WDeEJ6d9CMKq/cy/uVfYaGjB8W1x5zqrb1ma74QiIh0Bc3dU88wM69zLgicQe097TjOA8A5FzSz64AXgQzgXufcWjO7BVjhnHsGuM7MzgQCwAEgvmnPpF5lML6R4rPH+poMijuRT3p/mxP++TvGvvl7Pjn1emgyeUtloJIcb04Ko21ea5Oyut9FpLtrLjk/BLxmZvuBKuB1ADMbTaQLvkXOuaXA0ib7bmrw+juJBiyNVQerkz5335h55JTtZcQHj1DdYwDbT7ys0fG2bqm39vn4ikAFwXAw7tXvRES6mpj/93PO/cLMXgYGAn93ztU9Z+QhssiLdACtSeoA2068guyyfYx8/0EC2QXsPuGclNWdqERWZoulIlCR0Fr1IiJdSbNNGufc21H2bUhfOJKoVremzfjk1H/FW1PBmOV3EczM49NjTwOgOtS2ST0V98QP+w8rqYtIt9V+62tKSqSii9x5vKyb90MODRjP8a/9lj47VqSs7nilamCeBsuJSHempN7J1YRqUlJP2JvFmrP+jfLCkYx/5Vf03PNRm3a/p+oLRFuPAxAR6UiU1KVeKDOXNWffTHWP/kz6+3/Sc+9HbXbtqkBqknG8TwOIiHRFSuqdXKpHegeyC1h17s+pzu/H+Bduhq3/TGn9saQqGQfDQQJRnrsXEekOlNQ7uQxL/VzngZxerDr359T0GAh/vRS2vJryazSVym5ztdZFpLtSUu/k0vVMdiCnF5sv+C30GQUPXgYb/5GW69RJZVLXfXUR6a6U1Du5dE604nIL4ar/g75j4KHLYM3jabuWWuoiIq2npN7JZXoy01a3z+ODvEJY9BwMPQme+Bq8fWfKr+MP+QmGgymrTy11EemulNQ7uV7ZvdJWd/166tk94UtPwPHnwQs3wis/h/oJBlsvVY/l1dcXTG19IiKdhZJ6J9crq1daBssZRp+cPkd2+LLhkj/B1Cvhn7+Gp74FKUqeqU7q/rA/pfWJiHQWSuqdnMc8R1rUKVSQVRDpfm8owwtfuA3m/AQ+fBAeuAAqSlt9rVQ/guYPKamLSPekpN4FFOYUprzOvtl9ox8wgzk3wsJ7YNdKuPsMKGndcgCpblmHXZhAWM+qi0j3o6TeBRRmpz6pt/hFYeLFkQF0/nK4+0zY8GLS10p19zukvvUvItIZaOHpLsCX4WNKvymEXCgl9ZkZub7clgsOnQ5ffxke+Rd48FI4fTGcfiN4EvuumI4E7A/54/sMIiJdiJJ6F5HOUfDN6j0cvvYSPHsDvPbLSJf8RUsgN/77/OkY2JaO1r+ISEenpC6t58uBC34PQ4rg+RthyemRe+5DZ8R1ejwJePmGAE+8FaC03FGYbyyc5WP2WF/M8hoBLyLdke6pS2qYwfSvwVdfiLy/d37k0bdwy7cEWhqtvnxDgPuX+SktjzwbX1ruuH+Zn+UbYnfbawS8iHRHSuqSWkOK4Jo3YPyFkUlqHjgfDu2KWTwQDhB24WarfOKtAP4mE875g5H9sSipi0h3pKQuqZfdExbeDRf8AXa9D3+YBaseijoLXTyD5Opa6PHuB3W/i0j3pKQu6WEGU74I17wOx4yDp66Bhy6Hw3saFYvnfnphviW0H9RSF5HuSUld0qvw2Mjz7Gf/N2x5DX4/Ez74a32rPZ7ku3CWj8wmQzozvZH9sWj0u4h0R0rqkn6eDJj1rci99n4nwNPfgvvPg5JP4krqs8f6WDQ3s75lXphvLJqb2ezo92A42OK9ehGRrkaPtEnb6TsavvI8fPAAvHQz/OFk8ouuxHPC2YS9Wc2eOnts84+wReMP+cn2ZrcmYhGRTkUtdWlbHg9MWwTXrYCJF9P73XuY8fi36L/pVUhxy1rzv4tId5PWpG5m883sEzPbZGaLoxy/wczWmdlqM3vZzIanMx7pQPL7wYV3svnC2wnk9OSE125l6jM/pOfetSm7hAbLiUh3k7akbmYZwB3AOcA44AozG9ek2AdAkXNuEvA48D/pikc6poMDxrNywW/4+LTvkll1gBOf+wnjX/4lOU1GySdDLXUR6W7SeU99BrDJObcFwMweBs4H1tUVcM4ta1D+beBLaYxHOiB/yA/mYd+YuZSMnM3QNU8zbPUTFO54j71jzmD75Iup6XFMUnUrqYtId5POpD4Y2NngfTEws5nyXwOej3bAzK4GrgYYNmxYquKTDqBh4g17s9h+4qXsOe5Mhn/wKAM3vMSAjS8nndy1/KqIdDcdYvS7mX0JKAJOj3bcObcEWAJQVFQUexox6VRC4VDUx878uX3YePI17JhyMcM+fIKBn/ydARv+wd6xZ7Bj8sVU9+gfV/1qqYtId5POpL4LGNrg/ZDafY2Y2ZnAvwGnO+c0Y0g30lLSrcnry8bZ32TH5IUMXf0kg9a/yMAN/6BkxGx2TlhA2THHtap+EZGuJp1J/T1gjJmNJJLMLwe+2LCAmZ0I3AXMd859msZYpBOryevLpllXs2PSQoas/T8GffJ3jtn6BoeOOY7iCeezf/hJOE9Ge4fZqe2v2p+SenK9ueT6clNSl4gkLm1J3TkXNLPrgBeBDOBe59xaM7sFWOGcewb4NZAPPGZmADuccwvSFZN0LF5PYv/8/HmFbJmxiO1TLmXAxlcYsvb/GP/K/1Cdfwy7TjiHvWPmEcjplXT93VWZv4yP9n+Ukrp6ZfViyjFTUlKXiCTOXJSVszqyoqIit2LFivYOQ1LktZ2v4Ujy32A4RN8d7zHko6fptW8dYY+X/cNmsOe4z3Fg8GSG9BjG6N6jUxtwF/TJZ5+wp6L1jxDWmTFghlrrIilkZiudc0XxlFVTRtqV1+NN/t63J4P9I05i/4iTyD2wIzKgbtOrHLNtOVX5x1A96RKYcS30GtpyXd1UIBxgX+W+lNa5u3y3vkyJtBNNEyvtKlVd5JW9h7H5pK+z/Ir7WDfn+1T3GEDv5XfAbyfAvefAe3dDRWlKrtWV7KvYl/KFb/ZW7iUUDqW0ThGJj1rq0q5Sfd/bZfj49NjT+PTY0xjvyaPf5ldh9WPw3Pfh+Rth1FyYeDGMnQ8N7r93V6lupUNkhbzPqj+jX26/lNctIs1TUpd2lc7BbNZnJAyeDqf+APZ9BGseh4+egL99EzxeGHEqHP95OO5c6Dk4bXF0VGEXpiJQkZa6y/xlSuoi7UBJXdpVdkb6lkatr9sMBkyMbGfcDLtWwPpnYf1zsPQHkW3QiZHkPvoMGHhiZDW5Lq4yUJm2NefLAmVpqVdEmqekLu2qMKcwpSOv62RmZJKfmX/0AY8Hhs6IbGfdAiUbjiT4Zb+IbDl94Ni5cOwZcOw8KBiY8vg6gvJAefrq9qevbhGJTUld2lWvrF54zJPyFmNhdmF8BfuNhX43wKk3QHkJbFkGm16Gza9EuuoB+h0Pw2fD8JNh2Kwu01Vf5k9fazoQDlAdrCbbm76eGBE5mpK6tCuvx0vPrJ4cqD6Q0noLc+JM6g3l94NJl0Y25yL34Te9DNveiAy2W3FvpFyv4ZEEP3x2ZOszKtLF38mkexrdQDhANkrqIm1JSV3aXWF2YUqTusc89M7q3bpKGt6HP+W7EApGkvz25bD9TdjwAnz4YKRsbmHknvygE2HQVBg8FXoMaP0HSbMMS+/UuumuXySdQuEQ7+59t9G+QfmDGF4wvJ0iio+SurS7Ptl9Ulpfr6xeZKR6LvgMLwyaEtlmfSvSki/5BHYsh10rYdcHkS77utsIPQZFkvugKdB/AhwzDnoN61AtektzLOmuXySdakI11IQarzGWrqdFUklJXdpdri+XXG8ulcHKlNSXVNd7oszgmOMjW9FXI/v8lbB3Nex6H3Z/ALvfjwzCq5PZA445AfqPg2PG1/4cB7mp/VITL7XURWKrClYdta86WN0OkSRGSV06hGEFwzhUcygldfXN6ZuSehKWmQvDTopsdaoPQ8l62LcWPl0X+bn2KVh5/5EyOb2hcAwUjoa+o4+87jMKfOm7J62kLhJbtEZGqhoe6aSkLh3CgLwBDMjr+PehE5ZdcOQRujrOQdke2LcO9n8C+zdC6abIyPu6+/QAGBQMjnTb1229hx95XTAYMnxJhxb1kb8olm8I8MRbAUrLHYX5xsJZPmaPbf66ud7c1N8CEWlD0Z4OCYaDVAYqO/SCRUrqIm3NDAoGRbYxZzY+VlMGpZsjSX7/Rji4HQ7uiAzOW/PokXv2AJZRm/SHQo+Bkefpe9RuBYMig/V6DARvVtQwemb1bDHU5RsC3L/Mjz8YeV9a7rh/mR+g2cQeT90iHVmsRz7L/GVK6iISp6weRwbkNRX0w+FdRxL9wR1wYDsc2hmZJe/jPdBkYA8QmUynYBDkHwN5/SC3L+T1xZfXl0H+g5T7svBn9ySQXUDIl9NoMN8TbwXqE3odfzCyX0lduqrqYHXUe+oAB2sO0j+vfxtHFD8ldZHOwpsJfUZGtmicg6oDka79sj1weA+U7YWy3ZHXFZ9GegEq9kPtKN6xTaoIZWQSyC4gkN2TYFYet9Rkc8ibx2HyOOzyOEwuh1wehyrz6PFpIcGsfIKZeQSz8nEN5vFXUpfObHf57pjH9lXuY1SvUfg8yd/6SicldZGuwiwykj63D/Qf33xZfyVU7udA6SZ27nsfX/VhMqsORn5WH8JXdQivv5wTMkrJdxUUUEG2NZms5v8avw15swlm5hD25ZKT2xcy82u3vMiW1aP2df6Rn1n5jfd5s8GXc+RnRla3mIdfOo6wCzc7dXXYhdlbsZehPYa2YVTxU1IX6Y4ycyFzGAUFgzmY5Y05TW/De+pZ+Cmggr7eSi6fVsOJ/Wvw1pTj9Vfgq6nA6y8nw19JAR7AoKYcyveBvxz8FZH3/nLAJRZrRlbkKQBvTss/vVmRLcMHGZm1W0uv4yjr8YInI/LTMo6817P4Xc7eir0tzra4q3wXg/MH47GO94VTSV2kG8vwZFCYU0hJZUnU43X3zSOj3zMJ52dx6qz+DB/r47MYdU7rPy3yTH40zkGgqjbRl9cm+oojiT9YHTne8GewGgLVEKxq8rM6MrCwvKTxsVAwMrYgWEPCXyASZRlHkr3HC+Y58tqTEeWLgDfS81Bfvsn5noxIHeaJfGGofx1ta+l4iuqo/6wGWIMvMnZk31HHo+1rck7DL0QJnxPPdUjwHAi6EKWlH9PThRr+kQE4NGBc/Z7qYDU7y3Z2yNnllNRFurn+uf1jJnWIJPaWHmGrk+PNoUeshA6R/4Fm5kY2jkkw0iSEQxDy126BxF4H/Udeh0MQDoKr/RkOHdlXv7/hsdqfrkGZhue4JmWDNUfK4SJffly4ma2VxyUqLzAxxrFXv/Z0o/fbD2+nf27/DrdokZK6SDfXO6s3GZZBqFHrJDn9cvulIKIU8mSAJydyf16OaPFLQe0Xg3Ddv4naLxp1PR91r507+ni0ffXnEH89Rx1v7hwSqCfKcSJTwG48sLHRPnNHXjcVdmE2HdzEhL4TYpZpD0rqIt1chieDPjl9mm2tx6vdZvOTxJhFuv7RBEEQWbxl7b6VVOYk9uVvf9V+9lbs7VATZ3W8u/wi0ub6Zrc+GWdmZNLD10zXu0gHteXQlqSngN14YGPMZ9rbg5K6iNAnpw+GtVywGYXZhVqZTTqdz6o/Y1f5rqTPD7kQ6z9bj2umq74tqftdRPB5fPTM6snBmoNJ16Gud+lsAuEA6z9bH/N4vOseHKo5RHFZMUML2v/ZdbXURQSItLST5TEPvbN7pzAakfTbeGAj/pA/6rG6ORpKyyMt8Lp1D5ZviP4M+9bDWzvEeutpTepmNt/MPjGzTWa2OMrx08zsfTMLmtnF6YxFRJpXkFWQ9Ln5vvwOOTlsGPIAAAyOSURBVBGHSCwllSV8WvlpzOPNrXsQTdiF+bj043bvhk/bf4VmlgHcAZwDjAOuMLNxTYrtABYBDyIi7SrPl5f0ffV4l3EV6Sh2lO1o9nhdCz3e/QDlgXIO1BxoVVytlc6v1jOATc65Lc45P/AwcH7DAs65bc651YBmQxBpZ16PN+klJQsyk2/li7S1cn95zKVV6xTmR/+CG2t/nT3lseeNbwvpTOqDgZ0N3hfX7kuYmV1tZivMbEVJSeufpRWR6PJ9ybW4kz1PpD3srdjbYpmFs3xkNhlKnumN7G/O/qr9BELNzx2fTp3iJphzbolzrsg5V9SvXwebsUqkC0mmG92wpFv4Iu3BE8fKf7PH+lg0N7O+ZV6Ybyyam9nilMlm1q7jS9L5SNsuoOH4/iG1+0Skg8rxJj6dao43R4PkpFPpldWLHTR/Tx0SW/egTo/MHmR42m+mvnT+l/geMMbMRppZJnA58EwarycirZTrTbzFrVa6dDYFmQWtnmwplp6ZPdNSb7zSltSdc0HgOuBF4GPgUefcWjO7xcwWAJjZdDMrBi4B7jKztemKR0Ralu3NTvh/dsm07kXak9fjZWzvsSmvN8+X1+4T0KR1Rjnn3FJgaZN9NzV4/R6RbnkR6QA85iHbm53QXNbJtO5F2tvA/IH4w362HtqakvqyMrKY1G8SPk9i3fWppmliRaSR/Mz8hJJ6XmZeGqMRSZ/hBcOpCdXENRq+OV7zMrnfZLIyslIUWStiae8ARKRj6eHrQQnxPTpqmB5nk05tbO+xaemKby8asioijfTIjH/51Fxfrka+i3Qg+q9RRBpJ5Fl1zSQn0rGo+11EGvF5fIztPZaQC7VYtmdW+z6+IyKNKamLyFEG5Q9q7xBEJAnqfhcREekilNRFRES6CCV1ERGRLkJJXUREpItQUhcREekilNRFRES6CCV1ERGRLkJJXUREpItQUhcREekilNRFRES6CCV1ERGRLsKcc+0dQ0LMrATY3o4h9AX2t+P125o+b9fVnT4r6PN2ZV39sw53zvWLp2CnS+rtzcxWOOeK2juOtqLP23V1p88K+rxdWXf6rC1R97uIiEgXoaQuIiLSRSipJ25JewfQxvR5u67u9FlBn7cr606ftVm6py4iItJFqKUuIiLSRSipi4iIdBFK6jGY2Xwz+8TMNpnZ4ijHs8zskdrj75jZiLaPMjXMbKiZLTOzdWa21sy+E6XMHDM7ZPb/t3f+wVZVVRz/fHuKvITkVxpKE+E4FUjQEx01xqmYMSQVG23EocRoxjFwiiZzcLCGoZoismksx0xxBskxhwLEHxSEjZQjQiI/HkH8CqeIH6MQP0YIeKz+2Pvq4XLOe/e9e3nnvuP6zJy5++699j5r3b3PXXfvs+9ZWhOP7+Wha62QtEPS+mjL31LKJenB2L/rJDXloWe1SPpYos/WSDooaUqZTJfuW0mPS9orqTmR10fSUklb4mvvjLoToswWSRM6T+uOk2HvLEmb4lhdIKlXRt1Wx329kWHrdEk7E+N1TEbdVr/DC4uZ+VF2AA3ANmAQ0A1YCwwuk5kE/CqmxwFP5613Ffb2B5piuiewOcXezwDP5a1rDW3eAfRrpXwMsBgQcCXwat4618DmBmA34UEWhelb4BqgCWhO5P0EmBrTU4GZKfX6ANvja++Y7p23PR2091rgrJiemWZvLGt13NfbkWHrdOCeNuq1+R1e1MNn6ulcAWw1s+1mdgz4LTC2TGYsMCemfweMkqRO1LFmmNkuM1sd04eAjcBF+WqVO2OBJyywAuglqX/eSlXJKGCbmeX5RMaaY2bLgX1l2cnrcw5wU0rVzwNLzWyfme0HlgKjz5iiNSLNXjNbYmYn4tsVwIBOV+wMkNG3lVDJd3ghcaeezkXAvxLv/83pTu4dmXgxHQD6dop2Z5B4G+FTwKspxVdJWitpsaQhnapY7TFgiaTXJN2ZUl7JGOhqjAOeyigrUt8CXGBmu2J6N3BBikwR+xhgImGVKY22xn1X4e54q+HxjFsrRe3bNnGn7ryDpB7A74EpZnawrHg1Ydl2GPALYGFn61djRppZE3AdMFnSNXkrdCaR1A24EZiXUly0vj0FC+ux74n/7kqaBpwAnswQKcK4fxi4GBgO7AIeyFed+sKdejo7gQ8n3g+Ieakyks4CzgPe6hTtzgCSziY49CfNbH55uZkdNLPDMf0CcLakfp2sZs0ws53xdS+wgLBcl6SSMdCVuA5YbWZ7yguK1reRPaXbJfF1b4pMofpY0h3A9cD4+EPmNCoY93WPme0xsxYzOwk8SroNherb9uBOPZ1VwCWSPhpnOOOARWUyi4DSbtlbgBezLqR6J+4FmA1sNLOfZch8qLRnQNIVhLHTJX/ESDpXUs9SmrDJqLlMbBFwe9wFfyVwILGc2xW5jYyl9yL1bYLk9TkBeCZF5o/AtZJ6xyXca2Nel0PSaOBe4EYzeztDppJxX/eU7W35Iuk2VPIdXkzy3qlXrwdh9/Nmwg7KaTFvBuGiAehOWMrcCqwEBuWtcxW2jiQsT64D1sRjDHAXcFeUuRvYQNhFugK4Om+9q7B3ULRjbbSp1L9JewU8FPt/PTAib72rsPdcgpM+L5FXmL4l/FjZBRwn3Dv9GmF/yzJgC/AnoE+UHQE8lqg7MV7DW4Gv5m1LFfZuJdxDLl2/pX/mXAi8ENOp476ejwxb58Zrch3BUfcvtzW+P+07/L1w+GNiHcdxHKcg+PK74ziO4xQEd+qO4ziOUxDcqTuO4zhOQXCn7jiO4zgFwZ264ziO4xQEd+qO4ziOUxDcqTuO4zhOQXCn7jhVIKklxnRulvRsKY51In+tpNWSrs5b185A0uEat9dL0qTE+4HJ2Npt1G2U9JKkhnacr5uk5fHRz47T5XCn7jjVccTMhpvZpYQQkZPL8ocB9wE/yk3Drk0vYFKbUulMBOabWUulFSyE6VwG3NrBczpOrrhTd5za8Qrp4R0/AOwvz4zP4n4+zuabJd2aKPuypJVxtv9IabYpaZqkzZL+KukpSfeUz15j3vTW2op1Nkp6VNIGSUskNUb522NYy7WS5ralUxYdOO93Jf0jaRvwY+Di2Mas2HRDWv0UxpN45rukeZJ+Gdt/Q9JISXPj5zk7UW9hrOs4XQ536o5TA6KDG8W7QSMaoyPaBDwGfD+l2mjgP2Y2LM70/xDb+gRhpvhpMxsOtADjJV1GCEwxnPBc68sr0Cu1rVh8CfCQmQ0B/gvcrBBL/X7gc3GV4ZsVtFOL814O3AwMI0SUGxFlpwLb4qrHd7Lqp5y/GyEew45E9lBgu5mNBB4hBDG6FxgMfEHSOVGumQo+W8epR/y+keNUR6OkNYQZ+kZgacw/Ep0Zkq4CnpB0qZ0abGE98ICkmcBzZvaXmD8KuAxYFYOnNRJCh/YBFliMwiWpkqhTWW0tB/5pZmui3GvAQKA3MM/M3gQws31ttFOr8/YDnjGzo8BRSc+20nZa/XL6ERw+AJK6E5byfx6zDJhtMfKepBbgWLS5RdIxST3N7FArejhO3eFO3XGq44iZDZf0fkLYzsnAg0kBM3tFIT75B0k4QjPbLKmJMOv+gaRlZjaDECFujpndl2xH0pQMHU5w6qpb92S1jLYGAv9LZLUQHG8Wqe20V74D502jkvpHOPVzGEKIJ38yvh8GPBx1GkBYMUn+4DoHONpOvRwnd3z53XFqQJw9fwP4dvnOaUkfBxooi1Eu6ULgbTP7DTALaIpFy4BbJJ0f5fpI+ghhlntT3NXdE7ghyu8BzpfUNy4hX584TVZbWbwIfElS35J8B9tpr/zLwA2SukvqkbDhENCzlXqpmNl+wr33kmMfSgg5WuKThNCdEBx8KU20/U0zO97e8zpO3vhM3XFqhJm9LmkdcBvvLstDmLVOSNmFPRSYJekkIV7012M7f5d0P7BE0vti2WQzWyHpaYJz2gusivLHJc0AVgI7gU0JnVLbAnZn2LBB0g+Bl+KS9OvAHa2080ZGO+0976p4O2Ed4UfKeuCAmb0l6WWFjYCLCTHuK2UJMJIQT31o/HxKS/GN0fHDqQ4e4LPA8+04j+PUDR5P3XG6KAo73A+b2U/z1qUWSOphZofjrYzlwJ1mtrqK9pqAb5nZV9pZbz4w1cw2d/TcjpMXPlN3HKde+LWkwYR74XOqcegAZrZa0p8lNVT6X/W4a36hO3Snq+IzdcdxHMcpCL5RznEcx3EKgjt1x3EcxykI7tQdx3EcpyC4U3ccx3GcguBO3XEcx3EKgjt1x3EcxykI7tQdx3EcpyD8H2qUmI4UyyiwAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 576x360 with 1 Axes>"
      ]
     },
     "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()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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
}