QubitOperator support

projectq/backends/_qracksim/_cpp/qracksimulator.hpp

@@ -40,6 +40,10 @@ class QrackSimulator{
     using StateVector = std::vector<std::complex<float>, aligned_allocator<std::complex<float>,64>>;
     using Map = std::map<unsigned, unsigned>;
     using RndEngine = qrack_rand_gen;
+    using Term = std::vector<std::pair<unsigned, char>>;
+    using TermsDict = std::vector<std::pair<Term, calc_type>>;
+    using ComplexTermsDict = std::vector<std::pair<Term, std::complex<calc_type>>>;
+    using Matrix = std::vector<std::vector<std::complex<double>, aligned_allocator<std::complex<double>, 64>>>;
     enum Qrack::QInterfaceEngine QrackEngine = Qrack::QINTERFACE_QUNIT;
     enum Qrack::QInterfaceEngine QrackSubengine1 = Qrack::QINTERFACE_QFUSION;
 #if ENABLE_OPENCL
@@ -422,6 +426,31 @@ class QrackSimulator{
         delete[] substateVec;
     }
 
+    void apply_qubit_operator(ComplexTermsDict const& td, std::vector<unsigned> const& ids){
+        for (auto const& term : td){
+            apply_term(term.first, term.second, ids, {});
+        }
+    }
+
+    calc_type get_expectation_value(TermsDict const& td, std::vector<unsigned> const& ids){
+        calc_type expectation = 0;
+
+        std::size_t mask = 0;
+        for (unsigned i = 0; i < ids.size(); i++){
+            mask |= 1UL << map_[ids[i]];
+        }
+
+        run();
+
+        Qrack::QInterfacePtr qRegOrig = qReg->Clone();
+        for (auto const& term : td){
+            expectation += diagonalize(term.first, term.second, ids);
+            qReg = qRegOrig->Clone();
+        }
+
+        return expectation;
+    }
+
     std::tuple<Map, StateVector> cheat(){
         if (qReg == NULL) {
             StateVector vec(1, 0.0);
@@ -556,6 +585,51 @@ class QrackSimulator{
         }
     }
 
+    void apply_term(Term const& term, std::complex<calc_type> coeff, std::vector<unsigned> const& ids,
+                    std::vector<unsigned> const& ctrl){
+
+        std::complex<double> I(0., 1.);
+        Matrix X = {{0., 1.}, {1., 0.}};
+        Matrix Y = {{0., -I}, {I, 0.}};
+        Matrix Z = {{1., 0.}, {0., -1.}};
+        std::vector<Matrix> gates = {X, Y, Z};
+
+        for (auto const& local_op : term){
+            unsigned id = ids[local_op.first];
+            auto temp = gates[local_op.second - 'X'];
+            for (unsigned i = 0; i < 4; i++) {
+                temp[i / 2][i % 2] *= -coeff;
+            }
+            apply_controlled_gate(temp, {id}, ctrl);
+        }
+    }
+
+    calc_type diagonalize(Term const& term, std::complex<calc_type> coeff, std::vector<unsigned> const& ids){
+        calc_type expectation = 1;
+        calc_type angle = arg(coeff);
+        calc_type len = abs(coeff);
+        std::complex<double> phaseFac(cos(angle), sin(angle));
+        bitCapInt idPower;
+
+        std::complex<double> I(0., 1.);
+        Matrix X = {{M_SQRT1_2, M_SQRT1_2}, {M_SQRT1_2, -M_SQRT1_2}};
+        Matrix Y = {{M_SQRT1_2, -M_SQRT1_2 * I}, {M_SQRT1_2, M_SQRT1_2 * I}};
+        Matrix Z = {{1., 0.}, {0., 1.}};
+        std::vector<Matrix> gates = {X, Y, Z};
+
+        for (auto const& local_op : term){
+            unsigned id = ids[local_op.first];
+            auto temp = gates[local_op.second - 'X'];
+            for (unsigned i = 0; i < 4; i++) {
+                temp[i / 2][i % 2] *= phaseFac;
+            }
+            apply_controlled_gate(temp, {id}, {});
+            idPower = 1U << map_[id];
+            expectation *= (qReg->ProbMask(idPower, 0) - qReg->ProbMask(idPower, idPower));
+        }
+        return len * expectation;
+    }
+
     Map map_;
     std::shared_ptr<RndEngine> rnd_eng_;
     Qrack::QInterfacePtr qReg;

projectq/backends/_qracksim/_qracksim.cpp

@@ -51,6 +51,8 @@ PYBIND11_PLUGIN(_qracksim) {
         .def("apply_controlled_dec", &QrackSimulator::apply_controlled_dec)
         .def("apply_controlled_mul", &QrackSimulator::apply_controlled_mul)
         .def("apply_controlled_div", &QrackSimulator::apply_controlled_div)
+        .def("get_expectation_value", &QrackSimulator::get_expectation_value)
+        .def("apply_qubit_operator", &QrackSimulator::apply_qubit_operator)
         .def("get_probability", &QrackSimulator::get_probability)
         .def("get_amplitude", &QrackSimulator::get_amplitude)
         .def("set_wavefunction", &QrackSimulator::set_wavefunction)

projectq/backends/_qracksim/_simulator.py

@@ -1,4 +1,3 @@
-#   Copyright 2017 ProjectQ-Framework (www.projectq.ch)
 #
 #   Licensed under the Apache License, Version 2.0 (the "License");
 #   you may not use this file except in compliance with the License.
@@ -31,7 +30,8 @@
                           Deallocate,
                           UniformlyControlledRy,
                           UniformlyControlledRz,
-                          StatePreparation)
+                          StatePreparation,
+                          QubitOperator)
 from projectq.libs.math import (AddConstant,
                                 AddConstantModN,
                                 MultiplyByConstantModN)
@@ -119,6 +119,8 @@ def is_available(self, cmd):
             if (isinstance(cmd.gate, StatePreparation) and not cmd.control_qubits):
                 # Qrack has inexpensive ways of preparing a partial state, without controls.
                 return True
+            elif (isinstance(cmd.gate, QubitOperator) and not np.isclose(1, np.absolute(cmd.gate.coefficient))):
+                return True
         except:
             pass
 
@@ -154,12 +156,81 @@ def _convert_logical_to_mapped_qureg(self, qureg):
             return qureg
 
     def get_expectation_value(self, qubit_operator, qureg):
-        # To maintain compatibility with default Simulator, for the moment.
-        pass
+        """
+        Get the expectation value of qubit_operator w.r.t. the current wave
+        function represented by the supplied quantum register.
+
+        Args:
+            qubit_operator (projectq.ops.QubitOperator): Operator to measure.
+            qureg (list[Qubit],Qureg): Quantum bits to measure.
+
+        Returns:
+            Expectation value
+
+        Note:
+            Make sure all previous commands (especially allocations) have
+            passed through the compilation chain (call main_engine.flush() to
+            make sure).
+
+        Note:
+            If there is a mapper present in the compiler, this function
+            automatically converts from logical qubits to mapped qubits for
+            the qureg argument.
+
+        Raises:
+            Exception: If `qubit_operator` acts on more qubits than present in
+                the `qureg` argument.
+        """
+        qureg = self._convert_logical_to_mapped_qureg(qureg)
+        num_qubits = len(qureg)
+        for term, _ in qubit_operator.terms.items():
+            if not term == () and term[-1][0] >= num_qubits:
+                raise Exception("qubit_operator acts on more qubits than "
+                                "contained in the qureg.")
+        operator = [(list(term), coeff) for (term, coeff)
+                    in qubit_operator.terms.items()]
+        return self._simulator.get_expectation_value(operator,
+                                                     [qb.id for qb in qureg])
 
     def apply_qubit_operator(self, qubit_operator, qureg):
-        # To maintain compatibility with default Simulator, for the moment.
-        pass
+        """
+        Apply a (possibly non-unitary) qubit_operator to the current wave
+        function represented by the supplied quantum register.
+
+        Args:
+            qubit_operator (projectq.ops.QubitOperator): Operator to apply.
+            qureg (list[Qubit],Qureg): Quantum bits to which to apply the
+                operator.
+
+        Raises:
+            Exception: If `qubit_operator` acts on more qubits than present in
+                the `qureg` argument.
+
+        Warning:
+            This function allows applying non-unitary gates and it will not
+            re-normalize the wave function! It is for numerical experiments
+            only and should not be used for other purposes.
+
+        Note:
+            Make sure all previous commands (especially allocations) have
+            passed through the compilation chain (call main_engine.flush() to
+            make sure).
+
+        Note:
+            If there is a mapper present in the compiler, this function
+            automatically converts from logical qubits to mapped qubits for
+            the qureg argument.
+        """
+        qureg = self._convert_logical_to_mapped_qureg(qureg)
+        num_qubits = len(qureg)
+        for term, _ in qubit_operator.terms.items():
+            if not term == () and term[-1][0] >= num_qubits:
+                raise Exception("qubit_operator acts on more qubits than "
+                                "contained in the qureg.")
+        operator = [(list(term), coeff) for (term, coeff)
+                    in qubit_operator.terms.items()]
+        return self._simulator.apply_qubit_operator(operator,
+                                                    [qb.id for qb in qureg])
 
     def get_probability(self, bit_string, qureg):
         """
@@ -357,6 +428,9 @@ def _handle(self, cmd):
             self._simulator.apply_controlled_sqrtswap(ids1, ids2,
                                                       [qb.id for qb in
                                                        cmd.control_qubits])
+        elif isinstance(cmd.gate, QubitOperator):
+            ids = [qb.id for qb in cmd.qubits[0]]
+            self.apply_qubit_operator(cmd.gate, ids)
         elif isinstance(cmd.gate, AddConstant) or isinstance(cmd.gate, AddConstantModN):
             #Unless there's a carry, the only unitary addition is mod (2^len(ids))
             ids = [qb.id for qr in cmd.qubits for qb in qr]

projectq/backends/_qracksim/_simulator_test.py

@@ -45,6 +45,9 @@
 from projectq.backends import Simulator
 
 
+tolerance = 1e-6
+
+
 def test_is_qrack_simulator_present():
     _qracksim = pytest.importorskip("projectq.backends._qracksim._qracksim")
     import projectq.backends._qracksim._qracksim as _
@@ -640,3 +643,64 @@ def test_uniformly_controlled_r(sim, gate_classes):
     All(Measure) | correct_qb + correct_ctrl_qureg
     test_eng.flush(deallocate_qubits=True)
     correct_eng.flush(deallocate_qubits=True)
+
+def test_get_expectation_value(sim):
+    num_qubits = 2
+    test_eng = MainEngine(sim)
+    test_qureg = test_eng.allocate_qureg(num_qubits)
+    test_eng.flush()
+
+    qubit_op = QubitOperator("X0", 1)
+    test_eng.backend.set_wavefunction([1 / math.sqrt(2), 1 / math.sqrt(2), 0, 0],
+                                      test_qureg)
+    test_eng.flush()
+    assert(sim.get_expectation_value(qubit_op, test_qureg) == pytest.approx(1, rel=tolerance, abs=tolerance))
+    test_eng.backend.set_wavefunction([1 / math.sqrt(2), -1 / math.sqrt(2), 0, 0],
+                                      test_qureg)
+    test_eng.flush()
+    assert(sim.get_expectation_value(qubit_op, test_qureg) == pytest.approx(-1, rel=tolerance, abs=tolerance))
+
+    qubit_op = QubitOperator("Y0", 1)
+    test_eng.backend.set_wavefunction([1 / math.sqrt(2), 1j / math.sqrt(2), 0, 0],
+                                      test_qureg)
+    test_eng.flush()
+    assert(sim.get_expectation_value(qubit_op, test_qureg) == pytest.approx(1, rel=tolerance, abs=tolerance))
+    test_eng.backend.set_wavefunction([1 / math.sqrt(2), -1j / math.sqrt(2), 0, 0],
+                                      test_qureg)
+    test_eng.flush()
+    assert(sim.get_expectation_value(qubit_op, test_qureg) == pytest.approx(-1, rel=tolerance, abs=tolerance))
+
+    qubit_op = QubitOperator("Z0", 1)
+    test_eng.backend.set_wavefunction([1, 0, 0, 0],
+                                      test_qureg)
+    test_eng.flush()
+    assert(sim.get_expectation_value(qubit_op, test_qureg) == pytest.approx(1, rel=tolerance, abs=tolerance))
+    test_eng.backend.set_wavefunction([0, 1, 0, 0],
+                                      test_qureg)
+    test_eng.flush()
+    assert(sim.get_expectation_value(qubit_op, test_qureg) == pytest.approx(-1, rel=tolerance, abs=tolerance))
+
+    qubit_op = QubitOperator("Z0", 1 / 4)
+    test_eng.backend.set_wavefunction([1, 0, 0, 0],
+                                      test_qureg)
+    test_eng.flush()
+    assert(sim.get_expectation_value(qubit_op, test_qureg) == pytest.approx(1 / 4, rel=tolerance, abs=tolerance))
+    test_eng.backend.set_wavefunction([0, 1, 0, 0],
+                                      test_qureg)
+    test_eng.flush()
+    assert(sim.get_expectation_value(qubit_op, test_qureg) == pytest.approx(-1 / 4, rel=tolerance, abs=tolerance))
+
+    qubit_op = QubitOperator("Z0 Z1", 1)
+    test_eng.backend.set_wavefunction([1, 0, 0, 0],
+                                      test_qureg)
+    test_eng.flush()
+    assert(sim.get_expectation_value(qubit_op, test_qureg) == pytest.approx(1, rel=tolerance, abs=tolerance))
+    X | test_qureg[0]
+    test_eng.flush()
+    assert(sim.get_expectation_value(qubit_op, test_qureg) == pytest.approx(-1, rel=tolerance, abs=tolerance))
+    X | test_qureg[1]
+    test_eng.flush()
+    assert(sim.get_expectation_value(qubit_op, test_qureg) == pytest.approx(1, rel=tolerance, abs=tolerance))
+    X | test_qureg[0]
+    test_eng.flush()
+    assert(sim.get_expectation_value(qubit_op, test_qureg) == pytest.approx(-1, rel=tolerance, abs=tolerance))

projectq/setups/decompositions/qubitop2onequbit_test.py

@@ -18,7 +18,7 @@
 
 from projectq import MainEngine
 # Qrack simulator does not yet support QubitOperator, so always use the default simulator:
-from projectq.backends._sim import Simulator
+from projectq.backends import Simulator
 from projectq.cengines import (AutoReplacer, DecompositionRuleSet, DummyEngine,
                                InstructionFilter)
 from projectq.meta import Control