Add ZYZ and ABC decompositions for quantify-scheduler

opensquirrel/circuit.py

@@ -1,10 +1,12 @@
-from typing import Callable, Dict, List
+from typing import Callable, Dict
 
 import numpy as np
 
 import opensquirrel.parsing.antlr.squirrel_ir_from_string
 from opensquirrel import circuit_matrix_calculator, mckay_decomposer, merger, replacer, writer
 from opensquirrel.default_gates import default_gate_aliases, default_gate_set
+from opensquirrel.export import quantify_scheduler_exporter
+from opensquirrel.export_format import ExportFormat
 from opensquirrel.parsing.libqasm.libqasm_ir_creator import LibqasmIRCreator
 from opensquirrel.replacer import Decomposer
 from opensquirrel.squirrel_ir import Gate, SquirrelIR
@@ -119,3 +121,9 @@ def __repr__(self) -> str:
         """Write the circuit to a cQasm3 string."""
 
         return writer.squirrel_ir_to_string(self.squirrel_ir)
+
+    def export(self, format: ExportFormat):
+        if format == ExportFormat.QUANTIFY_SCHEDULER:
+            return quantify_scheduler_exporter.export(self.squirrel_ir)
+
+        raise Exception("Unknown export format")

opensquirrel/cnot_decomposer.py

@@ -0,0 +1,74 @@
+import math
+
+from opensquirrel.common import ATOL
+from opensquirrel.default_gates import cnot, ry, rz
+from opensquirrel.identity_filter import filter_out_identities
+from opensquirrel.replacer import Decomposer
+from opensquirrel.squirrel_ir import BlochSphereRotation, ControlledGate, Float, Gate
+from opensquirrel.zyz_decomposer import ZYZDecomposer, get_zyz_decomposition_angles
+
+
+class CNOTDecomposer(Decomposer):
+    """
+    Decomposes 2-qubit controlled unitary gates to CNOT + rz/ry.
+    Applying single-qubit gate fusion after this pass might be beneficial.
+
+    Source of the math: https://threeplusone.com/pubs/on_gates.pdf, chapter 7.5 "ABC decomposition"
+    """
+
+    @staticmethod
+    def decompose(g: Gate) -> [Gate]:
+        if not isinstance(g, ControlledGate):
+            # Do nothing:
+            # - BlochSphereRotation's are only single-qubit,
+            # - decomposing MatrixGate is currently not supported.
+            return [g]
+
+        if not isinstance(g.target_gate, BlochSphereRotation):
+            # Do nothing.
+            # ControlledGate's with 2+ control qubits are ignored.
+            return [g]
+
+        # Perform ZYZ decomposition on the target gate.
+        # This gives us an ABC decomposition (U = AXBXC, ABC = I) of the target gate.
+        # See https://threeplusone.com/pubs/on_gates.pdf
+        theta0, theta1, theta2 = get_zyz_decomposition_angles(g.target_gate.angle, g.target_gate.axis)
+        target_qubit = g.target_gate.qubit
+
+        # First try to see if we can get away with a single CNOT.
+        # FIXME: see https://github.com/QuTech-Delft/OpenSquirrel/issues/99 this could be extended, I believe.
+        if abs(abs(theta0 + theta2) % (2 * math.pi)) < ATOL and abs(abs(theta1 - math.pi) % (2 * math.pi)) < ATOL:
+            # g == rz(theta0) Y rz(theta2) == rz(theta0 - pi / 2) X rz(theta2 + pi / 2)
+            # theta0 + theta2 == 0
+
+            alpha0 = theta0 - math.pi / 2
+            alpha2 = theta2 + math.pi / 2
+
+            return filter_out_identities(
+                [
+                    rz(q=target_qubit, theta=Float(alpha2)),
+                    cnot(control=g.control_qubit, target=target_qubit),
+                    rz(q=target_qubit, theta=Float(alpha0)),
+                    rz(q=g.control_qubit, theta=Float(g.target_gate.phase - math.pi / 2)),
+                ]
+            )
+
+        A = [ry(q=target_qubit, theta=Float(theta1 / 2)), rz(q=target_qubit, theta=Float(theta2))]
+
+        B = [
+            rz(q=target_qubit, theta=Float(-(theta0 + theta2) / 2)),
+            ry(q=target_qubit, theta=Float(-theta1 / 2)),
+        ]
+
+        C = [
+            rz(q=target_qubit, theta=Float((theta0 - theta2) / 2)),
+        ]
+
+        return filter_out_identities(
+            C
+            + [cnot(control=g.control_qubit, target=target_qubit)]
+            + B
+            + [cnot(control=g.control_qubit, target=target_qubit)]
+            + A
+            + [rz(q=g.control_qubit, theta=Float(g.target_gate.phase))]
+        )

