Merge pull request #596 from mrava87/patch_numpyv2

Feature: migrate to numpy v2.0.0

Author: mrava87

.readthedocs.yaml

@@ -18,6 +18,6 @@ sphinx:
 # Declare the Python requirements required to build your docs
 python:
   install:
-    - requirements: requirements-dev.txt
+    - requirements: requirements-doc.txt
     - method: pip
       path: .

docs/source/installation.rst

@@ -321,6 +321,11 @@ In alphabetic order:
 
 dtcwt
 -----
+
+.. warning::
+
+   ``dtcwt`` is not yet supported with Numpy 2.
+
 `dtcwt <https://dtcwt.readthedocs.io/en/0.12.0/>`_ is a library used to implement the DT-CWT operators.
 
 Install it via ``pip`` with:
@@ -330,6 +335,7 @@ Install it via ``pip`` with:
    >> pip install dtcwt
 
 
+
 Devito
 ------
 `Devito <https://github.com/devitocodes/devito>`_ is a library used to solve PDEs via
@@ -468,6 +474,11 @@ or with ``pip`` via
 
 SPGL1
 -----
+
+.. warning::
+   
+   ``SPGL1`` is not yet supported with Numpy 2.
+
 `SPGL1 <https://spgl1.readthedocs.io/en/latest/>`_ is used to solve sparsity-promoting
 basis pursuit, basis pursuit denoise, and Lasso problems
 in :py:func:`pylops.optimization.sparsity.SPGL1` solver.

environment-dev-arm.yml

@@ -7,7 +7,7 @@ channels:
 dependencies:
   - python>=3.6.4
   - pip
-  - numpy>=1.21.0,<2.0.0
+  - numpy>=1.21.0
   - scipy>=1.11.0
   - pytorch>=1.2.0
   - cpuonly

environment-dev.yml

@@ -7,7 +7,7 @@ channels:
 dependencies:
   - python>=3.6.4
   - pip
-  - numpy>=1.21.0,<2.0.0
+  - numpy>=1.21.0
   - scipy>=1.11.0
   - pytorch>=1.2.0
   - cpuonly

environment.yml

@@ -3,5 +3,5 @@ channels:
   - defaults
 dependencies:
   - python>=3.6.4
-  - numpy>=1.21.0,<2.0.0
+  - numpy>=1.21.0
   - scipy>=1.14.0

pylops/basicoperators/restriction.py

@@ -1,12 +1,18 @@
 __all__ = ["Restriction"]
 
 import logging
-
 from typing import Sequence, Union
 
 import numpy as np
 import numpy.ma as np_ma
-from numpy.core.multiarray import normalize_axis_index
+
+# need to check numpy version since normalize_axis_index will be
+# soon moved from numpy.core.multiarray to from numpy.lib.array_utils
+np_version = np.__version__.split(".")
+if int(np_version[0]) < 2:
+    from numpy.core.multiarray import normalize_axis_index
+else:
+    from numpy.lib.array_utils import normalize_axis_index
 
 from pylops import LinearOperator
 from pylops.utils._internal import _value_or_sized_to_tuple
@@ -128,8 +134,13 @@ def __init__(
                 )
                 forceflat = None
 
-        super().__init__(dtype=np.dtype(dtype), dims=dims, dimsd=dimsd,
-                         forceflat=forceflat, name=name)
+        super().__init__(
+            dtype=np.dtype(dtype),
+            dims=dims,
+            dimsd=dimsd,
+            forceflat=forceflat,
+            name=name,
+        )
 
         iavareshape = np.ones(len(self.dims), dtype=int)
         iavareshape[axis] = len(iava)

pylops/linearoperator.py

@@ -1242,23 +1242,14 @@ def _get_dtype(
 ) -> DTypeLike:
     if dtypes is None:
         dtypes = []
-    opdtypes = []
     for obj in operators:
         if obj is not None and hasattr(obj, "dtype"):
-            opdtypes.append(obj.dtype)
-    return np.find_common_type(opdtypes, dtypes)
+            dtypes.append(obj.dtype)
+    return np.result_type(*dtypes)
 
 
 class _ScaledLinearOperator(LinearOperator):
