python 深度模仿 matlab 矩阵语法

直接上代码。MyMat (298行开始)类模仿matlab矩阵索引与赋值。从1开始计数，倒数从0开始。可以用单指标索引：1为第一个元素，2 为第一列第二行元素……以后还会加以改进。网友们也可以提出一些建议。（代码以pypi为准）

If you like matlab as well as python, then this is your choice. Class MyMat imitates the reference and assignment of matlab matrices. The index is form 1, instead of 0, and the last index is 0, instead of -1. I will add more functions or classes in the future to improve the modules. You can give me any advice under the article.

You can download the codes from pypi. https://pypi.python.org/pypi/mymat

截止10月1日（接近一个月），至少已有九千人下载代码。

</pre><p></p><p><pre name="code" class="python">import numpy as np
import numpy.linalg as LA

'''notation (abbreviation):
num: complex number
ind: index (int, slice or list, or tuple of index sometimes)
sglind: single index (int, slice or list)
vec: vector
colvec: column vector
slc: slice
ran: range
lst: list
mat: matrix
r(m): row
c(n): column
'''


# python for my matrix type
# special matrices

def O(m, n=None):
    # 0 matrix
    if n==None:
        n = m
    return MyMat(np.zeros((m, n)))


def One(m, n=None):
    # 1 matrix
    if n==None:
        n = m
    return MyMat(np.ones((m, n)))


def I(m, n=None):
    # identity matrix
    if n==None:
        n = m
    return MyMat(np.eye(m, n))


def E(i,j,m,n=None):
    if n==None: n = m
    A = O(m,n)
    A[i,j] = 1
    return A


def R(m, n=None, randfun='rand'):
    # random matrix
    if n==None:
        n = m
    return MyMat(getattr(np.random, randfun)(m, n))


def H(m, n=None):
    # Hilbert matrix
    if n==None:
        n = m
    return MyMat([[1 / (k + l - 1) for l in range(1, n+1)] for k in range(1, m+1)])


def TestMat(m, n=None):
    if n==None:
        n = m
    return MyMat([[l + (k-1)*n for l in range(1, n+1)] for k in range(1, m+1)])


# maps
def mat_map(f,*As):
    # matrices (np.matrix) in As have the same size
    rn, cn = As[0].shape    # size (shape) of matrix
    return np.mat([[f(*(A[k, l] for A in As)) for l in range(cn)] for k in range(rn)])
    

def mymat_map(f,*As):
    # matrices (MyMat) in As have the same size
    rn, cn = As[0].shape    # size (shape) of matrix
    return MyMat([[f(*(A[k, l] for A in As)) for l in range(1, cn+1)] for k in range(1, rn+1)])

# relations
# will be used in test
def equal(A, *Ms):
    if myMat(A).isempty():
        for M in Ms:
            if not MyMat(M).isempty():
                return False
        return True
    for M in Ms:
        M = np.mat(M)
        if A.shape != M.shape or (A != M).any():
            return False
    return True


def equal_tol(A, *Ms, tol=0.001):
    for k, M in enumerate(Ms):
        M = np.mat(M)
        for N in Ms[k+1:]:
            N = np.mat(N)
            if A.shape != M.shape or LA.norm(A - M)>=tol:
                return False
    return True

''' # Index class
# it may be used in the future
def indadd(ind, x):
    if isinstance(ind, int):
        return ind + x
    elif isinstance(ind, slice):
        return slice(ind.start + x, ind.stop + x, ind.step)
    else:
        return [a + x for a in ind]

def indinv(ind):
    if isinstance(ind, int):
        return ind
    elif isinstance(ind, slice):
        return slice(ind.stop, ind.start, -ind.step)
    else:
        ind.reverse()
        return ind

def indneg(ind):
    if isinstance(ind, int):
        return -ind
    elif isinstance(ind, slice):
        return slice(-ind.stop, -ind.start, -ind.step)
    else:
        return [-a for a in ind]
    
class Index:
    # value: int, slice, list
    def __init__(self, value=0):
        self.value=value

    def __add__(self, other):
        if isinstance(other, int):
            return Index(indadd(self.value, other))

    def __neg__(self):
        return Index(indneg(self.value))
    
    def __invert__(self):
        return Index(indinv(self.value))

    def __str__(self):
        return ','.join(map(str, ind2iter(self.value)))
'''

