Simulating impact of T1 on VQE.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pylab\n",
    "from qiskit import Aer, IBMQ\n",
    "from qiskit.aqua import QuantumInstance, aqua_globals\n",
    "from qiskit.aqua.algorithms import VQE \n",
    "from qiskit.aqua.algorithms import NumPyEigensolver\n",
    "from qiskit.aqua.components.optimizers import SPSA\n",
    "from qiskit.aqua.components.variational_forms import RY\n",
    "from qiskit.aqua.operators import WeightedPauliOperator\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from qiskit import execute, QuantumCircuit, QuantumRegister, ClassicalRegister\n",
    "from qiskit.quantum_info.operators import Kraus, SuperOp\n",
    "from qiskit.providers.aer import QasmSimulator\n",
    "from qiskit.tools.visualization import plot_histogram\n",
    "\n",
    "# Qiskit Aer noise module imports\n",
    "from qiskit.providers.aer.noise import NoiseModel\n",
    "from qiskit.providers.aer.noise.errors import QuantumError, ReadoutError\n",
    "from qiskit.providers.aer.noise.errors import pauli_error\n",
    "from qiskit.providers.aer.noise.errors import depolarizing_error\n",
    "from qiskit.providers.aer.noise.errors import thermal_relaxation_error\n",
    "from qiskit.providers.aer import noise"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/anaconda3/envs/QPT2/lib/python3.8/site-packages/qiskit/providers/ibmq/ibmqfactory.py:192: UserWarning: Timestamps in IBMQ backend properties, jobs, and job results are all now in local time instead of UTC.\n",
      "  warnings.warn('Timestamps in IBMQ backend properties, jobs, and job results '\n"
     ]
    }
   ],
   "source": [
    "provider = IBMQ.load_account()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Define a problem to be solved with VQE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of qubits: 2\n"
     ]
    }
   ],
   "source": [
    "pauli_dict = {\n",
    "    'paulis': [{\"coeff\": {\"imag\": 0.0, \"real\": -1.052373245772859}, \"label\": \"II\"},\n",
    "              {\"coeff\": {\"imag\": 0.0, \"real\": 0.39793742484318045}, \"label\": \"ZI\"},\n",
    "              {\"coeff\": {\"imag\": 0.0, \"real\": -0.39793742484318045}, \"label\": \"IZ\"},\n",
    "              {\"coeff\": {\"imag\": 0.0, \"real\": -0.01128010425623538}, \"label\": \"ZZ\"},\n",
    "              {\"coeff\": {\"imag\": 0.0, \"real\": 0.18093119978423156}, \"label\": \"XX\"}\n",
    "              ]\n",
    "}\n",
    "\n",
    "qubit_op = WeightedPauliOperator.from_dict(pauli_dict)\n",
    "num_qubits = qubit_op.num_qubits\n",
    "print('Number of qubits: {}'.format(num_qubits))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Solve it classically for comparison"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Reference value: [-1.85727503-3.77606111e-17j]\n"
     ]
    }
   ],
   "source": [
    "ee = NumPyEigensolver(qubit_op.copy()) \n",
    "result = ee.run()\n",
    "ref = result['eigenvalues'] \n",
    "print('Reference value: {}'.format(ref))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Performance *with* noise\n",
    "\n",
    "Now we will add noise. \n",
    "\n",
    "Note: simulation with noise takes significantly longer than without noise."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "NoiseModel:\n",
      "  Basis gates: ['cx', 'id', 'u2', 'u3']\n",
      "  Instructions with noise: ['measure', 'u3', 'cx', 'reset', 'u2']\n",
      "  Qubits with noise: [0, 1, 2, 3]\n",
      "  Specific qubit errors: [('reset', [0]), ('reset', [1]), ('reset', [2]), ('reset', [3]), ('measure', [0]), ('measure', [1]), ('measure', [2]), ('measure', [3]), ('u2', [0]), ('u2', [1]), ('u2', [2]), ('u2', [3]), ('u3', [0]), ('u3', [1]), ('u3', [2]), ('u3', [3]), ('cx', [0, 0]), ('cx', [0, 1]), ('cx', [0, 2]), ('cx', [0, 3]), ('cx', [1, 0]), ('cx', [1, 1]), ('cx', [1, 2]), ('cx', [1, 3]), ('cx', [2, 0]), ('cx', [2, 1]), ('cx', [2, 2]), ('cx', [2, 3]), ('cx', [3, 0]), ('cx', [3, 1]), ('cx', [3, 2]), ('cx', [3, 3])]\n"
     ]
    }
   ],
   "source": [
    "# T1 and T2 values for qubits 0-3\n",
    "T1s = np.random.normal(50e3, 10e3, 4) # Sampled from normal distribution mean 50 microsec\n",
    "T2s = np.random.normal(70e3, 10e3, 4)  # Sampled from normal distribution mean 50 microsec\n",
    "\n",
    "# Truncate random T2s <= T1s\n",
    "T2s = np.array([min(T2s[j], 2 * T1s[j]) for j in range(4)])\n",
    "\n",
    "# Instruction times (in nanoseconds)\n",
    "time_u1 = 0   # virtual gate\n",
    "time_u2 = 50  # (single X90 pulse)\n",
    "time_u3 = 100 # (two X90 pulses)\n",
    "time_cx = 300\n",
    "time_reset = 1000  # 1 microsecond\n",
    "time_measure = 1000 # 1 microsecond\n",
    "\n",
    "# QuantumError objects\n",
    "errors_reset = [thermal_relaxation_error(t1, t2, time_reset)\n",
    "                for t1, t2 in zip(T1s, T2s)]\n",
    "errors_measure = [thermal_relaxation_error(t1, t2, time_measure)\n",
    "                  for t1, t2 in zip(T1s, T2s)]\n",
    "errors_u1  = [thermal_relaxation_error(t1, t2, time_u1)\n",
    "              for t1, t2 in zip(T1s, T2s)]\n",
    "errors_u2  = [thermal_relaxation_error(t1, t2, time_u2)\n",
    "              for t1, t2 in zip(T1s, T2s)]\n",
    "errors_u3  = [thermal_relaxation_error(t1, t2, time_u3)\n",
    "              for t1, t2 in zip(T1s, T2s)]\n",
    "errors_cx = [[thermal_relaxation_error(t1a, t2a, time_cx).expand(\n",
    "             thermal_relaxation_error(t1b, t2b, time_cx))\n",
    "              for t1a, t2a in zip(T1s, T2s)]\n",
    "               for t1b, t2b in zip(T1s, T2s)]\n",
    "\n",
    "# Add errors to noise model\n",
    "noise_thermal = NoiseModel()\n",
    "for j in range(4):\n",
    "    noise_thermal.add_quantum_error(errors_reset[j], \"reset\", [j])\n",
    "    noise_thermal.add_quantum_error(errors_measure[j], \"measure\", [j])\n",
    "    noise_thermal.add_quantum_error(errors_u1[j], \"u1\", [j])\n",
    "    noise_thermal.add_quantum_error(errors_u2[j], \"u2\", [j])\n",
    "    noise_thermal.add_quantum_error(errors_u3[j], \"u3\", [j])\n",
    "    for k in range(4):\n",
    "        noise_thermal.add_quantum_error(errors_cx[j][k], \"cx\", [j, k])\n",
    "\n",
    "print(noise_thermal)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "VQE on Aer qasm simulator (with noise): (-1.8193398600806043+0j)\n",
      "Delta from reference: [0.03793517+3.77606111e-17j]\n"
     ]
    }
   ],
   "source": [
    "backend = Aer.get_backend('qasm_simulator')\n",
    "quantum_instance = QuantumInstance(backend=backend, seed_simulator=167, seed_transpiler=167,\n",
    "                                   noise_model=noise_thermal)\n",
    "\n",
    "counts1 = []\n",
    "values1 = []\n",
    "def store_intermediate_result1(eval_count, parameters, mean, std):\n",
    "    counts1.append(eval_count)\n",
    "    values1.append(mean)\n",
    "    \n",
    "\n",
    "aqua_globals.random_seed = 167\n",
    "optimizer = SPSA(maxiter=200) \n",
    "from qiskit.circuit.library import RealAmplitudes\n",
    "var_form = RealAmplitudes(num_qubits)\n",
    "vqe = VQE(qubit_op, var_form, optimizer, callback=store_intermediate_result1)\n",
    "vqe_result1 = vqe.run(quantum_instance)\n",
    "print('VQE on Aer qasm simulator (with noise): {}'.format(vqe_result1['eigenvalue'])) #'energy' deprecated\n",
    "print('Delta from reference: {}'.format(vqe_result1['eigenvalue']-ref)) #'energy' deprecated"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/anaconda3/envs/QPT2/lib/python3.8/site-packages/numpy/core/_asarray.py:85: ComplexWarning: Casting complex values to real discards the imaginary part\n",
      "  return array(a, dtype, copy=False, order=order)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.lines.Line2D at 0x7ffa60a8abe0>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmkAAAFNCAYAAABbpPhvAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nOzdeXxcdb3/8dd31uxJ27TpTlvK0rJrZUcQEFBkEcXrclF/ouj14n7V671eRe91A5ergrK5IIKKV1AQFARKKVuhQCnd931L2mZPZv3+/jjnTGYmk8xkmazv5+ORB8nMyZlv0tC+8/l8F2OtRURERERGFt9wD0BEREREulNIExERERmBFNJERERERiCFNBEREZERSCFNREREZARSSBMREREZgRTSRESGiTFmtjGm1Rjj7+Uaa4yZPwRjudUY81/Ffh0RKZxCmohkMMa83xiz3A0Pe40xfzPGnD3c4xqLrLU7rLUV1toEgDHmKWPMR4dpLJ+w1v73cLy2iOSmkCYiKcaYzwP/C3wbqANmAz8DrhjOcaUzxgSGewwiIkNBIU1EADDGVAPfBP7VWnu/tbbNWhuz1j5krf2ie03YGPO/xpg97tv/GmPC7nPnGWN2GWO+YIw54Fbh/p/73GnGmH3pbT1jzDuNMSvd933GmH83xmw2xhw0xtxnjJnoPjfHbflda4zZATxpjPEbY35gjGkwxmw1xlzvXhPwvhZjzC/cMew2xvyP99rGmA8bY54xxnzfGHPY/fy3pY1rojHmV+7Xd9gY8+e0595hjFlhjGk0xjxnjDmxh+/lN4wxP3XfDxpj2owxN7kflxpjOt3X8b62gDHmW8A5wM1uFfPmtFteaIzZ6L7uLcYY08Pr3uB+735jjGkxxqw2xixKe36BW61rdJ+7PO25Xxtj/sd9v9YY81f3ukPGmKXGGJ/73HRjzJ+MMfXu9+7Tvf1ciUj/KaSJiOcMoAR4oJdr/hM4HTgZOAk4Ffhq2vNTgWpgBnAtcIsxZoK1dhnQBpyfdu37gXvd9z8FXAmcC0wHDgO3ZL32ucAC4GLgY8Db3HG8wf3cdL8G4sB84BTgIiC9jXgasB6oBW4EfpEWfO4GyoDjgCnAjwCMMacAvwQ+DkwCbgMe9EJqliXAee77bwL2AW92Pz4DWG+tPZT+Cdba/wSWAte7LdDr055+h3ufE4H3uN+DnlwO/B6oAR4EbnbHHwQeAh5zv65PAfcYY47JcY8vALuAyTgV1f8ArBvUHgJew/kzvgD4rDGmt/GISD8ppImIZxLQYK2N93LNB4BvWmsPWGvrgW8A16Q9H3Ofj1lrHwFaAS8E/A54H4AxphJ4u/sYwCeA/7TW7rLWRoAbgHdntTZvcKt7HThB5cfu9YeB73oXGWPq3Ht/1r3+AE7Qem/avbZba+9w54LdBUwD6owx03DC3yestYfdr2OJ+znXAbdZa5dZaxPW2ruACE5ozfY8cJQxZhJOOPsFMMMYU4ETNpfk+JzefNda22it3QEsxgmnPXnGWvuI+7XdjROmccdZ4d4raq19Evgr7p9JlhjO9+QI93uw1DoHPb8JmGyt/aZ7jy3AHWR+b0VkkCikiYjnIFCbZ87XdGB72sfb3cdS98gKee04wQCcqtlVbuXpKuAVa613ryOAB9z2WiOwFkjgVHE8O7PGsbOH544AgsDetPvdhlM98uzz3rHWtrvvVgCzgENu8Mt2BPAF757ufWdlff3ePTuA5TiB7M04oew54Cz6F9L2pb2f/j0t5NoS9890OrDTWptMe347TkUs203AJuAxY8wWY8y/u48fAUzP+h78B5l/TiIySDQBV0Q8z+NUhq4E/q+Ha/bg/EO92v14tvtYXtbaNcaY7TiVqvRWJzgh6yPW2mezP88YM8e7RdrDe4GZaR/PyrpXBKjNUxXMZScw0RhTY61tzPHct6y13yrwXktw2runAC+5H1+M0yJ+uofPsT08Phj2ALOMMb60oDYb2NBtENa24LQ8v2CMOR5nHuBLON+Drdbao4o4ThFxqZImIgBYa5uAr+HMI7vSGFPmTnp/mzHmRvey3wFfNcZMNsbUutf/tg8vcy/wGZzq0h/THr8V+JYx5ggA9/69rSi9D/iMMWaGMaYG+HLa17EXZ97VD4wxVcZZlHCkMebcfINzP/dvwM+MMRPcr9+bS3YH8AnjLIIwxphyY8ylbus2lyXAB4E11too8BTOvLitbqs4l/3AvHzj7KdlOJW1L7lf13nAZTjz1zK4CyTmu/P0mnCqmkngRaDFGPNldwGE3xhzvDHmTUUas8i4ppAmIinW2h8An8dZDFCPUzm5HvBWOP4PThtvJfA68Ir7WKF+h9Pue9Ja25D2+I9xJrk/ZoxpAV7AmdzfkztwgthK4FXgEZyFAgn3+Q8CIWANziKE/8OZY1WIa3DmZK0DDgCfBbDWLsdZsHCze89NwId7uc9zQCldVbM1QCc9V9HA+T68211V+pMCx1sQNyhehlPJbMDZWuWD1tp1OS4/CngcZ07h88DPrLWL3Xlu78CZE7fVvc+dOItFRGSQGWcuqIjI6GWcLTRutdYeMdxjEREZLKqkicio47ba3u7uLzYD+Dq9bx0iIjLqqJImIqOOMaYMZ87XsUAH8DDwGWtt87AOTERkECmkiYiIiIxAaneKiIiIjEAKaSIiIiIj0JjczLa2ttbOmTNnuIchIiIiktfLL7/cYK2dnP34mAxpc+bMYfny5cM9DBEREZG83NNYulG7U0RERGQEUkgTERERGYEU0kRERERGIIU0ERERkRFIIU1ERERkBFJIExERERmBFNJERERERiCFNBEREZERSCFNREREZARSSBuAg60RVu1uGu5hiIiIyBikkDYAdyzdyod++eJwD0NERETGIIW0AWiLxGmNxId7GCIiIjIGKaQNQDyZJJ60wz0MERERGYMU0gYgnrAkkhZrFdRERERkcCmkDYBXRYslFNJERERkcCmkDUAskQSctqeIiIjIYFJIG4CEKmkiIiJSJAppA+CFs3hiaCppq3Y3ceUtz9IRTQzJ64mIiMjwUUgbgITb5kwM0QrP1XuaWLGzkf3NnUPyeiIiIjJ8FNIGILVwYIhCmvd6mgMnIiIy9imkDUBq4cAQtTu9il00rjlwIiIiY51C2gAM1sKB5zY38NU/v96H11MlTUREZKxTSBuA1MKBAbYfn1x7gHuX7ch7XULtThERkXFDIW0AUqFpgJW0tmicpIVknrltcbU7RURExg2FtAHw2o4DbT+2RZwtNRJ5jpdSu1NERGT8UEgbgK724wAraZF4xv3yv55CmoiIyFgXGO4BjGbxHJUtay3GmD7dpy0az7hft9dJJElYq3aniIjIOKJK2gB0bcHhhKb2aJy5X3mEW5ds7nbt+n0tLPivv7O7saPbc6l2Zw9z2256bD3vv2NZas6a2p0iIiJjn0LaAGS3H/e4Aey+l3Z2u3bbwTY6Ygn25gppbiUt1kMbc09jJ3saO3JW7kRERGRsUkgbAG8LDu+/TR0xAKpKg92ujcadYBXNEbDyzUmLxZPEkzZ1DNVAV5OKiIjIyKeQNgDZoelgaxSA6l5CWq6A1e62O3uakxZLJEkk0+akDbCSZq1lX5PO/xQRERnJFNL64e4XtvOXFbtTgctrdzb0FtJ62K7DWptqd/Y0Jy2aSBJPJAdtTtozmxo463tPKqiJiIiMYFrd2Q+/W7aD6TUlqTlksVQlLQJAVWmADftbCPl91FWVcKg9mqqkZR8h1RlL4hXQetpaI5Zw2p3xtM1ztza0MXtiGX5f31aSAtS3REgkLQ2tEaZWl/T580VERKT4VEnrh1DARySeTDtxwKukOSHNYPj3P63kW4+s5c6lW7ji5mfTQlpmEGt156NBL3PSEtadk+Y8v7epkwt/uITHVu/r1/i9MUTiWoAgIiIyUqmS1g+hgI9oPNm1cMANT167M5ZI0hZJEA7EaWiNcLAtkmp3ZlfL2qNdIS3fnDQvpNW3OpWwejcUFurXz27lqQ31XLCgDoBIPNGnzxcREZGho0paP4QDPjpjXQHHq6R5oSmaSBJNJOmMJ+iMJbG2q2IWy9qI1tsjDboqacmk5X23v8Dja/Y794tnLhzwVoN2RPsWsl7f3cxrOxuJxVVJExERGekU0voh5PfRHk0PaV4lLZL6OBpPEoklU9Wqlk5ne47svdDaclTSOmIJnt9ykJd3HM543GuZtnbGU9f1RWcs4bZO3ZDWx88XERGRoaOQ1g+hQGZI84JXQ4sT0mJuJS0ST6SqVc0dXiUtK6SlzUnzKnLevb3nuuaQuYEv0r+Q1hFLZLRpI/EkG/e3aHNcERGREUghrR9CAV9mBSxhicQTNHd2hapYIkkknkyFNK+Slj3vLL3dmaqkuSHNq5hltydbI869OvvY7uyIJtzw6LZnWyK87cdLeeCV3X26j4iIiBSfFg70QziQ3e5MpjayBYi67c5OXzI1d63FDVzZG9Gmhz1vTlp7zHnMm8cW9SpfsQG2O91KnDemAy0R4kmb8zxRERERGV7DUkkzxlxtjFltjEkaYxb1ct0lxpj1xphNxph/H8ox9sZb3emJJS2N7bHUx/FUJS2RVkmLu89lV9K6z0lLtTujududral2Z9/alF6FzltR2tjuBMvD7dEeP0dERESGx3C1O1cBVwFP93SBMcYP3AK8DVgIvM8Ys3Bohte7kN+f8XE8kcyoakXceV9OuzNr4UAiewuO9NWd7py0SFa7M2tfM29OWV9Xd3oVNO/+h91geahNIU1ERGSkGZaQZq1da61dn+eyU4FN1tot1too8HvgiuKPLr9QIPPbFkvY1ErJkmBXKzQaT9LpVru65qtlVtJaI5lz26Cr0tXabeFAZsDr7KHd2RFN8J1H1mbcG7rao16FrskNaaqkiYiIjDwjeeHADGBn2se73MeGXXZIiyeTqflelSXBjA1qmzucIOQFpnh2JS3HiQNemGqNxLHWdq3GzAplPYW0Zzc1cNvTW3huU0PG413tTue/jR1OODvUFmNbQxsHmnWWp4iIyEhRtJBmjHncGLMqx1tRqmHGmOuMMcuNMcvr6+uL8RIp4eyQlrCpilllSSBjxWZTRyzj2u7HQnVf3dm1BUcio/KWXUnraeHA1oY2oHuFrDOeucWHN4/ucFuU6+5ezrceWZvzfiIiIjL0ira601p74QBvsRuYlfbxTPexnl7vduB2gEWLFuU+X2mQhPzd251elaqyJMi+pq6KVHawiiWz250xQn4f0UTXWaBeiGqNxDOObio4pB10QtrBtLlmiaRNLXbwFis0ps1JO9gWobYi3OPXLCIiIkNrJLc7XwKOMsbMNcaEgPcCDw7zmIA87c5wIGMxQLZoPMl//XkVK3c1Ak6lbWJ5COiqsnWktUvTK3HZZ212RhN8/9H1vOOnS/n5U5tTj29zK2mH0rYFSW+NpubMua8XTTgLHZo7Yzy1/gDf/du6Xr9+ERERKb7h2oLjncaYXcAZwMPGmEfdx6cbYx4BsNbGgeuBR4G1wH3W2tXDMd5s3UNaZruzN00dMe5+YXvqXM7mjjgT3JCWXUkDMrb2yF500BFL8PuXdrJqdzN/emVX6nGv3Zm+arMjI6RlLijwNHfEeXjlXm5dsrlbm1ZERESG1rBsZmutfQB4IMfje4C3p338CPDIEA6tINntzniia9PaqpJgr5/rbavhbX/R3BljzqRy5z45Qlpv22O0RxMkks79dh1ux1onLO51262H2nNX0tLnzKVr7oylxrVyVyPnHDW5169FREREimcktztHrG6VNHcLDmOgLOzv4bMc3irPRrdSld7u3N3YwWU/fYZt7pwy6H17jEg8STxpmVFTSmcsycG2aMbnpge89JDW01y25o5YaoPbFTsae/06REREpLh0LFQ/dNsnLWnpiCUoCfi7VdmypUJae5Rk0tIaiadC2m9f2J46mcBzuICNZo+fUcXuxg52He5g56F2AI6uq+Bga5T33v48V79xFkfXVea9T9KSOiLqtV0KaSIiIsNJlbR+SA9p4YDPbXcmKQn6COYJaV4Ia2yP0dIZx1pSIa0y3JWZfcb57+H2/HPDjp9eDTgtz8fX7qemLMgZ8yaxu7GDF7Yc4hfPbO2xelYSzBzvPnevtBU7G7G2qItkRUREpBcKaf2Qvk9aacjv7pOWoCToJ+A3vX5ua6Rrl/9m96goL6SVpYW0Se52GI0FnAZw/AwnpG3c38o/1uznbcdPY0pVSer5NXubWb2nKefnTqksIeAzTChz5tJZCxXhAA2t0YwtPERERGRoKaT1Q0ZIC/qJJZN0xpOUBv15K2neKtCm9lhqBeXErNWdAFMqnZBWSCVt5oRSasqC/Pq5bbRHE1x20jQmufcEMAbuW74r5+eWhfx8/+qT+Nxbj049duxUpzWavt+biIiIDC2FtH5IP2C9JNhVSQsH889J87RE4qlK1YQyJ1Clb40xORXS8lezaspC1FWW0NQRY25tOafNnZQKfrMnlnHs1CrW7m3O+bmlIT9XnjKDk2fVpB47dpoT0vYqpImIiAwbhbR+SJ+TVhL0E3O34HDmpPXe7ky3w53kX10axO8zGVtvTCgLEfCZHkNa+vy1mrJg6hzOr1+2EL/PMKnCCWkLplUyo6Yk5z3AqQSCc1KCZ8G0KgD2NXUU/LWIiIjI4NLqzn4IZbQ7fTR3xp2QFvATKLCSBrDD3S6juqx7SCsN+akoCXC4LXe7s7osSEskTkU4QNDv46Z3n8TKXY2cd8wUoKs6t2BaFQdbe67GeSGtKm0T3iMnVxDwGVXSREREhpEqaf2QXUlLX92Zq91ZFsq9d5pXSasqCRD0mYw5aWVBPxXhQI+b2XohrMad8P/moydz/flHpZ4/YlI5Hz5zDu88ZQZTq3uupJWEulfSJpWHqKsq0Zw0ERGRYaSQ1g/pQaw06CfmzkkrDfkJBrq3O6tLc59CsP1gOz7jrKb0+zI/rywcoCIc6HHrDC+cef/N5vcZbrj8OI6YVM60XkKaV0kLBXyp96vLgkyrLlElTUREZBgppPVDt0qae8B6SSBzdadxc5d3VFR5VkVtx6F2qkqDGGO6tUnLQn6qcoQ7LyB6wc+rqPVmWnUpkLkq1QuF6VW+qlKn5VlTGmJqdQl7NSdNRERk2Cik9UN62AkHfe7qziThoJ+Ar+u5ipATerzwk74PGjhnb3phq1slLeTPmCfm8Taf9SpoPVXp0nmVtLKQn4D7OhXuWLzqGThhsjzkJxTwpSppyaTlmw+tYdOB1ryvIyIiIoNHIa0furc7k3RGndWdIbfdaQyUu0HIC1JeMPL7TCoseVW2QLeQFsiYJ+YpcUNVTalTQSukkubNSSsN+lNVQG8sJekhrTRIjXu/qdWlROJJXt15mF8+u5WP37087+uIiIjI4FFI6wdfWshy2p3WaXembWYb9PtSVa9Uu9M9fL0k4EttVlvrbpXhVdKmVIb56qULeOvCOipzVtLceWOlTtWrt0UB6Z8zoSxISahrfN69S9PanXVV4VTVra7KGd+uwx0Z4xMREZGhoS04+ikU8JGMJQj6fUTiSRJJS2lauzPs9xEO+PGlVdTK3PZnOOjn1x85lfX7WlKbyHqhLxz08dFz5gHkDGlee7Ik6OPBT53d66KAdNOqSzGGVEjL1e78xuXHE086JyKUu2NtdE88yD5UXkRERIpLIa2fwgEf8aQl6O/aOiO93RkM+AgHnaDmVdRSLcaAj6PrKjm6rjJ1P69Slb7wIHe70w2BAT9HTq4oeLxnH1VLWyTO4nUHnLGUdA9p3ikHzus4j3tHVxV6koKIiIgMDoW0fgoFnApa+kKBzHanocQNaOGAE3i8lZThYPd907zPC2WEtJ7bneFg30LTf7x9AQDn3rQYSAuMPezh5oXB5g5V0kRERIaD/uXtp1DAR8BnCKQdA5W+BUcorZIWzpqsH84ReLxKWvpz3ly2UMCHNyUsFdL6GZq6zUnLERiha65ac6cX0nJfJyIiIsWhSlo/hfw+gv7MszrDaWd3Bv0+wm5Q84JV+py0bN6ctPSKlRekQn4fWIgmkqlQlesehfBC2htmT6C2IsyZR07KeV2JG8qaO+JdYxAREZEho5DWT6GAH78vngpe4FSl0tuWR9VVYoyhsiSAMV1bcZT0UknLDGnO9UG/IWkN0UTXPSrD/fuj8+5fFgrwhYuO6fG67Epafyt3IiIi0j8Kaf0UCjiVtHOPnpx6LH1OWijg48uXHAtAezTO7IllbDvonNWZu5LWfU6at5lt0O8sUgBYNGcCFx9fxxuPmNC/cacqfb1vqeFV0rS6U0REZHjoX95+Cvt9+H2GWRPLUo+VZFXSPGWhAGfOr00Fo75X0nypdmjAbzj/2DqM6d++ZamFDXlCV0nIeV6rO0VERIaH/uXtp1DAl1o08P7TZgNOgEqfk5bNe6wkVyXN74W0rudSc9ICPvxupc3vG9gfWSqk5blPyO/DGDjcHk2NQURERIaO2p39FAr4UkHnG5cfx7lHT+aUWTW0RxOp57N5Aam31Z2ZFTg/fp8T/Lzwl318VF+lbxHSG2MMpUF/r19Ptj+8tIOF06o5YWb1gMYoIiIiCmn9duUpM2hoiQBO8Ln4uKmp99P/m86rluWspHlz0gJd4ckYQ0U4QNBtrQL4+tnm9KRvtptPSVpIKyQcfvlPrwOw7buXDmCEIiIiAgpp/Xb5SdNzPu5VqHJVy0K9VNICOSppAFWlAQLpc9IGWEkLFdjuhMw91JLWDuh1RUREpG800WiQGeMcvp6rndhbJc3v775wAKAyHCTkN6lKmj9PmzKfroUD+e+TfqpBUhlNRERkSKmSVgRBvy/nHK7U4eu9VdKynju6roKA30dLZzzjun6PLdBzOzabKmkiIiLDRyGtCJyJ/jnanYFeKmmpdmfmcz/6p5MxxnDpT5Y61w10TlqOLUJ6kj5OZTQREZGhpZBWBP/6lvmcPKum2+OpSlqOw9F7qqR5+6F5z/sHvLqza7+1fFRJExERGT4KaUXw8XOPzPl4ap+0HIeV+1OrO3NXuPy+wsNVb3pbfZqtJGNOWuEhLZm0+AYYJkVERMY7LRwYQhXueZveJrXpeqqkpZ53Q9VAt+DoW0jrCpOJZOGv0RqN93lcIiIikkkhbQjNnlTGXR85lQsX1nV7zquQhXsIT11bcAzsj8w7OL2QA9Mz56QVXklrcs/7FBERkf5Tu3OIpR/Ins4LYT1tjeEfpDlpV71hBjNqSnMuXsjWlzlp6SGuqSPGrP4PUURERFAlbcRIzUnz5w5PgUGakzalsoTLetiIN1tJH/ZJS3++uVOVNBERkYFSSBshhmpOWl/0pZKW/nxzh0KaiIjIQCmkjRD+fCFtkI6F6otwH/ZJywxpWjggIiIyUAppI0RPZ3d6BmtOWl+UZqzuzFNJS1v92aRKmoiIyIAppI0QPZ3d6RmsOWl9UdLPdqdCmoiIyMAppI0QXgjraWsMb07aQI+F6ovSUNdY+tLubOqIEY33YWM1ERER6UYhbYQI5DlxYLCOheqL9JMR8lfSut6/+4XtfOTXLxVrWCIiIuOCQtoI4bUxezoJwD9Im9n2RUmof/ukATyzqaEoYxIRERkvFNJGiNTh63kqaUOY0bIqab1f6z1//Vvm84HTZuMzzhmeIiIi0j/DEtKMMVcbY1YbY5LGmEU9XDPLGLPYGLPGvfYzQz3OoXThwil8+ZJjmVZdkvN5b07aUFbSSkOFHwvlrf6sqwozf0oFSQuNWkAgIiLSb8NVSVsFXAU83cs1ceAL1tqFwOnAvxpjFg7F4IbDlMoS/uW8IzE9LAwYjjlpC6dV8ZGz5jJ7YlneLTi8EOfzGSaWhwA41BYB4HBblH+67Xk2HWgt7oBFRETGkGEJadbatdba9Xmu2WutfcV9vwVYC8wYivGNRP5h2Mw2FPDxtcsWMqE8VHC702cMtRVhABpaowD8/qWdLNt6iN++sL2YwxURERlTRsWcNGPMHOAUYNnwjmT4TKoIU10axDeEIc3jM4Xvk+YzpFXSnJC2clcjALMnlhVxlCIiImNLoFg3NsY8DkzN8dR/Wmv/0of7VAB/Aj5rrW3u5brrgOsAZs+e3cfRjnz/fPps3nHitGF5bZ8xBe+TZoxhUoUT0g62Ou3OlbuaAIgmtHeaiIhIoYoW0qy1Fw70HsaYIE5Au8dae3+e17sduB1g0aJFY25ZYTjgp67Kn//CIiikkmbT2p0TytyQ1haloTXC7sYOANqjiaKOU0REZCwZse1O48yg/wWw1lr7w+Eez3hmjMkb0ryFBT7j7PVWUxbkYGuU13Y2pq5pj+jgdRERkUIN1xYc7zTG7ALOAB42xjzqPj7dGPOIe9lZwDXA+caYFe7b24djvOOds+dZ79d0zUlz5sxNLA9xqC3KjkPtABgD7TFV0kRERApVtHZnb6y1DwAP5Hh8D/B29/1ngKGfJS/d+H2GeJ75ZKnVne7ChtryMA2tEXYd7qA06Ke2MkSH2p0iIiIFG7HtThk5fHnanUs31hNzQ5y3+NSrpO081M6siaWUhwK0qd0pIiJSsGGppMno4sxJy/3c0o31XPOLF7n0BGflqdfunFQR4sVtUYJ+HzMnlHG4PUqH2p0iIiIFU0iTvHym52Oh9rgrN725Z14lbW5tOYfaorR2xnnTnAlE40lV0kRERPpA7U7Jy9dLJS2acJ4I+J105h1rdeaRte7zSWZOKKM05NcWHCIiIn2gkCZ59bZPWizuzEXzu+HMa3ceO7WSSe7JAzMnlFIW8qvdKSIi0gcKaZKXMabHA9a9UwT8qQPgncd9PsMZR04CYNbEMspCAdoiCmkiIiKFUkiTvPy9HAsVdStpQTedee1OgHecOJ2J5SHm1JY7lbSo5qSJiIgUSgsHJC+fr5d2Z1YlzZcW0i45fioXH1eHMYbykJ/2WAJrbUaQExERkdxUSZO8ejsWymt3Bv1eSOv+uQCloQDWQmdMh6yLiIgUQiFN8vIV0O7MVUlLVxZyDodvV8tTRESkIAppkldvqzu9kOaFs546mV0hTYsHRERECqGQJnn5jCGRZ05azN0vzd9jJc2Z/qiQJiIiUqiNc9cAACAASURBVBiFNMnLZwzJHqaSeeEs7l7gy56U5ioLq90pIiLSFwppkldvx0J57c64G9Z6yGiUBdXuFBER6QuFNMmrt2OhIl5IcytpPW2vUR5Wu1NERKQvFNIkr0L2SeuqpOUOaaVa3SkiItInCmmSl+mlkpZaOJDM0+7U6k4REZE+UUiTvAqZkxbL2oojm1Z3ioiI9E1BIc0Y8wNjzHHFHoyMTL1twZE9J62nkFbuVtLaImp3ioiIFKLQStpa4HZjzDJjzCeMMdXFHJSMLM4WHLlDWkfMqYzFvXZnDz9RAb+P0qCfls5YUcYoIiIy1hQU0qy1d1przwI+CMwBVhpj7jXGvKWYg5ORobdjoTq9kJZn4QBAZUmAlk5V0kRERApR8Jw0Y4wfONZ9awBeAz5vjPl9kcYmI0Rvx0J5B6bHE8nUtT2pUEgTEREpWKCQi4wxPwIuA54Avm2tfdF96nvGmPXFGpyMDD5fL/ukuZU0b3VnT/ukAVSWBGnRnDQREZGCFBTSgJXAV621bTmeO3UQxyMjkOmlkubNSfO24ujp7E6AqpKA5qSJiIgUqNCQ9hpwTFaVpAnYbq1tGvRRyYjinDjQPaTFE8nUgoFEgXPS9jZ1FmeQIiIiY0yhIe1nwBtwKmoGOB5YDVQbY/7FWvtYkcYnI4C/h81sO+Ndp67HUsdC9XyfynBQlTQREZECFbpwYA9wirV2kbX2jcApwBbgrcCNxRqcjAw9LRzwVnZC2urOXlYOeAsHPn/fCj72m+WDP1AREZExpNBK2tHW2tXeB9baNcaYY621W3qbKC5jg3G34LDWZiwMyAhpeY6FAqfd2R5NcP8ruwE40NzJlKqS4gxaRERklCu0krbGGPNzY8y57tvP3MfCgPpXY5w3zyy7mJYe0rKvzaWyJJjx8R9f3jXwwYmIiIxRhYa0DwGbgM+6b1uAD+MENG1oO8Z51bHslqe3R1rmtb0vHEj34Io9Ax+ciIjIGJW33eluYvuItfYtwA9yXNI66KOSEcWbZ5a9eCB3Ja3n+1SlhbSJ5SH2NnUMyvhERETGoryVNGttAkjqvM7xywxSJa0i3NXuXDitiubOeGp/NREREclU6MKBVuB1Y8w/gNSGttbaTxdlVDKieBvUdg9pfZ2T1vXjtmBaJc9sauBwW1SLB0RERHIoNKTd777JOOQzududuapgppfabGZIqwLgoEKaiIhITgWFNGvtXcaYUmC2tVZndY4zPbU7vfM6g35DrKATB5x2Z01ZkGnVpQAcaosO9nBFRETGhIJWdxpjLgNWAH93Pz7ZGPNgMQcmI0dqC46swlncraSFA/7UY72d3elV0uoqS5hUEQKcSlq63Y0dROOapyYiIlLoFhw34Byk3ghgrV0BzCvSmGSE6WkLDq/dWRLs+jHqbW/jkqCfkN9HXXUJE8udkHaoNZJxv4t/9DR3PrNlkEYuIiIyehUa0mI5DlJXuWOc8LbgSLghbe3eZu5bvjPV4gz5u36Memt3grP1xoyaUiaUhTAms915uD1KayTOi1sPDfaXICIiMuoUunBgtTHm/YDfGHMU8GngueINS0YSX9bqzrf9eCkAN1y2EIBw0J92be/3+sWHFzG5MozfZ5hQFspodx5ucw6veG1nY7cjqERERMabQitpnwKOAyLA74BmnJMHZBzo6Vgo77zOcKDwStpx06uZUums5pxYHmL1nmZ+8sRGkknLwTan9Xm4PcbOQ9roVkRExrdCV3e2A//pvsk409OctKi3cCCtktaX4tfE8hAvbj3Eip2NnH/slFQlDeC1XY3MnlTW/0GLiIiMcoWu7jzaGHO7MeYxY8yT3luxBycjQ0/7pHknDniVNJ+hTy3KqrQD11/f3cSh9q7W58pdjRnXHmqL8qFfvsjBtIUGIiIiY1mhc9L+CNwK3Al032ZexrTUPmlZKa0jGgecfdIgf6sz26rdXWtRXt/dRJ3bBp01sZS9TZ0Z167d28ySDfWs29fCWfPDfXodERGR0ajQOWlxa+3PrbUvWmtf9t76+6LGmKuNMauNMUljzKI81/qNMa8aY/7a39eTgUlfONAR7cro7dEEIb8v9XxfQ9rnLzqaCWVBFh0xgdW7mzjcHqWyJEBtRZimjljGtd78t+yWq4iIyFhVaEh7yBjzSWPMNGPMRO9tAK+7CrgKeLqAaz8DrB3Aa8kA+dyfkqSFhrR2Y3s0QcBv8LuT1vq6GPM9i2bx6tcu4pTZNazd18KBlk4mlYeoLg12D2nu/LdEds9VRERkjCq03fkh979fTHvM0s8Nba21ayH//CVjzEzgUuBbwOf781oycOmVtIbWrnlj7dE4AZ9JnTLQ10qa5/gZ1UTjSV7ceohZE8uoKQ2ypb4t4xpV0kREZLwpdHXn3GIPpAf/C3wJqBym1xfSt+DIDmkJgn5farPbfHuk9eSUWRMAaGiNctLMGqpLgzS2Zx4XFXc3zs1xpruIiMiY1Gu70xjzpbT3r8567tt5PvdxY8yqHG9XFDIwY8w7gAOFzn0zxlxnjFlujFleX19fyKdIgdJXd6avruzwQprJvK6vZk0sZXKlsxhgYnmI6rIQLZF4RmsznlS7U0RExpd8c9Lem/b+V7Keu6S3T7TWXmitPT7H218KHNtZwOXGmG3A74HzjTG/7eX1brfWLrLWLpo8eXKBLyGFSN8nLX1OWlvWnDRfP0tpxhgWHeFU0ya6c9KshZbOGC9vP8QND65OVdLU7hQRkfEiX0gzPbyf6+NBZa39irV2prV2Dk5YfNJa+8/FfE3JzZs7mEyS0e7sjHmVtIG1OwHe6Ia0CeUhakqd/dOaOmL8Y80Bfv3cttTGuQppIiIyXuQLabaH93N9XDBjzDuNMbuAM4CHjTGPuo9PN8Y80t/7SnH0VElLLRzwDWzhAMCiOc5i4ckVYardkNbYHqM14qzy7Iw5W3+o3SkiIuNFvoUDJxljmnGqZqXu+7gfl/T3Ra21DwAP5Hh8D/D2HI8/BTzV39eTgfFCWNJath1sI+T3EU0kaY8mmFQeTq3uHMiB6CfNrObm95/C+cdOYc0e58esqSNGW8QJZ97pBqqkiYjIeNFrSLPW+nt7XsYHr0IWS1g27m9l4fQqVuxsdFd3mlQ4G0i70xjDO06cDkBNmVtJ64jR0umcatCRqqT1/zVERERGk0I3s5VxzCuQbT/YRiSe5PgZVYDTegz4ffjdn6KBtDvTVXlz0tqjtEWckOa1O7OPphIRERmrFNIkLy98rdvXAsDx06tTz6XPSfMPpJSWpjpt4UCbez6odxxVQu1OEREZJxTSJC8vpK3d68wVWzi9KvVcKNC1unOQCmmEA35Kg35n4UBWu1Nz0kREZLxQSJO8vALZun0tzJ5YRlVJMPXcYK3uzFZT5pzf2ap2p4iIjFMKaZKXt0ltfUuEeZPLCQa6fmwCg7RPWrbq0iCNHbFuc9K0BYeIiIwXCmmSV3qFrCIcIOjv+jjoN2khbXAraYfaorS5c9FSqzuV0UREZJxQSJO80itk4YCfkL/rxyaYtrpzEDMatRVhth9sT33sLRxQu1NERMYLhTTJK32T2nDQRzAtpAV8vlQ7dLBWdwJMrgxnnG7Q4W5mq9WdIiIyXiikSV6ZlTQfgSFod06uDGd83KnVnSIiMs4opEle6eErFPAR9KUvHDCDcixUttqKzJCmdqeIiIw3CmmSV3pICwf8+HyGgFteC/q72p2DubqzWyUtrmOhRERkfFFIk7zSCmeE3e03vHlpQb8vVUkb1HZnD5U0zUkTEZHxQiFN8sqspHnhzHnM2czWu27wXjO7khaJOyU0tTtFRGS8UEiTvDJCWtAPOHPTwNnM1puL5hvElDaxPJRzSw9V0kREZLxQSJO8sld3grP1BkDIX5xjoYJ+HxPLQt0e1+pOEREZLxTSJC+Tq90ZcNudGXPSBvd1s1d4gtqdIiIyfiikSV7ZJw5A18KBgM+k2pyDuQUHdJ+XBlrdKSIi44dCmuSVfpKAV0kLpa3u9J4e7Era7Ell1FVlBjW1O0VEZLxQSJO8cq/uTNuCowhz0gC+fMmx3PPR0zMeS6jdKSIi40RguAcgI1969goHs7bg8BuSdvDP7gSoLg1SEsz8PUKrO0VEZLxQJU3yyj5xAJwFA+CENX+R5qQBGUdQAViFNBERGScU0iSvXO3OkL9rK45ire4EZ++19Puq3SkiIuOFQprklXt1Z9fZnSa1cKAIKY2uqh1odaeIiIwfCmmSV8Y+acHshQPpm9kW5/WDaTfW6k4RERkvFNIkr1xbcATTjoUq5pw07zU8aneKiMh4oZAmeaVXyELZ+6T5TKrN6S9SSAumhzRV0kREZJxQSJO80itkobSTBsCpqPlSB6wX5/W9+W+g1Z0iIjJ+KKRJXumVtNTWG4GusOYVuorX7uy6r9qdIiIyXiikSV65Vm1mHgtVnBMHPOl7pWl1p4iIjBcKaZJXrvCVfuJAsVd3plfStLpTRETGC4U0ySvXXLP0szt9RTq70xPwaXWniIiMPwppkleu8JWam+YbgnanKmkiIjIOKaRJXrnnpKW1O4t4LBRkbsGhkCYiIuOFQprklSt8ecEp4DepdmjxjoXS6k4RERl/FNIkr1xba0woCxHy+ygLBboqaUXbJy2tkqbVnSIiMk4EhnsAMjpdccp0TpldQ0U4UPxjodJKeTpxQERExgtV0qRfwgE/R9VVAl3hrHhbcGh1p4iIjD8KaTJgXiWteGd36lgoEREZfxTSZMC8cFasdqcOWBcRkfFIIU0GrOirO3UslIiIjEMKaTJgvqLvk5a2ma3mpImIyDih1Z1SkIuPq+OihVNzPpc6u7NIKc3bJy3oN2p3iojIuDEslTRjzNXGmNXGmKQxZlEv19UYY/7PGLPOGLPWGHPGUI5Tutx2zSLe9caZOZ/zpeakFee1vXZnyO9TJU1ERMaN4Wp3rgKuAp7Oc92Pgb9ba48FTgLWFntg0nf+Ih+w7rU7w0G/joUSEZFxY1jandbatdD7akBjTDXwZuDD7udEgegQDE/6yFvdWawtOLx90kJ+n9qdIiIybozkhQNzgXrgV8aYV40xdxpjyod7UNKdl82KfcB6OOjTsVAiIjJuFC2kGWMeN8asyvF2RYG3CABvAH5urT0FaAP+vZfXu84Ys9wYs7y+vn4QvgIpVLGPhQq69w/5fTpxQERExo2itTuttRcO8Ba7gF3W2mXux/9HLyHNWns7cDvAokWL9C/5ECr2nLRUuzOgdqeIiIwfI7bdaa3dB+w0xhzjPnQBsGYYhyQ9GKp90sIBre4UEZHxY7i24HinMWYXcAbwsDHmUffx6caYR9Iu/RRwjzFmJXAy8O2hH63k44Wzou2T5rU7Az6t7hQRkXFjuFZ3PgA8kOPxPcDb0z5eAfS4j5qMDEPV7gwH/JqTJiIi48aIbXfK6OHzFbfdWRL0A1Aa9KOMJiIi44WOhZIBKwv6OX5GFUdPrSzK/S86ro4bkyeyYX+LKmkiIjJuqJImAxbw+/jrp87hLcdMKcr9q0qCvOdNs/D7dHaniIiMHwppMmr4fEarO0VEZNxQSJNRw2+cSto5Nz7JH17aMdzDERERKSqFNBk1fD6DtbDzUAdbG9qHezgiIiJFpZAmo0b6Ae6xhA7xFBGRsU0hTUYNf9pPazSukCYiImObQpqMGukHuCukiYjIWKeQJqOGP2233KjanSIiMsYppMmo4VclTURExhGFNBk1fKqkiYjIOKKQJqOGP+1sUFXSRERkrFNIk1EjY06aQpqIiIxxCmkyamSs7lS7U0RExjiFNBk10itp2sxWRETGOoU0GTW0ulNERMYThTQZNXyakyYiIuOIQpqMGunHQkUU0kREZIxTSJNRw6eFAyIiMo4opMmokR7StHBARETGOoU0GTW0T5qIiIwnCmkyavi0ulNERMYRhTQZNdIrafGkJZm0wzgaERGR4lJIk1HDn/XTqsUDIiIylimkyaiR3u4EhTQRERnbFNJk1OgW0jQvTURExjCFNBk10uekgUKaiIiMbQppMmqokiYiIuOJQpqMGtmVNG1oKyIiY5lCmowa3upOL6vp/E4RERnLFNJk1PDanZUlQUCrO0VEZGxTSJNRoyukBQDNSRMRkbFNIU1GDW9OmldJ05w0EREZyxTSZNRQJU1ERMYThTQZNbxKWpVCmoiIjAMKaTJqeKs7tXBARETGA4U0GTVMVruzpy04rLV85vev8sTa/UM2NhERkcEWGO4BiBTKnxXSelo4cLAtyl9W7CGRtFywoG7IxiciIjKYVEmTUSPg9+akue3OHippG/e3ArBmT3OP93p5+2E6oolBHmGX1kicmx5dR2eseK8hIiJjm0KajBozakr5xuXHcdlJ04GeQ9qmeiekbWloozUS7/b8fct38q6fP8c9y7YXbaxPb6jnlsWbeWr9gaK9hoiIjG0KaTJqGGP40JlzmFwZBnoJaftbUu+v3ZtZTTvQ0snX/7IagP3NnUUaKRxw7/3y9sNFew0RERnbFNJk1An4DMb0vLpz44FWpleXALB6d1PGcyt2NNLhtiAPtkWLNsb61ggAy4cxpCWTVhv+ioiMYsMS0owxVxtjVhtjksaYRb1c9zn3ulXGmN8ZY0qGcpwyMhljCPp9PYa0TQdaOXN+LbUVIVZlzUvb2tAGwPwpFUWupDkhbdXupmGbl3bDQ6u59CdLh+W1RURk4IarkrYKuAp4uqcLjDEzgE8Di6y1xwN+4L1DMzwZ6cJ+X6rduWRDPR/+1Yus29dMQ2uEAy0R5k+pYOH0ala7Ie1Xz27lg798kQ37W6mtCHHUlAr2u0GqGA60OPeOJSyr0qp5Ow+1Z3xcLNsPtnHPsh1s2N/K4SJWDAt19wvbWfQ/j3PpT5aSSNrhHo6IDJNVu5u46mfP8vzmgwVdv7uxA2vH798ZwxLSrLVrrbXrC7g0AJQaYwJAGbCnuCOT0cLvN/zq2W386z2v8J1H1vLU+nouv/lZvv3IWgDOO2Yyx0+vYuP+Fu5+fhvfeGgNT2+o5x9r9jG3tpy6qpKMSloyaQd1tWd9S4STZ9UAmS3P//7rGq77zfJBe52e/PypzakwtHZfz6tch0JLZ4zv/W0dAKv3NPPqjpExT2/Jhnqu/fVLOrlCZAj96ZVdvLKjkQ/c+QLr8vzdtKexg3NvXMyPn9g4RKMbeUbsnDRr7W7g+8AOYC/QZK19bHhHJSNFY3sMgIdf38u6fS188eJjqAwHuP+V3Zw2dyLHTq3iuOnVxJOWG/++nnmTywFo7owzt7acKVVhWjrjtEed1Z/ff2w9b75pccEt0O0H2/jaX1ZxoCX39QdaIiyYVsnc2vKMxQOrdjexp6mTQ0WsbkXjSR5+fS/nHTMZgLV7W/J8Rpd1+5q56EdL+Mea7hsB92V+267D7Ty9oZ7v/G0tn/vDa7RG4vzkvScT8BkeXzv8K16bO2N88Y+v8cS6A7ywpbDf6IfaWK84bj/YxjZ3+kFPmtpjvLazMeMXqFW7m7hl8SY2u6u4e1LfEiESH9tb4MQSyaJO2yiG5zcf5Oi6CpIWnt3U9f9ePJHs9jP/3OaDxJOWnz21mZ2H2od6qCNC0UKaMeZxdy5Z9tsVBX7+BOAKYC4wHSg3xvxzL9dfZ4xZboxZXl9fPzhfhIxY71k0k4+fO49zjqplenUJ1549l69dthCAj5w9F4DjplcB0BKJ80+LZjGv1glqc2srqKt0pjd6c8eeWHuA+pYIn/rdqxn/ILRF4t3+4li/r4VLf/IMv3l+O/e9tLPb2OKJJAfbIkyuCPOG2RN4ZfthrLUcbouyp8n5C3Xd3sGvbllr+a8/r+JbD6+hpTPONacfQW1FmLV7m7HWsn5fC8le/uE/0NzJB+5Yxob9rfzkiY0ZLYY/Lt/J0V/9G++59XkOuosiovEkT2+oZ0vWP5ab61s558bFfPCXL/LLZ7byxLr9nD2/ljPn13Lq3IkZJ0G0R+PE3fBnre11fAOx81A7l/30mdSf1/cfXU9Da4RwwMejq/fx0rZD/H3VPu5cuqXHf/yttTy8cm9qXmNvNu5vKajN/MflO7nhwdUs33Yo4/HOWIJzvvckH/7Vi6nvd/ZYAJ5Yu59zb1rMd91KZX/FE0k++/tX+cFj3Rsca/Y005ZjK5t0WxvauO43y3lmYwMAB1sjPLp6H//z1zW87/YXeHl719dnreUL973GuTc9xfk/eIpbFm/qdr8n1u7npkfX8ZYfPMUVtzzLhT9cQlskzuu7mrj85me46dH13ca6cX8L3/v7OpraY+xt6uAt33+K7/2tkIZNl5e2HeKGB1d3+4Vk2ZaDPPL63kELfS9vP8yX/28lT64b2Kkov3hmK+d8bzGbDnT9IratoS31y2df9bWtuH5fS59+4TzcFmXdvhYuP2k6U6tKWLmrEXA6Gdf84kWuvOXZjL9/l205SGU4gM/A7U9v6dPYxoqinThgrb1wgLe4ENhqra0HMMbcD5wJ/LaH17sduB1g0aJFY/tXUOHGd58EOP+4tEUTlAT9XHHyDE6ZNYHZk8oAmD2xjMpwgJZInLcurGNLfRtbGtqYW1tGRdjZEHdfcyeVJQHW729h0RETeGnbIa6+7Tne9YaZXHHyDC7536e56Lg6/ufKE/jVs1t5dlMDAZ8Pg7P4YPH6eq4//yhufnIjD762h2nVpXz6gvlYC5OrSphaXcqfXtnFB3/5IpMrwqnxr9nbzJnza7t9XbsOt7P7cAfGGKy1nDZvEk0dMSKxBFOqnGDZ0BohaS1TKrvW0XTGEqze08zdLzh7v5WH/Jw1v5YF0ypZt6+ZO5Zu4duPrONT58/nCxcdk/N7etvTW2jsiPHhM+fw6+e2sWJnI0dMKqcjluC3y3ZQV1nCS9sPceczW+mIJnjg1d00dcSYM6mMJ75wHn6fs9nwE2v3Yy3c8cFFnD5vIkkL4YDz++BbF9bxjYfW8NrORqpKg1x+8zPEEkn+7aJj2FzfxrItB7n3Y6cztbqESDzBnUu3cvb8Wk5yW8e9aemMcbgtxrSaEoL+rt8/D7VFec9tz7O3qZNfP7eNo+oquPuF7XzojDkcaOnkdy/u4J5lO1LXL9lQz93Xntbt/mv3tvCv975CadDPd646gStPmQE4/7B5R5YB/PjxjfzvExuYUVPKO0+ZQcDn45ozjnD/XJ1fHlbvaeahlXu4bckWgn7D71/awZIvvoU698/4ibUH2NPUyb7mTr724Gpuef8bAOcfszuWbuHWJZt536mzuXXJZgAeem0PX77kmIxxpNvW0EZFSYDatJ/BdN9/bAN/XrEHY+DCBXVsOtDKY2v28Z5Fs7j2ruWUBH18+ZJjKQ36iSaSXHnKDLY1tPHzpzZzyfFT+c8HVtEaibNqdxOPf+FcPvuHFSzd2IDfZ5hQFuI9t73A3deeyqIjJvLbF7bzp1d28eEz59DQGuGmR9dzyqwaTp83ib+8tpuTZ03gk/e8QiSe5KRZNXzi3Hl8+5F1PLZmH39ZsYeq0iDnHT2Zv6/eR2skTkU4wPObD3LtXS/RHk2weN0BJlWEaI3E+dMru/jSJcdQEvSnvn9bGlqZP6Wy2/dg2ZaDXHvXclojcU6fN5FLjp+GtZZIPMnHf/syje0xjp9Rxf3/chahQGH1jVgiyU+f3MRLWw9x8uwaPnLWXCaVh/jqn1exdm8zf1i+k89ccBTPbGpgwbRKPnnefKbXlHa7zxNr93PH0i184aJjeNOcianH/7FmP9FEkuvvfZXJlWEisSQvbjvEnEll/OwDb2Sh+4tq+s/ObU9v4VtXHs/bTpiWeu7xNfv53H0r8BnDTe8+kec2H+Tdb5zJ8TOqe/zattS3ctnNzzB3Ujl/uf4swgEf1oLPl/kzeP8ru7hv+U6mV5eyYJozntPnTeL13U28uqOR//rzKpo6YjzvVrTfd8cLvHVhHZ8490iWbT3E6UdOwm8Mj67exzcuP441e5tZs6eZE2dV87W/rCaRtPzoPSen/t4/0NJJNJ5k5oSybmNevP4A8ydXMGti13PWWloi8dQm6SPNSD4WagdwujGmDOgALgCKP5lHRpWA30d1addfmN7/qOD8ZXHSrBoOtHQyb3IFZx9Vy30v7+SYqVUkks5vyq/vamLTAady8h+XLuBAc4T//usavvHQGu5cupUDLRHuXbaDmtIQN6f9xv+xc+ZSGvRz8+JN/OKZrXz/sQ28aY5TNXvXz58HYEplmLlu9W6pW2EAqAgHWLev6zffDftbWLmriVNm13DNnctS1TaAa04/gntf3EEiafn3tx3L1W+cyaU/WUo44Ofxz59LKODjryv38Pk/vMYRk8oIB3xMKg9xxpG1lAT9LJhWxR1Lt7BmTzM1ZUFuXryJ46ZXc/ZRtQR8hpKgn3X7mvn7qn3cu2wHl580nX+7+Bj+9MouvnL/6zR1xGjuiNEWTfDVSxewbOshfv6UEwwuO2k6cyaV8dMnN/HXlXu44mQntDy1vp5j6ip568LuR3JdvWgWtyzexNceXE0kliDgM5wwYwLf+du6VMXyg79cxl/+9Wz+++E13LtsBzc9up7vXHUC7zt1ds6fgUTS8vDre/nqA6+n2tm/+vCbmON+7294cDUNrRGuesMM7n9lN5/5/QrqKkv4t4uP4ZmN9Tzy+j4+fu483nHCdP62ai8/e2ozW+pbmTe5IvUzcsNDq5lUHiLgMxw/o4rP/mEFr+1qZEt9G5sOtPKba0/lyMkVPL/5ID96fAMXLazjpW2HuHnxJqyFHz2+ITXes+ZP4rnNB7EWLj1xGl9469Fc9KOn+fETG/nchUfzq2e3snz7YaZUhrlwYR0PvLKb5s4YNgm3Pr2Znz+1mQllQX721GZmTyzjvafO4sa/r2fHoXaOmOS19WP8fdU+Zk8sY2tDG19/cDUza0p5+NPnUBryZ3z/1u9r4banN/POU2aw4c5Q8QAAFPlJREFUdGMDV9/2fGqe3tKNDdRWhDhhRjXfeGhN6nN+/ew2JpaHWL79MH9btY95k8v51gXH85nfr+CLf1zJ0o0NXPfmeVz35nkE/T7eecuzfPp3K4gnkzS2xzhr/iS+9o6FRBNJVuxs5Jt/XcNHz5nHv/3xNapKAkTiSf72mXNYMK2KZNJy13Pb+f6jG9jd2MGXLzmWNx4xgT+v2MMTa/dz8XFT+fKfVlJXVcKnL5jPNx9aw7p9LZxzVC1LNzZwy+JN1FaEqa0Is2F/Cz9+YiN3feRUzj16Mp2xBB+/+2Ve3n6Y1kicGTWllIf93PviToJ+H1/6v5UcOaWCxvauX15++I8NfPbCo1LBb+WuRr7+4GrOnl/Lu984k9ZInOOmO+Hm9y/u4CdPbOTYqZXctmQzdz+/nTcfXcvavc3c+K4Tuf/VXfz4iY1UlwZZtbuJpzc0cPe1pzK9ppSHXtvD+v0tnDSzhk/e8woA//XnVTzy6XPw+QxN7TFe3XGYo6ZUsG5fC9ZCMGD46Nlz+evKvbz/zhe47+NncHSdE0hvXryJH/5jAwGf4Y8v70qFtETS8q1H1lJbEcYYuO7ulwF4YctB7vigs/lCeqiJJ5Lct3wX9yzbTsBnWL+/hQt+sISOWIJwwMeD15+d2svyhS0H+fx9r3Hk5HJe2d7I/a/uZuaEUk6cWcOJM2t4dPX+1C+WJ8+q4ao3zOCXz2zlpkfX8+qOw+w41M4Hz3A6An9fvY9XdzamAu706hLaYwlaOuPc8+J2vvK2BWxtaOPqW58nGk9w/yfPzAjjr+9q4v/96iWOqavkr58+m0TScs+yHfzm+W3sOtzBPR89jdPmTsQYQzJpsZD6xXM4meFYNWGMeSfwU2Ay0AissNZebIyZDtxprX27e903gH8C4sCrwEettXmX5C1atMguX648J7CvqZOEtcyoKcVay+b6NuZPqaC5M8aJN3RNcSwP+Vnx9YtSFZhvP7KW25/ewimza9i0v5WWSJwz5k1iwbQq/vDSDh793JvZ3xzhXT9/DoBT507k3o+expq9zVx+87MA3P/JMzl5Zg3//fAaqkqC/PiJjUyrLuGoukpW7mpkRk0pZSE/L21z5qz5DCQtfOeqE5hSGebGv69n/f4WTp83kbJQgCUb6plbW86W+laSFj585hxKQ35+89w2OuPOfI7LT5rOd991AgGfj1DAx9aGNu5YuoWKcIDr3jyPa+9azurdTfiMIeg3zK+rZNXuJhJJS1nIz4PXn8X8KZUs2VDPx+5aTmnITyjg42BrhBf+4wJ2Hmrn6luf5+PnHsmXLzmWZNLyth8vpaE1wrveOJPDbVH+vGI3HzlrLl95+4Kcfya/fWE7X/3zKirCAW75wBs4YUY1F/zgKQJ+H/99xXF84revcExdJev3t/CRs+ayek8TG/a38NFz5qW2NNnS0EZbJMHVi2by0Gt72HW4gxNnVnP1oln88LH1dMaSvOPEaXTGkzz02h4+/9ajufLkGbz5psUA3HbNG7n4uKmpn5Gp7r56B1o6Oeu7T3L2/FouOm4qnTGnYrhyl7Mi98IFU7j5/W/gK/e/zkOv7aEk6Hx/ktby1gV1rN7TzMG2CEu++BaaO2Mkk7DjUDvPbmrg2KmVvLTtML98divnHzuFb7/zBOqqwhhjuOHB1fz6uW3UVoRoaHXaRx85ay4XLJjCB+5cxrRqZ6FL0sL7T5vNly4+hu88so5rzjiCkqCfC3+4hPOOmUx5KMCc2jL+8NIuGtLapEfXVbBhfyuXnTSdj50zl+c3H+TVHU6raU9TB1vr23j6S29hd2MHf351N5Mqwry8/TCPr93P5y48muvPn8+tSzYzZ1I5JUEf197l/P16zelHkLCWT59/FFOrS/jaX1bxm+edf7yf+8r5qWrvip2NvPf25zlt7iTed+oszjtmSirkPLp6Hx+/+2V8BsIBPx2xBKfOmch9nzgjNf7vP7qemxdv4uz5tdz5oUWE/D7OuXEx7dE4/7+9O4+vqrzzOP75ZSMLWYCwQyBAkEVEtgAqRdDpoIKidSmilanL6KhDW3V0xNZaB0fHjrZ1QR2LtVWLGyrDKBVZREXZt7BKAgIJkEDIQpab5OaZP+5JSEKAQCG5we/79bqv3POcc0+ee353+d3nec55UtrFsnxnLm/cOpyLUhKprHRk5pXQMT6S0U8vJjOv5KjX4IgerXnrthHc/946Zq/O5PqhXejdPpYbhyfx8ucZ1QPVYyPDKCytoEdiDAvuG83UWWuZsy6LhOhwnrl+IC8uSmfld4doERaCr6KSiNAQQkOMX03ox4LN+1m9K49e7Vry9h0jyDhQxH/N28Ln23JIah3Nx/86itziMl5clM5PL0zmYJGPya8uo7jOCUyhIUavti254wc9uO/ddTxyRV9uHtmN91dl8vAHG3jvzpEkJ8bQpkYr6a6DxfzopaWUlPmZOKgT+/JLWbAlm6sGdiIhOoK3lu/il+P7sWRbDgcO+1izK4+XbhpM345xPPJhGn06xPI/X+yo3l/b2BakJrfmyWsG8OD76/l4wz5iI8N44uoBFJSW8/nWHBKiw/lobRapya15YfJgcgp93DJzOSFmzPvZKLLySknPOcwFPdsQGxnOV9sPMPnVZQzoHM9vrxtI29gWtI6JAOA/P97My0sySE6M4a3bhxMdEcbQ/5hP/07xrN2dR2R4CKXllTxz/UDmrt/LpqwCJqUm8eoXGYSHhRBigUs1TTivE2lZ+UwY2In3Vu5m6/5CSssrGZSUQHr2YQpKKxie3JpducWEhhjFZX5Geq18nRIief2nqbQIq/2j5kwxs1XOuaMuSdYkSdqZpiRNTsQ5x7Dpn1Hud/xj//Z0bRXNvZekVK/3Vfh5Zv42rhvSlSJfBcVlfkb0CPzKqupi8Vc6fj1nI11aRTF5RDdatgg0TA/6zaccKi7nywfH1Gpyn/LactrHRtIhPpLfL/iWwUkJ5JeUc83gLgzp1opHPkxjdO+2/HJ8YGzdzgNFzFqxm3vG9sJf6bj5j8sw4M7RPZn51Q5W7DxEeKgxoHM8j088l+cXbufesSm1ujjqKiwtZ9oHabSOiaC4rIKsvFL6dYrjrtE9iY8Kr9VVsTErn5iIMELM+C63iFEpgRMRcgp9JLaMqO5a27a/kPvfXceGzMD2h30VvHvnyFrdMnWP/Yqdhzi3cxzREYFjlpFzGDMjOTGmOmH5ychuPDqhP5v3FjD+uS8B6NYmmuiIMHq0jSGnwMfynbn06xjH3WN68cP+7QkPDWF3bjHPLfyWT9L2UVLm566LezL1khTCQkO44eWvaRvbgue97sP6vLh4O8/O30a5/8hn47VDujB79R5m3HQkuTvsq8CAvfkl/Pen21iafpD8knKmX30uk4d3O+b+N2UVkNK+Za0uWV+Fnyc/2cJHa7N47Mr+pGXmM+XC7rRt2YJh0z/jUHE51wzuTFR4KL+a0K/WF4dzjuFPLCC70Ff95XVRr0TuGduLIl8F8VHhDOyawIuL0vndgm1UfeQnJ8bgnGPnwWIeHNeHuy7uWaueuw4W8/SnW3n8qv4kREfUWvfrORv528Z9zP/F6OrXPQS69+5/dx0d4iKPStLLKiqP2U04Y3E6T83bwis3DyEtM58f9u9Qq6utoLScj9Zkct3QrtXJ3fbsQqZ9kMa+glJuTE3in0f3PGq/WXklZBf66JQQyV++/o75m/Yzpk87ZixOr+52/vmlvZl66ZH3fl5xGTMWp9OtTQzjB3bk2fnb+EHvtow5px3l/koWb81h2gcbyC70ERcZxt1jenHd0K788qM0Ssv8LN+RS6GvgoTocErL/cy6Y2T1md5w5ISQ+lppdh4oYuGWbApKy+nbMY5lGbm8tnQHb942nBHJbbjtzytZuCWbiLDAJYgSosNZOe1SwkKPPq67Dhbz6Jw0vsnIpWNCJKN6JfLgZX1YsyuPya8uAwLvp4SocDolRPHCjYOr3//OOX776VbCQkJIiA5n/Z58PliTSVxkGAWlFUy7vC+3jUo+qnt91vJdPDR7Ay3CQqiodMRHhTNzyrBaz79KSZmfu95cxdRLUhiU1KrWOucc+wt81T9iAJ78ZAsvfZ5OXGQYs+4YybIdB5lyQXfmrt/LvX9dA8C4/h14YNw5lFVU8ot31rF5bwGdE6LIzCshxODZG87nm4xc1u7Oo0+HWCalJpGa3Jp5afu4841VpLRryY4DRbSKiSCn0Me5nePo2iqaGTcNOar+p5uSNJE6CkrLiQ4PrfcD7u+RXVDK3zbu46YR3Wp9iFWNXSot95NT6KvVhVBzfUPkFPrYnn2YQUkJ1V9aTalq7I4ZZOQUVY89ORUV/kq27i+s7jICeH3pTmIjw7h6UOfqY1RZ6diYVUC/TnH1fuFV1anm8an6vDvRcc4uKCWvpJzC0grW78njlpHdySspr/6lXx9/pWN3bjHd2kQ3OI711bnuY99ftYfDvgpuuaD7MR8388sd7Mot5qHL+lBa7j8qqaqyMSufLXsLuaBXGzrGB8Y+5RaV0So6/KTr7Kvwn9ZWhkNFZbQ6zvE9XUrK/Dz+f5tYlnGQSalJ3Daqx0nvY93uPJ78ZAsPX96XAV1qj9ual7aP/12fxRMTBxAXFXbKrwU4kqxUtfRW+Ct5eUkGuUVl9OsYx3ld4klpf/T4uuMp91dy0VML6dWuJa9NSW3w+LpXlqTz/MLtTL96QPX8yfXZsCef91fvITYyjEmpSfWOsTtVCzbvJyIspPoHIwTG4/787bVc2rc9PxrSpbrcX+koKCknPiqcDZn5dG8TQ3z0scedbdtfSI/EGPYX+mgVHc47K3Yze00mYSHG7H+58LQ9h2NRkiYiIiLkl5QHzpo8yTFXlZXupB8jDXOsJC2YTxwQERGR0yw+6tTOZFSC1viC9mK2IiIiIt9nStJEREREgpCSNBEREZEgpCRNREREJAgpSRMREREJQkrSRERERILQWXkJjoycIm54+eumroaIiIjIKVNLmoiIiEgQ0owDIiIiIk3oWDMOqCVNREREJAidlS1pZpYDfHcG/0UicOAM7l/OPMWweVP8mj/FsPlTDE+fbs65tnULz8ok7Uwzs5X1NUtK86EYNm+KX/OnGDZ/iuGZp+5OERERkSCkJE1EREQkCClJOzWvNHUF5O+mGDZvil/zpxg2f4rhGaYxaSIiIiJBSC1pIiIiIkFISdpJMLNxZrbVzLab2UNNXR+pn5nNNLNsM0urUdbazOab2bfe31ZeuZnZH7yYrjezwU1Xc6liZl3NbJGZbTKzjWY21StXHJsBM4s0s+Vmts6L32NeebKZLfPi9LaZRXjlLbzl7d767k1ZfznCzELNbI2ZzfWWFcNGpCStgcwsFHgBuAzoB0wys35NWys5hj8B4+qUPQQscM6lAAu8ZQjEM8W73QHMaKQ6yvFVAPc55/oBI4C7vfeb4tg8+ICxzrmBwPnAODMbATwFPOuc6wUcAm71tr8VOOSVP+ttJ8FhKrC5xrJi2IiUpDVcKrDdOZfhnCsDZgFXNXGdpB7OuSVAbp3iq4DXvfuvAxNrlP/ZBXwDJJhZx8apqRyLc26vc261d7+QwJdEZxTHZsGLw2FvMdy7OWAs8J5XXjd+VXF9D7jEzKyRqivHYGZdgCuAV71lQzFsVErSGq4zsLvG8h6vTJqH9s65vd79fUB7777iGuS8bpNBwDIUx2bD6yZbC2QD84F0IM85V+FtUjNG1fHz1ucDbRq3xlKP3wH/BlR6y21QDBuVkjT53nGBU5p1WnMzYGYtgfeBnznnCmquUxyDm3PO75w7H+hCoCeiTxNXSU6CmY0Hsp1zq5q6Lt9nStIaLhPoWmO5i1cmzcP+qu4v72+2V664BikzCyeQoL3pnJvtFSuOzYxzLg9YBIwk0A0d5q2qGaPq+Hnr44GDjVxVqe1C4Eoz20lgeM9Y4Pcoho1KSVrDrQBSvDNbIoAfA3OauE7ScHOAW7z7twAf1Sj/iXd24Aggv0Z3mjQRbyzLH4HNzrlnaqxSHJsBM2trZgne/SjgHwiMK1wEXOttVjd+VXG9FljodBHPJuWc+3fnXBfnXHcC33cLnXOTUQwblS5mexLM7HICffShwEzn3PQmrpLUw8z+ClwMJAL7gUeBD4F3gCTgO+B651yulww8T+Bs0GLgn5xzK5ui3nKEmV0EfAFs4Mh4mIcJjEtTHIOcmZ1HYBB5KIHGgHecc78xsx4EWmVaA2uAm5xzPjOLBP5CYOxhLvBj51xG09Re6jKzi4H7nXPjFcPGpSRNREREJAipu1NEREQkCClJExEREQlCStJEREREgpCSNBEREZEgpCRNREREJAgpSRORZs/M/Ga2tsbtoRM/qt79LDazoae7fg34vxO9CeRFRKqFnXgTEZGgV+JNQdRcTQTmApuauiIiEjzUkiYiZyUzG2dm79ZYvtjM5nr3Z5jZSjPbaGaPNWBfw8xsqZmtM7PlZhZrZpFm9pqZbTCzNWY2xtt2ipk9X+Oxc72LgWJmh81surefb8ysvZldAFwJPO21AvY8zYdCRJopJWkicjaIqtPdeQPwGTDczGK8bW4gcKV0gGnOuaHAecBo7wr59fKmgXsbmOqcGwhcCpQAdxOY530AMAl43bvq+vHEAN94+1kC3O6cW0pgSp0HnHPnO+fST+H5i8hZSEmaiJwNSrwEp+r2tnOuApgHTPAmfL6CI/MMXm9mqwlMa9MfON54sHOAvc65FQDOuQJv3xcBb3hlWwhMU9X7BPUsI9CtCbAK6H6Sz1NEvkc0Jk1EzmazgHsIzCW40jlXaGbJwP3AMOfcITP7E3CiFrCTUUHtH8A1911eY9JpP/oMFpHjUEuaiJzNPgcGA7dzpKszDigC8s2sPXDZCfaxFehoZsMAvPFoYQQmgJ/slfUmMOn7VmAncL6ZhZhZVyC1AfUsBGJP4nmJyPeAkjQRORvUHZP2JIBzzk+ge/Ey7y/OuXUEujm3AG8BXx1vx865MgLj2Z4zs3XAfAKtYy8CIWa2gcCYtSnOOZ+3vx0EztT8A7C6AfWfBTzgnYCgEwdEBAA70vIuIiIiIsFCLWkiIiIiQUhJmoiIiEgQUpImIiIiEoSUpImIiIgEISVpIiIiIkFISZqIiIhIEFKSJiIiIhKElKSJiIiIBKH/B0DRgKg9q/hDAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=[10,5])\n",
    "plt.plot(counts1, values1)\n",
    "plt.xlabel('Eval count')\n",
    "plt.ylabel('Energy')\n",
    "plt.title('Convergence with noise');\n",
    "plt.axhline(ref)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Lower T1 values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "NoiseModel:\n",
      "  Basis gates: ['cx', 'id', 'u2', 'u3']\n",
      "  Instructions with noise: ['measure', 'u3', 'cx', 'reset', 'u2']\n",
      "  Qubits with noise: [0, 1, 2, 3]\n",
      "  Specific qubit errors: [('reset', [0]), ('reset', [1]), ('reset', [2]), ('reset', [3]), ('measure', [0]), ('measure', [1]), ('measure', [2]), ('measure', [3]), ('u2', [0]), ('u2', [1]), ('u2', [2]), ('u2', [3]), ('u3', [0]), ('u3', [1]), ('u3', [2]), ('u3', [3]), ('cx', [0, 0]), ('cx', [0, 1]), ('cx', [0, 2]), ('cx', [0, 3]), ('cx', [1, 0]), ('cx', [1, 1]), ('cx', [1, 2]), ('cx', [1, 3]), ('cx', [2, 0]), ('cx', [2, 1]), ('cx', [2, 2]), ('cx', [2, 3]), ('cx', [3, 0]), ('cx', [3, 1]), ('cx', [3, 2]), ('cx', [3, 3])]\n"
     ]
    }
   ],
   "source": [
    "# T1 and T2 values for qubits 0-3\n",
    "T1s = np.random.normal(5e3, 1e3, 4) # Sampled from normal distribution mean 50 microsec\n",
    "T2s = np.random.normal(7e3, 1e3, 4)  # Sampled from normal distribution mean 50 microsec\n",
    "\n",
    "# Truncate random T2s <= T1s\n",
    "T2s = np.array([min(T2s[j], 2 * T1s[j]) for j in range(4)])\n",
    "\n",
    "# Instruction times (in nanoseconds)\n",
    "time_u1 = 0   # virtual gate\n",
    "time_u2 = 50  # (single X90 pulse)\n",
    "time_u3 = 100 # (two X90 pulses)\n",
    "time_cx = 300\n",
    "time_reset = 1000  # 1 microsecond\n",
    "time_measure = 1000 # 1 microsecond\n",
    "\n",
    "# QuantumError objects\n",
    "errors_reset = [thermal_relaxation_error(t1, t2, time_reset)\n",
    "                for t1, t2 in zip(T1s, T2s)]\n",
    "errors_measure = [thermal_relaxation_error(t1, t2, time_measure)\n",
    "                  for t1, t2 in zip(T1s, T2s)]\n",
    "errors_u1  = [thermal_relaxation_error(t1, t2, time_u1)\n",
    "              for t1, t2 in zip(T1s, T2s)]\n",
    "errors_u2  = [thermal_relaxation_error(t1, t2, time_u2)\n",
    "              for t1, t2 in zip(T1s, T2s)]\n",
    "errors_u3  = [thermal_relaxation_error(t1, t2, time_u3)\n",
    "              for t1, t2 in zip(T1s, T2s)]\n",
    "errors_cx = [[thermal_relaxation_error(t1a, t2a, time_cx).expand(\n",
    "             thermal_relaxation_error(t1b, t2b, time_cx))\n",
    "              for t1a, t2a in zip(T1s, T2s)]\n",
    "               for t1b, t2b in zip(T1s, T2s)]\n",
    "\n",
    "# Add errors to noise model\n",
    "noise_thermal = NoiseModel()\n",
    "for j in range(4):\n",
    "    noise_thermal.add_quantum_error(errors_reset[j], \"reset\", [j])\n",
    "    noise_thermal.add_quantum_error(errors_measure[j], \"measure\", [j])\n",
    "    noise_thermal.add_quantum_error(errors_u1[j], \"u1\", [j])\n",
    "    noise_thermal.add_quantum_error(errors_u2[j], \"u2\", [j])\n",
    "    noise_thermal.add_quantum_error(errors_u3[j], \"u3\", [j])\n",
    "    for k in range(4):\n",
    "        noise_thermal.add_quantum_error(errors_cx[j][k], \"cx\", [j, k])\n",
    "\n",
    "print(noise_thermal)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "NoiseModel:\n",
      "  Basis gates: ['cx', 'id', 'u2', 'u3']\n",
      "  Instructions with noise: ['measure', 'u3', 'cx', 'u2', 'id']\n",
      "  Qubits with noise: [0, 1, 2, 3, 4]\n",
      "  Specific qubit errors: [('id', [0]), ('id', [1]), ('id', [2]), ('id', [3]), ('id', [4]), ('u2', [0]), ('u2', [1]), ('u2', [2]), ('u2', [3]), ('u2', [4]), ('u3', [0]), ('u3', [1]), ('u3', [2]), ('u3', [3]), ('u3', [4]), ('cx', [0, 1]), ('cx', [1, 0]), ('cx', [1, 2]), ('cx', [1, 3]), ('cx', [2, 1]), ('cx', [3, 1]), ('cx', [3, 4]), ('cx', [4, 3]), ('measure', [0]), ('measure', [1]), ('measure', [2]), ('measure', [3]), ('measure', [4])]\n",
      "VQE on Aer qasm simulator (with noise): (-1.5882408399504115+0j)\n",
      "Delta from reference: [0.26903419+3.77606111e-17j]\n"
     ]
    }
   ],
   "source": [
    "device = provider.get_backend('ibmq_valencia')\n",
    "coupling_map = device.configuration().coupling_map\n",
    "#noise_model = noise.device.basic_device_noise_model(device.properties())\n",
    "noise_model = NoiseModel.from_backend(device.properties()) # method above is deprecated\n",
    "basis_gates = noise_model.basis_gates\n",
    "\n",
    "print(noise_model)\n",
    "\n",
    "backend = Aer.get_backend('qasm_simulator')\n",
    "quantum_instance = QuantumInstance(backend=backend, seed_simulator=167, seed_transpiler=167,\n",
    "                                   noise_model=noise_thermal)\n",
    "\n",
    "counts = []\n",
    "values = []\n",
    "def store_intermediate_result1(eval_count, parameters, mean, std):\n",
    "    counts.append(eval_count)\n",
    "    values.append(mean)\n",
    "    \n",
    "\n",
    "aqua_globals.random_seed = 167\n",
    "#optimizer = SPSA(max_trials=200)\n",
    "optimizer = SPSA(maxiter=200) #max_trials property is deprecated\n",
    "#var_form = RY(num_qubits)\n",
    "var_form = RealAmplitudes(num_qubits) #RY is deprecated\n",
    "vqe = VQE(qubit_op, var_form, optimizer, callback=store_intermediate_result1)\n",
    "vqe_result1 = vqe.run(quantum_instance)\n",
    "print('VQE on Aer qasm simulator (with noise): {}'.format(vqe_result1['eigenvalue'])) #'energy' deprecated\n",
    "print('Delta from reference: {}'.format(vqe_result1['eigenvalue']-ref)) #'energy' deprecated"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.lines.Line2D at 0x7ffa609fa580>"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmkAAAFNCAYAAABbpPhvAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nOzdd3zdZd3/8dd1VvZq0jbde9FCKRRo2ZQhIOgtiiIqKgiIMhScN3iL3urtT0VARaWg7I3svQqU7kH3HmmTpm2SZo+zr98f33NOT5qkSdMmadr38/HIg+R7vud7rpym9J3P57qur7HWIiIiIiKHF1dPD0BEREREWlJIExERETkMKaSJiIiIHIYU0kREREQOQwppIiIiIochhTQRERGRw5BCmohIDzHGDDXG1Btj3Ps5xxpjRnfDWP5pjPlFV7+OiHScQpqINGOMudIYszgWHnYaY940xpze0+M6Ellrt1trM621EQBjzIfGmO/00Fi+a6393554bRFpnUKaiCQYY24F7gF+B/QHhgJ/Bz7fk+NKZozx9PQYRES6g0KaiABgjMkBfg1831r7grW2wVobsta+aq39ceycFGPMPcaY0tjHPcaYlNhjZxtjSowxtxljymJVuG/HHjvFGLMrua1njPmCMWZF7HOXMeZnxpjNxpg9xphnjTF9Yo8Nj7X8rjHGbAc+MMa4jTF3GWMqjDFbjTE3xs7xxL8XY8y/YmPYYYz5Tfy1jTHfMsZ8Yoz5kzGmKvb8i5LG1ccY81Ds+6syxryU9NglxphlxphqY8xcY8xxbbyXvzLG/DX2udcY02CM+WPs6zRjjD/2OvHvzWOM+S1wBvC3WBXzb0mXPM8YszH2uvcZY0wbr3tn7L171BhTZ4xZbYyZmvT4hFi1rjr22OeSHnvYGPOb2OcFxpjXYudVGmNmG2NcsccGGmP+Y4wpj713N+/v50pEOk8hTUTipgOpwIv7Oed2YBpwPDAZOBm4I+nxQiAHGARcA9xnjMmz1i4AGoAZSedeCTwZ+/wm4L+As4CBQBVw3z6vfRYwAfgMcC1wUWwcJ8Sem+xhIAyMBqYAFwDJbcRTgPVAAfAH4F9JwecxIB2YCPQD7gYwxkwB/g1cD+QD9wOvxEPqPj4Czo59fhKwCzgz9vV0YL21tjL5Cdba24HZwI2xFuiNSQ9fErvOccCXY+9BWz4HPA3kAq8Af4uN3wu8CrwT+75uAp4wxoxr5Rq3ASVAX5yK6n8DNhbUXgWW4/wZnwv8wBizv/GISCcppIlIXD5QYa0N7+ecrwG/ttaWWWvLgV8B30h6PBR7PGStfQOoB+Ih4CngqwDGmCzg4tgxgO8Ct1trS6y1AeBO4Ev7tDbvjFX3mnCCyr2x86uA38dPMsb0j137B7Hzy3CC1hVJ19pmrX0gNhfsEWAA0N8YMwAn/H3XWlsV+z4+ij3nOuB+a+0Ca23EWvsIEMAJrfuaB4wxxuTjhLN/AYOMMZk4YfOjVp6zP7+31lZba7cDs3DCaVs+sda+EfveHsMJ08TGmRm7VtBa+wHwGrE/k32EcN6TYbH3YLZ1bvR8EtDXWvvr2DW2AA/Q/L0VkUNEIU1E4vYABe3M+RoIbEv6elvsWOIa+4S8RpxgAE7V7LJY5ekyYKm1Nn6tYcCLsfZaNbAWiOBUceKK9xlHcRuPDQO8wM6k692PUz2K2xX/xFrbGPs0ExgCVMaC376GAbfFrxm77pB9vv/4NZuAxTiB7EycUDYXOI3OhbRdSZ8nv6cdOTc19mc6ECi21kaTHt+GUxHb1x+BTcA7xpgtxpifxY4PAwbu8x78N83/nETkENEEXBGJm4dTGfov4Pk2zinF+Yd6dezrobFj7bLWrjHGbMOpVCW3OsEJWVdba+fs+zxjzPD4JZIO7wQGJ309ZJ9rBYCCdqqCrSkG+hhjcq211a089ltr7W87eK2PcNq7U4BFsa8/g9Mi/riN59g2jh8KpcAQY4wrKagNBTa0GIS1dTgtz9uMMZNw5gEuwnkPtlprx3ThOEUkRpU0EQHAWlsD/A/OPLL/Msakxya9X2SM+UPstKeAO4wxfY0xBbHzHz+Al3kSuAWnuvRc0vF/Ar81xgwDiF1/fytKnwVuMcYMMsbkAj9N+j524sy7ussYk22cRQmjjDFntTe42HPfBP5ujMmLff/xuWQPAN81ziIIY4zJMMZ8Nta6bc1HwFXAGmttEPgQZ17c1liruDW7gZHtjbOTFuBU1n4S+77OBi7Fmb/WTGyBxOjYPL0anKpmFFgI1BljfhpbAOE2xkwyxpzURWMWOaoppIlIgrX2LuBWnMUA5TiVkxuB+ArH3+C08VYAK4GlsWMd9RROu+8Da21F0vF7cSa5v2OMqQPm40zub8sDOEFsBfAp8AbOQoFI7PGrAB+wBmcRwvM4c6w64hs4c7LWAWXADwCstYtxFiz8LXbNTcC39nOduUAae6tmawA/bVfRwHkfvhRbVfqXDo63Q2JB8VKcSmYFztYqV1lr17Vy+hjgPZw5hfOAv1trZ8XmuV2CMydua+w6D+IsFhGRQ8w4c0FFRHov42yh8U9r7bCeHouIyKGiSpqI9DqxVtvFsf3FBgG/ZP9bh4iI9DqqpIlIr2OMSceZ8zUeaAJeB26x1tb26MBERA4hhTQRERGRw5DanSIiIiKHIYU0ERERkcPQEbmZbUFBgR0+fHhPD0NERESkXUuWLKmw1vbd9/gRGdKGDx/O4sWLe3oYIiIiIu2K3Y2lBbU7RURERA5DCmkiIiIihyGFNBEREZHDkEKaiIiIyGFIIU1ERETkMKSQJiIiInIYUkgTEREROQwppImIiIgchhTSRERERA5DCmkHYU99gJUlNT09DBERETkCKaQdhAdmb+VbDy3sttdbU1rLl++fR1Mw0m2vKSIiIj1DIe0g1DSFqA+EO3TujuomLvv7HCobgp1+veUl1SzcWsmuWn+nryEiIiK9g0LaQQiEIoQi0WbH3l2zu9UgtnpHDUu3V7O5vL7TrxcMR5v9V0RERI5cCmkHIRCOErUQiVrAqZZd++hifvL88sQ5s9aV8dPnVxCIBatQKwGroj7A/C172n29eDir9Ye4+uFFbNhddyi+DRERETkMKaQdBH/ImRsWr6ZtLW8AoLZpbwt09sYK/rO0JBGwApGWIe0Lf5/DFTPnt/t6wfjrVDTwwboyFm6t7NS41+6s5YK7P6KmKdSp54uIiEjXU0g7CP6wE9Li4alojxPSBvdJS5wTCEcIRy2NQSe4tdaqLK5sAiDcSoADp9K2ubw+8dw6v3Ot+DUP1NqdtWzYXc/2PY2der6IiIh0PU9PD6A3C4SatzDj8836pPsS5+xtUbYd0hLnRqJ43C1z813vrGdRURXnH9PfuVasAtYQOLBVnmtKa9lcXp9ovTZ0MuSJiIhI11NI64Qbn1xK/+zURCUtFHHmpG2OtTuTFxME9ql+7TekhaMk5buE8rogNU2hg66kPTqviHfX7OamGaMBaEhamfrh+jL6ZPg4bnDuAV1TREREuobanZ2wvbKRTWX1+OOVtFgo2xKrpAWTQlryZP99H9tXWwGuIRAmEIq0uFbDAe6X1hiM0BSKJIJj8vYh//vaGv72waYDup6IiIh0HYW0TkjzumkKRgiE9y4caApG2FHtzC0LJIWteCiLtyj3DWLJQSnQRkhrDIYJhKOJMFgXC2mNHdyjLc4fC2iJdmdSu7QhEEmEPxEREel5and2QrrPTXl9IKmSZimuasQ6XU+C4SjPLykhM8WTCHJtzUkrS9qYtq2QVh9wQlqikhZbPXqglbSmUIRI1CbanMntzsZguMMb84qIiEjXUyWtE9JTPDQGI8224IjPEwMniP3rk608sWBbUrDa2+58c+VOyusCAOyuDTR7Xtwzi7ZTXOmsvoxXvOIhqi4Qq6QFw8zdVMGmso7tlxZf6LC3Xbp3zE2hCHX+MNv2NPDh+rIOXU9ERES6jkJaJ6Qn2p2xOwDE2p1xgXCUQChCILS3tRgPRrX+EN97cinPLNoOQFldciVtbxj76X9W8tySEmBvxSseBGuSVnf++PkV/OX9lnPJVu2owcZLezHxhQ57n+9cLxSJEopY6v1hZn68hZue+rRzb4yIiIgcMgppnZDuc9MQCCcqX6FwNFGVSve5CYaj+EMR/OFIixZlbVMIa2FP7NZRu5PanfFz48fq/WGstYlrxytp8Ws1BsPsaQi0uA3Vqh01XPLXT1psdhsPkvHn18cqdI2x43X+MFWNQer8Yd16SkREpIcppHVCms+TmGMGzpy0eADKTfMSjETxh6PNqm2JSlosINU0Ol83a3dGmoe0hkCYplCE2F2n9rY7Y9eqbAjiD0VbhLTS2AKG4qqmZsfbqqTFxx6MRBPjqW7q/I3gRURE5OAppHVCus/d7OtQJJqoRuWk+1qtpO27fUZ1LCjtrNkbpOLnxOer1e8zmT8ezuKhraLeCVLVjUFeXraDX7y0qtm1k1upAE3BfeakBVrutxafB1fVoJWeIiIiPUkhrRNaD2lO0MlL9xIIO4sK/Elz0uLiCwiqGp2AtXF3PUNit5EKhFtW0pK3yUiu3iWragxxy9PLeGz+NraU1yeqdPGwFxcI7VNJC8ZD2t7XKIs9Jz4+ERER6RkKaZ2Q1iKk2UTQyU330hBwWpTOvmTNt8mIB6SaRucOAlsrGjh2UA6wt5JWFms51vvDzbbJaGueWFMoQn6Gc6uCl5aVJl6jbJ+QFm937ntbqaZQy608qhXSREREepRCWie01e70uV2k+zyJdmIgFG0RrOLVsOqmEFsrGghHLRMHNg9pu+PtzkDzkLY/8ZWfL326I1EFi1fSnlq4nVnryxK3r4q3S/e2O1uGtD0NQZZur+rQa4uIiMihp81sOyHN2/xtc7bgCJPmc+PzuPbep7OVW0DFq1jVjUHW7aoFSFTS4lW3RLszGO7wTdCDkSgel2F7ZSPbY/PKKmIh7c/vbmBc/6wWz6lPLBxo+RrPLipmeUkNb//gTMYVtnyuiIiIdC1V0johI6VlJa0hGCHd58bn3v9bGo7urWZ9ur0at8swfoATggL7LBxoCEQS22R0xPFDnJujr9/lbG5bVhcgFIlSUR9oMT/NuX7blbSVO2oA2FrR0OHXFxERkUNHIa0TWrQ7Y9ttpPncpHg6/pYu2FrJsPx0slO9gFMNs9bu3SftANqdAMcNdkJaWVK7dNueBqxtudIToDEUYXN5faK6lyzeEo3fj1RERES6l9qdnbBvu9NZOBAmw+c5oJC2blctF00qTFTfguEo9YEwjcEIWSke6gJhqhv3vxWGx2US1bnJQ3ISx90uQyRqExWxqlauYy2ce9dHjO2f2ex4/LkAO6oU0kRERHqCKmmdsG8lLRhbOBCfk9ZR1sLEgTm4XAaPyxAIRymJhaLRseCUfEeCZGleZwxD+qQnjk0alIMxzufD853jK0pqWn2+x2USn28qqwdIhMUheWmJx3ZUN3b4+xEREZFDRyGtE/YNaeHYFhzpBxjSAI4ZmA1AisdFMBxlVazydfKIPkDrbUqAPrEtN/pmppCZ4lT2BuWmkZ+RAsCYfs48t5VthLS82PPBaW36PC6y05y264iCjMRjaneKiIj0DIW0Tmi5T1o00e5sb+HAviYOcEKaLxbSVpTUkJniSaz43F0baLWFmpfhBKqcdC95GV5y0ryket0U5jgh7ZiB2RjTdiUtPymkgRM8s1I9eFyGwXlOFa5fVoranSIiIj1EIa0T0n37zklLWjjgdbfxrL28bqfVWJCZQr/sVCAppO2oYdKgbLJiiwl21/pbBCqAvHTnWG6al7x0H/2znXBWmO20Kof2SeeYAdmtbgMCtAh+6V4npOWme8lLd1779NEFVDWGmt02SkRERLqHQlonuF0m0db0eVwED2ALDoCMFA+ZKR4mxlqdACkeN/XBMGt31nLc4FwyY9t8lNUGyE33JeaaxcXbnbnpXi6cVMglxw0ESFTSctK8nDa6oM0xTIhV8AoynfPTYpW0nDQvU4blMWlQNtNH5QNaPCAiItITtLqzk9J9boLhKJkpnkQlLd3n2e+ctPhKzFSPmy+eMiixZQY4YW/1jhqC4SjHDsohM2XvthyZKc6qUX8oSqrX+W9mioerpg/j/GMKE/PXAAbkOJW0nHQvp47KZ+bHW8hN9yZWiRrjLFg4ZmA2m393MX9+dz33zdpMus/DF6YMproxyDnj+nHOuH4sKqoEnHlpY1rZDFdERES6jkJaJ6V73TR5nFtBNQWjBCPRNhcOxLfTyEz1UN0YItXr4sefGd/sHJ/bxeZyZ5XliIKMZhvmjijIYN2uWvyhKFmpXvyhACkeN/9z6TEtXmt4fgbGQP/sVMYXZpHmdXPsoBxmb6xIjKXWHybF48LtMvRNqqR96cTBza4Vb6nWtLKPmoiIiHStHml3GmMuN8asNsZEjTFT93PehcaY9caYTcaYn3XnGNuTnuIh1evG4zaJENNWSIuvmsxKdTJxaivz1nweFza2gWxehi+xYhPg+KG5iblu8WukeFv/o7twUiGv33QGg3LTSPd5eO3m07n1/LGJx3NjwSvF41yvb1ZqYuwtxh17rfj9RkVERKT79NSctFXAZcDHbZ1gjHED9wEXAccAXzXGtCwd9ZD02N0FfG5XYsf+9DY2s40Hq6xYC7O1kJb8vLx0LxnJIW1IbuLxrNjxtjbNdbtMYlsPgFF9MxMtUHDmsCU/v29WSuL7aTlu59x6hTQREZFu1yPtTmvtWgCz72z45k4GNllrt8TOfRr4PLCmywfYAWleN6leN163q1klLR5+knftj9/2KTNRSWsZsJIXIqR53c3em7H9sxLXzdxPNa4t8ZDoc7sSm+DGK3H9YiFt37soxMfpdhnq/Gp3ioiIdLfDeXXnIKA46euS2LFWGWOuM8YsNsYsLi8v7/LBxQOZ17O33Znmc+NzOyEo3ioEyE7zNDu2v0pan3Rfi/DqrCZ1npPZTiWtrbG6XYYUryvx2nvbnW1X0owxZKV6qFMlTUREpNt1WSXNGPMeUNjKQ7dba18+1K9nrZ0JzASYOnWqPdTX39c3Tx1OVWOQR+dtS4S0jKTVnWleN42eCIFwNNE2jP831dP6nDTY244EuHnG6MRWGYl2Z2q8XdnxSlo8bHndrkQVL369jBQPg3LTGNInrdXnZqV6qD+Am7yLiIjIodFlIc1ae95BXmIHMCTp68GxY4eFs8f1A+CphcU0hSIAze7dmRprhwbC0UQFLV4Fa63dGQ9d8RWVALdeMC7pcVezaxxIJQ2csGUwiddJDnnv3XpWm1uHZKV41e4UERHpAYdzu3MRMMYYM8IY4wOuAF7p4TG1kLx5bXbq3oUDKV53Iox1aHVn7Drx2z3tq6OrO9uSneol1ZtUSUt6flqsHdqazFSPVneKiIj0gJ7aguMLxpgSYDrwujHm7djxgcaYNwCstWHgRuBtYC3wrLV2dU+Md3887r3hpk+GL6mS5sz/cpm9t5Ha36T/ve3OlreAgpaVtNZapvuTneolzetOqqR17I8+W3PSREREekRPre58EXixleOlwMVJX78BvNGNQztg3lgFzGWcgNUQu89lqsedCEUtts/Yz+rOPm2ENN++c9IOsJJ207mjCYajzNnkbGrb0TltWale6gN1B/RaIiIicvB0x4GDlGhTpvucVZjuvZW0FK8zR+2EYXmcMaYgsV9Za1WwlFYWDrT2+LGDcjhhaC7jDvA2TaeOcu7juWRbVbPrtUerO0VERHrG4TwnrVfwxtqd8RuepyQvHPC4SPG4OH5ILo9dc0pig9r9tTvz2mx3Os8ZlJfGC987jX7ZqZ0ab+o++6S1JzPFCWnWdvmCWREREUmiStpBirc78zOdcGWMU02Lr+5MXjW5d3uOttudbS4cSNrs9mBcdsIg+mWlJObJtScr1UskavGHoqS1speaiIiIdA1V0g6SNxaa8jNSEsd8HmcVZbrP3axqllxl21drW3A0ezwW7Lzu/d6loV0DctK4fOqQ9k+Mia8m1TYcIiIi3UuVtIMUn4MWb3cCnDO+H1OH9WFcYRZVjcHE8fj2F61VpNprd+am+Ujzuptt+dEdspJust4vu52TRURE5JBRSDtE4u1OgL9+dUqr54wsyOT2iycwY3y/Fo+dO74fO88exdA+6a0+9xvTh3HO+L7t3e/0kIuHNN11QEREpHsppB2keHjJz2i9ApbM5TJce+bIVh8bmJvGTy4c3+ZzM1M8jC/s/lJWfMsPtTtFRES6l+akHaR4eOmTNCftSLJ3TpoqaSIiIt1JIe0gxcNLTlrrqzJ7u/gdDlRJExER6V4KaQepNhZe4hWnI83edmeYP7+znsfmFfXoeERERI4WCmkHaeKAHAAKczq3uezhbm8lLcyrK3byzprdPTwiERGRo8ORWf7pRr/6/ES+MX0Y/Tt5B4DDndtlyPC5qfOHqfOHqG3Sj4yIiEh3UCXtIKV63UwalNPTw+hSzk3WQ9T6w1pAICIi0k0U0qRdWake9tQHCYajiTl4IiIi0rUU0qRdWakedlQ3AVDbpEqaiIhId1BIk3ZlpnopjYW0YCSKPxTp4RGJiIgc+RTSpF1ZqR5qk+aiqeUpIiLS9RTSpF3Z++wBp8UDIiIiXU8hTdoV39A2rrZJlTQREZGuppAm7YpvaBtXq0qaiIhIl1NIk3bte8sr3cdTRESk6ymkSbtatjtVSRMREelqCmnSrni70+MygFZ3ioiIdAeFNGlXfHVnv6wU3C6jdqeIiEg30N2ypV3xdmdWqpfGUETtThERkW6gkCbtii8cyEr10BTyqpImIiLSDRTSpF2ZzUJaRFtwiIiIdAOFNGnX3kqaF38oqs1sRUREuoEWDki7UjxufB4XWakeslI9ui2UiIhIN1BIkw657fyxXHbCYNJ9bppCkZ4ejoiIyBFPIU065PqzRnHisDzSfB4agwppIiIiXU0hTQ5ImtdNU1DtThERka6mkCYHJN3npjEUwVrb00MRERE5oimkyQFJ87mxFgLhaE8PRURE5IimkCYHJN3nBqBJ89JERES6lEKaHJB4SGvUCk8REZEupZAmByTN52xsq8UDIiIiXUshTQ5ImjdWSVO7U0REpEsppMkBSbQ7FdJERES6lEKaHJC0+MIBzUkTERHpUgppckC0ulNERKR79EhIM8ZcboxZbYyJGmOmtnHOEGPMLGPMmti5t3T3OKWldK+zcEDtThERka7VU5W0VcBlwMf7OScM3GatPQaYBnzfGHNMdwxO2pbqc35ktLpTRESka3l64kWttWsBjDH7O2cnsDP2eZ0xZi0wCFjTHWOU1qX7VEkTERHpDr1iTpoxZjgwBViwn3OuM8YsNsYsLi8v766hHXXiW3Bo4YCIiEjX6rJKmjHmPaCwlYdut9a+fADXyQT+A/zAWlvb1nnW2pnATICpU6fq7t9dxO0ypHhcWjggIiLSxbospFlrzzvYaxhjvDgB7Qlr7QsHPyo5FNJ9brU7RUREuthh2+40zoS1fwFrrbV/7unxyF7pPo9CmoiISBfrqS04vmCMKQGmA68bY96OHR9ojHkjdtppwDeAGcaYZbGPi3tivNJcqtdFU0irO0VERLpST63ufBF4sZXjpcDFsc8/Adpe/ik9Jt3n0Zw0ERGRLnbYtjvl8JXWzpy0VTtqmLOpohtHJCIicuRRSJMDlu5z73cLjnve28AdL63qxhGJiIgceRTS5IC1t7qzrC5AZUOwG0ckIiJy5FFIkwOW6nXvd05aeV2AWn+ISFTb1YmIiHSWQpocMKeS1vrqzmjUUlEfwFqoaQp188hERESOHAppcsBGFmRS1Rji9RU7+enzK2gI7A1sNU0hQhGnglbVqJaniIhIZymkyQG7fOpgslI8fP/JpTyzuJiFRZWJx8rrA4nPqxXSREREOk0hTQ5YVqqXr00blvh6c1l94vOKur0hrapB7U4REZHOUkiTTvnRBWN594dnkp/hY1NZfWKOWrNKmuakiYiIdFqP3HFAej+P28WY/lmM6pfJoqJKpv3ufc4/ppBxhZmJc9TuFBER6TyFNDkoY/plsnCrMyftP0tLAPB5XESiVgsHREREDoLanXJQxvRzKmfD8tM5b0I/ADJ8bnLTvFQ1qt0pIiLSWQppclBG98sC4LPHDuCGs0cDUNUYIjfd26l2Z1mdn2A4ekjHKCIi0hsppMlBOXFYHpefOJivTxvGicPymDwkl6tPG0Feuu+AV3eGI1HOvesjnlywrYtGKyIi0ntoTpoclDSfmz9ePjnx9cvfPw2A7zyymJKqxgO6Vk1TiDp/mNIa/yEdo4iISG+kSpp0ibx07wHfFqrW72zjUedv/ZZTIiIiRxOFNOkSfbNSKK8LUJG0b1p7amOhrj6gkCYiIqKQJl3iiycOJmotMz/ewubyek7833eZu7mCP729nvfX7gagIRBmU9LdCuKVtwaFNBEREc1Jk64xqm8mnz9+EI/OK2LJtir2NAS5+90NLCqq4pQRfTh3Qn/+8v5GHplXxKLbzyMr1UutP1ZJU7tTREREIU26zk8vHM/S7VUs2VZFZoqHRUVVAHxaXI0/FGHO5gr8oSjvrN5NYyiCtRaAOlXSREREOhbSjDF3Af+21q7u4vHIEaQwJ5VnrpvOc4uLGdE3gxuf/JTMFA/1gTAfbyhnTWktALc9txyAaSP7AFAf0Ca4IiIiHZ2TthaYaYxZYIz5rjEmpysHJUeOwpxUbjp3DBccU8jJw/vw689PxGXgvlmbiFoYmJOaOHdrRQOgdqeIiAh0MKRZax+01p4GXAUMB1YYY540xpzTlYOTI4fP4+LZ707nshMGM3V4H5aX1OBxGX5xyTHkZ/gA2F3rrAStD4QTrU8REZGjVYdXdxpj3MD42EcFsBy41RjzdBeNTY5Q9115Al+ZOoSrpg/nomMHsOQX55OdurfzHopYAro1lIiIHOU6OiftbuBS4H3gd9bahbGH/p8xZn1XDU6OTH2zUvh/Xzqu2bH+2anU+vdux1EfCJPqdXf30ERERA4bHa2krQAmW2uvTwpocScf4jHJUagwaW4aaF6aiIhIR7fgWA6MM8YkH6sBtm+ZykYAACAASURBVFlraw75qOSo0z/bCWmpXhf+UFR3HRARkaNeR0Pa34ETcCpqBpgErAZyjDE3WGvf6aLxyVGif3YKAANz09hS3qD7d4qIyFGvo+3OUmCKtXaqtfZEYAqwBTgf+ENXDU6OHoWxStqg3DQASqubqPWHeHnZDi6+dzbhiBYSiIjI0aWjlbSxyRvZWmvXGGPGW2u37NMCFemUfrGQNjDHCWm3Pbec81b1Z0ifNNbsrGXdrjomDdL2fCIicvToaCVtjTHmH8aYs2Iff48dSwG0PbwctHglbWCskgZQXNlIVUMQgEVFlT0yLhERkZ7S0ZD2TWAT8IPYxxbgWzgBTRvaykEb0z+TM8f25dwJ/RLH9jQEqWp0fgdYHLvvp4iIyNGi3XZnbBPbN6y15wB3tXJKfSvHRA5Ius/Do1ef3OxOA1WNQSqTKmnWWtReFxGRo0W7lTRrbQSI6n6d0h2SQ1gkatle2QhAWV2AXbX+nhqWiIhIt+vowoF6YKUx5l2gIX7QWntzl4xKjmqnjy5gS3k9pTV+appCFGansqvW+XxATlr7FxARETkCdDSkvRD7EOlyj3/nFGatL+PbDy0CYEifNHbV+nUXAhEROap0KKRZax8xxqQBQ621ulendLn8DF/i88F56SwqqtJdCERE5KjSodWdxphLgWXAW7GvjzfGvNKVA5OjW1763pA2JM9pcTYEIj01HBERkW7X0S047sS5kXo1gLV2GTCyi8YkQn5mUiWtTzoA9QFtySciIkePjoa0UCs3Utd9eqTLpHndpHicH88hefGQpkqaiIgcPToa0lYbY64E3MaYMcaYvwJzO/uixpjLjTGrjTFRY8zUds51G2M+Nca81tnXk97HGJOYlzakT7zdqTlpIiJy9OhoSLsJmAgEgKeAWpw7D3TWKuAy4OMOnHsLsPYgXkt6qbxYSCvITCHF49LCAREROap0dHVnI3B77OOgWWvXAu3uHm+MGQx8FvgtcOuheG3pPfpk+Ej3uUn1uslK9SikiYjIUaVDIc0YMxb4ETA8+TnW2hldM6yEe4CfAFld/DpyGBqYk8aO7CYAMlI82idNRESOKh3dzPY54J/Ag0CHZm8bY94DClt56HZr7csdeP4lQJm1dokx5uwOnH8dcB3A0KFDOzJEOcz9+MJx1DY5KzozUzyakyYiIkeVjoa0sLX2HwdyYWvteZ0YT7LTgM8ZYy4GUoFsY8zj1tqvt/F6M4GZAFOnTrWtnSO9S0FmCgWZKUCskqaQJiIiR5GOLhx41RjzPWPMAGNMn/hHVw7MWvtza+1ga+1w4Argg7YCmhz5MvcJaS8sLWFxUWUPjkhERKRrdTSkfRP4Mc62G0tiH4s7+6LGmC8YY0qA6cDrxpi3Y8cHGmPe6Ox15ciV3O4MR6Lc/uIqHppT1LODEhER6UIdXd054lC+qLX2ReDFVo6XAhe3cvxD4MNDOQbpXZLbnZvK62kKRahuCvbwqERERLrOfitpxpifJH1++T6P/a6rBiWyr+QtOJYXVwNQ3ajbRImIyJGrvXbnFUmf/3yfxy48xGMRaVOGz4M/FCUcibKs2LlDmUKaiIgcydoLaaaNz1v7WqTLZKS4AWgIRJIqaWp3iojIkau9kGbb+Ly1r0W6TFaqM32yvD7A+t11pHhcNAQjBMPRNp8TCEeIRPVjKiIivVN7IW2yMabWGFMHHBf7PP71sd0wPhHAWTgAsGDrHiJRy7SR+QD7XTzw9QcX8NvXddtXERHpnfYb0qy1bmtttrU2y1rriX0e/9rbXYMUyYyFtDmbKgA4c2xfAGramJdmrWVFSQ3b9jR0zwBFREQOsY7ukybSo0YUZADw7prdDMxJZWz/TACqm1oPaeV1AQLhKA1B3aVARER6J4U06RWG5Wdw2uh8QhHL5CG55KX7AKhqaL3dWVzVCDgLDURERHojhTTpNb4xbTgAxw3OJSfN6ba3VUkrrmwC0E3ZRUSk1+roDdZFetz5x/TnF5ccw38dPxCfx/n9oq05acWVsUqa2p0iItJLKaRJr+F2Ga453blDmbUWj8tQ1cZeaWp3iohIb6d2p/RKxhhy073N2p1ltX6ue3QxO2ua9rY7g2Gs1V5pIiLS+yikSa+Vk+bl1eWlzFpXBsCDn2zlnTW7eXphcaKSZi00hVRNExGR3kchTXotr9tFnT/M9Y8vYVeNn6cWbAfgiQXb2FHdRGF2KkDixuwiIiK9iUKa9Fq/+txEvn/OKILhKNc+upi6QJivnTKUivogfTNTuOHsUYDmpYmISO+khQPSa50yMp+ThvfhpU9LWbmjhvOP6c9PLhzPsuJqbj1/LOHYfTu1DYeIiPRGqqRJr+ZyGT53/EA8LsNPLxxPTpqX128+g3Mn9E/cSkohTUREeiNV0qTXu+XcMXzpxMGM6pvZ7Hj8puzaK01ERHojVdKk10v1ulsENIAMnxvQnDQREemdFNLkiJWhdqeIiPRiCmlyxIqHNG3BISIivZFCmhyx4u3OxmDLdmckanUnAhEROawppMkRy+N2keJxtWh31vlDTPn1O7yzZncPjUxERKR9CmlyRMtM8bRod24pb6DWH2ZNaW0PjUpERKR9CmlyRMtI8bRod26vdO7rWV4f6IkhiYiIdIhCmhzR0n3uFpW0REirU0gTEZHDl0KaHNEyUzwt5qQVK6SJiEgvoJAmR7Q+GT521fqbHVMlTUREegOFNDmiHT80ly3lDVQ1BBPHkuektbcNxx0vreRdrQIVEZEeoJAmR7Spw/oAsGRbFQChSJTS6iYyUzwEw1Fq/W1vdFvnD/H4/O28vGxHt4xVREQkmUKaHNGOG5yD121YtK0SgB1VTUQtTBmaC0B5ndMKfXD2Fk77/Qe8srw08dxNZfUAbC5vSBxbXlzNjU8uJRSJYq3l8fnbKK1u6q5vR0REjiIKaXJES/W6mTQoh/lbKrHW8p+lJQCcN6E/AGWxeWmz1pexo7qJW57+lLLYHLaNu52QtqW8nmjUaYs+vaiY11bsZO3OWkqqmrjjpVU8Mq+oe78pERE5KiikyRHvwomFLC+u5vrHlvDg7K1cctwAThudD8CVDyzgT2+vZ01pLeMLs7B2b2t0Y1kdAIFwlB2xatm8zRWAU1Fbs9PZDHfVjhqKKxvZVePf96VFREQ6zdPTAxDpatedOZKGQJgHZm8lxevitgvG0Sfdl3j8oTlbaQhGuHHGGP7fW+t4Y9UunllczPbKRjwuQzhqWVRUyc4aP0V7nEUHy4prGJznLEZYUVLD1x5cQGMwwis3nsbA3LQ2xxIMRwlHo6T7Do+/etZaQhGLz6Pf10REDjf6P7Mc8Ywx3HrBONb+74Us+58LGFGQQXba3pDUELsjwfFDcpg8OIdXl5fy4fpytpQ3MH2UU3G79dnlfPn+eQAM6ZPG8pK9lbQ6f5jtlY1U1Af46gPzmbd5T5tj+dFzy7n43tn4Qy1v+h4XjkRb7O3WlqKKBp5dXEwgHOE3r61hc3l9s8eD4SjFlY3U+kPc8dJKdtbsnT9X1RDkipnzufDejw/pzeattXywbjfBcDTxdSR66K5f1RBk/pa232MRkSOFQpoclYwxvH7z6Xxw21mxr2FcYTYnxlaDjijIAGDayPzEcyYPzmFU3wy+eMJgNpfXM3/zHkb3ywTA53bx0LdOAuCbDy2krNbP+l11XP7PuSzd7rRPt+1p4NUVpRTtaeSfH23GWktNU6jZuHZUN/H5++Zw4m/e5ecvrGDe5j28+GkJ9763MbGv286aJmoaQ7y7Zjfn/vkjfvL8Cn7+wkoe/GQrtz67vFkgemD2Fmbc9SF/eGsdj8/fzk+eX5EIZLc8s4wFWyvZUt7QbHFEe8KR6H4fX15Sw9UPL+bReUUUVzZywd0fc9uzy5qd01oo7EhQjEQt1z22mCtmzmdTrB19oB6fv40z/vAB5/zpQz5Y17Pbq3y6vYrlxdX7ffzfn2xN/Ay15VCG7J5yJHwPvcmqHTUtfqk7GlQ2BKnzh9o/8TBxePRcRHrAxIE5AEwalE1DIEJmioezx/XloTlb+fvXTmDbngamjcxnT32QNJ+LH39mPNGoZVlJNfe8t5G6QJjrjx/I32Zt4vTRBZwzvh/DCzKYcdeH/O/ra5m/ZQ/ldQFuf3EVr9x4Gve+txGPy3DKiHzueW8jD88toqYpxAs3nMrO2Hy2X7+6hoZAmAsnFvLC0h08tbA4Md7iqkaqGoK8v66MqcPyyM/00S8rBWvhhaU78HlcLC+u5tF5RXz7tBEAvL16F6GI5fH528lK8TB7YwX3zdrEOeP78fGGcr48dTDPLi7h/o82s2ZnLX/56hRG9c3koTlbyUnz4jKGPQ1BJg/OISPFw+yN5dz1zgZ+eP5YvnvWKD5cX0ZGioeThvchErW4XYaFW50q17OLi3l4bhElVU1s29PIbRc0srWigUA4yu0vruSBq6YyeUguDYEw33lkMfmZPv525Qn7/TN7aM5WFhU5geUPb60nM8XDzy4aT7/sVOc1FxVTFwhzzekjqGoI8vSiYr42bSjZqV4AVpRU88tXVjN5cA71gTDXPLKY/9xwKicMzWP7nkYG5aXhdplD90OWxB+K8NqKnVxy3ABSvW6qGoJc9e+FuIzhox+fzQfrynhz1S7+9/OTqGwI8vDcrTy72Fno4nO7ePw7p3DyiD7Nrrmrxs/VDy/ijDEF/PziCYnjoUiU5xaXUB8I8ZWpQwlFo6R53WSktP6//FAkyouf7mBYn3ROSfrFpCNKqhq55ell/PC8sZw+pqDF49ZajHHe07dW7WJPQwB/KMrxQ3I5cVgee+oD3PnqGt5fu5v7rjyBc8b3O6DXP1SWbKviqYXbGZyXxs0zxuBq4+dgZ00TG3fXc+bYvi0ei9/NZEif9C4bZ0V9gHSf+6CmTATDUb7574X4PC7euuVMUn0uUjzuQzjKrvHX9zfiD0f48WfGd/g5oUiU5cXVnDA0D5fLcOUD8+mblcJj15zShSM9dMyR+NvL1KlT7eLFi3t6GNJLrCmtxR+OcMLQPMD5S+1177/I/MG63dzz3kbuunwyu2sDDMtPT/yP+YbHl/Dmql0UZqfytVOGcte7G8hL91LVGOLaM0bwg/PG8vSiYlbvqOGNVTs5dlBOInjkpXt58tppTBiQTWVDkCXbqhien87Dc4t4YsF2ACYPyWV5cTXpPjeXHDeAQbnp3P3eBm48ZzSrS2uYs3kP00fmk+5z8+aqXWSneqj1h/nzlyfz0YZyXl5WSv/sFOr9Yeb+7FzOu/ujRJVu+sh8br1gLJf/c16b3/ug3DR2VDfxy0uP4a53NpCf6ePyEwdz93sbOXFoHm6XYV5SO/JHF4zlT+9sID/Dx56GIG6XIRK1TBvZh6eunca1jy7hvbVOReuF752a+HPY16ayOi7+yyecOaYv6T53YruUm88dw63nj2X9rjo++5fZhKOWe684nk+3V/Pw3CLG9Mvkmeun0xSKcOUD8/GHIrzzw7Nwuwzn3fUROWlepg7P44kF27n8xMH88fLJ7f7MVDc638ery3fy/trd/O3KE0jztfxHzh+KEI5aMlM83P3uBu59fyM3zRjN9WeN4s5XVvNCbLVxfmZK4s8gJ81LTVMIn9vFVdOH8fVpw7j6kUXUNoV4/7azyUlzAmcoEuWie2ezqayerBQPi+44j1SvM4b7Zm3ij2+vB2DG+H6sKKlhUG4qL3zvtBYhNByJ8sV/zGV5SQ256V4+uO1sVpRU8/ySEv50+eTENdty7aOLeXfNbvLSvbxxyxkMyNk7J/OtVbv41aurefCbU7EWLvnrJ4nHCjJ9fPyTc/jpf1by9qpd9M1KoT4Q5rtnjeKLJw6itNrPW6t2MSw/na+ePLTZa+6obuJ7TyzlmAFZXDp5IOt31fHN6cNbBKunF25nzc5azpvQPxGqrLW8vXoXnxZX85mJhYmft6/OnM/CokoiUcvMb5zIBRMLm13rleWlPDRnK6tLawmGo7x20+lMGpTDK8tLmTQwm/s/2sIzi4sZmJPKJz+d0WbISxb/xSZZVUOQiLUUZKa0OL+szs9n7v6Y4wbn8sjVJyeOVzYEueXpT1m7s5Y/Xj6Zc8btP+i+tqKUG5/8FACfx4XbGL4xfRjj+mdx1ri+7Kz2M64wKzFXdXetn1AkyuC8Qxc+A+EI768t46yxfRO/PGwqq+eJBdvwuV387KLxLNlWxYnD8jDGsKO6ibP+MAuABf99Lvmx9ycStdz/8WY+M7GQUX0zm73G1ooGvv7gAnZUN/GHLx3H9JH5nBG7xvu3ndXi/GTFlY34QxHG9M86ZN/z/hhjllhrp+573H3nnXd2ywC608yZM++87rrrenoY0kv0zUpp9g9LRyopIwoy+erJQ8nPTGFofnriH06A4wblkuJ18+evHM9ZY/uSk+Ylzevm2jNGct2ZI/F53JwwNI/PTCpk/e56Zq0vJ93n5k+XT+bGGWMYG/ufQprPzai+meRnpjBhQDaPz9/G+RP68z+XHMPj87cRiliuO3Mkl04eSEV9kO/PGMVnJhbywtIStlY0sHan0w7897dOYmRBBl89eRgXTiwkYi3VTSGumj6cU0cXsKKkmg276zl2UA5Ltlfx1spdZKd5+etXp3DFSUO47YJxnD6mgM9MLOTr04Zx+2cnMH9LJU8u3E4wEqWmKcTibVUMz89g5Y4aSqqaOGdcX0pr/Fxx0hB+/JlxPDp/G1WNocT78fnjB/LyslKaQhGeW1LCzTNGs2F3HWtL6xjTP4v/e3Mtd7y4ilF9M5n58RbcbsNvXl9LYyDMw1efxOTBuVTUB0jzullWXMU3pg3jhieW0hAMc8yAbJ5YsJ1VO2qYMiSXVaW1uF2G/3tjHVWNQR646iRGFGTg87gYlp/BEwu3s6y4mmMH5TBrfTlvrdpFnT/Mp9uruOHxpcxaX8aY/pn8+rU1uI2huLKRKx+Yz0Nzinhz1U62VDRQVNFAqtfFwFxnvuKe+iAFmSnc+uwy/uflVfTLSuXudzdgsSzZVs2/PtnK8pIavnnqcI4fmktRRSM/vXA83zl9BIuKKvn2qcP525VTuGBiIbnpPk4clse/5mylKRjh7Ng/wK8uL+XJhdv5+rShLN5WxbGDcxjdL5PiykZueupTZozvx2ePHcDTi4ppDEbYXRugMDuVYwfnUFTRwP97az2vLi9lVWktr63YyY3njGbu5j0s3lbJI3OLWF1ay6h+GUwYkJ342a5uDBKMRBNVl7dX7+Le9zdy5SlDWV1ay4bd9Xxu8kCMMUSjlhueWMr2ykY+3lDO8uIaahpDvPWDMzhnfD+eXFjM+l11vLVqFzfOGM0dl0zgzVW7eH3lTmatL+fhuUXM37qH99aWcd6E/olq6aodNVz5wHy272lkRUkNL3y6gw/Xl7O5vJ5zJ/TD43KxsqSGLeUNfO/JpawpreXV5TspyErh3dW7KKlq4pZnlrF0exXPLS5me2UTGT4Pf521ke+fM5pdNX5WltZyxUlDMMawYXcdO2v8XPfoErweF589dgBrd9XhMjAoL40rH1jAkwu3s2JHDScNz2P9bmcc/WPjbY21lnvf38j1jy3hxKF5ifATjkT5wt/ncO/7G+mbmcLEQTnNnnfL08tYs7OWbXsauXBSYSLI/XvOVp5aWIy1YIALJhZS3RikPhButeJ25yuriVq44exR+NxuRhSk8+ziEt5Zs5uZH2/hqYXbqWoMMmN8f+54aSU3P7WMpxcVc/nUwW1WY/2hCNsqG7HWcsdLq/jL+xu5650NVDWGmDI0l39+tJk3V+3ilBH5fLi+nGsfXcyTC7ZTHwhzzvh+NAUjXP7PeczeWM7ibVU0BsP87IWVjCjIoGhPA3e/u4HiqiYiUUu/rBROGOaE69dX7OSOl1Yxa30Z/3X8oGa/LP3+jXV8WlxFYU4q2/Y0kp7i5oN1ZRgDHpeLoX3SeXhOEScMy2VntZ+sFA/GGJZsq+Ty++fx6opSrjltRKIS3JV+9atf7bzzzjtn7ntclTSRHvT+2t1c88hivnXqcO783MT9nlta3UTfrBQ8LsNZf/yQ7ZWNLPzvcxP/eMXVNIbwegz//cJKVpXW8u4Pz9zv/2ReX7GT37+1lle+fzrPLi7mb7M2ccdnJ/CVk4a2+ZyVJTVc+rdPmDQom01l9fhDUZ66dhr/+mQL760t4w9fPI5TR+dTmJ2Kx+3iFy+tYlFRJa/ceDo+j4tAOMJlf5/L6tJaCjJ9zP7JDF5bUcpP/7OCqIWsFA8ZKZ7EfVfj1be7Lp/MF08cnBjH80tK+NFzy5k+Mp95W/Zw7xXHc+aYvnztwQWs3VXLOz84k1+/toZPNlVgLTx2zcmcMaZ5m8ofihCMRMnwefjt62tZVFTJyh01AJw4LI81pbU07bPQY0RBBrnpXkKRKKeOKmDmx1sApyoRXzBx9ri+zN5YgcuQWEH7wFVT+ePb65g4IIcvnzSEE4bmdvgfgF+8tIrHF2zjphljqPOHmLWuDI/bxZu3nMG0373P2P5Z/PJzx3DD40upqA/wxs1n0CfDxxf/MZfPHjuATzZVsLykmuvOHMXDc7YSjlVx6vxhJg7M5rWbTueRuUX89YNN+DwuPG5DVoqXs8f15dunjSDV6+LCe2ZTHwhz2wVjOWZANtc/toT+2am89P3TeHReEb95fS3fOX0E00fls6a0lrve3cB1Z47k+SUlVDYE+d7Zo/jJhU6r6lsPLeTD9eVMHpzD09dNT/zj+tGGcr790EL6ZPh46tppfPn+eWSmeijMTmXCgGyeW1xCnwwf//rWVP4+azNVjUFOGt6HP7+7ganD8pg0KIeH5xYBzi9gz393Ol/8x1wq6p3V2MbApIE5PH7NKdzz/gaeiYVYgA9uO4v5Wyr57xdXcunkgZTX+Zm/xdkI2+d28d6tZzE0P52bnvqU2RvL+cpJQ3jg4y2MK8xm+sh8bpwxmqm/eZcbZ4zhxGF5/N8ba8lK9XDlKUN5eE4R00bmMzgvjZmzt1Bc2US6z43bZfjBeWM5eXgf5m6u4P/eXMeYfplsLKvn5nPHMHFgNm+t2sW5E/px45Of8t2zRvHI3CJOG13AfV+bQjhiueje2QzMTSUnzcvq0lr+/OXjueHxJRgDM6+aytJtVczfUsnGsjp+eN5YfvDMMn50wVhunDEm8fPVFIywqayeWevLWLuzlrdW7+LGc0bz1w82cclxA3hz1S4umzKIzx8/iGkj+3D/x1s4Y0wBxw3OTfx8PrlwO5MH57BqRy2njs7HZQwfrCtr9nMc/zsyul8mQ/LSmL2xggeumsqrK0p5YekOHrhqKjc/9Wni71y8cg/wrVOHs6y4mj0NAZ69fjqF2alc+rdPqKgLUtkY5OJJhdxzxRSWbKvkow0VzPx4M1+YMojh+Rn835vrGF+YRZ0/zNTheby9ehcjCjJZu7OWcf2zWL+7jouPLeS8Cf2546VVRK3FH9pbMe1qbVXSFNJEelA4EuXhuUVcdsJg+mT42n9CzIOzt7BwayUzr2rxdzohvqrS007rtrXndSQ4vLq8lLH9s3h47lbW7arjhRtOZWNZPT9/YSX/+PoJ9MvaGx5bG0tJVSPf+NdCrj9zJFfE2lkLtuxhzuY9fOvU4VQ2BLju0SV8YcogHpi9hfEDsnnmumnNxtYQCHP1w4tYsLWSK08Zyu++cCwA9YEwRRUNTBqUw5srd3LDE0s5dVQ+T147rd3vKxSJ8p1HFrOnIcBz15/Kgq17+MNb67nzcxMprW4ixePijLF9yUzxJN6r4spGNpXX89bKXZw0og9rSmv595ytALxx8xnUNIUYXpDerGJ7oPyhCNc/toSPNpQn/qG7+yuT+cKUwTw6r4j/eXk1AOk+N49efTJTh/dJvPfGGCrqA3zl/nlsLm9gytBc/nLFFFwuw69eWc11Z45sdn4kanloThG/fWMtAF+YMgiPy/CfpSVMGpTDipKaxGu9+L3TGFeYRTgS5epHFvPxhvLEmKePzOexa04mGIny4fpyZozvl2if1vlDVNQHGZ6f3uLn7ZONFRTmpDC6XxbPLi7mnnc3kJ3mZd2uOi44pj+/+a9JLX45eXnZDu54aRV1/jCfmzyQVK+LS44byJlj+7K4qJIFWytpDIb550dbeOI7pyQWBZXV+vn+k0tJ9bp57JpTiEQtf353PffN2szgvDS+dsowdlQ3Mr4wm69PGwbA3E0VXPngAoyBaSPyeeq6vT9XX/rHXDaW1dMQCDM0Px1rnbZbnDFwwtA8PnvsAM4/pj+3PrssMd0B4IShuTxz/XR++MwyXluxs9nzslO9zP3ZDB6fv43/e3NdYosggHu+cjzVjUHufHUN6T43/bNTKa8LUB9bKT48P52qxhB1/hAuY5j7sxkt3sO4mqYQF93zMaU1fsb0y+S1m0/nly+v5ulFzhzZSYOyWbWjljSvm6nD85gyJJdH5m1LLIS65dwx/PD8sVhreXXFTrbvaWDq8D40hSK8s3o3Jw3P49LJA6luDHH+3R9R3RjCGGe7pJ9fNIEfPrPMmSOZn862PY0UZqfywY/OIt3nYcGWPVzzyGKyUj388tKJfPfxJfz+smMp2tPI/R9v5qlrp3HTU58mpg+8cfMZFGT5OPMPs/CHolw1fRg3nD2K8+76iIZghBOH5bFkWxWnjy5gzmbnl7lx/bO454rjueje2S3CbFc5rEKaMeZy4E5gAnCytbbVRGWMyQUeBCYBFrjaWtv2ZJkYhTSR7mOtxVo6NAented2JBBW1AfITPG0OT+qIRAmzetudQyhSJTfvr6WK08Zmmgld2RcUdux1ndrgmGnbdUvK4WHvn1y+0/ooEA4wqfbq5kyNBeDaba/3dxNFawqreHzxw9qs9VW3RikpKqJiQOz233f4wsdlhdX89j8bQB8/5xR/OiCcazcUUNxZRPTR+W3+OWizh9i/a46fB4XdcpJuwAADNVJREFUxw7KOaStovpAmMw22m3gvO+l1U0MayX4xdU0hZpNT4jb92exzh8iM9b+as1f39/IXe9u4PeXHZv4JQPghaUl/OX9jZw2uoCfXTQej8vFvz7ZwvRRBdz5ymoslueuP7VZW25TWT2byuqJWsvpYwrITvVirWV1aS07a/xUNwb58fMr+P45oxKT5t9bs5s5myvIz/DhD0W5+dwxbKmo58J7ZuN1G9685UzK6wLM2VTBl6cOYWh+Oo/NK+IXL6/mwomF/PMbJ+73vW4IhNlYVs+QvDTyM1OoaQrxzupdLNxayXNL/n979x50V1Wfcfz7JAESIqCCRAIUkHKV+wQFakck6Qy3GrBWcGAqHUdmGGupU7RQO9OBGZk67XhFGTMUpdUKFrW14NSCYsFyM9xvQpFLCUQIdxLI/dc/zk54jbm8IXnPWSf5fmbOvGevs9+918lvQh72WnuvOUzfd0cWLV3OI/Pm81R349Nnjt+PZxcs4pMz9l7nXMYVRs69XTH/694nX+KCq+7nwpMP4PRLbuVvTtyPEw+auvJ3bnv8ef7o4pvYasI4tp20BTd8+n0sXLKM3//cdbyyaCkTxoWvnnYYk7YYv3Iu4nPzF/HYc6+y79u3YfJWE7j67rlc/9A8LvzAgcx54VV2234yjz27gGdeWcSBO2/HpC3HM/OinzNh/Di+d9ZRo/ouG6K1kLYfsBz4OnDOWkLaZcANVXVJki2BratqzferdwxpkgZt4ZJljEuG/kHBLy9cwhmX3sqM/adw1nv37Mv8nGFQVTz09Hz2nvKmUf+ZLF66nHFhva9uA9z5xIu8c+q2a72pafny4vgv38AJB+7EJ6b/9tWfxUuXc/5/3Mdp796N/aduu5ojrNvSZcu5+p65HLPvjmwzcQsWLlnGiV/5OS++upibzpu+zpuuNpYVN6z81bH7ctbRewK9cHfN/U+zz9u34fgDd9rgc3z+moe47MbHuPm86au9MWhjaiqkrTx58jPWENKSbAfcCbyj1rOThjRJkvrjmZcX8vLCpSufG9kPjz67gC9d+xAXnHTAykfsbGzzFy1lqwnj+hI81xTSWn5O2h7APOAbSQ4GbgPOrqrRP3VTkiSNqR23nciOb+zC3Bu2xw6T+eKph47pOdY2tN4vYxYPk1yb5N7VvGaO8hATgMOAi6vqUGABcO5azndmktlJZs+bN29Nu0mSJA2FMYuJVTVjAw8xB5hTVbd021eylpBWVbOAWdAb7tzAc0uSJA1UszNaq+rXwBNJ9umapgP3D7BLkiRJfTOQkJbk5CRzgCOBq5P8uGufmuRHI3b9BPDtJHcDhwAX9r+3kiRJ/TeQWXFV9QPgB6tpfwo4fsT2ncCan9YpSZK0iWp2uFOSJGlzZkiTJElqkCFNkiSpQYY0SZKkBhnSJEmSGmRIkyRJapAhTZIkqUGGNEmSpAYZ0iRJkhpkSJMkSWqQIU2SJKlBhjRJkqQGGdIkSZIaZEiTJElqkCFNkiSpQYY0SZKkBhnSJEmSGmRIkyRJapAhTZIkqUGGNEmSpAYZ0iRJkhpkSJMkSWqQIU2SJKlBhjRJkqQGGdIkSZIaZEiTJElqkCFNkiSpQYY0SZKkBhnSJEmSGmRIkyRJapAhTZIkqUGGNEmSpAYZ0iRJkhpkSJMkSWqQIU2SJKlBhjRJkqQGGdIkSZIaZEiTJElqkCFNkiSpQYY0SZKkBhnSJEmSGjSQkJbkj5Pcl2R5kmlr2e+T3X73JvlOkon97KckSdKgDOpK2r3AB4Dr17RDkp2BPwemVdUBwHjg1P50T5IkabAmDOKkVfUAQJJ17ToBmJRkCbA18NQYd02SJKkJzc5Jq6ongX8A/g+YC7xUVf812F5JkiT1x5iFtCTXdnPJVn3NHOXvvwWYCewBTAUmJzl9LfufmWR2ktnz5s3bOF9CkiRpQMZsuLOqZmzgIWYAj1bVPIAk3weOAr61hvPNAmYBTJs2rTbw3JIkSQPV7HAnvWHOI5Jsnd7ktenAAwPukyRJUl8M6hEcJyeZAxwJXJ3kx1371CQ/AqiqW4ArgduBe7q+zhpEfyVJkvotVZveyOC0adNq9uzZg+6GJEnSOiW5rap+67mxLQ93SpIkbbYMaZIkSQ0ypEmSJDXIkCZJktQgQ5okSVKDDGmSJEkNMqRJkiQ1yJAmSZLUoDFbu3OQHpm3gFO+ftOguyFJkvSGeSVNkiSpQS4LJUmSNEAuCyVJkjREDGmSJEkNMqRJkiQ1yJAmSZLUIEOaJElSgwxpkiRJDTKkSZIkNciQJkmS1CBDmiRJUoMMaZIkSQ0ypEmSJDVok1y7M8k84PExPMUOwLNjeHyNPWs43Kzf8LOGw88abjy7VdXbVm3cJEPaWEsye3ULoWp4WMPhZv2GnzUcftZw7DncKUmS1CBDmiRJUoMMaW/MrEF3QBvMGg436zf8rOHws4ZjzDlpkiRJDfJKmiRJUoMMaeshybFJHkzycJJzB90frV6SS5M8k+TeEW1vTXJNkv/tfr6la0+SL3c1vTvJYYPruVZIsmuS65Lcn+S+JGd37dZxCCSZmOTWJHd19Tu/a98jyS1dna5IsmXXvlW3/XD3+e6D7L9el2R8kjuSXNVtW8M+MqSNUpLxwFeB44D9gQ8n2X+wvdIafBM4dpW2c4GfVNVewE+6bejVc6/udSZwcZ/6qLVbCvxlVe0PHAF8vPv7Zh2HwyLgmKo6GDgEODbJEcDngC9U1e8CLwAf7fb/KPBC1/6Fbj+14WzggRHb1rCPDGmj9y7g4ap6pKoWA5cDMwfcJ61GVV0PPL9K80zgsu79ZcBJI9r/qXpuBt6cZKf+9FRrUlVzq+r27v0r9P6R2BnrOBS6OszvNrfoXgUcA1zZta9avxV1vRKYniR96q7WIMkuwAnAJd12sIZ9ZUgbvZ2BJ0Zsz+naNBymVNXc7v2vgSnde+vauG7Y5FDgFqzj0OiGye4EngGuAX4FvFhVS7tdRtZoZf26z18Ctu9vj7UaXwQ+DSzvtrfHGvaVIU2bnerd0uxtzUMgyZuA7wF/UVUvj/zMOratqpZV1SHALvRGIvYdcJe0HpKcCDxTVbcNui+bM0Pa6D0J7Dpie5euTcPh6RXDX93PZ7p269qoJFvQC2jfrqrvd83WcchU1YvAdcCR9IahJ3QfjazRyvp1n28HPNfnruo3/R7w/iSP0ZvecwzwJaxhXxnSRu8XwF7dnS1bAqcCPxxwnzR6PwQ+0r3/CPDvI9r/pLs78AjgpRHDaRqQbi7LPwIPVNXnR3xkHYdAkrcleXP3fhLwB/TmFV4HfLDbbdX6rajrB4Gflg/xHKiqOq+qdqmq3en9e/fTqjoNa9hXPsx2PSQ5nt4Y/Xjg0qr67IC7pNVI8h3gaGAH4Gngb4F/A74L/A7wOPChqnq+CwMX0bsb9FXgT6tq9iD6rdcleQ9wA3APr8+H+Wt689KsY+OSHERvEvl4ehcDvltVFyR5B72rMm8F7gBOr6pFSSYC/0xv7uHzwKlV9chgeq9VJTkaOKeqTrSG/WVIkyRJapDDnZIkSQ0ypEmSJDXIkCZJktQgQ5okSVKDDGmSJEkNMqRJGnpJliW5c8Tr3HX/1mqP87Mk0zZ2/0Zx3pO6BeQlaaUJ695Fkpr3WrcE0bA6CbgKuH/QHZHUDq+kSdokJTk2yb+O2D46yVXd+4uTzE5yX5LzR3Gsw5PcmOSuJLcm2SbJxCTfSHJPkjuSvK/b94wkF4343au6h4GSZH6Sz3bHuTnJlCRHAe8H/r67CrjnRv6jkDSkDGmSNgWTVhnuPAW4Fnh3ksndPqfQe1I6wGeqahpwEPDe7gn5q9UtA3cFcHZVHQzMAF4DPk5vnfcDgQ8Dl3VPXV+bycDN3XGuBz5WVTfSW1LnU1V1SFX96g18f0mbIEOapE3Ba13AWfG6oqqWAv8J/GG34PMJvL7O4IeS3E5vWZt3AmubD7YPMLeqfgFQVS93x34P8K2u7Zf0lqnaex39XExvWBPgNmD39fyekjYjzkmTtCm7HPgzemsJzq6qV5LsAZwDHF5VLyT5JrCuK2DrYym/+T/AI4+9ZMSi08vwv8GS1sIraZI2Zf8NHAZ8jNeHOrcFFgAvJZkCHLeOYzwI7JTkcIBuPtoEegvAn9a17U1v0fcHgceAQ5KMS7Ir8K5R9PMVYJv1+F6SNgOGNEmbglXnpP0dQFUtoze8eFz3k6q6i94w5y+BfwH+Z20HrqrF9OazfSXJXcA19K6OfQ0Yl+QeenPWzqiqRd3xHqV3p+aXgdtH0f/LgU91NyB444AkAPL6lXdJkiS1witpkiRJDTKkSZIkNciQJkmS1CBDmiRJUoMMaZIkSQ0ypEmSJDXIkCZJktQgQ5okSVKD/h9d0IIkvrTs1QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=[10,5])\n",
    "plt.plot(counts, values)\n",
    "plt.xlabel('Eval count')\n",
    "plt.ylabel('Energy')\n",
    "plt.title('Convergence with noise');\n",
    "plt.axhline(ref)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "NoiseModel:\n",
      "  Basis gates: ['cx', 'id', 'u2', 'u3']\n",
      "  Instructions with noise: ['measure', 'u3', 'cx', 'u2', 'id']\n",
      "  Qubits with noise: [0, 1, 2, 3, 4]\n",
      "  Specific qubit errors: [('id', [0]), ('id', [1]), ('id', [2]), ('id', [3]), ('id', [4]), ('u2', [0]), ('u2', [1]), ('u2', [2]), ('u2', [3]), ('u2', [4]), ('u3', [0]), ('u3', [1]), ('u3', [2]), ('u3', [3]), ('u3', [4]), ('cx', [0, 1]), ('cx', [1, 0]), ('cx', [1, 2]), ('cx', [1, 3]), ('cx', [2, 1]), ('cx', [3, 1]), ('cx', [3, 4]), ('cx', [4, 3]), ('measure', [0]), ('measure', [1]), ('measure', [2]), ('measure', [3]), ('measure', [4])]\n",
      "VQE on Aer qasm simulator (with noise): (-1.7623537800696387+0j)\n",
      "Delta from reference: [0.09492125+3.77606111e-17j]\n"
     ]
    }
   ],
   "source": [
    "device = provider.get_backend('ibmq_valencia')\n",
    "coupling_map = device.configuration().coupling_map\n",
    "noise_model = NoiseModel.from_backend(device.properties())\n",
    "basis_gates = noise_model.basis_gates\n",
    "\n",
    "print(noise_model)\n",
    "\n",
    "backend = Aer.get_backend('qasm_simulator')\n",
    "quantum_instance = QuantumInstance(backend=backend, seed_simulator=167, seed_transpiler=167,\n",
    "                                   noise_model=noise_model)\n",
    "\n",
    "counts1 = []\n",
    "values1 = []\n",
    "def store_intermediate_result1(eval_count, parameters, mean, std):\n",
    "    counts1.append(eval_count)\n",
    "    values1.append(mean)\n",
    "    \n",
    "\n",
    "aqua_globals.random_seed = 167\n",
    "optimizer = SPSA(maxiter=200) \n",
    "from qiskit.circuit.library import RealAmplitudes\n",
    "var_form = RealAmplitudes(num_qubits)\n",
    "vqe = VQE(qubit_op, var_form, optimizer, callback=store_intermediate_result1)\n",
    "vqe_result1 = vqe.run(quantum_instance)\n",
    "print('VQE on Aer qasm simulator (with noise): {}'.format(vqe_result1['eigenvalue'])) #'energy' deprecated\n",
    "print('Delta from reference: {}'.format(vqe_result1['eigenvalue']-ref)) #'energy' deprecated"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.lines.Line2D at 0x7ffa6055c580>"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmkAAAFNCAYAAABbpPhvAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+j8jraAAAgAElEQVR4nOzdd3hc5Zn+8e87Vd29VzA2xnQwPZQktISwpC6kQpJfCJtCCunJZglLstk0SBZSHGoIJUBwQgvFpmNj446xcW+ybEtWb6Np7++Pc85oJI2ksWSN5PH9uS5fSDNH57ySDb55nrcYay0iIiIiMrT4BnsAIiIiItKVQpqIiIjIEKSQJiIiIjIEKaSJiIiIDEEKaSIiIiJDkEKaiIiIyBCkkCYiMkiMMVONMU3GGH8P11hjzFE5GMsfjTH/OdDPEZHsKaSJSAfGmE8YY5a54WGPMeZfxph3Dfa48pG1dqe1tsRamwAwxrxkjPl/gzSW66y1/z0YzxaRzBTSRCTFGPNN4FbgZ8A4YCrwe+CKwRxXOmNMYLDHICKSCwppIgKAMWYYcBPwZWvtY9baZmttzFr7hLX22+41YWPMrcaYCvfXrcaYsPveBcaYcmPMDcaYSrcK91n3vTOMMXvT23rGmA8ZY9a4H/uMMd8zxmwxxlQbYx42xox035vutvw+b4zZCbxgjPEbY35tjNlvjNlmjPmKe03A+16MMXe6Y9htjLnZe7Yx5hpjzGvGmF8ZY2rdr39f2rhGGmPudr+/WmPMP9Le+4AxZpUxps4Ys8gYc0I3P8ufGGP+z/04aIxpNsb80v280BgTcZ/jfW8BY8xPgXOB29wq5m1pt7zQGLPJfe7txhjTzXNvdH92fzHGNBpj3jbGzE17/xi3Wlfnvvdvae/dY4y52f14tDHmSfe6GmPMq8YYn/veRGPM340xVe7P7vqe/lyJSN8ppImI5yygAJjfwzU/BM4ETgJOBE4HfpT2/nhgGDAJ+DxwuzFmhLV2CdAMvCft2k8AD7gffxX4IHA+MBGoBW7v9OzzgWOAS4AvAO9zx3GK+7Xp7gHiwFHAycDFQHob8QxgAzAa+AVwZ1rwuQ8oAo4FxgK3ABhjTgbuAr4IjAL+BDzuhdROXgYucD8+DdgLnOd+fhawwVpbk/4F1tofAq8CX3FboF9Je/sD7n1OAP7d/Rl059+Ah4DhwOPAbe74g8ATwHPu9/VV4H5jzNEZ7nEDUA6Mwamo/gCwblB7AliN83v8XuDrxpiexiMifaSQJiKeUcB+a228h2s+Cdxkra201lYBPwE+nfZ+zH0/Zq19GmgCvBDwIPBxAGNMKfB+9zWA64AfWmvLrbVtwI3ARzu1Nm90q3utOEHlt+71tcDPvYuMMePce3/dvb4SJ2hdlXavHdbaP7tzwe4FJgDjjDETcMLfddbaWvf7eNn9mmuBP1lrl1hrE9bae4E2nNDa2WJgpjFmFE44uxOYZIwpwQmbL2f4mp783FpbZ63dCbyIE06785q19mn3e7sPJ0zjjrPEvVfUWvsC8CTu70knMZyfyTT3Z/CqdQ56Pg0YY629yb3HVuDPdPzZishBopAmIp5qYHQvc74mAjvSPt/hvpa6R6eQ14ITDMCpmn3YrTx9GFhhrfXuNQ2Y77bX6oD1QAKniuPZ1Wkcu7p5bxoQBPak3e9PONUjz17vA2tti/thCTAFqHGDX2fTgBu8e7r3ndLp+/fu2Qoswwlk5+GEskXAOfQtpO1N+zj9Z5rNtQXu7+lEYJe1Npn2/g6cilhnvwQ2A88ZY7YaY77nvj4NmNjpZ/ADOv4+ichBogm4IuJZjFMZ+iDwaDfXVOD8Rf22+/lU97VeWWvXGWN24FSq0lud4ISsz1lrX+/8dcaY6d4t0l7eA0xO+3xKp3u1AaN7qQpmsgsYaYwZbq2ty/DeT621P83yXi/jtHdPBt50P78Ep0X8SjdfY7t5/WCoAKYYY3xpQW0qsLHLIKxtxGl53mCMOQ5nHuCbOD+DbdbamQM4ThFxqZImIgBYa+uBH+PMI/ugMabInfT+PmPML9zLHgR+ZIwZY4wZ7V7/1wN4zAPA13CqS4+kvf5H4KfGmGkA7v17WlH6MPA1Y8wkY8xw4Ltp38cenHlXvzbGlBlnUcIMY8z5vQ3O/dp/Ab83xoxwv39vLtmfgeuMswjCGGOKjTGXua3bTF4GPgOss9ZGgZdw5sVtc1vFmewDjuxtnH20BKey9h33+7oAuBxn/loH7gKJo9x5evU4Vc0ksBRoNMZ8110A4TfGHGeMOW2AxixyWFNIE5EUa+2vgW/iLAaowqmcfAXwVjjejNPGWwO8BaxwX8vWgzjtvhestfvTXv8tziT354wxjcAbOJP7u/NnnCC2BlgJPI2zUCDhvv8ZIASsw1mE8CjOHKtsfBpnTtY7QCXwdQBr7TKcBQu3uffcDFzTw30WAYW0V83WARG6r6KB83P4qLuq9HdZjjcrblC8HKeSuR9na5XPWGvfyXD5TGABzpzCxcDvrbUvuvPcPoAzJ26be587cBaLiMhBZpy5oCIihy7jbKHxR2vttMEei4jIwaJKmogcctxW2/vd/cUmAf9Fz1uHiIgcclRJE5FDjjGmCGfO12ygFXgK+Jq1tmFQByYichAppImIiIgMQWp3ioiIiAxBCmkiIiIiQ1BebmY7evRoO3369MEehoiIiEivli9fvt9aO6bz63kZ0qZPn86yZcsGexgiIiIivXJPY+lC7U4RERGRIUghTURERGQIUkgTERERGYIU0kRERESGIIU0ERERkSFIIU1ERERkCFJIExERERmCFNJEREREhiCFNBEREZEhSCGtH6qb2li7u36whyEiIiJ5SCGtH+54bRtX37V0sIchIiIieUghrR+aInGa2uI5e159a4wHluzEWpuzZ4qIiMjgUEjrh1giSSyRzNnznl+3jx/Mf4vy2tacPVNEREQGh0JaP0QTSZIWEsncVLaicScQtsVzFwxFRERkcCik9UMsYd1/5iY0JZLJnD5PREREBo9CWj/E4l1DUyJpB2zOmBcKo6qkiYiI5D2FtH6IpypbTnhqbosz4wdP8/uXtmS8fn9TW8bXH1y6k4tvebnX53ltVVXSRERE8p9CWj9EO7U799RHAHh0eXmXa1fvquO0ny5g2/7mLu99/7G32LivqdfnxZOqpImIiBwuFNL6oXO7syESA6CsINDl2j31EaztvpoGvS9A8OakRVVJExERyXtd04RkzQtnXruztjkKQGlBsMu1XrDqqQoWSyTx+/xdXq+oa6WmOZqqpHnPExERkfylkNYPsU5zxKqbvJDW9cfqhbOeqmDdVdJ+u2ATb26v4f3HT+hwLxEREclfanf2g9fu9ELT/manlZkppLXFEx2+JpN4NxWypmiclmgirZKmkCYiIpLvVEnrgxsff5vRJaFUWPLCk1dJC/jbs28yaUlYmwpyPbUqY8nM4SsaTxJPJtvnpMWTtMUThPw+jDH9/4ZERERkyFElrQ+W7ahh5c66tDlpXrvTqaSltyPvXrSdi37zclpIO/B2p3P8lE2FwaqmNk6+6Xle2bS/z99Drk5JEBERkb5RSOuDcMBPWzzZfuKAG8Cq3YUDsUSSbz68ip8+tY6d1c3sqGlJHeXUeU5aeljqLsDFEkniiWTq2j31rbREE+yo7rqdRzaWba9hzo+f6XGlqYiIiAwuhbQ+CAd8tMUTqVDlBa/9Te0hbV1FA+v2NNAWT2Kts9Gt9166pkg89XF3c9JicaeK5oXCRvdrWqKJPo1/pxsa97r7uomIiMjQo5DWB+GAj2g82T4nzQ1PNc1eu9OZgxaJJYnEnCDV4AarzgsHvL3VoH1um3OPJEn387ZEknjSpuak9TekeeP2FjOIiIjI0KOQ1gehgK9juzORxFqbWjgQSyRpizsBLRLruNFt54UD9a3pIc25NpG0nPO/L/Dwsl3O18SdVqcXBr3qW2s0Tl94JyV4YxMREZGhR6s7+8CbkxZNa3c2tMY7HNvkhTSvWuVVvzrPSetQSfNCWFucqsY2trpHSLVXvjoGvgOtpEXdMXvVPK/KJyIiIkPPoFTSjDEfM8a8bYxJGmPm9nDdpcaYDcaYzcaY7+VyjD0JB3y0xRIdThzw9khzPne2yHDanc41XsWs80a0Da3xDl8HTkgDaHTDmBfsvFDVmKqkHVjI+t3CTXz0D4vS7qdKmoiIyFA1WO3OtcCHgVe6u8AY4wduB94HzAE+boyZk5vh9Swc9NESS2DdzmU8kUwtDAAnbHl7mUW8SlprLPVeuvRKmrd602tndp7H5t2rt0pabXOUz9/zJpWNHRcG7KptoaKuVZU0ERGRQ8CghDRr7Xpr7YZeLjsd2Gyt3WqtjQIPAVcM/Oh6F/L7O6zKjCWSqapWUchphba5CwfaUnPSMq/ubEibkxZLa3dCeovUed17hvd+Szcha+WuWha+U8mKHbUdXm9uixNL2NQYIvEEDZEY1mrPNBERkaFmKC8cmATsSvu83H1t0IWDvo4rMROWiFudKi0IpCpckViiS/Wr88KBhvQtOJLdtDvde3hz0rxM5S0cWLu7nk37GlP3qWxwWq/evm2eprY40UQyFfpqmqKc+bOFPLlmzwH+BERERGSgDVhIM8YsMMaszfBrQKphxphrjTHLjDHLqqqqBuIRKeFAxx9bPNG+1UZZQTAVsuJJS0ub83r6AevzV5an9ihraO26BYdXpWtMVd+81ZgdK2ct0QT3L9nBB/7vNb7wl2Wp16sanZBW09QxpDW3JUgkbeo+FfURWqIJ1u9pSF2zubKJZ9YqtImIiAy2AQtp1toLrbXHZfj1zyxvsRuYkvb5ZPe17p43z1o711o7d8yYMf0Zeq9CnUJaLD2kFQY7tELTt9gAJ5R942+rU9trZFrd2dypkhbrZqJ/azTBbS9sBqAu7TmVjZkrad59vX/Wuu+nb2p752vbuP6hVak92kRERGRwDOV255vATGPMEcaYEHAV8PggjwlwtuBIF020V6dKCwIdttlo7VT98tqbtS1OQKpriRHwOYekx92va0ybk5ZMtp/Z2Xnz2dZYInWflmgiNbfMWzBQk6HdCdDstkm99yvqW6lqbKMhEqOuJUo0nqSivjXbH4eIiIgMgMHaguNDxphy4CzgKWPMs+7rE40xTwNYa+PAV4BngfXAw9batwdjvJ11bnc6lTQnYJUVBHv8Wq+9Wdfi/LOmOcqY0rBzn07tzpZookPI61xJq22JEoklGVUcIhpPprbz8CppnUNac6cFCdXutiF76iNcfddSfvL4ulTo276/pcfvQ0RERAbWoGxma62dD8zP8HoF8P60z58Gns7h0LLS05y00oKef6ReC7PODUO1LVHGlobZUx9JHfvUnHaSQHrQ6jwnzQttM8eVUL21hqqmCMOKghkXDiSTlmZ3QYMX1lKVtLpWEklLSTiQar9uq27mXTNH9/yDEBERkQEzlNudQ1bXOWk2VfEqK+y5kuZVsbw5ZB0qaZ0OUPfe98S7mSc2c2wp4FTQrLXtCwfSNthN366j2V3MUNfavuI0aZ3KmjeHbltVc4/fh4iIiAwsHQvVB13npDntzpDfR0Gn9zpriLS3O2OJJI2ROGNKCwBn4cDa3fUdNsatTgta3Zk5rgRwVnU2tDrbbIQDPmqao1hrMcZ0uKc3N63z9mg1zdFUdW57tUKaiIjIYFJI64NwsFMlzT2nMxz0EQyYHr/WC0F1LdHUvLSxbiXt8dW7+cH8tzpcX91pG41MvEpaVWNbatHArHGlvLW7nsa2eIdtQaBjOzVdbUv7CtHt+xXSREREBpPanX3QZU6au/dYYdBPyJ/dj7S+NZaqko0tc0JaRV37VhjFIaci13nyfyZTRxURCvjckObcc/Z4J7htq2omkbQdKmnpH2cyrDDIzpqW1GpTERERyT2FtD7oHNKi7sKBgqCfYJYhLWlhR7WzgnJMiRPSygrbC5sThxcCXfc6y2REUZAxJWH2NkR4dHk5ACdMHgbAFbe/zifveCM1Tw06nnoQDvgIB3yMdsfgfW08adlT3/HsTxEREckdtTv7oPOctFg8ScIYCoK+LosKerLNbSmOLXPmpLWlbbExYXghmyqbum13+n2GRNISCvgoDPoZUxrmn6sqALj+vTM5ccrw1LVLt9WwqybzvmelBQFu/8QpVDa28dUHVwJw/KRhvLppP7tqWpgysijr70dEREQOHoW0PkivpIUDPmKJJBYoPIBKGrSvoPTmpHlnfgJMcINbTTcLB0YUBdnfFGV4YRBjTKpS9okzpvLNi2alThE4acpwCoN+Fm+tznifwpCfM44c1eHsz+MnOVW4XbXaK01ERGSwqN3ZB+mVtKKQn3jS0hpNEA76Cfp7XjiQbuv+JgBGlYQwhg6T+0cUhwgHfN22O0cUhQAYXuRs+XH12dM4elwp/3nZHADGDyvgvs+fzoNfOJNpo7qvhhWHnJw+sjiUem3OxDL8PtNt9U1EREQGnippfZDe0iwKBYjGk0TiSYYVBjMuHCgI+ojEkhjTcduLbfubKQkHCAf8BHwdt8koLQhQWhDstt3ZHtKcf1573gyuPW9Gh2vOnemcYTptVHG330uRu0BheFEoNb7RJWEmDCtQJU1ERGQQqZLWB+EOIc1PLJGkLZagMOjL2O4cXugEqZJQx0y8vynKiGKnEhbw+TpsVlsc8lNWEOh2nzTv64b3snku0HMlLeyMye8zjCwKEfQbikJ+powoYleNE9KefXsv9S2xbu8hIiIiB59CWh+k75OWane6qzvTq2zex8PcIFXiHhmVHvJGupWwQKc26fCiEKUFgS7ndXo6tzt7MrWHyf+FwfbW7cjiEMMKQxhjmDKykF21rVTUtfLF+5Zz/UMre32OiIiIHDwKaX2Q3tJMtTtjCQoC7QsHQn5fKgB5Ic2rWpWE2ytqx0woAyDgc0JaQdDHHz91KpceN57SDIe1e/f02pxeWOuJV0nzWpvpitPGMqokxAg39E0ZUdRh37WKOs1PExERySXNSeuDgN+X2gKjKOSnstE5Fqow5CfknjgQCvgoCPqob20/z9MLRAVBP99/32yMgavPnp66JziLEi49bjyQ+bD2opCf1liC0oIAX3n3UVw0Z1yv4y0tCDKqOERB0E9L1AlbAZ8h7o7f8+kzp6eOrfK23vBOHjiQrUVERESk/xTS+igc8NESTVAY8qcOWA+nzUkLB3wUdKqklYSdzwuCPr54fsdJ/kG3kpbeCs0U0rx7hgM+vvzuo7Ie77RRRcQSlt1uRay0IEBtS6xDJe2yEyakPvY21vVWlyqkiYiI5Jb+5u2jcMCHMU5oisaTRONJCgLtc9JCAeewdb/PpMJZSVolrbNUJS2YHtK6tju9ylc4wz168t1LZ/PdS2entgjx7p2pBZo+xvpWp7KW7XFXIiIicnCoktZHoYBTNQsFfKn9zQpD/k6VNB8FaRW14p5CmltJSw9D3bU7AQoOsLJ1xpGjUvePJRKpe3cX0ry5bw1eSFMlTUREJKf0N28fhQPOYeohf3tIKwj4UiErHPATDvopCPpTLUyvklaYsZJmUl/nyVRJS7U7D7CS5vHCVlmqkpY5pxeGOoU0VdJERERySn/z9lE44CPgN6kKGNDhgPWQW0ELB3ypQNVeSev6Yw/4MrU72wOU331OXytpHi+kefcuDvdcSatXJU1ERGRQ6G/ePvIWCQTTwovT7mxfAFCQFtSg5zlp6V/nKXODVNBvUu8X9nFOWvtzvJDWSyXNvX+dQpqIiMig0Jy0PvJaneknDIQ7LRz41JnTqGxsY+KwAo4cU5w6SD1TSPMqZaEM7c6g34ffGCIkOWpsKROG1XHk6O6Peupx3O74Tpg8jEgswSlTR2S8rsANg7Xe6k61O0VERHJKIa2PwgGnahbs0O7suAXHebPGpN574YYLeGbt3tR1nQXSvs5TlhbSvMdMH1XE4u+/t8/j9sLWiOIQt3/ylG6vUyVNRERkcCmk9VE46CPQud3ZaU5aZ95Gt5kWDmRqd3rzxkIBX+pgdr/PdPnaAxq3V+nz93yfoN9HwGeoa9E+aSIiIoNBIa2PjhxdQtDv6xCqCkPOvmh+n+mwStPjBbjM7c72VaGeVEjz+0i6Kc1bYNBXXtjKdBB8Z4VBP43uytVsrhcREZGDRyGtj358+RwAKhsj1DRHCfh8zHHP4Qx1Cm+eUA8hLehrP07K481JCwV8xJPOQev9raSlz5nrTUGoPaQlk7ZfzxUREZEDo5DWT2NLC7jh4qM7vBb0m4whyGuNZj5xoGu7MxRwwp7TCnVeD/Q3pPkPrJLmUUYTERHJLYW0AfDDy45hzoRhXV5vr6T1sHCg03ulBUGCfh/WJt3rDk4l7cBDmlKaiIhILimkDYArT5ua8fXUnLQM89W8Clm4U3gqKwgQCvhSlaz+z0lznp3NlhoFaUdGWYU0ERGRnFJIy6GpI4t4z+yxzJ3edW+y9hMHOga40oIAQb+PhJvS+j0nzWt3Bnq/T2FaVS+hkCYiIpJTCmk5VBjyc9c1p2V8L9MWHACfPms6fh/8ZfEO4GC0O437vCwqaZqTJiIiMmi0r8IQkWnhAMBHT53Mh06eTNCttB2sSlo27c70OWlqd4qIiOSWQtoQ4bU7u9sawwtxwYO0T1o2W3Ckh7SESmkiIiI5pZA2RKQWDmRYVADtqz8P1j5pWbU7Q/1rdz6/bh+n/XQBjZHYgX+xiIjIYU4hbYjIdHZnh/fdcNbvOWl+J3gFs7hPf7fgePbtvVQ1tvFWeT3zV5bTGk0c8D1EREQOV1o4MESkKmkZ9lBLf7+/lbTpo4sYWxrOeH5oZx3npGX/jETSEkskWba9BoA7X9vGwncqSX4MPnLq5AMes4iIyOFIIW2I8CpkXqWrM6892d8TB644aRJXnDQpq2sLQwc+J62yIcLVd79JU1uMXTWtACx8pxKAPfWtBzhaERGRw5dC2hAR7ObEAY8/1e7MXYe64ADbnTuqm/n0nUvZ1xChLe6ckDC+rIC9DREA9tRHBmagIiIieUhz0oaI9oUDPa/u7G8l7UAcaLvzi/ctpyES46Frz+R9x42nNBzgo2ntzX0N7SHte39fwzV3Lz2o4xUREcknqqQNEf5eVncerH3SDkRhqD0wZlNJq6hr5UMnT+LkqSO45cqT2NcQYXdtK7e9uJlxZeFUJS2WSPLQm7sGbNwiIiL5QJW0IcJrd/a2T9pgVdKymZNmLfjc8RUE/UwbVczZR43m1e+8m/fMHsteN6S9vnn/wAxYREQkjyikDREF7ly07lZdHqzVnQc2pgPbJy1pLT7TdXxTRhYxvqyQ6uYobfEET63Zk3ovEtO2HCIiIpmo3TlEXHrcBEIBH+OHFWR831swkM0mtAfLgR4LlbTQXYac4H5flQ1tHSppDa2xDmFQREREHKqkDRHDCoN86OTu9xDz2p25rKSVFLRn+EQWIS1hbard2ZkXPlftqqOiPsLJU4cDUN+q0whEREQyUUg7RKTanRnaiQNlzoQyfnvVSRw7sSyrdqftpt0J7ZW0Z9buBeCCWWMBaNCRUSIiIhkppB0iLjh6LJ86c2q3laqBYIzhipMmEQ74+t3u9CppL7xTScjv46wZowBV0kRERLozKCHNGPMxY8zbxpikMWZuN9dMMca8aIxZ5177tVyPcyg5bfpIbv7g8YPybJ8xWW3BkUjabit9pQVBPnzKJFpjCeZMLGN0SQjIPqS1xRPc9sIm2uJaaCAiIoeHwaqkrQU+DLzSwzVx4AZr7RzgTODLxpg5uRicdOQzptctOLxKm+mhHXvTFcdx7MQyLjl2PMMKgwA0tMazGsPy7bX86rmNLN9Rm+WoRUREDm2DsrrTWrseev4L3Vq7B9jjftxojFkPTALW5WKM0s7n630LDu/97uakAZSEAzx1/bmAs6EtZF9Ji7gVtGzPEBURETnUHRJz0owx04GTgSU9XHOtMWaZMWZZVVVVroZ2WPAZ0+ucNK8dmu2UuaDfR1HIn327M+aEurhCmoiIHCYGrJJmjFkAjM/w1g+ttf88gPuUAH8Hvm6tbejuOmvtPGAewNy5c/U3+UHkzEnr+RqvwnUgCxuGFQZpyHpOmhPSkgppIiJymBiwkGatvbC/9zDGBHEC2v3W2sf6PyrpC2N6bzPaLNqdnQ0rDB7QwgFQu1NERA4fQ7bdaZwJa3cC6621vxns8RzO/L6D3+4EKCvIPqRF3HanQpqIiBwuBmsLjg8ZY8qBs4CnjDHPuq9PNMY87V52DvBp4D3GmFXur/cPxngPd1m1O92QdiAnIpQVBmmIZLe6M1VJy2IrEBERkXwwWKs75wPzM7xeAbzf/fg1IHc7t0q3fNm0O51CV48rdjsbVhhk/Z5upxkCsLc+wuryutTCAVXSRETkcDFk250ydJgsNrPtU7uzMNBru/Mvi7fzH39dTktMc9JEROTwopAmvfIbQ29dxr60OyePKKKpLc6+hki311Q1tpG07fupKaSJiMjhQiFNeuVsZptdJe1A2p2nTB0OwIq0UwT21LemNroFqG6OAlDXEu3wHBERkXynkCa9Msb0OmG/fQuO7O977MRhhAK+1FFPrdEEF/76ZW5/cXPqmuqmNgDqWrxK2gEMXERE5BCmkCa9yqbd2T4nLfuUFgr4OHHyMFbsdELaO3sbaI4meHLNntQ1+5ucClptKqQppYmIyOFBIU165TO9txm9uWL+AwhpAKdMHcHa3Q20RhOsrXBWem6ubGJLVRPWWqqbvUpatMNzRERE8p1CmvTKl8XqTu/tA8xovGf2WKKJJA8s3cm6inoKg34A/vbmLpqjidQmtrVeSFNGExGRw8Sg7JMmhxZjDL11GfvS7gQ448hRnHPUKH7/4mbKCoOcMm04w4tCzHtlK5Vpqz7bTxxQu1NERA4PqqRJr/xZrO5MtTsPZOWA67uXzqa6Ocq2/c0cPa6M3111MufNGsM/VlVkeM4B315EROSQpJAmvcqm3ZnsY7sT4ITJw/ntVScBcPaMUfh9hnOPGt3Nc9TvFBGRw4PandIrk8XZnbaP7U7PFSdN4qI54ygKOX8kT5g8LON1cU1KExGRw4Qqadd/MbYAACAASURBVNIrn4FkLynNe7sv7U6PF9AAjps0LOOeazpgXUREDhcKadIrv6/3dqc3J60fGa2D4nCAmWNLKQr5O7zeW1gUERHJFwpp0itfFu3OvhwL1Zvzjx7DnAllFATb/5jGFdJEROQwoTlp0iuTRbuz/ViogxfSvnfpbADm/nQBkZjO7hQRkcOLQpr0yp/F6k5vrpj/INZmfW7v1NvgFnTigIiIHD7U7pRe+XyD0+70FIYU0kRE5PCjkCa9Mlmc3dnfLTh6okqaiIgcjhTSpFcHspntgR6wno0OlTTNSRMRkcOEQpr0yp/F6s6DvQVHuvRKmrbgEBGRw4VCmvTK10O7M5m03PP6Nprb4sDAzElL3ytNW3CIiMjhQqs7pVfGGKx15p11DmHPvL2XG59Yx+zxpUD/ThzojippIiJyOFIlTXrlLQbIlI+qm539y5qjcffag/98zUkTEZHDkUKa9Mrb+yxTyzMSTQAQDjhBakC24Aiq3SkiIocfhTTplUlV0roGpNaYF9KcP0oDUUlLn5OmdqeIiBwuFNKkV6l2Z7Lrey3RjiFtIOakFWgzWxEROQwppEmvemx3upW0UKqSNgCrO912ZzjgU0gTEZHDhkKa9MrXU7vTraQF3SQ3ABmNYUVB55+FQRLWcsPDq3l98/6D/yAREZEhJKuQZoz5tTHm2IEejAxNpofVnd6ctHjCO2D94Ke09x03gXs/dzoThxcSSyT5+4pyFm+pPujPERERGUqyraStB+YZY5YYY64zxgwbyEHJ0OLlLm/S/qubqvifp9cDaSHNnbA2EO3OgqCf82eNIeAztMWc58QyTZATERHJI1mFNGvtHdbac4DPANOBNcaYB4wx7x7IwcnQ4FXHvHbn8+v2cd8bO4D2OWnRxMAdC+Xx+QyRuPO8REJz00REJL9lPSfNGOMHZru/9gOrgW8aYx4aoLHJENG53RmJJYglnEqWNyctFh+4SprHbwwRt5Km/dJERCTfZXUslDHmFuByYCHwM2vtUvet/zXGbBiowcnQ4FXHrFtJi8SSxBKWZNKmtuAYyHanJ+A3qcpdXO1OERHJc9me3bkG+JG1tjnDe6cfxPHIEOQFr0QqpLnVs2Sy/eNUu3PgQpovrZKmrThERCTfZRvSVgNHdzrypx7YYa2tP+ijkiHF37nd6bY2YwmbWjgQdV8bwIyG32do6xQKRURE8lW2Ie33wCk4FTUDHAe8DQwzxvyHtfa5ARqfDAGm0+rOVPUsnqS5zTlY3Ws/DsQWHB5/+sIBVdJERCTPZbtwoAI42Vo711p7KnAysBW4CPjFQA1OhobOm9m2pVZ0JmmOdtwnbaAXDngVNC0cEBGRfJdtSJtlrX3b+8Rauw6Yba3dOjDDkqHElzoWyvmnNy+sMRJLVbSiCW/hwMCNI71KF09o4YCIiOS3bNud64wxfwC87TaudF8LA7EBGZkMGV0qaW7Lsaa5/bfe25LDN4ApLf3eqqSJiEi+y7aSdjWwGfi6+2srcA1OQNOGtnnOC2npW3AA1LZEU9fkot0ZSAtpmpMmIiL5rtdKmruJ7dPW2ncDv85wSdNBH5UMKaktONwOozd5v7Y5LaQlc3DiQFoAjKndKSIiea7XSpq1NgEkdV7n4cufmpPWcXVnTVolzWMGcuFA2p9WVdJERCTfZTsnrQl4yxjzPJDa0NZae/2AjEqGFJM2J81am2p31rV0nY44sFtwtKe0uPZJExGRPJdtSHvM/XVQGGM+BtwIHAOcbq1d1sO1fmAZsNta+4GDNQbJXvucNGiLt7cZG1q7hrSBXd3Z/rGOhRIRkXyXVUiz1t5rjCkEplprD8ZZnWuBDwN/yuLarwHrgbKD8FzpAy94JZKWtlh7OGp0N7LteO3A7pPmUbtTRETyXVarO40xlwOrgGfcz08yxjze14daa9dnE/aMMZOBy4A7+vos6T9v64uktalFAwCNESekhQPtf4wGNKSltTt1LJSIiOS7bLfguBHnIPU6AGvtKuDIARpTuluB7wDqbQ0iX9rZnd6iAWhvdxaHA2nXDtw4tHBAREQOJ9mGtFiGg9R7DE7GmAXGmLUZfl2RzQONMR8AKq21y7O8/lpjzDJjzLKqqqpsvkSy5AWv9EUDAE1uu7Mw6E+7Nleb2Sq3i4hIfst24cDbxphPAH5jzEzgemBRT19grb2wn2M7B/g3Y8z7gQKgzBjzV2vtp7p53jxgHsDcuXNVZjmI2vdJsx0qaY0Rp5JWFEoLaQNYSgvoxAERETmMZFtJ+ypwLNAGPAg04Jw8MGCstd+31k621k4HrgJe6C6gycDqrt3ZGIkTDvhS224MZKsTOi4c0BYcIiKS77IKadbaFmvtD621p1lr57ofR/r6UGPMh4wx5cBZwFPGmGfd1ycaY57u631lYHRod6ZtwdESTXQKaQOb0tTuFBGRw0lW7U5jzCzgW8D09K+x1r6nLw+11s4H5md4vQJ4f4bXXwJe6suzpP+8cJSwlra0ShpAQdCfakMOdEjTFhwiInI4yXZO2iPAH3G2wkj0cq3kmQ7tznjHClZB0J8Kcb5sm+d95PdrTpqIiBw+sg1pcWvtHwZ0JDJkeV3GpO24cACcPdK8ClcuK2makyYiIvku29rHE8aYLxljJhhjRnq/BnRkMmS0HwuVud3py1W7U3PSRETkMJJtJe1q95/fTnvNkpsNbWWQtW/BQWqfNJ9x2p/hgC9tTtrAjiM9pGlOmoiI5Ltsz+48YqAHIkOXN9csvd1ZEg7QEIlTEPTjFdAGco806BjSYgmLtRYzwNU7ERGRwdJju9MY8520jz/W6b2fDdSgZGhJb3dG4gn8PkNRyMn34YAv9f6Ab8HR6f4qpomISD7rbU7aVWkff7/Te5ce5LHIEOWFo+v+uoLbX9xCQcBHMOC8VhD052yftECnSl0soXlpIiKSv3oLaaabjzN9LnmqcxfT5zOE3NPOw8HcnTjQuZ2qeWkiIpLPegtptpuPM30ueapzOGqMxAl6IS3gH5QtOEB7pYmISH7rbeHAicaYBpyqWaH7Me7nBQM6MhkyMoWvUMAJaQU5rKQF/J1CmtqdIiKSx3oMadZaf64GIkNX5/B19LjSjpW01IkDuV04oHaniIjkswE+yEfyQXo4+uZFs3jm6+em5qR1rKTlbgsOULtTRETym0Ka9Cq9QlYcDmCMIRhor6S1b8ExsOPoEtJ0NJSIiOQxhTTpVXo2Kgo5HfCQ39uCw4dbVBv4zWy7LBzQnDQREclfCmnSq/Q2phfSOs5J83W5biB4lTTvMWp3iohIPlNIk16lh6/CoFtJ67C607tugMfhPsAbg9qdIiKSzxTSpFfp4as47CwIDqYWDuRun7RAp5Cm1Z0iIpLPFNKkVx0qaV3anb6ctTu9+xe4IS2mOWkiIpLHFNKkV5nmpLUvHPCnLRwY2HF4c9K8oKhKmoiI5DOFNOlVevgqCjrtzlCgvZLmy/E+aZqTJiIihwOFNOlVT+3OXM5J6xLS1O4UEZE8ppAmvep5Cw5fakL/gG9m681JC3khTZU0ERHJXwpp0qv0AlnXLTj8g9DudJ6tdqeIiOQzhTTpVfpxTF4g887uDAd97e3OgT5xoMsWHGp3iohI/lJIk15lqpAdPb6UmWNLGFYYxO/P1dmdzj8L1e4UEZHDQGCwByBDX6bwdd6sMTz/zfMBcrhwoL3FCmp3iohIflMlTXpleglfXhvSn6MD1ttXdyqkiYhI/lJIk37zKmi9hbn+GlMa5uSpwzlpynBAc9JERCS/qd0p/RbI0Zy0wpCf+V86h8qGCAAxtTtFRCSPqZIm/ebL0Zw0T8BdQaBjoUREJJ8ppEm/+XO0T1rn52lOmoiI5DOFNOm39tWduXmed8JBPKE5aSIikr8U0iRrR40tyfh6ritp3hw4VdJERCSfaeGAZOXFb13AqJJQxvdytQWHJ+DTsVAiIpL/FNIkK0eMLu72Pe84qBwV0lJtVW3BISIi+UztTum3QI7bncYYgn6jdqeIiOQ1hTTpNy+c5ard6T1LIU1ERPKZQpr0mz/H7U6AoM+nOWkiIpLXFNKk39y9ZXPW7gTw+43mpImISF5TSJN+87urLXPY7STgM8TU7hQRkTymkCb95h+EOWkBn4+E2p0iIpLHFNKk39xCGiaX7U6fIaZ2p4iI5DGFNOm3wCC0O4vDflraErl7oIiISI4NSkgzxnzMGPO2MSZpjJnbw3XDjTGPGmPeMcasN8aclctxSna8hQP+HFbSSguCNERiOXueiIhIrg1WJW0t8GHglV6u+y3wjLV2NnAisH6gByYHzlvVmct2Z1lBgMZIPGfPExERybVBORbKWrseev5L3RgzDDgPuMb9migQzcHw5ADl+oB1cCppW/c35+x5IiIiuTaU56QdAVQBdxtjVhpj7jDGdH+ApAya9gPWc/fMskJV0kREJL8N2F+rxpgFxpi1GX5dkeUtAsApwB+stScDzcD3enjetcaYZcaYZVVVVQfhO5BsDUYlrawgSENrDGu1DYeIiOSnAWt3Wmsv7OctyoFya+0S9/NH6SGkWWvnAfMA5s6dq7+5c8g/CHPSSguCxJOWSCxJYcifs+eKiIjkypBtd1pr9wK7jDFHuy+9F1g3iEOSbrRX0nL3zLJC5/8vGiIxtu1v5u7Xt+Xu4SIiIjkwWFtwfMgYUw6cBTxljHnWfX2iMebptEu/CtxvjFkDnAT8LPejld60z0nLbSUNoDES47EV5fzkiXVEYto3TURE8sdgre6cD8zP8HoF8P60z1cB3e6jJkPDYG3BAVDfGqeh1dkvrTWaoCCo1qeIiOSHIdvulEPHYLQ70ytpDe4qzxZV0kREJI8opEm/Bbx2Zw4racNSc9LiNEbaK2kiIiL5QiFN+s3nVdIGaU5aQ6tTSVNIExGRfKKQJv3WvgVH7p5Z5oa0htZ46gzPlqg2txURkfyhkCb95vfnvt1ZEPQR8BkaI7HUyQOtmpMmIiJ5RCFN+q0w6Gf6qCKOGJ27U7uMMZQVBmmMdFzdKSIiki8GZQsOyS9Bv4+Xvv3unD+3tCBAXWuMJrfN2aKQJiIieUSVNDlklRUE2VPXind8Z3q7s741RiyRHKSRiYiI9J9CmhyyhhcF2V7dkvo8vd158S0vc++i7YMwKhERkYNDIU0OWZNHFLK/qS31udfujCWS7Gtoo7y2dbCGJiIi0m8KaXLImjaq40IFr92Z+qfmqImIyCFMIU0OWdNHFXX4vDXacVPbZu2bJiIihzCFNDlkTR3ZsZLWGIlz76Lt2pJDRETygkKaHLKmpVXSRhQFeXljFf/1+NssfKcS0JYcIiJyaFNIk0NWcTjA6JIwAOPKCqhujgJQ2eAsJtAxUSIicihTSJND2vRRRRSF/JQWtO/LXNXkhTRV0kRE5NClEwfkkDZ7Qik1LVEKgv7Ua1WNEUAhTUREDm2qpMkh7buXzuaB/3cmRaH2kLa/yWl7qt0pIiKHMlXS5JBWWhCktCBIUSit3dmodqeIiBz6VEmTvJDe7qx3t+BoiydJJO1gDUlERKRfFNIkL6S3O9Op5SkiIocqhTTJC92FNG1oKyIihyqFNMkL6e3OdM0KaSIicohSSJO84FXSjOn4utqdIiJyqFJIk7xQ6FbSJg4r7PC62p0iInKoUkiTvHD85GEcO7GMEyYP6/C62p0iInKoUkiTvHDsxGE8df25jCsrAKAg6PzRblW7U0REDlEKaZJXSsLOprbeweva0FZERA5VCmmSV4rdkDaqOAQMfruzqU2VPBER6RuFNMkrJWFnAcEot5I2mO3OLVVNnHLT87yysSrnz16/p4HVu+py/lwRETl4FNIkr3SupA1mu3P+it1EE0lW7sx9WLrh4dV8/t5lxBLJnD+7O3vrI7y0oXKwhyEiQ4CO7MuOQprkFS+kFYcDhAO+jCEtEkvwH39dzsqdtb3eb19DhEjswIOetZZ/rt4NwOaqpgP++v7Y1xBh3Z4G9je1sWDdvpw+uzuRWIJr7l7K5+9dRmMkNtjD6dY7extoGMLjO9zEE8kuvx/NbXGWbK0e0n+OpGdvbK3m2P96hrW763u9Npm0PL66gr31kRyMbOhRSJO84i0cKAz5KQ4HaInGu/yHfsH6ffxr7V5+/dzG1Gt/fWMHH5/3BtF4e+WpLZ7gkltf4ef/eifjs/bWR7A28/8Nri6vZ1dNKyG/j82VBzektcUTPVbIXt7gtFdLwgEeWLrzoDyzvjXGP1ftpjmLOXYNkRgf+cMizvjZAo798TPc9MQ6fvSPtbyzt5FE0g5KZTGThkiMNeXOWGKJJN99dA2X3voqv3gm8++3tZZHl5fznUdXs3hLdbf3bYnGqW5q6/bPRrpXN1Vx85PreHJNRcbxJXupNtz2wiZ+/q932FXT0uuzBlL6vzfgjP2ZtXt5dVNVt39WX91UxS3Pb2RLD/8T873H3uKEG5/jmruXkkhaKhsjnPeLF7ly3hv8dsGmDte+tKGSy373Ks+9vbf/35AMGGstP3t6PZFYstffq2TS8rW/reL6B1fym+c35GiEQ4tCmuQVr5JWFPRTGPRT1djGlx9YwYk/eY6r71pKNJ5k/gqnwvXa5v2sq2gA4JHl5SzeWs1dr29L3euNrTXUtcR49u29Xf7CXb6jlrN+vpAF6zO375Zsdf4Sv+yECWytaiKRtFz/4Eq+8sCKrP7y7o61lqvvWsp5v3iRVd3MOVv4zj7GlxXw2XOm89rm/exvasv6/g2RGG9ur6Et3rF6+PsXN/O1h1Zx/i9fpKKutcd7LNlaw/IdtZw8ZQTnHDWau17fxqPLy/nCuUfgM7Bse03Gr4vGu1ZNDpZoPMn3H3uLl935gTurW7jitte54vbX2ba/mZufXMfflu1iTGmY1zdnDmCLt1bzrUdW88jycv63myDX3BbnrP95gVNvXsBH/7iYndXdh6dk0vKtR1Zzx2vb+PYjazoEnbqWKOf+74tc9ec3qGuJZvz6t8rr+dVzG/njy1v45sOrsv1RZFRe20J9a+affTJpM/6ZbYjEWL6jlnf2NnD8jc/yuXveZHddK4mk5ZN/XsJ1f13Op+9cyln/s5CfPrUuVflavauOO1/bxufvXcZvF27iittezxj+t1Y18diKco4cXcxLG6pYU17H35fvpro5yhGji3k9LShv2NvIl+9fwcZ9jVx733IWrh8aFeTuNLfFeXR5Of+3cBPX/mUZ9y/Z0W2YrWyIcN19y3lt0/4u77XFE2x1Q+7BnNrQGk30+N+p+pYD+/d0/Z6GVOfir2/sYE15PQVBHy9n+J7SrSqv44nVFZSEAyze2vHfy9rmKD+Y/xb//qee/z071PlvvPHGwR7DQTdv3rwbr7322sEehgyCprYYf1m8g/ceM5aywiD/WFXBlqpmzps1hlc2VlEcDvDg0p1cedpUNlc2krRw6tQR3PTk2wQDPpZsq2bWuFJmjC3hzle3saa8nqa2OJccO54xpWEWb6nmB/Pf4pWNVZTXtlIS9nPhMeOw1vKD+W8xf+VuLjt+AvNe2UoiabnipEk8+/Y+LjpmHP/91Do27G1kzsQyjhpb0mXsyaTl1gWb2LC3gdkTythe3UJB0Me8V7YyoijIqJIwL22s4vYXt2Ct5Yk1FXz6zGmEAj6qm9q467VtPLKsnMdX7+Gq06Zw+YkTuX/JTmaMKeaYCWX4Op+Z1cn6PQ1cfMsrPLBkJwY4e8bo1Li+/9hbTBlRxM7aFkrCAU4/YiT3L9nJXxbtIJ6wzBpXmrrPYyvKWbGzjse/8i4+dPIkWqJxzp4xmu9eOpvn1+2jpjnKR06d3OHZkViCT925hP9+Yj01LVHOnzUGYwzWWhaur2R4UYjCUMfzWaPxJH5fz99TMml5cUMlv39xC48sL2f1rjo+dcZUPnvPm+ypj9AWT/J2RT1Pr93LF849gvNnjeHx1RVcddpUSgoCHe713b+vIRJPct15M5i/cjeXnziRke7cR8/LG6t4ZHk5nzhjKos3V3P/kp1MGVnE0eNLSSQtizZXs2V/E1NHFvHm9hruXbSDD5wwgXV7Gjhv1hgmDXdOzLhv8Q6eX7+PysYIFXWtXHrchNT3U17byrDCILcu2MjmyiY+964jeGL1Hq48bQrhgD/1M/H+kr3ztW34fIbxwwoy/oxaonEu+NVL/GXxdkJ+H0WhALFEkpJwgHsXbecTdyzhdws3s3R7DUkL00cXsWJHLVf+6Q3uWbSd1zbtpzWWoKIuwhOrK9i6v5mF6yu5+YPH8fHTp1LfEuPvK3azcH0loYCPL963nJc2VHHE6GL+6/I5PL66grnTR1BWEOTKP73Bo8t3MbI4zJ9f3Up5bSsPfOFM/rpkB+OHFfDYit3MGlvKpceO54k1FVxz9nQKgn6+/tAqaltiPPO183jxnUoWb6nm4WW7iMaTnDR1ROp7rW+N4TMm9TMqr23hc/e8yVNv7eGMI0biM4Y3t9cwdWQR4ISVoN+HtZaHl+3ize01nJx2P3DmVhnApP37tW1/M396eStzp48g4HdqIQvX7yOetIwqCfO1h1by+5e2sHhrNc1tCR5fXUFrNMF5s8Z0+f255fmNPLy8nMdW7uaYCR3/23HL8xu5/sGVGGP4f/cuo7QgwENv7uKt8npmjClO/U8rOCvNb12wkdte2MwRo4uYOLzj6Sz1rTHueHUrDyzZyTf+tor61jjnH911PFuqmjj3Fy8ydWQR6/Y0pALn/Ut28u7ZY1Pfr2dvfYQP3v469y7ewcL1lTy8bBfnHz2G9x83gSfWVHD69JH8/qUt7K2PcMLk4ak/u8YYHn6znCXbqrnu/Bk8v24fHzllMj4Du2tb+dPLW3hw6S4q6loJB/2866jRGf98p3t8dQV+n+ny7621lrcrGhhVEsYAn7vnTV7dtJ9Ljh3f6z0Plp/85Cd7brzxxnmdXw9kuljkUDWyOEzAZxhTGubqs6cTTyRJWMsvP3oiF9/yMj//1zsUBH188bwjqWqM8OKGSs45ajRJC7f8+0n88eUtXHvfci6aM46VO2s5ddoIlu+o5Zm393LMhFJ+9vR63nLnUQR8htc2O/8n+NuFm3hw6S4AvntJCyt21nHezNGp/6DOe9UJbaNLQvxw/lsUhwKs21PPO3sa+czZ02mMxLh30XYWrK8k4DP8c3UFb5XX89FTJ/PQm7v41XMbuOezp/O7hZuYPKKQX3/sRK6c9wYPL9vFUWNL+NJfV9DYFscY+Pe5k/ne+2bj9xkmjyjkf5/ZwI/+sZbjJg1jfFkB15w9nTOOHAXAxn2NbNzXyAdOmMitCzbi9xlOmz6CB5fu5KvvmUko4GPlrlp217Vyw8WzeGzFbh5dXs4bW6tZsq2GUMDHG1urOWXacFbvquPiOeNZubOOYyaUpkLVDy+bk/r9mTttBA+9uYuPz3uD2pYo1793JnvrIzy8bBcb9jVyzozR3P36dj540iTKCoPc/OQ6Fr5TyUdPncyvPnYi4ASV//nXev6yeAf/9/GTubiH/5D+4tkN/PHlLQCcPHU4K3fW8aX7V7BqVx2//tiJPLaynNc3V3PkmGK+c+lsNuxtBGDJtmpaognmr9jNrz52Ire9uIlFW6r50WXHcPmJE7llwUZueX4jHzp5EqNLwpw4xfnL5aUNlRSH/Nx4+bF86YIZfO2hVXzjb6uYMKyAexZt58k1ewD44vlH0hSJUxj088PLjuGpt/bwu4WbiMQS/OcH5nDfGzs4ffpIZowt5vFVFWzf38y2/c08vrqC+St38+1Ljuafqyq4/ISJfPz0qcx7ZSuX/e41hhUGefbr5/GPlbu5+al1XDRnPH9fUc7okhD/+tp5jCkNd/kZPb9uH42ROLPHl3LTk+sACPl9nDtzNAvfqeS8WWM4akwJC9bv41uPrObu17dR3xqjOOxn0gjn9/36987k4jnj+PSdS3hgyU7ef/x4PnnGVIwxXDRnHK9uquLL96/gO4+u4YjRxfz5M6cydWQxSWsJB3y8umk/r2zcz7o9DUwbVcQX/rLM+Xfp0tnMHFfKCZOGcffr26lvjfHldx/FtFFF2Odh6bYafMb59/BHlx3DlJFFfPU9M7nhkdUAbK3awGUnTGRMaZjfPLeB372wmXDAxzcumsWnzpzGR/+wmOa2OElruXLeYk6ZOoJ/rqrglx89gR3VLcx7ZSsfOXUSLdEE/1zltKQnDS/k4mPHp4LbzU+tJ+j38d1Lj+bK06aSTFpueHgVK3bWYbF848JZ/PSp9dz3xg5KwgH+44IZ/GvtXv7jghlcd94MygoD/Ogfa7nz9W3sqm1hRFGI82aNYXhRkCVba3hw6U7ed9x4Nlc28ZvnNnLRMePw+QyxRJK/vVlO0sJvnt+Iz8CP//k2PgMWuP3FzVx3/gy+dcnR1DRH+eQdS9iwt4Gg38efXt7KvM+MTP0ZiMaTfOWBFby6aT/DCoOcNGU4d72+jW37m7hwzjg+cspkCtyj955cvYd40nLTk+uoaY7i95nUIoBYIsnvP3lKKrC+sbWanz29nmg8yZcumMGqXXV8+JTJ3HTFsazf08jvXtjMJ+5Ygs9A0kJFXSvDi4Lc9sJmvnbhLF7eWMkJk4bx/uPH87uFm1i8tZoX1leyYP0+fD7DR06ZxJ76CE+uqeA7lxyNMYbfPLeBVeX1fOHcIzh3ZnvIvOPVrdz81HrGloZ5/CvvYvywAu5dtJ07X9vGR06ZzC0LNvKlC2Zw9ozRvOhOGfncu6Zz7MSOp9jkmulP62Womjt3rl22bNlgD0MGyebKRqaPKu7yf3QPL9vFdx5dw2+vOokrTprEfW/s4D//sZYzjxzJmvJ6Vv34YpLW8oeXtvCXxduJJyy/ufIk7lm0jdc3V3POUaN4fXM115w9nYZIjKkji7h1wSauPmsa9y7ewXtnj2XhO5VcddoUHnpzFzd/8DguO34CJ//38/gMjCkNyWki2AAAE9BJREFU89fPn8HH/7wk1YIsDPppdRcmlIYDfObsadz12vbUawCnHzGSndUtlBQE2FzZxI8/MIfPvesIPvbHRayraCCWsBw5ppjbPnEKR4wu7lBduvnJddzx2jYuPGYsNc1Rdta0UtcS5aixJRw/aRgvbqhkf1OUa86ezj2LtvP1C2dy4pThfPbuNxlfVsCokhD1rTGqGttY9qMLWbi+kq//bRVFIT83XXEc1lq+/egapo0qYkd1C2cdOYq3dtfzoZMn8d8fPK7L780bW6u54eHVjC0L0xiJp+brnTh5GJ895wjec8xYTv/pAmaOLWXdngaCfsPUkUVU1EW49cqT8PsMS7fX8IeXtjC6JERDa5ybrjiWexZtpzkapy2WZFhhkDOOHMn6PY0s31HLx0+fynXnH8mk4YVc+JuX2V7dwrkzR3PvZ0/n8dUVfP1vq7j7mtN49+yxJJKWk256jvNmjmHJthr2N7VhDBjgi+fP4IaLZhHw+/jRP97ir2+0z/f7txMnMrwoyLNv7+XEycOZ95m5gNOSueTWV6hsdH6/v33J0Wza18iTa/ZgDFx2/ARuvepkLv+/11Lh33Pn1XMJ+H1cfddSSsKB1J5748rC7GtoY3RJiEeuO5sjRhfzb7e9xtrd9SQtXHjMWBasr6Q0HKCxLc6RY4opr20llkhy8pThfOvioykpCLBsey2NkTgvb6xkb32E1777HlaX17GzpoU/v7qVtbsb+I8LZqT+8rPW8vRbe/nKgyuwFh74whnMGFPC/W/s4AvnHUlpQZCWaJzWaIKRxaEOlSVwKjXzV5Tz7tljmTaqOPX6p+9cwoa9jdS2RPnoqVP44WXH8I2/reLocaXccPEsjDH8+rkN/N8LmznnqFHc89nTSVrLiT95DoBILMmk4YUsvOF8CoJ+4okkv39pC0eNLeGrD65kXGmYI8YU8/rmai45dhyJpDMvdcaYYrZUNfPIdWcRiyf5xB1LAAj6DbGETf25XF1ej8/ANy6cxbPr9rJxbxMXHTuOklCAvy3bxZlHjqQtnmRdRQMfPmUSi7ZUs6O6hVnjSthU2cTwwiC1LTGuPmsar2+pZnNlE6NLQrz07Xen5tA2tcX55B1LaIzEqGmOUpfWTvQZeOKr72Ljvka+8bfV/OnTp+Izhqff2sP8lbv5wrlH8Ob2Wv7r8jnc9fp2PnjSRGaMKeGXz23gqTV7ePS6s3hg6U6eWF3BHVefxmubqrhn0XZ+9bETaW5LsGjLfp5Zu5d40vLzDx/PVadPJZG0/O8z7/Dc23vZXt3CpOGFfPndR1EYcgLe7rpWGiNxjhxdzJSRRZQWBDhmQhm/fHYDN1w0iw37GvnUmdP45B1LGFEU5L8uP5bLT5zY4c+DtZYHlu6kOBTg3Jmj+dVzG3hw6S4CPsOwwiDVzU6b/6vvOYpvXDiL0366gMkji3irvI5QwEcsYXnuG+exfHst3/n7Gs6dOZojRxdz7+Idqf+ufuSUyeyqbeHiOeO4+an1nDtzNCt2OG3X044YySsbq/CmfQZ8BguMLysgkbQ0R+McPa6UW686ickjirr8t+xgM8Yst9bO7fK6QpocTiobI4wtddo+u2paOPcXLwJw1WlT+PlHTsj4NS3ROL9dsIl/rNpNUSjAM18/l3DAz+bKJi78zcsATnXl30/kqnlvsHxnLdbC09efy5yJZTy4dCf/8/R6PnnmNL576Wy27W/mtc37uWDWmNTk/hljSrjg6DEUBP08/OYuVpfXEfT7uGfRdu757Gks2eYEk6DfsOQHFzKyOMSqXXX86eUtjCoJ8a2Lj2Z4UajL2Jva4qzcWfv/27vz6KjKNI/j3ycJhLAY2XcEWRRcQEFApBsEPA3u9nFjtFt7GO2Z1hnbaXV0esa2nbZPz7HVAe32qO3uuLXd7cLYCyICioKoKJvscWFJAgESslWq6pk/7k0oICFBSeoW/j7n1Endt25dnuJJpZ563/fel/GDumBmlFbV8Ks/f8rWXZUsXLedvFbZ9OvclpVbShne92ienjGa9q1zuCqcqF0dTwYfUGcNYdzALsTiSR57ZxPfOaEHA7q0Y1dFjJG/eINE0pkytDsL1hYTSyS555LhBwxp7q+sqob/mr2K0wd25qJT9u574wvL+NNHmxneJ59HrhrF2m17uPLRxfs8d/rovtw6dSjTH3mPVVtL6dAmh28P7kpuThabdpSzvnAPA7u1Z8rQbvxwwkBahQX7zvIY5bE4vfLzyMoKCo+CHRUM6LK3aLjztVV1cxNnjB/Amm1l3PSd4xgR9pbV2ra7iq27K3ll2RaeW/I5WWZU1iTqPuhqLS0o4bWPt3D+iF6MPKYTO/ZUc979b3NKv4788qKTyG/bil+HPX73XDqcv60q5Ptjj2HMsZ2JxZOM+sUcSqvi3H7uMAZ1a8/xPTpw/5vr+fvxA+ri3rKrkup4kmufWsq6oj18a3AXHph+KjPnrmP66L6UVsWZv6aIZ5d8Ue8cxR9OOJbbpg3d5/dmxebdjA17XFO9+P4XFJVVcf2kwQfNb1M9smAjd72+muO6d+Dpfxhd9/5MVVhaxe8WbuT6SYPJz2sFBD0jq7eWMaJvPucP701+21YHPO+VZZv58/JtbN5VSd9Oedx76Qhyc7K449WVdV+sHr36NACuf/ZDFqwt5vEfjOaJRQVcOaZf3bD+wK7tOX1gZwpLq3jwrQ28vGwzuypquHpcf24/dxjbSquYcu98KmIJJgzpykm98/nHiQN5eMFGCraXc8moPnxrcFdi8SSbtpeTn9eqweHnZNL568ptlFTEOPfkXlTE4vTMzyOeSDJ15kJ2lscoraqhJuH06ZjHWzdNPOALKQR/tybfM5+aRJLte2L8aOJAbpl6PKu2lHL2rIV1++XmZPF3Y/oxflAXJg/tvs8x3J1FG3bw0z8tpyBl3td/nDOU6niSyUO7cXyPo4DgTNxpMxeyLvzilZuTRdKdt/9tUt1yfQdTHU9w8YPvUrC9nDd+MoG/rNjGrLnreHrGGIb1Ooqn3i3g9ldWYgZzbvw2YAzq1p7dFTWcPWshZvDlzkr6dMxj9j+P5+aXPmHOqsK6LzjH9+jAy9edwYbiPTy5qIBPvtxN76PzuHhkH+6ds5a7LjqJu15fjbvzr2cNoaismv98eQU5Wcai2ybX/d41FxVpIvW46LfvAPDcNWPruvMbUvteqe0hqP0meEyndpwxqDNmxrw1RTz+TgEj+uTz4ylDyAp7tWrnrWQ1Mocq1Z7qOPPXFHP2ST3YUFzOlHvnM+3EHjx45civ8EoPVLC9nIQ7Hdu2ZtGG7Uw9oUe9f+wb871HF7NySykLbzmTjz7fxUMLNnDfZSPqluY6VGsLy5g5dx0/O28Y3ToE32on3D2P9rk5DOt5FJ+VVPDMjDHktc6mpDzGL19fzfTR/Rh5TMfGD94E1fEEVzyymFgiySvXnXFAj1B93IOCdmnBTsYe26nR/8faOTe1KmMJCkur6J9SLNZ64p1NlMcSXHfmoEbjmLemiEcXbmLW9FMOmHcDwWT/xRtLqIjFGXtsZ9q2zq7rJe7Qpnk/hBpSVZNg3qdFTB7andY5LXMuWzyR5On3PmPqiT3omR/MzUoknV0VsboLYR9MdTzBhqJyhvbsUJfHxRt3kHQ4feCBhe3hsmZbGRf85m26tM/lD/80jg5tcmjbuuFZSx98VsIDb64nnnQe+t7Iun3vfG0VPfJzmXpCT9rlZjf6mqtqEnxeUsHawjJe+uBL7r54eL1D5+8XlDDzjXX07dSW55Z8ziUj+3B3OE2hKcqr4+yurDlgvhwE75k7Xl2JmXHH+SfU+/z3C0ro2j6X/l3akUw62/dUE0skmTV3HT+aOKje99fBbN5VyZJNO/b5EtlcVKSJ1KOqJkHr7KxDKp7S5cWlXzC6f6dD/kPT3IrKqqioTjRrXCXlMdrlZpObc/BC+nBJJp1YItlo4S7S0tYWlpGf16pJvVPpUlWT4L45a7lqXP96Cy45kIo0ERERkQhqqEjTddJEREREIkhFmoiIiEgEqUgTERERiaC0FGlmdomZrTSzpJkdMAabst+N4X4rzOw5M4vuTEkRERGRwyhdPWkrgO8CCxrawcx6A/8CjHL3E4Fs4PKWCU9EREQkvdKyLJS7rwaacv2hHCDPzGqAtsCWZg5NREREJBIiOyfN3TcDvwY+B7YCu939bw3tb2bXmtlSM1taXFzcUmGKiIiINItmK9LM7I1wLtn+twua+PyOwAXAAKAX0M7Mrmxof3d/2N1Hufuorl27NrSbiIiISEZotuFOd5/yNQ8xBdjk7sUAZvZHYBzwzNeNTURERCTqIjvcSTDMOdbM2loweW0ysDrNMYmIiIi0iHRdguMiM/sSOB34PzP7a9jey8xeB3D3xcBLwIfA8jDWh9MRr4iIiEhL09qdIiIiImmktTtFREREMoiKNBEREZEIUpEmIiIiEkEq0kREREQiSEWaiIiISASpSBMRERGJIBVpIiIiIhGkIk1EREQkgppt7c502lhczmUPvZvuMERERES+MvWkiYiIiESQloUSERERSSMtCyUiIiKSQVSkiYiIiESQijQRERGRCFKRJiIiIhJBKtJEREREIkhFmoiIiEgEqUgTERERiSAVaSIiIiIRpCJNREREJIJUpImIiIhEkIo0ERERkQg6ItfuNLNi4LNm/Ce6ANub8fjS/JTDzKb8ZT7lMPMph4fPMe7edf/GI7JIa25mtrS+hVAlcyiHmU35y3zKYeZTDpufhjtFREREIkhFmoiIiEgEqUj7ah5OdwDytSmHmU35y3zKYeZTDpuZ5qSJiIiIRJB60kREREQiSEXaITCzqWa2xszWm9mt6Y5H6mdmj5lZkZmtSGnrZGZzzGxd+LNj2G5mNivM6Sdmdmr6IpdaZtbXzOaZ2SozW2lmN4TtymMGMLM2ZrbEzD4O8/fzsH2AmS0O8/SCmbUO23PD7fXh4/3TGb/sZWbZZvaRmc0Ot5XDFqQirYnMLBv4DTANGAZMN7Nh6Y1KGvAEMHW/tluBue4+GJgbbkOQz8Hh7VrgwRaKUQ4uDvzE3YcBY4Hrwveb8pgZqoFJ7j4cGAFMNbOxwH8D97n7IGAnMCPcfwawM2y/L9xPouEGYHXKtnLYglSkNd1oYL27b3T3GPA8cEGaY5J6uPsCoGS/5guAJ8P7TwIXprQ/5YH3gKPNrGfLRCoNcfet7v5heL+M4EOiN8pjRgjzsCfcbBXeHJgEvBS275+/2ry+BEw2M2uhcKUBZtYHOAf4XbhtKIctSkVa0/UGvkjZ/jJsk8zQ3d23hve3Ad3D+8prxIXDJqcAi1EeM0Y4TLYMKALmABuAXe4eD3dJzVFd/sLHdwOdWzZiqcf/ALcAyXC7M8phi1KRJt84HpzSrNOaM4CZtQf+APzY3UtTH1Meo83dE+4+AuhDMBJxfJpDkkNgZucCRe7+Qbpj+SZTkdZ0m4G+Kdt9wjbJDIW1w1/hz6KwXXmNKDNrRVCg/a+7/zFsVh4zjLvvAuYBpxMMQ+eED6XmqC5/4eP5wI4WDlX2dQZwvpkVEEzvmQTMRDlsUSrSmu59YHB4Zktr4HLg1TTHJE33KnBVeP8q4JWU9u+HZweOBXanDKdJmoRzWR4FVrv7vSkPKY8ZwMy6mtnR4f084CyCeYXzgIvD3fbPX21eLwbedF3EM63c/TZ37+Pu/Qk+79509ytQDluULmZ7CMzsbIIx+mzgMXe/K80hST3M7DlgItAFKAR+BrwMvAj0Az4DLnX3krAYeIDgbNAK4AfuvjQdccteZjYeWAgsZ+98mH8nmJemPEacmZ1MMIk8m6Az4EV3v9PMjiXolekEfARc6e7VZtYGeJpg7mEJcLm7b0xP9LI/M5sI3OTu5yqHLUtFmoiIiEgEabhTREREJIJUpImIiIhEkIo0ERERkQhSkSYiIiISQSrSRERERCJIRZqIZDwzS5jZspTbrY0/q97jvGVmow53fE34dy8MF5AXEamT0/guIiKRVxkuQZSpLgRmA6vSHYiIRId60kTkiGRmU83s9ynbE81sdnj/QTNbamYrzeznTTjWaWa2yMw+NrMlZtbBzNqY2eNmttzMPjKzM8N9rzazB1KeOzu8GChmtsfM7gqP856ZdTezccD5wN1hL+DAw/xfISIZSkWaiBwJ8vYb7rwMeAMYY2btwn0uI7hSOsBP3X0UcDIwIbxCfr3CZeBeAG5w9+HAFKASuI5gnfeTgOnAk+FV1w+mHfBeeJwFwDXuvohgSZ2b3X2Eu2/4Cq9fRI5AKtJE5EhQGRY4tbcX3D0O/AU4L1zw+Rz2rjN4qZl9SLCszQnAweaDHQdsdff3Ady9NDz2eOCZsO1TgmWqhjQSZ4xgWBPgA6D/Ib5OEfkG0Zw0ETmSPQ9cT7CW4FJ3LzOzAcBNwGnuvtPMngAa6wE7FHH2/QKceuyalEWnE+hvsIgchHrSRORINh84FbiGvUOdRwHlwG4z6w5Ma+QYa4CeZnYaQDgfLYdgAfgrwrYhBIu+rwEKgBFmlmVmfYHRTYizDOhwCK9LRL4BVKSJyJFg/zlpvwJw9wTB8OK08Cfu/jHBMOenwLPAOwc7sLvHCOaz3W9mHwNzCHrHfgtkmdlygjlrV7t7dXi8TQRnas4CPmxC/M8DN4cnIOjEAREBwPb2vIuIiIhIVKgnTURERCSCVKSJiIiIRJCKNBEREZEIUpEmIiIiEkEq0kREREQiSEWaiIiISASpSBMRERGJIBVpIiIiIhH0/6V256yT8iASAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x360 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=[10,5])\n",
    "plt.plot(counts1, values1)\n",
    "plt.xlabel('Eval count')\n",
    "plt.ylabel('Energy')\n",
    "plt.title('Convergence with noise');\n",
    "plt.axhline(ref)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<h3>Version Information</h3><table><tr><th>Qiskit Software</th><th>Version</th></tr><tr><td>Qiskit</td><td>0.20.0</td></tr><tr><td>Terra</td><td>0.15.1</td></tr><tr><td>Aer</td><td>0.6.1</td></tr><tr><td>Ignis</td><td>0.4.0</td></tr><tr><td>Aqua</td><td>0.7.5</td></tr><tr><td>IBM Q Provider</td><td>0.8.0</td></tr><tr><th>System information</th></tr><tr><td>Python</td><td>3.8.2 (default, May  6 2020, 02:49:43) \n",
       "[Clang 4.0.1 (tags/RELEASE_401/final)]</td></tr><tr><td>OS</td><td>Darwin</td></tr><tr><td>CPUs</td><td>6</td></tr><tr><td>Memory (Gb)</td><td>16.0</td></tr><tr><td colspan='2'>Thu Oct 15 11:08:51 2020 EDT</td></tr></table>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<div style='width: 100%; background-color:#d5d9e0;padding-left: 10px; padding-bottom: 10px; padding-right: 10px; padding-top: 5px'><h3>This code is a part of Qiskit</h3><p>&copy; Copyright IBM 2017, 2020.</p><p>This code is licensed under the Apache License, Version 2.0. You may<br>obtain a copy of this license in the LICENSE.txt file in the root directory<br> of this source tree or at http:https://www.apache.org/licenses/LICENSE-2.0.<p>Any modifications or derivative works of this code must retain this<br>copyright notice, and modified files need to carry a notice indicating<br>that they have been altered from the originals.</p></div>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import qiskit.tools.jupyter\n",
    "%qiskit_version_table\n",
    "%qiskit_copyright"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}