Create config options for using numpy.matrix as state space return type

murrayrm · murrayrm · commit dc333a6e6ce6 · 2019-06-07T22:16:14.000-07:00
diff --git a/control/config.py b/control/config.py
@@ -7,13 +7,20 @@
 # files.  For now, you can just choose between MATLAB and FBS default
 # values.
 
+import warnings
+import numpy as np
+from .statesp import StateSpaceMatrix
+
 # Bode plot defaults
 bode_dB = False                 # Bode plot magnitude units
 bode_deg = True                 # Bode Plot phase units
 bode_Hz = False                 # Bode plot frequency units
 bode_number_of_samples = None   # Bode plot number of samples
 bode_feature_periphery_decade = 1.0  # Bode plot feature periphery in decades
 
+# State space return type (change to StateSpaceMatrix at next major revision)
+ss_return_type = np.matrix
+
 
 def reset_defaults():
     """Reset configuration values to their default values."""
@@ -22,6 +29,7 @@ def reset_defaults():
     global bode_Hz; bode_Hz = False
     global bode_number_of_samples; bode_number_of_samples = None
     global bode_feature_periphery_decade; bode_feature_periphery_decade = 1.0
+    global ss_return_type; ss_return_type = np.matrix   # TODO: update
 
 
 # Set defaults to match MATLAB
@@ -52,3 +60,49 @@ def use_fbs_defaults():
     global bode_deg; bode_deg = True
     global bode_Hz; bode_Hz = False
 
+#
+# State space function return type
+#
+# These functions are used to set the return type for state space functions
+# that return a matrix.  In the original version of python-control, these
+# functions returned a numpy.matrix object.  To handle the eventual
+# deprecation of the numpy.matrix type, the StateSpaceMatrix object type,
+# which implements matrix multiplications and related numpy.matrix operations
+# was created.  To avoid breaking existing code, the return type from
+# functions that used to return an np.matrix object now call the
+# get_ss_return_type function to obtain the return type.  The return type can
+# also be set using the `return_type` keyword, overriding the default state
+# space matrix return type.
+#
+
+
+# Set state space return type
+def set_ss_return_type(subtype, warn=True):
+    global ss_return_type
+    # If return_type is np.matrix, issue a pending deprecation warning
+    if (subtype is np.matrix and warn):
+        warnings.warn("Return type numpy.matrix is soon to be deprecated.",
+                      stacklevel=2)
+    ss_return_type = subtype
+
+
+# Get the state space return type
+def get_ss_return_type(subtype=None, warn=True):
+    global ss_return_type
+    return_type = ss_return_type if subtype is None else subtype
+    # If return_type is np.matrix, issue a pending deprecation warning
+    if (return_type is np.matrix and warn):
+        warnings.warn("Returning numpy.matrix, soon to be deprecated; "
+                      "make sure calling code can handle nparray.",
+                      stacklevel=2)
+    return return_type
+
+
+# Function to turn on/off use of np.matrix type
+def use_numpy_matrix(flag=True, warn=True):
+    if flag and warn:
+        warnings.warn("Return type numpy.matrix is soon to be deprecated.",
+                      stacklevel=2)
+        set_ss_return_type(np.matrix, warn=False)
+    else:
+        set_ss_return_type(StateSpaceMatrix)
diff --git a/control/modelsimp.py b/control/modelsimp.py
@@ -45,19 +45,19 @@
 
 # External packages and modules
 import numpy as np
-import warnings
 from .exception import ControlSlycot
 from .lti import isdtime, isctime
 from .statesp import StateSpace
 from .statefbk import gram
+from .config import get_ss_return_type
 
 __all__ = ['hsvd', 'balred', 'modred', 'era', 'markov', 'minreal']
 
 # Hankel Singular Value Decomposition
 #   The following returns the Hankel singular values, which are singular values
 #of the matrix formed by multiplying the controllability and observability
 #grammians
-def hsvd(sys, return_type=np.matrix):
+def hsvd(sys, return_type=None):
     """Calculate the Hankel singular values.
 
     Parameters
@@ -90,12 +90,6 @@ def hsvd(sys, return_type=np.matrix):
     >>> H = hsvd(sys)
 
     """
-    # If return_type is np.matrix, issue a pending deprecation warning
-    if (return_type is np.matrix):
-        warnings.warn("Returning numpy.matrix, soon to be deprecated; "
-                      "make sure calling code can handle nparray.",
-                      stacklevel=2)
-
     # TODO: implement for discrete time systems
     if (isdtime(sys, strict=True)):
         raise NotImplementedError("Function not implemented in discrete time")
