feat: refactor main and richardson to be modular

This commit is contained in:
Gromiusz 2024-10-19 23:24:26 +02:00
parent 6124e2db6a
commit 5e68bddc64
2 changed files with 201 additions and 0 deletions

101
code/main_abstract.py Normal file
View File

@ -0,0 +1,101 @@
import unittest
import numpy as np # For testing ONLY
from richardson_abstract import modified_richardson
class TestModifiedRichardson(unittest.TestCase):
def setUp(self):
self.A_2x2 = np.random.rand(2, 2).tolist()
self.b_2x2 = np.random.rand(2).tolist()
self.x0_2x2 = np.random.rand(2).tolist()
self.alpha_2x2 = 0.1
self.A_3x3 = np.random.rand(3, 3).tolist()
self.b_3x3 = np.random.rand(3).tolist()
self.x0_3x3 = np.random.rand(3).tolist()
self.alpha_3x3 = 0.15
def test_convergence_2x2(self):
print("Testing 2x2 Convergence")
print(f"A: {self.A_2x2}")
print(f"b: {self.b_2x2}")
print(f"x0: {self.x0_2x2}")
richardson = modified_richardson(self.A_2x2, self.b_2x2, self.x0_2x2, self.alpha_2x2)
result = richardson()
# result = modified_richardson(self.A_2x2, self.b_2x2, self.x0_2x2, self.alpha_2x2)
expected_solution = np.linalg.solve(np.array(self.A_2x2), np.array(self.b_2x2))
print(f"Result: {result}")
print(f"Expected: {expected_solution}")
for r, e in zip(result, expected_solution):
self.assertAlmostEqual(r, e, places=4)
def test_convergence_3x3(self):
print("Testing 3x3 Convergence")
print(f"A: {self.A_3x3}")
print(f"b: {self.b_3x3}")
print(f"x0: {self.x0_3x3}")
richardson = modified_richardson(self.A_3x3, self.b_3x3, self.x0_3x3, self.alpha_3x3)
result = richardson()
# result = modified_richardson(self.A_3x3, self.b_3x3, self.x0_3x3, self.alpha_3x3)
expected_solution = np.linalg.solve(np.array(self.A_3x3), np.array(self.b_3x3))
print(f"Result: {result}")
print(f"Expected: {expected_solution}")
for r, e in zip(result, expected_solution):
self.assertAlmostEqual(r, e, places=2)
def test_invalid_alpha(self):
richardson = modified_richardson(self.A_2x2, self.b_2x2, self.x0_2x2, alpha=-0.1)
with self.assertRaises(ValueError):
richardson()
# modified_richardson(self.A_2x2, self.b_2x2, self.x0_2x2, alpha=-0.1)
def test_non_square_matrix(self):
A = [[1, 2, 3], [4, 5, 6]] # Not a square matrix
b = [7, 8]
richardson = modified_richardson(A, b, self.x0_2x2, self.alpha_2x2)
with self.assertRaises(ValueError):
richardson()
# modified_richardson(A, b, self.x0_2x2, self.alpha_2x2)
def test_dimension_mismatch(self):
b = [1, 2, 3] # Length mismatch with A_2x2
richardson = modified_richardson(self.A_2x2, b, self.x0_2x2, self.alpha_2x2)
with self.assertRaises(ValueError):
richardson()
# modified_richardson(self.A_2x2, b, self.x0_2x2, self.alpha_2x2)
def test_zero_matrix(self):
A = [[0, 0], [0, 0]]
b = [0, 0]
richardson = modified_richardson(A, b, self.x0_2x2, self.alpha_2x2)
result = richardson()
# result = modified_richardson(A, b, self.x0_2x2, self.alpha_2x2)
# Solution should be [0, 0]
print("Testing Zero Matrix")
print(f"A: {A}")
print(f"b: {b}")
print(f"Result: {result}")
self.assertEqual(result, [0, 0])
def test_large_system(self):
# A large test case designed to take a long time to converge
size = 10 #1000
A = np.random.rand(size, size) + size * np.eye(size) # Large diagonally dominant matrix
b = np.random.rand(size)
x0 = np.random.rand(size)
alpha = 0.01 / size # Small alpha to ensure convergence
print("Testing Large System")
#print(f"A: {A}")
#print(f"b: {b}")
#print(f"x0: {x0}")
richardson = modified_richardson(A.tolist(), b.tolist(), x0.tolist(), alpha, tol=1e-6, max_iter=500000)
result = richardson()
# result = modified_richardson(A.tolist(), b.tolist(), x0.tolist(), alpha, tol=1e-6, max_iter=500000)
expected_solution = np.linalg.solve(A, b)
print(f"Result: {result}")
print(f"Expected: {expected_solution}")
for r, e in zip(result, expected_solution):
self.assertAlmostEqual(r, e, places=2)
if __name__ == '__main__':
unittest.main()

100
code/richardson_abstract.py Normal file
View File

@ -0,0 +1,100 @@
from abc import ABC, abstractmethod
class modified_richardson_base(ABC):
"""
Solves the system of linear equations Ax = b using the Modified Richardson iteration method.
Parameters:
A : list of lists
Coefficient matrix (n x n).
b : list
Right-hand side vector (n).
x0 : list
Initial guess for the solution (n).
alpha : float
Relaxation parameter (0 < alpha < 2 / max(eigenvalue(A))).
tol : float, optional
Tolerance for the stopping criterion (default is 1e-6).
max_iter : int, optional
Maximum number of iterations (default is 1000).
Returns:
x : list
Approximate solution to the system of equations.
"""
def __init__(self, A, b, x0, alpha, tol=1e-6, max_iter=1000):
self.A = A
self.b = b
self.x0 = x0
self.alpha = alpha
self.tol = tol
self.max_iter = max_iter
self.n = len(A)
self.x = self.x0[:]
def check_input_data(self):
if len(self.A) != len(self.A[0]):
raise ValueError("Matrix A must be square.")
if len(self.b) != self.n:
raise ValueError("Dimension mismatch between A and b.")
if self.alpha <= 0:
raise ValueError("Alpha must be greater than 0.")
@abstractmethod
def vector_norm(self, v):
pass
@abstractmethod
def mat_vec_mult(self, mat, vec):
pass
@abstractmethod
def vec_sub(self, v1, v2):
pass
@abstractmethod
def vec_add(self, v1, v2):
pass
@abstractmethod
def vec_scalar_mult(self, scalar, vec):
pass
def __call__(self):
self.check_input_data()
x = self.x
r = self.vec_sub(self.b, self.mat_vec_mult(self.A, x))
iteration = 0
while self.vector_norm(r) > self.tol and iteration < self.max_iter:
x = self.vec_add(x, self.vec_scalar_mult(self.alpha, r))
r = self.vec_sub(self.b, self.mat_vec_mult(self.A, x))
iteration += 1
if iteration == self.max_iter:
raise ValueError("Maximum number of iterations reached before convergence")
return x
class modified_richardson(modified_richardson_base):
def vector_norm(self, v):
return sum(vi ** 2 for vi in v) ** 0.5
def mat_vec_mult(self, mat, vec):
return [sum(mat[i][j] * vec[j] for j in range(len(vec))) for i in range(len(mat))]
def vec_sub(self, v1, v2):
return [v1[i] - v2[i] for i in range(len(v1))]
def vec_add(self, v1, v2):
return [v1[i] + v2[i] for i in range(len(v1))]
def vec_scalar_mult(self, scalar, vec):
return [scalar * vi for vi in vec]