{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Hadamard Test\n",
    "\n",
    "Consider the observable $O$ and two generic quantum states $\\ket{\\psi}$ and $\\ket{\\phi}$. We want to calculate the quantity\n",
    "$$\n",
    "\\bra{\\psi} O \\ket{\\phi}.\n",
    "$$\n",
    "where $O$ is a Pauli operator.\n",
    "\n",
    "First of all we shall prepare the states $\\ket{\\psi}$ and $\\ket{\\phi}$ using a quantum circuit for each of them. So we  have\n",
    "$$\n",
    "\\ket{\\psi} = U_{\\psi}\\ket{0} \\qquad \\ket{\\phi} = U_{\\phi}\\ket{0}\n",
    "$$\n",
    "\n",
    "Let's define an observable we want to use:\n",
    "$$\n",
    "O = X_1X_2\n",
    "$$\n",
    "\n",
    "Now we can evaluate the matrix element using the following fact:\n",
    "$$\n",
    "\\bra{\\psi}O\\ket{\\phi} = \\bra{0}U_\\psi^\\dagger O U_\\phi\\ket{0}\n",
    "$$\n",
    "This is just an expectation value which can be solved with a simple Hadamard test. The probability to measure $0$ or $1$ in the ancilla qubit is\n",
    "\n",
    "$$\n",
    "P(0) = \\frac{1}{2} \\left[ I + Re \\bra{\\psi} O \\ket{\\phi} \\right]\n",
    "$$\n",
    "\n",
    "$$\n",
    "P(1) = \\frac{1}{2} \\left[ I - Re \\bra{\\psi} O \\ket{\\phi} \\right]\n",
    "$$\n",
    "\n",
    "The difference between the probability of $0$ and $1$ gives \n",
    "\n",
    "$$\n",
    "P(0)-P(1) = Re \\bra{\\psi} O \\ket{\\phi}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### A- Numerical result as a reference: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Psi state:  (0.707107,0)\n",
      "(0,0)\n",
      "(0.707107,0)\n",
      "(0,0)\n",
      "\n",
      "Phi state:  (0,0)\n",
      "(1,0)\n",
      "(0,0)\n",
      "(0,0)\n",
      "\n",
      "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": [
    "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",
    "@cudaq.kernel\n",
    "def psi(num:int):\n",
    "    q = cudaq.qvector(num)\n",
    "    h(q[1])\n",
    "    \n",
    "@cudaq.kernel\n",
    "def phi(n:int):\n",
    "    q = cudaq.qvector(n)\n",
    "    x(q[0])\n",
    "\n",
    "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)\n",
    "\n",
    "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": [
    "### B- Using ``sample`` algorithmic primitive to sample the ancilla qubit and compute the expectation value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{ 0:85356 1:14644 }\n",
      "\n",
      "Observable QC:  0.70712 + - 0.0015811092713661505\n",
      "Numerical result 0.7071067690849304\n"
     ]
    }
   ],
   "source": [
    "import cudaq\n",
    "\n",
    "cudaq.set_target('nvidia')\n",
    "\n",
    "@cudaq.kernel\n",
    "def U_psi(q:cudaq.qview):\n",
    "    h(q[1])\n",
    "\n",
    "@cudaq.kernel\n",
    "def U_phi(q:cudaq.qview):\n",
    "    x(q[0])\n",
    "\n",
    "@cudaq.kernel  \n",
    "def ham_cir(q:cudaq.qview):\n",
    "    x(q[0])\n",
    "    x(q[1])\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)\n",
    "\n",
    "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))\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### C- Use multi-GPUs to compute the matrix elements"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Number of QPUs: 5\n",
      "0\n",
      "{ 0:63807 1:36193 }\n",
      "\n",
      "Observable QC:  0.27614 + - 0.0021491238917289066\n",
      "1\n",
      "{ 0:49929 1:50071 }\n",
      "\n",
      "Observable QC:  -0.00142 + - 0.0022360657230949183\n",
      "2\n",
      "{ 0:50041 1:49959 }\n",
      "\n",
      "Observable QC:  0.00082 + - 0.0022360672257336093\n",
      "3\n",
      "{ 0:50276 1:49724 }\n",
      "\n",
      "Observable QC:  0.00552 + - 0.0022360339102974265\n"
     ]
    }
   ],
   "source": [
    "import cudaq\n",
    "\n",
    "cudaq.set_target(\"nvidia-mqpu\")\n",
    "\n",
    "target = cudaq.get_target()\n",
    "qpu_count = target.num_qpus()\n",
    "print(\"Number of QPUs:\", qpu_count)\n",
    "\n",
    "@cudaq.kernel\n",
    "def U_psi(q:cudaq.qview, theta:float):\n",
    "    ry(theta, q[1])\n",
    "\n",
    "@cudaq.kernel\n",
    "def U_phi(q:cudaq.qview, theta: float):\n",
    "    rx(theta, q[0])\n",
    "\n",
    "@cudaq.kernel  \n",
    "def ham_cir(q:cudaq.qview):\n",
    "    x(q[0])\n",
    "    x(q[1])\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)\n",
    "    \n",
    "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, qubit_num, angle[i], theta[i], shots_count = shots, qpu_id = i%qpu_count)  \n",
    "    result.append(count)  \n",
    "\n",
    "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('Observable QC: ',  mean_val[i],'+ -', error)\n",
    "    i += 1\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Diagonalize the matrix using for example Numpy or CuPy. In this example, since we are having 2x2 matrix, we use numpy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 0.27614 -0.00142]\n",
      " [ 0.00082  0.00552]]\n",
      "Eigen values: \n",
      "[0.00551752 0.27614248]\n",
      "Eigenvector: \n",
      "[[ 0.00303004 -0.99999541]\n",
      " [-0.99999541 -0.00303004]]\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('Eigen values: ')\n",
    "print(E)\n",
    "\n",
    "print('Eigenvector: ')\n",
    "print(V)"
   ]
  }
 ],
 "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.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}