diff --git a/pylops/LinearOperator.py b/pylops/LinearOperator.py index ea07d8a3..04b859e0 100644 --- a/pylops/LinearOperator.py +++ b/pylops/LinearOperator.py @@ -11,8 +11,14 @@ from scipy.sparse.linalg import eigsh as sp_eigsh from scipy.sparse.linalg import lobpcg as sp_lobpcg from scipy.sparse.linalg import lsqr, spsolve -from scipy.sparse.linalg.interface import _ProductLinearOperator +# need to check scipy version since the interface submodule changed into +# _interface from scipy>=1.8.0 +sp_version = sp.__version__.split(".") +if int(sp_version[0]) <= 1 and int(sp_version[1]) < 8: + from scipy.sparse.linalg.interface import _ProductLinearOperator +else: + from scipy.sparse.linalg._interface import _ProductLinearOperator from pylops.optimization.solver import cgls from pylops.utils.backend import get_array_module, get_module, get_sparse_eye from pylops.utils.estimators import trace_hutchinson, trace_hutchpp, trace_nahutchpp @@ -298,7 +304,7 @@ def __truediv__(self, y, niter=100): else: if isinstance(y, np.ndarray): # numpy backend - xest = lsqr(self, y, iter_lim=niter)[0] + xest = lsqr(self, y, iter_lim=niter, atol=1e-8, btol=1e-8)[0] else: # cupy backend ncp = get_array_module(y) diff --git a/pylops/basicoperators/BlockDiag.py b/pylops/basicoperators/BlockDiag.py index 4a13a18d..e29ee827 100644 --- a/pylops/basicoperators/BlockDiag.py +++ b/pylops/basicoperators/BlockDiag.py @@ -1,8 +1,16 @@ import multiprocessing as mp import numpy as np +import scipy as sp from scipy.sparse.linalg.interface import LinearOperator as spLinearOperator -from scipy.sparse.linalg.interface import _get_dtype + +# need to check scipy version since the interface submodule changed into +# _interface from scipy>=1.8.0 +sp_version = sp.__version__.split(".") +if int(sp_version[0]) <= 1 and int(sp_version[1]) < 8: + from scipy.sparse.linalg.interface import _get_dtype +else: + from scipy.sparse.linalg._interface import _get_dtype from pylops import LinearOperator from pylops.basicoperators import MatrixMult diff --git a/pylops/basicoperators/HStack.py b/pylops/basicoperators/HStack.py index 0f1a9b09..beeac03e 100644 --- a/pylops/basicoperators/HStack.py +++ b/pylops/basicoperators/HStack.py @@ -1,8 +1,16 @@ import multiprocessing as mp import numpy as np +import scipy as sp from scipy.sparse.linalg.interface import LinearOperator as spLinearOperator -from scipy.sparse.linalg.interface import _get_dtype + +# need to check scipy version since the interface submodule changed into +# _interface from scipy>=1.8.0 +sp_version = sp.__version__.split(".") +if int(sp_version[0]) <= 1 and int(sp_version[1]) < 8: + from scipy.sparse.linalg.interface import _get_dtype +else: + from scipy.sparse.linalg._interface import _get_dtype from pylops import LinearOperator from pylops.basicoperators import MatrixMult diff --git a/pylops/basicoperators/VStack.py b/pylops/basicoperators/VStack.py index aaba01db..248a4249 100644 --- a/pylops/basicoperators/VStack.py +++ b/pylops/basicoperators/VStack.py @@ -1,8 +1,16 @@ import multiprocessing as mp import numpy as np +import scipy as sp from scipy.sparse.linalg.interface import LinearOperator as spLinearOperator -from scipy.sparse.linalg.interface import _get_dtype + +# need to check scipy version since the interface submodule changed into +# _interface from scipy>=1.8.0 +sp_version = sp.__version__.split(".") +if int(sp_version[0]) <= 1 and int(sp_version[1]) < 8: + from scipy.sparse.linalg.interface import _get_dtype +else: + from scipy.sparse.linalg._interface import _get_dtype from pylops import LinearOperator from pylops.basicoperators import MatrixMult diff --git a/pylops/utils/describe.py b/pylops/utils/describe.py index 2fef5108..b9ea6b81 100644 --- a/pylops/utils/describe.py +++ b/pylops/utils/describe.py @@ -1,12 +1,25 @@ import random import string -from scipy.sparse.linalg.interface import ( - _AdjointLinearOperator, - _ProductLinearOperator, - _SumLinearOperator, - _TransposedLinearOperator, -) +import scipy as sp + +# need to check scipy version since the interface submodule changed into +# _interface from scipy>=1.8.0 +sp_version = sp.__version__.split(".") +if int(sp_version[0]) <= 1 and int(sp_version[1]) < 8: + from scipy.sparse.linalg.interface import ( + _AdjointLinearOperator, + _ProductLinearOperator, + _SumLinearOperator, + _TransposedLinearOperator, + ) +else: + from scipy.sparse.linalg._interface import ( + _AdjointLinearOperator, + _ProductLinearOperator, + _SumLinearOperator, + _TransposedLinearOperator, + ) from pylops import LinearOperator from pylops.basicoperators import BlockDiag, HStack, VStack diff --git a/pytests/test_avo.py b/pytests/test_avo.py index a50301e1..b66c5a1f 100755 --- a/pytests/test_avo.py +++ b/pytests/test_avo.py @@ -15,12 +15,14 @@ from pylops.avo.prestack import AVOLinearModelling from pylops.utils import dottest +np.random.seed(0) + # Create medium parameters for single contrast vp1, vs1, rho1 = 2200.0, 1300.0, 2000 # upper medium vp0, vs0, rho0 = 2300.0, 1400.0, 2100 # lower medium # Create medium parameters for multiple contrasts -nt0 = 501 +nt0 = 201 dt0 = 0.004 t0 = np.arange(nt0) * dt0 vp = 1200 + np.arange(nt0) + filtfilt(np.ones(5) / 5.0, 1, np.random.normal(0, 80, nt0)) @@ -126,5 +128,7 @@ def test_AVOLinearModelling(par): ) assert dottest(AVOop, ntheta * nt0, 3 * nt0) - minv = lsqr(AVOop, AVOop * m, damp=1e-20, iter_lim=1000, show=0)[0] + minv = lsqr( + AVOop, AVOop * m, damp=1e-20, iter_lim=1000, atol=1e-8, btol=1e-8, show=0 + )[0] assert_array_almost_equal(m, minv, decimal=3) diff --git a/pytests/test_basicoperators.py b/pytests/test_basicoperators.py index 827a6b3a..2e8fae2a 100755 --- a/pytests/test_basicoperators.py +++ b/pytests/test_basicoperators.py @@ -44,7 +44,9 @@ def test_Regression(par): assert dottest(LRop, par["ny"], order + 1) x = np.array([1.0, 2.0, 0.0, 3.0, -1.0], dtype=np.float32) - xlsqr = lsqr(LRop, LRop * x, damp=1e-10, iter_lim=300, show=0)[0] + xlsqr = lsqr( + LRop, LRop * x, damp=1e-10, iter_lim=300, atol=1e-8, btol=1e-8, show=0 + )[0] assert_array_almost_equal(x, xlsqr, decimal=3) y = LRop * x diff --git a/pytests/test_combine.py b/pytests/test_combine.py index 402745a2..2e4c197a 100755 --- a/pytests/test_combine.py +++ b/pytests/test_combine.py @@ -57,7 +57,9 @@ def test_VStack(par): Vop, 2 * par["ny"], par["nx"], complexflag=0 if par["imag"] == 0 else 3 ) - xlsqr = lsqr(Vop, Vop * x, damp=1e-20, iter_lim=300, show=0)[0] + xlsqr = lsqr(Vop, Vop * x, damp=1e-20, iter_lim=300, atol=1e-8, btol=1e-8, show=0)[ + 0 + ] assert_array_almost_equal(x, xlsqr, decimal=4) # use numpy matrix directly in the definition of the operator @@ -87,7 +89,9 @@ def test_HStack(par): Hop, par["ny"], 2 * par["nx"], complexflag=0 if par["imag"] == 0 else 3 ) - xlsqr = lsqr(Hop, Hop * x, damp=1e-20, iter_lim=300, show=0)[0] + xlsqr = lsqr(Hop, Hop * x, damp=1e-20, iter_lim=300, atol=1e-8, btol=1e-8, show=0)[ + 0 + ] assert_array_almost_equal(x, xlsqr, decimal=4) # use numpy matrix directly in the definition of the operator @@ -126,7 +130,9 @@ def test_Block(par): Bop, 2 * par["ny"], 2 * par["nx"], complexflag=0 if par["imag"] == 0 else 3 ) - xlsqr = lsqr(Bop, Bop * x, damp=1e-20, iter_lim=500, show=0)[0] + xlsqr = lsqr(Bop, Bop * x, damp=1e-20, iter_lim=500, atol=1e-8, btol=1e-8, show=0)[ + 0 + ] assert_array_almost_equal(x, xlsqr, decimal=3) # use numpy matrix directly in the definition of the operator @@ -171,7 +177,9 @@ def test_BlockDiag(par): BDop, 2 * par["ny"], 2 * par["nx"], complexflag=0 if par["imag"] == 0 else 3 ) - xlsqr = lsqr(BDop, BDop * x, damp=1e-20, iter_lim=500, show=0)[0] + xlsqr = lsqr( + BDop, BDop * x, damp=1e-20, iter_lim=500, atol=1e-8, btol=1e-8, show=0 + )[0] assert_array_almost_equal(x, xlsqr, decimal=3) # use numpy matrix directly in the definition of the operator diff --git a/pytests/test_wavedecomposition.py b/pytests/test_wavedecomposition.py index ca5198ca..6a0de6f5 100755 --- a/pytests/test_wavedecomposition.py +++ b/pytests/test_wavedecomposition.py @@ -123,7 +123,7 @@ def test_WavefieldDecomposition2D(par, create_data2D): ntaper=ntaper, dottest=True, dtype="complex128", - **dict(damp=1e-10, iter_lim=10) + **dict(damp=1e-10, atol=1e-8, btol=1e-8, iter_lim=10) ) assert np.linalg.norm(p2d_minus_est - p2d_minus) / np.linalg.norm(p2d_minus) < 2e-1 assert np.linalg.norm(p2d_plus_est - p2d_plus) / np.linalg.norm(p2d_plus) < 2e-1 @@ -169,7 +169,7 @@ def test_WavefieldDecomposition3D(par, create_data3D): ntaper=ntaper, dottest=True, dtype="complex128", - **dict(damp=1e-10, iter_lim=10, show=2) + **dict(damp=1e-10, iter_lim=10, atol=1e-8, btol=1e-8, show=2) ) assert np.linalg.norm(p3d_minus_est - p3d_minus) / np.linalg.norm(p3d_minus) < 3e-1 assert np.linalg.norm(p3d_plus_est - p3d_plus) / np.linalg.norm(p3d_plus) < 3e-1