phasemapper: Cleaned up the two approaches (mx and my, m and beta)

pyramid/phasemapper.py

@@ -62,6 +62,7 @@ def phase_mag_real(res, projection, method, b_0=1, jacobi=None):
         the phasemap as a 2 dimensional array
     
     '''
+    # Function for creating the lookup-tables:
     def phi_lookup(method, n, m, res, b_0):
         if method == 'slab':    
             def F_h(n, m):
@@ -73,13 +74,6 @@ def phase_mag_real(res, projection, method, b_0=1, jacobi=None):
         elif method == 'disc':
             in_or_out = np.logical_not(np.logical_and(n == 0, m == 0))
             return coeff * m / (n**2 + m**2 + 1E-30) * in_or_out
-#    def phi_mag(i, j):  # TODO: rename
-#        return (np.cos(beta[j,i])*phi_u[v_dim-1-j:(2*v_dim-1)-j, u_dim-1-i:(2*u_dim-1)-i]
-#               -np.sin(beta[j,i])*phi_v[v_dim-1-j:(2*v_dim-1)-j, u_dim-1-i:(2*u_dim-1)-i])
-#                                          
-#    def phi_mag_deriv(i, j):  # TODO: rename
-#        return -(np.sin(beta[j,i])*phi_u[v_dim-1-j:(2*v_dim-1)-j, u_dim-1-i:(2*u_dim-1)-i]
-#                +np.cos(beta[j,i])*phi_v[v_dim-1-j:(2*v_dim-1)-j, u_dim-1-i:(2*u_dim-1)-i])
     # Process input parameters:
     v_dim, u_dim = np.shape(projection[0])
     v_mag, u_mag = projection 
@@ -90,22 +84,35 @@ def phase_mag_real(res, projection, method, b_0=1, jacobi=None):
     uu, vv = np.meshgrid(u, v)
     phi_u = phi_lookup(method, uu, vv, res, b_0)
     phi_v = phi_lookup(method, vv, uu, res, b_0)    
-    '''CALCULATE THE PHASE'''
+    # Calculation of the phase:
     phase = np.zeros((v_dim, u_dim))
-    for j in range(v_dim):  # TODO: only iterate over pixels that have a magn. > threshold (first >0)
-        for i in range(u_dim):
-            if (u_mag[j,i] != 0 and v_mag[j,i] != 0) or jacobi is not None: # TODO: same result with or without?
-                phase_u = phi_u[v_dim-1-j:(2*v_dim-1)-j, u_dim-1-i:(2*u_dim-1)-i]                
+    threshold = 0
+    if jacobi is not None:  # With Jacobian matrix (slower)
+        jacobi[:] = 0  # Jacobi matrix --> zeros
+        ############################### TODO: NUMERICAL CORE  #####################################
+        for j in range(v_dim):
+            for i in range(u_dim):
+                phase_u = phi_u[v_dim-1-j:(2*v_dim-1)-j, u_dim-1-i:(2*u_dim-1)-i]
+                jacobi[:,i+u_dim*j] = phase_u.reshape(-1)
+                if abs(u_mag[j, i]) > threshold:
+                    phase += u_mag[j,i] * phase_u
                 phase_v = phi_v[v_dim-1-j:(2*v_dim-1)-j, u_dim-1-i:(2*u_dim-1)-i]
-                phase += u_mag[j,i]*phase_u - v_mag[j,i]*phase_v
-                if jacobi is not None:
-                    jacobi[:,i+u_dim*j]             =  phase_u.reshape(-1)
-                    jacobi[:,u_dim*v_dim+i+u_dim*j] = -phase_v.reshape(-1)
-    
+                jacobi[:,u_dim*v_dim+i+u_dim*j] = -phase_v.reshape(-1)
+                if abs(v_mag[j, i]) > threshold:
+                    phase -= v_mag[j,i] * phase_v
+        ############################### TODO: NUMERICAL CORE  #####################################
+    else:  # Without Jacobi matrix (faster)
+        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]
+    # Return the phase:
     return phase
     
     
-def phase_mag_real_ANGLE(res, projection, method, b_0=1, jacobi=None):  # TODO: Modify
+def phase_mag_real_alt(res, projection, method, b_0=1, jacobi=None):  # TODO: Modify
     '''Calculate phasemap from magnetization data (real space approach).
     Arguments:
         res        - the resolution of the grid (grid spacing) in nm
