Bound the FiveTuple rescaling by the encoded depth, in SMT and MILP (#65)

src/gates.py

@@ -15,6 +15,12 @@
 
 self_adjoints = ['h', 'x', 'y', 'z', 'cx', 'cz', 'cy', 'ch', 'swap', 'ccx', 'ccz']
 
+# FiveTuple: how much a gate raises k (the sqrt(2)^k denominator), the other gates keep it
+k_increments = {'h': 1, 'ch': 1, 'sx': 2, 'sxdg': 2, 'csx': 2, 'sqrtswap': 2}
+
+def max_k_increment(gate_set) -> int:
+    return max((k_increments.get(gate, 0) for gate in gate_set), default=0)
+
 # OpenQASM 3 definitions of supported gates that are not part of stdgates.inc
 nonstandard_gate_definitions = {
     'sxdg': "gate sxdg a { inv @ sx a; }",

src/generator.py

@@ -751,8 +751,9 @@ def multiply_by_identity_unscaled(result_vec, source_vec):
             return self.And(*eqs)
 
         if self.logic == "QF_LIA" or self.logic == "QF_NIA":
-            if (d // 2) == 0:
-                d = 2
+            # the case split below covers k in [0, d); a larger k would leave r1 and r2 unconstrained, and any
+            # circuit would pass the equivalence check (#65)
+            self.add_assertion(self.LT(k, self.Int(d)))
             for i in range(d):
                 pow2 = self.Int(2 ** i)
                 r1_multiply_by_m = self.And(multiply_by_m_scaled(r1, v1, pow2), multiply_by_identity_unscaled(r2, v2))
@@ -767,14 +768,15 @@ def multiply_by_identity_unscaled(result_vec, source_vec):
         else:
             raise ValueError("rescaling with wrong logic")
 
-    def add_milp_rescaling(self, r1, r2, v1, v2, pair_idx, d):
-        # first encode k as the result of the operation floor(abs(k1 - k2)/2)
+    def add_milp_rescaling(self, r1, r2, v1, v2, pair_idx, max_gap, value_bigM):
+        # first encode k as the result of the operation floor(abs(k1 - k2)/2), where abs(k1 - k2) <= max_gap
+        # value_bigM bounds |v - r| for every entry of v1, v2 and r1, r2
         complex_representation = v1.element_representation
         if self.mode == "milp":
-            k = self.declare_integer(f"k{pair_idx}", lb=0, ub=d)
+            d = max_gap // 2
+            k = self.declare_integer(f"k{pair_idx}", lb=0)
             q = self.declare_bool(f"q{pair_idx}")
-
-            bigM = 2*d + 1
+            bigM = 2*max_gap + 1
             sleq = self.declare_bool(f"sleq{pair_idx}")
             self.add_assertion((v1.k - v2.k) - (2*k + q) <= bigM * sleq)
             self.add_assertion((2*k + q) - (v1.k - v2.k) <= bigM * sleq)
@@ -788,7 +790,7 @@ def add_milp_rescaling(self, r1, r2, v1, v2, pair_idx, d):
                 self.add_assertion(k == self.Sum([constants[i] * i for i in range(d+1)]))
             # now sleq, si, even encode all case splits needed for the rescaling
             # encode all combinations to assign to r1, r2
-            bigM = (2**(d+1)) + 1
+            bigM = value_bigM
             # for every si, I,M is multiplied by 2^i
             for i in range(len(v1)):
                 for j, s in enumerate(constants):

src/synth.py

@@ -19,14 +19,20 @@
 from pysmt.shortcuts import Portfolio, Symbol, Real, And, Equals, Plus, GT, LT, get_env, Int, Or, Not, Implies, GE, LE, Ite, Minus, Div, Times, write_smtlib, Solver
 from pysmt.typing import REAL, INT
 from pareto import Pareto
-from gates import GateSet, supported_gates, check_supported, Circuit
+from gates import GateSet, check_supported, Circuit, k_increments, max_k_increment
 from pulp import *
 from gurobipy import GRB
 from solvers import PortfolioSMTSolver, SMTSolver
 import time
 from logger import Logger
 import numpy as np
 
+def integral_k(k):
+    # a FiveTuple vector is divided by sqrt(2)^k, k has to be an integer (4.0 is accepted)
+    if k is None or int(k) != k:
+        raise ValueError(f"k of a FiveTuple vector must be an integer, got {k!r}")
+    return int(k)
+
 def rus_recovery_pairs(output_vectors, vector_pairs, measured_qubits):
     # (input, output) pairs for the RUS recovery operation: it maps the state after a failed
     # attempt (measured qubits projected to 1) back to the original input state
@@ -54,6 +60,8 @@ def __init__(self) -> None:
         self.circuit_q = None  # qubits the synthesized circuit acts on, the first circuit_q of the q
         self.d = None
         self.max_k = 0
+        self.pair_k = []
+        self.k_increment = 0
         self.encoding_method = None
         self.curr_depth = 0
         self.v = 0
@@ -185,7 +193,7 @@ def add_k_incr(sel, out_k, inp_k, increment):
                         propagate_identities[pos].append(bool_var)
                         propagate_identities[other].append(bool_var)
                     if self.complex_representation == FiveTuple:
-                        add_k_incr(bool_var, out.k, inp.k, 1)
+                        add_k_incr(bool_var, out.k, inp.k, k_increments[gate])
                 elif gate == 's':
                     for pos in range(size):
                         one_flag = (pos >> q) & 1
@@ -251,7 +259,7 @@ def add_k_incr(sel, out_k, inp_k, increment):
                         propagate_identities[pos].append(bool_var)
                         propagate_identities[other].append(bool_var)
                     if self.complex_representation == FiveTuple:
-                        add_k_incr(bool_var, out.k, inp.k, 2)
+                        add_k_incr(bool_var, out.k, inp.k, k_increments[gate])
                 elif gate == 'sxdg':
                     modified_positions = []
                     for pos in range(size):
@@ -264,7 +272,7 @@ def add_k_incr(sel, out_k, inp_k, increment):
                         propagate_identities[pos].append(bool_var)
                         propagate_identities[other].append(bool_var)
                     if self.complex_representation == FiveTuple:
-                        add_k_incr(bool_var, out.k, inp.k, 2)
+                        add_k_incr(bool_var, out.k, inp.k, k_increments[gate])
                 elif gate == 'y':
                     modified_positions = []
                     for pos in range(size):
@@ -363,7 +371,7 @@ def add_k_incr(sel, out_k, inp_k, increment):
                         elif cannot_propagate:
                             saved_operations[pos].append((bool_var, inp[pos]))
                     if self.complex_representation == FiveTuple:
-                        add_k_incr(bool_var, out.k, inp.k, 1)
+                        add_k_incr(bool_var, out.k, inp.k, k_increments[gate])
                 elif gate == 'csx':
                     modified_positions = []
                     for pos in range(size):
@@ -384,7 +392,7 @@ def add_k_incr(sel, out_k, inp_k, increment):
                         elif cannot_propagate:
                             saved_operations[pos].append((bool_var, inp[pos]))
                     if self.complex_representation == FiveTuple:
-                        add_k_incr(bool_var, out.k, inp.k, 2)
+                        add_k_incr(bool_var, out.k, inp.k, k_increments[gate])
                 elif gate == 'cy':
                     modified_positions = []
                     for pos in range(2**self.q):
@@ -500,7 +508,7 @@ def add_k_incr(sel, out_k, inp_k, increment):
                         elif cannot_propagate:
                             saved_operations[pos].append((bool_var, inp[pos]))
                     if self.complex_representation == FiveTuple:
-                        add_k_incr(bool_var, out.k, inp.k, 2)
+                        add_k_incr(bool_var, out.k, inp.k, k_increments[gate])
             elif len(rest) == 4: # three qubit
                 q1 = rest[1]
                 q2 = rest[2]
@@ -696,7 +704,7 @@ def encode_equivalence(self, vectors1, vectors2, depth, already_measured : bool
                         for i in range(8):
                             selectors.append(self.gen.declare_bool(f"global_phase_sel_{i}"))
                         self.gen.ExactlyOne(*selectors)
-                        vec = Vector(q=2**self.q, generator=self.gen, element_representation=self.complex_representation, k=0, n=self.q, name=f"Multiplied_by_global_phase_{vec2.name}", bound=2**self.max_k)
+                        vec = Vector(q=2**self.q, generator=self.gen, element_representation=self.complex_representation, k=0, n=self.q, name=f"Multiplied_by_global_phase_{vec2.name}", bound=vec2.bound)
                         vec_copy = vec2.copy()
                         for i in range(8):
                             # add selection shift
@@ -714,12 +722,21 @@ def encode_equivalence(self, vectors1, vectors2, depth, already_measured : bool
                         vec_target = vec
 
                     # fivetuple rescaling
-                    rescaled1 = Vector(q=2**self.q, generator=self.gen, element_representation=self.complex_representation, k=0, n=self.q, name=f"Rescaled1_{vec1.name}", bound=2**self.max_k)
-                    rescaled2 = Vector(q=2**self.q, generator=self.gen, element_representation=self.complex_representation, k=0, n=self.q, name=f"Rescaled2_{vec2.name}", bound=2**self.max_k)
+                    # after `depth` layers, k lies in [k_in, k_in + self.k_increment*depth], which bounds |k - k_target|
+                    # (self.max_k is fixed from the input and too small past its depth)
+                    k_in, k_target = self.pair_k[pair_idx]
+                    max_gap = max(k_in + self.k_increment * depth - k_target, k_target - k_in)
+                    # the rescaling multiplies by at most 2^(max_gap // 2), and an odd gap adds two components
+                    value_bound = max(vec1.bound, vec_target.bound)
+                    rescaled_bound = 2 * value_bound * 2**(max_gap // 2)
+                    rescaled1 = Vector(q=2**self.q, generator=self.gen, element_representation=self.complex_representation, k=0, n=self.q, name=f"Rescaled1_{vec1.name}", bound=rescaled_bound)
+                    # vec2 is the target, the same at every depth: without the depth in the name, MILP would get back the
+                    # first depth's PuLP variables, with that depth's (smallest) bound
+                    rescaled2 = Vector(q=2**self.q, generator=self.gen, element_representation=self.complex_representation, k=0, n=self.q, name=f"Rescaled2_{vec2.name}_{depth}", bound=rescaled_bound)
                     if self.gen.mode == "milp" or self.gen.mode == "gurobi":
-                        self.gen.add_milp_rescaling(rescaled1, rescaled2, vec1, vec_target, pair_idx, self.max_k)
+                        self.gen.add_milp_rescaling(rescaled1, rescaled2, vec1, vec_target, pair_idx, max_gap, value_bound + rescaled_bound)
                     else:
-                        self.gen.add_rescaling(rescaled1, rescaled2, vec1, vec_target, pair_idx, self.max_k)
+                        self.gen.add_rescaling(rescaled1, rescaled2, vec1, vec_target, pair_idx, max_gap // 2 + 1)
                     for i in range(2**self.q):
                         self.gen.add_assertion(self.gen.Equals(rescaled1[i], rescaled2[i]))
 
@@ -901,6 +918,9 @@ def synthesis(self, qasm_file=None, matrix=None, vectors=None, vector_pairs=None
             self.gate_set.add_gate("id", weight=0, qubits=1)
         self.circuit_q = self.q // 2 if vectors == "jamiolkowski" else self.q
         self.v = len(vector_pairs)
+        if self.complex_representation != Complex:
+            for state in (state for pair in vector_pairs for state in pair):
+                state.k = integral_k(state.k)
         self.curr_depth = self.d
         max = 0
         for pair_idx, (input_vector, output_vector) in enumerate(vector_pairs):
@@ -1005,6 +1025,10 @@ def _synth(self, vector_pairs, output_qasm="circuit.qasm") -> bool:
         # encode vector bounds
         inter_states = []
         target_states = []
+        # k of every (input, target) pair and its largest increase per layer, used to bound the FiveTuple k
+        if self.complex_representation != Complex:
+            self.pair_k = [(input_state.k, output_state.k) for input_state, output_state in vector_pairs]
+            self.k_increment = max_k_increment(self.gate_set)
         bounds = [1] * (self.d + 1)
         if self.complex_representation != Complex:
             for d in range(self.d + 1):
@@ -1029,7 +1053,7 @@ def _synth(self, vector_pairs, output_qasm="circuit.qasm") -> bool:
             for d in range(self.d + 1):
                 Inter = None
                 if self.complex_representation == FiveTuple:
-                    Inter = Vector(q=2**self.q, name=f"I_{pair_idx}_{d}", generator=self.gen, element_representation=self.complex_representation, k=input_state.k if d == 0 else 0, n = input_state.n, bound = bounds[d], k_bound = 2 * (d+1))
+                    Inter = Vector(q=2**self.q, name=f"I_{pair_idx}_{d}", generator=self.gen, element_representation=self.complex_representation, k=input_state.k if d == 0 else 0, n = input_state.n, bound = bounds[d], k_bound = input_state.k + self.k_increment * d)
                 else:
                     Inter = Vector(q=2**self.q, name=f"I_{pair_idx}_{d}", generator=self.gen, element_representation=self.complex_representation, bound = bounds[d], k_bound = 2 * (d+1))
 
@@ -1132,7 +1156,7 @@ def _synth(self, vector_pairs, output_qasm="circuit.qasm") -> bool:
                     for pair_idx in range(len(vector_pairs)):
                         if self.curr_depth > self.d:
                             if self.complex_representation == FiveTuple:
-                                state = Vector(q=2**self.q, name=f"I_{pair_idx}_{self.curr_depth}", generator=self.gen, element_representation=self.complex_representation, k=0, n = vector_pairs[pair_idx][0].n, bound = self.bound, k_bound = 2 * (self.curr_depth+1))
+                                state = Vector(q=2**self.q, name=f"I_{pair_idx}_{self.curr_depth}", generator=self.gen, element_representation=self.complex_representation, k=0, n = vector_pairs[pair_idx][0].n, bound = self.bound, k_bound = vector_pairs[pair_idx][0].k + self.k_increment * self.curr_depth)
                             else:
                                 state = Vector(q=2**self.q, name=f"I_{pair_idx}_{self.curr_depth}", generator=self.gen, element_representation=self.complex_representation, bound = self.bound, k_bound = 2 * (self.curr_depth+1))
                             inter_states[pair_idx].append(state)