MAINT: shims for array copy semantics (#4482)

* MAINT: shims for array copy semantics

* Fixes #4481

* this gets rid of segfaults/many test failurse against
latest NumPy `main` while still passing the full suite
against NumPy `1.26.4`; it is pretty hard to test since
needs special branch of SciPy (for the same reason, upstream
hasn't fully resolved copy semantics)

* I still see 2-3 test failures against `main` locally,
but likely still a large improvement; may want to let
this sit a little bit while the upstream storm settles,
but I suspect this is the gist of the shims we'd want

* MAINT: PR 4482 revisions

* abstract the array copy semantics (NumPy 1.x vs. 2.x) handling
to a utility function to reduce code duplication

* add a missing `dtype=None` to an `__array__` signature, for
more formal correctness with NumPy 2.x

* MAINT: PR 4482 revisions

* Update `vector_of_best_fit()` function in the helix analysis
code to avoid use of `np.linalg.linalg` namespace, which is
deprecated/private in NumPy `2.x`. Instead, just access `svd`
via the public `np.linalg` namespace.

* This fixes the last two remaining test failures for
MDAnalysis against NumPy `2.0.0rc1`.

* Switch the array `copy` shim to use the more appropriate
`NumpyVersion` rather than a raw string check for version `2`.

package/MDAnalysis/analysis/helix_analysis.py

@@ -121,7 +121,7 @@ def vector_of_best_fit(coordinates):
     """
     centered = coordinates - coordinates.mean(axis=0)
     Mt_M = np.matmul(centered.T, centered)
-    u, s, vh = np.linalg.linalg.svd(Mt_M)
+    u, s, vh = np.linalg.svd(Mt_M)
     vector = vh[0]
 
     # does vector face first local helix origin?

package/MDAnalysis/lib/formats/cython_util.pyx

@@ -26,6 +26,8 @@ cimport numpy as cnp
 from libc.stdlib cimport free
 from cpython cimport PyObject, Py_INCREF
 
+from MDAnalysis.lib.util import no_copy_shim
+
 cnp.import_array()
 
 
@@ -71,7 +73,7 @@ cdef class ArrayWrapper:
         self.data_type = data_type
         self.ndim = ndim
 
-    def __array__(self):
+    def __array__(self, dtype=None, copy=None):
         """ Here we use the __array__ method, that is called when numpy
             tries to get an array from the object."""
         ndarray = cnp.PyArray_SimpleNewFromData(self.ndim,
@@ -110,7 +112,7 @@ cdef cnp.ndarray ptr_to_ndarray(void* data_ptr, cnp.int64_t[:] dim, int data_typ
     array_wrapper = ArrayWrapper()
     array_wrapper.set_data(<void*> data_ptr, <int*> &dim[0], dim.size, data_type)
 
-    cdef cnp.ndarray ndarray = np.array(array_wrapper, copy=False)
+    cdef cnp.ndarray ndarray = np.array(array_wrapper, copy=no_copy_shim)
     # Assign our object to the 'base' of the ndarray object
     ndarray[:] = array_wrapper.__array__()
     # Increment the reference count, as the above assignement was done in

package/MDAnalysis/lib/transformations.py

@@ -170,6 +170,7 @@
 from numpy.linalg import norm
 
 from .mdamath import angle as vecangle
+from MDAnalysis.lib.util import no_copy_shim
 
 def identity_matrix():
     """Return 4x4 identity/unit matrix.
@@ -326,7 +327,7 @@ def rotation_matrix(angle, direction, point=None):
     M[:3, :3] = R
     if point is not None:
         # rotation not around origin
-        point = np.array(point[:3], dtype=np.float64, copy=False)
+        point = np.array(point[:3], dtype=np.float64, copy=no_copy_shim)
         M[:3, 3] = point - np.dot(R, point)
     return M
 
@@ -497,7 +498,7 @@ def projection_matrix(point, normal, direction=None,
 
     """
     M = np.identity(4)
-    point = np.array(point[:3], dtype=np.float64, copy=False)
+    point = np.array(point[:3], dtype=np.float64, copy=no_copy_shim)
     normal = unit_vector(normal[:3])
     if perspective is not None:
         # perspective projection
@@ -515,7 +516,7 @@ def projection_matrix(point, normal, direction=None,
         M[3, 3] = np.dot(perspective, normal)
     elif direction is not None:
         # parallel projection
-        direction = np.array(direction[:3], dtype=np.float64, copy=False)
+        direction = np.array(direction[:3], dtype=np.float64, copy=no_copy_shim)
         scale = np.dot(direction, normal)
         M[:3, :3] -= np.outer(direction, normal) / scale
         M[:3, 3] = direction * (np.dot(point, normal) / scale)
@@ -970,8 +971,8 @@ def superimposition_matrix(v0, v1, scaling=False, usesvd=True):
     True
 
     """
-    v0 = np.array(v0, dtype=np.float64, copy=False)[:3]
-    v1 = np.array(v1, dtype=np.float64, copy=False)[:3]
+    v0 = np.array(v0, dtype=np.float64, copy=no_copy_shim)[:3]
+    v1 = np.array(v1, dtype=np.float64, copy=no_copy_shim)[:3]
 
     if v0.shape != v1.shape or v0.shape[1] < 3:
         raise ValueError("vector sets are of wrong shape or type")
@@ -1314,7 +1315,7 @@ def quaternion_from_matrix(matrix, isprecise=False):
     True
 
     """
-    M = np.array(matrix, dtype=np.float64, copy=False)[:4, :4]
+    M = np.array(matrix, dtype=np.float64, copy=no_copy_shim)[:4, :4]
     if isprecise:
         q = np.empty((4, ), dtype=np.float64)
         t = np.trace(M)

package/MDAnalysis/lib/util.py

@@ -2552,3 +2552,11 @@ def wrapper(self, *args, **kwargs):
                         self._kwargs[key] = arg
         return func(self, *args, **kwargs)
     return wrapper
+
+
+def no_copy_shim():
+    if np.lib.NumpyVersion >= "2.0.0rc1":
+        copy = None
+    else:
+        copy = False
+    return copy