From 5e68bddc64b5700b8b3924fa54e7d2deb5323eb8 Mon Sep 17 00:00:00 2001 From: Gromiusz Date: Sat, 19 Oct 2024 23:24:26 +0200 Subject: [PATCH] feat: refactor main and richardson to be modular --- code/main_abstract.py | 101 ++++++++++++++++++++++++++++++++++++ code/richardson_abstract.py | 100 +++++++++++++++++++++++++++++++++++ 2 files changed, 201 insertions(+) create mode 100644 code/main_abstract.py create mode 100644 code/richardson_abstract.py diff --git a/code/main_abstract.py b/code/main_abstract.py new file mode 100644 index 00000000..3398ee28 --- /dev/null +++ b/code/main_abstract.py @@ -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() diff --git a/code/richardson_abstract.py b/code/richardson_abstract.py new file mode 100644 index 00000000..2abff903 --- /dev/null +++ b/code/richardson_abstract.py @@ -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] + + + + + \ No newline at end of file