From beb7bbed5a43099a1b2f41077d794ea89efca18e Mon Sep 17 00:00:00 2001 From: Ben Greiner Date: Thu, 4 Nov 2021 20:16:55 +0100 Subject: [PATCH 01/10] try least-squares if regular solve returns nan --- control/statesp.py | 26 +++++++++++++++----------- 1 file changed, 15 insertions(+), 11 deletions(-) diff --git a/control/statesp.py b/control/statesp.py index f04bc55b3..25dc6be09 100644 --- a/control/statesp.py +++ b/control/statesp.py @@ -53,10 +53,10 @@ import math import numpy as np -from numpy import any, array, asarray, concatenate, cos, delete, \ - dot, empty, exp, eye, isinf, ones, pad, sin, zeros, squeeze, pi +from numpy import any, asarray, concatenate, cos, delete, \ + dot, empty, exp, eye, isinf, ones, pad, sin, zeros, squeeze from numpy.random import rand, randn -from numpy.linalg import solve, eigvals, matrix_rank +from numpy.linalg import eigvals, lstsq, matrix_rank, solve from numpy.linalg.linalg import LinAlgError import scipy as sp from scipy.signal import cont2discrete @@ -887,6 +887,8 @@ def horner(self, x, warn_infinite=True): """ # Make sure the argument is a 1D array of complex numbers x_arr = np.atleast_1d(x).astype(complex, copy=False) + if len(x_arr.shape) > 1: + raise ValueError("input list must be 1D") # return fast on systems with 0 or 1 state if not config.defaults['statesp.use_numpy_matrix']: @@ -908,10 +910,6 @@ def horner(self, x, warn_infinite=True): # Fall back because either Slycot unavailable or cannot handle # certain cases. - # Make sure that we are operating on a simple list - if len(x_arr.shape) > 1: - raise ValueError("input list must be 1D") - # Preallocate out = empty((self.noutputs, self.ninputs, len(x_arr)), dtype=complex) @@ -919,10 +917,14 @@ def horner(self, x, warn_infinite=True): # TODO: can this be vectorized? for idx, x_idx in enumerate(x_arr): try: - out[:, :, idx] = np.dot( - self.C, - solve(x_idx * eye(self.nstates) - self.A, self.B)) \ - + self.D + xI_A = x_idx * eye(self.nstates) - self.A + xI_A_inv = solve(xI_A, self.B) + # gh-664: not singular but the underlying LAPACK routine + # was not satisfied with the condition. Try least-squares + # solver. + if np.any(np.isnan(xI_A_inv)): # pragma: no cover + xI_A_inv, _, _, _ = lstsq(xI_A, self.B, rcond=None) + out[:, :, idx] = np.dot(self.C, xI_A_inv) + self.D except LinAlgError: # Issue a warning messsage, for consistency with xferfcn if warn_infinite: @@ -936,6 +938,8 @@ def horner(self, x, warn_infinite=True): else: out[:, :, idx] = complex(np.inf, np.nan) + + return out def freqresp(self, omega): From 88f62c4c6abd2c94e39de81daca87b9b78bde9e7 Mon Sep 17 00:00:00 2001 From: Ben Greiner Date: Thu, 4 Nov 2021 20:36:29 +0100 Subject: [PATCH 02/10] DEBUG --- control/statesp.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/control/statesp.py b/control/statesp.py index 25dc6be09..5f7f7eca1 100644 --- a/control/statesp.py +++ b/control/statesp.py @@ -919,10 +919,11 @@ def horner(self, x, warn_infinite=True): try: xI_A = x_idx * eye(self.nstates) - self.A xI_A_inv = solve(xI_A, self.B) - # gh-664: not singular but the underlying LAPACK routine - # was not satisfied with the condition. Try least-squares - # solver. + # gh-664: x_IA not singular but the underlying LAPACK + # routine was not satisfied with the condition. + # Try least-squares solver. if np.any(np.isnan(xI_A_inv)): # pragma: no cover + print(f"DEBUG: {xI_A}, {self.B}") xI_A_inv, _, _, _ = lstsq(xI_A, self.B, rcond=None) out[:, :, idx] = np.dot(self.C, xI_A_inv) + self.D except LinAlgError: From e0509261cc279665fadb9b67320db08ad3f46e4b Mon Sep 17 00:00:00 2001 From: Ben Greiner Date: Thu, 4 Nov 2021 20:47:27 +0100 Subject: [PATCH 03/10] DEBUG --- control/statesp.py | 13 +++++++++---- 1 file changed, 9 insertions(+), 4 deletions(-) diff --git a/control/statesp.py b/control/statesp.py index 5f7f7eca1..42c1e71ea 100644 --- a/control/statesp.py +++ b/control/statesp.py @@ -919,12 +919,17 @@ def horner(self, x, warn_infinite=True): try: xI_A = x_idx * eye(self.nstates) - self.A xI_A_inv = solve(xI_A, self.B) - # gh-664: x_IA not singular but the underlying LAPACK - # routine was not satisfied with the condition. - # Try least-squares solver. + # gh-664: xI_A did not raise singular matrix error, + # but the underlying LAPACK routine was not satisfied + # with the condition. + print(f"DEBUG np.linalg.solve(\n," + f"{xI_A}\n" + f"{self.B} = \n" + f"{xI_A_inv}") if np.any(np.isnan(xI_A_inv)): # pragma: no cover - print(f"DEBUG: {xI_A}, {self.B}") xI_A_inv, _, _, _ = lstsq(xI_A, self.B, rcond=None) + print(f"DEBUG np.linalg.lstsq(...) = \n," + f"{xI_A_inv}") out[:, :, idx] = np.dot(self.C, xI_A_inv) + self.D except LinAlgError: # Issue a warning messsage, for consistency with xferfcn From c22d031dd244ebe5c76444dae664354f0bac1173 Mon Sep 17 00:00:00 2001 From: Ben Greiner Date: Thu, 4 Nov 2021 20:56:04 +0100 Subject: [PATCH 04/10] DEBUG2 --- control/statesp.py | 6 ------ 1 file changed, 6 deletions(-) diff --git a/control/statesp.py b/control/statesp.py index 42c1e71ea..f7f68b77c 100644 --- a/control/statesp.py +++ b/control/statesp.py @@ -920,16 +920,10 @@ def horner(self, x, warn_infinite=True): xI_A = x_idx * eye(self.nstates) - self.A xI_A_inv = solve(xI_A, self.B) # gh-664: xI_A did not raise singular matrix error, - # but the underlying LAPACK routine was not satisfied - # with the condition. print(f"DEBUG np.linalg.solve(\n," f"{xI_A}\n" f"{self.B} = \n" f"{xI_A_inv}") - if np.any(np.isnan(xI_A_inv)): # pragma: no cover - xI_A_inv, _, _, _ = lstsq(xI_A, self.B, rcond=None) - print(f"DEBUG np.linalg.lstsq(...) = \n," - f"{xI_A_inv}") out[:, :, idx] = np.dot(self.C, xI_A_inv) + self.D except LinAlgError: # Issue a warning messsage, for consistency with xferfcn From fb56d92d6090073f546f4d6066a3919169c76557 Mon Sep 17 00:00:00 2001 From: Ben Greiner Date: Thu, 4 Nov 2021 21:05:50 +0100 Subject: [PATCH 05/10] remove DEBUG --- control/statesp.py | 5 ----- 1 file changed, 5 deletions(-) diff --git a/control/statesp.py b/control/statesp.py index f7f68b77c..eb0d3e015 100644 --- a/control/statesp.py +++ b/control/statesp.py @@ -919,11 +919,6 @@ def horner(self, x, warn_infinite=True): try: xI_A = x_idx * eye(self.nstates) - self.A xI_A_inv = solve(xI_A, self.B) - # gh-664: xI_A did not raise singular matrix error, - print(f"DEBUG np.linalg.solve(\n," - f"{xI_A}\n" - f"{self.B} = \n" - f"{xI_A_inv}") out[:, :, idx] = np.dot(self.C, xI_A_inv) + self.D except LinAlgError: # Issue a warning messsage, for consistency with xferfcn From 5939d9644fca9127e11be9ef73326b4d2c98b709 Mon Sep 17 00:00:00 2001 From: Ben Greiner Date: Thu, 4 Nov 2021 21:11:19 +0100 Subject: [PATCH 06/10] Revert "remove DEBUG" This reverts commit fb56d92d6090073f546f4d6066a3919169c76557. --- control/statesp.py | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/control/statesp.py b/control/statesp.py index eb0d3e015..f7f68b77c 100644 --- a/control/statesp.py +++ b/control/statesp.py @@ -919,6 +919,11 @@ def horner(self, x, warn_infinite=True): try: xI_A = x_idx * eye(self.nstates) - self.A xI_A_inv = solve(xI_A, self.B) + # gh-664: xI_A did not raise singular matrix error, + print(f"DEBUG np.linalg.solve(\n," + f"{xI_A}\n" + f"{self.B} = \n" + f"{xI_A_inv}") out[:, :, idx] = np.dot(self.C, xI_A_inv) + self.D except LinAlgError: # Issue a warning messsage, for consistency with xferfcn From 478fffb6cd58ce237ceeeb693124b5b537ea1625 Mon Sep 17 00:00:00 2001 From: Ben Greiner Date: Thu, 4 Nov 2021 21:48:47 +0100 Subject: [PATCH 07/10] no conda mkl in GHA --- .github/workflows/python-package-conda.yml | 3 +++ 1 file changed, 3 insertions(+) diff --git a/.github/workflows/python-package-conda.yml b/.github/workflows/python-package-conda.yml index 67f782048..0618ad862 100644 --- a/.github/workflows/python-package-conda.yml +++ b/.github/workflows/python-package-conda.yml @@ -28,6 +28,9 @@ jobs: conda create -q -n test-environment python=${{matrix.python-version}} source $CONDA/bin/activate test-environment + # gh-664, https://github.com/numpy/numpy/issues/20233, https://github.com/conda/conda/issues/11023 -- disable MKL + conda install nomkl + # Set up (virtual) X11 sudo apt install -y xvfb From cc44d3ccb047c43a568aadef4824310ac4a4962e Mon Sep 17 00:00:00 2001 From: Ben Greiner Date: Thu, 4 Nov 2021 21:56:26 +0100 Subject: [PATCH 08/10] show np config --- .github/workflows/python-package-conda.yml | 3 +++ 1 file changed, 3 insertions(+) diff --git a/.github/workflows/python-package-conda.yml b/.github/workflows/python-package-conda.yml index 0618ad862..2b6c008e6 100644 --- a/.github/workflows/python-package-conda.yml +++ b/.github/workflows/python-package-conda.yml @@ -52,6 +52,9 @@ jobs: # Use xvfb-run instead of pytest-xvfb to get proper mpl backend # Use coverage instead of pytest-cov to get .coverage file # See https://github.com/python-control/python-control/pull/504 + + python3 -c 'import numpy; numpy.show_config()' + xvfb-run --auto-servernum coverage run -m pytest control/tests - name: Coveralls parallel From 8757bdb5a39b75e919d59fab5ab10fd737510f90 Mon Sep 17 00:00:00 2001 From: Ben Greiner Date: Thu, 4 Nov 2021 22:09:34 +0100 Subject: [PATCH 09/10] cleanup --- control/statesp.py | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/control/statesp.py b/control/statesp.py index f7f68b77c..09abbb706 100644 --- a/control/statesp.py +++ b/control/statesp.py @@ -53,10 +53,10 @@ import math import numpy as np -from numpy import any, asarray, concatenate, cos, delete, \ - dot, empty, exp, eye, isinf, ones, pad, sin, zeros, squeeze +from numpy import any, array, asarray, concatenate, cos, delete, \ + dot, empty, exp, eye, isinf, ones, pad, sin, zeros, squeeze, pi from numpy.random import rand, randn -from numpy.linalg import eigvals, lstsq, matrix_rank, solve +from numpy.linalg import solve, eigvals, matrix_rank from numpy.linalg.linalg import LinAlgError import scipy as sp from scipy.signal import cont2discrete @@ -938,8 +938,6 @@ def horner(self, x, warn_infinite=True): else: out[:, :, idx] = complex(np.inf, np.nan) - - return out def freqresp(self, omega): From a9ec2dc48f88bffe88d1d70f4c0abe13d1d2acb7 Mon Sep 17 00:00:00 2001 From: Ben Greiner Date: Thu, 4 Nov 2021 22:19:22 +0100 Subject: [PATCH 10/10] Revert "no conda mkl in GHA" This reverts commit 478fffb6cd58ce237ceeeb693124b5b537ea1625. --- .github/workflows/python-package-conda.yml | 3 --- 1 file changed, 3 deletions(-) diff --git a/.github/workflows/python-package-conda.yml b/.github/workflows/python-package-conda.yml index 2b6c008e6..3b04949b9 100644 --- a/.github/workflows/python-package-conda.yml +++ b/.github/workflows/python-package-conda.yml @@ -28,9 +28,6 @@ jobs: conda create -q -n test-environment python=${{matrix.python-version}} source $CONDA/bin/activate test-environment - # gh-664, https://github.com/numpy/numpy/issues/20233, https://github.com/conda/conda/issues/11023 -- disable MKL - conda install nomkl - # Set up (virtual) X11 sudo apt install -y xvfb