# MyMat Class

def isreal(x):
    return isinstance(x, (int, float))

def isnum(x):
    return np.isscalar(x)

def isscalar(x):
    return isnum(x) or isinstance(x, MyMat) and x.isscalar()

def issgl(x):
    # single index used in python
    return isinstance(x, (int, slice, list))

def isemp(x):
    # empty index
    return isinstance(x, list) and x==[] or isinstance(x, tuple) and (x[0]==[] or x[1]==[])

def iscrd(x):
    # coordinate type index
    # isinstance(x, tuple) and len(x) == 2
    return isinstance(x, tuple) and isinstance(x[0], list) and isinstance(x[1], list)

COLON = slice(None)

# transform of the indices
import itertools

def times(ind):
    '''
>>> times(([1,2],[3,4,6]))
([1, 1, 1, 2, 2, 2], [3, 4, 6, 3, 4, 6])
'''
    if iscrd(ind):
        lst = itertools.product(ind[0],ind[1])
        b=([], [])
        for a in lst:
            b[0].append(a[0])
            b[1].append(a[1])
        return b
    else:
        return ind

def ind2iter(ind, N):
    # slice to range, int and list to list
    if isinstance(ind, list):
        return ind
    elif isinstance(ind, slice): # slc is a slice
        a = ind.start if ind.start else 0
        b = ind.stop if ind.stop else N # b <= N
        return range(a, b, c) if ind.step else range(a,b)
    else:
        return [ind]

def ind2list(ind, N):
    # slice to range, int and list to list
    if isinstance(ind, list):
        return ind
    elif isinstance(ind, slice): # slc is a slice
        a = ind.start if ind.start else 0
        b = ind.stop if ind.stop else N # b <= N
        return list(range(a, b, c)) if ind.step else list(range(a,b))
    else:
        return [ind]
    
def ind2ind(ind):
    '''index of matlab type to index of python type
This is the main gimmick to define MyMat
also see slice2slice
'''
    if isinstance(ind, tuple):
        return tuple(ind2ind(a) for a in ind)
    if isinstance(ind, list):
        return [a - 1 for a in ind]
    elif isinstance(ind, slice):
        return slice2slice(ind)       # see its definition
    else:
        return ind - 1

def slice2slice(slc):
    '''slice of matlab type to slice of python type
for example: M[1:3] := M[0:3] 
>>> slice2slice(slice(1, 6, 2))
slice(0, 6, 2)
>>> slice2slice(slice(-3, 0))
slice(-4, None, None)
'''
    return slice(slc.start-1 if slc.start else None, None if slc.stop == 0 else slc.stop, slc.step)


def int2pair(ind, r, c=1):
    '''
>>> int2pair(ind2ind([3,4,7,10]),5,5)
([2, 3, 1, 4], [0, 0, 1, 1])
'''
    if isinstance(ind, int):
        return np.mod(ind, r), np.floor_divide(ind, r)
    else:
        ind=ind2iter(ind, r*c)
        return [np.mod(a, r) for a in ind], [np.floor_divide(a, r) for a in ind]

def pair2int(ind, r, c=None):
    '''inverse of int2pair
>>> pair2int(int2pair([4,7,12,22],5,5),5,5)
[4, 7, 12, 22]
'''
    if isinstance(ind[0], int) and isinstance(ind[1], int):
        return ind[0] + r * ind[1]
    else:
        if c==None:
            c=max(ind[1])
        return [a + r * b for a, b in zip(ind2iter(ind[0], r), ind2iter(ind[1], c))]

