Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Optimal Synthesis into RXX(theta) #6551

Merged
merged 20 commits into from
Nov 3, 2021
Merged
Changes from 16 commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
17 changes: 17 additions & 0 deletions qiskit/quantum_info/synthesis/xx_decompose/__init__.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,17 @@
# This code is part of Qiskit.
#
# (C) Copyright IBM 2021
#
# This code is licensed under the Apache License, Version 2.0. You may
# obtain a copy of this license in the LICENSE.txt file in the root directory
# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
#
# Any modifications or derivative works of this code must retain this
# copyright notice, and modified files need to carry a notice indicating
# that they have been altered from the originals.

"""
Submodule symbol exports for XX decomposition.
"""

from .decomposer import XXDecomposer
298 changes: 298 additions & 0 deletions qiskit/quantum_info/synthesis/xx_decompose/circuits.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,298 @@
# This code is part of Qiskit.
#
# (C) Copyright IBM 2021
#
# This code is licensed under the Apache License, Version 2.0. You may
# obtain a copy of this license in the LICENSE.txt file in the root directory
# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
#
# Any modifications or derivative works of this code must retain this
# copyright notice, and modified files need to carry a notice indicating
# that they have been altered from the originals.

"""
Tools for building optimal circuits out of XX interactions.
Inputs:
+ A set of native XX operations, described as strengths.
+ A right-angled path, computed using the methods in `paths.py`.
Output:
+ A circuit which implements the target operation (expressed exactly as the exponential of
`a XX + b YY + c ZZ`) using the native operations and local gates.
"""

from functools import reduce
import math
from operator import itemgetter

import numpy as np

from qiskit.circuit.quantumcircuit import QuantumCircuit
from qiskit.circuit.library.standard_gates import RXXGate, RYYGate, RZGate
from qiskit.exceptions import QiskitError

from .paths import decomposition_hop
from .utilities import EPSILON, safe_arccos
from .weyl import (
apply_reflection,
apply_shift,
canonical_rotation_circuit,
reflection_options,
shift_options,
)


# pylint:disable=invalid-name
def decompose_xxyy_into_xxyy_xx(a_target, b_target, a_source, b_source, interaction):
"""
Consumes a target canonical interaction CAN(a_target, b_target) and source interactions
CAN(a1, b1), CAN(a2), then manufactures a circuit identity of the form
CAN(a_target, b_target) = (Zr, Zs) CAN(a_source, b_source) (Zu, Zv) CAN(interaction) (Zx, Zy).
Returns the 6-tuple (r, s, u, v, x, y).
"""

cplus, cminus = np.cos(a_source + b_source), np.cos(a_source - b_source)
splus, sminus = np.sin(a_source + b_source), np.sin(a_source - b_source)
ca, sa = np.cos(interaction), np.sin(interaction)

uplusv = (
1
/ 2
* safe_arccos(
cminus ** 2 * ca ** 2 + sminus ** 2 * sa ** 2 - np.cos(a_target - b_target) ** 2,
2 * cminus * ca * sminus * sa,
)
)
uminusv = (
1
/ 2
* safe_arccos(
cplus ** 2 * ca ** 2 + splus ** 2 * sa ** 2 - np.cos(a_target + b_target) ** 2,
2 * cplus * ca * splus * sa,
)
)

u, v = (uplusv + uminusv) / 2, (uplusv - uminusv) / 2

# NOTE: the target matrix is phase-free
middle_matrix = reduce(
np.dot,
[
RXXGate(2 * a_source).to_matrix() @ RYYGate(2 * b_source).to_matrix(),
np.kron(RZGate(2 * u).to_matrix(), RZGate(2 * v).to_matrix()),
RXXGate(2 * interaction).to_matrix(),
],
)

phase_solver = np.array(
[
[
1 / 4,
1 / 4,
1 / 4,
1 / 4,
],
[
1 / 4,
-1 / 4,
-1 / 4,
1 / 4,
],
[
1 / 4,
1 / 4,
-1 / 4,
-1 / 4,
],
[
1 / 4,
-1 / 4,
1 / 4,
-1 / 4,
],
]
)
inner_phases = [
np.angle(middle_matrix[0, 0]),
np.angle(middle_matrix[1, 1]),
np.angle(middle_matrix[1, 2]) + np.pi / 2,
np.angle(middle_matrix[0, 3]) + np.pi / 2,
]
r, s, x, y = np.dot(phase_solver, inner_phases)

# If there's a phase discrepancy, need to conjugate by an extra Z/2 (x) Z/2.
generated_matrix = reduce(
np.dot,
[
np.kron(RZGate(2 * r).to_matrix(), RZGate(2 * s).to_matrix()),
middle_matrix,
np.kron(RZGate(2 * x).to_matrix(), RZGate(2 * y).to_matrix()),
],
)
if (abs(np.angle(generated_matrix[3, 0]) - np.pi / 2) < 0.01 and a_target > b_target) or (
abs(np.angle(generated_matrix[3, 0]) + np.pi / 2) < 0.01 and a_target < b_target
):
x += np.pi / 4
y += np.pi / 4
r -= np.pi / 4
s -= np.pi / 4

return r, s, u, v, x, y


def xx_circuit_step(source, strength, target, embodiment):
"""
Builds a single step in an XX-based circuit.
`source` and `target` are positive canonical coordinates; `strength` is the interaction strength
at this step in the circuit as a canonical coordinate (so that CX = RZX(pi/2) corresponds to
pi/4); and `embodiment` is a Qiskit circuit which enacts the canonical gate of the prescribed
interaction `strength`.
"""

permute_source_for_overlap, permute_target_for_overlap = None, None

# apply all possible reflections, shifts to the source
for source_reflection_name in reflection_options:
reflected_source_coord, source_reflection, reflection_phase_shift = apply_reflection(
source_reflection_name, source
)
for source_shift_name in shift_options:
shifted_source_coord, source_shift, shift_phase_shift = apply_shift(
source_shift_name, reflected_source_coord
)

# check for overlap, back out permutation
source_shared, target_shared = None, None
for i, j in [(0, 0), (0, 1), (0, 2), (1, 0), (1, 1), (1, 2), (2, 0), (2, 1), (2, 2)]:

if (
abs(np.mod(abs(shifted_source_coord[i] - target[j]), np.pi)) < EPSILON
or abs(np.mod(abs(shifted_source_coord[i] - target[j]), np.pi) - np.pi)
< EPSILON
):
source_shared, target_shared = i, j
break
if source_shared is None:
continue

# pick out the other coordinates
source_first, source_second = [x for x in [0, 1, 2] if x != source_shared]
target_first, target_second = [x for x in [0, 1, 2] if x != target_shared]

# check for arccos validity
r, s, u, v, x, y = decompose_xxyy_into_xxyy_xx(
float(target[target_first]),
float(target[target_second]),
float(shifted_source_coord[source_first]),
float(shifted_source_coord[source_second]),
float(strength),
)
if any(math.isnan(val) for val in (r, s, u, v, x, y)):
continue

# OK: this combination of things works.
# save the permutation which rotates the shared coordinate into ZZ.
permute_source_for_overlap = canonical_rotation_circuit(source_first, source_second)
permute_target_for_overlap = canonical_rotation_circuit(target_first, target_second)
break

if permute_source_for_overlap is not None:
break

if permute_source_for_overlap is None:
raise QiskitError(
f"Error during RZX decomposition: Could not find a suitable Weyl "
f"reflection to match {source} to {target} along {strength}."
)

prefix_circuit, affix_circuit = QuantumCircuit(2), QuantumCircuit(2)

