Added basic linalg functions and refactored gaussian elimination algorithm

benchmarks/gaussian_elimination_and_lu_decomposition_benchmark.py

@@ -0,0 +1,184 @@
+from __future__ import annotations
+
+import argparse
+import os
+import random
+import time
+from collections.abc import Callable
+
+import matplotlib.pyplot as plt
+
+from matrix_calculus.lu_decomposition import LUDecomposition
+from matrix_calculus.matrix_processor import MatrixProcessor
+
+
+def generate_matrix(n: int, seed: int = 0) -> list[list[float]]:
+    random.seed(seed)
+
+    A = [[random.uniform(-1, 1) for _ in range(n)] for _ in range(n)]
+
+    for i in range(n):
+        A[i][i] += n  # stabilność
+
+    return A
+
+
+def measure_time(func: Callable, *args) -> float:
+    start = time.perf_counter()
+    func(*args)
+    end = time.perf_counter()
+    return end - start
+
+
+def benchmark_gaussian_elimination(out_dir: str, n_max: int):
+    sizes = list(range(10, n_max + 1, 10))
+
+    times_no_pivot = []
+    times_pivot = []
+
+    for n in sizes:
+        M = generate_matrix(n, seed=42)
+
+        t1 = measure_time(MatrixProcessor.gaussian_elimination, M, False)
+        t2 = measure_time(MatrixProcessor.gaussian_elimination, M, True)
+
+        times_no_pivot.append(t1)
+        times_pivot.append(t2)
+
+        print(f"[Gauss] n={n} done")
+
+    path = os.path.join(out_dir, "gauss_benchmark.png")
+
+    plt.figure()
+    plt.plot(sizes, times_no_pivot, label="bez pivotingu")
+    plt.plot(sizes, times_pivot, label="z pivotingiem")
+    plt.xlabel("Rozmiar macierzy (n)")
+    plt.ylabel("Czas [s]")
+    plt.title("Eliminacja Gaussa – porównanie")
+    plt.legend()
+    plt.grid()
+    plt.savefig(path)
+
+    print(f"[SAVE] {path}")
+
+
+def benchmark_lu_decomposition(out_dir: str, n_max: int):
+    sizes = list(range(10, n_max + 1, 10))
+
+    times_no_pivot = []
+    times_pivot = []
+
+    for n in sizes:
+        M = generate_matrix(n, seed=123)
+
+        t1 = measure_time(LUDecomposition.decompose, M, False)
+        t2 = measure_time(LUDecomposition.decompose, M, True)
+
+        times_no_pivot.append(t1)
+        times_pivot.append(t2)
+
+        print(f"[LU] n={n} done")
+
+    path = os.path.join(out_dir, "lu_benchmark.png")
+
+    plt.figure()
+    plt.plot(sizes, times_no_pivot, label="bez pivotingu")
+    plt.plot(sizes, times_pivot, label="z pivotingiem")
+    plt.xlabel("Rozmiar macierzy (n)")
+    plt.ylabel("Czas [s]")
+    plt.title("LU faktoryzacja – porównanie")
+    plt.legend()
+    plt.grid()
+    plt.savefig(path)
+
+    print(f"[SAVE] {path}")
+
+
+def print_example_results():
+    M = [
+        [2.0, 1.0, 1.0],
+        [4.0, -6.0, 0.0],
+        [-2.0, 7.0, 2.0],
+    ]
+
+    print("\n====================")
+    print("MACIERZ WEJŚCIOWA")
+    print("====================")
+    for row in M:
+        print(row)
+
+    print("\n====================")
+    print("GAUSS BEZ PIVOTINGU")
+    print("====================")
+    result = MatrixProcessor.gaussian_elimination(M, pivoting=False)
+    for row in result:
+        print(row)
+
+    print("\n====================")
+    print("GAUSS Z PIVOTINGIEM (bez normalizacji przekątnej)")
+    print("====================")
+    result = MatrixProcessor.gaussian_elimination(
+        M,
+        pivoting=True,
+        normalize_diagonal=False,
+    )
+    for row in result:
+        print(row)
+
+    print("\n====================")
+    print("LU BEZ PIVOTINGU")
+    print("====================")
+    L, U = LUDecomposition.decompose(M, pivoting=False)
+    print("L:")
+    for row in L:
+        print(row)
+    print("U:")
+    for row in U:
+        print(row)
+
+    print("\n====================")
+    print("LU Z PIVOTINGIEM")
+    print("====================")
+    P, L, U = LUDecomposition.decompose(M, pivoting=True)
+    print("P:")
+    for row in P:
+        print(row)
+    print("L:")
+    for row in L:
+        print(row)
+    print("U:")
+    for row in U:
+        print(row)
+
+
+def benchmark_gaussian_elimination_and_lu_decomposition(out_dir: str, n_max: int):
+    os.makedirs(out_dir, exist_ok=True)
+
+    benchmark_gaussian_elimination(out_dir, n_max)
+    benchmark_lu_decomposition(out_dir, n_max)
+    print_example_results()
+
+
+if __name__ == "__main__":
+    parser = argparse.ArgumentParser()
+
+    parser.add_argument(
+        "--output",
+        type=str,
+        default="results",
+        help="folder na wykresy",
+    )
+
+    parser.add_argument(
+        "--n",
+        type=int,
+        default=150,
+        help="maksymalny rozmiar macierzy",
+    )
+
+    args = parser.parse_args()
+
+    benchmark_gaussian_elimination_and_lu_decomposition(
+        out_dir=args.output,
+        n_max=args.n,
+    )

matrix_calculus/linalg/__init__.py

