diff --git a/src/gates.py b/src/gates.py index 2e82e88..b2f203d 100644 --- a/src/gates.py +++ b/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; }", diff --git a/src/generator.py b/src/generator.py index d9c9acd..86a5e49 100644 --- a/src/generator.py +++ b/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): diff --git a/src/synth.py b/src/synth.py index 081ca93..c820ca4 100644 --- a/src/synth.py +++ b/src/synth.py @@ -19,7 +19,7 @@ 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 @@ -27,6 +27,12 @@ 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) diff --git a/tests/test_rescaling_bound.py b/tests/test_rescaling_bound.py new file mode 100644 index 0000000..8e8d767 --- /dev/null +++ b/tests/test_rescaling_bound.py @@ -0,0 +1,139 @@ +"""FiveTuple equivalence when the k gap between the output and the target exceeds the input's max_k (#65).""" + +import numpy as np +import pytest +from qiskit import QuantumCircuit +from qiskit.quantum_info import Operator + +from complex.fivetuple import FiveTuple +from complex.vector import Vector +from gates import GateSet, max_k_increment +from synth import Synthesizer + +SMT = dict(solving="smt", solver="yices2", complex_representation="FiveTuple") + + +def write_qasm(path, gates): + with open(path, "w") as f: + f.write('OPENQASM 3.0;\ninclude "stdgates.inc";\nqubit[1] q;\n') + f.writelines(f"{gate} q[0];\n" for gate in gates) + return str(path) + + +def operator(gates): + # reference unitary computed by Qiskit, independently of MILQ's simulator + qc = QuantumCircuit(1) + for gate in gates: + getattr(qc, gate)(0) + return Operator(qc).data + + +def synthesized_gates(circuit): + return [gate.name for gate in circuit.gates if gate is not None and gate.name != "id"] + + +def test_shorter_depth_bound_still_needs_an_exact_circuit(workdir): + # t x t x = w*I; with d=1 the search goes past the input depth, where the old bound left the equivalence + # unconstrained and h h h h (= I) was accepted + qasm_file = write_qasm(workdir / "in.qasm", ["t", "x", "t", "x"]) + res, circuit, _ = Synthesizer().synthesis(qasm_file=qasm_file, vectors="all", mode="incremental", d=1, + output_qasm="out.qasm", **SMT) + assert res + np.testing.assert_allclose(operator(synthesized_gates(circuit)), operator(["t", "x", "t", "x"]), atol=1e-9) + + +def test_unreachable_target_is_not_synthesized(workdir): + # X is not a product of H gates, at any depth + qasm_file = write_qasm(workdir / "in.qasm", ["x"]) + res, circuit, _ = Synthesizer().synthesis(qasm_file=qasm_file, vectors="all", mode="incremental", + gate_set=GateSet(["h"]), max_depth=6, output_qasm="out.qasm", **SMT) + assert (res, circuit) == (False, None) + + +def test_restricted_gate_set_gives_an_exact_circuit(workdir): + # Z = h sx sx h; the old bound accepted sx sx (= X) + qasm_file = write_qasm(workdir / "in.qasm", ["z"]) + res, circuit, _ = Synthesizer().synthesis(qasm_file=qasm_file, vectors="all", mode="incremental", + gate_set=GateSet(["h", "sx"]), output_qasm="out.qasm", **SMT) + assert res + gates = synthesized_gates(circuit) + assert len(gates) == 4 + np.testing.assert_allclose(operator(gates), operator(["z"]), atol=1e-9) + + +def zero_to_one(k_target): + # |0> -> |1>, with the target written as 4/sqrt(2)^k_target + inp = Vector(q=2, element_representation=FiveTuple, k=0) + out = Vector(q=2, element_representation=FiveTuple, k=k_target) + inp[0] = FiveTuple(a=1, b=0, c=0, d=0) + out[1] = FiveTuple(a=4, b=0, c=0, d=0) + return [(inp, out)] + + +@pytest.mark.parametrize("k_target", [4, 4.0]) +@pytest.mark.parametrize("mode", ["basic", "incremental"]) +def test_custom_target_with_a_large_k(workdir, mode, k_target): + # the old bound only covered the gap of the input depth, and the empty circuit was accepted + res, circuit, _ = Synthesizer().synthesis(vectors="custom", vector_pairs=zero_to_one(k_target), q=1, d=1, + gate_set=GateSet(["x", "h"]), mode=mode, output_qasm="out.qasm", **SMT) + assert res + assert synthesized_gates(circuit) == ["x"] + + +def test_float_k_reaches_pysmt_as_an_integer(workdir): + # yices2 tolerates a k of 4.0 in the formula; pySMT crashed on it, so the vectors must carry the int + res, circuit, _ = Synthesizer().synthesis(vectors="custom", vector_pairs=zero_to_one(4.0), q=1, d=1, + gate_set=GateSet(["x", "h"]), mode="basic", output_qasm="out.qasm", + solving="pysmt", solver="z3", complex_representation="FiveTuple") + assert res + assert synthesized_gates(circuit) == ["x"] + + +def test_custom_target_with_a_fractional_k(workdir): + with pytest.raises(ValueError, match="k of a FiveTuple vector must be an integer, got 4.5"): + Synthesizer().synthesis(vectors="custom", vector_pairs=zero_to_one(4.5), q=1, d=1, gate_set=GateSet(["x", "h"]), + mode="incremental", output_qasm="out.qasm", **SMT) + + +MILP = dict(solving="milp", solver="cbc", complex_representation="FiveTuple") + + +@pytest.mark.parametrize("target, gate_set, expected", [ + # past the input's max_k, the stale k bound and big-M made the correct circuit infeasible (#65 F15) + (["x"], ["sx", "z"], ["sx", "sx"]), + (["z"], ["h", "sx"], ["h", "sx", "sx", "h"]), +]) +def test_milp_finds_circuits_past_the_input_depth(workdir, target, gate_set, expected): + qasm_file = write_qasm(workdir / "in.qasm", target) + res, circuit, _ = Synthesizer().synthesis(qasm_file=qasm_file, vectors="all", mode="incremental", + gate_set=GateSet(gate_set), max_depth=6, output_qasm="out.qasm", **MILP) + assert res + gates = synthesized_gates(circuit) + assert len(gates) == len(expected) + np.testing.assert_allclose(operator(gates), operator(target), atol=1e-9) + + +def test_milp_custom_target_with_a_large_k(workdir): + res, circuit, _ = Synthesizer().synthesis(vectors="custom", vector_pairs=zero_to_one(4), q=1, d=1, + gate_set=GateSet(["x", "h"]), mode="basic", output_qasm="out.qasm", **MILP) + assert res + assert synthesized_gates(circuit) == ["x"] + + +@pytest.mark.parametrize("mode", ["basic", "incremental"]) +def test_milp_input_with_a_large_k(workdir, mode): + # the intermediate states' k was capped at 2*(depth+1), below the input's own k of 6 + inp = Vector(q=2, element_representation=FiveTuple, k=6) + out = Vector(q=2, element_representation=FiveTuple, k=6) + inp[0] = FiveTuple(a=1, b=0, c=0, d=0) + out[1] = FiveTuple(a=1, b=0, c=0, d=0) + res, circuit, _ = Synthesizer().synthesis(vectors="custom", vector_pairs=[(inp, out)], q=1, d=1, + gate_set=GateSet(["x", "h"]), mode=mode, output_qasm="out.qasm", **MILP) + assert res + assert synthesized_gates(circuit) == ["x"] + + +@pytest.mark.parametrize("gates, expected", [(["x", "z", "t"], 0), (["h", "t"], 1), (["h", "sx"], 2), (["cx", "csx"], 2)]) +def test_k_increment_follows_the_gate_set(gates, expected): + # the k bounds of the encoding come from this, a gate missing from k_increments would make them too small + assert max_k_increment(GateSet(gates)) == expected