# the basic formula we're trying to work with is:
# target^p_t_f_o =
# rs * (source^s_reflection * s_shift)^p_s_f_o * uv * operation * xy
# but we're rearranging it into the form
# target = affix source prefix
# and computing just the prefix / affix circuits.

# the outermost prefix layer comes from the (inverse) target permutation.
prefix_circuit.compose(permute_target_for_overlap.inverse(), inplace=True)
# the middle prefix layer comes from the local Z rolls.
prefix_circuit.rz(2 * x, [0])
prefix_circuit.rz(2 * y, [1])
prefix_circuit.compose(embodiment, inplace=True)
prefix_circuit.rz(2 * u, [0])
prefix_circuit.rz(2 * v, [1])
# the innermost prefix layer is source_reflection, shifted by source_shift,
# finally conjugated by p_s_f_o.
prefix_circuit.compose(permute_source_for_overlap, inplace=True)
prefix_circuit.compose(source_reflection, inplace=True)
prefix_circuit.global_phase += -np.log(reflection_phase_shift).imag
prefix_circuit.global_phase += -np.log(shift_phase_shift).imag

# the affix circuit is constructed in reverse.
# first (i.e., innermost), we install the other half of the source transformations and p_s_f_o.
affix_circuit.compose(source_reflection.inverse(), inplace=True)
affix_circuit.compose(source_shift, inplace=True)
affix_circuit.compose(permute_source_for_overlap.inverse(), inplace=True)
# then, the other local rolls in the middle.
affix_circuit.rz(2 * r, [0])
affix_circuit.rz(2 * s, [1])
# finally, the other half of the p_t_f_o conjugation.
affix_circuit.compose(permute_target_for_overlap, inplace=True)

return {"prefix_circuit": prefix_circuit, "affix_circuit": affix_circuit}


def canonical_xx_circuit(target, strength_sequence, basis_embodiments):
"""
Assembles a Qiskit circuit from a specified `strength_sequence` of XX-type interactions which
emulates the canonical gate at canonical coordinate `target`. The circuits supplied by
`basis_embodiments` are used to instantiate the individual XX actions.
NOTE: The elements of `strength_sequence` are expected to be normalized so that np.pi/2
corresponds to RZX(np.pi/2) = CX; `target` is taken to be a positive canonical coordinate;
and `basis_embodiments` maps `strength_sequence` elements to circuits which instantiate
these gates.
"""
# empty decompositions are easy!
if len(strength_sequence) == 0:
return QuantumCircuit(2)

# assemble the prefix / affix circuits
prefix_circuit, affix_circuit = QuantumCircuit(2), QuantumCircuit(2)
while len(strength_sequence) > 1:
source = decomposition_hop(target, strength_sequence)
strength = strength_sequence[-1]

preceding_prefix_circuit, preceding_affix_circuit = itemgetter(
"prefix_circuit", "affix_circuit"
)(xx_circuit_step(source, strength / 2, target, basis_embodiments[strength]))

prefix_circuit.compose(preceding_prefix_circuit, inplace=True)
affix_circuit.compose(preceding_affix_circuit, inplace=True, front=True)

target, strength_sequence = source, strength_sequence[:-1]

circuit = prefix_circuit

# lastly, deal with the "leading" gate.
if target[0] <= np.pi / 4:
circuit.compose(basis_embodiments[strength_sequence[0]], inplace=True)
else:
_, source_reflection, reflection_phase_shift = apply_reflection("reflect XX, YY", [0, 0, 0])
_, source_shift, shift_phase_shift = apply_shift("X shift", [0, 0, 0])

circuit.compose(source_reflection, inplace=True)
circuit.compose(basis_embodiments[strength_sequence[0]], inplace=True)
circuit.compose(source_reflection.inverse(), inplace=True)
circuit.compose(source_shift, inplace=True)
circuit.global_phase += -np.log(shift_phase_shift).imag
circuit.global_phase += -np.log(reflection_phase_shift).imag

circuit.compose(affix_circuit, inplace=True)

return circuit
278 changes: 278 additions & 0 deletions qiskit/quantum_info/synthesis/xx_decompose/decomposer.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,278 @@
# This code is part of Qiskit.
#
# (C) Copyright IBM 2021
#
# This code is licensed under the Apache License, Version 2.0. You may
# obtain a copy of this license in the LICENSE.txt file in the root directory
# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
#
# Any modifications or derivative works of this code must retain this
# copyright notice, and modified files need to carry a notice indicating
# that they have been altered from the originals.

"""
Driver for a synthesis routine which emits optimal XX-based circuits.
"""

import heapq
import math
from operator import itemgetter
from typing import Callable, Optional

import numpy as np

from qiskit.circuit.quantumcircuit import QuantumCircuit
from qiskit.circuit.library.standard_gates import RXXGate, RZXGate
from qiskit.exceptions import QiskitError
from qiskit.extensions import UnitaryGate
from qiskit.quantum_info.operators import Operator, average_gate_fidelity
from qiskit.quantum_info.synthesis.one_qubit_decompose import ONE_QUBIT_EULER_BASIS_GATES
from qiskit.quantum_info.synthesis.two_qubit_decompose import TwoQubitWeylDecomposition
from qiskit.transpiler.passes.optimization.optimize_1q_decomposition import (
Optimize1qGatesDecomposition,
)

from .circuits import apply_reflection, apply_shift, canonical_xx_circuit
from .utilities import EPSILON
from .polytopes import XXPolytope


def _average_infidelity(p, q):
"""
Computes the infidelity distance between two points p, q expressed in positive canonical
coordinates.
"""

a0, b0, c0 = p
a1, b1, c1 = q

return 1 - 1 / 20 * (
4
+ 16
* (
math.cos(a0 - a1) ** 2 * math.cos(b0 - b1) ** 2 * math.cos(c0 - c1) ** 2
+ math.sin(a0 - a1) ** 2 * math.sin(b0 - b1) ** 2 * math.sin(c0 - c1) ** 2
)
)


class XXDecomposer:
"""
A class for optimal decomposition of 2-qubit unitaries into 2-qubit basis gates of XX type
(i.e., each locally equivalent to CAN(alpha, 0, 0) for a possibly varying alpha).
Args:
euler_basis: Basis string provided to OneQubitEulerDecomposer for 1Q synthesis.
Defaults to "U".
embodiments: A dictionary mapping interaction strengths alpha to native circuits which
embody the gate CAN(alpha, 0, 0). Strengths are taken so that pi/2 represents the class
of a full CX.
backup_optimizer: If supplied, defers synthesis to this callable when XXDecomposer
has no efficient decomposition of its own. Useful for special cases involving 2 or 3
applications of XX(pi/2), in which case standard synthesis methods provide lower
1Q gate count.
NOTE: If embodiments is not passed, or if an entry is missing, it will be populated as needed
using the method _default_embodiment.
"""

def __init__(
self,
euler_basis: str = "U",
embodiments: Optional[dict] = None,
backup_optimizer: Optional[Callable] = None,
):
self._decomposer1q = Optimize1qGatesDecomposition(ONE_QUBIT_EULER_BASIS_GATES[euler_basis])
self.gate = RZXGate(np.pi / 2)
self.embodiments = embodiments if embodiments is not None else {}
self.backup_optimizer = backup_optimizer

self._check_embodiments()

@staticmethod
def _default_embodiment(strength):
"""
If the user does not provide a custom implementation of XX(strength), then this routine
defines a default implementation using RZX.
"""
xx_circuit = QuantumCircuit(2)

