kuhyx/WUT_Computer_Science
1
0
Fork 0
mirror of https://github.com/kuhyx/WUT_Computer_Science.git synced 2026-07-04 15:23:11 +02:00
Code Issues Actions Packages Projects Releases Wiki Activity

Revert "Merge branch 'main' of https://github.com/kuhyx/PORR"

Browse Source 
This reverts commit a7123acced, reversing
changes made to 4512cf9305.
This commit is contained in:
Krzysztof Rudnicki 2024-11-28 17:48:11 +01:00
parent a7123acced
commit 87d6bb0d97
7 changed files with 224 additions and 114 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

43
code/eigenvalue_methods.py Normal file
View File

@ -0,0 +1,43 @@
import numpy as np
import scipy
class EigenvalueMethods:
@staticmethod
def get_sing_vals(file_path):
mat_contents = scipy.io.loadmat(file_path)
A = mat_contents['S'][0][0] # Pobranie pierwszego elementu z pola 'S'
singular_values = A['s'].flatten()
return singular_values
@staticmethod
def power_method(LinAlgType, A, type, max_iter=100, tol=1e-6):
if (type == 'nemeth12'):
singular_vals = EigenvalueMethods.get_sing_vals("nemeth12_SVD.mat")
return np.max(singular_vals)
n = len(A)
x = [1] * n
lambda_old = 0
for _ in range(max_iter):
x = LinAlgType.matrix_vector_multiply(A, x)
lambda_new = LinAlgType.vector_norm(x)
x = LinAlgType.vector_scalar_divide(x, lambda_new)
if abs(lambda_new - lambda_old) < tol:
break
lambda_old = lambda_new
return lambda_new
@staticmethod
def inverse_power_method(LinAlgType, A, type, max_iter=100, tol=1e-6):
if (type == 'nemeth12'):
singular_vals = EigenvalueMethods.get_sing_vals("nemeth12_SVD.mat")
return np.min(singular_vals)
n = len(A)
I = [[1 if i == j else 0 for j in range(n)] for i in range(n)]
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(LinAlgType, A_inv, type, max_iter, tol)

141
code/linear_algebra_utils.py
View File

