Fix numpy 2 change and depreciated functions

.github/workflows/ci-ubuntu.yml

@@ -1,11 +1,13 @@
-# This workflow will install Python dependencies, run tests and lint with a variety of Python versions
+# This workflow will install Python dependencies, run tests with a variety of Python versions
 # For more information see: https://help.github.com/actions/language-and-framework-guides/using-python-with-github-actions
 
 name: CI-ubuntu
 
 on:
   push:
     branches: ['*']
+    paths-ignore:         # Skip CI in this case
+      - README.md    
   pull_request:
     branches: [ master ]
 
@@ -14,11 +16,11 @@ jobs:
     runs-on: ubuntu-latest
     strategy:
       matrix:
-        python-version: ["3.8", "3.9", "3.10", "3.11", "3.12"]
+        python-version: ["3.11", "3.12", "3.13", "3.14"]
     steps:
-    - uses: actions/checkout@v2
+    - uses: actions/checkout@main
     - name: Set up Python ${{ matrix.python-version }}
-      uses: actions/setup-python@v1
+      uses: actions/setup-python@main
       with:
         python-version: ${{ matrix.python-version }}
     - name: Install dependencies

.github/workflows/ci-windows.yml

@@ -11,8 +11,8 @@ jobs:
   test:
     runs-on: windows-latest
     steps:
-      - uses: actions/checkout@v2
-      - uses: conda-incubator/setup-miniconda@v2
+      - uses: actions/checkout@main
+      - uses: conda-incubator/setup-miniconda@main
         with:
           activate-environment: polze_env
           python-version: "3.11"

polze/_polze.py

@@ -71,7 +71,7 @@
 >>> z_ref = np.array([2, 3, 4])*np.pi
 >>> p.sort(); np.linalg.norm(p - p_ref) < 1e-10
 True
->>> z.sort(); np.linalg.norm(z - z_ref) < 1e-10
+>>> z.sort(); np.linalg.norm(z - z_ref) < 10e-10
 True
 
 The contour may be split into several annular regions if the upper bound for
