TB05AD wrapper #10
TB05AD wrapper #10
Conversation
|
The email address used in the commits begins with |
|
No, it was not intentional. There was an error in my .gitconfig file. According to this, I can change the commit email by rebasing. Is that the preferred method? |
|
See @slivingston's take on rebasing in pull requests here. In brief, rebasing in PRs is fine. You can't depend on python-control in Slycot or its tests; that's why Travis CI is failing. I'd like to do a more thorough review, but it will take some time. |
rabraker
force-pushed
the
master
branch
2 times, most recently
from
8f505a2 to
60d6c2c
Compare
|
Thanks for the link. I didn't know that a tool like travis CL even existed until today. After removing the dependence on python-control and making a couple tweaks to the transform.pyf file, the builds all succeed now. |
There was a problem hiding this comment.
Looks good, thank you. I've made a few suggestions, but nothing major.
I haven't looked at your proposed fix to TB05AD.f -- I'll need to study DGEBAL, and try to resurrect what little Fortran knowledge I have to do that.
I think @slivingston would prefer you open a new issue for this SLICOT bug, and remove the "Proposed fix" c.omment in this PR.
slycot/examples.py
Outdated
| @@ -240,3 +240,26 @@ def tb01pd_example(): | |||
| print('reduced order', out[-2]) | |||
| print(out) | |||
|
|
|||
|
|
|||
| def tb05ad_example(): | |||
| import numpy as np | |||
There was a problem hiding this comment.
A docstring wouldn't hurt here ("Demonstrate TB05AD with a second-order system with natural frequency 10r/s and damping 1.05").
slycot/src/transform.pyf
Outdated
| @@ -135,7 +135,97 @@ subroutine tb04ad_c(rowcol,n,m,p,a,lda,b,ldb,c,ldc,d,ldd,nr,index_bn,dcoeff,lddc | |||
| double precision intent(hide,cache),dimension(ldwork) :: dwork | |||
| integer optional :: ldwork = max(1,n*(n+1)+max(n*p+2*n+max(n,p),max(3*p,m))) | |||
| integer intent(out) :: info | |||
| end subroutine tb04ad_c | |||
| end subroutine tb04ad_c | |||
There was a problem hiding this comment.
assume you didn't mean to indent end subroutine here?
slycot/src/transform.pyf
Outdated
| end subroutine tb04ad_c | ||
| end subroutine tb04ad_c | ||
| subroutine tb05ad_ag(baleig,inita,n,m,p,freq,a,lda,b,ldb,c,ldc,rcond,g,ldg,evre,evim,hinvb,ldhinv,iwork,dwork,ldwork,zwork,lzwork,info) | ||
| fortranname tb05ad |
There was a problem hiding this comment.
A one-line comment (like you have for tb05ad_nh) would be good here.
From the TB05AD docs, 'C' and 'B' (and also 'E', which is the same as 'B'?) are other first letter options. Do you know if these other options are useful?
There was a problem hiding this comment.
After reading the documentation again I do think 'C' could be useful. The difference between 'N' and 'C' is that 'C' computes RCOND (an estimate of the reciprocal of the condition number of H), 'N' does not. Here,
H = (freq*I - A)
so if A has modes with zero or very little damping H will become poorly conditioned for frequencies near those modes. So 'C' could be useful in order to throw a warning that the calculation at a particular frequency may be inaccurate if the condition number is really bad. You could then either just throw that point away or use the information to adjust the points you're evaluating at.
I'd be willing to add option 'C' at some point in the future if there is interest in that.
The difference between 'A' and 'B' is that 'A' calculates eigenvalues, and RCOND while 'B' only computes eigenvalues. Both balance matrix A. It is my understanding that these two options are only to be used on the "initial" call to TB05AD. In that sense, I don't think its worth adding a different wrapper just to avoid estimating the condition number in a single call.
| integer optional,depend(n) :: lzwork = n*n+2*n | ||
| integer intent(out) :: info | ||
| end subroutine tb05ad_ag | ||
| subroutine tb05ad_ng(baleig,inita,n,m,p,freq,a,lda,b,ldb,c,ldc,rcond,g,ldg,evre,evim,hinvb,ldhinv,iwork,dwork,ldwork,zwork,lzwork,info) |
There was a problem hiding this comment.
Do we need NG? Aren't AG and NH sufficient?
[Later] oh, I see you have a comment about this in transform.py. Fair enough.
slycot/tests/test.py
Outdated
|
|
||
| class test_tb05ad(unittest.TestCase): | ||
|
|
||
| def setUp(self): |
There was a problem hiding this comment.
Omit setUp if it does nothing.
slycot/transform.py
Outdated
| used to transform A1 to Hessenberg form. See docs for lappack | ||
| DGEHRD(): | ||
| http://www.netlib.org/lapack/explore-3.1.1-html/dgehrd.f.html | ||
| However, it does not apparently to contain information on P, the |
There was a problem hiding this comment.
"apparently to contain" =>"apparently contain" ?
slycot/transform.py
Outdated
| e.info = out[-1] | ||
| raise e | ||
| if out[-1] == 1: | ||
| error_text = ("More than 30*"+str(n)+"iterations are required " |
There was a problem hiding this comment.
why 30*str(n) ? The docs just say 30 iterations.
| # Checks done, do computation. | ||
| if job == 'AG': | ||
| out = _wrapper.tb05ad_ag(n, m, p, jomega, A, B, C) | ||
| error_handler(out, arg_list) |
There was a problem hiding this comment.
This is good; it's tidier than the error handling elsewhere in Slycot.
It might be good to turn your error_handler into the constructor for (currently non-existant) SlycotValueError, inherited from ValueError, that handle alls the re-indexing, etc., in one place. However, that's for another PR.
slycot/transform.py
Outdated
| @@ -333,6 +333,254 @@ def tb04ad(n,m,p,A,B,C,D,tol1=0.0,tol2=0.0,ldwork=None): | |||
| kdcoef = max(index)+1 | |||
| return A[:Nr,:Nr],B[:Nr,:m],C[:p,:Nr],Nr,index,dcoeff[:porm,:kdcoef],ucoeff[:porm,:porp,:kdcoef] | |||
|
|
|||
|
|
|||
| def tb05ad(n, m, p, jomega, A, B, C, job='NG'): | |||
| """At, Bt, Ct, g_jw, rcond, ev, hinvb = tb05ad_a(n, m, p, jomega, A, B, C): | |||
There was a problem hiding this comment.
This doc string isn't right: what is returned depends on job (which isn't listed as an argument).
I'm not sure if we (Slycot and python-control) have a good convention for different return forms. Numpy appears to just omit the return list in the header (e.g., numpy.linalg.svd) and give details below.
tb05ad_a => tb05ad.
slycot/transform.py
Outdated
| e = ValueError("The shape of A is (" + str(A.shape[0]) + "," + | ||
| str(A.shape[1]) + "), but expected (" + str(n) + | ||
| "," + str(n) + ")") | ||
| e.info = -6 |
There was a problem hiding this comment.
I'm not sure it's worth faking the info attribute. You haven't done it for the "Unrecognized job" error at the end (where it's not clear what it would be...).
|
Thanks for the review. I have quite a few deadlines coming in the next couple of weeks, so it may be a bit before I can get to it though. |
|
I believe this have addressed all your comments. I made a few more minor changes to the docstring for tb05ad in transform.py that fixed a few other typos and made the wording between 'job' descriptions more consistent. I can rebase all of it into a single commit, but I thought having them separated might be useful for the moment. |
slycot/tests/test.py
Outdated
| @@ -64,6 +64,7 @@ def test_td04ad_static(self): | |||
| nr,a,b,c,d = transform.td04ad(rc,nin,nout,index,den,num) | |||
|
|
|||
|
|
|||
|
|
|||
There was a problem hiding this comment.
a nitpick: please remove this new newline, then this file won't be in the changeset at all.
|
I've removed that newline as requested. I should be able to file a bug report for TB05AD.f in the next few days. |
|
A Bug report has been submitted, as requested. |
|
OK, thanks for all the fixes. This looks good to me. @slivingston, what do you think? |
|
Last thing: could you please squash this to one commit? I have push rights now, will merge to master when that's done. |
…hon function tb05ad in transform.py.
|
It's all squashed now. |
This pull request is to address issue 116 in python-control. It provides a wrapper for the TB05AD slicot routine which computes the frequency response of a statespace system. The routine is mostly useful when computing the frequency response at many frequencies.
The wrapper implements three different modalities for TB05AD:
1). You pass in general (A, B, C) matrices. The A matrix is balanced, transformed to upper Hessenberg form, and the B, and C matrices are appropriately transformed. This returns the frequency response as well as the transformed (A, B, C) matrices.
2). The same as (1), except no balancing is performed.
3). You pass matrices which have already been transformed to upper Hessenberg form.
These three different jobs are wrapped as three different f2py functions in transform.pyf. In python, they are all called via the same function in transform.py 'tb05ad' by setting a different job flag. These are 'AG' (All, general), 'NG' (none, general), 'NH' (none, hessernberg) respectively.
Using job 'NG' (2) then job 'NH' (3) works perfectly as far as I can tell.
For certain A matrices, job 'AG' (1) gives the wrong the results. I believe this is a bug in TB05AD.f circa line 354. Above that line, I have commented out my proposed change which I believe fixes the issue. But I wanted somebody else to see it fail as well to make sure I'm not crazy :). When doing the balancing, TB05AD calls the lapack routine DGEBAL. From what I can tell, TB05AD misinterprets the permutation information provided by DGEBAL, but this misinterpretation is only a problem when the permutation corresponds to non-symmetric permutation matrix. This is similar to the scipy matrix_balance bug which I reported here. Before this was fixed, TB05AD and linalg.matrix_balance would yield the same transformed set (A, B, C), but you would get the wrong frequency response.
The failure is documented in the unit tests as a known failure.
I deviated some in writing the unit tests from what is currently implemented. The motivation was to test the three different routines on the same set of matrices while avoiding copy-and-paste coding. I'm sure there is a better way...
Let me know what you think.
Thanks!
Arnold