Update __array__ methods for Numpy 2.0 compatibility

As of Numpy 2.0, implementers of `__array__` are expected and required
to have a signature

    def __array__(self, dtype=None, copy=None): ...

In Numpys before 2.0, the `copy` argument will never be passed, and the
expected signature was

    def __array__(self, dtype=None): ...

Because of this, we have latitude to set `copy` in our implementations
to anything we like if we're running against Numpy 1.x, but we should
default to `copy=None` if we're running against Numpy 2.0.

The semantics of the `copy` argument to `np.array` changed in Numpy 2.0.
Now, `copy=False` means "raise a `ValueError` if a copy is required" and
`copy=None` means "copy only if required".  In Numpy 1.x, `copy=False`
meant "copy only if required".  In _both_ Numpy 1.x and 2.0,
`ndarray.astype` takes a `copy` argument, and in both, `copy=False`
means "copy only if required".  In Numpy 2.0 only, `np.asarray` gained a
`copy` argument with the same semantics as the `np.array` copy argument
from Numpy 2.0.

Further, the semantics of the `__array__` method changed in Numpy 2.0,
particularly around copying.  Now, Numpy will assume that it can pass
`copy=True` and the implementer will handle this.  If `copy=False` is
given and a copy or calculation is required, then the implementer is
required to raise `ValueError`.  We have a few places where the
`__array__` method may (or always does) calculate the array, so in all
these, we must forbid `copy=False`.

With all this in mind: this PR sets up all our implementers of
`__array__` to either default to `copy=None` if they will never actually
need to _use_ the `copy` argument within themselves (except perhaps to
test if it was set by Numpy 2.0 to `False`, as Numpy 1.x will never set
it), or to a compatibility shim `_numpy_compat.COPY_ONLY_IF_NEEDED` if
they do naturally want to use it with those semantics.  The pattern

    def __array__(self, dtype=None, copy=_numpy_compat.COPY_ONLY_IF_NEEDED):
        dtype = self._array.dtype if dtype is None else dtype
        return np.array(self._array, dtype=dtype, copy=copy)

using `array` instead of `asarray` lets us achieve all the desired
behaviour between the interactions of `dtype` and `copy` in a way that
is compatible with both Numpy 1.x and 2.x.

qiskit/__init__.py

@@ -53,6 +53,7 @@
 
 
 import qiskit._accelerate
+import qiskit._numpy_compat
 
 
 # Globally define compiled submodules. The normal import mechanism will not find compiled submodules

qiskit/_numpy_compat.py

@@ -0,0 +1,73 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2024.
+#
+# 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.
+
+"""Compatiblity helpers for the Numpy 1.x to 2.0 transition."""
+
+import re
+import typing
+import warnings
+
+import numpy as np
+
+# This version pattern is taken from the pypa packaging project:
+# https://github.com/pypa/packaging/blob/21.3/packaging/version.py#L223-L254 which is dual licensed
+# Apache 2.0 and BSD see the source for the original authors and other details.
+_VERSION_PATTERN = r"""
+    v?
+    (?:
+        (?:(?P<epoch>[0-9]+)!)?                           # epoch
+        (?P<release>[0-9]+(?:\.[0-9]+)*)                  # release segment
+        (?P<pre>                                          # pre-release
+            [-_\.]?
+            (?P<pre_l>(a|b|c|rc|alpha|beta|pre|preview))
+            [-_\.]?
+            (?P<pre_n>[0-9]+)?
+        )?
+        (?P<post>                                         # post release
+            (?:-(?P<post_n1>[0-9]+))
+            |
+            (?:
+                [-_\.]?
+                (?P<post_l>post|rev|r)
+                [-_\.]?
+                (?P<post_n2>[0-9]+)?
+            )
+        )?
+        (?P<dev>                                          # dev release
+            [-_\.]?
+            (?P<dev_l>dev)
+            [-_\.]?
+            (?P<dev_n>[0-9]+)?
+        )?
+    )
+    (?:\+(?P<local>[a-z0-9]+(?:[-_\.][a-z0-9]+)*))?       # local version
+"""
+
+VERSION = np.lib.NumpyVersion(np.__version__)
+VERSION_PARTS: typing.Tuple[int, ...]
+"""The numeric parts of the Numpy release version, e.g. ``(2, 0, 0)``.  Does not include pre- or
+post-release markers (e.g. ``rc1``)."""
+if match := re.fullmatch(_VERSION_PATTERN, np.__version__, flags=re.VERBOSE | re.IGNORECASE):
+    # Assuming Numpy won't ever introduce epochs, and we don't care about pre/post markers.
+    VERSION_PARTS = tuple(int(x) for x in match["release"].split("."))
+else:
+    # Just guess a version.  We know all existing Numpys have good version strings, so the only way
+    # this should trigger is from a new or a dev version.
+    warnings.warn(
+        f"Unrecognized version string for Numpy: '{np.__version__}'.  Assuming Numpy 2.0.",
+        RuntimeWarning,
+    )
+    VERSION_PARTS = (2, 0, 0)
+
+COPY_ONLY_IF_NEEDED = None if VERSION_PARTS >= (2, 0, 0) else False
+"""The sentinel value given to ``np.array`` and ``np.ndarray.astype`` (etc) to indicate that a copy
+should be made only if required."""

