Closed KybernetikJo closed 1 year ago
bnavigator
commented
1 year ago Hi, thanks for your contribution. I didn't check in detail yet, but if it is really true that there is a call by reference, then this would be the case in many more functions than ab13bd and we have a general problem in the whole package.
KybernetikJo
commented
1 year ago The following changes are actually enough:
integer check(n>=0) :: n
integer check(m>=0) :: m
integer check(p>=0) :: p
double precision intent(in,out,copy), dimension(n,n),depend(n) :: a
integer intent(hide),depend(a) :: lda = shape(a,0)
double precision intent(in,out,copy), dimension(n,m),depend(n,m) :: b
integer intent(hide),depend(b) :: ldb = shape(b,0)
double precision intent(in,out,copy), dimension(p,n),depend(n,p) :: c
integer intent(hide),depend(c) :: ldc = shape(c,0)
double precision intent(in,out,copy), dimension(p,m),depend(m,p) :: d
integer intent(hide),depend(d) :: ldd = shape(d,0)Intent(in,out,copy) copies the parameter. The function has no side effects on the state space matrices A,B,C,D
TODO:
bnavigator
commented
1 year ago Actually, according to https://numpy.org/doc/stable/f2py/signature-file.html, we should have the desired intent(in) by default.
inThe corresponding argument is considered to be input-only. This means that the value of the argument is passed to a Fortran/C function and that the function is expected to not change the value of this argument.The following rules apply:
- If none of intent(in | inout | out | hide) are specified, intent(in) is assumed.
I wonder if it the stray
causes the inputs to be altered?
Could you check if you get altered inputs in other wrappers as well?
bnavigator
commented
1 year ago Intent(in,out,copy) copies the parameter.
Yes, but we actually don't want the output anywhere.
KybernetikJo
commented
1 year ago Thanks for your reply. I appreciate that.
First, I am going to build up some theoretical understanding about f2py. Then I can run some experiments.
Intent(in,out,copy) copies the parameter.
Yes, but we actually don't want the output anywhere.
I am not sure. Why SLICOT returns them?
roryyorke
commented
1 year ago FWIW I can reproduce this with Slycot 0.5.2 on Ubuntu 20.04 in a Conda environment.
From the SLICOT docs for AB13MD, the A, B, C and D arguments are input/output; they contain the state-space matrices of "the numerator factor Q of the right coprime factorization with inner denominator of G". So I think we violate this expectation of f2py:
inThe corresponding argument is considered to be input-only. This means that the value of the argument is passed to a Fortran/C function and that the function is expected to not change the value of this argument.
Specifically, we violate "the function is expected to not change the value of this argument". Whoops! As @bnavigator noted, this could be a problem in quite a few places, we'll have to check.
roryyorke
commented
1 year ago See https://github.com/roryyorke/Slycot/tree/ab13md-alter-args, whence:
https://github.com/python-control/Slycot/commit/0277e72467233816866187cf120e3630c01da9e8 - add regression test from @KybernetikJo 's example; test fails
https://github.com/python-control/Slycot/commit/46993569110cce710ac078421111d76e6070d6de , remove (strange?) intent(out) :: ab13bd, test fails:
diff --git a/slycot/src/analysis.pyf b/slycot/src/analysis.pyf
index 633ff6ec..183a22f2 100644
--- a/slycot/src/analysis.pyf
+++ b/slycot/src/analysis.pyf
@@ -321,7 +321,6 @@ function ab13bd(dico,jobn,n,m,p,a,lda,b,ldb,c,ldc,d,ldd,nq,tol,dwork,ldwork,iwar
integer intent(hide),depend(n,m,p) :: ldwork = max(1,max(m*(n+m)+max(n*(n+5),max(m*(m+2),4*p)),n*(max(n,p)+4)+min(n,p)))
integer intent(out) :: iwarn
integer intent(out) :: info
- double precision intent(out) :: ab13bd
end function ab13bd
subroutine ab13dd(dico,jobe,equil,jobd,n,m,p,fpeak,a,lda,e,lde,b,ldb,c,ldc,d,ldd,gpeak,tol,iwork,dwork,ldwork,cwork,lcwork,info) ! in AB13DD.f
character intent(in) :: dicointent(in,out,copy) as suggested by @KybernetikJo , test passes:diff --git a/slycot/src/analysis.pyf b/slycot/src/analysis.pyf
index 183a22f2..999e3b65 100644
--- a/slycot/src/analysis.pyf
+++ b/slycot/src/analysis.pyf
@@ -307,13 +307,13 @@ function ab13bd(dico,jobn,n,m,p,a,lda,b,ldb,c,ldc,d,ldd,nq,tol,dwork,ldwork,iwar
integer check(n>=0) :: n
integer check(n>=0) :: m
integer check(n>=0) :: p
- double precision dimension(n,n),depend(n) :: a
+ double precision intent(in,out,copy),dimension(n,n),depend(n) :: a
integer intent(hide),depend(a) :: lda = shape(a,0)
- double precision dimension(n,m),depend(n,m) :: b
+ double precision intent(in,out,copy),dimension(n,m),depend(n,m) :: b
integer intent(hide),depend(b) :: ldb = shape(b,0)
- double precision dimension(p,n),depend(n,p) :: c
+ double precision intent(in,out,copy),dimension(p,n),depend(n,p) :: c
integer intent(hide),depend(c) :: ldc = shape(c,0)
- double precision dimension(p,m),depend(m,p) :: d
+ double precision intent(in,out,copy),dimension(p,m),depend(m,p) :: d
integer intent(hide),depend(d) :: ldd = shape(d,0)
integer intent(out) :: nq
double precision :: toland https://github.com/python-control/Slycot/commit/2625e1ffb512565db2aef883e4515e7a8c2e337f, intent(in,copy) also passes.
roryyorke
commented
1 year ago I don't find the f2py docs all that clear. I don't understand the difference between intent variants in,out, inout, and inplace; nor the difference between copy and overwrite.
For the record, intent(in, overwrite) fails the ab13bd test.
roryyorke
commented
1 year ago according to
grep -n "dimension" ~/src/slycot/slycot/src/*.pyf|grep -v "intent" >> audit.out
grep -n "dimension" ~/src/slycot/slycot/src/*.pyf|fgrep "intent(in)" >> audit.outthere are 139 array arguments with implicit or explicit intent(in) on master (including the 4 identified by @KybernetikJo ).
Some of these are fine, e.g., the first arg to AB05MD is marked SLICOT as "input".
We could wholesale change all of these to intent(in,copy), but it's probably best to audit all these.
Worst-case is that some dependency is relying on this bad behaviour, but I think that is unlikely.
ab13bd specifically is not used in python-control, at least.
roryyorke
commented
1 year ago I've compared each case from above with the comments in the associated SLICOT files, and marked each line OK or NOTOK. There are 19 "NOTOK" lines, including the 4 already found.
I'd appreciate spot checks of the results.
Next steps:
roryyorke
commented
1 year ago I think there no problems at python-control level.
Fix options:
intent(in,out) and update Slycot docstrings where needed.intent(in,out,copy), and returned the resulting matricesintent(in, copy): users expecting mutation will be surprised.
| Function | Slycot args altered | Python-control |
|---|---|---|
| ab13bd | Yes | Not used |
| sb03od | Yes | Used, but args local to func |
| sg03od | Yes | Not used |
| tc01od_r | Yes | Not used. |
| ab08nd | No - mistake in analysis | Used, but no problem |
| tb05ad_nh | No - case NG alters in SLICOT, is properly handled | Used, but no problem |
H2 / L2 norm
ABCD passed straight through
docstring does not indicate that arg will be changed
not used in python-control
Conclusions:
lyapunov equation
problem SLICOT args: a, q if FACT is not 'F'. FACT can be specified by caller, is default 'N'.
Advertised as "in-out" in doc string.
used in python-control statefbk.gram. Args A and Q are constructed inside the function.
Conclusions:
cholesky factor of generalized lyapunov equation solution
problem args: a, e, q, z. Passed straight through to SLICOT.
Not used in python-control
Conclusions:
problem args: pcoeff, qcoeff
This is in,out for tc01od_l
pcoeff and qcoeff are both mutated and returned.
Conclusions
mistake in initial analysis; args are input only in SLICOT
used in slycot.tb05ad
problem args a,b,c
Complex freq. resp for A hessenberg (job='NH') A,B,C passed straight through.
Ah - but only output if inita is "NG". For this case tb05ad_ng, with `intent(in,out,copy)` is used.
Conclusions:
KybernetikJo
commented
1 year ago I think the line
double precision intent(out) :: ab13bd https://github.com/python-control/Slycot/commit/46993569110cce710ac078421111d76e6070d6de
is there because of ab13bd is a function, and not a subroutine. Fortran functions do have a return value, subrountines do not.
I think ab13bd is the only function in SLICOT, I do not know why, but there might be a reason. I am not sure deleting the line, is the right way. It might be it though.
KybernetikJo
commented
1 year ago I was wrong. ab13cd is another function. There are a few of them in SLICOT.
It also has a return value.
AB13CD = ( GAMMAL + GAMMAU )/TWOFWIW, I have done some background reading about F2PY, and this is my current understanding. I am a newbie to fortran, therefore you should take this with a bit of toast.
The stable version of f2py stable do have some duplication in the documentation.
I fixed that, so the dev version f2py dev is a better read.
The online "book" python fortran can also be helpful.
The f2py was design with Fortran90/95 in mind.
The goal of F2PY is, to make the use of Fortran in Python easy, and to make the API PYTHONIC. The behavior of F2PY depends heavily on the data types used in the Fortran procedures, the data types used in Python and on the F2PY signature (none, infile.f or .pyf). Therefore, very simple conclusions cannot be drawn, some understanding of Fortran seems to be necessary and good templates might be helpful.
First, Fortran has two procedures, subroutines and functions. (This effects the F2PY signature file pyf) Both subroutines and functions have access to variables in the parent scope by argument association, this is similar to call by reference.
All arguments in Fortran77 can be seen as call by reference, by default.
In Fortran90/95 the intent attribute has been introduced intent(in), intent(out), intent(inout). This allows the compiler to check for unintentional errors and provides self-documentation (in Fortran77 the API was only documented in comments). In addition, some compilers can do code optimization based on intent, mabye?. Is there performance advantage of using intent(in) and intent(out)?
intent in Fortran90/95 Writing fast Fortran routines for Python
intent(in): The variable is an input to the subroutine only. Its value
must not be changed during the course of the subroutine.
intent(out): The variable is an output from the subroutine only. Its
input value is irrelevant. We must assign this variable a
value before exiting the subroutine.
intent(inout): The subroutine both uses and modifies the data in the
variable. Thus, the initial value is sent and we ultimately
make modifications base on it.
In Fortran there is also the concept of a pure function, which do not have side effects. What is a pure function?
In Fortran90/95 a new keyword pure was introduced.
In Fortran2003 a value attribute has been introduced. Now call by value is possible. F2PY do not have a good support of Fortran2003 features, so it can be mostly ignored.
KybernetikJo
commented
1 year ago There are many option how to design F2PY signature files. The intent attribute in signature files as tons of options.
===================
However, currently we need a solution for the case, where the parameter are labeled as input/output in SLICOT,
are arrays in fortran and numpy arrays in python. Mostly for vectors and matrices e.g state space matrices A,B,C,D.
https://github.com/python-control/Slycot/pull/200#issuecomment-1646797617
Fix options:
- API-preserving and easiest: change to
intent(in,out)and update Slycot docstrings where needed.- API-changing 1: change to
intent(in,out,copy), and returned the resulting matricesAPI-changing 2: change to
intent(in, copy): users expecting mutation will be surprised.
- removes capability for users wanting SLICOT results
intent(in, copy): users expecting mutation will be surprised.This option would potentially severely limit SLICOT procedures.
I think, if we have an INPUT/OUTPUT argument in SLICOT, we should use intent(in,out) or intent(in,out,copy).
intent(in,out) and update Slycot docstrings where needed.Here, we would have to deal with Slycot function side effects somewhere.
We could use shallow copies in the outer wrapper python files.
Or we could deal with the side effects in python-control.
This could lead to many bugs.
Having said that, there are inplace function in python
[].sort()I think in numpy there are not that common. For instance
a_out = np.sort(a)intent(in,out,copy), and returned the resulting matricesI still think this is the best option for us. This makes the slycot._wrapper function already more PYTHONIC.
===================
We could use even more intent attribute in the signature files to make the slycot._wrapper functions easier and save to use, and more PYTHONIC.
This would make massive changes in pyhton-control necessary.
The output
double precision intent(out) :: ab13bdin the .pyf file is there because ab13bd is function and not a subrountine. (see https://github.com/python-control/Slycot/pull/200#issuecomment-1644326484)
In the following example the optional attribute with default values is used heavily. The dimension n,m,p are determind by shape.
function ab13bd(dico,jobn,n,m,p,a,lda,b,ldb,c,ldc,d,ldd,nq,tol,dwork,ldwork,iwarn,info) ! in AB13BD.f
character optional :: dico = 'C'
character optional :: jobn = 'H'
integer optional :: n = shape(a,0)
integer optional :: m = shape(b,1)
integer optional :: p = shape(c,0)
double precision dimension(n,n), intent(in,out,copy) :: a
integer optional, depend(a) :: lda = shape(a,0)
double precision dimension(n,m), intent(in,out,copy) :: b
integer optional, depend(b) :: ldb = shape(b,0)
double precision dimension(p,n), intent(in,out,copy) :: c
integer optional, depend(c) :: ldc = shape(c,0)
double precision dimension(p,m), intent(in,out,copy) :: d
integer optional, depend(d) :: ldd = shape(d,0)
integer intent(out) :: nq
double precision optional :: tol = 0.0
double precision intent(hide,cache),dimension(ldwork),depend(ldwork) :: dwork
integer optional :: ldwork = max(1,max(m*(n+m)+max(n*(n+5),max(m*(m+2),4*p)),n*(max(n,p)+4)+min(n,p)))
integer intent(out) :: iwarn
integer intent(out) :: info
double precision intent(out) :: ab13bd
end function ab13bdThis would allow a minimal function call
slycot._wrapper.ab13bd(A,B,C,D)but also some optional attributes
out = slycot._wrapper.ab13bd(A, B, C, D, dico=dico, jobn=jobn, n=n, m=m, p=p, tol=tol)are possible.
KybernetikJo
commented
1 year ago To sum it up,
I would keep these commits and I would go for intent(in,out,copy) for the short term.
In the long run, we might reconsider the whole slycot API,
since the inner and outer wrappers for different procedures are quite heterogeneous.
Good templates could be a good solution.
KybernetikJo
commented
1 year ago FYI
A new api execution test with the optional attribute fac54678ca4990328138991920e0edf553fe3df9 7d7a5fae9a4bf6abc35478ab40853fe59ba60949 was successful.
This Draft/PR fixes the issues #198 and #199.
change parameter check #198
uses a shallow copy of state space matrices #199
Further testing is needed. Edge cases?
Unit tests (pytest) could be added?