{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Stablization via Symmetrization\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We've seen how we can implement the \"symmetrized tensor product\" of several qubits in terms of a quantum circuit: depending on the number of qubits we want to symmetrize, we adjoint a number of auxilliary qubits, prepare them in a special initial state, apply a sequence of controlled swaps on the auxilliaries and the original qubits, undo the initialization, and postselect on the auxilliary qubits all being $0/\\uparrow$. The key was that the algorithm requires post-selection, and so is only successful part of the time. By preparing our original qubits in the directions provided by a constellation of \"Majorana stars\", we were thus able to prepare spin-$j$ states in the guise of $2j$ symmetrized qubits. \n",
    "\n",
    "But nothing limits us to symmetrizing qubits. In fact, we can apply the same algorithm to symmetrize quantum systems consisting of several qubits: we simply repeat the controlled swaps on however many qubits comprise each subsystem. In fact, this played a role in the motivation behind the paper which provided our symmetrization algorithm: \"Stabilization Of Quantum Computations by Symmetrization\".\n",
    "\n",
    "The idea is the following: Suppose we have a quantum computer, and we're worried about noise. Instead of just running our quantum circuit, we could run $n$ copies of the circuit in parallel, while periodically using the symmetrization algorithm to project the $n$ copies as a whole into the permutation symmetric subspace. Suppose there were no errors: then naturally, each copy of the circuit would be precisely the same, and symmetrizing across them would have no effect. But if there are some errors that asymmetrically affect the different copies, then it's possible that projecting the copies into the permutation symmetric subspace will project off these errors! Indeed, the permutation symmetric subspace \"is the smallest subspace containing all possible error-free states. It thus corresponds to the 'most probing' test we can legitimately apply, which will be passed by all error-free states,\" since we can't just do any probe we like, since such poking around will generally disturb the quantum computation. Moreover, if we rapidly and periodically do this projection, which is probabilistic, then due to the \"quantum Zeno effect,\" we should be able to stay within the subspace with high probability.\n",
    "\n",
    "The problem, however, is that to do the symmetrization at all, we need to add in extra qubits and perform highly entangling operations, again and again. Given that NISQ-era quantum computers are notoriously unreliable about things like controlled swap operations, one might wonder if the errors introduced by the symmetrization algorithm itself would outweight the errors corrected by symmetrization! And in fact, that sadly seems to be the case--for now!\n",
    "\n",
    "But with an eye to the future, let's check out how it might work. \n",
    "\n",
    "(Incidentally, depending on one's particular quantum circuit, there may be other relevant subspaces that one could project into, without disturbing the computation, and which would provide some much needed stabliziation. The \"stabilization via symmetrization\" technique, however is nice insofar as it can help independently of the structure of one's circuits. Indeed, if such stabilization could be provided by the quantum computer itself that would be ideal.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr>\n",
    "\n",
    "Let's start off with a simple random circuit consisting of two qubits and two layers of gates. Note that we're using the library `pytket` for our circuits. `random_pytket_circuit` returns a dictionary specifying the circuit (for easy manipulation), and requires `build_pytket_circuit` to actually construct the circuit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div id=\"glowscript\" class=\"glowscript\"></div>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "application/javascript": [
       "if (typeof Jupyter !== \"undefined\") { window.__context = { glowscript_container: $(\"#glowscript\").removeAttr(\"id\")};}else{ element.textContent = ' ';}"
      ],
      "text/plain": [
       "<IPython.core.display.Javascript object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from spheres import *\n",
    "\n",
    "n_qubits = 2\n",
    "depth = 2\n",
    "\n",
    "circ_info = random_pytket_circuit(n_qubits=n_qubits, depth=depth)\n",
    "circ = build_pytket_circuit(circ_info)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The resulting circuit will look something like this (the image below was pregenerated, since `pytket` won't display its latex output to jupyter notebooks):\n",
    "\n",
    "<img src=\"img/circ_tex.png\">\n",
    "\n",
    "Each layer consists of some single qubit gates and a CNOT."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "sym_circ_info = symmetrize_pytket_circuit(circ_info,\n",
    "                                          n_copies=2,\n",
    "                                          every=1, \n",
    "                                          pairwise=False)\n",
    "sym_circ = sym_circ_info[\"circuit\"]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The symmetrized circuit might look something like this: \n",
    "\n",
    "<img src=\"img/sym_circ_tex.png\">\n",
    "\n",
    "We have two copies of the circuit running in parallel, and we do symmetrization after each of the two layers."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First, let's evaluate the original circuit analytically to see what the answers should be."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "analytical dist for original circuit:\n",
      "  00: 0.500000\n",
      "  01: 0.000000\n",
      "  10: 0.500000\n",
      "  11: 0.000000\n"
     ]
    }
   ],
   "source": [
    "basis = list(product([0,1], repeat=n_qubits))\n",
    "def display_dist(dist):\n",
    "    global basis\n",
    "    for bitstr in basis:\n",
    "        val = dist[bitstr] if bitstr in dist else 0\n",
    "        print(\"  %s: %f\" % (\"\".join([str(b) for b in bitstr]), val))\n",
    "\n",
    "circ_dist = eval_pytket_circuit_ibm(circ, analytic=True)\n",
    "print(\"analytical dist for original circuit:\")\n",
    "display_dist(circ_dist)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's add in some noise. We can choose between `thermal_noise_model`, `bitflip_noise_model` and `ibmq_16_melbourne_noise_model` model, or add in a new one!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "noisy dist for original circuit:\n",
      "  00: 0.512625\n",
      "  01: 0.002875\n",
      "  10: 0.482500\n",
      "  11: 0.002000\n"
     ]
    }
   ],
   "source": [
    "circ_noise_model = thermal_noise_model(n_qubits=n_qubits, on_qubits=list(range(n_qubits)))\n",
    "noisy_circ_counts = eval_pytket_circuit_ibm(circ.copy().measure_all(), noise_model=circ_noise_model)\n",
    "noisy_circ_dist = probs_from_counts(noisy_circ_counts)\n",
    "\n",
    "print(\"noisy dist for original circuit:\")\n",
    "display_dist(noisy_circ_dist)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We could try to do error mitigation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mitigated dist for original circuit:\n",
      "  00: 0.504593\n",
      "  01: 0.002901\n",
      "  10: 0.490434\n",
      "  11: 0.002073\n"
     ]
    }
   ],
   "source": [
    "circ_meas_filter = qiskit_error_calibration(n_qubits, circ_noise_model)    \n",
    "mitigated_circ_dist = probs_from_counts(qiskit_pytket_counts(circ_meas_filter.apply(pytket_qiskit_counts(noisy_circ_counts))))\n",
    "\n",
    "print(\"mitigated dist for original circuit:\")\n",
    "display_dist(mitigated_circ_dist)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's look at the symmetrized circuit. We evaluate the circuit, do the postselection, and average the results across the copies. Here we haven't included any noise, so this won't be any different from simply running the original circuit with a higher number of shots."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "clean sym circ dists:\n",
      "  postselected dists:\n",
      "    circuit 0:\n",
      "      00: 0.502000\n",
      "      10: 0.498000\n",
      "    circuit 1:\n",
      "      00: 0.514375\n",
      "      10: 0.485625\n",
      "  averaged dists:\n",
      "    00: 0.508188\n",
      "    10: 0.491812\n"
     ]
    }
   ],
   "source": [
    "def display_sym_dists(sym_dists):\n",
    "    global basis\n",
    "    exp_dists, avg_dist = sym_dists[\"exp_dists\"], sym_dists[\"avg_dist\"], \n",
    "    print(\"  postselected dists:\")\n",
    "    for i, dist in enumerate(exp_dists):\n",
    "        print(\"    circuit %d:\" % i)\n",
    "        for bitstr in basis:\n",
    "            print(\"      %s: %f\" % (\"\".join([str(b) for b in bitstr]), dist[bitstr])) if bitstr in dist else None\n",
    "    print(\"  averaged dists:\")\n",
    "    for bitstr in basis:\n",
    "        print(\"    %s: %f\" % (\"\".join([str(b) for b in bitstr]), avg_dist[bitstr])) if bitstr in avg_dist else None\n",
    "        \n",
    "clean_sym_circ_counts = eval_pytket_circuit_ibm(sym_circ)\n",
    "\n",
    "print(\"clean sym circ dists:\")\n",
    "display_sym_dists(process_symmetrized_pytket_counts(sym_circ_info, clean_sym_circ_counts))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let's look at the performance of the symmetrized circuits with some noise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "noisy sym circ dists:\n",
      "  postselected dists:\n",
      "    circuit 0:\n",
      "      00: 0.506192\n",
      "      01: 0.015902\n",
      "      10: 0.460597\n",
      "      11: 0.017309\n",
      "    circuit 1:\n",
      "      00: 0.506192\n",
      "      01: 0.020124\n",
      "      10: 0.458204\n",
      "      11: 0.015480\n",
      "  averaged dists:\n",
      "    00: 0.506192\n",
      "    01: 0.018013\n",
      "    10: 0.459401\n",
      "    11: 0.016395\n"
     ]
    }
   ],
   "source": [
    "sym_circ_noise_model = thermal_noise_model(n_qubits=len(sym_circ.qubits), on_qubits=list(range(len(sym_circ.qubits))))\n",
    "noisy_sym_circ_counts = eval_pytket_circuit_ibm(sym_circ, noise_model=sym_circ_noise_model)\n",
    "\n",
    "print(\"noisy sym circ dists:\")\n",
    "noisy_sym_circ_dists = process_symmetrized_pytket_counts(sym_circ_info, noisy_sym_circ_counts)\n",
    "display_sym_dists(noisy_sym_circ_dists)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We could try doing error mitigation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mitigated sym circ dists:\n",
      "  postselected dists:\n",
      "    circuit 0:\n",
      "      00: 0.494203\n",
      "      01: 0.015483\n",
      "      10: 0.472539\n",
      "      11: 0.017775\n",
      "    circuit 1:\n",
      "      00: 0.498741\n",
      "      01: 0.020111\n",
      "      10: 0.465265\n",
      "      11: 0.015883\n",
      "  averaged dists:\n",
      "    00: 0.496472\n",
      "    01: 0.017797\n",
      "    10: 0.468902\n",
      "    11: 0.016829\n"
     ]
    }
   ],
   "source": [
    "sym_circ_meas_filter = qiskit_error_calibration(len(sym_circ.qubits), sym_circ_noise_model)    \n",
    "mitigated_sym_circ_counts = qiskit_pytket_counts(sym_circ_meas_filter.apply(pytket_qiskit_counts(noisy_sym_circ_counts)))\n",
    "mitigated_sym_circ_dists = process_symmetrized_pytket_counts(sym_circ_info, mitigated_sym_circ_counts)\n",
    "\n",
    "print(\"mitigated sym circ dists:\")\n",
    "display_sym_dists(mitigated_sym_circ_dists)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And now let's look at how well we did."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "expected results: [0.5 0.  0.5 0. ]\n",
      "actual circ results: [0.50459295 0.00290095 0.49043356 0.00207255]\n",
      "actual sym_circ results: [0.496472   0.01779701 0.46890165 0.01682934]\n",
      "\n",
      "circ error: 0.011195\n",
      "sym circ error: 0.039743\n"
     ]
    }
   ],
   "source": [
    "expected_probs = np.array([circ_dist[b] if b in circ_dist else 0 for b in basis])\n",
    "actual_circ_probs = np.array([mitigated_circ_dist[b] if b in mitigated_circ_dist else 0 for b in basis])\n",
    "actual_sym_circ_probs = np.array([mitigated_sym_circ_dists[\"avg_dist\"][b] if b in mitigated_sym_circ_dists[\"avg_dist\"] else 0 for b in basis])\n",
    "\n",
    "print(\"expected results: %s\" % expected_probs)\n",
    "print(\"actual circ results: %s\" % actual_circ_probs)\n",
    "print(\"actual sym_circ results: %s\" % actual_sym_circ_probs)\n",
    "print()\n",
    "print(\"circ error: %f\" % np.linalg.norm(np.array(actual_circ_probs) - np.array(expected_probs)))\n",
    "print(\"sym circ error: %f\" % np.linalg.norm(np.array(actual_sym_circ_probs) - np.array(expected_probs)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To get a better sense of under what conditions the algorithm might be helpful, we can easily run experiments in batches. For example, below we evaluate circuits with two qubits, with two copies to symmetrize over, where we symmetrize every two layers, randomly pairwise instead of across all the circuits, with thermal noise, and where the errors apply to auxilliary qubits as much as the other qubits. We evaluate such circuits for depths from $1$ to $8$, averaging over many randomly generated circuits for each depth, and graph the results (we could have chosen another parameter to vary, besides depth)!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "depth: 1\n",
      "depth: 2\n",
      "depth: 3\n",
      "depth: 4\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA2IAAAEWCAYAAAAEgyRsAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABM3klEQVR4nO3dd5hU5fn/8fdNXRFEmokCCnZUilIFFCxEbJQoERtgiSXWxJ4YYyy/RM3XFo0GlWBvIIpYo4CNIkXQIKCAKEUF6Ujd5f798ZxlZ4fZndlld2Z29/O6rr125rS5z5x6n6eMuTsiIiIiIiKSPtUyHYCIiIiIiEhVo0RMREREREQkzZSIiYiIiIiIpJkSMRERERERkTRTIiYiIiIiIpJmSsRERERERETSLCsTMTNzM9u/iHFnm9m76Y6pIjGzt8xscKbjkJIzs73NbL2ZVS+DZc0ys547H1XlVVGPFTMbYmYfpzjtcDO7o7xjyjZm1t/MFkXH0+Fp/Ny0XaPM7FYzeyYdn1VS5R2bmR1lZnPLa/nlycwWmtnxKUzXIrofqpGOuMqKBf8xs1Vm9mk07FIz+zE6HhtlOkapXGLPNxXtuMnKRKw47v6su/8q/31xSVtJmdntZvaFmeWa2a1lscxMcPcT3f3J8vwMM7s7uslZa2bfmtkfy/Pz0s3MBpvZtGj9FkfrW+4Htbt/5+513T2vDJZ1qLuPL4OwEoouqPl/28xsY8z7s8vrc8tSOo4VyZh/AJdHx9Nn5fEBiS748deoMvysnma2uKyXW1G5+0fuflCm48hWqTysMbPaZjYsus79YGZ/SHHZR8Wd/9dHx8Fp0STdgV5AM3fvZGY1gXuBX0XH44okyzczu8vMVkR/d5mZxYwfamZzo+vOkATrdJ+ZLY0SwX9Fn58/Pj7uPDP7ZyrrvTPMrF10T7Eh+t+uvD9TksuG+5gKl4iVs3nA9cAbmQ6kAngCONjddwO6Ameb2a8zHFMhtnOlSnWAq4HGQGfgOODaMgir1KKLU9Ycs9EFta671wW+A06NGfZs/nQV5amUVGzR8dEkZtA+wKxMxVPR6Ditkm4FDiAcK8cA15tZ72QzRUlw7Pn/FGA98HY0yT7AQnf/OXr/CyCH1I/Hi4B+QFugDXAqcHHM+JnA74DpCea9EegAHAYcCBwB3BwTe2zcvwQ2Ai+nGFepmFkt4DXgGaAB8CTwWjQcM2sSm2imQ/z9UUmP/8pyvsiK+xh3L/YPWEi4Af0cWAO8COSkMN91wPfAUuB8wIH9o3HjgQtjph0CfBzz3oErgQXAT8A9QLX4aYEPo2l/JpwEziDcOI8BVgMrgY/y5031j3Cw3FrCecYDtwOfAOuAd4HGKczXHZgQxbsIGBINrw88BSwHviWcSGK/g0+Ah6JtMgc4Li6W2O/3fGA2sAp4B9gnGm7AfcAyYC3wBXBYSdY7Wk7TaN7rU5y+GuFkOR9YAbwENIzGvUV4ih07/Uzg19Hrg4H/Rtt2LvCbmOmGA48Ab0b7xHXAj0D1mGl+DcwsxTr+AXg9xWl7AouBP0b770Lg7JjxJwOfRd/5oth9DWgR7dM1YrblndH23gj8FvgiZvr/AlNi3n8E9Is5do+PXncCpkaf+SNwb8w8XWL2wZlAz1J8P7Gflb/+NwA/AE8Td4zHHOf554TahBKM76L4HgV2SfGzhwMPEx6grAMmA/vFjO8KTCEcK1OAromOFWB/4INoup+AF2OmK3K/SxLXvwj79PpoG/4SuJ9wLM4BDo+ZvlUUz2rCDUufmHGNgNHR9vuUcK75OJX4ojjuSHGfvYZwPvgeOK+E++x50bhVwCVAR8J1YzXwUNznJTwnleYP2Bf4a7QP3hDtS+spuDbMj9/f4r+XFNZ/F+D/COfiNcDH0bDvouWuj/6OZMfrWbL9L+l1A9iVcPxvi/msvQg30i8RrhXrov2mQ8x8ewEjCdeRb4ArY8bdCowgXO/WAhdG8dxBOB+sB14n7HvPRtNMAVrELOOBaJuvBaYBR8Ut/5kk2y5/3xkcfZc/AX+KGV+bcLwsjf7uB2rHbrOYaW8AlkTfw1yiayLFXG9SiCulfTr6jJuj/WNZtD3qx4w/Nxq3AvgThc+XxV0P8+OoUUSc+fOtA74E+secSzYBedF2XF3E/EsJpVT5728HXoj5PidTcC26lLB/7XD/B/wH+E/0+oK4z36ecBzmHydjUzhnTQAuinl/ATApwed+THTPFDNsKjAg5v1ZwKIi1n8w4T7TYoadAsyItvEEoE2Kx1Nx2/FXhH0z9nO+A3rHfNcLCeexljtxLizJ/dHxFJwzPwc2AzWAPtF2Xk04H7SKWcYO0xcTS3HXs+EUc80uZpkpnW9IctwUs/yFZOA+JtXAPo12wIaEi+clSebpHQVxGOEC8hwlT8TGRZ+3N/AVBTdLiaaNvbj+LVr5mtHfUUQ7P+Gm6F8prHNpE7H5hCcwu0Tv/55knn2infDMKNZGQLto3FOEJyj1op3qK+CCmO8gF/h9NN8ZhIt8w5hY8r+vvoSSvlaEg+xmYEI07gTCzrw7ISlrBewZjTsL+DxJ/DdScMOzgFANIZXv6ipgEtAs2nH/DTwfjRsEfBIz7SGEA7l2tC8tIlwcawCHEy7ch8Qc3GuAboSTYg7h4nRizPJGEW629o6Wu3eKMb+abHvGTNsz2j73RnH3IJz4DooZ3zqKsQ3hWOmX6AQSbcvvgEOjda5HuMg1jrb9j4QTfD3CfrcRaJTgpDIRODd6XRfoEr1uSrhonBTF0yt63yRmG48pxQksF7grWv9dSH4Cu4+QaDSM1uV14G8pft/Do5g7Rd/RsxTcTDQk3ESdG407M3rfKOb7zT9WnifcJOXvO92j4cXud0ni+gloHy1vLOHiPQioTrjhHRdNW5NwnP4RqAUcSzg35O8zLxAu7LsSzqtLKHgglcpxkUoilgvcFsVyErABaFCCffbRaD1/RdhHXwX2IOxjy4Aeyc5J0fgxwI1J4q0TfY/jom3/CNC5qP2riPfbv5cU1v/haF9pGm27roR9O3/da8Qsd0jMtkll/0vpukFc4hENuzX6rk+K4vob0Q1rtK2mAbcQ9ql9CefpE2Lm3UooeagW8/nzgP0IDwO/JFx7jo/if4rohjtaxjmE61YNwnn1B6IbdUqWiD0WfX5bwg1eq2j8bYRrxR5AE8KN8e3x3wdwEOEY2Ctmuft5kutNCnGluk+fH31v+xLOr68AT0fjDiFcJ4+OPv9ewr52fLL4SJ6IDSDcm1Uj3Af8TME1fAhx59y4eRtEy/5FzLDTiR70Rcv8kIJSs1XEPDiKmWdXwrmqZ6JjINF6kPyctYaY45lQwrUuwWcXlYjFJiBnR59dP8H8Yyn8UOnwaLt2JhxPgwnXttokP56K246/B96K++wxwDUx77sQzmMrCOe1c4E6xe2nCbZDSe+PFhKSzuaE4+/AaB/qRTgPXk/Yr2tFyyg0fTGxJLueDaeIa3aSdUzpfMOO+1tW38eksnEXAufEvL8beDTJPMOIuZhEGzc22PEkT8R6x7z/HfB+MdPGXlxvIyQw+ydbt2LiL20idnNczG8nmecmYFSC4dWBLcTc6BGK5cfHfAdLKfx05VMKbrS3f7+Ep/EXxExXjXCDsQ/h4PiKcAIoUalhzPKMcMD/FaiX4jyzKVyCtyfhhiA/0fiZglK7O4Fh0eszgI/ilvVv4C/R6+HAU3HjbwCejV43jNZ9zxKu4/mEJyNJSzij6XsSDuBdY4a9BPy5iOnvB+6LXrdgx0TstrjpPyKU7HUhPEF/ifDw4xhikmcKn1Q+jLZR47hl3UB0wxAz7B1gcAm/o9jP6hntvzkx44dQxAks2od+pnAp1pHANyl+9nDg8Zj3JwFzotfnAp/GTT+RgpLn8RQcK08BQ4l7oJBsv0sS12Mx768AZse8b030pJrwwOgHYo5DQmJ4K+F8sJVQFTh/3P+j4GY/leMilURsI4UTimVECXuK+2zTmPErgDNi3o8Ero5eF3lOSnF7P0544vsm8BuiEpIE08VfG+Lfb/9eilv/KL6NQNsEn5G/7kUlYqnsfyldNyg6EXsv5v0hwMbodWfgu7jpb6Kg5OJW4MO48eMpXCL1f8TcQBKqiM0oZtusyv+eKFki1ixm2KfAwOj1fOCkmHEnEKq8Ffo+COeRZYSEsWbcZxR5vUkhrlT36feB38WMO4iCa9otxNxkEm6Yt1BwvizuerjD/pXk+5wB9I3fD4uYtnm07NjzdK/87zfme1gZxXhTEcs5l/CAKfZ+pNBnx68Hyc9ZeRQ+3x0QzW9x8yRKxO4glDA3IdRAmBzNu2fcdPtEn9MyZtgjRIl+zLC5hIepyY6n4rbjn4lLNAjJx60Jvs/ahPPam9F3/3j8NEVsh9LcHy0Ezo95/2fgpZj31QgP/Xommr6YWIq8nsXEkvCaXZI/ijjfxO9vJVjeQjJwH5Nqe5MfYl5vIDzxKc5ehMw837cpfk6s+Pn3SnG+ewiZ+LtmtsDMbizFZ5dWSb+n5oQLTbz80o7Y7+1bwlO4fEs82tIx4xN9R/sAD5jZajNbTTiwjXCBGUuo3vgwsCxqALtbkpgL8eAzwo3KX1OcbR9gVExMswknxF+4+zpCcfXAaNozCSes/Pk6588XzXs24WSbL3a/gZBUn2pmuxJObh+5+/eprp+Z9SM8ZT7R3X9KdT5glRfUj4eY7WNmnc1snJktN7M1hCovjYtZVvw6fUA4SRwdvR5PuFD0iN4ncgHhgcgcM5tiZqdEw/cBBsR9p90JF5GdsdzdN6U4bRNCCce0mBjejoanqqhjby92PP/EH0v5riccG59a6HHy/Gh4KvtdUX6Meb0xwfvYOBe5+7YEcTYhXMyLOqfuTHyxVrh7bsz77d9jivtsquta5DkpxTgPI1wgZxAePGxOcb5kilr/xoSnx4nO1cmksv+V9LoRL37+nKg9wz7AXnH7xR8J7XXyxZ9bIPXtiJlda2azzWxNtPz6FH8uS3Udijp+E17n3H0eoU3vrYRr2Qtmlj9dkdebFOIqyfEbH2eN6DMK3Q9F14XYzipKHZ+ZDTKzGTHzHkbq3//66H/sNX83QslFfqwLCaUzLQj3CYkMJtzgexHjE0l2zlqfIK71KX7GnYRq1DMIJaivEhKiH+OmO5dwU/1NXFzXxMXVnLANkx1PxW3H+PXJX6d1ccOIzmefR/FvIWzTVJTm/ih+WKH9OLoeLaLw+SrRMuIVdz3LV+LzXhmeb1KVlvuY8mr4/z1h5823d9z4nwkB50t0wxA//9JUPtjd17n7Ne6+L6Gu6x/M7LhU5s2ARYQqIPF+Ipw49okZtjfhyUS+pnGNO4v6jhYBF7v77jF/u7j7BAB3f9Dd2xOepB5IaFdVGjWKWJdEFhESm9iYctw9f/2eB840syMJN0DjYub7IG6+uu5+acyyC52oo2VOJJQgnUuo55uSqNHyY4TGm1+kOl+kQZT85YvdPs8Riq+bu3t9QvWX4hrqxl984hOxD0iSiLn71+5+JqFazV3AiCi+RYQSsdjvdFd3/3vqq5pSzIWOeTOLPeZ/ItzUHBoTQ30PjWd31lIKH0ew47EUAnb/wd1/6+57EUqg/2WhR9ZU9ruyiLN5XGcs+XEuJ5SwFnVOTUd8Jd1ni1PsOSkZd+9CKP2tAYw1s0/N7HJL3iX2BpJfdxL5iVAtLdH5LdmNYcr7XwpKcqML4Xv+Ju57rufuJ+3EMrczs6MIDy9+Q6jCuTuh6lNZdjoQ//0VeS/g7s+5e/doeiec5yD59aa84swl3PwXuh8yszqE6lX5ShWfme1DuD5dTqjqujvwPwq+/2K3rbuvimJrGzO4LTEdapjZyYSn+u8THnLHx9CccC16qrjPSiDZOWtWcXEVx903uvvl7t40ug9cAUyLSwogVG2O7zF3EXBnXFx13P15kh9PxW3HWUCbuHu2NhT+rhtF57FPCVUmqwPHROe7VJT4/ijBsEL7cRRvcwqfr1I5ZxR3PSuVNJ1v4qXlPqa8ErGXgCFmdkh00vlL3PgZwK/NrE50o3NBgmVcZ2YNogP9KkInIYn8SKirC4CZnWJm+0c70BrCE4n4AzAhM6tpZjmE76WGmeXk9yxjBd0Ut0hlWSl6FjjezH5jZjWiA7Gdh67LXwLuNLN60Qn3D4TSnXx7AFdGMQ8gtLd4M8FnPArcZGaHRutRP5oeM+sYPemuSdjBNpHCd2Vm1czs4mj7mJl1Ai4jnKzzp1locd3KxsV0Z7Re+T0G9Y0Z/ybhZHAbocOE/JjGAAea2bnReteM1qFVkpCfIhzArQl195Mys2MJ2+c0d/80wfjhZjY8yWL+ama1ohPIKRT0zFQPWOnum6Lv7qxUYooxgVD1pROh2tMsoqdhhCqIidbnHDNrEn2Xq6PB2ygoMTzBzKpH+3xPM2tWwpiSmQkcaqEL3xzCk2tg+1O3x4D7zGyPKN6mZnZCTPxuZj1L8blvEvaZs6Jj7AzCQ4cx8ROa2YCY9V5FOAlvo/T7XUlMJiQK10fL70moBvZCdD54Bbg1OmceQngKnS8d8e3sPhuryHNSqtx9trtfT7hJ+CvhIcRCM0t0Lck3Azgr2s97R/Ok8lnbCNXt7zWzvaL5jzSz2oQkeRsx16A4Ke9/KfgRaGRm9VOc/lNgnZndYGa7RHEfZmYdS/HZidQjJBvLCdfLW9jxqf/Oeh64ObpGNCZU89vht8nM7CAzOzbaJpso6NgEkl9vyirO35tZSzOrS6g6/KKHEtYRwClm1t1CL3m3Ufjeq7Tx7Uo4Ry2P5juPwqUnPwLNos8sylOE77eBmR1M6AxqeLS8xoRqwBcSzjenmtlJcfOfS2jfWdLS4mTnrKcID9GbWijZvCY/rii2WtF1xICa0XWrWjSuaXScmpl1IVS3K3QPamZdCaUz8b0lPgZcEt0XmZntamYnm1k9kh9PxW3H8YR70SstdK9/eTR8bDTtBYRqcT0I57Pm7n6Du8+Oi3uImS0s5XeaipeAk83sOAv3hdcQ2mym9JAsRpHXsxIuJ1Y6zjfJ7NR9TFHKJRFz97cIbQjGEqoJjo2b5D5CkeuPhCcSz7Kj1wgNI2cQqqo9UcTH3Qo8aaEo8DeEusTvEYqCJxI65xgHYGaPmtmjxYT+GOEEfiah0f5GwokGwgX/W3Yio4/n7t8R6sZeQ6ieM4OCp0BXEJKjBYR60M8RbgbyTSas60+EovjTPcFvc7j7KMKTwRfMbC3hidmJ0ejdCOu8ioIene6B7T9KWtwTqP4U9Nb0DPDP6C+/q9ZGhIariTxAeLr+rpmti6brHBPzZsKN5/HReucPX0doND2Q8MTlBwoaUhZnFFG1AXffEMWY/8PJ8aW1+f5MKPZ+0wp+U+KtmPHNCfXQi/ID4XtdSti/L3H3OdG43wG3Ret+C+HklzIPVVumA7PcfUs0eCLwrbsvK2K23sAsM1tP+P4HRk8OFxE6T/gj4QS3iFAqmn9R+2PcepeKu39FuAF5D/iasE/HuoFwrpgU7afvEZLN/Keu6wg9c5b0c1cQkuBrCPv39cApnriaaUdgcvQdjQaucvcFO7HflSTOLYQL1YmEY/pfwKCYfeZyQtWNHwg3JP+Jmbfc42Mn99lYSc5JWPiR7ZR+l9Dd89z9DXcfQDjGJxYz+VWE73g1ocrOqyUI+1rC/jeFcK6+i9D+YQNRr6bRNajQ0+sS7n/FivaF54EF0WcVW10/SuBPAdoR2vD8RLixTjWRS+YdQtWbrwjXj02kVm2pJO4gdL7wOeH7nx4Ni1cb+DthHX8gPKi8KRpX7PWmjAwj1Lb4kPBdbyJcw4kelF1GuJZ9T7guxP4eXKnic/cvCW34JhLupVpT+Jo0llDi8oOZFbW//YVwHf+WUJviHnd/Oxo3FHjN3d+M9uMLgMetcMlzolKlpFI4Z/2b0NHBF4TzwxvRsHzvEu7PukZxbiTUEIFQcj2BcP/0JKHjn3fjQhgMvBLFERvXVEIy+hBhO80jtAtK5XgqcjtG5/d+hO9rNaHdeb+46/c+7j4gOp/lFfHVFXnfURbXAXefS+gQ45/R+p1KqBG0pdgZd1xOsutZaZT6fJMN9zHFxuclqtZbembmwAEe6nJXOGZ2M6G+6L+TTlz+sQwhdDDQPdOxJGJm3YHLPFSFywpmNp9QHeq9MlhWLcKTkTbuvjXB+J6ERqNlXapUJZnZOYTi/puSTiwiIlIJmdm7hAeEs5NOLBVGpfhBtnRw90RP4SQBd/+YHZ8UZIyZnUaowhFfMlsq0dOesqz2JcVw9x2qIomIiFQl7v6rTMcgZa/UVROjor71Cf52uvivMomq+CX6nlJqeCo7x8zGE7qkvSxBY12RnWKhd8VEx/fZmY4tns7ZkgnZeg3M1rhEKoKY5h2J/opq8pFsmUcVtcyyjj+bpK1qooiIiIiIiATl1WuiiIiIiIiIFEFtxEQkqcaNG3uLFi0yHYaISIUxbdq0n9y9JD9MLyJVjBIxEUmqRYsWTJ06NdNhiIhUGGb2baZjEJHspqqJIiIiIiIiaaZETEREREREJM2UiImIiIiIiKSZ2oiJSKls3bqVxYsXs2nTpkyHImmSk5NDs2bNqFmzZqZDERERqfCUiIlIqSxevJh69erRokULzCzT4Ug5c3dWrFjB4sWLadmyZabDERERqfBUNVFESmXTpk00atRISVgVYWY0atRIJaAiIiJlRImYiJSakrCqRdtbRESk7KhqooiISAZt2wb/+Q8sWwb16kHdukX/r1sXcnJAObGISMWnRExEJAVdu3ZlwoQJmQ5DKpmNG+Hcc2HkyNTnqV69+GStuCSuqHHVq5ffOoqISGJKxEREUpAoCcvNzaVGjbI5jcYvK9Vll2UMkl7Ll0OfPjB5Mtx7L1x6KaxfD+vWpfY/fth33xV+v2FD6rHsskvpk7hE/1VqJyKSnK7eIlJh9evXj0WLFrFp0yauuuoqLrroIh599FHmz5/PPffcA8Dw4cOZOnUqDz30ELfffjvPPPMMTZo0oXnz5rRv355rr7220DJ//PFHLrnkEhYsWADAI488QteuXalbty7r169n/Pjx/PnPf6ZBgwbMmTOH2bNnc8MNN/D2229TrVo1fvvb33LFFVcUWub8+fO57LLLWL58OXXq1OGxxx7j4IMPZsiQIeTk5PDZZ5/RrVs3Vq5cWej9oEGDuOSSS9iwYQP77bcfw4YNo0GDBvTs2ZN27drx8ccfc+aZZ3LNNdek5wuXMvPVV3DSSbBkCbz8Mpx2WhiekwONG5fNZ+Tlwc8/p5bEJfq/ahUsWlR4eG5uap9drVrpSu2KS/z0vEFEKhud1kRk5119NcyYUbbLbNcO7r+/2EmGDRtGw4YN2bhxIx07duS0007jtNNO48gjj9yeiL344ov86U9/YsqUKYwcOZKZM2eydetWjjjiCNq3b7/DMq+88kp69OjBqFGjyMvLY/369TtMM336dP73v//RsmVLHnnkERYuXMiMGTOoUaMGK1eu3GH6/ATxgAMOYPLkyfzud79j7NixQPgZgAkTJlC9enWGDBlS6H2bNm345z//SY8ePbjlllv461//yv3Rd7JlyxamTp1asu9UssLHH0PfviFZGTcOunQpn8+pXh122y38lZXNm3eu1G7x4sLvf/459c/Oydn5Kpix/3fZRaV2IpJZSsREpMJ68MEHGTVqFACLFi3i66+/pkuXLuy7775MmjSJAw44gDlz5tCtWzceeOAB+vbtS05ODjk5OZx66qkJlzl27FieeuopAKpXr079+vV3mKZTp07bf0vrvffe45JLLtlePbBhw4aFpl2/fj0TJkxgwIAB24dt3rx5++sBAwZQPaaBTv77NWvWsHr1anr06AHA4MGDCy3jjDPOSP2Lkqzx4osweDDsvTe89Rbst1+mIyqZ2rXDX6NGZbO8bdt2LLUrSaK3Zk0oVYydd8uW1D47v9SuuOStQQO4666yWVcRkXhKxERk5yUpuSoP48eP57333mPixInUqVOHnj17bv+Nq4EDB/LSSy9x8MEH079//zLvdn3XXXdNedpt27ax++67M6OIEsP4ZaW67JLEIJnnDnffDTfeCN27w6uvll0yU5FVqxYSnnr1YM89y2aZW7bsXKnd0qUF72vUUCImIuVHiZiIVEhr1qyhQYMG1KlThzlz5jBp0qTt4/r378+dd97JZ599xl3RXVS3bt24+OKLuemmm8jNzWXMmDFcdNFFOyz3uOOO45FHHuHqq6/eXjUxUalYvl69evHvf/+bY445ZnvVxNhSsd12242WLVvy8ssvM2DAANydzz//nLZt2xa7fvXr16dBgwZ89NFHHHXUUTz99NPbS8ekYsnNhcsvh3//G844A4YPD9XspHzUqgUNG4Y/EZFsph90FpEKqXfv3uTm5tKqVStuvPFGusQ0tGnQoAGtWrXi22+/pVOnTgB07NiRPn360KZNG0488URat26dMMF64IEHGDduHK1bt6Z9+/Z8+eWXxcZx4YUXsvfee9OmTRvatm3Lc889t8M0zz77LE888QRt27bl0EMP5bXXXktpHZ988kmuu+462rRpw4wZM7jllltSmk+yx7p1oWfEf/87lIY995ySMBERCczdMx2DiGS5Dh06eHzHELNnz6ZVq1YZiqh01q9fT926ddmwYQNHH300Q4cO5Ygjjsh0WBVKRdzumbJ0KZx8MnzxBfzrX5CgAFYqMTOb5u4dMh2HiGQvVU0UkSrjoosu4ssvv2TTpk0MHjxYSZiUmy++CN3Tr14Nr78OJ56Y6YhERCTbKBETkSojUbVBkbL23/+G3wWrVw8++ij8EoOIiEg8tRETyWJm1tvM5prZPDO7McH42mb2YjR+spm1iIY3MrNxZrbezB6Kmb6Omb1hZnPMbJaZ/T2NqyNS6Q0bFkrCWrSASZOUhImISNGUiIlkKTOrDjwMnAgcApxpZofETXYBsMrd9wfuA/I7Wt4E/Bm4NsGi/+HuBwOHA93MTJWmRHaSO/z5z3DBBXDMMaEkrHnzTEclIiLZTImYSPbqBMxz9wXuvgV4AegbN01f4Mno9QjgODMzd//Z3T8mJGTbufsGdx8Xvd4CTAealedKiFR2mzfDoEFwxx1w/vnwxhtQzC8eiIiIAErERLJZU2BRzPvF0bCE07h7LrAGSOlnYs1sd+BU4P0ixl9kZlPNbOry5ctLFrlIFbFqFfTuDc88A7ffDo8/DjVrZjoqERGpCJSIiVRBZlYDeB540N0XJJrG3Ye6ewd379CkSZP0BljFLF26lNNPPz3TYUgJLVwI3brBJ5/A00/DzTeDWaajEhGRikKJmEj2WgLEtjJpFg1LOE2UXNUHVqSw7KHA1+5+/86HKTtrr732YsSIETsMz83NLbPPiF9WqsvOy8srsxgqkylToHNn+P57ePddOOecTEckIiIVjRIxkew1BTjAzFqaWS1gIDA6bprRwODo9enAWE/yK+1mdgchYbu6bMNNv379+tG+fXsOPfRQhg4dCsCjjz7Kddddt32a4cOHc/nllwNw++23c9BBB9G9e3fOPPNM/vGPf+ywzJdffpnDDjuMtm3bcvTRRwNw9NFHM2PGjO3TdO/enZkzZ3LrrbcyePBgjjrqKPbZZx9eeeUVrr/+elq3bk3v3r3ZunXrDsufN28exx9/PG3btuWII45g/vz5LFy4kMMOO2x7vH369OHYY4/luOOOY/369Zx33nm0bt2aNm3aMHLkyB2WOW3aNHr06EH79u054YQT+P777wHo2bMnV199NR06dOCBBx7Y4f3777/P4YcfTuvWrTn//PPZvHkzAC1atOCGG27giCOO4OWXXy7NpqnURo+Gnj2hTh2YMCG8FhERKSn9jphIlnL3XDO7HHgHqA4Mc/dZZnYbMNXdRwNPAE+b2TxgJSFZA8DMFgK7AbXMrB/wK2At8CdgDjDdQj2qh9z98Z2J9eqrISZPKRPt2sH99xc/zbBhw2jYsCEbN26kY8eOnHbaaZx22mkceeSR3HPPPQC8+OKL/OlPf2LKlCmMHDmSmTNnsnXrVo444gjat2+/wzJvu+023nnnHZo2bcrq1asBuOCCCxg+fDj3338/X331FZs2baJt27aMGjWK+fPnM27cOL788kuOPPJIRo4cyd13303//v1544036NevX6Hln3322dx4443079+fTZs2sW3bNpYtW1ZomunTp/P555/TsGFDbrjhBurXr88XX3wBwKpVqwpNu3XrVq644gpee+01mjRpsn19hw0bBsCWLVuYOnUqAK+//vr295s2beKAAw7g/fff58ADD2TQoEE88sgjXH311QA0atSI6dOnJ9lKVc8//wlXXQXt28OYMfCLX2Q6IhERqahUIiaSxdz9TXc/0N33c/c7o2G3REkY7r7J3Qe4+/7u3im2vZe7t3D3hu5e192bufuX7r7Y3c3dW7l7u+hvp5KwTHrwwQdp27YtXbp0YdGiRXz99dc0adKEfffdl0mTJrFixQrmzJlDt27d+OSTT+jbty85OTnUq1ePU089NeEyu3XrxpAhQ3jssce2V8sbMGAAY8aMYevWrQwbNowhQ4Zsn/7EE0+kZs2atG7dmry8PHr37g1A69atWbhwYaFlr1u3jiVLltC/f38AcnJyqFOnzg4x9OrVi4YNGwLw3nvvcdlll20f16BBg0LTzp07l//973/06tWLdu3acccdd7B48eLt488444xC0+e/nzt3Li1btuTAAw8EYPDgwXz44YdFzlfV5eXBH/4AV14JffrA+PFKwkREZOeoRExEdlqykqvyMH78eN577z0mTpxInTp16NmzJ5s2hd76Bw4cyEsvvcTBBx9M//79sRL0oPDoo48yefJk3njjDdq3b8+0adNo1KgRvXr14rXXXuOll15i2rRp26evXbs2ANWqVaNmzZrbP6tatWqlbuO16667pjytu3PooYcyceLElJaV6rJLEkNlt2FDaAM2alRIxO69F6pXz3RUIiJS0alETEQqpDVr1tCgQQPq1KnDnDlzmDRp0vZx/fv357XXXuP5559n4MBQW7Nbt268/vrrbNq0ifXr1zNmzJiEy50/fz6dO3fmtttuo0mTJixaFH5B4MILL+TKK6+kY8eOO5RKpapevXo0a9aMV199FYDNmzezYcOGYufp1asXDz/88Pb38VUTDzroIJYvX749Edu6dSuzZs1KGstBBx3EwoULmTdvHgBPP/00PXr0KMnqVAnLlsGxx8Krr8J998EDDygJExGRsqFETEQqpN69e5Obm0urVq248cYb6dKly/ZxDRo0oFWrVnz77bd06tQJgI4dO9KnTx/atGnDiSeeSOvWramf4Fd3r7vuOlq3bs1hhx1G165dadu2LQDt27dnt91247zzztupuJ9++mkefPBB2rRpQ9euXfnhhx+Knf7mm29m1apV2zsQGTduXKHxtWrVYsSIEdxwww20bduWdu3aMWHChKRx5OTk8J///IcBAwbQunVrqlWrxiWXXLJT61bZzJ0LRx4JM2fCiBGhLaSIiEhZsSQdrImI0KFDB8/v8CHf7NmzadWqVYYiKp3169dTt25dNmzYwNFHH83QoUM54ogjUpp36dKl9OzZkzlz5lCtWtV9hlURt3tpfPQR9O0LNWrA66+HrupFSsLMprl7h0zHISLZq+reTYhIlXPRRRfRrl07jjjiCE477bSUk7CnnnqKzp07c+edd1bpJKyqeOEFOP54aNIEJk5UEiYiIuVDnXWISJXx3HPPlWq+QYMGMWjQoDKORrKNO9x1F9x0Exx1VGgXFnVeKSIiUub0aFdESk1Vm6uWyry9c3Ph4otDEjZwILz7rpIwEREpX0rERKRUcnJyWLFiRaW+OZcC7s6KFSvIycnJdChlbt06OPVUeOyxkIg9+yxUwtUUEZEso6qJIlIqzZo1Y/HixSxfvjzToUia5OTk0KxZs0yHUaaWLIGTT4b//Q+GDoXf/jbTEYmISFWhRExESqVmzZq0bNky02GIlNrnn4ckbPVqGDMGevfOdEQiIlKVqGqiiIhUOe++C927hw46Pv5YSZiIiKSfEjEREalShg2Dk06Cli1h0iSIfrNbREQkrZSIiYhIleAON98MF1wAxx0XfrS5kjV5ExGRCkRtxEREpNLbvBnOPx+eey4kYo88AjVrZjoqERGpypSIiYhIpbZqFfTvDx98AHfeGbqoN8t0VCIiUtUpERMRkUrrm29Ce7AFC8Lvg511VqYjEhERCZSIiYhIpfTpp+GHmrduDb0k9uiR6YhEREQKqLMOERGpdF57DXr2hDp1YMIEJWEiIpJ9lIiJiEil8uCDoU1Y69ahe/qDD850RCIiIjtSIiYiIpVCXh5cfTVcdRX06QPjxsEvfpHpqERERBJTIiYiIhXehg1w+unwwAMhERs5MlRLFBERyVbqrENERCq0ZctCpxxTpsD994dETEREJNspERMRkQprzpzQPf0PP8Arr0C/fpmOSEREJDVKxEREpEL68MOQeNWsCePHQ6dOmY5IREQkdWojJiIiFc5zz0GvXrDHHjBxopIwERGpeJSIiYhIheEO/+//wdlnQ5cu4TfC9t0301GJiIiUnKomiohIhbB1K/zud/D443DWWTBsGNSunemoRERESkclYiIikvXWrg09Iz7+OPzpT/DMM0rCRESkYlOJmIiIZLXFi+Hkk2HWLHjsMbjwwkxHJCIisvOUiImISNaaOTMkYWvXwptvwq9+lemIREREyoaqJoqISFZ65x3o3j28/ugjJWEiIlK5KBETEZGs8/jjoSRsv/1g8mRo2zbTEYmIiJQtJWIiIpI1tm0LnXH89rdw/PHhR5ubNs10VCIiImVPbcRERCQrbN4M550Hzz8fErGHH4aaNTMdlYiISPlQIiYiIhm3ciX06xfagv3tb3DDDWCW6ahERETKj6omimQxM+ttZnPNbJ6Z3ZhgfG0zezEaP9nMWkTDG5nZODNbb2YPxc3T3sy+iOZ50Ey3u5JZCxZA166hLdhzz8GNNyoJExGRyk+JmEiWMrPqwMPAicAhwJlmdkjcZBcAq9x9f+A+4K5o+Cbgz8C1CRb9CPBb4IDor3fZRy+Smk8/hS5dYNkyeO89OPPMTEckIiKSHkrERLJXJ2Ceuy9w9y3AC0DfuGn6Ak9Gr0cAx5mZufvP7v4xISHbzsz2BHZz90nu7sBTQL/yXAmRoowaBT17Qt26MHEiHHVUpiMSERFJHyViItmrKbAo5v3iaFjCadw9F1gDNEqyzMVJlgmAmV1kZlPNbOry5ctLGLpI8e6/H047Ddq0gUmT4KCDMh2RiIhIeikRE5GE3H2ou3dw9w5NmjTJdDhSSeTlwVVXwe9/HzrnGDsW9tgj01GJiIiknxIxkey1BGge875ZNCzhNGZWA6gPrEiyzGZJlilSLjZsCKVgDz4YErGXX4Y6dTIdlYiISGYoERPJXlOAA8yspZnVAgYCo+OmGQ0Mjl6fDoyN2n4l5O7fA2vNrEvUW+Ig4LWyD12ksB9/DO3BRo8Oidi990L16pmOSkREJHP0O2IiWcrdc83scuAdoDowzN1nmdltwFR3Hw08ATxtZvOAlYRkDQAzWwjsBtQys37Ar9z9S+B3wHBgF+Ct6E+k3MyeDSedFJKxUaOgb3yXMyIiIlWQFfPwXEQEgA4dOvjUqVMzHYZUQB98ENqC1aoFY8ZAx46ZjkgkPcxsmrt3yHQcIpK9VDVRRETKxbPPQq9esOeeoWdEJWEiIiIFlIiJiEiZcoc774RzzoGuXeGTT6Bly0xHJSIikl3URkxERMrM1q1w6aXwxBMhEXv8cahdO9NRiYiIZB+ViImISJlYuxZOPjkkYTffDE89pSRMRESkKCoRExGRnbZ4cegZcfbskIidf36mIxIREcluSsRERGSnzJgRSsLWrYM33wwddIiIiEjxVDVRRERK7e234aijoFq10CmHkjAREZHUKBETEZFSGToUTjkF9t8fJk+G1q0zHZGIiEjFoURMRERKZNs2uOkmuPhi+NWv4MMPYa+9Mh2ViIhIxaI2YiIikrJNm+C88+CFF0Ii9tBDUENXEhERkRLT5VNERFKyYgX06wcffwx//ztcfz2YZToqERGRikmJmIiIJDV/fuiefuHCUBp2xhmZjkhERKRiUyImIiLFmjQJ+vSBvDx4/33o3j3TEYmIiFR86qxDRESK9MorcMwxUK8eTJigJExERKSsKBETEZEduMN998Hpp0O7dqFU7KCDMh2ViIhI5aFETETKzYoVsHVrpqOQksrLg6uugj/8Afr3h7FjoUmTTEclIiJSuaiNmIiUm0MOgWXLoH59aNw4/DVqlPx/o0ZQq1amo6+afv4ZzjoLRo8Oidg990A1PbITEREpc0rERKTc3HYb/PhjKBn76afw/8cf4csvw/v164uet169kidvOTnpW7fK6Icf4NRTYfp0+Oc/4fLLMx2RiIhI5aVETETKzcUXFz9+8+aQnOUnavnJWqL/X30V/q9dW/Tydt21ZMlb48awyy5lu84V1ezZcOKJsHw5vPpqSMhERESk/CgRE5GMqV0b9tor/KVqyxZYubL4pC3///z54fXq1UUvb5ddSp681alTuX7IePz40Basdm344APo0CHTEYmIiFR+SsREpEKpVQt++cvwl6rc3KKTt/hh334b/q9cWfTycnJST9ry/9etm53J2zPPwPnnw/77w5tvQosWmY5IRESkalAiJiKVXo0asMce4S9VeXmwalXyKpMrVsDMmQVVLN0TL69WrdSTtvzXu+1WfsmbO9xxB9xyC/TsGX4vrEGD8vksERER2ZESMRGRBKpXL0iKUpWXF6pBJqsyuWIFzJoVXq9cGeZLpEaNkidv9esn7+Vw61a45BIYNgzOPRcef1y9VIqIiKSbEjERkTJSvXpBD44HHpjaPNu2wZo1RSdtsa/nzi14n5tbdAwNGxaftD3/PPz3v6E07NZbs7PKpIiISGWnRExEJIOqVQtVAhs0CO20UuEeeo9MpbfJ+fNh8uTwOv/HtWvUCKVh551XfuslIiIixVMiJlLOzMyAZu6+KNOxSOVgFqog1q8P++6b2jzu4Xfbfvop9BRZks5OREREpOwpERMpZ+7uZvYm0DrTsUjVZRZ+JLtevUxHIiIiIgBJmnSLSBmZbmYdMx2EiIiIiGQHlYiJpEdn4Gwz+xb4GTBCYVmbzIYlIiIiIpmgREwkPU7IdAAiIiIikj1UNVEkDdz9W2B34NTob/domIiIiIhUQUrERNLAzK4CngX2iP6eMbMrMhuViIiIiGSKqiaKpMcFQGd3/xnAzO4CJgL/zGhUIiIiIpIRKhETSQ8D8mLe50XDRERERKQKUomYSHr8B5hsZqOi9/2AJzIXjoiIiIhkkkrERMqZmVUDJgHnASujv/Pc/f4U5u1tZnPNbJ6Z3ZhgfG0zezEaP9nMWsSMuykaPtfMTogZ/nszm2Vm/zOz580spwxWU0RERERKQCViIuXM3beZ2cPufjgwPdX5zKw68DDQC1gMTDGz0e7+ZcxkFwCr3H1/MxsI3AWcYWaHAAOBQ4G9gPfM7EDgl8CVwCHuvtHMXoqmG77TKyoiIiIiKVOJmEh6vG9mp5lZSdqFdQLmufsCd98CvAD0jZumL/Bk9HoEcFz0GX2BF9x9s7t/A8yLlgfhAcwuZlYDqAMsLd0qiYiIiEhpKRETSY+LgZeBzWa21szWmdnaJPM0BRbFvF8cDUs4jbvnAmuARkXN6+5LgH8A3wHfA2vc/d3SrZKIiIiIlJYSMZFyFrUR6+3u1dy9lrvv5u713H23DMTSgFBa1pJQZXFXMzuniGkvMrOpZjZ1+fLl6QxTpOrZtAlyczMdhYiIpJESMZFy5u7bgIdKMesSoHnM+2bRsITTRFUN6wMripn3eOAbd1/u7luBV4CuRcQ91N07uHuHJk2alCJ8EUnJCy9Ao0bh75RT4P/+D6ZPh7y85POKiEiFpURMJD1K00ZsCnCAmbU0s1qETjVGx00zGhgcvT4dGOvuHg0fGPWq2BI4APiUUCWxi5nViWI5Dphd+tUSkVLLzYXrroMzz4TDDw//v/4arr0W2reHxo2hXz+4/36YORO2bct0xCIiUobUa6JIelwM/B7IM7NNhB9z9uKqJ7p7rpldDrwDVAeGufssM7sNmOruowm/Rfa0mc0jdIs/MJp3VtQj4pdALnCZu+cRfstsBKH3xlzgM2Bo+ayyiBRpxQoYOBDeew8uuwzuvRdq1Qrjli6F8eNh3Ljw99prYXijRtCjBxxzTPg75BAo0bMdERHJJhYenotIeYraiZ0NtHT328xsb2BPd5+c4dBS0qFDB586dWqmwxCpHGbOhP79YckSeOQROP/84qdftCgkZPnJ2cKFYXiTJtCzZ0FidtBBSsyyiJlNc/cOmY5DRLKXEjGRNDCzR4BtwLHu3irqNONdd++Y4dBSokRMpIy8+GJIvHbfHV55BTp3LvkyFi4sKC0bNw4WLw7Df/nLwonZ/vsrMcsgJWIikowSMZE0MLPp7n6EmX0W/bAzZjbT3dtmOrZUKBET2Ul5efDHP8Ldd0O3bjBiREicdpY7LFhQODH7/vswrmnTgqSsZ09o2VKJWRopERORZNRGTCQ9tppZdcABzKwJoYRMRCq7lStDRxzvvguXXAIPPFDQHmxnmcF++4W/Cy8MidlXXxVUY3z3XXjmmTDt3nsXJGbHHBPei4hIxqhETCQNzOxs4AzgCOBJQg+HN7v7yxkNLEUqERMppS++CD0fLl4MDz8ckqV0cofZswtKy8aPDx2FAOy7b+GqjE3jfy9edoZKxEQkGSViImliZgcTuos34H13rzDdxisREymFl1+GIUOgfv3QHqxLl0xHFLrAnzWrIDH74ANYtSqMO+CAgmqMxxxTNlUnqzAlYiKSjBIxEUlKiZhICeTlwc03w9//Dl27hvZge+6Z6agSy8uDzz8vKC374ANYuzaMO/jgwm3M9MPuJaJETESSUSImIkkpERNJ0apVcNZZ8PbbcPHF8OCDZdceLB3y8uCzzwpKzD76CNavD+MOPbQgMevRI/yumRRJiZiIJKNETESSUiImkoL//S+0B/vuO3joIbjookxHtPO2boXp0wsSs48/hg0bQichbdoUlJb16BG65JftlIiJSDJKxEQkKSViIkmMHAmDB0O9euF1166Zjqh8bNkCU6YU9Mr4ySewaVNIzA4/vKDE7KijYLfdMh1tRikRE5FklIiJSFJKxESKkJcHt9wC/+//hc44Ro6EvfbKdFTps3kzTJ5cUGI2cWJI1qpXh/btCzr+6N4d6tbNdLRppURMRJJRIiYiSSkRE0lg9erQHuytt0K39A89BLVrZzqqzNq4ESZNKkjMJk8O1Rtr1ICOHQuqMnbrBnXqZDracqVETESSUSImIkkpEROJ8+WXoT3YwoWhQ46LLw7V86Swn3+GCRMKqjJOmQK5uVCzJnTuXFCV8cgjIScn09GWKSViIpKMEjERSUqJmEiMUaNg0CDYdddQFbFbt0xHVHGsXx86/MjvLn/q1PDbZrVrh6qd+YlZ584VvnRRiZiIJKNETESSUiImQkgY/vIXuOMO6NQp/Ehz06aZjqpiW7s2dJGfX5Xxs8/AHXbZJXR4kl+VsWPHivUzACgRE5HklIiJSFJKxKTKW70azjkH3ngDzj8fHn640lWlywqrVsGHHxZUZZw5MwyvUyd0+JFfYta+fWh3lsWUiIlIMkrERCQpJWJSpc2eDX37wjffwAMPwKWXqj1YuqxYAR98UFBiNmtWGF63bugiPz8xO/zw0FNjFlEiJiLJZPfjJBERkUx67TU499xQVW7s2HDzL+nTqBH8+tfhD2DZssKJ2VtvheH168PRRxd0l9+2LVSrlrGwRURSoRIxEUlKJWJS5WzbBn/9K9x2W2ifNHIkNG+e6agk3vffh2qM+VUZv/46DG/QAHr0KCgxO/TQtCdmKhETkWSUiIlIUkrEpEpZsyaUgr3+OgwZAo88ovZgFcXixQVJ2fjxsGBBGN64ceHErFWrcq9eqkRMRJJRIiYiSSkRkypjzpzw+2Dz58N998Fll6k9WEX27bcFidm4cfDdd2H4HnsUVGM85hg48MAy385KxEQkGSViIpKUEjGpEkaPDj0j5uTAiBGhzZFUHu6hw5XYxGzJkjBuzz0LkrKePWG//XY6MVMiJiLJKBETkaSUiEmltm0b3H473Hpr6BZ91Ci1B6sK3GHevIKkbPx4+OGHMK5Zs4LE7JxzoGbNEi9eiZiIJKNETESSUiImldbatTBoUOgdcdAgePTR0EOiVD3uMHdu4cTMHX78sVQdfSgRE5Fk1H29iIhUTXPnhvZgX38dfh/siivUHqwqM4ODDw5/l14akrAlS9QNvoiUGyViIiJS9YwZA2efDbVqwXvvhXZBIrHMQhVFEZFyosc8IiJSdWzbBnfcAX36wP77w9SpSsJERCQjVCImIiJVw7p1MHhw6IzjnHNg6FC1BxMRkYxRIiYiIpXf11+H9mBz54bfB7vqKrUHExGRjFIiJiIildubb8JZZ0GNGvDuu3DssZmOSERERG3ERESkknKHO++EU06Bli1DezAlYSIikiVUIiYiIpXP+vUwZAiMHBlKwx57DOrUyXRUIiIi2ykRExGRymXevNAebPZs+L//g9//Xu3BREQk6ygRExGRyuPtt+HMM8OP8L7zDhx/fKYjEhERSUhtxEREpOJzh7//HU46CfbZJ7QHUxImIiJZTCViIiJSsa1fD+efDy+/DAMHwhNPqD2YiIhkPSViIiJScc2fH9qDffkl3HMPXHON2oOJiEiFoKqJIlnMzHqb2Vwzm2dmNyYYX9vMXozGTzazFjHjboqGzzWzE2KG725mI8xsjpnNNrMj07Q6ImXrnXegY0dYsiS0Dbv2WiVhIiJSYSgRE8lSZlYdeBg4ETgEONPMDomb7AJglbvvD9wH3BXNewgwEDgU6A38K1oewAPA2+5+MNAWmF3e6yJSptzh7rtDe7BmzUJ7sF69Mh2ViIhIiSgRE8lenYB57r7A3bcALwB946bpCzwZvR4BHGdmFg1/wd03u/s3wDygk5nVB44GngBw9y3uvrr8V0WkjPz8c2gHdsMNcPrpMHEi7LtvpqMSEREpMSViItmrKbAo5v3iaFjCadw9F1gDNCpm3pbAcuA/ZvaZmT1uZrsm+nAzu8jMpprZ1OXLl5fF+ojsnAULoGtXGDEC7roLXngBdk24+4qIiGQ9JWIiVUsN4AjgEXc/HPgZ2KHtGYC7D3X3Du7eoUmTJumMUWRH//1vaA/23Xfw5ptw/fVqDyYiIhWaEjGR7LUEaB7zvlk0LOE0ZlYDqA+sKGbexcBid58cDR9BSMxEspM7/OMf0Ls37LVXaA92wgnJ5xMREclySsREstcU4AAza2lmtQidb4yOm2Y0MDh6fTow1t09Gj4w6lWxJXAA8Km7/wAsMrODonmOA74s7xURKZUNG+Dss+G66+DXvw7twfbbL9NRiYiIlAn9jphIlnL3XDO7HHgHqA4Mc/dZZnYbMNXdRxM63XjazOYBKwnJGtF0LxGSrFzgMnfPixZ9BfBslNwtAM5L64qJpOKbb6B/f/j8c/jb30LnHKqKKCIilYiFh+ciIkXr0KGDT506NdNhSFXx/vvwm9/Atm3w3HNw4omZjkikxMxsmrt3yHQcIpK9VDVRRESygzvcey/86lew554wZYqSMBERqbSUiImISOZt2ADnnAPXXAP9+oX2YPvvn+moREREyo0SMRERyaxvv4Vu3eD55+HOO8PvhNWrl+moREREypU66xARkcwZOza0B8vNhTFj4KSTMh2RiIhIWqhETERE0s8d7r8/tAfbY4/QHkxJmIiIVCFKxEREJL02boRBg+D3v4c+fWDyZDjggExHJSIiklZKxEREJH2++w66d4dnn4Xbb1d7MBERqbLURkxERNJj/HgYMAC2bIHRo+GUUzIdkYiISMaoRExERMqXOzz4IBx/PDRuDJ9+qiRMRESqPCViIiJSfjZuhCFD4KqrQvI1eTIcdFCmoxIREck4JWIiIlI+Fi2Co46Cp56CW2+FV16B3XbLdFQiIiJZQW3ERESk7H3wQWgPtmkTvPZa6B1RREREtlOJmIiIlB13eOih0B6sYcPQHkxJmIiIyA6UiImISNnYtAnOPx+uuAJOPDG0Bzv44ExHJSIikpWUiImIyM5bvBiOPhqGD4e//AVefRXq1890VCIiIllLbcRERGTnfPQRnH46bNgAo0ZBv36ZjkhERCTrqURMRERKxx3+9S849ljYfffQHkxJmIiISEqUiImISMlt2gQXXgiXXQYnnBCSsFatMh2ViIhIhaFETERESmbJEujRA4YNg5tvhtGj1R5MRESkhNRGTEREUvfJJ3DaafDzzzByJPz615mOSEREpEJSiZiIiCTnDo8+CsccA/Xqha7plYSJiIiUmhIxEREp3ubNcNFFcOml4Yeap0yBQw7JdFQiIiIVmhIxEREp2tKl0LMnPP44/PGP8PrroYdEERER2SlqIyYiIolNmBDag61bByNGhNciIiJSJlQiJiIiOxo6NJSE7borTJqkJExERKSMKRETEZECmzfDxReHv2OPDe3BDjss01GJiIhUOkrEREQk+P770Cvi0KFw443wxhvQoEGmoxIREamU1EZMRERC9cNf/xrWrIGXXoIBAzIdkYiISKWmEjERkaru8cehRw/YZZeQkCkJExERKXdKxEREqqotW8Jvg/32t6FjjilToHXrTEclIiJSJSgRExGpin74IXTG8eijcP318Oab0LBhpqMSERGpMtRGTESkqpk8OXRHv2oVvPACnHFGpiMSERGpcpSIiYhURlu3wtKl8N134e/bbwtev/8+7LVX+MHmtm0zHamIiEiVpERMRKQiWr26ILGK/ctPuJYuhW3bCs/TuDHsvTeccw7cfTc0apSR0EVERESJmIhI9snNLbo0K/9v7drC89SsCc2bwz77wHHHhYQr/2+ffcK4OnUysz4iIiKyAyViIlnMzHoDDwDVgcfd/e9x42sDTwHtgRXAGe6+MBp3E3ABkAdc6e7vxMxXHZgKLHH3U9KwKhJr7drEyVX+sCVLdizNatQoJFX77Rd+dDk/wcpPtn7xC6im/pdEREQqCiViIlkqSpYeBnoBi4EpZjba3b+MmewCYJW7729mA4G7gDPM7BBgIHAosBfwnpkd6O550XxXAbOB3dK0OlVHbi58/33xpVlr1hSep0aNUGK1994FSVb83667ZmZ9REREpFwoERPJXp2Aee6+AMDMXgD6ArGJWF/g1uj1COAhM7No+Avuvhn4xszmRcubaGbNgJOBO4E/pGNFKpV16xInV/lJ15IlkJdXeJ6GDUMy1bJl+OHkRKVZ1atnZn1EREQkI5SIiWSvpsCimPeLgc5FTePuuWa2BmgUDZ8UN2/T6PX9wPVAveI+3MwuAi4C2HvvvUu1AhVOXl7y0qzVqwvPU6MGNGsWEqqjjy6cYO29dyjpqlfsVy0iIiJVkBIxkSrEzE4Blrn7NDPrWdy07j4UGArQoUMHL//o0mDdusQlWbGlWbm5hefZffeQXO2zDxx1VOEOMPbeG375S5VmiYiISIkpERPJXkuA5jHvm0XDEk2z2MxqAPUJnXYUNW8foI+ZnQTkALuZ2TPufk75rEIa5eXBDz8k7vwi/2/VqsLzVK9eUJrVvXvi0qzd1IxOREREyp4SMZHsNQU4wMxaEpKogcBZcdOMBgYDE4HTgbHu7mY2GnjOzO4ldNZxAPCpu08EbgKISsSurTBJ2Pr1sGhR0b0NLl68Y2lW/foFyVW3bjt26b7nnirNEhERkYxQIiaSpaI2X5cD7xC6rx/m7rPM7DZgqruPBp4Ano4641hJSNaIpnuJ0LFHLnBZTI+J2WfbtuSlWStXFp6nenVo2jQkVV277tgBRvPmIRETERERyULmXjmafohI+enQoYNPnTq19Av4+eeiS7O++y6M27q18Dz16hW0zUrUnftee4WOMkREspCZTXP3DpmOQ0Syl+5iRKR8uIeSqq+/hhUrCo+rVq2gNKtzZxgwYMcSLZVmiYiISCWmRExEyocZHHIItGuXuDSrZs1MRygiIiKSMUrERKT8PPFEpiMQERERyUrVMh2AiIiIiIhIVaNETEREREREJM2UiImIiIiIiKSZEjEREREREZE0UyImIiIiIiKSZkrERERERERE0kyJmIiIiIiISJopERMREREREUkzc/dMxyAiWc7MlgPflnL2xsBPZRhOJlWWdaks6wFal2xUWdYDdm5d9nH3JmUZjIhULkrERKRcmdlUd++Q6TjKQmVZl8qyHqB1yUaVZT2gcq2LiGQfVU0UERERERFJMyViIiIiIiIiaaZETETK29BMB1CGKsu6VJb1AK1LNqos6wGVa11EJMuojZiIiIiIiEiaqURMREREREQkzZSIiYiIiIiIpJkSMRHZaWY2zMyWmdn/ihhvZvagmc0zs8/N7Ih0x5iqFNalp5mtMbMZ0d8t6Y4xFWbW3MzGmdmXZjbLzK5KME2F2C4prktF2S45Zvapmc2M1uWvCaapbWYvRttlspm1yECoxUpxPYaY2fKYbXJhJmJNlZlVN7PPzGxMgnFZv01EpOKpkekARKRSGA48BDxVxPgTgQOiv87AI9H/bDSc4tcF4CN3PyU94ZRaLnCNu083s3rANDP7r7t/GTNNRdkuqawLVIztshk41t3Xm1lN4GMze8vdJ8VMcwGwyt33N7OBwF3AGZkIthiprAfAi+5+eQbiK42rgNnAbgnGVYRtIiIVjErERGSnufuHwMpiJukLPOXBJGB3M9szPdGVTArrUiG4+/fuPj16vY5wg9k0brIKsV1SXJcKIfqu10dva0Z/8b1m9QWejF6PAI4zM0tTiClJcT0qDDNrBpwMPF7EJFm/TUSk4lEiJiLp0BRYFPN+MRX0RjpyZFQl6y0zOzTTwSQTVaM6HJgcN6rCbZdi1gUqyHaJqsDNAJYB/3X3IreLu+cCa4BGaQ0yBSmsB8BpUbXXEWbWPL0Rlsj9wPXAtiLGV4htIiIVixIxEZGSmQ7s4+5tgX8Cr2Y2nOKZWV1gJHC1u6/NdDw7I8m6VJjt4u557t4OaAZ0MrPDMhxSqaSwHq8DLdy9DfBfCkqUsoqZnQIsc/dpmY5FRKoWJWIikg5LgNin4c2iYRWOu6/Nr5Ll7m8CNc2scYbDSihquzMSeNbdX0kwSYXZLsnWpSJtl3zuvhoYB/SOG7V9u5hZDaA+sCKtwZVAUevh7ivcfXP09nGgfZpDS1U3oI+ZLQReAI41s2fipqlQ20REKgYlYiKSDqOBQVEvfV2ANe7+faaDKg0z+2V+2xAz60Q4j2bdDVkU4xPAbHe/t4jJKsR2SWVdKtB2aWJmu0evdwF6AXPiJhsNDI5enw6Mdfesan+VynrEtTfsQ2jbl3Xc/SZ3b+buLYCBhO/7nLjJsn6biEjFo14TRWSnmdnzQE+gsZktBv5CaLyPuz8KvAmcBMwDNgDnZSbS5FJYl9OBS80sF9gIDMzSG7JuwLnAF1E7HoA/AntDhdsuqaxLRdkuewJPmll1QrL4kruPMbPbgKnuPpqQdD5tZvMIHccMzFy4RUplPa40sz6EXi9XAkMyFm0pVMBtIiIVjGXndUpERERERKTyUtVEERERERGRNFMiJiIiIiIikmZKxERERERERNJMiZiIiIiIiEiaKRETERERERFJMyViIiKSdczsVjO7thTztTOzk3Z2OSIiIuVNiZiIiFQm7Qi/jSYiIpLVlIiJiEhWMLM/mdlXZvYxcFA0bD8ze9vMppnZR2Z2cDR8uJk9amZTo3lOMbNawG3AGWY2w8zOiBZ9iJmNN7MFZnZlZtZORESksBqZDkBERMTM2gMDCSVaNYDpwDRgKHCJu39tZp2BfwHHRrO1ADoB+wHjgP2BW4AO7n55tNxbgYOBY4B6wFwze8Tdt6ZlxURERIqgRExERLLBUcAod98AYGajgRygK/CymeVPVztmnpfcfRvwtZktICRcibzh7puBzWa2DPgFsLgc1kFERCRlSsRERCRbVQNWu3u7IsZ7kvf5Nse8zkPXPhERyQJqIyYiItngQ6Cfme1iZvWAU4ENwDdmNgDAgrYx8wwws2pmth+wLzAXWEeogigiIpLVlIiJiEjGuft04EVgJvAWMCUadTZwgZnNBGYBfWNm+w74NJr+EnffRGgrdkhcZx0iIiJZx9yLqskhIiKSncxsODDG3UdkOhYREZHSUImYiIiIiIhImqlETEREREREJM1UIiYiIiIiIpJmSsRERERERETSTImYiIiIiIhImikRExERERERSTMlYiIiIiIiImn2/wFhIJXQg3KFagAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "from spheres import *\n",
    "\n",
    "params = {\"n_qubits\": 1,\\\n",
    "          #\"depth\": 2,\\\n",
    "          \"n_copies\": 3,\\\n",
    "          \"every\": 2,\\\n",
    "          \"pairwise\": True,\\\n",
    "          \"noise_model\": thermal_noise_model,\\\n",
    "          \"noise_model_name\": thermal_noise_model,\\\n",
    "          \"error_on_all\": True,\\\n",
    "          \"backend\": Aer.get_backend(\"qasm_simulator\"),\\\n",
    "          \"shots\": 8000}\n",
    "\n",
    "depths = list(range(1, 5))\n",
    "n_runs = 50\n",
    "experiments = []\n",
    "for depth in depths:\n",
    "    print(\"depth: %d\" % depth )\n",
    "    circ_errors, sym_circ_errors = [], []\n",
    "    for i in range(n_runs):\n",
    "        experiment = eval_pytket_symmetrization_performance(**params, depth=depth)\n",
    "        circ_errors.append(experiment[\"circ_error\"])\n",
    "        sym_circ_errors.append(experiment[\"sym_circ_error\"])\n",
    "    experiments.append({**params, **{\"depth\": depth,\\\n",
    "                                     \"circ_error\": sum(circ_errors)/len(circ_errors),\\\n",
    "                                     \"sym_circ_error\": sum(sym_circ_errors)/len(sym_circ_errors)}})\n",
    "\n",
    "plot_symmetrization_performance(\"depth\", experiments)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Some results, picked somewhat arbitrarily:\n",
    "\n",
    "![result1](img/sym_circ_result1.png)\n",
    "\n",
    "![result2](img/sym_circ_result2.png)\n",
    "\n",
    "![result3](img/sym_circ_result3.png)\n",
    "\n",
    "![result4](img/sym_circ_result4.png)\n",
    "\n",
    "<img style=\"display:block; margin:auto\" src=\"img/sym_circ_result5.png\" width=450>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr>\n",
    "\n",
    "Clearly, it's not a given that this method of \"stablization\" will deliver results at least at the moment. It would be interesting to find some set of conditions that would make this algorithm relevant. The difficulties are, of course, errors in the symmetrization procedure itself, and also the scaling of the auxilliary qubits. Unless intermediate measurements are possible, allowing one to reuse auxilliary qubits, the cost of applying the algorithm quickly becomes prohibitive--and there aren't enough gains at the small scale, it seems, to make it worth using. In the future, however, one could imagine having thousands of qubits, as well as reliable CSWAP gates, so that one could reliably symmetrize over hundreds copies of the circuit, etc. See the attached papers for some discussion of the types of errors that this method is best at dealing with."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr>\n",
    "\n",
    "## Bibliography \n",
    "\n",
    "[Error Symmetrization in Quantum Computers](https://arxiv.org/abs/quant-ph/9611046)\n",
    "\n",
    "[Stabilization of Quantum Computations By Symmetrization](https://heyredhat.github.io/refs/Stabilization_of_Quantum_Compu.pdf)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "VPython",
   "language": "python",
   "name": "vpython"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}