♻️ Consolidate DD circuit execution

CHANGELOG.md

@@ -79,6 +79,9 @@ releases may include breaking changes.
 
 #### Other additions
 
+- ✨ Add a low-level package-aware `dd::sample` overload that accepts a
+  caller-owned random-number generator and returns counts, the retained state,
+  and the execution count ([#2273]) ([**@simon1hofmann**])
 - 🐳 Add dev container configuration for a consistent local development
   environment ([#1786]) ([**@denialhaag**])
 
@@ -113,6 +116,9 @@ releases may include breaking changes.
 
 ### Removed
 
+- 💥 Move the high-level `mqt.core.dd.sample` and
+  `mqt.core.dd.simulate_statevector` helpers to MQT DDSIM and make virtual DD
+  execution internal ([#2273]) ([**@simon1hofmann**])
 - 💥 Remove `CircuitOptimizer`. Move circuit flattening and final-measurement
   removal to `QuantumComputation`, equivalence-checking transformations to MQT
   QCEC, and mapping transformations to MQT QMAP. Move single-qubit gate fusion
@@ -851,6 +857,7 @@ for previous changelogs._
 
 <!-- PR links -->
 
+[#2273]: https://github.com/munich-quantum-toolkit/core/pull/2273
 [#2262]: https://github.com/munich-quantum-toolkit/core/pull/2262
 [#2259]: https://github.com/munich-quantum-toolkit/core/pull/2259
 [#2258]: https://github.com/munich-quantum-toolkit/core/pull/2258

UPGRADING.md

@@ -75,6 +75,18 @@ factories. Move required implementations to the package that uses them. The
 `BUILD_MQT_CORE_BENCHMARKS` option and its legacy DD evaluation target were also
 removed.
 
+### Circuit simulation helpers
+
+The high-level Python `mqt.core.dd.sample` and
+`mqt.core.dd.simulate_statevector` helpers moved to `mqt.ddsim`. Import the same
+names from `mqt.ddsim>=2.6.0` instead. The package-aware C++ `dd::sample` and
+`dd::simulate` primitives and the Python `mqt.core.dd.simulate` binding remain
+available for callers that manage a DD package and input state directly.
+
+The public `dd::isExecutableVirtually` and `dd::applyVirtualOperation` helpers
+were removed. Virtual execution is an internal detail of Core's circuit
+simulation and has no public replacement.
+
 ### Python 3.11 and split-mode wheels
 
 MQT Core now requires Python 3.11 or newer. Upgrade the Python environment

bindings/dd/register_dd.cpp

@@ -13,7 +13,6 @@
 #include "dd/Node.hpp"
 #include "dd/Package.hpp"
 #include "dd/Simulation.hpp"
-#include "dd/StateGeneration.hpp"
 #include "ir/QuantumComputation.hpp"
 
 #include <nanobind/nanobind.h>
@@ -55,61 +54,6 @@ NB_MODULE(MQT_CORE_MODULE_NAME, m) {
   // DD Package
   registerDDPackage(m);
 
-  m.def(
-      "sample",
-      [](const qc::QuantumComputation& qc, const size_t shots = 1024U,
-         const size_t seed = 0U) { return dd::sample(qc, shots, seed); },
-      "qc"_a, "shots"_a = 1024U, "seed"_a = 0U,
-      R"pb(Sample from the output distribution of a quantum computation.
-
-This function classically simulates the quantum computation and repeatedly samples from the output distribution.
-It supports mid-circuit measurements, resets, and classical control.
-
-Args:
-    qc: The quantum computation.
-    shots: The number of samples to take.
-        If the quantum computation contains no mid-circuit measurements or resets, the circuit is simulated once and the samples are drawn from the final state.
-        Otherwise, the circuit is simulated once for each sample.
-        Defaults to 1024.
-    seed: The seed for the random number generator.
-        If set to a specific non-zero value, the simulation is deterministic.
-        If set to 0, the RNG is randomly seeded.
-        Defaults to 0.
-
-Returns:
-    A histogram of the samples.
-    Each sample is a bitstring representing the measurement outcomes of the qubits in the quantum computation.
-    The leftmost bit corresponds to the most significant qubit, that is, the qubit with the highest index (big-endian).
-    If the circuit contains measurements, only the qubits that are actively measured are included in the output distribution.
-    Otherwise, all qubits in the circuit are measured.)pb");
-
-  m.def(
-      "simulate_statevector",
-      [](const qc::QuantumComputation& qc) {
-        const auto dd = std::make_unique<dd::Package>(qc.getNqubits());
-        const auto in = makeZeroState(qc.getNqubits(), *dd);
-        const auto sim = simulate(qc, in, *dd);
-        return getVector(sim);
-      },
-      "qc"_a,
-      R"pb(Simulate the quantum computation and return the final state vector.
-
-This function classically simulates the quantum computation and returns the state vector of the final state.
-It does not support measurements, resets, or classical control.
-
-Since the state vector is guaranteed to be exponentially large in the number of qubits, this function is only suitable for small quantum computations.
-Consider using the :func:`~mqt.core.dd.simulate` or the :func:`~mqt.core.dd.sample` functions, which never explicitly construct the state vector, for larger quantum computations.
-
-Notes:
-    This function internally constructs a :class:`~mqt.core.dd.DDPackage`, creates the zero state, and simulates the quantum computation via the :func:`simulate` function.
-    The state vector is then extracted from the resulting DD via the :meth:`~mqt.core.dd.VectorDD.get_vector` method.
-
-Args:
-    qc: The quantum computation. Must only contain unitary operations.
-
-Returns:
-    The state vector of the final state.)pb");
-
   m.def(
       "build_unitary",
       [](const qc::QuantumComputation& qc) {
@@ -140,8 +84,7 @@ Consider using the :func:`~mqt.core.dd.build_functionality` function, which neve
         R"pb(Simulate a quantum computation.
 
 This function classically simulates a quantum computation for a given initial state and returns the final state (represented as a DD).
-Compared to the `sample` function, this function does not support measurements, resets, or classical control.
-It only supports unitary operations.
+This function only supports unitary operations; it does not support measurements, resets, or classical control.
 
 The simulation is effectively computed by sequentially applying the operations of the quantum computation to the initial state.
 

docs/dd_package.md

@@ -28,134 +28,22 @@ to work with decision diagrams in MQT Core from Python.
 
 ## Quickstart
 
-In its simplest use case, the MQT Core DD package can be used as a classical
-circuit simulator using the {py:func}`~mqt.core.dd.sample` function. The
-underlying simulation approach supports mid-circuit measurements, reset
-operations, as well as classically-controlled operations. For example, the
-following code snippet demonstrates how to simulate the iterative quantum phase
-estimation algorithm shown in
-[the MQT Core IR Quickstart guide](mqt_core_ir).
+MQT DDSIM 2.6 or newer owns high-level decision-diagram circuit simulation. Use
+its `mqt.ddsim.sample` helper for sampling—including dynamic circuits—and
+`mqt.ddsim.simulate_statevector` for a full state vector. MQT Core exposes the
+lower-level DD package operations used to build specialized tools.
 
-```{code-cell} ipython3
-from mqt.core.dd import sample
-from mqt.core.ir import QuantumComputation
-from mqt.core.ir.operations import OpType
-
-from math import pi
-
-theta = 3 * pi / 8
-precision = 3
-
-# Create an empty quantum computation
-qc = QuantumComputation()
-
-# Counting register
-q = qc.add_qubit_register(1, "q")
-
-# Eigenstate register
-psi = qc.add_qubit_register(1, "psi")
-
-# Classical register for the result, the estimated phase is `0.c_2 c_1 c_0 * pi`
-c = qc.add_classical_register(precision, "c")
-
-# Prepare psi in the eigenstate |1>
-qc.x(psi[0])
-
-for i in range(precision):
-  # Hadamard on the working qubit
-  qc.h(q[0])
-
-  # Controlled phase gate
-  qc.cp(2**(precision - i - 1) * theta, q[0], psi[0])
-
-  # Iterative inverse QFT
-  for j in range(i):
-    qc.if_(op_type=OpType.p, target=q[0], control_bit=c[j], params=[-pi / 2**(i - j)])
-  qc.h(q[0])
-
-  # Measure the result
-  qc.measure(q[0], c[i])
-
-  # Reset the qubit if not finished
-  if i < precision - 1:
-    qc.reset(q[0])
-
-# Run the simulation
-counts = sample(qc, 1024)
-```
-
-```{code-cell} ipython3
----
-tags: [remove-cell]
----
-from pathlib import Path
-from matplotlib import pyplot as plt
-
-def generate_plot(counts: dict[str, int], name: str, light: bool) -> None:
-    if light:
-        plt.style.use('default')
-    else:
-        plt.style.use('dark_background')
-
-    # Create the bar plot
-    fig, ax = plt.subplots()
-    bars = ax.bar(counts.keys(), counts.values(), color='#0065bd')
-
-    # Annotate counts above the bars
-    for bar in bars:
-        height = bar.get_height()
-        ax.annotate(f'{height}',
-                    xy=(bar.get_x() + bar.get_width() / 2, height),
-                    xytext=(0, 3),  # 3 points vertical offset
-                    textcoords="offset points",
-                    ha='center', va='bottom')
-
-    # Set background to transparent
-    fig.patch.set_alpha(0.0)
-    ax.patch.set_alpha(0.0)
-
-    # Remove top and right borders
-    ax.spines['top'].set_visible(False)
-    ax.spines['right'].set_visible(False)
-
-    plt.xlabel("Measurement Outcome")
-    plt.ylabel("Counts")
-
-    # export to SVG (ensure the directory exists)
-    Path("_build/html/_images").mkdir(parents=True, exist_ok=True)
-    filename = "_build/html/_images/fig-" + name + ("-light" if light else "-dark") + ".svg"
-    plt.savefig(filename, format="svg")
-
-name = 'qpe'
-generate_plot(counts, name, light=True)
-generate_plot(counts, name, light=False)
-```
-
-```{raw} html
-<img src="_images/fig-qpe-light.svg" class="only-light" style="display: block; margin: auto; width: 75%;" alt="QPE measurement counts">
-<img src="_images/fig-qpe-dark.svg" class="only-dark" style="display: block; margin: auto; width: 75%;" alt="QPE measurement counts">
-```
-
-The {py:func}`~mqt.core.dd.sample` function is a high-level interface to the
-decision diagram package that does not require any knowledge of the underlying
-data structure. In a similar fashion, the
-{py:func}`~mqt.core.dd.simulate_statevector` and
-{py:func}`~mqt.core.dd.build_unitary` functions can be used to perform
-statevector simulation or to construct the unitary matrix representation of a
-quantum circuit, respectively.
+The {py:func}`~mqt.core.dd.build_unitary` function constructs the unitary matrix
+representation of a circuit:
 
 ```{code-cell} ipython3
-from mqt.core.dd import simulate_statevector
-
 import numpy as np
 
+from mqt.core.ir import QuantumComputation
+
 qc = QuantumComputation(2)
 qc.h(0)
 qc.cx(0, 1)
-
-vec = np.array(simulate_statevector(qc), copy=False)
-with np.printoptions(precision=3, suppress=True):
-  print(vec)
 ```
 
 ```{code-cell} ipython3
@@ -166,15 +54,14 @@ with np.printoptions(precision=3, suppress=True):
   print(unitary)
 ```
 
-Both of these functions are inherently limited in their scalability due to the
-exponential growth of the resulting data structures. MQT Core also allows one to
-work with decision diagrams directly, which is particularly useful for larger
-quantum circuits. To this end, the {py:class}`~mqt.core.dd.DDPackage` class
-provides a low-level interface to the decision diagram package. An instance of
-this class can be used to simulate quantum circuits (see
-{py:func}`~mqt.core.dd.simulate`), construct unitary matrices (see
-{py:func}`~mqt.core.dd.build_functionality`), or perform other operations on
-decision diagrams.
+This function is inherently limited in its scalability due to the exponential
+growth of the resulting data structures. MQT Core also allows one to work with
+decision diagrams directly, which is particularly useful for larger quantum
+circuits. To this end, the {py:class}`~mqt.core.dd.DDPackage` class provides a
+low-level interface to the decision diagram package. An instance of this class
+can be used to simulate quantum circuits (see {py:func}`~mqt.core.dd.simulate`),
+construct unitary matrices (see {py:func}`~mqt.core.dd.build_functionality`), or
+perform other operations on decision diagrams.
 
 ```{code-cell} ipython3
 from mqt.core.dd import DDPackage, simulate

include/mqt-core/dd/Operations.hpp

@@ -206,32 +206,18 @@ VectorDD applyIfElseOperation(const qc::IfElseOperation& op, const VectorDD& in,
                               const std::vector<bool>& measurements,
                               const qc::Permutation& permutation = {});
 
-/**
- * @brief Check whether @p op is virtually executable.
- *
- * @param op The operation in question.
- * @return Whether @p op is virtually executable.
- */
-bool isExecutableVirtually(const qc::Operation& op) noexcept;
-
-/**
- * @brief Apply virtual operation @p op.
- *
- * @param op The virtual operation to apply.
- * @param permutation If suitable, the to be updated permutation.
- */
-void applyVirtualOperation(const qc::Operation& op,
-                           qc::Permutation& permutation) noexcept;
-
 /**
  * @brief Apply global phase to a given DD.
  *
+ * @details The registered root reference owned by @p in is transferred to this
+ * function. The returned DD owns exactly one registered root reference.
+ *
  * @param in The input DD
  * @param phase The phase to apply
  * @param dd The DD package to use
  * @return The output DD
  */
-VectorDD applyGlobalPhase(VectorDD& in, const fp& phase, Package& dd);
+VectorDD applyGlobalPhase(const VectorDD& in, const fp& phase, Package& dd);
 
 /**
  * @brief Change the permutation of a given DD.

include/mqt-core/dd/Simulation.hpp

@@ -18,6 +18,7 @@
 
 #include <cstddef>
 #include <map>
+#include <random>
 #include <string>
 
 namespace qc {
@@ -26,6 +27,13 @@ class QuantumComputation;
 
 namespace dd {
 
+/** @brief Result of sampling a quantum computation. */
+struct SamplingResult {
+  std::map<std::string, std::size_t> counts{};
+  VectorDD state{};
+  std::size_t executions = 0U;
+};
+
 /**
  * @brief Simulate a purely-quantum @ref qc::QuantumComputation on a given input
  * state using decision diagrams.
@@ -49,6 +57,32 @@ namespace dd {
 VectorDD simulate(const qc::QuantumComputation& qc, const VectorDD& in,
                   Package& dd);
 
+/**
+ * @brief Execute and sample a quantum computation using a caller-owned state,
+ * package, and random-number generator.
+ *
+ * @details This is the low-level execution primitive underlying the sampling
+ * convenience functions. Static circuits are executed once and sampled
+ * without collapsing the retained state. Dynamic circuits are executed once
+ * per shot and retain the final state from the last execution.
+ *
+ * The registered root reference owned by @p in is transferred to this
+ * function. The returned @ref SamplingResult::state owns exactly one
+ * registered root reference in @p dd and must eventually be passed to
+ * @ref Package::decRef. For a dynamic circuit with zero shots, the input state
+ * is returned unchanged and no execution is reported.
+ *
+ * @param qc The quantum computation to execute
+ * @param in The input state whose registered root reference is transferred
+ * @param dd The DD package to use for execution
+ * @param shots The number of samples to draw
+ * @param rng The random-number generator to use
+ * @return Counts, the retained final state, and the number of executions
+ */
+[[nodiscard]] SamplingResult sample(const qc::QuantumComputation& qc,
+                                    VectorDD in, Package& dd, std::size_t shots,
+                                    std::mt19937_64& rng);
+
 /**
  * @brief Sample from the output distribution of a quantum computation
  *
@@ -81,10 +115,11 @@ std::map<std::string, std::size_t> sample(const qc::QuantumComputation& qc,
  *
  * @details This is a more general version of @ref dd::sample that allows for
  * choosing the input state to simulate as well as the DD package to use for the
- * simulation.
+ * simulation. The registered root reference owned by @p in is transferred to
+ * this function.
  *
  * @param qc The quantum computation to simulate
- * @param in The input state to simulate. Represented as a vector DD.
+ * @param in The input state whose registered root reference is transferred.
  * @param dd The DD package to use for the simulation
  * @param shots The number of shots to sample
  * @param seed The seed for the random number generator