@ -6,7 +6,6 @@ from abc import ABC, abstractmethod
from concurrent.futures import ThreadPoolExecutor
from functools import partial
from time_measurement import time_measurement, time_accumulator
import numpy as np
class LinearAlgebraUtils(ABC):
@staticmethod
@ -59,6 +58,11 @@ class LinearAlgebraUtils(ABC):
def matrix_matrix_subtraction(A, B):
pass
@staticmethod
@abstractmethod
def gaussian_elimination(A, b):
pass
class SequentialLinearAlgebraUtils(ABC):
@staticmethod
@ -79,7 +83,7 @@ class SequentialLinearAlgebraUtils(ABC):
@staticmethod
def matrix_scalar_multiply(A, w):
return A * w
return [[w * A[i][j] for j in range(len(A[0]))] for i in range(len(A))]
@staticmethod
def vector_vector_subtraction(v1, v2):
@ -101,12 +105,39 @@ class SequentialLinearAlgebraUtils(ABC):
def matrix_matrix_subtraction(A, B):
return [[A[i][j] - B[i][j] for j in range(len(A[0]))] for i in range(len(A))]
@staticmethod
def gaussian_elimination(A, b):
n = len(A)
M = [row[:] for row in A]
for i in range(n):
M[i].append(b[i])
for k in range(n):
if M[k][k] == 0:
for i in range(k + 1, n):
if M[i][k] != 0:
M[k], M[i] = M[i], M[k]
break
for i in range(k + 1, n):
factor = M[i][k] / M[k][k]
for j in range(k, n + 1):
M[i][j] -= factor * M[k][j]
x = [0] * n
for i in range(n - 1, -1, -1):
x[i] = M[i][-1] / M[i][i]
for k in range(i - 1, -1, -1):
M[k][-1] -= M[k][i] * x[i]
return x
class ThreadsLinearAlgebraUtils(ABC):
NUM_THREADS = 4
@staticmethod
@time_measurement(time_accumulator)
def get_chunk_size(data):
num_elements = len(data)
num_threads = min(ThreadsLinearAlgebraUtils.NUM_THREADS, num_elements)
@ -116,7 +147,6 @@ class ThreadsLinearAlgebraUtils(ABC):
@staticmethod
@time_measurement(time_accumulator)
def divide_vectors_to_chunks(v1, v2):
chunk_size, num_threads, remainder = ThreadsLinearAlgebraUtils.get_chunk_size(v1)
@ -130,7 +160,6 @@ class ThreadsLinearAlgebraUtils(ABC):
return chunks
@staticmethod
@time_measurement(time_accumulator)
def divide_vector_or_matrix_to_chunks(v):
chunk_size, num_threads, remainder = ThreadsLinearAlgebraUtils.get_chunk_size(v)
@ -234,23 +263,8 @@ class ThreadsLinearAlgebraUtils(ABC):
@staticmethod
@time_measurement(time_accumulator)
def divide_matrixes_to_chunks(A, B):
num_rows = len(A)
num_threads = ThreadsLinearAlgebraUtils.NUM_THREADS
if num_threads > num_rows:
num_threads = num_rows
if num_rows == 0:
return []
chunk_size = num_rows // num_threads
remainder = num_rows % num_threads
chunks = []
start = 0
for _ in range(num_threads):
end = start + chunk_size + (1 if remainder > 0 else 0)
chunks.append((A[start:end], B[start:end]))
start = end
if remainder > 0:
remainder -= 1
return chunks
chunk_size = len(A) // ThreadsLinearAlgebraUtils.NUM_THREADS
return [(A[i:i + chunk_size], B[i:i + chunk_size]) for i in range(0, len(A), chunk_size)]
@staticmethod
@time_measurement(time_accumulator)
@ -265,6 +279,40 @@ class ThreadsLinearAlgebraUtils(ABC):
results = executor.map(subtract_chunk, chunks)
return [row for chunk in results for row in chunk]
@staticmethod
@time_measurement(time_accumulator)
def gaussian_elimination(A, b):
n = len(A)
M = [row[:] for row in A]
for i in range(n):
M[i].append(b[i])
for k in range(n):
# Pivoting
if M[k][k] == 0:
for i in range(k + 1, n):
if M[i][k] != 0:
M[k], M[i] = M[i], M[k]
break
# Threads
def eliminate_row(i):
factor = M[i][k] / M[k][k]
for j in range(k, n + 1):
M[i][j] -= factor * M[k][j]
with ThreadPoolExecutor(max_workers=ThreadsLinearAlgebraUtils.NUM_THREADS) as executor:
rows_to_eliminate = range(k + 1, n)
executor.map(eliminate_row, rows_to_eliminate)
x = [0] * n
for i in range(n - 1, -1, -1):
x[i] = M[i][-1] / M[i][i]
for k in range(i - 1, -1, -1):
M[k][-1] -= M[k][i] * x[i]
return x
@time_measurement(time_accumulator)
def process_row(params):
@ -282,7 +330,6 @@ def multiply_by_scalar(pair):
element, scalar = pair
return element * scalar
class ProcessLinearAlgebraUtils:
@staticmethod
@time_measurement(time_accumulator)
@ -375,26 +422,52 @@ class ProcessLinearAlgebraUtils:
result = pool.map(multiply_by_scalar, [(element, omega) for element in vector])
return list(result)
@staticmethod
@time_measurement(time_accumulator)
def sum_of_squares(row):
return sum(x ** 2 for x in row)
@staticmethod
@time_measurement(time_accumulator)
def matrix_norm(A):
with Pool() as pool:
row_sums = pool.map(ProcessLinearAlgebraUtils.sum_of_squares, A)
row_sums = pool.map(lambda row: sum(x ** 2 for x in row), A)
return math.sqrt(sum(row_sums))
@staticmethod
@time_measurement(time_accumulator)
def subtract_rows(row_from_A, row_from_B):
return [a - b for a, b in zip(row_from_A, row_from_B)]
def matrix_matrix_subtraction(A, B):
def subtract_rows(row_pair):
return [a - b for a, b in zip(*row_pair)]
with Pool() as pool:
result = pool.starmap(subtract_rows, zip(A, B))
return result
@staticmethod
@time_measurement(time_accumulator)
def matrix_matrix_subtraction(A, B):
def gaussian_elimination(A, b):
try:
n = len(A)
A = [list(row) + [b_i] for row, b_i in zip(A, b)]
for k in range(n):
# Pivoting
max_index = max(range(k, n), key=lambda x: abs(A[x][k]))
if A[max_index][k] == 0:
raise ValueError("Matrix is singular and cannot be solved.")
A[k], A[max_index] = A[max_index], A[k]
# Parallel row processing
with Pool() as pool:
result = pool.starmap(ProcessLinearAlgebraUtils.subtract_rows, zip(A, B))
return result
results = pool.map(process_row, [(A, k, i) for i in range(k + 1, n)])
# Update remaining rows in matrix
for i in range(k + 1, n):
A[i] = results[i - (k + 1)]
# Back substitution
x = [0] * n
for i in range(n - 1, -1, -1):
sum_ax = sum(A[i][j] * x[j] for j in range(i + 1, n))
x[i] = (A[i][-1] - sum_ax) / A[i][i]
return x
except Exception as e:
print(f"Error during Gaussian elimination: {e}")
return None

