docs: demonstrate the packaged phase estimation API

docs/sphinx/examples/python/08_quantum_phase_estimation.py

@@ -9,6 +9,21 @@
 # [Begin Documentation]
 """Example 8 — Quantum phase estimation with three energy oracles.
 
+The physical question is how to estimate a molecule's energy using a quantum
+circuit. Phase estimation extracts an energy-dependent phase into a counting
+register whose measurement produces a distribution of estimates. This example
+uses the same small hydrogen molecule and initial state throughout, so readers
+can follow how exact time evolution, approximate time evolution, and a
+qubitization walk each connect their phases to the same energy problem.
+
+You should learn how the initial state's overlap with energy eigenstates
+determines which energies are observed, and how counting-register resolution
+limits their precision. The explicit circuits and the packaged PhaseEstimation
+API show the roles of controlled operations, measured bits, and energy
+decoding. Comparing their distributions with exact references also explains
+why walk phases require a cosine-based decoder and mirror folding, and how to
+judge sampling variation separately from a systematic error.
+
 QPE turns the eigenphases of a unitary into measured bit strings. This
 example estimates ground-state energies of the H2/STO-3G Hamiltonian with
 one QPE shell around three different oracle constructions:
@@ -29,23 +44,33 @@
 Hartree-Fock has dominant ground-state overlap. Each modal estimate is
 checked to land within one phase cell of the exact ground energy.
 
+The final section executes PhaseEstimation's packaged textbook walk ladder
+on the same input. Its public decoders, exact phase marginal, sampled
+statistical check, and cosine phase cells are compared with the
+centered walk above. This is another implementation of method 3, not another
+physical oracle.
+
 The walk QPE core follows the centered exponent schedule of Berry et al.,
 "Rapid initial-state preparation for the quantum simulation of strongly
 correlated molecules", PRX Quantum 6, 020327 (2025).
 
 Prerequisite: PySCF (`pip install pyscf`). Run: python3 08_quantum_phase_estimation.py
+Add --verbose to show individual validation checks; failures are always reported.
 """
 
 from __future__ import annotations
 
+import argparse
 import math
 import os
+import time
 
 import cudaq
 import numpy as np
 import scipy
 
-from cudaq_algorithms import PauliLCU, chemistry, stateprep, trotter
+from cudaq_algorithms import PauliLCU, PhaseEstimation, chemistry, stateprep, trotter
+from cudaq_algorithms import sim_utils as sim
 
 # ----------------------------------------------------------------------
 # Classical decoding helpers: counts -> probabilities -> modal bin -> energy
