From 00821a80ec86dc3ba759c02fc69682ebdec48957 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Pablo=20Le=20H=C3=A9naff?= Date: Thu, 25 Jan 2024 15:15:42 +0100 Subject: [PATCH] Create merge pass - Circuit.merge_single_qubit_gates - mcKay decomposer is now only the decomposition - Circuit.decompose is introduced, pass it the McKay decomposer - Decomposer abstract base class - The decomposer/replacer now checks that every used decomposition is valid in terms of circuit matrix! --- opensquirrel/circuit.py | 40 +++++----- opensquirrel/common.py | 13 +++ opensquirrel/default_gates.py | 37 ++++++++- opensquirrel/mckay_decomposer.py | 112 +++++++------------------- opensquirrel/merger.py | 84 ++++++++++++++++++++ opensquirrel/replacer.py | 112 ++++++++++++++++++++++---- opensquirrel/squirrel_ir.py | 65 ++++++++++++--- test/ir_equality_test_base.py | 38 +++++++++ test/test_integration.py | 21 +++-- test/test_mckay_decomposer.py | 132 +++++++++++++------------------ test/test_merger.py | 108 +++++++++++++++++++++++++ test/test_replacer.py | 125 +++++++++++++++++++++++++++-- 12 files changed, 665 insertions(+), 222 deletions(-) create mode 100644 opensquirrel/merger.py create mode 100644 test/ir_equality_test_base.py create mode 100644 test/test_merger.py diff --git a/opensquirrel/circuit.py b/opensquirrel/circuit.py index 1e522a00..dd4676fa 100644 --- a/opensquirrel/circuit.py +++ b/opensquirrel/circuit.py @@ -1,11 +1,12 @@ -from typing import Callable, Dict +from typing import Callable, Dict, List import numpy as np import opensquirrel.parsing.antlr.squirrel_ir_from_string -from opensquirrel import circuit_matrix_calculator, mckay_decomposer, replacer, writer +from opensquirrel import circuit_matrix_calculator, mckay_decomposer, merger, replacer, writer from opensquirrel.default_gates import default_gate_aliases, default_gate_set from opensquirrel.parsing.libqasm.libqasm_ir_creator import LibqasmIRCreator +from opensquirrel.replacer import Decomposer from opensquirrel.squirrel_ir import Gate, SquirrelIR @@ -21,7 +22,7 @@ class Circuit: h q[0] - >>> c.decompose_mckay() + >>> c.decompose(decomposer=mckay_decomposer.McKayDecomposer) >>> c version 3.0 @@ -84,27 +85,25 @@ def number_of_qubits(self) -> int: def qubit_register_name(self) -> str: return self.squirrel_ir.qubit_register_name - def decompose_mckay(self): - """Perform gate fusion on all one-qubit gates and decompose them in the McKay style. - - * all one-qubit gates on same qubit are merged together, without attempting to commute any gate - * two-or-more-qubit gates are left as-is - * merged one-qubit gates are decomposed according to McKay decomposition, that is: - gate ----> Rz.Rx(pi/2).Rz.Rx(pi/2).Rz - * _global phase is deemed irrelevant_, therefore a simulator backend might produce different output - for the input and output circuit - those outputs should be equivalent modulo global phase. + def merge_single_qubit_gates(self): + """Merge all consecutive 1-qubit gates in the circuit. + Gates obtained from merging other gates become anonymous gates. """ - self.squirrel_ir = mckay_decomposer.decompose_mckay(self.squirrel_ir) # FIXME: inplace + merger.merge_single_qubit_gates(self.squirrel_ir) - def replace(self, gate_name: str, f): - """Manually replace occurrences of a given gate with a list of gates. + def decompose(self, decomposer: Decomposer): + """Generic decomposition pass. It applies the given decomposer function + to every gate in the circuit.""" + replacer.decompose(self.squirrel_ir, decomposer) - * this can be called decomposition - but it's the least fancy version of it - * function parameter gives the decomposition based on parameters of original gate + def replace(self, gate_generator: Callable[..., Gate], f): + """Manually replace occurrences of a given gate with a list of gates. + `f` is a callable that takes the arguments of the gate that is to be replaced + and returns the decomposition as a list of gates. """ - replacer.replace(self.squirrel_ir, gate_name, f) + replacer.replace(self.squirrel_ir, gate_generator, f) def test_get_circuit_matrix(self) -> np.ndarray: """Get the (large) unitary matrix corresponding to the circuit. @@ -117,9 +116,6 @@ def test_get_circuit_matrix(self) -> np.ndarray: return circuit_matrix_calculator.get_circuit_matrix(self.squirrel_ir) def __repr__(self) -> str: - """Write the circuit to a cQasm3 string. - - * comments are removed - """ + """Write the circuit to a cQasm3 string.""" return writer.squirrel_ir_to_string(self.squirrel_ir) diff --git a/opensquirrel/common.py b/opensquirrel/common.py index 0e071054..cb002587 100644 --- a/opensquirrel/common.py +++ b/opensquirrel/common.py @@ -41,3 +41,16 @@ def can1(axis: Tuple[float, float, float], angle: float, phase: float = 0) -> np ) return result + + +def are_matrices_equivalent_up_to_global_phase(matrix_a: np.ndarray, matrix_b: np.ndarray) -> bool: + first_non_zero = next( + (i, j) for i in range(matrix_a.shape[0]) for j in range(matrix_a.shape[1]) if abs(matrix_a[i, j]) > ATOL + ) + + if abs(matrix_b[first_non_zero]) < ATOL: + return False + + phase_difference = matrix_a[first_non_zero] / matrix_b[first_non_zero] + + return np.allclose(matrix_a, phase_difference * matrix_b) diff --git a/opensquirrel/default_gates.py b/opensquirrel/default_gates.py index 4c38b9aa..f34d3652 100644 --- a/opensquirrel/default_gates.py +++ b/opensquirrel/default_gates.py @@ -18,6 +18,11 @@ def x90(q: Qubit) -> Gate: return BlochSphereRotation(qubit=q, axis=(1, 0, 0), angle=math.pi / 2, phase=0) +@named_gate +def xm90(q: Qubit) -> Gate: + return BlochSphereRotation(qubit=q, axis=(1, 0, 0), angle=-math.pi / 2, phase=0) + + @named_gate def y(q: Qubit) -> Gate: return BlochSphereRotation(qubit=q, axis=(0, 1, 0), angle=math.pi, phase=math.pi / 2) @@ -76,5 +81,35 @@ def crk(control: Qubit, target: Qubit, k: Int) -> Gate: return ControlledGate(control, BlochSphereRotation(qubit=target, axis=(0, 0, 1), angle=theta, phase=theta / 2)) -default_gate_set = [h, x, x90, y, y90, z, z90, cz, cr, crk, cnot, rx, ry, rz, x] +@named_gate +def swap(q1: Qubit, q2: Qubit) -> Gate: + return MatrixGate( + np.array( + [ + [1, 0, 0, 0], + [0, 0, 1, 0], + [0, 1, 0, 0], + [0, 0, 0, 1], + ] + ), + [q1, q2], + ) + + +@named_gate +def sqrt_swap(q1: Qubit, q2: Qubit) -> Gate: + return MatrixGate( + np.array( + [ + [1, 0, 0, 0], + [0, (1 + 1j) / 2, (1 - 1j) / 2, 0], + [0, (1 - 1j) / 2, (1 + 1j) / 2, 0], + [0, 0, 0, 1], + ] + ), + [q1, q2], + ) + + +default_gate_set = [h, x, x90, xm90, y, y90, z, z90, cz, cr, crk, cnot, rx, ry, rz, x, swap, sqrt_swap] default_gate_aliases = {"X": x, "RX": rx, "RY": ry, "RZ": rz, "Hadamard": h, "H": h} diff --git a/opensquirrel/mckay_decomposer.py b/opensquirrel/mckay_decomposer.py index 60ebf2f0..01610d63 100644 --- a/opensquirrel/mckay_decomposer.py +++ b/opensquirrel/mckay_decomposer.py @@ -5,27 +5,36 @@ from opensquirrel.common import ATOL, normalize_angle from opensquirrel.default_gates import rz, x90 +from opensquirrel.replacer import Decomposer from opensquirrel.squirrel_ir import BlochSphereRotation, Float, Gate, Qubit, SquirrelIR -class _McKayDecomposerImpl: - def __init__(self, number_of_qubits: int, qubit_register_name: str): - self.output = SquirrelIR(number_of_qubits=number_of_qubits, qubit_register_name=qubit_register_name) - self.accumulated_1q_gates = {} +class McKayDecomposer(Decomposer): + @staticmethod + def decompose(g: Gate) -> [Gate]: + """Return the McKay decomposition of a 1-qubit gate as a list of gates. + gate ----> Rz.Rx(pi/2).Rz.Rx(pi/2).Rz - def decompose_and_add(self, qubit: Qubit, angle: float, axis: Tuple[float, float, float]): - if abs(angle) < ATOL: - return + The global phase is deemed _irrelevant_, therefore a simulator backend might produce different output + for the input and output - the results should be equivalent modulo global phase. + + Relevant literature: https://arxiv.org/abs/1612.00858 + """ + if not isinstance(g, BlochSphereRotation): + return [g] + + if abs(g.angle) < ATOL: + return [] # McKay decomposition - za_mod = sqrt(cos(angle / 2) ** 2 + (axis[2] * sin(angle / 2)) ** 2) - zb_mod = abs(sin(angle / 2)) * sqrt(axis[0] ** 2 + axis[1] ** 2) + za_mod = sqrt(cos(g.angle / 2) ** 2 + (g.axis[2] * sin(g.angle / 2)) ** 2) + zb_mod = abs(sin(g.angle / 2)) * sqrt(g.axis[0] ** 2 + g.axis[1] ** 2) theta = pi - 2 * atan2(zb_mod, za_mod) - alpha = atan2(-sin(angle / 2) * axis[2], cos(angle / 2)) - beta = atan2(-sin(angle / 2) * axis[0], -sin(angle / 2) * axis[1]) + alpha = atan2(-sin(g.angle / 2) * g.axis[2], cos(g.angle / 2)) + beta = atan2(-sin(g.angle / 2) * g.axis[0], -sin(g.angle / 2) * g.axis[1]) lam = beta - alpha phi = -beta - alpha - pi @@ -34,84 +43,19 @@ def decompose_and_add(self, qubit: Qubit, angle: float, axis: Tuple[float, float phi = normalize_angle(phi) theta = normalize_angle(theta) + decomposed_g = [] + if abs(lam) > ATOL: - self.output.add_gate(rz(qubit, Float(lam))) + decomposed_g.append(rz(g.qubit, Float(lam))) - self.output.add_gate(x90(qubit)) + decomposed_g.append(x90(g.qubit)) if abs(theta) > ATOL: - self.output.add_gate(rz(qubit, Float(theta))) + decomposed_g.append(rz(g.qubit, Float(theta))) - self.output.add_gate(x90(qubit)) + decomposed_g.append(x90(g.qubit)) if abs(phi) > ATOL: - self.output.add_gate(rz(qubit, Float(phi))) - - def flush(self, q): - if q not in self.accumulated_1q_gates: - return - p = self.accumulated_1q_gates.pop(q) - self.decompose_and_add(q, p["angle"], p["axis"]) - - def flush_all(self): - while len(self.accumulated_1q_gates) > 0: - self.flush(next(iter(self.accumulated_1q_gates))) - - def accumulate(self, qubit, bloch_sphere_rotation: BlochSphereRotation): - axis, angle, phase = bloch_sphere_rotation.axis, bloch_sphere_rotation.angle, bloch_sphere_rotation.phase - - if qubit not in self.accumulated_1q_gates: - self.accumulated_1q_gates[qubit] = {"angle": angle, "axis": axis, "phase": phase} - return - - existing = self.accumulated_1q_gates[qubit] - combined_phase = phase + existing["phase"] - - a = angle - l = axis - b = existing["angle"] - m = existing["axis"] - - combined_angle = 2 * acos(cos(a / 2) * cos(b / 2) - sin(a / 2) * sin(b / 2) * np.dot(l, m)) - - if abs(sin(combined_angle / 2)) < ATOL: - self.accumulated_1q_gates.pop(qubit) - return - - combined_axis = ( - 1 - / sin(combined_angle / 2) - * (sin(a / 2) * cos(b / 2) * l + cos(a / 2) * sin(b / 2) * m + sin(a / 2) * sin(b / 2) * np.cross(l, m)) - ) - - self.accumulated_1q_gates[qubit] = {"angle": combined_angle, "axis": combined_axis, "phase": combined_phase} - - def process_gate(self, gate: Gate): - qubit_arguments = [arg for arg in gate.arguments if isinstance(arg, Qubit)] - - if len(qubit_arguments) >= 2: - [self.flush(q) for q in qubit_arguments] - self.output.add_gate(gate) - return - - if len(qubit_arguments) == 0: - assert False, "Unsupported" - return - - assert isinstance(gate, BlochSphereRotation), f"Not supported for single qubit gate `{gate.name}`: {type(gate)}" - - self.accumulate(qubit_arguments[0], gate) - - -def decompose_mckay(squirrel_ir): - impl = _McKayDecomposerImpl(squirrel_ir.number_of_qubits, squirrel_ir.qubit_register_name) - - for statement in squirrel_ir.statements: - if not isinstance(statement, Gate): - continue - - impl.process_gate(statement) - - impl.flush_all() + decomposed_g.append(rz(g.qubit, Float(phi))) - return impl.output # FIXME: instead of returning a new IR, modify existing one + return decomposed_g diff --git a/opensquirrel/merger.py b/opensquirrel/merger.py new file mode 100644 index 00000000..17309904 --- /dev/null +++ b/opensquirrel/merger.py @@ -0,0 +1,84 @@ +from math import acos, atan2, cos, pi, sin, sqrt +from typing import List, Optional + +import numpy as np + +from opensquirrel.common import ATOL +from opensquirrel.squirrel_ir import BlochSphereRotation, Gate, Qubit, SquirrelIR + + +def compose_bloch_sphere_rotations(a: BlochSphereRotation, b: BlochSphereRotation) -> BlochSphereRotation: + """Computes the Bloch sphere rotation resulting from the composition of two Bloch sphere rotations. + The first rotation is applied and then the second. + The resulting gate is anonymous except if `a` is the identity and `b` is not anonymous, or vice versa. + + Uses Rodrigues' rotation formula, see for instance https://en.wikipedia.org/wiki/Rodrigues%27_rotation_formula. + """ + assert a.qubit == b.qubit, "Cannot merge two BlochSphereRotation's on different qubits" + + combined_angle = 2 * acos( + cos(a.angle / 2) * cos(b.angle / 2) - sin(a.angle / 2) * sin(b.angle / 2) * np.dot(a.axis, b.axis) + ) + + if abs(sin(combined_angle / 2)) < ATOL: + return BlochSphereRotation.identity(a.qubit) + + combined_axis = ( + 1 + / sin(combined_angle / 2) + * ( + sin(a.angle / 2) * cos(b.angle / 2) * a.axis + + cos(a.angle / 2) * sin(b.angle / 2) * b.axis + + sin(a.angle / 2) * sin(b.angle / 2) * np.cross(a.axis, b.axis) + ) + ) + + combined_phase = a.phase + b.phase + + generator = b.generator if a.is_identity() else a.generator if b.is_identity() else None + arguments = b.arguments if a.is_identity() else a.arguments if b.is_identity() else None + + return BlochSphereRotation( + qubit=a.qubit, + axis=combined_axis, + angle=combined_angle, + phase=combined_phase, + generator=generator, + arguments=arguments, + ) + + +def merge_single_qubit_gates(squirrel_ir: SquirrelIR): + accumulators_per_qubit: dict[Qubit, BlochSphereRotation] = { + Qubit(q): BlochSphereRotation.identity(Qubit(q)) for q in range(squirrel_ir.number_of_qubits) + } + + statement_index = 0 + while statement_index < len(squirrel_ir.statements): + statement = squirrel_ir.statements[statement_index] + + if not isinstance(statement, Gate): + # Skip, since statement is not a gate + statement_index += 1 + continue + + if isinstance(statement, BlochSphereRotation): + # Accumulate + already_accumulated = accumulators_per_qubit.get(statement.qubit) + + composed = compose_bloch_sphere_rotations(statement, already_accumulated) + accumulators_per_qubit[statement.qubit] = composed + + del squirrel_ir.statements[statement_index] + continue + + for qubit_operand in statement.get_qubit_operands(): + if not accumulators_per_qubit[qubit_operand].is_identity(): + squirrel_ir.statements.insert(statement_index, accumulators_per_qubit[qubit_operand]) + accumulators_per_qubit[qubit_operand] = BlochSphereRotation.identity(qubit_operand) + statement_index += 1 + statement_index += 1 + + for accumulated_bloch_sphere_rotation in accumulators_per_qubit.values(): + if not accumulated_bloch_sphere_rotation.is_identity(): + squirrel_ir.statements.append(accumulated_bloch_sphere_rotation) diff --git a/opensquirrel/replacer.py b/opensquirrel/replacer.py index 1d9bb085..d2564b93 100644 --- a/opensquirrel/replacer.py +++ b/opensquirrel/replacer.py @@ -1,31 +1,113 @@ -from typing import List +from abc import ABC, abstractmethod +from typing import Callable, List -from opensquirrel.squirrel_ir import Comment, Gate, SquirrelIR +from opensquirrel.circuit_matrix_calculator import get_circuit_matrix +from opensquirrel.common import are_matrices_equivalent_up_to_global_phase +from opensquirrel.squirrel_ir import ( + BlochSphereRotation, + ControlledGate, + Gate, + MatrixGate, + Qubit, + SquirrelIR, + SquirrelIRVisitor, +) -def replace(squirrel_ir: SquirrelIR, gate_name_to_replace: str, f): +class Decomposer(ABC): + @staticmethod + @abstractmethod + def decompose(gate: Gate) -> [Gate]: + raise NotImplementedError() + + +class _QubitReIndexer(SquirrelIRVisitor): + def __init__(self, mappings: List[Qubit]): + self.mappings = mappings + + def visit_bloch_sphere_rotation(self, g: BlochSphereRotation): + result = BlochSphereRotation( + qubit=Qubit(self.mappings.index(g.qubit)), angle=g.angle, axis=g.axis, phase=g.phase + ) + return result + + def visit_matrix_gate(self, g: MatrixGate): + mapped_operands = [Qubit(self.mappings.index(op)) for op in g.operands] + result = MatrixGate(matrix=g.matrix, operands=mapped_operands) + return result + + def visit_controlled_gate(self, controlled_gate: ControlledGate): + control_qubit = Qubit(self.mappings.index(controlled_gate.control_qubit)) + target_gate = controlled_gate.target_gate.accept(self) + result = ControlledGate(control_qubit=control_qubit, target_gate=target_gate) + return result + + +def get_replacement_matrix(replacement: List[Gate], qubit_mappings): + replacement_ir = SquirrelIR(number_of_qubits=len(qubit_mappings), qubit_register_name="q_temp") + qubit_remapper = _QubitReIndexer(qubit_mappings) + + for gate in replacement: + gate_with_remapped_qubits = gate.accept(qubit_remapper) + replacement_ir.add_gate(gate_with_remapped_qubits) + + return get_circuit_matrix(replacement_ir) + + +def check_valid_replacement(statement, replacement): + expected_qubit_operands = statement.get_qubit_operands() + replacement_qubit_operands = set() + [replacement_qubit_operands.update(g.get_qubit_operands()) for g in replacement] + + if set(expected_qubit_operands) != replacement_qubit_operands: + raise Exception(f"Replacement for gate {statement.name} does not seem to operate on the right qubits") + + replacement_matrix = get_replacement_matrix(replacement, expected_qubit_operands) + replaced_matrix = get_replacement_matrix([statement], expected_qubit_operands) + + if not are_matrices_equivalent_up_to_global_phase(replacement_matrix, replaced_matrix): + raise Exception(f"Replacement for gate {statement.name} does not preserve the quantum state") + + +def decompose(squirrel_ir: SquirrelIR, decomposer: Decomposer): + """Applies `decomposer` to every gate in the circuit, replacing each gate by the output of `decomposer`. + When `decomposer` decides to not decompose a gate, it needs to return a list with the intact gate as single element. + """ statement_index = 0 while statement_index < len(squirrel_ir.statements): statement = squirrel_ir.statements[statement_index] - if isinstance(statement, Comment): - statement_index += 1 - continue - if not isinstance(statement, Gate): - raise Exception("Unsupported") - - if statement.name != gate_name_to_replace: statement_index += 1 continue - # FIXME: handle case where if f is not a function but directly a list. + replacement: List[Gate] = decomposer.decompose(statement) + + check_valid_replacement(statement, replacement) - replacement: List[Gate] = f(*statement.arguments) squirrel_ir.statements[statement_index : statement_index + 1] = replacement statement_index += len(replacement) - # TODO: Here, check that the semantic of the replacement is the same! - # For this, need to update the simulation capabilities. - # TODO: Do we allow skipping the replacement, based on arguments? +class _GenericReplacer(Decomposer): + def __init__(self, gate_generator: Callable[..., Gate], replacement_function): + self.gate_generator = gate_generator + self.replacement_function = replacement_function + + def decompose(self, g: Gate) -> [Gate]: + if g.is_anonymous or g.generator != self.gate_generator: + return [g] + return self.replacement_function(*g.arguments) + + +def replace(squirrel_ir: SquirrelIR, gate_generator: Callable[..., Gate], f): + """Does the same as decomposer, but only applies to a given gate.""" + + def generic_replacer(g: Gate) -> [Gate]: + if g.is_anonymous or g.generator != gate_generator: + return [g] + return f(*g.arguments) + + generic_replacer = _GenericReplacer(gate_generator, f) + + decompose(squirrel_ir, generic_replacer) diff --git a/opensquirrel/squirrel_ir.py b/opensquirrel/squirrel_ir.py index 2d41e68e..3455a6d9 100644 --- a/opensquirrel/squirrel_ir.py +++ b/opensquirrel/squirrel_ir.py @@ -1,12 +1,14 @@ +import cmath import inspect -from abc import ABC +from abc import ABC, abstractmethod from dataclasses import dataclass from functools import wraps +from math import cos, sin from typing import Callable, List, Optional, Tuple import numpy as np -from opensquirrel.common import ATOL, normalize_angle, normalize_axis +from opensquirrel.common import ATOL, X, Y, Z, normalize_angle, normalize_axis class SquirrelIRVisitor(ABC): @@ -36,6 +38,7 @@ def visit_controlled_gate(self, controlled_gate: "ControlledGate"): class IRNode(ABC): + @abstractmethod def accept(self, visitor: SquirrelIRVisitor): pass @@ -79,8 +82,13 @@ class Statement(IRNode, ABC): class Gate(Statement, ABC): - arguments: Optional[Tuple[Expression, ...]] = None + # Note: two gates are considered equal even when their generators/arguments are different. generator: Optional[Callable[..., "Gate"]] = None + arguments: Optional[Tuple[Expression, ...]] = None + + def __init__(self, generator, arguments): + self.generator = generator + self.arguments = arguments @property def name(self) -> Optional[str]: @@ -90,20 +98,40 @@ def name(self) -> Optional[str]: def is_anonymous(self) -> bool: return self.arguments is None + @abstractmethod + def get_qubit_operands(self) -> List[Qubit]: + raise NotImplementedError + class BlochSphereRotation(Gate): generator: Optional[Callable[..., "BlochSphereRotation"]] = None - def __init__(self, qubit: Qubit, axis: Tuple[float, float, float], angle: float, phase: float = 0): + def __init__( + self, + qubit: Qubit, + axis: Tuple[float, float, float], + angle: float, + phase: float = 0, + generator=None, + arguments=None, + ): + Gate.__init__(self, generator, arguments) self.qubit: Qubit = qubit self.axis = normalize_axis(np.array(axis).astype(np.float64)) self.angle = normalize_angle(angle) self.phase = normalize_angle(phase) + @staticmethod + def identity(q: Qubit) -> "BlochSphereRotation": + return BlochSphereRotation(qubit=q, axis=(1, 0, 0), angle=0, phase=0) + def __repr__(self): return f"BlochSphereRotation({self.qubit}, axis={self.axis}, angle={self.angle}, phase={self.phase})" def __eq__(self, other): + if not isinstance(other, BlochSphereRotation): + return False + if self.qubit != other.qubit: return False @@ -120,11 +148,19 @@ def accept(self, visitor: SquirrelIRVisitor): visitor.visit_gate(self) return visitor.visit_bloch_sphere_rotation(self) + def get_qubit_operands(self) -> List[Qubit]: + return [self.qubit] + + def is_identity(self) -> bool: + # Angle and phase are already normalized. + return abs(self.angle) < ATOL and abs(self.phase) < ATOL + class MatrixGate(Gate): generator: Optional[Callable[..., "MatrixGate"]] = None - def __init__(self, matrix: np.ndarray, operands: List[Qubit]): + def __init__(self, matrix: np.ndarray, operands: List[Qubit], generator=None, arguments=None): + Gate.__init__(self, generator, arguments) assert len(operands) >= 2, "For 1q gates, please use BlochSphereRotation" assert matrix.shape == (1 << len(operands), 1 << len(operands)) @@ -133,6 +169,8 @@ def __init__(self, matrix: np.ndarray, operands: List[Qubit]): def __eq__(self, other): # TODO: Determine whether we shall allow for a global phase difference here. + if not isinstance(other, MatrixGate): + return False # FIXME: a MatrixGate can hide a ControlledGate. https://github.com/QuTech-Delft/OpenSquirrel/issues/88 return np.allclose(self.matrix, other.matrix) def __repr__(self): @@ -142,15 +180,21 @@ def accept(self, visitor: SquirrelIRVisitor): visitor.visit_gate(self) return visitor.visit_matrix_gate(self) + def get_qubit_operands(self) -> List[Qubit]: + return self.operands + class ControlledGate(Gate): generator: Optional[Callable[..., "ControlledGate"]] = None - def __init__(self, control_qubit: Qubit, target_gate: Gate): + def __init__(self, control_qubit: Qubit, target_gate: Gate, generator=None, arguments=None): + Gate.__init__(self, generator, arguments) self.control_qubit = control_qubit self.target_gate = target_gate def __eq__(self, other): + if not isinstance(other, ControlledGate): + return False # FIXME: a MatrixGate can hide a ControlledGate. https://github.com/QuTech-Delft/OpenSquirrel/issues/88 if self.control_qubit != other.control_qubit: return False @@ -163,16 +207,19 @@ def accept(self, visitor: SquirrelIRVisitor): visitor.visit_gate(self) return visitor.visit_controlled_gate(self) + def get_qubit_operands(self) -> List[Qubit]: + return [self.control_qubit] + self.target_gate.get_qubit_operands() + def named_gate(gate_generator: Callable[..., Gate]) -> Callable[..., Gate]: @wraps(gate_generator) - def wrapper(*args): + def wrapper(*args, **kwargs): for i, par in enumerate(inspect.signature(gate_generator).parameters.values()): if not issubclass(par.annotation, Expression): raise TypeError("Gate argument types must be expressions") - result = gate_generator(*args) - result.generator = gate_generator + result = gate_generator(*args, **kwargs) + result.generator = wrapper result.arguments = args return result diff --git a/test/ir_equality_test_base.py b/test/ir_equality_test_base.py new file mode 100644 index 00000000..84ce8f24 --- /dev/null +++ b/test/ir_equality_test_base.py @@ -0,0 +1,38 @@ +import unittest + +import numpy as np + +from opensquirrel import circuit_matrix_calculator +from opensquirrel.common import ATOL, are_matrices_equivalent_up_to_global_phase + + +class IREqualityTestBase(unittest.TestCase): + def check_equivalence_up_to_global_phase(self, matrix_a, matrix_b): + self.assertTrue(are_matrices_equivalent_up_to_global_phase(matrix_a, matrix_b)) + + def modify_ir_and_check(self, ir, action, expected_ir=None): + """ + Checks whether the IR action preserves: + - the number of qubits, + - the qubit register name(s), + - the circuit matrix up to a global phase factor. + """ + + # Store matrix before decompositions. + expected_matrix = circuit_matrix_calculator.get_circuit_matrix(ir) + + expected_number_of_qubits = ir.number_of_qubits + expected_qubit_register_name = ir.qubit_register_name + + action(ir) + + self.assertEqual(ir.number_of_qubits, expected_number_of_qubits) + self.assertEqual(ir.qubit_register_name, expected_qubit_register_name) + + # Get matrix after decompositions. + actual_matrix = circuit_matrix_calculator.get_circuit_matrix(ir) + + self.check_equivalence_up_to_global_phase(actual_matrix, expected_matrix) + + if expected_ir is not None: + self.assertEqual(ir, expected_ir) diff --git a/test/test_integration.py b/test/test_integration.py index 48a253d3..1bbcaedb 100644 --- a/test/test_integration.py +++ b/test/test_integration.py @@ -4,6 +4,7 @@ from opensquirrel.circuit import Circuit from opensquirrel.default_gates import * +from opensquirrel.mckay_decomposer import McKayDecomposer class IntegrationTest(unittest.TestCase): @@ -30,7 +31,7 @@ def test_simple(self): # myCircuit.replace( - "cnot", + cnot, lambda control, target: [ h(target), cz(control, target), @@ -40,7 +41,9 @@ def test_simple(self): # Do 1q-gate fusion and decompose with McKay decomposition. - myCircuit.decompose_mckay() + myCircuit.merge_single_qubit_gates() + + myCircuit.decompose(decomposer=McKayDecomposer) # Write the transformed circuit as a cQasm3 string. @@ -62,16 +65,16 @@ def test_simple(self): x90 qreg[1] rz qreg[1], 3.1415927 cz qreg[0], qreg[1] -rz qreg[1], 3.1415927 -x90 qreg[1] -rz qreg[1], 1.572389 -x90 qreg[1] -rz qreg[1], 3.1415927 rz qreg[0], 1.5707963 x90 qreg[0] rz qreg[0], 0.84159265 x90 qreg[0] rz qreg[0], 1.5707963 +rz qreg[1], 3.1415927 +x90 qreg[1] +rz qreg[1], 1.572389 +x90 qreg[1] +rz qreg[1], 3.1415927 """, ) @@ -107,7 +110,9 @@ def test_qi(self): """ ) - myCircuit.decompose_mckay() + myCircuit.merge_single_qubit_gates() + + myCircuit.decompose(decomposer=McKayDecomposer) output = str(myCircuit) expected = """version 3.0 diff --git a/test/test_mckay_decomposer.py b/test/test_mckay_decomposer.py index daee55d3..155ab1dc 100644 --- a/test/test_mckay_decomposer.py +++ b/test/test_mckay_decomposer.py @@ -1,91 +1,71 @@ import unittest +from test.ir_equality_test_base import IREqualityTestBase -from opensquirrel import circuit_matrix_calculator, mckay_decomposer from opensquirrel.default_gates import * +from opensquirrel.mckay_decomposer import McKayDecomposer from opensquirrel.squirrel_ir import Float, Qubit, SquirrelIR -class DecomposeMcKayTests(unittest.TestCase): - def check_equivalence_up_to_global_phase(self, matrix_a, matrix_b): - first_non_zero = next( - (i, j) for i in range(matrix_a.shape[0]) for j in range(matrix_a.shape[1]) if abs(matrix_a[i, j]) > ATOL +class DecomposeMcKayTests(IREqualityTestBase): + def test_ignores_2q_gates(self): + self.assertEqual(McKayDecomposer.decompose(cnot(Qubit(0), Qubit(1))), [cnot(Qubit(0), Qubit(1))]) + self.assertEqual( + McKayDecomposer.decompose(cr(Qubit(2), Qubit(3), Float(2.123))), [cr(Qubit(2), Qubit(3), Float(2.123))] ) - if abs(matrix_b[first_non_zero]) < ATOL: - return False - - phase_difference = matrix_a[first_non_zero] / matrix_b[first_non_zero] - - self.assertTrue(np.allclose(matrix_a, phase_difference * matrix_b)) - - def check_mckay_decomposition(self, squirrel_ir, expected_ir=None): - """ - Check whether the mcKay decomposition transformation applied to the input IR preserves the - circuit matrix up to a global phase factor. - """ - - # Store matrix before decompositions. - expected_matrix = circuit_matrix_calculator.get_circuit_matrix(squirrel_ir) - - output = mckay_decomposer.decompose_mckay(squirrel_ir) - - self.assertEqual(output.number_of_qubits, squirrel_ir.number_of_qubits) - self.assertEqual(output.qubit_register_name, squirrel_ir.qubit_register_name) - - if expected_ir is not None: - self.assertEqual(output, expected_ir) - - # Get matrix after decompositions. - actual_matrix = circuit_matrix_calculator.get_circuit_matrix(output) - - self.check_equivalence_up_to_global_phase(actual_matrix, expected_matrix) - - def test_one(self): - ir = SquirrelIR(number_of_qubits=2, qubit_register_name="squirrel") - - ir.add_gate(ry(Qubit(0), Float(23847628349.123))) - ir.add_gate(rx(Qubit(0), Float(29384672.234))) - ir.add_gate(rz(Qubit(0), Float(9877.87634))) - - self.check_mckay_decomposition(ir) - - def test_two(self): - ir = SquirrelIR(number_of_qubits=2, qubit_register_name="squirrel") - - ir.add_gate(ry(Qubit(0), Float(23847628349.123))) - ir.add_gate(cnot(Qubit(0), Qubit(1))) - ir.add_gate(rx(Qubit(0), Float(29384672.234))) - ir.add_gate(rz(Qubit(0), Float(9877.87634))) - ir.add_gate(cnot(Qubit(0), Qubit(1))) - ir.add_gate(rx(Qubit(0), Float(29384672.234))) - ir.add_gate(rz(Qubit(0), Float(9877.87634))) - - self.check_mckay_decomposition(ir) - - def test_small_random(self): - ir = SquirrelIR(number_of_qubits=4, qubit_register_name="q") + def test_identity_empty_decomposition(self): + self.assertEqual(McKayDecomposer.decompose(BlochSphereRotation.identity(Qubit(0))), []) + + def test_x(self): + self.assertEqual( + McKayDecomposer.decompose(x(Qubit(0))), + [ + # FIXME: we can do better here. See https://github.com/QuTech-Delft/OpenSquirrel/issues/89. + rz(Qubit(0), Float(-math.pi / 2)), + x90(Qubit(0)), + x90(Qubit(0)), + rz(Qubit(0), Float(-math.pi / 2)), + ], + ) - ir.add_gate(h(Qubit(2))) - ir.add_gate(cr(Qubit(2), Qubit(3), Float(2.123))) - ir.add_gate(h(Qubit(1))) - ir.add_gate(h(Qubit(0))) - ir.add_gate(h(Qubit(2))) - ir.add_gate(h(Qubit(1))) - ir.add_gate(h(Qubit(0))) - ir.add_gate(cr(Qubit(2), Qubit(3), Float(2.123))) + def test_y(self): + self.assertEqual( + McKayDecomposer.decompose(y(Qubit(0))), [rz(Qubit(0), Float(math.pi)), x90(Qubit(0)), x90(Qubit(0))] + ) - expected_ir = SquirrelIR(number_of_qubits=4, qubit_register_name="q") + def test_z(self): + self.assertEqual( + McKayDecomposer.decompose(z(Qubit(0))), + [ + rz(Qubit(0), Float(-math.pi / 2)), + x90(Qubit(0)), + rz(Qubit(0), Float(math.pi)), + x90(Qubit(0)), + rz(Qubit(0), Float(math.pi / 2)), + ], + ) - expected_ir.add_gate(x90(Qubit(2))) - expected_ir.add_gate(rz(Qubit(2), Float(1.5707963267948966))) - expected_ir.add_gate(x90(Qubit(2))) - expected_ir.add_gate(cr(Qubit(2), Qubit(3), Float(2.123))) - expected_ir.add_gate(x90(Qubit(2))) - expected_ir.add_gate(rz(Qubit(2), Float(1.5707963267948966))) - expected_ir.add_gate(x90(Qubit(2))) - expected_ir.add_gate(cr(Qubit(2), Qubit(3), Float(2.123))) + def test_hadamard(self): + self.assertEqual( + McKayDecomposer.decompose(h(Qubit(0))), + [ + BlochSphereRotation(Qubit(0), axis=(1, 0, 0), angle=math.pi / 2, phase=0.0), + BlochSphereRotation(Qubit(0), axis=(0, 0, 1), angle=math.pi / 2, phase=0.0), + BlochSphereRotation(Qubit(0), axis=(1, 0, 0), angle=math.pi / 2, phase=0.0), + ], + ) - self.check_mckay_decomposition(ir, expected_ir) + def test_arbitrary(self): + self.assertEqual( + McKayDecomposer.decompose(BlochSphereRotation(qubit=Qubit(0), angle=5.21, axis=(1, 2, 3), phase=0.324)), + [ + rz(Qubit(0), Float(0.018644578210707863)), + x90(Qubit(0)), + rz(Qubit(0), Float(2.520651583905213)), + x90(Qubit(0)), + rz(Qubit(0), Float(2.2329420137988887)), + ], + ) if __name__ == "__main__": diff --git a/test/test_merger.py b/test/test_merger.py new file mode 100644 index 00000000..71bcc7b5 --- /dev/null +++ b/test/test_merger.py @@ -0,0 +1,108 @@ +import unittest +from test.ir_equality_test_base import IREqualityTestBase + +from opensquirrel import merger +from opensquirrel.default_gates import * +from opensquirrel.merger import compose_bloch_sphere_rotations +from opensquirrel.squirrel_ir import Float, Qubit, SquirrelIR + + +class MergerTest(IREqualityTestBase): + def test_compose_bloch_sphere_rotations_same_axis(self): + q = Qubit(123) + a = BlochSphereRotation(qubit=q, axis=(1, 2, 3), angle=0.4) + b = BlochSphereRotation(qubit=q, axis=(1, 2, 3), angle=-0.3) + composed = compose_bloch_sphere_rotations(a, b) + self.assertEqual(composed, BlochSphereRotation(qubit=q, axis=(1, 2, 3), angle=0.1)) + + def test_compose_bloch_sphere_rotations_different_axis(self): + # Visualizing this in 3D is difficult... + q = Qubit(123) + a = BlochSphereRotation(qubit=q, axis=(1, 0, 0), angle=math.pi / 2) + b = BlochSphereRotation(qubit=q, axis=(0, 0, 1), angle=-math.pi / 2) + c = BlochSphereRotation(qubit=q, axis=(0, 1, 0), angle=math.pi / 2) + composed = compose_bloch_sphere_rotations(a, compose_bloch_sphere_rotations(b, c)) + self.assertEqual(composed, BlochSphereRotation(qubit=q, axis=(1, 1, 0), angle=math.pi)) + + def test_single_gate(self): + ir = SquirrelIR(number_of_qubits=2, qubit_register_name="q") + ir.add_gate(ry(Qubit(0), Float(1.2345))) + + expected_ir = SquirrelIR(number_of_qubits=2, qubit_register_name="q") + expected_ir.add_gate(ry(Qubit(0), Float(1.2345))) + + self.modify_ir_and_check(ir, action=merger.merge_single_qubit_gates, expected_ir=expected_ir) + + # Check that when no fusion happens, generator and arguments of gates are preserved. + self.assertEqual(ir.statements[0].generator, ry) + self.assertEqual(ir.statements[0].arguments, (Qubit(0), Float(1.2345))) + + def test_two_hadamards(self): + ir = SquirrelIR(number_of_qubits=4, qubit_register_name="q") + + ir.add_gate(h(Qubit(2))) + ir.add_gate(h(Qubit(2))) + + expected_ir = SquirrelIR(number_of_qubits=4, qubit_register_name="q") + + self.modify_ir_and_check(ir, action=merger.merge_single_qubit_gates, expected_ir=expected_ir) + + def test_two_hadamards_different_qubits(self): + ir = SquirrelIR(number_of_qubits=4, qubit_register_name="q") + ir.add_gate(h(Qubit(0))) + ir.add_gate(h(Qubit(2))) + + expected_ir = SquirrelIR(number_of_qubits=4, qubit_register_name="q") + expected_ir.add_gate(h(Qubit(0))) + expected_ir.add_gate(h(Qubit(2))) + + self.modify_ir_and_check(ir, action=merger.merge_single_qubit_gates, expected_ir=expected_ir) + + def test_merge_different_qubits(self): + ir = SquirrelIR(number_of_qubits=4, qubit_register_name="q") + ir.add_gate(ry(Qubit(0), Float(math.pi / 2))) + ir.add_gate(rx(Qubit(0), Float(math.pi))) + ir.add_gate(rz(Qubit(1), Float(1.2345))) + ir.add_gate(ry(Qubit(2), Float(1))) + ir.add_gate(ry(Qubit(2), Float(3.234))) + + expected_ir = SquirrelIR(number_of_qubits=4, qubit_register_name="q") + expected_ir.add_gate( + BlochSphereRotation(qubit=Qubit(0), axis=(1, 0, 1), angle=math.pi) + ) # This is hadamard with 0 phase... + expected_ir.add_gate(rz(Qubit(1), Float(1.2345))) + expected_ir.add_gate(ry(Qubit(2), Float(4.234))) + + self.modify_ir_and_check(ir, action=merger.merge_single_qubit_gates, expected_ir=expected_ir) + + self.assertTrue(ir.statements[0].is_anonymous) # When fusion happens, the resulting gate is anonymous. + self.assertEqual(ir.statements[1].generator, rz) # Otherwise it keeps the same generator and arguments. + self.assertEqual(ir.statements[1].arguments, (Qubit(1), Float(1.2345))) + self.assertTrue(ir.statements[2].is_anonymous) + + def test_merge_and_flush(self): + ir = SquirrelIR(number_of_qubits=4, qubit_register_name="q") + ir.add_gate(ry(Qubit(0), Float(math.pi / 2))) + ir.add_gate(rz(Qubit(1), Float(1.5))) + ir.add_gate(rx(Qubit(0), Float(math.pi))) + ir.add_gate(rz(Qubit(1), Float(-2.5))) + ir.add_gate(cnot(Qubit(0), Qubit(1))) + ir.add_gate(ry(Qubit(0), Float(3.234))) + + expected_ir = SquirrelIR(number_of_qubits=4, qubit_register_name="q") + expected_ir.add_gate( + BlochSphereRotation(qubit=Qubit(0), axis=(1, 0, 1), angle=math.pi) + ) # This is hadamard with 0 phase... + expected_ir.add_gate(rz(Qubit(1), Float(-1.0))) + expected_ir.add_gate(cnot(Qubit(0), Qubit(1))) + expected_ir.add_gate(ry(Qubit(0), Float(3.234))) + + self.modify_ir_and_check(ir, action=merger.merge_single_qubit_gates, expected_ir=expected_ir) + + self.assertTrue(ir.statements[0].is_anonymous) + self.assertEqual(ir.statements[3].generator, ry) + self.assertEqual(ir.statements[3].arguments, (Qubit(0), Float(3.234))) + + +if __name__ == "__main__": + unittest.main() diff --git a/test/test_replacer.py b/test/test_replacer.py index f3c67a13..13e0d617 100644 --- a/test/test_replacer.py +++ b/test/test_replacer.py @@ -2,24 +2,135 @@ from opensquirrel import replacer from opensquirrel.default_gates import * +from opensquirrel.replacer import Decomposer, check_valid_replacement from opensquirrel.squirrel_ir import Qubit, SquirrelIR -def hadamard_decomposition(q: Qubit): - return [ - y90(q), - x(q), - ] +class ReplacerTest(unittest.TestCase): + @staticmethod + def test_check_valid_replacement(): + check_valid_replacement(BlochSphereRotation.identity(Qubit(0)), [BlochSphereRotation.identity(Qubit(0))]) + check_valid_replacement( + BlochSphereRotation.identity(Qubit(0)), + [BlochSphereRotation.identity(Qubit(0)), BlochSphereRotation.identity(Qubit(0))], + ) + + check_valid_replacement(BlochSphereRotation.identity(Qubit(0)), [h(Qubit(0)), h(Qubit(0))]) + + check_valid_replacement(h(Qubit(0)), [h(Qubit(0)), h(Qubit(0)), h(Qubit(0))]) + + check_valid_replacement( + cnot(Qubit(0), Qubit(1)), [cnot(Qubit(0), Qubit(1)), BlochSphereRotation.identity(Qubit(0))] + ) + + # Arbitrary global phase change is not considered an issue. + check_valid_replacement( + cnot(Qubit(0), Qubit(1)), + [cnot(Qubit(0), Qubit(1)), BlochSphereRotation(Qubit(0), angle=0, axis=(1, 0, 0), phase=621.6546)], + ) + + def test_check_valid_replacement_wrong_qubit(self): + with self.assertRaisesRegex(Exception, r"Replacement for gate h does not seem to operate on the right qubits"): + check_valid_replacement(h(Qubit(0)), [h(Qubit(1))]) + + with self.assertRaisesRegex( + Exception, r"Replacement for gate cnot does not seem to operate on the right qubits" + ): + check_valid_replacement(cnot(Qubit(0), Qubit(1)), [cnot(Qubit(2), Qubit(1))]) + + with self.assertRaisesRegex(Exception, r"Replacement for gate cnot does not preserve the quantum state"): + check_valid_replacement(cnot(Qubit(0), Qubit(1)), [cnot(Qubit(1), Qubit(0))]) + + def test_check_valid_replacement_cnot_as_sqrt_swap(self): + # https://en.wikipedia.org/wiki/Quantum_logic_gate#/media/File:Qcircuit_CNOTsqrtSWAP2.svg + c = Qubit(0) + t = Qubit(1) + check_valid_replacement( + cnot(control=c, target=t), + [ + ry(t, Float(math.pi / 2)), + sqrt_swap(c, t), + z(c), + sqrt_swap(c, t), + rz(c, Float(-math.pi / 2)), + rz(t, Float(-math.pi / 2)), + ry(t, Float(-math.pi / 2)), + ], + ) + + with self.assertRaisesRegex(Exception, r"Replacement for gate cnot does not preserve the quantum state"): + check_valid_replacement( + cnot(control=c, target=t), + [ + ry(t, Float(math.pi / 2)), + sqrt_swap(c, t), + z(c), + sqrt_swap(c, t), + rz(c, Float(-math.pi / 2 + 0.01)), + rz(t, Float(-math.pi / 2)), + ry(t, Float(-math.pi / 2)), + ], + ) + + with self.assertRaisesRegex( + Exception, r"Replacement for gate cnot does not seem to operate on the right qubits" + ): + check_valid_replacement( + cnot(control=c, target=t), + [ + ry(t, Float(math.pi / 2)), + sqrt_swap(c, t), + z(c), + sqrt_swap(c, Qubit(2)), + rz(c, Float(-math.pi / 2 + 0.01)), + rz(t, Float(-math.pi / 2)), + ry(t, Float(-math.pi / 2)), + ], + ) + + def test_check_valid_replacement_large_number_of_qubits(self): + # If we were building the whole circuit matrix, this would run out of memory. + check_valid_replacement(h(Qubit(9234687)), [y90(Qubit(9234687)), x(Qubit(9234687))]) + + with self.assertRaisesRegex(Exception, r"Replacement for gate h does not seem to operate on the right qubits"): + check_valid_replacement(h(Qubit(9234687)), [y90(Qubit(698446519)), x(Qubit(9234687))]) + + with self.assertRaisesRegex(Exception, r"Replacement for gate h does not preserve the quantum state"): + check_valid_replacement(h(Qubit(9234687)), [y90(Qubit(9234687)), x(Qubit(9234687)), x(Qubit(9234687))]) + + def test_replace_generic(self): + squirrel_ir = SquirrelIR(number_of_qubits=3, qubit_register_name="test") + + squirrel_ir.add_gate(h(Qubit(0))) + squirrel_ir.add_gate(cnot(Qubit(0), Qubit(1))) + + # A stupid replacer function that adds identities before and after single-qubit gates. + class TestDecomposer(Decomposer): + def decompose(self, g: Gate) -> [Gate]: + if isinstance(g, BlochSphereRotation): + return [BlochSphereRotation.identity(g.qubit), g, BlochSphereRotation.identity(g.qubit)] + + return [g] + + replacer.decompose(squirrel_ir, decomposer=TestDecomposer()) + + expected_ir = SquirrelIR(number_of_qubits=3, qubit_register_name="test") + + expected_ir.add_gate(BlochSphereRotation.identity(Qubit(0))) + expected_ir.add_gate(h(Qubit(0))) + expected_ir.add_gate(BlochSphereRotation.identity(Qubit(0))) + expected_ir.add_gate(cnot(Qubit(0), Qubit(1))) + + self.assertEqual(expected_ir, squirrel_ir) -class ReplacerTest(unittest.TestCase): def test_replace(self): squirrel_ir = SquirrelIR(number_of_qubits=3, qubit_register_name="test") squirrel_ir.add_gate(h(Qubit(0))) squirrel_ir.add_comment(Comment("Test comment.")) - replacer.replace(squirrel_ir, "h", hadamard_decomposition) + replacer.replace(squirrel_ir, h, lambda q: [y90(q), x(q)]) expected_ir = SquirrelIR(number_of_qubits=3, qubit_register_name="test")