Repository navigation
SLICOT SB10YD wrapper #153
Description
Activity
Hello @JinmingChe,
You can find the documentation for SB10YD at http://slicot.org/objects/software/shared/doc/SB10YD.html and the source code along with the documentation header is at https://github.com/SLICOT/SLICOT-Reference/blob/main/src/SB10YD.f
Please feel free to create a pull request and we can guide you along with what needs to be improved in order to merge it into Slycot. You can use #118 as an example.
- Create a F2PY routine signature in
synthesis.pyf - Create a wrapper function in synthesis.py. Please convert the Fortran doc header into an equivalent Python docstring adhering to the numpydoc standard along with our special treatment of
WarnsandRaisesblocks - There is no Fortran example for
SB10YD. Please try to find or create one or more test for the new python function and add them totest_sb.py(or create a new separate test file)
... we should try to generalize these instructions and put them into the developer wiki.
- Create a F2PY routine signature in
Thanks very much for your reply and guides. I will follow your advises, learn to pull request and add the SB10YD.
Hi, I follow the instruction in F2PY, use the smart way create a signature file from SB10YD.f. And I modified the SB10YD.pyf. Next, should I copy the subroutine sb10yd to synthesis.pyf and then build the extension module by running synthesis.pyf ?
Next, should I copy the subroutine sb10yd to synthesis.pyf
Yes
and then build the extension module by running synthesis.pyf ?
No, build the full wrapper in a development environment. If you are using conda, there is a wikipage.
Alternatively use
pip install .orpython setup.py installin a venv. You need all build dependencies (fortran compiler, CMake, scikit-build, LAPACK/BLAS).Even without a venv,
python setup.py buildwill compile and create a full package in_skbuild/<your-platform>/cmake-install. You can put this into your PYTHONPATH and runimport slycot._wrapper print(slycot._wrapper.sb10yd.__doc__)
to see the automatically generated call signature of the sb10yd function.
Hi, I followed the approach mentioned above, got the following error:
ValueError: failed to create intent(cache|hide)|optional array-- must have defined dimensions but got (-1,1,) The above exception was the direct cause of the following exception: Traceback (most recent call last): File "/home/pi/Coding/Stability/TF_Estimate.py", line 56, in <module> A, B, C, D = sb10yd(0,0,9,real_H1_resp,imag_H1_resp,frequency,1,-1) File "/home/pi/Coding/Stability/venv/lib/python3.7/site-packages/slycot-0.4.0-py3.7-linux-armv7l.egg/slycot/synthesis.py", line 1766, in sb10yd A,B,C,D,info = _wrapper.sb10yd(discfl,flag,lendat,rfrdat,ifrdat,omega,n,tol,ldwork,lzwork) ValueError: failed in converting hidden `a' of _wrapper.sb10yd to C/Fortran arrayMy example is:
omega = 2*np.pi frequency = omega*np.arange(1,10,1) H1 = control.frd(H,frequency) real_H1_resp = [4.74192334e-03, 1.18688787e-03, 5.27621687e-04, 2.96810037e-04, 1.89965190e-04, 1.31922823e-04, 9.69240212e-05, 7.42080156e-05, 5.86337982e-05]#H1.fresp.real imag_H1_resp = [-0.11917753, -0.05965949, -0.03978174, -0.0298386, -0.02387173, -0.01989349, -0.01705176, -0.01492041, -0.01326265]#H1.fresp.imag print(H1.fresp) print(real_H1_resp) print(imag_H1_resp) print(frequency) A, B, C, D = sb10yd(0,0,9,real_H1_resp,imag_H1_resp,frequency,1,-1)My synthesis.pyf:
subroutine sb10yd(discfl,flag,lendat,rfrdat,ifrdat,omega,n,a,lda,b,c,d,tol,iwork,dwork,ldwork,zwork,lzwork,info) ! in SB10YD.f integer intent(in),check(discfl>=0 && discfl<=1) :: discfl integer intent(in),check(discfl>=0 && discfl<=1) :: flag integer intent(in),check(lendat>=2) :: lendat double precision intent(in),dimension(lendat),depend(lendat) :: rfrdat double precision intent(in),dimension(lendat),depend(lendat) :: ifrdat double precision intent(in),dimension(lendat),depend(lendat) :: omega integer intent(inout),check(n>0 && n<=lendat-1) :: n double precision intent(out),dimension(lda,n),depend(n) :: a integer intent(hide),check(shape(a,0)==lda),depend(a) :: lda=shape(a,0) double precision intent(out),dimension(n) :: b double precision intent(out),dimension(n) :: c double precision intent(out),dimension(1) :: d double precision :: tol integer intent(hide,cache),dimension(max(2,2*n)),depend(n) :: iwork double precision intent(hide,cache),dimension(ldwork) :: dwork integer optional,check(ldwork>=max(3,6*n)),depend(n) :: ldwork complex*16 intent(hide,cache),dimension(lzwork) :: zwork integer optional,:: lzwork integer intent(out):: info end subroutine sb10ydMy synthesis.py
def sb10yd(discfl,flag,lendat,rfrdat,ifrdat,omega,n,tol,ldwork=None,lzwork=None): """ A,B,C,D = sb10yd(discfl,flag,lendat,rfrdat,ifrdat,omega,n,tol,[ldwork,lzwork]) To fit frequency response data with a stable, minimum phase SISO system Parameters ---------- discfl : int Indicatres the type of the system, as follows: = 0: continuous-time system; = 1: discrete-time system. flag : int If flag = 0, then the system zeros and poles are not constrained. If flag = 1, then the system zeros and poles will have negative real parts in the continuous-time case, or moduli less than 1 in the discrete-time case. Consequently, flag must be equal to 1 in mu-synthesis routines. lendat : int The length of the vectors rfrdat, ifrdat and omega. length >= 2. rfrdat : double precision array, dimension (lendat) The real part of the frequency data to be fitted. ifrdat : double precision array, dimension (lendat) The imaginary part of the frequency data to be fitted. omega : double precision array, dimension (lendat) The frequencies corresponding to rfrdat and ifrdat. These value must be nonnegative and monotonically increasing. Additionally, for discrete-time systems they must be between 0 and PI. n : integer On entry, the desired order of the system to be fitted. n <= lendat-1. tol : int, optional The length of the cache array. ldwork >= max( 1, 2*n*n + 2*n + n*MAX( 5, n + m + np ) ). For good performance, ldwork must generally be larger. Returns ------- A : (n, n) double precision array The leading n-by-n part of this array contains the matrix A. B : (n) double precision array The computed vector B. C : (n) double precision array The computed vector C. D : (1) double precision array The computed scalar D. Raises ------ SlycotArithmeticError :info == 0: successful exit; :info < 0: if infor = -i, the i-th argument had an illegal value :info = 1: if the discret --> continous transformation cannot be made; :info = 2: if the system poles cannot be found; :info = 3: if the inverse system cannot be found, i.e., D is (close to) zero; :info = 4: if the system zeros cannot be found; :info = 5: if the state-space representation of the new transfer function T(s) cannot found; :info = 6: if the continous --> discrete transformation cannot be made. The iteration for computing singular value decomposition did not converge. """ hidden = ' (hidden by the wrapper)' arg_list = ['discfl', 'flag', 'lendat', 'rfrdat', 'ifrdat', 'omega', 'n', 'A', 'lda'+hidden, 'B', 'C', 'D', 'tol', 'iwork'+hidden, 'dwork'+hidden, 'ldwork', 'zwork'+hidden, 'lzwork', 'info'] if ldwork is None: lw1 = 2*lendat + 4*2048 lw2 = lendat + 6*2048 mn = min(2*lendat,2*n+1) if n > 0: lw3 = 2*lendat*(2*n+1) + max(2*lendat,2*n+1) + max(mn+6*n+4,2*mn+1) elif n == 0: lw3 = 4*lendat + 5 if flag == 1: lw4 = max(n*n+5*n,6*n+1+min(1,n)) elif flag == 0: lw4 = 0 ldwork = max(2, lw1, lw2, lw3, lw4) A,B,C,D,info = _wrapper.sb10yd(discfl,flag,lendat,rfrdat,ifrdat,omega,n,tol,ldwork,lzwork) raise_if_slycot_error(info, arg_list, sb10yd.__doc__) return A[:nsys,:nsys],B[:nsys],C[:nsys],D[:1]Hope you can help, thanks.
Bodouble precision intent(out),dimension(lda,n),depend(n) :: a integer intent(hide),check(shape(a,0)==lda),depend(a) :: lda=shape(a,0)That is kind of circular.
You can either do it like
ab08nd():double precision intent(in),dimension(lda,*),check(shape(a,1)>=n) :: a integer intent(hide),check(lda>=max(1,n)) :: lda=shape(a,0)Or restrict the shape of
Ato be exactly square, like all the other methods in synthesis.pyf do. No need to supply extra unused rows from Python, although the FORTRAN routines would allow it for memory management reasons.integer intent(in),check(n>=0),required :: n double precision dimension(n,n),depend(n) :: a integer intent(hide),depend(a) :: lda=shape(a,0)Please create a pull request!
Hi, Thanks so much. I am still learning Git, so what is a pull request?
HI, I got new error messages:
Traceback (most recent call last): File "/home/pi/Coding/Stability/TF_Estimate.py", line 56, in <module> A, B, C, D = sb10yd(0,0,9,real_H1_resp,imag_H1_resp,frequency,2,-1) File "/home/pi/Coding/Stability/venv/lib/python3.7/site-packages/slycot-0.4.0-py3.7-linux-armv7l.egg/slycot/synthesis.py", line 1768, in sb10yd raise_if_slycot_error(info, arg_list, sb10yd.__doc__) File "/home/pi/Coding/Stability/venv/lib/python3.7/site-packages/slycot-0.4.0-py3.7-linux-armv7l.egg/slycot/exceptions.py", line 236, in raise_if_slycot_error raise globals()[exception](fmessage, info) slycot.exceptions.SlycotArithmeticError: if infor = -i, the i-th argument had an illegal valueI am not confidence on the synthesis.py function, can you help me to check it out?
Bo
It is easier to review code in a pull request. Please follow:
https://github.com/python-control/python-control/wiki/How-to-contribute-with-a-pull-request (replace all
python-control/python-controlwithpython-control/Slycotand all remaining standalonepython-controlwithSlycot)Github help: https://docs.github.com/en/pull-requests/collaborating-with-pull-requests
This will also run the tests on our servers as soon as you commit code and we see the occuring errors immediately.
info < 1means that_wrapper.sb10yd()did not pass the correct arguments to the Fortran routine. If you remove the:info < 0:line from the docstring, theraise_if_slycot_error()function will tell you which argument was wrong. Double check the number and order of input and output parameters of_wrapper.sb10yd().Closed by #203
Hi, I am a student working on speech signal processing.
I want to implement the fitfrd function in python. It seems that sb10yd of slycot has implemented this function. But it is not currently wrapped. It would be so kind if you can give me some advices.
Sincerely.