51
code/matrix_generator.py
View File

@ -1,12 +1,20 @@
import numpy as np
import scipy.io
import os
class MatrixGenerator:
@staticmethod
def generate_spd_matrix(n: int) -> np.ndarray:
"""
Generates a random symmetric positive definite matrix of size n x n.
Parameters:
n (int): The size of the matrix to generate.
Returns:
np.ndarray: A symmetric positive definite matrix of size n x n.
"""
A = np.random.rand(n, n)
spd_matrix = np.dot(A, A.T) + n * MatrixGenerator.generate_identity_matrix(n) # Adding n*I ensures positive definiteness
spd_matrix = np.dot(A, A.T) + n * np.eye(n) # Adding n*I ensures positive definiteness
return spd_matrix
@staticmethod
@ -33,54 +41,17 @@ class MatrixGenerator:
if type == 'spd':
if size is None:
raise ValueError("Size must be provided for SPD matrix generation.")
try:
matrix, vector, lambda_min, lambda_max = MatrixGenerator.load_from_file("spd_"+str(size)+".npz")
except FileNotFoundError as e:
matrix = MatrixGenerator.generate_spd_matrix(size)
vector = np.random.uniform(-1, 1, size)
lambda_min, lambda_max = MatrixGenerator.calculate_eigenvalues(matrix)
MatrixGenerator.save_to_file(matrix, vector, lambda_min, lambda_max, "spd_"+str(size)+".npz")
elif type == 'nemeth12':
try:
matrix, vector, lambda_min, lambda_max = MatrixGenerator.load_from_file("nemeth12.npz")
except FileNotFoundError as e:
matrix = -1 * MatrixGenerator.get_matrix_from_file("nemeth12.mat", 1)
size = matrix.shape[0]
vector = MatrixGenerator.generate_alternating_vector(size)
lambda_min, lambda_max = MatrixGenerator.calculate_eigenvalues(matrix)
MatrixGenerator.save_to_file(matrix, vector, lambda_min, lambda_max, "nemeth12.npz")
elif type == 'poli3':
try:
matrix, vector, lambda_min, lambda_max = MatrixGenerator.load_from_file("poli3.npz")
except FileNotFoundError as e:
matrix = MatrixGenerator.get_matrix_from_file("poli3.mat", 2)
size = matrix.shape[0]
vector = MatrixGenerator.generate_alternating_vector(size)
lambda_min, lambda_max = MatrixGenerator.calculate_eigenvalues(matrix)
MatrixGenerator.save_to_file(matrix, vector, lambda_min, lambda_max, "poli3.npz")
else:
raise ValueError("Invalid type specified. Choose 'spd', 'nemeth12', or 'poli3'.")
return matrix, vector, lambda_min, lambda_max
@staticmethod
def calculate_eigenvalues(A):
eigenvalues = np.linalg.eigvals(A)
lambda_min = np.min(eigenvalues)
lambda_max = np.max(eigenvalues)
return lambda_min, lambda_max
@staticmethod
def save_to_file(matrix, vector, lambda_min, lambda_max, file_path):
np.savez(file_path, matrix=matrix, vector=vector, lambda_min=lambda_min, lambda_max=lambda_max)
def load_from_file(file_path):
if not os.path.exists(file_path):
raise FileNotFoundError(f"The file {file_path} does not exist.")
data = np.load(file_path)
matrix = data['matrix']
vector = data['vector']
lambda_min = data['lambda_min']
lambda_max = data['lambda_max']
return matrix, vector, lambda_min, lambda_max
return matrix, vector

