Refactor code for improved readability and maintainability

Algebra/Cholesky.py

@@ -1,24 +1,40 @@
 import numpy as np
 from math import sqrt
 
-def cholesky_decomposition(A):
+def cholesky_decomposition(A: np.ndarray) -> np.ndarray:
     """
     Perform Cholesky decomposition on a matrix A.
     A must be a symmetric, positive-definite matrix.
 
     Returns the lower triangular matrix L such that A = L * L^T.
+
+    Raises:
+        ValueError: If the input matrix is not square, not symmetric, 
+                    or not positive definite.
     """
-    n = len(A)
-    L = np.zeros((n, n))  # Initialize L with zeros
-
+    # Validate input matrix
+    if A.ndim != 2 or A.shape[0] != A.shape[1]:
+        raise ValueError("Input matrix must be square")
+
+    n = A.shape[0]
+    L = np.zeros_like(A, dtype=float)  # Initialize L with zeros
+
+    # Check matrix symmetry
+    if not np.allclose(A, A.T, atol=1e-10):
+        raise ValueError("Input matrix must be symmetric")
+
     for k in range(n):
-        # Calculate the diagonal element L[k][k]
-        sum_squares = sum(L[k][j] ** 2 for j in range(k))
-        L[k][k] = sqrt(A[k][k] - sum_squares)
-
-        # Calculate the off-diagonal elements L[i][k]
+        # Calculate diagonal element
+        diagonal_term = A[k, k] - np.sum(L[k, :k] ** 2)
+
+        # Check positive definiteness
+        if diagonal_term <= 0:
+            raise ValueError("Matrix is not positive definite")
+
+        L[k, k] = sqrt(diagonal_term)
+
+        # Calculate off-diagonal elements for row k
         for i in range(k + 1, n):
-            sum_products = sum(L[i][j] * L[k][j] for j in range(k))
-            L[i][k] = (A[i][k] - sum_products) / L[k][k]
-
+            L[i, k] = (A[i, k] - np.sum(L[i, :k] * L[k, :k])) / L[k, k]
+
     return L

Algebra/Gauss.py

@@ -1,57 +1,36 @@
-from numpy import zeros,array, eye, copy, dot
-
-def Gauss(A,b):
-    size = len(A)
-    A = array(A)
-    P = eye(size) 
-    for i in range(size-1):
-        maxRow = i
-        for j in range(i+1,size):
-            if abs(A[j][i]) > abs(A[maxRow][i]):
-                maxRow = j
-        if maxRow != i:
-            #### Swap
-            P[[maxRow,i]] = P[[i, maxRow]]
-            A[[maxRow,i]] = A[[i, maxRow]]
-        for k in range(i+1, size):
-            A[k][i] /= A[i][i]
-            for j in range(i+1, size):
-                A[k][j] -= A[k][i]*A[i][j]
-
-    L = eye(size)
-    for i in range(size):
-        for j in range(i):
-            L[i][j] = A[i][j]
-
-    U = zeros((size,size))
-    for i in range(size):
-        for j in range(i,size):
-            U[i][j] = A[i][j]
-
-    b = array(b)
-
-    y = copy(dot(P,b))
-
-    # forwards
-    for i in range(size):
-        if L[i][i] == 0: #Division by 0
-            y[i] = 0
-            continue
-        for j in range(i):
-            y[i] = (y[i] - L[i][j]*y[j])/L[i][i]
+import numpy as np
 
+def gauss_elimination(A, b):
+    A = np.asarray(A, dtype=float)
+    b = np.asarray(b, dtype=float)
 
