♻️ Delegate exact circuit execution to MQT Core

CHANGELOG.md

@@ -12,6 +12,9 @@ releases may include breaking changes.
 
 ### Changed
 
+- 🚚 Move the high-level `sample` and `simulate_statevector` helpers from MQT
+  Core to `mqt.ddsim`, and consolidate exact DD circuit execution in Core
+  ([#978]) ([**@simon1hofmann**])
 - ♻️ Migrate circuit normalization to the MQT Core v4 API ([#977])
   ([**@simon1hofmann**])
 - 💥 Drop support for x86 macOS and stop publishing the respective wheels
@@ -176,6 +179,7 @@ _📚 Refer to the [GitHub Release Notes] for previous changelogs._
 
 <!-- PR links -->
 
+[#978]: https://github.com/munich-quantum-toolkit/ddsim/pull/978
 [#977]: https://github.com/munich-quantum-toolkit/ddsim/pull/977
 [#976]: https://github.com/munich-quantum-toolkit/ddsim/pull/976
 [#975]: https://github.com/munich-quantum-toolkit/ddsim/pull/975

UPGRADING.md

@@ -6,6 +6,13 @@ of changes including minor and patch releases, please refer to the
 
 ## [Unreleased]
 
+### MQT Core v4 simulation helpers
+
+Import the high-level `sample` and `simulate_statevector` helpers from
+`mqt.ddsim` instead of `mqt.core.dd`. The package-aware Core DD primitives
+remain available for applications that manage their own `DDPackage` and input
+state.
+
 ### macOS support
 
 MQT DDSIM no longer supports x86 macOS. Use Apple silicon with macOS 13.3 or

docs/simulators/CircuitSimulator.md

@@ -77,6 +77,25 @@ print(result)
 As expected, the output distribution is approximately 50% for the states
 $|00\rangle$ and $|11\rangle$ each.
 
+For one-off simulations, the package also provides convenience functions:
+
+```{code-cell} ipython3
+from mqt.core.ir import QuantumComputation
+
+from mqt.ddsim import sample, simulate_statevector
+
+counts = sample(qc, shots=1024)
+
+unitary_qc = QuantumComputation(2)
+unitary_qc.h(0)
+unitary_qc.cx(0, 1)
+statevector = simulate_statevector(unitary_qc)
+```
+
+The `sample` helper supports terminal and mid-circuit measurements, resets, and
+classical control. The `simulate_statevector` helper is intended for unitary
+circuits and returns an independent NumPy array.
+
 If we would like to obtain the full state vector of the quantum circuit, we can
 query it after the simulation as follows:
 

include/CircuitSimulator.hpp

@@ -26,6 +26,7 @@
 #include <stdexcept>
 #include <string>
 #include <utility>
+#include <vector>
 
 struct ApproximationInfo {
   enum ApproximationStrategy : std::uint8_t { FidelityDriven, MemoryDriven };
@@ -130,10 +131,12 @@ class CircuitSimulator : public Simulator {
   std::size_t approximationRuns{0};
   long double finalFidelity{1.0L};
 
+  std::map<std::string, std::size_t> simulateWithHooks(std::size_t shots);
+
   struct CircuitAnalysis {
     bool isDynamic = false;
     bool hasMeasurements = false;
-    std::map<qc::Qubit, size_t> measurementMap;
+    std::vector<std::pair<qc::Qubit, std::size_t>> terminalMeasurements;
   };
 
   CircuitAnalysis analyseCircuit();

include/DeterministicNoiseSimulator.hpp

@@ -98,6 +98,10 @@ class DeterministicNoiseSimulator : public CircuitSimulator {
         shots);
   }
 
+  std::map<std::string, std::size_t> simulate(std::size_t shots) override {
+    return simulateWithHooks(shots);
+  }
+
   void initializeSimulation(std::size_t nQubits) override;
   char measure(dd::Qubit i) override;
   void reset(qc::NonUnitaryOperation* nonUnitaryOp) override;
@@ -131,4 +135,5 @@ class DeterministicNoiseSimulator : public CircuitSimulator {
   double measurementThreshold = 0.01;
   dd::ddsim::DensityDDPackage densityDD;
   dd::ddsim::DeterministicNoiseFunctionality deterministicNoiseFunctionality;
+  bool densityStateInitialized = false;
 };

python/mqt/ddsim/__init__.py

@@ -28,6 +28,7 @@ def _dll_patch() -> None:
     _dll_patch()
     del _dll_patch
 
+from ._simulation import sample, simulate_statevector
 from ._version import version as __version__
 from .provider import DDSIMProvider
 from .pyddsim import (
@@ -56,4 +57,6 @@ def _dll_patch() -> None:
     "UnitarySimulator",
     "UnitarySimulatorMode",
     "__version__",
+    "sample",
+    "simulate_statevector",
 ]

python/mqt/ddsim/_simulation.py

@@ -0,0 +1,66 @@
+# Copyright (c) 2023 - 2026 Chair for Design Automation, TUM
+# Copyright (c) 2025 - 2026 Munich Quantum Software Company GmbH
+# All rights reserved.
+#
+# SPDX-License-Identifier: MIT
+#
+# Licensed under the MIT License
+
+"""High-level decision-diagram simulation helpers."""
+
+from __future__ import annotations
+
+from typing import TYPE_CHECKING
+
+from mqt.core.dd import DDPackage
+from mqt.core.dd import simulate as simulate_dd
+
+from .pyddsim import CircuitSimulator
+
+if TYPE_CHECKING:
+    import numpy as np
+    from mqt.core.ir import QuantumComputation
+    from numpy.typing import NDArray
+
+_MAX_SEED = (1 << 63) - 1
+
+
+def sample(qc: QuantumComputation, shots: int = 1024, seed: int = 0) -> dict[str, int]:
+    """Sample from the output distribution of a quantum computation.
+
+    Circuits with only terminal measurements are executed once. Circuits with
+    mid-circuit measurements, resets, or classical control are executed once
+    per shot.
+
+    Args:
+        qc: The quantum computation to simulate.
+        shots: The number of samples to draw.
+        seed: The non-negative signed 64-bit random seed. Zero selects a random
+            seed for compatibility with the former `mqt.core.dd.sample` helper.
+
+    Returns:
+        A histogram mapping big-endian measurement bitstrings to counts.
+    """
+    if not 0 <= seed <= _MAX_SEED:
+        msg = f"seed must be between 0 and {_MAX_SEED}"
+        raise ValueError(msg)
+
+    simulator = CircuitSimulator(qc, seed=-1 if seed == 0 else seed)
+    return simulator.simulate(shots)
+
+
+def simulate_statevector(qc: QuantumComputation) -> NDArray[np.complex128]:
+    """Simulate a unitary quantum computation and return its state vector.
+
+    Args:
+        qc: The unitary quantum computation to simulate.
+
+    Returns:
+        A copy of the final state vector.
+    """
+    package = DDPackage(qc.num_qubits)
+    final_state = simulate_dd(qc, package.zero_state(qc.num_qubits), package)
+    try:
+        return final_state.get_vector()
+    finally:
+        package.dec_ref_vec(final_state)

src/CircuitSimulator.cpp

@@ -14,6 +14,7 @@
 #include "dd/FunctionalityConstruction.hpp"
 #include "dd/Node.hpp"
 #include "dd/Operations.hpp"
+#include "dd/Simulation.hpp"
 #include "dd/StateGeneration.hpp"
 #include "ir/Definitions.hpp"
 #include "ir/operations/IfElseOperation.hpp"
@@ -28,9 +29,34 @@
 #include <memory>
 #include <stdexcept>
 #include <string>
+#include <utility>
 
 std::map<std::string, std::size_t>
 CircuitSimulator::simulate(std::size_t shots) {
+  // Derived simulators customize the hook-based execution path. Approximation
+  // also depends on its established per-operation scheduling.
+  if (typeid(*this) != typeid(CircuitSimulator) ||
+      (approximationInfo.stepNumber > 0U &&
+       approximationInfo.stepFidelity < 1.0)) {
+    return simulateWithHooks(shots);
+  }
+
+  const auto& roots = dd->getRootSet<dd::vNode>();
+  const auto previousRoot = rootEdge;
+  const auto previousRootIsRegistered = roots.contains(previousRoot);
+
+  auto result =
+      dd::sample(*qc, dd::makeZeroState(qc->getNqubits(), *dd), *dd, shots, mt);
+  if (previousRootIsRegistered) {
+    dd->decRef(previousRoot);
+  }
+  rootEdge = result.state;
+  singleShots += result.executions;
+  return std::move(result.counts);
+}
+
+std::map<std::string, std::size_t>
+CircuitSimulator::simulateWithHooks(std::size_t shots) {
   const auto analysis = analyseCircuit();
 
   // easiest case: all gates are unitary --> simulate once and sample away on
@@ -53,7 +79,8 @@
     for (const auto& [bit_string, count] : measureAllNonCollapsing(shots)) {
       std::string resultString(qc->getNcbits(), '0');
 
-      for (auto const& [qubit_index, bitIndex] : analysis.measurementMap) {
+      for (const auto& [qubit_index, bitIndex] :
+           analysis.terminalMeasurements) {
         resultString[cbits - bitIndex - 1] =
             bit_string[qubits - qubit_index - 1];
       }
@@ -102,7 +129,8 @@
       }
 
       for (unsigned int i = 0; i < quantum.size(); ++i) {
-        analysis.measurementMap[quantum.at(i)] = classic.at(i);
+        analysis.terminalMeasurements.emplace_back(quantum.at(i),
+                                                   classic.at(i));
       }
     }
 
@@ -126,6 +154,10 @@
 }
 
 void CircuitSimulator::initializeSimulation(const std::size_t nQubits) {
+  const auto& roots = dd->getRootSet<dd::vNode>();
+  if (roots.contains(rootEdge)) {
+    dd->decRef(rootEdge);
+  }
   rootEdge = dd::makeZeroState(static_cast<dd::Qubit>(nQubits), *dd);
 }
 

src/DeterministicNoiseSimulator.cpp

@@ -10,6 +10,7 @@
 
 #include "DeterministicNoiseSimulator.hpp"
 
+#include "DensityNode.hpp"
 #include "Simulator.hpp"
 #include "dd/ComplexNumbers.hpp"
 #include "dd/DDDefinitions.hpp"
@@ -31,7 +32,12 @@ using CN = dd::ComplexNumbers;
 
 void DeterministicNoiseSimulator::initializeSimulation(
     const std::size_t nQubits) {
+  if (densityStateInitialized) {
+    dd::ddsim::DensityMatrixDD::alignDensityEdge(rootEdge);
+    densityDD.decRef(rootEdge);
+  }
   rootEdge = densityDD.makeZeroDensityOperator(static_cast<dd::Qubit>(nQubits));
+  densityStateInitialized = true;
 }
 
 void DeterministicNoiseSimulator::applyOperationToState(

test/CMakeLists.txt

@@ -9,6 +9,7 @@
 package_add_test(
   mqt-ddsim-test
   MQT::DDSim
+  test_circuit_execution.cpp
   test_circuit_sim.cpp
   test_shor_sim.cpp
   test_fast_shor_sim.cpp

test/python/simulator/test_standalone_simulator.py

@@ -12,10 +12,11 @@
 import unittest
 
 import numpy as np
+import pytest
 from mqt.core import load
 from mqt.core.ir import QuantumComputation
 
-from mqt.ddsim import CircuitSimulator
+from mqt.ddsim import CircuitSimulator, sample, simulate_statevector
 
 
 class MQTStandaloneSimulatorTests(unittest.TestCase):
@@ -56,6 +57,34 @@ def test_standalone_with_seed(self) -> None:
         assert "000" in result
         assert "111" in result
 
+    @staticmethod
+    def test_sampling_helper() -> None:
+        circ = QuantumComputation(1, 2)
+        circ.x(0)
+        circ.measure(0, 0)
+        circ.measure(0, 1)
+
+        assert sample(circ, shots=16, seed=1337) == {"11": 16}
+
+        with pytest.raises(ValueError, match="seed must be between 0 and"):
+            sample(circ, shots=16, seed=-1)
+
+    @staticmethod
+    def test_statevector_helper() -> None:
+        circ = QuantumComputation(1)
+        circ.x(0)
+        circ.gphase(np.pi / 2)
+
+        assert np.allclose(simulate_statevector(circ), np.array([0, 1j]))
+
+    @staticmethod
+    def test_statevector_helper_rejects_nonunitary_circuits() -> None:
+        circ = QuantumComputation(1, 1)
+        circ.measure(0, 0)
+
+        with pytest.raises(ValueError, match="DD for non-unitary operation not available"):
+            simulate_statevector(circ)
+
     def test_standalone_simple_approximation(self) -> None:
 
         # creates a state with <2% probability of measuring |1x>