"""
# int2pair(ind2ind) == ind2ind(int2pair2)
def int2pair2(ind, r, c=1):
    '''single index (int) => pair index (tuple)
example:
>>> int2pair(7, 3)
(1, 3)
    if isinstance(ind, int):
        return (r, np.floor_divide(ind, r)) if np.mod(ind, r)==0 else (np.mod(ind, r), np.floor_divide(ind, r)+1)
    else:  # extend to lists or slices
        ind=ind2iter(ind, r * c)
        p=([],[])
        for a in ind:
            if np.mod(a, r)==0:
                p[0].append(r); p[1].append(np.floor_divide(a, r))
            else:
                p[0].append(np.mod(a, r)); p[1].append(np.floor_divide(a, r)+1)
        return p
"""

def maxind(ind):
    # if issgl(ind):
    if isinstance(ind, int):
        return ind
    elif isinstance(ind, list):
        return max(ind)
    elif isinstance(ind, slice):
        return ind.stop if ind.stop else 0
    else:
        return maxind(ind[0]), maxind(ind[1])

# definition of MyMat Class
class MyMat(np.matrix):
    # python for imitating matrices (2D) in matlab
    # matrix on complex numbers
    def __new__(cls, data, *args, **kw):
        if isinstance(data, list):
            data=np.array(data)
        elif isinstance(data, str):
            data = np.array(np.matrix(data))
        elif hasattr(data, 'tolist'):
            data = np.array(data.tolist())
        return super(MyMat, cls).__new__(cls, data, *args, **kw)

    def __mul__(self, other):
        if isinstance(other, MyMat):
            if other.isscalar():
                return super(MyMat, self).__mul__(other.toscalar())
            else:
                return super(MyMat, self).__mul__(other)
        else:
            if self.isscalar():
                return MyMat([self[1,1] * other])
            else:
                return super(MyMat, self).__mul__(other)
            

    def __rmul__(self, other):
        return self * other   

    def __imul__(self, other):
        self[:,:] = self * other
        return self

    def __pow__(self, other):
        # A**B := A.*B
        return MyMat(np.multiply(self, other))

    def __rpow__(self, other):
        # B**A := B.*A
        return MyMat(np.multiply(other, self))

    def __ipow__(self, other):
        self[:,:] = self ** other
        return self

    def __floordiv__(self, other):
        # A//B := A./B
        return MyMat(np.divide(self, other))

    def __rfloordiv__(self, other):
        # B//A := B./A
        return MyMat(np.divide(other, self))

    def __ifloordiv__(self, other):
        self[:,:] = self // other
        return self

    def __truediv__(self, other):
        # A/B := A/B (A*inv(B))
        if isinstance(other, MyMat) and other.isscalar():
            return MyMat(np.divide(self, other.toscalar()))
        elif isnum(other):
            return self * (1 / other)  # = 1/other * self       
        return self * other.I

    def __rtruediv__(self, other):
        # r/A := r/A (r * inv(A))
        return other * self.I

    def __itruediv__(self, other):
        self[:,:] = self / other
        return self
    
    def __xor__(self, num):
        # A^n := A^n, n is int
        if isinstance(num, MyMat) and num.isscalar():
            return super(MyMat, self).__pow__(num.toscalar())
        return super(MyMat, self).__pow__(num)

    def __rxor__(self, num):
        pass

    def __ixor__(self, num):
        self[:,:] = self ^ other
        return self

    def __lshift__(self, other):
        # A<<B := A.^B
        return MyMat(np.power(self, other))

    def __rlshift__(self, other):
        # r<<A := r.^A
        return MyMat(np.power(other, self))

    def __ilshift__(self, other):
        self[:,:] = self << other
        return self

    @property
    def I(self):
        return self.getI()

    def __getitem__(self, ind):
        r, c = self.shape
        if iscrd(ind):
            return self.get((ind[0], COLON)).get((COLON, ind[1]))
        else:
            return self.get(ind)

    def __setitem__(self, ind, val):
        r, c = self.shape
        if issgl(ind):
            self.set(ind, val)
        else:
            if iscrd(ind):
                if isinstance(val, MyMat):
                    if val.isscalar():
                        self.set(times(ind), val.toscalar())
                    else:
                        self.set(times(ind), val.tovec(dim=2).tolist())
                elif isnum(val):
                    self.set(times(ind), val)
                else: # isinstance(val, list)
                    self.set(times(ind), val)
            else:
                self.set(ind, val)

    def get(self, ind):
        r, c = self.shape
        if issgl(ind):
            return super(MyMat, self).__getitem__(int2pair(ind2ind(ind), r, c))
        return super(MyMat, self).__getitem__(ind2ind(ind))

    def set(self, ind, val):
        r, c = self.shape
        if issgl(ind):   # single index
            mx=maxind(ind)
            if mx > r*c:
                mx1, mx2=int2pair(mx-1, r, c)
                self.just(max(r, mx1+1), max(c, mx2+1))
            if isnum(val):
                super(MyMat, self).__setitem__(int2pair(ind2ind(ind), r, c), val)
            elif isinstance(val, list):
                super(MyMat, self).__setitem__(int2pair(ind2ind(ind), r, c), val)
            elif isinstance(val, MyMat):
                if val.isscalar():
                    super(MyMat, self).__setitem__(int2pair(ind2ind(ind), r, c), val[1, 1])
                elif val.isvec():
                    super(MyMat, self).__setitem__(int2pair(ind2ind(ind), r, c), val.tolist())
                elif val.iscolvec():
                    super(MyMat, self).__setitem__(int2pair(ind2ind(ind), r, c), val.T.tolist())
            elif isinstance(val, np.matrix):
                self.set(ind, MyMat(val))
        else:      # double index
            mx1, mx2 = maxind(ind)
            if mx1 > r or mx2 > c:
                self.just(max(r, mx1), max(c, mx2))
            if isinstance(ind[0], slice) and isinstance(ind[0], slice):
                super(MyMat, self).__setitem__(ind, val)
            else:
                super(MyMat, self).__setitem__(ind2ind(ind), val)

    ''' def set_(self, ind, val):
        # extension of setitem
        r, c = self.shape
        if issgl(ind):
            mx=maxind(ind)
            if mx > r*c:
                mx1, mx2=int2pair(mx-1, r, c)
                self.just(max(r, mx1+1), max(c, mx2+1))
        else:
            mx1, mx2 = maxind(ind)
            if mx1 > r or mx2 > c:
                self.just(max(r, mx1), max(c, mx2))
        self.set(ind, val)
        return self
'''

    def update(self, other):
        self.resize(other.shape, refcheck=False)
        super(MyMat, self).__setitem__((COLON, COLON), other)

    def __call__(self, *ind):
        if len(ind)==1:
            return self[ind[0]]
        return self[ind]

    def __str__(self):
        return '['+self.join(ch2=';\n ')+']'

    def join(self, ch1=', ', ch2='; '):
        r, c = self.shape
        return ch2.join(ch1.join(str(super(MyMat, self).__getitem__((k, l))) for l in range(c)) for k in range(r))

    def cat(self, *others, dim=1):
        if self.isempty():
            return np.concatenate(others, axis=dim-1)
        elif not others:
            return self
        else:
            return np.concatenate((self,)+others, axis=dim-1)

    def __or__(self, other):
        return np.concatenate((self, other), 1)

    def __ior__(self, other):
        self.update(self | other)
        return self

    def __and__(self, other):
        return np.concatenate((self, other))

    def __iand__(self, other):
        self.update(self & other)
        return self

    def wequal(self, *others):
        # weakly equal
        if self.isempty():
            for other in others:
                if not MyMat(other).isempty():
                    return False
            return True
        for other in others:
            M = MyMat(other)
            if self.shape != M.shape or (self != M).any():
                return False
        return True

    def equal(self, *others):
        # others are MyMat
        if self.isempty():
            for other in others:
                if not other.isempty():
                    return False
            return True
        for other in others:
            if self.shape != other.shape or (self != M).any():
                return False
        return True
    
    def toscalar(self):
        return super(MyMat, self).__getitem__((0,0))
    
    def tolist(self):
        if self.isvec():
            return super(MyMat, self).tolist()[0]
        else:
            return super(MyMat, self).tolist()

    def toarray(self):
        if self.isvec():
            return self.A1
        else:
            return self.A

    def tomatrix(self):
        return np.matrix(self.tolist())

    def tomat(self):
        # == tomatrix
        return np.mat(self.tolist())

    def tovec(self, dim=1):
        if dim==1:
            return MyMat([[self[k, l]] for l in range(1, self.shape[1]+1) for k in range(1, self.shape[0]+1)])
        elif dim==2:
            lst=[]
            for k in range(1, self.shape[0]+1):
                lst.extend(self[k,:].tolist()) 
            return MyMat(lst)

    def just(self, m, n=None):
        if n == None: n = m
        r, c = self.shape
        if m>0 and n>0:
            if m<=r:
                temp = super(MyMat, self).__getitem__((slice(m), slice(r)))
            else:
                temp = self | O(m - r, c)
            if n<=c:
                temp = super(MyMat, temp).__getitem__((slice(c), slice(n)))
            else:
                temp = temp & O(m, n - c)
        elif m==0 or n==0:
            temp = Empty
        else:
            temp = O(abs(m), abs(n))
        self.update(temp)
        # return self

    def delete(self, ind1=[], ind2=[]):
        r, c = self.shape
        if not isemp(ind1):
            ind1=ind2list(ind2ind(ind1), r)
            if ind1 == list(range(r)):
                return Empty
            self=np.concatenate(tuple(super(MyMat, self).__getitem__((slice(a+1,b), COLON)) for a, b in zip([-1]+ind1, ind1+[r]) if b-a!=1))
        if not isemp(ind2):
            ind2=ind2list(ind2ind(ind2), c)
            if ind2 == list(range(c)):
                return Empty
            return np.concatenate(tuple(super(MyMat, self).__getitem__((COLON, slice(a+1,b))) for a, b in zip([-1]+ind2, ind2+[c]) if b-a!=1), 1)
        else:
            return self

    '''def delete(self, ind1, ind2):
        ......
        if not isemtpy(ind1):
            self=super(PyMat, self).__getitem__((slice(ind1[0]), COLON)).cat(
                *(tuple(super(PyMat, self).__getitem__((slice(a+1,b), COLON)) for a, b in zip(ind1[:-1], ind1[1:]) if b-a!=1)
                +(super(PyMat, self).__getitem__((slice(ind1[-1]+1,None), COLON)),)))
            if self.isempty():
                return Empty
        if not isemtpy(ind2):
            return super(PyMat, self).__getitem__((COLON, slice(ind2[0]))).cat(
                *(tuple(super(PyMat, self).__getitem__((COLON, slice(a+1,b))) for a, b in zip(ind2[:-1], ind2[1:]) if b-a!=1)
                +(super(PyMat, self).__getitem__((COLON, slice(ind2[-1]+1,None))),)), dim=2)
        else:
            return self'''
            
    def norm(self, p=2):
        if self.isvec():
            return LA.norm(self.toarray(), p)
        elif self.iscolvec():
            return LA.norm(self.T.toarray(), p)
        else:
            return LA.norm(self, p)

    def repmat(self, m=1):
        return np.tile(self, m-1)
        
    def apply(self, fun):
        return mymat_map(fun, self)

    def isempty(self):
        return self.shape[1]==0 or self.shape[0]==0

    def isscalar(self):
        return self.shape == (1,1)

    def isvec(self):
        return self.shape[0] == 1

    def iscolvec(self):
        return self.shape[1] == 1

    @property
    def size(self):
        return self.shape

    def poly(self, p):
        # self is a square matrix
        r, c = self.shape  #  r == c
        if len(p)==1:
            return p[0]*I(r)
        return p[0]*I(r)+self.poly(p[1:])*self

    def expm(self, N=10):
        p=[1/np.prod(range(1, n+1)) for n in range(N)]
        return self.poly(p)

# matrix instance
MMat = myMat = MyMat # alias
Empty = MyMat([])    # empty set
im = MyMat([[0, 1], [-1, 0]])

def c(a=0, b=0):
    # linear represenation of complex number
    return MyMat([[a, b], [-b, a]])

def q(a=0, b=0, c=0, d=0):
    # linear represenation of quaternions
    return MyMat([[a, b, c, d], [-b, a, -d, c], [-c, d, a, -b], [-d, -c, b, a]])