-    x = zeros((size,1))
-    x[size-1] = y[size-1]/U[size-1][size-1]
-
-    for i in range(size-1,-1,-1):
-        if U[i][i] == 0: #Division by 0
-            x[i] = 0
-            continue
-        xv = y[i]
-        for j in range(i+1,size):
-            xv -= U[i][j]*x[j]
-        xv /= U[i][i]
-        x[i] = xv
-
-    return x
+    if A.ndim != 2 or A.shape[0] != A.shape[1]:
+        raise ValueError("Coefficient matrix A must be square")
+
+    n = A.shape[0]
+
+    if b.ndim != 1 or len(b) != n:
+        raise ValueError("Right-hand side vector b must have same length as rows of A")
+
+    augmented = np.column_stack((A, b))
+
+    for i in range(n):
+        max_row = np.argmax(np.abs(augmented[i:, i])) + i
+        if max_row != i:
+            augmented[[i, max_row]] = augmented[[max_row, i]]
+
+        pivot = augmented[i, i]
+        if np.abs(pivot) < 1e-10:
+            raise ValueError("Matrix is nearly singular, cannot solve using Gaussian elimination")
+
+        augmented[i, i:] /= pivot
+
+        for j in range(i+1, n):
+            factor = augmented[j, i]
+            augmented[j, i:] -= factor * augmented[i, i:]
+
+    x = np.zeros(n)
+    for i in range(n-1, -1, -1):
+        x[i] = augmented[i, -1] - np.sum(augmented[i, i+1:-1] * x[i+1:])
+
+    return x

Algebra/function_to_graph.py

@@ -1,41 +1,60 @@
 import matplotlib.pyplot as plt
 import numpy as np
 
-
-xmin = -100
-xmax = 100
-ymin = -100
-ymax = 100
-points = xmax - xmin
-x = np.linspace(xmin, xmax, points)
-
-
-while True:
-    eq = input("Enter the equation for y in terms of x: ")
-
-    # Use eval() to evaluate the user's input as a mathematical expression
+def get_valid_expression(prompt: str) -> str:
+    while True:
+        expr = input(prompt).strip()
+        if is_safe_expression(expr):
+            return expr
+        print("Invalid or unsafe expression. Please try again.")
+
+def is_safe_expression(expr: str) -> bool:
+    allowed_names = {'sin': np.sin, 'cos': np.cos, 'tan': np.tan,
+                     'exp': np.exp, 'log': np.log, 'sqrt': np.sqrt,
+                     'pi': np.pi, 'e': np.e}
+    code = compile(expr, '<string>', 'eval')
+    for name in code.co_names:
+        if name not in allowed_names:
+            return False
+    return True
+
+def main():
+    xmin, xmax = -100, 100
+    ymin, ymax = -100, 100
+    points = 1000
+
+    x = np.linspace(xmin, xmax, points)
+
+    eq = get_valid_expression("Enter the equation for y in terms of x: ")
+
+    allowed_names = {'sin': np.sin, 'cos': np.cos, 'tan': np.tan,
+                     'exp': np.exp, 'log': np.log, 'sqrt': np.sqrt,
+                     'pi': np.pi, 'e': np.e, 'x': x}
     try:
-        y = eval(eq)  # Evaluate the expression with x as a variable
-        break
-    except:
-        print("Invalid input equation.")
-
-
-fig, ax = plt.subplots()
-plt.axis([xmin, xmax, ymin, ymax])
-plt.plot([xmin, xmax], [0, 0], "b")
-plt.plot([0, 0], [ymin, ymax], "b")
-
-ax.set_xlabel("x values")
-ax.set_ylabel("y values")
-ax.set_title("Equation Graph")
-ax.grid(True)
-
-ax.set_xticks(np.arange(xmin, xmax, 20))
-ax.set_yticks(np.arange(ymin, ymax, 20))
-
-
-plt.plot(x, y, label=f"y= {eq}")
-
-plt.legend()
-plt.show()
+        y = eval(eq, {"__builtins__": None}, allowed_names)
+    except Exception as e:
+        print(f"Error evaluating expression: {e}")
+        return
+
+    fig, ax = plt.subplots(figsize=(10, 8))
+    plt.axis([xmin, xmax, ymin, ymax])
+
+    plt.axhline(y=0, color='blue', linestyle='-', alpha=0.5)
+    plt.axvline(x=0, color='blue', linestyle='-', alpha=0.5)
+
+    ax.set_xlabel("x values", fontsize=12)
+    ax.set_ylabel("y values", fontsize=12)
+    ax.set_title("Equation Graph", fontsize=14)
+    ax.grid(True, linestyle='--', alpha=0.7)
+
+    ax.set_xticks(np.arange(xmin, xmax+1, 20))
+    ax.set_yticks(np.arange(ymin, ymax+1, 20))
+
+    plt.plot(x, y, label=f"y = {eq}", linewidth=2)
+
+    plt.legend(fontsize=12)
+    plt.tight_layout()
+    plt.show()
+
+if __name__ == "__main__":
+    main()