-    """
-    Sum Linear Operator
-
-    Modified version of scipy _ScaledLinearOperator which uses a modified
-    _get_dtype where the scalar and operator types are passed separately to
-    np.find_common_type. Passing them together does lead to problems when using
-    np.float32 operators which are cast to np.float64
-
-    """
+    """Scaled Linear Operator"""
 
     def __init__(
         self,
@@ -1269,7 +1260,15 @@ def __init__(
             raise ValueError("LinearOperator expected as A")
         if not np.isscalar(alpha):
             raise ValueError("scalar expected as alpha")
-        dtype = _get_dtype([A], [type(alpha)])
+        if isinstance(alpha, complex) and not np.iscomplexobj(
+            np.ones(1, dtype=A.dtype)
+        ):
+            # if the scalar is of complex type but not the operator, find out type
+            dtype = _get_dtype([A], [type(alpha)])
+        else:
+            # if both the scalar and operator are of real or complex type, use type
+            # of the operator
+            dtype = A.dtype
         super(_ScaledLinearOperator, self).__init__(dtype=dtype, shape=A.shape)
         self.args = (A, alpha)
 
@@ -1465,7 +1464,7 @@ def __init__(self, A: LinearOperator, p: int) -> None:
         if not isintlike(p) or p < 0:
             raise ValueError("non-negative integer expected as p")
 
-        super(_PowerLinearOperator, self).__init__(dtype=_get_dtype([A]), shape=A.shape)
+        super(_PowerLinearOperator, self).__init__(dtype=A.dtype, shape=A.shape)
         self.args = (A, p)
 
     def _power(self, fun: Callable, x: NDArray) -> NDArray:

pyproject.toml

@@ -30,7 +30,7 @@ classifiers = [
     "Topic :: Scientific/Engineering :: Mathematics",
 ]
 dependencies = [
-    "numpy >= 1.21.0 , < 2.0.0",
+    "numpy >= 1.21.0",
     "scipy >= 1.11.0",
 ]
 dynamic = ["version"]

pytests/test_dtcwt.py

@@ -3,6 +3,9 @@
 
 from pylops.signalprocessing import DTCWT
 
+# currently test only if numpy<2.0.0 is installed...
+np_version = np.__version__.split(".")
+
 par1 = {"ny": 10, "nx": 10, "dtype": "float64"}
 par2 = {"ny": 50, "nx": 50, "dtype": "float64"}
 
@@ -17,6 +20,8 @@ def sequential_array(shape):
 @pytest.mark.parametrize("par", [(par1), (par2)])
 def test_dtcwt1D_input1D(par):
     """Test for DTCWT with 1D input"""
+    if int(np_version[0]) >= 2:
+        return
 
     t = sequential_array((par["ny"],))
 
@@ -31,6 +36,8 @@ def test_dtcwt1D_input1D(par):
 @pytest.mark.parametrize("par", [(par1), (par2)])
 def test_dtcwt1D_input2D(par):
     """Test for DTCWT with 2D input (forward-inverse pair)"""
+    if int(np_version[0]) >= 2:
+        return
 
     t = sequential_array(
         (
@@ -50,6 +57,8 @@ def test_dtcwt1D_input2D(par):
 @pytest.mark.parametrize("par", [(par1), (par2)])
 def test_dtcwt1D_input3D(par):
     """Test for DTCWT with 3D input (forward-inverse pair)"""
+    if int(np_version[0]) >= 2:
+        return
 
     t = sequential_array((par["ny"], par["ny"], par["ny"]))
 
@@ -64,6 +73,9 @@ def test_dtcwt1D_input3D(par):
 @pytest.mark.parametrize("par", [(par1), (par2)])
 def test_dtcwt1D_birot(par):
     """Test for DTCWT birot (forward-inverse pair)"""
+    if int(np_version[0]) >= 2:
+        return
+
     birots = ["antonini", "legall", "near_sym_a", "near_sym_b"]
 
     t = sequential_array(

pytests/test_sparsity.py

@@ -5,6 +5,9 @@
 from pylops.basicoperators import FirstDerivative, Identity, MatrixMult
 from pylops.optimization.sparsity import fista, irls, ista, omp, spgl1, splitbregman
 
+# currently test spgl1 only if numpy<2.0.0 is installed...
+np_version = np.__version__.split(".")
+
 par1 = {
     "ny": 11,
     "nx": 11,
@@ -359,6 +362,9 @@ def test_ISTA_FISTA_multiplerhs(par):
 )
 def test_SPGL1(par):
     """Invert problem with SPGL1"""
+    if int(np_version[0]) >= 2:
+        return
+
     np.random.seed(42)
     Aop = MatrixMult(np.random.randn(par["ny"], par["nx"]))
 
@@ -412,6 +418,6 @@ def test_SplitBregman(par):
         x0=x0 if par["x0"] else None,
         restart=False,
         show=False,
-        **dict(iter_lim=5, damp=1e-3)
+        **dict(iter_lim=5, damp=1e-3),
     )
     assert (np.linalg.norm(x - xinv) / np.linalg.norm(x)) < 1e-1