Add measurement gate. Fixes.

apps/qsim_amplitudes.cc

@@ -172,7 +172,7 @@ int main(int argc, char* argv[]) {
     }
   };
 
-  using Runner = QSimRunner<IO, BasicGateFuser<GateQSim<float>>, Simulator>;
+  using Runner = QSimRunner<IO, BasicGateFuser<IO, GateQSim<float>>, Simulator>;
 
   Runner::Parameter param;
   param.seed = opt.seed;

apps/qsim_base.cc

@@ -112,7 +112,7 @@ int main(int argc, char* argv[]) {
   using Simulator = qsim::Simulator<For>;
   using StateSpace = Simulator::StateSpace;
   using State = StateSpace::State;
-  using Runner = QSimRunner<IO, BasicGateFuser<GateQSim<float>>, Simulator>;
+  using Runner = QSimRunner<IO, BasicGateFuser<IO, GateQSim<float>>, Simulator>;
 
   StateSpace state_space(circuit.num_qubits, opt.num_threads);
   State state = state_space.CreateState();

apps/qsim_von_neumann.cc

@@ -113,7 +113,7 @@ int main(int argc, char* argv[]) {
     IO::messagef("entropy=%g\n", entropy);
   };
 
-  using Runner = QSimRunner<IO, BasicGateFuser<GateQSim<float>>, Simulator>;
+  using Runner = QSimRunner<IO, BasicGateFuser<IO, GateQSim<float>>, Simulator>;
 
   Runner::Parameter param;
   param.seed = opt.seed;

lib/fuser_basic.h

@@ -24,7 +24,7 @@
 
 namespace qsim {
 
-template <typename Gate>
+template <typename IO, typename Gate>
 struct BasicGateFuser final {
   using GateFused = qsim::GateFused<Gate>;
 
@@ -103,8 +103,8 @@ struct BasicGateFuser final {
 
         if (gate.time < prev_time) {
           // This function assumes that gate times are ordered.
-          // Just stop sielently if this is not the case.
-          // TODO: report an error here.
+          // Just stop silently if this is not the case.
+          IO::errorf("gate times should be ordered.\n");
           gates_fused.resize(0);
           return gates_fused;
         }

lib/gate_appl.h

@@ -62,7 +62,8 @@ inline void CalcMatrix4(unsigned q0, unsigned q1,
 }
 
 /**
- * Applies the given gate to the simulator state.
+ * Applies the given gate to the simulator state. Measurement gates should not
+ * be applied in this function.
  * @param simulator Simulator object. Provides specific implementations for
  *   applying one- and two-qubit gates.
  * @param gate The gate to be applied.
@@ -85,7 +86,8 @@ inline void ApplyGate(
 }
 
 /**
- * Applies the given fused gate to the simulator state.
+ * Applies the given fused gate to the simulator state. Measurement gates
+ * should not be applied in this function.
  * @param simulator Simulator object. Provides specific implementations for
  *   applying one- and two-qubit gates.
  * @param gate The gate to be applied.

lib/hybrid.h

@@ -26,7 +26,8 @@
 namespace qsim {
 
 // Hybrid Feynmann-Schrodiner simulator.
-template <typename IO, typename GateT, template <typename> class FuserT,
+template <typename IO, typename GateT,
+          template <typename, typename> class FuserT,
           typename Simulator, typename For>
 struct HybridSimulator final {
  public:
@@ -64,7 +65,7 @@ struct HybridSimulator final {
   };
 
  public:
-  using Fuser = FuserT<GateHybrid>;
+  using Fuser = FuserT<IO, GateHybrid>;
   using GateFused = typename Fuser::GateFused;
 
   struct HybridData {

lib/parfor.h

@@ -18,16 +18,13 @@
 #include <omp.h>
 
 #include <cstdint>
+#include <utility>
 #include <vector>
 
 namespace qsim {
 
 template <uint64_t MIN_SIZE>
 struct ParallelForT {
-  static bool CanBeParallel(uint64_t size) {
-    return size >= MIN_SIZE;
-  }
-
   static uint64_t GetIndex0(
       uint64_t size, unsigned num_threads, unsigned thread_id) {
     return size >= MIN_SIZE ? size * thread_id / num_threads : 0;
@@ -63,7 +60,7 @@ struct ParallelForT {
 
   template <typename Function, typename Op, typename... Args>
   static std::vector<typename Op::result_type> RunReduceP(
-      unsigned num_threads, uint64_t size, Function&& func, const Op& op,
+      unsigned num_threads, uint64_t size, Function&& func, Op&& op,
       Args&&... args) {
     std::vector<typename Op::result_type> partial_results;
 
@@ -101,8 +98,9 @@ struct ParallelForT {
   template <typename Function, typename Op, typename... Args>
   static typename Op::result_type RunReduce(unsigned num_threads,
                                             uint64_t size, Function&& func,
-                                            const Op& op, Args&&... args) {
-    auto partial_results = RunReduceP(num_threads, size, func, op, args...);
+                                            Op&& op, Args&&... args) {
+    auto partial_results = RunReduceP(
+        num_threads, size, func, std::move(op), args...);
 
     typename Op::result_type result = 0;
 

lib/run_qsim.h

@@ -86,6 +86,9 @@ struct QSimRunner final {
 
     auto fused_gates = Fuser::FuseGates(circuit.num_qubits, circuit.gates,
                                         times_to_measure_at);
+    if (fused_gates.size() == 0 && circuit.gates.size() > 0) {
+      return false;
+    }
 
     unsigned cur_time_index = 0;
 
@@ -148,6 +151,9 @@ struct QSimRunner final {
 
     auto fused_gates = Fuser::FuseGates(circuit.num_qubits, circuit.gates,
                                         maxtime);
+    if (fused_gates.size() == 0 && circuit.gates.size() > 0) {
+      return false;
+    }
 
     // Apply fused gates.
     for (std::size_t i = 0; i < fused_gates.size(); ++i) {

lib/run_qsimh.h

@@ -82,7 +82,14 @@ struct QSimHRunner final {
     }
 
     auto fgates0 = Fuser::FuseGates(hd.num_qubits0, hd.gates0, maxtime);
+    if (fgates0.size() == 0 && hd.gates0.size() > 0) {
+      return false;
+    }
+
     auto fgates1 = Fuser::FuseGates(hd.num_qubits1, hd.gates1, maxtime);
+    if (fgates1.size() == 0 && hd.gates1.size() > 0) {
+      return false;
+    }
 
     rc = HybridSimulator::Run(
         param, hd, parts, fgates0, fgates1, bitstrings, results);

lib/seqfor.h

@@ -16,15 +16,12 @@
 #define SEQFOR_H_
 
 #include <cstdint>
+#include <utility>
 #include <vector>
 
 namespace qsim {
 
 struct SequentialFor {
-  static bool CanBeParallel(uint64_t size) {
-    return false;
-  }
-
   static uint64_t GetIndex0(
       uint64_t size, unsigned num_threads, unsigned thread_id) {
     return 0;
@@ -45,7 +42,7 @@ struct SequentialFor {
 
   template <typename Function, typename Op, typename... Args>
   static std::vector<typename Op::result_type> RunReduceP(
-      unsigned num_threads, uint64_t size, Function&& func, const Op& op,
+      unsigned num_threads, uint64_t size, Function&& func, Op&& op,
       Args&&... args) {
     typename Op::result_type result = 0;
 
@@ -59,8 +56,8 @@ struct SequentialFor {
   template <typename Function, typename Op, typename... Args>
   static typename Op::result_type RunReduce(unsigned num_threads,
                                             uint64_t size, Function&& func,
-                                            const Op& op, Args&&... args) {
-    return RunReduceP(num_threads, size, func, op, args...)[0];
+                                            Op&& op, Args&&... args) {
+    return RunReduceP(num_threads, size, func, std::move(op), args...)[0];
   }
 };
 

lib/statespace.h

@@ -35,7 +35,7 @@ inline void do_not_free(void*) noexcept {}
 }  // namespace detail
 
 // Routines for state-vector manipulations.
-template <typename For, typename FP>
+template <typename Impl, typename For, typename FP>
 class StateSpace {
  public:
   using fp_type = FP;
@@ -103,10 +103,8 @@ class StateSpace {
     For::Run(num_threads_, raw_size_, f, src, dest);
   }
 
- protected:
-  template <typename Func>
-  double Norm(uint64_t size, Func&& PartialNorms, const State& state) const {
-    auto partial_norms = PartialNorms(num_threads_, size, state);
+  double Norm(const State& state) const {
+    auto partial_norms = static_cast<const Impl&>(*this).PartialNorms(state);
 
     double norm = partial_norms[0];
     for (std::size_t i = 1; i < partial_norms.size(); ++i) {
@@ -116,10 +114,22 @@ class StateSpace {
     return norm;
   }
 
-  template <typename RGen, typename Func1, typename Func2>
-  MeasurementResult VirtualMeasure(
-      const std::vector<unsigned>& qubits, RGen& rgen, const State& state,
-      uint64_t size, Func1&& PartialNorms, Func2&& FindMeasuredBits) const {
+  template <typename RGen>
+  MeasurementResult Measure(const std::vector<unsigned>& qubits,
+                            RGen& rgen, State& state) const {
+    auto result =
+        static_cast<const Impl&>(*this).VirtualMeasure(qubits, rgen, state);
+
+    if (result.valid) {
+      static_cast<const Impl&>(*this).CollapseState(result, state);
+    }
+
+    return result;
+  }
+
+  template <typename RGen>
+  MeasurementResult VirtualMeasure(const std::vector<unsigned>& qubits,
+                                   RGen& rgen, const State& state) const {
     MeasurementResult result;
 
     result.valid = true;
@@ -134,7 +144,7 @@ class StateSpace {
       result.mask |= uint64_t{1} << q;
     }
 
-    auto partial_norms = PartialNorms(num_threads_, size, state);
+    auto partial_norms = static_cast<const Impl&>(*this).PartialNorms(state);
 
     for (std::size_t i = 1; i < partial_norms.size(); ++i) {
       partial_norms[i] += partial_norms[i - 1];
@@ -145,12 +155,12 @@ class StateSpace {
 
     unsigned m = 0;
     while (r > partial_norms[m]) ++m;
-    r -= m > 0 ? partial_norms[m - 1] : 0;
-
-    uint64_t k0 = For::GetIndex0(size, num_threads_, m);
-    uint64_t k1 = For::GetIndex1(size, num_threads_, m);
+    if (m > 0) {
+      r -= partial_norms[m - 1];
+    }
 
-    result.bits = FindMeasuredBits(k0, k1, r, result.mask, state);
+    result.bits = static_cast<const Impl&>(*this).FindMeasuredBits(
+        m, r, result.mask, state);
 
     result.bitstring.reserve(qubits.size());
     result.bitstring.resize(0);