qiskit/circuit/__init__.py

@@ -833,7 +833,7 @@
 ``__array__``.  This is used by :meth:`Gate.to_matrix`, and has the signature:
 
 .. currentmodule:: None
-.. py:method:: __array__(dtype=None)
+.. py:method:: __array__(dtype=None, copy=None)
 
     Return a Numpy array representing the gate.  This can use the gate's :attr:`~Instruction.params`
     field, and may assume that these are numeric values (assuming the subclass expects that) and not
@@ -875,7 +875,9 @@ def power(self, exponent: float):
             # Also we have an efficient representation of power.
             return RXZGate(exponent * self.params[0])
 
-        def __array__(self, dtype=None):
+        def __array__(self, dtype=None, copy=None):
+            if copy is False:
+                raise ValueError("cannot produce a matrix without calculation")
             cos = math.cos(0.5 * self.params[0])
             isin = 1j * math.sin(0.5 * self.params[0])
             return np.array([
@@ -1340,6 +1342,7 @@ def __array__(self, dtype=None):
 """
 
 from .exceptions import CircuitError
+from . import _utils
 from .quantumcircuit import QuantumCircuit
 from .classicalregister import ClassicalRegister, Clbit
 from .quantumregister import QuantumRegister, Qubit, AncillaRegister, AncillaQubit

qiskit/circuit/_utils.py

@@ -15,6 +15,8 @@
 
 import math
 import numpy
+
+from qiskit import _numpy_compat
 from qiskit.exceptions import QiskitError
 from qiskit.circuit.exceptions import CircuitError
 from .parametervector import ParameterVectorElement
@@ -117,8 +119,9 @@ def with_gate_array(base_array):
     nonwritable = numpy.array(base_array, dtype=numpy.complex128)
     nonwritable.setflags(write=False)
 
-    def __array__(_self, dtype=None):
-        return numpy.asarray(nonwritable, dtype=dtype)
+    def __array__(_self, dtype=None, copy=_numpy_compat.COPY_ONLY_IF_NEEDED):
+        dtype = nonwritable.dtype if dtype is None else dtype
+        return numpy.array(nonwritable, dtype=dtype, copy=copy)
 
     def decorator(cls):
         if hasattr(cls, "__array__"):
@@ -149,15 +152,21 @@ def matrix_for_control_state(state):
     if cached_states is None:
         nonwritables = [matrix_for_control_state(state) for state in range(2**num_ctrl_qubits)]
 
-        def __array__(self, dtype=None):
-            return numpy.asarray(nonwritables[self.ctrl_state], dtype=dtype)
+        def __array__(self, dtype=None, copy=_numpy_compat.COPY_ONLY_IF_NEEDED):
+            arr = nonwritables[self.ctrl_state]
+            dtype = arr.dtype if dtype is None else dtype
+            return numpy.array(arr, dtype=dtype, copy=copy)
 
     else:
         nonwritables = {state: matrix_for_control_state(state) for state in cached_states}
 
-        def __array__(self, dtype=None):
-            if (out := nonwritables.get(self.ctrl_state)) is not None:
-                return numpy.asarray(out, dtype=dtype)
+        def __array__(self, dtype=None, copy=_numpy_compat.COPY_ONLY_IF_NEEDED):
+            if (arr := nonwritables.get(self.ctrl_state)) is not None:
+                dtype = arr.dtype if dtype is None else dtype
+                return numpy.array(arr, dtype=dtype, copy=copy)
+
+            if copy is False and copy is not _numpy_compat.COPY_ONLY_IF_NEEDED:
+                raise ValueError("could not produce matrix without calculation")
             return numpy.asarray(
                 _compute_control_matrix(base, num_ctrl_qubits, self.ctrl_state), dtype=dtype
             )

qiskit/circuit/delay.py

@@ -17,9 +17,11 @@
 from qiskit.circuit.exceptions import CircuitError
 from qiskit.circuit.instruction import Instruction
 from qiskit.circuit.gate import Gate
+from qiskit.circuit import _utils
 from qiskit.circuit.parameterexpression import ParameterExpression
 
 
+@_utils.with_gate_array(np.eye(2, dtype=complex))
 class Delay(Instruction):
     """Do nothing and just delay/wait/idle for a specified duration."""
 
@@ -53,10 +55,6 @@ def duration(self, duration):
         """Set the duration of this delay."""
         self.params = [duration]
 
-    def __array__(self, dtype=None):
-        """Return the identity matrix."""
-        return np.array([[1, 0], [0, 1]], dtype=dtype)
-
     def to_matrix(self) -> np.ndarray:
         """Return a Numpy.array for the unitary matrix. This has been
         added to enable simulation without making delay a full Gate type.

qiskit/circuit/library/generalized_gates/pauli.py

@@ -63,13 +63,13 @@ def inverse(self, annotated: bool = False):
         r"""Return inverted pauli gate (itself)."""
         return PauliGate(self.params[0])  # self-inverse
 
-    def __array__(self, dtype=None):
+    def __array__(self, dtype=None, copy=None):
         """Return a Numpy.array for the pauli gate.
         i.e. tensor product of the paulis"""
         # pylint: disable=cyclic-import
         from qiskit.quantum_info.operators import Pauli
 
-        return Pauli(self.params[0]).__array__(dtype=dtype)
+        return Pauli(self.params[0]).__array__(dtype=dtype, copy=copy)
 
     def validate_parameter(self, parameter):
         if isinstance(parameter, str):

qiskit/circuit/library/generalized_gates/permutation.py

@@ -147,8 +147,11 @@ def __init__(
 
         super().__init__(name="permutation", num_qubits=num_qubits, params=[pattern])
 
-    def __array__(self, dtype=None):
+    def __array__(self, dtype=None, copy=None):
         """Return a numpy.array for the Permutation gate."""
+        if copy is False:
+            raise ValueError("cannot produce matrix without calculation")
+
         nq = len(self.pattern)
         mat = np.zeros((2**nq, 2**nq), dtype=dtype)
 

qiskit/circuit/library/generalized_gates/unitary.py

@@ -18,6 +18,7 @@
 import typing
 import numpy
 
+from qiskit import _numpy_compat
 from qiskit.circuit.gate import Gate
 from qiskit.circuit.controlledgate import ControlledGate
 from qiskit.circuit.annotated_operation import AnnotatedOperation, ControlModifier
@@ -118,10 +119,10 @@ def __eq__(self, other):
             return False
         return matrix_equal(self.params[0], other.params[0])
 
-    def __array__(self, dtype=None):
+    def __array__(self, dtype=None, copy=_numpy_compat.COPY_ONLY_IF_NEEDED):
         """Return matrix for the unitary."""
-        # pylint: disable=unused-argument
-        return self.params[0]
+        dtype = self.params[0].dtype if dtype is None else dtype
+        return numpy.array(self.params[0], dtype=dtype, copy=copy)
 
     def inverse(self, annotated: bool = False):
         """Return the adjoint of the unitary."""

qiskit/circuit/library/hamiltonian_gate.py

@@ -21,6 +21,7 @@
 from numbers import Number
 import numpy as np
 
+from qiskit import _numpy_compat
 from qiskit.circuit.gate import Gate
 from qiskit.circuit.quantumcircuit import QuantumCircuit
 from qiskit.circuit.quantumregister import QuantumRegister
@@ -92,18 +93,22 @@ def __eq__(self, other):
         times_eq = self.params[1] == other.params[1]
         return operators_eq and times_eq
 
-    def __array__(self, dtype=None):
+    def __array__(self, dtype=None, copy=None):
         """Return matrix for the unitary."""
-        # pylint: disable=unused-argument
         import scipy.linalg
 
+        if copy is False:
+            raise ValueError("cannot produce matrix without calculation")
         try:
-            return scipy.linalg.expm(-1j * self.params[0] * float(self.params[1]))
+            time = float(self.params[1])
         except TypeError as ex:
             raise TypeError(
                 "Unable to generate Unitary matrix for "
                 "unbound t parameter {}".format(self.params[1])
             ) from ex
+        arr = scipy.linalg.expm(-1j * self.params[0] * time)
+        dtype = complex if dtype is None else dtype
+        return np.array(arr, dtype=dtype, copy=_numpy_compat.COPY_ONLY_IF_NEEDED)
 
     def inverse(self, annotated: bool = False):
         """Return the adjoint of the unitary."""

qiskit/circuit/library/standard_gates/global_phase.py

@@ -69,10 +69,12 @@ def inverse(self, annotated: bool = False):
         """
         return GlobalPhaseGate(-self.params[0])
 
-    def __array__(self, dtype=complex):
+    def __array__(self, dtype=None, copy=None):
         """Return a numpy.array for the global_phase gate."""
+        if copy is False:
+            raise ValueError("cannot produce matrix without calculation")
         theta = self.params[0]
-        return numpy.array([[numpy.exp(1j * theta)]], dtype=dtype)
+        return numpy.array([[numpy.exp(1j * theta)]], dtype=dtype or complex)
 
     def __eq__(self, other):
         if isinstance(other, GlobalPhaseGate):

qiskit/circuit/library/standard_gates/p.py

@@ -140,8 +140,10 @@ def inverse(self, annotated: bool = False):
         """
         return PhaseGate(-self.params[0])
 
-    def __array__(self, dtype=None):
+    def __array__(self, dtype=None, copy=None):
         """Return a numpy.array for the Phase gate."""
+        if copy is False:
+            raise ValueError("cannot produce matrix without calculation")
         lam = float(self.params[0])
         return numpy.array([[1, 0], [0, exp(1j * lam)]], dtype=dtype)
 
@@ -280,8 +282,10 @@ def inverse(self, annotated: bool = False):
         r"""Return inverted CPhase gate (:math:`CPhase(\lambda)^{\dagger} = CPhase(-\lambda)`)"""
         return CPhaseGate(-self.params[0], ctrl_state=self.ctrl_state)
 
-    def __array__(self, dtype=None):
+    def __array__(self, dtype=None, copy=None):
         """Return a numpy.array for the CPhase gate."""
+        if copy is False:
+            raise ValueError("cannot produce matrix without calculation")
         eith = exp(1j * float(self.params[0]))
         if self.ctrl_state:
             return numpy.array(

qiskit/circuit/library/standard_gates/r.py

@@ -93,8 +93,10 @@ def inverse(self, annotated: bool = False):
         """
         return RGate(-self.params[0], self.params[1])
 
-    def __array__(self, dtype=None):
+    def __array__(self, dtype=None, copy=None):
         """Return a numpy.array for the R gate."""
+        if copy is False:
+            raise ValueError("cannot produce matrix without calculation")
         theta, phi = float(self.params[0]), float(self.params[1])
         cos = math.cos(theta / 2)
         sin = math.sin(theta / 2)

qiskit/circuit/library/standard_gates/rx.py

@@ -120,8 +120,10 @@ def inverse(self, annotated: bool = False):
         """
         return RXGate(-self.params[0])
 
-    def __array__(self, dtype=None):
+    def __array__(self, dtype=None, copy=None):
         """Return a numpy.array for the RX gate."""
+        if copy is False:
+            raise ValueError("cannot produce matrix without calculation")
         cos = math.cos(self.params[0] / 2)
         sin = math.sin(self.params[0] / 2)
         return numpy.array([[cos, -1j * sin], [-1j * sin, cos]], dtype=dtype)
@@ -264,8 +266,10 @@ def inverse(self, annotated: bool = False):
         """
         return CRXGate(-self.params[0], ctrl_state=self.ctrl_state)
 
-    def __array__(self, dtype=None):
+    def __array__(self, dtype=None, copy=None):
         """Return a numpy.array for the CRX gate."""
+        if copy is False:
+            raise ValueError("cannot produce matrix without calculation")
         half_theta = float(self.params[0]) / 2
         cos = math.cos(half_theta)
         isin = 1j * math.sin(half_theta)

qiskit/circuit/library/standard_gates/rxx.py

@@ -122,8 +122,10 @@ def inverse(self, annotated: bool = False):
         """
         return RXXGate(-self.params[0])
 
-    def __array__(self, dtype=None):
+    def __array__(self, dtype=None, copy=None):
         """Return a Numpy.array for the RXX gate."""
+        if copy is False:
+            raise ValueError("cannot produce matrix without calculation")
         theta2 = float(self.params[0]) / 2
         cos = math.cos(theta2)
         isin = 1j * math.sin(theta2)

qiskit/circuit/library/standard_gates/ry.py

@@ -119,8 +119,10 @@ def inverse(self, annotated: bool = False):
         """
         return RYGate(-self.params[0])
 
-    def __array__(self, dtype=None):
+    def __array__(self, dtype=None, copy=None):
         """Return a numpy.array for the RY gate."""
+        if copy is False:
+            raise ValueError("cannot produce matrix without calculation")
         cos = math.cos(self.params[0] / 2)
         sin = math.sin(self.params[0] / 2)
         return numpy.array([[cos, -sin], [sin, cos]], dtype=dtype)
@@ -259,8 +261,10 @@ def inverse(self, annotated: bool = False):
         ."""
         return CRYGate(-self.params[0], ctrl_state=self.ctrl_state)
 
-    def __array__(self, dtype=None):
+    def __array__(self, dtype=None, copy=None):
         """Return a numpy.array for the CRY gate."""
+        if copy is False:
+            raise ValueError("cannot produce matrix without calculation")
         half_theta = float(self.params[0]) / 2
         cos = math.cos(half_theta)
         sin = math.sin(half_theta)

qiskit/circuit/library/standard_gates/ryy.py

@@ -122,8 +122,10 @@ def inverse(self, annotated: bool = False):
         """
         return RYYGate(-self.params[0])
 
-    def __array__(self, dtype=None):
+    def __array__(self, dtype=None, copy=None):
         """Return a numpy.array for the RYY gate."""
+        if copy is False:
+            raise ValueError("cannot produce matrix without calculation")
         theta = float(self.params[0])
         cos = math.cos(theta / 2)
         isin = 1j * math.sin(theta / 2)