Matrices for the Python interpreter. An attempt to port https://github.com/jalawson/ulinalg
_B='//' _A='**' stypes=[bool,int] ddtype=int estypes=[] flt_eps=1 class matrix(object): def __init__(A,data,cstride=0,rstride=0,dtype=None): F=rstride;E=dtype;D=cstride;B=data if D!=0: if D==1:A.n=F;A.m=int(len(B)/A.n) else:A.m=D;A.n=int(len(B)/A.m) A.cstride=D;A.rstride=F;A.data=B else: A.n=1 if type(B)==int:A.m=1 else: A.m=len(B) if type(B[0])==list:A.n=len(B[0]) A.data=[B[C][G]for C in(range(A.m))for G in(range(A.n))];A.cstride=1;A.rstride=A.n if E is None:A.dtype=stypes[max([stypes.index(type(C))for C in(A.data)])] elif E in stypes:A.dtype=E else:raise TypeError('unsupported type',E) A.data=[A.dtype(C)for C in(A.data)] def __len__(A):return A.m def __eq__(A,other): B=other if A.shape==B.shape:D=all([A.data[C]==B.data[C]for C in(range(A.size()))]);return D and A.shape==B.shape else:raise ValueError('shapes not equal') def __ne__(A,other):return not __eq__(other) def __iter__(A): A.cur=0 if A.m==1:A.cnt_lim=A.n else:A.cnt_lim=A.m return A def __next__(A): if A.cur>=A.cnt_lim:raise StopIteration A.cur=A.cur+1 if A.m==1:return A.data[A.cur-1] else:return A[A.cur-1] def slice_to_offset(A,r0,r1,c0,c1):B=[A.data[C*A.rstride+D*A.cstride]for C in(range(r0,r1))for D in(range(c0,c1))];return matrix(B,cstride=1,rstride=c1-c0) def slice_indices(B,index,axis=0): C=axis;A=index if isinstance(A.start,type(None)):D=0 else:D=min(int(A.start),B.shape[C]) if isinstance(A.stop,type(None)):E=B.shape[C] else:E=min(int(A.stop),B.shape[C]) return D,E def __getitem__(D,index): A=index if type(A)==tuple: if isinstance(A[0],int):B=A[0];E=B+1 else:B,E=D.slice_indices(A[0],0) if isinstance(A[1],int):C=A[1];F=C+1 else:C,F=D.slice_indices(A[1],1) elif type(A)==list:raise NotImplementedError('Fancy indexing') else:B=A;E=B+1;C=0;F=D.n G=D.slice_to_offset(B,E,C,F) if E==B+1 and F==C+1:return G.data[0] else:return G def __setitem__(C,index,val): B=index;A=val if type(B)!=tuple:raise NotImplementedError('Need to use the slice [1,:] format.') if isinstance(B[0],int):D=B[0];H=D+1 else:D,H=C.slice_indices(B[0],0) if isinstance(B[1],int):E=B[1];I=E+1 else:E,I=C.slice_indices(B[1],1) if type(A)==matrix:A=A.data elif type(A)not in[list,tuple]:A=[A] if not all([type(F)in stypes for F in A]):raise ValueError('Non numeric entry') else: G=0 for F in range(D,H): for J in range(E,I):C.data[F*C.rstride+J*C.cstride]=C.dtype(A[G]);G=(G+1)%len(A) def __repr__(B): C=0 for D in B.data:C=max(C,len(repr(D))) A='mat([';E=0 for D in range(B.m): F=0;A=A+'[' for H in range(B.n): G=repr(B.data[E+F]);A=A+G+' '*(C-len(G)) if H<B.n-1:A=A+', ' F=F+B.cstride if D<B.m-1:A=A+'],\n ' else:A=A+']' E=E+B.rstride A=A+'])';return A def __neg__(A):B=[A.data[C]*-1 for C in(range(len(A.data)))];return matrix(B,cstride=A.cstride,rstride=A.rstride) def __do_op__(C,a,b,op): B='division by zero';A=op if A=='+':return a+b elif A=='-':return a-b elif A=='*':return a*b elif A==_A:return a**b elif A=='/': try:return a/b except ZeroDivisionError:raise ZeroDivisionError(B) elif A==_B: try:return a//b except ZeroDivisionError:raise ZeroDivisionError(B) else:raise NotImplementedError('Unknown operator ',A) def __OP__(A,a,op): G='could not be broadcast';E=op if type(a)in stypes:F=[A.__do_op__(A.data[B],a,E)for B in(range(len(A.data)))];return matrix(F,cstride=A.cstride,rstride=A.rstride) elif type(a)==list: if A.n==1 and len(a)==A.m:return A.__OP__(matrix([a]).T,E) elif len(a)==A.n:return A.__OP__(matrix([a]),E) else:raise ValueError(G) elif type(a)==matrix: if A.m==a.m and A.n==a.n:F=[A.__do_op__(A[(B,C)],a[(B,C)],E)for B in(range(A.m))for C in(range(A.n))];return matrix(F,cstride=1,rstride=A.n) elif A.m==a.m: D=A.copy() for B in range(A.n): for C in range(A.m):D[(C,B)]=A.__do_op__(D[(C,B)],a[(C,0)],E) return D elif A.n==a.n: D=A.copy() for B in range(A.m): for C in range(A.n):D[(B,C)]=A.__do_op__(D[(B,C)],a[(0,C)],E) return D else:raise ValueError(G) raise NotImplementedError('__OP__ matrix + ',type(a)) def __add__(A,a):return A.__OP__(a,'+') def __radd__(A,a):return A.__add__(a) def __sub__(A,a): if type(a)in estypes:return A.__add__(-a) raise NotImplementedError('__sub__ matrix -',type(a)) def __rsub__(A,a):A=-A;return A.__add__(a) def __mul__(A,a):return A.__OP__(a,'*') def __rmul__(A,a):return A.__mul__(a) def __truediv__(A,a):return A.__OP__(a,'/') def __rtruediv__(A,a):return A.__OP__(a,'/') def __floordiv__(A,a):return A.__OP__(a,_B) def __rfloordiv__(A,a):return A.__OP__(a,_B) def __pow__(A,a):return A.__OP__(a,_A) def __rpow__(A,a):return A.__OP__(a,_A) def copy(A):return matrix([B for B in(A.data)],cstride=A.cstride,rstride=A.rstride) def size(A,axis=0):return[A.m*A.n,A.m,A.n][axis] @property def shape(self):return self.m,self.n @shape.setter def shape(self,nshape): B=nshape;A=self if B[0]*B[1]==A.size():A.m,A.n=B;A.cstride=1;A.rstride=A.n else:raise ValueError('total size of new matrix must be unchanged') return A @property def is_square(self):return self.m==self.n def reshape(B,nshape):A=B.copy();A.shape=nshape;return A @property def T(self):return self.transpose() def transpose(A): B=matrix(A.data,cstride=A.rstride,rstride=A.cstride) if A.cstride==A.rstride:B.shape=A.n,A.m return B def reciprocal(A,n=1):return matrix([n/B for B in(A.data)],cstride=A.cstride,rstride=A.rstride) def apply(A,func,*B,**C):return matrix([func(D,*B,**C)for D in(A.data)],cstride=A.cstride,rstride=A.rstride) def matrix_isclose(x,y,rtol=1e-05,atol=flt_eps): for A in range(x.size()): try:B=[abs(x.data[A]-y.data[A])<=atol+rtol*abs(y.data[A])for A in(range(len(x.data)))] except (AttributeError,IndexError):B=[False for A in(range(len(x.data)))] return matrix(B,cstride=x.cstride,rstride=x.rstride,dtype=bool) def matrix_equal(x,y,tol=0): A=False if type(y)==matrix: if x.shape==y.shape:A=all([abs(x.data[B]-y.data[B])<=tol for B in(range(x.size()))]) return A def matrix_equiv(x,y): A=False if type(y)==matrix: if x.size()==y.size():A=all([x.data[B]==y.data[B]for B in(range(len(x.data)))]) return A def fp_eps(): A=1 while 1+A>1:A=A/2 return 2*A flt_eps=fp_eps() try:stypes.append(float);ddtype=float except:pass try:stypes.append(complex) except:pass estypes=[matrix] estypes.extend(stypes)