Merge pull request #326 from mrava87/master

Bugfix: Modifications to handle scipy 1.8.0

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)

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

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

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

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

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)

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

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

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