Add Stim simulator

.gitmodules

@@ -45,3 +45,6 @@
  path = tpls/Crow
  url = https://github.com/CrowCpp/Crow.git
  ignore = dirty
+[submodule "tpls/Stim"]
+ path = tpls/Stim
+ url = https://github.com/quantumlib/Stim

runtime/nvqir/CMakeLists.txt

@@ -38,6 +38,7 @@ install(TARGETS ${LIBRARY_NAME}
  INCLUDES DESTINATION include/nvqir)
 
 add_subdirectory(qpp)
+add_subdirectory(stim)
 
 if (CUSTATEVEC_ROOT AND CUDA_FOUND) 
  add_subdirectory(custatevec)

runtime/nvqir/stim/CMakeLists.txt

@@ -0,0 +1,42 @@
+# ============================================================================ #
+# Copyright (c) 2022 - 2024 NVIDIA Corporation & Affiliates. #
+# All rights reserved. #
+# #
+# This source code and the accompanying materials are made available under #
+# the terms of the Apache License 2.0 which accompanies this distribution. #
+# ============================================================================ #
+
+set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wno-ctad-maybe-unsupported")
+set(INTERFACE_POSITION_INDEPENDENT_CODE ON)
+
+set(STIM_SOURCE_DIR ${CMAKE_SOURCE_DIR}/tpls/Stim)
+set(STIM_BINARY_DIR ${CMAKE_BINARY_DIR}/tpls/Stim)
+
+add_subdirectory(${STIM_SOURCE_DIR} ${STIM_BINARY_DIR})
+
+macro (AddStimBackend LIBRARY_NAME SOURCE_FILE)
+ add_library(${LIBRARY_NAME} SHARED ${SOURCE_FILE})
+ set_property(GLOBAL APPEND PROPERTY CUDAQ_RUNTIME_LIBS ${LIBRARY_NAME})
+
+ set (STIM_DEPENDENCIES "")
+ list(APPEND STIM_DEPENDENCIES fmt::fmt-header-only cudaq-common)
+
+ target_include_directories(${LIBRARY_NAME}
+ PUBLIC
+ $<BUILD_INTERFACE:${CMAKE_CURRENT_SOURCE_DIR}>
+ $<BUILD_INTERFACE:${CMAKE_SOURCE_DIR}/runtime>
+ $<INSTALL_INTERFACE:include>)
+
+ target_link_libraries(${LIBRARY_NAME}
+ PUBLIC libstim
+ PRIVATE ${STIM_DEPENDENCIES})
+
+ set_target_properties(${LIBRARY_NAME}
+ PROPERTIES INSTALL_RPATH "${CMAKE_INSTALL_RPATH}:${LLVM_BINARY_DIR}/lib")
+
+ install(TARGETS ${LIBRARY_NAME} DESTINATION lib)
+endmacro()
+
+AddStimBackend(nvqir-stim StimCircuitSimulator.cpp)
+
+add_target_config(stim)

runtime/nvqir/stim/StimCircuitSimulator.cpp

