vqe_rcs_qiskit.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3061cde8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Installation of the requirements\n",
    "#!python -m pip install -r requirements.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "bf3f9f51",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SparsePauliOp(['IIII', 'IIIZ', 'IIZZ', 'IIZI', 'IZII', 'IZIZ', 'ZXIX', 'IXZX', 'ZXZX', 'IXIX', 'ZZZI', 'ZZZZ', 'IZZZ', 'ZZIZ', 'ZIZI'],\n",
      "              coeffs=[-0.81054798+0.j,  0.17218393+0.j, -0.22575349+0.j,  0.12091263+0.j,\n",
      "  0.17218393+0.j,  0.16892754+0.j,  0.0452328 +0.j, -0.0452328 +0.j,\n",
      " -0.0452328 +0.j,  0.0452328 +0.j, -0.22575349+0.j,  0.16614543+0.j,\n",
      "  0.16614543+0.j,  0.17464343+0.j,  0.12091263+0.j])\n"
     ]
    }
   ],
   "source": [
    "'''\n",
    "(C) Renata Wong 2023\n",
    "\n",
    "Qiskit code for testing fidelity of randomised classical shadow on the ground state energy of molecules.\n",
    "\n",
    "Procedure:\n",
    "1. Choose a variational ansatz with initial parameters selected at random.\n",
    "2. Generate a set of random basis change operators.\n",
    "3. Apply the random operators to change bases in the ansatz.\n",
    "4. Measure the ansatz in the Pauli Z basis and store the results as a shadow.\n",
    "5. Obtain the expectation value of the molecular Hamiltonian from the shadow.\n",
    "6. Optimize for minimum Hamiltonian expectation value. \n",
    "7. Feed the calculated angles/parameters back to the ansatz.\n",
    "8. Repeat steps 3-7 till the optimization is completed. \n",
    "9. Output the minimized expectation value of the molecular Hamiltonian and the mean-square-root-error. \n",
    "\n",
    "Note: Below we perform calculations on the molecular Hamiltonian of H_2.\n",
    "To perform calculations on other molecules, you will need to specify their geometry, charge and spin \n",
    "to replace the values in the driver. \n",
    "\n",
    "Note: predicting_quantum_properties module comes from https://github.com/hsinyuan-huang/predicting-quantum-properties\n",
    "'''\n",
    "\n",
    "from qiskit import QuantumCircuit, execute\n",
    "from qiskit_aer import QasmSimulator\n",
    "\n",
    "from qiskit_nature.units import DistanceUnit\n",
    "from qiskit_nature.second_q.drivers import PySCFDriver\n",
    "from qiskit_nature.second_q.mappers import BravyiKitaevMapper\n",
    "\n",
    "from qiskit.algorithms.minimum_eigensolvers import NumPyMinimumEigensolver\n",
    "from qiskit_nature.second_q.algorithms import GroundStateEigensolver\n",
    "\n",
    "from qiskit.circuit.library import EfficientSU2\n",
    "from qiskit.algorithms.optimizers import SLSQP, COBYLA, SPSA\n",
    "\n",
    "import time\n",
    "import numpy as np\n",
    "from functools import partial\n",
    "import matplotlib.pyplot as plt\n",
    "from collections import Counter\n",
    "\n",
    "from predicting_quantum_properties.data_acquisition_shadow import randomized_classical_shadow\n",
    "from predicting_quantum_properties.prediction_shadow import estimate_exp\n",
    "\n",
    "# handle deprecation issues\n",
    "import qiskit_nature\n",
    "qiskit_nature.settings.use_pauli_sum_op = False\n",
    "import h5py\n",
    "H5PY_DEFAULT_READONLY=1\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "\n",
    "# classically obtained ground state energy\n",
    "EXPECTED_EIGENVALUE = -1.86\n",
    "\n",
    "# specifying the geometry of the molecule in question\n",
    "driver = PySCFDriver(\n",
    "    atom=\"H 0 0 0; H 0 0 0.735\",\n",
    "    basis=\"sto3g\",\n",
    "    charge=0,\n",
    "    spin=0,\n",
    "    unit=DistanceUnit.ANGSTROM,\n",
    ")\n",
    "\n",
    "\n",
    "problem = driver.run()\n",
    "hamiltonian = problem.hamiltonian\n",
    "\n",
    "# electronic Hamiltonian of the system\n",
    "second_q_op = hamiltonian.second_q_op()\n",
    "\n",
    "# The Bravyi-Kitaev repserentation of the Fermionic Hamiltonian\n",
    "mapper = BravyiKitaevMapper()\n",
    "bkencoded_hamiltonian = mapper.map(second_q_op)\n",
    "\n",
    "print(bkencoded_hamiltonian)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "9b64bd5c",
   "metadata": {},
   "outputs": [],
   "source": [
    "def process_hamiltonian(hamiltonian, derandomize = False):\n",
    "\n",
    "    hamiltonian_observables = []\n",
    "    hamiltonian_coefficients = []\n",
    "\n",
    "    for observable in hamiltonian.paulis:\n",
    "        op_list = []\n",
    "        for op_index, pauli_op in enumerate(observable):\n",
    "            pauli_op = str(pauli_op)\n",
    "            if pauli_op == 'I' or pauli_op == 'X' or pauli_op == 'Y' or pauli_op == 'Z':\n",
    "                op_list.append((pauli_op, op_index))\n",
    "        hamiltonian_observables.append(op_list)\n",
    "    \n",
    "    hamiltonian_coefficients = hamiltonian.coeffs.real\n",
    "    \n",
    "    system_size = len(hamiltonian_observables[0])\n",
    "\n",
    "    # removing all occurrences of Pauli-I, for all-Pauli-I there is an empty list left\n",
    "    # these observables are needed for estimate_exp()\n",
    "    observables_xyze = []\n",
    "\n",
    "    for observable in hamiltonian_observables:\n",
    "        XYZE = []\n",
    "        for pauli in observable:\n",
    "            if pauli[0] != 'I':\n",
    "                XYZE.append(pauli)\n",
    "        observables_xyze.append(XYZE)   \n",
    "    \n",
    "    # derandomisation procedure requires that coefficients are non-negative\n",
    "    if derandomize == True:\n",
    "        absolute_coefficients = [abs(coeffcient) for coeffcient in hamiltonian_coefficients]\n",
    "\n",
    "        # removing the empty list as well\n",
    "        # these observables are needed for derandomisation procedure\n",
    "        observables_xyz = []\n",
    "        for idx, observable in enumerate(observables_xyze):\n",
    "            if observable:\n",
    "                observables_xyz.append(observable)\n",
    "            else:\n",
    "                absolute_coefficients.pop(idx)\n",
    "        \n",
    "        return system_size, observables_xyze, observables_xyz, hamiltonian_coefficients, absolute_coefficients\n",
    "    \n",
    "    return system_size, observables_xyze, hamiltonian_coefficients\n",
    "\n",
    "\n",
    "\n",
    "# process the Hamiltonian to obtain properly formatted data\n",
    "hamiltonian_data = process_hamiltonian(bkencoded_hamiltonian, derandomize = False)\n",
    "system_size, observables_xyze, hamiltonian_coefficients = hamiltonian_data\n",
    "\n",
    "#print('HAMILTONIAN\\n', observables_xyze)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "1d747421",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAbsAAADWCAYAAABIdCwuAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAvXElEQVR4nO3deVxUZf//8dfMALILiEqCC8hSopDiiguYlubXr+JtWmreaX7TwF+5ttyhlWHaz6y4v7/b7G4x7rtU7iRLM80WAffKShM3VFxSMfcFxQWG3x8TKAIzA8zMmXP6PB+PeSRnzlzzhs65PnOuuc45urKysjKEEEIIDdMrHUAIIYSwNyl2QgghNE+KnRBCCM2TYieEEELzpNgJIYTQPCl2QgghNE+KnRBCCM2TYieEEELzpNgJIYTQPCl2QgghNE+KnRBCCM2TYieEEELzpNgJIYTQPCl2QgghNE+KnRBCCM2TYieEEELzpNgJIYTQPCl2QgghNE+KnRBCCM2TYieEEELzpNgJIYTQPCl2QgghNE+KnRBCCM2TYieEEELzpNgJIYTQPCl2QgghNM9F6QDObt86uHxKmff2aQJR99XttUrlrk9m4TiyXdeObNfqJ8XOgsun4MIxpVPUnlpzC8dQ6/ah1txCeTKMKYQQQvOk2AkhhNA8GcYUqnf1Oly+DjrA1wPcXZVOJIRwNlLshOqUlcHhM7B5Pxw8Beeu3HpOBwT6QGQQdI+AZv6KxRRCOBEpdjYwbWEie45swWBwRa83EOQfysg+qSTEDlM6mllqzP37RcjcCofOVP98GXD6sumxaT9EB8OwzuDn6dCYmqDG7QPUm1vYlxQ7GxnVdyaj+s6gtLSEFZv/wdwlIwkPbk9wYLjS0cxSU+6tByDrRygxWv+aXcehYBWM7g5tgu2XTavUtH3cTq25hf3IBBUbMxhceLDLE5QaSzh4YrvScazm7Lk35UPm97UrdOWKb8L7uZAnU9brzNm3j5qoNbewPSl2Nnaz5AarNi8EICQwUuE01nPm3IfPmI7ozEkfZXrUxFgGH22Cs0W2zfZn4czbhzlqzS1sT4qdjSz57lWSZvox8AUPPlw7g6nD3iesWQwAx88cICU9jpslNwD4JOd1Mta+qGTcCuZyz1k8kq27V1Ws+1JGEtv2fe3QfCWlsGSL6bu4+rpeYvq+r8wWjdmIM2WpjmzXQiucutgZjUbmz59PREQE7u7uxMbGkpubS1RUFOPHj1c6XiUj+6TyedoFsl4+Q+e7B7DjQHbFc8GB4fRoN5TMdXMpPHeInO2ZjOyTqmDaW8zlTh6cTsbamRRfL2LDzuV4uTekY9QDDs338xE4dcl27e3/HQpO2669uvj9InzyAzz/CUxZAi8thzW/QtE1ZXNVR7ZroRVOXezGjRtHWloaEyZMYM2aNQwfPpwRI0ZQUFBAXFyc0vGq5ePpz9Rh7/P93i/ZnLeiYvnwxGfYumcVcxaPIHlQOm4uDRRMWVV1uf29mzCkxyQWrHiaJd/N5slBbzk816Z8dbRprT0n4PXVsGU/XLtpWnaxGNbuNC0/c1m5bObIdi3UzmmL3dKlS8nIyGDlypVMnz6d3r17k5qaSrdu3SgpKaFDhw5KR6yRr2cAQ3tOZdFXL2A0mmZUuBhcaRfWi6Li87QN7aFwwupVl7tfpzEcO51PUven8fUMcGieq9fhyFnbt7u3UJnhw6JrsGg9lBqrH5a9dA0+WO+8Q5uyXQs1c9piN2fOHPr3709CQkKl5eHh4bi6uhITYxp/P3z4MAkJCURGRtKuXTs2bNigRNwqhvScxLlLhXzz078BOHxyF7sOb6J9eF9Wf/+ewulqdmdugGaNwhWZsn3svH3avXqj8onojrL1INwsrfn7x7IyKLyg/DCrObJdC7XSlZU53+fIY8eO0bx5cz744AMef/zxSs+NGDGCvXv38ssvvwDQr18/Bg8eTEpKCps3b2bYsGEcOnQINzc3s++h0+msyjL/yWxiWyfW6fcoZzQamfZOAsmD0gkJjGTSgnjmjf8Wf5+mZl+342AO09/pXaf3tEXucvMyxzCgy/9Y9cm9Ppnv1KbXWO4fv6jSMnMzLs2ZvLjyz5++2ptje3Lq1lgdDU3NJjiqFzp9zZ8xy8qMbPviNTZ/Yt/vvmS7Vm67FrZlbQlzyiO7Y8dMJ0QFBQVVWl5cXExubm7FEOaZM2fYuHEj48aNAyA+Pp5mzZqRnZ2NM/liy0IiguOIDInD092HMf3SeHvlZKVjOT2d3qDKtmticGlgttABUFaG3mD+g5qzkO1aqIlTHtkdOHCAiIgI3nrrLSZPnlyxfNasWbz88sssWLCAlJQUfv75Zx566CEKCgoq1hk+fDh9+/a12WzNbZnK3T/LLwQ6PlK31yqVuz6Z77TjKHxo5ah0+RHfnUdwNZn2IDR38Fc1mVvh+4OWT6MY0RW6tLZvFtmua8eW27VQhlNeLiwsLIyYmBjmzJlDQEAAwcHBZGVlsXr1agCnnYkpbCvYThdxNujhrob2aduc7hGm7+3MaeAC7Vs6Jo8QfyZOOYyp1+tZtmwZ0dHRJCcnM3bsWAIDA5k4cSIGg6FickqLFi34/fffuX79esVrDx06RMuW0ltoQSNv8HW3fbvNA8DF8aOYNG8EPWu4iEf5N8jDOoObU34EFULdnLLYAURGRpKdnc2VK1c4evQoaWlp7Ny5kzZt2uDh4QFAYGAg3bt354MPPgBg8+bNHD9+nN695YtkLdDpoKsdJsvZo01rDekIA+8Fzzu+lmvkDY/3go6hisQSQvNU9Rly27ZtdO3atdKyd955hzFjxpCeno6bmxtLly61OBNTqEd8BKzbY7psmC34uEMHBQ/89TroGw2Jd8P0TNOyp++H0Mam4i6EsA/VFLuioiLy8/NJSUmptDwsLIz169crlMrkg9V/Y9fhTUS36k5I4ygys+cyeei7xLZO4JOc19m8awVN/VvyzMMZ3Cy5zrPv9iW4UTjPj/zYaXMDbNi5nIUrJrFkxm8UXy9SJLefJwyMhc9/tk17zjJMePswalgT5XKYU9P20S60J/My/8qpC0dxMbiR+mgmbi7uim7XNWWNbhXP1Ld7cejkTt6Zsp3gwHAKzx1iXuZf0aEjsGEIz434CIPewIxFAykqvkD6xI0Ozy/sz2mHMe/k7e1NaWkpTz31lNJRKjl0Mo8r1y7xZsp6Ll09y7UbVxiW8AyxrRM4X3SK7QezSZ+4kdC7YtiU9zkeDbxJHZWpdGyzuctt+DWLxn7NARTN3SsK7mlmfp3Jiy3PxIyPgJjmtsulZea2j4MntuPi4sabKevp12ks3/28WNHtw1xWg96FWWM+p2e7hyrW93b3Y/bYVbyZsp6ggFB+2Gua+Db78VU1vYXQANUUO2eVd2gjHSNNF5HtEHE/+tvO38r/bRuxYYl/PNeXPUe2KBGxWuZyA3y/ZzUdIvqi0ym/iej1MLYn3H1X3dvoHAYPdbRdJq0zt30ENgzGaDSNKxcVX8DXs5EiGcuZy6rT6aqc5O7j6Y+Xh2k6rovBFb1OgdlKwuGcYEBH3S5fPceqLe/w6Ya3KCq+QELscPy8TeNSV65dwNPdFwAv94YUXbugYNLKzOUG+Oanf/HcIx/x1Y+LzLTiOG4u8EQirNsNX+00XV/SGu6ukNTBdN6afCdmPXPbh69XINdvFvP46/dg0Bn4f0//4LRZzTlz8QQ/5X/DqD4zHJBSKE2KXT35eAbwWL9XiI8exNbdqzh98dYZr17uDTn9xxmwV69dwtvdT6GUVZnL/cuBdbRp2Q1XF+ea6GPQw/1tIbYFbNgHPxSY7lFXHe8GplmXPSJN3/uJ2jG3ffyU/zUNvRqz6Jk9rP81i2W58xl9v3L3sTOXtSY3Sq7z+n8eY+qw9zAYpBv8M1B+jErl2ob2YGeBaYLMjoM5FcM7AJHNO/FrQS4AP+//lntadq22DSWYy334ZB5bdq3kb+/158jvu/jwK+f65NvEF4Z2grShMLkfDL1tePLhLqaro7w8xDTFXwpd3ZjbPsrKyiruFNDQK5Ar1y4qkrGcuaw1Sc8az6D4ibRs2sbe8YSTkGJXT6FBbXExuDJtYSIuBlfc3bwqnvP3bkK7sF5MXtCDgye2Ex+dpFzQO5jLPaTH07z+5DrmPvEVLZtGM7b/bAWT1szNBVoFQs+oW8u6hSt30riWmNs+OkY+wNFTe5i2MJF/rX2RQd1SzLRkf+ayAqR9NJyf9n/NvMzH2Jy3gt2Ht7AxbznLN6QzbWEiG3d+plBy4Uhy/G4D4wbMrfj3+l+zyMx+jeDACGJbJ/BI7+d4pPdzFc8XXy/itaWPEtW8kxJRKzGXu1z5NGxnyi0cw9z28eJfsyqtq/T2YS7rzNGfVFl/5eyqd8mdsWggAb71mAUlnJpTXgjamcgFc2tHqQvmlp92UNdbAClBycyyXdeOXAha/eTIzgIfBU/4rc97K5Vbyb+XsJ5s1+p4X2E7UuwsiLpP6QR1o9bcwjHUun2oNbdQnkxQEUIIoXlS7IQQQmieFDshhBCaJ8VOCCGE5kmxE0IIoXlS7IQQQmieFDshhBCaJ8VOCCGE5kmxE0IIoXlS7IQQQmieFDshhBCaJ8VOCCGE5kmxE0IIoXly1wML9q2Dy6eUeW+fJnW/yrtSueuTWQhL1Lhdq7UP0RopdhZcPqXcTS7rQ625hTBHjdu1GjNrkQxjCiGE0DwpdkIIITRPhjGFcKBrN2FvIfx29tayjA3QzA9aBEJEUzDIR1AhbE6KnQ1MW5jIniNbMBhc0esNBPmHMrJPKgmxw5SOZpZac6vRhavwdR5sOwQ3Sio/t/2o6QHQ0APiI6D3PeAme2edqHG7VmNmtZHdyUZG9Z3JqL4zKC0tYcXmfzB3yUjCg9sTHBiudDSz1JpbTX4ogOXbTEd1llwshjW/moriyG4Q2tj++bRIjdu1GjOriQyY2JjB4MKDXZ6g1FjCwRPblY5jNbXmdnard8CSLdYVutudvgz/+BZ2/mafXH8Watyu1ZhZDaTY2djNkhus2rwQgJDASIXTWE+tuZ3Z+n2mocu6KjVCxkYoUOgcLS1Q43atxsxqIMXORpZ89ypJM/0Y+IIHH66dwdRh7xPWLAaA42cOkJIex82SGwB8kvM6GWtfVDJuBXO55yweydbdqyrWfSkjiW37vlYqqqr8fhFW/mx+nfRRpoc5pUZYsrXq93zCPDXuj7Iv2pdTFzuj0cj8+fOJiIjA3d2d2NhYcnNziYqKYvz48UrHq2Rkn1Q+T7tA1stn6Hz3AHYcyK54LjgwnB7thpK5bi6F5w6Rsz2TkX1SFUx7i7ncyYPTyVg7k+LrRWzYuRwv94Z0jHpAwbTqsfwnKDHapq0zl2HdHtu0VR83S03fJX7xi2l4tuAUlJUpnap6atwfZV+0L6cuduPGjSMtLY0JEyawZs0ahg8fzogRIygoKCAuLk7peNXy8fRn6rD3+X7vl2zOW1GxfHjiM2zds4o5i0eQPCgdN5cGCqasqrrc/t5NGNJjEgtWPM2S72bz5KC3FE6pDr9fgn2Ftm1z837TUZ5Sth+BFz+FjzfDd7tNw7P/+w3MX2Mqxs5Kjfuj7Iv24bTFbunSpWRkZLBy5UqmT59O7969SU1NpVu3bpSUlNChQwelI9bI1zOAoT2nsuirFzAaTT2Ui8GVdmG9KCo+T9vQHgonrF51uft1GsOx0/kkdX8aX88AhROqw0+HbN/mpWLIP2n7dq2x65jpu8PiaibZnLgA/+8buHzN4bGspsb9UfZF23PaYjdnzhz69+9PQkJCpeXh4eG4uroSE2May37xxReJjIxEr9eTlZWlRNRqDek5iXOXCvnmp38DcPjkLnYd3kT78L6s/v49hdPV7M7cAM0ahcv051o4ctbyOs7UrjllZbDiF9CZef5iMWzMd2isWlPj/ij7om055Xl2x44dIy8vjylTplR57ujRo0RHR9OggWnYoX///owZM4bHH3/c0TErvJGcU2WZl7svy185B5i+e/z78id5asgCQgIjmbQgnvjowfj7NHVw0sos5RZ1U3jBTu2et0+75hw5C6cuWV5v8354MMb+eayhxv1R9kX7c8oju2PHTJcIDwoKqrS8uLiY3NzcSkOY8fHxhIWF1fo9dDqdVY/c3Jx6/S4AX2xZSERwHJEhcXi6+zCmXxpvr5xs8XW5uTlW57RH7rqoT+b6PMop8d53Pk6fq1wdymdd3vmw9vlyK778yuG/y/0DHrLq//vla6DT6+2eR6n9UY37Yn1zq+VhLac8sgsMDAQgPz+fAQMGVCyfN28ehYWFTjs5pSaDu0+s9HP3tkl0b5ukTJg6ePaRDKUjqIrxjyntamnXnBvFF61a7+b1K847NfMOat4fZV+sO6csdmFhYcTExDBnzhwCAgIIDg4mKyuL1atXA9ik2JVZuWNuy1TuXlQJCYmULaxbB6JU7vpkro/Ji03/tfb/qz2lr4XDZ279XJ6tynqjzD9/p3GjBrHyDcf+fiWl8NJncOV6zevogJ5tvFjggL+9GrdrtfYhWuOUw5h6vZ5ly5YRHR1NcnIyY8eOJTAwkIkTJ2IwGCompwjhjJrbaaKcvdo1x8UAve+p+XkdoNdDwt0OiyREnTjlkR1AZGQk2dnZlZaNHj2aNm3a4OHhoVAqISyLaQEbbDw70dUAdzezbZvWuq+N6a4NG/NNxe324wSDHsb2hLv8lMkmhLWc8siuJtu2basyhDlz5kxCQkLYsmULEyZMICQkhIMHDyqUUAgIbwJNfW3bZlwr8HSzbZvW0uvgoU4wtT90um0u2IAYeDEJokOUySVEbTjtkd2dioqKyM/PJyUlpdLytLQ00tLSFEpl8sHqv7Hr8CaiW3UnpHEUmdlzmTz0XaJbxTP17V4cOrmTd6ZsJzgwnOLrRTz7bl+CG4Xz/MiPnTJ3bOsEBs9sSHiz9gC89NhyfD0DmLFoIEXFF0ifuFHR3M5Op4PBHeDdHNu018AF+rWzTVv10aKR6bZDPxSYfn7ACTJVx9x2/Y/Pn+JQ4U7uahTGlIfew6A3KLpd16bvuHjlDC9+OAiDwRUv94bMePQ/GI2lTtOfODvVHNl5e3tTWlrKU089pXSUSg6dzOPKtUu8mbKeS1fPcu3GFYYlPENs6wQMehdmjfmcnu1uTd/2aOBN6qhMBRObmMsNEBrUjjeSc3gjOafiag2zH19lrklxmzbB0Ln2Z8RUKykO/L1s05bWmduu9/32IyUlN3gjOYeWTaP5/o8LKyu1Xde27/D28OetlI28mZxLZHAcW3evcpr+RA1UU+ycVd6hjXSMNF2QtUPE/ej1horndDqd4ieO18RcboCjp/Yw5e2evL/6eaeY4ahGwzpDZFDNz09ebHkm5n1toGtr2+bSMnPbdeHZAkLvMk1ua93sXnYd2axIxnK17TsMegN6vanLLi0rJTgwwnFhNUCKXT1dvnqOf339EtMWJrLku1e5fFUdVzywlDvjuf28mbyeoqvn2bL7C4VSqpurAZ5IrNsRnkFvGgr973tNw6LCOua265DGUfxakAvA9gPruFJ8QaGUJnXpO/Ye/YGUv3dk+4F13BUQ6oCU2qGa7+yclY9nAI/1e4X46EFs3b2K0xcVOqGmlizlLh+6jG+bxIHjvxAfPUiJmKrnajB9zxXbAlb8bN2ltyKawl86ygzHujC3XYcH30uroLZMf6c3rYLa4qfwqEtd+o67W3Tm7UnbWJb7Bl/9uIihvapeUlFUT47s6qltaA92FqwHYMfBHIzGUoUTWcdc7uIbVyj94+ddhzfRrJGMo9VXdDD8bSBM7AM9IqBlI/BwMxVD7wam4c4+beD5gTCxrxS6urK0P46+/0XmP5mNr2cjutzzX0pErFDbvuPmbVfQ8XL3xc1VTsGqDTmyq6fQoLa4GFyZtjCRNi274e7mRanx1m2l0z4aTt7hjRw/s5+HE58lvu1gBdPeYi738dP7eWPZ43i4eRMUEMpfH5ilcFpt0OkgIsj0EPZhbrs2Go0888/70OsNtA/vwz0tujhtVqjadwT43sW7Xz6DXqfHxyOA50Z8pGB69ZFiZwPjBsyt+Pf6X7PIzH6N4MAIYlsnMHP0J5XWLb5exGtLHyWqeSdHx6zCXO6Fk3+usv6MRQMJ8L3LkRGFqDVz23V1dxdQcruuTd8B8GZybqWfnak/cXa6MplqZ5aS17XzC4GOj9TttUrlrk/m+iif1Vjd3QKEbSn5t1bjdq3WPkRr5MjOAp8m6nxvpXIr+fcS2qfG7VqtfYjWSLGzIOo+pRPUjVpzC2GOGrdrNWbWIpmNKYQQQvOk2AkhhNA8KXZCCCE0T4qdEEIIzZNiJ4QQQvOk2AkhhNA8KXZCCCE0T4qdEEIIzZNiJ4QQQvOk2AkhhNA8KXZCCCE0T4qdEEIIzZNiJ4QQQvPkrgcW7FsHl08p894+Tep+xXSlctcnsxBapMY+RI2ZLZFiZ8HlU8rdeLE+1JpbCK1R476oxsyWyDCmEEIIzZNiJ4QQQvOk2AkhrFJqhBPnb/18tgjKypTLI0RtyHd2Qoga3SyFX47A1gNw9CyUGG89l7YCvBpAVBD0iITQxqDTKZdVCHOk2NnAtIWJ7DmyBYPBFb3eQJB/KCP7pJIQO0zpaGapNbdwjN3H4T/fw8Ximte5ch1+PmJ6RAbBI10gwNtxGbVCrfuimnLLMKaNjOo7ky9eLWL5y2fp12kMc5eM5PiZA0rHskituYX9GMtg+TZ4N8d8obtT/kn4v1/CLo3N4nMUte6Lasktxc7GDAYXHuzyBKXGEg6e2K50HKupNbewrbIyyPoB1u+r2+uvl8AH62HXcdvm+jNR677o7Lml2NnYzZIbrNq8EICQwEiF01hPrbmFbf14CDZb+FCePsr0qImxDD7aBBeu2jbbn4Va90Vnz+3Uxc5oNDJ//nwiIiJwd3cnNjaW3NxcoqKiGD9+vNLxKlny3askzfRj4AsefLh2BlOHvU9YsxgAjp85QEp6HDdLbgDwSc7rZKx9Ucm4FczlnrN4JFt3r6pY96WMJLbt+1qpqMLOLhXDZ9ts09a1m/DJ97Zpyxau3YTLxaYZpc5K+hD7cupiN27cONLS0pgwYQJr1qxh+PDhjBgxgoKCAuLi4pSOV8nIPql8nnaBrJfP0PnuAew4kF3xXHBgOD3aDSVz3VwKzx0iZ3smI/ukKpj2FnO5kwenk7F2JsXXi9iwczle7g3pGPWAgmmFPW3aD8U3bdfe7hNw/Lzl9expzwlY8C08/wnMXA6pWfDZT7X7LtJRpA+xL6ctdkuXLiUjI4OVK1cyffp0evfuTWpqKt26daOkpIQOHTooHbFaPp7+TB32Pt/v/ZLNeSsqlg9PfIate1YxZ/EIkgel4+bSQMGUVVWX29+7CUN6TGLBiqdZ8t1snhz0lsIphb2UGmHLftu3u8kObVordy/8MxsO3HaNx2s3TcvfXGM6T9AZSR9iH05b7ObMmUP//v1JSEiotDw8PBxXV1diYmI4f/48AwcOJDIyktjYWB544AEOHFB+FpCvZwBDe05l0VcvYDSaxk1cDK60C+tFUfF52ob2UDhh9arL3a/TGI6dziep+9P4egYonFDYy+8X4dI127e7/6Tt27TGsXOmIzio/sT3S8Xw8WbHZqoN6UNszymL3bFjx8jLy2PYsKrnahw9epTo6GgaNGiATqdj8uTJ5Ofns2PHDgYOHMjYsWMVSFzVkJ6TOHepkG9++jcAh0/uYtfhTbQP78vq799TOF3N7swN0KxROMGB4QqmEvb22zn7tHv6MhTfsE/b5mzMB3Pnt5cBh04rP8xqjvQhtuWUJ5UfO2Y6UScoKKjS8uLiYnJzc3nwwQcB8PPzo2/fvhXPx8fHM2/ePKveQ2flpR7mP5lNbOtEs+u8kZxTZZmXuy/LXzH1IEajkb8vf5KnhiwgJDCSSQviiY8ejL9PU7Pt5ubm0GlEb6ty2iN3XdQnc31M+tj08d3a/6+iss5JM+j2UFqlZeZmXJp7fvLiyj8Hh97D+RN765Gu9sa8VUDDxqEW1xs4YhLb1/6vXbOosQ+xJrM1ueuitpnLrLxmnVMe2QUGBgKQn59fafm8efMoLCyscXJKeno6SUlJ9o5Xa19sWUhEcByRIXF4uvswpl8ab6+crHQsIW5jvw8JOju2XRO93mDVejor11Oa9CH1pyuztiw6kNFopH379hQWFjJ//nyCg4PJyspi9erVHD16lK1bt9KlS5dKr5k1axZr1qxh3bp1eHp62izLtkzl7uvkFwIdH6nba5XKXZ/M9VF+NGHpaERUb1M+LPvRunXL/8Z3HsHVZPZQ8HavW666+nA9/HrM8oWq/09fCDd/cFRvauxD1JjZEqc8stPr9Sxbtozo6GiSk5MZO3YsgYGBTJw4EYPBQExMTKX1Z8+ezapVq/jqq69sWuiE+LMIsdO8AT9Pxxc6gO6R5gudTgdNfKF1E8dlEspyyu/sACIjI8nOzq60bPTo0bRp0wYPD4+KZbNmzWL16tV88803+Pn5OTilENrQzB/cXU1T821JqWIS0dR0J4aN+VWf0+nARQ+PxstdGv5MnLbYVWfbtm107dq14uddu3bx8ssv07p1axITEyuWb9++3fHhhFAxVwN0CoMNdbwmZk26KTQBT6eDoR2hsQ9k76l86bK774KB90KwvzLZhDJUU+yKiorIz88nJSWlYll0dLTVM3GEEOb1ioTN+213Sa3mAcoOE+p0kHA39IyEqUtNy15KAn8v5TIJ5aim2Hl7e1NaWqp0jGp9sPpv7Dq8iehW3QlpHEVm9lwmD32XJv4tmJf5V3ToCGwYwnMjPsKgNzBj0UCKii+QPnGjU+aObZ3AT/nfkLluLsYyIxP++w2aN47i2Xf7EtwonOdHfqxobmEfjX2hfzv4ckf92zLoYURX5xgm1N82M8HZCl1N+2B0q3imvt2LQyd38s6U7RXnqH2S8zqbd62gqX9Lnnk4g5sl1x2+X9Y2M8CGnctZuGISS2b8RvH1IkX6EqecoKImh07mceXaJd5MWc+lq2e5duMKwxKeIbZ1At7ufsweu4o3U9YTFBDKD3tXAzD78VUWWrU/c7mv3yzmy63/5LXx3/BGcg6RIXF4NPAmdVSm0rGFnd3XxjTMZ87kxZZnYg6JM30PKGpmbh806F2YNeZzerZ7qGL980Wn2H4wm/SJGwm9K4ZNeZ87fL+sbeZyG37NorFfcwDF+hIpdvWUd2gjHSNNFzbtEHF/pfN7fDz98fJoCJgu9aPXOc85PeZy7z6yBZ1OzwvvP8hrS0dTfOOKUjGFgxn08HgvaNOsbq/XYSp0PZzvDi9Ox9w+qNPpqpwwnv/bNmLDEv9Yvy97jmxxWNZytc0M8P2e1XSI6ItOp2y5kWJXT5evnuNfX7/EtIWJLPnuVS5frXrlgDMXT/BT/jcVG4kzMJf7/OXfOXe5kDn/s4bolvF8ueWfCiYVjubmAv+TAEkdwKUWn88CfeCp+03fkwnLrOk7bnfl2gU83X0B8HJvSNG1Cw5IWVltMwN889O/6NPhUQekM08139k5Kx/PAB7r9wrx0YPYunsVpy9WPhPzRsl1Xv/PY0wd9h4Gg/P8uc3l9nJvSNtWPTDoDdwbfh/LcucrmFQoQa+HxHsgtoVp0sqWg1BUw4WiQwKgRwR0aGUqlMI6lvqOO3m5N+T0H2d6X712CW93PwekrKy2mX85sI42Lbvh6uLmoIQ1kyO7emob2oOdBesB2HEwB6Ox8iSa9KzxDIqfSMumbZSIVyNzuaOad+LoqT0AHDyxnaAAy9cYFNrk7wX/dS+88hdI/W94rAcM6wTDO8OE3qaro0x/ELqGS6GrLUt9x50im3fi14JcAH7e/y33tOxqdn17qG3mwyfz2LJrJX97rz9Hft/Fh1/NcETMasnmWU+hQW1xMbgybWEibVp2w93Ni1JjCQC7D29hY95yfj9/hOUb0hnSYxI92g1ROLGJudx+3o2JCUtg6tu9aODqyd9GLVE4rVCaXmeardnYV+kk2mFuHwRI+2g4eYc3cvzMfh5OfJb4toNpF9aLyQt60MSvBX/pOdnpMw/p8TRDejwNwOQFPRjbf7bDM5eTYmcD4wbMrfj3+l+zyMx+jeDACGJbJ7By9uUq689YNJAAXwtT3hzAXO6hvaYwtNeUiueLrxfx2tJHiWreSYmoQmiSuX1w5uhPqqz/SO/neKT3cxU/K7Ff1jZzufJTrZTqS5zyQtDORK0XRJULQQtxi5Lbhxr7EDVmtkSO7CzwUfAKEPV5b6VyK/n3EsIZqbEPUWNmS6TYWRB1n9IJ6katuYXQGjXui2rMbInMxhRCCKF5UuyEEEJonhQ7IYQQmifFTgghhOZJsRNCCKF5UuyEEEJonhQ7IYQQmifFTgghhOZJsRNCCKF5UuyEEEJonhQ7IYQQmifFTgghhOZJsRNCCKF5ctcDC/atg8unlHlvnyZ1v/q4Urnrk1kILVJjH6LGzJZIsbPg8inlbmJYH2rNLYTWqHFfVGNmS2QYUwghhOZJsRNCCKF5MowphNCkU5dg13H47eytZf/4Fpr5QatAaBsCbtID/mnI/2obmLYwkT1HtmAwuKLXGwjyD2Vkn1QSYocpHc0steYWwpzDZ2DNDth3supzB343PdbvAw836BYOD7QFd1fH57ydWvdFNeWWYmcjo/rOZFTfGZSWlrBi8z+Yu2Qk4cHtCQ4MVzqaWWrNLcSdSo3w5Q7I3g1lVqxffAPW7YZfjsCobhDe1O4RzVLrvqiW3PKdnY0ZDC482OUJSo0lHDyxXek4VlNrbiHAVOg+3GAqXtYUutudvwILv4Odv9klWq2pdV909txS7GzsZskNVm1eCEBIYKTCaayn1txCAGT9CHn1mCpfWgYZG+HoWcvr2pta90Vnzy3FzkaWfPcqSTP9GPiCBx+uncHUYe8T1iwGgONnDpCSHsfNkhsAfJLzOhlrX1QybgVzuecsHsnW3asq1n0pI4lt+75WKqoQ1dp9HLYcML9O+ijTw5xSIyzZAjdLbZetNqQPsS+nLnZGo5H58+cTERGBu7s7sbGx5ObmEhUVxfjx45WOV8nIPql8nnaBrJfP0PnuAew4kF3xXHBgOD3aDSVz3VwKzx0iZ3smI/ukKpj2FnO5kwenk7F2JsXXi9iwczle7g3pGPWAgmmrKr4BG/bd+vmHArhRolwe4VjGMli+zXbtnbwIm/Jt115tSB9iX05d7MaNG0daWhoTJkxgzZo1DB8+nBEjRlBQUEBcXJzS8arl4+nP1GHv8/3eL9mct6Ji+fDEZ9i6ZxVzFo8geVA6bi4NFExZVXW5/b2bMKTHJBaseJol383myUFvKZyysp8Pw4vL4dPbOrslW+Dlz2D/74rFEg60rxDOFNm2zU37TUVUKdKH2IfTFrulS5eSkZHBypUrmT59Or179yY1NZVu3bpRUlJChw4dlI5YI1/PAIb2nMqir17AaDQC4GJwpV1YL4qKz9M2tIfCCatXXe5+ncZw7HQ+Sd2fxtczQOGEt+wthI82VT/kVHwT/pkNJ847PpdwrJ8O277N05crn5unBOlDbM9pi92cOXPo378/CQkJlZaHh4fj6upKTIxpTDgpKYmYmBjat29P586d+fbbb5WIW8WQnpM4d6mQb376NwCHT+5i1+FNtA/vy+rv31M4Xc3uzA3QrFG4000jXrOj5ufKyqC0FL7b7bg8Qhn2mlDy2zn7tFsb0ofYllOeZ3fs2DHy8vKYMmVKleeOHj1KdHQ0DRqYDuEzMjLw8/MD4JdffiExMZFz585hMBgclveN5Jwqy7zcfVn+immPMRqN/H35kzw1ZAEhgZFMWhBPfPRg/H2UPbHHUm5ndfoyHLHQyZUB24/Aw13kKhlaVWqE05fs03bhBfu0WxPpQ+zPKY/sjh0zzSEOCgqqtLy4uJjc3NxKQ5jlhQ7g4sWL6HQ6ysosD7jrdDqrHrm5OfX+fb7YspCI4DgiQ+LwdPdhTL803l452eLrcnNzrM5pj9x1UZ/M1j5i47pblaW0DAKaBNs9jzyUebh7elU5p6581uWdD2ufL/feBx/aLKca+xCl+o+6ZLaWU37mDQwMBCA/P58BAwZULJ83bx6FhYVVJqdMnDiRNWvWcPHiRT799FNcXJzr1xrcfWKln7u3TaJ72yRlwtTBs49kKB2hkuKiM1atZzSWcP3qBfuGEYox/jENv6ysrFadnjVKblyzaXv1JX1I/enKrDkMcjCj0Uj79u0pLCxk/vz5BAcHk5WVxerVqzl69Chbt26lS5cuVV6Xm5vLlClTWL9+Pd7e3jbJsi1Tufs6+YVAx0fq9lqlctcns7XKymD+GtMElJo2Xh0Q0xzG9rJvFqGsV1eahrUtKT96m7zYunb/0hF6RdU91+3U2IeoMbMlTjmMqdfrWbZsGdHR0SQnJzN27FgCAwOZOHEiBoOhYnLKnRISEtDr9WzatMnBiYUj6XTQv535QqfTQZ9oR6YSSmhup8l99mpXKMe5xvtuExkZSXZ2dqVlo0ePpk2bNnh4eABQVFTE2bNnadmyJWCaoHLw4EHuueceh+cVjtWuuWnyybIfTOdElY9ilZWBqwEe6wktGimbUdhf+5bw8xHbtunvBS1l29Ecpy121dm2bRtdu3at+PnKlSs8/PDDFBUV4eLigru7Ox9//DEtWrRQMKVwlG7hEB0MWw+azovS6aB1E+gUBp5uSqcTjtAmGPw94fxV27XZPQL0TjnmJepDNcWuqKiI/Px8UlJSKpY1bdqUrVu3KpjK5IPVf2PX4U1Et+pOSOMoMrPnMnnou7QKiubFDwdhMLji5d6QGY/+B6OxlGff7Utwo3CeH/mxU+b2cm/Iwj9mep06f4QhPSfxl56TmbFoIEXFF0ifuFHR3Lfz9TDdj0z8ORn0kBRnuuOBLQR6Q08bfVdnjZr2wehW8Ux9uxeHTu7knSnbCQ4Mp6T0ZpVlxdeLHN6f1CYzwOCZDQlv1h6Alx5bjq9ngCJ9iWo+v3h7e1NaWspTTz2ldJRKDp3M48q1S7yZsp5LV89y7cYVhiU8Q2zrBLw9/HkrZSNvJucSGRzH1t2r8GjgTeqoTKVjm80dHnwvbyTn8EZyDqF3xdDlnoEAzH58lYVWhXC82BYQ18r8OpMXW56cotfBiG7QwEGHAOb2QYPehVljPqdnu4cq1q9umaP7k9pmBggNalfRn5RfQUWJvkQ1xc5Z5R3aSMdI04VNO0Tcj15/62R2g96A/o/xkNKyUoIDIxTJWB1zucsV37jC+csnFb/ygRCWPNIVIoMsr1cTHTCym2kY3FHM7YM6na7KCePVLXO02mYGOHpqD1Pe7sn7q5+36hxoe1HNMKazunz1HKu2vMOnG96iqPgCCbHD8fO+tcfsPfoD//tZCm4u7gzrNU3BpJVZyg3w4941dIzqr1BCIaznaoAnEuGzbbDZwu1+7uTjDiO6mr7/cyRr9kFnU5fMGc/tx8fDn79/+iRbdn9BfPQgB6WtTIpdPfl4BvBYv1eIjx7E1t2rOH2x8skpd7fozNuTtrEs9w2++nERQ3tVvQSaEizlBtiU9xnDE59VIJ0QtedqgOFdTMOaq3dYvqScqwE6hcJ/3QteCtxAwJp90NnUJXP50GV82yQOHP9FsWInw5j11Da0BzsL1gOw42AORuOty/CX32gRTNeLc3P1cHi+mpjLDVBSepOjp/bQulmsEvGEqLOou2BKf5jWH/q1gzbNoIkvNPKGEH/oHAbDO8OsIabiqEShA8v7oDOqbebiG1co/WOdXYc30axRa7tnrIkc2dVTaFBbXAyuTFuYSJuW3XB386LUaLp76MET23n3y2fQ6/T4eATw3IiPFE57i7ncAL8cWMe9re9TMKEQ9dO8kenhrCztg2kfDSfv8EaOn9nPw4nPEt92cLXLnDlzE/+WvLHscTzcvAkKCOWvD8xyaN7bSbGzgXED5lb8e/2vWWRmv0ZwYASxrRN4Mzm30rrF14t4bemjRDXv5OiYVZjL3SmqH52i+lVaf8aigQT43uXomEJolrl9cOboT6qsf+cyJfqT2mZeOPnnKsuU6Euc8tqYzkSt14jT8rUxhVATNfYhasxsiRzZWeCj4OSo+ry3UrmV/HsJ4YzU2IeoMbMlcmQnhBBC82Q2phBCCM2TYieEEELzpNgJIYTQPCl2QgghNE+KnRBCCM2TYieEEELzpNgJIYTQPCl2QgghNE+KnRBCCM2TYieEEELzpNgJIYTQPCl2QgghNE+KnRBCCM2TYieEEELzpNgJIYTQPCl2QgghNE+KnRBCCM2TYieEEELz/j+SiJIgyOB6dAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 568.197x264.88 with 1 Axes>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "'''\n",
    "VARIATIONAL ANSATZ\n",
    "Note that for molecules other than H_2 you may need to specify a different number of reps.\n",
    "'''\n",
    "\n",
    "reps = 1   \n",
    "ansatz = EfficientSU2(system_size, su2_gates=['rx', 'ry'], entanglement='circular', reps=reps, skip_final_rotation_layer=False)\n",
    "\n",
    "    \n",
    "ansatz.decompose().draw('mpl')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "c69ce62c",
   "metadata": {},
   "outputs": [],
   "source": [
    "'''\n",
    "COST FUNCTION\n",
    "'''\n",
    "\n",
    "\n",
    "def basis_change_circuit(pauli_op):\n",
    "\n",
    "    #   Generating circuit with just the basis change operators\n",
    "    #\n",
    "    #       pauli_op: n-qubit Pauli operator\n",
    "\n",
    "    basis_change = QuantumCircuit(ansatz.num_qubits, ansatz.num_qubits)\n",
    "    \n",
    "    for idx, op in enumerate(pauli_op):\n",
    "        if op == 'X':\n",
    "            basis_change.h(idx)\n",
    "        elif op == 'Y':\n",
    "            basis_change.h(idx)\n",
    "            basis_change.p(-np.pi/2, idx)\n",
    "\n",
    "    return basis_change\n",
    "\n",
    "\n",
    "\n",
    "def min_cost():\n",
    "    \n",
    "    return min(cost_history)\n",
    "\n",
    "\n",
    "def log_cost(cost):\n",
    "\n",
    "    global cost_history\n",
    "    cost_history.append(cost)\n",
    "    \n",
    "    \n",
    "\n",
    "\n",
    "def objective_function(operators, params):\n",
    "    \n",
    "    backend = QasmSimulator(method='statevector', shots=1)\n",
    "    \n",
    "    pauli_op_dict = Counter(tuple(x) for x in operators)    \n",
    "    \n",
    "    shadow = []\n",
    "    for pauli_op in pauli_op_dict:\n",
    "        \n",
    "        qc = ansatz.bind_parameters(params)\n",
    "        qc = qc.compose(basis_change_circuit(pauli_op))                                   \n",
    "        qc.measure(reversed(range(system_size)), range(system_size))\n",
    "        result = execute(qc, backend, shots=pauli_op_dict[pauli_op]).result()\n",
    "        counts = result.get_counts()    \n",
    "        \n",
    "        # generate a separate circuit for each basis change Pauli operator\n",
    "        output_str = list(list(counts.keys())[list(counts.values()).index(1)])\n",
    "        output = [int(i) for i in output_str]\n",
    "        eigenvals = [x+1 if x == 0 else x-2 for x in output]\n",
    "        snapshot = [(op, eigenval) for op, eigenval in zip(pauli_op, eigenvals)] \n",
    "                \n",
    "        shadow.append(snapshot)\n",
    "    \n",
    "    \n",
    "    cost = 0.0\n",
    "    for term, weight in zip(observables_xyze, hamiltonian_coefficients):\n",
    "        sum_product, match_count = estimate_exp(shadow, term)\n",
    "        if match_count != 0:\n",
    "            exp_val = sum_product / match_count\n",
    "            cost += weight * exp_val\n",
    "            \n",
    "    log_cost(cost)\n",
    "    \n",
    "    return cost\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ea2eb6b9",
   "metadata": {},
   "outputs": [],
   "source": [
    "'''\n",
    "RUNNING EXPERIMENTS\n",
    "''' \n",
    "\n",
    "start_time = time.time()\n",
    "rmse_errors = []\n",
    "\n",
    "print('NUMBER OF OPERATORS | RANDOMISED OPERATORS | AVERAGE RMSE ERROR\\n')\n",
    "measurement_range = [50, 100, 200, 600, 1000, 1600]\n",
    "\n",
    "for num_operators in measurement_range:   \n",
    "    \n",
    "    basis_change = randomized_classical_shadow(num_operators, system_size)\n",
    "\n",
    "    tuples = (tuple(pauli) for pauli in basis_change)\n",
    "    counts = Counter(tuples)\n",
    "\n",
    "    optimizer = SPSA(maxiter=2000) \n",
    "    cost_function = partial(objective_function, basis_change)\n",
    "    \n",
    "    expectation_values = []\n",
    "    num_experiments = 5\n",
    "\n",
    "\n",
    "    for iteration in range(num_experiments):\n",
    "        \n",
    "        cost_history = []\n",
    "        \n",
    "        params = np.random.rand(ansatz.num_parameters)\n",
    "        result = optimizer.minimize(fun=cost_function, x0=params)\n",
    "        minimal_cost = min_cost()\n",
    "        expectation_values.append(minimal_cost)\n",
    "                \n",
    "        print(\"EXPERIMENT {}: GROUND STATE ENERGY FOUND = {}\".format(iteration, minimal_cost))\n",
    "\n",
    "\n",
    "    rmse_randomised_cs = np.sqrt(np.sum([(expectation_values[i] - EXPECTED_EIGENVALUE)**2 \n",
    "                                           for i in range(num_experiments)])/num_experiments)\n",
    "\n",
    "    rmse_errors.append(rmse_randomised_cs)\n",
    "    \n",
    "    print('{} | {} | {}'.format(num_operators, counts, rmse_randomised_cs))\n",
    "\n",
    "\n",
    "elapsed_time = time.time() - start_time\n",
    "print(\"Execution time = \", time.strftime(\"%H:%M:%S\", time.gmtime(elapsed_time)))\n",
    "\n",
    "          \n",
    "points = measurement_range\n",
    "num_points = len(measurement_range)\n",
    "\n",
    "plt.plot([i for i in points], [rmse_errors[i] for i in range(num_points)], 'r', label='SPSA(maxiter=2000)')\n",
    "plt.xlabel('Number of measurements')\n",
    "plt.ylabel('Average RMSE error')\n",
    "plt.legend(loc=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "685faa90",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.8.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}