{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0191bb28-28d7-44e1-98c4-5c6cfe18b6f9",
   "metadata": {},
   "source": [
    "# Noisy Circuit Simulation\n",
    "## Overview\n",
    "In this notebook, we will show how to perform noisy circuit simulation using TenCirChem. The hydrogen molecule $\\rm{H}_2$ is used for the demonstration. Using parity transformation, the system corresponds to 2 qubits.\n",
    "We will firstly illustrate the workflow step by step, and then use well-wrapped high-level interfaces to complete the same task with much fewer lines of code.\n",
    "We will also show it is possible to use Qiskit ansatz directly in TenCirChem.\n",
    "In the end, we study the effect of gate noise and the number of measurement shots on the accuracy and uncertainty of VQE."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a4fb5490-5cb1-4bc7-bdf5-c245e39149f8",
   "metadata": {},
   "source": [
    "## Setup\n",
    "TenCirChem uses the `HEA` class for noisy circuit simulation, which has a different set of interfaces from the UCC classes. \n",
    "A primary reason is that *all* ideal UCC circuits are too large for both noisy circuit simulation and execution on state of the art quantum devices.\n",
    "Consequently, usually hardware-efficient ansatz (HEA) have to be adopted for such tasks.\n",
    "Another reason is that inherently the workflow for the noiseless UCC circuit simulation and noisy HEA circuit simulation are very different."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "db2a228b-9150-47cc-815c-07ebec4aa3b5",
   "metadata": {},
   "outputs": [],
   "source": [
    "from tencirchem import HEA\n",
    "from tencirchem.molecule import h2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a0401605-86d4-43ba-8c6f-9c19d234aaee",
   "metadata": {},
   "source": [
    "## Define the Hamiltonian\n",
    "The first difference between `UCC` and `HEA` is that the `HEA` class is not initialized by molecular inputs.\n",
    "Rather, `HEA` accepts the Hamiltonian in `openfermion.QubitOperator` form and the circuit ansatz for inputs.\n",
    "To facilitate defining the Hamiltonian, the `UCC` class provides several useful shortcuts via class attributes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "0e8345ef-a98a-4d6a-b5fc-e252dcfbb158",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.7141392859919029 [] +\n",
       "-1.2527052599711868 [0^ 0] +\n",
       "-0.48227117798977825 [0^ 1^ 0 1] +\n",
       "-0.6745650967143663 [0^ 2^ 0 2] +\n",
       "-0.18126641677772592 [0^ 2^ 1 3] +\n",
       "-0.6635375947675042 [0^ 3^ 0 3] +\n",
       "-0.18126641677772592 [0^ 3^ 1 2] +\n",
       "-0.47569770336145906 [1^ 1] +\n",
       "-0.18126641677772592 [1^ 2^ 0 3] +\n",
       "-0.6635375947675038 [1^ 2^ 1 2] +\n",
       "-0.18126641677772592 [1^ 3^ 0 2] +\n",
       "-0.6974673850129379 [1^ 3^ 1 3] +\n",
       "-1.2527052599711868 [2^ 2] +\n",
       "-0.48227117798977825 [2^ 3^ 2 3] +\n",
       "-0.47569770336145906 [3^ 3]"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from tencirchem import UCCSD\n",
    "\n",
    "uccsd = UCCSD(h2)\n",
    "# Hamiltonian as openfermion.FermionOperator\n",
    "h_fermion_op = uccsd.h_fermion_op\n",
    "h_fermion_op"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "6a36370f-8eb9-46fb-b4a7-2e4a0fad24b3",
   "metadata": {},
   "outputs": [],
   "source": [
    "from tencirchem import parity"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "1b1eb25d-35ad-42f2-b324-7d51c357e6a9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.33948675952516477 [] +\n",
       "0.18126641677772592 [X0 X1] +\n",
       "-0.39422935037950685 [Z0] +\n",
       "-0.01123932304807404 [Z0 Z1] +\n",
       "0.3942293503795067 [Z1]"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# use parity transformation for qubit Hamiltonian and remove 2 qubits\n",
    "h_qubit_op = parity(h_fermion_op, uccsd.n_qubits, uccsd.n_elec)\n",
    "h_qubit_op"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7de06ff3-aabf-4b41-a903-4d6f92d81fc1",
   "metadata": {},
   "source": [
    "An alternative approach is to use openfermion tools such as `openfermion.MolecularData`.\n",
    "\n",
    "## Define the ansatz\n",
    "Next, define the ansatz using TensorCircuit. A 1-layer $R_y$ ansatz composed of 4 parameters is constructed. A TensorCircuit ansatz is a function that has one argument (the circuit parameters) and returns the parameterized circuit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "f3a77d6c-5a76-41d8-9eaa-bf1a228b3c57",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAATEAAACuCAYAAABeIjpKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAUUElEQVR4nO3dfVRVdb7H8Tf4BPKgIDgQqPgAoSg4anh1rInCmUgtszFtzKnGcpoya3KkP2ZVM3dWdUm7a2XNjLqaGWu6MczV6irOjFbkGrOHwdBylJVCoh7gpEcQFZ9Azv3jjCSJeQ4czua3z+e1lgvZZz98Zf/4uPdv//beIW63242IiKFCrS5ARKQzFGIiYjSFmIgYTSEmIkZTiImI0RRiImI0hZiIGE0hJiJGU4iJiNEUYiJiNIWYiBhNISYiRlOIiYjRFGIiYjSFmIgYTSEmIkZTiImI0RRiImI0hZiIGE0hJiJGU4iJiNEUYiJiNIWYiBhNISYiRlOIiYjRFGIiYjSFmIgYTSEmIkZTiImI0RRiImI0hZiIGE0hJiJGU4iJiNEUYiJiNIWYiBhNISYiRlOIiYjRelpdgFzK7YaWJqur8E1oLwgJsboK+1Ab8J5CrBtqaYL3VlhdhW9yFkOP3lZXYR9qA97T6aSIGE0hJiJGU4iJiNEUYiJiNIWYiBhNISYiRlOIiYjRNE7MRj6t3MLPV+a0mRbWO4Lk+DRyx81n5ncepkcP7XI7C8Y2YK9/jQCQM/ZOstNvxo2b+hNO3v7kVVZueIyDh8v52Q9WW12eBEAwtQGFmA2lJo0jd/xdrd/PmPwgC55L52//fJl7b3qa/pHxFlYngRBMbUB9YkEgvHcE6UP+A7fbTc3RSqvLEQvYuQ0oxIJE7b8bbnTfWIsrEavYtQ3odNKGzjSdoqHRhdvt6Q/Z8OFKKqp3kD4om+T4NKvLkwAIpjZg+xBzuVw899xzvPHGGzgcDuLj45k1axbPPPMMixcv5g9/+AMvvvgiixYtsrpUv3l181O8uvmpNtOmjJ7Fw7f9xqKKrNd8Hg4chcaz0KsHJPSDmAirq+o6wdQGbB1iO3fuJC8vD6fTSUREBKNGjaKmpoYVK1ZQWVlJXV0dAGPHjrW2UD+bNnEh12XOprmlif21uyjaUoCrwUHvXmGt8zz92lxa3C08Mf8vrdOOn6rj/uUZLJy+nBvHzbOidL9rOA3b9sKHFXDizFfTQ4BRSXDt1ZCeaFl5XSaY2oBt+8RcLhczZszA6XSyZMkSamtrKSsrw+l0UlBQwMaNGyktLSUkJITMzEyry/WrpLhUxqXlkp2ex5ycfH597wY+d5TywroHWud5eNZv2V21jZIdha3TXnzzITKGTjGm8V6Jow6e/yts/lfbAANwA7urYWUJrC/zPITQToKpDdg2xBYvXozD4WDRokUsX76cqKio1s/y8/PJysqiubmZlJQUoqOjLay062WkTCZ33Hy2fFrE7qoPAE/n7pLZv+eltxbhaqjhH5+t5bPKLTw6a6XF1frHkRPwuxI4fubK85aUw993dX1NVrJzG7BliJWXl1NUVERcXBzPPvtsu/OMHz8egKysrDbT9+/fzy233EJUVBQxMTH86Ec/4ujRo11ec1ebl/sEoaE9eGXTk63Trkm/ie9m3kFB4V28+MaDPDb7ZaIjBlhYpf8U7/T0f3lr8y44erLLyukW7NoGbBlihYWFtLS0MG/ePCIjI9udJzw8HGgbYidOnCAnJweHw0FhYSGrV69m69atTJ8+nZaWloDU3lWS4kaQkzWXHRXvsuuLra3TF85YTvXRCq5Jz2PiyGkWVug/Dadg1yHflnEDH+zrknK6Dbu2AVuGWElJCQA5OTmXncfhcABtQ2z16tVUV1fz1ltvMX36dGbPns3rr7/ORx99xPr167u26AC488ZfEBoSyiubv/qfOLx3BImxwxiaMMbCyvzrkypo6UAfV+kXfi+l27FjG7Dl1ckDBw4AMGTIkHY/b25uZtu2bUDbECsuLmbKlCkMHjy4ddqkSZMYNmwYGzZsYObMmR2qZ8KECTidTq/n790znNWLfD8syBp+PW8vu/xv75BvjWTTc+d9Xq83UtNSOdd8ukvW7ausGb8idcoCn5c7fgYGDU7B3dLcBVX5JtjaQEJCAtu3b+/QsrYMscbGRgBOn27/B1pUVITL5SIqKoqhQ4e2Tt+zZw+zZ8++ZP6MjAz27NnT4XqcTifV1dVezx/Wq2+Ht2WV2poazjSdsroMAIafPNHhZaurHbhbuuaX3BdqA96zZYglJCRQX19PWVkZkyZNavNZbW0tS5cuBSAzM5OQi16UV19fT//+/S9ZX2xsLJ9//nmn6vFF757hHd6WVRKvuqrbHImFNh/v0HKnjzu5KtG3fdVVgq0N+Po7cjFbhlhubi7l5eUUFBQwdepU0tI8t1mUlpYyf/58XC4XELhBrr4eJp8/F9h3Dj7/0y2dXse+vfu6zXsnG07Dr970vV/slkkJrPp3X6nV1Aa8Z8uO/fz8fAYMGMChQ4fIyMhgzJgxpKamkp2dzbBhw7jhhhuAS4dXxMTEcOzYsUvWV1dXR2ysvW6atbN+4ZA5yLdlQkJg0oiuqUe6li1DLDk5ma1btzJt2jTCwsKoqqoiNjaWVatWsXHjRvbu3QtcGmIjR45st+9rz549jBw5MiC1i3/M+DZEhl15vgtuGgOx7Y/GkW7OliEGnkAqLi7mxIkTnDhxgo8//piFCxfS2NhIVVUVoaGhjB49us0y06dP5/33328dfgHw8ccfU1lZyYwZMwL9T5BOGBAJD97gOSq7kqkZ8L3RV55Puifbhtjl7N69G7fbTWpqKn37tr0CtHDhQhITE7n11lspLi5m7dq13HnnnWRnZ3PrrbdaVLF01FUx8PObIS+z/TDLHAQP3QjTxnpOJ8VMQRdiu3Z5bpL7+qkkQHR0NCUlJSQmJjJ37lzuu+8+Jk+eTHFxMaGhQfejsoWoMPj+GHhyJjzyPejb+6vpP74OUrvHxUjpBFtenfwm3xRiAMOHD6e4uDiQJQWM48g+lhXdTUOji4iwfiyds4aUhAyrywqIHqEwNN7zLDGA0CA88rLr/g+6w4srhZidvbDuJ9w8cSFrHt/LnJzHWVZ0j9UlSQDZdf8HXYiVlJTgdruZNs28G107o/7kYfY6tpM7zvMGnGvH3M6RY4eodlVYXJkEgp33f9CFWLA6cuwQsdGJrS9ODQkJYWDMYA4fO2hxZRIIdt7/CjERMZpCLEjE9x9E3fFazp/3PKHB7XZzuP4gA/sPvsKSYgd23v8KsSAREzmQEUnjeKfsNQC27lpHXP9kkuJ0r00wsPP+D7ohFsHs0dtXsazoHgpLnqFvWDRL7/ij1SVJANl1/yvEgsiggVez4uEPrS5DLGLX/a/TSRExmkJMRIymEBMRoynERMRoCjERMZquTnZDob0gZ7HVVfgmtJfVFdiL2oD3FGLdUEgI3ealG2INtQHv6XRSRIymEBMRoynERMRoCjERMZpCTESMphATEaMpxETEaAoxETGaQkxEjKYQExGjKcRExGgKMRExmkJMRIymEBMRoynERMRoCjERMZpCTESMpie7dkNuN7Q0WV2Fb0J7eZ5GKv6hNuA9hVg31NIE762wugrf5CzW45T9SW3AezqdFBGjKcRExGgKMRExmkJMRIymEBMRo+nqpNja8dNwqM7z5+hJOHXOM/30OfioEgbFQkI/6KH/zo2lEBPbaToPnx6E9/dClav9ec6dhz9/5Pl7394wcTh8JxXiogJXp/iHQkxsw+2G7fvh/8rg5Fnvlzt1Dt4r9/z59hC4fQJEhnVdneJfCjEb+bRyCz9fmdNmWljvCJLj08gdN5+Z33mYHj3sucsbTsNfPobd1Z1bz44DsM8JP8iGsYP9U1sgBWMbsNe/RgDIGXsn2ek348ZN/Qknb3/yKis3PMbBw+X87AerrS7P75wN8Lt3PUHmDyfPwpqt8P0xcNMYM2+nCqY2oBCzodSkceSOv6v1+xmTH2TBc+n87Z8vc+9NT9M/Mt7C6vzr8HF46R04ecb/6960y/M1L9P/6+5qwdQGdE0mCIT3jiB9yH/gdrupOVppdTl+c7YJVr/XNQF2waZdnn4209m1DYCOxIJG7b8bbnTfWIsr8Z8NO8F10rdlHrsJosM9Qy/+++/eLfPGdkhNgH7hPpfYrdixDUCQHIm5XC7y8/MZMWIEYWFhDBo0iEceeYTGxkYWLFhASEgIL730ktVl+s2ZplM0NLo4dvII+2t3seKNh6io3kH6oGyS49OsLs8vKr70DKHwVXQ49O/r+eqtU+fgf//p+7asFAxt4ALbH4nt3LmTvLw8nE4nERERjBo1ipqaGlasWEFlZSV1dXUAjB071tpC/ejVzU/x6uan2kybMnoWD9/2G4sq8r8L/VWB8i8HOOog2ZCDmGBoAxfY+kjM5XIxY8YMnE4nS5Ysoba2lrKyMpxOJwUFBWzcuJHS0lJCQkLIzDSw9/Yypk1cSMH9b/P0gr9y380FRPWNxdXgoHevrwY/Pf3aXH79pzvaLHf8VB1z/jORd8v+J9Al++TLBtj3ZeC3u21f4LfZUXZvAxezdYgtXrwYh8PBokWLWL58OVFRXw3Hzs/PJysri+bmZlJSUoiOjrawUv9KiktlXFou2el5zMnJ59f3buBzRykvrHugdZ6HZ/2W3VXbKNlR2DrtxTcfImPoFG4cN8+Ksr32QYU12/1kP5wx5Gmrdm8DF7NtiJWXl1NUVERcXBzPPvtsu/OMHz8egKysrNZpF0IvOzubPn36EGLiIKGvyUiZTO64+Wz5tIjdVR8Ans7dJbN/z0tvLcLVUMM/PlvLZ5VbeHTWSourvbIKC47CwHOr0sGj1my7s+zWBi5m2xArLCykpaWFefPmERkZ2e484eGe3t2LQ6yiooJ169aRkJDANddcE5BaA2Fe7hOEhvbglU1Ptk67Jv0mvpt5BwWFd/HiGw/y2OyXiY4YYGGVV9Z0HmqPWbd9R5112+4su7SBr7NtiJWUlACQk5Nz2XkcDgfQNsSuu+46amtrWb9+Pbm5uV1bZAAlxY0gJ2suOyreZdcXW1unL5yxnOqjFVyTnsfEkdMsrNA7tcegxW3d9g8ZHGJ2aQNfZ9urkwcOHABgyJAh7X7e3NzMtm3bgLYhFhrq/1yfMGECTqfT6/l79wxn9SL/9yLfeeMveG9nIa9sfpLlD7wHeAZBJsYOY2jCmE6tOzUtlXPNfrrv5xskXJ3DlB//qd3PLowB+ybRYV99/eVtl5/vcuPI/v7OP/jF3B96WW3HBVsbSEhIYPv27R1a1rYh1tjYCMDp0+3/UIuKinC5XERFRTF06NAurcXpdFJd7f2dyWG9+nZoO1nDr+ftZZc/TBnyrZFseu58h9Z9JbU1NZxpOtUl675Y7281XPazC2PAvBEa6v28F2tqdvu0LztKbcB7tg2xhIQE6uvrKSsrY9KkSW0+q62tZenSpQBkZmZ2eed9QkKCT/P37mne0PDEq64KyJFYTL/LP/DruBebjw7zBFhLCxz/htuVLreunj3cJCUlXXlDnRRsbcDX35GL2TbEcnNzKS8vp6CggKlTp5KW5hmlXFpayvz583G5PE/LC8QgV18Pk8+fM++dg/v27gvIOwe/bIBni9v/zJvbiH55m+cI7PgZ+OWbvm9/Zt71/PlXDt8X9JHagPdsG2L5+fm8/vrrHDp0iIyMDNLT0zlz5gwVFRXk5eWRkpLCpk2b2vSHBavnf7rF6hK8Fh8NfXrC2WZrtj/IkBH7vjKpDXydba9OJicns3XrVqZNm0ZYWBhVVVXExsayatUqNm7cyN69nhvvFGJmCQ2BpBjrtj/IrNEHQcG2R2IAI0eOpLj40nOPkydPUlVVRWhoKKNHj7agMumM0cnwxZHAb7dfuLUBKu2zdYhdzu7du3G73aSlpdG376VXgdauXQvAnj172nyfkpLChAkTAleotCt7GPz1U2huCex2J6XqrUjdUVCG2K5dnkcgXO5Ucvbs2e1+f/fdd7NmzZourU2uLDIMxg4J7MMKQ0Ng0vDAbU+8pxBrh9tt4ZBw8UpeJnx2CM4FqIP/hlHQr2NDt7oNx5F9LCu6m4ZGFxFh/Vg6Zw0pCRlWl9VpQXlwfKUQk+5vQCTc8u3AbCuhn+eFIaZ7Yd1PuHniQtY8vpc5OY+zrOgeq0vyi6A8ErtwX2WwOHn6GPc/P5qzTaeJ7zeIpvNncR79ghvHz2fJ7JetLq/DJqfCnmrYU+P9MhcGsXozMBagVw/44STo2cP3+rqT+pOH2evYzn/dvxmAa8fczktvLqLaVUFS3AiLq+ucoAyxYBMZ3p8bxv6Q8D5R3DX1CUo/30RhyTNGBxh4+qnuvhZWlsB+L69WevtcffB04t97LQy2wbCKI8cOERud2PrOyZCQEAbGDObwsYPGh1hQnk4Go4qanYxI8px/7XN8woirAnQu1sX69IQHcuDqRP+vd+H1MKrr7zCSTlKIBYkvvh5iSfYIMYA+veAn18PM8Z7Tv85KS4DHp/k/GK0U338QdcdrOX/ecyXE7XZzuP4gA/sb+Jrzr1GIBQFXQzWEhBDXz3NY8YXzs04/dqW7CQ2F69Mh/2YYk9yxt3YPiIQ5E+GnN0Bs+8/RNFZM5EBGJI3jnbLXANi6ax1x/ZONP5UE9YkFhYrqHW1OHyPD+rP+w98a3yfWnvhoWPBdqG+EDys8wzC+PA6XGzUT0QeGxXsuElyd6Olns6tHb1/FsqJ7KCx5hr5h0Sy9449Wl+QXIW4Niup2THyCQc5iLHmCgTfONkN1nedFu83nPR324b0hOQZiIjp21NbV1Aa8pyMxsb0+PWHYQM8fsR/1iYmI0RRiImI0hZiIGE0hJiJG09XJbsjthpYmq6vwTWiv7nmVz1RqA95TiImI0XQ6KSJGU4iJiNEUYiJiNIWYiBhNISYiRlOIiYjRFGIiYjSFmIgYTSEmIkZTiImI0RRiImI0hZiIGE0hJiJGU4iJiNEUYiJiNIWYiBhNISYiRlOIiYjRFGIiYjSFmIgYTSEmIkZTiImI0RRiImI0hZiIGE0hJiJG+3+C9Bi8ol4g7wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 370.906x200.667 with 1 Axes>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import numpy as np\n",
    "from tensorcircuit import Circuit\n",
    "\n",
    "# the ansatz\n",
    "n_qubits = 2\n",
    "init_guess = [0, np.pi, 0, 0]\n",
    "\n",
    "\n",
    "def get_circuit(params):\n",
    "    c = Circuit(n_qubits)\n",
    "    # the ansatz body\n",
    "    c.ry(0, theta=params[0])\n",
    "    c.ry(1, theta=params[1])\n",
    "    c.cnot(0, 1)\n",
    "    c.ry(0, theta=params[2])\n",
    "    c.ry(1, theta=params[3])\n",
    "    return c\n",
    "\n",
    "\n",
    "get_circuit(init_guess).draw(output=\"mpl\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4fbc9f10-10d2-4be5-a8dc-b71132eb66a3",
   "metadata": {},
   "source": [
    "## Noisy simulation with the `HEA` class\n",
    "\n",
    "The rest of the workflow is similar to `UCC`: create an `HEA` instance, run the kernel and print the summary.\n",
    "\n",
    "By using the `\"tensornetwork-noise\"` engine, a depolarizing noise with probability $p=0.02$ for the CNOT gate is added\n",
    "\n",
    "$$\n",
    "\\mathcal{E}(\\rho) = (1 - p)\\rho + \\frac{p}{15} \\sum_{i} P_j \\rho P_j\n",
    "$$\n",
    "\n",
    "Here $P_j$ are two qubit Pauli matrices except $I$. An alternative form is\n",
    "\n",
    "$$\n",
    "\\mathcal{E}(\\rho) = \\frac{16p}{15} \\frac{I}{4} + (1 - \\frac{16p}{15}) \\rho\n",
    "$$\n",
    "\n",
    "When $p=0$, no error is introduced. When $p=\\frac{15}{16}$, this is a completely depolarizing channel. If $\\rho$ is a pure state, the gate fidelity is $1-\\frac{4}{5}p$,\n",
    "which is is approximately 98% if $p=0.02$.\n",
    "\n",
    "The gradient of the parameters is evaluated by the parameter-shift rule and not by auto-differentiation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "8ac5b872-9e63-4698-add2-baf2ea43c8a3",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "############################### Circuit ###############################\n",
      "   #qubits  #gates  #CNOT  #multicontrol  depth #FLOP\n",
      "0        2       5      1              0      3    96\n",
      "######################### Optimization Result #########################\n",
      "            e: -1.1202549357480136\n",
      "          fun: array(-1.12025494)\n",
      "     hess_inv: <4x4 LbfgsInvHessProduct with dtype=float64>\n",
      "   init_guess: array([0.        , 3.14159265, 0.        , 0.        ])\n",
      "          jac: array([5.16941739e-10, 5.55111512e-17, 0.00000000e+00, 0.00000000e+00])\n",
      "      message: 'CONVERGENCE: NORM_OF_PROJECTED_GRADIENT_<=_PGTOL'\n",
      "         nfev: 5\n",
      "          nit: 3\n",
      "         njev: 5\n",
      "     opt_time: 0.25640320777893066\n",
      " staging_time: 1.1920928955078125e-06\n",
      "       status: 0\n",
      "      success: True\n",
      "            x: array([-2.25973122e-01,  3.14159265e+00, -6.82347069e-17,  7.19412730e-17])\n"
     ]
    }
   ],
   "source": [
    "hea = HEA(h_qubit_op, get_circuit, init_guess, engine=\"tensornetwork-noise\")\n",
    "hea.kernel()\n",
    "hea.print_summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "13e86d30-4039-4479-b259-fbb40211fced",
   "metadata": {},
   "source": [
    "Different noise models can be specified using the `NoiseConf` class from TensorCircuit. Here we show a customized noise model with $p=0.25$ and gate fidelity of 80%. The energy obtained is higher than the $p=0.02$ result."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "b1c42f33-1ba7-435c-93ac-76db2959b42b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "############################### Circuit ###############################\n",
      "   #qubits  #gates  #CNOT  #multicontrol  depth #FLOP\n",
      "0        2       5      1              0      3    96\n",
      "######################### Optimization Result #########################\n",
      "            e: -0.9245310332616318\n",
      "          fun: array(-0.92453103)\n",
      "     hess_inv: <4x4 LbfgsInvHessProduct with dtype=float64>\n",
      "   init_guess: array([0.        , 3.14159265, 0.        , 0.        ])\n",
      "          jac: array([ 3.87354315e-10, -5.55111512e-17, -5.55111512e-17,  0.00000000e+00])\n",
      "      message: 'CONVERGENCE: NORM_OF_PROJECTED_GRADIENT_<=_PGTOL'\n",
      "         nfev: 5\n",
      "          nit: 3\n",
      "         njev: 5\n",
      "     opt_time: 0.25928783416748047\n",
      " staging_time: 9.5367431640625e-07\n",
      "       status: 0\n",
      "      success: True\n",
      "            x: array([-2.25973122e-01,  3.14159265e+00, -1.83807579e-16, -9.76918418e-17])\n"
     ]
    }
   ],
   "source": [
    "from tensorcircuit.noisemodel import NoiseConf\n",
    "from tensorcircuit.channels import isotropicdepolarizingchannel\n",
    "\n",
    "engine_conf = NoiseConf()\n",
    "# larger noise, corresponding to 80% fidelity\n",
    "channel = isotropicdepolarizingchannel(p=0.25, num_qubits=2)\n",
    "engine_conf.add_noise(\"cnot\", channel)\n",
    "\n",
    "hea = HEA(h_qubit_op, get_circuit, init_guess, engine=\"tensornetwork-noise\", engine_conf=engine_conf)\n",
    "hea.kernel()\n",
    "hea.print_summary()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "48608036-ae47-4e77-9375-899b186ef088",
   "metadata": {},
   "source": [
    "The `\"tensornetwork-noise\"` engine does not consider measurement uncertainty. In other words, the engine assumes that infinite shots are taken for each term in the Hamiltonian.\n",
    "In order to take measurement uncertainty into account, use the `\"tensornetwork-noise&shot\"` engine. \n",
    "By default the number of shots for each term is 4096.\n",
    "Note that the engine is significantly slower to run due to the computational overhead.\n",
    "Similar to the `UCC` class, the `HEA` class supports switching engine at runtime.\n",
    "The following example shows that by increasing measurement shots the energy uncertainty is suppressed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "7c21a332-926b-430c-8675-d7cdcc6fe2bd",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "shots number: 4096\n",
      "-0.9187451689899536\n",
      "-0.9285959899718158\n",
      "-0.9285225235212169\n",
      "-0.9349106044977904\n",
      "-0.9264353453475235\n",
      "shots number: 4096 * 128\n",
      "-0.9250756148331946\n",
      "-0.9245772630966033\n",
      "-0.9238598336661094\n",
      "-0.924377692124025\n",
      "-0.923522568588401\n"
     ]
    }
   ],
   "source": [
    "# default shots is 4096\n",
    "print(\"shots number: 4096\")\n",
    "for i in range(5):\n",
    "    print(hea.energy(engine=\"tensornetwork-noise&shot\"))\n",
    "print(\"shots number: 4096 * 128\")\n",
    "hea.shots = 4096 * 128\n",
    "for i in range(5):\n",
    "    print(hea.energy(engine=\"tensornetwork-noise&shot\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "93ec8bbd-d955-44cb-9419-885a5cb2e789",
   "metadata": {},
   "source": [
    "Finally, noiseless results are available by using the `\"tensornetwork\"` engine, in which case the energy is in exact agreement with the FCI energy for the hydrogen molecule system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "093c63e4-6a5d-4cc3-8e48-965482c91a78",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-1.1372744055294381, -1.1372744055294384)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "hea.energy(engine=\"tensornetwork\"), uccsd.e_fci"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a58aea25-6fda-410d-ac9e-2bea15c60d4c",
   "metadata": {
    "tags": []
   },
   "source": [
    "## Built-in ansatz, UCC ansatz and Qiskit ansatz\n",
    "TenCirChem has implemented the $R_y$ ansatz.\n",
    "We can rebuild the `HEA` instance using the `HEA.ry` function, which builds the qubit Hamiltonian and the $R_y$ ansatz automatically,\n",
    "and then run the kernel. Note that the default initial guess is randomly generated. The result is exactly the same as what we've done step by step in the above code."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "89cc0336-9831-4edf-a4c5-eff91344d6ae",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-1.1202549357480136"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "hea = HEA.ry(uccsd.int1e, uccsd.int2e, uccsd.n_elec, uccsd.e_core, n_layers=1, engine=\"tensornetwork-noise\")\n",
    "hea.init_guess = init_guess\n",
    "hea.kernel()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f6c59f5a-0cbc-47b3-8993-899b376335e2",
   "metadata": {},
   "source": [
    "If desired, the UCC ansatz defined in the `UCC` class can be fed into the `HEA` class for noisy circuit simulation. Note that TenCirChem by default uses a multi-qubit controlled $R_y$ gate for UCC circuit simulation, which has to be decomposed into elementary gates first for actual execution on quantum computers. The gradient is turned off because parameter-shift rule is not applicable to the circuit. The energy obtained is significantly higher than HEA ansatz, because there are so many CNOT gates in the circuit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "3e71b357-13bd-4356-9b65-f65575a3037c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-0.8208126798813514, 18)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from functools import partial\n",
    "\n",
    "get_circuit = partial(uccsd.get_circuit, decompose_multicontrol=True)\n",
    "hea = HEA(uccsd.h_qubit_op, get_circuit, uccsd.init_guess, engine=\"tensornetwork-noise\")\n",
    "hea.grad = \"free\"\n",
    "hea.kernel(), get_circuit().gate_count([\"cnot\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8fec7993-23cd-4baa-8c4d-7b0e87ec8c68",
   "metadata": {},
   "source": [
    "`HEA` also directly accepts qiskit parameterized circuit as input. While qiskit offers a rich library of circuits, TenCirChem offers much faster simulation speed and friendly user interfaces."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "e9793a5e-551d-4701-aa0c-01266635f7dc",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAATEAAACuCAYAAABeIjpKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAXTUlEQVR4nO3de1xUdf7H8RfDRZCLgpogKKCCIArkBdcyk8JVUsvcKE2t1LSb2abJ1q9tq+1nRl72t2atsq3bbVPKzFXcUje1TLtgXjLxkUJichltBJWLCMj8/pgkCdQZnJnD98zn+Xjw0DnzPed8xvPl7TlnvuccN7PZbEYIIRRl0LoAIYS4GhJiQgilSYgJIZQmISaEUJqEmBBCaRJiQgilSYgJIZQmISaEUJqEmBBCaRJiQgilSYgJIZQmISaEUJqEmBBCaRJiQgilSYgJIZQmISaEUJqEmBBCaRJiQgilSYgJIZQmISaEUJqEmBBCaRJiQgilSYgJIZQmISaEUJqEmBBCaRJiQgilSYgJIZQmISaEUJqEmBBCaRJiQgilSYgJIZQmISaEUJqEmBBCaRJiQgilSYgJIZQmISaEUJqEmBBCaR5aFyCaMpuhvlbrKmxj8AQ3N62r0A/pA9aTEGuF6mth6xKtq7BN8ixw99K6Cv2QPmA9OZwUQihNQkwIoTQJMSGE0iTEhBBKkxATQihNQkwIoTQJMSGE0mScmI7sy9/GE8uSG03z9vIlrFM0Kf0mM/b6R3F3l02uZ67YB/T1aQQAyYkTSIq5BTNmysqNbP7mLZatn82PJw7y+B2ZWpcnnMCV+oCEmA5FhfYjpf+khtdjrnuYaS/H8NHXrzNl5Dza+3XSsDrhDK7UB+ScmAvw8fIlJvw3mM1mik/ma12O0ICe+4CEmIso+bnjBrQN0rgSoRW99gE5nNSh6toqTleaMJst50PWf7GMvKI9xHRNIqxTtNblCSdwpT6g+xAzmUy8/PLLrFmzhsLCQjp16sS4ceN48cUXmTVrFitWrOCVV15h5syZWpdqN29tepa3Nj3baNqQPuN49PZXNapIe3Xn4ehJqDwHnu4Q3A4CfbWuynFcqQ/oOsT27t1LamoqRqMRX19fevfuTXFxMUuWLCE/P5/S0lIAEhMTtS3UzkYNmsHQ+DTq6ms5UrKfrG0ZmE4X4uXp3dBm3jvjqTfX88zk9xqmnakqZfrCOGaMXsjN/SZqUbrdnT4LOw7BF3lQXv3LdDegdyjc0AtiQjQrz2FcqQ/o9pyYyWRizJgxGI1G5syZQ0lJCbt378ZoNJKRkcGGDRvIycnBzc2N+Ph4rcu1q9COUfSLTiEpJpW7ktN5Ycp6vi/M4a8fPNjQ5tFxr3GgYAdb9qxsmPbKh48QFzlEmc57JYWlsOg/sOm7xgEGYAYOFMGyLbBut+UmhHriSn1AtyE2a9YsCgsLmTlzJgsXLsTf37/hvfT0dBISEqirqyMiIoKAgAANK3W8uIjrSOk3mW37sjhQsBOwnNydk/YPlq6diel0MZ99u5pv87fx+3HLNK7WPn4qh79tgTPVV2675SB8vN/xNWlJz31AlyF28OBBsrKy6NixI/Pnz2+2Tf/+/QFISEhoNP3IkSPceuut+Pv7ExgYyD333MPJkycdXrOjTUx5BoPBnTc3/qlh2sCYkdwYfycZKyfxypqHmZ32OgG+HTSs0n6y91rOf1lr0344WeGwcloFvfYBXYbYypUrqa+vZ+LEifj5+TXbxsfHB2gcYuXl5SQnJ1NYWMjKlSvJzMxk+/btjB49mvr6eqfU7iihHXuSnDCePXmfsP+H7Q3TZ4xZSNHJPAbGpDIodpSGFdrP6SrYf8y2eczAzsMOKafV0Gsf0GWIbdmyBYDk5ORLtiksLAQah1hmZiZFRUWsXbuW0aNHk5aWxrvvvsuXX37JunXrHFu0E0y4+WkMbgbe3PTL/8Q+Xr6EBHUnMrivhpXZ1zcFUN+Cc1w5P9i9lFZHj31Al99OHj16FIDw8PBm36+rq2PHjh1A4xDLzs5myJAhdOvWrWHa4MGD6d69O+vXr2fs2LEtqmfAgAEYjUar23t5+JA50/bdgoQew9i84NK/veGdY9n48nmbl2uNqOgoaurOOmTZtkoY8zxRQ6bZPN+ZaujaLQJzfZ0DqrKNq/WB4OBgdu3a1aJ5dRlilZWVAJw92/w/aFZWFiaTCX9/fyIjIxum5+bmkpaW1qR9XFwcubm5La7HaDRSVFRkdXtvz7YtXpdWSoqLqa6t0roMAHpUlLd43qKiQsz1jvklt4X0AevpMsSCg4MpKytj9+7dDB48uNF7JSUlzJ07F4D4+HjcLnpQXllZGe3bt2+yvKCgIL7//vurqscWXh4+LV6XVkK6dGk1e2KGujMtmu/sGSNdQmzbVo7ian3A1t+Ri+kyxFJSUjh48CAZGRkMHz6c6GjLZRY5OTlMnjwZk8kEOG+Qq627yedrnPvMwUUPbbvqZRw+dLjVPHfy9Fl4/kPbz4vdOjiY5T+fK9Wa9AHr6fLEfnp6Oh06dODYsWPExcXRt29foqKiSEpKonv37tx0001A0+EVgYGBnDp1qsnySktLCQrS10WzetbOB+K72jaPmxsM7umYeoRj6TLEwsLC2L59O6NGjcLb25uCggKCgoJYvnw5GzZs4NChQ0DTEIuNjW323Fdubi6xsbFOqV3Yx5hrwc/7yu0uGNkXgpofjSNaOV2GGFgCKTs7m/LycsrLy/nqq6+YMWMGlZWVFBQUYDAY6NOnT6N5Ro8ezeeff94w/ALgq6++Ij8/nzFjxjj7I4ir0MEPHr7Jsld2JcPj4Ld9rtxOtE66DbFLOXDgAGazmaioKNq2bfwN0IwZMwgJCeG2224jOzub1atXM2HCBJKSkrjttts0qli0VJdAeOIWSI1vPsziu8IjN8OoRMvhpFCTy4XY/v2Wi+R+fSgJEBAQwJYtWwgJCWH8+PHcf//9XHfddWRnZ2MwuNw/lS74e8OIvvCnsfDYb6Gt1y/Tpw6FqNbxZaS4Crr8dvJyLhdiAD169CA7O9uZJQkncDdAZCfLvcQADLLnpRsSYjqSX7yPv6yeTtW5cjq3D+cPE97m6PED/M/rqYR16sVLMzYR6HcN1TVVLHp/GoeO5eDmZmBq6osMjb8DgMzsuWzbl0VUaD+ev2+tth9I2MTa7f+Pj/6HHfvX4OnRBnd3T6aMnMfAXiMA+OCzv7Bu56t4e/mxfPZebT+QlVwuxC5cV6lHC7Lu44k7/0nP0EQ+/noFmdlPMGLgFMI69WrUId//dCGe7m1488k8SkqPMGvJIBJ7JBPg24EZoxcQ3jmOnQfWavY5RMtYu/37Rt7ApJRnaOPpQ37xPmb/bSirninGx8uX3w19nJ6h1/Lav3+v2eewlZzo0Ym8oj34tPGjZ2giAMMH3MsXueuoratp0vbTfVmMHmy5OV5IUCTxPYbx+XcfOrNcYWe2bP+kmFTaeFq+6YgM7gtmM6crfnJmuXblcntielVSeoQjJft5YHFiw7RzNVWYzjS9ZvPEqR/pHPjLxfHBgRGcOPWjM8oUDmLL9r/Yxl3/JDioe6P+oBoJMR2J6TaIl6ZvbHh9x3P6eUCquDJbt//uw5/w9ubnyZi+udE1xKqRw0mdCAnq3mhvqrL6DNU1lXQMCG3S9pr23ThedrThtbGsgGvad2vSTqjDlu0PsC//Uxa+N4UXpqyn6zW9nFWmQ0iI6UTP0EQ8DJ58c2gzAOt3vsaNCXfh6dH0ityh8Wlkf2G5j3pJ6RG+zd/G9X3GOrNcYWe2bP9vf/iMjFWT+fN9/6ZHF/W/pZfDSR156u5/seC9KSxZ8xBdOvTkybvfocD4XZN2acPmsui9qdwzvwcGgzszb19KO9+OGlQs7Mna7b/o/WnU1p1jQdaUhmlPTnibyBC5s6vQWGRIX1577Mq3/fHx8uWPk7KcUJFwJmu3/5t/0NfDBORwUuc83L0orzrJA4sTKas4ccX2mdlzWbV1Pn4+gU6oTjiardv/g8/+wpI1Dyu1Z+5mNuvtsaHqc/YN8ewheRat5qaIl/PsGstNE9v5wPPjtK7m0qQPWE/2xIQQSpMQE0IoTU7st0IGT8uuuUoMnlpXoC/SB6wnIdYKubmpcX5JOI70AevJ4aQQQmkSYkIIpUmICSGUJiEmhFCahJgQQmkSYkIIpUmICSGUJiEmhFCahJgQQmkSYkIIpUmICSGUJiEmhFCahJgQQmkSYkIIpUmICSGUJiEmhFCahJgQQmlyZ9dWyGyG+lqtq7CNwdNyN1JhH9IHrCch1grV18rjulyd9AHryeGkEEJpEmJCCKVJiAkhlCYhJoRQmoSYEEJp8u2k0LUzZ+FYqeXnZAVU1Vimn62BL/OhaxAEtwN3+e9cWRJiQndqz8O+H+HzQ1Bgar5NzXlY9aXl7229YFAPuD4KOvo7r05hHxJiQjfMZth1BP69GyrOWT9fVQ1sPWj5uTYcfjcA/LwdV6ewLwkxHdmXv40nliU3mubt5UtYp2hS+k1m7PWP4u6uz01++iy89xUcKLq65ew5CoeNcEcSJHazT23O5Ip9QF+fRgCQnDiBpJhbMGOmrNzI5m/eYtn62fx44iCP35GpdXl2ZzwNf/vEEmT2UHEO3tgOI/rCyL5qXk7lSn1AQkyHokL7kdJ/UsPrMdc9zLSXY/jo69eZMnIe7f06aVidfZ04A0v/CxXV9l/2xv2WP1Pj7b9sR3OlPiDfybgAHy9fYsJ/g9lspvhkvtbl2M25Wsjc6pgAu2Djfst5NtXptQ+A7Im5jJKfO25A2yCNK7Gf9XvBVGHbPLNHQoCPZejF4o+tm2fNLogKhnY+NpfYquixD4CL7ImZTCbS09Pp2bMn3t7edO3alccee4zKykqmTZuGm5sbS5cu1bpMu6mureJ0pYlTFT9xpGQ/S9Y8Ql7RHmK6JhHWKVrr8uwi77hlCIWtAnygfVvLn9aqqoH3v7Z9XVpyhT5wge73xPbu3UtqaipGoxFfX1969+5NcXExS5YsIT8/n9LSUgASExO1LdSO3tr0LG9terbRtCF9xvHo7a9qVJH9XThf5SzfFUJhKYQpshPjCn3gAl3viZlMJsaMGYPRaGTOnDmUlJSwe/dujEYjGRkZbNiwgZycHNzc3IiPV/Ds7SWMGjSDjOmbmTftP9x/Swb+bYMwnS7Ey/OXwU/z3hnPC2/f2Wi+M1Wl3PXnED7Z/S9nl2yT46fh8HHnr3fHYeevs6X03gcupusQmzVrFoWFhcycOZOFCxfi7//LcOz09HQSEhKoq6sjIiKCgIAADSu1r9COUfSLTiEpJpW7ktN5Ycp6vi/M4a8fPNjQ5tFxr3GgYAdb9qxsmPbKh48QFzmEm/tN1KJsq+3M02a93xyBakXutqr3PnAx3YbYwYMHycrKomPHjsyfP7/ZNv379wcgISGhYdqF0EtKSqJNmza4qThI6FfiIq4jpd9ktu3L4kDBTsBycndO2j9YunYmptPFfPbtar7N38bvxy3TuNory9NgLwwslyr9eFKbdV8tvfWBi+k2xFauXEl9fT0TJ07Ez8+v2TY+PpazuxeHWF5eHh988AHBwcEMHDjQKbU6w8SUZzAY3Hlz458apg2MGcmN8XeSsXISr6x5mNlprxPg20HDKq+s9jyUnNJu/YWl2q37aumlD/yabkNsy5YtACQnJ1+yTWFhIdA4xIYOHUpJSQnr1q0jJSXFsUU6UWjHniQnjGdP3ifs/2F7w/QZYxZSdDKPgTGpDIodpWGF1ik5BfVm7dZ/TOEQ00sf+DXdfjt59OhRAMLDw5t9v66ujh07dgCNQ8xgsH+uDxgwAKPRaHV7Lw8fMmfa/yzyhJufZuvelby56U8sfHArYBkEGRLUncjgvle17KjoKGrq7HTdz2UE90pmyNS3m33vwhiwywnw/uXP526/dLtLjSP7+L+f8fT4u62stuVcrQ8EBweza9euFs2r2xCrrKwE4OzZ5v9Rs7KyMJlM+Pv7ExkZ6dBajEYjRUXWX5ns7dm2RetJ6DGMzQsuvZsS3jmWjS+fb9Gyr6SkuJjq2iqHLPtiXp1PX/K9C2PArGEwWN/2YrV1Zpu2ZUtJH7CebkMsODiYsrIydu/ezeDBgxu9V1JSwty5cwGIj493+Mn74OBgm9p7eag3NDykSxen7IkFtrv0Db/OWLH6AG9LgNXXw5nLXK50qWV5uJsJDQ298oqukqv1AVt/Ry6m2xBLSUnh4MGDZGRkMHz4cKKjLaOUc3JymDx5MiaT5W55zhjkautu8vka9Z45ePjQYac8c/D4aZif3fx71lxG9Nztlj2wM9Xw3Ie2r39s6jBWPV9o+4w2kj5gPd2GWHp6Ou+++y7Hjh0jLi6OmJgYqqurycvLIzU1lYiICDZu3NjofJirWvTQNq1LsFqnAGjjAefqtFl/V0VG7NtKpT7wa7r9djIsLIzt27czatQovL29KSgoICgoiOXLl7NhwwYOHbJceCchphaDG4QGarf+rmqNPnAJut0TA4iNjSU7u+mxR0VFBQUFBRgMBvr06aNBZeJq9AmDH35y/nrb+WgboKJ5ug6xSzlw4ABms5no6Gjatm36LdDq1asByM3NbfQ6IiKCAQMGOK9Q0ayk7vCffVBX79z1Do6SpyK1Ri4ZYvv3W26BcKlDybS0tGZf33vvvbzxxhsOrU1cmZ83JIY792aFBjcY3MN56xPWkxBrhtms4ZBwYZXUePj2GNQ46QT/Tb2hXcuGbgkHc8md4yuFmKryi/cxc0kSUxfE8tTfR3Kq4if25W9j1FM+PLA4kbKKEwB8/PUKpi/qy4g/eLBm+/81WkZm9lzunteNZ98Y6/wPYIMOfnDrtc5ZV3A7ywNDWjtrt/+Kj55m+qK+PLA4kQcWJ7J176qGZaiy/S/mkntiF66r1JsFWffxxJ3/pGdoIh9/vYLM7CcYMXAKYZ16sXz23oZ2UWH9+eOk91i1pendPWaMXkB45zh2HljrvMJb6LooyC2C3GLr57kwiNWagbEAnu5w92DwcLe9PmezdvvfOWwuU1PnAWA6XcS0BbH0i0qhnW9Hpbb/BS65J6ZHeUV78GnjR8/QRACGD7iXL3LXUVtX06Rtjy4JhHeOxc1N7c1vcIN7b4BIGx7cs/hjyyBXawbGuhtgyg3QTYFhFbZsfz+f9g1/P3uuAjNm6s1O/pbEjlxyT0yPSkqPcKRkPw8sTmyYdq6mCtMZx1/np6U2HvBgMqzYDt+X2He5U4dCrxD7LdORbN3+H36+hHU7X8V0qpDH014n0O8aJ1VqfxJiOhLTbRAvTd/Y8PqO5/TzbMHLaeMJDwyDzw7Bhr2We45djehgGD8Igpq/DV2rZcv2v33ILG4fMov84n28tHISA6J/q9x9xC5Q+3hCNAgJ6s6JUz82vK6sPkN1TSUdAxx/sXJrYDDAsBhIvwX6hrXsqd0d/OCuQfDQTeoFWEu3f48uCXQMCGVf/jYHV+g4siemEz1DE/EwePLNoc30jx7O+p2vcWPCXXh6aHBFroY6BcC0G6GsEr7IswzDOH4GLjVqxrcNdO9k+ZKgV4jlPJuKbNn+R4/nEt65NwDFpnzyivfQ7efXKpIQ05Gn7v4XC96bwpI1D9GlQ0+evPsdCozfNWm3MecN3tj4Ryqqyth5YC3vf7qQF6asp2eok8YsOEGgL9ySYPk5VwdFpZYH7dadt5yw9/GCsEBLOx08RgGwfvv/fUM6xtIjuBs8cXf3YObYpYR3jtWgYvuQENORyJC+vPbYlW/7M2LgfYwYeJ/jC2ol2nhA92ssP3pm7fb/36mXuJeRouScmM55uHtRXnWy0WDHy8nMnsuqrfPx85ErnfXAFba/m1musWl1VLwhXvIsNLkhnl5JH7Ce7IkJIZQmISaEUJocTrZCZjPU12pdhW0Mnvr5lq81kD5gPQkxIYTS5HBSCKE0CTEhhNIkxIQQSpMQE0IoTUJMCKE0CTEhhNIkxIQQSpMQE0IoTUJMCKE0CTEhhNIkxIQQSpMQE0IoTUJMCKE0CTEhhNIkxIQQSpMQE0IoTUJMCKE0CTEhhNIkxIQQSpMQE0IoTUJMCKE0CTEhhNIkxIQQSpMQE0IoTUJMCKG0/wcT2IJXR//e4gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 370.906x200.667 with 1 Axes>"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from qiskit import QuantumCircuit\n",
    "from qiskit.circuit.library import TwoLocal\n",
    "\n",
    "ansatz = QuantumCircuit(2)\n",
    "ansatz = ansatz.compose(TwoLocal(2, rotation_blocks=\"ry\", entanglement_blocks=\"cx\", reps=1).decompose())\n",
    "ansatz.draw(output=\"mpl\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "efd1b20d-101a-484e-be43-6f1be5e942e5",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-1.1202549357480136"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "hea = HEA(h_qubit_op, ansatz, init_guess, engine=\"tensornetwork-noise\")\n",
    "hea.kernel()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ef4cd7bf-2405-48d7-885a-47c698cf6d71",
   "metadata": {},
   "source": [
    "## Effect of gate noise on VQE energy\n",
    "Here we show how the CNOT depolarizing error affects the optimized VQE energy.\n",
    "If a one-layer ansatz is adopted, the energy increases linearly with the error probability up to $p=0.8$. This is a consequence of the special form of the ansatz: there's only one CNOT gate in the circuit.\n",
    "Increasing the number of layers does not lead to more accurate energy estimation. \n",
    "Rather, the noise caused by more CNOT gates deteriorates the result."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "b39aa3a3-5d11-4685-bba2-884ae5a59c9f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "14.445068836212158"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from tencirchem import get_noise_conf, set_backend\n",
    "\n",
    "# cnot error probabilities\n",
    "error_probs = [0.001, 0.05, 0.10, 0.20, 0.40, 0.60, 0.80]\n",
    "# number of layers in the ansatz\n",
    "n_layers = [1, 2, 3]\n",
    "e_list = []\n",
    "# record run time\n",
    "import time\n",
    "\n",
    "time1 = time.time()\n",
    "for n_layer in n_layers:\n",
    "    hea = HEA.ry(uccsd.int1e, uccsd.int2e, uccsd.n_elec, uccsd.e_core, n_layers=n_layer, engine=\"tensornetwork-noise\")\n",
    "    # set HF init state\n",
    "    hea.init_guess = np.zeros_like(hea.init_guess)\n",
    "    hea.init_guess[1] = np.pi\n",
    "    for error_prob in error_probs:\n",
    "        hea.engine_conf = get_noise_conf(error_prob)\n",
    "        e = hea.kernel()\n",
    "        e_list.append(e)\n",
    "e_array = np.array(e_list).reshape(3, -1)\n",
    "time.time() - time1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "51e17cb8-3690-4f8a-87f4-0c58dcfe8a6c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhgAAAEbCAYAAACY+P0jAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAABrOklEQVR4nO3deVzb9f3A8Vc4WqAHIb2LrSVUe3g2gPWqOhu8rdOGss7pvIC6Td2hRHTuN+cmhqnbdB5QN53HKhCvanVK6lXrVUjrUbVa0tpae4f0gpYj398fX5MSSIBALuD9fDz6aPL9fr7f7zsB8n3nc2oURVEQQgghhAihuGgHIIQQQoiBRxIMIYQQQoScJBhCCCGECLmEaAcQi4455hgyMzODPm7Lli2kp6f36dqhOEesnSeWYgnVeSSW8J4nlmIJ1XmiHUt9fT1r167t8/WF6DFFdHLxxRdH9LhQnyPWzhNLsYTqPBJLeM8TS7GE6jzRjiVU1xeip6SJJMYsXLhwQJ4nFGLpNcn7Ev7zhEIsvaZYel+EiARJMGJMLH0ghvI8oRBLr0nel/CfJxRi6TXF0vsiRCRIghFC8gHin7wvgcl745+8L4HJeyP6C0kwQkj+8P2T9yUweW/8k/clgPVVLDz4W1hfHe1IhOiWJBhCCNEfNO6At4ugaRu8U6g+FyKGSYIhhBCxTlHgnUXQsk993rwP3r0+ujEJ0Q1JMIQQItatr4INL4DSpj5X2sDxvLpdiBglCYYQQsSw73Z9Sttb19JxVUoFDa1vF/CdUybPErFJZvIUQohoUhRo2glDRkJCkrrN8QJ8/hCtDV8z8cBmv98ENSjQvJdVr8zFnfcpk5PHRjRsIbojCYYQQkTCge9hRy3sdcDeDbBvg/r/3g3QegAueQvSz1LLHtwF3y3v9gM6Abh0/3a+2FkLky8I8wsQIjiSYAghRF+1tcD+TYeTB08CMcsMY2apZTa+Au8UBTiBBhq3Hn6afjac/R/WxcdTv/Imzmnc7ffDuhUNLw4bg157VKhfkRB9JgmGH1u2bGHevHksXLhQxuMLIdRmjMZtauKQmgkp49TtG5bCihvgwHeguDsfN/n8wwlG2nQYkwUjM2BEBozUq49HZsCIIyF+6OHjUjMhNZMDe7/hqrHHsm7TSlLdrT5NJW5gb1w8vxgzg/91EfqSJUtYsmQJW7Zs6eObIERwJMHwIz09naVLl0Y7DCEGp/VVsOImmPMATM2L/PX3OGDDi52bMtoOqvvP/g9Mv1J9HJ+k1lx4HnuThx8SiLE5h8878QzIqw06nJ0JQ1g0ZgaV2z/z2R4HLBozg50JQ7o83vNFad68eUFfW4i+kARDCBE7PJNJNbvUyaQmngkpIeq82HYI9n3r24ThSSIMt0KmSS235xt4/3edj9fEwfBJh4eKAoybDZeuVBOKlPGg0fQpRLfiZtWedTy34z1uPtLk3V41fBwL9m/jkgM7SQBa0PDSsDFUjxjfp+sJEU6SYAghYkOgyaTOe66Hx7vVjpSeBGL0iTD6eHXfd2/CUiN0Guz5A+fawwlG2nSYmu9bEzEyA4ZPhvhE3+OGpsKEU4N9pT7alDZWNHzO8zve4/ntK9lyaBcAR6ccgWHkVLWQRsP1Y2ZydtNKtO5W9v3QNCJELJMEQwgRGzyTSXm0n0xq6gI1AVHaIO6Hj629G2B1WbvaiI3gbj58/El/OpxgDEsHFEhIOZwwtO8HMfrEw8eNOBLOeTbMLxY2H9zBXfXP8OKO99nZsufw5eNTuGjMbKYNO8Kn/M6EIRSNmcE/dq3jxjHTu20aESLaJMEQQkRf4w619gINnWoZbD+DVXeqHSlPvAWyf69ubzsEax/1LauJhxGT1eSh/Q06dSpctR2Sx/S5GaO3mtoOsfWQE33KBACS44byry2v48ZNWsIILhl7CvPHnY5RZyApXk0eNjXtICkukYPuFgCqR4zv1CySFJfI6MTUyL4YIXpAEgwhRHT5NI34acJwt0DDF+rjvY7D20dMgazb1WQiVa/+P/yIwzUc7cXFh64vRxD2th7g1Z0f89yO93h158dkjTyKd0+6H4DRQ1L569HXcfwIPWemHU+in7gnJ49l3WmPs6tdDUdHoxNTZZItEZMkwRBCRI/ihs0236aRQM5/CY4wHn6ekASz/xy+2Hppd/NeXt75Ic9tX8Ebu+00Ky3efZsP7qSp7RDJPwxJ/e0UU6DTeE1OHisJhOiXJMEQQkTeIRd89QR8/hC42yDjx7DxZd8RGh6aeMi4BDL6xzDLa9bex9KdH3ifH51yBPPHnc5lY08na+RRaKLURCNEpPW7BKOsrAytVguAy+WiuLg4qONzc3OpqakJQ2RCiG7t/lxNKtY9pU6PDTAkFU74HWx5G5r34NtMolHX6DjjkSgE27VNTTt+GPnxHv859hYyfuhbcenY09jYtP2HpOI0jhk+RZIKMSj1qwTDk1wUFhYCYLPZKCoqory8vEfHW61WbDZbOEMUQviz/SP4sAS2vHV4m+4YOO4GOPpySBwOZz4KNT/pcKCibo9C/wl/vjmwhed2rOC57e9Ru/dr7/YXdqz0Nnf8fGIuV6WfE60QhYgZ/Wq59tLSUhYsWOB9bjQaqaio6NGxLpcLp9MZrtCEEF1xt6rJhSYe9PPVhb3yP4NjitTkAtShqBmXqmXgh7KXqduj7Mv9mzju/UKOXnk1Jd/8m9q9X6NBwxlpx/GPadeTP/4sb1mprYgch8OB2WwmLS2NzMzMTvs9+3Jzc7Hb7VGIMDTsdjtZWVl+X2Ms6zcJhsPhwOVyeZtH2utJrURVVZVPciKECJMddbD8Knj/lsPbxp8Kp/0NfrYBzrOqq4Z2vBFrNGptReII9XmUmkYUReHjPV+xfPdq77ZJSWNY3/g9CZp4zh2VTfnMm9h65rO8k3MfNx55KelJoyMeZyxoc7tZ8eV2qj/YyIovt9Pm9rMeSxjp9XosFgslJSU4nU7MZrPPfovFQmFhITU1NRgMhh6d02w2k5cX/inqHQ4HFRUV5ObmdlvWYDBgsVjCHlOo9ZsmEofD4Xe7VqvF5XJ1eazNZsNoNHZZpj3PYmcesuiZEN1oa4b6avjsn7D9Q3Vb4nDI+T/1f40GTvh19+dJGQtnlR9eiyRCTSNtShsrG9by/I6VPL/jPTYf3MkJI/SsOUWdZ2N4QjLLZt3FrJFTSfMkQDHOs8iZR6gXO3tp1WbMz9Sxxdno3ZauS8FyeRaX5EwK6bW6o9Vqqa6uJjc3l6KiIvR6vXdfsN/6c3Nzu72n9JXdbqe2tjaomnWdThfWmMKh3yQYgeh0um5/QC6XC71e3+NfGlnsTIge2r8FviiHtRXQtF3dFpeoTrV93K8gYVjw55y6IGLNIm87P2HJtrd4ccf77Gh2ebcPj09m+rBJHHI3MzROnfTq7FGzIhJTqHT8YhTKxc5eWrWZKx5c0WnWku+djVzx4AqeumFOxJMMo9GI0WgkLy+Purq6Pp0n3AwGAwaDAavVGvZrRVPUEgyr1UplZWW35UpKSrqs2uouuaioqPB2ChVChNhn/4TV96iPh02EY66HmQWHlzOPMQfbmhkal+jtJ7H4u1f57za142lawgjmjT2Z+WNPJ3dUlnc2zcHgwKHWgPviNRqShsR7n+9tauaWp2v9ruqioM7Fan6mjouy0jnYErjJJE4DyUNCewuqrq4mLS2t28/9srIyby2Hw+Hwjka02+0UFBTgcrmor6/3lq+oqPB+SXU4HN5jzWYzTqeT6upqjEYjDoeDrKwssrOzKS8v96lJCQeXy+WNraamhqKiIu/90mq19ig+z3vheV0mk8l7rMFgoKioyDvyMthmmqglGCaTCZOp+0lmPAL9oDy1E/7Y7Xays7N7FZ8QooOWRvjmv5A2Ayacpm479nrY/gEc+0t1LouOi4HFgH2tjby2axXPbX+PZbs+4v2T/s7xI9TPjJ9NmMvwhGTmjz2dH+lO9Dub5mAwvqAq4L5zTpjIc787y/tc/8vnOdQaOHFQgC3ORt5ft5Mr/vkeu/cd8lvOkKHjnTvP623Ifmm1WiwWC2azmQULFvjts5eXl0dRUZG3psLhcHinL/D0dSgqKvKW99QytC9vs9koLCxEq9X6nEuv11NSUhL09Am9VVpa6m0SMplMZGZmUldXh1arxWQydRtfXl4e+fn53ntxbm6u91ye8+t0OvLz83tUIdBRv+nkqdfr0Wq1fvtiBKrScjqd2Gw2ysrKKCsr83YAKisrG/BVU0KEzB6H2mHzySPg7QKw33N434jJ8OO3YWpeTCUXDS37+M+WN7hk9f8x5u088j/9C1Xb3+FA20Fe27XKW+78MSdRPvPXnDM6e9AmF8FSlAAr0nawzdUU5kj8Ky4uRqfTUVBQ0Gmf3W7v1CdPr9d77xXgv69DdXW1t4ldr9d7v7gajUacTqfPCJVw11q050l22l+7/fOu4nM4HFitVp8v+rm5ud5pHzz3W09zTm86mfarv6iSkhJv5ghqZtm+GszzhnmyM0+bnIfdbqeioiJi2aUQ/ZZnCu/PHoRvl+Gd/GpkBhxxtrp+SIwOx1y1Zx2nfvxrWtvNCjo1ZSLzx87hsnGnkTNyWhSji03bFgfu8xLf4ee85DdnMv/et7s953htMmvvvyTg/rgw/vpUV1eTlZXVaWhqbW2t3wTA08Tg78uqyWSivLyctLQ0DAYD+fn5PveQwsJCKisrMRgMQQ8o6Kvq6moAb9ON0+ns1G0gUHw2mw2tVuuTkOzevdvnS3xfk6V+lWAUFxdTVlbmnfuivr7eZ5Itm81GeXm53wSifZ+PoqIi8vLyIvqLIES/8r/5sOHFw88nnaNOijX5fHXhsBDb1LSjVwt6bT64g+e3ryQpLpGiSRcBcMIIPcPjkzkiaTTzx57O/HFzOFZm0+zSsKE9vxXMPXY86boUvnc2+u2HoQEm6lI4ddoY4uOiU0luMBgoLCwkLy/PZ+hqb0eH1NTUYLfbqayspLS0FMB7n8nPz2fu3LlYLBYcDkef7ytZWVk+N/mumlzsdjulpaXk5uZiNBr9JgSB4vN0L2gfb6jvif0qwQC6rH0oLCwM2LEn2D4fQgwqzi/V5o7EH0Z9TD4PvlsO069S+1ekhe9b/6amHUxbebV3SXJ/kuISWXfa40xOHsv6xi08t12dovvjvesAyEgeT+ERF6LRaBgSl8jXp/+bMUO0YYt5MIuPi8NyeRZXPLgCDZ0mdgfAcnlW1JILD4vFQkZGhncuDFBvoJ4EoT2Hw0F+fr7f83g6jHqaCkpKSpg7d673XmQwGNDpdFit1pAMJe3pCBiXy8XcuXOpq6vzJhaeBKp9R9RA8RkMBr/vRaD5pnqj3/TBEEKEmLsNNrwES43w7Ez4+unD+6ZdCT/fos5FEcbkAmBXy54ukwuAg+4W/rn5JU54fxFHvXc1t37zLz7euw4NGuZoj+XGyT/2aRKR5CK8LsmZxFM3zGGiLsVn+0RdSlSGqLYf8eGh1WpZvHixT5OBwWAgOzvbp1nA04wS6AuoZ6RGex1rCoqKiigoKAj6S2xfZpf2TD7ZPhbP+To2DfmLz2g0kp2d3ak/YiiX0+h3NRhCiD46uBu+/Bd8/jDs+1bdpomDPesPl0lIjk5sXXA0buPT/Q7iNXGcrTuR+WNP55KxpzJ+aP+bgGgguCRnEhdlpfP+up1sczUxXpsc8WYRz1ThVqsVl8uFxWLx+fZtMpk6jX6oqanBbDZ7myHq6+u9tQZ2u927z2w2+5zPcyN2OBwsXrzY55wmk4lVq1bRUw6Hg/Lycmw2m/eamZmZAWvg/cVlMBgoLi7GbDaTk5MDqH0yzGZzp9qYQPF53gun0+mt3TCZTNhsNm+Titls7jR5WU9plJ52CR5E5s2bJxNtiYFHccPbRWpNRdtBddtQnTpvxbHXw4gjoxKWfe83ZH34y27LPX2cmVZ3GxePPRld4sgIRBZb2tzuPt3M5XMtfFwuF7W1tTHbry9a8UkNhhADmeJWaydA/b9xm5pcjJ6ldto86icxWVvhz4xhkzGMPCraYURFLE3LLVTt+ypUVVXF3ISOsRCf9MEQYiA6sBU+/iM8OflwMwjASX+CS1dCXh3MuDomkosdh1zRDiGmeablbp9cwOFpuV9atTlKkQ1uZrPZ2zcjFtcJiYX4JMHww7PYWfuFgoSIeYoCW1fCGwvhqclQeycc2AJfPXG4zJhZMOHUmJnDwq24uf7LB6IdRsxqc7sxP1MXcFpuUKfl7moV0yVLljBv3ryQL3Y22Hlm+6yoqIjJEYqxEJ/0wfBD2ipFv9J2CL5+Rl0XZNfhJcYZf5q64Jj+MoixdTUURfHOS3GP41lK1v+722PqTn5o0DWRrPhyOxeULu+23Kslc5kzo+v1X+RzTUSa1GAI0d+1NcN7v1aTi/gkmH4N5NnhsvfUPhYxlFx827Sd+Wv+xBPfv+HdtnD8j0iK63qa8aS4REYnpoY7vJjT0+m2ozUttxBdkU6eQvQnigJb3oSNr8Bp96tNHUNGgOFWiEuAGddC0qhoR9nJwbZm7v22mrsdz9LkPsQHe77g8glnMyQukSNTxrHutMd7NZPnQDde27M+Mj0tJ0QkSYIhRH/QvA++fkptBmn4Ut2WOR8mnK4+zroterF145WdH3LTV4/gaNoKwJlpx/Pg9F8ypF2txeTksYMygQjkf2u2oB02hFOnjenxtNxCxBpJMISIlvVVsOImdbbMqXn+yzSsg88fUjtqtuxTtyUOh2k/h2HpEQu1NxyNW7nxq4dZtusjACYOHcV9RxeSP/4sWRckgB17mih+uo7nPtrE1PEj+ODPF/SLabmF8EcSDCGioXGHOulVswveKYSJZ0JKh2/w2z6A5089/Fw7Te20Oe1KGBL7E03taHaxbNdHJGoS+M2Rl/F7/U8ZkZDS/YGDkKIoPPWug98/u5qGA83EaTScPysdBcU7LXfHeTAmyjwYIsZJgiFEpCkKvLPocI1E8z5493o4a7Ha/DHhNHX72JMg9ShIm6EmFkfMPTxpVgxSFIWvDmxmxvDJAJysncHfpi3ivNHZTB82OcrRxa712/Zy0+OrePfL7QCccGQaD14zm1kZh+cuiIVpuYUIlvx2ChFp66tgwwvgWZxLaQPH8/DEePjfZeqwU1CXRc//BC54CSblxnRy8eX+TZxTdyuzPrweR+NW7/ZfH3mZJBdd+GrLHk65/TXe/XI7yUPiuSv/RN7+47k+yYVHfFwcc2aMI++UKcyZMU6SCw6vRZKWlkZmZman/Z59ubm5nRYA60/sdjtZWVl+X2MskxoMISKpcYdae9GpRR1wt0DSaNi/GVKnqttiYKbNruxrbeQuxzP87dvnaVXaGBqXyEd7vkKfMiHaofUL0yaO5PTpY2lzK/z9qhz040ZEO6Qe2dS0IyZG/ej1eiwWC6NGjaK0tNS7EJiH53H7bd3xLCpWXV0d8njbKysrAw6vBFteXh6wrMFgwGKxeCfP6i8kwRAiUnyaRvyMCdDEgXb64eQihimKwpJtb3Hz1xVsPaQuEX3xmJP527RFZKZMjHJ0sWv/wRbuf+ULbjh/BmnDhqDRaHjyV6czPCmh33R83dS0g2krr+aguyVgmaS4RNad9njERgZptVqqq6vJzc3ttPJnsN/6c3NzcblcIY7QV8dEqKioiNzcXGpqagIeE4vTkXdHEgwhIsW5Vm0aCURxw4bn1XK6YyIXV5AUReEC++38b3ctAJnJE/nH9Ou5cMzsKEcW29745Ht+859VbNp1gB17DvLPa9X3a0Ry15OMxZpdLXu6TC4ADrpb2NWyJ6JDj41GI0ajkby8PO/y6709Tzi5XC7sdrvPYmRFRUVkZWXhcDh6tSx6rJJGPCEiRXcMZFwKmnj/+zXx6rTeMZxcAGg0Gk7RziQ5bih/mXo1n59aIclFF3buPcg1D69k/n1vs2nXASaPHsa87Nga+XGgtSngv4NtzT5lmzx9hPp43mDO01PV1dXY7XbvIl+BlJWVYbVasVqt3qYKCNzXoaKiApvN5i3vOTYzM5O0tDRsNhug9gnx9PlwOBwBr19bW+uz35NUBFtz4nK5vPEUFRX59DPpaXye4z3/tz82Ly8Pm82G2WzGbDYHFRtIDYZfnsXOFi5cyMKFC6MdjhgIDnyv1lCc+ShseQua99BpVoMhI+GMR6IVYUBuxc2T39uYNuwITtHOBKB4ygKumniOTI7VBUVReOa9Ddz2X7t36On15xzN7+cfz/CkyNVaLFmyhCVLlnS52NnwNy8JuO+C0SexzPBn73Njbc9vNFNWXBmwr0b2yKNZdfI/e3yuntBqtVgsFsxmMwsWLPDWELSXl5dHUVGRt6bC4XB4myf89XXw3HTbl7fZbBQWFqLVan3OpdfrKSkpobi4uMsYGxoafLZ5EoBgay9KS0u9TUImk4nMzEzq6urQarWYTKZu48vLyyM/P9+7GFpubq73XJ7z63Q68vPzqaysDCo2kBoMv9LT01m6dKkkFyI0nF/Cc6fAK+ep64Kc+Sid+2Ao6vaOc2FEWd3erznt499w9dp7+cWXD9L2w8iXpPghklx04x+vfsn1iz+k4UAzx03W8ub/ncM9l2dFNLkAWLhwIUuXLiU9PbYnZguV4uJidDodBQUFnfbZ7XZsNptPM4her8fpdHpv8v76OlRXV3trF/R6PdnZ2YCadDidTp+ag940cZSWllJeXu43IeqKJ9lpf+32z7uKz+FwYLVafVZazc3N9XY21Wq1OBwODAaDN/EKltRgCBFOW1fCqxfDoQbQHg3Ne2HqAlhfCRuXqkNUNfGQcYm6PUbsbt7L79c/Qfl3y1BQGB6fzOUTzsatKMT3j76IUXfFGZlU2L7murlHc8N500lMiN3vc/vPfingvvgOTXq2bAunr/ptj867cc6TAffFhXHYdXV1NVlZWZ2GptbW1vpNAPR6PTU1NX77X5hMJsrLy0lLS8NgMJCfn+9TQ1FYWEhlZSUGg6FT8tITZrOZoqIiCgsLgzoO8I50cblcOBwOnE4nTqfTp0yg+Gw2G1qt1ich2b17t9+mm96K3d94Ifo7xwuw1KgmF+NOhktXwojJ6gJlZz4KiT8MSYyhppE2pY3yza9w9MqrefS7V1BQ+On4H7HutH9z85Q8EuPkO0kgqzc4uX3JahRFrZ0aNWIoq8su5rcXzYzp5AJgWEJywH9JHVbjTY4fGpLzBnOeYBkMBgoLC8nL852Cv7ejQ2pqaqirq8NoNFJaWurTbyM/P9/b58PhcARVC+Hp69A+ucjKyiItLc37r/21OrLb7eTl5VFVVYVWq/WbEASKz+VyodfrvZ1jjUYjFouly5EswYrt33oh+qvPH4b/zYe2gzBlHsxbDsmjD+9PGQtnlUPyeDizPGaaRl7e+SGLvnwAZ8s+jhuewTvZ9/LM8SVMjMEVWmPFgUOtlPzXzll/fJ0HXvuSF1dt9u4bmhigQ68IO4vFgtPp9KnaNxqNfjtfOhwOcnJy/J7Hc3P2NBNs2LDBpz+CwWBAp9NhtVqDGkrqqTnwJBeeWoi6ujoaGhq8/wL153C5XMydOxeLxUJhYSF6vd6bQLV/jYHiMxgMft+LUA7RlQRDiFD77CF495eAAjML4bznINHPGhxTF8DVWwMvdBYhbsXtfXzJmFO5cPRsHpj+C+wnP8wZuuOjGFnsW/7ZVmaXLOOf//sKt6KQd/KRnD49NpLFcBmdmEpSXNf9SJLiEhmdmBqhiA5PVtWeVqtl8eLFPk0GBoOB7Oxsn2YBTzNK+74I7blcrk6jUjrWFBQVFVFQUBDwHB3Z7Xbsdrv3Jm+3270dKnvK4XB4ayE8PK+1Y9OQv/iMRiPZ2dneTqwe7d+bvpL6TiFCLePHsLoMZl4HWb9Xm0RiUKu7jYc3L+XR75bx4ex/MDJhGBqNhpdn/anfTPoULTv3HuS2/9p59v2NAEwalcLfrsrh3BMGfkfKycljWXfa4zExk6dnqnCr1YrL5cJisfg0UZhMpk6jH2pqaryzdYKanHjmzbDb7d59nsmwPOfz3IgdDgeLFy/2OafJZGLVqlU9itlT8+ByuToN/QzUkdJfXAaDgeLiYsxms7f2pbq6GrPZTH5+fo/i87wXTqfTm9yYTCZsNhsWi8V7vY6Tl/WURvE0GAqvefPmsXTp0miHIfoTdyu075/QvDemVzx9x/kpv/rqn3y+fyMA9x5dyO+m9Ozb12CnKApz//QGq+p3o9HAL86ZFvGhp70hn2vh43K5qK2tDfskXb0VrfikiUSIvmrcBs+dDN88e3hbjCYXWw7u4qeflnJW7c18vn8jusQRlM+8iV8feWm0Q+s3NBoNfzCdwLGTtLz5h+gMPRXR176vQlVVVcwlF7EQnyQYQvSF62t1jouddfD+LdDaFO2I/FIUhb9uqGL6ymtZsu0tNGhYdMRFfH3a4xQecWGnoYjisNY2N39f9gX/fvMb77azjhnPyrvOJztzdBdHioHMbDZ7+2bE4johsRCfJBhC9Na2D+H5U2HfRhiZCT9+K2ZXP9VoNHyy38H+tiZOTp1B7cn/5JGZNzIqRmtaYsWajU7O+uPr3FG5htuWrGZrQ6N3X1yc9FMZzDyzfVZUVPS4c2ckxUJ80slTiN7Y+DK8ka/WWIzJhguXxcxQU49vm7aToIknPUn9ll12VAFGnYErJxrDOsnRQNB4qJW/PP8ZD73+FW1uhbRhQ/jzT2YxXhubCaSIPM8Ml7EqFuKTBEOIYK2tgHevV9cWmXw+nFsFicOjHZXXwbZm/rqxitINlVw45iSqT7gDgIlJo7gq/ZwoRxf73vx8Kzc9/jEbdx4AYP7syZT9LIuxqZJcCBEMSTCECNbeejW5mH6NOiNnfOx08Htl54fc9NUjOJq2ArCzeQ9NbYfCOmviQPLd7gPMv+9tWtsU0nUp/O3nOZw/a+APPRUiHIIaprpmzRpqa2txuVzs3r0bgFGjRqHX69Hr9Zx44onhijOisrKySE9Pl9VUhX+KGxzPg35+zMxxUd/4PTd99QjLdn0EwMSho7jv6ELyx58lc1oE6c7qT9h/sIU/mE5gRHLsJI+91X41Vc98D0JEQrcJxsaNG7nnnnuoqqpi1KhRzJo1C51O5zOfudPpxOFwsGfPHkwmE0VFRUyZMiUC4YeHjBcXPloOgL1UnTQrISna0XRi223nQvsdNCstJGoS+M2Rl3GH/nKGx2iH01jy7c793PxULf9nOoFjJ6cB6oibgZiUyeeaiLQum0gWLVrEhg0bKCws5NFHH+3RCZcvX05hYSFZWVmUlpaGJEghoqZxB7x6EexYBfs2gTHw6pDRckrqDMYN1TJ92CQemP4Lpg+bHO2QYl5rm5tH3ljHn5/7lMbmNvY2tfD67bkAAzK5ECIa/CYYe/bsoaCggNtuuy3oZo+5c+cyd+5cVq9eTX5+PosXL2bkSBkKJ/qhPfXwynmwZz0kjYJjr492RAB8dWATj25exv3TiojTxDEsIZmPZj/A+CE6uTn2wKffNvCrf33E6o3qug2nTx/LA1efFOWohBh4/CYYVVVVVFVV9enEs2bNorKyksWLF1NQUNCncwkRcdtXwasXQtNOGDEFLvofpE2Lakj7Whu5y/EMf/v2eVqVNo4ZfiQFR1wAwIShstppdxoPtVL64mc8+Jo69FSbkshdP5nFlWdkypwWQoSB38HwoUwIJLkQ/c63r8FLZ6nJxehZMP+DqCYXiqKwZOtbTF95LX/dWE2r0sZFo2fzI90JUYupP6p8fyN/X/YlbW6FS0+aTO09F3HVWVMlueir9VXw+ARYXx3xS3sW40pLSyMzM7PTfs++3NzcTiuM9id2u52srCy/rzGWyTBVIdpr2Q9v/hxaG2HSOXCuFYaMiFo4n+3bwK+++ifvNnwGQGbyRP4+fREXjTk5ajH1J+07bF55ph7bZ1u5fE4GF8w6IsqRDRCNO+DtImh2wTuFMPHMiE44p9frsVgsjBo1itLSUu9Kox6ex4FWKfXHs2ppdXX4EiaXy+VtJaivr/eu0Np+Jdj2DAYDFovFOztnfxH0dH6PPfYYU6dOJT4+nsceewxQO3aWlJSEPDghIi5xOJz3PMy4Di54OarJhaIoLPryH7zb8BnJcUP589Sr+PzUCkkuekBRFKre30jun2toam4FID4ujmdunCPJRagoCryzCFr2qc+b96kT0EWBVqulurqasrIy7zLsHsF+68/Nze203Hmomc1msrOzKSwsxGKxoNPpyMvL6/KYWFzvpDtB1WAsXryY8vJy75vjGVM9d+5csrOzeeyxx7juuuvCEqhHWVmZzxDZ4uLiHh1nNpu9v2g6nS4m544XUeJuBdc60B2jPp9wuvovDDY17WBXy56A+3UJI5iYNIohcYloNBr+Pu16yjZWcd/RRUxOjq2pyGPVpl0H+M0TH/PGp+pkYxW2b7jpghlRjmoAWl8FG144/FxpU+eHWV8FUxdEPByj0YjRaCQvL69P831EYtVRh8OBzWbzTuWdmZnZ536PsSioBMOzprxH+x9iamoqDQ0NoYvMD09yUVhYCIDNZqOoqIjy8vKAx7hcLubOncvy5cvRarXetqwg5hcTA1lLI9T8BL5/B378LowOX7+GTU07mLbyag66WwKW0aDhl5Mu5sEZvwIgJ3Wad6pv0bU2t5tH3/iau577lAOHWhmSEEfxvGO4/pyjox1a7Gs5EHifJt53/peWA9C0A94pAjRA+89Sjbp94llqU0mX540L+eKA1dXVpKWlUVFR4b1P+FNWVoZerwfUm73ni6rdbqegoACXy0V9fb23fEVFBXq9HpfLhcPh8B5rNptxOp1UV1djNBpxOBxkZWWRnZ1NeXm5t1xHNTU1Ps9XrVrVq8TG5XJ5Y6upqaGoqMibtFit1h7F53kvPK/LZDJ5jzUYDBQVFXnjDaapCYJMMNLS0rrcH+4hcqWlpWzYsMH73Gg0kpub22WCUVpaSn5+vrfWw2AwdPrhikGqaRe8ejFs/xDih8L+TWFNMHa17OkyuQBQUHhq63JKj7pWJsoKwmebGrjh3x9R51CHnp46bQwPXH0S0yamRjmyfmJxF2vpTL4ALlp2+Pm/x0BbU4DCyuGmkvOeg6emwMFd/ouOyYa8Vb2N2C+tVovFYsFsNrNgwQK/fRry8vIoKiry3tAdDge5ubnU1NT47etgtVoBfMrbbDYKCwvRarU+59Lr9ZSUlPS4Zt1zfpfL1as+H6WlpRQVFXkTg8zMTOrq6tBqtZhMpm7jy8vLIz8/31ujn5ub6z2X5/w6nY78/HwqKyuDji+oPhjr16/3ed6+FmDv3r2d9oeSw+HA5XL5/YWx2WwBjysrK/Nmbp5ykagCEzFu7wZ44TQ1uRiaBvNsMOXiaEcFQOVxt0tyEaS/PP8ZdQ4nI5MT+cdVObxWYpTkImy6qf31NJU410YmnA6Ki4vR6XR+RzDa7XZsNpvPPUCv1+N0Or33B399Haqrq3G5XN7y2dnZgHovcTqdPiNUAtVa+FNRUYHD4ej1qqft72uea7d/3lV8DocDq9Xq012g/Rd2rVbrjc2TeAUrqBqM/Px8jjrqKG699VaysrLYs2ePd30Si8US1l63HTvueGi1Wu8PPtAxnqofvV5PUVEReXl5kmQMZjtXw7ILoHEbDJ+sznGhi502+jFD5cbYE21uN/Fx6nekv/4si+FJCdyVfyIT0lKiHFk/VLA/8D5NvO/zq3eC7XL4dpmaTPgrn3GJ2qfpio1dnDfoMQY9Vl1dTVZWVqehqbW1tX4TAE8Tg7/7gslkory8nLS0NAwGA/n5+T41FIWFhVRWVmIwGDolL93xNONUVFSQkZHBhg0bAo4k8cdzz/U03TidTpxOZ6dr+IvPZrOh1Wp9EpLdu3f73GuDSZb8CSrBmDVrFo888giLFi3yCSItLY2qqqqoLHam0+k6vaEenhi1Wq03Q7RYLGRkZHTZX2TLli3MmzfP+1wWPRtAdn0CL56hDkcddTxc9BoMmxjtqEQQnPsPcfuS1QA8UqCOqJk0ehiPLTo1mmHFHM8iZx5btmwJXDhxWM9PPGQ4nLUY/jsNmvfQqQ/GkJFwxiPBnzeEDAYDhYWF5OXlYTabvdsDfRntTk1NDXa7ncrKSu8SGJ4kIz8/n7lz52KxWHA4HN0mGC6Xi9LSUkpKSrzJhNFoxOVyYbPZKC0t9bm/dtXkYrfbKS0tJTc3F6PR6DchCBSfy+VCr9f7xBvqL95Bz4NhNBpZv349q1ev9tYMzJo1K+gLW63WHrXplJSUdFl9FCi5gMNVXZ7qLDhc49FVppmeni6LAg1UaTNg3CmgtMJ5L4DUFvQbiqLw3EffUvy0nZ17D6LRwG8vmslRE2QpAn86fjFq/6Wpz1LGwpmPqh2kfSjq9gjOhRGI58ukxWLx1hQYjUa/a2Q5HI6AQ1M9HUY9TQUlJSXMnTvXe9M3GAzodDqsVmuPhpI6HA7KysooKiryGREJ6v2ppyNgPAMY6urqvImF5zztO6IGis9gMPh9LwJ1ReiNoBOMxx57zJsJlZeXM2vWLJYvX+7NvHrKZDIFNVQ0UFWNJwvr6hin0+nzhnnalsQgoSiAolbJxg+B85+HuES1Y2cEWbetiOj1+qM2t5v31+1km6uJ8dpkTp02hvi4ODbvOsBv/rOK1z/5HoBpE0fy4DWzJbmIpqkLYH0lbFyqNpV4mkaiMES1/YgPD61W22mpCoPBQHZ2ts8XTE8zSqD7kWekRvtRKR3vOUVFRRQUFPRoJKXBYKC4uNjnHJ4mjGBqEDz9Etufx/OF2263+2z3F5/RaCQ7O7tTPwybzRayaRx6NQ9GcXFxxOfB0Ov13sSg4w830A9Fq9X6DL/xcLlcPrUaYgBzt8HKX6uPT38ANBp1Mq1IhqC4Kf76Me771hrR6/Y3L63ajPmZOrY4G73bJqYlM/fY8Tz/8WYOHGolMT6OW+Ydw28vmsnQxPguzibCTqNRayu2vKXO5Nm+aSRCPFOFe0ZiWCwWny+TJpOpU015TU2Nd7ZOUJMTz73Mbrd793lmBfWczzOaxDPrZnsmk4lVq3o+IqakpISysjLvc5fLxfLlywOW9xeXJ1Exm83k5OQAap8Ms9ncqTYmUHye98LpdHprN0wmEzabzVuRYDabvSNVgqVRgpgQ4q9//Su33HKL93nHhKLj/lDrOA+G1WqlpqbG2+vV0yu2fXuV1Wpl1apV3h6wVquV8vLyLoeqzps3T5pIBoLWJrD9TO3RjgZMH8PYyCaWTW2HuOIzC8/teA+ABE08rf46xv0gKS6Rdac9Pugm1Xpp1WaueHBFp/EJ7WdZOPmoMTx4zUlMT5dmrd4I2+fa+ipYcRPMeQCmdj0b5UDlmSMqVgcPRCu+fjUPRnFxMWVlZVRUVABq5tl+DgybzeatYfEwmUw4nU5vtrh7926ZB2MwOOiE1y6Bre9B3BAwPhXx5GLHoQYuWfNHPtzzJUM0iTxx7M2cpj2my5k8RyemDrrkos3txvxMnd/Bj55t2pQhvFpyNokJUmsRc6YuiEqzSLS176tQVVXV5cRe0RAL8QWVYERzHgyPriYwKSws9PsmxtoPXoTZvk3wynnQ8CUMSYXzX4T0syIawroDmznffjsbmraRljCCl2b9kTlpxwEMugSiO++v2+nTLOKPq7GZD7/ZxZwZ4yIUlRBdM5vNZGVlUVhYGJPrhMRCfEENRPbMg/Gvf/2LNWvWeOfBeOyxx8jKymLRokXhilOIntn1KTx/ippcDEuHS1dEPLkAsGyoYkPTNvTJE/hg9t+9yYXobJsr0KyQvSsnRCR4ZvusqKiIybWtYiG+Xs+DsWHDBm8NRjTnwRDCx74N6gRaumPgwtdgxKSohPHPGb8kKS6RO6deyZgh2qjE0B8oisI3W/f2qOx4rcxuKmKHZ9hqrIqF+KI2D4YQYZFxCZxrVRdbSuq6z1AoKYrC0p0fMG/MKWg0GlLik3h45o0Ru35/9N3uA/z2yVpeW93FBFCoHT0n6lI4ddqYyAQmhAiJXs/VOmvWLObPny/JhYi+zx+Gfd8efq6/NKLJRYu7leu+uJ8fr/kjd9Y/FbHr9ncbd+7ntdVbSIyP49KTJqFBTSba8zy3XJ7lnRZcCNE/BPUXu2jRIo466qhwxSJEcBQ3rPwtvPtLePm8rpeGDpM9LQe4wH47/97yOnHEMVaaQ7q0p7HZ+/j06eP4y8JZrLzrPJ781RyeumEOE3W+64hM1KXw1A1zuCQnOk1dQojeC7qJ5J577gm4b+/evYwcKTPriQhoOwS2K6G+Sn0+41pIiOwiV5uadnCB/XbWHviWYfFJVB5/OxeOmR3RGPqLg81t/HXp51TYvmbFXeczZYw62dmN5x9eZO6SnElclJXudyZPIUT/E9RfblZWFrm5uQH3t19Upj/zLHbWfqEgEUMOueDlc9XkIi4RjM/ArJvVmQUjpG7v18z+6EbWHviWCUN1rMi5X5KLAN77agen3vEaZUvX4mpsofqDjQHLxsfFMWfGOPJOmcKcGeMkuQiBJUuWMG/evK4XOxMiDIKayXPNmjXeaVczMzPR6XQ+i7WYzWa++eabsAQaSTKTZwzb/x28cj44P4fEEXD+C3DE3IiGsKflAFNWXIGrdT/HDp/Cq4Y/MylJ5rboyHWgmT9UreHxt9T5ccalJnHfldnMy54U9kn5RGfyuSYiLagmkrPPPjvgSmsul0s+NET4vftLNblImQAXvQqjT4x4CKmJw/j7tEX8d9tbVB3/e1KjtCR1LFtau5nfPVnrnbviqrMyuSt/FtphQ6IcmRAiUoJKMPR6PcuXLyc11f9aAAsWDL7pYkWEnVUOb7bAmY/AiCMjdlm34mbboQYmJo0C4Ofp53DFRCNxGqnC9+eTjU62uZqYOn4ED15zEqdPlxk4Reh5VvWuqKhAp9N5J5cCdSkJzxTZnrWoPDzN+aNGjUKr1aLT6TCZTN6FxOx2O5WVld7zmkymXi/4NagpQbDZbH3a319cfPHF0Q5BtNfwdVQvf6C1Sbls9Z3KlHd+pmw76IxqLLGqrc2tbHc1ep83HWpV/vHqF0rTodYoRiXaG8ifawaDQSksLOy0va6uzmd7XV2dYjAYlJqaGp9y9fX1islkUvR6fY/OK3omqK9fc+d23dbd3X4hgvbpg7BkOqyLzvwSOw41cHZtMc/veI/vDzmx7+3/fYxC7astezjnLzXMv+9tWtvcACQNiefG82eQNEQWJxPhF2itDYPBQGZmpvd5Xl4eFoul06qier3ep/aju/OKngmqicTfMNTFixezZ88eDAYDZ599dkiDE4OY4oYPS2C1ugouO2ph2hURDeGrA5u4wP57NjRtQ5c4ghdP/KOsKdLOoZY27nt5Lfe+/AUtbW6GJyWwdrOLE6bIh3K/09jFHDLx8TA0qWdl4+IgKbl3ZcPE02fQ0ywSaMlyo9EoTSAhFlSCYTabeeSRR3y2FRQUALB69Wruvfdebr755tBFJwantmZ48xr45hn1+ey/gKEkoiG84/yUS9fcSUPrPvTJE3jV8GemDZPJnjzeX7eDGx//mHXfq+uInD8rnfuvzOaIUdLhtV+aPTzwvjkXwMPLDj8/ayw0BVj9NvtMePztw8/PmwINu/yXPSYbnl0VbKQ9YrPZ0Ov16PV672raVqs1YHLhMVCmWogVQTWRKF2MaJUpw0VINO+FZReoyUVcApz9BGTdFtE5Ll7fVUtu3a00tO7j5NQZfDj7H5Jc/KCpuZVfP/Ex5/7Fxrrv9zI2NYknf3U6lb8+Q5ILETOqq6s7bXM4HD7NJf50l4CI4HRZg7FhwwZWr17t8/yFF17wm2isWrUKu90uNRii91qb4MUzYdcaSBgG5z0Hk8+NeBg5qUejTx7PcSMyePLYYpLjh0Y8hlg1JCGOzza5APj5mZnc9ZNZpMnQ0/7vo/2B98V36Efz9o7AZTtOjPa/jT0v20e1tbWUlZWxe/duqqqqpDYiBnSZYGRkZOByubDZbJjNZjQaDTU1NZ3KabVacnJyePTRR8MWqBgEEpLhyAvhwPdw4aswNitil25T2ojXqB+kusSRvJtzH6OHpMowVGBrQyMjU4YwbGgC8XFxPHTtbHbuPcicGTL0dMBICaL2KVxl+yg7O5vi4mIAcnJyOu3X6/XU19d3eQ7PCuEiNLr99Jw1axa33HILtbW1FBQU4Ha7O/1zOp28/vrrZGRkRCJmMdC0rxE76S7I/ySiyYWrZT/n1t3Gw5sOz3I4dmjaoE8u3G6Ff735Ddm3LuPu5z/zbp+enirJhYhpRqOx0wgQk8mEzWbr8rju9ovg9PgT1GAwcM4554QzFjEYOV6EpUZo+aHTmEYDKeMjdvlvm7Zz+se/YblzNbetfxxny96IXTuWfbVlD+fdbePXT6xib1MLH63f6R2CKkSs02q1nWac9ky2FSiJcLlcMiw1xIL6ivbGG28MiuXaZbGzCPn8UXh9Pmx5Ez57MOKXr9v7NSd/dBNrD3zLxKGjeCv7r+gSB/dqwIda2ih94TNOu+M1Pvh6J8OGJlD2syxev91IQvzgrtHprwbDYmdOp7NH5aqrqzGbzZ2SDJfLRUVFBSaTqVfnFf7Jcu1+pKeny6JA4aQo8PEdUPcX9fnMAjjxdxEN4eUdH/CTT++m0X2I44ZnsMxw16BfsGztZhc/f+g979DTc06YyN9/nsOk0TI6pD9buHAhCxcuZN68edEOJeQ8U4U7HA4cDgdarZbc3NyAo0EMBgN1dXWYzWZqamq8U4UD3v4bgHeqcIfD4V3Is6SkxO86XCKwoFZTXbx4Mfn5+QGTiOuvv77TPBn9kaw6GELrq2DFTTDnAZiaB20t8HYBrPuPuv+kP0HW7yM6DPWhTUu58auHcePmnFFZVJ/we0YmyE30e2cj2be+QtKQBMp+ZmD+7CNlAcMBRD7XRKQFVYORk5NDaWkp4H+5dukgI3w07oC3i6DZBe8UwthseOd62Pw6aOLhzHKYeW3Ew2pyH8KNm+vSz+fhGTeQGBd0Rd6AsXqDk1kZarvzRF0KS359BsdNTkM3XIbmCiH6RpZrF+GhKPDOImjZpz5v3qcmFztrISEFzqmCKRdGJbTfHWni2OFTOHdU9qD9nd3mauLmp2p5adVmXrj5LIzHTwTgzJmR62ArhBjYguq1pdfraWhowOl0dvrndruZP39+uOIU/c36KtjwAiht6nOlTa25OOE3cMlbEU0uth9q4JrP72Nvq7ougkaj4bzROYMyuXC7Ff791nqyb32Fl1ZtJj5Ow1ffy8gZIUToBVWDYbFYSE1NDbjf32p0YhBq3KHWXqAB2nfx0cCae2FGQcRCab9gWZP7EEuOvy1i1441X2/dy43//piV69SZGA0ZOh68ZjbHH5kW5ciEEANRUAlGoOXYPSuqdhziIwYhn6aRjv2HFbWp5N3r1WnAw+wd56f8eM0fcbXuJzN5IndmXhn2a8aq8pp13LZkNc2tblKGxHOH6QSuP+do4kM8XbMQQniEpHebZ0XVkpISbydQMUg516pNI4EobeB4Xi2nOyZsYTz9vY1r1t5Pi9LKKakzeWnWHxkzRBu268W68doUmlvd5B4/gb/9PIcjx3SxeqYQQoRA0AnG888/T01NTacJSAJ1/hSDjO4YyLgUNrwE+Jn5URMPGZeELblQFIU/O57hD/VPApA37gz+c+wtg27Bsn1NLXz1/R5yMkcDcEnOJF4tmcvp08cOyr4nQojIC6p+9LnnnuO6665j9+7dNDQ0oCgKaWlpKIpCamoqlZWV4YpT9BcaDcy4Fr/JBRoYMhLOCN9cKc6WfTz63TIAiqcs4Nnjbxt0ycWrq78jp2QZpvveYefeg97tc2aMk+RCCBExQdVg2Gw2b83F6tWraWho4OyzzwZgz549PP/881x22WWhj1L0Hwed8O4vAuxU4MxHISV8M2aOGjKSZbP+zMd7v6LwiOgMg42W7a4mbnm6jhc+3gRAxtjhbG1oYszIpChHJoQYjIKqwTAYDN7Her2e8vJy7/PU1FR2794dushE//TZg7B/E4zMhCkXq00ioP6vvwymLgj5Jb9t2s6rOz/2Pj9xZOagSi4UReGJt9Whpy98vIn4OA2/vnAGH/7lAhkhIoSImqBqMNpXr6amplJfX8+3337LkUceCai1GAOBZ7Ezzxz+IghZvwc0MOUiGHYE/HeaOpNnmJpGavd8zcWr/0BD6z7eyv4rp2hnhvwasayl1c2P//oW7365HYBZU3Q8eM1JnDBFVoUUqiVLlrBkyZIBvdiZiE1BJRiKorBgwQJWr17NN998w6233orBYKCsrAxFUVi1alW44owoWeysD+LiIecPh5+fVX54LZIQN420X7Ds+OF6JiWNCen5+4PEhDiOmjCC2vpd3D7/eH5xzjRZ9VT4GMiLnYnYFlSCUVBQgE6nIz8/HwCTyYTT6aSgoACNRkNdXV1YghQxrnkvrLkPDCWQ0KG9f+qCsDSL/HPTS9z01SO4cXPuqGyqTrh90CxYtqp+F6NHJJExVh1qeueCE/n1hTOZIkNPhRAxJOhhqh2nAy8sLKSwsDBkAYl+RlHUBc3WPwu71sAFL4X1cm1KGzevq+Dvm9S5NgrSz+ehQbJg2b6mFv5k/YRy29ecOWMcS81no9FoSE0ZQmrKkGiHJ4QQPgb+p7IIry//pSYXmniYZQ775Z7+frk3ubjnqGspnrJgUAy9/N+aLfzmiVV852wEYEJaMgdb2kgeIn/CQojYFLCx9t577w36ZL05RvRjuz+H925UH8/+C0w4NeyXvGKikcsnnM2zx9+GOSN/wCcXO/Y0cdVD75F3/zt852xkyphhvHjLj6goOlWSCyFETAuYYNTU1AR9st4cI/qplgPwRj60NsGkc2HWLWG71PrGLRxyNwMQp4nj6eNuJX/8WWG7XqxYs9FJ9q3LeO6jTcRpNNx4/gw+vPtC5h43IdqhCSFEtwJ+BaqpqWHUqFFBnczlcvU1HtFfrLgRGr6AlPEw90nQhGfkwtvOT7h0zZ1cMDqHp4+7dcDXWLQ3Iz2VsalJTB49jAevmc2sDBl6KoToP7qsY83IyECn8/1Qq62tJTs7u1PZ3bt3s2bNmpAGJ2LU/i2w4XlAA8ZnwjYz51Pf27j2hwXLNjRtZ39bEyMSUsJyrUhrc7t5f91OtrmaGK9N5tRpY3C74akVDq6YoycxIY6hifG8cPOPmJCWLENPhRD9TsAEw2QyUVVV1Wn7rbfeyj333OP3mAULQj8csaOysjLvomoul4vi4uJuj6moqPAuxlZfX09JSYkszNYXw9NhwRr4bjkccXbIT68oCnc5nuH/BuiCZS+t2oz5mTq2/NBhE2DMyKEkD4ln065GXAea+e1F6oRhk0YPjqG3QoiBJ2CC4ZnroqOuqqgDHRMqnuTCMyzWZrNRVFTkM2W5v2MKCwt9kpKCggKqq6vDGuuAN+JImHFNyE/b7G6h8Iu/85/v1f485in53H3U1cSFqQkm0l5atZkrHlyB0mH7zr2HABielMARuoFRSyOEGNwCfmp3nO+iJ3pzTDBKS0t9akmMRiMVFRVdHlNZWelTW6HVaqWvSG99/AfYuCysl7j8s3v4z/c1xGvieHTGjdxz9LUDJrloc7sxP1PXKblob0RSIvNPnhyxmIQQIlwCfnLv3bvX73ZFCfzxGOiYUHA4HN5mjo5sNlvA4/R6Pbm5ud6kwuFwoNfrwxTlAOZ4AWrvglcvhoZ1YbvMLyfNY1TiSF4+8U8UTboobNeJhvfX7fRpFvFnq6uJ99ftjFBEQggRPgETDLPZ/6RJXTWRBDomFBwOh9/t3dVILF68GKfTSVpaGmazGZvN1mWTivBj37fw1g/NISfeDGnTQnr6g23N3sdn6U5gw5wnOX/MSSG9RizY5moKaTkhhIhlXQ5TfeGFF0hNTfXZbrfbeeuttzrVZLhcri5rEsJFp9PhdDoD7tdqtRQVFVFTU0NZWRlGo5EFCxZ02cnTs5qqx6BeVbWtBd74CRxywbjZ6oRaIbR0xwdc/+UD1GTdw8zh6qq8A2WkSHuO7ftYvSHw72l747XJYY5GDAaeVVQ9ZDVVEWkBEwyHw4HJZOq0XVEUv4mEoihBzVFgtVqprKzstlxJSQkGgyHg/q6SC1BrVXJzc6mursbhcJCXl0dWVhb19fUBj5HVVNv5+Pew/UMYkgq5z0J8YshO/eCmF7npq0dQULj/2+d47JjfhuzcsaKl1c2D//uK0hc+41BLG2NGDmXX3kN++2FogIm6FE6dNvhWhRWh1/GLkaymKiItYIKh1+uxWCw9Hs7Z0NBASUlJjy9sMpn8JjBdxeOPy+UKuM/Tb8NoNHrPUVdXR25uLlarNajrD0qb/gery9THP/o3jJwSktMGWrBsoLE7dvOrf3/EZ5tcAJw1cxwXZ0/i5idr0YBPkuFJzS2XZxEfNzA6tQohBreACYbRaAx6VEg4m0j0ej1ardZvJ01PAtGRw+HwmyDl5eWFI8SB57s31f+P/SVkXtbjwzY17WBXyx6/+5ramvmT42ne2F0HDMwFyw4cauUu6yc88sbXuBWFtGFDuPunBi4/PQONRsO41ORO82BM1KVguTyLS3ImRTFyIYQInYAJhsViCfpkvTkmGCUlJdhsNu88GFar1WepeIfDgdVq9U6+ZTQasVgsnUaf1NXVSUfPnji1DCacDpPO6fEhm5p2MG3l1Rx0t3RZbogmgSePKx5wa4q43Qrn3FXDp5saAMg7+UgsP8tizMgkb5lLciZxUVZ6p5k8peZCCDGQBEwwOnbu7IneHBOM4uJiysrKvHNf1NfX+yQKnhEi7Wf3rK6uprS0lFGjRnlHnIQ7ERpQMoJrt93Vsqfb5ALgkRk3DrjkAiAuTsO1c4/i3qWf87ercjj3hHS/5eLj4pgzY1yEoxNCiMjRKH4mtti4cSNTpkwJyQVCea5ImTdv3uDt5LnlbVj9Vzj78V6tMWLf+w1ZH/6y23J1Jz+EYeRRvQgwtiiKwrPvb2TsyCTvKqdut0JTSxvDhspy6iJ2DOrPNREVfutk6+rquPfee/t88sWLF2O32/t8HhEhTTuh5qew6VWwl0Y7mpi3Ycd+fvzXtygs/4Ab/v0R+w+qNTdxcRpJLoQQg57fBGP+/PnMmjWL/Px83nzzzaBPunz5cvLz8xk1ahSXXdbzzoEiihQ3LL8SGrdC2gyY/edoRxSzWtvc/OPVL5l92zLe/HwbQxPjuObsoxiaEB/t0IQQImYE/Jo1d+5csrOzMZvNmM1msrOzyc3NxWAwoNPpGDlyJKBOD+5wOHA4HLzxxhvU1taSk5NDRUVF2PtkiBBac686LDU+Cc6pgkRZxdOfNRud3PDvj1izUe3EecaMcfzj6hymjh8Z5ciEECK2dFmPm5qayqOPPgrAc889x7PPPsvdd9+Nw+Fgzx51GKJWqyUtLY2srCwWLFjgLS/6kW0fwEe3q4/nPACjjo1uPDFq3fd7OOuPr9PmVoee/vkns7jiDP2AGmIrhBCh0uOG4vnz54d9tVQRBQcb1KnA3a0wNR9mXNen09U3fh+iwGLPtImpXJx1BPFxGsp+lsXYVJnSWwghApGeaINd0w6IHwIjM+GsCujDt/F9rY3c8NXD3ZZLiktkdGLsN5/t2neQPz/3Kbddepw3mfjX9acyRPpaCCFEt/wOUx3ssrKySE9PHzyLnDXvgwPfh2SV1Me+e437v32Oh2b8itQE//04RiemMjk5+CGwkaIoClUfbMT8jJ3d+w5hOvlIHv/FadEOS4he8Sx6tmXLFurq6qIdjhhEJMHwY1CMF29rCenCZe21uFtJjOuflWPf7tzPr59Yhe2zrQAcM0nLg9ecRE7m6ChHJkTfDIrPNRFTZG7iwah5H1TPgk8fgD7ml/tbmyhY+zd2Nru82/pjctHa5ubB177kpJJl2D7bytDEOP5gOoEVd54nyYUQQvRC/7sTiL5RFHhnETjXqkNTp/0chvauP8TBtmZ+vOaPLHeu5ssDm1iRc3+/HVHx4P++4g+VawA4ffpYHrj6JI6aIENPhRCityTBGGy+ehy++S9o4iF3Sa+Ti1Z3Gws/u5vlztUMj0/m/mlF/Ta5ALju7KOofH8ji3KP5sozMomL67+vRQghYoEkGIOJcy2s+JX6ePafYULvOi66FTfXfXE/L+54n6FxiSyddScnpU4PYaDh99bn26j8YCOPXDcbjUbDiORE3r/rfEkshBAiRCTBGCxaGuGNfGhtUpdfn1Xc/TF+KIrCb9Y9yn++ryFeE0fV8b/nR7oTQxtrGO3ed4jbltj573sbADhjxlh+eroeQJILIYQIoYAJxqhRo9DpdBiNRjIzMzEajZx44okRDE2E1Hs3qTUYKeNh7pOg6V3/3vu/fY4HNr0IwBPH3My8saeEMMjwURSF6g++xfxMHbv2HUKjgcK5R3NR1qRohyaEEANSwARDURSqqqqYNWuW3/179uyhoqKCzMxMWdAs1ikKaI+GuCFgfAZSxvX6VJeNPY1HNr/Cr4+8lJ9NNIYwyPDZtOsAv3niY974VB16OiM9lQevOYnZR42JcmRCCDFwBUwwsrOzAyYXoK5Tcsstt7B69WoWLFhAZmYmJSUl3kXQRBStr4IVN6nrikzNU2fnnHULHH05DJvYp1NnpEzgk1MeYVhC/5km++qHV/Lx+l0MSYjDfMmx/PrCGTIbpxBChFnAenK9Xt+jE8yaNYuqqirWr1+PTqcLWWCilxp3wNtF0LQN3i6EvRsO7+tlcvHC9vdYuuODw6fpR8kFwD2XG5gzYywf/OUCii85VpILIYSIgIAJRschh8899xyLFi2ipKSEN998s1P5xYsX43a7Qx+h6DnPHBct+9TnzXug8gTY/VmvT2nbbecnn5Zy2Sd3srJhbYgCDZ+m5lb+r2oN97/yhXdbTuZolt06l6NlXgshhIiYLvtgtOdZTXXq1Km4XC4cDgdGo5EpU6YA6rLtRmP/aJMfsNZXwYYX2m1Q1GTj62fglHuCPt0Hri/48Zo/0qy0YBo3h5O1sT0U9Z0vtnHT4x9Tv30/QxPj+MmpU5ioSwE6J8xCCCHCK2CCEegD2Wg0YrFY/Pa16GmzSqzbsmUL8+bN61+LnTXuUGsv0AAdpv/+ohxO+C2k9HyBsU/3ObjA/nsOtB3knFFZPH2cmXhNbDYtOPcf4vYlq3l6hQOAiWnJ3Hdljje5EGIwa7/YmRCRFLCJxGaz8cILL7B3716f7WlpaQE7cmq12pAGFy3p6eksXbq0/yQXPk0jftYWad4H717f49Otb9zCOXUluFr3c6p2Js+f8AeGxg0JXbwhoigK1g83kn3rMp5e4UCjgYK5R7Hqnou4KOuIaIcnRExYuHAhS5cuJT09PdqhiEEmYA1GfX09JpMJUGsmcnNzMRqNOByOgCeTaugoca7t0DTSgdIGjufVcrpjujzVzmYXubW3sr25geOH63ll1l0x26nz+4YmFi3+kEMtbqb/MPT0ZBl6KoQQMSFgguFpCqmpqcFms/Hoo4/y6KOPAmC32zEajZxzzjnMnTvXW6PhcrkiErToQHcMZFwKG5eqyURHmnjIuKTb5AJAlziC80bnYHPaeSOrlLTEEWEIuPcURfEmsum6FO6YfwJNza385sKZDE2MzSYcIYQYjDRKx96cP1i8eDEFBQU+21avXo3NZqOmpoba2lpcLhcajcbbwdNut/PNN99EJPBwmjdvHkuXLo12GMFp3AHPHAUtezvs0MBQLSz8qsd9MBRFwdmyj1FDYmvUxdrNLm58/GPuXjhLJskSIkj98nNN9GsB+2B0TC5AnfPilltu4Y033sDpdFJfX88jjzzC3Llzqamp6bL5RIRZylg4q8LPDgXOfLTL5OJgWzOWDZW0uFsBtakrlpKLg81t3Fn9Caf/4TU+Xr+L3z+7JtohCSGE6EafFjvLyMigsLCQwsJCQJ39U0SBoqizdU5dAOsrDzeVeJpGpi4IeGiLu5UFn/6Zl3d+yKf7HDxzfEkEA+/eii+3c8O/P6Z+uzq3x8VZR3DvFfJ7JoQQsS6kq6nKPBhRUvsnaNwOJ9+t1lZseQuaXTBkJJzxSMDD3Iqbq9fey8s7PyQpbghFR1wYuZi70XCgmd8/u5on36kHYLw2mfuuzGZetixOJoQQ/UFIE4x77gl+MifRRw3roO5ucDfDEXMhcz6cVX54LZIATSOKonDDVw/xzNY3SdDEYz3hDs7QHR/h4AN7bfV33uTi2rOncueCE0lNib2hskIIIfwLaYIhIkxR4N1fqMnF5PNB/8OqtlMXdNksAnDH+id4ePPLaNDw5LHFXDhmdgQC7lqb2018nNotaOFpGXzw9U5+cloGp03r+QRhQgghYkPATp6iH/jmv7DlTYhPgjn/VPth9MD9G638ZcMSAB6ecQMLJ/wonFF2q83tprxmHbNve5W9TS2A2tH0wWtmS3IhhBD9lNRg9FcHG2Dlb9XH2XdAas+naT9hRCbD4pO4PWMhiyZdFKYAe+aL71z86l8fsap+NwCPv7Wemy6YEdWYhBBC9J0kGP3VR7dB0w5ImwEn3hzUoXNHzeKLUx9jUlL05pI42NzGX5d+zt+WfUlLm5sRSQncueBErj37qKjFJIQQInQkwfAj5hc7O9gA9Vb18ZmPQnz3nR9tu+2kDx3NjOGTAZicHP6mhza3m/fX7WSbq4nx2mROnTaG+Lg43vtqBzc+/jHfbFUnBbvQcAT3XZlNuixOJkTIyWJnIloCzuQ5mPWLGe+adsKGpTDz2m6LrmxYS27draTED2XlSX9j2rDwD/V8adVmzM/UscXZ6N2WrkvBcnkWL63aRPWH3zIuNck79FTWsREivPrF55oYUKQGo79KHuOTXGxq2sGulj2diq078B2FX/yNJvchztIdT0by+LCH9tKqzVzx4IpO67p+72zkigdX8NC1JzF65FBuu/R4tMNk6KkQQgxEkmD0J/s2w8460P/YZ/Omph1MW3k1B90tAQ/VoOFv0xYxJC4xrCG2ud2Yn6nzt2g8CqAB/vLC56y9f553SKoQQoiBRz7h+5P3boT/XQof3eGzeVfLni6TCwAFhQNtB8MZHQDvr9vp0yzSOQ7Y4mzk/XU7wx6LEEKI6JEEo7/YsBQ2vAhxCTA1P9rRBLTN1RTSckIIIfonSTD6g5YDsOIG9fEJv4NRx0Y3ngAOtbRh+/T7HpUdr00OczRCCCGiSfpg9Ae1f4L9m2DEkeqkWjGqpc3Nu19u77KMBpioS+HUadGbg0MIIUT4SQ1GrNv9GXxyv/p4zkOQOCy68XSwt6kFt1vt0jk8KZGHrjuZG86fjgY1mWjP89xyeZZ08BRCiAGu39VgOBwObDYb1dXV1NTU9OiYsrIytFotAC6Xi+Li4jBGGEKKAu9cD+5WdSGzKbGznDqow1FvfqqWkh8fyzU/zMB59rETOPvYCcyeOqbTPBgTf5gH45IcWXJdCCEGun6VYNjtdmpra3G5XDidzh4d40kuCgsLAbDZbBQVFVFeXh7OUENDo4Gs2+HDW+H0fwQs9r7riwgGBVsbGvndk7W8XPcdAE++W8/VP5rqM1nWJTmTuCgr3e9MnkIIIQa+fpVgGAwGDAYDVqu1x8eUlpayYcMG73Oj0Uhubm7/SDAAjjwfJp8XcKXUz/dt4Nav/9XtaZLiEhmdmNqnUNxuhcffXs8fKtewt6mFhHgNv7lwJsXzjvU7E2d8XBxzZozr0zWFEEL0T/0qwQiWw+HA5XJ5m0fas9lsGI3GyAfVlfVVsOImmPMATD4XhoxUtwdILrYdcnLh6js44D7I7NTp/GPa9STG+f+Rjk5M7dP6I99s3csv//URH3ytzl+RnTmKf14zm2MmaXt9TiGEEAPXgE8w/NFqtbhcroDHeRY784jIomeNO+DtImh2wVvXABo45R449hf+i7cd5JLV/8emgzs4KiWdZbP+zChPQhIG+5pa+OibXQwbmsD/5Z1AofEoae4QIoZ5FjnzkMXORKQN6AQjEJ1O12UfjvT09MguCqQo8M4iaNmnPm/Zr/6/+7OAh7y+q5aP965DlzgibMmFp+8EgEE/in9eexJnzRzPpNGxNZJFCNFZxy9G7b80CREJUUswrFYrlZWV3ZYrKSnBYDCE9No97SAaMeurYMMLnbePyQl4yKXjTufZ429jwhAdRw1LD2k4exqb+WPVJzy1op73/nQ+09PVvhtXnJEZ0usIIYQYuKKWYJhMJkwmU1ivodfr/W53uVwB90Vc4w619gINdFwi7IPfwZSLIOVw3wlFUbwdKvPHnxXycF6u28zvnqxla4M6lffrn3zvTTCEEEKInhrQjeh6vR6tVuu3L0ZMdPD0aRrxs/5o8z5493rv07ecazjl45vYcnBXyEPZ2tDI5Q+s4Kf/WMHWhib0Y4fzyq1nc9MFM0J+LSGEEANfv0wwAjVxOBwOysrKfLaVlJRgs9m8z61Wq3dOjKhzrlWbRpQ2//uVNnA8D861rDuwmcvW/ImP9nzFPRu6b1oKxtMrHOSULGNp7Wbi4zT89qKZfHj3BZw5c3xIryOEEGLw6FcJhsPhwGw2U15ejt1ux2w2U1FR4d1vs9k6zW9RXFyMy+WioqKCiooKVq1aFTtzYOiOgYxLQRPvf78mHvSXsWv4EVxovwNX635OSZ3JX48uCGkYO/ceZE9jC1l6HSv+dB53LjiR5CGDsv+vEEKIENEoiuKnbn5wmzdvXuRGkTTugP9Og+Y9+DaTaGColkM/+ZS5a+9jpWstGcnj+fCkfzB2aFqfLtnc2sbWhiaOHDMcgNY2N8++v5GFp02RoadCDFAR/VwTgn5WgzEgpYyFMx+lcx8MBeWMR7im/ilWutaSmjCMV2bd1efk4qNvdnL6Hf/jsnvf5lCL2jSTEB/Hz+boJbkQQggRMlIPHm2OFw43lWxcqva70MRDxiXcm6Dhv9veIkETj/WEO5g5/MhuT9fmdvtd/2NfUwt3Vn9CxfKvURQYPWIoX2/dy3GT+5awCCGEEP5IghFNB3fjfvMqNC37qT/zEaZ8ZyO+ZR9ticP4wnALMzUapqVM4ndT5mMc1f1cIC+t2txpBdN0XQr5p06h8v2N3u0/PT2DuxcaGDViaNhemhBCiMFNEowo2vthCSOb97JmyHCyNlczXzeFf+xax426KVjX/AGAoXGJnDsqu9tzvbRqM1c8uKJTQ8sWZyP3v6Kutpoxdjj/uOokfnSsjA4RQggRXpJgRMue9Qz/6gkAbh59NG6NhuoR46ke4XvzP+RuYVfLni4XKmtzuzE/U+dvJg2v4UkJrLzrPEYkDwlB8EIIIUTXpFdftHxYQpy7hddSRrE8ZVSfTvX+up0+zSL+7D/YypqNDX26jhBCCNFTkmD44VlNtf1KhCG17QOot6Jo4rhl1NF9P52rKaTlhBADx5IlS5g3b56spioiTppI/AjraqqKAit/B8Bu/WWs1ezp8yldjc09KudZGVUIMXh4VlWV1VRFpEkNRqQpbpiaB8OOYOsJN/XpVPuaWrjlqVp++5/aLstpUEeTnDptTJ+uJ4QQQvSUJBiRFhcPJ/wGfuagJWVcr0+zeoOTk0qW8WjN1wCcPn0MGtRkoj3Pc8vlWTKRlhBCiIiRO060xCf26fAjRqXQ2NzGlDHDePGWH/Habbk8dcMcJupSfMpN1KXw1A1zuCRnUp+uJ4QQQgRD+mBEyiEXLLsIZhXDlItBo2F0YipJcYkcdLcEPCwpLpHRiakoisKbn2/j7GPHo9FoGDMyiRduPotp6akMG6r+GC/JmcRFWel+Z/IUQgghIkkSjEipuxu2rYQPS+DIC3ATx+Tksaw77XF2teyhze3mk29d7N53kFEjkjjhSC3xcXGMTkyleU8SF/39Td79cjtP33A6l+RMBsCg7zy8NT4ujjkzet/0IoQQQoSCJBiRsHcjfPaA+viUMv71fQ3/3fYWlcffzuTksaz+/JDfKb7vXjiLN3fsxPLi5xxsaSN5SDzO/T0bMSKEEEJEkyQYkfDR7dB2CNLP5oPUKVxfewstSivPbH2TKTuyA07x/fOHVnqfn33seP5+1UlkjB0e2diFEEKIXpAEI9x21MI3/wVgZ9btXPbJXbQorcwfezq/PGIex/7t5S6n+NZo4OFrZ3P5HD0aTccxIkIIIURskt5/4aQo8P7NALQd9VMu3vIi25qdHDt8Ck8cewsffL2r2ym+FQWOHDNckgshhBD9iiQY4bR1BXz/Dkr8UMxpGXy05yvSEkbw4ol/ZHhCskzxLYQQYsCSJpJwmjAHLljKe9++zH3Oj4kjjsoTbiMzZSIA41KTenQameJbCCFEfyM1GH6EbLEzjQamXEyaoQR98gTKjr6O3FFZADi276Pspc+7PhyZ4lsI0Tey2JmIFo2iKF31MRyU5s2b17fFzlr2Q1szJOm8m/a1NjI8Ppk2t8KD//uKu5//jIMtbSTGx9HS5kYDPp09PT0uZBZOIUQo9PlzTYggSQ1GONjvQXk6k/rV93g3jUhIobVNYe6f3uAPlWs42NLGWTPHUXvPhTwtU3wLIYQYYKQPRqisr4IVN0HO/6F8cj+a1ibMG6u4eKyBn6efA0BiQhynTx/Hhh37Kf2pgZ+enoFGo0E/boRM8S2EEGJAkQSjDzY17WBXyx40jTuZufxahrTtx73iBuKVVt5L0vJiyhhO3ZrIl4qLGUdoAbj9suP49YUzGDPSt4OnTPEthBBiIJEEo5c2Ne1g2sqrOdjWjHXbJxzXth8NEKe0AvC70UfjRsMdj33DKxOaqbkjl/i4OFKGJpAyVN52IYQQA5vUwffSrpY9HHS3sGD/duYf2OnN1DydM49saULRuNEkHeKkqaNpaZW+tEIIIQYP+SrdS21uN2Nam3l055e48c3U3MCjO7/k7WQd9xaczJXHZEUpSiGEECI6JMHopU82NvDIzi8Y4W49nFy0qP/FASNppXzLF+wZdQAaD0B8PAxt1++i8UDgk8fFQVJy78o2Narzi/uj0UBySu/KHmwCtztwHCnDelf20EFoawtN2eQUNW6A5kPQ2hqasknJ6vsM0NIMLS2hKTs0Sf29CLpsi1o+kCFDISEh+LKtrep7EUjiEEhMDL5sW5v6swtYNlEtH2xZt1v9XQtF2YQE9b0A9W+iqYsp/IMpG8zffU/Ltv+bECKGSYLRS5pdXzL/wE7fjY8dfpgAXMpO4EfqhjkXwMPLDhc4a2zgD6bsM+Hxtw8/P28KNOzyX/aYbHh21eHnP54J33/rv2zmTHhx7eHnC3Og/gv/ZSceCa9vPPz8qjNgba3/smmj4d1278X150PtO/7LJqfAx+0+OH8zH1a86r8swGftEqCSK6DGGrjsR/sPf/jeWQRL/xO47Ds7QPfDBGZlv4XKhwOX/d8GSJ+iPn7gdnji3sBlX/gcph6jPl58NzxyZ+CySz6GY3PUx0//A+4vDlz2329BzlnqY2sF3P2rwGUfegXOuFB9vOwZuOPqwGXvrYJz89THy1+AmxcELnvX4/Djq9TH778Ov7wocNnb/gkLf6k+tq+Aa34UuOxvy+DqW9THX9ph4UmBy17/f/CLP6qPHV/CpccGLnvVzfC7v6qPt26C8zICl83/Bfz+IfVxwy44c2zgsvN+Dn95Qn3c1Aizu1jhONcE91cfft5V2Z5+Rnwmza2if5A+GL2kjJ7Bc8PG0MX3XiGEEGLQkpk8/ejJjHerXOu48P0i1m1aSaqnmeSHGm434IqLxzDpFJ6b/QBZqUdLE0mgstJE0ouy0kQCSBNJkGQmTxFp0kTSS/FxcexMGMKiMTOo3P6ZuvGHz9M4YNG4mXybkoxm2HD/HwjBfEgEUzY5pfsyvSmbFMSCa8GUHdqzBd+CLjtk6OGbQCjLJg45fNOKWtnEwzfvUJZNSDicbISybHx8z3+HgykbFxeeshpNeMpC+MoKEYOkicSPnix2NjoxlaS4RKqGj/NpKmlBg3XYWKpHjCcpLpHRiamRCVoIIfyQxc5EtEgTiR89rUr0mcnzpXMZ0raf5vgRfHHJ/1BSxjA6MZXJyV10FhNCiAiRJhIRadJE0geTk8eqCcTIo2Duv2DFTQyd8wCzxp8a7dCEEEKIqJImklCZuoAlSffD1LxoRxJzumpqGuzkvfFP3pfA5L0R/YUkGCEkf/j+yfsSmLw3/sn7Epi8N6K/kAQjxoTqwyPWzhMKsfSa5H0J/3lCIZZeUyy9L0JEgiQYMSaWPhBDeZ5QiKXXJO9L+M8TCrH0mmLpfREiEmQUiR/HHHMMmZmZQR+3ZcsW0tPT+3TtUJwj1s4TS7GE6jwSS3jPE0uxhOo80Y6lvr6etWvXdl9QiBCRBEMIIYQQISdNJEIIIYQIOUkwhBBCCBFykmAIIYQQIuQkwRBCCCFEyMlU4T1UVlaGVqsFwOVyUVxcHJZj+qPevjeg9mwHKC8vD1t80dLXn39ubi41NTVhiCz6evvemM1m7wgvnU6HyWQKV4hR05v3pqKiApfLhVarpb6+npKSEu85hIgaRXTLYrEo5eXl3uc1NTVKYWFhyI/pj3rzOouLi32eFxYWKkajMSzxRUtff/7V1dXKQP3z7M1709DQoBgMBqWhoUFRFEWpq6sbkO9Pbz9rPO+LoqjvlclkCleIQvTYwPsLDQOtVuvzB6woSrcfbr05pj8K9nU2NDQoRqPR5xjPzaK+vj5MUUZeX37+DQ0NSnl5+YD8fVGU3r03xcXFisVi8dlWU1MT6tCirjfvjcFg6LRtoCXson+SPhjdcDgc3qrHjmw2W8iO6Y96+zpra2txOBze53q9HlCrgweCvv78q6qqWLBgQRgii77evjdlZWUYjUYcDoe3nNFoDFeYUdHb90av15Obm+v9+3E4HN6/KSGiSRKMbrS/Eban1WoD3hB7c0x/1JvXqdVqaWhowGAweLd5PjwHyodiX37+NpttwN042+vL35PnBqzX6ykqKhpQyTr0/vdm8eLFOJ1O0tLSMJvN2Gy2AdmnSfQ/kmD0kk6nw+l0hv2Y/ijY11laWkp5efmA75TWk/fFcwMdbLp6bzw3Xq1Wi8FgQK/XY7FYyMvLi2SIUdPd741Wq6WoqAiTyURZWRnV1dUD6ouM6L8kweil3iQKgyG5gOBep9lspqioiMLCwjBGFBu6e18qKioG5KiInujqvdHpdABkZ2d7t3m+1Q+0Wgx/uvu9MZvN6PV6qqurqa+vx+l0kpWVFaHohAhMEoxuBPo22dU3zd4c0x/19XVarVYyMzMHXHLRm/fFbrf73EAHqr78PXW80Wq12oDNCv1Rb94bT7ORp1lNr9dTV1eHXq/HarWGLVYhekISjG7o9fqAH2SB2sp7c0x/1JfX6fnm6UkuXC7XgLlZ9OZ9cTqd2Gw2ysrKKCsrw2w2A2rnxoF0o+jNe6PVatHr9Z2OcblcAyop681743A4/DYtDpbmIxHbJMHogZKSEp+qWKvV6vOt2+FweCeO6ukxA0Vv3hu73Y7dbsdgMOBwOLDb7ZSWlnqrwgeCYN8Xo9FIcXGx919RUREAxcXFA67ZpDe/MxaLxWfSMavVitFo9OksPBD05vfGbrd36nNRV1c34H5vRP8jy7X3UPvZ9err67FYLN59FRUVWCwW76yUPTlmIAnmvXG5XGRkZPjthDbQfhV78zsD6k2lsrLSe3PJy8sbUDVf0Lv3xjNbJcDu3bvl7+kHLpeL0tJSRo0a5e2bUlhYOOA7TYvYJwmGEEIIIUJOmkiEEEIIEXKSYAghhBAi5CTBEEIIIUTISYIhhBBCiJCTBEMIIYQQIScJhhBCCCFCThIMIYQQQoScJBhCCCGECDlJMAYgs9lMbm4ueXl5mM1mysrKvDMgeqagdrlc5OXlkZWVhUaj8W5vz2azkZubi0ajITMzs9P0zR4Oh4OioiLvtTwrpAZaWyQtLY3MzExyc3MpKiqiqKiItLQ00tLSyMvLo6ioiLy8PDIzM0lLSwvNmyK8bDYbWVlZpKWlhWw10mDPabfbSUtLw263e597jq+oqAhYTgjRjyhiwKipqVH0er1SXl7eaZ/FYlEsFovS8UdeU1OjGI1GBVBqamr8ntdoNCoNDQ1+95WXlyt6vV6pr6/32d7Q0KAYjcZOsTQ0NCh6vb7T+QwGg6LX63tUVvRdQ0NDlz/zcJ+zrq5O0ev1Sl1dnc92rVbr8zsTqJyiKJ1+54QQsSUhuumNCBVPbUNNTY3fdSuKi4sDrrBosVi8NQf+1sbIzc31u66BzWbzHtNxOWmtVkt1dTUZGRnA4VVTnU4nRUVFPVonQavVYjabcTqdsq5CiGm12pC/p8Gc02Aw+P1d67jgXaBydrsdh8MRcBlzIUT0SRPJAJGXl4fJZOpyUazFixf73e5JBhwOh3eZ8J5es7i4OOCHvFarxWKxUFRU5G2icblcQa2AmZ2dPWCWcRehE6qmHSFE+EiCMQB4+lj460fRnlarDXhzNxgMFBYWUlZW1qP2bs818/Pzuyy3YMECAJ/EJTs7u9vze+j1+gG1jLvoO7vdHlQiLISIDmkiGQAqKysBerSkd6BaDIDy8nKqqqooKCigrq6uR9fsrjZCq9Wi1+upqqqivLw8qNoLz/E9OcZut1NZWUlmZiYulwutVuttlrHZbJjNZhwOB8uXL6e2ttb7+srLy7vd7+FwOCgvLyczMxPovJR2T8/TXvtjFi9e7K2t2b17Ny6XC4vF4m12CFWc7blcLp9OlXV1dT7X9OhYJjc3F5PJ1KtzOhwO8vLycDgcWCwW78+pI3/lrFYrNTU13tfseWyxWLy/Y3a7HYPBwOLFizEYDLhcLrKysnA6nSxYsCDgz8LDZrPhcDioqamhuroaq9WK0+mkvr6enJycgK9bCNFBtDuBiL7TarWdOm/2VE1NjU9nuZqaGgVQLBaLd1v7x725psFg6Lasv06ePeXp3NrxfO07C3o6IBYXFysNDQ1KdXW1zzHd7a+rq1MMBoPPNerr6zt1Qu3uPIEAnTrTVldXK1qtNqjz9zRORVF/hkaj0WdbXV1dp2tWV1crJpOp07HV1dWdXkdPz+kp27ETsL9Oyv7KAX6v73l/Ou6rrq72W94fz++7wWBQTCaTT9yAdC4VooekiUT4MBqNmEwm7zflUPD0vwgXzxDZjtvaf2tv3wFRq9ViMpl8Og92t9/TCbY9vV6PwWDwuXZ35wlEq9WSl5fnU3NgMpnQ6/VBnb+ncbYv357BYCA7O7tT2Y7NZkaj0Vt70Ntz9rTpK5gmMk/NVcdaCofD0aOaB7vd7q0JtNvt5Ofn+/xM9Hq99P8QoockwRgAQt2TfvHixWi12i77dARzzXCOAvGMJujYryNQ51BPs0Eg/vZ7ruGvCSo3N5eqqqqgr+OPvxup0Wjs8fl7E6e/axoMBp+y7ZMYl8uF3W7H5XLhdDp7/Do6njOcioqKvM0cHsGMbjEYDN5jO76XDocj7AmzEAOFJBgDgKejZU++WQWaLKs9z+gPm82G1Wrt0zVdLhcul6tH/UN6w3MjsNlsVFRUeP/V1tb6bWvv7tuwv/2eawQ61vMag7lOT40aNcrvDS1Ucfrj6cfSntVqJSsri4KCgl4ljP7OGS6eJMHz86+oqPB2Nu4pm82GwWDweZ2Bkg4hhH+SYAwAhYWFaLXabjuvgdp5sKfnNBqNFBQU+K3i7+k1S0tLAQJ2MuwrT02K0WiksLCw07+Oursx+tvvuYa/b+yem2bH40JVY7N7926/5wpVnP7U19f7lKuoqKCgoIDq6mqqq6sxGo1BJ1Adzxlq7TuVglqL4dnm6fQbDH/zyVit1h53OhZCSIIxIGi1WhYvXozVau1yiKnVavU7rDTQN8vy8nJcLpffqu2eXNPlclFWVobFYgnbhEieb5m1tbWd9oWqrdxzDX/nW7VqVchGFfhLDKxWa4+/ffcmTn/XtNvtPtf09Gdp/zNsf1zHm3tPztlXHROGjr/DhYWF3qHbvalxsNls5OTk+GwrLy8PW6IsxEAkCcYAYTKZqKmpIS8vr9MHPuBt6uj47cvTbu+PXq/v8gPVZDJRXV1NXl5epyTD4XAwd+5cLBYLxcXFwb6coFRXV2OxWHxuMi6Xq1MnxO6aCLraX11d7U24POx2O3a7vdPQ3542Rfi7Rnuen2PHn0Go4gQ6DUf29F3o6ufeftK09v/35pz+XkdPtmVnZ7Nq1SoA77DUjgoLC71NHcHw9LNo/3dRVlbmnStGCNEzGkVRlGgHIULLbDZjt9vR6/XezoBGo9Hng9blclFQUOCt9jUajZ1ucB55eXkB93nO5RkhkJmZ6Z3DwWw2d1lzYbfbKS8vx+FweL91G41G9Ho9eXl5QX3z9JzL83o7zoNRXl6O1WpFr9djNBo7zS/R1f5A19i9ezclJSVBn8eftLQ0b5Kk1WpxuVzU19f7NEGFKk4Ps9mMxWKhrKzMe83du3d3SgTsdjulpaXk5OSg1WrR6XSYTCby8vLQ6/UUFRV5f849OafnfJ7X4allaL/NZDKRn5/fqZwnWfUssJebmwvgN4m12+3U1tYGnRRUVFRgsVi8v5seklwIERxJMISIAWlpaSxevFgmcQohT+fOYPtf5OXlodPpetSnSQgRmDSRCCEGhPY1YdC7zp1weOl5IUTfSIIhRAzoal4J0TNWq9Vb62C1WnvVpOHpfyFDUYXoO0kwhIgim83m7UdgsVhkEa8+KCwsRK/XU1FR4TPjaU/ZbDbvLKSefkxCiN6TPhhCCCGECDmpwRBCCCFEyEmCIYQQQoiQkwRDCCGEECH3/2l+YzbkZBw+AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 400x300 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "from matplotlib import pyplot as plt\n",
    "import mpl_config\n",
    "\n",
    "\n",
    "plt.plot(error_probs, e_array[0], marker=\"o\", linestyle=\"--\", label=\"Noisy 1-layer\")\n",
    "plt.plot(error_probs, e_array[1], marker=\"s\", linestyle=\"--\", label=\"Noisy 2-layer\")\n",
    "plt.plot(error_probs, e_array[2], marker=\"d\", linestyle=\"--\", label=\"Noisy 3-layer\")\n",
    "plt.plot(error_probs, [uccsd.e_fci] * len(error_probs), linestyle=\"--\", label=\"FCI\")\n",
    "\n",
    "plt.xlabel(\"CNOT error probability $p$\")\n",
    "plt.ylabel(\"$E$ (Hartree)\")\n",
    "plt.legend(loc=(1, 0.5))\n",
    "plt.savefig(\"error-probability.pdf\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4fd7b167-5c84-4025-b481-fc6991d921a2",
   "metadata": {},
   "source": [
    "## Effect of measurement shots on VQE energy uncertainty\n",
    "Here we show how measurement shots affect the standard deviation of the optimized VQE energy.\n",
    "As expected, the standard deviation is in linear relation with $\\sqrt{\\frac{1}{N}}$.\n",
    "In general, increasing $p$ causes a larger standard deviation of the energy, yet the effect is far less significant than adjusting $N$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "aefacfff-3d2f-4082-a295-17b96f8fa6b9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "11.673604011535645"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "hea = HEA.ry(uccsd.int1e, uccsd.int2e, uccsd.n_elec, uccsd.e_core, n_layers=1, engine=\"tensornetwork-noise\")\n",
    "hea.kernel()\n",
    "# cnot error probabilities\n",
    "error_probs = [0.10, 0.20, 0.40]\n",
    "# number of measurement shots\n",
    "shots_list = [2**i for i in range(8, 14)]\n",
    "e_list = []\n",
    "# record run time\n",
    "time1 = time.time()\n",
    "for error_prob in error_probs:\n",
    "    for shots in shots_list:\n",
    "        hea.engine_conf = get_noise_conf(error_prob)\n",
    "        hea.shots = shots\n",
    "        e = []\n",
    "        for i in range(64):\n",
    "            e.append(hea.energy(engine=\"tensornetwork-noise&shot\"))\n",
    "        e_list.append(e)\n",
    "e_array = np.array(e_list).reshape(len(error_probs), -1, 64)\n",
    "time.time() - time1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "f2ecbfce-1205-49d8-95df-83c85cf90e08",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAX8AAAEiCAYAAAARVNJOAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/NK7nSAAAACXBIWXMAAA9hAAAPYQGoP6dpAABVxklEQVR4nO3deXxU5b348c9kA8KSyYQ9smSCoqCVbO4skkStC6gkUKt4b3ubRG9re1s1A+1te9vai5Pb5XfrbSWJXVExyQDirpmooLiQZEDrhjATFtkCJJMAScgy5/fHYYZMMpNMJjOZLN/368VLzplzzjzHhO8885zn+X41iqIoCCGEGFHCQt0AIYQQA0+CvxBCjEAS/IUQYgSS4C+EECNQRKgbMNjMnz+fxMTEHo85fPgw8fHxPl3P12ODcU05Vo6VYwfu2FC/f0/HWq1WPv30U/edinBzxx13BOSYvh4bjGvKsXKsHDtwx4b6/Xs61tN+GfYJsnvuuSegxwXz2L4YDO2Ve+u7wdDe4XpvQ+2+pOffRaB7/kPJcL0vRZF7G6qG670N9H1Jzz9AgvUJH2rD9b5A7m2oGq735td97SuFv06DfWUBaYNGUWSFb2fLli3jhRdeCHUzhBDigqZaeHYutNphlBbu2QPRk30+3VNck56/EEIMZooC2x6AttPqdutp2P5gvy8rwV8IIQazfaVQswWUDnVb6QDbZnV/P0jwF0KIwaqpVu31o+nygga25amv+0mCvxBCDEZuwz1dH80q/R7+keDfxeHDh1m2bBkbN24MdVOEECNZ3afuwz1dOYd/6j71/DqwceNGli1bxuHDh7u9JrN9upDZPkKIQUFR4LW7oWYr3Xv+gCYcEpbDLZt6vZTM9hFCiKGivRnam/AY+NFA1ARY9KTfl5fEbkIIMdi01MHLt8Px90ETCUpblwMUWLy+T3P9u5KevxBCDCYtp2DLDWrgH6WF5RWQcJc6zAPqf/V3w5yV/Xob6fmLQaOgoACtVguA3W4nPz+/13NsNhtms5mysjLKy8sDck0hQmpULExcAK2NcPtrEHc5aOfC4bfUFb79HO5xkuA/SHQ4HLy35wTH7M1M1Y7hurmTCA8bOV/MnEE6NzcXALPZTF5eHoWFhV7PsVgsVFVVYbfbqaurC8g1hQg5TRgs/Rs0n4Bx53PzR0+GJYXwzg9g4R/6NdzjehuZ7eMuFLN9tlYewvBMNYfrmlz74nXRGO9NYXnajAFtS6jExsZSU1Pj6qUDaDQafPn1NJlMrFu3jurq6oBdU4gBVfOCOq3zxj+rwT/AZLbPILS18hCrn3jHLfADHKlrYvUT77C18tCAtqeoqAiDwYDFYsFkMmEymTAYDNhstqC9p81mw263uwVpJ7PZPGiuKURQfPYUvHYXfPE3+PwvA/a2MuwTQGfPtXt9LVyjYXRUuNuxHQ4Hjz5d5XEil4K6oNvwTDU3Xj7F6xBQmAbGRF34MTadayd6lH8/VrvdTmpqKjqdjuzsbKqrq9FqtWi1WgwGA2VlgUkl25W3DxatVovdbh801xQioBQFqh+DnT9Tty/9Nlz6rwP29hL8A2hqjvdESzddOZ1NDy9xbeu/u4mmVi8r985TgMN1TVz2H1tpbO461UuVnKBj2y9ucW2nrX2ZT3+3vE/tdrLZbCQnJ1NSUkJWVpar16zT6XrsLefl5fl0/ZSUFNf4uy90Op3Hsfz+CMY1hegzRwe8+3345E/qdspP4KpfgaZrDp/gkeA/BDgGaIw6OTkZUB+kGgwG136bzYZer/d6XrAeoAYjSEvgFyHX3gLm+8C2CdCoD3Cv+N6AN0OCfwAdK/Y+7za8yye67Y8r2LGnlhW/ebvX6/79ezdw/VzPT/fDunQUKtfd1uv1elNVVUVqaqpru7y83G070Lx9sNjt9h4/dAb6mkIERN0ncOAlCIuCjKdhTnZImiHBP4DG9mGsfeyoCNIvn0q8LpojdU3eFnAzXRdN+uVTfZ726e94v5PNZkOn07k9KC0tLe02k6az/g776PV6tFqtx28YGRkZvjV8AK4pREBMToXMjeoCrvgbQ9YMCf5dOLN63nPPPUGvHxoeFobx3hRWP/EOGtwzeDg79MZ7UwZ0vn/Xsf2CggJyc3ODPuyzdu1azGaz68PBZDK5fVDYbDZMJpPHRVrehnJ6u6YQA6Z+Dzja1AVbAPq7BuRtN27cyMaNGyWrpy9G+jz/vLw8tFotaWlprqA6UAGz82pcq9WK0Wh0vVZUVITRaMRqtbr22Ww2CgsLMZvNWCwW8vPzSUxMdGtvT9cUYkAc/xBevg3CR8Hd78H4WQPeBE9xTYJ/F6FK6TxYVvgmJiZSVlbmevgrhOiHA6/C61lqds7JaXDbyzBm0oA3w1Nck2GfQSI8LIyFl00JdTNc0z2FEP30xd/hrX9Ti67MuBluMUHkuFC3ykVW+AoXs9ksM2GE6C9FAYsR3vxXNfBfch/c+uKgCvwgwV+cZ7FYXOPhRUVFIW6NEEPYp4XwwRr17wsegfS/Q3hkaNvkgQz7CEBd4OUpJbIQoo8uvgc+K1J7/At+FOrWeCXBXwgh+qu9BSJGq38fFQMrPoDwqNC2qRcy7COEEP3RdBw2Xwe7fnNh3yAP/DAIe/7+VF7q6Ry73U5pqZpwzWq1YrPZKC4u9pjqVwgh+qRhH7x4MzTa4OwRmPcddeXuEDCoev6dKy/l5uaSnJzca+qA3s4xGAykpqaSm5uL0Wh0pSsWQoh+qa2GzdergX+CHu5+d8gEfujHIq/du3e7SuidOnUKgLi4OPR6PXq9ngULFvT5mv5UXurtnMzMTDIzM13fBgoKCli3bh319fUerxeqRV5CiCHkUDm8dje0nYGJSXD7KxA9NdSt8qrfi7z279/P448/TmlpKXFxcSQlJbklAdu3bx87d+7EZrPR0NBAVlYWeXl5zJ49u9dr91Z5yVMyLl/O6TqDpbKyUhJ7CSH8t/c5qLhfzdVzUTrcslktqj7E+Bz8H3jgAWpqasjNzWX9+vU+nVNRUUFubi4pKSmsW7eux2P9qbzU13NMJhN2uz1oFamEECNAy0k18M9ZdX4O/6hQt8gvvQb/hoYGcnJy+PGPf9znoZz09HTS09PZtWsXq1atori4mAkT+vYJ6U/lJU/nFBUVYbfbe01d4Mzq6TQQ2T2Fyt+H/YAr4VvXDKP+XFOIHl3xPXWMf+YtQSm2HgjObJ5OnrJ6ovSiqKiot0N81tO1ysvLFU/N0Wq1SmFhYcDOKSwsVLRarVJfX+/x9TvuuMNrG4PhQNNxpbrhS69/DjQdH9D2hIrRaHT7mZWXlyu5ubk9npOfn++2nZubq2RkZPTrmkJ0096qKB/8RFGaT4W6JX7zFNd6Df4DxWq1egzkgFJeXu7XOfX19Up+fr5boHeeU1ZW5vGaAxn8DzQdV0aX36rweqbXP6PLbx0RHwCePpB76pvU19crGRkZbudUV1crgGK1Wv26phDdtJ5WlBdvUZQ/oihbFiuKwxHqFvnFU1zr13eWp556ijlz5hAeHs5TTz0FqOP8a9eu7fO1Olde6srbA9rezrHZbBQUFLgNATmfBQyGef4n2xpocXguzO7U4mjjZFvDALVIHR4zGAxYLBZMJhMmkwmDweD1+Uog9Pbg3puqqiq3djmT0tntdr+vKYRL8wnYuhQOvgYRY2DBowNaYD3Y/A7+xcXFrF+/HoPBQFVVlWt/eno6a9ascX0Y9IWz8pKTp2pOzjFeX85JTk4mPz/fLVNlSUkJycnJQZnxc7a92euflo7Wbsc2d5zz6brNHee8XrfrNZo6Wvxuv91uJzU1lbS0NLKzs8nIyCArK4vMzEy3gu6B5s/Dfq1WS319vdszHOfvgV6v9+uaQrg07lfn8NdWwigdLHsTZve/PvZg4vcKX7vd7hb0O9d4jYmJ8TqPvif5+fkUFBS4skparVa3B3hms5nCwkK3h3a9nbN27Vq3Dwy73U5FRUWf2+aLcW8u9/rarROv4uXkx1zbk99eSZPDt+B/+66fYW8/4/G11AmXUHnN/7m25+3IYf+iDT622J0zl39JSQlZWVmuXrNOp+uxt9zfGr7e9PVh/7p16ygsLOzxW50/EwjECHPyI3jpFmg6BuNmwh2vQ+yloW5VwPkd/GNjY3t8XePn16OeZmM4V/H25RytViszPHzk7EVbLBa3nr6nIuidBaKGryd9CdIGg4G8vLxeP1wk8IseKQ6o+Bc18Osuh9tfg3HxoW5VUPgd/Pft2+e2rXRahdvY2Njt9ZHgzNKtXl8L14S7bdcuKWX3aSs3VPae8vWlpF+yYHyix9fCukw1++z6Yh9a2rOqqipSU1Nd2+Xl5W7bgebtg8Vut/tUXMZkMnWr3dvfa4oRShMGN5XAB2th6V+GVLqGvvI7+K9atYqLL76YNWvWkJKSQkNDgyvlg9FoHJELqcZGjOnTsWN8XBwyJnyUz9eODh/tcxs8sdlsbqu2AUpLS92G9brq77BP5wf3XQNzb89mnMNRzuva7Xbq6ur6dU0xAjXsg5g56t9j58LXN4e2PQPA7+CflJTEk08+6Vr56+z5x8bGUlpa6lduHxF6Xcf2CwoKyM3NDfqwj/PBvTOIe3rYbzKZ3IbwLBYLFouFrKws1+yekpIS12yz3q4pBIoCO38Gu4xqcfUZmaFu0YDxO7FbZxaLhZqaGvR6PUlJSYFoV8gMZGK3g821zN3xrR6ne44Oi2TP9X9l5pjJA9KmvLw8tFotaWlprvHxgQqYnVfjWq1WV1lJUKegGo1G10peu91OQkKCx5k7nX+le7qmGOEc7bDtQfj8/MzEax6H5ODNagslT3EtIMF/OBnorJ4Hm2t7nMc/MTJmwAI/QGJiImVlZb2mwRBiSGtrgvJ7YP8L6jj/oidh/vD9Vugprg2aRV4j1cwxk0mecLHXPwMZ+OHCdE8hhq2WOngxUw384aPh5k3DOvB7M6gWeYnQMpvNMhNGDG8tdbDlBjj2njqTZ1k56O8MdatCwu/g71zklZOT022c399FXoOBM6tn54x4I4HFYnGNhzsXzAkx7IyKhclpMDYe7nwHpt0Q6hYF1caNG1m2bJnHrJ6DbpFXqMXHx4/ISl7JycndCt8IMexoNLDkKWg5AWOnh7o1QedMSd85Tb2T3z1/WeQlhBgSal6A8vvA0aFuh0eOiMDfG7+Dv3OR15///Gd2797tWuT11FNPkZKSwgMPPBDIdgohRN999hS8dhfsfQY+k+HMzmSRlxBi+FEUqP417Pypun3pt2FeTmjbNMj4HfxBXSa/b98+du3a5VpGP9QXeQkhhjhHB7z7ffjkT+p2yk/gql8Nq1z8gdCv4O+UlJQkQV8IEXrtLWC+D2ybAA0s/INac1d04/eY/wMPPMDFF18cyLYIIUT/2L+Ag69AWBTc9JwE/h70q+f/+OOPe32tsbGRCRMm9OfyQgjRNxMXwE2lEDkW4m8MdWsGNb+Df0pKCpmZ3jPgGQwGnnzySX8vPzLtK4V3fqB+VZ2THerWDLjOSdjsdnufi/BkZmZ2W6vQ32uKIaB+D3S0wMQr1e3Zt4e2PUOE38E/LS2NdevWAWoysM454O12uxTJ7qumWng7D1rtsC0Xpi+G6IHN6xNKziDtzCBqNpvJy8vzOV20yWTymI66P9cUQ8DxD+Hl2yAsEu5+DyYkhLpFQ4bfWT11Oh12u91jvVS73Y5Go6Gjo6O/7RtwA53VE1Cnpb22Qk00pXSAJhwSlsMtmwa2HSEUGxtLTU2N2++TRqPBl19Pu91OaWkpeXl5bsf355piCDjwKryeBe1NasqG216GMZNC3apBKaBZPfV6PfX19dTV1XX743A4WLFiRb8bPGLsK4WaLWrgB/W/ts3q/gFWVFSEwWDAYrFgMpkwmUwYDAZsNlvQ3tNZiMVTR8KXb5ClpaWsXLkyoNcUg9wX/4BX7lAD/4ybYfmbEvj7yO9hH6PRSExMjNfXfS3tN6y0nfX+miYcIkZ3P7a5FrblARqgc49Uo+6fvkR9eOX1umHQucRjWxNERve97ag96NTUVHQ6HdnZ2VRXV6PVatFqtRgMhqCV5vT2waLVaj0Wa+nMbDZ7LMvYn2uKQUxRYPf/wPvni65cch/c+GcIjwptu4agfo35d1VcXExDQwPJycmkp6f3q2Gh4szq6UyI1CfF47y/NvNWuP3lC9t/naz2WrxSoPU0bH8QjmyHlpOeD5uUCtmVF7afmwer9/el1S7OXP4lJSVkZWW5es06na7H3nJ/a/h6o9PpXNXEvHEWZPc1oPtyTTGIff7UhcC/4BG41qh2gIRHGzduZOPGjYHN6ulpNk9Ojrp8eteuXfzmN7/hkUce8ffyITOosno6h3+itAPyds4iLhaLBYPhQjk7T0XQOwvWA9TegnRRUVGfS0xK4B/i5nwDPiuGOatgwcOhbs2g11NWT7+Df08PzZKSkqioqPD30kNXzhnvr2nC3be/Vat+hTXfCwdevjDe3/WchOWQ/o8ertul1/ONz3xvrxdVVVWkpqa6tsvLy922A83bB4uzV++JxWLpsU3+XFMMUu0tED5KTc8QNR7u2qFm5hT94nPwr6mpYdeuXW7bW7Zs8fghUFlZicViGZI9/37paWze27FLiuHZudDaQLcx/6gJam3RPl3Xv/F+J5vN5jZtF9QHqtXV1V7P6e+wj16vR6vVevyG4Wk8H9QevMVicQ1HOQu7FxQUoNfrXcNWfbmmGISajqtTOfV3Q8qP1X0S+AND6QOLxaIUFBQoGo1GCQsLUzQaTbc/sbGxyk033aTYbLa+XHrQuOOOOwb+Tb98TlH+SPc/e0sGvCmFhYWKXq93bRuNRiU/Pz/o72s0GpXCwkLXdllZmZKbm+vatlqtitFo9Hp+dXW10vXXubdrikHOvk9RNiSq/xb+MklRmk+FukVDlqe45tc8f4vFQlFREevXrw/Cx1FojfR5/nl5eWi1WtLS0lzj430dV/dX59W4VqvVVVYS1PF9o9Ho6uF3ZjKZKCkpwWQykZubS3Z2tqt339M1xSB2wgIvfV2dDTchAW5/HbSSS8xfnuKa34u8Nm3aNCzn8ock+IO6wvfZueoK31GxcM8XIVnhm5iYSFlZmevhrxAD7pBZLcDSdkbN1XP7qxA9NdStGtICushrOAb+kIqeDEsKYcxUWFwYstQOzumeQoTE3ufg5VvVwB+/FO7cJoE/SCSl82AyZyV862jIkrqZzWaZCSNCq7URHG2QuBJuf0Wd9CCCol+rI3pL6SyGDovF4hoPLyqSWqciRObnqsM8N21Up3eKoPE7+PuS0lkMHcnJyZSXl2O1WgfsAa8QdLTBBz+B5k4r2GfeIqt2B4CkdBZChEbbWXg9Gw6+Cke2wV3bJegPIL+D/9KlS3tN6SyEEB41n1QXb9XuVBMTJq+VwD/AJKVzF87Ebhs3bgx1U4QYnhr3w5br1cA/SgfLKmD2baFu1bC0ceNGli1b5jGxm9/z/CsqKnrM3Nnb64NVyOb5CzESnPwYXroFmo7CuBnq4i3dZaFu1bDnKa75PezjLbA70zpnZWX5e2khxHCkKPDWt9TAr7scbn8NxsWHulUjVsAH2XJycnjkkUekTqoQwp1GAzeVqEna7twugT/E/O75A2zevJny8vJuOdK9PQgWQoxA9i9Be4n695g5I6o29WDmd89/06ZNfOc73+HUqVPU19ejKAqxsbEoikJMTAwlJSWBbKcQYqhRFPjwp/DcfDjwSqhbI7rwu+dvNptdPf5du3ZRX1/P0qVLAWhoaGDz5s3cfffdgWmlEGJocbTDtgfVsosApz6GWbeGtk3Cjd89/87Jv/R6vdsYf0xMDKdOnepfy4QQQ1N7s5qi/POn1Ln7iwsheU2oWyW68Dv4d17EFRMTg9Vq5cCBA659DQ0N/WuZEGLoaamDFzLU2hTho+DmTWq+HjHo9KuG78qVK9m1axd79+5lzZo1JCcnU1BQgKIoVFZWBrKdQojB7pwdtiyE+s9glBZufRGm3RDqVgkv/A7+OTk56HQ6Vq1aBUBWVhZ1dXXk5OSg0Wh6rPkqhBiGomJg2vVqPerbX4O4y0PdoiHvYHMtJ9u8j6JMjIxh5hj/an/4vcJ3uJIVvkL0g6Mdmk/A2GmhbsmQd7C5lrk7vkWLo83rMaPDItlz/V97/QAIaCUvIYSg5gV4fZUa9AHCIiTwB8jJtoYeAz9Ai6Otx28GPfEp+P/mN7/p84X9OUcIMYR89pRaa9daCp+uD3VrRB/5FPzLy8v7fGF/zhkMJKunEL1QFKh6DN7OAcUBl34L5j8Q6lYJD3rK6unTA9/y8nLi4uL69KZ2u71Pxw8W8fHxMuYvhDeODnj3+/DJn9Tt5B/D1Y+peXtEQNW3ne73Ne655x7uueceli1b1u01n2f7JCQkoNPp3PZVVVWRmpra7dhTp06xe/fuvrdUCDF47CuFd34AC/8Ac7KhvQUqVoPVBGjghv+Frz0U6lYOKy0drbx44gP+caScV08Gd7q8T8E/KyuL0tLSbvvXrFnjtYj7ypUr+9cyIUToNNXC23nQaodtuTB9MTQfhwOvQlgUZGyAOfJvPNDeqtvNyo8fG5D38in4O+fyd9VTqUZv5wghBjlFgW0PgHPYofU0bH9QzcZ5yyYIi4SLloa2jcPA3rOH2XDUTGzkOH44S618mBmXQuqES8iMSyZl/MVkffyroL2/T8Hfn5KMQ7WMoxAj3r5SqNlyYVvpANtmdb/09vvlVGsjJcfeZsPRCj5o+ByA+FET+f7MOwnXhBMRFk7lNf8HqPP8R4dF9jrPf2JkjF9t8Sn4NzY2MmHChG77e1of5u0cIcQg1lSr9vrRAJ3/fWtgWx5MXwLR/q0oHclePbGTosOv8PKJnbQp6pqIMMK4eWIKq6dloCio/8s7mTlmMnuu/2vQVvj6FPwNBgNPPvlkt/09Dft4O0cIMUi5Dfd07dgp7sM/okfOjrEzRr508kOer30PgKTxc1g9PZ17pt7I1FE6r9cA9QPA3+DeG5+nem7ZsoWYGPevFxaLhbfeeqvbNwC73Y7ZbA5cK4UQwVf3qftwT1fO4Z+6T0E3f+DaNYRYm47w9NEKNhyp4K+XP8zC2CsA+Lf4WxgXPobV09K5fHxCiFup8in422w2jwXZFUXxGOQVRenxW0FPCgoKXCUg7XY7+fn5/T6noKAAAKvVCiD1hYXoSlFgz9M9H6MJh4TlEvi7qG87TemxbWw4WsEO+6eu/c8efcsV/JMnXEzyhItD1USPfAr+er0eo9Hoc13e+vp61q5d2+fGOIN4bq6a/9tsNpOXl9djsO7tHIPBgNFodB2fl5dHZmbmkF2BLERQaDTQfkb9e1jk+Vw9Xcb8oybAIhnKdapvO03Op7/nxRMf0qqoD2XDCCM9bgH3T8vgzsnXh7iFvVB8kJeX58th/T5Hq9Uq9fX1bvt6a2JP59TX1ysZGRlur1dXVyuAYrVaPV7vjjvu6HO7hRiSmk8pSuPBC9vnGhRl/8uK8uVzivJHuv/ZWxK6tg4CDodDOdh83LXd4ehQZm27T+H1TOWKHbnK/9SUKoebT4awhd55ims+9fw795x91ddzbDYbdrvd47cLs9lMRkZGn89JTU2lqqoKm83mKjup1+uBoZt+Qoh+UxSwlsE7D0HsPFj+ptrzj5qg1tlVFNhXolbjUjouDPeM0Gme+5uPucbxT7Y1cHTxc0SFRRKmCWP9vO8zbZSOK8cnhrqZfeZT8O/6oDcY59hsNo/7tVqt10Dd2zlarZb6+nq315zPKJwfAkKMKKcPwfZ/hwMvqduj49SVu9FTLxyj0cDi9XD4LXWF7wgc7mloO0vZ8e1sOGpme/0/XfvHhI3in2dqSJlwCQC3TEwLVRP7rdfgv3//fmbPnh2QN/PnWjqdjrq6uoCds27dOgoLC70+v3Bm9XRyJkYSYkhTHPDJk/DBGmg7o47rp/xELawePqr78dGTYUnhhdw+I2hu/zNHK/i3T3/HufOLqzRoWKpbwOpp6dw95QbGR0SHuIW927hxo1tmYr+yelZXV2MymXjkkUf61Zji4mLi4uL6HPz7Gvh7OsdgMJCXl+d6OOyJZPUUw87Zo/DaCjj+vro99TpYUgy6eT2fN2flsB/qURQFy+m9jA6LYv642QBcOU7POUcb88bO4v7pGdw7bSkXjZ4U2ob2UddOq6esnr3m81+xYgVJSUmsWrWKN998s8+NqKioYNWqVcTFxXH33Xd7Pc7bMIzdbvf6Wl/OMZlMJCYm9hj4hRiWRsepC7cix8OiP8Jd7/Qe+Ie5g821rLNtZP57OaR+8D1+bbvQS758fAKfXFfEJ9cVYUhYNeQCv698GvNPT08nNTUVg8GAwWAgNTWVzMxMkpOT0el0rjQOjY2N2Gw2bDYbb7zxBlVVVaSlpVFUVNTrMwC9Xo9Wq8Vms3UL3J4e9vblHOc4vzPw2+126urqZNxfDF/Hd8LEJAiPhPAoyNyoFlgfPyPULQuZxvazbDr+LhuOmHm7/mOU81NZR4dFMTosym19kvNbwHDmcw3fmJgY1q9fT2VlJRkZGTz33HNkZWUxe/ZswsPDCQ8PJyEhgaysLEpKSsjMzKSqqoonn3zS54e/a9eudVs0ZjKZ3HrqNpvNtWDL13MsFgsWi4Xk5GRsNhsWi4V169Z1q00gxLDQ2gjbvwubrobd/3Nhf9zlIzrwA6RXGfj2p7/lrfqPUFBYHPs1/jz/RxxbXMJfLn/Y74WpA6HD4eCdz49T9v5+3vn8OB0OR7+vqVGUHrKzhUDn1bpWq9VtymhRURFGo9G1Ure3c+x2OwkJCR5nC3m7bU9V7oUYEmpeUGfynD3/cG/+g7D4T6FtUwgoisLu01aePfom/5W4mrERYwAw1pTw18Ovc//0TO6dtpRZY6aEuKW+2Vp5CMMz1Ryua3Lti9dFY7w3heVpvn2ge4prgy74h5oEfzHkNB2Dd76vzt0HmJCoztS5KD207Rpgh1tO8szRN9lw1MwnZ/YD8PQVBu6dpv5/aHO0E6EJH9Q9/K62Vh5i9RPvdEuz57yDDQ8t9OkDwFNc87mMoxBiEDrwKpi/Cefs6mKsBQ9D6s8hcvBPRwyEs+3NbKpVx/Er6na7xvFHhUWybNK16MdMcx0bGTa0wl2Hw4HhmepugR/UxBsawPBMNbenxBMe5vMIvsvQ+r8hhHAXkwjtzTApGZY8BZOSQt2iAVXbaudfPrnwbGOh9nJWT88ge8oitJHjQtiy/ntvzwm3oZ6uFOBwXRPv7TnBwsv6PoQlwV+IoaSjDY5sgxnnZ7NpL4E7t8GkFBhiPdu++ufpGv5xpJyG9rMUzf8hAAnR0/iX6ZkkjpnGfdPSSYie1stVho5/Hqzv/SDgmL3Zr+sP798WIYaT2mp4+ztw8iO4+z2Yeo26f8rVoW1XEB09d4pnj77FhqNmPjqtpnOJ0ITz64u/xaQoLQB/u/zRELYwsJrOtfNC1SE2bLex/fPjPp0zVTvGr/eS4C/EYNfWBJU/h49+p6ZpGBWr5uMZxt44WcXvDmym/JQFB+q0xkhNBLdPuprV09KJiRgb4hYGh6XmFDmF77u2R0WEca7d87RODTBdF811c/1bhOYx+IeHh9PR0eHXBYUQAXSoXK2d21ijbs/5Btzw/yB6aExT9JVDcdCudBAVFgnA3qYjvH6qCoBrY+Zx//QMVk5dhC5y+NQFr21o5rn39qMo8INbLwPg+rmTWTJvCtdfOplv3qBnV00dq594B+hWXQEA470pfj3sBS/BX2Z/CjEIvPsf8PH/qn8fexEsfhJm3x7SJgXaZ2cO8I8jZp45WsF/6r9J3gz1/lZNXUxtq537pqVz8dj4ELcycNo7HLzx8RE2bLfx2u7DtHcoxI6N4oHMSxgVqU5DfXHNhSm6MyeOZcNDC7vN85/ex3n+nnic5z+Se/4pKSnEx8dLNk8Rep//Bd76DlzxPbj61xA1PtQtCojac/VsPPY2/zhSjuX0Ptf+m+NSeS3lv0PYsuD58mgj/9hmZeOOGmobWlz7UxPjWL1Qzzdv0DM6Ktzr+R0OB+/tOcExezNTtWO4bu4kn3r8zuyehw8fprq62u01v4P/7t27qaqqwm63c+rUKQDi4uLQ6/Xo9XoWLFjQa8MGI1nkJULm9EF1de7Ua9VtRYFT/4SJXwttuwLEoTi4e/cveOnkh3Qo6jh2hCacWydexf3TM7ht4tWMDo8KcSuD4z+f28X/vvI5ABPHj+KeGxJYvVDPZRdpB+T9+73Ia//+/Tz++OOUlpYSFxdHUlISOp3OlVph37597Ny5E5vNRkNDA1lZWeTl5QWsHoAQw5KjAz75E3z4Y4jSwj2fqgVUNJohHfgdioOPT9ewYIJa5SpME0ar0k6H4uCqCXO5f3omq6YuZmJU34tFDVYOh8KOPbVs2G7jmzcksGS+WiRn9SI9e482snpRIjdfOZ3ICP/G6QPJ5+D/wAMPUFNTQ25uLuvXr/fpnIqKCnJzc0lJSWHdunV+N1KIYavuU3Vo5/gH6nbclepq3aih+2Bzz9lDbDhi5umjb3KwpZb9Czcwc4xaDGbdxd/m93MfYO7Y4ZVk7qtTZ3n23RqefsdGTe0ZAM61dbiC/9zpMZT8cHEom9hNr8G/oaGBnJwcfvzjH/d5KCc9PZ309HR27drFqlWrKC4udqV/FmJE6zgH1b8Gy+PgaFNz7V9rhPl5oAl9r7CvTrTaee7Y22w4UkFl4x7X/gkR0fzzTI0r+A/FWrfeOBwKW3Ye5Ol3bFR8chTnAPr40RGsuGYW9y8e3Pfaa/AvLS2ltLS0X2+SlJRESUkJxcXF5OTk9OtaQgx55xpg0zVg/0LdTlgOC/8Pxl0U2nahFjk52dbg9fWJkTGuQO5UcWoXt1h+TLuiPicM14RxS1wa90/P4I5J1zDGU5nIYUCjAePWT/j8sPr/a+Flk1m9KJFlqTMYO2rwL6HqtYWBDNYS+IUARsXAxCvV4ugL/w/0d6uRJMQONtcyd8e3aDlfu9aT0WGRPH35GiZERpMZlwLAVTFzidJEcOV4PfdPy+AbU5cweVTsQDV7QJw83ULZ+wfY/OFBns+/kbGjItBoNHz/1svYX3uGexfqSZg8tHIJ+fTxFBcXh06nIyMjg8TERDIyMobsbB4hQqLmBZicBmPP555Z+H9qFs7RgydInmxr6DHwA7Q42sj6+FdcOV7P7mvV4D8+Ipp9C//GtFFxA9HMAdPhcFDxz2Ns2G7lZcth2jrUGUrP7zzIvQvVKoD3LRy61QB9Cv6KolBaWkpSkueMgQ0NDRQVFZGYmNhjnV4hRpyzR+Hd74PVBIlZcPP5nPtjJoa2Xf0wJmwUSePn0NLR6pqaOZwC/3F7M0+W72HjuzUcqb+QNG3B7FjuW6jn60nDY9GZT8E/NTXVa+AHtcTjo48+yq5du1i5ciWJiYmsXbtWHu6KkUtR4PM/w3uPQGuD2sufMEed1hnmfTHPUFCe8jjXx84PdTOCprmtg9+++BkAunGj+MZ1s7lvkZ4rZg6eb2mB4FPw97XQeVJSEqWlpWRnZ6PT6Whvb+9X44QYkux74e1cOPK2uj0pBW58CiYuCGWrenWi1fuD3s7GDJOFWIqi8MHek2zYbuVcWwd/fvB6AGZPGscjd8znylmxfD0pnlGRQ/vD2hufgn/XsmebNm2ivLyc2NhYMjMzWbp0qdvrxcXFbNq0KXCtFGKoOGSGV25Xp3JGjIGrHoOvfX9Q59o/0nIK4/4S1h96KdRNGRBH65tcc/L3HTsNQHiYhv++J5kp59Mj/zz7ylA2cUD4PObf2YoVK1ixYgVz5szBbrdjs9nIyMhwreTVarVkZGQEvLFCDHpTr4ExU0A7Fxavh5jB+0DwcMtJjPtLKPrqFc718qB3ONj+2XH+8OrnlH98FMf5mDZ2VAR3XjWT1Yv0TI4ZHeIWDiy/ev5OGRkZGI1Gj2P7vg4VCTGktZ1VE7Bd8V11cVbkOLXQytjpg2L6pieHW07yeE0JxYcvBP3rtfO5d9pS/v3zJ0LcusByOBTCwtSfw5dHG3n9oyMAXHPxJFYv0nPXVTMZPyYylE0MGZ+Cv9lsZsuWLaSnp7sF+tjYWK8PdZ35foaaw4cPs2zZMsnqKXp38A011/7p/RA+CubnqvvHDe7ZIOtqnuOPh9QkXwu1l/PzxNUs1S3gUMsJRoet73We/8TIwZ2Lx362FdMHB9iw3cq/LE7k20svBiDrmlkcOnWWexfquWTayJiM0jmrZ1c+BX+r1UpWVhag9ugzMzPJyMjAZrN5Pcfbt4XBLj4+XrJ6ip61nIIdP4I9/1C3x82A8bNC26YeHGqppdXRTmL0dAAMCSv54uwhfqK/hyWxV7r+rc4cM5k91/+1zyt8BwOHQ2HbZ8fZsN3KC9WHONd2vvpXRJgr+GvHRvGLlQtC2MqB5+zELlu2rNtrPgV/5/BOeXk5ZrOZ9evXu5K7WSwWMjIyuOmmm9y+Gdjt9sDdgRCDgaLA3o2w4z+g+QSggSsegqsfG5S59g8217Ku5jn+fPg1bp2YxvNJvwBgxujJmFONHs+ZOWbyoAzu3iiKQsHWT/j7NiuHTl0odjJ/hpbVC/WsvG526Bo3yPkU/LOzs0lKSiIpKYn8/HwAdu3ahdlspry8nJKSEgoLC9FoNK6HvRaLJagNF2LAvfcwfPR79e+6y2FJ8YUi6oPIgebjrKt5jr8cfp02RZ1u3djexDlHK6PChv40zdb2DqIi1OmXGo2GSuspDp1qQhsdSfa1s1m9KJEFs2OH7OjDQAlYJa+amhrXNwOz2UxDQ8OQrAYmxVyEV8crYetiSP4xJOXDIJvvfqD5OP9ds5G/Hn7DFfRvjL2SnyeuZrFu6NYFALWHX207xYbtNjZ/eIAdj93KzIlqEfcde2o5UtfE7SkXMSZq8E6pDaV+F3PpSUJCArm5ueTmqg+9UlNTA3VpIYJnXym88wNY+AeYk+3+2qlP4EQ1XPov6vaUNLj/EIwenKkMnq99j6KvXgFgqW4BP9ffx6IhHvSdRc43bLfxxeELzyK27DzoVvRc9F3QPiZlnr8Y9Jpq4e08NbvmtlyYvhiiJ0N7i5prf9fjQBhMvgp0aqAZTIHf1nSU2lY712jVtuVedCvvN3zG92Ys54bYy0Pcuv45dPIs+c9Uu4qcA4yODOfOtBmsXpTIDZdKwO+voAX/xx9/PFiXFqL/FAW2PQBt6gpPWk/D9gfha/8Bb+eA/XxBkoTlagrmQcTWdJRf1zzL34+Uc0n0RfzzukLCNeGMCR/Fc1/7Saib57fTzW2uOfex46J465NjtHcoapHzRYmsuHomMdGDa6htKJMBMjEy7SuFmi0XtpUOsG1W/wBETx1UufYBrE1H+LVtI/84Wu4qgD5j9CTq2k4zKUob2sb5qbG5jU3n5+Q3tXbw/mNfR6PRMG50JH/8t6uZd1HMgBU5H2l6Df779+8PWAH2QF5LCL811aq9fjRAt/kOcMl9sPAJGKUd4IZ5Zms6yi9tT/P00QpX0L85LpWfJ97Htdp5IW5d310ocm7l+cpDNLeerwAWpsFWe4bEKeq02RXXDN61E8NBr8G/uroak8nEI4880q83Ki4uJi4uToK/CC234R4PgV8TBu1NgybwA+xtOszfj5QDcEtcKj9PXO0a5x9qtlYe4qclu1xFzgEumTaB1Yv03HN9giuxmgi+XoP/ihUrqKioYNWqVeTl5XXL4NmbiooKioqKWLVqlRR6EaFX96n7cE9XikMd+qn7FHShyVm/9+xhPj27nzsnqymGb4pL4eFZWWRPWcjVQyzot7R2cK69wzVWPzoqjJraM64i56sXJZKWGCdz8kPApzH/9PR0UlNTMRgMGAwGUlNTyczMJDk5GZ1O51rV29jYiM1mw2az8cYbb1BVVUVaWhpFRUXExAyuh2ZihBo9EaKnQ9MRz69rwtWHvCEI/HvOHuIx27M8e/QtxkeMYf/CK9FGjkOj0fCbubkD3p7++Gh/HRu22yh9fz/funGOK61C+uXTeOqBa7k9ZWgUOR/O/FrktWnTJkpKSlyBvqFBnX+r1WqJjY0lJSWFlStXsmLFiuC1PEhSUlKIj4+XxG7DjaKcn9P/XTU3j0cadbjnni/UKZ8D5IuzB3nM9iwbj76NA3VM//aJV/N/l32PWWOmDFg7+uvU6XOUvr+fDdut/POg3bU/NTGOt35+c+gaNoJ1TuxWXV3t9lrAVvgOF7LCdxhqqoXt/w628wWG4q6Ei++BD9Z0P/amEpizckCadaD5OD/e+xc2Hnsb5fzzh2WTruVnifeSMuGSAWlDoDz0lw959t0aWtvVD6+oiDDuSLmI1YsSWTJ/CuFhYSFu4cgW1BW+QgxKjTVgugpaTqrVtJJ/Aik/hrBIOP4h7H9BnebpHO4ZoMAP0K50UHJ8GwoKyyddy88S7yN5wsUD9v496XA4eG/PCY7Zm5mqHcN1cye5BXDb8dMkTB7nGquPDA+jtd3BgtmxrF6USNY1s9CNGxWq5gsf+B38H3jgASoqKti7d28g2yNEYI2fDZOSoek4pP/NvY7u4vVw+C11hW/UBFj0ZFCb8tmZA1TU7eKhmXcCkBg9nf+d++9cp51H0oQ5QX3vvthaeQjDM9UcrruQJTNeF80vVl5JW4fChu1W3ttzgtd/ksF151MrfP/Wy/jXJXP42qzhVeR8OOtXz7+nVbyNjY1eC70IEVTWzXBRuroyV6OBzGchcnz3RGzRk2FJ4YXcPkEa5//kdA2/sj1L2fHtAKTrkpg3Tp3D/t2Z3fOsh9LWykOsfuKdbpNgD9c18Z3177u2wzQaLDV1ruA/e9K4AWylCAS/g39KSgqZmZleXzcYDDz5ZHB7UkK4aapVH+haTTAvB5YUqft7ysczZ2XQhnr+ebqGX9mecQV9gBWTbyBqkBZz73A4MDxT7Wn1g0t4mIYf330F992gZ7ouesDaJgLP79/CtLQ01q1bB0BiYiI6nc5VutFut2M2mwPSQCF8sq8Utn/3wth+9DR1hk8I5o8fPXeK73/xJ0zH33Hty5qykJ/p7+OK8QkD3h5fPbdjv9tQjycdDoVrL54kgX8Y8Dv4L126FLvd7rFWr91ul0UbYmB07u0DxH0Nlv4NJiWFrEnjw6N5q+4jNGjInrKIn+q/yeWDMOgrisKnh+xsrTrE85WH3FIm9+SYvTnILRMDwe/gr9frqaio8Lp4a+XKgZs1IUaoozvg1Ts7zeT5MaT8ZMCLrOxutLLhqJnfXJKrJiWLGMNf5j9MYvQ05o+bPaBt8cVnX9kpfX8/z+88hPX4adf+8DANHY6eBn1UUyUFw7Dgd/A3Go09rtrNy8vz99JC+Cbm/LTIEPX2dzXu4xfWDWw9oT4IXRR7BcsnXwfAssnXDmhbeuJwKCgorqmamz44wG9f/AyAUZFhZFwxnTvTZpD5telc/9NXOVLX5HHcXwNM10Vz3dxJA9d4ETR+B//09HSP+4uLi2loaCArK8vvRgnh1dEdME3NeUP0ZFj+JmjnDmhv39K4l19Yn+aF80Ffg4ZvTF3CpWNnDFgbetPhcPDB3pNsrTzIC1Vf8ZvVqdyechEAd101ky+PNnJn2kxuunK6K4c+gPHeFFY/8U63fKeaTq/Lgq3hIeDTDnJycgBYu3at64GwEP3WfEJdpWs1wU2lF0ouxl0xYE2wt53h/k8KePHEB4Aa9O+ZuoT/1N/LZeNmDlg7vGnvcLBjTy3P7zzEi9WHON7Q4nrtZctXruB/+cxYNjy00OM1lqfNYMNDC7vN85+ui8Z4bwrL0wbPB5zon34F/82bN1NeXk5dXZ3bfm8PgoXwy74yNfC3nFRX4p45GJJmxESM5VDLCcII455pS/hP/Te5dGzogz7AydMtpK55mVOnz7n2xURHcmvSRdx51QyWzp/m87WWp83g9pT4Hlf4iqHP7+C/adMmcnJyyMjIcAV7nU5HXV0dMTExlJSUBLKdYiRqPqFO37SWqdu6K9RVupOS+3XZg821nGzzPrNlYmQMM8dMZmfDF/z+wGaK5/2QcRFj0Gg0FM77ATERY5kbwiGec20dvPnJMQ6ePENe5ly1zeNHM1U7BkWB21Mu4s60GSyeN4WoiHC/3iM8LIyFlw2dpHKi7/wO/maz2dXj37VrF/X19a5c/w0NDWzevHlI5u8/fPgwy5Ytk6yeobb/ZXjzXy/09pPXQupP+z22f7C5lrk7vkWLo83rMVGaCK7TzuPt+o8BWDA+EUPCKgCuirm0X+/vr+bWdso/PsoLVYd4dddhGpvbGBMVzn2LEl2pkUt/uJjpsWOICJceulB1zurZld/BPzn5Qu9Lr9eTm5vrCv4xMTGcOuUtbe7gFh8fL1k9B4OwSDXwB6i373SyraHHwA/QqrTzdv3HhGvCuG9aOndPviEg7+2PbZ8d4y9v7uP1j45w9ly7a/+02DEsT51B07l2V/CfOXFsqJopBilnJ3bZsu5pRPwO/p0XccXExGC1Wjlw4ACzZqk5S5w5/oXw2emDMP78GPrMm+Drz8PMrw/4vH2AOyZew+8uzWNOdPyAvm9DUyuR4WFEnw/ou2rq2LxTfcYxIy6a5WkzuTNtBmmJEwkLk4WUwn9+B39FUVi5ciW7du1i7969rFmzhuTkZAoKClAUhcrKykC2UwxnzSdg+/fg4KvwjU8ufAAkLA9Zk/5rzuoBC/x1Z87xyq7DPL/zIG99eownvn0V37xBD8CdV82k7kwrd6bNIClBJyvnRcD4FPz379/frfB6Tk4OOp2OVavUsdCsrCzq6urIyclBo9F0qxojhEfWTbD9QfUDQBMOR7bB3NVBe7uWjtagXbsvTjS28FL1V2ytPMi2z4/T3nFhVn3lvlOu4D970jh+uWpBiFophjOfgn9eXh6vv/56t/1dyzTm5uaSmzu0ao2KEGk+eX4mT6m6rbsclv4VJqcG5e2On6vnT4de5A8Hnw/K9fui7sw5LvnBFreAP3+GljvTZnBn2kwujZd61yL4fAr+5eXl/Pa3v+Xhhx8OdnvESGDdDNsfuNDbT15zfiZP4Cs/fXK6ht8f3MzTR96kVen5QW8wfHXqLC9UHeKruib++x71obVu3CjSEifS3NrO8rSZLE+dwcXTpPaFGFg+BX+tVktMTAxr1qxBo9GwatUqFixYEOSmiWHr+Adq4NfNV3PyBKm33+po48aqfNec/qtjLuWuSdezZt+fg/J+TvtPnGFr5SGerzxIlVWd9RYepuHh2+cTN179gHv+0RtdD3WFCAWffvuMRiPf+c53XNvFxcWsX7+eOXPmkJubKxW7RO/azkDk+WpPV/0SoqfAFd8LaG+/ueMcW2p38I2pSwjThBEVFsl3Z9zBp2cP8KNZK7hWO4+DzbX8l+0fPU73HB0WycTIvg+9bP7wAL9/+TN276937dNo4NpLJrE8dQYR4Rce1krgF6GmURSlWwK/8PBwOjo6ej25pqYGk8mEzWYjOzvbNc+/PwoKCtyKwuTn5/f7HJvNhtlspqysjPLy8h6v5anKveiH5pPwzkNwugbu2gFh/q047YlzPP9Ph17kZFsDryQ9xtcnXQWos9K6zpDxdYVvb7443MDkmNGuQuV/eXMvP/hbJWEaDQsvm8zy1BnckTpDUiCLkPMU1/rV/UhISODRRx8FoKKigjVr1jBx4kSysrK6zQ7yhTOIOx8am81m8vLyKCws9Psci8VCVVUVdru9Ww4iEWTWzedn8tSqY/vH3oPpnhOK+eOfp2v4/YHNPHP0wnj+rNFTaHFcmNHjaWrkzDGTfQruXSmKwieH7K4hnT1HGvnt/ankZlwCwB2pM9BoNNyechGTJoz2866EGBj96vl3tWnTJtatW8euXbvIyMjwOEOoJ7GxsdTU1LglhdNoNHhoYp/PMZlMrFu3rtcpqNLzDwBnb3/fc+p2gMf269tOs+rjX1N+yuLad23MPH40627unHw9ET5+u+hwOHpNXqYoCrtq6thadYitlQexHj/jei0yPIwf3HoZP8++MiD3JUSwBLznD7B7924KCwspLS3FbrcTExPDo48+2udiLjabzWs2ULPZTEZGRkDOEQGyrxTe+QEs/MOF9MrQvbefZIC0n/V7bN+hOAjTqIFZGzGO4+fshBFG1pSF/HDW3VyjvaxP19taeahb2uJ4D2mLG5rayPhVOW0dDkAtfpL5teksT53BLQvi0Y4d+NXHQgSCX4u8GhsbKSoqorCwEJvNhqIoZGRkkJeX123uv69sNpvH/VqtFrvdHrBzRAA01cLbedBqh225MH2xWljF0QGW/1YDf4B6+8fO1fGnQy+y8dhb7LrmSVd2zeL5/8HkKC2zx0zt8zW3Vh5i9RPvdKtWdbiuifueeIdUfRxv/dfNAGjHRnFbcjwajYblqTO6FT8RYqjyKfgbDAZKSkp48803MRqNmM1mFEVBr9fz+OOPk5ub22NJx/5wpokO9jlOzqyeTpLdswtFgW0PQNv52q+tp2Hbg/D1TerD3KV/g70b+93b9zSev/HYW+RcdCvgf3bNDocDwzPVHssUOlXZTrG/9gyzJ6uzk/7xvRskrYIYUpzZPJ38zupZVlaG2Wymvl6dwpabm0teXh5JScGvmepPEO/Pg13J6tmLfaVQs+XCttIBNZvV/XNWQtzlEPdrvy7tUBy8frKK3x3YhLlul2t/5/H8/npvzwm3oR5vrMdPu4K/BH4x1HTttPYrq2dCQgJFRUV+D+v0Rq/Xe9xvt9u9vubPOaIfmmrVXn+3Cq/A27kwfYk6/OOnA83HuW3XT1FQ+jWe35WzvOHFUydwzN7s0zl1Z871fpAQQ5hPwT8jI4M33ngjqA3R6/VotVpsNlu3wO3twa0/5wg/KQqUfxNaG+gW+EFdxLX9Qbhlk8+XPHaujrfqPuKeaTcCkBA9jfunZxAXOYGHZi73azzfyRnwt+w8yNbKQ5w8fY6fZV3JNRdP9Ol8mZsvhjuPwb/rNMns7GxPhwXc2rVrMZvNrjn7JpPJLVGczWbDZDK5LeLq7RwnmePfT7t/C4crvL+udIBtM9R9qj7s7cHHp238/sBmnj36Fu1KB1fHXIo+Wq0x+7fLH/W7iR0OBzu+OMHmnQd4oeorTjReKGAee35WznVzJxGvi+ZIXZPHcX8NarHy6+ZO8rsdQgwFHuf5GwwG3nzzTdauXTvgpRg7r9a1Wq0YjUbXa0VFRRiNRqxWq8/n2Gw2CgsLMZvNWCwW8vPzSUxM9Jp9VOb5A21NsH8rRE6A2bep+5rr4G+TQXHgseevCVfz73vp+TsUB6+dH8+v6DSef512Hv936fdImjCn381ubm0n4bubXRWvYsdGcUfqDO6+aiaLLptCZIQ6VdQ524cud+Ic2d/w0EK36Z5CDHWe4prH4A9qJa5169ZRUVERkg+BUBmxwV9xqLn09/xDzbHfdhqmXA0rPrhwTON+KE3yMPSjgVFauOcLj2P+n5yuIfvjx/ji7CEAwjXnx/Nn3s3VfoznO3v4W3YeZM/RBl5Ze2GI70d/r+Rcu4O70maweN5UV8Dvytd5/kIMB31a5BUTE8Pjjz9OQ0MDRUVFpKWlkZeX55bgTQwD9V+oAf/Lp+HMoQv7JyTAjJvVufvnV8wejIym46qfk/DuD7tcRMF21c+J0MD5Glx0KB2Ea9TzZo2ZwtFzdUyIiCYn/lYemrmcWWOm9KmZztW4mz88yAvVh6htuDCk8/lXdi67SAvA7/4lzafrLU+bwe0p8b2u8BViuPLa8/fkf/7nfygtLWXVqlU88sgjwWxXyIy4nv8ry2D/i+rfo2LU6Zpz74ep16spKc872FzL3B3foqWjFdOxj1h+9gQRQBsato6dRPa0KxkdFsmWBf/Fc8fe5rMzB/nw6j+4pkm+W/8JV47XMz4ius9N3Lijhv98bpdbwI8dG8XtKRdx91Uze+zhCyECkN7h0Ucf5dFHH6W4uJjU1FQyMzNZu3atpHQeCtpb4MBLai//hv9Ve/YAl35Lnckz936YfQdEeE5IdrKtQU2DrNHw4KR5LG3egdbRzumwcP59kjp00+Jo4+uWn7jO2dnwhWtY54bYy31qprOHf1HcWBLOz7PXjYuitqHFFfDvumomSyTgC9EvfuX2ycnJIScnh02bNrF06VL5EBisFAWO7YA9G2BfyfmxemBymlo5C0B/l/qnD05ERJE36TL+9+Qevj/pUk5EXMhvE4aG7KmL+jSe7wz4W3YeZGuVOqTzw9vmuWrX3jh/KpsfWcLieVOIigh8SmghRqJ+JXZbsWIFK1asoKKigqVLl5KWlobBYPArnbPohbdEap6cs8NH/w++3ACNnfIfjb1ILY4+Z1Wf3/50u/viqLLxUykb330e/gtJv+S2SVf3ej2HQ+G9L9V5+M9XHuo2pBPZqfBJVEQ4mV+b3uc2CyG8C0g5ofT0dKqqqqioqCA3N5fY2FiMRqN8CPSTs+hIRMsp5r31HcLbTtPx9nf4bMIM2kfHuRcd6fRglrBI+Oi3F6pnJWbBJashfgloug+VdC540uZo53cHNnGwpZYDzbXqf1uO09jee0oEgGmjdD4d51AU7nviXU6dVlfSxo6N4rbk80M686WHL0SwBbSWXHp6Ounp6ezatYvc3Fw0Gg1Go1Hq/fqh6wPWeW2n1XnorY18+epysqddyXhNONYZ32DS/hfVSllZVepD2sixaqnEMVMgYTlKRDSfnz3IwVPVbkHd+fe0mEsou1IdBorQhPNL69M0OQKT3qDD4eD9L0+yZecBqm2nePNnNxMWpiEiPIz7Fuo5dfqcBHwhQiAohUSTkpJ44403qKmpwWAw0NDQgMFgCEiZx2BzZvUMdTZP5wPWlWeOs+LsCdf+CCDrbC2vHraQeq6RiXtfc732+pd/55NR4zjYUsvU8VNZqz/ffkUh9YPv0ewloOsix7v+rtFo+O7MZURpIpg1ZgozR09m5ujJnGprZGHlj3xquzPgP7/zIM9XHuR4pyGdndaTXHOxunr2sW8EPzGgECOZM7unp6yefZrq6S/nh8D+/fsxGo3ceOONwX5Lvw2WqZ6Wxr3c8m4Oew7uIMbRjrd5LUfDR/HM+KlsGD+Nj0ddCOJJ4+dgufZPru1rPvw+TR3nmDl6MrPGTHYF9VljJjNr9BTiR/ec88bSuJeUD77ba7t/O/6nPLnxpFsCtZjoSG5PmcFdV83gxvlTpYcvxAALSiWv3uzfv5+CggJMJhMA+fn5VFZWBvtthz5F4ckTnzHeS+B3ADtGxXDjRakomgimj9Jx3fmgPmv0FC4d675K9YOr/9Cv5sRGjEfTEY4S7r28p6YjnBnjdByzHyImOpLbki/i7qtnSsAXYhAKWvDvWvglKyuLtWvXDkgNgOFgdMNet+GersKAhecaMCd+m+v1K4kMC+7n+FcHNIzZehfKaO/PAjQto5j8UBybHl4iY/hCDHIBjRidyztarVa0Wm3QK30NV+fGTOWr8CjiO1rxVEqkHQ3Pj52EfnJa0AM/wDF7M2FN46BpXI/H1Ta2kH3t7KC3RwjRPwGJGs6x/KKiIhRFITk5mbKysqAVfhnumo7vZNJLt3BRRysK6hBP56EfB9B4fmXta54vEXAbtlt7PwjJgy/EUBGQ9fH5+fkUFhayYsUKqqurqaqqksDvr8//SsSW67noXD0HI0bzn7rEbj+kMOCBSZe5rawNNIdDcavrkDan5wfCGtSsmJIHX4ihISDBv7S0lPz8fEpLS2VM31/tzfDmv8Fb3ybK0c62cdO5eub1/HdsApvGTqL9/GFtaDCNnUzZ+KmMDotkYmRgh9Na2zv4xzYrqWtf5o2Pj7j2f/fmS/nt/alooNswlHPbeG+KZMUUYogI2GBxXl4eu3fvlgVdfjq+86dM+eIv6grctF9yQ3I+H56rd63wZetNKG2n0USN55Kvb6W66wrffjrd3MZf397HH1/7giP16jTNIvNebr4yHgDduFHkZlzClJgx3fLgT5c8+EIMOQEL/gkJCYG61IjS7ujgR1+u56nGjzgw9TomXfVLuCidcGDmmMlqcJ9wMdz4FLzzAyIW/oGvTb4mYO9/orGFJ9/YQ7H5S+xNbQBMix3D9265lG8t6V5dS/LgCzE8BH+aiPCso40z/3yCu9uOUl7/EYRp2JD8Q350Ubrn4+esVP8E2Mrfb6PKegqAi6dN4D9uvYxV181mVKT3aZrhYWEsvKxvxViEEIOLBP9QOPMVTa/eybgT1dwQm8D7k69gwxX53Dn5+qC/9Uf760icOp5xoyMBeDBzLk8qe/jh7fO4PfkiwsI8TSwVQgw3EvwH2qFyzr2+kuhWO/awCL6aMIv3rvp/XDE+eMNmiqKw/fPj/P6lz6j45BiPfzOZ795yKQBZ18wi+9pZrqyeQoiRQYL/QHF0QPVjKJW/YBQKllHj+e0lq/jfq3/DxKjgLIDrcDh4qforfv/yZ1Tb6gAI02g4XH/hYa309IUYmST4dxGUrJ7NJ8B8Hxx6Aw2wbdo1bLl4JX+b91DQVudu2G7ldy99xr5jpwEYHRnO/Yv1PPT1y5g9qedVukKI4aGnrJ4S/LuIj48PeFbPY/WfMuXIO2gixsCi9dww914Wa4Kb9+b1j46w79hpYsdGkZN+MQ/cNJdJEzzX5xVCDE/OTuyyZcu6vSbBPxB6KLG4o/5T7v7yTzyU+HXWLvgp4RMX0New76xx621q5XF7M396Yw//sjgR/RQ1rfOjd8zn6jkT+dclcxg/JrK/dyiEGGYk+PdXUy28nQetdtiWC9MXQ/go2PYAL05NY8VxM21KO6bxevImzKavyQ+2Vh7qtqgq/vyiqstnavnDK5/zzLs2zrU5aGxq4/f/mgbAlbN1XDnbt5KKQoiRR4K/nw4213Ky1U7C9u+hPV9iUWk9zenX7iTqzBFGnznAvP1bUWZdx4opi/j75Y8yNqJvSc+2Vh5i9RPv0LXazuG6Ju574h31Pc/vu2rORG5eIEXOhRC+keDvB2d93WUNhyg5/k/Xfo3SwYRj7wNwIGI035x6BT+d86/8p/6bhHkonN6TDocDwzPV3QJ/ZwqQ+bVpPHzHfK67ZJJM1xRC+EyCvx9OtjUwvvUs60983i3dMqjJ126ankTO/B/ySEK2p0v06r09J9yGerz54W3zuH5uYPL7CCFGDgn+/uilxKIGhV+fsqKPW+DX5c+1dfBS9Vc+Hdu5Vq4QQvhKgr8feiuxGAFkna3lM/teNSmbDxRF4aMD9Ty93Ubp+/upP9vq03lSPEUI4Q8J/n5oibmYTWMnsfzsCY//A10lFrW+Bf76s63c+t9mPjlkd+2bph3N6ZZ2zra0exz316CmUpbiKUIIf0geXn9oNDw4aR6nwyJwdHmpc4lFb9raHXy0v861rY1W5+GPigxjxdUz2fLIEj7/f3eyPuda9e26vv35/0rxFCGEv6Tn76cTEVE8MOkyt9k+0HOJxc++srNhu42S9/bTdK6dfU/cxbjRkWg0GorzriU+biyxYy+ctzxtBhseWijFU4QQASfBvx9Kx01h5ZljruGfNjRsHTuJsvFTXcfUnTmH6YMDPL3dxq5Ovf1JE0az50gjKfo4AC6fGevxPaR4ihAiGCT4+2FiZAyRRNKmaePBSfNY2rwDraOd052GeyKJpPqzM2QWbaG1XR0cigjX8PUF8dy3SE/mFdOJjPAtgEvxFCFEoEnw78KXrJ7xoyYy481vcLTFzhngoXEJ/GbyJh6uzeLMF0lEA1NHa8n80cU83PE5V8zUct9CPdnXzpbkakKIAdNTVk+Noig9LSIdcZYtW9ZrVs93Pj/Oresqer3WK2vTma6LJvF8sjUhhAgFT3FNBo794OvCqmP2Zgn8QohBSYK/H3xdWCULsIQQg5UEfz9cN3cS8brobvPvnTSoaZdlAZYQYrCS4O+H8LAwjPemALIASwgxNEl08pNzAdZ0XbTb/um6aDY8tFAWYAkhBjWZ6tkPsgBLCDFUSfDvJ1mAJYQYiqSL6oeNGzeGuglBMVzvC+Tehqrhem+D4b4k+PthMPzggmG43hfIvQ1Vw/XeBsN9SfAPMl9/yH35ZQjWsX0xGNor99Z3g6G9w/Xehtp9SfAPMgn+g+PYvhgM7ZV767tQt3eo3Zfk9uli/vz5JCYm9njM4cOHiY+P9+l6vh4bjGvKsXKsHDtwx4b6/Xs61mq18umnn7rtk+AvhBAjkAz7CCHECCTBXwghRiAJ/kIIMQKN+BW+BQUFaLVaAOx2O/n5+f0+x2azYTabKSsro7y8PNBN9lkw7q2goABQHyABFBYWBq7BfRDoe7Pb7ZSWlgLqvdlsNoqLi13HD5Rg/Mw6y8zMDNnvZKDvzWw2U1hYSGZmJnq9nvLyctLS0sjKygpG8/1uZ3/OMRgMrgkoOp0usPemjGBGo1EpLCx0bZeXlyu5ubn9Oqe6ulopLCxUjEajkpycHPhG+ygY95afn+92fG5urpKRkRGgFvsuGPeWm5urVFdXu20P9L0F4746KysrU0L1Tz4Y91ZWVqbo9XoFUPR6vduxAykY91ZfX68kJycr9fX1iqKocSXQP7sRHfy1Wq3rf65Tb/+DfT2nrKwspME/0PdWX1+vZGRkuL3u/IW0Wq2BaLLPgvFzy8jIUIxGo2vbaDQqWq22323ti2D+PtbX1yuFhYUhC/7BuLeysrJur4dCMO4tPz/f7fdRUdQPiEAascHfarV6/AEBXv8n9+WcUAb/YNxbfX29otVq3XrH9fX1CuC2L9iC/XNzysrKUrKysvrX2D4I9n0VFha6fl4DLVj3NhiCf7Duzfnvymq1BjzoO43YMX+bzeZxv1arxW63B+ycUAjGvWm1Wurr691eM5vNAOj1ev8b20cD8XMzmUzY7XbKysr8bWafBfO+zGYzGRkZ/W2i34J5b6Wlpeh0OgAqKysxGo39amtfBePenK/bbDb0ej16vZ68vDyys7MD+nMcscHfG51OR11dXdDPCYVA39u6desoLCwc8IeingTq3oqKirDb7SQnJweyeX4LxH3Z7Xb0ev2g6qBA/+9Nr9eTnJzs6nzU1dWRnZ09oB/a3vTn3pzBX6vVun4PjUYjCQkJ3Tpg/SFTPbvwJ4gPhcAPgb03g8FAXl4eubm5/W1WQATq3nJzc8nPzycxMZGEhISQB8z+3ldRUVFIZr/4or/31jnwA6xcudL1rS3U+nNvzm8yqamprtec3wqc37YDYcQGf29DFc5eUqDOCYVg35vJZCIxMTEkgT8Y92a32zEYDG5BIyMjI+D/2HoSjPuyWCxuASRUgvX7aDKZ3F5zfgP1NqwSDMG4t87fZDrTarWBvbegPEkYIrRabbeZKr39L/H1nMEw2ycY91ZeXq6UlZW5tuvr60My2yeQ9+Zp1pJzX7AetvW1jf6cU15erhiNRtef3NxcBVCMRqPbz3AgBPrenA+vO7/u3DfQD4GD8W9Nr9d3+90jwJMrRmzPH2Dt2rVuPTuTyeTWm7XZbK5FTb6e4xTqoaBg3JvFYsFisZCcnIzNZsNisbBu3TrX19SBEuh7S05OJj8/361HVlJSQnJy8oA+KA30fWVkZJCfn+/6k5eXB0B+fv6ADwUF+t60Wi2FhYVuPzPnENdAP4MKxr81o9HothjPZDKRkZER0GdRIz6rZ+dVdlar1W22QFFREUaj0bWa1ZdzbDYbhYWFmM1mLBaLa/w4FEMkgbw3u93udQw8FL9Cgf652e12ioqKXNvO10O5wjcQ9+VkMpkoKSlxBZlAzxzxRbB/ZqdOnRrw2T5Owfi5OScfQHDubcQHfyGEGIlG9LCPEEKMVBL8hRBiBJLgL4QQI5AEfyGEGIEk+AshxAgkwV8IIUYgCf5CCDECSfAXQogRSIK/EF0YDAY0Go380WhC/aMQQST5/IXoxJnJUxa+i+FOev5CdLJu3TrWrl0b6mYIEXQS/IU4z9nrD1TGS5vNhsFgCMi1hAg0Cf5CnFdUVOS1128ymcjOzu7T9QoLC8nMzOz2HtnZ2Wg0GlJSUrqdY7fbyc7OJjY2lsTERFcaZiECTbJ6CnFeSkoK1dXVrm1nhS9Qi6Dr9Xq3HOt9vZ6Tsw6CyWSiurraY4727OxsiouLB0V9ZDE8Sc9fCNTc6l172c6CIYWFhX0uotFTCUWz2UxxcTGgfjvwJC0tTQK/CCoJ/kKgVu4KZMGdwsLCHodstFotubm5bsVInOx2uwR+EXQS/MWIV1RUxKpVqwJ6zaqqql6/LTg/HLp+AJjN5gGvsiVGHgn+YsQrLCwkPz8/YNfrKXg7ayCDWjtYr9d3G/qx2WxutWmFCAYJ/mJEcxbGDqSehnyqqqrcngXk5eVhsViw2WwBbYMQvZHgL0a0YCzqstvtXnvuXcfznc8ZnL1/X8b7zWYziYmJHp8XCOErCf5iWMrOzu61N+0cngnkw9W+rgfQarVkZWW5AnnXbwaeZGRk9Hshmt1u79f5YuiT4C+GFeeqWpPJhNFo7PFYg8EQ8F5/YWEhK1eu9No2T98I8vLysNvtmEwmt2cCPYmLi/O7jTabjdLSUr/PF8ODBH8xrNjtdoxGo6s37a2HazabSU1NDWiv3/le3q7p7UGw89uHtzn/gdbbh6IYGSSrpxhWnL3mtWvXYjKZKCoq8jiTx2g0BjzYlpaW9ji3v6fx/NzcXAoKCrqlg7Db7RQVFZGcnIzdbqeystIVvJ3fFkBdp1BWVuY6z2KxuFYl22w2srKy0Ov1mM1mqqqqqKurA9QPHp1O5/U9xDCmCDFMJScnK3q9vtv+6upqJSsrq0/XysrKUjIyMno8pqfX6+vre3xPq9WqAEp1dbXbfqPRqJSXl7u2CwsLXfs7Xy8jI8N1rtVq7daW5ORkpb6+XlEURcnPz3ddp6f3EMObDPuIYWvt2rXYbDbMZrPbfoPB4HPPNi8vj+zsbEwmE2azmczMTNcYfWc2m81rrz47O5uEhATMZrPXh8F6vZ7c3Nxu4/1ZWVlkZ2eTkpJCQUGB2/OEtLQ019+1Wq2rN+8pHYVer/c6zt/Te4jhS4Z9xLCVlZWFVqvFaDS6xtqdM4B8XUTl69CQyWTyOuTTeTimr++l0+mor6/HYrFQUlJCdnZ2n5LL9cZutwf9PcTgJD1/MaytXbsWs9nsCvp96fX3RUlJSVBSMqxbtw6bzUZycjJGo9GnB9SrVq3q9m3HYrF47NGbzWa/3kMMfZLSWQxrdrud2NhYcnNzMRgM5OXlBbxXa7FYXNk/A62goACtVotOp6Ourg6dToderycnJweA4uJi1/RWZ/B2Pti1WCzo9XoqKytZtWqVayjIZrNhNBpJSUkhIyMDk8nU7T0CVdBGDF4S/MWw5xyzz8rKIi8vL+A9dIPB4BZchRgKZNhHDHvOhVw2my0oQzPOIRMhhhLp+YsRISUlxe3BrxAjnQR/IYQYgWTYRwghRiAJ/kIIMQJJ8BdCiBFIgr8QQoxA/x+74CSHd/w1/wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 400x300 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(np.array(shots_list) ** (-0.5), e_array.std(axis=2)[0], marker=\"o\", label=\"$p=0.10$\", linestyle=\"--\")\n",
    "plt.plot(np.array(shots_list) ** (-0.5), e_array.std(axis=2)[1], marker=\"s\", label=\"$p=0.20$\", linestyle=\"--\")\n",
    "plt.plot(np.array(shots_list) ** (-0.5), e_array.std(axis=2)[2], marker=\"d\", label=\"$p=0.40$\", linestyle=\"--\")\n",
    "plt.xlabel(r\"$\\sqrt{1 / N_{\\rm{shots}}}$\")\n",
    "plt.ylabel(r\"$\\sqrt{\\textrm{Var}(E)}$ (Hartree)\")\n",
    "plt.legend()\n",
    "plt.savefig(\"uncertainty.pdf\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f8279834-6f87-456c-a1db-7caf482b237a",
   "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.7.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}