Pre-Trajectory Sampling with Batch Execution

.github/pre-commit/spelling_allowlist.txt

@@ -281,6 +281,7 @@ macOS
 makefiles
 merchantability
 mps
+multinomial
 multithreaded
 mutex
 namespace
@@ -319,6 +320,7 @@ qumode
 qumodes
 reStructuredText
 reconfigurable
+reproducibility
 runtime
 runtimes
 rvalue

.github/pre-commit/spelling_allowlist_cxx.txt

@@ -1,4 +1,5 @@
 clang
 gcc
 nullary
+pre
 typedef

.github/workflows/config/gitlab_commits.txt

@@ -1,2 +1,2 @@
 nvidia-mgpu-repo: cuda-quantum/cuquantum-mgpu.git
-nvidia-mgpu-commit: e41bb18a5e5ff5be44090dfc4cb6c59fb511e9c7
+nvidia-mgpu-commit: 69fa05df00a069c3bc5c040294678014a396fa68

docs/sphinx/api/languages/cpp_api.rst

@@ -120,6 +120,10 @@ Common
 .. doxygenclass:: cudaq::sample_result
     :members:
 
+.. _classcudaq_1_1ptsbe_1_1sample__result:
+
+**cudaq::`ptsbe`::sample_result** — Result type returned by ``ptsbe::sample()``, extending `cudaq::sample_result` with optional execution data (trace and per-trajectory info). See ``PTSBESampleResult.h`` and ``PTSBEExecutionData.h``.
+
 .. doxygenstruct:: cudaq::sample_options
     :members:
 

docs/sphinx/examples/python/ptsbe_end_to_end_workflow.ipynb