# NOTE: One could branch here on `strength == np.pi / 2` and decide to use a CX-based
# circuit in this one case where it's available.
xx_circuit.h(0)
xx_circuit.rzx(strength, 0, 1)
xx_circuit.h(0)

return xx_circuit

def _check_embodiments(self):
"""
Checks that `self.embodiments` is populated with legal circuit embodiments: the key-value
pair (angle, circuit) satisfies Operator(circuit) approx RXX(angle).to_matrix().
"""

for angle, embodiment in self.embodiments:
ecpeterson marked this conversation as resolved.
Show resolved Hide resolved
actual = Operator(RXXGate(angle))
purported = Operator(embodiment)
if average_gate_fidelity(actual, purported) < 1 - EPSILON:
raise QiskitError(
f"RXX embodiment provided for angle {angle} disagrees with RXXGate({angle})"
)

@staticmethod
def _best_decomposition(canonical_coordinate, available_strengths):
"""
Finds the cheapest sequence of `available_strengths` which supports the best approximation
to `canonical_coordinate`. Returns a dictionary with keys "cost", "point", and "operations".
NOTE: `canonical_coordinate` is a positive canonical coordinate. `strengths` is a dictionary
mapping the available strengths to their (infidelity) costs, with the strengths
themselves normalized so that pi/2 represents CX = RZX(pi/2).
"""
best_point, best_cost, best_sequence = [0, 0, 0], 1.0, []
priority_queue = []
heapq.heappush(priority_queue, (0, []))
canonical_coordinate = np.array(canonical_coordinate)

while True:
sequence_cost, sequence = heapq.heappop(priority_queue)

strength_polytope = XXPolytope.from_strengths(*[x / 2 for x in sequence])
candidate_point = strength_polytope.nearest(canonical_coordinate)
candidate_cost = sequence_cost + _average_infidelity(
canonical_coordinate, candidate_point
)

if candidate_cost < best_cost:
best_point, best_cost, best_sequence = candidate_point, candidate_cost, sequence

if strength_polytope.member(canonical_coordinate):
break

for strength, extra_cost in available_strengths.items():
if len(sequence) == 0 or strength <= sequence[-1]:
heapq.heappush(
priority_queue, (sequence_cost + extra_cost, sequence + [strength])
)

return {"point": best_point, "cost": best_cost, "sequence": best_sequence}

def num_basis_gates(self, unitary):
"""
Counts the number of gates that would be emitted during re-synthesis.
NOTE: Used by ConsolidateBlocks.
"""
strengths = self._strength_to_infidelity(1.0)

# get the associated _positive_ canonical coordinate
weyl_decomposition = TwoQubitWeylDecomposition(unitary)
target = [getattr(weyl_decomposition, x) for x in ("a", "b", "c")]
if target[-1] < -EPSILON:
target = [np.pi / 2 - target[0], target[1], -target[2]]

best_sequence = self._best_decomposition(target, strengths)["sequence"]
return len(best_sequence)

@staticmethod
def _strength_to_infidelity(basis_fidelity, approximate=False):
"""
Converts a dictionary mapping ZX strengths to fidelities to a dictionary mapping ZX
strengths to infidelities. Also supports one of the other formats Qiskit uses: if only a
lone float is supplied, it extends it from CX over CX/2 and CX/3 by linear decay.
"""

if isinstance(basis_fidelity, float):
if not approximate:
slope, offset = 1e-10, 1e-12
else:
slope, offset = (1 - basis_fidelity) / 2, (1 - basis_fidelity) / 2
return {
strength: slope * strength / (np.pi / 2) + offset
for strength in [np.pi / 2, np.pi / 4, np.pi / 6]
}
elif isinstance(basis_fidelity, dict):
return {
strength: (1 - fidelity if approximate else 1e-12 + 1e-10 * strength / (np.pi / 2))
for (strength, fidelity) in basis_fidelity.items()
}

raise TypeError("Unknown basis_fidelity payload.")

def __call__(self, unitary, basis_fidelity=1.0, approximate=True):
"""
Fashions a circuit which (perhaps `approximate`ly) models the special unitary operation
`unitary`, using the circuit templates supplied at initialization as `embodiments`. The
routine uses `basis_fidelity` to select the optimal circuit template, including when
performing exact synthesis; the contents of `basis_fidelity` is a dictionary mapping
interaction strengths (scaled so that CX = RZX(pi/2) corresponds to pi/2) to circuit
fidelities.
Args:
unitary (Operator or ndarray): 4x4 unitary to synthesize.
basis_fidelity (dict or float): Fidelity of basis gates. Can be either (1) a dictionary
mapping XX angle values to fidelity at that angle; or (2) a single float f,
interpreted as {pi: f, pi/2: f/2, pi/3: f/3} .
approximate (bool): Approximates if basis fidelities are less than 1.0 .
Returns:
QuantumCircuit: Synthesized circuit.
"""
strength_to_infidelity = self._strength_to_infidelity(
basis_fidelity, approximate=approximate
)

# get the associated _positive_ canonical coordinate
weyl_decomposition = TwoQubitWeylDecomposition(unitary)
target = [getattr(weyl_decomposition, x) for x in ("a", "b", "c")]
if target[-1] < -EPSILON:
target = [np.pi / 2 - target[0], target[1], -target[2]]

# scan for the best point
best_point, best_sequence = itemgetter("point", "sequence")(
self._best_decomposition(target, strength_to_infidelity)
)
# build the circuit building this canonical gate
embodiments = {
k: self.embodiments.get(k, self._default_embodiment(k))
for k, v in strength_to_infidelity.items()
}
ecpeterson marked this conversation as resolved.
Show resolved Hide resolved
circuit = canonical_xx_circuit(best_point, best_sequence, embodiments)

if best_sequence == [np.pi / 2, np.pi / 2, np.pi / 2] and self.backup_optimizer is not None:
return self.backup_optimizer(unitary, basis_fidelity=basis_fidelity)

# change to positive canonical coordinates
if weyl_decomposition.c >= -EPSILON:
# if they're the same...
corrected_circuit = QuantumCircuit(2)
corrected_circuit.rz(np.pi, [0])
corrected_circuit.compose(circuit, [0, 1], inplace=True)
corrected_circuit.rz(-np.pi, [0])
circuit = corrected_circuit
else:
# else they're in the "positive" scissors part...
corrected_circuit = QuantumCircuit(2)
_, source_reflection, _ = apply_reflection("reflect XX, ZZ", [0, 0, 0])
_, source_shift, _ = apply_shift("X shift", [0, 0, 0])

corrected_circuit.compose(source_reflection.inverse(), inplace=True)
corrected_circuit.rz(np.pi, [0])
corrected_circuit.compose(circuit, [0, 1], inplace=True)
corrected_circuit.rz(-np.pi, [0])
corrected_circuit.compose(source_shift.inverse(), inplace=True)
corrected_circuit.compose(source_reflection, inplace=True)
corrected_circuit.global_phase += np.pi / 2

circuit = corrected_circuit

circ = QuantumCircuit(2, global_phase=weyl_decomposition.global_phase)

circ.append(UnitaryGate(weyl_decomposition.K2r), [0])
circ.append(UnitaryGate(weyl_decomposition.K2l), [1])
circ.compose(circuit, [0, 1], inplace=True)
circ.append(UnitaryGate(weyl_decomposition.K1r), [0])
circ.append(UnitaryGate(weyl_decomposition.K1l), [1])

circ = self._decomposer1q(circ)

return circ
395 changes: 395 additions & 0 deletions qiskit/quantum_info/synthesis/xx_decompose/paths.py

Large diffs are not rendered by default.

