Feature/isometries

projectq/backends/_sim/_cppkernels/simulator.hpp

@@ -208,6 +208,58 @@ class Simulator{
             fused_gates_ = fused_gates;
     }
 
+    template <class M>
+    void apply_uniformly_controlled_gate(std::vector<M> &unitaries,
+                                         unsigned target_id,
+                                         std::vector<unsigned> choice_ids,
+                                         std::vector<unsigned> ctrl_ids){
+        run();
+        std::size_t n = vec_.size();
+        std::size_t dist = 1UL << map_[target_id];
+
+        auto mask = get_control_mask(ctrl_ids);
+
+        #pragma omp parallel for collapse(2) schedule(static)
+        for(std::size_t high = 0; high < n; high += 2*dist){
+            for(std::size_t low = 0; low < dist; ++low){
+                std::size_t entry = high+low;
+                if((entry&mask) == mask) {
+                    unsigned u = 0;
+                    for(std::size_t i = 0; i < choice_ids.size(); ++i)
+                        u |= ((entry >> map_[choice_ids[i]]) & 1) << i;
+
+                    auto &m = unitaries[u];
+                    std::complex<double> v[2];
+                    v[0] = vec_[entry];
+                    v[1] = vec_[entry + dist];
+                    vec_[entry]        = v[0]*m[0][0] + v[1]*m[0][1];
+                    vec_[entry + dist] = v[0]*m[1][0] + v[1]*m[1][1];
+                }
+            }
+        }
+    }
+
+    template <class M>
+    void apply_diagonal_gate(std::vector<calc_type> angles,
+                             std::vector<unsigned> ids,
+                             std::vector<unsigned> ctrl_ids)
+    {
+        run();
+        std::size_t n = vec_.size();
+        complex_type I(0., 1.);
+        auto mask = get_control_mask(ctrl_ids);
+
+        #pragma omp parallel for schedule(static)
+        for(std::size_t entry = 0; entry < n; ++entry) {
+            if((entry&mask) == mask) {
+                unsigned a = 0;
+                for(std::size_t i = 0; i < ids.size(); ++i)
+                    a |= ((entry >> map_[ids[i]]) & 1) << i;
+                vec_[entry] *= std::exp(I * angles[a]);
+            }
+        }
+    }
+
     template <class F, class QuReg>
     void emulate_math(F const& f, QuReg quregs, std::vector<unsigned> ctrl,
                       unsigned num_threads=1){

projectq/backends/_sim/_cppsim.cpp

@@ -49,6 +49,8 @@ PYBIND11_PLUGIN(_cppsim) {
         .def("is_classical", &Simulator::is_classical)
         .def("measure_qubits", &Simulator::measure_qubits_return)
         .def("apply_controlled_gate", &Simulator::apply_controlled_gate<MatrixType>)
+        .def("apply_uniformly_controlled_gate", &Simulator::apply_uniformly_controlled_gate<MatrixType>)
+        .def("apply_diagonal_gate", &Simulator::apply_diagonal_gate<MatrixType>)
         .def("emulate_math", &emulate_math_wrapper<QuRegs>)
         .def("get_expectation_value", &Simulator::get_expectation_value)
         .def("apply_qubit_operator", &Simulator::apply_qubit_operator)

projectq/backends/_sim/_pysim.py

@@ -20,6 +20,7 @@
 
 import random
 import numpy as _np
+import cmath
 
 
 class Simulator(object):
@@ -397,6 +398,43 @@ def apply_controlled_gate(self, m, ids, ctrlids):
             pos = [self._map[ID] for ID in ids]
             self._multi_qubit_gate(m, pos, mask)
 
+    def apply_uniformly_controlled_gate(self, unitaries, target_id,
+                                        choice_ids, ctrl_ids):
+        choice_pos = [self._map[ID] for ID in choice_ids]
+        pos = self._map[target_id]
+        mask = self._get_control_mask(ctrl_ids)
+
+        def kernel(u, d, m):
+            return u * m[0][0] + d * m[0][1], u * m[1][0] + d * m[1][1]
+
+        dist = 1 << pos
+        n = len(self._state)
+        for high in range(0, n, 2*dist):
+            for low in range(0, dist):
+                entry = high+low
+                if (entry & mask) == mask:
+                    u = 0
+                    for i in range(len(choice_pos)):
+                        u |= ((entry >> choice_pos[i]) & 1) << i
+                    id1 = entry
+                    id2 = entry + dist
+                    self._state[id1], self._state[id2] = kernel(
+                        self._state[id1],
+                        self._state[id2],
+                        unitaries[u])
+
+    def apply_diagonal_gate(self, angles, ids, ctrlids):
+        pos = [self._map[ID] for ID in ids]
+        mask = self._get_control_mask(ctrlids)
+
+        n = len(self._state)
+        for entry in range(n):
+            if (entry & mask) == mask:
+                a = 0
+                for i in range(len(pos)):
+                    a |= ((entry >> pos[i]) & 1) << i
+                self._state[entry] *= cmath.exp(1j*angles[a])
+
     def _single_qubit_gate(self, m, pos, mask):
         """
         Applies the single qubit gate matrix m to the qubit at position `pos`

projectq/backends/_sim/_simulator.py

@@ -19,19 +19,23 @@
 """
 
 import math
+import cmath
 import random
+from projectq.types import WeakQubitRef
 from projectq.cengines import BasicEngine
 from projectq.meta import get_control_count, LogicalQubitIDTag
-from projectq.ops import (NOT,
+from projectq.ops import (All,
+                          NOT,
                           H,
                           R,
                           Measure,
                           FlushGate,
                           Allocate,
                           Deallocate,
                           BasicMathGate,
-                          TimeEvolution)
-from projectq.types import WeakQubitRef
+                          TimeEvolution,
+                          UniformlyControlledGate,
+                          DiagonalGate)
 
 try:
     from ._cppsim import Simulator as SimulatorBackend
@@ -104,7 +108,9 @@ def is_available(self, cmd):
         if (cmd.gate == Measure or cmd.gate == Allocate or
                 cmd.gate == Deallocate or
                 isinstance(cmd.gate, BasicMathGate) or
-                isinstance(cmd.gate, TimeEvolution)):
+                isinstance(cmd.gate, TimeEvolution) or
+                isinstance(cmd.gate, UniformlyControlledGate) or
+                isinstance(cmd.gate, DiagonalGate)):
             return True
         try:
             m = cmd.gate.matrix
@@ -399,6 +405,26 @@ def _handle(self, cmd):
             qubitids = [qb.id for qb in cmd.qubits[0]]
             ctrlids = [qb.id for qb in cmd.control_qubits]
             self._simulator.emulate_time_evolution(op, t, qubitids, ctrlids)
+        elif isinstance(cmd.gate, UniformlyControlledGate):
+            choice_ids = [qb.id for qb in cmd.qubits[0]]
+            ctrl_ids = [qb.id for qb in cmd.control_qubits]
+            target_id = cmd.qubits[1][0].id
+            unitaries = [gate.matrix.tolist() for gate in cmd.gate.gates]
+            assert len(unitaries) == 2**len(choice_ids)
+            self._simulator.apply_uniformly_controlled_gate(unitaries,
+                                                            target_id,
+                                                            choice_ids,
+                                                            ctrl_ids)
+            if cmd.gate.up_to_diagonal:
+                angles = [-cmath.phase(p) for p in cmd.gate.decomposition[1]]
+                ids = [target_id]+choice_ids
+                self._simulator.apply_diagonal_gate(angles, ids, [])
+        elif isinstance(cmd.gate, DiagonalGate):
+            ids = [q.id for qr in cmd.qubits for q in qr]
+            ctrlids = [qb.id for qb in cmd.control_qubits]
+            angles = cmd.gate.angles
+            assert len(angles) == 2**len(ids)
+            self._simulator.apply_diagonal_gate(angles, ids, ctrlids)
         elif len(cmd.gate.matrix) <= 2 ** 5:
             matrix = cmd.gate.matrix
             ids = [qb.id for qr in cmd.qubits for qb in qr]

projectq/backends/_sim/_simulator_test.py

@@ -19,6 +19,7 @@
 
 import copy
 import math
+import cmath
 import numpy
 import pytest
 import random
@@ -31,7 +32,8 @@
                                LocalOptimizer, NotYetMeasuredError)
 from projectq.ops import (All, Allocate, BasicGate, BasicMathGate, CNOT,
                           Command, H, MatrixGate, Measure, QubitOperator,
-                          Rx, Ry, Rz, S, TimeEvolution, Toffoli, X, Y, Z)
+                          Rx, Ry, Rz, S, T, TimeEvolution, Toffoli, X, Y, Z,
+                          DiagonalGate, UniformlyControlledGate)
 from projectq.meta import Control, Dagger, LogicalQubitIDTag
 from projectq.types import WeakQubitRef
 
