Merge branch 'main' into replace_npfloat128

Author: NoureldinYosri

qualtran/rotation_synthesis/README.md

@@ -13,9 +13,11 @@
 #  limitations under the License.
 
 from qualtran.rotation_synthesis.math_config import NumpyConfig, with_dps
-from qualtran.rotation_synthesis.protocols.clifford_t_synthesis import (
+from qualtran.rotation_synthesis.matrix import to_cirq, to_quirk, to_sequence
+from qualtran.rotation_synthesis.protocols import (
     diagonal_unitary_approx,
     fallback_protocol,
+    magnitude_approx,
     mixed_diagonal_protocol,
     mixed_fallback_protocol,
 )

qualtran/rotation_synthesis/__init__.py

@@ -15,12 +15,15 @@
 from __future__ import annotations
 
 import abc
-from typing import Optional, Sequence
+from typing import Optional, Sequence, Union
 
 import attrs
+import cirq
+import numpy as np
 
 import qualtran.rotation_synthesis._typing as rst
 import qualtran.rotation_synthesis.math_config as mc
+import qualtran.rotation_synthesis.matrix.clifford_t_repr as ctr
 import qualtran.rotation_synthesis.rings as rings
 from qualtran.rotation_synthesis.matrix import su2_ct
 from qualtran.rotation_synthesis.rings import zsqrt2
@@ -33,7 +36,7 @@ def expected_num_ts(self, config: mc.MathConfig) -> rst.Real: ...
 
     @abc.abstractmethod
     def diamond_norm_distance_to_rz(self, theta: rst.Real, config: mc.MathConfig) -> rst.Real:
-        r"""Returns the diamond norm distance to $Rz(2\theta)$."""
+        r"""Returns the diamond norm distance to $e^{i\theta Z}$."""
 
 
 @attrs.frozen
@@ -110,6 +113,79 @@ def from_sequence(cls, seq: Sequence[str], twirl: bool = False) -> UnitaryChanne
         n = sum(g.startswith("T") for g in seq)
         return UnitaryChannel(u.matrix[0, 0], u.matrix[1, 0], n, twirl)
 
+    def to_cirq(self, fmt: str = "xz", qs: Optional[Sequence[cirq.Qid]] = None) -> cirq.Circuit:
+        """Retruns a representation of the channel as a cirq circuit.
+
+        Args:
+            fmt: The gates to use (see the documentation of to_sequence).
+            qs: Optional qubits to operate on.
+        Returns:
+            A cirq circuit
+        Raises:
+            ValueError: If twirl=True
+        """
+        if self.twirl:
+            raise ValueError("to_cirq is not supported when twirl=True")
+        if qs:
+            (q,) = qs
+        else:
+            q = cirq.q(0)
+        return cirq.Circuit(ctr.to_cirq(self.to_matrix(), fmt, q))
+
+    def to_quirk(self, fmt: str = "xz") -> str:
+        """Retruns a quirk link representing the channel operation.
+
+        Args:
+            fmt: The gates to use (see the documentation of to_sequence).
+        Returns:
+            A quirk link.
+        Raises:
+            ValueError: If twirl=True
+        """
+        if self.twirl:
+            raise ValueError("to_quirk is not supported when twirl=True")
+        gates = ctr.to_quirk(self.to_matrix(), fmt)
+        cols = '[' + ','.join(f'[{g}]' for g in gates) + ']'
+        return "https://algassert.com/quirk#circuit={\"cols\":%s}" % cols
+
+    def diamond_norm_distance_to_unitary(
+        self, unitary: np.ndarray, config: mc.MathConfig
+    ) -> rst.Real:
+        r"""Returns the diamond norm distance between self and the given unitary.
+
+        From Theorem B.1 of arxiv:2203.10064, the diamond norm distance between two untiaries
+        $U, V$ is $|v_1 - v_0|$ where $v_i$ are the eigen values of $V^\dagger U$. Geometrically
+        this is the diameter of the smallest disc in the complex plane that contains both
+        eigenvalues.
+        """
+        # W = V^\dagger U
+        w = self.to_matrix().adjoint().numpy(config) @ unitary
+        # Compute the eigen values of W
+        a = config.one
+        b = -w[0, 0] - w[1, 1]
+        c = w[0, 0] * w[1, 1] - w[0, 1] * w[1, 0]
+        d = config.sqrt(b**2 - 4 * a * c)
+        eigv0, eigv1 = [(-b - d) / (2 * a), (-b + d) / (2 * a)]
+        # Compute the norm of the difference.
+        diameter_vec = eigv1 - eigv0
+        return config.sqrt(diameter_vec.real**2 + diameter_vec.imag**2)
+
+    @classmethod
+    def from_unitaries(
+        cls, *unitaries: Union[UnitaryChannel, su2_ct.SU2CliffordT]
+    ) -> UnitaryChannel:
+        if not unitaries:
+            raise ValueError('at least one unitary should be provided')
+
+        unitary = su2_ct.ISqrt2
+        for u in unitaries:
+            if isinstance(u, UnitaryChannel):
+                unitary = unitary @ u.to_matrix()
+            else:
+                unitary = unitary @ u
+        unitary = unitary.rescale()
+        return UnitaryChannel(unitary.matrix[0, 0], unitary.matrix[1, 0], unitary.num_t_gates())
+
 
 @attrs.frozen
 class ProjectiveChannel(Channel):
@@ -173,6 +249,71 @@ def diamond_norm_distance_to_rz(self, theta: rst.Real, config: mc.MathConfig) ->
             theta - self.failure_angle(config), config
         )
 