@@ -105,7 +130,136 @@ def decode_walk_qpe(probabilities: np.ndarray,
     return mode, -alpha * math.cos(phase)
 
 
+def check(label, error, tolerance=1e-10, *, verbose=False):
+    passed = np.isfinite(error) and error <= tolerance
+    if verbose or not passed:
+        print(f"CHECK {label}: {error:.3e} <= {tolerance:.3e} -> "
+              f"{'PASS' if passed else 'FAIL'}")
+    if not passed:
+        raise AssertionError(label)
+
+
+def packaged_walk_comparison(encoding,
+                             prepare_hartree_fock,
+                             centered_kernel,
+                             centered_sampled,
+                             counting_qubits,
+                             shots,
+                             exact_ground_energy,
+                             *,
+                             verbose=False):
+    print("\nPackaged PhaseEstimation: the same walk, a textbook power ladder")
+    qpe = PhaseEstimation(encoding)
+    measured = qpe.kernel(counting_qubits,
+                          state_prep=prepare_hartree_fock,
+                          measure=True)
+    cudaq.set_random_seed(109)
+    counts = cudaq.sample(measured, shots_count=shots)
+    ranked = qpe.energies_from_counts(counts, counting_qubits, top=3)
+    sampled = qpe.phase_probabilities_from_counts(counts, counting_qubits)
+    # The public decoder also accepts full-register counts, whose layout is
+    # system | phase | ancilla. No manual string slicing is needed here.
+    check("only phase bits measured",
+          max(abs(len(str(bits)) - counting_qubits) for bits in counts),
+          0,
+          verbose=verbose)
+    exact = sim.qpe.phase_probabilities(qpe,
+                                        counting_qubits,
+                                        state_prep=prepare_hartree_fock)
+    bins = 1 << counting_qubits
+    check("exact phase marginal normalized",
+          abs(float(exact.sum()) - 1),
+          verbose=verbose)
+    # A priori union/Hoeffding bound over phase-bin subsets, delta=1e-6.
+    # Compare unfolded distributions before identifying conjugate energies.
+    delta = 1e-6
+    tv_bound = min(
+        1., np.sqrt((bins * np.log(2) + np.log(2 / delta)) / (2 * shots)))
+    tv = lambda a, b: 0.5 * float(np.sum(np.abs(a - b)))
+    check("sampled versus exact unfolded phase TV",
+          tv(sampled, exact),
+          tv_bound,
+          verbose=verbose)
+
+    # Same physical oracle, different schedule and allocation order. The
+    # centered schedule applies W^(n - M/2) instead of W^n; the extra
+    # W^(-M/2) is a phase on each walk eigencomponent, so the phase-register
+    # marginal agrees bin by bin before any folding. Disable measurement in
+    # the centered shell to obtain its exact phase marginal.
+    centered_state = np.asarray(
+        cudaq.get_state(centered_kernel, counting_qubits, encoding.num_system,
+                        encoding.num_ancilla, False))
+    centered_exact = np.sum(np.abs(
+        centered_state.reshape(1 << encoding.num_system,
+                               1 << encoding.num_ancilla, bins))**2,
+                            axis=(0, 1))
+    check("centered versus textbook exact phase marginals",
+          tv(centered_exact, exact),
+          verbose=verbose)
+    check("centered sampled versus its exact unfolded phase TV",
+          tv(centered_sampled, centered_exact),
+          tv_bound,
+          verbose=verbose)
+    check("centered versus textbook sampled phase TV",
+          tv(centered_sampled, sampled),
+          2 * tv_bound,
+          verbose=verbose)
+    decoded = np.array(
+        [qpe.phase_to_energy(j, counting_qubits) for j in range(bins)])
+    expected = -encoding.alpha * np.cos(2 * np.pi * np.arange(bins) / bins)
+    check("public cosine decoder on every phase bin",
+          float(np.max(np.abs(decoded - expected))),
+          verbose=verbose)
+
+    mode, _ = decode_walk_qpe(exact, encoding.alpha)
+    sampled_mode, sampled_energy = decode_walk_qpe(sampled, encoding.alpha)
+    centered_mode, _ = decode_walk_qpe(centered_sampled, encoding.alpha)
+    check("sampled dominant energy matches public ranking",
+          abs(sampled_energy - ranked[0][0]),
+          verbose=verbose)
+    check("centered and textbook dominant energy cell",
+          abs(centered_mode - mode),
+          0,
+          verbose=verbose)
+    check("sampled and exact dominant energy cell",
+          abs(sampled_mode - mode),
+          0,
+          verbose=verbose)
+    lower = -encoding.alpha * math.cos(
+        2 * math.pi * max(0., mode - 0.5) / bins)
+    upper = -encoding.alpha * math.cos(
+        2 * math.pi * min(bins / 2, mode + 0.5) / bins)
+    check("exact ground energy lies in cosine phase cell",
+          max(lower - exact_ground_energy, exact_ground_energy - upper, 0.),
+          1e-9,
+          verbose=verbose)
+    print(
+        f"m={counting_qubits}, shots={shots}, alpha={encoding.alpha:.9f}, "
+        f"total qubits={encoding.num_system+counting_qubits+encoding.num_ancilla}, "
+        f"textbook controlled-walk steps={bins-1}")
+    print(f"sampling seed=109; unfolded TV={tv(sampled, exact):.6f}, "
+          f"a priori bound={tv_bound:.6f}, delta={delta}")
+    for value, probability in ranked:
+        print(f"  folded E={value:+.12f} Ha, probability={probability:.6f}")
+    print(f"ground-energy cell: [{lower:+.9f}, {upper:+.9f}] Ha")
+    print("The encoded Hamiltonian already includes core/nuclear energy once.")
+    print("PhaseEstimation builds circuits and decodes; the caller executes.")
+    print(
+        "QPE weights reflect input overlaps; it does not select the ground state."
+    )
+    print(
+        "Long simulated ladders (from m=13) need sim.run_on_large_stack, "
+        "including when cudaq.sample is the executor; they are not run here.")
+
+
 def main() -> int:
+    parser = argparse.ArgumentParser(description=__doc__.splitlines()[0])
+    parser.add_argument(
+        "--verbose",
+        action="store_true",
+        help="show each validation check (failures are always reported)")
+    verbose = parser.parse_args().verbose
+    started = time.perf_counter()
     cudaq.set_target(os.environ.get("CUDAQ_DEFAULT_SIMULATOR", "qpp-cpu"))
 
     try:
@@ -292,7 +446,8 @@ def qpe_trotter(n_counting_q: int, n_system_q: int,
     controlled_walk_dagger = encoding.controlled_adjoint_walk_step_kernel()
 
     @cudaq.kernel
-    def qpe_walk(n_counting_q: int, n_system_q: int, n_signal_q: int):
+    def qpe_walk(n_counting_q: int, n_system_q: int, n_signal_q: int,
+                 measure_phase: bool):
         counting_and_signal = cudaq.qvector(n_counting_q + n_signal_q)
         counting = counting_and_signal.front(n_counting_q)
         signal = counting_and_signal.back(n_signal_q)
@@ -327,7 +482,8 @@ def qpe_walk(n_counting_q: int, n_system_q: int, n_signal_q: int):
                 swap(counting[bit], counting[n_counting_q - 1])
 
         iqft(counting)
-        mz(counting)
+        if measure_phase:
+            mz(counting)
 
     # 6. Execute all three circuits: the same phase register, the same
     #    Hamiltonian, the same Hartree-Fock state, and 20,000 shots each.
@@ -360,6 +516,7 @@ def qpe_walk(n_counting_q: int, n_system_q: int, n_signal_q: int):
                                counting_qubits,
                                system_qubit_count,
                                encoding.num_ancilla,
+                               True,
                                shots_count=shots)
     walk_sampled = sampled_probabilities(walk_counts, counting_qubits, shots)
 
@@ -401,6 +558,15 @@ def qpe_walk(n_counting_q: int, n_system_q: int, n_signal_q: int):
     )
     assert walk_lower - 1e-9 <= exact_ground_energy <= walk_upper + 1e-9
 
