Qiskit-Pauli_decomposing.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 67,
   "outputs": [
    {
     "data": {
      "text/plain": "<IPython.core.display.HTML object>",
      "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.19.2</td></tr><tr><td><code>qiskit-aer</code></td><td>0.10.3</td></tr><tr><td><code>qiskit-ignis</code></td><td>0.7.0</td></tr><tr><td><code>qiskit-ibmq-provider</code></td><td>0.18.3</td></tr><tr><td><code>qiskit</code></td><td>0.34.2</td></tr><tr><td><code>qiskit-nature</code></td><td>0.3.1</td></tr><tr><td><code>qiskit-optimization</code></td><td>0.3.1</td></tr><tr><td><code>qiskit-machine-learning</code></td><td>0.3.1</td></tr><tr><th>System information</th></tr><tr><td>Python version</td><td>3.8.12</td></tr><tr><td>Python compiler</td><td>GCC 7.5.0</td></tr><tr><td>Python build</td><td>default, Oct 12 2021 13:49:34</td></tr><tr><td>OS</td><td>Linux</td></tr><tr><td>CPUs</td><td>24</td></tr><tr><td>Memory (Gb)</td><td>62.71316146850586</td></tr><tr><td colspan='2'>Fri Jun 10 14:12:10 2022 +03</td></tr></table>"
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import qiskit.tools.jupyter\n",
    "%qiskit_version_table"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "# Decomposing Hamiltonian to Pauli observables\n",
    "\n",
    "The $n$-fold Pauli operator set is defined as $G_n=\\{I,X,Y,Z \\}^{\\otimes n}$, that is as the set containing all the possible tensor products between $n$ Pauli matrices. This set of tensor products form an orthogonal basis for the vector space of $2^n \\times 2^n$ complex matrices. To see this first we notice that the space has a dimension of $4^n$ and we also have $4^n$ vectors ( the vectors are operators in this case). So we only need to show that they are linearly independent. We can actually show something stronger. It can be easily seen that the members of the Pauli group are orthogonal under the Hilbert-Schmidt inner product. The H-S inner product of two matrices is defined as $\\operatorname{Tr}(AB^\\dagger)$. We can easily verify from the definition that the Pauli group is a mutually orthogonal set under this inner product. We simply have to use the elementary property $\\operatorname{Tr}(C \\otimes D) = \\operatorname{Tr}(C)\\operatorname{Tr}(D)$, which means $(\\mathbb{C}^{2\\times 2})^{\\otimes n}=\\mathbb{C}^{2^n\\times 2^n}$\n",
    "\n",
    "\n",
    "In summary, products of Pauli matrices form a basis for the set of Hermitian matrices of dimensions that are powers of $2$.\n",
    "\n",
    "More specifically, fix an integer $n$ and let $N\\equiv 2^n$, define $\\boldsymbol \\sigma\\equiv(\\sigma_x,\\sigma_y,\\sigma_y)$, and $\\tilde \\boldsymbol \\sigma\\equiv (I,\\sigma_x,\\sigma_y,\\sigma_y)$. Moreover, let $J\\in\\{0,1,2,3\\}^{n}$ be a tuple of $n$ integers with each $J_i\\in\\{0,1,2,3\\}$. Consider the matrices of the form\n",
    "\n",
    "$$\n",
    "\\tilde \\boldsymbol \\sigma_J\\equiv\\prod_{k=1}^n \\tilde \\boldsymbol \\sigma^{(k)}_{J_k},\n",
    "$$\n",
    "\n",
    "where $\\tilde \\boldsymbol \\sigma_j^{(k)}$ denotes the Pauli matrix $\\tilde \\boldsymbol \\sigma_j$ applied on the $k$-th qubit. You can check that these are all Hermitians. Moreover, $\\tilde \\boldsymbol \\sigma_J$ is traceless for all $J\\neq(0,...,0)$ (for which $\\tilde \\boldsymbol \\sigma_{(0,...,0)}=I$).\n",
    "\n",
    "There are $4^n=N^2$ matrices of this form (one for each possible choice of $J$), and for any pair of tuples $J,K$ we have $\\operatorname{Tr}(\\tilde \\boldsymbol \\sigma_J\\tilde \\boldsymbol \\sigma_K)=N\\delta_{JK}$. Moreover, the space of Hermitian $N\\times N$ matrices also has dimension $N^2$. It follows that $\\{\\tilde \\boldsymbol \\sigma_J\\}_J$ is a basis for this space. Explicitly, you can decompose an arbitrary Hermitian matrix $H$ as\n",
    "\n",
    "$$\n",
    "H=\\frac{1}{N}\\sum_J \\operatorname{Tr}(\\tilde \\boldsymbol \\sigma_J H)\\tilde \\boldsymbol \\sigma_J.\n",
    "$$\n",
    "\n",
    "Note that the coefficients in any such expansion are always real. This is not by chance. Indeed, being more careful, we should state that the set of real linear combinations of products of Pauli matrices give the set of Hermitian matrices (notice that if $A$ is Hermitian then $\\alpha A$ is Hermitian iff $\\alpha\\in\\mathbb R$, so this is not surprising).\n",
    "\n",
    "More general matrices can be generated if we allow for complex coefficients in the expansion. Indeed, note that\n",
    "\n",
    "$$\n",
    "\\frac{I+Z}{2}=\\begin{pmatrix}1&0 \\\\ 0 & 0\\end{pmatrix},\n",
    "\\qquad \\frac{I-Z}{2}=\\begin{pmatrix}0&0 \\\\ 0 & 1\\end{pmatrix}, \\\\\n",
    "\\frac{X+iY}{2}=\\begin{pmatrix}0 & 1 \\\\ 0 & 0\\end{pmatrix}\n",
    "\\qquad\n",
    "\\frac{X-iY}{2}=\\begin{pmatrix}0 & 0 \\\\ 1 & 0\\end{pmatrix}.\n",
    "$$\n",
    "\n",
    "Therefore the complex span of Pauli matrices can be used to generate arbitrary $2\\times 2$ matrices. This then translates into the same result for arbitrary $2^n$-dimensional spaces, as if $\\mathcal V\\equiv\\{v_k\\}$ is a basis for $V$, then the sets of tensor products of elements of $\\mathcal V$ form a basis for $V^{\\otimes n}$.\n",
    "\n",
    "This is a special instance of the more general fact that any matrix can be decomposed uniquely as sum of a Hermitian and a skew-Hermitian matrix, as\n",
    "\n",
    "$$\n",
    "A=\\frac{A+A^\\dagger}{2}+i\\frac{A-A^\\dagger}{2i},\n",
    "$$\n",
    "\n",
    "and the fact that products of Pauli matrices give you a basis for the set of Hermitian matrices.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## Operator Flow\n",
    "\n",
    "Qiskit provides classes representing states and operators and sums, tensor products, and compositions thereof. These algebraic constructs allow us to build expressions representing operators.\n",
    "\n",
    "We introduce expressions by building them from Pauli operators. In subsequent sections we explore in more detail operators and states, how they are represented, and what we can do with them. In the last section we construct a state, evolve it with a Hamiltonian, and compute expectation values of an observable.\n",
    "\n",
    "### Pauli operators, sums, compositions, and tensor products\n",
    "\n",
    "The most important base operators are the Pauli operators.\n",
    "The Pauli operators are represented like this.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "I X Y Z\n"
     ]
    }
   ],
   "source": [
    "from qiskit.opflow import I, X, Y, Z\n",
    "print(I, X, Y, Z)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "These operators may also carry a coefficient."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.5 * I\n",
      "2.5 * X\n"
     ]
    }
   ],
   "source": [
    "print(1.5 * I)\n",
    "print(2.5 * X)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "These coefficients allow the operators to be used as terms in a sum."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.0 * X\n",
      "+ 2.0 * Y\n"
     ]
    }
   ],
   "source": [
    "print(X + 2.0 * Y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "Tensor products are denoted with a caret, like this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "XYZ\n"
     ]
    }
   ],
   "source": [
    "print(X^Y^Z)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "Composition is denoted by the `@` symbol."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "iI\n"
     ]
    }
   ],
   "source": [
    "print(X @ Y @ Z)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "In the preceding two examples, the tensor product and composition of Pauli operators were immediately reduced to the equivalent (possibly multi-qubit) Pauli operator. If we tensor or compose more complicated objects, the result is objects representing the unevaluated operations. That is, algebraic expressions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "For example, composing two sums gives"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1j * Z\n",
      "+ -1j * Y\n",
      "+ 1.0 * I\n",
      "+ 1j * X\n"
     ]
    }
   ],
   "source": [
    "print((X + Y) @ (Y + Z))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "And tensoring two sums gives"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "scrolled": true,
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.0 * XY\n",
      "+ 1.0 * XZ\n",
      "+ 1.0 * YY\n",
      "+ 1.0 * YZ\n"
     ]
    }
   ],
   "source": [
    "print((X + Y) ^ (Y + Z))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "Let's take a closer look at the types introduced above. First the Pauli operators."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": "(PauliOp(Pauli('I'), coeff=1.0), PauliOp(Pauli('X'), coeff=1.0))"
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(I, X)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "Each Pauli operator is an instance of `PauliOp`, which wraps an instance of `qiskit.quantum_info.Pauli`, and adds a coefficient `coeff`. In general, a `PauliOp` represents a weighted tensor product of Pauli operators."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": "PauliOp(Pauli('XYZ'), coeff=2.0)"
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "2.0 * X^Y^Z"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "For the encoding of the Pauli operators as pairs of Boolean values, see the documentation for `qiskit.quantum_info.Pauli`.\n",
    "\n",
    "All of the objects representing operators, whether as \"primitive\"s such as `PauliOp`, or algebraic expressions carry a coefficient"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.2 * (\n",
      "  1.1 * XY\n",
      "  + 1.4300000000000002 * XZ\n",
      ")\n"
     ]
    }
   ],
   "source": [
    "print(1.1 * ((1.2 * X)^(Y + (1.3 * Z))))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "In the following we take a broader and deeper look at Qiskit's operators, states, and the building blocks of quantum algorithms."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "\n",
    "## Part I: State Functions and Measurements\n",
    "\n",
    "Quantum states are represented by subclasses of the class `StateFn`. There are four representations of quantum states: `DictStateFn` is a sparse representation in the computational basis, backed by a `dict`. `VectorStateFn` is a dense representation in the computational basis backed by a numpy array. `CircuitStateFn` is backed by a circuit and represents the state obtained by executing the circuit on the all-zero computational-basis state. `OperatorStateFn` represents mixed states via a density matrix. (As we will see later, `OperatorStateFn` is also used to represent observables.)\n",
    "\n",
    "Several `StateFn` instances are provided for convenience. For example `Zero, One, Plus, Minus`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "from qiskit.opflow import (StateFn, Zero, One, Plus, Minus, H,\n",
    "                           DictStateFn, VectorStateFn, CircuitStateFn, OperatorStateFn)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "`Zero` and `One` represent the quantum states $|0\\rangle$ and $|1\\rangle$. They are represented via `DictStateFn`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "DictStateFn({'0': 1}) DictStateFn({'1': 1})\n"
     ]
    }
   ],
   "source": [
    "print(Zero, One)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "`Plus` and `Minus`, representing states $(|0\\rangle + |1\\rangle)/\\sqrt{2}$ and $(|0\\rangle - |1\\rangle)/\\sqrt{2}$ are represented via circuits. `H` is a synonym for `Plus`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CircuitStateFn(\n",
      "   ┌───┐\n",
      "q: ┤ H ├\n",
      "   └───┘\n",
      ") CircuitStateFn(\n",
      "   ┌───┐┌───┐\n",
      "q: ┤ X ├┤ H ├\n",
      "   └───┘└───┘\n",
      ")\n"
     ]
    }
   ],
   "source": [
    "print(Plus, Minus)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "Indexing into quantum states is done with the `eval` method. These examples return the coefficients of the `0` and `1` basis states. (Below, we will see that the `eval` method is used for other computations, as well.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.0\n",
      "0.0\n",
      "1.0\n",
      "(0.7071067811865475+0j)\n",
      "(-0.7071067811865475+8.7e-17j)\n"
     ]
    }
   ],
   "source": [
    "print(Zero.eval('0'))\n",
    "print(Zero.eval('1'))\n",
    "print(One.eval('1'))\n",
    "print(Plus.eval('0'))\n",
    "print(Minus.eval('1'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "The dual vector of a quantum state, that is the *bra* corresponding to a *ket* is obtained via the `adjoint` method. The `StateFn` carries a flag `is_measurement`, which is `False` if the object is a ket and `True` if it is a bra."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "Here, we construct $\\langle 1 |$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": "DictStateFn({'1': 1}, coeff=1.0, is_measurement=True)"
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "One.adjoint()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "For convenience, one may obtain the dual vector with a tilde, like this"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": "DictStateFn({'1': 1}, coeff=1.0, is_measurement=True)"
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "~One"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "### Algebraic operations and predicates\n",
    "\n",
    "Many algebraic operations and predicates between `StateFn`s are supported, including:\n",
    "\n",
    "* `+` - addition\n",
    "* `-` - subtraction, negation (scalar multiplication by -1)\n",
    "* `*` - scalar multiplication\n",
    "* `/` - scalar division\n",
    "* `@` - composition\n",
    "* `^` - tensor product or tensor power (tensor with self n times)\n",
    "* `**` - composition power (compose with self n times)\n",
    "* `==` - equality\n",
    "* `~` - adjoint, alternating between a State Function and Measurement\n",
    "\n",
    "Be very aware that these operators obey the [Python rules for operator precedence](https://docs.python.org/3/reference/expressions.html#operator-precedence), which might not be what you expect mathematically.  For example, `I^X + X^I` will actually be parsed as `I ^ (X + X) ^ I == 2 * (I^X^I)` because Python evaluates `+` before `^`.  In these cases, you can use the methods (`.tensor()`, etc) or parentheses."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "`StateFn`s carry a coefficient. This allows us to multiply states by a scalar, and so to construct sums."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "Here, we construct $(2 + 3i)|0\\rangle$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": "DictStateFn({'0': 1}, coeff=(2+3j), is_measurement=False)"
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(2.0 + 3.0j) * Zero"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "Here, we see that adding two `DictStateFn`s returns an object of the same type. We construct $|0\\rangle + |1\\rangle$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "DictStateFn({'0': 1.0, '1': 1.0})\n"
     ]
    }
   ],
   "source": [
    "print(Zero + One)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "Note that you must normalize states by hand. For example, to construct $(|0\\rangle + |1\\rangle)/\\sqrt{2}$, we write"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "DictStateFn({'0': 1.0, '1': 1.0}) * 0.7071067811865475\n"
     ]
    }
   ],
   "source": [
    "import math\n",
    "\n",
    "v_zero_one = (Zero + One) / math.sqrt(2)\n",
    "print(v_zero_one)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "In other cases, the result is a symbolic representation of a sum. For example, here is a representation of $|+\\rangle + |-\\rangle$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SummedOp([\n",
      "  CircuitStateFn(\n",
      "     ┌───┐\n",
      "  q: ┤ H ├\n",
      "     └───┘\n",
      "  ),\n",
      "  CircuitStateFn(\n",
      "     ┌───┐┌───┐\n",
      "  q: ┤ X ├┤ H ├\n",
      "     └───┘└───┘\n",
      "  )\n",
      "])\n"
     ]
    }
   ],
   "source": [
    "print(Plus + Minus)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "The composition operator is used to perform an inner product, which by default is held in an unevaluated form. Here is a representation of $\\langle 1 | 1 \\rangle$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ComposedOp([\n",
      "  DictMeasurement({'1': 1}),\n",
      "  DictStateFn({'1': 1})\n",
      "])\n"
     ]
    }
   ],
   "source": [
    "print(~One @ One)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "Note that the `is_measurement` flag causes the (bra) state `~One` to be printed `DictMeasurement`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "Symbolic expressions may be evaluated with the `eval` method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": "1.0"
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(~One @ One).eval()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": "0.9999999999999998"
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(~v_zero_one @ v_zero_one).eval()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "Here is $\\langle - | 1 \\rangle = \\langle (\\langle 0| - \\langle 1|)/\\sqrt{2} | 1\\rangle$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": "(-0.7071067811865475-8.7e-17j)"
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(~Minus @ One).eval()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "The composition operator `@` is equivalent to calling the `compose` method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ComposedOp([\n",
      "  DictMeasurement({'1': 1}),\n",
      "  DictStateFn({'1': 1})\n",
      "])\n"
     ]
    }
   ],
   "source": [
    "print((~One).compose(One))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "Inner products may also be computed using the `eval` method directly, without constructing a `ComposedOp`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": "1.0"
     },
     "execution_count": 26,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(~One).eval(One)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "Symbolic tensor products are constructed as follows. Here is $|0\\rangle \\otimes |+\\rangle$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "TensoredOp([\n",
      "  DictStateFn({'0': 1}),\n",
      "  CircuitStateFn(\n",
      "     ┌───┐\n",
      "  q: ┤ H ├\n",
      "     └───┘\n",
      "  )\n",
      "])\n"
     ]
    }
   ],
   "source": [
    "print(Zero^Plus)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "This may be represented as a simple (not compound) `CircuitStateFn`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {
    "scrolled": false,
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CircuitStateFn(\n",
      "     ┌───┐\n",
      "q_0: ┤ H ├\n",
      "     └───┘\n",
      "q_1: ─────\n",
      "          \n",
      ")\n"
     ]
    }
   ],
   "source": [
    "print((Zero^Plus).to_circuit_op())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "Tensor powers are constructed using the caret `^` as follows. Here are $600 (|11111\\rangle + |00000\\rangle)$, and $|10\\rangle^{\\otimes 3}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "DictStateFn({'11111': 1.0, '00000': 1.0}) * 600.0\n",
      "DictStateFn({'101010': 1})\n"
     ]
    }
   ],
   "source": [
    "print(600 * ((One^5) + (Zero^5)))\n",
    "print((One^Zero)^3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "The method `to_matrix_op` converts to `VectorStateFn`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "VectorStateFn(Statevector([ 0.25-6.1e-17j, -0.25+6.1e-17j,  0.25-6.1e-17j,\n",
      "             -0.25+6.1e-17j, -0.25+6.1e-17j,  0.25-6.1e-17j,\n",
      "             -0.25+6.1e-17j,  0.25-6.1e-17j,  0.25-6.1e-17j,\n",
      "             -0.25+6.1e-17j,  0.25-6.1e-17j, -0.25+6.1e-17j,\n",
      "             -0.25+6.1e-17j,  0.25-6.1e-17j, -0.25+6.1e-17j,\n",
      "              0.25-6.1e-17j],\n",
      "            dims=(2, 2, 2, 2)))\n",
      "CircuitStateFn(\n",
      "     ┌───┐\n",
      "q_0: ┤ X ├\n",
      "     ├───┤\n",
      "q_1: ┤ H ├\n",
      "     ├───┤\n",
      "q_2: ┤ X ├\n",
      "     ├───┤\n",
      "q_3: ┤ H ├\n",
      "     └───┘\n",
      ")\n",
      "{'1101': 0.2529296875, '1111': 0.2529296875, '0111': 0.248046875, '0101': 0.24609375}\n"
     ]
    }
   ],
   "source": [
    "print(((Plus^Minus)^2).to_matrix_op())\n",
    "print(((Plus^One)^2).to_circuit_op())\n",
    "print(((Plus^One)^2).to_matrix_op().sample())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "Constructing a StateFn is easy. The `StateFn` class also serves as a factory, and can take any applicable primitive in its constructor and return the correct StateFn subclass. Currently the following primitives can be passed into the constructor, listed alongside the `StateFn` subclass they produce:\n",
    "\n",
    "* str (equal to some basis bitstring) -> DictStateFn\n",
    "* dict  -> DictStateFn\n",
    "* Qiskit Result object -> DictStateFn\n",
    "* list -> VectorStateFn\n",
    "* np.ndarray -> VectorStateFn\n",
    "* Statevector -> VectorStateFn\n",
    "* QuantumCircuit -> CircuitStateFn\n",
    "* Instruction -> CircuitStateFn\n",
    "* OperatorBase -> OperatorStateFn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "DictStateFn({'0': 1})\n",
      "True\n",
      "VectorStateFn(Statevector([0.+0.j, 1.+0.j, 1.+0.j, 0.+0.j],\n",
      "            dims=(2, 2)))\n",
      "CircuitStateFn(\n",
      "     ┌──────────────────────────────────────────────────────────┐\n",
      "q_0: ┤0                                                         ├\n",
      "     │  RealAmplitudes(θ[0],θ[1],θ[2],θ[3],θ[4],θ[5],θ[6],θ[7]) │\n",
      "q_1: ┤1                                                         ├\n",
      "     └──────────────────────────────────────────────────────────┘\n",
      ")\n"
     ]
    }
   ],
   "source": [
    "print(StateFn({'0':1}))\n",
    "print(StateFn({'0':1}) == Zero)\n",
    "\n",
    "print(StateFn([0,1,1,0]))\n",
    "\n",
    "from qiskit.circuit.library import RealAmplitudes\n",
    "print(StateFn(RealAmplitudes(2)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## Part II: `PrimitiveOp`s\n",
    "\n",
    "The basic Operators are subclasses of `PrimitiveOp`. Just like `StateFn`, `PrimitiveOp` is also a factory for creating the correct type of `PrimitiveOp` for a given primitive. Currently, the following primitives can be passed into the constructor, listed alongside the `PrimitiveOp` subclass they produce:\n",
    "\n",
    "* Terra's Pauli -> PauliOp\n",
    "* Instruction -> CircuitOp\n",
    "* QuantumCircuit -> CircuitOp\n",
    "* 2d List -> MatrixOp\n",
    "* np.ndarray -> MatrixOp\n",
    "* spmatrix -> MatrixOp\n",
    "* Terra's quantum_info.Operator -> MatrixOp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "from qiskit.opflow import X, Y, Z, I, CX, T, H, S, PrimitiveOp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "### Matrix elements"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "The `eval` method returns a column from an operator. For example, the Pauli $X$ operator is represented by a `PauliOp`. Asking for a column returns an instance of the sparse representation, a `DictStateFn`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": "PauliOp(Pauli('X'), coeff=1.0)"
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "DictStateFn({'1': (1+0j)})\n"
     ]
    }
   ],
   "source": [
    "print(X.eval('0'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "It follows that indexing into an operator, that is obtaining a matrix element, is performed with two calls to the `eval` method."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "We have $X = \\left(\\begin{matrix} 0 & 1 \\\\\n",
    "                            1 & 0\n",
    "                            \\end{matrix} \\right)$. And the matrix element $\\left\\{X \\right\\}_{0,1}$ is"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": "(1+0j)"
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X.eval('0').eval('1')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "Here is an example using the two qubit operator `CX`, the controlled `X`, which is represented by a circuit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "          \n",
      "q_0: ──■──\n",
      "     ┌─┴─┐\n",
      "q_1: ┤ X ├\n",
      "     └───┘\n",
      "[[1. 0. 0. 0.]\n",
      " [0. 0. 0. 1.]\n",
      " [0. 0. 1. 0.]\n",
      " [0. 1. 0. 0.]]\n"
     ]
    }
   ],
   "source": [
    "print(CX)\n",
    "print(CX.to_matrix().real) # The imaginary part vanishes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": "VectorStateFn(Statevector([0.+0.j, 0.+0.j, 0.+0.j, 1.+0.j],\n            dims=(2, 2)), coeff=1.0, is_measurement=False)"
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "CX.eval('01')  # 01 is the one in decimal. We get the first column."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": "(1+0j)"
     },
     "execution_count": 38,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "CX.eval('01').eval('11')  # This returns element with (zero-based) index (1, 3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "### Applying an operator to a state vector"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "Applying an operator to a state vector may be done with the `compose` method (equivalently, `@` operator). Here is a representation of $X | 1 \\rangle = |0\\rangle$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ComposedOp([\n",
      "  X,\n",
      "  DictStateFn({'1': 1})\n",
      "])\n"
     ]
    }
   ],
   "source": [
    "print(X @ One)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "A simpler representation, the `DictStateFn` representation of $|0\\rangle$, is obtained with `eval`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": "DictStateFn({'0': (1+0j)}, coeff=(1+0j), is_measurement=False)"
     },
     "execution_count": 40,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(X @ One).eval()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "The intermediate `ComposedOp` step may be avoided by using `eval` directly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": "DictStateFn({'0': (1+0j)}, coeff=(1+0j), is_measurement=False)"
     },
     "execution_count": 41,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X.eval(One)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "Composition and tensor products of operators are effected with `@` and `^`. Here are some examples."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(1+0j)\n",
      "          ┌───┐          ┌───┐     \n",
      "q_0: ──■──┤ H ├───────■──┤ H ├─────\n",
      "     ┌─┴─┐└───┘┌───┐┌─┴─┐└───┘┌───┐\n",
      "q_1: ┤ X ├──■──┤ H ├┤ X ├──■──┤ H ├\n",
      "     └───┘┌─┴─┐├───┤└───┘┌─┴─┐├───┤\n",
      "q_2: ──■──┤ X ├┤ H ├──■──┤ X ├┤ H ├\n",
      "     ┌─┴─┐└───┘├───┤┌─┴─┐└───┘├───┤\n",
      "q_3: ┤ X ├──■──┤ H ├┤ X ├──■──┤ H ├\n",
      "     └───┘┌─┴─┐├───┤└───┘┌─┴─┐├───┤\n",
      "q_4: ─────┤ X ├┤ H ├─────┤ X ├┤ H ├\n",
      "          └───┘└───┘     └───┘└───┘\n",
      "CircuitStateFn(\n",
      "     ┌───┐┌───┐     ┌───┐          ┌───┐     \n",
      "q_0: ┤ X ├┤ H ├──■──┤ H ├───────■──┤ H ├─────\n",
      "     ├───┤├───┤┌─┴─┐└───┘┌───┐┌─┴─┐└───┘┌───┐\n",
      "q_1: ┤ X ├┤ H ├┤ X ├──■──┤ H ├┤ X ├──■──┤ H ├\n",
      "     ├───┤├───┤└───┘┌─┴─┐├───┤└───┘┌─┴─┐├───┤\n",
      "q_2: ┤ X ├┤ H ├──■──┤ X ├┤ H ├──■──┤ X ├┤ H ├\n",
      "     ├───┤├───┤┌─┴─┐└───┘├───┤┌─┴─┐└───┘├───┤\n",
      "q_3: ┤ X ├┤ H ├┤ X ├──■──┤ H ├┤ X ├──■──┤ H ├\n",
      "     ├───┤├───┤└───┘┌─┴─┐├───┤└───┘┌─┴─┐├───┤\n",
      "q_4: ┤ X ├┤ H ├─────┤ X ├┤ H ├─────┤ X ├┤ H ├\n",
      "     └───┘└───┘     └───┘└───┘     └───┘└───┘\n",
      ")\n",
      "CircuitStateFn(\n",
      "     ┌─────────────┐     \n",
      "q_0: ┤0            ├─────\n",
      "     │             │     \n",
      "q_1: ┤1 Pauli(XII) ├─────\n",
      "     │             │┌───┐\n",
      "q_2: ┤2            ├┤ H ├\n",
      "     └─────────────┘└───┘\n",
      ")\n"
     ]
    }
   ],
   "source": [
    "print(((~One^2) @ (CX.eval('01'))).eval())\n",
    "\n",
    "print(((H^5) @ ((CX^2)^I) @ (I^(CX^2)))**2)\n",
    "print((((H^5) @ ((CX^2)^I) @ (I^(CX^2)))**2) @ (Minus^5))\n",
    "print(((H^I^I)@(X^I^I)@Zero))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ComposedOp([\n",
      "  DictMeasurement({'1': 1}),\n",
      "  CircuitStateFn(\n",
      "     ┌───┐┌───┐\n",
      "  q: ┤ X ├┤ H ├\n",
      "     └───┘└───┘\n",
      "  )\n",
      "])\n"
     ]
    }
   ],
   "source": [
    "print(~One @ Minus)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## Part III `ListOp` and subclasses"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "### `ListOp`\n",
    "\n",
    "`ListOp` is a container for effectively vectorizing operations over a list of operators and states."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ComposedOp([\n",
      "  ListOp([\n",
      "    DictMeasurement({'1': 1}),\n",
      "    DictMeasurement({'0': 1})\n",
      "  ]),\n",
      "  ListOp([\n",
      "    DictStateFn({'1': 1}),\n",
      "    DictStateFn({'0': 1})\n",
      "  ])\n",
      "])\n"
     ]
    }
   ],
   "source": [
    "from qiskit.opflow import ListOp\n",
    "\n",
    "print((~ListOp([One, Zero]) @ ListOp([One, Zero])))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "For example, the composition above is distributed over the lists (`ListOp`) using the simplification method `reduce`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ListOp([\n",
      "  ListOp([\n",
      "    ComposedOp([\n",
      "      DictMeasurement({'1': 1}),\n",
      "      DictStateFn({'1': 1})\n",
      "    ]),\n",
      "    ComposedOp([\n",
      "      DictMeasurement({'1': 1}),\n",
      "      DictStateFn({'0': 1})\n",
      "    ])\n",
      "  ]),\n",
      "  ListOp([\n",
      "    ComposedOp([\n",
      "      DictMeasurement({'0': 1}),\n",
      "      DictStateFn({'1': 1})\n",
      "    ]),\n",
      "    ComposedOp([\n",
      "      DictMeasurement({'0': 1}),\n",
      "      DictStateFn({'0': 1})\n",
      "    ])\n",
      "  ])\n",
      "])\n"
     ]
    }
   ],
   "source": [
    "print((~ListOp([One, Zero]) @ ListOp([One, Zero])).reduce())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "### `ListOp`s: `SummedOp`, `ComposedOp`, `TensoredOp`\n",
    "\n",
    "`ListOp`, introduced above, is useful for vectorizing operations. But, it also serves as the superclass for list-like composite classes.\n",
    "If you've already played around with the above, you'll notice that you can easily perform operations between `OperatorBase`s which we may not know how to perform efficiently in general (or simply haven't implemented an efficient procedure for yet), such as addition between `CircuitOp`s. In those cases, you may receive a `ListOp` result (or subclass thereof) from your operation representing the lazy execution of the operation. For example, if you attempt to add together a `DictStateFn` and a `CircuitStateFn`, you'll receive a `SummedOp` representing the sum of the two. This composite State function still has a working `eval` (but may need to perform a non-scalable computation under the hood, such as converting both to vectors).\n",
    "\n",
    "These composite `OperatorBase`s are how we construct increasingly complex and rich computation out of `PrimitiveOp` and `StateFn` building blocks.\n",
    "\n",
    "Every `ListOp` has four properties:\n",
    "\n",
    "* `oplist` - The list of `OperatorBase`s which may represent terms, factors, etc.\n",
    "* `combo_fn` - The function taking a list of complex numbers to an output value which defines how to combine the outputs of the `oplist` items. For broadcasting simplicity, this function is defined over NumPy arrays.\n",
    "* `coeff` - A coefficient multiplying the primitive. Note that `coeff` can be int, float, complex or a free `Parameter` object (from `qiskit.circuit` in Terra) to be bound later using `my_op.bind_parameters`.\n",
    "* `abelian` - Indicates whether the Operators in `oplist` are known to mutually commute (usually set after being converted by the `AbelianGrouper` converter).\n",
    "\n",
    "Note that `ListOp` supports typical sequence overloads, so you can use indexing like `my_op[4]` to access the `OperatorBase`s in `oplist`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "### `OperatorStateFn`\n",
    "\n",
    "We mentioned above that `OperatorStateFn` represents a density operator. But, if the `is_measurement` flag is `True`, then `OperatorStateFn` represents an observable. The expectation value of this observable can then be constructed via `ComposedOp`. Or, directly, using `eval`. Recall that the `is_measurement` flag (property) is set via the `adjoint` method."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "Here we construct the observable corresponding to the Pauli $Z$ operator. Note that when printing, it is called `OperatorMeasurement`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "OperatorMeasurement(Z)\n"
     ]
    },
    {
     "data": {
      "text/plain": "OperatorStateFn(PauliOp(Pauli('Z'), coeff=1.0), coeff=1.0, is_measurement=True)"
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "print(StateFn(Z).adjoint())\n",
    "StateFn(Z).adjoint()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "Here, we compute $\\langle 0 | Z | 0 \\rangle$, $\\langle 1 | Z | 1 \\rangle$, and $\\langle + | Z | + \\rangle$, where $|+\\rangle = (|0\\rangle + |1\\rangle)/\\sqrt{2}$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(1+0j)\n",
      "(-1+0j)\n",
      "0j\n"
     ]
    }
   ],
   "source": [
    "print(StateFn(Z).adjoint().eval(Zero))\n",
    "print(StateFn(Z).adjoint().eval(One))\n",
    "print(StateFn(Z).adjoint().eval(Plus))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "## Part IV: Converters\n",
    "\n",
    "Converters are classes that manipulate operators and states and perform building blocks of algorithms. Examples include changing the basis of operators and Trotterization.\n",
    "Converters traverse an expression and perform a particular manipulation or replacement, defined by the converter's `convert()` method, of the Operators within. Typically, if a converter encounters an `OperatorBase` in the recursion which is irrelevant to its conversion purpose, that `OperatorBase` is left unchanged."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from qiskit.opflow import I, X, Y, Z, H, CX, Zero, ListOp, PauliExpectation, PauliTrotterEvolution, CircuitSampler, MatrixEvolution, Suzuki\n",
    "from qiskit.circuit import Parameter\n",
    "from qiskit import Aer"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "### Evolutions, `exp_i()`, and the `EvolvedOp`\n",
    "\n",
    "Every `PrimitiveOp` and `ListOp` has an `.exp_i()` function such that `H.exp_i()` corresponds to $e^{-iH}$. In practice, only a few of these Operators have an efficiently computable exponentiation (such as MatrixOp and the PauliOps with only one non-identity single-qubit Pauli), so we need to return a placeholder, or symbolic representation, (similar to how `SummedOp` is a placeholder when we can't perform addition). This placeholder is called `EvolvedOp`, and it holds the `OperatorBase` to be exponentiated in its `.primitive` property.\n",
    "\n",
    "Qiskit operators fully support parameterization, so we can use a `Parameter` for our evolution time here. Notice that there's no \"evolution time\" argument in any function. The Operator flow exponentiates whatever operator we tell it to, and if we choose to multiply the operator by an evolution time, $e^{-iHt}$, that will be reflected in our exponentiation parameters."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "#### Weighted sum of Pauli operators\n",
    "A Hamiltonian expressed as a linear combination of multi-qubit Pauli operators may be constructed like this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "two_qubit_H2 =  (-1.0523732 * I^I) + \\\n",
    "                (0.39793742 * I^Z) + \\\n",
    "                (-0.3979374 * Z^I) + \\\n",
    "                (-0.0112801 * Z^Z) + \\\n",
    "                (0.18093119 * X^X)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "Note that `two_qubit_H2` is represented as a `SummedOp` whose terms are `PauliOp`s."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-1.0523732 * II\n",
      "+ 0.39793742 * IZ\n",
      "- 0.3979374 * ZI\n",
      "- 0.0112801 * ZZ\n",
      "+ 0.18093119 * XX\n"
     ]
    }
   ],
   "source": [
    "print(two_qubit_H2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "Next, we multiply the Hamiltonian by a `Parameter`. This `Parameter` is stored in the `coeff` property of the `SummedOp`. Calling `exp_i()` on the result wraps it in `EvolvedOp`, representing exponentiation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "e^(-i*1.0*θ * (\n",
      "  -1.0523732 * II\n",
      "  + 0.39793742 * IZ\n",
      "  - 0.3979374 * ZI\n",
      "  - 0.0112801 * ZZ\n",
      "  + 0.18093119 * XX\n",
      "))\n",
      "EvolvedOp(PauliSumOp(SparsePauliOp(['II', 'IZ', 'ZI', 'ZZ', 'XX'],\n",
      "              coeffs=[-1.0523732 +0.j,  0.39793742+0.j, -0.3979374 +0.j, -0.0112801 +0.j,\n",
      "  0.18093119+0.j]), coeff=1.0*θ), coeff=1.0)\n"
     ]
    }
   ],
   "source": [
    "evo_time = Parameter('θ')\n",
    "evolution_op = (evo_time*two_qubit_H2).exp_i()\n",
    "print(evolution_op) # Note, EvolvedOps print as exponentiations\n",
    "print(repr(evolution_op))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "We construct `h2_measurement`, which represents `two_qubit_H2` as an observable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "OperatorMeasurement(-1.0523732 * II\n",
      "+ 0.39793742 * IZ\n",
      "- 0.3979374 * ZI\n",
      "- 0.0112801 * ZZ\n",
      "+ 0.18093119 * XX)\n"
     ]
    }
   ],
   "source": [
    "h2_measurement = StateFn(two_qubit_H2).adjoint()\n",
    "print(h2_measurement)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "We construct a Bell state $|\\Phi_+\\rangle$ via $\\text{CX} (H\\otimes I) |00\\rangle$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CircuitStateFn(\n",
      "     ┌───┐     \n",
      "q_0: ┤ H ├──■──\n",
      "     └───┘┌─┴─┐\n",
      "q_1: ─────┤ X ├\n",
      "          └───┘\n",
      ")\n"
     ]
    }
   ],
   "source": [
    "bell = CX @ (I ^ H) @ Zero\n",
    "print(bell)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "Here is the expression $H e^{-iHt} |\\Phi_+\\rangle$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ComposedOp([\n",
      "  OperatorMeasurement(-1.0523732 * II\n",
      "  + 0.39793742 * IZ\n",
      "  - 0.3979374 * ZI\n",
      "  - 0.0112801 * ZZ\n",
      "  + 0.18093119 * XX),\n",
      "  e^(-i*1.0*θ * (\n",
      "    -1.0523732 * II\n",
      "    + 0.39793742 * IZ\n",
      "    - 0.3979374 * ZI\n",
      "    - 0.0112801 * ZZ\n",
      "    + 0.18093119 * XX\n",
      "  )),\n",
      "  CircuitStateFn(\n",
      "       ┌───┐     \n",
      "  q_0: ┤ H ├──■──\n",
      "       └───┘┌─┴─┐\n",
      "  q_1: ─────┤ X ├\n",
      "            └───┘\n",
      "  )\n",
      "])\n"
     ]
    }
   ],
   "source": [
    "evo_and_meas = h2_measurement @ evolution_op @ bell\n",
    "print(evo_and_meas)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "Typically, we want to approximate $e^{-iHt}$ using two-qubit gates. We achieve this with the `convert` method of `PauliTrotterEvolution`, which traverses expressions applying trotterization to all `EvolvedOp`s encountered. Although we use `PauliTrotterEvolution` here, there are other possibilities, such as `MatrixEvolution`, which performs the exponentiation exactly."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ComposedOp([\n",
      "  OperatorMeasurement(-1.0523732 * II\n",
      "  + 0.39793742 * IZ\n",
      "  - 0.3979374 * ZI\n",
      "  - 0.0112801 * ZZ\n",
      "  + 0.18093119 * XX),\n",
      "  CircuitStateFn(\n",
      "       ┌───┐     ┌────────────────────────────────────────┐\n",
      "  q_0: ┤ H ├──■──┤0                                       ├\n",
      "       └───┘┌─┴─┐│  exp(-i II + IZ + ZI + ZZ + XX)(1.0*θ) │\n",
      "  q_1: ─────┤ X ├┤1                                       ├\n",
      "            └───┘└────────────────────────────────────────┘\n",
      "  )\n",
      "])\n"
     ]
    }
   ],
   "source": [
    "trotterized_op = PauliTrotterEvolution(trotter_mode=Suzuki(order=2, reps=1)).convert(evo_and_meas)\n",
    "# We can also set trotter_mode='suzuki' or leave it empty to default to first order Trotterization.\n",
    "print(trotterized_op)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "`trotterized_op` contains a `Parameter`. The `bind_parameters` method traverses the expression binding values to parameter names as specified via a `dict`. In this case, there is only one parameter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "bound = trotterized_op.bind_parameters({evo_time: .5})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "`bound` is a `ComposedOp`. The second factor is the circuit. Let's draw it to verify that the binding has taken place."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": "     ┌───┐     ┌──────────────────────────────────────┐\nq_0: ┤ H ├──■──┤0                                     ├\n     └───┘┌─┴─┐│  exp(-i II + IZ + ZI + ZZ + XX)(0.5) │\nq_1: ─────┤ X ├┤1                                     ├\n          └───┘└──────────────────────────────────────┘",
      "text/html": "<pre style=\"word-wrap: normal;white-space: pre;background: #fff0;line-height: 1.1;font-family: &quot;Courier New&quot;,Courier,monospace\">     ┌───┐     ┌──────────────────────────────────────┐\nq_0: ┤ H ├──■──┤0                                     ├\n     └───┘┌─┴─┐│  exp(-i II + IZ + ZI + ZZ + XX)(0.5) │\nq_1: ─────┤ X ├┤1                                     ├\n          └───┘└──────────────────────────────────────┘</pre>"
     },
     "execution_count": 57,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "bound[1].to_circuit().draw()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "### Expectations\n",
    "\n",
    "`Expectation`s are converters that enable the computation of expectation values of observables. They traverse an Operator tree, replacing `OperatorStateFn`s (observables) with equivalent instructions which are more amenable to\n",
    "computation on quantum or classical hardware. For example, if we want to measure the expectation value of an Operator `o` expressed as a sum of Paulis with respect to some state function, but can only access diagonal measurements on quantum hardware, we can create an observable `~StateFn(o)` and use a ``PauliExpectation`` to convert it to a diagonal measurement and circuit pre-rotations to append to the state.\n",
    "\n",
    "Another interesting `Expectation` is the `AerPauliExpectation`, which converts the observable into a `CircuitStateFn` containing a special expectation snapshot instruction which `Aer` can execute natively with high performance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SummedOp([\n",
      "  ComposedOp([\n",
      "    OperatorMeasurement(-1.0523732 * II),\n",
      "    II\n",
      "  ]),\n",
      "  ComposedOp([\n",
      "    OperatorMeasurement(0.39793742 * IZ),\n",
      "    II\n",
      "  ]),\n",
      "  ComposedOp([\n",
      "    OperatorMeasurement(-0.3979374 * ZI),\n",
      "    II\n",
      "  ]),\n",
      "  ComposedOp([\n",
      "    OperatorMeasurement(-0.0112801 * ZZ),\n",
      "    II\n",
      "  ]),\n",
      "  ComposedOp([\n",
      "    OperatorMeasurement(0.18093119 * ZZ),\n",
      "         ┌───┐\n",
      "    q_0: ┤ H ├\n",
      "         ├───┤\n",
      "    q_1: ┤ H ├\n",
      "         └───┘\n",
      "  ])\n",
      "])\n"
     ]
    }
   ],
   "source": [
    "# Note that XX was the only non-diagonal measurement in our H2 Observable\n",
    "print(PauliExpectation(group_paulis=False).convert(h2_measurement))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "By default `group_paulis=True`, which will use the `AbelianGrouper` to convert the `SummedOp` into groups of mutually qubit-wise commuting Paulis. This reduces circuit execution overhead, as each group can share the same circuit execution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SummedOp([\n",
      "  ComposedOp([\n",
      "    OperatorMeasurement(0.18093119 * ZZ\n",
      "    - 1.0523732 * II),\n",
      "         ┌───┐\n",
      "    q_0: ┤ H ├\n",
      "         ├───┤\n",
      "    q_1: ┤ H ├\n",
      "         └───┘\n",
      "  ]),\n",
      "  ComposedOp([\n",
      "    OperatorMeasurement(0.39793742 * IZ\n",
      "    - 0.3979374 * ZI\n",
      "    - 0.0112801 * ZZ),\n",
      "    II\n",
      "  ])\n",
      "])\n"
     ]
    }
   ],
   "source": [
    "print(PauliExpectation().convert(h2_measurement))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "Note that converters act recursively, that is, they traverse an expression applying their action only where possible. So we can just convert our full evolution and measurement expression. We could have equivalently composed the converted `h2_measurement` with our evolution `CircuitStateFn`. We proceed by applying the conversion on the entire expression."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {
    "scrolled": true,
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SummedOp([\n",
      "  ComposedOp([\n",
      "    OperatorMeasurement(0.18093119 * ZZ\n",
      "    - 1.0523732 * II),\n",
      "    CircuitStateFn(\n",
      "         ┌───┐     ┌────────────────────────────────────────┐┌───┐\n",
      "    q_0: ┤ H ├──■──┤0                                       ├┤ H ├\n",
      "         └───┘┌─┴─┐│  exp(-i II + IZ + ZI + ZZ + XX)(1.0*θ) │├───┤\n",
      "    q_1: ─────┤ X ├┤1                                       ├┤ H ├\n",
      "              └───┘└────────────────────────────────────────┘└───┘\n",
      "    )\n",
      "  ]),\n",
      "  ComposedOp([\n",
      "    OperatorMeasurement(0.39793742 * IZ\n",
      "    - 0.3979374 * ZI\n",
      "    - 0.0112801 * ZZ),\n",
      "    CircuitStateFn(\n",
      "         ┌───┐     ┌────────────────────────────────────────┐\n",
      "    q_0: ┤ H ├──■──┤0                                       ├\n",
      "         └───┘┌─┴─┐│  exp(-i II + IZ + ZI + ZZ + XX)(1.0*θ) │\n",
      "    q_1: ─────┤ X ├┤1                                       ├\n",
      "              └───┘└────────────────────────────────────────┘\n",
      "    )\n",
      "  ])\n",
      "])\n"
     ]
    }
   ],
   "source": [
    "diagonalized_meas_op = PauliExpectation().convert(trotterized_op)\n",
    "print(diagonalized_meas_op)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "Now we bind multiple parameter values into a `ListOp`, followed by `eval` to evaluate the entire expression. We could have used `eval` earlier if we bound earlier, but it would not be efficient. Here, `eval` will convert our `CircuitStateFn`s to `VectorStateFn`s through simulation internally."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [],
   "source": [
    "evo_time_points = list(range(8))\n",
    "h2_trotter_expectations = diagonalized_meas_op.bind_parameters({evo_time: evo_time_points})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "Here are the expectation values $\\langle \\Phi_+| e^{iHt} H e^{-iHt} |\\Phi_+\\rangle$ corresponding to the different values of the parameter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "data": {
      "text/plain": "array([-0.88272211+0.0e+00j, -0.88272211+0.0e+00j, -0.88272211+0.0e+00j,\n       -0.88272211+0.0e+00j, -0.88272211+0.0e+00j, -0.88272211+0.0e+00j,\n       -0.88272211+0.0e+00j, -0.88272211+5.6e-17j])"
     },
     "execution_count": 62,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "h2_trotter_expectations.eval()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "### Executing `CircuitStateFn`s with the `CircuitSampler`\n",
    "\n",
    "The `CircuitSampler` traverses an Operator and converts any `CircuitStateFn`s into approximations of the resulting state function by a `DictStateFn` or `VectorStateFn` using a quantum backend. Note that in order to approximate the value of the `CircuitStateFn`, it must 1) send the state function through a depolarizing channel, which will destroy all phase information and 2) replace the sampled frequencies with **square roots** of the frequency, rather than the raw probability of sampling (which would be the equivalent of sampling the **square** of the state function, per the Born rule)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sampled Trotterized energies:\n",
      " [-0.88272211 -0.88272211 -0.88272211 -0.88272211 -0.88272211 -0.88272211\n",
      " -0.88272211 -0.88272211]\n"
     ]
    }
   ],
   "source": [
    "sampler = CircuitSampler(backend=Aer.get_backend('aer_simulator'))\n",
    "# sampler.quantum_instance.run_config.shots = 1000\n",
    "sampled_trotter_exp_op = sampler.convert(h2_trotter_expectations)\n",
    "sampled_trotter_energies = sampled_trotter_exp_op.eval()\n",
    "print('Sampled Trotterized energies:\\n {}'.format(np.real(sampled_trotter_energies)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "pycharm": {
     "name": "#%% md\n"
    }
   },
   "source": [
    "Note again that the circuits are replaced by dicts with ***square roots*** of the circuit sampling probabilities. Take a look at one sub-expression before and after the conversion:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {
    "pycharm": {
     "name": "#%%\n"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Before:\n",
      "\n",
      "ComposedOp([\n",
      "  OperatorMeasurement(0.18093119 * ZZ\n",
      "  - 1.0523732 * II),\n",
      "  CircuitStateFn(\n",
      "       ┌───┐     ┌────────────────────────────────────┐┌───┐\n",
      "  q_0: ┤ H ├──■──┤0                                   ├┤ H ├\n",
      "       └───┘┌─┴─┐│  exp(-i II + IZ + ZI + ZZ + XX)(0) │├───┤\n",
      "  q_1: ─────┤ X ├┤1                                   ├┤ H ├\n",
      "            └───┘└────────────────────────────────────┘└───┘\n",
      "  )\n",
      "])\n",
      "\n",
      "After:\n",
      "\n",
      "ComposedOp([\n",
      "  OperatorMeasurement(0.18093119 * ZZ\n",
      "  - 1.0523732 * II),\n",
      "  DictStateFn({'11': 0.7126096406869612, '00': 0.701560760020114})\n",
      "])\n"
     ]
    }
   ],
   "source": [
    "print('Before:\\n')\n",
    "print(h2_trotter_expectations.reduce()[0][0])\n",
    "print('\\nAfter:\\n')\n",
    "print(sampled_trotter_exp_op[0][0])"
   ]
  }
 ],
 "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.9.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}