Added probability and amplitude functions to simulator. (#111)

projectq/backends/_sim/_cppkernels/simulator.hpp

@@ -264,6 +264,41 @@ class Simulator{
         return expectation;
     }
 
+    calc_type get_probability(std::vector<bool> const& bit_string,
+                              std::vector<unsigned> const& ids){
+        run();
+        for (auto id : ids)
+            if (map_.count(id) == 0)
+                throw(std::runtime_error("get_probability(): Unknown qubit id. Please make sure you have called eng.flush()."));
+        std::size_t mask = 0, bit_str = 0;
+        for (unsigned i = 0; i < ids.size(); ++i){
+            mask |= 1UL << map_[ids[i]];
+            bit_str |= (bit_string[i]?1UL:0UL) << map_[ids[i]];
+        }
+        calc_type probability = 0.;
+        #pragma omp parallel for reduction(+:probability) schedule(static)
+        for (std::size_t i = 0; i < vec_.size(); ++i)
+            if ((i & mask) == bit_str)
+                probability += std::norm(vec_[i]);
+        return probability;
+    }
+
+    complex_type const& get_amplitude(std::vector<bool> const& bit_string,
+                                      std::vector<unsigned> const& ids){
+        run();
+        std::size_t chk = 0;
+        std::size_t index = 0;
+        for (unsigned i = 0; i < ids.size(); ++i){
+            if (map_.count(ids[i]) == 0)
+                break;
+            chk |= 1UL << map_[ids[i]];
+            index |= (bit_string[i]?1UL:0UL) << map_[ids[i]];
+        }
+        if (chk + 1 != vec_.size())
+            throw(std::runtime_error("The second argument to get_amplitude() must be a permutation of all allocated qubits. Please make sure you have called eng.flush()."));
+        return vec_[index];
+    }
+
     void emulate_time_evolution(TermsDict const& tdict, calc_type const& time,
                                 std::vector<unsigned> const& ids,
                                 std::vector<unsigned> const& ctrl){

projectq/backends/_sim/_cppsim.cpp

@@ -52,6 +52,8 @@ PYBIND11_PLUGIN(_cppsim) {
         .def("emulate_math", &emulate_math_wrapper<QuRegs>)
         .def("get_expectation_value", &Simulator::get_expectation_value)
         .def("emulate_time_evolution", &Simulator::emulate_time_evolution)
+        .def("get_probability", &Simulator::get_probability)
+        .def("get_amplitude", &Simulator::get_amplitude)
         .def("run", &Simulator::run)
         .def("cheat", &Simulator::cheat)
         ;

projectq/backends/_sim/_pysim.py

@@ -248,6 +248,65 @@ def get_expectation_value(self, terms_dict, ids):
             self._state = _np.copy(current_state)
         return expectation
 
+    def get_probability(self, bit_string, ids):
+        """
+        Return the probability of the outcome `bit_string` when measuring
+        the qubits given by the list of ids.
+
+        Args:
+            bit_string (list[bool|int]): Measurement outcome.
+            ids (list[int]): List of qubit ids determining the ordering.
+
+        Returns:
+            Probability of measuring the provided bit string.
+
+        Raises:
+            RuntimeError if an unknown qubit id was provided.
+        """
+        for i in range(len(ids)):
+            if ids[i] not in self._map:
+                raise RuntimeError("get_probability(): Unknown qubit id. "
+                                   "Please make sure you have called "
+                                   "eng.flush().")
+        mask = 0
+        bit_str = 0
+        for i in range(len(ids)):
+            mask |= (1 << self._map[ids[i]])
+            bit_str |= (bit_string[i] << self._map[ids[i]])
+        probability = 0.
+        for i in range(len(self._state)):
+            if (i & mask) == bit_str:
+                e = self._state[i]
+                probability += e.real**2 + e.imag**2
+        return probability
+
+    def get_amplitude(self, bit_string, ids):
+        """
+        Return the probability amplitude of the supplied `bit_string`.
+        The ordering is given by the list of qubit ids.
+
+        Args:
+            bit_string (list[bool|int]): Computational basis state
+            ids (list[int]): List of qubit ids determining the
+                ordering. Must contain all allocated qubits.
+
+        Returns:
+            Probability amplitude of the provided bit string.
+
+        Raises:
+            RuntimeError if the second argument is not a permutation of all
+            allocated qubits.
+        """
+        if not set(ids) == set(self._map):
+            raise RuntimeError("The second argument to get_amplitude() must"
+                               " be a permutation of all allocated qubits. "
+                               "Please make sure you have called "
+                               "eng.flush().")
+        index = 0
+        for i in range(len(ids)):
+            index |= (bit_string[i] << self._map[ids[i]])
+        return self._state[index]
+
     def emulate_time_evolution(self, terms_dict, time, ids, ctrlids):
         """
         Applies exp(-i*time*H) to the wave function, i.e., evolves under

projectq/backends/_sim/_simulator.py

@@ -126,6 +126,40 @@ def get_expectation_value(self, qubit_operator, qureg):
         return self._simulator.get_expectation_value(operator,
                                                      [qb.id for qb in qureg])
 
+    def get_probability(self, bit_string, qureg):
+        """
+        Return the probability of the outcome `bit_string` when measuring
+        the quantum register `qureg`.
+
+        Args:
+            bit_string (list[bool|int]|string[0|1]): Measurement outcome.
+            qureg (Qureg|list[Qubit]): Quantum register.
+
+        Returns:
+            Probability of measuring the provided bit string.
+        """
+        bit_string = [bool(int(b)) for b in bit_string]
+        return self._simulator.get_probability(bit_string,
+                                               [qb.id for qb in qureg])
+
+    def get_amplitude(self, bit_string, qureg):
+        """
+        Return the probability amplitude of the supplied `bit_string`.
+        The ordering is given by the quantum register `qureg`, which must
+        contain all allocated qubits.
+
+        Args:
+            bit_string (list[bool|int]|string[0|1]): Computational basis state
+            qureg (Qureg|list[Qubit]): Quantum register determining the
+                ordering. Must contain all allocated qubits.
+
+        Returns:
+            Probability amplitude of the provided bit string.
+        """
+        bit_string = [bool(int(b)) for b in bit_string]
+        return self._simulator.get_amplitude(bit_string,
+                                             [qb.id for qb in qureg])
+
     def cheat(self):
         """
         Access the ordering of the qubits and the state vector directly.

projectq/backends/_sim/_simulator_test.py

@@ -16,6 +16,7 @@
 """
 
 import copy
+import math
 import numpy
 import pytest
 import random
@@ -38,7 +39,8 @@
                           BasicGate,
                           BasicMathGate,
                           QubitOperator,
-                          TimeEvolution)
+                          TimeEvolution,
+                          All)
 from projectq.meta import Control
 
 from projectq.backends import Simulator
@@ -204,6 +206,70 @@ def test_simulator_emulation(sim):
     Measure | (qubit1 + qubit2 + qubit3)
 
 
+def test_simulator_probability(sim):
+    eng = MainEngine(sim)
+    qubits = eng.allocate_qureg(6)
+    All(H) | qubits
+    eng.flush()
+    bits = [0, 0, 1, 0, 1, 0]
+    for i in range(6):
+        assert (eng.backend.get_probability(bits[:i], qubits[:i])
+                == pytest.approx(0.5**i))
+    extra_qubit = eng.allocate_qubit()
+    with pytest.raises(RuntimeError):
+        eng.backend.get_probability([0], extra_qubit)
+    del extra_qubit
+    All(H) | qubits
+    Ry(2 * math.acos(math.sqrt(0.3))) | qubits[0]
+    eng.flush()
+    assert eng.backend.get_probability([0], [qubits[0]]) == pytest.approx(0.3)
+    Ry(2 * math.acos(math.sqrt(0.4))) | qubits[2]
+    eng.flush()
+    assert eng.backend.get_probability([0], [qubits[2]]) == pytest.approx(0.4)
+    assert (eng.backend.get_probability([0, 0], qubits[:3:2])
+            == pytest.approx(0.12))
+    assert (eng.backend.get_probability([0, 1], qubits[:3:2])
+            == pytest.approx(0.18))
+    assert (eng.backend.get_probability([1, 0], qubits[:3:2])
+            == pytest.approx(0.28))
+    Measure | qubits
+
+
+def test_simulator_amplitude(sim):
+    eng = MainEngine(sim)
+    qubits = eng.allocate_qureg(6)
+    All(X) | qubits
+    All(H) | qubits
+    eng.flush()
+    bits = [0, 0, 1, 0, 1, 0]
+    assert eng.backend.get_amplitude(bits, qubits) == pytest.approx(1. / 8.)
+    bits = [0, 0, 0, 0, 1, 0]
+    assert eng.backend.get_amplitude(bits, qubits) == pytest.approx(-1. / 8.)
+    bits = [0, 1, 1, 0, 1, 0]
+    assert eng.backend.get_amplitude(bits, qubits) == pytest.approx(-1. / 8.)
+    All(H) | qubits
+    All(X) | qubits
+    Ry(2 * math.acos(0.3)) | qubits[0]
+    eng.flush()
+    bits = [0] * 6
+    assert eng.backend.get_amplitude(bits, qubits) == pytest.approx(0.3)
+    bits[0] = 1
+    assert (eng.backend.get_amplitude(bits, qubits)
+            == pytest.approx(math.sqrt(0.91)))
+    Measure | qubits
+    # raises if not all qubits are in the list:
+    with pytest.raises(RuntimeError):
+        eng.backend.get_amplitude(bits, qubits[:-1])
+    # doesn't just check for length:
+    with pytest.raises(RuntimeError):
+        eng.backend.get_amplitude(bits, qubits[:-1] + [qubits[0]])
+    extra_qubit = eng.allocate_qubit()
+    eng.flush()
+    # there is a new qubit now!
+    with pytest.raises(RuntimeError):
+        eng.backend.get_amplitude(bits, qubits)
+
+
 def test_simulator_expectation(sim):
     eng = MainEngine(sim, [])
     qureg = eng.allocate_qureg(3)