Parameter_Shift_Rule.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Information about the repository QOSQ Previous Challenges.\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Parameter Shift Rule Derivation using Qiskit and validate with Pennylane \n",
    "\n",
    "## Find the lowest eigenvalue of the followinf matrix:\n",
    "\n",
    "I will use Qiskit for the implementation.\n",
    "At the beggining I just need to calculate the gradient descent in the circuit and compare them with the pennylane library.\n",
    "\n",
    "## Background:\n",
    "\n",
    "First we need to understand the problem in order to tackle it in teh best way. I will begin showing what it is a VQE (Variational Quantum Eigensolver). How we can implement parameter quantum circuits and expectations values. At the end, how can we optimze the parameters in order to find the value or the expected behaviour.  \n",
    "\n",
    "Variational Theorem:  \n",
    "If you are familiar with hamiltonians and expectations values. You will find that the eigenvalue attach many interest in many research areas, in particular we are interested into the lowest value in the hamiltonian. The variational theorem states that:  \n",
    "\n",
    "$$  \\langle \\psi | \\hat{H} | \\psi \\rangle  = \\lambda_1 \\ge \\lambda_0 $$\n",
    "\n",
    "This equation tell us that any evaluation (expectation value) of a random state $\\psi$ in our hamiltonian will be greater or equal to the lowest eigenvalue possible of my hamiltonian. Lets assume that we parametrized the state $ \\psi $, now the equation has the following form:\n",
    "\n",
    "$$  \\langle \\psi(\\theta) | \\hat{H} | \\psi(\\theta) \\rangle  = \\lambda_1 \\ge \\lambda_0 $$\n",
    "\n",
    "In principle, the problem is still the same as before, but now we can modify the parameter $\\theta$ in order to find the eigenvalue. Now, we need to introduce the concept of parametrized quantum circuit. As you know (or maybe not) the quantum circuits can be parametrized in order to modify some gates. The perfect gates to do this are the $RX$, $RY$, $RZ$ and, the control version of this three gates which are $CRX$, $CRY$, $CRZ$ or any linear combination of theses. This will help us to created our $\\theta$ state.  \n",
    "\n",
    "The second step is to evaluate the quantum state with our hamiltonian. The hamiltonian is consider that main part of our research, becouse we need to deal with it in a very carefull way to reduce some complexity. For our problem our hamiltonian is very simple but take care when you are planing to deal with a bigger one.  \n",
    "\n",
    "This is our hamiltonian to find the lowest eigenvalue.  \n",
    "\n",
    "$$\n",
    "\n",
    "\\begin{equation}\n",
    "\\begin{bmatrix}\n",
    "    0 & 0 & 0 & 0 \\\\\n",
    "    0 & -1 & 1 & 0 \\\\\n",
    "    0 & 1 & -1 & 0 \\\\\n",
    "    0 & 0 & 0 & 0 \\\\\n",
    "\\end{bmatrix}  \n",
    "\\end{equation}\n",
    "\n",
    "$$\n",
    "\n",
    "If you carefully this is a non-unitary matrix. This can be anoying at the moment to insert into a quantum circuit.Becouse we work with unitary transformations. Nevertheless, a very interesting properties about the expectation value is that this is linear, such as:\n",
    "\n",
    "$$  \\langle \\psi(\\theta) | \\hat{H} | \\psi(\\theta) \\rangle  = \\langle \\psi(\\theta) | \\hat{H}_1+\\hat{H}_2 | \\psi(\\theta) \\rangle =  \\langle \\psi(\\theta) | \\hat{H}_1 | \\psi(\\theta) \\rangle +  \\langle \\psi(\\theta) | \\hat{H}_2 | \\psi(\\theta) \\rangle $$\n",
    "\n",
    "Now we can convert our non-unitary matrix into two separete unitary matrices. Caviat; this is not the only way to convert this matrix into a unitary.\n",
    "\n",
    "$$\n",
    "\n",
    "\\begin{equation}\n",
    "\\begin{bmatrix}\n",
    "    0 & 0 & 0 & 0 \\\\\n",
    "    0 & -1 & 1 & 0 \\\\\n",
    "    0 & 1 & -1 & 0 \\\\\n",
    "    0 & 0 & 0 & 0 \\\\\n",
    "\\end{bmatrix} =\n",
    "\n",
    "\\begin{bmatrix}\n",
    "    1 & 0 & 0 & 0 \\\\\n",
    "    0 & -1 & 0 & 0 \\\\\n",
    "    0 & 0 & -1 & 0 \\\\\n",
    "    0 & 0 & 0 & 1 \\\\\n",
    "\\end{bmatrix} +\n",
    "\n",
    "\\begin{bmatrix}\n",
    "    -1 & 0 & 0 & 0 \\\\\n",
    "    0 & 0 & 1 & 0 \\\\\n",
    "    0 & 1 & 0 & 0 \\\\\n",
    "    0 & 0 & 0 & -1 \\\\\n",
    "\\end{bmatrix}  \n",
    "\\end{equation}\n",
    "\n",
    "$$\n",
    "\n",
    "Perfect, now our matrix in divide into two unitary matrices. Now in order to optimize the quantum parameters to find the lowest eigenvalue we need to find the gradient of our quantum circuit with respect to the parameters. Maybe you are wonder, how can I achieve this? In the resource section you will find a more detail explanation [1], here a I will give you the short part.  \n",
    "\n",
    "We will build something call the Parameter-Shift Rule in order to find the derivate in our circuit with certain parameters. So, lets assume we have $N$ parameters in our quantum circuit and we want to find the derivate/gradient in the $N_i$. With the parameter-shift rule we need only evaluate our circuit in the following way:\n",
    "\n",
    "$$  \\nabla_i  \\langle \\psi(\\theta_i) | \\hat{H} | \\psi(\\theta_i) \\rangle = \\frac{1}{2} [ \\langle \\psi(\\theta_i+\\frac{\\pi}{2}) | \\hat{H} | \\psi(\\theta_i+\\frac{\\pi}{2}) \\rangle - \\langle \\psi(\\theta_i-\\frac{\\pi}{2}) | \\hat{H} | \\psi(\\theta_i-\\frac{\\pi}{2}) \\rangle ] $$\n",
    "\n",
    "As you can see we are evaluating the quantum circuit and calculating the expectation value just moving one parameter at the time and the other are constant. With this tecnique we just need to iterate for all the parameters and we will get all the gradients for all the points. With the gradient we can apply a simply gradient descent algorithm to find the lowest value in our hamiltonian as follow:\n",
    "\n",
    "$$ \\theta_{new} = \\theta_{old} - l*\\nabla_i \\langle \\hat{H} \\rangle $$\n",
    "\n",
    "Finally, just iterate many times to find the lowest eigenvalue!\n",
    "\n",
    "## Resources:\n",
    "[1] Parameter Shift-Rule: https://pennylane.ai/qml/demos/tutorial_backprop.html\n",
    "\n",
    "## Comments:\n",
    "- Remember that the just need the gradient in the circuit to calculate the values.\n",
    "- Remember the order of the operators.\n",
    "\n",
    "### Importing\n",
    "\n",
    "We will work with Qiskit and Pennylane. Pennylane is a very focus framework to work with Quantum Machine Learning projects.in the other hand, Qiskit is more friendly and it is used for most student when are involved in the quantum computing world. We will use qiskit to build the circuit a test the problem and pennylane to corroborate our problem and results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from qiskit import QuantumCircuit, Aer, execute, QuantumRegister, transpile\n",
    "from qiskit.circuit import Parameter, ParameterVector\n",
    "from qiskit.visualization import plot_histogram\n",
    "from qiskit.quantum_info import Operator\n",
    "from qiskit.opflow import CircuitStateFn, CircuitOp\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We need a quantum circuit with just 2 qubits. Why? because our operator (matrix) is a 4x4 matrix which is $2^n$ x $2^n$ where $n$ is the number of qubits. Solving this we have that $n$ must be equal to 2. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAukAAAEfCAYAAAD851iSAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAA9zUlEQVR4nO3dd3hUddrG8XsS0klIAiQBEhOaoYQiTQUURAKIjWIDF1ZdsayIDXfVfXdFXURXoqisrLIqCvaAKAjSe5MWOoSAQOg1jRAIybx/zDpkTJ1kJnNm8v1cVy5zZk55gjzMnXN+53dMZrPZLAAAAACG4eXqAgAAAADYIqQDAAAABkNIBwAAAAyGkA4AAAAYDCEdAAAAMBhCOgAAAGAwhHQAAADAYAjpAAAAgMEQ0gEAAACDIaQDAAAABkNIBwAAAAyGkA4AAAAYDCEdAAAAMBhCOgAAAGAwhHQAAADAYAjpAAAAgMEQ0gEAAACDIaQDAAAABkNIBwAAAAyGkA4AAAAYDCEdAAAAMBhCOgAAAGAwhHQAAADAYAjpAAAAgMEQ0gEAAACDIaQDAAAABkNIBwAAAAyGkA4AAAAYDCEdAAAAMBhCOgAAAGAwhHQAAADAYAjpAAAAgMEQ0gEAAACDIaQDAAAABkNIBwAAAAyGkA4AAAAYDCEdAAAAMBhCOgAAAGAwhHQAAADAYAjpAAAAgMEQ0gEAAACDIaQDAAAABkNIBwAAAAyGkA4AAAAYDCEdAAAAMBhCOgAAAGAwhHQAAADAYAjpAAAAgMEQ0gEAAACDIaQDAAAABkNIBwAAAAyGkA4AAAAYDCEdAAAAMBhCOgAAAGAwhHQAAADAYAjpAAAAgMEQ0gEAAACDIaQDAAAABlPL1QUAFbH+K+litqurcC9+wVLnIa6uAp6GXrQfvQhnoR/t4269SEiHW7iYLeVluboKAPQiYBz0o2djuAsAAABgMIR0AAAAwGAI6QAAAIDBENIBAAAAgyGkAwAAAAZDSAcAAAAMhpAOAAAAGAwhHQAAADAYQjoAAABgMIR0AAAAwGBquboAwNWem9RTW/cvK/E9L5OXAv1DFFo7Qs0adVDX1neqe8Ig+dTyreYqgZqBfgSMgV50PUI6UIZCc6FyLmQo50KGDp9K1dKUr9WwXjO9NPQrxcd0cnV5QI1CPwLGQC9WD0I6UER8TGfFx3SxLpvNhco8f1o7D67W6cwjkqSjp9P01496692RaxQb2dJVpQIej34EjIFedA1COlBElxb9NbzPmGKvFxQWaNaaSZr0w1MqNBfqfF6mJv34tN4YMa/6iwRqCPoRMAZ60TW4cRSoAG8vbw3oNlL39PyL9bWNqfN1MiPdhVUBNRP9CBgDvehchHTADnd0fcJmefuvK11UCQD6ETAGetE5COmAHeqHRis4IMy6fDbrmAurAWo2+hEwBnrROQjpgJ18fQKs31+6nOfCSgDQj4Ax0IuOR0gH7HC5IF9ZuWesy6G167uwGqBmox8BY6AXnYOQDthh268rlH/5onW5Rcy1LqwGqNnoR8AY6EXnIKQDFXQpP0+fzHnRunx1dCc1adjWhRUBNRf9CBgDveg8zJMOlKGwsFBZuWe0bf9yfbHon9p3NEWSFBJYV6Pv/dS1xQE1DP0IGAO9WD0I6UARUxe8oqkLXilznVax1+vJgR+ocVRCNVUF1Ez0I2AM9KJrENIBO8RGttILQ6apQd0mpa6zdf9yjf5PT3mZvDVx1Ho1a9S+xPV2HFitZz+4QTKZ9N7ItYqP6eSkqgHPVF4/jp12n5Zu+Ubtm96kfz26SCaTqcT1zGaz/vLhzUrZt0Q929+nv93/lTPLBjxORT4bn5vUU1v3L7Nrv+MfW6J2TXtWsTr3RUivZgUFBVq9erX27t2rkydPKjw8XDExMbrxxhsVFBTk6vJqvPiYzoqP6WJdPn8hQ0fOpGn3oXWSpIMndmrU+9dp/ONLFRvZqsR9tG1yo2677jHNWjNJSd89pIlP/iJvb9tWu5Sfp6Tv/qRCc6Hu6fEXAno1KjRLqcel3Uel9LPSiUwpv0Dy9pLCgqSYcKlJhNT+KsmXfyFdqqr9+MSA97U5bZFS9i3R7LUf6vbrHyvxOD+t/Ugp+5YoNKi+Rg543zk/DEp0LENKOSSln5GOZkgXLkkmk1TbX4oOk2LrSR1ipTqBrq60ZnPEZ2PjBm1UUHi53GMdOZWqjPOnJElB/qFVL96N8RFUTTIyMjR27Fh99tlnOnXqVLH3Q0JCNGLECI0bN04+Pj4aM2aMXnnlyqWlHTt2qFWrkv/iw3G6tOiv4X3GFHt9/9GtGvflUB04sUMZ50/plc8Ha9LTm+RXZF7Yoh7u/6bW7pqttCOb9e2ytzSk14s273++YIzST+5WdP2r9cc+ZV9ChGMUFkqr06Slu6TTOSWvk3tJOnJOWrtPmrFBurap1DdBCvSr3lphUdV+DK1dXyMHTNTYL+7Tf3/6i65teasiQmNs1jmVcViT51geaT5y4ETVCarntJ8HV+w5Js3bJu0v/nEoScrLl05nWwL8rM1SmxjplrZSVJ3qrRMWjvhsrMgvwIdO7tbI9zpLkgbf+GypV6JrCmZ3qQZz5sxR8+bNNX78+BIDuiRlZWUpKSlJ/fv318WLF7Vlyxbre0FBQWrRokV1lYsSNGnYVmMfnqvgwHBJUvrJ3Zq28LVS1w/0D9bTgz+UZBnLd+jkbut7qYc36rtl42UymfTsXf+Vr4+/c4uHTmVJ7y+UkteXHtB/Ly9fWrZbevMnaccR59YH+9jTjz3b36turQco92K23vluRLH3J0x/VLl5WeqWMFA92t3j1Lph6auv10qTFpce0H+v0CxtOSSNnyMt3GH5hRvGYO9nY1kuXDqvVz8frAsXc9Q6rptG9H/TkaW6JUK6k33yySe64447dPr0aetrXl5e6tSpkwYPHqxBgwapadOm1vcWLlyoDz74QCkpKdbX2rVrJy8v/le5WkRojB7sN9a6PHPFuzqXfaLU9bu0uEW9OwxT/uWLSvr2TyosLNTlgnyN/+ZBFRYW6Pbr/6w2TW6ojtJrtF9PSUk/W/5bGZkXpMlLpSW7HFoWqsiefhw1aJKCA8K0IXWefl5/ZeaJBRun6pfdcxQcGK5RAz9wes01XfYF6d35litVlXG5UJqdIk1ZKV0ucGhpqAJ7PxtL827yozp4YqdCa0fo//7wbbFhojURyc+JZs+erREjRqigwPKviZeXl5599lkdOnRI69evV3JysqZPn660tDTNmjVLYWFhkqTXXntNBw4csO6nY8eOrigfJejf5WFFhTeWJOXl5+rrJW+Uuf7jd05QaO0I7Ty4WjNXva8vF43Vr8e3KTIsVg/3L3tbVF36Wek/iy1n76rqh03Sij1V3w8cp6L9GB4SpcfvmCBJ+nDWszqdeVTnsk9o0o9PS5Iev2OCwkOiqqPkGiv3kvTBIssY9Krami5NW205ww5jsPez8fdmrZ6kRZu/kJfJSy8O/VL16jR0Rpluh5DuJOnp6Ro2bJgK/3ddzt/fX/PmzVNSUpIaNWpUbP3bbrtN33//vUwmk86dO2fzHiHdOLy9a+nenn+1Ls9ZN1lZ58+Uun5IYLhGDpgoSfp07kv6avHrkqSnBn+oAL/azi22hrt0WfpspXSx/PuUKuz7jdLhs47bH6rGnn5M7DRcXVr0V86FDL074zG9P3OksnPP6tqWtyqx47DqKrnGmrFeOpbpuP2lHJJWpTpuf6gaez8bi0o9vFGTZj0jSRre91V1aH6zU2p0R4R0J3n00UeVkZFhXf7iiy/Uu3fvMrfp0aOH+vbtW+z1Dh06OLo8VEGfzg+obojlt/y8S+c1Y8WEMtfv0e5udWs9QHn5ubpckK8+nR5Q5/ji/5/hWD9tsdx45kiFZunLNVIBY2INw55+fHrwhwr0D9HanbO0YmuygvzrWO8dgfNsPyxtOOD4/c7a7PgeR+XZ+9koSdm55/Ta53cp//JFdWnRX0N7veTkKt0LId0J5s2bp7lz51qXhwwZokGDBlVo21tuucVmOSAggFldDMa3lp8G3/isdfmH1ROVm1f6J0VhYaHO5VwZnxcf09mp9UHKuuC8oSlHMyw3scEY7OnH+qHRevS28dblR24br3p1il/ZhOOYzdKcLeWvVxmXCqRFO52zb9jP3s9Gs9msN78epuPnDigyLFZ/HTK11GcZ1FSEdCd4440rY7F8fHw0fvz4Mta21bx5c5vldu3aydvbu9T1jx07ppdeeklt2rRRSEiIQkJC1KZNG7300ks6fvy4/cWjQm677lHr3ew5FzL045rSbzqbueo97Ty4xjol1SdzX9LpzKPVUmdNtTbNueNVV3KZ3VDs6ccOzROLfF/21U1U3YHTll9snWXjr5a51WEM9vTiV4vHad2un+Tj7au/D/tOIf/bDlcQ0h0sNTVVS5cutS4PGDBADRtW/AaIgADbuUXLGuqycOFCJSQkaNy4cdq+fbuys7OVnZ2t7du3a9y4cUpISNCiRYvs/hlQvgC/2hrYbZR1ecbyd3QpP6/YesfO7Nenc/8mSXphyBfq0DxR5/MyNfH7J6qt1ppo00Hn7n//KSkj17nHQMVVtB9R/TYdcO7+LxVYhtPAGCraiylpS/TZ/H9Isty4zRXmkhHSHSw5Odlmefjw4XZtX3Qcu1T6TaM7duzQoEGDdPbsWXl7e+uZZ57R0qVLtXLlSr3yyisKCAjQmTNnNHDgQO3YscOuGlAxA7qPst78eS7nhOb+8nGxdd5OHqG8/Fzd0PYudW8zUM/c9ZH8fYO0asdMLd+aXGx9VF1evuUpos52qGL3RKGaVKQfUf2qo0/SuZnbUMrrxTNZx/T6F0NUWFigXtcM1e1dH3dFmW7BZDabmcTIgRITE7Vw4UJJlikXz549qzp1Kv6ItMmTJ+uRRx6xLqekpKhdu3ZlHmfatGm6//77bd5ftmyZevXqpcLCQiUmJmr+/PmV+XEqzNnjyKa++KuiwuOcegxHm7Nust5JfkTBgeH6ePROhQVHSpJmrJigST8+o/DgKP139E4FB4Y55fjHzx7QsHGNnbJvI2vQvKvueXmV04+z7vtXtXb6y04/jtG4Yy8WVbQvqutnqam9KJNJT3ycq1q+zn1g25E9K5T82o1OPYZRuVs/FhRc1vMf3axt+5crNrKV3h/1iwJ8g6rt+K7oxarEbM6kO9jmzZut38fHx9sV0CVpw4YN1u/9/PzUunXrYuts2rTJGtBvvvnmYgFdsswU88ADD0iSFixYYPNwJDjf6cwj+mj285Kkx29/xxrQJWlAt1FqcdW1Opt9XB/Ofs5VJXos/+C61XOcIMZPAmXx8Q10ekCX6EV38sncl7Rt/3IF+NXWP4ZPr9aA7o4I6Q6UlZWlM2euXNuzZyz6b1asWGH9vm3btqpVq/gTt2bMmGH9/k9/+lOp+3rooYes30+fPt3uWlB5705/TOfzMtU5vp8SO9kOefLy8tJzd38sH29fzVv/qTbt5b4Bh6qmi4NmcRESMAIGBLiH1dt/0LfL3pIkPXPXZF0V0cLFFRkfz1x1oMxM24GwoaGhdm2/bt067dp15dnjpY1HLxrke/ToUer+unTposDAQOXm5mrlypV21YLKW7z5S63dNVuBfsGlzsEcF9Va9/V6UVMXvKIJyY/oo+e2yd83sJor9Ux5OdUzQPViNR0HcFf5l3J1+dIF1fINKH/lKrh4nl50B299+4AkyaeWn35YNVE/rJpY7jb9Oj+kfl0eKnc9T0VId6ILFy7Ytf6ECRNslkub2WXnTsvEsHXq1CnzbL2Pj4+aNm2qbdu2WbeBc2XknNIHPzwlSXrolnGKCLuq1HWH9HpJK7Ym68CJHZoy7+967Pak6irTo50+tEXmwkKZvJx7ofDkwc3lrwTUZGazTh3aogbNrnPqYehF95BzIUOSlH/5onYcqNh9QzV9mlRuHHWg7OxshYSEWJebN2+u1NSKTai8ePFi3Xyz7aNwN27cWCyoX7x4Uf7+ljF+rVu31vbt28vc76233qo5c+ZYt/X19a1QPUaz8iMpL8vVVbgX/xCp+yPlr+eJ3pgtHXfyDC9jBkqhNfDiB71ov5rci9PXSyuc/FyB+6+XOjdx7jGMin60j7v1ImPSHSg4OFhRUVHW5bS0NO3fv7/c7Y4cOaIHH3zQ5jVfX18lJCQUWzc7+8rTu2rXrl3uvouuU3RbwJN1iHPu/ptG1MyADtjL2b3oW0tKiHbuMQBXIaQ7WPfu3a3fm81mvfrqq2Wun56err59++rQIdvnjCckJJR41rvoEJqKnBX38/MrcVvAk13XVPJ24r9u3ZqXvw4AKa6e1NA5s8xKkjrFSQHueYEYKBch3cF+m/bwN5999plGjx5dLCDn5uZq8uTJateunXbs2KG6dW2njSvtptGiTyS9dKn8ZyFfvHixxG0BTxYSIN1wtXP23TBMalf6rQYAijCZpP5tnbNv31rSzcVnKQY8BiHdwfr3719sxpWkpCQ1aNBA/fv319ChQ9WzZ081aNBAjzzyiM6dO6eoqCiNGzfOZpvSQnpwcLD1+5ycnHLrKbpO0W0BT9e/nVTfwX/lvUzS0Ouce5Ye8DQJ0VInJzw/5vb2Ut3yR30CbouPGgczmUyaNm2amjVrZvN6Zmam5s6dq6+++krLli1TVpblTo8OHTpoyZIlxaZvLG1mFz8/P9WrV0+SdPjw4XLrSU9PlyRFRES47U2jQGX41pL+2F3y93HcPgd1kqJ5bgpgt0GdpIahjttf+6ukbk66WgYYBSHdCaKjo7Vq1SoNGzZMXqVMAxcfH6+kpCStW7dOLVq00MaNG63v+fj4qG3b0q8PtmrVSpIl+B89erTU9fLz87Vv3z6bbYCaJDpceqyXY8asDuwodScUAJUS6Cs9frNjxqe3v0r6Q1fLlS3AkzFPupNERETo888/V1JSkhYtWmQ9692oUSO1bNlS7du3t1m/aEhv3bq1zQ2fv3fDDTdo+fLlkqRly5ZpyJAhJa63fv165ebmSrK9oRWoSeLqSc/2k75aI+0/Zf/2oYHSvddKLe1/gDCAIoL9pVGJ0o+bpdV77d/ex9syjK1HvOTkxyAAhkBId7L69evrvvvuK3OdrKwspaWlWZdLG+rym0GDBmns2LGSpI8//rjUkP7JJ59Yvx88eHBFSwY8Tv1gaWSitCZNWrpLOlWB2UgDfaVrm0p9Epg9AnAUfx/pni6Ws+ELtkt7T5S/jZfJcrN2v7ZSZEj56wOegpBuAJs2bVLRZ0qVdtPobzp06KDevXtr4cKFWrRokb788ksNHTrUZp3ly5drypQpkqTExMRiZ+6BmsbLZJk6sWszSzDYdVQ6fLbkkPCHrlLbGMu4dgCOd3WU5et4ppRyyNKL20u4zWpAR6lDrGXGJqCm4YKRARQd6iKVH9IlacKECdbZWoYPH67nnntOy5cv1+rVq/Xqq6/qlltuUUFBgYKDg/XOO+84pW5PcfzsASU+b9K/vn6gwtvMWz9Fic+bNG/9FJvXE5832XydzTpepdr2pG+w2d8fXo+r0v5gmRLu6ijpzg7SE6U8cbpTYwK6Kxi5F5dvTbbZ33OTelZpf7CIqiP1ayM93KPk93u2IKC7gpF78bul4232Z0+N7oaPIQMoGtK9vb3LvGn0N61bt9aMGTN077336uzZs3r77bf19ttv26xTt25dffPNN2rduuZNJJv4vO0dRV4mL9UOCFOTBm11y7UPq9c1Q0vZsuoiw2LVp9MDkqQAv+Lzg53KOKzP5v9D6/f8rOzzZxQe0kBdWw/QsMSXFRxoe1dV3ZCGGpb4siRpxooJTqsZcBYj9mLW+TNauf17/bLrJ/16fJtOZx5RrVq+ahzVRn07P6i+nR4sdtN/bGQray9OXfCK02oGnMWIvViShRun6c2vh0mSnrlrsvpf+7DN+63iumpY4svKuZCh71e+65R6jYKQbgBFQ3qrVq0q/NCh3r17a/v27Xrvvfc0a9YsHTx4UJIUGxur22+/XU899ZSioqKcUrO7+O1D9XJBvtJP7tbqnT8oZd8SpaZv0GN3vF3O1qXrnjBQLWOvU93gBsXeiwyL0/A+Y0rc7ujpfXrq312VkXNSXVvfqZiIFtpz6Bd9v/JdbdjzsyY8sUohQVcebFWvTkPrvuZvmFLpegFXM1IvLtv6nd6b8bjCQxqofdObFBF6lc7lnNDKbTP09ncPa/3uufr7sO9kMl0JNbGRraz7IqTDnRmpF3/vZEa6Js4cqQC/2rpwseRnwbSO66rWcV11/OwBQjqcKysrS3v3XrnNvSJDXYpq0KCBxo0bV+xhSLD4/T8Km/Yu0guTEzVj5QQN6D5KUeFxldpvUEAdBQXUsXu7977/szJyTuqJO9/TgO5PWl//z4/PavqKd/TJz3/T04P/U6maACMzUi9G179arz74o65tcavNGfOH+r2uke930Ypt07Vy2wzd0JYb7uF5jNSLRZnNZo3/5kGFBNZV9zaD9N2y8ZXel6dgTLqLbd682eam0fJmdkHVdGh+s2Lqt5DZbNae9PXF3j9+9oDGTrtPg1+up/4v+uvP73bS2p2zi61X2ti7shw9vU8bU+crKixOd3R9wua94X1ekb9vkBZtnKoLl87b/XMB7saVvXhNs166vtXtxYa0hIdE6bbrHpMkbdm/1K6fB3BXruzFor5f+Z5S9i3W6Hs/lb9vUKX24WkI6S7Wo0cPmc1m69eTTz5Z/kaoErMsvxQVvZQtSSczDurJ97ro+LkD6t1xmHq2u1cHjm/Xy1PuVErakiofN2WfZR8dr+5TLBwE+gerdVw35eXnatfBtVU+FuAOXNWLZanlbXlErbcXF5pRc7i6Fw+e2KWP57yggd2fUtsmNzpsv+6Of4VQo2xKXajDp/bIZDIpPqazzXtb9i3V8MQxGtbnZetrN10zVC/9t5++XfaW2je7qUrHPnxqjySpUf2SH1vZqF5zbUydryOnUtWh+c1VOhZgdK7sxdIUFFzWgo2fS5I6xfdzyjEAo3F1LxYUXNabXw9TRNhVeuiW16u8P09CSIdH+3z+GEmWG2QOn9qjVTtmymw2a/ANzygyLNZm3ciwWA3t/X82r3WO76uI0Ku059AvVa7lfF6mJCnIv+Qxe7+9npOXUeVjAUZjpF4szX/nvKADx7erS4v+6hzf12nHAVzJaL04deGr2ndks97580r5+TDfZlGEdHi032ZhMJlMqu0fqjaNb1C/zn9S745/KLZu04bt5e3lXez1+qEx2nVwjdNrBTyZ0Xvx+5XvKXl5kmIiWuivQ6Y65RiAERipF3cdWqevFr+uu258Tq3irq/y/jwNIR0ebcFb5vJX+p8g/9ASX/f2qqVCc2GVa/ntTPlvZ9R/77fXa5dSB+DOjNSLvzdz1UR98MNTio1spX89skghgeEOPwZgFEbpxYKCy/rX18MVXe9q/bHfa1Xal6cipAPVJLp+vCTpyKnUEt8/ctoyFWdpY9YBON6MFRM06cdnFBeVoH89ukhhtSNcXRJQI1y4lKPD//s8vPVF/xLXeSd5hN5JHqGB3Z/Sn++cUI3VGQMhHagm7ZtabrDZmDpfhYWFNjO85OZla8eBVfL3CVTL2OtcVSJQo3y95E19POcFNW3YXm8+skB1guq5uiSgxvCp5ad+Xf5U4ntpRzYp7chmJTTuruj68WoVWzOHwhDSgWrSsF5Tdby6jzamztePq/9t8zCjz+e/rLxL53XrdY8qgPlhAaebtuA1fTb/H2oe3VFvjJjPEBegmvn5BOi5u/9b4nufzx+jtCObldjxj+p/7cPVXJlxENKBajRq4Ad66t9d9e8fRmlz2iJdFdFSuw+tU8q+JYquf7Ue6jfW1SUCHm/+hs/02fx/yMvLW20a36CZK98rtk5kWJz6dn6g+osDgP8hpAPVqGG9pvr3Uxv02bx/aMOen/XL7jkKD26ggd2f0rDElxUcGObqEgGPd/zsr5KkwsICzVgxocR12jbpQUgH4FKEdHgke+5ejwqPK3P9pMeXFnutb+cHKv0BHhEao+fv/bRS2wLuxoi9OLzPGA3vM8aubQB3Z8ReLA09akFIBxxs6/5lSnze8mjlb/5+TOEhUZXe1570DRr53pUnwP3+QRMASufIXly+NVmvTb3bUaUBNYoje/G7peP10U/PO6o0QyOkAw40LPFlm+UAv9pV2l/dkIY2+6wdEFql/QE1haN7MTaylc0+I8PiqrQ/oKZwdC+2iutqs8+mDdtXaX9GRkgHHMjRl+fq1WnIJT+gEhzdN7GRrehFoBIc3Tet47qqdVxXh+7TqLzKXwUAAABAdSKkAwAAAAZDSAcAAAAMhpAOAAAAGAwhHQAAADAYQjoAAABgMEzBCLfgF+zqCtwPf2ZwBv5e2Y8/MzgLf7fs425/XoR0uIXOQ1xdAQCJXgSMhH70bAx3AQAAAAyGkA4AAAAYDCEdAAAAMBhCOgAAAGAwhHQAAADAYAjpAAAAgMEQ0gEAAACDIaQDAAAABkNIBwAAAAyGkA4AAAAYDCEdAAAAMBhCOgAAAGAwhHQAAADAYAjpAAAAgMEQ0gEAAACDIaQDAAAABkNIBwAAAAyGkA4AAAAYDCEdAAAAMBhCOgAAAGAwhHQAAADAYAjpAAAAgMEQ0gEAAACDIaQDAAAABkNIBwAAAAyGkA4AAAAYDCEdAAAAMBhCOgAAAGAwhHQAAADAYAjpAAAAgMEQ0gEAAACDIaQDAAAABkNIBwAAAAyGkA4AAAAYDCEdAAAAMJhari4AqIj1X0kXs11dhXvxC5Y6D3F1FfA09KL96EU4C/1oH3frRUI63MLFbCkvy9VVAKAXAeOgHz0bw10AAAAAgyGkAwAAAAZDSAcAAAAMhpAOAAAAGAwhHQAAADAYQjoAAABgMIR0AAAAwGAI6QAAAIDB8DAjADVGQaF0IlNKPysdPitl55W83rLdUky41ChM8vOp3hqBmsBsls6et/Ri+hnL9yWZv02KDpdi6krB/tVbI+BqhHQAHu3wWemX/dLBM9LRc1J+QfnbfL/R8l+TpMg6lpDQ/iqpVUPJi+uPQKVk5Eq/7JP2nbSE89xL5W8zZ+uV70MDLb88xzeQOjWW/PkFGh6OkI4a77lJPbV1/7IS3/MyeSnQP0ShtSPUrFEHdW19p7onDJJPLd9qrhL2yC+QthySVqRKB09Xfj9mScczLV8bfpXCAqWuzaXrmnFWz1noR89iNkt7T0irUqVth6VCc+X3lZFr+dp2WJq1WercWOp2tdQg1GHlogh60fUI6UAZCs2FyrmQoZwLGTp8KlVLU75Ww3rN9NLQrxQf08nV5eF38vKlhTukNWnS+YuO3/+5XOmnLdLP26R2MVK/tlJEiOOPg5LRj+6j0Cyt2yct2SWdzHL8/i9ellbutXw1jZASE6QWDRx/HJSMXqwehHSgiPiYzoqP6WJdNpsLlXn+tHYeXK3TmUckSUdPp+mvH/XWuyPXKDaypatKxe/sOSZ9vdYSpJ2toFDadNByRu/WdtKN8QyDcQb60T2dzpa+WmsZ1lId9p2U9i2Wrm0q3dlBCuRkrsPRi65BSAeK6NKiv4b3GVPs9YLCAs1aM0mTfnhKheZCnc/L1KQfn9YbI+ZVf5GwkZcv/bhJWp1W/cfOL5BmbpK2pEtDruOsuqPRj+6l0CytTJVmb5YuVeDeD0dbt0/afVS67zqpZcPqP74noxddg3M/QAV4e3lrQLeRuqfnX6yvbUydr5MZ6S6sCvtPSm/+5JqAXtSvp6S35kgr9ri2jpqCfjSejFzp3wulGRtcE9B/k3lB+nCJ5apaRW4SR9XQi85FSAfscEfXJ2yWt/+60kWVYMdhadJi6VwpU7dVt/wCafoG6cfNlpvl4Hz0ozGczpbem199w1sqYu0+6aMllittcD560TkI6YAd6odGKzggzLp8NuuYC6upuTYdkD5ebswzZYt3St/9QlCvDvSj6x3PtAT00uY5d6W9J6QPFlVsqkdUDb3oHIR0wE6+PgHW7y9dLuVpOHCaHYelaaurNpWbs61Os4xVJ6g7H/3oOmdypEmLpCwD/7EfOiNNXmqZDQbORS86HiEdsMPlgnxl5Z6xLofWru/Camqe/SelKSuNHdB/s2y3ZTpIOA/96Do5eZaAnnnB1ZWU79dT0qfLpcJCV1fiuehF5yCkA3bY9usK5V++MgF3i5hrXVhNzZKXL01dbcwhLqWZs7VqD1NC2ehH10leL53OcXUVFbf7mLSMG7udhl50DqZgBCroUn6ePpnzonX56uhOatKwrQsrqllmbXbOTaIT7i/+2tNfOGbfZrP05RppdH/Jx9sx+4QF/eg6KQellEOO368ze1GS5myRWjWSIpkq1aHoRefhTDpQhsLCQmXknNKKrdM1auL12p3+iyQpJLCuRt/7qYurqzlSj0ur9rq6iso5kSX9vNXVVXgG+tH1svOk79a7uorKyS+QvlrDsBdHoBerB2fSgSKmLnhFUxe8UuY6rWKv15MDP1DjqIRqqqpmy8u3PL3QnS3eJbWNkWLruboS90I/Gk/yeun8xfLXM6oDp6Wlu6VerVxdiXuhF12DkA7YITaylV4YMk0N6jYpdZ2t+5dr9H96ysvkrYmj1qtZo/YlrrfjwGo9+8ENksmk90auVXxMJydV7d4W7TDOXOiVZTZbws2z/SSTydXVeI7y+nHstPu0dMs3at/0Jv3r0UUylfKHbzab9ZcPb1bKviXq2f4+/e3+r5xZtttKPS5tccIwl+o2d6vUuYkU7O/qSjxHRT4bn5vUU1v3L7Nrv+MfW6J2TXtWsTr3RUh3gYKCAq1evVp79+7VyZMnFR4erpiYGN14440KCgpydXk1WnxMZ8XHdLEun7+QoSNn0rT70DpJ0sETOzXq/es0/vGlio0s+VRM2yY36rbrHtOsNZOU9N1DmvjkL/L2tm21S/l5SvruTyo0F+qeHn8hoJficoG0xsVPE3WU9LOW6eA4m15xVe3HJwa8r81pi5Syb4lmr/1Qt1//WInH+WntR0rZt0ShQfU1csD7zvlhPICnPFE3v0Bat0/q3drVlbgPR3w2Nm7QRgWF5c+FeeRUqjLOn5IkBfmHVr14N0ZIr0YZGRkaO3asPvvsM506darY+yEhIRoxYoTGjRsnHx8fjRkzRq+8cuXy0o4dO9SqFdfonKlLi/4a3mdMsdf3H92qcV8O1YETO5Rx/pRe+XywJj29SX5F5oUt6uH+b2rtrtlKO7JZ3y57S0N6vWjz/ucLxij95G5F179af+xT9iXEmmzLISnHjS+t/97KvYR0e1S1H0Nr19fIARM19ov79N+f/qJrW96qiNAYm3VOZRzW5DmWR5qPHDhRdYL4H1SSc+el7UdcXYXjrNor9WopeXFnXoU44rOxIr8AHzq5WyPf6yxJGnzjs6Veia4p+OtZTebMmaPmzZtr/PjxJQZ0ScrKylJSUpL69++vixcvasuWLdb3goKC1KJFi+oqF7/TpGFbjX14roIDwyVJ6Sd3a9rC10pdP9A/WE8P/lCSZSzfoZO7re+lHt6o75aNl8lk0rN3/Ve+PlxzLc1KN71ZtDSbD7j3eF6jsKcfe7a/V91aD1DuxWy9892IYu9PmP6ocvOy1C1hoHq0u8epdbuzNWme9XCuc+elXUddXYX7s/ezsSwXLp3Xq58P1oWLOWod100j+r/pyFLdEiG9GnzyySe64447dPr0lQmTvby81KlTJw0ePFiDBg1S06ZNre8tXLhQH3zwgVJSUqyvtWvXTl78yu9SEaExerDfWOvyzBXv6lz2iVLX79LiFvXuMEz5ly8q6ds/qbCwUJcL8jX+mwdVWFig26//s9o0uaE6SndLR85ZHkLiSS4XWi6zo+rs6cdRgyYpOCBMG1Ln6ef1V2aeWLBxqn7ZPUfBgeEaNfADp9fsrjxp2FlRnnYSwFXs/WwszbvJj+rgiZ0KrR2h//vDt8WGidZEpD4nmz17tkaMGKGCAssTWLy8vPTss8/q0KFDWr9+vZKTkzV9+nSlpaVp1qxZCgsLkyS99tprOnDggHU/HTt2dEX5+J3+XR5WVHhjSVJefq6+XvJGmes/fucEhdaO0M6DqzVz1fv6ctFY/Xp8myLDYvVw/7K3renW73d1Bc7xi4f+XK5Q0X4MD4nS43dMkCR9OOtZnc48qnPZJzTpx6clSY/fMUHhIVHVUbJb2nPMMvWip9l9VMpygyemugN7Pxt/b9bqSVq0+Qt5mbz04tAvVa9OQ2eU6XYI6U6Unp6uYcOGqfB/k7L6+/tr3rx5SkpKUqNGjYqtf9ttt+n777+XyWTSuXPnbN4jpBuDt3ct3dvzr9blOesmK+v8mVLXDwkM18gBEyVJn859SV8tfl2S9NTgDxXgV9u5xbq5g6X/sbq1E5mWaSVRdfb0Y2Kn4erSor9yLmTo3RmP6f2ZI5Wde1bXtrxViR2HVVfJbslTe9EsKd1Df7bqZu9nY1Gphzdq0qxnJEnD+76qDs1vdkqN7oiQ7kSPPvqoMjIyrMtffPGFevfuXeY2PXr0UN++fYu93qFDB0eXh0rq0/kB1Q2x/Jafd+m8ZqyYUOb6PdrdrW6tBygvP1eXC/LVp9MD6hxf/P8xrigslI6cdXUVzmGWZSgPHMOefnx68IcK9A/R2p2ztGJrsoL861jvHUHp0j20FyXP/tmqm72fjZKUnXtOr31+l/IvX1SXFv01tNdLTq7SvRDSnWTevHmaO3eudXnIkCEaNGhQhba95ZZbbJYDAgKY1cVAfGv5afCNz1qXf1g9Ubl52aWuX1hYqHM5V8bnxcd0dmp9nuBktnSpwNVVOA9n7xzHnn6sHxqtR28bb11+5Lbxqlen+FVNXGE2e3aQ9eSfrbrZ+9loNpv15tfDdPzcAUWGxeqvQ6aW+iyDmoqQ7iRvvHFlPJaPj4/Gjx9fxtq2mjdvbrPcrl07eXt7F1svIyNDixYt0htvvKG77rpLcXFxMplMMplMiouLq3TtKN9t1z1qvZs950KGflxT+k1nM1e9p50H11inpPpk7ks6ncm0AmU55OEhlmDgWPb0Y4fmiUW+L/vKJqSMXCnHA8ej/4ZedCx7evGrxeO0btdP8vH21d+HfaeQ/22HKwjpTpCamqqlS5dalwcMGKCGDSt+E0RAgO38oqUNdbnmmmvUu3dvvfjii5o+fboOHjxYqXphvwC/2hrYbZR1ecbyd3Qpv/gn2bEz+/Xp3L9Jkl4Y8oU6NE/U+bxMTfz+iWqr1R0d9vAPTk//+apbRfsR9vP0v6tZF7h51JEq2ospaUv02fx/SLLcuM0V5pIR0p0gOTnZZnn48OF2bV90HLtU+k2j5iKT1tatW1d9+vRR7drcjFhdBnQfZb3581zOCc395eNi67ydPEJ5+bm6oe1d6t5moJ656yP5+wZp1Y6ZWr41udj6sPD0ucQ96QFNRlGRfoT9asLfVU++UuAK5fXimaxjev2LISosLFCva4bq9q6Pu6JMt8AklE6wZMkS6/deXl664Qb75sL+/cOOSgvpI0eOVGxsrDp16qTGjS1TH8XFxSknJ8fOiqvO2ePIpr74q6LC45yy76THl1Zqu+DAMP34z9LH281ZN1kpaYsVHBiuJ/83w0tUeJwe7PdPTfrxGf175pO6ptnNCg4Mq9Txy3PgwAHdYGrslH07262jktWsy+AKrz/hfscev7L7e/qLiq13LjNHJlNw5Q7iYs7sRcl5/ehK7tyLbROf0E1/nFjh9Y3Si1LF+7FDp846sX9D5Q/kQu722VhQcFljv7hP53JOKDaylZ6+66MqVGg/V/SiuQpPAeNMuhNs3rzZ+n18fLzq1Klj1/YbNlz5x8LPz0+tW7cucb3Ro0fr7rvvtgZ0GMfpzCP6aPbzkqTHb39HYcGR1vcGdBulFlddq7PZx/Xh7OdcVaKhefv4uboEp6rFU2bhJrxreXYvSpI3/VhtPpn7krbtX64Av9r6x/DpCvANcnVJhkZId7CsrCydOXPlrjd7xqL/ZsWKFdbv27Ztq1q1uODhbt6d/pjO52Wqc3w/JXayHe7k5eWl5+7+WD7evpq3/lNt2rvIRVUaV0GBZ08kXnD5kqtLACqk0MN7UaIfq8vq7T/o22VvSZKeuWuyropo4eKKjI/052CZmZk2y6GhoXZtv27dOu3atcu6zEOM3M/izV9q7a7ZCvQLLnUO5rio1rqv14uauuAVTUh+RB89t03+voHVXKlxFXj4TX+XL3GnGtyDp/eiVDN+RiN469sHJEk+tfz0w6qJ+mFV+cOo+nV+SP26POTkyoyLkO5kFy7Y92E8YcIEm2UeYuReMnJO6YMfnpIkPXTLOEWEXVXqukN6vaQVW5N14MQOTZn3dz12e1J1lWl4ly5kuboEp/L0nw+eoyb8Xa0JP6MR5FzIkCTlX76oHQdWVWibmj5NqslclRHtKCY7O1shISHW5ebNmys1NbVC2y5evFg332z7ONyNGzfaFdTj4uJ08OBBxcbG6sCBAxXezuhWfiTl8e+oXfxDpO6PuLqKylmZKiWvr55jlXRjWkVvOKushGjp4R7OPYaz0Iv2c+dePJ4pvTG7eo7lil7095Fev1vyctNn6NCP9nG3XmRMuoMFBwcrKirKupyWlqb9+/eXu92RI0f04IMP2rzm6+urhIQEh9cIGF20hz/TIsbDfz54johgydeDr7lHh7tvQIfnI6Q7Qffu3a3fm81mvfrqq2Wun56err59++rQoUM2ryckJMjX19cpNQJG1jDUsz84CelwF15eUrRzZok1BHoRRkZId4IHHnjAZvmzzz7T6NGji41Pz83N1eTJk9WuXTvt2LFDdevWtXmfm0ZRU/nWkqLsm7nUrRAM4E48+cqWJ/9scH+EdCfo37+/evSwHXCalJSkBg0aqH///ho6dKh69uypBg0a6JFHHtG5c+cUFRWlcePG2WxDSEdNFlO3/HXcUZ0AKTjA1VUAFXeVh/aixC/MMDZCuhOYTCZNmzZNzZo1s3k9MzNTc+fO1VdffaVly5YpK8tyt0eHDh20ZMmSYtM3MrMLarL2pU+M49bax7q6AsA+rRpKPt6ursLxGoZJ9d3zwb+oIQjpThIdHa1Vq1Zp2LBh8vIq+Y85Pj5eSUlJWrdunVq0aKGNGzda3/Px8VHbtm2rq1zAcOIbSHVru7oKx+vW3NUVAPYJ9JM6xLm6Csfr3lwyefC9L3B/HnzPtutFRETo888/V1JSkhYtWqTDhw9Lkho1aqSWLVuqffv2NusXDemtW7eWn5/nP44ZKI2XyRJof9zs6kocJz5Kiggpfz3AaLo3l9btc3UVjuPvI3WMc3UVQNkI6dWgfv36uu+++8pcJysrS2lpadZlhroA0rVNpTlbpMuFrq7EMbpd7eoKgMqJqWsZm37ojKsrcYzOTSQ/H1dXAZSN4S4GsWnTJhV9rhQ3jQJSkJ90TZyrq3CM0ECpdSNXVwFUXncP+iWTYWdwB5xJN4iiQ12kioX0lJQUpaSk2LyWk5Nj/e+UKVNs3uvevXuxm1khHT97QMPGNVZixz/qL/dNqdA289ZP0fhvH9Toez5V384PWF9PfN52gOM3fz+m8JAoVdbyrcl6berd1uW2TXoo6fGlld6fO+qbIG05JF267OpKqua29pI3p0XKRC8aW4dYaelu6eg5V1dSNV2aePYUr45ALxoDId0gioZ0b2/vCt00OnPmTL3yyislvnfmzJliTzD99NNPa0xI//0/Cl4mL9UOCFOTBm11y7UPq9c1Q5127MiwWPXp9IAkKcDP9s7H5VuTtXX/Mu07mqL9R7co92K2br7mfr0wdFqJ+4qNbKVhiS9LkqYuKPn/taerFyzd3l6avsHVlVReQnTNHf9q1F6c/NNflXp4g46cSlXm+dPy8wlQRFisurUeoDu7jVRIkO28g/SiVMtbGnqd9PbPUqG5/PWNqE6ANKCGXqg2ai/+4fU4nTh3sMTtwmpH6tuXj9u8VpN6kZBuEEVDeqtWrRQQwETKjvBbI18uyFf6yd1avfMHpexbotT0DXrsjrcrvd/uCQPVMvY61Q1uUOy9yLA4De8zpsTtvlj4T+0/tkUBfrVVr060ck/uLvM4sZGtrPvy9H+MytLtamlLupR2wtWV2C/QV7q7C7NIGK0XZ6x4R80adVCHqxMVWjtCeZfOa9fBtfp8wRj9tO4jvffkWkWExljXpxctosOlPgnSz9tcXUnl3HutpSdrMqP1oiQF+dfRoBueLvZ6gG/xKb5qUi8S0g0gKytLe/futS5XdDz6mDFjNGbMGCdV5Rl+/4/Cpr2L9MLkRM1YOUEDuo9SVHhcpfYbFFBHQQH2Xy99/I53VK9OtBrVa6at+5dp9H9uqtTxaxovk3TftdK/5rjfsJdBnSxn72o6o/XiD69lydfHv9jrn8z9m75a/Lq+XjxOowZ9UKmaPF3v1tLWw+437KVLE6kV94UYrhclqXZAaJkhvqZihKQBbN682eamUWZ2cZ4OzW9WTP0WMpvN2pO+vtj7x88e0Nhp92nwy/XU/0V//fndTlq7c3ax9eatn6LE502at36KXcdv3+wmRddvLlNNP61aCfWCpYFOuEz99BfFvxylXUzNHeZSHlf3YkkBXZJ6tLtHknTk9N4S34dl2Mv91zv+AUfO7MW6tWvuMJfyuLoXUTpCugH06NFDZrPZ+vXkk0+6uiSPZpblF6LfB+WTGQf15HtddPzcAfXuOEw9292rA8e36+UpdyolbYkrSsXvXN/McqndHTSNkO7vyjCXshixF9funCVJatyAh8mVpVGY9NCN7nEzdIi/9HgvhrmUxdW9mH/5ohZunKYvF72uGSveVUraEhUUFjhs/+6K4S6oUTalLtThU3tkMpkUH9PZ5r0t+5ZqeOIYDevzsvW1m64Zqpf+20/fLntL7ZsxNMUIbmkr5eVLy/e4upLSXVVXGtFT8uVf2FIZpRe/WzpeFy7l6HxeplIPb9D2X1eqSYO2uu+mFxx2DE/VsqE0rJv0+Urj3kga5Cc91styJQ4lM0Ivns0+rje/HmbzWlR4Y42+51O1a9rDIcdwR3yEwKN9Pn+MJMsNModP7dGqHTNlNps1+IZnFBkWa7NuZFishvb+P5vXOsf3VUToVdpz6JfqKhnlMJksw158a0kLd7i6muKaRlgCuj8PSrFh1F78btl4ncu5ckdy5/h+ev7eKQqtXd+hx/FU7a+Sat0oTVlhvIeO1QmQHr+Z6RZ/z2i92LfTg0pocoPiIlsrwC9Yx87u1w+rJmrOuo/0t49v0bsj16hpw3YOOZa7IaTDo/1257fJZFJt/1C1aXyD+nX+k3p3/EOxdZs2bC9vr+KDLOuHxmjXwTVOrxUVZzJZ5h0P8pNmbTbOWbx2MZYhLpxBL86ovfjb9G7nsk9ox4HV+njuC3pswjX654Oz1Tya+4MqIiFaevQmacpK6fxFV1djEVXH8sty3eKTg9R4RuvFomfpJalxVIKeHvwfBfjWVvLyJE2dP0ZjHvjeIcdyN3yUwKMteKvi6S3IP7TE1729aqnQbLBTRJAk3dRSah4pfblGOprhujoCfS2zuHSMYwx6aYzei2HBkereZqCaR3fQg29erX99PVyTR293yrE8UfMo6YXbpORfLNOluorJJPVqKfVr6/gbWz2F0XvxN7dd/5iSlydp66/LnXocIyOkA3Br0eHSs/2kBTukBdur/6x6QrR0TxcphGkWPUJkWKyuimylfUdTlHn+tOoE1XN1SW4j2F968EZp80EpeX31n1WPqiMNvd5yTwjcX2iQZchZ3qXzLq7EdQjpANxeLW/LDaVtoqXp66VfTzv/mGGB0q3tOXvuic5kHZUkeZk4FVsZ18RKzSKlHzZKmw46/xdnv1rSjfFSnzacPfckOw+tlSQ1CG/i4kpch5AOwGNEh0tP9ZUOn5VW7ZU2/ipdcvAsXi0aSN2vllo1lLzcYPo5FHf4VKrCakcWe/BKYWGhpsz7uzJyTqpVbFcFB4a5qEL3F+wv/aGbdNs10po0y1fWBcceo0Edy9OIOzXmRm13dfDELkWEXaUA3yCb14+fPaCJM0dKkm7uUHysfE1BSAeq0artM7Vq+0xJ0rlsyw1rOw+u0b++fkCSVCeonh69fbyLqvMc0eGWx3/ffo20fr/l62hG5c/ohQVK7WOlrs2l+kzl5vZ+2T1HH895UQmNuysqvLFCAuvqXM4Jbd23TMfO7ld4cJSevXuyq8v0CKGBlqtcfRKkbYelNXulX09V/pfnQF8p/n+/KDepz1Usd7dsyzdKXp6kNo1vVGRYrGV2lzP7tG7XT7p0OU9dWvTX3T1Gu7pMlyGkA9Vo39EULdj4mc1rx87u17Gz+yVZxsMS0h0n0Ffq0cLydemydOSc5Sx7+lnLf3PyLGHhcoHloSw+3pYzcg3DpJhwS9iPCZdql/xwSripa5r3Vr8uadr+60qlHdmsnLwM+fsGKbre1erdcZgGdB+lkMBwV5fpUby9LNM1tr9KKiiUTmZZ+jD9rHT4jHT2vJRfYPkySfKpZenHyJD/9WFdSy+GBxHMPUm7pjcp/dQe7TuyWTsOrFLepfOqHRCqhMbd1bvDMPXuOKxGP6HbZC76PHrAoFZ+JOVlubqK8iU+b1LbJj2U9PhSl+/bP0Tq/ojDy0ANRy/SizAOd+hHerHyOJMOONjW/cuU+LzlN/9v/n5M4SFRld7X8q3Jem3q3Y4qDahR6EXAGOjFyiGkAw40LNH2oQwBflV7kkZsZCubfUaGxVVpf0BNQS8CxkAvVh7DXeAW3OGSntG422U9uAd60X70IpyFfrSPu/UiE4gBAAAABkNIBwAAAAyGkA4AAAAYDCEdAAAAMBhCOgAAAGAwhHQAAADAYJgnHW7BL9jVFbgf/szgDPy9sh9/ZnAW/m7Zx93+vJgnHQAAADAYhrsAAAAABkNIBwAAAAyGkA4AAAAYDCEdAAAAMBhCOgAAAGAwhHQAAADAYAjpAAAAgMEQ0gEAAACDIaQDAAAABkNIBwAAAAyGkA4AAAAYDCEdAAAAMBhCOgAAAGAwhHQAAADAYAjpAAAAgMEQ0gEAAACDIaQDAAAABkNIBwAAAAyGkA4AAAAYDCEdAAAAMBhCOgAAAGAwhHQAAADAYAjpAAAAgMEQ0gEAAACDIaQDAAAABkNIBwAAAAyGkA4AAAAYDCEdAAAAMBhCOgAAAGAwhHQAAADAYAjpAAAAgMEQ0gEAAACDIaQDAAAABkNIBwAAAAyGkA4AAAAYzP8DRRW0NRJsABIAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 968.63x361.2 with 1 Axes>"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "phi = ParameterVector('Phi',6)\n",
    "\n",
    "qc = QuantumCircuit(2)\n",
    "qc.rx(phi[0],0)\n",
    "qc.rx(phi[1],1)\n",
    "qc.cnot(0,1)\n",
    "qc.ry(phi[2],0)\n",
    "qc.ry(phi[3],1)\n",
    "qc.cnot(0,1)\n",
    "qc.rz(phi[4],0)\n",
    "qc.rz(phi[5],1)\n",
    "\n",
    "qc.draw('mpl', scale=2.5)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In principle, your circuit can have the quantum gates that you want as long as has parameters and also a varierity of them. Depending of the problem itself the number of the parameters can increased, but careful to use to many parameters which this will be the quantum circuit more complex and slowly to evaluate expectation values.  \n",
    "The next function will calculate the expectation value given a initial parameters and operators for a given quantum circuit. I will explain each line in detail."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def ShiftRuleV3(qc, initial_values, operators : list):\n",
    "\n",
    "    num_para = qc.num_parameters\n",
    "    num_qubits = qc.num_qubits\n",
    "\n",
    "    operators_list = []\n",
    "\n",
    "    ### Convert the matrices into operators using qiskit functions.\n",
    "    for k in range(len(operators)):\n",
    "        Opt = Operator(operators[k])\n",
    "        operators_list.append(CircuitOp(Opt.to_instruction()))\n",
    "\n",
    "    ### Expectation Value Calculation ###\n",
    "    # First we assign the parameters and calculate the expectation value.\n",
    "    qc_normal = qc.bind_parameters({phi:initial_values})\n",
    "    psi_normal = CircuitStateFn(qc_normal)\n",
    "\n",
    "    exp_value = 0.0\n",
    "    # Here we iterate over the different operators.\n",
    "    for k in range(len(operators)):\n",
    "        exp_value += psi_normal.adjoint().compose(operators_list[k]).compose(psi_normal).eval().real\n",
    "\n",
    "    gradients = []\n",
    "    ### Gradient Calculation ###\n",
    "    for i in range(num_para):\n",
    "\n",
    "        exp_value_plus = 0.0\n",
    "        exp_value_minus = 0.0\n",
    "\n",
    "        actual_parameter = initial_values\n",
    "        actual_parameter[i] = actual_parameter[i] + np.pi/2\n",
    "\n",
    "        qc_plus = qc.bind_parameters({phi:actual_parameter})\n",
    "\n",
    "        actual_parameter[i] = actual_parameter[i] - np.pi\n",
    "\n",
    "        qc_minus = qc.bind_parameters({phi:actual_parameter})\n",
    "\n",
    "        psi_plus = CircuitStateFn(qc_plus)\n",
    "        psi_minus = CircuitStateFn(qc_minus)\n",
    "\n",
    "        for k in range(len(operators)):\n",
    "            exp_value_plus += psi_plus.adjoint().compose(operators_list[k]).compose(psi_plus).eval().real\n",
    "            exp_value_minus += psi_minus.adjoint().compose(operators_list[k]).compose(psi_minus).eval().real\n",
    "\n",
    "        gradients.append(0.5*(exp_value_plus-exp_value_minus))\n",
    "\n",
    "        actual_parameter[i] = actual_parameter[i] + np.pi/2\n",
    "\n",
    "    return exp_value, np.array(gradients)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Solving the problem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Eigenvalues: [ 0. -2.  0.  0.]\n",
      "Initial Parameters: [0.37454012 0.95071431 0.67454012 0.45071431 0.25071431 0.75071431]\n",
      "Initial Expectation Value: -0.09125348550048962\n",
      "Initial Gradient Values: [-0.05948042 -0.07455684  0.06623967 -0.17180918  0.18112438 -0.18112438]\n"
     ]
    }
   ],
   "source": [
    "initial = np.array([0.37454012, 0.95071431, 0.67454012, 0.45071431, 0.25071431, 0.75071431])\n",
    "#initial = np.array([0.5,0.5,0.5,0.5,0.5,0.5])\n",
    "#initial = np.random.random(6)\n",
    "\n",
    "I = np.array([[1,0],[0,1]])\n",
    "X = np.array([[0,1],[1,0]])\n",
    "Y = np.array([[0, -1j],[1j, 0]])\n",
    "Z = np.array([[1,0],[0,-1]])\n",
    "\n",
    "#U = np.kron(np.kron(Z,I),Y)\n",
    "\n",
    "U1 = np.array([[1,0,0,0],\n",
    "     [0,-1,0,0],\n",
    "     [0,0,-1,0],\n",
    "     [0,0,0,1]])\n",
    "\n",
    "U2 = np.array([[-1,0,0,0],\n",
    "     [0,0,1,0],\n",
    "     [0,1,0,0],\n",
    "     [0,0,0,-1]])\n",
    "\n",
    "R1 = np.kron(I,I)\n",
    "R2 = np.kron(X,X)\n",
    "R3 = np.kron(Y,Y)\n",
    "R4 = np.kron(Z,Z)\n",
    "\n",
    "eigen_val = np.linalg.eigvals(U1+U2)\n",
    "print(\"Eigenvalues: {}\".format(eigen_val))\n",
    "\n",
    "H = [U1, U2]\n",
    "#H = [-R1, R2, R3, R4]\n",
    "\n",
    "exp_value, gradient = ShiftRuleV3(qc, initial,H)\n",
    "print(\"Initial Parameters: {}\".format(initial))\n",
    "print(\"Initial Expectation Value: {}\".format(exp_value))\n",
    "print(\"Initial Gradient Values: {}\".format(gradient))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Optimizer(qc, initial_values, operators, lr, epochs):\n",
    "\n",
    "    initial = initial_values\n",
    "    exp_values = []\n",
    "\n",
    "    for i in range(epochs):\n",
    "        exp, grad = ShiftRuleV3(qc, initial,operators)\n",
    "        exp_values.append(exp)\n",
    "        initial = initial - lr*grad\n",
    "        \n",
    "    return exp_values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Lowest Value: -1.9999981834066167\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAEWCAYAAABIVsEJAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAA1UElEQVR4nO3dd3gc1dn38e9PzXKvwsi9gjHYGBCmhmoTSoKB0EMntMAT0gikPSGEJJA8tBSKIYADvEBCb6H3YsAGYxuDuw3uVbblIkvW/f4xR84iVFZlNbvS/bmYa2fPlL13ZObeOWfmHJkZzjnnXH1lxR2Ac865zOQJxDnnXIN4AnHOOdcgnkCcc841iCcQ55xzDeIJxDnnXIN4AnF1ktRPUomk7AZuXyJpUDrFlAqp+J7NSdLVku4P82l3fF368QTSAkk6R9I0SZskLZN0m6Qu9dh+gaQxle/N7Asz62Bm2xoST9h2XkO2TVVMqdAU37MmkoZKekjSSknrJc2W9FdJfVLxeU15fCW9Lul7tSwfIMlCwiqRtFzSM5LGNvazU0XSvZKujTuOuHkCaWEk/QS4HrgC6AzsC/QHXpKUF2dsrmEkDQHeB5YAe5hZJ+AAYC5wYA3b5DRfhE2mi5l1AHYHXgIel3ROvCG5WpmZTy1kAjoBJcDJVco7ACuB88L7q4FHgIeBDcBHwO5h2X1ABbA57OtnwADAgJywzuvAtcC7YZ2nge7AA8B64ENgQMLnGzAE6BXWr5w2Rf8EDWAw8CqwGlgV9tWlHjH1Ap4C1gBzgAsSPv9q4F/AP8P3/RQoquEYfmW/Cd/3e2F+CPAGsC7E+XDV7xnmu4fjUnk8rgXerrLuxcBsoBj4O6AaYrofeLqOv/0hwCLgSmBZOGZdgWfC335tmO+TsM3A8F02EJ2w/wbcX91xIPox8g9gKbA4fJ/ssOwc4G3g/8LnzAeOCst+D2wDtoS/3d+SOeah/KfAciAr4W/8aPg+84EfJKw7GpgUjvdy4MaEZQcS/VstBr4EzgnlbULMX4RtbgfaVjmePwFWhO99blh2IVAGbA3fqda/TUueYg/Apyb8Y8KRQHnV/xHDsgnAg2H+6vA/wIlAbvgfdT6QG5YvAMYkbFv1ZPI60Ul6cDixzABmAWOAHKIT9T0J228/sVaJ6YGEmIYAY8P/1AXAm8DNCevWFdObwK1APjAqnGQOS/i+W4CjgWzgj8DEGo7h105mfDWBPAj8kujqPR84sLrvCTwUpnbA8HDiqppAngG6AP1CvEfWENMywkmvlr/9IeFvf304hm2Jkth3QgwdgX8DTyRs8x5wY1j/IKJEUlMCeRy4A2gP7AB8AFwUlp1D9O/pgnB8LyG6WlLV45fsMQ/lg0L5LuF4Twb+F8gLy+YB30z4LmeG+Q7AvmG+f/hepxH9W+8OjArLbiL60dEtHJ+ngT9WOZ7XhO2OJvrB0zUsvxe4Nu7/5+OevAqrZekBrDKz8mqWLQ3LK002s0fMrIzoJJJPVN2VrHvMbK6ZrQP+A8w1s5fDZ/8b2KO2jSVdCQwDzgMwszlm9pKZlZrZyhDTwckEIqkvUZXOlWa2xcymAHcBZyWs9raZPWdRnf59RNUkDVFGdFLqFT7r7WriySY6cf/GzDaZ2QyiBF7VdWZWbGZfAK8RJb7q9CBKIpX7v0xScWgvuDNhvYrwmaVmttnMVpvZoyGGDURXAweHffQD9gZ+HdZ/k+gE+jWSehKdQH9oZhvNbAXRyffUhNUWmtmd4fhOAAqBnjV8n2QtCa/dQqwFZnaNmW21qK3pzoQYyoAhknqYWYmZTQzlpwMvm9mDZlYWjskUSSK6kviRma0Jx+cPVb5TGXBN2O45oquNnRv5nVoUTyAtyyqgRw3134VheaUvK2fMrILocr1XPT5recL85mred6hpQ0lHAZcDx5nZ5lDWMzQSL5a0nqjapkdN+6iiF1B5Eqi0EOid8H5ZwvwmIL+B7QQ/AwR8IOlTSedVs04B0ZXYlwllX1azXtWYajpmq4n+fgCY2d/MrAtwM9Gv40orzWxL5RtJ7STdIWlhOKZvAl1CgusFrDWzjQnbL6zh8/uHz1kaElcx0dXIDtV9FzPbFGZr/DeQpMq/35oQQ6/Kzw8x/IL/JqnzgZ2AzyV9KOlbobwvUVtRVQVEV2aTE/b3fCivtLrKj7Ha/katkieQluU9oBQ4IbFQUgfgKOCVhOK+CcuzgD789xdfyrpolrQz0S/Uk80s8aT6h/C5IyxqJD6D6ERdqbaYlgDdJHVMKOtHVFdfX5Un1HYJZTtuD8JsmZldYGa9gIuAW0Mjd6KVRNUfiXdI9aXhXqHK37QGVY/RT4h+Me8TjulBoVxEV6RdJbVPWL9fDfv9kujfVQ8z6xKmTma2a5LxN/Tf0/FE7Q8zQwzzEz6/i5l1NLOjAcxstpmdRpTUrgceCd/tS6Kq1qpWEf3Q2TVhf50tasRP5XdqUTyBtCChOum3wF8lHSkpV9IAogbkRURVN5X2knRC+BX+Q6ITROVl/3KiOuYmJakT8CTwy2qqfjoSVRGsk9Sb6C6yRDXGFBLRu8AfJeVLGkn0i/T++sYYqs8WA2dIyg5XGNtPQJJOSrh1di3RiaSiyj62AY8BV4ergGF8tTqtvq4GviHpxnBskNSDqG2gNh2JTpLFkroBv0mIcSFRo/NvJeVJOhD4dnU7MbOlwIvADZI6ScqSNFhSUlWM1PPfU7gavSzE+/NwhfwBsEHSlZLahr/NbpL2DtucIakgrFscdlVB1M42RtLJknIkdZc0Kqx3J3CTpB3CPnpL+mYqvlNL5QmkhTGzPxFd2v8f0R0p7xP9CjvczEoTVn0SOIXoJHgmcEJoD4GokflX4dL+p00Y3p5Ev4hvSrjnvyQs+21Yvg54lugEnKiumE4jaoxdQtTg+xsze7mBcV5AlMBWA7sSJadKewPvh7ifAi636p/9uIzoBoPKO6IeJErS9WZms4B9iK5oPpG0AXiH6Lv+upZNbyZqTF9F9OPg+SrLTw/7XUN0sv5nLfs6i6jxegbRv5lHSKhWq8MtwImS1kr6Sy3rFUvaCEwjanM5yczuhu1J+VtE7UTzw3e6i+gYQ3QDyafh73ILcGpoB/oi7Osn4XtO4b/tX1cS3QwyMVTxvUzybRz/AIaHf49PJLlNi1N5l4RrRSRdTXS30Blxx9JaSLoe2NHMzo47Fueail+BOJcCkoZJGqnIaKIqtcfjjsu5ppSJT6s6lwk6ElVb9SKqL7+BqNrQuRYj9iosSUcS1VlmA3eZ2XVVlrchqpvdi6hO+hQzWxCW/Zzol902oqdSX2jG0J1zrlWLtQor3I/+d6JbTIcDp0kaXmW184nuVx9C9PDS9WHb4UQP/exK1IB2q7znUOecazZxV2GNBuZU3sUi6SFgHNGdHpXGEd3GCNGdH38LT5GOAx4KdxbNlzQn7O+9mj6sR48eNmDAgKb+Ds4516JNnjx5lZkVVC2PO4H05qtP6C4iuq2w2nXMrFzSOqL+bHrz3+cWKrftXWVbJF1I1GUB/fr1Y9KkSU0WvHPOtQaSqu2loMXfhWVm482syMyKCgq+lkCdc841UNwJZDFf7eKhD1/vfmL7OuGp6c5EjenJbOuccy5F4k4gHwJDJQ1UNNjRqURP9yZ6Cqh8+OpE4FWLbh17CjhVUhtJA4GhRN0dOOecawaxtoGENo3LgBeIbuO928w+lXQNMMnMniLqMuC+0Ei+htDdcljvX0QN7uXApZZGw5s651xLF/tzIM2pqKjIvBHdOefqR9JkMyuqWh53FZZzzrkM5QnEOedcg3gCScIH89dw11vzqKhoPdV9zjlXF08gSXj6kyVc++xnnDfhQ1aVNGhIB+eca3E8gSThmnG78rtxu/Lu3NUcdctbvD17Vd0bOedcC+cJJAmSOHO/ATx12QF0bpvLmXe/z/XPf07Ztoq6N3bOuRbKE0g9DNuxE09fdiCn7t2X216fy0m3v8eXazbFHZZzzsXCE0g9tc3L5o8njOTvp+/J3JUlHH3LWzz1yZK4w3LOuWbnCaSBjhlZyHM/+AZDe3bgBw9+zM8e+YRNW8vjDss555qNJ5BG6NutHQ9ftB+XHjqYf09exEm3v+dJxDnXangCaaTc7Cyu+OYwxp9ZxIyl6/nFY9NoTd3DOOdaL08gTWTs8J78eMxOPDFlCfdNrHbsFeeca1E8gTShSw8dwuHDduB3z8xg8sK1cYfjnHMp5QmkCWVliRtPHkVh57Z8/4HJrNzgT60751ouTyBNrHO7XG4/Yy+KN5XxPw9+RLk/bOica6E8gaTA8F6d+MPxI5g4bw1/fmFm3OE451xKeAJJke/s1Ycz9u3HHW/O4z/TlsYdjnPONblYEoikbpJekjQ7vHatZp1Rkt6T9KmkqZJOSVh2r6T5kqaEaVSzfoEk/fpbwxnVtwtXPDKVOStK4g7HOeeaVFxXIFcBr5jZUOCV8L6qTcBZZrYrcCRws6QuCcuvMLNRYZqS6oAbok1ONrd+d0/ycrK4+P7JbCz1hwydcy1HXAlkHDAhzE8Ajqu6gpnNMrPZYX4JsAIoaK4Am0qvLm3562l7MG9lCVc+OtUfMnTOtRhxJZCeZlbZMLAM6FnbypJGA3nA3ITi34eqrZsktall2wslTZI0aeXKlY0OvCEOGNKDK745jGemLuWedxbEEoNzzjW1lCUQSS9Lml7NNC5xPYt+ktf4s1xSIXAfcK6ZVd4T+3NgGLA30A24sqbtzWy8mRWZWVFBQXwXMBcfPIjDhu3ADS/OZN2mstjicM65ppKyBGJmY8xst2qmJ4HlITFUJogV1e1DUifgWeCXZjYxYd9LLVIK3AOMTtX3aCqSuOKbO7Nx6zbum7gg7nCcc67R4qrCego4O8yfDTxZdQVJecDjwD/N7JEqyyqTj4jaT6anMtimskthJw7ZuYB73lnAlrJtcYfjnHONElcCuQ4YK2k2MCa8R1KRpLvCOicDBwHnVHO77gOSpgHTgB7Atc0afSNccvBgVm/cyr8nL4o7FOecaxS1pruCioqKbNKkSbHGYGaccNu7rCop5bWfHEJOtj/L6ZxLb5Imm1lR1XI/ezUzSVx88GC+XLOZ56Yvizsc55xrME8gMRi7S08GF7Tn9tfn+nMhzrmM5QkkBllZ4qKDBzNj6Xremr0q7nCcc65BPIHE5LhRvdmxUz63vT637pWdcy4NeQKJSV5OFucfOJD35q3mky+L4w7HOefqzRNIjE7bpx+d8nO4/Q2/CnHOZR5PIDHq0CaHM/frz/OfLmPuSu/u3TmXWTyBxOyc/QeSl53FnW/OizsU55yrF08gMSvo2IaTivrw2EeLWbF+S9zhOOdc0jyBpIELvzGY8ooK/vHO/LhDcc65pHkCSQP9urfj6BGFPDDxC9Zt9q7enXOZwRNImrj44MGUlJbzwPsL4w7FOeeS4gkkTezWuzPfGNqDu9/2rt6dc5nBE0gaueTgwawqKeXRj7yrd+dc+vMEkkb2G9ydkX06M+HdBd7JonMu7XkCSSOSOKmoL7OWl/DZ0g1xh+Occ7WKLYFI6ibpJUmzw2vXGtbbljAi4VMJ5QMlvS9pjqSHwxC4Ge+YEYXkZIknP1kcdyjOOVerOK9ArgJeMbOhwCvhfXU2m9moMB2bUH49cJOZDQHWAuenNtzm0a19HgftVMDTU5ZQUeHVWM659BVnAhkHTAjzE4Djkt1QkoDDgEcasn26GzeqF0vWbeHDBWviDsU552oUZwLpaWZLw/wyoGcN6+VLmiRpoqTjQll3oNjMysP7RUDv1IXavMYO70nb3Gye/GRJ3KE451yNUppAJL0saXo107jE9Sy65aim+pr+YTD304GbJQ2uZwwXhgQ0aeXKlQ37Is2sXV4OR+zak+emLWVreUXc4TjnXLVSmkDMbIyZ7VbN9CSwXFIhQHhdUcM+FofXecDrwB7AaqCLpJywWh+g2lZnMxtvZkVmVlRQUNCk3y+Vxo3qRfGmMt6clRlJzznX+sRZhfUUcHaYPxt4suoKkrpKahPmewAHADPCFctrwIm1bZ/JvjG0gK7tcr0ayzmXtuJMINcBYyXNBsaE90gqknRXWGcXYJKkT4gSxnVmNiMsuxL4saQ5RG0i/2jW6FMsNzuLY0YW8tKMZZSUlte9gXPONbOculdJDTNbDRxeTfkk4Hth/l1gRA3bzwNGpzLGuI0b1Zv7J37BSzOWcfwefeIOxznnvsKfRE9je/XrSu8ubXlyildjOefSjyeQNJaVJY4d1Yu3Zq9iVUlp3OE459xXeAJJc+NG9WJbhfHctKV1r+ycc83IE0iaG7ZjJ4bt2NGrsZxzaccTSAY4dlQvJi9cy5drNsUdinPObecJJAMcu3svAJ7yZ0Kcc2nEE0gG6NO1HXsP6MoTHy/2gaacc2nDE0iGOHZUb2av8IGmnHPpI6kEIqm/pDFhvq2kjqkNy1XlA00559JNnQlE0gVE427cEYr6AE+kMCZXDR9oyjmXbpK5ArmUqBPD9QBmNhvYIZVBuer5QFPOuXSSTAIpNbOtlW9CF+r+EzgGPtCUcy6dJJNA3pD0C6CtpLHAv4GnUxuWq44PNOWcSyfJJJCrgJXANOAi4DngV6kMytXMB5pyzqWLOrtzN7MK4M4wuZhVDjT1+MeLGTO8pmHknXMu9epMIJLmU02bh5kNSklErla52VmcuFcf7n5nAYvWbqJP13Zxh+Sca6WSqcIqAvYO0zeAvwD3pzIoV7tzDxiIgH+8PT/uUJxzrVidCcTMVidMi83sZuCYxnyopG6SXpI0O7x2rWadQyVNSZi2SDouLLtX0vyEZaMaE0+m6dWlLcfu3ouHP/ySdZvK4g7HOddKJfMg4Z4JU5Gki2n8ULhXAa+Y2VDglfD+K8zsNTMbZWajgMOATcCLCatcUbnczKY0Mp6Mc8FBg9i0dRv3v78w7lCcc61UMonghoT5cmABcHIjP3cccEiYnwC8DlxZy/onAv8xM+/PPNilsBMH7VTAPe8s4PwDB5Kfmx13SM65ViaZKqxDE6axZnaBmc1s5Of2NLPKIfaWAXXdTnQq8GCVst9LmirpJkltatpQ0oWSJkmatHJly7r19aKDBrGqpJQnPvb+sZxzzU81dQ8u6ce1bWhmN9a6Y+llYMdqFv0SmGBmXRLWXWtmX2sHCcsKgalALzMrSyhbBuQB44G5ZnZNbfEAFBUV2aRJk+paLWOYGd/669tsLtvGyz86mKwsxR2Sc64FkjTZzIqqltdWhdWoHnfNbEwtwSyXVGhmS0MyWFHLrk4GHq9MHmHflVcvpZLuAX7amFgzlSQuPGgQlz80hVc+X8FYfy7EOdeMakwgZvbbFH7uU8DZwHXh9cla1j0N+HliQULyEXAcMD1Fcaa9Y0YU8qfnZzL+zbmeQJxzzSqZu7DyJV0q6VZJd1dOjfzc64CxkmYDY8J7wl1edyV89gCgL/BGle0fkDSNqHuVHsC1jYwnY+VkZ3H+gQP5cMFaPvpibdzhOOdakWQeJLyPqC3jm0Qn8j5Ao4bFC8+UHG5mQ81sjJmtCeWTzOx7CestMLPeoTuVxO0PM7MRZrabmZ1hZiWNiSfTnbJ3Xzq3zWX8G/PiDsU514okk0CGmNmvgY1mNoHoIcJ9UhuWq4/2bXI4Y99+vDBjGfNXbYw7HOdcK5FMAqlsvC6WtBvQGR9QKu2cvf8AcrOyuOstvwpxzjWPZBLI+NDVyK+JGr9nANenNCpXbzt0zOeEPXvzyORFrCopjTsc51wrUGMCkTRD0q+A18xsrZm9YWaDzGwHM7ujpu1cfL73jUGUllfwz/e8exPnXOrVdgVyGtAeeFHSB5J+FJ7ZcGlqyA4dGLNLT+57bwGbt26LOxznXAtXYwIxs0/M7OdmNhj4AdAPeF/Sa5IuaLYIXb1cdPAg1m4q49+Tv4w7FOdcC5dMGwhmNtHMfgScBXQB/pbKoFzDFfXvyh79unDXW/PZVlF9NzXOOdcUknmQcG9JN0paCFwN3AH0SnVgrmEkcdFBg/hizSaen74s7nCccy1YjV2ZSPoDcAqwBngIOMDMFjVXYK7hxg7fkUEF7fn5Y1Pp0i6XA4b0iDsk51wLVNsVyBbgSDPb28xu8OSRObKzxIRzR7Nj53zOuvsD7p/od2U555pebY3o15jZ7OYMxjWdvt3a8egl+3PQ0B786onp/ObJ6ZRvq6h7Q+ecS1JSjeguM3XMz+Wus/fmewcOZMJ7Czn33g9Zt9nHUHfONQ1PIC1cdpb41beGc/13RjBx3mqOv/Ud7y/LOdckkkogknpL2l/SQZVTqgNzTeuUvftx//n7sHbjVo77+zu8O2dV3CE55zJcMrfxXg+8A/wKuCJMrXIEwEy3z6DuPHnpgezQsQ1n3f0BD7zvjevOuYarcUz07StIM4GRZpbxPfS1tDHRG2rDljL+58GPeX3mSvYb1J0Dh/Zg30HdGdmnM7nZXqvpnPuqhoyJXmkekAtkfAJxkY75ufzj7L259bU5PDN1KX9+YSYA7fKyKRrQjX0HdWPfQd0Z0dsTinOuZslcgTwK7A68QkISMbMfNOqDpZOInmzfBRhtZtVeGkg6ErgFyAbuMrPK4W8HEj3g2B2YDJxpZltr+0y/Aqne6pJS3p+/honzVjNx3mpmLY8GeGwfEsouhZ3o07Utvbu2pW/XtvTu0o62edkxR+2cay6NuQJ5KkxNbTpwAlHXKNWSlA38HRgLLAI+lPSUmVWOSXKTmT0k6XbgfOC2FMTZ4nXv0IajRxRy9Iios+VVJaW8P28N781bxfvz1vDu3FWUbfvqD43u7fO2J5U+XduxQ8c2dO+QR/f2bejWPo8eHaLXvBy/gnGupaozgZjZBEl5wE6haKaZNfphAjP7DKK+m2oxGphjZvPCug8B4yR9BhwGnB7Wm0B0NeMJpAn06NCGY0YWcszIKKFsqzBWbihl0dpNLC7ezKK1m1m0dhOL1m7m86UbePmzFWwtr/4hxU75OXTv0Ibu7fPo0i6Pzm1z6dIu9yuv0Xze9vmO+TledeZcBqgzgUg6hOgEvQAQ0FfS2Wb2Zkoji/QGEvslX0Q0Hnt3oNjMyhPKe1e3A0kXAhcC9OvXL3WRtmDZWWLHzvns2Dmfr13DAmbG+i3lrC4pZc3Grawq2crqjaWsLtka3kfzi4s3M2PJOtZtLmNjHeOVtM3NplPbHDrl59KpbS6d8nPCa26V8uh95+3znoCcay7JVGHdABxhZjMBJO0EPAjsVdeGkl4Gdqxm0S/N7Mn6BNpQZjYeGA9RG0hzfGZrI2n71cOgguS22VpewfotZRRvKmPd5jLWbd5K8aYyNmwpZ/3mMtZvKWP95vLodUsZqzduZf6qjazbXMb6LeV1dlXfoU2UVLq2z6VL2zw6t8ula7tovvLKp0eHNvToEKreOuTRJsfbdZyrj2QSSG5l8gAws1mScpPZuZmNaXBkkcVA34T3fULZaqCLpJxwFVJZ7jJEXk7W9hN4fZkZm7Zu+2qS2Rwloui1nHWbyygOSal401aWrNvMuk1lFG8uqzH5dMzPiRJKaMPp3iGPws75FHZuG712iV7zcz3ROAfJJZBJku4C7g/vvws0161MHwJDwx1Xi4FTgdPNzCS9BpxIdCfW2UCzXNG4+EmifZsc2rfJobBz/bY1MzaUllO8sYzVG0uj6raSUlaVhPmNW1m1oZR5q0p4f34pazd9vbmvW/s8duyUT68u+fTq0pYB3dszqKA9g3p0oHfXtmRn1dqu51yLkcxtvG2AS4EDQ9FbwK2NfbBQ0vHAX4ECoBiYYmbflNSL6Hbdo8N6RwM3E93Ge7eZ/T6UDyJKHt2Aj4Ez6orJb+N19bWlbBtL121hafHm6HXdZpYkvF+8djMbSsu3r5+XnUX/7u0YVNCegT06MKigPYMLOjC8sJPf+uwyVk238daZQFoSTyCuqZkZqzduZd7KjcxfVcK8lRuZG+a/WLNp++3P2Vli6A4d2L1PF0b27czufbqw844dvbHfZYR6Pwci6V9mdrKkacDXsoyZjWziGJ3LOJK2t+WMHtjtK8vKt1WwaO1mZi3fwLTF6/hk0TpemLGMhydFNxbm5WQxvLATu/fpzB79uvKNoT3o3oA2IefiUuMViKRCM1sqqX91y80s43ri8ysQFzcz44s1m/hk0TqmflnM1EXrmL5kHZu2bkOCkX26cNjOO3DYsB3YtVcnsrw9xaWBBldhSbrezK6sqywTeAJx6WhbhfHpknW89vlKXp25gqmLijGDgo5tOHTnAg4btgMHDi2gQ5tk7nlxruk1JoF8ZGZ7VimbmolVWJ5AXCZYVVLKGzOjZPLmrJVs2FJObrbYd1B3vrtPP8bs0pMcbztxzajeCUTSJcD3gUHA3IRFHYF3zOyMVASaSp5AXKYp21bB5IVree3zFTwzdSmLizdT2Dmf00f349TR/Sjo6G0mLvUakkA6A12BPwJXJSzaYGZrUhJlinkCcZlsW4XxymfLuW/iQt6avYrcbHH0iELO2m8Ae/brUle/cs41WKNv45W0A5Bf+d7Mvmi68JqHJxDXUsxdWcJ97y3k0cmL2FBazq69OnHWfv05dvfe/ryJa3KNaQP5NnAj0AtYAfQHPjOzXVMRaCp5AnEtzcbSch7/eDH3vbeQmcs30L19HlcdNYwT9+rjVySuydSUQJJpibsW2BeYZWYDgcOBiU0cn3OuAdq3yeGMffvz/A+/wYMX7Ev/7u244pGpnHzHe3y+bH3c4bkWLpkEUmZmq4EsSVlm9hpU26u3cy4mkthvcHceuXh/rv/OCOasKOGYv7zN75+dQUlCVyvONaVkEkixpA7Am8ADkm4BNqY2LOdcQ2RliVP27serPzmEk4v6cOdb8xlzwxs8N20pranbItc8kkkg44BNwI+A54lu6f1WKoNyzjVO1/Z5/PGEkTx6yf50a5/H9x/4iLPv+ZAFq/y3n2s6ySSQ/zWzCjMrN7MJZvYXIOOeQneuNdqrf1eeuuwAfvPt4Xy0cC1H3Pwmf31lNhV1DMjlXDKSSSBjqyk7qqkDcc6lRk52FuceMJBXf3IwRwzvyQ0vzeLyh6dQWl77sMLO1aW23ngrn0QfLGlqwqKOwLupDsw517R26JTPX0/bg916d+a6/3zOivVbGH9WEZ3bJjXAqHNfU9sVyP8Dvk000t+3E6a9zOy7zRCbc66JSeLigwdzy6mj+OiLtZx0+7ssLt4cd1guQ9WYQMxsnZktAG4B1pjZwtCFe7mkfZorQOdc0xs3qjcTzhvN0uItnHDrO8xY4s+MuPpLpg3kNqAk4X1JKGswSSdJ+lRShaRqnymR1FfSa5JmhHUvT1h2taTFkqaE6ejGxONca7T/4B78+5L9yJI4+Y73eGfOqrhDchkmmQQiS7iB3MwqqKXtJEnTgROIni2pSTnwEzMbTvQk/KWShicsv8nMRoXpuUbG41yrNGzHTjz2/f3p07UtZ9/9AY9/vCjukFwGSSaBzJP0A0m5YbocmNeYDzWzz8xsZh3rLDWzj8L8BuAzoHdjPtc593WFndvyr4v3Y/TAbvzo4U/4+2tz/KFDl5RkEsjFwP7AYmARsA9wYSqDqkrSAGAP4P2E4sskTZV0t6SutWx7oaRJkiatXLky1aE6l5E65edy77mjOW5UL/78wkx++/QMTyKuTnVWRZnZCuDU+u5Y0svAjtUs+qWZPVmP/XQAHgV+aGaVLX23Ab8DLLzeAJxX3fZmNh4YD1FvvEl/AedambycLG48eRQ9OrThrrfnM6B7O845YGDcYbk0VmcCkbQT0Qm7p5ntJmkkcKyZXVvbdmY2prHBScolSh4PmNljCftenrDOncAzjf0s51zUl9Yvjt6FhWs28btnP2NYYSf2HdQ97rBcmkqmCutO4OdAGYCZTaUBVyT1pWgwg38QjT1yY5VlhQlvjydqlHfONYGsLHHjybszoHs7Ln3gI5b4cyKuBskkkHZm9kGVskb1Dy3peEmLgP2AZyW9EMp7Saq8o+oA4EzgsGpu1/2TpGnhCflDiTp6dM41kY75uYw/q4jS8gouvn8yW8q82xP3dcncjrtK0mCi9gYknQgsbcyHmtnjwOPVlC8Bjg7zbwPVDqlmZmc25vOdc3UbXNCBm04ZxQX/nMSvnpjOn08c6aMcuq9I5grkUuAOYJikxcAPie7Mcs61cGOH9+Tyw4fyyORF3DdxYdzhuDSTzF1Y84AxktoDWeGZDOdcK3H54UP5dMk6rnl6Bjv37Mg+3qjugjqvQCR1l/QX4C3gdUm3SPJ/Qc61EllZ4sZTRtGvWzsu/X8fsXSdN6q7SDJVWA8BK4HvACeG+YdTGZRzLr10ys9l/Fl7saWsgovv/8gb1R2QXAIpNLPfmdn8MF0L9Ex1YM659DJkh47ccPLufPJlMf/75HR/Ut0llUBelHSqpKwwnQy8kOrAnHPp55u77sgPDhvCvyYt4v73v4g7HBezZBLIBUSDS5WG6SHgIkkbJPkgAs61Mj8csxOH7lzAtc/M4Ms1m+IOx8WozgRiZh3NLMvMcsOUFco6mlmn5gjSOZc+srLEH04YQXaW+O3TM+IOx8Uombuwzq/yPlvSb1IXknMu3RV2bssPDh/Ky58t59XPl9e9gWuRkqnCOlzSc5IKJe0GTAQ6pjgu51yaO++AgQwuaM/VT83wu7JaqWSqsE4HJgDTgGeJulX/aaoDc86lt7ycLK4ZtxtfrNnEHW80aow5l6GSqcIaClxO1K36QuBMSe1SHZhzLv0dMKQHx4ws5NbX53iDeiuUTBXW08Cvzewi4GBgNvBhSqNyzmWMXx2zizeot1LJJJDRZvYKgEVuIBqDwznnvEG9FasxgUj6GYCZrZd0UpXF56QyKOdcZvEG9daptiuQxFEHf15l2ZEpiMU5l6G8Qb11qi2BqIb56t4751o5b1BvfWpLIFbDfHXv60XSSZI+lVQhqaiW9RaEoWunSJqUUN5N0kuSZofXro2JxznXNLxBvXWpLYHsLmm9pA3AyDBf+X5EIz93OnAC8GYS6x5qZqPMLDHRXAW8YmZDgVfCe+dczAo7t+Vyb1BvNWpMIGaWbWadQp9XOWG+8n1uYz7UzD4zs5mN2MU4oocbCa/HNSYe51zTOdcb1FuNZG7jjZMRdSc/WdKFCeU9zWxpmF9GLeOTSLpQ0iRJk1auXJnKWJ1zeIN6a5KyBCLpZUnTq5nG1WM3B5rZnsBRwKWSDqq6gkWj2tTYJmNm482syMyKCgoK6v9FnHP1dsCQHhwzopDb35jLqpLSuMNxKZKyBGJmY8xst2qmJ+uxj8XhdQXwODA6LFouqRAgvK5o6vidc43zkyN2Yuu2Cm57fW7cobgUSdsqLEntJXWsnAeOIGp8B3gKODvMnw0knZScc81jUEEHvrNnb+6buJCl6zbHHY5LgVgSiKTjJS0C9gOelfRCKO8l6bmwWk/gbUmfAB8Az5rZ82HZdcBYSbOBMeG9cy7N/M9hQzEz/vbqnLhDcSmQE8eHmtnjRFVSVcuXAEeH+XnA7jVsvxo4PJUxOucar2+3dpy6dz8e/OALLjpoMP26e0feLUnaVmE551qGyw4bQnaWuOWV2XGH4pqYJxDnXEr17JTPWfv15/GPFzFnxYa4w3FNyBOIcy7lLj54MG1zs7npZb8KaUk8gTjnUq57hzacd+BAnp26lE+XrIs7HNdEPIE455rF974xiE75Odz00qy4Q3FNxBOIc65ZdG6by0UHD+blz1bw8Rdr4w7HNQFPIM65ZnPO/gPo3j6PG170q5CWwBOIc67ZtG+TwyWHDObtOat4b+7quMNxjeQJxDnXrM7Ytz89O7XhhhdnEvWF6jKVJxDnXLPKz83mssOGMmnhWt6Y5UMsZDJPIM65ZndKUV/6dG3LDS/O8quQDOYJxDnX7PJysrj88KFMW7yOFz71oW8zlScQ51wsjt+jN4MK2vPnFz6nbFtF3OG4BvAE4pyLRU52FlceOYy5Kzfy0AdfxB2OawBPIM652BwxvCf7DurGjS/NYt3msrjDcfXkCcQ5FxtJ/OqY4RRvLuNvr3pHi5nGE4hzLla79e7MiXv24d53F7Bw9ca4w3H1ENeQtidJ+lRShaSiGtbZWdKUhGm9pB+GZVdLWpyw7Ohm/QLOuSZ1xTd3Jjc7iz8+93ncobh6iOsKZDpwAvBmTSuY2UwzG2Vmo4C9gE18dRjcmyqXm9lz1e7EOZcRduiUzyUHD+b5T5cxcZ53cZIpYkkgZvaZmc2sxyaHA3PNbGGqYnLOxeuCgwbRq3M+1z47g4oKf7gwE2RKG8ipwINVyi6TNFXS3ZK61rShpAslTZI0aeVK7zbBuXSVn5vNz44cxvTF63ns48Vxh+OSkLIEIullSdOrmcbVcz95wLHAvxOKbwMGA6OApcANNW1vZuPNrMjMigoKCur/RZxzzebY3Xuxe98u/PmFz9m0tTzucFwdUpZAzGyMme1WzfRkPXd1FPCRmW3v78DMlpvZNjOrAO4ERjdl7M65eGRlif/91i4sX1/KHW/MizscV4dMqMI6jSrVV5IKE94eT9Qo75xrAfbq341jRhZyx5tzWbpuc9zhuFrEdRvv8ZIWAfsBz0p6IZT3kvRcwnrtgbHAY1V28SdJ0yRNBQ4FftRMoTvnmsFVRw6jwuDPz9fnXhvX3HLi+FAze5yv3pJbWb4EODrh/UagezXrnZnSAJ1zserbrR3nHziQ216fyzkHDGBkny5xh+SqkQlVWM65Vuj7hwymR4c8fvfMDB8zJE15AnHOpaWO+bn8eOzOfLhgLf+ZvizucFw1PIE459LWKXv3ZZfCTvzy8WksWOX9ZKUbTyDOubSVnSVu++6eAJx374es3bg15ohcIk8gzrm0NqBHe8afVcSitZu56P7JlJZvizskF3gCcc6lvb0HdOPPJ43kg/lruOrRad6oniZiuY3XOefqa9yo3ny5ZhP/9+Is+nVrx4/G7hR3SK2eJxDnXMa49NAhLFi9iVtemU2/bu34zl594g6pVfME4pzLGJL4w/EjWFK8masem0qvLm3Zb/DXnjV2zcTbQJxzGSUvJ4vbztiL/t3bc9F9k5izoiTukFotTyDOuYzTuW0u95yzN7nZWZx374esLimNO6RWyROIcy4j9e3WjrvOLmL5+i1c8M9JbCnz23ubmycQ51zG2qNfV24+ZRQffVHMaXdOZPLCNXGH1Kp4AnHOZbSjRhRyy6mjWLR2M9+57T2+N2ESs5ZviDusVsETiHMu440b1Zs3rjiEnx6xE+/PW82RN7/JT//9CYuLfUCqVFJreqKzqKjIJk2aFHcYzrkUWrtxK39/bQ7/fG8hCM7atz+XHjqEru3z4g4tY0mabGZFXyv3BOKca4kWF2/mppdm8dhHi2ifl8NFBw/i9H36080TSb2lXQKR9Gfg28BWYC5wrpkVV7PekcAtQDZwl5ldF8oHAg8RjVg4GTjTzGrtqtMTiHOtz6zlG/jT8zN5+bPlAOzYKZ/hvToxvLDT9td+3dqRlaWYI01f6ZhAjgBeNbNySdcDmNmVVdbJBmYRjYu+CPgQOM3MZkj6F/CYmT0k6XbgEzO7rbbP9ATiXOs1dVEx789bw4yl65mxZD1zVpawrSI6/7XPy2ZYYZRMduycT8f8HDq0CVN+Dp3yc7fPd2iTQ5ucLKTWk3BqSiCxdWViZi8mvJ0InFjNaqOBOWY2D0DSQ8A4SZ8BhwGnh/UmAFcDtSYQ51zrNbJPl6+Mrb6lbBuzl5cwY+k6Plu6gRlL1vPEx4vZUFqe1P6yBDlZWWRnafuUkzCfFRKMFCYUXqMuWQTRmyqqFjVVovrD8SMYPbBbk+yrUrr0hXUe8HA15b2BLxPeLwL2Iaq2Kjaz8oTy3tXtWNKFwIUA/fr1a6p4nXMZLj83mxF9OjOiT+ftZWZGaXkFJaXllGwpp6S0nA3htaS0jJIt5WwoLaes3NhWUUF5hbGtwhJeK6LXbUaFgWGE/zCz8Prf91V9raQJK4jat8luup0FKU0gkl4Gdqxm0S/N7Mmwzi+BcuCBVMRgZuOB8RBVYaXiM5xzLYMk8nOzyc/NpkeHNnGHk/ZSmkDMbExtyyWdA3wLONyqb4xZDPRNeN8nlK0GukjKCVchleXOOeeaSWwPEoa7q34GHGtmm2pY7UNgqKSBkvKAU4GnQrJ5jf+2m5wNPJnqmJ1zzv1XnE+i/w3oCLwkaUq4kwpJvSQ9BxCuLi4DXgA+A/5lZp+G7a8EfixpDlGbyD+a+ws451xrFuddWENqKF8CHJ3w/jnguWrWm0d0l5ZzzrkYeF9YzjnnGsQTiHPOuQbxBOKcc65BPIE455xrkFbVG6+klcDCBm7eA1jVhOGkWibFm0mxQmbFm0mxQmbFm0mxQuPi7W9mBVULW1UCaQxJk6rrTCxdZVK8mRQrZFa8mRQrZFa8mRQrpCZer8JyzjnXIJ5AnHPONYgnkOSNjzuAesqkeDMpVsiseDMpVsiseDMpVkhBvN4G4pxzrkH8CsQ551yDeAJxzjnXIJ5AkiDpSEkzJc2RdFXc8dRG0gJJ00IPx2k3ALykuyWtkDQ9oaybpJckzQ6vXeOMsVINsV4taXE4vlMkHV3bPpqTpL6SXpM0Q9Knki4P5Wl3fGuJNS2Pr6R8SR9I+iTE+9tQPlDS++Hc8HAYdiJdY71X0vyEYzuq0Z/lbSC1k5QNzALGEg2d+yFwmpnNiDWwGkhaABSZWVo+4CTpIKAE+KeZ7RbK/gSsMbPrQoLuamZXxhlniKu6WK8GSszs/+KMrTqSCoFCM/tIUkdgMnAccA5pdnxrifVk0vD4KhqYvL2ZlUjKBd4GLgd+DDxmZg+FISk+MbPb0jTWi4FnzOyRpvosvwKp22hgjpnNM7OtwEPAuJhjylhm9iawpkrxOGBCmJ9AdCKJXQ2xpi0zW2pmH4X5DURj6PQmDY9vLbGmJYuUhLe5YTLgMKDyhJwux7amWJucJ5C69Qa+THi/iDT+h070D+VFSZMlXRh3MEnqaWZLw/wyoGecwSThMklTQxVX7NVB1ZE0ANgDeJ80P75VYoU0Pb6SsiVNAVYALwFzgeIw8B2k0bmhaqxmVnlsfx+O7U2SGj3ouyeQludAM9sTOAq4NFTDZIwwXHE616veBgwGRgFLgRtijaYakjoAjwI/NLP1icvS7fhWE2vaHl8z22Zmo4A+RDUTw+KNqGZVY5W0G/Bzopj3BroRjeraKJ5A6rYY6Jvwvk8oS0tmtji8rgAeJzNGbVwe6sQr68ZXxBxPjcxsefifswK4kzQ7vqHO+1HgATN7LBSn5fGtLtZ0P74AZlYMvAbsB3SRVDmya9qdGxJiPTJUG5qZlQL30ATH1hNI3T4Ehoa7LfKAU4GnYo6pWpLahwZJJLUHjgCm175VWngKODvMnw08GWMstao8EQfHk0bHNzSe/gP4zMxuTFiUdse3pljT9fhKKpDUJcy3Jbqp5jOik/OJYbV0ObbVxfp5wo8IEbXVNPrY+l1YSQi3Et4MZAN3m9nv442oepIGEV11QDTe/f9Lt1glPQgcQtS19HLgN8ATwL+AfkTd7Z9sZrE3XtcQ6yFE1SsGLAAuSmhfiJWkA4G3gGlARSj+BVHbQlod31piPY00PL6SRhI1kmcT/fD+l5ldE/6fe4ioSuhj4IzwCz82tcT6KlAACJgCXJzQ2N6wz/IE4pxzriG8Css551yDeAJxzjnXIJ5AnHPONYgnEOeccw3iCcQ551yDeAJxrglI2pbQy+kUNWGvzZIGKKFHYOfSRU7dqzjnkrA5dB3hXKvhVyDOpZCi8Vn+pGiMlg8kDQnlAyS9Gjq2e0VSv1DeU9LjYSyHTyTtH3aVLenOML7Di+EJYyT9QNGYGlMlPRTT13StlCcQ55pG2ypVWKckLFtnZiOAvxH1aADwV2CCmY0EHgD+Esr/ArxhZrsDewKfhvKhwN/NbFegGPhOKL8K2CPs5+LUfDXnqudPojvXBCSVmFmHasoXAIeZ2bzQeeAyM+suaRXRgEploXypmfWQtBLok9gdRuju/CUzGxreXwnkmtm1kp4nGvTqCeCJxnZN4Vx9+BWIc6lnNczXR2L/Stv4b/vlMcDfia5WPkzoGda5lPME4lzqnZLw+l6Yf5eoZ2eA7xJ1LAjwCnAJbB8UqHNNO5WUBfQ1s9eIxnboDHztKsi5VPFfK841jbZhBLhKz5tZ5a28XSVNJbqKOC2U/Q9wj6QrgJXAuaH8cmC8pPOJrjQuIRpYqTrZwP0hyQj4Sxj/wblm4W0gzqVQaAMpMrNVccfiXFPzKiznnHMN4lcgzjnnGsSvQJxzzjWIJxDnnHMN4gnEOedcg3gCcc451yCeQJxzzjXI/wcKV8lmodo/2AAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "epochs = 35\n",
    "\n",
    "epochs_line = np.linspace(0,epochs, epochs)\n",
    "lr = 0.5\n",
    "\n",
    "exp_value = Optimizer(qc, initial,H, lr, epochs)\n",
    "\n",
    "print(\"Lowest Value: {}\".format(exp_value[-1]))\n",
    "plt.plot(epochs_line, exp_value)\n",
    "plt.ylabel(\"Expectation Value\")\n",
    "plt.xlabel(\"Epochs\")\n",
    "plt.title(\"Optimization usign Gradient Descent\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Comparison with Pennylane"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Initial Parameters: [0.5 0.5 0.5 0.5 0.5 0.5]\n",
      "Initial Expectation Value: 0.6758712218347047\n"
     ]
    }
   ],
   "source": [
    "# here we are just using the same value for all the parameters.\n",
    "initial = np.array([0.5,0.5,0.5,0.5,0.5,0.5])\n",
    "\n",
    "# MEasure the Z @ Z operator for simplicity.\n",
    "H = [R4]\n",
    "\n",
    "exp_value, gradient = ShiftRuleV3(qc, initial,H)\n",
    "print(\"Initial Parameters: {}\".format(initial))\n",
    "print(\"Initial Expectation Value: {}\".format(exp_value))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Pennylane implementations:\n",
    "from pennylane import numpy as np\n",
    "import pennylane as qml\n",
    "\n",
    "dev = qml.device(\"default.qubit\", wires=2)\n",
    "\n",
    "@qml.qnode(dev)\n",
    "def circuit(params):\n",
    "    qml.RX(params[0], wires=0)\n",
    "    qml.RX(params[0], wires=1)\n",
    "\n",
    "    qml.CNOT(wires=[0,1])\n",
    "\n",
    "    qml.RY(params[0], wires=0)\n",
    "    qml.RY(params[0], wires=1)\n",
    "\n",
    "    qml.CNOT(wires=[0,1])\n",
    "\n",
    "    qml.RZ(params[0], wires=0)\n",
    "    qml.RZ(params[0], wires=1)\n",
    "\n",
    "    return qml.expval(qml.PauliZ(0) @ qml.PauliZ(1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "tensor(0.67587122, requires_grad=True)"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "par = np.array([0.5], requires_grad=True)\n",
    "\n",
    "circuit(par)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Excelent! There are the same values. With this code you can generalize the parameter shift rule for any number of qubits."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "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.9.5"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "b692c5112b84cb0f0146006066a71ec03118a3195a234c45307942b05862c083"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}