Updating pseudopy to the recent python ecosystem evolutions

.github/workflows/ci-ubuntu.yml

@@ -0,0 +1,30 @@
+# This workflow will install Python dependencies and run tests with a variety of Python versions
+
+name: CI-ubuntu
+
+on:
+  push:
+    branches: ['*']
+    paths-ignore:         # Don't trigger on files that are updated by the CI
+      - README.md    
+  pull_request:
+    branches: [ master ]
+
+jobs:
+  build:
+    runs-on: ubuntu-latest
+    strategy:
+      matrix:
+        python-version: ["3.11", "3.12", "3.13", "3.14"]
+    steps:
+    - uses: actions/checkout@main
+    - name: Set up Python ${{ matrix.python-version }}
+      uses: actions/setup-python@main
+      with:
+        python-version: ${{ matrix.python-version }}
+    - name: Install the package
+      run: |       
+        pip install .[test]        
+    - name: Run test
+      run: |
+        python -m pseudopy.tests.test

.gitignore

@@ -1,5 +1,8 @@
-*.pyc
-.*.swp
-build
-dist
-MANIFEST
+# Byte-code
+__pycache__/
+*.py[cod]
+
+# Distribution / packaging
+*.egg-info/
+build/
+dist/

README.md

@@ -25,16 +25,21 @@ pyplot.show()
 ```
 
 ## Installation
-### Dependencies
-PseudoPy depends on numpy, scipy, matplotlib and shapely. If you are on Debian/Ubuntu, you can install these dependencies with
+
+To install the package from pypi, run
+
+`pip install pseudopy`
+
+It will install all dependencies like numpy, scipy, matplotlib and shapely. **Consider to use virtual environments to ensure better isolation**.
+To have the latest version, you can also use pip to install the package from the sources, available on the repository.
+
+To run the test suite, run
 ```