opensquirrel/default_gates.py

@@ -43,6 +43,11 @@ def z90(q: Qubit) -> Gate:
     return BlochSphereRotation(qubit=q, axis=(0, 0, 1), angle=math.pi / 2, phase=0)
 
 
+@named_gate
+def zm90(q: Qubit) -> Gate:
+    return BlochSphereRotation(qubit=q, axis=(0, 0, 1), angle=-math.pi / 2, phase=0)
+
+
 @named_gate
 def rx(q: Qubit, theta: Float) -> Gate:
     return BlochSphereRotation(qubit=q, axis=(1, 0, 0), angle=theta.value, phase=0)
@@ -111,5 +116,10 @@ def sqrt_swap(q1: Qubit, q2: Qubit) -> Gate:
     )
 
 
-default_gate_set = [h, x, x90, xm90, y, y90, z, z90, cz, cr, crk, cnot, rx, ry, rz, x, swap, sqrt_swap]
+@named_gate
+def ccz(control1: Qubit, control2: Qubit, target: Qubit) -> Gate:
+    return ControlledGate(control1, cz(control2, target))
+
+
+default_gate_set = [h, x, x90, xm90, y, y90, z, z90, zm90, cz, cr, crk, cnot, rx, ry, rz, x, swap, sqrt_swap, ccz]
 default_gate_aliases = {"X": x, "RX": rx, "RY": ry, "RZ": rz, "Hadamard": h, "H": h}

opensquirrel/export/quantify_scheduler_exporter.py

@@ -0,0 +1,90 @@
+import math
+
+from opensquirrel.common import ATOL
+from opensquirrel.default_gates import x, z
+from opensquirrel.squirrel_ir import (
+    BlochSphereRotation,
+    ControlledGate,
+    MatrixGate,
+    Qubit,
+    SquirrelIR,
+    SquirrelIRVisitor,
+)
+
+try:
+    import quantify_scheduler
+    import quantify_scheduler.operations.gate_library as quantify_scheduler_gates
+except Exception as e:
+    pass
+
+
+_unsupported_gates_exception = Exception(
+    "Cannot export circuit: it contains unsupported gates - decompose them to the "
+    "Quantify-scheduler gate set first (rxy, rz, cnot, cz)"
+)
+
+
+class _ScheduleCreator(SquirrelIRVisitor):
+    def _get_qubit_string(self, q: Qubit) -> str:
+        return f"{self.qubit_register_name}[{q.index}]"
+
+    def __init__(self, qubit_register_name: str):
+        self.qubit_register_name = qubit_register_name
+        self.schedule = quantify_scheduler.Schedule(f"Exported OpenSquirrel circuit")
+
+    def visit_bloch_sphere_rotation(self, g: BlochSphereRotation):
+        if abs(g.axis[2]) < ATOL:
+            # Rxy rotation.
+            theta: float = g.angle
+            phi: float = math.atan2(g.axis[1], g.axis[0])
+            self.schedule.add(quantify_scheduler_gates.Rxy(theta=theta, phi=phi, qubit=self._get_qubit_string(g.qubit)))
+            return
+
+        if abs(g.axis[0]) < ATOL and abs(g.axis[1]) < ATOL:
+            # Rz rotation.
+            self.schedule.add(quantify_scheduler_gates.Rz(theta=g.angle, qubit=self._get_qubit_string(g.qubit)))
+            return
+
+        raise _unsupported_gates_exception
+
+    def visit_matrix_gate(self, g: MatrixGate):
+        raise _unsupported_gates_exception
+
+    def visit_controlled_gate(self, g: ControlledGate):
+        if not isinstance(g.target_gate, BlochSphereRotation):
+            raise _unsupported_gates_exception
+
+        if g.target_gate == x(g.target_gate.qubit):
+            self.schedule.add(
+                quantify_scheduler_gates.CNOT(
+                    qC=self._get_qubit_string(g.control_qubit), qT=self._get_qubit_string(g.target_gate.qubit)
+                )
+            )
+            return
+
+        if g.target_gate == z(g.target_gate.qubit):
+            self.schedule.add(
+                quantify_scheduler_gates.CZ(
+                    qC=self._get_qubit_string(g.control_qubit), qT=self._get_qubit_string(g.target_gate.qubit)
+                )
+            )
+            return
+
+        raise _unsupported_gates_exception
+
+
+def export(squirrel_ir: SquirrelIR):
+    if "quantify_scheduler" not in globals():
+
+        class QuantifySchedulerNotInstalled:
+            def __getattr__(self, attr_name):
+                raise Exception("quantify-scheduler is not installed, or cannot be installed on your system")
+
+        global quantify_scheduler
+        quantify_scheduler = QuantifySchedulerNotInstalled()
+        global quantify_scheduler_gates
+        quantify_scheduler_gates = QuantifySchedulerNotInstalled()
+
+    schedule_creator = _ScheduleCreator(squirrel_ir.qubit_register_name)
+    squirrel_ir.accept(schedule_creator)
+    return schedule_creator.schedule

