Added initcython.py (import this befor importing .pyx files to use Cython-Code)

scripts/get_jacobi.py

@@ -48,6 +48,8 @@ def phase_from_mag():
     np.savetxt('../output/jacobi.npy', jacobi)
     print 'Time for Real Space Approach with Jacobi-Matrix (Slab): ' + str(toc - tic)
     
+    
+    
     return jacobi
     
     

src/pyramex/SConstruct

@@ -6,11 +6,17 @@ PYTHON_LIBPATH = distutils.sysconfig.get_python_lib()
 PYTHON_INCPATH = distutils.sysconfig.get_python_inc()
 PYTHON_LIBRARY = "python" + sys.version[:3]
 
-cythonBuilder = Builder(action ='cython $SOURCE',
+CythonBuilder = Builder(action ='cython $SOURCE',
                         suffix = '.c',
                         src_suffix = '.pyx')
 
+def PydBuilder(env,source):
+    cCode = env.Py2C(source)
+    env.SharedLibrary(target='karl.pyd', source=cCode, LIBPREFIX='')
+
 env = Environment(ENV      =   os.environ,
-                  BUILDERS = {'Py2C':cythonBuilder})
+                  BUILDERS = {'Py2C':CythonBuilder})
+
+env.AddMethod(CythonBuilder,'Cython')
 
-env.Py2C('c1.c', 'c1.pyx')
\ No newline at end of file
+env.Cython('c1')
\ No newline at end of file

src/pyramex/initcython.py

+# -*- coding: utf-8 -*-
+"""Script to initialize Cythons pyximport to ensure compatibility with MinGW compiler and NumPy."""
+
+
+import os
+import numpy
+import pyximport
+
+
+if os.name == 'nt':
+    if os.environ.has_key('CPATH'):
+        os.environ['CPATH'] = os.environ['CPATH'] + numpy.get_include()
+    else:
+        os.environ['CPATH'] = numpy.get_include()
+#    # XXX: assuming that MinGW is installed in C:\MinGW (default)
+#    #      for PythonXY the default is C:\MinGW-xy
+#    if os.environ.has_key('PATH'):
+#        os.environ['PATH'] = os.environ['PATH'] + ';C:\MinGW\bin'
+#    else:
+#        os.environ['PATH'] = 'C:\MinGW\bin'
+
+    mingw_setup_args = { 'options': { 'build_ext': { 'compiler': 'mingw32' } } }
+    pyximport.install(setup_args=mingw_setup_args)
+
+
+elif os.name == 'posix':
+    if os.environ.has_key('CFLAGS'):
+        os.environ['CFLAGS'] = os.environ['CFLAGS'] + ' -I' + numpy.get_include()
+    else:
+        os.environ['CFLAGS'] = ' -I' + numpy.get_include()
+        
+    pyximport.install()
\ No newline at end of file

src/pyramex/setup.py

