''' Part of the micro-linalg project to provide a small matrix / linear algebra package for MicroPython (Python3) The MIT License (MIT) Copyright (c) 2015 Jamie Lawson Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions: The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software. THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE. ''' import sys stypes = [bool, int] ddtype = int estypes = [] flt_eps = 1 class matrix(object): def __init__(self, data, cstride=0, rstride=0, dtype=None): ''' Builds a matrix representation of 'data'. 'data' can be a list (columns) of lists (rows) [[1,2,3],[4,5,6]] or a simple list organized as determined by rstride and cstride: [1,2,3,4,5,6] cstride=1, rstride=3. Elements will be of highest type included in 'data' or 'dtype' can be used to force the type. ''' if cstride != 0: if cstride == 1: self.n = rstride self.m = int(len(data) / self.n) else: self.m = cstride self.n = int(len(data) / self.m) self.cstride = cstride self.rstride = rstride self.data = data else: # else determine shape from list passed in self.n = 1 if type(data) == int: self.m = 1 else: # it is a list self.m = len(data) # is data[0] a list if (type(data[0]) == list): self.n = len(data[0]) self.data = [data[i][j] for i in range(self.m) for j in range(self.n)] self.cstride = 1 self.rstride = self.n # ensure all elements are of the same type if dtype is None: self.dtype = stypes[max([stypes.index(type(i)) for i in self.data])] else: if dtype in stypes: self.dtype = dtype else: raise TypeError('unsupported type', dtype) self.data = [self.dtype(i) for i in self.data] def __len__(self): return self.m def __eq__(self, other): if self.shape == other.shape: res = all([self.data[i] == other.data[i] for i in range(self.size())]) return res and (self.shape == other.shape) else: raise ValueError('shapes not equal') def __ne__(self, other): return not __eq__(other) def __iter__(self): self.cur = 0 # determine proper axis if self.m == 1: self.cnt_lim = self.n else: self.cnt_lim = self.m return self def __next__(self): ''' Returns a matrix if m > 1 else the next numeric element of the vector. (Numpy returns vectors if selected via slice) ''' if self.cur >= self.cnt_lim: raise StopIteration self.cur = self.cur + 1 if self.m == 1: return self.data[self.cur - 1] else: return self[self.cur - 1] def slice_to_offset(self, r0, r1, c0, c1): # check values and limit them nd = [self.data[i * self.rstride + j * self.cstride] for i in range(r0, r1) for j in range(c0, c1)] return matrix(nd, cstride=1, rstride=(c1 - c0)) def slice_indices(self, index, axis=0): # handles the unsupported slice.indices() method in uPy. # If implemented: # midx = index.indices(self.m) # should work. if isinstance(index.start, type(None)): s0 = 0 else: s0 = min(int(index.start), self.shape[axis]) if isinstance(index.stop, type(None)): p0 = self.shape[axis] else: p0 = min(int(index.stop), self.shape[axis]) return (s0, p0) def __getitem__(self, index): if type(index) == tuple: # int and int # int and slice # slice and int # slice and slice if isinstance(index[0], int): s0 = index[0] p0 = s0 + 1 else: # row slice s0, p0 = self.slice_indices(index[0], 0) if isinstance(index[1], int): s1 = index[1] p1 = s1 + 1 else: # column slice s1, p1 = self.slice_indices(index[1], 1) elif type(index) == list: # list of indices etc raise NotImplementedError('Fancy indexing') else: # type is int? This will default to returning a row s0 = index p0 = s0 + 1 s1 = 0 p1 = self.n # resultant matrix z = self.slice_to_offset(s0, p0, s1, p1) # if it's a single entry then return that entry as int, float etc. if (p0 == s0 + 1) and (p1 == s1 + 1): return z.data[0] else: return z def __setitem__(self, index, val): if type(index) != tuple: # need to make it a slice without the slice function raise NotImplementedError('Need to use the slice [1,:] format.') # int and int => single entry gets changed # combinations of int and slice => row and columns take on elements from val if isinstance(index[0], int): s0 = index[0] p0 = s0 + 1 else: # slice s0, p0 = self.slice_indices(index[0], 0) if isinstance(index[1], int): s1 = index[1] p1 = s1 + 1 else: # slice s1, p1 = self.slice_indices(index[1], 1) if type(val) == matrix: val = val.data elif type(val) not in [list, tuple]: val = [val] if not all([type(i) in stypes for i in val]): raise ValueError('Non numeric entry') else: # assign list values wrapping as necessary to fill destination k = 0 for i in range(s0, p0): for j in range(s1, p1): self.data[i * self.rstride + j * self.cstride] = (self.dtype(val[k])) k = (k + 1) % len(val) # there is also __delitem__ # def __str__(self): def __repr__(self): # things that use __str__ will fallback to __repr__ # find max string field size for formatting l = 0 for i in self.data: l = max(l, len(repr(i))) s = 'mat([' r = 0 for i in range(self.m): c = 0 s = s + '[' for j in range(self.n): s1 = repr(self.data[r + c]) s = s + s1 + ' ' * (l - len(s1)) if (j < (self.n - 1)): s = s + ', ' c = c + self.cstride if (i < (self.m - 1)): s = s + '],\n ' else: s = s + ']' r = r + self.rstride s = s + '])' return s # Reflected operations are not yet implemented in MicroPython # __rmul__ for example will not be invoked def __neg__(self): ndat =[self.data[i] * (-1) for i in range(len(self.data))] return matrix(ndat, cstride=self.cstride, rstride=self.rstride) def __do_op__(self, a, b, op): if op == '+': return (a + b) elif op == '-': return (a - b) elif op == '*': return (a * b) elif op == '**': return (a ** b) elif op == '/': try: return (a / b) except ZeroDivisionError: raise ZeroDivisionError('division by zero') elif op == '//': try: return (a // b) except ZeroDivisionError: raise ZeroDivisionError('division by zero') else: raise NotImplementedError('Unknown operator ', op) def __OP__(self, a, op): if type(a) in stypes: # matrix - scaler elementwise operation ndat = [self.__do_op__(self.data[i], a, op) for i in range(len(self.data))] return matrix(ndat, cstride=self.cstride, rstride=self.rstride) elif (type(a) == list): # matrix - list elementwise operation # hack - convert list to matrix and resubmit then it gets handled below # if self.n = 1 try transpose otherwise broadcast error to match numpy if (self.n == 1) and (len(a) == self.m): return self.__OP__(matrix([a]).T, op) elif len(a) == self.n: return self.__OP__(matrix([a]), op) else: raise ValueError('could not be broadcast') elif (type(a) == matrix): if (self.m == a.m) and (self.n == a.n): # matrix - matrix elementwise operation # use matrix indices to handle views ndat = [self.__do_op__(self[i, j], a[i, j], op) for i in range(self.m) for j in range(self.n)] return matrix(ndat, cstride=1, rstride=self.n) # generalize the following two elif for > 2 dimensions? elif (self.m == a.m): # m==m n!=n => column-wise row operation Y = self.copy() for i in range(self.n): # this call _OP_ once for each row and __do_op__ for each element for j in range(self.m): Y[j, i] = self.__do_op__(Y[j, i], a[j, 0], op) return Y elif (self.n == a.n): # m!=m n==n => row-wise col operation Y = self.copy() for i in range(self.m): # this call _OP_ once for each col and __do_op__ for each element for j in range(self.n): Y[i, j] = self.__do_op__(Y[i, j], a[0, j], op) return Y else: raise ValueError('could not be broadcast') raise NotImplementedError('__OP__ matrix + ', type(a)) def __add__(self, a): ''' matrix - scaler elementwise addition''' return self.__OP__(a, '+') def __radd__(self, a): ''' scaler - matrix elementwise addition''' ''' commutative ''' return self.__add__(a) def __sub__(self, a): ''' matrix - scaler elementwise subtraction ''' if type(a) in estypes: return self.__add__(-a) raise NotImplementedError('__sub__ matrix -', type(a)) def __rsub__(self, a): ''' scaler - matrix elementwise subtraction ''' self = -self return self.__add__(a) def __mul__(self, a): ''' matrix scaler elementwise multiplication ''' return self.__OP__(a, '*') def __rmul__(self, a): ''' scaler * matrix elementwise multiplication commutative ''' return self.__mul__(a) def __truediv__(self, a): ''' matrix / scaler elementwise division ''' return self.__OP__(a, '/') def __rtruediv__(self, a): ''' scaler / matrix elementwise division ''' return self.__OP__(a, '/') def __floordiv__(self, a): ''' matrix // scaler elementwise integer division ''' return self.__OP__(a, '//') def __rfloordiv__(self, a): ''' scaler // matrix elementwise integer division ''' return self.__OP__(a, '//') def __pow__(self, a): ''' matrix ** scaler elementwise power ''' return self.__OP__(a, '**') def __rpow__(self, a): ''' scaler ** matrix elementwise power ''' return self.__OP__(a, '**') def copy(self): """ Return a copy of matrix, not just a view """ return matrix([i for i in self.data], cstride=self.cstride, rstride=self.rstride) def size(self, axis=0): """ 0 entries 1 rows 2 columns """ return [self.m * self.n, self.m, self.n][axis] @property def shape(self): return (self.m, self.n) @shape.setter def shape(self, nshape): """ check for proper length """ if (nshape[0] * nshape[1]) == self.size(): self.m, self.n = nshape self.cstride = 1 self.rstride = self.n else: raise ValueError('total size of new matrix must be unchanged') return self @property def is_square(self): return self.m == self.n def reshape(self, nshape): """ check for proper length """ X = self.copy() X.shape = nshape return X @property def T(self): return self.transpose() def transpose(self): """ Return a view """ X = matrix(self.data, cstride=self.rstride, rstride=self.cstride) if self.cstride == self.rstride: # handle column vector X.shape = (self.n, self.m) return X def reciprocal(self, n=1): return matrix([n / i for i in self.data], cstride=self.cstride, rstride=self.rstride) def apply(self, func, *args, **kwargs): """ call a scalar function on each element, returns a new matrix passes *args and **kwargs to func unmodified note: this is not useful for matrix-matrix operations e.g. y = x.apply(math.sin) y = x.apply(lambda a,b: a>b, 5) # equivalent to y = x > 5 y = x.apply(operators.gt, 5) # equivalent to y = x > 5 (not in micropython) """ return matrix([func(i, *args, **kwargs) for i in self.data], cstride=self.cstride, rstride=self.rstride) def matrix_isclose(x, y, rtol=1E-05, atol=flt_eps): ''' Returns a matrix indicating equal elements within tol''' for i in range(x.size()): try: data = [abs(x.data[i] - y.data[i]) <= atol+rtol*abs(y.data[i]) for i in range(len(x.data))] except (AttributeError, IndexError): data = [False for i in range(len(x.data))] return matrix(data, cstride=x.cstride, rstride=x.rstride, dtype=bool) def matrix_equal(x, y, tol=0): ''' Matrix equality test with tolerance same shape''' res = False if type(y) == matrix: if x.shape == y.shape: res = all([abs(x.data[i] - y.data[i]) <= tol for i in range(x.size())]) return res def matrix_equiv(x, y): ''' Returns a boolean indicating if X and Y share the same data and are broadcastable''' res = False if type(y) == matrix: if x.size() == y.size(): res = all([x.data[i] == y.data[i] for i in range(len(x.data))]) return res def fp_eps(): ''' Determine floating point resolution ''' e = 1 while 1 + e > 1: e = e / 2 return 2 * e flt_eps = fp_eps() try: if sys.implementation.name == 'micropython' and sys.platform == 'linux': # force this as there seems to be some interaction with # some operations done using the C library with a smaller epsilon (doubles) flt_eps = 1.19E-7 # single precision IEEE 2**-23 double 2.22E-16 == 2**-52 except: pass # Determine supported types try: stypes.append(float) ddtype = float except: pass try: stypes.append(complex) except: pass # extended types estypes = [matrix] estypes.extend(stypes)