Introduce general ABAdecomposer class

juanboschero · juanboschero · commit 918bc2a8712e · 2024-06-14T15:33:34.000+02:00
diff --git a/opensquirrel/decomposer/aba_decomposer.py b/opensquirrel/decomposer/aba_decomposer.py
@@ -0,0 +1,142 @@
+from __future__ import annotations
+
+import math
+from abc import ABC, abstractmethod
+from collections.abc import Callable, Iterable
+
+import numpy as np
+
+from opensquirrel.common import ATOL
+from opensquirrel.decomposer.general_decomposer import Decomposer
+from opensquirrel.default_gates import Rx, Ry, Rz
+from opensquirrel.ir import Axis, AxisLike, BlochSphereRotation, Float, Gate, Qubit
+from opensquirrel.utils.identity_filter import filter_out_identities
+
+
+class ABADecomposer(Decomposer, ABC):
+    @property
+    @abstractmethod
+    def ra(self) -> Callable[..., BlochSphereRotation]: ...
+
+    @property
+    @abstractmethod
+    def rb(self) -> Callable[..., BlochSphereRotation]: ...
+
+    _gate_list: list[str] = [Rx.__name__, Ry.__name__, Rz.__name__]
+
+    def __init__(self) -> None:
+        self.index_a = self._gate_list.index(self.ra.__dict__["__wrapped__"].__name__)
+        self.index_b = self._gate_list.index(self.rb.__dict__["__wrapped__"].__name__)
+
+    def get_decomposition_angles(self, alpha: float, axis: AxisLike) -> tuple[float, float, float]:
+        """
+        Gives the angles used in the A-B-A 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 = Ra(theta3) Rb(theta2) Ra(theta1)
+
+        """
+        axis = Axis(axis)
+
+        if not (-math.pi + ATOL < alpha <= math.pi + ATOL):
+            ValueError("Angle needs to be normalized")
+
+        if abs(alpha - math.pi) < ATOL:
+            # alpha == pi, math.tan(alpha / 2) is not defined.
+
+            p: float
+            if abs(axis[self.index_a]) < ATOL:
+                theta2 = math.pi
+                p = 0
+                m = 2 * math.acos(axis[self.index_b])
+
+            else:
+                p = math.pi
+                theta2 = 2 * math.acos(axis[self.index_a])
+
+                if abs(axis[self.index_a] - 1) < ATOL or abs(axis[self.index_a] + 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(axis[self.index_b] / math.sqrt(1 - axis[self.index_a] ** 2))
+
+        else:
+            p = 2 * math.atan2(axis[self.index_a] * math.sin(alpha / 2), math.cos(alpha / 2))
+
+            acos_argument = math.cos(alpha / 2) * math.sqrt(1 + (axis[self.index_a] * 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 = axis[self.index_b] * 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
+
+    def decompose(self, g: Gate) -> list[Gate]:
+        if not isinstance(g, BlochSphereRotation):
+            # Only decomposer single-qubit gates.
+            return [g]
+
+        theta1, theta2, theta3 = self.get_decomposition_angles(g.angle, g.axis)
+        a1 = self.ra.__dict__["__wrapped__"](g.qubit, Float(theta1))
+        b = self.rb.__dict__["__wrapped__"](g.qubit, Float(theta2))
+        a2 = self.ra.__dict__["__wrapped__"](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([a1, b, a2])
+
+
+class ZYZDecomposer(ABADecomposer):
+    ra = Rz
+    rb = Ry
+
+
+class XYXDecomposer(ABADecomposer):
+    ra = Rx
+    rb = Ry
+
+
+class YZYDecomposer(ABADecomposer):
+    ra = Ry
+    rb = Rz
+
+
+class XZXDecomposer(ABADecomposer):
+    ra = Rx
+    rb = Rz
+
+
+class YXYDecomposer(ABADecomposer):
+    ra = Ry
+    rb = Rx
+
+
+class ZXZDecomposer(ABADecomposer):
+    ra = Rz
+    rb = Rx
diff --git a/opensquirrel/decomposer/xyx_decomposer.py b/opensquirrel/decomposer/xyx_decomposer.py
@@ -80,7 +80,7 @@ def decompose(g: Gate) -> [Gate]:
             # Only decomposer single-qubit gates.
             return [g]
 
-        theta1, theta2, theta3 = get_xyx_decomposition_angles(g.angle, g.axis)
+        theta1, theta2, theta3 = XYXDecomposer().get_decomposition_angles(g.angle, g.axis)
 
         x1 = Rx(g.qubit, Float(theta1))
         y = Ry(g.qubit, Float(theta2))
diff --git a/opensquirrel/decomposer/zyz_decomposer.py b/opensquirrel/decomposer/zyz_decomposer.py
@@ -79,7 +79,7 @@ def decompose(self, g: Gate) -> list[Gate]:
             # Only decomposer single-qubit gates.
             return [g]
 
-        theta1, theta2, theta3 = get_zyz_decomposition_angles(g.angle, g.axis)
+        theta1, theta2, theta3 = ZYZDecomposer().get_decomposition_angles(g.angle, g.axis)
 
         z1 = Rz(g.qubit, Float(theta1))
         y = Ry(g.qubit, Float(theta2))
diff --git a/test/decomposer/test_xyx_decomposer.py b/test/decomposer/test_xyx_decomposer.py
@@ -0,0 +1,42 @@
+import math
+
+import pytest
+
+from opensquirrel.decomposer.aba_decomposer import XYXDecomposer
+from opensquirrel.default_gates import CNOT, CR, H, I, Rx, Ry, S, X, Y
+from opensquirrel.ir import BlochSphereRotation, Float, Gate, Qubit
+
+
+@pytest.fixture(name="decomposer")
+def decomposer_fixture() -> XYXDecomposer:
+    return XYXDecomposer()
+
+
+@pytest.mark.parametrize(
+    "gate, expected_result",
+    [
+        (CNOT(Qubit(0), Qubit(1)), [CNOT(Qubit(0), Qubit(1))]),
+        (CR(Qubit(2), Qubit(3), Float(2.123)), [CR(Qubit(2), Qubit(3), Float(2.123))]),
+        (I(Qubit(0)), []),
+        (
+            S(Qubit(0)),
+            [Rx(Qubit(0), Float(math.pi / 2)), Ry(Qubit(0), Float(math.pi / 2)), Rx(Qubit(0), Float(-math.pi / 2))],
+        ),
+        (Y(Qubit(0)), [Ry(Qubit(0), Float(math.pi))]),
+        (Ry(Qubit(0), Float(0.9)), [Ry(Qubit(0), Float(0.9))]),
+        (X(Qubit(0)), [Rx(Qubit(0), Float(math.pi))]),
+        (Rx(Qubit(0), Float(0.123)), [Rx(Qubit(0), Float(0.123))]),
+        (H(Qubit(0)), [Rx(Qubit(0), Float(math.pi)), Ry(Qubit(0), Float(math.pi / 2))]),
+        (
+            BlochSphereRotation(qubit=Qubit(0), angle=5.21, axis=(1, 2, 3), phase=0.324),
+            [
+                Rx(Qubit(0), Float(0.8251439260060653)),
+                Ry(Qubit(0), Float(-1.030183660156084)),
+                Rx(Qubit(0), Float(-1.140443520488592)),
+            ],
+        ),
+    ],
+    ids=["CNOT", "CR", "I", "S", "Y", "Ry", "X", "Rx", "H", "arbitrary"],
+)
+def test_xyx_decomposer(decomposer: XYXDecomposer, gate: Gate, expected_result: list[Gate]) -> None:
+    assert decomposer.decompose(gate) == expected_result
diff --git a/test/decomposer/test_xzx_decomposer.py b/test/decomposer/test_xzx_decomposer.py
@@ -0,0 +1,42 @@
+import math
+
+import pytest
+
+from opensquirrel.decomposer.aba_decomposer import XZXDecomposer
+from opensquirrel.default_gates import CNOT, CR, H, I, Rx, Rz, S, X, Z
+from opensquirrel.ir import BlochSphereRotation, Float, Gate, Qubit
+
+
+@pytest.fixture(name="decomposer")
+def decomposer_fixture() -> XZXDecomposer:
+    return XZXDecomposer()
+
+
+@pytest.mark.parametrize(
+    "gate, expected_result",
+    [
+        (CNOT(Qubit(0), Qubit(1)), [CNOT(Qubit(0), Qubit(1))]),
+        (CR(Qubit(2), Qubit(3), Float(2.123)), [CR(Qubit(2), Qubit(3), Float(2.123))]),
+        (I(Qubit(0)), []),
+        (S(Qubit(0)), [Rz(Qubit(0), Float(math.pi / 2))]),
+        (Z(Qubit(0)), [Rz(Qubit(0), Float(math.pi))]),
+        (Rz(Qubit(0), Float(0.9)), [Rz(Qubit(0), Float(0.9))]),
+        (X(Qubit(0)), [Rx(Qubit(0), Float(math.pi))]),
+        (Rx(Qubit(0), Float(0.123)), [Rx(Qubit(0), Float(0.123))]),
+        (
+            H(Qubit(0)),
+            [Rx(Qubit(0), Float(math.pi / 2)), Rz(Qubit(0), Float(math.pi / 2)), Rx(Qubit(0), Float(math.pi / 2))],
+        ),
+        (
+            BlochSphereRotation(qubit=Qubit(0), angle=5.21, axis=(1, 2, 3), phase=0.324),
+            [
+                Rx(Qubit(0), Float(0.43035280630630446)),
+                Rz(Qubit(0), Float(-1.030183660156084)),
+                Rx(Qubit(0), Float(-0.7456524007888308)),
+            ],
+        ),
+    ],
+    ids=["CNOT", "CR", "I", "S", "Y", "Ry", "X", "Rx", "H", "arbitrary"],
+)
+def test_xzx_decomposer(decomposer: XZXDecomposer, gate: Gate, expected_result: list[Gate]) -> None:
+    assert decomposer.decompose(gate) == expected_result
diff --git a/test/decomposer/test_yxy_decomposer.py b/test/decomposer/test_yxy_decomposer.py
@@ -0,0 +1,45 @@
+import math
+
+import pytest
+
+from opensquirrel.decomposer.aba_decomposer import YXYDecomposer
+from opensquirrel.default_gates import CNOT, CR, H, I, Rx, Ry, S, X, Y
+from opensquirrel.ir import BlochSphereRotation, Float, Gate, Qubit
+
+
+@pytest.fixture(name="decomposer")
+def decomposer_fixture() -> YXYDecomposer:
+    return YXYDecomposer()
+
+
+@pytest.mark.parametrize(
+    "gate, expected_result",
+    [
+        (CNOT(Qubit(0), Qubit(1)), [CNOT(Qubit(0), Qubit(1))]),
+        (CR(Qubit(2), Qubit(3), Float(2.123)), [CR(Qubit(2), Qubit(3), Float(2.123))]),
+        (I(Qubit(0)), []),
+        (
+            S(Qubit(0)),
+            [Ry(Qubit(0), Float(math.pi / 2)), Rx(Qubit(0), Float(math.pi / 2)), Ry(Qubit(0), Float(-math.pi / 2))],
+        ),
+        (Y(Qubit(0)), [Ry(Qubit(0), Float(math.pi))]),
+        (Ry(Qubit(0), Float(0.9)), [Ry(Qubit(0), Float(0.9))]),
+        (X(Qubit(0)), [Rx(Qubit(0), Float(math.pi))]),
+        (Rx(Qubit(0), Float(0.123)), [Rx(Qubit(0), Float(0.123))]),
+        (
+            H(Qubit(0)),
+            [Ry(Qubit(0), Float(math.pi / 4)), Rx(Qubit(0), Float(math.pi)), Ry(Qubit(0), Float(-math.pi / 4))],
+        ),
+        (
+            BlochSphereRotation(qubit=Qubit(0), angle=5.21, axis=(1, 2, 3), phase=0.324),
+            [
+                Ry(Qubit(0), Float(0.9412144817800217)),
+                Rx(Qubit(0), Float(-0.893533136099803)),
+                Ry(Qubit(0), Float(-1.5568770630164868)),
+            ],
+        ),
+    ],
+    ids=["CNOT", "CR", "I", "S", "Y", "Ry", "X", "Rx", "H", "arbitrary"],
+)
+def test_yxy_decomposer(decomposer: YXYDecomposer, gate: Gate, expected_result: list[Gate]) -> None:
+    assert decomposer.decompose(gate) == expected_result
diff --git a/test/decomposer/test_yzy_decomposer.py b/test/decomposer/test_yzy_decomposer.py
@@ -0,0 +1,52 @@
+import math
+
+import pytest
+
+from opensquirrel.decomposer.aba_decomposer import YZYDecomposer
+from opensquirrel.default_gates import CNOT, CR, H, I, Rx, Ry, Rz, S, X, Y
+from opensquirrel.ir import BlochSphereRotation, Float, Gate, Qubit
+
+
+@pytest.fixture(name="decomposer")
+def decomposer_fixture() -> YZYDecomposer:
+    return YZYDecomposer()
+
+
+@pytest.mark.parametrize(
+    "gate, expected_result",
+    [
+        (CNOT(Qubit(0), Qubit(1)), [CNOT(Qubit(0), Qubit(1))]),
+        (CR(Qubit(2), Qubit(3), Float(2.123)), [CR(Qubit(2), Qubit(3), Float(2.123))]),
+        (I(Qubit(0)), []),
+        (S(Qubit(0)), [Rz(Qubit(0), Float(math.pi / 2))]),
+        (Y(Qubit(0)), [Ry(Qubit(0), Float(math.pi))]),
+        (Ry(Qubit(0), Float(0.9)), [Ry(Qubit(0), Float(0.9))]),
+        (
+            X(Qubit(0)),
+            [Ry(Qubit(0), Float(math.pi / 2)), Rz(Qubit(0), Float(math.pi)), Ry(Qubit(0), Float(-math.pi / 2))],
+        ),
+        (
+            Rx(Qubit(0), Float(0.123)),
+            [
+                Ry(Qubit(0), Float(math.pi / 2)),
+                Rz(Qubit(0), Float(0.12300000000000022)),
+                Ry(Qubit(0), Float(-math.pi / 2)),
+            ],
+        ),
+        (
+            H(Qubit(0)),
+            [Ry(Qubit(0), Float(math.pi / 4)), Rz(Qubit(0), Float(math.pi)), Ry(Qubit(0), Float(-math.pi / 4))],
+        ),
+        (
+            BlochSphereRotation(qubit=Qubit(0), angle=5.21, axis=(1, 2, 3), phase=0.324),
+            [
+                Ry(Qubit(0), Float(0.013919263778408464)),
+                Rz(Qubit(0), Float(-0.893533136099803)),
+                Ry(Qubit(0), Float(-0.6295818450148737)),
+            ],
+        ),
+    ],
+    ids=["CNOT", "CR", "I", "S", "Y", "Ry", "X", "Rx", "H", "arbitrary"],
+)
+def test_yzy_decomposer(decomposer: YZYDecomposer, gate: Gate, expected_result: list[Gate]) -> None:
+    assert decomposer.decompose(gate) == expected_result
diff --git a/test/decomposer/test_zxz_decomposer.py b/test/decomposer/test_zxz_decomposer.py
diff --git a/test/decomposer/test_zyz_decomposer.py b/test/decomposer/test_zyz_decomposer.py
diff --git a/test/test_xyx_decomposer.py b/test/test_xyx_decomposer.py