@@ -5,32 +5,12 @@ Created on Fri May 03 10:27:04 2013
 @author: Jan
 """
 
-#from distutils.core import setup, Extension
-#
-#module1 = Extension('c1',
-#                    sources = ['c1.pyx'])
-#
-#setup (name = 'pyramex',
-#       version = '1.0',
-#       description = 'This is a demo package',
-#       ext_modules = [module1])
-
-
+#call with: python setup.py build_ext --inplace --compiler=mingw32
 
 from distutils.core import setup
-from distutils.extension import Extension
-from Cython.Distutils import build_ext
-
-def get_extensions():
-	return [
-		Extension('hello', ['hello.pyx']),
-		Extension('c1',    ['c1.pyx']),
-		Extension('c2',    ['c2.pyx']),
-		Extension('c3',    ['c3.pyx'])
-		]
+from Cython.Build import cythonize
 
 setup(
-  name = 'pyramex',
-  cmdclass = {'build_ext': build_ext},
-  ext_modules = get_extensions()
+    name = "My hello app",
+    ext_modules = cythonize('hello.pyx'), # accepts a glob pattern
 )
\ No newline at end of file

src/pyramid/phasemap.py

@@ -112,34 +112,44 @@ def real_space(mag_data, method, b_0=1, jacobi=None):
     phi_sin = phi_pixel(method, yy_big, xx_big, res, b_0)
             
     def phi_mag(i, j):  # TODO: rename
-        return (np.cos(beta[j,i])*phi_cos[y_dim-1-j:(2*y_dim-1)-j, 
-                                          x_dim-1-i:(2*x_dim-1)-i]
-               -np.sin(beta[j,i])*phi_sin[y_dim-1-j:(2*y_dim-1)-j,
-                                          x_dim-1-i:(2*x_dim-1)-i])
+        return (np.cos(beta[j,i])*phi_cos[y_dim-1-j:(2*y_dim-1)-j, x_dim-1-i:(2*x_dim-1)-i]
+               -np.sin(beta[j,i])*phi_sin[y_dim-1-j:(2*y_dim-1)-j, x_dim-1-i:(2*x_dim-1)-i])
                                           
     def phi_mag_deriv(i, j):  # TODO: rename
-        return -(np.sin(beta[j,i])*phi_cos[y_dim-1-j:(2*y_dim-1)-j, 
-                                           x_dim-1-i:(2*x_dim-1)-i]
-                +np.cos(beta[j,i])*phi_sin[y_dim-1-j:(2*y_dim-1)-j,
-                                           x_dim-1-i:(2*x_dim-1)-i])   
+        return -(np.sin(beta[j,i])*phi_cos[y_dim-1-j:(2*y_dim-1)-j, x_dim-1-i:(2*x_dim-1)-i]
+                +np.cos(beta[j,i])*phi_sin[y_dim-1-j:(2*y_dim-1)-j, x_dim-1-i:(2*x_dim-1)-i])
+                                           
+    def phi_mag_fd(i, j, h):  # TODO: rename
+        return ((np.cos(beta[j,i]+h) - np.cos(beta[j,i])) / h 
+                      * phi_cos[y_dim-1-j:(2*y_dim-1)-j, x_dim-1-i:(2*x_dim-1)-i]
+               -(np.sin(beta[j,i]+h) - np.sin(beta[j,i])) / h 
+                      * phi_sin[y_dim-1-j:(2*y_dim-1)-j, x_dim-1-i:(2*x_dim-1)-i])
     
     '''CALCULATE THE PHASE'''
     phase = np.zeros((y_dim, x_dim))
     
     # TODO: only iterate over pixels that have a magn. > threshold (first >0)
+    jacobi_fd = jacobi.copy()
+    h = 0.0001
+    
     for j in range(y_dim):
         for i in range(x_dim):
-            if (mag[j, i] != 0 ):#or jacobi is not None): # TODO: same result with or without?
+            #if (mag[j,i] != 0 ):#or jacobi is not None): # TODO: same result with or without?
                 phi_mag_cache = phi_mag(i, j)
                 phase += mag[j,i] * phi_mag_cache
                 if jacobi is not None:
                     jacobi[:,i+x_dim*j] = phi_mag_cache.reshape(-1)
-                    jacobi[:,x_dim*y_dim+i+x_dim*j] = (mag[j,i]*phi_mag_deriv(i, j)).reshape(-1)         
+                    jacobi[:,x_dim*y_dim+i+x_dim*j] = (mag[j,i]*phi_mag_deriv(i,j)).reshape(-1)
+                    
+                    jacobi_fd[:,i+x_dim*j] = phi_mag_cache.reshape(-1)
+                    jacobi_fd[:,x_dim*y_dim+i+x_dim*j] = (mag[j,i]*phi_mag_fd(i,j,h)).reshape(-1)  
+                    
+                    
     
-    return phase
-
+    jacobi_diff = jacobi_fd - jacobi
+    assert (np.abs(jacobi_diff) < 1.0E-8).all(), 'jacobi matrix is not the same'
     
-
+    return phase
     
 
 def phase_elec(mag_data, b_0=1, v_0=0, v_acc=30000):