✨ Add QCO density-matrix DD simulation

CHANGELOG.md

@@ -29,10 +29,15 @@ releases may include breaking changes.
   [**@denialhaag**], [**@taminob**], [**@DRovara**], [**@li-mingbao**],
   [**@Ectras**], [**@MatthiasReumann**], [**@simon1hofmann**], [**@J4MMlE**])
 - ✨ Add decision diagram-based construction, simulation, and sampling for QCO
-  programs, including mid-circuit `measure` / `reset`, concrete `if` /
-  `index_switch` / `scf.for` / `func.call`, classical SSA and CBit registers,
-  dense multi-wire embedding, output-aware multi-shot sampling, and Python
-  bindings ([#1915], [#1973], [#2077]) ([**@simon1hofmann**])
+  programs, including mid-circuit `measure` / `reset`, concrete QCO, SCF, and
+  multi-block CFG control flow, non-recursive calls, bound parameters, classical
+  integer, floating-point, and common math operations, CBit registers,
+  one-dimensional memrefs, dynamic quantum allocation and separable
+  deallocation, qtensors, dense multi-wire embedding, Python input bindings,
+  optional sampling input states, output-aware multi-shot sampling,
+  density-matrix simulation and sampling with physical partial trace, and Python
+  bindings ([#1915], [#1973], [#2077], [#2078], [#2079], [#2080])
+  ([**@simon1hofmann**], [**@burgholzer**])
 - ✨ Add immutable MLIR compiler targets, QDMI device integration, and target
   compilation through C++, Python, and `mqt-cc` ([#1687], [#1993], [#1999],
   [#2049]) ([**@MatthiasReumann**], [**@simon1hofmann**], [**@burgholzer**])
@@ -131,8 +136,10 @@ releases may include breaking changes.
 - 💥 Remove `MQT::CoreAlgorithms`, its fixed-circuit factories, and the legacy
   DD package evaluation. MQT Core provides no direct replacement ([#2214])
   ([**@burgholzer**])
-- 💥 Remove the unowned decision-diagram approximation algorithm and
-  density-matrix support from MQT Core ([#1466], [#2154]) ([**@burgholzer**])
+- 💥 Remove the unowned decision-diagram approximation algorithm and the legacy
+  raw density-matrix and noise APIs. Compiler-backed QCO density simulation
+  remains available through generic matrix DDs ([#1466], [#2154])
+  ([**@burgholzer**])
 - 💥 Make `nlohmann_json` an implementation detail and replace JSON-typed
   decision-diagram statistics APIs with strings and streams ([#2138])
   ([**@denialhaag**])
@@ -914,6 +921,9 @@ for previous changelogs._
 [#2105]: https://github.com/munich-quantum-toolkit/core/pull/2105
 [#2084]: https://github.com/munich-quantum-toolkit/core/pull/2084
 [#2082]: https://github.com/munich-quantum-toolkit/core/pull/2082
+[#2080]: https://github.com/munich-quantum-toolkit/core/pull/2080
+[#2079]: https://github.com/munich-quantum-toolkit/core/pull/2079
+[#2078]: https://github.com/munich-quantum-toolkit/core/pull/2078
 [#2077]: https://github.com/munich-quantum-toolkit/core/pull/2077
 [#2074]: https://github.com/munich-quantum-toolkit/core/pull/2074
 [#2066]: https://github.com/munich-quantum-toolkit/core/pull/2066

UPGRADING.md

@@ -105,16 +105,19 @@ The Python bindings depend on `nanobind-backend`, which supplies the
 interpreter-specific nanobind runtime. This dependency does not change the C++
 API or the Python import paths.
 
-### Removal of DD approximation and density-matrix support
+### Removal of DD approximation and legacy density/noise APIs
 
 MQT Core no longer provides the decision-diagram approximation algorithm. The
 algorithm had no production owner in the MQT ecosystem. Remove uses of the
 `dd/Approximation.hpp` header, the `dd::ApproximationMetadata` type, and the
 `dd::approximate` function. MQT Core does not provide a replacement.
 
-MQT Core also no longer provides density-matrix decision diagrams or the noise
-operations that depended on them. Consumers must provide this functionality or
-use another implementation.
+MQT Core also no longer provides the legacy raw density-matrix DD types or the
+noise operations that depended on them. Compiler-backed QCO density simulation
+remains available through generic matrix DDs: use
+`mqt.core.mlir.make_density_matrix`, then call `QCOProgram.simulate_density` or
+`QCOProgram.sample_density`. Consumers of the removed raw APIs must migrate to
+QCO or another implementation.
 
 ### Private `nlohmann_json` dependency
 

bindings/mlir/register_mlir.cpp

@@ -8,12 +8,14 @@
  * Licensed under the MIT License
  */
 
+#include "dd/Edge.hpp"
 #include "dd/Node.hpp"
 #include "dd/Package.hpp"
 #include "mlir/Compiler/Programs.h"
 #include "mlir/Compiler/QDMIAdapter.h"
 #include "mlir/Compiler/Target.h"
 #include "mlir/Dialect/MQT/IR/MQTDialect.h"
+#include "mlir/Dialect/QCO/IR/QCODialect.h"
 #include "mlir/Dialect/QCO/Utils/DDFunctionality.h"
 #include "qdmi/Client.hpp"
 #include "qdmi/driver/SessionConfig.hpp"
@@ -22,6 +24,8 @@
 #include <llvm/Support/Error.h>
 #include <llvm/Support/raw_ostream.h>
 #include <mlir/Dialect/Func/IR/FuncOps.h>
+#include <mlir/IR/BuiltinAttributes.h>
+#include <mlir/IR/BuiltinTypes.h>
 #include <mlir/IR/Diagnostics.h>
 #include <mlir/IR/MLIRContext.h>
 #include <mlir/Support/LogicalResult.h>
@@ -39,6 +43,7 @@
 #include <cstddef>
 #include <cstdint>
 #include <filesystem>
+#include <map>
 #include <memory>
 #include <optional>
 #include <random>
@@ -48,6 +53,7 @@
 #include <string_view>
 #include <system_error>
 #include <utility>
+#include <variant>
 #include <vector>
 
 namespace mqt {
@@ -139,6 +145,66 @@ entryFunc(const mlir::QCOProgram& program) {
   return func;
 }
 
+using QCODDBindingValue = std::variant<bool, int64_t, double>;
+using QCODDBindingMap = std::map<size_t, QCODDBindingValue>;
+
+[[nodiscard]] static mlir::qco::DDBindings
+makeQCODDBindings(mlir::func::FuncOp func,
+                  const QCODDBindingMap& pythonBindings) {
+  mlir::qco::DDBindings bindings;
+  for (const auto& [index, binding] : pythonBindings) {
+    if (index >= func.getNumArguments()) {
+      throw nb::value_error("QCO DD binding argument index is out of range");
+    }
+
+    mlir::Value argument = func.getArgument(static_cast<unsigned>(index));
+    const mlir::Type type = argument.getType();
+    mlir::Attribute attribute;
+    if (type.isInteger(1)) {
+      if (const auto* value = std::get_if<bool>(&binding)) {
+        attribute = mlir::BoolAttr::get(func.getContext(), *value);
+      }
+    } else if (mlir::isa<mlir::IndexType, mlir::IntegerType>(type)) {
+      if (const auto* value = std::get_if<int64_t>(&binding)) {
+        attribute = mlir::IntegerAttr::get(type, *value);
+      }
+    } else if (const auto floatType = mlir::dyn_cast<mlir::FloatType>(type)) {
+      if (const auto* value = std::get_if<double>(&binding)) {
+        attribute = mlir::FloatAttr::get(floatType, *value);
+      }
+    } else if (const auto tensorType =
+                   mlir::dyn_cast<mlir::RankedTensorType>(type);
+               tensorType && tensorType.getRank() == 1 &&
+               tensorType.isDynamicDim(0) &&
+               mlir::isa<mlir::qco::QubitType>(tensorType.getElementType())) {
+      if (const auto* value = std::get_if<int64_t>(&binding);
+          value != nullptr && *value >= 0) {
+        attribute = mlir::IntegerAttr::get(
+            mlir::IndexType::get(func.getContext()), *value);
+      }
+    }
+
+    if (!attribute) {
+      throw nb::value_error(
+          "QCO DD binding value does not match the entry argument type");
+    }
+    bindings[argument] = attribute;
+  }
+  return bindings;
+}
+
+template <class Node>
+static void requireLiveReference(const dd::Edge<Node>& state,
+                                 dd::Package& ddPackage,
+                                 const char* argumentName = "initial_state") {
+  if (dd::Edge<Node>::trackingRequired(state) &&
+      !ddPackage.getRootSet<Node>().contains(state)) {
+    std::string message(argumentName);
+    message.append(" must have a live reference in dd_package");
+    throw nb::value_error(message.c_str());
+  }
+}
+
 [[nodiscard]] static std::mt19937_64 makeRng(const uint64_t seed) {
   if (seed == 0) {
     return std::mt19937_64(std::random_device{}());
@@ -993,20 +1059,24 @@ LLVM bitcode.)pb");
 
   qcoProgram.def(
       "build_functionality",
-      [](const mlir::QCOProgram& program, dd::Package& ddPackage) {
+      [](const mlir::QCOProgram& program, dd::Package& ddPackage,
+         const QCODDBindingMap& pythonBindings) {
         auto func = entryFunc(program);
+        const auto bindings = makeQCODDBindings(func, pythonBindings);
         return takeFailureOr(
             func.getContext(),
-            "cannot build DD functionality for this QCO program",
-            [&] { return mlir::qco::buildFunctionality(func, ddPackage); });
+            "cannot build DD functionality for this QCO program", [&] {
+              return mlir::qco::buildFunctionality(func, ddPackage, bindings);
+            });
       },
-      "dd_package"_a,
+      "dd_package"_a, nb::kw_only(), "bindings"_a = QCODDBindingMap{},
       // Keep the DD package alive while the returned matrix DD is alive.
       nb::keep_alive<0, 2>(),
       R"pb(Build a matrix DD for a static unitary QCO program.
 
 Args:
     dd_package: DD package with enough qubits for the program.
+    bindings: Concrete entry-argument values keyed by zero-based argument index.
 
 Returns:
     Matrix DD of the program functionality.
@@ -1017,20 +1087,20 @@ LLVM bitcode.)pb");
   qcoProgram.def(
       "simulate",
       [](const mlir::QCOProgram& program, const dd::VectorDD& initialState,
-         dd::Package& ddPackage, const uint64_t seed) {
-        if (dd::VectorDD::trackingRequired(initialState) &&
-            !ddPackage.getRootSet<dd::vNode>().contains(initialState)) {
-          throw nb::value_error(
-              "initial_state must have a live reference in dd_package");
-        }
+         dd::Package& ddPackage, const uint64_t seed,
+         const QCODDBindingMap& pythonBindings) {
+        requireLiveReference(initialState, ddPackage);
         auto func = entryFunc(program);
+        const auto bindings = makeQCODDBindings(func, pythonBindings);
         auto rng = makeRng(seed);
         return takeFailureOr(
             func.getContext(), "cannot simulate this QCO program", [&] {
-              return mlir::qco::simulate(func, initialState, ddPackage, rng);
+              return mlir::qco::simulate(func, initialState, ddPackage, rng,
+                                         bindings);
             });
       },
-      "initial_state"_a, "dd_package"_a, "seed"_a = 0U,
+      "initial_state"_a, "dd_package"_a, "seed"_a = 0U, nb::kw_only(),
+      "bindings"_a = QCODDBindingMap{},
       // Keep the DD package alive while the returned vector DD is alive.
       nb::keep_alive<0, 3>(),
       R"pb(Simulate a QCO program on a DD state.
@@ -1042,6 +1112,7 @@ LLVM bitcode.)pb");
     dd_package: DD package with enough qubits for the program.
     seed: RNG seed. ``0`` (default) selects nondeterministic seeding. Any other
         value produces reproducible measurement and reset results.
+    bindings: Concrete entry-argument values keyed by zero-based argument index.
 
 Returns:
     Output state DD.
@@ -1050,31 +1121,149 @@ LLVM bitcode.)pb");
     ValueError: When ``initial_state`` has no live reference in ``dd_package``,
         has too few qubits, or the program is unsupported for simulation.)pb");
 
+  m.def(
+      "make_density_matrix",
+      [](const dd::VectorDD& state, const size_t numQubits,
+         dd::Package& ddPackage) {
+        requireLiveReference(state, ddPackage, "state");
+        try {
+          return mlir::qco::makeDensityMatrix(state, numQubits, ddPackage);
+        } catch (const std::invalid_argument& error) {
+          throw nb::value_error(error.what());
+        }
+      },
+      "state"_a, "num_qubits"_a, "dd_package"_a, nb::keep_alive<0, 3>(),
+      R"pb(Construct ``|psi><psi|`` from a pure DD state.
+
+The input vector reference remains owned by the caller. The returned matrix DD
+has a live reference in ``dd_package``.
+
+Args:
+    state: Pure state with a live reference in ``dd_package``.
+    num_qubits: Number of active qubits represented by ``state``.
+    dd_package: DD package that owns ``state`` and has enough qubits.
+
+Returns:
+    Density-matrix DD for the pure state.
+
+Raises:
+    ValueError: When ``state`` has no live reference in ``dd_package`` or
+        ``num_qubits`` does not cover the state or exceeds the DD package
+        capacity.)pb");
+
+  qcoProgram.def(
+      "simulate_density",
+      [](const mlir::QCOProgram& program, const dd::MatrixDD& initialState,
+         dd::Package& ddPackage, const uint64_t seed,
+         const QCODDBindingMap& pythonBindings) {
+        requireLiveReference(initialState, ddPackage);
+        auto func = entryFunc(program);
+        const auto bindings = makeQCODDBindings(func, pythonBindings);
+        auto rng = makeRng(seed);
+        return takeFailureOr(
+            func.getContext(), "cannot density-simulate this QCO program", [&] {
+              return mlir::qco::simulateDensity(func, initialState, ddPackage,
+                                                rng, bindings);
+            });
+      },
+      "initial_state"_a, "dd_package"_a, "seed"_a = 0U, nb::kw_only(),
+      "bindings"_a = QCODDBindingMap{}, nb::keep_alive<0, 3>(),
+      R"pb(Simulate a QCO program on a density-matrix DD.
+
+Args:
+    initial_state: Input density-matrix DD with a live reference in
+        ``dd_package``. It represents exactly the program's inferred initial
+        quantum register; skipped DD levels denote identity factors within
+        that register. A valid input reference is consumed.
+    dd_package: DD package with enough qubits for the program.
+    seed: RNG seed. ``0`` (default) selects nondeterministic seeding. Any other
+        value produces reproducible measurement and reset results.
+    bindings: Concrete entry-argument values keyed by zero-based argument index.
+
+Returns:
+    Output density-matrix DD.
+
+Raises:
+    ValueError: When ``initial_state`` has no live reference in ``dd_package``
+        or the program is unsupported for simulation.)pb");
+
   qcoProgram.def(
       "sample",
       [](const mlir::QCOProgram& program, dd::Package& ddPackage,
-         const size_t shots, const uint64_t seed) {
+         const size_t shots, const uint64_t seed,
+         const std::optional<dd::VectorDD>& initialState,
+         const QCODDBindingMap& pythonBindings) {
         auto func = entryFunc(program);
+        const auto bindings = makeQCODDBindings(func, pythonBindings);
         auto rng = makeRng(seed);
         return takeFailureOr(
-            func.getContext(), "cannot sample this QCO program",
-            [&] { return mlir::qco::sample(func, ddPackage, shots, rng); });
+            func.getContext(), "cannot sample this QCO program", [&] {
+              if (initialState) {
+                requireLiveReference(*initialState, ddPackage);
+                return mlir::qco::sample(func, *initialState, ddPackage, shots,
+                                         rng, bindings);
+              }
+              return mlir::qco::sample(func, ddPackage, shots, rng, bindings);
+            });
       },
-      "dd_package"_a, "shots"_a = 1024U, "seed"_a = 0U,
+      "dd_package"_a, "shots"_a = 1024U, "seed"_a = 0U, nb::kw_only(),
+      "initial_state"_a = nb::none(), "bindings"_a = QCODDBindingMap{},
       R"pb(Sample the declared outputs of a QCO program.
 
 Args:
     dd_package: DD package with enough qubits for the program.
     shots: Number of shots (default 1024).
     seed: RNG seed. ``0`` (default) selects nondeterministic seeding. Any other
         value produces reproducible results.
+    initial_state: Optional input state with a live reference in ``dd_package``.
+        A valid input reference is consumed.
+    bindings: Concrete entry-argument values keyed by zero-based argument index.
+
+Returns:
+    Histogram of returned CBit registers in return order, each MSB first. If
+    no CBit result exists, final ``measureAll`` bitstrings instead.
+
+Raises:
+    ValueError: When ``initial_state`` has no live reference in ``dd_package``
+        or the program is unsupported for sampling.)pb");
+
+  qcoProgram.def(
+      "sample_density",
+      [](const mlir::QCOProgram& program, const dd::MatrixDD& initialState,
+         dd::Package& ddPackage, const size_t shots, const uint64_t seed,
+         const QCODDBindingMap& pythonBindings) {
+        requireLiveReference(initialState, ddPackage);
+        auto func = entryFunc(program);
+        const auto bindings = makeQCODDBindings(func, pythonBindings);
+        auto rng = makeRng(seed);
+        return takeFailureOr(
+            func.getContext(), "cannot density-sample this QCO program", [&] {
+              return mlir::qco::sampleDensity(func, initialState, ddPackage,
+                                              shots, rng, bindings);
+            });
+      },
+      "initial_state"_a, "dd_package"_a, "shots"_a = 1024U, "seed"_a = 0U,
+      nb::kw_only(), "bindings"_a = QCODDBindingMap{},
+      R"pb(Sample the declared outputs of a QCO program from a density-matrix DD.
+
+Args:
+    initial_state: Input density-matrix DD with a live reference in
+        ``dd_package``. It represents exactly the program's inferred initial
+        quantum register; skipped DD levels denote identity factors within
+        that register. A valid input reference is consumed.
+    dd_package: DD package with enough qubits for the program.
+    shots: Number of shots (default 1024).
+    seed: RNG seed. ``0`` (default) selects nondeterministic seeding. Any other
+        value produces reproducible results.
+    bindings: Concrete entry-argument values keyed by zero-based argument index.
 
 Returns:
     Histogram of returned CBit registers in return order, each MSB first. If
     no CBit result exists, final ``measureAll`` bitstrings instead.
 
 Raises:
-    ValueError: When the program is unsupported for sampling.)pb");
+    ValueError: When ``initial_state`` has no live reference in ``dd_package``
+        or the program is unsupported for sampling.)pb");
 
   m.def("compile_program", &compileProgram, "program"_a, nb::kw_only(),
         "output"_a = mlir::ProgramFormat::QC, "inplace"_a = false,