add return_type keyword to functions that previously returned np.matrix

murrayrm · murrayrm · commit fbc58ecc5e7e · 2019-04-16T21:52:58.000-07:00
diff --git a/control/modelsimp.py b/control/modelsimp.py
@@ -45,9 +45,10 @@
 
 # External packages and modules
 import numpy as np
+import warnings
 from .exception import ControlSlycot
 from .lti import isdtime, isctime
-from .statesp import StateSpace, ssmatrix
+from .statesp import StateSpace
 from .statefbk import gram
 
 __all__ = ['hsvd', 'balred', 'modred', 'era', 'markov', 'minreal']
@@ -56,7 +57,7 @@
 #   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):
+def hsvd(sys, return_type=np.matrix):
     """Calculate the Hankel singular values.
 
     Parameters
@@ -66,9 +67,12 @@ def hsvd(sys):
 
     Returns
     -------
-    H : Matrix
+    H : matrix
         A list of Hankel singular values
 
+    return_type: nparray subtype, optional (default = numpy.matrix)
+        Set the ndarray subtype for the return value
+
     See Also
     --------
     gram
@@ -86,6 +90,12 @@ def hsvd(sys):
     >>> 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")
@@ -96,11 +106,12 @@ def hsvd(sys):
     w, v = np.linalg.eig(WoWc)
 
     hsv = np.sqrt(w)
-    hsv = np.array(hsv, ndmin=2)                # was np.matrix(hsv)
+    hsv = np.array(hsv, ndmin=2)        # was np.matrix(hsv)
     hsv = np.sort(hsv)
     hsv = np.fliplr(hsv)
-    # Return the Hankel singular values
-    return hsv
+
+    # Return the Hankel singular values (casting type, if needed)
+    return hsv.view(type=return_type)
 
 def modred(sys, ELIM, method='matchdc'):
     """
diff --git a/control/statefbk.py b/control/statefbk.py
@@ -42,6 +42,7 @@
 # External packages and modules
 import numpy as np
 import scipy as sp
+import warnings
 from . import statesp
 from .exception import ControlSlycot, ControlArgument, ControlDimension
 
@@ -219,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):
+def acker(A, B, poles, return_type=np.matrix):
     """Pole placement using Ackermann method
 
     Call:
@@ -238,12 +239,18 @@ def acker(A, B, poles):
         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 = np.array(A)
-    b = np.array(B)
+    a = statesp.ssmatrix(A)
+    b = statesp.ssmatrix(B)
 
     # Make sure the system is controllable
-    ct = ctrb(A, B)
+    ct = ctrb(A, B, return_type=np.ndarray)
     if np.linalg.matrix_rank(ct) != a.shape[0]:
         raise ValueError("System not reachable; pole placement invalid")
 
@@ -258,7 +265,7 @@ def acker(A, B, poles):
     K = np.linalg.solve(ct, pmat)
 
     K = K[-1][:]                # Extract the last row
-    return K
+    return K.view(type=return_type)
 
 def lqr(*args, **keywords):
     """lqr(A, B, Q, R[, N])
@@ -368,14 +375,18 @@ def lqr(*args, **keywords):
 
     return K, S, E
 
-def ctrb(A,B):
+
+def ctrb(A, B, return_type=np.matrix):
     """Controllabilty matrix
 
     Parameters
     ----------
     A, B: array_like or string
         Dynamics and input matrix of the system
 
+    return_type: nparray subtype, optional (default = numpy.matrix)
+        Set the ndarray subtype for the return value
+
     Returns
     -------
     C: matrix
@@ -386,25 +397,36 @@ def ctrb(A,B):
     >>> 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 = np.array(A)
-    bmat = np.array(B)
+    amat = statesp.ssmatrix(A)
+    bmat = statesp.ssmatrix(B)
     n = np.shape(amat)[0]
     # Construct the controllability matrix
     ctrb = bmat
     for i in range(1, n):
         ctrb = np.hstack((ctrb, np.dot(np.linalg.matrix_power(amat, i), bmat)))
-    return ctrb
 
-def obsv(A, C):
+    # Return the observability matrix in the desired type
+    return ctrb.view(type=return_type)
+
+
+def obsv(A, C, return_type=np.matrix):
     """Observability matrix
 
     Parameters
     ----------
     A, C: array_like or string
         Dynamics and output matrix of the system
 
+    return_type: nparray subtype, optional (default = numpy.matrix)
+        Set the ndarray subtype for the return value
+
     Returns
     -------
     O: matrix
@@ -414,20 +436,28 @@ def obsv(A, C):
     --------
     >>> 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 = np.array(A)
-    cmat = np.array(C)
+    amat = statesp.ssmatrix(A)
+    cmat = statesp.ssmatrix(C)
     n = np.shape(amat)[0]
 
     # Construct the controllability matrix
     obsv = cmat
     for i in range(1, n):
         obsv = np.vstack((obsv, np.dot(cmat, np.linalg.matrix_power(amat, i))))
-    return obsv
 
-def gram(sys,type):
+    # Return the observability matrix in the desired type
+    return obsv.view(type=return_type)
+
+
+def gram(sys, type, return_type=np.matrix):
     """Gramian (controllability or observability)
 
     Parameters
@@ -436,8 +466,9 @@ def gram(sys,type):
         State-space system to compute Gramian for
     type: String
         Type of desired computation.
-        `type` is either 'c' (controllability) or 'o' (observability). To compute the
-        Cholesky factors of gramians use 'cf' (controllability) or 'of' (observability)
+        `type` is either 'c' (controllability) or 'o' (observability). To
+        compute the Cholesky factors of gramians use 'cf' (controllability) or
+        'of' (observability)
 
     Returns
     -------
@@ -463,6 +494,11 @@ def gram(sys,type):
     >>> 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):
@@ -501,7 +537,7 @@ def gram(sys,type):
         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
+        return gram.view(type=return_type)
 
     elif type=='cf' or type=='of':
         #Compute cholesky factored gramian from slycot routine sb03od
@@ -524,4 +560,4 @@ def gram(sys,type):
             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
+        return gram.view(type=return_type)
diff --git a/control/tests/modelsimp_test.py b/control/tests/modelsimp_test.py
@@ -5,6 +5,7 @@
 
 import unittest
 import numpy as np
+import warnings
 from control.modelsimp import *
 from control.matlab import *
 from control.exception import slycot_check
@@ -17,10 +18,20 @@ def testHSVD(self):
         C = np.array([[6., 8.]])
         D = np.array([[9.]])
         sys = ss(A,B,C,D)
-        hsv = hsvd(sys)
+        hsv = hsvd(sys, return_type=np.ndarray)
         hsvtrue = np.array([[24.42686, 0.5731395]]) # from MATLAB
         np.testing.assert_array_almost_equal(hsv, hsvtrue)
 
+        # Make sure default type values are correct
+        self.assertTrue(isinstance(hsv, np.ndarray))
+
+        # Check that default type generates a warning
+        # TODO: remove this check with matrix type is deprecated
+        with warnings.catch_warnings(record=True) as w:
+            hsv = hsvd(sys)
+            self.assertTrue(issubclass(w[-1].category, UserWarning))
+            self.assertTrue(isinstance(hsv, np.ndarray))
+
     def testMarkov(self):
         U = np.array([[1.], [1.], [1.], [1.], [1.]])
         Y = U
diff --git a/control/tests/statefbk_test.py b/control/tests/statefbk_test.py