@@ -111,7 +105,7 @@ def hsvd(sys, return_type=np.matrix):
     hsv = np.fliplr(hsv)
 
     # Return the Hankel singular values (casting type, if needed)
-    return hsv.view(type=return_type)
+    return hsv.view(type=get_ss_return_type(return_type))
 
 def modred(sys, ELIM, method='matchdc'):
     """
diff --git a/control/statefbk.py b/control/statefbk.py
@@ -42,8 +42,8 @@
 # External packages and modules
 import numpy as np
 import scipy as sp
-import warnings
 from . import statesp
+from .config import get_ss_return_type
 from .exception import ControlSlycot, ControlArgument, ControlDimension
 
 __all__ = ['ctrb', 'obsv', 'gram', 'place', 'place_varga', 'lqr', 'acker']
@@ -220,7 +220,7 @@ def place_varga(A, B, p, dtime=False, alpha=None):
     return -F
 
 # Contributed by Roberto Bucher <roberto.bucher@supsi.ch>
-def acker(A, B, poles, return_type=np.matrix):
+def acker(A, B, poles, return_type=None):
     """Pole placement using Ackermann method
 
     Call:
@@ -239,12 +239,6 @@ def acker(A, B, poles, return_type=np.matrix):
         Gains such that A - B K has given eigenvalues
 
     """
-    # If return_type is np.matrix, issue a pending deprecation warning
-    if (return_type is np.matrix):
-        warnings.warn("Returning numpy.matrix, soon to be deprecated; "
-                      "make sure calling code can handle nparray.",
-                      stacklevel=2)
-
     # Convert the inputs to matrices (arrays)
     a = statesp.ssmatrix(A)
     b = statesp.ssmatrix(B)
@@ -265,7 +259,7 @@ def acker(A, B, poles, return_type=np.matrix):
     K = np.linalg.solve(ct, pmat)
 
     K = K[-1][:]                # Extract the last row
-    return K.view(type=return_type)
+    return K.view(type=get_ss_return_type(return_type))
 
 def lqr(*args, **keywords):
     """lqr(A, B, Q, R[, N])
@@ -375,7 +369,7 @@ def lqr(*args, **keywords):
 
     return K, S, E
 
-def ctrb(A, B, return_type=np.matrix):
+def ctrb(A, B, return_type=None):
     """Controllabilty matrix
 
     Parameters
@@ -396,12 +390,6 @@ def ctrb(A, B, return_type=np.matrix):
     >>> C = ctrb(A, B)
 
     """
-    # If return_type is np.matrix, issue a pending deprecation warning
-    if (return_type is np.matrix):
-        warnings.warn("Returning numpy.matrix, soon to be deprecated; "
-                      "make sure calling code can handle nparray.",
-                      stacklevel=2)
-
     # Convert input parameters to matrices (if they aren't already)
     amat = statesp.ssmatrix(A)
     bmat = statesp.ssmatrix(B)
@@ -412,10 +400,10 @@ def ctrb(A, B, return_type=np.matrix):
                                for i in range(1, n)])
 
     # Return the observability matrix in the desired type
-    return ctrb.view(type=return_type)
+    return ctrb.view(type=get_ss_return_type(return_type))
 
 
-def obsv(A, C, return_type=np.matrix):
+def obsv(A, C, return_type=None):
     """Observability matrix
 
     Parameters
@@ -436,12 +424,6 @@ def obsv(A, C, return_type=np.matrix):
     >>> O = obsv(A, C)
 
     """
-    # If return_type is np.matrix, issue a pending deprecation warning
-    if (return_type is np.matrix):
-        warnings.warn("Returning numpy.matrix, soon to be deprecated; "
-                      "make sure calling code can handle nparray.",
-                      stacklevel=2)
-
     # Convert input parameters to matrices (if they aren't already)
     amat = statesp.ssmatrix(A)
     cmat = statesp.ssmatrix(C)
@@ -452,10 +434,10 @@ def obsv(A, C, return_type=np.matrix):
                                for i in range(1, n)])
 
     # Return the observability matrix in the desired type