+    def to_cirq(self, fmt: str = "xz", qs: Optional[Sequence[cirq.Qid]] = None) -> cirq.Circuit:
+        """Retruns a representation of the channel as a cirq circuit.
+
+        Args:
+            fmt: The gates to use (see the documentation of to_sequence).
+            qs: Optional qubits to operate on.
+        Returns:
+            A cirq circuit
+        Raises:
+            ValueError: If the correction channel is not a UnitaryChannel.
+        """
+        if qs:
+            q0, q1 = qs
+        else:
+            q0, q1 = cirq.LineQubit.range(2)
+        correction = self.correction
+        if not isinstance(correction, UnitaryChannel):
+            raise ValueError('to_cirq does not support a non unitary correction')
+        return cirq.Circuit(
+            cirq.CNOT(q0, q1),
+            cirq.CircuitOperation(self.rotation.to_cirq(fmt, (q0,)).freeze()),
+            cirq.CNOT(q0, q1),
+            cirq.measure(q1, key='m'),
+            cirq.CircuitOperation(correction.to_cirq(fmt, (q0,)).freeze()).with_classical_controls(
+                'm'
+            ),
+        )
+
+    def to_quirk(self, fmt: str = "xz") -> str:
+        """Retruns a quirk link representing the channel operation.
+
+        Args:
+            fmt: The gates to use (see the documentation of to_sequence).
+        Returns:
+            A quirk link.
+        """
+        correction = self.correction
+        if not isinstance(correction, UnitaryChannel):
+            raise ValueError(f"to_quirk is not supported for correction of type {type(correction)}")
+        rot = ctr.to_quirk(self.rotation.to_matrix(), fmt)
+        cor = ctr.to_quirk(correction.to_matrix(), fmt)
+        first_row = []
+        second_row = []
+        # CNOT
+        first_row.append("\"•\"")
+        second_row.append("\"X\"")
+        # rotation
+        first_row.extend(rot)
+        second_row.extend("1" for _ in rot)
+        # CNOT
+        first_row.append("\"•\"")
+        second_row.append("\"X\"")
+        # measure
+        first_row.append("1")
+        second_row.append("\"Measure\"")
+        # correction
+        first_row.extend(cor)
+        second_row.extend("\"•\"" for _ in cor)
+        cols = (
+            '['
+            + ','.join(f'[{g1},{g2}]' for g1, g2 in zip(first_row, second_row, strict=True))
+            + ']'
+        )
+        return "https://algassert.com/quirk#circuit={\"cols\":%s}" % cols
+
 
 @attrs.frozen
 class ProbabilisticChannel(Channel):

qualtran/rotation_synthesis/channels/channel.py

@@ -12,11 +12,13 @@
 #  See the License for the specific language governing permissions and
 #  limitations under the License.
 
