boost of phi_2d calculations

epileptor/phi_2d.py

@@ -1,4 +1,5 @@
 import numpy as np
+from epileptor.phi_2d_iterative import calculate
 
 
 class PhiSolver:
@@ -11,14 +12,17 @@ class PhiSolver:
         return sigma
 
     def __init__(self, x_nh, y_nh, sigma):
-        self.phi_prev = None
+        self.phi_prev = np.zeros((y_nh, x_nh))
+        # make on item nonzero so we don't divide on zero further
+        self.phi_prev[1, 1] = 1
         self.x_nh = x_nh
         self.y_nh = y_nh
         self.sigma = sigma
         self.A = self.inverse_matrix_coef_explicit()
 
     def ode_step(self, nu):
-        self.phi_prev, _ = self.iterative_step(nu, 1e-8, self.phi_prev)
+        # self.phi_prev, _ = self.iterative_step(nu, 1e-8, self.phi_prev)
+        self.phi_prev, _ = calculate(self.sigma, nu, 1e-8, self.phi_prev)
         return self.phi_prev
         # return self.inverse_step(nu)
 
@@ -304,7 +308,7 @@ if __name__ == "__main__":
     import matplotlib.pyplot as plt
     import timeit
 
-    timeit_num = 1000
+    timeit_num = 10
     x_nh = 40
     y_nh = 40
     y_h = x_h = 3
@@ -326,21 +330,27 @@ if __name__ == "__main__":
 
     # inverse
     phi_inverse = solver.inverse_step(nu)
-    print('inverse time: {}'.format(timeit.timeit('solver.inverse_step(nu)', number=timeit_num, globals=globals())))
+    print('inverse time: {}'.format(
+        timeit.timeit('solver.inverse_step(nu)', number=timeit_num, globals=globals())
+    ))
 
     # iterative
     phi_iterative, _ = solver.iterative_step(nu, 1e-8)
-    print('iterative time: {}'.format(timeit.timeit('solver.iterative_step(nu, 1e-8)', number=timeit_num, globals=globals())))
+    print('iterative time: {}'.format(
+        timeit.timeit('solver.iterative_step(nu, 1e-8)', number=timeit_num, globals=globals())
+    ))
 
-    # iterative inverse
-    phi_iterative_inverse, _ = solver.iterative_inverse_step(nu, 1e-6)
-    # print('iterative inverse time: {}'.format(timeit.timeit()))
+    # iterative plain
+    phi_iterative_plain, _ = calculate(sigma, nu, 1e-8, np.ones(nu.shape))
+    print('iterative plain time: {}'.format(
+        timeit.timeit('calculate(sigma, nu, 1e-8, np.ones(nu.shape))', number=timeit_num, globals=globals())
+    ))
 
-    print(np.allclose(phi_inverse, phi_iterative_inverse, atol=1e-3))
+    print(np.allclose(phi_inverse, phi_iterative_plain, atol=1e-3))
     # idx = 10
-    # print(np.hstack((phi_inverse[:, idx][:, np.newaxis], phi_iterative[:, idx][:, np.newaxis], phi_iterative_inverse[:, idx][:, np.newaxis])))
-    # print(np.hstack((phi_inverse[idx, :][:, np.newaxis], phi_iterative[idx, :][:, np.newaxis], phi_iterative_inverse[idx, :][:, np.newaxis])))
-    print(np.sum(phi_inverse), np.sum(phi_iterative), np.sum(phi_iterative_inverse))
+    # print(np.hstack((phi_inverse[:, idx][:, np.newaxis], phi_iterative[:, idx][:, np.newaxis], phi_iterative_plain[:, idx][:, np.newaxis])))
+    # print(np.hstack((phi_inverse[idx, :][:, np.newaxis], phi_iterative[idx, :][:, np.newaxis], phi_iterative_plain[idx, :][:, np.newaxis])))
+    print(np.sum(phi_inverse), np.sum(phi_iterative), np.sum(phi_iterative_plain))
 
     x = np.linspace(0, x_h, x_nh)
     y = np.linspace(0, y_h, y_nh)
@@ -355,15 +365,15 @@ if __name__ == "__main__":
     ax_phi1 = plt.subplot(222, projection='3d')
     ax_phi1.plot_surface(x, y, phi_inverse, cmap='plasma')
     # ax_phi1.axis([x_h / 2 - 2 * x_h, x_h / 2 + 2 * x_h, 0, y_h])
-    ax_phi1.set_title('Phi1')
+    ax_phi1.set_title('Phi inverse')
 
     ax_phi2 = plt.subplot(223, projection='3d')
     ax_phi2.plot_surface(x, y, phi_iterative, cmap='plasma')
     # ax_phi2.axis([x_h / 2 - 2 * x_h, x_h / 2 + 2 * x_h, 0, y_h])
-    ax_phi2.set_title('Phi2')
+    ax_phi2.set_title('Phi iterative numpy')
 
     ax_phi3 = plt.subplot(224, projection='3d')
