Repository: TeamGraphix/graphix

Issue:
Measurement pattern for Toffoli gate 
Implement measurement pattern for Toffoli gate in `graphix.Circuit` class for transpilation. Please consider:
- minimize the number of commands by choosing an appropriate decomposition before translating into a measurement pattern.
- Make sure the `Circuit.simulate_statevector()` runs with Toffoli gate.
- Please also implement 'opt' version (see, e.g. `Circuit._rz_command_opt` and `_rzz_command_opt):
    -  By performing Pauli measurement preprocessing for gate-level measurement patterns (eg `Circuit._rz_command`), it is possible to implement a significantly shorter measurement pattern for a gate, with the cost of added local Clifford commands (see `Circuit._rz_command_opt`).
    - Make sure you only perform Pauli measurements on the intermediate qubits (not input or output qubits of the gate-level pattern), so that subsequent pattern optimization works as intended.
    - Please add tests (in `tests/test_transpiler.py` see also `tests/test_pattern.py`) and make sure the pattern simulation results are as expected even after various forms of pattern optimizations (see `tests/test_pattern.py`). 

Please see our [contribution guide](https://github.com/TeamGraphix/graphix/blob/master/CONTRIBUTING.md) before starting, and do not forget to comment on this issue to show your interest. We are always happy to answer questions, so do not hesitate to ask on this page.

Thank you for your interest, and good luck! 

Discussion:
shinich1: @king-p3nguin hi, just want to see how it's going as UnitaryHack finishes early next week.
king-p3nguin: @shinich1 Thank you for the reminder. I am a bit busy and unsure I can meet the deadline, but I will work on this issue even if the UnitaryHack is finished.
natestemen: Nice work with this hack @king-p3nguin! Just a reminder that PRs are eligible for unitaryHACK bounties as long as the issue is closed before June 20th! The hackathon finishes EOD June 13th [AoE](https://time.is/Anywhere_on_Earth), and as long as the PR is up by this deadline, it is eligible! We understand it can take some time to get reviewed and make changes.
shinich1: @king-p3nguin Good that it works now! could you modify test cases to include ccz gate (perhaps add to random_circuit.py for now), and also add lines to docs (.rst files) where necessary?
also - perhaps you can use `-l 120` option for black to save up on line numbers in transpilar.py (see above)?
king-p3nguin: @shinich1 
Okay, I think I can do that. Unfortunately, there still seems to be some issue in the ccx commands...
I thought simply replacing the commands from the original pattern with each variable would work, but it did not. I wrote the following code to replace the commands. Maybe you can spot the bug.

```python
from graphix.transpiler import Circuit
import numpy as np

circuit = Circuit(3)
# toffoli
circuit.h(2)
circuit.cnot(1, 2)
circuit.rz(2, -np.pi / 4)
circuit.cnot(0, 2)
circuit.rz(2, np.pi / 4)
circuit.cnot(1, 2)
circuit.rz(1, np.pi / 4)
circuit.rz(2, -np.pi / 4)
circuit.cnot(0, 2)
circuit.cnot(0, 1)
circuit.rz(2, np.pi / 4)
circuit.rz(0, np.pi / 4)
circuit.rz(1, -np.pi / 4)
circuit.h(2)
circuit.cnot(0, 1)

# transpile
pattern = circuit.transpile()
pattern.standardize()
pattern.shift_signals()
pattern.perform_pauli_measurements()

dictonary = {  # for opt = True
    0: "ancilla[0]",
    1: "ancilla[1]",
    5: "ancilla[2]",
    9: "ancilla[3]",
    13: "ancilla[4]",
    21: "ancilla[5]",
    25: "control_node1",
    28: "target_node",
    29: "ancilla[6]",
    30: "control_node2",
}
# dictonary = { # for opt = False
#     0: "ancilla[0]",
#     1: "ancilla[1]",
#     2: "ancilla[2]",
#     3: "ancilla[3]",
#     4: "ancilla[4]",
#     5: "ancilla[5]",
#     6: "ancilla[6]",
#     7: "ancilla[7]",
#     8: "ancilla[8]",
#     9: "ancilla[9]",
#     10: "ancilla[10]",
#     11: "ancilla[11]",
#     12: "ancilla[12]",
#     13: "ancilla[13]",
#     14: "ancilla[14]",
#     15: "ancilla[15]",
#     16: "ancilla[16]",
#     17: "ancilla[17]",
#     18: "ancilla[18]",
#     19: "ancilla[19]",
#     20: "ancilla[20]",
#     21: "ancilla[21]",
#     22: "ancilla[22]",
#     23: "ancilla[23]",
#     24: "ancilla[24]",
#     26: "ancilla[25]",
#     27: "ancilla[26]",
#     29: "ancilla[27]",
#     25: "control_node1",
#     28: "target_node",
#     30: "control_node2",
# }
import numpy as np
for i in range(len(pattern.seq)):
    if pattern.seq[i][0] == "E":
        print(
            f'seq.append(["E", ({dictonary[pattern.seq[i][1][0]]}, {dictonary[pattern.seq[i][1][1]]})])'
        )
    elif pattern.seq[i][0] == "M":
        s_domain = pattern.seq[i][4]
        new_s_domain = []
        for sidx in s_domain:
            try:
                new_s_domain.append(dictonary[sidx])
            except KeyError:
                continue
        t_domain = pattern.seq[i][5]
        new_t_domain = []
        for tidx in t_domain:
            try:
                new_t_domain.append(dictonary[tidx])
            except KeyError:
                continue
        try:
            print(
                f'seq.append(["M", {dictonary[pattern.seq[i][1]]}, "{pattern.seq[i][2]}", {pattern.seq[i][3]}, [{", ".join(new_s_domain)}], [{", ".join(new_t_domain)}], {pattern.seq[i][6]}])'
            )
        except IndexError:
            print(
                f'seq.append(["M", {dictonary[pattern.seq[i][1]]}, "{pattern.seq[i][2]}", {pattern.seq[i][3]}, [{", ".join(new_s_domain)}], [{", ".join(new_t_domain)}]])'
            )
    elif pattern.seq[i][0] == "X":
        # remove duplicates
        _domain = np.array(pattern.seq[i][2])
        uind = np.unique(_domain)
        unique_domain = []
        for ind in uind:
            if np.mod(np.count_nonzero(_domain == ind), 2) == 1:
                unique_domain.append(ind)

        new_domain = []
        for idx in unique_domain:
            try:
                new_domain.append(dictonary[idx])
            except KeyError:
                continue
        if new_domain == []:
            continue
        print(
            f'seq.append(["X", {dictonary[pattern.seq[i][1]]}, [{", ".join(new_domain)}]])'
        )
    elif pattern.seq[i][0] == "Z":
        # remove duplicates
        _domain = np.array(pattern.seq[i][2])
        uind = np.unique(_domain)
        unique_domain = []
        for ind in uind:
            if np.mod(np.count_nonzero(_domain == ind), 2) == 1:
                unique_domain.append(ind)

        new_domain = []
        for idx in unique_domain:
            try:
                new_domain.append(dictonary[idx])
            except KeyError:
                continue
        if new_domain == []:
            continue
        print(
            f'seq.append(["Z", {dictonary[pattern.seq[i][1]]}, [{", ".join(new_domain)}]])'
        )
```
king-p3nguin: Also, I was a bit confused about how the "X" and "Z" commands work. `pattern.print_pattern()` seems to remove duplicate nodes in the list, but `StatevectorBackend.correct_byproduct(cmd)` does not?
shinich1: https://github.com/TeamGraphix/graphix/pull/56#issuecomment-1626050694

What's this code supposed to do? and what exactly is the issue? could you elaborate? 
shinich1: > Also, I was a bit confused about how the "X" and "Z" commands work. `pattern.print_pattern()` seems to remove duplicate nodes in the list, but `StatevectorBackend.correct_byproduct(cmd)` does not?

the results are the same, we are thinking about a way to make sure we simplify signal lists in an efficient way, not just when printing
king-p3nguin: > [#56 (comment)](https://github.com/TeamGraphix/graphix/pull/56#issuecomment-1626050694)
> 
> What's this code supposed to do? and what exactly is the issue? could you elaborate?

Sorry, I did not explain it enough. What I was doing is that I generated the commands from a circuit that has only one Toffoli gate, and replaced the node number in the commands with `ancilla[i]`, `control_node1`, `control_node2`, and `target_node` so that I can simply copy & paste directly into the code.
shinich1: I see, it might be worth looking into graph state for toffoli gate and check their LC form after pauli measurement, e.g. look for isolated nodes. Just to check, have you modified `pattern.perform_pauli_measurements()` so it leave input qubits unmeasured? 
you may need to find alternative way and perhaps try out ZX calculus to get correct pattern fragment - rather than copy and paste. 
king-p3nguin: Yes, I also tried skipping the measurement of input qubits and the results were the same.
I'm kind of stuck with this method, so will check ZX calculus.
king-p3nguin: By using PyZX and `graphix.generator.generate_from_graph()`, I was able to create commands for ccx gate (only for `opt = False`).
Commands I just created from `pattern.perform_pauli_measurements(leave_input=True)` does not seem to be working, so I will work on it later. (I might be missing something)

It took me a while to figure out how to convert a ZX-diagram into graph state. The process was more complicated than I expected. It might be a good idea to add to the example how to create a measurement pattern from a ZX-diagram created by pyzx. (It would also be a good use case for `graphix.generator.generate_from_graph()`)
king-p3nguin: I was not able to reduce the nodes a lot for `opt=True`, but it works for now.

@shinich1 Will you review this PR?
shinich1: @king-p3nguin found a better pattern, please take a look. what worked was to manually perform Pauli measurements of the pattern using graph state sim while making sure no local Hadamard remain on either input or output nodes (because these cannot be commuted through next or previous commands), sometimes interleaving with equivalent_fill_nodes of the graph state simulator. 

the script I used (works also after perform_pauli_measurements)
```python
circuit = Circuit(3)
circuit.ccx(0, 1, 2)
pattern = circuit.transpile()
pattern.standardize()
pattern.shift_signals()
pattern.draw_graph()
nodes, edges = pattern.get_graph()
vops = pattern.get_vops()
print(nodes)
print(edges)
from graphix import GraphState
nodes2 = [i for i in range(21)]
mapping = {nodes[i]: nodes2[i] for i in range(21)}
edges2 = [(mapping[i], mapping[j]) for i,j in edges]
vop2 = dict()
for i in nodes:
    if i in vops.keys():
        vop2[mapping[i]] = vops[i]
    else:
        vop2[mapping[i]] = 0
g=GraphState(nodes=nodes2, edges=edges2, vops=vop2)

results = dict()
results[3] = g.measure_x(3)
g.equivalent_fill_node(2)
results[4] = g.measure_x(4)
results[6] = g.measure_x(6)
results[12] = g.measure_x(12)
g.equivalent_fill_node(1)
results[11] = g.measure_x(11)
results[19] = g.measure_x(19)
results[9] = g.measure_x(9)
# 17, 14
g.draw()
print(results)

```
shinich1: @king-p3nguin would you be happy to add type hints?
shinich1: @king-p3nguin looks great, thanks! if it's ready on your side, please squash and merge.

Patch:
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 2ffe97c..33f45ce 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -13,6 +13,7 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
   - Only `networkx` backend was available for pattern optimization.
   By setting the `use_rustworkx` option to True while using `Pattern.perform_pauli_measurements()`,
   graphix will run pattern optimization using `rustworkx` (#98)
+- Added `.ccx()` method to `graphix.Circuit`.
 
 ### Fixed
 
diff --git a/graphix/gflow.py b/graphix/gflow.py
index 908c45c..740169d 100644
--- a/graphix/gflow.py
+++ b/graphix/gflow.py
@@ -416,7 +416,7 @@ def find_odd_neighbor(graph, candidate, vertices):
     """
     out = []
     for c in candidate:
-        if np.mod(len(set(graph.neighbor(c)) ^ vertices), 2) == 1:
+        if np.mod(len(set(graph.neighbors(c)) ^ vertices), 2) == 1:
             out.append(c)
     return out
 
diff --git a/graphix/ops.py b/graphix/ops.py
index 25d824e..305a0a5 100644
--- a/graphix/ops.py
+++ b/graphix/ops.py
@@ -26,6 +26,18 @@ class Ops:
     cz = np.array([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, -1]])
     cnot = np.array([[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 0, 1], [0, 0, 1, 0]])
     swap = np.array([[1, 0, 0, 0], [0, 0, 1, 0], [0, 1, 0, 0], [0, 0, 0, 1]])
+    ccx = np.array(
+        [
+            [1, 0, 0, 0, 0, 0, 0, 0],
+            [0, 1, 0, 0, 0, 0, 0, 0],
+            [0, 0, 1, 0, 0, 0, 0, 0],
+            [0, 0, 0, 1, 0, 0, 0, 0],
+            [0, 0, 0, 0, 1, 0, 0, 0],
+            [0, 0, 0, 0, 0, 1, 0, 0],
+            [0, 0, 0, 0, 0, 0, 0, 1],
+            [0, 0, 0, 0, 0, 0, 1, 0],
+        ]
+    )
 
     @staticmethod
     def Rx(theta):
diff --git a/graphix/transpiler.py b/graphix/transpiler.py
index 58f0ec7..5340de6 100644
--- a/graphix/transpiler.py
+++ b/graphix/transpiler.py
@@ -3,7 +3,10 @@
 accepts desired gate operations and transpile into MBQC measurement patterns.
 
 """
+from __future__ import annotations
+
 from copy import deepcopy
+from typing import Optional, Sequence
 
 import numpy as np
 
@@ -25,7 +28,7 @@ class Circuit:
         List containing the gate sequence applied.
     """
 
-    def __init__(self, width):
+    def __init__(self, width: int):
         """
         Parameters
         ----------
@@ -35,7 +38,7 @@ class Circuit:
         self.width = width
         self.instruction = []
 
-    def cnot(self, control, target):
+    def cnot(self, control: int, target: int):
         """CNOT gate
 
         Parameters
@@ -50,7 +53,22 @@ class Circuit:
         assert control != target
         self.instruction.append(["CNOT", [control, target]])
 
-    def h(self, qubit):
+    def swap(self, qubit1: int, qubit2: int):
+        """SWAP gate
+
+        Parameters
+        ---------
+        qubit1 : int
+            first qubit to be swapped
+        qubit2 : int
+            second qubit to be swapped
+        """
+        assert qubit1 in np.arange(self.width)
+        assert qubit2 in np.arange(self.width)
+        assert qubit1 != qubit2
+        self.instruction.append(["SWAP", [qubit1, qubit2]])
+
+    def h(self, qubit: int):
         """Hadamard gate
 
         Parameters
@@ -61,7 +79,7 @@ class Circuit:
         assert qubit in np.arange(self.width)
         self.instruction.append(["H", qubit])
 
-    def s(self, qubit):
+    def s(self, qubit: int):
         """S gate
 
         Parameters
@@ -83,7 +101,7 @@ class Circuit:
         assert qubit in np.arange(self.width)
         self.instruction.append(["X", qubit])
 
-    def y(self, qubit):
+    def y(self, qubit: int):
         """Pauli Y gate
 
         Parameters
@@ -94,7 +112,7 @@ class Circuit:
         assert qubit in np.arange(self.width)
         self.instruction.append(["Y", qubit])
 
-    def z(self, qubit):
+    def z(self, qubit: int):
         """Pauli Z gate
 
         Parameters
@@ -105,7 +123,7 @@ class Circuit:
         assert qubit in np.arange(self.width)
         self.instruction.append(["Z", qubit])
 
-    def rx(self, qubit, angle):
+    def rx(self, qubit: int, angle: float):
         """X rotation gate
 
         Parameters
@@ -118,7 +136,7 @@ class Circuit:
         assert qubit in np.arange(self.width)
         self.instruction.append(["Rx", qubit, angle])
 
-    def ry(self, qubit, angle):
+    def ry(self, qubit: int, angle: float):
         """Y rotation gate
 
         Parameters
@@ -131,7 +149,7 @@ class Circuit:
         assert qubit in np.arange(self.width)
         self.instruction.append(["Ry", qubit, angle])
 
-    def rz(self, qubit, angle):
+    def rz(self, qubit: int, angle: float):
         """Z rotation gate
 
         Parameters
@@ -144,7 +162,7 @@ class Circuit:
         assert qubit in np.arange(self.width)
         self.instruction.append(["Rz", qubit, angle])
 
-    def rzz(self, control, target, angle):
+    def rzz(self, control: int, target: int, angle: float):
         r"""ZZ-rotation gate.
         Equivalent to the sequence
         CNOT(control, target),
@@ -156,7 +174,9 @@ class Circuit:
 
         Parameters
         ---------
-        qubit : int
+        control : int
+            control qubit
+        target : int
             target qubit
         angle : float
             rotation angle in radian
@@ -165,7 +185,24 @@ class Circuit:
         assert target in np.arange(self.width)
         self.instruction.append(["Rzz", [control, target], angle])
 
-    def i(self, qubit):
+    def ccx(self, control1: int, control2: int, target: int):
+        r"""CCX (Toffoli) gate.
+
+        Prameters
+        ---------
+        control1 : int
+            first control qubit
+        control2 : int
+            second control qubit
+        target : int
+            target qubit
+        """
+        assert control1 in np.arange(self.width)
+        assert control2 in np.arange(self.width)
+        assert target in np.arange(self.width)
+        self.instruction.append(["CCX", [control1, control2, target]])
+
+    def i(self, qubit: int):
         """identity (teleportation) gate
 
         Parameters
@@ -176,7 +213,7 @@ class Circuit:
         assert qubit in np.arange(self.width)
         self.instruction.append(["I", qubit])
 
-    def transpile(self, opt=False):
+    def transpile(self, opt: bool = False):
         """gate-to-MBQC transpile function.
 
         Parameters
@@ -201,6 +238,8 @@ class Circuit:
                 )
                 pattern.seq.extend(seq)
                 Nnode += 2
+            elif instr[0] == "SWAP":
+                out[instr[1][0]], out[instr[1][1]] = out[instr[1][1]], out[instr[1][0]]
             elif instr[0] == "I":
                 pass
             elif instr[0] == "H":
@@ -262,13 +301,34 @@ class Circuit:
                         "YZ-plane measurements not accepted and Rzz gate\
                         cannot be directly transpiled"
                     )
+            elif instr[0] == "CCX":
+                if opt:
+                    ancilla = [Nnode + i for i in range(11)]
+                    (
+                        out[instr[1][0]],
+                        out[instr[1][1]],
+                        out[instr[1][2]],
+                        seq,
+                    ) = self._ccx_command_opt(out[instr[1][0]], out[instr[1][1]], out[instr[1][2]], ancilla)
+                    pattern.seq.extend(seq)
+                    Nnode += 11
+                else:
+                    ancilla = [Nnode + i for i in range(18)]
+                    (
+                        out[instr[1][0]],
+                        out[instr[1][1]],
+                        out[instr[1][2]],
+                        seq,
+                    ) = self._ccx_command(out[instr[1][0]], out[instr[1][1]], out[instr[1][2]], ancilla)
+                    pattern.seq.extend(seq)
+                    Nnode += 18
             else:
                 raise ValueError("Unknown instruction, commands not added")
         pattern.output_nodes = out
         pattern.Nnode = Nnode
         return pattern
 
-    def standardize_and_transpile(self, opt=True):
+    def standardize_and_transpile(self, opt: bool = True):
         """gate-to-MBQC transpile function.
         Commutes all byproduct through gates, instead of through measurement
         commands, to generate standardized measurement pattern.
@@ -305,6 +365,9 @@ class Circuit:
                 self._instr.append(["XC", instr[1][1], seq[7][2]])
                 self._instr.append(["ZC", instr[1][1], seq[8][2]])
                 self._instr.append(["ZC", instr[1][0], seq[9][2]])
+            elif instr[0] == "SWAP":
+                out[instr[1][0]], out[instr[1][1]] = out[instr[1][1]], out[instr[1][0]]
+                self._instr.append(instr)
             elif instr[0] == "I":
                 pass
             elif instr[0] == "H":
@@ -461,7 +524,20 @@ class Circuit:
         pattern.seq = command_seq
         return pattern
 
-    def _commute_with_cnot(self, target):
+    def _commute_with_swap(self, target: int):
+        assert self._instr[target][0] in ["XC", "ZC"]
+        assert self._instr[target + 1][0] == "SWAP"
+        if self._instr[target][1] == self._instr[target + 1][1][0]:
+            self._instr[target][1] = self._instr[target + 1][1][1]
+            self._commute_with_following(target)
+        elif self._instr[target][1] == self._instr[target + 1][1][1]:
+            self._instr[target][1] = self._instr[target + 1][1][0]
+            self._commute_with_following(target)
+        else:
+            self._commute_with_following(target)
+        return target
+
+    def _commute_with_cnot(self, target: int):
         assert self._instr[target][0] in ["XC", "ZC"]
         assert self._instr[target + 1][0] == "CNOT"
         if self._instr[target][0] == "XC" and self._instr[target][1] == self._instr[target + 1][1][0]:  # control
@@ -478,7 +554,7 @@ class Circuit:
             self._commute_with_following(target)
         return target
 
-    def _commute_with_H(self, target):
+    def _commute_with_H(self, target: int):
         assert self._instr[target][0] in ["XC", "ZC"]
         assert self._instr[target + 1][0] == "H"
         if self._instr[target][1] == self._instr[target + 1][1]:
@@ -491,7 +567,7 @@ class Circuit:
         else:
             self._commute_with_following(target)
 
-    def _commute_with_S(self, target):
+    def _commute_with_S(self, target: int):
         assert self._instr[target][0] in ["XC", "ZC"]
         assert self._instr[target + 1][0] == "S"
         if self._instr[target][1] == self._instr[target + 1][1]:
@@ -503,7 +579,7 @@ class Circuit:
         self._commute_with_following(target)
         return target
 
-    def _commute_with_Rx(self, target):
+    def _commute_with_Rx(self, target: int):
         assert self._instr[target][0] in ["XC", "ZC"]
         assert self._instr[target + 1][0] == "Rx"
         if self._instr[target][1] == self._instr[target + 1][1]:
@@ -516,7 +592,7 @@ class Circuit:
         else:
             self._commute_with_following(target)
 
-    def _commute_with_Ry(self, target):
+    def _commute_with_Ry(self, target: int):
         assert self._instr[target][0] in ["XC", "ZC"]
         assert self._instr[target + 1][0] == "Ry"
         if self._instr[target][1] == self._instr[target + 1][1]:
@@ -526,7 +602,7 @@ class Circuit:
         else:
             self._commute_with_following(target)
 
-    def _commute_with_Rz(self, target):
+    def _commute_with_Rz(self, target: int):
         assert self._instr[target][0] in ["XC", "ZC"]
         assert self._instr[target + 1][0] == "Rz"
         if self._instr[target][1] == self._instr[target + 1][1]:
@@ -539,7 +615,7 @@ class Circuit:
         else:
             self._commute_with_following(target)
 
-    def _commute_with_Rzz(self, target):
+    def _commute_with_Rzz(self, target: int):
         assert self._instr[target][0] in ["XC", "ZC"]
         assert self._instr[target + 1][0] == "Rzz"
         if self._instr[target][0] == "XC":
@@ -550,7 +626,7 @@ class Circuit:
                 self._M[self._instr[target + 1][3]][4].extend(self._instr[target][2])
         self._commute_with_following(target)
 
-    def _commute_with_following(self, target):
+    def _commute_with_following(self, target: int):
         """Internal method to perform the commutation of
         two consecutive commands that commutes.
         commutes the target command with the following command.
@@ -564,7 +640,7 @@ class Circuit:
         self._instr.pop(target + 1)
         self._instr.insert(target, A)
 
-    def _find_byproduct_to_move(self, rev=False, skipnum=0):
+    def _find_byproduct_to_move(self, rev: bool = False, skipnum: int = 0):
         """Internal method for reordering commands
         Parameters
         ----------
@@ -602,6 +678,8 @@ class Circuit:
                 continue
             if self._instr[target + 1][0] == "CNOT":
                 target = self._commute_with_cnot(target)
+            elif self._instr[target + 1][0] == "SWAP":
+                target = self._commute_with_swap(target)
             elif self._instr[target + 1][0] == "H":
                 self._commute_with_H(target)
             elif self._instr[target + 1][0] == "S":
@@ -620,7 +698,7 @@ class Circuit:
             target += 1
 
     @classmethod
-    def _cnot_command(self, control_node, target_node, ancilla):
+    def _cnot_command(self, control_node: int, target_node: int, ancilla: Sequence[int]):
         """MBQC commands for CNOT gate
 
         Parameters
@@ -654,7 +732,7 @@ class Circuit:
         return control_node, ancilla[1], seq
 
     @classmethod
-    def _h_command(self, input_node, ancilla):
+    def _h_command(self, input_node: int, ancilla: int):
         """MBQC commands for Hadamard gate
 
         Parameters
@@ -678,7 +756,7 @@ class Circuit:
         return ancilla, seq
 
     @classmethod
-    def _s_command(self, input_node, ancilla):
+    def _s_command(self, input_node: int, ancilla: Sequence[int]):
         """MBQC commands for S gate
 
         Parameters
@@ -706,7 +784,7 @@ class Circuit:
         return ancilla[1], seq
 
     @classmethod
-    def _x_command(self, input_node, ancilla):
+    def _x_command(self, input_node: int, ancilla: Sequence[int]):
         """MBQC commands for Pauli X gate
 
         Parameters
@@ -734,7 +812,7 @@ class Circuit:
         return ancilla[1], seq
 
     @classmethod
-    def _y_command(self, input_node, ancilla):
+    def _y_command(self, input_node: int, ancilla: Sequence[int]):
         """MBQC commands for Pauli Y gate
 
         Parameters
@@ -767,7 +845,7 @@ class Circuit:
         return ancilla[3], seq
 
     @classmethod
-    def _z_command(self, input_node, ancilla):
+    def _z_command(self, input_node: int, ancilla: Sequence[int]):
         """MBQC commands for Pauli Z gate
 
         Parameters
@@ -795,7 +873,7 @@ class Circuit:
         return ancilla[1], seq
 
     @classmethod
-    def _rx_command(self, input_node, ancilla, angle):
+    def _rx_command(self, input_node: int, ancilla: Sequence[int], angle: float):
         """MBQC commands for X rotation gate
 
         Parameters
@@ -825,7 +903,7 @@ class Circuit:
         return ancilla[1], seq
 
     @classmethod
-    def _ry_command(self, input_node, ancilla, angle):
+    def _ry_command(self, input_node: int, ancilla: Sequence[int], angle: float):
         """MBQC commands for Y rotation gate
 
         Parameters
@@ -860,7 +938,7 @@ class Circuit:
         return ancilla[3], seq
 
     @classmethod
-    def _rz_command(self, input_node, ancilla, angle):
+    def _rz_command(self, input_node: int, ancilla: Sequence[int], angle: float):
         """MBQC commands for Z rotation gate
 
         Parameters
@@ -890,7 +968,7 @@ class Circuit:
         return ancilla[1], seq
 
     @classmethod
-    def _rz_command_opt(self, input_node, ancilla, angle):
+    def _rz_command_opt(self, input_node: int, ancilla: int, angle: float):
         """optimized MBQC commands for Z rotation gate
 
         Parameters
@@ -916,7 +994,7 @@ class Circuit:
         return input_node, seq
 
     @classmethod
-    def _rzz_command_opt(self, control_node, target_node, ancilla, angle):
+    def _rzz_command_opt(self, control_node: int, target_node: int, ancilla: int, angle: float):
         """Optimized MBQC commands for ZZ-rotation gate
 
         Parameters
@@ -946,17 +1024,204 @@ class Circuit:
         return control_node, target_node, seq
 
     @classmethod
-    def _sort_outputs(self, pattern, output_nodes):
-        """Sort the node indices of ouput qubits.
+    def _ccx_command(self, control_node1: int, control_node2: int, target_node: int, ancilla: Sequence[int]):
+        """MBQC commands for CCX gate
 
         Parameters
         ---------
-        input_node : int
-            input node index
-        ancilla : list of two ints
+        control_node1 : int
+            first control node on graph
+        control_node2 : int
+            second control node on graph
+        target_node : int
+            target node on graph
+        ancilla : list of int
             ancilla node indices to be added to graph
-        angle : float
-            measurement angle in radian
+
+        Returns
+        ---------
+        control_out1 : int
+            first control node on graph after the gate
+        control_out2 : int
+            second control node on graph after the gate
+        target_out : int
+            target node on graph after the gate
+        commands : list
+            list of MBQC commands
+        """
+        assert len(ancilla) == 18
+        seq = [["N", ancilla[i]] for i in range(18)]  # assign new qubit labels
+        seq.append(["E", (target_node, ancilla[0])])
+        seq.append(["E", (ancilla[0], ancilla[1])])
+        seq.append(["E", (ancilla[1], ancilla[2])])
+        seq.append(["E", (ancilla[1], control_node2)])
+        seq.append(["E", (control_node1, ancilla[14])])
+        seq.append(["E", (ancilla[2], ancilla[3])])
+        seq.append(["E", (ancilla[14], ancilla[4])])
+        seq.append(["E", (ancilla[3], ancilla[5])])
+        seq.append(["E", (ancilla[3], ancilla[4])])
+        seq.append(["E", (ancilla[5], ancilla[6])])
+        seq.append(["E", (control_node2, ancilla[6])])
+        seq.append(["E", (control_node2, ancilla[9])])
+        seq.append(["E", (ancilla[6], ancilla[7])])
+        seq.append(["E", (ancilla[9], ancilla[4])])
+        seq.append(["E", (ancilla[9], ancilla[10])])
+        seq.append(["E", (ancilla[7], ancilla[8])])
+        seq.append(["E", (ancilla[10], ancilla[11])])
+        seq.append(["E", (ancilla[4], ancilla[8])])
+        seq.append(["E", (ancilla[4], ancilla[11])])
+        seq.append(["E", (ancilla[4], ancilla[16])])
+        seq.append(["E", (ancilla[8], ancilla[12])])
+        seq.append(["E", (ancilla[11], ancilla[15])])
+        seq.append(["E", (ancilla[12], ancilla[13])])
+        seq.append(["E", (ancilla[16], ancilla[17])])
+        seq.append(["M", target_node, "XY", -0.0, [], []])
+        seq.append(["M", ancilla[0], "XY", -0.0, [target_node], []])
+        seq.append(["M", ancilla[1], "XY", -0.0, [ancilla[0]], []])
+        seq.append(["M", control_node1, "XY", -0.0, [], []])
+        seq.append(["M", ancilla[2], "XY", -1.75, [ancilla[1], target_node], []])
+        seq.append(["M", ancilla[14], "XY", -0.0, [control_node1], []])
+        seq.append(["M", ancilla[3], "XY", -0.0, [ancilla[2], ancilla[0]], []])
+        seq.append(["M", ancilla[5], "XY", -0.25, [ancilla[3], ancilla[1], ancilla[14], target_node], []])
+        seq.append(["M", control_node2, "XY", -0.25, [], []])
+        seq.append(["M", ancilla[6], "XY", -0.0, [ancilla[5], ancilla[2], ancilla[0]], []])
+        seq.append(["M", ancilla[9], "XY", -0.0, [control_node2, ancilla[0], ancilla[5], ancilla[2], ancilla[0]], []])
+        seq.append(["M", ancilla[7], "XY", -1.75, [ancilla[6], ancilla[3], ancilla[1], ancilla[14], target_node], []])
+        seq.append(["M", ancilla[10], "XY", -1.75, [ancilla[9], ancilla[14]], []])
+        seq.append(["M", ancilla[4], "XY", -0.25, [ancilla[14]], []])
+        seq.append(["M", ancilla[8], "XY", -0.0, [ancilla[7], ancilla[5], ancilla[2], ancilla[0]], []])
+        seq.append(
+            [
+                "M",
+                ancilla[11],
+                "XY",
+                -0.0,
+                [ancilla[10], control_node2, ancilla[0], ancilla[5], ancilla[2], ancilla[0]],
+                [],
+            ]
+        )
+        seq.append(
+            [
+                "M",
+                ancilla[12],
+                "XY",
+                -0.25,
+                [ancilla[8], ancilla[14], ancilla[6], ancilla[3], ancilla[1], ancilla[14], target_node],
+                [],
+            ]
+        )
+        seq.append(
+            [
+                "M",
+                ancilla[16],
+                "XY",
+                -0.0,
+                [
+                    ancilla[4],
+                    control_node1,
+                    ancilla[2],
+                    control_node2,
+                    ancilla[7],
+                    ancilla[10],
+                    ancilla[0],
+                    ancilla[0],
+                    ancilla[5],
+                    ancilla[2],
+                    ancilla[0],
+                    ancilla[5],
+                    ancilla[2],
+                    ancilla[0],
+                    control_node2,
+                    ancilla[0],
+                    ancilla[5],
+                    ancilla[2],
+                    ancilla[0],
+                ],
+                [],
+            ]
+        )
+        seq.append(["X", ancilla[17], [ancilla[14], ancilla[16]]])
+        seq.append(["X", ancilla[15], [ancilla[9], ancilla[11]]])
+        seq.append(["X", ancilla[13], [ancilla[0], ancilla[2], ancilla[5], ancilla[7], ancilla[12]]])
+        seq.append(["Z", ancilla[17], [ancilla[4], ancilla[5], ancilla[7], ancilla[10], control_node1]])
+        seq.append(["Z", ancilla[15], [control_node2, ancilla[2], ancilla[5], ancilla[10]]])
+        seq.append(["Z", ancilla[13], [ancilla[1], ancilla[3], ancilla[6], ancilla[8], target_node]])
+        return ancilla[17], ancilla[15], ancilla[13], seq
+
+    @classmethod
+    def _ccx_command_opt(self, control_node1: int, control_node2: int, target_node: int, ancilla: Sequence[int]):
+        """Optimized MBQC commands for CCX gate
+
+        Parameters
+        ---------
+        control_node1 : int
+            first control node on graph
+        control_node2 : int
+            second control node on graph
+        target_node : int
+            target node on graph
+        ancilla : list of int
+            ancilla node indices to be added to graph
+
+        Returns
+        ---------
+        control_out1 : int
+            first control node on graph after the gate
+        control_out2 : int
+            second control node on graph after the gate
+        target_out : int
+            target node on graph after the gate
+        commands : list
+            list of MBQC commands
+        """
+        assert len(ancilla) == 11
+        seq = [["N", ancilla[i]] for i in range(11)]  # assign new qubit labels
+        seq.append(["E", (control_node1, ancilla[8])])
+        seq.append(["E", (control_node2, ancilla[4])])
+        seq.append(["E", (control_node2, ancilla[5])])
+        seq.append(["E", (control_node2, ancilla[2])])
+        seq.append(["E", (control_node2, ancilla[0])])
+        seq.append(["E", (target_node, ancilla[6])])
+        seq.append(["E", (ancilla[0], ancilla[6])])
+        seq.append(["E", (ancilla[1], ancilla[10])])
+        seq.append(["E", (ancilla[2], ancilla[10])])
+        seq.append(["E", (ancilla[2], ancilla[6])])
+        seq.append(["E", (ancilla[3], ancilla[6])])
+        seq.append(["E", (ancilla[3], ancilla[10])])
+        seq.append(["E", (ancilla[4], ancilla[10])])
+        seq.append(["E", (ancilla[5], ancilla[9])])
+        seq.append(["E", (ancilla[6], ancilla[7])])
+        seq.append(["E", (ancilla[8], ancilla[10])])
+        seq.append(["M", target_node, "XY", -0.0, [], []])
+        seq.append(["M", control_node1, "XY", -0.0, [], []])
+        seq.append(["M", ancilla[0], "XY", -1.75, [target_node], [], 6])
+        seq.append(["M", ancilla[8], "XY", -0.0, [control_node1], []])
+        seq.append(["M", ancilla[2], "XY", -0.25, [target_node, ancilla[8]], [], 6])
+        seq.append(["M", control_node2, "XY", -0.25, [], []])
+        seq.append(["M", ancilla[3], "XY", -1.75, [ancilla[8], target_node], [], 6])
+        seq.append(["M", ancilla[4], "XY", -1.75, [ancilla[8]], [], 6])
+        seq.append(["M", ancilla[1], "XY", -0.25, [ancilla[8]], [], 6])
+        seq.append(["M", ancilla[5], "XY", -0.0, [control_node2, ancilla[0], ancilla[2], ancilla[4]], []])
+        seq.append(["M", ancilla[6], "XY", -0.25, [target_node], []])
+        seq.append(["X", ancilla[10], [ancilla[8]]])
+        seq.append(["X", ancilla[9], [ancilla[5]]])
+        seq.append(["X", ancilla[7], [ancilla[0], ancilla[2], ancilla[3], ancilla[6]]])
+        seq.append(["Z", ancilla[10], [control_node1, ancilla[1], ancilla[2], ancilla[3], ancilla[4]]])
+        seq.append(["Z", ancilla[9], [control_node2, ancilla[0], ancilla[2], ancilla[4]]])
+        seq.append(["Z", ancilla[7], [target_node]])
+
+        return ancilla[10], ancilla[9], ancilla[7], seq
+
+    @classmethod
+    def _sort_outputs(self, pattern: Pattern, output_nodes: Sequence[int]):
+        """Sort the node indices of ouput qubits.
+
+        Parameters
+        ---------
+        pattern : :meth:`~graphix.pattern.Pattern`
+            pattern object
+        output_nodes : list of int
+            output node indices
 
         Returns
         ---------
@@ -979,12 +1244,16 @@ class Circuit:
             elif pattern.seq[i][1] in old_out:
                 pattern.seq[i][1] = output_nodes[old_out.index(pattern.seq[i][1])]
 
-    def simulate_statevector(self, input_state=None):
+    def simulate_statevector(self, input_state: Optional[Statevec] = None):
         """Run statevector simultion of the gate sequence, using graphix.Statevec
 
+        Parameters
+        ----------
+        input_state : :meth:`~graphix.Statevec`
+
         Returns
         -------
-        stete : graphix.Statevec
+        stete : :meth:`~graphix.Statevec`
             output state of the statevector simulation.
         """
 
@@ -996,6 +1265,8 @@ class Circuit:
         for i in range(len(self.instruction)):
             if self.instruction[i][0] == "CNOT":
                 state.CNOT((self.instruction[i][1][0], self.instruction[i][1][1]))
+            elif self.instruction[i][0] == "SWAP":
+                state.swap((self.instruction[i][1][0], self.instruction[i][1][1]))
             elif self.instruction[i][0] == "I":
                 pass
             elif self.instruction[i][0] == "S":
@@ -1016,5 +1287,9 @@ class Circuit:
                 state.evolve_single(Ops.Rz(self.instruction[i][2]), self.instruction[i][1])
             elif self.instruction[i][0] == "Rzz":
                 state.evolve(Ops.Rzz(self.instruction[i][2]), [self.instruction[i][1][0], self.instruction[i][1][1]])
+            elif self.instruction[i][0] == "CCX":
+                state.evolve(Ops.ccx, [self.instruction[i][1][0], self.instruction[i][1][1], self.instruction[i][1][2]])
+            else:
+                raise ValueError(f"Unknown instruction: {self.instruction[i][0]}")
 
         return state


Test patch:
diff --git a/tests/random_circuit.py b/tests/random_circuit.py
index 2c2dc6a..71980e8 100644
--- a/tests/random_circuit.py
+++ b/tests/random_circuit.py
@@ -1,5 +1,5 @@
 import numpy as np
-
+from copy import deepcopy
 from graphix.transpiler import Circuit
 
 GLOBAL_SEED = None
@@ -71,7 +71,21 @@ def genpair(n_qubits, count, rng):
     return pairs
 
 
-def get_rand_circuit(nqubits, depth, use_rzz=False, seed=None):
+def gentriplet(n_qubits, count, rng):
+    triplets = []
+    for i in range(count):
+        choice = [j for j in range(n_qubits)]
+        x = rng.choice(choice)
+        choice.pop(x)
+        y = rng.choice(choice)
+        locy = np.where(y == np.array(deepcopy(choice)))[0][0]
+        choice.pop(locy)
+        z = rng.choice(choice)
+        triplets.append((x, y, z))
+    return triplets
+
+
+def get_rand_circuit(nqubits, depth, use_rzz=False, use_ccx=False, seed=None):
     rng = get_rng(seed)
     circuit = Circuit(nqubits)
     gate_choice = [0, 1, 2, 3, 4, 5, 6, 7]
@@ -79,8 +93,13 @@ def get_rand_circuit(nqubits, depth, use_rzz=False, seed=None):
         for j, k in genpair(nqubits, 2, rng):
             circuit.cnot(j, k)
         if use_rzz:
-            for j, k in genpair(nqubits, 1, rng):
+            for j, k in genpair(nqubits, 2, rng):
                 circuit.rzz(j, k, np.pi / 4)
+        if use_ccx:
+            for j, k, l in gentriplet(nqubits, 2, rng):
+                circuit.ccx(j, k, l)
+        for j, k in genpair(nqubits, 4, rng):
+            circuit.swap(j, k)
         for j in range(nqubits):
             k = rng.choice(gate_choice)
             if k == 0:
diff --git a/tests/test_tnsim.py b/tests/test_tnsim.py
index ddb4b39..290a05a 100644
--- a/tests/test_tnsim.py
+++ b/tests/test_tnsim.py
@@ -302,6 +302,21 @@ class TestTN(unittest.TestCase):
         value2 = tn_mbqc.expectation_value(random_op2, [0, 1])
         np.testing.assert_almost_equal(value1, value2)
 
+    def test_ccx(self):
+        nqubits = 4
+        depth = 6
+        for _ in range(10):
+            circuit = rc.get_rand_circuit(nqubits, depth)
+            circuit.ccx(0, 1, 2)
+            pattern = circuit.transpile()
+            pattern.minimize_space()
+            state = circuit.simulate_statevector()
+            tn_mbqc = pattern.simulate_pattern(backend="tensornetwork")
+            random_op3 = random_op(3)
+            value1 = state.expectation_value(random_op3, [0, 1, 2])
+            value2 = tn_mbqc.expectation_value(random_op3, [0, 1, 2])
+            np.testing.assert_almost_equal(value1, value2)
+
     def test_with_graphtrans(self):
         nqubits = 4
         depth = 6
diff --git a/tests/test_transpiler.py b/tests/test_transpiler.py
index 8a47290..0225cf3 100644
--- a/tests/test_transpiler.py
+++ b/tests/test_transpiler.py
@@ -93,8 +93,31 @@ class TestTranspiler_UnitGates(unittest.TestCase):
         state_mbqc = pattern.simulate_pattern()
         np.testing.assert_almost_equal(np.abs(np.dot(state_mbqc.flatten().conjugate(), state.flatten())), 1)
 
+    def test_ccx(self):
+        nqubits = 4
+        depth = 6
+        for _ in range(10):
+            circuit = rc.get_rand_circuit(nqubits, depth, use_ccx=True)
+            pattern = circuit.transpile()
+            pattern.minimize_space()
+            state = circuit.simulate_statevector()
+            state_mbqc = pattern.simulate_pattern()
+            np.testing.assert_almost_equal(np.abs(np.dot(state_mbqc.flatten().conjugate(), state.flatten())), 1)
+
 
 class TestTranspiler_Opt(unittest.TestCase):
+    def test_ccx_opt(self):
+        nqubits = 4
+        depth = 6
+        for _ in range(10):
+            circuit = rc.get_rand_circuit(nqubits, depth, use_ccx=True)
+            circuit.ccx(0, 1, 2)
+            pattern = circuit.transpile(opt=True)
+            pattern.minimize_space()
+            state = circuit.simulate_statevector()
+            state_mbqc = pattern.simulate_pattern()
+            np.testing.assert_almost_equal(np.abs(np.dot(state_mbqc.flatten().conjugate(), state.flatten())), 1)
+
     def test_transpile_opt(self):
         nqubits = 2
         depth = 1