BIN
code/nemeth12_SVD.mat Normal file
View File

Binary file not shown.

BIN
code/poli3_SVD.mat Normal file
View File

Binary file not shown.

21
code/richardson_method.py
View File

@ -1,14 +1,14 @@
import linear_algebra_utils as linAlg
from eigenvalue_methods import EigenvalueMethods
from matrix_generator import MatrixGenerator
from processing_type import ProcessingType
from time_measurement import time_measurement, time_accumulator, tests_time
import time
import gc
import numpy as np
class RichardsonMethod:
@time_measurement(time_accumulator)
def __init__(self, method: ProcessingType, A, b, lambda_min, lambda_max, max_iterations, size: int, x0=None, tol=1e-5):
def __init__(self, method: ProcessingType, A, type, b, max_iterations, size: int, x0=None, tol=1e-5):
self.LinAlg = self.assign_LinAlgType(method)
self.A = A
self.b = b
@ -16,21 +16,24 @@ class RichardsonMethod:
self.max_iterations = max_iterations
self.tol = tol
# self.I = MatrixGenerator.generate_identity_matrix(size)
self.lambda_min = lambda_min
self.lambda_max = lambda_max
self.lambda_min, self.lambda_max = RichardsonMethod.calculate_eigenvalues(self.LinAlg, self.A, type)
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(LinAlgType, A, type):
return EigenvalueMethods.inverse_power_method(LinAlgType, A, type), EigenvalueMethods.power_method(LinAlgType, A, type)
@staticmethod
def calculate_omega(lambda_min, lambda_max):
return 2 / (lambda_min + lambda_max)
@staticmethod
def convergence_norm(LinAlgType, A, omega, I) -> bool:
wA = LinAlgType.matrix_scalar_multiply(A, omega)
IMinuswA = LinAlgType.matrix_matrix_subtraction(I, wA)
norm = LinAlgType.matrix_norm(IMinuswA)
wA = LinAlgType.LinAlg.matrix_scalar_multiply(A, omega)
IMinuswA = LinAlgType.LinAlg.matrix_matrix_subtraction(I, wA)
norm = LinAlgType.LinAlg.matrix_norm(IMinuswA)
return norm
@staticmethod
@ -51,8 +54,8 @@ class RichardsonMethod:
time_accumulator.total_time = 0
start = time.perf_counter()
x = self.x0[:]
# if RichardsonMethod.convergence_norm(self.LinAlg, self.A, self.omega, self.I) >= 1:
# return RichardsonMethod.convergence_norm(self.LinAlg, self.A, self.omega, self.I), "Richardson method for those values will NOT converge",
#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 = self.LinAlg.matrix_vector_multiply(self.A, x)

62
code/tests.py
View File

