# 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 __future__ import annotations from copy import copy from dataclasses import dataclass, field from itertools import combinations 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`. """ rng = np.random.default_rng(seed) 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.") rng.shuffle(paragraphs) for convex_subpolytope in paragraphs: sentences = convex_subpolytope.inequalities + convex_subpolytope.equalities if len(sentences) == 0: continue dimension = len(sentences[0]) - 1 rng.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`. """ 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. # A1, A1inv, A2, A2inv, A3, A3inv # 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`. """