@@ -0,0 +1,6 @@
+from __future__ import annotations
+
+from .matrix_properties import det
+from .matrix_properties import rank
+
+__all__ = ["det", "rank"]

matrix_calculus/linalg/matrix_properties.py

@@ -0,0 +1,32 @@
+from __future__ import annotations
+
+from ..matrix_processor import MatrixProcessor
+
+
+def det(matrix: list[list[float]]) -> float:
+    """Calculate the determinant of a square matrix."""
+    processed_matrix, stats = MatrixProcessor.gaussian_elimination(
+        matrix,
+        stats=True,
+        normalize_diagonal=False,
+    )
+
+    if not stats["reversible"]:
+        return 0.0
+
+    determinant = (-1) ** stats["rows_swaps"]
+    for i in range(len(processed_matrix)):
+        determinant *= processed_matrix[i][i]
+
+    return determinant
+
+
+def rank(matrix: list[list[float]]) -> float:
+    """Calculate the rank of a matrix."""
+    _, stats = MatrixProcessor.gaussian_elimination(
+        matrix,
+        stats=True,
+        normalize_diagonal=False,
+    )
+
+    return stats["rank"]

matrix_calculus/matrix_processor.py

@@ -1,12 +1,61 @@
 from __future__ import annotations
 
+from typing import Any
+from typing import Literal
+from typing import overload
+
 
 class MatrixProcessor:
+    @overload
+    @staticmethod
+    def gaussian_elimination(
+        matrix: list[list[float]],
+        pivoting: bool = ...,
+        stats: Literal[True] = ...,
+        normalize_diagonal: bool = ...,
+    ) -> tuple[list[list[float]], dict[str, Any]]:
+        pass
+
+    @overload
     @staticmethod
-    def gaussian_elimination(matrix: list[list[float]]) -> list[list[float]]:
+    def gaussian_elimination(
+        matrix: list[list[float]],
+        pivoting: bool = ...,
+        stats: Literal[False] = ...,
+        normalize_diagonal: bool = ...,
+    ) -> list[list[float]]:
+        pass
+
+    @staticmethod
+    def gaussian_elimination(
+        matrix: list[list[float]],
+        pivoting: bool = True,
+        stats: bool = False,
+        normalize_diagonal: bool = True,
+    ) -> list[list[float]] | tuple[list[list[float]], dict]:
+        if pivoting:
+            result = MatrixProcessor.gaussian_elimination_with_pivoting(
+                matrix,
+                normalize_diagonal,
+            )
+        else:
+            result = MatrixProcessor.gaussian_elimination_without_pivoting(
+                matrix,
+                normalize_diagonal,
+            )
+
+        return result[0] if not stats else result
+
+    @staticmethod
+    def gaussian_elimination_without_pivoting(
+        matrix: list[list[float]],
+        normalize_diagonal: bool = True,
+    ) -> tuple[list[list[float]], dict]:
         M = [row[:] for row in matrix]
         n = len(M)
         eps = 1e-12
+        rows_swaps = 0
+        rank = n
 
         def swap_rows(row1, row2):
             for col in range(n):
@@ -22,25 +71,37 @@ def swap_rows(row1, row2):
                         diagonal_element_value = M[row][diagonal_element]
 
                         swap_rows(diagonal_element, row)
+                        rows_swaps += 1
 
             if abs(diagonal_element_value) < eps:
-                raise ValueError("Matrix is singular")
+                rank -= 1
+                continue
 
-            for col in range(diagonal_element, n):
-                M[diagonal_element][col] /= diagonal_element_value
+            if normalize_diagonal:
+                for col in range(diagonal_element, n):
+                    M[diagonal_element][col] /= diagonal_element_value
 
             for row in range(diagonal_element + 1, n):
-                factor = M[row][diagonal_element]
+                factor = (
+                    M[row][diagonal_element]
+                    if normalize_diagonal
+                    else M[row][diagonal_element] / diagonal_element_value
+                )
                 for col in range(diagonal_element, n):
                     M[row][col] -= factor * M[diagonal_element][col]
 
-        return M
+        return M, {"rows_swaps": rows_swaps, "rank": rank, "reversible": rank == n}
 
     @staticmethod
-    def gaussian_elimination_pivoting(matrix: list[list[float]]) -> list[list[float]]:
+    def gaussian_elimination_with_pivoting(
+        matrix: list[list[float]],
+        normalize_diagonal: bool = True,
+    ) -> tuple[list[list[float]], dict]:
         M = [row[:] for row in matrix]
         n = len(M)
         eps = 1e-12
+        rows_swaps = 0
+        rank = n
 
         def swap_rows(row1, row2):
             for col in range(n):
@@ -57,14 +118,24 @@ def swap_rows(row1, row2):
 
             if highest_pivot_row != diagonal_element:
                 swap_rows(diagonal_element, highest_pivot_row)
+                rows_swaps += 1
 
             pivot_value = M[diagonal_element][diagonal_element]
             if abs(pivot_value) < eps:
-                raise ValueError("Matrix is singular")
+                rank -= 1
+                continue
+
+            if normalize_diagonal:
+                for col in range(diagonal_element, n):
+                    M[diagonal_element][col] /= pivot_value
 
             for row in range(diagonal_element + 1, n):
-                factor = M[row][diagonal_element] / pivot_value
+                factor = (
+                    M[row][diagonal_element]
+                    if normalize_diagonal
+                    else M[row][diagonal_element] / pivot_value
+                )
                 for col in range(diagonal_element, n):
                     M[row][col] -= factor * M[diagonal_element][col]
 
-        return M
+        return M, {"rows_swaps": rows_swaps, "rank": rank, "reversible": rank == n}