{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "40e4ae41-82bd-4c33-8855-d3c596d29d49",
   "metadata": {},
   "source": [
    "# Example 4 (Advanced options): PEC of H$_{2}$ using Qiskit IBM Runtime Service.\n",
    "\n",
    "In this section, we demonstrate how to use `qc2` alongside `Qiskit` and the `IBM Quantum Runtime Service` for calculating the potential energy curve of the hydrogen molecule. We utilize the `EstimatorV2` primitive as implemented in the `qiskit-ibm-runtime` package, in conjunction with the `ibmq_qasm_simulator`, to simulate a real quantum backend."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b1f5d21a-fbab-40f6-98bc-e7cc625551de",
   "metadata": {},
   "source": [
    "### Import required packages"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "9a6303c0-283f-42e9-9c87-8b21eb994f08",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "# ASE Atoms object\n",
    "from ase import Atoms\n",
    "\n",
    "# IBMQ runtime service packages\n",
    "from qiskit_ibm_runtime import (\n",
    "    QiskitRuntimeService,\n",
    "    EstimatorV2 as Estimator,\n",
    "    EstimatorOptions,\n",
    "    Session\n",
    ")\n",
    "\n",
    "# Qiskit-related packages\n",
    "from qiskit_nature.second_q.circuit.library import HartreeFock, UCCSD\n",
    "from qiskit_nature.second_q.mappers import JordanWignerMapper\n",
    "from qiskit_algorithms.optimizers import COBYLA\n",
    "from qiskit.transpiler.preset_passmanagers import generate_preset_pass_manager\n",
    "\n",
    "# ignore package import warnings\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "\n",
    "# qc2 packages\n",
    "from qc2.data import qc2Data\n",
    "from qc2.ase import PySCF"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3d0474e0-ef1a-4438-a1a9-7d1406a63551",
   "metadata": {},
   "source": [
    "### Initiate IBM runtime service\n",
    "\n",
    "In the following steps, we assume that you have already created an __[IBM Quantum account](https://quantum.ibm.com/)__ and saved your accound credentials (via your personal token) on disk; for further info click __[here](https://github.com/Qiskit/qiskit-ibm-runtime/blob/main/README.md)__."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "6fc4456c-96d2-4d89-ad06-5452fe58ad43",
   "metadata": {},
   "outputs": [],
   "source": [
    "# instantiating `QiskitRuntimeService`\n",
    "service = QiskitRuntimeService()\n",
    "\n",
    "# run on simulator\n",
    "backend = service.backend(\"ibmq_qasm_simulator\")\n",
    "# Use the following line if you want to run on a real quantum backend\n",
    "# backend = service.least_busy(simulator=False)\n",
    "\n",
    "options = EstimatorOptions()\n",
    "\n",
    "# set seed for reproducible runs\n",
    "options.seed_estimator = 42\n",
    "\n",
    "# error supression options\n",
    "options.optimization_level = 1\n",
    "\n",
    "# error mitigation options\n",
    "options.resilience_level = 0  # 2 = ZNE\n",
    "\n",
    "# execution options\n",
    "options.default_shots = 6000\n",
    "\n",
    "# instantiate the Estimator to calculate expectation values\n",
    "estimator = Estimator(backend=backend, options=options)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6c57a140-6a65-4b95-840c-5ea8c27fa8a0",
   "metadata": {},
   "source": [
    "### Set up calculation details"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "c2bbf6f2-9472-459a-983f-698e9bc6e85e",
   "metadata": {},
   "outputs": [],
   "source": [
    "# define activate space\n",
    "n_active_electrons = (1, 1)\n",
    "n_active_spatial_orbitals = 2\n",
    "\n",
    "# define the type of fermionic-to-qubit transformation\n",
    "mapper = JordanWignerMapper()\n",
    "\n",
    "# define file to write the calculated data to\n",
    "output_file = 'h2_pec_ibmq_qasm_simulator.data'\n",
    "file = open(output_file, \"w\", encoding='utf-8')\n",
    "\n",
    "# set up equilibrium geometry and initial circuit parameters in rads\n",
    "req = 0.737166\n",
    "req_params = [1.57079357, 1.57087253, 1.45852109]\n",
    "\n",
    "# define vector of internuclear distances to loop over\n",
    "r_distances = [r*req for r in np.arange(0.5, 5, 0.5)]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c4e640c6-8c67-46ad-9e24-103b96969ee1",
   "metadata": {},
   "source": [
    "### Run VQE for each selected geometry"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "e263182f-a351-4816-b30a-2af628d77159",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "* Reference energy (Hartree): -0.8323959925858835\n",
      "* Saving qchem data in h2.hdf5\n",
      "\n",
      ">>> Final ground state energy (hartree): -0.8366257350658743\n",
      "\n",
      "* Reference energy (Hartree): -1.11690055771897\n",
      "* Saving qchem data in h2.hdf5\n",
      "\n",
      ">>> Final ground state energy (hartree): -1.1349489643349955\n",
      "\n",
      "* Reference energy (Hartree): -1.0347677703871998\n",
      "* Saving qchem data in h2.hdf5\n",
      "\n",
      ">>> Final ground state energy (hartree): -1.0739508097501043\n",
      "\n",
      "* Reference energy (Hartree): -0.9186044326889025\n",
      "* Saving qchem data in h2.hdf5\n",
      "\n",
      ">>> Final ground state energy (hartree): -1.0002826958494946\n",
      "\n",
      "* Reference energy (Hartree): -0.8185133354304801\n",
      "* Saving qchem data in h2.hdf5\n",
      "\n",
      ">>> Final ground state energy (hartree): -0.9588193874386988\n",
      "\n",
      "* Reference energy (Hartree): -0.7444704693706197\n",
      "* Saving qchem data in h2.hdf5\n",
      "\n",
      ">>> Final ground state energy (hartree): -0.939401703419809\n",
      "\n",
      "* Reference energy (Hartree): -0.6935395490285441\n",
      "* Saving qchem data in h2.hdf5\n",
      "\n",
      ">>> Final ground state energy (hartree): -0.9350495720648158\n",
      "\n",
      "* Reference energy (Hartree): -0.6597173498273768\n",
      "* Saving qchem data in h2.hdf5\n",
      "\n",
      ">>> Final ground state energy (hartree): -0.9335246431974368\n",
      "\n",
      "* Reference energy (Hartree): -0.6376873315318315\n",
      "* Saving qchem data in h2.hdf5\n",
      "\n",
      ">>> Final ground state energy (hartree): -0.9329587957364209\n",
      "\n"
     ]
    }
   ],
   "source": [
    "for r in r_distances:\n",
    "\n",
    "    # set Atoms object\n",
    "    mol = Atoms(\n",
    "        'H2',\n",
    "        positions=[[0.,  0.,  0.],\n",
    "                   [0.,  0.,  r]]\n",
    "    )\n",
    "\n",
    "    # instantiate qc2Data class\n",
    "    qc2data = qc2Data(\n",
    "        molecule=mol,\n",
    "        filename='h2.hdf5',\n",
    "        schema='qcschema'\n",
    "    )\n",
    "\n",
    "    # specify the qchem calculator\n",
    "    qc2data.molecule.calc = PySCF(method='scf.RHF', basis='sto-3g')\n",
    "\n",
    "    # run calculation and save qchem data in the hdf5 file\n",
    "    qc2data.run()\n",
    "\n",
    "    # set up qubit Hamiltonian and core energy based on given activate space\n",
    "    e_core, qubit_op = qc2data.get_qubit_hamiltonian(\n",
    "        n_active_electrons,\n",
    "        n_active_spatial_orbitals,\n",
    "        mapper,\n",
    "        format='qiskit'\n",
    "    )\n",
    "\n",
    "    # build up the initial quantum circuit\n",
    "    reference_state = HartreeFock(\n",
    "        n_active_spatial_orbitals,\n",
    "        n_active_electrons,\n",
    "        mapper,\n",
    "    )\n",
    "\n",
    "    # set up the ansatz using unitary CC as variational form\n",
    "    ansatz = UCCSD(\n",
    "        n_active_spatial_orbitals,\n",
    "        n_active_electrons,\n",
    "        mapper,\n",
    "        initial_state=reference_state,\n",
    "    )\n",
    "\n",
    "    # create initial callback dictionary\n",
    "    callback_dict = {\n",
    "            \"backend\": None,\n",
    "            \"iters\": 0,\n",
    "            \"cost_history\": [],\n",
    "            \"parameters\": [req_params],\n",
    "            \"gs_energy_history\": [],\n",
    "            \"variance\": []\n",
    "    }\n",
    "\n",
    "    def cost_func(params: list) -> float:\n",
    "        \"\"\"Return electronic energy from Estimator for given set of parameters.\n",
    "\n",
    "        Args:\n",
    "            params: list\n",
    "                List of VQE parameters.\n",
    "\n",
    "        Returns:\n",
    "            energy: float\n",
    "                Calculated electronic ground-state energy.\n",
    "        \"\"\"\n",
    "        # perform transpilation\n",
    "        pm = generate_preset_pass_manager(backend=backend, optimization_level=0)\n",
    "        isa_circuit = pm.run(ansatz)\n",
    "        isa_observable = qubit_op.apply_layout(isa_circuit.layout)\n",
    "\n",
    "        # run job\n",
    "        job = estimator.run(\n",
    "            [(isa_circuit, isa_observable, params)]\n",
    "        )\n",
    "        result = job.result()\n",
    "        energy = result[0].data.evs\n",
    "        variance = result[0].data.stds\n",
    "        parameters = list(params)\n",
    "\n",
    "        callback_dict[\"iters\"] += 1\n",
    "        callback_dict[\"cost_history\"].append(energy)\n",
    "        callback_dict[\"backend\"] = job.backend()\n",
    "        callback_dict[\"gs_energy_history\"].append(energy+e_core)\n",
    "        callback_dict[\"variance\"].append(variance)\n",
    "        callback_dict[\"parameters\"].append(parameters)\n",
    "\n",
    "        # print intermediate data if needed\n",
    "        # print(\"=== COST FUNCTION SUMMARY ===\")\n",
    "        # print(f\">>> Job ID: {job.job_id()}\")\n",
    "        # print(f\">>> Job Status: {job.status()}\")\n",
    "        # print(\"=============================\")\n",
    "        # print(f\">>> Job Input: {job.inputs}\")\n",
    "        # print(\"=============================\")\n",
    "        # print(f\">>> Backend: {job.backend()}\")\n",
    "        # print(f\">>> {result}\")\n",
    "        # print(\"=============================\")\n",
    "        # print(f\">>> Expectation value (Hartree): {energy}\")\n",
    "        # print(f\">>> Total ground state energy (Hartree): {energy+e_core}\")\n",
    "        # print(\"=============================\\n\")\n",
    "        return energy\n",
    "\n",
    "    # initiate Qiskit Runtime Section\n",
    "    with Session(service=service, backend=backend) as session:\n",
    "\n",
    "        # initial circuit parameters\n",
    "        initial_theta = np.array(callback_dict['parameters'][-1])\n",
    "\n",
    "        # set up the classical optimizer\n",
    "        optimizer = COBYLA(maxiter=300)\n",
    "        \n",
    "        # find the best circuit parameters that minimizes the electronic energy\n",
    "        res = optimizer.minimize(\n",
    "            cost_func,\n",
    "            x0=initial_theta\n",
    "        )\n",
    "        print(f\">>> Final ground state energy (hartree): {res.fun+e_core}\\n\")\n",
    "\n",
    "        # write relevant data to file\n",
    "        data_to_write = f\"{r} {callback_dict['gs_energy_history'][-1]} \" \\\n",
    "                        f\"{callback_dict['variance'][-1]}\\n\"\n",
    "        file.write(data_to_write)\n",
    "\n",
    "        # use the optimized circuit parameters as guess for the next geometry\n",
    "        req_params = callback_dict['parameters'][-1]\n",
    "\n",
    "        session.close()\n",
    "\n",
    "file.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5d3c7be0-962e-4417-9bc5-de161fd015ce",
   "metadata": {},
   "source": [
    "### Calculate PEC using PySCF for comparison\n",
    "\n",
    "With our VQE points done, let's calculate the PEC using a classical qchem backend to compare with, but now using a finer grid of internuclear distances."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "7bdfa00e-420a-4864-be4a-ad8b089a711d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "* Energy (hartree) -0.841556 :: r (angstrom) 0.369\n",
      "* Energy (hartree) -0.924999 :: r (angstrom) 0.405\n",
      "* Energy (hartree) -0.987489 :: r (angstrom) 0.442\n",
      "* Energy (hartree) -1.034243 :: r (angstrom) 0.479\n",
      "* Energy (hartree) -1.068956 :: r (angstrom) 0.516\n",
      "* Energy (hartree) -1.094322 :: r (angstrom) 0.553\n",
      "* Energy (hartree) -1.112352 :: r (angstrom) 0.590\n",
      "* Energy (hartree) -1.124586 :: r (angstrom) 0.627\n",
      "* Energy (hartree) -1.132223 :: r (angstrom) 0.663\n",
      "* Energy (hartree) -1.136210 :: r (angstrom) 0.700\n",
      "* Energy (hartree) -1.137302 :: r (angstrom) 0.737\n",
      "* Energy (hartree) -1.136103 :: r (angstrom) 0.774\n",
      "* Energy (hartree) -1.133097 :: r (angstrom) 0.811\n",
      "* Energy (hartree) -1.128672 :: r (angstrom) 0.848\n",
      "* Energy (hartree) -1.123139 :: r (angstrom) 0.885\n",
      "* Energy (hartree) -1.116748 :: r (angstrom) 0.921\n",
      "* Energy (hartree) -1.109701 :: r (angstrom) 0.958\n",
      "* Energy (hartree) -1.102166 :: r (angstrom) 0.995\n",
      "* Energy (hartree) -1.094279 :: r (angstrom) 1.032\n",
      "* Energy (hartree) -1.086157 :: r (angstrom) 1.069\n",
      "* Energy (hartree) -1.077900 :: r (angstrom) 1.106\n",
      "* Energy (hartree) -1.069594 :: r (angstrom) 1.143\n",
      "* Energy (hartree) -1.061316 :: r (angstrom) 1.179\n",
      "* Energy (hartree) -1.053131 :: r (angstrom) 1.216\n",
      "* Energy (hartree) -1.045097 :: r (angstrom) 1.253\n",
      "* Energy (hartree) -1.037263 :: r (angstrom) 1.290\n",
      "* Energy (hartree) -1.029672 :: r (angstrom) 1.327\n",
      "* Energy (hartree) -1.022360 :: r (angstrom) 1.364\n",
      "* Energy (hartree) -1.015354 :: r (angstrom) 1.401\n",
      "* Energy (hartree) -1.008678 :: r (angstrom) 1.437\n",
      "* Energy (hartree) -1.002347 :: r (angstrom) 1.474\n",
      "* Energy (hartree) -0.996374 :: r (angstrom) 1.511\n",
      "* Energy (hartree) -0.990763 :: r (angstrom) 1.548\n",
      "* Energy (hartree) -0.985517 :: r (angstrom) 1.585\n",
      "* Energy (hartree) -0.980631 :: r (angstrom) 1.622\n",
      "* Energy (hartree) -0.976101 :: r (angstrom) 1.659\n",
      "* Energy (hartree) -0.971916 :: r (angstrom) 1.695\n",
      "* Energy (hartree) -0.968065 :: r (angstrom) 1.732\n",
      "* Energy (hartree) -0.964534 :: r (angstrom) 1.769\n",
      "* Energy (hartree) -0.961307 :: r (angstrom) 1.806\n",
      "* Energy (hartree) -0.958367 :: r (angstrom) 1.843\n",
      "* Energy (hartree) -0.955696 :: r (angstrom) 1.880\n",
      "* Energy (hartree) -0.953277 :: r (angstrom) 1.917\n",
      "* Energy (hartree) -0.951092 :: r (angstrom) 1.953\n",
      "* Energy (hartree) -0.949123 :: r (angstrom) 1.990\n",
      "* Energy (hartree) -0.947353 :: r (angstrom) 2.027\n",
      "* Energy (hartree) -0.945764 :: r (angstrom) 2.064\n",
      "* Energy (hartree) -0.944341 :: r (angstrom) 2.101\n",
      "* Energy (hartree) -0.943069 :: r (angstrom) 2.138\n",
      "* Energy (hartree) -0.941933 :: r (angstrom) 2.175\n",
      "* Energy (hartree) -0.940921 :: r (angstrom) 2.211\n",
      "* Energy (hartree) -0.940019 :: r (angstrom) 2.248\n",
      "* Energy (hartree) -0.939218 :: r (angstrom) 2.285\n",
      "* Energy (hartree) -0.938506 :: r (angstrom) 2.322\n",
      "* Energy (hartree) -0.937875 :: r (angstrom) 2.359\n",
      "* Energy (hartree) -0.937315 :: r (angstrom) 2.396\n",
      "* Energy (hartree) -0.936819 :: r (angstrom) 2.433\n",
      "* Energy (hartree) -0.936380 :: r (angstrom) 2.470\n",
      "* Energy (hartree) -0.935991 :: r (angstrom) 2.506\n",
      "* Energy (hartree) -0.935648 :: r (angstrom) 2.543\n",
      "* Energy (hartree) -0.935345 :: r (angstrom) 2.580\n",
      "* Energy (hartree) -0.935077 :: r (angstrom) 2.617\n",
      "* Energy (hartree) -0.934841 :: r (angstrom) 2.654\n",
      "* Energy (hartree) -0.934633 :: r (angstrom) 2.691\n",
      "* Energy (hartree) -0.934450 :: r (angstrom) 2.728\n",
      "* Energy (hartree) -0.934289 :: r (angstrom) 2.764\n",
      "* Energy (hartree) -0.934147 :: r (angstrom) 2.801\n",
      "* Energy (hartree) -0.934022 :: r (angstrom) 2.838\n",
      "* Energy (hartree) -0.933912 :: r (angstrom) 2.875\n",
      "* Energy (hartree) -0.933816 :: r (angstrom) 2.912\n",
      "* Energy (hartree) -0.933732 :: r (angstrom) 2.949\n",
      "* Energy (hartree) -0.933658 :: r (angstrom) 2.986\n",
      "* Energy (hartree) -0.933594 :: r (angstrom) 3.022\n",
      "* Energy (hartree) -0.933537 :: r (angstrom) 3.059\n",
      "* Energy (hartree) -0.933488 :: r (angstrom) 3.096\n",
      "* Energy (hartree) -0.933445 :: r (angstrom) 3.133\n",
      "* Energy (hartree) -0.933407 :: r (angstrom) 3.170\n",
      "* Energy (hartree) -0.933374 :: r (angstrom) 3.207\n",
      "* Energy (hartree) -0.933346 :: r (angstrom) 3.244\n",
      "* Energy (hartree) -0.933321 :: r (angstrom) 3.280\n",
      "* Energy (hartree) -0.933300 :: r (angstrom) 3.317\n",
      "* Energy (hartree) -0.933281 :: r (angstrom) 3.354\n",
      "* Energy (hartree) -0.933265 :: r (angstrom) 3.391\n",
      "* Energy (hartree) -0.933251 :: r (angstrom) 3.428\n",
      "* Energy (hartree) -0.933239 :: r (angstrom) 3.465\n",
      "* Energy (hartree) -0.933228 :: r (angstrom) 3.502\n",
      "* Energy (hartree) -0.933219 :: r (angstrom) 3.538\n",
      "* Energy (hartree) -0.933211 :: r (angstrom) 3.575\n",
      "* Energy (hartree) -0.933204 :: r (angstrom) 3.612\n",
      "* Energy (hartree) -0.933199 :: r (angstrom) 3.649\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "import pyscf\n",
    "from pyscf import gto, scf, fci\n",
    "\n",
    "# set up equilibrium geometry and vector of internuclear distances\n",
    "req = 0.737166\n",
    "r_distances = [r*req for r in np.arange(0.5, 5, 0.05)]\n",
    "\n",
    "# define file to write data to\n",
    "output_file = 'pyscf_fci_sto3g_h2.data'\n",
    "\n",
    "with open(output_file, \"w\", encoding='utf-8') as file:\n",
    "\n",
    "    for r in r_distances:\n",
    "\n",
    "        # run HF\n",
    "        mol_h2 = gto.M(atom=f'H 0 0 0; H 0 0 {r}', basis='sto-3g', verbose=0)\n",
    "        scf_h2 = scf.RHF(mol_h2)\n",
    "        hf_energy = scf_h2.kernel()\n",
    "\n",
    "        # run FCI\n",
    "        hf_h2 = mol_h2.RHF().run()\n",
    "        fci_energy = fci.FCI(hf_h2).kernel()[0]\n",
    "\n",
    "        print(f'* Energy (hartree) {fci_energy:.6f} :: r (angstrom) {r:.3f}')\n",
    "\n",
    "        # write data to file\n",
    "        data_to_write = f\"{r} {hf_energy} {fci_energy}\\n\"\n",
    "        file.write(data_to_write)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6774ab15-1f7e-47e4-a900-a0ef4b424d89",
   "metadata": {},
   "source": [
    "### Plot PECs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "1daf832b-298a-4bd9-996c-17c41af484ac",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAksAAAHJCAYAAACVLX1MAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/H5lhTAAAACXBIWXMAAA9hAAAPYQGoP6dpAAB2SklEQVR4nO3dd3gU1fs28Hs3vTdSIaQQSEIPBJCiVAnFSC/SEWmCKE1BFARFQBFRVLB9KQKCVBGRXqSEFoi0ECEQQEhCCUnIpmfP+wdv5pclm2U32c1mk/tzXXOZOXNm5plx2Tw5c+YcmRBCgIiIiIjUkhs7ACIiIqKKjMkSERERkQZMloiIiIg0YLJEREREpAGTJSIiIiINmCwRERERacBkiYiIiEgDJktEREREGjBZIiIiItKAyRIRERGRBkyWiIiIiDRgskREWlm1ahVkMhkSEhLKdV+qGs6cOYNWrVrBzs4OMpkMMTExxg6JSMJkicgACpODwsXa2hp16tTBxIkTkZycrPPxTpw4gY8++gipqan6D9aI5yICgLy8PPTr1w8pKSn48ssv8csvv8DPz09t3cJ/W2fPnlW7vV27dqhfv74hw6UqiMkSkQHNmzcPv/zyC7755hu0atUKy5cvR8uWLZGZmanTcU6cOIG5c+eWW7Kk7lxDhw5FVlZWib/EiEorPj4et27dwrRp0zBmzBgMGTIELi4uxg6LSGJu7ACIKrOuXbsiPDwcAPDGG2/Azc0NS5Yswe+//47XXnvNyNHpxszMDGZmZsYOo1wpFArY2dkZOwy9qajXc//+fQCAs7OzcQMhKgFblojKUYcOHQAAN2/elMrOnz+Prl27wtHREfb29ujYsSNOnjwpbf/oo48wffp0AEBAQID0aK+w/8/du3fx+uuvw9PTE1ZWVqhXrx7+97//qZz3o48+gkwmw/Xr1zFixAg4OzvDyckJI0eOVGnl0nQudf2Obt26hTfffBPBwcGwsbGBm5sb+vXrV6a+Sfq8ntIc88qVKxg0aBBcXFzQpk0bafvhw4cRHh4Oa2tr1KpVC99//720T6FDhw5BJpNh27ZtxWJYv349ZDIZoqKinnv9o0aNgo+PD6ysrBAQEIDx48cjNzcXADBixAj4+/sX2+/ZWEq6ns2bN0Mmk+HIkSPFjvH9999DJpPh0qVLOt07TZ73+R4xYgTatm0LAOjXrx9kMhnatWun9fGJygNblojKUXx8PADAzc0NAHD58mW8+OKLcHR0xLvvvgsLCwt8//33aNeuHY4cOYIWLVqgd+/e+Pfff/Hrr7/iyy+/RLVq1QAA7u7uSE5OxgsvvACZTIaJEyfC3d0df/31F0aNGoX09HS88847Kufv378/AgICsGDBApw7dw4//fQTPDw8sGjRIgDQeC51zpw5gxMnTmDgwIGoUaMGEhISsHz5crRr1w5XrlyBra2tTvdH39dTmmP269cPtWvXxqeffgohBICnv/C7dOkCb29vzJ07FwUFBZg3b16x+9KuXTv4+vpi3bp16NWrl8q2devWoVatWmjZsmWJ13/v3j00b94cqampGDNmDEJCQnD37l1s3rwZmZmZsLS01Ol+qrue7t27w97eHr/99puUpBTauHEj6tWrJ/X50fXePUubz/fYsWNRvXp1fPrpp5g0aRKaNWsGT0/P515XWloaHj58WKw8Ly9P+5tDpC1BRHq3cuVKAUDs379fPHjwQNy5c0ds2LBBuLm5CRsbG/Hff/8JIYTo2bOnsLS0FPHx8dK+9+7dEw4ODuKll16Syj7//HMBQNy8eVPlPKNGjRLe3t7i4cOHKuUDBw4UTk5OIjMzUwghxJw5cwQA8frrr6vU69Wrl3Bzc1MpK+lchddUtLzw+EVFRUUJAGLNmjUa91XHENej6zFfe+21YnFFRkYKW1tbcffuXans2rVrwtzcXDz7NTpz5kxhZWUlUlNTpbL79+8Lc3NzMWfOHI3XP2zYMCGXy8WZM2eKbVMqlUIIIYYPHy78/PyKbS+M/9l1ddfz2muvCQ8PD5Gfny+VJSYmCrlcLubNmyeVaXvvSqLt5/vQoUMCgNi0aZPG4wnxf58lTUu9evWeexwiXfAxHJEBderUCe7u7vD19cXAgQNhb2+Pbdu2oXr16igoKMDevXvRs2dPBAYGSvt4e3tj0KBBOHbsGNLT00s8thACW7ZsQWRkJIQQePjwobREREQgLS0N586dU9ln3LhxKusvvvgiHj16pPE8mtjY2Eg/5+Xl4dGjRwgKCoKzs3Oxcz+PIa5HH8csKCjA/v370bNnT/j4+EjlQUFB6Nq1a7HrGDZsGHJycrB582apbOPGjcjPz8eQIUNKvH6lUont27cjMjJS6udWVNFHbLp49noAYMCAAbh//z4OHz4slW3evBlKpRIDBgwAULp7V1RZP9/P8+2332Lfvn3FloYNG5b6mEQl4WM4IgP69ttvUadOHZibm8PT0xPBwcGQy5/+jfLgwQNkZmYiODi42H6hoaFQKpW4c+cO6tWrp/bYDx48QGpqKn744Qf88MMPausUdpwtVLNmTZX1wjeOHj9+DEdHR52vLysrCwsWLMDKlStx9+5d6bEV8PQxiS4McT2lOWZAQECx7VlZWQgKCiq2r7qykJAQNGvWDOvWrcOoUaMAPH0E98ILL6itX+jBgwdIT0/X+2vvz14PAHTp0gVOTk7YuHEjOnbsCOBpQte4cWPUqVNHikfXe1dUWT/fz9O8eXO1SaWLi4vK47mcnByMHz8e+/fvR2pqKurWrYsvv/xS4+NQomcxWSIyoJK+0PVBqVQCAIYMGYLhw4errfPsX9klvc1WNMnRxVtvvYWVK1finXfeQcuWLeHk5ASZTIaBAwdK8WnLENdTmmMWbS0rrWHDhuHtt9/Gf//9h5ycHJw8eRLffPNNmY8LlNzCVFBQoLZc3fVYWVmhZ8+e2LZtG7777jskJyfj+PHj+PTTT6U6pbl3FVF+fj78/f1x7Ngx1KhRA7/99hsiIyORkJAAe3t7Y4dHJoLJEpGRuLu7w9bWFnFxccW2Xb16FXK5HL6+vgDU/4J0d3eHg4MDCgoK0KlTJ73Fpcvjns2bN2P48OH44osvpLLs7OxSjQdliOvRxzE9PDxgbW2N69evF9umrgwABg4ciClTpuDXX39FVlYWLCwspMdbmmJ1dHRUeRNNHRcXF7X399atWxr3e9aAAQOwevVqHDhwALGxsRBCqMRY1nuny+fbkOzs7DB79mxpvfD/TVxcHJo2bWrw81PlwD5LREZiZmaGzp074/fff1d51T45ORnr169HmzZtpEdjhWPjFP0laWZmhj59+mDLli1qf8E+ePCgVHGpO5ema3i2VWrZsmUltnI871j6vh59HNPMzAydOnXC9u3bce/ePan8+vXr+Ouvv9TuU61aNXTt2hVr167FunXr0KVLF+nNwpLI5XL07NkTf/zxh9rRqQvvc61atZCWloYLFy5I2xITE9UOV6BJp06d4Orqio0bN2Ljxo1o3ry5yiO7st47XT7f5enatWtISUnR+EiU6FlsWSIyok8++QT79u1DmzZt8Oabb8Lc3Bzff/89cnJy8Nlnn0n1Cv8CnjVrFgYOHAgLCwtERkZi4cKFOHToEFq0aIHRo0ejbt26SElJwblz57B//36kpKToHFNJ51LnlVdewS+//AInJyfUrVsXUVFR2L9/vzQ0gq4McT36OOZHH32EvXv3onXr1hg/fjwKCgrwzTffoH79+iXOYTZs2DD07dsXAPDxxx9rFeunn36KvXv3om3bthgzZgxCQ0ORmJiITZs24dixY3B2dsbAgQPx3nvvoVevXpg0aRIyMzOxfPly1KlTR6dO9RYWFujduzc2bNgAhUKBxYsXF6tT1nun7ee7vGRlZWHIkCGYOXMmnJycyv38ZMKM8QoeUWVX+HqzulfAn3Xu3DkREREh7O3tha2trWjfvr04ceJEsXoff/yxqF69upDL5Sqv4ScnJ4sJEyYIX19fYWFhIby8vETHjh3FDz/8IO1b+Br5gwcP1Mb57Cv96s6lru7jx4/FyJEjRbVq1YS9vb2IiIgQV69eFX5+fmL48OHPPY86hrieshyz0IEDB0RYWJiwtLQUtWrVEj/99JOYOnWqsLa2Vls/JydHuLi4CCcnJ5GVlfXc6y5069YtMWzYMOHu7i6srKxEYGCgmDBhgsjJyZHq7N27V9SvX19YWlqK4OBgsXbt2hKHDijpeoQQYt++fQKAkMlk4s6dO2rraHPvNNHm812aoQNK+rfVtm1btUMH5Obmiu7du4tBgwZJwzAQaUsmRCl7dhIRVXE9e/bE5cuXce3atWLb8vPz4ePjg8jISPz8889GiI4KKZVKDBo0CAqFAtu2bYO5OR+qkG7YZ4mISAtZWVkq69euXcOuXbtKnJpj+/btePDgAYYNG1YO0ZEmY8eOlR5nMlGi0mDLEhGRFry9vTFixAgEBgbi1q1bWL58OXJycnD+/HnUrl1bqnfq1ClcuHABH3/8MapVq6bz4JykX7du3YK/vz+sra1Vhpr466+/8OKLLxoxMjIlTLGJiLTQpUsX/Prrr0hKSoKVlRVatmyJTz/9VCVRAoDly5dj7dq1aNy4MVatWmWcYEni5+dX6nHEiAqxZYmIiIhIA/ZZIiIiItKAyRIRERGRBuyzpAdKpRL37t2Dg4NDqWcGJyIiovIlhMCTJ0/g4+MjTXKuDpMlPbh37165zHFERERE+nfnzh3UqFGjxO1MlvTAwcEBwNObbYy5joiIiEh36enp8PX1lX6Pl4TJkh4UPnpzdHRkskRERGRinteFhh28iYiIiDRgskRERESkAZMlIiIiIg3YZ4mI6P8rKChAXl6escMgIj2xsLBQmROwtJgsEVGVJ4RAUlISUlNTjR0KEemZs7MzvLy8yjQOIpMlIqryChMlDw8P2NracnBZokpACIHMzEzcv38fAODt7V3qYzFZIqIqraCgQEqU3NzcjB0OEemRjY0NAOD+/fvw8PAo9SM5dvAmoiqtsI+Sra2tkSMhIkMo/Lddlv6ITJaIiPD8QemIyDTp4982kyUiIiIiDZgsEREREWnAZImIyESNGDECMpms2HL9+nUAT9/ye+uttxAYGAgrKyv4+voiMjISBw4ckI7h7++PpUuXFjt2QEAA9u/fr/H8hw8fhkwm09uQC6+++ipq1qwJa2treHt7Y+jQobh3755ejk1UFkyWiIhMWJcuXZCYmKiyBAQEICEhAU2bNsXBgwfx+eef4+LFi9i9ezfat2+PCRMmaDzmhQsX8PjxY7Rt27acruKp9u3b47fffkNcXBy2bNmC+Ph49O3bt1xjIFLHZJKllJQUDB48GI6OjnB2dsaoUaOQkZGhcZ+kpCQMHToUXl5esLOzQ5MmTbBly5YyH5eIqKKwsrKCl5eXymJmZoY333wTMpkMp0+fRp8+fVCnTh3Uq1cPU6ZMwcmTJzUe8/fff0eXLl1gYWGBW7duITIyEi4uLrCzs0O9evWwa9cuJCQkoH379gAAFxcXyGQyjBgxAgCQk5ODSZMmwcPDA9bW1mjTpg3OnDnz3GuZPHkyXnjhBfj5+aFVq1aYMWMGTp48qfIW048//ghfX1/Y2tqiV69eWLJkCZydnUt9/4i0YTLjLA0ePBiJiYnYt28f8vLyMHLkSIwZMwbr168vcZ9hw4YhNTUVO3bsQLVq1bB+/Xr0798fZ8+eRVhYWKmPWx4UCgXs7e0BABkZGbCzszNqPERVSXZ2tjSQnTEUJhmllZKSgt27d2P+/Plqvzuel1zs2LEDU6ZMAQBMmDABubm5+Pvvv2FnZ4crV67A3t4evr6+2LJlC/r06YO4uDg4OjpKY9q8++672LJlC1avXg0/Pz989tlniIiIwPXr1+Hq6qr1Naxbtw6tWrWChYUFAOD48eMYN24cFi1ahFdffRX79+/Hhx9+qMOdISodk0iWYmNjsXv3bpw5cwbh4eEAgGXLlqFbt25YvHgxfHx81O534sQJLF++HM2bNwcAfPDBB/jyyy8RHR2NsLCwUh+XiCq3+/fvY+XKlUY7/8iRI1GzZk2t6u7cuVP6wwoAunbtiunTp0MIgZCQEJ3PfffuXVy4cAFdu3YFANy+fRt9+vRBgwYNAACBgYFS3cLEx8PDQ0rAFAoFli9fjlWrVknH+PHHH7Fv3z78/PPPmD59usbzv/fee/jmm2+QmZmJF154ATt37pS2LVu2DF27dsW0adMAAHXq1MGJEydU6hAZgkk8houKioKzs7OU0ABAp06dIJfLcerUqRL3a9WqFTZu3IiUlBQolUps2LAB2dnZaNeuXZmOm5OTg/T0dJWFiMgY2rdvj5iYGGn5+uuvIYQo9fF27NiBNm3aSMnPpEmT8Mknn6B169aYM2cOLly4oHH/+Ph45OXloXXr1lKZhYUFmjdvjtjYWADAuHHjYG9vLy1FTZ8+HefPn8fevXthZmaGYcOGSdcTFxcn/fFb6Nl1IkMwiWQpKSkJHh4eKmXm5uZwdXVFUlJSifv99ttvyMvLg5ubG6ysrDB27Fhs27YNQUFBZTruggUL4OTkJC2+vr5luDoiotKzs7NDUFCQtHh7e6N27dqQyWS4evWqzsfbsWMHXn31VWn9jTfewI0bNzB06FBcvHgR4eHhWLZsWZlinjdvnkqCV1S1atVQp04dvPzyy9iwYQN27dr13D5WRIZm1MdwM2bMwKJFizTWKfxLpDQ+/PBDpKamYv/+/ahWrRq2b9+O/v374+jRo1KTcmnMnDlTep4PAOnp6UyYiCoRDw8PjBw50qjnLwtXV1dERETg22+/xaRJk4r1W0pNTVXbbykjIwOHDh3C8uXLVcp9fX0xbtw4jBs3DjNnzsSPP/6It956C5aWlgCezq9XqFatWrC0tMTx48fh5+cH4Ok0E2fOnME777wjXZ8216hUKgE8bc0HgODg4GIdxbXpOE5UVkZNlqZOnSq9PVGSwMBAeHl5FetsmZ+fj5SUFHh5eandLz4+Ht988w0uXbqEevXqAQAaNWqEo0eP4ttvv8WKFStKdVzg6dsnVlZWWlwhEZkia2trrfsMVVTffvstWrdujebNm2PevHlo2LAh8vPzsW/fPixfvlztH6K7d+9GnTp14O/vL5W988476Nq1K+rUqYPHjx/j0KFDCA0NBQD4+flBJpNh586d6NatG2xsbGBvb4/x48dj+vTpcHV1Rc2aNfHZZ58hMzMTo0aNKjHeU6dO4cyZM2jTpg1cXFwQHx+PDz/8ELVq1ULLli0BAG+99RZeeuklLFmyBJGRkTh48CD++usvTlVDBmfUx3Du7u4ICQnRuFhaWqJly5ZITU1FdHS0tO/BgwehVCrRokULtcfOzMwEAMjlqpdoZmYm/bVSmuMSEZmCwMBAnDt3Du3bt8fUqVNRv359vPzyyzhw4ECxlqNCv//+u8ojOOBpq9GECRMQGhqKLl26oE6dOvjuu+8AANWrV8fcuXMxY8YMeHp6YuLEiQCAhQsXok+fPhg6dCiaNGmC69evY8+ePXBxcSkxXltbW2zduhUdO3ZEcHAwRo0ahYYNG+LIkSPSH6etW7fGihUrsGTJEjRq1Ai7d+/G5MmTy/TmIJE2ZKIsPQHLUdeuXZGcnIwVK1ZIr/iHh4dLr/jfvXsXHTt2xJo1a9C8eXPk5eWhbt268Pb2xuLFi+Hm5obt27dj+vTp0l9B2hxXG+np6XByckJaWhocHR31cr0cOoCofGRnZ+PmzZsICAio0r908/Pz4enpib/++sukOk2PHj0aV69exdGjR40dClVQmv6Na/v72yQ6eAPAunXrEBISgo4dO6Jbt25o06YNfvjhB2l7Xl4e4uLipBYlCwsL7Nq1C+7u7oiMjETDhg2xZs0arF69WkqUtDkuEVFVkJKSgsmTJ6NZs2bGDkWjxYsX459//sH169exbNkyrF69GsOHDzd2WFTJmUzLUkXGliUi08WWJdPSv39/HD58GE+ePEFgYCDeeustjBs3zthhUQWmj5YlkxiUkoiICHg6JAxReTOZx3BERERExsBkiYiIiEgDJktEREREGjBZIiIiItKAyRIRERGRBnwbjoioFIQQ0rhuhRQKBTw9PQEAycnJxYb8sLW15dQcRCaILUtERKWQmZkJe3t7laUwUQIAT0/PYtufTa6qooSEBMhkMsTExAAADh8+DJlMhtTUVKPGVRW0a9dOmszYkEaMGIGePXsa/DzlickSEVEV9c8//+C1116Dr68vbGxsEBoaiq+++sqg5/T19UViYiLq169v0PNQcVu3bsXHH39s7DCKeTaBroj4GI6IqIzUPXIrVPTRXEUTHR0NDw8PrF27Fr6+vjhx4gTGjBkDMzMzaVLc50lISEBAQAC0nQzCzMwMXl5eZQmbSsnV1dXYIRhcXl4eLCws9H5ctiwREZWRnZ2dxsVQFAoFhg0bBnt7e3h7e+OLL74o9qglJycH7733Hnx9fWFlZYWgoCD8/PPPAIDXX38dX331Fdq2bYvAwEAMGTIEI0eOxNatWw0Wc0mtCMePH0fDhg1hbW2NF154AZcuXZK2rVq1Cs7Ozti5cyeCg4Nha2uLvn37IjMzE6tXr4a/vz9cXFwwadIkFBQUSPvdv38fkZGRsLGxQUBAANatWwd/f38sXbpUq1ivXbuGl156CdbW1qhbty727dsHmUyG7du3S3Xee+891KlTB7a2tggMDMSHH36IvLw8afs///yD9u3bw8HBAY6OjmjatCnOnj1bpuvS5LvvvkPt2rVhbW0NT09P9O3bV9r27GfD398fn3zyifQZ8vPzw44dO/DgwQP06NED9vb2aNiwoRQvAHz00Udo3LixyjmXLl0Kf3//EmPavXs32rRpA2dnZ7i5ueGVV15BfHy8tD0gIAAAEBYWBplMhnbt2gEAlEol5s2bhxo1asDKygqNGzfG7t27pf0KP0sbN25E27ZtYW1tjXXr1ml1n3TFliUiIhM1ffp0HDlyBL///js8PDzw/vvv49y5cyq/zIYNG4aoqCh8/fXXaNSoEW7evImHDx+WeMy0tDSjtEBMnz4dX331Fby8vPD+++8jMjIS//77r9RKkJmZia+//hobNmzAkydP0Lt3b/Tq1QvOzs7YtWsXbty4gT59+qB169YYMGAAgKd9Z+7du4dDhw7BwsICkyZNwv3797WKR6lUonfv3vD09MSpU6eQlpamtr+Pg4MDVq1aBR8fH1y8eBGjR4+Gg4MD3n33XQDA4MGDERYWhuXLl8PMzAwxMTEqLR+lua6SnD17FpMmTcIvv/yCVq1aISUlBUePHtW4z5dffolPP/0UH374Ib788ksMHToUrVq1wuuvv47PP/8c7733HoYNG4bLly+X+uUEhUKBKVOmoGHDhsjIyMDs2bPRq1cvxMTEQC6X4/Tp02jevDn279+PevXqwdLSEgDw1Vdf4YsvvsD333+PsLAw/O9//8Orr76Ky5cvo3bt2tLxZ8yYgS+++AJhYWGGm99RUJmlpaUJACItLU1vx8zIyBAABACRkZGht+MSkaqsrCxx5coVkZWVpdN+2v4bNdS/5SdPnghLS0vx22+/SWWPHj0SNjY24u233xZCCBEXFycAiH379ml1zOPHjwtzc3OxZ88ereO4efOm0OVXSWH98+fPCyGEOHTokAAgNmzYUOw6Nm7cKIQQYuXKlQKAuH79ulRn7NixwtbWVjx58kQqi4iIEGPHjhVC/N+1nz59WtoeGxsrAIgvv/zyuXHu2bNHmJubi7t370plf/31lwAgtm3bVuJ+n3/+uWjatKm07uDgIFatWqW2bmmuS5MtW7YIR0dHkZ6ernZ727Ztpc+GEEL4+fmJIUOGSOuJiYkCgPjwww+lsqioKAFAJCYmCiGEmDNnjmjUqJHKcb/88kvh5+cnrQ8fPlz06NGjxDgfPHggAIiLFy8KIYp/Jgr5+PiI+fPnq5Q1a9ZMvPnmmyr7LV26tMRzCaH537i2v7/5GI6IyATFx8cjNzcXLVq0kMpcXV0RHBwsrcfExMDMzAxt27Z97vEuXbqEHj16YM6cOejcubPGuvXq1ZPe8KtXrx4AqLz117VrV52vp2XLlsWuIzY2ViqztbVFrVq1pHVPT0/4+/vD3t5epayw5Sg2Nhbm5uZo2rSptD0kJATOzs5axRMbGwtfX1/4+PiojbHQxo0b0bp1a3h5ecHe3h4ffPABbt++LW2fMmUK3njjDXTq1AkLFy5UefxUmuvS5OWXX4afnx8CAwMxdOhQrFu37rlvYDZs2FDlPADQoEGDYmXatsipc+3aNbz22msIDAyEo6Oj9Miu6H16Vnp6Ou7du4fWrVurlLdu3VrlcwEA4eHhpY5NW0yWiIgqKRsbG63qXblyBR07dsSYMWPwwQcfPLf+rl27EBMTg5iYGOzatQsApPWYmBj89NNPZYpbnWc77cpkMrVlSqVS7+cuSVRUFAYPHoxu3bph586dOH/+PGbNmoXc3FypzkcffYTLly+je/fuOHjwIOrWrYtt27ZJ2/V5XQ4ODjh37hx+/fVXeHt7Y/bs2WjUqJHGYRmKnqvwMZu6ssLzy+XyYp35i/bRUicyMhIpKSn48ccfcerUKZw6dQoAVO5TWRiyX2AhJktERCaoVq1asLCwkH7xAMDjx4/x77//SusNGjSAUqnEkSNHSjzO5cuX0b59ewwfPhzz58/X6tx+fn4ICgpCUFAQ/Pz8AEBaDwoKQvXq1XW+npMnTxa7jtDQUJ2PUygkJAT5+fmIjo6WyuLi4rQezyk0NBR37txBYmKi2hgB4MSJE/Dz88OsWbMQHh6O2rVr49atW8WOVadOHUyePBl79+5F7969sXLlytJdlBbMzc3RqVMnfPbZZ7hw4QISEhJw8OBBvR3f3d0dSUlJKgmTplf+Hz16hLi4OHzwwQfo2LEjQkND8fjxY5U6hX2UinZid3R0hI+PD44fP65S9/jx46hbt64erkQ37OBNRFRGCoWiVNvKwt7eHqNGjcL06dPh5uYGDw8PzJo1C3L5//0N7O/vj+HDh+P111+XOnjfunUL9+/fR//+/XHp0iV06NABERERmDJlCpKSkgA8fb3f3d3dIHGXZN68eXBzc4OnpydmzZqFatWqlWlgw+DgYHTp0gVjx47F8uXLYW5ujnfeeUfr1rZOnTqhTp06GD58OD7//HOkp6dj1qxZKnVq166N27dvY8OGDWjWrBn+/PNPlVajrKwsTJ8+HX379kVAQAD+++8/nDlzBn369Cn1dWmyc+dO3LhxAy+99BJcXFywa9cuKJVKlUezZdWuXTs8ePAAn332Gfr27Yvdu3fjr7/+gqOjo9r6Li4ucHNzww8//ABvb2/cvn0bM2bMUKnj4eEBGxsb7N69GzVq1IC1tTWcnJwwffp0zJkzB7Vq1ULjxo2xcuVKxMTEGOyNN03YskREVEbqRutWN6q3vn3++ed48cUXERkZiU6dOqFNmzYqfXQAYPny5ejbty/efPNNhISEYPTo0VICt3nzZjx48ABr166Ft7e3tDRr1sxgMZdk4cKFePvtt9G0aVMkJSXhjz/+kFocSmvlypXw8fFB27Zt0bt3b4wZMwYeHh5a7SuXy7Ft2zZkZWWhefPmeOONN4q1vL366quYPHkyJk6ciMaNG+PEiRP48MMPpe1mZmZ49OgRhg0bhjp16qB///7o2rUr5s6dW6brKomzszO2bt2KDh06IDQ0FCtWrMCvv/4q9SvTh9DQUHz33Xf49ttv0ahRI5w+fRrTpk0rsb5cLseGDRsQHR2N+vXrY/Lkyfj8889V6pibm+Prr7/G999/Dx8fH/To0QMAMGnSJEyZMgVTp05FgwYNsHv3buzYsUPlTbjyIhPPPnwknaWnp8PJyQlpaWklZteaCM4xRWQ02dnZuHnzJgICAnR67VihUKh0wtVGRkaGwftXtGvXDo0bN9Z6LKGqxt/fH++8806pp/2QyWTYtm1bpZvOozLT9G9c29/ffAxXARTOMVUSdX+ZlseXLhGVzNbWFhkZGSpl2vyRQ0Smh8kSEVEpyGSyYsmQnZ2d1tN+kHGtW7cOY8eOVbvNz88Ply9fLueInu/o0aMah2V4Nnkn/WGyVMGY6hxTRFQxHD582NghVGgJCQkAgCdPnqiMUVWUprnFjJkMh4eHV+jJZiszJksVjKHnkiIioqdjEjk4OBg7DJ3Y2NggKCjI2GFUSXwbjogIxm0xICLD0ce/bSZLRFSlFZ2olYgqn8J/25oerz4PH8MRUZVmZmYGZ2dnae4rDstBVDkUDstz//59ODs7w8zMrNTHYrJERFWel5cXgLJNFkpEFZOzs7P0b7y0mCwRUZUnk8ng7e0NDw+P504KSkSmw8LCokwtSoWYLFUwxphjioieMjMz08sXKxFVLkyWKhiOo0RERFSx8G04IiIiIg2YLFUAhXNMFV1OnDghbZ82bRquX7+usp1zTBEREZUPPoarANTNMeXq6ir9bGlpCblczpG9iYiIjMBkWpZSUlIwePBgODo6wtnZGaNGjXrupIFJSUkYOnQovLy8YGdnhyZNmmDLli0qdfz9/SGTyVSWhQsXGvJStGJtba2ynpWVZaRIiIiIqjaTaVkaPHgwEhMTsW/fPuTl5WHkyJEYM2YM1q9fX+I+w4YNQ2pqKnbs2IFq1aph/fr16N+/P86ePYuwsDCp3rx58zB69GhpvSLMF2RjY6OyzmSJiIjIOEyiZSk2Nha7d+/GTz/9hBYtWqBNmzZYtmwZNmzYgHv37pW434kTJ/DWW2+hefPmCAwMxAcffABnZ2dER0er1HNwcICXl5e0VITHXWxZIiIiqhhMIlmKioqCs7MzwsPDpbJOnTpBLpfj1KlTJe7XqlUrbNy4ESkpKVAqldiwYQOys7PRrl07lXoLFy6Em5sbwsLC8PnnnyM/P19jPDk5OUhPT1dZ9O3ZsV6ys7P1fg4iIiJ6PpN4DJeUlAQPDw+VMnNzc7i6uiIpKanE/X777TcMGDAAbm5uMDc3h62tLbZt24agoCCpzqRJk9CkSRO4urrixIkTmDlzJhITE7FkyZISj7tgwQLMnTu37BemA7YsERERGYdRW5ZmzJhRrHP1s8vVq1dLffwPP/wQqamp2L9/P86ePYspU6agf//+uHjxolRnypQpaNeuHRo2bIhx48bhiy++wLJly5CTk1PicWfOnIm0tDRpuXPnTqlj1BZbloiIiIzDqC1LU6dOxYgRIzTWCQwMhJeXV7EJLvPz85GSklLi5Hjx8fH45ptvcOnSJdSrVw8A0KhRIxw9ehTffvstVqxYoXa/Fi1aID8/HwkJCQgODlZbx8rKClZWVs+5Ov1iyxIREZFxGDVZcnd3h7u7+3PrtWzZEqmpqYiOjkbTpk0BAAcPHoRSqUSLFi3U7pOZmQkAkMtVG8/MzMygVCpLPFdMTAzkcnmxx37GxmSJiIjIOEyig3doaCi6dOmC0aNH4/Tp0zh+/DgmTpyIgQMHwsfHBwBw9+5dhISE4PTp0wCAkJAQBAUFYezYsTh9+jTi4+PxxRdfYN++fejZsyeApx3Hly5din/++Qc3btzAunXrMHnyZAwZMgQuLi7Guly1mCwREREZh0l08AaAdevWYeLEiejYsSPkcjn69OmDr7/+Wtqel5eHuLg4qUXJwsICu3btwowZMxAZGYmMjAwEBQVh9erV6NatG4Cnj9M2bNiAjz76CDk5OQgICMDkyZMxZcoUo1yjJuyzREREZBwmkyy5urpqHIDS398fQgiVstq1axcbsbuoJk2a4OTJk3qL0ZBycnKgVCqLPVYkIiIiw+JvXhPC1iUiIqLyx2TJhLDfEhERUfljsmRCmCwRERGVPyZLJoTJEhERUfljsmRC2GeJiIio/DFZMiFsWSIiIip/TJZMCJMlIiKi8sdkyYQwWSIiIip/TJZMCPssERERlT8mSyaELUtERETlj8mSCWHLEhERUfljsmRC2LJERERU/pgsmRAmS0REROWPyZIJYbJERERU/pgsmZCCggLk5eUZOwwiIqIqhcmSiWHrEhERUflismRi+EYcERFR+WKyVEHZ2dlBCIGEhARYWlpK5WxZIiIiKl9Mlio4GxsblXUmS0REROWLyVIFZ21trbLOZImIiKh8MVmq4J5tWWKfJSIiovLFZKmCs7CwgJmZmbTOliUiIqLyxWTJBBRtXWKyREREVL6YLJmAoskSH8MRERGVLyZLJoAtS0RERMbDZMkEFH0jjskSERFR+WKyZAL4GI6IiMh4mCyZALYsERERGQ+TJRPwbMuSEMKI0RAREVUtTJZMQNFkSQiBnJwcI0ZDRERUtTBZMgGcH46IiMh4mCyZACZLRERExsNkyQQ8O5ku34gjIiIqPyaTLKWkpGDw4MFwdHSEs7MzRo0ahYyMDI37xMfHo1evXnB3d4ejoyP69++P5OTkMh+3vLFliYiIyHhMJlkaPHgwLl++jH379mHnzp34+++/MWbMmBLrKxQKdO7cGTKZDAcPHsTx48eRm5uLyMhIKJXKUh/XGJgsERERGY9MmMB76LGxsahbty7OnDmD8PBwAMDu3bvRrVs3/Pfff/Dx8Sm2z969e9G1a1c8fvwYjo6OAIC0tDS4uLhg79696NSpU6mOq056ejqcnJyQlpYmnUuflEolPv74Y2m9Q4cOePHFF/V+HiIioqpE29/fJtGyFBUVBWdnZymhAYBOnTpBLpfj1KlTavfJycmBTCaDlZWVVGZtbQ25XI5jx46V+riFx05PT1dZDEkul6tcB/ssERERlR+TSJaSkpLg4eGhUmZubg5XV1ckJSWp3eeFF16AnZ0d3nvvPWRmZkKhUGDatGkoKChAYmJiqY8LAAsWLICTk5O0+Pr6lvEKn4+T6RIRlT+FQgGZTAaZTAaFQmHscMqkMl1LeTNqsjRjxgzpf1xJy9WrV0t1bHd3d2zatAl//PEH7O3t4eTkhNTUVDRp0gRyedkue+bMmUhLS5OWO3fulOl42ij6RhxbloiIiMqPuTFPPnXqVIwYMUJjncDAQHh5eeH+/fsq5fn5+UhJSYGXl1eJ+3bu3Bnx8fF4+PAhzM3N4ezsDC8vLwQGBgJAqY9rZWWl8lisPLBliYhMhUKhgL29PQAgIyMDdnZ2Ro6IqGyMmiy5u7vD3d39ufVatmyJ1NRUREdHo2nTpgCAgwcPQqlUokWLFs/dv1q1atI+9+/fx6uvvqqX45YnJktERETGYRJ9lkJDQ9GlSxeMHj0ap0+fxvHjxzFx4kQMHDhQemPt7t27CAkJwenTp6X9Vq5ciZMnTyI+Ph5r165Fv379MHnyZAQHB2t93Iqi6GM4JktElRP7lBBVTEZtWdLFunXrMHHiRHTs2BFyuRx9+vTB119/LW3Py8tDXFwcMjMzpbK4uDjMnDkTKSkp8Pf3x6xZszB58mSdjltRFG1ZYp8lIiKi8mMyyZKrqyvWr19f4nZ/f388O2TUwoULsXDhwjIdt6Iomizl5uaioKAAZmZmRoyIiIioajCJx3DE+eGIiIiMhcmSieCUJ0RERMbBZMlEMFkiIiIyDiZLJoLJEhERkXEwWTIR7LNERERkHEyWTARbloiIiIyDyZKJsLS0hEwmk9aZLBEREZUPJksmQiaTccoTIiIiI2CyZEKYLBEREZU/JksmpOjM3Zw3ioiIqHwwWTIh9vb20s8ZGRlGjISo4uDks0RkaEyWTAhbloiIiMofkyUTUjRZyszMhFKpNGI0REREVQOTJRNSNFkSQrCTNxERUTlgsmRCivZZAthviYiIqDyY61I5NTUV27Ztw9GjR3Hr1i1kZmbC3d0dYWFhiIiIQKtWrQwVJ0G1ZQlgvyUiIn0SQiAzM1OlrOj3rLrvXFtbW5UBgyuKynQtFYFWydK9e/cwe/ZsrFu3Dj4+PmjevDkaN24MGxsbpKSk4NChQ1i8eDH8/PwwZ84cDBgwwNBxV0lMloiIDCczM7NYC35Rnp6excoyMjKKfTdXBJXpWioCrZKlsLAwDB8+HNHR0ahbt67aOllZWdi+fTuWLl2KO3fuYNq0aXoNlPgYjogqHrZgUFWgVbJ05coVuLm5aaxjY2OD1157Da+99hoePXqkl+BIlaWlJSwsLJCXlweALUtEZHyVtQUjOTm5xBgVCoXa66qoKtO1GItWydLzEqWy1ift2dnZITU1FQCTJSIiQ7Gzs6vwCZ22KtO1GItOHbwL/fLLL1ixYgVu3ryJqKgo+Pn5YenSpQgICECPHj30HSMVwWSJqHKojI+vKkoLhhACSqUSBQUFKCgoUPuzUqmUlsL10nRtOHnyJKytraFUKiGEkM5d+LOmsmeXwti1WX/ez9nZ2Tpfy5o1a2BlZSUdp/BYz/tZlzJ1dbTZBgBmZmYYO3asxjqGonOytHz5csyePRvvvPMO5s+fj4KCAgCAs7Mzli5dymTJwIo2dzNZIjJdlfHxlbYtGHfu3IG5uTny8vKKLfn5+dJ/S1oKCgpUfi5cCtdLO2Bvbm6uzvscPHgQlpaWpTqfIZXmWv77778KeS2FzM1L1b6jn3PrusOyZcvw448/omfPnli4cKFUHh4ezk7d5cDW1lb6mR28iciQlEolsrOzkZWVhaysLGRnZxdbClu6dfHLL79U6F/KRM/SOVm6efMmwsLCipVbWVmxpaMcPNuyJISo0M3yRPR85fX4SqlUQqFQICMjAwqFQuXnrKwsZGZmQqFQIDMzU0qOnqc0LRiVmUwmg1wulyZ3Vrf+vKXwOJrKiq4X/lz0v6X5/1KrVi1YW1urHEfbn3Up0/Q7S9M2MzOzErcZms7JUkBAAGJiYuDn56dSvnv3boSGhuotMFKv6BdqQUEBcnJypA83EZkmfXTAzcvLQ1paGtLS0pCeni4tT548wZMnT6Sk6Hn9QozNzMwMFhYWMDc3lxYLCwup3MzMDObm5ir/Lbo8WyaXy4v9LJfL1f6cnZ2NTz/9VKd43333XTg4OKgkLhWBQqHAuHHjdNpnwIABFfoxrzHpnCxNmTIFEyZMQHZ2NoQQOH36NH799VcsWLAAP/30kyFipCLUDUzJZImoanj8+DGSk5Px+PFjpKSkIDU1VVqe7SxeHiwtLUv1OG3IkCFwcnKChYWFymJubg653HizcD2vg726bcaOWRvaXguVTOdk6Y033oCNjQ0++OADZGZmYtCgQfDx8cFXX32FgQMHGiJGKkLdwJQcqoGoali+fLlB+vpYWlrCzs4Otra20mJjYyP9t+hibW0Na2trWFlZQS6XQ6FQ4P3339fpfDVr1qzwLRiVaeyhynQtxlKqruWDBw/G4MGDkZmZiYyMDHh4eOg7LioBpzwhqhzy8/MNfg4HBwdpsbe3V/m58NGfnZ0dLCws9HI+tmBQZVWqZCk/Px+HDx9GfHw8Bg0aBODp/HGOjo4aX4WlsmOyRGR68vPzkZycjLt37+LevXtISkrC3bt3y3RMMzMzuLi4wNnZGU5OTnB2dpZ+LvwuLu8OsabegmFra1vsLeOiHezVdcQv+oZyRVKZrqUi0DlZunXrFrp06YLbt28jJycHL7/8MhwcHLBo0SLk5ORgxYoVhoiT/j8bGxvI5XJpHBEOH0BUsQghkJaWhjt37uDOnTv477//kJycXGzsn9KMBdS9e3dUr14dLi4uUqdi0h+ZTKbx8aApjYRdma6lItA5WXr77bcRHh6Of/75R6WvTK9evTB69Gi9BkfFFf4DePLkCQC2LBEZmxACjx49QkJCAhISEnD79m3p36e2tH181ahRowr3C44tGFQV6JwsHT16FCdOnCjWydDf37/MzcqkHSZLRMaVnp6O69ev4+bNm0hISNCphdfGxgY+Pj5wcnKSykz58RVbMKgq0DlZKpxH51n//fcfHBwc9BKUOikpKXjrrbfwxx9/QC6Xo0+fPvjqq6809pGKj4/HtGnTcOzYMeTk5KBLly5YtmyZyheTv78/bt26pbLfggULMGPGDINdS1kV/eJhskRkeAUFBUhISMD169cRHx+PBw8eaLWfmZkZqlevLi0+Pj5wdnaGTCbjv10iE6JzstS5c2csXboUP/zwA4Cnf1VkZGRgzpw56Natm94DLDR48GAkJiZi3759yMvLw8iRIzFmzBisX79ebX2FQoHOnTujUaNGOHjwIADgww8/RGRkJE6ePKkyLsa8efNUHiEaMunTh6IJIvssERlGVlYWrl27hn///RfXrl3TakRkGxsb+Pv7w9fXF76+vvD29i6xkzUfXxGZDp2TpcWLF6NLly6oW7cusrOzMWjQIFy7dg3VqlXDr7/+aogYERsbi927d+PMmTMIDw8H8HSOum7dumHx4sXw8fEpts/x48eRkJCA8+fPw9HREQCwevVquLi44ODBg+jUqZNU18HBAV5eXgaJ3RDYskRkGJmZmbh69SouX76MmzdvPne0a2tra/j7+0uLh4eH1p2u+fiKyHTonCz5+vrin3/+wcaNG/HPP/8gIyMDo0aNwuDBg2FjY2OIGBEVFQVnZ2cpUQKATp06QS6X49SpU+jVq1exfXJyciCTyWBlZSWVWVtbQy6X49ixYyrJ0sKFC/Hxxx+jZs2aGDRoECZPnmzU2Y2fp+gXaG5uLvLy8vQ2TgpRVZOTk4MrV67gypUruHHjhsa31GQyGapXr45atWohKCgIPj4+FX70ZiIqO50ygry8PISEhGDnzp3SwJTlISkpqdjAl+bm5nB1dUVSUpLafV544QXY2dnhvffew6effgohBGbMmIGCggIkJiZK9SZNmoQmTZrA1dUVJ06cwMyZM5GYmIglS5aUGE9OTg5ycnKk9fT09DJeoW7UjbXk7OxcrjEQmTKlUokbN27gwoULiI2N1ThApKWlJYKCghAcHIygoCA+CiOqgnRKliwsLLSahVpbM2bMwKJFizTWiY2NLdWx3d3dsWnTJowfPx5ff/015HI5XnvtNTRp0kTlL8EpU6ZIPzds2BCWlpYYO3YsFixYoNIqVdSCBQswd+7cUsWlD892ameyRFWFEKLYHGjPm8/L1tZWejT2+PFjREdH48KFCxpf77exsUFoaChCQ0Ph7+9foVuaicjwdP4GmDBhAhYtWoSffvqpzF8gU6dOxYgRIzTWCQwMhJeXF+7fv69Snp+fj5SUFI19jTp37oz4+Hg8fPgQ5ubmcHZ2hpeXFwIDA0vcp0WLFsjPz0dCQgKCg4PV1pk5c6ZKkpWeng5fX1+N16FPz7YssZM3VRWZmZka34BV9wp+eno67t69i+joaFy/fr3Efa2trREaGop69erB39+/3Ee/JqKKS+ds58yZMzhw4AD27t2LBg0aFPvFvXXrVq2P5e7uDnd39+fWa9myJVJTUxEdHY2mTZsCAA4ePAilUokWLVo8d/9q1apJ+9y/fx+vvvpqiXVjYmIgl8s1zndnZWVVYqtTeeCUJ0Ta++abb0p8k00ul6N27dpo1KgRateuzRYkIlJL528GZ2dn9OnTxxCxlCg0NBRdunTB6NGjsWLFCuTl5WHixIkYOHCg9Cbc3bt30bFjR6xZswbNmzcHAKxcuRKhoaFwd3dHVFQU3n77bUyePFlqMYqKisKpU6fQvn17ODg4ICoqCpMnT8aQIUPg4uJSrteoC7YsEal/tb5Q0VfwMzIyig2i6+HhgSZNmqB+/fp844yInkvnZGnlypWGiOO51q1bh4kTJ6Jjx47SoJRff/21tD0vLw9xcXEq/Rni4uIwc+ZMpKSkwN/fH7NmzcLkyZOl7VZWVtiwYQM++ugj5OTkICAgAJMnT1Z5xFYRyeVy2NjYICsrCwBblqhq0vXVejMzM9SvXx9NmzZFjRo1OK8aEWlN52SpQ4cO2Lp1a7EOxenp6ejZs6c0AKS+ubq6ljgAJfB0JO5nx0RZuHAhFi5cWOI+TZo0wcmTJ/UWY3myt7dnskSlplAopL4/GRkZlbp1xdHRES+++CLCwsIMNrwJEVVuOidLhw8fVvv8Pzs7G0ePHtVLUPR8dnZ20pQLTJaISvbmm29W+FH5iahi0zpZunDhgvTzlStXVMY3KigowO7du1G9enX9RkclKvpGEJMlopJx0EgiKiutk6XGjRtDJpNBJpOhQ4cOxbbb2Nhg2bJleg2OSlZ0YDx28Kaq4OHDh9i1a5exwyCiKkjrZKlwnqTAwECcPn1a5ZV/S0tLeHh4cFySclS0ZSkrKwsFBQW8/1QpZWZm4tChQ4iOjlYZOZ+IqLxonSz5+fkhLy8Pw4cPh5ubG/z8/AwZFz3Hsx1yMzMz2S+DKpWCggKcOXMGR44cUTtzgKbHz3w0TUT6pPN0J9u2bcPs2bMNFQ9pSd3AlEyWqLKIj4/HX3/9hUePHpVYR91o3UREhqBzz8cePXpg+/btBgiFdPHslA/st0SVwZMnT7B582asXbu2WKJUrVo1DBgwwEiREVFVpvPQAbVr18a8efNw/PhxNG3atFgLx6RJk/QWHJWMU55QZaJUKnH27FkcPHiwWL8ka2trtGvXDuHh4ZDL5cX+MCg6Wre6Ub2LvgxBRFQaOidLP//8M5ydnREdHY3o6GiVbTKZjMlSOWGyRJXF/fv3sWPHDty9e7fYtvDwcLRv314l4dE0gKauo3oTPY+dnV2xAY9NVWW6lvKmc7J08+ZNQ8RBOrKwsIClpaU0QCgfw5GpUSqVOH78OI4cOYKCggKVbd7e3njllVekuR/JtPCXMlU2nGLbhNnb2yMlJQUAVObEI6roHjx4gO3bt+PevXsq5ZaWlmjfvj2aN2/OwSSJqMIoVbL033//YceOHbh9+3axqU+WLFmil8Do+ezs7KRkiS1LZAqEEIiKisLBgweLtSbVqVMH3bt3h6Ojo5GiIyJST+dk6cCBA3j11VcRGBiIq1evon79+khISIAQAk2aNDFEjFSCon0z2GeJKronT55g+/btuHHjhkq5lZUVunTpgkaNGkEmkxkpOiKikumcLM2cORPTpk3D3Llz4eDggC1btsDDwwODBw9Gly5dDBEjlYDJEpmKf//9F7///nuxx8VBQUGIjIxkaxIRVWg6J0uxsbH49ddfn+5sbo6srCzY29tj3rx56NGjB8aPH6/3IEm9ZyfTFULwL3OqUPLz87F//36cOnVKpdzCwgIRERFo0qQJP7NEVOHpnCzZ2dlJ/ZS8vb0RHx+PevXqAXg60SWVn6ItS0qlEllZWRxThiqM1NRUbNq0qVgnbi8vL/Tp0wfVqlUzUmQVF98iI6qYdE6WXnjhBRw7dgyhoaHo1q0bpk6diosXL2Lr1q144YUXDBEjleDZ8WQyMjKYLFGFcP36dWzduhVZWVkq5S+88AI6duwIc3O+iEtEpkPnb6wlS5ZIb17NnTsXGRkZ2LhxI2rXrs034crZs/080tPT4eHhYaRoiJ6+7fb333/j8OHDKuW2trbo2bMnateubZzAiIjKQOdkKTAwUPrZzs4OK1as0GtApD0nJyeV9bS0NCNFQgRkZ2dj69atuHbtmkp59erV0a9fv2KfVyIiU1HqtvDc3Fzcv38fSqVSpbxmzZplDoq0Y2dnB7lcLv0/SE9PN3JEVFWlpKTg119/LdZvMTw8HBEREXzsRkQmTedvsH///RejRo3CiRMnVMoL38R6dqA5Mhy5XA5HR0ekpqYCYLJExpGQkIDffvtNpX+Subk5IiMj0bBhQyNGRkSkHzonSyNHjoS5uTl27twJb29vvvZrZEWTJT6Go/IWHR2NXbt2qbQwOzk5YeDAgfDy8jJiZERE+qNzshQTE4Po6GiEhIQYIh7SUdF+IEyWqLwIIbB3716cPHlSpdzX1xf9+/dXGQOMiMjU6Zws1a1bl+MpVSBF34hLT0/nwJRkcPn5+di+fTsuX76sUt6oUSO88sor7J9ERJWOVt9qRfvCLFq0CO+++y4+/fRTNGjQABYWFip1OW1B+SraspSfn4/MzMxi4y8R6UtWVhY2btyIW7duqZR36tQJrVq1YqJORJWSVsmSs7OzypegEAIdO3ZUqcMO3sahbvgAJktkCGlpaVi3bh0ePHgglZmZmaF3796oW7euESMjIjIsrZKlQ4cOGToOKiV1A1P6+PgYKRqqrB48eIBffvkFT548kcqsra0xcOBA+Pn5GTEyIiLD0ypZatu2raHjoFLiwJRkaImJiVi7di0yMzOlMicnJwwePBju7u5GjIyIqHzItal0+/ZtnQ569+7dUgVDurO2tlbpN8ZkifTp9u3bWL16tUqi5OnpiVGjRjFRIqIqQ6tkqVmzZhg7dizOnDlTYp20tDT8+OOPqF+/PrZs2aK3AEkzmUym0rrEgSlJX+Lj4/HLL78gJydHKvP19cWIESPg4OBgxMiIiMqXVo/hrly5gvnz5+Pll1+GtbU1mjZtCh8fH1hbW+Px48e4cuUKLl++jCZNmuCzzz5Dt27dDB03FeHk5CQN58CWJdKHq1evYvPmzSovbAQGBmLAgAGwtLQ0YmREROVPq2TJzc0NS5Yswfz58/Hnn3/i2LFjuHXrFrKyslCtWjUMHjwYERERqF+/vqHjJTWeHWuJqCghhMpjNABQKBRqfwaeJkp//vknhBBSWUhICPr06cMxlIioStLpm8/GxgZ9+/ZF3759DRUPlULRx3BPnjyBUqmEXK7VE1aqAjIzMzWOqO3p6Vms7P3335dakBo2bIgePXrwM0VEVZbJfPvNnz8frVq1gq2tLZydnbXaRwiB2bNnw9vbGzY2NujUqROuXbumUiclJQWDBw+Go6MjnJ2dMWrUKGRkZBjgCgynaMuSEELl9W6isggLC0PPnj2ZKBFRlWYy34C5ubno168fxo8fr/U+n332Gb7++musWLECp06dgp2dHSIiIpCdnS3VGTx4MC5fvox9+/Zh586d+PvvvzFmzBhDXILBcPgA0lZycjIyMjLULsnJySp1w8LCEBkZyVG5iajKM5kOCHPnzgUArFq1Sqv6QggsXboUH3zwAXr06AEAWLNmDTw9PbF9+3YMHDgQsbGx2L17N86cOYPw8HAAwLJly9CtWzcsXrzYZAZ3ZLJE2rKzs9NqhPdGjRoxUSIi+v9MpmVJVzdv3kRSUhI6deoklTk5OaFFixaIiooCAERFRcHZ2VlKlICnc1zJ5XKcOnWqxGPn5OQgPT1dZTEmdaN4E5VFt27dTCZRsrOzgxACQghO9UNEBqFzsvTsmzMVVVJSEoDinVc9PT2lbUlJSfDw8FDZbm5uDldXV6mOOgsWLICTk5O0+Pr66jl63VhYWMDW1lZaZ8sSlZWpJEpEROVB52TJ09MTr7/+Oo4dO1bmk8+YMQMymUzjcvXq1TKfR99mzpyJtLQ0ablz546xQ+LwAURERAaic5+ltWvXYtWqVejQoQP8/f3x+uuvY9iwYaXq3zN16lSMGDFCY53AwECdjwsAXl5eAJ52aPX29pbKk5OT0bhxY6nO/fv3VfbLz89HSkqKtL86VlZWsLKyKlVchuLk5CS1hrFliYiISH90blnq2bMntm/fjrt372LcuHFYv349/Pz88Morr2Dr1q3Iz8/X+lju7u4ICQnRuJR2tOCAgAB4eXnhwIEDUll6ejpOnTqFli1bAgBatmyJ1NRUREdHS3UOHjwIpVKJFi1alOq8xlK0ZYnJEhERkf6UuoO3u7s7pkyZggsXLmDJkiXYv38/+vbtCx8fH8yePbvYiMFldfv2bcTExOD27dsoKChATEwMYmJiVMZECgkJwbZt2wA87XPxzjvv4JNPPsGOHTtw8eJFqQWsZ8+eAIDQ0FB06dIFo0ePxunTp3H8+HFMnDgRAwcONJk34QoVfSMuKysLeXl5RoyGKiqFQqFxISKi4ko9dEBycjJWr16NVatW4datW+jbty9GjRqF//77D4sWLcLJkyexd+9evQU6e/ZsrF69WloPCwsDABw6dAjt2rUDAMTFxam0qrz77rtQKBQYM2YMUlNT0aZNG+zevRvW1tZSnXXr1mHixIno2LEj5HI5+vTpg6+//lpvcZcXdcMHVKtWzUjRUEWlbrRuIiLSTCaKTgClha1bt2LlypXYs2cP6tatizfeeANDhgxRGVU7Pj4eoaGhyM3N1Xe8FVJ6ejqcnJyQlpZW7DX+8nL79m2sXLlSWh86dGip+3tR5fLkyROdP5cZGRl8DZ+IKj1tf3/r3LI0cuRIDBw4EMePH0ezZs3U1vHx8cGsWbN0PTSVAQemJHWEEDh06BDef/99qczGxgZ9+vRB3bp1ATxtJX42MSo6FAURUVWnc7KUmJj43C9SGxsbzJkzp9RBke4cHBwgk8mkmeKZLBEA/P333zh//rz0ooSFhQWGDRsGFxcXqY62o3oTEVVVOidL+fn5asfxkclksLKyKvXba1Q2crkcDg4O0v8bJkt04cIFHD58WFqXyWTo168fatSowc7cREQ60DlZcnZ21ji6b40aNTBixAjMmTOHM5WXMycnJylZ4sCUVdutW7ewY8cOlbLIyEjUrl3bSBEREZkunZOlVatWYdasWRgxYgSaN28OADh9+jRWr16NDz74AA8ePMDixYthZWWl0k+CDM/JyUkaTZwtS1XXo0ePsHHjRhQUFEhlL774ovQGKRER6UbnZGn16tX44osv0L9/f6ksMjISDRo0wPfff48DBw6gZs2amD9/PpOlcvbslCdCCM7xVcVkZmZi/fr1yMrKksrq16+P9u3bGzEqIiLTpvNzshMnTqj9CzUsLAxRUVEAgDZt2uD27dtlj450UjRZysvLQ3Z2thGjofJWUFCA3377DSkpKVKZr68vevTowaSZiKgMdE6WfH198fPPPxcr//nnn+Hr6wvg6WOAom/bUPng8AFVlxACf/75J27duiWVubi4YMCAATA3L/XYs0REhFI8hlu8eDH69euHv/76Sxpn6ezZs7h69So2b94MADhz5gwGDBig30jpudQlS5omBKbK4/Tp0zh//ry0bm1tjUGDBnFIACIiPdA5WXr11VcRFxeH77//HnFxcQCArl27Yvv27fD39wcAjB8/Xq9BknbYslQ1xcfHY8+ePdJ64RABnO6GiEg/dEqW8vLy0KVLF6xYsQILFiwwVExUSjY2NjA3N0d+fj4ADh9QFTx69AibN29G0VmLIiIiONUNEZEe6dRnycLCAhcuXDBULFRGMplMpXWJLUuVW3Z2NjZs2KDSkb9JkybSkB5ERKQfOnfwHjJkiNoO3lQxFE2WHj9+bMRIyJCEENi2bRsePnwoldWsWRPdunXjm29ERHpWqulO/ve//2H//v1o2rRpsQ6kS5Ys0VtwpDtXV1fcuHEDwNNHNBxrqXI6cuQI/v33X2ndyckJ/fv3h5mZmRGjIiKqnHROli5duoQmTZoAgMqXNQD+Uq4A3NzcpJ+zs7ORlZXFGeQrmX///RdHjhyR1s3NzTFw4EC++UZEZCA6J0uHDh0yRBykJ0WTJeBp6xKTpcrj0aNH2Lp1q0rZq6++yiEiiIgMqNQz3V6/fh179uyRplUo+jYOGc+zr4sX7dNCpi03NxcbN25ETk6OVNaiRQs0aNDAiFEREVV+OidLjx49QseOHVGnTh1069YNiYmJAIBRo0Zh6tSpeg+QdOPk5KTSb+XRo0dGjIb0RQiBHTt24MGDB1KZn58fXn75ZSNGRURUNeicLE2ePBkWFha4ffu2yuOdAQMGYPfu3XoNjnQnl8vh6uoqrTNZqhxOnTqFy5cvS+sODg7o27cvO3QTEZUDnfss7d27F3v27EGNGjVUymvXrq0yLxUZj5ubm9QCwWTJ9N25cwf79u2T1s3MzNC/f3/Y29sbMSoioqpD55YlhUKhtsNwSkoKrKys9BIUlU3RlqWUlBQolUojRkNloVAosHnzZpX/hxEREcX+WCEiIsPROVl68cUXsWbNGmldJpNBqVTis88+Q/v27fUaHJVO0U7eBQUFHMnbRCmVSmzbtk1l2poGDRogPDzciFEREVU9Oj+G++yzz9CxY0ecPXsWubm5ePfdd3H58mWkpKTg+PHjhoiRdKRu+AAXFxcjRVO5KBQK6fFXRkaGQcc2+vvvvxEfHy+tu7u745VXXuF4ZkRE5UznlqX69evj33//RZs2bdCjRw8oFAr07t0b58+fR61atQwRI+lIXbJEpiU+Pl5l4EkLCwv069cPlpaWRoyKiKhq0rllCXj6evqsWbP0HQvpia2tLaytraUJVpksmZYnT56oHXjS3d3dSBEREVVtpUqWUlNTcfr0ady/f79Y5+Fhw4bpJTAqPZlMBjc3N9y9excAkyVTolQqsXXrVmRmZkpl4eHhqF+/vhGjIiKq2nROlv744w8MHjwYGRkZcHR0VOk/IZPJmCxVEEyWTNPff/+NhIQEad3b2xsRERHGC4iIiHTvszR16lS8/vrryMjIQGpqKh4/fiwtKSkphoiRSqFov6W0tDTk5eUZMRrSxs2bN1X6KVlaWqJv374wNy9VAzAREemJzsnS3bt3MWnSJE7OWsE928mbiWzFplAoivVTioyMVBkzi4iIjEPnZCkiIgJnz541RCykR3wjznQIIbBt2zZkZGRIZU2aNGE/JSKiCkLn9v3u3btj+vTpuHLlCho0aAALCwuV7a+++qregqPSe7ZFgslSxXX8+HGV8ZQ8PDzQpUsXI0ZERERF6ZwsjR49GgAwb968YttkMhkKCgrKHhWVmaWlJRwdHaXRn5ksVUx3797FoUOHpHULCwv07du32B8hRERkPDonS5xnzHS4ubkxWarAcnJysGXLFpV/U127di2X8ZTs7OwghDD4eYiIKgOd+ywZy/z589GqVSvY2trC2dlZq32EEJg9eza8vb1hY2ODTp064dq1ayp1/P39IZPJVJaFCxca4ArKX9F+S0yWKhYhBP788088fvxYKqtfvz4aN25svKCIiEgtrZOlbt26qUzIunDhQqSmpkrrjx49Qt26dfUaXFG5ubno168fxo8fr/U+n332Gb7++musWLECp06dgp2dHSIiIqSRrQvNmzcPiYmJ0vLWW2/pO3yjKJosZWVlqQx0SMZ14cIFXLx4UVp3dnZG9+7dOe8bEVEFpHWytGfPHuTk5Ejrn376qcrr6Pn5+YiLi9NvdEXMnTsXkydPRoMGDbSqL4TA0qVL8cEHH6BHjx5o2LAh1qxZg3v37mH79u0qdR0cHODl5SUthpwctTzxjbiK6dGjR/jzzz+ldZlMhj59+sDa2tqIURERUUm0Tpae7d9Q0fs73Lx5E0lJSejUqZNU5uTkhBYtWiAqKkql7sKFC+Hm5oawsDB8/vnnyM/P13jsnJwcpKenqywVEZOliqegoABbtmxRGSS0ffv2qFGjhhGjIiIiTSrt0MBJSUkAAE9PT5VyT09PaRsATJo0CU2aNIGrqytOnDiBmTNnIjExEUuWLCnx2AsWLMDcuXMNE7geOTs7Qy6XSx2IHz58aOSI6NChQ0hMTJTW/f390bp1ayNGREREz6N1y1Jh5+dny8pixowZxTpXP7tcvXq1TOd4nilTpqBdu3Zo2LAhxo0bhy+++ALLli1TeeT4rJkzZyItLU1a7ty5Y9AYS0sul6uMt8RRvI0rISEBx48fl9ZtbGzQq1cvyOUm854FEVGVpHXLkhACI0aMgJWVFQAgOzsb48aNk/r3aEouSjJ16lSMGDFCY53AwECdjwsAXl5eAIDk5GR4e3tL5cnJyRrfOGrRogXy8/ORkJCA4OBgtXWsrKyk+1DRubm5SS1KbFkynqysLGzbtk2lLDIyEo6OjkaKiIiItKV1sjR8+HCV9SFDhhSrM2zYMJ1O7u7ubrAxZQICAuDl5YUDBw5IyVF6ejpOnTql8Y26mJgYyOVyeHh4GCSu8la031JKSgqUSiVbMsqZEAI7d+5U6dvWpEkThIaGGjEqIiLSltbJ0sqVKw0Zx3Pdvn0bKSkpuH37NgoKChATEwMACAoKgr29PQAgJCQECxYsQK9evSCTyfDOO+/gk08+Qe3atREQEIAPP/wQPj4+6NmzJwAgKioKp06dQvv27eHg4ICoqChMnjwZQ4YMgYuLi5GuVL+KJksFBQVIT0/Xepwq0o9//vkHV65ckdbd3NwQERFhxIiIiEgXJtPBe/bs2Vi9erW0HhYWBuBph9l27doBAOLi4lTGgnr33XehUCgwZswYpKamok2bNti9e7f0iraVlRU2bNiAjz76CDk5OQgICMDkyZMxZcqU8rswA6tWrZrK+oMHD5gslaOUlBT89ddf0rpcLkfv3r1haWlpxKiIiEgXMlHRxwAwAenp6XByckJaWlqF64OSnZ2NRYsWSevt27fHSy+9ZMSITJtCoZBaMjMyMjSOyaVUKrFy5Ur8999/UlnHjh3Rpk0bg8dJRETPp+3vb3ZeqeSsra1VHikWfW2dDOvo0aMqiZKfnx9atWplxIiIiKg0mCxVAUXfBmSyVD7u3r2LI0eOSOtWVlYcJoCIyETxm7sKKJospaWlcY44A8vNzcXWrVtVRrnv3r07nJycjBgVERGVFpOlKqBosgSwdcnQ9uzZozIAaP369bWe05CIiCoeJktVAJOl8hMXF4dz585J646OjujWrZsRIyIiorJislQF2NraqjwCYrJkGBkZGdixY4dKWc+ePWFjY2OkiIiISB+YLFUR7ORtWEII/PHHHyr9wVq2bImAgAAjRkVERPpgMoNSUtl4e3tLkxI/fvwYWVlZbPF4DiFEsc7wCoVC7c/nz5/HpUuXYGFhAZlMBk9PT3To0KHcYiUiIsNhslRFPNtvKSkpia0ez5GZmSkNQKmOp6dnsbL3338fNjY26NWrF8zN+c+LiKgy4GO4KoKdvMtPhw4d1CZSRERkmvinbxVhb28PBwcHPHnyBACTJV0lJyeXOLWJQqGQkiM/Pz+0bNmyPEMjIiIDY7JUhfj4+CAuLg4AkyVd2dnZaZwHrtArr7wCmUxWDhEREVF54WO4KsTLy0v6+dGjR8jJyTFiNJUTR+kmIqp8mCxVIT4+PirrSUlJRoqEiIjIdDBZqkKe7eR97949I0VCRERkOpgsVSEODg4qr8KzZYmIiOj5mCxVMUVbl9iyRERE9Hx8G66K8fb2xrVr1wAADx8+RG5uLiwtLY0cVcVXdLRuXbYREZHpY7JUxagbybtmzZpGisZ0cJBJIqKqi4/hqhiO5E1ERKQbtixVMY6OjrC1tZUmiGWyVDJbW1tkZGTgxIkTOHz4MAAgNzcXixcvBqB+VG9bW9vyDpOIiAyMyVIVI5PJ4O3tjfj4eADA3bt3jRxRxSWTyZCeno6TJ09K/bqK9u/SdlRvIiIybXwMVwXVqFFD+vnhw4fSfHGkKi8vD1u3boVSqZTKOnfubMSIiIjIGJgsVUGBgYEq6zdu3DBSJBXb/v378fDhQ2k9NDQUDRo0MGJERERkDEyWqqDq1aurPE66efOmEaOpmOLj43H69Glp3d7enpPkEhFVUUyWqiAzMzP4+flJ6zdu3IAQwogRVSwKhQLbt29XKevRowc7bxMRVVFMlqqooo/injx5gkePHhkxmopDCIEdO3YgIyNDKmvevDmCgoKMGBURERkTk6UqKiAgQGWd/ZaeOnv2LP79919p3cPDAy+//LIRIyIiImNjslRFeXh4qLz2bsh+SwqFAjKZDDKZrEJPDfLgwQPs3btXWjczM0Pv3r1hbs4RNoiIqjImS1WUTCZTeRR38+ZNlVfkq5r8/Hxs2bIF+fn5UtnLL7/MaU6IiIjJUlVW9FFcTk5OlR7N+8CBA0hOTpbWg4KC0Lx5cyNGREREFQWTpSqM4y09df36dZw8eVJat7W1RY8ePThMABERAWCyVKU5OTnB1dVVWq+KyVJGRobaYQLs7e2NExAREVU4JpMszZ8/H61atYKtrS2cnZ212mfr1q3o3Lkz3NzcIJPJEBMTU6xOdnY2JkyYADc3N9jb26NPnz4qj2Mqu6KtS3fu3EFeXp4RoylfQghs27ZNpdN5s2bNUKdOHSNGRUREFY3JJEu5ubno168fxo8fr/U+CoUCbdq0waJFi0qsM3nyZPzxxx/YtGkTjhw5gnv37qF37976CNkkFO23VFBQgNu3bxsxmvJ1/PhxldY0T09Pzv1GRETFmMw70XPnzgUArFq1Sut9hg4dCgBISEhQuz0tLQ0///wz1q9fjw4dOgAAVq5cidDQUJw8eRIvvPBCmWI2BerGW6pVq5aRoik/d+7cwcGDB6V1CwsL9O3bl8MEEBFRMSbTsmQI0dHRyMvLQ6dOnaSykJAQ1KxZE1FRUSXul5OTg/T0dJXFVNnY2MDHx0darwrzxGVnZ2PLli0qU7x069YN1apVM2JURERUUVXpZCkpKQmWlpbF+kB5enoiKSmpxP0WLFgAJycnafH19TVwpIZVtHUpMTERmZmZRozGsAqnM0lLS5PKGjRogEaNGhkxKiIiqsiMmizNmDFDGtm5pOXq1avGDFGtmTNnIi0tTVru3Llj7JDK5NkhBOLj440UieGdOnUKsbGx0rqLiwu6d+/OYQKIiKhERu2gMXXqVIwYMUJjnWd/keuTl5cXcnNzkZqaqtK6lJycDC8vrxL3s7KygpWVlcHiKm++vr6wsLCQ3oS7dOkSGjRoYOSo9O+///7Dvn37pHUzMzP07du3Uv2/JCIi/TNqsuTu7g53d3ejnb9p06awsLDAgQMH0KdPHwBAXFwcbt++jZYtWxotrvJmYWGB0NBQXLhwAcDTQRoVCoXK3HGmLjMzE5s2bVKZ0iUiIkKlvxYREZE6JtNn6fbt24iJicHt27dRUFCAmJgYxMTEICMjQ6oTEhKCbdu2SespKSmIiYnBlStXADxNhGJiYqT+SE5OThg1ahSmTJmCQ4cOITo6GiNHjkTLli2rxJtwRTVs2FD6WalU4tKlS0aMRr8Kx1Mq2hG/fv36CA8PN2JURERkKkwmWZo9ezbCwsIwZ84cZGRkICwsDGFhYTh79qxUJy4uTqXj7o4dOxAWFobu3bsDAAYOHIiwsDCsWLFCqvPll1/ilVdeQZ8+ffDSSy/By8sLW7duLb8LqyACAgLg4OAgrRe2MlUGR48exfXr16V1Nzc3vPLKK+ynREREWpGJou9PU6mkp6fDyckJaWlpcHR0NHY4pbZv3z6cOHFCWn/zzTf18phUoVBI04dkZGSU6+O9GzduYO3atdIwAebm5hg9ejQ8PDxKdTxjXgsREemXtr+/TaZliQyv6KM4wPRblx4/fozNmzerjKf0yiuvlDpRIiKiqonJEkk8PT3h6ekprV+8eBGm2vCYm5uLDRs2ICsrSypr0qQJx1MiIiKdMVkiFUVbl9LS0nDr1i0jRlM6Qgj8/vvvuH//vlTm6+uLbt26GTEqIiIyVUyWSEWDBg1UOj7/888/RoymdI4dOya9AQkADg4O6N+/P8zMzIwYFRERmSomS6TCwcFBZSDQK1euSINVmoJ///1XZYJcMzMzDBgwQOqUTUREpCsmS1RM0Udxubm5iIuLM2I02ktKSsKWLVtUyrp3747q1asbKSIiIqoMmCxRMSEhIbC0tJTWTeFRXFpaGtavX4/c3FyprHnz5ggLCzNiVEREVBkwWaJiLC0tERoaKq1fv34d9+7dM2JEmmVnZ2P9+vV48uSJVFarVi107txZ7+eys7ODEAJCCI6xRERURTBZIrWaNWumsl60H1BFUlBQgE2bNqm8+ebp6Yl+/fqxQzcREekFkyVSq3r16ggJCZHW4+PjkZCQYLyA1BBCYOfOnbhx44ZU5ujoiEGDBsHKysqIkRERUWXCZIlK1L59e5X1gwcPVphBKoUQ2L9/P2JiYqQyKysrDBo0yKSnnCEiooqHyRKVyMPDQ2XE6zt37uDatWtGjOj/HDp0SGUeO7lcjv79+6uMQE5ERKQPTJZIo7Zt20Iu/7+PSUVoXTpy5AiOHj2qUvbqq6+qjA9FRESkL0yWSCMXFxc0bdpUWk9OTsalS5eMFs+xY8dw+PBhlbLIyEjO+UZERAbDZIme68UXX4S5ubm0fvjwYRQUFJR7HCdOnMCBAwdUyrp27YomTZqUeyxERFR1MFmi53JwcECLFi2k9ZSUFPz999/ldn4hBPbs2YN9+/aplHfu3BnNmzcvtziIiKhqYrJEWmndurXK6/h///23ymS1hpKXl4dNmzbh5MmTKuUdO3ZEy5YtDX5+IiIiJkukFRsbG3Tv3l2lbPv27UhKSlIpE0JAoVAUWwqp21ZSh/HMzEysWbMGsbGxKuURERFo06aNnq6MiIhIM/PnVyF6qkGDBkhOTsbx48cBPG312bBhA0aPHi1N/ZGZmQl7e/sSj6Hu1f6MjIxiU4ckJSVh06ZNSElJkcrMzc3Ru3dvlalYiIiIDI0tS6STDh06oHbt2tJ6WloaNm3apLcO30qlEkePHsWPP/6okijZ2Nhg2LBhTJSIiKjcMVkincjlcvTu3RvVqlWTym7duoWtW7ciOztbpW5ycjIyMjLULsnJycWO/ejRI6xcuRIHDx6EUqmUyl1dXTFq1Cj4+voa7sKIiIhKwMdwpDNra2sMHDgQP/30k5QgXblyBXfu3EHHjh2lenZ2dsUer6mTlZWF6OhoHD16FPn5+SrbQkJCEBkZCVtbW/1eBBERkZaYLFGpuLm5oW/fvli/fr3UCvTkyRP89ttvOh/rm2++gUwmUymzsrJC165d0bBhw2LbiIiIyhMfw1Gp1apVCyNGjICrq2uZjpOXl6eyHhgYiPHjx6NRo0ZMlIiIyOjYskRl4uvri3HjxuHgwYPFxkLSlZOTE1588UU0adKESRIREVUYTJaozCwsLBAREYG6deuW6jGcj48P2rVrh9DQUJVJe4mIiCoCmTD2FPKVQHp6OpycnJCWlgZHR0djh2NUqampcHFxAaB+/KRCCoVCGo9JUz0iIiJD0fb3N1uWSK8sLCykn4uO3P0sTduIiIgqEiZLZDDqRusmIiIyNewgQkRERKQBW5ZIr2xtbZGRkaFSplAopFam5OTkYv2TOOAkERFVZEyWSK9kMpnGztrajupNRERUUfAxHBEREZEGJpMszZ8/H61atYKtrS2cnZ212mfr1q3o3Lkz3NzcIJPJEBMTU6xOu3btIJPJVJZx48bpN3giIiIyWSaTLOXm5qJfv34YP3681vsoFAq0adMGixYt0lhv9OjRSExMlJbPPvusrOESERFRJWEyfZbmzp0LAFi1apXW+wwdOhQAkJCQoLGera0tvLy8ShsaERERVWIm07JkSOvWrUO1atVQv359zJw5E5mZmRrr5+TkID09XWUhIiKiyslkWpYMZdCgQfDz84OPjw8uXLiA9957D3Fxcdi6dWuJ+yxYsEBq6SIiIqLKzagtSzNmzCjWufrZ5erVqwaNYcyYMYiIiECDBg0wePBgrFmzBtu2bUN8fHyJ+8ycORNpaWnScufOHYPGSERERMZj1JalqVOnYsSIERrrBAYGlk8w/1+LFi0AANevX0etWrXU1rGysoKVlVV5hkVERERGYtRkyd3dHe7u7sYMoZjC4QW8vb2NGwgRERFVCCbTZ+n27dtISUnB7du3UVBQICU1QUFBsLe3BwCEhIRgwYIF6NWrFwBI9e/duwcAiIuLAwB4eXnBy8sL8fHxWL9+Pbp16wY3NzdcuHABkydPxksvvYSGDRuW/0USERFRhWMyydLs2bOxevVqaT0sLAwAcOjQIbRr1w7A02QoLS1NqrNjxw6MHDlSWh84cCAAYM6cOfjoo49gaWmJ/fv3Y+nSpVAoFPD19UWfPn3wwQcflMMVERERkSmQCSGEsYMwdenp6XByckJaWhocHR2NHU6Fo1AopNa/jIwMzg1HREQVgra/vznOEhEREZEGTJaIiIiINGCyRERERKQBkyUiIiIiDZgsEREREWnAZImIiIhIAyZLRERERBowWSIiIiLSgMkSERERkQZMloiIiIg0YLJEREREpAGTJSIiIiINmCwRERERacBkiYiIiEgDJktEREREGjBZIiIiItKAyRIRERGRBkyWiIiIiDRgskRERESkgbmxA6DKz87ODkIIY4dBRERUKmxZIiIiItKAyRIRERGRBkyWiIiIiDRgskRERESkAZMlIiIiIg2YLBERERFpwGSJiIiISAMmS0REREQaMFkiIiIi0oDJEhEREZEGTJaIiIiINGCyRERERKQBkyUiIiIiDZgsEREREWnAZImIiIhIA3NjB1AZCCEAAOnp6UaOhIiIiLRV+Hu78Pd4SZgs6cGTJ08AAL6+vkaOhIiIiHT15MkTODk5lbhdJp6XTtFzKZVK3Lt3Dw4ODpDJZMYOp5j09HT4+vrizp07cHR0NHY4FQLvSXG8J8Xxnqji/SiO96Q4U7onQgg8efIEPj4+kMtL7pnEliU9kMvlqFGjhrHDeC5HR8cK/8Etb7wnxfGeFMd7oor3ozjek+JM5Z5oalEqxA7eRERERBowWSIiIiLSgMlSFWBlZYU5c+bAysrK2KFUGLwnxfGeFMd7oor3ozjek+Iq4z1hB28iIiIiDdiyRERERKQBkyUiIiIiDZgsEREREWnAZImIiIhIAyZLlcS3334Lf39/WFtbo0WLFjh9+nSJdVetWgWZTKayWFtbl2O0hvf3338jMjISPj4+kMlk2L59+3P3OXz4MJo0aQIrKysEBQVh1apVBo+zvOh6Pw4fPlzsMyKTyZCUlFQ+AZeDBQsWoFmzZnBwcICHhwd69uyJuLi45+63adMmhISEwNraGg0aNMCuXbvKIVrDK839qOzfJcuXL0fDhg2lwRVbtmyJv/76S+M+lfXzUUjXe1JZPiNMliqBjRs3YsqUKZgzZw7OnTuHRo0aISIiAvfv3y9xH0dHRyQmJkrLrVu3yjFiw1MoFGjUqBG+/fZbrerfvHkT3bt3R/v27RETE4N33nkHb7zxBvbs2WPgSMuHrvejUFxcnMrnxMPDw0ARlr8jR45gwoQJOHnyJPbt24e8vDx07twZCoWixH1OnDiB1157DaNGjcL58+fRs2dP9OzZE5cuXSrHyA2jNPcDqNzfJTVq1MDChQsRHR2Ns2fPokOHDujRowcuX76stn5l/nwU0vWeAJXkMyLI5DVv3lxMmDBBWi8oKBA+Pj5iwYIFauuvXLlSODk5lVN0xgdAbNu2TWOdd999V9SrV0+lbMCAASIiIsKAkRmHNvfj0KFDAoB4/PhxucRUEdy/f18AEEeOHCmxTv/+/UX37t1Vylq0aCHGjh1r6PDKnTb3o6p9lwghhIuLi/jpp5/UbqtKn4+iNN2TyvIZYcuSicvNzUV0dDQ6deoklcnlcnTq1AlRUVEl7peRkQE/Pz/4+vo+96+CqiAqKkrlHgJARESExntYFTRu3Bje3t54+eWXcfz4cWOHY1BpaWkAAFdX1xLrVKXPiTb3A6g63yUFBQXYsGEDFAoFWrZsqbZOVfp8ANrdE6ByfEaYLJm4hw8foqCgAJ6enirlnp6eJfYvCQ4Oxv/+9z/8/vvvWLt2LZRKJVq1aoX//vuvPEKukJKSktTew/T0dGRlZRkpKuPx9vbGihUrsGXLFmzZsgW+vr5o164dzp07Z+zQDEKpVOKdd95B69atUb9+/RLrlfQ5qUx9uQDt70dV+C65ePEi7O3tYWVlhXHjxmHbtm2oW7eu2rpV5fOhyz2pLJ8Rc2MHQOWvZcuWKn8FtGrVCqGhofj+++/x8ccfGzEyqiiCg4MRHBwsrbdq1Qrx8fH48ssv8csvvxgxMsOYMGECLl26hGPHjhk7lApB2/tRFb5LgoODERMTg7S0NGzevBnDhw/HkSNHSkwOqgJd7kll+YwwWTJx1apVg5mZGZKTk1XKk5OT4eXlpdUxLCwsEBYWhuvXrxsiRJPg5eWl9h46OjrCxsbGSFFVLM2bN6+UycTEiROxc+dO/P3336hRo4bGuiV9TrT9t2YKdLkfz6qM3yWWlpYICgoCADRt2hRnzpzBV199he+//75Y3arw+QB0uyfPMtXPCB/DmThLS0s0bdoUBw4ckMqUSiUOHDig8RlyUQUFBbh48SK8vb0NFWaF17JlS5V7CAD79u3T+h5WBTExMZXqMyKEwMSJE7Ft2zYcPHgQAQEBz92nMn9OSnM/nlUVvkuUSiVycnLUbqvMnw9NNN2TZ5nsZ8TYPcyp7DZs2CCsrKzEqlWrxJUrV8SYMWOEs7OzSEpKEkIIMXToUDFjxgyp/ty5c8WePXtEfHy8iI6OFgMHDhTW1tbi8uXLxroEvXvy5Ik4f/68OH/+vAAglixZIs6fPy9u3bolhBBixowZYujQoVL9GzduCFtbWzF9+nQRGxsrvv32W2FmZiZ2795trEvQK13vx5dffim2b98url27Ji5evCjefvttIZfLxf79+411CXo3fvx44eTkJA4fPiwSExOlJTMzU6rz7L+d48ePC3Nzc7F48WIRGxsr5syZIywsLMTFixeNcQl6VZr7Udm/S2bMmCGOHDkibt68KS5cuCBmzJghZDKZ2Lt3rxCian0+Cul6TyrLZ4TJUiWxbNkyUbNmTWFpaSmaN28uTp48KW1r27atGD58uLT+zjvvSHU9PT1Ft27dxLlz54wQteEUvvr+7FJ4H4YPHy7atm1bbJ/GjRsLS0tLERgYKFauXFnucRuKrvdj0aJFolatWsLa2lq4urqKdu3aiYMHDxoneANRdz8AqPx/f/bfjhBC/Pbbb6JOnTrC0tJS1KtXT/z555/lG7iBlOZ+VPbvktdff134+fkJS0tL4e7uLjp27CglBUJUrc9HIV3vSWX5jMiEEKL82rGIiIiITAv7LBERERFpwGSJiIiISAMmS0REREQaMFkiIiIi0oDJEhEREZEGTJaIiIiINGCyRERERKQBkyUiMhh/f38sXbq0xO0jRoxAz549yy0eTRISEiCTyRATE6PzvgcOHEBoaCgKCgr0H5gJe/jwITw8PExuhnmiZzFZIqrERowYAZlMJi1ubm7o0qULLly4YOzQjErfSdq7776LDz74AGZmZno7ZlmVJfnTl2rVqmHYsGGYM2eO0WIg0gcmS0SVXJcuXZCYmIjExEQcOHAA5ubmeOWVV4wdVqVx7NgxxMfHo0+fPsYOpVRyc3MNevyRI0di3bp1SElJMeh5iAyJyRJRJWdlZQUvLy94eXmhcePGmDFjBu7cuYMHDx5IdS5evIgOHTrAxsYGbm5uGDNmDDIyMqTthS0xixcvhre3N9zc3DBhwgTk5eVJde7fv4/IyEjY2NggICAA69at0zlWpVKJBQsWICAgADY2NmjUqBE2b94sbT98+DBkMhkOHDiA8PBw2NraolWrVoiLi1M5zieffAIPDw84ODjgjTfewIwZM9C4cWMAwEcffYTVq1fj999/l1rcDh8+LO1748YNtG/fHra2tmjUqBGioqI0xrxhwwa8/PLLsLa2lsri4+PRo0cPeHp6wt7eHs2aNcP+/ftV9vP398enn36K119/HQ4ODqhZsyZ++OEHlTonTpxA48aNYW1tjfDwcGzfvl2ltejx48cYPHgw3N3dYWNjg9q1a2PlypUAgICAAABAWFgYZDIZ2rVrB+D//l/Onz8fPj4+CA4OBqD9Z+DTTz+Fp6cnnJ2dMW/ePOTn52P69OlwdXVFjRo1pPMXqlevHnx8fLBt2zaN95GoImOyRFSFZGRkYO3atQgKCoKbmxsAQKFQICIiAi4uLjhz5gw2bdqE/fv3Y+LEiSr7Hjp0CPHx8Th06BBWr16NVatWYdWqVdL2ESNG4M6dOzh06BA2b96M7777Dvfv39cpvgULFmDNmjVYsWIFLl++jMmTJ2PIkCE4cuSISr1Zs2bhiy++wNmzZ2Fubo7XX39d2rZu3TrMnz8fixYtQnR0NGrWrInly5dL26dNm4b+/furtLi1atVK5djTpk1DTEwM6tSpg9deew35+fklxnz06FGEh4erlGVkZKBbt244cOAAzp8/jy5duiAyMhK3b99WqffFF18gPDwc58+fx5tvvonx48dLiV96ejoiIyPRoEEDnDt3Dh9//DHee+89lf0//PBDXLlyBX/99RdiY2OxfPlyVKtWDQBw+vRpAMD+/fuRmJiIrVu3SvsdOHAAcXFx2LdvH3bu3Kn1Z+DgwYO4d+8e/v77byxZsgRz5szBK6+8AhcXF5w6dQrjxo3D2LFji/VRat68OY4ePVriPSSq8Iw9ky8RGc7w4cOFmZmZsLOzE3Z2dgKA8Pb2FtHR0VKdH374Qbi4uIiMjAyp7M8//xRyuVwkJSVJx/Hz8xP5+flSnX79+okBAwYIIYSIi4sTAMTp06el7bGxsQKA+PLLLzXG16NHDyGEENnZ2cLW1lacOHFCpc6oUaPEa6+9JoQQ4tChQwKA2L9/v0qsAERWVpYQQogWLVqICRMmqByjdevWolGjRmrPW+jmzZsCgPjpp5+kssuXLwsAIjY2tsRrcHJyEmvWrClxe6F69eqJZcuWSet+fn5iyJAh0rpSqRQeHh5i+fLlQgghli9fLtzc3KTrEkKIH3/8UQAQ58+fF0IIERkZKUaOHKn2fIXXU1i30PDhw4Wnp6fIycmRynT5DBQUFEh1goODxYsvviit5+fnCzs7O/Hrr7+qnHPy5MmiXbt2Gu8PUUXGliWiSq59+/aIiYlBTEwMTp8+jYiICHTt2hW3bt0CAMTGxqJRo0aws7OT9mndujWUSqXK46169eqpdGD29vaWWo5iY2Nhbm6Opk2bSttDQkLg7OysdZzXr19HZmYmXn75Zdjb20vLmjVrEB8fr1K3YcOGKnEAkGKJi4tD8+bNVeo/u66JpmOrk5WVpfIIDnjasjRt2jSEhobC2dkZ9vb2iI2NLdayVPRcMpkMXl5eKtfRsGFDlWM/ex3jx4/Hhg0b0LhxY7z77rs4ceKEVtfYoEEDWFpaSuu6fAbk8v/7teHp6YkGDRpI62ZmZnBzcyt2v2xsbJCZmalVbEQVkbmxAyAiw7Kzs0NQUJC0/tNPP8HJyQk//vgjPvnkE62PY2FhobIuk8mgVCr1Fmdh/5g///wT1atXV9lmZWVVYiwymQwA9BaLrseuVq0aHj9+rFI2bdo07Nu3D4sXL0ZQUBBsbGzQt2/fYp2py3pPC5PeXbt2Yd++fejYsSMmTJiAxYsXa9yvaFKkC3XxanMNKSkpcHd3L9U5iSoCtiwRVTEymQxyuRxZWVkAgNDQUPzzzz9QKBRSnePHj0Mul0udf58nJCQE+fn5iI6Olsri4uKQmpqqdVx169aFlZUVbt++jaCgIJXF19dX6+MEBwfjzJkzKmXPrltaWuptTKSwsDBcuXJFpez48eMYMWIEevXqhQYNGsDLywsJCQk6HTc4OBgXL15ETk6OVPbsdQCAu7s7hg8fjrVr12Lp0qVSJ/HCliNtrlMfnwFNLl26hLCwsDIfh8hYmCwRVXI5OTlISkpCUlISYmNj8dZbbyEjIwORkZEAgMGDB8Pa2hrDhw/HpUuXcOjQIbz11lsYOnQoPD09tTpHcHAwunTpgrFjx+LUqVOIjo7GG2+8ARsbG63jdHBwwLRp0zB58mSsXr0a8fHxOHfuHJYtW4bVq1drfZy33noLP//8M1avXo1r167hk08+wYULF6RWIuDpm2gXLlxAXFwcHj58qPJWn64iIiJw7NgxlbLatWtj69atiImJwT///INBgwbp3PJVuM+YMWMQGxuLPXv2SC1Ghdcye/Zs/P7777h+/TouX76MnTt3IjQ0FADg4eEBGxsb7N69G8nJyUhLSyvxXPr4DJQkMzMT0dHR6Ny5c5mOQ2RMTJaIKrndu3fD29sb3t7eaNGihfS2U+Gr5La2ttizZw9SUlLQrFkz9O3bFx07dsQ333yj03lWrlwJHx8ftG3bFr1798aYMWPg4eGh0zE+/vhjfPjhh1iwYAFCQ0PRpUsX/Pnnn9Jr8NoYPHgwZs6ciWnTpqFJkya4efMmRowYodL3Z/To0QgODkZ4eDjc3d1x/PhxneJ89nyXL19W6duzZMkSuLi4oFWrVoiMjERERASaNGmi03EdHR3xxx9/ICYmBo0bN8asWbMwe/ZsAJCuxdLSEjNnzkTDhg3x0ksvwczMDBs2bAAAmJub4+uvv8b3338PHx8f9OjRo8Rz6eszoM7vv/+OmjVr4sUXXyzzsYiMRSaEEMYOgojIkF5++WV4eXnhl19+Mcjxp0+fjvT0dHz//fcGOX6hdevWYeTIkUhLS9Op1c6YXnjhBUyaNAmDBg0ydihEpcYO3kRUqWRmZmLFihWIiIiAmZkZfv31V+zfvx/79u0z2DlnzZqF7777DkqlUuVtsbJas2YNAgMDUb16dfzzzz9477330L9/f5NJlB4+fIjevXvjtddeM3YoRGXCliUiqlSysrIQGRmJ8+fPIzs7G8HBwfjggw/Qu3dvY4ems88++wzfffcdkpKS4O3tLY28bWtra+zQiKoUJktEREREGrCDNxEREZEGTJaIiIiINGCyRERERKQBkyUiIiIiDZgsEREREWnAZImIiIhIAyZLRERERBowWSIiIiLSgMkSERERkQb/D+vAGN/y3qyZAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "from scipy.interpolate import splrep, splev\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# step 1: read the data from the file\n",
    "data_pyscf = np.loadtxt('pyscf_fci_sto3g_h2.data')  #, skiprows=1)  # Assuming the first row is a header\n",
    "data_ibmq = np.loadtxt('h2_pec_ibmq_qasm_simulator.data')\n",
    "\n",
    "# extract the columns\n",
    "r = data_pyscf[:, 0]\n",
    "energy = data_pyscf[:, 2]\n",
    "\n",
    "r_ibmq = data_ibmq[:, 0]\n",
    "energy_ibmq = data_ibmq[:, 1]\n",
    "\n",
    "# standard deviation \n",
    "std_ibmq = np.sqrt((data_ibmq[:, 2]))\n",
    "\n",
    "# step 2: perform spline interpolation using SciPy\n",
    "spline = splrep(r, energy, s=0)  # Adjust s for smoothing factor\n",
    "\n",
    "# Evaluate the spline at more points for a smoother curve\n",
    "r_interp = np.linspace(r.min(), r.max(), 100)\n",
    "energy_interp = splev(r_interp, spline)\n",
    "\n",
    "# step 3: plot the interpolated data using Matplotlib\n",
    "plt.figure()\n",
    "plt.plot(r_interp, energy_interp, label='FCI/sto-3g', linewidth=2.5, color='gray')\n",
    "\n",
    "plt.errorbar(r_ibmq, energy_ibmq, std_ibmq, label='qc2 + ibmq_qasm_simulator',\n",
    "             barsabove=False, ecolor='black', elinewidth=1.5,\n",
    "             linestyle='None', marker='s', mfc='white',\n",
    "             mec='black', ms=10*0.75, mew=2*0.75)\n",
    "\n",
    "plt.xlabel('Bond length (angstrom)')\n",
    "plt.ylabel('Energy (hartree)')\n",
    "plt.title('Potential energy curve of H$_2$')\n",
    "plt.legend()\n",
    "plt.grid(False)\n",
    "plt.show()"
   ]
  }
 ],
 "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.11.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}