demos/qiskit_alt_demo.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "3e23cbd3",
   "metadata": {},
   "source": [
    "# qiskit_alt\n",
    "\n",
    "This python package provides a thin wrapper around some features of Qiskit that have been (re-)implemented in Julia and provides a Python interface. The input and output of the Python interface are the same as the input and output to Python qiskit. At present, we have prepared two high level demonstrations\n",
    "\n",
    "* Performing the Jordan-Wigner transform from a Fermionic operator to a Pauli operator.\n",
    "\n",
    "* Computing the Fermionic operator from integrals computed by `pyscf`.\n",
    "\n",
    "In both cases, we will see that the `qiskit_alt` implementation is much more performant.\n",
    "\n",
    "We have also prepared some lower-level demonstrations of performance gains\n",
    "\n",
    "* Converting an operator from the computational basis to the Pauli basis.\n",
    "\n",
    "* Creating a `SparsePauliOp` from a list of strings"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "86f7fd8f",
   "metadata": {},
   "source": [
    "The Python package has been installed in a virtual environment created with `python -m venv ./env`. The Julia packages have been installed in a local environment in the standard way, via a spec in `Project.toml` file.\n",
    "\n",
    "When we import the package `qiskit_alt`, the Julia environment is also activated.\n",
    "There are two options for communcating with Julia: `PyCall.jl/pyjulia` and `PythonCall.jl/juliacall`.\n",
    "Here we use the first by passing `calljulia=\"pyjulia\"` when initializing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "29b7099b",
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import qiskit_alt\n",
    "qiskit_alt.project.ensure_init(calljulia=\"pyjulia\", compile=False, depot=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "cb29172e",
   "metadata": {},
   "source": [
    "We assume that no one is familiar with Julia, much less with `pyjulia`, the package we use to call Julia from Python. So, we inject a bit of tutorial.\n",
    "\n",
    "The default `Module` in Julia `Main` is available. You can think of it as a namespace. And, as always, objects from the `Module` `Base` have been imported into `Main`.\n",
    "\n",
    "As an example of how `pyjulia` works, we create an `Array` of `Float64` zeros on the Julia side, and they are automatically copied to a numpy array when returned to Python."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "0f936bc8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0.])"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "julia = qiskit_alt.project.julia\n",
    "julia.Main.zeros(10)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b733f29d",
   "metadata": {},
   "source": [
    "There are several ways to call Julia from Python and vice versa, and to specifiy features such as the copying vs. sharing semantics. We won't go into much of this in this demo."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dea87cf9",
   "metadata": {},
   "source": [
    "## Electronic structure\n",
    "\n",
    "Part of a workflow for, say, VQE involves using qiskit-nature to do the following:\n",
    "* Get a description of a model Hamiltonian from the package `pyscf` by passing it a description of the geometry of a molecule.\n",
    "* Convert that description of a Hamiltonian to a qiskit-native Fermionic operator.\n",
    "* Convert the Fermionic operator to a qubit operator expressed in the Pauli basis.\n",
    "\n",
    "The last step above may be done in several ways, one of which is known as the Jordan-Wigner transform. It is this step that we will benchmark here."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c4ca7000",
   "metadata": {},
   "source": [
    "### qiskit-nature\n",
    "\n",
    "First, we see how this is done in qiskit-nature. We need to specify the geometry of the molecule and the\n",
    "[basis set](https://en.wikipedia.org/wiki/Basis_set_(chemistry)). We choose `sto3g`, one of the smallest, simplest, basis sets."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "ca1047f9",
   "metadata": {},
   "outputs": [],
   "source": [
    "from qiskit_nature.drivers import UnitsType, Molecule\n",
    "from qiskit_nature.drivers.second_quantization import ElectronicStructureDriverType, ElectronicStructureMoleculeDriver\n",
    "\n",
    "# Specify the geometry of the H_2 molecule\n",
    "geometry = [['H', [0., 0., 0.]],\n",
    "            ['H', [0., 0., 0.735]]]\n",
    "\n",
    "basis = 'sto3g'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c6726e7d",
   "metadata": {},
   "source": [
    "Then, we compute the fermionic Hamiltonian like this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "fc1697b6",
   "metadata": {},
   "outputs": [],
   "source": [
    "molecule = Molecule(geometry=geometry,\n",
    "                     charge=0, multiplicity=1)\n",
    "driver = ElectronicStructureMoleculeDriver(molecule, basis=basis, driver_type=ElectronicStructureDriverType.PYSCF)\n",
    "\n",
    "from qiskit_nature.problems.second_quantization import ElectronicStructureProblem\n",
    "from qiskit_nature.converters.second_quantization import QubitConverter\n",
    "from qiskit_nature.mappers.second_quantization import JordanWignerMapper\n",
    "\n",
    "es_problem = ElectronicStructureProblem(driver)\n",
    "second_q_op = es_problem.second_q_ops()\n",
    "\n",
    "fermionic_hamiltonian = second_q_op[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f6000f0a",
   "metadata": {},
   "source": [
    "The Jordan-Wigner transform is performed like this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "1d41dd67",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "SparsePauliOp(['YYYY', 'XXYY', 'YYXX', 'XXXX', 'IIII', 'ZIII', 'IZII', 'ZZII', 'IIZI', 'ZIZI', 'IZZI', 'IIIZ', 'ZIIZ', 'IZIZ', 'IIZZ'],\n",
       "              coeffs=[ 0.0452328 +0.j,  0.0452328 +0.j,  0.0452328 +0.j,  0.0452328 +0.j,\n",
       " -0.81054798+0.j, -0.22575349+0.j,  0.17218393+0.j,  0.12091263+0.j,\n",
       " -0.22575349+0.j,  0.17464343+0.j,  0.16614543+0.j,  0.17218393+0.j,\n",
       "  0.16614543+0.j,  0.16892754+0.j,  0.12091263+0.j])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "qubit_converter = QubitConverter(mapper=JordanWignerMapper())\n",
    "nature_qubit_op = qubit_converter.convert(fermionic_hamiltonian)\n",
    "nature_qubit_op.primitive"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f00c5538",
   "metadata": {},
   "source": [
    "### qiskit_alt\n",
    "\n",
    "The only high-level code in `qiskit_alt` was written to support this demo. So doing the JW-transform is less verbose."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "7ecebbe7",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "SparsePauliOp(['IIII', 'ZIII', 'IZII', 'ZZII', 'IIZI', 'ZIZI', 'IZZI', 'XXXX', 'YYXX', 'XXYY', 'YYYY', 'IIIZ', 'ZIIZ', 'IZIZ', 'IIZZ'],\n",
       "              coeffs=[-0.09057899+0.j, -0.22575349+0.j,  0.17218393+0.j,  0.12091263+0.j,\n",
       " -0.22575349+0.j,  0.17464343+0.j,  0.16614543+0.j,  0.0452328 +0.j,\n",
       "  0.0452328 +0.j,  0.0452328 +0.j,  0.0452328 +0.j,  0.17218393+0.j,\n",
       "  0.16614543+0.j,  0.16892754+0.j,  0.12091263+0.j])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import qiskit_alt.electronic_structure\n",
    "fermi_op = qiskit_alt.electronic_structure.fermionic_hamiltonian(geometry, basis)\n",
    "pauli_op = qiskit_alt.electronic_structure.jordan_wigner(fermi_op)\n",
    "pauli_op.simplify() # The Julia Pauli operators use a different sorting convention; we sort again for comparison. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3d6b0ffa",
   "metadata": {},
   "source": [
    "Note that the constant term differs. The qiskit-nature version ignores the nuclear-repulsion term. I need to open an issue about whether and how to handle it."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f32494c5",
   "metadata": {},
   "source": [
    "### Benchmarking\n",
    "\n",
    "Computing the Hamiltonian for a larger molecule or a larger basis set takes more time and produces a Hamiltonian with more factors and terms. Here we compare the performance of `qiskit_alt` and `qiskit-nature` on combinations of $\\text{H}_2$ and $\\text{H}_2\\text{O}$ molecules for several basis sets."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9c118b30",
   "metadata": {},
   "source": [
    "First we benchmark qiskit-nature, and record the times in `nature_times`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "83562a9c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "geometry=h2_geometry, basis='sto3g' 3.37 ms\n",
      "geometry=h2_geometry, basis='631g' 39.25 ms\n",
      "geometry=h2_geometry, basis='631++g' 209.80 ms\n",
      "geometry=h2o_geometry, basis='sto3g' 246.12 ms\n",
      "geometry=h2o_geometry, basis='631g' 3004.95 ms\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[3.371502598747611,\n",
       " 39.25394311081618,\n",
       " 209.79926399886608,\n",
       " 246.11722519621253,\n",
       " 3004.945662105456]"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%run ../bench/jordan_wigner_nature_time.py\n",
    "nature_times"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c6f6e210",
   "metadata": {},
   "source": [
    "Next we benchmark qiskit_alt, and record the times in `alt_times`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "56f3abc4",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "geometry=h2_geometry, basis='sto3g' 1.30 ms\n",
      "geometry=h2_geometry, basis='631g' 3.44 ms\n",
      "geometry=h2_geometry, basis='631++g' 24.52 ms\n",
      "geometry=h2o_geometry, basis='sto3g' 18.70 ms\n",
      "geometry=h2o_geometry, basis='631g' 469.48 ms\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[1.3044350082054734,\n",
       " 3.4391670022159815,\n",
       " 24.519264418631792,\n",
       " 18.69659200310707,\n",
       " 469.4802389945835]"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%run ../bench/jordan_wigner_alt_time.py\n",
    "alt_times"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "74a38fe5",
   "metadata": {},
   "source": [
    "We compare the relative performance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "8269b475",
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[2.5846459022790462,\n",
       " 11.413793830169755,\n",
       " 8.556507259632268,\n",
       " 13.163747979060131,\n",
       " 6.400579646420698]"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "[t_nature / t_qk_alt for t_nature, t_qk_alt in zip(nature_times, alt_times)]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e7e0524b",
   "metadata": {},
   "source": [
    "We see that\n",
    "* qiskit_alt is at least ten times faster\n",
    "* The relative performance increases as the problem in some sense gets larger.\n",
    "\n",
    "In fact, another problem, not shown here, finishes in 18s with qiskit_alt and in 5730s in qiskit-nature.\n",
    "In this case, `qiskit_alt` is 320 times faster than `qiskit-nature`. I don't have an idea about the origin of this scaling."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f6ca3e40",
   "metadata": {},
   "source": [
    "### Computing the Fermonic operator\n",
    "\n",
    "Computing the Fermionic operator from the output of `pyscf` is also much more efficient in `qiskit_alt`."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f2dd8cfa",
   "metadata": {},
   "source": [
    "We benchmark qiskit-nature computing the fermionic Hamiltonian"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "1eb2798b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "geometry=h2_geometry, basis='sto3g' 65.00 ms\n",
      "geometry=h2_geometry, basis='631g' 105.01 ms\n",
      "geometry=h2_geometry, basis='631++g' 298.17 ms\n",
      "geometry=h2o_geometry, basis='sto3g' 438.85 ms\n",
      "geometry=h2o_geometry, basis='631g' 24030.27 ms\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[65.00064351130277,\n",
       " 105.00940189231187,\n",
       " 298.17473497241735,\n",
       " 438.84517257101834,\n",
       " 24030.2699431777]"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%run ../bench/fermionic_nature_time.py\n",
    "nature_times"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a56c2fd1",
   "metadata": {},
   "source": [
    "We benchmark qiskit_alt computing the fermionic Hamiltonian."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "b6650d69",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "geometry=h2_geometry, basis='sto3g' 59.39 ms\n",
      "geometry=h2_geometry, basis='631g' 62.03 ms\n",
      "geometry=h2_geometry, basis='631++g' 67.03 ms\n",
      "geometry=h2o_geometry, basis='sto3g' 77.55 ms\n",
      "geometry=h2o_geometry, basis='631g' 193.82 ms\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[59.38674798235297,\n",
       " 62.029914115555584,\n",
       " 67.0257663121447,\n",
       " 77.55265950690955,\n",
       " 193.81799241527915]"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%run ../bench/fermionic_alt_time.py\n",
    "alt_times"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0c8e4ac9",
   "metadata": {},
   "source": [
    "We compare the relative performance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "e197f5fb",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1.094531115437033,\n",
       " 1.6928832385079522,\n",
       " 4.448658350040971,\n",
       " 5.658673414441957,\n",
       " 123.9836902844905]"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "[t_nature / t_qk_alt for t_nature, t_qk_alt in zip(nature_times, alt_times)]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ef5d18bf",
   "metadata": {},
   "source": [
    "We see again that, as the problem size increases, `qiskit_alt` is increasingly more performant."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bdfd094a",
   "metadata": {},
   "source": [
    "## Discussion\n",
    "\n",
    "The Julia implemenation consists of these packages\n",
    "\n",
    "* [`QuantumOps.jl`](https://github.com/Qiskit-Extensions/QuantumOps.jl) implementing Fermionic and Pauli operators and calculations using them.\n",
    "\n",
    "* [`ElectronicStructure.jl`](https://github.com/Qiskit-Extensions/ElectronicStructure.jl) provides an interface to electronic structure packages.\n",
    "\n",
    "* [`ElectronicStructurePySCF.jl`](https://github.com/Qiskit-Extensions/ElectronicStructurePySCF.jl) provides an interface to `pyscf`\n",
    "\n",
    "* [`QiskitQuantumInfo.jl`](https://github.com/Qiskit-Extensions/QiskitQuantumInfo.jl) provides data structures that mirror Python Qiskit data structures. These are used as intermedidate structures for converting from `QuantumOps` and `ElectronicStructure` to Python Qiskit. In the future these might be used directly for calculations.\n",
    "\n",
    "\n",
    "The Python interface is a Python package `qiskit_alt`. This could contain a mixture of Julia and Python code. Or all the Julia code might be moved to the Julia packages.\n",
    "\n",
    "### Implementation\n",
    "\n",
    "In the examples above, the following happens.\n",
    "\n",
    "* Julia code calls `pyscf` and stores the results in Julia data structures.\n",
    "\n",
    "* These data are used to construct a Fermionic operator as a data structure defined in `QuantumOps`.\n",
    "\n",
    "* The Jordan-Wigner transform, implemented in `QuantumOps` is used to compute a Pauli operator.\n",
    "\n",
    "* The Pauli operator (as a structure in `QuantumOps`) is converted to a Qiskit-like operator defined in `QiskitQuantumInfo.jl`.\n",
    "\n",
    "* The operator defined in `QiskitQuantumInfo.jl` is sent to Python and converted to numpy arrays, which are then used to construct native Qiskit types. The conversion to numpy arrays is provided by `pyjulia`.\n",
    "\n",
    "### Complexity, dynamism\n",
    "\n",
    "* It is worth noting that operators in `QuantumOps` are *not* highly optimized implementations. In fact, much of the code for the two types of operators is shared, they inherit from a parent class. There are other implementations of Pauli operators in Julia that are much more efficient for instance in [`QuantumClifford.jl`](https://github.com/Krastanov/QuantumClifford.jl).\n",
    "\n",
    "* [Issue](https://github.com/Qiskit-Extensions/QuantumOps.jl/issues/17) for improving performance of Jordan-Wigner in `QuantumOps`.\n",
    "    * Precompute one and two-body terms\n",
    "    * Use @threads\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "af37da81",
   "metadata": {},
   "source": [
    "# More demos\n",
    "\n",
    "Here are some smaller scale demonstrations.\n",
    "\n",
    "## Converting a matrix to the Pauli basis\n",
    "\n",
    "Here we convert a matrix representing an operator in the computational basis to the Pauli basis.\n",
    "In this case, `qiskit_alt` is much more performant than `qiskit.quantum_info`.\n",
    "This is how it is done in `QuantumOps`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "9bda9915",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "SparsePauliOp(['III', 'XII', 'YII', 'ZII', 'IXI', 'XXI', 'YXI', 'ZXI', 'IYI', 'XYI', 'YYI', 'ZYI', 'IZI', 'XZI', 'YZI', 'ZZI', 'IIX', 'XIX', 'YIX', 'ZIX', 'IXX', 'XXX', 'YXX', 'ZXX', 'IYX', 'XYX', 'YYX', 'ZYX', 'IZX', 'XZX', 'YZX', 'ZZX', 'IIY', 'XIY', 'YIY', 'ZIY', 'IXY', 'XXY', 'YXY', 'ZXY', 'IYY', 'XYY', 'YYY', 'ZYY', 'IZY', 'XZY', 'YZY', 'ZZY', 'IIZ', 'XIZ', 'YIZ', 'ZIZ', 'IXZ', 'XXZ', 'YXZ', 'ZXZ', 'IYZ', 'XYZ', 'YYZ', 'ZYZ', 'IZZ', 'XZZ', 'YZZ', 'ZZZ'],\n",
       "              coeffs=[ 0.47494296+0.00000000e+00j,  0.40823308+0.00000000e+00j,\n",
       "  0.        +8.72050729e-03j,  0.08853865+0.00000000e+00j,\n",
       "  0.48022917+0.00000000e+00j,  0.50427456+0.00000000e+00j,\n",
       "  0.        +1.86142504e-01j,  0.06825203+0.00000000e+00j,\n",
       "  0.        -1.20594387e-01j,  0.        +1.63950963e-01j,\n",
       "  0.0933073 +0.00000000e+00j,  0.        +3.95169450e-02j,\n",
       " -0.08316598+0.00000000e+00j, -0.00272138+0.00000000e+00j,\n",
       "  0.        -7.34809038e-02j,  0.06674014+0.00000000e+00j,\n",
       "  0.58852484+0.00000000e+00j,  0.51606481+0.00000000e+00j,\n",
       "  0.        +8.79202746e-02j,  0.03380507+0.00000000e+00j,\n",
       "  0.52470279+0.00000000e+00j,  0.57627299+0.00000000e+00j,\n",
       "  0.        -3.41978024e-02j,  0.20037087+0.00000000e+00j,\n",
       "  0.        +4.90473087e-02j,  0.        -1.21577947e-01j,\n",
       "  0.04581339+0.00000000e+00j,  0.        -1.73195154e-01j,\n",
       "  0.0348815 +0.00000000e+00j,  0.16682054+0.00000000e+00j,\n",
       "  0.        +3.50230791e-02j,  0.00527009+0.00000000e+00j,\n",
       "  0.        +1.79330324e-01j,  0.        -8.29205184e-02j,\n",
       "  0.13369856+0.00000000e+00j,  0.        +9.74266701e-02j,\n",
       "  0.        +5.62465253e-03j,  0.        -1.35618682e-01j,\n",
       "  0.14589336+0.00000000e+00j,  0.        -1.29918791e-01j,\n",
       " -0.06372318-0.00000000e+00j, -0.1256793 -0.00000000e+00j,\n",
       "  0.        +1.78768198e-01j,  0.0769969 +0.00000000e+00j,\n",
       "  0.        -1.21335084e-01j,  0.        +1.96184160e-01j,\n",
       "  0.05176287+0.00000000e+00j,  0.        +8.92921828e-02j,\n",
       " -0.0737263 +0.00000000e+00j,  0.06150848+0.00000000e+00j,\n",
       "  0.        +1.37063194e-01j,  0.04306766+0.00000000e+00j,\n",
       " -0.26085734+0.00000000e+00j,  0.14409961+0.00000000e+00j,\n",
       "  0.        +7.95247641e-02j, -0.0744302 +0.00000000e+00j,\n",
       "  0.        -3.14516650e-02j,  0.        -1.15053673e-02j,\n",
       "  0.02982605+0.00000000e+00j,  0.        -1.85210774e-04j,\n",
       " -0.09773203+0.00000000e+00j,  0.03884013+0.00000000e+00j,\n",
       "  0.        +5.20931164e-03j, -0.31915262+0.00000000e+00j])"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from qiskit_alt.pauli_operators import QuantumOps, PauliSum_to_SparsePauliOp\n",
    "import numpy as np\n",
    "m = np.random.rand(2**3, 2**3) # 3-qubit operator\n",
    "pauli_sum = QuantumOps.PauliSum(m) # This is a wrapped Julia object\n",
    "PauliSum_to_SparsePauliOp(pauli_sum) # Convert to qiskit.quantum_info.SparsePauliOp"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0bfdfacb",
   "metadata": {},
   "source": [
    "Note that the `numpy` matrix was automatically converted to a Julia `Matrix`."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9268c181",
   "metadata": {},
   "source": [
    "We run benchmarks of conversion of matrices to the Pauli basis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "a6945b8f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "nqubits=2, 1.23 ms\n",
      "nqubits=3, 5.00 ms\n",
      "nqubits=4, 24.66 ms\n",
      "nqubits=5, 142.61 ms\n",
      "nqubits=6, 1053.21 ms\n",
      "nqubits=7, 8999.82 ms\n",
      "nqubits=8, 107648.63 ms\n"
     ]
    }
   ],
   "source": [
    "%run ../bench/from_matrix_quantum_info.py"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "465a7bf5",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "nqubits=2, 1.29 ms\n",
      "nqubits=3, 1.36 ms\n",
      "nqubits=4, 11.33 ms\n",
      "nqubits=5, 27.75 ms\n",
      "nqubits=6, 27.85 ms\n",
      "nqubits=7, 107.10 ms\n",
      "nqubits=8, 567.22 ms\n"
     ]
    }
   ],
   "source": [
    "%run ../bench/from_matrix_alt.py"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "94fd017c",
   "metadata": {},
   "source": [
    "Here are the ratios of the times for `qiskit.quantum_info` to those for `qiskit_alt`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "217f2f7a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0.9556047308329915,\n",
       " 3.67816523510366,\n",
       " 2.1765062775946946,\n",
       " 5.138368881804313,\n",
       " 37.8119442420315,\n",
       " 84.0323151825261,\n",
       " 189.7838578070803]"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "[t_pyqk / t_qk_alt for t_pyqk, t_qk_alt in zip(pyqk_times, qk_alt_times)]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "40fb432c",
   "metadata": {},
   "source": [
    "Again, the performance gain increases with the problem size."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0cd10e98",
   "metadata": {},
   "source": [
    "## Creating a `SparsePauliOp` from a list of strings\n",
    "\n",
    "\n",
    "Here, we create a `SparsePauliOp` from a list of `n` strings, each with `k` single-Pauli factors, and simplify the result."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fbfeab47",
   "metadata": {},
   "source": [
    "First, using `qiskit.quantum_info`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "39ad43f3",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "k=10, n=10, 0.5427255993708968 ms\n",
      "k=10, n=100, 1.8118745530955493 ms\n",
      "k=10, n=1000, 16.81242184713483 ms\n",
      "k=10, n=5000, 83.37073109578341 ms\n",
      "k=10, n=10000, 167.40924795158207 ms\n",
      "k=10, n=100000, 1695.4269821289927 ms\n",
      "k=100, n=10, 0.336431001778692 ms\n",
      "k=100, n=100, 1.897125702816993 ms\n",
      "k=100, n=1000, 17.57452649762854 ms\n",
      "k=100, n=5000, 87.75187629507855 ms\n",
      "k=100, n=10000, 176.62636120803654 ms\n",
      "k=100, n=100000, 1799.1630749311298 ms\n"
     ]
    }
   ],
   "source": [
    "%run ../bench/pauli_from_list_qinfo.py"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bc7f9d07",
   "metadata": {},
   "source": [
    "Now, using `qiskit_alt`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "a4b9eeda",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "k=10, n=10, 1.2364529073238373 ms\n",
      "k=10, n=100, 1.3390239444561303 ms\n",
      "k=10, n=1000, 2.580544597003609 ms\n",
      "k=10, n=5000, 8.644396206364036 ms\n",
      "k=10, n=10000, 25.033850746694952 ms\n",
      "k=10, n=100000, 228.47767625935376 ms\n",
      "k=100, n=10, 5.032895249314606 ms\n",
      "k=100, n=100, 1.4746043947525322 ms\n",
      "k=100, n=1000, 3.6469499580562115 ms\n",
      "k=100, n=5000, 17.004536895547062 ms\n",
      "k=100, n=10000, 33.38807285763323 ms\n",
      "k=100, n=100000, 370.6562086008489 ms\n"
     ]
    }
   ],
   "source": [
    "%run ../bench/pauli_from_list_alt.py"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a1372173",
   "metadata": {},
   "source": [
    "The results were written to lists `quantum_info_times` and `qkalt_times`. We compare the performance:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "e7bfa245",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0.43893754154015063,\n",
       " 1.3531308088978768,\n",
       " 6.515067349216332,\n",
       " 9.644482865605532,\n",
       " 6.687315093691049,\n",
       " 7.420536701381927,\n",
       " 0.06684641446183648,\n",
       " 1.286531973977582,\n",
       " 4.818965628745723,\n",
       " 5.16049786207691,\n",
       " 5.2901035037623005,\n",
       " 4.853994168133864]"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "[x / y for x,y in zip(quantum_info_times, qkalt_times)]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7c28aadc",
   "metadata": {},
   "source": [
    "We see that the performance improvement in `qiskit_alt` is significant, but does not increase with the number of terms `n`. Further benchmarks show that the time required to convert the strings from Python to Julia takes all the time.\n",
    "\n",
    "We see this in the following."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "43d88903",
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Initializing Julia interpreter. This may take some time...\n"
     ]
    }
   ],
   "source": [
    "%load_ext julia.magic"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0b82727a",
   "metadata": {},
   "source": [
    "Generate `1000` ten-qubit Pauli strings."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "22ac98aa",
   "metadata": {},
   "outputs": [],
   "source": [
    "%julia using Random: randstring\n",
    "%julia pauli_strings = [randstring(\"IXYZ\", 10) for _ in 1:1000]\n",
    "None;"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "77f21d2d",
   "metadata": {},
   "source": [
    "Benchmark converting these to a `QuantumOps.PauliSum`. Note that as the sums are always sorted."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "4c1dade6",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "    Updating registry at `~/julia_qiskit_repos/qiskit_alt/venv/julia_project/qiskit_alt-1.7.2/depot/registries/QuantumRegistry`\n",
      "    Updating git-repo `https://github.com/Qiskit-Extensions/QuantumRegistry`\n",
      "    Updating registry at `~/julia_qiskit_repos/qiskit_alt/venv/julia_project/qiskit_alt-1.7.2/depot/registries/General.toml`\n",
      "   Resolving package versions...\n",
      "  No Changes to `~/julia_qiskit_repos/qiskit_alt/venv/julia_project/qiskit_alt-1.7.2/Project.toml`\n",
      "  No Changes to `~/julia_qiskit_repos/qiskit_alt/venv/julia_project/qiskit_alt-1.7.2/Manifest.toml`\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  341.921 μs (1012 allocations: 134.17 KiB)\n"
     ]
    }
   ],
   "source": [
    "%julia import Pkg; Pkg.add(\"BenchmarkTools\")\n",
    "%julia using BenchmarkTools: @btime\n",
    "%julia @btime QuantumOps.PauliSum($pauli_strings)\n",
    "None;"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5b1facf8",
   "metadata": {},
   "source": [
    "Check that we are actually getting the right result."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "ff503a66",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "999\n",
      "10-factor QuantumOps.PauliTerm{Vector{QuantumOps.Paulis.Pauli}, Complex{Int64}}:\n",
      "IIIIXZZIZZ * (1 + 0im)\n"
     ]
    }
   ],
   "source": [
    "%julia pauli_sum = QuantumOps.PauliSum(pauli_strings);\n",
    "%julia println(length(pauli_sum))\n",
    "%julia println(pauli_sum[1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "6255a2ef",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "46.06705539358601"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "6.9 * 2.29 / .343  # Ratio of time to construct PauliSum via qiskit_alt to time in pure Julia"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d2c36826",
   "metadata": {},
   "source": [
    "So the pure Julia code is `46` times faster than the qiskit.quantum_info.\n",
    "**But, the `qiskit.quantum_info` is also extracting a possible phase !**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "f3d77af7",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/lapeyre/julia_qiskit_repos/qiskit_alt/venv/lib/python3.9/site-packages/qiskit/tools/jupyter/__init__.py:133: RuntimeWarning: matplotlib can't be found, ensure you have matplotlib and other visualization dependencies installed. You can run '!pip install qiskit-terra[visualization]' to install it from jupyter\n",
      "  warnings.warn(\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<h3>Version Information</h3><table><tr><th>Qiskit Software</th><th>Version</th></tr><tr><td><code>qiskit-terra</code></td><td>0.20.0</td></tr><tr><td><code>qiskit_alt</code></td><td>0.1.13</td></tr><tr><td><code>qiskit-nature</code></td><td>0.3.1</td></tr><tr><th>System information</th></tr><tr><td>Python version</td><td>3.9.7</td></tr><tr><td>Python compiler</td><td>GCC 11.1.0</td></tr><tr><td>Python build</td><td>default, Oct 10 2021 15:13:22</td></tr><tr><td>OS</td><td>Linux</td></tr><tr><td>CPUs</td><td>12</td></tr><tr><td>Memory (Gb)</td><td>62.76226043701172</td></tr><tr><td colspan='2'>Tue Apr 05 11:13:16 2022 EDT</td></tr></table>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import qiskit.tools.jupyter\n",
    "d = qiskit.__qiskit_version__._version_dict\n",
    "d['qiskit_alt'] = qiskit_alt.__version__\n",
    "%qiskit_version_table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fc683f8e",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "qiskit_env_march2022",
   "language": "python",
   "name": "qiskit_env_march2022"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}