WIP: check whether new workflow is ok

Author: lorenzotomada

.github/workflows/tests.yml

@@ -5,31 +5,27 @@ on:
     branches:
       - main
       - final_project
+      - gpu_branch
 
   pull_request:
     branches:
       - main
       - final_project
+      - gpu_branch
 
 jobs:
-  docs:
+  test:
     name: Testing using pytest
     runs-on: ubuntu-latest
+    container:
+      image: cupy/cupy:v13.4.0
+
     steps:
     - uses: actions/checkout@v3
 
-    - name: Set up Python
-      uses: actions/setup-python@v3
-      with:
-        python-version: 3.10.16
-    
-    - name: Set up Python environment
-      run: |
-        python -m venv venv
-        source venv/bin/activate
-
-    - name: Install dependencies
+    - name: Install additional dependencies
       run: |
+        python3 -m venv venv
         source venv/bin/activate
         python -m pip install .
 

shell/submit.sbatch

@@ -6,7 +6,7 @@
 #SBATCH --gpus=1
 #SBATCH --ntasks-per-node=1
 #SBATCH --gpus-per-task=1
-#SBATCH --mem=32000
+#SBATCH --mem=2000
 #SBATCH --time=01:00:00
 #SBATCH --mail-user=ltomada@sissa.it
 #SBATCH --output=%x.o%j.%N
@@ -24,11 +24,13 @@ echo '------------------------------------------------------'
 #
 
 cd $SLURM_SUBMIT_DIR           
-export SLURM_NTASKS_PER_NODE=2  # due to Ulysses's bug
+export SLURM_NTASKS_PER_NODE=1  # due to Ulysses's bug
 
 module load cuda/12.1
 conda init
+source ~/.bashrc
 conda activate ~/miniconda3/envs/devtools_scicomp
-
+pip install cupy-cuda12x
+python -m pip freeze > requirements.txt
 # Run the script
-python scripts/run.py fit --config=experiments/config.yaml
+#python scripts/run.py fit --config=experiments/config.yaml

src/pyclassify/__init__.py

@@ -1,13 +1,15 @@
 __all__ = [
     "eigenvalues_np",
     "eigenvalues_sp",
+    "eigenvalues_cp",
     "power_method",
     "power_method_numba",
 ]
 
 from .eigenvalues import (
     eigenvalues_np,
     eigenvalues_sp,
+    eigenvalues_cp,
     power_method,
     power_method_numba,
-)
+)

src/pyclassify/eigenvalues.py

@@ -1,16 +1,21 @@
 import numpy as np
 import scipy.sparse as sp
+import scipy.linalg as spla
 from line_profiler import profile
 from numpy.linalg import eig, eigh
-from pyclassify.utils import check_A_square_matrix, power_method_numba_helper
+from pyclassify.utils import check_square_matrix, check_symm_square, power_method_numba_helper
+import cupy as cp
+cp.cuda.Device(0)
+import cupy.linalg as cpla
+import cupyx.scipy.sparse as cpsp
 
 
 @profile
 def eigenvalues_np(A, symmetric=True):
     """
     Compute the eigenvalues of a square matrix using NumPy's `eig` or `eigh` function.
 
-    This function checks if the input matrix is square (and is actually a matrix) using 'check_A_square_matrix', and then computes its eigenvalues.
+    This function checks if the input matrix is square (and is actually a matrix) using 'check_square_matrix', and then computes its eigenvalues.
     If the matrix is symmetric, it uses `eigh` (which is more efficient for symmetric matrices).
     Otherwise, it uses `eig`.
 
@@ -27,7 +32,7 @@ def eigenvalues_np(A, symmetric=True):
         TypeError: If the input is not a NumPy array or a SciPy sparse matrix.
         ValueError: If number of rows != number of columns.
     """
-    check_A_square_matrix(A)
+    check_square_matrix(A)
     eigenvalues, _ = eigh(A) if symmetric else eig(A)
     return eigenvalues
 
