Speed up ODGP solver by using an iterator to generate solutions dynamically

LICENSE

@@ -18,4 +18,4 @@ FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
 OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
-SOFTWARE.
+SOFTWARE.

pygridsynth/gridsynth.py

@@ -5,6 +5,7 @@
 from .diophantine import NO_SOLUTION, diophantine_dyadic
 from .loop_controller import LoopController
 from .mymath import solve_quadratic, sqrt
+from .quantum_gate import Rz
 from .region import ConvexSet, Ellipse
 from .ring import DRootTwo
 from .synthesis_of_cliffordT import decompose_domega_unitary
@@ -17,14 +18,14 @@ class EpsilonRegion(ConvexSet):
     def __init__(self, theta, epsilon):
         self._theta = theta
         self._epsilon = epsilon
-        self._d = 1 - epsilon**2 / 2
+        self._d = sqrt(1 - epsilon**2 / 4)
         self._z_x = mpmath.cos(-theta / 2)
         self._z_y = mpmath.sin(-theta / 2)
         D_1 = mpmath.matrix([[self._z_x, -self._z_y], [self._z_y, self._z_x]])
-        D_2 = mpmath.matrix([[4 * (1 / epsilon) ** 4, 0], [0, (1 / epsilon) ** 2]])
+        D_2 = mpmath.matrix([[64 * (1 / epsilon) ** 4, 0], [0, 4 * (1 / epsilon) ** 2]])
         D_3 = mpmath.matrix([[self._z_x, self._z_y], [-self._z_y, self._z_x]])
         p = mpmath.matrix([self._d * self._z_x, self._d * self._z_y])
-        ellipse = Ellipse(D_1 * D_2 * D_3, p)
+        ellipse = Ellipse(D_1 @ D_2 @ D_3, p)
         super().__init__(ellipse)
 
     @property
@@ -82,31 +83,22 @@ def generate_complex_unitary(sol):
     )
 
 
-def generate_target_Rz(theta):
-    return mpmath.matrix(
-        [[mpmath.exp(-1.0j * theta / 2), 0], [0, mpmath.exp(1.0j * theta / 2)]]
-    )
-
-
 def error(theta, gates):
-    Rz = generate_target_Rz(theta)
-    U = DOmegaUnitary.from_gates(gates).to_complex_matrix
-    E = U - Rz
-    return sqrt(mpmath.fabs(E[0, 0] * E[1, 1] - E[0, 1] * E[1, 0]))
+    u_target = Rz(theta)
+    u_approx = DOmegaUnitary.from_gates(gates).to_complex_matrix
+    return 2 * sqrt(
+        mpmath.re(1 - mpmath.re(mpmath.conj(u_target[0, 0]) * u_approx[0, 0]) ** 2)
+    )
 
 
 def check(theta, gates):
     t_count = gates.count("T")
     h_count = gates.count("H")
-    U_decomp = DOmegaUnitary.from_gates(gates)
-    # Rz = generate_target_Rz(theta)
-    # U = U_decomp.to_complex_matrix
+    u_approx = DOmegaUnitary.from_gates(gates)
     e = error(theta, gates)
     print(f"{gates=}")
     print(f"{t_count=}, {h_count=}")
-    # print(f"{Rz=}")
-    print(f"U_decomp={U_decomp.to_matrix}")
-    # print(f"{U=}")
+    print(f"u_approx={u_approx.to_matrix}")
     print(f"{e=}")
 
 