+import cirq
 import numpy as np
 import pytest
 
 import qualtran.rotation_synthesis.channels as ch
 import qualtran.rotation_synthesis.math_config as mc
+import qualtran.rotation_synthesis.matrix.clifford_t_repr as ctr
 import qualtran.rotation_synthesis.matrix.su2_ct as su2_ct
 
 
@@ -46,6 +48,10 @@ def test_unitary_from_sequence(gates):
 def test_diamond_distance_for_unitary(gates, theta, distance):
     c = ch.UnitaryChannel.from_sequence(gates)
     np.testing.assert_allclose(c.diamond_norm_distance_to_rz(theta, mc.NumpyConfig), distance)
+    u = np.zeros((2, 2), complex)
+    u[1, 1] = np.exp(-1j * theta)
+    u[0, 0] = u[1, 1].conjugate()
+    np.testing.assert_allclose(c.diamond_norm_distance_to_unitary(u, mc.NumpyConfig), distance)
 
 
 @pytest.mark.parametrize(
@@ -306,3 +312,38 @@ def test_diamond_distance_for_mixed_fallback(
     np.testing.assert_allclose(
         float(c.diamond_norm_distance_to_rz(theta, mc.NumpyConfig)), distance, atol=3e-5
     )
+
+
+@pytest.mark.parametrize(
+    "gates",
+    [["I", "Z", "I", "Tz"], ["I", "S", "Tz"], ["I", "S", "Tz"], ["I", "Z", "S", "Tz"], ["I", "S"]],
+)
+@pytest.mark.parametrize("fmt", ["xz", "xyz"])
+def test_unitary_to_cirq(gates, fmt):
+    u = ch.UnitaryChannel.from_sequence(gates)
+    assert u.to_cirq(fmt) == cirq.Circuit(ctr.to_cirq(u.to_matrix(), fmt))
+
+
+@pytest.mark.parametrize(
+    "gates1",
+    [["I", "Z", "I", "Tz"], ["I", "S", "Tz"], ["I", "S", "Tz"], ["I", "Z", "S", "Tz"], ["I", "S"]],
+)
+@pytest.mark.parametrize(
+    "gates2",
+    [["I", "Z", "I", "Tz"], ["I", "S", "Tz"], ["I", "S", "Tz"], ["I", "Z", "S", "Tz"], ["I", "S"]],
+)
+@pytest.mark.parametrize("fmt", ["xz", "xyz"])
+def test_fallback_to_cirq(gates1, gates2, fmt):
+    c = ch.ProjectiveChannel(
+        ch.UnitaryChannel.from_sequence(gates1), ch.UnitaryChannel.from_sequence(gates2)
+    )
+    q0, q1 = cirq.LineQubit.range(2)
+    assert c.to_cirq(fmt) == cirq.Circuit(
+        cirq.CNOT(q0, q1),
+        cirq.CircuitOperation(cirq.FrozenCircuit(ctr.to_cirq(c.rotation.to_matrix(), fmt, q0))),
+        cirq.CNOT(q0, q1),
+        cirq.measure(q1, key='m'),
+        cirq.CircuitOperation(
+            cirq.FrozenCircuit(ctr.to_cirq(c.correction.to_matrix(), fmt, q0))  # type: ignore[attr-defined]
+        ).with_classical_controls('m'),
+    )

qualtran/rotation_synthesis/channels/channel_test.py

@@ -213,7 +213,7 @@ def plot(self, ax: Optional[plt.Axes] = None, add_label: bool = True, **patch_ar
             fig, ax = plt.subplots(1)
 
         theta = float(self.tilt(mc.NumpyConfig))
-        e = self.rotate(-theta, mc.NumpyConfig)
+        e = self.rotate(theta, mc.NumpyConfig)
         w = 2 / np.sqrt(float(e.D[0, 0]))
         h = 2 / np.sqrt(float(e.D[1, 1]))
         c = self.center.astype(float).tolist()

qualtran/rotation_synthesis/lattice/geometry.py

@@ -40,6 +40,7 @@ class MathConfig:
     _ceil: Callable[[Real], int] = math.ceil
     _arctan2: Callable[[Real, Real], Real] = math.atan2
     _arcsin: Callable[[Real], Real] = math.asin
+    _arccos: Callable[[Real], Real] = math.acos
     _number: Callable[[Real], Real] = float
 
     def number(self, x) -> Real:
@@ -75,6 +76,9 @@ def __hash__(self) -> int:
     def arcsin(self, x: Real) -> Real:
         return self._arcsin(x)
 
+    def arccos(self, x: Real) -> Real:
+        return self._arccos(x)
+
 
 NumpyConfig = MathConfig(
     "numpy",
@@ -91,6 +95,7 @@ def arcsin(self, x: Real) -> Real:
     lambda x: int(np.ceil(x)),
     np.arctan2,
     np.arcsin,
+    np.arccos,
     np.longdouble,
 )
 
@@ -124,5 +129,6 @@ def with_dps(dps: int) -> MathConfig:
         lambda x: int(mpmath.ceil(x)),
         mpmath.atan2,
         mpmath.asin,
+        mpmath.acos,
         mpmath.mpf,
     )

qualtran/rotation_synthesis/math_config.py

@@ -14,8 +14,10 @@
 
 r"""A submodule for compiling $\mathbb{Z}[e^{i \pi/4}]$ matrices to Clifford+T as well as generating them."""
 
+from qualtran.rotation_synthesis.matrix.clifford_t_repr import to_cirq, to_quirk, to_sequence
 from qualtran.rotation_synthesis.matrix.generation import (
+    generate_cliffords,
     generate_rotations,
     generate_rotations_iter,
 )
-from qualtran.rotation_synthesis.matrix.su2_ct import generate_cliffords, SU2CliffordT
+from qualtran.rotation_synthesis.matrix.su2_ct import SU2CliffordT

qualtran/rotation_synthesis/matrix/__init__.py

@@ -0,0 +1,56 @@
+#  Copyright 2025 Google LLC
+#
+#  Licensed under the Apache License, Version 2.0 (the "License");
+#  you may not use this file except in compliance with the License.
+#  You may obtain a copy of the License at
+#
+#      https://www.apache.org/licenses/LICENSE-2.0
+#
+#  Unless required by applicable law or agreed to in writing, software
+#  distributed under the License is distributed on an "AS IS" BASIS,
+#  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+#  See the License for the specific language governing permissions and
+#  limitations under the License.
+
+import numpy as np
+
+import qualtran.rotation_synthesis._typing as rst
+import qualtran.rotation_synthesis.math_config as mc
+
+
+def su_unitary_to_zxz_angles(
+    u: np.ndarray, config: mc.MathConfig
+) -> tuple[rst.Real, rst.Real, rst.Real]:
+    det = u[0, 0] * u[1, 1] - u[0, 1] * u[1, 0]
+    # print(det)
+    # assert config.isclose(det, 1)
+
+    cos_phi = (u[1, 1] * u[0, 0] + u[1, 0] * u[0, 1]).real
+    # clip to the range [-1, 1] to ensure numerical stability.
+    cos_phi = min(1, max(-1, cos_phi))
+    phi = config.arccos(cos_phi)
+
+    sum_half_theta = config.arctan2(u[1, 1].imag, u[1, 1].real)
+    diff_half_theta = config.arctan2(u[1, 0].imag, u[1, 0].real) - 1.5 * config.pi
+
+    theta1 = sum_half_theta + diff_half_theta
+    theta2 = sum_half_theta - diff_half_theta
+
+    return theta1, phi, theta2
+
+
+def rx(phi: rst.Real, config: mc.MathConfig) -> np.ndarray:
+    c = config.cos(phi / 2)
+    s = config.sin(phi / 2)
+    return np.array([[c, -1j * s], [-1j * s, c]])
+
+
+def rz(theta: rst.Real, config: mc.MathConfig) -> np.ndarray:
+    v = config.cos(theta / 2) + 1j * config.sin(theta / 2)
+    return np.diag([v.conjugate(), v])
+
+
+def unitary_from_zxz(
+    theta1: rst.Real, phi: rst.Real, theta2: rst.Real, config: mc.MathConfig
+) -> np.ndarray:
+    return rz(theta1, config) @ rx(phi, config) @ rz(theta2, config)