opensquirrel/export_format.py

@@ -0,0 +1,5 @@
+from enum import Enum
+
+
+class ExportFormat(Enum):
+    QUANTIFY_SCHEDULER = 0

opensquirrel/identity_filter.py

@@ -0,0 +1,5 @@
+from opensquirrel.squirrel_ir import Gate
+
+
+def filter_out_identities(l: [Gate]):
+    return [g for g in l if not g.is_identity()]

opensquirrel/squirrel_ir.py

@@ -92,7 +92,7 @@ def __init__(self, generator, arguments):
 
     @property
     def name(self) -> Optional[str]:
-        return self.generator.__name__ if self.generator else None
+        return self.generator.__name__ if self.generator else "<anonymous>"
 
     @property
     def is_anonymous(self) -> bool:
@@ -102,6 +102,10 @@ def is_anonymous(self) -> bool:
     def get_qubit_operands(self) -> List[Qubit]:
         raise NotImplementedError
 
+    @abstractmethod
+    def is_identity(self) -> bool:
+        raise NotImplementedError
+
 
 class BlochSphereRotation(Gate):
     generator: Optional[Callable[..., "BlochSphereRotation"]] = None
@@ -183,6 +187,9 @@ def accept(self, visitor: SquirrelIRVisitor):
     def get_qubit_operands(self) -> List[Qubit]:
         return self.operands
 
+    def is_identity(self) -> bool:
+        return np.allclose(self.matrix, np.eye(2 ** len(self.operands)))
+
 
 class ControlledGate(Gate):
     generator: Optional[Callable[..., "ControlledGate"]] = None
@@ -210,17 +217,29 @@ def accept(self, visitor: SquirrelIRVisitor):
     def get_qubit_operands(self) -> List[Qubit]:
         return [self.control_qubit] + self.target_gate.get_qubit_operands()
 
+    def is_identity(self) -> bool:
+        return self.target_gate.is_identity()
+
 
 def named_gate(gate_generator: Callable[..., Gate]) -> Callable[..., Gate]:
     @wraps(gate_generator)
     def wrapper(*args, **kwargs):
-        for i, par in enumerate(inspect.signature(gate_generator).parameters.values()):
+        result = gate_generator(*args, **kwargs)
+        result.generator = wrapper
+
+        all_args = []
+        arg_index = 0
+        for par in inspect.signature(gate_generator).parameters.values():
             if not issubclass(par.annotation, Expression):
                 raise TypeError("Gate argument types must be expressions")
 
-        result = gate_generator(*args, **kwargs)
-        result.generator = wrapper
-        result.arguments = args
+            if par.name in kwargs:
+                all_args.append(kwargs[par.name])
+            else:
+                all_args.append(args[arg_index])
+                arg_index += 1
+
+        result.arguments = tuple(all_args)
         return result
 
     return wrapper

