{
  "cells": [
    {
      "cell_type": "markdown",
      "id": "frontmatter",
      "metadata": {},
      "source": [
        "---\n",
        "title: \"Quickstart\"\n",
        "description: \"Quickstart for the latest version of Sample-based quantum diagonalization (SQD)\"\n",
        "---\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "intro",
      "metadata": {},
      "source": [
        "# Quickstart\n",
        "\n",
        "{/* cspell:ignore mcscf, ncas, nelecas, hcore, binom */}\n",
        "\n",
        "This guide demonstrates a minimal working example of the `qiskit-addon-sqd` package. We use sample-based quantum diagonalization (SQD) to refine noisy samples taken from a quantum circuit that approximates the ground state of a chemistry Hamiltonian. Ground state energies are approximated by diagonalizing the Hamiltonian in subspaces formed by the refined samples.\n",
        "\n",
        "For more explanatory end-to-end workflows that run on quantum hardware, check out the tutorials:\n",
        "\n",
        "* [Sample-based quantum diagonalization](/docs/tutorials/sample-based-quantum-diagonalization): chemistry Hamiltonian workflow\n",
        "* [Sample-based Krylov quantum diagonalization](/docs/tutorials/sample-based-krylov-quantum-diagonalization): fermionic lattice model workflow\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "prepare-header",
      "metadata": {},
      "source": [
        "## 1. Prepare the inputs for SQD\n",
        "\n",
        "First, we define the chemistry Hamiltonian and quantum circuit that approximates its ground state. The system studied in this example is N₂ at equilibrium geometry in STO-3G basis with two core orbitals frozen. This system has `8` spatial orbitals and `10` electrons, which means the quantum circuit is defined on `16` qubits, where qubits `0-7` represent spin-up orbitals and qubits `8-15` represent spin-down orbitals. Since this is a closed-shell system with `10` electrons, we know each physically-valid sample should have hamming weight `10`, where the hamming weight of both the left and right halves of the bitstring is `5`.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 1,
      "id": "load-molecule",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Spatial orbitals: 8\n",
            "Qubits: 16\n",
            "Electrons (alpha, beta): (5, 5)\n"
          ]
        }
      ],
      "source": [
        "from pyscf import ao2mo, gto, mcscf, scf\n",
        "\n",
        "# Construct the Hamiltonian, pull out integrals and other information about the system\n",
        "mol = gto.M(\n",
        "    atom=\"N 0 0 0; N 0 0 1.09768\", basis=\"sto-3g\", spin=0, charge=0, verbose=0\n",
        ")\n",
        "mf = scf.RHF(mol).run()\n",
        "cas = mcscf.CASCI(mf, ncas=8, nelecas=10)\n",
        "hcore, nuclear_repulsion_energy = cas.get_h1eff()\n",
        "eri = cas.get_h2eff()\n",
        "num_orbitals = cas.ncas\n",
        "eri = ao2mo.restore(1, eri, num_orbitals)\n",
        "n_alpha, n_beta = cas.nelecas\n",
        "nelec = (n_alpha, n_beta)\n",
        "\n",
        "print(\n",
        "    f\"Spatial orbitals: {num_orbitals}\\nQubits: {num_orbitals * 2}\\nElectrons (alpha, beta): {nelec}\"\n",
        ")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "samples-explanation",
      "metadata": {},
      "source": [
        "## 2. Generate noisy samples\n",
        "\n",
        "In a real workflow you would design an ansatz circuit, optimize it for hardware, and sample it on the QPU. Here we skip that and draw `10,000` uniformly random bitstrings instead. While our samples contain no information about the actual ground state, we will see that the SQD protocol can still recover the true ground state by iteratively refining these samples, diagonalizing, and using the energies at each iteration to inform the next round of bitstring refinement. We can eventually recover the exact ground state in this manner because the subspace containing the ground state is small enough to calculate exactly. Specifically, since we know each half of each bitstring (`8` bits) should have hamming weight `5`, we know the full Hilbert subspace of the system is size $\\binom{8}{5}^2=3136$.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 2,
      "id": "generate-samples",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Generated 10000 uniformly-random, 16-qubit samples.\n"
          ]
        }
      ],
      "source": [
        "import numpy as np\n",
        "from qiskit_addon_sqd.counts import generate_bit_array_uniform\n",
        "\n",
        "# Generate some uniformly-random bitstrings to simulate QPU samples\n",
        "rng = np.random.default_rng(24)\n",
        "bit_array = generate_bit_array_uniform(\n",
        "    10_000, num_orbitals * 2, rand_seed=rng\n",
        ")\n",
        "print(\n",
        "    f\"Generated {bit_array.num_shots} uniformly-random, {bit_array.num_bits}-qubit samples.\"\n",
        ")"
      ]
    },