From 9da0fd4313c6931acc7e5327705e738b0ac320ef Mon Sep 17 00:00:00 2001 From: mrava87 Date: Fri, 26 Jul 2024 22:53:38 +0300 Subject: [PATCH 01/12] build: remove forcing numpy<2 --- environment-dev-arm.yml | 2 +- environment-dev.yml | 2 +- environment.yml | 2 +- pyproject.toml | 2 +- requirements-dev.txt | 2 +- 5 files changed, 5 insertions(+), 5 deletions(-) diff --git a/environment-dev-arm.yml b/environment-dev-arm.yml index a84f7199..413a9759 100755 --- a/environment-dev-arm.yml +++ b/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 diff --git a/environment-dev.yml b/environment-dev.yml index 5922a184..ef51f696 100755 --- a/environment-dev.yml +++ b/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 diff --git a/environment.yml b/environment.yml index 43df259a..e09650de 100755 --- a/environment.yml +++ b/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 diff --git a/pyproject.toml b/pyproject.toml index 05604374..6144f6e2 100644 --- a/pyproject.toml +++ b/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"] diff --git a/requirements-dev.txt b/requirements-dev.txt index ef8bbc55..703b377f 100644 --- a/requirements-dev.txt +++ b/requirements-dev.txt @@ -1,4 +1,4 @@ -numpy>=1.21.0,<2.0.0 +numpy>=1.21.0 scipy>=1.11.0 --extra-index-url https://download.pytorch.org/whl/cpu torch>=1.2.0 From 3f9b941884e14e0865b1fc595ef6634bd132c3d9 Mon Sep 17 00:00:00 2001 From: mrava87 Date: Fri, 26 Jul 2024 23:04:31 +0300 Subject: [PATCH 02/12] bug: fix linearoperator for numpy2 --- pylops/linearoperator.py | 27 +++++++++++++-------------- 1 file changed, 13 insertions(+), 14 deletions(-) diff --git a/pylops/linearoperator.py b/pylops/linearoperator.py index 0a719cf3..44e561cd 100644 --- a/pylops/linearoperator.py +++ b/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: From 28cc528abeed5a07d1906e692be4dc4b7e868fff Mon Sep 17 00:00:00 2001 From: mrava87 Date: Fri, 26 Jul 2024 23:17:55 +0300 Subject: [PATCH 03/12] test: add safeguards to dtcwt and spgl1 for numpy2 --- pytests/test_dtcwt.py | 80 +++++++++++++++++++++++-------------------- 1 file changed, 42 insertions(+), 38 deletions(-) diff --git a/pytests/test_dtcwt.py b/pytests/test_dtcwt.py index b0cf2b61..1bbba529 100644 --- a/pytests/test_dtcwt.py +++ b/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,66 +20,67 @@ 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: + t = sequential_array((par["ny"],)) - t = sequential_array((par["ny"],)) - - for level in range(1, 10): - Dtcwt = DTCWT(dims=t.shape, level=level, dtype=par["dtype"]) - x = Dtcwt @ t - y = Dtcwt.H @ x + for level in range(1, 10): + Dtcwt = DTCWT(dims=t.shape, level=level, dtype=par["dtype"]) + x = Dtcwt @ t + y = Dtcwt.H @ x - np.testing.assert_allclose(t, y) + np.testing.assert_allclose(t, y) @pytest.mark.parametrize("par", [(par1), (par2)]) def test_dtcwt1D_input2D(par): """Test for DTCWT with 2D input (forward-inverse pair)""" - - t = sequential_array( - ( - par["ny"], - par["ny"], + if int(np_version[0]) < 2: + t = sequential_array( + ( + par["ny"], + par["ny"], + ) ) - ) - for level in range(1, 10): - Dtcwt = DTCWT(dims=t.shape, level=level, dtype=par["dtype"]) - x = Dtcwt @ t - y = Dtcwt.H @ x + for level in range(1, 10): + Dtcwt = DTCWT(dims=t.shape, level=level, dtype=par["dtype"]) + x = Dtcwt @ t + y = Dtcwt.H @ x - np.testing.assert_allclose(t, y) + np.testing.assert_allclose(t, y) @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: + t = sequential_array((par["ny"], par["ny"], par["ny"])) - t = sequential_array((par["ny"], par["ny"], par["ny"])) + for level in range(1, 10): + Dtcwt = DTCWT(dims=t.shape, level=level, dtype=par["dtype"]) + x = Dtcwt @ t + y = Dtcwt.H @ x - for level in range(1, 10): - Dtcwt = DTCWT(dims=t.shape, level=level, dtype=par["dtype"]) - x = Dtcwt @ t - y = Dtcwt.H @ x - - np.testing.assert_allclose(t, y) + np.testing.assert_allclose(t, y) @pytest.mark.parametrize("par", [(par1), (par2)]) def test_dtcwt1D_birot(par): """Test for DTCWT birot (forward-inverse pair)""" - birots = ["antonini", "legall", "near_sym_a", "near_sym_b"] - - t = sequential_array( - ( - par["ny"], - par["ny"], + if int(np_version[0]) < 2: + birots = ["antonini", "legall", "near_sym_a", "near_sym_b"] + + t = sequential_array( + ( + par["ny"], + par["ny"], + ) ) - ) - for _b in birots: - print(f"birot {_b}") - Dtcwt = DTCWT(dims=t.shape, biort=_b, dtype=par["dtype"]) - x = Dtcwt @ t - y = Dtcwt.H @ x + for _b in birots: + print(f"birot {_b}") + Dtcwt = DTCWT(dims=t.shape, biort=_b, dtype=par["dtype"]) + x = Dtcwt @ t + y = Dtcwt.H @ x - np.testing.assert_allclose(t, y) + np.testing.assert_allclose(t, y) From 59de6504375487e18857a0bd41d5ff0464360af7 Mon Sep 17 00:00:00 2001 From: mrava87 Date: Fri, 26 Jul 2024 23:19:36 +0300 Subject: [PATCH 04/12] doc: force readthedocs to use numpy=1.21.0,<2.0.0 +scipy>=1.11.0 +--extra-index-url https://download.pytorch.org/whl/cpu +torch>=1.2.0 +numba +pyfftw +PyWavelets +spgl1 +scikit-fmm +sympy +devito +dtcwt +matplotlib +ipython +pytest +pytest-runner +setuptools_scm +docutils<0.18 +Sphinx +pydata-sphinx-theme +sphinx-gallery +numpydoc +nbsphinx +image +pre-commit +autopep8 +isort +black +flake8 +mypy From 85dfc1b48720c8e2cca17962a9cd52c715dcc5b4 Mon Sep 17 00:00:00 2001 From: mrava87 Date: Fri, 26 Jul 2024 23:30:40 +0300 Subject: [PATCH 05/12] doc: force also test_sparsity --- pytests/test_sparsity.py | 38 +++++++++++++++++++++----------------- 1 file changed, 21 insertions(+), 17 deletions(-) diff --git a/pytests/test_sparsity.py b/pytests/test_sparsity.py index b4ef5a30..c87bc690 100644 --- a/pytests/test_sparsity.py +++ b/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,24 +362,25 @@ def test_ISTA_FISTA_multiplerhs(par): ) def test_SPGL1(par): """Invert problem with SPGL1""" - np.random.seed(42) - Aop = MatrixMult(np.random.randn(par["ny"], par["nx"])) - - x = np.zeros(par["nx"]) - x[par["nx"] // 2] = 1 - x[3] = 1 - x[par["nx"] - 4] = -1 - - x0 = ( - np.random.normal(0, 10, par["nx"]) - + par["imag"] * np.random.normal(0, 10, par["nx"]) - if par["x0"] - else None - ) - y = Aop * x - xinv = spgl1(Aop, y, x0=x0, iter_lim=5000)[0] + if int(np_version[0]) < 2: + np.random.seed(42) + Aop = MatrixMult(np.random.randn(par["ny"], par["nx"])) + + x = np.zeros(par["nx"]) + x[par["nx"] // 2] = 1 + x[3] = 1 + x[par["nx"] - 4] = -1 + + x0 = ( + np.random.normal(0, 10, par["nx"]) + + par["imag"] * np.random.normal(0, 10, par["nx"]) + if par["x0"] + else None + ) + y = Aop * x + xinv = spgl1(Aop, y, x0=x0, iter_lim=5000)[0] - assert_array_almost_equal(x, xinv, decimal=1) + assert_array_almost_equal(x, xinv, decimal=1) @pytest.mark.parametrize("par", [(par1), (par2), (par1j), (par2j)]) From 599c5ab608e029ba4ab97aa91bed32ef87ca1cbe Mon Sep 17 00:00:00 2001 From: mrava87 Date: Sat, 27 Jul 2024 12:21:17 +0300 Subject: [PATCH 06/12] tests: disable tests for TorchOperator on darwin os --- pytests/test_torchoperator.py | 79 ++++++++++++++++++++--------------- 1 file changed, 46 insertions(+), 33 deletions(-) diff --git a/pytests/test_torchoperator.py b/pytests/test_torchoperator.py index 38246a20..a74bfaa2 100755 --- a/pytests/test_torchoperator.py +++ b/pytests/test_torchoperator.py @@ -1,3 +1,5 @@ +import platform + import numpy as np import pytest import torch @@ -17,53 +19,64 @@ def test_TorchOperator(par): must equal the adjoint of operator applied to the same vector, the two results are also checked to be the same. """ - Dop = MatrixMult(np.random.normal(0.0, 1.0, (par["ny"], par["nx"]))) - Top = TorchOperator(Dop, batch=False) + # temporarily, skip tests on mac as torch seems not to recognized + # numpy when v2 is installed + if platform.system is not "Darwin": + Dop = MatrixMult(np.random.normal(0.0, 1.0, (par["ny"], par["nx"]))) + Top = TorchOperator(Dop, batch=False) - x = np.random.normal(0.0, 1.0, par["nx"]) - xt = torch.from_numpy(x).view(-1) - xt.requires_grad = True - v = torch.randn(par["ny"]) + x = np.random.normal(0.0, 1.0, par["nx"]) + xt = torch.from_numpy(x).view(-1) + xt.requires_grad = True + v = torch.randn(par["ny"]) - # pylops operator - y = Dop * x - xadj = Dop.H * v + # pylops operator + y = Dop * x + xadj = Dop.H * v - # torch operator - yt = Top.apply(xt) - yt.backward(v, retain_graph=True) + # torch operator + yt = Top.apply(xt) + yt.backward(v, retain_graph=True) - assert_array_equal(y, yt.detach().cpu().numpy()) - assert_array_equal(xadj, xt.grad.cpu().numpy()) + assert_array_equal(y, yt.detach().cpu().numpy()) + assert_array_equal(xadj, xt.grad.cpu().numpy()) @pytest.mark.parametrize("par", [(par1)]) def test_TorchOperator_batch(par): """Apply forward for input with multiple samples (= batch) and flattened arrays""" - Dop = MatrixMult(np.random.normal(0.0, 1.0, (par["ny"], par["nx"]))) - Top = TorchOperator(Dop, batch=True) + # temporarily, skip tests on mac as torch seems not to recognized + # numpy when v2 is installed + if platform.system is not "Darwin": + Dop = MatrixMult(np.random.normal(0.0, 1.0, (par["ny"], par["nx"]))) + Top = TorchOperator(Dop, batch=True) - x = np.random.normal(0.0, 1.0, (4, par["nx"])) - xt = torch.from_numpy(x) - xt.requires_grad = True + x = np.random.normal(0.0, 1.0, (4, par["nx"])) + xt = torch.from_numpy(x) + xt.requires_grad = True - y = Dop.matmat(x.T).T - yt = Top.apply(xt) + y = Dop.matmat(x.T).T + yt = Top.apply(xt) - assert_array_equal(y, yt.detach().cpu().numpy()) + assert_array_equal(y, yt.detach().cpu().numpy()) @pytest.mark.parametrize("par", [(par1)]) def test_TorchOperator_batch_nd(par): """Apply forward for input with multiple samples (= batch) and nd-arrays""" - Dop = MatrixMult(np.random.normal(0.0, 1.0, (par["ny"], par["nx"])), otherdims=(2,)) - Top = TorchOperator(Dop, batch=True, flatten=False) - - x = np.random.normal(0.0, 1.0, (4, par["nx"], 2)) - xt = torch.from_numpy(x) - xt.requires_grad = True - - y = (Dop @ x.transpose(1, 2, 0)).transpose(2, 0, 1) - yt = Top.apply(xt) - - assert_array_equal(y, yt.detach().cpu().numpy()) + # temporarily, skip tests on mac as torch seems not to recognized + # numpy when v2 is installed + if platform.system is not "Darwin": + Dop = MatrixMult( + np.random.normal(0.0, 1.0, (par["ny"], par["nx"])), otherdims=(2,) + ) + Top = TorchOperator(Dop, batch=True, flatten=False) + + x = np.random.normal(0.0, 1.0, (4, par["nx"], 2)) + xt = torch.from_numpy(x) + xt.requires_grad = True + + y = (Dop @ x.transpose(1, 2, 0)).transpose(2, 0, 1) + yt = Top.apply(xt) + + assert_array_equal(y, yt.detach().cpu().numpy()) From 24848c49d9e5d3fc99ce141b1499cb653f8c2c43 Mon Sep 17 00:00:00 2001 From: mrava87 Date: Sat, 27 Jul 2024 12:29:46 +0300 Subject: [PATCH 07/12] minor: change numpy import to suppress depracation warning --- pylops/basicoperators/restriction.py | 12 ++++++++---- pytests/test_torchoperator.py | 6 +++--- 2 files changed, 11 insertions(+), 7 deletions(-) diff --git a/pylops/basicoperators/restriction.py b/pylops/basicoperators/restriction.py index c2e51a31..865e155a 100644 --- a/pylops/basicoperators/restriction.py +++ b/pylops/basicoperators/restriction.py @@ -1,12 +1,11 @@ __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 +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 +127,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) diff --git a/pytests/test_torchoperator.py b/pytests/test_torchoperator.py index a74bfaa2..bc022a9b 100755 --- a/pytests/test_torchoperator.py +++ b/pytests/test_torchoperator.py @@ -21,7 +21,7 @@ def test_TorchOperator(par): """ # temporarily, skip tests on mac as torch seems not to recognized # numpy when v2 is installed - if platform.system is not "Darwin": + if platform.system() != "Darwin": Dop = MatrixMult(np.random.normal(0.0, 1.0, (par["ny"], par["nx"]))) Top = TorchOperator(Dop, batch=False) @@ -47,7 +47,7 @@ def test_TorchOperator_batch(par): """Apply forward for input with multiple samples (= batch) and flattened arrays""" # temporarily, skip tests on mac as torch seems not to recognized # numpy when v2 is installed - if platform.system is not "Darwin": + if platform.system() != "Darwin": Dop = MatrixMult(np.random.normal(0.0, 1.0, (par["ny"], par["nx"]))) Top = TorchOperator(Dop, batch=True) @@ -66,7 +66,7 @@ def test_TorchOperator_batch_nd(par): """Apply forward for input with multiple samples (= batch) and nd-arrays""" # temporarily, skip tests on mac as torch seems not to recognized # numpy when v2 is installed - if platform.system is not "Darwin": + if platform.system() != "Darwin": Dop = MatrixMult( np.random.normal(0.0, 1.0, (par["ny"], par["nx"])), otherdims=(2,) ) From 45b5fbc6739304225de5ac8093b34223527120d4 Mon Sep 17 00:00:00 2001 From: mrava87 Date: Sat, 27 Jul 2024 12:40:32 +0300 Subject: [PATCH 08/12] minor: different import in Restriction based on numpy version Allow Restriction using normalize_axis_index differently based on numpy version to enable using numpy2 whilst still using numpy1 for rtd. --- pylops/basicoperators/restriction.py | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/pylops/basicoperators/restriction.py b/pylops/basicoperators/restriction.py index 865e155a..fc81a252 100644 --- a/pylops/basicoperators/restriction.py +++ b/pylops/basicoperators/restriction.py @@ -5,7 +5,14 @@ import numpy as np import numpy.ma as np_ma -from numpy.lib.array_utils 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 From b9618f547baf133374fd8c5e5cccf4c7ff327a66 Mon Sep 17 00:00:00 2001 From: Carlos da Costa Date: Sun, 4 Aug 2024 18:38:49 -0700 Subject: [PATCH 09/12] chore: simplify testing logic --- pytests/test_dtcwt.py | 86 +++++++++++++++++++++++-------------------- 1 file changed, 47 insertions(+), 39 deletions(-) diff --git a/pytests/test_dtcwt.py b/pytests/test_dtcwt.py index 1bbba529..979a7f76 100644 --- a/pytests/test_dtcwt.py +++ b/pytests/test_dtcwt.py @@ -20,67 +20,75 @@ 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: - t = sequential_array((par["ny"],)) + if int(np_version[0]) >= 2: + return - for level in range(1, 10): - Dtcwt = DTCWT(dims=t.shape, level=level, dtype=par["dtype"]) - x = Dtcwt @ t - y = Dtcwt.H @ x + t = sequential_array((par["ny"],)) - np.testing.assert_allclose(t, y) + for level in range(1, 10): + Dtcwt = DTCWT(dims=t.shape, level=level, dtype=par["dtype"]) + x = Dtcwt @ t + y = Dtcwt.H @ x + + np.testing.assert_allclose(t, y) @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: - t = sequential_array( - ( - par["ny"], - par["ny"], - ) + if int(np_version[0]) >= 2: + return + + t = sequential_array( + ( + par["ny"], + par["ny"], ) + ) - for level in range(1, 10): - Dtcwt = DTCWT(dims=t.shape, level=level, dtype=par["dtype"]) - x = Dtcwt @ t - y = Dtcwt.H @ x + for level in range(1, 10): + Dtcwt = DTCWT(dims=t.shape, level=level, dtype=par["dtype"]) + x = Dtcwt @ t + y = Dtcwt.H @ x - np.testing.assert_allclose(t, y) + np.testing.assert_allclose(t, y) @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: - t = sequential_array((par["ny"], par["ny"], par["ny"])) + if int(np_version[0]) >= 2: + return + + t = sequential_array((par["ny"], par["ny"], par["ny"])) - for level in range(1, 10): - Dtcwt = DTCWT(dims=t.shape, level=level, dtype=par["dtype"]) - x = Dtcwt @ t - y = Dtcwt.H @ x + for level in range(1, 10): + Dtcwt = DTCWT(dims=t.shape, level=level, dtype=par["dtype"]) + x = Dtcwt @ t + y = Dtcwt.H @ x - np.testing.assert_allclose(t, y) + np.testing.assert_allclose(t, y) @pytest.mark.parametrize("par", [(par1), (par2)]) def test_dtcwt1D_birot(par): """Test for DTCWT birot (forward-inverse pair)""" - if int(np_version[0]) < 2: - birots = ["antonini", "legall", "near_sym_a", "near_sym_b"] - - t = sequential_array( - ( - par["ny"], - par["ny"], - ) + if int(np_version[0]) >= 2: + return + + birots = ["antonini", "legall", "near_sym_a", "near_sym_b"] + + t = sequential_array( + ( + par["ny"], + par["ny"], ) + ) - for _b in birots: - print(f"birot {_b}") - Dtcwt = DTCWT(dims=t.shape, biort=_b, dtype=par["dtype"]) - x = Dtcwt @ t - y = Dtcwt.H @ x + for _b in birots: + print(f"birot {_b}") + Dtcwt = DTCWT(dims=t.shape, biort=_b, dtype=par["dtype"]) + x = Dtcwt @ t + y = Dtcwt.H @ x - np.testing.assert_allclose(t, y) + np.testing.assert_allclose(t, y) From 734776d16d4947098b5d74792b9ae1adbc160ae6 Mon Sep 17 00:00:00 2001 From: Carlos da Costa Date: Sun, 4 Aug 2024 18:41:08 -0700 Subject: [PATCH 10/12] chore: simplify testing logic --- pytests/test_sparsity.py | 40 +++++++++++++++++++++------------------- 1 file changed, 21 insertions(+), 19 deletions(-) diff --git a/pytests/test_sparsity.py b/pytests/test_sparsity.py index c87bc690..00c7d944 100644 --- a/pytests/test_sparsity.py +++ b/pytests/test_sparsity.py @@ -362,25 +362,27 @@ def test_ISTA_FISTA_multiplerhs(par): ) def test_SPGL1(par): """Invert problem with SPGL1""" - if int(np_version[0]) < 2: - np.random.seed(42) - Aop = MatrixMult(np.random.randn(par["ny"], par["nx"])) - - x = np.zeros(par["nx"]) - x[par["nx"] // 2] = 1 - x[3] = 1 - x[par["nx"] - 4] = -1 - - x0 = ( - np.random.normal(0, 10, par["nx"]) - + par["imag"] * np.random.normal(0, 10, par["nx"]) - if par["x0"] - else None - ) - y = Aop * x - xinv = spgl1(Aop, y, x0=x0, iter_lim=5000)[0] + if int(np_version[0]) >= 2: + return - assert_array_almost_equal(x, xinv, decimal=1) + np.random.seed(42) + Aop = MatrixMult(np.random.randn(par["ny"], par["nx"])) + + x = np.zeros(par["nx"]) + x[par["nx"] // 2] = 1 + x[3] = 1 + x[par["nx"] - 4] = -1 + + x0 = ( + np.random.normal(0, 10, par["nx"]) + + par["imag"] * np.random.normal(0, 10, par["nx"]) + if par["x0"] + else None + ) + y = Aop * x + xinv = spgl1(Aop, y, x0=x0, iter_lim=5000)[0] + + assert_array_almost_equal(x, xinv, decimal=1) @pytest.mark.parametrize("par", [(par1), (par2), (par1j), (par2j)]) @@ -416,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 From f8e332301ca59d09988a8ef4682c22575f1d95cd Mon Sep 17 00:00:00 2001 From: Carlos da Costa Date: Sun, 4 Aug 2024 18:42:34 -0700 Subject: [PATCH 11/12] chore: simplify testing logic --- pytests/test_torchoperator.py | 74 ++++++++++++++++++----------------- 1 file changed, 39 insertions(+), 35 deletions(-) diff --git a/pytests/test_torchoperator.py b/pytests/test_torchoperator.py index bc022a9b..43f33e3f 100755 --- a/pytests/test_torchoperator.py +++ b/pytests/test_torchoperator.py @@ -21,25 +21,27 @@ def test_TorchOperator(par): """ # temporarily, skip tests on mac as torch seems not to recognized # numpy when v2 is installed - if platform.system() != "Darwin": - Dop = MatrixMult(np.random.normal(0.0, 1.0, (par["ny"], par["nx"]))) - Top = TorchOperator(Dop, batch=False) + if platform.system() == "Darwin": + return - x = np.random.normal(0.0, 1.0, par["nx"]) - xt = torch.from_numpy(x).view(-1) - xt.requires_grad = True - v = torch.randn(par["ny"]) + Dop = MatrixMult(np.random.normal(0.0, 1.0, (par["ny"], par["nx"]))) + Top = TorchOperator(Dop, batch=False) - # pylops operator - y = Dop * x - xadj = Dop.H * v + x = np.random.normal(0.0, 1.0, par["nx"]) + xt = torch.from_numpy(x).view(-1) + xt.requires_grad = True + v = torch.randn(par["ny"]) - # torch operator - yt = Top.apply(xt) - yt.backward(v, retain_graph=True) + # pylops operator + y = Dop * x + xadj = Dop.H * v - assert_array_equal(y, yt.detach().cpu().numpy()) - assert_array_equal(xadj, xt.grad.cpu().numpy()) + # torch operator + yt = Top.apply(xt) + yt.backward(v, retain_graph=True) + + assert_array_equal(y, yt.detach().cpu().numpy()) + assert_array_equal(xadj, xt.grad.cpu().numpy()) @pytest.mark.parametrize("par", [(par1)]) @@ -47,18 +49,20 @@ def test_TorchOperator_batch(par): """Apply forward for input with multiple samples (= batch) and flattened arrays""" # temporarily, skip tests on mac as torch seems not to recognized # numpy when v2 is installed - if platform.system() != "Darwin": - Dop = MatrixMult(np.random.normal(0.0, 1.0, (par["ny"], par["nx"]))) - Top = TorchOperator(Dop, batch=True) + if platform.system() == "Darwin": + return + + Dop = MatrixMult(np.random.normal(0.0, 1.0, (par["ny"], par["nx"]))) + Top = TorchOperator(Dop, batch=True) - x = np.random.normal(0.0, 1.0, (4, par["nx"])) - xt = torch.from_numpy(x) - xt.requires_grad = True + x = np.random.normal(0.0, 1.0, (4, par["nx"])) + xt = torch.from_numpy(x) + xt.requires_grad = True - y = Dop.matmat(x.T).T - yt = Top.apply(xt) + y = Dop.matmat(x.T).T + yt = Top.apply(xt) - assert_array_equal(y, yt.detach().cpu().numpy()) + assert_array_equal(y, yt.detach().cpu().numpy()) @pytest.mark.parametrize("par", [(par1)]) @@ -66,17 +70,17 @@ def test_TorchOperator_batch_nd(par): """Apply forward for input with multiple samples (= batch) and nd-arrays""" # temporarily, skip tests on mac as torch seems not to recognized # numpy when v2 is installed - if platform.system() != "Darwin": - Dop = MatrixMult( - np.random.normal(0.0, 1.0, (par["ny"], par["nx"])), otherdims=(2,) - ) - Top = TorchOperator(Dop, batch=True, flatten=False) + if platform.system() == "Darwin": + return + + Dop = MatrixMult(np.random.normal(0.0, 1.0, (par["ny"], par["nx"])), otherdims=(2,)) + Top = TorchOperator(Dop, batch=True, flatten=False) - x = np.random.normal(0.0, 1.0, (4, par["nx"], 2)) - xt = torch.from_numpy(x) - xt.requires_grad = True + x = np.random.normal(0.0, 1.0, (4, par["nx"], 2)) + xt = torch.from_numpy(x) + xt.requires_grad = True - y = (Dop @ x.transpose(1, 2, 0)).transpose(2, 0, 1) - yt = Top.apply(xt) + y = (Dop @ x.transpose(1, 2, 0)).transpose(2, 0, 1) + yt = Top.apply(xt) - assert_array_equal(y, yt.detach().cpu().numpy()) + assert_array_equal(y, yt.detach().cpu().numpy()) From 65102078793e0a7cc86700d1d49ee4108bd5b0d7 Mon Sep 17 00:00:00 2001 From: Carlos da Costa Date: Sun, 4 Aug 2024 19:25:14 -0700 Subject: [PATCH 12/12] fix: add warning that NumPy 2 is not supported with dtcwt and spgl1 --- docs/source/installation.rst | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/docs/source/installation.rst b/docs/source/installation.rst index 8de4c0be..3b6cd49d 100755 --- a/docs/source/installation.rst +++ b/docs/source/installation.rst @@ -321,6 +321,11 @@ In alphabetic order: dtcwt ----- + +.. warning:: + + ``dtcwt`` is not yet supported with Numpy 2. + `dtcwt `_ 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 `_ 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 `_ is used to solve sparsity-promoting basis pursuit, basis pursuit denoise, and Lasso problems in :py:func:`pylops.optimization.sparsity.SPGL1` solver.