Select Git revision
Jack Eckert authored
linAlg.py 22.79 KiB
# Jack Eckert
'''
This is NOT an alternative to numpy. This is purely an exploratory program that I made
in order to teach myself classes.
This program is a shoddy attempt at a custom-built vector/matrix calculator. It contains
two classes, vector and matrix that have all of the standard operations (besides one
that allows you to calculate the eigenvalues/eigenvectors of a matrix). Documentation
on this program may be quite scarce, as it was a very early project in my development
journey. The efficiency of the program also isn't very scalable, and it only really works
well for small matrices.
Comments made in docstrings will represents comments on my old code that explain major
differences there would be between programs if I were to remake this program in the modern day.
'''
# dependencies
import math
import random
# MATRIX CLASS
'''Having each entry be a seperate parameter acted as a neat way to practice the *args keyword,
but a more efficient program might implement the "numpy approach" of structuring entries of list
within lists.'''
# first parameter is the number of rows, second is the number of columns, and the remaining parameters
# are the entries of the matrix, from left to right, top to bottom
class matrix():
def __init__(self, rows, columns, *entries):
self.numRows = rows
self.numCols = columns
'''This should be a private attribute. I meant it to be used as one but I didn't know how to implement them
at the time. I will simply write "private" above all methods/attributes that were supposed to be private.'''
self.entryList = list(entries)
self.pivList = []
for entry in entries:
if not(isinstance(entry, int) or isinstance(entry, float)):
raise TypeError("Please enter numeric values")
if len(entries) != rows * columns:
raise ValueError("Incorrect number of entries")
for row in range(1, rows + 1):
for col in range(1, columns + 1):
self.__dict__[f'{row} {col}'] = entries[((row - 1) * columns + col) - 1]
'''Having the key for indexing an item as a string is awkward considering the string just
contains numerals and whitespace and you have to convert integers to strings every time.
A better approach would likely be to have the key be a number which indexes the row as a
list, and then a second key could be added to index the column. E.g. instead of something
like matrix["3 4"] it would be matrix[3][4].'''
# gets an item where the key is a string that contains the row number (starting at 1)
# followed by a space followed by the column number (starting at 1)
def __getitem__(self, x):
return self.__dict__[x]
# sets an item based on the key
def __setitem__(self, d, x):
setattr(self, d, x)
'''I was not aware of the justification specifiers for f-strings when I wrote this code. In the
modern day, I would likely use those as sometimes columns aren't aligned. This doens't affect the
actual math, though.'''
# converts a matrix to a string, in a typical form similar to how it is written on paper
def __str__(self):
i = [len(str(round(entry, 3))) for entry in self.entryList]
i = max(i)
self.str = ""
for row in range(1, self.numRows + 1):
self.str += "|"
for col in range(1, self.numCols + 1):
strAdd = (f'{self.getEntry(row, col):.03f}')
counter = 0
for char in strAdd:
if char.isdigit():
counter += int(char)
if counter == 0 and strAdd[0] == "-":
strAdd = strAdd[1:]
self.str += ((" " * (i - len(strAdd))) + f"{strAdd} ")
self.str = self.str[:-1] + "|\n"
return self.str[:-1]
# allows to check if two matrices are the same by every entry
def __eq__(self, m):
if self.entryList == m.entryList:
return True
else:
return False
# allows adding two matrices entry-by-entry
def __add__(self, m):
entries = []
for i in range(len(self.entryList)):
entries.append(self.entryList[i] + m.entryList[i])
return matrix(self.numRows, self.numCols, *entries)
# allows subtracting two matrices entry-by-entry
def __sub__(self, m):
entries = []
for i in range(len(self.entryList)):
entries.append(self.entryList[i] - m.entryList[i])
return matrix(self.numRows, self.numCols, *entries)
# either multiplies a matrix by a constant (if m is an integer or float) or by
# another matrix (if m is another matrix)
def __mul__(self, m):
if isinstance(m, int) or isinstance(m, float):
entries = [entry * m for entry in self.entryList]
return matrix(self.numRows, self.numCols, *entries)
if isinstance(m, vector):
m = m.matrix()
if not(isinstance(m, matrix)):
raise TypeError('A matrix must be multiplied by a integer, float, vector, or matrix.')
if self.numCols != m.numRows:
raise ValueError('The amount of columns in the first matrix must equal the amount of rows in the second')
if isinstance(m, matrix):
entries = []
for row in range(1, self.numRows + 1):
for col in range(1, m.numCols + 1):
eSum = 0
for i in range(1, self.numCols + 1):
eSum += (getattr(self, f'{row} {i}') * getattr(m, f'{i} {col}'))
entries.append(eSum)
return matrix(self.numRows, m.numCols, *entries)
# handles pre-multiplication by a constant
def __rmul__(self, x):
if not(isinstance(x, float) or isinstance(x, int)):
raise TypeError('Must be an integer or float')
entries = [entry * x for entry in self.entryList]
return matrix(self.numRows, self.numCols, *entries)
# divides all entries of a matrix by specified value
def __truediv__(self, x):
if not(isinstance(x, float) or isinstance(x, int)):
raise TypeError('Must be an integer or float')
entries = [entry / x for entry in self.entryList]
return matrix(self.numRows, self.numCols, *entries)
# all entries of the matrix are inverted and multiplied by a specified value
def __rtruediv__(self, x):
if not(isinstance(x, float) or isinstance(x, int)):
raise TypeError('Must be an integer or float')
entries = [x / entry for entry in self.entryList]
return matrix(self.numRows, self.numCols, *entries)
# negates all entries in a matrix
def __neg__(self):
return self * -1
# creates a copy of a matrix
def copy(self):
return matrix(self.numRows, self.numCols, *self.entryList)
'''private'''
# returns the row of a given index of the entry list
def getRow(self, i):
return (i // self.numCols) + 1
'''private'''
# returns the row of a given index of the entry list
def getCol(self, i):
if (i + 1) % self.numCols == 0:
return self.numCols
else:
return (i + 1) % self.numCols
'''private'''
# returns the index of an entry in the entry list given a row and column
def getInd(self, row, col):
return ((row - 1) * self.numCols + col) - 1
# returns a specified row in list format
def listRow(self, n):
return self.entryList[self.getInd(n, 1):self.getInd(n, self.numCols) + 1]
# returns a specified column in list format
def listCol(self, n):
if n < 1 or n > self.numCols:
raise TypeError(f"{n} is not a valid column")
putList = []
for i in range(len(self.entryList)):
if self.getCol(i) == n:
putList.append(self.entryList[i])
return putList
# returns an entry given its correspondig row and column
def getEntry(self, row, col):
return self.entryList[self.getInd(row, col)]
# calculates the determinant of a matrix
def det(self):
if self.numRows != self.numCols:
raise ValueError("Matrix must have equal number of rows and columns")
if self.numRows == 1:
return self['1 1']
if self.numRows == 2:
return self['1 1'] * self['2 2'] - self['1 2'] * self['2 1']
else:
sign = 1
entries = []
det = 0
for col in range(1, self.numCols + 1):
colC = col
entries.clear()
for i in range(1, len(self.entryList) + 1):
if colC == self.numCols:
colC = 0
elif i > self.numCols and i % self.numCols != colC:
entries.append(self.entryList[i - 1])
det += self[f'1 {col}'] * sign * matrix(self.numRows - 1, self.numCols - 1, *entries).det()
sign *= -1
return det
# calculates the transpose of a matrix
def transp(self):
entries = [0] * len(self.entryList)
for i in range(len(self.entryList)):
entries[i] = self.entryList[self.getInd(self.getCol(i), self.getRow(i))]
return matrix(self.numCols, self.numRows, *entries)
# calculates the minor of a given entry of a matrix
def minor(self, row, col):
entries = []
for i in range(len(self.entryList)):
if self.getRow(i) != row and self.getCol(i) != col:
entries.append(self.entryList[i])
return matrix(self.numRows - 1, self.numCols - 1, *entries).det()
# calculates the cofactor matrix
def cofact(self):
entries = []
for row in range(1, self.numRows + 1):
for col in range(1, self.numCols + 1):
if (row + col) % 2 == 0:
entries.append(self.minor(row, col))
else:
entries.append(-self.minor(row, col))
return matrix(self.numRows, self.numCols, *entries)
# calculates the adjunct matrix
def adj(self):
cof = self.cofact()
return cof.transp()
# calculates the inverse matrix
def inv(self):
return (self.adj() / self.det())
# swaps two given rows
def swapRow(self, row1, row2):
copy = self.copy()
for col in range(1, self.numCols + 1):
self[f'{row1} {col}'] = copy[f'{row2} {col}']
self[f'{row2} {col}'] = copy[f'{row1} {col}']
self.entryList[self.getInd(row1, col)] = copy[f'{row2} {col}']
self.entryList[self.getInd(row2, col)] = copy[f'{row1} {col}']
# multiplies a given row by a constant
def mulRow(self, row, x):
for col in range(1, self.numCols + 1):
self[f'{row} {col}'] *= x
self.entryList[self.getInd(row, col)] *= x
# given row is added to by a constant times a different row
def sumRow(self, row1, x, row2):
if row1 == row2:
self.mulRow(row1, x)
else:
for col in range(1, self.numCols + 1):
self[f'{row1} {col}'] += x * self.entryList[self.getInd(row2, col)]
self.entryList[self.getInd(row1, col)] += x * self.entryList[self.getInd(row2, col)]
# adds a vector or matrix onto the right of a given matrix, and returns that new matrix
def augment(self, vect):
if isinstance(vect, vector):
vect = vect.matrix()
if not(isinstance(vect, matrix)):
raise TypeError("Must augment with a matrix or vector")
if self.numRows != vect.numRows:
raise ValueError("Must have an equal number of rows")
entries = []
for row in range(1, self.numRows + 1):
entries.extend(self.listRow(row))
entries.extend(vect.listRow(row))
return matrix(self.numRows, self.numCols + vect.numCols, *entries)
# gives the row echelon form of a matrix
def ref(self):
copyM = self.copy()
pivRow = 1
self.pivList.clear()
for col in range(1, self.numCols + 1):
for row in range(pivRow, self.numRows + 1):
if round(copyM.getEntry(row, col), 5) != 0:
copyM.sumRow(pivRow, 1 / copyM.getEntry(row, col), row)
self.pivList.append((pivRow, col))
if row == pivRow:
for r in range(row + 1, self.numRows + 1):
copyM.sumRow(r, -copyM.getEntry(r, col), pivRow)
else:
for r in range(row, self.numRows + 1):
copyM.sumRow(r, -copyM.getEntry(r, col), pivRow)
pivRow += 1
if pivRow > copyM.numRows:
return copyM
return copyM
# gives the reduced row echelon form of a matrix
def rref(self):
copyM = self.ref()
for coord in self.pivList[::-1]:
for row in range(coord[0] - 1, 0, -1):
copyM.sumRow(row, -copyM.getEntry(row, coord[1]), coord[0])
return copyM
# gives the value of all of the variables for a specific matrix
# DOES NOT ACCOUNT FOR FREE VARIABLES
def solve(self):
copyM = self.rref()
putStr = ""
for coord in self.pivList:
putStr += f"x{coord[1]} = {round(copyM.getEntry(coord[0], copyM.numCols), 5)}, "
return putStr[:-2]
# converts matrix to a column
def vector(self):
if self.numCols == 1:
return vector(self.entryList)
else:
raise TypeError("Must have only one column")
# gives the matrix that can be postmultiplied by the "left" matrix to result in self
def factorR(self, left):
return left.inv() * self
# gives the matrix that can be premultiplied by the "right" matrix to result in self
def factorL(self, right):
return self * right.inv()
# returns a new matrix in upper triangular form
def upper(self):
outputMatrix = self.copy()
pivRow = 1
for col in range(1, self.numCols):
for row in range(pivRow + 1, self.numRows + 1):
r = pivRow + 1
while round(outputMatrix.getEntry(pivRow, col), 5) == 0 and r <= outputMatrix.numRows:
outputMatrix.swapRow(r, pivRow)
if outputMatrix.getEntry(pivRow, col) != 0:
outputMatrix.sumRow(row, -outputMatrix.getEntry(row, col) / outputMatrix.getEntry(pivRow, col), pivRow)
pivRow += 1
return outputMatrix
# returns a new matrix in lower triangular form
def lower(self):
return self.factorL(self.upper())
# calculates the trace
def trace(self):
if self.numRows != self.numCols:
raise ValueError("Must be a square matrix")
output = 0
for i in range(1, self.numRows + 1):
output += self.getEntry(i, i)
return output
# calculates the orthogonal matrix
def orthog(self):
outputMatrix = matrix(self.numRows, 1, *self.listCol(1))
for col in range(2, self.numCols + 1):
v = vector(*self.listCol(col))
for k in range(1, col):
if self.numRows != outputMatrix.listCol(k).count(0):
u = vector(*(outputMatrix.listCol(k)))
v -= u.proj(v)
outputMatrix = outputMatrix.augment(v.unit())
col1Mag = vector(*outputMatrix.listCol(1)).magnitude
for row in range(1, self.numRows + 1):
outputMatrix[f'{row} {1}'] /= col1Mag
outputMatrix.entryList[outputMatrix.getInd(row, 1)] /= col1Mag
return outputMatrix
# returns the matrix in upper hessenburg form
def upHess(self):
output = self.copy()
for col in range(1, self.numCols):
for row in range(col + 2, self.numRows + 1):
output = givRot(self.numRows, row, col + 1, output.getEntry(col + 1, col), output.getEntry(row, col)) * output
return output
# returns a tuple that represents the QR decomposition
def decomposeQR(self):
return -self.orthog(), -self.factorR(self.orthog())
#TODO
def eigenvalues(self, precision=100):
pass
# returns an identity matrix of the given size
def getIdentity(d):
entries = []
diagonal = 0
for i in range(d ** 2):
if i == diagonal:
entries.append(1)
diagonal += d + 1
else:
entries.append(0)
return matrix(d, d, *entries)
# performs a given rotation
def givRot(d, p, q, a, b):
put = getIdentity(d)
entries = put.entryList
for i in range(len(entries)):
if i == put.getInd(p, p) or i == put.getInd(q, q):
entries[i] = a / (math.sqrt(a ** 2 + b ** 2))
if i == put.getInd(p, q):
entries[i] = -b / (math.sqrt(a ** 2 + b ** 2))
if i == put.getInd(q, p):
entries[i] = b / (math.sqrt(a ** 2 + b ** 2))
return matrix(d, d, *entries)
# gives a matrix containing random elements
def randMat(rows, cols, low, high):
entries = []
for i in range(rows * cols):
entries.append(random.uniform(low, high))
return matrix(rows, cols, *entries)
# END OF MATRIX
# VECTOR CLASS
# parameters are the components of the vector
class vector():
def __init__(self, *comps):
self.dimension = 0
for comp in comps:
self.dimension += 1
self.__dict__[f'x{self.dimension}'] = comp
self.magnitude = math.sqrt(sum(((getattr(self, f'x{i}')) ** 2) for i in range(1, self.dimension + 1)))
# getter method
def __getitem__(self, x):
return self.__dict__[x]
# setter method
def __setitem__(self, d, x):
setattr(self, d, x)
# allows you to iterate through the vector
def __iter__(self):
self.a = 0
return self
def __next__(self):
if self.a < self.dimension:
self.a += 1
return self.__dict__[f'x{self.a}']
else:
raise StopIteration
# gives the dimension of the vector when the len() function is called on it
def __len__(self):
return self.dimension
# gives the magnitude of the vector when the abs() function is called on it
def __abs__(self):
return self.magnitude
# returns the vector as a string
def __str__(self):
self.str = "<"
for x in self:
self.str += f'{x}, '
if len(self.str) > 1:
return self.str[:-2] + ">"
else:
return "< >"
# allows adding two vectors component-wise
def __add__(self, v):
if len(self) != len(v):
raise ValueError("Vectors must have the same dimension")
components_sum = [getattr(self, f'x{i}') + getattr(v, f'x{i}') for i in range(1, len(self) + 1)]
return vector(*components_sum)
# allows subtracting two vectors component-wise
def __sub__(self, v):
if len(self) != len(v):
raise ValueError("Vectors must have the same dimension")
components_dif = [getattr(self, f'x{i}') - getattr(v, f'x{i}') for i in range(1, len(self) + 1)]
return vector(*components_dif)
# allows multiplying the vector by a constant
def __mul__(self, x):
scalProd = [getattr(self, f'x{i}') * x for i in range(1, len(self) + 1)]
return vector(*scalProd)
def __rmul__(self, x):
scalProd = [getattr(self, f'x{i}') * x for i in range(1, len(self) + 1)]
return vector(*scalProd)
# allows dividing the vector by a constant
def __truediv__(self, x):
return self * (1 / x)
def __rtruediv__(self, x):
scaleQuot = [x / getattr(self, f'x{i}') for i in range(1, len(self) + 1)]
return vector(*scaleQuot)
# negates all components of the vector
def __neg__(self):
components_neg = [0 - getattr(self, f'x{i}') for i in range(1, len(self) + 1)]
return vector(*components_neg)
# allows checking if two vectors are equivalent
def __eq__(self, v):
return all([getattr(self, f'x{i}') == getattr(v, f'x{i}') for i in range(1, len(self) + 1)])
# returns the vector in column matrix form
def matrix(self):
entries = [getattr(self, f'x{i}') for i in range(1, len(self) + 1)]
return matrix(len(self), 1, *entries)
# returns the dot product of two vectors
def dot(self, v):
if len(self) != len(v):
raise ValueError(f"Vectors must be the same dimension | Dim1 = {self.dimension} Dim2 = {v.dimension}")
dotProd = sum(getattr(self, f'x{i}') * getattr(v, f'x{i}') for i in range(1, len(self) + 1))
return dotProd
# returns the cross product of a variable amount of vectors based on their dimension
def cross(self, *vectors):
if len(self) < 2 or len(vectors) + 2 < 2:
raise ValueError("Vectors must both have a dimension of 2 or higher")
if len(self) != len(vectors) + 2:
raise ValueError("There must be d - 1 vectors, where d is the dimension of the vectors")
vList = [1, 1]
putList = []
vList.extend(list(self))
for vect in vectors:
vList.extend(list(vect))
vList.insert(0, 1)
crossMatrix = matrix(len(vectors) + 2, len(vectors) + 2, *vList)
sign = 1
for i in range(1, len(self) + 1):
putList.append(crossMatrix.minor(1, i) * sign)
sign *= -1
return vector(*putList)
'''the default vector should be variable dimensions based on the dimension of the instance vector.'''
# gets the angle between two vectors, the unit vector (1, 0) by default
def angle(self, v=None):
if v is None:
v = vector(1, 0)
return math.acos((self.dot(v)) / (self.magnitude * v.magnitude))
# gets the angle in degrees
def angleDeg(self, v=None):
if v is None:
v = vector(1, 0)
return math.acos((self.dot(v)) / (self.magnitude * v.magnitude)) * 180 / math.pi
# returns a vector with length one in the same direction of the given vector
def unit(self):
return self / self.magnitude
# returns a normal vector of length one
'''Should change based on the dimension of the vector. This code assumes a three-dimensional
vector is being input. This is likely because the cross-product was only implemented for 3 dimensions
until near the end of this program's development, and I forgot to change it to support other dimensions'''
def norm(self, v):
return self.cross(v).unit()
# the magnitude of the vector projected on the given vector
def comp(self, v):
return self.dot(v) / self.magnitude
# the vector that represents the vector projected on the given vector
def proj(self, v):
return self.unit() * self.comp(v)
# END OF VECTOR