@@ -53,22 +58,46 @@ def eigenvalues_sp(A, symmetric=True):
         TypeError: If the input is not a NumPy array or a SciPy sparse matrix.
         ValueError: If number of rows != number of columns.
     """
-    check_A_square_matrix(A)
+    check_square_matrix(A)
     eigenvalues, _ = (
-        sp.linalg.eigsh(A, k=A.shape[0] - 1)
+        spla.eigsh(A, k=A.shape[0] - 1)
         if symmetric
-        else sp.linalg.eigs(A, k=A.shape[0] - 1)
+        else spla.eigs(A, k=A.shape[0] - 1)
     )
     return eigenvalues
 
 
+@profile
+def eigenvalues_cp(A):
+    """
+    Compute the eigenvalues of a sparse matrix using CuPy's `eigsh` function.
+
+    This function checks if the input matrix is square and symmetric, then computes its eigenvalues using
+    CupY's sparse linear algebra solvers. For symmetric matrices, it uses `eigsh` for
+    more efficient computation.
+
+    Args:
+        A (cpsp.spmatrix): A square sparse matrix whose eigenvalues are to be computed.
+
+    Returns:
+        np.ndarray: An array containing the eigenvalues of the sparse matrix `A`.
+
+    Raises:
+        TypeError: If the input is not a CuPy sparse symmetric matrix.
+        ValueError: If number of rows != number of columns.
+    """
+    check_symm_square(A)
+    eigenvalues, _ = cpla.eigsh(A, k=A.shape[0] - 1)
+    return eigenvalues
+
+
 @profile
 def power_method(A, max_iter=500, tol=1e-4, x=None):
     """
     Compute the dominant eigenvalue of a square matrix using the power method.
 
     Args:
-        A (np.ndarray or sp.spmatrix): A square matrix whose dominant eigenvalue is to be computed.
+        A (np.ndarray or sp.spmatrix or cpsp.spmatrix): A square matrix whose dominant eigenvalue is to be computed.
         max_iter (int, optional): Maximum number of iterations to perform (default is 500).
         tol (float, optional): Tolerance for convergence based on the relative change between iterations
                                (default is 1e-4).
@@ -81,19 +110,19 @@ def power_method(A, max_iter=500, tol=1e-4, x=None):
         TypeError: If the input is not a NumPy array or a SciPy sparse matrix.
         ValueError: If number of rows != number of columns.
     """
-    check_A_square_matrix(A)
+    check_square_matrix(A)
     if x is None:
         x = np.random.rand(A.shape[0])
-    x /= np.linalg.norm(x)
+    x /= spla.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 /= spla.norm(x)
+        update_norm = spla.norm(x - x_old) / spla.norm(x_old)
         x_old = x.copy()
         iteration += 1
 
@@ -121,6 +150,44 @@ def power_method_numba(A):
     return power_method_numba_helper(A)
 
 
+@profile
+def power_method_cp(A, max_iter=500, tol=1e-4, x=None):
+    """
+    Compute the dominant eigenvalue of a square matrix using the power method.
+
+    Args:
+        A (cp.spmatrix): A square matrix whose dominant eigenvalue is to be computed.
+        max_iter (int, optional): Maximum number of iterations to perform (default is 500).
+        tol (float, optional): Tolerance for convergence based on the relative change between iterations
+                               (default is 1e-4).
+        x (cp.ndarray, optional): Initial guess for the eigenvector. If None, a random vector is generated.
+
+    Returns:
+        float: The approximated dominant eigenvalue of the matrix `A`, computed as the Rayleigh quotient x @ A @ x.
+
+    Raises:
+        TypeError: If the input is not a NumPy array or a SciPy sparse matrix.
+        ValueError: If number of rows != number of columns.
+    """
+    check_square_matrix(A)
+    if x is None:
+        x = cp.random.rand(A.shape[0])
+    x /= cpla.norm(x)
+    x_old = x
+
+    iteration = 0
+    update_norm = tol + 1
+
+    while iteration < max_iter and update_norm > tol:
+        x = A @ x
+        x /= cpla.norm(x)
+        update_norm = cpla.norm(x - x_old) / cpla.norm(x_old)
+        x_old = x.copy()
+        iteration += 1
+
+    return x @ A @ x
+
+
 def Lanczos_PRO(A, q, m=None, toll=np.sqrt(np.finfo(float).eps)):
     """
     Perform the Lanczos algorithm for symmetric matrices.
@@ -223,6 +290,5 @@ def QR_method(A_copy, tol=1e-10, max_iter=100):
             Q = Q@R
         A=A@Q
         iter+=1
-
 
     return np.diag(A), Q