@ -1,52 +1,59 @@
import pytest
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
from time_measurement import time_measurement, tests_time
def calculate_norm_numpy(I, w, A):
# Calculate the difference between I and w * A
difference = I - w * A
# Calculate the Euclidean norm of the difference
norm = np.linalg.norm(difference)
return norm
def calculate_eigenvalues(A):
# Calculate the eigenvalues of matrix A
eigenvalues = np.linalg.eigvals(A)
# Find the minimum and maximum eigenvalues
lambda_min = np.min(eigenvalues)
lambda_max = np.max(eigenvalues)
return lambda_min, lambda_max
def calcualte_norm_from_matrix_numpy(A, n):
def calcualte_norm_from_matrix_numpy(A, n, max_iterations):
lambda_min, lambda_max = calculate_eigenvalues(A)
omega = 2 / (lambda_min + lambda_max)
I = np.eye(n)
return calculate_norm_numpy(I, omega, A)
@time_measurement(tests_time)
def solution_lib(A, b):
return np.linalg.solve(A, b)
@pytest.mark.parametrize("n", [2, 5, 10, 50, 100, 300, 500, 750, 1000, 5000, 10000])
@pytest.mark.parametrize("n", [2, 3, 4, 5, 10, 20, 50, 100])
@pytest.mark.parametrize("processing_type", [ProcessingType.SEQUENTIAL, ProcessingType.THREADS, ProcessingType.PROCESSES])
@pytest.mark.parametrize("matrix_type", ["spd", "nemeth12", "poli3"])
@pytest.mark.parametrize("matrix_type", ["spd", "nemeth12"])#, "poli3"])
@time_measurement(tests_time)
def test_richardson_vs_cg(n: int, processing_type: ProcessingType, matrix_type: str, capsys):
print("matrix type: ", matrix_type)
print("matrix size: ", n if matrix_type == "spd" else "fixed")
tolerance = 8e-3
tolerance = 7e-3
max_iterations=100
if matrix_type in ["nemeth12", "poli3"] and n != 2:
pytest.skip("Fixed matrix size for nemeth12 and poli3, skipping redundant runs.")
if matrix_type == "spd":
A, b, lambda_min, lambda_max = MatrixGenerator.generate_matrix_and_vector('spd', size=n)
A, b = MatrixGenerator.generate_matrix_and_vector('spd', size=n)
elif matrix_type == "poli3":
A, b, lambda_min, lambda_max = MatrixGenerator.generate_matrix_and_vector('poli3')
A, b = MatrixGenerator.generate_matrix_and_vector('poli3')
elif matrix_type == "nemeth12":
A, b, lambda_min, lambda_max = MatrixGenerator.generate_matrix_and_vector('nemeth12')
A, b = MatrixGenerator.generate_matrix_and_vector('nemeth12')
else:
raise ValueError("Invalid matrix type specified. Choose 'spd', 'poli3', or 'nemeth12'.")
richardson_solver = RichardsonMethod(processing_type, A, b, lambda_min, lambda_max, max_iterations, size=A.shape[0], tol=1e-7)
richardson_solver = RichardsonMethod(processing_type, A, matrix_type, b, max_iterations, size=A.shape[0], tol=1e-7)
# solution_richardson, info_richardson = richardson_solver.solve()
solution_richardson, info_richardson = None, None
with capsys.disabled():
@ -56,24 +63,37 @@ def test_richardson_vs_cg(n: int, processing_type: ProcessingType, matrix_type:
captured = capsys.readouterr()
print("Captured output:", captured.out)
solution = solution_lib(A,b)
solution_cg, info = cg(A, b, atol=0.)
assert_converged(solution_richardson, info_richardson, solution, tolerance, A, n)
if info == 0: # SciPy CG converged
assert_scipy_converged(solution_richardson, info_richardson, solution_cg, tolerance, A, b, max_iterations, n)
else: # SciPy CG did not converge
assert_scipy_not_converged(solution_richardson, info_richardson, A, b)
def assert_converged(solution_richardson, info_richardson, solution, tolerance, A, n):
def assert_scipy_converged(solution_richardson, info_richardson, solution_cg, tolerance, A, b, max_iterations, n):
if info_richardson == "Richardson method for those values will NOT converge":
numpy_norm = calcualte_norm_from_matrix_numpy(A, n)
print("Richardson did not converge, while SciPy did")
numpy_norm = calcualte_norm_from_matrix_numpy(A, n, max_iterations)
print("Numpy norm: ", numpy_norm, " Richardson norm: ", solution_richardson)
assert False, "Richardson did not converge"
assert False, "Richardson did not converge, while SciPy did"
else:
difference = np.linalg.norm(solution_richardson - solution)
print(f"Difference between Richardson and numpy solutions: {difference:.8f}")
difference = np.linalg.norm(solution_richardson - solution_cg)
print(f"Difference between Richardson and CG solutions: {difference:.8f}")
if difference < tolerance:
print("Both Richardson and numpy method converged and calculated correct values.")
print("Both Richardson and CG converged and calculated correct values.")
else:
print("Solution numpy:\n", solution)
print("Solution CG:\n", solution_cg)
print("Solution Richardson:\n", solution_richardson)
assert difference < tolerance, f"The solutions are different! Difference: {difference:.8f}"
def assert_scipy_not_converged(solution_richardson, info_richardson, A, b):
if info_richardson == "Richardson method for those values will NOT converge":
print("Richardson and SciPy did not converge")
else:
print("Richardson converged while SciPy did not:", solution_richardson)
print("Matrix A:\n", A)
print("Vector b:\n", b)
assert False, "Richardson converged while SciPy did not"
if __name__ == "__main__":
pytest.main()