BUG: avoid soon-to-be-deprecated direct mutations of ndarray.shape attributes

Author: neutrinoceros

CONTRIBUTING.rst

@@ -720,9 +720,7 @@ Source code style guide
 * Do not use nested classes unless you have a very good reason to, such as
   requiring a namespace or class-definition modification.  Classes should live
   at the top level.  ``__metaclass__`` is exempt from this.
-* Avoid copying memory when possible. For example, don't do
-  ``a = a.reshape(3, 4)`` when ``a.shape = (3, 4)`` will do, and ``a = a * 3``
-  should be ``np.multiply(a, 3, a)``.
+* Avoid copying memory when possible.
 * In general, avoid all double-underscore method names: ``__something`` is
   usually unnecessary.
 * When writing a subclass, use the super built-in to access the super class,

yt/_maintenance/numpy2_compat.py

@@ -1,9 +1,37 @@
 # avoid deprecation warnings in numpy >= 2.0
+from dataclasses import dataclass
+from importlib.metadata import version
+from typing import Literal, TypeAlias
 
 import numpy as np
+from packaging.version import Version
 
 if hasattr(np, "trapezoid"):
     # np.trapz is deprecated in numpy 2.0 in favor of np.trapezoid
     trapezoid = np.trapezoid
 else:
     trapezoid = np.trapz  # type: ignore # noqa: NPY201
+
+
+NUMPY_VERSION = Version(version("numpy"))
+
+CopyFunc: TypeAlias = Literal["reshape"]
+MIN_NUMPY_VERSION: dict[CopyFunc, Version] = {
+    "reshape": Version("2.1.0"),
+}
+
+
+@dataclass(frozen=True, slots=True)
+class CopyKwarg:
+    value: bool | None
+
+    def get(self, f: CopyFunc, /) -> dict[Literal["copy"], bool | None]:
+        if NUMPY_VERSION >= MIN_NUMPY_VERSION[f]:
+            return {"copy": self.value}
+        else:
+            return {}
+
+
+COPY_NONE = CopyKwarg(None)
+COPY_TRUE = CopyKwarg(True)
+COPY_FALSE = CopyKwarg(False)

yt/fields/field_detector.py

@@ -276,18 +276,19 @@ def fcoords(self):
             np.mgrid[0 : 1 : self.nd * 1j, 0 : 1 : self.nd * 1j, 0 : 1 : self.nd * 1j]
         )
         if self.flat:
-            fc.shape = (self.nd * self.nd * self.nd, 3)
+            fc = fc.reshape(self.nd * self.nd * self.nd, 3)
         else:
             fc = fc.transpose()
         return self.ds.arr(fc, units="code_length")
 
     @property
     def fcoords_vertex(self):
-        rng = np.random.default_rng()
-        fc = rng.random((self.nd, self.nd, self.nd, 8, 3))
         if self.flat:
-            fc.shape = (self.nd * self.nd * self.nd, 8, 3)
-        return self.ds.arr(fc, units="code_length")
+            shape = (self.nd * self.nd * self.nd, 8, 3)
+        else:
+            shape = (self.nd, self.nd, self.nd, 8, 3)
+        rng = np.random.default_rng()
+        return self.ds.arr(rng.random(shape), units="code_length")
 
     @property
     def icoords(self):
@@ -297,21 +298,23 @@ def icoords(self):
             0 : self.nd - 1 : self.nd * 1j,
         ]
         if self.flat:
-            ic.shape = (self.nd * self.nd * self.nd, 3)
+            return ic.reshape(self.nd * self.nd * self.nd, 3)
         else:
-            ic = ic.transpose()
-        return ic
+            return ic.transpose()
 
     @property
     def ires(self):
-        ir = np.ones(self.nd**3, dtype="int64")
-        if not self.flat:
-            ir.shape = (self.nd, self.nd, self.nd)
-        return ir
+        if self.flat:
+            shape = (self.nd**3,)
+        else:
+            shape = (self.nd, self.nd, self.nd)
+        return np.ones(shape, dtype="int64")
 
     @property
     def fwidth(self):