-    ax_phi3.plot_surface(x, y, phi_iterative_inverse, cmap='plasma')
-    ax_phi3.set_title('Phi3')
+    ax_phi3.plot_surface(x, y, phi_iterative_plain, cmap='plasma')
+    ax_phi3.set_title('Phi iterative plain')
 
     plt.show()

epileptor/phi_2d_iterative.py

@@ -0,0 +1,58 @@
+# for boosting of phi_2d computation this module does phi_2d.iterative_step in raw cycles and applies JIT Numba.
+# Pythran gives the same acceleration as Numba hence it is not used as it requries additional compilation step.
+from numba import jit
+import numpy as np
+
+
+# pythran export calculate(float64, float64[][], float64, float64[][])
+@jit
+def calculate(sigma, nu, l2_target, phi_init):
+    l2_norm = 1
+    iterations = 0
+    b = nu * sigma
+    # initial guess
+    phi = phi_init.copy()
+    phi_prev = np.empty_like(phi)
+
+    m, n = nu.shape
+    sigma4 = 1. / (4. + sigma)
+
+    while l2_norm > l2_target:
+        for i in range(0, m):
+            for j in range(0, n):
+                phi_prev[i, j] = phi[i, j]
+        # omp parallel for schedule(guided)
+        for i in range(1, m - 1):
+            # i - rows
+            for j in range(1, n - 1):
+                # j - columns
+                phi[i, j] = sigma4 * (
+                        phi_prev[i - 1, j] + phi_prev[i + 1, j] + phi_prev[i, j - 1] + phi_prev[i, j + 1] +
+                        b[i, j]
+                )
+        for i in range(0, m):
+            phi[i, 0] = phi[i, 1]
+            phi[i, n - 1] = phi[i, n - 2]
+        for j in range(0, n):
+            phi[0, j] = phi[1, j]
+            phi[m - 1, j] = phi[m - 2, j]
+
+        phi[0, 0] = phi[1, 1]
+        phi[-1, 0] = phi[-2, 1]
+        phi[0, -1] = phi[1, -2]
+        phi[-1, -1] = phi[-2, -2]
+
+        l2_norm_num = 0
+        l2_norm_div = 0
+        # omp parallel for schedule(guided)
+        for i in range(0, m):
+            # i - rows
+            for j in range(0, n):
+                l2_norm_num += (phi[i, j] - phi_prev[i, j]) ** 2
+                l2_norm_div += phi_prev[i, j] ** 2
+        l2_norm_div = max(l2_norm_div, 1e-6)
+        l2_norm = l2_norm_num ** 0.5 / l2_norm_div ** 0.5
+
+        iterations += 1
+
+    return phi, iterations

epileptor/phi_2d_iterative.pyx

@@ -0,0 +1,57 @@
+# experimental cython module.
+# WARNING! Do not use it.
+
+import numpy as np
+cimport numpy as cnp
+
+def calculate(double sigma, cnp.ndarray[cnp.float64_t, ndim=2] nu, double l2_target, cnp.ndarray[cnp.float64_t, ndim=2] phi_init):
+    cdef double l2_norm = 1
+    cdef unsigned int iterations = 0
+    cdef double[:,:] b = np.multiply(nu, sigma)
+    cdef double[:,:] phi = np.empty_like(nu)
+    # phi[:,:] = phi_init[:,:]
+    cdef double[:,:] phi_prev = np.empty_like(phi)
+
+    cdef unsigned int i, j, m, n
+    cdef double l2_norm_num, l2_norm_div
+    m = len(nu)
+    n = len(nu[0])
+    cdef double sigma4 = 1. / (4. + sigma)
+
+    while l2_norm > l2_target:
+        for i in range(0, m):
+            for j in range(0, n):
+                phi_prev[i, j] = phi[i, j]
+        for i in range(1, m - 1):
+            # i - rows
+            for j in range(1, n - 1):
+                # j - columns
+                phi[i, j] = sigma4 * (
+                        phi_prev[i - 1, j] + phi_prev[i + 1, j] + phi_prev[i, j - 1] + phi_prev[i, j + 1] +
+                        b[i, j]
+                )
+        for i in range(0, m):
+            phi[i, 0] = phi[i, 1]
+            phi[i, n - 1] = phi[i, n - 2]
+        for j in range(0, n):
+            phi[0, j] = phi[1, j]
+            phi[m - 1, j] = phi[m - 2, j]
+
+        phi[0, 0] = phi[1, 1]
+        phi[-1, 0] = phi[-2, 1]
+        phi[0, -1] = phi[1, -2]
+        phi[-1, -1] = phi[-2, -2]
+
+        l2_norm_num = 0
+        l2_norm_div = 0
+        for i in range(0, m):
+            # i - rows
+            for j in range(0, n):
+                l2_norm_num += (phi[i, j] - phi_prev[i, j]) ** 2
+                l2_norm_div += phi_prev[i, j] ** 2
+        l2_norm_div = max(l2_norm_div, 1e-6)
+        l2_norm = l2_norm_num ** 0.5 / l2_norm_div ** 0.5
+
+        iterations += 1
+
+    return phi, iterations