N_oscillators_case_scalability.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/NSPFD/quantum-computing/blob/main/N_oscillators_case_scalability.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "97fb99d0-2ebc-4143-a4ba-ea1e3c874b4a",
      "metadata": {
        "id": "97fb99d0-2ebc-4143-a4ba-ea1e3c874b4a"
      },
      "source": [
        "## N-Harmonic oscillators\n",
        "\n",
        "### Model description\n",
        "This notebook implements the same algorithmic implementation as the one of the $N=2$ harmonic oscillators extending our computations from $N=2$ to larger numbers of oscillators"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "0c3f0d87-55a5-4db3-8e32-e0a8df957aeb",
      "metadata": {
        "id": "0c3f0d87-55a5-4db3-8e32-e0a8df957aeb"
      },
      "source": [
        "$\\underline{\\textbf{Important comment}}$: Here, we fix all the masses to be equal to 1000 and all spring constants to be equal to 1 and we vary the number of oscillators N. We study two different connectivity topologies: 1) linear chain, i.e. each oscillator is connected with only two others (one on the left and on on the right) except the ones at the boundaries of the chain, 2) full connectivity where each oscillator is coupled to all the other $N-1$ oscillators."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "415dc16a-cfe4-4b92-bd99-a7c6f102a94c",
      "metadata": {
        "id": "415dc16a-cfe4-4b92-bd99-a7c6f102a94c",
        "outputId": "46f44114-3eeb-41b7-f1bd-87b30d7d3aec"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "mass matrix: (M):\n",
            "[[10  0  0  0]\n",
            " [ 0 10  0  0]\n",
            " [ 0  0 10  0]\n",
            " [ 0  0  0 10]]\n",
            "\n",
            "(Spring Constants Matrix) K:\n",
            "[[1. 1. 0. 0.]\n",
            " [1. 1. 1. 0.]\n",
            " [0. 1. 1. 1.]\n",
            " [0. 0. 1. 1.]]\n",
            "\n",
            " Sq_K matrix:\n",
            "[[1. 1. 0. 0.]\n",
            " [1. 1. 1. 0.]\n",
            " [0. 1. 1. 1.]\n",
            " [0. 0. 1. 1.]]\n",
            "\n",
            " F matrix:\n",
            "[[ 2. -1. -0. -0.]\n",
            " [-1.  3. -1. -0.]\n",
            " [-0. -1.  3. -1.]\n",
            " [-0. -0. -1.  2.]]\n",
            "\n",
            " A matrix:\n",
            "[[ 0.2 -0.1  0.   0. ]\n",
            " [-0.1  0.3 -0.1  0. ]\n",
            " [ 0.  -0.1  0.3 -0.1]\n",
            " [ 0.   0.  -0.1  0.2]]\n",
            "\n",
            "Q Matrix:\n",
            "[[ 1.  1.  0.  0.  0.  0.  0.  0.  0.  0.]\n",
            " [ 0. -1.  0.  0.  1.  1.  0.  0.  0.  0.]\n",
            " [ 0.  0. -0.  0.  0. -1.  0.  1.  1.  0.]\n",
            " [ 0.  0.  0. -0.  0.  0. -0.  0. -1.  1.]]\n",
            "\n",
            "B Matrix:\n",
            "[[ 0.31622777  0.31622777  0.          0.          0.          0.\n",
            "   0.          0.          0.          0.        ]\n",
            " [ 0.         -0.31622777  0.          0.          0.31622777  0.31622777\n",
            "   0.          0.          0.          0.        ]\n",
            " [ 0.          0.          0.          0.          0.         -0.31622777\n",
            "   0.          0.31622777  0.31622777  0.        ]\n",
            " [ 0.          0.          0.          0.          0.          0.\n",
            "   0.          0.         -0.31622777  0.31622777]]\n",
            "\n",
            "Padded B Matrix (B_pad):\n",
            "[[ 0.31622777  0.31622777  0.          0.          0.          0.\n",
            "   0.          0.          0.          0.          0.          0.\n",
            "   0.          0.          0.          0.        ]\n",
            " [ 0.         -0.31622777  0.          0.          0.31622777  0.31622777\n",
            "   0.          0.          0.          0.          0.          0.\n",
            "   0.          0.          0.          0.        ]\n",
            " [ 0.          0.          0.          0.          0.         -0.31622777\n",
            "   0.          0.31622777  0.31622777  0.          0.          0.\n",
            "   0.          0.          0.          0.        ]\n",
            " [ 0.          0.          0.          0.          0.          0.\n",
            "   0.          0.         -0.31622777  0.31622777  0.          0.\n",
            "   0.          0.          0.          0.        ]\n",
            " [ 0.          0.          0.          0.          0.          0.\n",
            "   0.          0.          0.          0.          0.          0.\n",
            "   0.          0.          0.          0.        ]\n",
            " [ 0.          0.          0.          0.          0.          0.\n",
            "   0.          0.          0.          0.          0.          0.\n",
            "   0.          0.          0.          0.        ]\n",
            " [ 0.          0.          0.          0.          0.          0.\n",
            "   0.          0.          0.          0.          0.          0.\n",
            "   0.          0.          0.          0.        ]\n",
            " [ 0.          0.          0.          0.          0.          0.\n",
            "   0.          0.          0.          0.          0.          0.\n",
            "   0.          0.          0.          0.        ]\n",
            " [ 0.          0.          0.          0.          0.          0.\n",
            "   0.          0.          0.          0.          0.          0.\n",
            "   0.          0.          0.          0.        ]\n",
            " [ 0.          0.          0.          0.          0.          0.\n",
            "   0.          0.          0.          0.          0.          0.\n",
            "   0.          0.          0.          0.        ]\n",
            " [ 0.          0.          0.          0.          0.          0.\n",
            "   0.          0.          0.          0.          0.          0.\n",
            "   0.          0.          0.          0.        ]\n",
            " [ 0.          0.          0.          0.          0.          0.\n",
            "   0.          0.          0.          0.          0.          0.\n",
            "   0.          0.          0.          0.        ]\n",
            " [ 0.          0.          0.          0.          0.          0.\n",
            "   0.          0.          0.          0.          0.          0.\n",
            "   0.          0.          0.          0.        ]\n",
            " [ 0.          0.          0.          0.          0.          0.\n",
            "   0.          0.          0.          0.          0.          0.\n",
            "   0.          0.          0.          0.        ]\n",
            " [ 0.          0.          0.          0.          0.          0.\n",
            "   0.          0.          0.          0.          0.          0.\n",
            "   0.          0.          0.          0.        ]\n",
            " [ 0.          0.          0.          0.          0.          0.\n",
            "   0.          0.          0.          0.          0.          0.\n",
            "   0.          0.          0.          0.        ]]\n",
            "\n",
            "Hamiltonian (H):\n",
            "[[ 0.          0.          0.          0.          0.31622777  0.31622777\n",
            "   0.          0.          0.          0.          0.          0.\n",
            "   0.          0.        ]\n",
            " [ 0.          0.          0.          0.          0.         -0.31622777\n",
            "   0.          0.          0.31622777  0.31622777  0.          0.\n",
            "   0.          0.        ]\n",
            " [ 0.          0.          0.          0.          0.          0.\n",
            "   0.          0.          0.         -0.31622777  0.          0.31622777\n",
            "   0.31622777  0.        ]\n",
            " [ 0.          0.          0.          0.          0.          0.\n",
            "   0.          0.          0.          0.          0.          0.\n",
            "  -0.31622777  0.31622777]\n",
            " [ 0.31622777  0.          0.          0.          0.          0.\n",
            "   0.          0.          0.          0.          0.          0.\n",
            "   0.          0.        ]\n",
            " [ 0.31622777 -0.31622777  0.          0.          0.          0.\n",
            "   0.          0.          0.          0.          0.          0.\n",
            "   0.          0.        ]\n",
            " [ 0.          0.          0.          0.          0.          0.\n",
            "   0.          0.          0.          0.          0.          0.\n",
            "   0.          0.        ]\n",
            " [ 0.          0.          0.          0.          0.          0.\n",
            "   0.          0.          0.          0.          0.          0.\n",
            "   0.          0.        ]\n",
            " [ 0.          0.31622777  0.          0.          0.          0.\n",
            "   0.          0.          0.          0.          0.          0.\n",
            "   0.          0.        ]\n",
            " [ 0.          0.31622777 -0.31622777  0.          0.          0.\n",
            "   0.          0.          0.          0.          0.          0.\n",
            "   0.          0.        ]\n",
            " [ 0.          0.          0.          0.          0.          0.\n",
            "   0.          0.          0.          0.          0.          0.\n",
            "   0.          0.        ]\n",
            " [ 0.          0.          0.31622777  0.          0.          0.\n",
            "   0.          0.          0.          0.          0.          0.\n",
            "   0.          0.        ]\n",
            " [ 0.          0.          0.31622777 -0.31622777  0.          0.\n",
            "   0.          0.          0.          0.          0.          0.\n",
            "   0.          0.        ]\n",
            " [ 0.          0.          0.          0.31622777  0.          0.\n",
            "   0.          0.          0.          0.          0.          0.\n",
            "   0.          0.        ]]\n",
            "\n",
            "Padded Hamiltonian (H_pad):\n",
            "[[0. 0. 0. ... 0. 0. 0.]\n",
            " [0. 0. 0. ... 0. 0. 0.]\n",
            " [0. 0. 0. ... 0. 0. 0.]\n",
            " ...\n",
            " [0. 0. 0. ... 0. 0. 0.]\n",
            " [0. 0. 0. ... 0. 0. 0.]\n",
            " [0. 0. 0. ... 0. 0. 0.]]\n"
          ]
        }
      ],
      "source": [
        "import numpy as np\n",
        "\n",
        "# The number N of oscillators needs to be a power of 2. Specify first the number of qubits n the system will correspond to.\n",
        "n=2\n",
        "N=2**n\n",
        "\n",
        "# Create the mass matrix M as an N x N diagonal matrix with all diagonal elements equal to 1000\n",
        "M = np.diag([10] * N)\n",
        "\n",
        "# Create the spring constants matrix K as an N x N symmetric matrix with all elements equal to 1\n",
        "#K = np.ones((N, N)) #This command is for full oscillators connectivity\n",
        "\n",
        "#The following command is for oscillators in a linear chain arrangement\n",
        "# Initialize the K matrix with zeros\n",
        "K = np.zeros((N, N))\n",
        "\n",
        "# Fill the K matrix according to the specified pattern\n",
        "for i in range(N):\n",
        "    if i > 0:\n",
        "        K[i, i-1] = 1  # Element (i, i-1)\n",
        "    K[i, i] = 1  # Element (i, i)\n",
        "    if i < N-1:\n",
        "        K[i, i+1] = 1  # Element (i, i+1)\n",
        "\n",
        "# Compute the inverse square root of the mass matrix\n",
        "inv_sqrt_M = np.diag(1 / np.sqrt(np.diag(M)))\n",
        "\n",
        "# Sq_K matrix whose elements are the square root of the elements of the K matrix\n",
        "Sq_K = np.sqrt(K)\n",
        "\n",
        "# F-matrix: Constructed using the spring constants\n",
        "F = np.zeros((N, N))\n",
        "for j in range(N):\n",
        "    F[j, j] = np.sum(K[j, :])  # Diagonal elements\n",
        "    for k in range(N):\n",
        "        if j != k:\n",
        "            F[j, k] = -K[j, k]  # Off-diagonal elements\n",
        "\n",
        "#A-matrix\n",
        "A = inv_sqrt_M @ F @ inv_sqrt_M\n",
        "\n",
        "# Initialize the Q matrix (N x (N*(N+1)/2))\n",
        "Q = np.zeros((N, N * (N + 1) // 2))\n",
        "\n",
        "# Fill the Q matrix based on the formula\n",
        "column = 0\n",
        "for j in range(N):\n",
        "    for k in range(j, N):\n",
        "        if j == k:\n",
        "            Q[j, column] = Sq_K[j, k]\n",
        "        else:\n",
        "            Q[j, column] = Sq_K[j, k]\n",
        "            Q[k, column] = -Sq_K[j, k]\n",
        "        column += 1\n",
        "\n",
        "# Compute the B matrix\n",
        "B = inv_sqrt_M @ Q\n",
        "\n",
        "# Create the padded B matrix B_pad\n",
        "B_pad = np.zeros((N**2, N**2))\n",
        "B_pad[:N, :B.shape[1]] = B\n",
        "\n",
        "# Construct the Hamiltonian H\n",
        "size_H = N + N * (N + 1) // 2\n",
        "H = np.zeros((size_H, size_H))\n",
        "H[:N, N:] = B\n",
        "H[N:, :N] = B.T\n",
        "\n",
        "# Construct the padded Hamiltonian H_pad\n",
        "size_H_pad = 2 * N**2\n",
        "H_pad = np.zeros((size_H_pad, size_H_pad))\n",
        "H_pad[:N**2, N**2:] = B_pad\n",
        "H_pad[N**2:, :N**2] = B_pad.T\n",
        "\n",
        "# Print the matrices\n",
        "print(\"mass matrix: (M):\")\n",
        "print(M)\n",
        "print(\"\\n(Spring Constants Matrix) K:\")\n",
        "print(K)\n",
        "print(\"\\n Sq_K matrix:\")\n",
        "print(Sq_K)\n",
        "print(\"\\n F matrix:\")\n",
        "print(F)\n",
        "print(\"\\n A matrix:\")\n",
        "print(A)\n",
        "print(\"\\nQ Matrix:\")\n",
        "print(Q)\n",
        "print(\"\\nB Matrix:\")\n",
        "print(B)\n",
        "print(\"\\nPadded B Matrix (B_pad):\")\n",
        "print(B_pad)\n",
        "print(\"\\nHamiltonian (H):\")\n",
        "print(H)\n",
        "print(\"\\nPadded Hamiltonian (H_pad):\")\n",
        "print(H_pad)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "9e5f8eed-e560-4e46-bb47-4bf15dd1e415",
      "metadata": {
        "id": "9e5f8eed-e560-4e46-bb47-4bf15dd1e415"
      },
      "source": [
        "### LCU Pauli Strings decomposition of the Hamiltonian matrix $\\bar{\\textbf{H}}$\n",
        "\n",
        "In the notebook for the two oscillators, we discussed how to decompose a Hamiltonian matrix to the basis of the Pauli matrices together with the identity. For the case of $N=2$ oscillators, the padded Hamiltonian, performing the evolution, is $8\\times 8$. For the general case of $N$ oscillators, the padded Hamiltonian is $2N^2\\times 2N^2$, where $N=2^n$, i.e. is a power of 2. In other words, the padded Hamiltonian is a $(2^{2n+1})\\times (2^{2n+1})$\n",
        "\n",
        "$\\textbf{Important clarification point}$: The padded Hamiltonian $\\bar{\\textbf{H}}$ has dimensionality $2N^2\\times 2N^2$, where $N=2^n$, i.e. is a power of 2. In other words, the padded Hamiltonian is a $(2^{2n+1})\\times (2^{2n+1})$. This is the dimensionality of a Hamiltonian acting on the Hilbert space of 3 qubits. For the case of the 2 oscillators, the Hamiltonian is equivalent to one acting on the space of 3 qubits. Hence, when the Hamiltonian is decomposed to a series of Pauli strings, we need strings made of 3 Pauli matrices.\n",
        "\n",
        "Any hermitian operator acting on a Hilbert space $\\mathcal{H}$ can be decomposed to a sum of the identity and Pauli matrices tensor products as\n",
        "\n",
        "$\\textbf{H}=\\sum_{i_1=0}^3\\sum_{i_2=0}^3,...,\\sum_{i_{2n+1}=0}^3c_{i_1 i_2...i_{2n+1}}\\,\\textbf{P}_{i_1}\\otimes \\textbf{P}_{i_2} \\otimes \\cdot \\cdot \\cdot \\otimes \\textbf{P}_{i_{2n+1}},$\n",
        "\n",
        "where $\\textbf{P}_{i_1}\\otimes \\textbf{P}_{i_2} \\otimes \\cdot \\cdot \\cdot \\otimes \\textbf{P}_{i_{2n+1}}$ with $\\textbf{P}_i=\\left(\\textbf{I}, \\textbf{X}, \\textbf{Y}, \\textbf{Z}\\right)$ forms a basis with an orthonormal product given by\n",
        "\n",
        "\n",
        "For the $N$-oscillators, resulting to the $2n+1$-qubits case, as we explained in the 2_oscilators notebook, the orthonormality reads\n",
        "\n",
        "$\n",
        "\\langle \\textbf{P}_{i_1}\\otimes \\textbf{P}_{i_2} \\otimes \\cdot \\cdot \\cdot \\otimes \\textbf{P}_{i_{2n+1}}, \\textbf{P}_{i_1^\\prime}\\otimes \\textbf{P}_{i_2^\\prime} \\otimes \\cdot \\cdot \\cdot \\otimes \\textbf{P}_{i_{2n+1}^\\prime} \\rangle=\\frac{1}{2^{2n+1}}\\text{Tr}\\left[\\left(\\textbf{P}_{i_1}\\otimes \\textbf{P}_{i_2} \\otimes \\cdot \\cdot \\cdot \\otimes \\textbf{P}_{i_{2n+1}}\\right)\\left(\\textbf{P}_{i_1^\\prime}\\otimes \\textbf{P}_{i_2^\\prime} \\otimes \\cdot \\cdot \\cdot \\otimes \\textbf{P}_{i_{2n+1}^\\prime}\\right)\\right]=\\delta_{i_1 i_1^\\prime}\\delta_{i_2 i_2^\\prime}\\cdot \\cdot \\cdot\\delta_{i_{2n+1} i_{2n+1}^\\prime}.\n",
        "$\n",
        "\n",
        "The (padded) $2N^2\\times N^2$ Hamiltonian of the $N$ oscillators that maps the evolution to that of a $2n+1$-qubit system is decomposed to the LCU Paili basis $(\\textbf{P}_{i_1}\\otimes \\textbf{P}_{i_2} \\otimes \\cdot \\cdot \\cdot \\otimes \\textbf{P}_{i_{2n+1}})$ with its expansion coefficients given by\n",
        "\n",
        "$\n",
        "c_{i_1 i_2...i_{2n+1}}=\\frac{1}{2^{2n+1}}\\text{Tr}\\left[\\left(\\textbf{P}_{i_1}\\otimes \\textbf{P}_{i_2} \\otimes \\cdot \\cdot \\cdot \\otimes \\textbf{P}_{i_{2n+1}}\\right)\\textbf{H}\\right].\n",
        "$\n",
        "\n",
        "The code below computes the $c_{i_1 i_2...i_{2n+1}}$ coefficients from which it constructs the LCU decomposition of the Hamiltonian."
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "2b6e3f61-0c0c-41c5-87d9-557000c42dbc",
      "metadata": {
        "id": "2b6e3f61-0c0c-41c5-87d9-557000c42dbc",
        "outputId": "f9313753-be95-4262-9a08-5255826f95d2"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Decomposition of H_pad into Pauli strings:\n",
            "0.02 * XIIIX\n",
            "0.04 * XIIIZ\n",
            "0.02 * XIIZX\n",
            "0.04 * XIIZZ\n",
            "0.02 * XIXII\n",
            "0.04 * XIXIX\n",
            "-0.02 * XIXIZ\n",
            "-0.02 * XIXXX\n",
            "-0.02 * XIXYY\n",
            "0.02 * XIXZI\n",
            "-0.02 * XIXZZ\n",
            "0.02 * XIYXY\n",
            "-0.02 * XIYYX\n",
            "0.04 * XIYZY\n",
            "0.02 * XIZIX\n",
            "0.04 * XIZIZ\n",
            "0.02 * XIZZX\n",
            "0.04 * XIZZZ\n",
            "0.04 * XXIXI\n",
            "-0.02 * XXIXX\n",
            "0.02 * XXIYY\n",
            "0.04 * XXZXI\n",
            "-0.02 * XXZXX\n",
            "0.02 * XXZYY\n",
            "-0.02 * XYIXY\n",
            "0.04 * XYIYI\n",
            "-0.02 * XYIYX\n",
            "-0.02 * XYZXY\n",
            "0.04 * XYZYI\n",
            "-0.02 * XYZYX\n",
            "0.02 * XZIIX\n",
            "0.04 * XZIIZ\n",
            "0.02 * XZIZX\n",
            "0.04 * XZIZZ\n",
            "0.02 * XZXII\n",
            "0.04 * XZXIX\n",
            "-0.02 * XZXIZ\n",
            "-0.02 * XZXXX\n",
            "-0.02 * XZXYY\n",
            "0.02 * XZXZI\n",
            "-0.02 * XZXZZ\n",
            "0.02 * XZYXY\n",
            "-0.02 * XZYYX\n",
            "0.04 * XZYZY\n",
            "0.02 * XZZIX\n",
            "0.04 * XZZIZ\n",
            "0.02 * XZZZX\n",
            "0.04 * XZZZZ\n",
            "-0.02 * YIIIY\n",
            "-0.02 * YIIZY\n",
            "0.02 * YIXXY\n",
            "-0.02 * YIXYX\n",
            "0.04 * YIXZY\n",
            "-0.02 * YIYII\n",
            "-0.04 * YIYIX\n",
            "0.02 * YIYIZ\n",
            "0.02 * YIYXX\n",
            "0.02 * YIYYY\n",
            "-0.02 * YIYZI\n",
            "0.02 * YIYZZ\n",
            "-0.02 * YIZIY\n",
            "-0.02 * YIZZY\n",
            "-0.02 * YXIXY\n",
            "0.04 * YXIYI\n",
            "-0.02 * YXIYX\n",
            "-0.02 * YXZXY\n",
            "0.04 * YXZYI\n",
            "-0.02 * YXZYX\n",
            "-0.04 * YYIXI\n",
            "0.02 * YYIXX\n",
            "-0.02 * YYIYY\n",
            "-0.04 * YYZXI\n",
            "0.02 * YYZXX\n",
            "-0.02 * YYZYY\n",
            "-0.02 * YZIIY\n",
            "-0.02 * YZIZY\n",
            "0.02 * YZXXY\n",
            "-0.02 * YZXYX\n",
            "0.04 * YZXZY\n",
            "-0.02 * YZYII\n",
            "-0.04 * YZYIX\n",
            "0.02 * YZYIZ\n",
            "0.02 * YZYXX\n",
            "0.02 * YZYYY\n",
            "-0.02 * YZYZI\n",
            "0.02 * YZYZZ\n",
            "-0.02 * YZZIY\n",
            "-0.02 * YZZZY\n",
            "------------------------------------------------------------------------------------------------------------------------------------------\n",
            "Verification\n",
            "\n",
            "Original Hamiltonian H_pad:\n",
            "[[0. 0. 0. ... 0. 0. 0.]\n",
            " [0. 0. 0. ... 0. 0. 0.]\n",
            " [0. 0. 0. ... 0. 0. 0.]\n",
            " ...\n",
            " [0. 0. 0. ... 0. 0. 0.]\n",
            " [0. 0. 0. ... 0. 0. 0.]\n",
            " [0. 0. 0. ... 0. 0. 0.]]\n",
            "\n",
            "Reconstructed Hamiltonian H_reconstructed:\n",
            "[[0.+0.j 0.+0.j 0.+0.j ... 0.+0.j 0.+0.j 0.+0.j]\n",
            " [0.+0.j 0.+0.j 0.+0.j ... 0.+0.j 0.+0.j 0.+0.j]\n",
            " [0.+0.j 0.+0.j 0.+0.j ... 0.+0.j 0.+0.j 0.+0.j]\n",
            " ...\n",
            " [0.+0.j 0.+0.j 0.+0.j ... 0.+0.j 0.+0.j 0.+0.j]\n",
            " [0.+0.j 0.+0.j 0.+0.j ... 0.+0.j 0.+0.j 0.+0.j]\n",
            " [0.+0.j 0.+0.j 0.+0.j ... 0.+0.j 0.+0.j 0.+0.j]]\n",
            "\n",
            "Difference (H_pad - H_reconstructed):\n",
            "[[0.+0.j 0.+0.j 0.+0.j ... 0.+0.j 0.+0.j 0.+0.j]\n",
            " [0.+0.j 0.+0.j 0.+0.j ... 0.+0.j 0.+0.j 0.+0.j]\n",
            " [0.+0.j 0.+0.j 0.+0.j ... 0.+0.j 0.+0.j 0.+0.j]\n",
            " ...\n",
            " [0.+0.j 0.+0.j 0.+0.j ... 0.+0.j 0.+0.j 0.+0.j]\n",
            " [0.+0.j 0.+0.j 0.+0.j ... 0.+0.j 0.+0.j 0.+0.j]\n",
            " [0.+0.j 0.+0.j 0.+0.j ... 0.+0.j 0.+0.j 0.+0.j]]\n",
            "\n",
            "Reconstruction successful?: True\n"
          ]
        }
      ],
      "source": [
        "import numpy as np\n",
        "from itertools import product\n",
        "\n",
        "# Number of qubits\n",
        "k=2*(int(np.log2(N)))+1\n",
        "\n",
        "# Define Pauli matrices and identity matrix\n",
        "I = np.array([[1, 0], [0, 1]])\n",
        "X = np.array([[0, 1], [1, 0]])\n",
        "Y = np.array([[0, -1j], [1j, 0]])\n",
        "Z = np.array([[1, 0], [0, -1]])\n",
        "\n",
        "pauli_matrices = [I, X, Y, Z]\n",
        "pauli_labels = ['I', 'X', 'Y', 'Z']\n",
        "\n",
        "def generate_pauli_strings(k):\n",
        "    \"\"\"Generate all tensor products of k Pauli matrices\"\"\"\n",
        "    pauli_strings = []\n",
        "    labels = []\n",
        "    for prods in product(range(4), repeat=k):\n",
        "        ps = pauli_matrices[prods[0]]\n",
        "        label = pauli_labels[prods[0]]\n",
        "        for i in range(1, k):\n",
        "            ps = np.kron(ps, pauli_matrices[prods[i]])\n",
        "            label += pauli_labels[prods[i]]\n",
        "        pauli_strings.append(ps)\n",
        "        labels.append(label)\n",
        "    return pauli_strings, labels\n",
        "\n",
        "def decompose_hamiltonian(H):\n",
        "    \"\"\"Decompose Hamiltonian into Pauli strings\"\"\"\n",
        "    w = H.shape[0]\n",
        "    pauli_strings, labels = generate_pauli_strings(k)\n",
        "\n",
        "    #coeffs = [np.trace(np.dot(ps, H)) / w for ps in pauli_strings]\n",
        "    coeffs = [np.real(np.trace(np.dot(ps, H)) / w) for ps in pauli_strings]\n",
        "\n",
        "    decomposition = {label: coef for label, coef in zip(labels, coeffs)}\n",
        "\n",
        "    return decomposition, pauli_strings, coeffs\n",
        "\n",
        "def reconstruct_hamiltonian(pauli_strings, coefficients):\n",
        "    \"\"\"Reconstruct Hamiltonian from Pauli strings and coefficients\"\"\"\n",
        "    H_reconstructed = sum(c * ps for c, ps in zip(coefficients, pauli_strings))\n",
        "    return H_reconstructed\n",
        "\n",
        "# Decompose the Hamiltonian\n",
        "decomposition, pauli_strings, coeffs = decompose_hamiltonian(H_pad)\n",
        "\n",
        "# Print the decomposition\n",
        "print(\"Decomposition of H_pad into Pauli strings:\")\n",
        "for label, coef in decomposition.items():\n",
        "    if not np.isclose(coef, 0):\n",
        "        print(f\"{coef:.2f} * {label}\")\n",
        "print(\"------------------------------------------------------------------------------------------------------------------------------------------\")\n",
        "\n",
        "# Reconstruct the Hamiltonian\n",
        "H_reconstructed = reconstruct_hamiltonian(pauli_strings, coeffs)\n",
        "\n",
        "# Verification part\n",
        "print(\"Verification\")\n",
        "print(\"\\nOriginal Hamiltonian H_pad:\")\n",
        "print(H_pad)\n",
        "print(\"\\nReconstructed Hamiltonian H_reconstructed:\")\n",
        "print(H_reconstructed)\n",
        "print(\"\\nDifference (H_pad - H_reconstructed):\")\n",
        "print(H_pad - H_reconstructed)\n",
        "print(\"\\nReconstruction successful?:\", np.allclose(H_pad, H_reconstructed))\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "719eb699-87ce-47d6-911f-d1fc21be84f5",
      "metadata": {
        "id": "719eb699-87ce-47d6-911f-d1fc21be84f5",
        "outputId": "20dd026a-8bdc-47eb-91da-e19bacf8b29f"
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "[0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.01976423537605237,\n",
              " 0.0,\n",
              " 0.03952847075210474,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.01976423537605237,\n",
              " 0.0,\n",
              " 0.03952847075210474,\n",
              " 0.01976423537605237,\n",
              " 0.03952847075210474,\n",
              " 0.0,\n",
              " -0.01976423537605237,\n",
              " 0.0,\n",
              " -0.01976423537605237,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " -0.01976423537605237,\n",
              " 0.0,\n",
              " 0.01976423537605237,\n",
              " 0.0,\n",
              " 0.0,\n",
              " -0.01976423537605237,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.01976423537605237,\n",
              " 0.0,\n",
              " 0.0,\n",
              " -0.01976423537605237,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.03952847075210474,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.01976423537605237,\n",
              " 0.0,\n",
              " 0.03952847075210474,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.01976423537605237,\n",
              " 0.0,\n",
              " 0.03952847075210474,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.03952847075210474,\n",
              " -0.01976423537605237,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.01976423537605237,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.03952847075210474,\n",
              " -0.01976423537605237,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.01976423537605237,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " -0.01976423537605237,\n",
              " 0.0,\n",
              " 0.03952847075210474,\n",
              " -0.01976423537605237,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " -0.01976423537605237,\n",
              " 0.0,\n",
              " 0.03952847075210474,\n",
              " -0.01976423537605237,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.01976423537605237,\n",
              " 0.0,\n",
              " 0.03952847075210474,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.01976423537605237,\n",
              " 0.0,\n",
              " 0.03952847075210474,\n",
              " 0.01976423537605237,\n",
              " 0.03952847075210474,\n",
              " 0.0,\n",
              " -0.01976423537605237,\n",
              " 0.0,\n",
              " -0.01976423537605237,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " -0.01976423537605237,\n",
              " 0.0,\n",
              " 0.01976423537605237,\n",
              " 0.0,\n",
              " 0.0,\n",
              " -0.01976423537605237,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.01976423537605237,\n",
              " 0.0,\n",
              " 0.0,\n",
              " -0.01976423537605237,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.03952847075210474,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.01976423537605237,\n",
              " 0.0,\n",
              " 0.03952847075210474,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.01976423537605237,\n",
              " 0.0,\n",
              " 0.03952847075210474,\n",
              " 0.0,\n",
              " 0.0,\n",
              " -0.01976423537605237,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " -0.01976423537605237,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.01976423537605237,\n",
              " 0.0,\n",
              " 0.0,\n",
              " -0.01976423537605237,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.03952847075210474,\n",
              " 0.0,\n",
              " -0.01976423537605237,\n",
              " -0.03952847075210474,\n",
              " 0.0,\n",
              " 0.01976423537605237,\n",
              " 0.0,\n",
              " 0.01976423537605237,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.01976423537605237,\n",
              " 0.0,\n",
              " -0.01976423537605237,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.01976423537605237,\n",
              " 0.0,\n",
              " 0.0,\n",
              " -0.01976423537605237,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " -0.01976423537605237,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " -0.01976423537605237,\n",
              " 0.0,\n",
              " 0.03952847075210474,\n",
              " -0.01976423537605237,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " -0.01976423537605237,\n",
              " 0.0,\n",
              " 0.03952847075210474,\n",
              " -0.01976423537605237,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " -0.03952847075210474,\n",
              " 0.01976423537605237,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " -0.01976423537605237,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " -0.03952847075210474,\n",
              " 0.01976423537605237,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " -0.01976423537605237,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " -0.01976423537605237,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " -0.01976423537605237,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.01976423537605237,\n",
              " 0.0,\n",
              " 0.0,\n",
              " -0.01976423537605237,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.03952847075210474,\n",
              " 0.0,\n",
              " -0.01976423537605237,\n",
              " -0.03952847075210474,\n",
              " 0.0,\n",
              " 0.01976423537605237,\n",
              " 0.0,\n",
              " 0.01976423537605237,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.01976423537605237,\n",
              " 0.0,\n",
              " -0.01976423537605237,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.01976423537605237,\n",
              " 0.0,\n",
              " 0.0,\n",
              " -0.01976423537605237,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " -0.01976423537605237,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " 0.0,\n",
              " ...]"
            ]
          },
          "execution_count": 6,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "coeffs"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "39d9922d-f94d-447d-9fc7-c55d12c9da19",
      "metadata": {
        "id": "39d9922d-f94d-447d-9fc7-c55d12c9da19"
      },
      "source": [
        "### Block Encoding of the (padded) Hamiltonian matrix"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "df5fe8f5-1b72-4f77-98b5-e551d0567fc8",
      "metadata": {
        "id": "df5fe8f5-1b72-4f77-98b5-e551d0567fc8"
      },
      "outputs": [],
      "source": [
        "import matplotlib.pyplot as plt\n",
        "import numpy as np\n",
        "import scipy\n",
        "\n",
        "from classiq import *\n",
        "from classiq.execution import (\n",
        "    ClassiqBackendPreferences,\n",
        "    ClassiqSimulatorBackendNames,\n",
        "    ExecutionPreferences,\n",
        ")\n",
        "from classiq.qmod.symbolic import pi\n",
        "\n",
        "def get_normalized_lcu_coef(lcu_coef):\n",
        "\n",
        "    normalization_factor = sum(lcu_coef)\n",
        "    prepare_prob = [c / normalization_factor for c in lcu_coef]\n",
        "    coef_size = int(np.ceil(np.log2(len(prepare_prob))))\n",
        "    prepare_prob += [0] * (2**coef_size - len(prepare_prob))\n",
        "\n",
        "    print(\"The size of the block encoding:\", coef_size)\n",
        "    print(\"\\n The normalized coefficients:\", prepare_prob)\n",
        "    print(\"\\n The normalization factor:\", normalization_factor)\n",
        "\n",
        "    return normalization_factor, coef_size, prepare_prob\n",
        "\n",
        "@qfunc\n",
        "def apply_pauli_term(pauli_string: PauliTerm, x: QArray[QBit]):\n",
        "    repeat(\n",
        "        count=x.len,\n",
        "        iteration=lambda index: switch(\n",
        "            pauli_string.pauli[index],\n",
        "            [\n",
        "                lambda: IDENTITY(x[pauli_string.pauli.len - index - 1]),\n",
        "                lambda: X(x[pauli_string.pauli.len - index - 1]),\n",
        "                lambda: Y(x[pauli_string.pauli.len - index - 1]),\n",
        "                lambda: Z(x[pauli_string.pauli.len - index - 1]),\n",
        "            ],\n",
        "        ),\n",
        "    )\n",
        "\n",
        "\n",
        "@qfunc\n",
        "def lcu_paulis(\n",
        "    pauli_terms_list: CArray[PauliTerm],\n",
        "    probs: CArray[CReal],\n",
        "    block: QNum,\n",
        "    data: QArray[QBit],\n",
        "):\n",
        "    within_apply(\n",
        "        lambda: inplace_prepare_state(probs, 0.0, block),\n",
        "        lambda: repeat(\n",
        "            count=pauli_terms_list.len,\n",
        "            iteration=lambda i: control(\n",
        "                block == i, lambda: apply_pauli_term(pauli_terms_list[i], data)\n",
        "            ),\n",
        "        ),\n",
        "    )\n",
        "\n",
        "## fix the execution preferences for this tutorial\n",
        "Execution_Prefs = ExecutionPreferences(\n",
        "    num_shots=1,\n",
        "    backend_preferences=ClassiqBackendPreferences(\n",
        "        backend_name=ClassiqSimulatorBackendNames.SIMULATOR_STATEVECTOR\n",
        "    ),\n",
        ")\n",
        "\n",
        "def get_projected_state_vector(\n",
        "    execution_result,\n",
        "    measured_var: str,\n",
        "    projections: dict,\n",
        ") -> np.ndarray:\n",
        "    \"\"\"\n",
        "    This function returns a reduced statevector from execution results.\n",
        "    measured var: the name of the reduced variable\n",
        "    projections: on which values of the other variables to project, e.g., {\"ind\": 1}\n",
        "    \"\"\"\n",
        "    projected_size = len(execution_result[0].value.output_qubits_map[measured_var])\n",
        "    proj_statevector = np.zeros(2**projected_size).astype(complex)\n",
        "    for sample in execution_result[0].value.parsed_state_vector:\n",
        "        if all(sample.state[key] == projections[key] for key in projections.keys()):\n",
        "            proj_statevector[int(sample.state[measured_var])] += sample.amplitude\n",
        "    return proj_statevector"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "ad209f08-5f6c-4ba2-ac66-e5dd81f48271",
      "metadata": {
        "id": "ad209f08-5f6c-4ba2-ac66-e5dd81f48271"
      },
      "outputs": [],
      "source": [
        "\"\"\"\n",
        "Generating the HAMILTONIAN variable that is input to the function lcu_paulis() connecting our LCU decomposition\n",
        "to the next steps of block encoding\n",
        "\"\"\"\n",
        "HAMILTONIAN = [\n",
        "    PauliTerm(pauli=[getattr(Pauli, label) for label in term], coefficient=coeff)\n",
        "    for term, coeff in zip(decomposition.keys(), coeffs) if not np.isclose(coeff, 0)\n",
        "]"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "7a080854-70ab-4ae6-8974-54b44ea38690",
      "metadata": {
        "id": "7a080854-70ab-4ae6-8974-54b44ea38690",
        "outputId": "4112727f-f619-4c57-c2d0-b5b26fe6504c"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "The size of the block encoding: 7\n",
            "\n",
            " The normalized coefficients: [0.008928571428571444, 0.017857142857142887, 0.008928571428571444, 0.017857142857142887, 0.008928571428571444, 0.017857142857142887, 0.008928571428571444, 0.008928571428571444, 0.008928571428571444, 0.008928571428571444, 0.008928571428571444, 0.008928571428571444, 0.008928571428571444, 0.017857142857142887, 0.008928571428571444, 0.017857142857142887, 0.008928571428571444, 0.017857142857142887, 0.017857142857142887, 0.008928571428571444, 0.008928571428571444, 0.017857142857142887, 0.008928571428571444, 0.008928571428571444, 0.008928571428571444, 0.017857142857142887, 0.008928571428571444, 0.008928571428571444, 0.017857142857142887, 0.008928571428571444, 0.008928571428571444, 0.017857142857142887, 0.008928571428571444, 0.017857142857142887, 0.008928571428571444, 0.017857142857142887, 0.008928571428571444, 0.008928571428571444, 0.008928571428571444, 0.008928571428571444, 0.008928571428571444, 0.008928571428571444, 0.008928571428571444, 0.017857142857142887, 0.008928571428571444, 0.017857142857142887, 0.008928571428571444, 0.017857142857142887, 0.008928571428571444, 0.008928571428571444, 0.008928571428571444, 0.008928571428571444, 0.017857142857142887, 0.008928571428571444, 0.017857142857142887, 0.008928571428571444, 0.008928571428571444, 0.008928571428571444, 0.008928571428571444, 0.008928571428571444, 0.008928571428571444, 0.008928571428571444, 0.008928571428571444, 0.017857142857142887, 0.008928571428571444, 0.008928571428571444, 0.017857142857142887, 0.008928571428571444, 0.017857142857142887, 0.008928571428571444, 0.008928571428571444, 0.017857142857142887, 0.008928571428571444, 0.008928571428571444, 0.008928571428571444, 0.008928571428571444, 0.008928571428571444, 0.008928571428571444, 0.017857142857142887, 0.008928571428571444, 0.017857142857142887, 0.008928571428571444, 0.008928571428571444, 0.008928571428571444, 0.008928571428571444, 0.008928571428571444, 0.008928571428571444, 0.008928571428571444, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]\n",
            "\n",
            " The normalization factor: 2.213594362117862\n"
          ]
        }
      ],
      "source": [
        "\"\"\"\n",
        "For the case of the two oscillators N=2, we also construct the HAMILTONIAN matrix manually to compare it with the automated construction\n",
        "of the previous cell and verify that the previous code is correct\n",
        "\"\"\"\n",
        "HAMILTONIAN_manually = [\n",
        "    PauliTerm(pauli=[Pauli.X, Pauli.I, Pauli.I], coefficient=0.25),\n",
        "    PauliTerm(pauli=[Pauli.X, Pauli.I, Pauli.X], coefficient=0.50),\n",
        "    PauliTerm(pauli=[Pauli.X, Pauli.I, Pauli.Z], coefficient=1.25),\n",
        "    PauliTerm(pauli=[Pauli.X, Pauli.X, Pauli.X], coefficient=0.75),\n",
        "    PauliTerm(pauli=[Pauli.X, Pauli.Y, Pauli.Y], coefficient=0.75),\n",
        "    PauliTerm(pauli=[Pauli.X, Pauli.Z, Pauli.I], coefficient=0.25),\n",
        "    PauliTerm(pauli=[Pauli.X, Pauli.Z, Pauli.X], coefficient=0.50),\n",
        "    PauliTerm(pauli=[Pauli.X, Pauli.Z, Pauli.Z], coefficient=1.25),\n",
        "    PauliTerm(pauli=[Pauli.Y, Pauli.I, Pauli.Y], coefficient=-0.50),\n",
        "    PauliTerm(pauli=[Pauli.Y, Pauli.X, Pauli.Y], coefficient=0.75),\n",
        "    PauliTerm(pauli=[Pauli.Y, Pauli.Y, Pauli.X], coefficient=-0.75),\n",
        "    PauliTerm(pauli=[Pauli.Y, Pauli.Z, Pauli.Y], coefficient=-0.50)\n",
        "]\n",
        "\n",
        "\n",
        "\"\"\"\n",
        "We store the coefficients of the expansion in lcu_pauli_coef to imput them in the get_normalized_lcu_coef() function and\n",
        "obtain the normalization factor, the normalized Hamiltonian, and its preparation probabilities.\n",
        "\"\"\"\n",
        "lcu_pauli_coef = [c.real for c in coeffs if c.real!=0]\n",
        "\n",
        "normalization_factor = sum(np.abs(lcu_pauli_coef))\n",
        "lcu_pauli_coef2=[]\n",
        "lcu_pauli_coef2 = [abs(c) / normalization_factor for c in lcu_pauli_coef]\n",
        "\n",
        "normalization_ham, lcu_size_ham, prepare_probs_ham = get_normalized_lcu_coef(\n",
        "    #lcu_pauli_coef\n",
        "    np.abs(lcu_pauli_coef).tolist()\n",
        "    #lcu_pauli_coef2\n",
        ")"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "6bdd8bc8-e6d2-4769-8bbe-893c3c6ee74a",
      "metadata": {
        "id": "6bdd8bc8-e6d2-4769-8bbe-893c3c6ee74a",
        "outputId": "ffcea3f5-d8c7-4817-cb44-34f4b8bb9227"
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "[PauliTerm(pauli=[<Pauli.X: 1>, <Pauli.I: 0>, <Pauli.I: 0>, <Pauli.I: 0>, <Pauli.X: 1>], coefficient=0.01976423537605237),\n",
              " PauliTerm(pauli=[<Pauli.X: 1>, <Pauli.I: 0>, <Pauli.I: 0>, <Pauli.I: 0>, <Pauli.Z: 3>], coefficient=0.03952847075210474),\n",
              " PauliTerm(pauli=[<Pauli.X: 1>, <Pauli.I: 0>, <Pauli.I: 0>, <Pauli.Z: 3>, <Pauli.X: 1>], coefficient=0.01976423537605237),\n",
              " PauliTerm(pauli=[<Pauli.X: 1>, <Pauli.I: 0>, <Pauli.I: 0>, <Pauli.Z: 3>, <Pauli.Z: 3>], coefficient=0.03952847075210474),\n",
              " PauliTerm(pauli=[<Pauli.X: 1>, <Pauli.I: 0>, <Pauli.X: 1>, <Pauli.I: 0>, <Pauli.I: 0>], coefficient=0.01976423537605237),\n",
              " PauliTerm(pauli=[<Pauli.X: 1>, <Pauli.I: 0>, <Pauli.X: 1>, <Pauli.I: 0>, <Pauli.X: 1>], coefficient=0.03952847075210474),\n",
              " PauliTerm(pauli=[<Pauli.X: 1>, <Pauli.I: 0>, <Pauli.X: 1>, <Pauli.I: 0>, <Pauli.Z: 3>], coefficient=-0.01976423537605237),\n",
              " PauliTerm(pauli=[<Pauli.X: 1>, <Pauli.I: 0>, <Pauli.X: 1>, <Pauli.X: 1>, <Pauli.X: 1>], coefficient=-0.01976423537605237),\n",
              " PauliTerm(pauli=[<Pauli.X: 1>, <Pauli.I: 0>, <Pauli.X: 1>, <Pauli.Y: 2>, <Pauli.Y: 2>], coefficient=-0.01976423537605237),\n",
              " PauliTerm(pauli=[<Pauli.X: 1>, <Pauli.I: 0>, <Pauli.X: 1>, <Pauli.Z: 3>, <Pauli.I: 0>], coefficient=0.01976423537605237),\n",
              " PauliTerm(pauli=[<Pauli.X: 1>, <Pauli.I: 0>, <Pauli.X: 1>, <Pauli.Z: 3>, <Pauli.Z: 3>], coefficient=-0.01976423537605237),\n",
              " PauliTerm(pauli=[<Pauli.X: 1>, <Pauli.I: 0>, <Pauli.Y: 2>, <Pauli.X: 1>, <Pauli.Y: 2>], coefficient=0.01976423537605237),\n",
              " PauliTerm(pauli=[<Pauli.X: 1>, <Pauli.I: 0>, <Pauli.Y: 2>, <Pauli.Y: 2>, <Pauli.X: 1>], coefficient=-0.01976423537605237),\n",
              " PauliTerm(pauli=[<Pauli.X: 1>, <Pauli.I: 0>, <Pauli.Y: 2>, <Pauli.Z: 3>, <Pauli.Y: 2>], coefficient=0.03952847075210474),\n",
              " PauliTerm(pauli=[<Pauli.X: 1>, <Pauli.I: 0>, <Pauli.Z: 3>, <Pauli.I: 0>, <Pauli.X: 1>], coefficient=0.01976423537605237),\n",
              " PauliTerm(pauli=[<Pauli.X: 1>, <Pauli.I: 0>, <Pauli.Z: 3>, <Pauli.I: 0>, <Pauli.Z: 3>], coefficient=0.03952847075210474),\n",
              " PauliTerm(pauli=[<Pauli.X: 1>, <Pauli.I: 0>, <Pauli.Z: 3>, <Pauli.Z: 3>, <Pauli.X: 1>], coefficient=0.01976423537605237),\n",
              " PauliTerm(pauli=[<Pauli.X: 1>, <Pauli.I: 0>, <Pauli.Z: 3>, <Pauli.Z: 3>, <Pauli.Z: 3>], coefficient=0.03952847075210474),\n",
              " PauliTerm(pauli=[<Pauli.X: 1>, <Pauli.X: 1>, <Pauli.I: 0>, <Pauli.X: 1>, <Pauli.I: 0>], coefficient=0.03952847075210474),\n",
              " PauliTerm(pauli=[<Pauli.X: 1>, <Pauli.X: 1>, <Pauli.I: 0>, <Pauli.X: 1>, <Pauli.X: 1>], coefficient=-0.01976423537605237),\n",
              " PauliTerm(pauli=[<Pauli.X: 1>, <Pauli.X: 1>, <Pauli.I: 0>, <Pauli.Y: 2>, <Pauli.Y: 2>], coefficient=0.01976423537605237),\n",
              " PauliTerm(pauli=[<Pauli.X: 1>, <Pauli.X: 1>, <Pauli.Z: 3>, <Pauli.X: 1>, <Pauli.I: 0>], coefficient=0.03952847075210474),\n",
              " PauliTerm(pauli=[<Pauli.X: 1>, <Pauli.X: 1>, <Pauli.Z: 3>, <Pauli.X: 1>, <Pauli.X: 1>], coefficient=-0.01976423537605237),\n",
              " PauliTerm(pauli=[<Pauli.X: 1>, <Pauli.X: 1>, <Pauli.Z: 3>, <Pauli.Y: 2>, <Pauli.Y: 2>], coefficient=0.01976423537605237),\n",
              " PauliTerm(pauli=[<Pauli.X: 1>, <Pauli.Y: 2>, <Pauli.I: 0>, <Pauli.X: 1>, <Pauli.Y: 2>], coefficient=-0.01976423537605237),\n",
              " PauliTerm(pauli=[<Pauli.X: 1>, <Pauli.Y: 2>, <Pauli.I: 0>, <Pauli.Y: 2>, <Pauli.I: 0>], coefficient=0.03952847075210474),\n",
              " PauliTerm(pauli=[<Pauli.X: 1>, <Pauli.Y: 2>, <Pauli.I: 0>, <Pauli.Y: 2>, <Pauli.X: 1>], coefficient=-0.01976423537605237),\n",
              " PauliTerm(pauli=[<Pauli.X: 1>, <Pauli.Y: 2>, <Pauli.Z: 3>, <Pauli.X: 1>, <Pauli.Y: 2>], coefficient=-0.01976423537605237),\n",
              " PauliTerm(pauli=[<Pauli.X: 1>, <Pauli.Y: 2>, <Pauli.Z: 3>, <Pauli.Y: 2>, <Pauli.I: 0>], coefficient=0.03952847075210474),\n",
              " PauliTerm(pauli=[<Pauli.X: 1>, <Pauli.Y: 2>, <Pauli.Z: 3>, <Pauli.Y: 2>, <Pauli.X: 1>], coefficient=-0.01976423537605237),\n",
              " PauliTerm(pauli=[<Pauli.X: 1>, <Pauli.Z: 3>, <Pauli.I: 0>, <Pauli.I: 0>, <Pauli.X: 1>], coefficient=0.01976423537605237),\n",
              " PauliTerm(pauli=[<Pauli.X: 1>, <Pauli.Z: 3>, <Pauli.I: 0>, <Pauli.I: 0>, <Pauli.Z: 3>], coefficient=0.03952847075210474),\n",
              " PauliTerm(pauli=[<Pauli.X: 1>, <Pauli.Z: 3>, <Pauli.I: 0>, <Pauli.Z: 3>, <Pauli.X: 1>], coefficient=0.01976423537605237),\n",
              " PauliTerm(pauli=[<Pauli.X: 1>, <Pauli.Z: 3>, <Pauli.I: 0>, <Pauli.Z: 3>, <Pauli.Z: 3>], coefficient=0.03952847075210474),\n",
              " PauliTerm(pauli=[<Pauli.X: 1>, <Pauli.Z: 3>, <Pauli.X: 1>, <Pauli.I: 0>, <Pauli.I: 0>], coefficient=0.01976423537605237),\n",
              " PauliTerm(pauli=[<Pauli.X: 1>, <Pauli.Z: 3>, <Pauli.X: 1>, <Pauli.I: 0>, <Pauli.X: 1>], coefficient=0.03952847075210474),\n",
              " PauliTerm(pauli=[<Pauli.X: 1>, <Pauli.Z: 3>, <Pauli.X: 1>, <Pauli.I: 0>, <Pauli.Z: 3>], coefficient=-0.01976423537605237),\n",
              " PauliTerm(pauli=[<Pauli.X: 1>, <Pauli.Z: 3>, <Pauli.X: 1>, <Pauli.X: 1>, <Pauli.X: 1>], coefficient=-0.01976423537605237),\n",
              " PauliTerm(pauli=[<Pauli.X: 1>, <Pauli.Z: 3>, <Pauli.X: 1>, <Pauli.Y: 2>, <Pauli.Y: 2>], coefficient=-0.01976423537605237),\n",
              " PauliTerm(pauli=[<Pauli.X: 1>, <Pauli.Z: 3>, <Pauli.X: 1>, <Pauli.Z: 3>, <Pauli.I: 0>], coefficient=0.01976423537605237),\n",
              " PauliTerm(pauli=[<Pauli.X: 1>, <Pauli.Z: 3>, <Pauli.X: 1>, <Pauli.Z: 3>, <Pauli.Z: 3>], coefficient=-0.01976423537605237),\n",
              " PauliTerm(pauli=[<Pauli.X: 1>, <Pauli.Z: 3>, <Pauli.Y: 2>, <Pauli.X: 1>, <Pauli.Y: 2>], coefficient=0.01976423537605237),\n",
              " PauliTerm(pauli=[<Pauli.X: 1>, <Pauli.Z: 3>, <Pauli.Y: 2>, <Pauli.Y: 2>, <Pauli.X: 1>], coefficient=-0.01976423537605237),\n",
              " PauliTerm(pauli=[<Pauli.X: 1>, <Pauli.Z: 3>, <Pauli.Y: 2>, <Pauli.Z: 3>, <Pauli.Y: 2>], coefficient=0.03952847075210474),\n",
              " PauliTerm(pauli=[<Pauli.X: 1>, <Pauli.Z: 3>, <Pauli.Z: 3>, <Pauli.I: 0>, <Pauli.X: 1>], coefficient=0.01976423537605237),\n",
              " PauliTerm(pauli=[<Pauli.X: 1>, <Pauli.Z: 3>, <Pauli.Z: 3>, <Pauli.I: 0>, <Pauli.Z: 3>], coefficient=0.03952847075210474),\n",
              " PauliTerm(pauli=[<Pauli.X: 1>, <Pauli.Z: 3>, <Pauli.Z: 3>, <Pauli.Z: 3>, <Pauli.X: 1>], coefficient=0.01976423537605237),\n",
              " PauliTerm(pauli=[<Pauli.X: 1>, <Pauli.Z: 3>, <Pauli.Z: 3>, <Pauli.Z: 3>, <Pauli.Z: 3>], coefficient=0.03952847075210474),\n",
              " PauliTerm(pauli=[<Pauli.Y: 2>, <Pauli.I: 0>, <Pauli.I: 0>, <Pauli.I: 0>, <Pauli.Y: 2>], coefficient=-0.01976423537605237),\n",
              " PauliTerm(pauli=[<Pauli.Y: 2>, <Pauli.I: 0>, <Pauli.I: 0>, <Pauli.Z: 3>, <Pauli.Y: 2>], coefficient=-0.01976423537605237),\n",
              " PauliTerm(pauli=[<Pauli.Y: 2>, <Pauli.I: 0>, <Pauli.X: 1>, <Pauli.X: 1>, <Pauli.Y: 2>], coefficient=0.01976423537605237),\n",
              " PauliTerm(pauli=[<Pauli.Y: 2>, <Pauli.I: 0>, <Pauli.X: 1>, <Pauli.Y: 2>, <Pauli.X: 1>], coefficient=-0.01976423537605237),\n",
              " PauliTerm(pauli=[<Pauli.Y: 2>, <Pauli.I: 0>, <Pauli.X: 1>, <Pauli.Z: 3>, <Pauli.Y: 2>], coefficient=0.03952847075210474),\n",
              " PauliTerm(pauli=[<Pauli.Y: 2>, <Pauli.I: 0>, <Pauli.Y: 2>, <Pauli.I: 0>, <Pauli.I: 0>], coefficient=-0.01976423537605237),\n",
              " PauliTerm(pauli=[<Pauli.Y: 2>, <Pauli.I: 0>, <Pauli.Y: 2>, <Pauli.I: 0>, <Pauli.X: 1>], coefficient=-0.03952847075210474),\n",
              " PauliTerm(pauli=[<Pauli.Y: 2>, <Pauli.I: 0>, <Pauli.Y: 2>, <Pauli.I: 0>, <Pauli.Z: 3>], coefficient=0.01976423537605237),\n",
              " PauliTerm(pauli=[<Pauli.Y: 2>, <Pauli.I: 0>, <Pauli.Y: 2>, <Pauli.X: 1>, <Pauli.X: 1>], coefficient=0.01976423537605237),\n",
              " PauliTerm(pauli=[<Pauli.Y: 2>, <Pauli.I: 0>, <Pauli.Y: 2>, <Pauli.Y: 2>, <Pauli.Y: 2>], coefficient=0.01976423537605237),\n",
              " PauliTerm(pauli=[<Pauli.Y: 2>, <Pauli.I: 0>, <Pauli.Y: 2>, <Pauli.Z: 3>, <Pauli.I: 0>], coefficient=-0.01976423537605237),\n",
              " PauliTerm(pauli=[<Pauli.Y: 2>, <Pauli.I: 0>, <Pauli.Y: 2>, <Pauli.Z: 3>, <Pauli.Z: 3>], coefficient=0.01976423537605237),\n",
              " PauliTerm(pauli=[<Pauli.Y: 2>, <Pauli.I: 0>, <Pauli.Z: 3>, <Pauli.I: 0>, <Pauli.Y: 2>], coefficient=-0.01976423537605237),\n",
              " PauliTerm(pauli=[<Pauli.Y: 2>, <Pauli.I: 0>, <Pauli.Z: 3>, <Pauli.Z: 3>, <Pauli.Y: 2>], coefficient=-0.01976423537605237),\n",
              " PauliTerm(pauli=[<Pauli.Y: 2>, <Pauli.X: 1>, <Pauli.I: 0>, <Pauli.X: 1>, <Pauli.Y: 2>], coefficient=-0.01976423537605237),\n",
              " PauliTerm(pauli=[<Pauli.Y: 2>, <Pauli.X: 1>, <Pauli.I: 0>, <Pauli.Y: 2>, <Pauli.I: 0>], coefficient=0.03952847075210474),\n",
              " PauliTerm(pauli=[<Pauli.Y: 2>, <Pauli.X: 1>, <Pauli.I: 0>, <Pauli.Y: 2>, <Pauli.X: 1>], coefficient=-0.01976423537605237),\n",
              " PauliTerm(pauli=[<Pauli.Y: 2>, <Pauli.X: 1>, <Pauli.Z: 3>, <Pauli.X: 1>, <Pauli.Y: 2>], coefficient=-0.01976423537605237),\n",
              " PauliTerm(pauli=[<Pauli.Y: 2>, <Pauli.X: 1>, <Pauli.Z: 3>, <Pauli.Y: 2>, <Pauli.I: 0>], coefficient=0.03952847075210474),\n",
              " PauliTerm(pauli=[<Pauli.Y: 2>, <Pauli.X: 1>, <Pauli.Z: 3>, <Pauli.Y: 2>, <Pauli.X: 1>], coefficient=-0.01976423537605237),\n",
              " PauliTerm(pauli=[<Pauli.Y: 2>, <Pauli.Y: 2>, <Pauli.I: 0>, <Pauli.X: 1>, <Pauli.I: 0>], coefficient=-0.03952847075210474),\n",
              " PauliTerm(pauli=[<Pauli.Y: 2>, <Pauli.Y: 2>, <Pauli.I: 0>, <Pauli.X: 1>, <Pauli.X: 1>], coefficient=0.01976423537605237),\n",
              " PauliTerm(pauli=[<Pauli.Y: 2>, <Pauli.Y: 2>, <Pauli.I: 0>, <Pauli.Y: 2>, <Pauli.Y: 2>], coefficient=-0.01976423537605237),\n",
              " PauliTerm(pauli=[<Pauli.Y: 2>, <Pauli.Y: 2>, <Pauli.Z: 3>, <Pauli.X: 1>, <Pauli.I: 0>], coefficient=-0.03952847075210474),\n",
              " PauliTerm(pauli=[<Pauli.Y: 2>, <Pauli.Y: 2>, <Pauli.Z: 3>, <Pauli.X: 1>, <Pauli.X: 1>], coefficient=0.01976423537605237),\n",
              " PauliTerm(pauli=[<Pauli.Y: 2>, <Pauli.Y: 2>, <Pauli.Z: 3>, <Pauli.Y: 2>, <Pauli.Y: 2>], coefficient=-0.01976423537605237),\n",
              " PauliTerm(pauli=[<Pauli.Y: 2>, <Pauli.Z: 3>, <Pauli.I: 0>, <Pauli.I: 0>, <Pauli.Y: 2>], coefficient=-0.01976423537605237),\n",
              " PauliTerm(pauli=[<Pauli.Y: 2>, <Pauli.Z: 3>, <Pauli.I: 0>, <Pauli.Z: 3>, <Pauli.Y: 2>], coefficient=-0.01976423537605237),\n",
              " PauliTerm(pauli=[<Pauli.Y: 2>, <Pauli.Z: 3>, <Pauli.X: 1>, <Pauli.X: 1>, <Pauli.Y: 2>], coefficient=0.01976423537605237),\n",
              " PauliTerm(pauli=[<Pauli.Y: 2>, <Pauli.Z: 3>, <Pauli.X: 1>, <Pauli.Y: 2>, <Pauli.X: 1>], coefficient=-0.01976423537605237),\n",
              " PauliTerm(pauli=[<Pauli.Y: 2>, <Pauli.Z: 3>, <Pauli.X: 1>, <Pauli.Z: 3>, <Pauli.Y: 2>], coefficient=0.03952847075210474),\n",
              " PauliTerm(pauli=[<Pauli.Y: 2>, <Pauli.Z: 3>, <Pauli.Y: 2>, <Pauli.I: 0>, <Pauli.I: 0>], coefficient=-0.01976423537605237),\n",
              " PauliTerm(pauli=[<Pauli.Y: 2>, <Pauli.Z: 3>, <Pauli.Y: 2>, <Pauli.I: 0>, <Pauli.X: 1>], coefficient=-0.03952847075210474),\n",
              " PauliTerm(pauli=[<Pauli.Y: 2>, <Pauli.Z: 3>, <Pauli.Y: 2>, <Pauli.I: 0>, <Pauli.Z: 3>], coefficient=0.01976423537605237),\n",
              " PauliTerm(pauli=[<Pauli.Y: 2>, <Pauli.Z: 3>, <Pauli.Y: 2>, <Pauli.X: 1>, <Pauli.X: 1>], coefficient=0.01976423537605237),\n",
              " PauliTerm(pauli=[<Pauli.Y: 2>, <Pauli.Z: 3>, <Pauli.Y: 2>, <Pauli.Y: 2>, <Pauli.Y: 2>], coefficient=0.01976423537605237),\n",
              " PauliTerm(pauli=[<Pauli.Y: 2>, <Pauli.Z: 3>, <Pauli.Y: 2>, <Pauli.Z: 3>, <Pauli.I: 0>], coefficient=-0.01976423537605237),\n",
              " PauliTerm(pauli=[<Pauli.Y: 2>, <Pauli.Z: 3>, <Pauli.Y: 2>, <Pauli.Z: 3>, <Pauli.Z: 3>], coefficient=0.01976423537605237),\n",
              " PauliTerm(pauli=[<Pauli.Y: 2>, <Pauli.Z: 3>, <Pauli.Z: 3>, <Pauli.I: 0>, <Pauli.Y: 2>], coefficient=-0.01976423537605237),\n",
              " PauliTerm(pauli=[<Pauli.Y: 2>, <Pauli.Z: 3>, <Pauli.Z: 3>, <Pauli.Z: 3>, <Pauli.Y: 2>], coefficient=-0.01976423537605237)]"
            ]
          },
          "execution_count": 12,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "HAMILTONIAN"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "3d5cc028-c214-4094-b04b-0afa35dcebfb",
      "metadata": {
        "id": "3d5cc028-c214-4094-b04b-0afa35dcebfb"
      },
      "source": [
        "### State preparation\n",
        "\n",
        "Now, we construct the initial state:\n",
        "\n",
        "$\n",
        "\\ket{\\psi_o}=\n",
        "\\begin{pmatrix}\n",
        "\\dot{\\vec{y}}(0)\\\\\n",
        "i\\textbf{B}^\\dagger\\vec{y}(0)\n",
        "\\end{pmatrix}, \\quad (22)\n",
        "$\n",
        "\n",
        "Note that the vector (22) is $(N+M)$-dimensional where $M=N*(N+1)/2$. We want to work with the power of 2 dimensions. This can be achieved by working with the padded $\\bar{\\textbf{B}}$ matrix. Then, the resulting initial quantum state where the initial positions and velocities are encoded is:\n",
        "\n",
        "$\n",
        "\\ket{\\psi_o}=\n",
        "\\begin{pmatrix}\n",
        "\\dot{\\vec{y}}(0)\\\\\n",
        "i\\bar{\\textbf{B}}^\\dagger\\vec{y}(0)\n",
        "\\end{pmatrix}, \\quad (23)\n",
        "$\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "5c7f6066-2552-465a-bb39-7a0f11b99b5d",
      "metadata": {
        "id": "5c7f6066-2552-465a-bb39-7a0f11b99b5d",
        "outputId": "1ab1c976-5b75-4925-bea2-24070c5ee274"
      },
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "ydot_pad_o:\n",
            "[ 0.65513657 -0.9809765  -0.80818839 -0.00617164  0.          0.\n",
            "  0.          0.          0.          0.          0.          0.\n",
            "  0.          0.          0.          0.        ]\n",
            "\n",
            "y_pad_o:\n",
            "[ 0.22867976  0.1912002  -0.6153515   0.87272955  0.          0.\n",
            "  0.          0.          0.          0.          0.          0.\n",
            "  0.          0.          0.          0.        ]\n",
            "\n",
            "psi_pad_o:\n",
            "[ 0.65513657+0.j         -0.9809765 +0.j         -0.80818839+0.j\n",
            " -0.00617164+0.j          0.        +0.j          0.        +0.j\n",
            "  0.        +0.j          0.        +0.j          0.        +0.j\n",
            "  0.        +0.j          0.        +0.j          0.        +0.j\n",
            "  0.        +0.j          0.        +0.j          0.        +0.j\n",
            "  0.        +0.j          0.        +0.07231489j  0.        +0.01185208j\n",
            "  0.        +0.j          0.        +0.j          0.        +0.06046281j\n",
            "  0.        +0.25505404j  0.        +0.j         -0.        -0.19459123j\n",
            " -0.        -0.47057255j  0.        +0.27598132j  0.        +0.j\n",
            "  0.        +0.j          0.        +0.j          0.        +0.j\n",
            "  0.        +0.j          0.        +0.j        ]\n",
            "\n",
            "Initial amplitudes without imaginary phase (amps):\n",
            "[ 0.41818684 -0.62617701 -0.51588289 -0.00393948  0.          0.\n",
            "  0.          0.          0.          0.          0.          0.\n",
            "  0.          0.          0.          0.          0.04616005  0.00756542\n",
            "  0.          0.          0.03859463  0.16280612  0.         -0.12421149\n",
            " -0.30037591  0.17616442  0.          0.          0.          0.\n",
            "  0.          0.        ]\n",
            "\n",
            "U_matrix:\n",
            "[[1.+0.j 0.+0.j 0.+0.j ... 0.+0.j 0.+0.j 0.+0.j]\n",
            " [0.+0.j 1.+0.j 0.+0.j ... 0.+0.j 0.+0.j 0.+0.j]\n",
            " [0.+0.j 0.+0.j 1.+0.j ... 0.+0.j 0.+0.j 0.+0.j]\n",
            " ...\n",
            " [0.+0.j 0.+0.j 0.+0.j ... 0.+1.j 0.+0.j 0.+0.j]\n",
            " [0.+0.j 0.+0.j 0.+0.j ... 0.+0.j 0.+1.j 0.+0.j]\n",
            " [0.+0.j 0.+0.j 0.+0.j ... 0.+0.j 0.+0.j 0.+1.j]]\n"
          ]
        }
      ],
      "source": [
        "import numpy as np\n",
        "\n",
        "\"\"\"\n",
        "Generate random N-dimensional lists ydot_o and y_o with elements between [-1, 1]\n",
        "\"\"\"\n",
        "ydot_o = np.random.uniform(-1, 1, N)\n",
        "y_o = np.random.uniform(-1, 1, N)\n",
        "\n",
        "\n",
        "\"\"\"\n",
        "Extend ydot_o and y_o to form ydot_pad_o and y_pad_o, respectively, with dimension N^2\n",
        "\"\"\"\n",
        "ydot_pad_o = np.pad(ydot_o, (0, N**2 - N), mode='constant', constant_values=0)\n",
        "y_pad_o = np.pad(y_o, (0, N**2 - N), mode='constant', constant_values=0)\n",
        "\n",
        "\"\"\"\n",
        "Calculate jB_pad^\\dagger y_pad_o (with j = 1j representing the imaginary unit)\n",
        "\"\"\"\n",
        "B_pad_dagger_y_pad_o = 1j * np.dot(B_pad.conj().T, y_pad_o)\n",
        "\n",
        "\"\"\"\n",
        "Construct the column vector psi_pad_o\n",
        "\"\"\"\n",
        "psi_pad_o = np.concatenate((ydot_pad_o, B_pad_dagger_y_pad_o))\n",
        "\n",
        "\"\"\"\n",
        "Store the real part of psi_pad_o to the variable amps\n",
        "\"\"\"\n",
        "#amps = np.real(psi_pad_o)\n",
        "amps = np.concatenate((np.real(psi_pad_o[:N**2]), np.imag(psi_pad_o[N**2:])))\n",
        "norm=np.linalg.norm(amps)\n",
        "amps = amps/np.linalg.norm(amps)\n",
        "\n",
        "\"\"\"\n",
        "Define the matrix U_matrix\n",
        "\"\"\"\n",
        "identity = np.eye(N**2)\n",
        "identity_j = 1j * np.eye(N**2)\n",
        "U_matrix = np.block([\n",
        "    [identity, np.zeros((N**2, N**2))],\n",
        "    [np.zeros((N**2, N**2)), identity_j]\n",
        "])\n",
        "\n",
        "\"\"\"\n",
        "Print the results\n",
        "\"\"\"\n",
        "print(\"ydot_pad_o:\")\n",
        "print(ydot_pad_o)\n",
        "print(\"\\ny_pad_o:\")\n",
        "print(y_pad_o)\n",
        "print(\"\\npsi_pad_o:\")\n",
        "print(psi_pad_o)\n",
        "print(\"\\nInitial amplitudes without imaginary phase (amps):\")\n",
        "print(amps)\n",
        "print(\"\\nU_matrix:\")\n",
        "print(U_matrix)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "a4cf0f2e-252f-484d-8fe1-51ec10c42810",
      "metadata": {
        "id": "a4cf0f2e-252f-484d-8fe1-51ec10c42810",
        "outputId": "648e0957-862f-42e9-e941-bf3587353289"
      },
      "outputs": [
        {
          "data": {
            "text/plain": [
              "array([ 0.41818684, -0.62617701, -0.51588289, -0.00393948,  0.        ,\n",
              "        0.        ,  0.        ,  0.        ,  0.        ,  0.        ,\n",
              "        0.        ,  0.        ,  0.        ,  0.        ,  0.        ,\n",
              "        0.        ,  0.04616005,  0.00756542,  0.        ,  0.        ,\n",
              "        0.03859463,  0.16280612,  0.        , -0.12421149, -0.30037591,\n",
              "        0.17616442,  0.        ,  0.        ,  0.        ,  0.        ,\n",
              "        0.        ,  0.        ])"
            ]
          },
          "execution_count": 15,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "amps"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "27a86cea-ecee-4455-8da7-cce9b2c47ac3",
      "metadata": {
        "id": "27a86cea-ecee-4455-8da7-cce9b2c47ac3"
      },
      "outputs": [],
      "source": [
        "\"\"\"\n",
        "Finally, we are ready to define the function that creates and prepares the initial state in the actual quantum algorithm in the qmod environment.\n",
        "\"\"\"\n",
        "@qfunc\n",
        "def my_state_prep(x: Output[QArray[QBit]]):\n",
        "    prepare_amplitudes(amplitudes=amps.tolist(), bound=0.01, out=x)\n",
        "    unitary(U_matrix.tolist(), x)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "4b952cb2-70f7-433c-b00b-b5ce3640e12c",
      "metadata": {
        "id": "4b952cb2-70f7-433c-b00b-b5ce3640e12c"
      },
      "outputs": [],
      "source": [
        "@qfunc\n",
        "def main(data: Output[QNum], block: Output[QNum]):\n",
        "    allocate(lcu_size_ham, block)\n",
        "    my_state_prep(data)\n",
        "    lcu_paulis(HAMILTONIAN, prepare_probs_ham, block, data)\n",
        "    #lcu_paulis(HAMILTONIAN, np.abs(prepare_probs_ham).tolist(), block, data)\n",
        "\n",
        "\n",
        "qmod = create_model(main, execution_preferences=Execution_Prefs)\n",
        "qprog = synthesize(qmod)\n",
        "show(qprog)\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "2b356731-0ff0-47ae-875e-4abb66c92bf7",
      "metadata": {
        "id": "2b356731-0ff0-47ae-875e-4abb66c92bf7"
      },
      "outputs": [],
      "source": [
        "res = execute(qprog).result()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "3da5c798-4a9f-4ac5-b143-ff69fad02350",
      "metadata": {
        "id": "3da5c798-4a9f-4ac5-b143-ff69fad02350"
      },
      "outputs": [],
      "source": [
        "state_result = get_projected_state_vector(res, \"data\", {\"block\": 0.0})"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "b92b9836-3b7a-4ba0-8b84-490c34bb35d1",
      "metadata": {
        "id": "b92b9836-3b7a-4ba0-8b84-490c34bb35d1"
      },
      "outputs": [],
      "source": [
        "PAULI_MATRICES_DICT = {\n",
        "    Pauli.I: np.array([[1, 0], [0, 1]], dtype=np.complex128),\n",
        "    Pauli.Z: np.array([[1, 0], [0, -1]], dtype=np.complex128),\n",
        "    Pauli.X: np.array([[0, 1], [1, 0]], dtype=np.complex128),\n",
        "    Pauli.Y: np.array([[0, -1j], [1j, 0]], dtype=np.complex128),\n",
        "}\n",
        "\n",
        "\n",
        "def pauli_list_to_mat(pauli_list: list) -> np.ndarray:\n",
        "    real_matrix = 0\n",
        "    for term in pauli_list:\n",
        "        single_str = PAULI_MATRICES_DICT[term.pauli[0]]\n",
        "        for pauli in term.pauli[1:]:\n",
        "            single_str = np.kron(single_str, PAULI_MATRICES_DICT[pauli])\n",
        "        real_matrix += term.coefficient * single_str\n",
        "\n",
        "    assert np.allclose(\n",
        "        np.transpose(np.conjugate(real_matrix)), real_matrix\n",
        "    ), \"matrix not Hermitian\"\n",
        "    assert np.allclose(np.imag(real_matrix), 0.0), \"matrix is not real-valued\"\n",
        "    return np.real(real_matrix)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "65445112-3ff4-4ff9-a10a-2ae7ef706126",
      "metadata": {
        "id": "65445112-3ff4-4ff9-a10a-2ae7ef706126"
      },
      "outputs": [],
      "source": [
        "matrix = pauli_list_to_mat(HAMILTONIAN)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "8eacb79c-0dc2-437a-b7b3-3688c880db4a",
      "metadata": {
        "id": "8eacb79c-0dc2-437a-b7b3-3688c880db4a"
      },
      "outputs": [],
      "source": [
        "normalization_ham"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "5e7c4ada-89fc-4676-b5ea-06488adbc8d9",
      "metadata": {
        "id": "5e7c4ada-89fc-4676-b5ea-06488adbc8d9"
      },
      "outputs": [],
      "source": [
        "matrix"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "0ad06c36-f54d-43a2-869b-050050a30144",
      "metadata": {
        "id": "0ad06c36-f54d-43a2-869b-050050a30144"
      },
      "outputs": [],
      "source": [
        "expected_state = (matrix /  normalization_ham) @ (U_matrix @ amps)\n",
        "relative_phase = np.angle(expected_state[0] / state_result[0])\n",
        "state_result = state_result * np.exp(1j * relative_phase)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "c7c8a055-7b28-4834-b7f8-534e529159cf",
      "metadata": {
        "id": "c7c8a055-7b28-4834-b7f8-534e529159cf"
      },
      "outputs": [],
      "source": [
        "print(\"The resulting state:\", state_result)\n",
        "print(\"The expected state:\", expected_state)\n",
        "assert (np.allclose(state_result, expected_state) < 1e-16)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "012b5610-5a2c-4afb-b668-292506dc7952",
      "metadata": {
        "id": "012b5610-5a2c-4afb-b668-292506dc7952"
      },
      "source": [
        "#### Success!!!"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "96eb48b2-2fe0-4fed-a87f-1d459555fda9",
      "metadata": {
        "id": "96eb48b2-2fe0-4fed-a87f-1d459555fda9"
      },
      "outputs": [],
      "source": [
        "EVOLUTION_TIME = 10\n",
        "EPS = 0.1    # Precision"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "e2771293-a4e8-4d52-939e-f547e2654831",
      "metadata": {
        "id": "e2771293-a4e8-4d52-939e-f547e2654831"
      },
      "outputs": [],
      "source": [
        "@qfunc\n",
        "def my_walk_operator(block: QArray[QBit], data: QArray[QBit]) -> None:\n",
        "    lcu_paulis(HAMILTONIAN, prepare_probs_ham, block, data)  # Unitary implementing H\n",
        "    reflect_about_zero(block) # Π_|0>M operator\n",
        "    RY(2 * pi, block[0])      # for the negative sign found in the unitary construction of W"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "d19d6c68-d3e8-4c28-8fc4-5531a2fbb0ee",
      "metadata": {
        "id": "d19d6c68-d3e8-4c28-8fc4-5531a2fbb0ee"
      },
      "outputs": [],
      "source": [
        "from scipy.special import eval_chebyt, jv\n",
        "\n",
        "def get_cheb_coef(epsilon, t):\n",
        "    poly_degree = int(\n",
        "        np.ceil(\n",
        "            t\n",
        "            + np.log(epsilon ** (-1)) / np.log(np.exp(1) + np.log(epsilon ** (-1)) / t)\n",
        "        )\n",
        "    )\n",
        "    cos_coef = [jv(0, t)] + [\n",
        "        2 * jv(2 * k, t) * (-1) ** k for k in range(1, poly_degree // 2 + 1)\n",
        "    ]\n",
        "    sin_coef = [\n",
        "        -2 * jv(2 * k - 1, t) * (-1) ** k for k in range(1, poly_degree // 2 + 1)\n",
        "    ]\n",
        "    return cos_coef, sin_coef"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "da1071e5-a22d-4a0e-8bfa-d8883431cf28",
      "metadata": {
        "id": "da1071e5-a22d-4a0e-8bfa-d8883431cf28"
      },
      "outputs": [],
      "source": [
        "normalized_time = normalization_ham * EVOLUTION_TIME   # because block encoding for H stores H/alpha_bar (normalization constant)\n",
        "\n",
        "cos_coef, sin_coef = get_cheb_coef(EPS, normalized_time)\n",
        "\n",
        "combined_sin_cos_coef = []\n",
        "for k in range(len(cos_coef) - 1):\n",
        "    combined_sin_cos_coef.append(cos_coef[k])\n",
        "    combined_sin_cos_coef.append(sin_coef[k])\n",
        "combined_sin_cos_coef.append(cos_coef[-1])\n",
        "if len(sin_coef) == len(cos_coef):\n",
        "    combined_sin_cos_coef.append(sin_coef[-1])"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "55ff76f0-9222-460a-80d2-642350a6f488",
      "metadata": {
        "id": "55ff76f0-9222-460a-80d2-642350a6f488"
      },
      "outputs": [],
      "source": [
        "signs_cheb_coef = np.sign(combined_sin_cos_coef).tolist()\n",
        "generalized_signs = [\n",
        "    (1 - signs_cheb_coef[s]) + (s) % 2 for s in range(len(signs_cheb_coef))\n",
        "]\n",
        "positive_cheb_lcu_coef = np.abs(combined_sin_cos_coef)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "940ed3b3-baea-4ce2-9246-bedd51ec303c",
      "metadata": {
        "id": "940ed3b3-baea-4ce2-9246-bedd51ec303c"
      },
      "outputs": [],
      "source": [
        "normalization_exp, lcu_size_exp, prepare_probs_exp = get_normalized_lcu_coef(\n",
        "    positive_cheb_lcu_coef\n",
        ")"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "f01f8283-f5e3-434f-bacf-6d7d9d07c47f",
      "metadata": {
        "id": "f01f8283-f5e3-434f-bacf-6d7d9d07c47f"
      },
      "outputs": [],
      "source": [
        "@qfunc\n",
        "def lcu_cheb(\n",
        "    coef: CArray[CReal],\n",
        "    generalized_signs: CArray[CInt],\n",
        "    walk_operator: QCallable[QNum, QArray],\n",
        "    walk_block: QNum,\n",
        "    walk_data: QArray,\n",
        "    cheb_block: QNum,\n",
        "):\n",
        "\n",
        "    within_apply(\n",
        "        lambda: inplace_prepare_state(coef, 0.0, cheb_block),\n",
        "        lambda: repeat(\n",
        "            generalized_signs.len,\n",
        "            lambda k: control(\n",
        "                cheb_block == k,\n",
        "                lambda: (\n",
        "                    U(0, 0, 0, pi / 2 * generalized_signs[k], walk_data[0]),\n",
        "                    power(k, lambda: walk_operator(walk_block, walk_data)),\n",
        "                ),\n",
        "            ),\n",
        "        ),\n",
        "    )"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "4775a70f-93f2-4fc7-a435-bafa3d484514",
      "metadata": {
        "id": "4775a70f-93f2-4fc7-a435-bafa3d484514"
      },
      "outputs": [],
      "source": [
        "@qfunc\n",
        "def main(ham_block: Output[QNum], data: Output[QNum], exp_block: Output[QNum]):\n",
        "    allocate(lcu_size_exp, exp_block)\n",
        "    allocate(lcu_size_ham, ham_block)\n",
        "    my_state_prep(data)  # Initial state prep block\n",
        "    lcu_cheb(\n",
        "        prepare_probs_exp,\n",
        "        generalized_signs,\n",
        "        lambda x, y: my_walk_operator(x, y),\n",
        "        ham_block,\n",
        "        data,\n",
        "        exp_block,\n",
        "    )\n",
        "\n",
        "\n",
        "#qmod = create_model(main, execution_preferences=Execution_Prefs)\n",
        "#write_qmod(qmod, f'hamiltonian_simulation_qubitization_{N}_oscillators', decimal_precision=12)\n",
        "#qprog = synthesize(qmod)\n",
        "#show(qprog)\n",
        "\n",
        "qmod_with_constraints = set_constraints(\n",
        "    qmod, Constraints(optimization_parameter=\"depth\", max_width=100)\n",
        ")\n",
        "\n",
        "qprog = synthesize(qmod_with_constraints)\n",
        "show(qprog)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "e8932937-2511-48b4-9262-f5f6f6f46791",
      "metadata": {
        "id": "e8932937-2511-48b4-9262-f5f6f6f46791"
      },
      "outputs": [],
      "source": [
        "results = execute(qprog).result()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "5f4f52bc-4882-429c-8a9f-523eea9cc037",
      "metadata": {
        "id": "5f4f52bc-4882-429c-8a9f-523eea9cc037"
      },
      "outputs": [],
      "source": [
        "state_result = get_projected_state_vector(\n",
        "    results, \"data\", {\"exp_block\": 0.0, \"ham_block\": 0.0}\n",
        ")"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "264eae37-5210-40bc-bcf7-51c5ab29f3d5",
      "metadata": {
        "id": "264eae37-5210-40bc-bcf7-51c5ab29f3d5"
      },
      "outputs": [],
      "source": [
        "expected_state = (\n",
        "    1 / normalization_exp * scipy.linalg.expm(1j * matrix * EVOLUTION_TIME) @ (U_matrix @ amps)\n",
        ")\n",
        "relative_phase = np.angle(expected_state[0] / state_result[0])\n",
        "state_result = state_result * np.exp(\n",
        "    1j * relative_phase\n",
        ")  # rotate according to a global phase\n",
        "print(\"expected state:\", expected_state)\n",
        "print(\"resulting state:\", state_result)\n",
        "#assert np.linalg.norm(state_result - expected_state) < EPS"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "fd4b83b7-48ee-4ae0-bb26-adf7bf057b27",
      "metadata": {
        "id": "fd4b83b7-48ee-4ae0-bb26-adf7bf057b27"
      },
      "outputs": [],
      "source": [
        "print(\n",
        "    \"overlap between expected and resulting state:\",\n",
        "    np.abs(np.vdot(state_result, expected_state))\n",
        "    * normalization_exp\n",
        "    / np.linalg.norm(state_result),\n",
        ")"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "460ec011-b4ae-4152-b8b7-7ae260088f7b",
      "metadata": {
        "id": "460ec011-b4ae-4152-b8b7-7ae260088f7b"
      },
      "outputs": [],
      "source": []
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "2f5099b3-181c-4b76-bd51-11e682f861d2",
      "metadata": {
        "id": "2f5099b3-181c-4b76-bd51-11e682f861d2"
      },
      "outputs": [],
      "source": []
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "b50d1397-6b01-4a46-b03a-581b7daceb0c",
      "metadata": {
        "id": "b50d1397-6b01-4a46-b03a-581b7daceb0c"
      },
      "outputs": [],
      "source": []
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "id": "846b05d1-be6b-4197-9dd9-cba85e553c27",
      "metadata": {
        "id": "846b05d1-be6b-4197-9dd9-cba85e553c27"
      },
      "outputs": [],
      "source": []
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3 (ipykernel)",
      "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.11.7"
    },
    "colab": {
      "provenance": [],
      "include_colab_link": true
    }
  },
  "nbformat": 4,
  "nbformat_minor": 5
}