From cf60075ce74cd0fb9f0a3bd6123ac9946ec365cd Mon Sep 17 00:00:00 2001 From: iabdalkader Date: Sat, 22 Dec 2018 20:37:06 +0200 Subject: [PATCH] Add micro-linalg library. --- scripts/libraries/ulinalg.py | 173 +++++++++++++ scripts/libraries/umatrix.py | 473 +++++++++++++++++++++++++++++++++++ 2 files changed, 646 insertions(+) create mode 100644 scripts/libraries/ulinalg.py create mode 100644 scripts/libraries/umatrix.py diff --git a/scripts/libraries/ulinalg.py b/scripts/libraries/ulinalg.py new file mode 100644 index 000000000..7ad79c911 --- /dev/null +++ b/scripts/libraries/ulinalg.py @@ -0,0 +1,173 @@ +''' + +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 math +import umatrix + + +def zeros(m, n, dtype=umatrix.ddtype): + return umatrix.matrix([[0 for i in range(n)] for j in range(m)], dtype=dtype) + + +def ones(m, n, dtype=umatrix.ddtype): + return zeros(m, n, dtype) + 1 + + +def eye(m, dtype=umatrix.ddtype): + Z = zeros(m, m, dtype=dtype) + for i in range(m): + Z[i, i] = 1 + return Z + +def det_inv(x): + ''' Return (det(x) and inv(x)) + + Operates on a copy of x + Using elementary row operations convert X to an upper matrix + the product of the diagonal = det(X) + Continue to convert X to the identity matrix + All the operation carried out on the original identity matrix + makes it the inverse of X + ''' + if not x.is_square: + raise ValueError('Matrix must be square') + else: + # divide each row element by [0] to give a one in the first position + # (may have to find a row to switch with if first element is 0) + x = x.copy() + inverse = eye(len(x), dtype=float) + sign = 1 + factors = [] + p = 0 + while p < len(x): + d = x[p, p] + if abs(d) < umatrix.flt_eps: + # pivot == 0 need to swap a row + # check if swap row also has a zero at the same position + np = 1 + while (p + np) < len(x) and abs(x[p + np, p]) < umatrix.flt_eps: + np += 1 + if (p + np) == len(x): + # singular + return [0, []] + # swap rows + z = x[p + np] + x[p + np, :] = x[p] + x[p, :] = z + # do identity + z = inverse[p + np] + inverse[p + np, :] = inverse[p] + inverse[p, :] = z + # change sign of det + sign = -sign + continue + factors.append(d) + # change target row + for n in range(p, len(x)): + x[p, n] = x[p, n] / d + # need to do the entire row for the inverse + for n in range(len(x)): + inverse[p, n] = inverse[p, n] / d + # eliminate position in the following rows + for i in range(p + 1, len(x)): + # multiplier is that column entry + t = x[i, p] + for j in range(p, len(x)): + x[i, j] = x[i, j] - (t * x[p, j]) + for j in range(len(x)): + inverse[i, j] = inverse[i, j] - (t * inverse[p, j]) + p = p + 1 + s = sign + for i in factors: + s = s * i # determinant + # travel through the rows eliminating upper diagonal non-zero values + for i in range(len(x) - 1): + # final row should already be all zeros + # except for the final position + for p in range(i + 1, len(x)): + # multiplier is that column entry + t = x[i, p] + for j in range(i + 1, len(x)): + x[i, j] = x[i, j] - (t * x[p, j]) + for j in range(len(x)): + inverse[i, j] = inverse[i, j] - (t * inverse[p, j]) + return (s, inverse) + + +def pinv(X): + ''' Calculates the pseudo inverse Adagger = (A'A)^-1.A' ''' + Xt = X.transpose() + d, Z = det_inv(dot(Xt, X)) + return dot(Z, Xt) + + +def dot(X, Y): + ''' Dot product ''' + if X.size(2) == Y.size(1): + Z = [] + for k in range(X.size(1)): + for j in range(Y.size(2)): + Z.append(sum([X[k, i] * Y[i, j] for i in range(Y.size(1))])) + return umatrix.matrix(Z, cstride=1, rstride=Y.size(2)) + else: + raise ValueError('shapes not aligned') + + +def cross(X, Y, axis=1): + ''' Cross product + axis=1 Numpy default + axis=0 MATLAB, Octave, SciLab default + ''' + if axis == 0: + X = X.T + Y = Y.T + if (X.n in (2, 3)) and (Y.n in (2, 3)): + if X.m == Y.m: + Z = [] + for k in range(min(X.m, Y.m)): + z = X[k, 0] * Y[k, 1] - X[k, 1] * Y[k, 0] + if (X.n == 3) and (Y.n == 3): + Z.append([X[k, 1] * Y[k, 2] - X[k, 2] * Y[k, 1], + X[k, 2] * Y[k, 0] - X[k, 0] * Y[k, 2], z]) + else: + Z.append([z]) + if axis == 0: + return umatrix.matrix(Z).T + else: + return umatrix.matrix(Z) + else: + raise ValueError('shape mismatch') + else: + raise ValueError('incompatible dimensions for cross product' + ' (must be 2 or 3)') + +def eps(x = 0): + # ref. numpy.spacing(), Octave/MATLAB eps() function + if x: + return 2**(math.floor(math.log(abs(x))/math.log(2)))*umatrix.flt_eps + else: + return umatrix.flt_eps diff --git a/scripts/libraries/umatrix.py b/scripts/libraries/umatrix.py new file mode 100644 index 000000000..e7597864e --- /dev/null +++ b/scripts/libraries/umatrix.py @@ -0,0 +1,473 @@ +''' + +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)