Updating improvements

Author: lorenzotomada

scripts/profiling.py

@@ -0,0 +1,114 @@
+import cupyx.scipy.sparse as cpsp
+import cupy as cp
+from pyclassify import (
+    eigenvalues_np,
+    eigenvalues_sp,
+    eigenvalues_cp,
+    power_method,
+    power_method_numba,
+    power_method_cp,
+    QR,
+    QR_cp,
+)
+from pyclassify.utils import make_symmetric, read_config
+import numpy as np
+import scipy.sparse as sp
+import random
+import argparse
+import cProfile
+import os
+import csv
+
+
+seed = 8422
+random.seed(seed)
+cp.random.seed(seed)
+np.random.seed(seed)
+
+
+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["dim"]
+density = kwargs["density"]
+tol = kwargs["tol"]
+max_iter = kwargs["max_iter"]
+
+
+matrix = sp.random(dim, dim, density=density, format="coo")
+matrix = make_symmetric(matrix)
+cp_symm_matrix = cpsp.csr_matrix(matrix)
+
+
+logs = "./logs"
+iteration_factor = 300
+
+
+def profile_with_cprofile(func_name, func, *args, **kwargs):
+    log_file = "./logs/timings.csv"
+    with open(log_file, "a", newline="") as csvfile:
+        writer = csv.writer(csvfile)
+
+        profile_output = cProfile.Profile()
+        profile_output.enable()
+
+        result = func(*args, **kwargs)
+
+        profile_output.disable()
+
+        stats = profile_output.getstats()
+        total_time = sum(stat.totaltime for stat in stats)
+
+        print(f"{func_name}: {total_time:.6f} s")
+
+        writer.writerow([func_name, dim, total_time])
+
+    return result
+
+
+profile_with_cprofile(
+    "eigenvalues_np", eigenvalues_np, matrix.toarray(), symmetric=True
+)
+profile_with_cprofile("eigenvalues_sp", eigenvalues_sp, matrix, symmetric=True)
+profile_with_cprofile("power_method", power_method, matrix)
+profile_with_cprofile("power_method_numba", power_method_numba, matrix.toarray())
+profile_with_cprofile("QR", QR, matrix.toarray(), max_iter=iteration_factor * dim)
+
+
+def profile_with_cupy_profiler(func_name, func, *args, **kwargs):
+    log_file = "./logs/timings.csv"
+    with open(log_file, "a", newline="") as csvfile:
+        writer = csv.writer(csvfile)
+        start = cp.cuda.Event()
+        end = cp.cuda.Event()
+
+        start.record()
+        result = func(*args, **kwargs)
+        end.record()
+        end.synchronize()
+
+        result = func(*args, **kwargs)
+        elapsed_time = cp.cuda.get_elapsed_time(start, end) / 1000
+        print(f"{func_name}: {elapsed_time:.6f} s")
+
+        writer.writerow([func_name, dim, elapsed_time])
+
+    return result
+
+
+profile_with_cupy_profiler("eigenvalues_cp", eigenvalues_cp, cp_symm_matrix)
+profile_with_cupy_profiler("power_method_cp", power_method_cp, cp_symm_matrix)
+profile_with_cupy_profiler(
+    "QR_cp",
+    QR_cp,
+    cp_symm_matrix,
+    q0=cp.random.rand(dim),
+    tol=1e-3,
+    max_iter=iteration_factor * dim,
+)

scripts/run.py

@@ -7,15 +7,23 @@
     power_method,
     power_method_numba,
     power_method_cp,
+    QR,
+    QR_cp,
 )
 from pyclassify.utils import make_symmetric, read_config
 import numpy as np
 import scipy.sparse as sp
 import random
 import argparse
+import cProfile
+import cupyx.profiler as profiler
+import os
+import csv
 
 
 random.seed(226)
+cp.random.seed(226)
+np.random.seed(226)
 
 
 parser = argparse.ArgumentParser()
@@ -40,23 +48,45 @@
 cp_symm_matrix = cpsp.csr_matrix(matrix)
 
 
