[CQT-434] Implementation of transitions between different two-qubit gate semantics

opensquirrel/ir/expression.py

@@ -7,6 +7,7 @@
 import numpy as np
 from numpy.typing import ArrayLike, DTypeLike, NDArray
 
+from opensquirrel.common import ATOL
 from opensquirrel.ir.ir import IRNode, IRVisitor
 
 
@@ -220,7 +221,7 @@ def __eq__(self, other: Any) -> bool:
         """
         if not isinstance(other, self.__class__):
             return False
-        return np.array_equal(self, other)
+        return np.allclose(self, other, atol=ATOL)
 
     @overload
     def __getitem__(self, i: int, /) -> np.float64: ...

opensquirrel/ir/semantics/canonical_gate.py

@@ -1,14 +1,31 @@
-from typing import Any
+from __future__ import annotations
+
+from typing import TYPE_CHECKING, Any
 
 import numpy as np
 from numpy.typing import NDArray
 
-from opensquirrel.ir import AxisLike, IRVisitor
 from opensquirrel.ir.expression import BaseAxis
 from opensquirrel.ir.semantics.gate_semantic import GateSemantic
 
+if TYPE_CHECKING:
+    from opensquirrel.ir import AxisLike, IRVisitor
+    from opensquirrel.ir.semantics import BlochSphereRotation
+
 
 class CanonicalAxis(BaseAxis):
+    restrict: bool = True
+
+    @staticmethod
+    def in_weyl_chamber(tx: float, ty: float, tz: float) -> bool:
+        """Checks if the canonical axis is in the Weyl chamber.
+
+        Returns:
+            True if the canonical axis is in the Weyl chamber, False otherwise.
+
+        """
+        return 1 / 2 >= tx >= ty >= tz >= 0 or 1 / 2 >= (1 - tx) >= ty >= tz > 0
+
     @staticmethod
     def parse(axis: AxisLike) -> NDArray[np.float64]:
         """Parse and validate an `AxisLike`.
@@ -27,6 +44,7 @@ def parse(axis: AxisLike) -> NDArray[np.float64]:
             ValueError: If the axis cannot be flattened to length 3.
 
         """
+
         if isinstance(axis, CanonicalAxis):
             return axis.value
 
@@ -39,8 +57,9 @@ def parse(axis: AxisLike) -> NDArray[np.float64]:
         if len(axis) != 3:
             msg = f"axis has size {len(axis)!r}: requires an ArrayLike of length 3"
             raise ValueError(msg)
-
-        return CanonicalAxis.restrict_to_weyl_chamber(axis)
+        if CanonicalAxis.restrict:
+            return CanonicalAxis.restrict_to_weyl_chamber(axis)
+        return axis
 
     @staticmethod
     def restrict_to_weyl_chamber(axis: NDArray[np.float64]) -> NDArray[np.float64]:
@@ -92,8 +111,12 @@ def accept(self, visitor: IRVisitor) -> Any:
 
 
 class CanonicalGateSemantic(GateSemantic):
-    def __init__(self, axis: AxisLike) -> None:
+    def __init__(self, axis: AxisLike, rotations: list[BlochSphereRotation] | None = None) -> None:
         self.axis = CanonicalAxis(axis)
+        self.rotations = rotations
+        if self.rotations and len(self.rotations) != 4:
+            msg = f"invalid number of rotations, expected 4 but got {len(self.rotations)}"
+            raise ValueError(msg)
 
     def accept(self, visitor: IRVisitor) -> Any:
         """Accepts visitor and processes this IR node."""
@@ -106,7 +129,10 @@ def is_identity(self) -> bool:
             True if the canonical gate semantic represents an identity operation, False otherwise.
 
         """
-        return self.axis == CanonicalAxis((0, 0, 0))
+        if not self.rotations:
+            return self.axis == CanonicalAxis((0, 0, 0))
+        return self.axis == CanonicalAxis((0, 0, 0)) and all(rotation.is_identity() for rotation in self.rotations)
 
     def __repr__(self) -> str:
-        return f"CanonicalGateSemantic(axis={self.axis})"
+        rotation_str = f", rotations={self.rotations}" if self.rotations else ""
+        return f"CanonicalGateSemantic(axis={self.axis}{rotation_str})"

opensquirrel/ir/two_qubit_gate.py

@@ -5,6 +5,7 @@
 
 from opensquirrel.ir import Gate, IRVisitor, Qubit, QubitLike
 from opensquirrel.ir.semantics import CanonicalGateSemantic, ControlledGateSemantic, MatrixGateSemantic
