examples/notebooks/3_Expectation_value_of_observables.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "fc24b8dc",
   "metadata": {},
   "source": [
    "# Computing expectation values of observables\n",
    "\n",
    "In addition of sampling the circuit in a basis or extracting the state-vector in\n",
    "the ideal case, MPQP also allows you to compute expectation values of observables\n",
    "with respect to the state generated by the circuit.\n",
    "\n",
    "To demonstrate this, let us create a circuit generating a 2-qubit state."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "a3c9a6e8",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "     ┌───┐      ┌─────────┐ \n",
      "q_0: ┤ H ├──■───┤ Ry(2.6) ├─\n",
      "     └───┘┌─┴─┐┌┴─────────┴┐\n",
      "q_1: ─────┤ X ├┤ Ry(-0.87) ├\n",
      "          └───┘└───────────┘\n"
     ]
    }
   ],
   "source": [
    "from mpqp import QCircuit\n",
    "from mpqp.gates import *\n",
    "circuit = QCircuit([H(0), CNOT(0,1), Ry(2.6, 0), Ry(-0.87, 1)])\n",
    "print(circuit)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "27afecad",
   "metadata": {},
   "source": [
    "We then import the two needed objects: ``Observable`` and ``ExpectationMeasure``. The first takes as parameter a hermitian matrix (as a numpy array), and the second is the measure to be added to the circuit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "d97f41ba",
   "metadata": {},
   "outputs": [],
   "source": [
    "from mpqp.measures import Observable, ExpectationMeasure"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "df526c3f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import numpy as np\n",
    "from mpqp.tools.maths import is_hermitian\n",
    "\n",
    "matrix = np.array([[4,  2,  3, 8],\n",
    "                   [2, -3,  1, 0],\n",
    "                   [3,  1, -1, 5],\n",
    "                   [8,  0,  5, 2]])\n",
    "is_hermitian(matrix)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "444d0fb2",
   "metadata": {},
   "outputs": [],
   "source": [
    "obs = Observable(matrix)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "61271501",
   "metadata": {},
   "source": [
    "The ``ExpectationMeasure`` takes as parameter the list of qubits corresponding to the observable. The indices of qubits can be given non-ordered, non-contigous and restricted to some qubits, and MPQP will automatically operate on the circuit and the observable to handle that. One has also to precise the number of shots, when sampling the circuit to compute the expectation. If the number of shots is zero (by default), the exact value is returned."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "b6730e48",
   "metadata": {},
   "outputs": [],
   "source": [
    "circuit.add(ExpectationMeasure(obs))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "b1cb8ad4",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "BatchResult: 2 results\n",
      "Result: circuit 1, ATOSDevice, MYQLM_PYLINALG\n",
      " Expectation value: -0.12200000000000011\n",
      " Error/Variance: 0.09998468351171062\n",
      "Result: circuit 1, IBMDevice, AER_SIMULATOR\n",
      " Expectation value: 1.9081958235744878e-17\n",
      " Error/Variance: 0.0\n",
      "expectation_value: -0.12200000000000011\n",
      "expectation_value: 1.9081958235744878e-17\n"
     ]
    }
   ],
   "source": [
    "from mpqp.execution import run, ATOSDevice, IBMDevice\n",
    "results = run(circuit, [ATOSDevice.MYQLM_PYLINALG, IBMDevice.AER_SIMULATOR])\n",
    "print(results)\n",
    "for result in results:\n",
    "    print(f\"expectation_value:\", result.expectation_value)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "2a099233",
   "metadata": {},
   "outputs": [],
   "source": [
    "circuit = circuit.without_measurements()\n",
    "circuit.add(ExpectationMeasure(obs, shots=2000))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "148e5900",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "BatchResult: 2 results\n",
      "Result: circuit 1, ATOSDevice, MYQLM_PYLINALG\n",
      " Expectation value: -3.722\n",
      " Error/Variance: 0.10571682351161849\n",
      "Result: circuit 1, IBMDevice, AER_SIMULATOR\n",
      " Expectation value: -3.5935083233096723\n",
      " Error/Variance: 0.0\n"
     ]
    }
   ],
   "source": [
    "results = run(circuit, [ATOSDevice.MYQLM_PYLINALG, IBMDevice.AER_SIMULATOR])\n",
    "print(results)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "057b9856",
   "metadata": {},
   "source": [
    "## Pauli String representation\n",
    "\n",
    "An observable can also be represented by a Pauli string. This section will demonstrate how to create and manipulate PauliString objects."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "4ac9c7e0",
   "metadata": {},
   "outputs": [],
   "source": [
    "from mpqp.measures import I, X, Y, Z"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aed703a9",
   "metadata": {},
   "source": [
    "> ⚠ **pauli_string import**: pauli atoms are named `I`, `X`, `Y`, and `Z`. If\n",
    "> you have conflicts with the gates of the same name, you can:\n",
    "> \n",
    "> - **Rename Import:**\n",
    "> ```python \n",
    "> from mpqp.measures import X as Pauli_X,  Y as Pauli_Y\n",
    "> ps = Pauli_X + Pauli_Y/2\n",
    "> ```\n",
    "> \n",
    "> - **Import the whole module:**\n",
    "> ```python\n",
    "> from mpqp.measures import pauli_string\n",
    "> ps = pauli_string.X + pauli_string.Y/2\n",
    "> ```"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "30d41d0f",
   "metadata": {},
   "source": [
    "A `PauliString` is based on the following hierarchy:\n",
    "- an *atom* is the most elemental building brick of pauli string, it it either\n",
    "  `I`, `X`, `Y` or `Z`;\n",
    "- a *monomial* is a tensor product of *atoms* multiplied by a real coefficient,\n",
    "  for instance `0.3 * I⊗Z⊗Y`. In MPQP, the tensor product is denoted as `@`, so\n",
    "  the previous example would be expressed as `0.3 * I@Z@Y`;\n",
    "- a *string* is a sum of *monomials*, for instance `0.3 * I@Z@Y + X@X@X`.\n",
    "\n",
    "In practice, you never need to handle these types independently, you can express\n",
    "everything in term of expression based on the fundamental atoms. Let's see how."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7d07ec4b",
   "metadata": {},
   "source": [
    "### Creating and Manipulating PauliString Objects\n",
    "\n",
    "Let's illustrate how to create and manipulate PauliString objects:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "bc5ed517",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ps_1=1*I@Z - 3*X@Y\n"
     ]
    }
   ],
   "source": [
    "ps_1 = I@Z - 3 * X@Y\n",
    "print(f\"{ps_1=}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "937f1a84",
   "metadata": {},
   "source": [
    "#### Using Symbolic Variables in a PauliString\n",
    "\n",
    "We can also define `PauliString` objects with symbolic coefficients using `sympy`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "032ff198",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(θ)*I@X + (k)*Z@Y\n"
     ]
    }
   ],
   "source": [
    "from sympy import symbols\n",
    "\n",
    "theta, k = symbols(\"θ k\")\n",
    "ps_symbol = theta * I @ X + k * Z @ Y\n",
    "print(ps_symbol)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d153fda4",
   "metadata": {},
   "source": [
    "We can substitute numerical values for the symbolic coefficients:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "5dd3c74b",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3.1416*I@X - 1*Z@Y\n"
     ]
    }
   ],
   "source": [
    "ps_symbol = ps_symbol.subs({theta: np.pi, k: -1})\n",
    "print(ps_symbol)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7824baf5",
   "metadata": {},
   "source": [
    "#### Simplifying and Rounding PauliStrings\n",
    "\n",
    "When dealing with PauliStrings, simplification and rounding are common operations:\n",
    "\n",
    "- **Simplify**: We combine like terms and eliminate those with zero coefficients;\n",
    "- **Round**: We can round the coefficients to a specified number of decimals (5 by default).\n",
    "\n",
    "> `str` on a `PauliString` will call both methods: `round()` and `simplify()`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "06784ef3",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ps_2 = 1*I@Z + 2.555555555*Y@I + 1*X@Z - 1*X@Z\n",
      "     = 1*I@Z + 2.555555555*Y@I\n",
      "    ~= 1*I@Z + 2.6*Y@I + 1*X@Z - 1*X@Z\n",
      "    ~= 1*I@Z + 2.5556*Y@I\n"
     ]
    }
   ],
   "source": [
    "ps_2 = I@Z + 2.555555555*Y@I + X@Z - X@Z\n",
    "print(\"ps_2 =\",repr(ps_2))\n",
    "print(\"     =\",repr(ps_2.simplify()))\n",
    "print(\"    ~=\",repr(ps_2.round(1)))\n",
    "print(\"    ~=\",ps_2)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bb3a48ec",
   "metadata": {},
   "source": [
    "#### Arithmetic Operations\n",
    "\n",
    "We can perform various arithmetic operations on PauliString objects, including addition, subtraction, scalar multiplication, scalar division, and matrix multiplication:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "b9cadb08",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Addition:\n",
      "(1*I@Z - 3*X@Y) + (1*I@Z + 2.6*Y@I) = 2*I@Z - 3*X@Y + 2.6*Y@I\n",
      "\n",
      "Subtraction:\n",
      "(1*I@Z - 3*X@Y) - (1*I@Z + 2.6*Y@I) = -3*X@Y - 2.6*Y@I\n",
      "\n",
      "Scalar product:\n",
      "2 * (1*I@Z - 3*X@Y) = 2*I@Z - 6*X@Y\n",
      "\n",
      "Scalar division:\n",
      "(1*I@Z + 2.6*Y@I) / 3 ~= 0.3333*I@Z + 0.8667*Y@I\n",
      "\n",
      "Tensor product:\n",
      "(1*I@Z - 3*X@Y) @ Z = 1*I@Z@Z - 3*X@Y@Z\n",
      "(1*I@Z - 3*X@Y) @ (1*I@Z + 2.6*Y@I) = 1*I@Z@I@Z + 2.6*I@Z@Y@I - 3*X@Y@I@Z - 7.8*X@Y@Y@I\n"
     ]
    }
   ],
   "source": [
    "ps_2 = ps_2.round(1).simplify()\n",
    "print(f\"\"\"Addition:\n",
    "({ps_1}) + ({ps_2}) = {ps_1 + ps_2}\n",
    "\n",
    "Subtraction:\n",
    "({ps_1}) - ({ps_2}) = {ps_1 - ps_2}\n",
    "\n",
    "Scalar product:\n",
    "2 * ({ps_1}) = {2 * ps_1}\n",
    "\n",
    "Scalar division:\n",
    "({ps_2}) / 3 ~= {ps_2 / 3}\n",
    "\n",
    "Tensor product:\n",
    "({ps_1}) @ Z = {ps_1@Z}\n",
    "({ps_1}) @ ({ps_2}) = {ps_1@ps_2}\"\"\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "59c34f41",
   "metadata": {},
   "source": [
    "### Create an Observable with a Pauli string"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "5b0537a1",
   "metadata": {},
   "outputs": [],
   "source": [
    "obs1 = Observable(ps_1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3b411669",
   "metadata": {},
   "source": [
    "Since there is an equivalence between definition from a Pauli string or from a\n",
    "Hermitian matrix, both these observable wan can also be expressed in term of the\n",
    "mean through which it was not defined (Pauli for matrix and vice versa):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "1732e099",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "`obs` created with matrix:\n",
      "matrix:\n",
      "[[4, 2 , 3 , 8],\n",
      " [2, -3, 1 , 0],\n",
      " [3, 1 , -1, 5],\n",
      " [8, 0 , 5 , 2]]\n",
      "Pauli string:\n",
      "0.5*I@I + 3.5*I@X + 1*I@Z + 1.5*X@I + 4.5*X@X + 1.5*X@Z - 3.5*Y@Y - 1.5*Z@X + 2.5*Z@Z\n",
      "\n",
      "\n",
      "`obs1` created with Pauli string:\n",
      "Pauli string:\n",
      "1*I@Z - 3*X@Y\n",
      "matrix:\n",
      "[[1 , 0  , 0 , -3j],\n",
      " [0 , -1 , 3j, 0  ],\n",
      " [0 , -3j, 1 , 0  ],\n",
      " [3j, 0  , 0 , -1 ]]\n"
     ]
    }
   ],
   "source": [
    "from mpqp.tools import pprint\n",
    "\n",
    "\n",
    "print(\"`obs` created with matrix:\")\n",
    "print(\"matrix:\")\n",
    "pprint(obs.matrix)\n",
    "print(\"Pauli string:\")\n",
    "print(obs.pauli_string)\n",
    "\n",
    "\n",
    "print(\"\\n\\n`obs1` created with Pauli string:\")\n",
    "print(\"Pauli string:\")\n",
    "print(obs1.pauli_string)\n",
    "print(\"matrix:\")\n",
    "pprint(obs1.matrix)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a63b436e",
   "metadata": {},
   "source": [
    "### ExpectationMeasure\n",
    "\n",
    "Next, we measure the expectation value of the observable by adding it to the circuit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "892f0b6e",
   "metadata": {},
   "outputs": [],
   "source": [
    "circuit = circuit.without_measurements()\n",
    "circuit.add(ExpectationMeasure(obs1, shots=1000))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "745aa999",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "BatchResult: 3 results\n",
      "Result: circuit 1, ATOSDevice, MYQLM_PYLINALG\n",
      " Expectation value: -0.031999999999999945\n",
      " Error/Variance: 0.100047316133245\n",
      "Result: circuit 1, IBMDevice, AER_SIMULATOR\n",
      " Expectation value: 1.9081958235744878e-17\n",
      " Error/Variance: 0.0\n",
      "Result: circuit 1, ATOSDevice, MYQLM_CLINALG\n",
      " Expectation value: -0.29800000000000015\n",
      " Error/Variance: 0.0996658079447173\n"
     ]
    }
   ],
   "source": [
    "results = run(\n",
    "    circuit,\n",
    "    [\n",
    "        ATOSDevice.MYQLM_PYLINALG,\n",
    "        IBMDevice.AER_SIMULATOR,\n",
    "        ATOSDevice.MYQLM_CLINALG,\n",
    "    ],\n",
    ")\n",
    "print(results)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ad9fed49",
   "metadata": {},
   "outputs": [],
   "source": [
    "#TODO: show examples with multi observable expectation value\n",
    "\n",
    "#TODO: maybe show also the set of commuting pauli string when those observables are decomposed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a6f5f27c",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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": 5
}