Added draft for a numcore routine performing a faster real space forward calculation.

pyramid/magdata.py

@@ -5,7 +5,6 @@
 import numpy as np
 import tables.netcdf3 as nc
 import matplotlib.pyplot as plt
-from mayavi import mlab
 
 
 class MagData:
@@ -164,6 +163,8 @@ class MagData:
             None
 
         '''
+        from mayavi import mlab
+
         res = self.res
         dim = self.dim
         # Create points and vector components as lists:

pyramid/numcore/__init__.py

+from phase_mag_real import *

pyramid/numcore/phase_mag_real.pyx

+import numpy as np
+import math
+cimport cython
+cimport numpy as np
+
+
+@cython.boundscheck(False)
+@cython.wraparound(False)
+def phase_mag_real_helper_1(
+        unsigned int v_dim, unsigned int u_dim,
+        double[:, :] phi_u,
+        double[:, :] phi_v,
+        double[:, :] u_mag,
+        double[:, :] v_mag,
+        double[:, :] phase, float threshold):
+    cdef unsigned int i, j, ii, jj, iii, jjj
+    cdef double u, v
+    for j in range(v_dim):
+        for i in range(u_dim):
+            u = u_mag[j, i]
+            v = v_mag[j, i]
+            iii = u_dim - 1 - i
+            jjj = v_dim - 1 - j
+            if abs(u) > threshold:
+               for jj in range(phase.shape[0]):
+                    for ii in range(phase.shape[1]):
+                        phase[jj, ii] += u * phi_u[jjj + jj, iii + ii]
+            if abs(v) > threshold:
+               for jj in range(phase.shape[0]):
+                    for ii in range(phase.shape[1]):
+                        phase[jj, ii] -= v * phi_v[jjj + jj, iii + ii]
+

pyramid/phasemapper.py

@@ -3,7 +3,8 @@
 
 
 import numpy as np
-
+import numcore
+import time
 
 # Physical constants
 PHI_0 = -2067.83    # magnetic flux in T*nm²
@@ -102,12 +103,19 @@ def phase_mag_real(res, projection, method, b_0=1, jacobi=None):
                     phase -= v_mag[j, i] * phase_v
         ############################### TODO: NUMERICAL CORE  #####################################
     else:  # Without Jacobi matrix (faster)
+##        phasecopy = phase.copy()
+##        start_time = time.time()
+##        numcore.phase_mag_real_helper_1(v_dim, u_dim, phi_u, phi_v, u_mag, v_mag, phasecopy, threshold)
+##        print time.time() - start_time
+##        start_time = time.time()
         for j in range(v_dim):
             for i in range(u_dim):
                 if abs(u_mag[j, i]) > threshold:
                     phase += u_mag[j, i] * phi_u[v_dim-1-j:(2*v_dim-1)-j, u_dim-1-i:(2*u_dim-1)-i]
                 if abs(v_mag[j, i]) > threshold:
                     phase -= v_mag[j, i] * phi_v[v_dim-1-j:(2*v_dim-1)-j, u_dim-1-i:(2*u_dim-1)-i]
+##        print time.time() - start_time
+##        print ((phase - phasecopy) ** 2).sum()
     # Return the phase:
     return phase
 

setup.py

@@ -15,10 +15,10 @@ from Cython.Build import cythonize
 
 setup(
     name = 'Pyramid',
-    version = '0.1',    
+    version = '0.1',
     description = 'PYthon based Reconstruction Algorithm for MagnetIc Distributions',
     author = 'Jan Caron',
     author_email = 'j.caron@fz-juelich.de',
-    packages = ['pyramid'],
-    ext_modules = cythonize(glob.glob(os.path.join('pyramid','numcore','*.pyx')))
+    packages = ['pyramid', 'pyramid.numcore'],
+    ext_modules = cythonize(glob.glob(os.path.join('pyramid', 'numcore', '*.pyx')))
 )