Performance improvement attempt

ELViM.py

@@ -11,11 +11,12 @@
 from numba import njit, prange
 import os
 import shutil
+import time
 import mdtraj as md
 import force_scheme as force
 import sys
 
-@njit(parallel=True, fastmath=False)         ####dissimilarity defautl, generalizar para diferentes sigma_0 e epsilon
+@njit(parallel=True, fastmath=False, cache=True)         ####dissimilarity defautl, generalizar para diferentes sigma_0 e epsilon
 def square_to_condensed(sqdmat):
    """
    Convert the symmetrical square matrix to the required condensed form.
@@ -24,20 +25,20 @@ def square_to_condensed(sqdmat):
    total = size*(size+1)//2
    d=np.zeros(total, dtype=np.float32) 
    for k in prange (size):
-      for l in prange (k, size):
+      for l in range (k, size):
          d[int(total - ((size - k) * (size - k + 1) / 2) + (l - k))] = sqdmat[k,l]
    return d 
 
 
-@njit(parallel=True, fastmath=False)         ####dissimilarity defautl, generalizar para diferentes sigma_0 e epsilon
+@njit(parallel=True, fastmath=False, cache=True)         ####dissimilarity defautl, generalizar para diferentes sigma_0 e epsilon
 def dissimilarity(coords, d, k, s0,size, atoms, total):
     epsilon=0.15
     sigma_0=s0
     for l in prange (k, size):
        N=0
        ep=0
-       for i in prange(atoms):
-          for j in prange (i+1,atoms):
+       for i in range(atoms):
+          for j in range (i+1,atoms):
              xijk = (coords[k][i][0] - coords[k][j][0])
              yijk = (coords[k][i][1] - coords[k][j][1])
              zijk = (coords[k][i][2] - coords[k][j][2])
@@ -53,6 +54,31 @@ def dissimilarity(coords, d, k, s0,size, atoms, total):
        d[int(total - ((size - k) * (size - k + 1) / 2) + (l - k))] = 1-ep/N
     pass   
 
+@njit(parallel=True, fastmath=False, cache=True)         ####dissimilarity defautl, generalizar para diferentes sigma_0 e epsilon
+def dissimilarity_no_verbose(coords, d, s0,size, atoms, total):
+    epsilon=0.15
+    sigma_0=s0
+    for k in prange (size):
+        for l in range (k, size):
+            N=0
+            ep=0
+            for i in range(atoms):
+                for j in range (i+1,atoms):
+                    xijk = (coords[k][i][0] - coords[k][j][0])
+                    yijk = (coords[k][i][1] - coords[k][j][1])
+                    zijk = (coords[k][i][2] - coords[k][j][2])
+                    rijk = np.sqrt(xijk*xijk + yijk*yijk + zijk*zijk)
+                    xijl = (coords[l][i][0] - coords[l][j][0])
+                    yijl = (coords[l][i][1] - coords[l][j][1])
+                    zijl = (coords[l][i][2] - coords[l][j][2])
+                    rijl = np.sqrt(xijl*xijl + yijl*yijl + zijl*zijl)
+                    dr = (rijk-rijl)**2
+                    sigma = sigma_0*abs(i-j)**epsilon 
+                    ep = ep + np.exp(-dr/(2*sigma**2))
+                    N = N + 1
+            d[int(total - ((size - k) * (size - k + 1) / 2) + (l - k))] = 1-ep/N
+    return d
+
 def calc_dmat(coords, s0, outdm, verbose):
     """
       Calculates the dissimilarity matrix using coordinates of C-alpha carbons.
@@ -83,9 +109,8 @@ def calc_dmat(coords, s0, outdm, verbose):
             dissimilarity(coords, dmat, k, s0, size, atoms, total)
 
     else:
-        for k in range (size):
-            dissimilarity(coords, dmat, k, s0, size, atoms, total)
-
+        dmat = dissimilarity_no_verbose(coords, dmat, s0, size, atoms, total)
+
     if outdm!=None:
         if os.path.exists(outdm):
             backup_file(outdm)

force_scheme.py

@@ -1,18 +1,18 @@
-import numpy as np
+import time
 import math
+import numpy as np
 from numba import njit, prange
-import mdtraj as md
 
 
 
-@njit(parallel=True)
+@njit(parallel=True, cache=True)
 def stress(distance_matrix, projection):
     size = len(projection)
     total = len(distance_matrix)
     den = 0
     num = 0
     for i in prange(size):
-        for j in prange(size):
+        for j in range (size):
             dr2 = math.sqrt((projection[i][0] - projection[j][0]) * (projection[i][0] - projection[j][0]) +
                             (projection[i][1] - projection[j][1]) * (projection[i][1] - projection[j][1]))
 
@@ -24,7 +24,7 @@ def stress(distance_matrix, projection):
     return math.sqrt(num / den)
 
 
-@njit(parallel=True, fastmath=False)
+@njit(parallel=True, fastmath=False, cache=True)
 def move(ins1, distance_matrix, projection, learning_rate):
     size = len(projection)
     total = len(distance_matrix)
@@ -61,6 +61,39 @@ def iteration(index, distance_matrix, projection, learning_rate):
         error += move(ins1, distance_matrix, projection, learning_rate)
     return error / size 
 
+@njit(parallel=False, fastmath=False, cache=True)
+def iteration_no_verbose(index, distance_matrix, projection, learning_rate):
+    size = len(projection)
+    error = 0
+
+    for i in prange(size):
+        ins1 = index[i]
+        size = len(projection)
+        total = len(distance_matrix)
+        error_l = 0
+        for ins2 in range(size):
+            if ins1 != ins2:
+                x1x2 = projection[ins2][0] - projection[ins1][0]
+                y1y2 = projection[ins2][1] - projection[ins1][1]
+                dr2 = max(math.sqrt(x1x2 * x1x2 + y1y2 * y1y2), 0.0001)
+
+                # getting te index in the distance matrix and getting the value
+                r = (ins1 + ins2 - math.fabs(ins1 - ins2)) / 2  # min(ins1,ins2)
+                s = (ins1 + ins2 + math.fabs(ins1 - ins2)) / 2  # max(ins1,ins2)
+                drn = distance_matrix[int(total - ((size - r) * (size - r + 1) / 2) + (s - r))]
+
+                # calculate the movement
+                delta = (drn - dr2) #* math.fabs(drn - dr2)
+
+                error_l += math.fabs(delta)
+
+                # moving
+                projection[ins2][0] += learning_rate * delta * (x1x2 / dr2)
+                projection[ins2][1] += learning_rate * delta * (y1y2 / dr2)
+        error += error_l / (size-1)
+
+    return error / size 
+
 #@njit(fastmath=False)
 def execute(distance_matrix, projection, max_it, verbose, learning_rate0=0.5, lrmin= 0.05, decay=0.95, tolerance=0.0000001):
     nr_iterations = 0
@@ -84,14 +117,14 @@ def execute(distance_matrix, projection, max_it, verbose, learning_rate0=0.5, lr
     else:
         for k in range(max_it):
             learning_rate = max(learning_rate0 * math.pow((1 - k / max_it), decay), lrmin)
-            new_error = iteration(index, distance_matrix, projection, learning_rate)
+            new_error = iteration_no_verbose(index, distance_matrix, projection, learning_rate)
             if math.fabs(new_error - error) < tolerance:
                 break
             error = new_error
             p_error[k] = error
             kstress[k] = stress(distance_matrix, projection)
             nr_iterations = k + 1  
-    
+
     # setting the min to (0,0)
     min_x = min(projection[:, 0])
     min_y = min(projection[:, 1])