Merge pull request #136 from C2QA/feature/two-mode-sum-gate

initial commit of two-mode sum gate

c2qa/circuit.py

@@ -359,7 +359,7 @@ def cv_sq2(self, theta, qumode_a, qumode_b, duration=100, unit="ns"):
             ),
             qargs=qumode_a + qumode_b,
         )
-    
+
     def cv_sq3(self, theta, qumode_a, qumode_b, qumode_c, duration=100, unit="ns"):
         """Three-mode squeezing gate
 
@@ -788,6 +788,61 @@ def cv_testqubitorderf(self, phi, qubit_1, qubit_2, duration=100, unit="ns"):
             qargs=[qubit_1] + [qubit_2],
         )
 
+    def cv_sum(self, scale, qumode_a, qumode_b, duration=100, unit="ns"):
+        """Two-mode sum gate.
+
+        Args:
+            scale (real): arbitrary real scale factor
+            qumode_a (list): list of qubits representing qumode
+            qumode_b (list): list of qubits representing qumode
+
+        Returns:
+            Instruction: QisKit instruction
+        """
+        return self.append(
+            ParameterizedUnitaryGate(
+                self.ops.sum,
+                [scale],
+                cutoffs=[
+                    QumodeRegister.calculate_cutoff(len(qumode_a)),
+                    QumodeRegister.calculate_cutoff(len(qumode_b)),
+                ],
+                num_qubits=len(qumode_a) + len(qumode_b),
+                label="sum",
+                duration=duration,
+                unit=unit,
+            ),
+            qargs=qumode_a + qumode_b,
+        )
+
+    def cv_c_sum(self, scale, qumode_a, qumode_b, qubit, duration=100, unit="ns"):
+        """Conditional two-mode sum gate.
+
+        Args:
+            scale (real): arbitrary real scale factor
+            qumode_a (list): list of qubits representing qumode
+            qumode_b (list): list of qubits representing qumode
+            qubit (Qubit): control Qubit
+
+        Returns:
+            Instruction: QisKit instruction
+        """
+        return self.append(
+            ParameterizedUnitaryGate(
+                self.ops.csum,
+                [scale],
+                cutoffs=[
+                    QumodeRegister.calculate_cutoff(len(qumode_a)),
+                    QumodeRegister.calculate_cutoff(len(qumode_b)),
+                ],
+                num_qubits=len(qumode_a) + len(qumode_b) + 1,
+                label="cSum",
+                duration=duration,
+                unit=unit,
+            ),
+            qargs=qumode_a + qumode_b + [qubit],
+        )
+
     def measure_z(self, qubit, cbit, duration=100, unit="ns"):
         """Measure qubit in z using probe qubits
 

c2qa/operators.py

@@ -31,11 +31,11 @@ def get_a2(self, cutoff_a: int, cutoff_b: int):
     def get_b1(self, cutoff_a: int, cutoff_b: int, cutoff_c: int):
         kron_ab = scipy.sparse.kron(self.get_a(cutoff_a), self.get_eye(cutoff_b))
         return scipy.sparse.kron(kron_ab, self.get_eye(cutoff_c)).tocsc()
-    
+
     def get_b2(self, cutoff_a: int, cutoff_b: int, cutoff_c: int):
         kron_ab = scipy.sparse.kron(self.get_eye(cutoff_a), self.get_a(cutoff_b))
         return scipy.sparse.kron(kron_ab, self.get_eye(cutoff_c)).tocsc()
-    
+
     def get_b3(self, cutoff_a: int, cutoff_b: int, cutoff_c: int):
         kron_ab = scipy.sparse.kron(self.get_eye(cutoff_a), self.get_eye(cutoff_b))
         return scipy.sparse.kron(kron_ab, self.get_a(cutoff_c)).tocsc()
@@ -44,7 +44,11 @@ def get_a12(self, cutoff_a: int, cutoff_b: int):
         return self.get_a1(cutoff_a, cutoff_b) * self.get_a2(cutoff_a, cutoff_b)
 
     def get_b123(self, cutoff_a: int, cutoff_b: int, cutoff_c: int):
-        return self.get_b1(cutoff_a, cutoff_b, cutoff_c) * self.get_b2(cutoff_a, cutoff_b, cutoff_c) * self.get_b3(cutoff_a, cutoff_b, cutoff_c)
+        return (
+            self.get_b1(cutoff_a, cutoff_b, cutoff_c)
+            * self.get_b2(cutoff_a, cutoff_b, cutoff_c)
+            * self.get_b3(cutoff_a, cutoff_b, cutoff_c)
+        )
 
     def get_a_dag(self, cutoff: int):
         """Creation operator"""
@@ -65,18 +69,22 @@ def get_a2_dag(self, cutoff_a: int, cutoff_b: int):
 
     def get_b1_dag(self, cutoff_a: int, cutoff_b: int, cutoff_c: int):
         return self.get_b1(cutoff_a, cutoff_b, cutoff_c).conjugate().transpose().tocsc()
-    
+
     def get_b2_dag(self, cutoff_a: int, cutoff_b: int, cutoff_c: int):
         return self.get_b2(cutoff_a, cutoff_b, cutoff_c).conjugate().transpose().tocsc()
-    
+
     def get_b3_dag(self, cutoff_a: int, cutoff_b: int, cutoff_c: int):
         return self.get_b3(cutoff_a, cutoff_b, cutoff_c).conjugate().transpose().tocsc()
 
     def get_a12_dag(self, cutoff_a: int, cutoff_b: int):
         return self.get_a1_dag(cutoff_a, cutoff_b) * self.get_a2_dag(cutoff_a, cutoff_b)
