SLICOT SB10YD wrapper

[JinmingChe]

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.

[bnavigator]

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.

1. Create a F2PY routine signature in synthesis.pyf
2. 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 Warns and Raises blocks
3. 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 to test_sb.py (or create a new separate test file)

... we should try to generalize these instructions and put them into the developer wiki.

[JinmingChe]

Thanks very much for your reply and guides. I will follow your advises, learn to pull request and add the SB10YD.

[JinmingChe]

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 ?

[bnavigator]

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 . or python setup.py install in a venv. You need all build dependencies (fortran compiler, CMake, scikit-build, LAPACK/BLAS).

Even without a venv, python setup.py build will compile and create a full package in _skbuild/<your-platform>/cmake-install. You can put this into your PYTHONPATH and run

import slycot._wrapper
print(slycot._wrapper.sb10yd.__doc__)

to see the automatically generated call signature of the sb10yd function.

[wenboklose]

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 array

My 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 sb10yd

My 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.
Bo

[bnavigator]

double 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 A to 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!

[wenboklose]

Hi, Thanks so much. I am still learning Git, so what is a pull request?

[wenboklose]

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 value

I am not confidence on the synthesis.py function, can you help me to check it out?

Bo

[bnavigator]

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-control with python-control/Slycot and all remaining standalone python-control with Slycot)

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.

[bnavigator]

info < 1 means that _wrapper.sb10yd() did not pass the correct arguments to the Fortran routine. If you remove the :info < 0: line from the docstring, the raise_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().

[bnavigator]

Closed by #203