@@ -0,0 +1,192 @@
+/*******************************************************************************
+ * Copyright (c) 2022 - 2024 NVIDIA Corporation & Affiliates. *
+ * All rights reserved. *
+ * *
+ * This source code and the accompanying materials are made available under *
+ * the terms of the Apache License 2.0 which accompanies this distribution. *
+ ******************************************************************************/
+
+#include "nvqir/CircuitSimulator.h"
+#include "nvqir/Gates.h"
+#include "stim.h"
+
+#include <bit>
+#include <iostream>
+#include <set>
+#include <span>
+
+using namespace cudaq;
+
+namespace nvqir {
+
+/// @brief The StimCircuitSimulator implements the CircuitSimulator
+/// base class to provide a simulator delegating to the Stim library from
+/// https://github.com/quantumlib/Stim.
+class StimCircuitSimulator : public nvqir::CircuitSimulatorBase<double> {
+protected:
+ stim::Circuit stimCircuit;
+ std::mt19937_64 randomEngine;
+
+ /// @brief Grow the state vector by one qubit.
+ void addQubitToState() override { addQubitsToState(1); }
+
+ /// @brief Override the default sized allocation of qubits
+ /// here to be a bit more efficient than the default implementation
+ void addQubitsToState(std::size_t qubitCount,
+ const void *stateDataIn = nullptr) override {
+ return;
+ }
+
+ /// @brief Reset the qubit state.
+ void deallocateStateImpl() override { stimCircuit.clear(); }
+
+ /// @brief Apply the noise channel on \p qubits
+ void applyNoiseChannel(const std::string_view gateName,
+ const std::vector<std::size_t> &qubits) override {
+ // Do nothing if no execution context
+ if (!executionContext)
+ return;
+
+ // Do nothing if no noise model
+ if (!executionContext->noiseModel)
+ return;
+
+ // Get the name as a string
+ std::string gName(gateName);
+
+ // Cast size_t to uint32_t
+ std::vector<std::uint32_t> stimTargets;
+ stimTargets.reserve(qubits.size());
+ for (auto q : qubits)
+ stimTargets.push_back(static_cast<std::uint32_t>(q));
+
+ // Get the Kraus channels specified for this gate and qubits
+ auto krausChannels =
+ executionContext->noiseModel->get_channels(gName, qubits);
+
+ // If none, do nothing
+ if (krausChannels.empty())
+ return;
+
+ // TODO
+ return;
+ }
+
+ void applyGate(const GateApplicationTask &task) override {
+ std::string gateName(task.operationName);
+ std::transform(gateName.begin(), gateName.end(), gateName.begin(),
+ ::toupper);
+ std::vector<std::uint32_t> stimTargets;
+
+ // These CUDA-Q rotation gates have the same name as Stim "reset" gates.
+ // Stim is a Clifford simulator, so it doesn't actually support rotational
+ // gates. Throw exceptions if they are encountered here.
+ if (gateName == "RX" || gateName == "RY" || gateName == "RZ")
+ throw std::runtime_error(fmt::format(
+ "Gate not supported by simulator: {}", task.operationName));
+
+ if (task.controls.size() > 1)
+ throw std::runtime_error(
+ "Gates with >1 controls not supported by stim simulator");
+ if (task.controls.size() >= 1)
+ gateName = "C" + gateName;
+ for (auto c : task.controls)
+ stimTargets.push_back(c);
+ for (auto t : task.targets)
+ stimTargets.push_back(t);
+ try {
+ stimCircuit.safe_append_u(gateName, stimTargets);
+ } catch (std::out_of_range &e) {
+ throw std::runtime_error(
+ fmt::format("Gate not supported by simulator: {}", e.what()));
+ }
+ }
+
+ /// @brief Set the current state back to the |0> state.
+ void setToZeroState() override { return; }
+
+ /// @brief Override the calculateStateDim because this is not a state vector
+ /// simulator.
+ std::size_t calculateStateDim(const std::size_t numQubits) override {
+ return 0;
+ }
+
+ /// @brief Measure the qubit and return the result. Collapse the
+ /// state vector.
+ bool measureQubit(const std::size_t index) override { return false; }
+
+ QubitOrdering getQubitOrdering() const override { return QubitOrdering::msb; }
+
+public:
+ StimCircuitSimulator() {
+ // Populate the correct name so it is printed correctly during
+ // deconstructor.
+ summaryData.name = name();
+ }
+ virtual ~StimCircuitSimulator() = default;
+
+ void setRandomSeed(std::size_t seed) override {
+ randomEngine = std::mt19937_64(seed);
+ }
+
+ bool canHandleObserve() override { return false; }
+
+ /// @brief Reset the qubit
+ /// @param index 0-based index of qubit to reset
+ void resetQubit(const std::size_t index) override {
+ flushGateQueue();
+ stimCircuit.safe_append_u(
+ "R", std::vector<std::uint32_t>{static_cast<std::uint32_t>(index)});
+ }
+
+ /// @brief Sample the multi-qubit state.
+ cudaq::ExecutionResult sample(const std::vector<std::size_t> &qubits,
+ const int shots) override {
+ std::vector<std::uint32_t> stimTargetQubits;
+ for (auto q : qubits)
+ stimTargetQubits.push_back(static_cast<std::uint32_t>(q));
+ stimCircuit.safe_append_u("M", stimTargetQubits);
+ if (false) {
+ std::stringstream ss;
+ ss << stimCircuit << '\n';
+ cudaq::log("Stim circuit is\n{}", ss.str());
+ }
+ auto ref_sample = stim::TableauSimulator<
+ stim::MAX_BITWORD_WIDTH>::reference_sample_circuit(stimCircuit);
+ stim::simd_bit_table<stim::MAX_BITWORD_WIDTH> sample =
+ stim::sample_batch_measurements(stimCircuit, ref_sample, shots,
+ randomEngine, false);
+ size_t bits_per_sample = stimCircuit.count_measurements();
+ std::vector<std::string> sequentialData;
+ sequentialData.reserve(shots);
+ // Only retain the final "qubits.size()" measurements. All other
+ // measurements were mid-circuit measurements that have been previously
+ // accounted for and saved.
+ assert(bits_per_sample >= qubits.size());
+ std::size_t first_bit_to_save = bits_per_sample - qubits.size();
+ CountsDictionary counts;
+ for (std::size_t shot = 0; shot < shots; shot++) {
+ std::string aShot(qubits.size(), '0');
+ for (std::size_t b = first_bit_to_save; b < bits_per_sample; b++) {
+ aShot[b - first_bit_to_save] = sample[b][shot] ? '1' : '0';
+ }
+ counts[aShot]++;
+ sequentialData.push_back(std::move(aShot));
+ }
+ ExecutionResult result(counts);
+ result.sequentialData = std::move(sequentialData);
+ return result;
+ }
+
+ bool isStateVectorSimulator() const override { return false; }
+
+ std::string name() const override { return "stim"; }
+ NVQIR_SIMULATOR_CLONE_IMPL(StimCircuitSimulator)
+};
+
+} // namespace nvqir
+
+#ifndef __NVQIR_QPP_TOGGLE_CREATE
+/// Register this Simulator with NVQIR.
+NVQIR_REGISTER_SIMULATOR(nvqir::StimCircuitSimulator, stim)
+#endif

