From 9eec3f506f3463ae83f7ebfc301f18baa866451d Mon Sep 17 00:00:00 2001 From: Ondra Lengal Date: Tue, 13 Oct 2026 10:37:25 +0200 Subject: [PATCH 1/4] fix: bound the FiveTuple rescaling by the depth, not by the input's max_k add_rescaling covered k = floor(|k_out - k_target| / 2) only below Synthesizer.max_k, which is computed once from the input. Past the input depth (incremental search, -d, restricted gate sets, custom vector pairs) the rescaled vectors were left unconstrained, so any circuit passed the equivalence check: t;x;t;x with -d 1 gave h;h;h;h. encode_equivalence now covers every gap the encoded depth can produce (each layer raises k by at most 2), and add_rescaling asserts k < d, so a gap outside the case split is UNSAT instead of unconstrained. Part of #65 (the MILP counterpart, F15, is not addressed here). Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_0177AhE9SXm3iJ1KP8GTSk6o --- src/generator.py | 5 ++- src/synth.py | 12 +++++- tests/test_rescaling_bound.py | 74 +++++++++++++++++++++++++++++++++++ 3 files changed, 88 insertions(+), 3 deletions(-) create mode 100644 tests/test_rescaling_bound.py diff --git a/src/generator.py b/src/generator.py index 25d02a7..806571d 100644 --- a/src/generator.py +++ b/src/generator.py @@ -752,8 +752,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)) diff --git a/src/synth.py b/src/synth.py index d985e0c..238b148 100644 --- a/src/synth.py +++ b/src/synth.py @@ -26,6 +26,9 @@ from logger import Logger import numpy as np +# the largest increase of k (FiveTuple) caused by one layer: sx, sxdg, csx and sqrtswap add 2, h and ch add 1 +MAX_K_INCREMENT = 2 + 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 @@ -52,6 +55,7 @@ def __init__(self) -> None: self.q = None self.d = None self.max_k = 0 + self.pair_k = [] self.encoding_method = None self.curr_depth = 0 self.v = 0 @@ -788,7 +792,11 @@ def encode_equivalence(self, vectors1, vectors2, depth, already_measured : bool 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) else: - self.gen.add_rescaling(rescaled1, rescaled2, vec1, vec_target, pair_idx, self.max_k) + # after `depth` layers, k lies in [k_in, k_in + MAX_K_INCREMENT*depth]; cover every possible + # floor(|k - k_target| / 2) (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 + MAX_K_INCREMENT * depth - k_target, k_target - k_in) + 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])) @@ -1043,6 +1051,8 @@ def _synth(self, vector_pairs, output_qasm="circuit.qasm") -> bool: # encode vector bounds inter_states = [] target_states = [] + # k of every (input, target) pair, used to bound the FiveTuple rescaling + self.pair_k = [(input_state.k, output_state.k) for input_state, output_state in vector_pairs] bounds = [1] * (self.d + 1) if self.complex_representation != Complex: for d in range(self.d + 1): diff --git a/tests/test_rescaling_bound.py b/tests/test_rescaling_bound.py new file mode 100644 index 0000000..34a09df --- /dev/null +++ b/tests/test_rescaling_bound.py @@ -0,0 +1,74 @@ +"""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 +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) + + +@pytest.mark.parametrize("mode", ["basic", "incremental"]) +def test_custom_target_with_a_large_k(workdir, mode): + # |0> -> |1>, with the target written as 4/sqrt(2)^4: the old bound only covered the gap of the input depth + inp = Vector(q=2, element_representation=FiveTuple, k=0) + out = Vector(q=2, element_representation=FiveTuple, k=4) + inp[0] = FiveTuple(a=1, b=0, c=0, d=0) + out[1] = FiveTuple(a=4, 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", **SMT) + assert res + assert synthesized_gates(circuit) == ["x"] From 18e94ba631010673a9a72fbf8c2d15c72235dee0 Mon Sep 17 00:00:00 2001 From: Ondra Lengal Date: Tue, 13 Oct 2026 11:59:05 +0200 Subject: [PATCH 2/4] fix: reject a non-integral FiveTuple k in vector pairs The rescaling bound is computed from the pairs' k, so a custom pair with k=4.0 raised a TypeError in range(). Integral floats are now accepted and anything else is rejected with a ValueError. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_0177AhE9SXm3iJ1KP8GTSk6o --- src/synth.py | 9 ++++++++- tests/test_rescaling_bound.py | 22 +++++++++++++++++----- 2 files changed, 25 insertions(+), 6 deletions(-) diff --git a/src/synth.py b/src/synth.py index 238b148..ea6c517 100644 --- a/src/synth.py +++ b/src/synth.py @@ -29,6 +29,12 @@ # the largest increase of k (FiveTuple) caused by one layer: sx, sxdg, csx and sqrtswap add 2, h and ch add 1 MAX_K_INCREMENT = 2 +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 @@ -1052,7 +1058,8 @@ def _synth(self, vector_pairs, output_qasm="circuit.qasm") -> bool: inter_states = [] target_states = [] # k of every (input, target) pair, used to bound the FiveTuple rescaling - self.pair_k = [(input_state.k, output_state.k) for input_state, output_state in vector_pairs] + if self.complex_representation != Complex: + self.pair_k = [(integral_k(input_state.k), integral_k(output_state.k)) for input_state, output_state in vector_pairs] bounds = [1] * (self.d + 1) if self.complex_representation != Complex: for d in range(self.d + 1): diff --git a/tests/test_rescaling_bound.py b/tests/test_rescaling_bound.py index 34a09df..ee490ba 100644 --- a/tests/test_rescaling_bound.py +++ b/tests/test_rescaling_bound.py @@ -61,14 +61,26 @@ def test_restricted_gate_set_gives_an_exact_circuit(workdir): np.testing.assert_allclose(operator(gates), operator(["z"]), atol=1e-9) -@pytest.mark.parametrize("mode", ["basic", "incremental"]) -def test_custom_target_with_a_large_k(workdir, mode): - # |0> -> |1>, with the target written as 4/sqrt(2)^4: the old bound only covered the gap of the input depth +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=4) + 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) - res, circuit, _ = Synthesizer().synthesis(vectors="custom", vector_pairs=[(inp, out)], q=1, d=1, + 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_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) From f269ec7bebf9b77d84fef563f4e611ba498ec857 Mon Sep 17 00:00:00 2001 From: Jakub Havlik Date: Sat, 17 Oct 2026 19:32:19 +0200 Subject: [PATCH 3/4] fix: integral_k check before synthesis even begins -- crashed pysmt before --- src/synth.py | 5 ++++- tests/test_rescaling_bound.py | 9 +++++++++ 2 files changed, 13 insertions(+), 1 deletion(-) diff --git a/src/synth.py b/src/synth.py index ea6c517..c0cf523 100644 --- a/src/synth.py +++ b/src/synth.py @@ -957,6 +957,9 @@ def synthesis(self, qasm_file=None, matrix=None, vectors=None, vector_pairs=None self.gate_set = gate_set self.qubits_to_measure = [] 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): @@ -1059,7 +1062,7 @@ def _synth(self, vector_pairs, output_qasm="circuit.qasm") -> bool: target_states = [] # k of every (input, target) pair, used to bound the FiveTuple rescaling if self.complex_representation != Complex: - self.pair_k = [(integral_k(input_state.k), integral_k(output_state.k)) for input_state, output_state in vector_pairs] + self.pair_k = [(input_state.k, output_state.k) for input_state, output_state in vector_pairs] bounds = [1] * (self.d + 1) if self.complex_representation != Complex: for d in range(self.d + 1): diff --git a/tests/test_rescaling_bound.py b/tests/test_rescaling_bound.py index ee490ba..046d198 100644 --- a/tests/test_rescaling_bound.py +++ b/tests/test_rescaling_bound.py @@ -80,6 +80,15 @@ def test_custom_target_with_a_large_k(workdir, mode, k_target): 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"]), From 0133d27c1a47430eda0f1ec982d42e7a201af341 Mon Sep 17 00:00:00 2001 From: Jakub Havlik Date: Sun, 18 Oct 2026 00:44:42 +0200 Subject: [PATCH 4/4] fix: bound the MILP FiveTuple rescaling by the encoded depth (#65) add_milp_rescaling took its k bound and both big-Ms from Synthesizer.max_k, which is fixed from the input. Past the input depth, or with custom pairs whose k is not reduced, the correct circuit was infeasible: CBC found no sx;sx for X from {sx, z}, nor h;sx;sx;h for Z from {h, sx}. - The MILP rescaling uses the same per-depth gap bound as the SMT one, with a parity big-M of 2*max_gap + 1 and a value big-M from the vectors' bounds. - k of the rescaling has no ub: incremental search gets back the first depth's PuLP variable with its stale bound, the one-hot bounds k anyway. - The intermediate states' k bound includes the input's k. - The per-layer increase of k comes from a k_increments table in gates.py, used by encode_layer and by the bounds, instead of a constant 2. --- src/gates.py | 6 +++++ src/generator.py | 13 +++++----- src/synth.py | 44 +++++++++++++++++---------------- tests/test_rescaling_bound.py | 46 ++++++++++++++++++++++++++++++++++- 4 files changed, 81 insertions(+), 28 deletions(-) diff --git a/src/gates.py b/src/gates.py index 5d3e45b..67cda08 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 806571d..e36331a 100644 --- a/src/generator.py +++ b/src/generator.py @@ -769,14 +769,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) @@ -790,7 +791,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 c0cf523..f0592f3 100644 --- a/src/synth.py +++ b/src/synth.py @@ -18,7 +18,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 @@ -26,9 +26,6 @@ from logger import Logger import numpy as np -# the largest increase of k (FiveTuple) caused by one layer: sx, sxdg, csx and sqrtswap add 2, h and ch add 1 -MAX_K_INCREMENT = 2 - 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: @@ -62,6 +59,7 @@ def __init__(self) -> None: 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 +183,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 self.basis == "pauli": raise NotImplementedError("Pauli basis not supported for H gate") elif gate == 's': @@ -269,7 +267,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 self.basis == "pauli": raise NotImplementedError("Pauli basis not supported for SX gate") elif gate == 'sxdg': @@ -285,7 +283,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 self.basis == "pauli": raise NotImplementedError("Pauli basis not supported for SXDG gate") elif gate == 'y': @@ -402,7 +400,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 self.basis == "pauli": raise NotImplementedError("Pauli basis not supported for CH gate") elif gate == 'csx': @@ -426,7 +424,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 self.basis == "pauli": raise NotImplementedError("Pauli basis not supported for CSX gate") elif gate == 'cy': @@ -563,7 +561,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 self.basis == "pauli": raise NotImplementedError("Pauli basis not supported for SQRTSWAP gate") elif len(rest) == 4: # three qubit @@ -775,7 +773,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 @@ -793,15 +791,18 @@ 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) + 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=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: - # after `depth` layers, k lies in [k_in, k_in + MAX_K_INCREMENT*depth]; cover every possible - # floor(|k - k_target| / 2) (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 + MAX_K_INCREMENT * depth - k_target, k_target - k_in) 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])) @@ -1060,9 +1061,10 @@ def _synth(self, vector_pairs, output_qasm="circuit.qasm") -> bool: # encode vector bounds inter_states = [] target_states = [] - # k of every (input, target) pair, used to bound the FiveTuple rescaling + # 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): @@ -1088,7 +1090,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 = class_to_use(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 = class_to_use(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 = class_to_use(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)) @@ -1189,7 +1191,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 = class_to_use(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 = class_to_use(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 = class_to_use(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 index 046d198..8e8d767 100644 --- a/tests/test_rescaling_bound.py +++ b/tests/test_rescaling_bound.py @@ -7,7 +7,7 @@ from complex.fivetuple import FiveTuple from complex.vector import Vector -from gates import GateSet +from gates import GateSet, max_k_increment from synth import Synthesizer SMT = dict(solving="smt", solver="yices2", complex_representation="FiveTuple") @@ -93,3 +93,47 @@ 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