Discrete-time additions to scipy.signal
![](https://secure.gravatar.com/avatar/c748a4ea7c8f04551253d59fcb96683a.jpg?s=120&d=mm&r=g)
Hi folks, I've been working on some additions to the scipy.signal module. A significant portion of my work involves dealing with linear systems in the discrete-time domain, but I found that most of the functionality I would need in scipy is implemented only in the continuous domain. I've therefore added some code to handle the discrete time cases. My fork is visible at: https://github.com/ArmstrongJ/scipy Additionally, I presented at PyCon 2011, and I briefly discussed a discrete algebraic Riccati equation solver I had implemented in Python. I've included that code into scipy.signal along with a continuous implementation. I'd appreciate some comments and review of the code. Everything seems to be working for me at this point, but any constructive criticism is greatly appreciated. -Jeff Jeff Armstrong - jba@sdf.lonestar.org SDF Public Access UNIX System - http://sdf.lonestar.org
![](https://secure.gravatar.com/avatar/3d3176cf99cae23d0ac119d1ea6c4d11.jpg?s=120&d=mm&r=g)
Hi Jeffrey, On Tue, Apr 26, 2011 at 2:40 PM, Jeffrey Armstrong <jba@sdf.lonestar.org> wrote:
Hi folks, I've been working on some additions to the scipy.signal module. A significant portion of my work involves dealing with linear systems in the discrete-time domain, but I found that most of the functionality I would need in scipy is implemented only in the continuous domain. I've therefore added some code to handle the discrete time cases. My fork is visible at:
https://github.com/ArmstrongJ/scipy
Additionally, I presented at PyCon 2011, and I briefly discussed a discrete algebraic Riccati equation solver I had implemented in Python. I've included that code into scipy.signal along with a continuous implementation.
I'd appreciate some comments and review of the code. Everything seems to be working for me at this point, but any constructive criticism is greatly appreciated.
I looked through your changes and at first glance it looks pretty good to me. It's not so easy to review however, because it's all in your master branch, it's a lot of code, and there is quite some code added that is later deleted again. It would be easier if you would create separate branches for separate features (not your master branch, keep that a clean mirror of the numpy master branch). There are at least two, the discrete versions of features already present for the continuous domain, and your additions of Riccati, Lyapunov and Sylvester solvers. Perhaps the sort kw for decomp_schur is a third. The current names of modules are not very descriptive. They could be changed to something like: c2d -> cont2discrete care+dare -> riccati lyap+dlyap -> lyapunov Or maybe put it in a "control" module. Can you explain the relation to pydare a bit? Is this code all from pydare and are you relicensing it as BSD and proposing it for inclusion, or is part of this new code? Cheers, Ralf
![](https://secure.gravatar.com/avatar/c748a4ea7c8f04551253d59fcb96683a.jpg?s=120&d=mm&r=g)
On Tue, 26 Apr 2011, Ralf Gommers wrote:
Hi Jeffrey,
I looked through your changes and at first glance it looks pretty good to me. It's not so easy to review however, because it's all in your master branch, it's a lot of code, and there is quite some code added that is later deleted again. It would be easier if you would create separate branches for separate features (not your master branch, keep that a clean mirror of the numpy master branch). There are at least two, the discrete versions of features already present for the continuous domain, and your additions of Riccati, Lyapunov and Sylvester solvers. Perhaps the sort kw for decomp_schur is a third.
Ralf, I apologize for the lack of branches. I must confess that this is my first large-scale experience in a git environment so I'm still getting the hang of it. You've summarized the changes pretty well: 1. Duplicating most of the continuous functionality in the discrete domain (dltisys.py, c2d.py) 2. Adding Lyapunov and algebraic Riccati equation solvers (lyap.py, dlyap.py, dare.py, care.py, decomp_schur.py, trsyl flapack wrapper). The sort kw in decomp_schur is necessary to properly solve the algebraic Riccati equation.
The current names of modules are not very descriptive. They could be changed to something like: c2d -> cont2discrete care+dare -> riccati lyap+dlyap -> lyapunov Or maybe put it in a "control" module.
The module names could be more descriptive, I agree. While I think the addition of a "control" module in scipy would perhaps be a good idea, I didn't have the audacity to make such a suggestion from the start.
Can you explain the relation to pydare a bit? Is this code all from pydare and are you relicensing it as BSD and proposing it for inclusion, or is part of this new code?
The pydare package (http://code.google.com/p/pydare/) is a small package providing discrete algebraic Riccati and Lyapunov equation solvers, and I am the sole author. The code was originally licensed under GPLv3 because the discrete Lyapunov equation solver utilized code that had been directly translated from some Octave code into Python, which was also GPL. The pydare package can also exploit SLICOT, another GPL'd package. To contribute the code to SciPy, I do indeed relicense the entirety of my code as BSD-licensed, and I'll assign copyright or whatever as necessary. The GPL-based code in the discrete Lyapunov equation solver has been redacted in its entirety, replaced by original code that relies on the solution of a Sylvester equation instead (requiring an inversion, which is less than ideal). All references to SLICOT have also been removed, and there is not any code present any longer that is not original. Jeff Armstrong - jba@sdf.lonestar.org SDF Public Access UNIX System - http://sdf.lonestar.org
![](https://secure.gravatar.com/avatar/3d3176cf99cae23d0ac119d1ea6c4d11.jpg?s=120&d=mm&r=g)
On Wed, Apr 27, 2011 at 1:28 PM, Jeffrey Armstrong <jba@sdf.lonestar.org> wrote:
On Tue, 26 Apr 2011, Ralf Gommers wrote:
Hi Jeffrey,
I looked through your changes and at first glance it looks pretty good to me. It's not so easy to review however, because it's all in your master branch, it's a lot of code, and there is quite some code added that is later deleted again. It would be easier if you would create separate branches for separate features (not your master branch, keep that a clean mirror of the numpy master branch). There are at least two, the discrete versions of features already present for the continuous domain, and your additions of Riccati, Lyapunov and Sylvester solvers. Perhaps the sort kw for decomp_schur is a third.
Ralf,
I apologize for the lack of branches. I must confess that this is my first large-scale experience in a git environment so I'm still getting the hang of it.
Are you trying to do this? If you need any help let me know (offline perhaps).
You've summarized the changes pretty well:
1. Duplicating most of the continuous functionality in the discrete domain (dltisys.py, c2d.py)
2. Adding Lyapunov and algebraic Riccati equation solvers (lyap.py, dlyap.py, dare.py, care.py, decomp_schur.py, trsyl flapack wrapper).
The sort kw in decomp_schur is necessary to properly solve the algebraic Riccati equation.
(1) looks like a good addition, and fits in well in scipy.signal. (2) could use some discussion, if and where it belongs in scipy. It looks like these equations are fairly widely used in control theory (which I don't know all that much about), but they're a bit more exotic/inaccessible than your average scipy function. Cheers, Ralf
The current names of modules are not very descriptive. They could be changed to something like: c2d -> cont2discrete care+dare -> riccati lyap+dlyap -> lyapunov Or maybe put it in a "control" module.
The module names could be more descriptive, I agree. While I think the addition of a "control" module in scipy would perhaps be a good idea, I didn't have the audacity to make such a suggestion from the start.
Can you explain the relation to pydare a bit? Is this code all from pydare and are you relicensing it as BSD and proposing it for inclusion, or is part of this new code?
The pydare package (http://code.google.com/p/pydare/) is a small package providing discrete algebraic Riccati and Lyapunov equation solvers, and I am the sole author. The code was originally licensed under GPLv3 because the discrete Lyapunov equation solver utilized code that had been directly translated from some Octave code into Python, which was also GPL. The pydare package can also exploit SLICOT, another GPL'd package.
To contribute the code to SciPy, I do indeed relicense the entirety of my code as BSD-licensed, and I'll assign copyright or whatever as necessary. The GPL-based code in the discrete Lyapunov equation solver has been redacted in its entirety, replaced by original code that relies on the solution of a Sylvester equation instead (requiring an inversion, which is less than ideal). All references to SLICOT have also been removed, and there is not any code present any longer that is not original.
Jeff Armstrong - jba@sdf.lonestar.org SDF Public Access UNIX System - http://sdf.lonestar.org _______________________________________________ SciPy-Dev mailing list SciPy-Dev@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-dev
![](https://secure.gravatar.com/avatar/c748a4ea7c8f04551253d59fcb96683a.jpg?s=120&d=mm&r=g)
On Thu, 28 Apr 2011, Ralf Gommers wrote:
Are you trying to do this? If you need any help let me know (offline perhaps).
I think I've successfully split the changes into two branches: discrete time additions: https://github.com/ArmstrongJ/scipy/tree/discrete-time-signal Lyapunov and algebraic Riccati solvers: https://github.com/ArmstrongJ/scipy/tree/riccati-lyapunov There could possibly be some minor mess with re-committing __init__.py in scipy.signal, but the final product is ok. Git has taken a bit of time to understand.
(1) looks like a good addition, and fits in well in scipy.signal.
Great! I'm glad to hear it.
(2) could use some discussion, if and where it belongs in scipy. It looks like these equations are fairly widely used in control theory (which I don't know all that much about), but they're a bit more exotic/inaccessible than your average scipy function.
A controls module might indeed be a good place for these to live. Octave (and MATLAB for that matter) keeps them in a controls toolbox, but the linear system functionality in scipy.signal is also present in these toolboxes. Maybe these solvers along with the scipy.signal functionality could provide an interesting place to start on a controls module (scipy.control). I do, however, think that my change introducing the 'sort' keyword to scipy.linalg.schur is a valuable addition. Being able to sort the eigenvalues in some manner is quite useful in a number of situations. This change is present in my riccati-lyapunov branch. -Jeff Jeff Armstrong - jba@sdf.lonestar.org SDF Public Access UNIX System - http://sdf.lonestar.org
![](https://secure.gravatar.com/avatar/ad13088a623822caf74e635a68a55eae.jpg?s=120&d=mm&r=g)
On Fri, Apr 29, 2011 at 8:39 AM, Jeffrey Armstrong <jba@sdf.lonestar.org> wrote:
On Thu, 28 Apr 2011, Ralf Gommers wrote:
Are you trying to do this? If you need any help let me know (offline perhaps).
I think I've successfully split the changes into two branches:
discrete time additions: https://github.com/ArmstrongJ/scipy/tree/discrete-time-signal
Lyapunov and algebraic Riccati solvers: https://github.com/ArmstrongJ/scipy/tree/riccati-lyapunov
There could possibly be some minor mess with re-committing __init__.py in scipy.signal, but the final product is ok. Git has taken a bit of time to understand.
(1) looks like a good addition, and fits in well in scipy.signal.
Great! I'm glad to hear it.
(2) could use some discussion, if and where it belongs in scipy. It looks like these equations are fairly widely used in control theory (which I don't know all that much about), but they're a bit more exotic/inaccessible than your average scipy function.
A controls module might indeed be a good place for these to live. Octave (and MATLAB for that matter) keeps them in a controls toolbox, but the linear system functionality in scipy.signal is also present in these toolboxes. Maybe these solvers along with the scipy.signal functionality could provide an interesting place to start on a controls module (scipy.control).
I do, however, think that my change introducing the 'sort' keyword to scipy.linalg.schur is a valuable addition. Being able to sort the eigenvalues in some manner is quite useful in a number of situations. This change is present in my riccati-lyapunov branch.
I browsed the source a bit but haven't tried it. a few comments: pep8 : space after comma in function arguments, class CareSolver: should subclass object for python 2.4, 2.5 compatibility Is mixing arrays and matrices a good idea? In general it's confusing.
From reading the code CareSolver uses matrices, but I didn't see whether it also uses arrays. The continuous and discrete system should use the same matrix/array pattern, and as far as I remember ltisys uses arrays only.
Since this has been a wished for feature on the mailing list for a long time, it is nice to see a concrete proposal to enhance scipy. Josef
-Jeff
Jeff Armstrong - jba@sdf.lonestar.org SDF Public Access UNIX System - http://sdf.lonestar.org _______________________________________________ SciPy-Dev mailing list SciPy-Dev@scipy.org http://mail.scipy.org/mailman/listinfo/scipy-dev
![](https://secure.gravatar.com/avatar/3d3176cf99cae23d0ac119d1ea6c4d11.jpg?s=120&d=mm&r=g)
On Fri, Apr 29, 2011 at 2:39 PM, Jeffrey Armstrong <jba@sdf.lonestar.org> wrote:
On Thu, 28 Apr 2011, Ralf Gommers wrote:
Are you trying to do this? If you need any help let me know (offline perhaps).
I think I've successfully split the changes into two branches:
discrete time additions: https://github.com/ArmstrongJ/scipy/tree/discrete-time-signal
Lyapunov and algebraic Riccati solvers: https://github.com/ArmstrongJ/scipy/tree/riccati-lyapunov
There could possibly be some minor mess with re-committing __init__.py in scipy.signal, but the final product is ok. Git has taken a bit of time to understand.
That was a good start. I've taken your discrete-time-signal branch, and did some more reorganizing to only have a few commits with no renames and deletions: https://github.com/rgommers/scipy/tree/armstrong-discrete. I then did some cleanups of cont2discrete in the last commit, to make the code cleaner. Some comments on your code: - function signatures with *args and **kwargs are usually a bad idea, it is better to have separate functions or some other way to make clear what the expected inputs are. (this is why I split up your cont2discrete function). - like Josef said, it should conform to PEP8. - I also agree with him that it's better to not mix arrays and matrices, and preferably avoid matrices completely. - the standard numpy/scipy conventions for imports, docstrings etc. should be followed, like "import numpy as np". See also https://github.com/numpy/numpy/blob/master/doc/HOWTO_DOCUMENT.rst.txt - the permissions of your test files was wrong, they were executable. Nose will ignore such files, so your tests were never executed. - the caps in the test files are a little jarring to read. They can be lowercase, and moved into the test class. Only the continuous constants are reused, the discrete ones can be in the actual test functions where they are used. In the cont2discrete file, it would also be helpful if you added some comments (or references) to the least transparent parts. For example, I can't tell what the like with all the multiple vstack/hstack's is doing. I propose you pull my branch of your work and make changes to that. Later the commits can be rewritten again to keep it easy to review. Cheers, Ralf
![](https://secure.gravatar.com/avatar/c748a4ea7c8f04551253d59fcb96683a.jpg?s=120&d=mm&r=g)
On Sat, 30 Apr 2011, Ralf Gommers wrote:
That was a good start. I've taken your discrete-time-signal branch, and did some more reorganizing to only have a few commits with no renames and deletions: https://github.com/rgommers/scipy/tree/armstrong-discrete. I then did some cleanups of cont2discrete in the last commit, to make the code cleaner.
I'll pull from here for any future work. The cont2discrete function needs some updates. I have to admit I don't remember why I chose to use *args on this function, but, for consistency with other LTI functions, it should probably be combined into a single function again, but accept a tuple as the system. The tuple's length would then be checked to see if it's a state-space model or a transfer function. In regards to the permissions, I've been working on Windows, so the tests were in fact executing on my end. After reading your email, it did occur to me that git was marking new files a "755", which was unexpected. I'll make sure it doesn't happen again. As far as styling goes, I'll do my best to stick with best practices. Sorry for any inconvenience. -Jeff Jeff Armstrong - jba@sdf.lonestar.org SDF Public Access UNIX System - http://sdf.lonestar.org
![](https://secure.gravatar.com/avatar/3d3176cf99cae23d0ac119d1ea6c4d11.jpg?s=120&d=mm&r=g)
On Mon, May 2, 2011 at 2:03 PM, Jeffrey Armstrong <jba@sdf.lonestar.org> wrote:
On Sat, 30 Apr 2011, Ralf Gommers wrote:
That was a good start. I've taken your discrete-time-signal branch, and did some more reorganizing to only have a few commits with no renames and deletions: https://github.com/rgommers/scipy/tree/armstrong-discrete. I then did some cleanups of cont2discrete in the last commit, to make the code cleaner.
I'll pull from here for any future work. The cont2discrete function needs some updates. I have to admit I don't remember why I chose to use *args on this function, but, for consistency with other LTI functions, it should probably be combined into a single function again, but accept a tuple as the system. The tuple's length would then be checked to see if it's a state-space model or a transfer function.
Not sure I understand that reasoning - changing a transfer function or a state space model look like separate operations to me.
In regards to the permissions, I've been working on Windows, so the tests were in fact executing on my end. After reading your email, it did occur to me that git was marking new files a "755", which was unexpected. I'll make sure it doesn't happen again.
As far as styling goes, I'll do my best to stick with best practices. Sorry for any inconvenience.
No problem at all. Cheers, Ralf
![](https://secure.gravatar.com/avatar/c748a4ea7c8f04551253d59fcb96683a.jpg?s=120&d=mm&r=g)
On Mon, 2 May 2011, Ralf Gommers wrote:
I'll pull from here for any future work. The cont2discrete function needs some updates. I have to admit I don't remember why I chose to use *args on this function, but, for consistency with other LTI functions, it should probably be combined into a single function again, but accept a tuple as the system. The tuple's length would then be checked to see if it's a state-space model or a transfer function.
Not sure I understand that reasoning - changing a transfer function or a state space model look like separate operations to me.
While the two operations are mildly different, most of the LTI functions in scipy.signal accept a tuple and determine whether they are dealing with a transfer function or a state-space model internally. If you look at scipy.signal.lsim, for example, you'll see the first input is the "system," which can be one of three types. For consistency, I would say that the cont2discrete function should behave similarly. I'm not implying its the correct route to take, only that it seems that it was handled this way in the past. -Jeff Jeff Armstrong - jba@sdf.lonestar.org SDF Public Access UNIX System - http://sdf.lonestar.org
participants (3)
-
Jeffrey Armstrong
-
josef.pktd@gmail.com
-
Ralf Gommers