feat: add clifford simulator(#1)

.github/workflows/wheels.yml

@@ -27,9 +27,9 @@ jobs:
                 -   os-arch: win_amd64
                     os: windows-2019
                 -   os-arch: macosx_x86_64
-                    os: macos-11
+                    os: macos-13
                 -   os-arch: macosx_arm64
-                    os: macos-11
+                    os: macos-13
         runs-on: ${{ matrix.os }}
 
         env:

CMakeLists.txt

@@ -9,8 +9,9 @@ if(NOT CMAKE_BUILD_TYPE)
     set(CMAKE_BUILD_TYPE Release)
 endif()
 
-set(CMAKE_CXX_STANDARD 17)
+set(CMAKE_CXX_STANDARD 20)
 set(CMAKE_CUDA_ARCHITECTURES 70;75;80;90)
+set(CMAKE_EXPORT_COMPILE_COMMANDS ON)
 if(SKBUILD)
 
 execute_process(
@@ -46,7 +47,7 @@ include(ExternalProject)
 ExternalProject_Add(Eigen3
 		PREFIX             ${EIGEN3_ROOT}
 		GIT_REPOSITORY     https://gitlab.com/libeigen/eigen.git
-		GIT_TAG            3.3.9
+		GIT_TAG            3.4
 
 		CONFIGURE_COMMAND ""
 		BUILD_COMMAND ""
@@ -56,7 +57,6 @@ ExternalProject_Add(Eigen3
 )
 list (APPEND PRJ_INCLUDE_DIRS ${EIGEN3_INCLUDE_DIR})
 
-
 find_package(pybind11 CONFIG)
 list (APPEND PRJ_INCLUDE_DIRS ${PYBIND11_INCLUDE_DIR})
 
@@ -74,7 +74,7 @@ if(CMAKE_HOST_SYSTEM_PROCESSOR STREQUAL "x86_64" OR CMAKE_HOST_SYSTEM_PROCESSOR
 	endif()
 endif()
 
-list (APPEND PRJ_INCLUDE_DIRS src/qfvm)
+list (APPEND PRJ_INCLUDE_DIRS src/qfvm src/qfvm_clifford)
 pybind11_add_module(${PROJECT_NAME} MODULE src/${PROJECT_NAME}/${PROJECT_NAME}.cpp)
 add_dependencies(${PROJECT_NAME} Eigen3) #must add dependence for ninja
 target_compile_options(${PROJECT_NAME} PUBLIC ${PRJ_COMPILE_OPTIONS})

pyproject.toml

@@ -18,6 +18,7 @@ archs = ["x86_64"]
 
 
 [tool.cibuildwheel.macos]
+environment = {MACOSX_DEPLOYMENT_TARGET = "13.6"}
 archs = ["x86_64", "arm64"]
 
 repair-wheel-command = [

quafu/results/results.py

@@ -91,14 +91,21 @@ class SimuResult(Result):
     """
 
     def __init__(self, input, input_form, count_dict: dict = None):
-        self.num = int(np.log2(input.shape[0]))
+        if input_form != "count_dict":
+            self.num = int(np.log2(input.shape[0]))
+        else:
+            # input is num qubits
+            self.num = input
         if input_form == "density_matrix":
             self.rho = np.array(input)
             self.probabilities = np.diag(input)
         elif input_form == "probabilities":
             self.probabilities = input
         elif input_form == "state_vector":
             self.state_vector = input
+        elif input_form == "count_dict":
+            # do nothing, only count dict
+            pass
         # come form c++ simulator
         # TODO: add count for py_simu
         if count_dict is not None:

quafu/simulators/simulator.py

@@ -86,6 +86,7 @@ def simulate(
         if use_gpu:
             if qc.executable_on_backend == False:
                 raise QuafuError("classical operation only support for `qfvm_qasm`")
+
             if use_custatevec:
                 try:
                     from .qfvm import simulate_circuit_custate
@@ -101,6 +102,14 @@ def simulate(
         else:
             count_dict, psi = simulate_circuit(qc, psi, shots)
 
+    elif simulator == "qfvm_clifford":
+        try:
+            from .qfvm import simulate_circuit_clifford
+        except ImportError:
+            raise QuafuError("you are not using the clifford version of pyquafu")
+
+        count_dict = simulate_circuit_clifford(qc, shots)
+
     elif simulator == "py_simu":
         if qc.executable_on_backend == False:
             raise QuafuError("classical operation only support for `qfvm_qasm`")
@@ -129,6 +138,9 @@ def simulate(
     elif output == "state_vector":
         return SimuResult(psi, output, count_dict)
 
+    elif output == "count_dict":
+        return SimuResult(max(qc.used_qubits) + 1, output, count_dict)
+
     else:
         raise ValueError(
             "output should in be 'density_matrix', 'probabilities', or 'state_vector'"

src/qfvm/qfvm.cpp

@@ -3,6 +3,8 @@
 #include <pybind11/numpy.h>
 #include <pybind11/pybind11.h>
 #include <random>
+#include <tuple>
+
 #ifdef _USE_GPU
 #include <cuda_simulator.cuh>
 #endif
@@ -11,6 +13,12 @@
 #include <custate_simu.cuh>
 #endif
 
+#ifdef USE_SIMD
+constexpr size_t _word_size = 256;
+#else
+constexpr size_t _word_size = 64;
+#endif
+
 namespace py = pybind11;
 
 template <typename T>
@@ -115,6 +123,82 @@ simulate_circuit(py::object const& pycircuit,
     return std::make_pair(outcount, np_inputstate);
 }
 
+std::map<uint, uint> simulate_circuit_clifford(py::object const& pycircuit,
+                                               const int& shots) {
+
+  auto circuit = Circuit(pycircuit);
+
+  // If measure all at the end, simulate once
+  uint actual_shots = shots;
+
+  // qbit, cbit
+  vector<std::pair<uint, uint>> measures = circuit.measure_vec();
+  std::map<uint, bool> cbit_measured;
+  for (auto& pair : measures) {
+    cbit_measured[pair.second] = true;
+  }
+
+  // Store outcome's count
+  std::map<uint, uint> outcount;
+
+  circuit_simulator<_word_size> cs(circuit.qubit_num());
+
+  for (uint i = 0; i < actual_shots; i++) {
+
+    simulate(circuit, cs);
+    uint outcome = 0;
+
+    if (!circuit.final_measure()) {
+      // qubit, cbit, measure result
+      auto measure_results = cs.current_measurement_record();
+
+      // make sure the order is the same with other simulators
+      std::sort(
+          measure_results.begin(), measure_results.end(),
+          [](auto& a, auto& b) { return std::get<1>(a) < std::get<1>(b); });
+
+      for (auto& measure_result : measure_results) {
+        outcome *= 2;
+        outcome += std::get<2>(measure_result);
+      }
+
+    } else if (circuit.final_measure() && !measures.empty()) {
+      for (auto& measure : measures) {
+        cs.do_circuit_instruction(
+            {"measure", std::vector<size_t>{measure.first},
+             std::vector<double>{static_cast<double>(measure.second)}});
+      }
+
+      // qubit, cbit, measure result
+      auto measure_results = cs.current_measurement_record();
+
+      // make sure the order is the same with other simulators
+      std::sort(
+          measure_results.begin(), measure_results.end(),
+          [](auto& a, auto& b) { return std::get<1>(a) < std::get<1>(b); });
+
+      for (auto& measure_result : measure_results) {
+        outcome *= 2;
+        outcome += std::get<2>(measure_result);
+      }
+    }
+
+    if (measures.empty()) {
+      continue;
+    }
+
+    if (outcount.find(outcome) != outcount.end())
+      outcount[outcome]++;
+    else
+      outcount[outcome] = 1;
+
+    cs.reset_tableau();
+    cs.sim_record.clear();
+  }
+
+  return outcount;
+}
+
 #ifdef _USE_GPU
 py::object simulate_circuit_gpu(py::object const& pycircuit,
                                 py::array_t<complex<double>>& np_inputstate) {
@@ -164,6 +248,9 @@ PYBIND11_MODULE(qfvm, m) {
         py::arg("circuit"),
         py::arg("inputstate") = py::array_t<complex<double>>(0),
         py::arg("shots"));
+  m.def("simulate_circuit_clifford", &simulate_circuit_clifford,
+        "Simulate with circuit using clifford", py::arg("circuit"),
+        py::arg("shots"));
 
 #ifdef _USE_GPU
   m.def("simulate_circuit_gpu", &simulate_circuit_gpu, "Simulate with circuit",

src/qfvm/simulator.hpp

@@ -1,7 +1,12 @@
 #pragma once
 
 #include "circuit.hpp"
+#include "clifford_simulator.h"
+#include "qasm.hpp"
 #include "statevector.hpp"
+#include "types.hpp"
+#include <cstddef>
+#include <vector>
 
 void apply_op(QuantumOperator& op, StateVector<data_t>& state) {
   bool matched = false;
@@ -114,3 +119,44 @@ void simulate(Circuit const& circuit, StateVector<data_t>& state) {
     apply_op(op, state);
   }
 }
+
+template <size_t word_size>
+void apply_measure(circuit_simulator<word_size>& cs, const vector<pos_t>& qbits,
+                   const vector<pos_t>& cbits) {
+  for (size_t i = 0; i < qbits.size(); i++) {
+    cs.do_circuit_instruction(
+        {"measure", std::vector<size_t>{qbits[i]},
+         std::vector<double>{static_cast<double>(cbits[i])}});
+  }
+}
+
+template <size_t word_size>
+void apply_op(QuantumOperator& op, circuit_simulator<word_size>& cs) {
+  // TODO: support args
+  switch (OPMAP[op.name()]) {
+  case Opname::measure:
+    apply_measure(cs, op.qbits(), op.cbits());
+    break;
+  case Opname::reset:
+    for (auto qubit : op.qbits()) {
+      cs.do_circuit_instruction(
+          {"reset", std::vector<size_t>{static_cast<size_t>(qubit)}});
+    }
+    break;
+  default:
+    auto qubits = op.positions();
+    cs.do_circuit_instruction(
+        {op.name(), std::vector<size_t>(qubits.begin(), qubits.end())});
+  }
+}
+
+template <size_t word_size>
+void simulate(Circuit const& circuit, circuit_simulator<word_size>& cs) {
+  // skip measure and handle it in qfvm.cpp
+  bool skip_measure = circuit.final_measure();
+  for (auto op : circuit.instructions()) {
+    if (skip_measure == true && op.name() == "measure")
+      continue;
+    apply_op(op, cs);
+  }
+}

src/qfvm_clifford/bit.h

@@ -0,0 +1,55 @@
+#ifndef BIT_H_
+#define BIT_H_
+
+#include <cstddef>
+#include <cstdint>
+
+// bit in byte
+struct bit {
+  uint8_t* byte;
+  uint8_t byte_index;
+
+  bit(void* ptr, size_t offset)
+      : byte(((uint8_t*)ptr + (offset / 8))), byte_index(offset & 7) {}
+
+  // copy assignment for bit in byte
+  inline bit& operator=(bool value) {
+    // make bit be 0
+    *byte &= ~((uint8_t)1 << byte_index);
+    // assignment
+    *byte |= uint8_t(value) << byte_index;
+    return *this;
+  }
+
+  inline bit& operator=(const bit& other) {
+    *this = bool(other);
+    return *this;
+  }
+
+  // bit operator
+  inline bit& operator^=(bool value) {
+    *byte ^= uint8_t(value) << byte_index;
+    return *this;
+  }
+
+  inline bit& operator&=(bool value) {
+    *byte &= (uint8_t(value) << byte_index) | ~(uint8_t(1) << byte_index);
+    return *this;
+  }
+
+  inline bit& operator|=(bool value) {
+    *byte |= uint8_t(value) << byte_index;
+    return *this;
+  }
+
+  // conversion operator
+  inline operator bool() const { return (*byte >> byte_index) & 1; }
+
+  void swap(bit other) {
+    bool b = bool(other);
+    other = bool(*this);
+    *this = b;
+  }
+};
+
+#endif