# 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. """Piecewise-polynomially-controlled Pauli rotations.""" from __future__ import annotations import warnings from typing import List, Optional import numpy as np from qiskit.circuit import QuantumRegister, AncillaRegister, QuantumCircuit, Gate from qiskit.circuit.exceptions import CircuitError from qiskit.utils.deprecation import deprecate_func from .functional_pauli_rotations import FunctionalPauliRotations from .polynomial_pauli_rotations import PolynomialPauliRotations, PolynomialPauliRotationsGate from .integer_comparator import IntegerComparator, IntegerComparatorGate class PiecewisePolynomialPauliRotations(FunctionalPauliRotations): r"""Piecewise-polynomially-controlled Pauli rotations. This class implements a piecewise polynomial (not necessarily continuous) function, :math:`f(x)`, on qubit amplitudes, which is defined through breakpoints and coefficients as follows. Suppose the breakpoints :math:`(x_0, ..., x_J)` are a subset of :math:`[0, 2^n-1]`, where :math:`n` is the number of state qubits. Further on, denote the corresponding coefficients by :math:`[a_{j,1},...,a_{j,d}]`, where :math:`d` is the highest degree among all polynomials. Then :math:`f(x)` is defined as: .. math:: f(x) = \begin{cases} 0, x < x_0 \\ \sum_{i=0}^{i=d}a_{j,i}/2 x^i, x_j \leq x < x_{j+1} \end{cases} where if given the same number of breakpoints as polynomials, we implicitly assume :math:`x_{J+1} = 2^n`. .. note:: Note the :math:`1/2` factor in the coefficients of :math:`f(x)`, this is consistent with Qiskit's Pauli rotations. Examples: >>> from qiskit import QuantumCircuit >>> from qiskit.circuit.library.arithmetic.piecewise_polynomial_pauli_rotations import\ ... PiecewisePolynomialPauliRotations >>> qubits, breakpoints, coeffs = (2, [0, 2], [[0, -1.2],[-1, 1, 3]]) >>> poly_r = PiecewisePolynomialPauliRotations(num_state_qubits=qubits, ...breakpoints=breakpoints, coeffs=coeffs) >>> >>> qc = QuantumCircuit(poly_r.num_qubits) >>> qc.h(list(range(qubits))); >>> qc.append(poly_r.to_instruction(), list(range(qc.num_qubits))); >>> qc.draw() ┌───┐┌──────────┐ q_0: ┤ H ├┤0 ├ ├───┤│ │ q_1: ┤ H ├┤1 ├ └───┘│ │ q_2: ─────┤2 ├ │ pw_poly │ q_3: ─────┤3 ├ │ │ q_4: ─────┤4 ├ │ │ q_5: ─────┤5 ├ └──────────┘ References: [1] Haener, T., Roetteler, M., & Svore, K. M. (2018). Optimizing Quantum Circuits for Arithmetic. `arXiv:1805.12445 `_ [2] Carrera Vazquez, A., Hiptmair, R., & Woerner, S. (2022). Enhancing the Quantum Linear Systems Algorithm using Richardson Extrapolation. `ACM Transactions on Quantum Computing 3, 1, Article 2 `_ """ @deprecate_func( since="2.2", additional_msg="Use the class PiecewisePolynomialPauliRotationsGate instead.", removal_timeline="in Qiskit 3.0", ) def __init__( self, num_state_qubits: Optional[int] = None, breakpoints: Optional[List[int]] = None, coeffs: Optional[List[List[float]]] = None, basis: str = "Y", name: str = "pw_poly", ) -> None: """ Args: num_state_qubits: The number of qubits representing the state. breakpoints: The breakpoints to define the piecewise-linear function. Defaults to ``[0]``. coeffs: The coefficients of the polynomials for different segments of the piecewise-linear function. ``coeffs[j][i]`` is the coefficient of the i-th power of x for the j-th polynomial. Defaults to linear: ``[[1]]``. basis: The type of Pauli rotation (``'X'``, ``'Y'``, ``'Z'``). name: The name of the circuit. """ # store parameters self._breakpoints = breakpoints if breakpoints is not None else [0] self._coeffs = coeffs if coeffs is not None else [[1]] # store a list of coefficients as homogeneous polynomials adding 0's where necessary self._hom_coeffs = [] self._degree = len(max(self._coeffs, key=len)) - 1 for poly in self._coeffs: self._hom_coeffs.append(poly + [0] * (self._degree + 1 - len(poly))) super().__init__(num_state_qubits=num_state_qubits, basis=basis, name=name) @property def breakpoints(self) -> List[int]: """The breakpoints of the piecewise polynomial function. The function is polynomial in the intervals ``[point_i, point_{i+1}]`` where the last point implicitly is ``2**(num_state_qubits + 1)``. Returns: The list of breakpoints. """ if ( self.num_state_qubits is not None and len(self._breakpoints) == len(self.coeffs) and self._breakpoints[-1] < 2**self.num_state_qubits ): return self._breakpoints + [2**self.num_state_qubits] return self._breakpoints @breakpoints.setter def breakpoints(self, breakpoints: List[int]) -> None: """Set the breakpoints. Args: breakpoints: The new breakpoints. """ self._invalidate() self._breakpoints = breakpoints if self.num_state_qubits and breakpoints: self._reset_registers(self.num_state_qubits) @property def coeffs(self) -> List[List[float]]: """The coefficients of the polynomials. Returns: The polynomial coefficients per interval as nested lists. """ return self._coeffs @coeffs.setter def coeffs(self, coeffs: List[List[float]]) -> None: """Set the polynomials. Args: coeffs: The new polynomials. """ self._invalidate() self._coeffs = coeffs # update the homogeneous polynomials and degree self._hom_coeffs = [] self._degree = len(max(self._coeffs, key=len)) - 1 for poly in self._coeffs: self._hom_coeffs.append(poly + [0] * (self._degree + 1 - len(poly))) if self.num_state_qubits and coeffs: self._reset_registers(self.num_state_qubits) @property def mapped_coeffs(self) -> List[List[float]]: """The coefficients mapped to the internal representation, since we only compare x>=breakpoint. Returns: The mapped coefficients. """ return _map_coeffs(self._hom_coeffs) @property def contains_zero_breakpoint(self) -> bool | np.bool_: """Whether 0 is the first breakpoint. Returns: True, if 0 is the first breakpoint, otherwise False. """ return np.isclose(0, self.breakpoints[0]) def evaluate(self, x: float) -> float: """Classically evaluate the piecewise polynomial rotation. Args: x: Value to be evaluated at. Returns: Value of piecewise polynomial function at x. """ y = 0 for i, breakpt in enumerate(self.breakpoints): y = y + (x >= breakpt) * (np.poly1d(self.mapped_coeffs[i][::-1])(x)) return y def _check_configuration(self, raise_on_failure: bool = True) -> bool: """Check if the current configuration is valid.""" valid = True if self.num_state_qubits is None: valid = False if raise_on_failure: raise AttributeError("The number of qubits has not been set.") if self.num_qubits < self.num_state_qubits + 1: valid = False if raise_on_failure: raise CircuitError( "Not enough qubits in the circuit, need at least " f"{self.num_state_qubits + 1}." ) if len(self.breakpoints) != len(self.coeffs) + 1: valid = False if raise_on_failure: raise ValueError("Mismatching number of breakpoints and polynomials.") return valid def _reset_registers(self, num_state_qubits: Optional[int]) -> None: """Reset the registers.""" self.qregs = [] if num_state_qubits: qr_state = QuantumRegister(num_state_qubits) qr_target = QuantumRegister(1) self.qregs = [qr_state, qr_target] # Calculate number of ancilla qubits required num_ancillas = num_state_qubits + 1 if self.contains_zero_breakpoint: num_ancillas -= 1 if num_ancillas > 0: qr_ancilla = AncillaRegister(num_ancillas) self.add_register(qr_ancilla) def _build(self): """If not already built, build the circuit.""" if self._is_built: return super()._build() circuit = QuantumCircuit(*self.qregs, name=self.name) qr_state = circuit.qubits[: self.num_state_qubits] qr_target = [circuit.qubits[self.num_state_qubits]] # Ancilla for the comparator circuit qr_ancilla = circuit.qubits[self.num_state_qubits + 1 :] # apply comparators and controlled linear rotations for i, point in enumerate(self.breakpoints[:-1]): if i == 0 and self.contains_zero_breakpoint: # the class itself is deprecated, no need to raise additional warnings during runtime with warnings.catch_warnings(): warnings.filterwarnings("ignore", category=DeprecationWarning, module="qiskit") poly_r = PolynomialPauliRotations( num_state_qubits=self.num_state_qubits, coeffs=self.mapped_coeffs[i], basis=self.basis, ) circuit.append(poly_r.to_gate(), qr_state[:] + qr_target) else: # the class itself is deprecated, no need to raise additional warnings during runtime with warnings.catch_warnings(): warnings.filterwarnings("ignore", category=DeprecationWarning, module="qiskit") comp = IntegerComparator(num_state_qubits=self.num_state_qubits, value=point) qr_state_full = qr_state[:] + [qr_ancilla[0]] # add compare qubit qr_remaining_ancilla = qr_ancilla[1:] # take remaining ancillas circuit.append( comp.to_gate(), qr_state_full[:] + qr_remaining_ancilla[: comp.num_ancillas] ) # the class itself is deprecated, no need to raise additional warnings during runtime with warnings.catch_warnings(): warnings.filterwarnings("ignore", category=DeprecationWarning, module="qiskit") poly_r = PolynomialPauliRotations( num_state_qubits=self.num_state_qubits, coeffs=self.mapped_coeffs[i], basis=self.basis, ) circuit.append( poly_r.to_gate().control(), [qr_ancilla[0]] + qr_state[:] + qr_target ) # uncompute comparator circuit.append( comp.to_gate().inverse(), qr_state_full[:] + qr_remaining_ancilla[: comp.num_ancillas], ) self.append(circuit.to_gate(), self.qubits) class PiecewisePolynomialPauliRotationsGate(Gate): r"""Piecewise-polynomially-controlled Pauli rotations. This class implements a piecewise polynomial (not necessarily continuous) function, :math:`f(x)`, on qubit amplitudes, which is defined through breakpoints and coefficients as follows. Suppose the breakpoints :math:`(x_0, ..., x_J)` are a subset of :math:`[0, 2^n-1]`, where :math:`n` is the number of state qubits. Further on, denote the corresponding coefficients by :math:`[a_{j,1},...,a_{j,d}]`, where :math:`d` is the highest degree among all polynomials. Then :math:`f(x)` is defined as: .. math:: f(x) = \begin{cases} 0, x < x_0 \\ \sum_{i=0}^{i=d}a_{j,i}/2 x^i, x_j \leq x < x_{j+1} \end{cases} where if given the same number of breakpoints as polynomials, we implicitly assume :math:`x_{J+1} = 2^n`. .. note:: Note the :math:`1/2` factor in the coefficients of :math:`f(x)`, this is consistent with Qiskit's Pauli rotations. Examples: >>> from qiskit import QuantumCircuit >>> from qiskit.circuit.library.arithmetic.piecewise_polynomial_pauli_rotations import\ ... PiecewisePolynomialPauliRotationsGate >>> qubits, breakpoints, coeffs = (2, [0, 2], [[0, -1.2],[-1, 1, 3]]) >>> poly_r = PiecewisePolynomialPauliRotationsGate(num_state_qubits=qubits, ...breakpoints=breakpoints, coeffs=coeffs) >>> >>> qc = QuantumCircuit(poly_r.num_qubits) >>> qc.h(list(range(qubits))); >>> qc.append(poly_r, list(range(qc.num_qubits))); >>> qc.draw() ┌───┐┌──────────┐ q_0: ┤ H ├┤0 ├ ├───┤│ │ q_1: ┤ H ├┤1 ├ └───┘│ │ q_2: ─────┤2 ├ │ pw_poly │ q_3: ─────┤3 ├ │ │ q_4: ─────┤4 ├ │ │ q_5: ─────┤5 ├ └──────────┘ References: [1] Haener, T., Roetteler, M., & Svore, K. M. (2018). Optimizing Quantum Circuits for Arithmetic. `arXiv:1805.12445 `_ [2] Carrera Vazquez, A., Hiptmair, R., & Woerner, S. (2022). Enhancing the Quantum Linear Systems Algorithm using Richardson Extrapolation. `ACM Transactions on Quantum Computing 3, 1, Article 2 `_ """ def __init__( self, num_state_qubits: int, breakpoints: list[int] | None = None, coeffs: list[list[float]] | None = None, basis: str = "Y", label: str | None = None, ) -> None: """ Args: num_state_qubits: The number of qubits representing the state. breakpoints: The breakpoints to define the piecewise-linear function. Defaults to ``[0]``. coeffs: The coefficients of the polynomials for different segments of the piecewise-linear function. ``coeffs[j][i]`` is the coefficient of the i-th power of x for the j-th polynomial. Defaults to linear: ``[[1]]``. basis: The type of Pauli rotation (``'X'``, ``'Y'``, ``'Z'``). label: An optional label for the gate. """ if coeffs is None: degree = 0 coeffs = [[1]] else: # store a list of coefficients as homogeneous polynomials adding 0's where necessary degree = len(max(coeffs, key=len)) - 1 coeffs = [poly + [0] * (degree + 1 - len(poly)) for poly in coeffs] # ensure the breakpoint contains 2 ** num_state_qubits breakpoints = breakpoints.copy() if breakpoints is None: breakpoints = [0, 2**num_state_qubits] elif breakpoints[-1] < 2**num_state_qubits: breakpoints.append(2**num_state_qubits) self.coeffs = coeffs self.breakpoints = breakpoints self.basis = basis self.num_compare = int(len(self.breakpoints) > 2) super().__init__( "PiecewisePolyPauli", num_state_qubits + self.num_compare + 1, [], label=label ) def evaluate(self, x: float) -> float: """Classically evaluate the piecewise polynomial rotation. Args: x: Value to be evaluated at. Returns: Value of piecewise polynomial function at x. """ mapped_coeffs = _map_coeffs(self.coeffs) y = 0 for i, breakpt in enumerate(self.breakpoints): y = y + (x >= breakpt) * (np.poly1d(mapped_coeffs[i][::-1])(x)) return y def _define(self): num_state_qubits = self.num_qubits - self.num_compare - 1 circuit = QuantumCircuit(self.num_qubits, name=self.name) qr_state = circuit.qubits[:num_state_qubits] if len(self.breakpoints) > 2: qr_target = [circuit.qubits[-2]] qr_compare = [circuit.qubits[-1]] else: qr_target = [circuit.qubits[-1]] qr_compare = [] # apply comparators and controlled linear rotations contains_zero_breakpoint = np.isclose(self.breakpoints[0], 0) mapped_coeffs = _map_coeffs(self.coeffs) for i, point in enumerate(self.breakpoints[:-1]): if i == 0 and contains_zero_breakpoint: # apply rotation poly_r = PolynomialPauliRotationsGate( num_state_qubits=num_state_qubits, coeffs=mapped_coeffs[i], basis=self.basis, ) circuit.append(poly_r, qr_state[:] + qr_target) else: # apply Comparator comp = IntegerComparatorGate(num_state_qubits=num_state_qubits, value=point) qr_state_full = qr_state[:] + qr_compare # add compare qubit circuit.append(comp, qr_state_full[:]) # apply controlled rotation poly_r = PolynomialPauliRotationsGate( num_state_qubits=num_state_qubits, coeffs=mapped_coeffs[i], basis=self.basis, ) circuit.append(poly_r.control(), qr_compare + qr_state[:] + qr_target) # uncompute comparator circuit.append(comp, qr_state_full[:]) self.definition = circuit def _map_coeffs(coeffs): mapped_coeffs = [] mapped_coeffs.append(coeffs[0]) degree = len(coeffs[0]) - 1 # all coeffs should have the same length by now for i in range(1, len(coeffs)): mapped_coeffs.append([]) for j in range(0, degree + 1): mapped_coeffs[i].append(coeffs[i][j] - coeffs[i - 1][j]) return mapped_coeffs