{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Using the Hadamard Test to Determine Quantum Krylov Subspace Decomposition Matrix Elements\n",
    "\n",
    "The Hadamard test, is a quantum algorithm for estimating the expectation values and is  a useful subroutine for a number of quantum applications ranging from estimation of molecular ground state energies to quantum semidefinite programming. This tutorial will briefly introduce the Hadamard test, and demonstrate how it can be implemented in CUDA-Q and then parallelized for a Quantum Krylov Subspace Diagonalization application.\n",
    "\n",
    "The Hadamard test is performed using a register with an ancilla qubit in the $\\ket{0}$ state and a prepared quantum state $\\ket{\\psi}$, the following circuit can be used to extract the expectation value from measurement of the ancilla.\n",
    "\n",
    "\n",
    "![Htest](./images/htest.png)\n",
    "\n",
    "The key insight is to note that $$P(0) = \\frac{1}{2} \\left[ I + Re \\bra{\\psi} O \\ket{\\phi} \\right]$$ and $$P(1) = \\frac{1}{2} \\left[ I - Re \\bra{\\psi} O \\ket{\\phi} \\right]$$ so their difference is equal to $$P(0)-P(1) = Re \\bra{\\psi} O \\ket{\\phi}$$\n",
    "\n",
    "\n",
    "More details and a short derivation can be found [here](https://en.wikipedia.org/wiki/Hadamard_test)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What if you want to perform the Hadamard test to compute an expectation value like $\\bra{\\psi} O \\ket{\\phi}$, where $\\ket{\\psi}$ and $\\ket{\\phi}$ are different states and $O$ is a Pauli Operator? This is a common subroutine for the QKSD, where matrix elements are determined by computing expectation values between different states.\n",
    "\n",
    "Defining $O$ as \n",
    "$$O = X_1X_2$$\n",
    "\n",
    "and given the fact that\n",
    "$$\\ket{\\psi} = U_{\\psi}\\ket{0} \\qquad \\ket{\\phi} = U_{\\phi}\\ket{0}$$\n",
    "\n",
    "We can combine the state preparation steps into the operator resulting in\n",
    "$$\\bra{\\psi}O\\ket{\\phi} = \\bra{0}U_\\psi^\\dagger O U_\\phi\\ket{0}$$\n",
    "Which corresponds to the following circuit\n",
    "![Htest2](./images/htestfactored.png)\n",
    "\n",
    "By preparing this circuit, and repeatedly measuring the ancilla qubit, allows for estimation of the expectation value as $$P(0)-P(1) = Re \\bra{\\psi} O \\ket{\\phi}$$\n",
    "\n",
    "\n",
    "The following sections demonstrate how this can be performed in CUDA-Q."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Numerical result as a reference: "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Before performing the Hadamard test, lets determine the exact expectation value by performing the matrix multiplications explicitly. The code below builds two CUDA-Q kernels corresponding to $\\ket{\\psi} = \\frac{1}{\\sqrt{2}}\\begin{pmatrix}1 \\\\ 0 \\\\ 1 \\\\ 0\\end{pmatrix}$ and $\\ket{\\phi} = \\begin{pmatrix}0 \\\\ 1 \\\\ 0 \\\\ 0\\end{pmatrix}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import cudaq\n",
    "import numpy as np\n",
    "from functools import reduce\n",
    "\n",
    "cudaq.set_target('nvidia')\n",
    "\n",
    "qubit_num = 2\n",
    "\n",
    "\n",
    "@cudaq.kernel\n",
    "def psi(num: int):\n",
    "    q = cudaq.qvector(num)\n",
    "    h(q[1])\n",
    "\n",
    "\n",
    "@cudaq.kernel\n",
    "def phi(n: int):\n",
    "    q = cudaq.qvector(n)\n",
    "    x(q[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The state vectors can be accessed using the `get_state` command and printed as numpy arrays"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Psi state:  SV: [(0.707107,0), (0,0), (0.707107,0), (0,0)]\n",
      "\n",
      "Phi state:  SV: [(0,0), (1,0), (0,0), (0,0)]\n",
      "\n"
     ]
    }
   ],
   "source": [
    "psi_state = cudaq.get_state(psi, qubit_num)\n",
    "print('Psi state: ', psi_state)\n",
    "\n",
    "phi_state = cudaq.get_state(phi, qubit_num)\n",
    "print('Phi state: ', phi_state)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Hamiltonian operator ($O$ in the derivation above) is defined as a CUDA-Q spin operator and converted to a matrix withe the `to_matrix`. The following line of code performs the explicit matrix multiplications to produce the exact expectation value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "hamiltonian:  [[0.+0.j 0.+0.j 0.+0.j 1.+0.j]\n",
      " [0.+0.j 0.+0.j 1.+0.j 0.+0.j]\n",
      " [0.+0.j 1.+0.j 0.+0.j 0.+0.j]\n",
      " [1.+0.j 0.+0.j 0.+0.j 0.+0.j]] \n",
      "\n",
      "Numerical expectation value:  (0.7071067690849304+0j)\n"
     ]
    }
   ],
   "source": [
    "ham = cudaq.spin.x(0) * cudaq.spin.x(1)\n",
    "ham_matrix = ham.to_matrix()\n",
    "print('hamiltonian: ', np.array(ham_matrix), '\\n')\n",
    "\n",
    "exp_val = reduce(np.dot, (np.array(psi_state).conj().T, ham_matrix, phi_state))\n",
    "\n",
    "print('Numerical expectation value: ', exp_val)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Using `Sample` to perform the Hadamard test"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Three CUDA-Q kernels are constructed below corresponding to $\\ket{\\psi}$, $\\ket{\\phi}$, and the Hamiltonian. A fourth kernel constructs the Hadamard test circuit and completes with a measurement of the ancilla qubit in the computational basis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "import cudaq\n",
    "\n",
    "cudaq.set_target('nvidia')\n",
    "\n",
    "\n",
    "@cudaq.kernel\n",
    "def U_psi(q: cudaq.qview):\n",
    "    h(q[1])\n",
    "\n",
    "\n",
    "@cudaq.kernel\n",
    "def U_phi(q: cudaq.qview):\n",
    "    x(q[0])\n",
    "\n",
    "\n",
    "@cudaq.kernel\n",
    "def ham_cir(q: cudaq.qview):\n",
    "    x(q[0])\n",
    "    x(q[1])\n",
    "\n",
    "\n",
    "@cudaq.kernel\n",
    "def kernel(n: int):\n",
    "    ancilla = cudaq.qubit()\n",
    "    q = cudaq.qvector(n)\n",
    "    h(ancilla)\n",
    "    cudaq.control(U_phi, ancilla, q)\n",
    "    cudaq.control(ham_cir, ancilla, q)\n",
    "    cudaq.control(U_psi, ancilla, q)\n",
    "\n",
    "    h(ancilla)\n",
    "\n",
    "    mz(ancilla)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The CUDA-Q Sample method computes 100000 sample ancilla measurements and uses them to estimate the expectation value. The standard error is provided as well. Try increasing the sample size and note the convergence of the expectation value and the standard error towards the numerical result."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{ 0:85281 1:14719 }\n",
      "\n",
      "Observable QC:  0.70562 + - 0.0015844563982640861\n",
      "Numerical result 0.7071067690849304\n"
     ]
    }
   ],
   "source": [
    "shots = 100000\n",
    "qubit_num = 2\n",
    "count = cudaq.sample(kernel, qubit_num, shots_count=shots)\n",
    "print(count)\n",
    "\n",
    "mean_val = (count['0'] - count['1']) / shots\n",
    "error = np.sqrt(2 * count['0'] * count['1'] / shots) / shots\n",
    "print('Observable QC: ', mean_val, '+ -', error)\n",
    "print('Numerical result', np.real(exp_val))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Multi-GPU evaluation of QKSD matrix elements using the Hadamard Test\n",
    "\n",
    "This example is small, but a more practical application of the Hadamard test such as QKSD will require much larger circuits. The QKSD method works by reducing the exponential $2^N$ Hilbert space into an exponentially smaller subspace using a set of non-orthogonal states which are easy to prepare on a quantum computer. The Hadamard test is used to compute the matrix elements of this smaller subspace which is then diagonalized using a classical method to produce the eigenvalues. [This paper](https://www.osti.gov/servlets/purl/1962060) described the method in more detail and is the source of the figure below.\n",
    "\n",
    "![Htest3](./images/QKSD.png)\n",
    "\n",
    "This method can be easily parallelized and multiple QPUs, if available could compute the matrix elements asynchronously.  The CUDA-Q `mqpu` backend allows you to simulate a computation across multiple simulated QPUs. The code below demonstrates how.\n",
    "\n",
    "First, the Hadamard test circuit is defined, but this time the $\\ket{\\psi}$ and $\\ket{\\phi}$ states contain parameterized rotations so that multiple states can be quickly generated for the sake of example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "import cudaq\n",
    "\n",
    "\n",
    "@cudaq.kernel\n",
    "def U_psi(q: cudaq.qview, theta: float):\n",
    "    ry(theta, q[1])\n",
    "\n",
    "\n",
    "@cudaq.kernel\n",
    "def U_phi(q: cudaq.qview, theta: float):\n",
    "    rx(theta, q[0])\n",
    "\n",
    "\n",
    "@cudaq.kernel\n",
    "def ham_cir(q: cudaq.qview):\n",
    "    x(q[0])\n",
    "    x(q[1])\n",
    "\n",
    "\n",
    "@cudaq.kernel\n",
    "def kernel(n: int, angle: float, theta: float):\n",
    "    ancilla = cudaq.qubit()\n",
    "    q = cudaq.qvector(n)\n",
    "    h(ancilla)\n",
    "    cudaq.control(U_phi, ancilla, q, theta)\n",
    "    cudaq.control(ham_cir, ancilla, q)\n",
    "    cudaq.control(U_psi, ancilla, q, angle)\n",
    "\n",
    "    h(ancilla)\n",
    "\n",
    "    mz(ancilla)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, the `nvidia-mqpu` backend is specified and the number of GPUs available is determined."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of QPUs: 5\n"
     ]
    }
   ],
   "source": [
    "cudaq.set_target(\"nvidia-mqpu\")\n",
    "\n",
    "target = cudaq.get_target()\n",
    "qpu_count = target.num_qpus()\n",
    "print(\"Number of QPUs:\", qpu_count)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `sample_async` command is then used to distribute the Hadamard test computations across multiple simulated QPUs. The results are saved in a list and accessed below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "shots = 100000\n",
    "angle = [0.0, 1.5, 3.14, 0.7]\n",
    "theta = [0.6, 1.2, 2.2, 3.0]\n",
    "qubit_num = 2\n",
    "\n",
    "result = []\n",
    "for i in range(4):\n",
    "    count = cudaq.sample_async(kernel,\n",
    "                               qubit_num,\n",
    "                               angle[i],\n",
    "                               theta[i],\n",
    "                               shots_count=shots,\n",
    "                               qpu_id=i % qpu_count)\n",
    "    result.append(count)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The four matrix elements are shown below and now be classically processed to produce the eigenvalues."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0\n",
      "{ 0:49865 1:50135 }\n",
      "\n",
      "QKSD Matrix Element:  -0.0027 + - 0.0022360598270171573\n",
      "1\n",
      "{ 0:49796 1:50204 }\n",
      "\n",
      "QKSD Matrix Element:  -0.00408 + - 0.002236049366181346\n",
      "2\n",
      "{ 0:49695 1:50305 }\n",
      "\n",
      "QKSD Matrix Element:  -0.0061 + - 0.002236026375068058\n",
      "3\n",
      "{ 0:49972 1:50028 }\n",
      "\n",
      "QKSD Matrix Element:  -0.00056 + - 0.002236067626884303\n"
     ]
    }
   ],
   "source": [
    "mean_val = np.zeros(len(angle))\n",
    "i = 0\n",
    "for count in result:\n",
    "    print(i)\n",
    "    i_result = count.get()\n",
    "    print(i_result)\n",
    "    mean_val[i] = (i_result['0'] - i_result['1']) / shots\n",
    "    error = np.sqrt(2 * i_result['0'] * i_result['1'] / shots) / shots\n",
    "    print('QKSD Matrix Element: ', mean_val[i], '+ -', error)\n",
    "    i += 1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Classically Diagonalize the Subspace Matrix\n",
    "\n",
    "For a problem of this size, numpy can be used to diagonalize the subspace and produce the eigenvalues and eigenvectors in the basis of non-orthogonal states. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[-0.0027  -0.00408]\n",
      " [-0.0061  -0.00056]]\n",
      "Eigenvalues: \n",
      "[-0.00782313  0.00456313]\n",
      "Eigenvector: \n",
      "[[-0.76575845  0.64312829]\n",
      " [-0.64312829 -0.76575845]]\n"
     ]
    }
   ],
   "source": [
    "import numpy as np\n",
    "\n",
    "my_mat = np.zeros((2, 2), dtype=float)\n",
    "m = 0\n",
    "for k in range(2):\n",
    "    for j in range(2):\n",
    "        my_mat[k, j] = mean_val[m]\n",
    "        m += 1\n",
    "\n",
    "print(my_mat)\n",
    "\n",
    "E, V = np.linalg.eigh(my_mat)\n",
    "\n",
    "print('Eigenvalues: ')\n",
    "print(E)\n",
    "\n",
    "print('Eigenvector: ')\n",
    "print(V)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}