
    {
      "cell_type": "markdown",
      "id": "postprocess-header",
      "metadata": {},
      "source": [
        "## 3. Post-process with SQD\n",
        "\n",
        "The objective of the SQD algorithm is to refine noisy samples taken from a QPU and ultimately diagonalize the Hamiltonian in a subspace spanned by those refined samples. This is done iteratively as follows:\n",
        "\n",
        "* Refine samples by postselecting only those which have correct particle number for both the spin-up and spin-down subsystems\n",
        "* Diagonalize in small subspaces spanned by the refined samples and collect information about the subspaces that yield the the best (lowest) energies.\n",
        "* For each iteration:\n",
        "  * Refine the noisy samples by probabilistically flipping bits in samples with incorrect particle number. The probabilities are calculated using the subspace information from the most recent diagonalization.\n",
        "  * Diagonalize in small subspaces spanned by the refined samples\n",
        "  * Collect information about the subspaces yielding the best (lowest) energies\n",
        "\n",
        "We can see from the output that the algorithm converges to the exact ground state energy after `6` iterations and that the diagonalization is performed over the full Hilbert subspace of dimension $\\binom{8}{5}^2=3136$. Of course, diagonalizing over the full subspace is usually not feasible, and therefore, systems with ground states having sparse support are the most likely to benefit from subspace diagonalization techniques such as SQD.\n",
        "\n",
        "For more information on the diagonalize\\_fermionic\\_hamiltonian function, visit the [API docs](/docs/api/qiskit-addon-sqd/fermion#diagonalize_fermionic_hamiltonian).\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 3,
      "id": "run-sqd",
      "metadata": {},
      "outputs": [
        {
          "name": "stdout",
          "output_type": "stream",
          "text": [
            "Exact (FCI) reference energy: -107.652521 Ha\n",
            "\n",
            "Iteration 1\n",
            "  Subsample 0\n",
            "    Energy: -106.908481\n",
            "    Subspace dimension: 784\n",
            "    Error vs exact: 0.744040 Ha\n",
            "Iteration 2\n",
            "  Subsample 0\n",
            "    Energy: -106.934921\n",
            "    Subspace dimension: 1764\n",
            "    Error vs exact: 0.717600 Ha\n",
            "Iteration 3\n",
            "  Subsample 0\n",
            "    Energy: -106.953430\n",
            "    Subspace dimension: 2209\n",
            "    Error vs exact: 0.699091 Ha\n",
            "Iteration 4\n",
            "  Subsample 0\n",
            "    Energy: -107.651946\n",
            "    Subspace dimension: 2916\n",
            "    Error vs exact: 0.000575 Ha\n",
            "Iteration 5\n",
            "  Subsample 0\n",
            "    Energy: -107.652003\n",
            "    Subspace dimension: 3025\n",
            "    Error vs exact: 0.000518 Ha\n",
            "Iteration 6\n",
            "  Subsample 0\n",
            "    Energy: -107.652521\n",
            "    Subspace dimension: 3136\n",
            "    Error vs exact: 0.000000 Ha\n",
            "Iteration 7\n",
            "  Subsample 0\n",
            "    Energy: -107.652521\n",
            "    Subspace dimension: 3136\n",
            "    Error vs exact: 0.000000 Ha\n"
          ]
        }
      ],
      "source": [
        "from functools import partial\n",
        "\n",
        "from pyscf import fci\n",
        "from qiskit_addon_sqd.fermion import (\n",
        "    SCIResult,\n",
        "    diagonalize_fermionic_hamiltonian,\n",
        "    solve_sci_batch,\n",
        ")\n",
        "\n",
        "exact_energy, _ = fci.direct_spin1.kernel(\n",
        "    hcore, eri, int(num_orbitals), nelec\n",
        ")\n",
        "exact_energy += nuclear_repulsion_energy\n",
        "print(f\"Exact (FCI) reference energy: {exact_energy:.6f} Ha\\n\")\n",
        "\n",
        "sci_solver = partial(solve_sci_batch, spin_sq=0.0, max_cycle=200)\n",
        "\n",
        "result_history = []\n",
        "\n",
        "\n",
        "def callback(results: list[SCIResult]):\n",
        "    result_history.append(results)\n",
        "    iteration = len(result_history)\n",
        "    print(f\"Iteration {iteration}\")\n",
        "    for i, result in enumerate(results):\n",
        "        error = result.energy + nuclear_repulsion_energy - exact_energy\n",
        "        print(f\"  Subsample {i}\")\n",
        "        print(f\"    Energy: {result.energy + nuclear_repulsion_energy:.6f}\")\n",
        "        print(\n",
        "            f\"    Subspace dimension: {np.prod(result.sci_state.amplitudes.shape)}\"\n",
        "        )\n",
        "        print(f\"    Error vs exact: {error:.6f} Ha\")\n",
        "\n",
        "\n",
        "samples_per_batch = 20\n",
        "result = diagonalize_fermionic_hamiltonian(\n",
        "    hcore,\n",
        "    eri,\n",
        "    bit_array,\n",
        "    samples_per_batch=samples_per_batch,\n",
        "    norb=num_orbitals,\n",
        "    nelec=nelec,\n",
        "    occupancies_tol=1e-7,\n",
        "    max_iterations=30,\n",
        "    sci_solver=sci_solver,\n",
        "    symmetrize_spin=True,\n",
        "    callback=callback,\n",
        "    seed=np.random.default_rng(32),\n",
        ")"
      ]
    },
    {
      "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
}