88 lines
3.1 KiB
Python
Executable File
88 lines
3.1 KiB
Python
Executable File
#!/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
|
|
from numpy import array
|
|
|
|
|
|
def leadidx(row):
|
|
"""Return index of leading entry in row."""
|
|
for index, entry in enumerate(row):
|
|
if entry: return index
|
|
else:
|
|
return inf
|
|
|
|
|
|
def gauss(mat, log=False):
|
|
"""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),
|
|
key=(lambda t: leadidx(t[1])))
|
|
if leadidx(ar) >= leadidx(matrix[i]): break
|
|
operations.append('R{} <-> R{}'.format(i + 1, eye + 1))
|
|
matrix[i], matrix[eye] = ar, matrix[i]
|
|
if log: print(operations[-1], *matrix, sep='\n')
|
|
|
|
row = matrix[i]
|
|
j = leadidx(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]
|
|
if log: print(operations[-1], *matrix, sep='\n')
|
|
|
|
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
|
|
if log: print(operations[-1], *matrix, sep='\n')
|
|
return matrix, operations
|
|
|
|
|
|
def gauss_jordan(mat, log=False):
|
|
"""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, log)
|
|
m = len(mat)
|
|
for i in range(1, m):
|
|
j = leadidx(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]
|
|
if log: print(operations[-1], *matrix, sep='\n')
|
|
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')
|