266 changes: 266 additions & 0 deletions qiskit/quantum_info/synthesis/xx_decompose/polytopes.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,266 @@
# This code is part of Qiskit.
#
# (C) Copyright IBM 2021
#
# This code is licensed under the Apache License, Version 2.0. You may
# obtain a copy of this license in the LICENSE.txt file in the root directory
# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
#
# Any modifications or derivative works of this code must retain this
# copyright notice, and modified files need to carry a notice indicating
# that they have been altered from the originals.

"""
Defines bare dataclasses which house polytope information, as well as a specialized data structure
which describes those two-qubit programs accessible to a given sequence of XX-type interactions.
"""

from copy import copy
from dataclasses import dataclass, field
ecpeterson marked this conversation as resolved.
Show resolved Hide resolved
from itertools import combinations
import random
from typing import List

import numpy as np

from qiskit.exceptions import QiskitError

from .utilities import EPSILON


@dataclass
class ConvexPolytopeData:
"""
The raw data underlying a ConvexPolytope. Describes a single convex
polytope, specified by families of `inequalities` and `equalities`, each
entry of which respectively corresponds to
inequalities[j][0] + sum_i inequalities[j][i] * xi >= 0
and
equalities[j][0] + sum_i equalities[j][i] * xi == 0.
"""

inequalities: List[List[int]]
equalities: List[List[int]] = field(default_factory=list)
name: str = ""


@dataclass
class PolytopeData:
"""
The raw data of a union of convex polytopes.
"""

convex_subpolytopes: List[ConvexPolytopeData]


def polytope_has_element(polytope, point):
"""
Tests whether `polytope` contains `point.
"""
return all(
-EPSILON <= inequality[0] + sum(x * y for x, y in zip(point, inequality[1:]))
for inequality in polytope.inequalities
) and all(
abs(equality[0] + sum(x * y for x, y in zip(point, equality[1:]))) <= EPSILON
for equality in polytope.equalities
)


def manual_get_vertex(polytope, seed=42):
"""
Returns a single random vertex from `polytope`.
"""
random.seed(seed)
ecpeterson marked this conversation as resolved.
Show resolved Hide resolved

if isinstance(polytope, PolytopeData):
paragraphs = copy(polytope.convex_subpolytopes)
elif isinstance(polytope, ConvexPolytopeData):
paragraphs = [polytope]
else:
raise TypeError(f"{type(polytope)} is not polytope-like.")

random.shuffle(paragraphs)
for convex_subpolytope in paragraphs:
sentences = convex_subpolytope.inequalities + convex_subpolytope.equalities
if len(sentences) == 0:
continue
dimension = len(sentences[0]) - 1
random.shuffle(sentences)
for inequalities in combinations(sentences, dimension):
matrix = np.array([x[1:] for x in inequalities])
b = np.array([x[0] for x in inequalities])
try:
vertex = np.linalg.inv(-matrix) @ b
if polytope_has_element(convex_subpolytope, vertex):
return vertex
except np.linalg.LinAlgError:
pass

raise QiskitError(f"Polytope has no feasible solutions:\n{polytope}")


@dataclass
class XXPolytope:
"""
Describes those two-qubit programs accessible to a given sequence of XX-type interactions.
NOTE: Strengths are normalized so that CX corresponds to pi / 4, which differs from Qiskit's
conventions around RZX elsewhere.
"""

# NOTE: This is _not_ a subclass of PolytopeData, because we're never going to call slow,
# generic PolytopeData functions on it.

total_strength: float = 0.0
max_strength: float = 0.0
place_strength: float = 0.0

@classmethod
def from_strengths(cls, *strengths):
"""
Constructs an XXPolytope from a sequence of strengths.
"""
total_strength, max_strength, place_strength = 0, 0, 0
for strength in strengths:
total_strength += strength
if strength >= max_strength:
max_strength, place_strength = strength, max_strength
elif strength >= place_strength:
place_strength = strength

return XXPolytope(
total_strength=total_strength, max_strength=max_strength, place_strength=place_strength
)

def add_strength(self, new_strength: float = 0.0):
"""
Returns a new XXPolytope with one new XX interaction appended.
"""
return XXPolytope(
total_strength=self.total_strength + new_strength,
max_strength=max(self.max_strength, new_strength),
place_strength=(
self.max_strength
if new_strength > self.max_strength
else new_strength
if new_strength > self.place_strength
else self.place_strength
),
)

@property
def _offsets(self):
"""
Returns b with A*x + b ≥ 0 iff x belongs to the XXPolytope.
"""
return np.array(
[
0,
0,
0,
np.pi / 2,
self.total_strength,
self.total_strength - 2 * self.max_strength,
self.total_strength - self.max_strength - self.place_strength,
]
)

def member(self, point):
"""
Returns True when `point` is a member of `self`.
"""
global A # pylint:disable=global-statement

reflected_point = point.copy().reshape(-1, 3)
rows = reflected_point[:, 0] >= np.pi / 4 + EPSILON
reflected_point[rows, 0] = np.pi / 2 - reflected_point[rows, 0]
reflected_point = reflected_point.reshape(point.shape)

return np.all(
self._offsets + np.einsum("ij,...j->...i", A, reflected_point) >= -EPSILON, axis=-1
)

def nearest(self, point):
"""
Finds the nearest point (in Euclidean or infidelity distance) to `self`.
"""
# pylint:disable=invalid-name

# NOTE: A CAS says that there are no degenerate double intersections, and the only
# degenerate triple intersections are
#
# (1, -1, 0), (0, 0, 1), (-1, 1, 1) and (1, 1, 0), (0, 0, 1), (1, 1, 1).
#
# Skipping this pair won't save much work, so we don't bother.

global A1, A1inv, A2, A2inv, A3, A3inv # pylint:disable=global-statement
# These global variables contain projection matrices, computed once-and-for-all, which
# produce the Euclidean-nearest projection.

if isinstance(point, np.ndarray) and len(point.shape) == 1:
y0 = point.copy()
elif isinstance(point, list):
y0 = np.array(point)
else:
raise TypeError(f"Can't handle type of point: {point} ({type(point)})")

reflected_p = y0[0] > np.pi / 4 + EPSILON
if reflected_p:
y0[0] = np.pi / 2 - y0[0]

# short circuit in codimension 0
if self.member(y0):
if reflected_p:
y0[0] = np.pi / 2 - y0[0]
return y0

# codimension 1
b1 = self._offsets.reshape(7, 1)
A1y0 = np.einsum("ijk,k->ij", A1, y0)
nearest1 = np.einsum("ijk,ik->ij", A1inv, b1 + A1y0) - y0

# codimension 2
b2 = np.array([*combinations(self._offsets, 2)])
A2y0 = np.einsum("ijk,k->ij", A2, y0)
nearest2 = np.einsum("ijk,ik->ij", A2inv, b2 + A2y0) - y0

# codimension 3
b3 = np.array([*combinations(self._offsets, 3)])
nearest3 = np.einsum("ijk,ik->ij", A3inv, b3)

# pick the nearest
nearest = -np.concatenate([nearest1, nearest2, nearest3])
nearest = nearest[self.member(nearest)]
smallest_index = np.argmin(np.linalg.norm(nearest - y0, axis=1))

if reflected_p:
nearest[smallest_index][0] = np.pi / 2 - nearest[smallest_index][0]
return nearest[smallest_index]