@@ -150,7 +142,6 @@ def gridsynth(
         if measure_time:
             time_of_solve_TDGP += time.time() - start
             start = time.time()
-
         for z in sol:
             if (z * z.conj).residue == 0:
                 continue
@@ -183,9 +174,11 @@ def gridsynth(
         k += 1
 
 
-def gridsynth_gates(
+def gridsynth_circuit(
     theta,
     epsilon,
+    wires=[0],
+    decompose_phase_gate=True,
     loop_controller=None,
     verbose=False,
     measure_time=False,
@@ -202,8 +195,33 @@ def gridsynth_gates(
     )
 
     start = time.time() if measure_time else 0.0
-    gates = decompose_domega_unitary(u_approx)
+    circuit = decompose_domega_unitary(
+        u_approx, wires=wires, decompose_phase_gate=decompose_phase_gate
+    )
     if measure_time:
         print(f"time of decompose_domega_unitary: {(time.time() - start) * 1000} ms")
         print(f"total time: {(time.time() - start_total) * 1000} ms")
-    return gates
+
+    return circuit
+
+
+def gridsynth_gates(
+    theta,
+    epsilon,
+    decompose_phase_gate=True,
+    loop_controller=None,
+    verbose=False,
+    measure_time=False,
+    show_graph=False,
+):
+    circuit = gridsynth_circuit(
+        theta=theta,
+        epsilon=epsilon,
+        wires=[0],
+        decompose_phase_gate=decompose_phase_gate,
+        loop_controller=loop_controller,
+        verbose=verbose,
+        measure_time=measure_time,
+        show_graph=show_graph,
+    )
+    return circuit.to_simple_str()

pygridsynth/normal_form.py

@@ -1,5 +1,7 @@
 from enum import Enum
 
+from .quantum_gate import HGate, QuantumCircuit, SGate, SXGate, TGate, WGate
+
 
 class Axis(Enum):
     I = 0
@@ -127,6 +129,19 @@ def from_str(cls, g):
         else:
             raise ValueError
 
+    @classmethod
+    def from_gate(cls, g):
+        if isinstance(g, HGate):
+            return CLIFFORD_H
+        elif isinstance(g, SGate):
+            return CLIFFORD_S
+        elif isinstance(g, SXGate):
+            return CLIFFORD_X
+        elif isinstance(g, WGate):
+            return CLIFFORD_W
+        else:
+            raise ValueError
+
     def __eq__(self, other):
         if isinstance(other, self.__class__):
             return (
@@ -182,15 +197,28 @@ def decompose_tconj(self):
         axis, c1, d1 = self.__class__._tconj(self._a, self._b)
         return axis, self.__class__(0, self._b, c1 + self._c, d1 + self._d)
 
-    def to_gates(self):
+    def to_circuit(self, wires):
         axis, c = self.decompose_coset()
-        return ("" if axis == Axis.I else axis.name) + "X" * c.b + "S" * c.c + "W" * c.d
+        circuit = QuantumCircuit()
+        if axis == Axis.H:
+            circuit.append(HGate(target_qubit=wires[0]))
+        elif axis == Axis.SH:
+            circuit.append(SGate(target_qubit=wires[0]))
+            circuit.append(HGate(target_qubit=wires[0]))
+        for _ in range(c.b):
+            circuit.append(SXGate(target_qubit=wires[0]))
+        for _ in range(c.c):
+            circuit.append(SGate(target_qubit=wires[0]))
+        for _ in range(c.d):
+            circuit.append(WGate())
+        return circuit
 
 
 class NormalForm:
-    def __init__(self, syllables, c):
+    def __init__(self, syllables, c, phase=0):
         self._syllables = syllables
         self._c = c
+        self._phase = phase
 
     @property
     def syllables(self):
@@ -204,13 +232,21 @@ def c(self):
     def c(self, c):
         self._c = c
 
+    @property
+    def phase(self):
+        return self._phase
+
+    @phase.setter
+    def phase(self, phase):
+        self._phase = phase
+
     def __repr__(self):
-        return f"NormalForm({repr(self._syllables)}, {repr(self._c)})"
+        return (
+            f"NormalForm({repr(self._syllables)}, {repr(self._c)}, {repr(self._phase)})"
+        )
 
     def _append_gate(self, g):
-        if g in ["H", "S", "X", "W"]:
-            self.c *= Clifford.from_str(g)
-        elif g == "T":
+        if isinstance(g, TGate):
             axis, new_c = self.c.decompose_tconj()
             if axis == Axis.I:
                 if len(self._syllables) == 0:
@@ -230,21 +266,30 @@ def _append_gate(self, g):
             elif axis == Axis.SH:
                 self._syllables.append(Syllable.SHT)
                 self.c = new_c
+        else:
+            self.c *= Clifford.from_gate(g)
 
     @classmethod
-    def from_gates(cls, gates):
-        normal_form = NormalForm([], CLIFFORD_I)
-        for g in gates:
+    def from_circuit(cls, circuit):
+        normal_form = NormalForm([], CLIFFORD_I, phase=circuit.phase)
+        for g in circuit:
             normal_form._append_gate(g)
         return normal_form
 
-    def to_gates(self):
-        gates = ""
+    def to_circuit(self, wires):
+        circuit = QuantumCircuit(phase=self.phase)
         for syllable in self._syllables:
-            if syllable != Syllable.I:
-                gates += syllable.name
-        gates += self._c.to_gates()
-        return "I" if gates == "" else gates
+            if syllable == Syllable.T:
+                circuit.append(TGate(target_qubit=wires[0]))
+            elif syllable == Syllable.HT:
+                circuit.append(HGate(target_qubit=wires[0]))
+                circuit.append(TGate(target_qubit=wires[0]))
+            elif syllable == Syllable.SHT:
+                circuit.append(SGate(target_qubit=wires[0]))
+                circuit.append(HGate(target_qubit=wires[0]))
+                circuit.append(TGate(target_qubit=wires[0]))
+        circuit += self._c.to_circuit(wires=wires)
+        return circuit
 
 
 CLIFFORD_I = Clifford(0, 0, 0, 0)

pygridsynth/odgp.py

@@ -1,72 +1,71 @@
+import itertools
+
 from .mymath import SQRT2, ceil, floor, floorlog, pow_sqrt2
 from .ring import LAMBDA, DRootTwo, ZRootTwo
 
 
 def _solve_ODGP_internal(I, J):
     if I.width < 0 or J.width < 0:
-        return []
+        return iter([])
     elif I.width > 0 and J.width <= 0:
         sol = _solve_ODGP_internal(J, I)
-        return [beta.conj_sq2 for beta in sol]
+        return map(lambda beta: beta.conj_sq2, sol)
     else:
         (n, _) = (0, 0) if J.width <= 0 else floorlog(J.width, LAMBDA.to_real)
         if n == 0:
-            sol = []
             a_min = ceil((I.l + J.l) / 2)
             a_max = floor((I.r + J.r) / 2)
-            for a in range(a_min, a_max + 1):
+
+            def gen_sol(a):
                 b_min = ceil(SQRT2() * (a - J.r) / 2)
                 b_max = floor(SQRT2() * (a - J.l) / 2)
-                for b in range(b_min, b_max + 1):
-                    sol.append(ZRootTwo(a, b))
-            return sol
+                return map(lambda b: ZRootTwo(a, b), range(b_min, b_max + 1))
+
+            return itertools.chain.from_iterable(map(gen_sol, range(a_min, a_max + 1)))
         else:
             lambda_n = LAMBDA**n
             lambda_inv_n = LAMBDA**-n
             lambda_conj_sq2_n = LAMBDA.conj_sq2**n
             sol = _solve_ODGP_internal(
                 I * lambda_n.to_real, J * lambda_conj_sq2_n.to_real
             )
-            sol = [beta * lambda_inv_n for beta in sol]
-            return sol
+            return map(lambda beta: beta * lambda_inv_n, sol)
 
 
 def solve_ODGP(I, J):
     if I.width < 0 or J.width < 0:
-        return []
+        return iter([])
 
     a = floor((I.l + J.l) / 2)
     b = floor(SQRT2() * (I.l - J.l) / 4)
     alpha = ZRootTwo(a, b)
     sol = _solve_ODGP_internal(I - alpha.to_real, J - alpha.conj_sq2.to_real)
-    sol = [beta + alpha for beta in sol]
-    sol = [
-        beta
-        for beta in sol
-        if I.within(beta.to_real) and J.within(beta.conj_sq2.to_real)
-    ]
+    sol = map(lambda beta: beta + alpha, sol)
+    sol = filter(
+        lambda beta: I.within(beta.to_real) and J.within(beta.conj_sq2.to_real), sol
+    )
     return sol
 
 
 def solve_ODGP_with_parity(I, J, beta):
     p = beta.parity
     sol = solve_ODGP((I - p) * SQRT2() / 2, (J - p) * (-SQRT2()) / 2)
-    sol = [alpha * ZRootTwo(0, 1) + p for alpha in sol]
+    sol = map(lambda alpha: alpha * ZRootTwo(0, 1) + p, sol)
     return sol
 
 
 def solve_scaled_ODGP(I, J, k):
     scale = pow_sqrt2(k)
     sol = solve_ODGP(I * scale, -J * scale if k & 1 else J * scale)
-    return [DRootTwo(alpha, k) for alpha in sol]
+    return map(lambda alpha: DRootTwo(alpha, k), sol)
 
 
 def solve_scaled_ODGP_with_parity(I, J, k, beta):
     if k == 0:
         sol = solve_ODGP_with_parity(I, J, beta.renew_denomexp(0))
-        return [DRootTwo.from_zroottwo(alpha) for alpha in sol]
+        return map(lambda alpha: DRootTwo.from_zroottwo(alpha), sol)
     else:
         p = beta.renew_denomexp(k).parity
         offset = DRootTwo.from_int(0) if p == 0 else DRootTwo.power_of_inv_sqrt2(k)
         sol = solve_scaled_ODGP(I - offset.to_real, J - offset.conj_sq2.to_real, k - 1)
-        return [alpha + offset for alpha in sol]
+        return map(lambda alpha: alpha + offset, sol)