-eigs_np = eigenvalues_np(matrix.toarray(), symmetric=True)
-eigs_sp = eigenvalues_sp(matrix, symmetric=True)
-eigs_cp = eigenvalues_cp(cp_symm_matrix)
+logs = "./logs"
+iteration_factor = 300
 
 
-index_np = np.argmax(np.abs(eigs_np))
-index_sp = np.argmax(np.abs(eigs_sp))
-index_cp = np.argmax(np.abs(eigs_cp))
+cProfile.run(
+    "eigenvalues_np(matrix.toarray(), symmetric=True)",
+    os.path.join(logs, f"eigenvalues_np_{dim}.prof"),
+)
+cProfile.run(
+    "eigenvalues_np(matrix, symmetric=True)",
+    os.path.join(logs, f"eigenvalues_sp_{dim}.prof"),
+)
+cProfile.run("power_method(matrix)", os.path.join(logs, f"power_method_{dim}.prof"))
+cProfile.run(
+    "power_method(matrix.toarray())",
+    os.path.join(logs, f"power_method_numba{dim}.prof"),
+)
+cProfile.run(
+    "QR(matrix.toarray(), max_iter=iteration_factor*d)",
+    os.path.join(logs, f"QR_{dim}.prof"),
+)
+
+
+def profile_function(func_name, func, *args, **kwargs):
+    cupyx.profiler.start()
+    result = func(*args, **kwargs)
+    cupyx.profiler.stop()
+
+    profiler_data = cupyx.profiler.get_profiler()
 
+    output_file = os.path.join(logs, f"{func_name}_{dim}.prof")
+    with open(output_file, "w") as f:
+        f.write(str(profiler_data))
 
-biggest_eigenvalue_np = eigs_np[index_np]
-biggest_eigenvalue_sp = eigs_sp[index_sp]
-biggest_eigenvalue_cp = eigs_cp[index_cp]
+    elapsed_time = profiler_data[0]["time"]
+    print(f"{func_name}: {elapsed_time/1000} s")
+    return result
 
-biggest_eigenvalue_pm = power_method(matrix)
-# biggest_eigenvalue_pm_numba = power_method_numba(matrix.toarray())
-biggest_eigenvalue_cp = power_method_cp(cp_symm_matrix)
 
-_, __ = QR(A)
-_, __ = QR_cp(A)
+profile_function("eigenvalues_cp", eigenvalues_cp, cp_symm_matrix)
+profile_function("power_method_cp", power_method_cp, cp_symm_matrix)
+profile_function("QR_cp", QR_cp, cp_symm_matrix, max_iter=iteration_factor * d)

shell/submit.sbatch

@@ -31,34 +31,24 @@ conda init
 source ~/.bashrc
 conda activate /scratch/ltomada/miniconda3/envs/dtsc
 
-mkdir -p tmp_data experiments logs
-touch tmp_data/timings.csv
+#touch log/timings.csv
 
-export LINE_PROFILE=0
-
-density=0.1
+density=0.2
 tol=0.0001
 max_iter=1000
 
+echo "function,dim,time">logs/timings.csv
 
-for dim in 100 500 1000 5000; do
+for dim in 100 500 1000; do
     echo "d: $dim" # for readability
     echo "dim: $dim" > experiments/config.yaml    
     echo "density: $density" >> experiments/config.yaml
     echo "tol: $tol" >> experiments/config.yaml
     echo "max_iter: $max_iter" >> experiments/config.yaml
 
-    python -m kernprof -l -o tmp_data/profile.dat scripts/run.py --config=experiments/config
-
-    if [ ! -f tmp_data/timings.csv ]; then
-        echo "dim,function_name,time" > tmp_data/timings.csv
-    fi
-
-    python -m line_profiler -rmt tmp_data/profile.dat | awk -v dim="$dim" '
-        /Function: / {func_name=$2}
-        /Total time:/ {if (func_name != "") print dim, func_name, $3}' OFS=',' >> tmp_data/timings.csv
+    python scripts/profiling.py
 
-    echo "dim: $dim - times saved in tmp_data/timings.csv" 
+    echo "dim: $dim - times saved in logs/timings.csv" 
     echo '--------------------------------------------------'
 done
 