A = np.array(
[
[1, -1, 0], # a ≥ b
[0, 1, -1], # b ≥ c
[0, 0, 1], # c ≥ 0
[-1, -1, 0], # pi/2 ≥ a + b
[-1, -1, -1], # strength
[1, -1, -1], # slant
[0, 0, -1], # frustrum
]
)
A1 = A.reshape(-1, 1, 3) # pylint:disable=too-many-function-args
A1inv = np.linalg.pinv(A1)
A2 = np.array([np.array([x, y], dtype=float) for (x, y) in combinations(A, 2)])
A2inv = np.linalg.pinv(A2)
A3 = np.array([np.array([x, y, z], dtype=float) for (x, y, z) in combinations(A, 3)])
A3inv = np.linalg.pinv(A3)
"""
These globals house matrices, computed once-and-for-all, which project onto the Euclidean-nearest
planes spanned by the planes/points/lines defined by the faces of any XX polytope.
See `XXPolytope.nearest`.
"""
39 changes: 39 additions & 0 deletions qiskit/quantum_info/synthesis/xx_decompose/utilities.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,39 @@
# This code is part of Qiskit.
#
# (C) Copyright IBM 2021
#
# This code is licensed under the Apache License, Version 2.0. You may
# obtain a copy of this license in the LICENSE.txt file in the root directory
# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
#
# Any modifications or derivative works of this code must retain this
# copyright notice, and modified files need to carry a notice indicating
# that they have been altered from the originals.

"""
Depository for generic utility snippets.
"""

import warnings

import numpy as np


EPSILON = 1e-6 # Fraction(1, 1_000_000)


# TODO: THIS IS A STOPGAP!!!
def safe_arccos(numerator, denominator):
"""
Computes arccos(n/d) with different (better?) numerical stability.
"""
threshold = 0.005

if abs(numerator) > abs(denominator) and abs(numerator - denominator) < threshold:
return 0.0
elif abs(numerator) > abs(denominator) and abs(numerator + denominator) < threshold:
return np.pi
else:
with warnings.catch_warnings():
warnings.filterwarnings("ignore", category=RuntimeWarning)
return np.arccos(numerator / denominator)
132 changes: 132 additions & 0 deletions qiskit/quantum_info/synthesis/xx_decompose/weyl.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,132 @@
# This code is part of Qiskit.
#
# (C) Copyright IBM 2021
#
# This code is licensed under the Apache License, Version 2.0. You may
# obtain a copy of this license in the LICENSE.txt file in the root directory
# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
#
# Any modifications or derivative works of this code must retain this
# copyright notice, and modified files need to carry a notice indicating
# that they have been altered from the originals.

"""
Simple circuit constructors for Weyl reflections.
"""

import numpy as np

from qiskit.circuit.quantumcircuit import QuantumCircuit
from qiskit.circuit.library.standard_gates import RXGate, RYGate, RZGate


reflection_options = {
"no reflection": ([1, 1, 1], 1, []),
"reflect XX, YY": ([-1, -1, 1], 1, [RZGate]),
"reflect XX, ZZ": ([-1, 1, -1], 1, [RYGate]),
"reflect YY, ZZ": ([1, -1, -1], 1, [RXGate]),
}
"""
A table of available reflection transformations on canonical coordinates.
Entries take the form
readable_name: (reflection scalars, global phase, [gate constructors]),
where reflection scalars (a, b, c) model the map (x, y, z) |-> (ax, by, cz),
global phase is a complex unit, and gate constructors are applied in sequence
and by conjugation to the first qubit and are passed pi as a parameter.
"""

shift_options = {
"no shift": ([0, 0, 0], 1, []),
"Z shift": ([0, 0, 1], 1j, [RZGate]),
"Y shift": ([0, 1, 0], -1j, [RYGate]),
"Y,Z shift": ([0, 1, 1], 1, [RYGate, RZGate]),
"X shift": ([1, 0, 0], -1j, [RXGate]),
"X,Z shift": ([1, 0, 1], 1, [RXGate, RZGate]),
"X,Y shift": ([1, 1, 0], -1, [RXGate, RYGate]),
"X,Y,Z shift": ([1, 1, 1], -1j, [RXGate, RYGate, RZGate]),
}
"""
A table of available shift transformations on canonical coordinates. Entries
take the form
readable name: (shift scalars, global phase, [gate constructors]),
where shift scalars model the map
(x, y, z) |-> (x + a pi / 2, y + b pi / 2, z + c pi / 2) ,
global phase is a complex unit, and gate constructors are applied to the first
and second qubits and are passed pi as a parameter.
"""


def apply_reflection(reflection_name, coordinate):
"""
Given a reflection type and a canonical coordinate, applies the reflection
and describes a circuit which enacts the reflection + a global phase shift.
"""
reflection_scalars, reflection_phase_shift, source_reflection_gates = reflection_options[
reflection_name
]
reflected_coord = [x * y for x, y in zip(reflection_scalars, coordinate)]
source_reflection = QuantumCircuit(2)
for gate in source_reflection_gates:
source_reflection.append(gate(np.pi), [0])

return reflected_coord, source_reflection, reflection_phase_shift


def apply_shift(shift_name, coordinate):
"""
Given a shift type and a canonical coordinate, applies the shift and
describes a circuit which enacts the shift + a global phase shift.
"""
shift_scalars, shift_phase_shift, source_shift_gates = shift_options[shift_name]
shifted_coord = [np.pi / 2 * x + y for x, y in zip(shift_scalars, coordinate)]

source_shift = QuantumCircuit(2)
for gate in source_shift_gates:
source_shift.append(gate(np.pi), [0])
source_shift.append(gate(np.pi), [1])

return shifted_coord, source_shift, shift_phase_shift


def canonical_rotation_circuit(first_index, second_index):
"""
Given a pair of distinct indices 0 ≤ (first_index, second_index) ≤ 2,
produces a two-qubit circuit which rotates a canonical gate
a0 XX + a1 YY + a2 ZZ
into
a[first] XX + a[second] YY + a[other] ZZ .
"""
conj = QuantumCircuit(2)

if (0, 1) == (first_index, second_index):
pass # no need to do anything
elif (0, 2) == (first_index, second_index):
conj.rx(-np.pi / 2, [0])
conj.rx(np.pi / 2, [1])
elif (1, 0) == (first_index, second_index):
conj.rz(-np.pi / 2, [0])
conj.rz(-np.pi / 2, [1])
elif (1, 2) == (first_index, second_index):
conj.rz(np.pi / 2, [0])
conj.rz(np.pi / 2, [1])
conj.ry(np.pi / 2, [0])
conj.ry(-np.pi / 2, [1])
elif (2, 0) == (first_index, second_index):
conj.rz(np.pi / 2, [0])
conj.rz(np.pi / 2, [1])
conj.rx(np.pi / 2, [0])
conj.rx(-np.pi / 2, [1])
elif (2, 1) == (first_index, second_index):
conj.ry(np.pi / 2, [0])
conj.ry(-np.pi / 2, [1])