@@ -521,6 +523,92 @@ def build_matrix(list_single_matrices):
                           init_wavefunction)
 
 
+def test_simulator_apply_diagonal_gate(sim):
+    eng = MainEngine(sim)
+    qureg = eng.allocate_qureg(4)
+    eng.flush()
+
+    target_1 = qureg[0]
+    control = qureg[1]
+    target_0 = qureg[2]
+    empty = qureg[3]
+
+    wf = [1./4.]*(1 << 4)
+    eng.backend.set_wavefunction(wf, qureg)
+
+    D = DiagonalGate(angles=range(4))
+    with Control(eng, control):
+        D | (target_0, target_1)
+
+    eng.flush()
+    qbit_to_bit_map, final_wavefunction = copy.deepcopy(eng.backend.cheat())
+    All(Measure) | qureg
+
+    desired_state = [1./4.*cmath.exp(1j*i) for i in
+                     [0, 0, 0, 2, 0, 0, 1, 3, 0, 0, 0, 2, 0, 0, 1, 3]]
+    assert numpy.allclose(final_wavefunction, desired_state)
+
+
+def test_simulator_apply_uniformly_controlled_gate():
+    eng = MainEngine()
+    qureg = eng.allocate_qureg(3)
+    eng.flush()
+    wf = [math.sqrt(1./8.)]*(1 << 3)
+    eng.backend.set_wavefunction(wf, qureg)
+
+    A = Rx(numpy.pi/5)
+    B = H
+    C = Rz(numpy.pi/5)
+    D = Ry(numpy.pi/3)
+    U = UniformlyControlledGate([A, B, C, D])
+    with Dagger(eng):
+        U | ([qureg[0], qureg[2]], qureg[1])
+    eng.flush()
+    qbit_to_bit_map, final_wavefunction = copy.deepcopy(eng.backend.cheat())
+    vec = numpy.array([final_wavefunction]).T
+    vec[[1, 2]] = vec[[2, 1]]  # reorder basis
+    vec[[5, 6]] = vec[[6, 5]]
+    reference = numpy.matrix(scipy.linalg.block_diag(A.matrix, B.matrix,
+                                                     C.matrix, D.matrix))
+    assert numpy.allclose(reference*vec, wf)
+
+
+def test_simulator_apply_uniformly_controlled_gate_with_control(sim):
+    eng = MainEngine(sim)
+    qureg = eng.allocate_qureg(5)
+    eng.flush()
+
+    control = qureg[0]
+    choice_1 = qureg[1]
+    target = qureg[2]
+    empty = qureg[3]
+    choice_0 = qureg[4]
+
+    gates = [X, H, S, T]
+    U = UniformlyControlledGate(gates)
+
+    id = Rz(0.0)
+    gates_equiv = [id, X, id, S, id, X, id, S,
+                   id, H, id, T, id, H, id, T]
+    U_equiv = UniformlyControlledGate(gates_equiv)
+
+    All(H) | qureg
+    with Control(eng, control):
+        U | ([choice_0, choice_1], target)
+
+    with Dagger(eng):
+        All(H) | qureg
+        U_equiv | (qureg[0:2]+qureg[3:5], target)
+
+    eng.flush()
+    qbit_to_bit_map, final_wavefunction = copy.deepcopy(eng.backend.cheat())
+    All(Measure) | qureg
+
+    print(final_wavefunction)
+    desired_state = [1.0] + [0.0]*31
+    assert numpy.allclose(final_wavefunction, desired_state)
+
+
 def test_simulator_set_wavefunction(sim, mapper):
     engine_list = [LocalOptimizer()]
     if mapper is not None:

projectq/libs/isometries/__init__.py

@@ -0,0 +1,7 @@
+from .decompositions import _decompose_diagonal_gate
+from .decompositions import _decompose_uniformly_controlled_gate
+from .decompositions import _decompose_isometry
+from .single_qubit_gate import _SingleQubitGate
+from .apply_decompositions import (_apply_isometry,
+                                   _apply_diagonal_gate,
+                                   _apply_uniformly_controlled_gate)