MFE: more clean up

benchmarks/magnetic_field_extrapolation/python_cython/magnetic_field_extrapolation.pyx

@@ -38,10 +38,12 @@ cdef calcB(np.ndarray B_x0, T_t alpha=0.0,
     cdef T_t [:, :] C = 4.0 / (n-1.0)**2 * np.sum(np.sum((B_x0 * np.sin(np.pi/y_max * u * y[:,None])[:,:,None])[:,None] * np.sin(np.pi/z_max * u * z[:,None])[:,None],-1),-1)
     cdef T_t [:, :] l = np.pi**2 * ((u[:]*u[:] / y_max)[:,None] + (u*u / z_max))
     l[0,0] = 1.0
-    cdef T_t [:, :] r = np.empty((n,n),dtype=T)
+    r_np = np.empty((n,n),dtype=T)
+    cdef T_t [:, :] r = r_np
     for i in range(n):
         for j in range(n):
             r[i,j] = np.sqrt(l[i,j] - alpha**2)
+    cdef T_t [:, :, :] exprx = np.exp((-r_np * x[:, None, None]))
 
     # Calculating B
     cdef T_t [:, :, :] Bx = np.empty((n,n,n), dtype=T)
@@ -66,10 +68,9 @@ cdef calcB(np.ndarray B_x0, T_t alpha=0.0,
             for k in range(n):
                 for m in range(n):
                     for q in range(n):
-                        exprx = exp(-r[m, q] * x[k])
-                        Bx[k, i, j] += temp_x[m, q] * exprx
-                        By[k, i, j] += temp_y[m, q] * exprx
-                        Bz[k, i, j] += temp_z[m, q] * exprx
+                        Bx[k, i, j] += temp_x[m, q] * exprx[k, m, q]
+                        By[k, i, j] += temp_y[m, q] * exprx[k, m, q]
+                        Bz[k, i, j] += temp_z[m, q] * exprx[k, m, q]
     return (Bx, By, Bz)
 
 
@@ -90,13 +91,9 @@ def main():
 
     B.start()
     for _ in range(B.size[2]):
-        Rx, Ry, Rz = calcB(window(B_x0.copy()))
+        calcB(window(B_x0.copy()))
     B.stop()
     B.pprint()
 
-    if B.outputfn:
-        R = Rx+Ry+Rz
-        B.tofile(B.outputfn, {'res': R, 'res_x': Rx, 'res_y': Ry, 'res_z': Rz})
-
 if __name__ == '__main__':
     main()

benchmarks/magnetic_field_extrapolation/python_cython/readme.rst

@@ -1 +1,5 @@
 Magnetic Field Extrapolation
+
+Compile with (in this folder)::
+
+    python setup.py build_ext --inplace -v

benchmarks/magnetic_field_extrapolation/python_numba/magnetic_field_extrapolation.py

@@ -34,6 +34,7 @@ def calcB_loop(B_x0, alpha=0.0,
     l = pi**2 * ((u**2 / y_max)[:,None] + (u**2 / z_max))
     l[0,0] = 1.0
     r = np.sqrt(l - alpha**2)
+    exprx = np.exp((-r * x[:, None, None]))
 
     # Calculating B
     Bx = np.empty((n,n,n),dtype=B_x0.dtype)
@@ -57,16 +58,15 @@ def calcB_loop(B_x0, alpha=0.0,
             for k in range(n):
                 for m in range(n):
                     for q in range(n):
-                        exprx = np.exp(-r[m, q] * x[k])
-                        Bx[k, i, j] += temp_x[m, q] * exprx
-                        By[k, i, j] += temp_y[m, q] * exprx
-                        Bz[k, i, j] += temp_z[m, q] * exprx
+                        Bx[k, i, j] += temp_x[m, q] * exprx[k, m, q]
+                        By[k, i, j] += temp_y[m, q] * exprx[k, m, q]
+                        Bz[k, i, j] += temp_z[m, q] * exprx[k, m, q]
     return (Bx, By, Bz)
 
 
 #@cuda.jit()
 @numba.jit(nopython=True)
-def loop(n, Bx, By, Bz, u, y, y_max, z, z_max, alpha, C, l, r, x, temp_x, temp_y, temp_z):
+def loop(n, Bx, By, Bz, u, y, y_max, z, z_max, alpha, C, l, r, temp_x, temp_y, temp_z, exprx):
     for i in range(n):
         for j in range(n):
             for k in range(n):
@@ -82,10 +82,9 @@ def loop(n, Bx, By, Bz, u, y, y_max, z, z_max, alpha, C, l, r, x, temp_x, temp_y
             for k in range(n):
                 for m in range(n):
                     for q in range(n):
-                        exprx = exp(-r[m, q] * x[k])
-                        Bx[k, i, j] += temp_x[m, q] * exprx
-                        By[k, i, j] += temp_y[m, q] * exprx
-                        Bz[k, i, j] += temp_z[m, q] * exprx
+                        Bx[k, i, j] += temp_x[m, q] * exprx[k, m, q]
+                        By[k, i, j] += temp_y[m, q] * exprx[k, m, q]
+                        Bz[k, i, j] += temp_z[m, q] * exprx[k, m, q]
 
 
 def calcB_numba(B_x0, alpha=0.0,
@@ -105,6 +104,7 @@ def calcB_numba(B_x0, alpha=0.0,
     l = pi**2 * ((u**2 / y_max)[:,None] + (u**2 / z_max))
     l[0,0] = 1.0
     r = np.sqrt(l - alpha**2)
+    exprx = np.exp((-r * x[:, None, None]))
 
     # Calculating B
     Bx = np.empty((n,n,n),dtype=B_x0.dtype)
@@ -113,7 +113,7 @@ def calcB_numba(B_x0, alpha=0.0,
     temp_x = np.empty((n, n), dtype=B_x0.dtype)
     temp_y = np.empty((n, n), dtype=B_x0.dtype)
     temp_z = np.empty((n, n), dtype=B_x0.dtype)
-    loop(n, Bx, By, Bz, u, y, y_max, z, z_max, alpha, C, l, r, x, temp_x, temp_y, temp_z)
+    loop(n, Bx, By, Bz, u, y, y_max, z, z_max, alpha, C, l, r, temp_x, temp_y, temp_z, exprx)
     return (Bx, By, Bz)
 
 
@@ -134,10 +134,12 @@ def main():
     B.start()
     for _ in range(B.size[2]):
         Rx, Ry, Rz = calcB_numba(window(B_x0.copy()))
+        """
         Rx1, Ry1, Rz1 = calcB_loop(window(B_x0.copy()))
         assert np.allclose(Rx, Rx1)
         assert np.allclose(Ry, Ry1)
         assert np.allclose(Rz, Rz1)
+        """
     B.stop()
     B.pprint()
 

benchmarks/magnetic_field_extrapolation/python_numpy/magnetic_field_extrapolation.py

@@ -66,7 +66,7 @@ def main():
 
     B.start()
     for _ in range(B.size[2]):
-        Rx, Ry, Rz = calcB(window(B_x0))
+        calcB(window(B_x0))
     B.stop()
     B.pprint()