Add sketch of run.py file, add script for profiling, numba for power method

Author: lorenzotomada

Makefile

@@ -19,5 +19,4 @@ format:
 	python -m black .
 
 all: tests docs format
-	@echo "Tests completed, documentation generated, and code formatted."
-
+	@echo "Tests completed, documentation generated, and code formatted."

experiments/config.yaml

@@ -0,0 +1,3 @@
+dimension: 1000
+
+density: 0.1

requirements.txt

@@ -30,6 +30,7 @@ Pygments==2.19.1
 pyparsing==3.2.1
 pytest==8.3.4
 python-dateutil==2.9.0.post0
+PyYAML==6.0.2
 requests==2.32.3
 scipy==1.15.2
 six==1.17.0

scripts/run.py

@@ -0,0 +1,45 @@
+from pyclassify import (
+    eigenvalues_np,
+    eigenvalues_sp,
+    power_method,
+    power_method_numba,
+)
+from pyclassify.utils import make_symmetric, read_config
+import numpy as np
+import scipy.sparse as sp
+import random
+import argparse
+
+
+random.seed(226)
+
+
+parser = argparse.ArgumentParser()
+parser.add_argument("--config", type=str, required=False, help="config file:")
+
+
+args = parser.parse_args()
+filename = (
+    args.config if args.config else "./experiments/config"
+)  # automatic choice if no argument is passed
+
+
+kwargs = read_config(filename)
+dim = kwargs["dimension"]
+density = kwargs["density"]
+
+
+matrix = sp.random(dim, dim, density=density, format="csr")
+matrix = make_symmetric(matrix)
+
+eigs_np = eigenvalues_np(matrix.toarray(), symmetric=True)
+eigs_sp = eigenvalues_sp(matrix, symmetric=True)
+
+index_np = np.argmax(np.abs(eigs_np))
+index_sp = np.argmax(np.abs(eigs_sp))
+
+biggest_eigenvalue_np = eigs_np[index_np]
+biggest_eigenvalue_sp = eigs_sp[index_sp]
+
+biggest_eigenvalue_pm = power_method(matrix)
+biggest_eigenvalue_pm_numba = power_method_numba(matrix.toarray())

shell/profile.sh

@@ -0,0 +1,3 @@
+#!/bin/bash
+python -m kernprof -l -o logs/profile_eigenvalues.dat scripts/run.py --config=experiments/config
+python -m line_profiler -rmt "logs/profile_eigenvalues.dat" > logs/eigenvalues.txt

src/pyclassify/__init__.py

@@ -0,0 +1,8 @@
+__all__ = [
+   'eigenvalues_np',
+   'eigenvalues_sp',
+    'power_method',
+    'power_method_numba',
+   ]
+
+from .eigenvalues import eigenvalues_np, eigenvalues_sp, power_method, power_method_numba

src/pyclassify/eigenvalues.py

@@ -0,0 +1,49 @@
+import numpy as np
+import scipy.sparse as sp
+from line_profiler import profile
+from numpy.linalg import eig, eigh
+from pyclassify.utils import check_A_square_matrix, power_method_numba_helper
+
+@profile
+def eigenvalues_np(A, symmetric=True):
+    check_A_square_matrix(A)
+    eigenvalues, _ = eigh(A) if symmetric else eig(A)
+    return eigenvalues
+
+
+@profile
+def eigenvalues_sp(A, symmetric=True):
+    check_A_square_matrix(A)
+    eigenvalues, _ = (
+        sp.linalg.eigsh(A, k=A.shape[0] - 1)
+        if symmetric
+        else sp.linalg.eigs(A, k=A.shape[0] - 1)
+    )
+    return eigenvalues
+
+
+@profile
+def power_method(A, max_iter=500, tol=1e-4, x=None):
+    check_A_square_matrix(A)
+    if x is None:
+        x = np.random.rand(A.shape[0])
+    x /= np.linalg.norm(x)
+    x_old = x
+
+    iteration = 0
+    update_norm = tol + 1
+
+    while iteration < max_iter and update_norm > tol:
+        x = A @ x
+        x /= np.linalg.norm(x)
+        update_norm = np.linalg.norm(x - x_old) / np.linalg.norm(x_old)
+        x_old = x.copy()
+        iteration += 1
+
+    return x @ A @ x
+
+
+
+@profile
+def power_method_numba(A):
+    return power_method_numba_helper(A)

src/pyclassify/utils.py

@@ -1,44 +1,33 @@
 import numpy as np
 import scipy.sparse as sp
-from numpy.linalg import eig, eigh
+import numba
+import os
+import yaml
+from line_profiler import profile
 
+# from numba.pycc import CC
 
 
-def make_symmetric(A):
-    if not sp.isspmatrix(A):
-        raise TypeError("Input must be a scipy sparse matrix")
-    
-    A_sym = (A + A.T) / 2
-    return A_sym
-
-
-
-def eigenvalues_np(A, symmetric=True):
-    if A.shape[0] != A.shape[1]:
-        raise ValueError("Matrix must be square!")
+def check_A_square_matrix(A):
     if not isinstance(A, (np.ndarray, sp.spmatrix)):
         raise TypeError("Input matrix must be a NumPy array or a SciPy sparse matrix!")