runtime/nvqir/stim/stim.yml

@@ -0,0 +1,13 @@
+# ============================================================================ #
+# Copyright (c) 2022 - 2024 NVIDIA Corporation & Affiliates. #
+# All rights reserved. #
+# #
+# This source code and the accompanying materials are made available under #
+# the terms of the Apache License 2.0 which accompanies this distribution. #
+# ============================================================================ #
+
+name: stim
+description: "Stim-based CPU-only backend target"
+config:
+ nvqir-simulation-backend: stim
+ preprocessor-defines: ["-D CUDAQ_SIMULATION_SCALAR_FP64"]

unittests/CMakeLists.txt

@@ -79,6 +79,9 @@ macro (create_tests_with_backend NVQIR_BACKEND EXTRA_BACKEND_TESTER)
  if (${NVQIR_BACKEND} STREQUAL "dm")
  target_compile_definitions(${TEST_EXE_NAME} PRIVATE -DCUDAQ_BACKEND_DM -DCUDAQ_SIMULATION_SCALAR_FP64)
  endif()
+ if (${NVQIR_BACKEND} STREQUAL "stim")
+ target_compile_definitions(${TEST_EXE_NAME} PRIVATE -DCUDAQ_BACKEND_STIM -DCUDAQ_SIMULATION_SCALAR_FP64)
+ endif()
  if (${NVQIR_BACKEND} STREQUAL "tensornet")
  target_compile_definitions(${TEST_EXE_NAME} PRIVATE -DCUDAQ_BACKEND_TENSORNET -DCUDAQ_SIMULATION_SCALAR_FP64)
  set(TEST_LABELS "gpu_required")
@@ -102,6 +105,7 @@ endmacro()
 # We will always have the QPP backend, create a tester for it
 create_tests_with_backend(qpp backends/QPPTester.cpp)
 create_tests_with_backend(dm backends/QPPDMTester.cpp)
+create_tests_with_backend(stim "")
 
 if (CUSTATEVEC_ROOT AND CUDA_FOUND) 
  create_tests_with_backend(custatevec-fp32 "")

unittests/integration/adjoint_tester.cpp