+from opensquirrel.ir.semantics.bsr import bsr_from_matrix
 from opensquirrel.ir.semantics.gate_semantic import GateSemantic
 from opensquirrel.utils import get_matrix
 
@@ -31,18 +32,35 @@ def matrix(self) -> MatrixGateSemantic:
         if self._matrix:
             return self._matrix
 
-        if self.controlled:
+        if self._controlled:
             self._matrix = MatrixGateSemantic(get_matrix(self, 2))
             return self._matrix
 
-        if self.canonical:
-            from opensquirrel.utils.matrix_expander import can2
+        if self._canonical:
+            from opensquirrel.utils.matrix_expander import can1, can2
+
+            if self._canonical.rotations:
+                k1, k2, k3, k4 = (
+                    can1(rotation.axis, rotation.angle, rotation.phase) for rotation in self._canonical.rotations
+                )
+                return MatrixGateSemantic(np.kron(k3, k4) @ can2(self._canonical.axis) @ np.kron(k1, k2))
+            return MatrixGateSemantic(can2(self._canonical.axis))
 
-            return MatrixGateSemantic(can2(self.canonical.axis))
-        return MatrixGateSemantic(np.eye(4))
+        msg = f"invalid gate semantic: {self.gate_semantic}"
+        raise ValueError(msg)
 
     @cached_property
-    def canonical(self) -> CanonicalGateSemantic | None:
+    def canonical(self) -> CanonicalGateSemantic:
+        if not self._canonical:
+            from opensquirrel.utils.matrix_expander import canonical_decomposition
+
+            k1, k2, k3, k4, axis = canonical_decomposition(np.array(self.matrix))
+
+            bsr1 = bsr_from_matrix(k1)
+            bsr2 = bsr_from_matrix(k2)
+            bsr3 = bsr_from_matrix(k3)
+            bsr4 = bsr_from_matrix(k4)
+            self._canonical = CanonicalGateSemantic(axis, [bsr1, bsr2, bsr3, bsr4])
         return self._canonical
 
     @cached_property

opensquirrel/passes/mapper/qgym_mapper.py

@@ -132,8 +132,7 @@ def _ir_to_graph(ir: IR) -> nx.Graph:
         for statement in ir.statements:
             if not isinstance(statement, Instruction):
                 continue
-            instruction = cast("Instruction", statement)  # type: ignore[redundant-cast]
-            qubit_indices = instruction.qubit_indices
+            qubit_indices = statement.qubit_indices
             for q_index in qubit_indices:
                 interaction_graph.add_node(q_index)
             if len(qubit_indices) >= 2:

opensquirrel/reader/libqasm_parser.py