-    
-    eigenvalues, _ = eigh(A) if symmetric else eig(A)
-    return eigenvalues
-
-
-
-def eigenvalues_sp(A, symmetric=True):
     if A.shape[0] != A.shape[1]:
         raise ValueError("Matrix must be square!")
-    if not isinstance(A, (np.ndarray, sp.spmatrix)):
-        raise TypeError("Input matrix must be a NumPy array or a SciPy sparse matrix!")
-
-    eigenvalues, _ = sp.linalg.eigsh(A, k=A.shape[0] - 1) if symmetric else sp.linalg.eigs(A, k=A.shape[0] - 1)
-    return eigenvalues
 
 
+@profile
+def make_symmetric(A):
+    check_A_square_matrix(A)
+    A_sym = (A + A.T) / 2
+    return A_sym
 
-def power_method(A, max_iter=1000, tol=1e-4):
-    assert A.shape[0] == A.shape[1], "Matrix must be square"
 
-    x = np.random.rand(A.shape[0])
+@numba.njit(nogil=True, parallel=True)
+def power_method_numba_helper(A, max_iter=500, tol=1e-4, x=None):
+    if A.shape[0] != A.shape[1]:
+        raise ValueError("Matrix must be square!")  # explain why re-written
+    if x is None:
+        x = np.random.rand(A.shape[0])
     x /= np.linalg.norm(x)
     x_old = x
 
@@ -48,8 +37,22 @@ def power_method(A, max_iter=1000, tol=1e-4):
     while iteration < max_iter and update_norm > tol:
         x = A @ x
         x /= np.linalg.norm(x)
-        update_norm = np.linalg.norm(x - x_old)/np.linalg.norm(x_old)
+        update_norm = np.linalg.norm(x - x_old) / np.linalg.norm(x_old)
         x_old = x.copy()
         iteration += 1
 
-    return x @ A @ x
+    return x @ A @ x
+
+
+def read_config(file: str) -> dict:
+    """
+    To read the desired configuration file, passed in input as a string
+    Input:
+        file: str (representing the location of file)
+    Returns:
+        dict (containing the configuration parameters needed)
+    """
+    filepath = os.path.abspath(f"{file}.yaml")
+    with open(filepath, "r") as stream:
+        kwargs = yaml.safe_load(stream)
+    return kwargs

test/test_.py

@@ -1,32 +1,47 @@
 import numpy as np
 import scipy.sparse as sp
 import pytest
-from pyclassify.utils import eigenvalues_np, eigenvalues_sp, power_method, make_symmetric
+from pyclassify import (
+    eigenvalues_np,
+    eigenvalues_sp,
+    power_method,
+    power_method_numba,
+)
+from pyclassify.utils import check_A_square_matrix, make_symmetric
 
 
 np.random.seed(42)
-sizes = [10, 100, 1000]
-densities = [0.1, 0.5]
 
 
+def test_checks_on_A():
+    ugly_nonsquare_matrix = np.random.rand(5, 3)
+
+    with pytest.raises(ValueError):
+        check_A_square_matrix(ugly_nonsquare_matrix)
+    with pytest.raises(TypeError):
+        check_A_square_matrix("definitely_not_a_matrix")
+
+
+sizes = [10, 100, 1000]
+densities = [0.1, 0.3]
+
 
 @pytest.mark.parametrize("size", sizes)
 @pytest.mark.parametrize("density", densities)
 def test_make_symmetric(size, density):
     matrix = sp.random(size, size, density=density, format="csr")
     symmetric_matrix = make_symmetric(matrix)
-    
+
     assert np.allclose(symmetric_matrix.toarray(), symmetric_matrix.toarray().T)
     assert symmetric_matrix.shape == matrix.shape
 
     with pytest.raises(TypeError):
         _ = make_symmetric("banana")
 
 
-
 @pytest.mark.parametrize("size", sizes)
 @pytest.mark.parametrize("density", densities)
-def test_power_method(size, density):
+def test_implementations_power_method(size, density):
     matrix = sp.random(size, size, density=density, format="csr")
     matrix = make_symmetric(matrix)
 
@@ -40,14 +55,8 @@ def test_power_method(size, density):
     biggest_eigenvalue_sp = eigs_sp[index_sp]
 
     biggest_eigenvalue_pm = power_method(matrix)
+    biggest_eigenvalue_pm_numba = power_method_numba(matrix.toarray())
 
     assert np.isclose(biggest_eigenvalue_np, biggest_eigenvalue_sp, rtol=1e-4)
     assert np.isclose(biggest_eigenvalue_pm, biggest_eigenvalue_sp, rtol=1e-4)
-
-    ugly_nonsquare_matrix = np.random.rand(5, 3)
-
-    with pytest.raises(ValueError):
-        _ = eigenvalues_np(eigenvalues_np(ugly_nonsquare_matrix, symmetric=True))
-
-    with pytest.raises(ValueError):
-        _ = eigenvalues_sp(eigenvalues_sp(ugly_nonsquare_matrix, symmetric=True))
+    assert np.isclose(biggest_eigenvalue_pm_numba, biggest_eigenvalue_sp, rtol=1e-4)