@@ -101,22 +101,27 @@ CUDAQ_TEST(AdjointTester, checkSimple) {
  EXPECT_EQ(1, counts2.size());
  EXPECT_TRUE(counts2.begin()->first == "00000");
 
+#ifndef CUDAQ_BACKEND_STIM
  auto counts3 = cudaq::sample(rotation_adjoint_test{});
  counts3.dump();
  EXPECT_EQ(1, counts3.size());
  EXPECT_TRUE(counts3.begin()->first == "0");
+#endif
 
  auto counts4 = cudaq::sample(twoqbit_adjoint_test{});
  counts4.dump();
  EXPECT_EQ(1, counts4.size());
  EXPECT_TRUE(counts4.begin()->first == "00");
 
+#ifndef CUDAQ_BACKEND_STIM
  auto counts5 = cudaq::sample(test_cudaq_adjoint{});
  counts5.dump();
  EXPECT_EQ(1, counts5.size());
  EXPECT_TRUE(counts5.begin()->first == "101");
+#endif
 }
 
+#ifndef CUDAQ_BACKEND_STIM
 CUDAQ_TEST(AdjointTester, checkNestedAdjoint) {
 
  struct xxxh_gates {
@@ -218,6 +223,7 @@ CUDAQ_TEST(AdjointTester, checkNestedAdjoint) {
  // ctrl ry pi / 4 1 2
  // }
 }
+#endif
 
 #ifndef CUDAQ_BACKEND_DM
 
@@ -251,12 +257,14 @@ static __qpu__ void bar() {
  cudaq::adjoint(foo, q);
 }
 
+#ifndef CUDAQ_BACKEND_STIM
 CUDAQ_TEST(AdjointTester, checkEvenAdjointNesting) {
  auto result = cudaq::get_state(bar);
  std::array<std::complex<double>, 2> expected = {1., 0};
  EXPECT_TRUE(essentially_equal(expected[0], result[0]));
  EXPECT_TRUE(essentially_equal(expected[1], result[1]));
 }
+#endif
 
 static __qpu__ void zaz(cudaq::qubit &q) { rz<cudaq::adj>(M_PI_2, q); }
 
@@ -268,11 +276,13 @@ static __qpu__ void bar_2() {
  cudaq::adjoint(foo_2, q);
 }
 
+#ifndef CUDAQ_BACKEND_STIM
 CUDAQ_TEST(AdjointTester, checkOddAdjointNesting) {
  auto result = cudaq::get_state(bar_2);
  std::array<std::complex<double>, 2> expected = {1., 0};
  EXPECT_TRUE(essentially_equal(expected[0], result[0]));
  EXPECT_TRUE(essentially_equal(expected[1], result[1]));
 }
+#endif
 
 #endif

unittests/integration/async_tester.cpp

@@ -11,6 +11,7 @@
 
 #ifndef CUDAQ_BACKEND_DM
 
+#ifndef CUDAQ_BACKEND_STIM
 CUDAQ_TEST(AsyncTester, checkObserveAsync) {
 
  using namespace cudaq::spin;
@@ -45,6 +46,7 @@ CUDAQ_TEST(AsyncTester, checkObserveAsync) {
  i++;
  }
 }
+#endif
 
 CUDAQ_TEST(AsyncTester, checkSampleAsync) {
  struct ghz {
@@ -71,6 +73,7 @@ CUDAQ_TEST(AsyncTester, checkSampleAsync) {
  cc3.get().dump();
 }
 
+#ifndef CUDAQ_BACKEND_STIM
 CUDAQ_TEST(AsyncTester, checkGetStateAsync) {
  struct ghz {
  auto operator()(int NQubits) __qpu__ {
@@ -107,3 +110,4 @@ CUDAQ_TEST(AsyncTester, checkGetStateAsync) {
  }
 }
 #endif
+#endif

unittests/integration/bug67_vqe_then_sample.cpp

@@ -12,7 +12,8 @@
 #include <cudaq/optimizers.h>
 #include <cudaq/platform.h>
 
-#if !defined CUDAQ_BACKEND_DM && !defined CUDAQ_BACKEND_TENSORNET
+#if !defined(CUDAQ_BACKEND_DM) && !defined(CUDAQ_BACKEND_TENSORNET) && \
+ !defined(CUDAQ_BACKEND_STIM)
 
 CUDAQ_TEST(VqeThenSample, checkBug67) {
 

unittests/integration/bug77_vqe_with_shots.cpp

@@ -12,7 +12,7 @@
 #include <cudaq/optimizers.h>
 #include <cudaq/platform.h>
 
-#ifndef CUDAQ_BACKEND_DM
+#if !defined(CUDAQ_BACKEND_DM) && !defined(CUDAQ_BACKEND_STIM)
 CUDAQ_TEST(VqeWithShots, checkBug77) {
 
  struct ansatz {