diff --git a/.github/workflows/ci-ubuntu.yml b/.github/workflows/ci-ubuntu.yml index 6377771..f779ef8 100644 --- a/.github/workflows/ci-ubuntu.yml +++ b/.github/workflows/ci-ubuntu.yml @@ -1,4 +1,4 @@ -# 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 @@ -6,6 +6,8 @@ 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 diff --git a/.github/workflows/ci-windows.yml b/.github/workflows/ci-windows.yml index ae277a6..9399dc9 100644 --- a/.github/workflows/ci-windows.yml +++ b/.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" diff --git a/polze/_polze.py b/polze/_polze.py index ff77094..a1536fa 100755 --- a/polze/_polze.py +++ b/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': diff --git a/polze/test.py b/polze/test.py index aeafd21..fc21d9a 100644 --- a/polze/test.py +++ b/polze/test.py @@ -19,10 +19,10 @@ """ import unittest import doctest -import polze -import numpy as np import cmath import sys +import polze +import numpy as np # Numpy 2.0 change default printing options making doctest failing. # https://numpy.org/neps/nep-0051-scalar-representation.html @@ -35,7 +35,6 @@ def f(z): return cmath.tan(z) # define the function - def df(z): return cmath.tan(z)**2 + 1 # and its derivative [optional] @@ -44,15 +43,17 @@ def df(z): def f_np(z): return np.tan(z) + def df_np(z): return np.tan(z)**2 + 1 + class TestBasic(unittest.TestCase): - """ Test suite. + """Test suite. """ def test_non_vectorized(self): - """ Test with non vectorized function. + """Test with non vectorized function. """ # Reuse a non vectorized function pz = polze.PZ((f, df), Rmax=5, Npz=10, Ni=1024, @@ -69,14 +70,12 @@ def test_non_vectorized(self): self.assertTrue(np.linalg.norm(z - z_ref) < 1e-12) def test_iterative_ref(self): - """ Test iterative refinement for simple roots. + """Test iterative refinement for simple roots. """ tol = 1e-12 # Define a vectorized function - f = lambda z: np.tan(z) # define the function - df = lambda z: np.tan(z)**2 + 1 # and its derivative [optional] # Low quality initial solution (small Ni) - pz = polze.PZ((f, df), Rmax=5, Npz=10, Ni=100, + pz = polze.PZ((f_np, df_np), Rmax=5, Npz=10, Ni=100, options={'_vectorized': True, '_tol': 1e-3, '_NR_tol': tol}) @@ -101,15 +100,13 @@ def test_iterative_ref(self): self.assertFalse(np.linalg.norm(z0 - z_ref) < tol) def test_iterative_ref_no_df(self): - """ Test iterative refinement for simple roots without df. + """Test iterative refinement for simple roots without df. """ tol = 1e-10 - # Define a vectorized function - f = lambda z: np.tan(z) # define the function # Low quality initial solution (small Ni) # Here we prescribe a big '_clean_tol' to avoid spurious root that # break the number of poles in the test - pz = polze.PZ(f, Rmax=5, Npz=10, Ni=100, + pz = polze.PZ(f_np, Rmax=5, Npz=10, Ni=100, options={'_vectorized': True, '_tol': 1e-3, '_NR_tol': tol, @@ -133,16 +130,13 @@ def test_iterative_ref_no_df(self): self.assertTrue(info['ErrTot'] < tol) # initial sols are not accurate enougth self.assertFalse(np.linalg.norm(z0 - z_ref) < tol) - + def test_contour_ref(self): - """ Test contour refinement for simple roots. + """Test contour refinement for simple roots. """ tol = 1e-12 - # Define a vectorized function - f = lambda z: np.tan(z) # define the function - df = lambda z: np.tan(z)**2 + 1 # and its derivative [optional] # Low quality initial solution (small Ni) - pz = polze.PZ((f, df), Rmax=5, Npz=10, Ni=100, + pz = polze.PZ((f_np, df_np), Rmax=5, Npz=10, Ni=100, options={'_vectorized': True, '_tol': 1e-3, '_NR_tol': tol}) @@ -165,14 +159,11 @@ def test_contour_ref(self): self.assertFalse(np.linalg.norm(z0 - z_ref) < tol) def test_with_finite_difference_on_circle(self): - """ Test using the finite difference on the circle for a + """Test using the finite difference on the circle for a non-vectorized function. """ # Here tol bigger than when the derivative are given. test_tol = 1e-9 - # define a non vectorized function - f = lambda z: cmath.tan(z) - pz = polze.PZ(f, Rmax=5, Npz=10, Ni=1024, options={'_vectorized': False}) _ = pz.solve() @@ -187,13 +178,13 @@ def test_with_finite_difference_on_circle(self): self.assertTrue(np.linalg.norm(z - z_ref) < test_tol) def test_zeros_only(self): - """ Test with function without poles. + """Test with function without poles. """ # define a non vectorized function - f = lambda z: np.sin(z) # define the function - df = lambda z: np.cos(z) # and its derivative [optional] + def f2(z): return np.sin(z) # define the function + def df2(z): return np.cos(z) # and its derivative [optional] # Npz is too small and will be corrected by `estimate_Npz` method - pz = polze.PZ((f, df), Rmax=5, Npz=1, Ni=2048, + pz = polze.PZ((f2, df2), Rmax=5, Npz=1, Ni=2048, options={'_zeros_only': True}) _ = pz.solve() p, z = pz.dispatch() @@ -206,7 +197,7 @@ def test_zeros_only(self): self.assertTrue(np.linalg.norm(z - z_ref) < 1e-12) def test_rational_fraction(self): - """ Test with a rational fraction. + """Test with a rational fraction. """ tol = 1e-10 # Define the rational fraction @@ -214,14 +205,14 @@ def test_rational_fraction(self): z_ref = np.array([1. + 0.2j, 2. - 0.3j, -1 + 0.4j, 7.]) p_ref = np.array([.5 + 0.2j, 3 + 0.3j, -7, 5.]) - def f(z): - """ Define the rational fraction from roots. + def fr(z): + """Define the rational fraction from roots. """ n = pol.polyval(z, pol.polyfromroots(z_ref)) d = pol.polyval(z, pol.polyfromroots(p_ref)) return n/d - pz = polze.PZ(f, Rmax=8, Npz=10, Ni=2048) + pz = polze.PZ(fr, Rmax=8, Npz=10, Ni=2048) _ = pz.solve() p, z = pz.dispatch() z.sort() @@ -238,14 +229,14 @@ def test_rational_fraction_multiple_roots(self): z_ref = np.array([1. + 0.2j, 1. + 0.2j, -1 + 0.4j, 7.]) p_ref = np.array([.5 + 0.2j, -7, -7, -7]) - def f(z): + def fr(z): """ Define the rational fraction from roots. """ n = pol.polyval(z, pol.polyfromroots(z_ref)) d = pol.polyval(z, pol.polyfromroots(p_ref)) return n/d - pz = polze.PZ(f, Rmax=8, Npz=10, Ni=2048) + pz = polze.PZ(fr, Rmax=8, Npz=10, Ni=2048) _ = pz.solve() (p, pm), (z, zm) = pz.dispatch(multiplicities=True) z.sort() @@ -278,7 +269,6 @@ def test_parallel(self): z.sort() self.assertTrue(np.linalg.norm(z - z_ref) < 1e-12) - def test_parallel_vectorized(self): """ Test accuracy with parallel execution with vectorized (numpy) func. @@ -308,7 +298,7 @@ def test_parallel_vectorized(self): suite = loadTestsFromTestCase(TestBasic) # Add doctest from _polze module suite.addTest(doctest.DocTestSuite(polze._polze, - optionflags=(doctest.ELLIPSIS|doctest.NORMALIZE_WHITESPACE))) + optionflags=doctest.ELLIPSIS | doctest.NORMALIZE_WHITESPACE)) # define the runner runner = unittest.TextTestRunner(verbosity=3) # run all the suite