Add BitMeasurementRegister to Core.

Add bitIndex to Measure instruction.
Add bit register size to QuantumState.
Add bit measurement register to QuantumState.
Add SparseArray.accumulate.
Add boost/1.85.0 to conanfile.py.

Change OperandsHelper.get_register_operand to work with bit operands.
Change GateConvertor.addGates to accept bit indices in a Measure instruction.
Refactor QuantumState.measure.
Change CMakeLists.txt to work with Boost.
Update CHANGELOG.md with minor aesthetic changes.

Remove QuantumState.{measureAll, prep}.
Remove InstructionExecutor.{MeasureAll, PrepZ}.

CHANGELOG.md

@@ -74,8 +74,8 @@ This project adheres to [Semantic Versioning](http://semver.org/).
 
 ### Changed
 
-- Removed submodules and `deps` folder, replaced with CMake FetchContent.
 - C++23.
+- Manage dependencies via CMake FetchContent.
 - gtest instead of doctest.
 - Various fixes to GitHub actions.
 - Fixed bug in bitset due to bool automatic casting to 32 bits.
@@ -84,7 +84,7 @@ This project adheres to [Semantic Versioning](http://semver.org/).
 
 ### Removed
 
--
+- Submodules and `deps` folder.
 
 ## [ 0.6.2 ] - [ 2023-02-21 ]
 

CMakeLists.txt

@@ -111,6 +111,7 @@ set(CMAKE_WINDOWS_EXPORT_ALL_SYMBOLS ON)
 #=============================================================================#
 
 find_package(absl)
+find_package(Boost REQUIRED)
 find_package(libqasm REQUIRED CONFIG)
 
 
@@ -138,8 +139,9 @@ target_compile_features(qx PRIVATE
     cxx_std_23
 )
 
-target_include_directories(qx
-    PUBLIC "${CMAKE_CURRENT_SOURCE_DIR}/include/"
+target_include_directories(qx PUBLIC
+    "${CMAKE_CURRENT_SOURCE_DIR}/include/"
+    "${Boost_INCLUDE_DIRS}"
 )
 
 if(CMAKE_COMPILER_IS_GNUCXX)
@@ -207,8 +209,12 @@ target_link_libraries(qx PUBLIC
 # Executables                                                                 #
 #=============================================================================#
 
-add_executable("qx-simulator" "${CMAKE_CURRENT_SOURCE_DIR}/src/qx-simulator/Simulator.cpp")
-target_link_libraries(qx-simulator qx)
+add_executable(qx-simulator
+    "${CMAKE_CURRENT_SOURCE_DIR}/src/qx-simulator/Simulator.cpp"
+)
+target_link_libraries(qx-simulator PRIVATE
+    qx
+)
 
 #=============================================================================#
 # Testing                                                                     #
@@ -232,6 +238,17 @@ if(QX_BUILD_PYTHON)
 endif()
 
 
+#=============================================================================#
+# Debug info                                                                  #
+#=============================================================================#
+
+message(STATUS
+    "[${PROJECT_NAME}] Target include directories:\n"
+    "      CMAKE_CURRENT_SOURCE_DIR/include/: ${CMAKE_CURRENT_SOURCE_DIR}/include/\n"
+    "      Boost_INCLUDE_DIRS: ${Boost_INCLUDE_DIRS}\n"
+)
+
+
 #=============================================================================#
 # Installation                                                                #
 #=============================================================================#

conanfile.py

@@ -48,6 +48,7 @@ def build_requirements(self):
 
     def requirements(self):
         self.requires("abseil/20230125.3", transitive_headers=True)
+        self.requires("boost/1.85.0")
         self.requires("fmt/10.2.1", transitive_headers=True)
         self.requires("libqasm/0.6.6", transitive_headers=True)
         self.requires("range-v3/0.12.0", transitive_headers=True)

include/qx/Circuit.hpp

@@ -15,6 +15,7 @@ namespace qx {
 class Circuit {
 public:
     struct Measure {
+        core::BitIndex bitIndex{};
         core::QubitIndex qubitIndex{};
     };
     struct MeasureAll {

include/qx/Core.hpp

@@ -1,6 +1,7 @@
 #pragma once
 
-#include <cstdint>  // size_t
+#include <boost/dynamic_bitset/dynamic_bitset.hpp>
+#include <cstdint>  // size_t, uint32_t
 #include <cstdlib>  // abs
 #include <complex>
 
@@ -19,6 +20,7 @@ struct BitIndex {
 };
 
 using BasisVector = utils::Bitset<config::MAX_QUBIT_NUMBER>;
+using BitMeasurementRegister = boost::dynamic_bitset<uint32_t>;
 
 inline constexpr bool isNotNull(std::complex<double> c) {
 #if defined(_MSC_VER)

include/qx/QuantumState.hpp

@@ -11,6 +11,7 @@
 #include "qx/Core.hpp"  // BasisVector, QubitIndex
 #include "qx/CompileTimeConfiguration.hpp"  // MAX_QUBIT_NUMBER
 #include "qx/DenseUnitaryMatrix.hpp"
+#include "qx/RegisterManager.hpp"
 #include "qx/SparseArray.hpp"
 
 
@@ -25,7 +26,7 @@ template <std::size_t NumberOfOperands>
 void applyImpl(DenseUnitaryMatrix<1 << NumberOfOperands> const &matrix,
                std::array<QubitIndex, NumberOfOperands> const &operands,
                BasisVector index,
-               SparseComplex sparseComplex,
+               SparseComplex const &sparseComplex,
                SparseArray::MapBasisVectorToSparseComplex &storage) {
     utils::Bitset<NumberOfOperands> reducedIndex;
     for (std::size_t i = 0; i < NumberOfOperands; ++i) {
@@ -47,8 +48,9 @@ void applyImpl(DenseUnitaryMatrix<1 << NumberOfOperands> const &matrix,
 
 class QuantumState {
 public:
-    explicit QuantumState(std::size_t n);
+    explicit QuantumState(std::size_t qubit_register_size, std::size_t bit_register_size);
     [[nodiscard]] std::size_t getNumberOfQubits() const;
+    [[nodiscard]] std::size_t getNumberOfBits() const;
     void reset();
     void testInitialize(std::initializer_list<std::pair<std::string, std::complex<double>>> values);
 
@@ -75,92 +77,27 @@ class QuantumState {
     }
 
     [[nodiscard]] const BasisVector &getMeasurementRegister() const;
+    [[nodiscard]] const BitMeasurementRegister &getBitMeasurementRegister() const;
+    [[nodiscard]] double getProbabilityOfMeasuringOne(QubitIndex qubitIndex);
+    [[nodiscard]] double getProbabilityOfMeasuringZero(QubitIndex qubitIndex);
+    void collapseQubit(QubitIndex qubitIndex, bool measuredState, double probabilityOfMeasuringOne);
 
     template <typename F>
-    void measure(QubitIndex qubitIndex, F &&randomGenerator) {
-        auto rand = randomGenerator();
-        double probabilityOfMeasuringOne = 0.;
-        data.forEach([qubitIndex, &probabilityOfMeasuringOne](auto const &kv) {
-            auto const &[basisVector, sparseComplex] = kv;
-            if (basisVector.test(qubitIndex.value)) {
-                probabilityOfMeasuringOne += std::norm(sparseComplex.value);
-            }
-        });
-        if (rand < probabilityOfMeasuringOne) {
-            data.eraseIf([qubitIndex](auto const &kv) {
-                auto const &[basisVector, _] = kv;
-                return !basisVector.test(qubitIndex.value);
-            });
-            data *= std::sqrt(1 / probabilityOfMeasuringOne);
-            measurementRegister.set(qubitIndex.value, true);
-        } else {
-            data.eraseIf([qubitIndex](auto const &kv) {
-                auto const &[basisVector, _] = kv;
-                return basisVector.test(qubitIndex.value);
-            });
-            data *= std::sqrt(1 / (1 - probabilityOfMeasuringOne));
-            measurementRegister.set(qubitIndex.value, false);
-        }
-    }
-
-    template <typename F>
-    void measureAll(F &&randomGenerator) {
-        auto rand = randomGenerator();
-        double probability = 0.;
-        auto [basisVector, sparseComplex] = std::invoke([this, &probability, rand] {
-            for (auto const &[bv, sc] : data) {  // does this work with non-ordered iteration?
-                probability += std::norm(sc.value);
-                if (probability > rand) {
-                    return SparseElement{ bv, sc };
-                }
-            }
-            throw std::runtime_error{ "Vector was not normalized at measurement location (a bug)" };
-        });
-        data.clear();
-        data[basisVector] = SparseComplex{ sparseComplex.value / std::abs(sparseComplex.value) };
-        measurementRegister = basisVector;
+    void measure(BitIndex bitIndex, QubitIndex qubitIndex, F &&randomGenerator) {
+        auto probabilityOfMeasuringOne = getProbabilityOfMeasuringOne(qubitIndex);
+        // measuredState will be true if we measured a 1, or false if we measured a 0
+        auto measuredState = (randomGenerator() < probabilityOfMeasuringOne);
+        collapseQubit(qubitIndex, measuredState, probabilityOfMeasuringOne);
+        measurementRegister.set(qubitIndex.value, measuredState);
+        bitMeasurementRegister.set(bitIndex.value, measuredState);
     }
 
-    template <typename F>
-    void prep(QubitIndex qubitIndex, F &&randomGenerator) {
-        // Measure + conditional X, and reset the measurement register.
-        auto rand = randomGenerator();
-        double probabilityOfMeasuringOne = 0.;
-        data.forEach([qubitIndex, &probabilityOfMeasuringOne](auto const &kv) {
-            auto const &[basisVector, sparseComplex] = kv;
-            if (basisVector.test(qubitIndex.value)) {
-                probabilityOfMeasuringOne += std::norm(sparseComplex.value);
-            }
-        });
-        if (rand < probabilityOfMeasuringOne) {
-            data.eraseIf([qubitIndex](auto const &kv) {
-                auto const &[basisVector, _] = kv;
-                return !basisVector.test(qubitIndex.value);
-            });
-            SparseArray::MapBasisVectorToSparseComplex newData;
-            for (auto [basisVector, sparseComplex] : data) {
-                auto newKey = basisVector;
-                newKey.set(qubitIndex.value, false);
-                newData.insert(std::make_pair(
-                    newKey,
-                    SparseComplex{ sparseComplex.value * std::sqrt(1 / probabilityOfMeasuringOne) }
-                ));
-            }
-            data = newData;  // could fix the interface
-        } else {
-            data.eraseIf([qubitIndex](auto const &kv) {
-                auto const &[basisVector, _] = kv;
-                return basisVector.test(qubitIndex.value);
-            });
-            data *= std::sqrt(1 / (1 - probabilityOfMeasuringOne));
-        }
-        measurementRegister.set(qubitIndex.value, false);
-    };
-
 private:
     std::size_t const numberOfQubits = 1;
+    std::size_t const numberOfBits = 1;
     SparseArray data;
     BasisVector measurementRegister{};
+    BitMeasurementRegister bitMeasurementRegister{};
 };
 
 }  // namespace qx::core

include/qx/SparseArray.hpp

@@ -4,6 +4,7 @@
 #include <algorithm>  // for_each, sort
 #include <complex>
 #include <cstdint>  // size_t, uint64_t
+#include <numeric>  // accumulate
 #include <stdexcept>  // runtime_error
 #include <utility>  // pair
 #include <vector>
@@ -60,6 +61,12 @@ class SparseArray {
     [[nodiscard]] std::size_t size() const;
     [[nodiscard]] std::vector<std::complex<double>> toVector() const;
 
+    template <typename T, typename F>
+    T accumulate(T init, F &&f) {
+        cleanUpZeros();
+        return std::accumulate(data_.begin(), data_.end(), init, f);
+    }
+
     template <typename F>
     void forEach(F &&f) {
         cleanUpZeros();

src/qx/Circuit.cpp

@@ -14,13 +14,7 @@ struct InstructionExecutor {
     explicit InstructionExecutor(core::QuantumState &s) : quantumState(s){};
 
     void operator()(Circuit::Measure const &m) {
-        quantumState.measure(m.qubitIndex, &random::randomZeroOneDouble);
-    }
-    void operator()(Circuit::MeasureAll const &) {
-        quantumState.measureAll(&random::randomZeroOneDouble);
-    }
-    void operator()(Circuit::PrepZ const &r) {
-        quantumState.prep(r.qubitIndex, &random::randomZeroOneDouble);
+        quantumState.measure(m.bitIndex, m.qubitIndex, &random::randomZeroOneDouble);
     }
     void operator()(Circuit::MeasurementRegisterOperation const &op) {
         op.operation(quantumState.getMeasurementRegister());
@@ -78,10 +72,6 @@ void Circuit::execute(core::QuantumState &quantumState, error_models::ErrorModel
         // std::visit(instructionExecutor, instruction);
         if (auto *measure = std::get_if<Circuit::Measure>(&instruction)) {
             instruction_executor(*measure);
-        } else if (auto *measureAll = std::get_if<Circuit::MeasureAll>(&instruction)) {
-            instruction_executor(*measureAll);
-        } else if (auto *prepZ = std::get_if<Circuit::PrepZ>(&instruction)) {
-            instruction_executor(*prepZ);
         } else if (auto *classicalOp = std::get_if<Circuit::MeasurementRegisterOperation>(&instruction)) {
             instruction_executor(*classicalOp);
         } else if (auto *instruction1 = std::get_if<Circuit::Unitary<1>>(&instruction)) {

src/qx/GateConvertor.cpp

@@ -131,10 +131,16 @@ void GateConvertor::addGates(const V3Instruction &instruction) {
         // A measure statement has the following syntax: b = measure q
         // The left-hand side operand, b, is the operand 0
         // The right-hand side operand, q, is the operand 1
-        for (const auto &q: operands_helper.get_register_operand(1)) {
-            auto controlBits = std::make_shared<std::vector<core::QubitIndex>>();
-            circuit_.add_instruction(Circuit::Measure{ core::QubitIndex{ static_cast<std::size_t>(q->value) } },
-                                     controlBits);
+        const auto &bit_indices = operands_helper.get_register_operand(0);
+        const auto &qubit_indices = operands_helper.get_register_operand(1);
+        auto controlBits = std::make_shared<std::vector<core::QubitIndex>>();
+        for (size_t i{ 0 }; i < bit_indices.size(); ++i) {
+            circuit_.add_instruction(
+                Circuit::Measure{
+                    core::BitIndex{ static_cast<std::size_t>(bit_indices[i]->value) },
+                    core::QubitIndex{ static_cast<std::size_t>(qubit_indices[i]->value) }
+                },
+                controlBits);
         }
     } else if (name == "CR") {
         addGates<2>(

src/qx/OperandsHelper.cpp

@@ -15,25 +15,38 @@ OperandsHelper::OperandsHelper(const V3Instruction &instruction, const RegisterM
 
 [[nodiscard]] V3Many<V3ConstInt> OperandsHelper::get_register_operand(int id) const {
     if (auto variable_ref = instruction_.operands[id]->as_variable_ref()) {
-        if (!is_qubit_variable(*variable_ref->variable)) {
+        auto ret = V3Many<V3ConstInt>{};
+        if (is_qubit_variable(*variable_ref->variable)) {
+            auto qubit_range = register_manager_.get_qubit_range(variable_ref->variable->name);
+            ret.get_vec().resize(qubit_range.size);
+            std::generate_n(ret.get_vec().begin(), qubit_range.size, [qubit_range, i=0]() mutable {
+                return v3_tree::make<V3ConstInt>(static_cast<v3_primitives::Int>(qubit_range.first + i++));
+            });
+        } else if (is_bit_variable(*variable_ref->variable)) {
+            auto bit_range = register_manager_.get_bit_range(variable_ref->variable->name);
+            ret.get_vec().resize(bit_range.size);
+            std::generate_n(ret.get_vec().begin(), bit_range.size, [bit_range, i=0]() mutable {
+                return v3_tree::make<V3ConstInt>(static_cast<v3_primitives::Int>(bit_range.first + i++));
+            });
+       } else {
             return {};
         }
-        auto qubit_range = register_manager_.get_qubit_range(variable_ref->variable->name);
-        auto ret = V3Many<V3ConstInt>{};
-        ret.get_vec().resize(qubit_range.size);
-        std::generate_n(ret.get_vec().begin(), qubit_range.size, [qubit_range, i=0]() mutable {
-            return v3_tree::make<V3ConstInt>(static_cast<v3_primitives::Int>(qubit_range.first + i++));
-        });
         return ret;
     } else if (auto index_ref = instruction_.operands[id]->as_index_ref()) {
-        if (!is_qubit_variable(*index_ref->variable)) {
+        auto ret = index_ref->indices;
+        if (is_qubit_variable(*index_ref->variable)) {
+            auto qubit_range = register_manager_.get_qubit_range(index_ref->variable->name);
+            std::for_each(ret.get_vec().begin(), ret.get_vec().end(), [qubit_range](const auto &index) {
+                index->value += qubit_range.first;
+            });
+        } else if (is_bit_variable(*index_ref->variable)) {
+            auto bit_range = register_manager_.get_bit_range(index_ref->variable->name);
+            std::for_each(ret.get_vec().begin(), ret.get_vec().end(), [bit_range](const auto &index) {
+                index->value += bit_range.first;
+            });
+        } else {
             return {};
         }
-        auto qubit_range = register_manager_.get_qubit_range(index_ref->variable->name);
-        auto ret = index_ref->indices;
-        std::for_each(ret.get_vec().begin(), ret.get_vec().end(), [qubit_range](const auto &index) {
-            index->value += qubit_range.first;
-        });
         return ret;
     }
     assert(false && "operand is neither a variable reference nor an index reference");