{
  "cells": [
    {
      "cell_type": "markdown",
      "id": "frontmatter",
      "metadata": {},
      "source": [
        "---\n",
        "title: \"Quickstart\"\n",
        "description: \"Quickstart for the latest version of Pauli propagation\"\n",
        "---\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "a9320fa8-31c5-4248-96b0-3549a13dda6f",
      "metadata": {},
      "source": [
        "---\n",
        "title: Quickstart\n",
        "description: A quickstart guide for the pauli-prop Qiskit addons package\n",
        "---\n",
        "\n",
        "# Quickstart\n",
        "\n",
        "In this guide we use the `pauli-prop` package to classically simulate the time dynamics of a 10-qubit kicked Ising model on a 1D spin chain.\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "3b5bf7dc-1cd8-41d3-b8ce-cc75a66fed8d",
      "metadata": {},
      "source": [
        "## Prepare the inputs for Pauli propagation\n",
        "\n",
        "The Hamiltonian considered is:\n",
        "\n",
        "$H = -J\\sum\\limits_{\\langle i,j \\rangle} Z_iZ_j + h\\sum\\limits_iX_i$\n",
        "\n",
        "where $J>0$ describes the coupling of nearest-neighbor spins, $i<j$, and $h$ is the global transverse field. A first-order Trotter decomposition of the time-evolved operator will be implemented as a quantum circuit, $U$, over $20$ Trotter steps. The coupling constant, $J$, will be fixed at $J=-\\frac{\\pi}{2}$, and $h$ will be fixed at $\\frac{\\pi}{6}$. The $ZZ$ interactions will be implemented using Clifford gates ($CX$, $Sdg$, $\\sqrt{Y}$).\n",
        "\n",
        "We implement the Trotterized time evolution as a quantum circuit and use $\\frac{\\pi}{6}$ for the non-Clifford rotations about the x-axis. The further these angles are from Clifford angles (for example, $\\theta=n\\frac{\\pi}{2}, n \\in \\mathbb{Z}$), the more difficult the system will be to simulate for Pauli propagation methods.\n",
        "\n",
        "For the choice of observable, we consider the average single-site magnetization, $\\frac{1}{N} \\sum_{i=1}^{N} \\langle z_i \\rangle$, where $N$ is the number of spins.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "id": "15d5eccf-eef0-4435-b67a-98a8038c400e",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "<Image src=\"/docs/images/addons/pauli-prop/guides/quickstart/extracted-outputs/15d5eccf-eef0-4435-b67a-98a8038c400e-0.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "execution_count": 1,
          "metadata": {},
          "output_type": "execute_result"
        }
      ],
      "source": [
        "import numpy as np\n",
        "from qiskit import QuantumCircuit\n",
        "from qiskit.quantum_info import SparsePauliOp\n",
        "from qiskit.transpiler import CouplingMap\n",
        "\n",
        "num_qubits = 10\n",
        "coupling_map = CouplingMap.from_line(num_qubits, bidirectional=False)\n",
        "\n",
        "# Num Trotter steps\n",
        "num_steps = 20\n",
        "theta_rx = np.pi / 6\n",
        "\n",
        "# Average single-site magnetization\n",
        "observable = (\n",
        "    SparsePauliOp(\n",
        "        [\n",
        "            \"I\" * iq + \"Z\" + \"I\" * (num_qubits - iq - 1)\n",
        "            for iq in range(num_qubits)\n",
        "        ]\n",
        "    )\n",
        "    / num_qubits\n",
        ")\n",
        "\n",
        "# Create the Trotter circuit\n",
        "num_qubits = 10\n",
        "num_steps = 20\n",
        "theta_rx = np.pi / 6\n",
        "circuit = QuantumCircuit(num_qubits)\n",
        "edges = CouplingMap.from_line(num_qubits, bidirectional=False).get_edges()\n",
        "for _ in range(num_steps):\n",
        "    circuit.rx(theta_rx, [i for i in range(num_qubits)])\n",
        "    for edge in edges:\n",
        "        circuit.sdg(edge)\n",
        "        circuit.ry(np.pi / 2, edge[1])\n",
        "        circuit.cx(edge[0], edge[1])\n",
        "        circuit.ry(-np.pi / 2, edge[1])\n",
        "circuit.draw(\"mpl\", fold=-1)"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "c05fa0f8-605e-4172-8880-0de89360f27d",
      "metadata": {},
      "source": [
        "## Simulate the time evolution of the system with Pauli propagation\n",
        "\n",
        "Once we have our circuit, $U$, and observable, $O$, we can easily simulate the system in a few steps:\n",
        "\n",
        "* Separate $U$, into its Clifford, $C$, and non-Clifford, $P$, parts such that $U=PC$ using `evolve_through_cliffords`\n",
        "* Propagate $O$ through $P$, resulting in a new operator, $O^\\prime$, using `pauli_prop.propagate_through_circuit`\n",
        "* Evolve $O^\\prime$ through the Clifford part of the circuit using the Clifford evolution support built into Qiskit\n",
        "* Approximate the expectation value as $\\langle0|O^\\prime|0\\rangle \\approx \\langle0|U^\\dagger OU|0\\rangle$ by summing the coefficients in $O^\\prime$ associated with fully-diagonal Pauli terms (Pauli terms containing either `I` or `Z` on all qubits). Remember, this is an approximation because we truncated terms from $O^\\prime$ as we propagated it through the non-Clifford part of the circuit.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 2,
      "id": "c8a27469-ab72-4d50-9544-d446af847da5",
      "metadata": {},
      "outputs": [],
      "source": [
        "import time\n",
        "\n",
        "from pauli_prop import evolve_through_cliffords, propagate_through_circuit\n",
        "\n",
        "cliff, non_cliff = evolve_through_cliffords(circuit)\n",
        "\n",
        "max_terms_list = [10**i for i in range(8)]\n",
        "approx_evs = []\n",
        "durations = []\n",
        "for max_terms in max_terms_list:\n",
        "    st = time.perf_counter()\n",
        "    evolved_obs = propagate_through_circuit(\n",
        "        observable, non_cliff, max_terms=max_terms, atol=1e-12, frame=\"h\"\n",
        "    )[0]\n",
        "    evolved_obs.paulis = evolved_obs.paulis.evolve(cliff, frame=\"h\")\n",
        "    durations.append(time.perf_counter() - st)\n",
        "    approx_evs.append(\n",
        "        float(evolved_obs.coeffs[~evolved_obs.paulis.x.any(axis=1)].sum())\n",
        "    )"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "7fb3e72f-20fb-4fd0-a301-2fcc6c0ddad5",
      "metadata": {},
      "source": [
        "As we run larger calculations, the expectation value approximations become more accurate. In this example, we saturate the full Pauli space at around $4^{10}\\approx10^6$, which is reflected in the curve flattening out between the final two points.\n",
        "\n",
        "While the plot below shows a monotonic convergence, Pauli propagation simulations do not generally converge monotonically. It is not unusual to see \"bumpy\" behavior in these types of plots.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 3,
      "id": "f12410aa-ef3b-4ee6-ab35-f7e9001a92b8",
      "metadata": {},
      "outputs": [
        {
          "data": {
            "text/plain": [
              "Text(0.5, 1.0, 'Simulating 20-step 1D Ising Model')"
            ]
          },
          "execution_count": 3,
          "metadata": {},
          "output_type": "execute_result"
        },
        {
          "data": {
            "text/plain": [
              "<Image src=\"/docs/images/addons/pauli-prop/guides/quickstart/extracted-outputs/f12410aa-ef3b-4ee6-ab35-f7e9001a92b8-1.avif\" alt=\"Output of the previous code cell\" />"
            ]
          },
          "metadata": {},
          "output_type": "display_data"
        }
      ],
      "source": [
        "import matplotlib.pyplot as plt\n",
        "from qiskit_aer import AerSimulator\n",
        "\n",
        "sim_circ = circuit.copy()\n",
        "sim_circ.save_statevector()\n",
        "backend = AerSimulator(method=\"statevector\")\n",
        "psi = backend.run(sim_circ).result().data()[\"statevector\"]\n",
        "exact_ev = psi.expectation_value(observable)\n",
        "\n",
        "ax1 = plt.gca()\n",
        "ax1.plot(max_terms_list, approx_evs, marker=\"o\", label=\"Approximate\")\n",
        "ax1.axhline(exact_ev, linestyle=\"--\", color=\"green\", label=\"Exact\")\n",
        "ax1.set_xscale(\"log\")\n",
        "ax1.set_xlabel(\"# terms kept\")\n",
        "ax1.set_ylabel(r\"$\\frac{1}{N} \\sum_{i=1}^{N} \\langle z_i \\rangle$\")\n",
        "\n",
        "ax2 = ax1.twinx()\n",
        "ax2.plot(\n",
        "    max_terms_list, durations, marker=\".\", label=\"Runtime\", color=\"orange\"\n",
        ")\n",
        "ax2.set_ylabel(\"Runtime (s)\", color=\"orange\")\n",
        "ax2.set_yscale(\"log\")\n",
        "\n",
        "handles1, labels1 = ax1.get_legend_handles_labels()\n",
        "handles2, labels2 = ax2.get_legend_handles_labels()\n",
        "ax1.legend(handles1 + handles2, labels1 + labels2, loc=\"lower right\")\n",
        "\n",
        "plt.title(f\"Simulating {num_steps}-step 1D Ising Model\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "id": "a1b8767d",
      "source": "© IBM Corp., 2017-2026"
    }
  ],
  "metadata": {
    "kernelspec": {
      "display_name": "Python 3",
      "language": "python",
      "name": "python3"
    },
    "language_info": {
      "codemirror_mode": {
        "name": "ipython",
        "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 5
}