-sudo apt-get install python-numpy python-scipy python-matplotlib python-shapely
+pip install pseudopy[test]
+python -m pseudopy.tests.test
 ```
+The first line will install the extra dependencies required for tests.
 
-### pip
-```pip install pseudopy```
-
-Note that you may need to add `sudo` if you want to install it system-wide.
 
 ## License
 PseudoPy is free software licensed under the [MIT License](http://opensource.org/licenses/mit-license.php).

pseudopy/nonnormal.py

@@ -6,6 +6,27 @@
 
 from .utils import Path, Paths, plot_finish
 
+# Numpy 2 change inf name
+if numpy.lib.NumpyVersion(numpy.__version__) >= '2.0.0b1':
+    inf = numpy.inf
+else:
+    inf = numpy.Inf
+
+
+def contour_to_paths(contours):
+    """Extract the different paths from a contour.
+    """
+    paths = Paths()
+    # New matplotlib version (>3.8) flatten all paths in get_path
+    # the call of to_polygons split the list for each contour
+    for path in contours.get_paths():
+        # But by default to_polygons simplify the path
+        path.should_simplify = False
+        # polygons contains vertices
+        for p in path.to_polygons():
+            paths.append(Path(p[:, 0] + 1j*p[:, 1]))
+    return paths
+
 
 def inv_resolvent_norm(A, z, method='svd'):
     r'''Compute the reciprocal norm of the resolvent
@@ -38,10 +59,10 @@ def matvec(x):
             x1 = x[:m]
             x2 = x[m:]
             ret1 = AH.dot(x2) - numpy.conj(z)*x2
-            ret2 = numpy.array(A.dot(x1), dtype=numpy.complex)
+            ret2 = numpy.array(A.dot(x1), dtype=complex)
             ret2[:n] -= z*x1
             return numpy.c_[ret1, ret2]
-        AH_A = LinearOperator(matvec=matvec, dtype=numpy.complex,
+        AH_A = LinearOperator(matvec=matvec, dtype=complex,
                               shape=(m+n, m+n))
 
         evals = eigsh(AH_A, k=2, tol=1e-6, which='SM', maxiter=m+n+1,
@@ -59,7 +80,7 @@ def __init__(self, A, points, method='svd'):
         Stores result in self.vals and points in self.points
         '''
         self.points = points
-        if method == 'lanczosinv':
+        if method == 'lanczos':
             self.vals = []
 
             # algorithm from page 375 of Trefethen/Embree 2005
@@ -85,7 +106,7 @@ def matvec(x):
                         trans=2,
                         check_finite=False
                         )
-                MH_M = LinearOperator(matvec=matvec, dtype=numpy.complex,
+                MH_M = LinearOperator(matvec=matvec, dtype=complex,
                                       shape=(n, n))
 
                 evals = eigsh(MH_M, k=1, tol=1e-3, which='LM',
@@ -113,7 +134,7 @@ def __init__(self, A,
 
         # call super constructor
         super(NonnormalMeshgrid, self).__init__(
-            A, self.Real.flatten() + 1j*self.Imag.flatten())
+            A, self.Real.flatten() + 1j*self.Imag.flatten(), method=method)
         self.Vals = numpy.array(self.vals).reshape((imag_n, real_n))
 
     def plot(self, epsilons, **kwargs):
@@ -130,11 +151,9 @@ def contour_paths(self, epsilon):
         ax = figure.gca()
         contours = ax.contour(self.Real, self.Imag, self.Vals,
                               levels=[epsilon])
-        paths = Paths()
-        if len(contours.collections) == 0:
-            return paths
-        for path in contours.collections[0].get_paths():
-            paths.append(Path(path.vertices[:, 0] + 1j*path.vertices[:, 1]))
+        if len(contours.get_paths()) == 0:
+            return Paths()
+        paths = contour_to_paths(contours)
         pyplot.close(figure)
         return paths
 
@@ -154,11 +173,9 @@ def contour_paths(self, epsilon):
         '''Extract the polygon patches for the provided epsilon'''
         figure = pyplot.figure()
         contours = pyplot.tricontour(self.triang, self.vals, levels=[epsilon])
-        paths = Paths()
-        if len(contours.collections) == 0:
-            return paths
-        for path in contours.collections[0].get_paths():
-            paths.append(Path(path.vertices[:, 0] + 1j*path.vertices[:, 1]))
+        if len(contours.get_paths()) == 0:
+            return Paths()
+        paths = contour_to_paths(contours)
         pyplot.close(figure)
         return paths
 
@@ -272,12 +289,12 @@ def sort(lambd):
                     g_demmel1 = kappa = p + r_norm
 
                     # Demmel 2
-                    g_demmel2 = numpy.Inf
+                    g_demmel2 = inf
                     if radii[-1] <= sep_min/(2*kappa):
                         g_demmel2 = p + r_norm**2 * radii[-1]/(0.5*sep_min - p*radii[-1])
 
                     # Michael Karow bound (personal communication)
-                    g_mika = numpy.Inf
+                    g_mika = inf
                     if radii[-1] <= sep_min/(2*kappa):
                         eps_sep = radii[-1]/sep_min
                         g_mika = (p - eps_sep)/(
@@ -322,3 +339,4 @@ def sort(lambd):
                     points.append(midpoint + radius*numpy.exp(1j*(rand+arg)))
         points = numpy.concatenate(points)
         super(NonnormalAuto, self).__init__(A, points, **kwargs)
+

pseudopy/normal.py

@@ -1,7 +1,7 @@
 import numpy
 from matplotlib import pyplot
 import shapely.geometry as geom
-from shapely.ops import cascaded_union
+from shapely.ops import unary_union
 
 from .utils import get_paths, plot_finish
 
@@ -35,7 +35,7 @@ def contour_paths(self, epsilon):
                    .buffer(epsilon) for lamda in self.evals]
 
         # pseudospectrum is union of circles
-        pseudospec = cascaded_union(circles)
+        pseudospec = unary_union(circles)
 
         return get_paths(pseudospec)
 

pseudopy/tests/test.py

@@ -1,9 +1,11 @@
+from itertools import product
+import numpy
+import pseudopy
+import shapely.geometry as geom
+from shapely.ops import unary_union
 import matplotlib
 matplotlib.use('Agg')
 
-import numpy
-import pseudopy
-from itertools import product
 
 def dict_merge(*dicts):
     items = []
@@ -47,7 +49,7 @@ def test():
         pseudopy.NonnormalPoints: [{'A': A, 'points': points}],
         pseudopy.Normal: [{'A': A}],
         pseudopy.NormalEvals: [{'evals': evals}]
-        }
+    }
 
     # define epsilons
     epsilons = [0.2, 0.7, 1.1]
@@ -57,7 +59,7 @@ def test():
             pseudo = cls(**param)
 
             # test plot
-            #yield run_plot, pseudo, epsilons
+            # yield run_plot, pseudo, epsilons
 
             # test contour_paths
             for epsilon in epsilons:
@@ -76,16 +78,14 @@ def run_contour_paths(pseudo, epsilon, evals):
     paths = pseudo.contour_paths(epsilon)
 
     # check if pseudospectrum is correct by matching the parts of it
-    import shapely.geometry as geom
-    from shapely.ops import cascaded_union
     # create circles
     circles = [geom.Point(numpy.real(lamda), numpy.imag(lamda))
                .buffer(epsilon) for lamda in evals]
-    exact_pseudo = cascaded_union(circles)
+    exact_pseudo = unary_union(circles)
     exact_paths = pseudopy.utils.get_paths(exact_pseudo)
 
     N = len(paths)
-    assert(N == len(exact_paths))
+    assert (N == len(exact_paths))
 
     # create polygons
     polys = [geom.Polygon([(numpy.real(z), numpy.imag(z))
@@ -100,9 +100,9 @@ def run_contour_paths(pseudo, epsilon, evals):
     for (i, j) in product(range(N), range(N)):
         M[i, j] = exact_polys[i].symmetric_difference(polys[j]).area
     for i in range(N):
-        assert(numpy.min(M[i, :]) < 0.1*epsilon)
+        assert (numpy.min(M[i, :]) < 0.1*epsilon)
 
 
 if __name__ == '__main__':
-    import nose
-    nose.main()
+    import nose2
+    nose2.main()

pseudopy/utils.py

@@ -1,6 +1,7 @@
 import numpy
 from matplotlib import pyplot
-
+import matplotlib.ticker as ticker
+import shapely.geometry as geom
 
 class Path(object):
     def __init__(self, vertices):
@@ -33,11 +34,10 @@ def _get_points(c):
                 for sub in [polygon.exterior]+list(polygon.interiors)]
 
     paths = Paths()
-    import shapely.geometry as geom
     if isinstance(obj, geom.polygon.Polygon):
         paths += _get_polygon_paths(obj)
     elif isinstance(obj, geom.multipolygon.MultiPolygon):
-        for polygon in obj:
+        for polygon in obj.geoms:
             paths += _get_polygon_paths(polygon)
     return paths
 
@@ -47,16 +47,19 @@ def plot_finish(contours, spectrum=None, contour_labels=True, autofit=True):
     if spectrum is not None:
         pyplot.plot(numpy.real(spectrum), numpy.imag(spectrum), 'o')
 
+
     if autofit:
         vertices = []
-        for collection in contours.collections:
-            for path in collection.get_paths():
-                vertices.append(path.vertices[:, 0] + 1j*path.vertices[:, 1])
+        # After matplotlib 3.8
+        for path in contours.get_paths():
+            vertices.append(path.vertices[:, 0] + 1j*path.vertices[:, 1])
         vertices = numpy.concatenate(vertices)
         pyplot.xlim(numpy.min(vertices.real), numpy.max(vertices.real))
         pyplot.ylim(numpy.min(vertices.imag), numpy.max(vertices.imag))
 
     # plot contour labels?
-    from matplotlib.ticker import LogFormatterMathtext
     if contour_labels:
-        pyplot.clabel(contours, inline=1, fmt=LogFormatterMathtext())
+        fmt = ticker.LogFormatterMathtext()
+        fmt.create_dummy_axis()
+        pyplot.clabel(contours, inline=1, fmt=fmt)        
+

pyproject.toml

@@ -0,0 +1,38 @@
+[build-system]
+requires = ["hatchling"]
+build-backend = "hatchling.build"
+
+[project]
+name = 'pseudopy'
+version = '1.3.0'
+description = "Compute and visualize pseudospectra of matrices (like eigtool)"
+readme = 'README.md'
+license = {file = 'LICENSE'}
+requires-python = '>=3.11'  # PEP 518, shapely
+authors = [
+  {name = "André Gaul", email = 'gaul@web-yard.de'},
+]
+dependencies = [
+  'matplotlib>=3.8', 
+  'numpy>=1.23',
+  'scipy>=0.12',
+  'shapely>=2.0'
+]
+classifiers = [
+ 'Development Status :: 4 - Beta',
+ 'Intended Audience :: Science/Research',
+ 'License :: OSI Approved :: MIT License',
+ 'Operating System :: OS Independent',
+ 'Programming Language :: Python',
+ 'Programming Language :: Python :: 3',
+ 'Topic :: Scientific/Engineering :: Mathematics'
+]
+
+[project.urls]
+homepage = 'https://github.com/andrenarchy/pseudopy'
+
+[project.optional-dependencies]
+# Define a "test" group for test dependencies
+test = [
+    "nose2",
+]