+    packaged_walk_comparison(encoding,
+                             prepare_hartree_fock,
+                             qpe_walk,
+                             walk_sampled,
+                             counting_qubits,
+                             shots,
+                             exact_ground_energy,
+                             verbose=verbose)
+    print(f"Total example elapsed: {time.perf_counter()-started:.2f}s")
     print("\nOK — all three QPE oracles land within one phase cell of the "
           "exact ground energy.")
     return 0

docs/sphinx/examples_rst/getting_started.rst

@@ -104,6 +104,24 @@ Requires ``qsppack`` for the phase factors, like example 2.
 Quantum phase estimation with three energy oracles
 --------------------------------------------------
 
+The physical question is how to estimate a molecule's energy using a quantum
+circuit. Phase estimation extracts an energy-dependent phase into a counting
+register whose measurement produces a distribution of estimates. This example
+uses the same small hydrogen molecule and initial state throughout, so readers
+can follow how exact time evolution, approximate time evolution, and a
+qubitization walk each connect their phases to the same energy problem.
+
+You should learn how the initial state's overlap with energy eigenstates
+determines which energies are observed, and how counting-register resolution
+limits their precision. The explicit circuits and the packaged PhaseEstimation
+API show the roles of controlled operations, measured bits, and energy
+decoding. Comparing their distributions with exact references also explains
+why walk phases require a cosine-based decoder and mirror folding, and how to
+judge sampling variation separately from a systematic error.
+
+Individual validation checks are shown with ``--verbose``. All checks still
+run by default, and failures are always reported.
+
 Quantum phase estimation reads the eigenphases of a unitary into a counting
 register: for :math:`U|E_j\rangle = e^{2\pi i \phi_j}|E_j\rangle`, an
 :math:`m`-qubit register and the controlled powers :math:`U^{2^k}` build a
