"""matrix.py The discussion about matrix indexing has been interminible and for the most part pretty pointless IMO. However, it does point out one thing: the interaction between the matrix and array classes is still pretty klunky despite a fair amount of effort trying to make them interoperate. Here's one possible alternative approach to enabling matrix operations within the context of numpy arrays. I started with a few requirements: 1. Matrices should be proper subclasses of arrays in the sense that one should be able to use matrices wherever one uses arrays with no change in the result. 2. Indexing into matrices should produce row or column vectors where appropriate. 3. There should be some syntax for matrix multiplication. I know that there are other operations defined on the matrix class, but I don't consider them worth the headache of having two class hierarchies. This is in part inspired by some comments Bill Spotz, but after listening to his later comments, I don't think that this is really what he had in mind. So, please don't blame him for any deficiencies you find herein. The usual caveats apply: this is just a sketch, it's probably buggy, use at your own risk, etc. However, you should be able to save this mail as "matrix.py" and execute it to run all of the examples below. I enforce point #1 by not overriding any ndarray methods except for __array_finalize__ and __getitem__ and they are overridden in a way that shouldn't effect code unless it is specifically testing for types. With that summarily taken care of, let's move onto point #2: indexing. First let's create a matrix. Scratch that, let's create an array of matrices. One things that this approach does is let you create arrays of matrices and vectors in a natural way. The matrix below has a shape of (3,2,2), which correspons to an array of 3 2x2 matrices.
m = matrix([ [[1, 0],[0,1]], [[2, 0],[0,2]], [[3, 0],[0,3]] ])
If we index on the first axis, we are selecting a single matrix, so we would expect an ndmatrix back.
m[0] ndmatrix([[1, 0], [0, 1]])
On the other hand, if we index along the second axis, we are indexing on the matrices columns and thus expect to get an array of row matrices.
m[:,0] ndrow([[1, 0], [2, 0], [3, 0]])
Finally if we index on the last index, we get an array of column matrices.
m[...,0] ndcolumn([[1, 0], [2, 0], [3, 0]])
Indexing of row and column matrices works similarly. It's not shown here and thus will probably turn out to be buggy. For matrix multiplication, I chose to abuse call. Perhaps someday we'll be able to convinve the powers that be to free up another operator, but for the time being this is better than "*" since that forces matrix and array into separate camps and arguably more readable than most of the new syntax proposals that I've seen. Here we multiply our set of three matrices 'm', with a single 2x2 matrix 'n':
n = matrix([[0, 1], [3, 0]]) (m)(n) ndmatrix([[[0, 1], [3, 0]], <BLANKLINE> [[0, 2], [6, 0]], <BLANKLINE> [[0, 3], [9, 0]]])
Things behave similarly when multiplying a matrix with a row vector or a row vector with a column vector.
c = column([3, 5]) (m)(c) ndcolumn([[ 3, 5], [ 6, 10], [ 9, 15]]) (n)(c) ndcolumn([5, 9])
r = row([[2, 1], [3, 4]]) (r)(c) array([11, 29]) (r[0])(c) 11
Note, however that you can't (for instance) multiply column vector with a row vector:
(c)(r) Traceback (most recent call last): ... TypeError: Cannot matrix multiply columns with anything
Finally, you'd like to be able to transpose a matrix or vector. Transposing a vector swaps the type from ndrow and ndcolumn while transposing a matrix transposes the last two axes of the array it is embedded in. Note that I use ".t" for this instead of ".T" to avoid interfering with the existing ".T". (Frankly, if it were up to me, I'd just deprecate .T).
r ndrow([[2, 1], [3, 4]]) r.t ndcolumn([[2, 1], [3, 4]]) n.t ndmatrix([[0, 3], [1, 0]])
That's about it. Enjoy (or not) as you will. """ from numpy import * _ULT = 1 # Utlimate AKA last _PEN = 2 # Penultimate AKA second to last # I'm sure there's a better way to do this and this is probably buggy... def _reduce(key, ndim): if isinstance(key, tuple): m = len(key) if m == 0: return 0 elif m == 1: key = key[0] else: if Ellipsis in key: key = (None,)*(ndim-m) + key if len(key) == ndim: return isinstance(key[-1], int) * _ULT + isinstance(key[-2], int) * _PEN if len(key) == ndim-1: return isinstance(key[-1], int) * _PEN if isinstance(key, int): if ndim == 1: return _ULT if ndim == 2: return _PEN return 0 class ndmatrix(ndarray): def __array_finalize__(self, obj): if obj.ndim < 2: raise ValueError("matrices must have dimension of at least 2") def __getitem__(self, key): rmask = _reduce(key, self.ndim) value = ndarray.__getitem__(self, key) if isinstance(value, ndmatrix): if rmask == _ULT | _PEN: return value.view(ndarray) if rmask == _ULT: return value.view(ndcolumn) if rmask == _PEN: return value.view(ndrow) return value def __call__(self, other): if isinstance(other, ndcolumn): return sum(self * other[...,newaxis,:], -1).view(ndcolumn) if isinstance(other, ndmatrix): return sum(self[...,newaxis] * other[...,newaxis,:], -2).view(ndmatrix) else: raise TypeError("Can only matrix multiply matrices by matrices or vectors") def transpose_vector(self): return self.swapaxes(-1,-2) t = property(transpose_vector) class ndcolumn(ndarray): def __array_finalize__(self, obj): if obj.ndim < 1: raise ValueError("vectors must have dimension of at least 1") def __getitem__(self, key): rmask = _reduce(key, self.ndim) value = ndarray.__getitem__(self, key) if isinstance(value, ndmatrix): if rmask == _ULT: return value.view(ndarray) return value def __call__(self, other): raise TypeError("Cannot matrix multiply columns with anything") def transpose_vector(self): return self.view(ndrow) t = property(transpose_vector) class ndrow(ndarray): def __array_finalize__(self, obj): if obj.ndim < 1: raise ValueError("vectors must have dimension of at least 1") def __getitem__(self, key): rmask = _reduce(key, self.ndim) value = ndarray.__getitem__(self, key) if isinstance(value, ndmatrix): if rmask == _ULT: return value.view(ndarray) return value def __call__(self, other): if isinstance(other, ndcolumn): return sum(self * other, -1).view(ndarray) if isinstance(othe, ndmatrix): raise notimplemented else: raise TypeError("Can only matrix multiply rows with matrices or columns") def transpose_vector(self): return self.view(ndcolumn) t = property(transpose_vector) def matrix(*args, **kwargs): m = array(*args, **kwargs) return m.view(ndmatrix) def column(*args, **kwargs): v = array(*args, **kwargs) return v.view(ndcolumn) def row(*args, **kwargs): v = array(*args, **kwargs) return v.view(ndrow) if __name__ == "__main__": import doctest, matrix doctest.testmod(matrix)