toric-code-example.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Toric code example\n",
    "\n",
    "In this example, we'll use PyMatching to estimate the threshold of the toric code \n",
    "under an independent noise model with perfect syndrome measurements.\n",
    "The decoding problem for the toric code is identical for $X$-type and $Z$-type \n",
    "errors, so we will only simulate decoding $Z$-type errors using $X$-type stabilisers in this example.\n",
    "\n",
    "First, we will construct a check matrix $H_X$ corresponding to the $X$-type stabilisers.\n",
    "Each element $H_X[i,j]$ will be 1 if the $i$th $X$ stabiliser acts non-trivially on \n",
    "the $j$th qubit, and is 0 otherwise.\n",
    "\n",
    "\n",
    "We will construct $H_X$ by taking the [hypergraph product](https://arxiv.org/abs/0903.0566) of two repetition codes.\n",
    "The hypergraph product code construction $HGP(H_1,H_2)$ takes as input the parity check \n",
    "matrices of two linear codes $C_1:=\\ker H_1$ and $C_2:= \\ker H_2$. The code $HGP(H_1,H_2)$ \n",
    "is a CSS code with the check matrix for the $X$ stabilisers given by\n",
    "\n",
    "$$H_X=[H_1\\otimes I_{n_2},I_{r_1}\\otimes H_2^T]$$\n",
    "\n",
    "and with the check matrix for the $Z$ stabilisers given by \n",
    "\n",
    "$$H_Z=[I_{n_1}\\otimes H_2,H_1^T\\otimes I_{r_2}]$$\n",
    "\n",
    "where $H_1$ has dimensions $r_1\\times n_1$, $H_2$ has dimensions $r_2\\times n_2$ and $I_l$ denotes the $l\\times l$ identity matrix.\n",
    "\n",
    "Since we only need the $X$ stabilisers of the toric code with lattice size L, we only need to construct $H_X$, using the check matrix of a repetition code with length L for both $H_1$ and $H_2$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.sparse import hstack, kron, eye, csr_matrix, block_diag\n",
    "\n",
    "\n",
    "def repetition_code(n):\n",
    "    \"\"\"\n",
    "    Parity check matrix of a repetition code with length n.\n",
    "    \"\"\"\n",
    "    row_ind, col_ind = zip(*((i, j) for i in range(n) for j in (i, (i+1)%n)))\n",
    "    data = np.ones(2*n, dtype=np.uint8)\n",
    "    return csr_matrix((data, (row_ind, col_ind)))\n",
    "\n",
    "\n",
    "def toric_code_x_stabilisers(L):\n",
    "    \"\"\"\n",
    "    Sparse check matrix for the X stabilisers of a toric code with \n",
    "    lattice size L, constructed as the hypergraph product of \n",
    "    two repetition codes.\n",
    "    \"\"\"\n",
    "    Hr = repetition_code(L)\n",
    "    H = hstack(\n",
    "            [kron(Hr, eye(Hr.shape[1])), kron(eye(Hr.shape[0]), Hr.T)],\n",
    "            dtype=np.uint8\n",
    "        )\n",
    "    H.data = H.data % 2\n",
    "    H.eliminate_zeros()\n",
    "    return csr_matrix(H)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "From the [K&uuml;nneth theorem](https://en.wikipedia.org/wiki/K%C3%BCnneth_theorem), the $X$ logical operators of the toric code are given by\n",
    "\n",
    "$$L_X=\\left(\n",
    "  \\begin{array}{cc}\n",
    "  \\mathcal{H}^1\\otimes \\mathcal{H}^0 & 0 \\\\\n",
    "  0 & \\mathcal{H}^0\\otimes \\mathcal{H}^1\n",
    "\\end{array} \\right)$$\n",
    "\n",
    "where $\\mathcal{H}^0$ and $\\mathcal{H}^1$ are the zeroth and first cohomology groups of the length-one chain complex that has the repetition code parity check matrix as its boundary operator. We can construct this matrix with the following function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "def toric_code_x_logicals(L):\n",
    "    \"\"\"\n",
    "    Sparse binary matrix with each row corresponding to an X logical operator \n",
    "    of a toric code with lattice size L. Constructed from the \n",
    "    homology groups of the repetition codes using the Kunneth \n",
    "    theorem.\n",
    "    \"\"\"\n",
    "    H1 = csr_matrix(([1], ([0],[0])), shape=(1,L), dtype=np.uint8)\n",
    "    H0 = csr_matrix(np.ones((1, L), dtype=np.uint8))\n",
    "    x_logicals = block_diag([kron(H1, H0), kron(H0, H1)])\n",
    "    x_logicals.data = x_logicals.data % 2\n",
    "    x_logicals.eliminate_zeros()\n",
    "    #print(H1.todense())\n",
    "    #print(H0)\n",
    "    return csr_matrix(x_logicals)\n",
    "\n",
    "#toric_code_x_logicals(L=5)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that we have the $X$ check matrix and $X$ logicals of the toric code, we can use PyMatching to simulate its performance using the minimum-weight perfect matching decoder and an error model of our choice.\n",
    "\n",
    "To do so, we first import the Matching class from PyMatching, and use it to construct a Matching object from the check matrix of the stabilisers:\n",
    "```\n",
    "from pymatching import Matching\n",
    "matching=Matching(H)\n",
    "```\n",
    "Constructing the Matching object, while efficient, is often slower than the decoding step itself. As a result, it's best to construct the Matching object only at the beginning of the experiment, and not before every use of the decoder, in order to obtain the best performance.\n",
    "\n",
    "We also choose a number of trials, `num_trials`.  For each trial, we simulate a $Z$ error under an independent noise model, in which each qubit independently suffers a $Z$ error with probability $p$:\n",
    "```\n",
    "noise = np.random.binomial(1, p, H.shape[1])\n",
    "```\n",
    "Here, `noise` is a binary vector and `noise[i]` is 1 if qubit $i$ suffers a $Z$ error, and 0 otherwise.\n",
    "\n",
    "The syndrome of the $X$ stabilisers is then calculated from the dot product (modulo 2) with the $X$ check matrix $H$:\n",
    "```\n",
    "syndrome = H@noise % 2\n",
    "```\n",
    "\n",
    "We can now use PyMatching to infer the most probable individual error given the syndrome:\n",
    "```\n",
    "prediction = matching.decode(syndrome)\n",
    "```\n",
    "\n",
    "We use this to predict which logical X operators have been flipped:\n",
    "\n",
    "```\n",
    "predicted_logicals_flipped = logicals@prediction % 2\n",
    "```\n",
    "\n",
    "The actual logicals that were flipped are:\n",
    "\n",
    "```\n",
    "actual_logicals_flipped = logicals@prediction % 2\n",
    "```\n",
    "\n",
    "Our decoder was successful if `actual_logical_observables` equals `predicted_logical_observables`.\n",
    "\n",
    "Taken together, we obtain the following function `num_decoding_failures` that returns the number of logical errors after `num_trials` Monte Carlo trials, simulating an independent error model with error probability `p`, with the $X$ stabiliser check matrix `H` and $X$ logical matrix `logicals`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pymatching import Matching\n",
    "\n",
    "def num_decoding_failures_via_physical_frame_changes(H, logicals, p, num_trials):\n",
    "    matching = Matching.from_check_matrix(H, weights=np.log((1-p)/p))\n",
    "    num_errors = 0\n",
    "    for i in range(num_trials):\n",
    "        noise = np.random.binomial(1, p, H.shape[1])\n",
    "        syndrome = H@noise % 2\n",
    "        prediction = matching.decode(syndrome)\n",
    "        predicted_logicals_flipped = logicals@prediction % 2\n",
    "        actual_logicals_flipped = logicals@noise % 2\n",
    "        if not np.array_equal(predicted_logicals_flipped, actual_logicals_flipped):\n",
    "            num_errors += 1\n",
    "    return num_errors"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can speed this up slightly by telling PyMatching about the logical operators matrix when we create the `pymatching.Matching` object, \n",
    "using the `faults_matrix` argument. \n",
    "By doing this, `pymatching.Matching.decode` directly predicts which logicals have been flipped. This is a bit faster, as it allows the decoder to make some more optimisations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def num_decoding_failures(H, logicals, p, num_trials):\n",
    "    matching = Matching.from_check_matrix(H, weights=np.log((1-p)/p), faults_matrix=logicals)\n",
    "    num_errors = 0\n",
    "    for i in range(num_trials):\n",
    "        noise = np.random.binomial(1, p, H.shape[1])\n",
    "        syndrome = H@noise % 2\n",
    "        predicted_logicals_flipped = matching.decode(syndrome)\n",
    "        actual_logicals_flipped = logicals@noise % 2\n",
    "        if not np.array_equal(predicted_logicals_flipped, actual_logicals_flipped):\n",
    "            num_errors += 1\n",
    "    return num_errors"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using this function, we can now estimate the threshold of the toric code by varying the error rate $p$, for a range of lattice sizes $L$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Simulating L=11...Simulating L=17...Simulating L=9...Simulating L=13...\n",
      "\n",
      "Simulating L=15...Simulating L=19...\n",
      "\n",
      "\n",
      "\n",
      "Finish simulating L=9...\n",
      "Finish simulating L=11...\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Process ForkPoolWorker-6:\n"
     ]
    },
    {
     "ename": "KeyboardInterrupt",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
      "File \u001b[0;32m<timed exec>:30\u001b[0m\n",
      "File \u001b[0;32m/opt/conda/lib/python3.10/multiprocessing/pool.py:367\u001b[0m, in \u001b[0;36mPool.map\u001b[0;34m(self, func, iterable, chunksize)\u001b[0m\n\u001b[1;32m    362\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mmap\u001b[39m(\u001b[38;5;28mself\u001b[39m, func, iterable, chunksize\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mNone\u001b[39;00m):\n\u001b[1;32m    363\u001b[0m     \u001b[38;5;124;03m'''\u001b[39;00m\n\u001b[1;32m    364\u001b[0m \u001b[38;5;124;03m    Apply `func` to each element in `iterable`, collecting the results\u001b[39;00m\n\u001b[1;32m    365\u001b[0m \u001b[38;5;124;03m    in a list that is returned.\u001b[39;00m\n\u001b[1;32m    366\u001b[0m \u001b[38;5;124;03m    '''\u001b[39;00m\n\u001b[0;32m--> 367\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_map_async\u001b[49m\u001b[43m(\u001b[49m\u001b[43mfunc\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43miterable\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mmapstar\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mchunksize\u001b[49m\u001b[43m)\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mget\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\n",
      "File \u001b[0;32m/opt/conda/lib/python3.10/multiprocessing/pool.py:768\u001b[0m, in \u001b[0;36mApplyResult.get\u001b[0;34m(self, timeout)\u001b[0m\n\u001b[1;32m    767\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mget\u001b[39m(\u001b[38;5;28mself\u001b[39m, timeout\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mNone\u001b[39;00m):\n\u001b[0;32m--> 768\u001b[0m     \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mwait\u001b[49m\u001b[43m(\u001b[49m\u001b[43mtimeout\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m    769\u001b[0m     \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39mready():\n\u001b[1;32m    770\u001b[0m         \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mTimeoutError\u001b[39;00m\n",
      "File \u001b[0;32m/opt/conda/lib/python3.10/multiprocessing/pool.py:765\u001b[0m, in \u001b[0;36mApplyResult.wait\u001b[0;34m(self, timeout)\u001b[0m\n\u001b[1;32m    764\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mwait\u001b[39m(\u001b[38;5;28mself\u001b[39m, timeout\u001b[38;5;241m=\u001b[39m\u001b[38;5;28;01mNone\u001b[39;00m):\n\u001b[0;32m--> 765\u001b[0m     \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_event\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mwait\u001b[49m\u001b[43m(\u001b[49m\u001b[43mtimeout\u001b[49m\u001b[43m)\u001b[49m\n",
      "File \u001b[0;32m/opt/conda/lib/python3.10/threading.py:607\u001b[0m, in \u001b[0;36mEvent.wait\u001b[0;34m(self, timeout)\u001b[0m\n\u001b[1;32m    605\u001b[0m signaled \u001b[38;5;241m=\u001b[39m \u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_flag\n\u001b[1;32m    606\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m signaled:\n\u001b[0;32m--> 607\u001b[0m     signaled \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43m_cond\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mwait\u001b[49m\u001b[43m(\u001b[49m\u001b[43mtimeout\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m    608\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m signaled\n",
      "File \u001b[0;32m/opt/conda/lib/python3.10/threading.py:320\u001b[0m, in \u001b[0;36mCondition.wait\u001b[0;34m(self, timeout)\u001b[0m\n\u001b[1;32m    318\u001b[0m \u001b[38;5;28;01mtry\u001b[39;00m:    \u001b[38;5;66;03m# restore state no matter what (e.g., KeyboardInterrupt)\u001b[39;00m\n\u001b[1;32m    319\u001b[0m     \u001b[38;5;28;01mif\u001b[39;00m timeout \u001b[38;5;129;01mis\u001b[39;00m \u001b[38;5;28;01mNone\u001b[39;00m:\n\u001b[0;32m--> 320\u001b[0m         \u001b[43mwaiter\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43macquire\u001b[49m\u001b[43m(\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m    321\u001b[0m         gotit \u001b[38;5;241m=\u001b[39m \u001b[38;5;28;01mTrue\u001b[39;00m\n\u001b[1;32m    322\u001b[0m     \u001b[38;5;28;01melse\u001b[39;00m:\n",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m: "
     ]
    }
   ],
   "source": [
    "%%time\n",
    "\n",
    "def simulate(L:int):\n",
    "    print(\"Simulating L={}...\".format(L))\n",
    "    Hx = toric_code_x_stabilisers(L)\n",
    "    logX = toric_code_x_logicals(L)\n",
    "    log_errors = []\n",
    "    for p in ps:\n",
    "        #print(p,)\n",
    "        num_errors = num_decoding_failures(Hx, logX, p, num_trials)\n",
    "        log_errors.append(num_errors/num_trials)\n",
    "#    log_errors_all_L.append(np.array(log_errors))\n",
    "    print(\"Finish simulating L={}...\".format(L))\n",
    "    return log_errors\n",
    "\n",
    "num_trials = 3000000\n",
    "Ls = range(9,21,2) #change from 8 to 9\n",
    "ps = np.linspace(0.01, 0.2, 9)\n",
    "ps = np.linspace(0.03, 0.15, 5)\n",
    "ps=[0.1,\n",
    " 0.07142857142857144,\n",
    " 0.051020408163265314,\n",
    " 0.036443148688046656,\n",
    " 0.0260308204914619,\n",
    " 0.018593443208187073,\n",
    " 0.013281030862990767,\n",
    " 0.009486450616421976]\n",
    "np.random.seed(2)\n",
    "log_errors_all_L = []\n",
    "from multiprocessing import Pool\n",
    "with Pool(processes=len(Ls)) as pool:\n",
    "        log_errors_all_L = pool.map(simulate,Ls)\n",
    "print(log_errors_all_L)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, let's plot the results! We expect to see a threshold of around 10.3%, although a precise estimate requires using more trials, larger lattice sizes and scanning more values of $p$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {},
   "outputs": [
    {
     "ename": "TypeError",
     "evalue": "unsupported operand type(s) for -: 'int' and 'list'",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mTypeError\u001b[0m                                 Traceback (most recent call last)",
      "Cell \u001b[0;32mIn [72], line 6\u001b[0m\n\u001b[1;32m      4\u001b[0m plt\u001b[38;5;241m.\u001b[39myscale(\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mlog\u001b[39m\u001b[38;5;124m'\u001b[39m),plt\u001b[38;5;241m.\u001b[39mxscale(\u001b[38;5;124m'\u001b[39m\u001b[38;5;124mlog\u001b[39m\u001b[38;5;124m'\u001b[39m)\n\u001b[1;32m      5\u001b[0m \u001b[38;5;28;01mfor\u001b[39;00m L, logical_errors \u001b[38;5;129;01min\u001b[39;00m \u001b[38;5;28mzip\u001b[39m(Ls, log_errors_all_L):\n\u001b[0;32m----> 6\u001b[0m     std_err \u001b[38;5;241m=\u001b[39m (logical_errors\u001b[38;5;241m*\u001b[39m(\u001b[38;5;241;43m1\u001b[39;49m\u001b[38;5;241;43m-\u001b[39;49m\u001b[43mlogical_errors\u001b[49m)\u001b[38;5;241m/\u001b[39mnum_trials)\u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39m\u001b[38;5;241m0.5\u001b[39m\n\u001b[1;32m      7\u001b[0m     plt\u001b[38;5;241m.\u001b[39merrorbar(ps, logical_errors, yerr\u001b[38;5;241m=\u001b[39mstd_err, label\u001b[38;5;241m=\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mL=\u001b[39m\u001b[38;5;132;01m{}\u001b[39;00m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;241m.\u001b[39mformat(L)\u001b[38;5;241m+\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124m,n=\u001b[39m\u001b[38;5;132;01m{}\u001b[39;00m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;241m.\u001b[39mformat(\u001b[38;5;241m2\u001b[39m\u001b[38;5;241m*\u001b[39mL\u001b[38;5;241m*\u001b[39mL))\n\u001b[1;32m      8\u001b[0m plt\u001b[38;5;241m.\u001b[39mplot(ps, ps, \u001b[38;5;124m'\u001b[39m\u001b[38;5;124m--\u001b[39m\u001b[38;5;124m'\u001b[39m,label\u001b[38;5;241m=\u001b[39m\u001b[38;5;124m'\u001b[39m\u001b[38;5;124my=x\u001b[39m\u001b[38;5;124m'\u001b[39m)  \u001b[38;5;66;03m#add y=x by weilei\u001b[39;00m\n",
      "\u001b[0;31mTypeError\u001b[0m: unsupported operand type(s) for -: 'int' and 'list'"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAicAAAGhCAYAAAC6URSFAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAAARRklEQVR4nO3dT2iUd/7A8c+oGC+Z2aZCaP4oHkpxWhCMIyh4yCWSgkKl0JN40IOwsEgOheKhUISeKh4cBemhuzdhoR6KIDnIKgjbrdVDj4IQaypiYWdiDpXq8zssDb+stptpZ/J8Jnm9YA7zzDjPJzHwffOdf5WiKIoAAEhiQ9kDAAD8f+IEAEhFnAAAqYgTACAVcQIApCJOAIBUxAkAkMqmsgfo1IsXL2J+fj4GBwejUqmUPQ4AsAJFUcTCwkKMjIzEhg2/vTfSd3EyPz8f4+PjZY8BAPwODx48iLGxsd+8T9/FyeDgYET854erVqslTwMArES73Y7x8fGldfy39F2c/PJUTrVaFScA0GdW8pIML4gFAFIRJwBAKuIEAEhFnAAAqYgTACAVcQIApCJOAIBUxAkAkIo4AQBSKSVO3nvvvXjttdfi/fffL+P0AEBipcTJX/7yl/jb3/5WxqkBgORKiZPJyckVffEPALD+dBwnN27ciEOHDsXIyEhUKpW4cuXKS/e5cOFC7NixI7Zs2RITExNx8+bNbswKAKwDHcfJ4uJi7Nq1K86fP//K2y9fvhynTp2K06dPx507d+LAgQMxPT0dc3Nzv2vAn376Kdrt9rILALB2dRwn09PTcebMmThy5Mgrbz979mwcP348Tpw4ETt37oxz587F+Ph4XLx48XcN+Omnn0atVlu6jI+P/67HAQD6Q1dfc/Ls2bO4fft2TE1NLTs+NTUVt27d+l2P+dFHH0Wr1Vq6PHjwoBujAgBJbermgz158iSeP38ew8PDy44PDw/Ho0ePlq4fPHgwvv3221hcXIyxsbH48ssvo9FovPIxBwYGYmBgoJtjAgCJdTVOflGpVJZdL4pi2bFr16714rQAwBrQ1ad1tm7dGhs3bly2SxIR8fjx45d2UwAAXqWrcbJ58+aYmJiI2dnZZcdnZ2dj//793TwVALBGdfy0ztOnT+PevXtL1+/fvx93796NoaGh2LZtW8zMzMTRo0djz549sW/fvrh06VLMzc3FyZMnuzo4ALA2dRwn33zzTUxOTi5dn5mZiYiIY8eOxRdffBEffPBB/Pjjj/HJJ5/EDz/8EO+8805cvXo1tm/f3r2pAYA1q1IURVH2EJ1ot9tRq9Wi1WpFtVotexwAYAU6Wb9L+W4dAIBfI04AgFTECQCQijgBAFIRJwBAKuIEAEilb+Kk2WxGvV7/1S8IBADWBp9zAgD0nM85AQD6ljgBAFIRJwBAKuIEAEhFnAAAqYgTACAVcQIApCJOAIBUxAkAkIo4AQBSEScAQCriBABIRZwAAKmIEwAgFXECAKTSN3HSbDajXq9Ho9EoexQAoIcqRVEUZQ/RiXa7HbVaLVqtVlSr1bLHAQBWoJP1u292TgCA9UGcAACpiBMAIBVxAgCkIk4AgFTECQCQijgBAFIRJwBAKuIEAEhFnAAAqYgTACAVcQIApCJOAIBUxAkAkIo4AQBSEScAQCriBABIRZwAAKn0TZw0m82o1+vRaDTKHgUA6KFKURRF2UN0ot1uR61Wi1arFdVqtexxAIAV6GT97pudEwBgfRAnAEAq4gQASEWcAACpiBMAIBVxAgCkIk4AgFTECQCQijgBAFIRJwBAKuIEAEhFnAAAqYgTACAVcQIApCJOAIBUxAkAkIo4AQBSEScAQCriBABIRZwAAKn0TZw0m82o1+vRaDTKHgUA6KFKURRF2UN0ot1uR61Wi1arFdVqtexxAIAV6GT97pudEwBgfRAnAEAq4gQASEWcAACpiBMAIBVxAgCkIk4AgFTECQCQijgBAFIRJwBAKuIEAEhFnAAAqYgTACAVcQIApCJOAIBUxAkAkIo4AQBSEScAQCriBABIRZwAAKmIEwAgFXECAKQiTgCAVMQJAJCKOAEAUumbOGk2m1Gv16PRaJQ9CgDQQ5WiKIqyh+hEu92OWq0WrVYrqtVq2eMAACvQyfrdNzsnAMD6IE4AgFTECQCQijgBAFIRJwBAKuIEAEhFnAAAqYgTACAVcQIApCJOAIBUxAkAkIo4AQBSEScAQCriBABIRZwAAKmIEwAgFXECAKQiTgCAVMQJAJCKOAEAUhEnAEAq4gQASEWcAACpiBMAIBVxAgCkIk4AgFTECQCQijgBAFIRJwBAKuIEAEhFnAAAqfRNnDSbzajX69FoNMoeBQDooUpRFEXZQ3Si3W5HrVaLVqsV1Wq17HEAgBXoZP3um50TAGB9ECcAQCriBABIRZwAAKmIEwAgFXECAKQiTgCAVMQJAJCKOAEAUhEnAEAq4gQASEWcAACpiBMAIBVxAgCkIk4AgFTECQCQijgBAFIRJwBAKuIEAEhFnAAAqYgTACAVcQIApCJOAIBUxAkAkIo4AQBSEScAQCriBABIRZwAAKmIEwAgFXECAKQiTgCAVMQJAJCKOAEAUhEnAEAq4gQASEWcAACpiBMAIBVxAgCkIk4AgFTECQCQSt/ESbPZjHq9Ho1Go+xRAIAeqhRFUZQ9RCfa7XbUarVotVpRrVbLHgcAWIFO1u++2TkBANYHcQIApCJOAIBUxAkAkIo4AQBSEScAQCriBABIRZwAAKmIEwAgFXECAKQiTgCAVMQJAJCKOAEAUhEnAEAq4gQASEWcAACpiBMAIBVxAgCkIk4AgFTECQCQijgBAFIRJwBAKuIEAEhFnAAAqYgTACAVcQIApCJOAIBUxAkAkIo4AQBSEScAQCriBABIRZwAAKmIEwAgFXECAKQiTgCAVMQJAJCKOAEAUhEnAEAq4gQASEWcAACpiBMAIBVxAgCkIk4AgFTECQCQijgBAFIRJwBAKuIEAEhFnAAAqYgTACAVcQIApCJOAIBUxAkAkIo4AQBSEScAQCriBABIRZwAAKmIEwAglb6Jk2azGfV6PRqNRtmjAAA9VCmKoih7iE602+2o1WrRarWiWq2WPQ4AsAKdrN99s3MCAKwP4gQASEWcAACpiBMAIBVxAgCkIk4AgFTECQCQijgBAFIRJwBAKuIEAEhFnAAAqYgTACAVcQIApCJOAIBUxAkAkIo4AQBSEScAQCriBABIRZwAAKmIEwAgFXECAKQiTgCAVMQJAJCKOAEAUhEnAEAq4gQASEWcAACpiBMAIBVxAgCkIk4AgFTECQCQijgBAFIRJwBAKuIEAEhFnAAAqYgTACAVcQIApCJOAIBUxAkAkIo4AQBSEScAQCriBABIRZwAAKmIEwAgFXECAKQiTgCAVMQJAJCKOAEAUhEnAEAq4gQASEWcAACpiBMAIBVxAgCkIk4AgFTECQCQijgBAFIRJwBAKuIEAEhFnAAAqYgTACAVcQIApCJOAIBUxAkAkIo4AQBSEScAQCriBABIRZwAAKmIEwAgFXECAKQiTgCAVMQJAJCKOAEAUhEnAEAq4gQASEWcAACpiBMAIBVxAgCkIk4AgFTECQCQijgBAFIRJwBAKuIEAEhFnAAAqYgTACAVcQIApCJOAIBUxAkAkIo4AQBSEScAQCriBABIRZwAAKmIEwAgFXECAKQiTgCAVEqJk6+++ireeuutePPNN+Pzzz8vYwQAIKlNq33Cn3/+OWZmZuL69etRrVZj9+7dceTIkRgaGlrtUQCAhFZ95+Trr7+Ot99+O0ZHR2NwcDDefffduHbt2mqPAQAk1XGc3LhxIw4dOhQjIyNRqVTiypUrL93nwoULsWPHjtiyZUtMTEzEzZs3l26bn5+P0dHRpetjY2Px8OHD3zc9ALDmdBwni4uLsWvXrjh//vwrb798+XKcOnUqTp8+HXfu3IkDBw7E9PR0zM3NRUREURQv/ZtKpfKr5/vpp5+i3W4vuwAAa1fHcTI9PR1nzpyJI0eOvPL2s2fPxvHjx+PEiROxc+fOOHfuXIyPj8fFixcjImJ0dHTZTsn3338fb7zxxq+e79NPP41arbZ0GR8f73RkAKCPdPU1J8+ePYvbt2/H1NTUsuNTU1Nx69atiIjYu3dvfPfdd/Hw4cNYWFiIq1evxsGDB3/1MT/66KNotVpLlwcPHnRzZAAgma6+W+fJkyfx/PnzGB4eXnZ8eHg4Hj169J8TbtoUn332WUxOTsaLFy/iww8/jNdff/1XH3NgYCAGBga6OSYAkFhP3kr8368hKYpi2bHDhw/H4cOHe3FqAKDPdfVpna1bt8bGjRuXdkl+8fjx45d2UwAAXqWrcbJ58+aYmJiI2dnZZcdnZ2dj//793TwVALBGdfy0ztOnT+PevXtL1+/fvx93796NoaGh2LZtW8zMzMTRo0djz549sW/fvrh06VLMzc3FyZMnuzo4ALA2dRwn33zzTUxOTi5dn5mZiYiIY8eOxRdffBEffPBB/Pjjj/HJJ5/EDz/8EO+8805cvXo1tm/f3r2pAYA1q1K86lPREmu321Gr1aLVakW1Wi17HABgBTpZv0v5VmIAgF8jTgCAVMQJAJCKOAEAUhEnAEAq4gQASKVv4qTZbEa9Xo9Go1H2KABAD/Xd55y0Wq3405/+FA8ePPA5JwDQJ9rtdoyPj8e///3vqNVqv3nfnnwrcS8tLCxERMT4+HjJkwAAnVpYWPifcdJ3OycvXryI+fn5GBwcjL1798a//vWvrj32L1VnV4ZXaTQaXf17W2/W+u+v336+TPOWOctqnbuX5+nmY/dyHSyKIhYWFmJkZCQ2bPjtV5X03c7Jhg0bYmxsLCIiNm7c2JOIqFar4oSX9Orvbb1Y67+/fvv5Ms1b5iyrde5enqcXj92rdfB/7Zj8om9eEPsqf/7zn8segXXE39sfs9Z/f/3282Wat8xZVuvcvTxPpv/Lbum7p3V6yZcKArCeZVkH+3rnpNsGBgbi448/joGBgbJHAYBVl2UdtHMCAKRi5wQASEWcAACpiBMAIBVxAgCkIk4AgFTEyQp99dVX8dZbb8Wbb74Zn3/+ednjAMCqe++99+K1116L999/v6fn8VbiFfj555+jXq/H9evXo1qtxu7du+Of//xnDA0NlT0aAKya69evx9OnT+Ovf/1r/P3vf+/ZeeycrMDXX38db7/9doyOjsbg4GC8++67ce3atbLHAoBVNTk5GYODgz0/z7qIkxs3bsShQ4diZGQkKpVKXLly5aX7XLhwIXbs2BFbtmyJiYmJuHnz5tJt8/PzMTo6unR9bGwsHj58uBqjA0BX/NG1cDWtizhZXFyMXbt2xfnz5195++XLl+PUqVNx+vTpuHPnThw4cCCmp6djbm4uIv7zNc//rVKp9HRmAOimP7oWrqZ1ESfT09Nx5syZOHLkyCtvP3v2bBw/fjxOnDgRO3fujHPnzsX4+HhcvHgxIiJGR0eX7ZR8//338cYbb6zK7ADQDX90LVxN6yJOfsuzZ8/i9u3bMTU1tez41NRU3Lp1KyIi9u7dG9999108fPgwFhYW4urVq3Hw4MEyxgWArlvJWriaNq36GZN58uRJPH/+PIaHh5cdHx4ejkePHkVExKZNm+Kzzz6LycnJePHiRXz44Yfx+uuvlzEuAHTdStbCiIiDBw/Gt99+G4uLizE2NhZffvllNBqNrs+z7uPkF//9GpKiKJYdO3z4cBw+fHi1xwKAVfO/1sLVeqfqun9aZ+vWrbFx48ZlZRgR8fjx45cKEgDWomxr4bqPk82bN8fExETMzs4uOz47Oxv79+8vaSoAWD3Z1sJ18bTO06dP4969e0vX79+/H3fv3o2hoaHYtm1bzMzMxNGjR2PPnj2xb9++uHTpUszNzcXJkydLnBoAuqev1sJiHbh+/XoRES9djh07tnSfZrNZbN++vdi8eXOxe/fu4h//+Ed5AwNAl/XTWui7dQCAVNb9a04AgFzECQCQijgBAFIRJwBAKuIEAEhFnAAAqYgTACAVcQIApCJOAIBUxAkAkIo4AQBSEScAQCr/B16cGhwSf15gAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "plt.figure()\n",
    "plt.yscale('log'),plt.xscale('log')\n",
    "for L, logical_errors in zip(Ls, log_errors_all_L):\n",
    "    std_err = (logical_errors*(1-logical_errors)/num_trials)**0.5\n",
    "    plt.errorbar(ps, logical_errors, yerr=std_err, label=\"L={}\".format(L)+\",n={}\".format(2*L*L))\n",
    "plt.plot(ps, ps, '--',label='y=x')  #add y=x by weilei\n",
    "plt.plot(ps, 40000000*ps**8, 'r--',label='y=x^8')\n",
    "plt.plot(ps, 400000000*ps**8, 'b--',label='spc y=x^8')\n",
    "plt.xlabel(\"Physical error rate\")\n",
    "plt.ylabel(\"Logical error rate\")\n",
    "plt.legend(loc=0);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "#save simulation result into json file\n",
    "#dictionary={L:{'p_qubit':ps.tolist(),'p_block':logical_errors.tolist()} for L, logical_errors in zip(Ls, log_errors_all_L)}\n",
    "dictionary={L:{'p_qubit':ps,'p_block':logical_errors} for L, logical_errors in zip(Ls, log_errors_all_L)}\n",
    "\n",
    "import json\n",
    "# Serializing json\n",
    "json_object = json.dumps(dictionary, indent=4)\n",
    "#print(json_object)\n",
    "\n",
    "filename=\"result/toric2.json\"\n",
    "#filename='tmp.json'\n",
    "with open(filename, \"w\") as outfile:\n",
    "    outfile.write(json_object)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "range(9, 21, 2)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "log_errors_all_L=[[0.22706666666666667, 0.05316866666666667, 0.009749666666666667, 0.00131, 0.00012, 2e-05, 0.0, 0.0], [0.22396733333333332, 0.04183933333333333, 0.00486, 0.00027, 2e-05, 0.0, 0.0, 0.0], [0.222357, 0.0316, 0.00239, 8e-05, 0.0, 0.0, 0.0, 0.0], [0.22052866666666668, 0.023869666666666668, 0.00122, 0.0001, 2e-05, 0.0, 0.0, 0.0], [0.22166766666666668, 0.01793, 0.00059, 0.0, 0.0, 0.0, 0.0, 0.0], [0.21693733333333334, 0.013439333333333333, 0.00037, 1e-05, 0.0, 0.0, 0.0, 0.0]]\n",
    "Ls"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkgAAAG8CAYAAAAl9Vo9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAACFv0lEQVR4nOzdeXxU1f3/8dcsmck6gYTsJGGHQEiAAIrsW0KUTalaRQRF/VLtz4Xa1t1qK1hrlaoFpbYiaCu1ipWITKKsQlUIhC3sBANZCSSZyT6Zub8/JpnskIQkk8Dn+XjMo5lz79z5zGDJm3POPUelKIqCEEIIIYRwUDu7ACGEEEKIzkYCkhBCCCFEPRKQhBBCCCHqkYAkhBBCCFGPBCQhhBBCiHokIAkhhBBC1CMBSQghhBCiHglIQgghhBD1SEASQgghhKhHApIQQgghRD0SkIQQQggh6rluA9K5c+eYNGkSgwcPJioqik8//dTZJQkhhBCik1Bdr5vVZmVlkZOTw7Bhw8jNzWXEiBEcP34cDw8PZ5cmhBBCCCfTOrsAZwkKCiIoKAgAf39/fHx8uHTpkgQkIYQQQnTdIbYdO3Ywa9YsgoODUalUfPHFFw3OWblyJb1798bV1ZWYmBh27tzZ6LX27t2LzWYjNDS0nasWQgghRFfQZXuQiouLiY6O5r777mPevHkNjq9fv57HH3+clStXMnbsWN577z3i4+NJTU0lLCzMcd7Fixe59957ef/991v0/jabjczMTLy8vFCpVFf9eYQQQgjR/hRFwWw2ExwcjFp9mX4i5RoAKBs2bKjTNnr0aGXJkiV12gYNGqQ89dRTjudlZWXK+PHjlbVr117xPcrKypTCwkLHIzU1VQHkIQ95yEMe8pBHF3ycO3fusr/3u2wP0uVUVFSQnJzMU089Vac9NjaW3bt3A6AoCosWLWLKlCksWLDgitdcvnw5L730UoP2c+fOYTAY2qZwIYQQQrQrk8lEaGgoXl5elz3vmgxIeXl5WK1WAgIC6rQHBASQnZ0NwK5du1i/fj1RUVGO+Uvr1q1j6NChjV7z6aefZunSpY7n1V+wwWCQgCSEEEJ0MVeaHnNNBqRq9T+8oiiOtnHjxmGz2Zp9Lb1ej16vb9P6hBBCCNE5ddm72C6nR48eaDQaR29Rtdzc3Aa9SkIIIYQQ9V2TAUmn0xETE0NSUlKd9qSkJG666SYnVSWEEEKIrqLLDrEVFRVx6tQpx/O0tDRSUlLw8fEhLCyMpUuXsmDBAkaOHMmYMWNYvXo16enpLFmyxIlVCyGEEKIr6LIBae/evUyePNnxvHoC9cKFC1mzZg133nknFy9e5OWXXyYrK4vIyEg2bdpEeHi4s0oWQgghRBdx3e7FdrVMJhPe3t4UFhbKXWxCCCFEF9Hc39/X5BwkIYQQQoirIQFJCCGEEKIeCUhCCCGEEPVIQBJCCCGEqEcCkhBCCCFEPRKQhBBCCCHqkYAkhBBCCFGPBCQhhBBCOI2iKJSdOMGFt9/BajY7uxyHLruSthBCCCG6JkVRKD92DJPRiNmYSEVaGgC68DC8Z892cnV2EpCEEEII0SEsubnkr/sIU6IRy0/pjnaVTofHuHFo/QOcWF1dEpCEEEII0S4URcFWWIimW7eqBrj4t78BoNLr8ZwwHq+4GXhOmojG09N5hTZCApIQQggh2oxis1GacgCz0YgpKRF97z6E/f19AFwC/PH9xRJcBwzAc8IE1B4eTq62aRKQhBBCCHFVFKuV0v37MRkTMScmUpmT4zhmKyrGVlaG2tUVAP/HHnNWmS0iAUkIIYQQVyXjV09i3rzZ8Vzt4YHnlCkY4mLxGDfOEY66EglIQgghhGgWpbKSkj17MBmN9Hj4YVz8/QHwuPFGinfvxmvKFLziYvEYOxa1Tufkaq+OBCQhhBBCNEmxWCj+/gfMiUbMSd9gLSgAQN+/Pz7z5wPgPXcO3W67FVUXD0W1SUASQgghRAOW7GwuvPU25m+/xVZY6GjXdO+O17RpuA0d6mjrikNoVyIBSQghhBDYysupvHABXc+eAKjd3CjcuBEsFjS+vnhNn4YhLg73UaNQaa/9+HDtf0IhhBBCNMpWVkbRzp2YjYkUbd2Krl9feq9fD4DG25uAp36Lvl9/3EfGoNJonFxtx5KAJIQQQlxHbKWlFG3fYZ9TtG07SkmJ41hlTi7WomI0nvb1iarnGF2PJCAJIYQQ15HMZ57B/HXNLfna4CAMsXF4xcXiFh2NSi372IMEJCGEEOKaZDWbKdq2DZPRSMCTT6Lr1QsAr6nTKDt8BENcLF5xcbhGRqJSqZxbbCckAUkIIYS4RlhNJsxbtmDebKR41y4UiwUAt6FR9Pi/hwAwxM/AcMvNEoquQAKSEEII0cVZsrLIevFFiv/3PVSFIgBd7954zYjDa9pUR9v1Ntm6tSQgCSGEEF1MZX4+loxM3CKHAPa1iUr2JoPFgr5/P7xi4zDMiEPXr5/0FLWSBCQhhBCiC6jMy8P8zTeYjEZKftyDrlcv+n6VANgXagxevhx9v77o+/Z1cqXXBglIQgghRCdlyc3FnJSE2ZhIyd69YLM5jqn0OqwmExqDAQBDXKyzyrwmSUASQgghOqkLb66gcMMGx3PXoUPtd5/FxqILC3NiZdc+CUhCCCGEk1kyMjAlJmE2Ggl47jnH3CLDjDjKz5zGEDfDHop6hji50uuHBCQhhBDCCSrOncOcmIhps5GyQ4cc7Waj0RGQPCdOxHPiRGeVeF2TgCSEEEJ0IEtWFucf+SVlqak1jSoV7iNH4hUXh9f06c4rTjhIQBJCCCHaUfmZNCyZmXiOGwuA1s8PS1YWqNW4jx5tn1M0bRpaPz8nVypqk4AkhBBCtLHyU6cwbTZiNhopP3kSbWAg/bZ8i0qtRqXV0vOtv6Dr2xetj4+zSxVNkIAkhBBCtIHyU6cwbdqEyZhIxenTNQdcXNAP6I+1sBBt9+4AuI8a5aQqRXNJQBJCCCFaQVEUAMdK1fnr/03+unX2NhcXPMaNwysuFq/Jk9F4ezutTtE6EpCEEEKIZlIUhbLDhzEbjZiMiQT9/vd43HgDAIb4eCxZmRji4vCcNAmNl5eTqxVXQwKSEEIIcRmKzUbZwYOYjImYjUYsmZmOY+bEREdAch8xHPcR7zirTNHGJCAJIYQQTbDk5HD2zp9TmZ3taFO5ueE5aaK9p2j8eCdWJ9qTBCQhhBACUKxWSvftw5KdjfesWQBo/f1RqdWo3d3xnDwZr7hYPMePR+3m5uRqRXu7rgNSQkICv/rVr7DZbPz2t7/lgQcecHZJQgghOpBSWUnJ3r2YjEbMSd9gzctD7e2NYcYMVC4uqFQqQle/h0toKGq93tnlig503QakyspKli5dytatWzEYDIwYMYLbbrsNH1mTQgghrnkl+/ZTuGED5m++wZqf72hXGwx4TZmCtajIcUu+vl8/Z5UpnOi6DUg//vgjQ4YMISTEvvHfzTffjNFo5K677nJyZUIIIdqaUlEBKhUqFxcAinbuoODTTwHQdOuG1/RpeMXG4XHDaFQ6nTNLFZ2E2tkFtNaOHTuYNWsWwcHBqFQqvvjiiwbnrFy5kt69e+Pq6kpMTAw7d+50HMvMzHSEI4CePXuSkZHREaULIYToALaKCsxbt5L51NOcGDeeoh07HMcM8fF0u/NOwv7xd/rv3EHQ73+P5/hxEo6EQ5ftQSouLiY6Opr77ruPefPmNTi+fv16Hn/8cVauXMnYsWN57733iI+PJzU1lbCwMMcCX7VVL/YlhBCia7KVlVG8axcmo5GiLVuxFRU5jhXt2InX1KkAuA4YQNBLv3NSlaIr6LIBKT4+nvj4+CaPv/HGGyxevNgx8XrFihUYjUZWrVrF8uXLCQkJqdNjdP78eW644YYmr1deXk55ebnjuclkaoNPIYQQoq1U5uVxOjYOW0mJo03r749XbCyGuFjcRoxwYnWiq+myAelyKioqSE5O5qmnnqrTHhsby+7duwEYPXo0hw8fJiMjA4PBwKZNm3jhhReavOby5ct56aWX2rVuIYQQzWMrKaFoxw4sWdn43rcIAG2PHrj07InVZMIQF4tXXBxuw4ahUnfZ2STCia7JgJSXl4fVaiUgIKBOe0BAANlVi31ptVr+/Oc/M3nyZGw2G7/5zW/w9fVt8ppPP/00S5cudTw3mUyEhoa2zwcQQgjRgLWoiKJt2zEbjRTt3IlSVobK1ZXud96B2t0dgND3/4a2Rw8JReKqXZMBqVr9OUWKotRpmz17NrNnz27WtfR6PXpZA0MIITpc0c6d5P/rE4q/+85+N1oVl7AwDHFx9raqgOTi7++sMkVrWCuh3ARlBVBW6HjYSgsoN+dDv8m49Yx2SmnXZEDq0aMHGo3G0VtULTc3t0GvkhBCiM7FWliISqdzrFZdfuIERVu2AKDr1QuvGXEY4uLQDxokN9c4WxMBh7JClNICKooLqCjKx1qSj63U3q4uL0RrMeNiMaO3lTR6WTXgBuzNszDyZxKQ2oxOpyMmJoakpCRuvfVWR3tSUhJz5sxxYmVCCCEaU5mfT9G332IyJlL8v/8R/Mof8K76+9orbga24hK84uLQD+gvoagtOQJOYaMPpayAypICLMX5WEsKUEoLUZUXoi434WIxo7MWN3lpFaCvelxJkeKKCXdMikfV/7pjwgM3bXBbfdIW67IBqaioiFOnTjmep6WlkZKSgo+PD2FhYSxdupQFCxYwcuRIxowZw+rVq0lPT2fJkiVOrFoIIUS1yosXMX/zLWbjZop/+BGsVsex0oOHHAFJ1zMEv0f/n7PK7Nxs1ibDTe2HtbSAyuJ8rKWFUFaAqtyEtsKES2XTAQfsIcel6nE5xYoeEx5VwaZu0ClSe2LRGrDqvLDpvVG5eaN264bWoxs6j+64eXbD4OmOt5sL3m4u+Ljp6F31s07rvLlkXTYg7d27l8mTJzueV0+gXrhwIWvWrOHOO+/k4sWLvPzyy2RlZREZGcmmTZsIDw93VslCCCGqWAsKODlxElRWOtr0EREY4uLwio1F36e384rrDErzIXM/ZKaAOdsxL6eyTi9OIdorBJxqmqpHUy4XcEx4UKzywOLihVVnQHH1Ru3qjca9Gy4e3dF7dMerVsDp5u5CaNXP3m4uuLpc7p07L5XS2IqJ4opMJhPe3t4UFhZiMBicXY4QQnRalpxczElJWLIyCfj1rx3tZ+ffg1JejldcLIbYWHTX6z9gK0og+yBk7KPspz0oGftwM59t0SWuFHBMijtmPBy9OLh6o3Lrhta9GzrPbhjcawKOoSrkeNcKOe46zTUztNnc399dtgdJCCFE52XJysKcmIjJmEjpvn32Ro0G3wcecGwCG/b391G7ujqxSiewWiA3lYr0vZjP/Igmcx8G82nU2IcXa38bP9n8Oaj04SclgELFAxMemHHH4uKForcPU6nduuHi0Q1Pd7cGoSbAzYX+1c/dXfDUaVGrr42Q0xEkIAkhhGgzJmMiF//xd8oOHKzT7jZiBIa4WFTaml8713w4stmw5Z3i4snvKUn7AV3OAXyLjqNTKtABtVfey1W6ccDWl4NKH7I9BmMLGkbPkJ4MCvRiQnc3+9CVmw5PVy0aCTkdQgKSEEKIVqtIT0fj7Y3G2xsAa36+PRypVLjHxOAVG4tXXCwu1/oSK4pCflYaOcf+R8VPe3DPO0BgyXE8lWL86p1qUtw5YOvDSW1/CroPhZAR9Azty8AgA78I8MRdJ7+aOwP5UxBCCNEi5WlpmI2JmBKNlKceJeD55/CZPx8Ar9jpKDYrXtOmXbOLNpZZrKSlnyPvxP+wnduL4dIhwsqO40sB3eufq7iQSm/OuQ2iyDcKl7CRBPUezMAgb8Z56q+ZeT3XIglIQgghrqj89GlMRiPmzUbKT5yoOaDRYMnIdDzV+vjgc/fdTqiw7dlsCufzSzlxPouC03tRZ+6je8Eh+lpOEKG60OD8SkXNGXU42Z6DKfePxjV8FMEDhhPt340RMizW5UhAEkIIcVlWs5kzc28Fi8XeoNXiceONGGbE4Tl1qmPSdVeWX1zBsWwzJzMvYvopBZfs/fibjzBYOc1kVQYaVa0bvquyTqYmhDzDECyBw/HqM5qQiNEM8DQwwDkfQbQxCUhCCCEA+36V5ceOYTIasZw7T8ifXwdA4+WF58QJYKnEKy4OrymT0XTr5txiW6m80sqp3CKOZ5s5nlVA4fmjuOfup1f5CaLUp7lTlY5eVbU2kwpHGMrX+lHQPRIleATd+91At36jCXbrjvPWeRbtTdZBaiVZB0kIcS1QFIWyI6mYjZsxGROxpKc7jvXb8i0uwcGO87rSfBlFsQ+PHc82cyzbxLEsEwVZp+lecIhIThOtPkOkKg1PVVmD15ZqDZh9otCGxuDddzSa0JHgFeiETyHag6yDJIQQ4rIK//tfLrz1NpaMDEebSq/Hc8J4vOJm1Okl6szhqLDEwrFsE8dzzBzLNnMsy8TFnPP0tZwgWn2GKNVp7lKfwVdlbvBbz6Jxo6zHUPRhMejCR0HwCNy698KtE39e0TEkIAkhxHVAsdkoTTmAS0gILgFVd5dptFgyMlC5ueE5cSKGuFg8J0xA7eHh3GKbUFFp4/QF+/DY0WyTfZgs20xR4SUi1WlEq04zTn2Gh9Vn6KnKA13d19vULlT2iMAldCSqnjEQPAIXv4G4qLvmVhiifckQWyvJEJsQorNTrFZK9+/HZEzEnJhIZU4Ofk88QY//ewgAa1Exxbt24Tl+HGp3dydXW0NRFDILyziWZeJYVQg6nm3m9IUiNLZyBqt+Ikp9hij1aYapTtNXndXwGqigxwBUISMgxB6GCBgCLtf44pTiimSITQghrkOK1UrJ3mTMRiOmpESsF/Icx9SenigVFY7nGk8PDHGxzijTwVRmqZonZOZ4toljWWaO55gxl1WiwUp/VQZR6tPcqzpDtPY0A9XncKnalqMO7zAIGe4IQ6qgaHCVf7yK1pOAJIQQ1xDFYuH8L36BraQEALXBgNeUKXjFxeIxdixqne4KV2gfFquNMxeK7ROma/UKZRSUVldOuCqHaNUZZqhPM0x3hkj1WVwpb3gx9x72IBQywt4zFDwcPOuvVy06qwprBXmleZRVllFaWUpJZYnj5zJrGRE+EfTv3t/ZZUpAEkKIrkixWCj+/gdMxs1UnD1Lr48+Auz7m3nPnYOtvBxDXBweN96Iygmh6IK5nG+O5vDDmYscqxoes1hrZnT4k0+0+jQ/157hBpc0BqvO4GkzN7yQzguCh9WEoZAR4B0KMom6zVltVhQUtGp7NCgsL+Ss6WxNeKn63+rHuJBxDPQZCEDqxVTeP/R+nfNqh57HRzzOrf1vBWBf7j4eTHywyTqWxiyVgCSEEKL5bBUVFO/ejdmYiHnLFmyFhY5j5SdPou9v/6US+MILTqkv/WIJxiPZJKZms/enfKpnuBoo4gZ1GqP0aYxx/YmB1pN4V9YM/aFUPTQ6CIyqG4Z8+4Na7YyP06koikKlrZKSyhJKK0vppu+Gq9Y+nyqrKIvUS6mNhpPSylJu7X8rA7rbl6/cnbmbdw+86zindu+NxWbhj+P/yM19bgbgx+wfWbptaZM1eeu9HQGpoLyApJ+SmjzXVGFy/OyudUev0eOqdcVN64arxv6/1Y8gj6Cr/r7aggQkIYToAvI//ZTc1/6EzVzTy6Lx9cVr+jQMM2ag6927w2tSFIXULBPGIzkkHsnmWHZNbYNU6Tzc7X+MZx/dy87VvKh6xEylBr8I+7yh6jDkPwS0zhkCvFqKolBuLa8bTqyl9Db0xt3FPgH++KXjHLhwoE4vTO1A83/R/0ffbn0BSDiTwMqUlXWCjlWpmXu1cupKxvccD8D3Wd/zwu6mQ/Ew/2GOgGSuMLM/d3+T55ZZa9aF8tZ5E+IZUhNiXGrCjKvWlVCvUMe5fb378swNzzjOc9O4Oc5z07rh716zL1+UXxR779nbkq/XKSQgCSFEJ2MrK6No505c+/dH16sXAC7+/tjMZrR+fnjFxuIVF4t7TAwqTcfeom61Kew9e4nE1ByMR7I5n1/qOOatLuVR/xTm2LbQw3QEaq/B2L133Z6hoGjQddxyAlablTJrWZPhZFTgKEeQ2ZO9hx+zf2wwpFT9/OWbXibUYA8Haw6vYeUBe5BRaHhT+Lr4dQzzHwbYg8zre19vssa5/eY6AlJZZRnnzOcaPU+r0mKxWRzP/d39ifaLdoQRN42bPaRUBZtehl6Oc4f5DePNSW/W9N7Ufo3WDY9afyajg0azed7mZn2/AR4B3DXormad21VIQBJCiE7AVlJC0Y6dmBONmLdtRykpwffBB/D/1a8A8BgzhvCP1uE2YgSqDh5yKrNY2XUqj8QjOXxzNIeLxTV3wrm6qHigZxa3a7YRlp2EqnrStdoFBsZD9F0QdiO4+zTrvQrKCigoL2gyzMzqOwu9Rg/A5rTN/JD9Q5NDS3+P+7uj5+LPyX9mXeq6Jt/3v3P+S59ufQD70NK7B95tusbyAkKp6T0prSytc1yn1jl6W2rr7d2bqWFTa0JJVYBxd3HHTetGuHe449xJoZPo161f3WGoqtDjonapc92xIWMZGzL2cl+rQ4BHAAEeAc0693onAUkIIZxEsVgwJyVhMiZStGMHSmnNL1qX4GA03Wo2gVXpdLiPHNlhtZnKLGw9lkvikRy2Hc+luKJmeMfbzYUZvUq5VfMtUTmJuGfVrENUYAgkNWwkN854E7WnPZx8lPpRk0GmrLKMpNuT8HCx91y8kfwGG05taLKuiT0n4uduv2Mt5UIK/znxnybPLbYUO36uDisqVA0CipvWrc5K4UN7DOXOgXfW6WGpPU+mp1dPx7m39r+VqeFTHcf0Gr1jknN9E3pOYELPCU3WW1sPtx70cOvRrHNF+5CAJIQQHUiprESlrfqrV6Ui+/d/wJqfD4BLz54YZsThFReHa2Rks7f3qD18ZLFaCPKsmeS6P3c/OSU5lFpKHedUhxObYuPXo37tOPflnSvYeuwSeXnBmE3BoNQM36m0BSwcOYQ7DUcYmPUFnExCXTWkVKxS8bWnOxs8PTmod4Hig/zo6oVb1WuPXjrKtnPbmqy/tLLUEZA8XDzwdPGsG05qDQPV/k7Gh4ynm75bnQm+tQNNoEfN/mkPRT3EQ1EPodfor/i9tiTIeOu98dZ7N+tc0bXIStqtJCtpCyGaq6zgIgXfJFFkTKQy7SzKJ++ASkWEbwQXVq5EKSvnyFADGcE6ymzllFhK6oQZV60rL4550XG9pduWcijvkON4ubVmrSCDzsCuu3Y5nj9gfIAfsn9otC6tWsvn8d/Z7zw7ks2+9Hwc29cDal0OWq9UBrjt4WcVp1lo1aMqrrn77IiHgW09wtjvF45G71VnIu9vRv3GMafnh6wfOGc+V6fnpvYjxDOkyV4XIdqarKQthBDNVG4td/Sw1Fm0rrIMnUbHqMBRjnM/PvoxeaV5dW6Trv7Z392fV8a9AoC1oIA/vXE74cmZRKbZ0Npq3u+p9++gvH9PNs/bjN/DDwPwy413cHTv0Ubr667vXicg5Zflk12c3azPNtBnIFbFWqtnxY2SYh/OZ/fg7HkDk1/fVutsFX0C1NzQ15WJfbQMy0+le+qP6LMOVB03g4c/DLsLhi9gSI/+DGlGDTcE3cANQTc0q14hOgsJSEKILuO8+Xydxedqr8Lr4+pTZ1jktT2vYSo31Tm3+rX9uvXjTxP/5Dg39j+xXCq71Oh7DvIZxKezPnU8//jox03eXRTmFQbApX/+k5xly5lTWek4lt4DfozQkDLEHVOQO8GuvnVee1PwTYQbwhsMLblr3fF08axz7jM3PEOFtaLBua4a1wbDR78e9WsqrTZ+PHuJxCM5JKXm1Fq9GrRqFWP6+hI7OIDpEQEEmg7C/rXw3w1QPYdHpYH+sTBigf1/NXUnCQtxLZKAJIS4KpW2SqyK1XFnkdVm5VDeIXsosZTWCSillaWEeoUS1yvOce6T25+se7eSte6t129MesPxXnO+mEOFraLROkYFjqoTkDae3khBeUGj5+o0ddfacdPaZ8u4qF0aDAP1NtRdX2hWn1mYKkyOUOJVZMVvTxq2iD64D4kEwHXAAKisRDOgL7qpE/GMnU6/AUOIVWubnP/yeMzjTXzDDTVnleEyi5WdJ/MwHsnm26M55JfU3Bbu5qJh4gA/4iIDmDIwAG9bPhz4BD5aB3knai7i09ceiqLvAq/ARt5FiGuXBCQhrmGKolBhq3D0nOg0Onxc7bdbl1WWsStjl70HpirM1A4nkT0iuaXPLYB9cbmHv3m4JuhY7IvglVaWUmmrZFafWSwbvwyASqWSBV8vaLKmyaGTHQFJrVKz9dzWOgvg1VZ79V2Abq7dsFgtjU7gHdh9YJ1zHxj6AJW2SkcvTO3zu+u71zn389mfo9PomjUP5hfDfoElJxfzN0mYjYmU7N0LNhvd7vo5QbFLAHAbPpy+m792rGHUUQpLLGw5noPxcA7bT1yg1FLzvXZ3d2FqRABxQwIZ378HrmoFTn0DX/4WTmwGW1Vvl4s7DJ5rD0ZhY2RLD3HdkoAkhBPZFFvdIaOqkOLr6uu4E6mwvJCv075udC+kMmsZE3pOYHbf2QBkFmVyv/H+Ote0KTWTX+4adBfP3PAMAEWWIh7f9niTtc3pO8cRkDQqDSkXUpo8t/bquzq1jjCvMPRafYPVdF21rgz2Hew4V6VS8fyNz6NVa+ucVx1oDLq6Eyi/vf3b5n2xwMIhC5t9bvVk4stRLBby//UvTMZESvftg1r3t7hGReE6KMLxXKXRdFg4yi4sIyk1G+ORHL4/c5FKW01dId3cmD7YHopG9eqOVqOGi6dh+7tw4F9grrk9n5CR9lA05DZwlRtPhJCAJMRl2BQbJZaSOhNya/ey9PHuQ5jBPu8kqyiL/5z8T8MF66z2Hpd5A+Y5gszRi0e59+t76wSL2h6Keoj/N/z/AXCp7BKv/PBKkzX6ufk5rqtRacgoymj0PK1aS+2bVt217gzzG9Zgvkv18yG+NdNvXbWurJi0ouFt11UPd21NwFCpVHx121fN+XoBmDdgXrPP7WhWkwlN9V0uWi2XPlyLJcP+/boNG4ZXXByG2Om4hIR0aF2nLxRhPGIPRQfOFdQ5NjDAi9gh9lA0JNhgH9KrKIHD/4Z96+Cn72pOdveFqJ/D8HsgYDBCiBoSkMQ1ocRSQk5JTuPhpLKUmIAY+njbV8k9fuk4/zr2rwa7TdfeD2lmn5mA/fbkh5IeavJ9fxXzKxZFLgLgQukFVh9c3eS5te/i0Wl0DcKRq8a1Tg9KNW+9N9PDp9fZA6n2/0b41PRc+Lj5sC5+XYN1YVy1rg1W33V3cWfdzU2vLFybWqVmavjUZp3b1VWcO4c5MRHTZiOW9HT6f7cTlYsLKpUK3wcfQKmowCs2FpfAjpuToygKB88XVoWibE5fqFkAUaWC4aHdiBsSSOyQQHr38Kh+EWTus4eiw59BefVwpQr6TYXhC2DgzV127zMh2psEJNFuqoePXDQujl/OeaV5pBWmNb41gLWUuF5xjiCzJ3sPHx75sNFdp8usZbx000uOuSy7Mndddtfp3435neO6eaV5fHbysybPvVRaczdT9eRdtUpdZzXd6t6W7q41c1n83f35+cCfN7pRo6vWlf7daibWhnmFsXne5jrXU6sa3z7Cx9WnzkTly3FRuzj2fRLNV/HTT5g2GzEbjZSlptYcUKkoO3IEt2HDAOj+8593WE0Wq40f0y5VrVGUQ7apJlC7aFSM6duDuCH2O8/8DbW2tCi5BAfX24NR7pGa9m5h9lA07G7w7okQ4vIkIF3H6m/e6O/u77gT6WzhWY5dOlbnNuraYebeIfcSbrDvG2Q8a+Qfh//RcAuBqh6Sd6e969gnaMf5Hby4+8XGCwL6ePdxBJmLpRfZfn57k+eWWEocP3toPfBy8Wq0h8VV61pn76Fe3r345bBfNhgmqn5efas22LccSL4nGRe1yxVX3w30COTZG5+97DnVXDQuhHh27LCMaNyldR+R80qtIUy1GvcbRmOIi8Nr6lS0fn4dVktJRSU7TuSReCSbb4/lUlhac+eZh07DpIH+xA4JYPIgfwyutXoEbVY4sw32r4NjX4G16k4/jR4Gz7YHo17joYP3cBOiK5OA1EkpioLFZmm467S1jAifCFy19n8xHrhwgIMXDjY6ebfUUspvRv/G8Yv446Mf88HhDxzXq3+79Ec3f0S0XzQA289vv+yu01PDpzoCkqnCROrF1CbPLaus+Zdvd313env3rrOvUe2AUjs0RPaI5KWbXmowtFT9qL1P0U0hN7H77t3N+m5DPEP4v+j/a9a5GrUGDR27W7poP+UnT2IyJuI+ciQeN9qHPN1HjwKNBo8bb8RrRlUo8mnexqptIb+4gm+P5WI8ks3Okxcos9RMqvf10DEtIoC4yABu6tsDV5d6/y3m/wQpH0PKP6Gw1tpMgVEw4l4Y+jNwq3vHnhCieSQgdUJvJr/Jh0c+bPLW5//O/a+jl2Xn+Z28d/C9Jq/1wNAHHKGj3FpOTklOg3NUqHDTulFpq1nULsQzhNGBoxvdqLF6a4BqY4PH8tepf20wBFU7+FSbHDaZyWGTm/U99PTqWWdTSCFaSlEUyk+cwGw0YjImUnH6NACGWbMcAUk/YAADdn2Hplu3Dqsrs6CUxKpJ1j+evYS11p1nPbu7ETckkLghgcSEd0ejrtdzaSmDYwn23qIz26FqPzRcvWHoHfY70YKiO+yzCHGtkoDUCWlUmjrhSKuqewt07TuRBvoMJL5XvGP/o/rhpPamlTP7zOSGoBscmz66ad1wc3FDp9Y1GD6aFj6NaeHTmlVvsGcwwZ7BV/mphWg7SmUlF956G7PRSMVPPznaVS4ueIwbh9eUmpCuUqnaPRwpisKp3Jo7zw5lFNY5PijQyxGKIoK8Gh/OzT5kn1d0cD2UFdS0955oH0KLmAkubg1fJ4RoFdmstpXac7PawvJCyq3ljqBT/+4jIURdiqJgOXcOXVjN/LEzc+ZSfvw4Kp0OjwnjMcTF4TlpEhovrw6pyWZTSDlfgPFINklHcjiTV/fOs5Hh3e13ng0OJMy3iXWYSgvg0Kew/yPISqlpN4TAsPkwfD5079WeH0OIa45sVtuFeeu9nV2CEJ2eYrNRdvAgJmMiZqORyosXGbB7F2oP+23uPR5+GKXSgufESWg8PTqkpopKG9+fuWgPRak55JrLHcd0GjVj+/kSNySQqREB+HnpG7+IzWZfq2jfOjj6JVTP4VO7wKCbYfi90HcyqGVunBDtSQKSEKLLUGw2SlNSMG3ejDkxicrsmh3tVW5ulJ04gfvw4QAY4mI7pKbi8kq2n7iA8Ug2W47lYi6rmcvnqdcyeZA/cUMCmDjADy/Xy/QGmzLtE673fwT5Z2va/SLs84qi7gSPHk2+XAjRtiQgCSG6jPyP/1nnlny1hweekyfjFReL57hxqN06Zg7OxaJyvj1adefZqTwqKmvuPOvhqa/a3iOAMX190Wsv09NTWWHfB23/Ovu+aNXbwui8YOg8e29RyAjZD00IJ5CAJITodJTKSkr27sW0eTMeN92EIdbeG+Q5eRIX3noLrymT8YqbgcfYm1DrmxiqamPnLpWQmJqD8Ug2e89eotaNZ4T7uldNsg5gWGgjd57Vl3vMHooOfAIleTXtYTfZe4sGzwFdxwwLCiEaJwFJCNEpKBYLxT/8iNloxPzNN1jz8wGozMl1BCRdz54M2L0LlUv737igKArHc8wYD+eQmJrNkUxTneNDgg1V23sEMDCgiTvPais3w+HP7cHo/J6ads8A++rWw+6BHv3a4ZMIIVrjug1I586dY8GCBeTm5qLVann++ee5/fbbnV2WENcdxWYj+8UXMScmYS2suf1d060bntOm4n3zzXXOb89wZLMp7EvPt2/vkZrDTxdrVmtXq2BULx/ihgQyfXAAoT5N3HlWm6LAuR/sE66PbABL1Z1sKg0MmGHvLeo3HTTX7V/FQnRa1+3/K7VaLStWrGDYsGHk5uYyYsQIbr75Zjw8pFtbiPZkq6ig7MgRx2RqlVpN+ekzWAsL0fj44DV9Ooa4WNxHj0albf+/osorrew+fZHEIzkkpeaQV1TrzjOtmgn9exA7OJCpEf74ejZzOM+cAwf+ZZ9wffFkTbtvP/uaRdF3gVdA068XQjjddRuQgoKCCAqyL6Lo7++Pj48Ply5dkoAkRDuwlZVRvGsXJqORoi1bsZWW0v+7nWi727fB8Hv0UQDcR41EpWn/29eLyivZeiyXxNQcth7Lpai85s4zL1ctUwf5EzskkIkD/PDQN/OvSWslnEqy9xad2AzVi726uMOQW+3BKOxGmXAtRBfRaQPSjh07+NOf/kRycjJZWVls2LCBuXPn1jln5cqV/OlPfyIrK4shQ4awYsUKxo8f3+L32rt3LzabjdDQ0DaqXghhKy2laMdOzEYjRdu2YSupGa7S+vtTcfasIyBVb/vRni6Yy/n2qH2S9a5TF6mw1tx55u9VfedZIDf28UWnbcGmrhdP2+cVpfwLimqWHaDnKHsoirwN9B2zOKUQou102oBUXFxMdHQ09913H/PmzWtwfP369Tz++OOsXLmSsWPH8t577xEfH09qaiphVavpxsTEUF5e3uC1iYmJBAfbt8a4ePEi9957L++//377fiAhrjOFX24k+8UXHc+1QUEYYmPxiovDbVg0qg7YWT79YgmJqdn2O89+yqf2vgG9e3gQO8Qeiob17Ib6Snee1VZRDKn/tfcWpdfaJNnd1z58Nvwe8I9ouw8ihOhwXWKrEZVK1aAH6YYbbmDEiBGsWrXK0RYREcHcuXNZvnx5s65bXl7O9OnTefDBB1mwYMEVz60dtkwmE6Ghoe2y1YgQXYm1qIiirdswJxrxGD+e7nfcAUDlxYucvetuvKZNwzAjDtehQ698p9dVUhSF1CwTiUfsPUXHss11jg8N8SauKhT18/dsWT2KAhn7YP9aOPQZVFRdW6WGvlPtE64HxINW14afSAjR1q7prUYqKipITk7mqaeeqtMeGxvL7t27m3hVXYqisGjRIqZMmXLFcASwfPlyXnrppVbVK8S1xmoyUbR1KyZjIsXffYdSUWFvzy9wBCStry99jZvbPRRZbQp7z15yrFF0Pr/UcUyjVjG6lw9xQwKIHRJIcLdWLCRZfNG+Qez+dZCbWtPevZe9pyj6bvAOufoPIoToVLpkQMrLy8NqtRIQUPcukICAALJrbT1wObt27WL9+vVERUXxxRdfALBu3TqGDh3a6PlPP/00S5cudTyv7kES4nqiKAoZjz6Kedt2sFgc7bpevfCaEYchLq7O+e0VjsosVnadyiPxSA7fHM3hYnGF45heq2bCAD/7nmeD/Onu0YoeHZsVTm+1h6JjX4Gt6rNqXSFitr23KHwcdMAwoRDCObpkQKpW/y9fRVGa/RfyuHHjsNlsVz6xil6vR99BK/YK0VlU5udTsncvhunTAfv/52wVFWCxoOvXF0NsHF5xcegH9G/3niJTmcV+59mRHLYdz6W4wuo45u3m4rjzbMKAHrjrWvlXW/5Z2P8xpPwTTOdr2oOi7ROuh94Obt2u6nMIIbqGLhmQevTogUajadBblJub26BXSQjRMpUXL2L+5lvMxs0U//AjWK24ffsNLiH2YST/xx5D9etfo+/X/qs+55rKSDqag/FIDv87nYfFWjNlMtDg6phkPbq3Dy6aVvbmWMrgWALsWwtp22vaXbvZN4gdfg8ERV3dBxFCdDldMiDpdDpiYmJISkri1ltvdbQnJSUxZ84cJ1YmRNdUefEiJqMRszGRkj17oFbvqn5wBJUXLzoCkuvgwe1aS3ZhGV8dymLToSz2pde986yvn0fVnmeBRPX0vrpeq6yD9iG0g/+GsoKa9j6T7L1Fg2aCi2vrry+E6NI6bUAqKiri1KlTjudpaWmkpKTg4+NDWFgYS5cuZcGCBYwcOZIxY8awevVq0tPTWbJkiROrFqLrqD0kXfy/78l5+feOY66RkRhmxOEVG4uuatmM9nTBXM7mw1lsPJDFnp8u1QlFQ0O8iR8aSOxg+51nV6U0Hw79x95blH2wpt3QE4bPh2HzoXv41b2HEOKa0GkD0t69e5k8ebLjefUE6YULF7JmzRruvPNOLl68yMsvv0xWVhaRkZFs2rSJ8HD5y02IplgyMzElJmI2JuI5aRI9/u8hADwnTcJtZAxeU6baQ1HP9r8rK7+4gs1Hskk4mMn/Tl/EVisUxYR3Z2ZUEPGRQQR6X2Uvjs0GZ3fae4uOboTKMnu7RgeDbrEPofWZDOr2X8FbCNF1dIl1kDqj5q6jIISzVZw/j9mYiCnRSNmBml4T/eAI+nz+eYfWUlhqIfFINgkHs9h1Ko/KWqkouqc3M6OCuSUqqHW34zd4s/P21a33r4OCn2ra/YfY70Ibegd4+F79+wghupRreh0kIcSVKYpC+sJFlPz4Y02jSoV7TAxecXF4xU7vkDqKyiv59mgOGw9ksuNEXp0tPgYHGZgZHcTMocGE+bpf/ZtVVsDxTfZQdHoLKFXvpTdA5Dx7MAoeIfuhCSGuSAKSENeI8rQ0inftpvv8u1GpVKhUKrQ9eoBajfuoUXjFxWKYPh2tn1+711JaYWXLsVwSDmay5Vgu5ZU1oai/vyczo4KZGR1EX7+rnFNULfeofduPg59AycWa9vCx9gnXg+eArg0CmBDiuiFDbK0kQ2yiMyg/dcpx91n5iRMA9Nn4Jfr+/QGoOHcOtZubPSi1szKLle0nLpBwMItvj+ZQUmudot49PJgZFcTMqGAGBrbRxq1lJjjyuT0YZeytafcMhGF32+cW+fZtm/cSQlwzZIhNiGtUxfnzFH6+AZPRSMXp0zUHtFo8xoyxL+RYRdfOq71XVNrYdSqPjQczSTqSg7m80nGsZ3c3e09RVBBDgg1ts5CkokD6/+yhKPULsJTY29VaGDDD3lvUbxpo5K82IcTVkb9FhOjkFEVBqahAXbWSe8Xp0+StXGk/6OKC50034TVjBl5TJqPx9m73eiqtNv535iIJB7LYfCSbwtKaLUcCDa7cEhXEzKgghoV2a7vVtc05cOCfsP8juFiz/Ae+/e3ziqLvAk//tnkvIYRAApIQnZKiKJQdPoI50YjJmIghLhb/X/0KAI8xYzDcHI/npEl4Tp6MxquNhqwuw2pT2HP2EgkHM/n6UHadvc96eOq5ZWggM6ODiQnrjlrdRqHIWgknE+0Trk8YQakasnPxgMhbYfi9EDpaJlwLIdqFzEFqJZmDJNqaoiiUHTyIyZiI2WjEkpHhOKbv358+G7/s0HpsNoX95/LZeMC+qnWuudxxrLu7C/FD7T1FN/T2RdNWoQjg0hn7EFrKx1CUU9Pec7S9t2jIraBv/1AohLg2yRwkIbqYs7ffQdnhw47nKjc3PCdOxBAXi+eECR1Sg6IoHDxfSMLBTL46mEVmYZnjmMFVS9wQe0/RTX19W7/3WWMqy+HYV5C8pu5+aAA3/B/E3A/+g9ru/YQQ4gokIAnRwRSrldL9+ynasRO/xx9DpbYHDdfBgyk/cwavSZPwiovDc/w41O7tf2u6oigczTKTcDCThINZpF8qcRzz1GuZPjiAmVFBjO/vh07bhqEIIO+kPRQd+Fet2/NV0HcKxCyEAfGg1bXtewohRDPIEFsryRCbaAmlspKSvcn2OUVJSVgv5AEQ/q9/4j58OACVly6hdndH7doxG6SezDGz8WAWCQczOXOh2NHu5qJhSoQ/s6KCmDTQH1eXNt6Cw1IKqV/Cvg/hp1017V5B9lvzhy+Q/dCEEO2mQ4bYysrKcO2gv8yF6IrKT5/m0tp1mJOSsF665GhXGwx4TZmC2sPD0ab18Wn3es7mFZNwMJONB7I4nmN2tOu0aiYP9GNmVDBTI/xx17VD53JOqj0UHfgEygrsbSo19I+FmEXQb7rcni+E6DRa/LeRzWbjlVde4d133yUnJ4cTJ07Qp08fnn/+eXr16sXixYvbo04hugTFYsFWXIymWzcArAUFFKxfD4DG2xvPaVMxzJiBxw03oNJ1zNDRuUslfHXI3lN0OMPkaHfRqJjQ34+Z0UFMiwjAy9Wl7d+8ohiObIDkD+F8rS1PvENhxL0wbD54t//GuEII0VItDkh/+MMf+PDDD3nttdd48MEHHe1Dhw7lzTfflIAkrju2igqKd+/GbEzEvGULhpvjCXrxRQDchg/HZ+FCPCaMx2P0aFQu7RBCGpFdWOaYU5RyrsDRrlGruKmvL7OigokbEoi3ezvVk3XAHooOfQrlVaGsejHHmPug72RQt/HQnRBCtKEWB6S1a9eyevVqpk6dypIlSxztUVFRHDt2rE2LE6KzspWXU7xrF2ajEfOWrdjMNcNVpfv2O35WqdUEPP1Uh9R0wVzO14ezSDiQxZ6fLlE9u1Clght7+zIzOogZQwLx9dS3TwHlZjj0H/swWmbNd0D33jW9RV4B7fPeQgjRxlockDIyMujXr1+DdpvNhsViaeQVQlx7zt5xJ+XHjzuea/388IqNxSsuFveYmA6r41JxBZsPZ5NwMJPvz1zEVuuWi5Hh3ZkZFcTNQ4PwN7TTXEFFgYx9kPwBHP4cLFWTvdUuEDHLfidarwmgbuO734QQop21OCANGTKEnTt3Eh5e9y6TTz/9lOFVd+MIca2wlZRQtGMnRdu2EfTyS455Qx5jx2ItLMQQF4tXXBxuw4Y5btdvb4WlFoxHskk4mMWuU3lYa6Wi6NBuzKoKRcHd3NqviNIC+/BZ8hrIqVm7Cd/+9lAUfRd4tP8GuUII0V5aHJBefPFFFixYQEZGBjabjc8//5zjx4+zdu1aEhIS2qNGITqUtaiYou3bMBsTKdqxA6XMvlii4ZZb8Bw/DgC///dL/J/8VYeFoqLySr5JzSHhYCY7TuRRYbU5jg0OMjAr2r4pbKhPO66bpChw7gd7KDryBVSW2ts1ehgyF0YshPCbZOsPIcQ1ocUBadasWaxfv55ly5ahUql44YUXGDFiBBs3bmT69OntUaMQHaLs2DEuvP0OxTt3olTU7DXmEhqKIS4WXVioo03t1o69M1VKK6x8eyyHhANZbD2eS3llTSgaEODJzCh7KOrj59m+hZRcsi/kmPwh5NUMK+I/2H57ftQd4Na9fWsQQogO1qpFR+Li4oiLi2vrWoToUNbCQmzFxbgEB9sbVCqKvv0WAJfwMAxxMzDMiEMfEdF2u9JfQZnFyvYTF0g4mMU3qTmUWqyOY316eDAzKoiZ0cEMCGjnvcgUBc7utIeio1+CtSowurhD5G0wYhH0HCm9RUKIa1aLA1KfPn3Ys2cPvr6+ddoLCgoYMWIEZ86cabPihGhrlfn5FG3ZgslopPh/32OYMYOQP70GgH7AAPx/8xs8xt6EfsCADgtFFZU2vjt1gYQDWSSl5mAur3Qc69ndzdFTNCTY0P41FV2wbxK7by1cOl3THhhl7y0a+jNw9W7fGoQQohNocUA6e/YsVqu1QXt5eTkZtXYfF6KzqLx0CfM332A2JlL8ww9QWRNALOfPoygKKpUKlUqF7/33dUxNVhv/O3ORhANZbD6STWFpzR2gQd6u3DLU3lMU3dO7/UORzQZnttpvzz/2Fdiqvh+dlz0QxSyEYLkBQwhxfWl2QPryyy8dPxuNRry9a/4VabVa+fbbb+nVq1ebFidEWzj3wIOUpaY6nusHDbLffRYbi75v3w6rw2pT+DHtEgkHM9l8OJuLxTXznPy89PZQFBXEiLDuqNUd0HtlyoKUj2DfOij4qaY9ZKQ9FA25DfTtPL9JCCE6qWYHpLlz5wKgUqlYuHBhnWMuLi706tWLP//5z21anBAtYcnJxZyURNGWLfR8+y3HPmde06eBSoVXXByG2OnoOjDI22wK+9LzSTiYxaZDWeSayx3HfDx0zIgMZGZUEDf09kXTEaHIZoVT39jvRDthBKWqN1jvDdF32u9EC4xs/zqEEKKTa3ZAstnsd9D07t2bPXv20KOHrHEinM+SlYU5KQnTZiOl+/dTvXx00Y4dGOLjAfD9v/+jxy9+0WE1KYrCwfOFJBzM5KuDWWQWljmOGVy1VaEomJv6+qLVdNACigXnYP9HsH8dmGoNhYeNsYeiwXNA145LBAghRBfT4jlIaWlp7VGHEC1SeugwOX/4A6UHDtRpdxs2zL5w44gRjraOWKtIURRSs0wkHMziq4NZpF8qcRzz1GuZPjiAmVFBjO/vh07bQaHIarH3EiWvsfcaUbWgpJuPfSHHmIXgN7BjahFCiC6mVbf5FxcXs337dtLT06motV4MwKOPPtomhQlRW8W5c9hKS3EdMAAATfdu9nCkUuE2YoRjTpFLYGCH1nUyx8zGA/ZNYc/kFTva3Vw0TI3wZ2ZUMJMG+uHq0oEbs15Ks9+FlvIxFOXUtPcab78TLWIWaNtpPzYhhLhGqBRFUa58Wo39+/dz8803U1JSQnFxMT4+PuTl5eHu7o6/v/91c5u/yWTC29ubwsJCDAaDs8u5JlWcPYvJmIjZaKQsNRXPyZMJXbXScbxwYwLuo0fjEuDfoXWl5RWTUBWKjufUbFKr06qZMtCfmdFBTBnkj7uuVf/+aJ3KCjiWYL8T7cy2mnYPP/smsSPuBd+Om5AuhBCdVXN/f7f4b/AnnniCWbNmsWrVKrp168b333+Pi4sL99xzD4899thVFS1E+ZkzmI1GTMZEyo8dqzmgVqPYrI5b8gG8Z83ssLrOXSoh4WAWCQczOZJpcrS7aFRM6O/HzOggpkUE4OXq0mE1AZB3CvatgZR/QUleVaMK+k6xD6ENiAetrmNrEkKIa0CLA1JKSgrvvfceGo0GjUZDeXk5ffr04bXXXmPhwoXcdttt7VGnuE5kPf8CpcnJ9icaDR433ohXXCxe06ah9fHp2FoKS/nqYBYJB7NIOVfgaNeoVYzt14OZUUHEDQ7E272DQ5GlzL66dfKH8NN3Ne1eQTD8Hhi+ALqHN/16IYQQV9TigOTi4uL4F3xAQADp6elERETg7e1Nenp6mxcorj2KolB+/Dgmo5Gib74h7MMPHeHHe+YtqD09MMTG4TllMtruHbvHV665jK8PZZNwMJM9Z/Md7WoV3NDbl1nRwcyIDMTHwwm9MrlH7aHowL+grMDeplJD/1j7nWj9Y0HTgcN6QghxDWvx36bDhw9n7969DBgwgMmTJ/PCCy+Ql5fHunXrGDp0aHvUKK4BiqJQlpqK2ZiIybgZy081Ydr87bd0v/12ALrfdRfd77qrQ2u7VFzB14ezSDiQxQ9pF7HVmpU3qld3ZkYFEz80EH8v1w6tC4CKEjiywX4n2vkfa9q9Q+09RcPvAe+Qjq9LCCGucS0OSMuWLcNstk9M/f3vf8/ChQv5xS9+Qb9+/fjggw/avEDR9ZUePEjGr57Ecu6co02l0+ExYTyGuDg8J03q8JoKSywYU7NJOJjFrlN5WGulomGh3ZgZFcQtUUEEebt1eG0AZB20h6JDn0J51ZwnlQYGxtvvROs7BdQdeGecEEJcZ1oUkBRFwc/PjyFDhgDg5+fHpk2b2qUw0TUpNhtlBw+iWK24x8QA4NKzJ5aMDFSurnhOmIBhRhweEyai8fTo0NrMZRa+OZpDwoEsdpy8gMVaE4qGBBscm8KG+jhpwcRyMxz+zB6MMvfXtHfvZb8Lbdh88OrYZQyEEOJ61eKA1L9/f44cOUL//v3bqybRxSg2G6X792MyGjEnJlGZnY37yJGEf7QOAK2PD2H/+AduUUNRu3ds+CipqGTLsVw2Hshk6/ELVFTaHMcGBng5eor6+DlpzzFFgcx9Vb1Fn4Glai0ltQtEzLT3FvWaAB2w2KUQQogaLQpIarWa/v37c/HiRQlIgpLkZEybvsacmEjlhQuOdrWHBy4hwShWKyqNfRjI48YbOqyuMouVbccvkHAwk2+P5lJqsTqO9enhwcyoIGZGBzMgwKvDamqgtMA+fJb8IeQcqmn37WcPRdF3gYds5yOEEM7S4jlIr732Gr/+9a9ZtWoVkZGyqeX1pHbgAbi4+m8Ubd8OgNrLC68pU/CKi8Nj7E2o9R27UnNFpY3vTl1g44EsklJzKCqvdBwL9XFzDJ8NDjI47sLscIoC536wh6IjG6Cy1N6u0cOQufY70cJvAmfVJ4QQwqHFAemee+6hpKSE6OhodDodbm51J7FeunSpzYoTzqdYLBT/8CNmoxHzt9/S+9N/4xJiv2vKe+4cNL6+GOJi8RgzBpWuY299r7Ta2H36IgkHMzEeyaGw1OI4Fuztyi1RQcyMCiaqp7fzQhFAySU48Il9lesLtRa/9B9sD0VRd4B7x67xJIQQ4vJaHJBWrFjRDmWIzkSpqKD4++8xbbaHIlthoeOY+dtv8bn3XgAM8fEY4uM7tDarTeGHtIskHMxi8+FsLhXX7AXo56XnlqFBzIoOYnhod9RqJ4YiRYGz39lDUeqXYC23t7u4w5Db7Ktc9xwlvUVCCNFJtTggLVy4sD3qEJ1E6YEDpD/4EDZTzXYaGh8fvKZPxxAXi/vo0R1ek82msC89n40HMtl0OJsL5nLHMR8PHfGRgcyMCmZ0bx80zgxFAEUX4MA/7ZvFXjxV0x4YZQ9FQ28HV2/n1SeEEKJZrvtld0tKSoiIiOD222/n9ddfd3Y5HcpWVkbxrl2gUuE1ZQoAur79UMrL0fj1wDB9Ol5xM3AfGVNn7lFHUBSFA+cLSTiQyVeHssgqLHMcM7hqmREZyKzoYMb08UWrcfIdXjYbpG2z34l2bBPYqob6dJ4w9Gf2SdfBw51YoBBCiJa67gPSK6+8wg03dNwdVs5mKy2laMdOzEYjRdu2YSspQR8R4QhIGk8Pen/2H3S9ezslFB3JNJFwMIuvDmVy7lKp45inXkvs4ABmRgcxrp8fOm0nuO3dnA37P7L3FhX8VNMeEmOfWxQ5D/ROWj5ACCHEVbmuA9LJkyc5duwYs2bN4vDhw84up12Zv/2Wwi83UrRjB0ppTfDQBgXhMXo0SmUlKq39Pwd9v34dWtuJHDMJBzJJOJjFmbxiR7ubi4ZpgwOYGRXExAF+uLp0gpWjbVY49a29t+jEZlCqlhDQe0P0nfZgFCh3dwohRFfXaQPSjh07+NOf/kRycjJZWVls2LCBuXPn1jln5cqV/OlPfyIrK4shQ4awYsUKxo8f3+z3ePLJJ/nTn/7E7t2727h657MWFaP2cHfcvWX6ahNmoxEAl5AQvGbEYYiLw3XoUKfc4XXmQhEJB7NIOJjJiZwiR7teq2byQH9mRQczZZA/brpOEIoACs7Ze4v2fwSm8zXtYWPsoWjwHNA5aQVuIYQQba5FAamyshJXV1dSUlLafQ2k4uJioqOjue+++5g3b16D4+vXr+fxxx9n5cqVjB07lvfee4/4+HhSU1MJCwsDICYmhvLy8gavTUxMZM+ePQwYMIABAwZcMwHJajJRtHUrJmMixd99R6//fIrrgAEAeM+7zR6M4uJwHTLYKaHo3KUSEg5msfFAJqlZNZPAXTQqJg7wY2ZUMNMGB+Cp7yS53WqBE0b7nWgnk4CqrUncukP03fbtP/wHObVEIYQQ7aNFv4m0Wi3h4eFYrdYrn3yV4uPjib/MLeRvvPEGixcv5oEHHgDsyw8YjUZWrVrF8uXLAUhOTm7y9d9//z2ffPIJn376KUVFRVgsFgwGAy+88EKj55eXl9cJW6Zad3k5k7WgAPOWrZiMmyne/T+w1KwFVLxrtyMgeY4di+fYsR1eX2ZBKZsOZbHxYBYHzhU42jVqFeP69WBmVBCxQwLxdnPp8NqalH/WPq9o/8dQlF3T3mu8fcL1oJng4uqs6oQQQnSAFv9T/bnnnuPpp5/mo48+wsfHOYvbVVRUkJyczFNPPVWnPTY2ttm9QcuXL3cEqTVr1nD48OEmw1H1+S+99FLri24HpQcOcHb+PVBZs2q0rl9fDLFxeMXFoR/gnO1gcs1lbDqYRcLBLPb+lO9oV6vgxj6+zIwKZkZkID4eHbuw5GVVVsDxr+yrXJ/ZWtPu4QfD7rYPo/n2dV59QgghOlSLA9Jbb73FqVOnCA4OJjw8HA+Pujuy79u3r82Ka0peXh5Wq5WAgIA67QEBAWRnZzfxqqvz9NNPs3TpUsdzk8lEaGhou7xXYyovXsSc9A0qFy3dqoYcXSMiULu54RIUhFdcLIa4uA6fYF3tUnEFXx+2D5/9kHYJRak5NrqXDzOjg5gRGYi/Vyfreck7ZR9CS/knlORVNaqg72R7b9GAeNB2oiAnhBCiQ7Q4INWfKO1M9efRKIrSqrk1ixYtuuI5er0efQfvL2bJzcX8zTeYNxsp2bsXbDZcQkPxvu02VCoVKp2Ovpu/Ruvr26F1VSsssWA8ks3Gg5nsPn0Rq60mFQ0L7cbMqCBuiQoiyNvtMldxAksZHN1ovxPtp+9q2j0DYfg9MGIBdO/lrOqEEEJ0Ai0OSC+++GJ71NEiPXr0QKPRNOgtys3NbdCr1BUVfL6Bgs8/ozR5H7W7YlwjI/GKi7XPM6ra96yjw5G5zMI3R3PYeCCLnScvYLHW1BcZYmBmVDC3DA0i1KcT3tGVe8zeW3TgX1BaNfSnUkP/WPsQWv9Y0HSSCeJCCCGcqtW/DZKTkzl69CgqlYrBgwczfHjHrRSs0+mIiYkhKSmJW2+91dGelJTEnDlzOqyO9lJ64ACle+0TzN2io/GKi8MrNhZdzxCn1FNSUcm3R3NJOJjJ1uMXqKi0OY4NCvSq6ikKpncPj8tcxUkqSiD1C3tv0bkfatoNPe13oQ2/B7yd870KIYTovFockHJzc/n5z3/Otm3b6NatG4qiUFhYyOTJk/nkk0/w8/Nrk8KKioo4dapmL6u0tDRSUlLw8fEhLCyMpUuXsmDBAkaOHMmYMWNYvXo16enpLFmypE3e35m6/Wwe+j698Zo+HZfgYKfUUGaxsu14LhsPZrHlaC6llpo7F/v4eTAzKphZUUH0D/BySn1XlHXQ3lt08FMor9psV6WBgfH2uUV9p4C6k6yxJIQQotNRKUrt6bRXduedd3L69GnWrVtHREQEAKmpqSxcuJB+/frxr3/9q00K27ZtG5MnT27QvnDhQtasWQPYF4p87bXXyMrKIjIykjfffJMJEya0yftficlkwtvbm8LCQgwGQ4e8Z3urqLSx8+QFEg5mkZSaQ1F5zd1xYT7uzIwKYmZUMBFBXk5ZR+mKys1w+DP7nWiZtW4W6N7L3ls0bD54BTqtPCGEEM7X3N/fLQ5I3t7efPPNN4waNapO+48//khsbCwFBQWtKriruVYCksVqY/fpiyQcyMR4JBtTWU0oCvZ25ZaoIGZFBzM0xLtzhiJFsYeh5A/t4aiialVutQtEzLTPLeo9EdSdYO82IYQQTtfc398tHmKz2Wy4uDRc1M/FxQWbzdbIK0RnY7Up/HDmIhsPZrH5cBb5JTWLS/p76bl5aBCzooMYHtodtboThiKAskI4+G/7MFr2oZp23372UDTsbvDo4bz6hBBCdGktDkhTpkzhscce41//+hfBVfNjMjIyeOKJJ5g6dWqbFyjahs2mkJyeT8KBTDYdzuaCuWZVcF8PHfFDA5kZFcyoXj5oOmsoUhQ496M9FB3+HCqrNt3V6O17ocUsgvCboDP2dAkhhOhSWhyQ3nnnHebMmUOvXr0IDQ1FpVKRnp7O0KFD+eijj9qjRtFKiqKQcq6AhINZbDqURVZhmeOYt5sLM4YEMjM6iDF9fNFqOvEQVMklOLjePox24WhNu1+EPRRF3QHuzlnVXQghxLWpxQEpNDSUffv2kZSUxLFjx1AUhcGDBzNt2rT2qE+0kKIoHMk0kXAwi4SDmZzPL3Uc89JrmT4kgFlRwYzt1wOdthOHIkWBn3bZQ1Hqf8Fa1eOldYPI2+zBqOco6S0SQgjRLloUkCorK3F1dSUlJYXp06czffr09qpLtNDxbDMJBzNJOJhFWl6xo91dp2FaRAAzo4KYMMAPV5dOfmt7cZ592499H8LFmmUeCBxqD0VDbwdXb6eVJ4QQ4vrQooCk1WoJDw/HarVe+WTR7k5fKCLhgL2n6GRukaNdr1UzZZA/M6OCmTLIHzddJw9FNhukbbP3Fh37CmxVk8Z1njD0Z/ZJ18HDpbdICCFEh2nxENtzzz3H008/zUcffYSPj8z7aEslFZUMfsEIQOrLcbjrGv7xpF8sIeFQJgkHskjNMjnadRo1Ewb4MSs6iKkRAXjqu8CWGeZs2P8R7FsLBT/VtIfE2ENR5DzQezqvPiGEENetFv8Wfeuttzh16hTBwcGEh4fj4VF3e4l9+/Y18UrRWpkFpXxVNafowPlCR7tWrWJsvx7MjAoidkgg3m4Nl1/odGxWOPWtfQjt+NegVPVG6r3tk61jFtqH04QQQggnanFAmjt3bjuUIeq7YC5ny7FzJBzMIvmnfEe7WgVj+voyMyqYGUMC6e6hc2KVLVB4vqq3aB2Yzte0h95on1s0eA7oOuEGt0IIIa5LLZ6kDXD//fcTGhraLgUJu0mvb6N6jXOVCkaF+zArOogZkUH4eemdW1xzWSvhpNE+t+hUEihVC4m6dYfou+zDaP6DnFujEEII0YgWT9J+/fXXWbhwYXvVI6ooCgwP68bMqGBuGRpEoLers0tqvvyz9p6i/R9BUXZNe6/x9t6iQTPBpQt9HiGEENedFg+xTZ06lW3btrFo0aJ2KEdUS1o6gf7+Xs4uo/kqK+D4JkheA2e2AVXdX+49YPh8e2+Rb18nFiiEEEI0X4sDUnx8PE8//TSHDx8mJiamwSTt2bNnt1lx17OQbm7OLqF5Lp62T7je/zGU5NW0951iD0UDbwZtF5knJYQQQlRpcUD6xS9+AcAbb7zR4JhKpZI1kq4HljI4lmDvLTq7s6bdMxCG3wMjFkD3Xs6qTgghhLhqLQ5INputPeoQXUHuMXtv0YF/QWnVnXUqNfSbbr89v38caLrA+ktCCCHEFVzVb7OysjJcXWWy7TWtogRSv7DfiXbu+5p2Q097T9Hwe8C7p9PKE0IIIdpDiwOS1Wpl2bJlvPvuu+Tk5HDixAn69OnD888/T69evVi8eHF71HldcNdpOfvqLc4uwy77kD0UHfw3lFctTqnSwMB4+51ofaeAupNvYSKEEEK0UosD0iuvvMKHH37Ia6+9xoMPPuhoHzp0KG+++aYEpK6svAgOf2YfRstIrmnvFm4fQhs2H7wCnVefEEII0UFaHJDWrl3L6tWrmTp1KkuWLHG0R0VFcezYsTYtTnQARYHM/fZQdOg/UFG16a3aBSJm2u9E6z0R1Grn1imEEEJ0oBYHpIyMDPr169eg3WazYbFY2qQo0QHKCuHQp/Y70bIP1bT79rOHoui7wNPPaeUJIYQQztTigDRkyBB27txJeHh4nfZPP/2U4cOHt1lhoh0oCpzfY59bdORzsJTY2zV6+15oMQshfKx9bxMhhBDiOtbigPTiiy+yYMECMjIysNlsfP755xw/fpy1a9eSkJDQHjWKq1VyyT7ZOnkNXDha0+4XYQ9FUXeCu4/TyhNCCCE6mxYHpFmzZrF+/XqWLVuGSqXihRdeYMSIEWzcuJHp06e3R42iNRQFftpl7y1K/S9Yy+3tWjeIvM0+jBY6WnqLhBBCiEaoFKV6z3jREiaTCW9vbwoLCzEYDG1z0YpiWBZs//mZTNB5XP78xhTnQco/Yd9auHiypj1wqP32/KG3g6t3m5QrhBBCdDXN/f0tyx5fC2w2SNtuvxPtaALYqibL6zwhcp49GAUPl94iIYQQopkkIHVl5mxI+djeW5R/tqY9eIQ9FEXeBnovZ1UnhBBCdFkSkLoamxVOb7FPuD7+NShVmwPrvSHqDvuk68ChTi1RCCGE6OokIHUVhRmw/yPYvw4Kz9W0h95oD0WD54LO3WnlCSGEENcSCUidmbUSTiba5xadTATFZm93625fyHHEveAf4dwahRBCiGtQswLS0qVLm33BN954o9XFiFq2vwYH14M5q6at13j77fkRs8DF1Xm1CSGEENe4ZgWk/fv3N+tiKrlLqu3sWmH/X/ceMOxuezDq0XCLFyGEEEK0vWYFpK1bt7Z3HaK+3hNg5P0w8BbQ6pxdjRBCCHFdkTlIndVdn7RuoUghhBBCXLVWBaQ9e/bw6aefkp6eTkVFRZ1jn3/+eZsUJoQQQgjhLOqWvuCTTz5h7NixpKamsmHDBiwWC6mpqWzZsgVvb9nCQgghhBBdX4sD0rJly3jzzTdJSEhAp9Pxl7/8haNHj3LHHXcQFhbWHjUKIYQQQnSoFgek06dPc8sttwCg1+spLi5GpVLxxBNPsHr16jYvUAghhBCio7U4IPn4+GA2mwEICQnh8OHDABQUFFBSUtK21QkhhBBCOEGLJ2mPHz+epKQkhg4dyh133MFjjz3Gli1bSEpKYurUqe1RoxBCCCFEh2pxQHrnnXcoKysD4Omnn8bFxYXvvvuO2267jeeff77NC2xPaWlp3H///eTk5KDRaPj+++/x8JBb64UQQojrnUpRFMXZRTjLxIkT+cMf/sD48eO5dOkSBoMBrbZ5mdFkMuHt7U1hYSEGg6FtCqoohmXB9p+fyZR1kIQQQog21tzf3y2eg7Rp0yaMRmOD9sTERL7++uuWXs5pjhw5gouLC+PHjwfsc6uaG46EEEIIcW1rcUB66qmnsFqtDdptNhtPPfVUmxQFsGPHDmbNmkVwcDAqlYovvviiwTkrV66kd+/euLq6EhMTw86dO5t9/ZMnT+Lp6cns2bMZMWIEy5Yta7PahRBCCNG1tbjL5OTJkwwePLhB+6BBgzh16lSbFAVQXFxMdHQ09913H/PmzWtwfP369Tz++OOsXLmSsWPH8t577xEfH09qaqpjPaaYmBjKy8sbvDYxMRGLxcLOnTtJSUnB39+fGTNmMGrUKKZPn95mn0EIIYQQXVOLA5K3tzdnzpyhV69eddpPnTrVphOc4+PjiY+Pb/L4G2+8weLFi3nggQcAWLFiBUajkVWrVrF8+XIAkpOTm3x9z549GTVqFKGhoQDcfPPNpKSkNBmQysvL64Qtk8nU4s8khBBCiK6hxUNss2fP5vHHH+f06dOOtlOnTvGrX/2K2bNnt2lxTamoqCA5OZnY2Ng67bGxsezevbtZ1xg1ahQ5OTnk5+djs9nYsWMHERERTZ6/fPlyvL29HY/qYCWEEEKIa0+LA9Kf/vQnPDw8GDRoEL1796Z3795ERETg6+vL66+/3h41NpCXl4fVaiUgIKBOe0BAANnZ2c26hlarZdmyZUyYMIGoqCj69+/PzJkzmzz/6aefprCw0PE4d+7cVX0GIYQQQnRerRpi2717N0lJSRw4cAA3NzeioqKYMGFCe9R3WSqVqs5zRVEatF3OlYbxatPr9ej1+hbVJ4QQQoiuqVX3tatUKmJjYxsMcXWUHj16oNFoGvQW5ebmNuhVEkIIIYRoqWYFpLfeeouHHnoIV1dX3nrrrcue++ijj7ZJYZej0+mIiYkhKSmJW2+91dGelJTEnDlz2v39hRBCCHFta1ZAevPNN5k/fz6urq68+eabTZ6nUqnaLCAVFRXVWTYgLS2NlJQUfHx8CAsLY+nSpSxYsICRI0cyZswYVq9eTXp6OkuWLGmT9xdCCCHE9atZASktLa3Rn9vT3r17mTx5suP50qVLAVi4cCFr1qzhzjvv5OLFi7z88stkZWURGRnJpk2bCA8P75D6hBBCCHHtuq73YrsashebEEII0fU09/d3iydpV/fk1KdSqXB1daVfv37MmTMHHx+fll5aCCGEEKJTaHFA2r9/P/v27cNqtTJw4EAUReHkyZNoNBoGDRrEypUr+dWvfsV3333X6JYkQgghhBCdXYsXipwzZw7Tpk0jMzOT5ORk9u3bR0ZGBtOnT+euu+4iIyODCRMm8MQTT7RHvUIIIYQQ7a7Fc5BCQkJISkpq0Dt05MgRYmNjycjIYN++fcTGxpKXl9emxXYmMgdJCCGE6Hqa+/u7xT1IhYWF5ObmNmi/cOGCYwPXbt26UVFR0dJLCyGEEEJ0Cq0aYrv//vvZsGED58+fJyMjgw0bNrB48WLmzp0LwI8//siAAQPaulYhhBBCiA7R4kna7733Hk888QQ///nPqaystF9Eq2XhwoWORSQHDRrE+++/37aVCiGEEEJ0kFavg1RUVMSZM2dQFIW+ffvi6enZ1rV1ajIHSQghhOh62m0dpGqenp74+PigUqmuu3AkhBBCiGtbi+cg2Ww2Xn75Zby9vQkPDycsLIxu3brx+9//HpvN1h41CiGEEEJ0qBb3ID377LP8/e9/59VXX2Xs2LEoisKuXbv43e9+R1lZGa+88kp71CmEEEII0WFaHJA+/PBD3n//fWbPnu1oi46OJiQkhIcfflgCkhBCCCG6vBYPsV26dIlBgwY1aB80aBCXLl1qk6KEEEIIIZypxQEpOjqad955p0H7O++8Q3R0dJsUJYQQQgjhTC0eYnvttde45ZZb+OabbxgzZgwqlYrdu3dz7tw5Nm3a1B41CiGEEEJ0qBb3IE2cOJETJ05w6623UlBQwKVLl7jttts4fvw448ePb48ahRBCCCE6VKvWQQoODm4wGfvcuXPcf//9/OMf/2iTwoQQQgghnKXFPUhNuXTpEh9++GFbXU4IIYQQwmnaLCAJIYQQQlwrJCAJIYQQQtQjAUkIIYQQop5mT9K+7bbbLnu8oKDgamsRQgghhOgUmh2QvL29r3j83nvvveqChBBCCCGcrdkB6YMPPmjPOoQQQgghOg2ZgySEEEIIUY8EJCGEEEKIeiQgCSGEEELUIwFJCCGEEKIeCUhCCCGEEPVIQBJCCCGEqEcCkhBCCCFEPRKQhBBCCCHqkYAkhBBCCFGPBCQhhBBCiHokIAkhhBBC1CMBSQghhBCiHglIQgghhBD1SEASQgghhKhHApIQQgghRD3XdUB68803GTJkCIMHD+bRRx9FURRnlySEEEKITuC6DUgXLlzgnXfeITk5mUOHDpGcnMz333/v7LKEEEII0QlonV2AM1VWVlJWVgaAxWLB39/fyRUJIYQQojPotD1IO3bsYNasWQQHB6NSqfjiiy8anLNy5Up69+6Nq6srMTEx7Ny5s9nX9/Pz48knnyQsLIzg4GCmTZtG37592/ATCCGEEKKr6rQBqbi4mOjoaN55551Gj69fv57HH3+cZ599lv379zN+/Hji4+NJT093nBMTE0NkZGSDR2ZmJvn5+SQkJHD27FkyMjLYvXs3O3bs6KiPJ4QQQohOrNMOscXHxxMfH9/k8TfeeIPFixfzwAMPALBixQqMRiOrVq1i+fLlACQnJzf5+k8//ZR+/frh4+MDwC233ML333/PhAkTGj2/vLyc8vJyx3OTydTizySEEEKIrqHT9iBdTkVFBcnJycTGxtZpj42NZffu3c26RmhoKLt376asrAyr1cq2bdsYOHBgk+cvX74cb29vxyM0NPSqPoMQQgghOq8uGZDy8vKwWq0EBATUaQ8ICCA7O7tZ17jxxhu5+eabGT58OFFRUfTt25fZs2c3ef7TTz9NYWGh43Hu3Lmr+gxCCCGE6Lw67RBbc6hUqjrPFUVp0HY5r7zyCq+88kqzztXr9ej1+hbVJ4QQQoiuqUv2IPXo0QONRtOgtyg3N7dBr5IQQgghREt1yYCk0+mIiYkhKSmpTntSUhI33XSTk6oSQgghxLWi0w6xFRUVcerUKcfztLQ0UlJS8PHxISwsjKVLl7JgwQJGjhzJmDFjWL16Nenp6SxZssSJVQshhBDiWtBpA9LevXuZPHmy4/nSpUsBWLhwIWvWrOHOO+/k4sWLvPzyy2RlZREZGcmmTZsIDw93VslCCCGEuEaoFNmhtVVMJhPe3t4UFhZiMBiaPM9qtWKxWJp30YoSWF21DtNDO0Dn3gaVis7KxcUFjUbj7DKEEOK60tzf3522B6mrUxSF7OxsCgoKWvAiG4z9s/3n81mg6pJTxEQLdOvWjcDAwBbdfSmEEKL9SUBqJ9XhyN/fH3d39+b9ArRZIa9qte4evUAtvQvXKkVRKCkpITc3F4CgoCAnVySEEKI2CUjtwGq1OsKRr69v819os4K2Kki5ukpAusa5ubkB9uUp/P39ZbhNCCE6ERnDaQfVc47c3WUOkbi86v9Gmj1PTQghRIeQgNSOZF6JuBL5b0QIITonCUhCCCGEEPVIQBJCCCGEqEcCkhBCCCFEPRKQhBBCCCHqkYAkHNauXYuvry/l5eV12ufNm8e9997bomtduHCBwMBAli1b5mj74Ycf0Ol0JCYmtkm9QgghRHuRgNSBSiwlTT7KreX1zi1t8tyyyrJmXbelbr/9dqxWK19++aWjLS8vj4SEBO677z527tyJp6fnZR/VgcjPz49//OMf/O53v2Pv3r0UFRVxzz338PDDDxMbG9uKb08IIYToOLJQZAe64Z83NHlsfMh4Vk552/F80qeTKbWWNXruyICRfDDjA8fzGZ/NIL88v8F5hxYealF9bm5u3H333XzwwQfcfvvtAHz88cf07NmTSZMmUVZWRkpKymWv4ePj4/j55ptv5sEHH2T+/PmMGjUKV1dXXn311RbVJIQQQjiDBCRRx4MPPsioUaPIyMggJCSEDz74gEWLFqFSqXBzc6Nfv34tut7rr79OZGQk//73v9m7dy+urq7tVLkQQgjRdiQgdaAf7v6hyWOaetuKbLt9K6gbHwFV19vEdvO8zVdfXJXhw4cTHR3N2rVriYuL49ChQ2zcuBGAnTt3Eh8ff9nXP/PMMzzzzDOO52fOnCEzMxObzcZPP/1EVFRUm9UqhBBCtBcJSB3I3eUKW4/YrLXOdWv2XmxXvG4LPfDAA7z55ptkZGQwbdo0QkNDARg5cmSLhtgqKiqYP38+d955J4MGDWLx4sUcOnSIgICANq1XCCGEaGsSkEQD8+fP58knn+Rvf/sba9eudbS3dIjt2WefpbCwkLfeegtPT0++/vprFi9eTEJCQnuULYQQQrQZuYtNNGAwGJg3bx6enp7MnTu3VdfYtm0bK1asYN26dRgMBtRqNevWreO7775j1apVbVuwEEII0cakB0k0Kisri/nz56PX61v1+kmTJjXYoT4sLIyCgoI2qE4IIYRoXxKQRB2XLl0iMTGRLVu28M477zi7HCGEEMIpJCCJOkaMGEF+fj5//OMfGThwoLPLEUIIIZxCApKo4+zZs84uQQghhHA6maQthBBCCFGPBCQhhBBCiHokIAkhhBBC1CMBSQghhBCiHglIQgghhBD1SEASQgghhKhHApLoFE6ePElAQADu7u7s2rXL2eUIIYS4zklAEk6XmZlJbGws48aNY/HixcycOZNDhw41OK+oqIhf/vKX9OzZEzc3NyIiImRfNyGEEO1CFooUTpWfn+8IR2vWrEGj0eDl5UVcXBzfffcdffr0cZz7xBNPsHXrVj766CN69epFYmIiDz/8MMHBwcyZM8eJn0IIIcS1RnqQhMPatWvx9fWlvLy8Tvu8efO49957W3StCxcuEBgYyLJlyxxtP/zwAzqdjsTERABKSkq4+eabGTduHB9++CEajQaAZcuW8cgjjxAbG0t2drbj9f/73/9YuHAhkyZNolevXjz00ENER0ezd+/e1n5kIYQQolEqRVEUZxfRFZlMJry9vSksLMRgMNQ5VlZWRlpaGr1798bV1dXRbispafqCGg1qFy1kH7Sfa+gHak3j56rVqJtxXbW7ezM/jV1paSlBQUH87W9/4/bbbwcgLy+PkJAQNm/ejFarJT4+/rLXeOaZZ3jmmWcA2LRpE3PnzmX37t0MGjSI4cOHc8stt7BixYoW1VVtyZIlJCcn88UXXxAcHMy2bduYPXs2X3/9NePGjWvVNZ2tqf9WhBBCtI/L/f6uTYbYOtDxETFNHvOYOIGwVSsdz0+Mn4hSWtroue6jRhG+bq3j+amp07Dm5zc4L+LY0RbV5+bmxt13380HH3zgCEgff/wxPXv2ZNKkSZSVlZGSknLZa/j4+Dh+vvnmm3nwwQeZP38+o0aNwtXVlVdffbVFNdX21ltv8eCDD9KzZ0+0Wi1qtZr333+/y4YjIYQQnZcEJFHHgw8+yKhRo8jIyCAkJIQPPviARYsWoVKpcHNzo1+/fi263uuvv05kZCT//ve/2bt371X1krz11lt8//33fPnll4SHh7Njxw4efvhhgoKCmDZtWquvK4QQQtQnAakDDdyX3PRBTd3htAE7t192iK22ft9+c7WlOQwfPpzo6GjWrl1LXFwchw4dYuPGjQDs3LmzRUNsAGfOnCEzMxObzcZPP/1EVFRUq+oqLS3lmWeeYcOGDdxyyy0AREVFkZKSwuuvvy4BSQghRJuSgNSBrjgnyGate25TAaml122hBx54gDfffJOMjAymTZtGaGgoACNHjmzREFtFRQXz58/nzjvvZNCgQSxevJhDhw4REBDQ4posFgsWiwV1vXCo0Wiw2Wwtvp4QQghxORKQRAPz58/nySef5G9/+xtr19bMdWrpENuzzz5LYWEhb731Fp6ennz99dcsXryYhISEFtdkMBiYOHEiv/71r3FzcyM8PJzt27ezdu1a3njjjRZfTwghhLgcuc1fNGAwGJg3bx6enp7MnTu3VdfYtm0bK1asYN26dRgMBtRqNevWreO7775r9eKOn3zyCaNGjWL+/PkMHjyYV199lVdeeYUlS5a06npCCCFEU6QHSTQqKyuL+fPno9frW/X6SZMmYbFY6rSFhYVRUFDQ6poCAwP54IMPWv16IYQQormuix6kW2+9le7du/Ozn/2swbGEhAQGDhxI//79ef/9951QXedy6dIlPvnkE7Zs2cIjjzzi7HKEEEIIp7guepAeffRR7r//fj788MM67ZWVlSxdupStW7diMBgYMWIEt912W52JxtebESNGkJ+fzx//+EcGDhzo7HKEEEIIp7guAtLkyZPZtm1bg/Yff/yRIUOGEBISAtgXNjQajdx1110dXGHncfbsWWeXIIQQQjid04fYduzYwaxZswgODkalUvHFF180OGflypWOrRhiYmLYuXNnm7x3ZmamIxwB9OzZk4yMjDa5thBCCCG6LqcHpOLiYqKjo3nnnXcaPb5+/Xoef/xxnn32Wfbv38/48eOJj48nPT3dcU5MTAyRkZENHpmZmZd978a2oVOpVFf3gYQQQgjR5Tl9iC0+Pv6yqzO/8cYbLF68mAceeACAFStWYDQaWbVqFcuXLwcgOfkyK1RfRkhISJ0eo/Pnz3PDDTc0em55eXmdXe5NJlOr3lMIIYQQnZ/Te5Aup6KiguTkZGJjY+u0x8bGsnv37qu+/ujRozl8+DAZGRmYzWY2bdpEXFxco+cuX74cb29vx6N6dWkhhBBCXHs6dUDKy8vDarU22JoiICCA7OzsZl8nLi6O22+/nU2bNtGzZ0/27NkDgFar5c9//jOTJ09m+PDh/PrXv8bX17fRazz99NMUFhY6HufOnWv9BxNCCCFEp+b0IbbmqD8vSFGUFs0VMhqNTR6bPXs2s2fPvuI19Hp9qxdNbK2SikoGv2CvPfXlONx1XeKPSwghhOjyOnUPUo8ePdBoNA16i3Jzc1u14akQQgghRHN06oCk0+mIiYkhKSmpTntSUhI33XSTk6q6ti1atKjV+6/V99hjjxETE4Ner2fYsGGNnqMoCq+//joDBgxAr9cTGhrKsmXL2uT9m2vNmjWoVKoGj7KyMsc5V1qOwmKx8Nvf/pahQ4fi4eFBcHAw99577xXvpBRCCNE5OX3MpqioiFOnTjmep6WlkZKSgo+PD2FhYSxdupQFCxYwcuRIxowZw+rVq0lPT5cNSrsARVG4//77+eGHHzh48GCj5zz22GMkJiby+uuvM3ToUAoLC8nLy+vgSu0b9B4/frxOm6urq+Pn6uUo7rvvPubNm9fg9SUlJezbt4/nn3+e6Oho8vPzefzxx5k9ezZ79+5t9/qFEEK0LacHpL179zJ58mTH86VLlwKwcOFC1qxZw5133snFixd5+eWXycrKIjIykk2bNhEeHu6skkUzvfXWWwBcuHCh0YB09OhRVq1axeHDh5u9rcnZs2fp3bs3n332GW+//TY//PAD/fv3591332XMmDGtrlWlUhEYGNjk8SstR+Ht7d2gp/Ptt99m9OjRpKenExYW1urahBBCdDynB6RJkyY1umBjbQ8//DAPP/xwB1XU9hRFodRivfKJNitYbPafKyopqaz5XkoqKlv13m4umjZb/DI+Pv6Kq5gXFRU1+3obN26kT58+JCQkMGPGDBRFYdq0abz22mtX3A/v2Wef5fXXX6d///48++yz3HXXXZw6dQqtVkt6ejqDBw++7Ovvuece3n333Tp1h4eHY7VaGTZsGL///e8ZPnx4sz9LYwoLC1GpVHTr1u2qriOEEKLjOT0gXQ9KLVbH3WjNV3di+sg/fNuq927Lu9/ef/99SktL2+RaAGfOnOGnn37i008/Ze3atVitVp544gl+9rOfsWXLlsu+9sknn+SWW24B4KWXXmLIkCGcOnWKQYMGERwcTEpKymVfbzAYHD8PGjSINWvWMHToUEwmE3/5y18YO3YsBw4coH///q36bGVlZTz11FPcfffddd5LCCFE1yABSTRb7X3r2oLNZqO8vJy1a9cyYMAAAP7+978TExPD8ePHLzvsFhUV5fg5KCgIsN/dOGjQILRaLf369Wt2HTfeeCM33nij4/nYsWMZMWIEb7/9tmOYsCUsFgs///nPsdlsrFy5ssWvF0II4XwSkDqAm4uG1JcbX6G7DpsVcg7bfw6IpKRScfQc7X1uaqt6gtxcNC1+TVPaeogtKCgIrVbrCEcAERERAKSnp182ILm4uDh+rh5CtNlsjte2dIitNrVazahRozh58mTzPkgtFouFO+64g7S0NLZs2SK9R0II0UVJQOoAKpWqeeHGpgKXqpUXdFpQ18xBctdpnb5QZFsPsY0dO5bKykpOnz5N3759AThx4gTAVU3Cb+kQW32KopCSksLQoUNb9L7V4ejkyZNs3bq1yVXZhRBCdH4SkEQDhYWFDQJG9bILLXHq1CmKiorIzs6mtLTUcc3Bgwej0+mYNm0aI0aM4P7772fFihXYbDYeeeQRpk+fXqdXqaVaOsT20ksvceONN9K/f39MJhNvvfUWKSkp/PWvf3Wcc6XlKCorK/nZz37Gvn37SEhIwGq1OhY49fHxQafTtfrzCCGE6HgSkEQD27Zta3AHV/WyCy3xwAMPsH37dsfz6mumpaXRq1cv1Go1Gzdu5P/9v//HhAkT8PDwID4+nj//+c+O11Tf1r9161YmTZrU6s90OQUFBTz00ENkZ2fj7e3N8OHD2bFjB6NHj3acc6XlKM6fP8+XX34J0GBRzPasXQghRPtQKVe6x140ymQy4e3tTWFhYYPhmrKyMtLS0ujdu3edxQavyGaF7Kr1ggKjKKlUrvu92LZt28att97KmTNn6N69u7PLaXOt/m9FCCFEq1zu93dtnXqrESE2b97MM888c02GIyGEEJ3X9dcl0YW467ScffUWZ5fhVK+++qqzSxBCCHEdkh4kIYQQQoh6JCAJIYQQQtQjAUkIIYQQoh4JSEIIIYQQ9UhAEkIIIYSoRwKSEEIIIUQ9EpA6s4pi+J23/VFR7OxqhBBCiOuGBCQhhBBCiHokIIk6Fi1axNy5c6/6OgcOHOCuu+4iNDQUNzc3IiIi+Mtf/tLgvEOHDjFx4kTc3NwICQnh5Zdfxhm736xYsYKBAwfi5uZGaGgoTzzxBGVlZXXOWblypWNLkJiYGHbu3FnnuKIo/O53vyM4OBg3NzcmTZrEkSNHOvJjCCGEaCOykrZoF8nJyfj5+fHRRx8RGhrK7t27eeihh9BoNPzyl78E7PvhTJ8+ncmTJ7Nnzx5OnDjBokWL8PDw4Fe/+lWH1frxxx/z1FNP8Y9//IObbrrJUQfAm2++CcD69et5/PHHWblyJWPHjuW9994jPj6e1NRUwsLCAHjttdd44403WLNmDQMGDOAPf/gD06dP5/jx43h5eXXY5xFCCNEGFNEqhYWFCqAUFhY2OFZaWqqkpqYqpaWlLbuotVJRMvbZH9ZKRSkvUpQXDfZHeVEbVX55CxcuVObMmdMu13744YeVyZMnO56vXLlS8fb2VsrKyhxty5cvV4KDgxWbzdboNdLS0hRA+eyzz5RJkyYpbm5uSlRUlLJ79+5W1/XII48oU6ZMqdO2dOlSZdy4cY7no0ePVpYsWVLnnEGDBilPPfWUoiiKYrPZlMDAQOXVV191HC8rK1O8vb2Vd999t8n3bvV/K0IIIVrlcr+/a5Mhto6gKPZJ1s15WErtj4piqCipuUZFSfOvUfvRhsNV8fHxeHp6XvZxOYWFhfj4+Die/+9//2PixIno9XpHW1xcHJmZmZw9e/ay13r22Wd58sknSUlJYcCAAdx1111UVlYCkJ6efsU6lyxZ4rjWuHHjSE5O5scffwTgzJkzbNq0iVtuse+DV1FRQXJyMrGxsXVqiI2NZffu3QCkpaWRnZ1d5xy9Xs/EiRMd5wghhOg6ZIitI1hKYFnw1V3j9X6te90zmaDzuLr3rvL+++9TWlraqtf+73//49///jdfffWVoy07O5tevXrVOS8gIMBxrHfv3k1e78knn3QEmJdeeokhQ4Zw6tQpBg0aRHBwMCkpKZetx2AwOH7++c9/zoULFxg3bhyKolBZWckvfvELnnrqKQDy8vKwWq2O2mrXmp2d7ai3dv21z/npp58uW4sQQojORwKSaLaQkJBWve7IkSPMmTOHF154genTp9c5plKp6jxXqnq86rfXFxUV5fg5KCgIgNzcXAYNGoRWq6Vfv+YHym3btvHKK6+wcuVKbrjhBk6dOsVjjz1GUFAQzz///GVrrd/WnHOEEEJ0fhKQOoKLu70n50psVsg5bP85IBIqy2t6jp48BTr31r13G4mPj29w51Z9RUVFdZ6npqYyZcoUHnzwQZ577rk6xwIDAx09L9Vyc3OBhj0x9bm4uDh+rg4gNpsNsA+xDR48+LKvv+eee3j33XcBeP7551mwYAEPPPAAAEOHDqW4uJiHHnqIZ599lh49eqDRaBqttbrOwMBAwN6TVB3Y6p8jhBCi65CA1BFUquYNc9ms4OJm/1nnAWpNzTGde5sNlbVWS4fYjhw5wpQpU1i4cCGvvPJKg+NjxozhmWeeoaKiAp1OB0BiYiLBwcENht5aoqVDbCUlJajVdafjaTQaFEVBURR0Oh0xMTEkJSVx6623Os5JSkpizpw5APTu3ZvAwECSkpIYPnw4YJ+7tH37dv74xz+2+rMIIYRwDglInYlaA8HDnV0FhYWFDQKGj4+P43b25jhy5AiTJ08mNjaWpUuXOnpfNBoNfn5+ANx999289NJLLFq0iGeeeYaTJ0+ybNkyXnjhhasalmrpENusWbN44403GD58uGOI7fnnn2f27NloNPaQunTpUhYsWMDIkSMZM2YMq1evJj093THZW6VS8fjjj7Ns2TL69+9P//79WbZsGe7u7tx9992t/ixCCCGcQwKSaGDbtm2OXpBqCxcuZM2aNc2+xqeffsqFCxf4+OOP+fjjjx3t4eHhjjvUvL29SUpK4pFHHmHkyJF0796dpUuXsnTpUsf5Z8+epXfv3mzdupVJkyZdzcdq0nPPPYdKpeK5554jIyMDPz8/Zs2aVafX68477+TixYu8/PLLZGVlERkZyaZNmwgPD3ec85vf/IbS0lIefvhh8vPzueGGG0hMTJQ1kIQQogtSKYoTli2+BphMJry9vSksLKwzXANQVlZGWlqaY9XlVqsorrn7rQ3vRutKtm3bxq233sqZM2fo3r27s8tpc23234oQQohmudzv79pkHSTRqW3evJlnnnnmmgxHQgghOi8ZYuvMdB7wu0JnV+FUr776qrNLEEIIcR2SHiQhhBBCiHokIAkhhBBC1CMBqR3J/HdxJfLfiBBCdE4SkNpB9SrPJSUlVzhTXO+q/xupvTK4EEII55NJ2u1Ao9HQrVs3x7YZ7u7ush+XqENRFEpKSsjNzaVbt26OBSmFEEJ0DhKQ2kn13lzVIUmIxnTr1s3x34oQQojOQwJSO1GpVAQFBeHv74/FYnF2OaITcnFxkZ4jIYTopCQgtTONRiO/BIUQQoguRiZpCyGEEELUIwFJCCGEEKIeCUhCCCGEEPXIHKRWql7gz2QyObkSIYQQQjRX9e/tKy3UKwGplcxmMwChoaFOrkQIIYQQLWU2m/H29m7yuEqRvQ5axWazkZmZiZeXV5svAjlq1Cj27NnTpte83lzP3+G19Nm70mfpjLV2hpqcUYPJZCI0NJRz585hMBg69L1F56coCmazmeDgYNTqpmcaSQ9SK6nVanr27Nku19ZoNPJ/6qt0PX+H19Jn70qfpTPW2hlqcmYNBoPB6Z9fdE6X6zmqJpO0O6FHHnnE2SV0edfzd3gtffau9Fk6Y62doabOUIMQrSFDbEIIIa4pJpMJb29vCgsLpQdJtJr0IAkhhLim6PV6XnzxRfR6vbNLEV2Y9CAJIYQQQtQjPUhCCCGEEPVIQBJCCCGEqEcC0nXm3LlzTJo0icGDBxMVFcWnn37q7JK6HPkOuz75MxRCXInMQbrOZGVlkZOTw7Bhw8jNzWXEiBEcP34cDw8PZ5fWZch32PXJn6EQ4kpkocjrTFBQEEFBQQD4+/vj4+PDpUuX5BdDC8h32PXJn6EQ4kpkiK2ZVq1aRVRUlGNl1jFjxvD111+36Xvs2LGDWbNmERwcjEql4osvvmj0vJUrV9K7d29cXV2JiYlh586drXq/vXv3YrPZnLKf3PLly1GpVDz++ONtet3r6Tt0loyMDO655x58fX1xd3dn2LBhJCcnt9n15c9QdKRbb72V7t2787Of/czZpYhORgJSM/Xs2ZNXX32VvXv3snfvXqZMmcKcOXM4cuRIo+fv2rULi8XSoP3YsWNkZ2c3+pri4mKio6N55513mqxj/fr1PP744zz77LPs37+f8ePHEx8fT3p6uuOcmJgYIiMjGzwyMzMd51y8eJF7772X1atXN/craDN79uxh9erVREVFXfY8+Q47n/z8fMaOHYuLiwtff/01qamp/PnPf6Zbt26Nni9/hqKze/TRR1m7dq2zyxCdkSJarXv37sr777/foN1qtSrR0dHKz372M6WystLRfvz4cSUwMFD54x//eMVrA8qGDRsatI8ePVpZsmRJnbZBgwYpTz31VLPrLisrU8aPH6+sXbu22a9pK2azWenfv7+SlJSkTJw4UXnssccaPU++w87pt7/9rTJu3LhmnSt/hqKr2Lp1qzJv3jxnlyE6GelBagWr1conn3xCcXExY8aMaXBcrVazadMm9u/fz7333ovNZuP06dNMmTKF2bNn85vf/KZV71tRUUFycjKxsbF12mNjY9m9e3ezrqEoCosWLWLKlCksWLCgVXVcjUceeYRbbrmFadOmXfY8+Q47py+//JKRI0dy++234+/vz/Dhw/nb3/7W6LnyZyiuVnOGW9tqqFWI+iQgtcChQ4fw9PREr9ezZMkSNmzYwODBgxs9Nzg4mC1btrBr1y7uvvtupkyZwtSpU3n33Xdb/f55eXlYrVYCAgLqtAcEBDQ5XFHfrl27WL9+PV988QXDhg1j2LBhHDp0qNU1tcQnn3zCvn37WL58ebPOl++w8zlz5gyrVq2if//+GI1GlixZctkhCvkzFFfjSsOtbTXUKkRj5C62Fhg4cCApKSkUFBTw2WefsXDhQrZv395kSAoLC2Pt2rVMnDiRPn368Pe//x2VSnXVddS/hqIozb7uuHHjsNlsV11DS507d47HHnuMxMREXF1dm/06+Q47F5vNxsiRI1m2bBkAw4cP58iRI6xatYp777230dfIn6Forfj4eOLj45s8/sYbb7B48WIeeOABAFasWIHRaGTVqlWOf4i15Q0E4voiPUgtoNPp6NevHyNHjmT58uVER0fzl7/8pcnzc3JyeOihh5g1axYlJSU88cQTV/X+PXr0QKPRNPhXcm5uboN/TXc2ycnJ5ObmEhMTg1arRavVsn37dt566y20Wi1Wq7XR18l32LkEBQU1+AdBREREnX+x1yd/hqI9tMVQqxCXIwHpKiiKQnl5eaPH8vLymDp1KhEREXz++eds2bKFf//73zz55JOtfj+dTkdMTAxJSUl12pOSkrjppptafd2OMHXqVA4dOkRKSorjMXLkSObPn09KSgoajabBa+Q77HzGjh3L8ePH67SdOHGC8PDwRs+XP0PRXtpiqBUgLi6O22+/nU2bNtGzZ0/27NnT1qWKLkqG2JrpmWeeIT4+ntDQUMxmM5988gnbtm1j8+bNDc612WzMmDGD8PBw1q9fj1arJSIigm+++YbJkycTEhLS6L+ii4qKOHXqlON5WloaKSkp+Pj4EBYWBsDSpUtZsGABI0eOZMyYMaxevZr09HSWLFnSfh++DXh5eREZGVmnzcPDA19f3wbtIN9hZ/XEE09w0003sWzZMu644w5+/PFHVq9e3eht8vJnKDrC1Qy1AhiNxrYuSVwrnHcDXddy//33K+Hh4YpOp1P8/PyUqVOnKomJiU2en5iYqJSWljZo379/v5Kent7oa7Zu3aoADR4LFy6sc95f//pXRy0jRoxQtm/fflWfzVkud5u/osh32Flt3LhRiYyMVPR6vTJo0CBl9erVTZ4rf4airVBvyYfy8nJFo9Eon3/+eZ3zHn30UWXChAkdXJ24FslebEIIITo9lUrFhg0bmDt3rqPthhtuICYmhpUrVzraBg8ezJw5c5p9t6wQTZEhNiGEEJ3SlYZbZahVtCfpQRJCCNEpbdu2jcmTJzdoX7hwIWvWrAHsC0W+9tprZGVlERkZyZtvvsmECRM6uFJxLZKAJIQQQghRj9zmL4QQQghRjwQkIYQQQoh6JCAJIYQQQtQjAUkIIYQQoh4JSEIIIYQQ9UhAEkIIIYSoRwKSEEIIIUQ9EpCEEEIIIeqRgCSEaLVt27ahUqkoKChol+ufPXsWlUpFSkpKm11TpVLxxRdftNn1hBDXJglIQogmLVq0CJVKhUqlwsXFhT59+vDkk09SXFzcIe8fGhrq2EJCtE57h1ghrlWyWa0Q4rJmzJjBBx98gMViYefOnTzwwAMUFxezatWqdn9vjUZDYGBgu7/P1VAUBavVilZb96/TiooKdDpdi6/X3Ne19vpCiOaRHiQhxGXp9XoCAwMJDQ3l7rvvZv78+Q2GqJKTkxk5ciTu7u7cdNNNHD9+HLAPkanVavbu3Vvn/Lfffpvw8HAURSE/P5/58+fj5+eHm5sb/fv354MPPnC8vv4Q25EjR7jlllswGAx4eXkxfvx4Tp8+DcCePXuYPn06PXr0wNvbm4kTJ7Jv374WfV5FUXjttdfo06cPbm5uREdH85///MdxvLpHxmg0MnLkSPR6PTt37mTSpEn88pe/ZOnSpfTo0YPp06cDsH37dkaPHo1erycoKIinnnqKyspKx/Wael19ixYtYu7cuSxfvpzg4GAGDBgAwEcffcTIkSPx8vIiMDCQu+++m9zcXMf3V73Za/fu3VGpVCxatKhZn1OI650EJCFEi7i5uWGxWOq0Pfvss/z5z39m7969aLVa7r//fgB69erFtGnTHIGn2gcffOAYvnv++edJTU3l66+/5ujRo6xatYoePXo0+t4ZGRlMmDABV1dXtmzZQnJyMvfff78jcJjNZhYuXMjOnTv5/vvv6d+/PzfffDNms7nZn++5557jgw8+YNWqVRw5coQnnniCe+65h+3bt9c57ze/+Q3Lly/n6NGjREVFAfDhhx+i1WrZtWsX7733HhkZGdx8882MGjWKAwcOsGrVKv7+97/zhz/8oc616r+uKd9++y1Hjx4lKSmJhIQEwN6T9Pvf/54DBw7wxRdfkJaW5ghBoaGhfPbZZwAcP36crKws/vKXv7Tocwpx3VKEEKIJCxcuVObMmeN4/sMPPyi+vr7KHXfcoSiKomzdulUBlG+++cZxzldffaUASmlpqaIoirJ+/Xqle/fuSllZmaIoipKSkqKoVColLS1NURRFmTVrlnLfffc1+v5paWkKoOzfv19RFEV5+umnld69eysVFRXNqr+yslLx8vJSNm7c6GgDlA0bNjR6flFRkeLq6qrs3r27TvvixYuVu+66q85n/uKLL+qcM3HiRGXYsGF12p555hll4MCBis1mc7T99a9/VTw9PRWr1drk6xqzcOFCJSAgQCkvL7/seT/++KMCKGazuU69+fn5LfqcQlzvpAdJCHFZCQkJeHp64urqypgxY5gwYQJvv/12nXOqe1AAgoKCABzDPHPnzkWr1bJhwwYA/vGPfzB58mR69eoFwC9+8Qs++eQThg0bxm9+8xt2797dZC0pKSmMHz8eFxeXRo/n5uayZMkSBgwYgLe3N97e3hQVFZGent6sz5qamkpZWRnTp0/H09PT8Vi7dq1jGK/ayJEjG7y+ftvRo0cZM2YMKpXK0TZ27FiKioo4f/78Za/VmKFDhzaYd7R//37mzJlDeHg4Xl5eTJo0CeCyn7kln1OI65VM0hZCXNbkyZNZtWoVLi4uBAcHNxpOardVhwGbzQaATqdjwYIFfPDBB9x2223885//ZMWKFY7z4+Pj+emnn/jqq6/45ptvmDp1Ko888givv/56g/dxc3O7bK2LFi3iwoULrFixgvDwcPR6PWPGjKGioqJZn7W65q+++oqQkJA6x/R6fZ3nHh4eDV5fv01RlDrhqLoNqNPe2LUaU/+84uJiYmNjiY2N5aOPPsLPz4/09HTi4uIu+5lb8jmFuF5JQBJCXJaHhwf9+vW7qms88MADREZGsnLlSiwWC7fddlud435+fixatIhFixYxfvx4fv3rXzcakKKiovjwww+xWCyNBrWdO3eycuVKbr75ZgDOnTtHXl5es+scPHgwer2e9PR0Jk6c2MJP2fj1PvvsszpBaffu3Xh5eTUIJq1x7Ngx8vLyePXVVwkNDQVoMCG+usfJarXWqastP6cQ1yIZYhNCtLuIiAhuvPFGfvvb33LXXXfV6Ql64YUX+O9//8upU6c4cuQICQkJRERENHqdX/7yl5hMJn7+85+zd+9eTp48ybp16xx3zfXr149169Zx9OhRfvjhB+bPn3/FXqfavLy8ePLJJ3niiSf48MMPOX36NPv37+evf/0rH374YYs/98MPP8y5c+f4f//v/3Hs2DH++9//8uKLL7J06VLU6qv/6zcsLAydTsfbb7/NmTNn+PLLL/n9739f55zw8HBUKhUJCQlcuHCBoqKiNv+cQlyLJCAJITrE4sWLqaiocNzhVk2n0/H0008TFRXFhAkT0Gg0fPLJJ41ew9fXly1btlBUVMTEiROJiYnhb3/7m6M36R//+Af5+fkMHz6cBQsW8Oijj+Lv79+iOn//+9/zwgsvsHz5ciIiIoiLi2Pjxo307t27xZ85JCSETZs28eOPPxIdHc2SJUtYvHgxzz33XIuv1Rg/Pz/WrFnDp59+yuDBg3n11Vcb9LyFhITw0ksv8dRTTxEQEMAvf/nLNv+cQlyLVEr1gLgQQrSjV155hU8++YRDhw45uxQhhLgi6UESQrSroqIi9uzZw9tvv82jjz7q7HKEEKJZJCAJIdrVL3/5S8aNG8fEiRP/f3t2TAMAAMMwjD/r/RmBHjaKqH33GsAqFxsAQFiQAABCIAEAhEACAAiBBAAQAgkAIAQSAEAIJACAEEgAACGQAADiAHlmwqHPkzhYAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "plt.figure()\n",
    "plt.yscale('log'),plt.xscale('log')\n",
    "for L, logical_errors in zip(Ls, log_errors_all_L):\n",
    "    std_err = (logical_errors*(1-logical_errors)/num_trials)**0.5\n",
    "    plt.errorbar(ps, logical_errors, yerr=std_err, label=\"L={}\".format(L)+\",n={}\".format(2*L*L))\n",
    "plt.plot(ps, ps, '--',label='y=x')  #add y=x by weilei\n",
    "plt.plot(ps, 100000000*ps**8, '--',label='y=x^8')  #add y=x^8 by weilei\n",
    "plt.xlabel(\"Physical error rate\")\n",
    "plt.ylabel(\"Logical error rate\")\n",
    "plt.legend(loc=0);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Noisy syndromes\n",
    "\n",
    "In the presence of measurement errors, each syndrome measurement is repeated $O(L)$ times, and decoding instead takes place over a 3D matching graph with an additional time dimension (see Section IV B of [this paper](https://arxiv.org/abs/quant-ph/0110143)). The time dimension can be added to the matching graph by specifying the number of repetitions when constructing the matching object:\n",
    "```\n",
    "matching = Matching(H, repetitions=T)\n",
    "```\n",
    "where here $T$ is the number of repetitions. For decoding, the difference syndrome should be supplied as an $r\\times T$ binary numpy matrix, where $r$ is the number of checks (rows in $H$). The difference syndrome in time step $t$ is the difference (modulo 2) between the syndrome measurement in time step $t$ and $t-1$, and ensures that any single measurement error results in two syndrome defects (at the endpoints of a timelike edge in the matching graph). The last round of syndrome measurements should be free of measurement errors to ensure that the overall syndrome has even parity: when qubits are measured individually at the end of a computation then the final round of syndrome measurement is indeed error-free (stabilisers can be determined exactly in post-processing), however the [overlapping recovery method](https://arxiv.org/abs/quant-ph/0110143) should be implemented when decoding must be completed before all qubits are measured.\n",
    "\n",
    "The following example demonstrates decoding in the presence of measurement errors using a phenomenological error model. In this error model, in each round of measurements each qubit suffers an error with probability $p$, and each syndrome is measured incorrectly with probability $q$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def num_decoding_failures_noisy_syndromes(H, logicals, p, q, num_trials, repetitions):\n",
    "    matching = Matching(H, weights=np.log((1-p)/p),\n",
    "                repetitions=repetitions, timelike_weights=np.log((1-q)/q), faults_matrix=logicals)\n",
    "    num_stabilisers, num_qubits = H.shape\n",
    "    num_errors = 0\n",
    "    for i in range(num_trials):\n",
    "        noise_new = (np.random.rand(num_qubits, repetitions) < p).astype(np.uint8)\n",
    "        noise_cumulative = (np.cumsum(noise_new, 1) % 2).astype(np.uint8)\n",
    "        noise_total = noise_cumulative[:,-1]\n",
    "        syndrome = H@noise_cumulative % 2\n",
    "        syndrome_error = (np.random.rand(num_stabilisers, repetitions) < q).astype(np.uint8)\n",
    "        syndrome_error[:,-1] = 0 # Perfect measurements in last round to ensure even parity\n",
    "        noisy_syndrome = (syndrome + syndrome_error) % 2\n",
    "        # Convert to difference syndrome\n",
    "        noisy_syndrome[:,1:] = (noisy_syndrome[:,1:] - noisy_syndrome[:,0:-1]) % 2\n",
    "        predicted_logicals_flipped = matching.decode(noisy_syndrome)\n",
    "        actual_logicals_flipped = noise_total@logicals.T % 2\n",
    "        if not np.array_equal(predicted_logicals_flipped, actual_logicals_flipped):\n",
    "            num_errors += 1\n",
    "    return num_errors"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We'll now simulate the performance of the decoder for a range of lattice sizes $L$ and physical error rate $p$ (taking $q=p$) and estimate the threshold:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Simulating L=8...\n",
      "Simulating L=10...\n",
      "Simulating L=12...\n",
      "CPU times: user 3min 58s, sys: 1.64 ms, total: 3min 58s\n",
      "Wall time: 3min 58s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "\n",
    "num_trials = 30000\n",
    "Ls = range(8,13,2)\n",
    "ps = np.linspace(0.02, 0.04, 7)\n",
    "log_errors_all_L = []\n",
    "for L in Ls:\n",
    "    print(\"Simulating L={}...\".format(L))\n",
    "    Hx = toric_code_x_stabilisers(L)\n",
    "    logX = toric_code_x_logicals(L)\n",
    "    log_errors = []\n",
    "    for p in ps:\n",
    "        num_errors = num_decoding_failures_noisy_syndromes(Hx, logX, p, p, num_trials, L)\n",
    "        log_errors.append(num_errors/num_trials)\n",
    "    log_errors_all_L.append(np.array(log_errors))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plotting the results, we find a threshold of around 3%, consistent with the threshold of 2.9% found in [this paper](https://arxiv.org/abs/quant-ph/0207088):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkQAAAGwCAYAAABIC3rIAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAA9hAAAPYQGoP6dpAAB7cUlEQVR4nO3dd3xN9x/H8dfN3pFIIkMksQWx96YoapQubRGjvw6rtFVao2jRVtVoaIsaVdVhlFZrx97EHglBhAgie957z++PKyGVkhs3uTfJ5/l49NF7zz33ez5f103ezvme71elKIqCEEIIIUQpZmbsAoQQQgghjE0CkRBCCCFKPQlEQgghhCj1JBAJIYQQotSTQCSEEEKIUk8CkRBCCCFKPQlEQgghhCj1LIxdgCnTarXcuHEDR0dHVCqVscsRQgghRD4oikJSUhLe3t6YmeXv3I8EojyEhIQQEhJCZmYmly5dMnY5QgghhCiAqKgoypcvn699VTJT9X9LSEigTJkyREVF4eTkZOxyhBBCCJEPiYmJ+Pr6Eh8fj7Ozc77eI2eIHiP7MpmTk5MEIiGEEKKY0We4iwyqFkIIIUSpJ4FICCGEEKWeBCIhhBBClHoyhsgANBoNWVlZxi6jxLC0tMTc3NzYZQghhChFJBA9BUVRiImJIT4+3tillDhlypTB09NT5n8SQghRJCQQPYXsMOTh4YGdnZ388jYARVFITU0lNjYWAC8vLyNXJIQQojSQQFRAGo0mJwyVLVvW2OWUKLa2tgDExsbi4eEhl8+EEEIUOhlUnYeQkBACAwNp1KjRf+6TPWbIzs6uqMoqVbL/XGVslhBCiKIggSgPQ4cO5ezZsxw+fPiJ+8plssIhf65CCCGKkgQiE5CaqcZ/7F/4j/2L1Ey1scsRQgghSh0JREIIIYQo9SQQCSGEEKLUk0BUCgUHB9OrVy+DtHX48GE6dOhAmTJlcHFxoVOnToSFhRmkbSGEECWfoijGLgGQQCSeQlJSEp07d6ZChQocPHiQPXv24OTkROfOneXuMCGEEP8pKjGKledW8s7Wdxi8ebCxywFkHiKDUhSFtCyN3u97eCB1QQdV21qaF/mdWRcuXODevXtMmTIFX19fACZNmkRQUBDXrl2jUqVKRVqPEEII05SuTufIrSPsid7Dnug9XE28muv16KRofBx9jFSdjgQiA0rL0hA4cdNTtdHw020Fet/ZKZ2xs3r6j7NLly7s3r37sfskJycDUK1aNdzc3Fi8eDEfffQRGo2GxYsXU7NmTfz8/J66FiGEEMXXtcRr7I7ezZ7oPRyJOUK6Jj3nNQuVBfXK1aOlT0ta+rTE28HbiJXer8nYBQjTsmjRItLS0vK1r6OjI6GhofTs2ZOpU6cCULVqVTZt2oSFhfzVEkKI0iRdnc7hmMM5Z4GuJV3L9Xo5u3K09GlJK59WNPFqgoOVg5EqzZv81jIgW0tzzk7prPf7UjPVOWeGjozvUKAzPbaWhlnewscn/6cs09LSGDRoEC1atODnn39Go9Ewc+ZMunbtyuHDh3OW4BBCCFEyXU28yp7oPeyO3s2RmCNkaDJyXrMws6C+R/2cs0CVy1Q26Ul3JRDlISQkhJCQEDQa/cYDqVSqp75sZWdlYZBLXwWlzyWzlStXcuXKFfbv34+ZmVnONhcXF/744w9eeeWVQq9XCCFE0UlTp+U6CxSVFJXrdU97z5wA1NSrKfaW9kaqVH8SiPIwdOhQhg4dSmJiIs7OzsYup0jpc8ksNTUVMzOzXIk/+7lWqy2sEoUQQhQRRVFyzgLtid7D4ZjDZGozc163MLOggUeDnBBUqUwlkz4L9DgSiEqphISER+YLcnV1pUKFCvluo2PHjnzwwQcMHTqU4cOHo9VqmTFjBhYWFrRr187AFQshhCgK2WeBdl/XDYi+nnw91+te9l45AaiJV5NidRbocSQQlVKhoaHUq1cv17YBAwawdOnSfLdRvXp1NmzYwOTJk2nWrBlmZmbUq1ePf/75By8vLwNXLIQQojAoisKVxCs5Z4GOxBzJdRbI0syS+uXq08qnFS19WlLRuWKxPQv0OBKISqGlS5fqFXwep2PHjnTs2NEgbQkhhCgaqVmpurNA92+Lj06OzvW6t713rrNAdpZ2Rqq06EggMgF2VhZcmdHN2GUIIYQooRRFITIxkj3X758FunWELO2DFQUszSxpUE43FqhV+VYEOAUU3VkgTRaYWxbNsR5DApEQQghRAqVmpXIo5lDOpbB/nwXycfDJOQvU2LNx0Z0FUhSIOQnhmyF8C6gz4M2dRXPsx5BAJIQQQpQAiqIQmRCZcxns6K2jj5wFaliuoS4ElW9ZtGeBMpLg0o4HISg5JvfrcVfA1b9oavkPEoiEEEKIYio1K5UDNw+wJ3oPe6P3ciPlRq7Xs88CtfJpRSPPRkV7FuhOOIRv0oWgq/vhoXCGpR1UbAtVOkLljlDGt2jqegwJREIIIUQxoSgKl+Iv5VwGOxp7FLX2waLgVmZWNPJslHMpzM/Jr+jOAmWlQeTu+2eBNkN87gVcca0IVTrrQpBfC7C0KZq68kkCkRBCCGHCUrJScp0FuplyM9frvo6+OQGokWcjbC2KcNmke1cfBKDIXaB+sIAr5lbg3xKqdNL9V7ZS0dVVABKIhBBCCBOiKAoR8RE5Z4GOxR579CyQV6OceYH8nPyKrjh1JkQd0AWgi5vhzoXcrzv5PAhAAa3B2rQWcH0cCUSmIDMFpnnrHn90A6xKxqyfQggh8ic5M5mDNw+yO3o3e2/sJSYl96Bjo54FSorRDYQO36wbGJ2Z9OA1lTn4NoGq90OQRyAU00kbJRAJIYQQRUxRFMLjw3POAh2/dRy18uAskLW5dc5YoFY+rajglP9llZ6aVgPRR++fBdqku0X+YXZuunFAVTpBpXZg61J0tRUiCUR5KOhq98VFcHAw8fHxrFu37qnbGjlyJHv27OH06dPUqFHjkfXRAE6dOsWwYcM4dOgQrq6uvPnmm0yYMKFETv0uhBD/JTkzOWcs0J7oPdxKvZXrdT8nv5yzQA3LNcTGoggHHafGQcQ2XQiK2Appcblf966vC0BVO4FXPTAzK7raiogEojyU5tXu9aUoCoMGDeLgwYOcPHnykdcTExPp2LEj7dq14/Dhw1y8eJHg4GDs7e157733jFCxEEIUDUVRuHjvYk4ACosNy3UWyMbcJtdZIF+nIrz1/N+TI14/DIr2wevWzlC5vS4EVX4GHDyKrjYjkUAknsrcuXMBuH37dp6B6KeffiI9PZ2lS5dibW1NrVq1uHjxIrNmzWL06NFylkgIUaKotWp2Xd/Fzus72RO9h9jU2Fyv+zv555wFalCuQdGeBXrS5IgeNXWXwqp2hvKNwbx0RYTS1dvCpiiQlar/+zJT836sD0s7gwxk69KlC7t3737sPsnJyflub//+/bRp0wZra+ucbZ07d2bcuHFcuXKFgICAAtcqhBCmIjUrlbURa1l2Zlmu2+JtzG1o7NVYF4K8Wxb9WaD8To5YpRM4ly+62h6SodZgbWFulGM/TAKRIWWlPrhbrKBmVi7Y+wx0d9qiRYtIS0t76nayxcTE4O/vn2tbuXLlcl6TQCSEKM7upN1h5bmV/HLhFxIzEwFwsXahW8VutPJpRQPPBlibWz+hFQPKSoMre3SDofOcHLHS/dvijTc5YpZGy7Gr99gVfptdF++g1ir8PbJVkdfxbxKIRC4+Pj4Gb/Pfl8UURclzuxBCFBeXEy6z/Mxy1l9an7NeWAXHCvQP7E+Pyj1kcsR/iYpLZefF2+y6eJt9l+6SnKH+1+sp+Load8oZCUSGZGmnO1Ojr8zUB2eG3o8AqwKsNWOg9WkMfcnM09OTmJjc16ljY3XX1LPPFAkhRHGgKArHYo+x9PRSQq+H5mwPcg9iYM2BtPNth7lZEVz6KQaTI6Zmqjlw+S67Lt5h58XbRN5JyfW6q70Vraq40aaqOy2ruOHhaPxlPCQQGZJK9fSXrazsjDoxo6EvmTVr1oyPPvqIzMxMrKysANi8eTPe3t6PXEoTQghTpNFq2HZtG0vPLOXUnVMAqFDR1rctA2sNpJ5HvcIvwsQnR1QUhfMxSey6eJudF29z5Mo9MjUP7lqzMFNRv4ILbaq507qKOzW9nTAzM62rBBKISqmEhIRH5gxydXWlQgX9Jv+KiIggOTmZmJgY0tLSctoMDAzEysqKV199lcmTJxMcHMxHH31EeHg406ZNY+LEiXLJTAhh0tLUafwR8QfLzy4nKikK0C2b0aNyD/oH9ifAuRDHQD5pckR7d90q8VU6Gm1yxLiUTHbfHwe0O/w2sUkZuV4v72JLm6rutK7qTvNKZXG0sSzyGvUhgaiUCg0NpV693P+qGTBgAEuXLtWrnSFDhrBz586c59ltRkZG4u/vj7OzM1u2bGHo0KE0bNgQFxcXRo8ezejRo5+6D0IIURjupt1l1YVVrDq/iviMeACcrZ15udrL9K3eFzdbt8I5cH4mR6x6f7V4I0yOqNZoOR4Vz677Y4FORidwf0goALaW5jSrVJbWVdxoXdWdADf7YvUPXwlEpdDSpUv1Dj7/JTQ09In71K5dm127dhnkeEIIUViuJl5l2ZllrL+0ngyN7myHj4MP/QP706tyL+wMNFYzR74nR+wMlTsYZXLE6/dS2XXxDrsu3mZvxB2S/jUYurqnY85ZoIb+LiZx+3xBSSASQghRqoXFhrH0zFK2X9uOgu6UR62ytQiuFUyHCh2wMDPgr8r8TI6YPRbICJMjpmVqOBB5N2cs0OXbuQdDu9hZ0rKKuy4EVXHDw8n4g6ENRQKRKbCyh08SjF2FEEKUGhqthtCoUJaeWUrY7bCc7W3Kt2FAzQE0LNfQMJd7THxyREVRuHgrmZ0XY9l18Q6HrsSRqX5wlsrcTEU93zI5Z4Fq+ThjbmKDoQ1FApEQQohSI12dzvpL61l+djlXE3WTFlqaWfJcxecYUHMAlcoYYJ4eE58cMT41k93hustgu8PvEJOYnut1nzK2tK7qTpuqbjSr5IazrWkPhjYUCURCCCFKvHvp93IGSsel6wYrO1o58nK1l3m1+qu427k//UFunYED38Lx5bm3G3lyRLVGy4nr8ey8Pxbo5PV4tA8NhraxNKNpxbK0rqI7C1TJvXgNhjYUCURCCCFKrKjEKJafXc66iHWka3RnQrztvekX2I/nqzyPveVTzvumzoTzG+DQIri2L/dr9fpBta5GmRzxRnya7m6w8NvsCb9DYnruwdDVyjnSuqrubrBG/q7YWBbfwdCGIoFICCFEiXPy9kmWnlnKtmvb0N6/c6uGaw2CawbTyb/T0w+UToiGo0vh2DJIvqXbpjKHGs9BoyHg36pIJ0dMz9JwMDKOnRd0ISgiNveKAs62lrSs4kabKu60quqGl3MRLi1STEggEkIIUSJoFS27ru9iyeklHIs9lrO9hU8LBtYcSGPPxk93KUhRIHInHF4E5zeCotFtd/CEBsHQYAA4PeUC3/kuRSEiNpmd9+8GOxQZR8ZDg6HNVFDXtwxtqnrQuqobQeXLlNjB0IYigSgPISEhhISEoNFojF2KEEKIJ8jQZPDnpT9ZdnYZkQmRAFiYWdA1oCsDag6gqkvVpztAegKE/awLQnfDH2z3awmNh0D158C88AceJ6RmsSfiTs6lsJsJuQdDeznb5NwN1qKSG852pWMwtKFIIMrD0KFDGTp0KImJiTg7Oxf68VKzUmmysgkAB189aPjJv4QQogRKyEjg1wu/8tO5n7ibfhcAB0sHXqz2Iq9Vf41y9k+5gHTMKV0IOvkrZKXqtlk5QJ1XdJfFPGo8ZQ8eT6NVOHH9wczQYVG5B0NbW5jRpKJuZug2Vd2p7OFQKgdDG4oEIiGEEMXK9aTrrDi3gjXha0hT6xajLmdXjn6B/ehTpQ8OVk8xgFmdCefWw6GFuhXls7nXgEaDdWHI2vEpe/DfYhLSdZMi3h8MnZCWlev1Kh4OtL5/FqhJgAyGNiQJRKVQcHAw8fHxrFu37qnbGjlyJHv27OH06dPUqFHjkQVjQ0ND+frrrzl06BCJiYlUqVKFDz74gNdee+2pjy2EKF3O3D3D0tNL2Xx1c85A6aouVQmuGcyzAc9iafYUl4jiox4Mkk65rdtmZgE1uuvOBvm1KJRB0ulZGg5ficuZGfrirdyDoZ1sLGhZxS3nlnjvMjIYurBIIBJPRVEUBg0axMGDBzl58uQjr+/bt4+goCA+/PBDypUrx19//UX//v1xcnKie/fuRqhYCFGcKIrC7ujdLDuzjEMxh3K2N/NqRnDNYJp5Nyv4ZSKtFiJDdbfMX/z7wTpijl7QYCDU7w9OXk/fiYcoisKl2ynsvH8Z7GDkXdKzHgyGVqmgTvky9ydGdKdOeWcszIt2EdfSSgKReCpz584F4Pbt23kGoo8++ijX8xEjRrBp0ybWrl0rgUgI8Z8yNZlsjNzIsjPLiIiPAMBcZc6zAc8SXDOY6q7VC954WjyErYQji+FuxIPt/q2g8Ru6uYMMOEg6IS2LfRF32BV+m10X7xAdn5br9XJO1rSu4k6bau60rOxGGTsrgx1b5J8EIgNSFCXnerY+Hn5PQd4PYGtha5DBdF26dGH37t2P3Sc5Ofmxrz9JQkICNWoU7mBEIUTxlJiZyG8XfuOncz9xO0136crOwo4Xqr7A6zVex8vhKc7Y3DwJhxfCyd8g+2etlSPU7au7LOZezQA90A2GPhWdkDMY+nhUPJqHRkNbWZjRJMA15zJY1XIyGNoUSCAyoDR1Ws7dYgXV9te2BXqfoe5OW7RoEWlpBQtl+fH7779z+PBhvvvuu0I7hhCi+LmZfJMfz/3I6ourSVXr7ujysPXgtcDXeKHqCzhZORWsYXUGnFmnu1vs+oNLbnjU1A2SDnrZILNIa7UKh6/EsfZ4NKsORz3yeiV3+5zB0E0DymJrJYOhTY0EIpGLj49PobUdGhpKcHAwCxcupGbNmoV2HCFE8XE+7jxLzyzln8h/0Nyf6LBymcoE1wyma0BXLAt66Sr+Ghz5AY79CKl3dNvMLCCwp+5sUIVmTz1IWlF0Z4LWh93gz5M3H1kktWOgB+2qlaN1VTfKu8h0KqZOApEB2VrYcvDVg3q/L02dlnNmKPSlUGwt9L+LoCDvyUthXTLbuXMn3bt3Z9asWfTv37+g5QkhSgBFUdh/Yz9LzizhwM0Ht7Y39mxMcM1gWvq0LNglJK0WLm/XDZIO3/RgkLSTz4NB0o5POTcREH4rifUnbrDhxA2u3E3N2e5oY8GzNT3pXsebZpXKYimDoYsVCUQGpFKpnvqyla2FrVEnZiyMS2ahoaE899xzfP755/zvf/8zaNtCiOIjS5vFP5H/sPTMUi7euwjoBkp38uvEgFoDqFm2gGeOU+MeDJKOu/xge0Ab3SDpql3A/Ol+3V27m8qGk7oQdD4mKWe7jaUZz9QoR4863rSp5o61hVwKK64kEJVSCQkJj8wZ5OrqSoUKFfRqJyIiguTkZGJiYkhLS8tpMzAwECsrK0JDQ+nWrRsjR46kT58+xMTEAGBlZYWrq6shuiKEMHHJmcn8fvF3Vpxbwa1U3UKotha29KnSh9cDX8fHoYCX6m8c140NOvU7qO9frrJ2grqvQsPB4P50S3bcSkznr5M3WX/iBmFR8TnbLc1VtKnqTvc63jxToxz21vKrtCSQT7GUCg0NpV69erm2DRgwgKVLl+rVzpAhQ9i5c2fO8+w2IyMj8ff3Z+nSpaSmpjJ9+nSmT5+es1+bNm0IDQ0tcP1CCNMXkxLDynMr+e3ibyRn6S61l7Upy2s1XuOlai/hbF2ApZGy0uHMWl0Qij7yYHu5WrqxQUEvgZV9gWu+l5LJ36dj2HDiBgci76LcvznMTAXNKpWlRx1vnq3pJeuElUAqRVGUJ+9WOmWvZZaQkICTU+47HNLT04mMjCQgIAAbG5unOo6sZfYoQ/75CiGK1oW4Cyw/u5yNlzeiVtQABDgHEFwzmG4Vu2Ftbq1/o/euPBgknRan22ZmCTV7QaM3wLdxgQdJJ2eo2XI2hvVhN9gdfgf1Q7fIN/BzoXuQF12DvPBwlJ9FxcXjfn//FzlDJIQQ4qkpisLBmIMsPb2UvTf25mxvUK4BA2sOpFX5Vpip9BxkrNXCpW26dcXCNwP3g4pTeWh4f5C0g0eB6k3P0hB6IZb1J26w7VwsGeoHs0UHejnRo6433Wp74esq/0AtLSQQmQA7SztODThl7DKEEEJvWdosNl/ZzLIzyzgXdw4AM5UZHSp0YGDNgdR2r61/o6lxcHyFbpD0vSsPtldqr7ssVqVzgQZJZ2m07Im4w4YTN9h85hbJGeqc1yq62dO9jjfd63hT2ePp5yUSxY8EIiGEEHpLyUphTfgafjz7IzdTbgJgY25Dr8q96B/YH18nX/0bjT4KhxfD6dUPDZJ2hnqv6QZJu1XWu0mtVuHQlTjWn7jB36duci/1werx3s42OSGopreTzBZdykkgEkIIkW+xqbGsPLeSXy/+SlKm7vZzVxtX+lbvyyvVXqGMTRn9GsxKg9NrdIOkbxx7sN0zSHfLfK0+eg+SVhSFk9cTWH/iBn+evMGtxIyc19wcrOha24sedbypX8EFMzMJQUJHAtFTkjHphUP+XIUwLZfiL7H0zFL+vPwnaq3uUpOfkx/9A/vTo1IPbCz0HHAcF6m7JHZ8BaTd020zt4Kaz+sGSZdvqPcg6Yu3klgfdoMNJ29wNY8JE3vU9aZZxbKyerzIkwSiArK01N1ymZqaiq2tYWaJFg+kpup+mGX/OQshjCM5M5mZR2ayOnx1zra67nUJrhVMO992+g2U1mogYqtukHTEVnIGSTtX0A2SrtcPHNz1qi97wsT1YTe4cOvBhIm2luY8E1iO7kFeMmGiyBcJRAVkbm5OmTJliI2NBcDOzk6uPxuAoiikpqYSGxtLmTJlMDeXH2JCGMuu67uYemAqMSkxOdsWdlxIU++m+jWUcheO/6i7bT7+6oPtlZ+5P0i6E5jl/7t+KzGdDSdusOHkTU48MmGiBz3qevNMDQ/srORXnMg/+dvyFDw9PQFyQpEwnDJlyuT8+Qohita99Ht8cfgL/rz8JwA+Dj580vwTmnrpEYQU5f4g6UW6MUKa++N4bMpAvdeh4SAoWyn/NaVksvH0TTacuMHByLhcEyY2r+RGjzredK7pKRMmigKTQPQUVCoVXl5eeHh4kJWV9eQ3iHyxtLSUM0NCGIGiKGy6uonpB6cTlx6HChWvB77OsLrD8j9hbGaq7i6xw4vgZtiD7V51dYOka/YGq/y1lZSexZazt1h/4gZ78pgwsUcdb7rU9pQJE4VBSCAyAHNzc/kFLoQo1mJTY/n0wKfsiNoBQCXnSkxuMZk67nXy18DdS7pLYsdXQHq8bpu5NdTqrRsk7VM/X4Ok07M07DivmzBx+/ncEybW9Haiex1vngvyoryLTJgoDEsCkRBClGKKorA2Yi0zD88kKSsJC5UFbwS9wZDaQ7Ayt3r8m7Ua3QzShxbqZpTOVqaCbt6gev3AvuwTa8jSaNkTfn/CxLP/mjDR3Z4edbx5LkgmTBSFSwKREEKUUlFJUUzeN5mDMQcBqFm2JlNaTKGqyxNWiU+5A8eWw5ElkHDt/kYVVOmoGyRd+ZknDpLWaBUORd6fMPH0TeIfmjDRp4wtz9XxonuQTJgoio4EojyEhIQQEhKCRqMxdilCCGFwGq2Gn879xLzj80jXpGNtbs3wesN5rcZrWJj9x68FRYHrh3Vng86uA02mbruti+5MUMOB4FrxscdVFIUT1xNYH3aDv049OmFit9pedJcJE4WRyGr3j1GQ1XKFEMKURdyLYNK+SZy8cxKARp6N+KTZJ1RwqpD3GxQFzv8Fu76AmycebPeurzsbVKs3WD5+LrYLMUmsPxHNhhM3uRaXe8LELrU86VHHh6YVXWXCRGEwstq9EEKIPGVpslh0ehHfn/wetVaNg6UDoxuOpk+VPv89ueLlUNg2RXf7fLagl6DJW+DT4LHHu3o3hQ0nbrD+xA0u3krO2Z49YWKPOt60ruomEyYKkyGBSAghSrjTd04zcd9Ewu+FA9CmfBvGNx2Pp/1/zPV1/ShsmwyRO3XPLe2h6dvQbCjYuf7ncWIS0vnz5A02nLjBiesJOdtlwkRRHMjfSiGEKKHS1GnMD5vP8rPL0SpaXKxdGNt4LF0CuuQ9UDn2HGz/FM7rJmTE3Eo3gWKr98DBI89jxKVksvGUbsLEQ1dyT5jYorIb3YNkwkRRPEggEkKIEuhwzGEm7ZtEVFIUAF0DuvJh4w9xtcnjDM+9K7BjOpz8BVBAZQZ1+kLbsbpb6P8lKT2LzWfuT5gYcQfNQxMmNvRzoUddb7rU8sLd0bqQeieE4UkgEkKIEiQpM4lZR2fx+8XfAfCw82Bi04m08W2Tx84xsGsmHF0K2vu3vdfoAe3Hg3u1XLumZ2nYfj6W9WE32H4hlsx/TZjYo443z9XxxqeMLHYtiicJREIIUULsjNrJlANTiE3Vra/4YtUXGdVgFI5Wjrl3TLsHe+fAgW9BnabbVqk9tJ+gm1H6IZF3Uli0+zI/HbyWa3v2hInd63hTyV0mTBTFnwQiIYQo5uLS45hxaAZ/R/4NgK+jL5ObT6aRZ6PcO2amwIEFsHcuZNwf9Fy+EXSYBAGtcnZTFIU9EXdYsvcKOy7E8vDkLINb+tO7fnkCvWTCRFGySCASQohiSlEU/o78mxmHZnAv4x5mKjP6B/bnnbrvYGvx0KUrdQYcXQa7voQU3dkjPGpChwlQ9dmcNcbSMjWsPR7N0n2RuW6Vb1fNnYEtAmhZ2U0mTBQllgQiIYQohmJSYvj0wKfsvK67Nb6KSxWmNJ9CLbdaD3bSanQDpUOnQ/z9S14u/tDuY6jVJ2d5jZsJaSzff5WfD13LWULDzsqcFxuUZ0BzfyrKJTFRCkggEkKIYkSraPn94u/MOjqLlKwULMwseDPoTQbXGoyl+f1b2xUFzm3Q3UJ/54Jum4MntBkD9fuDuSWKonDs6j1+2BvJP6djcu4UK+9iS3Bzf15s6IuzrdwqL0oPCURCCFFMXEu8xif7P+FwzGEAgtyCmNx8MpVdKj/Y6dIO3ezSN47pntuUgVajodEbYGVHplrLxuPRLNkbmWvyxCYBrgxqGcAzNcphLpfFRCkkgUgIIUycWqtmxdkVfBP2DRmaDGwtbBlebzivVn8V8+xV5aMOw/YpELlL99zSHpq9A82Hg40zd5MzWLk7nB8PXCU2SbeoqpW5GT3rehPcwp+a3s5G6p0QpkECkRBCmLCL9y4yae8kTt89DUATryZMajYJX0df3Q63zuoujV34S/fc3AoaDr4/u7Q7Z28ksmTvCf44cSNn7iB3R2v6NfXj1SYVcHOQyROFAAlEQghhkjI1mSw8tZBFJxehVtQ4WjryfqP3eb7y87rb3eMidYOlT/5KzuzSdV+FNh+icfJl67lbLNm7nwOX43LaDCrvzMAW/nSr7Y2VhawsL8TDJBAJIYSJOXn7JJP2TSIiPgKAdr7tGN90PB52HrrZpXd+AceWgVate0NgT2g3nkTHAH49HMWy/TuIitNNuGhupuLZWp4MauFP/QouMneQEP9BApEQQpiI1KxUvgn7hhVnV6Cg4GrjykdNPqKTXydUafdgyyQ4+N1Ds0t3gA4TiLSqyrJ9V/jtyDZSMjUAONta0rdxBfo388NbltMQ4okkEAkhhAk4cPMAn+z7hOjkaAC6V+zOmEZjKKOygN0zYe+8h2aXbozSYSJ71TX4YXMkOy6E5swmXcXDgYEtAni+ng+2VuZG6o0QxY8EIiGEMKLEzERmHZnF6vDVAHjaezKx6URaeTbWLbq660tIua3buVwtMlp/zOqkmixdd4WLtw7mtNO+ugcDW/jTsrKbXBYTogAkEAkhhJFsv7adTw98yu00XeB5udrLvFt3OA7n/oTf3oGE7NmlA7jXdAzfx9Xl59XXiU/V3XEms0kLYTgSiIQQoojdTbvLjEMz+OfKPwD4OfkxudknNIiLhkUdc2aXVhy9uFZrGDPvNGbjH3fQaCMBmU1aiMIggUgIIYqIoij8eflPPj/8OQkZCZirzAmuGcxbjtWxWf8+3Diu28/WhTMVBzM5pjmHd6QDujNIMpu0EIVHApEQQhSBm8k3mXJgCnui9wBQzaUakyu9TM3Dy+DKJAAUSzsOe77KhzdbE3nUAkjHysKMnnVkNmkhCpsEIiGEKERaRctvF35j1tFZpKpTsTSz5O1KvQm+chrL34bo9jGzYq9LTz649Qwx4Y6AzCYtRFGTQCSEEIXkSsIVJu2bxLFY3UKrdV2qMznDmoqbZwIKCmbssOvI+Lhu3Eh1A3SzSQ9qEUDX2l4ym7QQRUgCkRBCGJhaq2b52eXMD5uvW4zV3IaR1hXoG7Yds/uzS283b86nqb25nO6NuZmKbjKbtBBGJYFICCEM6ELcBSbsncC5uHMANLP2YFLkOXwyLgKwW6nLjMwXOaME4GxryVsym7QQJkECkRBCGECmJpPvTn7HD6d+0C3GambFmLvx9Iw/ggo4oq3KF1kvc0ipQRUPB6bJbNJCmBQJREII8ZTCYsOYuG8ikQm6eYKeSVfzcex13DRazmr9+FL9Eju0dWlfvRw/ymzSQpgkCURCCFFAqVmpzD0+l5XnVqKgUFar8PHtO3RMTSNSW47J6hfZbtGCF5pUYLvMJi2ESZNAJIQQBbDvxj6m7J9MdPINAHomJfNBXDypmjKMU7/KPsdn6deiMp818sXJRmaTFsLUSSASQgg9JGQkMPPwTNZdWgeAd5aaSXfjqJ5qwTx1X877vkS/VtX5VGaTFqJYkUAkhBD5tO3qNqbuncjdrERUisKrickMikvnF203NgYOoW/rmjKbtBDFlAQiIYR4gjtpd/h48/vsiz8KQEBmFh/fTiBC3Z4/moykb6s6Mpu0EMWcBCIhhPgPao2Wr7fN5vfoZaSaaTFXFAbGJ1EprRHxbcfyYuP6Mpu0ECWEBCIhhPiXxPQsftqxmc1XpxJhmwJmUCMjk1dSK1K1/VRq1m4gt80LUcKUikD0/PPPExoaSocOHfj999+NXY4QwkRF3knhl53HuXHlUw64RpFma4aVVuGVdGdebz8Tr2rNjF2iEKKQlIpANGLECAYNGsSyZcuMXYoQwsQoisLeiLus3H0G96jviPA8xkk3S8CM2moLJjX6mGpBLxi7TCFEISsVgahdu3aEhoYauwwhhAlJz9Kw9ng0K/ZcoNHdtVRz28zyCjZkmllip8C7FZ/n5ZaTMDOTpTWEKA2MPhpw165ddO/eHW9vb1QqFevWrXtkn/nz5xMQEICNjQ0NGjRg9+7dRV+oEKJESM/S8P2uSwRO2MixP+YxJu1tzvhtY5GrLZlmKlo4VWZdn3/o23qKhCEhShGjnyFKSUmhTp06DBw4kD59+jzy+i+//MK7777L/PnzadGiBd999x1dunTh7NmzVKhQAYAGDRqQkZHxyHs3b96Mt7d3ofdBCGH60rM0/HI4ivmhEVgmXWeF9Tccdo1lpLMTGpUKZ3MbPmzyEc9V7iUDpoUohYweiLp06UKXLl3+8/VZs2YxePBghgwZAsDs2bPZtGkTCxYsYPr06QAcPXrUILVkZGTkClaJiYkGaVcIYTwZ6vtBaMclYhLT6WG2lxF2yxhfzp4z1rpJFDtXeIaxTT/GzdbNyNUKIYzF6IHocTIzMzl69Chjx47Ntb1Tp07s27fP4MebPn06kydPNni7Qoiil6HW8OvhKOaHXuJmQjoOpLLA7kfKWh1gkIcb98zNMVMUZrT+ki4V//sfZUKI0sGkA9GdO3fQaDSUK1cu1/Zy5coRExOT73Y6d+7MsWPHSElJoXz58qxdu5ZGjRo9st+4ceMYPXp0zvPExER8fX0L3gEhRJHLUGv49ch15u+I4GZCOgDPOFxlluU3rLVMZJyrBxqVihou1fm6/Wx8HHyMXLEQwhSYdCDK9u/r+Yqi6HWNf9OmTfnaz9raGmtrmX5fiOIoU63l1yNRzN8RwY37QcjL0YL5fqFUu/w9nzg487eDCwDdK3ZnYrOJ2FjYGLNkIYQJMelA5Obmhrm5+SNng2JjYx85aySEKJ0y1Vp+P3qdkB0RRMenAVDOyZoPmtjR+8pkoq8cpp+nOxetrbBQmfN+ow94tfqrMnBaCJGLSQciKysrGjRowJYtW3j++edztm/ZsoWePXsasTIhhLFlqrWsPnadb7Y/CEIejta83bYSr9kfxurvQew1y2CMjxeJZma42rjyVZuvaOjZ0MiVCyFM0VMFovT0dGxsnu6Uc3JyMhERETnPIyMjCQsLw9XVlQoVKjB69Gj69etHw4YNadasGd9//z3Xrl3jrbfeeqrjCiGKpyyNltVHr/PNjgiu39MFIXdHa95uU4lX67pgs2UsyomfWeTsxFwXDxQVBLkFMavtLMrZy5llIUTe9A5EWq2Wzz77jG+//ZZbt25x8eJFKlasyIQJE/D392fw4MF6tXfkyBHatWuX8zx7UPOAAQNYunQpL7/8Mnfv3mXKlCncvHmTWrVqsXHjRvz8/PQtPd9CQkIICQlBo9EU2jGEEPrJ0mhZc+w687Y/CEJuDvfPCDWpgM2t47C4FynxVxnv4cZWezsA+lTpw0dNPsLK3MqY5QshTJxKURRFnzdMmTKFZcuWMWXKFN544w1Onz5NxYoV+fXXX/n666/Zv39/YdVa5BITE3F2diYhIQEnJydjlyNEqZSl0bL2WDTzdoQTFfcgCL3VpiKvNfHD1gLYPQtCpxNpoeJdLy8um6uwNLPkoyYf8UJVWYdMiNKmIL+/9T5DtHz5cr7//ns6dOiQ67JVUFAQ58+f17c5IYTIU5ZGy9rj0XyzPYJrcakAuDlY8VabSrogZGUO8VGw5n9wbR877Gz5qJwnyWjwsPVgVrtZ1HGvY+ReCCGKC70DUXR0NJUrV35ku1arJSsryyBFCSFKL3V2ENoRwdW7D4LQm60r8XrT+0EI4PRq2DAKbUYCC8q6862TLaChvkd9vmr7lcw6LYTQi96BqGbNmuzevfuRMTy//fYb9erVM1hhQojSRa3Rsi7sBvO2h+cEobL2VrzZpiKvN/XDzur+j6uMJNg4Bk6sJNFMxTi/yuwyywTg1eqv8n6j97E0szRWN4QQxZTegWjSpEn069eP6OhotFota9as4cKFCyxfvpw///yzMGoUQpRgao2WP+4HoSv3g5CrvRVvtq5Iv2YPBSGA60dh9WC4F0mElRUj/SpzTZ2Mtbk1E5tNpEelHkbqhRCiuNN7UDXoZn6eNm0aR48eRavVUr9+fSZOnEinTp0Ko0ajkUHVQhQetUbL+hM3mLc9gsg7KYAuCP2vdUX6NfXD3vqhIKTVwJ5ZsGM6KBo2ufkywdmGNG0GXvZezG43m8CygUbqiRDC1BTk93eBAlFpIYFICMPTaBXWn4hm3rYILt8PQi52lvyvdSX6N/tXEAJIuK4bOH11LxpgTpXGLFHrZq9v4tWEL1t/iYuNSxH3QghhyorkLrOKFSty+PBhypYtm2t7fHw89evX5/Lly/o2aXJkHiIhDE+jVdhw4gZzt4XnCkJvtK7IgGb+jwYhgDNrYcNISE8g3saRD6rU40DyFQAG1hzIiPojsDAz6Qn3hRDFhN5niMzMzIiJicHDwyPX9lu3blGhQgUyMjIMWqAxyRkiIZ6eRqvw58kbzNkWzuXbuiBUxs6SN1pVZEBzfxzyCkIZyfD3hxC2AoBz5YMY5WRFdFostha2TGk+hWcDni3KbgghipFCPUO0fv36nMebNm3C2dk557lGo2Hbtm34+/vnv1ohRImWHYTmbgvnUn6DEED0UVg9BOIuAyo21H+ByYlhZKTF4+voy+x2s6nqUrXoOiKEKBXyHYh69eoFgEqlYsCAAbles7S0xN/fn6+++sqgxQkhih+NVuGvUzeZuy2ciNhkAJxtLXmjVQADmvvjaPMft8RrNbB3NuyYBlo1WU7l+SqoAz9F7wCgpU9LZrSagbO1c97vF0KIp5DvQKTVagEICAjg8OHDuLnJpGdCiAe0DwWh8PtByMnGgjdaVSS4xWOCENwfOP0mXN0DwJ3A53jfQcXR+2Hof0H/450672BuZl7o/RBClE56j0aMjIwsjDqEEMWUVquw8bQuCF289SAIDbkfhJweF4QAzqy7P3A6HiztOdX2Xd69sYnY27HYW9rzWcvP6FChQ6H3QwhRuhXo9oyUlBR27tzJtWvXyMzMzPXaiBEjDFKYEMK0abUKf5+OYc62izlByNHGgiEtdUHI2fYJQSgjGf75EI7rBk7jXZ81jV7m09Pfk6XNwt/Jnznt51DRuWIh90QIIQoQiI4fP07Xrl1JTU0lJSUFV1dX7ty5g52dHR4eHhKIhCjhtFqFf87EMGdrOBduJQG6IDS4ZQADWwQ8OQgBRB+7P3D6EqAis8UIZtib8dvJEADa+7bns5af4WDlUIg9EUKIB/QORKNGjaJ79+4sWLCAMmXKcODAASwtLXn99dcZOXJkYdRY5GQeIiEepdUqbDoTw5xt4ZyPuR+ErC0Y1DKAQS3zGYS0Gtg7B3Z8Blo1OPkQ2/VzRl1excnwk6hQMazeMIbUHoKZyqyQeySEEA/oPQ9RmTJlOHjwINWqVaNMmTLs37+fGjVqcPDgQQYMGMD58+cLq9YiJ/MQCaELQpvPxjB7a+4gNLBlAINbBOBsl8+FVBOiYe2bcGW37nlgT441HcTo/Z9wN/0ujlaOfN7qc1qVb1VIPRFClBZFMlO1paUlKpUKgHLlynHt2jVq1KiBs7Mz165d07c5IYSJUhSFzWdvMXtrOOduJgLgYG3BoBb+DGoZQBk7q/w3dvYPWD8iZ+C08uwMVtlZ8kXoKNSKmiouVZjTdg6+Tr6F0xkhhHgCvQNRvXr1OHLkCFWrVqVdu3ZMnDiRO3fu8OOPP1K7du3CqFEIUYQURWHL/SB09qEgNLCFP4P1DUKZKfDPWDi2XPfcux7pveYz9eJPrD+tm+z1Wf9nmdx8MnaWdobuihBC5JvegWjatGkkJelOm0+dOpUBAwbw9ttvU7lyZZYsWWLwAoUQRUNRFLaei2X21oucuaELQvZW5gxsEcDglgG42OsRhABuHNcNnL4bAaig5bvcaDiQUXvGcPbuWcxUZoxuMJr+gf1zzjoLIYSx6BWIFEXB3d2dmjVrAuDu7s7GjRsLpTAhRNFQFIVt52KZve0ip6MfBKHgFv4MaVlR/yCk1cK+ubD9U9BmgaM39P6OgzbWfPBPP+5l3KOMdRm+bPMlTb2aFkKPhBBCf3oHoipVqnDmzBmqVKlSWDUJIYqAoihsPx/L7K3hnIpOAMDOypzg5v4MaVURV32DEEDiDd3A6chduuc1eqA8N5vlV/5k1u5ZaBUtNVxrMLvdbLwdvA3YGyGEeDp6BSIzMzOqVKnC3bt3JRAJUUwpisKOC7ogdPL6gyA0oLk/bxQ0CAGc2wDrh0PaPbC0gy6fk1qrD5/sn8zfV/4GoEelHkxoOgEbCxtDdUcIIQxC7zFEX3zxBR988AELFiygVq1ahVGTEKIQKIpC6IXbzN56kRMPBaH+zfx5o1UAZR2sC9ZwZgr8Mw6OLdM996oLfRYTZWXNyH/6EX4vHAuVBR80+oC+1fvKeCEhhEnSex4iFxcXUlNTUavVWFlZYWtrm+v1uLg4gxZoTDIPkSgJFEUh9OJtZm8N50RUPAC2lub0b+7H/1pVLHgQArgRdn/gdDigghYjoN149tw6xJhdY0jKTKKsTVm+avsVDco1MER3hBDiiYpkHqLZs2fr+5ZiR2aqFiWBoijsvB+Ewh4OQs38eKN1RdyeJghptbB/Hmyben/gtBc8/x1KQGsWnVrEvOPzUFAIcg9iVptZlLMvZ5hOCSFEIdH7DFFpImeIRHGkKAq7wu8we+tFjl+LB8DG0ox+Tf14s02lpwtCcH/g9FsQuVP3vPpz0GMeyRZWjN87nm3XtgHwYtUXGdt4LFbmBRyTJIQQBVQkZ4iEEKYpNVNN4MRNubbZWJrxehNdEHJ3fMogBHDuT1g/7MHA6WenQ/0BRCZeYeTmkUQmRGJpZsnHTT6mT9U+T388IYQoIhKIhCgBktKz+PSvcznPrS3MeL2pH2+2qYiHowHu6MpMgU0fwdGluudedaDPYnCrwvZr2/loz0ekZKXgYefB122/Jsg96OmPKYQQRUgCkRDF3LZztxi/7jQ3E9Jztm0e1Rq/svaGOcDNE/D74EcGTmvNLZh//Bu+O/kdAPU96vNV269ws3UzzHGFEKIISSASopi6k5zB5A1n2XDiBgB+Ze2Y/nxtmlc2UCDRauFACGydnGvgNBXbkJCRwLjQceyO1q1c/1qN13iv4XtYmlka5thCCFHE9ApEarUaGxsbwsLCZA4iIYxEURTWHo9myp9niU/NwkwFb7SqyLvPVMXWytwwB0m8Cevehss7dM/vD5zGzpXwe+GM3DGSqKQorM2tmdRsEt0rdTfMcYUQwkj0CkQWFhb4+fnJ7ehCGMn1e6l8tPY0uy7eBqCGlxNf9Amidnlnwx3k/Eb4YyikxYGFrW7gdINgUKn458o/TNw7kTR1Gt723nzd7msCywYa7thCCGEkel8yGz9+POPGjWPFihW4uroWRk1CiH/RaBWW77/Cl5sukJqpwcrCjJEdqvC/1hWxNDczzEEyU2Hzx3DkB91zzyDdwGn3qqi1auYencuSM0sAaOrVlC9af4GLjYthji2EEEamdyCaO3cuEREReHt74+fnh7197oGbx44dM1hxQggIv5XEmNUnc+YUauzvyvQ+tank7mC4g9w8qZtx+s4F3fPmw6H9BLCw5l76PT7Y9QEHbx4EYGCtgYyoNwILMxmCKIQoOfT+idarV69CKEMI8W+Zai3zQyMI2RFBlkbBwdqCsV2q82rjCpiZGWg9MK0WDsyHbZNBkwkOnvD8t1CpHQBn755l1I5R3Ei5ga2FLVNaTOFZ/2cNc2whhDAhMlN1Hh5euuPixYsyU7Uocseu3WPs6pNcvJUMQIfqHnz6fC28nG2f8E49JMXoBk5f2q57Xq2bbuC0fVkANlzawOT9k8nQZODr6MucdnOo4lLFcMcXQohCUpCZqgsciI4ePcq5c+dQqVQEBgZSr169gjRj0mTpDlHUUjLUzNx8gaX7rqAoUNbeikk9atI9yMuwq8Rf+Fs3cDr17v2B09OgwUBQqcjSZjHz8ExWnl8JQCufVsxoPQMnK/kOCCGKhyJZuiM2NpZXXnmF0NBQypQpg6IoJCQk0K5dO1atWoW7u7vehQshYNfF24xbc4ro+DQAetfzYcJzgbjYG3AtsMxU2DwejizWPfesfX/gdDUA7qTd4b3Q9zgWqxsL+Fadt3i7ztuYqQw0cFsIIUyU3j/lhg8fTmJiImfOnCEuLo579+5x+vRpEhMTGTFiRGHUKESJdi8lk9G/htH/h0NEx6fhU8aWpQMbMevluoYNQzGn4Pu2D8JQs2EwZFtOGDpx+wQvb3iZY7HHsLe0Z067OQytO1TCkBCiVND7kpmzszNbt26lUaNGubYfOnSITp06ER8fb8j6jEoumYnCpCgKf568yeQNZ7iTnIlKBQOa+fNB52rYWxvwDi6tFg4ugK2fPDRwegFUap+zy+8Xf2fawWlkabMIcA5gTrs5BDgHGK4GIYQoQkVyyUyr1WJp+ej0/JaWlmi1Wn2bE6JUupmQxoR1p9l6LhaAKh4OzOgTRAM/A8/rk3Tr/sDpbbrn1bpCj29yBk5najKZfmg6v1/8HYAOFTrwaYtPcbAy4C39QghRDOgdiNq3b8/IkSP5+eef8fb2BiA6OppRo0bRoUMHgxcoREmi1Sr8fPgaMzaeJylDjaW5infaVuaddpWwtjDQshvZLvxzf+D0Hd3A6c6fQcNBcH9w9q2UW4zeOZqTt0+iQsXwesMZXHuwXCITQpRKegeib775hp49e+Lv74+vry8qlYpr165Ru3ZtVqxYURg1ClEiXL6dzNg1pzgUGQdAXd8yfPFCEFXLORr2QFlpsHkCHF6oe16uNrzwYOA0wNFbR3kv9D3upt/F0cqRL1p/QUufloatQwghihG9A5Gvry/Hjh1jy5YtnD9/HkVRCAwM5JlnnimM+oQo9rI0WhbuvszsreFkqrXYWprzQedqDGjuj7mhJlgEyEyBad65tzUdCs9MAgtrQDduaeX5lcw8PBO1oqaKSxXmtJ2Dr5Ov4eoQQohiqMCr3Xfs2JGOHTsWVl1ClAinrifw4eqTnL2ZCECrKm5Me742vq52hj9YxNYHj+09dDNOV35wGTtdnc7UA1NZf2k9AF38u/BJ80+wsyyEWoQQopiR1e6FKARpmRpmb73Ioj2RaLQKZewsmdAtkN71fQw7wSKAosDBb2HTRw+2DdkGLhVynt5IvsG7O97lXNw5zFRmjG4wmv6B/Q1fixBCFFOy2r0QBrbv0h3GrTnF1bupADwX5MWk7jVxd7Q2/ME0WfD3hw/mFqrfH7rNAvMHd4IeuHmAD3Z+QHxGPC7WLnzZ5kuaeDUxfC1CCFGMyWr3QhhIQloWM/4+x8+HogDwdLJhaq9adAwsVzgHTIuH34Lh8g5ABZ2m6iZbvH/WR1EUlp1ZxtfHvkaraAksG8jstrPxcvAqnHqEEKIYk9XuhTCATWdimLDuNLFJGQC81qQCH3apjpPNo3N2GURcJKx8Ge5cAEs76LMIqnfLeTk1K5VJ+ybxz5V/AOhZqSfjm47HxsKmcOoRQohiTu9B1QCDBg3C17fk3pXy8Gr3QjxObFI6n6w/w8ZTMQBUdLNneu/aNKlYtvAOenU//PKabmFWR294dRV41cl5+VriNUbuGElEfAQWKgs+bPwhL1d7WcYLCSHEY+i9dIejoyOnTp3C39+/kEoyHbJ0h/gviqLw29HrfPbXORLSsjA3U/Fm64qM6FAFG0sDT7D4sBO/wPphuiU4vOpC31XgpLsElpqVSpOVD8YGudm68VWbr6hfrn7h1SOEECaoSJbu6NChA6GhoQQHB+v7ViFKhGt3U/lo7Sn2RNwBoJaPE5/3CaKmt3PhHVSrhdBpsOtL3fMa3eH578DqwRi+3dd35zyu7Vab2e1m42HnUXg1CSFECaJ3IOrSpQvjxo3j9OnTNGjQ4JFB1T169DBYcUKYEo1WYcneSL7afJG0LA3WFmaM7liVwS0DsDAvxOUustJ065GdWat73nIUtJ8IZg+Oue3qNsbuGZvzfH6H+ZSxKVN4NQkhRAmj9yUzM7P//sGvUqlK1LgbuWQmsp27mcjY1Sc5cT0BgKYVXZnROwh/N/snvPMpJd2CVX0h+iiYWUL32VDv9Vy7bL6ymQ93fYhaUdPZvzPTW03H0qyQBnMLIUQxUGSr3QtRWmSoNXyzPYIFoZdQaxUcbSz4qGsNXmnkW/iDlGNO6+4kS7wOti7w8grwz73e2N+RfzNu9zg0ioauAV35rOVnWJjp/bUWQohS76l+cqanp2NjI7fxipLpyJU4Plx9kku3UwDoFFiOqb1qUc6pCP7OX9wEvw+CzGQoWxle/RXKVsq1y4ZLGxi/dzxaRUuPSj2Y0nwK5maFOKBbCCFKML0HPmg0GqZOnYqPjw8ODg5cvnwZgAkTJrB48WKDFyhEUUvOUDPxj9O8+N1+Lt1Owc3Bmvmv1ee7fg0KPwwpChxYAD+/ogtDAa1hyNZHwtC6iHV8vOdjtIqW3lV6M7XFVAlDQgjxFPQORJ999hlLly7liy++wMrKKmd77dq1WbRokUGLE6Ko7TgfS6dZO1m+/yqKAi82KM/W0a3pWtur8C+RabLgr9Hwz1hQtLplOF5fo7tc9pDVF1czce9EFBReqvoSk5pNwkxViIO6hRCiFND7ktny5cv5/vvv6dChA2+99VbO9qCgIM6fP2/Q4oQoKneTM5jy51n+CLsBgK+rLdOfD6JlFbeiKeAJy3Bk++X8L3x68FMA+lbvy7jG42TCRSGEMAC9A1F0dDSVK1d+ZLtWqyUrK8sgRQlRVBRF4Y+wG0z58yxxKZmYqWBQiwBGd6qKnVURDU5+wjIc2X469xMzDs0AoF9gPz5o+IGEISGEMBC9f+LXrFmT3bt34+fnl2v7b7/9Rr169QxWmBCFLTo+jfFrT7Hjwm0Aqns68nmfIOr4lim6Iq7uh1WvQlpcnstwZFt2Zhkzj8wEYGCtgYyqP0rCkBBCGJDegWjSpEn069eP6OhotFota9as4cKFCyxfvpw///yzMGoUwqC0WoUfD1zli3/Ok5KpwcrcjOHtK/Nmm0pYWRThWJwTq2D98DyX4XjY4lOLmX1sNgBv1H6D4fWGSxgSQggD0zsQde/enV9++YVp06ahUqmYOHEi9evXZ8OGDXTs2LEwahTCYCJik/hw9SmOXr0HQEM/F2b0qU1lD8eiK0KrhR2fwW7dGR/dMhzfg5XdI7t+d+I7vgn7BoC367zN23XeljAkhBCFQO+ZqksTmam65MhUa/l25yW+2R5BpkaLvZU5H3apzutN/DAzK8KAkY9lOEA3tmnBiQUsOLEAgGF1h/FmnTeLrk4hhCjGimSmaiGKm7CoeMauPsn5mCQA2lVz59Pna+NTxrZoC3lkGY45UO+1R3ZTFIV5x+ex8NRCAN6t/y6Daw8u2lqFEKKUkUAkSqzUTDVfbb7Ikr2RaBVwtbdiUvdAetTxLvrLTvlYhgN0YejrY1+z5PQSAN5v+D4Dag4o2lqFEKIUkkAkSqQ94XcYt/YkUXFpAPSq682E5wIp62Bd9MXkYxkO0IWhL498yY9nfwRgbOOxvFbj0TNIQgghDE8CUR5CQkIICQlBo9EYuxShp4TULD796yy/Hb0OgLezDZ89X5t21T2KvpjsZTg2f6ybeTqgNby0/JGZp3W7Ksw4NIOV51cCML7JeF6u/nJRVyyEEKWWDKp+DBlUXXwoisLfp2OY+McZ7iRnoFJBv6Z+jHm2Og7WRsj9miz4ewwc+UH3vP4A6PYVmFs+sqtW0fLZgc/49eKvqFAxsdlEXqj6QhEXLIQQJUehDaoePXp0vouYNWtWvvcVwhBuJaYzYd1pNp+9BUAld3s+7xNEQ39X4xSUz2U4QBeGpuyfwurw1ahQMaXFFHpV7lXEBQshhMhXIDp+/Hi+GpP5UURR0moVfjkSxbSN50hKV2NhpuKdtpUY2r4y1hZGWvk91zIc9veX4eia564arYaJ+yay/tJ6zFRmfNriU7pX6l7EBQshhIB8BqIdO3YUdh1C6OXKnRTGrjnJgctxANQp78znLwRR3dOIlzbzuQwHgFqrZvze8fx1+S/MVeZMazmNrhXzDk5CCCEKnwyqFsWKWqNl0Z5Ivt5ykQy1FhtLM97vVI2BLQIwL8oJFv8tn8twgC4MfbT7I/6+8jcWKgtmtJ5BZ//ORVuvEEKIXAoUiA4fPsxvv/3GtWvXyMzMzPXamjVrDFKYEP925kYCH64+yenoRABaVnZj2vO1qVD20SUviowey3AAZGmz+HDXh2y5ugULMwtmtp5JB78ORViwEEKIvOgdiFatWkX//v3p1KkTW7ZsoVOnToSHhxMTE8Pzzz9fGDWKUiw1U03gxE0AmJup0GgVnGwsGP9cIC82KG/ccWtZabD2LTi7Tve85WhoP+GRZThydtdk8f7O99ketR1LM0tmtZ1FW9+2RVauEEKI/6Z3IJo2bRpff/01Q4cOxdHRkTlz5hAQEMCbb76Jl1felwiEKKjwW0k5jzVaha61PfmkR008HG2MWBX5XoYjW6Ymk9Gho9l5fSdWZlbMbjebVuVbFWHBQgghHifvf8o+xqVLl+jWrRsA1tbWpKSkoFKpGDVqFN9//73BCxSlk6IorDx4jZe+O5CzbW7fusx/rYHxw1DMaVjYXheGbF2g/x+PDUPp6nRG7BjBzus7sTa3Zl77eRKGhBDCxOh9hsjV1ZWkJN2/2n18fDh9+jS1a9cmPj6e1NRUgxcoSp/E9CzGrTnFXydv5tr+TI1yRqroIflchiNbmjqNEdtHcODmAWwtbJnXfh5NvJoUYcFCCCHyQ+9A1KpVK7Zs2ULt2rV56aWXGDlyJNu3b2fLli106CCDQ8XTOREVz/Cfj3MtLhULMxUfdK7GG60qYmbMO8hAr2U4sqVmpTJ8+3AOxRzC1sKW+R3m09CzYREWLYQQIr/0DkTffPMN6enpAIwbNw5LS0v27NlD7969mTBhgsELFKWDoigs3hPJ5/+cJ0uj4FPGlnmv1qN+hf8OHEVGkwUbP4CjuhXoH7cMR7aUrBSGbhvK0VtHsbOw49uO31LPo14RFSyEEEJfspbZY8haZkUjLiWTD347wbbzsQA8W9OTz/sE4Wz334GjyDyyDMen0GxonstwZEvOTObtrW8TdjsMB0sHvu34LXXc856gUQghhOEV2lpmD9u4cSPm5uZ07px7IrnNmzej0Wjo0qWLvk2KUuxQZBwjfj5OTGI6VhZmTOhWg9eb+pnGMjBxl+8vw3HxictwZEvMTOTtLW9z8s5JHK0c+b7j99Ryq1VEBQshhCgove8yGzt2LBqN5pHtWq2WsWPHGqQoUfJptApzt4Xzyvf7iUlMp6KbPWvfaU6/Zv6mEYau7oeFHXRhyNEbBv3zxDCUkJHA/zb/j5N3TuJs7cyiToskDAkhRDGh9xmi8PBwAgMDH9levXp1IiIiDFKUKNliE9N595cw9l26C0Dv+j5M7VkLe2sTWUnm4WU4vOvBKz//5zIc2eLT43ljyxucjzuPi7ULCzstpJprtSIqWAghxNPS+zeQs7Mzly9fxt/fP9f2iIgI7O3tDVWXKKF2XrzN6F/CuJuSiZ2VOVN71qJPg/LGLktHz2U4ssWlxzFk8xDC74XjauPKok6LqOJSpQgKFkIIYSh6XzLr0aMH7777LpcuXcrZFhERwXvvvUePHj0MWpwoObI0Wmb8fZ4BPxzibkom1T0dWT+spemEocxU+H3ggzDUcjS8uPyJYehO2h0GbxpM+L1w3GzdWNJ5iYQhIYQohvQ+Q/Tll1/y7LPPUr16dcqX1/0yu379Oq1atWLmzJkGL1AUf9fvpTLi5+McuxYPwOtNKzC+WyA2lubGLSxb0i34+RW4cUy3DEePuVD31Se+7XbqbQZvHkxkQiQeth4s6ryIAOeAIihYCCGEoRXoktm+ffvYsmULJ06cwNbWlqCgIFq3bl0Y9Yli7p/TMYz5/QSJ6WocbSz4ok8QXWqb0Jp3Madg5SuQeF03yeLLP4F/iye+7VbKLYZsHsKVxCuUsyvHD51/oIJThSIoWAghRGGQeYgeQ+YhKrj0LA3TNp5j+f6rANT1LcO8vvXwdX38JagideEfWD34/jIcVeDVXx67DEe2m8k3Gbx5MFFJUXjZe7G482J8HX2LoGAhhBD5UWjzEM2dO5f//e9/2NjYMHfu3MfuO2LEiHwdWJRcl28nM2zlcc7eTATgzdYVeb9zNSzN9R6yVjgKsAxHtujkaAZvGkx0cjQ+Dj780PkHvB28i6BoIYQQhSlfZ4gCAgI4cuQIZcuWJSDgv8dIqFQqLl++bNACjUnOEOlv7fHrfLz2NKmZGlztrfjqpTq0q+Zh7LIeKMAyHNmiEqMYvHkwN1NuUsGxAos7L8bT3rOQCxZCCKGvQjtDFBkZmedjIbKlZqqZ+McZfj96HYCmFV2Z80o9yjnZGLmyh6TFw28D4HIo+V2GI9vVxKsM2jSI2NRY/J38WdRpEeXsyxV2xUIIIYqIicyEZ1pCQkIICQnJc0Zu8ahzNxMZtvIYl26nYKaCER2qMLx9FcyNvUL9wwqwDEe2ywmXGbJpCLfTblPRuSKLOy/GzdatkAsWQghRlPQeVD169Oi8G1KpsLGxoXLlyvTs2RNXV1eDFGhMcsns8RRF4aeD15jy51ky1VrKOVkz++V6NKtU1til5XZ1H6x6DdLiwMkH+q4Cr6B8vfVS/CUGbxrM3fS7VC5TmUWdFlHW1sT6J4QQIpeC/P7WOxC1a9eOY8eOodFoqFatGoqiEB4ejrm5OdWrV+fChQuoVCr27NmT5xIfxYkEov+WmJ7FuNWn+OvUTQDaVnPnqxfrUNbB2siV/UvYz7BhxINlOPquAsf8jfu5eO8ib2x+g7j0OKq5VGNhp4W42Dx54LUQQgjjKsjvb71v++nZsyfPPPMMN27c4OjRoxw7dozo6Gg6duxI3759iY6OpnXr1owaNUrvDojiISwqnm5zd/PXqZtYmKn4qGt1fhjQyLTCkFYL26bCurd0YahGDwjemO8wdD7uPIM3DSYuPY4arjVY3HmxhCEhhCjB9D5D5OPjw5YtWx45+3PmzBk6depEdHQ0x44do1OnTty5c8egxRY1OUOUm1arsHhPJJ//cx61VqG8iy3z+tajXgUTCwqZqbDubTi7Tve85WhoPwHM8pf/z9w9w/82/4/EzERqla3Ftx2/xdnaufDqFUIIYVCFdpfZwxISEoiNjX0kEN2+fZvERN28M2XKlCEzM1PfpoUJi0vJ5L1fw9hx4TYAXWt7Mr13EM62T75dvUglxcDPffVehiPbqduneHPLmyRlJRHkHsS3z3yLo5VjIRYshBDCFOgdiHr27MmgQYP46quvaNSoESqVikOHDvH+++/Tq1cvAA4dOkTVqlUNXaswkgOX7zJy1XFuJWZgZWHGxOcCea1JBVT5uF29SBVwGY5sYbFhvL31bZKzkqnnUY/5HebjYOVQiAULIYQwFXpfMktOTmbUqFEsX74ctVoNgIWFBQMGDODrr7/G3t6esLAwAOrWrWvoeotUab9kptEqfLM9gjnbLqJVoKK7Pd/0rU+gtwn+WRRwGY5sR28d5Z2t75CqTqVhuYaEdAjBztKElhkRQgiRb0Vyl1m25ORkLl++jKIoVKpUCQeHkvcv6dIciG4lpvPuqjD2X74LQJ/65ZnSsyb21iY2dZWiwIH5sOljQIGANvDSsnwtw5HtcMxhhm4bSpo6jSaeTZjbfq6EISGEKMaKZAxRNgcHB1xdXVGpVCUyDJVmoRdiee/XE9xNycTOypxPe9Wid/3yxi7rUf9ehqNBMHSdma9lOLIduHmA4duGk65Jp7l3c+a0m4ONhQnNri2EEKJI6H3bvVarZcqUKTg7O+Pn50eFChUoU6YMU6dORavVFkaNoohkabRM33iO4CWHuZuSSQ0vJzYMb2maYSgtHn564X4YUkGnz+C52XqFob3Rexm2bRjpmnRa+bRibvu5EoaEEKKU0vsM0ccff8zixYuZMWMGLVq0QFEU9u7dyyeffEJ6ejqfffZZYdQpCllUXCrDfz5OWFQ8AP2b+fFR1xrYWJobt7C8/HsZjhcWQ7UuejWx6/ou3t3xLlnaLNr6tuWrNl9hZW5VSAULIYQwdXqPIfL29ubbb7+lR48eubb/8ccfvPPOO0RHRxu0QGMqLWOI/j51kzGrT5KUrsbRxoIvXwji2Vpexi4rb0+xDEe2Hdd2MHrnaNRaNc9UeIYvWn+BpR5nloQQQpi2IhlDFBcXR/Xq1R/ZXr16deLi4vRtThhRepaGz/46x48HrgJQr0IZ5r5SD19XExpQnJkC07x1j7vPgb/eB22W3stwZNt6dSsf7PwAtaKmk18nZrSegaWZhCEhhCjt9B5DVKdOHb755ptHtn/zzTfUqVPHIEWJwnfpdjLPz9+XE4bebFORX99sZlph6N82jNSFIT2X4cj2z5V/eH/n+6gVNV0CuvB5688lDAkhhAAKcIboiy++oFu3bmzdupVmzZqhUqnYt28fUVFRbNy4sTBqFAa2+uh1JvxxmtRMDWXtrfjqpTq0reZh7LLyplXnft7qPWg3Pt/LcGT76/JffLTnI7SKlu4VuzO1xVTMzUxwfJQQQgij0DsQtWnThosXLxISEsL58+dRFIXevXvzzjvv4O3tXRg1CgNJyVAz4Y/TrDmmG+fVrGJZZr9Sl3JOJnpnlVajOyuU7bnZ0HCg3s2sv7SeCXsnoFW09Krci0+afSJhSAghRC4FmofI29v7kbvJoqKiGDRoED/88INBChOGdfZGIsN+Psbl2ymYqeDdZ6oytF1lzM1MbPmNbFoN/DEUzqx9sC3oJb2bWRu+lkn7JqGg8ELVF5jQdAJmKr2vFAshhCjhDPabIS4ujmXLlhmqOWEgiqLw44Gr9Jq/l8u3UyjnZM3PbzRlRIcqJhyGtLBhBJz4GVQFP5Pz28XfmLhvIgoKL1d7WcKQEEKI/2Ri6zAIQ0pIy2Ls6pP8fToGgPbVPZj5Yh1c7U14vh1Fgb9Gw/EVoDKDPguhVh+9m/n5/M9MOzgNgNdrvM6YRmNMbzFaIYQQJkMCUQl1/No9hv98nOv30rA0V/Hhs9UZ3DLAtEOBosDfYx7MPt3r2wKFoRVnV/D54c8BCK4ZzOgGo02730IIIYxOAlEJo9UqLNx9mS83XUCtVfB1teWbvvWp41vG2KU9nqLoFmg99D2ggp4hUOdlvZtZenopXx39CoAhtYcwot4ICUNCCCGeKN+BqHfv3o99PT4+/mlrEU/pbnIG7/12gtALtwHoFuTF9N61cbIx8bl2FAW2ToIDIbrn3WdDvdf0bmbRqUXMOTYHgDeD3mRo3aEShoQQQuRLvgORs7PzE1/v37//UxckCmb/pbu8+8txbiVmYG1hxqTuNenb2Ld4BIIdn8FeXZCh60zdqvV6WnBiAfPD5gPwTt13eLvO2wYsUAghREmX70C0ZMmSwqxDFJBGqzB3WzjztoejVaCSuz0hr9WnumcxWXst9HPY9aXu8bOfQ+M39Hq7oiiEhIXw3cnvABhZfyRDag8xdJVCCCFKOBlDVIzFJKQzctVxDkbq1pB7sUF5JvesiZ1VMflYd38Fobo7wej0KTR9S6+3K4rC3ONzWXRqEQDvNXiP4FrBBi5SCCFEaVBMfnOKf9txPpb3fjtBXEom9lbmfPZ8bXrV8zF2Wfm3dy5sm6J73GESNB+u19sVRWHW0VksPbMUgDGNxtAvsJ+BixRCCFFaSCAqZjLVWr7cdJ6FuyMBqOntxLy+9ajo7mDkyvRwYAFsmaB73PYjaDVar7crisIXh79gxbkVAHzU5CP6Vu9r6CqFEEKUIhKIipGouFSG/XycE1HxAAQ392dc1+pYWxSjdbkOLYR/xuoet/4A2n6o19u1ipbpB6ez6sIqACY2m8iLVV80dJVCCCFKGQlExcTGUzf5cPVJktLVONlY8MULdXi2lqexy9LPkSWw8X3d4xbvQruP9Xq7VtEy9cBUfr/4OypUTG4+meerPG/4OoUQQpQ6EohMXHqWhql/nuWng9cAqF+hDHP71qO8i52RK9PT8RXw57u6x82GwTOfQD6mBEjNSqXJyiYAdK/YnQ2XN2CmMmNqi6n0qNSj8OoVQghRqkggMmERsckMW3mM8zFJALzdthKjO1bF0ryYLVB64hf4Y5juceM3dXeUFWB+pOwwNK3lNLpV7GbgIoUQQpRmEohM1O9HrzNh3WnSsjS4OVgx66W6tK7qbuyy9Hfqd1j3FqBAw0HQ5XO9wpCiKDmPzVXmzGg9g2f9ny2EQoUQQpRmEohMTHKGmonrTrPmeDQALSqX5euX6uLhZGPkygrgzDpY8z9QtFC/P3T9Su8zQ8vOLst5/GmLTyUMCSGEKBQSiEzImRsJDF95nMt3UjBTweiOVXm7bWXMzYrB8hv/dv4vWD0YFA3UeRWemwNm+l3q23Z1GwtOLMh53r5Ce0NXKYQQQgASiEyCoij8eOAqn/51jky1Fi9nG+a8Uo/GAa7GLq1gLm6CXweAVg21X4Se3+gdhi7EXWDcnnGFVKAQQgiRmwQiI0jNVBM4cRMA+8e1Z/L6s/xzJgaADtU9mPliHVzsrYxZYsFFbIVfXgdtFtR8Hnp9C2b6zZN0J+0Ow7YPI02dRmPPxhyKOVRIxQohhBA6EoiMrM+CfdyIT8fSXMXYLjUY1MK/eKxQn5fLobDqNdBkQo3u0HshmOv3VyxTk8moHaOISYnB38mfWW1n4WztXDj1CiGEEPcVs/u39RcVFUXbtm0JDAwkKCiI3377zdglodU+uHPqRnw6FVztWP12cwa3DCi+YejKHlj5CqjToWoX6PMDmFvq1YSiKEzeP5mw22E4Wjkyr/08CUNCCCGKRIk/Q2RhYcHs2bOpW7cusbGx1K9fn65du2Jvb2+0ms7FJOY87lLLk89fCMLJRr/wYFKu7oefXgJ1GlTuCC8tAwv9L/ktPbOU9ZfWY64y56s2X+Hv7G/4WoUQQog8lPhA5OXlhZeXFwAeHh64uroSFxdn1EBU0/vBWY+ZLwZhb12Mw1DUYfjpBchKgYrt4OUVYGGtdzM7o3by9dGvAfiw8Yc0825m6EqFEEKI/2T0S2a7du2ie/fueHt7o1KpWLdu3SP7zJ8/n4CAAGxsbGjQoAG7d+8u0LGOHDmCVqvF19f3Kas2nGJ7iQwg+his6A2ZyeDfCl5ZCZb6z5cUfi+cMbvGoKDwUtWXeKXaK4VQrBBCCPHfjH6GKCUlhTp16jBw4ED69OnzyOu//PIL7777LvPnz6dFixZ89913dOnShbNnz1KhQgUAGjRoQEZGxiPv3bx5M97e3gDcvXuX/v37s2jRov+sJSMjI1c7iYmJ/7lvqXfzBPzYCzISoUJzePUXsNJ/fbW49DiGbx9OqjqVxp6NGdtkbPEOiUIIIYollfLw2ghGplKpWLt2Lb169crZ1qRJE+rXr8+CBQ8m6KtRowa9evVi+vTp+Wo3IyODjh078sYbb9CvX7//3O+TTz5h8uTJj2xPSEjAyckp/x0p6WJOw7LnIO0e+DaB11eDtaPezWRpshiyeQjHYo/h6+jLyq4rKWNTxvD1CiGEKFUSExNxdnbW6/e30S+ZPU5mZiZHjx6lU6dOubZ36tSJffv25asNRVEIDg6mffv2jw1DAOPGjSMhISHnv6ioqALXXmLFnoPlPXRhyKcBvPZ7gcKQoih8evBTjsUew8HSgW/afyNhSAghhNEY/ZLZ49y5cweNRkO5cuVybS9XrhwxMTH5amPv3r388ssvBAUF5YxP+vHHH6ldu/Yj+1pbW2Ntrf+A4FLj9kVY1gNS74JXXXh9DdgU7MzZinMrWBO+BjOVGV+2+ZKKZSoatlYhhBBCDyYdiLL9e0yJoij5HmfSsmVLtFptYZRVuty9BMu6Q0oslKsN/daCbZkCNbX7+m5mHpkJwPsN36elT0sDFiqEEELoz6Qvmbm5uWFubv7I2aDY2NhHzhqJQhQXqQtDyTHgEQj9/wC7gq2zdjn+MmN2jUGraOldpTev13jdwMUKIYQQ+jPpQGRlZUWDBg3YsmVLru1btmyhefPmRqqqlLl3VReGEqPBrRr0Xw/2ZQvUVHx6PMO2DyM5K5n6HvUZ32S83FEmhBDCJBj9kllycjIRERE5zyMjIwkLC8PV1ZUKFSowevRo+vXrR8OGDWnWrBnff/89165d46233jJi1aVEwnVdGEqIgrKVYcB6cHAvUFNZ2ize2/keUUlR+Dj48HW7r7HUc2kPIYQQorAYPRAdOXKEdu3a5TwfPXo0AAMGDGDp0qW8/PLL3L17lylTpnDz5k1q1arFxo0b8fPzK7SaQkJCCAkJQaPRFNoxTF7iDV0Yir8KLgEwYAM4ehaoKUVRmHFwBodiDmFnYce89vNwtSnYJTchhBCiMJjUPESmpiDzGJQISbdgaVe4GwFlKkDwRihT8Nm9fz7/M9MOTkOFinnt59HGt40BixVCCCFyK3HzEAkjSL6tm2fobgQ4+8KAP58qDO27sY/PD30OwKgGoyQMCSGEMEkSiMQDKXdheU+4fR4cvXVjhlwKfmnySsIV3t/5PhpFQ49KPQiuGWy4WoUQQggDkkAkdFLj4MeeEHsGHDx1Y4ZcCz5ZYkJGAsO3DycpM4k67nWY1GyS3FEmhBDCZEkgEpAWDz8+DzGnwN5dF4bcKhe4ObVWzQc7P+BK4hW87L2Y3W42VuZWhqtXCCGEMDAJRKVdeiKs6AM3w8CurC4MuVd9qia/PPwl+2/ux9bClnnt5+Fm62aYWoUQQohCIoEoDyEhIQQGBtKoUSNjl1K4MpLhpxch+gjYuuhmoPao8VRN/nrhV1aeXwnA9JbTqeZazRCVCiGEEIVKbrt/jBJ9231mii4MXd0LNs66Gai96z5Vk4duHuLNLW+iVtSMqDeCN4LeMEytQgghhB7ktnuRP5mp8PMrujBk7aRbqPUpw9C1xGuM3jkataKma0BXhtQeYphahRBCiCIggai0yUqHX16DyF1g5QCvrwafBk/VZFJmEsO3DychI4HabrWZ3Hyy3FEmhBCiWJFAVJqoM+DXfnBpO1jawWu/g2/jp2pSo9UwZtcYLidcxsPOgznt5mBjYWOggoUQQoiiIYGotFBnwm/BEL4ZLGzh1V/Br9lTNzvr6Cz2RO/BxtyGue3n4m5XsMVfhRBCCGOSQFQaaLJg9SC4sBEsbKDvzxDQ6qmbXRO+huVnlwPwactPqVm25lO3KYQQQhiDBKKSTqOGNf+DcxvA3Ape+QkqtXvqZo/eOsrUA1MBeKfOO3T27/zUbQohhBDGIoGoJNNqYN3bcGYNmFnCSz9C5WeeutnrSdcZtWMUaq2azv6deavOWwYoVgghhDAeCUR5KBETM2q1sH44nPoVzCzgxaVQ7dmnbjYlK4Xh24dzL+MegWUDmdpiqtxRJoQQotiTiRkfo9hOzKjVwp8j4dhyUJnDCz9AzV5P3axGq+HdHe8Sej0Ud1t3fu72M+Xsyz19vUIIIYQBycSMAhQFNr5/PwyZQe/vDRKGAOYen0vo9VCsza2Z026OhCEhhBAlhgSikkRR4J+xcGQxoIJeC6D2CwZpev2l9fxw+gcApjSfQm332gZpVwghhDAFEohKCkWBzePh4Le65z3mQZ1XDNJ0WGwYn+z7BIA3ar9B14pdDdKuEEIIYSokEJUEigLbJsP+b3TPn5sN9fsZpOmbyTcZuWMkWdosOlTowLB6wwzSrhBCCGFKJBCVBKHTYc/XusddZ0LDgQZpNjUrleHbhxOXHkc1l2pMazkNM5X8lRFCCFHyyG+34m7nl7Dzc93jztOh8RsGaVaraPloz0dcuHcBVxtX5rWfh52lnUHaFkIIIUyNBKLibM/XsONT3eOOU6DZOwZrOiQshG3XtmFpZsmcdnPwcvAyWNtCCCGEqZFAlIdiMTHjvm9g6ye6x+0nQIuRBmt64+WNfH/yewA+af4JdT3qGqxtIYQQwhTJxIyPYbITMx78Dv4eo3vcdhy0HWuwpk/dPkXwP8FkajMZWGsgoxuMNljbQgghRFGQiRlLg8OLH4ShVu9Dmw8N1nRMSgwjdowgU5tJ2/JtGVnPcGedhBBCCFMmgag4OboM/rp/xqb5CGg/Hgy0jliaOo2RO0ZyJ+0OlctUZkbrGZibmRukbSGEEMLUSSAqLsJWwob7Z2yavqMbRG2gMKQoCuP3jOfs3bO4WLswr/087C3tDdK2EEIIURxIICoOTv4G694BFGj8P+g8zWBhCODbE9+y+epmLMws+Lrd15R3LG+wtoUQQojiQAKRqTu9Btb+D1CgwUDo8oVBw9CmK5uYf2I+ABOaTqBBuQYGa1sIIYQoLiQQmbKz62H1EFC0UO916DbLoGHozN0zjN8zHoB+gf3oXaW3wdoWQgghihMJRKbqwt/w+0BQNBD0CnSfC2aG+7hup95mxPYRpGvSaeHTgvcavGewtoUQQojiRgKRKQrfAr/2B60aar0AveaDAe/4SlenM3LHSGJTY6noXJEvW38pd5QJIYQo1SQQmZqIbbDqNdBkQmBPeP47g4YhRVGYtG8Sp+6cwtnamXnt5+Fo5Wiw9oUQQojiSAKRKbm8E1a9CpoMqP4c9FkM5hYGPcSiU4vYGLkRC5UFs9rMooJTBYO2L4QQQhRHEojyYJS1zK7shZ9fAXU6VH0WXlgC5pYGPcS2q9uYe3wuAOOajKOxV2ODti+EEEIUV7KW2WMU2lpmmSkwzVv3+KMbEHMafnweslKg8jPwykqwsDbc8YALcRfo93c/0tRpvFr9VcY1GWfQ9oUQQghTIWuZFUfRx2BFH10YqtgWXl5h8DB0J+0Ow7YPI02dRjOvZnzQ6AODti+EEEIUdxKIjG3Vq5CZBP6t4JWfwdLWoM1najIZtWMUMSkx+Dv582WbL7EwM+y4JCGEEKK4k0BkbBmJUKEZ9F0FVnYGbVpRFCbvn0zY7TAcrRyZ134eztbOBj2GEEIIURJIIDKGO+EPHvs0gNd+A2sHgx9m6ZmlrL+0HnOVOTPbzMTf2d/gxxBCCCFKAglExuDo+eDxyz+BteHnAdoZtZOvj34NwJhGY2ju3dzgxxBCCCFKCglExvBwALIx4N1r94XfC2fMrjEoKLxU9SX6Vu9r8GMIIYQQJYkEohImLj2O4duHk6pOpbFnY8Y2GYvKgAvCCiGEECWRBKISJEuTxagdo4hOjsbX0Zev2nyFpZlhJ3cUQgghSiIJRCWEoih8evBTjsUew8HSgW/af0MZmzLGLksIIYQoFmRCGmOwsodPEgza5IpzK1gTvgYzlRlftP6CimUqGrR9IYQQoiSTM0QlwO7ru5l5ZCYA7zV4j1blWxm5IiGEEKJ4kUBUzF2Ov8yYXWPQKlp6V+lNv8B+xi5JCCGEKHYkEOXBKKvdF0B8ejzDtg8jOSuZ+h71Gd9kvNxRJoQQQhSArHb/GIW22r0BZGmzeGvLWxyKOYSPgw8ru63E1cbV2GUJIYQQRier3ZcSiqIw4+AMDsUcws7Cjnnt50kYEkIIIZ6CBKJiaNWFVfx68VdUqPi89edUcali7JKEEEKIYk0CUTGz78Y+Pj/0OQDvNniXtr5tjVuQEEIIUQJIICpGriRc4f2d76NRNPSo1IOBNQcauyQhhBCiRJBAVEwkZCQwfPtwkjKTqONeh4nNJsodZUIIIYSBSCAqBtRaNR/s/IAriVfwtPdkdrvZWJtbG7ssIYQQosSQQFQMfHn4S/bf3I+thS3z2s/DzdbN2CUJIYQQJYoEIhP364VfWXl+JQDTW06numt1I1ckhBBClDwSiEzYoZuHmH5wOgDD6w2ng18HI1ckhBBClEwSiEzUtcRrjN45GrWipmtAV96o/YaxSxJCCCFKLAlEJigpM4nh24eTkJFAbbfaTG4+We4oE0IIIQqRBCITo9FqGLNrDJcTLuNh58GcdnOwsbAxdllCCCFEiSaByMTMOjqLPdF7sDG3YW77ubjbuRu7JCGEEKLEk0BkQtaEr2H52eUAfNryU2qWrWnkioQQQojSQQKREaRmpVJ7WW1qL6tNalYqAEdvHWXqgakAvFPnHTr7dzZmiUIIIUSpIoHIBFxPus6oHaNQa9V08uvEm3XeNHZJQgghRKkigSgPISEhBAYG0qhRo0I/VkpWCsO3D+dexj1quNbg05afYqaSj0UIIYQoSipFURRjF2GqEhMTcXZ2JiEhAScnJ4O1m5qVSpOVTQBo5dOK3dG7cbd1Z2W3lXjaexrsOEIIIURpVJDf33Iqwsh2R+/GysyKOe3mSBgSQgghjEQCkQmY0mIKtd1rG7sMIYQQotSSQGQEF+9dzHk8sOZAulXsZsRqhBBCCCGByAgCnAJyHv8v6H9GrEQIIYQQIIHIKCzNLXMeyx1lQgghhPHJb2MhhBBClHoSiIQQQghR6kkgEkIIIUSpJ4FICCGEEKWeBCIhhBBClHoWxi6gNLKztOPUgFPGLkMIIYQQ98kZIiGEEEKUehKIhBBCCFHqSSASQgghRKkngUgIIYQQpZ4EIiGEEEKUehKIhBBCCFHqSSASQgghRKkngUgIIYQQpZ4EIiGEEEKUehKIhBBCCFHqSSASQgghRKkngUgIIYQQpZ4EojyEhIQQGBhIo0aNjF2KEEIIIYqASlEUxdhFmKrExEScnZ1JSEjAycnJ2OUIIYQQIh8K8vtbzhAJIYQQotSzMHYBpiz75FliYqKRKxFCCCFEfmX/3tbnIpgEosdISkoCwNfX18iVCCGEEEJfSUlJODs752tfGUP0GFqtlhs3buDo6IhKpTJo24mJifj6+hIVFVUixyeV9P5Bye+j9K/4K+l9LOn9g5Lfx8Lqn6IoJCUl4e3tjZlZ/kYHyRmixzAzM6N8+fKFegwnJ6cS+Zc8W0nvH5T8Pkr/ir+S3seS3j8o+X0sjP7l98xQNhlULYQQQohSTwKREEIIIUo9CURGYm1tzaRJk7C2tjZ2KYWipPcPSn4fpX/FX0nvY0nvH5T8PppS/2RQtRBCCCFKPTlDJIQQQohSTwKREEIIIUo9CURCCCGEKPUkEAkhhBCi1JNAlE/z588nICAAGxsbGjRowO7dux+7/86dO2nQoAE2NjZUrFiRb7/9NtfrCxcupFWrVri4uODi4sIzzzzDoUOH9D6uoih88skneHt7Y2trS9u2bTlz5kyx6N/06dNp1KgRjo6OeHh40KtXLy5cuJBrn+DgYFQqVa7/mjZtWiz698knnzxSu6enZ659DPX5GauP/v7+j/RRpVIxdOjQnH1M9TNcs2YNDRs2pEyZMtjb21O3bl1+/PFHvY9rqt/B/PSvKL+DxupjUX4PjdG/4vwdfNiqVatQqVT06tVL7+Ma7OeoIp5o1apViqWlpbJw4ULl7NmzysiRIxV7e3vl6tWree5/+fJlxc7OThk5cqRy9uxZZeHChYqlpaXy+++/5+zz6quvKiEhIcrx48eVc+fOKQMHDlScnZ2V69ev63XcGTNmKI6Ojsrq1auVU6dOKS+//LLi5eWlJCYmmnz/OnfurCxZskQ5ffq0EhYWpnTr1k2pUKGCkpycnLPPgAEDlGeffVa5efNmzn93797Nd9+M2b9JkyYpNWvWzFV7bGxsrmMZ4vMzZh9jY2Nz9W/Lli0KoOzYsSNnH1P9DHfs2KGsWbNGOXv2rBIREaHMnj1bMTc3V/755x+9jmuq38H89K+ovoPG7GNRfQ+N1b/i/B3MduXKFcXHx0dp1aqV0rNnT72Pa6ifoxKI8qFx48bKW2+9lWtb9erVlbFjx+a5/5gxY5Tq1avn2vbmm28qTZs2/c9jqNVqxdHRUVm2bFm+j6vVahVPT09lxowZOa+np6crzs7Oyrfffpu/zuXjOP9mqP79W2xsrAIoO3fuzNk2YMCAR74g+jJW/yZNmqTUqVPnP99jqM9PUUznMxw5cqRSqVIlRavV5mwrLp+hoihKvXr1lPHjx+f7uMXpO6goj/bv3wrrO6goxutjUX0PTeUzLG7fQbVarbRo0UJZtGhRnnUW1XdQURRFLpk9QWZmJkePHqVTp065tnfq1Il9+/bl+Z79+/c/sn/nzp05cuQIWVlZeb4nNTWVrKwsXF1d833cyMhIYmJicu1jbW1NmzZt/rM2U+lfXhISEgAe2Sc0NBQPDw+qVq3KG2+8QWxs7BP7lc3Y/QsPD8fb25uAgABeeeUVLl++nPOaIT4/U+jjw3WsWLGCQYMGPbIYsql/hoqisG3bNi5cuEDr1q3zfdzi8h3Mq395KYzvIBi/j4X9PTR2/x6uo7h9B6dMmYK7uzuDBw8u0HEN9XMUZAzRE925cweNRkO5cuVybS9XrhwxMTF5vicmJibP/dVqNXfu3MnzPWPHjsXHx4dnnnkm38fN/r8+tZlK//5NURRGjx5Ny5YtqVWrVs72Ll268NNPP7F9+3a++uorDh8+TPv27cnIyDD5/jVp0oTly5ezadMmFi5cSExMDM2bN+fu3bs5x8luO7+1mVofH7Zu3Tri4+MJDg7Otd2UP8OEhAQcHBywsrKiW7duzJs3j44dO+b7uKb+HXxc//6tsL6Dxu5jUXwPTeUzLG7fwb1797J48WIWLlxY4OMa6ucoyGr3+fbvtK0oyiPbnrR/XtsBvvjiC37++WdCQ0OxsbHR+7j61pbfeouif9mGDRvGyZMn2bNnT67tL7/8cs7jWrVq0bBhQ/z8/Pjrr7/o3bv34zv1hHoLu39dunTJeVy7dm2aNWtGpUqVWLZsGaNHjy5wbfrUXJSf4eLFi+nSpQve3t65tpvyZ+jo6EhYWBjJycls27aN0aNHU7FiRdq2bavXcU31O5if/mUr7O/gf9Vc2H0syu+hsT/D4vQdTEpK4vXXX2fhwoW4ubk99XEN8flJIHoCNzc3zM3NH0masbGxjyTSbJ6ennnub2FhQdmyZXNtnzlzJtOmTWPr1q0EBQXpddzsOyViYmLw8vLKV22m0r+HDR8+nPXr17Nr1y7Kly//2Hq9vLzw8/MjPDz8SV0DTKN/2ezt7aldu3ZO7Yb4/MA0+nj16lW2bt3KmjVrnlivKX2GZmZmVK5cGYC6dety7tw5pk+fTtu2bUvEd/Bx/XtYYX4HwTT6mK0wvoem0L/i9h08c+YMV65coXv37jmva7VaACwsLLhw4QK+vr5F8h3MJpfMnsDKyooGDRqwZcuWXNu3bNlC8+bN83xPs2bNHtl/8+bNNGzYEEtLy5xtX375JVOnTuWff/6hYcOGeh83ICAAT0/PXPtkZmayc+fO/6zNVPoHugQ/bNgw1qxZw/bt2wkICHhivXfv3iUqKirXX/zHMWb//i0jI4Nz587l1G6Iz89U+rhkyRI8PDzo1q3bE+s1pc/w3xRFybmMUBK+g4/rX/bzwv4OgnH7+G+F8T00hf4Vt+9g9erVOXXqFGFhYTn/9ejRg3bt2hEWFoavr2+RfQdz6DUEu5TKvu1v8eLFytmzZ5V3331Xsbe3V65cuaIoiqKMHTtW6devX87+2bcbjho1Sjl79qyyePHiR243/PzzzxUrKyvl999/z3UrZFJSUr6Pqyi62w2dnZ2VNWvWKKdOnVL69u1b4NtFi7p/b7/9tuLs7KyEhobm2ic1NVVRFEVJSkpS3nvvPWXfvn1KZGSksmPHDqVZs2aKj49Psejfe++9p4SGhiqXL19WDhw4oDz33HOKo6OjwT8/Y/ZRURRFo9EoFSpUUD788MNH6jLlz3DatGnK5s2blUuXLinnzp1TvvrqK8XCwkJZuHBhvo+rKKb7HcxP/4rqO2jMPhbV99BY/VOU4vsd/Le87jIrqu+gosht9/kWEhKi+Pn5KVZWVkr9+vUfuS21TZs2ufYPDQ1V6tWrp1hZWSn+/v7KggULcr3u5+enAI/8N2nSpHwfV1F0txxOmjRJ8fT0VKytrZXWrVsrp06dKhb9y+t1QFmyZImiKIqSmpqqdOrUSXF3d1csLS2VChUqKAMGDFCuXbtWLPqXPReGpaWl4u3trfTu3Vs5c+ZMrnYM9fkZq4+KoiibNm1SAOXChQuP1GTKn+HHH3+sVK5cWbGxsVFcXFyUZs2aKatWrdLruIpiut/B/PSvKL+DxupjUX4PjfV3tLh+B//tv6YHKKrvoEpR7o9yEkIIIYQopWQMkRBCCCFKPQlEQgghhCj1JBAJIYQQotSTQCSEEEKIUk8CkRBCCCFKPQlEQgghhCj1JBAJIYQQotSTQCSEEEKIUk8CkRDiiUJDQ1GpVMTHxxdK+1euXEGlUhEWFmawNlUqFevWrTNYe0KIkk0CkRCC4OBgVCoVKpUKS0tLKlasyPvvv09KSkqRHN/X15ebN29Sq1atIjleSVTYoVWIks7C2AUIIUzDs88+y5IlS8jKymL37t0MGTKElJQUFixYUOjHNjc3x9PTs9CP8zQURUGj0WBhkfvHZmZmJlZWVnq3l9/3FbR9IYR+5AyREAIAa2trPD098fX15dVXX+W111575JLT0aNHadiwIXZ2djRv3pwLFy4AukteZmZmHDlyJNf+8+bNw8/PD0VRuHfvHq+99hru7u7Y2tpSpUoVlixZkvP+f18yO3PmDN26dcPJyQlHR0datWrFpUuXADh8+DAdO3bEzc0NZ2dn2rRpw7Fjx/Tqr6IofPHFF1SsWBFbW1vq1KnD77//nvN69hmXTZs20bBhQ6ytrdm9ezdt27Zl2LBhjB49Gjc3Nzp27AjAzp07ady4MdbW1nh5eTF27FjUanVOe//1vn8LDg6mV69eTJ8+HW9vb6pWrQrAihUraNiwIY6Ojnh6evLqq68SGxub8+fXrl07AFxcXFCpVAQHB+ern0IIHQlEQog82drakpWVlWvbxx9/zFdffcWRI0ewsLBg0KBBAPj7+/PMM8/kBJxsS5YsybkcN2HCBM6ePcvff//NuXPnWLBgAW5ubnkeOzo6mtatW2NjY8P27ds5evQogwYNygkYSUlJDBgwgN27d3PgwAGqVKlC165dSUpKynf/xo8fz5IlS1iwYAFnzpxh1KhRvP766+zcuTPXfmPGjGH69OmcO3eOoKAgAJYtW4aFhQV79+7lu+++Izo6mq5du9KoUSNOnDjBggULWLx4MZ9++mmutv79vv+ybds2zp07x5YtW/jzzz8B3ZmiqVOncuLECdatW0dkZGRO6PH19WX16tUAXLhwgZs3bzJnzhy9+ilEqacIIUq9AQMGKD179sx5fvDgQaVs2bLKSy+9pCiKouzYsUMBlK1bt+bs89dffymAkpaWpiiKovzyyy+Ki4uLkp6eriiKooSFhSkqlUqJjIxUFEVRunfvrgwcODDP40dGRiqAcvz4cUVRFGXcuHFKQECAkpmZma/61Wq14ujoqGzYsCFnG6CsXbs2z/2Tk5MVGxsbZd++fbm2Dx48WOnbt2+uPq9bty7XPm3atFHq1q2ba9tHH32kVKtWTdFqtTnbQkJCFAcHB0Wj0fzn+/IyYMAApVy5ckpGRsZj9zt06JACKElJSbnqvXfvnl79FELoyBkiIQQAf/75Jw4ODtjY2NCsWTNat27NvHnzcu2TfYYEwMvLCyDnsk2vXr2wsLBg7dq1APzwww+0a9cOf39/AN5++21WrVpF3bp1GTNmDPv27fvPWsLCwmjVqhWWlpZ5vh4bG8tbb71F1apVcXZ2xtnZmeTkZK5du5avvp49e5b09HQ6duyIg4NDzn/Lly/PuSyXrWHDho+8/9/bzp07R7NmzVCpVDnbWrRoQXJyMtevX39sW3mpXbv2I+OGjh8/Ts+ePfHz88PR0ZG2bdsCPLbP+vRTiNJOBlULIQBo164dCxYswNLSEm9v7zzDyMPbsn/5a7VaAKysrOjXrx9Lliyhd+/erFy5ktmzZ+fs36VLF65evcpff/3F1q1b6dChA0OHDmXmzJmPHMfW1vaxtQYHB3P79m1mz56Nn58f1tbWNGvWjMzMzHz1Nbvmv/76Cx8fn1yvWVtb53pub2//yPv/vU1RlFxhKHsbkGt7Xm3l5d/7paSk0KlTJzp16sSKFStwd3fn2rVrdO7c+bF91qefQpR2EoiEEIDul3DlypWfqo0hQ4ZQq1Yt5s+fT1ZWFr179871uru7O8HBwQQHB9OqVSs++OCDPANRUFAQy5YtIysrK89gtnv3bubPn0/Xrl0BiIqK4s6dO/muMzAwEGtra65du0abNm307GXe7a1evTpXMNq3bx+Ojo6PBJGCOH/+PHfu3GHGjBn4+voCPDKAPfuMkkajyVWXIfspREkml8yEEAZTo0YNmjZtyocffkjfvn1znemZOHEif/zxBxEREZw5c4Y///yTGjVq5NnOsGHDSExM5JVXXuHIkSOEh4fz448/5tzVVrlyZX788UfOnTvHwYMHee211554Vulhjo6OvP/++4waNYply5Zx6dIljh8/TkhICMuWLdO73++88w5RUVEMHz6c8+fP88cffzBp0iRGjx6NmdnT/5itUKECVlZWzJs3j8uXL7N+/XqmTp2aax8/Pz9UKhV//vknt2/fJjk52eD9FKIkk0AkhDCowYMHk5mZmXMHWjYrKyvGjRtHUFAQrVu3xtzcnFWrVuXZRtmyZdm+fTvJycm0adOGBg0asHDhwpyzRT/88AP37t2jXr169OvXjxEjRuDh4aFXnVOnTmXixIlMnz6dGjVq0LlzZzZs2EBAQIDeffbx8WHjxo0cOnSIOnXq8NZbbzF48GDGjx+vd1t5cXd3Z+nSpfz2228EBgYyY8aMR86s+fj4MHnyZMaOHUu5cuUYNmyYwfspREmmUrIvdAshhAF89tlnrFq1ilOnThm7FCGEyDc5QySEMIjk5GQOHz7MvHnzGDFihLHLEUIIvUggEkIYxLBhw2jZsiVt2rR55HKZEEKYOrlkJoQQQohST84QCSGEEKLUk0AkhBBCiFJPApEQQgghSj0JREIIIYQo9SQQCSGEEKLUk0AkhBBCiFJPApEQQgghSj0JREIIIYQo9f4Pc9wbQfPJDWEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "plt.figure()\n",
    "for L, logical_errors in zip(Ls, log_errors_all_L):\n",
    "    std_err = (logical_errors*(1-logical_errors)/num_trials)**0.5\n",
    "    plt.errorbar(ps, logical_errors, yerr=std_err, label=\"L={}\".format(L))\n",
    "plt.yscale(\"log\")\n",
    "plt.xlabel(\"Physical error rate\")\n",
    "plt.ylabel(\"Logical error rate\")\n",
    "plt.legend(loc=0);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.6"
  },
  "vscode": {
   "interpreter": {
    "hash": "397704579725e15f5c7cb49fe5f0341eb7531c82d19f2c29d197e8b64ab5776b"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}