diff --git a/slycot/tests/CMakeLists.txt b/slycot/tests/CMakeLists.txt index 80d751da..e397b0de 100644 --- a/slycot/tests/CMakeLists.txt +++ b/slycot/tests/CMakeLists.txt @@ -14,8 +14,16 @@ set(PYSOURCE test_transform.py test_sg02ad.py test_sg03ad.py + test_tb01id.py + test_tb01pd.py + test_tb03ad.py + test_tb04ad.py test_tb05ad.py + test_tc01od.py + test_tc04ad.py test_td04ad.py + test_tf01md.py + test_tf01rd.py test_tg01ad.py test_tg01fd.py ) diff --git a/slycot/tests/test_tb01id.py b/slycot/tests/test_tb01id.py new file mode 100644 index 00000000..efb6bcb7 --- /dev/null +++ b/slycot/tests/test_tb01id.py @@ -0,0 +1,15 @@ +import warnings + +import pytest + +from slycot import transform as tf +from slycot.exceptions import SlycotArithmeticError + +from .test_exceptions import assert_docstring_parse + + +class Test_tb01id: + + @pytest.mark.skip(reason="not implemented") + def test_tb01id(self): + pass diff --git a/slycot/tests/test_tb01pd.py b/slycot/tests/test_tb01pd.py new file mode 100644 index 00000000..f4018b8a --- /dev/null +++ b/slycot/tests/test_tb01pd.py @@ -0,0 +1,15 @@ +import warnings + +import pytest + +from slycot import transform as tf +from slycot.exceptions import SlycotArithmeticError + +from .test_exceptions import assert_docstring_parse + + +class Test_tb01pd: + + @pytest.mark.skip(reason="not implemented") + def test_tb01pd(self): + pass \ No newline at end of file diff --git a/slycot/tests/test_tb03ad.py b/slycot/tests/test_tb03ad.py new file mode 100644 index 00000000..a9cb6024 --- /dev/null +++ b/slycot/tests/test_tb03ad.py @@ -0,0 +1,21 @@ +import warnings + +import pytest + +from slycot import transform as tf +from slycot.exceptions import SlycotArithmeticError + +from .test_exceptions import assert_docstring_parse + + +class Test_tb03ad: + + @pytest.mark.parametrize( + 'fun, exception_class, erange, checkvars', + ((tf.tb03ad, SlycotArithmeticError, 2, {}),)) + def test_tb03ad_docparse(self, fun, exception_class, erange, checkvars): + assert_docstring_parse(fun.__doc__, exception_class, erange, checkvars) + + @pytest.mark.skip(reason="not implemented") + def test_tb03ad(self): + pass \ No newline at end of file diff --git a/slycot/tests/test_tb04ad.py b/slycot/tests/test_tb04ad.py new file mode 100644 index 00000000..7d241333 --- /dev/null +++ b/slycot/tests/test_tb04ad.py @@ -0,0 +1,15 @@ +import warnings + +import pytest + +from slycot import transform as tf +from slycot.exceptions import SlycotArithmeticError + +from .test_exceptions import assert_docstring_parse + + +class Test_tb04ad: + + @pytest.mark.skip(reason="not implemented") + def test_tb04ad(self): + pass \ No newline at end of file diff --git a/slycot/tests/test_tb05ad.py b/slycot/tests/test_tb05ad.py index 900a49a2..df5996d1 100644 --- a/slycot/tests/test_tb05ad.py +++ b/slycot/tests/test_tb05ad.py @@ -1,6 +1,3 @@ -# =================================================== -# tb05ad tests - import sys import numpy as np @@ -9,8 +6,11 @@ from scipy.linalg import eig, matrix_balance from slycot import transform +from slycot import transform as tf from slycot.exceptions import SlycotArithmeticError, SlycotParameterError +from .test_exceptions import assert_docstring_parse + # set the random seed so we can get consistent results. np.random.seed(40) CASES = {} @@ -36,216 +36,226 @@ 'B': np.random.randn(n, m), 'C': np.random.randn(p, n)} - -def test_tb05ad_ng(): - """ - Test that tb05ad with job 'NG' computes the correct - frequency response. - """ - for key in CASES: - sys = CASES[key] - check_tb05ad_AG_NG(sys, 10*1j, 'NG') - - -def test_tb05ad_ag(): - """ - Test that tb05ad with job 'AG' computes the correct - frequency response. - """ - for key in CASES: - sys = CASES[key] - check_tb05ad_AG_NG(sys, 10*1j, 'AG') - - -def test_tb05ad_nh(): - """Test that tb05ad with job = 'NH' computes the correct - frequency response after conversion to Hessenberg form. - - First call tb05ad with job='NH' to transform to upper Hessenberg - form which outputs the transformed system. - Subsequently, call tb05ad with job='NH' using this transformed system. - """ - jomega = 10*1j - for key in CASES: - sys = CASES[key] - sys_transformed = check_tb05ad_AG_NG(sys, jomega, 'NG') - check_tb05ad_NH(sys_transformed, sys, jomega) - - -def test_tb05ad_errors(): - """ - Test tb05ad error handling. We give wrong inputs and - and check that this raises an error. - """ - check_tb05ad_errors(CASES['pass1']) - - -def check_tb05ad_AG_NG(sys, jomega, job): - """ - Check that tb05ad computes the correct frequency response when - running jobs 'AG' and/or 'NG'. - - Inputs - ------ - - sys: A a dict of system matrices with keys 'A', 'B', and 'C'. - jomega: A complex scalar, which is the frequency we are - evaluating the system at. - job: A string, either 'AG' or 'NH' - - Returns - ------- - sys_transformed: A dict of the system matrices which have been - transformed according the job. - """ - n, m = sys['B'].shape - p = sys['C'].shape[0] - result = transform.tb05ad(n, m, p, jomega, - sys['A'], sys['B'], sys['C'], job=job) - g_i = result[3] - hinvb = np.linalg.solve(np.eye(n) * jomega - sys['A'], sys['B']) - g_i_solve = sys['C'].dot(hinvb) - assert_almost_equal(g_i_solve, g_i) - sys_transformed = {'A': result[0], 'B': result[1], 'C': result[2]} - return sys_transformed - - -def check_tb05ad_NH(sys_transformed, sys, jomega): - """ - Check tb05ad, computes the correct frequency response when - job='NH' and we supply system matrices 'A', 'B', and 'C' - which have been transformed by a previous call to tb05ad. - We check we get the same result as computing C(sI - A)^-1B - with the original system. - - Inputs - ------ - - sys_transformed: A a dict of the transformed (A in upper - hessenberg form) system matrices with keys - 'A', 'B', and 'C'. - - sys: A dict of the original un-transformed system matrices. - - jomega: A complex scalar, which is the frequency to evaluate at. - - """ - - n, m = sys_transformed['B'].shape - p = sys_transformed['C'].shape[0] - result = transform.tb05ad(n, m, p, jomega, sys_transformed['A'], - sys_transformed['B'], sys_transformed['C'], - job='NH') - g_i = result[0] - hinvb = np.linalg.solve(np.eye(n) * jomega - sys['A'], sys['B']) - g_i_solve = sys['C'].dot(hinvb) - assert_almost_equal(g_i_solve, g_i) - - -def check_tb05ad_errors(sys): - """ - Check the error handling of tb05ad. We give wrong inputs and - and check that this raises an error. - """ - n, m = sys['B'].shape - p = sys['C'].shape[0] - jomega = 10*1j - # test error handling - # wrong size A - with pytest.raises(SlycotParameterError) as cm: - transform.tb05ad( - n+1, m, p, jomega, sys['A'], sys['B'], sys['C'], job='NH') - assert cm.value.info == -7 - # wrong size B - with pytest.raises(SlycotParameterError) as cm: - transform.tb05ad( - n, m+1, p, jomega, sys['A'], sys['B'], sys['C'], job='NH') - assert cm.value.info == -9 - # wrong size C - with pytest.raises(SlycotParameterError) as cm: - transform.tb05ad( - n, m, p+1, jomega, sys['A'], sys['B'], sys['C'], job='NH') - assert cm.value.info == -11 - # unrecognized job - with pytest.raises(SlycotParameterError) as cm: - transform.tb05ad( - n, m, p, jomega, sys['A'], sys['B'], sys['C'], job='a') - assert cm.value.info == -1 - - -def test_tb05ad_resonance(): - """ Test tb05ad resonance failure. - - Actually test one of the exception messages. These - are parsed from the docstring, tests both the info index and the - message - """ - A = np.array([[0, -1], - [1, 0]]) - B = np.array([[1], - [0]]) - C = np.array([[0, 1]]) - jomega = 1j - with pytest.raises( - SlycotArithmeticError, - match=r"Either `freq`.* is too near to an eigenvalue of A,\n" - r"or `rcond` is less than the machine precision EPS.") as cm: - transform.tb05ad(2, 1, 1, jomega, A, B, C, job='NH') - assert cm.value.info == 2 - - -def test_tb05ad_balance(): - """Test balancing in tb05ad. - - Tests for the cause of the problem reported in issue #11 - balancing permutations were not correctly applied to the - C and D matrix. - """ - - # find a good test case. Some sparsity, - # some zero eigenvalues, some non-zero eigenvalues, - # and proof that the 1st step, with dgebal, does some - # permutation and some scaling - crit = False - n = 8 - while not crit: - A = np.random.randn(n, n) - A[np.random.uniform(size=(n, n)) > 0.35] = 0.0 - - Aeig = eig(A)[0] - neig0 = np.sum(np.abs(Aeig) == 0) - As, T = matrix_balance(A) - nperm = np.sum(np.diag(T == 0)) - nscale = n - np.sum(T == 1.0) - crit = nperm < n and nperm >= n//2 and \ - neig0 > 1 and neig0 <= 3 and nscale > 0 - - # print("number of permutations", nperm, "eigenvalues=0", neig0) - B = np.random.randn(8, 4) - C = np.random.randn(3, 8) - - # do a run - jomega = 1.0 - At, Bt, Ct, rcond, g_jw, ev, hinvb, info = transform.tb05ad( - 8, 4, 3, jomega, A, B, C, job='AG') - - # remove information on Q, in lower sub-triangle part of A - At = np.triu(At, k=-1) - - # now after the balancing in DGEBAL, and conversion to - # upper Hessenberg form: - # At = Q^T * (P^-1 * A * P ) * Q - # with Q orthogonal - # Ct = C * P * Q - # Bt = Q^T * P^-1 * B - # so test with Ct * At * Bt == C * A * B - # and verify that eigenvalues of both A matrices are close - assert_almost_equal(np.dot(np.dot(Ct, At), Bt), - np.dot(np.dot(C, A), B)) - # uses a sort, there is no guarantee on the order of eigenvalues - eigAt = eig(At)[0] - idxAt = np.argsort(eigAt) - eigA = eig(A)[0] - idxA = np.argsort(eigA) - assert_almost_equal(eigA[idxA], eigAt[idxAt]) +class Test_tb05ad: + + @pytest.mark.parametrize( + 'fun, exception_class, erange, checkvars', + ((tf.tb05ad, SlycotArithmeticError, 2, {'n30': 90, + 'jomega': 2.0, + 'rcond': 1e-12}),)) + def test_tb05ad_docparse(self, fun, exception_class, erange, checkvars): + assert_docstring_parse(fun.__doc__, exception_class, erange, checkvars) + + + def test_tb05ad_ng(self): + """ + Test that tb05ad with job 'NG' computes the correct + frequency response. + """ + for key in CASES: + sys = CASES[key] + self.check_tb05ad_AG_NG(sys, 10*1j, 'NG') + + + def test_tb05ad_ag(self): + """ + Test that tb05ad with job 'AG' computes the correct + frequency response. + """ + for key in CASES: + sys = CASES[key] + self.check_tb05ad_AG_NG(sys, 10*1j, 'AG') + + + def test_tb05ad_nh(self): + """Test that tb05ad with job = 'NH' computes the correct + frequency response after conversion to Hessenberg form. + + First call tb05ad with job='NH' to transform to upper Hessenberg + form which outputs the transformed system. + Subsequently, call tb05ad with job='NH' using this transformed system. + """ + jomega = 10*1j + for key in CASES: + sys = CASES[key] + sys_transformed = self.check_tb05ad_AG_NG(sys, jomega, 'NG') + self.check_tb05ad_NH(sys_transformed, sys, jomega) + + + def test_tb05ad_errors(self): + """ + Test tb05ad error handling. We give wrong inputs and + and check that this raises an error. + """ + self.check_tb05ad_errors(CASES['pass1']) + + + def check_tb05ad_AG_NG(self, sys, jomega, job): + """ + Check that tb05ad computes the correct frequency response when + running jobs 'AG' and/or 'NG'. + + Inputs + ------ + + sys: A a dict of system matrices with keys 'A', 'B', and 'C'. + jomega: A complex scalar, which is the frequency we are + evaluating the system at. + job: A string, either 'AG' or 'NH' + + Returns + ------- + sys_transformed: A dict of the system matrices which have been + transformed according the job. + """ + n, m = sys['B'].shape + p = sys['C'].shape[0] + result = transform.tb05ad(n, m, p, jomega, + sys['A'], sys['B'], sys['C'], job=job) + g_i = result[3] + hinvb = np.linalg.solve(np.eye(n) * jomega - sys['A'], sys['B']) + g_i_solve = sys['C'].dot(hinvb) + assert_almost_equal(g_i_solve, g_i) + sys_transformed = {'A': result[0], 'B': result[1], 'C': result[2]} + return sys_transformed + + + def check_tb05ad_NH(self, sys_transformed, sys, jomega): + """ + Check tb05ad, computes the correct frequency response when + job='NH' and we supply system matrices 'A', 'B', and 'C' + which have been transformed by a previous call to tb05ad. + We check we get the same result as computing C(sI - A)^-1B + with the original system. + + Inputs + ------ + + sys_transformed: A a dict of the transformed (A in upper + hessenberg form) system matrices with keys + 'A', 'B', and 'C'. + + sys: A dict of the original un-transformed system matrices. + + jomega: A complex scalar, which is the frequency to evaluate at. + + """ + + n, m = sys_transformed['B'].shape + p = sys_transformed['C'].shape[0] + result = transform.tb05ad(n, m, p, jomega, sys_transformed['A'], + sys_transformed['B'], sys_transformed['C'], + job='NH') + g_i = result[0] + hinvb = np.linalg.solve(np.eye(n) * jomega - sys['A'], sys['B']) + g_i_solve = sys['C'].dot(hinvb) + assert_almost_equal(g_i_solve, g_i) + + + def check_tb05ad_errors(self, sys): + """ + Check the error handling of tb05ad. We give wrong inputs and + and check that this raises an error. + """ + n, m = sys['B'].shape + p = sys['C'].shape[0] + jomega = 10*1j + # test error handling + # wrong size A + with pytest.raises(SlycotParameterError) as cm: + transform.tb05ad( + n+1, m, p, jomega, sys['A'], sys['B'], sys['C'], job='NH') + assert cm.value.info == -7 + # wrong size B + with pytest.raises(SlycotParameterError) as cm: + transform.tb05ad( + n, m+1, p, jomega, sys['A'], sys['B'], sys['C'], job='NH') + assert cm.value.info == -9 + # wrong size C + with pytest.raises(SlycotParameterError) as cm: + transform.tb05ad( + n, m, p+1, jomega, sys['A'], sys['B'], sys['C'], job='NH') + assert cm.value.info == -11 + # unrecognized job + with pytest.raises(SlycotParameterError) as cm: + transform.tb05ad( + n, m, p, jomega, sys['A'], sys['B'], sys['C'], job='a') + assert cm.value.info == -1 + + + def test_tb05ad_resonance(self): + """ Test tb05ad resonance failure. + + Actually test one of the exception messages. These + are parsed from the docstring, tests both the info index and the + message + """ + A = np.array([[0, -1], + [1, 0]]) + B = np.array([[1], + [0]]) + C = np.array([[0, 1]]) + jomega = 1j + with pytest.raises( + SlycotArithmeticError, + match=r"Either `freq`.* is too near to an eigenvalue of A,\n" + r"or `rcond` is less than the machine precision EPS.") as cm: + transform.tb05ad(2, 1, 1, jomega, A, B, C, job='NH') + assert cm.value.info == 2 + + + def test_tb05ad_balance(self): + """Test balancing in tb05ad. + + Tests for the cause of the problem reported in issue #11 + balancing permutations were not correctly applied to the + C and D matrix. + """ + + # find a good test case. Some sparsity, + # some zero eigenvalues, some non-zero eigenvalues, + # and proof that the 1st step, with dgebal, does some + # permutation and some scaling + crit = False + n = 8 + while not crit: + A = np.random.randn(n, n) + A[np.random.uniform(size=(n, n)) > 0.35] = 0.0 + + Aeig = eig(A)[0] + neig0 = np.sum(np.abs(Aeig) == 0) + As, T = matrix_balance(A) + nperm = np.sum(np.diag(T == 0)) + nscale = n - np.sum(T == 1.0) + crit = nperm < n and nperm >= n//2 and \ + neig0 > 1 and neig0 <= 3 and nscale > 0 + + # print("number of permutations", nperm, "eigenvalues=0", neig0) + B = np.random.randn(8, 4) + C = np.random.randn(3, 8) + + # do a run + jomega = 1.0 + At, Bt, Ct, rcond, g_jw, ev, hinvb, info = transform.tb05ad( + 8, 4, 3, jomega, A, B, C, job='AG') + + # remove information on Q, in lower sub-triangle part of A + At = np.triu(At, k=-1) + + # now after the balancing in DGEBAL, and conversion to + # upper Hessenberg form: + # At = Q^T * (P^-1 * A * P ) * Q + # with Q orthogonal + # Ct = C * P * Q + # Bt = Q^T * P^-1 * B + # so test with Ct * At * Bt == C * A * B + # and verify that eigenvalues of both A matrices are close + assert_almost_equal(np.dot(np.dot(Ct, At), Bt), + np.dot(np.dot(C, A), B)) + # uses a sort, there is no guarantee on the order of eigenvalues + eigAt = eig(At)[0] + idxAt = np.argsort(eigAt) + eigA = eig(A)[0] + idxA = np.argsort(eigA) + assert_almost_equal(eigA[idxA], eigAt[idxAt]) diff --git a/slycot/tests/test_tc01od.py b/slycot/tests/test_tc01od.py new file mode 100644 index 00000000..9b246ed7 --- /dev/null +++ b/slycot/tests/test_tc01od.py @@ -0,0 +1,15 @@ +import warnings + +import pytest + +from slycot import transform as tf +from slycot.exceptions import SlycotArithmeticError + +from .test_exceptions import assert_docstring_parse + + +class Test_tc01od: + + @pytest.mark.skip(reason="not implemented") + def test_tc01od(self): + pass \ No newline at end of file diff --git a/slycot/tests/test_tc04ad.py b/slycot/tests/test_tc04ad.py new file mode 100644 index 00000000..c8aa4fb0 --- /dev/null +++ b/slycot/tests/test_tc04ad.py @@ -0,0 +1,21 @@ +import warnings + +import pytest + +from slycot import transform as tf +from slycot.exceptions import SlycotArithmeticError + +from .test_exceptions import assert_docstring_parse + + +class Test_tc04ad: + + @pytest.mark.parametrize( + 'fun, exception_class, erange, checkvars', + ((tf.tc04ad, SlycotArithmeticError, 1, {'leri': 'L'}),)) + def test_tc04ad_docparse(self, fun, exception_class, erange, checkvars): + assert_docstring_parse(fun.__doc__, exception_class, erange, checkvars) + + @pytest.mark.skip(reason="not implemented") + def test_tc04ad(self): + pass \ No newline at end of file diff --git a/slycot/tests/test_td04ad.py b/slycot/tests/test_td04ad.py index 57c557d0..c1947a6a 100644 --- a/slycot/tests/test_td04ad.py +++ b/slycot/tests/test_td04ad.py @@ -1,220 +1,230 @@ -# -# test_td04ad.py - test suite for tf -> ss conversion -# RvP, 04 Jun 2018 - import numpy as np +import pytest from slycot import transform - - -def test_td04ad_c(): - """td04ad: Convert with 'C' option""" - - # for octave: - """ - num = { [0.0, 0.0, 1.0 ], [ 1.0, 0.0 ]; - [3.0, -1.0, 1.0 ], [ 0.0, 1.0 ]; - [0.0, 0.0, 1.0], [ 0.0, 2.0 ] }; - den = { [1.0, 0.4, 3.0], [ 1.0, 1.0 ]; - [1.0, 0.4, 3.0], [ 1.0, 1.0 ]; - [1.0, 0.4, 3.0], [ 1.0, 1.0 ]}; - """ - - m = 2 - p = 3 - d = 3 - num = np.array([ - [ [0.0, 0.0, 1.0], [1.0, 0.0, 0.0] ], - [ [3.0, -1.0, 1.0], [0.0, 1.0, 0.0] ], - [ [0.0, 0.0, 1.0], [0.0, 2.0, 0.0] ] ]) - - numc = np.zeros((max(1, m, p), max(1, m, p), d), dtype=float) - numc[:p,:m,:] = num - denc = np.array( - [ [1.0, 0.4, 3.0], [ 1.0, 1.0, 0.0 ] ]) - indc = np.array( - [ 2, 1 ], dtype=int) - - nref = 3 - Aref = np.array([ [-1, 0, 0], - [ 0, -0.4, -0.3], - [ 0, 10, 0] ]) - Bref = np.array([ [0, -1], - [1, 0], - [0, 0] ]) - Cref = np.array([ [1, 0, 0.1], - [-1, -2.2, -0.8], - [-2, 0, 0.1] ]) - Dref = np.array([ [0, 1], - [3, 0], - [0, 0] ]) - - nr, A, B, C, D = transform.td04ad('C', m, p, indc, denc, numc) - #print('A=\n', A, '\nB=\n', B, '\nC=\n', C, '\nD=\n', D) - np.testing.assert_equal(nref, nr) - # the returned state space representation is not guaranteed to - # be of one form for all architectures, so we transform back - # to tf and check for equality then - _, _, _, _, _, dcoeff, ucoeff = transform.tb04ad( - nr, m, p, A, B, C, D) - _, _, _, _, _, dcoeffref, ucoeffref = transform.tb04ad( - nref, m, p, Aref, Bref, Cref, Dref) - np.testing.assert_array_almost_equal(dcoeff,dcoeffref) - np.testing.assert_array_almost_equal(ucoeff,ucoeffref) - -def test_td04ad_r(): - """td04ad: Convert with 'R' option - - example program from - http://slicot.org/objects/software/shared/doc/TD04AD.html - """ - - m = 2 - p = 2 - rowcol = 'R' - index = [3, 3] - dcoeff = np.array([ [1.0, 6.0, 11.0, 6.0], [1.0, 6.0, 11.0, 6.0] ]) - - ucoeff = np.array([ [[1.0, 6.0, 12.0, 7.0], [0.0, 1.0, 4.0, 3.0]], - [[0.0, 0.0, 1.0, 1.0], [1.0, 8.0, 20.0, 15.0]] ]) - - nref = 3 - - Aref = np.array([ [ 0.5000, -0.8028, 0.9387], - [ 4.4047, -2.3380, 2.5076], - [-5.5541, 1.6872, -4.1620] ]) - Bref = np.array([ [-0.2000, -1.2500], - [ 0.0000, -0.6097], - [ 0.0000, 2.2217] ]) - Cref = np.array([ [0.0000, -0.8679, 0.2119], - [0.0000, 0.0000, 0.9002] ]) - Dref = np.array([ [1.0000, 0.0000], - [0.0000, 1.0000] ]) - - nr, A, B, C, D = transform.td04ad(rowcol, m, p, index, dcoeff, ucoeff) - #print('A=\n', A, '\nB=\n', B, '\nC=\n', C, '\nD=\n', D) - np.testing.assert_equal(nref, nr) - # order of states is not guaranteed, so we reorder the reference - rindex = np.flip(np.argsort(np.diag(A))) - Arref = Aref[rindex, :][:, rindex] - Brref = Bref[rindex, :] - Crref = Cref[:, rindex] - Drref = Dref - np.testing.assert_array_almost_equal(A, Arref,decimal=4) - np.testing.assert_array_almost_equal(B, Brref,decimal=4) - np.testing.assert_array_almost_equal(C, Crref,decimal=4) - np.testing.assert_array_almost_equal(D, Drref,decimal=4) - - -def test_staticgain(): - """td04ad: Convert a transferfunction to SS with only static gain""" - - # 2 inputs, 3 outputs? columns share a denominator - num = np.array([ [ [1.0], [2.0] ], - [ [0.2], [4.3] ], - [ [1.2], [3.2] ] ]) - p, m, d = num.shape - numc = np.zeros((max(1, m, p), max(1, m, p), d), dtype=float) - numc[:p,:m,:] = num - - # denc, columns share a common denominator - denc = np.array([ [ 1.0], [0.5] ]) - Dc = (num / denc).reshape((3,2)) - idxc = np.zeros((2,), dtype=int) - - # denr, rows share a common denominator - denr = np.array([ [1.0], [0.5], [3.0] ]) - idxr = np.zeros((3,), dtype=int) - Dr = (num / denr[:, np.newaxis]).reshape((3,2)) - - # fails with: - # On entry to TB01XD parameter number 5 had an illegal value - - n, A, B, C, D = transform.td04ad('C', 2, 3, idxc, denc, numc) - #print('A=\n', A, '\nB=\n', B, '\nC=\n', C, '\nD=\n', D) - assert A.shape == (0,0) - assert B.shape == (0,2) - assert C.shape == (3,0) - np.testing.assert_array_almost_equal(D, Dc) - - n, A, B, C, D = transform.td04ad('R', 2, 3, idxr, denr, num) - #print('A=\n', A, '\nB=\n', B, '\nC=\n', C, '\nD=\n', D) - assert A.shape == (0,0) - assert B.shape == (0,2) - assert C.shape == (3,0) - np.testing.assert_array_almost_equal(D, Dr) - -def test_td04ad_static(): - """Regression: td04ad (TFM -> SS transformation) for static TFM""" - from itertools import product - for nout, nin, rc in product(range(1, 6), range(1, 6), ['R', 'C']): - Dref = np.zeros((nout, nin)) - if rc == 'R': - num = np.reshape(np.arange(nout * nin), (nout, nin, 1)) - den = np.reshape(np.arange(1, 1 + nout), (nout, 1)) - index = np.repeat(0, nout) - Dref = num[:nout, :nin, 0] / np.broadcast_to(den, (nout, nin)) - else: - maxn = max(nout, nin) - num = np.zeros((maxn, maxn, 1)) - num[:nout, :nin, 0] = np.reshape( - np.arange(nout * nin), (nout, nin)) - den = np.reshape(np.arange(1, 1 + nin), (nin, 1)) - index = np.repeat(0, nin) - Dref = num[:nout, :nin, 0] / np.broadcast_to(den.T, (nout, nin)) - nr, A, B, C, D = transform.td04ad(rc, nin, nout, index, den, num) - np.testing.assert_equal(nr, 0) - for M in [A, B, C]: - np.testing.assert_equal(M, np.zeros_like(M)) - np.testing.assert_almost_equal(D, Dref) - -def test_mixfeedthrough(): - """Test case popping up from control testing - - a mix of feedthrough and dynamics. The problem from the control - package was somewhere else - """ - num = np.array([ [ [ 0.0, 0.0 ], [ 0.0, -0.2 ] ], - [ [ -0.1, 0.0 ], [ 0.0, 0.0 ] ] ]) - p, m, d = num.shape - numc = np.zeros((max(1, m, p), max(1, m, p), d), dtype=float) - numc[:p,:m,:] = num - denc = np.array([[1.0, 1.1], - [1.0, 0.0]]) - idxc = np.array([1, 0]) - n, A, B, C, D = transform.td04ad('C', 2, 2, idxc, denc, numc) - np.testing.assert_array_almost_equal(D, np.array([[0, 0],[-0.1, 0]])) - -def test_toandfrom(): - A = np.array([[-3.0]]) - B = np.array([[0.1, 0.0]]) - C = np.array([[1.0], - [0.0]]) - D = np.array([[0.0, 0.0], - [0.0, 1.0]]) - - tfout = transform.tb04ad(1, 2, 2, A, B, C, D) - - num = tfout[6] - den = tfout[5] - idxc = np.array([1, 0]) - n, At, Bt, Ct, Dt = transform.td04ad('R', 2, 2, idxc, den, num) - np.testing.assert_array_almost_equal(D, Dt) - np.testing.assert_array_almost_equal(A, At) - -def test_tfm2ss_6(): - """Python version of Fortran test program from - -- Bug in TD04AD when ROWCOL='C' #6 - This bug was fixed in PR #27""" - m = 1 - p = 1 - index = np.array([0]) - dcoeff = np.array([[0.5]]) - ucoeff = np.array([[[32]]]) - n, A, B, C, D = transform.td04ad('R', m, p, index, dcoeff, ucoeff) - assert n == 0 - np.testing.assert_array_almost_equal(D, np.array([[64]])) - n, A, B, C, D = transform.td04ad('C', m, p, index, dcoeff, ucoeff) - assert n == 0 - np.testing.assert_array_almost_equal(D, np.array([[64]])) +from slycot import transform as tf +from slycot.exceptions import SlycotArithmeticError, SlycotParameterError + +from .test_exceptions import assert_docstring_parse + + +class Test_tb04ad: + + @pytest.mark.parametrize( + 'fun, exception_class, erange, checkvars', + ((tf.td04ad, SlycotArithmeticError, (3,), {}),)) + def test_td04ad_docparse(self, fun, exception_class, erange, checkvars): + assert_docstring_parse(fun.__doc__, exception_class, erange, checkvars) + + + def test_td04ad_c(self): + """td04ad: Convert with 'C' option""" + + # for octave: + """ + num = { [0.0, 0.0, 1.0 ], [ 1.0, 0.0 ]; + [3.0, -1.0, 1.0 ], [ 0.0, 1.0 ]; + [0.0, 0.0, 1.0], [ 0.0, 2.0 ] }; + den = { [1.0, 0.4, 3.0], [ 1.0, 1.0 ]; + [1.0, 0.4, 3.0], [ 1.0, 1.0 ]; + [1.0, 0.4, 3.0], [ 1.0, 1.0 ]}; + """ + + m = 2 + p = 3 + d = 3 + num = np.array([ + [ [0.0, 0.0, 1.0], [1.0, 0.0, 0.0] ], + [ [3.0, -1.0, 1.0], [0.0, 1.0, 0.0] ], + [ [0.0, 0.0, 1.0], [0.0, 2.0, 0.0] ] ]) + + numc = np.zeros((max(1, m, p), max(1, m, p), d), dtype=float) + numc[:p,:m,:] = num + denc = np.array( + [ [1.0, 0.4, 3.0], [ 1.0, 1.0, 0.0 ] ]) + indc = np.array( + [ 2, 1 ], dtype=int) + + nref = 3 + Aref = np.array([ [-1, 0, 0], + [ 0, -0.4, -0.3], + [ 0, 10, 0] ]) + Bref = np.array([ [0, -1], + [1, 0], + [0, 0] ]) + Cref = np.array([ [1, 0, 0.1], + [-1, -2.2, -0.8], + [-2, 0, 0.1] ]) + Dref = np.array([ [0, 1], + [3, 0], + [0, 0] ]) + + nr, A, B, C, D = transform.td04ad('C', m, p, indc, denc, numc) + #print('A=\n', A, '\nB=\n', B, '\nC=\n', C, '\nD=\n', D) + np.testing.assert_equal(nref, nr) + # the returned state space representation is not guaranteed to + # be of one form for all architectures, so we transform back + # to tf and check for equality then + _, _, _, _, _, dcoeff, ucoeff = transform.tb04ad( + nr, m, p, A, B, C, D) + _, _, _, _, _, dcoeffref, ucoeffref = transform.tb04ad( + nref, m, p, Aref, Bref, Cref, Dref) + np.testing.assert_array_almost_equal(dcoeff,dcoeffref) + np.testing.assert_array_almost_equal(ucoeff,ucoeffref) + + def test_td04ad_r(self): + """td04ad: Convert with 'R' option + + example program from + http://slicot.org/objects/software/shared/doc/TD04AD.html + """ + + m = 2 + p = 2 + rowcol = 'R' + index = [3, 3] + dcoeff = np.array([ [1.0, 6.0, 11.0, 6.0], [1.0, 6.0, 11.0, 6.0] ]) + + ucoeff = np.array([ [[1.0, 6.0, 12.0, 7.0], [0.0, 1.0, 4.0, 3.0]], + [[0.0, 0.0, 1.0, 1.0], [1.0, 8.0, 20.0, 15.0]] ]) + + nref = 3 + + Aref = np.array([ [ 0.5000, -0.8028, 0.9387], + [ 4.4047, -2.3380, 2.5076], + [-5.5541, 1.6872, -4.1620] ]) + Bref = np.array([ [-0.2000, -1.2500], + [ 0.0000, -0.6097], + [ 0.0000, 2.2217] ]) + Cref = np.array([ [0.0000, -0.8679, 0.2119], + [0.0000, 0.0000, 0.9002] ]) + Dref = np.array([ [1.0000, 0.0000], + [0.0000, 1.0000] ]) + + nr, A, B, C, D = transform.td04ad(rowcol, m, p, index, dcoeff, ucoeff) + #print('A=\n', A, '\nB=\n', B, '\nC=\n', C, '\nD=\n', D) + np.testing.assert_equal(nref, nr) + # order of states is not guaranteed, so we reorder the reference + rindex = np.flip(np.argsort(np.diag(A))) + Arref = Aref[rindex, :][:, rindex] + Brref = Bref[rindex, :] + Crref = Cref[:, rindex] + Drref = Dref + np.testing.assert_array_almost_equal(A, Arref,decimal=4) + np.testing.assert_array_almost_equal(B, Brref,decimal=4) + np.testing.assert_array_almost_equal(C, Crref,decimal=4) + np.testing.assert_array_almost_equal(D, Drref,decimal=4) + + + def test_staticgain(self): + """td04ad: Convert a transferfunction to SS with only static gain""" + + # 2 inputs, 3 outputs? columns share a denominator + num = np.array([ [ [1.0], [2.0] ], + [ [0.2], [4.3] ], + [ [1.2], [3.2] ] ]) + p, m, d = num.shape + numc = np.zeros((max(1, m, p), max(1, m, p), d), dtype=float) + numc[:p,:m,:] = num + + # denc, columns share a common denominator + denc = np.array([ [ 1.0], [0.5] ]) + Dc = (num / denc).reshape((3,2)) + idxc = np.zeros((2,), dtype=int) + + # denr, rows share a common denominator + denr = np.array([ [1.0], [0.5], [3.0] ]) + idxr = np.zeros((3,), dtype=int) + Dr = (num / denr[:, np.newaxis]).reshape((3,2)) + + # fails with: + # On entry to TB01XD parameter number 5 had an illegal value + + n, A, B, C, D = transform.td04ad('C', 2, 3, idxc, denc, numc) + #print('A=\n', A, '\nB=\n', B, '\nC=\n', C, '\nD=\n', D) + assert A.shape == (0,0) + assert B.shape == (0,2) + assert C.shape == (3,0) + np.testing.assert_array_almost_equal(D, Dc) + + n, A, B, C, D = transform.td04ad('R', 2, 3, idxr, denr, num) + #print('A=\n', A, '\nB=\n', B, '\nC=\n', C, '\nD=\n', D) + assert A.shape == (0,0) + assert B.shape == (0,2) + assert C.shape == (3,0) + np.testing.assert_array_almost_equal(D, Dr) + + def test_td04ad_static(self): + """Regression: td04ad (TFM -> SS transformation) for static TFM""" + from itertools import product + for nout, nin, rc in product(range(1, 6), range(1, 6), ['R', 'C']): + Dref = np.zeros((nout, nin)) + if rc == 'R': + num = np.reshape(np.arange(nout * nin), (nout, nin, 1)) + den = np.reshape(np.arange(1, 1 + nout), (nout, 1)) + index = np.repeat(0, nout) + Dref = num[:nout, :nin, 0] / np.broadcast_to(den, (nout, nin)) + else: + maxn = max(nout, nin) + num = np.zeros((maxn, maxn, 1)) + num[:nout, :nin, 0] = np.reshape( + np.arange(nout * nin), (nout, nin)) + den = np.reshape(np.arange(1, 1 + nin), (nin, 1)) + index = np.repeat(0, nin) + Dref = num[:nout, :nin, 0] / np.broadcast_to(den.T, (nout, nin)) + nr, A, B, C, D = transform.td04ad(rc, nin, nout, index, den, num) + np.testing.assert_equal(nr, 0) + for M in [A, B, C]: + np.testing.assert_equal(M, np.zeros_like(M)) + np.testing.assert_almost_equal(D, Dref) + + def test_mixfeedthrough(self): + """Test case popping up from control testing + + a mix of feedthrough and dynamics. The problem from the control + package was somewhere else + """ + num = np.array([ [ [ 0.0, 0.0 ], [ 0.0, -0.2 ] ], + [ [ -0.1, 0.0 ], [ 0.0, 0.0 ] ] ]) + p, m, d = num.shape + numc = np.zeros((max(1, m, p), max(1, m, p), d), dtype=float) + numc[:p,:m,:] = num + denc = np.array([[1.0, 1.1], + [1.0, 0.0]]) + idxc = np.array([1, 0]) + n, A, B, C, D = transform.td04ad('C', 2, 2, idxc, denc, numc) + np.testing.assert_array_almost_equal(D, np.array([[0, 0],[-0.1, 0]])) + + def test_toandfrom(self): + A = np.array([[-3.0]]) + B = np.array([[0.1, 0.0]]) + C = np.array([[1.0], + [0.0]]) + D = np.array([[0.0, 0.0], + [0.0, 1.0]]) + + tfout = transform.tb04ad(1, 2, 2, A, B, C, D) + + num = tfout[6] + den = tfout[5] + idxc = np.array([1, 0]) + n, At, Bt, Ct, Dt = transform.td04ad('R', 2, 2, idxc, den, num) + np.testing.assert_array_almost_equal(D, Dt) + np.testing.assert_array_almost_equal(A, At) + + def test_tfm2ss_6(self): + """Python version of Fortran test program from + -- Bug in TD04AD when ROWCOL='C' #6 + This bug was fixed in PR #27""" + m = 1 + p = 1 + index = np.array([0]) + dcoeff = np.array([[0.5]]) + ucoeff = np.array([[[32]]]) + n, A, B, C, D = transform.td04ad('R', m, p, index, dcoeff, ucoeff) + assert n == 0 + np.testing.assert_array_almost_equal(D, np.array([[64]])) + n, A, B, C, D = transform.td04ad('C', m, p, index, dcoeff, ucoeff) + assert n == 0 + np.testing.assert_array_almost_equal(D, np.array([[64]])) diff --git a/slycot/tests/test_tf01md.py b/slycot/tests/test_tf01md.py new file mode 100644 index 00000000..d349f5f5 --- /dev/null +++ b/slycot/tests/test_tf01md.py @@ -0,0 +1,15 @@ +import warnings + +import pytest + +from slycot import transform as tf +from slycot.exceptions import SlycotArithmeticError + +from .test_exceptions import assert_docstring_parse + + +class Test_tf01md: + + @pytest.mark.skip(reason="not implemented") + def test_tf01md(self): + pass \ No newline at end of file diff --git a/slycot/tests/test_tf01rd.py b/slycot/tests/test_tf01rd.py new file mode 100644 index 00000000..f67226eb --- /dev/null +++ b/slycot/tests/test_tf01rd.py @@ -0,0 +1,15 @@ +import warnings + +import pytest + +from slycot import transform as tf +from slycot.exceptions import SlycotArithmeticError + +from .test_exceptions import assert_docstring_parse + + +class Test_tf01rd: + + @pytest.mark.skip(reason="not implemented") + def test_tf01rd(self): + pass \ No newline at end of file diff --git a/slycot/tests/test_tg01ad.py b/slycot/tests/test_tg01ad.py index 864d41be..64f4b84c 100644 --- a/slycot/tests/test_tg01ad.py +++ b/slycot/tests/test_tg01ad.py @@ -1,6 +1,3 @@ -# =================================================== -# tg01ad tests - import numpy as np from numpy.testing import assert_almost_equal, assert_equal, assert_raises diff --git a/slycot/tests/test_tg01fd.py b/slycot/tests/test_tg01fd.py index 902d1ef2..1168cc80 100644 --- a/slycot/tests/test_tg01fd.py +++ b/slycot/tests/test_tg01fd.py @@ -1,6 +1,3 @@ -# =================================================== -# tg01fd tests - import numpy as np from numpy.testing import assert_almost_equal, assert_equal, assert_raises diff --git a/slycot/tests/test_transform.py b/slycot/tests/test_transform.py deleted file mode 100644 index 928ac5d4..00000000 --- a/slycot/tests/test_transform.py +++ /dev/null @@ -1,22 +0,0 @@ -# -# test_transform.py - generic tests for transform routines -# repagh