diff --git a/.pylintrc b/.pylintrc index b11fca0a..a5fb49a1 100644 --- a/.pylintrc +++ b/.pylintrc @@ -33,7 +33,7 @@ unsafe-load-any-extension=no # A comma-separated list of package or module names from where C extensions may # be loaded. Extensions are loading into the active Python interpreter and may # run arbitrary code -extension-pkg-whitelist=retworkx, numpy, tweedledum +extension-pkg-whitelist=rustworkx, numpy, tweedledum [MESSAGES CONTROL] diff --git a/qiskit_qec/__init__.py b/qiskit_qec/__init__.py index 3727eb86..5810137a 100644 --- a/qiskit_qec/__init__.py +++ b/qiskit_qec/__init__.py @@ -15,6 +15,7 @@ from logging import NullHandler from . import ( + arithmetic, circuits, codes, decoders, diff --git a/qiskit_qec/arithmetic/__init__.py b/qiskit_qec/arithmetic/__init__.py new file mode 100644 index 00000000..e14a7298 --- /dev/null +++ b/qiskit_qec/arithmetic/__init__.py @@ -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. + +from .modn import gcd_ext, quo, div, ann, stab, unit diff --git a/qiskit_qec/arithmetic/modn.py b/qiskit_qec/arithmetic/modn.py new file mode 100644 index 00000000..739ea004 --- /dev/null +++ b/qiskit_qec/arithmetic/modn.py @@ -0,0 +1,444 @@ +# This code is part of Qiskit. +# +# (C) Copyright IBM 2017, 2020 +# +# 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. +# +# This code is adapted from XPFpackage: +# https://github.com/m-webster/XPFpackage, originally developed by Mark +# Webster. The original code is licensed under the GNU General Public +# License v3.0 and Mark Webster has given permission to use the code under +# the Apache License v2.0. + +"""Modular arithmetic Z/nZ.""" + +from typing import Tuple +import math +import numpy as np +from qiskit import QiskitError + + +def gcd_ext(a: int, b: int, n: int) -> Tuple[int, int, int, int, int]: + """Implements the extended Euclidean algorithm in the ring Z/nZ: for any + two integers a & b, find g, s, t, u, v that satisfy: + - g = gcd_n(a, b) = s * a + t * b, where gcd_n stands for the greatest + common divisor in the ring Z/nZ, s & t are called the Bezout coefficients + for a & b; + - (u * a + v * b) mod n = 0; + - (s * v - t * u) mod n = 1. + + Args: + a: input integer + b: input integer + n: modulus + + Returns: + g, s, t, u, v: g is the greatest common divisor of (a mod n) and (b mod n); + s, t, u, and v are integer coefficients satisfying (s*a + t*b) mod n = g, + (u*a + v*b) mod n = 0, (s*v - t*u) mod n = 1 + + Raises: + QiskitError: Input must be integers. + QiskitError: n must be a positive integer. + + Examples: + >>> gcd_ext(15, 6, 4) + (1, 1, -1, -2, 3) + + See Also: + _gcd_ext + """ + if not isinstance(a, (int, np.integer)) or not isinstance(b, (int, np.integer)): + raise QiskitError("Input must be integers.") + if not n > 0 or not isinstance(n, (int, np.integer)): + raise QiskitError("n must be a positive integer") + + return _gcd_ext(a, b, n) + + +def _gcd_ext(a: int, b: int, n: int) -> Tuple[int, int, int, int, int]: + """Implements the extended Euclidean algorithm in the ring Z/nZ: for any + two integers a & b, find g, s, t, u, v that satisfy: + - g = gcd_n(a, b) = s * a + t * b, where gcd_n stands for the greatest + common divisor in the ring Z/nZ, s & t are called the Bezout coefficients + for a & b; + - (u * a + v * b) mod n = 0; + - (s * v - t * u) mod n = 1. + + Args: + a: input integer + b: input integer + n: modulus + + Returns: + g_, s_, t_, u_, v_: g_ is the greatest common divisor of (a mod n) and + (b mod n); s_, t_, u_, and v_ are integer coefficients satisfying (s_ * + a + t_ * b) mod n = g_, (u_ * a + v_ * b) mod n = 0, (s_ * v_ - t_ * u_) + mod n = 1 + + Examples: + >>> _gcd_ext(15, 6, 4) + (1, 1, -1, -2, 3) + + See Also: + gcd_ext + """ + old_s_, s_ = 1, 0 + old_t_, t_ = 0, 1 + old_r_, r_ = a % n, b % n + + while r_ != 0: + q_ = old_r_ // r_ + old_r_, r_ = r_, old_r_ - q_ * r_ + old_s_, s_ = s_, old_s_ - q_ * s_ + old_t_, t_ = t_, old_t_ - q_ * t_ + + sgn = int(math.copysign(1, t_ * old_s_ - s_ * old_t_)) + u_, v_ = sgn * s_, sgn * t_ + g_, s_, t_ = old_r_, old_s_, old_t_ + + return (g_, s_, t_, u_, v_) + + +def quo(a: int, b: int, n: int) -> int: + """Computes the quotient of a/b in the ring Z/nZ, i.e. returns integer q + such that a = (b * q) mod n. Returns None if b mod n = 0. + + Args: + a: numerator + b: denominator + n: modulus + + Returns: + quotient of a/b in the ring Z/nZ + + Raises: + QiskitError: Input must be integers. + QiskitError: n must be a positive integer. + + Examples: + >>> quo(25, 5, 4) + 1 + + >>> quo(25, 8, 4) + None + + See Also: + _quo + """ + if not isinstance(a, (int, np.integer)) or not isinstance(b, (int, np.integer)): + raise QiskitError("Input must be integers.") + if not n > 0 or not isinstance(n, (int, np.integer)): + raise QiskitError("n must be a positive integer") + + return _quo(a, b, n) + + +def _quo(a: int, b: int, n: int) -> int: + """Computes the quotient of a/b in the ring Z/nZ, i.e. returns integer q + such that a = (b * q) mod n. Returns None if b mod n = 0. + + Args: + a: numerator + b: denominator + n: modulus + + Returns: + quotient of a/b in the ring Z/nZ + + Examples: + >>> _quo(25, 5, 4) + 1 + + >>> _quo(25, 8, 4) + None + + See Also: + quo + """ + a, b = a % n, b % n + if b == 0: + return None + return (a // b) % n + + +def div(a: int, b: int, n: int) -> int: + """Computes the divisor of a/b in the ring Z/nZ, i.e., returns integer d + such that b * d = a mod n. Returns None if no such d exists. + + Args: + a: numerator + b: denominator + n: modulus + + Returns: + divisor of a/b in the ring Z/nZ + + Raises: + QiskitError: Input must be integers. + QiskitError: n must be a positive integer. + + Examples: + >>> div(24, 8, 5) + 3 + + >>> div(24, 8, 4) + None + + >>> div(23, 10, 4) + None + + See Also: + _div + """ + if not isinstance(a, (int, np.integer)) or not isinstance(b, (int, np.integer)): + raise QiskitError("Input must be integers.") + if not n > 0 or not isinstance(n, (int, np.integer)): + raise QiskitError("n must be a positive integer") + + return _div(a, b, n) + + +def _div(a: int, b: int, n: int) -> int: + """Computes the divisor of a/b in the ring Z/nZ, i.e., returns integer d + such that b * d = a mod n. Returns None if no such d exists. + + Args: + a: numerator + b: denominator + n: modulus + + Returns: + divisor of a/b in the ring Z/nZ + + Examples: + >>> _div(24, 8, 5) + 3 + + >>> _div(24, 8, 4) + None + + >>> _div(23, 10, 4) + None + + See Also: + div + """ + a, b = a % n, b % n + if b == 0: + return None + gcd = math.gcd(b, n) + if a % gcd != 0: + return None + else: + r = a % b + while r > 0: + a += n + r = a % b + return a // b % n + + +def ann(a: int, n: int) -> int: + """Computes the annihilator of a in the ring Z/nZ, i.e., returns integer b + such that (b * a) mod N = 0. + + Args: + a: input integer + n: modulus + + Returns: + annihilator of a in the ring Z/nZ + + Raises: + QiskitError: Input must be integers. + QiskitError: n must be a positive integer. + + Examples: + >>> ann(10, 5) + 1 + + >>> ann(3, 5) + 0 + + See Also: + _ann + """ + if not isinstance(a, (int, np.integer)): + raise QiskitError("Input must be integers.") + if not n > 0 or not isinstance(n, (int, np.integer)): + raise QiskitError("n must be a positive integer") + + return _ann(a, n) + + +def _ann(a: int, n: int) -> int: + """Computes the annihilator of a in the ring Z/nZ, i.e., returns integer b + such that (b * a) mod n = 0. + + Args: + a: input integer + n: modulus + + Returns: + annihilator of a in the ring Z/nZ + + Examples: + >>> _ann(10, 5) + 1 + + >>> _ann(3, 5) + 0 + + See Also: + ann + """ + a = a % n + if a == 0: + return 1 + b = n // math.gcd(a, n) + return b % n + + +def stab(a: int, b: int, n: int) -> int: + """Returns a ring element c such that gcd(a + b * c, n) = gcd(a, b, n) in + the ring Z/nZ. + + Args: + a: input integer + b: input integer + n: modulus + + Returns: + ring element c such that gcd(a+b*c, N) = gcd(a, b, N) + + Raises: + QiskitError: Input must be integers. + QiskitError: n must be a positive integer. + + Examples: + >>> stab(25, 8, 6) + 0 + + >>> stab(24, 8, 6) + 1 + + See Also: + _stab + """ + if not isinstance(a, (int, np.integer)) or not isinstance(b, (int, np.integer)): + raise QiskitError("Input must be integers.") + if not n > 0 or not isinstance(n, (int, np.integer)): + raise QiskitError("n must be a positive integer") + + return _stab(a, b, n) + + +def _stab(a: int, b: int, n: int) -> int: + """Returns a ring element c such that gcd(a + b * c, n) = gcd(a, b, n) in + the ring Z/nZ. + + Args: + a: input integer + b: input integer + n: modulus + + Returns: + ring element c such that gcd(a+b*c, N) = gcd(a, b, N) + + Examples: + >>> _stab(25, 8, 6) + 0 + + >>> _stab(24, 8, 6) + 1 + + See Also: + stab + """ + a, b = a % n, b % n + gcd = math.gcd(math.gcd(a, b), n) + n_old = n + a, n = a // gcd, n // gcd + if n == 0: + c = 0 + else: + a = a % n + if a == 0: + c = 1 + else: + r = int(math.ceil(math.log2(math.log2(n)))) if n > 1 else 1 + for _ in range(r): + a = a * a % n + c = n // math.gcd(a, n) + return c % n_old + + +def unit(a: int, n: int) -> int: + """Computes a unit u such that for element a in the ring Z/nZ, i.e., + (a * u) mod n = gcd(a, n). + + Args: + a: input integer + n: modulus + + Returns: + unit of a in the ring Z/nZ + + Raises: + QiskitError: Input must be integers. + QiskitError: n must be a positive integer. + + Examples: + >>> unit(10, 5) + 1 + + >>> unit(3, 5) + 2 + + See Also: + _unit + """ + if not isinstance(a, (int, np.integer)): + raise QiskitError("Input must be integers.") + if not n > 0 or not isinstance(n, (int, np.integer)): + raise QiskitError("n must be a positive integer") + + return _unit(a, n) + + +def _unit(a: int, n: int) -> int: + """Computes a unit u such that for element a in the ring Z/nZ, i.e., + (a * u) mod n = gcd(a, n). + + Args: + a: input integer + n: modulus + + Returns: + unit of a in the ring Z/nZ + + Examples: + >>> _unit(10, 5) + 1 + + >>> _unit(3, 5) + 2 + + See Also: + unit + """ + a = a % n + if a == 0: + return 1 + gcd = math.gcd(a, n) + d = div(gcd, a, n) + if gcd == 1: + return d + c = stab(d, n // gcd, n) + return (d + c * n // gcd) % n diff --git a/qiskit_qec/circuits/repetition_code.py b/qiskit_qec/circuits/repetition_code.py index bbfc96be..9c73e261 100644 --- a/qiskit_qec/circuits/repetition_code.py +++ b/qiskit_qec/circuits/repetition_code.py @@ -18,7 +18,7 @@ from typing import List, Optional, Tuple import numpy as np -import retworkx as rx +import rustworkx as rx from qiskit import ClassicalRegister, QuantumCircuit, QuantumRegister, transpile from qiskit.circuit.library import XGate, RZGate diff --git a/qiskit_qec/decoders/base_matcher.py b/qiskit_qec/decoders/base_matcher.py index d2d80398..75bc8811 100644 --- a/qiskit_qec/decoders/base_matcher.py +++ b/qiskit_qec/decoders/base_matcher.py @@ -1,7 +1,7 @@ """Base matching object.""" from abc import ABC, abstractmethod from typing import List, Tuple, Dict, Set -import retworkx as rx +import rustworkx as rx class BaseMatcher(ABC): diff --git a/qiskit_qec/decoders/circuit_matching_decoder.py b/qiskit_qec/decoders/circuit_matching_decoder.py index 3cf9bd11..fbde457e 100644 --- a/qiskit_qec/decoders/circuit_matching_decoder.py +++ b/qiskit_qec/decoders/circuit_matching_decoder.py @@ -7,12 +7,12 @@ from typing import Dict, List, Tuple from sympy import Poly, Symbol, symbols -import retworkx as rx +import rustworkx as rx from qiskit import QuantumCircuit from qiskit_qec.analysis.faultenumerator import FaultEnumerator from qiskit_qec.decoders.decoding_graph import CSSDecodingGraph, DecodingGraph from qiskit_qec.decoders.pymatching_matcher import PyMatchingMatcher -from qiskit_qec.decoders.retworkx_matcher import RetworkXMatcher +from qiskit_qec.decoders.rustworkx_matcher import RustworkxMatcher from qiskit_qec.decoders.temp_code_util import temp_gauge_products, temp_syndrome from qiskit_qec.exceptions import QiskitQECError from qiskit_qec.noise.paulinoisemodel import PauliNoiseModel @@ -21,7 +21,7 @@ class CircuitModelMatchingDecoder(ABC): """Matching decoder for circuit noise.""" - METHOD_RETWORKX: str = "retworkx" + METHOD_RETWORKX: str = "rustworkx" METHOD_PYMATCHING: str = "pymatching" AVAILABLE_METHODS = {METHOD_RETWORKX, METHOD_PYMATCHING} @@ -64,7 +64,7 @@ def __init__( blocks : number of measurement blocks method : matching implementation uniform : use same edge weight everywhere? - annotate : for retworkx method, compute self.matcher.annotated_graph + annotate : for rustworkx method, compute self.matcher.annotated_graph """ self.n = n self.css_x_gauge_ops = css_x_gauge_ops @@ -94,7 +94,7 @@ def __init__( if self.method == self.METHOD_PYMATCHING: self.matcher = PyMatchingMatcher() else: - self.matcher = RetworkXMatcher(annotate) + self.matcher = RustworkxMatcher(annotate) self.z_gauge_products = temp_gauge_products(self.css_z_stabilizer_ops, self.css_z_gauge_ops) self.x_gauge_products = temp_gauge_products(self.css_x_stabilizer_ops, self.css_x_gauge_ops) @@ -136,7 +136,7 @@ def __init__( self.edge_weight_polynomials = {} self.symbols = None if not self.uniform: - fe = FaultEnumerator(circuit, order=1, method="stabilizer", model=self.model) + fe = FaultEnumerator(circuit, order=1, method="propagator", model=self.model) self.event_map = self._enumerate_events( self.css_x_gauge_ops, self.css_x_stabilizer_ops, @@ -304,18 +304,9 @@ def _revise_decoding_graph( for source, target in graph.edge_list(): edge_data = graph.get_edge_data(source, target) if "weight_poly" not in edge_data and edge_data["weight"] != 0: - if ( - "is_boundary" in graph.nodes()[source] and graph.nodes()[source]["is_boundary"] - ) or ( - "is_boundary" in graph.nodes()[target] and graph.nodes()[target]["is_boundary"] - ): - # TODO: we could consider removing the edge - edge_data["weight"] = 1000000 - logging.info("increase edge weight (%d, %d)", source, target) - else: - # Otherwise, remove the edge - remove_list.append((source, target)) - logging.info("remove edge (%d, %d)", source, target) + # Remove the edge + remove_list.append((source, target)) + logging.info("remove edge (%d, %d)", source, target) graph.remove_edges_from(remove_list) return graph diff --git a/qiskit_qec/decoders/decoding_graph.py b/qiskit_qec/decoders/decoding_graph.py index 1cbcf1e7..c26958a8 100644 --- a/qiskit_qec/decoders/decoding_graph.py +++ b/qiskit_qec/decoders/decoding_graph.py @@ -22,7 +22,7 @@ from typing import List, Tuple import numpy as np -import retworkx as rx +import rustworkx as rx from qiskit_qec.analysis.faultenumerator import FaultEnumerator @@ -50,17 +50,18 @@ def _make_syndrome_graph(self): S = rx.PyGraph(multigraph=False) self.hyperedges = [] - # get the circuit used as the base case - if isinstance(self.code.circuit, dict): - if "base" not in dir(self.code): - base = "0" + if self.code is not None: + + # get the circuit used as the base case + if isinstance(self.code.circuit, dict): + if "base" not in dir(self.code): + base = "0" + else: + base = self.code.base + qc = self.code.circuit[base] else: - base = self.code.base - qc = self.code.circuit[base] - else: - qc = self.code.circuit + qc = self.code.circuit - if self.code is not None: fe = FaultEnumerator(qc, method="stabilizer") blocks = list(fe.generate_blocks()) fault_paths = list(itertools.chain(*blocks)) diff --git a/qiskit_qec/decoders/pymatching_matcher.py b/qiskit_qec/decoders/pymatching_matcher.py index 82167c5f..197b0e65 100644 --- a/qiskit_qec/decoders/pymatching_matcher.py +++ b/qiskit_qec/decoders/pymatching_matcher.py @@ -3,7 +3,7 @@ from typing import List, Tuple, Dict, Set import logging -import retworkx as rx +import rustworkx as rx from pymatching import Matching from qiskit_qec.exceptions import QiskitQECError @@ -15,7 +15,7 @@ class PyMatchingMatcher(BaseMatcher): """Matching subroutines using PyMatching. - The input retworkx graph is expected to have the following properties: + The input rustworkx graph is expected to have the following properties: edge["weight"] : real edge weight edge["qubits"] : list of qubit ids associated to edge vertex["is_boundary"] : bool, true if boundary node diff --git a/qiskit_qec/decoders/retworkx_matcher.py b/qiskit_qec/decoders/retworkx_matcher.py index e345ea62..d3deb026 100644 --- a/qiskit_qec/decoders/retworkx_matcher.py +++ b/qiskit_qec/decoders/retworkx_matcher.py @@ -4,14 +4,14 @@ from copy import deepcopy from typing import Dict, List, Set, Tuple -import retworkx as rx +import rustworkx as rx from qiskit_qec.decoders.base_matcher import BaseMatcher -class RetworkXMatcher(BaseMatcher): - """Matching subroutines using retworkx. +class RustworkxMatcher(BaseMatcher): + """Matching subroutines using rustworkx. - The input retworkx graph is expected to have the following properties: + The input rustworkx graph is expected to have the following properties: edge["weight"] : real edge weight edge["measurement_error"] : bool, true if edge corresponds to measurement error edge["qubits"] : list of qubit ids associated to edge @@ -98,7 +98,7 @@ def _error_chain_from_vertex_path( Examine the edges along the path to extract the error chain. Store error chains as sets and merge using symmetric difference. - The vertex_path is a list of retworkx node indices. + The vertex_path is a list of rustworkx node indices. """ qubit_errors = set([]) measurement_errors = set([]) diff --git a/qiskit_qec/decoders/rustworkx_matcher.py b/qiskit_qec/decoders/rustworkx_matcher.py new file mode 100644 index 00000000..51486518 --- /dev/null +++ b/qiskit_qec/decoders/rustworkx_matcher.py @@ -0,0 +1,179 @@ +"""rustworkx matching object.""" + +import logging +from copy import deepcopy +from typing import Dict, List, Set, Tuple + +import rustworkx as rx +from qiskit_qec.decoders.base_matcher import BaseMatcher + + +class RustworkxMatcher(BaseMatcher): + """Matching subroutines using rustworkx. + + The input rustworkx graph is expected to have the following properties: + edge["weight"] : real edge weight + edge["measurement_error"] : bool, true if edge corresponds to measurement error + edge["qubits"] : list of qubit ids associated to edge + vertex["qubits"] : qubit ids involved in gauge measurement + vertex["time"] : integer time step of gauge measurement + + The annotated graph will also have "highlighted" properties on edges and vertices. + """ + + def __init__(self, annotate: bool = False): + """Create the matcher.""" + self.length = {} + self.path = {} + self.annotate = annotate + self.annotated_graph = None + super().__init__() + + def preprocess(self, graph: rx.PyGraph): + """Compute shortest paths between vertex pairs in decoding graph. + + Updates sets self.length and self.path. + """ + edge_cost_fn = lambda edge: edge["weight"] + length = rx.all_pairs_dijkstra_path_lengths(graph, edge_cost_fn) + self.length = {s: dict(length[s]) for s in length} + path = rx.all_pairs_dijkstra_shortest_paths(graph, edge_cost_fn) + self.path = {s: {t: list(path[s][t]) for t in path[s]} for s in path} + + def find_errors( + self, + graph: rx.PyGraph, + idxmap: Dict[Tuple[int, List[int]], int], + highlighted: List[Tuple[int, Tuple[int]]], + ) -> Tuple[Set[int], Set[Tuple[int, Tuple[int]]]]: + """Process a set of highlighted vertices and return error locations. + + Be sure to have called recompute_paths if needed. + """ + matching = self._compute_matching(idxmap, highlighted) + logging.info("process: matching = %s", matching) + qubit_errors, measurement_errors = self._compute_error_correction( + graph, idxmap, matching, highlighted + ) + logging.info("process: qubit_errors = %s", qubit_errors) + logging.debug("process: measurement_errors = %s", measurement_errors) + return qubit_errors, measurement_errors + + def _compute_matching( + self, + idxmap: Dict[Tuple[int, List[int]], int], + highlighted: List[Tuple[int, Tuple[int]]], + ) -> Set[Tuple[int, int]]: + """Compute a min. weight perfect matching of highlighted vertices. + + highlighted is a list of highlighted vertices given as tuples + (t, qubit_set). + Return the matching. + """ + gm = rx.PyGraph(multigraph=False) # matching graph + idx = 0 # vertex index in matching graph + midxmap = {} # map from (t, qubit_tuple) to vertex index + for v in highlighted: + gm.add_node({"dvertex": v}) + midxmap[v] = idx + idx += 1 + for i, high_i in enumerate(highlighted): + for j in range(i + 1, len(highlighted)): + vi = midxmap[high_i] + vj = midxmap[highlighted[j]] + vip = idxmap[high_i] + vjp = idxmap[highlighted[j]] + gm.add_edge(vi, vj, {"weight": -self.length[vip][vjp]}) + + def weight_fn(edge): + return int(edge["weight"]) + + matching = rx.max_weight_matching(gm, max_cardinality=True, weight_fn=weight_fn) + return matching + + def _error_chain_from_vertex_path( + self, graph: rx.PyGraph, vertex_path: List[int] + ) -> Tuple[Set[int], Set[Tuple[int, Tuple[int]]]]: + """Return a chain of qubit and measurement errors from a vertex path. + + Examine the edges along the path to extract the error chain. + Store error chains as sets and merge using symmetric difference. + The vertex_path is a list of rustworkx node indices. + """ + qubit_errors = set([]) + measurement_errors = set([]) + logging.debug("_error_chain_from_vertex_path %s", vertex_path) + for i in range(len(vertex_path) - 1): + v0 = vertex_path[i] + v1 = vertex_path[i + 1] + if graph.get_edge_data(v0, v1)["measurement_error"] == 1: + measurement_errors ^= set( + [(graph.nodes()[v0]["time"], tuple(graph.nodes()[v0]["qubits"]))] + ) + qubit_errors ^= set(graph.get_edge_data(v0, v1)["qubits"]) + logging.debug( + "_error_chain_for_vertex_path q = %s, m = %s", + qubit_errors, + measurement_errors, + ) + return qubit_errors, measurement_errors + + def _compute_error_correction( + self, + graph: rx.PyGraph, + idxmap: Dict[Tuple[int, List[int]], int], + matching: Set[Tuple[int, int]], + highlighted: List[Tuple[int, Tuple[int]]], + ) -> Tuple[Set[int], Set[Tuple[int, Tuple[int]]]]: + """Compute the qubit and measurement corrections. + + graph : the decoding graph + idxmap : maps (t, qubit_idx) to vertex index + matching : perfect matching computed by _compute_matching + highlighted : list of highlighted vertices + + Returns a tuple of sets, (qubit_errors, measurement_errors) where + qubit_errors contains the indices of qubits with errors and + measurement_errors contains tuples (t, qubit_set) indicating + failed measurements. + """ + used_paths = [] + qubit_errors = set([]) + measurement_errors = set([]) + for p in matching: + v0 = idxmap[highlighted[p[0]]] + v1 = idxmap[highlighted[p[1]]] + # Use the shortest paths between the matched vertices to + # identify all of the qubits in the error chains + path = self.path[v0][v1] + q, m = self._error_chain_from_vertex_path(graph, path) + # Add the error chains modulo two to get the total correction + # (uses set symmetric difference) + qubit_errors ^= q + measurement_errors ^= m + used_paths.append(path) + if self.annotate: + self.annotated_graph = self._make_annotated_graph(graph, used_paths) + return qubit_errors, measurement_errors + + def _make_annotated_graph(self, gin: rx.PyGraph, paths: List[List[int]]) -> rx.PyGraph: + """Highlight the vertex paths and return annotated graph. + + gin : decoding graph + paths : list of vertex paths, each given as a list of + vertex indices in the decoding graph. + """ + graph = deepcopy(gin) + for path in paths: + # Highlight the endpoints of the path + for i in [0, -1]: + graph.nodes()[path[i]]["highlighted"] = True + # Highlight the edges along the path + for i in range(len(path) - 1): + try: + idx = list(graph.edge_list()).index((path[i], path[i + 1])) + except ValueError: + idx = list(graph.edge_list()).index((path[i + 1], path[i])) + edge = graph.edges()[idx] + edge["highlighted"] = True + return graph diff --git a/qiskit_qec/decoders/temp_graph_util.py b/qiskit_qec/decoders/temp_graph_util.py index cd74bed4..15f341d9 100644 --- a/qiskit_qec/decoders/temp_graph_util.py +++ b/qiskit_qec/decoders/temp_graph_util.py @@ -1,11 +1,11 @@ """Temporary module with methods for graphs.""" import json import networkx as nx -import retworkx as rx +import rustworkx as rx def ret2net(graph: rx.PyGraph): - """Convert retworkx graph to equivalent networkx graph.""" + """Convert rustworkx graph to equivalent networkx graph.""" nx_graph = nx.Graph() for j, node in enumerate(graph.nodes()): nx_graph.add_node(j) diff --git a/qiskit_qec/linear/__init__.py b/qiskit_qec/linear/__init__.py index 6e0500a9..78194bc4 100644 --- a/qiskit_qec/linear/__init__.py +++ b/qiskit_qec/linear/__init__.py @@ -62,7 +62,16 @@ from .bit.bit import Bit from .smatrix_api.smatrix import SMatrix -from .matrix import create_lambda_matrix, augment_mat, rref, rank, rref_complete +from .matrix import ( + create_lambda_matrix, + augment_mat, + rref, + rank, + rref_complete, + do_row_op, + howell, + howell_complete, +) from .symplectic import ( all_commute, symplectic_product, diff --git a/qiskit_qec/linear/matrix.py b/qiskit_qec/linear/matrix.py index 4198f9d1..606b22cc 100644 --- a/qiskit_qec/linear/matrix.py +++ b/qiskit_qec/linear/matrix.py @@ -12,13 +12,14 @@ """Matrix ops.""" -from typing import List, Tuple +from typing import List, Tuple, Union import numpy as np from qiskit import QiskitError +from qiskit_qec.arithmetic.modn import gcd_ext, quo, ann, unit def create_lambda_matrix(n: int) -> np.ndarray: - r"""Creates the GF(2) Lambda matrix + r"""Creates the GF(2) Lambda matrix. [0 I_n] [I_n 0] @@ -52,7 +53,7 @@ def create_lambda_matrix(n: int) -> np.ndarray: def _create_lambda_matrix(n: int): - r"""Creates the GF(2) Lambda matrix + r"""Creates the GF(2) Lambda matrix. [0 I_n] [I_n 0] @@ -82,8 +83,8 @@ def _create_lambda_matrix(n: int): def augment_mat(matrix: np.ndarray, pos: str = "right") -> np.ndarray: """Augments a given matrix with the identify matrix. The position of the - identity matrix is given by the optional position argument - left -> [M | I] + identity matrix is given by the optional position argument left -> [M | + I] right -> [I | M] @@ -135,7 +136,7 @@ def augment_mat(matrix: np.ndarray, pos: str = "right") -> np.ndarray: def _augment_mat(matrix: np.array, pos: str) -> np.ndarray: """Augments a given matrix with the identify matrix. The position of the - identity matrix is given by the optional position argument + identity matrix is given by the optional position argument. left -> [M | I] @@ -250,7 +251,7 @@ def _rref(matrix: np.array) -> np.ndarray: def rank(matrix: np.ndarray) -> int: - """Computes the rank of the GF(2) matrix + """Computes the rank of the GF(2) matrix. Args: matrix: Input GF(2) matrix @@ -275,7 +276,7 @@ def rank(matrix: np.ndarray) -> int: def _rank(matrix: np.ndarray) -> int: - """Computes the rank of the GF(2) matrix + """Computes the rank of the GF(2) matrix. Args: matrix: Input GF(2) matrix @@ -299,8 +300,8 @@ def _rank(matrix: np.ndarray) -> int: def rref_complete(matrix: np.ndarray) -> Tuple[List[int], np.ndarray, np.ndarray, int]: - """Computes the Row Reduced Echelon Form for a GF(2) matrix as well - as pivots, transformation matrix and rank. + """Computes the Row Reduced Echelon Form for a GF(2) matrix as well as + pivots, transformation matrix and rank. Args: matrix : Input GF(2) matrix. @@ -354,8 +355,8 @@ def rref_complete(matrix: np.ndarray) -> Tuple[List[int], np.ndarray, np.ndarray def _rref_complete(matrix) -> Tuple[List[int], np.ndarray, np.ndarray, int]: - """Computes the Row Reduced Echelon Form for a GF(2) matrix as well - as pivots, transformation matrix and rank. + """Computes the Row Reduced Echelon Form for a GF(2) matrix as well as + pivots, transformation matrix and rank. Args: matrix : Input GF(2) matrix. @@ -493,3 +494,419 @@ def istack(mat: np.ndarray, size: int, interleave: bool = False) -> np.ndarray: if interleave: return np.hstack(size * [mat]).reshape((size * len(mat),) + mat.shape[1:]) return np.vstack(size * [mat]).reshape((size * len(mat),) + mat.shape[1:]) + + +# --------------------------------------------------------------- +# Span-preserving row operations and Howell matrix form +# This code is adapted from XPFpackage: +# https://github.com/m-webster/XPFpackage, originally developed by Mark +# Webster. The original code is licensed under the GNU General Public +# License v3.0 and Mark Webster has given permission to use the code under +# the Apache License v2.0. + + +def do_row_op(mat: np.ndarray, row_op: Tuple[str, List[int], List[int]], n: int) -> np.ndarray: + """Performs span-preserving row operations on a matrix in the ring Z/NZ. + + These include: + - Swap two rows mat[i] and mat[j]; + - Multiply a row by a scalar c (valid only when c is a unit); + - Add a multiple of one row to another; + - Append the product of a row by a scalar c to the end of matrix (used + when c is a zero divisor); + - Update two rows by multiplying by a full-rank 2x2 matrix. + + Args: + mat: input matrix + row_op: tuple (op, rows, coeff), where op is the operation to be + performed, rows is a list of row indices, and coeff is a list of + coefficients for the operation (empty list if op is 'swap') + n: modulus + + Returns: + matrix after performing the row operation + + Raises: + QiskitError: Modulus N must be a positive integer. + QiskitError: Input matrix must be a 2D array. + QiskitError: Row operation must be a valid operation ("swap", "unit", "add", "append", "update"). + QiskitError: Row operation must involve valid row indices for the input matrix. + QiskitError: Swap operation must involve two rows. + QiskitError: Unit operation must involve one row and one coefficient. + QiskitError: Add operation must involve two rows and one coefficient. + QiskitError: Append operation must involve one row and one coefficient. + QiskitError: Update operation must involve two rows and one 2x2 matrix (four coefficients). + + Examples: + >>> mat = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]]) + >>> do_row_op(mat, ('swap', [0, 2], []), 4) + array([[7, 8, 9], + [4, 5, 6], + [1, 2, 3]]) + + >>> do_row_op(mat, ('unit', [0], [2]), 4) + array([[2, 0, 2], + [4, 5, 6], + [7, 8, 9]]) + """ + if not isinstance(n, (int, np.integer)) or n <= 0: + raise QiskitError("Modulus N must be a positive integer.") + mat = np.array(mat, dtype=int) + if not mat.ndim == 2: + raise QiskitError("Input matrix must be a 2D array") + + # get number of rows + nrows = mat.shape[0] + op, rows, coeff = row_op + + if op not in ["swap", "unit", "add", "append", "update"]: + raise QiskitError( + 'Row operation must be a valid operation ("swap", "unit", "add", "append", "update").' + ) + if any((row >= nrows for row in rows)): + raise QiskitError("Row operation must involve valid row indices for the input matrix.") + if op == "swap": + if len(rows) != 2: + raise QiskitError("Swap operation must involve two rows.") + mat = _swap_rows(mat, rows) + elif op == "unit": + if len(rows) != 1 or len(coeff) != 1: + raise QiskitError("Unit operation must involve one row and one coefficient.") + mat = _multiply_unit(mat, rows[0], coeff[0], n) + elif op == "add": + if len(rows) != 2 or len(coeff) != 1: + raise QiskitError("Add operation must involve two rows and one coefficient.") + mat = _add_rows(mat, rows, coeff[0], n) + elif op == "append": + if len(rows) != 1 or len(coeff) != 1: + raise QiskitError("Append operation must involve one row and one coefficient.") + mat = _append_row(mat, rows[0], coeff[0], n) + elif op == "update": + if len(rows) != 2 or len(coeff) != 4: + raise QiskitError( + "Update operation must involve two rows and one 2x2 matrix (four coefficients)." + ) + mat = _update_rows(mat, rows, coeff, n) + + return mat + + +def _swap_rows(mat: np.ndarray, rows: Union[list, np.ndarray]) -> np.ndarray: + """Swaps two rows of a matrix. + + Args: + mat: input matrix + rows: list of indices of the two rows to swap + + Returns: + matrix with rows swapped + + Examples: + >>> mat = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]]) + >>> _swap_rows(mat, [0, 2]) + array([[7, 8, 9], + [4, 5, 6], + [1, 2, 3]]) + """ + mat[[rows[0], rows[1]]] = mat[[rows[1], rows[0]]] + return mat + + +def _multiply_unit(mat: np.ndarray, row: int, c: int, n: int) -> np.ndarray: + """Multiplies a row of a matrix by a scalar c (valid only when c is a unit) + in the ring Z/nZ. + + Args: + mat: input matrix + row: index of row to multiply + c: scalar to multiply by + n: modulus + + Returns: + matrix with row multiplied by scalar + + Examples: + >>> mat = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]]) + >>> _multiply_unit(mat, 0, 2, 4) + array([[2, 0, 2], + [4, 5, 6], + [7, 8, 9]]) + """ + mat[row] = np.mod(c * mat[row], n) + return mat + + +def _add_rows(mat: np.ndarray, rows: Union[list, np.ndarray], c: int, n: int) -> np.ndarray: + """Adds a multiple of one row to another of a matrix in the ring Z/nZ. + + Args: + mat: input matrix + rows: list of indices of the two rows in action + c: scalar to multiply by + n: modulus + + Returns: + matrix with rows added + + Examples: + >>> mat = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]]) + >>> _add_rows(mat, [0, 2], 2, 4) + array([[3, 2, 1], + [4, 5, 6], + [7, 8, 9]]) + """ + mat[rows[0]] = np.mod(mat[rows[0]] + c * mat[rows[1]], n) + return mat + + +def _append_row(mat: np.ndarray, row: int, c: int, n: int) -> np.ndarray: + """Appends the product of a row by a scalar c to the end of matrix (used + when c is a zero divisor) in the ring Z/nZ. + + Args: + mat: input matrix + row: index of row to multiply + c: scalar to multiply by + n: modulus + + Returns: + matrix with row multiplied by scalar and appended + + Examples: + >>> mat = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]]) + >>> _append_row(mat, 0, 2, 4) + array([[1, 2, 3], + [4, 5, 6], + [7, 8, 9], + [2, 0, 2]]) + """ + mat = np.vstack((mat, np.mod(c * mat[row], n))) + return mat + + +def _update_rows( + mat: np.ndarray, rows: Union[list, np.ndarray], c: Union[list, np.ndarray], n: int +) -> np.ndarray: + """Updates two rows by multiplying by a full-rank 2x2 matrix in the ring + Z/nZ. + + Args: + mat: input matrix + rows: list of indices of the two rows in action + c: components of the 2x2 matrix to multiply by + n: modulus + + Returns: + matrix with rows updated + + Examples: + >>> mat = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]]) + >>> _update_rows(mat, [0, 1], [2, 1, 3, 4], 4) + array([[2, 1, 0], + [3, 2, 1], + [7, 8, 9]]) + """ + r1 = np.mod(c[0] * mat[rows[0]] + c[1] * mat[rows[1]], n) + r2 = np.mod(c[2] * mat[rows[0]] + c[3] * mat[rows[1]], n) + mat[rows[0]], mat[rows[1]] = r1, r2 + return mat + + +def howell(mat: np.ndarray, n: int) -> np.ndarray: + """Computes the Howell form of a matrix in the ring Z/NZ. + + Args: + mat: input matrix + n: modulus + + Returns: + Howell form of input matrix + + Examples: + >>> mat = numpy.array([[8, 5, 5], + [0, 9, 8], + [0, 0, 10]]) + >>> n = 12 + >>> howell_mat = howell(mat, n) + >>> howell_mat + array([[4, 1, 0], + [0, 3, 0], + [0, 0, 1]]) + + See Also: + _howell, howell_complete, _howell_complete + """ + howell_mat, _, _ = howell_complete(mat, n) + return howell_mat + + +def _howell(mat: np.ndarray, n: int) -> np.ndarray: + """Computes the Howell form of a matrix in the ring Z/nZ. + + Args: + mat: input matrix + n: modulus + + Returns: + Howell form of input matrix + + Examples: + >>> mat = numpy.array([[8, 5, 5], + [0, 9, 8], + [0, 0, 10]]) + >>> n = 12 + >>> howell_mat = _howell(mat, n) + >>> howell_mat + array([[4, 1, 0], + [0, 3, 0], + [0, 0, 1]]) + + See Also: + howell, howell_complete, _howell_complete + """ + howell_mat, _, _ = _howell_complete(mat, n) + return howell_mat + + +def howell_complete(mat: np.ndarray, n: int) -> Tuple[np.ndarray, np.ndarray, np.ndarray]: + """Computes the Howell form of a matrix in the ring Z/nZ, the corresponding + transformation matrix, and the kernel matrix. + + Args: + mat: input matrix + n: modulus + + Returns: + howell_mat: Howell form of mat + transform_mat: transformation matrix (transform_mat @ mat = howell_mat) + kernel_mat: kernel of mat (kernel_mat @ mat = 0) + + Raises: + QiskitError: Modulus N must be a positive integer + QiskitError: Input matrix must be a 2D array + + Examples: + >>> mat = numpy.array([[8, 5, 5], + [0, 9, 8], + [0, 0, 10]]) + >>> n = 12 + >>> howell_mat, transform_mat, kernel_mat = howell_complete(mat, n) + >>> howell_mat + array([[4, 1, 0], + [0, 3, 0], + [0, 0, 1]]) + >>> transform_mat + array([[8, 1, 0], + [0, 7, 4], + [9, 3, 4]]) + >>> kernel_mat + array([[6, 6, 6], + [0, 4, 4]]) + + See Also: + _howell_complete, howell, _howell + """ + if not n > 0 or not isinstance(n, (int, np.integer)): + raise QiskitError("Modulus N must be a positive integer") + mat = np.array(mat, dtype=int) + if not mat.ndim == 2: + raise QiskitError("Input matrix must be a 2D array") + + return _howell_complete(mat, n) + + +def _howell_complete(mat: np.ndarray, n: int) -> Tuple[np.ndarray, np.ndarray, np.ndarray]: + """Computes the Howell form of a matrix in the ring Z/nZ, the corresponding + transformation matrix, and the kernel. + + Args: + mat: input matrix + n: modulus + + Returns: + howell_mat: Howell form of mat + transform_mat: transformation matrix (transform_mat @ mat = howell_mat) + kernel_mat: kernel of mat (kernel_mat @ mat = 0) + + Examples: + >>> mat = numpy.array([[8, 5, 5], + [0, 9, 8], + [0, 0, 10]]) + >>> n = 12 + >>> howell_mat, transform_mat, kernel_mat = _howell_complete(mat, n) + >>> howell_mat + array([[4, 1, 0], + [0, 3, 0], + [0, 0, 1]]) + >>> transform_mat + array([[8, 1, 0], + [0, 7, 4], + [9, 3, 4]]) + >>> kernel_mat + array([[6, 6, 3], + [0, 4, 4]]) + + See Also: + howell_complete, howell, _howell + """ + howell_mat = mat.copy() + transform_mat = np.eye(howell_mat.shape[0], dtype=int) + n_row, n_col = howell_mat.shape + + # set row index to 0 + r = 0 + # going through each column + for c in range(n_col): + # find j such that howell_mat[j, c] > 0 + j = r + while j < n_row and howell_mat[j, c] == 0: + j += 1 + if j < n_row: + # found j: if j > r, swap rows r and j + if j > r: + howell_mat = _swap_rows(howell_mat, [r, j]) + transform_mat = _swap_rows(transform_mat, [r, j]) + + # multiply row r by a unit to ensure that howell_mat[r, c] is a minimal + # representative + x = unit(howell_mat[r, c], n) + if x > 1: + howell_mat = _multiply_unit(howell_mat, r, x, n) + transform_mat = _multiply_unit(transform_mat, r, x, n) + + # eliminate entries in column c below row r + for i in range(r + 1, n_row): + if howell_mat[i, c] % n > 0: + _, s_, t_, u_, v_ = gcd_ext(howell_mat[r, c], howell_mat[i, c], n) + howell_mat = _update_rows(howell_mat, [r, i], [s_, t_, u_, v_], n) + transform_mat = _update_rows(transform_mat, [r, i], [s_, t_, u_, v_], n) + + # ensure entries in column c above row r are less than + # howell_mat[r, c] + b = howell_mat[r, c] + for i in range(r): + if howell_mat[i, c] >= b: + x = quo(howell_mat[i, c], b, n) + howell_mat = _add_rows(howell_mat, [i, r], -1 * x, n) + transform_mat = _add_rows(transform_mat, [i, r], -1 * x, n) + + # if b = howell_mat[r, c] is a zero divisor, find the annihilator x that + # eliminates howell_mat[r, c] and append a new row x * + # howell_mat[r] + x = ann(b, n) + if x > 0: + howell_mat = _append_row(howell_mat, r, x, n) + transform_mat = _append_row(transform_mat, r, x, n) + n_row = len(howell_mat) + r += 1 + + # remove rows of zeros + howell_mat = howell_mat[howell_mat.any(axis=1)] + + # compute the transformation matrix and kernel matrix + k = len(howell_mat) + kernel_mat = transform_mat[k:, :] + kernel_mat = kernel_mat[kernel_mat.any(axis=1)] + transform_mat = transform_mat[:k, :] + + return howell_mat, transform_mat, kernel_mat diff --git a/qiskit_qec/operators/__init__.py b/qiskit_qec/operators/__init__.py index f58c9258..75417c02 100644 --- a/qiskit_qec/operators/__init__.py +++ b/qiskit_qec/operators/__init__.py @@ -27,8 +27,12 @@ PauliList Pauli BasePauli + BaseXPPauli """ from .base_pauli import BasePauli from .pauli import Pauli from .pauli_list import PauliList +from .base_xp_pauli import BaseXPPauli +from .xp_pauli import XPPauli +from .xp_pauli_list import XPPauliList diff --git a/qiskit_qec/operators/base_pauli.py b/qiskit_qec/operators/base_pauli.py index 54429362..47127478 100644 --- a/qiskit_qec/operators/base_pauli.py +++ b/qiskit_qec/operators/base_pauli.py @@ -74,7 +74,7 @@ def __init__( Args: matrix: Input GF(2) symplectic matrix - phase_exp (optional): Phase exponent vector for imput matrix. A value of None will + phase_exp (optional): Phase exponent vector for input matrix. A value of None will result in an a complex coefficients of 1 for each Pauli operator. Defaults to None. order: Set to 'xz' or 'zx'. Defines which side the x and z parts of the input matrix @@ -773,7 +773,7 @@ def to_label( squeeze: bool = True, index_str: str = "", ) -> Union[str, List[str]]: - """Returns the string representatiojn for a Pauli or Paulis. + """Returns the string representation for a Pauli or Paulis. Args: output_pauli_encoding (optional): Encoding used to represent phases. diff --git a/qiskit_qec/operators/base_xp_pauli.py b/qiskit_qec/operators/base_xp_pauli.py new file mode 100644 index 00000000..d003ea17 --- /dev/null +++ b/qiskit_qec/operators/base_xp_pauli.py @@ -0,0 +1,1658 @@ +# This code is part of Qiskit. +# +# (C) Copyright IBM 2017, 2020 +# +# 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. + +# Part of the QEC framework +# +# This code is based on the paper: "The XP Stabiliser Formalism: a +# Generalisation of the Pauli Stabiliser Formalism with Arbitrary Phases", Mark +# A. Webster, Benjamin J. Brown, and Stephen D. Bartlett. Quantum 6, 815 +# (2022). +"""Module for base XP pauli""" + +import numbers +from typing import List, Optional, Union +import warnings + +import numpy as np + +# Must be imported as follows to avoid circular import errors +import qiskit.quantum_info.operators.symplectic.clifford +from qiskit import QiskitError +from qiskit.circuit import QuantumCircuit +from qiskit.circuit.barrier import Barrier +from qiskit.circuit.delay import Delay +from qiskit.circuit.instruction import Instruction +from qiskit.quantum_info.operators.base_operator import BaseOperator +from qiskit.quantum_info.operators.mixins import AdjointMixin, MultiplyMixin +from qiskit_qec.utils import xp_pauli_rep + + +# pylint: disable=no-member +class BaseXPPauli(BaseOperator, AdjointMixin, MultiplyMixin): + r"""Symplectic representation of a list of N-qubit XP operators with phases + using numpy arrays for generalized symplectic matrices and phase vectors. + + Base class for XPPauli and XPPauliList. + """ + + # External string formats used when displaying Pauli's as strings + EXTERNAL_TENSOR_ENCODING = xp_pauli_rep.DEFAULT_EXTERNAL_TENSOR_ENCODING + EXTERNAL_PHASE_ENCODING = xp_pauli_rep.DEFAULT_EXTERNAL_PHASE_ENCODING + EXTERNAL_XP_PAULI_ENCODING = EXTERNAL_PHASE_ENCODING + EXTERNAL_TENSOR_ENCODING + + EXTERNAL_SYNTAX = xp_pauli_rep.PRODUCT_SYNTAX + EXTERNAL_QUBIT_ORDER = xp_pauli_rep.DEFAULT_QUBIT_ORDER + + PRINT_PHASE_ENCODING = None + + # pylint: disable=unused-argument + def __init__( + self, + matrix: Union[np.ndarray, None] = None, + phase_exp: Union[None, np.ndarray, np.integer] = None, + precision: int = None, + order: str = "xz", + ) -> None: + # TODO: Need to update this docstring once the representation to be + # used for XP operators has been decided. In addition, + # (-i)^phase_exp Z^z X^x and GF(2) need to be updated to XP operators. + """A BaseXPPauli object represents a list N-qubit XP Pauli operators with phases. + Numpy arrays are used to represent the generalized symplectic matrix represention of these + XP operators. The phases of the XP operators are stored encoded. The phases of the XP operators + operators are internally encoded in the '-iZX' Pauli encoding (See the xp_pauli_rep + module for more details). That is a XP operator is represented as generalized symplectic + vector V and a phase exponent phase_exp such that: + + (-i)^phase_exp Z^z X^x + + where V = [x, z] and phase_exp is a vector of Z_2N elements (0,1,2,...,2N-1). A list + of XP operators is represented as a generalized symplectic matrix S and a phase exponent + vector phase_exp such that the rows or S are the generalized symplectic vector representations + of the XP operators and the phase_exp vector store the phase exponent of each + associated XP Operator. + + Args: + matrix: Input GF(2) symplectic matrix + phase_exp (optional): Phase exponent vector for input matrix. A value of None will + result in an a complex coefficients of 1 for each XP operator. Defaults to None. + precision: Precision of XP operators. Must be an integer greater than or equal to 2. + order: Set to 'xz' or 'zx'. Defines which side the x and z parts of the input matrix + + Raises: + QiskitError: matrix and phase_exp sizes are not compatible, or if precision is less than 2. + + Examples: + >>> matrix = numpy.array([[1,1,0,0],[0,1,0,1]]) + >>> base_xp_pauli = BaseXPPauli(matrix) + + See also: + XPPauli, XPPauliList + """ + if not (isinstance(precision, int) and (precision > 1)): + raise QiskitError( + "Precision of XP operators must be an integer greater than or equal to 2." + ) + + if matrix is None or matrix.size == 0: + matrix = np.empty(shape=(0, 0), dtype=np.int64) + phase_exp = np.empty(shape=(0,), dtype=np.int64) + matrix = np.atleast_2d(matrix) + num_qubits = matrix.shape[1] >> 1 + + if order == "zx": + nmatrix = np.empty(shape=matrix.shape, dtype=matrix.dtype) + nmatrix[:, :num_qubits] = matrix[:, num_qubits:] + nmatrix[:, num_qubits:] = matrix[:, :num_qubits] + matrix = nmatrix + + self.matrix = matrix + self.precision = precision + self._num_xppaulis = self.matrix.shape[0] + if phase_exp is None: + self._phase_exp = np.zeros(shape=(self.matrix.shape[0],), dtype=np.int64) + else: + self._phase_exp = np.atleast_1d(phase_exp) + if self._phase_exp.shape[0] != self.matrix.shape[0]: + raise QiskitError("matrix and phase_exp sizes are not compatible") + + super().__init__(num_qubits=num_qubits) + + # --------------------------------------------------------------------- + # Properties + # --------------------------------------------------------------------- + + @property + def x(self): + """ + Returns the X component of symplectic representation as a 2d matrix. + + Note: The XPPauli class over writes this method to return + a 1d array instead of a 2d array. Use the self._x method + if a 2d array is needed as _x method is markeded as @final + + Examples: + >>> matrix = np.array([[1,0,0,0],[0,1,1,1]], dtype=np.int64) + >>> phase_exp = np.array([0,1]) + >>> precision = 4 + >>> base_xp_pauli = BaseXPPauli(matrix, phase_exp, precision) + >>> base_xp_pauli.x + array([[1, 0], + [0, 1]]) + + See Also: + _x, z, _z + """ + return self.matrix[:, : self.num_qubits] + + @x.setter + def x(self, val: np.ndarray): + """ + Sets the X component of symplectic representation + + Args: + val: integer matrix used to set the X component of the symplectic representation + + Examples: + >>> matrix = np.array([[1,0,0,0],[0,1,1,1]], dtype=np.int64) + >>> phase_exp = np.array([0,1]) + >>> precision = 4 + >>> base_xp_pauli = BaseXPPauli(matrix, phase_exp, precision) + >>> base_xp_pauli.x = np.array([[1,1],[0,0]], dtype=np.int64) + >>> base_xp_pauli.x + array([[1, 1], + [0, 0]]) + + See Also: + x, z, _z + """ + self.matrix[:, : self.num_qubits] = val + + # @final Add when python >= 3.8 + @property + def _x(self): # pylint: disable=invalid-name + """ + Returns the X component of symplectic representation as a 2d matrix. + + Note: The XPPauli class over writes this method to return + a 1d array instead of a 2d array. Use the self._x method + if a 2d array is needed as _x method is markeded as @final + + Examples: + >>> matrix = np.array([[1,0,0,0],[0,1,1,1]], dtype=np.int64) + >>> phase_exp = np.array([0,1]) + >>> precision = 4 + >>> base_xp_pauli = BaseXPPauli(matrix, phase_exp, precision) + >>> base_xp_pauli._x + array([[1, 0], + [0, 1]]) + + See Also: + x, z, _z + """ + return self.matrix[:, : self.num_qubits] + + # @final Add when python >= 3.8 + @_x.setter + def _x(self, val: np.ndarray): # pylint: disable=invalid-name + """ + Sets the X component of symplectic representation + + Args: + val: integer matrix used to set the X component of the symplectic representation + + Examples: + >>> matrix = np.array([[1,0,0,0],[0,1,1,1]], dtype=np.int64) + >>> phase_exp = np.array([0,1]) + >>> precision = 4 + >>> base_xp_pauli = BaseXPPauli(matrix, phase_exp, precision) + >>> base_xp_pauli._x = np.array([[1,1],[0,0]], dtype=np.int64) + >>> base_xp_pauli._x + array([[1, 1], + [0, 0]]) + + See Also: + _x, z, _z + """ + self.matrix[:, : self.num_qubits] = val + + @property + def z(self): + """ + Returns the Z component of symplectic representation as a 2d matrix. + + Note: The XPPauli class over writes this method to return + a 1d array instead of a 2d array. Use the self._z method + if a 2d array is needed as _z method is markeded as @final + + Examples: + >>> matrix = np.array([[1,0,0,0],[0,1,1,1]], dtype=np.int64) + >>> phase_exp = np.array([0,1]) + >>> precision = 4 + >>> base_xp_pauli = BaseXPPauli(matrix, phase_exp, precision) + >>> base_xp_pauli.z + array([[0, 0], + [1, 1]]) + + See Also: + _z, x, _x + """ + return self.matrix[:, self.num_qubits :] + + @z.setter + def z(self, val: np.ndarray): + """ + Sets the Z component of symplectic representation + + Args: + val: integer matrix used to set the Z component of the symplectic representation + + Examples: + >>> matrix = np.array([[1,0,0,0],[0,1,1,1]], dtype=np.int64) + >>> phase_exp = np.array([0,1]) + >>> precision = 4 + >>> base_xp_pauli = BaseXPPauli(matrix, phase_exp, precision) + >>> base_xp_pauli.z = np.array([[1,1],[0,0]], dtype=np.int64) + >>> base_xp_pauli.z + array([[1, 1], + [0, 0]]) + + See Also: + z, x, _x + """ + self.matrix[:, self.num_qubits :] = val + + # @final Add when python >= 3.8 + @property + def _z(self): # pylint: disable=invalid-name + """ + Returns the Z component of symplectic representation as a 2d matrix. + + Note: The XPPauli class over writes this method to return + a 1d array instead of a 2d array. Use the self._z method + if a 2d array is needed as _z method is markeded as @final + + Examples: + >>> matrix = np.array([[1,0,0,0],[0,1,1,1]], dtype=np.int64) + >>> phase_exp = np.array([0,1]) + >>> precision = 4 + >>> base_xp_pauli = BaseXPPauli(matrix, phase_exp, precision) + >>> base_xp_pauli._z + array([[0, 0], + [1, 1]]) + + See Also: + z, x, _x + """ + return self.matrix[:, self.num_qubits :] + + # @final Add when python >= 3.8 + @_z.setter + def _z(self, val: np.ndarray): # pylint: disable=invalid-name + """ + Sets the Z component of symplectic representation + + Args: + val: integer matrix used to set the Z component of the symplectic representation + + Examples: + >>> matrix = np.array([[1,0,0,0],[0,1,1,1]], dtype=np.int64) + >>> phase_exp = np.array([0,1]) + >>> precision = 4 + >>> base_xp_pauli = BaseXPPauli(matrix, phase_exp, precision) + >>> base_xp_pauli._z = np.array([[1,1],[0,0]], dtype=np.int64) + >>> base_xp_pauli._z + array([[1, 1], + [0, 0]]) + + See Also: + _z, x, _x + """ + self.matrix[:, self.num_qubits :] = val + + @property + def num_y(self): + """_summary_""" + pass + + @property + def tensor_encoding(self): + """_summary_""" + pass + + @classmethod + def set_tensor_encoding(cls, encoding: str = xp_pauli_rep.DEFAULT_EXTERNAL_TENSOR_ENCODING): + """_summary_""" + pass + + @property + def phase_encoding(self): + """_summary_""" + pass + + @classmethod + def set_phase_encoding(cls, encoding: str = xp_pauli_rep.DEFAULT_EXTERNAL_PHASE_ENCODING): + """_summary_""" + pass + + @property + def pauli_encoding(self): + """_summary_""" + pass + + @classmethod + def set_pauli_encoding(cls, encoding: str = xp_pauli_rep.DEFAULT_EXTERNAL_XP_PAULI_ENCODING): + """_summary_""" + pass + + @property + def syntax(self): + """_summary_""" + pass + + @classmethod + def set_syntax(cls, syntax_code: Optional[int] = None, syntax_str: Optional[str] = "Product"): + """_summary_""" + pass + + @property + def print_phase_encoding(self): + """_summary_""" + pass + + @classmethod + def set_print_phase_encoding(cls, phase_encoding: Optional[str] = None): + """_summary_""" + pass + + @property + def qubit_order(self): + """_summary_""" + pass + + @classmethod + def set_qubit_order(cls, qubit_order: Optional[str] = None): + """_summary_""" + pass + + # --------------------------------------------------------------------- + # Magic Methods + # --------------------------------------------------------------------- + + def __imul__(self, other: "BaseXPPauli") -> "BaseXPPauli": + """_summary_""" + pass + + def __neg__(self) -> "BaseXPPauli": + """_summary_""" + pass + + # --------------------------------------------------------------------- + + def copy(self) -> "BaseXPPauli": + """_summary_""" + pass + + # --------------------------------------------------------------------- + # BaseOperator methods + # --------------------------------------------------------------------- + # Needed by GroupMixin class from BaseOperator class + + # pylint: disable=arguments-differ + def compose( + self, + other: "BaseXPPauli", + qargs: Optional[list] = None, + front: bool = False, + inplace: bool = False, + ) -> "BaseXPPauli": + r"""Return the composition of XPPaulis lists + + To be consistent with other compose functions in Qiskit, composition is defined via + left multiplication. That is + + A.compose(B) = B.A = B.dot(A) = A.compose(B, front=False) + + where . is the Pauli group multiplication and so B is applied after A. Likewise + + A.compose(B, front=True) = A.B = A.dot(B) + + That is B is applied first or at the front. + + This compose is: + + [A_1,A_2,...,A_k].compose([B_1,B_2,...,B_k]) = [A_1.compose(B_1),...,A_k.compose(B_k)] + + or + + [A].compose([B_1,B_2,...,B_k])) = [A.compose(B_1),...,A.compose(B_k)] + + Note: + This method does compose coordinate wise (which is different from the PauliTable compose + which should be corrected at some point). + This method is adapted from method XPMul from XPFpackage: + https://github.com/m-webster/XPFpackage, originally developed by + Mark Webster. The original code is licensed under the GNU General + Public License v3.0 and Mark Webster has given permission to use + the code under the Apache License v2.0. + + Args: + other: BaseXPPauli + front (bool): (default: False) + qargs (list or None): Optional, qubits to apply compose on + on (default: None->All). + inplace (bool): If True update in-place (default: False). + + Returns: + BaseXPPauli : Compositon of self and other + + Raises: + QiskitError: if number of qubits of other does not match qargs, or + if precision of other does not match precision of self. + + Examples: + >>> a = BaseXPPauli(matrix=np.array([0, 1, 0, 0, 2, 0], dtype=np.int64), + ... phase_exp=6, precision=4) + >>> b = BaseXPPauli(matrix=np.array([1, 1, 1, 3, 3, 0], dtype=np.int64), + ... phase_exp=2, precision=4) + >>> value = a.compose(b) + >>> value.matrix + array([[1, 0, 1, 3, 3, 0]], dtype=np.int64) + >>> value._phase_exp + array([6]) + + See also: + _compose + """ + + # Validation + if qargs is None and other.num_qubits != self.num_qubits: + raise QiskitError(f"other {type(self).__name__} must be on the same number of qubits.") + + if qargs and other.num_qubits != len(qargs): + raise QiskitError( + f"Number of qubits of the other {type(self).__name__} does not match qargs." + ) + + if other._num_xppaulis not in [1, self._num_xppaulis]: + raise QiskitError( + "Incompatible BaseXPPaulis. Second list must " + "either have 1 or the same number of XPPaulis." + ) + + if self.precision != other.precision: + raise QiskitError( + "Precision of the two BaseXPPaulis to be multiplied must be the same." + ) + + return self._compose(self, other, qargs=qargs, front=front, inplace=inplace) + + @staticmethod + def _compose( + a: "BaseXPPauli", + b: "BaseXPPauli", + qargs: Optional[list] = None, + front: bool = False, + inplace: bool = False, + ) -> "BaseXPPauli": + """Returns the composition of two BaseXPPauli objects. + + Note: + This method is adapted from method XPMul from XPFpackage: + https://github.com/m-webster/XPFpackage, originally developed by + Mark Webster. The original code is licensed under the GNU General + Public License v3.0 and Mark Webster has given permission to use + the code under the Apache License v2.0. + + Args: + a : BaseXPPauli object + b : BaseXPPauli object + qargs (Optional[list], optional): _description_. Defaults to None. + front (bool, optional): _description_. Defaults to False. + inplace (bool, optional): _description_. Defaults to False. + + Returns: + BaseXPPauli: _description_ + + See also: + antisymmetric_op, unique_vector_rep, _unique_vector_rep + """ + + if qargs is not None: + qargs = list(qargs) + [item + a.num_qubits for item in qargs] + amat = a.matrix[:, qargs] + else: + amat = a.matrix + bmat = b.matrix + + # Calculate the sum of generalized symplectic matrix for the composition, excluding D + x = amat[:, : a.num_qubits] + bmat[:, : b.num_qubits] + z = amat[:, a.num_qubits :] + bmat[:, b.num_qubits :] + mat = np.concatenate((x, z), axis=-1) + + # Calculate the phase of the composition, excluding D + phase_exp = a._phase_exp + b._phase_exp + # Calculate the antisymmetric operator, i.e. D + if front: + dinput = 2 * np.multiply(b.x, a.z) + d = b.antisymmetric_op(dinput) + else: + dinput = 2 * np.multiply(a.x, b.z) + d = a.antisymmetric_op(dinput) + + if qargs is None: + if not inplace: + result_x = x + d.x + result_z = z + d.z + result_phase_exp = phase_exp + d._phase_exp + result_mat = np.concatenate((result_x, result_z), axis=-1) + return BaseXPPauli( + matrix=result_mat, phase_exp=result_phase_exp, precision=a.precision + )._unique_vector_rep() + # Inplace update + a.x = x + d.x + a.z = z + d.z + a._phase_exp = phase_exp + d._phase_exp + return a._unique_vector_rep() + + # Qargs update + ret = a if inplace else a.copy() + ret.matrix[:, qargs] = mat + ret._phase_exp = phase_exp + d._phase_exp + ret = ret._unique_vector_rep() + return ret + + # --------------------------------------------------------------------- + + def tensor(self, other): + """_summary_""" + pass + + @staticmethod + def _tensor(a: "BaseXPPauli", b: "BaseXPPauli"): + """_summary_""" + pass + + # --------------------------------------------------------------------- + + def expand(self, other): + """_summary_""" + pass + + # --------------------------------------------------------------------- + # Needed by MultiplyMixin class + # --------------------------------------------------------------------- + + def _multiply(self, phase, roundit=True) -> "BaseXPPauli": # pylint: disable=arguments-renamed + """_summary_""" + pass + + # Needed by AdjointMixin class + + def transpose(self, inplace: bool = False) -> "BaseXPPauli": + """_summary_""" + pass + + def commutes(self, other: "BaseXPPauli", qargs: List = None) -> np.ndarray: + """_summary_""" + pass + + def _commutes(self, other: "BaseXPPauli", qargs: List = None) -> np.ndarray: + """_summary_""" + pass + + # --------------------------------------------------------------------- + # Extra all_commutes method + # --------------------------------------------------------------------- + + def all_commutes(self, other: "BaseXPPauli") -> np.ndarray: + """_summary_""" + pass + + # --------------------------------------------------------------------- + # Evolve Methods + # --------------------------------------------------------------------- + + def evolve(self, other: "BaseXPPauli", qargs: Union[None, List, int] = None, frame: str = "h"): + """_summary_""" + pass + + def _evolve_clifford( + self, + other: qiskit.quantum_info.operators.symplectic.clifford.Clifford, + qargs: Union[None, List, int] = None, + frame: str = "h", + ): + """_summary_""" + pass + + # --------------------------------------------------------------------- + # Helper Metods + # --------------------------------------------------------------------- + + def _eq(self, other): + """_summary_""" + pass + + # --------------------------------------------------------------------- + # Class methods for creating labels -> Uses xp_pauli_rep suite of methods + # --------------------------------------------------------------------- + + def to_label( + self, + output_pauli_encoding: Optional[str] = None, + no_phase: bool = False, + return_phase: bool = False, + syntax: Optional[int] = None, + qubit_order: Optional[str] = None, + index_start: int = 0, + squeeze: bool = True, + index_str: str = "", + ) -> Union[str, List[str]]: + """_summary_""" + pass + + # --------------------------------------------------------------------- + # The methods below should be deprecated eventually as they simply refer + # to the newer versions that are now elsewhere. + # --------------------------------------------------------------------- + + def _count_y(self) -> np.ndarray: + """_summary_""" + pass + + @staticmethod + def _stack(array: np.ndarray, size: int, vertical: bool = True) -> np.ndarray: + """_summary_""" + pass + + @staticmethod + def _phase_from_complex(coeff: numbers.Complex) -> np.ndarray: + """_summary_""" + pass + + # --------------------------------------------------------------------- + # Apply Clifford to BaseXPPauli + # --------------------------------------------------------------------- + + def _append_circuit( + self, + circuit: Union[Barrier, Delay, QuantumCircuit, Instruction], + qargs: Optional[List] = None, + ) -> "BaseXPPauli": + """_summary_""" + pass + + # --------------------------------------------------------------------- + # BaseXPPauli methods for XP arithmetic + # --------------------------------------------------------------------- + + def unique_vector_rep(self) -> "BaseXPPauli": + """Convert the BaseXPPauli operator into unique vector form, ie + phase_exp in Z modulo 2*precision, x in Z_2, z in Z modulo + precision. + + Note: + This method is adapted from method XPRound from XPFpackage: + https://github.com/m-webster/XPFpackage, originally developed by + Mark Webster. The original code is licensed under the GNU General + Public License v3.0 and Mark Webster has given permission to use + the code under the Apache License v2.0. + + Returns: + BaseXPPauli: Unique vector representation of self + + Examples: + >>> a = BaseXPPauli(matrix=np.array([0, 3, 1, 6, 4, 3], dtype=np.int64), + ... phase_exp=11, precision=4) + >>> a = a.unique_vector_rep() + >>> a.matrix + np.array([[0, 1, 1, 2, 0, 3]], dtype=np.int64) + >>> a._phase_exp + array([3]) + + See also: + _unique_vector_rep + """ + return self._unique_vector_rep() + + def _unique_vector_rep(self) -> "BaseXPPauli": + """Convert the BaseXPPauli operator into unique vector form, ie + phase_exp in Z modulo 2*precision, x in Z_2, z in Z modulo + precision. + + Note: + This method is adapted from method XPRound from XPFpackage: + https://github.com/m-webster/XPFpackage, originally developed by + Mark Webster. The original code is licensed under the GNU General + Public License v3.0 and Mark Webster has given permission to use + the code under the Apache License v2.0. + + Returns: + BaseXPPauli: Unique vector representation of self + """ + matrix = np.empty(shape=np.shape(self.matrix), dtype=np.int64) + + phase_exp = np.mod(self._phase_exp, 2 * self.precision) + matrix[:, : self.num_qubits] = np.mod(self.x, 2) + matrix[:, self.num_qubits :] = np.mod(self.z, np.expand_dims(self.precision, axis=-1)) + + return BaseXPPauli(matrix, phase_exp, self.precision) + + def rescale_precision(self, new_precision: int, inplace: bool = False) -> "BaseXPPauli": + """Rescale the generalized symplectic vector components of BaseXPPauli + operator to the new precision. + + Note: + This method is adapted from method XPSetN from XPFpackage: + https://github.com/m-webster/XPFpackage, originally developed by + Mark Webster. The original code is licensed under the GNU General + Public License v3.0 and Mark Webster has given permission to use + the code under the Apache License v2.0. + + Args: + new_precision: The target precision in which BaseXPPauli is to be expressed + inplace: If True, rescale BaseXPPauli in place, otherwise return a new + BaseXPPauli object. Defaults to False + + Returns: + BaseXPPauli: Resultant of rescaling the precision of BaseXPPauli + + Raises: + QiskitError: If it is not possible to express BaseXPPauli in new_precision + + Examples: + >>> a = BaseXPPauli( + ... matrix=np.array([1, 1, 1, 0, 0, 0, 0, 0, 0, 4, 0, 0, 0, 0], dtype=np.int64), + ... phase_exp=12, precision=8) + >>> a = a.rescale_precision(new_precision=2) + >>> a.matrix + array([[1, 1, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0]], dtype=int64) + >>> a._phase_exp + array([3]) + + See also: + _rescale_precision + """ + unique_xp_op = self.unique_vector_rep() + old_precision = self.precision + if new_precision < old_precision: + scale_factor = old_precision // new_precision + if (new_precision > old_precision) and (new_precision % old_precision > 0): + raise QiskitError("XP Operator cannot be expressed in new_precision.") + if (new_precision < old_precision) and ( + (old_precision % new_precision > 0) + or (np.sum(np.mod(unique_xp_op._phase_exp, scale_factor)) > 0) + or (np.sum(np.mod(unique_xp_op.z, scale_factor)) > 0) + ): + raise QiskitError("XP Operator cannot be expressed in new_precision.") + + return self._rescale_precision(new_precision, inplace) + + def _rescale_precision(self, new_precision: int, inplace: bool = False) -> "BaseXPPauli": + """Rescale the generalized symplectic vector components + of BaseXPPauli operator to the new precision. Returns None if the + rescaling is not possible, else returns the rescaled BaseXPPauli object. + + Note: + This method is adapted from method XPSetN from XPFpackage: + https://github.com/m-webster/XPFpackage, originally developed by + Mark Webster. The original code is licensed under the GNU General + Public License v3.0 and Mark Webster has given permission to use + the code under the Apache License v2.0. + + Args: + new_precision: The target precision in which BaseXPPauli is to be expressed + inplace: If True, rescale BaseXPPauli in place, else return a new + BaseXPPauli. Defaults to False + + Returns: + BaseXPPauli: Resultant of rescaling the precision of BaseXPPauli + + Warning: + If precision rescaling is not possible, a warning is raised and the unique vector + representation of the original XP operator is returned. + + See also: + unique_vector_rep + """ + unique_xp_op = self.unique_vector_rep() + old_precision = unique_xp_op.precision + matrix = np.empty(shape=np.shape(unique_xp_op.matrix), dtype=np.int64) + phase_exp = np.empty(shape=np.shape(unique_xp_op._phase_exp)) + + if new_precision > old_precision: + if np.mod(new_precision, old_precision) > 0: + warnings.warn( + "Precision rescaling is not possible. Returning the unique vector representation of the original XP operator." + ) + return unique_xp_op + scale_factor = new_precision // old_precision + phase_exp = scale_factor * unique_xp_op._phase_exp + matrix[:, unique_xp_op.num_qubits :] = scale_factor * np.atleast_2d(unique_xp_op.z) + + else: + scale_factor = old_precision // new_precision + if ( + (old_precision % new_precision > 0) + or (np.sum(np.mod(unique_xp_op._phase_exp, scale_factor)) > 0) + or (np.sum(np.mod(unique_xp_op.z, scale_factor)) > 0) + ): + warnings.warn( + "Precision rescaling is not possible. Returning the unique vector representation of the original XP operator." + ) + return unique_xp_op + phase_exp = unique_xp_op._phase_exp // scale_factor + matrix[:, unique_xp_op.num_qubits :] = np.atleast_2d(unique_xp_op.z) // scale_factor + + matrix[:, 0 : unique_xp_op.num_qubits] = unique_xp_op.x + + if not inplace: + return BaseXPPauli(matrix, phase_exp, new_precision) + self.matrix = matrix + self._phase_exp = phase_exp + self.precision = new_precision + return self + + def weight(self) -> Union[int, np.ndarray]: + """Return the weight, i.e. count of qubits where either z or x component is nonzero. + + Note: + This method is adapted from method XPDistance from XPFpackage: + https://github.com/m-webster/XPFpackage, originally developed by + Mark Webster. The original code is licensed under the GNU General + Public License v3.0 and Mark Webster has given permission to use + the code under the Apache License v2.0. + + Returns: + Union[int, np.ndarray]: Weight of BaseXPPauli + + Examples: + >>> a = BaseXPPauli( + ... matrix=np.array([1, 1, 1, 0, 0, 1, 0, 0, 3, 4, 0, 0, 0, 1], dtype=np.int64), + ... phase_exp = 12, precision = 8) + >>> a.weight() + array([5]) + + See also: + _weight + """ + return self._weight() + + def _weight(self) -> Union[int, np.ndarray]: + """Return the weight, i.e. count of qubits where either z or x component is nonzero. + + Note: + This method is adapted from method XPDistance from XPFpackage: + https://github.com/m-webster/XPFpackage, originally developed by + Mark Webster. The original code is licensed under the GNU General + Public License v3.0 and Mark Webster has given permission to use + the code under the Apache License v2.0. + + Returns: + Union[int, np.ndarray]: Weight of BaseXPPauli + """ + return np.sum(np.logical_or(self.x, self.z), axis=-1) + + def is_diagonal(self) -> np.ndarray: + """Return True if the XP operator is diagonal. + + Note: + This method is adapted from method XPisDiag from XPFpackage: + https://github.com/m-webster/XPFpackage, originally developed by + Mark Webster. The original code is licensed under the GNU General + Public License v3.0 and Mark Webster has given permission to use + the code under the Apache License v2.0. + + Returns: + np.ndarray: True, if BaseXPPauli is diagonal + + Examples: + >>> a = BaseXPPauli( + ... matrix=np.array([0, 0, 0, 0, 0, 0, 0, 0, 1, 3, 3, 3, 3, 3], dtype=np.int64), + ... phase_exp=0, precision=8) + >>> a.is_diagonal() + array([True]) + + >>> b = BaseXPPauli( + ... matrix=np.array([0, 1, 0, 1, 0, 0, 0, 0, 1, 3, 3, 3, 3, 3], dtype=np.int64), + ... phase_exp=12, precision=8) + >>> b.is_diagonal() + array([False]) + + See also: + _is_diagonal + """ + return self._is_diagonal() + + def _is_diagonal(self) -> np.ndarray: + """Return True if the XP operator is diagonal. + + Note: + This method is adapted from method XPisDiag from XPFpackage: + https://github.com/m-webster/XPFpackage, originally developed by + Mark Webster. The original code is licensed under the GNU General + Public License v3.0 and Mark Webster has given permission to use + the code under the Apache License v2.0. + + Returns: + np.ndarray: True, if BaseXPPauli is diagonal + """ + return np.where(np.sum(self.x, axis=-1) == 0, True, False) + + def antisymmetric_op(self, int_vec: np.ndarray) -> "BaseXPPauli": + """Return the antisymmetric operator corresponding to an integer vector, + with precision specified by BaseXPPauli. + + Note: + This method is adapted from method XPD from XPFpackage: + https://github.com/m-webster/XPFpackage, originally developed by + Mark Webster. The original code is licensed under the GNU General + Public License v3.0 and Mark Webster has given permission to use + the code under the Apache License v2.0. + + Args: + int_vec: An integer vector + + Returns: + BaseXPPauli: The antisymmetric operator + + Raises: + QiskitError: Input vector must be an integer array + + Examples: + >>> a = BaseXPPauli( + ... matrix=np.array([0, 0, 0, 0, 0, 0, 0, 0, 1, 2, 3, 3, 3, 3], dtype=np.int64), + ... phase_exp=0, precision=8) + >>> value = a.antisymmetric_op(a.z) + >>> value.matrix + array([0, 0, 0, 0, 0, 0, 0, 0, -1, -2, -3, -3, -3, -3], dtype=np.int64) + >>> value._phase_exp + array([15]) + + See also: + _antisymmetric_op + """ + # Check if int_vec is a valid integer vector + int_vec = np.atleast_2d(int_vec) + if int_vec.dtype != np.int64: + raise TypeError("Input vector must be an integer array.") + + return self._antisymmetric_op(int_vec, self.precision) + + @staticmethod + def _antisymmetric_op(int_vec: np.ndarray, precision: int) -> "BaseXPPauli": + """Return the antisymmetric operator of specified precision + corresponding to an integer vector. + + Note: + This method is adapted from method XPD from XPFpackage: + https://github.com/m-webster/XPFpackage, originally developed by + Mark Webster. The original code is licensed under the GNU General + Public License v3.0 and Mark Webster has given permission to use + the code under the Apache License v2.0. + + Args: + int_vec: An integer vector + precision: Precision of the antisymmetric operator + + Returns: + BaseXPPauli: The antisymmetric operator of specified precision + """ + phase_exp = np.sum(int_vec, axis=-1, dtype=np.int64) + x = np.zeros(np.shape(int_vec), dtype=np.int64) + z = -int_vec + matrix = np.concatenate((x, z), axis=-1) + + return BaseXPPauli(matrix=matrix, phase_exp=phase_exp, precision=precision) + + def inverse(self) -> "BaseXPPauli": + """Return the inverse of the XP operator. + + Note: + This method is adapted from method XPInverse from XPFpackage: + https://github.com/m-webster/XPFpackage, originally developed by + Mark Webster. The original code is licensed under the GNU General + Public License v3.0 and Mark Webster has given permission to use + the code under the Apache License v2.0. + + Returns: + BaseXPPauli: Inverse of BaseXPPauli + + Examples: + >>> a = BaseXPPauli( + ... matrix=np.array([0, 0, 0, 1, 0, 1, 1, 5, 5, 6, 1, 1, 4, 0], dtype=np.int64), + ... phase_exp=1, precision=8) + >>> value = a.inverse() + >>> value.matrix + array([0, 0, 0, 1, 0, 1, 1, 3, 3, 2, 1, 7, 4, 0], dtype=np.int64) + >>> value._phase_exp + array([5]) + + See also: + _inverse + """ + return self._inverse() + + def _inverse(self) -> "BaseXPPauli": + """Return the inverse of the XP operator. + + Note: + This method is adapted from method XPInverse from XPFpackage: + https://github.com/m-webster/XPFpackage, originally developed by + Mark Webster. The original code is licensed under the GNU General + Public License v3.0 and Mark Webster has given permission to use + the code under the Apache License v2.0. + + Returns: + BaseXPPauli: Inverse of BaseXPPauli + + See also: + _antisymmetric_op, _unique_vec_rep + """ + phase_exp = -self._phase_exp + matrix = np.concatenate((self.x, -self.z), axis=-1) + first = BaseXPPauli(matrix=matrix, phase_exp=phase_exp, precision=self.precision) + + dinput = -2 * np.multiply(self.x, self.z) + second = self._antisymmetric_op(dinput, self.precision) + + product = BaseXPPauli( + matrix=first.matrix + second.matrix, + phase_exp=first._phase_exp + second._phase_exp, + precision=self.precision, + ) + + return product._unique_vector_rep() + + def power(self, n: Union[int, list, np.ndarray]) -> "BaseXPPauli": + """Return the XP operator of specified precision raised to the power n. + + For a list of XP operators, power is performed element-wise: + + [A_1, ..., A_k].power([n_1, ..., n_k]) = [A_1.power(n_1), ..., A_k.power(n_k)]. + + Note: + This method is adapted from method XPPower from XPFpackage: + https://github.com/m-webster/XPFpackage, originally developed by + Mark Webster. The original code is licensed under the GNU General + Public License v3.0 and Mark Webster has given permission to use + the code under the Apache License v2.0. + + Args: + n: The power to which BaseXPPauli is to be raised + + Returns: + BaseXPPauli: BaseXPPauli raised to the power n + + Raises: + QiskitError: The number of powers in the array n must match the + number of XP operators + + Examples: + >>> a = BaseXPPauli( + ... matrix=np.array([1, 0, 1, 1, 5, 3, 5, 4], dtype=np.int64), + ... phase_exp=4, precision=6) + >>> value = a.power(n=5) + >>> value.matrix + array([1, 0, 1, 1, 5, 3, 5, 4], dtype=np.int64) + >>> value._phase_exp + array([4]) + + See also: + _power + """ + if isinstance(n, list): + n = np.array(n, dtype=np.int64) + if isinstance(n, np.ndarray) and n.shape != self.matrix.shape[:-1]: + raise QiskitError( + "The number of powers in the array n must match the number of XP operators." + ) + return self._power(n) + + def _power(self, n: Union[int, np.ndarray]) -> "BaseXPPauli": + """Return the XP operator of specified precision raised to the power n. + + Note: + This method is adapted from method XPPower from XPFpackage: + https://github.com/m-webster/XPFpackage, originally developed by + Mark Webster. The original code is licensed under the GNU General + Public License v3.0 and Mark Webster has given permission to use + the code under the Apache License v2.0. + + Args: + n: The power to which BaseXPPauli is to be raised + + Returns: + BaseXPPauli: BaseXPPauli raised to the power n + + See also: + _antisymmetric_op, _unique_vector_rep + """ + if isinstance(n, int): + n = np.array([n] * self.matrix.shape[0], dtype=np.int64) + + a = np.mod(n, 2) + + x = self.x * a[:, None] + z = self.z * n[:, None] + phase_exp = np.multiply(self._phase_exp, n) + matrix = np.concatenate((x, z), axis=-1) + first = BaseXPPauli(matrix=matrix, phase_exp=phase_exp, precision=self.precision) + + dinput = np.multiply(self.x, self.z) * (n - a)[:, None] + second = self._antisymmetric_op(dinput, self.precision) + + product = BaseXPPauli( + matrix=first.matrix + second.matrix, + phase_exp=first._phase_exp + second._phase_exp, + precision=self.precision, + ) + + return product._unique_vector_rep() + + def conjugate( + self, other: "BaseXPPauli", front: bool = True, inplace: bool = False + ) -> "BaseXPPauli": + """Return the conjugation of two BaseXPPauli operators. + + For single XP operators, this means + + A.conjugate(B, front=True) = A . B . A^{-1}, + + where . is the XP Pauli multiplication and A^{-1} is the inverse of A. + + Likewise, + + A.conjugate(B, front=False) = B . A . B^{-1}. + + For a list of XP operators, conjugation is performed element-wise: + + [A_1, ..., A_k].conjugate([B_1, ..., B_k]) = [A_1.conjugate(B_1), ..., A_k.conjugate(B_k)]. + + TODO: This method currently only supports conjugation of two XP operator + lists of the same length. + + Note: + This method is adapted from method XPConjugate from XPFpackage: + https://github.com/m-webster/XPFpackage, originally developed by + Mark Webster. The original code is licensed under the GNU General + Public License v3.0 and Mark Webster has given permission to use + the code under the Apache License v2.0. + + Args: + other: BaseXPPauli object + front (bool, optional): Whether to conjugate in front (True) or + behind (False), defaults to True + inplace (bool, optional): Whether to perform the conjugation in + place (True) or to return a new BaseXPPauli (False), defaults to + False + + Returns: + BaseXPPauli: Conjugated XP operator + + Raises: + QiskitError: Other BaseXPPauli must be on the same number of qubits + QiskitError: Incompatible BaseXPPaulis. Second list must either have + 1 or the same number of XPPaulis + QiskitError: Precision of the two BaseXPPaulis to be conjugated must + be the same + + Examples: + >>> a = BaseXPPauli(matrix=np.array([1, 0, 1, 1, 5, 3, 5, 4], dtype=np.int64), + ... phase_exp=4, precision=6) + >>> b = BaseXPPauli(matrix=np.array([1, 0, 0, 1, 4, 1, 0, 1], dtype=np.int64), + ... phase_exp=11, precision=6) + >>> value = a.conjugate(b) + >>> value.matrix + array([[1, 0, 0, 1, 0, 1, 0, 1]], dtype=np.int64) + >>> value._phase_exp + array([3]) + + See also: + _conjugate + """ + # Validation + if other.num_qubits != self.num_qubits: + raise QiskitError(f"Other {type(self).__name__} must be on the same number of qubits.") + + if other._num_xppaulis not in [1, self._num_xppaulis]: + raise QiskitError( + "Incompatible BaseXPPaulis. Second list must " + "either have 1 or the same number of XPPaulis." + ) + + if self.precision != other.precision: + raise QiskitError( + "Precision of the two BaseXPPaulis to be conjugated must be the same." + ) + + return self._conjugate(self, other, front=front, inplace=inplace) + + @staticmethod + def _conjugate( + a: "BaseXPPauli", b: "BaseXPPauli", front: bool = True, inplace: bool = False + ) -> "BaseXPPauli": + """Return the conjugation of two BaseXPPauli operators. + + Note: + This method is adapted from method XPConjugate from XPFpackage: + https://github.com/m-webster/XPFpackage, originally developed by + Mark Webster. The original code is licensed under the GNU General + Public License v3.0 and Mark Webster has given permission to use + the code under the Apache License v2.0. + + Args: + a: BaseXPPauli object + b: BaseXPPauli object + front (bool, optional): Whether to conjugate in front (True) or + behind (False), defaults to True + inplace (bool, optional): Whether to perform the conjugation in + place (True) or to return a new BaseXPPauli (False), defaults to + False + + Returns: + BaseXPPauli: Conjugation of XP operators a and b + + See also: + antisymmetric_op, _unique_vector_rep + """ + if front: + dinput = ( + 2 * np.multiply(a.x, b.z) + + 2 * np.multiply(a.z, b.x) + - 4 * np.multiply(np.multiply(a.x, b.x), a.z) + ) + d = b.antisymmetric_op(dinput) + product = BaseXPPauli( + matrix=b.matrix + d.matrix, + phase_exp=b._phase_exp + d._phase_exp, + precision=b.precision, + ) + else: + dinput = ( + 2 * np.multiply(b.x, a.z) + + 2 * np.multiply(b.z, a.x) + - 4 * np.multiply(np.multiply(b.x, a.x), b.z) + ) + d = a.antisymmetric_op(dinput) + product = BaseXPPauli( + matrix=a.matrix + d.matrix, + phase_exp=a._phase_exp + d._phase_exp, + precision=a.precision, + ) + if not inplace: + return product._unique_vector_rep() + else: + a.matrix = product.matrix + a._phase_exp = product._phase_exp + return a + + def commutator( + self, other: "BaseXPPauli", front: bool = True, inplace: bool = False + ) -> "BaseXPPauli": + """Return the commutator of two BaseXPPauli operators. + + For single XP operators, this means + + A.commutator(B, front=True) = [A, B] = A . B . A^{-1} . B^{-1}, + + where . is the XP Pauli multiplication and A^{-1} is the inverse of A. + + Likewise, + + A.commutator(B, front=False) = [B, A] = B . A . B^{-1}. A^{-1}. + + For a list of XP operators, commutator is computed element-wise: + + [A_1, ..., A_k].commutator([B_1, ..., B_k]) = [A_1.commutator(B_1), ..., A_k.commutator(B_k)]. + + TODO: This method currently only supports commutator of two XP operator + lists of the same length. + + Note: + This method is adapted from method XPCommutator from XPFpackage: + https://github.com/m-webster/XPFpackage, originally developed by + Mark Webster. The original code is licensed under the GNU General + Public License v3.0 and Mark Webster has given permission to use + the code under the Apache License v2.0. + + Args: + other: BaseXPPauli object + front (bool, optional): Whether self is the first element in the + commutator (True) or second (False), defaults to True + inplace (bool, optional): Whether to compute the commutator in + place (True) or to return a new BaseXPPauli (False), defaults to + False + + Returns: + BaseXPPauli: Commutator of XP operators + + Raises: + QiskitError: Other BaseXPPauli must be on the same number of qubits + QiskitError: Incompatible BaseXPPaulis. Second list must either have + 1 or the same number of XPPaulis + QiskitError: Precision of the two BaseXPPaulis in a commutator must + be the same + + Examples: + >>> a = BaseXPPauli(matrix=np.array([1, 0, 1, 1, 5, 3, 5, 4], dtype=np.int64), + ... phase_exp=4, precision=6) + >>> b = BaseXPPauli(matrix=np.array([1, 0, 0, 1, 4, 1, 0, 1], dtype=np.int64), + ... phase_exp=11, precision=6) + >>> value = a.commutator(b) + >>> value.matrix + array([[0, 0, 0, 0, 4, 0, 0, 0]], dtype=np.int64) + >>> value._phase_exp + array([8]) + + See also: + _commutator + """ + # Validation + if other.num_qubits != self.num_qubits: + raise QiskitError(f"Other {type(self).__name__} must be on the same number of qubits.") + + if other._num_xppaulis not in [1, self._num_xppaulis]: + raise QiskitError( + "Incompatible BaseXPPaulis. Second list must " + "either have 1 or the same number of XPPaulis." + ) + + if self.precision != other.precision: + raise QiskitError("Precision of the two BaseXPPaulis in a commutator must be the same.") + + return self._commutator(self, other, front=front, inplace=inplace) + + @staticmethod + def _commutator( + a: "BaseXPPauli", b: "BaseXPPauli", front: bool = True, inplace: bool = False + ) -> "BaseXPPauli": + """Return the commutator of two BaseXPPauli operators. + + Note: + This method is adapted from method XPCommutator from XPFpackage: + https://github.com/m-webster/XPFpackage, originally developed by + Mark Webster. The original code is licensed under the GNU General + Public License v3.0 and Mark Webster has given permission to use + the code under the Apache License v2.0. + + Args: + a: BaseXPPauli object + b: BaseXPPauli object + front (bool, optional): Whether self is the first element in the + commutator (True) or second (False), defaults to True + inplace (bool, optional): Whether to compute the commutator in + place (True) or to return a new BaseXPPauli (False), defaults to + False + + Returns: + BaseXPPauli: Commutator of XP operators a and b + + See also: + antisymmetric_op, _unique_vector_rep + """ + if front: + dinput = ( + 2 * np.multiply(a.x, b.z) + - 2 * np.multiply(a.z, b.x) + + 4 * np.multiply(np.multiply(a.x, b.x), a.z) + - 4 * np.multiply(np.multiply(a.x, b.x), b.z) + ) + result = a.antisymmetric_op(dinput) + else: + dinput = ( + 2 * np.multiply(b.x, a.z) + - 2 * np.multiply(b.z, a.x) + + 4 * np.multiply(np.multiply(b.x, a.x), b.z) + - 4 * np.multiply(np.multiply(b.x, a.x), a.z) + ) + result = b.antisymmetric_op(dinput) + if not inplace: + return result._unique_vector_rep() + else: + a.matrix = result.matrix + a._phase_exp = result._phase_exp + return a + + def degree(self) -> np.ndarray: + """Return the degree of the XP operator. + + Note: + This method is adapted from method XPDegree from XPFpackage: + https://github.com/m-webster/XPFpackage, originally developed by + Mark Webster. The original code is licensed under the GNU General + Public License v3.0 and Mark Webster has given permission to use + the code under the Apache License v2.0. + + Returns: + np.ndarray: Degree of BaseXPPauli + + Examples: + >>> a = BaseXPPauli(matrix=np.array([0, 0, 0, 2, 1, 0], dtype=np.int64), + ... phase_exp=2, precision=4) + >>> a.degree() + array([4]) + + See also: + _degree + """ + return self._degree() + + def _degree(self) -> np.ndarray: + """Return the degree of the XP operator. + + Note: + This method is adapted from method XPDegree from XPFpackage: + https://github.com/m-webster/XPFpackage, originally developed by + Mark Webster. The original code is licensed under the GNU General + Public License v3.0 and Mark Webster has given permission to use + the code under the Apache License v2.0. + + Returns: + np.ndarray: Degree of BaseXPPauli + + See also: + is_diagonal + """ + + gcd = np.gcd(self.z, self.precision) + precision_by_gcd = np.atleast_2d(np.floor_divide(self.precision, gcd)) + lcm = np.atleast_2d(precision_by_gcd)[:, 0] + for i, val in enumerate(precision_by_gcd): + for j in val: + lcm[i] = np.lcm(lcm[i], j) + + square = self.compose(self) + if not isinstance(square, type(self)): + square = type(self)(square) + gcd_square = np.gcd(square.z, square.precision) + precision_by_gcd_square = np.atleast_2d(np.floor_divide(square.precision, gcd_square)) + lcm_square = np.atleast_2d(precision_by_gcd_square)[:, 0] + for i, val in enumerate(precision_by_gcd_square): + for j in val: + lcm_square[i] = np.lcm(lcm_square[i], j) + + lcm_square = 2 * lcm_square + + # Do not modify the logic of this function. Naively, it looks like the + # algorithm that is used when the XP operator is non-diagonal (which + # involves squaring) gives the correct output for diagonal XP operators + # as well. However, that is not true. Counter example given by Mark + # Webster is the operator -I, where the faulty method would give the + # degree 2, while the actual degree is 1. + return np.where(self.is_diagonal(), lcm, lcm_square) + + def fundamental_phase(self) -> np.ndarray: + """Return the fundamental phase of the XP operator. + + Note: + This method is adapted from method XPFundamentalPhase from + XPFpackage: https://github.com/m-webster/XPFpackage, originally + developed by Mark Webster. The original code is licensed under the + GNU General Public License v3.0 and Mark Webster has given + permission to use the code under the Apache License v2.0. + + Returns: + np.ndarray: Fundamental phase of BaseXPPauli + + Examples: + >>> a = BaseXPPauli(matrix=np.array([1, 0, 1, 1, 5, 3, 5, 4], dtype=np.int64), + ... phase_exp=4, precision=6) + >>> a.fundamental_phase() + array([0]) + + See also: + _fundamental_phase + """ + return self._fundamental_phase() + + def _fundamental_phase(self) -> np.ndarray: + """Return the fundamental phase of the XP operator. + + Note: + This method is adapted from method XPFundamentalPhase from + XPFpackage: https://github.com/m-webster/XPFpackage, originally + developed by Mark Webster. The original code is licensed under the + GNU General Public License v3.0 and Mark Webster has given + permission to use the code under the Apache License v2.0. + + See also: + _degree, _power + """ + deg = self._degree() + return self._power(deg)._phase_exp + + def reset_eigenvalue(self, inplace: bool = False) -> "BaseXPPauli": + """Returns the adjusted XP operator such that +1 is an eigenvalue of it. + + Note: + This method is adapted from method XPSetEval from XPFpackage: + https://github.com/m-webster/XPFpackage, originally developed by + Mark Webster. The original code is licensed under the GNU General + Public License v3.0 and Mark Webster has given permission to use + the code under the Apache License v2.0. + + Args: + inplace: If True, adjust BaseXPPauli in place, else return a new + BaseXPPauli. Defaults to False + + Returns: + BaseXPPauli: XP operator with +1 as an eigenvalue + + Examples: + >>> a = BaseXPPauli(matrix=np.array([1, 0, 1, 1, 0, 1, 0, 4], dtype=np.int64), + ... phase_exp=4, precision=6) + >>> value = a.reset_eigenvalue() + >>> value.matrix + array([[1, 0, 1, 1, 0, 1, 0, 4]], dtype=np.int64) + >>> value._phase_exp + array([1]) + + See also: + _reset_eigenvalue + """ + return self._reset_eigenvalue(inplace) + + def _reset_eigenvalue(self, inplace: bool = False) -> "BaseXPPauli": + """Returns the adjusted XP operator such that +1 is an eigenvalue of it. + + Note: + This method is adapted from method XPSetEval from XPFpackage: + https://github.com/m-webster/XPFpackage, originally developed by + Mark Webster. The original code is licensed under the GNU General + Public License v3.0 and Mark Webster has given permission to use + the code under the Apache License v2.0. + Args: + inplace: If True, adjust BaseXPPauli in place, else return a new + BaseXPPauli. Defaults to False + + Returns: + BaseXPPauli: XP operator with +1 as an eigenvalue + + See also: + _fundamental_phase, _degree + """ + fphase = self._fundamental_phase() + deg = self._degree() + new_phase = np.mod(self._phase_exp - np.floor_divide(fphase, deg), 2 * self.precision) + if not inplace: + return BaseXPPauli(matrix=self.matrix, phase_exp=new_phase, precision=self.precision) + self._phase_exp = new_phase + return self + + +# --------------------------------------------------------------------- +# Evolution by Clifford gates +# --------------------------------------------------------------------- + + +# pylint: disable=unused-argument +def _evolve_h(base_xp_pauli: "BaseXPPauli", qubit: int) -> "BaseXPPauli": + """_summary_""" + pass + + +# pylint: disable=unused-argument +def _evolve_s(base_xp_pauli: "BaseXPPauli", qubit: int) -> "BaseXPPauli": + """_summary_""" + pass + + +# pylint: disable=unused-argument +def _evolve_sdg(base_xp_pauli: "BaseXPPauli", qubit: int) -> "BaseXPPauli": + """_summary_""" + pass + + +# pylint: disable=unused-argument +def _evolve_i(base_xp_pauli: "BaseXPPauli", qubit: int) -> "BaseXPPauli": + """_summary_""" + pass + + +def _evolve_x(base_xp_pauli: "BaseXPPauli", qubit: int) -> "BaseXPPauli": + """_summary_""" + pass + + +def _evolve_y(base_xp_pauli: "BaseXPPauli", qubit: int) -> "BaseXPPauli": + """_summary_""" + pass + + +def _evolve_z(base_xp_pauli: "BaseXPPauli", qubit: int) -> "BaseXPPauli": + """_summary_""" + pass + + +def _evolve_cx(base_xp_pauli: "BaseXPPauli", qctrl: int, qtrgt: int) -> "BaseXPPauli": + """_summary_""" + pass + + +def _evolve_cz( # pylint: disable=invalid-name + base_xp_pauli: "BaseXPPauli", q1: int, q2: int # pylint: disable=invalid-name +) -> "BaseXPPauli": + """_summary_""" + pass + + +def _evolve_cy(base_xp_pauli: "BaseXPPauli", qctrl: int, qtrgt: int) -> "BaseXPPauli": + """_summary_""" + pass + + +def _evolve_swap( # pylint: disable=invalid-name + base_xp_pauli: "BaseXPPauli", q1: int, q2: int # pylint: disable=invalid-name +) -> "BaseXPPauli": + """_summary_""" + pass diff --git a/qiskit_qec/operators/pauli_list.py b/qiskit_qec/operators/pauli_list.py index b90aeb61..96f6be49 100644 --- a/qiskit_qec/operators/pauli_list.py +++ b/qiskit_qec/operators/pauli_list.py @@ -10,13 +10,18 @@ # copyright notice, and modified files need to carry a notice indicating # that they have been altered from the originals. # Part of the QEC framework -"""Module fo Pauli List""" +# +# This code is based on the paper: "The XP Stabiliser Formalism: a +# Generalisation of the Pauli Stabiliser Formalism with Arbitrary Phases", Mark +# A. Webster, Benjamin J. Brown, and Stephen D. Bartlett. Quantum 6, 815 +# (2022). +"""Module for Pauli List""" import numbers from collections import defaultdict from typing import Iterable, List, Tuple, Union import numpy as np -import retworkx as rx +import rustworkx as rx from qiskit.exceptions import QiskitError from qiskit.quantum_info.operators.custom_iterator import CustomIterator from qiskit.quantum_info.operators.mixins import GroupMixin, LinearMixin diff --git a/qiskit_qec/operators/random.py b/qiskit_qec/operators/random.py index 14d8a704..3975c8ea 100644 --- a/qiskit_qec/operators/random.py +++ b/qiskit_qec/operators/random.py @@ -13,6 +13,7 @@ Random symplectic operator functions """ +from typing import Union import numpy as np from numpy.random import default_rng @@ -22,6 +23,7 @@ from qiskit_qec.operators.pauli import Pauli from qiskit_qec.operators.pauli_list import PauliList +from qiskit_qec.operators.xp_pauli import XPPauli def random_pauli(num_qubits, group_phase=False, seed=None): @@ -196,6 +198,53 @@ def random_clifford(num_qubits, seed=None): return Clifford(StabilizerTable(table, phase)) +def random_xppauli( + num_qubits: int, + precision: int = None, + seed: Union[int, np.random.Generator, None] = None, +) -> XPPauli: + """Return a random XPPauli. + + Note: + This method is adapted from XPFpackage: + https://github.com/m-webster/XPFpackage, originally developed by + Mark Webster. The original code is licensed under the GNU General + Public License v3.0 and Mark Webster has given permission to use + the code under the Apache License v2.0. + + Args: + num_qubits (int): the number of qubits. + precision (int): Precision of XP operators. Must be an integer + greater than or equal to 2. + seed (int or np.random.Generator): Optional. Set a fixed seed or + generator for RNG. + + Examples: + >>> value = random_xppauli(num_qubits=3, precision=8) + + Returns: + XPPauli: a random XPPauli + """ + if seed is None: + rng = np.random.default_rng() + elif isinstance(seed, np.random.Generator): + rng = seed + else: + rng = default_rng(seed) + z = rng.integers(precision, size=num_qubits, dtype=np.int64) + x = rng.integers(2, size=num_qubits, dtype=bool) + # TODO: Need to decide whether we will add an argument group_phase in + # analogy with random_pauli. If yes, its implementation goes here. + + # Mark's code randomizes phase modulo 2*precision. + phase = rng.integers(2 * precision, dtype=np.int64) + xppauli = XPPauli(data=np.concatenate((z, x)), phase_exp=phase, precision=precision) + return xppauli + + +# TODO: def random_xppauli_list(): + + def _sample_qmallows(n, rng=None): """Sample from the quantum Mallows distribution""" diff --git a/qiskit_qec/operators/xp_pauli.py b/qiskit_qec/operators/xp_pauli.py new file mode 100644 index 00000000..9cdf048f --- /dev/null +++ b/qiskit_qec/operators/xp_pauli.py @@ -0,0 +1,511 @@ +# This code is part of Qiskit. +# +# (C) Copyright IBM 2017, 2020 +# +# 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. +# Part of the QEC framework +"""Module for XPPauli""" +from typing import Any, List, Optional, Union + +import numpy as np +from qiskit.exceptions import QiskitError +from qiskit.quantum_info.operators.mixins import generate_apidocs +from qiskit_qec.operators.base_xp_pauli import BaseXPPauli + +# from qiskit_qec.utils import xp_pauli_rep + + +class XPPauli(BaseXPPauli): + """`XPPauli` inherits from `BaseXPPauli`""" + + # Set the max XPPauli string size before truncation + _truncate__ = 50 + + # pylint: disable=unused-argument + def __init__( + self, + data: Any, + *, + x: Union[List, np.ndarray, None] = None, + z: Union[List, np.ndarray, None] = None, + phase_exp: Union[int, np.ndarray, None] = None, + precision: int = None, + input_xppauli_encoding: str = BaseXPPauli.EXTERNAL_XP_PAULI_ENCODING, + ): + """XPPauli Init + + Args: + data (str): Still in progress + x ([type], optional): [description]. Defaults to None. + z ([type], optional): [description]. Defaults to None. + phase_exponent ([type], optional): [description]. Defaults to None. + precision: Precision of XP operators. Must be an integer greater than or equal to two. + + Raises: + QiskitError: Something went wrong. + + See also: + BaseXPPauli, XPPauliList + """ + if isinstance(data, np.ndarray): + matrix = np.atleast_2d(data) + if phase_exp is None: + phase_exp = 0 + elif isinstance(data, BaseXPPauli): + matrix = data.matrix[:, :] + phase_exp = data._phase_exp[:] + precision = data.precision + # TODO: elif isinstance(data, (tuple, list)), isinstance(data, str), + # isinstance(data, ScalarOp) may be implemented later, like Pauli class + else: + raise QiskitError("Invalid input data for XPPauli.") + + # Initialize BaseXPPauli + if matrix.shape[0] != 1: + raise QiskitError("Input is not a single XPPauli") + + super().__init__(matrix, phase_exp, precision) + # TODO check if this is needed + # self.vlist = self.matrix[0].tolist() + + # --------------------------------------------------------------------- + # Property Methods + # --------------------------------------------------------------------- + + @property + def name(self): + """Unique string identifier for operation type.""" + return "xppauli" + + # TODO: several @property methods exist in Pauli class, analogous methods + # may be added here later. + + @property + def phase_exp(self): + """Return the group phase exponent for the Pauli.""" + # Convert internal Pauli encoding to external Pauli encoding + # return pauli_rep.change_pauli_encoding( + # self._phase_exp, self.num_y, output_pauli_encoding=BasePauli.EXTERNAL_PAULI_ENCODING + # ) + pass + + @phase_exp.setter + def phase_exp(self, value): + # Convert external Pauli encoding to the internal Pauli Encoding + # self._phase_exp[:] = pauli_rep.change_pauli_encoding( + # value, + # self.num_y, + # input_pauli_encoding=BasePauli.EXTERNAL_PAULI_ENCODING, + # output_pauli_encoding=pauli_rep.INTERNAL_PAULI_ENCODING, + # same_type=False, + # ) + pass + + @property + def phase(self): + """Return the complex phase of the Pauli""" + # return pauli_rep.exp2cpx(self.phase_exp, input_encoding=BasePauli.EXTERNAL_PHASE_ENCODING) + pass + + @property + def x(self): + """The x vector for the XPPauli.""" + return self.matrix[:, : self.num_qubits][0] + + @x.setter + def x(self, val): + self.matrix[:, : self.num_qubits][0] = val + + @property + def z(self): + """The z vector for the XPPauli.""" + return self.matrix[:, self.num_qubits :][0] + + @z.setter + def z(self, val): + self.matrix[:, self.num_qubits :][0] = val + + # --------------------------------------------------------------------- + # BaseOperator methods + # --------------------------------------------------------------------- + + def compose( + self, + other: Union["XPPauli", BaseXPPauli], + qargs: Optional[List[int]] = None, + front: bool = False, + inplace: bool = False, + ) -> "XPPauli": + """Return the operator composition with another XPPauli. + + Note: + This method is adapted from method XPMul from XPFpackage: + https://github.com/m-webster/XPFpackage, originally developed by + Mark Webster. The original code is licensed under the GNU General + Public License v3.0 and Mark Webster has given permission to use + the code under the Apache License v2.0. + + Args: + other (XPPauli): a XPPauli object. + qargs (list or None): Optional, qubits to apply dot product + on (default: None). + front (bool): If True compose using right operator multiplication, + instead of left multiplication [default: False]. + inplace (bool): If True update in-place (default: False). + + Returns: + XPPauli: The composed XPPauli. + + Raises: + QiskitError: if other cannot be converted to an operator, or has + incompatible dimensions for specified subsystems, or + if precision of other does not match precision of self. + + .. note:: + Composition (``&``) by default is defined as `left` matrix multiplication for + matrix operators, while :meth:`dot` is defined as `right` matrix + multiplication. That is that ``A & B == A.compose(B)`` is equivalent to + ``B.dot(A)`` when ``A`` and ``B`` are of the same type. + + Setting the ``front=True`` kwarg changes this to `right` matrix + multiplication and is equivalent to the :meth:`dot` method + ``A.dot(B) == A.compose(B, front=True)``. + + Examples: + >>> a = XPPauli(data=np.array([0, 1, 0, 0, 2, 0], dtype=np.int64), phase_exp=6, precision=4) + >>> b = XPPauli(data=np.array([1, 1, 1, 3, 3, 0], dtype=np.int64), phase_exp=2, precision=4) + >>> value = a.compose(b) + >>> value.matrix + array([[1, 0, 1, 3, 3, 0]], dtype=np.int64) + >>> value._phase_exp + array([6]) + + See also: + _compose + """ + if qargs is None: + qargs = getattr(other, "qargs", None) + if not isinstance(other, XPPauli): + other = XPPauli(other) + if self.precision != other.precision: + raise QiskitError("Precision of the two XPPaulis to be multiplied must be the same.") + + return XPPauli(super().compose(other, qargs=qargs, front=front, inplace=inplace)) + + # --------------------------------------------------------------------- + + def unique_vector_rep(self) -> "XPPauli": + """Convert the XPPauli operator into unique vector form, ie + phase_exp in Z modulo 2*precision, x in Z_2, z in Z modulo + precision. + + Note: + This method is adapted from method XPRound from XPFpackage: + https://github.com/m-webster/XPFpackage, originally developed by + Mark Webster. The original code is licensed under the GNU General + Public License v3.0 and Mark Webster has given permission to use + the code under the Apache License v2.0. + + Returns: + XPPauli: Unique vector representation of self + + Examples: + >>> a = XPPauli(matrix=np.array([0, 3, 1, 6, 4, 3], dtype=np.int64), + ... phase_exp=11, precision=4) + >>> a = a.unique_vector_rep() + >>> a.matrix + np.array([[0, 1, 1, 2, 0, 3]], dtype=np.int64) + >>> a._phase_exp + array([3]) + + See also: + _unique_vector_rep + """ + return XPPauli(super().unique_vector_rep()) + + def rescale_precision(self, new_precision: int, inplace: bool = False) -> "XPPauli": + """Rescale the generalized symplectic vector components + of XPPauli operator to the new precision. Returns the rescaled XPPauli object. + + Note: + This method is adapted from method XPSetN from XPFpackage: + https://github.com/m-webster/XPFpackage, originally developed by + Mark Webster. The original code is licensed under the GNU General + Public License v3.0 and Mark Webster has given permission to use + the code under the Apache License v2.0. + + Args: + new_precision: The target precision in which XPPauli is to be expressed + inplace: If True, rescale XPPauli in place, else return a new XPPauli. + Defaults to False + + Returns: + XPPauli: Resultant of rescaling the precision of XPPauli + + Raises: + QiskitError: If it is not possible to express XPPauli in new_precision + + Examples: + >>> a = XPPauli( + ... data=np.array([1, 1, 1, 0, 0, 0, 0, 0, 0, 4, 0, 0, 0, 0], dtype=np.int64), + ... phase_exp=12, precision=8) + >>> a = a.rescale_precision(new_precision=2) + >>> a.matrix + array([[1, 1, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0]], dtype=np.int64) + >>> a._phase_exp + array([3]) + + See also: + _rescale_precision + """ + return XPPauli(super().rescale_precision(new_precision, inplace)) + + def antisymmetric_op(self, int_vec: np.ndarray) -> "XPPauli": + """Return the antisymmetric operator corresponding to an integer vector, + with precision specified by the XP operator. + + Note: + This method is adapted from method XPD from XPFpackage: + https://github.com/m-webster/XPFpackage, originally developed by + Mark Webster. The original code is licensed under the GNU General + Public License v3.0 and Mark Webster has given permission to use + the code under the Apache License v2.0. + + Args: + int_vec: An integer vector + + Returns: + XPPauli: The antisymmetric operator + + Examples: + >>> a = XPPauli( + ... data=np.array([0, 0, 0, 0, 0, 0, 0, 0, 1, 2, 3, 3, 3, 3], dtype=np.int64), + ... phase_exp=0, precision=8) + >>> value = a.antisymmetric_op(data.z) + >>> value.matrix + array([0, 0, 0, 0, 0, 0, 0, 0, -1, -2, -3, -3, -3, -3], dtype=np.int64) + >>> value._phase_exp + array([15]) + + See also: + _antisymmetric_op + """ + return XPPauli(super().antisymmetric_op(int_vec)) + + def inverse(self) -> "XPPauli": + """Return the inverse of the XP operator. + + Note: + This method is adapted from method XPInverse from XPFpackage: + https://github.com/m-webster/XPFpackage, originally developed by + Mark Webster. The original code is licensed under the GNU General + Public License v3.0 and Mark Webster has given permission to use + the code under the Apache License v2.0. + + Returns: + XPPauli: Inverse of XPPauli + + Examples: + >>> a = BaseXPPauli( + ... matrix=np.array([0, 0, 0, 1, 0, 1, 1, 5, 5, 6, 1, 1, 4, 0], dtype=np.int64), + ... phase_exp=1, precision=8) + >>> value = a.inverse() + >>> value.matrix + array([0, 0, 0, 1, 0, 1, 1, 3, 3, 2, 1, 7, 4, 0], dtype=np.int64) + >>> value._phase_exp + array([5]) + + See also: + _inverse + """ + return XPPauli(super().inverse()) + + def power(self, n: Union[int, list, np.ndarray]) -> "XPPauli": + """Return the XP operator of specified precision raised to the power n. + + Note: + This method is adapted from method XPPower from XPFpackage: + https://github.com/m-webster/XPFpackage, originally developed by + Mark Webster. The original code is licensed under the GNU General + Public License v3.0 and Mark Webster has given permission to use + the code under the Apache License v2.0. + + Args: + n: The power to which XPPauli is to be raised + + Returns: + XPPauli: XPPauli raised to the power n + + Examples: + >>> a = XPPauli( + ... data=np.array([1, 0, 1, 1, 5, 3, 5, 4], dtype=np.int64), + ... phase_exp=4, precision=6) + >>> value = a.power(n=5) + >>> value.matrix + array([1, 0, 1, 1, 5, 3, 5, 4], dtype=np.int64) + >>> value._phase_exp + array([4]) + + See also: + _power + """ + return XPPauli(super().power(n)) + + def conjugate( + self, other: Union["XPPauli", BaseXPPauli], front: bool = True, inplace: bool = False + ) -> "XPPauli": + """Return the conjugation of two XP operators. + + For single XP operators, this means + + A.conjugate(B, front=True) = A . B . A^{-1}, + + where . is the XP Pauli multiplication and A^{-1} is the inverse of A. + + Likewise, + + A.conjugate(B, front=False) = B . A . B^{-1}. + + For a list of XP operators, conjugation is performed element-wise: + + [A_1, ..., A_k].conjugate([B_1, ..., B_k]) = [A_1.conjugate(B_1), ..., A_k.conjugate(B_k)]. + + TODO: This method currently only supports conjugation of two XP operator + lists of the same length. + + Note: + This method is adapted from method XPConjugate from XPFpackage: + https://github.com/m-webster/XPFpackage, originally developed by + Mark Webster. The original code is licensed under the GNU General + Public License v3.0 and Mark Webster has given permission to use + the code under the Apache License v2.0. + + Args: + other: an XP operator + front (bool, optional): Whether to conjugate in front (True) or + behind (False), defaults to True + inplace (bool, optional): Whether to perform the conjugation in + place (True) or to return a new BaseXPPauli (False), defaults to + False + + Returns: + XPPauli: Conjugated XP operator + + Examples: + >>> a = XPPauli(data=np.array([1, 0, 1, 1, 5, 3, 5, 4], dtype=np.int64), + ... phase_exp=4, precision=6) + >>> b = XPPauli(data=np.array([1, 0, 0, 1, 4, 1, 0, 1], dtype=np.int64), + ... phase_exp=11, precision=6) + >>> value = a.conjugate(b) + >>> value.matrix + array([[1, 0, 0, 1, 0, 1, 0, 1]], dtype=np.int64) + >>> value._phase_exp + array([3]) + + See also: + _conjugate + """ + if not isinstance(other, XPPauli): + other = XPPauli(other) + + return XPPauli(super().conjugate(other, front=front, inplace=inplace)) + + def commutator( + self, other: Union["XPPauli", BaseXPPauli], front: bool = True, inplace: bool = False + ) -> "XPPauli": + """Return the commutator of two XP operators. + + For single XP operators, this means + + A.commutator(B, front=True) = [A, B] = A . B . A^{-1} . B^{-1}, + + where . is the XP Pauli multiplication and A^{-1} is the inverse of A. + + Likewise, + + A.commutator(B, front=False) = [B, A] = B . A . B^{-1}. A^{-1}. + + For a list of XP operators, commutator is computed element-wise: + + [A_1, ..., A_k].commutator([B_1, ..., B_k]) = [A_1.commutator(B_1), ..., A_k.commutator(B_k)]. + + TODO: This method currently only supports commutator of two XP operator + lists of the same length. + + Note: + This method is adapted from method XPCommutator from XPFpackage: + https://github.com/m-webster/XPFpackage, originally developed by + Mark Webster. The original code is licensed under the GNU General + Public License v3.0 and Mark Webster has given permission to use + the code under the Apache License v2.0. + + Args: + other: an XP operator + front (bool, optional): Whether self is the first element in the + commutator (True) or second (False), defaults to True + inplace (bool, optional): Whether to compute the commutator in + place (True) or to return a new BaseXPPauli (False), defaults to + False + + Returns: + XPPauli: Commutator of XP operators + + Examples: + >>> a = XPPauli(data=np.array([1, 0, 1, 1, 5, 3, 5, 4], dtype=np.int64), + ... phase_exp=4, precision=6) + >>> b = XPPauli(data=np.array([1, 0, 0, 1, 4, 1, 0, 1], dtype=np.int64), + ... phase_exp=11, precision=6) + >>> value = a.commutator(b) + >>> value.matrix + array([[0, 0, 0, 0, 4, 0, 0, 0]], dtype=np.int64) + >>> value._phase_exp + array([8]) + + See also: + _commutator + """ + if not isinstance(other, XPPauli): + other = XPPauli(other) + + return XPPauli(super().commutator(other, front=front, inplace=inplace)) + + def reset_eigenvalue(self, inplace: bool = False) -> "XPPauli": + """Returns the adjusted XP operator such that +1 is an eigenvalue of it. + + Note: + This method is adapted from method XPSetEval from XPFpackage: + https://github.com/m-webster/XPFpackage, originally developed by + Mark Webster. The original code is licensed under the GNU General + Public License v3.0 and Mark Webster has given permission to use + the code under the Apache License v2.0. + + Args: + inplace: If True, adjust XPPauli in place, else return a new XPPauli. + Defaults to False + + Returns: + XPPauli: XP operator with +1 as an eigenvalue + + Examples: + >>> a = XPPauli(data=np.array([1, 0, 1, 1, 0, 1, 0, 4], dtype=np.int64), + ... phase_exp=4, precision=6) + >>> value = a.reset_eigenvalue() + >>> value.matrix + array([[1, 0, 1, 1, 0, 1, 0, 4]], dtype=np.int64) + >>> value._phase_exp + array([1]) + + See also: + _reset_eigenvalue + """ + return XPPauli(super().reset_eigenvalue(inplace)) + + +# Update docstrings for API docs +generate_apidocs(XPPauli) diff --git a/qiskit_qec/operators/xp_pauli_list.py b/qiskit_qec/operators/xp_pauli_list.py new file mode 100644 index 00000000..996be7bc --- /dev/null +++ b/qiskit_qec/operators/xp_pauli_list.py @@ -0,0 +1,673 @@ +# This code is part of Qiskit. +# +# (C) Copyright IBM 2017, 2020 +# +# 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. +# Part of the QEC framework +# +# This code is based on the paper: "The XP Stabiliser Formalism: a +# Generalisation of the Pauli Stabiliser Formalism with Arbitrary Phases", Mark +# A. Webster, Benjamin J. Brown, and Stephen D. Bartlett. Quantum 6, 815 +# (2022). +"""Module for XPPauli List""" +import numbers +from typing import Iterable, List, Tuple, Union, Optional + +import numpy as np +from qiskit.exceptions import QiskitError +from qiskit.quantum_info.operators.mixins import GroupMixin, LinearMixin +from qiskit_qec.operators.base_xp_pauli import BaseXPPauli +from qiskit_qec.utils import xp_pauli_rep + + +# pylint: disable=unused-argument +# pylint: disable=no-member +class XPPauliList(BaseXPPauli, LinearMixin, GroupMixin): + """`XPPauliList` inherits from `BaseXPPauli`""" + + # Set the max number of qubits * xppaulis before string truncation + _truncate__ = 2000 + + def __init__( + self, + data: Union[BaseXPPauli, np.ndarray, Tuple[np.ndarray], Iterable, None] = None, + phase_exp: Union[int, np.ndarray, None] = None, + precision: Union[int, np.ndarray] = None, + *, + input_pauli_encoding: str = BaseXPPauli.EXTERNAL_XP_PAULI_ENCODING, + input_qubit_order: str = "right-to-left", + tuple_order: str = "zx", + ) -> None: + """Inits a XPPauliList + + Args: + data (str): List of XPPauli Operators. + phase_exp (int, optional): i**phase_exp. Defaults to 0. + input_qubit_order (str, optional): Order to read pdata. Defaults to "right-to-left". + precision: Precision of XP operators. Must be an integer/array of + integers greater than or equal to 2. + + Raises: + QiskitError: Something went wrong. + + See also: + BaseXPPauli, XPPauli + """ + if data is None: + matrix = np.empty(shape=(0, 0), dtype=np.bool_) + phase_exp = np.empty(shape=(0,), dtype=np.int8) + elif isinstance(data, BaseXPPauli): + matrix = data.matrix + phase_exp = data._phase_exp + precision = data.precision + # TODO elif isinstance(data, StabilizerTable), elif isinstance(data, PauliTable) + elif isinstance(data, np.ndarray): + if data.size == 0: + matrix = np.empty(shape=(0, 0), dtype=np.bool_) + phase_exp = np.empty(shape=(0,), dtype=np.int8) + # TODO elif isinstance(data[0], str): + else: + if phase_exp is None: + phase_exp = 0 + matrix, phase_exp, precision = xp_pauli_rep.from_array( + data, phase_exp, precision, input_pauli_encoding=input_pauli_encoding + ) + # TODO elif isinstance(data, tuple) + else: + # TODO Conversion as iterable of Paulis + pass + + super().__init__(matrix, phase_exp, precision) + + # TODO + # self.paulis = [ + # Pauli(self.matrix[i], phase_exp=self._phase_exp[i]) for i in range(self.matrix.shape[0]) + # ] + + # --------------------------------------------------------------------- + # Init Methods + # --------------------------------------------------------------------- + + # --------------------------------------------------------------------- + # Property Methods + # --------------------------------------------------------------------- + + @property + def phase(self): + """Return the phase vector of the XPPauliList. + + Note: This is different from the quantum_info phase property which + instead returns the phase_exp + """ + # TODO + pass + + @phase.setter + def phase(self, phase): + """Set the phase vector of the XPPauliList + + Args: + phase (numpy.ndarray or complex numbers): Array of phases, + phases must be one of [1,-1, 1j, -1j] + """ + # TODO + pass + + @property + def shape(self): + """The full shape of the :meth:`array`""" + return self._num_xppaulis, self.num_qubits + + @property + def size(self): + """The number of XPPauli rows in the table.""" + return self._num_xppaulis + + @property + def num_xppaulis(self) -> int: + """Returns the number of XPPauli's in List""" + return self._num_xppaulis + + @property + def phase_exp(self): + """Return the phase exponent vector of the XPPauliList""" + # TODO + pass + + @phase_exp.setter + def phase_exp(self, phase_exp, input_phase_encoding=BaseXPPauli.EXTERNAL_PHASE_ENCODING): + """Set the phase exponent vector of the XPPauliList. Note that this method + converts the phase exponents directly and does not take into account the + number of Y paulis in the representation. + + Args: + phase_exp (_type_): _description_ + input_phase_encoding (_type_, optional): _description_. Defaults to + BaseXPPauli.EXTERNAL_PHASE_ENCODING. + """ + # TODO + pass + + @property + def settings(self): + """Return settings.""" + return {"data": self.to_labels()} + + # --------------------------------------------------------------------- + # Magic Methods and related methods + # --------------------------------------------------------------------- + + def __getitem__(self, index): + """Return a view of the XPPauliList.""" + # Returns a view of specified rows of the XPPauliList + # This supports all slicing operations the underlying array supports. + # TODO + pass + + def getaslist(self, slc: Union[numbers.Integral, slice]) -> List["XPPauli"]: + """_summary_ + + Returns: + _type_: _description_ + """ + # TODO + pass + + def __setitem__(self, index, value): + """Update XPPauliList.""" + # TODO + pass + + def __repr__(self): + """Display representation.""" + return self._truncated_str(True) + + def __str__(self): + """Print representation.""" + return self._truncated_str(False) + + def _truncated_str(self, show_class): + # TODO + pass + + def __array__(self, dtype=None): + """Convert to numpy array""" + # pylint: disable=unused-argument + shape = (len(self),) + 2 * (2**self.num_qubits,) + ret = np.zeros(shape, dtype=complex) + for i, mat in enumerate(self.matrix_iter()): + ret[i] = mat + return ret + + def __eq__(self, other): + """Entrywise comparison of XPPauli equality.""" + if not isinstance(other, XPPauliList): + other = XPPauliList(other) + if not isinstance(other, BaseXPPauli): + return False + return self._eq(other) + + def __len__(self): + """Return the number of XPPauli rows in the table.""" + return self._num_xppaulis + + def _add(self, other, qargs=None): + """summary""" + pass + + # ---- + # + # ---- + + # --------------------------------------------------------------------- + # BaseOperator methods + # --------------------------------------------------------------------- + + def tensor(self, other): + """Return the tensor product with each XPPauli in the list. + + Args: + other (XPPauliList): another XPPauliList. + + Returns: + XPPauliList: the list of tensor product XPPaulis. + + Raises: + QiskitError: if other cannot be converted to a XPPauliList, does + not have either 1 or the same number of XPPaulis as + the current list. + """ + # TODO + pass + + def compose( + self, + other: "XPPauliList", + qargs: Optional[list] = None, + front: bool = False, + inplace: bool = False, + ) -> "XPPauliList": + """Return the composition self∘other for each XPPauli in the list. + + Note: + This method is adapted from the method XPMul from XPFpackage: + https://github.com/m-webster/XPFpackage, originally developed by + Mark Webster. The original code is licensed under the GNU General + Public License v3.0 and Mark Webster has given permission to use + the code under the Apache License v2.0. + + Args: + other (XPPauliList): another XPPauliList. + qargs (None or list): qubits to apply dot product on (Default: None). + front (bool): If True use `dot` composition method [default: False]. + inplace (bool): If True update in-place (default: False). + + Returns: + XPPauliList: the list of composed XPPaulis. + + Raises: + QiskitError: if other cannot be converted to a XPPauliList, does + not have either 1 or the same number of XPPaulis as + the current list, has the wrong number of qubits + for the specified qargs, or if precision of other + does not match precision of self. + + Examples: + >>> a = XPPauliList( + ... data=np.array([[0, 1, 0, 0, 2, 0], [0, 1, 0, 0, 2, 0]], dtype=np.int64), + ... phase_exp=np.array([6, 6]), precision=4) + >>> b = XPPauliList( + ... data=np.array([[1, 1, 1, 3, 3, 0], [1, 1, 1, 3, 3, 0]], dtype=np.int64), + ... phase_exp=np.array([2, 2]), precision=4) + >>> value = a.compose(b) + >>> value.matrix + array([[1, 0, 1, 3, 3, 0], [1, 0, 1, 3, 3, 0]], dtype=np.int64) + >>> value._phase_exp + array([6, 6]) + + See also: + _compose + """ + if qargs is None: + qargs = getattr(other, "qargs", None) + if not isinstance(other, XPPauliList): + other = XPPauliList(other) + if len(other) not in [1, len(self)]: + raise QiskitError( + "Incompatible XPPauliLists. Other list must " + "have either 1 or the same number of XPPaulis." + ) + if self.precision != other.precision: + raise QiskitError( + "Precision of the two XPPauliLists to be multiplied must be the same." + ) + + return XPPauliList(super().compose(other, qargs=qargs, front=front, inplace=inplace)) + + def rescale_precision(self, new_precision: int, inplace: bool = False) -> "XPPauliList": + """Rescale the generalized symplectic vector components + of XPPauli operator to the new precision. Returns the rescaled XPPauli object. + + Note: + This method is adapted from method XPSetN from XPFpackage: + https://github.com/m-webster/XPFpackage, originally developed by + Mark Webster. The original code is licensed under the GNU General + Public License v3.0 and Mark Webster has given permission to use + the code under the Apache License v2.0. + + Args: + new_precision: The target precision in which XPPauli is to be expressed + inplace: If True, rescale XPPauliList in place, else return a new + XPPauliList. Defaults to False + + Returns: + XPPauliList: Resultant of rescaling the precision of XPPauliList + + Raises: + QiskitError: If it is not possible to express XPPauliList in new_precision + + Examples: + >>> matrix1 = np.array([1, 1, 1, 0, 0, 0, 0, 0, 0, 4, 0, 0, 0, 0], dtype=np.int64) + >>> phase_exp1 = 12 + >>> matrix2 = np.array([1, 1, 0, 1, 0, 0, 0, 0, 0, 4, 0, 0, 0, 6], dtype=np.int64) + >>> phase_exp2 = 8 + >>> precision = 8 + >>> new_precision = 4 + >>> matrix = np.array([matrix1, matrix2]) + >>> phase_exp = np.array([phase_exp1, phase_exp2]) + >>> xppauli_list = XPPauliList(data=matrix, phase_exp=phase_exp, precision=precision) + >>> rescaled_xppauli_list = xppaulilist.rescale_precision(new_precision=new_precision) + >>> rescaled_xppauli_list.matrix + array([[1, 1, 1, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0], + [1, 1, 0, 1, 0, 0, 0, 0, 0, 2, 0, 0, 0, 3]] dtype=np.int64) + >>> rescaled_xppauli_list._phase_exp + array([6, 4]) + + See also: + _rescale_precision + """ + return XPPauliList(super().rescale_precision(new_precision, inplace)) + + def antisymmetric_op(self, int_vec: np.ndarray) -> "XPPauliList": + """Return the antisymmetric operators corresponding to the list of + integer vectors, with precision specified by BaseXPPauli. + + Note: + This method is adapted from method XPD from XPFpackage: + https://github.com/m-webster/XPFpackage, originally developed by + Mark Webster. The original code is licensed under the GNU General + Public License v3.0 and Mark Webster has given permission to use + the code under the Apache License v2.0. + + Args: + int_vec (np.ndarray): Array containing integer vectors + + Returns: + XPPauliList: The antisymmetric operators corresponding to the input vectors + + Examples: + >>> matrix = np.array([[0, 0, 0, 0, 0, 0, 0, 0, 1, 2, 3, 3, 3, 3], + ... [0, 0, 0, 0, 0, 0, 0, 3, 1, 2, 3, 7, 6, 3]], dtype=np.int64) + >>> phase_exp = np.array([0, 0]) + >>> precision = 8 + >>> xppauli_list = XPPauliList(data=matrix, phase_exp=phase_exp, precision=precision) + >>> value = xppauli_list.antisymmetric_op(xppauli_list.z) + >>> value.matrix + array([[0, 0, 0, 0, 0, 0, 0, 0, -1, -2, -3, -3, -3, -3], + [0, 0, 0, 0, 0, 0, 0, -3, -1, -2, -3, -7, -6, -3]], dtype=np.int64) + >>> value._phase_exp + array([15, 25]) + + See also: + _antisymmetric_op + """ + return XPPauliList(super().antisymmetric_op(int_vec)) + + def inverse(self) -> "XPPauli": + """Return the inverse of the list of XP operators. + + Note: + This method is adapted from method XPInverse from XPFpackage: + https://github.com/m-webster/XPFpackage, originally developed by + Mark Webster. The original code is licensed under the GNU General + Public License v3.0 and Mark Webster has given permission to use + the code under the Apache License v2.0. + + Returns: + XPPauliList: Inverse of XPPauliList + + Examples: + >>> matrix = np.array([[1, 1, 0, 1, 1, 0, 1, 2, 4, 4, 3, 1, 6, 1], + ... [0, 1, 0, 0, 1, 0, 1, 7, 7, 3, 4, 6, 2, 7]], dtype=np.int64) + >>> phase_exp = np.array([1, 0]) + >>> precision = 8 + >>> xppauli_list = XPPauliList(data=matrix, phase_exp=phase_exp, precision=precision) + >>> value = xppauli_list.inverse() + >>> value.matrix + array([[1, 1, 0, 1, 1, 0, 1, 2, 4, 4, 3, 1, 2, 1], + [0, 1, 0, 0, 1, 0, 1, 1, 7, 5, 4, 6, 6, 7]], dtype=np.int64) + >>> value._phase_exp + np.array([9, 8]) + + See also: + _inverse + """ + return XPPauliList(super().inverse()) + + def power(self, n: Union[int, list, np.ndarray]) -> "XPPauliList": + """Return the XP operators of specified precision raised to the power n. + + Note: + This method is adapted from method XPPower from XPFpackage: + https://github.com/m-webster/XPFpackage, originally developed by + Mark Webster. The original code is licensed under the GNU General + Public License v3.0 and Mark Webster has given permission to use + the code under the Apache License v2.0. + + Args: + n: The power to which XPPauliList is to be raised + + Returns: + XPPauliList: XPPauliList raised to the power n + + Examples: + >>> matrix = np.array([[1, 0, 1, 1, 5, 3, 5, 4], + ... [1, 0, 1, 1, 5, 4, 1, 5]], dtype=np.int64) + >>> phase_exp = np.array([4, 3]) + >>> precision = 6 + >>> n = np.array([5, 3]) + >>> xppauli_list = XPPauliList(data=matrix, phase_exp=phase_exp, precision=precision) + >>> value = xppauli_list.power(n=n) + >>> value.matrix + array([[1, 0, 1, 1, 5, 3, 5, 4], + [1, 0, 1, 1, 5, 0, 1, 5]], dtype=np.int64) + >>> value._phase_exp + np.array([4, 7]) + + See also: + _power + """ + return XPPauliList(super().power(n)) + + def conjugate( + self, other: "XPPauliList", front: bool = True, inplace: bool = False + ) -> "XPPauliList": + """Return the conjugation of two XP operators. + + For single XP operators, this means + + A.conjugate(B, front=True) = A . B . A^{-1}, + + where . is the XP Pauli multiplication and A^{-1} is the inverse of A. + + Likewise, + + A.conjugate(B, front=False) = B . A . B^{-1}. + + For a list of XP operators, conjugation is performed element-wise: + + [A_1, ..., A_k].conjugate([B_1, ..., B_k]) = [A_1.conjugate(B_1), ..., A_k.conjugate(B_k)]. + + TODO: This method currently only supports conjugation of two XP operator + lists of the same length. + + Note: + This method is adapted from method XPConjugate from XPFpackage: + https://github.com/m-webster/XPFpackage, originally developed by + Mark Webster. The original code is licensed under the GNU General + Public License v3.0 and Mark Webster has given permission to use + the code under the Apache License v2.0. + + Args: + other: List of XP operators to be conjugated with self + front (bool, optional): Whether to conjugate in front (True) or + behind (False), defaults to True + inplace (bool, optional): Whether to perform the conjugation in + place (True) or to return a new BaseXPPauli (False), defaults to + False + + Returns: + XPPauliList: List of conjugated XP operators + + Raises: + QiskitError: Other list must have either 1 or the same number of + XPPaulis + + Examples: + >>> a = XPPauliList(data=np.array([[1, 0, 1, 1, 5, 3, 5, 4], + ... [1, 0, 1, 0, 1, 5, 2, 0]], dtype=np.int64), + ... phase_exp=np.array([4, 7]), precision=6) + >>> b = XPPauliList(data=np.array([[1, 0, 0, 1, 4, 1, 0, 1], + ... [0, 1, 1, 0, 1, 3, 0, 5]], dtype=np.int64), + ... phase_exp=np.array([11, 2]), precision=6) + >>> value = a.conjugate(b) + >>> value.matrix + array([[1, 0, 0, 1, 0, 1, 0, 1], + [0, 1, 1, 0, 5, 5, 4, 5]], dtype=np.int64) + >>> value._phase_exp + array([3, 10]) + + See also: + _conjugate + """ + if not isinstance(other, XPPauliList): + other = XPPauliList(other) + if len(other) not in [1, len(self)]: + raise QiskitError( + "Incompatible XPPauliLists. Other list must " + "have either 1 or the same number of XPPaulis." + ) + + return XPPauliList(super().conjugate(other, front=front, inplace=inplace)) + + def commutator( + self, other: "XPPauliList", front: bool = True, inplace: bool = False + ) -> "XPPauliList": + """Return the commutator of two XP operators. + + For single XP operators, this means + + A.commutator(B, front=True) = [A, B] = A . B . A^{-1} . B^{-1}, + + where . is the XP Pauli multiplication and A^{-1} is the inverse of A. + + Likewise, + + A.commutator(B, front=False) = [B, A] = B . A . B^{-1}. A^{-1}. + + For a list of XP operators, commutator is computed element-wise: + + [A_1, ..., A_k].commutator([B_1, ..., B_k]) = [A_1.commutator(B_1), ..., A_k.commutator(B_k)]. + + TODO: This method currently only supports commutator of two XP operator + lists of the same length. + + Note: + This method is adapted from method XPCommutator from XPFpackage: + https://github.com/m-webster/XPFpackage, originally developed by + Mark Webster. The original code is licensed under the GNU General + Public License v3.0 and Mark Webster has given permission to use + the code under the Apache License v2.0. + + Args: + other: List of XP operators to be in the commutator with self + front (bool, optional): Whether self is the first element in the + commutator (True) or second (False), defaults to True + inplace (bool, optional): Whether to compute the commutator in + place (True) or to return a new BaseXPPauli (False), defaults to + False + + Returns: + XPPauliList: List of commutators of XP operators + + Raises: + QiskitError: Other list must have either 1 or the same number of + XPPaulis + + Examples: + >>> a = XPPauliList(data=np.array([[1, 0, 1, 1, 5, 3, 5, 4], + ... [1, 0, 1, 0, 1, 5, 2, 0]], dtype=np.int64), + ... phase_exp=np.array([4, 7]), precision=6) + >>> b = XPPauliList(data=np.array([[1, 0, 0, 1, 4, 1, 0, 1], + ... [0, 1, 1, 0, 1, 3, 0, 5]], dtype=np.int64), + ... phase_exp=np.array([11, 2]), precision=6) + >>> value = a.commutator(b) + >>> value.matrix + array([[0, 0, 0, 0, 4, 0, 0, 0], + [0, 0, 0, 0, 4, 4, 2, 0]], dtype=np.int64) + >>> value._phase_exp + array([8, 8]) + + See also: + _commutator + """ + if not isinstance(other, XPPauliList): + other = XPPauliList(other) + if len(other) not in [1, len(self)]: + raise QiskitError( + "Incompatible XPPauliLists. Other list must " + "have either 1 or the same number of XPPaulis." + ) + + return XPPauliList(super().commutator(other, front=front, inplace=inplace)) + + def reset_eigenvalue(self, inplace: bool = False) -> "XPPauliList": + """Returns the list of adjusted XP operators such that +1 is an eigenvalue of them. + + Note: + This method is adapted from method XPSetEval from XPFpackage: + https://github.com/m-webster/XPFpackage, originally developed by + Mark Webster. The original code is licensed under the GNU General + Public License v3.0 and Mark Webster has given permission to use + the code under the Apache License v2.0. + + Args: + inplace: If True, adjust XPPauliList in place, else return a new + XPPauliList. Defaults to False + + Returns: + XPPauliList: XP operators with +1 as an eigenvalue + + Examples: + >>> a = XPPauliList(data=np.array([[0, 0, 1, 1, 1, 0, 4, 2], + ... [1, 1, 0, 1, 0, 1, 0, 4]], dtype=np.int64), + ... phase_exp=np.array([7, 4]), precision=6) + >>> value = a.reset_eigenvalue() + >>> value.matrix + array([[0, 0, 1, 1, 1, 0, 4, 2], + [1, 1, 0, 1, 0, 1, 0, 4]], dtype=np.int64) + >>> value._phase_exp + array([6, 1]) + + See also: + _reset_eigenvalue + """ + return XPPauliList(super().reset_eigenvalue(inplace)) + + # def transpose(self): + # """Return the transpose of each XPPauli in the list.""" + # # TODO + # pass + + # def adjoint(self): + # """Return the adjoint of each XPPauli in the list.""" + # # TODO + # pass + + # --------------------------------------------------------------------- + # Utility methods + # --------------------------------------------------------------------- + + # --------------------------------------------------------------------- + # Custom Iterators + # --------------------------------------------------------------------- + + # --------------------------------------------------------------------- + # Class methods + # --------------------------------------------------------------------- + + @classmethod + def from_symplectic(cls, z, x, phase_exp=0): + """Construct a XPPauliList from a symplectic data. + + Args: + z (np.ndarray): 2D boolean Numpy array. + x (np.ndarray): 2D boolean Numpy array. + phase_exp (np.ndarray or None): Optional, 1D integer array from Z_4. + + Returns: + XPPauliList: the constructed XPPauliList. + + Note: Initialization this way will copy matrices and not reference them. + + TODO: Fix this method to be more general and not in old form only + (i.e. include matrix inputs ...) + """ + # TODO + pass diff --git a/qiskit_qec/utils/__init__.py b/qiskit_qec/utils/__init__.py index 52bfd589..353cb46b 100644 --- a/qiskit_qec/utils/__init__.py +++ b/qiskit_qec/utils/__init__.py @@ -28,6 +28,7 @@ indexer pauli_rep + xp_pauli_rep """ -from . import indexer, pauli_rep +from . import indexer, pauli_rep, xp_pauli_rep diff --git a/qiskit_qec/utils/pauli_rep.py b/qiskit_qec/utils/pauli_rep.py index 1c49812e..126e0091 100644 --- a/qiskit_qec/utils/pauli_rep.py +++ b/qiskit_qec/utils/pauli_rep.py @@ -102,13 +102,13 @@ # Different string syntax formats are available. The they are "product" syntax and # "index" syntax. "product" syntax represents a Pauli operator of the form # :math: $p * T_1 \otimes T_2 \otimes ... \otimes T_n$ as :math" $pT1T2T3...Tn$. See the -# following exmaples: +# following examples: # # -iX \otimes Y \otimes Z -> -iXYZ # X \otimes Y \otimes Z \otimes I \otimes I \otimes I -> XYZII # # The index syntax only represents the non identity Paulis. Index syntax specifically -# indiciates the index that the Pauli's are acting on. Following Qiskit's current internal +# indicates the index that the Pauli's are acting on. Following Qiskit's current internal # indexing: # # -iX \otimes Y \otimes Z -> -iZ0Y1X2 @@ -256,7 +256,7 @@ def _is_pattern(string, pattern): def get_phase_encodings() -> List[str]: - """Returns the availble phase encodings + """Returns the available phase encodings Returns: encoding: List of available phase encodings @@ -328,7 +328,7 @@ def split_pauli_enc(encoding: str) -> Tuple[str, str]: """Splits the Pauli encoding into the phase and tensor encodings Args: - encoding: Pauli encpoding + encoding: Pauli encoding Raises: QiskitError: Encoding not valid @@ -336,7 +336,7 @@ def split_pauli_enc(encoding: str) -> Tuple[str, str]: Returns: phase_enc, tensor_enc: phase encoding and tensor encoding - Exampes: + Examples: >>> encoding = "iXZ' >>> split_pauli_encoding(encoding) ('i', 'XZ') @@ -359,7 +359,7 @@ def _split_pauli_enc(encoding: str) -> Tuple[str, str]: def get_phase_enc(encoding: str) -> str: - """Returns the phase encodeing part of the Pauli encoding string + """Returns the phase encoding part of the Pauli encoding string Args: encoding: Pauli encoding string @@ -372,7 +372,7 @@ def get_phase_enc(encoding: str) -> str: def get_tensor_enc(encoding: str) -> str: - """Returns the tensor encodeing part of the Pauli encoding string + """Returns the tensor encoding part of the Pauli encoding string Args: encoding: Pauli encoding string @@ -716,7 +716,7 @@ def squeeze(array_: Any, scalar: bool = False) -> bool: >>> squeeze(numpy.array([[[[1,2,3,4,5]]]])) array([1,2,3,4,5]) """ - array_ = np.squeeze(array_) + array_ = np.squeeze(np.array(array_, dtype=object)) if array_.shape == () and scalar is True: return array_.item() else: @@ -1238,7 +1238,7 @@ def exp2expstr( different encodings have a specific syntaxs. Args: - phase_exp: Phase encosings to convert to string representations + phase_exp: Phase encodings to convert to string representations input_encoding: Encoding of the input phase exponents same_type (optional): Scalar/Vector return flag. Defaults to True. diff --git a/qiskit_qec/utils/xp_pauli_rep.py b/qiskit_qec/utils/xp_pauli_rep.py new file mode 100644 index 00000000..7461b75b --- /dev/null +++ b/qiskit_qec/utils/xp_pauli_rep.py @@ -0,0 +1,1127 @@ +# This code is part of Qiskit. +# +# (C) Copyright IBM 2017, 2020 +# +# 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. + +# Part of the QEC framework +# +# This code is based on the paper: "The XP Stabiliser Formalism: a +# Generalisation of the Pauli Stabiliser Formalism with Arbitrary Phases", Mark +# A. Webster, Benjamin J. Brown, and Stephen D. Bartlett. Quantum 6, 815 +# (2022). +""" +N-qubit XPPauli Representation Encodings and Conversion Module +""" + +# pylint: disable=invalid-name,anomalous-backslash-in-string +# pylint: disable=bad-docstring-quotes # for deprecate_function decorator + +import numbers +import re +from typing import Any, Iterable, List, Optional, Tuple, Union + +import numpy as np +from qiskit.circuit import Gate +from qiskit.quantum_info.operators.scalar_op import ScalarOp +from qiskit.exceptions import QiskitError +from scipy.sparse import csr_matrix + + +# ------------------------------------------------------------------------------- +# Module Variables/States +# ------------------------------------------------------------------------------- + +# Set the internal encodings +# The internal encoding cannot be changed by changing this constant +# These constants are for reference only and do not change the behavior of +# the XPPauli methods. See [ref] for details on the different encodings +# TODO: Include ref for above. + +INTERNAL_TENSOR_ENCODING = "XP" +INTERNAL_PHASE_ENCODING = "w" +INTERNAL_XP_PAULI_ENCODING = INTERNAL_PHASE_ENCODING + INTERNAL_TENSOR_ENCODING + +DEFAULT_EXTERNAL_TENSOR_ENCODING = "XP" +DEFAULT_EXTERNAL_PHASE_ENCODING = "w" +DEFAULT_EXTERNAL_XP_PAULI_ENCODING = ( + DEFAULT_EXTERNAL_PHASE_ENCODING + DEFAULT_EXTERNAL_TENSOR_ENCODING +) + +# Set the external encodings +# The external encodings may be changed via the phase_rep_format, +# symp_rep_format, or pauli_rep_format methods. Any method changing +# these values must make sure to update the others +# Tensor encodings are: 'XZ', 'XZY', 'ZX', 'YZX' +# Phase encodings are: 'i', '-i', 'is', '-is' +# See [ref] for details on the different encodings +# TODO: Include ref for above. +# TODO update this comment after formats have been finalized. + +# w is exp(pi*i/N), as defined in XP Formalism paper. If the phase exponent is +# p, then the operator's phase is w**p, where p is an integer between and +# including 0 and 2N-1. P is diag(1,w**2). +PHASE_ENCODINGS = ["w"] +PHASE_ENCODINGS_IMI = [] +PHASE_ENCODINGS_ISMIS = [] +TENSOR_ENCODINGS = ["XP"] +Y_TENSOR_ENCODINGS = [] +XP_PAULI_ENCODINGS = [ + "wXP", +] +XP_PAULI_ENCODINGS_SPLIT = { + "wXP": ("w", "XP"), +} + +# Different string syntax formats are available. The they are "product" syntax and +# "index" syntax. "product" syntax represents a XPPauli operator of the form +# :math: $p * T_1 \otimes T_2 \otimes ... \otimes T_n$ as :math" $pT1T2T3...Tn$. See the +# following exmaples: +# +# -iX \otimes Y \otimes Z -> -iXYZ +# X \otimes Y \otimes Z \otimes I \otimes I \otimes I -> XYZII +# +# The index syntax only represents the non identity XPPaulis. Index syntax specifically +# indiciates the index that the XPPauli's are acting on. Following Qiskit's current internal +# indexing: +# +# -iX \otimes Y \otimes Z -> -iZ0Y1X2 +# X \otimes Y \otimes Z \otimes I \otimes I \otimes I -> Z3Y4X5 +# TODO update this comment/following REGEX after formats have been finalized. + +# -i, -1j, +1, ... +COMPLEX_REGEX = r"[\-+]?1?[ij]?" +# (i,2), (1j, 0), ( +i , 2 ), ... +ENC_I_REGEX = r"\(\s*[+]?[ij]\s*\,\s*[0123]\s*\)" +# (-i,2), (-1j, 0), ( -i , 2 ), ... +ENC_MI_REGEX = r"\(\s*\-1?[ij]\s*\,\s*[0123]\s*\)" +# (i,0)(-1,1), (1j,1)(-1 , 0), ... +ENC_IS_REGEX = r"\(\s*1?[ij]\s*\,\s*[01]\s*\)\s*\(\s*\-1\s*\,\s*[01]\s*\)" +# (-i,0)(-1,1), (-1j,1)(-1, 0), ... +ENC_MIS_REGEX = r"\(\s*\-1?[ij]\s*\,\s*[01]\s*\)\s*\(\s*\-1\s*\,\s*[01]\s*\)" + +CAP_STAND_ENC_I_REGEX = r"\(i\,([0123])\)" +CAP_STAND_ENC_MI_REGEX = r"\(-i\,([0123])\)" +CAP_STAND_ENC_IS_REGEX = r"\(i\,([01])\)\(-1\,([01])\)" +CAP_STAND_ENC_MIS_REGEX = r"\(-i\,([01])\)\(-1\,([01])\)" + +CAP_STAND_I_PATTERN = re.compile(f"^{CAP_STAND_ENC_I_REGEX}$") +CAP_STAND_MI_PATTERN = re.compile(f"^{CAP_STAND_ENC_MI_REGEX}$") +CAP_STAND_IS_PATTERN = re.compile(f"^{CAP_STAND_ENC_IS_REGEX}$") +CAP_STAND_MIS_PATTERN = re.compile(f"^{CAP_STAND_ENC_MIS_REGEX}$") + +CAP_PHASE_PATTERN = { + "i": CAP_STAND_I_PATTERN, + "-i": CAP_STAND_MI_PATTERN, + "is": CAP_STAND_IS_PATTERN, + "-is": CAP_STAND_MIS_PATTERN, +} + +COMPLEX_PATTERN = re.compile(f"^({COMPLEX_REGEX})$") + +I_PATTERN = re.compile(f"^({ENC_I_REGEX})$") +MI_PATTERN = re.compile(f"^({ENC_MI_REGEX})$") +IS_PATTERN = re.compile(f"^({ENC_IS_REGEX})$") +MIS_PATTERN = re.compile(f"^({ENC_MIS_REGEX})$") + +PHASE_PATTERN = { + "complex": COMPLEX_PATTERN, + "i": I_PATTERN, + "-i": MI_PATTERN, + "is": IS_PATTERN, + "-is": MIS_PATTERN, +} + +XP_PAULI_START_REGEX = r"\(?[IXZY].*" +SPLIT_PATTERN = re.compile(f"^(.*?)({XP_PAULI_START_REGEX})") + +PHASE_REGEX = r"[\-+]?1?[ij]?" +XP_PAULI_REGEX = r"[IXZY]" + +INDEX_REGEX = r".*?[0-9].*" +INDEX_INDICATE_PATTERN = re.compile(f"^{INDEX_REGEX}$") + +LATEX_REGEX = r".*?_\{[0-9].*" +LATEX_INDICATE_PATTERN = re.compile(f"^{LATEX_REGEX}$") + + +ENC_INDEX_XZ_REGEX = r"(\((X[0-9]+|Z[0-9]+|X([0-9]+)Z\3)\))+" +ENC_INDEX_XZY_REGEX = r"([XZY][0-9]+)+" +ENC_INDEX_ZX_REGEX = r"(\((X[0-9]+|Z[0-9]+|Z([0-9]+)X\3)\))+" +ENC_INDEX_YZX_REGEX = r"([XZY][0-9]+)+" + + +INDEX_XZ_PATTERN_CAP = re.compile(r"\((X[0-9]+|Z[0-9]+|X(?:[0-9]+)Z[0-9]+)\)") +INDEX_ZX_PATTERN_CAP = re.compile(r"\((Z[0-9]+|X[0-9]+|Z(?:[0-9]+)X[0-9]+)\)") + +INDEX_XZ_PATTERN = re.compile(f"^{ENC_INDEX_XZ_REGEX}$") +INDEX_XZY_PATTERN = re.compile(f"^{ENC_INDEX_XZY_REGEX}$") +INDEX_ZX_PATTERN = re.compile(f"^{ENC_INDEX_ZX_REGEX}$") +INDEX_YZX_PATTERN = re.compile(f"^{ENC_INDEX_YZX_REGEX}$") + +TENSOR_INDEX_PATTERN = { + "XZ": INDEX_XZ_PATTERN, + "XZY": INDEX_XZY_PATTERN, + "ZX": INDEX_ZX_PATTERN, + "YZX": INDEX_YZX_PATTERN, +} + +ENC_LATEX_XZ_REGEX = r"(\((X_\{[0-9]+\}|Z_\{[0-9]+\}|X_\{([0-9]+\})Z\3)\))+" +ENC_LATEX_XZY_REGEX = r"([XZY]_\{[0-9]+\})+" +ENC_LATEX_ZX_REGEX = r"(\((X_\{[0-9]+\}|Z_\{[0-9]+\}|Z_\{([0-9]+\})X\3)\))+" +ENC_LATEX_YZX_REGEX = r"([XZY]_\{[0-9]+\})+" + + +LATEX_XZ_PATTERN_CAP = re.compile(r"\((X_\{[0-9]+\}|Z_\{[0-9]+\}|X(?:_\{[0-9]+\})Z_\{[0-9]+\})\)") +LATEX_ZX_PATTERN_CAP = re.compile(r"\((Z_\{[0-9]+\}|X_\{[0-9]+\}|Z(?:_\{[0-9]+\})X_\{[0-9]+\})\)") + +LATEX_XZ_PATTERN = re.compile(f"^{ENC_LATEX_XZ_REGEX}$") +LATEX_XZY_PATTERN = re.compile(f"^{ENC_LATEX_XZY_REGEX}$") +LATEX_ZX_PATTERN = re.compile(f"^{ENC_LATEX_ZX_REGEX}$") +LATEX_YZX_PATTERN = re.compile(f"^{ENC_LATEX_YZX_REGEX}$") + +TENSOR_LATEX_PATTERN = { + "XZ": LATEX_XZ_PATTERN, + "XZY": LATEX_XZY_PATTERN, + "ZX": LATEX_ZX_PATTERN, + "YZX": LATEX_YZX_PATTERN, +} + + +ENC_PRODUCT_XZ_CAP = r"\(([XZ]|XZ|I)\)" +ENC_PRODUCT_ZX_CAP = r"\(([ZX]|ZX|I)\)" + +PRODUCT_XZ_PATTERN_CAP = re.compile(f"{ENC_PRODUCT_XZ_CAP}") +PRODUCT_ZX_PATTERN_CAP = re.compile(f"{ENC_PRODUCT_ZX_CAP}") + +ENC_PRODUCT_XZ_REGEX = r"(\(([XZ]|XZ|I)\))+" +ENC_PRODUCT_XZY_REGEX = r"[XZYI]+" +ENC_PRODUCT_ZX_REGEX = r"(\(([ZX]|ZX|I)\))+" +ENC_PRODUCT_YZX_REGEX = r"[XZYI]+" + +PRODUCT_XZ_PATTERN = re.compile(f"^{ENC_PRODUCT_XZ_REGEX}$") +PRODUCT_XZY_PATTERN = re.compile(f"^{ENC_PRODUCT_XZY_REGEX}$") +PRODUCT_ZX_PATTERN = re.compile(f"^{ENC_PRODUCT_ZX_REGEX}$") +PRODUCT_YZX_PATTERN = re.compile(f"^{ENC_PRODUCT_YZX_REGEX}$") + +TENSOR_PRODUCT_PATTERN = { + "XZ": PRODUCT_XZ_PATTERN, + "XZY": PRODUCT_XZY_PATTERN, + "ZX": PRODUCT_ZX_PATTERN, + "YZX": PRODUCT_YZX_PATTERN, +} + +PRODUCT_SYNTAX = 0 +INDEX_SYNTAX = 1 +LATEX_SYNTAX = 2 +XP_SYMPLECTIC_SYNTAX = 3 +DEFAULT_SYNTAX = 0 +SYNTAX_TO_TEXT = ["product", "index", "latex", "XP symplectic"] + +DEFAULT_QUBIT_ORDER = "right-to-left" +QUBIT_ORDERS = ["right-to-left", "left-to-right"] + + +def _is_pattern(string, pattern): + """_summary_""" + return bool(pattern.search(string)) + + +# ------------------------------------------------------------------------------- +# Encoding lists +# ------------------------------------------------------------------------------- + + +def get_phase_encodings() -> List[str]: + """Returns the available phase encodings + + Returns: + encoding: List of available phase encodings + + Examples: + >>> get_phase_encodings() + ['w'] + + See Also + get_tensor_encodings, get_pauli_encodings + """ + return PHASE_ENCODINGS + + +def get_tensor_encodings() -> List[str]: + """Returns the available tensor encodings + + Returns: + encoding: List of available tensor encodings + + Examples: + >>> get_tensor_encodings() + ['XP'] + + See Also: + get_phase_encodings, get_pauli_encodings + """ + return TENSOR_ENCODINGS + + +def get_xp_pauli_encodings() -> List[str]: + """Returns the available XPPauli encodings + + Returns: + encodings : List of available XPPauli encodings + + Example: + >>> get_xp_pauli_encodings() + ['wXP'] + + See Also: + get_phase_encodings, get_tensor_encodings + + """ + return XP_PAULI_ENCODINGS + + +# ------------------------------------------------------------------------------- +# Encoding Methods and Conversions +# ------------------------------------------------------------------------------- + + +def split_xp_pauli_enc(encoding: str) -> Tuple[str, str]: + """Splits the XPPauli encoding into the phase and tensor encodings + + Args: + encoding: XPPauli encoding + + Raises: + QiskitError: Encoding not valid + + Returns: + phase_enc, tensor_enc: phase encoding and tensor encoding + + Examples: + >>> encoding = "wXP" + >>> split_xp_pauli_encoding(encoding) + ('w', 'XP') + + See Also: + _split_xp_pauli_encoding + """ + if encoding not in XP_PAULI_ENCODINGS: + raise QiskitError(f"Encoding not valid: {encoding}") + return _split_xp_pauli_enc(encoding) + + +def _split_xp_pauli_enc(encoding: str) -> Tuple[str, str]: + """Splits the XPPauli encoding into the phase and tensor encodings + + Args: + encoding: XPPauli encoding string + """ + return XP_PAULI_ENCODINGS_SPLIT[encoding] + + +def get_phase_enc(encoding: str) -> str: + """Returns the phase encoding part of the XPPauli encoding string + + Args: + encoding: XPPauli encoding string + + Returns: + phase_enc: phase encoding + + Examples: + >>> encoding = "wXP" + >>> get_phase_enc(encoding) + 'w' + """ + phase_part, _ = split_xp_pauli_enc(encoding) + return phase_part + + +def get_tensor_enc(encoding: str) -> str: + """Returns the tensor encoding part of the XPPauli encoding string + + Args: + encoding: XPPauli encoding string + + Returns: + tensor_enc: tensor encoding + + Examples: + >>> encoding = "wXP" + >>> get_tensor_enc(encoding) + 'XP' + """ + _, tensor_part = split_xp_pauli_enc(encoding) + return tensor_part + + +# TODO depending on what encoding formats are decided, these need to be +# implemented or removed. +# pylint: disable=unused-argument +def change_xp_pauli_encoding( + phase_exp: Any, + y_count: Union[np.array, int] = 0, + *, + input_xp_pauli_encoding: str = INTERNAL_XP_PAULI_ENCODING, + output_xp_pauli_encoding: str = DEFAULT_EXTERNAL_XP_PAULI_ENCODING, + same_type=True, +) -> Any: + """_summary_""" + pass + + +# pylint: disable=unused-argument +def _change_xp_pauli_encoding( + phase_exponent: np.ndarray, + y_count: np.ndarray, + input_xp_pauli_encoding: str, + output_xp_pauli_encoding: str, +) -> Any: + """_summary_""" + pass + + +# TODO depending on what encoding formats are decided, these need to be +# implemented or removed. +# pylint: disable=unused-argument +def stand_phase_str( + phase_str: str, same_type: bool = True, imaginary: str = "i" +) -> Union[np.ndarray, str]: + """_summary_""" + pass + + +# pylint: disable=unused-argument +def _stand_phase_str(phase_string: np.ndarray, imaginary: str) -> np.ndarray: + """_summary_""" + pass + + +# pylint: disable=unused-argument +def is_scalar(obj: Any, scalar_pair: bool = False): + """_summary_""" + pass + + +# pylint: disable=unused-argument +def squeeze(array_: Any, scalar: bool = False) -> bool: + """_summary_""" + pass + + +# pylint: disable=unused-argument +def is_exp_type(phase_exp: Any, input_phase_encoding: str) -> bool: + """_summary_""" + pass + + +# Conversion function between XPPauli phases as + +# complex string: "0-1j", ... +# complex type: 0-1j, ... +# encoded string: "(-i, 1)", ... +# encoded type: (-i,0)(-1,1), ... + +# Methods convering between type do not change the encoding scheme if present +# Encoding schemes can be change via the exp2exp method + +# This method are all vector enabled and have two forms: method and _method. +# The method versions have all the checks and conversions in them. The raw +# methods assume correct input formats and generally do very little checking +# or conversions. +# +# The same_type parameter is used to allow scalar values to be output +# when scalar values are input + +# ---------------------------------------------------------------------- +# pylint: disable=unused-argument +def cpxstr2cpx( + cpx_str: Union[str, np.ndarray, List[str]], + same_type: bool = True, +) -> Union[np.ndarray, numbers.Complex]: + """_summary_""" + pass + + +# pylint: disable=unused-argument +def _cpxstr2cpx(cpx_string: np.ndarray) -> np.ndarray: + """_summary_""" + pass + + +# ---------------------------------------------------------------------- +# pylint: disable=unused-argument +def cpx2cpxstr( + cpx: Union[numbers.Complex, np.ndarray], same_type: bool = True, ones: bool = False +) -> Union[str, np.ndarray]: + """_summary_""" + pass + + +# pylint: disable=unused-argument +def _cpx2cpxstr(cpx: np.ndarray, one_str: str) -> Union[str, np.ndarray]: + """_summary_""" + pass + + +# ---------------------------------------------------------------------- + + +# pylint: disable=unused-argument +def exp2cpx( + phase_exp: Any, input_encoding: str, same_type: bool = True +) -> Union[np.ndarray, numbers.Complex]: + """_summary_""" + pass + + +# pylint: disable=unused-argument +def _exp2cpx(phase_exp: np.ndarray, input_encoding: str) -> np.ndarray: + """_summary_""" + pass + + +# ---------------------------------------------------------------------- + + +# pylint: disable=unused-argument +def cpx2exp( + cpx: numbers.Complex, output_encoding: str, same_type: bool = True, roundit: bool = True +) -> Union[np.ndarray, Tuple[numbers.Integral, numbers.Integral], numbers.Integral]: + """_summary_""" + pass + + +# pylint: disable=unused-argument +def _cpx2exp(cpx: numbers.Complex, encoding: str) -> np.ndarray: + """_summary_""" + pass + + +# ---------------------------------------------------------------------- + + +# pylint: disable=unused-argument +def expstr2exp(exp_str: Union[np.ndarray, str], same_type: bool = True) -> Any: + """_summary_""" + pass + + +# pylint: disable=unused-argument +def _expstr2exp(exp_string, encoding: str) -> np.ndarray: + """_summary_""" + pass + + +# ---------------------------------------------------------------------- + + +# pylint: disable=unused-argument +def exp2expstr( + phase_exp: Any, + input_encoding: str = DEFAULT_EXTERNAL_XP_PAULI_ENCODING, + same_type: bool = True, +) -> Union[np.ndarray, str]: + """Converts encoded phases (exponents) to their string representations + + Note: This method does more than apply str method as the string representations of the + different encodings have a specific syntaxs. + + Args: + phase_exp: Phase encodings to convert to string representations + input_encoding: Encoding of the input phase exponents. Defaults to + DEFAULT_EXTERNAL_XP_PAULI_ENCODING. + same_type (optional): Scalar/Vector return flag. Defaults to True. + + Raises: + QiskitError: Invalid phase exponent encoding + + Returns: + exp_str: string representations of given phase exponents + + Examples: + TODO + + See Also: + _exp2expstr + """ + if input_encoding not in PHASE_ENCODINGS: + raise QiskitError(f"Invalid phase exponent encoding: {input_encoding}") + + phase_exp = np.atleast_1d(phase_exp) + + return _exp2expstr(phase_exp, input_encoding) + + +def _exp2expstr(phase_exp: np.ndarray, encoding: str) -> np.ndarray: + """Converts encoded phases (exponents) to their string representations + + Note: This method does more than apply str as the string representations of the + different encodings have a specific syntax. + + Args: + phase_exp: Phase encosings to convert to string representations + encoding: Encoding of the input phase exponents + + Raises: + QiskitError: Invalid phase exponent encoding + + Returns: + exp_str: The encoding is not supported + + Examples: + TODO + + See Also: + exp2expstr + """ + if encoding == "w": + return np.array(["(w," + str(item) + ")" for item in phase_exp]) + else: + raise QiskitError(f"The encoding {encoding} is not supported.") + + +# ---------------------------------------------------------------------- + + +# pylint: disable=unused-argument +def exp2exp( + phase_exp: Union[np.ndarray, Any], + input_encoding=INTERNAL_PHASE_ENCODING, + output_encoding=DEFAULT_EXTERNAL_PHASE_ENCODING, + same_type: bool = True, +): + """_summary_""" + pass + + +# pylint: disable=unused-argument +def _exp2exp(phase_exp, input_encoding, output_encoding): + """_summary_""" + pass + + +# ---------------------------------------------------------------------- + + +# pylint: disable=unused-argument +def cpxstr2expstr( + cpx_str: Union[np.ndarray, str], encoding: str, same_type: bool = True +) -> Union[np.ndarray, str]: + """_summary_""" + pass + + +# pylint: disable=unused-argument +def expstr2cpxstr( + exp_str: Union[np.ndarray, str], encoding: str, same_type: bool = True, ones: bool = False +) -> Union[str, np.ndarray]: + """_summary_""" + pass + + +# pylint: disable=unused-argument +def cpxstr2exp(cpx_str: Union[np.ndarray, str], encoding: str, same_type: bool = True) -> Any: + """_summary_""" + pass + + +# pylint: disable=unused-argument +def exp2cpxstr( + phase_exp: Any, encoding: str, same_type: bool = True, ones: bool = False +) -> Union[str, np.ndarray]: + """_summary_""" + pass + + +# pylint: disable=unused-argument +def expstr2cpx( + phase_str: Union[np.ndarray, str], encoding: str, same_type: bool = True +) -> Union[np.ndarray, numbers.Complex]: + """_summary_""" + pass + + +# pylint: disable=unused-argument +def cpx2expstr( + cpx: Union[np.ndarray, numbers.Complex], encoding: str, same_type: bool = True +) -> Union[np.ndarray, str]: + """_summary_""" + pass + + +# ---------------------------------------------------------------------- + + +# pylint: disable=unused-argument +def str2exp( + phase_str: Union[np.ndarray, str], encoding: str, same_type: bool = True +) -> Union[np.ndarray, str]: + """_summary_""" + pass + + +# ---------------------------------------------------------------------- + + +# pylint: disable=unused-argument +def _str2exp(phase_str: np.ndarray, encoding: str) -> np.ndarray: + """_summary_""" + pass + + +# --------------------------------------------------------------------- +# Label parsing helper functions +# --------------------------------------------------------------------- + + +# pylint: disable=unused-argument +def split_pauli(pauli_str: str, same_type: bool = True) -> Tuple[str, str]: + """_summary_""" + pass + + +# pylint: disable=unused-argument +def _split_pauli(pauli_str: str): + """_summary_""" + + def _split(p_str): + pass + + pass + + +# ---------------------------------------------------------------------- + + +# pylint: disable=unused-argument +def encode_of_phase_str(phase_str: str, same_type: bool = True) -> Union[np.ndarray, str]: + """_summary_""" + pass + + +# pylint: disable=unused-argument +def _encode_of_phase_str(phase_str: str) -> np.ndarray: + """_summary_""" + + # pylint: disable=unused-variable + def _find_encode(p_str): + pass + + pass + + +# ---------------------------------------------------------------------- + + +# pylint: disable=unused-argument +def encode_of_tensor_str( + tensor_str: str, encoded: bool = True, same_type: bool = True +) -> List[Tuple[List, Union[str, int]]]: + """_summary_""" + pass + + +# pylint: disable=unused-argument +def _encode_of_tensor_str( + tensor_str: np.ndarray, encoded: bool +) -> List[Tuple[List, Union[str, int]]]: + """_summary_""" + + # pylint: disable=unused-variable + def _find_encode(t_str: str, ecode: bool): + pass + + pass + + +# ---------------------------------------------------------------------- + + +# pylint: disable=unused-argument +def str2symplectic( + pauli_str: Union[np.ndarray, str], + qubit_order: str = "right-to-left", + output_encoding: Optional[str] = INTERNAL_XP_PAULI_ENCODING, + index_start: int = 0, + same_type: bool = True, +) -> Tuple[np.ndarray, Union[np.array, Any]]: + """_summary_""" + pass + + +# pylint: disable=unused-argument +def _str2symplectic( + pauli_str: np.ndarray, + qubit_order: str, + output_encoding: str, + index_start: int, +) -> Tuple[np.ndarray, np.ndarray]: + """_summary_""" + pass + + +# ---------------------------------------------------------------------- +def xp_symplectic2str( + matrix: np.ndarray, + phase_exp: Any = None, + precision: int = None, + input_encoding: str = INTERNAL_XP_PAULI_ENCODING, + output_phase_encoding: str = None, + no_phase: bool = False, + output_tensor_encoding: str = DEFAULT_EXTERNAL_TENSOR_ENCODING, + syntax: str = INDEX_SYNTAX, + qubit_order: str = "right-to-left", + index_start: int = 0, + same_type: bool = True, + index_str="", +) -> Union[np.ndarray, str]: + """Converts a symplectic matrix and phase to string representations + + Args: + matrix: Generalized symplectic matrix for XP operator + phase_exp (optional): Phase exponent(s) for matrix. A value of + None will lead to unity phases. Defaults to None. + precision: Precision of XP operator. + input_encoding (optional): XPPauli encoding of phase relative to + matrix. Defaults to INTERNAL_XP_PAULI_ENCODING. + output_phase_encoding (optional): Encoding used to represent phases. + A value of None will result in complex phases notation. Defaults + to None. + no_phase (optional): When set to True, no phase will appear no matter + what encoding is selected. + output_tensor_encoding (optional): Encoding of XPPauli tensor + (without phase). Defaults to DEFAULT_EXTERNAL_TENSOR_ENCODING. + syntax (optional): Syntax of pauli tensor. Values are + PRODUCT_SYNTAX = 0, INDEX_SYNTAX=1, LATEX_SYNTAX=2 and XP_SYMPLECTIC_SYNTAX=3. + Defaults to INDEX_SYNTAX. + qubit_order (optional): Order in which qubits are read. options are + "right-to-left" and "left-to-right". Defaults to "right-to-left". + index_start (optional): Lowest value for index in index syntax tensors. + Defaults to 0 + same_type (optional): Scalar/Vector return flag. Defaults to True. + index_str (optional): String that get inserted between operator and numbers in + index format. Default is "". + + Raises: + QiskitError: Unsupport syntax + + Returns: + xp_pauli_str: XPPauli strings + + Examples: + >>> precision = 8 + >>> matrix1 = np.array([1, 1, 1, 0, 0, 0, 0, 0, 0, 4, 0, 0, 0, 0], dtype=np.int64) + >>> phase_exp1 = 12 + >>> matrix2 = np.array([1, 1, 1, 0, 0, 0, 0, 0, 0, 2, 3, 0, 0, 0], dtype=np.int64) + >>> phase_exp2 = 2 + >>> matrix = np.array([matrix1, matrix2]) + >>> phase_exp = np.array([phase_exp1, phase_exp2]) + >>> xp_symplectic2str(matrix, phase_exp, precision) + np.array(["XP8((w,12)(X(P,4))2(X)1(X)0)", "XP8((w,2)(P,3)3(X(P,2))2(X)1(X)0)"]) + + >>> xp_symplectic2str(matrix, phase_exp, precision, qubit_order="left-to-right") + np.array(["XP8((w,12)(X)0(X)1(X(P,4))2)", "XP8((w,2)(X)0(X)1(X(P,2))2(P,3)3)"]) + + >>> xp_symplectic2str(matrix, phase_exp, precision, syntax=XP_SYMPLECTIC_SYNTAX) + np.array(["XP8(12|1 1 1 0 0 0 0|0 0 4 0 0 0 0)", "XP8(2|1 1 1 0 0 0 0|0 0 2 3 0 0 0)"]) + + >>> xp_symplectic2str(matrix, phase_exp, precision, no_phase=True) + np.array(["XP8((X(P,4))2(X)1(X)0)", "XP8((P,3)3(X(P,2))2(X)1(X)0)"]) + + >>> xp_symplectic2str(matrix, phase_exp, precision, syntax=PRODUCT_SYNTAX) + np.array(["XP8((w,12)(I)(I)(I)(I)(X(P,4))(X)(X))", "XP8((w,2)(I)(I)(I)(P,3)(X(P,2))(X)(X))"]) + + >>> xp_symplectic2str(matrix, phase_exp, precision, syntax=LATEX_SYNTAX) + np.array( + [ + "XP_{8}((w,12)(XP^{4})_{2}(X)_{1}(X)_{0})", + "XP_{8}((w,2)(P^{3})_{3}(XP^{2})_{2}(X)_{1}(X)_{0})", + ] + ) + """ + matrix = np.atleast_2d(matrix) + num_qubits = matrix.shape[1] >> 1 + matrix[:, 0:num_qubits] = np.mod(matrix[:, 0:num_qubits], 2) + matrix[:, num_qubits:] = np.mod(matrix[:, num_qubits:], precision) + if no_phase: + phase_str = np.full((matrix.shape[0],), "") + else: + if phase_exp is None: + phase_exp = np.zeros(shape=(matrix.shape[0],), dtype=np.int64) + else: + phase_exp = np.atleast_1d(phase_exp) + phase_exp = np.mod(phase_exp, 2 * precision) + + # If multiple phase/tensor encodings are implemented, the conversion + # needs to go here. + + if output_phase_encoding is None: + if syntax != XP_SYMPLECTIC_SYNTAX: + phase_str = exp2expstr(phase_exp, "w") + + tensor_str = [] + + _XPENC = ["(I)", "(X)", "(P,{zexp})", "(X(P,{zexp}))"] + _XP_LATEX_ENC = ["(I)", "(X)", "(P{zexp})", "(XP{zexp})"] + _ENC = {"XP": _XPENC} + + if syntax == PRODUCT_SYNTAX: + for xppauli in matrix: + tmp_tensor_str = "" + for index in range(num_qubits): + tmp_enc = "" + rep = "" + if xppauli[index + num_qubits] > 1: + rep = str(xppauli[index + num_qubits]) + if xppauli[index + num_qubits] == 0: + tmp_enc = _ENC[output_tensor_encoding][xppauli[index]] + else: + tmp_enc = _ENC[output_tensor_encoding][2 + xppauli[index]].replace( + "{zexp}", rep + ) + tmp_enc = tmp_enc.replace("X(P,)", "XP") + tmp_enc = tmp_enc.replace("(P,)", "(P)") + if tmp_enc: + if qubit_order == "left-to-right": + tmp_tensor_str += tmp_enc + else: + tmp_tensor_str = tmp_enc + tmp_tensor_str + + tensor_str.append(tmp_tensor_str) + elif syntax == XP_SYMPLECTIC_SYNTAX: + for i, xppauli in enumerate(matrix): + tmp_tensor_str = "XP" + str(precision) + tmp_tensor_str += ( + "(" + + str(phase_exp[i]) + + "|" + + str(xppauli[:num_qubits])[1:-1] + + "|" + + str(xppauli[num_qubits:])[1:-1] + + ")" + ) + + tensor_str.append(tmp_tensor_str) + elif syntax in (INDEX_SYNTAX, LATEX_SYNTAX): + if syntax == LATEX_SYNTAX: + ind_str_repr = _ind_to_latex_repr + sup_str_repr = _sup_to_latex_repr + _TMP_ENC = {"XP": _XP_LATEX_ENC} + else: + ind_str_repr = str + sup_str_repr = str + _TMP_ENC = _ENC + + for xppauli in matrix: + tmp_tensor_str = "" + for index in range(num_qubits): + if output_tensor_encoding == "XP": + tmp_enc = "" + rep = "" + if xppauli[index + num_qubits] > 1: + rep = sup_str_repr(xppauli[index + num_qubits]) + if xppauli[index] == 1 and xppauli[index + num_qubits] == 0: + tmp_enc = _TMP_ENC[output_tensor_encoding][1] + elif xppauli[index + num_qubits] > 0: + tmp_enc = _TMP_ENC[output_tensor_encoding][2 + xppauli[index]].replace( + "{zexp}", rep + ) + + tmp_enc = tmp_enc.replace("X(P,)", "XP") + tmp_enc = tmp_enc.replace("(P,)", "(P)") + if tmp_enc: + if qubit_order == "left-to-right": + tmp_tensor_str += ( + tmp_enc + index_str + ind_str_repr(index + index_start) + ) + else: + tmp_tensor_str = ( + tmp_enc + + index_str + + ind_str_repr(index + index_start) + + tmp_tensor_str + ) + + tensor_str.append(tmp_tensor_str) + + else: + raise QiskitError(f"Unsupported syntax: {syntax}") + + if syntax != XP_SYMPLECTIC_SYNTAX: + if syntax in (PRODUCT_SYNTAX, INDEX_SYNTAX): + result = [ + "XP" + str(precision) + "(" + p_str + t_str + ")" + for p_str, t_str in zip(phase_str, tensor_str) + ] + else: + result = [ + "XP" + _ind_to_latex_repr(precision) + "(" + p_str + t_str + ")" + for p_str, t_str in zip(phase_str, tensor_str) + ] + else: + result = tensor_str + if matrix.shape[0] == 1 and same_type: + return result[0] + else: + return np.array(result) + + +# ---------------------------------------------------------------------- +# pylint: disable=unused-argument +def str2str( + pauli_str: Union[np.ndarray, str], + syntax_output: str, + phase_encoding_output_string: str = None, + tensor_encoding_output_string: str = DEFAULT_EXTERNAL_TENSOR_ENCODING, + qubit_order_input: str = "right-to-left", + qubit_order_output: str = "right-to-left", + index_start_input: int = 0, + index_start_output: int = 0, +) -> Union[np.ndarray, str]: + """_summary_""" + pass + + +# ---------------------------------------------------------------------- +# Array(s) to Symplectic +# ---------------------------------------------------------------------- + + +def from_array( + matrix: Union[List, Tuple, np.ndarray], + phase_exp: Union[int, List, Tuple, np.ndarray] = None, + precision: Union[int, List, Tuple, np.ndarray] = None, + input_pauli_encoding: Optional[str] = None, +) -> Tuple[np.ndarray, np.ndarray, np.ndarray]: + """Convert array and phase_exp to the matrix, phase_exp and precision in BaseXPPauli's internal + XPPauli encoding (xp_pauli_rep.INTERNAL_XP_PAULI_ENCODING) + Args: + matrix (_type_): _description_ + phase_exp (_type_): _description_ + precision (_type_): Precision of XP operator + input_pauli_encoding: input XPPauli encoding + Returns: + _type_: _description_ + """ + if input_pauli_encoding is None: + input_pauli_encoding = DEFAULT_EXTERNAL_XP_PAULI_ENCODING + if isinstance(matrix, np.ndarray) and matrix.dtype == np.int64: + matrix_data = matrix + else: + matrix_data = np.asarray(matrix, dtype=np.int64) + matrix_data = np.atleast_2d(matrix_data) + # TODO + # if not is_symplectic_matrix_form(matrix_data): + # raise QiskitError("Input matrix not a symplectic matrix or symplectic vector") + if phase_exp is None or (isinstance(phase_exp, numbers.Integral) and phase_exp == 0): + phase_exp = np.zeros(shape=(matrix_data.shape[0],)) + # TODO may need to implement change_pauli_encoding + return matrix_data, phase_exp, precision + + +# pylint: disable=unused-argument +def from_split_array( + x: Union[List, Tuple, np.ndarray], + z: Union[List, Tuple, np.ndarray], + phase_exp: Union[int, List, Tuple, np.ndarray], + input_pauli_encoding: Optional[str] = None, +) -> Tuple[np.ndarray, np.ndarray]: + """_summary_""" + pass + + +# ---------------------------------------------------------------------- +# Symplectic to Complex Matrix conversion +# ---------------------------------------------------------------------- + + +# pylint: disable=unused-argument +def to_cpx_matrix( + matrix: np.ndarray, + phase_exp: Optional[np.ndarray], + pauli_encoding: str = INTERNAL_XP_PAULI_ENCODING, + sparse: bool = False, +) -> Union[np.ndarray, csr_matrix]: + """_summary_""" + pass + + +# pylint: disable=unused-argument +def _to_cpx_matrix( + matrix: np.ndarray, phase_exp: int, num_qubits: int, sparse: bool = False +) -> Union[np.ndarray, csr_matrix]: + """_summary_""" + pass + + +# ---------------------------------------------------------------------- +# Utility functions +# ---------------------------------------------------------------------- + + +# pylint: disable=unused-argument +def indices_to_boolean(indices: Iterable[int], dim: int) -> np.ndarray: + """_summary_""" + pass + + +# ---------------------------------------------------------------------- +def _ind_to_latex_repr(index: int) -> str: + """Adds curly braces and an underscore to an index. + + Args: + index: An integer + + Returns: + str: string in LaTeX syntax + """ + return f"_{{{index}}}" + + +# ---------------------------------------------------------------------- +def _sup_to_latex_repr(superscript: int) -> str: + """Adds curly braces and a caret to a superscript. + + Args: + superscript: An integer + + Returns: + str: string in LaTeX syntax + """ + return f"^{{{superscript}}}" + + +# pylint: disable=unused-argument +def boolean_to_indices(booleans: Iterable[bool]) -> np.ndarray: + """_summary_""" + pass + + +# pylint: disable=unused-argument +def scalar_op2symplectic( + op: ScalarOp, output_encoding: str = DEFAULT_EXTERNAL_PHASE_ENCODING +) -> Tuple[np.ndarray, Union[np.array, Any]]: + """_summary_""" + pass + + +# ---------------------------------------------------------------------- + + +# pylint: disable=unused-argument +def gate2symplectic( + gate: Gate, encoding: str = INTERNAL_XP_PAULI_ENCODING +) -> Tuple[np.ndarray, Union[np.array, Any]]: + """_summary_""" + pass + + +# ---------------------------------------------------------------------- diff --git a/requirements.txt b/requirements.txt index 7177ba80..65a5c2c0 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,8 +1,8 @@ qiskit-terra>=0.21.2 qiskit-aer>=0.11.0 pybind11<=2.9.1 -PyMatching>=0.6.0 -retworkx>=0.11.0 +PyMatching>=0.6.0,!=2.0.0 +rustworkx>=0.11.0 networkx>=2.6.3 sympy>=1.9 numpy>=1.21.0 diff --git a/tests/heavy_hex_codes/test_heavy_hex_decoder.py b/tests/heavy_hex_codes/test_heavy_hex_decoder.py index 7188d4bc..fd93f317 100644 --- a/tests/heavy_hex_codes/test_heavy_hex_decoder.py +++ b/tests/heavy_hex_codes/test_heavy_hex_decoder.py @@ -149,7 +149,7 @@ def test_d3_2(self): self.correct_all_1(c, circ, dec, self.model) def test_d3_3(self): - """Check 3, zx, z, retworkx.""" + """Check 3, zx, z, rustworkx.""" blocks = 3 round_schedule = "zx" basis = "z" @@ -185,7 +185,7 @@ def test_d3_3(self): basis=basis, round_schedule=round_schedule, blocks=blocks, - method="retworkx", + method="rustworkx", uniform=False, ) # self.no_faults_success(c, circ, dec, self.model) diff --git a/tests/linear/test_matrix.py b/tests/linear/test_matrix.py index 9cfc40d1..589ef2f5 100644 --- a/tests/linear/test_matrix.py +++ b/tests/linear/test_matrix.py @@ -24,6 +24,9 @@ augment_mat, rref_complete, rank, + do_row_op, + _howell_complete, + howell_complete, ) @@ -332,3 +335,92 @@ def test_rank(self): dtype=np.bool_, ) self.assertEqual(rank(matrix), 4) + + def test_do_row_op(self): + """Tests do row operation.""" + matrix = np.array( + [[1, 2, 3], [4, 5, 6], [7, 8, 9]], + ) + n = 4 + + swap_mat = do_row_op(matrix, ("swap", [0, 2], []), n) + expected_swap_mat = np.array( + [[7, 8, 9], [4, 5, 6], [1, 2, 3]], + ) + self.assertTrue(np.equal(swap_mat, expected_swap_mat).all()) + + unit_mat = do_row_op(matrix, ("unit", [0], [2]), n) + expected_unit_mat = np.array( + [[2, 0, 2], [4, 5, 6], [7, 8, 9]], + ) + self.assertTrue(np.equal(unit_mat, expected_unit_mat).all()) + + add_mat = do_row_op(matrix, ("add", [0, 2], [2]), n) + expected_add_mat = np.array( + [[3, 2, 1], [4, 5, 6], [7, 8, 9]], + ) + self.assertTrue(np.equal(add_mat, expected_add_mat).all()) + + append_mat = do_row_op(matrix, ("append", [0], [2]), n) + expected_append_mat = np.array( + [[1, 2, 3], [4, 5, 6], [7, 8, 9], [2, 0, 2]], + ) + self.assertTrue(np.equal(append_mat, expected_append_mat).all()) + + update_mat = do_row_op(matrix, ("update", [0, 1], [2, 1, 3, 4]), n) + expected_update_mat = np.array( + [[2, 1, 0], [3, 2, 1], [7, 8, 9]], + ) + self.assertTrue(np.equal(update_mat, expected_update_mat).all()) + + def test_invalid_do_row_op(self): + """Tests invalid do row operation.""" + matrix = np.array( + [[1, 2, 3], [4, 5, 6], [7, 8, 9]], + ) + n = 4 + + fraction_n = 1.1 + negative_n = -2 + invalid_matrix = np.array([[[1, 1], [2, 1]]]) + invalid_op = "invalid" + invalid_row_idx = 3 + + self.assertRaises(QiskitError, do_row_op, matrix, ("unit", [0], [1]), fraction_n) + self.assertRaises(QiskitError, do_row_op, matrix, ("unit", [0], [1]), negative_n) + self.assertRaises(QiskitError, do_row_op, invalid_matrix, ("unit", [0], [1]), n) + self.assertRaises(QiskitError, do_row_op, matrix, (invalid_op, [0], [1]), n) + self.assertRaises(QiskitError, do_row_op, matrix, ("unit", [invalid_row_idx], [1]), n) + + self.assertRaises(QiskitError, do_row_op, matrix, ("swap", [0, 1, 2], []), n) + self.assertRaises(QiskitError, do_row_op, matrix, ("unit", [0], [1, 2]), n) + self.assertRaises(QiskitError, do_row_op, matrix, ("add", [0], [2]), n) + self.assertRaises(QiskitError, do_row_op, matrix, ("append", [0, 2], [1]), n) + self.assertRaises(QiskitError, do_row_op, matrix, ("update", [0, 2], [2, 1, 3]), n) + + def test_howell_complete(self): + """Tests howell.""" + matrix = np.array( + [[8, 5, 5], [0, 9, 8], [0, 0, 10]], + ) + n = 12 + + howell_mat, transform_mat, kernel_mat = _howell_complete(matrix, n) + expected_howell_mat = np.array([[4, 1, 0], [0, 3, 0], [0, 0, 1]]) + expected_transform_mat = np.array([[8, 1, 0], [0, 7, 4], [9, 3, 4]]) + expected_kernel_mat = np.array([[6, 6, 3], [0, 4, 4]]) + self.assertTrue(np.equal(howell_mat, expected_howell_mat).all()) + self.assertTrue(np.equal(transform_mat, expected_transform_mat).all()) + self.assertTrue(np.equal(kernel_mat, expected_kernel_mat).all()) + + def test_invalid_howell_complete(self): + """Tests invalid howell.""" + valid_array = np.array([[1, 2], [3, 4]]) + valid_n = 8 + + negative_n = -2 + self.assertRaises(QiskitError, howell_complete, valid_array, negative_n) + fraction_n = 1.5 + self.assertRaises(QiskitError, howell_complete, valid_array, fraction_n) + invalid_array = np.array([[[1, 1], [2, 1]]]) + self.assertRaises(QiskitError, howell_complete, invalid_array, valid_n) diff --git a/tests/matching/test_circuitmatcher.py b/tests/matching/test_circuitmatcher.py index 001bf36b..abeb85db 100644 --- a/tests/matching/test_circuitmatcher.py +++ b/tests/matching/test_circuitmatcher.py @@ -63,7 +63,7 @@ def setUp(self) -> None: self.x_logical = x_logical def test_no_errors(self): - """Test the case with no errors using retworkx.""" + """Test the case with no errors using rustworkx.""" shots = 100 seed = 100 dec = ThreeBitDecoder( @@ -79,7 +79,7 @@ def test_no_errors(self): "z", "z", 2, - "retworkx", + "rustworkx", False, ) result = execute( @@ -139,7 +139,7 @@ def test_no_errors_pymatching(self): self.assertEqual(failures, 0) def test_correct_single_errors(self): - """Test the case with single faults using retworkx.""" + """Test the case with single faults using rustworkx.""" dec = ThreeBitDecoder( 3, self.x_stabilizers, @@ -153,7 +153,7 @@ def test_correct_single_errors(self): "z", "z", 2, - "retworkx", + "rustworkx", False, ) dec.update_edge_weights(self.pnm) @@ -165,7 +165,7 @@ def test_correct_single_errors(self): self.assertEqual(fail[0], 0) def test_correct_single_errors_uniform(self): - """Test the case with single faults using retworkx.""" + """Test the case with single faults using rustworkx.""" dec = ThreeBitDecoder( 3, self.x_stabilizers, @@ -179,7 +179,7 @@ def test_correct_single_errors_uniform(self): "z", "z", 2, - "retworkx", + "rustworkx", True, ) dec.update_edge_weights(self.pnm) @@ -243,7 +243,7 @@ def test_correct_single_errors_pymatching_uniform(self): self.assertEqual(fail[0], 0) def test_error_pairs(self): - """Test the case with two faults using retworkx.""" + """Test the case with two faults using rustworkx.""" dec = ThreeBitDecoder( 3, self.x_stabilizers, @@ -257,7 +257,7 @@ def test_error_pairs(self): "z", "z", 2, - "retworkx", + "rustworkx", False, ) dec.update_edge_weights(self.pnm) @@ -271,7 +271,7 @@ def test_error_pairs(self): self.assertEqual(failures, 128) def test_error_pairs_uniform(self): - """Test the case with two faults using retworkx.""" + """Test the case with two faults using rustworkx.""" dec = ThreeBitDecoder( 3, self.x_stabilizers, @@ -285,7 +285,7 @@ def test_error_pairs_uniform(self): "z", "z", 2, - "retworkx", + "rustworkx", True, ) dec.update_edge_weights(self.pnm) @@ -352,7 +352,12 @@ def test_error_pairs_propagator_pymatching_uniform(self): corrected_outcomes = dec.process(outcome) fail = temp_syndrome(corrected_outcomes, self.z_logical) failures += fail[0] - self.assertEqual(failures, 156) + # For pymatching v0.x, there are 168 failures, whereas for pymatching >v2.0.0 there are 156. + # The reason for the difference is that many of these test cases have degenerate solutions + # (Both versions of pymatching are giving valid minimum-weight perfect matching solutions, but + # the predictions they make are not always the same when there is more than one valid + # minimum-weight solution.) + self.assertTrue(failures in {156, 168}) if __name__ == "__main__": diff --git a/tests/matching/test_pymatchingmatcher.py b/tests/matching/test_pymatchingmatcher.py index c0cb9c7a..f0653668 100644 --- a/tests/matching/test_pymatchingmatcher.py +++ b/tests/matching/test_pymatchingmatcher.py @@ -1,8 +1,8 @@ -"""Tests for the retworkx matcher subroutines.""" +"""Tests for the rustworkx matcher subroutines.""" import unittest from typing import Dict, Tuple -import retworkx as rx +import rustworkx as rx from qiskit_qec.decoders.pymatching_matcher import PyMatchingMatcher diff --git a/tests/matching/test_repetitionmatcher.py b/tests/matching/test_repetitionmatcher.py index d307cb8b..a8920df0 100644 --- a/tests/matching/test_repetitionmatcher.py +++ b/tests/matching/test_repetitionmatcher.py @@ -37,8 +37,8 @@ def setUp(self) -> None: self.code_circuit_5 = RepetitionCodeCircuit(5, 2) self.z_logical_5 = self.code_circuit.css_z_logical - def test_no_errors(self, method="retworkx"): - """Test the case with no errors using retworkx.""" + def test_no_errors(self, method="rustworkx"): + """Test the case with no errors using rustworkx.""" def gint(c): """Casts to int if possible""" @@ -74,8 +74,8 @@ def test_no_errors_pymatching(self): """Test the case with no errors using pymatching.""" self.test_no_errors(method="pymatching") - def test_correct_single_errors(self, method="retworkx"): - """Test the case with single faults using retworkx.""" + def test_correct_single_errors(self, method="rustworkx"): + """Test the case with single faults using rustworkx.""" for logical in ["0", "1"]: dec = RepetitionDecoder(self.code_circuit, self.pnm, method, False, logical) qc = self.code_circuit.circuit[logical] @@ -91,8 +91,8 @@ def test_correct_single_errors_pymatching(self): """Test the case with two faults using pymatching.""" self.test_correct_single_errors(method="pymatching") - def test_error_pairs(self, dec_method="retworkx", fe_method="stabilizer"): - """Test the case with two faults on a d=5 code using retworkx.""" + def test_error_pairs(self, dec_method="rustworkx", fe_method="stabilizer"): + """Test the case with two faults on a d=5 code using rustworkx.""" expected_failures = {"0": 0, "1": 0} for logical in ["0", "1"]: dec = RepetitionDecoder(self.code_circuit_5, self.pnm, dec_method, False, logical) diff --git a/tests/matching/test_retworkxmatcher.py b/tests/matching/test_retworkxmatcher.py index 4f411bbd..6aca8ec6 100644 --- a/tests/matching/test_retworkxmatcher.py +++ b/tests/matching/test_retworkxmatcher.py @@ -1,13 +1,13 @@ -"""Tests for the retworkx matcher subroutines.""" +"""Tests for the rustworkx matcher subroutines.""" import unittest from typing import Dict, Tuple -import retworkx as rx -from qiskit_qec.decoders.retworkx_matcher import RetworkXMatcher +import rustworkx as rx +from qiskit_qec.decoders.rustworkx_matcher import RustworkxMatcher -class TestRetworkXMatcher(unittest.TestCase): - """Tests for the retworkx matcher subroutines.""" +class TestRustworkxMatcher(unittest.TestCase): + """Tests for the rustworkx matcher subroutines.""" def make_test_graph(self) -> Tuple[rx.PyGraph, Dict[Tuple[int, Tuple[int]], int]]: """Make a basic decoding graph. @@ -29,7 +29,7 @@ def make_test_graph(self) -> Tuple[rx.PyGraph, Dict[Tuple[int, Tuple[int]], int] return graph, idxmap def setUp(self) -> None: - self.rxm = RetworkXMatcher(annotate=True) + self.rxm = RustworkxMatcher(annotate=True) def test_preprocess(self): """Test preprocessing example.""" diff --git a/tests/operators/test_xp_pauli.py b/tests/operators/test_xp_pauli.py new file mode 100644 index 00000000..73155e4a --- /dev/null +++ b/tests/operators/test_xp_pauli.py @@ -0,0 +1,290 @@ +# This code is part of Qiskit. +# +# (C) Copyright IBM 2017, 2020. +# +# 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. + +# pylint: disable=invalid-name + +"""Tests for Pauli operator class.""" + +import unittest + +import numpy as np +from ddt import ddt +from qiskit.test import QiskitTestCase +from qiskit_qec.operators.xp_pauli import XPPauli + +# TODO from qiskit_qec.utils.pauli_rep import split_pauli, cpxstr2exp + +# from qiskit.quantum_info.operators.symplectic.pauli import _split_pauli_label, _phase_from_label + + +@ddt +class TestXPPauliInit(QiskitTestCase): + """Tests for XPPauli initialization.""" + + def test_array_init(self): + """Test array initialization.""" + # Test case taken from Mark's paper on XPF + matrix = np.array([1, 1, 1, 0, 0, 0, 0, 0, 0, 4, 0, 0, 0, 0], dtype=np.int64) + phase_exp = 12 + precision = 8 + xppauli = XPPauli(data=matrix, phase_exp=phase_exp, precision=precision) + np.testing.assert_equal(xppauli.matrix, np.atleast_2d(matrix)) + np.testing.assert_equal(xppauli._phase_exp, phase_exp) + np.testing.assert_equal(xppauli.precision, precision) + + +@ddt +class TestXPPauli(QiskitTestCase): + """Tests for XPPauli operator class.""" + + def test_precision_rescale(self): + """Test precision rescaling method.""" + # Test case taken from Mark's paper on XPF + matrix = np.array([1, 1, 1, 0, 0, 0, 0, 0, 0, 4, 0, 0, 0, 0], dtype=np.int64) + phase_exp = 12 + precision = 8 + new_precision = 2 + xppauli = XPPauli(data=matrix, phase_exp=phase_exp, precision=precision) + rescaled_xppauli = xppauli.rescale_precision(new_precision=new_precision) + target_matrix = np.array([1, 1, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0], dtype=np.int64) + target_phase_exp = 3 + target_xppauli = XPPauli( + data=target_matrix, phase_exp=target_phase_exp, precision=new_precision + ) + np.testing.assert_equal(target_xppauli.matrix, rescaled_xppauli.matrix) + np.testing.assert_equal(target_xppauli._phase_exp, rescaled_xppauli._phase_exp) + np.testing.assert_equal(target_xppauli.precision, rescaled_xppauli.precision) + + def test_weight(self): + """Test weight method.""" + matrix = np.array([1, 1, 1, 0, 0, 1, 0, 0, 3, 4, 0, 0, 0, 1], dtype=np.int64) + phase_exp = 12 + precision = 8 + xppauli = XPPauli(data=matrix, phase_exp=phase_exp, precision=precision) + value = xppauli.weight() + target = 5 + self.assertEqual(target, value) + + def test_diagonal(self): + """Test is_diagonal method.""" + # Test case taken from Mark's paper, Table 5. + matrix = np.array([0, 0, 0, 0, 0, 0, 0, 3, 3, 3, 3, 3, 3, 3], dtype=np.int64) + phase_exp = 0 + precision = 8 + xppauli = XPPauli(data=matrix, phase_exp=phase_exp, precision=precision) + value = xppauli.is_diagonal() + target = np.array([True]) + self.assertEqual(target, value) + + # Test case taken from Mark's paper, Table 5. + matrix = np.array([0, 0, 0, 0, 0, 0, 0, 0, 1, 3, 3, 3, 3, 3], dtype=np.int64) + phase_exp = 0 + precision = 8 + xppauli = XPPauli(data=matrix, phase_exp=phase_exp, precision=precision) + value = xppauli.is_diagonal() + target = np.array([True]) + self.assertEqual(target, value) + + matrix = np.array([0, 1, 0, 1, 0, 0, 0, 0, 1, 3, 3, 3, 3, 3], dtype=np.int64) + phase_exp = 12 + precision = 8 + xppauli = XPPauli(data=matrix, phase_exp=phase_exp, precision=precision) + value = xppauli.is_diagonal() + target = np.array([False]) + self.assertEqual(target, value) + + def test_antisymmetric_op(self): + """Test antisymmetric_op method.""" + + matrix = np.array([0, 0, 0, 0, 0, 0, 0, 0, 1, 2, 3, 3, 3, 3], dtype=np.int64) + phase_exp = 0 + precision = 8 + xppauli = XPPauli(data=matrix, phase_exp=phase_exp, precision=precision) + dinput = np.array(xppauli.z) + value = xppauli.antisymmetric_op(dinput) + + target_matrix = np.array([0, 0, 0, 0, 0, 0, 0, 0, -1, -2, -3, -3, -3, -3], dtype=np.int64) + target_phase_exp = 15 + target_precision = 8 + target = XPPauli(data=target_matrix, phase_exp=target_phase_exp, precision=target_precision) + np.testing.assert_equal(target.matrix, value.matrix) + np.testing.assert_equal(target._phase_exp, value._phase_exp) + np.testing.assert_equal(target.precision, value.precision) + + def test_inverse(self): + """Test inverse method.""" + + matrix = np.array([0, 0, 0, 1, 0, 1, 1, 5, 5, 6, 1, 1, 4, 0], dtype=np.int64) + phase_exp = 1 + precision = 8 + xppauli = XPPauli(data=matrix, phase_exp=phase_exp, precision=precision) + value = xppauli.inverse() + + target_matrix = np.array([0, 0, 0, 1, 0, 1, 1, 3, 3, 2, 1, 7, 4, 0], dtype=np.int64) + target_phase_exp = 5 + target_precision = 8 + target = XPPauli(data=target_matrix, phase_exp=target_phase_exp, precision=target_precision) + np.testing.assert_equal(target.matrix, value.matrix) + np.testing.assert_equal(target._phase_exp, value._phase_exp) + np.testing.assert_equal(target.precision, value.precision) + + def test_power(self): + """Test power method.""" + matrix = np.array([1, 1, 1, 0, 0, 1, 0, 0, 3, 4, 0, 0, 0, 1], dtype=np.int64) + phase_exp = 12 + precision = 8 + n = np.array([5]) + xppauli = XPPauli(data=matrix, phase_exp=phase_exp, precision=precision) + value = xppauli.power(n=n) + + target_matrix = np.array([1, 1, 1, 0, 0, 1, 0, 0, 3, 4, 0, 0, 0, 5], dtype=np.int64) + target_phase_exp = 8 + target_precision = 8 + target = XPPauli(data=target_matrix, phase_exp=target_phase_exp, precision=target_precision) + np.testing.assert_equal(target.matrix, value.matrix) + np.testing.assert_equal(target._phase_exp, value._phase_exp) + np.testing.assert_equal(target.precision, value.precision) + + def test_multiplication(self): + """Test compose method.""" + # Test case taken from Mark's code. + a_matrix = np.array([0, 1, 0, 0, 2, 0], dtype=np.int64) + a_phase_exp = 6 + a_precision = 4 + a = XPPauli(data=a_matrix, phase_exp=a_phase_exp, precision=a_precision) + b_matrix = np.array([1, 1, 1, 3, 3, 0], dtype=np.int64) + b_phase_exp = 2 + b_precision = 4 + b = XPPauli(data=b_matrix, phase_exp=b_phase_exp, precision=b_precision) + value = a.compose(b) + + target_matrix = np.array([1, 0, 1, 3, 3, 0], dtype=np.int64) + target_phase_exp = 6 + target_precision = 4 + target = XPPauli(data=target_matrix, phase_exp=target_phase_exp, precision=target_precision) + np.testing.assert_equal(target.matrix, value.matrix) + np.testing.assert_equal(target._phase_exp, value._phase_exp) + np.testing.assert_equal(target.precision, value.precision) + + def test_conjugate(self): + """Test conjugate method.""" + a_matrix = np.array([1, 0, 1, 1, 5, 3, 5, 4], dtype=np.int64) + a_phase_exp = 4 + a_precision = 6 + a = XPPauli(data=a_matrix, phase_exp=a_phase_exp, precision=a_precision) + b_matrix = np.array([1, 0, 0, 1, 4, 1, 0, 1], dtype=np.int64) + b_phase_exp = 11 + b_precision = 6 + b = XPPauli(data=b_matrix, phase_exp=b_phase_exp, precision=b_precision) + value_front = a.conjugate(b, front=True) + value_back = a.conjugate(b, front=False) + + target_matrix_front = np.array([1, 0, 0, 1, 0, 1, 0, 1], dtype=np.int64) + target_phase_exp_front = 3 + target_precision_front = 6 + target_front = XPPauli( + data=target_matrix_front, + phase_exp=target_phase_exp_front, + precision=target_precision_front, + ) + target_matrix_back = np.array([1, 0, 1, 1, 3, 3, 5, 4], dtype=np.int64) + target_phase_exp_back = 0 + target_precision_back = 6 + target_back = XPPauli( + data=target_matrix_back, + phase_exp=target_phase_exp_back, + precision=target_precision_back, + ) + np.testing.assert_equal(target_front.matrix, value_front.matrix) + np.testing.assert_equal(target_front._phase_exp, value_front._phase_exp) + np.testing.assert_equal(target_front.precision, value_front.precision) + np.testing.assert_equal(target_back.matrix, value_back.matrix) + np.testing.assert_equal(target_back._phase_exp, value_back._phase_exp) + np.testing.assert_equal(target_back.precision, value_back.precision) + + def test_commutator(self): + """Test commutator method.""" + a_matrix = np.array([1, 0, 1, 1, 5, 3, 5, 4], dtype=np.int64) + a_phase_exp = 4 + a_precision = 6 + a = XPPauli(data=a_matrix, phase_exp=a_phase_exp, precision=a_precision) + b_matrix = np.array([1, 0, 0, 1, 4, 1, 0, 1], dtype=np.int64) + b_phase_exp = 11 + b_precision = 6 + b = XPPauli(data=b_matrix, phase_exp=b_phase_exp, precision=b_precision) + value_front = a.commutator(b, front=True) + value_back = a.commutator(b, front=False) + + target_matrix_front = np.array([0, 0, 0, 0, 4, 0, 0, 0], dtype=np.int64) + target_phase_exp_front = 8 + target_precision_front = 6 + target_front = XPPauli( + data=target_matrix_front, + phase_exp=target_phase_exp_front, + precision=target_precision_front, + ) + target_matrix_back = np.array([0, 0, 0, 0, 2, 0, 0, 0], dtype=np.int64) + target_phase_exp_back = 4 + target_precision_back = 6 + target_back = XPPauli( + data=target_matrix_back, + phase_exp=target_phase_exp_back, + precision=target_precision_back, + ) + np.testing.assert_equal(target_front.matrix, value_front.matrix) + np.testing.assert_equal(target_front._phase_exp, value_front._phase_exp) + np.testing.assert_equal(target_front.precision, value_front.precision) + np.testing.assert_equal(target_back.matrix, value_back.matrix) + np.testing.assert_equal(target_back._phase_exp, value_back._phase_exp) + np.testing.assert_equal(target_back.precision, value_back.precision) + + def test_degree(self): + """Test degree method.""" + matrix = np.array([0, 0, 0, 2, 1, 0], dtype=np.int64) + phase_exp = 2 + precision = 4 + xppauli = XPPauli(data=matrix, phase_exp=phase_exp, precision=precision) + value = xppauli.degree() + + target = 4 + self.assertEqual(target, value) + + def test_fundamental_phase(self): + """Test fundamental_phase method.""" + matrix = np.array([1, 0, 1, 1, 5, 3, 5, 4], dtype=np.int64) + phase_exp = 4 + precision = 6 + xppauli = XPPauli(data=matrix, phase_exp=phase_exp, precision=precision) + value = xppauli.fundamental_phase() + + target = np.array([0]) + self.assertEqual(target, value) + + def test_reset_eigenvalue(self): + """Test reset_eigenvalue method.""" + matrix = np.array([1, 1, 0, 1, 0, 1, 0, 4], dtype=np.int64) + phase_exp = 4 + precision = 6 + xppauli = XPPauli(data=matrix, phase_exp=phase_exp, precision=precision) + value = xppauli.reset_eigenvalue() + + target_matrix = np.array([1, 1, 0, 1, 0, 1, 0, 4], dtype=np.int64) + target_phase_exp = 1 + target_precision = 6 + target = XPPauli(data=target_matrix, phase_exp=target_phase_exp, precision=target_precision) + np.testing.assert_equal(target.matrix, value.matrix) + np.testing.assert_equal(target._phase_exp, value._phase_exp) + np.testing.assert_equal(target.precision, value.precision) + + +if __name__ == "__main__": + unittest.main() diff --git a/tests/operators/test_xp_pauli_list.py b/tests/operators/test_xp_pauli_list.py new file mode 100644 index 00000000..6477b5cf --- /dev/null +++ b/tests/operators/test_xp_pauli_list.py @@ -0,0 +1,344 @@ +# This code is part of Qiskit. +# +# (C) Copyright IBM 2017, 2020. +# +# 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 XPPauliList class.""" + +import unittest +import numpy as np +from ddt import ddt +from qiskit.test import QiskitTestCase +from qiskit_qec.operators.xp_pauli_list import XPPauliList + + +class TestXPPauliListInit(QiskitTestCase): + """Tests for XPPauliList initialization.""" + + def test_array_init(self): + """Test array initialization.""" + # Matrix array initialization + precision = 8 + matrix1 = np.array([1, 1, 1, 0, 0, 0, 0, 0, 0, 4, 0, 0, 0, 0], dtype=np.int64) + phase_exp1 = 12 + matrix2 = np.array([1, 1, 1, 0, 0, 0, 0, 0, 0, 2, 3, 0, 0, 0], dtype=np.int64) + phase_exp2 = 2 + matrix = np.array([matrix1, matrix2]) + phase_exp = np.array([phase_exp1, phase_exp2]) + xppaulilist = XPPauliList(data=matrix, phase_exp=phase_exp, precision=precision) + np.testing.assert_equal(xppaulilist.matrix, matrix) + np.testing.assert_equal(xppaulilist._phase_exp, phase_exp) + np.testing.assert_equal(xppaulilist.precision, precision) + + +@ddt +class TestXPPauliListOperator(QiskitTestCase): + """Tests for XPPauliList base operator methods.""" + + def test_precision_rescale(self): + """Test precision rescaling method.""" + matrix1 = np.array([1, 1, 1, 0, 0, 0, 0, 0, 0, 4, 0, 0, 0, 0], dtype=np.int64) + phase_exp1 = 12 + matrix2 = np.array([1, 1, 0, 1, 0, 0, 0, 0, 0, 4, 0, 0, 0, 6], dtype=np.int64) + phase_exp2 = 8 + precision = 8 + new_precision = 4 + matrix = np.array([matrix1, matrix2]) + phase_exp = np.array([phase_exp1, phase_exp2]) + + xppaulilist = XPPauliList(data=matrix, phase_exp=phase_exp, precision=precision) + rescaled_xppaulilist = xppaulilist.rescale_precision(new_precision=new_precision) + target_matrix1 = np.array([1, 1, 1, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0], dtype=np.int64) + target_phase_exp1 = 6 + target_matrix2 = np.array([1, 1, 0, 1, 0, 0, 0, 0, 0, 2, 0, 0, 0, 3], dtype=np.int64) + target_phase_exp2 = 4 + target_matrix = np.array([target_matrix1, target_matrix2]) + target_phase_exp = np.array([target_phase_exp1, target_phase_exp2]) + target_xppaulilist = XPPauliList( + data=target_matrix, phase_exp=target_phase_exp, precision=new_precision + ) + np.testing.assert_equal(target_xppaulilist.matrix, rescaled_xppaulilist.matrix) + np.testing.assert_equal(target_xppaulilist._phase_exp, rescaled_xppaulilist._phase_exp) + np.testing.assert_equal(target_xppaulilist.precision, rescaled_xppaulilist.precision) + + def test_weight(self): + """Test weight method.""" + precision = 8 + matrix1 = np.array([1, 1, 1, 0, 0, 1, 0, 0, 3, 4, 0, 0, 0, 1], dtype=np.int64) + phase_exp1 = 12 + matrix2 = np.array([1, 1, 1, 0, 0, 0, 0, 0, 0, 2, 0, 0, 0, 4], dtype=np.int64) + phase_exp2 = 3 + matrix = np.array([matrix1, matrix2]) + phase_exp = np.array([phase_exp1, phase_exp2]) + xppaulilist = XPPauliList(data=matrix, phase_exp=phase_exp, precision=precision) + value = xppaulilist.weight() + target = np.array([5, 4]) + np.testing.assert_equal(target, value) + + def test_diagonal(self): + """Test is_diagonal method.""" + precision = 8 + matrix1 = np.array([0, 0, 0, 0, 0, 0, 0, 3, 3, 3, 3, 3, 3, 3], dtype=np.int64) + phase_exp1 = 0 + + matrix2 = np.array([0, 0, 0, 0, 0, 0, 0, 0, 1, 3, 3, 3, 3, 3], dtype=np.int64) + phase_exp2 = 0 + + matrix3 = np.array([0, 1, 0, 1, 0, 0, 0, 0, 1, 3, 3, 3, 3, 3], dtype=np.int64) + phase_exp3 = 12 + matrix = np.array([matrix1, matrix2, matrix3]) + phase_exp = np.array([phase_exp1, phase_exp2, phase_exp3]) + xppaulilist = XPPauliList(data=matrix, phase_exp=phase_exp, precision=precision) + + value = xppaulilist.is_diagonal() + target = np.array([True, True, False]) + np.testing.assert_equal(target, value) + + def test_antisymmetric_op(self): + """Test antisymmetric_op method.""" + + matrix = np.array( + [ + [0, 0, 0, 0, 0, 0, 0, 0, 1, 2, 3, 3, 3, 3], + [0, 0, 0, 0, 0, 0, 0, 3, 1, 2, 3, 7, 6, 3], + ], + dtype=np.int64, + ) + phase_exp = np.array([0, 0]) + precision = 8 + xppauli_list = XPPauliList(data=matrix, phase_exp=phase_exp, precision=precision) + dinput = xppauli_list.z + value = xppauli_list.antisymmetric_op(dinput) + + target_matrix = np.array( + [ + [0, 0, 0, 0, 0, 0, 0, 0, -1, -2, -3, -3, -3, -3], + [0, 0, 0, 0, 0, 0, 0, -3, -1, -2, -3, -7, -6, -3], + ], + dtype=np.int64, + ) + target_phase_exp = np.array([15, 25]) + target_precision = 8 + target = XPPauliList( + data=target_matrix, phase_exp=target_phase_exp, precision=target_precision + ) + np.testing.assert_equal(target.matrix, value.matrix) + np.testing.assert_equal(target._phase_exp, value._phase_exp) + np.testing.assert_equal(target.precision, value.precision) + + def test_inverse(self): + """Test inverse method.""" + + matrix = np.array( + [ + [1, 1, 0, 1, 1, 0, 1, 2, 4, 4, 3, 1, 6, 1], + [0, 1, 0, 0, 1, 0, 1, 7, 7, 3, 4, 6, 2, 7], + ], + dtype=np.int64, + ) + phase_exp = np.array([1, 0]) + precision = 8 + xppauli_list = XPPauliList(data=matrix, phase_exp=phase_exp, precision=precision) + value = xppauli_list.inverse() + + target_matrix = np.array( + [ + [1, 1, 0, 1, 1, 0, 1, 2, 4, 4, 3, 1, 2, 1], + [0, 1, 0, 0, 1, 0, 1, 1, 7, 5, 4, 6, 6, 7], + ], + dtype=np.int64, + ) + target_phase_exp = np.array([9, 8]) + target_precision = 8 + target = XPPauliList( + data=target_matrix, phase_exp=target_phase_exp, precision=target_precision + ) + np.testing.assert_equal(target.matrix, value.matrix) + np.testing.assert_equal(target._phase_exp, value._phase_exp) + np.testing.assert_equal(target.precision, value.precision) + + def test_power(self): + """Test power method.""" + matrix = np.array( + [ + [1, 1, 1, 0, 0, 1, 0, 0, 3, 4, 0, 0, 0, 1], + [1, 1, 1, 0, 0, 1, 0, 0, 3, 4, 0, 0, 0, 1], + ], + dtype=np.int64, + ) + phase_exp = np.array([12, 12]) + precision = 8 + n = np.array([5, 5]) + xppauli_list = XPPauliList(data=matrix, phase_exp=phase_exp, precision=precision) + value = xppauli_list.power(n=n) + + target_matrix = np.array( + [ + [1, 1, 1, 0, 0, 1, 0, 0, 3, 4, 0, 0, 0, 5], + [1, 1, 1, 0, 0, 1, 0, 0, 3, 4, 0, 0, 0, 5], + ], + dtype=np.int64, + ) + target_phase_exp = np.array([8, 8]) + target_precision = 8 + target = XPPauliList( + data=target_matrix, phase_exp=target_phase_exp, precision=target_precision + ) + np.testing.assert_equal(target.matrix, value.matrix) + np.testing.assert_equal(target._phase_exp, value._phase_exp) + np.testing.assert_equal(target.precision, value.precision) + + def test_multiplication(self): + """Test compose method.""" + a_matrix = np.array([[0, 1, 0, 0, 2, 0], [0, 1, 0, 0, 2, 0]], dtype=np.int64) + a_phase_exp = np.array([6, 6]) + a_precision = 4 + a = XPPauliList(data=a_matrix, phase_exp=a_phase_exp, precision=a_precision) + b_matrix = np.array([[1, 1, 1, 3, 3, 0], [1, 1, 1, 3, 3, 0]], dtype=np.int64) + b_phase_exp = np.array([2, 2]) + b_precision = 4 + b = XPPauliList(data=b_matrix, phase_exp=b_phase_exp, precision=b_precision) + value = a.compose(b) + + target_matrix = np.array([[1, 0, 1, 3, 3, 0], [1, 0, 1, 3, 3, 0]], dtype=np.int64) + target_phase_exp = np.array([6, 6]) + target_precision = 4 + target = XPPauliList( + data=target_matrix, phase_exp=target_phase_exp, precision=target_precision + ) + np.testing.assert_equal(target.matrix, value.matrix) + np.testing.assert_equal(target._phase_exp, value._phase_exp) + np.testing.assert_equal(target.precision, value.precision) + + def test_conjugate(self): + """Test conjugate method.""" + a_matrix = np.array([[1, 0, 1, 1, 5, 3, 5, 4], [1, 0, 1, 0, 1, 5, 2, 0]], dtype=np.int64) + a_phase_exp = np.array([4, 7]) + a_precision = 6 + a = XPPauliList(data=a_matrix, phase_exp=a_phase_exp, precision=a_precision) + b_matrix = np.array([[1, 0, 0, 1, 4, 1, 0, 1], [0, 1, 1, 0, 1, 3, 0, 5]], dtype=np.int64) + b_phase_exp = np.array([11, 2]) + b_precision = 6 + b = XPPauliList(data=b_matrix, phase_exp=b_phase_exp, precision=b_precision) + value_front = a.conjugate(b, front=True) + value_back = a.conjugate(b, front=False) + + target_matrix_front = np.array( + [[1, 0, 0, 1, 0, 1, 0, 1], [0, 1, 1, 0, 5, 5, 4, 5]], dtype=np.int64 + ) + target_phase_exp_front = np.array([3, 10]) + target_precision_front = 6 + target_front = XPPauliList( + data=target_matrix_front, + phase_exp=target_phase_exp_front, + precision=target_precision_front, + ) + target_matrix_back = np.array( + [[1, 0, 1, 1, 3, 3, 5, 4], [1, 0, 1, 0, 5, 1, 4, 0]], dtype=np.int64 + ) + target_phase_exp_back = np.array([0, 11]) + target_precision_back = 6 + target_back = XPPauliList( + data=target_matrix_back, + phase_exp=target_phase_exp_back, + precision=target_precision_back, + ) + np.testing.assert_equal(target_front.matrix, value_front.matrix) + np.testing.assert_equal(target_front._phase_exp, value_front._phase_exp) + np.testing.assert_equal(target_front.precision, value_front.precision) + np.testing.assert_equal(target_back.matrix, value_back.matrix) + np.testing.assert_equal(target_back._phase_exp, value_back._phase_exp) + np.testing.assert_equal(target_back.precision, value_back.precision) + + def test_commutator(self): + """Test commutator method.""" + a_matrix = np.array([[1, 0, 1, 1, 5, 3, 5, 4], [1, 0, 1, 0, 1, 5, 2, 0]], dtype=np.int64) + a_phase_exp = np.array([4, 7]) + a_precision = 6 + a = XPPauliList(data=a_matrix, phase_exp=a_phase_exp, precision=a_precision) + b_matrix = np.array([[1, 0, 0, 1, 4, 1, 0, 1], [0, 1, 1, 0, 1, 3, 0, 5]], dtype=np.int64) + b_phase_exp = np.array([11, 2]) + b_precision = 6 + b = XPPauliList(data=b_matrix, phase_exp=b_phase_exp, precision=b_precision) + value_front = a.commutator(b, front=True) + value_back = a.commutator(b, front=False) + + target_matrix_front = np.array( + [[0, 0, 0, 0, 4, 0, 0, 0], [0, 0, 0, 0, 4, 4, 2, 0]], dtype=np.int64 + ) + target_phase_exp_front = np.array([8, 8]) + target_precision_front = 6 + target_front = XPPauliList( + data=target_matrix_front, + phase_exp=target_phase_exp_front, + precision=target_precision_front, + ) + target_matrix_back = np.array( + [[0, 0, 0, 0, 2, 0, 0, 0], [0, 0, 0, 0, 2, 2, 4, 0]], dtype=np.int64 + ) + target_phase_exp_back = np.array([4, 4]) + target_precision_back = 6 + target_back = XPPauliList( + data=target_matrix_back, + phase_exp=target_phase_exp_back, + precision=target_precision_back, + ) + np.testing.assert_equal(target_front.matrix, value_front.matrix) + np.testing.assert_equal(target_front._phase_exp, value_front._phase_exp) + np.testing.assert_equal(target_front.precision, value_front.precision) + np.testing.assert_equal(target_back.matrix, value_back.matrix) + np.testing.assert_equal(target_back._phase_exp, value_back._phase_exp) + np.testing.assert_equal(target_back.precision, value_back.precision) + + def test_degree(self): + """Test degree method.""" + matrix = np.array([[1, 0, 1, 1, 5, 3, 5, 4], [1, 0, 1, 1, 5, 4, 1, 5]], dtype=np.int64) + phase_exp = np.array([4, 3]) + precision = 6 + xppauli_list = XPPauliList(data=matrix, phase_exp=phase_exp, precision=precision) + value = xppauli_list.degree() + + target = np.array([2, 6]) + np.testing.assert_equal(target, value) + + def test_fundamental_phase(self): + """Test fundamental_phase method.""" + matrix = np.array([[1, 0, 1, 1, 5, 3, 5, 4], [1, 0, 1, 1, 5, 4, 1, 5]], dtype=np.int64) + phase_exp = np.array([4, 3]) + precision = 6 + xppauli_list = XPPauliList(data=matrix, phase_exp=phase_exp, precision=precision) + value = xppauli_list.fundamental_phase() + + target = np.array([0, 0]) + np.testing.assert_equal(target, value) + + def test_reset_eigenvalue(self): + """Test reset_eigenvalue method.""" + matrix = np.array([[0, 0, 1, 1, 1, 0, 4, 2], [1, 1, 0, 1, 0, 1, 0, 4]], dtype=np.int64) + phase_exp = np.array([7, 4]) + precision = 6 + xppauli_list = XPPauliList(data=matrix, phase_exp=phase_exp, precision=precision) + value = xppauli_list.reset_eigenvalue() + + target_matrix = np.array( + [[0, 0, 1, 1, 1, 0, 4, 2], [1, 1, 0, 1, 0, 1, 0, 4]], dtype=np.int64 + ) + target_phase_exp = np.array([6, 1]) + target_precision = 6 + target = XPPauliList( + data=target_matrix, phase_exp=target_phase_exp, precision=target_precision + ) + np.testing.assert_equal(target.matrix, value.matrix) + np.testing.assert_equal(target._phase_exp, value._phase_exp) + np.testing.assert_equal(target.precision, value.precision) + + +if __name__ == "__main__": + unittest.main() diff --git a/tests/repetition_codes/test_codes.py b/tests/repetition_codes/test_codes.py index 2262150b..3d9d13ed 100644 --- a/tests/repetition_codes/test_codes.py +++ b/tests/repetition_codes/test_codes.py @@ -387,6 +387,10 @@ def test_transpilation(self): "Error: Wrong number of cx gates after transpilation.", ) + def test_empty_decoding_grapg(self): + """Test initializtion of decoding graphs with None""" + DecodingGraph(None) + if __name__ == "__main__": unittest.main() diff --git a/tests/utils/test_xp_pauli_rep.py b/tests/utils/test_xp_pauli_rep.py new file mode 100644 index 00000000..b8d01ca4 --- /dev/null +++ b/tests/utils/test_xp_pauli_rep.py @@ -0,0 +1,95 @@ +"""Test xp pauli rep.""" + +from unittest import TestCase +import numpy as np + +from qiskit_qec.utils.xp_pauli_rep import ( + xp_symplectic2str, + PRODUCT_SYNTAX, + LATEX_SYNTAX, + XP_SYMPLECTIC_SYNTAX, +) + + +class TestXPPauliRep(TestCase): + """Test xp pauli rep.""" + + def test_xp_symplectic2str(self): + """Tests xp_symplectic2str function.""" + + precision = 8 + matrix1 = np.array([1, 1, 1, 0, 0, 0, 0, 0, 0, 4, 0, 0, 0, 0], dtype=np.int64) + phase_exp1 = 12 + matrix2 = np.array([1, 1, 1, 0, 0, 0, 0, 0, 0, 2, 3, 0, 0, 0], dtype=np.int64) + phase_exp2 = 2 + matrix = np.array([matrix1, matrix2]) + phase_exp = np.array([phase_exp1, phase_exp2]) + + np.testing.assert_equal( + xp_symplectic2str(matrix, phase_exp, precision), + np.array(["XP8((w,12)(X(P,4))2(X)1(X)0)", "XP8((w,2)(P,3)3(X(P,2))2(X)1(X)0)"]), + ) + np.testing.assert_equal( + xp_symplectic2str(matrix, phase_exp, precision, qubit_order="left-to-right"), + np.array(["XP8((w,12)(X)0(X)1(X(P,4))2)", "XP8((w,2)(X)0(X)1(X(P,2))2(P,3)3)"]), + ) + np.testing.assert_equal( + xp_symplectic2str(matrix, phase_exp, precision, syntax=XP_SYMPLECTIC_SYNTAX), + np.array(["XP8(12|1 1 1 0 0 0 0|0 0 4 0 0 0 0)", "XP8(2|1 1 1 0 0 0 0|0 0 2 3 0 0 0)"]), + ) + np.testing.assert_equal( + xp_symplectic2str(matrix, phase_exp, precision, no_phase=True), + np.array(["XP8((X(P,4))2(X)1(X)0)", "XP8((P,3)3(X(P,2))2(X)1(X)0)"]), + ) + np.testing.assert_equal( + xp_symplectic2str(matrix, phase_exp, precision, syntax=PRODUCT_SYNTAX), + np.array( + ["XP8((w,12)(I)(I)(I)(I)(X(P,4))(X)(X))", "XP8((w,2)(I)(I)(I)(P,3)(X(P,2))(X)(X))"] + ), + ) + np.testing.assert_equal( + xp_symplectic2str(matrix, phase_exp, precision, syntax=LATEX_SYNTAX), + np.array( + [ + "XP_{8}((w,12)(XP^{4})_{2}(X)_{1}(X)_{0})", + "XP_{8}((w,2)(P^{3})_{3}(XP^{2})_{2}(X)_{1}(X)_{0})", + ] + ), + ) + + # Tests conversion to unique vector format for different formats + precision = 4 + matrix1 = np.array([1, 2, 3, 0, 0, 0, 0, 0, 0, 5, 0, 0, 0, 0], dtype=np.int64) + matrix2 = np.array([1, 2, 3, 0, 0, 0, 0, 0, 0, 2, 6, 0, 0, 0], dtype=np.int64) + matrix = np.array([matrix1, matrix2]) + + np.testing.assert_equal( + xp_symplectic2str(matrix, phase_exp, precision), + np.array(["XP4((w,4)(XP)2(X)0)", "XP4((w,2)(P,2)3(X(P,2))2(X)0)"]), + ) + np.testing.assert_equal( + xp_symplectic2str(matrix, phase_exp, precision, qubit_order="left-to-right"), + np.array(["XP4((w,4)(X)0(XP)2)", "XP4((w,2)(X)0(X(P,2))2(P,2)3)"]), + ) + np.testing.assert_equal( + xp_symplectic2str(matrix, phase_exp, precision, syntax=XP_SYMPLECTIC_SYNTAX), + np.array(["XP4(4|1 0 1 0 0 0 0|0 0 1 0 0 0 0)", "XP4(2|1 0 1 0 0 0 0|0 0 2 2 0 0 0)"]), + ) + np.testing.assert_equal( + xp_symplectic2str(matrix, phase_exp, precision, no_phase=True), + np.array(["XP4((XP)2(X)0)", "XP4((P,2)3(X(P,2))2(X)0)"]), + ) + np.testing.assert_equal( + xp_symplectic2str(matrix, phase_exp, precision, syntax=PRODUCT_SYNTAX), + np.array( + ["XP4((w,4)(I)(I)(I)(I)(XP)(I)(X))", "XP4((w,2)(I)(I)(I)(P,2)(X(P,2))(I)(X))"] + ), + ) + np.testing.assert_equal( + xp_symplectic2str(matrix, phase_exp, precision, syntax=LATEX_SYNTAX), + np.array( + ["XP_{4}((w,4)(XP)_{2}(X)_{0})", "XP_{4}((w,2)(P^{2})_{3}(XP^{2})_{2}(X)_{0})"] + ), + ) + + # TODO add scalar return test (single operator, same_type=True/False)