intel-qs/notebooks/noise_via_chi_matrix_example.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "----\n",
    "# Example using the chi-matrix to simulate noisy circuits\n",
    "----\n",
    "\n",
    "As a first example, we apply three known noise channels to an initial pure state. Each channel is applied multiple times sequentially and we compute how the squared overlap changes depending on the number of applications.\n",
    "\n",
    "When a single channel corresponds to a noisy evolution for a certain fixed duration, the number of channel applications plays the role of time.\n",
    "\n",
    "We consider the following channels:\n",
    "* depolarizing\n",
    "* dephasing\n",
    "* amplitude damping\n",
    "\n",
    "At the end of the notebook we include some quick test to confirm the implementation of class `ChiMatrix` and its methods.\n",
    "\n",
    "----\n",
    "\n",
    "### Import Intel QS library\n",
    "\n",
    "We start by importing the Python library with the class and methods defined in the C++ implementation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Import the Python library with the C++ class and methods of Intel Quantum Simulator.\n",
    "# Since the library is not contained in the same folder of this notebook, its path has to be added.\n",
    "import sys\n",
    "sys.path.insert(0, '../build/lib')\n",
    "import intelqs_py as iqs\n",
    "\n",
    "import numpy as np              # Import NumPy library\n",
    "import matplotlib.pyplot as plt # Import graphical library for plots\n",
    "import pandas as pd             # Import library to store results in dataframes\n",
    "import pickle                   # Import library to save dataframes\n",
    "\n",
    "# Name of files and directories\n",
    "data_dir = './noise_via_chi_matrix/'\n",
    "file_df1 = data_dir + 'study_TYPE-channel_qQUBITS.pkl'\n",
    "num_qubits = 4                  # number of qubits"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## [0] Simulate repetead applications of the same channel\n",
    "\n",
    "For all channels, the initial state is $|-\\rangle_3 |+\\rangle_2 |1\\rangle_1 |0\\rangle_0$, with qubits ordered from right to left.\n",
    "\n",
    "### Depolarizing channel\n",
    "\n",
    "The input state $\\rho$ is transformed to the output state $\\rho^\\prime$ according to: \n",
    "$$\\rho^\\prime = (1-p) \\rho + \\frac{p}{3} ( X.\\rho.X + Y.\\rho.Y + Z.\\rho.Z )$$\n",
    "\n",
    "Interpretation: nothing happens with probability $(1-p)$, otherwise each of the Pauli errors happens with equal probability $p/3$. Convergence towards the maximally mixed state happens with (single qubit) rate $\\lambda=4 \\frac{p}{3}$.\n",
    "\n",
    "\n",
    "### Dephasing channel\n",
    "\n",
    "The input state $\\rho$ is transformed to the output state $\\rho^\\prime$ according to: \n",
    "$$\\rho^\\prime = (1-p) \\rho + p Z.\\rho.Z$$\n",
    "\n",
    "Interpretation: nothing happens with probability $(1-p)$, otherwise the Pauli Z error happens with probability $p$.\n",
    "\n",
    "### Amplitude-damping channel\n",
    "\n",
    "The Kraus operators are $M_0=\\begin{bmatrix} 1 & 0 \\\\ 0 & \\sqrt{1-p}\\end{bmatrix}$\n",
    "and $M_1=\\begin{bmatrix} 0 & \\sqrt{p} \\\\ 0 & 0 \\end{bmatrix}$.\n",
    "\n",
    "The decomposition in Pauli matrices leads to:\n",
    "$M_0 = \\frac{1+\\sqrt{1-p}}{2} I + \\frac{1-\\sqrt{1-p}}{2} Z$ and \n",
    "$M_1 = \\frac{\\sqrt{p}}{2} (X  + i Y)$.\n",
    "\n",
    "The result is:\n",
    "$$\\rho^\\prime = M_0 \\rho M_0^\\dagger + M_1 \\rho M_1^\\dagger = \\sum_{i,j} \\chi_{i,j} \\sigma_i \\rho \\sigma_j$$\n",
    "with\n",
    "$$\\chi = \\frac{1}{4} \\begin{bmatrix} (1+\\sqrt{1-p})^2 & 0 & 0 & p \\\\\n",
    "                                    0 & p^2 & -i p^2 & 0 \\\\\n",
    "                                    0 & i p^2 & p^2 & 0 \\\\\n",
    "                                    p & 0 & 0 & (1-\\sqrt{1-p})^2 \\end{bmatrix}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-- state 0 of the ensemble\n",
      "-- state 100 of the ensemble\n",
      "-- state 200 of the ensemble\n",
      "-- state 300 of the ensemble\n",
      "-- state 400 of the ensemble\n",
      "data saved in file  ./noise_via_chi_matrix/example_amplitude-damping-channel_q4.pkl\n"
     ]
    }
   ],
   "source": [
    "assert True and \"data already saved in pandas dataframes\"\n",
    "\n",
    "#channel_type = 'depolarizing'\n",
    "#channel_type = 'dephasing'\n",
    "channel_type = 'amplitude-damping'\n",
    "\n",
    "# Prapare the initial state\n",
    "assert num_qubits == 4                                        # number of qubits\n",
    "index = 0*1 + 1*2 + 0*4 + 1*8\n",
    "psi = iqs.QubitRegister(num_qubits, 'base', index, 0) # |psi> = |1010> = |1>_3 .|0>_2 .|1>_1 .|0>_0\n",
    "psi.ApplyHadamard(2)                                  # |psi> = |1+10>\n",
    "psi.ApplyHadamard(3)                                  # |psi> = |-+10>\n",
    "\n",
    "# Associate a random-number-generator to the state\n",
    "rng = iqs.RandomNumberGenerator();\n",
    "rng_seed = 7777;\n",
    "rng.SetSeedStreamPtrs( rng_seed );\n",
    "\n",
    "# Define the chi-matrix\n",
    "p=0.01\n",
    "chi = iqs.CM4x4()\n",
    "if channel_type == 'depolarizing':\n",
    "    chi[0,0] = 1-p\n",
    "    chi[1,1] = p/3\n",
    "    chi[2,2] = p/3\n",
    "    chi[3,3] = p/3\n",
    "elif channel_type == 'dephasing':\n",
    "    chi[0,0] = 1-p\n",
    "    chi[3,3] = p\n",
    "elif channel_type == 'amplitude-damping':\n",
    "    chi[0,0] = (1+np.sqrt(1-p))**2 /4\n",
    "    chi[0,3] = p/4\n",
    "    chi[3,0] = p/4\n",
    "    chi[3,3] = (1-np.sqrt(1-p))**2 /4\n",
    "    chi[1,1] = p**2 /4\n",
    "    chi[1,2] = -1j * p**2 /4\n",
    "    chi[2,1] = +1j * p**2 /4\n",
    "    chi[2,2] = p**2 /4\n",
    "else:\n",
    "    assert False\n",
    "chi.SolveEigenSystem()\n",
    "chi.Print(True)\n",
    "\n",
    "# parameter of the evolution\n",
    "num_ensemble_states = 500\n",
    "num_time_steps = 200\n",
    "collective_list = []\n",
    "for s in range(num_ensemble_states):\n",
    "    if s%100==0:\n",
    "        print(\"-- state\", s, \"of the ensemble\")\n",
    "    psi_s = iqs.QubitRegister(psi)\n",
    "    psi_s.SetRngPtr(rng);\n",
    "    for t in range(0,num_time_steps):\n",
    "        sq_overlap = np.absolute(psi_s.ComputeOverlap(psi))**2 # the overlap is computed as if at the begin of time step t\n",
    "        if sq_overlap>1:\n",
    "            print(s, t, sq_overlap)\n",
    "            psi_s.Print('large overlap?')\n",
    "        if psi_s.GetOverallSignOfChannels() < 0:\n",
    "            sq_overlap = sq_overlap * (-1)\n",
    "            print('possible error: overall negative sign')\n",
    "        for q in range(num_qubits):\n",
    "            psi_s.ApplyChannel(q, chi)         # then the channel is implemented\n",
    "        collective_list.append([s, t, sq_overlap])\n",
    "df = pd.DataFrame(collective_list, columns =['s', 't', 'sq_overlap'])\n",
    "filename = file_df1.replace('TYPE', channel_type)\n",
    "filename = filename.replace('QUBITS', str(num_qubits))\n",
    "df.to_pickle(filename)\n",
    "print('data saved in file ', filename)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## [1] Plot the results (one type of channel at a time)\n",
    "\n",
    "Depolarizing:  \n",
    "$|0\\rangle\\langle 0| \\rightarrow I$  \n",
    "$|1\\rangle\\langle 1| \\rightarrow I$\n",
    "\n",
    "Dephasing:  \n",
    "$|0\\rangle\\langle 0| \\rightarrow |0\\rangle\\langle 0|$  \n",
    "$|1\\rangle\\langle 1| \\rightarrow |1\\rangle\\langle 1|$\n",
    "\n",
    "Amplitude-damping:  \n",
    "$|0\\rangle\\langle 0| \\rightarrow |0\\rangle\\langle 0|$  \n",
    "$|1\\rangle\\langle 1| \\rightarrow |0\\rangle\\langle 0|$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<bound method NDFrame.head of          s    t  sq_overlap\n",
      "0        0    0    1.000000\n",
      "1        0    1    0.980309\n",
      "2        0    2    0.961031\n",
      "3        0    3    0.942155\n",
      "4        0    4    0.923673\n",
      "...    ...  ...         ...\n",
      "99995  499  195    0.032765\n",
      "99996  499  196    0.032267\n",
      "99997  499  197    0.031777\n",
      "99998  499  198    0.031295\n",
      "99999  499  199    0.030821\n",
      "\n",
      "[100000 rows x 3 columns]>\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA4MAAAGBCAYAAADcwGTYAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAACQ0ElEQVR4nOzdd3xUVfrH8c+Z9E4SCB1CU4pIFQEVgkgRFVZdV10bu+u6q9tQFzsrLuq6trXr4qrgz127IlhAUKIgoGCX3jsBkpDeZub8/phk0mECk0zK9/16jXPLufc+cw05eeace46x1iIiIiIiIiItiyPQAYiIiIiIiEjDUzIoIiIiIiLSAikZFBERERERaYGUDIqIiIiIiLRASgZFRERERERaICWDIiIiIiIiLZCSQRERERERkRZIyaCIiEgzYYwZZYyZb4zZa4yxxpipPhzT3xjzmTGmoPS4vxljTAOEKyIiAaZkUEREpPmIBn4C/gIUHKuwMSYWWAykAacBfwamAzfVY4wiItJIGGttoGMQERERPzPG5AJ/tNbOOUqZ64F/Am2ttQWl2+4Crgc6Wf2RICLSrKllUEREpOUaASwrSwRLLQI6AMkBiUhERBpMcKADqE+tW7e2ycnJJ3SOvLw8oqKi/BNQC6d76T+6l/6je+k/J3Ivv/7668PW2jZ+DkmOrR2wp8q2tAr7th/t4BOtZ/Xvz0P3oZzuhYfug4fuQ7n6qmObdTKYnJzMmjVrTugcqamppKSk+CegFk730n90L/1H99J/TuReGmN2+jcaqYOqXUFNLds9O425DrgOoG3btjz88MPHfeHc3Fyio6OP+/jmQvehnO6Fh+6Dh+5DuRO5F2PGjKm1jm3WyaCIiIgc1QE8LYAVJZW+p1EDa+1sYDbA0KFD7Yl8maIvYzx0H8rpXnjoPnjoPpSrr3uhZwZFRERarpXAWcaY8ArbxgH7gB0BiUhERBqMkkEREZFmwhgTbYwZaIwZiKeO71K63qV0/z+MMZ9UOOR/QD4wxxhzijHmIuA24FGNJCoi0vwpGRQREWk+hgLflr4igHtKl/9eur890KOssLU2C09LYAdgDfA08AjwaMOFLCIigaJnBkVERJoJa20q5QPA1LR/ag3bfgRG1V9UIiLSWKllUEREREREpAVSy6BIFdnZ2Rw8eJCSkpJAh1KruLg41q9fH+gwmgXdS/+p6V6GhISQlJREbGxsgKISERGR2igZFKkgOzubtLQ0OnbsSEREBMbU2tsqoHJycoiJiQl0GM2C7qX/VL2X1loKCgrYu3cvgBJCERGRRqbBu4kaY0YZY+YbY/YaY6wxZqoPx/Q3xnxmjCkoPe5vprH+lS5N2sGDB+nYsSORkZGNNhEUaSqMMURGRtKxY0cOHjwY6HBERESkikA8MxgN/AT8BSg4VmFjTCywGM/kt6cBfwamAzfVY4zSQpWUlBARERHoMESalYiIiEbd7VpERKSlavBuotbaD4EPAYwxc3w45AogErjGWlsA/GSM6QPcZIyp13mQSlxuNmS4GF7iIjwkqL4uI42MWgRF/Ev/pkRERBqnpvDM4AhgWWkiWGYRMAtIBrbXx0XvfX8dr6/eTU6Rk76nZDDqpDb1cRkREZEWKfm2DzwLCz8IbCCNhe5DOd0LD90HD90Hrx0p/j9nU0gG2wF7qmxLq7CvUjJojLkOuA6gbdu2pKamHtdFd+4uIqfICcD/ffIN7n1hx3UeKZebm3vc/z8aSlxcHDk5OYEO45hcLleTiLOqSZMmAfDhhx8CsGzZMs477zw++OADzjrrLACefvppOnfuzOTJk/167fvvv58HHniA7OzsStv9fS9///vfs3z5cn766Se/ndOf/vvf/3L99dfz448/0rVrV7+e+2j3srCwsNH/+xcREWlpmkIyCFC1K6ipZTvW2tnAbIChQ4falJSU47qgo8MhFr/4FQBb88I43vNIudTU1EZ/H9evX98kRpZsqiNgBgV5uluXxX7mmWeycuVK+vbt69323HPPceaZZ3LFFVf49dphYWGVrl3G3/cyJCQEY0yj/f9z8cUXM3DgQHr16uW9J/5ytHsZHh7OoEGD/Ho9EREROTFNIRk8gKcFsKKk0vc06smwbglEhARRUOJi2+E8dqbn0TUxqr4uJ9IixcbGMnz48ECH0aK0adOGNm3U7V0Cb8cD5zWJLwkbgu5DOd0LD90HD92HcvXVuyYQo4nW1UrgLGNMeIVt44B9wI76umh4SBAjeyR611M3HqqvS4nUmy1btnDVVVfRrVs3IiIi6N69O9dffz2ZmZmVyk2dOpVOnTqxZs0aRo4cSUREBCeffDIffODpp//oo4+SnJxMbGwsU6ZM4dChyv8ejDHceeed3HfffXTq1ImIiAhGjRrFd999d9T4UlNTMcZ4f8ElJyezc+dO/vvf/2KMwRjD1KlTvTEmJydXO0dKSkq1iuLbb7/lrLPOIjw8nI4dOzJr1ixqGmvK6XTyyCOP0Lt3b8LCwujQoQM333wzhYWFR427zCeffMLgwYMJDw+nR48e/Pvf/66x3N13383gwYOJi4ujdevWnH322axatarGezFv3jx+97vfkZCQQHx8PDfeeCMul4vVq1dz5plnEhUVRb9+/Vi0aFGl48v+H65YsYLTTjuN8PBwkpOTefLJJyuVmzNnDsYYduzY4d2WnJzMlVdeyWuvvUafPn2Iiopi6NChLF++vNpnefzxx0lOTiY8PJxhw4axYsUKkpOTvf+fREREpOlo8JZBY0w00LN01QF0McYMBDKstbuMMf8Ahllrx5aW+R9wNzDHGHMvcBJwG3BPfY4kCpBychs+2eCZGyt140GuGZlcn5cT8bt9+/bRqVMnHnvsMeLj49m2bRv3338/kyZNYuXKlZXKZmdnc/XVV/PXv/6VDh06cN9993HxxRfzhz/8gU2bNvH000+TlpbGtGnT+MMf/sAbb7xR6fiXX36ZLl268NRTT1FUVMTf/vY3xo4dy+bNm0lISPAp3nfffZdJkyYxYMAAZs6cCVDnVqzDhw9z9tln065dO+bOnUtYWBgPPfQQu3btqlb2yiuvZMGCBdx6662MHDmS9evXM2PGDHbs2MHbb7991OusX7+eSZMmMXToUF577TWKioqYOXMmubm53u6wZfbu3cuNN95Ip06dyMvL45VXXmHUqFGsWbOGU089tVLZadOmcdFFF/H666/z+eefc++99+J0OlmyZAnTp0+nY8eO3HvvvVx00UXs3LmT1q1be4/Nzs7m0ksv5dZbb6Vnz5689tpr/PnPfyYmJuaYydqyZcvYuHEjs2bNIjw8nBkzZnD++eezY8cOWrVqBcB//vMfpk2bxm9+8xsuueQStm7dyi9/+UuOHDly1HOLiIhI4xSIbqJDgaUV1u8pfc0FpgLtgR5lO621WcaYccDTwBogE3gEeLS+A005OQlYC8CKrekUaoqJFsk74l0jsOOB8+pUftSoUYwaNcq7PnLkSHr27MlZZ53Ft99+W+kZrpycHJ577jlv+Q4dOjBgwADef/991q1b501wfvrpJ5588klcLlelpKegoICPP/6YqChPd+rTTz+dXr168a9//YtZs2b5FO+gQYMICwujdevWx9199F//+hd5eXksWrSILl26ADBu3Lhqg6UsW7aM119/neeee47f/e53AJxzzjkkJCRw5ZVX8t133zFw4MBar3PvvfcSExNT6TOPHDmSHj160KFDh0pl//Of/3iXXS4XEydOpF+/frzwwgs8/vjjlcqeffbZPProo964P/jgA5566imWLVvGmWeeCUD79u0ZMGAAH3zwAddcc4332JycHGbPns1ll10GwMSJE9m7dy93330311xzzVGneMjOzua7774jPj4egHbt2nHaaafx4Ycf8stf/hK3280999zDueeeW+nztGvXjosvvrjW84qIiEjj1eDdRK21qdZaU8Nraun+qdba5CrH/GitHWWtDbfWtrfW1nurIEDnhEjaR3n+eCpyulm1Lb2+LyniV8XFxdx///307t2biIgIQkJCvKN2bty4sVLZqKioSolj7969AU+CVDHp6927N06nk/3791c6ftKkSd6kCDxdD4cPH16tBbK+rVy5kuHDh3sTQfB8tgsuuKBSuYULFxIaGsqUKVNwOp3e1/jx4wH4/PPPAU/yVnF/2a+elStXVvvMnTt35owzzqgW05IlSxgzZgyJiYkEBwcTEhLCpk2bqv0/ADj33HMrrffu3ZuoqChvIli2DWD37t2VygYFBVVLzC677DJ27drF3r17a7ljHiNGjPAmggD9+/cH8Lao7tmzhz179nDJJZdUOm7KlCkEBzeFx89FRESkqqbwzGBAndq6/I9gPTcoTc3tt9/OzJkzufLKK/nggw/46quveOeddwCqPRdX1hWwTGhoKEClBKHi9qrHt23bttr127Zte8wkxN/2799faywVHTx4kOLiYtq3b09ISIj3lZTkGZ8qPd3z5U+PHj0q7Z87d26drvPNN98wadIkoqOjeeGFF1i1ahWrV69mwIABNT6bWNP9ru3/TdXj4+PjCQkJqTGeY/1/qNqVt2yk0bJrlCX/ZfenTFBQUKWuqiIiItJ06OvcYzi1TTCLdnrmG0zdeBDoF9iApMHVtWtmY/Laa69x9dVXc9ddd3m35ebm1su10tKqD+6blpZGx44d/XL+8PBwiouLq21PT08nMbF8sKf27dvXGktFiYmJhIeHs3Dhwkqte2XKunouWLCAoqIi7/Zu3brV6Tpvv/02wcHBvPPOO5UStczMzGpJ3onKzMykpKSk0nXK4jnR/w/t27cHPEl0RS6Xi8OHD5/QuUVERCQw1DJ4DCclOIgM9bQO7kjPZ/vhvABHJOK7/Pz8ai1FL730Ur1c68MPPyQvr/zfx44dO1i1ahUjRoyo03nCwsIoKCiotr1r166kpaVVSjy2bt1aravliBEjWLVqVaUulHl5eSxYsKBSuYkTJ1JYWEh2djZDhw6t9ipLBvv3719pe1niOWLEiGqfeffu3XzxxReVrpOfn09QUFCl5/U+/fTTGge0OVEul6vawDevvfYaXbp0OeFksFOnTnTq1Ik333yz0vZ58+bhdDpP6NwiIiISGEoGjyHEYRjZo7wLlKd1UKRpmDhxInPnzuWZZ57h448/5ve//z0rVqyol2tFREQwfvx45s2bx+uvv87EiROJjY3lxhtvrNN5+vbty7Jly3j//fdZs2aNdwqESy65BGMMV1xxBYsWLeK///0vU6ZMqdZF8cYbbyQqKorx48fz+uuvM2/ePMaPH09ERESlcikpKVx++eVcffXVzJo1i0WLFrF48WKef/55LrzwQjZt2nTUOO+66y6ys7O9n/mNN95g/Pjx1bqJTpw4kdzcXKZOnconn3zCs88+y5VXXum3FtOKYmJiuOWWW3jqqadYtGgRU6dOZcmSJfz9738/6uAxvnA4HNx999189NFHXHvttSxatIhnn32Wm266ibi4OBwOVSciIiJNjWpvH6ScXD60/VI9NyhNyJNPPsnkyZO58847ufTSS8nJyeHVV1+tl2tdffXVnHfeefzxj3/kmmuuoU2bNnzyySc+TytR5h//+Acnn3wyv/jFLzjttNO8U0z07NmTt956i7179/Kzn/2MBx98kEcffZSTTjqp0vGtW7fmk08+oXXr1lxzzTX84Q9/YOLEifz617+udq1XXnmF2267jbfeeospU6bw85//nKeeeopevXrV+DxgRX369OHDDz8kPz+fSy+9lNtuu41p06YxduzYSuUmTJjAE088wRdffMH555/Piy++yMsvv0zPnj1rOfPxi42N5bXXXmPu3LlMmTKFpUuX8vjjj1cacfREXHvttfzrX/9i8eLFTJkyhRdeeME7J2RcXJxfriEiIiINxzTAoJwBM3ToULtmzZoTOkdqaio9BwzjzH96ZsMIDXbw/d/GExGqKSbqKjU1tdrk4I3N+vXr6dOnT6DDOKacnBxiYmICHYZX2aTz9957b6BDqbPGdi+PV1kr4J49exr0uqtXr2bYsGG8/PLL/OxnP6v1Xh7r35Yx5mtr7dD6ilPqx4nWs02hXmgIug/ldC88dB88dB/Knci9OFodqwFkfNApPpJeSdFsPphLsdPNym2HObv30VsNRESam+3bt/P0009z1llnERsby/r167n//vvp1q0bF198MS6XK9AhioiISB0oGfTRmN5JbD7oGYUxdeMhJYMi0uJERETw008/8fLLL5OZmUl8fDznnHMODzzwAJGRkeTk5AQ6RBEREakDJYM+SjmpDbM/3wZ4kkFr7QkPyCDSXDTn7uZNxZw5c+r9Gu3atWPhwoX1fh0RERFpGBpAxkdDkxOIKn1OcFdGPts0xYSIiIiIiDRhSgZ9FBrs4IyeFaeY0KiiIiIiIiLSdCkZrIOUk5O8y5pvsPlSl0cR/9K/KRERkcZJyWAdVJxv8MttGeQXOwMYjdSHkJAQCgoKAh2GSLNSUFBASEhIoMMQERGRKpQM1kGHVhGc3NYzh1axy82KLekBjkj8LSkpib1795Kfn6/WDJETZK0lPz+fvXv3kpSUdOwDREREpEFpNNE6Sundho1pnuHTUzcd5Jy+mmKiOYmNjQVg3759lJSUBDia2hUWFhIeHh7oMJoF3Uv/qelehoSE0LZtW++/LREREWk8lAzWUcpJSfz7M88UE0s3aIqJ5ig2NrbR/+GamprKoEGDAh1Gs6B76T+6lyIiIk2LuonW0dDkeKLDPDn03iMFbErLDXBEIiIiIiIidadksI5CghyMPql8IJlFaw8EMBoREREREZHjo2TwOIzvV/6coJJBERERERFpipQMHoezeycRGuS5dWv3ZbM7Iz/AEYmIiIiIiNSNksHjEBMewsieid71j9elBTAaERERERGRulMyeJwm9GvnXVZXURERERERaWqUDB6nc/q0pWxGiTU7MkjPLQpsQCIiIiIiInWgZPA4tYkJY2jXeADcFpasV1dRERERERFpOpQMnoDKXUWVDIqIiIiISNOhZPAEVEwGl28+TE5hSQCjERERERER8Z2SwRPQOSGSvu1jASh2uUndeCjAEYmIiIiIiPhGyeAJ0qiiIiIiIiLSFCkZPEETTmnrXU7deIgipyuA0YiIiIiIiPhGyeAJOrltDF0TIwHILXKyYkt6gCMSERERERE5NiWDJ8gYo66iIiIiIiLS5CgZ9IMJ/cq7ii5el4bLbQMYjYiIiIiIyLEpGfSDQZ3jaRMTBkB6XjFfbc8IcEQiIiIiIiJHp2TQDxwOw8QKXUU//HF/AKMRERERERE5NiWDfnLeqe29yx/9dEBdRUVEREREpFFTMugnpyUn0Dra01X0cG6RuoqKiIiIiEijpmTQT4IchnNPUVdRERERERFpGpQM+tGk/uoqKiIiIiIiTYOSQT8a1q1yV9HVO9RVVEREREREGiclg35UtavoBz+oq6iIiIiIiDROSgb9TF1FRURERESkKQgOdADNTVlX0cO5Rd6uosO7JwY6LBERkcYjcycc2khC+g+wqTjQ0QSc7kO5Jnsv4jpC236BjkKkzpQM+lmQwzDxlLa8smoX4BlVVMmgiIhIBZsWwUfTORXgx0AHE3i6D+Wa9L0Y8isYcBlYC9jK79ZdfRsWLDWWSTz8I6zPPcZ5qOF87pqvX2OZmo6vy3mOsa9inFS5nneZKuWoVK7Hnt1QsLCGckc7h6/XshXefD3fiV6rwues47USIoYDKfibksF6cF7/DhWSwQPcfUE/ghwmwFGJiIiISL35+iXPyw/6A/zkl1M1aZ0B9gQ6isYhvFePejmvksF6oK6iIiIiR9GqM/QcR3pGBokJCYGOJuB0H8o1uXvhLoFtqYGOQuS4KRmsB+oqKiIichQnnwsnn8uPqamkpKQEOpqA030o1yTvRdo6WHI3FGQCBoyp8u4oXaaGfVXLeLYdTk+ndes2FcpQrUy1d+OoZR8+lKntPLXE7Mt5Kl7bex4qn5sq16hSbsvWrfTs2fOY5SqvH9+1fC5X6z5fyx3ftdK3HaE+KBmsJ5P6t1dXURERaXDGmBuA6UB7YC0wzVq77CjlJwAzgVOAIuALYLq1dlP9RyvSDLTtC1e86ddT/tQUk+J6sKc4lZ4jUgIdRqNQtC+1Xs6rqSXqyendEmkdHQpoAnoREWkYxphLgceB+4FBwArgI2NMl1rKdwPeA5aVlj8HiAA+bJCARUQkoAKSDBpjbjDGbDfGFBpjvjbGnHWM8hOMMSuNMTnGmMPGmPeMMSc1VLzHw9NVtHwC+g9/1AT0IiJS724C5lhrn7fWrrfW/gnYD1xfS/khQAhwu7V2i7X2O+AfQA9jTOsGiVhERAKmwZPBlvStpSagFxGRhmKMCcWT3H1cZdfHwMhaDlsDlADXGmOCjDExwDXAamvt4XoLVkREGoVAtAy2mG8tK3YVPZSjrqIiIlKvWgNBQFqV7WlAu+rFwVq7AxgH3IPnecEsPKPan19vUYqISKPRoAPIVPjW8uEqu3z91vI/QCRN5FvLsq6iGlVUREQaUNVuKKaGbZ4dxrQDXgBeBl4FYoC/A28YY8621rprOOY64DqAtm3bkpqaetyB5ubmntDxzYXuQzndCw/dBw/dh3L1dS8aejTRo31reU5NB1hrdxhjxgFvAk/jac38Fji3HuP0m4qjin7ww35mnN+XkCCN2yMiIn53GHBRvRUwier1bpk/AHnW2lvKNhhjrgR24/mSdnnVA6y1s4HZAEOHDrUnMuJhqkZMBHQfKtK98NB98NB9KFdf9yJQU0vU27eW/vzGEk48C3dbS3yYIbPIkp5XzFNvf8qgpJY5o4e+3fEf3Uv/0b30H93LwLLWFhtjvsbT7bPiOPfjgLdrOSwSTwJZUdm6vrkUEWnmGjorqfdvLf35jSX4Jwu/tHADz322FYDNJfHcmDLkhM7XVOnbHf/RvfQf3Uv/0b1sFB4F/s8Y8xWe+QJ/D3QAngMwxvwDGGatHVta/gPgRmPM3cD/8Hzhej+eOvbrBo5dREQaWIN+62etLcZTuYyrsmscnlFFa9Lkv7W8eHBH7/KSdQfJyi8JYDQiItJcWWtfB6YBdwHfAWcCk6y1O0uLtAd6VCj/KfBLYAqeRzAW4XlOf6K1Nq/BAhcRkYAIRDL1KDDVGHOtMaaPMeZxqnxraYz5pEL5D4DBxpi7jTG9jDGDgZdoQt9a9mobw6md4gAodrl5/8d9AY5IRESaK2vtM9baZGttmLV2iLX28wr7plprk6uUf81aO9haG22tbWOtvcBau67BAxcRkQbX4MlgS/3W8qJB5a2Db3+9J4CRiIiIiIiIBKibZUv81vKCAR0IdhgAvtl1hO2Hm0weKyIiIiIizVCTeOauOUiMDiPl5CTv+rvfqHVQREREREQCR8lgA/r5kPKuou98uxe3u8bZNEREREREROqdksEGNKZ3EnERIQDsySxg9Y6MAEckIiIiIiItlZLBBhQWHMQFA9p719/5Zm8AoxERERERkZZMyWADu2hwJ+/yBz/up6C46hSKIiIiIiIi9U/JYAMb1LkV3VpHAZBb5OTjdQcCHJGIiIiIiLRESgYbmDGm8pyD6ioqIiIiIiIBoGQwAC4c3BHjmXKQZZsPse9IQWADEhERERGRFkfJYAB0io/kjB6tAbAW3v5acw6KiIiIiEjDUjIYIJcMLR9I5s2v92jOQRERERERaVBKBgNkQr92xIYHA7ArI58vt2vOQRERERERaThKBgMkPCSIyQM7eNff/VZdRUVEREREpOEoGQygCweVdxX96McDFJZozkEREREREWkYSgYDaHCXVnRJiAQgp8jJ0g0HAxyRiIiIiIi0FEoGA8gYw5QKXUXnfac5B0VEREREpGEoGQywKQPLJ6BfuuEQmXnFAYxGRERERERaCiWDAdYzKZoBneIAKHa5eUtzDoqIiIiISANQMtgI/PL0Lt7l/365U3MOioiIiIhIvVMy2AhcMKADMaVzDu5Iz2fF1vQARyQiIiIiIs2dz8mgMaaVMeYeY8zHxpi1pe8zjTGt6jG+FiEyNJiLB5dPM/HKqp0BjEZERBqa6lgREQkEn5JBY8wAYDNwOxAOrCt9vwPYZIzpX28RthBXVOgqunh9GgeyCgMYjYiINBTVsSIiEii+tgw+AaQDvay1o6y1l1hrRwEnARnAk/UVYEvRq20Mp3dLAMDltry+eneAIxIRkQaiOlZERALC12TwNGCGtbZS/0Vr7Q7gbmCYn+Nqka4Y3tW7/OpXu3C63AGMRkREGojqWBERCQhfk8F0oKiWfYWl++UETezXjtbRoQAcyC7kkw0HAxyRiIg0ANWxIiISEL4mg88C040x4RU3GmMigL8CT/s7sJYoNNjBL4Z29q7/98tdAYxGREQaiOpYEREJiGAfy0UCXYFdxpgPgTSgLTAJKACijDF/Ly1rrbV3+z3SFuLyYV149rOtWAufbzrEzvQ8uiZGBTosERGpP6pjRUQkIHxNBu+osHx1DfvvrLBs8TzjIMehc0IkKSe1YenGQwD878td3D6pT4CjEhGReqQ6VkREAsKnbqLWWkcdXkH1HXRzd2WFgWTeWLObwhJXAKMREZH6pDpWREQCxedJ56XhpJycRMdWEQBk5pcw/7t9AY5IRERERESaGyWDjVCQw3D1iPLWwRe/2I61NoARiYiIiIhIc+NzMmiMuc4Y860xJt8Y46r6qs8gW6LLTutCRIinN9CGAzms3KaRxUVEmivVsSIiEgg+JYPGmKuBJ4HVQDjwEvAKkA1sBf5e+9FyPOIiQ7h4SEfv+ktf7AhcMCIiUm9Ux4qISKD42jI4DfgHcH3p+jPW2muA7niGvVazVT2YOrKbd3nJ+jR2pucFMBoREakn01AdKyIiAeBrMtgL+Bxwl75CAay1mcB9wF/qJboWrmdSNKNPagOAtTB3xc4ARyQiIvVAdayIiASEr8lgAeCwnlFMDuD5trJMLtDB34GJx6/OSPYuv7FmNzmFJYELRkRE6kOLq2M/2fkJv170a5448AS/XvRr/vHlP8guzg50WCIiLY6vk87/CPQElgDLgDuMMdsBJzAT2FAv0QmjerWhe5soth3KI7fIyVtf7+FXZ3Q79oEiItJUtLg6Ni0/jdUHVntWDsDqA6tpH9WeqadMDWhcItK0WGtxWzdu3OXL1u3Zh2fdYr2j8ldcr/juPcZa77mqlnPjBlv9vLW9l52nUixVy5Vdu8p5vbFWKJ9Zklkv99DXZHA25d9UzsBTYS0vXc8BfubfsKSMw2H41RndmDHvJwDmrNjBNSOScThMgCMTERE/UR0L7MvTnLpNmbUWl3VR4i7xvFwl5cul6063s9I2p9uJy+3CZV04rWf5p9yfyNyc6V13WZennHXhcrsqba+67mu5sm1liYPLuiq9+7q94nrriNbcdtptnJRwUnkZt9t7vUrnc7tqPH/F95/yfiJna06tsVQ8R6VrHOVzuayrUsJUlnyUJR1uKixXSHaqbq+aeFVMdKrtr6Xs0a5ZcX+JswTHfx21HteSXJpwab2c16dk0Fr7eoXlLcaYfsAIIBJYYa09XC/RCQAXD+7IQws3kF3oZGd6Pp9uOMg5fdsGOiwREfGDlljHnt3lbHq26smTy5/ku/zvAHC5NYPGiXJbNwXOAvJL8sl35lPgLKDQWUiRq6jyy1m+XOgqpNhVTKGz9N1V/l5WrmoSV1OiV+Iu8bZonLAV/jlNQzqQd4BpqdP8e9Llxy7SIjgDHUDj4Ld/X1X42jJYibU2D883l9IAIkODuXxYF/79+TbAMwm9kkERkeapJdSx7aLa0S6qHb3De5cng7blJoNu6yarKIvsomyyS7LJLsompziH7OLy99ziXG+CV5bslb3nleRR4CygwFkQ6I8i0uAcxoEDB8YYDAaHKV/2vpctVygDlJfHgMG7XOlYU6EMFcqYyuf3nrNiHGXnpYaYjMGBA0+RY5dPKEqol/tXazJojOlSlxNZa3edeDhSm6tGdOX5ZdtwW1ixNZ21+7Lo1yEu0GGJiMhxUB3rUfbHEzSvZNBt3aQXpJNemE5GQYbnvTDDu55R6HllFmZ6Er2SXGgG/4eDTBAhjhDPKyiEYBNMSJBnPdgR7N1eVibIEUSI8bwHmSCCHEEcPniYju07EmSCCHYEe7cHm+BK5Squ11qudFvVMmX7HMaBwziqLzs8yUXZcsX9FcuXvX+d9jUPr3mY3JLcavscxkGwI7jG47xlHNW3pR9Kp13bduXbHTUfW3H/0fZVjaks+ShLnBx4lssSkarbyxKisiTlWOepmEAdbX9t5ytLglZ8sYJRZ4066nlaitTU1Ho579FaBndAndojg04sFDmaTvGRnHtKez74cT8Az3++jccuGxTgqERE5DjtQHWs51vxUk2pm2iJq4TdubvZnb2b/Xn7OZB3gAP5B9ifu5+0/DTS8tNwugPTty0iOILI4EgiQyIJDw4nPCicsKAwwoLDCHN43sODwgkNCvXsCw7z7A+qsD043LstNCiU0KDQ8iSvSkJXMfELcpz4j2lqaiopZ6Sc+I1oQGd0PIMzOp7h13OmpqaSMirFr+dsiiIcEUSFRAU6jGbtaMngr6lbRSX17LpR3b3J4IIf9nPLxN50aBUR4KhEROQ4qI6lSjLYCFsG80vy2ZS5iU2Zm9iRvYMdWTvYmb2Tvbl7/RqvwRAdGk1saCyxobHEhMYQExrjXY4NjSU6NNqb5NX2Hh4U7peETERajlqTQWvtnLJlY0wcUGitLWqIoKRmAzq34vRuCXy5PQOX2/Li8u3cdX7fQIclIiJ1pDrWwx/dRHdl7+L1ja+zJ2cP1w+8nt4JvY/rPCXuEtanr+fbg9+yNn0tGzI2sCNrx3EP2tAqrBWtI1qTGJFIQngCieGl76XrCeEJxIfHExcWx5ov1nD2mLOP6zoiIifimAPIGGOCgXTgQmBBvUckR/W70d35cnsGAK9+tYs/je1FXERIgKMSEZHj0dLr2Iotg3UdJn539m6e/f5ZFmwrv215zjz+M/4/Ph3vtm7Wp6/n872f8/WBr/nh8A91GoClfVR7usZ2pUN0B8+AOJHtvAPjtItqR0Sw7z13KibFIiIN6ZjJoLXWaYxJAxpf/40WKOWkJHomRbPlYC55xS7+9+Uurk/pEeiwRETkOLT0OrZiEuTrM3YlrhKe//F5nv/x+WrH7Mo++igsTreTL/d/yaIdi/h8z+ekF6YfM75usd3ondibHnE96BrbleS4ZLrEdCE8ONyneEVEGjNfp5Z4BbgW+LAeYxEfOByG687qzi1v/wDAS19s5zdndiM0WN8qiog0US22jq1ry+D2rO3c/NnNbM7cXOP+9IJ0rLXVRhjcemQr725+lw+2f8DhgtqnbewY3ZHBSYMZ0GYAfRL70Cu+V51a+EREmhpfk8EdwC+NMauB94D9VHnw3Vr7on9Dk9pMGdSBhz/eyMGcIg7mFPHed3u5ZGjnQIclIiLHZwcttI6t1DJoj94yuGr/Km5aehM5JTnebae2OZW/Dv0rv1v8OwqcBRS7i8ktySUmNAZrLV8d+Io5a+ewfG/Ns3fHh8VzVqezGNlhJEPaDqFdVDv/fDARkSbC12Tw6dL3jsCQGvZboFlWVI1RWHAQU89I5sGFGwF4ftk2fj6kU4uaa0VEpBlpsXWsr1NLrNi3gj9/+meKXJ4xdsKCwpg2eBqX976cIEcQCeEJ7M3dC0BGYQY7s3fy0OqH+ObgN9XOlRieyKTukxjfdTz9W/fX6Jsi0qL5mgx2q9copM6uOL0rT3+6hbxiF5vSckndeIgxvZMCHZaIiNRdi61jfekmuiFjA3/59C/eRDApMoknz36Svonlo2knhid6k8HrPr6OfXn7Kp3DYDi7y9lc3OtiRnQYQbDD1z9/RESaN59+G1prd/rzosaYG4DpQHtgLTDNWrvsKOUN8Bfg93gqzQxgrrX2Nn/G1ZTERYRw2bAuvLB8OwCzP9+mZFBEpAnydx3blBxrAJmsoiymLZ1GoasQ8Izg+cKEF+gcU/nRiISIBO9yxUQw2BHMRT0v4pp+19Altou/wxcRafLq9NWYMeZUYBSQCPzbWnvAGNMTSLPW5hz9aO85LgUeB24Alpe+f2SM6WutrW0YsEeA8/EkkD8CcXgSyRbtV2ckM2fFDlxuy8pt6fy4J4v+neICHZaIiBwHf9SxTU0Q5V00a2oZfGj1Q94Wv6iQKJ4b91y1RBA8LYNVndPlHG4ccqOSQBGRo/ApGTTGhOEZ7ewiwOB5fmEBcAB4ENgE+NpKdxMwx1r7fOn6n4wxE4HrgdtruPbJwJ+AU6216yvs+tbH6zVbneIjOf/U9rz3nedb0H9/vpWnfjk4wFGJiEhd+LmObVIqPuteddL51QdW897W97zr955xL93jutd4nl7xvbzLSZFJzBo5i5EdR/o5WhGR5sfX+QjuA84BrgLa4qmsynwETPDlJMaYUDwPx39cZdfHQG2/tacA24CJxphtxpgdxpi5xhj1iQSuG1VeMX744352Z+QHMBoRETkOfqljm6KKLYMVk0FrLQ+tfsi7Pq7rOM7pek6t57mo10Vcd+p13DDgBt6d8q4SQRERH/naTfRy4C5r7f+MMVWH3doOJPt4ntZAEJBWZXsanoqwJt2BrsBlwFQ835g+DCwwxoywtnK/EmPMdcB1AG3btiU1NdXH0GqWm5t7wueob/0SHaxNd+O2MPO1ZVzVNyzQIdWoKdzLpkL30n90L/1H9/K4+auObXIqtQxWGE105b6VrM/wdAYKDwrn1tNuPep5IoIj+NOgP9VPkCIizZivyWAisL6WfQ6grtmHrbJuathW9fxXWWs3ARhjrgI2AqcBX1Y6sbWzgdkAQ4cOtSkpKXUMrbLU1FRO9Bz1zdHhEFe/+BUAy/a5ue+K4bSLCw9wVNU1hXvZVOhe+o/upf/oXh43f9exTUZtLYMv/lQ+k8aFvS6kbVTbBo1LRKSl8LWb6HZgRC37huFJzHxxGHABVWd1TaJ6a2GZ/YCzLBEstRlwAnoqHDirV2sGdm4FQLHTzTOpWwIbkIiI1IW/6tgmx1D9mcHtWdv58oDne94gE8Q1/a4JSGwiIi2Br8ngy8BtxpgrgNDSbdYYMwa4ER8nw7XWFgNfA+Oq7BoHrKjlsC+AYGNMjwrbuuNp1Wyxw3FXZIxh2jnlD8+/9tVu9h0pCGBEIiJSB36pY5uioAq9Ysu6iX68o3xYgdGdRtMxumODxyUi0lL4mgw+CHwA/B+eOf7AMy3EEmChtfbJOlzzUWCqMeZaY0wfY8zjQAfgOQBjzD+MMZ9UKL8E+AZ40RgzyBgzCE/F+CWwpg7XbdZGn9SGQV1aAVDsUuugiEgT4s86tkmpOOl8WcvgxzvLk8EJyc127BwRkUbBp2TQWuuy1l4GjMYz599/gCeAs621V9Tlgtba14FpwF3Ad8CZwKQKk+62B3pUKO/GM8fgQeBzYBGwB5hSdfCYlswYw43nnORdf331bvaqdVBEpNHzZx0LYIy5wRiz3RhTaIz52hhz1jHKG2PMNGPMBmNMkTFmvzHmgeP6MHVUcdJ5l3WxPWs7mzI9T4WEOkIZ3Xl0Q4QhItJi1WnSeWvtMmDZiV7UWvsM8Ewt+6bWsG0/cMmJXre5O6tXa4Z0jefrnZmUuCxPL93C/Rf2D3RYIiLiA3/UscaYS4HHgRvwtC7eAHxkjOlrrd1Vy2GP4PnSdTrwIxCH54vZelexZdDtdrN091Lv+pkdzyQqJKohwhARabF8ahk0xnxT+q2hhvNqxKq2Dr65Zjd7MjXvoIhIY+bnOvYmYI619nlr7Xpr7Z/wDMR2fS3XPhn4E57eNu9Za7dZa7+11n7oh1iOqWLLoNM6WX1gtXddrYIiIvXP12cG0/A807DbGPOhMeYyY0zjm7tAOKNnIqclxwN4WwdFRKRR80sda4wJBYYAH1fZ9TFQ2yzsU4BtwERjzDZjzA5jzFxjTFJdr388KrYMlrhK+CbtG+/6ae1Oa4gQRERaNJ+6iVprzy2tGH4JXAn8D8gxxrwFvGKtXXrUE0iDKWsd/OV/PMNyv7lmDzek9KRzQmSAIxMRkZr4sY5tDQRRfaqmNOCcWo7pDnQFLgOm4pnz92FggTFmRE3P5htjrgOuA2jbti2pqak+hlddQV75s+05JTne5VZBrdiyZgtbzdbjPndTkpube0L3sTnRvfDQffDQfShXX/fC52cGrbUHgceAx4wxfYCr8FRcU40xe6y1Xf0enRyXET0SGZacwFc7MnC6LU98spmHLhkQ6LBERKQWfq5jbZV1U8O2MmWT2l9VNp+vMeYqPHMbnoZn5O6qsc4GZgMMHTrUpqSk1CG0yhZ+uhAyq28f2nEoY8aMOe7zNjWpqamcyH1sTnQvPHQfPHQfytXXvfC1m2gl1tr1wN+BO4F9QCd/BiUnxhjDjePKnx18+5s9bE7LOcoRIiLSWJxAHXsYcAHtqmxPonprYZn9gLMsESy1GXACXXyN+Xg5avkzpG9i3/q+tIiIcBzJoDHmbGPMS3gqlpfxTPPwJ38HJidmRI9EzurVGgC3hYcWbQxwRCIiciwnUsdaa4uBr4FxVXaNA1bUctgXQLAxpkeFbd3x9BzaWfMh/lNxAJmKlAyKiDQMX0cTPcUY84AxZhewGM9cSI8Dva21I0qnipBG5taJvb3LH69L4+udNfTFERGRgPJzHfsonq6l1xpj+hhjHgc6AM+VXusfxphPKpRfAnwDvGiMGWSMGQS8iKd76JoT/3RHV1vLYK/4XvV9aRERwfdnBn8AsoA3gf8rnQtJGrlTOsZxwYAOLPh+HwD/XLiB168bjjEmwJGJiEgFfqtjrbWvG2MSgbvwzBX4EzDJWlvWytce6FGhvNsYcz6eSe4/BwrwJKQ31TR4jL8ZqtdHUSFRtIloU9+XFhERfE8GLwXmW2uL6jMY8b+bx53ERz/ux+m2fLU9g9SNhxjTu0FGDBcREd/4tY4tbUmssTXRWju1hm37gUv8ce26MsYQbIJxWqd3W3Jssr60FBFpID51E7XWvllWSRljoo0xnY0xUfUbmvhDcusoLh9WPgbAPxduwOWubVA5ERFpaC29jq363GC3uG4BikREpOXxeQAZY8wEY8wa4AiwA8gyxnxljKn6oLo0Mn8a25PI0CAANhzI4b3v9gY4IhERqagl17FBjqBK651jOgcoEhGRlsfXAWQmAB8A0cAs4AbgXiAG+LAlVFZNWVJMONeeWf5N6yMfb6LI6QpgRCIiUqal17FBpnIy2D6qfYAiERFpeXxtGZwJfAz0tdbeY639t7V2JtAPz4Pm99RPeOIvvx3VnYSoUAD2HinglVW7AhyRiIiUmkkLrmOrtgy2i6o6TaKIiNQXX5PBAcDTVUcWK11/Bhjo57jEz2LCQ/jDmJ7e9Sc+2cyR/OIARiQiIqVadB2rlkERkcDxNRksAmJr2RdTul8auSuHd6FrYiQAWQUlPLZkc4AjEhERVMdWopZBEZGG42symArMMsZUGuLLGNMFT/eWpf4NS+pDWHAQt5/bx7v+f6t2suVgTgAjEhERWngdm1WUVWk9PDg8QJGIiLQ8viaDtwJxwEZjzOfGmNeNMZ8Bm4FWpfulCZjQry3DuycA4HJb7vtgfYAjEhFp8Vp0Heuy5QOaRQRHBDASEZGWx9d5BjcBpwJPAGHAYCAceBwYaK1Vf8MmwhjDjPP7Ujaf79KNh/hs06HABiUi0oKpji0XG1pbb1kREakPwb4WtNbuB/5aj7FIA+nXIY5fDOnM62t2A3Dv++s44y9nERzk87STIiLiR6pjPeLC4gIdgohIi6K//luomyecRFTpRPSbD+by6leaakJERAJLyaCISMNSMthCJcWEc0OFqSYeXbyJrPySAEYkIiItXVyokkERkYakZLAF+82Z3ejYyvOwfmZ+CY99sinAEYmISEsWG6ZnBkVEGpKSwRYsPCSIOyaVTzXx8sqdrN+fHcCIRESkJVPLoIhIw1Iy2MJN6t+OEd0TAc9UEzPm/YTbbQMclYiItEQnJ5wc6BBERFoUJYMtnDGGWT/rR0iQZ66JNTszefubPQGOSkREWopHUx4lKiSKkR1Gcm63cwMdjohIi1Lr1BLGmBfrcB5rrf2NH+KRAOiZFMNvzuzOc59tBeCBjzYwvm874iJDAhyZiEjzpDq23Liu4xjTeQzBDp9nuxIRET852m/es4GK/QVbAXGAE0gHEkuPzwIy6yk+aSB/HtuT+d/tZV9WIel5xTz08Qbu/Vn/QIclItJcqY6tQImgiEhg1NpN1FqbbK3tZq3tBlwF5AKXARHW2vZABHB56fYrGyJYqT+RocH87YK+3vX/frmLH/YcCVxAIiLNmOpYERFpDHx9ZvBR4B/W2jestS4Aa63LWvs68ADwWD3FJw1oQr92jD6pDQDWwox5P+HSYDIiIvVNdayIiASEr8lgf2BLLfs2A6f4JxwJJGMM90zuR2iw58fi+z1ZvPrVrgBHJSLS7KmOFRGRgPA1GTwA/KKWfZcBaf4JRwItuXUUvx/dw7v+z482cCCrMIARiYg0e6pjRUQkIHxNBh8Dfm2M+cAYM9UYc27p+4fANXi6uEgzcUNKD7q1jgIgp8jJjPd+wlp1FxURqSePoTpWREQCwKfhu6y1jxtjcoG7gYqTAO0GfmutrcsQ2dLIhYcE8Y+L+nPZ7FUALF6Xxkc/HWBS//YBjkxEpPlRHSsiIoHi86Tz1toXgK6lr+Gl78mqpJqn4d0TuXxYZ+/6395by5H84gBGJCLSfKmOFRGRQPA5GQTPrLfW2t3W2q9K39V3sBm77dw+JMWEAXA4t4j7Plgf4IhERJov1bEiItLQfE4GjTGDjDHvGGMOG2OcxpjBpdvvN8ZMrL8QJVDiIkKY9bPyQeze/HoPyzcfDmBEIiLNk+pYEREJBJ+SQWPMmcBKoDfwvyrHuYHf+z80aQwm9GvHpP7tvOt3vPsjBcWuAEYkItK8qI4VEZFA8bVl8AFgEdAPuKnKvm+Awf4MShqXmZP7ERvuGWtoV0Y+Dy3aGOCIRESaFdWxIiISEL4mg4OBZ0ufX6j6DMNhoI1fo5JGJSkmnLvO6+tdf2nFdlZtSw9gRCIizYrqWBERCQhfk8FCILKWfe2BLP+EI43VJUM7Mfokz98j1sJf3/ye3CJngKMSEWkWVMeKiEhA+JoMLgemGWOCKmwr+/byN8Cnfo1KGh1jDP+8+FRvd9E9mQXc+/66AEclItIsqI4VEZGA8DUZnIGnG8v3pcsWuMYYsxTPfEj31E940pi0iwuvNLroa6t38+mGtABGJCLSLKiOFRGRgPApGbTWfg+cBaQBdwIG+GPp7tHWWo0o0kJMHtCh0uiit779I5l5moxeROR4qY4VEZFAOWYyaIwJMcZMAY5Ya8cCMUAnINZaO8Za+219BymNhzGGe3/Wn9bRnsnoD+UUMeO9nwIclYhI06Q6VkREAumYyaC1tgR4A0guXS+01u6z1ubXc2zSSCVEhfLARf296+//sJ/53+8LYEQiIk2T6lgREQkkX58Z3AYk1Wcg0rSc07ctlwzp5F2/890f2Z2hv11ERI6D6lgREQkIX5PBB4E7jTGa60i8/nZBXzonRACQU+jkL699i9PlDnBUIiJNjupYEREJiGAfy50NJADbjTGrgP1UnhjXWmuv8Xdw0rjFhIfwxGWDuOS5lTjdlm92HeHxTzZz8/iTAx2aiEhTojpWREQCwteWwTOBEuAQ0KN0/awqL58ZY24wxmw3xhQaY742xvh0vDGmlzEmxxiTW5frSf0Z1CWem8af5F1/aukWVm5ND2BEIiJNjl/rWBEREV/51DJore3mrwsaYy4FHgduwDPR7g3AR8aYvtbaXUc5LhR4DfgcGO2veOTE/X5UD77YcpgvtqRjLdz4+nd89JeziI8KDXRoIiKNnj/rWBERkbrwtWXQn24C5lhrn7fWrrfW/glPl5jrj3HcP4EfgDfrO0CpG4fD8OgvBhIfGQLAgexCbnn7B6y1xzhSREREREQCpc7JoDEmyRjTperLx2NDgSHAx1V2fQyMPMpx5wHnA3+ua7zSMNrGhvPwJQO864vXpTF3xY7ABSQi0gSdSB0rIiJSVz51EzXGOIB7gd8BrWopFuTDqVqXlkursj0NOKeWa7cHngcustbmGGOOFet1wHUAbdu2JTU11Yewapebm3vC52gpgoBxXYNZvNMJwKz31+E8tI2erTw/GrqX/qN76T+6l/6je3l8/FjHioiI1Imvo4lOA/6Ap6vmvcB9gBu4ovT9gTpet2r/QVPDtjKvAM9aa1f5dGJrZwOzAYYOHWpTUlLqGFplqampnOg5WpLhZ7i45LmV/Lg3C5eFF9bD+38aQWJ0mO6lH+le+o/upf/oXh63afi3jhUREfGJr91EfwX8HU9FBfCutfZuoA+wF/C1C8thwAW0q7I9ieqthWXOBu42xjiNMU7gBSCqdP06H68rDSQ8JIhnrhhMXITn+cH9WYVMe/07XG49PygiUgt/1bEiIiJ14msy2B1YY611AU4gAsBaWwI8Bvzal5NYa4uBr4FxVXaNA1bUclh/YGCF19+AgtJlDSbTCHVOiORfl5Y/P7hs82Ee/2RzACMSEWnU/FLHioiI1JWvyWAWEF66vA+oOKt4MJ7Jcn31KDDVGHOtMaaPMeZxoAPwHIAx5h/GmE/KCltrf6r4wvMtqbt0PbMO15UGdHbvtvxxTE/v+hOfbOaHQ84ARiQi0mj5s44VERHxma/PDH4L9AUWlb7uMcYU4PkG8z7gG18vaK193RiTCNwFtAd+AiZZa3eWFmmPZ9JdaeJuHHcS3+7O5Istnkno//1DEReOzadzQmSAIxMRaVT8VseKiIjUha/J4GN4urEA3A0MBv5bur4T+GNdLmqtfQZ4ppZ9U49x7BxgTl2uJ4ER5DA8cdkgzntiOQeyC8krgd+/8jVv/X4kEaEaGE9EpNRj+LGOFZHmJzs7m4MHD1JSUhLoUBpUXFwc69evD3QYjUJN9yIkJISkpCRiY2OP+7w+JYPW2sUVlg8YY4bhab2LBNaXPtcgUk1idBhPXzGYy2avpMRlWbsvm1ve/oEnLhvIsaYJERFpCVTHisjRZGdnk5aWRseOHYmIiGhRfz/l5OQQExMT6DAahar3wlpLQUEBe/fuBTjuhLDOk86XXtxaa7dYa39QJSXHMqRrPPdMPsW7vuD7ffz7820BjEhEpPFSHSsiFR08eJCOHTsSGRnZohJBOTpjDJGRkXTs2JGDBw8e93l8nXR+1LHKWGs/P+4opNn75eld+HjNelJ3ewaR+efCDfRuF0PKyUkBjkxEJLBUx4rI0ZSUlBARERHoMKSRioiIOKHuw74+M5hK7ZPCl9FDYHJUV/YJJS8ohtU7MrEW/vzqt7z3xzPp1joq0KGJiARSKqpjReQo1CIotTnRnw1fk8ExNWxLBM4HRqOH28UHwQ7DM1cMYfJTy9mfVUh2oZPrXl7Du384g+gwX38URUSaHdWxIiISED49M2it/ayG1zvW2l8D84EL6jdMaS7axITx76uGEBrs+dHbfDCXG1//Drf7WF+Ki4g0T6pjRaQlmTdvHo8++milbampqRhjSE1NDUxQdeByuXj44Yc5++yzadu2LTExMQwePJgXXngBt9sd6PDq7LgGkKniA+AXfjiPtBCndmrFAxf1964vXpfG459sDmBEIiKNlupYEWlWakoGBw8ezMqVKxk8eHCAovJdQUEB9957L6eccgqzZ89m3rx5jBkzht/+9rfceuutgQ6vzvzRN+9koOmlwRJQFw3uxNp92bywfDsAj3+ymT7tY5l4SrsARyYi0qiojhWRZi82Npbhw4cHOgyfREREsG3bNhISErzbxo4dS2ZmJk8++SR///vfm9SAPz61DBpjrq7hda0x5jHgAeCjeo1SmqXbz+3NmT1be9dvfuM7NqXlBDAiEZGGpzpWRFqKqVOnMnfuXPbu3YsxBmMMycnJNXYTTUlJYfz48SxcuJCBAwcSERHBoEGD+PLLL3E6ndxxxx20b9+ehIQEpk6dSl5eXqVr5efnc+utt9KtWzdCQ0Pp1q0b99133wl35QwKCqqUCJY57bTTKCoq4vDhwyd0/obma8vgnFq2FwGvA3/xSzTSogQHOXjy8kFMfno5uzMKyCt28es5q3nn+pEkxYYHOjwRkYYyp5btqmNFpFmZMWMGhw4dYvXq1cyfPx+AsLAwsrKyaiy/bds2pk+fzp133kl0dDS33HILkydPZvLkyTidTubMmcP69euZPn06SUlJPPjggwA4nU4mTJjAunXrmDFjBv3792fVqlXMmjWLjIwMHnnkEe81nE6nT7EHBx89bfrss89o1aoV7du39+l8jYWvyWC3GrYVWmvT/BmMtDzxUaHMvmooFz2zgoISF3syC7j6xa94/XcjiIsICXR4IiINwa91rDHmBmA60B5YC0yz1i7z4bhewDeAsdZGH8+1RUSOpkePHrRp04bQ0NBK3UJrGzgmIyODlStX0r17dwDcbjdTpkxh+/btLFmyBIAJEybw+eef8+abb3qTwVdffZXly5fz2WefMWqUZyrXsWPHAnDPPfdw6623kpTkmes6JMS3vzetrX2ww0WLFvHGG28wa9asYyaNjY1P0Vprd9Z3INJy9Wkfy9NXDOK3L3+Ny23ZcCCH3768hpd/PYzwEE2tJSLNmz/rWGPMpcDjwA3A8tL3j4wxfa21u45yXCjwGvA5nuksRKQRS77tg0CH4LXjgfPq7dw9e/b0JoIAvXv3BjwJYEW9e/dmwYIFWGsxxrBw4UK6du3KyJEjK7X8jR8/nrvuuotVq1YxefJkAFavXn1CMa5bt47LL7+clJSUFjuAjMgJO7t3W/558an89c3vAfhqewZ/ee1bnrliCEEOTbQqIuKjm4A51trnS9f/ZIyZCFwP3H6U4/4J/AB8hpJBEWkkWrVqVWk9NDQUgPj4+GrbnU4nLpeL4OBgDh48yM6dO2tt9UtPT/cuDxw48Ljj27ZtG+PGjaNbt27MmzevybUKgo/JoDHGDfg6EZy11ja9OyEB9/MhnTicW8QDH20AYNHaNO6a9xP3X3gKxighFJHmyV91bGnr3hDg4Sq7PgZGHuX65+GZ4H4wcLGPcYiINFqJiYl069aNN954o8b9ycnJ3uXj7Sa6Z88exo4dS2xsLAsXLiQ2Nva44w0kX5O2WcCvgAhgAZAGtMNTeeQDL+F7RSZSq9+N6s6hnCLvlBOvfrWL1tGh3Dz+5ABHJiJSb/xVx7YGgkqPrygNOKemA4wx7YHngYustTm+fPFmjLkOuA6gbdu2JzRJdG5ubpOYZLq+6T6U073wqHgf4uLiyMmpPNr6j3eOCkBUNasa27E4HA7y8/MrHZefn+99L9vucrmw1lYql5ubC0BhYWGl7UVFRd5YgoODSUlJ4e2338YYw0knnXTUuH39eat4vcOHDzNhwgSstbz77ruEh4fX+T7UlcvlqvUahYWFx/3vxtdksATYCUyw1uaXbTTGRAGLgBJr7X3HFYFIBcYY7pzUh/TcIuZ9tw+AJz/dQkRoEDek9AxwdCIi9cLfdWzVxNHUsK3MK8Cz1tpVPp/c2tnAbIChQ4falJSUOoRWWWpqKidyfHOh+1BO98Kj4n1Yv349MTExgQ3IjwYMGMBLL73EK6+8wtChQwkPDycyMhKAyMhI72cNCgrC6XRW+uzR0Z6xrcLDwyttDwsLAyAmJobg4GB+85vf8OqrrzJ58mRuvvlmBgwYQHFxMVu3bmX+/PnMmzfPe83Ro+vWM76goICLL76YXbt28eKLL3LkyBGOHDni3d+3b996aSXMycmp9ecgPDycQYMGHdd5fU0Gfwf8oWIlBWCtzTPGPAw8CSgZFL9wOAwP/nwAmfklfLbpEAAPLtxIeHAQvz6zpkH3RESaNH/VsYcBF55WxYqSqN5aWOZsYLQx5u7SdQM4jDFO4IbSxE9ExG+uvfZaVq1axR133MGRI0fo2rUrc+bM8es1QkJCWLRoEQ888ACzZ89m+/btREVF0aNHD8477zzvs4fHIy0tjW+//RaAK664otr+pUuXNqkvNHxNBlsDtd21UCDRP+GIeIQGO3juyiH8as5XrNqWAcDf319HeEgQvzy9S4CjExHxK7/UsdbaYmPM18A44M0Ku8YBb9dyWP8q61OAO4FhwF5frisiUhdRUVG8+uqr1bZXfSYvNTW1WrfI5OTkGqd4mDlzJjNnzqy0LTw8vMbtJ6q2GJoqh4/l1gD3GGM6VtxYuj4TOLExWUVqEBEaxAvXnMaQruUjRt0570fe+WZPAKMSEfE7f9axjwJTjTHXGmP6GGMeBzoAz5We8x/GmE/KCltrf6r4wpMAukvXM0/sY4mISGPna8vgn4FPga3GmFV4upu0BYbjebj9l/UTnrR0UWHBvPSr07ji+S/5cW8W1sJf3/ye0GAH55/aIdDhiYj4g9/qWGvt68aYROAuPJPO/wRMqjCXYXughx9jPy7ZCxeS/tJLxGdns/3ZZwk/6SSSbr6ZoCrDyIuISP3yddL5b40xPfHMXzQcT7eS/XiGr/6XtTb9aMeLnIjY8BBe/vUwLn9+FRsO5OC28JfXvsPltkwZ2PHYJxARacT8Xcdaa58Bnqll39RjHDsHmFOX6x0P5+F0Cr//gVCgECj8/gdCOneh9XW/re9Ly1FYa6GkBOt0Yp1OHOHhmBN4tiqQrLXgdoPLha3wbp1OcLuxLpfn3ekCd+my213p3bNsK+z3LHv2ly1bsG5COnUirFv9jWvg7ZZYsXtiTcvWlo8WVcv+6tu8/yndVKULZNVtNZY5xjrVu4Ee8xhrMUVFOIuLjxmHrWGb77EevUzl3TV0D63p+KNtP45jLGAiIqAeBhLyeT7A0sroTr9HIOKD+KhQXrn2dC7990q2HsrD5bbc+LonIbxocKdAhycickJUx0LJ/n2BDqHBWZcLd0EBjqwsinfswF1QgLugEFtcjC0uwl1UhC0q9q7boiLcxcWebUVFnjJV1m1xeTLneZVAiRNbUsv2CttwOivFZ0JCiL/ySqLOPAPKypVUOL62bSVVrlO2rVIMpdd3ubAuJ7jcWLeL+IxMtj/zbJUkzgXOykkdTmfldZfLk+BVWG9owe3aeVq33W6wZcmj53NRtmzd4CpNNCsmrGXLpduTSkpYbwy43ZQ88TiFAfg8jUEQnuGWBfBxPsS68nXSeQfgsNY6K2ybAJwCfGqt/bZeohOpoHV0GK9dN4JfPr+KzQdzcVu4+c3vcbosvzitc6DDExE5Li2xjo2dMJ6IU/qx/okniVixAgBb0nT+5LPW4s7NxZWVhTs7G1fpy7OcgysnG3dWNq6cHNy5ubjz8z2JXn4eNr/Au25L50ZrA2wN7EeqkS0pIeOll8h46aUGu2ZZa3FT5DxwAOeBA34517Fn/BTxD19bBl8FioCrAYwxv6e8C0qJMeY8a+2SeohPpJI2MWG8dt1wrvjPl2w4kIO1cMvbP+B0W40yKiJNVYurY4PbtCG4TRuKT+rlTQY5gWSwZO9enJlHCO/XF2NO7M9oV3Y2JXv3UrJvHyX79uM8fBhn+mFch9Nxpqd7l21Z17XmJigIExzsTVSbvOBgjMPh+VxBQZ73snWHA4KDMI4gcDgqbDdgSpeNAYcDghwY4ygv5yhfLtm7l+KdO48diz9U/Pk2pjRpNOXZY5X91LLf1Lbfl/Wq/8a8cdR2fA3HHPOanm0lJSWEhoYctcwxz3OUaxsfytR+nmPsO0r5Wn9PHeWY/HpqHfY1GRwO3FphfTrwH+BmPBPP3gk0q4pKGq/E6DBe/e1wrnzhS9buywbgjnd/pNjpYuoZmodQRJqcFlvH2qAg77K7DsmVdTrJX72a7EWLyFu2nJK9nlkwEn/7W5JuvumYx7sLCijato3iLVso2rKFom3bKdmzh5J9+3BXGcq+vpnISJzBwYTHxeGIjMSEh+EIDcOEeV6OsFBMSKh33YSG4AgLw4RWWQ8L85QLDcWEBGOCPS+CgzHBIeXbQkJKt4dUKuct6/AMNF+8ezcHH3wQV+YRCCk9R1nZkArnrW1bSJXrl1272vag8kTN4eC7H35g8NChULpuak3igo++XvbeQEr27cOZkelJJB1VkkZjvJ/Pm1w6gmop68A4DJ8vX87olBRwONiwaRPhffqc8BcdTVFRTg4h9fCcXJNUT7+bfE0Gkyidb6j0IfduwFPW2hxjzEvA/+olOpFaxEeF8r9rh3P1i1/y/Z4sAGYuWMeRghL+MrZXi/yFKSJNVsutY4PL/wzxpZtoyYEDZL7yCkfefgdXZvWZL7IXLqyWDFqnk6JNm8j/7jsKv/+egu++p3jXrhoHuKgLExFBUHwrgmJiCYqNxREbS1BMDEFxsTi822IIio7GERnpSfQiInFEReKIiMAREYEJD8c4HKSmpja6SapDO3em05NPNvh1S/LziRg4sMGve6JCOnQgpIMfRzkPDS0fvMcY/V0j9cbXZDCb8klvU4DD1tofStddQLif4xI5prjIEP7v2tOZ+uJXfLPrCACPLdlMZl4xd1/QD4dDvzhFpElosXWsDfItGXRlZ3Po8SfIfP31aoOcVFSyZw/uwkJc2dnkLVtGbupn5K1YgTsvz+eYTHi45w/7jh0J6dDB06W1dSJBiYkEJ7YmuHUiwYmJOKKifD6niEhj5WsyuAK4zRjjBKYBH1bY1xPQLOASELHhIbxy7en8/pVv+HzTIQDmrtxJZn4Jj/xiACFBDddFRETkOLXcOrZiy2At3URzPvmE/X+7G1d65Rk2gpOSiJkwgZhx57Dvtttw7tsP1rLzl1dQuG7d0a/rcBDapQthvXoS2rMnYT16EtqlMyEdOxKUkKBWGBFpMXxNBm8BPgDmA9uAmRX2XQqs9G9YIr6LDA3mP1cP5aY3vuP9H/YDMP/7fWQXlvDsFUOICA06xhlERAKqxdaxNrj893PVlkFrLYeffobDTz1VaXvEkCEk/ubXRI8e7XkOCwjvdRK5+zy//2tKBIPbtiVi0CAiBg4gYsAAwvv2xREW5u+PIyI+mDdvHtu2beOmm8q7dKempjJmzBiWLl3a6LpM12T58uX85z//Yc2aNWzYsIFOnTqxY8eOGsuuXbuWG2+8kRUrVhAWFsbkyZN55JFHSEhIaNiga+HrpPObgZOMMYk1TH77F8A/4+iKHKfQYAePXzaIuIgQ/vvlLgBSNx7i8udX8cI1Q0mMVqUvIo1Ti65jj/LM4OGnnubw00+XF23Xjra33UbMhPHVWu7CTjqJ3M8+K9/gcBA5dCjRo0cTPXoUoT16qLVPpJGYN28eS5YsqZQMDh48mJUrV9K3b98ARua7Tz75hGXLljF06FCMMeTUMrjLvn37SElJoXfv3rz11lscOXKE6dOnc/7557N8+XIcDTjIUW18nnQevJPiVt32o//CETl+QQ7DvT87hcSoUJ74dAsA3+0+wkXPruClqafRvU10gCMUEaldS6xjbS3J4JF58yolglEjR9Dx0Uc9E3rXIP6yS8n76ksocRI76VxiL5hMSNukeotbRPwrNjaW4cOHBzoMn82YMYO7774bgCuvvJLly5fXWO6hhx6ipKSEBQsW0Kr091eHDh0YPXo08+bN46KLLmqokGsV+HRUxI+MMdw0/mRmXtDXO1XLzvR8Lnp2BWt2ZAQ2OBERqaTi1BJlzwwW79lD2qx7vdujzjqLTs89V2siCBDSsSPdXn+dbu+8TeK11yoRFGmkpk6dyty5c9m7dy+mdJTU5ORkUlNTMcaQmprqLZuSksL48eNZuHAhAwcOJCIigkGDBvHll1/idDq54447aN++PQkJCUydOpW8KgNF5efnc+utt9KtWzdCQ0Pp1q0b9913H263+4Q/h68tevPnz+e8887zJoIAo0aNokuXLrz33nsnHIc/1KllUKSpmHpGN9q3iuAvr31LYYmbI/kl/PI/X/LoLwZw/ql+HPpZRESOXw0tgwf+/nfv6J+hXbvS6fHHcJQNsS8iTdqMGTM4dOgQq1evZv78+QCEhYWRlZVVY/lt27Yxffp07rzzTqKjo7nllluYPHkykydPxul0MmfOHNavX8/06dNJSkriwQcfBMDpdDJhwgTWrVvHjBkz6N+/P6tWrWLWrFlkZGTwyCOPeK/hPMoIxRUFB9ctbSooKGD79u1ce+211fb169ePdcca6KqBKBmUZmtCv3a8dt0Irp27msO5xRQ73fzxf9+yMz2fG1L0/IiISKBV7SZa8OOP5H2+zLPB4aDDQw/iiIwMUHQi4m89evSgTZs2hIaGVuoWWrFFsKKMjAxWrlxJ9+7dAXC73UyZMoXt27ezZMkSACZMmMDnn3/Om2++6U0GX331VZYvX85nn33GqFGjABg7diwA99xzD7feeitJSZ4eBCEhIT7Fbus4N2lmZibWWuLj46vtS0hIYOPGjXU6X31RMijN2sDOrXjn+jOYOucrth3yfNP80KKNbErL4Z8Xn0p4iEYaFREJmEpTS5Rw+NnnvOuxkyYRceqpgYhKpPGbGRfoCMrNrLlVzx969uzpTQQBevfuDXgSwIp69+7NggULsNZijGHhwoV07dqVkSNHVmr5Gz9+PHfddRerVq1i8uTJAKxevbpeYi9LHmtqfKhrYlmflAxKs9clMZJ3rh/J7/7va77c7nlu8L3v9rH9cB6zrxpKu7hmO5+ziEijVvGZQdfhw+R++ql3vfXvrgtESCLSiLSq8qxwaGmX8aqtbaGhoTidTlwuF8HBwRw8eJCdO3fW2uqXXmHe0oEDB/o15jLx8fEYY8jIqD5mRWZmZtOaWkKkqWsVGcr//eZ07p6/lle/8kw98cOeLCY/tZx/XzWEQV2qN+GLiEg9q+UZnKhRZxHWq1cDByMizUViYiLdunXjjTfeqHF/cnKyd7m+uolGRkaSnJzM2rVrq+1bt24do0ePrtP56ouSQWkxQoMd3H/hKfRtH8PMBetwuS0Hc4q4dPYq7r+wPz8f0inQIYqItCi2lmQwbvKUBo5EpImpx66Z9S0sLIyCgoJ6vcbEiRN5++23iY6O9nYtrU19dRMFmDx5MnPnziUrK4u4OE/X3uXLl7Nz505vN9VAUzIoLYoxhqtGJNOjTTQ3/O8bjuSXUOx089c3v+frnZncfUFfPUcoItJQgmr4fRscTHRK4/jGXET8r2/fvmRkZPDss88ydOhQwsP9/7jOFVdcwUsvvcTYsWO5+eabGTBgAMXFxWzdupX58+czb948IksHpxo6dGidz3/o0CE+++wzAHbt2kV+fj5vvfUW4Pl8ffv2BWD69Om88sorTJ48mdtvv52srCxuueUWhg0bxoUXXuinT3tilAxKizSyZ2ve+8MZ/PblNWxKywXg1a92sXZfFk//cjCdEzR6nYhIvTMGQkKgwoTzkQMHEhQdHcCgRKQ+XXvttaxatYo77riDI0eO0LVrV+bMmePXa4SEhLBo0SIeeOABZs+ezfbt24mKiqJHjx6cd9553mcPj9fatWu55JJLKm0rW7/77ruZOXMmAB07dmTp0qXcdNNNXHzxxYSGhjJlyhQeeeQRn+cqrG9KBqXF6poYxbs3nMFt7/zIgu/3AZ7nCC94ajmPXTqQlJM1abGISH1zhITgrpgMVhhuXkSan6ioKF599dVq26s+k5eamkpOTk6lbcnJyTU+uzdz5kxvAlYmPDy8xu3+kJKS4vMzhP3792fx4sV+j8FfGkdKKhIgUWHBPHHZQGZe0Jdgh2fo3yP5Jfxqzmoe+XgjTpc7wBGKiDRvpsrgDRGn9g9QJCIiLY+SQWnxjDFMPaMbr/9uBO1iPf3WrYUnP93CL5//kn1H6vchZxERKRfeX8mgiEhDUTIoUmpI13je//OZnNEz0bvtqx0ZTHpiGYvXpQUwMhGR5suVVXlUxOB4TfUjItJQlAyKVNA6OoyXf306fx1/EqW9RjmSX8JvX17DPQvWUuR0BTZAEZFmLKRjx0CHICLSoigZFKkiyGH449m9eP13I2gfVz7c8Utf7ODiZ1ew/XBeAKMTEWm+Qjp0CHQIIiItipJBkVqclpzAh38+i3P6tPVu+2lvNuc/sYz3vtsbwMhERJonJYMiIg1LyaDIUcRHhfL81UO4+4K+hAZ5/rnkFbv4y2vfcctb35Nf7AxwhCIizUd43z6BDkFEpEUJSDJojLnBGLPdGFNojPnaGHPWUcqmGGPeM8bsN8bkG2N+MMb8uiHjlZbNGMOvzujGOzeMJDmxfDL6N9bsYfJTX7DhQHYAoxMRadra3nE7ACGdOtHq0ksDHI2ISMvS4MmgMeZS4HHgfmAQsAL4yBjTpZZDRgI/Aj8HTgGeBWYbY37ZAOGKeJ3SMY73/3wWPxtY3o1py8Fcpjz1Bf/9cqfPk4+KiEi5hKuvpsfHi+j+4Qc4wsOPfYCIiPhNIFoGbwLmWGuft9aut9b+CdgPXF9TYWvt/dbau6y1X1hrt1lrnwXeAS5uwJhFAIgOC+Zflw7koZ+fSkRIEABFTjd3vvsTf/zft2QVlAQ4QhGRpie0SxccoaGBDkNEpMVp0GTQGBMKDAE+rrLrYzwtgL6KBTL9FZdIXRhjuGRoZxb86Qx6t4vxbv/gx/1MfOxzFv50QK2EIiIiIjWYN28ejz76aKVtqampGGNITU0NTFB14HK5ePjhhzn77LNp27YtMTExDB48mBdeeAG3212t/Nq1axk/fjzR0dEkJibyq1/9ioyMjABEXrOGbhlsDQQBVWfwTgPa+XICY8z5wFhgtn9DE6mbnkkxzPvDGVw5vLyH8/6sQn7/ytf8Zu4admfkBzA6ERERkcanpmRw8ODBrFy5ksGDBwcoKt8VFBRw7733csoppzB79mzmzZvHmDFj+O1vf8utt95aqey+fftISUmhoKCAt956i6effpolS5Zw/vnn15g4BkJwgK5btdnE1LCtGmPMGcD/gD9ba7+qpcx1wHUAbdu2PeFvGHJzc5vEtxRNQXO9l+e0griBYfzfuiKyiz3bPt1wkOWbDjK5RwgTu4UQXDaDvZ8013sZCLqX/qN7KSIixyM2Npbhw4cHOgyfREREsG3bNhISErzbxo4dS2ZmJk8++SR///vfiYiIAOChhx6ipKSEBQsW0KpVKwA6dOjA6NGjmTdvHhdddFEgPkIlDZ0MHgZcVG8FTKJ6a2ElxpgzgQ+Bv5U+N1gja+1sSlsNhw4dalNSUk4kXlJTUznRc4hHc76XKcB1+SU8uGgD//tqF9ZCsRve2lzCd1lh3PuzUxjePdFv12vO97Kh6V76j+6liIgczdSpU5k7dy7geewGoGvXrsyZM4cxY8awdOlSbz2SkpJCUVERd999N7fddhsbN26kd+/ePPfccwwZMoS//e1vvPTSSxQVFTF58mSefvppoqKivNfKz8/nnnvu4Y033mDv3r107NiRa6+9lttvvx2H4/g7RwYFBVVKBMucdtppvPTSSxw+fJjOnTsDMH/+fM477zxvIggwatQounTpwnvvvdfykkFrbbEx5mtgHPBmhV3jgLdrO84YMwr4AJhprX2sXoMUOU5xkSHcd2F/fj6kE3e++xPr9numnNhyMJfLZq9iysAO3HZub9rHRQQ4UhEREZGGN2PGDA4dOsTq1auZP38+AGFhYWRlZdVYftu2bUyfPp0777yT6OhobrnlFiZPnszkyZNxOp3MmTOH9evXM336dJKSknjwwQcBcDqdTJgwgXXr1jFjxgz69+/PqlWrmDVrFhkZGTzyyCPeazidvs0ZHRx89LTps88+o1WrVrRv3x7wdCfdvn071157bbWy/fr1Y926dT5dt74Fopvoo8D/GWO+Ar4Afg90AJ4DMMb8AxhmrR1bup6CJxF8BvivMaasVdFlrT3UsKGLHNugLvHM/+MZzF25k0c/3khesQuA977bx8dr07g+pQfXjepOeOlopCIiIiItQY8ePWjTpg2hoaGVuoXW9ohBRkYGK1eupHv37gC43W6mTJnC9u3bWbJkCQATJkzg888/58033/Qmg6+++irLly/ns88+Y9SoUYCnKyfAPffcw6233kpSUhIAISEhPsV+tMEBFy1axBtvvMGsWbO8SWNmZibWWuLj46uVT0hIYOPGjT5dt741eDJorX3dGJMI3AW0B34CJllrd5YWaQ/0qHDIVCAS+Gvpq8xOILm+4xU5HsFBDn5zZjfO69+eWe+v44Mf9wNQUOLi0cWbeH31bu6Y1IdJ/dt5u0mIiIiI+KL/3P6BDsHrx2t+rLdz9+zZ05sIAvTu3RvwJIAV9e7dmwULFmCtxRjDwoUL6dq1KyNHjqzU8jd+/HjuuusuVq1axeTJkwFYvXr1CcW4bt06Lr/8clJSUioNIFOWPNb0d15jGnU+IAPIWGufwdPSV9O+qTWsT62prEhj1y4unKevGMxV29K5Z8E61pd2Hd17pIA//O8bhiUn8LcL+nJKx7gARyoiIiLSuFR81g4gtHQ+0qqtbaGhoTidTlwuF8HBwRw8eJCdO3fW2uqXnp7uXR44cOBxx7dt2zbGjRtHt27dmDdvXqWupPHx8RhjapxGIjMzs8bnDgMhUKOJirQow7sn8v6fzuT11bt5+OONZOR5hh39akcGFzy1nEuHdubm8SfTJiYswJGKiIiING2JiYl069aNN954o8b9ycnJ3uXj7Sa6Z88exo4dS2xsLAsXLiQ2NrbS/sjISJKTk1m7dm21c61bt47Ro0f7dN36pmRQpIEEOQy/PL0L553anic+2czcFTtwui3WwmurdzP/+3386oxkrjurB3GRvv1iEhERkZanPrtm1rewsDAKCgrq9RoTJ07k7bffJjo62tu1tDbH00300KFDnHPOOQAsXryYNm3a1Fhu8uTJzJ07l6ysLOLiPL3Ali9fzs6dO73dVANNyaBIA4uLCGHG+X25fFgX7vtgHUs3esZByi928fTSrby8cifXndWdX53Zjegw/RMVERGR5qNv375kZGTw7LPPMnToUMLDw/1+jSuuuIKXXnqJsWPHcvPNNzNgwACKi4vZunUr8+fPZ968eURGRgIwdOjQOp27oKCACRMmsGPHDl588UX27NnDnj17vPv79u3rbSWcPn06r7zyCpMnT+b2228nKyuLW265hWHDhnHhhRf67wOfAP2lKRIgPZOieelXw1i64SD/XLiBDQdyAMgpdPLI4k28tGIH14/uwVUjumrkUREREWkWrr32WlatWsUdd9zBkSNHvPMM+lNISAiLFi3igQceYPbs2Wzfvp2oqCh69OjBeeed53328HikpaXx7bffAp6ks6qKcyV27NiRpUuXctNNN3HxxRcTGhrKlClTeOSRR05orkN/UjIoEmBjeicx+qQ2vP/jfh5bvIlth/MAyMgr5r4P1/P8sm38bnQPLh/WmchQ/ZMVERGRpisqKopXX3212vaqz+SlpqaSk5NTaVtycnKNI3HOnDmTmTNnVtoWHh5e4/YTVVsMtenfvz+LFy/2awz+1DhSUpEWzuEwTB7QgY9vHMVDPz+VTvHlE9MfzCli1vvrOPOfS3nq081kFZQEMFIRERERaS6UDIo0IsFBDi4Z2plPb05h1s9OoW1s+eiiGXnFPPzxJs584FPe2lRMem5RACMVERERkaZOyaBIIxQa7OCq4V35bPoYZv3sFDq2Km8pzCly8v62Es7456fcs2AtuzPyAxipiIiIiDRVSgZFGrHwkCCuGt6V1OkpPHzJALq3ifLuKyxx89IXOxj90FL++L9v+H73kcAFKiIiIiJNjpJBkSYgJMjBz4d0YvGNo3n6l4PpHFP+T9dt4f0f9jPl6S/4xb9XsnhdGm637w82i4iIiEjLpKEJRZqQIIfhvFPbE5m+ATr04z/LtvHFlnTv/q+2Z/DV9gy6t47iN2d14+LBnTQthYiIiIjUSMmgSBNkjCHl5CTGnJzET3uzeGH5dhZ8vw9naYvgtsN53PnuTzzy8SZ+MbQzvxzWhS6JkQGOWkREREQaE3UTFWniTukYx78uHciyW8fwu1HdiQkr/44nI6+Y5z7byqiHlnLVC1+y8Kf9lLjcAYxWRERERBoLtQyKNBPt4yK4fVIf/jS2F6+v3s2Ly7ez90iBd/+yzYdZtvkwSTFhXHpaZy49rTOd4tVaKCIiItJSKRkUaWaiw4L5zZnduGZEV1I3HuK/X+4kddMhbOmYMgdzinjy0y08tXQLI3skcvHgTkw8pR2Rofp1ICIiItKSqJuoSDMVHOTgnL5teelXw1h2yxj+OKYnbWLKJ7G3Fr7Yks5Nb3zPafcuYfqb37NqW7pGIhUREZF6M2/ePB599NFK21JTUzHGkJqaGpig6mjmzJkYY6q9fvazn1Uru3btWsaPH090dDSJiYn86le/IiMjo+GDroWaAkRagE7xkfx1wsn85ZxeLFmXxv++2sXyLYe9rYV5xS7e/HoPb369h84JEVw0qBMXD+6kQWdERETEr+bNm8eSJUu46aabvNsGDx7MypUr6du3bwAjq7vly5cTFFQ+antCQkKl/fv27SMlJYXevXvz1ltvceTIEaZPn87555/P8uXLcTgC3y6nZFCkBQkJcnBu//ac2789+7MKePfbvbz99R62HsrzltmdUcDjn2zm8U82M7BzK84/tT3nndqe9nERAYxcREREmqvY2FiGDx8e6DDq7PTTTyc4uPZ06qGHHqKkpIQFCxbQqlUrADp06MDo0aOZN28eF110UQNFWrvAp6MiEhDt4yK4IaUnS24azbs3jOTK4V2IiwipVOa73Ue494P1jPjHp/z82RXM+WI7B3MKAxSxiIiINGVTp05l7ty57N2719u1Mjk5ucZuoikpKYwfP56FCxcycOBAIiIiGDRoEF9++SVOp5M77riD9u3bk5CQwNSpU8nLy6t0rfz8fG699Va6detGaGgo3bp147777sPtbrhR1efPn895553nTQQBRo0aRZcuXXjvvfcaLI6jUcugSAtnjGFQl3gGdYlnxvl9+WT9Qd76eg+fbzrknbcQYM3OTNbszOSe99dxercEzj+1AxP6tav0HKKIiIhIbWbMmMGhQ4dYvXo18+fPByAsLIysrKway2/bto3p06dz5513Eh0dzS233MLkyZOZPHkyTqeTOXPmsH79eqZPn05SUhIPPvggAE6nkwkTJrBu3TpmzJhB//79WbVqFbNmzSIjI4NHHnnEew2n0+lT7DW1AHbu3JmDBw/SqVMnLrvsMmbOnElEhKcnVUFBAdu3b+faa6+tdly/fv1Yt26dT9etb0oGRcQrLDiISf3bM6l/e47kF7No7QHe/2E/K7am4ypNDK2FVdsyWLUtgxnv/cSgzq0Y17cd4/u1pUeb6AB/AhEREWmsevToQZs2bQgNDa3ULbS2gWMyMjJYuXIl3bt3B8DtdjNlyhS2b9/OkiVLAJgwYQKff/45b775pjcZfPXVV1m+fDmfffYZo0aNAmDs2LEA3HPPPdx6660kJSUBEBISUvWyNbK2/Avynj178sADDzBo0CCMMXz88cf861//4ptvvmHx4sUAZGZmYq0lPj6+2rkSEhLYuHGjT9etb0oGRaRGrSJDufS0Llx6WhfSc4v46KcDvP/DPr7cnuEdeMZa+GbXEb7ZdYR/LtxA9zZRjOvblvF92zKoczwOhwnshxAREWmG1vfuE+gQvPpsWF9v5+7Zs6c3EQTo3bs34EkAK+rduzcLFizAWosxhoULF9K1a1dGjhxZqeVv/Pjx3HXXXaxatYrJkycDsHr16jrHdeWVV1ZaHzduHJ06dWLatGksWbKEc845x5s8GlP9b6GKiWWgKRkUkWNKjA7jyuFduXJ4Vw5mF/LRTwf48Mf9rN6RQcWZKLYdyuPfn23j359to3V0GGf3bkPKyUmc0bN1tecRRaR+GGNuAKYD7YG1wDRr7bJayqYANwLDgDhgC/CYtfbFBglWROQoKj5rBxAaGgpQrbUtNDQUp9OJy+UiODiYgwcPsnPnzlpb/dLT073LAwcO9Eusl19+OdOmTWP16tWcc845xMfHY4ypcRqJzMzMaiOPBoqSQRGpk6TYcK4Zmcw1I5PJzCvm0w0H+XjdAT7fdJiCEpe33OHcIt5Ys4c31uwhyGEY1LkVKSe3YfRJSfTrEKtWQ5F6YIy5FHgcuAFYXvr+kTGmr7V2Vw2HjAR+BB4E9gMTgNnGmEJr7f8aKGwREb9KTEykW7duvPHGGzXuT05O9i4fTzfRoylrCYyMjCQ5OZm1a9dWK7Nu3TpGjx7t0/nqm5JBETlu8VGhXDykExcP6URhiYsvthxm8bo0lqxP43Busbecy229A9A8/PEmWkeHMqpXG0ad1IaRPRJJig0P4KcQaVZuAuZYa58vXf+TMWYicD1we9XC1tr7q2x61hgzBrgYUDIo0kjVZ9fM+hYWFkZBQUG9XmPixIm8/fbbREdHe7uW1uZ4uonW5L///S/gmW6izOTJk5k7dy5ZWVnExcUBnrkJd+7c6e2mGmhKBkXEL8JDghjbpy1j+7TF5bZ8t/sIn206xGcbD/LD3iwqfqF2OLeYd77dyzvf7gWgV1I0I3skMqJHa4Z3T6BVZGiAPoVI02WMCQWGAA9X2fUxnhZAX8UCe/wVl4hIRX379iUjI4Nnn32WoUOHEh7u/y+Er7jiCl566SXGjh3LzTffzIABAyguLmbr1q3Mnz+fefPmERkZCcDQoUPrfP5BgwZx9dVXc/LJJ2OMYfHixTz55JNMnDiRMWPGeMtNnz6dV155hcmTJ3P77beTlZXFLbfcwrBhw7jwwgv99nlPhJJBEfG7IIdhSNd4hnSN56ZxJ5GeW8SyzYf5bNMhPt90iPS84krlNx/MZfPBXOau3Ikx0K9DLCN7tGZEj0SGJScQFaZfVSI+aA0EAWlVtqcB5/hyAmPM+cBY4IyjlLkOuA6gbdu2tY4C6Ivc3NwTOr650H0op3vhUfE+xMXFkZOTE9iA/OjSSy9l2bJl3HHHHRw5coQuXbrw7LPPAp65Acs+q8vlwlpb6bPn5uYCUFhYWGl7UVERADk5Od4pIN566y0effRRnnvuOXbu3ElkZCTdunVjwoQJFBUV4XKVP9pSV927d+eJJ54gLS0Nl8tFt27duPXWW5k2bVqluGJjY3n//fe54447uPjiiwkNDWXSpEncd9991eZFPBaXy1Xrz0FhYeFx/7sxjWk0G38bOnSoXbNmzQmdIzU1lZSUFP8E1MLpXvpPU76Xbrflp31ZpG48xBdbDvPtriMUu2qfADbYYRjYuRUjeyQyrFsiA7u0ItqPyWFTvpeNzYncS2PM19baun89K17GmA7AXmBUxQFjjDF3A5dba4/aV8oYcwbwEXCrtfZZX655ovWs/v156D6U073wqHgf1q9fT58+jWf00IaUk5NDTExMoMNoFI52L471M3K0OlZft4tIg3I4DKd2asWpnVrx57G9KCh28fXOTFZsPcyKren8sOdIpRFKnRWeN4QtOAyc3C6WwV1aeVsfuyRE1jh0s0gLcxhwAe2qbE+iemthJcaYM4EPgb/5mgiKiEjTp2RQRAIqIjSIM3u15sxerQHILixh9fYMvtiSzoqth9lwoHKXCLeF9fuzWb8/m/9+6RkcMTEqlMGlieHgLvH07xhHRGhQg38WkUCy1hYbY74GxgFvVtg1Dni7tuOMMaOAD4CZ1trH6jVIERFpVJQMikijEhse4h2IBiA9t4hV2zJYue0wX+88wsYD2ZVaDgHS84pZvC6Nxes8jR9BDkOvpGj6d4zj1E5x9O/Uit7tYggPUYIozd6jwP8ZY74CvgB+D3QAngMwxvwDGGatHVu6noInEXwG+K8xpqxV0WWtPdSwoYuISENTMigijVpidBjnndqe805tD0BukZPvdx/h652ZfLMrk292ZpJd6Kx0jMtt2XAghw0Hcnjza8+giMEOw8ntYjzJYcdW9O8Yx0ntohv884jUJ2vt68aYROAuPJPO/wRMstbuLC3SHuhR4ZCpQCTw19JXmZ1Acn3HKyIigaVkUESalOiwYM7o2Zozenq6lbrdlq2HcvlmV2ZpgniErYdyqTo2ltNtWbsvm7X7snmV3YAnQWwXCUMOfEvvdrH0bh9D3/axJMWE6RlEabKstc/gaemrad/UGtan1lRWRESaPyWDItKkORyGXm1j6NU2hktP6wJ4Wg/X7s3ix71Z/LDH8779cPUhnJ1uy55c2PPdPt5jn3d7fGSINzns0z6WPu1i6dU2Wt1MRUREpFlRMigizU50WDCnd0/k9O6J3m1ZBSWs3ZfFj3uy+GGv531XRn6Nx2fml7ByWzort6V7tzkMdEmIpGdSND2SounZJppebWPo0SaKmPCQev9MIiIiIv6mZFBEWoS4iBBG9mjNyB6tvdtyi5y89uFnRHToyfr92WzY73nOMLfIWe14t4Ud6fnsSM9nyfqDlfa1iw2nZ1J0tVdiVKi6m4qIiEijpWRQRFqs6LBgesYHkXJ6V+82ay17Mgs8yeGBHDYcyGb9/hx2pOdVew6xzIHsQg5kF7J8y+FK2+MjQ0huHUVyYhRdEiJJbh1J18QouiZEkqBEUURERAJMyaCISAXGGDonRNI5IZLx/crn7i4scbHtUB6bD+aw9WAuWw7lsuVgLtsP51HiqjlLzMwvIXPXEb7ddaTavpiwYLokRnoSxcRIkhMj6ZIQRaf4CNrFhRMS5KivjygiIiICKBkUEfFJeEgQfTvE0rdDbKXtTpebnRn5bDnoSQ4rJor5xa5az5dT5PSOblqVw3i6nnaMj6BTfCQdW0XQMT6Cjq0i6BQfQYdWERrMRkREmqR58+axbds2brrpJu+21NRUxowZw9KlS0lJSQlccD5avnw5//nPf1izZg0bNmygU6dO7Nixo8aya9eu5cYbb2TFihWEhYUxefJkHnnkERISEiqV2717NzfeeCOLFy/GWss555zDY489RpcuXer1sygZFBE5AcFBDnq0iaZHm2gm9Cvfbq1lf1YhO9Pz2Zmex470fHZl5LHjcD67MvJrfC6xjNvCvqxC9mUVsnpHZo1lWkeHeZLFVhG0jwunXVw4bWM97+1iw0mKDSMsWAmjiIg0LvPmzWPJkiWVksHBgwezcuVK+vbtG8DIfPfJJ5+wbNkyhg4dijGGnJycGsvt27ePlJQUevfuzVtvvcWRI0eYPn06559/PsuXL8fh8PQCys/P5+yzzyYsLIy5c+dijOGuu+5izJgx/PDDD0RFRdXbZ1EyKCJSD4wxdGjlacUb0SOx0j5rLel5xexMz2Nn6aA0u9Lz2JWRz94jBaRlFx3z/IdzizicW8T3u4/UWiYhKpS2seG0L0sUY8NpFxdWKWmMiwjRs4vS4NK2Z7Pzp8Mc3OHmq5xtJHSIpvvA1jjUPVqkRYqNjWX48OGBDsNnM2bM4O677wbgyiuvZPny5TWWe+ihhygpKWHBggW0atUKgA4dOjB69GjmzZvHRRddBMDzzz/Ptm3b2LhxIz179gTg1FNPpVevXvz73/+ulDj7m5JBEZEGZoyhdXQYraPDGNI1odr+IqeL/UcK2XukgD2Z+ezNLGDPkQL2ZBawN7OAA9mFuNy1jGZTQUZeMRl5xazfX70rapnQIAeto0NpHRNWGlMobbzLYd7lNtFhxEYEK3EUv0jbkc3qD3YAcGit533sNX3oPaJ94IKSRslaCxbc1oLbs24tWLctX7YW6676XnFf+TIW3G7Puy/ly69V4Zru8tjK1sOjQ4iODyvdQaXr2dLRx6wbLLbK/rJr4V3GQu5+y8616Vi3xel2UZRfUuW+1HbDPP+pttvWWAxjwBFkPAulgoKP8nu+luse63plfnPtr/m/V14G8NYnXbt05fnnX2D8hHP4eNESRo8aDcC48WMpKi7mzjvu4q677mDT5k2cfNLJPP30MwwZMoSZM+9m7stzKSoq4vzzL+DJx5+s1IKWl5/PrHv/zttvvcXefXvp2KEjv/rVr7n1ltu8LXJlgdZ6P2v4POWLnh8El9ONtVT7fwTw3nvzmTjhXMJDoijM8+wfNmQEnTt34Z133uW8cycDMG/eewwbdjqdO3T1nqdD206MGDGSd9+dxx9+/yfctYxPcKKUDIqINDJhwUGeUUhb19wtxOlyk5ZTxJ7SlsQD2YWkZRWWjmpaRFpWIQdzCvEhX6TY5fZ2ST2WqoljYlQoCdGhnveoMNz57rp+VBGv/VuzfEoGrbUU5paQmZZP1sF8cjOLyDtShMNhGDiuC7GtIxog2lpic1ucTjeuYjclxS6cxS6cxW6cJe7SZVeFZU8ZV4mbtK1uVqRvweVy43Ja3E43LmfpsqvCcoXtAEldY0jsGI3L5cbtsqWv0uWyY92l252Vy7gqli19ucrKuC3WZQmPDqFN1xiCghyebe7yfRXX3a7SpMhVed3tqnBMWRxVz2MpP87tiWH9W0u9SVptSUVLsPOz7wE47cp4sg4VBDga//jjb29k/94DfPvDN7z8/GsAhIWGkp3h+dIyN6OQI2meOYCdxW62bd3GLbfcwrQ/3ExUVDSz/vE3fnbhz5hwzrm4nC4e++fTbNqyib//Ywaxka342+2zPMc6nVx0+fls2ryBG/90C3169+Xrb1dz//33sX/PQe656z5vTE5n7Y9tVBQcXHPaVFLkwrpstf9HBYUF7Nixnct/fiXZhyvvO6nHyaz9cS1ZBz2fde1Pa5k4blK1c/RMPokFH84j61ABIZE+hVlnSgZFRJqY4CCHZ1CZVrX/0etyWw7nFnGgLEksfS9PGj3LeUcZ5KaqYyWOvzkltM6fRVqmtsmxnHZeMt9/voPi0kdtigtq/oPM5XSzb/MR9m85wv6tWRzalUNRfs1l92w6wthr+hAVF1beSlNHbrelKK+EgpwSCnKLPe85xRTll1BU4KK4wElxodPzXuCkuNCzrajASUmh7/+eqjq8bledj0nfm3vc1/NFTkYhh3bV/CxUfaqhXUuaieSu3UlMaE1oSChDB5/m3f7FymU1ls88ksH773xMcpduALjdbq757eXs2r2Tt/47H4Axo89h1VdfMP+Ded5k8N35b/Hl6pXMe/1DRpx+BgCjzkgB4JHH/8kffz+NNq3bANCxZyK+SNuRVafPmpV1BGstreJaVdvXqlU8W7Zt9q4fycokrpZyR7KO1Om6daVkUESkGQpyGNrGep4VHHCUcgXFLg7nFnEot4hDOZ7nEA/nFHMot5DDOcXefYdzio6ZOMaGqQup+KZtt1jadoslLWcnuz73/OFfVCUZTN+byw9L97D1m4O1Jn9VZe7P460H1hAU7GDKtIGERYbgLHHRpkuMt0taQU4xmWn55GYUkpNRSE5GkXe5IKeYwtySY3cZk4ZlPF0KjaP03btecbnCPocBAw6Hp/ujL+Vr3lbaldEYHA7Kz1XxPMbTeJl1sABnsavC9TyBm9LeiGVlK+73nLu8u6SnrGd/ZmYmiYkJYCA41BAaEVzxdvDK31Y15P+Bo7rsb8MAT1xl/+Uo1YEJMqWfyzPImTEQHOq5UcGhDkLCSrc7oEf3nvTq1dN7zj59egNwdso5BAWXP2N8Uq+T+fiThTiCPf8vly77hM6dujBi+Ai8zcsGxp49jgcevpfvfvyacyecB8CSjz4/VshgICQ8uMYyjiADDir9PwK8nyMkLJiwyIr7PD9PxhhCw8u3B4cEedYrXMQR5FkJjQgGh2+/B+tKyaCISAsWERrknVfxWMoSx4OlSWPZM4npucVk5BWRFFHzyKcitQmq0Jhc1jJYmFfC8jc2s/HLA7UeFxwWRKukCFq1jSQ2MZxvFlVuVXM53bzz8DeVtrXvGUfm/nzvczv1JSjYQXCog+AQB8GhQaUvR/l7SBAhFdaDgh3s2rOLnr264wgyBAU7Sl8Gh3fZgSO4dF+QZ/nQrhwO7871JD1BxrM9yHhfQcFl647S/eXL3u3BVbc7cDg8+91uy8Ed2RTmlWAcxru9bNmUrnuXHQYTZHCY0vey8qb8OOOg1vOUlVu+fBmjR4+GKolfS5OamkpKykAA1q9fT6ukeuoj6AeJHaLrVD4sIhiHw5DQvvxRiJiECO97fDvP9uDQIOITWnnXAdoUtAKgU3I7EjuWXzcuMRqn00mrthEEBweTlZPB7j27aNs1vsYYim2e9/qj24zwKe7auomGhns+T9X/R6HRHTDGUFCSS1ybyvtyC3Jo3SaRVm092+Pj48kvyvGulykoziU+Pp5WSZG1jlh6opQMioiIT46VOKampjZsQNLkOULKl4sLnKTvy2XBE9+Td6TyiLrRCWEk929N+55xtOsWR0xieKUEoc/IDsx/4jty0mt/9nX/Ft+7eIVFBhMeHUJkTCjh0SFExIYSHhVCWEQwoRHBhIYHlb6Xrkd4vtEPDQ86rhFRC1N3Mzila52OadM5ps7XqauKf6w3FEewIShEo8rKiUlMTKRbt2688cYbNe5PTk72LoeEhNRYpipbxy4DkZGRJCcns3bt2mr71q1b5/nSo1S/fv1qLVff020oGRQREZGACKrwN1jmgXzee+w7CrKLvduST23NoHFdaN8jztP1rxat2kZy9X0jsdbyxZtb+P7T3bWWDQ51EN8uitjW4UQnhBMTH05MQjjRCWFExYURHh1SqfuZSGP0h+fODnQIxy0sLIyCgvodEGfixIm8/fbbREdH07t376OWXb16db3FMXnyZObOnUtWVhZxcXGAZ8L6nTt3Mnny5Erl/vrXv7Jt2za6d+8OwI4dO/jiiy944IEH6i0+UDIoIiIiAeKo8oV8WSIYEh7EOVP70n1gmzqdzxjDGZf05NSxnbwtdRu/PEBRvpP49lHEt4skJj78qImliNSvvn37kpGRwbPPPsvQoUMJDw/3+zWuuOIKXnrpJcaOHcvNN9/MgAEDKC4uZuvWrcyfP5958+YRGenp5TJ06NA6n//QoUN89tlnAOzatYv8/HzeeustwPP5ylrzpk+fziuvvMLkyZO5/fbbycrK4pZbbmHYsGFceOGF3vP99re/5amnnmLKlCnce++9GGOYMWMGnTt35ne/+92J3o6jCkgyaIy5AZgOtAfWAtOstTUPI+Qp3x94ChgGZAD/BmbZurbXioiISKPhCPYMxmErzIMSHOLggj8OoH3PVsd1TmMMsYnlI+32GdnhRMMUET+69tprWbVqFXfccQdHjhyha9euzJkzx6/XCAkJYdGiRTzwwAPMnj2b7du3ExUVRY8ePTjvvPMIDT2x0a/Xrl3LJZdcUmlb2frdd9/NzJkzAejYsSNLly7lpptu4uKLLyY0NJQpU6bwyCOPVJjrEKKiovj000+58cYbueqqq7DWMnbsWB577DGio+v2TGZdNXgyaIy5FHgcuAFYXvr+kTGmr7W22rjKxphYYDHwOXAacDIwB8gDHmmgsEVERMTPjDEEBYOzuDwZHP6zHsedCIpI4xcVFcWrr75abXvVNp7U1NRqg6YkJyfX+OzezJkzvQlYmfDw8Bq3+0NKSorPzxD279+fxYsXH7Ncly5dePvtt080tDoLRKf4m4A51trnrbXrrbV/AvYD19dS/gogErjGWvuTtfZt4J/ATaYlDi8lIiLSjDiL3ZXWTxndMUCRiIi0PA3aMmiMCQWGAA9X2fUxMLKWw0YAy6y1FZ80XQTMApKB7X4Os5KB06ZBq1aVN/7iF3DDDZCfD5MmVT9o6lTP6/Bh+PnPq++//nq49FLYvRuuuqr6/ptvhgsugI0boaZ+wnfdBeecA999B9OmVd9///0wciSsWAF33FF9/2OPwcCBsGQJ3Htv9f3//jecfDIsWACP1ND4+n//B507w+uvw7PPVt//1lvQujXMmeN5lRp45IjnXn74IURGwjPPQE2jPJWNSPjww/D++5X3RUTARx95lmfNgk8+qbw/MRHKvlW5/XZYubLy/k6d4JVXPMvTpnnuYUUnnQSzZ3uWr7sONm2qvH/gQM/9A7jyStizp/L+ESPgH//wLF98MaSnV94/dizMmOFZPvdcqPoA9fnnw1//6llOSaGa0p89R2Fhzfv1s1fjz55XDT973p9L0M+eDz97tf7e00ii4gc9BrXR4C0iIg2oobuJtgaCgLQq29OAc2o5ph1Q5a8e7/HtqJIMGmOuA64DaNu27QkPdd7f5eLIkSOVth3ctIl9qak4Cgs5tco+gAMbNnAgNZWQrCz61bB/79q1HEpNJezgQfrUsH/3jz+SHhNDxK5dnFzD/p3ff09mcDDRW7bQs4b92775huziYmL/v717D9aivu84/v4EjlyqGAUBBYP3hohWG5KoA8hYmURNoonTWJvU0ElsRk0vMTNJrE2lzXhrEypNYnPrlKZjRmOTeGu9xMuRNIIJWBQISLmJcOTmJXDkwEH89o/fHliW59yfc56Hs5/XzM7h2d9vn/09X37P7n6f3f3t0qWcVKF81cKFNL/xBkc9/zwTKpS/+OyztLzyCiOXLOH4CuXL589n9+rVHLNsGeMqlC/75S/Zc+SRjF2xgrG58r1ZLF+YN4+3hw7luJUrGV1h+cXZ/9nxq1czslC+t6WFJVn5hLVrOapQvuftt1mWlZ+4fj1HFsp3NzSwPCs/ZcMGDi+U72xqYmVWflpTE8ML5c0bNrAqK5+4eTNDCuW/Xb+etVn56Vu30rB9+wHlr69dy0tZ+Rmvvcag3QcOn/7q6tW8nJWfVSE2bX1vZ3PzQf0S3Pfa63ttKvW9vbnvuPteKu+o77W33Vvc2Ehzc7MfL2HdNun8cSx9eiODhwzivMtPqXVzzMxKRf05Bouk44CNwLT8gDGSbgKujIiDxn6V9BjwckR8JjdvArAOODciFrS3vsmTJ8fChQt71eb00M/pvXoPSxzL6nEsq8exrJ7exFLSoojo/pBuVlO93c82NjYy5byprHpuC6PGH86o8X3/7Lx65O3Qfo5Fko/D8uXLmThxYm0bVCM7duzgiCPKuV0o6igWnfWRjvax/X1mcBuwl3RGL280B58tbLOpnfp0sIyZmZkdAgYfNoh3n3NsrZthZlZK/XphfkS0AouAGYWiGcAz7Sw2H5gqaWihfhPp7KCZmZmZmZl1Uy3u0p4NzJT0WUkTJc0BjgO+AyDpVkn5kRl+BOwE5kqaJOnjwFeA2X7OoJmZmZkNdD7ktfb0tm/0+3MGI+IeSSOBvyE9dH4pcHFEvJRVORY4OVf/t5JmAN8GFgKvk54vOLtfG25mZmZm1s8aGhpoaWlh+PDhtW6K1aGWlhYaGhp6vHy/J4MAEXEncGc7ZTMrzFsCTOvjZpmZmZmZ1ZXRo0ezceNGxo0bx7Bhw/Bjtg3SGcGWlhY2btzImDFjevw+NUkGzczMzMyscyNGjACgqamJPXv21Lg1/WvXrl0MHTq084olUCkWDQ0NjBkzZl8f6Qkng2ZmZmZmdWzEiBG9OuA/VDU2NnL22WfXuhl1oa9iUYsBZMzMzMzMzKzGnAyamZmZmZmVkJNBMzMzMzOzEnIyaGZmZmZmVkJOBs3MzMzMzErIyaCZmZmZmVkJKSJq3YY+I2kr8FIv32YUsK0KzTHHspocy+pxLKunN7GcEBHHVLMx1veqsJ/19y9xHPZzLBLHIXEc9uuTfeyATgarQdLCiJhc63YMBI5l9TiW1eNYVo9jad3lPpM4Dvs5FonjkDgO+/VVLHyZqJmZmZmZWQk5GTQzMzMzMyshJ4Od+16tGzCAOJbV41hWj2NZPY6ldZf7TOI47OdYJI5D4jjs1yex8D2DZmZmZmZmJeQzg2ZmZmZmZiXkZNDMzMzMzKyEnAy2Q9K1ktZK2iVpkaSptW5TvZM0S1IUpk25cmV1miS1SGqUdHot21wvJE2T9ICkjVncZhbKO42dpCGSvilpm6Q3s/cb368fpE50IZ5zK/TVBYU6pY+npBsk/VrSdklbJT0oaVKhjvum9chA3896u+5tSJ6k6yS9kMViu6T5ki7JlZciDkWS/jr7fnwrN2/Ax0JVOGauVgycDFYg6QpgDnALcDbwDPCwpHfVtGGHhheBY3PTGbmyLwFfBP4ceB+wBfi5pCP6u5F16HBgKfCXQEuF8q7E7g7gcuBKYCowAnhI0qC+a3bd6iyeAI9zYF+9uFB+B47ndOBO4DzgAuAt4HFJR+fquG9at5VkP+vturcheRuALwO/D0wGngTuk3RmVl6WOOwj6RzgauCFQlFZYtHbY+Y7qEYMIsJTYQKeBb5fmPd/wK21bls9T8AsYGk7ZQJeAW7MzRsG7AA+V+u219MENAMzuxM74EigFfhkrs7xwNvAB2v9meopntm8ucBDHSzjeFaOy+HAXuAj2Wv3TU897Uul2s96u76v/d6GHBiP14DPlTEO2edZTfqRoBH4Vpn6BL08Zq5mDHxmsEDSYcB7gccKRY+Rftmyjp2UXRKzVtLdkk7K5p8IjCUX14hoAebhuHamK7F7L9BQqPMysBzHtz1TJG2RtFLS9yWNzpU5npUdQbqi5PXstfumdZv3s0B5vzvehgCSBkn6I1Jy/AzljMP3gP+MiCcL88sUi94cM1ctBk4GDzYKGARsLszfTPqPsfY9C8wELiKd9h8LPCNpJPtj57h2X1diN5b0a+u2DurYfo8AVwF/QLoM4/3Ak5KGZOWOZ2VzgMXA/Oy1+6b1hPez5f3ulHobIukMSc3AbuA7wMciYgnli8PVwCnAVysUlyUWvT1mrloMBnencskUH8CoCvMsJyIezr9WGpBjDfBpoG1wDse153oSO8e3goi4O/dyiaRFwEvAJcBPO1i0tPGUNBuYAkyJiL2FYvdN6wnvD0r03fE2BEj3iJ0FvJN0r9e/S5qeKx/wcZD0u6R7hadGRGsHVQd0LPrwmLnbMfCZwYNtI2Xaxax6NAdn6NaBiGgGlgGnAm0jJDmu3deV2G0i/dI+qoM61o6IaCLd3H9qNsvxzJH0T6Qb1C+IiDW5IvdN6wnvZ0v23fE2JImI1ohYFRELI+IG0lnSL1CuOJxL+gxLJb0l6S3gfODa7N+vZvXKEIt9enDMXLUYOBksyH6lWATMKBTNIF3XbV0kaSjwbtJNsGtJHXdGoXwqjmtnuhK7RcCeQp3xwEQc305JGgWMI/VVcDz3kTQH+GPSQdyKQrH7pnWb97NAib473oZ06B3AEMoVh/tIo2aelZsWAndn/15JeWKxTw+OmasXg1qPplOPE3AFaYSez2ZBnUMaCWxCrdtWzxPwddKvOycCHwAeAra3xY00pPJ24OPAJNIXvwk4otZtr/VEuon8rGzaCfxt9u93dTV2wL8AG4ELSUO1P0X61XFQrT9fPcUzK/s66dfJE0hDn88nnRl0PA+M47ezfncB6RfKtunwXB33TU896VsDfj/r7bq3IYVY3EY6mD+BlAzdShr58aIyxaGd2DSSjSZallhQhWPmasWg5sGo1wm4FlhHusl3ETCt1m2q9ynXUVuzzvkT4D25cpGG0n0F2AU8DUyqdbvrYSIlJFFhmtvV2AFDgW+SLrHYCTwIHF/rz1Zv8SQNz/wo6Zk9raR7BecWY+V4Bu3EMIBZuTrum556NA30/ay3696GFD7D3Gx/szvb/zxO7hEAZYlDO7Fp5MBkcMDHgiocM1crBsrezMzMzMzMzErE9wyamZmZmZmVkJNBMzMzMzOzEnIyaGZmZmZmVkJOBs3MzMzMzErIyaCZmZmZmVkJORk0MzMzMzMrISeDZjmSLpN0fYX50yWFpOn936rukzRL0gW1boeZmVm9kjQz27eHpNMqlE/PlV9Yizaa9TUng2YHugw4KBkEngPOzf4eCm4CnAyamZl1bgfwJxXmX5WVmQ1YTgbNuiAitkfEgojYXuu2mJmZWVX9FPiUJLXNkDQMuBz4Sc1aZdYPnAyaZSTNBT4NjMtdFrIuKzvoMlFJjZL+R9KHJC2W1CLpfyV9QNJgSbdIekXSa5LmSvqdwvqGS7pd0lpJrdnfGyV1+L3M3vtrklZL2iVpW9aOKVl5ZFVvzH2OWbnlz5f0hKQdkt6U9KikSYV1tH22SyUtlbRb0gpJnyjUO03SzyRtydqyXtK9kgZ3L/pmZmY18x/ABGBKbt7HgEE4GbQBzgdsZvt9DTgGeB/w0Wze7k6WOQX4R+BmoBn4B+CBbBoMzAQmZnW2AF+ClNABjwLvyda7BDgH+CpwNPDFDtb5ZeALwI3AYmAEMDlbDtLlrPOBucB3s3kbsvVeAtwP/Bfwqdz7/ULSmRHxcuGz/TMwK2v7NcDdkrZGxFNZnYeAN7KybcA44GL8Q5OZmR06XgLmkS4V/UU27yrgZ6R9u9mA5WTQLBMRqyVtBVojYkEXFxsJnBcRawCys3r3AydGRNvN5o9Kmgb8IVkyCFxJ+gXy/IiYl817IrtC5SZJt0fElnbWeS7wWETMyc17MPc5FmTvs7HC55gDPB0Rl7bNkPQUsIaUgP5Vru4Y4Ny295D0CLAM+HtgqqRRwKnApRHxQG65H7XTbjMzs3r1Q+Abkv4COAq4ELiotk0y63v+9d6sd1a2JYKZFdnfRwv1VgDjc/cjfIj0S+Qz2WWfg7OzhY8BDaSzhO35NXCxpJslTZF0WFcaKulU4GTgrsI6d5LOJE4rLPJyPpmMiL3AvcD7s6T3VVISeZukq7P3NzMzOxTdCwwBPgJ8EtgEPFHTFpn1AyeDZr3zeuF1awfzB5PuPwAYTbo/YU9h+lVWPrKDdd5CGi30o6TLWV6V9G/ZmbqOjM7+/muF9X64wjo3V3iPzcBhwDEREcAMYCFwK7BS0hpJ13TSDjMzs7oSETuA+0iXil4F3BURb9e0UWb9wJeJmtXGq8Ba4BPtlK9rb8GI2APcDtwuaSwpkZsNDAeu6GSdADcAj1coby28HlOhzpis3tasLWuAq7Iznr8HfB64U9K6iHi4g7aYmZnVmx+S7ql/B+l2DrMBz8mg2YF2A8P6YT2PkIasbo6IFZ1Vbk9EbAJ+IOliID8iaCsHf44XSUnm6RFxWxfe/nhJ5+TuGRxEuu/xV8VfS7OzhIslXQ98JmuLk0EzMzuU/Bz4MfBGRCyrdWPM+oOTQbMD/QY4OrvUcSGwKyKW9MF67gL+lDRozDeA50mXX55MuvzzsojYWWlBSfdn9Z8jXY56NukexO/mqv0GuCQb9OV1oCkimiRdB9yf3Wf4Y9IIoGOA84D1ETE79x6bgXsk3UQ6E3gNcFr2F0lnkgakuQdYRboEdibwFvBkjyNjZmZWA9m98T4jaKXiZNDsQD8gDd5yC/BO0iAvJ1R7JRGxR9IHga8AfwacCLwJrCZdolK8ZDNvHukM3XWkS0PXkx5pcXOuzudJj4V4kHRD/N8BsyLiv7ORTW8kfdZhpJvkF5CSurxV2fveQho1dB1wZe6xEpuydV8PjAd2kR6R8eGIWNT1aJiZmZlZLShd3WVmtp+kRmBwREzprK6ZmZmZHZo8mqiZmZmZmVkJORk0MzMzMzMrIV8mamZmZmZmVkI+M2hmZmZmZlZCTgbNzMzMzMxKyMmgmZmZmZlZCTkZNDMzMzMzKyEng2ZmZmZmZiXkZNDMzMzMzKyE/h9pkuhrH7UNcwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1080x432 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "figure saved as file  ./noise_via_chi_matrix/study_amplitude-damping-channel_q4_t200_s500.png\n"
     ]
    }
   ],
   "source": [
    "#channel_type = 'depolarizing'\n",
    "#channel_type = 'dephasing'\n",
    "channel_type = 'amplitude-damping'\n",
    "\n",
    "filename = file_df1.replace('TYPE', channel_type)\n",
    "filename = filename.replace('QUBITS', str(num_qubits))\n",
    "df = pd.read_pickle(filename)\n",
    "print(df.head)\n",
    "\n",
    "fig, axs = plt.subplots(1,2, sharex=False, figsize=(15, 6))\n",
    "fig.subplots_adjust(wspace=0.3)\n",
    "axs[0].tick_params(axis='both', labelsize=14)\n",
    "axs[1].tick_params(axis='both', labelsize=14)\n",
    "\n",
    "dfs = df.groupby(['t'])['sq_overlap'].mean().reset_index()\n",
    "axs[0].plot(dfs.t, dfs.sq_overlap, linewidth=3, label=channel_type)\n",
    "if channel_type == 'depolarizing':\n",
    "    #axs[0].plot(dfs.t, 1/2+1/2*np.exp(-rate*dfs.t))\n",
    "    axs[0].plot(dfs.t, [1/2**num_qubits]*len(dfs.t), 'r--')\n",
    "elif channel_type == 'dephasing':\n",
    "    #axs[0].plot(dfs.t, 1/4+3/4*np.exp(-p*dfs.t))\n",
    "    axs[0].plot(dfs.t, [1/2**(num_qubits/2)]*len(dfs.t), 'r--')\n",
    "elif channel_type == 'amplitude-damping':\n",
    "    axs[0].plot(dfs.t, [0]*len(dfs.t), 'r--')\n",
    "\n",
    "# Fix time step, compute the convergence when average is over more ensemble states.\n",
    "for t in [2, 10, 20, 50, 100]:\n",
    "    # If we do not take a copy, pandas would complain that:\n",
    "    # 'A value is trying to be set on a copy of a slice from a DataFrame.'\n",
    "    dft = df.loc[(df['t'] == t)].copy()\n",
    "    a = dft['sq_overlap'].to_numpy()\n",
    "    a = np.cumsum(a)\n",
    "    b = np.ones(num_ensemble_states)\n",
    "    b = np.cumsum(b)\n",
    "    c = np.divide(a, b)\n",
    "    dft['avg_sq_overlap'] = c\n",
    "    axs[1].plot(dft.s, dft.avg_sq_overlap, linewidth=3, label='time={}'.format(t))\n",
    "\n",
    "axs[0].set_ylabel('squared overlap', fontsize=16)\n",
    "axs[0].set_xlabel('time steps', fontsize=16)\n",
    "axs[0].grid()\n",
    "axs[0].legend(loc='upper right', fontsize=16)\n",
    "axs[1].set_xlabel('num states in the average', fontsize=16)\n",
    "axs[1].set_xlabel('M', fontsize=16)\n",
    "axs[1].set_ylabel('squared overlap', fontsize=16)\n",
    "axs[1].legend(loc='lower right', fontsize=16)\n",
    "axs[1].grid()\n",
    "\n",
    "filename = data_dir + 'study_{}-channel_q{}_t{}_s{}.png'.format(\n",
    "    channel_type, num_qubits, num_time_steps, num_ensemble_states)\n",
    "fig.savefig(filename)\n",
    "plt.show()\n",
    "print('figure saved as file ', filename)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### UTILITY CHECKS\n",
    "\n",
    "Confirm the correct behavior of ChiMatrix methods by considering the 2x2 symmetric matrix:\n",
    "$$M = \\begin{bmatrix} 1 & 3 \\\\ 3 & 7 \\end{bmatrix}$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 80,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Original matrix:\n",
      " [[1 3]\n",
      " [3 7]] \n",
      "\n",
      "Eigenvalues: [-0.24264069  8.24264069]\n",
      "Eigenvectors: [-0.92387953  0.38268343] and [-0.38268343 -0.92387953] \n",
      "\n",
      "Before standardization:  M =?=\n",
      " [[1. 3.]\n",
      " [3. 7.]] \n",
      "\n",
      "After standardization, eigenvectors: [[-0.92387953  0.38268343]] and [[-0.38268343 -0.92387953]]\n",
      "After standardization:  M =?=\n",
      " [[1. 3.]\n",
      " [3. 7.]] \n",
      "\n",
      "After normalization, eigenvectors: [[-2.69121547  1.11473795]] and [[-1.11473795 -2.69121547]]\n"
     ]
    }
   ],
   "source": [
    "from numpy import linalg as la\n",
    "\n",
    "M = np.array([[1,3],[3,7]])\n",
    "print('Original matrix:\\n', M, '\\n')\n",
    "w, v = la.eig(M)\n",
    "print('Eigenvalues:', w)\n",
    "print('Eigenvectors:', v[:,0], 'and', v[:,1], '\\n')\n",
    "\n",
    "# Eigenvectors need to be 'standardized' so that: M = w0 |v0><v0| + w1 |v1><v1|\n",
    "v0 = np.array(v[:,0]); v0.shape = (2,1)\n",
    "v1 = np.array(v[:,1]); v1.shape = (2,1)\n",
    "assert v0.conj().transpose().dot(v1) == 0\n",
    "M_current = w[0] * v0.dot(v0.conj().transpose()) + w[1] * v1.dot(v1.conj().transpose())\n",
    "# Is this currently the case?\n",
    "print('Before standardization:  M =?=\\n', M_current, '\\n')\n",
    "\n",
    "# If not, standardize the eigenvectors.\n",
    "# (this is done automatically inside class ChiMatrix when method SolveEigensystem is called)\n",
    "a0 = v0.conj().transpose().dot(M.dot(v0))\n",
    "s0 = v0 * np.sqrt(abs(w[0]/a0))\n",
    "a1 = v1.conj().transpose().dot(M.dot(v1))\n",
    "s1 = v1 * np.sqrt(abs(w[1]/a1))\n",
    "M_current = w[0] * s0.dot(s0.conj().transpose()) + w[1] * s1.dot(s1.conj().transpose())\n",
    "np.set_printoptions(linewidth=150)\n",
    "print('After standardization, eigenvectors:', s0.transpose(), 'and', s1.transpose())\n",
    "print('After standardization:  M =?=\\n', M_current, '\\n')\n",
    "\n",
    "# Normalization of the eigenvectors.\n",
    "p_norm = abs(w[0]) + abs(w[1])\n",
    "p0 = abs(w[0]) / p_norm\n",
    "p1 = abs(w[1]) / p_norm\n",
    "n0 = s0 * np.sqrt(p_norm)\n",
    "n1 = s1 * np.sqrt(p_norm)\n",
    "print('After normalization, eigenvectors:', n0.transpose(), 'and', n1.transpose())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "----\n",
    "## END\n",
    "----"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}