Skip to content

Change to code for ctrb and obsv #300

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom
Merged
merged 2 commits into from
May 18, 2019
Merged

Change to code for ctrb and obsv #300

merged 2 commits into from
May 18, 2019

Conversation

billtubbs
Copy link
Contributor

@billtubbs billtubbs commented May 18, 2019 •
edited
Loading

This is not a bug fix. It's a minor speed improvement to statefbk.ctrb and statefbk.obsv.

[It's the first time I have contributed to someone else's project so welcome feedback and advice if I'm not doing this correctly]

Summary

I noticed that you can avoid multiple hstack operations by creating a list first.

I've done some testing and I think the behaviour of these methods is unchanged. See my test outputs below.

I realise this is a very minor improvement so feel free to ignore.

1. Testing ctrb()

Python 3.7.3 | packaged by conda-forge | (default, Mar 27 2019, 15:43:19) 
Type 'copyright', 'credits' or 'license' for more information
IPython 7.5.0 -- An enhanced Interactive Python. Type '?' for help.

In [1]: import numpy as np                                                      

In [2]: import control                                                          

In [3]: np.random.seed(1)                                                       

In [4]: cum = 0                                                                 

In [5]: for i in range(10): 
   ...:     A = np.random.randn(i, i) 
   ...:     B = np.random.randn(i, 1) 
   ...:     cum += control.ctrb(A, B).sum() 
   ...:                                                                         

In [6]: cum                                                                     
Out[6]: -10066.24518631401

In [7]: control.__file__                                                        
Out[7]: '/Users/billtubbs/python-control/python-control/control/__init__.py'

Results with official control package for comparison

Python 3.6.7 | packaged by conda-forge | (default, Feb 28 2019, 02:16:08) 
Type 'copyright', 'credits' or 'license' for more information
IPython 7.3.0 -- An enhanced Interactive Python. Type '?' for help.

In [1]: import numpy as np                                                                         

In [2]: import control                                                                             

In [3]: np.random.seed(1)                                                                          

In [4]: cum = 0.                                                                                   

In [5]: for i in range(10): 
   ...:     A = np.random.randn(i, i) 
   ...:     B = np.random.randn(i, 1) 
   ...:     cum += control.ctrb(A, B).sum() 
   ...:                                                                                            

In [6]: cum                                                                                        
Out[6]: -10066.245186314005

In [7]: control.__file__                                                                           
Out[7]: '/Users/billtubbs/anaconda3/envs/torch/lib/python3.6/site-packages/control/__init__.py'

2. Testing obsv()

In [1]: import numpy as np                                                      

In [2]: import control                                                          

In [3]: control.__file__                                                        
Out[3]: '/Users/billtubbs/python-control/python-control/control/__init__.py'

In [4]: np.random.seed(1)                                                       

In [5]: cum = 0.                                                                

In [6]: for i in range(10): 
   ...:     A = np.random.randn(i, i) 
   ...:     C = np.random.randn(1, i) 
   ...:     cum += control.obsv(A, C).sum() 
   ...:                                                                         

In [7]: cum                                                                     
Out[7]: -7589.3515984634905

Results with official control package for comparison

In [12]: np.random.seed(1)                                                                         

In [13]: cum = 0.                                                                                  

In [14]: for i in range(10): 
    ...:     A = np.random.randn(i, i) 
    ...:     C = np.random.randn(1, i) 
    ...:     cum += control.obsv(A, C).sum() 
    ...:                                                                                           

In [15]: cum                                                                                       
Out[15]: -7589.351598463492

In [16]: control.__file__                                                                          
Out[16]: '/Users/billtubbs/anaconda3/envs/torch/lib/python3.6/site-packages/control/__init__.py'

3. Speed Tests (In Python 3.7 on a 2013 MacBook Pro)

In [8]: def ctrb1(A, B):  
   ...:      amat = np.mat(A)  
   ...:      bmat = np.mat(B)  
   ...:      n = np.shape(amat)[0]  
   ...:      ctrb = bmat  
   ...:      for i in range(1, n):  
   ...:          ctrb = np.hstack((ctrb, amat**i*bmat))  
   ...:      return ctrb 
   ...:                                                                         

In [9]: def ctrb2(A, B):  
   ...:      amat = np.mat(A)  
   ...:      bmat = np.mat(B)  
   ...:      n = np.shape(amat)[0]  
   ...:      ctrb = np.hstack([bmat] + [amat**i*bmat for i in range(1, n)])  
   ...:      return ctrb 
   ...:                                                                         

In [11]: np.random.seed(1)                                                      

In [12]: A = np.random.randn(8, 8)                                              

In [13]: B = np.random.randn(8, 1)                                              

In [14]: %timeit ctrb1(A, B)                                                    
456 µs ± 18 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

In [15]: %timeit ctrb2(A, B)                                                    
378 µs ± 13.2 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

In [16]: 1-378/456                                                              
Out[16]: 0.17105263157894735

@murrayrm
Copy link
Member

murrayrm commented May 18, 2019 •
edited
Loading

@billtubbs Thanks for the contribution to python-control!

Your proposed changes are not passing the Travis CI unit tests. Can you fix those and then resubmit? The problem could be in your code or could be in the unit tests, but something is not right.

Also, I wonder whether it might not be faster by using the fact that each iteration of the computation just needs to pre-/post-multiple the last col/row of the matrix by A, rather than computing A**i.

@billtubbs
Copy link
Contributor Author

I see I made an error. I checked out your unit tests and realised my mistake. I will make the change shortly and test it with your tests this time...

Sorry about that.

@billtubbs
Copy link
Contributor Author

I fixed the error and it seems to pass the unit tests now:

In [1]: from control.tests.statefbk_test import TestStatefbk                    

In [2]: tests = TestStatefbk()                                                  

In [3]: tests.testCtrbSISO()                                                    

In [4]: tests.testCtrbMIMO()                                                    

In [5]: tests.testObsvSISO()                                                    

In [6]: tests.testObsvMIMO()

@coveralls
Copy link

Coverage Status

Coverage decreased (-0.02%) to 78.205% when pulling 9309259 on billtubbs:bt-statefbk into c0092e0 on python-control:master.

@murrayrm
Copy link
Member

Changes look good. More pythonic and if we get some speed up as well, that's great.

@murrayrm murrayrm merged commit ace1683 into python-control:master May 18, 2019
@billtubbs
Copy link
Contributor Author

Cool, thanks. It's only a minor speed improvement of 13-17% on something that wasn't exactly computationally demanding...

Anyhow, I am interested in this project. If there is any way I can help let me know. I'd say I'm an intermediate/good python programmer and currently learning control engineering.

@billtubbs billtubbs deleted the bt-statefbk branch May 19, 2019 00:03
@billtubbs
Copy link
Contributor Author

Also, I wonder whether it might not be faster by using the fact that each iteration of the computation just needs to pre-/post-multiple the last col/row of the matrix by A, rather than computing A**i.

Ah yes, good idea. I'll look into that.

@murrayrm
Copy link
Member

Places that could use help: https://github.com/python-control/python-control/issues

@billtubbs
Copy link
Contributor Author

billtubbs commented May 19, 2019 •
edited
Loading

I worked on some faster versions. Turns out converting the arguments to np.mat was wasting a lot of time. However, the fastest two solutions are not that pretty and hardly as readable as ctrb2.

What do you think?

Setup:

In [1]: import numpy as np                                                      

In [2]: np.random.seed(1)                                                       

In [3]: A = np.random.randn(8, 8)                                               

In [4]: B = np.random.randn(8, 1)

In [47]: def ctrb1(A, B):  # Original
    ...:     amat = np.mat(A) 
    ...:     bmat = np.mat(B) 
    ...:     n = np.shape(amat)[0] 
    ...:     ctrb = bmat 
    ...:     for i in range(1, n): 
    ...:         ctrb = np.hstack((ctrb, amat**i*bmat)) 
    ...:     return ctrb 
    ...:                                                                        

In [48]: %timeit ctrb1(A, B)   # Most recent improvement                                          
227 µs ± 24.8 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

In [49]: def ctrb2(A, B):  
    ...:      amat = np.mat(A)  
    ...:      bmat = np.mat(B)  
    ...:      n = np.shape(amat)[0]  
    ...:      ctrb = np.hstack([bmat] + [amat**i*bmat for i in range(1, n)])  
    ...:      return ctrb 
    ...:                                                                        

In [50]: %timeit ctrb2(A, B)                                                    
175 µs ± 783 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)