@@ -117,78 +124,58 @@ def phase_mag_real_ANGLE(res, projection, method, b_0=1, jacobi=None):  # TODO:
         the phasemap as a 2 dimensional array
     
     '''
+    # Function for creating the lookup-tables:
     def phi_lookup(method, n, m, res, b_0):
         if method == 'slab':    
             def F_h(n, m):
                 a = np.log(res**2 * (n**2 + m**2))
                 b = np.arctan(n / m)
                 return n*a - 2*n + 2*m*b            
-            coeff = -b_0 * res**2 / ( 2 * PHI_0 )
             return coeff * 0.5 * ( F_h(n-0.5, m-0.5) - F_h(n+0.5, m-0.5) 
                                   -F_h(n-0.5, m+0.5) + F_h(n+0.5, m+0.5) ) 
         elif method == 'disc':
-            coeff = - b_0 * res**2 / ( 2 * PHI_0 )
             in_or_out = np.logical_not(np.logical_and(n == 0, m == 0))
             return coeff * m / (n**2 + m**2 + 1E-30) * in_or_out
-    
+    # Function for the phase contribution of one pixel:
+    def phi_mag(i, j):
+        return (np.cos(beta[j,i])*phi_u[v_dim-1-j:(2*v_dim-1)-j, u_dim-1-i:(2*u_dim-1)-i]
+               -np.sin(beta[j,i])*phi_v[v_dim-1-j:(2*v_dim-1)-j, u_dim-1-i:(2*u_dim-1)-i])
+    # Function for the derivative of the phase contribution of one pixel:                                   
+    def phi_mag_deriv(i, j):
+        return -(np.sin(beta[j,i])*phi_u[v_dim-1-j:(2*v_dim-1)-j, u_dim-1-i:(2*u_dim-1)-i]
+                +np.cos(beta[j,i])*phi_v[v_dim-1-j:(2*v_dim-1)-j, u_dim-1-i:(2*u_dim-1)-i])
+    # Process input parameters:
     v_dim, u_dim = np.shape(projection[0])
     v_mag, u_mag = projection
-    
     beta = np.arctan2(v_mag, u_mag)
-    mag = np.hypot(u_mag, v_mag) 
-                
-    '''CREATE COORDINATE GRIDS'''
-    u = np.linspace(0,(u_dim-1),num=u_dim)
-    v = np.linspace(0,(v_dim-1),num=v_dim)
-    uu, vv = np.meshgrid(u,v)
-     
-    u_lookup = np.linspace(-(u_dim-1), u_dim-1, num=2*u_dim-1)
-    v_lookup = np.linspace(-(v_dim-1), v_dim-1, num=2*v_dim-1)
-    uu_lookup, vv_lookup = np.meshgrid(u_lookup, v_lookup)
-    
-    phi_cos = phi_lookup(method, uu_lookup, vv_lookup, res, b_0)
-    phi_sin = phi_lookup(method, vv_lookup, uu_lookup, res, b_0)
-            
-    def phi_mag(i, j):  # TODO: rename
-        return (np.cos(beta[j,i])*phi_cos[v_dim-1-j:(2*v_dim-1)-j, u_dim-1-i:(2*u_dim-1)-i]
-               -np.sin(beta[j,i])*phi_sin[v_dim-1-j:(2*v_dim-1)-j, u_dim-1-i:(2*u_dim-1)-i])
-                                          
-    def phi_mag_deriv(i, j):  # TODO: rename
-        return -(np.sin(beta[j,i])*phi_cos[v_dim-1-j:(2*v_dim-1)-j, u_dim-1-i:(2*u_dim-1)-i]
-                +np.cos(beta[j,i])*phi_sin[v_dim-1-j:(2*v_dim-1)-j, u_dim-1-i:(2*u_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[v_dim-1-j:(2*v_dim-1)-j, u_dim-1-i:(2*u_dim-1)-i]
-               -(np.sin(beta[j,i]+h) - np.sin(beta[j,i])) / h 
-                      * phi_sin[v_dim-1-j:(2*v_dim-1)-j, u_dim-1-i:(2*u_dim-1)-i])
-    
-    '''CALCULATE THE PHASE'''
+    mag  = np.hypot(u_mag, v_mag) 
+    coeff = -b_0 * res**2 / ( 2 * PHI_0 )            
+    # Create lookup-tables for the phase of one pixel:
+    u = np.linspace(-(u_dim-1), u_dim-1, num=2*u_dim-1)
+    v = np.linspace(-(v_dim-1), v_dim-1, num=2*v_dim-1)
+    uu, vv = np.meshgrid(u, v)
+    phi_u = phi_lookup(method, uu, vv, res, b_0)
+    phi_v = phi_lookup(method, vv, uu, res, b_0)    
+    # Calculation of the phase:
     phase = np.zeros((v_dim, u_dim))
-    
-    # TODO: only iterate over pixels that have a magn. > threshold (first >0)
-    if jacobi is not None:
-        jacobi_fd = jacobi.copy()
-    h = 0.0001
-    
-    for j in range(v_dim):
-        for i in range(u_dim):
-            #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+u_dim*j] = phi_mag_cache.reshape(-1)
+    threshold = 0
+    if jacobi is not None:  # With Jacobian matrix (slower)
+        jacobi[:] = 0  # Jacobi matrix --> zeros
+        ############################### TODO: NUMERICAL CORE  #####################################
+        for j in range(v_dim):
+            for i in range(u_dim):
+                phase_cache = phi_mag(i, j)
+                jacobi[:,i+u_dim*j] = phase_cache.reshape(-1)
+                if mag[j, i] > threshold:
+                    phase += mag[j,i]*phase_cache
                     jacobi[:,u_dim*v_dim+i+u_dim*j] = (mag[j,i]*phi_mag_deriv(i,j)).reshape(-1)
-                    
-                    jacobi_fd[:,i+u_dim*j] = phi_mag_cache.reshape(-1)
-                    jacobi_fd[:,u_dim*v_dim+i+u_dim*j] = (mag[j,i]*phi_mag_fd(i,j,h)).reshape(-1)  
-                    
-                    
-    
-    if jacobi is not None:
-        jacobi_diff = jacobi_fd - jacobi
-        assert (np.abs(jacobi_diff) < 1.0E-8).all(), 'jacobi matrix is not the same'
-    
+        ############################### TODO: NUMERICAL CORE  #####################################
+    else:  # Without Jacobi matrix (faster)
+        for j in range(v_dim):
+            for i in range(u_dim):
+                if abs(mag[j, i]) > threshold:
+                    phase += mag[j,i] * phi_mag(i, j)
+    # Return the phase:
     return phase
     
 

pyramid/test/run_tests.py

@@ -31,7 +31,7 @@ def run():
     suite.addTest(loader.loadTestsFromTestCase(TestCaseAnalytic))
     suite.addTest(loader.loadTestsFromTestCase(TestCaseReconstructor))
     runner.run(suite)
-
+    
 
 if __name__ == '__main__':
     run()
\ No newline at end of file

pyramid/test/test_compliance.py

@@ -51,10 +51,15 @@ class TestCaseCompliance(unittest.TestCase):
         self.assertTrue(self.pep8 is not None,
                         msg="Install Python pep8 module to fully execute testbench!")
         errors = 0
-        for dir_name in ["test", os.environ["BUILDDIR"]]:
+        for dir_name in "./":
             errors += self.checkDirectory(dir_name)
         self.assertEqual(errors, 0)
 
+
+if __name__ == '__main__':
+    suite = unittest.TestLoader().loadTestsFromTestCase(TestCaseCompliance)
+    unittest.TextTestRunner(verbosity=2).run(suite)
+
 #    def test_variables(self):
 #        """
 #        Function for checking that all attributes are present.

scripts/get_jacobi.py

@@ -10,12 +10,13 @@ import pyramid.magcreator as mc
 import pyramid.projector as pj
 import pyramid.phasemapper as pm
 from pyramid.magdata import MagData
+from pyramid.phasemap import PhaseMap
 import time
 import pdb, traceback, sys
 from numpy import pi
 
 
-def phase_from_mag():
+def get_jacobi():
     '''Calculate and display the phase map from a given magnetization.
     Arguments:
         None