src/pyclassify/utils.py

@@ -1,29 +1,32 @@
 import numpy as np
 import scipy.sparse as sp
+import cupy as cp
+import cupyx.scipy.sparse as cpsp
 import numba
 import os
 import yaml
 from line_profiler import profile
 
+
 # from numba.pycc import CC
 
 
-def check_A_square_matrix(A):
+def check_square_matrix(A):
     """
     Checks if the input matrix is a square matrix of type NumPy ndarray or SciPy sparse matrix.
     This is done to ensure that the input matrix `A` is both:
-    1. Of type `np.ndarray` (NumPy array) or `scipy.sparse` (SciPy sparse matrix).
+    1. Of type `np.ndarray` (NumPy array) or `scipy.sparse.spmatrix` (SciPy sparse matrix) or 'cupyx.scipy.sparse.spmatrix'.
     2. A square matrix.
 
     Args:
-        A (np.ndarray or sp.spmatrix): The matrix to be checked.
+        A (np.ndarray or sp.spmatrix or cpsp.spmatrix): The matrix to be checked.
 
     Raises:
         TypeError: If the input is not a NumPy array or a SciPy sparse matrix.
         ValueError: If number of rows != number of columns.
     """
-    if not isinstance(A, (np.ndarray, sp.spmatrix)):
-        raise TypeError("Input matrix must be a NumPy array or a SciPy sparse matrix!")
+    if not isinstance(A, (np.ndarray, sp.spmatrix, cpsp.spmatrix)):
+        raise TypeError("Input matrix must be a NumPy array or a SciPy/CuPy sparse matrix!")
     if A.shape[0] != A.shape[1]:
         raise ValueError("Matrix must be square!")
 
@@ -33,24 +36,46 @@ def make_symmetric(A):
     """
     Ensures the input matrix is symmetric by averaging it with its transpose.
 
-    This function first checks if the matrix is square using the `check_A_square_matrix` function.
+    This function first checks if the matrix is square using the `check_square_matrix` function.
     Then, it makes the matrix symmetric by averaging it with its transpose.
 
     Args:
-        A (np.ndarray or sp.spmatrix): The input square matrix to be made symmetric.
+        A (np.ndarray or sp.spmatrix or cpsp.spmatrix): The input square matrix to be made symmetric.
 
     Returns:
-        np.ndarray or sp.spmatrix: The symmetric version of the input matrix.
+        np.ndarray or sp.spmatrix or cpsp.spmatrix: The symmetric version of the input matrix.
 
     Raises:
-        TypeError: If the input matrix is not a NumPy array or SciPy sparse matrix.
+        TypeError: If the input matrix is not a NumPy array or SciPy or CuPy sparse matrix.
         ValueError: If the input matrix is not square.
     """
-    check_A_square_matrix(A)
+    check_square_matrix(A)
     A_sym = (A + A.T) / 2
     return A_sym
 
 
+def check_symm_square(A):
+    """
+    Checks if the input matrix is a square symmetric matrix of type SciPy/CuPy sparse matrix.
+    This is done to ensure that the input matrix `A` is all of the following:
+    1. A scipy sparse matrix or CuPy sparse matrix.
+    2. A square matrix.
+    3. Symmetric.
+
+    Args:
+        A (sp.spmatrix or cpsp.spmatrix): The matrix to be checked.
+
+    Raises:
+        TypeError: If the input is not a SciPy or CuPy sparse matrix.
+        ValueError: If number of rows != number of columns or the matrix is not symmetric.
+    """
+    check_square_matrix(A)
+    if isinstance(A, sp.spmatrix) and not np.allclose(A.toarray(), A.toarray().T):
+        raise ValueError("Matrix must be symmetric!")
+    elif isinstance(A, cpsp.spmatrix) and not cp.allclose(A.get(), A.get().T):
+        raise ValueError("Matrix must be symmetric!")
+
+
 @numba.njit(nogil=True, parallel=True)
 def power_method_numba_helper(A, max_iter=500, tol=1e-4, x=None):
     """
@@ -109,4 +134,4 @@ def read_config(file: str) -> dict:
     filepath = os.path.abspath(f"{file}.yaml")
     with open(filepath, "r") as stream:
         kwargs = yaml.safe_load(stream)
-    return kwargs
+    return kwargs