{
  "cells": [
    {
      "cell_type": "markdown",
      "id": "frontmatter",
      "metadata": {},
      "source": [
        "---\n",
        "title: \"Projeção do operador de Pauli de referência\"\n",
        "description: \"Projeção do operador de Pauli de referência para a versão mais recente da diagonalização quântica baseada em amostras (SQD)\"\n",
        "---\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "55900f6c-2fb6-4d2c-8745-29bad8b66d9f",
      "metadata": {},
      "source": [
        "<span id=\"benchmark-pauli-operator-projection\" />\n",
        "\n",
        "# Projeção do operador de Pauli de referência\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "b5caa5b4",
      "metadata": {},
      "source": [
        "<span id=\"pauli-string-action-on-a-computational-basis-state\" />\n",
        "\n",
        "### Ação da corda de Pauli em um estado fundamental computacional\n",
        "\n",
        "A ação de uma corda de Pauli sobre um estado de base computacional é bastante trivial e consiste, ela própria, em um único estado de base computacional. Essa é uma consequência direta da estrutura das matrizes de Pauli, que possuem apenas um único elemento diferente de zero em cada uma de suas linhas. Consequentemente, sua ação sobre um qubit é:\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 da sequência de bits que identifica a base computacional será designado por $x \\in \\{0, 1 \\}$ Para manter a implementação o mais leve possível, vamos\n",
        "representar as sequências de bits com `bool` as variáveis: $0\\rightarrow \\textrm{False}$ e\n",
        "$1\\rightarrow \\textrm{True}$.\n",
        "\n",
        "Para representar a ação de cada operador de Pauli em um estado de base computacional,\n",
        "atribuiremos três variáveis a ele: `diag`, `sign`, `imag`.\n",
        "\n",
        "* `diag` indica se o operador é 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` Identifica se há uma mudança de sinal no elemento da matriz associado\n",
        "  a 0 ou a 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` Identifica se há um componente complexo no elemento da matriz:\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 um operador de Pauli arbitrário como $\\sigma \\in \\{ I, \\sigma_x, \\sigma_y \\sigma_z\\}$. A ação do operador de Pauli sobre um estado em base computacional\n",
        "pode, então, ser representada pela operação 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",
        "Isso se generaliza diretamente para um número arbitrário de qubits.\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "a8fd7e11",
      "metadata": {},
      "source": [
        "Vamos verificar se isso 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": [
        "Geramos um grande número de sequências de bits (50 milhões) para um 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",
        "### Funções de projeção de Pauli no modelo SQD de referência\n",
        "\n",
        "A corda de Pauli em questão é $\\sigma_z \\otimes ... \\otimes \\sigma_z$.\n",
        "\n",
        "Diferentes dimensões do subespaço são consideradas por meio do corte da matriz de sequências de bits. Medimos o tempo de projeção do subespaço para os diferentes tamanhos de subespaço.\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": [
        "Agora faremos o mesmo para 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
}