From cec34bee27b184769f81e473e95a49421e9eb0a4 Mon Sep 17 00:00:00 2001 From: Nguyễn Gia Phong Date: Wed, 17 Oct 2018 22:14:52 +0700 Subject: MATH1.1 --- toys/gauss.py | 83 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 83 insertions(+) create mode 100755 toys/gauss.py (limited to 'toys') diff --git a/toys/gauss.py b/toys/gauss.py new file mode 100755 index 0000000..793be06 --- /dev/null +++ b/toys/gauss.py @@ -0,0 +1,83 @@ +#!/usr/bin/env python3 +"""This module provide Gaussian and Gauss-Jordan eliminations on numpy +array of floats since they are not provided by numpy. + +Note that matrix operations are straightforward with numpy.array, namely +add (+ -), scalar ops (+ - * / **) and mul (@), so try to make use of +the library while doing your matrix homework ;-P +""" + +from itertools import islice +from math import inf, isclose +from numpy import array + + +def non0(row): + """Return index of first nonzero entry in row.""" + for index, entry in enumerate(row): + if entry: return index + else: + return inf + + +def gauss(mat): + """Return a matrix in row-echelon form, which is row-equivalent to mat, + along with shorthand method of notation of row operations performed. + """ + if mat.size == 0: return mat, [] + m, matrix, operations = len(mat), mat.__copy__(), [] + # Fuck optimization: we won't have 1000x1000 matrices for homework anyway. + for i in range(m): + while i + 1 < m: + eye, ar = min(islice(enumerate(matrix), i + 1, m), + default=inf, key=(lambda t: non0(t[1]))) + if non0(ar) >= non0(matrix[i]): break + operations.append('R{} <-> R{}'.format(i + 1, eye + 1)) + matrix[i], matrix[eye] = ar, matrix[i] + + row = matrix[i] + j = non0(row) + if j == inf: return matrix, operations + if row[j] != 1: # floating point accuracy is a total disaster tho + operations.append('({})R{} -> R{}'.format(1 / row[j], i+1, i+1)) + row /= row[j] + + for l in range(i + 1, m): + k = -matrix[l][j] + if k: # OMG I want Python 3.7 so much + operations.append('R{} + {}R{} -> R{}'.format( + l + 1, '' if k == 1 else '({})'.format(k), i + 1, l + 1)) + matrix[l] += k * row + return matrix, operations + + +def gauss_jordan(mat): + """Return a matrix in reduced row-echelon form, which is + row-equivalent to mat, along with shorthand method of notation of + row operations performed. + """ + if mat.size == 0: return mat, [] + matrix, operations = gauss(mat) + m = len(mat) + for i in range(1, m): + j = non0(matrix[i]) + if j == inf: return matrix, operations + k = -matrix[i - 1][j] + if k: + operations.append('R{} + {}R{} -> R{}'.format( + i, '' if k == 1 else '({})'.format(k), i + 1, i)) + matrix[i - 1] += k * matrix[i] + return matrix, operations + + +if __name__ == '__main__': + A = array([[1, -2, 3, 9], + [-1, 3, 0, -4], + [2, -5, 5, 17]], dtype=float) + print(*A, sep='\n', end='\n\n') + mat, ops = gauss(A) + print(*ops, sep='\n') + print(*mat, sep='\n', end='\n\n') + mat, ops = gauss_jordan(A) + print(*ops, sep='\n') + print(*mat, sep='\n', end='\n\n') -- cgit 1.4.1