@@ -0,0 +1,217 @@
+{
+  "cells": [
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "## PTSBE end-to-end workflow\n",
+        "\n",
+        "**PTSBE** (Pre-Trajectory Sampling with Batch Execution) is a method for sampling from **noisy** quantum circuits efficiently. Instead of simulating the full density matrix and sampling once per shot, PTSBE:\n",
+        "\n",
+        "1. **Traces** the kernel to get the gate sequence and qubit layout.\n",
+        "2. **Extracts noise sites** by matching the noise model to the trace (each noisy gate becomes a *noise site* with a set of Kraus outcomes, e.g. I, X, Y, Z for depolarization).\n",
+        "3. **Generates trajectories** — each trajectory is one possible *realization* of noise (one Kraus outcome per site). A **sampling strategy** decides which trajectories to use (e.g. **Exhaustive**: all combinations, or **Probabilistic**: sample by probability).\n",
+        "4. **Allocates shots** across trajectories (e.g. **proportional** to trajectory probability, or **uniform**).\n",
+        "5. **Runs batches** — for each trajectory, the circuit is run as a **noiseless** circuit with that trajectory's outcomes applied; results are collected.\n",
+        "6. **Aggregates** all per-trajectory counts into a single `SampleResult`.\n",
+        "\n",
+        "You get the same statistics as density-matrix sampling (in the limit of many shots), but with the ability to batch many shots per trajectory and to control cost via the number of trajectories. This notebook runs the full workflow with a single API call: `cudaq.ptsbe.sample()`."
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "### 1. Set up the environment\n",
+        "\n",
+        "Use the density-matrix simulator target (required for PTSBE). Set a random seed for reproducibility."
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {},
+      "outputs": [],
+      "source": [
+        "import cudaq\n",
+        "\n",
+        "cudaq.set_target(\"density-matrix-cpu\")\n",
+        "cudaq.set_random_seed(42)"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "### 2. Define the circuit and noise model\n",
+        "\n",
+        "Define a kernel and attach a noise model. Each gate you add to the noise model becomes a **noise site** when that gate appears in the circuit. For **single-qubit** gates (e.g. H) use **DepolarizationChannel** with one qubit; for **CNOTs** pass the **qubit pair** `[control, target]` and use **Depolarization2** (two-qubit depolarization). Here we use a small Bell-style circuit with depol on the Hadamard and on the controlled-X (qubit pair [0, 1])."
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {},
+      "outputs": [],
+      "source": [
+        "@cudaq.kernel\n",
+        "def bell_with_noise():\n",
+        "    q = cudaq.qvector(2)\n",
+        "    h(q[0])\n",
+        "    x.ctrl(q[0], q[1])\n",
+        "    mz(q)\n",
+        "\n",
+        "noise = cudaq.NoiseModel()\n",
+        "noise.add_channel(\"h\", [0], cudaq.DepolarizationChannel(0.05))\n",
+        "noise.add_channel(\"x\", [0, 1], cudaq.Depolarization2(0.03))"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "### 3. Run PTSBE sampling\n",
+        "\n",
+        "Call `cudaq.ptsbe.sample()` with the kernel, noise model, and shot count. Optional arguments:\n",
+        "\n",
+        "- **sampling_strategy** — how trajectories are chosen:\n",
+        "  - **ExhaustiveSamplingStrategy()**: use all possible trajectories (every combination of Kraus outcomes per noise site).\n",
+        "  - **ProbabilisticSamplingStrategy(seed=...)**: sample trajectories randomly according to their probabilities; use a seed for reproducibility.\n",
+        "  - **OrderedSamplingStrategy()**: use the top-$k$ trajectories by probability (highest first), up to `max_trajectories`.\n",
+        "- **shot_allocation** — how shots are split across the chosen trajectories:\n",
+        "  - **PROPORTIONAL** (default): allocate shots in proportion to each trajectory’s probability.\n",
+        "  - **UNIFORM**: give each trajectory the same number of shots.\n",
+        "  - **LOW_WEIGHT_BIAS**: bias more shots toward low-weight (fewer errors) trajectories; optional `bias_strength` (default 2.0).\n",
+        "  - **HIGH_WEIGHT_BIAS**: bias more shots toward high-weight trajectories; optional `bias_strength` (default 2.0).\n",
+        "  Example: `ShotAllocationStrategy(type=cudaq.ptsbe.ShotAllocationType.UNIFORM)` or `ShotAllocationStrategy(type=cudaq.ptsbe.ShotAllocationType.LOW_WEIGHT_BIAS, bias_strength=5.0)`.\n",
+        "- **max_trajectories**: cap the number of trajectories (useful for large shot counts).\n",
+        "- **return_execution_data** (bool): If ``True``, the result includes trace instructions and per-trajectory data (``result.ptsbe_execution_data``); see section 6 at the end."
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {},
+      "outputs": [],
+      "source": [
+        "shots = 1000000\n",
+        "\n",
+        "strategy = cudaq.ptsbe.ProbabilisticSamplingStrategy(seed=42)\n",
+        "result = cudaq.ptsbe.sample(\n",
+        "    bell_with_noise,\n",
+        "    noise_model=noise,\n",
+        "    shots_count=shots,\n",
+        "    sampling_strategy=strategy,\n",
+        ")\n",
+        "\n",
+        "print(\"PTSBE sample result:\")\n",
+        "print(result)\n",
+        "print(f\"Total shots: {result.get_total_shots()}\")"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "### 4. Compare with standard (density-matrix) sampling\n",
+        "\n",
+        "To verify that PTSBE matches the usual noisy simulation, run standard `cudaq.sample()` with the same kernel and noise model. With enough shots, the two outcome distributions should be close (see the **PTSBE accuracy validation** example for a Hellinger fidelity comparison)."
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {},
+      "outputs": [],
+      "source": [
+        "result_standard = cudaq.sample(bell_with_noise, noise_model=noise, shots_count=shots)\n",
+        "print(\"Standard (density-matrix) sample result:\")\n",
+        "print(result_standard)\n",
+        "print(f\"Total shots: {result_standard.get_total_shots()}\")"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "### 5. Return execution data\n",
+        "\n",
+        "Pass `return_execution_data=True` to get the PTSBE execution data: the trace (gate, noise, measurement instructions) and the list of trajectories with their probabilities and shot allocations. Use `result.ptsbe_execution_data` and `result.has_execution_data()`."
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {},
+      "outputs": [],
+      "source": [
+        "result_with_data = cudaq.ptsbe.sample(\n",
+        "    bell_with_noise,\n",
+        "    noise_model=noise,\n",
+        "    shots_count=shots,\n",
+        "    sampling_strategy=strategy,\n",
+        "    return_execution_data=True,\n",
+        ")\n",
+        "assert result_with_data.has_execution_data()\n",
+        "data = result_with_data.ptsbe_execution_data\n",
+        "\n",
+        "Gate = cudaq.ptsbe.TraceInstructionType.Gate\n",
+        "Noise = cudaq.ptsbe.TraceInstructionType.Noise\n",
+        "Measurement = cudaq.ptsbe.TraceInstructionType.Measurement\n",
+        "print(\"Execution data:\")\n",
+        "print(f\"  Instructions: {len(data.instructions)} total\")\n",
+        "n_gate = sum(1 for i in data.instructions if i.type == Gate)\n",
+        "n_noise = sum(1 for i in data.instructions if i.type == Noise)\n",
+        "n_meas = sum(1 for i in data.instructions if i.type == Measurement)\n",
+        "print(f\"    Gates: {n_gate}, Noise: {n_noise}, Measurements: {n_meas}\")\n",
+        "print(f\"  Trajectories: {len(data.trajectories)}\")\n",
+        "total_traj_shots = sum(t.num_shots for t in data.trajectories)\n",
+        "print(f\"  Sum of trajectory shots: {total_traj_shots} (expected {shots})\")\n",
+        "print(\"  Trace instructions (first 5):\")\n",
+        "for i, inst in enumerate(data.instructions[:5]):\n",
+        "    print(f\"    [{i}] type={inst.type}, name={inst.name}, targets={list(inst.targets)}\")\n",
+        "if data.trajectories:\n",
+        "    t0 = data.trajectories[0]\n",
+        "    print(f\"  Example trajectory: id={t0.trajectory_id}, probability={t0.probability:.6f}, num_shots={t0.num_shots}\")\n",
+        "    print(\"  Selected trajectory (kraus_selections):\")\n",
+        "    for sel in t0.kraus_selections:\n",
+        "        print(f\"    circuit_location={sel.circuit_location}, kraus_operator_index={sel.kraus_operator_index}, is_error={sel.is_error}\")"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "### 6. Two API options:\n",
+        "\n",
+        "We have two ways to sample from a noisy circuit:\n",
+        "\n",
+        "1. **`cudaq.sample(kernel, noise_model=noise, use_ptsbe=..., ...)`** — The main sample API with a **`use_ptsbe`** option. When `use_ptsbe=False` (or omitted), sampling uses the standard density-matrix path; when `use_ptsbe=True`, sampling uses the PTSBE (trajectory) path.\n",
+        "2. **`cudaq.ptsbe.sample(kernel, noise_model=noise, ...)`** — The dedicated PTSBE API. There is no `use_ptsbe` option; sampling always uses the PTSBE path (trajectories, strategies, execution data, etc.).\n",
+        "\n",
+        "**When to choose which:** Use **`cudaq.sample`** with `use_ptsbe=True` when you want PTSBE from the main sample API (one call, one place for all sampling). Use **`cudaq.ptsbe.sample`** when you want the dedicated PTSBE API (explicit PTSBE entry point, no flag). Both provide the same PTSBE path (large shot counts, batching, trajectory strategies, execution data)."
+      ]
+    }
+  ],
+  "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.12.3"
+    }
+  },
+  "nbformat": 4,
+  "nbformat_minor": 4
+}

docs/sphinx/using/examples/examples.rst

@@ -19,6 +19,7 @@ Examples that illustrate how to use CUDA-Q for application development are avail
       Multi-GPU Workflows <multi_gpu_workflows.rst>
       Optimizers & Gradients  <../../examples/python/optimizers_gradients.ipynb>
       Noisy Simulations <../../examples/python/noisy_simulations.ipynb>
+      PTSBE End-to-End Workflow <../../examples/python/ptsbe_end_to_end_workflow.ipynb>
       Constructing Operators <operators.rst>
       Performance Optimizations <../../examples/python/performance_optimizations.ipynb>
       Using Quantum Hardware Providers <hardware_providers.rst>

python/cudaq/__init__.py

@@ -141,6 +141,7 @@ def _configure_cuda_library_paths() -> None:
 from .runtime.observe import observe_async
 from .runtime.run import run
 from .runtime.run import run_async
+from .runtime import ptsbe
 from .runtime.translate import translate
 from .runtime.state import (get_state, get_state_async, to_cupy)
 from .runtime.draw import draw