-        fw = np.ones((self.nd**3, 3), dtype="float64") / self.nd
-        if not self.flat:
-            fw.shape = (self.nd, self.nd, self.nd, 3)
+        if self.flat:
+            shape = (self.nd**3, 3)
+        else:
+            shape = (self.nd, self.nd, self.nd, 3)
+        fw = np.full(shape, 1 / self.nd, dtype="float64")
         return self.ds.arr(fw, units="code_length")

yt/fields/geometric_fields.py

@@ -1,5 +1,6 @@
 import numpy as np
 
+from yt._maintenance.numpy2_compat import COPY_FALSE
 from yt.utilities.lib.geometry_utils import compute_morton
 from yt.utilities.math_utils import (
     get_cyl_r,
@@ -128,14 +129,17 @@ def _morton_index(data):
         LE = data.ds.domain_left_edge - eps * uq
         RE = data.ds.domain_right_edge + eps * uq
         # .ravel() only copies if it needs to
-        morton = compute_morton(
-            data["index", "x"].ravel(),
-            data["index", "y"].ravel(),
-            data["index", "z"].ravel(),
-            LE,
-            RE,
+        morton = np.reshape(
+            compute_morton(
+                data["index", "x"].ravel(),
+                data["index", "y"].ravel(),
+                data["index", "z"].ravel(),
+                LE,
+                RE,
+            ),
+            data["index", "x"].shape,
+            **COPY_FALSE.get("reshape"),
         )
-        morton.shape = data["index", "x"].shape
         return morton.view("f8")
 
     registry.add_field(

yt/frontends/amrvac/datfile_utils.py

@@ -157,6 +157,7 @@ def get_single_block_field_data(istream, byte_offset, block_shape, field_idx):
     byte_size_field = struct.calcsize(fmt)
 
     istream.seek(byte_offset + byte_size_field * field_idx)
-    data = np.fromfile(istream, "=f8", count=np.prod(field_shape))
-    data.shape = field_shape[::-1]
+    data = np.fromfile(
+        istream, dtype=("=f8", field_shape[::-1]), count=np.prod(field_shape)
+    )
     return data.T

yt/frontends/amrvac/io.py

@@ -108,8 +108,7 @@ def _read_data(self, fid, grid, field):
         byte_size_field = count * 8  # size of a double
 
         fid.seek(offset + byte_size_field * field_idx)
-        data = np.fromfile(fid, "=f8", count=count)
-        data.shape = field_shape[::-1]
+        data = np.fromfile(fid, dtype=("=f8", field_shape[::-1]), count=count)
         data = data.T
         # Always convert data to 3D, as grid.ActiveDimensions is always 3D
         while len(data.shape) < 3:

yt/frontends/chimera/data_structures.py

@@ -10,6 +10,7 @@
 
 import numpy as np
 
+from yt._maintenance.numpy2_compat import COPY_FALSE
 from yt.data_objects.index_subobjects.unstructured_mesh import SemiStructuredMesh
 from yt.data_objects.static_output import Dataset
 from yt.geometry.api import Geometry
@@ -103,7 +104,9 @@ def _initialize_mesh(self):
                     mylog.warning(
                         "Yin-Yang File Detected; This data is not currently supported."
                     )
-                coords.shape = (nxd * nyd * nzd, 3)
+                coords = np.reshape(
+                    coords, (nxd * nyd * nzd, 3), **COPY_FALSE.get("reshape")
+                )
                 # Connectivity is an array of rows, each of which corresponds to a grid cell.
                 # The 8 elements of each row are integers representing the cell vertices.
                 # These integers reference the numerical index of the element of the

yt/frontends/gadget/io.py

@@ -462,10 +462,9 @@ def _yield_coordinates(self, data_file, needed_ptype=None):
                 if needed_ptype is not None and ptype != needed_ptype:
                     continue
                 # The first total_particles * 3 values are positions
-                pp = np.fromfile(f, dtype=dt, count=count * 3).astype(
+                pp = np.fromfile(f, dtype=(dt, (count, 3)), count=count * 3).astype(
                     dt_native, copy=False
                 )
-                pp.shape = (count, 3)
                 yield ptype, pp
 
     def _get_smoothing_length(self, data_file, position_dtype, position_shape):

yt/frontends/http_stream/io.py

@@ -1,5 +1,6 @@
 import numpy as np
 
+from yt._maintenance.numpy2_compat import COPY_FALSE
 from yt.funcs import mylog
 from yt.utilities.io_handler import BaseParticleIOHandler
 from yt.utilities.on_demand_imports import _requests as requests
@@ -36,26 +37,35 @@ def _read_particle_coords(self, chunks, ptf):
         for data_file in self._sorted_chunk_iterator(chunks):
             for ptype in ptf:
                 s = self._open_stream(data_file, (ptype, "Coordinates"))
-                c = np.frombuffer(s, dtype="float64")
-                c.shape = (c.shape[0] / 3.0, 3)
+                c_raw = np.frombuffer(s, dtype="float64")
+                c = np.reshape(
+                    c_raw, (c_raw.shape[0] // 3, 3), **COPY_FALSE.get("reshape")
+                )
                 yield ptype, (c[:, 0], c[:, 1], c[:, 2]), 0.0
 
     def _read_particle_fields(self, chunks, ptf, selector):
         # Now we have all the sizes, and we can allocate
         for data_file in self._sorted_chunk_iterator(chunks):
             for ptype, field_list in sorted(ptf.items()):
                 s = self._open_stream(data_file, (ptype, "Coordinates"))
-                c = np.frombuffer(s, dtype="float64")
-                c.shape = (c.shape[0] / 3.0, 3)
+                c_raw = np.frombuffer(s, dtype="float64")
+                c = np.reshape(
+                    c_raw, c_raw.shape[0] // 3, 3, **COPY_FALSE.get("reshape")
+                )
                 mask = selector.select_points(c[:, 0], c[:, 1], c[:, 2], 0.0)
                 del c
+                del c_raw
                 if mask is None:
                     continue
                 for field in field_list:
                     s = self._open_stream(data_file, (ptype, field))
-                    c = np.frombuffer(s, dtype="float64")
+                    c_raw = np.frombuffer(s, dtype="float64")
                     if field in self._vector_fields:
-                        c.shape = (c.shape[0] / 3.0, 3)
+                        c = np.reshape(
+                            c_raw,
+                            (c_raw.shape[0] // 3, 3),
+                            **COPY_FALSE.get("reshape"),
+                        )
                     data = c[mask, ...]
                     yield (ptype, field), data
 

yt/frontends/open_pmd/io.py

@@ -2,6 +2,7 @@
 
 import numpy as np
 
+from yt._maintenance.numpy2_compat import COPY_FALSE
 from yt.frontends.open_pmd.misc import get_component, is_const_component
 from yt.geometry.selection_routines import GridSelector
 from yt.utilities.io_handler import BaseIOHandler
@@ -189,13 +190,15 @@ def _read_fluid_selection(self, chunks, selector, fields, size):
                         continue
                     component = fname.replace("_", "/").replace("-", "_")
                     if component.split("/")[0] not in grid.ftypes:
-                        data = np.full(grid.ActiveDimensions, 0, dtype=np.float64)
+                        data = np.full_like(mask, 0)
                     else:
-                        data = get_component(ds, component, grid.findex, grid.foffset)
+                        data = np.reshape(  # Workaround - casts a 2D (x,y) array to 3D (x,y,1)
+                            get_component(ds, component, grid.findex, grid.foffset),
+                            mask.shape,
+                            **COPY_FALSE.get("reshape"),
+                        )
+
                     # The following is a modified AMRGridPatch.select(...)
-                    data.shape = (
-                        mask.shape
-                    )  # Workaround - casts a 2D (x,y) array to 3D (x,y,1)
                     count = grid.count(selector)
                     rv[field][ind[field] : ind[field] + count] = data[mask]
                     ind[field] += count