[CQT-43] Take global phases into account when verifying decompositions

opensquirrel/circuit.py

@@ -2,11 +2,15 @@
 
 from collections.abc import Callable
 from typing import TYPE_CHECKING, Any
+from dataclasses import dataclass
+
+import numpy as np
 
 from opensquirrel.passes.exporter.export_format import ExportFormat
 
 if TYPE_CHECKING:
-    from opensquirrel.ir import IR, Gate
+    from numpy.typing import ArrayLike, NDArray
+    from opensquirrel.ir import IR, Gate, QubitLike
     from opensquirrel.passes.decomposer import Decomposer
     from opensquirrel.passes.mapper import Mapper
     from opensquirrel.register_manager import RegisterManager
@@ -99,7 +103,7 @@ def decompose(self, decomposer: Decomposer) -> None:
         """
         from opensquirrel.passes.decomposer import general_decomposer
 
-        general_decomposer.decompose(self.ir, decomposer)
+        general_decomposer.decompose(self, decomposer)
 
     def map(self, mapper: Mapper) -> None:
         """Generic qubit mapper pass.
@@ -116,7 +120,7 @@ def replace(self, gate_generator: Callable[..., Gate], f: Callable[..., list[Gat
         """
         from opensquirrel.passes.decomposer import general_decomposer
 
-        general_decomposer.replace(self.ir, gate_generator, f)
+        general_decomposer.replace(self, gate_generator, f)
 
     def export(self, fmt: ExportFormat | None = None) -> Any:
         if fmt == ExportFormat.QUANTIFY_SCHEDULER:
@@ -129,3 +133,31 @@ def export(self, fmt: ExportFormat | None = None) -> Any:
             return cqasmv1_exporter.export(self)
         msg = "unknown exporter format"
         raise ValueError(msg)
+
+
+@dataclass(init=False)
+class PhaseMap(Circuit):
+
+    def __init__(self, phase_map: ArrayLike[np.complex128]):
+        """Initialize a PhaseMap object."""
+        self.phase_map = phase_map
+
+    def __contains__(self, qubit: QubitLike) -> bool:
+        """Checks if qubit is in the phase map."""
+        if qubit in self.phase_map:
+            return True
+
+        return False
+
+    def set_phase_map(self, phase_map: NDArray[np.complex128]) -> None:
+        self.qubit_phase_map = phase_map
+
+    def add_qubit_phase(self, qubit: QubitLike, phase: np.complex128) -> None:
+        from opensquirrel.ir import Qubit
+
+        self.qubit_phase_map[Qubit(qubit).index] += phase
+
+    def get_qubit_phase(self, qubit: QubitLike) -> np.complex128:
+        from opensquirrel.ir import Qubit
+
+        return np.complex128(self.qubit_phase_map[Qubit(qubit).index])

opensquirrel/common.py

@@ -39,13 +39,39 @@ def are_matrices_equivalent_up_to_global_phase(
     Returns:
         Whether two matrices are equivalent up to a global phase.
     """
+    phase_difference = calculate_phase_difference(matrix_a, matrix_b)
+
+    if not phase_difference:
+        return False
+
+    return np.allclose(matrix_a, phase_difference * matrix_b)
+
+
+def calculate_phase_difference(matrix_a: NDArray[np.complex128], matrix_b: NDArray[np.complex128]) -> np.complex128:
+    """Calculates the phase difference between two matrices.
+
+    Args:
+        matrix_a: first matrix.
+        matrix_b: second matrix.
+
+    Returns:
+        The phase difference between the two matrices.
+    """
     first_non_zero = next(
         (i, j) for i in range(matrix_a.shape[0]) for j in range(matrix_a.shape[1]) if abs(matrix_a[i, j]) > ATOL
     )
 
     if abs(matrix_b[first_non_zero]) < ATOL:
-        return False
+        return np.complex128(1)
 
-    phase_difference = matrix_a[first_non_zero] / matrix_b[first_non_zero]
+    return np.complex128(matrix_a[first_non_zero] / matrix_b[first_non_zero])
 
-    return np.allclose(matrix_a, phase_difference * matrix_b)
+
+def to_euler_form(scalar: np.complex128) -> np.complex128:
+    """ " Derives the Euler rotation angle from a scalar.
+    Args:
+        scalar: scalar to convert.
+    Returns:
+        Euler phase angle of scalar.
+    """
+    return np.complex128(-1j * np.log(scalar))

opensquirrel/ir.py

@@ -187,6 +187,12 @@ def __init__(self, index: QubitLike) -> None:
     def __hash__(self) -> int:
         return hash(str(self.__class__) + str(self.index))
 
+    def __eq__(self, other: Any) -> bool:
+        """Compare two qubits."""
+        if not isinstance(other, Qubit):
+            return False
+        return self.__hash__() == other.__hash__()
+
     def __repr__(self) -> str:
         return f"Qubit[{self.index}]"
 

opensquirrel/passes/decomposer/general_decomposer.py

@@ -2,21 +2,44 @@
 
 from abc import ABC, abstractmethod
 from collections.abc import Callable, Iterable
+from typing import TYPE_CHECKING
+
+import numpy as np
 
 from opensquirrel.circuit_matrix_calculator import get_circuit_matrix
-from opensquirrel.common import are_matrices_equivalent_up_to_global_phase
-from opensquirrel.ir import IR, Gate
+from opensquirrel.common import (
+    ATOL,
+    are_matrices_equivalent_up_to_global_phase,
+    calculate_phase_difference,
+    to_euler_form,
+)
+from opensquirrel.default_gates import Rz
+from opensquirrel.ir import Float, Gate
 from opensquirrel.reindexer import get_reindexed_circuit
 
+if TYPE_CHECKING:
+    from opensquirrel.circuit import Circuit
+
 
 class Decomposer(ABC):
     @abstractmethod
     def decompose(self, gate: Gate) -> list[Gate]:
         raise NotImplementedError()
 
 
-def check_gate_replacement(gate: Gate, replacement_gates: Iterable[Gate]) -> None:
+def check_gate_replacement(gate: Gate, replacement_gates: Iterable[Gate], circuit: Circuit | None = None) -> list[Gate]:
+    """
+    Verifies the replacement gates against the given gate.
+    Args:
+        gate: original gate
+        replacement_gates: gates replacing the gate
+        qc: circuit to verify
+
+    Returns:
+        Returns verified list of replacement gates with possible correction.
+    """
     gate_qubit_indices = [q.index for q in gate.get_qubit_operands()]
+    qubit_list = gate.get_qubit_operands()
     replacement_gates_qubit_indices = set()
     for g in replacement_gates:
         replacement_gates_qubit_indices.update([q.index for q in g.get_qubit_operands()])
@@ -32,25 +55,37 @@ def check_gate_replacement(gate: Gate, replacement_gates: Iterable[Gate]) -> Non
         msg = f"replacement for gate {gate.name} does not preserve the quantum state"
         raise ValueError(msg)
 
+    if qc is not None:
+        phase_difference = calculate_phase_difference(replaced_matrix, replacement_matrix)
+        euler_phase = to_euler_form(phase_difference)
+        [qc.add_qubit_phase(q, euler_phase) for q in gate.get_qubit_operands()]
+
+        if len(gate_qubit_indices) > 1:
+            relative_phase = float(np.real(qc.get_qubit_phase(qubit_list[1]) - qc.get_qubit_phase(qubit_list[0])))
+            if abs(relative_phase) > ATOL:
+                list(replacement_gates).append(Rz(gate.get_qubit_operands()[0], Float(-1 * relative_phase)))
+
+    return list(replacement_gates)
+
 
-def decompose(ir: IR, decomposer: Decomposer) -> None:
+def decompose(circuit: Circuit, decomposer: Decomposer) -> None:
     """Applies `decomposer` to every gate in the circuit, replacing each gate by the output of `decomposer`.
     When `decomposer` decides to not decomposer a gate, it needs to return a list with the intact gate as single
     element.
     """
     statement_index = 0
-    while statement_index < len(ir.statements):
-        statement = ir.statements[statement_index]
+    while statement_index < len(circuit.ir.statements):
+        statement = circuit.ir.statements[statement_index]
 
         if not isinstance(statement, Gate):
             statement_index += 1
             continue
 
         gate = statement
         replacement_gates: list[Gate] = decomposer.decompose(statement)
-        check_gate_replacement(gate, replacement_gates)
+        replacement_gates = check_gate_replacement(gate, replacement_gates, circuit)
 
-        ir.statements[statement_index : statement_index + 1] = replacement_gates
+        circuit.ir.statements[statement_index : statement_index + 1] = replacement_gates
         statement_index += len(replacement_gates)
 
 
@@ -66,8 +101,8 @@ def decompose(self, g: Gate) -> list[Gate]:
         return self.replacement_function(*arguments)
 
 
-def replace(ir: IR, gate_generator: Callable[..., Gate], f: Callable[..., list[Gate]]) -> None:
+def replace(circuit: Circuit, gate_generator: Callable[..., Gate], f: Callable[..., list[Gate]]) -> None:
     """Does the same as decomposer, but only applies to a given gate."""
     generic_replacer = _GenericReplacer(gate_generator, f)
 
-    decompose(ir, generic_replacer)
+    decompose(circuit, generic_replacer)

test/test_integration.py

@@ -113,6 +113,75 @@ def test_Spin2_backend() -> None:  # noqa: N802
     )
 
 
+def test_integration_global_phase() -> None:
+    qc = Circuit.from_string(
+        """
+        version 3.0
+
+        qubit[3] q
+
+        H q[0:2]
+        Ry(1.5789) q[0]
+        H q[0]
+        CNOT q[1], q[0]
+        Ry(3.09) q[0]
+        Ry(0.318) q[1]
+        Ry(0.18) q[2]
+        CNOT q[2], q[0]
+        """,
+    )
+
+    # Decompose 2-qubit gates to a decomposition where the 2-qubit interactions are captured by CNOT gates
+    qc.decompose(decomposer=CNOTDecomposer())
+
+    # Replace CNOT gates with CZ gates
+    qc.replace(
+        CNOT,
+        lambda control, target: [
+            H(target),
+            CZ(control, target),
+            H(target),
+        ],
+    )
+
+    # Merge single-qubit gates and decompose with McKay decomposition.
+    qc.merge_single_qubit_gates()
+    qc.decompose(decomposer=McKayDecomposer())
+
+    assert (
+        str(qc)
+        == """version 3.0
+
+qubit[3] q
+
+Rz(1.5707963) q[1]
+X90 q[1]
+Rz(1.5707963) q[1]
+Rz(3.1415927) q[0]
+X90 q[0]
+Rz(0.0081036221) q[0]
+X90 q[0]
+Rz(3.1415927) q[0]
+CZ q[1], q[0]
+X90 q[2]
+Rz(1.3907963) q[2]
+X90 q[2]
+X90 q[0]
+Rz(0.051592695) q[0]
+X90 q[0]
+Rz(3.1415927) q[0]
+CZ q[2], q[0]
+Rz(1.5707963) q[0]
+X90 q[0]
+Rz(1.5707963) q[0]
+Rz(3.1415927) q[1]
+X90 q[1]
+Rz(2.8235927) q[1]
+X90 q[1]
+"""
+    )
+
+
 def test_hectoqubit_backend() -> None:
     qc = Circuit.from_string(
         """

test/test_replacer.py

@@ -62,7 +62,7 @@ def decompose(self, g: Gate) -> list[Gate]:
                     return [I(g.qubit), g, I(g.qubit)]
                 return [g]
 
-        decompose(circuit.ir, decomposer=TestDecomposer())
+        decompose(circuit, decomposer=TestDecomposer())
 
         builder2 = CircuitBuilder(3)
         builder2.I(0)
@@ -78,7 +78,7 @@ def test_replace(self) -> None:
         builder1.H(0)
         circuit = builder1.to_circuit()
 
-        replace(circuit.ir, H, lambda q: [Y90(q), X(q)])
+        replace(circuit, H, lambda q: [Y90(q), X(q)])
 
         builder2 = CircuitBuilder(3)
         builder2.Y90(0)