VQE_Tutorial.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Variational Quantum Eigensolver Tutorial"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So you want to find the lowest eigenvalue/eigenvector of a very large matrix. Maybe that matrix represents the energy operator of some quantum system whose ground state you're interested in. Maybe you have other motives. Either way, since the matrix is very large, it's going to be inconvenient to directly solve for its eigenvalues and eigenvectors. On the other hand, if we actually had a quantum system, we could evaluate expectation values like $\\langle \\psi \\mid H \\mid \\psi \\rangle$. Could we use that to help us? \n",
    "\n",
    "Suppose we had some ansatz for the eigenvector that takes some parameters, for example, an angle $\\theta$: $\\mid \\alpha(\\theta) \\rangle$. We could search over $\\langle \\alpha(\\theta) \\mid H \\mid \\alpha(\\theta) \\rangle$ until we find a value of $\\theta$ that minimizes the expectation value of the operator, and then we've solved our problem.\n",
    "\n",
    "In order to do this, however, we have to figure out how to actually measure $\\langle \\psi \\mid H \\mid \\psi \\rangle$ on our quantum computer. All we have are $n$ qubits, the ability to implement unitary gates, and the ability to measure each of the qubits along the $Z$ axis, obtaining a bit: $0$ or $1$. So what do we do?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr>\n",
    "\n",
    "For example, suppose we have $H = \\begin{pmatrix} 4 & -2i \\\\ 2i & 2 \\end{pmatrix}$.\n",
    "\n",
    "Let's stare at the Pauli matrices:\n",
    "\n",
    "$I = \\begin{pmatrix} 1 & 0 \\\\ 0 & 1 \\end{pmatrix}, X = \\begin{pmatrix} 0 & 1 \\\\ 1 & 0 \\end{pmatrix}, Y = \\begin{pmatrix} 0 & -i \\\\ i & 0 \\end{pmatrix}, Z = \\begin{pmatrix} 1 & 0 \\\\ 0 & -1 \\end{pmatrix} $\n",
    "\n",
    "After a while, we might realize that $H$ can be written $3 + 2Y + Z$.\n",
    "\n",
    "$\\begin{pmatrix} 4 & -2i \\\\ 2i & 2 \\end{pmatrix} = 3\\begin{pmatrix} 1 & 0 \\\\ 0 & 1 \\end{pmatrix} + 2\\begin{pmatrix} 0 & -i \\\\ i & 0 \\end{pmatrix} + \\begin{pmatrix} 1 & 0 \\\\ 0 & -1 \\end{pmatrix} $\n",
    "\n",
    "By the laws of expectation values, it follows that:\n",
    "\n",
    "$\\langle H \\rangle = 3+ 2\\langle Y \\rangle + \\langle Z \\rangle$\n",
    "\n",
    "In other words, instead of measuring $\\langle H \\rangle$ directly, we could measure $\\langle Y \\rangle$ and $\\langle Z \\rangle$ separately (on separate preparations of the system), and add up the results, appropriately weighted, to obtain $\\langle H \\rangle$. In other words, we can split our expectation value into a weighted sum of expectation values of Pauli's.\n",
    "\n",
    "Indeed, $I, X, Y, Z$ form a basis for $2 \\times 2$ Hermitian operators, and so we can decompose any $2\\times 2$ Hermitian operator into a sum as we did above: \n",
    "\n",
    "$H = \\frac{tr(H)}{2}I + \\frac{tr(XH)}{2}X + \\frac{tr(YH)}{2}Y + \\frac{tr(ZH)}{2}Z$\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "our matrix:\n",
      "Quantum object: dims = [[2], [2]], shape = (2, 2), type = oper, isherm = True\n",
      "Qobj data =\n",
      "[[4.+0.j 0.-2.j]\n",
      " [0.+2.j 2.+0.j]]\n",
      "\n",
      "pauli decomposition:\n",
      "I: 3.00\n",
      "Y: 2.00\n",
      "Z: 1.00\n",
      "\n",
      "still our matrix:\n",
      "Quantum object: dims = [[2], [2]], shape = (2, 2), type = oper, isherm = True\n",
      "Qobj data =\n",
      "[[4.+0.j 0.-2.j]\n",
      " [0.+2.j 2.+0.j]]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "import qutip as qt\n",
    "\n",
    "########################################################################\n",
    "\n",
    "def construct_pauli_basis():\n",
    "    return {\"I\": qt.identity(2),\\\n",
    "            \"X\": qt.sigmax(),\\\n",
    "            \"Y\": qt.sigmay(),\\\n",
    "            \"Z\": qt.sigmaz()}\n",
    "\n",
    "def op_pauli(O, basis):\n",
    "    return dict([(s, (o.dag()*O).tr()/np.sqrt(len(basis)))\\\n",
    "                    for s, o in basis.items()])\n",
    "\n",
    "def pauli_op(P, basis):\n",
    "    return sum([P[s]*o for s, o in basis.items()])\n",
    "\n",
    "def print_pauli(P):\n",
    "    for s, v in P.items():\n",
    "        if not np.isclose(v, 0):\n",
    "            print(\"%s: %.2f\" % (s, v.real))\n",
    "\n",
    "########################################################################\n",
    "\n",
    "basis = construct_pauli_basis()\n",
    "H = qt.Qobj(np.array([[4, -2j],\\\n",
    "                      [2j, 2]]))\n",
    "H_pauli = op_pauli(H, basis)\n",
    "H2 = pauli_op(H_pauli, basis)\n",
    "\n",
    "print(\"our matrix:\")\n",
    "print(H)\n",
    "print(\"\\npauli decomposition:\")\n",
    "print_pauli(H_pauli)\n",
    "print(\"\\nstill our matrix:\")\n",
    "print(H2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can check that $\\langle H \\rangle = 3+ 2\\langle Y \\rangle + \\langle Z \\rangle$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<H> calculated the usual way: 3.6952\n",
      "<H> calculated the pauli way: 3.6952\n"
     ]
    }
   ],
   "source": [
    "# We take a sum over each operator and non-zero weight in the \n",
    "# Pauli decomposition. If we're dealing with the identity matrix, \n",
    "# we just add the weight to the sum; otherwise, we take the expectation value \n",
    "# of that Pauli operator, multiplied by its weight, and add that to the pile.\n",
    "def expect(H, psi, basis):\n",
    "    H_pauli = op_pauli(H, basis)\n",
    "    return sum([v*(psi.dag()*basis[s]*psi\\\n",
    "                       if s.count(\"I\") != len(s)\\\n",
    "                   else 1)\\\n",
    "                if not np.isclose(v, 0) else 0\n",
    "                    for s, v in H_pauli.items()])[0,0].real\n",
    "    \n",
    "psi = qt.rand_ket(2)\n",
    "print(\"<H> calculated the usual way: %.4f\" % (psi.dag()*H*psi)[0,0].real)\n",
    "print(\"<H> calculated the pauli way: %.4f\" % expect(H, psi, basis))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr>\n",
    "\n",
    "So we've reduced the problem of measuring $\\langle H \\rangle$ to that of measuring the $\\langle X \\rangle$/$\\langle Y \\rangle$/$\\langle Z \\rangle$'s separately (the constant term doesn't require any actual measurement). $\\langle Z \\rangle$ presents no problem for our quantum computer--it's the default: we simply prepare and measure some large number $N$ of identical copies of the initial qubit: we get $\\uparrow$ or $\\downarrow$ each time. Over the course of the experiments, we accumulate counts, the number $N_{\\uparrow}$ and the number $N_{\\downarrow}$, which we divide by the total number $N$ of experiments, to obtain the probabilities $p_{\\uparrow} = \\frac{N_{\\uparrow}}{N}, p_{\\downarrow} = \\frac{N_{\\downarrow}}{N}$. We associate $\\uparrow$ with eigenvalue $\\lambda_{\\uparrow} = 1$ and $\\downarrow$ with eigenvalue $\\lambda_{\\downarrow} = -1$, and sum them, weighted by their probabilities. This is the expectation value $\\langle \\psi \\mid Z \\mid \\psi \\rangle = (1)p_{\\uparrow} + (-1)p_{\\uparrow}$. You could just imagine: we write down a $1$ or a $-1$ depending on the outcome, and them sum them up and divide by $N$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<Z> calculated the usual way: -0.1198\n",
      "<Z> estimated stochastically: -0.1152\n"
     ]
    }
   ],
   "source": [
    "def expectZ(psi, n_shots=5000):\n",
    "    dm = psi*psi.dag()\n",
    "    return sum([np.random.choice([1, -1],\\\n",
    "                        p=[dm[0,0].real, dm[1,1].real])\\\n",
    "                    for i in range(n_shots)])/n_shots\n",
    "\n",
    "psi = qt.rand_ket(2)\n",
    "print(\"<Z> calculated the usual way: %.4f\" % qt.expect(qt.sigmaz(), psi))\n",
    "print(\"<Z> estimated stochastically: %.4f\" % (expectZ(psi)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr>\n",
    "\n",
    "But what about $\\langle X \\rangle$ or $\\langle Y \\rangle$?\n",
    "\n",
    "Recall that the probabilities are related to the complex amplitudes of the state vector: $p_{\\uparrow} = |\\alpha_{\\uparrow}|^{2}$ and $p_{\\downarrow} = |\\alpha_{\\downarrow}|^{2}$, for a qubit state $\\begin{pmatrix} \\alpha_{\\uparrow} \\\\ \\alpha_{\\downarrow} \\end{pmatrix}$ aka $\\alpha_{\\uparrow}\\mid \\uparrow \\rangle +\\alpha_{\\uparrow}\\mid \\downarrow \\rangle$, where $\\uparrow$ and $\\downarrow$ are in the $Z$ direction. The probabilities also appear as the diagonal entries of the density matrix, which for a pure state is $\\mid \\psi \\rangle \\langle \\psi \\mid$.\n",
    "\n",
    "Such a state can be interpreted in terms of a point on the sphere, given by $(<X>, <Y>, <Z>)$, and a complex phase. The two complex orthogonal eigenstates of a Hermitian matrix (like $X$, $Y$, and $Z$) correspond to two antipodal points on the sphere.\n",
    "\n",
    "Another way of thinking about this is in terms of the <a href=\"https://en.wikipedia.org/wiki/Stereographic_projection\">stereographic projection</a>:\n",
    "\n",
    "<a href=\"https://mathematica.stackexchange.com/questions/23793/stereographic-projection\"><img style=\"margin: auto; width: 500px; display: block\" src=\"https://i.stack.imgur.com/6zeHc.jpg\"></a>\n",
    "\n",
    "We determine a point on the plane by the ratio $\\frac{\\alpha_{\\downarrow}}{\\alpha_{\\uparrow}}$: the $(x, y)$ coordinates are the real and imaginary parts of the ratio. We then project this point onto the unit sphere. Each point on the plane is mapped to a unique location on the sphere. If ${\\alpha_{\\uparrow}} = 0$, then we have an extra point, not on the plane, which corresponds to the pole of projection on the sphere. (The standard thing is to do the projection from the South Pole, despite the picture.)\n",
    "\n",
    "The point is we can think of our qubit as picking out a rotation axis *built out of expectation values*. And indeed, if the qubit is measured along the axis it's pointed in, it remains in exactly that same state. Moreover, under the action of rotations, given by complex exponentials of Paulis, this expected axis rotates just as you'd expect (although, of course, there's more sublety regarding the phase!)\n",
    "\n",
    "<hr>\n",
    "\n",
    "So, if you think about it, we just have to do a rotation before our measurement. Suppose our qubit is pointing up in the $Y^{+}$ direction, so a $Y$ measurement should give $\\uparrow$ every time. But we're only allowed $Z$ measurements. We can get the same effect as a $Y$ measurement by first doing a rotation of $90^{\\circ}$ aka $\\frac{\\pi}{2}$ around the $X$ axis, which rotates $Y^{+}$ into $Z^{+}$ (and $Y^{-}$ into $Z^{-}$). So we can turn a $Z$ measurement into a $Y$ measurement by first applying $R_{x}(\\frac{\\pi}{2}) = e^{-\\frac{\\pi}{4}iX}$ to our qubit $\\mid \\psi \\rangle$, and then measuring.\n",
    "\n",
    "Similarly, if we had to measure $\\langle X \\rangle$, we could apply $R_{y}(-\\frac{\\pi}{2}) = e^{\\frac{\\pi}{4}}iY$, rotating the opposite way a quarter turn around the $Y$ axis, to rotate $X^{+}$ into $Z^{+}$, etc."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "X+: [0.9999999999999998, 0.0, 0.0]\n",
      "rotates to: [0.0, 0.0, 1.0]\n",
      "\n",
      "Y+: [0.0, 0.9999999999999998, 0.0]\n",
      "rotates to: [0.0, 0.0, 1.0]\n"
     ]
    }
   ],
   "source": [
    "def qubit_xyz(qubit):\n",
    "    return [qt.expect(qt.sigmax(), qubit),\\\n",
    "            qt.expect(qt.sigmay(), qubit),\\\n",
    "            qt.expect(qt.sigmaz(), qubit)]\n",
    "\n",
    "# If you need <X> -> Ry\n",
    "Xp = (qt.basis(2,0)+qt.basis(2,1)).unit()\n",
    "print(\"X+: %s\" % qubit_xyz(Xp))\n",
    "print(\"rotates to: %s\" % qubit_xyz((1j*qt.sigmay()*np.pi/4).expm()*Xp))\n",
    "print()\n",
    "# If you need <Y> -> Rx\n",
    "Yp = (qt.basis(2,0)+1j*qt.basis(2,1)).unit()\n",
    "print(\"Y+: %s\" % qubit_xyz(Yp))\n",
    "print(\"rotates to: %s\" % qubit_xyz((-1j*qt.sigmax()*np.pi/4).expm()*Yp))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<X> calculated the usual way: -0.9694\n",
      "<X> estimated stochastically: -0.9696\n",
      "\n",
      "<Y> calculated the usual way: -0.2143\n",
      "<Y> estimated stochastically: -0.2100\n"
     ]
    }
   ],
   "source": [
    "print(\"<X> calculated the usual way: %.4f\" % qt.expect(qt.sigmax(), psi))\n",
    "print(\"<X> estimated stochastically: %.4f\" % expectZ((1j*qt.sigmay()*np.pi/4).expm()*psi))\n",
    "print()\n",
    "print(\"<Y> calculated the usual way: %.4f\" % qt.expect(qt.sigmay(), psi))\n",
    "print(\"<Y> estimated stochastically: %.4f\" % expectZ((-1j*qt.sigmax()*np.pi/4).expm()*psi))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Putting it together, we decompose our matrix into a sum of Pauli's, and then for each term in this sum, we do a separate experiment (state preparation, pre-measurement rotation, measurement) many times, to accumulate the expectation value of that Pauli. We add 'em all up, each with their weight, which gives us the expectation value of the original matrix on that state. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<H> calculated the usual way: 2.4516\n",
      "<H> estimated stochastically: 2.4448\n"
     ]
    }
   ],
   "source": [
    "def expect(H, psi, basis):\n",
    "    Rmap = {\"I\": qt.identity(2),\\\n",
    "            \"X\": (1j*qt.sigmay()*np.pi/4).expm(),\\\n",
    "            \"Y\": (-1j*qt.sigmax()*np.pi/4).expm(),\\\n",
    "            \"Z\": qt.identity(2)}\n",
    "    H_pauli = op_pauli(H, basis)\n",
    "    return sum([v*(expectZ(Rmap[s]*psi)\\\n",
    "                       if s.count(\"I\") != len(s)\\\n",
    "                   else 1)\\\n",
    "                if not np.isclose(v, 0) else 0\n",
    "                    for s, v in H_pauli.items()])\n",
    "\n",
    "print(\"<H> calculated the usual way: %.4f\" % qt.expect(H, psi))\n",
    "print(\"<H> estimated stochastically: %.4f\" % expect(H, psi, basis))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr>\n",
    "\n",
    "So all that's great, but what if we have multiple qubits, say $n$?\n",
    "\n",
    "Well, it turns out that strings of Paulis of length $n$ also form a basis for Hermitian matrices, in this case, $2^{n} \\times 2^{n}$ dim matrices. By a \"string,\" I mean, a tensor product of Paulis, e.g.:\n",
    "\n",
    "$XYYZIIIX$ means: $X \\otimes Y \\otimes Y \\otimes Z \\otimes I \\otimes I \\otimes I \\otimes X$.\n",
    "\n",
    "For two qubits, we have the 16 matrices:\n",
    "\n",
    "$\\begin{matrix} II & IX & IY & IZ \\\\\n",
    "                XI & XX & XY & XZ \\\\\n",
    "                YI & YX & YY & YZ \\\\\n",
    "                ZI & ZX & XY & ZZ \\end{matrix}$\n",
    "\n",
    "Thus we can decompose an $n$-qubit operator just as easily as before:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "our matrix:\n",
      "Quantum object: dims = [[2, 2], [2, 2]], shape = (4, 4), type = oper, isherm = True\n",
      "Qobj data =\n",
      "[[ 1.  0.  0.  0.]\n",
      " [ 0.  0. -1.  0.]\n",
      " [ 0. -1.  0.  0.]\n",
      " [ 0.  0.  0.  1.]]\n",
      "\n",
      "pauli decomposition:\n",
      "II: 0.50\n",
      "XX: -0.50\n",
      "YY: -0.50\n",
      "ZZ: 0.50\n",
      "\n",
      "still our matrix:\n",
      "Quantum object: dims = [[2, 2], [2, 2]], shape = (4, 4), type = oper, isherm = True\n",
      "Qobj data =\n",
      "[[ 1.  0.  0.  0.]\n",
      " [ 0.  0. -1.  0.]\n",
      " [ 0. -1.  0.  0.]\n",
      " [ 0.  0.  0.  1.]]\n"
     ]
    }
   ],
   "source": [
    "from itertools import product\n",
    "\n",
    "def construct_pauli_basis(n):\n",
    "    IXYZ = {\"I\": qt.identity(2),\\\n",
    "            \"X\": qt.sigmax(),\\\n",
    "            \"Y\": qt.sigmay(),\\\n",
    "            \"Z\": qt.sigmaz()}\n",
    "    return dict([(\"\".join(p),\\\n",
    "                  qt.tensor(*[IXYZ[s] for s in p]))\\\n",
    "                for p in product(IXYZ.keys(), repeat=n)])\n",
    "\n",
    "n = 2\n",
    "basis = construct_pauli_basis(n)\n",
    "H = qt.Qobj(np.array([[1,0,0,0],\\\n",
    "                      [0,0,-1,0],\\\n",
    "                      [0,-1,0,0],\\\n",
    "                      [0,0,0,1]]))\n",
    "H.dims = [[2]*n, [2]*n]\n",
    "H_pauli = op_pauli(H, basis)\n",
    "H2 = pauli_op(H_pauli, basis)\n",
    "\n",
    "print(\"our matrix:\")\n",
    "print(H)\n",
    "print(\"\\npauli decomposition:\")\n",
    "print_pauli(H_pauli)\n",
    "print(\"\\nstill our matrix:\")\n",
    "print(H2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In other words, we could perform the following decomposition:\n",
    "\n",
    "$H = \\begin{pmatrix}1 & 0 & 0 & 0\\\\ \n",
    "0 & 0 & -1 & 0\\\\ \n",
    "0 & -1 & 0 & 0\\\\ \n",
    "0 & 0 & 0 & 1\\end{pmatrix} = \\frac{1}{2} -\\frac{1}{2}XX -\\frac{1}{2}YY + \\frac{1}{2}ZZ $\n",
    "\n",
    "The expectation value becomes:\n",
    "\n",
    "$\\langle H \\rangle = \\frac{1}{2} -\\frac{1}{2}\\langle XX \\rangle -\\frac{1}{2} \\langle YY \\rangle + \\frac{1}{2} \\langle ZZ \\rangle$\n",
    "\n",
    "<hr>\n",
    "\n",
    "Now it is a fact that $A \\otimes B \\mid \\alpha \\rangle \\mid \\beta \\rangle = A\\mid \\alpha \\rangle \\otimes B \\mid \\beta \\rangle$ for unitaries $A$ and $B$. In other words, $A \\otimes B$ on a separable state just transforms $\\mid \\alpha \\rangle$ by $A$ and $\\mid \\beta \\rangle$ by $B$. One thing this implies is that the tensor products of $A$ and $B$'s eigenvectors are also eigenvectors of $A \\otimes B$. Moreover, the corresponding eigenvalues of $A \\otimes B$ are given by the products of $A$ and $B$ eigenvalues."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "XY eigenvalues: [1.0, -1.0, -1.0, 1.0]\n",
      "\n",
      "XY eigenvector / XY on eigenvector:\n",
      "[[ 0.5+0.j   0.5+0.j ]\n",
      " [ 0. -0.5j  0. -0.5j]\n",
      " [-0.5+0.j  -0.5+0.j ]\n",
      " [ 0. +0.5j  0. +0.5j]]\n",
      "\n",
      "XY eigenvector / XY on eigenvector:\n",
      "[[ 0.5+0.j  -0.5+0.j ]\n",
      " [ 0. +0.5j  0. -0.5j]\n",
      " [-0.5+0.j   0.5+0.j ]\n",
      " [ 0. -0.5j  0. +0.5j]]\n",
      "\n",
      "XY eigenvector / XY on eigenvector:\n",
      "[[-0.5+0.j   0.5+0.j ]\n",
      " [ 0. +0.5j  0. -0.5j]\n",
      " [-0.5+0.j   0.5+0.j ]\n",
      " [ 0. +0.5j  0. -0.5j]]\n",
      "\n",
      "XY eigenvector / XY on eigenvector:\n",
      "[[-0.5+0.j  -0.5+0.j ]\n",
      " [ 0. -0.5j  0. -0.5j]\n",
      " [-0.5+0.j  -0.5+0.j ]\n",
      " [ 0. -0.5j  0. -0.5j]]\n",
      "\n"
     ]
    }
   ],
   "source": [
    "XY = qt.tensor(qt.sigmax(), qt.sigmay())\n",
    "Xl, Xv = qt.sigmax().eigenstates()\n",
    "Yl, Yv = qt.sigmay().eigenstates()\n",
    "\n",
    "l, v = [], []\n",
    "for i in range(2):\n",
    "    for j in range(2):\n",
    "        l.append(Xl[i]*Yl[j]) # product of eigenvalues\n",
    "        v.append(qt.tensor(Xv[i],Yv[j])) # tensor product of eigenvectors\n",
    "\n",
    "print(\"XY eigenvalues: %s\\n\" % l)\n",
    "for v_ in v:\n",
    "    print(\"XY eigenvector / XY on eigenvector:\")\n",
    "    print(np.concatenate([v_.full(), (XY*v_).full()], axis=1))\n",
    "    print()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You might notice, however, that:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(array([-1., -1.,  1.,  1.]),\n",
       " array([Quantum object: dims = [[2, 2], [1, 1]], shape = (4, 1), type = ket\n",
       " Qobj data =\n",
       " [[0.70710678+0.j        ]\n",
       "  [0.        +0.j        ]\n",
       "  [0.        +0.j        ]\n",
       "  [0.        -0.70710678j]],\n",
       "        Quantum object: dims = [[2, 2], [1, 1]], shape = (4, 1), type = ket\n",
       " Qobj data =\n",
       " [[0.        +0.j        ]\n",
       "  [0.        -0.70710678j]\n",
       "  [0.70710678+0.j        ]\n",
       "  [0.        +0.j        ]],\n",
       "        Quantum object: dims = [[2, 2], [1, 1]], shape = (4, 1), type = ket\n",
       " Qobj data =\n",
       " [[-0.70710678+0.j        ]\n",
       "  [ 0.        +0.j        ]\n",
       "  [ 0.        +0.j        ]\n",
       "  [ 0.        -0.70710678j]],\n",
       "        Quantum object: dims = [[2, 2], [1, 1]], shape = (4, 1), type = ket\n",
       " Qobj data =\n",
       " [[0.        +0.j        ]\n",
       "  [0.        +0.70710678j]\n",
       "  [0.70710678+0.j        ]\n",
       "  [0.        +0.j        ]]], dtype=object))"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "XY.eigenstates()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "gives different eigenstates!\n",
    "\n",
    "Yet:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "XY eigvec (1.00) / in alternate basis:\n",
      "[[ 0.5       +0.j   0.        +0.j ]\n",
      " [ 0.        -0.5j  0.        +0.j ]\n",
      " [-0.5       +0.j  -0.70710678+0.j ]\n",
      " [ 0.        +0.5j -0.70710678+0.j ]]\n",
      "\n",
      "XY eigvec (-1.00) / in alternate basis:\n",
      "[[ 0.5       +0.j   0.70710678+0.j ]\n",
      " [ 0.        +0.5j -0.70710678+0.j ]\n",
      " [-0.5       +0.j   0.        +0.j ]\n",
      " [ 0.        -0.5j  0.        +0.j ]]\n",
      "\n",
      "XY eigvec (-1.00) / in alternate basis:\n",
      "[[-0.5       +0.j  -0.70710678+0.j ]\n",
      " [ 0.        +0.5j -0.70710678+0.j ]\n",
      " [-0.5       +0.j   0.        +0.j ]\n",
      " [ 0.        +0.5j  0.        +0.j ]]\n",
      "\n",
      "XY eigvec (1.00) / in alternate basis:\n",
      "[[-0.5       +0.j   0.        +0.j ]\n",
      " [ 0.        -0.5j  0.        +0.j ]\n",
      " [-0.5       +0.j   0.70710678+0.j ]\n",
      " [ 0.        -0.5j -0.70710678+0.j ]]\n",
      "\n"
     ]
    }
   ],
   "source": [
    "XYl, XYv = XY.eigenstates()\n",
    "for i, v_ in enumerate(v):\n",
    "    print(\"XY eigvec (%.2f) / in alternate basis:\" % (l[i]))\n",
    "    print(np.concatenate([v_.full(), (v_.transform(XYv)).full()], axis=1))\n",
    "    print()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the alternate basis, the components are associated with the eigenvalues in order $\\{-1, -1, 1, 1\\}$.\n",
    "\n",
    "Each of the tensor product eigenstates, in the basis provided by the eigenstates calculated by qutip, lives either in the $1$ or $-1$ sector: in other words, can be written as a sum of eigenstates in the other basis with the same eigenvalue, indeed, the same eigenvalue associated with the tensor product state. So the ambiguity is due to the degeneracy of the eigenvalues. \n",
    "\n",
    "In any case, we're good to go!\n",
    "\n",
    "<hr>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we have multiple qubits, for each Pauli string in the decomposition of our operator, e.g., $XYY$, we want to follow the same rule as before, just on each of the qubits. Recall:\n",
    "\n",
    "$X \\rightarrow R_{y}(-\\frac{\\pi}{2}) = e^{\\frac{\\pi}{4}}iY$\n",
    "\n",
    "$Y \\rightarrow R_{x}(\\frac{\\pi}{2}) = e^{-\\frac{\\pi}{4}iX}$\n",
    "\n",
    "So for $XYY$,  we'd do a $Z$ measurement on each of the qubits in turn, but first: a rotation $R_{y}(-\\frac{\\pi}{2})$ on the first qubit, and a rotation $R_{x}(\\frac{\\pi}{2})$ on the second two.\n",
    "\n",
    "Upon measurement, we get three bits, e.g., $1, -1, 1$, representing the results of a $Z$ measurement on each of the three qubits. We multiply these three eigenvalues together, to get the eigenvalue of the whole state: $-1$. We add it to the pile, and then repeat the experiment thousands of times, eventually dividing by the number of experiments. And then we have the expectation value of that Pauli operator, which then is weighted in the sum which gives the expectation value of the original operator."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<H> calculated the usual way: 0.4919\n",
      "<H> estimated stochastically: 0.4873\n"
     ]
    }
   ],
   "source": [
    "# Tensors an operator appropriately with the identity so that it acts\n",
    "# on the i-th of n subsystems.\n",
    "def upgrade(op, i, n):\n",
    "    return qt.tensor(*[op if i == j \\\n",
    "                          else qt.identity(op.shape[0])\\\n",
    "                                for j in range(n)])\n",
    "\n",
    "# Expectation value of Z measurements on n qubits.\n",
    "# If stochastic == False, we simply calculate the expectation value \n",
    "# of the ...ZZZ... operator on the state in the usual way. \n",
    "# Otherwise: For n_shots, we iterate over the qubits, calculating \n",
    "# the density matrix of each qubit by a partial trace,\n",
    "# extracting the probabilities for a Z measurement  \n",
    "# from its diagonal entries. We flip a coin and randomly select\n",
    "# an eigenvalue {1, -1}. We project that qubit into the\n",
    "# corresponding Z eigenstate, and move onto the next.\n",
    "# Once we've gone through all the qubits, we multiply\n",
    "# all the eigenvalues together to get the eigenvalue\n",
    "# of the whole state, and add that to our running total,\n",
    "# which, after n_shots, we divide by n_shots to get our\n",
    "# expectation value.\n",
    "def expectZ(qubits, n_shots=1024, stochastic=True):\n",
    "    n = len(qubits.dims[0])\n",
    "    if not stochastic:\n",
    "        return qt.expect(qt.tensor(*[qt.sigmaz()]*n), qubits)\n",
    "    running = 0\n",
    "    Zp = [v*v.dag() for v in qt.sigmaz().eigenstates()[1]]\n",
    "    for i in range(n_shots):\n",
    "        choices = []\n",
    "        state = qubits.copy()\n",
    "        for j in range(n):\n",
    "            dm = state.ptrace(j)\n",
    "            choice = np.random.choice([1, -1], p=[dm[0,0].real, dm[1,1].real])\n",
    "            choices.append(choice)\n",
    "            state = (upgrade(Zp[0 if choice == -1 else 1], j, n)*state).unit()\n",
    "        running += np.prod(choices)\n",
    "    running = running/n_shots\n",
    "    return running\n",
    "\n",
    "def expect(H, psi, basis, stochastic=True):\n",
    "    Rmap = {\"I\": qt.identity(2),\\\n",
    "            \"X\": (1j*qt.sigmay()*np.pi/4).expm(),\\\n",
    "            \"Y\": (-1j*qt.sigmax()*np.pi/4).expm(),\\\n",
    "            \"Z\": qt.identity(2)}\n",
    "    H_pauli = op_pauli(H, basis)\n",
    "    return sum([v*(expectZ(\\\n",
    "                        qt.tensor(*[Rmap[o] for o in s])*psi,\\\n",
    "                        stochastic=stochastic)\\\n",
    "                            if s.count(\"I\") != len(s)\\\n",
    "                   else 1)\\\n",
    "                if not np.isclose(v, 0) else 0\n",
    "                    for s, v in H_pauli.items()])\n",
    "\n",
    "psi = qt.rand_ket(4)\n",
    "psi.dims = [[2,2],[1,1]]\n",
    "print(\"<H> calculated the usual way: %.4f\" % qt.expect(H, psi))\n",
    "print(\"<H> estimated stochastically: %.4f\" % expect(H, psi, basis))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr>\n",
    "\n",
    "So we can take expectation values to our heart's content. Now we turn to issue of parameterizing the state.\n",
    "\n",
    "As before, we're working with:\n",
    "\n",
    "$H = \\begin{pmatrix}1 & 0 & 0 & 0\\\\ \n",
    "0 & 0 & -1 & 0\\\\ \n",
    "0 & -1 & 0 & 0\\\\ \n",
    "0 & 0 & 0 & 1\\end{pmatrix} = \\frac{1}{2} -\\frac{1}{2}XX -\\frac{1}{2}YY + \\frac{1}{2}ZZ $\n",
    "\n",
    "We'll take as our ansatz for the state: \n",
    "\n",
    "$\\mid \\psi(\\theta) \\rangle = (R_{x}(\\theta) \\otimes I)(CNOT)(h \\otimes I)\\mid \\uparrow \\uparrow \\rangle$\n",
    "\n",
    "Here, $h = \\frac{1}{\\sqrt{2}}\\begin{pmatrix} 1 & 1 \\\\ 1 & -1 \\end{pmatrix}$, the Hadamard gate, aka the Fourier Transform operator ($h^{\\dagger}Xh = Z$).\n",
    "\n",
    "$CNOT = \\begin{pmatrix} 1 & 0 & 0 & 0 \\\\ 0 & 1 & 0 & 0 \\\\ 0 & 0 & 0 & 1 \\\\ 0 & 0 & 1 & 0 \\end{pmatrix}$ is the controlled-$X$ or controlled-$NOT$ gate, which flips the second qubit if the first qubit is $\\mid \\downarrow \\rangle$. Acting on $\\mid \\uparrow \\uparrow \\rangle$, $h$ takes the first qubit to an equal superposition of $\\uparrow$ and $\\downarrow$: in other words, from $Z^{+}$ to $X^{+}$. Then the $CNOT$ conditions the flippage of the second qubit on the first qubit which is in a superposition. This leads to the maximally entangled state: $\\frac{1}{\\sqrt{2}} (\\mid \\uparrow \\uparrow \\rangle + \\mid \\downarrow \\downarrow \\rangle)$. Finally, we perform a rotation of the first qubit around the $X$ axis by some angle $\\theta$, which allows us to parameterize all manner of maximally entangled states, from correlated to anticorrelated.\n",
    "\n",
    "Then we can then iterate over values of $\\theta$ to try and find our ground state.\n",
    "\n",
    "<hr>\n",
    "\n",
    "Putting it all together:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "our matrix:\n",
      "Quantum object: dims = [[2, 2], [2, 2]], shape = (4, 4), type = oper, isherm = True\n",
      "Qobj data =\n",
      "[[ 1.  0.  0.  0.]\n",
      " [ 0.  0. -1.  0.]\n",
      " [ 0. -1.  0.  0.]\n",
      " [ 0.  0.  0.  1.]]\n",
      "\n",
      "compare:\n",
      "lowest eigval =  -1.0000\n",
      "eigenvector = \n",
      "Quantum object: dims = [[2, 2], [1, 1]], shape = (4, 1), type = ket\n",
      "Qobj data =\n",
      "[[0.        ]\n",
      " [0.70710678]\n",
      " [0.70710678]\n",
      " [0.        ]]\n",
      "\n",
      "theta = 0.0000, <H> = 1.0000\n",
      "theta = 0.3307, <H> = 0.9458\n",
      "theta = 0.6614, <H> = 0.7891\n",
      "theta = 0.9921, <H> = 0.5469\n",
      "theta = 1.3228, <H> = 0.2455\n",
      "theta = 1.6535, <H> = -0.0826\n",
      "theta = 1.9842, <H> = -0.4017\n",
      "theta = 2.3149, <H> = -0.6773\n",
      "theta = 2.6456, <H> = -0.8795\n",
      "theta = 2.9762, <H> = -0.9864\n",
      "theta = 3.3069, <H> = -0.9864\n",
      "theta = 3.6376, <H> = -0.8795\n",
      "theta = 3.9683, <H> = -0.6773\n",
      "theta = 4.2990, <H> = -0.4017\n",
      "theta = 4.6297, <H> = -0.0826\n",
      "theta = 4.9604, <H> = 0.2455\n",
      "theta = 5.2911, <H> = 0.5469\n",
      "theta = 5.6218, <H> = 0.7891\n",
      "theta = 5.9525, <H> = 0.9458\n",
      "theta = 6.2832, <H> = 1.0000\n",
      "\n",
      "final theta = 3.3069, <H> = -0.9864\n",
      "final vector =\n",
      "Quantum object: dims = [[2, 2], [1, 1]], shape = (4, 1), type = ket\n",
      "Qobj data =\n",
      "[[-0.05839242+0.j        ]\n",
      " [ 0.        -0.70469165j]\n",
      " [ 0.        -0.70469165j]\n",
      " [-0.05839242+0.j        ]]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, '<H>')"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgQAAAF3CAYAAADaXFNkAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAABCP0lEQVR4nO3dd3zV9d3+8dc7ixAIIxAgkEDCHgoIEWQpDtyjrjpaB6hUra0dtrcdd62977u3v7Z22KpVERStVXHUhXsyFAjKniGDJEAIK4RAyPr8/uDYO6UBE0jO54zr+XicxznnO865zmHkynd8vuacQ0RERKJbjO8AIiIi4p8KgYiIiKgQiIiIiAqBiIiIoEIgIiIiqBCIiIgIEOc7gE9du3Z1mZmZvmOIiIgExdKlS3c451IbmxfVhSAzM5OcnBzfMURERILCzAqPNE+7DERERESFQERERFQIREREBBUCERERQYVAREREUCEQERERVAhEREQEFQIRERFBhUBEREQIsUJgZjPNbLuZrTrCfDOzB8ws18xWmNmoBvPONbP1gXl3By+1iIhI+AupQgA8AZx7lPnnAQMCt+nAwwBmFgs8GJg/FLjGzIa2alIREZEIElKFwDn3CbDrKItcAsx2h3wGdDKzNGAMkOucy3POVQPPBpYNmoIdlSzK24lzLphvKyIiEWrzzv18umln0N4v3C5u1AsoavC8ODCtseljG3sBM5vOoa0L9O7du8WCPflpAbMWFNCnSxJXjErn8tHp9OzUtsVeX0REIt/+6lrmrtzGnJwiFuXvok+XJD66azJm1urvHW6FoLFvxB1l+r9PdO5R4FGA7OzsFvt1/kfnDOKEnh2Zs7SI+9/dwO/f28DE/l25YnQ65wzrQWJ8bEu9lYiIRBDnHDmFu5mTU8QbK7ZSWV1HZpck7jp7IJeNSg9KGYDwKwTFQEaD5+nAFiDhCNODJikhjstHH9oysHnnfl74vJgXlxZz57PLSE6M4+IRPbkyO4MR6R2D9ocrIiKha2v5AV76vIQXlhaTv6OSpIRYLjgxjSuzMzg5s3PQf1ZYqO3zNrNM4HXn3AmNzLsAuAM4n0O7BB5wzo0xszhgA3AmUAIsAa51zq0+2ntlZ2e7nJycFv4E/6e+3vFp3k7m5BTx5qptHKytZ0C39lyZnc7XTupFt+TEVntvEREJPVU1dby7ppQ5S4uZv7GMegdjs1K4MjuD807oQbs2rft7upktdc5lNzovlAqBmf0dmAx0BUqBe4B4AOfcX+1QXfoLh85E2A9Mdc7lBNY9H/gjEAvMdM79z1e9X2sXgob2VtXwxoqtzMkp4vPNe4iNMSYPTOXK7HTOGNydhLiQOr5TRERaiHOOFcXlzFlaxKvLtrC3qpZendpy+aheXD46nT5d2gUtS9gUgmALZiFoKHf7Pl5YWsxLnxezveIgKe0SuGRkT64cncHQnh2CnkdERFpeWcVB/vFFCXOWFrGhdB9t4mI494QeXDk6g/H9uhATE/zdxyoER+CrEHyptq6eeRt3MGdpEe+t2U51XT3DenbgytHpXDKyF53bJXjLJiIizVddW8+H67czJ6eYj9Zvp7beMTKjE1dmp3Ph8J50bBvvNZ8KwRH4LgQN7a6s5pVlJcxZWszqLXtJiI3hrKHduHJ0BpMHpepARBGRELaxtIJnlxTxjy9K2FlZTWpyGy4b1YsrR6fTv1uy73j/pEJwBKFUCBpas2Uvc5YW8cqyLeyqrObG8Zncc9FQlQIRkRD0wbpSps9eihmcNaQ7V2anc+qAVOJiQ+/YsKMVgnA77TAqDO3ZgXt6DuMn5w3hf99cy6wFBaQmt+Hbp/f3HU1ERBpYWriL2//2OYPTknli6hi6tm/jO9IxUyEIYQlxMfznBUPZs7+G3769npR2CVwzpuVGVxQRkWO3flsFU2ctIa1j27AvA6BCEPJiYozfXDGc3fur+dnLK+mcFM+5J6T5jiUiEtWKd+/n+pmLSIyPZfa08C8DEGIXN5LGxcfG8NA3RjEyoxPf/fsyFm7a4TuSiEjU2rnvINc/vpgD1XXMvmkMGSlJviO1CBWCMJGUEMfMG0+mT5ckps9eyqqSct+RRESizr6DtUx9Ygklew7w+I0nM7hH5Iwdo0IQRjolJTD7pjF0SIzjxlmLKdhR6TuSiEjUOFhbx61PLWX1lr08eO0oTs5M8R2pRakQhJm0jm2ZfdNY6uod189czPa9Vb4jiYhEvLp6xw+fX8783B3cd9mJnDW0u+9ILU6FIAz179aeWVPHsGPfQW6YtYTyAzW+I4mIRCznHPe+tprXV2zlJ+cN5srsjK9eKQypEISpkRmdeOS60eRur+CW2TlU1dT5jiQiEpH+/EEusz8tZPqpffnWaf18x2k1KgRhbNKAVO7/+kiWFOziO3//gtq6et+RREQiytOfFfL7dzdw2ahe3H3uYN9xWpUKQZi7eERP7rlwKO+uKeVnL68imoeiFhFpSXNXbuU/X1nFGYO78f8uH+7l6oTBpIGJIsCNE7LYVVnNAx/k0qV9Aj+O8BYrItLaFubu4HvPLmNU7848eO0o4kPwugQtTYUgQnx/ykB2VFbz0EebSGmXwM2T+vqOJCISllYWl3PL7ByyurZj5g0n0zYh1nekoFAhiBBmxn9dcgK7K6v57zfW0qV9ApeelO47lohIWMnfUcmNsxbTKSmBJ6eNoWNSvO9IQRP520CiSGyM8YerRjKubxd+NGcFH67f7juSiEjYKN1bxXWPL8IBs28aQ4+Oib4jBZUKQYRJjI/l0etHM6hHMrc9vZSlhbt9RxIRCXnl+2u4/vHF7K6s5ompJ9Mvtb3vSEGnQhCBkhPjeWLqGHp0SGTaE0vYUFrhO5KISMiqqqnj5tlLyNuxj0euy2Z4eiffkbxQIYhQqclteOqmsSTExXD944sp2XPAdyQRkZBTW1fPHc98Tk7hbv5w1UgmDujqO5I3KgQRLCMlidnTxlBZXct1jy9iV2W170giIiHDOcfdL63kvbXb+dXFw7hweE/fkbxSIYhwQ9I68PgNJ1Oy+wBTZy2m8mCt70giIiHhvrfW8cLSYu48cwDXjcv0Hcc7FYIoMCYrhb9cO4pVW/Zy69NLqa7VEMciEt0e+ySPRz7O45un9OZ7Zw3wHSckqBBEiSlDu/O/l57IvI07+OGc5dTXa4hjEYlOLy4t5n/mruX8E3tw78UnYBbZQxI3lQYmiiJfPzmDnZXV/L+31tGrU1vuPk9DHItIdPl0005+/OIKxvfrwh+uGklshF+foDm0hSDK3HpaX67KzuDRTzaxdute33FERIKmuraen/1jJemd2/Lo9dm0iYuOIYmbSoUgypgZd583mA5t47nn1dW6OqKIRI1ZC/LJK6vknouG0r6NNpAfLqQKgZmda2brzSzXzO5uZP6PzGxZ4LbKzOrMLCUwr8DMVgbm5QQ/ffjo3C6Bu84exOL8Xby2YqvvOCIira50bxUPvL+RMwZ344zB3X3HCUkhUwjMLBZ4EDgPGApcY2ZDGy7jnPutc26kc24k8BPgY+fcrgaLnB6Ynx2s3OHqmjG9GdazA79+Y61ORRSRiHffm+uoqXP84sKhX71wlAqZQgCMAXKdc3nOuWrgWeCSoyx/DfD3oCSLQLExxr0XD2Pb3ioe/DDXdxwRkVazpGAXL39Rwi2nZpHZtZ3vOCErlApBL6CowfPiwLR/Y2ZJwLnAiw0mO+AdM1tqZtNbLWUEyc5M4bKTevHYvDzyd1T6jiMi0uLq6h33vLKatI6JfPv0/r7jhLRQKgSNnftxpCPeLgIWHLa7YIJzbhSHdjl828xObfRNzKabWY6Z5ZSVlR1f4ghw93mDaRMXy69eW+07iohIi3tm8WbWbN3LT88fQlKCDiQ8mlAqBMVARoPn6cCWIyx7NYftLnDObQncbwde5tAuiH/jnHvUOZftnMtOTU097tDhrluHRO48cwAfri/j/bWlvuOIiLSY3ZXV3P/Oek7pm8KFw9N8xwl5oVQIlgADzCzLzBI49EP/1cMXMrOOwGnAKw2mtTOz5C8fA2cDq4KSOgLcMD6Tfqnt+NXra6iqqfMdR0SkRfzunfVUVNVqNMImCplC4JyrBe4A3gbWAs8751ab2a1mdmuDRS8F3nHONdzp3R2Yb2bLgcXAG865t4KVPdwlxMXwy4uHUbhzPzPm5fmOIyJy3FaVlPPM4s1cd0ofBvVI9h0nLFg0D0yTnZ3tcnI0ZMGXvvVUDp9s2MH7PzyNnp3a+o4jInJMnHNc8ddPKdhRyQd3TaZj23jfkUKGmS090qn5IbOFQPz7+QVDqXeO/5m71ncUEZFj9vIXJSwt3M1/nDtYZaAZVAjknzJSkrhtcj/eWLGVhZt2+I4jItJsFVU1/O+b6xiR0YkrRqf7jhNWVAjkX9x6Wj/SO7fl3lfXUFNX7zuOiEiz/PmDXMoqDnLvxcOI0ZUMm0WFQP5FYnwsP79gKOtLK3jq00LfcUREmix3ewUz5+fz9ex0RmZ08h0n7KgQyL85Z1h3Jg3oyh/e28COfQd9xxER+UrOOe59bQ1tE2L58bmDfccJSyoE8m/MjHsuGsaB6jp+89Y633FERL7S26tLmbdxBz+YMpCu7dv4jhOWVAikUf27tWfaxCyezylmWdEe33FERI6oqqaO/3p9DYO6J3PdKX18xwlbKgRyRN85oz/dkttwzyurqK+P3vEqRCS0PfzRJkr2HOCXFw8jLlY/1o6Vvjk5ouTEeH5y/mCWF5czZ2nRV68gIhJkRbv289ePN3Hh8DTG9eviO05YUyGQo/rayF5k9+nMb95aT/mBGt9xRET+xX+/sYYYM352wRDfUcKeCoEclZnxy4uHsWt/NX94d4PvOCIi//TJhjLeXl3KHWf0J62jhls/XioE8pVO6NWRb4ztzVOfFbJu217fcUREqK6t55evrSazSxI3T8ryHSciqBBIk/xwyiCSE+O455XVRPMFsUQkNDyxMJ+8skp+cdFQ2sTF+o4TEVQIpEk6t0vgrrMHsSh/F6+v2Oo7johEse17q/jTexs5Y3A3zhjc3XeciKFCIE12zZjeDOvZgV/PXUvlwVrfcUQkSt335jpq6hy/uHCo7ygRRYVAmiw2xrj34mFsLa/iwQ9zfccRkSiUU7CLl74o4ZZTs8js2s53nIiiQiDNkp2ZwmUn9WLGvHwKdlT6jiMiUaSu3vGLV1aT1jGRb5/e33eciKNCIM1293mDSYiL4Vevr/EdRUSiyN8Xb2bN1r389PwhJCXE+Y4TcVQIpNm6dUjkzjMH8MG67by/ttR3HBGJArsrq/ndO+s5pW8KFw5P8x0nIqkQyDG5YXwm/VLb8avX11BVU+c7johEuN+9s56KqlruvfgEzMx3nIikQiDHJCEuhl9ePIzCnft5fH6+7zgiEsFWlZTzzOLNXHdKHwb1SPYdJ2KpEMgxmzQglXOGdecvH+SyZc8B33FEJAI55/jlq6tJSUrg+1MG+o4T0VQI5Lj8/IKh1DvHr+eu9R1FRCLQP5aVkFO4m/84dzAd28b7jhPRVAjkuGSkJHHb5H68vmIrCzft8B1HRCJIRVUNv567jhEZnbhidLrvOBFPhUCO262n9SO9c1vufXUNdfW6zoGItIy/fJBLWcVB7r14GDExOpCwtakQyHFLjI/lJ+cNYX1pBXNX6joHInL8duw7yJOfFnDZSb0YmdHJd5yooEIgLeK8E3rQL7UdD320SVdDFJHjNmtBPgdr6/n2GRqRMFhUCKRFxMQYt03uz9qte/lofZnvOCISxvZW1TD708LALxrtfceJGioE0mIuGdmTXp3a6sJHInJcnv6skIqqWm6frK0DwRRShcDMzjWz9WaWa2Z3NzJ/spmVm9mywO0XTV1XWl98bAzTT+1LTuFuFufv8h1HRMJQVU0dM+fnc+rAVE7o1dF3nKgSMoXAzGKBB4HzgKHANWbW2MWu5znnRgZuv2rmutLKvp6dQZd2CdpKICLH5PmcInbsq+b2yf18R4k6IVMIgDFArnMuzzlXDTwLXBKEdaUFtU2IZdrELD7eUMaqknLfcUQkjNTU1fPIx3mM7tOZsVkpvuNEnVAqBL2AogbPiwPTDjfOzJab2ZtmNqyZ62Jm080sx8xyysp08FtruG5cH5LbxPHwR5t8RxGRMPLqsi2U7DnA7ZP76QJGHoRSIWjsT//w89c+B/o450YAfwb+0Yx1D0107lHnXLZzLjs1NfVYs8pRdEiM57pxfZi7ait5Zft8xxGRMFBf73j4400M7pHMGYO7+Y4TlUKpEBQDGQ2epwNbGi7gnNvrnNsXeDwXiDezrk1ZV4Jr2sQsEmJj+OvH2kogIl/tnTWl5G7fx23aOuBNKBWCJcAAM8syswTgauDVhguYWQ8L/E0xszEcyr+zKetKcHVt34arT87g5S9KdCVEETkq5xwPf5RLny5JXHBimu84UStkCoFzrha4A3gbWAs875xbbWa3mtmtgcWuAFaZ2XLgAeBqd0ij6wb/U0hDt5zaF+fgsXl5vqOISAhbkLuT5cXlfOvUfsTFhsyPpagT5ztAQ4HdAHMPm/bXBo//AvylqeuKX+mdk7hkZC+eXVzEHaf3p0v7Nr4jiUgIeuijXLolt+Hy0Y0eCy5Boiomreq2yX2pqq3jiYUFvqOISAj6YvNuFm7ayS2T+tImLtZ3nKimQiCtqn+3ZM4Z2oMnFxZQUVXjO46IhJiHPtpEx7bxXDu2t+8oUU+FQFrd7af3Y29VLX9btNl3FBEJIeu3VfDumlJuHJ9JuzYhtQc7KqkQSKsbnt6JSQO6MmNePlU1db7jiEiI+OvHm0hKiOXG8Zm+owgqBBIkt03ux459B5mztNh3FBEJAUW79vPq8i1cO6Y3ndsl+I4jqBBIkIzr24WTenfikY83UVtX7zuOiHj2yCebiDXj5kl9fUeRABUCCQoz4/bJ/SnefYDXVmgQSZFotr2iiudzirl8dDo9Oib6jiMBKgQSNGcO7sag7sk8/NEm6usbvdSEiESBx+fnU1tXz62naetAKFEhkKCJiTFum9yPDaX7eG9tqe84IuJB+f4anv60kAuH96RPl3a+40gDKgQSVBcOTyMjpS0PfrQJ57SVQCTazP60gMrqOm6b3M93FDmMCoEEVVxsDN86tR/Li/bw6aadvuOISBDtr65l5oJ8zhzcjSFpHXzHkcOoEEjQXTE6ndTkNjz0kS6NLBJNnl1cxO79Ndx+urYOhCIVAgm6xPhYbp6YxfzcHSwv2uM7jogEQXVtPY/Ny2NsVgqj+6T4jiONUCEQL75xSh86JMbx0Ee5vqOISBD844sStpZXcfvp/X1HkSNQIRAv2reJ48bxmby9upSNpRW+44hIK6qrdzz88SZO6NWBUwd09R1HjkCFQLy5cUIWbeNjefhjHUsgEsneWrWN/B2V3D65P2bmO44cgQqBeJPSLoFrxvTmlWVbKNq133ccEWkFzjke/DCXvqntOGdYD99x5ChUCMSrW07NIsbgsXl5vqOISCv4eEMZa7bu5dbT+hEbo60DoUyFQLxK69iWy05K57klRZRVHPQdR0Ra2EMfbqJnx0S+NrKX7yjyFVQIxLtvndaX6rp6Zi7I9x1FRFrQkoJdLC7YxS2n9iUhTj9uQp3+hMS7vqntOf/ENJ76tJDyAzW+44hIC3now1xS2iVw9cm9fUeRJlAhkJBw22n92Hewlqc/K/QdRURawOot5Xy4voxpEzJpmxDrO440gQqBhIQTenVk8qBUHp+fz4HqOt9xROQ4PfzRJtq3ieO6cZm+o0gTqRBIyLh9cn92VVbz3JLNvqOIyHHI31HJ3JVb+eYpfejYNt53HGkiFQIJGWOyUjg5szOPfpJHdW297zgicowe+XgT8bEx3DQxy3cUaQYVAgkpt0/uz5byKl5ZVuI7iogcg23lVbz4eTFfz84gNbmN7zjSDCoEElImD0plSFoHHv54E3X1znccEWmmx+blUe9g+ql9fUeRZlIhkJBiZtw+uR95ZZW8vXqb7zgi0gy7Kqt5ZtFmLhnRk4yUJN9xpJlCqhCY2blmtt7Mcs3s7kbmf8PMVgRuC81sRIN5BWa20syWmVlOcJNLSzr/xDQyuyTx0Ee5OKetBCLh4omFBRyoqeO2yf18R5FjEDKFwMxigQeB84ChwDVmNvSwxfKB05xzw4H/Ah49bP7pzrmRzrnsVg8srSY2xrj1tH6sKtnLJxt3+I4jIk2w72AtTyzI5+yh3RnQPdl3HDkGIVMIgDFArnMuzzlXDTwLXNJwAefcQufc7sDTz4D0IGeUILl0VC96dEjkoQ9zfUcRkSZ4ZlEhe6tquf30/r6jyDEKpULQCyhq8Lw4MO1IbgLebPDcAe+Y2VIzm36klcxsupnlmFlOWVnZcQWW1tMmLpabJ2WxKH8XSwt3+Y4jIkdxsLaOx+blM6F/F0ZmdPIdR45RKBWCxq6L2egOZDM7nUOF4D8aTJ7gnBvFoV0O3zazUxtb1zn3qHMu2zmXnZqaeryZpRVdO7Y3HdvGM2OeLnokEspeXbaFsoqD3Hqajh0IZ6FUCIqBjAbP04Ethy9kZsOBGcAlzrmdX053zm0J3G8HXubQLggJY0kJcVw9JoO3V2+jaNd+33FEpBHOOWYuKGBQ92Qm9u/qO44ch1AqBEuAAWaWZWYJwNXAqw0XMLPewEvAdc65DQ2mtzOz5C8fA2cDq4KWXFrN9eMyMTOe0kWPRELSovxdrN26l6kTDv1blfAVMoXAOVcL3AG8DawFnnfOrTazW83s1sBivwC6AA8ddnphd2C+mS0HFgNvOOfeCvJHkFbQq1Nbzj2hB39fvJnKg7W+44jIYWbOz6dzUjxfO+loh3xJOIjzHaAh59xcYO5h0/7a4PHNwM2NrJcHjDh8ukSGaRMyeWPFVl78vJjrdeU0kZCxeed+3l1byu2T+5EYr0sch7uQ2UIgciSjendmRHpHnlhQQL2GMxYJGU9+WkCsGdedkuk7irQAFQIJeWbGtIlZ5O2o5OMNOlVUJBRUVNXw3JIiLhieRo+Oib7jSAtQIZCwcN4JaXRLbsPMBToFUSQUvLC0mH0Ha5k6QZc4jhQqBBIWEuJiuH5cH+Zt3MGG0grfcUSiWn2944mFBYzq3UkDEUUQFQIJG9eM6U2buBhmLSjwHUUkqn2wbjuFO/czbaK2DkQSFQIJG13at+FrI3vx8hfF7K6s9h1HJGrNWphPWsdEzhnWw3cUaUEqBBJWpk7MpKqmnr8v2ew7ikhUWrdtLwtyd3L9uEziY/UjJJLoT1PCyuAeHZjQvwuzFxZSU1fvO45I1Jk1v4DE+BiuGZPx1QtLWFEhkLAzdXwW2/ZW8daqbb6jiESVnfsO8vKyEi4blU6npATfcaSFqRBI2DljcDf6dEnSKYgiQfb3xZuprq1n6vhM31GkFagQSNiJiTGmjs/ki817+GLzbt9xRKJCdW09sz8t5NSBqQzonuw7jrQCFQIJS1dkZ5DcJk6nIIoEyZurtrK94iBTJ2T6jiKtRIVAwlL7NnF8/eQM5q7cytbyA77jiEQ05xyPz8+nb2o7ThuQ6juOtBIVAglbN47PpN45nvq00HcUkYj2+ebdrCguZ+qELGJizHccaSUqBBK2MlKSOGtId/6+eDMHqut8xxGJWDMXFNAhMY7LR/XyHUVakQqBhLVpE7PYvb+Gfywr8R1FJCKV7DnAW6u2cc2Y3iQlxPmOI61IhUDC2tisFIamdWDWgnycc77jiESc2Z8WAHC9TjWMeCoEEtbMjKkTMtlQuo8FuTt9xxGJKPura3l2cRHnDOtOr05tfceRVnZMhcDMUszsF2amK1uIdxeN6EnX9gkaqEikhb30eQnlB2qYNkFXNYwGx7qF4BvAPcCNLRdF5NgkxsfyjbF9+GDddvLK9vmOIxIR6usdsxbkMzy9I6P7dPYdR4LgWAvB9cCiwL2Id984pTfxscaTCwt8RxGJCPNyd7CprJKpEzIx06mG0aDZhcDMhgLDgWuBDDM7ucVTiTRTt+RELhrRkzlLiyk/UOM7jkjYmzk/n9TkNlxwYk/fUSRIjmULwQ3AW865AuAfaLeBhIhpE7LYX13H80uKfEcRCWu52/fx8YYyrj+lDwlxOvY8WjTrT9rMYjh0/MBTgUlPA1eZWXxLBxNprhN6dWRMZgpPLCygtq7edxyRsPXEwnwS4mK4dmxv31EkiJpb/c4CkoBXA8/fBaqBC1sylMixmjYxk5I9B3hvbanvKCJhac/+al5cWsLXRvakS/s2vuNIEDW3EFwPzHHOVQM45+qBZ9BuAwkRU4b2IL1zW2bOL/AdRSQsPbukiAM1dUzVqYZRp8mFwMySgUv5v90FX/obcJ6ZdW3JYCLHIjbGuGFcJosLdrGqpNx3HJGwUltXz+yFBYzr24UhaR18x5Ega84WghjgPOfc/IYTnXNfAGcAurqMhISvn5xBUkKsBioSaaa3V5eypbyKaRO1dSAaNbkQOOfKnXOfHGHefOfc7uMNY2bnmtl6M8s1s7sbmW9m9kBg/gozG9XUdSV6dGwbz5Wj03lt+Ra2V1T5jiMSNmYuyKdPlyTOGNzNdxTxIGTOJzGzWOBB4DxgKHBNYMyDhs4DBgRu04GHm7GuRJEbxmdSU+f422ebfUcRCQvLi/awtHA3N4zLJDZGAxFFoyZdy9LM8oGvvJScc67vcWQZA+Q65/IC7/kscAmwpsEylwCz3aHL2n1mZp3MLA3IbMK6EkX6prbnjMHd+NuiQm6b3I/E+FjfkURC2qwF+bRvE8eV2em+o4gnTb249V8aPDbgv4HfAy15ebleQMMRZYqBsU1YplcT15UoM21CFt98fBGvLd/CldkZvuOIhKzSvVW8vmIr14/LJDlRw8pEqyYVAufc/Q2fm9k9wIwvfyNvIY1tozp8q8SRlmnKuodewGw6h3Y30Lu3Bt2IZBP6d2Fg9/bMWlDAFaPTNR67yBE8/Vkhdc5x4/hM31HEo5A5hoBDv9U3/DUuHdjSxGWasi4AzrlHnXPZzrns1NTU4w4tocvMmDohizVb97Iof5fvOCIhqaqmjr8t2sxZQ7rTu0uS7zjiUSgVgiXAADPLMrME4Gr+b0TEL70KXB842+AUoNw5t7WJ60oUuvSkXnROimfmfJ2CKNKYV5aVsKuymmkaiCjqhUwhcM7VAncAbwNrgeedc6vN7FYzuzWw2FwgD8gFHgNuP9q6Qf4IEoIS42O5Zkxv3l1byuad+33HEQkpzjlmLShgcI9kTumb4juOeNbUswx+0Mh6N5nZvxxU6Jz7/fGEcc7N5dAP/YbT/trgsQO+3dR1RQCuG9eHRz/J48lPC/jPC3U2qsiXPt20k3XbKvjNFcN1jI00+SyD7xz2fBtw7WHTHIfOPBAJKWkd23L+iWk8v6SI708ZSPs2Tf1rLxLZZi7Ip0u7BC4e0dN3FAkBTdpl4JzLasLteMYgEGlVUydkUnGwlhdyir56YZEoULCjkvfXbecbY3trnA4Bmndxow5mNukI8yaYWeeWiyXSsk7q3ZmTenfiiYUF1Nd/5RhbIhHviYUFxMUY3zylj+8oEiKac1ChA940swkNJ5rZSOADQBVTQtq0CVkU7NzPh+u3+44i4lVFVQ0vLC3mouE96dYh0XccCRHNubhRBfAKcP1hs74JvO2c29GSwURa2rkn9KBHh0RdBVGi3vM5xew7WMtUnWooDTT3tMPZwJVmFg9gZjEcOrjwiRbOJdLi4mNjuH58Hxbk7mTdtr2+44h4UVfveGJhPidndubE9I6+40gIaW4heBfYD1wUeH4mkAC81pKhRFrLNSf3JjE+hlnzC3xHEfHi/bWlFO06oK0D8m+aVQicc/XA3/i/3QbXAc8552paOphIa+jcLoHLRqXz8rISyioO+o4jEnQz5uXTq1Nbzh7a3XcUCTHHMlLhbOBcM8sALgWebNlIIq1r2oQsqmvrefqzQt9RRIJqedEeFhfsYuqETOJiQ2agWgkRzf4bERgSeCXwDFDsnFvc4qlEWlH/bu05c3A3nvqskKqaOt9xRIJmxvx8ktvEcdXJuhy4/LtjrYhPARM4tLVAJOzcPKkvuyqreenzEt9RRIKiZM8B5q7cyjVje5OcGO87joSgYy0ETwP3ArNaMItI0JzSN4UTenVgxvw8DVQkUWFW4IqfN4zP9BtEQtYxFQLn3C7n3L3OuW0tHUgkGMyMmyf2Ja+sUgMVScTbW1XDs0uKuODENHp1aus7joQoHVUiUeuC4WmkdUxkxjwNVCSR7fklRew7WMstk3TJGTkyFQKJWvGxMdw4PpNP83ayqqTcdxyRVlFbV8+sBQWMzUrRQERyVCoEEtWuHtObdgmxzJiX5zuKSKuYu2obJXsOaOuAfCUVAolqHdvG8/WTM3h9xVa2lh/wHUekRTnnmDEvj75d23HG4G6+40iIUyGQqDdtQhb1zvHEwgLfUURa1JKC3awoLmfaxCxiYsx3HAlxKgQS9TJSkjjvhDSeWbSZfQdrfccRaTGPzcujc1I8l49K9x1FwoAKgQhw86QsKqpqeX5Jke8oIi0if0cl760t5bpT+tA2IdZ3HAkDKgQiwEm9O5PdpzMzF+RTW1fvO47IcXt8fh7xMTF8c1wf31EkTKgQiATcPCmL4t0HeHt1qe8oIsdld2U1Lywt5msn9aRbcqLvOBImVAhEAqYM7UGfLknMmK9TECW8/W1RIVU19dysUw2lGVQIRAJiY4xpE7L4YvMelhbu8h1H5JgcrK3jyU8LOW1gKgO7J/uOI2FEhUCkgSuz0+nYNp7HPtFwxhKeXlm2hbKKg9w8Kct3FAkzKgQiDSQlxPGNsb15e802CndW+o4j0izOOR6fl8/gHslM7N/VdxwJMyoEIoe5YXwmcTHGrAUFvqOINMu8jTtYX1rBzZP6YqaBiKR5VAhEDtO9QyIXjejJ8zlFlO+v8R1HpMkem5dHt+Q2XDyip+8oEoZUCEQacfPEvuyvruNviwt9RxFpknXb9jJv4w5uGJ9JQpz+a5fmC4m/NWaWYmbvmtnGwH3nRpbJMLMPzWytma02szsbzPulmZWY2bLA7fzgfgKJNEN7dmBi/648ubCA6loNVCShb8a8fNrGx/KNsb19R5EwFRKFALgbeN85NwB4P/D8cLXAD51zQ4BTgG+b2dAG8//gnBsZuM1t/cgS6W6elEXp3oO8tnyL7ygiR7V9bxWvLCvhyux0OiUl+I4jYSpUCsElwJOBx08CXzt8AefcVufc54HHFcBaoFewAkr0OW1gKgO6tWfG/Hycc77jiBzR7E8Lqa13TJugUw3l2IVKIejunNsKh37wA0e9cLeZZQInAYsaTL7DzFaY2czGdjk0WHe6meWYWU5ZWVkLRJdIZWbcPCmLtVv3snDTTt9xRBq1v7qWpxcVcvbQ7mR2bec7joSxoBUCM3vPzFY1crukma/THngR+J5zbm9g8sNAP2AksBW4/0jrO+cedc5lO+eyU1NTj+3DSNS4ZGQvurZP4LF5Gs5YQtOLS4vZs79GwxTLcYsL1hs558460jwzKzWzNOfcVjNLA7YfYbl4DpWBvznnXmrw2qUNlnkMeL3lkks0S4yP5fpxmfz+3Q1sLK1ggIaClRBSV+94fH4+IzI6kd3niBtGRZokVHYZvArcEHh8A/DK4QvYoVE2HgfWOud+f9i8tAZPLwVWtVJOiULfPKUPbeJieHy+hjOW0PL+2lIKdu7nlklZGohIjluoFIL7gClmthGYEniOmfU0sy/PGJgAXAec0cjphb8xs5VmtgI4Hfh+kPNLBEtpl8Dlo9N56YsSyioO+o4j8k8z5uXTq1Nbzh3Ww3cUiQBB22VwNM65ncCZjUzfApwfeDwfaLQCO+eua9WAEvVumpjFM4s289RnhfxgykDfcURYXrSHxQW7+PkFQ4iLDZXf7SSc6W+RSBP0S23PWUO68fRnhVTV1PmOI8Jj8/JIbhPHVSdn+I4iEUKFQKSJbp7Ul12V1bz4ebHvKBLlinfv581V27hmbG+SE+N9x5EIoUIg0kRjs1I4oVcHHp+fT329BioSf55YUIABN47P9B1FIogKgUgTmRm3TOpLXlklH65v9MxYkVa3t6qGZ5cUccHwNHp2aus7jkQQFQKRZjj/xDTSOiZqoCLx5rnFRew7WMvNEzUQkbQsFQKRZoiPjWHqhEw+y9vFqpJy33EkytTU1TNrQT5js1I4Mb2j7zgSYVQIRJrp6jG9aZcQywxtJZAge3PVNraUV3GLhimWVqBCINJMHRLjuerk3ry+Yitbyw/4jiNRwjnHjHl59O3ajjMGH/X6byLHRIVA5BhMnZBJvXM8saDAdxSJEovzd7GiuJxpE7OIidEwxdLyVAhEjkFGShLnnZjGM4s3s+9gre84EgUem5dP56R4Lh+V7juKRCgVApFjdMukvlRU1fLckiLfUSTC5ZXt4/11pVx3Sh/aJsT6jiMRSoVA5BiNzOjEyZmdmbUgn9q6et9xJILNXJBPfGwM143L9B1FIpgKgchxuGliX4p3H+Dt1aW+o0iE2l1ZzQtLi7l0ZC9Sk9v4jiMRTIVA5DhMGdqdPl2SeGxeHs5pOGNpeYcuqFXPTZOyfEeRCKdCIHIcYmOMmyZmsaxoD0sLd/uOIxGmqqaOJz8t5LSBqQzsnuw7jkQ4FQKR43TF6HQ6to1nxrx831Ekwry6fAs79h3UQEQSFCoEIscpKSGOb57Sm7fXbGNT2T7fcSRC1NU7Hvskj8E9kpnQv4vvOBIFVAhEWsDUCVm0jY/lj+9t9B1FIsRry7ewcfs+7jijP2YaiEhanwqBSAvo2r4NUydk8tryLazdutd3HAlzNXX1/OG9DQxJ68D5J6T5jiNRQoVApIVMn9SP5MQ47n9ng+8oEuZeWFpM4c793HX2QA1TLEGjQiDSQjomxfOtU/vy3tpSvtisMw7k2FTV1PHA+xs5qXcnXcRIgkqFQKQF3Tghi5R2CdpKIMfsmUWb2VpexY/OHqRjBySoVAhEWlD7NnHcPrkf83N38Ommnb7jSJjZX13LQx/lMr5fF8b37+o7jkQZFQKRFvbNU/rQvUMbfvfOeo1eKM0ya0EBO/ZV88OzB/mOIlFIhUCkhSXGx/KdMwawtHA3H60v8x1HwkT5gRoe+XgTZwzuxug+nX3HkSikQiDSCr6enUFGSlt+98566uu1lUC+2ox5eeytquWHZw/0HUWilAqBSCtIiIvhe2cOZPWWvby1epvvOBLidu47yMz5+VxwYhrDenb0HUeilAqBSCv52km96N+tPb9/dwN12kogR/HwR5s4UFPH96do64D4ExKFwMxSzOxdM9sYuG90B5qZFZjZSjNbZmY5zV1fJJhiY4wfTBlI7vZ9/OOLEt9xJERtK69i9meFXHpSOv27tfcdR6JYSBQC4G7gfefcAOD9wPMjOd05N9I5l32M64sEzbnDejCsZwf++P4GqmvrfceREPTnDzbinON7Zw3wHUWiXKgUgkuAJwOPnwS+FuT1RVpFTIxx19mDKNp1gDlLi3zHkRBTtGs/zy0p4qqTM8hISfIdR6JcqBSC7s65rQCB+yON1+mAd8xsqZlNP4b1MbPpZpZjZjllZTolTFrf5EGpjO7TmT+/n0tVTZ3vOBJC/vjeRmJjjO+coa0D4l/QCoGZvWdmqxq5XdKMl5ngnBsFnAd828xObW4O59yjzrls51x2ampqc1cXaTYz44dnD2Tb3iqe/qzQdxwJEbnbK3j5i2KuH9eH7h0SfccRIS5Yb+ScO+tI88ys1MzSnHNbzSwN2H6E19gSuN9uZi8DY4BPgCatL+LL+H5dmdC/Cw9/tIlrxvSmXZug/dOTEPWHdzfSNj6W2yb39x1FBAidXQavAjcEHt8AvHL4AmbWzsySv3wMnA2saur6Ir7ddfYgdlZWM2tBvu8o4tmqknLeWLmVmyYeuhiWSCgIlUJwHzDFzDYCUwLPMbOeZjY3sEx3YL6ZLQcWA28459462voioeSk3p05a0g3Hvkkj/L9Nb7jiEe/f3cDHdvGc/OpfX1HEfmnkNhu6ZzbCZzZyPQtwPmBx3nAiOasLxJqfjBlEOc/MI9H523iR+cM9h1HPFhauJsP1m3nx+cOokNivO84Iv8UKlsIRKLC0J4duHB4WuCqdgd9xxEPfvf2erq2T+DG8Zm+o4j8CxUCkSD7/pSBVNXU8dCHm3xHkSBbkLuDT/N2cvvk/iQlhMQGWpF/UiEQCbJ+qe25fFQ6Ty8qZGv5Ad9xJEicc/z27fWkdUzk2rG9fccR+TcqBCIefPfMATjneOD9XN9RJEjeX7udZUV7+O6ZA0iMj/UdR+TfqBCIeJCRksQ1Y3ozJ6eIwp2VvuNIK6uvd/zunfX06ZLEFaPTfccRaZQKgYgnd5zen9gY40/vbfQdRVrZ3FVbWbetgu+fNZD4WP23K6FJfzNFPOnWIZEbxmfy8rISNpZW+I4jraS2rp7fv7uBgd3bc9GInr7jiByRCoGIR7ee1o92CXH8/t0NvqNIK3n5ixLyyir5wZRBxMaY7zgiR6RCIOJRSrsEpk3M4s1V21hVUu47jrSw6tp6/vT+Roand+ScYd19xxE5KhUCEc9unpRFx7bx/O6d9b6jSAt7bslmincf4IdnD8JMWwcktKkQiHjWITGeW0/rx0fry8gp2OU7jrSQA9V1/PmDXMZkpnDqgK6+44h8JRUCkRBww/g+dG3fht++vR7nnO840gKe+qyA7RUHuescbR2Q8KBCIBICkhLiuOP0fizK38X83B2+48hxqqiq4eGPNnHqwFTGZKX4jiPSJCoEIiHimrG96dWpLb/TVoKwN3N+Abv313DX2QN9RxFpMhUCkRDRJi6W757Zn+XF5by7ptR3HDlGe/ZXM2NeHucM687w9E6+44g0mQqBSAi5fFQ6WV3b8ft3N1Bfr60E4eivH+exr7qWH0wZ5DuKSLOoEIiEkLjYGL531gDWbavgtRVbfMeRZtpeUcUTC/O5eERPBvVI9h1HpFlUCERCzEXDezKoezJ/fG8jtXX1vuNIMzz04SZq6hzfP0vHDkj4USEQCTExMcYPzh5I/o5KXvq8xHccaaKSPQd4ZtFmrhydTmbXdr7jiDSbCoFICDp7aHdGpHfkT+9v5GBtne840gR/fv/QVSu/c+YAz0lEjo0KgUgIMjN+ePYgSvYc4MmFBb7jyFdYv62COUuLuTZw6qhIOFIhEAlRkwZ05awh3fjdOxtYt22v7zhyBAdr67jz2S/onBTPd87o7zuOyDFTIRAJUWbGfZcPp0NiHN97dhlVNdp1EIp++9Z61m2r4LdXjKBL+za+44gcMxUCkRDWtX0bfnvFiEM/cN7W1RBDzfyNO5gxP5/rTunD6YO7+Y4jclxUCERC3OmDu3H9uD48Pj+fTzaU+Y4jAbsrq/nhnGX0S23HT88f4juOyHFTIRAJAz89fwj9u7XnrjnL2VVZ7TtO1HPO8ZOXVrKrspo/XX0SbRNifUcSOW4qBCJhIDE+lj9dPZLd+6v5yUsrdPEjz+bkFPPW6m3cdfYgTujV0XcckRahQiASJob17MiPzhnE26tLeT6nyHecqFWwo5JfvraacX27cMukvr7jiLSYkCgEZpZiZu+a2cbAfedGlhlkZssa3Paa2fcC835pZiUN5p0f9A8hEgQ3T+zL+H5duPe1NeTvqPQdJ+rU1NXzveeWERdj3P/1EcTEmO9IIi0mJAoBcDfwvnNuAPB+4Pm/cM6td86NdM6NBEYD+4GXGyzyhy/nO+fmBiO0SLDFBH4QxcfG8L3nllGjax0E1Z8/yGVZ0R5+fdmJ9NQARBJhQqUQXAI8GXj8JPC1r1j+TGCTc66wNUOJhKK0jm359aUnsrxoDw8EhsuV1pdTsIu/fLCRy0b14sLhPX3HEWlxoVIIujvntgIE7r/qhN6rgb8fNu0OM1thZjMb2+UgEkkuGJ7G5aPSefDDXJYU7PIdJ+JVVNXwveeW0atzW+69eJjvOCKtImiFwMzeM7NVjdwuaebrJAAXA3MaTH4Y6AeMBLYC9x9l/elmlmNmOWVlOqdbwtcvLx5Keuckvv/cMvZW1fiOE9HueXU1W/Yc4I9XjSQ5Md53HJFWEbRC4Jw7yzl3QiO3V4BSM0sDCNxvP8pLnQd87pwrbfDapc65OudcPfAYMOYoOR51zmU757JTU1Nb5sOJeJCcGM8frhrJ1vIqfvnKat9xItbrK7bw0ucl3HHGAEb3SfEdR6TVhMoug1eBGwKPbwBeOcqy13DY7oIvy0TApcCqFk0nEqJG9+nMHaf356UvSnht+RbfcSLOlj0H+OlLKxmZ0Ynv6sJFEuFCpRDcB0wxs43AlMBzzKynmf3zjAEzSwrMf+mw9X9jZivNbAVwOvD94MQW8e87Z/TnpN6d+NnLKynZc8B3nIhRV+/4wfPLqK13/PGqkcTFhsp/lyKtIyT+hjvndjrnznTODQjc7wpM3+KcO7/Bcvudc12cc+WHrX+dc+5E59xw59zFXx6gKBIN4mJj+ONVIw/9AHtuGXX1GsWwJTw2L4/P8nbxy4uGkdm1ne84Iq0uJAqBiByfPl3acc/Fw1iUv4tHP8nzHSfsrSop5/531nPusB5cmZ3uO45IUKgQiESIK0enc/6JPfj9u+tZVVL+1StIow5U13Hns1+Q0i6B/73sRMw0GqFEBxUCkQhhZvz60hPp0q4N3332Cw5U1/mOFJZ+PXctm8oquf/KkXRul+A7jkjQqBCIRJBOSQnc//UR5JVV8j9z1/iOE3Y+WFfKU58VcvPELCYO6Oo7jkhQqRCIRJgJ/btyy6Qsnv5sM++vLf3qFQSAsoqD/PiFFQzukcyPzh3kO45I0KkQiESgu84ZxJC0Dvz4hRWUVRz0HSfkOef48QvL2VtVywPXnESbuFjfkUSCToVAJAK1iYvlT1ePZN/BWn78wnKc06mIR/P0Z4V8uL6Mn5w3mIHdk33HEfFChUAkQg3snsxPzhvMh+vLeOozXRj0SHK3V/Dfb6zl1IGp3Dg+03ccEW9UCEQi2A3jMzltYCr/88ZaNpZW+I4Tcqpr67nz2WW0axPH764YrlMMJaqpEIhEMDPjt1cOp12bOO58dhkHa3UqYkP3v7ue1Vv28v8uH063Dom+44h4pUIgEuG6JSfym8uHs2brXu5/Z4PvOCFj4aYdPPpJHteM6c2Uod19xxHxToVAJAqcNbQ7147tzWPz8liYu8N3HO/K99fww+eXk9WlHf954RDfcURCggqBSJT4+QVDyOrSjh88v5w9+6t9x/HGOcdP/7GSsoqD/PHqkSQlxPmOJBISVAhEokRSQhx/uvokduw7yPefW0blwVrfkYLOOceMefm8sWIr358ykOHpnXxHEgkZKgQiUeTE9I7cc9FQPtpQxsV/mc+6bXt9RwqavVU13PH3L/ifuWs5a0h3bj2tn+9IIiFFhUAkylw3LpO/3TSW8gO1XPKXBTy3ZHPED1y0qqSci/48n7dWbePH5w7i0etGExujUwxFGlIhEIlC4/t3Ze6dE8nO7Mx/vLiS7z+3jH0RuAvBOceTCwu47KGFVNfW8+z0U7h9cn9iVAZE/o2OphGJUt2SE5k9bSwPfpjLH9/bwIricv5y7SiG9uzgO1qLKD9Qw90vruDNVds4Y3A3fnflCFJ0OWORI9IWApEoFhtjfPfMATxzyynsO1jL1x5awDOLwn8XworiPVz453m8s6aUn54/mBnXZ6sMiHwFFQIR4ZS+XZh75yTGZqXw05dX8t1nl1FRVeM7VrM555g5P5/LH15IXZ3j+W+NY/qp/bSLQKQJtMtARADo2r4NT04dw8Mfb+L+d9azsngPf7l2FCf06ug7WpOU76/hRy8s5501pZw15NAugk5J2iog0lTaQiAi/xQTY3z79P48O30cVTX1XPbQQp76rDDkdyEsK9rD+Q/M44N12/n5BUN47PpslQGRZlIhEJF/MyYrhbl3TmJ8/y785z9WccczX7A3BHchHBpoKI8rHl4IwJxbx3HzpL66aqHIMdAuAxFpVEq7BGbecDKPfJLH795Zz8qSch68dhQnpofGLoQ9+6u5a85y3lu7nbOHdue3V4ygY1K871giYUtbCETkiGJijNsm9+O56adQU1fP5Q8v5MmFBd53ISwt3M0FD8zn4w1l3HPRUB65brTKgMhxUiEQka+UnZnC3O9OYtKArtzz6mpue/pzyg8EfxdCfb3jkY83cdUjnxITAy/cOp6pE7K0i0CkBagQiEiTdG6XwIwbsvnZ+UN4b20pF/55HsuL9gTt/XdVVnPz7Bz+9811TBnande/M4kRGZ2C9v4ikU6FQESazMy45dS+PH/rOOrr4Yq/LuTx+fmtvgthScEuLnhgHvM37uBXlwzjoW+MomNb7SIQaUkqBCLSbKN6d+aN707ktIHd+K/X1zD9qaWU72/5XQj19Y6HPsrl6kc/IyEuhpduH8/14zK1i0CkFZjvg4MAzOxK4JfAEGCMcy7nCMudC/wJiAVmOOfuC0xPAZ4DMoEC4OvOud1f9b7Z2dkuJ6fRtxKRJnDOMXNBAfe9uZaObRPISGnboq+/90ANm8oquWB4GvdddiLJidoqIHI8zGypcy67sXmhctrhKuAy4JEjLWBmscCDwBSgGFhiZq8659YAdwPvO+fuM7O7A8//o/Vji0Q3M+OmiVlk9+nMXz7MpaqmrkVfPzkxnlsm9eWqkzO0VUCklYVEIXDOrQW+6h/8GCDXOZcXWPZZ4BJgTeB+cmC5J4GPUCEQCZoRGZ147PpGf+kQkTARTscQ9AKKGjwvDkwD6O6c2woQuO92pBcxs+lmlmNmOWVlZa0WVkREJJwEbQuBmb0H9Ghk1s+cc6805SUamdbsAyCcc48Cj8KhYwiau76IiEgkClohcM6ddZwvUQxkNHieDmwJPC41szTn3FYzSwO2H+d7iYiIRJVw2mWwBBhgZllmlgBcDbwamPcqcEPg8Q1AU7Y4iIiISEBIFAIzu9TMioFxwBtm9nZgek8zmwvgnKsF7gDeBtYCzzvnVgde4j5giplt5NBZCPcF+zOIiIiEs5AYh8AXjUMgIiLR5GjjEITEFgIRERHxS4VAREREVAhEREREhUBERERQIRARERFUCERERAQVAhERESHKxyEwszKgsAVfsiuwowVfLxxF+3cQ7Z8f9B2AvoNo//wQut9BH+dcamMzoroQtDQzyznSgA/RItq/g2j//KDvAPQdRPvnh/D8DrTLQERERFQIRERERIWgpT3qO0AIiPbvINo/P+g7AH0H0f75IQy/Ax1DICIiItpCICIiIioELcbMzjWz9WaWa2Z3+84TbGY208y2m9kq31l8MLMMM/vQzNaa2Wozu9N3pmAzs0QzW2xmywPfwb2+M/lgZrFm9oWZve47iw9mVmBmK81smZlF3fXlzayTmb1gZusC/x+M852pqbTLoAWYWSywAZgCFANLgGucc2u8BgsiMzsV2AfMds6d4DtPsJlZGpDmnPvczJKBpcDXouzvgAHtnHP7zCwemA/c6Zz7zHO0oDKzHwDZQAfn3IW+8wSbmRUA2c65UDwHv9WZ2ZPAPOfcDDNLAJKcc3s8x2oSbSFoGWOAXOdcnnOuGngWuMRzpqByzn0C7PKdwxfn3Fbn3OeBxxXAWqCX31TB5Q7ZF3gaH7hF1W8cZpYOXADM8J1Fgs/MOgCnAo8DOOeqw6UMgApBS+kFFDV4XkyU/TCQ/2NmmcBJwCLPUYIusLl8GbAdeNc5F23fwR+BHwP1nnP45IB3zGypmU33HSbI+gJlwKzAbqMZZtbOd6imUiFoGdbItKj6zUgOMbP2wIvA95xze33nCTbnXJ1zbiSQDowxs6jZfWRmFwLbnXNLfWfxbIJzbhRwHvDtwO7EaBEHjAIeds6dBFQCYXNMmQpByygGMho8Twe2eMoingT2m78I/M0595LvPD4FNpN+BJzrN0lQTQAuDuxDfxY4w8ye9hsp+JxzWwL324GXObRLNVoUA8UNtoy9wKGCEBZUCFrGEmCAmWUFDiK5GnjVcyYJosABdY8Da51zv/edxwczSzWzToHHbYGzgHVeQwWRc+4nzrl051wmh/4P+MA5903PsYLKzNoFDqolsKn8bCBqzjxyzm0DisxsUGDSmUDYHFgc5ztAJHDO1ZrZHcDbQCww0zm32nOsoDKzvwOTga5mVgzc45x73G+qoJoAXAesDOxDB/ipc26uv0hBlwY8GTjrJgZ43jkXlafeRbHuwMuH+jFxwDPOubf8Rgq67wB/C/xymAdM9ZynyXTaoYiIiGiXgYiIiKgQiIiICCoEIiIiggqBiIiIoEIgIiIiqBCISJCZ2fTAFfEqzex1M+vpO5OIqBCISBCZ2SXAH4AfcGgENwe8amb6v0jEM41DICJBY2ZLgI+dc3cFnqdwaJjvy6JsECeRkKNWLiJBYWadgWzgjS+nOed2AUuBKb5yicghKgQiEixZgfvcw6ZvBPoFOYuIHEbXMhCRYPnyuvBrA2PdfykB+CT4cUSkIRUCEQmWysD9mcDOBtN/DSQGP46INKRCICLBkhe4r3DO/XO3gZm1ATb5iSQiX9IxBCISFM65PcASYNKX0wJl4FTgHU+xRCRAWwhEJJj+C5hpZps5dDDhLwL3b3lNJSIah0BEgsvMbgF+DqQC7wPfcs5t8ZtKRFQIRERERMcQiIiIiAqBiIiIoEIgIiIiqBCIiIgIKgQiIiKCCoGIiIigQiAiIiKoEIiIiAgqBCIiIgL8fz7/Q4XRv25HAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 576x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "from qutip.qip.operations import cnot, hadamard_transform\n",
    "\n",
    "########################################################################\n",
    "\n",
    "def ansatz(theta):\n",
    "    vac = qt.basis(4,0)\n",
    "    vac.dims = [[2,2], [1,1]]\n",
    "    return qt.tensor((-1j*qt.sigmax()*theta/2).expm(), qt.identity(2))*\\\n",
    "           cnot()*\\\n",
    "           qt.tensor(hadamard_transform(), qt.identity(2))*\\\n",
    "           vac\n",
    "\n",
    "########################################################################\n",
    "\n",
    "print(\"our matrix:\")\n",
    "H = qt.Qobj(np.array([[1,0,0,0],\\\n",
    "                      [0,0,-1,0],\\\n",
    "                      [0,-1,0,0],\\\n",
    "                      [0,0,0,1]]))\n",
    "H.dims = [[2]*n, [2]*n]\n",
    "print(H)\n",
    "\n",
    "Hl, Hv = H.eigenstates()\n",
    "print(\"\\ncompare:\")\n",
    "print(\"lowest eigval =  %.4f\" % (Hl[0]))\n",
    "print(\"eigenvector = \\n%s\" % (Hv[0]))\n",
    "print()\n",
    "\n",
    "########################################################################\n",
    "\n",
    "basis = construct_pauli_basis(2)\n",
    "\n",
    "samples = 20\n",
    "stochastic = False\n",
    "verbose = True\n",
    "\n",
    "expectation_values = []\n",
    "thetas = np.linspace(0, 2*np.pi, samples)\n",
    "for theta in thetas:\n",
    "    expectation_values.append(expect(H, ansatz(theta), basis, stochastic=stochastic))\n",
    "    print(\"theta = %.4f, <H> = %.4f\" % (theta, expectation_values[-1])) if verbose else None\n",
    "least_index = np.argmin(expectation_values)\n",
    "\n",
    "final_theta = thetas[least_index]\n",
    "final_vector = ansatz(final_theta)\n",
    "print(\"\\nfinal theta = %.4f, <H> = %.4f\" % (final_theta, expectation_values[least_index]))\n",
    "print(\"final vector =\\n%s\" % ansatz(final_theta))\n",
    "\n",
    "########################################################################\n",
    "\n",
    "fig = plt.figure(figsize=(8,6))\n",
    "ax = fig.add_subplot(111)\n",
    "ax.plot(thetas, expectation_values)\n",
    "ax.set_xlabel('θ', fontsize=14)\n",
    "ax.set_ylabel('<H>', fontsize=14)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr>\n",
    "\n",
    "That's all well and good, but let's actually run this baby on a quantum computer.\n",
    "\n",
    "We use IBM's framekwork qiskit. The algorithm directly translates over, although we construct all the circuits ahead of time, transpile/simplify them, and then send them in bulk to the quantum computer (or classical simulator), divvied up into several \"jobs.\"\n",
    "\n",
    "Make sure to have an <a href=\"https://quantum-computing.ibm.com/\">account</a> with IBM, and replace ```YOUR_TOKEN``` to run things non-locally."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "our matrix:\n",
      "Quantum object: dims = [[2, 2], [2, 2]], shape = (4, 4), type = oper, isherm = True\n",
      "Qobj data =\n",
      "[[ 1.  0.  0.  0.]\n",
      " [ 0.  0. -1.  0.]\n",
      " [ 0. -1.  0.  0.]\n",
      " [ 0.  0.  0.  1.]]\n",
      "\n",
      "compare:\n",
      "lowest eigval =  -1.0000\n",
      "eigenvector = \n",
      "Quantum object: dims = [[2, 2], [1, 1]], shape = (4, 1), type = ket\n",
      "Qobj data =\n",
      "[[0.        ]\n",
      " [0.70710678]\n",
      " [0.70710678]\n",
      " [0.        ]]\n",
      "\n",
      "laying out circuits...\n",
      "evaluating...\n",
      "\n",
      "theta = 0.0000, <H> = 1.00000000\n",
      "theta = 0.1065, <H> = 0.99340820\n",
      "theta = 0.2130, <H> = 0.97875977\n",
      "theta = 0.3195, <H> = 0.95410156\n",
      "theta = 0.4260, <H> = 0.90258789\n",
      "theta = 0.5325, <H> = 0.87036133\n",
      "theta = 0.6390, <H> = 0.80346680\n",
      "theta = 0.7455, <H> = 0.72900391\n",
      "theta = 0.8520, <H> = 0.64843750\n",
      "theta = 0.9585, <H> = 0.57299805\n",
      "theta = 1.0649, <H> = 0.47924805\n",
      "theta = 1.1714, <H> = 0.38647461\n",
      "theta = 1.2779, <H> = 0.29516602\n",
      "theta = 1.3844, <H> = 0.20361328\n",
      "theta = 1.4909, <H> = 0.08154297\n",
      "theta = 1.5974, <H> = -0.03222656\n",
      "theta = 1.7039, <H> = -0.12109375\n",
      "theta = 1.8104, <H> = -0.24145508\n",
      "theta = 1.9169, <H> = -0.36254883\n",
      "theta = 2.0234, <H> = -0.43310547\n",
      "theta = 2.1299, <H> = -0.52905273\n",
      "theta = 2.2364, <H> = -0.63378906\n",
      "theta = 2.3429, <H> = -0.69677734\n",
      "theta = 2.4494, <H> = -0.77612305\n",
      "theta = 2.5559, <H> = -0.84375000\n",
      "theta = 2.6624, <H> = -0.88989258\n",
      "theta = 2.7689, <H> = -0.92382812\n",
      "theta = 2.8754, <H> = -0.96533203\n",
      "theta = 2.9819, <H> = -0.98632812\n",
      "theta = 3.0883, <H> = -0.99902344\n",
      "theta = 3.1948, <H> = -0.99853516\n",
      "theta = 3.3013, <H> = -0.98901367\n",
      "theta = 3.4078, <H> = -0.96508789\n",
      "theta = 3.5143, <H> = -0.92456055\n",
      "theta = 3.6208, <H> = -0.88183594\n",
      "theta = 3.7273, <H> = -0.83032227\n",
      "theta = 3.8338, <H> = -0.77270508\n",
      "theta = 3.9403, <H> = -0.69238281\n",
      "theta = 4.0468, <H> = -0.61157227\n",
      "theta = 4.1533, <H> = -0.53295898\n",
      "theta = 4.2598, <H> = -0.44287109\n",
      "theta = 4.3663, <H> = -0.33837891\n",
      "theta = 4.4728, <H> = -0.26196289\n",
      "theta = 4.5793, <H> = -0.13403320\n",
      "theta = 4.6858, <H> = -0.03930664\n",
      "theta = 4.7923, <H> = 0.08374023\n",
      "theta = 4.8988, <H> = 0.19897461\n",
      "theta = 5.0052, <H> = 0.28466797\n",
      "theta = 5.1117, <H> = 0.40942383\n",
      "theta = 5.2182, <H> = 0.48901367\n",
      "theta = 5.3247, <H> = 0.59033203\n",
      "theta = 5.4312, <H> = 0.65795898\n",
      "theta = 5.5377, <H> = 0.73120117\n",
      "theta = 5.6442, <H> = 0.80786133\n",
      "theta = 5.7507, <H> = 0.85571289\n",
      "theta = 5.8572, <H> = 0.90991211\n",
      "theta = 5.9637, <H> = 0.94653320\n",
      "theta = 6.0702, <H> = 0.97900391\n",
      "theta = 6.1767, <H> = 0.99584961\n",
      "theta = 6.2832, <H> = 1.00000000\n",
      "\n",
      "final_theta: 3.0883, <H> = -0.9990\n",
      "final vector = \n",
      "Quantum object: dims = [[4], [1]], shape = (4, 1), type = ket\n",
      "Qobj data =\n",
      "[[ 1.88235512e-02+0.j        ]\n",
      " [ 4.32824586e-17-0.70685619j]\n",
      " [-4.32824586e-17-0.70685619j]\n",
      " [ 1.88235512e-02+0.j        ]]\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, '<H>')"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgQAAAF3CAYAAADaXFNkAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAABBOElEQVR4nO3dd3yV5d3H8c8vE8IKIQFCEnYYYcgIoIiKIgouXDhq1VZba9U+tnbZ5agdPh1WrdYWV7XuLU5QRBFFIOwNYSeMhJUAIWRdzx85+qQYIIHkXOfkfN+vV17nnHuE70mAfHPd933d5pxDREREIluU7wAiIiLinwqBiIiIqBCIiIiICoGIiIigQiAiIiKoEIiIiAgQ4zuAT8nJya5r166+Y4iIiATFvHnzdjjnUmpbF9GFoGvXruTk5PiOISIiEhRmtvFw63TIQERERFQIRERERIVAREREUCEQERERVAhEREQEFQIRERFBhUBERERQIRARERFUCERERIQQKwRm9oSZFZjZ0sOsNzN70MxyzWyxmQ2psW6cma0KrLs9eKlFRETCX0gVAuDfwLgjrB8PZAY+bgAeATCzaODhwPos4Eozy2rUpCIiIk1ISN3LwDk3w8y6HmGTCcDTzjkHfGFmiWaWCnQFcp1z6wDM7IXAtssbOfJX1hbuY3tRKV2SW5DauhlRURasP1pERJoI5xzFpRXs2HeQHXsPAjCie7ug/NkhVQjqIA3YXON1XmBZbctH1PYJzOwGqkcX6Ny5c4MFe3VeHv/4eC0AcTFRdElKoEu7FnRLrn7M6tSaQemJKgoiIgJA/p4DvLEgn/kbd1MYKAA79pdRVlH11TZ9U1vz3q2nBCVPuBWC2n6auiMs//pC5yYBkwCys7Nr3eZYfPvkbozqmcyGnSVs2LmfDTv2s3FnCZ+uKeRg4Jub0iqesVkdOLtfR07q3o64mFA7YiMiIo1pb2k57y3ZxmsL8vhi3S4AendoRcc2zchs34rkVnGktIwnuWU87VrGkdqmWdCyhVshyAMyarxOB7YAcYdZHjQpreJJaRXPyJ7/vbyqyrGtuJS5G3Yxddl23liQz3OzN9GqWQxn9GnP2f06Mrp3Cglx4fatEBGRuqiscnyyuoDX5ufzwfLtHKyooltyC24b24sLB6XRuV2C74hA+BWCycAtgXMERgBFzrmtZlYIZJpZNyAfuAL4hsecX4mKMjolNmfCoDQmDEqjtLySz3J3MGXZNj5cUcCbC7eQ1CKOm0/vyTdP7Ex8TLTvyCIi0kA27SzhRy8tZN7G3bRNiOXyYRlcNDiNQRmJmIXWIeSQKgRm9jwwGkg2szzgTiAWwDn3T+Bd4BwgFygBvh1YV2FmtwBTgGjgCefcsqC/gTpoFhvNmL4dGNO3AxWVVczZsIt/TF/LPW8v54mZ6/nR2F5cNDiNaJ1rICIStpxzvDwvj7snLyMqyvjTpQO5cFBaSB8qtuoT9iNTdna2y8nJ8R0DgJlrdvCnKStZnFdEZvuW/PTs3ozN6hByDVJERI5s1/4yfvHaYqYs286J3ZP462WDSEts7jsWAGY2zzmXXes6FYLQKARQ3SjfW7qNv0xZxbod+xncOZF7JvSnf1ob39FERKQOpq8q4GevLGZPSRk/Pbs33xnVPaSuLlMhOIxQKwRfqqis4pV5edz3wWoOVlTx6vdPomf7Vr5jiYjIYZSWV/KHd1fw9KyN9OrQkvsvH0xWp9a+Y33NkQpB6B7MiGAx0VFcMbwzr9w4ktjoKK55fA5biw74jiUiIrXYvKuESx75nKdnbeT6Ud2YfMuokCwDR6NCEMI6t0vg398eRnFpBd96Yi5FJeW+I4mISA2frinkgodmsmlXCY9fm81vzsuiWWx4Xi2mQhDi+qe1YdLVQ1m3Yx/ffTqH0vJK35FERCKec45HPl7LtU/MoX2rZrx1yyjG9O3gO9ZxUSEIAyN7JnPfZYOYs2EXt76wgMqqyD3vQ0TEt30HK7jp2fn87/srOWdAKq/dNJKuyS18xzpuKgRh4vwTOnHHeVlMWbadO95cSiSfDCoi4svawn1c+PBnTF2+nV+f25e/XzmYFvEhNaXPMWsa7yJCXDeqG9v3lvKvT9bRoXUz/mdMpu9IIiIRY3HeHr7x6GziY6L4z/XDGdkj2XekBqVCEGZuH9eHwr0Hue+D1VQ5xy2n9yQmWgM9IiKNKX/PAa5/KofEhFhe/N5JITPRUENSIQgzZsb/XjKQqirH/R+uYcbqQu67bFCTOH4lIhKK9paWc92Tcyktr+S574xokmUAdA5BWIqNjuL+KwbzwBWDyC3YxzkPfsrzczbpvAIRkQZWUVnFLc8tYG3hPh65aiiZHZruJHEqBGFswqA03v/hqQzunMgvXlvCd57KoXDvQd+xRESaBOccd721jE9WF/K7C/szKrNpnTNwKBWCMNcpsTn/uW4Ed5yXxczcHZx9/wymLNvmO5aISNh7fOZ6nvliE987rTtXDO/sO06jUyFoAqKijOtGdePtH4witU0zvvefeTw7e6PvWCIiYWvqsm38/t0VjO/fkZ+f3cd3nKBQIWhCMju04vWbTuaUzGTueXs56wr3+Y4kIhJ2luQVcesLCxmYnsh9lw0KqbsVNiYVgiYmLiaKP196AvEx0fzopUVUVFb5jiQiEjY27tzP9U/NJalFHI9eM5TmceF5X4JjoULQBHVs04zfXdifRZv38PD0tb7jiIiEhbzdJXzj0dmUV1bx5LeH0b5VM9+RgkqFoIk6/4ROTBjUiQc/WsPivD2+44iIhLRtRaVc9dhs9paW85/rR9CrCV9eeDgqBE3Yby/oT0rLeH704kIOlOkuiSIitSnce5BvPPYFO/eV8dR1w+mf1sZ3JC9UCJqwNgmx/GXiCawt3M//vr/SdxwRkZCza38Z33xsNlv3lPLkt4cxuHNb35G8USFo4kZlJvOtkV359+cb+HRNoe84IiIho6iknG8+NpsNO/fz+LXZDOua5DuSVyoEEeD28X3okdKCn768mKKSct9xRES821tazjVPzCa3YB//unooI3s27VkI60KFIAI0i43m/ssHs2PfQX7z5lLfcUREvKqsclz/VA7LthTz8FVDGN27ve9IIUGFIEIMSG/DrWMymbxoC+8v3eo7joiIN8/N3sic9bv448UDGJvVwXeckKFCEEG+P7oHWamtuXPyMvaW6tCBiESeHfsO8ucpqxjZox2XDk33HSekqBBEkJjoKP5w8QAK9h7kr1NX+44jIhJ09763kgPllfx2Qn/MImNK4rpSIYgwgzIS+eaILjw9a4MmLBKRiJKzYRevzMvjO6d0p2f7lr7jhBwVggj003G9adcynl++vkT3OhCRiFBRWcWv31hKWmJzfnBGT99xQlJIFQIzG2dmq8ws18xur2X9T81sYeBjqZlVmllSYN0GM1sSWJcT/PTho3WzWO48P4ul+cU8PUu3SRaRpu+pWRtZuW0vvzkvi4S4GN9xQlLIFAIziwYeBsYDWcCVZpZVcxvn3J+dc4Occ4OAXwCfOOd21djk9MD67GDlDlfnDkhldO8U/jp1FVuLDviOIyLSaLYXl/K3D1ZzWq8Uzu6nqwoOJ2QKATAcyHXOrXPOlQEvABOOsP2VwPNBSdYEmRn3TOhPpXPcNXmZ7zgiIo3m9++soKyyirsv6KcTCY8glApBGrC5xuu8wLKvMbMEYBzwao3FDphqZvPM7IZGS9mEZCQlcOuYXkxZtp0Plm/3HUdEpMF9nruDyYu2cONpPeia3MJ3nJAWSoWgttrmDrPt+cBnhxwuONk5N4TqQw43m9mptf4hZjeYWY6Z5RQWam7/75zSjd4dWnHnm0vZf7DCdxwRkQZTVlHFb95cSkZSc24a3cN3nJAXSoUgD8io8Tod2HKYba/gkMMFzrktgccC4HWqD0F8jXNuknMu2zmXnZKSctyhw11sdBR/uLg/W4pKuf9DzU0gIk3H4zPXs7ZwP3df0I9msdG+44S8UCoEc4FMM+tmZnFU/9CffOhGZtYGOA14s8ayFmbW6svnwFmAJu2vo6FdkrhyeGee+GwDK7cV+44jInLcNu8q4cFpazizbwfO6KMTCesiZAqBc64CuAWYAqwAXnLOLTOzG83sxhqbXgRMdc7tr7GsAzDTzBYBc4B3nHPvByt7U/Dzcb1p1SyGe95ejnOHO1IjIhL6nHP8+o2lRBn8dkI/33HCRkhdjOmcexd495Bl/zzk9b+Bfx+ybB1wQiPHa9ISE+L40Zm9uHPyMj5Yvp2z+nX0HUlE5JhMXrSFT1YXcuf5WXRKbO47TtgImREC8e+qEZ3JbN+S37+7goMVlb7jiIjU256SMn771nJOSG/DNSd19R0nrKgQyFdioqO44/wsNu4s4cnPNviOIyJSb79/ZwV7DpTzx4sHEh2lOQfqQ4VA/sspmSmc2bc9D32US8HeUt9xRETq7PPcHbw8L48bTu1OVqfWvuOEHRUC+ZpfnZvFwYpK/jJlle8oIiJ1UlpeyS9fX0KXdgncOibTd5ywpEIgX9MtuQXfPrkbL8/LY0leke84IiJH9dBHuWzYWcLvLxygOQeOkQqB1OqWM3qSlBDH3W8t02WIIhLSVm3byz8/WcvFQ9IYlZnsO07YUiGQWrVuFstPzu5NzsbdvL14q+84IiK1qqpy3P7aYlo1i+HX52YdfQc5LBUCOazLsjPom9qae99byYEyXYYoIqHnmdkbWbBpD785L4ukFnG+44Q1FQI5rOgo487zs8jfc4BJM9b5jiMi8l827tzPve+t5JTMZC4aXOvNcaUeVAjkiE7s3o5zBnTkkU9yydtd4juOiAgAlVWOH7+0iOgo438vGYiZ5hw4XioEclS/OjcLw7hr8nLfUUREAHjs03XkbNzN3Rf00/TEDUSFQI4qLbE5Pzwzkw9XbGfqsm2+44hIhFu5rZi/Tl3N2f066FBBA1IhkDq5blQ3endoxd1vLaekrMJ3HBGJUGUVVdz24iJaN4/hDxcN0KGCBqRCIHUSGx3F7y7qT/6eAzwwbY3vOCISoR6ctoblW4v5w0UDaNcy3necJkWFQOpsWNckLstO5/FP17Nq217fcUQkwszftJt/fJzLpUPTdYv2RqBCIPVy+/i+1ROAvLGEqirNYCgiwXGgrJIfv7SI1DbNueN8TUDUGFQIpF6SWsTxi/F9mbthN6/Mz/MdR0QixL3vrWD9jv38eeJAWjeL9R2nSVIhkHq7dGg62V3a8sd3V7B7f5nvOCLSxH2Wu4OnZm3k2yd3ZWQP3augsagQSL1FRRm/u6g/e0sruPe9lb7jiEgT5pzj3vdW0jkpgZ+P6+M7TpOmQiDHpE/H1lw/qhsv5mwmZ8Mu33FEpImatW4nS/KLuPG0HrqtcSNTIZBjduuZmaQlNufXbyzVCYYi0igmzVhHcss4Lh6iCYgamwqBHLOEuBh+Nq43K7ftZYpmMBSRBrZyWzEfryrkWyO7anQgCFQI5LicN7AT3VNa8MC0NRolEJEGNWnGOhLiovnmiV18R4kIKgRyXKKjjB+c0ZOV2/Yydfl233FEpInYWnSAyQu3cFl2BokJcb7jRAQVAjlu5w/sRLfkFjw4bQ3OaZRARI7fk59twAHXj+rmO0rEUCGQ4xYTHcXNp/dk+dZiPlxR4DuOiIS54tJynpu9iXMGpJKRlOA7TsRQIZAGceGgTnRpl8AD01ZrlEBEjsvzszex72AF3zu1u+8oEUWFQBrEl6MES/OLmb5KowQicmzKKqp44rP1nNyzHf3T2viOE1FCqhCY2TgzW2VmuWZ2ey3rR5tZkZktDHzcUdd9pfFdNDiNjKTmPPChziUQkWPz5sJ8thcf5IZTe/iOEnFCphCYWTTwMDAeyAKuNLPabmn1qXNuUODjt/XcVxpRbHQUN4/uyaK8Ij5ZXeg7joiEGeccj366jj4dW3Fqpu5ZEGwhUwiA4UCuc26dc64MeAGYEIR9pQFdPCSdtMTmPKArDkSknj5eVcjq7fu44dTumJnvOBEnlApBGrC5xuu8wLJDnWRmi8zsPTPrV899pZHFxURx0+k9WLBpD5+u2eE7joiEkX/NWEtqm2acf0In31EiUigVgtrq4KG/Ys4HujjnTgD+DrxRj32rNzS7wcxyzCynsFDD2o3h0qHpdGrTTKMEIlJni/P28MW6XVw/qhux0aH0oylyhNJXPQ/IqPE6HdhScwPnXLFzbl/g+btArJkl12XfGp9jknMu2zmXnZKS0pD5JSA+Jprvj+7BvI27+XztTt9xRCTEOef43/dX0qZ5LFcM7+w7TsQKpUIwF8g0s25mFgdcAUyuuYGZdbTAgSUzG051/p112VeC67JhGXRs3UxXHIjIUU1Ztp3Pcndy29hetIyP8R0nYoVMIXDOVQC3AFOAFcBLzrllZnajmd0Y2OxSYKmZLQIeBK5w1WrdN/jvQr705SjBnA27mLVOowQiUrvS8kp+985yendoxVUjNDrgk0Xyb2/Z2dkuJyfHd4wmq7S8ktP+PJ0u7Vrw4g0n6qxhEfmav09bw18/WM1z3xnByJ661LCxmdk851x2betCZoRAmp5msdHcNLonc9ZrlEBEvm5r0QH+8fFaxvXrqDIQAlQIpFFdPiyDDq3juV/nEojIIf747kqqnONX5/b1HUVQIZBGplECEanN3A27mLxoC987tbvuaBgiVAik0WmUQERqqqxy3DV5GaltmnHjaN2zIFSoEEij0yiBiNT04tzNLNtSzC/O6UtCnC4zDBUqBBIUGiUQEYCiknL+MnUVw7smcf7AVN9xpAYVAgkKjRKICMD901azp6SMOy/I0qXIIUaFQIJGowQikS23YB9Pz9rIFcM7069TG99x5BAqBBI0GiUQiWzPfLGRaDN+PLaX7yhSCxUCCSqNEohEprKKKt5cmM/YrA60axnvO47UQoVAgkqjBCKRafqqAnaXlHPJ0DTfUeQwVAgk6L4cJfjbB6s1SiASIV6dl0dyy3hOzdRt50OVCoEEXbPYaG4+vSdzN+zms1yNEog0dbv2lzF9VQEXDe5ETLR+7IQqfWfEi8uHZZDaphn3fbBKowQiTdzkhfmUVzouGZruO4ocgQqBeBEfE80tZ/Rk/qY9fLK60HccEWlEr87Pp1+n1vTp2Np3FDkCFQLxZuLQDNISm+tcApEmbNW2vSzJL+KSIRodCHUqBOJNXEwU/zOmJ4vyivhoZYHvOCLSCF6dn0dMlDFhUCffUeQoVAjEq4uHpNM5KYH7NEog0uRUVFbx+oJ8Rvdur7kHwoAKgXgVGx3F/4zJZNmWYqYu3+47jog0oE9zd1C49yCX6mTCsKBCIN5dOKgT3ZJb8LcPVlNVpVECkabi1Xl5tE2I5Yw+7X1HkTpQIRDvYqKjuHVMJiu37eX9Zdt8xxGRBlB0oJypy7dzwQmdiIvRj5pwoO+ShITzT+hEj5TqUYJKjRKIhL13Fm+lrKJKcw+EERUCCQnRUcYPz+zFmoJ9vLNkq+84InKcXp2fR2b7lgxI022Ow4UKgYSMcwek0qtDS+7/UKMEIuFsXeE+5m3czSVD0zEz33GkjlQIJGRERRk/OrMX6wr38+bCfN9xROQYvTY/nyiDiwbrzobhRIVAQsrZ/TrSN7U1D05bQ0Vlle84IlJPVVWO1xfkc0pmCh1aN/MdR+pBhUBCSlSU8cMzM9mws4TXF2iUQCTcTFtZQP6eAzqZMAypEEjIOSurA/06tebvH+VSrlECkbBRWeX485SVdEtuwfj+HX3HkXoKqUJgZuPMbJWZ5ZrZ7bWsv8rMFgc+PjezE2qs22BmS8xsoZnlBDe5NCQz47axvdi0q4TX5uf5jiMidfTqvDxWb9/HT8/uTWx0SP14kToIme+YmUUDDwPjgSzgSjPLOmSz9cBpzrmBwD3ApEPWn+6cG+Scy270wNKozujTnhPS2/DgtFzKKjRKIBLqSssrue+D1QzKSNToQJgKmUIADAdynXPrnHNlwAvAhJobOOc+d87tDrz8AtBBqibKzPjh2F7k7znAK/M0SiAS6p78bAPbiku5fXwfXWoYpkKpEKQBm2u8zgssO5zrgfdqvHbAVDObZ2Y3NEI+CbLRvVIY3DmRhz5aw8GKSt9xROQw9pSU8Y+PczmjT3tO7N7Odxw5RqFUCGqrlLXOTmNmp1NdCH5eY/HJzrkhVB9yuNnMTj3MvjeYWY6Z5RQWFh5vZmlEX55LsKWolJfmbj76DiLixcPTc9l3sIKfj+vjO4och1AqBHlARo3X6cCWQzcys4HAY8AE59zOL5c757YEHguA16k+BPE1zrlJzrls51x2SkpKA8aXxjCqZzLDurbl4elrKS3XKIFIqMnbXcJTn2/kkiHp9O7YynccOQ6hVAjmAplm1s3M4oArgMk1NzCzzsBrwNXOudU1lrcws1ZfPgfOApYGLbk0GrPq2Qu3FZfywpxNvuOIyCHu+2A1GNw2tpfvKHKcQqYQOOcqgFuAKcAK4CXn3DIzu9HMbgxsdgfQDvjHIZcXdgBmmtkiYA7wjnPu/SC/BWkkJ/Vox4huSTz8sUYJRELJ8i3FvL4gn2+P7EqnxOa+48hxMuci9yYy2dnZLidHUxaEgy/W7eSKSV/wm/OyuH5UN99xRAT41pNzWLBpDzN+ejptEmJ9x5E6MLN5h7s0P2RGCESO5MTu7RjZox2PfJxLSVmF7zgiEe/ztTv4eFUhN5/eQ2WgiVAhkLBx29he7NhXxss5mpdAxCfnHPe+t5JObZpxzUldfceRBqJCIGEju2sS/dNa81KOLkEU8WnW2p0szivi1jMzaRYb7TuONBAVAgkrE4dmsGxLMcu2FPmOIhKxnp29icSEWCYMOtLccRJuVAgkrEwY1Im46CgdNhDxpHDvQaYs28YlQ9I1OtDEqBBIWElMiGNsvw68uTBfNz0S8eDleZupqHJcObyz7yjSwFQIJOxMHJrO7pJypq3Y7juKSESpqnK8MGczI7ol0bN9S99xpIGpEEjYOSUzhY6tm/Gy7oIoElQzc3ewaVcJ3xih0YGmSIVAwk50lHHxkDQ+XlXA9uJS33FEIsZzszeR1CKOcf07+o4ijUCFQMLSxOwMqhy8Nj/fdxSRiFBQXMoHK7Zz6dB04mN0MmFTpEIgYalbcguGdW3Ly/M2E8nTb4sEy0s5m6nUyYRNmgqBhK2JQzNYV7if+Zv2+I4i0qRVVjmen7OZkT3a0S25he840khUCCRsnTMwleax0bysmQtFGtWMNYXk7zmgkwmbOBUCCVst42M4d2Aqby/eqhseiTSi52ZvIrllHGdl6WTCpkyFQMLaxKHp7DtYwftLt/mOItIkbSsq5aOVBVw6NIO4GP3IaMr03ZWwNrxbEl3aJWgqY5FG8uLcL08mzPAdRRqZCoGENTPj0iHpzFq3k007S3zHEWlSKqscL87dxCmZyXRpp5MJmzoVAgl7lwxNxwxema9RApGG9PGqArYUlfINXWoYEVQIJOx1SmzOqJ7JvDovj6oqzUkg0lCem72JlFbxnJnVwXcUCYJjKgRmlmRmd5iZTjmVkDAxO4P8PQf4bO0O31FEmoRtRaVMX1XAxKHpxEbrd8dIcKzf5auAO4FvNVwUkWN3VlYH2ibE8tzsTb6jiDQJry3Io8rBZdk6mTBSHGshuAaYHXgU8a5ZbDSXZWcwdfl23fBI5Dg553hlXh7Duralq2YmjBj1LgRmlgUMBL4BZJjZsAZPJXIMrhzeOXBWtGYuFDke8zftYV3hfiYO1ehAJDmWEYJrgfedcxuAN9BhAwkRXZNbcEpmMs/P2URFZZXvOCJh65V5eTSPjeacgam+o0gQ1asQmFkU1ecP/Cew6BngcjOLbehgIsfiqhFd2FpUyvRVhb6jiISlA2WVvL1oC+MHdKRlfIzvOBJE9R0hOBNIACYHXn8AlAHnNWQokWM1pm97OrSO59nZG31HEQlLU5dvY+/BCi4dmu47igRZfQvBNcDLzrkyAOdcFfAcOmwgISI2OorLh3Xmk9WFbN6lmQtF6uvlnDzS2zbnxG7tfEeRIKtzITCzVsBF/P/hgi89C4w3s+SGDCZyrK4cnoEBz83RJYgi9fHlXB6XDEknKsp8x5Egq88IQRQw3jk3s+ZC59wC4Ayg8njDmNk4M1tlZrlmdnst683MHgysX2xmQ+q6r0SO1DbNGdO3Ay/N3UxZhU4uFKmr1+bl4Rw6XBCh6lwInHNFzrkZh1k30zm3+3iCmFk08DAwHsgCrgxc4ljTeCAz8HED8Eg99pUIctWIzuzcX8aUZbotskhdOOd4ZX4eJ3ZPIiMpwXcc8SCU5qMcDuQ659YFzlF4AZhwyDYTgKddtS+ARDNLreO+EkFOzUwhI6k5z3yhkwtF6mLuht1s3FmiuQciWJ2uKTGz9cBR7xrjnOt+HFnSgJozyuQBI+qwTVod95UIEhVlfGN4F/73/ZXkFuylZ/tWviOJhLRX5m2mRVw04wfoFjWRqq4XmT5U47kBvwPuA3Y2YJbazmA5tIQcbpu67Fv9CcxuoPpwA50765aeTdnE7HTu+2AVz3yxibsu6Oc7jkjIKimr4J3FWzl3YCoJcZp7IFLV6TvvnPtrzddmdifwmHNuXQNmyQNqjlWlA1vquE1cHfYFwDk3CZgEkJ2drXvlNmHJLeMZ1z+VV+fn8fNxfWgeF+07kkhIem/JNvaXVXKpDhdEtFA6h2AukGlm3cwsDriC/58A6UuTgWsCVxucCBQ557bWcV+JQN8c0Zm9pRW8tbjWfigiVE9V3KVdAsO6tvUdRTwKmULgnKsAbgGmACuAl5xzy8zsRjO7MbDZu8A6IBd4FLjpSPsG+S1ICBreLYme7VvyrE4uFKnV5l0lzFq3k0uHpGOmuQciWUgdLHLOvUv1D/2ay/5Z47kDbq7rviJmxlUjOnP3W8vJ2bCL7K5JviOJhJRX5+dhBhdr7oGIV9erDG6rZb/rzey/Tip0zt3XUMFEGsrlwzJ4eHouf5m6iue/e6J+CxIJKK+s4uWcPE7ukUxaYnPfccSzuo4Q/OCQ19uAbxyyzFF95YFISEmIi+GW03ty11vL+Sx3J6MyNcu2CMDr8/PJ33OA313Y33cUCQF1vcqgW2MHEWlMV47ozKOfrufPU1Zycs+TNUogEa+isoqHpucyIK0No3un+I4jIaA+NzdqbWanHGbdyWam01MlZMXHRHPrmEwW5RXxwfLtvuOIeDd50RY27SrhljN6qiALUL+rDBzwnpmdXHOhmQ0CPgJ0kbeEtIuHpNE9uQV/nbqayipNQSGRq7LK8dD0XPp0bMXYvh18x5EQUZ+bG+0F3gSuOWTVN4EpzrkdDRlMpKHFREfxo7G9WLV9L29rXgKJYO8u2cq6wv384IxM3eZYvlLfeQieBiaaWSyAmUVRfXLhvxs4l0ijOHdAKn1TW3PfB6spr9StkSXyVFU5/v7RGnq2b8n4/rpvgfy/+haCD4AS4PzA6zFUTxv8VkOGEmksUVHGT87qxcadJbwyL893HJGgm7p8G6u37+MHZ/TU6ID8l3oVAudcFfAs/3/Y4GrgRedceUMHE2ksZ/Rpz5DOiTzw4RpKyyt9xxEJGuccD07LpVtyC84b2Ml3HAkxxzJ18dPAODPLAC4CnmrYSCKNy8z4ydm92VZcyjOa0lgiyLQVBSzfWsxNo3sQrdEBOUS9C0HgHgFLgOeAPOfcnAZPJdLIRvZIZlTPZB75eC37Dlb4jiPS6JyrPncgI6k5Fw5O8x1HQtCx3tzoP8DJVI8WiISln5zdm537y3hy5nrfUUQa3Yw1O1iUV8RNo3sSGx0y97WTEHKsfyueAe4GnmzALCJBNSgjkbFZHZj06Tr2a5RAmrDqcwfW0KlNMy4ZopsYSe2OqRA453Y55+52zm1r6EAiwXTjad3ZW1rBO4u3+o4i0mhmrd3JvI27+f7oHsTFaHRAaqe/GRLRhnRuS4+UFryYs9l3FJFG88gna2nfKp6J2Rm+o0gIUyGQiGZmXD4sg3kbd5NbsM93HJEGt62olJm5O7hieGeaxWqGeTk8FQKJeBcNTicmynhZowTSBL25MB/n4CJdWSBHoUIgES+lVTxn9GnPq/PzNJ2xNDmvL8hnUEYi3ZJb+I4iIU6FQAS4fFgGO/aVMX1lge8oIg1m+ZZiVm7by8VDNDogR6dCIAKc1iuF9q3ieUmHDaQJeWNhPjFRpmmKpU5UCESovjXyJUPTmb6qkILiUt9xRI5bZZXjzYX5jO7dnqQWcb7jSBhQIRAJmDg0ncoqx6vz831HETlun6/dwfbigzpcIHWmQiAS0D2lJcO7JvFyzmacc77jiByX1+fn06pZDGf0ae87ioQJFQKRGi4blsG6HfvJ2bjbdxSRY1ZSVsH7y7Zx7oBUzT0gdaZCIFLDOQM60jI+hhfn6uRCCV9Tl22npKxScw9IvagQiNSQEBfD+Sek8s7irewtLfcdR+SYvLYgn7TE5gzrmuQ7ioQRFQKRQ1yWncGB8krd8EjCUkFxKTPXFHLR4DSiosx3HAkjKgQihxiUkUhm+5a64ZGEpcmLtlDl4CJdXSD1FBKFwMySzOwDM1sTeGxbyzYZZjbdzFaY2TIzu7XGurvMLN/MFgY+zgnuO5Cm5MsbHi3YtIc12/f6jiNSL68vyOeE9Db0SGnpO4qEmZAoBMDtwDTnXCYwLfD6UBXAj51zfYETgZvNLKvG+r855wYFPt5t/MjSlF00OI2YKNPMhRJWVm/fy7ItxVyokwnlGIRKIZgAPBV4/hRw4aEbOOe2OufmB57vBVYA+lsvjaJdy3jO7NuB1+bnU1ahGx5JeHhtfj7RUcb5J2iqYqm/UCkEHZxzW6H6Bz9wxJk0zKwrMBiYXWPxLWa22MyeqO2Qg0h9XT4sg537y/hINzySMFAVmKr4tF4pJLeM9x1HwlDQCoGZfWhmS2v5mFDPz9MSeBX4oXOuOLD4EaAHMAjYCvz1CPvfYGY5ZpZTWFh4bG9GIsIpmcl0aB3PyzpsIGHgi/U72VpUqrkH5JjFBOsPcs6debh1ZrbdzFKdc1vNLBWo9VcyM4ulugw865x7rcbn3l5jm0eBt4+QYxIwCSA7O1vz08phxURHccmQdP75yVq2F5fSoXUz35FEDus/szbSMj6GsVkdfEeRMBUqhwwmA9cGnl8LvHnoBmZmwOPACufcfYesS63x8iJgaSPllAgzMTuDKld9bFYkVH26ppD3lm7ju6d011TFcsxCpRDcC4w1szXA2MBrzKyTmX15xcDJwNXAGbVcXvgnM1tiZouB04EfBTm/NFHdklvohkcS0g5WVHLHm8vo2i6B753W3XccCWNBO2RwJM65ncCYWpZvAc4JPJ8J1DrtlnPu6kYNKBFtYnY6P31lMfM27iZbU8FKiHl0xjrW79jPU9cN1+iAHJdQGSEQCVnnDEilRVy05iSQkLN5Vwl//yiXcwZ05LReKb7jSJhTIRA5ihbxMZw3sBNvL97K/oMVvuOIfOXut5YRHWX85ryso28schQqBCJ1cNmwdErKKnlniW54JKHhg+Xb+XBFAbeOySS1TXPfcaQJUCEQqYMhndvSPaWF5iSQkHCgrJK7Ji8js31LrhvVzXccaSJUCETqwMy4LDuDuRt2s65wn+84EuEenp5L/p4D3HNhf2Kj9d+4NAz9TRKpo4sHpxEdZbw8L893FIlgawv3MWnGOi4anMaJ3dv5jiNNiAqBSB21b92M03un8Oq8PCoqdcMjCT7nHHe+uYz42Ch+cU4f33GkiVEhEKmHidkZFOw9yIw1ug+GBN/7S7cxM3cHPzmrN+1baSptaVgqBCL1cEaf9iS3jOPFuTq5UILLOcfDH+fSPaUF3zyxi+840gSpEIjUQ2x0FBcNTmPaigJ27DvoO45EkNnrd7E0v5jvjOpOdFStk7aKHBcVApF6mpidQUWV440FuuGRBM/jM9fTNiGWi4fo9sbSOFQIROqpV4dWDMpI5N+fb6CkTDMXSuNbv2M/H67YzjdP7KL7FUijUSEQOQa3j+9D3u4D/GXKat9RJAI8+dl6YqOiuPoknTsgjUeFQOQYnNi9HVef2IUnP19PzoZdvuNIE1ZUUs7LOXlcMKiTriyQRqVCIHKMfj6+D53aNOdnryymtLzSdxxpop6bs4kD5ZVcrymKpZGpEIgco5bxMfzx4gGs27Gf+z9c4zuONEHllVU89fkGTu7Zjr6prX3HkSZOhUDkOJzaK4XLszOYNGMtizbv8R1Hmph3l2xlW3Ep3xnV3XcUiQAqBCLH6Zfn9iWlVTw/e2UxZRWa0lgahnOOxz5dT4+UFpzWK8V3HIkAKgQix6lN81j+cNEAVm3fy0PTc33HkSZi7obdLMkv4rpR3YjSREQSBCoEIg1gTN8OXDioE/+YnsvyLcW+40gT8Nin66onIhqc7juKRAgVApEGcuf5/UhMiOWnryyiXHdDlOOwYcd+PghMRNQ8ThMRSXCoEIg0kLYt4rhnQn+WbSnmsU/X+44jYezfn28gJsq4WjcxkiBSIRBpQOMHpDKmT3se+TiXvaXlvuNIGCoqKeelnM1ccEIa7VtrIiIJHhUCkQZ265mZFJdW8J8vNvqOImHo6VkbKCnTREQSfCoEIg1sYHoip/ZK4fFP13OgTDMYSt0VlZQz6dN1jM3qQFYnTUQkwaVCINIIfnBGT3buL+P5OZt8R5Ew8tjMdewtreC2sb18R5EIpEIg0giGdU1iRLck/jVjLQcrNEogR7dz30GemLmecwemappi8UKFQKSR/OCMTLYXH+SVeXm+o0gY+NeMdRwor+RHZ2b6jiIRSoVApJGc3LMdJ2Qk8sjHazUvgRxRQXEpT32+gQsHp9GzfSvfcSRChUQhMLMkM/vAzNYEHtseZrsNZrbEzBaaWU599xcJJjPjB6f3JG/3ASYv3OI7joSwh6fnUlHluHWMRgfEn5AoBMDtwDTnXCYwLfD6cE53zg1yzmUf4/4iQTOmb3v6prbm4Y9zqaxyvuNICMrfc4Dn52zmsux0urRr4TuORLBQKQQTgKcCz58CLgzy/iKNwsy45fSerCvcz3tLt/qOIyHooY/WAHDLGRodEL9CpRB0cM5tBQg8tj/Mdg6YambzzOyGY9gfM7vBzHLMLKewsLCB4osc3rj+Heme0oKHPsrFOY0SyP/bsGM/L+XkceXwDNISm/uOIxEuaIXAzD40s6W1fEyox6c52Tk3BBgP3Gxmp9Y3h3NuknMu2zmXnZKie4xL44uOMm4e3ZOV2/YybUWB7zgSQh6ctoaYKOPm03v6jiISvELgnDvTOde/lo83ge1mlgoQeKz1f03n3JbAYwHwOjA8sKpO+4v4csGgTmQkNefv0zVKINVyC/by+sJ8rh3ZVfcskJAQKocMJgPXBp5fC7x56AZm1sLMWn35HDgLWFrX/UV8io2O4sbTerBo8x4+XbPDdxwJAX/7cA0JsdF879TuvqOIAKFTCO4FxprZGmBs4DVm1snM3g1s0wGYaWaLgDnAO86594+0v0gouXRoOh1axzNpxjrfUcSzpflFvLN4K98+uRvtWsb7jiMCQIzvAADOuZ3AmFqWbwHOCTxfB5xQn/1FQkl8TDTXnNSVP09ZxYqtxZqeNkI557hr8jKSWsTxXY0OSAgJlRECkYhw1YjONI+N5omZ631HEU/eWJhPzsbd/Ozs3rRpHus7jshXVAhEgigxIY5Lhqbx5sItFO496DuOBNm+gxX88d2VDExvw2XZGb7jiPwXFQKRIPv2yd0oq6zimS82+o4iQfb3aWso2HuQuy/oR1SU+Y4j8l9UCESCrEdKS8b0ac8zX2yktFy3Ro4Uawv38cRn65k4NJ3BnXW7FQk9KgQiHlw/qhs795fx5sJ831EkCJxz3P3WcprFRPOzcX18xxGplQqBiAcn9WhH39TWPD5zvSYqigAfrihgxupCfji2FymtdJmhhCYVAhEPzIzrR3Vj9fZ9mqioiSstr+Set5eT2b4l15zUxXcckcNSIRDx5PwTUkluGc/jugSxSXt0xjo27Srh7gv6ERut/3IldOlvp4gn8THRXHtSFz5ZXcia7Xt9x5FGkL/nAA9/nMs5Azoysmey7zgiR6RCIOLRVSd2IT4miic+0yhBU/SHd1YA8KtzszwnETk6FQIRj5JaxHHxkDRem5/Prv1lvuNIA5q1difvLNnKTaN7kpbY3HcckaNSIRDx7LqTu3GwoopnNVFRk1FZ5fjt28tJS2zODbpfgYQJFQIRzzI7tOK0Xik8/cVGDlZooqKm4KWczazYWswvz+lLs9ho33FE6kSFQCQEfOeUbhTuPcgLczb7jiLHqbi0nL9MWcXwrkmcM6Cj7zgidaZCIBICRvVM5pTMZP743gpW64qDsPbQR7nsKinjjvOzMNP9CiR8qBCIhAAz477LBtEyPpabn53PgTIdOghH63fs58nA/Qr6p7XxHUekXlQIREJESqt4/nb5Cawp2Mdv317mO44cg9+/s4L4mGh+cnZv31FE6k2FQCSEnJKZwvdH9+D5OZt5a9EW33GkHmau2cGHK7Zz8+k9ad+qme84IvWmQiASYm4b24shnRP5xWtL2LSzxHccqYOKyip++/YyOiclcN2orr7jiBwTFQKREBMbHcWDVw4myuCW5+dTVlHlO5IcxfNzNrF6+z5+eU5f4mN0maGEJxUCkRCU3jaBP106kMV5Rfx5ykrfceQIikrKue+D1ZzYPYmz+3XwHUfkmKkQiISocf1TufrELjz66Xo+Wrnddxw5jPunraboQDl3nNdPlxlKWFMhEAlhvzq3L306tuLHLy2ioLjUdxw5RG7BXv4zayOXD+tMVqfWvuOIHBcVApEQ1iw2moe+MYT9ZZX84d0VvuNIDc45fvv2CprHRfOTs3r5jiNy3FQIREJcz/YtueGU7ryxcAs5G3b5jiMB01cVMGN1IbeOyaRdy3jfcUSOmwqBSBi46fQedGzdjLveWkZllfMdJ+KVVVTxu7dX0D2lBdec1NV3HJEGoUIgEgYS4mL4xTl9WJpfzMs5ugGSb0/P2sC6Hfv5zblZxMXov1FpGvQ3WSRMXHBCJ4Z1bcufpqyi6EC57zgRa8e+gzzw4RpG907h9D7tfccRaTAhUQjMLMnMPjCzNYHHtrVs09vMFtb4KDazHwbW3WVm+TXWnRP0NyHSyMyMuy7ox56SMh74cI3vOBHrr1NXcaC8kl+fm+U7ikiDColCANwOTHPOZQLTAq//i3NulXNukHNuEDAUKAFer7HJ375c75x7NxihRYKtX6c2XDG8M0/N2sAa3SY56JbmF/HC3M1cO7IrPdu39B1HpEGFSiGYADwVeP4UcOFRth8DrHXObWzMUCKh6Cdn9aZFXDR3v7Uc53SCYbBUX2a4nLYJcfzPmEzfcUQaXKgUgg7Oua0AgcejHZi7Anj+kGW3mNliM3uitkMOXzKzG8wsx8xyCgsLjy+1iAdJLeL48Vm9mZm7g6nLNYNhsLy7ZBtz1u/ix2f1ok3zWN9xRBpc0AqBmX1oZktr+ZhQz88TB1wAvFxj8SNAD2AQsBX46+H2d85Ncs5lO+eyU1JS6v9GRELAVSM607tDK373znJKyyt9x2nySsurJ4bq07EVVwzr7DuOSKMIWiFwzp3pnOtfy8ebwHYzSwUIPBYc4VONB+Y757761cg5t905V+mcqwIeBYY35nsR8S0mOoo7z89i864DPPbpOt9xmrxJM9aRv+cAd57fj+go3a9AmqZQOWQwGbg28Pxa4M0jbHslhxwu+LJMBFwELG3QdCIhaGTPZM4Z0JGHp69l864S33GarAWbdvPQR7mcM6AjJ/Vo5zuOSKMJlUJwLzDWzNYAYwOvMbNOZvbVFQNmlhBY/9oh+//JzJaY2WLgdOBHwYkt4tevzs0iOsr46SuLqNIMhg2uoLiUG5+ZR4c28fz+wgG+44g0qhjfAQCcczupvnLg0OVbgHNqvC4BvlbRnXNXN2pAkRCVlticO87P4mevLObJzzdw/ahuviM1GQcrKrnxmXkUH6jgtZtG0rZFnO9IIo0qVEYIROQYTRyazpl92/On91eSW6C5CRrKXZOXM3/THv4y8QT6purWxtL0qRCIhDkz4w8XDyAhLprbXlpEeWWV70hh79nZG3l+ziZuGt2DcwemHn0HkSZAhUCkCWjfqhl/uGgAi/OKeHh6ru84YW3uhl3cNXkZo3un8OOzevuOIxI0KgQiTcT4AalcOKgTD32Uy+K8Pb7jhKWtRQf4/jPzSW+bwANXDNYlhhJRVAhEmpC7L+hPcst4bntpkSYsqqfS8kpu/M88DpRVMOnqoZqNUCKOCoFIE9ImIZY/XTqQ3IJ9/GXKKt9xwoZzjl++toRFeUXcd/kgMju08h1JJOhUCESamFN7pXD1iV14/LP1zFq703ecsPDnKat4bUE+t43txdn9OvqOI+KFCoFIE/SLc/rQJSmBn7y8iILiUt9xQtq/P1vPPz5eyzdGdOYHZ/T0HUfEGxUCkSYoIS6G+68YzO6SMq6Y9AXbVQpq9fbiLdz99nLOyurAPRP6Y6aTCCVyqRCINFGDMhJ56rrhbC8u5UqVgq/5fO0ObntxEdld2vLglbqiQESFQKQJG9Y16atScMWkL9hWpFIAsGxLETc8PY+uyQk8ds0wmsVG+44k4p0KgUgTlx0oBQXFpVwxaRZbiw74juTV5l0lfOvJubRqFsNT1w2nTYIuLxQBFQKRiJDdNYmnrx/Ojn3V5xREainYue8g1zwxh7KKKp6+bjipbZr7jiQSMlQIRCLE0C7VIwU7A6Vgy57IKgUFe0u56rHZbNlzgMevzdZcAyKHUCEQiSBDu7Tl6euHs2tfGRP/OYsFm3b7jhQUm3aWcOkjs9i0q4THrs0mu2uS70giIUeFQCTCDOnclme/OwKAif+cxT8/WUtVlfOcqvGs2FrMJf/8nOLScp79zghOyUzxHUkkJKkQiESggemJvPs/pzA2qwP3vreSa5+cQ+Heg75jNbi5G3Zx2b9mEW3Gy987icGd2/qOJBKyVAhEIlSbhFj+cdUQfndhf+as38X4Bz5lxupC37EazPSVBVz9+GxSWsbzyvdP0jkDIkehQiASwcyMb57Yhcm3jKJtQizXPDGHe99bSXllle9ox+WNBfl89+kcerZvyUs3nkR62wTfkURCngqBiNC7Yysm3zKKK4dn8M9P1nLlpC/Yvb/Md6xj8vjM9fzwxYVkd23L8989keSW8b4jiYQFFQIRAaB5XDR/vHggD1wxiMV5RVzyz8/ZvKvEd6w6q6xy3DV5Gfe8vZxx/Try728Pp1UzTTokUlcqBCLyXyYMSuPp64dTuPcglzzyOcu3FPuOdFQHyiq56dl5/PvzDVx3cjcevmqIpiMWqScVAhH5mhO7t+OVG0cSZcbl/5rF52t3+I50WDv3HeTKR79g6vLt3Hl+Fnecn6UbFYkcAxUCEalV746teO2mkXRs04xvPTGXtxZt8R3pa9YV7uPiRz5nxdZiHrlqKN8+uZvvSCJhS4VARA6rU2JzXrlxJIMyEvnB8wt4fOZ635G+Mm/jLi555HP2llbw/A0nMq5/R9+RRMKaCoGIHFGbhFievn444/p15J63l3PbiwvJLdjnLU9JWQX/+DiXKx+dTWJCHK/fNJIhmnBI5LjF+A4gIqGvWWw0D181hD9PWcWTn63ntQX5nNm3Pd89pTvDuyVh1vjH7EvLK3l29iYe+TiXHfvKGNOnPX+eeAJJLeIa/c8WiQTmXNOdw/xosrOzXU5Oju8YImFlx76D/GfWRv7zxUZ27S/jhPQ2fPfU7ozr15GY6IYfdCyrqOKlnM089FEu24pLGdmjHT8+qxdDu+gGRSL1ZWbznHPZta4LhUJgZhOBu4C+wHDnXK0/pc1sHPAAEA085py7N7A8CXgR6ApsAC5zzh31Nm4qBCLH7kBZJa/Oz+PxmetZv2M/6W2bc/HgNE7tlcIJGYnEHmc52LnvINNWFvD3j9awedcBhnZpy4/P6sXIHskN9A5EIk84FIK+QBXwL+AntRUCM4sGVgNjgTxgLnClc265mf0J2OWcu9fMbgfaOud+frQ/V4VA5PhVVTk+XLGdx2euZ+6GXVQ5aBUfw0k92nFKrxROzUymS7sWR/wcxaXlLM0rYlFeEUvy97BocxH5ew4A0D+tNT8+qzeje6UE5dCESFN2pEIQEucQOOdWAEf7xz4cyHXOrQts+wIwAVgeeBwd2O4p4GPgqIVARI5fVJRxVr+OnNWvI0Ul5Xy+dgcz1uxgxupCpi7fDkB62+a0axFHlaueUbDKVX9UVjlKy6u++uEP0DkpgcGdE7l2ZBcGd25Ldpe2KgIiQRAShaCO0oDNNV7nASMCzzs457YCOOe2mln7w30SM7sBuAGgc+fOjRRVJDK1SYhl/IBUxg9IxTnH+h37+XTNDr5Yt5OSskqio4wogyizwHMjJtq4sn0GA9MTGZjehsQEnSQo4kPQCoGZfQjUdqHwr5xzb9blU9SyrN7HO5xzk4BJUH3IoL77i0jdmBndU1rSPaUl147s6juOiBxF0AqBc+7M4/wUeUBGjdfpwJdTp203s9TA6EAqUHCcf5aIiEhECaeJieYCmWbWzczigCuAyYF1k4FrA8+vBeoy4iAiIiIBIVEIzOwiM8sDTgLeMbMpgeWdzOxdAOdcBXALMAVYAbzknFsW+BT3AmPNbA3VVyHcG+z3ICIiEs5C4rJDX3TZoYiIRJIjXXYYEiMEIiIi4pcKgYiIiKgQiIiIiAqBiIiIoEIgIiIiqBCIiIgIKgQiIiKCCoGIiIigQiAiIiJE+EyFZlYIbGzAT5kM7GjAzxeOIv1rEOnvH/Q1AH0NIv39Q+h+Dbo451JqWxHRhaChmVnO4aaEjBSR/jWI9PcP+hqAvgaR/v4hPL8GOmQgIiIiKgQiIiKiQtDQJvkOEAIi/WsQ6e8f9DUAfQ0i/f1DGH4NdA6BiIiIaIRAREREVAgajJmNM7NVZpZrZrf7zhNsZvaEmRWY2VLfWXwwswwzm25mK8xsmZnd6jtTsJlZMzObY2aLAl+Du31n8sHMos1sgZm97TuLD2a2wcyWmNlCM8vxnSfYzCzRzF4xs5WB/w9O8p2prnTIoAGYWTSwGhgL5AFzgSudc8u9BgsiMzsV2Ac87Zzr7ztPsJlZKpDqnJtvZq2AecCFEfZ3wIAWzrl9ZhYLzARudc594TlaUJnZbUA20No5d57vPMFmZhuAbOdcKF6D3+jM7CngU+fcY2YWByQ45/Z4jlUnGiFoGMOBXOfcOudcGfACMMFzpqByzs0AdvnO4Ytzbqtzbn7g+V5gBZDmN1VwuWr7Ai9jAx8R9RuHmaUD5wKP+c4iwWdmrYFTgccBnHNl4VIGQIWgoaQBm2u8ziPCfhjI/zOzrsBgYLbnKEEXGC5fCBQAHzjnIu1rcD/wM6DKcw6fHDDVzOaZ2Q2+wwRZd6AQeDJw2OgxM2vhO1RdqRA0DKtlWUT9ZiTVzKwl8CrwQ+dcse88weacq3TODQLSgeFmFjGHj8zsPKDAOTfPdxbPTnbODQHGAzcHDidGihhgCPCIc24wsB8Im3PKVAgaRh6QUeN1OrDFUxbxJHDc/FXgWefca77z+BQYJv0YGOc3SVCdDFwQOIb+AnCGmT3jN1LwOee2BB4LgNepPqQaKfKAvBojY69QXRDCggpBw5gLZJpZt8BJJFcAkz1nkiAKnFD3OLDCOXef7zw+mFmKmSUGnjcHzgRWeg0VRM65Xzjn0p1zXan+P+Aj59w3PccKKjNrETiplsBQ+VlAxFx55JzbBmw2s96BRWOAsDmxOMZ3gKbAOVdhZrcAU4Bo4Ann3DLPsYLKzJ4HRgPJZpYH3Omce9xvqqA6GbgaWBI4hg7wS+fcu/4iBV0q8FTgqpso4CXnXEReehfBOgCvV/djYoDnnHPv+40UdD8Ang38crgO+LbnPHWmyw5FREREhwxEREREhUBERERQIRARERFUCERERAQVAhEREUGFQESCzMxuCNwRb7+ZvW1mnXxnEhEVAhEJIjObAPwNuI3qGdwcMNnM9H+RiGeah0BEgsbM5gKfOOd+EnidRPU03xdH2CROIiFHrVxEgsLM2gLZwDtfLnPO7QLmAWN95RKRaioEIhIs3QKPuYcsXwP0CHIWETmE7mUgIsHy5X3hVwTmuv9SHDAj+HFEpCYVAhEJlv2BxzHAzhrL/wA0C34cEalJhUBEgmVd4HGvc+6rwwZmFg+s9RNJRL6kcwhEJCicc3uAucApXy4LlIFTgameYolIgEYIRCSY7gGeMLNNVJ9MeEfg8X2vqURE8xCISHCZ2XeBXwMpwDTge865LX5TiYgKgYiIiOgcAhEREVEhEBEREVQIREREBBUCERERQYVAREREUCEQERERVAhEREQEFQIRERFBhUBERESA/wNRaEcovr3AdgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import qutip as qt\n",
    "from math import pi\n",
    "from itertools import *\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from qiskit import QuantumCircuit, execute\n",
    "from qiskit import Aer, IBMQ, transpile\n",
    "from qiskit.providers.ibmq.managed import IBMQJobManager\n",
    "\n",
    "#IBMQ.save_account('YOUR_TOKEN')\n",
    "\n",
    "samples = 60\n",
    "n_shots = 4096\n",
    "classical = False # if True, uses \"ibmq_qasm_simulator\" \n",
    "quantum = False # if True, uses \"ibmq_vigo\"\n",
    "\n",
    "# Or: set local to False, and specify particular backend.\n",
    "local = True \n",
    "backend_name = \"ibmq_5_yorktown\" #\"ibmq_16_melbourne\", etc.\n",
    "\n",
    "##############################################################################\n",
    "\n",
    "# For each Pauli string of the form \"XIXXY\", of length n,\n",
    "# construct the tensor product of the corresponding\n",
    "# Pauli operators. These form a basis for 2^n x 2^n Hermitian matrices.\n",
    "def construct_pauli_basis(n):\n",
    "    IXYZ = {\"I\": qt.identity(2),\\\n",
    "            \"X\": qt.sigmax(),\\\n",
    "            \"Y\": qt.sigmay(),\\\n",
    "            \"Z\": qt.sigmaz()}\n",
    "    return dict([(\"\".join(p),\\\n",
    "                  qt.tensor(*[IXYZ[s] for s in p]))\\\n",
    "                for p in product(IXYZ.keys(), repeat=n)])\n",
    "\n",
    "# Find the components of an operator in the Pauli basis.\n",
    "def op_pauli(O, basis):\n",
    "    return dict([(s, (o.dag()*O).tr()/np.sqrt(len(basis)))\\\n",
    "                    for s, o in basis.items()])\n",
    "\n",
    "# Recovers the original matrix from the Pauli components.\n",
    "def pauli_op(P, basis):\n",
    "    return sum([P[s]*o for s, o in basis.items()])\n",
    "\n",
    "# Pretty print for the Pauli decomposition.\n",
    "def print_pauli(P):\n",
    "    for s, v in P.items():\n",
    "        if not np.isclose(v, 0):\n",
    "            print(\"%s: %.2f\" % (s, v))\n",
    "\n",
    "##############################################################################\n",
    "\n",
    "if quantum:\n",
    "    local = False\n",
    "    backend_name = \"ibmq_vigo\"\n",
    "    \n",
    "if classical:\n",
    "    local = False\n",
    "    backed_name = \"ibmq_qasm_simulator\"\n",
    "\n",
    "##############################################################################\n",
    "\n",
    "print(\"our matrix:\")\n",
    "H = qt.Qobj(np.array([[1,0,0,0],\\\n",
    "                      [0,0,-1,0],\\\n",
    "                      [0,-1,0,0],\\\n",
    "                      [0,0,0,1]]))\n",
    "H.dims = [[2]*2, [2]*2]\n",
    "print(H)\n",
    "Hl, Hv = H.eigenstates()\n",
    "\n",
    "# Let's make sure we know the right answer.\n",
    "print(\"\\ncompare:\")\n",
    "print(\"lowest eigval =  %.4f\" % (Hl[0]))\n",
    "print(\"eigenvector = \\n%s\" % (Hv[0]))\n",
    "print()\n",
    "\n",
    "basis = construct_pauli_basis(2)\n",
    "H_pauli = op_pauli(H, basis)\n",
    "\n",
    "##############################################################################\n",
    "\n",
    "def ansatz(theta):\n",
    "    circ = QuantumCircuit(2)\n",
    "    circ.h(0)\n",
    "    circ.cx(0,1)\n",
    "    circ.rx(theta, 0)\n",
    "    return circ\n",
    "\n",
    "##############################################################################\n",
    "\n",
    "if local:\n",
    "    backend = Aer.get_backend('qasm_simulator')\n",
    "else:    \n",
    "    provider = IBMQ.load_account()\n",
    "    job_manager = IBMQJobManager()\n",
    "    backend = provider.get_backend(backend_name)\n",
    "\n",
    "print(\"laying out circuits...\")\n",
    "# Creates circuits: for each theta, we run through all the operators \n",
    "# in the Pauli decomposition, and create a circuit for each one\n",
    "# that starts with the ansatz, then consists of pre-measurement\n",
    "# rotations on each of the qubits, and then a full measurement.\n",
    "circs = []\n",
    "circ_indices = []\n",
    "thetas = np.linspace(0, 2*np.pi, samples)\n",
    "for t in thetas:\n",
    "    for n, v in H_pauli.items():\n",
    "        if not np.isclose(v, 0):\n",
    "            circ = ansatz(t)\n",
    "            for i, o in enumerate(n):\n",
    "                if o == \"X\":\n",
    "                    circ.ry(-pi/2, i)\n",
    "                elif o == \"Y\":\n",
    "                    circ.rx(pi/2, i)\n",
    "            circ.measure_all()\n",
    "            circs.append(circ)\n",
    "            circ_indices.append((t, n))\n",
    "\n",
    "print(\"evaluating...\\n\")\n",
    "# We optimize the circuits, and send them off to be evaluated.\n",
    "circs = transpile(circs, backend=backend)\n",
    "jb = execute(circs, backend, shots=n_shots) if local \\\n",
    "        else job_manager.run(circs, backend=backend, name='vqe', shots=n_shots) \n",
    "jb_results = jb.result() if local else jb.results()\n",
    "\n",
    "# We run through the results. For each Pauli operator P in the\n",
    "# decomposition, we get the counts which, when converted into probabilities,\n",
    "# weight the eigenvalues, which are the product of the {-1, 1}'s of each\n",
    "# qubit. We weight each of these ...Z... expectation values by\n",
    "# the component of P in the Pauli decomposition, and sum to get the \n",
    "# full expectation value.\n",
    "# We do this for each theta, and find the theta with the least <H>.\n",
    "expectation_values = []\n",
    "for t in thetas:\n",
    "    values = []\n",
    "    for n, v in H_pauli.items():\n",
    "        if not np.isclose(v, 0):\n",
    "            if n.count(\"I\") == len(n):\n",
    "                values.append(H_pauli[n])\n",
    "            else:\n",
    "                counts = jb_results.get_counts(circ_indices.index((t, n)))\n",
    "                values.append(\\\n",
    "                    H_pauli[n]*\\\n",
    "                        sum([(c/n_shots)*\\\n",
    "                                np.prod([1 if s_=='1' \\\n",
    "                                            else -1 for s_ in s])\\\n",
    "                        for s, c in counts.items()])\\\n",
    "                    )\n",
    "    e = sum(values)\n",
    "    print(\"theta = %.4f, <H> = %.8f\" % (t, e))\n",
    "    expectation_values.append(e)\n",
    "min_result = np.argmin(expectation_values)\n",
    "final_theta = thetas[min_result]\n",
    "print(\"\\nfinal_theta: %.4f, <H> = %.4f\" % (final_theta, expectation_values[min_result]))\n",
    "\n",
    "# Use a statevector simulator to look at our final state.\n",
    "state_backend = Aer.get_backend('statevector_simulator')\n",
    "state = execute(ansatz(final_theta), backend=state_backend).result().get_statevector()\n",
    "print(\"final vector = \\n%s\" % qt.Qobj(state))\n",
    "\n",
    "# Plot!\n",
    "fig = plt.figure(figsize=(8,6))\n",
    "ax = fig.add_subplot(111)\n",
    "ax.plot(thetas, expectation_values)\n",
    "ax.set_xlabel('θ', fontsize=14)\n",
    "ax.set_ylabel('<H>', fontsize=14)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can check out the transpiled circuits:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[        ┌───┐      ░ ┌─┐   \n",
       "    q_0: ┤ H ├──■───░─┤M├───\n",
       "         └───┘┌─┴─┐ ░ └╥┘┌─┐\n",
       "    q_1: ─────┤ X ├─░──╫─┤M├\n",
       "              └───┘ ░  ║ └╥┘\n",
       " meas: 2/══════════════╩══╩═\n",
       "                       0  1 ,\n",
       "         ┌───┐     ┌────────────┐ ░ ┌─┐   \n",
       "    q_0: ┤ H ├──■──┤ U2(pi,-pi) ├─░─┤M├───\n",
       "         └───┘┌─┴─┐├────────────┤ ░ └╥┘┌─┐\n",
       "    q_1: ─────┤ X ├┤ U2(pi,-pi) ├─░──╫─┤M├\n",
       "              └───┘└────────────┘ ░  ║ └╥┘\n",
       " meas: 2/════════════════════════════╩══╩═\n",
       "                                     0  1 ,\n",
       "         ┌───┐     ┌────────────────┐ ░ ┌─┐   \n",
       "    q_0: ┤ H ├──■──┤ U2(-pi/2,pi/2) ├─░─┤M├───\n",
       "         └───┘┌─┴─┐├────────────────┤ ░ └╥┘┌─┐\n",
       "    q_1: ─────┤ X ├┤ U2(-pi/2,pi/2) ├─░──╫─┤M├\n",
       "              └───┘└────────────────┘ ░  ║ └╥┘\n",
       " meas: 2/════════════════════════════════╩══╩═\n",
       "                                         0  1 ,\n",
       "         ┌───┐      ░ ┌─┐   \n",
       "    q_0: ┤ H ├──■───░─┤M├───\n",
       "         └───┘┌─┴─┐ ░ └╥┘┌─┐\n",
       "    q_1: ─────┤ X ├─░──╫─┤M├\n",
       "              └───┘ ░  ║ └╥┘\n",
       " meas: 2/══════════════╩══╩═\n",
       "                       0  1 ,\n",
       "         ┌───┐     ┌────────────────────────┐ ░ ┌─┐   \n",
       "    q_0: ┤ H ├──■──┤ U3(0.10649,-pi/2,pi/2) ├─░─┤M├───\n",
       "         └───┘┌─┴─┐└────────────────────────┘ ░ └╥┘┌─┐\n",
       "    q_1: ─────┤ X ├───────────────────────────░──╫─┤M├\n",
       "              └───┘                           ░  ║ └╥┘\n",
       " meas: 2/════════════════════════════════════════╩══╩═\n",
       "                                                 0  1 ,\n",
       "         ┌───┐     ┌────────────────┐ ░ ┌─┐   \n",
       "    q_0: ┤ H ├──■──┤ U2(-3.0351,pi) ├─░─┤M├───\n",
       "         └───┘┌─┴─┐└─┬────────────┬─┘ ░ └╥┘┌─┐\n",
       "    q_1: ─────┤ X ├──┤ U2(pi,-pi) ├───░──╫─┤M├\n",
       "              └───┘  └────────────┘   ░  ║ └╥┘\n",
       " meas: 2/════════════════════════════════╩══╩═\n",
       "                                         0  1 ,\n",
       "         ┌───┐     ┌───────────────────────┐ ░ ┌─┐   \n",
       "    q_0: ┤ H ├──■──┤ U3(1.6773,-pi/2,pi/2) ├─░─┤M├───\n",
       "         └───┘┌─┴─┐└───┬────────────────┬──┘ ░ └╥┘┌─┐\n",
       "    q_1: ─────┤ X ├────┤ U2(-pi/2,pi/2) ├────░──╫─┤M├\n",
       "              └───┘    └────────────────┘    ░  ║ └╥┘\n",
       " meas: 2/═══════════════════════════════════════╩══╩═\n",
       "                                                0  1 ,\n",
       "         ┌───┐     ┌────────────────────────┐ ░ ┌─┐   \n",
       "    q_0: ┤ H ├──■──┤ U3(0.10649,-pi/2,pi/2) ├─░─┤M├───\n",
       "         └───┘┌─┴─┐└────────────────────────┘ ░ └╥┘┌─┐\n",
       "    q_1: ─────┤ X ├───────────────────────────░──╫─┤M├\n",
       "              └───┘                           ░  ║ └╥┘\n",
       " meas: 2/════════════════════════════════════════╩══╩═\n",
       "                                                 0  1 ,\n",
       "         ┌───┐     ┌────────────────────────┐ ░ ┌─┐   \n",
       "    q_0: ┤ H ├──■──┤ U3(0.21299,-pi/2,pi/2) ├─░─┤M├───\n",
       "         └───┘┌─┴─┐└────────────────────────┘ ░ └╥┘┌─┐\n",
       "    q_1: ─────┤ X ├───────────────────────────░──╫─┤M├\n",
       "              └───┘                           ░  ║ └╥┘\n",
       " meas: 2/════════════════════════════════════════╩══╩═\n",
       "                                                 0  1 ,\n",
       "         ┌───┐     ┌────────────────┐ ░ ┌─┐   \n",
       "    q_0: ┤ H ├──■──┤ U2(-2.9286,pi) ├─░─┤M├───\n",
       "         └───┘┌─┴─┐└─┬────────────┬─┘ ░ └╥┘┌─┐\n",
       "    q_1: ─────┤ X ├──┤ U2(pi,-pi) ├───░──╫─┤M├\n",
       "              └───┘  └────────────┘   ░  ║ └╥┘\n",
       " meas: 2/════════════════════════════════╩══╩═\n",
       "                                         0  1 ]"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "[circ.draw() for circ in circs[:10]]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<hr>\n",
    "\n",
    "So in the end, the answer is: when $\\theta = \\pi$, the ansatz is $-\\frac{i}{\\sqrt{2}} \\begin{pmatrix} 0 \\\\ 1 \\\\ 1 \\\\ 0 \\end{pmatrix}$, which corresponds (up to phase) to the lowest eigenstate of $H$ with eigenvalue $-1$. We could imagine more sophisticated schemes, involving more elaborate ansatze, using fancier algorithms to search through the parameters, from gradient descent onward, generating thousands upon thousands of circuits that when implemented thousands upon thousands of times allow one to estimate the lowest eigenvalue/eigenvector of a matrix faster than a classical computer.\n",
    "\n",
    "Special thanks to <a href=\"https://qosf.org/\">qosf</a>, and <a href=\"https://www.mustythoughts.com/variational-quantum-eigensolver-explained\">this</a>."
   ]
  }
 ],
 "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": 4
}