@@ -239,7 +239,7 @@ def _get_registry(
         registry = OrderedDict()
         for variable in filter(type_check, ast.variables):
             registry[variable.name] = register_cls(variable.typ.size, variable.name)
-        return registry or OrderedDict({register_cls.default_name: register_cls(0)})
+        return cast("Registry", registry or OrderedDict({register_cls.default_name: register_cls(0)}))
 
     def _create_register_manager(self, ast: Any) -> RegisterManager:
         qubit_registry = self._get_registry(ast, QubitRegister, LibQasmParser._is_qubit_type)

opensquirrel/utils/matrix_expander.py

@@ -1,6 +1,7 @@
 from __future__ import annotations
 
 import cmath
+import itertools
 import math
 from collections.abc import Iterable
 from math import pi
@@ -9,6 +10,7 @@
 import numpy as np
 from numpy.typing import NDArray
 
+from opensquirrel.common import ATOL
 from opensquirrel.ir import (
     Axis,
     AxisLike,
@@ -140,9 +142,7 @@ def _controlled_gate(self, gate: TwoQubitGate) -> NDArray[np.complex128]:
             msg = f"gate {gate!r} does not have a controlled gate semantic"
             raise ValueError(msg)
 
-        qubit_operands = list(reversed(gate.qubit_operands))
-
-        if any(q.index >= self.qubit_register_size for q in qubit_operands):
+        if any(q.index >= self.qubit_register_size for q in gate.qubit_operands):
             msg = "index out of range"
             raise IndexError(msg)
 
@@ -310,6 +310,139 @@ def can2(canonical_axis: AxisLike) -> NDArray[np.complex128]:
     )
 
 
+def nearest_kronecker_product(c: NDArray[np.complex128]) -> tuple[NDArray[np.complex128], NDArray[np.complex128]]:
+    if c.shape != (4, 4):
+        msg = f"c has to have the shape (4, 4), but has shape {c.shape} instead."
+        raise ValueError(msg)
+
+    c = c.reshape(2, 2, 2, 2).swapaxes(1, 2).reshape(4, 4)
+
+    u, sv, vh = np.linalg.svd(c)
+    a = np.sqrt(sv[0]) * u[:, 0].reshape(2, 2)
+    b = np.sqrt(sv[0]) * vh[0, :].reshape(2, 2)
+    return a, b
+
+
+def global_phase_and_su4(unitary: NDArray[np.complex128]) -> tuple[float, NDArray[np.complex128]]:
+    """Extract global phase so that det(u_su) = 1.  U = e^{i alpha} u_su."""
+    d = np.linalg.det(unitary)
+    alpha = np.angle(d) / 4.0
+    u_su = unitary * np.exp(-1j * alpha)
+    return alpha, u_su
+
+
+def _orthogonal_diagonalization(
+    matrix: NDArray[np.complex128],
+) -> tuple[NDArray[np.complex128], NDArray[np.complex128]]:
+    """
+    Adapted from: https://github.com/gecrooks/quantumflow/blob/master/quantumflow/decompositions.py#L324
+    """
+    rng = np.random.default_rng()
+    max_attempts = 16
+    for _ in range(max_attempts):
+        c = rng.uniform(0, 1)
+        m = c * matrix.real + (1 - c) * matrix.imag
+        _, eigvecs = np.linalg.eigh(m)
+        eigvecs = np.array(eigvecs, dtype=complex)
+        eigvals = np.diag(eigvecs.transpose() @ matrix @ eigvecs)
+
+        # Finish if we got a correct answer.
+        reconstructed = eigvecs @ np.diag(eigvals) @ eigvecs.transpose()
+        if np.allclose(matrix, reconstructed):
+            return eigvals, eigvecs
+
+    msg = "matrix is not orthogonally diagonalizable"
+    raise np.linalg.LinAlgError(msg)
+
+
+def lambdas_to_coords(lambdas: NDArray[np.complex128]) -> tuple[float, float, float]:
+    l1, l2, _, l4 = lambdas
+    c1 = np.real(1j * np.log(l1 * l2))
+    c2 = np.real(1j * np.log(l2 * l4))
+    c3 = np.real(1j * np.log(l1 * l4))
+    coords = np.asarray((c1, c2, c3)) / np.pi
+
+    coords[np.abs(coords - 1) < ATOL] = -1
+    if all(coords < 0):
+        coords += 1
+
+    if np.abs(coords[0] - coords[1]) < ATOL:
+        coords[1] = coords[0]
+
+    if np.abs(coords[1] - coords[2]) < ATOL:
+        coords[2] = coords[1]
+
+    if np.abs(coords[0] - coords[1] - 1 / 2) < ATOL:
+        coords[1] = coords[0] - 1 / 2
+
+    coords[np.abs(coords) < ATOL] = 0
+    return coords
+
+
+def constrain_to_weyl_chamber(
+    lambdas: NDArray[np.complex128],
+) -> tuple[tuple[float, float, float], NDArray[np.int8], NDArray[np.int8]]:
+    for permutation in itertools.permutations(range(4)):
+        for signs in ([1, 1, 1, 1], [1, 1, -1, -1], [-1, 1, -1, 1], [1, -1, -1, 1]):
+            signed_lambdas = lambdas * np.asarray(signs)
+            lambdas_perm = signed_lambdas[list(permutation)]
+
+            coords = lambdas_to_coords(lambdas_perm)
+
+            if CanonicalAxis.in_weyl_chamber(*coords):
+                return coords, np.asarray(signs), np.asarray(permutation)
+
+    msg = "could not find a permutation and signs to put lambdas in the Weyl chamber."
+    raise ValueError(msg)
+
+
+def canonical_decomposition(
+    unitary: NDArray[np.complex128],
+) -> tuple[
+    NDArray[np.complex128],
+    NDArray[np.complex128],
+    NDArray[np.complex128],
+    NDArray[np.complex128],
+    CanonicalAxis,
+]:
+    """
+    This implentation of the canonical decomposition is heavily based on:
+    https://github.com/gecrooks/quantumflow/blob/master/quantumflow/decompositions.py
+    """
+
+    _, u_su = global_phase_and_su4(unitary)
+
+    m = (1 / np.sqrt(2)) * np.array([[1, 0, 0, 1j], [0, 1j, 1, 0], [0, 1j, -1, 0], [1, 0, 0, -1j]], dtype=np.complex128)
+    m_dag = m.conj().T
+
+    u_m = m_dag @ u_su @ m
+    q = u_m.T @ u_m
+
+    eig_vals, eig_vecs = _orthogonal_diagonalization(q)
+    lambdas = np.sqrt(eig_vals)
+    det_d = np.prod(lambdas)
+    if det_d.real < 0:
+        lambdas[0] = -lambdas[0]
+
+    (tx, ty, tz), signs, perm = constrain_to_weyl_chamber(lambdas)
+    lambdas = (signs * lambdas)[perm]
+    o2 = (np.diag(signs) @ eig_vecs.T)[perm]
+    d = np.diag(lambdas)
+    o1 = u_m @ o2.T @ d.conj()
+
+    if np.linalg.det(o2).real < 0:
+        o2[0, :] *= -1
+        o1[:, 0] *= -1
+
+    q1 = m @ o1 @ m_dag
+    q2 = m @ o2 @ m_dag
+
+    k1, k2 = nearest_kronecker_product(q2)
+    k3, k4 = nearest_kronecker_product(q1)
+
+    return k1, k2, k3, k4, CanonicalAxis(tx, ty, tz)
+
+
 def get_matrix(gate: Gate, qubit_register_size: int) -> NDArray[np.complex128]:
     """Compute the unitary matrix corresponding to the gate applied to those qubit operands, taken
     among any number of qubits. This can be used for, e.g.,

pyproject.toml

@@ -47,6 +47,7 @@ dev = [
 export = [
     "pyqt5-qt5==5.15.2; sys_platform != 'darwin'",
     "quantify-scheduler==0.28.0; sys_platform != 'darwin'",
+    "qcodes<0.57.0; sys_platform != 'darwin'",
 ]
 qgym_mapper = [
     "qgym==0.3.1",

tests/ir/semantics/test_canonical_gate_semantic.py

@@ -4,7 +4,7 @@
 from numpy.typing import NDArray
 
 from opensquirrel.ir import GateSemantic
-from opensquirrel.ir.semantics import CanonicalAxis, CanonicalGateSemantic
+from opensquirrel.ir.semantics import BlochSphereRotation, CanonicalAxis, CanonicalGateSemantic
 
 
 class TestCanonicalAxis:
@@ -30,10 +30,34 @@ class TestCanonicalGateSemantic:
     def semantic(self) -> CanonicalGateSemantic:
         return CanonicalGateSemantic((0, 0, 0))
 
-    def test_eq(self, semantic: CanonicalGateSemantic) -> None:
-        assert semantic.is_identity()
+    @pytest.fixture
+    def semantic_with_rotations(self) -> CanonicalGateSemantic:
+        """Fixture for a CanonicalGateSemantic with rotations."""
+        rotations = [
+            BlochSphereRotation(axis=(1, 0, 0), angle=0.5, phase=0.1),
+            BlochSphereRotation(axis=(0, 1, 0), angle=1.0, phase=0.2),
+            BlochSphereRotation(axis=(0, 0, 1), angle=1.5, phase=0.3),
+            BlochSphereRotation(axis=(1, 0, 0), angle=0.75, phase=0.4),
+        ]
+        return CanonicalGateSemantic((0.25, 0.25, 0.25), rotations)
 
     def test_init(self, semantic: CanonicalGateSemantic) -> None:
         assert isinstance(semantic, GateSemantic)
         assert hasattr(semantic, "axis")
         assert isinstance(semantic.axis, CanonicalAxis)
+
+    def test_is_identity_with_zero_axis(self, semantic: CanonicalGateSemantic) -> None:
+        assert semantic.is_identity()
+
+    def test_is_identity_with_non_zero_axis(self, semantic_with_rotations: CanonicalGateSemantic) -> None:
+        assert not semantic_with_rotations.is_identity()
+
+    def test_rotations_attribute_list(self, semantic_with_rotations: CanonicalGateSemantic) -> None:
+        assert semantic_with_rotations.rotations is not None
+        assert len(semantic_with_rotations.rotations) == 4
+        assert all(isinstance(rot, BlochSphereRotation) for rot in semantic_with_rotations.rotations)
+
+    def test_invalid_number_of_rotations(self) -> None:
+        rotations = [BlochSphereRotation(axis=(1, 0, 0), angle=0.5, phase=0.1)]
+        with pytest.raises(ValueError, match="invalid number of rotations, expected 4 but got 1"):
+            CanonicalGateSemantic((0, 0, 0), rotations)