src/pyclassify/eigenvalues.py

@@ -196,7 +196,7 @@ def power_method_cp(A, max_iter=500, tol=1e-4, x=None):
     return x @ A @ x
 
 
-@jit(nopython=True)
+# @jit(nopython=True)
 def Lanczos_PRO(A, q=None, m=None, tol=1e-8):
     r"""
     Perform the Lanczos algorithm for symmetric matrices.
@@ -267,8 +267,8 @@ def Lanczos_PRO(A, q=None, m=None, tol=1e-8):
     return Q, np.array(alpha), np.array(beta[:-1])
 
 
-@jit(nopython=True)
-def QR_method(diag, off_diag, tol=1e-8, max_iter=100):
+# @jit(nopython=True)
+def QR_method(diag, off_diag, tol=1e-8, max_iter=1000):
     """
     Compute the eigenvalues of a tridiagonal matrix using the QR algorithm.
 
@@ -366,6 +366,8 @@ def QR_method(diag, off_diag, tol=1e-8, max_iter=100):
                         x_new = x_0 - np.cos(x_0) * np.cos(x_0) * (
                             np.tan(x_0) + b_2 / d
                         )
+                        if np.isclose(x_new, 0):
+                            break
                         err_rel = np.abs((x_new - x_0) / x_new)
                         x_0 = x_new
                         iter_newton += 1
@@ -394,7 +396,7 @@ def QR_method(diag, off_diag, tol=1e-8, max_iter=100):
 
 
 @profile
-def QR(A, q0=None, tol=1e-8, max_iter=100):
+def QR(A, q0=None, tol=1e-8, max_iter=1000):
     """
     Compute the eigenvalues of a square matrix using the QR algorithm.
     Done using the Lanczos algorithm to compute the tridiagonal matrix and then the QR
@@ -480,12 +482,11 @@ def Lanczos_PRO_cp(A, q=None, m=None, tol=1e-8):
         beta.append(cp.linalg.norm(r))
 
         if cp.abs(beta[j]) < 1e-15:
-
-            return cp.array(alpha), cp.array(beta[:-1])
+            return Q, cp.array(alpha), cp.array(beta[:-1])
     return Q, cp.array(alpha), cp.array(beta[:-1])
 
 
-def QR_method_cp(diag, off_diag, tol=1e-8, max_iter=100):
+def QR_method_cp(diag, off_diag, tol=1e-8, max_iter=1000):
     """
     Compute the eigenvalues of a tridiagonal matrix using the QR algorithm.
 
@@ -568,7 +569,7 @@ def QR_method_cp(diag, off_diag, tol=1e-8, max_iter=100):
                         s = -cp.sqrt(2) / 2
                     else:
                         c = s = cp.sqrt(2) / 2
-
+                # if off diagonal ... do nothing
                 else:
                     b_2 = off_diag[0]
                     if off_diag[0] * d > 0:
@@ -582,6 +583,8 @@ def QR_method_cp(diag, off_diag, tol=1e-8, max_iter=100):
                         x_new = x_0 - cp.cos(x_0) * cp.cos(x_0) * (
                             cp.tan(x_0) + b_2 / d
                         )
+                        if cp.isclose(cp.linalg.norm(x_new), 0):
+                            break
                         err_rel = cp.abs((x_new - x_0) / x_new)
                         x_0 = x_new
                         iter_newton += 1
@@ -610,7 +613,7 @@ def QR_method_cp(diag, off_diag, tol=1e-8, max_iter=100):
 
 
 @profile
-def QR_cp(A, q0=None, tol=1e-8, max_iter=100):
+def QR_cp(A, q0=None, tol=1e-8, max_iter=1000):
     """
     Compute the eigenvalues of a square matrix using the QR algorithm.
     Done using the Lanczos algorithm to compute the tridiagonal matrix and then the QR

test/test_.py

@@ -158,7 +158,7 @@ def test_QR_method(size):
 
 
 @pytest.mark.parametrize("size", sizes)
-def test_qr(size):
+def test_QR(size):
     eig = np.arange(1, size + 1)
     A = np.diag(eig)
     U = scipy.stats.ortho_group.rvs(size)