-    
+
     def get_b123_dag(self, cutoff_a: int, cutoff_b: int, cutoff_c: int):
-        return self.get_b1_dag(cutoff_a, cutoff_b, cutoff_c) * self.get_b2_dag(cutoff_a, cutoff_b, cutoff_c) * self.get_b3_dag(cutoff_a, cutoff_b, cutoff_c)
+        return (
+            self.get_b1_dag(cutoff_a, cutoff_b, cutoff_c)
+            * self.get_b2_dag(cutoff_a, cutoff_b, cutoff_c)
+            * self.get_b3_dag(cutoff_a, cutoff_b, cutoff_c)
+        )
 
     def get_a12dag(self, cutoff_a: int, cutoff_b: int):
         return self.get_a1(cutoff_a, cutoff_b) * self.get_a2_dag(cutoff_a, cutoff_b)
@@ -159,9 +167,10 @@ def s3(self, theta, cutoff_a, cutoff_b, cutoff_c):
             csc_matrix: operator matrix
         """
 
-        arg = (numpy.conjugate(theta * 1j) * self.get_b123_dag(cutoff_a, cutoff_b, cutoff_c)) - (
-            theta * 1j * self.get_b123(cutoff_a, cutoff_b, cutoff_c)
-        )
+        arg = (
+            numpy.conjugate(theta * 1j)
+            * self.get_b123_dag(cutoff_a, cutoff_b, cutoff_c)
+        ) - (theta * 1j * self.get_b123(cutoff_a, cutoff_b, cutoff_c))
 
         return scipy.sparse.linalg.expm(arg)
 
@@ -497,3 +506,42 @@ def gate_from_matrix(self, matrix):
             csc_matrix: operator matrix
         """
         return scipy.sparse.csc_matrix(matrix)
+
+    def sum(self, scale, cutoff_a, cutoff_b):
+        """Two-qumode sum gate
+
+        Args:
+            scale (real): arbitrary scale factor
+
+        Returns:
+            csc_matrix: operator matrix
+        """
+        # TODO verify below implementation
+        #     equation 205 from https://arxiv.org/pdf/2407.10381
+        #     vs equation 4 from https://journals.aps.org/prl/pdf/10.1103/PhysRevLett.88.097904
+        #   Equation 205 with matrix multiplication is not unitary, how to handle different cutoffs
+        a_mat = self.get_a(cutoff_a) + self.get_a_dag(cutoff_a)
+        b_mat = self.get_a_dag(cutoff_b) - self.get_a(cutoff_b)
+        arg = (scale / 2) * (scipy.sparse.kron(a_mat, b_mat))
+
+        return scipy.sparse.linalg.expm(arg)
+
+    def csum(self, scale, cutoff_a, cutoff_b):
+        """Conditional two-qumode sum gate
+
+        Args:
+            scale (real): arbitrary scale factor
+
+        Returns:
+            csc_matrix: operator matrix
+        """
+        # TODO verify below implementation
+        #     equation 205 from https://arxiv.org/pdf/2407.10381
+        #     vs equation 4 from https://journals.aps.org/prl/pdf/10.1103/PhysRevLett.88.097904
+        #   Equation 205 with matrix multiplication is not unitary, how to handle different cutoffs
+        a_mat = self.get_a(cutoff_a) + self.get_a_dag(cutoff_a)
+        b_mat = self.get_a_dag(cutoff_b) - self.get_a(cutoff_b)
+        arg = (scale / 2) * (scipy.sparse.kron(a_mat, b_mat))
+        arg = scipy.sparse.kron(zQB, arg)
+
+        return scipy.sparse.linalg.expm(arg)

tests/test_gates.py

@@ -359,6 +359,26 @@ def test_multiboson_sampling(capsys):
         assert result.success
 
 
+def test_two_mode_sum():
+    circuit, qmr = create_unconditional()
+
+    z = random.random()
+    circuit.cv_sum(z, qmr[0], qmr[1])
+
+    state, result, fock_counts = c2qa.util.simulate(circuit)
+    assert_changed(state, result)
+
+
+def test_conditional_two_mode_sum():
+    circuit, qmr, qbr = create_conditional()
+
+    z = random.random()
+    circuit.cv_c_sum(z, qmr[0], qmr[1], qbr[0])
+
+    state, result, fock_counts = c2qa.util.simulate(circuit)
+    assert_changed(state, result)
+
+
 def random_real_and_complex():
     return [
         random.random(),

tests/test_operators.py

@@ -46,6 +46,16 @@ def test_s(self):
     def test_s2(self):
         assert is_unitary_matrix(self.ops.s2(random.random(), self.cutoff, self.cutoff))
 
+    def test_sum(self):
+        assert is_unitary_matrix(
+            self.ops.sum(random.random(), self.cutoff, self.cutoff)
+        )
+
+    def test_csum(self):
+        assert is_unitary_matrix(
+            self.ops.csum(random.random(), self.cutoff, self.cutoff)
+        )
+
 
 class TestMatrices:
     """Test that the operators produce the values we expect.
@@ -128,3 +138,15 @@ def test_s2(self):
         rand = self.ops.s2(random.random(), self.cutoff, self.cutoff)
 
         assert not allclose(one, rand)
+
+    def test_sum(self):
+        one = self.ops.sum(1, self.cutoff, self.cutoff)
+        rand = self.ops.sum(random.random(), self.cutoff, self.cutoff)
+
+        assert not allclose(one, rand)
+
+    def test_csum(self):
+        one = self.ops.csum(1, self.cutoff, self.cutoff)
+        rand = self.ops.csum(random.random(), self.cutoff, self.cutoff)
+
+        assert not allclose(one, rand)