@@ -133,8 +151,18 @@ three different constructions of the unitary whose phase reveals energy:
      - :math:`E = -\alpha\cos\theta`
 
 The walk case is what `cudaq_algorithms.PhaseEstimation` packages (see
-:doc:`../guide/qubitization_qsvt`); the example keeps its own shell so the
-three oracles share one circuit and one readout.
+:doc:`../guide/qubitization_qsvt`). After the three-oracle comparison, the
+example executes that packaged interface with the same Hamiltonian and
+Hartree-Fock preparation: build a measured kernel, call ``cudaq.sample``,
+and decode with ``energies_from_counts``. It is another implementation
+of the walk approach, rather than a fourth physical oracle.
+
+The packaged section checks that the textbook and centered exact phase
+marginals agree bin by bin, compares seeded sampled distributions with
+their exact marginals under a declared statistical bound, and verifies
+the nonlinear cosine phase-cell edges. Only phase bits
+are measured; public decoders handle register interpretation and mirror
+folding. The encoded identity/core energy is included once.
 
 The product-formula oracle shows a subtlety of composing library evolution
 under control: `Trotter` evolves the nonidentity part of the Hamiltonian, so
@@ -155,7 +183,8 @@ cheaper classically. The walk is here because fault-tolerant algorithms
 commonly assume structured PREPARE/SELECT access, and this shell is exactly
 how such an oracle drops into QPE. Six counting qubits keep the script fast
 on the CPU simulator; every modal estimate is asserted to land within one
-phase cell (about 34 mHa here) of the exact ground energy. Resolution scales
+phase cell of the exact ground energy. The time-oracle half-cell is about
+34 mHa here; walk cells use their own asymmetric cosine edges. Resolution scales
 with the register: at 12 counting qubits the same circuits reach roughly
 0.1--0.5 mHa.
 

docs/sphinx/guide/qubitization_qsvt.rst

@@ -183,6 +183,31 @@ puts a hand-written QPE shell around three oracles — exact time evolution,
 a product formula and this walk — to compare their decoders; `PhaseEstimation`
 packages the walk case.
 
+The worked example also executes the packaged API on the same H2 Hamiltonian
+and Hartree-Fock preparation, alongside the three-oracle comparison.
+This is another implementation of the walk approach, not a fourth physical
+oracle. Core/nuclear energy is already encoded once.
+
+It checks that the exact centered and textbook phase marginals agree bin by
+bin (the centered schedule's extra :math:`W^{-M/2}` is a phase on each walk
+eigencomponent, so the readout distribution is unchanged), then compares
+seeded 20,000-shot histograms to the exact *unfolded* marginals using an
+a priori total-variation bound. With M=2^m bins and failure budget delta,
+the bound is
+
+.. math::
+
+   t=\min\left(1,\sqrt{\frac{M\log 2+\log(2/\delta)}{2S}}\right).
+
+At m=6, S=20,000 and delta=1e-6, this is about 0.0384. The sampling check
+uses unfolded distributions: folding mirror bins coarsens a distribution and
+can only lower total variation, so a folded check is weaker. Energy
+acceptance uses the cosine of
+the neighboring half-bin edges, clipped to the folded endpoints, rather
+than a uniform energy tolerance or fitted offset. QPE retains the input's
+spectral overlap weights; a high-overlap Hartree-Fock state explains the
+dominant ground-energy cell.
+
 State preparation injection
 ---------------------------