-    return obsv.view(type=return_type)
+    return obsv.view(type=get_ss_return_type(return_type))
 
 
-def gram(sys, type, return_type=np.matrix):
+def gram(sys, type, return_type=None):
     """Gramian (controllability or observability)
 
     Parameters
@@ -492,12 +474,6 @@ def gram(sys, type, return_type=np.matrix):
     >>> Ro = gram(sys,'of'), where Wo=Ro'*Ro
 
     """
-    # If return_type is np.matrix, issue a pending deprecation warning
-    if (return_type is np.matrix):
-        warnings.warn("Returning numpy.matrix, soon to be deprecated; "
-                      "make sure calling code can handle nparray.",
-                      stacklevel=2)
-
     #Check for ss system object
     if not isinstance(sys,statesp.StateSpace):
         raise ValueError("System must be StateSpace!")
@@ -535,7 +511,7 @@ def gram(sys, type, return_type=np.matrix):
         A = np.array(sys.A)         # convert to NumPy array for slycot
         X,scale,sep,ferr,w = sb03md(n, C, A, U, dico, job='X', fact='N', trana=tra)
         gram = X
-        return gram.view(type=return_type)
+        return gram.view(type=get_ss_return_type(return_type))
 
     elif type=='cf' or type=='of':
         #Compute cholesky factored gramian from slycot routine sb03od
@@ -558,4 +534,4 @@ def gram(sys, type, return_type=np.matrix):
             C[0:n,0:m] = sys.C.transpose()
             X,scale,w = sb03od(n, m, A, Q, C.transpose(), dico, fact='N', trans=tra)
         gram = X
-        return gram.view(type=return_type)
+        return gram.view(type=get_ss_return_type(return_type))
diff --git a/control/tests/statefbk_test.py b/control/tests/statefbk_test.py
@@ -11,6 +11,8 @@
 from control.matlab import *
 from control.exception import slycot_check, ControlDimension
 from control.mateqn import care, dare
+from control.config import use_numpy_matrix
+from control.statesp import StateSpaceMatrix
 
 class TestStatefbk(unittest.TestCase):
     """Test state feedback functions"""
@@ -29,11 +31,15 @@ def testCtrbSISO(self):
         A = np.array([[1., 2.], [3., 4.]])
         B = np.array([[5.], [7.]])
         Wctrue = np.array([[5., 19.], [7., 43.]])
-        Wc = ctrb(A, B, return_type=np.ndarray)
+
+        use_numpy_matrix(False)
+        Wc = ctrb(A, B)
         np.testing.assert_array_almost_equal(Wc, Wctrue)
+        self.assertTrue(isinstance(Wc, StateSpaceMatrix))
 
-        # Check that default type generates a warning
+        # Check that default using np.matrix generates a warning
         # TODO: remove this check with matrix type is deprecated
+        use_numpy_matrix(True)
         with warnings.catch_warnings(record=True) as w:
             Wc = ctrb(A, B)
             self.assertTrue(issubclass(w[-1].category, UserWarning))
@@ -53,6 +59,7 @@ def testObsvSISO(self):
         A = np.array([[1., 2.], [3., 4.]])
         C = np.array([[5., 7.]])
         Wotrue = np.array([[5., 7.], [26., 38.]])
+        use_numpy_matrix(True)  # set default, but override in next statement
         Wo = obsv(A, C, return_type=np.ndarray)
         np.testing.assert_array_almost_equal(Wo, Wotrue)
 
@@ -61,6 +68,7 @@ def testObsvSISO(self):
 
         # Check that default type generates a warning
         # TODO: remove this check with matrix type is deprecated
+        use_numpy_matrix(True)
         with warnings.catch_warnings(record=True) as w:
             Wo = obsv(A, C)
             self.assertTrue(issubclass(w[-1].category, UserWarning))
@@ -95,6 +103,7 @@ def testGramWc(self):
 
         # Check that default type generates a warning
         # TODO: remove this check with matrix type is deprecated
+        use_numpy_matrix(True)
         with warnings.catch_warnings(record=True) as w:
             Wc = gram(sys, 'c')
             self.assertTrue(issubclass(w[-1].category, UserWarning))
@@ -359,7 +368,10 @@ def test_dare(self):
 
 
 def test_suite():
-   return unittest.TestLoader().loadTestsFromTestCase(TestStatefbk)
+    
+   status1 = unittest.TestLoader().loadTestsFromTestCase(TestStatefbk)
+   status2 = unittest.TestLoader().loadTestsFromTestCase(TestStatefbk)
+   return status1 and status2
 
 if __name__ == '__main__':
     unittest.main()
