Coverage increased (+0.3%) to 81.228% when pulling a5126f8 on murrayrm:np.matrix into 8a11cd3 on python-control:master.

One overall comment: If control functions, methods, and attributes now return numpy arrays instead of numpy matrices, this change could break lots of code in the wild. If a user gets these matrix outputs and then operates on them with operations that assume a matrix type (e.g. overloaded *), then their code will now break or return incorrect results. It would probably be best to deprecate this in the control package so that users see warnings that output types will change in the next release. A 6 month deprecation time might be the minimum for this to register for average users.

@moorepants Good point. Most of these changes are internal and so shouldn't affect any user outputs. But I'll go through and check to see whether the output of any function actually changed pre/post PR. The only think that might have changed are some of the functions that return matrices (eg, Gramians).

Note also that the StateSpaceMatrix class does define the "standard" overloaded operators (*,^, [,]) so user code operating on state space matrices should function normally. But I completely agree this is something we want to be careful about.

I checked through the code. Most functions that return matrices were already returning ndarrays (typically via a call to scipy.linalg.*). The one set that I found that could have compatibility problems is:

• ctrb and obsv used to return matrix objects and now return ndarray objects.

Thoughts on the best approach to handling this? We could:

• Keep the return type as matrix for now, with a deprecation warning. Change to ndarray in a future release.

• Return a StateSpaceMatrix object, which would be compatible for most (all?) operations and allow matrix-like manipulation of the observability or controllability gramians.

• Return an ndarray object and let people rewrite their code as they discover the change.

We could also allow a subtype keyword for the function and return an ndarray of the listed subtype (with a default to either matrix, StateSpaceMatrix, or ndarray, depending on how we resolve the issue).

Updated the code so that all functions that used to return an ndarray of subtype numpy.matrix return that type by default, with a user warning that this will be deprecated in the future. This can be overridden using the return_type keyword.

Once we are OK with the code changing (perhaps in the next major or minor release), we can change the code to return numpy.ndarray by default (user will still be able to override using return_type keyword).

Merged changes from current version of master (including implementing updated ctrb and obsv computations in PR #300).

@billtubbs FYI. Perhaps worth waiting until this PR is merged before making any more changes to ctrb and obsv. I took your recent update and applied that same change here (but without using the matrix class).

StateSpaceMatrix should probably have an __imul__ method.
Using __new__ here seems unusual to me; is it obvious why it's necessary?

Putting return_type=np.matrix arguably makes a long-term commitment to this interface (i.e., the "return_type" argument will be expected in the future); an alternative might be a config option for return types of the relevant functions. It also suggests more generality than is on offer; that is, a naive user might expect return_type=StateSpaceMatrix to work (it might, though I assume it's not what's intended) Perhaps a flag instead, "return_npmatrix=True"?

The code for turning Matlab-style matrix strings into nd.arrays is in numpy.matrixlib.defmatrix._convert_from_string -- it's not all that long or complicated, one could copy it if necessary.

Less importantly, if we didn't have to support Python 2.7 this refactoring could have used @ to be a bit tidier. Pity there's no tidy "matrix-power" operator.

Thanks for the comments @roryyorke. Responses (no actions yet):

• The code for __new__ is based on the numpy.matrix.__new__ function. I'm pretty sure we need to use __new__ instead of __init__ because we are creating the object (via a call to np.ndarray.__new__) with the desired storage type, shape, etc. If we just call __init__ then that information will have already been created by numpy.ndarray (since __init__ gets called after __new__).

• I'm pretty sure that return_type=StateSpaceMatrix will actually work. You can also use return_type=numpy.ndarray, if you wanted to make sure you got a NumPy array and not a StateSpace object. But I agree that we may want the simpler behavior of just turning off the matrix class. Need to think on that one.

• I hadn't thought about what other operators to implement. In addition to __imul__ there are things like __iadd__ and __isub__ as well. Implementing all of them is tedious, but not implementing them means that you get the default numpy.ndarray behavior, which could be confusing. Another alternative would be to define all of those to raise a NotImplementedError exception, so at least you couldn't "accidentally" use it.

np.matrix inherits most numeric dunder methods from np.ndarray, except mul and pow (see below). Raising NotImplementedError for unimplemented mul and pow methods is an easy, safe option.

prefixes = ['','r','i']                                                                                       
bases = ['add','sub','mul','matmul','truediv','floordiv','mod','divmod','pow','lshift','rshift','and','or','xor']
makemeth = lambda p, b: '__' + p + b + '__'
# there's no such thing as __idivmod__; otherwise all exist
methods = [makemeth(p,b) for p,b in itertools.product(prefixes, bases) if hasattr(np.ndarray, makemeth(p,b))]
matonly = [m for m in methods if getattr(np.matrix, m) is not getattr(np.ndarray, m)]                          
print(matonly)
# output: ['__mul__', '__pow__', '__rmul__', '__rpow__', '__imul__', '__ipow__']

Committed a new set of changes to address the issues raised by @roryyorke. Most recent updates:

• The use of numpy.matrix as a return type is now set as a configuration option. The default is for functions in statesp.py and (hsvd in modelsimp.py) to return numpy.matrix objects (for backward compatibility), but the call use_numpy_matrix(False) will change the return type to StateSpaceMatrix. There is also a function set_ss_return_type that can set the return type to any subclass of numpy.ndarray (including ndarray itself).

• New unit test code to test the functionality of use_numpy_matrix(False) and set_ss_return_type.

• Added code to raise theNotImplementedError exception for dunder methods ipow and imul, so that they don't fall back to the (incompatible) numpy.ndarray methods.

In a future release, we can change the defaults in config.py to return either StateSpaceMatrix objects or numpy.ndarray objects and the dependence on numpy.matrix will be removed.

I decided not to copy the string processing code from numpy.matrixintopython-control though we may eventually want to put that into the matlab submodule. I think the main module should follow whatever Numpy does in terms of allowing string initialization (my guess is that it will go away).

In looking through the code just now, I see what we have some inconsistencies in how we return matrix objects. In some place we are returning matrix objects (now convertible to StateSpaceMatrix via changes in this PR) but in other place we return ndarray objects. For example, place and lqr return ndarray objects, but ctrb, obsv, and acker return matrix-like objects (either numpy.matrix or StateSpaceMatrix, depending on settings).

On the one hand, it seems useful to have functions that return a gain matrix sending back something that can be used like the numpy.matrix class, so that you can say (for example)

K = lqr(A, B, Q, R)
Acl = A - B * K

(note that this won't actually work in the current code, tthough if you used acker instead of lqr then it would!).

On the other hand, at some point we should probably just bite the bullet and use ndarray everywhere, so that we are compatible with Numpy. If we did that, we could just get rid of the new StateSpaceMatrix class in this PR completely.

Any thoughts?

Will this work?

K = lqr(A, B, Q, R)
Acl = A - B @ K

Using @ should be what we move towards.

Yes, that will work (though not in python 2.7, of course).

Following the discussion in issue #233, this PR is being abandoned in favor of PR #314, which is aligned with the approach in issue #233 (and is much more straightforward).