@@ -177,15 +177,16 @@
 |";-:.
 ```
 """
+from collections import namedtuple
+from concurrent.futures import ProcessPoolExecutor
+import logging
+import sys
 import scipy.linalg as spl
 from scipy.spatial import KDTree
 import numpy as np
 import matplotlib.pyplot as plt
-from collections import namedtuple
-from concurrent.futures import ProcessPoolExecutor
+from numpy.lib import NumpyVersion
 
-import logging
-import sys
 
 # Create console logger handler
 _logging_level_default = logging.INFO
@@ -201,7 +202,14 @@
     logger.addHandler(ch)
 
 
-class PZ(object):
+# Define custom behavior depending on numpy version
+if NumpyVersion(np.__version__) < '2.0.0':
+    trapz = np.trapz
+else:
+    trapz = np.trapezoid
+
+
+class PZ():
     """Poles-zeros solver class.
 
     Attributes
@@ -226,17 +234,28 @@ class PZ(object):
         parameters (see below).
     """
 
-    _Npz_limit = 5;       """Number of Pole and Zeros max befor splitting path"""
-    _tol = 1e-5;          """General tolerance"""
-    _clean_tol = 1e-5;    """Spurious roots multiplicity tolerance"""
-    _NR_tol = 1e-12;      """Newton Raphson tolerance"""
-    _NiterMax  = 15;       """Max number of Newton Raphson iterations"""
-    _vectorized = False;  """Vectorized function evaluation"""
-    _RiShift = 1.02;      """Scaling factor of the radius when s0 is not an integer"""
-    _zeros_only = False;  """Use 0-moment for Npz estimate"""
-    _Niter_for_int = 5;   """Max Number of iterations to estimate the moments"""
-    _parallel = False;    """Split evaluation loop in a process pool executor"""
-    _max_workers = 1;     """Number of workers in the process pool executor"""
+    _Npz_limit = 5
+    """Number of Pole and Zeros max befor splitting path"""
+    _tol = 1e-5
+    """General tolerance"""
+    _clean_tol = 1e-5
+    """Spurious roots multiplicity tolerance"""
+    _NR_tol = 1e-12
+    """Newton Raphson tolerance"""
+    _NiterMax = 15
+    """Max number of Newton Raphson iterations"""
+    _vectorized = False
+    """Vectorized function evaluation"""
+    _RiShift = 1.02
+    """Scaling factor of the radius when s0 is not an integer"""
+    _zeros_only = False
+    """Use 0-moment for Npz estimate"""
+    _Niter_for_int = 5
+    """Max Number of iterations to estimate the moments"""
+    _parallel = False
+    """Split evaluation loop in a process pool executor"""
+    _max_workers = 1
+    """Number of workers in the process pool executor"""
 
     def __init__(self, fdf, Rmax, Npz, Ni=1024, R0=0., split=0, options=None):
 
@@ -252,7 +271,7 @@ def __init__(self, fdf, Rmax, Npz, Ni=1024, R0=0., split=0, options=None):
             self.df = None
 
         # dict for caching some data
-        self._h_cached = dict()
+        self._h_cached = {}
         self.Rmax = Rmax
 
         self.R0 = R0
@@ -311,7 +330,8 @@ def _eval(self, f, z):
             if self._parallel:
                 # Split the array for each worker
                 z_split = np.array_split(z, max_workers)
-                logger.debug('Parallel run with {} workers.'.format(max_workers))
+                logger.debug(
+                    f'Parallel run with {max_workers} workers.')
                 # array for storing each process output
                 fz_split = np.empty((max_workers,), dtype=object)
                 with ProcessPoolExecutor(max_workers=max_workers) as executor:
@@ -326,7 +346,8 @@ def _eval(self, f, z):
             fz = np.zeros(z.shape, dtype=complex)
             if self._parallel:
                 chunksize = z.size // max_workers
-                logger.debug('Parallel run with {} workers.'.format(max_workers))
+                logger.debug(
+                    f'Parallel run with {max_workers} workers.')
                 with ProcessPoolExecutor(max_workers=max_workers) as executor:
                     for i, fzi in enumerate(executor.map(f, z, chunksize=chunksize)):
                         fz[i] = fzi
@@ -361,11 +382,11 @@ def _eval_or_cache(self, Ri, z):
         """
 
         # Look for cached value
-        if Ri in self._h_cached.keys():
+        if Ri in self._h_cached:
             h = self._h_cached[Ri]
             # Check only if both have the same number of integration points
             if h.size == z.size:
-                logger.debug('Used cached value for Ri ={}'.format(Ri))
+                logger.debug(f'Used cached value for Ri ={Ri}')
                 return h
         else:
             # Compute f with loop or vectorized computation
@@ -417,12 +438,12 @@ def moment(self, K, Ri):
         z_rho = (z)/self._rho
         for k in range(0, K):
             # TODO try to limit numerical roundoff
-            s[k] = alpha * np.trapz(h * z_rho**k, z)
+            s[k] = alpha * trapz(h * z_rho**k, z)
 
         # Need modify the contour in this case.
         if abs(s[0]-np.round(s[0])) > self._tol:
-            logger.info('The first moment is not an integer (Ri={}, s0={}, tol={}).'.format(Ri, s[0], self._tol)
-                  + ' Some roots may lie on the integration path or increase Ni ({}).'.format(self.Ni))
+            logger.info(f'The first moment is not an integer (Ri={Ri}, s0={s[0]}, tol={self._tol}).'
+                        + f' Some roots may lie on the integration path or increase Ni ({self.Ni}).')
             return None
         else:
             return s.copy()
@@ -490,7 +511,8 @@ def estimate_Npz(self):
                     Ri = self._RiShift*Ri
                     self._Ri[i] = Ri
                     Niter += 1
-                    logger.info('  Try new radius Ri={} (Niter={})'.format(Ri, Niter))
+                    logger.info(
+                        f'  Try new radius Ri={Ri} (Niter={Niter})')
                     if Niter > Niter_for_int:
                         raise RuntimeError(('Max number of iteration is '
                                             'reached. The moment s0 is still '
@@ -542,7 +564,8 @@ def solve(self):
                     Ri = self._RiShift*Ri
                     self._Ri[i] = Ri
                     Niter += 1
-                    logger.info('  Try new radius Ri={} (Niter={})'.format(Ri, Niter))
+                    logger.info(
+                        f'  Try new radius Ri={Ri} (Niter={Niter})')
                     if Niter > Niter_for_int:
                         raise RuntimeError(('  Max number of iteration is '
                                             'reached. The moment s0 is still '
@@ -553,7 +576,8 @@ def solve(self):
                 s -= self._s[i-1]
             if np.round(s[0]) > Npz:
                 # necessary but not sufficient condition, ex Np ~ Nz -> s0~0
-                raise ValueError('The Npz upperbound ({}) if not sufficient (s0={})'.format(Npz, s[0]))
+                raise ValueError(
+                    f'The Npz upperbound ({Npz}) if not sufficient (s0={s[0]})')
             # build Hankel matrix
             H = spl.hankel(s[0:Npz], s[Npz-1:2*Npz-1])
             H2 = spl.hankel(s[1:Npz+1], s[Npz:2*Npz])
@@ -566,9 +590,7 @@ def solve(self):
                 D = D[~isinf]
                 NpzFinite = len(D)
                 logger.warning('Found inf in eig at ' +
-                               'Ri={}. Keep {} finite eigenvalue on {}.'.format(Ri,
-                                                                                NpzFinite,
-                                                                                Npz))
+                               f'Ri={Ri}. Keep {NpzFinite} finite eigenvalue on {Npz}.')
             else:
                 NpzFinite = Npz
 
@@ -619,11 +641,13 @@ def iterative_ref(self):
         f = self.f
         if self.df is None:
             logger.info('NR requires the derivative. Use FD approximation...')
+
             def df(z):
                 """ Five order finite difference approximation of df.
                 """
                 # eps = np.finfo(float).eps
-                h = 1e-4  # good compromize with cancellation error for f~1 with O(h**4)
+                # good compromize with cancellation error for f~1 with O(h**4)
+                h = 1e-4
                 # if f has zeros at z, Taylor series converge
                 if mu_ > 0:
                     fm2, f2 = f(z-2*h), f(z+2*h)
@@ -633,7 +657,7 @@ def df(z):
                 else:
                     # if f has pole at z, Taylor series DO NOT converge
                     # but 1/f Taylor series does...
-                    g = lambda x: 1/f(x)
+                    def g(x): return 1/f(x)
                     gm2, g2 = g(z-2*h), g(z+2*h)
                     g1, gm1 = g(z+h), g(z-h)
                     dgz = ((gm2 - g2)/12. + (g1 - gm1)*(2./3.))/h
@@ -672,7 +696,8 @@ def df(z):
 
         ErrTot = np.sum(np.abs(Err))
         if (np.abs(Err) > self._NR_tol).any():
-            logger.info('Refine step finalized. Some roots have not converged. Please check `info`.')
+            logger.info(
+                'Refine step finalized. Some roots have not converged. Please check `info`.')
 
         # package output
         info = {'Err': Err, 'Niter': Niter, 'ErrTot': ErrTot}
@@ -739,11 +764,12 @@ def contour_ref(self, Ni=500, scaling=0.01):
             pz_, mus_ = pz_ref.clean()
             # Expect only one value
             if pz_.size > 1:
-                logger.warning('Several roots have been found during contour refinement. Expect one.')
+                logger.warning(
+                    'Several roots have been found during contour refinement. Expect one.')
             # store solution and Error for check
-            pz_raf[n] = pz_
-            mus_raf[n] = mus_
-            Err[n] = abs(pz_ - pz)
+            pz_raf[n] = pz_.item()
+            mus_raf[n] = mus_.item()
+            Err[n] = abs(pz_ - pz).item()
 
         # Package output
         info = {'deviation': Err, 'radii': radii}
@@ -777,9 +803,7 @@ def clean(self, tol=None, split=False):
             dist_mult_to_int = np.abs(m[ind] - np.round(m[ind].real))
             if (dist_mult_to_int > tol).any():
                 logger.warning('Some multiplicities are far from integer ' +
-                               'at Ri={}. (Max dist = {}, tol = {}).'.format(self._Ri[i],
-                                                                             dist_mult_to_int.max(),
-                                                                             tol))
+                               f'at Ri={self._Ri[i]}. (Max dist = {dist_mult_to_int.max()}, tol = {tol}).')
             # TODO add a sort option
             pz.append(self.pz[i][ind])
             mu.append(m[ind])
@@ -804,9 +828,9 @@ def display(self, clean=True):
         print('-----------------------------------------------------------------------------------------------------------------------------')
         for i, (p, mu) in enumerate(zip(pz, m)):
             if np.round(mu) > 0:
-                print('                                             {}    {}'.format(p, mu))
+                print(f'                                             {p}    {mu}')
             else:
-                print('{}  |                                            {}'.format(p, mu))
+                print(f'{p}  |                                            {mu}')
         print('-----------------------------------------------------------------------------------------------------------------------------')
 
     def plot(self, ax=None, variable='\\nu'):
@@ -835,8 +859,10 @@ def plot(self, ax=None, variable='\\nu'):
 
         for pzi, mi in zip(pz, m):
             ind = mi.real > 0
-            ax.plot(pzi[ind].real, pzi[ind].imag, '.', color=col[col_index])  # zeros
-            ax.plot(pzi[~ind].real, pzi[~ind].imag, '+', color=col[col_index])  # poles
+            ax.plot(pzi[ind].real, pzi[ind].imag, '.',
+                    color=col[col_index])  # zeros
+            ax.plot(pzi[~ind].real, pzi[~ind].imag,
+                    '+', color=col[col_index])  # poles
             # color alternance
             col_index *= -1
 
@@ -877,7 +903,8 @@ def mapz(self, NgridPoint=51, variable='z', bounds=None, ax=None,
         """
 
         if bounds is None:
-            bounds = (-self.Rmax*(1+1j) + self.R0, self.Rmax * (1+1j) + self.R0)
+            bounds = (-self.Rmax*(1+1j) + self.R0,
+                      self.Rmax * (1+1j) + self.R0)
 
         z_real, z_imag = np.meshgrid(np.linspace(bounds[0].real, bounds[1].real, NgridPoint),
                                      np.linspace(bounds[0].imag, bounds[1].imag, NgridPoint))
@@ -886,7 +913,7 @@ def mapz(self, NgridPoint=51, variable='z', bounds=None, ax=None,
         FZ = self._eval(self.f, Z.ravel()).reshape(Z.shape)
 
         # plot
-        if not(ax):
+        if not ax:
             fig = plt.figure()
             ax = plt.gca()
         if part == 'modulus':