Project_QC.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "concrete-musician",
   "metadata": {},
   "source": [
    "# Import some libraries."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "cardiovascular-contract",
   "metadata": {},
   "outputs": [],
   "source": [
    "from qiskit.opflow import Z, I, X, Y\n",
    "from qiskit.opflow import CircuitStateFn, StateFn, CircuitSampler, PauliExpectation, ListOp, AerPauliExpectation\n",
    "from qiskit.providers.aer import AerSimulator\n",
    "from qiskit.providers.aer.noise import NoiseModel\n",
    "from qiskit.utils import QuantumInstance\n",
    "from qiskit.circuit import QuantumCircuit, QuantumRegister, Parameter, ParameterVector\n",
    "from qiskit.algorithms.optimizers import ADAM, SPSA, GradientDescent\n",
    "from qiskit.circuit.library import TwoLocal, EfficientSU2\n",
    "from qiskit.circuit.library.n_local import RealAmplitudes\n",
    "from qiskit import Aer, execute, transpile, IBMQ\n",
    "from qiskit.providers.aer.noise import NoiseModel\n",
    "from qiskit.providers.aer.noise import QuantumError\n",
    "from qiskit.providers.aer.noise import amplitude_damping_error\n",
    "from qiskit.providers.aer.noise import depolarizing_error\n",
    "from qiskit.algorithms import VQE, NumPyMinimumEigensolver\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "handy-uncertainty",
   "metadata": {},
   "source": [
    "Let's define the number of qubits $N$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "ordered-separation",
   "metadata": {},
   "outputs": [],
   "source": [
    "num_qub = 4"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "planned-scottish",
   "metadata": {},
   "source": [
    "Let's define the transverse field Ising Hamiltonian on a 1D lattice: \n",
    "\\begin{equation}\n",
    "H=J_z\\sum_{i=1}^{N} Z_{i} Z_{i+1} + h\\sum_{i=1}^{N} X_{i},\n",
    "\\end{equation}\n",
    "where $Z_{i}$ and $X_{i}$ are respectively the $Z$ and $X$ Pauli operators acting on the $i$-th qubit. The transverse-field Ising model is the quantum version of the classical Ising model. The Hamiltonian is composed by two types of terms: \n",
    "\n",
    "- nearest neighbour interactions determined by the alignment or anti-alignment of spin projections along the $z$ axis. The energy of these interactions is quantified by $J$. If $J<0$ we have the tendency of alignment between spins and so ferromagnetic coupling, while if $J>0$ the tendency is of anti-alignment and thus the coupling is anti-ferromagnetic. \n",
    "- an external magnetic field along the x axis which creates an energetic bias for the alignment of spins projections along the $x$ axis. The strenght of the magnetic field is determined by $h$. \n",
    "\n",
    "We have considered periodic boundary conditions in the system, namely $N+1 \\equiv 1$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "decent-relay",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.0 * ZZII\n",
      "+ 1.0 * IZZI\n",
      "+ 1.0 * IIZZ\n",
      "+ 1.0 * ZIIZ\n",
      "+ 1.0 * XIII\n",
      "+ 1.0 * IXII\n",
      "+ 1.0 * IIXI\n",
      "+ 1.0 * IIIX\n"
     ]
    }
   ],
   "source": [
    "J = np.ones(num_qub)\n",
    "h = np.ones(num_qub)\n",
    "\n",
    "def create_Hamiltonian(num_qub, J, h):\n",
    "    H = 0\n",
    "    \n",
    "    if(num_qub == 2):\n",
    "        return J[0]*(Z^Z) + h[0]*(X^I) + h[1]*(I^X)\n",
    "    \n",
    "    for i in range(num_qub):\n",
    "        if(i==0 or i==num_qub-1):\n",
    "            pauli_str = J[i]*Z\n",
    "        else: \n",
    "            pauli_str = J[i]*I\n",
    "        for j in range(num_qub):\n",
    "            if(j!=i):\n",
    "                if(j==i+1 or j==i-1):\n",
    "                    pauli_str ^= Z\n",
    "                else:\n",
    "                    pauli_str ^= I\n",
    "        H += pauli_str\n",
    "        \n",
    "    for i in range(num_qub):\n",
    "        if(i==0):\n",
    "            pauli_str = h[i]*X\n",
    "        else: \n",
    "            pauli_str = h[i]*I\n",
    "        for j in range(num_qub-1):\n",
    "            if(j==i-1 and i!=0):\n",
    "                pauli_str ^= X\n",
    "            else:\n",
    "                pauli_str ^= I\n",
    "        H += pauli_str\n",
    "\n",
    "    return H\n",
    "\n",
    "H = create_Hamiltonian(num_qub, J, h)\n",
    "print(H)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "interstate-cassette",
   "metadata": {},
   "source": [
    "The ansatz considered is an hardware-efficient ansatz and consists in the alternation of rotations around the $Y$ axis on all the qubits followed by $CX$ (to create entanglement) on all possible qubit pairs. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "lasting-glory",
   "metadata": {},
   "outputs": [],
   "source": [
    "num_lay = 2\n",
    "params = ParameterVector('θ', num_qub*num_lay)\n",
    "entanglement = 'linear'\n",
    "ansatz_gates = ['ry', 'rz']\n",
    "\n",
    "def create_ansatz(num_qub, num_lay, params):\n",
    "    ansatz = QuantumCircuit(num_qub)\n",
    "    for i in range(num_lay):\n",
    "        for j in range(num_qub):\n",
    "            ansatz.ry(params[j+i*num_qub], j)\n",
    "        ansatz.barrier()\n",
    "        if(num_qub > 2):\n",
    "            for j in range(num_qub-1):\n",
    "                    #if(j == num_qub-1):\n",
    "                    #    ansatz.cx(0, j)\n",
    "                    #else:\n",
    "                    ansatz.cx(j, j+1)\n",
    "        else: \n",
    "            ansatz.cx(0, 1)\n",
    "        ansatz.barrier()\n",
    "    return ansatz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "level-canada",
   "metadata": {},
   "outputs": [],
   "source": [
    "#ansatz = EfficientSU2(num_qubits = num_qub, su2_gates = ansatz_gates, entanglement = entanglement, reps = num_lay, insert_barriers = True)\n",
    "#ansatz = ansatz.decompose()\n",
    "#ansatz.draw('mpl', filename='./ansatz.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "infrared-vessel",
   "metadata": {},
   "source": [
    "Define a class which implement the VQE. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "intended-flood",
   "metadata": {},
   "outputs": [],
   "source": [
    "provider         = IBMQ.load_account()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "loose-photographer",
   "metadata": {},
   "outputs": [],
   "source": [
    "class VQE:\n",
    "    def __init__(self, ansatz, H, noise = 'no', prob = 0):\n",
    "        self.H = H\n",
    "        self.ansatz = ansatz\n",
    "        if(noise == 'no'):\n",
    "            noise_model = None \n",
    "        if(noise == 'yes'):\n",
    "            backend = provider.get_backend('ibmq_manila')\n",
    "            noise_model = NoiseModel.from_backend(backend)\n",
    "        if(noise == 'model'):\n",
    "            noise_model = NoiseModel()\n",
    "            err_1gate = amplitude_damping_error(prob)\n",
    "            err_2gates = err_1gate.tensor(err_1gate)\n",
    "            noise_model.add_all_qubit_quantum_error(err_1gate, [\"u1\", \"u2\", \"u3\"])\n",
    "            noise_model.add_all_qubit_quantum_error(err_2gates, [\"cx\"])\n",
    "            \n",
    "        self.simulator = Aer.get_backend('qasm_simulator')\n",
    "        self.shots = 10000\n",
    "        #self.shots = 1\n",
    "        self.lr = 0.1\n",
    "        \n",
    "        self.quantum_instance = QuantumInstance(backend = self.simulator, shots = self.shots, noise_model = noise_model)\n",
    "        #self.optimizer = ADAM(maxiter=1, tol=1e-06, lr=self.lr, beta_1=0.9, beta_2=0.99, amsgrad=False)\n",
    "        self.optimizer = GradientDescent(maxiter=1, learning_rate=self.lr)\n",
    "        \n",
    "        self.expectationer = PauliExpectation()\n",
    "        #self.expectationer = AerPauliExpectation()\n",
    "        \n",
    "\n",
    "    def expectation_value(self, parameters):\n",
    "        \n",
    "        ansatz_dum = self.ansatz.bind_parameters(parameters)\n",
    "        #ansatz_dum = create_ansatz(num_qub, num_lay, parameters)\n",
    "        ansatz_state = CircuitStateFn(ansatz_dum)\n",
    "        \n",
    "        measurable_expression = StateFn(self.H, is_measurement=True).compose(ansatz_state) \n",
    "        \n",
    "        expectation = self.expectationer.convert(measurable_expression)  \n",
    "\n",
    "        sampler = CircuitSampler(self.quantum_instance).convert(expectation)\n",
    "\n",
    "        return sampler.eval().real\n",
    "            \n",
    "    \n",
    "    def gradient(self, parameters):\n",
    "        \n",
    "        gradients = np.zeros_like(parameters)\n",
    "        for i in range(len(parameters)):\n",
    "           \n",
    "            shift = np.zeros_like(parameters)\n",
    "            shift[i] = np.pi/2\n",
    "            gradients[i] = (self.expectation_value(parameters + shift) - self.expectation_value(parameters - shift))/2\n",
    "        \n",
    "        return gradients\n",
    "    \n",
    "    def update(self,parameters):\n",
    "\n",
    "        parameters, loss, it = self.optimizer.optimize(parameters.size,\n",
    "                        lambda param: self.expectation_value(param),\n",
    "                        gradient_function= lambda param: self.gradient(param), initial_point=parameters)\n",
    "        \n",
    "        #loss = self.expectation_value(parameters)\n",
    "        #gradients = self.gradient(parameters)\n",
    "        #parameters -= self.lr*gradients \n",
    "        \n",
    "        return loss, parameters\n",
    "  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "graphic-reserve",
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib notebook\n",
    "plt.ion()\n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(111)\n",
    "fig.show()\n",
    "fig.canvas.draw()\n",
    "\n",
    "num_lay = 4 # 1 for 2 qubits, 2 for 3 qubits and 4 for 4 qubits\n",
    "entanglement = 'linear'\n",
    "ansatz_gates = ['ry', 'rz']\n",
    "ansatz = EfficientSU2(num_qubits = num_qub, su2_gates = ansatz_gates, entanglement = entanglement, reps = num_lay)\n",
    "np.random.seed(1234)\n",
    "parameters = np.pi * np.random.normal(size = len(ansatz.parameters))\n",
    "\n",
    "noise = 'yes'\n",
    "prob = 0.05\n",
    "vqe = VQE(ansatz, H, noise = noise, prob = prob)\n",
    "epochs = 500\n",
    "losses = []\n",
    "exact_energy = min(np.real(np.linalg.eig(H.to_matrix())[0])) \n",
    "print(\"Exact ground energy: \", exact_energy)\n",
    "\n",
    "for i in range(epochs):\n",
    "    \n",
    "    loss, parameters = vqe.update(parameters)\n",
    "    losses.append(loss)\n",
    "    ax.clear()\n",
    "    ax.plot(losses, 'b.', label = 'VQE energy')\n",
    "    ax.plot(exact_energy*np.ones_like(range(epochs)),'k--', label = 'Exact ground energy')\n",
    "    ax.legend()\n",
    "    ax.set_title(\"VQE energy vs epochs and exact ground energy\")\n",
    "    ax.set_ylabel(\"Energy (a.u.)\")\n",
    "    ax.set_xlabel(\"Epoch\")\n",
    "    fig.canvas.draw()\n",
    "   \n",
    "print(\"Ground energy from VQE: \", loss)\n",
    "print(\"Percentual accuracy: \", 100*np.abs((loss-exact_energy)/exact_energy), \"%\")\n",
    "np.savetxt('ground_%iqub_%ilay_ryrz_%s_%sNoise_losses' % (num_qub, num_lay, entanglement, noise), losses)\n",
    "np.savetxt('ground_%iqub_%ilay_ryrz_%s_%sNoise_fin_params'  % (num_qub, num_lay, entanglement, noise), parameters)\n",
    "plt.savefig(\"ground_%iqub_%ilay_ryrz_%s_%sNoise.pdf\"  % (num_qub, num_lay, entanglement, noise))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "banned-insulin",
   "metadata": {},
   "source": [
    "Plotting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "retired-grave",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(111)\n",
    "num_lay = 4\n",
    "entanglement = 'linear'\n",
    "noise = 'yes'\n",
    "losses = np.loadtxt('Project_QC/Ising_model/shot-based_backend/with_noise/ground_%iqub_%ilay_ryrz_%s_%sNoise_losses' % (num_qub, num_lay, entanglement, noise))\n",
    "exact_energy = min(np.real(np.linalg.eig(H.to_matrix())[0]))\n",
    "ax.plot(losses, 'b.', label = 'VQE energy')\n",
    "ax.plot(exact_energy*np.ones_like(range(len(losses))),'k--', label = 'Exact ground energy')\n",
    "ax.legend()\n",
    "ax.set_title(\"VQE energy vs epochs and exact ground energy\")\n",
    "ax.set_ylabel(\"Energy (a.u.)\")\n",
    "ax.set_xlabel(\"Epoch\")\n",
    "fig.canvas.draw()\n",
    "print(\"Exact ground energy: \", exact_energy)\n",
    "print(\"VQE ground energy: \",losses[-1])\n",
    "print(\"Accuracy: \", 100*np.abs((losses[-1]-exact_energy)/exact_energy), \"% \\n\")\n",
    "plt.savefig('Project_QC/Ising_model/shot-based_backend/with_noise/ground_%iqub_%ilay_ryrz_%s_%sNoise.pdf'  % (num_qub, num_lay, entanglement, noise))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "trying-strap",
   "metadata": {},
   "source": [
    "## Overlap based method for excited states and excited energies\n",
    "Once we have obtained an approximation of the ground state we apply the overlap based method to find some low excited energies and corresponding states. \n",
    "### $1^{st}$ excited energy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "catholic-recall",
   "metadata": {},
   "outputs": [],
   "source": [
    "num_lay = 4\n",
    "entanglement = 'linear'\n",
    "ansatz_gates = ['ry', 'rz']\n",
    "noise = 'yes'\n",
    "ansatz = EfficientSU2(num_qubits = num_qub, su2_gates = ansatz_gates, entanglement = entanglement, reps = num_lay)\n",
    "parameters = np.loadtxt('Project_QC/Ising_model/shot-based_backend/with_noise/ground_%iqub_%ilay_ryrz_%s_%sNoise_fin_params' % (num_qub, num_lay, entanglement, noise))\n",
    "\n",
    "ground_circ = ansatz.bind_parameters(parameters)\n",
    "ground_state = CircuitStateFn(ground_circ)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "earlier-production",
   "metadata": {},
   "outputs": [],
   "source": [
    "class first_exc_VQE(VQE):\n",
    "    \n",
    "    def expectation_value(self, parameters):\n",
    "        alpha = 2\n",
    "        ansatz_dum = self.ansatz.bind_parameters(parameters)\n",
    "\n",
    "        ansatz_state = CircuitStateFn(ansatz_dum)\n",
    "        measurable_expression = StateFn(self.H, is_measurement=True).compose(ansatz_state) \n",
    "        overlap = alpha * np.abs((~ground_state @ ansatz_state).eval())\n",
    "        expectation = self.expectationer.convert(measurable_expression)  \n",
    "\n",
    "        sampler = CircuitSampler(self.quantum_instance).convert(expectation)\n",
    "        result = sampler.eval().real + overlap\n",
    "        return result\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "frequent-clinic",
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib notebook\n",
    "plt.ion()\n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(111)\n",
    "fig.show()\n",
    "fig.canvas.draw()\n",
    "\n",
    "num_lay = 4\n",
    "entanglement = 'linear'\n",
    "ansatz_gates = ['ry', 'rz']\n",
    "ansatz = EfficientSU2(num_qubits = num_qub, su2_gates = ansatz_gates, entanglement = entanglement, reps = num_lay)\n",
    "parameters = np.pi * np.random.normal(size = len(ansatz.parameters))\n",
    "\n",
    "noise = 'yes'\n",
    "vqe = first_exc_VQE(ansatz, H, noise = noise)\n",
    "epochs = 500\n",
    "losses = []\n",
    "exact_energy = np.sort(np.real(np.linalg.eig(H.to_matrix())[0]))[1]\n",
    "\n",
    "print(\"Exact 1st excited energy: \", exact_energy)\n",
    "\n",
    "for i in range(epochs):\n",
    "    \n",
    "    loss, parameters = vqe.update(parameters)\n",
    "    losses.append(loss)\n",
    "    ax.clear()\n",
    "    ax.plot(losses, 'b.', label = 'VQE energy')\n",
    "    ax.plot(exact_energy*np.ones_like(range(epochs)),'k--', label = r'Exact $1^{st}$ excited energy')\n",
    "    ax.legend()\n",
    "    ax.set_title(r\"VQE energy vs epochs and exact $1^{st}$ excited energy\")\n",
    "    ax.set_ylabel(\"Energy (a.u.)\")\n",
    "    ax.set_xlabel(\"Epoch\")\n",
    "    fig.canvas.draw()\n",
    "    \n",
    "print(\"1st excited energy from VQE: \", loss)\n",
    "print(\"Percentual accuracy: \", 100*np.abs((loss-exact_energy)/exact_energy), \" %\")\n",
    "np.savetxt('Project_QC/Ising_model/shot-based_backend/with_noise/1st_exc_%iqub_%ilay_ryrz_%s_%sNoise_losses' % (num_qub, num_lay, entanglement, noise), losses)\n",
    "np.savetxt('Project_QC/Ising_model/shot-based_backend/with_noise/1st_exc_%iqub_%ilay_ryrz_%s_%sNoise_fin_params'  % (num_qub, num_lay, entanglement, noise), parameters)\n",
    "plt.savefig(\"Project_QC/Ising_model/shot-based_backend/with_noise/1st_exc_%iqub_%ilay_ryrz_%s_%sNoise.pdf\"  % (num_qub, num_lay, entanglement, noise))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "accepting-detail",
   "metadata": {},
   "source": [
    "Plotting"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "powered-seminar",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(111)\n",
    "num_lay = 4\n",
    "entanglement = 'linear'\n",
    "ansatz_gates = ['ry', 'rz']\n",
    "noise = 'yes'\n",
    "losses = np.loadtxt('Project_QC/Ising_model/shot-based_backend/with_noise/1st_exc_%iqub_%ilay_ryrz_%s_%sNoise_losses' % (num_qub, num_lay, entanglement, noise))\n",
    "exact_energy = np.sort(np.real(np.linalg.eig(H.to_matrix())[0]))[1]\n",
    "ax.plot(losses, 'b.', label = 'VQE energy')\n",
    "ax.plot(exact_energy*np.ones_like(range(len(losses))),'k--', label = r'Exact $1^{st}$ excited energy')\n",
    "ax.legend()\n",
    "ax.set_title(r\"VQE energy vs epochs and exact $1^{st}$ excited energy\")\n",
    "ax.set_ylabel(\"Energy (a.u.)\")\n",
    "ax.set_xlabel(\"Epoch\")\n",
    "fig.canvas.draw()\n",
    "print(\"Exact 1st excited ground energy: \", exact_energy)\n",
    "print(\"VQE 1st excited energy: \",losses[-1])\n",
    "print(\"Accuracy: \", 100*np.abs((losses[-1]-exact_energy)/exact_energy), \"% \\n\")\n",
    "plt.savefig(\"Project_QC/Ising_model/shot-based_backend/with_noise/1st_exc_%iqub_%ilay_ryrz_%s_%sNoise.pdf\"  % (num_qub, num_lay, entanglement, noise))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "caring-dining",
   "metadata": {},
   "source": [
    "### $2^{nd}$ excited energy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "adapted-bunny",
   "metadata": {},
   "outputs": [],
   "source": [
    "num_lay = 4\n",
    "entanglement = 'linear'\n",
    "ansatz_gates = ['ry', 'rz']\n",
    "noise = 'yes'\n",
    "ansatz = EfficientSU2(num_qubits = num_qub, su2_gates = ansatz_gates, entanglement = entanglement, reps = num_lay)\n",
    "parameters = np.loadtxt('Project_QC/Ising_model/shot-based_backend/with_noise/1st_exc_%iqub_%ilay_ryrz_%s_%sNoise_fin_params'  % (num_qub, num_lay, entanglement, noise))\n",
    "first_exc_circ = ansatz.bind_parameters(parameters)\n",
    "first_exc_state = CircuitStateFn(first_exc_circ)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "living-funds",
   "metadata": {},
   "outputs": [],
   "source": [
    "class sec_exc_VQE(VQE):\n",
    "    \n",
    "    def expectation_value(self, parameters):\n",
    "        alpha = 5\n",
    "        ansatz_dum = self.ansatz.bind_parameters(parameters)\n",
    "\n",
    "        ansatz_state = CircuitStateFn(ansatz_dum)\n",
    "        measurable_expression = StateFn(self.H, is_measurement=True).compose(ansatz_state) \n",
    "        overlap = alpha * (np.abs((~ground_state @ ansatz_state).eval()) + np.abs((~first_exc_state @ ansatz_state).eval()))\n",
    "        expectation = self.expectationer.convert(measurable_expression)  \n",
    "\n",
    "        sampler = CircuitSampler(self.quantum_instance).convert(expectation)\n",
    "        result = sampler.eval().real + overlap\n",
    "        return result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "celtic-building",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "%matplotlib notebook\n",
    "plt.ion()\n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(111)\n",
    "fig.show()\n",
    "fig.canvas.draw()\n",
    "\n",
    "num_lay = 6\n",
    "entanglement = 'linear'\n",
    "ansatz_gates = ['ry', 'rz']\n",
    "ansatz = EfficientSU2(num_qubits = num_qub, su2_gates = ansatz_gates, entanglement = entanglement, reps = num_lay)\n",
    "parameters = np.pi * np.random.normal(size = len(ansatz.parameters))\n",
    "\n",
    "noise = 'yes'\n",
    "vqe = sec_exc_VQE(ansatz, H, noise = noise)\n",
    "epochs = 500\n",
    "losses = []\n",
    "exact_energy = np.sort(np.real(np.linalg.eig(H.to_matrix())[0]))[2]\n",
    "\n",
    "print(\"Exact 2nd excited energy: \", exact_energy)\n",
    "\n",
    "for i in range(epochs):\n",
    "    \n",
    "    loss, parameters = vqe.update(parameters)\n",
    "    losses.append(loss)\n",
    "    ax.clear()\n",
    "    ax.plot(losses, 'b.', label = 'VQE energy')\n",
    "    ax.plot(exact_energy*np.ones_like(range(epochs)),'k--', label = r'Exact $2^{nd}$ excited energy')\n",
    "    ax.legend()\n",
    "    ax.set_title(r\"VQE energy vs epochs and exact $2^{nd}$ excited energy\")\n",
    "    ax.set_ylabel(\"Energy (a.u.)\")\n",
    "    ax.set_xlabel(\"Epoch\")\n",
    "    fig.canvas.draw()\n",
    "\n",
    "print(\"2nd excited energy from VQE: \", loss)\n",
    "print(\"Percentual accuracy: \", 100*np.abs((loss-exact_energy)/exact_energy), \" %\")\n",
    "np.savetxt('Project_QC/Ising_model/shot-based_backend/with_noise/2nd_exc_%iqub_%ilay_ryrz_%s_%sNoise_losses' % (num_qub, num_lay, entanglement, noise), losses)\n",
    "np.savetxt('Project_QC/Ising_model/shot-based_backend/with_noise/2nd_exc_%iqub_%ilay_ryrz_%s_%sNoise_fin_params'  % (num_qub, num_lay, entanglement, noise), parameters)\n",
    "plt.savefig(\"Project_QC/Ising_model/shot-based_backend/with_noise/2nd_exc_%iqub_%ilay_ryrz_%s_%sNoise.pdf\"  % (num_qub, num_lay, entanglement, noise))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "thick-maryland",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(np.mean(100*np.abs((losses[-30:]-exact_energy)/exact_energy)))\n",
    "print(losses[-30:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "sonic-office",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(111)\n",
    "num_lay = 6\n",
    "entanglement = 'linear'\n",
    "ansatz_gates = ['ry', 'rz']\n",
    "noise = 'yes'\n",
    "losses = np.loadtxt('Project_QC/Ising_model/shot-based_backend/with_noise/2nd_exc_%iqub_%ilay_ryrz_%s_%sNoise_losses' % (num_qub, num_lay, entanglement, noise))\n",
    "exact_energy = np.sort(np.real(np.linalg.eig(H.to_matrix())[0]))[2]\n",
    "ax.plot(losses, 'b.', label = 'VQE energy')\n",
    "ax.plot(exact_energy*np.ones_like(range(len(losses))),'k--', label = r'Exact $2^{nd}$ excited energy')\n",
    "ax.legend()\n",
    "ax.set_title(r\"VQE energy vs epochs and exact $2^{nd}$ excited energy\")\n",
    "ax.set_ylabel(\"Energy (a.u.)\")\n",
    "ax.set_xlabel(\"Epoch\")\n",
    "fig.canvas.draw()\n",
    "print(\"Exact 2nd excited ground energy: \", exact_energy)\n",
    "print(\"VQE 2nd excited energy: \",losses[-1])\n",
    "print(\"Accuracy: \", 100*np.abs((losses[-1]-exact_energy)/exact_energy), \"% \\n\")\n",
    "plt.savefig(\"Project_QC/Ising_model/shot-based_backend/with_noise/2nd_exc_%iqub_%ilay_ryrz_%s_%sNoise.pdf\"  % (num_qub, num_lay, entanglement, noise))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "behind-clone",
   "metadata": {},
   "source": [
    "## Three lowest energies and states"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "accessory-lying",
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib notebook\n",
    "plt.ion()\n",
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(111)\n",
    "fig.show()\n",
    "fig.canvas.draw()\n",
    "state_labels = ['ground', '1st_exc', '2nd_exc']\n",
    "nums_lay = [4, 4, 6]\n",
    "colors = ['blue', 'green', 'red']\n",
    "noise = 'no'\n",
    "for i in range(len(state_labels)): \n",
    "    losses = np.loadtxt('Project_QC/Ising_model/shot-based_backend/without_noise/%s_%iqub_%ilay_ryrz_%s_%sNoise_losses' % (state_labels[i], num_qub, nums_lay[i], entanglement, noise))\n",
    "    exact_energy = np.sort(np.real(np.linalg.eig(H.to_matrix())[0]))[i]\n",
    "    print(\"Exact %s energy: \" % (state_labels[i]), exact_energy)\n",
    "    print(\"VQE %s energy: \" % (state_labels[i]),np.mean(losses[-50:]))\n",
    "    print(\"Accuracy: \", 100*np.abs((np.mean(losses[-50:])-exact_energy)/exact_energy), \"% \\n\")\n",
    "    ax.plot(losses, 'b.', label = 'VQE energy', color = colors[i])\n",
    "    ax.plot(exact_energy*np.ones_like(range(len(losses))),'k--', label = r'Exact %s energy' % (state_labels[i]))\n",
    "    ax.legend()\n",
    "    ax.set_xlim([-2, 500])\n",
    "    title = r\"VQE energies vs epochs and exact energies\"\n",
    "    if(noise == 'yes'):\n",
    "        title += ' with noise'\n",
    "    ax.set_title(title)\n",
    "    ax.set_ylabel(\"Energy (a.u.)\")\n",
    "    ax.set_xlabel(\"Epoch\")\n",
    "    ax.set_ylim(top=2)\n",
    "fig.canvas.draw()\n",
    "plt.savefig(\"Project_QC/Ising_model/shot-based_backend/without_noise/3_energies_%sNoise.pdf\" % noise)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "geological-scheduling",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Exact ground energy:  -2.2360679774997876\n",
      "VQE ground energy:  -2.2363799999999996\n",
      "Accuracy:  0.013954070419669765 % \n",
      "\n",
      "VQE ground energy:  -2.044988\n",
      "Accuracy:  8.545356376573114 % \n",
      "\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZQAAAEYCAYAAAB7twADAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABYaklEQVR4nO2deXxU1d3/39+ZLIQshYQlQtBADShbQgjEaIRYXBC1WvelLlWLVm1rbd1+trV9tNqqrdatSBeXPrbWp1YfraJWHgJEIorByKIsApaALBJo9m3y/f1xZ+ZOwkwyITeZMHPer9d9zcy9Z8753s89937v2UVVMRgMBoOht7gibYDBYDAYogPjUAwGg8HgCMahGAwGg8ERjEMxGAwGgyMYh2IwGAwGRzAOxWAwGAyOYByKwWAwGBzBOBSDwWAwOIJxKFGMiDwjIvdG2o5IICLrRKSki+PbROTk/rPIWbo7v4FGT/JiV9dGRCaIyGoRqRWR7zlr5eFJd3m5P/NK1DgUEXlLRP4ryP6zRWSXiMR5f18lImtEpMG7/0kR+UpA+G0i0igidQHb4/15Lobeo6qTVLUUnHEe3eWbgHCXisgqb775QkQWiUhxb9IORuD5xRi3AaWqmqqqj0bamGCISKKI/FFEPvc6vtUicnqk7HH6XuiKqHEowDPA5SIinfZfDjyvqm0i8kPgV8CtwFeA44Bs4G0RiQ/4z1mqmhKw3dRb43wOLZaIlnMON9+IyC3AI8B9wEjgSOBJ4Oz+tTiqOQpY19M/9XNejAO2A7Ox8stPgBdFJLsfbYgMqhoVG5AE/AeYFbBvKNAE5AJpQB1wYaf/pQB7gCu9v7cBJ/cg3XxgNVAL/A/wN+DegLhuBz4GmrEy2rFAKXAA68b4ekBcChwd8PsZX1wB8f3IG99/vGkNCjg+Dajw2vI34IXA/3eyexTwErAX2Ap8rwfpdPffzufclUa3Ai91su0x4JEgNn8LeC3g92bgxYDf24G8wOsI/BloBxq91/+2cM4xIM5w881XvOEu6EHeCWlDV/mkcz716r3Dq+8GYE5316kLe2712lMP/BHLMS7yxv0OMDRM+0Lmxe7sIsQ9CPwf4MG6p+uA8V3ZQZC8GCTO64HXgSeAL4GdwCl98Hz6GDivi+Nd6RXOc+FOYD2wH3iajvfrNrq+F4Lmn0M6T6eFi+QG/B74Q8Dv64CPvN/nAm0hMtWzWKWYkJk5RHoJwOfA94F44FyghY4O5SNgDJbDi8d6CP4/73+/5r2IE3qQcd733pDpwCfA9Z1s+YE3nfOBVoI4FKyS6YfAT73/GwdsAU4LI51w/ht4zt1pdATWw2uI93cc1oN6ehC7x2E9OFze/30O7Ag4th9wdb6Owa5pV+fYKVy4+SZkuC7yT1AbussngecETMBypKO8+7OBr3Z3nbqw5z0sJzLaex0qsB52iVgP9Lu7s48u8mI4dgW7XgHHSoFrvd+7s2MbAXkxRHxPAtXAaV7b7gbeCRLun1h5L9j2z26u80gsJ3hMN8+RoPcu4T0X1nrPMx14N8jxoPcCIfLPoT6Do6nKC6wb/AIRSfL+vsK7D2AY8KWqtgX53xfA8IDfr4jIgYDt2yHSOw7rAfioqraq6j+wHhCBPKqq21W10Rs+Bfilqrao6v9hZdRLenCOj6rqTlWtBl4D8gJsicd6s29V1b8DH4SIYwYwXFX/y2vHFixnfHEY6YT738BzDqmRqn4BLAMu8O6ai3WdPuxstDetWq8ts4G3gB0icoz393JVbQ8t3UGEOsdAws03GV2E66kNPcknHqyH/UQRiVfVbar6GeFdp2A8pqq7VXUHsBxYqaqrVbUZeBnLuXRnX1d58VDtCkY4OgXmxWBM9f7/LW/eWR8skKqeqapDQmxnhjLQWyX6PPCsqn7axXmEe++G4nHveVYDvyD8Z0qo/HNIREUdtw9VLRORvcDZIvI+VuY913v4S2CYiMQFuemPwCp++zhHVd8JI8lRWG/IgWsAbO8UZnun8Ns7PfQ+x3obDJddAd8bvHGGsuXzEHEcBYwSkQMB+9xYD5Du0gnnv53PuTuNngW+g/Vg+SZW0TwUS4ES4Gjv9wNYzqTI+7snhDrHQMLNN/u6CNdTG8LOJ6q6WURuBn4GTBKRt4BbCO86BWN3wPfGIL9TwrCvq7x4qHYFIxydOue1zkwBAl8YJxPCqfQUEXFh5eUWoKt22J7cu6EIPM/PCZ6XDyJU/lHVnT1MH4iuRnkfz2GVTC4H3lZV3w1RjlWPem5gYBFJBk6n5w8jsN5QR3fqCDCmU5jATLITGOPNaD6OxKq/BOuBMjjgWGYvbTkyRNjtwNZOb1mpqjovjHTC+W/gOYej0SvAVBGZDJyJ9UYXCp9DOdH7fSmWQ5lN6GuoIfaHQ7j5phyrWuOcXqTlo7t80gFV/YuqFmM9rBWrA0FvrnFv7esqLzppVzg6hbz2IjIW66V6Q8DuaVjVZJ3DLurU8zNwWxQkvGC3QZ2nqq1dnEd39244z4XAe+pILG2CcZAeIfLPIRGtDuVkrLcOX3UXqvof4OfAYyIyV0Tivb0u/gfrLbSrh1goyrGKjDeJSJyInA3M7CL8Sqz2gtu86ZcAZ2E1wIGVkS8VEbeIzMV6SPbEljbge15bzu3ClveBGhG5XUSSvOlNFpEZYaTT0/92q5GqNgF/B/4CvK+q/+4i/aXASVh14lVYb7ZzsaqcVof4z26suvoeE26+8Yb7KfCEiJwjIoO9YU8XkQd6mGx3+cSPWOMyviYiiVgOrRFL795c497a11VedNKusHUKwVRgTacSzjSgsnNAVT1dO/b8DNyCdQn+HVaHgbO6qG7z0d29+xHdPxduFJEsEUnHalP6W4i0OtwLXeSfQyLqHIqqbgNWAMnAq52OPYAl9kNYdfFbsTz/yapaHxD0tU5vIC+HSKsF6831Gqyql29i1eE2dxH+61hvtl9iNQheEVC3+n2sG+IAcBnWm3u45+2z5SqsxumLgH+ECOvxppOHpcGXwB+weip1l06P/tsDjZ7Fqn7oqroLVd2I1UNlufd3DVaj7rte24JxP/Bjb3vYj7o8weBphpVvVPU3WNVNP8aqCtuOVdXxSg/T6y6fBJII/NIbbhcwAvh/vbnGvbWvq7zopF091CkYUwkojYhIBtbb/9qe2hKIiByF1SEoD9gV8By5LFj4MO7dcJ4LfwHexroXtmB1gAhG53shaP4J70wPRjpW28UWInI11tvnCd28FfckzpXAAlV92on4opFgGonIkcCnQKbXSQxY+iLfGAyBiMgzQJWq/jjStvSEqGqU7ymq+icRaQWOBw7pwSAis7HqYL/EenuYCrzpmJFRQHcaeevAbwFeGOjOBJzJNwZDNBIxhyIiY7DaOzKxBtssVNXfdgojwG+BeVgNU1epaoWTdqhql1UsYTABeBGr98tnwPlqdYU12ITUyNu4vRurZ8rciFnYQxzINwZD1BGxKi8ROQI4QlUrRCQVa7DTOaq6PiDMPOC7WA6lEPitqhZGxGCDwWAwdEnEGuVV9QtfaUNVa7FGCXfuZ3828JxavAcM8Toig8FgMAwwBkQbircb5jSsboCBjKbjgJ0q776DqpREZD4wH2Dw4MHTR48ejYiQmJiI2+2moaEBgLi4OJKSkqitrfX9j5SUFBoaGvB4rE5CycnJtLa20tLSAsCgQYMQERobrd5/8fHxJCYmUldXB4DL5SI5OdmROOrr62lvt3oxpqSk0NzcTGur1YU9KSkJVaWpqQmAhIQE4uPjqa+3Ohq53W4GDx58UBwNDQ3+3+HGUVdX55uagdTUVBobG2lrs8brDR48GI/HQ3Oz1VGrO42diMOJ6+SzY6Bep0OJ41Cvk6oSHx8/IK9Tf99PjY2NiMiAvE7BNO6r67Rx48YvVTVwxpAeE3GHIiIpWBPF3RykQbbzzMEQYqCSqi4EFgIUFBToQw89RElJiZOmHraUlpYaLTA6BGK0sDFaWIhIT0fnH0REx6GINc/NS1gT7AUbM1FFxxGgWYQeAdqB3Nzc3hsYJRgtLIwONkYLG6OFc0TMoQRMTfCJd0BYMF4FrhCL44D/hNuDKrB4GOsYLSyMDjZGCxujhXNEsoRyAtZ8W18TkY+82zwRuV5ErveGeQNr1OdmrIkDbwg38i1btjhu8OGK0cLC6GBjtLAxWjhHxNpQVLWM4G0kgWEUuLF/LDIY+obW1laqqqr8DbcDga985St88sknkTZjQBBrWgwaNIisrCzi4+O7D9xDIt4o31dkZ2dH2oQBg9HCIlI6VFVVkZqaSnZ2NnLQCtWRobm5mcTExEibMSCIJS1UlX379lFVVcXYsWMdjz/qJof0kZ6eHmkTBgxGC4tI6dDU1ERGRsaAcSZgdSM1WMSSFiJCRkZGn5WWo9ahVFQ4OkPLYY3RwiKSOgwkZwL4xycYYk+LvsyLUetQDAaDwdC/RK1DGTp0aKRNGDAYLSxiVYeSkhLeeuutDvt+97vfccMNVqfJdevW8bWvfY3x48fz1a9+lbvvvts/MvyZZ55h+PDh5OXl+bf16x1ZIXdAMG/ePGprazlw4ABPPvmkf39paSlnnhlyqfgeUVpayooVKxyJC+D44493LC6niVqHYgYr2RgtLGJVh0suuYQXXui4iOE//vEPLrnkEhobG/n617/OHXfcwcaNG1mzZg3vv/8+v/2tPfH3RRddxEcffeTfJk6c2C92+6YV6UveeOMNRo0adZBDcRKnHYqTcTlN1DqUpUsPZYn46MRoYXE46VBeDvffb332lvPPP59//vOf/nmetm3bxo4dOyguLuYvf/kLJ5xwAqeeeipgzQn1+OOP8+CDD/Yojf/+7/9m5syZ5OXlcd111/nnh0pJSeGuu+4iNzeX4447jt27dwOwd+9ezjvvPGbMmMGMGTN49913AfjZz37G/PnzOfXUU7niiivYu3cvp5xyCvn5+Vx33XUcddRRfPnll/zkJz/p4PTuuusuHn300Q42PfDAA/59P/jBD/ja174GwOLFi/nmN78JWD3/tm3bxh133MFnn31GXl4et956KwB1dXWcf/75HHPMMVx22WX+ubgWL17MtGnTmDJlCldffbVf1+zsbL788ksAVq1aRUlJCdu2bWPBggU8/PDD5OXlsXz58g42/uxnP+Pqq6+mpKSEcePGdTiH3/zmN0yePJnJkyfzyCOP+PenpKQA8MUXXzBr1izy8vKYPHmyP+63336boqIi8vPzueCCCzrMZdfnqGrUbdOnT9clS5aowcJoYREpHdavX9+j8CtWqCYlqbrd1ueKFb23Yd68efrKK6+oqur999+v3/ve91RV9Qc/+IE+8sgjB4UfMmSI7t+/X59++mkdNmyY5ubm+reGhoYOYdevX69nnnmmtrS0qKrqd77zHX322WdVVRXQV199VVVVb731Vr3nnntUVfWSSy7R5cuXq6rq559/rsccc4yqqt59992an5/vT+PGG2/U++67T1VVFy1apIDu3btXt27dqtOmTVNVVY/Ho+PGjdMvv/yyg13l5eV6/vnnq6pqcXGxzpgxQ1taWvRnP/uZLliwQFVVjzrqKN26datu3bpVJ02a5P/vkiVLNC0tTbdv364ej0ePO+44Xb58uTY2NmpWVpZu2LBBVVUvv/xyffjhh/1x7d27V1VVP/jgA509e7b/nB588MGg1+Xuu+/WoqIibWpq0r1792p6erq2tLToqlWrdPLkyVpXV6e1tbU6ceJEraioUFXV5ORkVVV96KGH9N5771VV1ba2Nq2pqdG9e/fqiSeeqHV1daqq+stf/lJ//vOfH5RusDwJrNJePnujtr/cQOtVE0mMFhaHiw6lpdDSAh6P9VlaCkVFvYvTV+119tln88ILL/DYY48B1gtlMF00YJ2kiy66iMcffzxk3IsXL+bDDz9kxowZADQ2NjJixAjAmn3X1xYxffp0/vWvfwHwzjvvdGiLqamp8U+B8vWvf52kpCQAysrKePnllwGYO3euvx0sOzubjIwMVq9eze7du5k2bRoZGRkd7Jo+fToffvghtbW1JCYmkp+fz6pVq1i+fPlBpZlgzJw5k6ysLADy8vLYtm0bqampjB07lvHjxwNw5ZVX8sQTT3DzzTd3G18ozjjjDBITE0lMTGTEiBHs3r2bsrIyvvGNb5CcnAzAueeey/Lly5k2bZr/fzNmzODqq6+mtbWVc845h7y8PJYuXcr69es54YQTAGhpaaGot5mnB0StQ5k9e3akTRgwGC0sDhcdSkogIcFyJgkJ1u/ecs4553DLLbdQUVFBY2MjJ554IgCTJk1i2bJlHcJu2bKFYcOGMWTIkLDiVlWuvPJK7r///oOOxcfH+x2W2+32t4u0t7dTXl7udxyB+B6ivrhDce211/LMM8+wa9curr766qBpZ2dn8/TTT3P88cczdepUlixZwmeffcaxxx7rD5eSkhK0WihwsKPP9q7siYuL83dm6Mk4j56m42PWrFksW7aM119/ncsvv5xbb72VoUOHcsopp/DXv/417PSdJGrbUCorKyNtwoDBaGFxuOhQVASLF8M991ifTrxgpqSkUFJSwtVXX80ll1ziH3tx2WWXUVZWxjvvvANYpYvvfe97/PznPw877jlz5vD3v/+dPXv2AFBdXc3nn3c9E/qpp57aodTz0UcfBQ1XXFzMiy++CFhtA/v37/cf+8Y3vsGbb77JBx98wGmnnRb0/7NmzeKhhx5i1qxZnHjiiSxYsIC8vLwOpbKGhgZSU1PDmiTymGOOYdu2bWzevBmAP//5z/4XlezsbD788EMAXnrpJf9/wo27s92vvPIKDQ0N1NfX8/LLL/tfAnx8/vnnjBgxgm9/+9tcc801VFRUcNxxx/Huu+/67WtoaGDjxo09Srs3RK1DCcx4sY7RwuJw0qGoCO680xln4uOSSy6hsrKSiy++2N9onpSUxKuvvsovfvELxo8fz7BhwzjhhBO47LLL/P/729/+1qHbcOdeRhMnTuTee+/l1FNPZerUqZxyyil88UXXk4I/+uijrFq1iqlTpzJx4kQWLFgQNNzdd9/N22+/TX5+PosWLeKII44gNTUVsKrTTjrpJC688ELcbnfQ/5944ol88cUXFBUVMXLkSAYNGnTQg7m9vZ2MjAxOOOEEJk+e7G+UD8agQYN4+umnueCCC5gyZQoul4vrr7/eb+v3v/99TjzxxA72nHXWWbz88stBG+VDkZ+fz1VXXcXMmTMpLCzk2muv7VDdBVbvsby8PKZNm8ZLL73E97//fYYPH84zzzzDJZdcwtSpUznuuOP49NNPw0rTEXrbCDMQN9Mo3xGjhcXh0ijfH9TU1ATd//LLL+vYsWN127Zt/WxRcJqamrS1tVVVVVesWKG5ubn+Yx6PR3Nzc3Xjxo29SiOUFtGMaZTvIfn5+ZE2YcBgtLAwOtgMHjw46P5zzjmHc845p3+N6YJ///vfXHjhhbS3t5OQkMDvf/97ANavX8+ZZ57JN77xDXJycnqVRigtDD0nah1KdXU1aWlpkTZjQGC0sDA62LS1tYWsJhpI5OTksHr16oP2T5w40bF1TA4XLQ4HorYNZdu2bZE2YcBgtLAwOti0tLRE2oQBg9HCOSK9pvyfRGSPiKwNcbxERP4TsKLjT/vbRoPBYDCER6SrvJ4BHgee6yLMclXt8Sxt48aNO1Sbog6jhYXRwSYhISHSJgwYjBbOEdESiqouA6r7Im5f10KD0cKH0cHGtBnYGC2cI9IllHAoEpFKYCfwI1VdFyyQiMwH5gOMGjWKd999l5SUFMaNG0dqaqp/UFtGRkaH0cFxcXEUFxdTUVFBTU0NAAUFBezevZvt27cDVsNgYmIia9daNXMjRoxg/PjxlJWVAdZI16KiIlatWuUfcVtYWEhVVRU7duwAYMKECbjdbv90E5mZmYwdO5Zy7+x/SUlJFBYWsnLlShobG60TLypi69at7Nq1C7AaIj0eDxs2bABg9OjRZGVlsXLlSsAavFZQUEB5ebl/wrri4mJWrlzJoEGDAJg8eTLNzc1s2rQJgDFjxjBy5EhWrVoFQFpaGvn5+ZSVlflHNc+aNYt169axb98+wJq1t7a21t8omp2dTXp6un8Bq6FDh5Kbm8vSpUv9U3vMnj2byspK/1iQ/Px8qqur/e0a/XGd6uvrOeOMM/r9OqWnp9Pe3k59fT0ALpeL5ORk6urq/COiU1JSaGpq8ms+aNAgVNV/HePj44mPj/cPSPTFEThgLiUlhcbGxg5jTDweDyeffDK33HIL8+bNIy4ujoaGBh5//HG2bt3KwoULef/997n11lvZsWMHqspFF13E7bffjsvl4n/+53+4/fbbOeKIIwBr+po///nP/uVj3W43SUlJHUaap6amUl9f7x81PnjwYFpbW2ltbQWs+0VE/KPJ4+LiGDRokD8OESElJaVDHMnJybS0tHQZR2Jiol9jXxyBGicnJ9Pc3NxB47POOounnnoKl8vFP/7xD7773e9SX1/P8uXLeeyxx1i0aJGj16miooK///3v/OpXvzroOvnachISEvzXyafx4MGDO1zr1NRUGhoa/HEMHjyYtra2DnG43W5/Hu18nXy6Bd5PjtDbfse93YBsYG2IY2lAivf7PGBTOHGacSgdMVpYxOo4lAULFuhVV13VYV9BQYEuW7ZMGxoadNy4cfrWW2+pqmp9fb3OnTtXf/Ob36iq6tNPP6033nhjv9usqv7xJ31NTU1N0MkhzzjjjH5JPxL01TiUAd3LS1VrVLXO+/0NIF5EhoXz384TxcUyRguLw0mH8vJy7r//fn/JqDcEm75+165dZvp6rNL1/v37ezR9fSAlJSXcfvvtzJw5k/Hjx/tHwjc1NfGtb32LKVOmMG3aNJYsWQJ0XLhr6dKl/tkHpk2b5i+BPPjgg8yYMYOpU6dy99139+g6RJzeeqTebnRdQskExPt9JvBv3++utunTp6vH4+mx145WjBYWkdKh59PXr9CkpCR1u92alJSkKxyYv77z9PU//OEPVdVMX3/UUUfpnj17wp6+vjOzZ8/WW265RVVVX3/9dZ0zZ46qWlPL+0qFn3zyiY4ZM0YbGxs7lHzOPPNMLSsrU1XV2tpabW1t1bfeeku//e1va3t7u3o8Hj3jjDN06dKlB6XbW6KyhCIifwXKgQkiUiUi14jI9SJyvTfI+cBabxvKo8DF3hPvls4zqMYyRguLw0WH0tJSWlpa/PXqpaWlvY4zcNVG3zT2EP709YErNnaeIThw+vq8vDwWL17sb1/rPH29r83snXfe4aabbiIvL4+vf/3rXU5ff/HFFwOhp69/++23w5q+3tfOuXz58g7zefnaXTrjm77e5XL5p68PxrnnnnvQ+ZWVlXH55ZcD1oSSRx111EGTNJ5wwgnccsstPProoxw4cIC4uDjefvtt//nk5+fz6aef+ts7Dwci2iivqpd0c/xxrG7FBkPMUFJSQkJCAi0tLSQkJFDiwPz1naevz8vLA8z09V0RbFr5rsIFhgnnvfeOO+7gjDPO4I033uC4447jnXfeQVW58847ue6667r9/0BkQLeh9Ia4uMOhA1v/YLSwOFx0KCoqYvHixdxzzz0sXrzYkQWSOk9f78NMX29xKFPMd8WsWbN4/vnnAdi4cSP//ve/mTBhQocwn332GVOmTOH222+noKCATz/9lNNOO40//elP/t5YO3bs8Ot6OBC1DqW4uDjSJgwYjBYWh5MORUVF3HnnnY6uthc4fb1vTI6Zvt5ytuFOXx8uN9xwAx6PhylTpnDRRRfxzDPPdCjxADzyyCNMnjyZ3NxckpKSOP300zn11FO59NJLKSoqYsqUKZx//vmOOro+p7eNMANxmz59un744YfhtE3FBEYLi0jpEOluw8HwrTnemVicvj6UFtGMmb6+hzg6WOcwx2hhYXSw8Q0Y7EwsTl8fSgtDz4lah2IwGA5/+mP6eoNzRG0bSkFBQaRNGDAYLSyMDjZmUSkbo4VzRK1D8Y3INRgtfBgdbHxzYhmMFk4StQ7FN2GgwWjhw+hgYx6iNkYL54hah2IwGAyG/iVqHUpve35EE0YLC6ODTecxEbGM0cI5otahmExiY7SwMDrYBJu/K1YxWjhH1DoU3yJLBqOFj1jXwe12+0e7T58+nV/+8peOxHvgwAGefPLJkMefeuopjjjiiA6j7desWeNI2k7gW2zK0HvMOBSDIUZISkryz5lVW1vr2JLIPodyww03BD3+8ccfc++993LNNdc4kp5h4BK1JZQRI0ZE2oQBg9HCwuhg45so84MPPmDq1Kk0NTVRX1/PpEmT/CW5c845h+nTpzNp0iQWLlzo/+9zzz3H1KlTyc3N5fLLLw+6OFUga9as8c9uPBA5XCYNPRzwLV4VVRQUFOh7771nMoqXtrY2owWR0+GTTz7pMF16sOnoL7zwQm644QYaGhqYN2/eQcevuuoqrrrqKr788kvOP//8DsfCXS/F7XYzZcoU/+8777yTiy66iB//+Mc0NTXR2NhIVlYWd955J2DNGpyenk5jYyMzZsxg6dKl7Nq1i3PPPZd3332XYcOGUV1dTU1NDWeeeWbIKsWMjAxGjx6Ny2W9v95www3Mnz8/LJv7Aw2xJkw00zlPAojIh6raq9G/UfuUKSsrc2QdiWjAaGER6zqEqvL66U9/yowZMxg0aFCHZXQfffRRXn75ZcAaw7Np0yY++OADzj//fIYNs1biTk9P73KOtO3btzNixAg+/vjjg47V19dzww03+Nd8CZzhuD+pq6tzrPov1omoQxGRPwFnAntUdXKQ4wL8FpgHNABXqWpF/1ppMDhLVyWKwYMHd3l82LBhjqzgGEh1dTV1dXW0trbS1NREcnIypaWlvPPOO5SXlzN48GBKSkpoamrq8dv8xx9/zDHHHBP02D/+8Q/OP/98zjrrLC666KKIORSDc0S6DeUZYG4Xx08HcrzbfOB34UZsuojaGC0sjA42gU5h/vz53HPPPVx22WXcfvvtAPznP/9h6NChDB48mE8//ZT33nsPsBbTevHFF9m3bx9gOaOuFqdas2ZNSIdSVVXFmDFjAEKuZ9IfxFp1V18S6SWAl4lIdhdBzgae887V/56IDBGRI1S169V7wNGFiQ53jBYWsa5D4NK/YK3RPnHiROLi4rj00kvxeDwcf/zx/N///R9z585lwYIFTJ06lQkTJnDccccB1pLBd911F7Nnz8btdjNt2jSeeeYZ/+JUp59+Og8++KA/jTVr1rB06VIWLVoEWA/v5cuXk5KSQlZWFlVVVeTl5UV0CvmUlJSIpR1tRLxR3utQ/hmiyuufwC9Vtcz7ezFwu6quChJ2PlYphlGjRk3//e9/z+DBgxk3bhypqalUVlYCVgNh4DracXFxFBcXU1FR4a8LLigoYPfu3f65n3JyckhMTPQ3Oo4YMYLx48dTVlYGWG++RUVFrFq1yr90Z2FhIVVVVezYsQOACRMm4Ha7Wb9+PQCZmZmMHTuW8vJywKrfLiwsZOXKlTQ2NgLWA3Dr1q3s2rULsKbs9ng8bNiwAYDRo0eTlZXFypUrAevGKCgooLy8nObmZsBapXDp0qX+N8DJkyfT3NzMpk2bABgzZgwjR45k1SpL0rS0NPLz8ykrK/Ovjz1r1izWrVvnfyvNzc2ltrbWP314dnY26enpVFRYtZFDhw4lNzeXpUuX+qtIZs+eTWVlpX8J1/z8fKqrq9m2bRtAv1ynpqYm5s6d2+/XKT09ncmTJ1NfXw+Ay+UiOTmZuro6/9rjKSkpNDU1+TUfNGgQquq/jvHx8cTHx9PQ0NAhjsCSQUpKCo2NjXg8Hr+tHo+HlpYWwFrhMC4ujoaGBtrb24mPj2fw4MEd4khNTaWhocEfx+DBg2lra+sQh9vt9p+72+0mKSnJr6cvjvr6er+TGDx4MK2trf45sxITExER9u3bx49+9COSkpKYPXs2Z599NmA5nZSUlA5xJCcn09LSclAcvjEkcXFxJCYm+jX2xRGocXJyMs3NzQdp3NjYiMvlIj4+noSEhAF1nXwaO32dNm/ezLRp0zrcTyeddFKvG+UHukN5Hbi/k0O5TVU/7CrOgoICfeihh2K6ATaQ0tJSowWR0yFYj5pI4+Q4lMOdWNSir3p5RboNpTuqgDEBv7OAnRGyxWAwGAxdMNAdyqvAFWJxHPCfcNpPwKrKMFgYLSyMDjbJycmRNmHAYLRwjog6FBH5K1AOTBCRKhG5RkSuF5HrvUHeALYAm4HfA8HndghCVVWV4/YerhgtLIwONr66doPRwkki3cvrkm6OK3DjocS9Y8cOM125F6OFRSR1GGijsVtbWxk0aFCkzRgQxJoWfdluPtCrvAyGw55Bgwaxb9++Pr2RDYZwUFX27dvXZw40aqdemTBhQqRNGDAYLSwipYNvvMXevXsjkn4wPB5PRAcTDiRiTYtBgwaRlZXVJ3FHrUOJpQzSHUYLi0jpEB8fz9ixYyOSdij27NljZl/2YrRwjqit8vINTDMYLXwYHWyMFjZGC+eIWodiMBgMhv4lah1KZmZmpE0YMBgtLIwONkYLG6OFc0StQxloddaRxGhhYXSwMVrYGC2cI2odim8yP4PRwofRwcZoYWO0cI6odSgGg8Fg6F+i1qEkJSVF2oQBg9HCwuhgY7SwMVo4R8Snr+8LCgoK1Le+h8FgMBi6Jxamrz9kfItOGYwWPowONkYLG6OFc0StQ/GtVGYwWvgwOtgYLWyMFs4RtQ7FYDAYDP1L1LahvPvuuyQmJkbalAFBc3Oz0QKjQyBGCxujhcVh34YiInNFZIOIbBaRO4IcLxGR/4jIR97tp+HGvXXrVmeNPYwxWlgYHWyMFjZGC+eImEMRETfwBHA6MBG4REQmBgm6XFXzvNt/hRv/rl27HLL08MdoYWF0sDFa2BgtnCOSJZSZwGZV3aKqLcALwNkRtMdgMBgMvSCSDmU0sD3gd5V3X2eKRKRSRBaJyKRwI584MVhhJzYxWlgYHWyMFjZGC+eI5AJbwRbY7txDoAI4SlXrRGQe8AoQdFFwEZkPzAcYNWoUlZWVxMfHM27cOFJTU6msrAQgIyODSZMmsWzZMgDi4uIoLi6moqKCmpoaAAoKCti9ezfbt1v+Licnh8TERNauXQvAiBEjGD9+PGVlZQAkJiZSVFTEqlWrqKurA6CwsJCqqip27NgBWKsFut1u/9oLmZmZjB071j+PUFJSEoWFhaxcudLfjbGoqIitW7f6i+QTJ07E4/GwYcMGAEaPHk1WVpa/H31KSgoFBQWUl5fT3NwMQHFxMdu2bfOnO3nyZJqbm9m0aRMAY8aMYeTIkfgGgqalpZGfn09ZWRltbW0AzJo1i3Xr1rFv3z4AcnNzqa2tZcuWLQBkZ2eTnp5ORUUFAEOHDiU3N5elS5f611KfPXs2lZWV7N+/H4D8/Hyqq6vZtm0bQL9cJ9++gXqdNm7cyJ49e/rlOrW2tjJixIgBeZ36+35au3Yt8fHxA/I6Qf/dT04QsV5eIlIE/ExVT/P+vhNAVe/v4j/bgAJV/bKruAsKCvShhx6ipKTEOYMPY0pLS40WGB0CMVrYGC0sDvdeXh8AOSIyVkQSgIuBVwMDiEimiIj3+0wse/f1u6UGg8Fg6Jawq7xEJBloUlWPEwmrapuI3AS8BbiBP6nqOhG53nt8AXA+8B0RaQMagYs1zCLV6NHBmmNiE6OFhdHBxmhhY7RwjpBVXiLiwio1XAbMAJqBRGAv8AawUFU39ZOdPaKgoECXL19uZhH10tjYaLTA6BCI0cLGaGHR11VeS4CvAncCmao6RlVHACcC7wG/FJFv9ibxvsRM+GZjtLAwOtgYLWyMFs7RVZXXyara2nmnqlYDLwEviUh8n1lmMBgMhsOKkCWUYM7kUMJEipSUlEibMGAwWlgYHWyMFjZGC+c4pG7DIvJPVT2zD+xxBLPAlsFgMPSMSHYb/nZvEu0PfAOcDEYLH0YHG6OFjdHCOQ7JoajqF04b4jS+ka0Go4UPo4ON0cLGaOEc3Y5DEZGtHDwlCqo6rk8sMhgMBsNhSbdtKCKSEfBzEHABkK6qYa9N0t8UFBToe++9R1xcJKcqGzi0tbUZLTA6BGK0sDFaWPRLG4qq7gvYdqjqI8DXepNof7Bx48ZImzBgMFpYGB1sjBY2Rgvn6NahiEh+wFbgnRoltR9s6xW+2UANRgsfRgcbo4WN0cI5winn/TrgexuwFbiwb8wxGAwGw+FKtw5FVU/qD0OcZvLkyZE2YcBgtLAwOtgYLWyMFs5xSN2GRSTfaUOcxnQFtDFaWBgdbIwWNkYL5zjUgY3fcdSKPsC3gprBaOHD6GBjtLAxWjjHoQ5sHPAj5Q0Gg8HQv4TV+VpEhmKt5T7It09Vl/WVUU4wZsyYSJswYDBaWBgdbIwWNkYL5win2/C1wDKslRV/7v38mROJi8hcEdkgIptF5I4gx0VEHvUe/7gnbTcjR450wsSowGhhYXSwMVrYGC2cI5wqr+9jrdj4ubfH1zSsVRt7hYi4gSeA04GJwCUiMrFTsNOxSkY5wHzgd+HGb2YbtjFaWBgdbIwWNkYL5wjHoTSpahOAiCSq6qfABAfSnglsVtUtqtoCvACc3SnM2cBzavEeMEREjnAgbYPBYDA4TDhtKFUiMgR4BfiXiOwHdjqQ9mhge2A6QGEYYUYDB812LCLzsUoxjBo1ioaGBkpLSxk3bhypqalUVlYCkJGRwaRJk1i2zGoCiouLo7i4mIqKCmpqagAoKChg9+7dbN9uJZ2Tk0NiYiJr164FYMSIEYwfP56ysjIAEhMTKSoqYtWqVdTV1QFQWFhIVVUVO3bsAGDChAm43W7Wr18PQGZmJmPHjvVPnZ2UlERhYSErV66ksbERgKKiIrZu3cquXbsAmDhxIh6Phw0bNljijB5NVlaWfwnTlJQUCgoKKC8v93eFLC4upr29ndLSUsDqc9/c3Ozv2TJmzBhGjhzpf0tLS0sjPz+fsrIy2traAJg1axbr1q1j3759AOTm5lJbW8uWLVsAyM7OJj09nYqKCgCGDh1Kbm4uS5cuRVUREWbPnk1lZSX79+8HID8/n+rqarZt2wbQL9fJp8lAvU4bN270j9ru6+vU0NBAZWXlgLxO/X0/+Z4VA/E6Qf/dT07QowW2RGQ28BXgTW+p4tATFrkAOE1Vr/X+vhyYqarfDQjzOnC/qpZ5fy8GblPVD7uK+9hjj9UrrriCkpISioqKemOmwWAwxAR9OjmkiBy0LqaqLlXVV33OJFiYHlAFBHavyOLgkk84YQ5iw4YN/PjHP2bOnDlm8Rzwv/nFOkYHG6OFjdHCObpqQ/lfEfm1iMwSkWTfThEZJyLXiMhbwNxepP0BkCMiY0UkAbgYeLVTmFeBK7y9vY4D/hPO4l6qSnt7Oy0tLf6qnljGV9SOdYwONkYLG6OFc4RsQ1HVOSIyD7gOOME7FqUN2AC8DlypqrsONWFVbRORm7C6IbuBP6nqOu9sxqjqAuANYB6wGWgAvhVO3CKCiJCQkEBJScmhmmgwGAyGHtCjNpTDhWOPPVYvv/xyTjrpJNOGArS3t+NyHeosO9GD0cHGaGFjtLDolwW2Dk+S2b37asA4E4B169ZF2oQBgdHBxmhhY7Rwjqh0KBs2wOOPj2TOHDBt8vi7J8Y6Rgcbo4WN0cI5otKhqEJ7u9DSAqZN3mAwGPqHcObyekhEJvWHMU4hAi7XCkTuJyPDFFFyc3MjbcKAwOhgY7SwMVo4RzgllE+BhSKyUkSuF5Gv9LVRvWXMmHpcrpNpb/8JN99sxqLU1tZG2oQBgdHBxmhhY7Rwjm4diqr+QVVPAK4AsoGPReQvIjJglwb2eGppb2+mvd1jxqKAf0qHWMfoYGO0sDFaOEdYbSjemYGP8W5fApXALSLyQh/adsikpqYSHx+P2+02Y1EMBoOhn+h2ckgR+Q3wdWAxcJ+qvu899CsR2dCXxh0qycnJ/OUvf2HDhg1mPi+syeYMRodAjBY2RgvnCGe24bXAj1W1IcixmQ7b4xgnn3wy5557bqTNGBCkp6dH2oQBgdHBxmhhY7RwjnCqvD4CjhGR/IDtqyISp6r/6WP7Dhnf1M8Go4UPo4ON0cLGaOEc4ZRQngTygY8BASZ7v2eIyPWq+nYf2nfIrFuXRnk5lJRAjNd4GQwGQ78QjkPZBlyjqusAvMv03grcA/wDGHAOpb4efvjDPNraICEBFi+ObacydOjQSJswIDA62BgtbIwWzhFOldcxPmcCoKrrgWmqOmD72tXWQlubC48HM1oeM3DLh9HBxmhhY7RwjnAcykYR+Z2IzPZuT3r3JQKtfWzfIZGaCi5XGSL343aXE+u9hpcuXRppEwYERgcbo4WN0cI5wqnyuhK4AbgZqw2lDPgRljMZoIMb64GTgTastbsWE8szD0fjEgWHgtHBxmhhY7Rwji4dindA42uqejLw6yBB6vrEql5SW1uLx9OKajttbdZI+VgeiyIikTZhQGB0sDFa2BgtnKPbBbZE5FXgcie7CItIOvA3rKlctgEXqur+IOG2AbWAB2gLd/GXY489Vj///HNaWlpISEhg8eLFMe1QDAaDoTv6a4GtJmCNiPxRRB71bb1JFLgDWKyqOVj1UXd0EfYkVc3ryYkmJyezYMEC7rnnHuNMgMrKykibMCAwOtgYLWyMFs4RThvK697NSc4GSrzfnwVKgdudTODII4/kiiuucDLKw5b9+w8q/MUkRgcbo4WN0cI5unUoqvqsiCQBR6qqU3N3jVTVL7zxfyEiI0IlD7wtIgo8paoLQ0UoIvOB+QCjRo2irq6O0tJSxo0bR2pqqv8tJCMjg0mTJrFs2TIA4uLiKC4upqKigpqaGgAKCgrYvXs327dvByAnJ4fExETWrl0LwIgRIxg/fjxlZWUAJCYmUlRUxKpVq6irs5qVCgsLqaqqYseOHQBMmDABt9vN+vXrAcjMzGTs2LH+qfWTkpIoLCxk5cqVNDY2AlBUVMTWrVvZtWsXABMnTsTj8bBhg3UZRo8eTVZWFitXrgQgJSWFgoICysvLaW5uBqC4uJimpib/jMuTJ0+mubmZTZs2ATBmzBhGjhzJqlWrAEhLSyM/P5+ysjLa2toAmDVrFuvWrfOvbJebm0ttba1/ltbs7GzS09P9I46HDh1Kbm4uS5cuRVUREWbPnk1lZaX/5s3Pz6e6uppt27YB9Mt1qq+vBxiw12njxo3s2bOnX65TXV0dlZWVA/I69ff95HtWDMTrBP13PzlBOG0oZwEPAQmqOlZE8oD/UtWvd/O/d4DMIIfuAp5V1SEBYfer6kGji0RklKru9DqcfwHfVdVl3ZwTBQUF+stfLuWDD5LNSHmgpqaGtLS0SJsRcYwONkYLG6OFhRNtKOFUef0MaxLIUgBV/UhExnb3J2/PsKCIyG4ROcJbOjkC2BMijp3ezz0i8rLXjm4dSn09nHVWEi0t5bhcpTzxRAnz58euV6murjY3DEaHQIwWNkYL5winUb4tSA+v3nbcfhVrfAvez//tHEBEkkUk1fcdOBVr5uNuqa2F5ub3aG+fQ1vbT7jppthetdFXDI51jA42Rgsbo4VzhONQ1orIpYBbRHJE5DFgRS/T/SVwiohsAk7x/kZERonIG94wI4EyEakE3gdeV9U3w4ncGim/BGgBPHg8ZtVGg8Fg6GvCcSjfBSYBzcBfgRqsUfOHjKruU9U5qprj/az27t+pqvO837eoaq53m6Sqvwg3/uRkuPfeXOLjE3C53CQmxvaqjePGjYu0CQMCo4ON0cLGaOEc4fTyasBqSL+r781xjuuuO4HZsxdTWloa86s2pqamRtqEAYHRwcZoYWO0cI5uSygiMl5EForI2yLyf76tP4zrDVZ3uSLgTmJ5Hi8wA7d8GB1sjBY2RgvnCKeX1/8AC4A/YE2Bcliwbl0at94Kzc2mp5fBYDD0B+E4lDZV/V2fW+IwGzeOorm5nPb2ObS3t3DTTQlMmRKb07BkZGRE2oQBgdHBxmhhY7RwjnAa5V8TkRtE5AgRSfdtfW5ZL7nwwhG4XKWYnl4wadKkSJswIDA62BgtbIwWzhGOQ7kSa8nfFcCH3m1VXxrlBK2ty3jiiRLT0wv80y3EOkYHG6OFjdHCOcLp5dXtqPiByvz5RUyZYnp6GQwGQ38QsoQiIrcFfL+g07H7+tIoJ4iLC6d5KDYwWlgYHWyMFjZGC+cIOTmkiFSoan7n78F+DzQKCgp01apVlJeXM2fOHJqbm3G5XDzxxBPMnz8/0uYZDAbDgKOvF9iSEN+D/R5wVFRUUFpaSnNzM+3t7bS1tXHTTTfF5JxevmmwYx2jg43RwsZo4RxdORQN8T3Y7wFHTU0NJSUliNin6PF4YrKnl5PrHRzOGB1sjBY2Rgvn6KryMFdEarBKI0ne73h/D+pzyxyhCJEngJsAD/HxiTHb08tgMBj6mpAORVXd/WmI0xQUFPDYY6A6H5iCSCnf+lZs9vQqKOhVtWjUYHSwMVrYGC2cI5xxKIclu3fvpqQEEhLA7S4iPr4EKI3JNpTdu3dH2oQBgdHBxmhhY7Rwjqh1KNu3b6eoCB55BKZPL0d1Dr///U+YMyf2FtvyreUd6xgdbIwWNkYL54hahwJQXg433wwffFBKa2sLHo+HlpbYnYLFYDAY+pKIOBQRuUBE1olIu4iErMAUkbkiskFENovIHT1JIycnh9JSaGkB1RIgARE3CQmxNwVLTk5OpE0YEBgdbIwWNkYL54hUCWUtcC4QchIdEXEDTwCnAxOBS0RkYrgJJCYmUlICbjdY66E8Aszhu999JOYa5hMTEyNtwoDA6GBjtLAxWjhHRByKqn6iqhu6CTYT2OxdCrgFeAE4O9w01q5dS1ERXH01QDlwM6rv8NBDN7Jw4cJDtv1wZO3atZE2YUBgdLAxWtgYLZxjIE9iMxoIbC2rAgpDBRaR+cB8gFGjRlFXV0dpaSmnnnoMf/jDEtramoF22tvbuemmm2htbWXSpEnExcVRXFxMRUWFf4BTQUEBu3fv9jfW5eTkkJiY6M94I0aMYPz48ZSVlQHWG05RURGrVq2irq4OgMLCQqqqqtixYwcAEyZMwO12s379egAyMzMZO3asv4NAUlIShYWFrFy5ksbGRgCKiorYunUru3btAmDixIl4PB42bLB88ejRo8nKymLlypUApKSkUFBQQHl5Oc3NzQAUFxfT1NTkbzeaPHkyzc3NbNq0CYAxY8YwcuRIVq2yJpBOS0sjPz+fsrIy2traAJg1axbr1q1j3759AOTm5lJbW8uWLVsAyM7OJj093T/ieOjQoeTm5rJ06VJUFRFh9uzZVFZWsn//fgDy8/Oprq5m27ZtgLWud2pqqn/1vIyMDCZNmuSfCdaJ61RfXw8wYK/Txo0b2bNnT79cp7q6OiorKwfkderv+8n3rBiI1wn6735yBFXtkw14B6tqq/N2dkCYUqAgxP8vAP4Q8Pty4LFw0p4+fbquW7dOfTz11AoViVOsEf7qcrn0vvvu01ghUItYxuhgY7SwMVpYAKu0l8/9PiuhqOrJvYyiChgT8DsL2Bnun8ePH+//PmVKEW73E7S1xeaI+UAtYhmjg43RwsZo4RwDudvwB0COiIwVkQTgYuDVcP/sKz4DlJb6RswvBeaTm3ulw6YObAK1iGWMDjZGCxujhXNEqtvwN0SkCqv71esi8pZ3/ygReQNAVduwJuF6C/gEeFFV1x1Ker4R8y4XwLOsWvX7mBzgaDAYDH1JpHp5vayqWaqaqKojVfU07/6dqjovINwbqjpeVb+qqr/oSRqBXQF9I+bHjStFpIX29tga4Gi6RVoYHWyMFjZGC+cIucDW4YxvgS0f5eUwZw40N5fT3l6CSCvx8fGUlpbG3JgUg8FgCEZfL7B1WBPoUHwj5tvbAQRVEBnwa4Q5RqAWsYzRwcZoYWO0cI6odSi+/utgt6GIlAJtgNLW1hYzVV6BWsQyRgcbo4WN0cI5otahBOJrQ3G5SoAErNMWMjIyImqXwWAwRBNR61AKCzsOqrcGpfrm9HLj8bRz8803x0RPr85axCpGBxujhY3Rwjmi1qFUVVV1+G1PFLkPaAfaaW6OjZ5enbWIVYwONkYLG6OFc0StQ/HN+ePDniiyBKvay43bHRtT2XfWIlYxOtgYLWyMFs4RtQ4lGFdcAUlJRbhci3G7v80ZZ8TWiHmDwWDoS6LWoUyYMOGgfb7G+YICcLme5bXXYmPEfDAtYhGjg43RwsZo4RxR61DcVoNJB2J1SeBgWsQiRgcbo4WN0cI5otah+NZJCKTjksBuQHC73VHfjhJMi1jE6GBjtLAxWjhH1DqUYHScJFIQia0R8waDwdCXRK1DyczMPGhfUREsXgwFBaWItKEaGyPmg2kRixgdbIwWNkYL54hahzJ27NiQxz76qATV2Bkx35UWsYTRwcZoYWO0cI6odSihem6VloLHY4+Yb2+P/hHz0XxuPcHoYGO0sDFaOEekFti6QETWiUi7iIScLllEtonIGhH5SEQcmRLUnijSGjGv2h4TPb0MBoOhr4lUCWUtcC6wLIywJ6lqXk/n6U9KSgq639eOcvbZJbjdCbhcbhISonvEfCgtYg2jg43RwsZo4RxxkUhUVT+Bvu1h1d2Eb2+9VYTqYkSe47TT+syMAYGZ/M7C6GBjtLAxWjhHRBxKD1DgbRFR4ClVXRgqoIjMB+YDjBo1itdff53k5GTGjRtHamoqlZWVAGRkZLBkySSam8W74Naz/O//trBo0dP8+te/ZtKkSRQUFLB79262b98OQE5ODomJiaxduxaAESNGMH78eMrKygBrCdGioiJWrVrlX1uhsLCQqqoq/zxBEyZMwO12+/u8Z2ZmMnbsWH/9bVJSEoWFhaxcuZLGxkYAioqK2Lp1K7t27QJg4sSJeDweNmzYAMDo0aPJyspi5cqVAKSkpFBQUEB5eTnNzc0AFBcXU1paSlycdaknT55Mc3MzmzZtAmDMmDGMHDnSv8hQWloa+fn5lJWV0dbWBsCsWbNYt24d+6wpm8nNzaW2tpYtW7YAkJ2dTXp6OhUVFQAMHTqU3Nxcli5diqoiIsyePZvKykr2798PQH5+PtXV1Wzbtg0g6HWaNGkSy5ZZhdi4uDiKi4upqKigpqYGoMfXqbGxkdNPP33AXqeNGzeyZ8+efrlO9fX1ZGVlDcjr1N/3U0VFBcnJyQPyOkH/3U+OoKp9sgHvYFVtdd7ODghTChR0Ecco7+cIoBKYFU7a06dP1yVLlmgoVqxQTUhQhfsU3Aqo2+3W++67L+R/Dme60iKWMDrYGC1sjBYWwCrt5XO/z0ooqnqyA3Hs9H7uEZGXgZmE1+7SLVZtWwmQgMvVEvXtKAaDwdDXDNhuwyKSLCKpvu/AqVglnLAoKioKeay0FKzSZxEijzBu3BweeeSRLv9zOBOt59VTjA42Rgsbo4VzRKrb8DdEpAprCcXXReQt7/5RIvKGN9hIoExEKoH3gddV9c1w09i6dWvIY/YULOWo3syWLYujeixKV1rEEkYHG6OFjdHCOSLiUFT1ZVXNUtVEVR2pqqd59+9U1Xne71tUNde7TVLVX/QkDV/DWzB8XYdPPrkUl6uF9nYPTU1NPPfcc706r4FKV1rEEkYHG6OFjdHCOQZslVdfU1QE551Xgog1dbWq8vTTT0dtKcVgMBj6mqh1KBMnTuzyuLU2ShHt7VcD1niYaJ0osjstYgWjg43RwsZo4RxR61A8Hk+Xx+21Ua4ABiESvSPmu9MiVjA62BgtbIwWzhG1DsU3qCwUdsN8EW73Yk488dtceWV0rjHfnRaxgtHBxmhhY7Rwjqh1KN3hW1/e7Yb2dli27E889dRTlJSUmHYUg8FgOASi1qGMHj262zD79lnORPU5oAVVpaWlJep6e4WjRSxgdLAxWtgYLZwjah1KVlZWt2FKSqwSSmeirRthOFrEAkYHG6OFjdHCOaLWofgmeOsOawqWK4B4/75FixZFVbVXuFpEO0YHG6OFjdHCOaLWoYRDxylYrsHXfbi5uZkHHngggpYZDAbD4UfUOpSUlJRuw/h6erndEB9/Be6A+q9XXnmFhQtDzpZ/WBGOFrGA0cHGaGFjtHAOsWYtji4KCgrUtx5Bd5SXg68NftmyQtavf99/bObMmaY4bDAYYgIR+VB7uDJuZ6K2hNKTNpBnn4Xf/x42brymw/7Vq1dHRVtKNJyDExgdbIwWNkYL54hah+JbYa07fCPmPR7weOaTnX2O/1hra2tUtKWEq0W0Y3SwMVrYGC2cI2odSrgEdh1WhR07bsPtttcde+WVV7j99tsjY5zBYDAcRkRtG8p7773nX0e9O77zHViwwPrucsGoUYVUVdltKSLCggULmD9/fl+Y2+e0tbWFrUU0Y3SwMVrYGC0sTBtKF2zcuDHssNOm2d/b22H27I5tKarK9ddff9j2+uqJFtGM0cHGaGFjtHCOSK3Y+KCIfCoiH4vIyyIyJES4uSKyQUQ2i8gdPUljz549YYfdt88qmVhpQmrqfG677bYOYVSV6667jtmzZ/Od73znsGrI64kW0YzRwcZoYWO0cI5IlfP+Bdypqm0i8ivgTqBDQ4VYK189AZwCVAEfiMirqrreaWNKSiAuzjedPTz9NCxZ8is2btzIK6+80iHssmXLWLZsGQsXLqSgoICSkhJqamrYtWsXmZmZXHHFFWaNaoPBEJNExKGo6tsBP98Dzg8SbCawWVW3AIjIC8DZQFgOZfLkyWHbU1QEV18NTz1lOZS2Nqv312233cbrr79Oa2vrQf9pb2/n/fff5/333++wf+HChfzud79jypQplJaWUlJS0sHBlJeXd9jf+XdgmIyMDPbt23dQHKHiCrWvJ1p0R7D4I0FP7SgvL+edd94hMTGxx3aHk1Z3YcK1NzAc0G1eOVSbfXmip3H6Jk71vTiFyoPdhevPfNRVWr3JF71JN9TxYNodStzhpuX0NYh4o7yIvAb8TVX/u9P+84G5qnqt9/flQKGq3hQinvnAfIBRo0ZNf+aZZ4iPj2fcuHGkpqZSWVkJQEZGBpMmTWLZsmUAxMXFUVxczDPPbOD664+muVkQEQoKGrnyyk+Bcn7729+yadOmHp2X2+2mvb0dEWHMmDFkZmYiIrz//vu0t7fjdruZNGkSa9asQVUREU488URqamr4+OOPaW9v950XLpeLCy64gISEBFSVwsJCAH70ox/R3NyMqjJ16lTOOOMMHnnkEZqbm4mPj+fGG2+krq6OTz/9FI/HQ3p6OiUlJXz00Uds2bKF9PR0CgsL2bJlCx988AGtra2MHTuWffv20dzczJFHHklOTg719fVUV1czZswYnnjiCdra2oiPj+ehhx5i0qRJZGdns3nzZv72t79RV1fH2rVrSU9PZ8iQIRx99NFs3ryZAwcOUF1dzeTJk6mvr/fP7FxfX096ejrHH388GzduZOPGjaSnp3PhhRdy6aWXcuutt/LGG28wfPhwfvWrX5GYmEh5eTk//OEPaWtrw+12M2PGDGpqanC73Zx44ok0NjaSlpZGe3s748ePZ9GiRbz55pt4PB7i4uIoKiqiurqampoapk6dysyZM3nzzTfZt28fc+bMoba2ltWrV9Pa2srRRx/d4b9nnnkmc+fO5c033/RXlYgIK1asoL29HZfLRW5uLuPGjWPdunUkJCSQkJDA6tWr/cePPPJIxowZQ3p6OuPGjWPTpk24XC6qqqr84eLi4nC5XLS2tiIiTJ8+nY8//pjW1lbcbjeTJ0+mqamJhIQEUlJSSE9P54wzzmDx4sXs3LnTn89EhKOPPpq8vDzGjBlDa2srb775JvHx8cTHx1NZWem3a/z48Zx11lk0NjYSHx9PTU0NU6ZMoaqqiu3bt/Piiy/686Xb7eb73/8+Tz75JC0tLYgIF1xwATt27ODdd9/1h3O5XBQVFfHee+/h8XgQEUaPHs2uXbtos+Y94qtf/Srf+MY3OHDgAHl5eeTn51NUVMSqVav461//yhtvvEF8fDwNDQ0kJyeTmppKamoqI0aMYNiwYWzatImGhgbGjRtHbW0tmzdv5uSTT+b+++/nzjvv5MEHH/RfvyuvvJI9e/ZQXV1NWloaixcvpq2tDRFh0qRJnHzyyWzatInq6mpGjhzJaaedxqJFi9i7d6//fHbs2IGIEB8fz7e//W2OPvpoXnvtNdLS0lizZg3//ve/qaurY+fOnQDEx8dz7LHHcuDAAb7yla9QVFTkn9XcZ9e1117LmjVrWLFiBb5ncnx8POedd16H++nUU09l69atPProo379Ro4cyfDhw3G5XBx77LHU19fT3NzMqFGjSE9PZ9OmTbz11lu0tbWhqowcOZKUlBS2bt2Kx+PB5XJx/PHHU1ZW9qmqHtujB11nVLVPNuAdYG2Q7eyAMHcBL+N1bJ3+fwHwh4DflwOPhZP29OnTdcmSJdpTbrtN1SqjWFt8vOqKFaorVqzQxMREBQbMlpGREXEbhg8frrNmzdJZs2apy+VyPP709PQOv10ul9522206c+bMiJ+72dBBgwY5Hqfb7dZjjz1W8/Ly9KijjupVXJmZmX2SLztvItKvuvdheu29fe73WZWXqp7c1XERuRI4E5ijwYtJVcCYgN9ZwE7nLDyYjz7q+Lu11ar6uvPOIpYsWUJpaSkHDhzg4YcfDloN1p/s27cvoukD7N271//m1hdUV1d3+N3e3h4VA02jhaamJsfj9Hg8fPLJJ47E1V/LUPR3LU8fpie9jSAibSgiMherEX62qjaECPYBkCMiY4EdwMXApeGmMWbMmO4DdeK88+DtgNad+HirwR6gqKjIX894zjnn+Ns4Vq9eDcC0adN4/vnn/VVpBoPBEHP0VZVXN9Vhm4HtwEfebYF3/yjgjYBw84CNwGfAXeHGP336dK2tre1xlZeq6lNPqc6cqZqXp3rssarnnGNVe4XLbbfd1qNitsvl0ry8vF4VzTtXDR3q1t9F90hvTunm0y7YNeyppj0Jn56efsj55lCvtYiE9d9QWgRubrdbc3JywoovMzNThw8fHvE8A1ZVb05OTpe25uXlHWRvamqq47a4XC7Nzs4OW0dAhwwZEjJsb5/tEW+U7wsKCgr0oYce8veU6SkLF8J119m/XS743e8g3IHynXtpZWRkcPPNN9PS0oLb7ebqq69m2rRpHXpw+Xp3PP300/5Gwvz8/IO6JaelpVFaWkpLSwsJCQlcc801zJ8/n4ULF/LHP/6RQYMGkZ6eTnV1NU1NTZSUlFBdXU16ejoPP/ywvyH7lltuOai7M9i9igAeeOABXnvtNX+nAV+mcblcTJkyhf3797N9+3Z/o25xcTETJ070x/Xcc88dFL+vB4vv/H2lPF+4adOmsXr1atavX8/evXuZMGECp59+OjfeeKO/EdKHzyaAxMREHn30UX9cwdi7dy/Dhw/v0L174cKFvPTSS5x33nlMmTKlQw+bQHvT0tL46KOPGD58OHv37iUvL48hQ4Z06Inn0y/UPl9pNpg+gfnBp/3OnTspKSlhyJAhHDhwoEP65513HvPnz/fnNd/xvLw8/3X1sWjRIn/j/llnncVtt91GRUUFNTU1B5WyFy1axGuvveZv0A94ucPtdvPkk08yZcoUHnjgATZs2MDw4cOZOHGiX59AXW6++Waam5txuVzccsstDBkyJGTvNV9eC2y49zVqJyYmsnjxYgDmzJnjv4/mzZsH0CHf+M7j+eefp6ysDFUNmt/T0tL49a9/jcfjOShPud1uLrroIl588UXa2tpwuVy4XC7a29tJSEhg8eLFHezeuXMnOTk5Ha6L7zlw0kkn+e/VJUuWHHTd09LSDrovfde7tLSUUaNG+cfEPfDAA7z66qv+c7r22ms79AQL7CHm02PXrl1UV1fz7rvvoqodtOz8jLruuut2qGqvlq80DiUIp53WseoLeu5UOnMo3Uad6srni+9Q4u6qG2tf2RvKjq4cVDjjfwLPI5YIdo260qLzNe+Jxt2l25Pwve3qGk633eeee441a9YwZcqUoC95XeV9pzRw8px68z8npl6JWoeycOFC8vPzD+n/nUsoPkQgNxeam2H4cJg4Ea64whrHMpCpqKg4ZC2iCaODjdHCxmhhYRxKCHqywFYoFi6E++6Dzz/vOpzbDWedZX3PzLTmBdu3z2rMH+iOxgnKy62ecL05XyfiGAh0Po/+OC9fGhkZsZXvDM7jhEOJ2ik2y8rKKC4uPuT/z59v3aB33QVd+VyPBzrNzgJYPcTOOCO8tAId0YEDVvfl4cNh716r55mvmi1wdcnOJSPfMV/VeWamFWbNGnjkkRomTEjjttvsB11gPBD6wdc53kB7Fy2C116zJtR0uSzH6kvD999QD7vAYzffbJX6XC544omO1YqdH8oLF8JLL9m6BEvDdz4ZGbB6tWV7Zibk5lZy/fW5B51bZz07x9nZ/mBap6XBww9bsyy43XDLLfDYY9Z5iVjanH568PiCae27fsHOJfAafO97Vhpg6ZeYCIsX2/EG6jVlih2/6hfcfvsRQOg81d01DEYoJ9r5unX3/2DaB7M13PBdvegtWFDJW2/lsmGDdd+lp9sad04n1Pn31O7utIOD0wvMr77z6SqPdrZn9WpYvx6amuCaazreP07VBEdtCaU3bSg+ysthzhxobHTGrkMlOxtOOAFeeMFyYIH7hwyB/fvh3/8+2PG5XNbD3urAYXUxz8qCHTs6hhWxfrvdVjrp6VBdbZXOuiuhhbI3Ph4++8yXvkVcHFx0EVRUwKefWmn60g5k1izrc+9e2LzZikMEhg3r6Njy8mDdOmu8UCDB4rRoJyfHxdChlm7/+pcdTgROPNE690WL7Id0IC4XFBTABx90/ZIRLr40L7sMnn8eli8/ON7Q5xKarCy49FJ47z0I7MXeMS7F7Rba2ztqcMop1ktNS4uVr7Zv73gNXS7wvaft3Ws5sOZm69MX3hcuP996cH32GQQOH8rMtLaEBMjJgU2bQqfns8vttvYHHsvMhN27e6ZPYiIsWWK9aP3xj1aamzbZ90dnXC6YOhU+/tjOh0cdBUce2fE+CWZ3KHz3bWDVeW0t/PWvdhq+8wUrrSFDbBu6QgRGjoQ9e7oPm5dnx2k9K8YfUN04NLyzCJF+tDqUtrY2hgwZ0mH/hRdeyA033EBDQ4O/h0ggV111FVdddRVffvkl559vTS9WU2M9xKy3uu8gchHDhm1n797Lg6T8Q+AsYAMQpBGGHwMnY/WUvjnI8fuA44EVwP8LcvwRIA9rEoJ7gxx/CpgAvAb8OsjxP2ONFf0b8Lsgx/8ODAOe8W6deQMYDDwJvBjkeKn38yHgn52OJQGLvN/vARZ3Op4BvOT9fifQeTbnLMA3O8/NWBoGMh7wLS8wH6u3eSB5wMNYD45vYo2bDaQIuN/7/Tyg88DROcBPvN9PBzq/ZZwJ/Mj7vYSDuRC4AWjA6g3fmau825cEn9ruO8BFWL3tTd47mFLvZ/d5z+1eTMfOXf2R9x7xfh/Iea8A1VW9GtwYtVVenZ3JoZKWZm2ZmTBhglXFkJUFJ58MZhmFw41eDwQ2RAGdegobHCRqSyhPP/00U6ZM6dN0Auu8Fy2yq158s4ENHEIX6WMLo4ON0cLG0iI93aoCG1j3bn9iSigh6Y+5roqKDm7IDdUY2Llx7Omn7QbcwNq3wAb6jIyD69Z9vcqqqzvWYfvqYn2N5YEdBS6++N+cdNJR/PGPVl11QoJlZ02N1Ui3d69V2gp8cwusK/bF6+s00LnxeciQjo1+PrsC68g//LBj/Hl51mfg/GmBdbqhmDULjjsOHnqoYzgRqx54yhQ73b17rTiHDLHsfvXVWjyeNH/bjssFl1xi2ddpFQLAPm/fdGXvvNPxf5MmHdxYHvhiccwxcOaZdvqlpVb7kcdjxTF8+MGdHfbts+vRc3IObiQOvA6lpTBqFIwfb18T6PhA9Nn6wgtWvHFxVp4NbFuZNQvGjOlYh3/00VbYwPSrq6GsrGNbgq8twNeGEqwtz8dll0F9PWzYYNnquw6B6SUmWu0bvv2dX8465/9A+wKprrYanwPzebA2qsxMGDKkhh/84Cv+RurnnrPuT9/6SD4dfW1Hn39un2fgfQLw7rt2HjnhBGtfZzt890dnvURgxoyObXS+NpHdu207oGMagecf+FwI1NFHXJwVty+PBrYdvvJKda8n5ovaEooTjfJ9SU+6lHbVuysUgb1qxo/vfkBfT3rzHEp32FDnEKrXVqBjDuz11FWPl+7ssQd4Hty9d84c6wEiYjcmd+6NFM55h9sbqLteaz3t/hvqhSZUL6Uf/nADa9dOOKgXoRPnF3hdFi2CnTt7pme4eaCnlJfDHXdYThHsHnHNzQffH93dD+HYHo6dwe6LUPu66gEWKu5g/4Hg96IZhxKCgoIC/de//sXQob3qsBA17N+/32hB1zpEy1iYcInlPNH5WseyFoGYcShdUFtbazKJF6OFRVc6BFZfxgKxnCc6X+tY1sJpXJE2oK/YsmVLpE0YMBgtLIwONkYLG6OFc0StQzEYDAZD/xKVbSgisheoxxolZrBGjBktjA6BGC1sjBYWE1Q1tTcRRGUbiqoOF5FVvW1gihaMFhZGBxujhY3RwkJEejejLqbKy2AwGAwOYRyKwWAwGBwhmh3Kwu6DxAxGCwujg43RwsZoYdFrHaKyUd5gMBgM/U80l1AMBoPB0I8Yh2IwGAwGR4g6hyIic0Vkg4hsFpE7Im1PXyMifxKRPSKyNmBfuoj8S0Q2eT+HBhy706vNBhE5LTJWO4+IjBGRJSLyiYisE5Hve/fHohaDROR9Ean0avFz7/6Y08KHiLhFZLWI/NP7Oya1EJFtIrJGRD7ydRN2VAtVjZoNcAOfAeOABKASmBhpu/r4nGcB+cDagH0PAHd4v98B/Mr7faJXk0RgrFcrd6TPwSEdjgDyvd9TsZbNmxijWgiQ4v0eD6wEjotFLQI0uQX4C/BP7++Y1ALYBgzrtM8xLaKthDIT2KyqW1S1BXgBODvCNvUpqroMqO60+2zgWe/3Z4FzAva/oKrNqroV2Iyl2WGPqn6hqhXe77XAJ8BoYlMLVdU6789476bEoBYAIpIFnAH8IWB3TGoRAse0iDaHMhpr0W0fVd59scZIVf0CrActMMK7Pyb0EZFsYBrWm3lMauGt4vkI2AP8S1VjVgusBd1vAwKXbotVLRR4W0Q+FBHfCjWOaRFtU68EW77S9Iu2iXp9RCQFeAm4WVVrRIKdshU0yL6o0UJVPUCeiAwBXhaRyV0Ej1otRORMYI+qfigiJeH8Jci+qNDCywmqulNERgD/EpFPuwjbYy2irYRSBYwJ+J0F7IyQLZFkt4gcAeD93OPdH9X6iEg8ljN5XlX/4d0dk1r4UNUDQCkwl9jU4gTg6yKyDasK/Gsi8t/Ephao6k7v5x7gZawqLMe0iDaH8gGQIyJjRSQBuBh4NcI2RYJXgSu9368E/jdg/8UikigiY4EcIMhq6ocfYhVF/gh8oqq/CTgUi1oM95ZMEJEk4GTgU2JQC1W9U1WzVDUb63nwf6r6TWJQCxFJFpFU33fgVGAtTmoR6V4HfdCLYR5WD5/PgLsibU8/nO9fgS+AVqw3imuADGAxsMn7mR4Q/i6vNhuA0yNtv4M6FGMVxz8GPvJu82JUi6nAaq8Wa4GfevfHnBaddCnB7uUVc1pg9X6t9G7rfM9HJ7UwU68YDAaDwRGircrLYDAYDBHCOBSDwWAwOIJxKAaDwWBwBONQDAaDweAIxqEYDAaDwRGMQzEYDgER8XhnbPVtjs1sLSLZgbNHGwyHC9E29YrB0F80qmpepI0wGAYSpoRiMDiId72JX3nXI3lfRI727j9KRBaLyMfezyO9+0eKyMvetUsqReR4b1RuEfm9dz2Tt70j3g2GAY1xKAbDoZHUqcrrooBjNao6E3gca6ZbvN+fU9WpwPPAo979jwJLVTUXa12bdd79OcATqjoJOACc16dnYzA4gBkpbzAcAiJSp6opQfZvA76mqlu8k1XuUtUMEfkSOEJVW737v1DVYSKyF8hS1eaAOLKxppzP8f6+HYhX1Xv74dQMhkPGlFAMBufREN9DhQlGc8B3D6a903AYYByKweA8FwV8lnu/r8Ca7RbgMqDM+30x8B3wL4qV1l9GGgxOY956DIZDI8m7IqKPN1XV13U4UURWYr2wXeLd9z3gTyJyK7AX+JZ3//eBhSJyDVZJ5DtYs0cbDIcdpg3FYHAQbxtKgap+GWlbDIb+xlR5GQwGg8ERTAnFYDAYDI5gSigGg8FgcATjUAwGg8HgCMahGAwGg8ERjEMxGAwGgyMYh2IwGAwGR/j/YDwm+Ww9B1QAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig = plt.figure()\n",
    "ax = fig.add_subplot(111)\n",
    "fig.canvas.draw()\n",
    "num_qub = 2\n",
    "num_lay = 1\n",
    "H = create_Hamiltonian(num_qub, J, h)\n",
    "exact_energy = np.sort(np.real(np.linalg.eig(H.to_matrix())[0]))[0]\n",
    "print(\"Exact ground energy: \", exact_energy)\n",
    "\n",
    "losses = np.loadtxt('Project_QC/Ising_model/shot-based_backend/without_noise/%i_qub/ground_%iqub_%ilay_ryrz_linear_noNoise_losses' % (num_qub, num_qub, num_lay))\n",
    "print(\"VQE ground energy: \", np.mean(losses[-50:]))\n",
    "print(\"Accuracy: \", 100*np.abs((np.mean(losses[-50:])-exact_energy)/exact_energy), \"% \\n\")\n",
    "ax.plot(losses, 'b.', label = 'VQE energy without noise', color = 'blue')\n",
    "\n",
    "losses = np.loadtxt('Project_QC/Ising_model/shot-based_backend/with_noise/%i_qub/ground_%iqub_%ilay_ryrz_linear_yesNoise_losses' % (num_qub, num_qub, num_lay))\n",
    "print(\"VQE ground energy: \", np.mean(losses[-50:]))\n",
    "print(\"Accuracy: \", 100*np.abs((np.mean(losses[-50:])-exact_energy)/exact_energy), \"% \\n\")\n",
    "ax.plot(losses, 'b.', label = 'VQE energy with noise', color = 'black')\n",
    "ax.plot(exact_energy*np.ones_like(range(len(losses))),'k--', label = r'Exact $E_{0}$')\n",
    "\n",
    "ax.legend()\n",
    "ax.set_xlim([-2, 500])\n",
    "title = r\"VQE ground energy with QC noise model for $n = %i$ qubits\" % num_qub\n",
    "ax.set_title(title)\n",
    "ax.set_ylabel(\"Energy (a.u.)\")\n",
    "ax.set_xlabel(\"Epoch\")\n",
    "ax.set_ylim(top=2)\n",
    "#ax.set_facecolor('#e4e6e8')\n",
    "#plt.grid(color='white', linestyle = '--', linewidth = 2)\n",
    "plt.grid(linestyle = '--')\n",
    "fig.canvas.draw()\n",
    "plt.savefig(\"with_without_noise_shot-based_%i_qub.pdf\" % num_qub)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "gross-marina",
   "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}