opensquirrel/zyz_decomposer.py

@@ -0,0 +1,97 @@
+import math
+from typing import Tuple
+
+from opensquirrel.common import ATOL
+from opensquirrel.default_gates import ry, rz
+from opensquirrel.identity_filter import filter_out_identities
+from opensquirrel.replacer import Decomposer
+from opensquirrel.squirrel_ir import BlochSphereRotation, Float, Gate
+
+
+def get_zyz_decomposition_angles(alpha: float, axis: Tuple[float, float, float]):
+    """
+    Gives the angles used in the Z-Y-Z decomposition of the Bloch sphere rotation
+    characterized by a rotation around `axis` of angle `alpha`.
+
+    Parameters:
+      alpha: angle of the Bloch sphere rotation
+      axis: _normalized_ axis of the Bloch sphere rotation
+
+    Returns:
+      a triple (theta1, theta2, theta3) corresponding to the decomposition of the
+      arbitrary Bloch sphere rotation into U = rz(theta3) ry(theta2) rz(theta1)
+    """
+
+    nx, ny, nz = axis
+
+    assert abs(nx**2 + ny**2 + nz**2 - 1) < ATOL, "Axis needs to be normalized"
+
+    assert -math.pi + ATOL < alpha <= math.pi + ATOL, "Angle needs to be normalized"
+
+    if abs(alpha - math.pi) < ATOL:
+        # alpha == pi, math.tan(alpha / 2) is not defined.
+
+        if abs(nz) < ATOL:
+            theta2 = math.pi
+            p = 0
+            m = 2 * math.acos(ny)
+
+        else:
+            p = math.pi
+            theta2 = 2 * math.acos(nz)
+
+            if abs(nz - 1) < ATOL or abs(nz + 1) < ATOL:
+                m = p  # This can be anything, but setting m = p means theta3 == 0, which is better for gate count.
+            else:
+                m = 2 * math.acos(ny / math.sqrt(1 - nz**2))
+
+    else:
+        p = 2 * math.atan2(nz * math.sin(alpha / 2), math.cos(alpha / 2))
+
+        acos_argument = math.cos(alpha / 2) * math.sqrt(1 + (nz * math.tan(alpha / 2)) ** 2)
+
+        # This fixes float approximations like 1.0000000000002 which acos doesn't like.
+        acos_argument = max(min(acos_argument, 1.0), -1.0)
+
+        theta2 = 2 * math.acos(acos_argument)
+        theta2 = math.copysign(theta2, alpha)
+
+        if abs(math.sin(theta2 / 2)) < ATOL:
+            m = p  # This can be anything, but setting m = p means theta3 == 0, which is better for gate count.
+        else:
+            acos_argument = ny * math.sin(alpha / 2) / math.sin(theta2 / 2)
+
+            # This fixes float approximations like 1.0000000000002 which acos doesn't like.
+            acos_argument = max(min(acos_argument, 1.0), -1.0)
+
+            m = 2 * math.acos(acos_argument)
+
+    theta1 = (p + m) / 2
+
+    theta3 = p - theta1
+
+    return theta1, theta2, theta3
+
+
+class ZYZDecomposer(Decomposer):
+    @staticmethod
+    def decompose(g: Gate) -> [Gate]:
+        if not isinstance(g, BlochSphereRotation):
+            # Only decompose single-qubit gates.
+            return [g]
+
+        theta1, theta2, theta3 = get_zyz_decomposition_angles(g.angle, g.axis)
+
+        z1 = rz(g.qubit, Float(theta1))
+        y = ry(g.qubit, Float(theta2))
+        z2 = rz(g.qubit, Float(theta3))
+
+        # Note: written like this, the decomposition doesn't preserve the global phase, which is fine
+        # since the global phase is a physically irrelevant artifact of the mathematical
+        # model we use to describe the quantum system.
+
+        # Should we want to preserve it, we would need to use a raw BlochSphereRotation, which would then
+        # be an anonymous gate in the resulting decomposed circuit:
+        # z2 = BlochSphereRotation(qubit=g.qubit, angle=theta3, axis=(0, 0, 1), phase = g.phase)
+
+        return filter_out_identities([z1, y, z2])