@@ -27,22 +28,21 @@ def phase_from_mag():
     b_0 = 1.0  # in T    
     dim = (1, 3, 3)  # in px (y,x)
     res = 10.0  # in nm
-    beta = 0*pi/4    
+    beta = pi/4    
     
     center = (0, 1, 1)  # in px (y,x) index starts with 0!
     width  = (0, 1, 1)  # in px (y,x)
-    mag_shape = mc.Shapes.slab(dim, center, width)
 
-    mag_data = MagData(res, mc.create_mag_dist(mag_shape, beta))
-    
+    mag_data = MagData(res, mc.create_mag_dist(mc.Shapes.slab(dim, center, width), beta))
     projection = pj.simple_axis_projection(mag_data)
     
     '''NUMERICAL SOLUTION'''
     # numerical solution Real Space (Slab):
     jacobi = np.zeros((dim[2]*dim[1], 2*dim[2]*dim[1]))
     tic = time.clock()
-    pm.phase_mag_real(res, projection, 'slab', b_0, jacobi=jacobi)
+    phase_map = PhaseMap(res, pm.phase_mag_real(res, projection, 'slab', b_0, jacobi=jacobi))
     toc = time.clock()
+    phase_map.display()
     np.savetxt('../output/jacobi.npy', jacobi)
     print 'Time for Real Space Approach with Jacobi-Matrix (Slab): ' + str(toc - tic)
     