In [51]: def ctrb3(A, B):   # Eliminating np.mat conversion
    ...:     n = A.shape[0] 
    ...:     ctrb = np.hstack([B] + [np.linalg.matrix_power(A, i).dot(B) for i i
    ...: n range(1, n)]) 
    ...:     return ctrb 
    ...:                                                                        

In [52]: %timeit ctrb3(A, B)                                                    
70.5 µs ± 690 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)

In [53]: def ctrb4(A, B):    # Eliminating the matrix power operation
    ...:     n = A.shape[0] 
    ...:     items = [B] 
    ...:     x = np.eye(n) 
    ...:     for i in range(1, n): 
    ...:         x = x.dot(A) 
    ...:         items.append(x.dot(B)) 
    ...:     ctrb = np.hstack(items) 
    ...:     return ctrb 
    ...:                                                                        

In [54]: %timeit ctrb4(A, B)                                                    
26.6 µs ± 212 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)

In [56]: np.isclose(ctrb1(A, B), ctrb3(A, B)).all()                             
Out[56]: True

In [57]: np.isclose(ctrb1(A, B), ctrb4(A, B)).all()                             
Out[57]: True

In [58]: def g(A, B):   # Using a generator function
    ...:     n = A.shape[0] 
    ...:     x = np.eye(n) 
    ...:     yield B 
    ...:     for i in range(1, n):  
    ...:         x = x.dot(A) 
    ...:         yield x.dot(B) 
    ...:                                                                        