return conj
5 changes: 4 additions & 1 deletion qiskit/transpiler/layout.py
Original file line number Diff line number Diff line change
@@ -308,7 +308,10 @@ def from_intlist(int_list, *qregs):
num_qubits = sum(reg.size for reg in qregs)
# Check if list is too short to cover all qubits
if len(int_list) != num_qubits:
raise LayoutError("Integer list length must equal number of qubits in circuit.")
raise LayoutError(
f"Integer list length ({len(int_list)}) must equal number of qubits "
f"in circuit ({num_qubits}): {int_list}."
)
out = Layout()
main_idx = 0
for qreg in qregs:
Original file line number Diff line number Diff line change
@@ -56,11 +56,7 @@ def __init__(self, kak_basis_gate=None, force_consolidate=False, basis_gates=Non
if kak_basis_gate is not None:
self.decomposer = TwoQubitBasisDecomposer(kak_basis_gate)
elif basis_gates is not None:
kak_basis_gate = unitary_synthesis._choose_kak_gate(basis_gates)
if kak_basis_gate is not None:
self.decomposer = TwoQubitBasisDecomposer(kak_basis_gate)
else:
self.decomposer = None
self.decomposer = unitary_synthesis._basis_gates_to_decomposer_2q(basis_gates)
else:
self.decomposer = TwoQubitBasisDecomposer(CXGate())

39 changes: 31 additions & 8 deletions qiskit/transpiler/passes/synthesis/unitary_synthesis.py
Original file line number Diff line number Diff line change
@@ -23,8 +23,16 @@
from qiskit.dagcircuit.dagcircuit import DAGCircuit
from qiskit.extensions.quantum_initializer import isometry
from qiskit.quantum_info.synthesis import one_qubit_decompose
from qiskit.quantum_info.synthesis.xx_decompose import XXDecomposer
from qiskit.quantum_info.synthesis.two_qubit_decompose import TwoQubitBasisDecomposer
from qiskit.circuit.library.standard_gates import iSwapGate, CXGate, CZGate, RXXGate, ECRGate
from qiskit.circuit.library.standard_gates import (
iSwapGate,
CXGate,
CZGate,
RXXGate,
RZXGate,
ECRGate,
)
from qiskit.transpiler.passes.synthesis import plugin
from qiskit.providers.models import BackendProperties

@@ -38,6 +46,7 @@ def _choose_kak_gate(basis_gates):
"iswap": iSwapGate(),
"rxx": RXXGate(pi / 2),
"ecr": ECRGate(),
"rzx": RZXGate(pi / 4), # typically pi/6 is also available
}

kak_gate = None
@@ -78,6 +87,23 @@ def _choose_bases(basis_gates, basis_dict=None):
return out_basis


def _basis_gates_to_decomposer_2q(basis_gates, pulse_optimize=None):
kak_gate = _choose_kak_gate(basis_gates)
euler_basis = _choose_euler_basis(basis_gates)

if isinstance(kak_gate, RZXGate):
backup_optimizer = TwoQubitBasisDecomposer(
CXGate(), euler_basis=euler_basis, pulse_optimize=pulse_optimize
)
return XXDecomposer(euler_basis=euler_basis, backup_optimizer=backup_optimizer)
elif kak_gate is not None:
return TwoQubitBasisDecomposer(
kak_gate, euler_basis=euler_basis, pulse_optimize=pulse_optimize
)
else:
return None


class UnitarySynthesis(TransformationPass):
"""Synthesize gates according to their basis gates."""

@@ -332,15 +358,12 @@ def run(self, unitary, **options):
qubits = options["coupling_map"][1]

euler_basis = _choose_euler_basis(basis_gates)
kak_gate = _choose_kak_gate(basis_gates)

decomposer1q, decomposer2q = None, None
if euler_basis is not None:
decomposer1q = one_qubit_decompose.OneQubitEulerDecomposer(euler_basis)
if kak_gate is not None:
decomposer2q = TwoQubitBasisDecomposer(
kak_gate, euler_basis=euler_basis, pulse_optimize=pulse_optimize
)
else:
decomposer1q = None

decomposer2q = _basis_gates_to_decomposer_2q(basis_gates, pulse_optimize=pulse_optimize)

synth_dag = None
wires = None
26 changes: 23 additions & 3 deletions qiskit/transpiler/passes/utils/gate_direction.py
Original file line number Diff line number Diff line change
@@ -20,7 +20,7 @@

from qiskit.circuit import QuantumRegister
from qiskit.dagcircuit import DAGCircuit
from qiskit.circuit.library.standard_gates import RYGate, HGate, CXGate, ECRGate
from qiskit.circuit.library.standard_gates import RYGate, HGate, CXGate, ECRGate, RZXGate


class GateDirection(TransformationPass):
@@ -39,6 +39,12 @@ class GateDirection(TransformationPass):
│ ECR │ = └┬──────────┤│ ECR │├───┤
q_1: ┤1 ├ q_1: ─┤ RY(pi/2) ├┤0 ├┤ H ├
└──────┘ └──────────┘└──────┘└───┘
┌──────┐ ┌───┐┌──────┐┌───┐
q_0: ┤0 ├ q_0: ┤ H ├┤1 ├┤ H ├
│ RZX │ = ├───┤│ RZX │├───┤
q_1: ┤1 ├ q_1: ┤ H ├┤0 ├┤ H ├
└──────┘ └───┘└──────┘└───┘
"""

def __init__(self, coupling_map):
@@ -69,6 +75,18 @@ def __init__(self, coupling_map):
self._ecr_dag.apply_operation_back(HGate(), [qr[0]], [])
self._ecr_dag.apply_operation_back(HGate(), [qr[1]], [])

@staticmethod
def _rzx_dag(parameter):
_rzx_dag = DAGCircuit()
qr = QuantumRegister(2)
_rzx_dag.add_qreg(qr)
_rzx_dag.apply_operation_back(HGate(), [qr[0]], [])
_rzx_dag.apply_operation_back(HGate(), [qr[1]], [])
_rzx_dag.apply_operation_back(RZXGate(parameter), [qr[1], qr[0]], [])
_rzx_dag.apply_operation_back(HGate(), [qr[0]], [])
_rzx_dag.apply_operation_back(HGate(), [qr[1]], [])
return _rzx_dag

def run(self, dag):
"""Run the GateDirection pass on `dag`.
@@ -113,10 +131,12 @@ def run(self, dag):
dag.substitute_node_with_dag(node, self._cx_dag)
elif node.name == "ecr":
dag.substitute_node_with_dag(node, self._ecr_dag)
elif node.name == "rzx":
dag.substitute_node_with_dag(node, self._rzx_dag(*node.op.params))
else:
raise TranspilerError(
"Flipping of gate direction is only supported "
"for CX and ECR at this time."
f"Flipping of gate direction is only supported "
f"for CX, ECR, and RZX at this time, not {node.name}."
)

return dag
Original file line number Diff line number Diff line change
@@ -0,0 +1,8 @@
---
features:
- |
This release introduces a decomposition method for two-qubit gates which
targets user-defined sets of RZX gates. Transpiler users can enable
decomposition for {RZX(pi/2), RZX(pi/4), and RZX(pi/6)} specifically by including
`'rzx'` in their `basis_gates` list. Quantum information package users can
find the method itself under the `XXDecomposer` class.
13 changes: 13 additions & 0 deletions test/python/quantum_info/xx_decompose/__init__.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,13 @@
# This code is part of Qiskit.
#
# (C) Copyright IBM 2021
#
# This code is licensed under the Apache License, Version 2.0. You may
# obtain a copy of this license in the LICENSE.txt file in the root directory
# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
#
# Any modifications or derivative works of this code must retain this
# copyright notice, and modified files need to carry a notice indicating
# that they have been altered from the originals.

"""Qiskit integration tests for XX-based synthesis."""
142 changes: 142 additions & 0 deletions test/python/quantum_info/xx_decompose/test_circuits.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,142 @@
# This code is part of Qiskit.
#
# (C) Copyright IBM 2021
#
# This code is licensed under the Apache License, Version 2.0. You may
# obtain a copy of this license in the LICENSE.txt file in the root directory
# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
#
# Any modifications or derivative works of this code must retain this
# copyright notice, and modified files need to carry a notice indicating
# that they have been altered from the originals.

