{
  "cells": [
    {
      "cell_type": "markdown",
      "id": "frontmatter",
      "metadata": {},
      "source": [
        "---\n",
        "title: \"Proyección del operador de Pauli de referencia\"\n",
        "description: \"Proyección del operador de Pauli de referencia para la última versión de la diagonalización cuántica basada en muestras (SQD)\"\n",
        "---\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "55900f6c-2fb6-4d2c-8745-29bad8b66d9f",
      "metadata": {},
      "source": [
        "<span id=\"benchmark-pauli-operator-projection\" />\n",
        "\n",
        "# Proyección del operador de Pauli de referencia\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "b5caa5b4",
      "metadata": {},
      "source": [
        "<span id=\"pauli-string-action-on-a-computational-basis-state\" />\n",
        "\n",
        "### Acción de la cuerda de Pauli en un estado fundamental computacional\n",
        "\n",
        "La acción de una cuerda de Pauli sobre un estado de base computacional es bastante trivial, y constituye en sí misma un único estado de base computacional. Esto es una consecuencia directa de la estructura de las matrices de Pauli, que solo tienen un único elemento distinto de cero en cada una de sus filas. Por lo tanto, su acción sobre un qubit es:\n",
        "\n",
        "***\n",
        "\n",
        "$$\n",
        "\\sigma_x |0 \\rangle = |1 \\rangle\n",
        "$$\n",
        "\n",
        "$$\n",
        "\\sigma_x |1 \\rangle = |0 \\rangle\n",
        "$$\n",
        "\n",
        "***\n",
        "\n",
        "$$\n",
        "\\sigma_y |0 \\rangle = i|1 \\rangle\n",
        "$$\n",
        "\n",
        "$$\n",
        "\\sigma_y |1 \\rangle = -i|0 \\rangle\n",
        "$$\n",
        "\n",
        "***\n",
        "\n",
        "$$\n",
        "\\sigma_z |0 \\rangle = |0 \\rangle\n",
        "$$\n",
        "\n",
        "$$\n",
        "\\sigma_z |1 \\rangle = -|1 \\rangle\n",
        "$$\n",
        "\n",
        "***\n",
        "\n",
        "$$\n",
        "I |0 \\rangle = |0 \\rangle\n",
        "$$\n",
        "\n",
        "$$\n",
        "I |1 \\rangle = |1 \\rangle\n",
        "$$\n",
        "\n",
        "***\n",
        "\n",
        "Cada bit de la cadena de bits que identifica la base computacional se etiquetará con un $x \\in \\{0, 1 \\}$ e. Para que la implementación sea lo más ligera posible,\n",
        "representaremos las cadenas de bits con `bool` las variables: $0\\rightarrow \\textrm{False}$ y\n",
        "$1\\rightarrow \\textrm{True}$.\n",
        "\n",
        "Para representar la acción de cada operador de Pauli en un estado de base computacional,\n",
        "le asignaremos tres variables: `diag`, `sign`, `imag`.\n",
        "\n",
        "* `diag` indica si el operador es diagonal:\n",
        "\n",
        "  * $\\textrm{diag}(I) = \\textrm{True}$\n",
        "  * $\\textrm{diag}(\\sigma_x) = \\textrm{False}$\n",
        "  * $\\textrm{diag}(\\sigma_y) = \\textrm{False}$\n",
        "  * $\\textrm{diag}(\\sigma_z) = \\textrm{True}$\n",
        "\n",
        "* `sign` Determina si se produce un cambio de signo en el elemento de la matriz asociado\n",
        "  al 0 o al 1:\n",
        "\n",
        "  * $\\textrm{sign}(I) = \\textrm{False}$\n",
        "  * $\\textrm{sign}(\\sigma_x) = \\textrm{False}$\n",
        "  * $\\textrm{sign}(\\sigma_y) = \\textrm{True}$\n",
        "  * $\\textrm{sign}(\\sigma_z) = \\textrm{True}$\n",
        "\n",
        "* `imag` Indica si el elemento de la matriz contiene un componente complejo:\n",
        "\n",
        "  * $\\textrm{imag}(I) = \\textrm{False}$\n",
        "  * $\\textrm{imag}(\\sigma_x) = \\textrm{False}$\n",
        "  * $\\textrm{imag}(\\sigma_y) = \\textrm{True}$\n",
        "  * $\\textrm{imag}(\\sigma_z) = \\textrm{False}$\n",
        "\n",
        "Denominamos « $\\sigma \\in \\{ I, \\sigma_x, \\sigma_y \\sigma_z\\}$ » a un operador de Pauli arbitrario. La acción del operador de Pauli sobre un estado de base computacional\n",
        "puede representarse entonces mediante la operación lógica:\n",
        "\n",
        "$$\n",
        "\\sigma |x \\rangle = |x == \\textrm{diag}(\\sigma) \\rangle (-1)^{x\\textrm{ and sign}(\\sigma)}\n",
        "(i)^{\\textrm{imag}(\\sigma)}.\n",
        "$$\n",
        "\n",
        "Esto se generaliza directamente a un número arbitrario de qubits.\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "a8fd7e11",
      "metadata": {},
      "source": [
        "Comprobemos si esto funciona:\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "id": "dcb15308",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "-------------------\n",
            "I\n",
            "|False> -->  |False>    ME:(1+0j)\n",
            "|True> -->  |True>    ME:(1+0j)\n",
            "-------------------\n",
            "SX\n",
            "|False> -->  |True>    ME:(1+0j)\n",
            "|True> -->  |False>    ME:(1+0j)\n",
            "-------------------\n",
            "SZ\n",
            "|False> -->  |False>    ME:(1+0j)\n",
            "|True> -->  |True>    ME:(-1+0j)\n",
            "-------------------\n",
            "SY\n",
            "|False> -->  |True>    ME:1j\n",
            "|True> -->  |False>    ME:(-0-1j)\n"
          ]
        }
      ],
      "source": [
        "def connected_element_and_amplitude_bool(\n",
        "    x: bool, diag: bool, sign: bool, imag: bool\n",
        ") -> tuple[bool, complex]:\n",
        "    \"\"\"\n",
        "    Finds the connected element to computational basis state |x> under\n",
        "    the action of the Pauli operator represented by (diag, sign, imag).\n",
        "\n",
        "    Args:\n",
        "        x: Value of the bit, either True or False.\n",
        "        diag: Whether the Pauli operator is diagonal (I, Z)\n",
        "        sigma: Whether the Pauli operator's rows differ in sign (Y, Z)\n",
        "        imag: Whether the Pauli operator is purely imaginary (Y)\n",
        "\n",
        "    Returns:\n",
        "        A length-2 tuple:\n",
        "            - The connected element to x, either False or True\n",
        "            - The matrix element\n",
        "    \"\"\"\n",
        "    return x == diag, (-1) ** (x and sign) * (1j) ** (imag)\n",
        "\n",
        "\n",
        "sigma_indices = [0, 1, 2, 3]\n",
        "sigma_string = [\"I\", \"SX\", \"SZ\", \"SY\"]\n",
        "sigma_diag = [True, False, True, False]\n",
        "sigma_sign = [False, False, True, True]\n",
        "sigma_imag = [False, False, False, True]\n",
        "qubit_values = [False, True]\n",
        "\n",
        "for xi in sigma_indices:\n",
        "    print(\"-------------------\")\n",
        "    print(sigma_string[xi])\n",
        "    for x in qubit_values:\n",
        "        x_p, matrix_element = connected_element_and_amplitude_bool(\n",
        "            x, sigma_diag[xi], sigma_sign[xi], sigma_imag[xi]\n",
        "        )\n",
        "        print(\n",
        "            \"|\"\n",
        "            + str(x)\n",
        "            + \"> -->  |\"\n",
        "            + str(x_p)\n",
        "            + \">    ME:\"\n",
        "            + str(matrix_element)\n",
        "        )"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "f40a3463",
      "metadata": {},
      "source": [
        "Generamos un gran número de cadenas de bits (50 M) para un sistema de 40 qubits:\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "74c16c91-cc5c-46ce-aff8-17d0e71ac50f",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Total number of unique bitstrings: 49998839\n"
          ]
        }
      ],
      "source": [
        "import numpy as np\n",
        "from qiskit_addon_sqd.qubit import sort_and_remove_duplicates\n",
        "\n",
        "rand_seed = 22\n",
        "np.random.seed(rand_seed)\n",
        "\n",
        "# Generate some random bitstrings for testing\n",
        "\n",
        "\n",
        "def random_bitstrings(n_samples, n_qubits):\n",
        "    return (\n",
        "        np.round(np.random.rand(n_samples, n_qubits))\n",
        "        .astype(\"int\")\n",
        "        .astype(\"bool\")\n",
        "    )\n",
        "\n",
        "\n",
        "n_qubits = 40\n",
        "bts_matrix = random_bitstrings(50_000_000, n_qubits)\n",
        "\n",
        "# We need to sort the bitstrings and only keep the unique ones\n",
        "# NOTE: It is essential for the projection code to have the bitstrings sorted!\n",
        "bts_matrix = sort_and_remove_duplicates(bts_matrix).astype(\"bool\")\n",
        "\n",
        "# Final subspace dimension after getting rid of duplicated bitstrings\n",
        "d = bts_matrix.shape[0]\n",
        "\n",
        "print(\"Total number of unique bitstrings: \" + str(d))"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "184bf287",
      "metadata": {},
      "source": [
        "<span id=\"benchmark-sqd-pauli-projection-functions\" />\n",
        "\n",
        "### Funciones de proyección de Pauli del modelo de referencia SQD\n",
        "\n",
        "La cuerda de Pauli que nos ocupa es $\\sigma_z \\otimes ... \\otimes \\sigma_z$.\n",
        "\n",
        "Se tienen en cuenta diferentes dimensiones de subespacios mediante la división de la matriz de cadenas de bits. Mide el tiempo que tarda la proyección en el subespacio para los distintos tamaños de subespacio.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 3,
      "id": "8fe182bc",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Iteration 0 took 0.201246s\n",
            "Iteration 1 took 0.348222s\n",
            "Iteration 2 took 0.576333s\n",
            "Iteration 3 took 0.78356s\n",
            "Iteration 4 took 1.016162s\n",
            "Iteration 5 took 1.305325s\n",
            "Iteration 6 took 1.392751s\n",
            "Iteration 7 took 1.632433s\n",
            "Iteration 8 took 1.826521s\n",
            "Iteration 9 took 2.02903s\n",
            "Iteration 10 took 2.297458s\n",
            "Iteration 11 took 2.588042s\n",
            "Iteration 12 took 2.738746s\n",
            "Iteration 13 took 2.906144s\n",
            "Iteration 14 took 3.148833s\n",
            "Iteration 15 took 3.323253s\n",
            "Iteration 16 took 3.664171s\n",
            "Iteration 17 took 3.680663s\n",
            "Iteration 18 took 4.008313s\n",
            "Iteration 19 took 4.173532s\n"
          ]
        }
      ],
      "source": [
        "import time\n",
        "\n",
        "from qiskit.quantum_info import Pauli\n",
        "from qiskit_addon_sqd.qubit import matrix_elements_from_pauli\n",
        "\n",
        "pauli = Pauli(\"Z\" * n_qubits)\n",
        "\n",
        "# Different subspace sizes to test\n",
        "d_list = np.linspace(d / 1000, d, 20).astype(\"int\")\n",
        "\n",
        "# To store the walltime\n",
        "time_array = np.zeros(20)\n",
        "\n",
        "for i in range(20):\n",
        "    int_bts_matrix = bts_matrix[: d_list[i], :]\n",
        "    time_1 = time.time()\n",
        "    _ = matrix_elements_from_pauli(int_bts_matrix, pauli)\n",
        "    time_array[i] = time.time() - time_1\n",
        "    print(f\"Iteration {i} took {round(time_array[i], 6)}s\")"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 4,
      "id": "9abb110f",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<Image src=\"/docs/images/addons/qiskit-addon-sqd/guides/benchmark-pauli-projection/extracted-outputs/9abb110f-0.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "import matplotlib.pyplot as plt\n",
        "\n",
        "# Data for energies plot\n",
        "x1 = d_list\n",
        "y1 = time_array\n",
        "\n",
        "# Plot energies\n",
        "plt.title(\"Runtime vs subspace dimension 40 qubits\")\n",
        "plt.xlabel(\"Subspace dimension (millions)\")\n",
        "plt.ylabel(\"Wall time [s]\")\n",
        "plt.xticks([1e7, 2e7, 3e7, 4e7, 5e7], [str(i) for i in [10, 20, 30, 40, 50]])\n",
        "plt.plot(x1, y1, marker=\".\", markersize=20)\n",
        "plt.tight_layout()\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "8c486bc7",
      "metadata": {},
      "source": [
        "Ahora hacemos lo mismo con 60 qubits:\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 5,
      "id": "359ed3f3",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Total number of unique bitstrings: 50000000\n"
          ]
        }
      ],
      "source": [
        "n_qubits = 60\n",
        "bts_matrix = random_bitstrings(50_000_000, n_qubits)\n",
        "\n",
        "# We need to sort the bitstrings and just keep the unique ones\n",
        "bts_matrix = sort_and_remove_duplicates(bts_matrix).astype(\"bool\")\n",
        "\n",
        "# Final subspace dimension after getting rid of duplicated bitstrings\n",
        "d = bts_matrix.shape[0]\n",
        "\n",
        "print(\"Total number of unique bitstrings: \" + str(d))"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 6,
      "id": "7bb0d8d8",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Iteration 0 took 0.236567s\n",
            "Iteration 1 took 0.424116s\n",
            "Iteration 2 took 0.673399s\n",
            "Iteration 3 took 0.905164s\n",
            "Iteration 4 took 1.168936s\n",
            "Iteration 5 took 1.454204s\n",
            "Iteration 6 took 1.74778s\n",
            "Iteration 7 took 1.920795s\n",
            "Iteration 8 took 2.259994s\n",
            "Iteration 9 took 2.550674s\n",
            "Iteration 10 took 2.681287s\n",
            "Iteration 11 took 3.04411s\n",
            "Iteration 12 took 3.293262s\n",
            "Iteration 13 took 3.471247s\n",
            "Iteration 14 took 3.726639s\n",
            "Iteration 15 took 4.072854s\n",
            "Iteration 16 took 4.221037s\n",
            "Iteration 17 took 4.498535s\n",
            "Iteration 18 took 4.741108s\n",
            "Iteration 19 took 5.159038s\n"
          ]
        }
      ],
      "source": [
        "pauli = Pauli(\"Z\" * n_qubits)\n",
        "\n",
        "# Different subspace sizes to test\n",
        "d_list = np.linspace(d / 1000, d, 20).astype(\"int\")\n",
        "\n",
        "# It is better to do this once\n",
        "row_array = np.arange(d)\n",
        "\n",
        "# To store the walltime\n",
        "time_array = np.zeros(20)\n",
        "\n",
        "for i in range(20):\n",
        "    int_bts_matrix = bts_matrix[: d_list[i], :]\n",
        "    int_row_array = row_array[: d_list[i]]\n",
        "    time_1 = time.time()\n",
        "    _ = matrix_elements_from_pauli(int_bts_matrix, pauli)\n",
        "    time_array[i] = time.time() - time_1\n",
        "    print(f\"Iteration {i} took {round(time_array[i], 6)}s\")"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 7,
      "id": "6b961c81",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<Image src=\"/docs/images/addons/qiskit-addon-sqd/guides/benchmark-pauli-projection/extracted-outputs/6b961c81-0.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "# Data for energies plot\n",
        "x1 = d_list\n",
        "y1 = time_array\n",
        "\n",
        "fig, axs = plt.subplots(1, 1, figsize=(6, 6))\n",
        "\n",
        "# Plot energies\n",
        "axs.plot(x1, y1, marker=\".\", markersize=20)\n",
        "axs.set_title(\"Runtime vs subspace dimension 60 qubits\")\n",
        "axs.set_xlabel(\"Subspace dimension (millions)\")\n",
        "plt.xticks([1e7, 2e7, 3e7, 4e7, 5e7], [str(i) for i in [10, 20, 30, 40, 50]])\n",
        "axs.set_ylabel(\"Wall time [s]\")\n",
        "\n",
        "plt.tight_layout()\n",
        "plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "id": "a1b8767d",
      "source": "© IBM Corp., 2017-2026"
    }
  ],
  "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"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 5
}