[CQT-243] Implementation of global phase

opensquirrel/circuit.py

@@ -5,6 +5,7 @@
 
 from opensquirrel.ir import IR, AsmDeclaration, Gate
 from opensquirrel.passes.exporter import ExportFormat
+from opensquirrel.phase_map import PhaseMap
 
 if TYPE_CHECKING:
     from opensquirrel.passes.decomposer import Decomposer
@@ -43,6 +44,7 @@ def __init__(self, register_manager: RegisterManager, ir: IR) -> None:
         """Create a circuit object from a register manager and an IR."""
         self.register_manager = register_manager
         self.ir = ir
+        self.phase_map = PhaseMap(self.register_manager.qubit_register.register_size)
 
     def __repr__(self) -> str:
         """Write the circuit to a cQASM 3 string."""
@@ -102,7 +104,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 export(self, fmt: ExportFormat | None = None) -> Any:
         if fmt == ExportFormat.QUANTIFY_SCHEDULER:
@@ -139,7 +141,7 @@ def replace(self, gate: type[Gate], replacement_gates_function: Callable[..., li
         """
         from opensquirrel.passes.decomposer import general_decomposer
 
-        general_decomposer.replace(self.ir, gate, replacement_gates_function)
+        general_decomposer.replace(self, gate, replacement_gates_function)
 
     def validate(self, validator: Validator) -> None:
         """Generic validator pass. It applies the given validator to the circuit."""

opensquirrel/common.py

@@ -39,6 +39,7 @@ def are_matrices_equivalent_up_to_global_phase(
     Returns:
         Whether two matrices are equivalent up to a global phase.
     """
+
     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
     )
@@ -51,6 +52,24 @@ def are_matrices_equivalent_up_to_global_phase(
     return np.allclose(matrix_a, phase_difference * matrix_b, atol=ATOL)
 
 
+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 np.complex128(1)
+
+    return np.complex128(matrix_a[first_non_zero] / matrix_b[first_non_zero])
+
+
 def is_identity_matrix_up_to_a_global_phase(matrix: NDArray[np.complex128]) -> bool:
     """Checks whether matrix is an identity matrix up to a global phase.
 
@@ -60,3 +79,13 @@ def is_identity_matrix_up_to_a_global_phase(matrix: NDArray[np.complex128]) -> b
         Whether matrix is an identity matrix up to a global phase.
     """
     return are_matrices_equivalent_up_to_global_phase(matrix, np.eye(matrix.shape[0], dtype=np.complex128))
+
+
+def get_phase_angle(scalar: np.complex128) -> np.float64:
+    """Derives the Euler rotation angle from a scalar.
+    Args:
+        scalar: scalar to convert.
+    Returns:
+        Euler phase angle of scalar.
+    """
+    return np.float64(-1j * np.log(scalar))

opensquirrel/ir.py

@@ -275,6 +275,12 @@ def __init__(self, index: QubitLike) -> None:
             msg = "index must be a QubitLike"
             raise TypeError(msg)
 
+    def __eq__(self, other: Any) -> bool:
+        """Compare two qubits."""
+        if not isinstance(other, Qubit):
+            return False
+        return self.__hash__() == other.__hash__()
+
     def __hash__(self) -> int:
         return hash(str(self.__class__) + str(self.index))
 

opensquirrel/passes/decomposer/general_decomposer.py

@@ -2,14 +2,25 @@
 
 from abc import ABC, abstractmethod
 from collections.abc import Callable, Iterable
-from typing import Any
+from typing import TYPE_CHECKING, Any
+
+import numpy as np
 
 from opensquirrel.circuit_matrix_calculator import get_circuit_matrix
-from opensquirrel.common import are_matrices_equivalent_up_to_global_phase, is_identity_matrix_up_to_a_global_phase
+from opensquirrel.common import (
+    ATOL,
+    are_matrices_equivalent_up_to_global_phase,
+    calculate_phase_difference,
+    get_phase_angle,
+    is_identity_matrix_up_to_a_global_phase,
+)
 from opensquirrel.default_instructions import is_anonymous_gate
-from opensquirrel.ir import IR, Gate
+from opensquirrel.ir import Float, Gate, Rz
 from opensquirrel.reindexer import get_reindexed_circuit
 
+if TYPE_CHECKING:
+    from opensquirrel.circuit import Circuit
+
 
 class Decomposer(ABC):
     def __init__(self, **kwargs: Any) -> None: ...
@@ -19,13 +30,23 @@ 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
+        circuit: 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()]
     replacement_gates_qubit_indices = set()
     replaced_matrix = get_circuit_matrix(get_reindexed_circuit([gate], gate_qubit_indices))
+    qubits = gate.get_qubit_operands()
 
     if is_identity_matrix_up_to_a_global_phase(replaced_matrix):
-        return
+        return []
 
     for g in replacement_gates:
         replacement_gates_qubit_indices.update([q.index for q in g.get_qubit_operands()])
@@ -40,25 +61,40 @@ 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 circuit is not None:
+        phase_difference = calculate_phase_difference(replaced_matrix, replacement_matrix)
+        euler_phase = get_phase_angle(phase_difference)
+        for qubit in qubits:
+            circuit.phase_map.add_qubit_phase(qubit, euler_phase)
+
+        if len(qubits) > 1:
+            relative_phase = float(
+                np.real(circuit.phase_map.get_qubit_phase(qubits[0]) - circuit.phase_map.get_qubit_phase(qubits[1]))
+            )
+            if abs(relative_phase) > ATOL:
+                list(replacement_gates).append(Rz(qubits[0], Float(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)
 
 
@@ -73,8 +109,8 @@ def decompose(self, gate: Gate) -> list[Gate]:
         return self.replacement_gates_function(*gate.arguments)
 
 
-def replace(ir: IR, gate: type[Gate], replacement_gates_function: Callable[..., list[Gate]]) -> None:
+def replace(circuit: Circuit, gate: type[Gate], replacement_gates_function: Callable[..., list[Gate]]) -> None:
     """Does the same as decomposer, but only applies to a given gate."""
     generic_replacer = _GenericReplacer(gate, replacement_gates_function)
 
-    decompose(ir, generic_replacer)
+    decompose(circuit, generic_replacer)

opensquirrel/phase_map.py

@@ -0,0 +1,20 @@
+import numpy as np
+
+from opensquirrel.ir import Qubit, QubitLike
+
+
+class PhaseMap:
+    def __init__(self, qubit_register_size: int) -> None:
+        """Initialize a PhaseMap object."""
+
+        self.qubit_phase_map = np.zeros(qubit_register_size, dtype=np.float64)
+
+    def __contains__(self, qubit: QubitLike) -> bool:
+        """Checks if qubit is in the phase map."""
+        return qubit in self.qubit_phase_map
+
+    def add_qubit_phase(self, qubit: QubitLike, phase: np.float64) -> None:
+        self.qubit_phase_map[Qubit(qubit).index] += phase
+
+    def get_qubit_phase(self, qubit: QubitLike) -> np.float64:
+        return np.float64(self.qubit_phase_map[Qubit(qubit).index])

tests/integration/test_integration.py

@@ -123,9 +123,11 @@ def test_spin2plus_backend() -> None:
 X90 q[0]
 Rz(-1.5707964) q[0]
 b[0] = measure q[0]
-Rz(2.3561946) q[1]
+Rz(0.78539824) q[1]
 X90 q[1]
-Rz(-3.1415927) q[1]
+Rz(1.5707963) q[1]
+X90 q[1]
+Rz(1.5707963) q[1]
 b[2] = measure q[1]
 Rz(1.5707963) q[1]
 X90 q[1]

tests/test_phase_map.py

@@ -0,0 +1,77 @@
+from __future__ import annotations
+
+import math
+
+import numpy as np
+
+from opensquirrel import Circuit, CircuitBuilder
+from opensquirrel.passes.decomposer import CNOT2CZDecomposer, CNOTDecomposer, McKayDecomposer, XYXDecomposer
+from opensquirrel.passes.merger import SingleQubitGatesMerger
+
+
+def test_integration_global_phase() -> None:
+    circuit = 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
+    circuit.decompose(decomposer=CNOTDecomposer())
+
+    # Replace CNOT gates with CZ gates
+    circuit.decompose(decomposer=CNOT2CZDecomposer())
+    # Merge single-qubit gates and decompose with McKay decomposition.
+    circuit.merge(merger=SingleQubitGatesMerger())
+    circuit.decompose(decomposer=McKayDecomposer())
+    assert (
+        str(circuit)
+        == """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]
+CZ q[1], q[0]
+X90 q[2]
+Rz(1.3907963) q[2]
+X90 q[2]
+Rz(3.1415927) q[0]
+X90 q[0]
+Rz(0.051592654) q[0]
+X90 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_phase_map_circuit_builder() -> None:
+    circuit_builder = CircuitBuilder(4)
+    circuit_builder.H(0).H(1).H(2).H(3).Rz(0, math.pi / 8).Ry(1, math.pi / 4).Rx(2, 3 * math.pi / 8).Ry(
+        3, math.pi / 2
+    ).CNOT(1, 2).CNOT(0, 3).Rz(0, math.pi / 8).Ry(1, math.pi / 4).Rx(2, 3 * math.pi / 8).Ry(3, math.pi / 2)
+    circuit = circuit_builder.to_circuit()
+    circuit.merge(SingleQubitGatesMerger())
+    circuit.decompose(XYXDecomposer())
+    assert np.allclose(circuit.phase_map.qubit_phase_map, [math.pi / 2, math.pi / 2, -math.pi / 2, math.pi / 2])

tests/test_replacer.py

@@ -61,7 +61,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)
@@ -77,7 +77,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)