In [61]: def ctrb5(A, B): 
    ...:     return np.hstack(tuple(g(A, B))) 
    ...:                                                                        

In [63]: %timeit ctrb5(A, B)                                                    
27.1 µs ± 147 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)

In [64]: np.isclose(ctrb1(A, B), ctrb5(A, B)).all() 
Out[64]: True

ctrb4() looks like the best for speed.

@murrayrm
Copy link
Member

The numpy matrix class is being deprecated, so eventually the np.mat will have to go away. Check out PR #292 for a refactored version and perhaps build off of that?

The PR is now merged, so if you want to make additional changes it will have to go in a new PR.

@billtubbs
Copy link
Contributor Author

billtubbs commented May 19, 2019 •
edited
Loading

Understood. Seems like a good idea to replace np.mat but looks like a big job...

I'll have to add things like the following to ctrb and obsv:

def ctrb(A, B):
    amat = np.array(A, ndmin=2)
    bmat = np.array(B, ndmin=2) 
    ...

or even:

def ctrb(A, B):
    if isinstance(A, str):
        A = np.mat(A)
    if isinstance(B, str):
        B = np.mat(B)
    amat = np.array(A, ndmin=2, dtype=float)
    bmat = np.array(B, dtype=float) 
    ...

I'll maybe work on this next weekend and make sure it's working before starting a new PR.

@murrayrm
Copy link
Member

Current changes are better than what was there. OK to keep tweaking, but there are lots of other things to work on that could use work, if you have time.

@murrayrm murrayrm mentioned this pull request May 19, 2019
Closed
@murrayrm murrayrm added this to the 0.8.3 milestone Jan 4, 2020
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Reviewers
No reviews
Assignees
No one assigned
Labels
None yet
Projects
None yet
Milestone
0.8.3
Development

Successfully merging this pull request may close these issues.
3 participants
@billtubbs @murrayrm @coveralls