Add fix and test for emulation of controlled time evolution.

projectq/backends/_sim/_cppkernels/simulator.hpp

@@ -335,14 +335,15 @@ class Simulator{
         unsigned s = std::abs(time) * op_nrm + 1.;
         complex_type correction = std::exp(-time * I * tr / (double)s);
         auto output_state = vec_;
+        auto ctrlmask = get_control_mask(ctrl);
         for (unsigned i = 0; i < s; ++i){
             calc_type nrm_change = 1.;
             for (unsigned k = 0; nrm_change > 1.e-12; ++k){
                 auto coeff = (-time * I) / double(s * (k + 1));
                 auto current_state = vec_;
                 auto update = StateVector(vec_.size(), 0.);
                 for (auto const& tup : td){
-                    apply_term(tup.first, ids, ctrl);
+                    apply_term(tup.first, ids, {});
                     #pragma omp parallel for schedule(static)
                     for (std::size_t j = 0; j < vec_.size(); ++j){
                         update[j] += vec_[j] * tup.second;
@@ -354,14 +355,17 @@ class Simulator{
                 for (std::size_t j = 0; j < vec_.size(); ++j){
                     update[j] *= coeff;
                     vec_[j] = update[j];
-                    output_state[j] += update[j];
-                    nrm_change += std::norm(update[j]);
+                    if ((j & ctrlmask) == ctrlmask){
+                        output_state[j] += update[j];
+                        nrm_change += std::norm(update[j]);
+                    }
                 }
                 nrm_change = std::sqrt(nrm_change);
             }
             #pragma omp parallel for schedule(static)
             for (std::size_t j = 0; j < vec_.size(); ++j){
-                output_state[j] *= correction;
+                if ((j & ctrlmask) == ctrlmask)
+                    output_state[j] *= correction;
                 vec_[j] = output_state[j];
             }
         }

projectq/backends/_sim/_pysim.py

@@ -351,6 +351,7 @@ def emulate_time_evolution(self, terms_dict, time, ids, ctrlids):
         s = int(op_nrm + 1.)
         correction = _np.exp(-1j * time * tr / float(s))
         output_state = _np.copy(self._state)
+        mask = self._get_control_mask(ctrlids)
         for i in range(s):
             j = 0
             nrm_change = 1.
@@ -359,16 +360,20 @@ def emulate_time_evolution(self, terms_dict, time, ids, ctrlids):
                 current_state = _np.copy(self._state)
                 update = 0j
                 for t, c in terms_dict:
-                    self._apply_term(t, ids, ctrlids)
+                    self._apply_term(t, ids)
                     self._state *= c
                     update += self._state
                     self._state = _np.copy(current_state)
                 update *= coeff
                 self._state = update
-                output_state += update
+                for i in range(len(update)):
+                    if (i & mask) == mask:
+                        output_state[i] += update[i]
                 nrm_change = _np.linalg.norm(update)
                 j += 1
-            output_state *= correction
+            for i in range(len(update)):
+                if (i & mask) == mask:
+                    output_state[i] *= correction
             self._state = _np.copy(output_state)
 
     def apply_controlled_gate(self, m, ids, ctrlids):

projectq/backends/_sim/_simulator_test.py

@@ -419,7 +419,7 @@ def test_simulator_applyqubitoperator(sim):
 
 
 def test_simulator_time_evolution(sim):
-    N = 9  # number of qubits
+    N = 8  # number of qubits
     time_to_evolve = 1.1  # time to evolve for
     eng = MainEngine(sim, [])
     qureg = eng.allocate_qureg(N)
@@ -435,10 +435,13 @@ def test_simulator_time_evolution(sim):
     op += 1.1 * Qop(())
     op += -1.4 * Qop("Y0 Z1 X3 Y5")
     op += -1.1 * Qop("Y1 X2 X3 Y4")
-    TimeEvolution(time_to_evolve, op) | qureg
+    ctrl_qubit = eng.allocate_qubit()
+    H | ctrl_qubit
+    with Control(eng, ctrl_qubit):
+        TimeEvolution(time_to_evolve, op) | qureg
     eng.flush()
     qbit_to_bit_map, final_wavefunction = copy.deepcopy(eng.backend.cheat())
-    Measure | qureg
+    Measure | (qureg, ctrl_qubit)
     # Check manually:
 
     def build_matrix(list_single_matrices):
@@ -465,7 +468,13 @@ def build_matrix(list_single_matrices):
     init_wavefunction = numpy.array(init_wavefunction, copy=False)
     final_wavefunction = numpy.array(final_wavefunction, copy=False)
     res = scipy.sparse.linalg.expm_multiply(res_matrix, init_wavefunction)
-    assert numpy.allclose(res, final_wavefunction)
+
+    half = int(len(final_wavefunction) / 2)
+    hadamard_f = 1. / math.sqrt(2.)
+    # check evolution and control
+    assert numpy.allclose(hadamard_f * res, final_wavefunction[half:])
+    assert numpy.allclose(final_wavefunction[:half], hadamard_f *
+                          init_wavefunction)
 
 
 def test_simulator_set_wavefunction(sim):