@@ -51,7 +51,7 @@ def phase_from_mag():
     
 if __name__ == "__main__":
     try:
-        jacobi = phase_from_mag()
+        jacobi = get_jacobi()
     except:
         type, value, tb = sys.exc_info()
         traceback.print_exc()

scripts/invert2.py

+# -*- coding: utf-8 -*-
+"""Create random magnetic distributions."""
+
+import random as rnd
+import pdb, traceback, sys
+import numpy as np
+from numpy import pi
+import pylab
+from copy import deepcopy
+
+import pyramid.magcreator as mc
+from pyramid.magdata import MagData
+import pyramid.phasemapper as pm
+import pyramid.projector as pj
+import pyramid.holoimage as hi
+from pyramid.phasemap import PhaseMap
+from scipy.optimize import leastsq
+
+
+def create_random_distribution():
+    '''Calculate, display and save a random magnetic distribution to file.
+    Arguments:
+        None
+    Returns:
+        None
+
+    '''
+    # Input parameters:
+    count = 200
+    dim = (1, 32, 32)
+    b_0 = 1    # in T
+    res = 10 # in nm
+    rnd.seed(12)
+    # Create lists for magnetic objects:
+    mag_shape_list = np.zeros((count,) + dim)
+    beta_list      = np.zeros(count)
+    magnitude_list = np.zeros(count)
+    for i in range(count):
+        pixel = (rnd.randrange(dim[0]), rnd.randrange(dim[1]), rnd.randrange(dim[2]))
+        mag_shape_list[i,...] = mc.Shapes.pixel(dim, pixel)
+        beta_list[i] = 2*pi*rnd.random()
+        magnitude_list[i] = rnd.random()
+    # Create magnetic distribution
+    magnitude = mc.create_mag_dist_comb(mag_shape_list, beta_list, magnitude_list)
+    mag_data = MagData(res, magnitude)
+#    mag_data.quiver_plot()
+#    mag_data.save_to_llg('output/mag_dist_random_pixel.txt')
+
+
+    magx, magy, magz = mag_data.magnitude
+    maskx, masky, maskz = magx != 0, magy != 0, magz !=0
+    x_t = np.concatenate([magx[maskx], magy[masky], magz[maskz]])
+    print "x_t", x_t
+
+    def F(x):
+        mag_data_temp = MagData(deepcopy(mag_data.res), deepcopy(mag_data.magnitude))
+        magx, magy, magz = mag_data_temp.magnitude
+        maskx, masky, maskz = magx != 0, magy != 0, magz !=0
+#        print maskx.sum() + masky.sum() + maskz.sum()
+        assert len(x) == maskx.sum() + masky.sum() + maskz.sum()
+        magx[maskx] = x[:maskx.sum()]
+        magy[masky] = x[maskx.sum():maskx.sum() + masky.sum()]
+        magz[maskz] = x[maskx.sum() + masky.sum():]
+        projection = pj.simple_axis_projection(mag_data_temp)
+        phase_map_slab = PhaseMap(res, pm.phase_mag_real_ANGLE(res, projection, 'slab', b_0))
+        return phase_map_slab.phase.reshape(-1)
+
+    y_m = F(x_t)
+    print "y_m", y_m
+
+    lam = 1e-6
+
+    def J(x_i):
+#        print "x_i", x_i
+        y_i = F(x_i)
+#        dd1 = np.zeros(mx.shape)
+#        dd2 = np.zeros(mx.shape)
+        term1 = (y_i - y_m)
+        term2 = lam * x_i
+#        dd1[:, :-1, :] += np.diff(mx, axis=1)
+#        dd1[:, -1, :] += np.diff(mx, axis=1)[:, -1, :]
+#        dd1[:, :, :-1] += np.diff(mx, axis=2)
+#        dd1[:, :, -1] += np.diff(mx, axis=2)[:, :, -1]
+
+#        dd2[:, :-1, :] += np.diff(my, axis=1)
+#        dd2[:, -1, :] += np.diff(my, axis=1)[:, -1, :]
+#        dd2[:, :, :-1] += np.diff(my, axis=2)
+#        dd2[:, :, -1] += np.diff(my, axis=2)[:, :, -1]
+
+#        result = np.concatenate([term1, np.sqrt(abs(dd1.reshape(-1))), np.sqrt(abs(dd2.reshape(-1)))])
+#        result = np.concatenate([term1, np.sqrt(abs(dd1.reshape(-1)))])
+#        print result
+        return np.concatenate([term1, term2])
+
+    x_f, _ = leastsq(J, np.zeros(x_t.shape))
+    y_f = F(x_f)
+#    print "y_m", y_m
+#    print "y_f", y_f
+#    print "dy", y_f - y_m
+#    print "x_t", x_t
+#    print "x_f", x_f
+#    print "dx", x_f - x_t
+#    pylab.show()
+    projection = pj.simple_axis_projection(mag_data)
+    phase_map  = PhaseMap(res, pm.phase_mag_real(res, projection, 'slab'))
+    hi.display(hi.holo_image(phase_map, 10))
+
+
+if __name__ == "__main__":
+    try:
+        create_random_distribution()
+    except:
+        type, value, tb = sys.exc_info()
+        traceback.print_exc()
+        pdb.post_mortem(tb)