"""
Tests for qiskit-terra/qiskit/quantum_info/synthesis/xx_decompose/circuits.py .
"""

from operator import itemgetter
import random
import unittest

import ddt
import numpy as np

from qiskit.circuit import QuantumCircuit
from qiskit.circuit.library import RZGate
import qiskit.quantum_info.operators
from qiskit.quantum_info.synthesis.weyl import weyl_coordinates
from qiskit.quantum_info.synthesis.xx_decompose.circuits import (
decompose_xxyy_into_xxyy_xx,
xx_circuit_step,
)

from .utilities import canonical_matrix

EPSILON = 0.001


@ddt.ddt
class TestMonodromyCircuits(unittest.TestCase):
"""Check circuit synthesis step routines."""

def __init__(self, *args, seed=42, **kwargs):
super().__init__(*args, **kwargs)
random.seed(seed)
np.random.seed(seed)

def _generate_xxyy_test_case(self):
"""
Generates a random pentuple of values (a_source, b_source), beta, (a_target, b_target) s.t.
CAN(a_source, b_source) * exp(-i(s ZI + t IZ)) * CAN(beta)
=
exp(-i(u ZI + v IZ)) * CAN(a_target, b_target) * exp(-i(x ZI + y IZ))
admits a solution in (s, t, u, v, x, y).
Returns (source_coordinate, interaction, target_coordinate).
"""
source_coordinate = [random.random(), random.random(), 0.0]
source_coordinate = [
source_coordinate[0] * np.pi / 8,
source_coordinate[1] * source_coordinate[0] * np.pi / 8,
0.0,
]
interaction = [random.random() * np.pi / 8]
z_angles = [random.random() * np.pi / 8, random.random() * np.pi / 8]
prod = (
canonical_matrix(*source_coordinate)
@ np.kron(RZGate(2 * z_angles[0]).to_matrix(), RZGate(2 * z_angles[1]).to_matrix())
@ canonical_matrix(interaction[0], 0.0, 0.0)
)
target_coordinate = weyl_coordinates(prod)

self.assertAlmostEqual(target_coordinate[-1], 0.0, delta=EPSILON)

return source_coordinate, interaction, target_coordinate

# pylint:disable=invalid-name
def test_decompose_xxyy(self):
"""
Test that decompose_xxyy_into_xxyy_xx correctly recovers decompositions.
"""

for _ in range(100):
source_coordinate, interaction, target_coordinate = self._generate_xxyy_test_case()

r, s, u, v, x, y = decompose_xxyy_into_xxyy_xx(
target_coordinate[0],
target_coordinate[1],
source_coordinate[0],
source_coordinate[1],
interaction[0],
)

prod = (
np.kron(RZGate(2 * r).to_matrix(), RZGate(2 * s).to_matrix())
@ canonical_matrix(*source_coordinate)
@ np.kron(RZGate(2 * u).to_matrix(), RZGate(2 * v).to_matrix())
@ canonical_matrix(interaction[0], 0.0, 0.0)
@ np.kron(RZGate(2 * x).to_matrix(), RZGate(2 * y).to_matrix())
)
expected = canonical_matrix(*target_coordinate)
self.assertTrue(np.all(np.abs(prod - expected) < EPSILON))

def test_xx_circuit_step(self):
"""
Test that `xx_circuit_step` correctly generates prefix/affix circuits relating source
canonical coordinates to target canonical coordinates along prescribed interactions, all
randomly selected.
"""

for _ in range(100):
source_coordinate, interaction, target_coordinate = self._generate_xxyy_test_case()

source_embodiment = qiskit.QuantumCircuit(2)
source_embodiment.append(
qiskit.extensions.UnitaryGate(canonical_matrix(*source_coordinate)), [0, 1]
)
interaction_embodiment = qiskit.QuantumCircuit(2)
interaction_embodiment.append(
qiskit.extensions.UnitaryGate(canonical_matrix(*interaction)), [0, 1]
)

prefix_circuit, affix_circuit = itemgetter("prefix_circuit", "affix_circuit")(
xx_circuit_step(
source_coordinate, interaction[0], target_coordinate, interaction_embodiment
)
)

target_embodiment = QuantumCircuit(2)
target_embodiment.compose(prefix_circuit, inplace=True)
target_embodiment.compose(source_embodiment, inplace=True)
target_embodiment.compose(affix_circuit, inplace=True)
self.assertTrue(
np.all(
np.abs(
qiskit.quantum_info.operators.Operator(target_embodiment).data
- canonical_matrix(*target_coordinate)
)
< EPSILON
)
)
109 changes: 109 additions & 0 deletions test/python/quantum_info/xx_decompose/test_decomposer.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,109 @@
# This code is part of Qiskit.
#
# (C) Copyright IBM 2021
#
# This code is licensed under the Apache License, Version 2.0. You may
# obtain a copy of this license in the LICENSE.txt file in the root directory
# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
#
# Any modifications or derivative works of this code must retain this
# copyright notice, and modified files need to carry a notice indicating
# that they have been altered from the originals.

"""
Tests for qiskit-terra/qiskit/quantum_info/synthesis/xx_decompose/qiskit.py .
"""

import random
from statistics import mean
import unittest

import ddt
import numpy as np
from scipy.stats import unitary_group

from qiskit.quantum_info.operators import Operator
from qiskit.quantum_info.synthesis.xx_decompose.decomposer import (
XXDecomposer,
TwoQubitWeylDecomposition,
)

from .utilities import canonical_matrix

EPSILON = 0.001


@ddt.ddt
class TestMonodromyQISKit(unittest.TestCase):
"""Check QISKit routines."""

decomposer = XXDecomposer(euler_basis="PSX")

def __init__(self, *args, seed=42, **kwargs):
super().__init__(*args, **kwargs)
random.seed(seed)
np.random.seed(seed)

def test_random_compilation(self):
"""Test that compilation gives correct results."""
for _ in range(100):
unitary = unitary_group.rvs(4)
unitary /= np.linalg.det(unitary) ** (1 / 4)

# decompose into CX, CX/2, and CX/3
circuit = self.decomposer(unitary, approximate=False)
decomposed_unitary = Operator(circuit).data

self.assertTrue(np.all(unitary - decomposed_unitary < EPSILON))

def test_compilation_determinism(self):
"""Test that compilation is stable under multiple calls."""
for _ in range(10):
unitary = unitary_group.rvs(4)
unitary /= np.linalg.det(unitary) ** (1 / 4)

# decompose into CX, CX/2, and CX/3
circuit1 = self.decomposer(unitary, approximate=False)
circuit2 = self.decomposer(unitary, approximate=False)

self.assertEqual(circuit1, circuit2)

@ddt.data(np.pi / 3, np.pi / 5, np.pi / 2)
def test_default_embodiment(self, angle):
"""Test that _default_embodiment actually does yield XX gates."""
embodiment = self.decomposer._default_embodiment(angle)
embodiment_matrix = Operator(embodiment).data
self.assertTrue(np.all(canonical_matrix(angle, 0, 0) - embodiment_matrix < EPSILON))

def test_compilation_improvement(self):
"""Test that compilation to CX, CX/2, CX/3 improves over CX alone."""
slope, offset = (64 * 90) / 1000000, 909 / 1000000 + 1 / 1000
strength_table = self.decomposer._strength_to_infidelity(
basis_fidelity={
strength: 1 - (slope * strength / (np.pi / 2) + offset)
for strength in [np.pi / 2, np.pi / 4, np.pi / 6]
},
approximate=True,
)
limited_strength_table = {np.pi / 2: strength_table[np.pi / 2]}

