Add Interface for Extending Parallel Processing Methods (#4)

* Add choice between sequential or parallel(threads) processing

The LinearAlgebraUtils class has been changed to abstract class.

* remove needless imports and move the functionality from constructor to the function

* add ProcessingType class for a more transparent selection of methods

* the refactor of the assign_method

* rename the assign_method method

code/eigenvalue_methods.py

@@ -1,16 +1,14 @@
-from linear_algebra_utils import LinearAlgebraUtils
-
 class EigenvalueMethods:
     @staticmethod
-    def power_method(A, max_iter, tol=1e-6):
+    def power_method(LinAlgType, A, max_iter, tol=1e-6):
         n = len(A)
         x = [1] * n
         lambda_old = 0
         
         for _ in range(max_iter):
-            x = LinearAlgebraUtils.matrix_vector_multiply(A, x)
+            x = LinAlgType.matrix_vector_multiply(A, x)
-            lambda_new = LinearAlgebraUtils.vector_norm(x)
+            lambda_new = LinAlgType.vector_norm(x)
-            x = LinearAlgebraUtils.vector_scalar_divide(x, lambda_new)
+            x = LinAlgType.vector_scalar_divide(x, lambda_new)
             if abs(lambda_new - lambda_old) < tol:
                 break
             lambda_old = lambda_new
@@ -18,10 +16,10 @@ class EigenvalueMethods:
         return lambda_new
 
     @staticmethod
-    def inverse_power_method(A, max_iter, tol=1e-6):
+    def inverse_power_method(LinAlgType, A, max_iter, tol=1e-6):
         n = len(A)
         I = [[1 if i == j else 0 for j in range(n)] for i in range(n)]
-        A_inv = [LinearAlgebraUtils.gaussian_elimination(A.tolist(), I_col) for I_col in I]
+        A_inv = [LinAlgType.gaussian_elimination(A.tolist(), I_col) for I_col in I]
         A_inv = list(map(list, zip(*A_inv)))
 
-        return 1 / EigenvalueMethods.power_method(A_inv, max_iter, tol)
+        return 1 / EigenvalueMethods.power_method(LinAlgType, A_inv, max_iter, tol)

code/linear_algebra_utils.py

@@ -1,13 +1,71 @@
 import math
+from abc import ABC, abstractmethod
 
-class LinearAlgebraUtils:
+class LinearAlgebraUtils(ABC):
+    @staticmethod
+    @abstractmethod
+    def dot_product(v1, v2):
+        pass
+
+    @staticmethod
+    @abstractmethod
+    def matrix_vector_multiply(A, x):
+        pass
+    
+    @staticmethod
+    @abstractmethod
+    def vector_norm(v):
+        pass
+
+    @staticmethod
+    @abstractmethod
+    def vector_scalar_divide(x, scalar):
+        pass
+
+    @staticmethod
+    @abstractmethod
+    def matrix_scalar_multiply(A, w):
+        pass
+
+    @staticmethod
+    @abstractmethod
+    def vector_vector_subtraction(v1, v2):
+        pass
+
+    @staticmethod
+    @abstractmethod
+    def vector_vector_addition(v1, v2):
+        pass
+
+    @staticmethod
+    @abstractmethod
+    def scalar_matrix_multiply(omega, vector):
+        pass
+
+    @staticmethod
+    @abstractmethod
+    def matrix_norm(A):
+        pass
+
+    @staticmethod
+    @abstractmethod
+    def matrix_matrix_subtraction(A, B):
+        pass
+
+    @staticmethod
+    @abstractmethod
+    def gaussian_elimination(A, b):
+        pass
+
+
+class SequentialLinearAlgebraUtils(ABC):
     @staticmethod
     def dot_product(v1, v2):
         return sum(x*y for x, y in zip(v1, v2))
 
     @staticmethod
     def matrix_vector_multiply(A, x):
-        return [LinearAlgebraUtils.dot_product(row, x) for row in A]
+        return [SequentialLinearAlgebraUtils.dot_product(row, x) for row in A]
     
     @staticmethod
     def vector_norm(v):
@@ -69,3 +127,49 @@ class LinearAlgebraUtils:
                 M[k][-1] -= M[k][i] * x[i]
         
         return x
+    
+
+class ThreadsLinearAlgebraUtils(ABC):
+    @staticmethod
+    def dot_product(v1, v2):
+        pass
+
+    @staticmethod
+    def matrix_vector_multiply(A, x):
+        pass
+    
+    @staticmethod
+    def vector_norm(v):
+        pass
+
+    @staticmethod
+    def vector_scalar_divide(x, scalar):
+        pass
+
+    @staticmethod
+    def matrix_scalar_multiply(A, w):
+        pass
+
+    @staticmethod
+    def vector_vector_subtraction(v1, v2):
+        pass
+
+    @staticmethod
+    def vector_vector_addition(v1, v2):
+        pass
+
+    @staticmethod
+    def scalar_matrix_multiply(omega, vector):
+        pass
+
+    @staticmethod
+    def matrix_norm(A):
+        pass
+
+    @staticmethod
+    def matrix_matrix_subtraction(A, B):
+        pass
+
+    @staticmethod
+    def gaussian_elimination(A, b):
+        pass

code/processing_type.py

@@ -0,0 +1,5 @@
+from enum import Enum, auto
+
+class ProcessingType(Enum):
+    SEQUENTIAL = auto()
+    THREADS = auto()

code/richardson_method.py

@@ -1,43 +1,59 @@
 from linear_algebra_utils import LinearAlgebraUtils
+from linear_algebra_utils import SequentialLinearAlgebraUtils
+from linear_algebra_utils import ThreadsLinearAlgebraUtils
 from eigenvalue_methods import EigenvalueMethods
 from matrix_generator import MatrixGenerator
+from processing_type import ProcessingType
 
 class RichardsonMethod:
-    def __init__(self, A, b, max_iterations, size: int, x0=None, tol=1e-5):
+    def __init__(self, method: ProcessingType, A, b, max_iterations, size: int, x0=None, tol=1e-5):
+        self.LinAlg = self.assign_LinAlgType(method)
         self.A = A
         self.b = b
         self.x0 = x0 if x0 is not None else [0.0] * len(b)
         self.max_iterations = max_iterations
         self.tol = tol
         self.I = MatrixGenerator.generate_identity_matrix(size)
-        self.lambda_min, self.lambda_max = RichardsonMethod.calculate_eigenvalues(self.A, max_iterations)
+        self.lambda_min, self.lambda_max = RichardsonMethod.calculate_eigenvalues(self.LinAlg, self.A, max_iterations)
         if self.lambda_min < 0:
             raise ValueError("Matrix A is not positive semi-definite.")
         self.omega = RichardsonMethod.calculate_omega(self.lambda_min, self.lambda_max)
-
+        
     @staticmethod
-    def calculate_eigenvalues(A, max_iterations):
+    def calculate_eigenvalues(LinAlgType, A, max_iterations):
-        return EigenvalueMethods.inverse_power_method(A, max_iterations), EigenvalueMethods.power_method(A, max_iterations)
+        return EigenvalueMethods.inverse_power_method(LinAlgType, A, max_iterations), EigenvalueMethods.power_method(LinAlgType, A, max_iterations)
 
     @staticmethod
     def calculate_omega(lambda_min, lambda_max):
         return 2 / (lambda_min + lambda_max)
     
     @staticmethod
-    def convergence_norm(A, omega, I) -> bool:
+    def convergence_norm(LinAlgType, A, omega, I) -> bool:
-        wA = LinearAlgebraUtils.matrix_scalar_multiply(A, omega)
+        wA = LinAlgType.LinAlg.matrix_scalar_multiply(A, omega)
-        IMinuswA = LinearAlgebraUtils.matrix_matrix_subtraction(I, wA)
+        IMinuswA = LinAlgType.LinAlg.matrix_matrix_subtraction(I, wA)
-        norm = LinearAlgebraUtils.matrix_norm(IMinuswA)
+        norm = LinAlgType.LinAlg.matrix_norm(IMinuswA)
         return norm
+    
+    @staticmethod
+    def assign_LinAlgType(method):
+        metody = {
+            ProcessingType.SEQUENTIAL: SequentialLinearAlgebraUtils,
+            ProcessingType.THREADS: ThreadsLinearAlgebraUtils
+        }
+        
+        try:
+            return metody[method]
+        except KeyError:
+            raise ValueError("Unknown method, please use 'SEQUENTIAL' or 'THREADS'.")
 
     def solve(self):
         x = self.x0[:]
-        #if RichardsonMethod.convergence_norm(self.A, self.omega, self.I) >= 1:
+        #if RichardsonMethod.convergence_norm(self.LinAlg, self.A, self.omega, self.I) >= 1:
         #    return RichardsonMethod.convergence_norm(self.A, self.omega, self.I), "Richardson method for those values will NOT converge", 
 
         for iteration in range(self.max_iterations):
-            Ax = LinearAlgebraUtils.matrix_vector_multiply(self.A, x)
+            Ax = self.LinAlg.matrix_vector_multiply(self.A, x)
-            residual = LinearAlgebraUtils.vector_vector_subtraction(self.b, Ax)
+            residual = self.LinAlg.vector_vector_subtraction(self.b, Ax)
-            x = LinearAlgebraUtils.vector_vector_addition(x, LinearAlgebraUtils.scalar_matrix_multiply(self.omega, residual))
+            x = self.LinAlg.vector_vector_addition(x, self.LinAlg.scalar_matrix_multiply(self.omega, residual))
 
         return x, 0

code/tests.py

@@ -3,6 +3,7 @@ import numpy as np
 from scipy.sparse.linalg import cg
 from matrix_generator import MatrixGenerator
 from richardson_method import RichardsonMethod
+from processing_type import ProcessingType
 
 def calculate_norm_numpy(I, w, A):
     # Calculate the difference between I and w * A
@@ -37,7 +38,7 @@ def test_richardson_vs_cg(n: int):
     tolerance = 1e-5
     max_iterations=1000
     A, b = MatrixGenerator.generate_random_matrix_and_vector(n)
-    richardson_solver = RichardsonMethod(A, b, max_iterations, size=n, tol=1e-7)
+    richardson_solver = RichardsonMethod(ProcessingType.SEQUENTIAL , A, b, max_iterations, size=n, tol=1e-7)
     solution_richardson, info_richardson = richardson_solver.solve()
     
     solution_cg, info = cg(A, b)