clever_costs = []
naive_costs = []
for _ in range(200):
unitary = unitary_group.rvs(4)
unitary /= np.linalg.det(unitary) ** (1 / 4)

weyl_decomposition = TwoQubitWeylDecomposition(unitary)
target = [getattr(weyl_decomposition, x) for x in ("a", "b", "c")]
if target[-1] < -EPSILON:
target = [np.pi / 2 - target[0], target[1], -target[2]]

# decompose into CX, CX/2, and CX/3
clever_costs.append(self.decomposer._best_decomposition(target, strength_table)["cost"])
naive_costs.append(
self.decomposer._best_decomposition(target, limited_strength_table)["cost"]
)

# the following are taken from Fig 14 of the XX synthesis paper
self.assertAlmostEqual(mean(clever_costs), 1.445e-2, delta=5e-3)
self.assertAlmostEqual(mean(naive_costs), 2.058e-2, delta=5e-3)
61 changes: 61 additions & 0 deletions test/python/quantum_info/xx_decompose/test_polytopes.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,61 @@
# This code is part of Qiskit.
#
# (C) Copyright IBM 2021
#
# This code is licensed under the Apache License, Version 2.0. You may
# obtain a copy of this license in the LICENSE.txt file in the root directory
# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
#
# Any modifications or derivative works of this code must retain this
# copyright notice, and modified files need to carry a notice indicating
# that they have been altered from the originals.

"""
Tests for qiskit-terra/qiskit/quantum_info/synthesis/xx_decompose/xx_polytope.py .
"""

import random
import unittest

import ddt
import numpy as np
from numpy import pi

from qiskit.quantum_info.synthesis.xx_decompose.polytopes import XXPolytope

EPSILON = 0.001


@ddt.ddt
class TestMonodromyXXPolytope(unittest.TestCase):
"""Check specialized XX polytope routines."""

def __init__(self, *args, seed=42, **kwargs):
super().__init__(*args, **kwargs)
random.seed(seed)
np.random.seed(seed)

@ddt.data(
([pi / 4, pi / 4, pi / 4], [0.52359878, 0.39269908, 0.31415927]),
([pi / 6, pi / 8, pi / 10], [pi / 6, pi / 8, pi / 10]),
([pi / 4, pi / 4, 0], [0.615228561, 0.615228561, 0]),
([pi / 4, pi / 8, 0], [pi / 4, pi / 8, 0]),
)
@ddt.unpack
def test_nearest(self, offbody, expected):
"""Check that the nearest point calculator recovers some known cases."""
polytope = XXPolytope.from_strengths(pi / 6, pi / 8, pi / 10)
result = polytope.nearest(np.array(offbody))
print(offbody, result)
self.assertTrue(np.all(np.abs(np.array(expected) - result) < EPSILON))

def test_add_strengths(self):
"""
Tests that adding new strengths to an existing XXPolytope is equivalent to forming the
appropriate XXPolytope from scratch.
"""
for _ in range(100):
strengths = [random.random() for _ in range(4)]
small_polytope = XXPolytope.from_strengths(*strengths[:-1])
large_polytope = XXPolytope.from_strengths(*strengths)
self.assertEqual(small_polytope.add_strength(strengths[-1]), large_polytope)
113 changes: 113 additions & 0 deletions test/python/quantum_info/xx_decompose/test_weyl.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,113 @@
# This code is part of Qiskit.
#
# (C) Copyright IBM 2021
#
# This code is licensed under the Apache License, Version 2.0. You may
# obtain a copy of this license in the LICENSE.txt file in the root directory
# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
#
# Any modifications or derivative works of this code must retain this
# copyright notice, and modified files need to carry a notice indicating
# that they have been altered from the originals.

"""
Tests for qiskit-terra/qiskit/quantum_info/synthesis/xx_decompose/weyl.py .
"""

from itertools import permutations
import random
import unittest

import ddt
import numpy as np

from qiskit.quantum_info.operators import Operator
from qiskit.quantum_info.synthesis.xx_decompose.weyl import (
apply_reflection,
apply_shift,
canonical_rotation_circuit,
reflection_options,
shift_options,
)

from .utilities import canonical_matrix

EPSILON = 0.001


@ddt.ddt
class TestMonodromyWeyl(unittest.TestCase):
"""Check Weyl action routines."""

def __init__(self, *args, seed=42, **kwargs):
super().__init__(*args, **kwargs)
random.seed(seed)
np.random.seed(seed)

def test_reflections(self):
"""Check that reflection circuits behave as expected."""
for name in reflection_options:
coordinate = [random.random() for _ in range(3)]
reflected_coordinate, reflection_circuit, reflection_phase = apply_reflection(
name, coordinate
)
original_matrix = canonical_matrix(*coordinate)
reflected_matrix = canonical_matrix(*reflected_coordinate)
reflect_matrix = Operator(reflection_circuit).data
with np.printoptions(precision=3, suppress=True):
print(name)
print(reflect_matrix.conjugate().transpose(1, 0) @ original_matrix @ reflect_matrix)
print(reflected_matrix * reflection_phase)
self.assertTrue(
np.all(
np.abs(
reflect_matrix.conjugate().transpose(1, 0)
@ original_matrix
@ reflect_matrix
- reflected_matrix * reflection_phase
)
< EPSILON
)
)

def test_shifts(self):
"""Check that shift circuits behave as expected."""
for name in shift_options:
coordinate = [random.random() for _ in range(3)]
shifted_coordinate, shift_circuit, shift_phase = apply_shift(name, coordinate)
original_matrix = canonical_matrix(*coordinate)
shifted_matrix = canonical_matrix(*shifted_coordinate)
shift_matrix = Operator(shift_circuit).data
with np.printoptions(precision=3, suppress=True):
print(name)
print(original_matrix @ shift_matrix)
print(shifted_matrix * shift_phase)
self.assertTrue(
np.all(
np.abs(original_matrix @ shift_matrix - shifted_matrix * shift_phase) < EPSILON
)
)

def test_rotations(self):
"""Check that rotation circuits behave as expected."""
for permutation in permutations([0, 1, 2]):
coordinate = [random.random() for _ in range(3)]
rotation_circuit = canonical_rotation_circuit(permutation[0], permutation[1])
original_matrix = canonical_matrix(*coordinate)
rotation_matrix = Operator(rotation_circuit).data
rotated_matrix = canonical_matrix(
coordinate[permutation[0]],
coordinate[permutation[1]],
coordinate[permutation[2]],
)
self.assertTrue(
np.all(
np.abs(
rotation_matrix.conjugate().transpose(1, 0)
@ original_matrix
@ rotation_matrix
- rotated_matrix
)
< EPSILON
)
)
26 changes: 26 additions & 0 deletions test/python/quantum_info/xx_decompose/utilities.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,26 @@
# This code is part of Qiskit.
#
# (C) Copyright IBM 2021
#
# This code is licensed under the Apache License, Version 2.0. You may
# obtain a copy of this license in the LICENSE.txt file in the root directory
# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
#
# Any modifications or derivative works of this code must retain this
# copyright notice, and modified files need to carry a notice indicating
# that they have been altered from the originals.

"""
Helper functions for the XXDecomposer test suite.
"""

from qiskit.circuit.library import RXXGate, RYYGate, RZZGate


def canonical_matrix(a=0.0, b=0.0, c=0.0):
"""
Produces the matrix form of a "canonical operator"
exp(-i (a XX + b YY + c ZZ)) .
"""
return RXXGate(2 * a).to_matrix() @ RYYGate(2 * b).to_matrix() @ RZZGate(-2 * c).to_matrix()