Jacobi-matrix calculation via new Kernel helper class!

phasemapper: added Kernel helper class with methods to calculate the jacobi
matrix or to multiply a vector with the (transposed) jacobi matrix without
explicitly calculating the full matrix
magcreator: added Shape.ellipsoid
numcore: added (experimental) speed-up for jacobi-matrix (not faster)
scripts: modified get_jacobi to test the new Kernel class,
tests: test_compliance now also searches and creates a list of TODOs

pyramid/magcreator.py

@@ -137,6 +137,36 @@ class Shapes:
                             for z in range(dim[0])])
         return mag_shape
 
+    @classmethod
+    def ellipsoid(cls, dim, center, width):
+        '''Create the shape of a sphere.
+
+        Parameters
+        ----------
+        dim : tuple (N=3)
+            The dimensions of the grid `(z, y, x)`.
+        center : tuple (N=3)
+            The center of the sphere in pixel coordinates `(z, y, x)`.
+        radius : float
+            The radius of the sphere in pixel coordinates.
+
+        Returns
+        -------
+        mag_shape : :class:`~numpy.ndarray` (N=3)
+            The magnetic shape as a 3D-array with values between 1 and 0.
+
+        '''
+        assert np.shape(dim) == (3,), 'Parameter dim has to be a a tuple of length 3!'
+        assert np.shape(center) == (3,), 'Parameter center has to be a a tuple of length 3!'
+#        assert radius > 0 and np.shape(radius) == (), 'Radius has to be a positive scalar value!'  # TODO: change
+        mag_shape = np.array([[[np.sqrt((x-center[2])**2/(width[2]/2)**2
+                                      + (y-center[1])**2/(width[1]/2)**2
+                                      + (z-center[0])**2/(width[0]/2)**2) <= 1
+                            for x in range(dim[2])]
+                            for y in range(dim[1])]
+                            for z in range(dim[0])])
+        return mag_shape
+
     @classmethod
     def filament(cls, dim, pos, axis='y'):
         '''Create the shape of a filament.

pyramid/numcore/phase_mag_real.pyx

@@ -58,3 +58,50 @@ def phase_mag_real_core(
                 for q in range(v_dim):
                     for p in range(u_dim):
                         phase[q, p] -= v_m * v_phi[q_c + q, p_c + p]
+
+
+@cython.boundscheck(False)
+@cython.wraparound(False)
+def get_jacobi_core(
+    unsigned int v_dim, unsigned int u_dim,
+    double[:, :] v_phi, double[:, :] u_phi,
+    double[:, :] jacobi):
+    '''DOCSTRING!'''
+    # TODO: Docstring!!!
+    cdef unsigned int i, j, p, q, p_min, p_max, q_min, q_max, u_column, v_column, row
+    for j in range(v_dim):
+        for i in range(u_dim):
+            q_min = (v_dim-1) - j
+            q_max = (2*v_dim-1) - j
+            p_min = (u_dim-1) - i
+            p_max = (2*u_dim-1) - i
+            v_column = i + u_dim*j + u_dim*v_dim
+            u_column = i + u_dim*j
+            for q in range(q_min, q_max):
+                for p in range(p_min, p_max):
+                    q_c = q - (v_dim-1-j)  # start at zero for jacobi column
+                    p_c = p - (u_dim-1-i)  # start at zero for jacobi column
+                    row = p_c + u_dim*q_c
+                    # u-component:
+                    jacobi[row, u_column] =  u_phi[q, p]
+                    # v-component (note the minus!):
+                    jacobi[row, v_column] = -v_phi[q, p]
+
+
+@cython.boundscheck(False)
+@cython.wraparound(False)
+def get_jacobi_core(
+    unsigned int v_dim, unsigned int u_dim,
+    np.ndarray v_phi, np.ndarray u_phi,
+    np.ndarray jacobi):
+    '''DOCSTRING!'''
+    # TODO: Docstring!!!
+    cdef unsigned int i, j, p, q, p_min, p_max, q_min, q_max, u_column, v_column, row
+    for j in range(v_dim):
+        for i in range(u_dim):
+            # v_dim*u_dim columns for the u-component:
+            jacobi[:, i+u_dim*j] = \
+                np.reshape(u_phi[v_dim-1-j:(2*v_dim-1)-j, u_dim-1-i:(2*u_dim-1)-i], -1)
+            # v_dim*u_dim columns for the v-component (note the minus!):
+            jacobi[:, u_dim*v_dim+i+u_dim*j] = \
+                -np.reshape(v_phi[v_dim-1-j:(2*v_dim-1)-j, u_dim-1-i:(2*u_dim-1)-i], -1)

pyramid/phasemapper.py

@@ -26,55 +26,99 @@ Q_E = 1.602E-19    # electron charge in C
 C = 2.998E8        # speed of light in m/s
 
 
-def phase_mag_fourier(res, projection, padding=0, b_0=1):
-    '''Calculate the magnetic phase from magnetization data (Fourier space approach).
-
-    Parameters
-    ----------
-    res : float
-        The resolution of the grid (grid spacing) in nm.
-    projection : tuple (N=3) of :class:`~numpy.ndarray` (N=2)
-        The in-plane projection of the magnetization as a tuple, storing the `u`- and `v`-component
-        of the magnetization and the thickness projection for the resulting 2D-grid.
-    padding : int, optional
-        Factor for the zero padding. The default is 0 (no padding). For a factor of n the number
-        of pixels is increase by ``(1+n)**2``. Padded zeros are cropped at the end.
-    b_0 : float, optional
-        The magnetic induction corresponding to a magnetization `M`\ :sub:`0` in T.
-        The default is 1.
-
-    Returns
-    -------
-    phase : :class:`~numpy.ndarray` (N=2)
-        The phase as a 2-dimensional array.
+# TODO: Kernel class? Creator, transform to fourier, get_jacobi?
 
-    '''
-    v_dim, u_dim = np.shape(projection[0])
-    v_mag, u_mag = projection[:-1]
-    # Create zero padded matrices:
-    u_pad = u_dim/2 * padding
-    v_pad = v_dim/2 * padding
-    u_mag_big = np.zeros(((1 + padding) * v_dim, (1 + padding) * u_dim))
-    v_mag_big = np.zeros(((1 + padding) * v_dim, (1 + padding) * u_dim))
-    u_mag_big[v_pad:v_pad+v_dim, u_pad:u_pad+u_dim] = u_mag
-    v_mag_big[v_pad:v_pad+v_dim, u_pad:u_pad+u_dim] = v_mag
-    # Fourier transform of the two components:
-    u_mag_fft = np.fft.fftshift(np.fft.rfft2(u_mag_big), axes=0)
-    v_mag_fft = np.fft.fftshift(np.fft.rfft2(v_mag_big), axes=0)
-    # Calculate the Fourier transform of the phase:
-    nyq = 0.5 / res  # nyquist frequency
-    f_u = np.linspace(0, nyq, u_mag_fft.shape[1])
-    f_v = np.linspace(-nyq, nyq, u_mag_fft.shape[0], endpoint=False)
-    f_uu, f_vv = np.meshgrid(f_u, f_v)
-    coeff = (1j*b_0) / (2*PHI_0)
-    phase_fft = coeff * res * (u_mag_fft*f_vv - v_mag_fft*f_uu) / (f_uu**2 + f_vv**2 + 1e-30) #* (8*(res/2)**3)*np.sinc(res/2*f_uu)*np.sinc(res/2*f_vv) * np.exp(2*pi*1j*())
-    # Transform to real space and revert padding:
-    phase_big = np.fft.irfft2(np.fft.ifftshift(phase_fft, axes=0))
-    phase = phase_big[v_pad:v_pad+v_dim, u_pad:u_pad+u_dim]
-    return phase
-
-
-def phase_mag_real(res, projection, method='discs', b_0=1, jacobi=None):
+class Kernel:
+    # TODO: Docstrings!!!
+    
+    def __init__(self, dim, res, geometry='disc', b_0=1):
+        # TODO: Docstring!!!
+        def get_elementary_phase(geometry, n, m, res):
+            if geometry == '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
+                return 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 geometry == 'disc':
+                in_or_out = np.logical_not(np.logical_and(n == 0, m == 0))
+                return m / (n**2 + m**2 + 1E-30) * in_or_out
+        self.dim = dim  # !!! size of the FOV, kernel is bigger!
+        self.res = res
+        self.geometry = geometry
+        self.b_0 = b_0
+        coeff = -res**2 / (2*PHI_0)
+        v_dim, u_dim = dim
+        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)
+        self.u = coeff * get_elementary_phase(geometry, uu, vv, res)
+        self.v = coeff * get_elementary_phase(geometry, vv, uu, res)
+        size = 3*np.array(dim) - 1  # dim + (2*dim - 1) magnetisation + kernel
+        fsize = 2 ** np.ceil(np.log2(size)).astype(int)  # next multiple of 2
+        self.u_fft = np.fft.rfftn(self.u, fsize)
+        self.v_fft = np.fft.rfftn(self.v, fsize)
+
+    def get_jacobi(self):
+        # TODO: Docstring!!!
+        '''CAUTIOUS! Just use for small dimensions or computer will explode violently!'''
+        v_dim, u_dim = self.dim
+        jacobi = np.zeros((v_dim*u_dim, 2*v_dim*u_dim))  
+    #    nc.get_jacobi_core(dim[0], dim[1], v_phi, u_phi, jacobi)
+    #    return jacobi
+        for j in range(v_dim):
+            for i in range(u_dim):
+                u_column = i + u_dim*j
+                v_column = i + u_dim*j + u_dim*v_dim
+                u_min = (u_dim-1) - i
+                u_max = (2*u_dim-1) - i
+                v_min = (v_dim-1) - j
+                v_max = (2*v_dim-1) - j
+                # u_dim*v_dim columns for the u-component:
+                jacobi[:, u_column] = self.u[v_min:v_max, u_min:u_max].reshape(-1)
+                # u_dim*v_dim columns for the v-component (note the minus!):
+                jacobi[:, v_column] = -self.v[v_min:v_max, u_min:u_max].reshape(-1)
+        return jacobi
+
+    def multiply_jacobi(self, vector):
+        # TODO: Docstring!!!
+        # vector: v_dim*u_dim elements for u_mag and v_dim*u_dim elements for v_mag
+        v_dim, u_dim = self.dim
+        size = v_dim * u_dim
+        assert len(vector) == 2*size, 'vector size not compatible!'
+        result = np.zeros(size)
+        for s in range(size):  # column-wise (two columns at a time, u- and v-component)
+            i = s % u_dim
+            j = int(s/u_dim)
+            u_min = (u_dim-1) - i
+            u_max = (2*u_dim-1) - i
+            v_min = (v_dim-1) - j
+            v_max = (2*v_dim-1) - j
+            result += vector[s]*self.u[v_min:v_max, u_min:u_max].reshape(-1)  # u
+            result += vector[s+size]*-self.v[v_min:v_max, u_min:u_max].reshape(-1)  # v        
+        return result
+
+    def multiply_jacobi_T(self, vector):
+        # TODO: Docstring!!!
+        # vector: v_dim*u_dim elements for u_mag and v_dim*u_dim elements for v_mag
+        v_dim, u_dim = self.dim
+        size = v_dim * u_dim
+        assert len(vector) == size, 'vector size not compatible!'
+        result = np.zeros(2*size)
+        for s in range(size):  # row-wise (two rows at a time, u- and v-component)
+            i = s % u_dim
+            j = int(s/u_dim)
+            u_min = (u_dim-1) - i
+            u_max = (2*u_dim-1) - i
+            v_min = (v_dim-1) - j
+            v_max = (2*v_dim-1) - j
+            result[s] = np.sum(vector*self.u[v_min:v_max, u_min:u_max].reshape(-1))  # u
+            result[s+size] = np.sum(vector*-self.v[v_min:v_max, u_min:u_max].reshape(-1))  # v        
+        return result
+
+
+def phase_mag_real(res, projection, geometry='disc', b_0=1, jacobi=None):
     '''Calculate the magnetic phase from magnetization data (real space approach).
 
     Parameters
@@ -84,7 +128,7 @@ def phase_mag_real(res, projection, method='discs', b_0=1, jacobi=None):
     projection : tuple (N=3) of :class:`~numpy.ndarray` (N=2)
         The in-plane projection of the magnetization as a tuple, storing the `u`- and `v`-component
         of the magnetization and the thickness projection for the resulting 2D-grid.
-    method : {'disc', 'slab'}, optional
+    geometry : {'disc', 'slab'}, optional
         Specifies the elemental geometry to use for the pixel field.
         The default is 'disc', because of the smaller computational overhead.
     b_0 : float, optional
@@ -104,9 +148,10 @@ def phase_mag_real(res, projection, method='discs', b_0=1, jacobi=None):
     dim = np.shape(projection[0])
     v_mag, u_mag = projection[:-1]
 
-    # Get lookup-tables for the phase of one pixel:
-    u_phi = get_kernel(method, 'u', dim, res, b_0)
-    v_phi = get_kernel(method, 'v', dim, res, b_0)
+    # Create kernel (lookup-tables for the phase of one pixel):
+    kernel = Kernel(dim, res, geometry, b_0)
+    u_phi = kernel.u
+    v_phi = kernel.v
 
     # Calculation of the phase:
     phase = np.zeros(dim)
@@ -211,41 +256,99 @@ def phase_mag_real_fast(res, projection, kernels_fourier, b_0=1):
     u_phase = np.fft.irfftn(u_mag_f * u_kern_f, fsize)[fslice].copy()
     v_phase = np.fft.irfftn(v_mag_f * v_kern_f, fsize)[fslice].copy()
     return u_phase - v_phase
-    
-    
-def get_kernel_fourier(method, orientation, dim, res, b_0=1):
 
-    kernel = get_kernel(method, orientation, dim, res, b_0)
- 
+
+def phase_mag_fourier(res, projection, padding=0, b_0=1):
+    '''Calculate the magnetic phase from magnetization data (Fourier space approach).
+
+    Parameters
+    ----------
+    res : float
+        The resolution of the grid (grid spacing) in nm.
+    projection : tuple (N=3) of :class:`~numpy.ndarray` (N=2)
+        The in-plane projection of the magnetization as a tuple, storing the `u`- and `v`-component
+        of the magnetization and the thickness projection for the resulting 2D-grid.
+    padding : int, optional
+        Factor for the zero padding. The default is 0 (no padding). For a factor of n the number
+        of pixels is increase by ``(1+n)**2``. Padded zeros are cropped at the end.
+    b_0 : float, optional
+        The magnetic induction corresponding to a magnetization `M`\ :sub:`0` in T.
+        The default is 1.
+
+    Returns
+    -------
+    phase : :class:`~numpy.ndarray` (N=2)
+        The phase as a 2-dimensional array.
+
+    '''
+    v_dim, u_dim = np.shape(projection[0])
+    v_mag, u_mag = projection[:-1]
+    # Create zero padded matrices:
+    u_pad = u_dim/2 * padding
+    v_pad = v_dim/2 * padding
+    u_mag_big = np.zeros(((1 + padding) * v_dim, (1 + padding) * u_dim))
+    v_mag_big = np.zeros(((1 + padding) * v_dim, (1 + padding) * u_dim))
+    u_mag_big[v_pad:v_pad+v_dim, u_pad:u_pad+u_dim] = u_mag
+    v_mag_big[v_pad:v_pad+v_dim, u_pad:u_pad+u_dim] = v_mag
+    # Fourier transform of the two components:
+    u_mag_fft = np.fft.fftshift(np.fft.rfft2(u_mag_big), axes=0)
+    v_mag_fft = np.fft.fftshift(np.fft.rfft2(v_mag_big), axes=0)
+    # Calculate the Fourier transform of the phase:
+    nyq = 0.5 / res  # nyquist frequency
+    f_u = np.linspace(0, nyq, u_mag_fft.shape[1])
+    f_v = np.linspace(-nyq, nyq, u_mag_fft.shape[0], endpoint=False)
+    f_uu, f_vv = np.meshgrid(f_u, f_v)
+    coeff = (1j*b_0) / (2*PHI_0)
+    phase_fft = coeff * res * (u_mag_fft*f_vv - v_mag_fft*f_uu) / (f_uu**2 + f_vv**2 + 1e-30) #* (8*(res/2)**3)*np.sinc(res/2*f_uu)*np.sinc(res/2*f_vv) * np.exp(2*pi*1j*())
+    # Transform to real space and revert padding:
+    phase_big = np.fft.irfft2(np.fft.ifftshift(phase_fft, axes=0))
+    phase = phase_big[v_pad:v_pad+v_dim, u_pad:u_pad+u_dim]
+    return phase
+
+
+def phase_mag(res, projection, kernel, b_0=1):
+    '''Calculate the magnetic phase from magnetization data (real space approach).
+
+    Parameters
+    ----------
+    res : float
+        The resolution of the grid (grid spacing) in nm.
+    projection : tuple (N=3) of :class:`~numpy.ndarray` (N=2)
+        The in-plane projection of the magnetization as a tuple, storing the `u`- and `v`-component
+        of the magnetization and the thickness projection for the resulting 2D-grid.
+    method : {'disc', 'slab'}, optional
+        Specifies the elemental geometry to use for the pixel field.
+        The default is 'disc', because of the smaller computational overhead.
+    b_0 : float, optional
+        The magnetic induction corresponding to a magnetization `M`\ :sub:`0` in T.
+        The default is 1.
+    jacobi : :class:`~numpy.ndarray` (N=2), optional
+        Specifies the matrix in which to save the jacobi matrix. The jacobi matrix will not be
+        calculated, if no matrix is specified (default), resulting in a faster computation.
+
+    Returns
+    -------
+    phase : :class:`~numpy.ndarray` (N=2)
+        The phase as a 2-dimensional array.
+
+    '''  # TODO: Docstring!!!
+    # Process input parameters:
+    v_mag, u_mag = projection[:-1]
+    dim = np.shape(u_mag)
+
     size = 3*np.array(dim) - 1  # dim + (2*dim - 1) magnetisation + kernel
     fsize = 2 ** np.ceil(np.log2(size)).astype(int)
+    fslice = tuple([slice(0, int(sz)) for sz in size])
 
-    return np.fft.rfftn(kernel, fsize)
+    u_mag_f = np.fft.rfftn(u_mag, fsize)
+    v_mag_f = np.fft.rfftn(v_mag, fsize)
 
+    v_kern_f, u_kern_f = kernels_fourier
 
-def get_kernel(method, orientation, dim, res, b_0=1):
-    
-    def get_elementary_phase(method, n, m, res):
-        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
-            return 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':
-            in_or_out = np.logical_not(np.logical_and(n == 0, m == 0))
-            return m / (n**2 + m**2 + 1E-30) * in_or_out
-    
-    coeff = -b_0 * res**2 / (2*PHI_0)
-    v_dim, u_dim = dim
-    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)
-    if orientation == 'u':
-        return coeff * get_elementary_phase(method, uu, vv, res)
-    elif orientation == 'v':
-        return coeff * get_elementary_phase(method, vv, uu, res)
+    fslice = (slice(dim[0]-1, 2*dim[0]-1), slice(dim[1]-1, 2*dim[1]-1))
+    u_phase = np.fft.irfftn(u_mag_f * u_kern_f, fsize)[fslice].copy()
+    v_phase = np.fft.irfftn(v_mag_f * v_kern_f, fsize)[fslice].copy()
+    return u_phase - v_phase
 
 
 def phase_elec(res, projection, v_0=1, v_acc=30000):

scripts/compare methods/compare_pixel_fields.py

@@ -33,7 +33,7 @@ def compare_pixel_fields():
     # Input parameters:
     res = 1.0  # in nm
     phi = 0  # in rad
-    dim = (1, 32, 32)
+    dim = (1, 64, 64)
     pixel = (0,  int(dim[1]/2),  int(dim[2]/2))
     limit = 0.35
 
@@ -57,37 +57,37 @@ def compare_pixel_fields():
     projection = pj.simple_axis_projection(mag_data)
     # Kernel of a disc in real space:
     phase_map_disc = PhaseMap(res, pm.phase_mag_real(res, projection, 'disc'), 'mrad')
-    phase_map_disc.display('Phase of one Pixel (Disc)', limit=limit)
+#    phase_map_disc.display('Phase of one Pixel (Disc)', limit=limit)
     # Kernel of a slab in real space:
     phase_map_slab = PhaseMap(res, pm.phase_mag_real(res, projection, 'slab'), 'mrad')
-    phase_map_slab.display('Phase of one Pixel (Slab)', limit=limit)
+#    phase_map_slab.display('Phase of one Pixel (Slab)', limit=limit)
     # Kernel of the Fourier method:
     phase_map_fft = PhaseMap(res, pm.phase_mag_fourier(res, projection, padding=0), 'mrad')
     phase_map_fft.display('Phase of one Pixel (Fourier)', limit=limit)
     
     # Kernel of the Fourier method, calculated directly:
     phase_map_fft_kernel = PhaseMap(res, get_fourier_kernel(), 'mrad')
-    phase_map_fft_kernel.display('Phase of one Pixel (Fourier Kernel)', limit=limit)
+#    phase_map_fft_kernel.display('Phase of one Pixel (Fourier Kernel)', limit=limit)
     # Kernel differences:
     print 'Fourier Kernel, direct and indirect method are identical:', \
           np.all(phase_map_fft_kernel.phase - phase_map_fft.phase) == 0
     phase_map_diff = PhaseMap(res, phase_map_disc.phase - phase_map_fft.phase, 'mrad')
-    phase_map_diff.display('Phase difference of one Pixel (Disc - Fourier)')
+#    phase_map_diff.display('Phase difference of one Pixel (Disc - Fourier)')
 
     phase_inv_fft = np.abs(np.fft.ifftshift(np.fft.rfft2(phase_map_fft.phase), axes=0))**2
     phase_map_inv_fft = PhaseMap(res, phase_inv_fft)
-    phase_map_inv_fft.display('FT of the Phase of one Pixel (Fourier, Power)')
+#    phase_map_inv_fft.display('FT of the Phase of one Pixel (Fourier, Power)')
 
     phase_inv_disc = np.abs(np.fft.ifftshift(np.fft.rfft2(phase_map_disc.phase), axes=0))**2
     phase_map_inv_disc = PhaseMap(res, phase_inv_disc)
-    phase_map_inv_disc.display('FT of the Phase of one Pixel (Disc, Power)')
+#    phase_map_inv_disc.display('FT of the Phase of one Pixel (Disc, Power)')
     
     import matplotlib.pyplot as plt
     from matplotlib.ticker import IndexLocator 
     
     fig = plt.figure()
     axis = fig.add_subplot(1, 1, 1)
-    x = range(dim[1])
+    x = np.linspace(-dim[1]/res/2, dim[1]/res/2-1, dim[1])
     y_ft = phase_map_fft.phase[:, dim[1]/2]
     y_re = phase_map_disc.phase[:, dim[1]/2]
     axis.axhline(0, color='k')
@@ -97,8 +97,9 @@ def compare_pixel_fields():
     axis.grid()
     axis.legend()
     axis.set_title('Real Space Kernel')
-    axis.set_xlim(0, dim[1]-1)
+    axis.set_xlim(-dim[1]/2, dim[1]/2-1)
     axis.xaxis.set_major_locator(IndexLocator(base=dim[1]/8, offset=0))
+
     
     fig = plt.figure()
     axis = fig.add_subplot(1, 1, 1)

scripts/create distributions/create_core_shell_disc.py

@@ -34,7 +34,7 @@ def create_core_shell_disc():
     filename = directory + '/mag_dist_core_shell_disc.txt'
     res = 1.0  # in nm
     density = 1
-    dim = (64, 64, 64)
+    dim = (32, 32, 32)
     center = (dim[0]/2-0.5, int(dim[1]/2)-0.5, int(dim[2]/2)-0.5)
     radius_core = dim[1]/8
     radius_shell = dim[1]/4

scripts/reconstruction/get_jacobi.py

@@ -44,18 +44,56 @@ def get_jacobi():
 
     mag_data = MagData(res, mc.create_mag_dist_homog(mc.Shapes.slab(dim, center, width), phi))
     projection = pj.simple_axis_projection(mag_data)
+    print 'Projection calculated!'
 
     '''NUMERICAL SOLUTION'''
-    # numerical solution Real Space (Slab):
-    jacobi = np.zeros((dim[2]*dim[1], 2*dim[2]*dim[1]))
+    # numerical solution Real Space:
+    dim_proj = np.shape(projection[0])
+    size = np.prod(dim_proj)
+    kernel = pm.Kernel(dim_proj, res, 'disc')
+
+    tic = time.clock()
+    kernel.multiply_jacobi(np.ones(2*size))  # column-wise
+    toc = time.clock()
+    print 'Time for one multiplication with the Jacobi-Matrix:      ', toc - tic
+
+    jacobi = np.zeros((size, 2*size))
     tic = time.clock()
-    phase_map = PhaseMap(res, pm.phase_mag_real(res, projection, 'slab', b_0, jacobi=jacobi))
+    phase_map = PhaseMap(res, pm.phase_mag_real(res, projection, 'disc', b_0, jacobi=jacobi))
     toc = time.clock()
     phase_map.display()
     np.savetxt(filename, jacobi)
-    print 'Time for Real Space Approach with Jacobi-Matrix (Slab): ' + str(toc - tic)
+    print 'Time for Jacobi-Matrix during phase calculation:         ', toc - tic
+    
+    tic = time.clock()
+    jacobi_test = kernel.get_jacobi()
+    toc = time.clock()
+    print 'Time for Jacobi-Matrix from the Kernel:                  ', toc - tic
+    
+    unity = np.eye(2*size)
+    jacobi_test2 = np.zeros((size, 2*size))
+    tic = time.clock()
+    for i in range(unity.shape[1]):
+        jacobi_test2[:, i] = kernel.multiply_jacobi(unity[:, i])  # column-wise
+    toc = time.clock()
+    print 'Time for getting the Jacobi-Matrix (vector-wise):        ', toc - tic
+
+    unity_transp = np.eye(size)
+    jacobi_transp = np.zeros((2*size, size))
+    tic = time.clock()
+    for i in range(unity_transp.shape[1]):
+        jacobi_transp[:, i] = kernel.multiply_jacobi_T(unity_transp[:, i])  # column-wise
+    toc = time.clock()
+    print 'Time for getting the transp. Jacobi-Matrix (vector-wise):', toc - tic
 
-    return jacobi
+    print 'Methods (during vs kernel) in accordance?                ', \
+        np.logical_not(np.all(jacobi-jacobi_test))
+    print 'Methods (during vs vector-wise) in accordance?           ', \
+        np.logical_not(np.all(jacobi-jacobi_test2))
+    print 'Methods (transponed Jacobi) in accordance?               ', \
+        np.logical_not(np.all(jacobi.T-jacobi_transp))
+    
+    return (jacobi, jacobi_test, jacobi_test2, jacobi_transp)
 
 
 if __name__ == "__main__":

scripts/reconstruction/reconstruct_core_shell.py

@@ -89,6 +89,5 @@ mask = mag_data.get_mask()
 dim = mag_data.dim
 res = mag_data.res
 mag_data_reconstruct = MagData(res, (np.zeros(dim),)*3)
-print np.shape(mag_data_reconstruct.magnitude)
 
 mag_vector = mag_data_reconstruct.get_vector(mask)

scripts/test_h5py.py

+# -*- coding: utf-8 -*-
+"""
+Created on Thu Nov 21 17:35:37 2013
+
+@author: Jan
+"""
+
+import h5py
+import numpy as np
+
+with h5py.File('testfile.hdf5', 'w') as f:
+
+    dset = f.create_dataset('testdataset', (100,), dtype='i')
+    
+    print 'dset.shape:', dset.shape
+    
+    print 'dset.dtype:', dset.dtype
+    
+    dset[...] = np.arange(100)
+    
+    print 'dset[0]:', dset[0]
+
+    print 'dset[10]:', dset[10]
+    
+    print 'dset[0:100:10]:', dset[0:100:10]
+    
+    print 'dset.name:', dset.name
+    
+    print 'f.name:', f.name
+    
+    grp = f.create_group('subgroup')
+    
+    dset_big = grp.create_dataset('another_set', (1000, 1000, 1000), dtype='f')
+    
+    for i in range(dset_big.shape[0]):
+        print 'i:', i
+        dset_big[i, ...] = np.ones(dset_big.shape[1:])

scripts/test_methods.py

@@ -35,7 +35,7 @@ def compare_methods():
     geometry = 'sphere'
     if geometry == 'slab':
         center = (dim[0]/2-0.5, dim[1]/2-0.5, dim[2]/2.-0.5)  # in px (z, y, x) index starts at 0!
-        width = (dim[0]/2, dim[1]/2., dim[2]/2.)  # in px (z, y, x)
+        width = (dim[0]/3., dim[1]/4., dim[2]/2.)  # in px (z, y, x)
         mag_shape = mc.Shapes.slab(dim, center, width)
     elif geometry == 'disc':
         center = (dim[0]/2-0.5, dim[1]/2.-0.5, dim[2]/2.-0.5)  # in px (z, y, x) index starts at 0!
@@ -58,7 +58,6 @@ def compare_methods():
     import time
     start = time.time()
     projection = pj.single_tilt_projection(mag_data, tilt)
-    pj.quiver_plot(projection)
     print 'Total projection time:', time.time() - start
     # Construct phase maps:
     phase_map_mag = PhaseMap(res, pm.phase_mag_fourier(res, projection, padding=1))

scripts/test_vadim.py

@@ -47,7 +47,7 @@ if not os.path.exists(directory):
 res = 1.0  # in nm
 phi = pi/4
 density = 30
-dim = (64, 256, 256)  # in px (z, y, x)
+dim = (256, 256, 64)  # in px (z, y, x)
 
 l = dim[2]/4.
 r = dim[2]/4. + dim[2]/2.
@@ -59,99 +59,118 @@ center = (dim[0]/2.-0.5, dim[1]/2.-0.5, dim[2]/2.-0.5)  # in px (z, y, x) index
 radius = dim[1]/8  # in px
 height = dim[0]/4  # in px
 width = (dim[0]/4, dim[1]/4, dim[2]/4)
-mag_shape_sphere = mc.Shapes.sphere(dim, center, radius)
-mag_shape_disc = mc.Shapes.disc(dim, center, radius, height)
-mag_shape_slab = mc.Shapes.slab(dim, center, (height, radius*2, radius*2))
+#mag_shape_sphere = mc.Shapes.sphere(dim, center, radius)
+#mag_shape_disc = mc.Shapes.disc(dim, center, radius, height)
+#mag_shape_slab = mc.Shapes.slab(dim, center, (height, radius*2, radius*2))
+mag_shape_ellipse = mc.Shapes.ellipsoid(dim, (center[0], 0*center[1], center[2]), (radius*5, radius*10, radius*2))
 # Create magnetization distributions:
-mag_data_sphere = MagData(res, mc.create_mag_dist_homog(mag_shape_sphere, phi))
-mag_data_disc = MagData(res, mc.create_mag_dist_homog(mag_shape_disc, phi))
-mag_data_vortex = MagData(res, mc.create_mag_dist_vortex(mag_shape_disc))
-mag_data_slab = MagData(res, mc.create_mag_dist_homog(mag_shape_slab, phi))
+#mag_data_sphere = MagData(res, mc.create_mag_dist_homog(mag_shape_sphere, phi))
+#mag_data_disc = MagData(res, mc.create_mag_dist_homog(mag_shape_disc, phi))
+#mag_data_vortex = MagData(res, mc.create_mag_dist_vortex(mag_shape_disc))
+#mag_data_slab = MagData(res, mc.create_mag_dist_homog(mag_shape_slab, phi))
+mag_data_ellipse = MagData(res, mc.create_mag_dist_homog(mag_shape_ellipse, phi))
 # Project the magnetization data:
-projection_sphere = pj.simple_axis_projection(mag_data_sphere)
-projection_disc = pj.simple_axis_projection(mag_data_disc)
-projection_vortex = pj.simple_axis_projection(mag_data_vortex)
-projection_slab = pj.simple_axis_projection(mag_data_slab)
-# Magnetic phase map of a homogeneously magnetized disc:
-phase_map_mag_disc = PhaseMap(res, pm.phase_mag_real_conv(res, projection_disc))
-phase_map_mag_disc.save_to_txt(directory+'phase_map_mag_disc.txt')
-axis, _ = hi.display_combined(phase_map_mag_disc, density)
+#projection_sphere = pj.simple_axis_projection(mag_data_sphere)
+#projection_disc = pj.simple_axis_projection(mag_data_disc)
+#projection_vortex = pj.simple_axis_projection(mag_data_vortex)
+#projection_slab = pj.simple_axis_projection(mag_data_slab)
+projection_ellipse = pj.simple_axis_projection(mag_data_ellipse)
+## Magnetic phase map of a homogeneously magnetized disc:
+#phase_map_mag_disc = PhaseMap(res, pm.phase_mag_real_conv(res, projection_disc))
+#phase_map_mag_disc.save_to_txt(directory+'phase_map_mag_disc.txt')
+#axis, _ = hi.display_combined(phase_map_mag_disc, density)
+#axis.axvline(l, b/dim[1], t/dim[1], linewidth=2, color='g')
+#axis.axvline(r, b/dim[1], t/dim[1], linewidth=2, color='g')
+#axis.axhline(b, l/dim[2], r/dim[2], linewidth=2, color='g')
+#axis.axhline(t, l/dim[2], r/dim[2], linewidth=2, color='g')
+#axis.arrow(l+(r-l)/2, b, 0, t-b, length_includes_head=True,
+#              head_width=(r-l)/10, head_length=(r-l)/10, fc='g', ec='g')
+#plt.savefig(directory+'phase_map_mag_disc.png')
+#x, y = calculate_charge_batch(phase_map_mag_disc)
+#np.savetxt(directory+'charge_integral_mag_disc.txt', np.array([x, y]).T)
+#plt.figure()
+#plt.plot(x, y)
+#plt.savefig(directory+'charge_integral_mag_disc.png')
+## Magnetic phase map of a vortex state disc:
+#phase_map_mag_vortex = PhaseMap(res, pm.phase_mag_real_conv(res, projection_vortex))
+#phase_map_mag_vortex.save_to_txt(directory+'phase_map_mag_vortex.txt')
+#axis, _ = hi.display_combined(phase_map_mag_vortex, density)
+#axis.axvline(l, b/dim[1], t/dim[1], linewidth=2, color='g')
+#axis.axvline(r, b/dim[1], t/dim[1], linewidth=2, color='g')
+#axis.axhline(b, l/dim[2], r/dim[2], linewidth=2, color='g')
+#axis.axhline(t, l/dim[2], r/dim[2], linewidth=2, color='g')
+#axis.arrow(l+(r-l)/2, b, 0, t-b, length_includes_head=True,
+#              head_width=(r-l)/10, head_length=(r-l)/10, fc='g', ec='g')
+#plt.savefig(directory+'phase_map_mag_vortex.png')
+#x, y = calculate_charge_batch(phase_map_mag_vortex)
+#np.savetxt(directory+'charge_integral_mag_vortex.txt', np.array([x, y]).T)
+#plt.figure()
+#plt.plot(x, y)
+#plt.savefig(directory+'charge_integral_mag_vortex.png')
+## MIP phase of a slab:
+#phase_map_mip_slab = PhaseMap(res, pm.phase_elec(res, projection_slab, v_0=1, v_acc=300000))
+#phase_map_mip_slab.save_to_txt(directory+'phase_map_mip_slab.txt')
+#axis, _ = hi.display_combined(phase_map_mip_slab, density)
+#axis.axvline(l, b/dim[1], t/dim[1], linewidth=2, color='g')
+#axis.axvline(r, b/dim[1], t/dim[1], linewidth=2, color='g')
+#axis.axhline(b, l/dim[2], r/dim[2], linewidth=2, color='g')
+#axis.axhline(t, l/dim[2], r/dim[2], linewidth=2, color='g')
+#axis.arrow(l+(r-l)/2, b, 0, t-b, length_includes_head=True,
+#              head_width=(r-l)/10, head_length=(r-l)/10, fc='g', ec='g')
+#plt.savefig(directory+'phase_map_mip_slab.png')
+#x, y = calculate_charge_batch(phase_map_mip_slab)
+#np.savetxt(directory+'charge_integral_mip_slab.txt', np.array([x, y]).T)
+#plt.figure()
+#plt.plot(x, y)
+#plt.savefig(directory+'charge_integral_mip_slab.png')
+## MIP phase of a disc:
+#phase_map_mip_disc = PhaseMap(res, pm.phase_elec(res, projection_disc, v_0=1, v_acc=300000))
+#phase_map_mip_disc.save_to_txt(directory+'phase_map_mip_disc.txt')
+#axis, _ = hi.display_combined(phase_map_mip_disc, density)
+#axis.axvline(l, b/dim[1], t/dim[1], linewidth=2, color='g')
+#axis.axvline(r, b/dim[1], t/dim[1], linewidth=2, color='g')
+#axis.axhline(b, l/dim[2], r/dim[2], linewidth=2, color='g')
+#axis.axhline(t, l/dim[2], r/dim[2], linewidth=2, color='g')
+#axis.arrow(l+(r-l)/2, b, 0, t-b, length_includes_head=True,
+#              head_width=(r-l)/10, head_length=(r-l)/10, fc='g', ec='g')
+#plt.savefig(directory+'phase_map_mip_disc.png')
+#x, y = calculate_charge_batch(phase_map_mip_disc)
+#np.savetxt(directory+'charge_integral_mip_disc.txt', np.array([x, y]).T)
+#plt.figure()
+#plt.plot(x, y)
+#plt.savefig(directory+'charge_integral_mip_disc.png')
+## MIP phase of a sphere:
+#phase_map_mip_sphere = PhaseMap(res, pm.phase_elec(res, projection_sphere, v_0=1, v_acc=300000))
+#phase_map_mip_sphere.save_to_txt(directory+'phase_map_mip_sphere.txt')
+#axis, _ = hi.display_combined(phase_map_mip_sphere, density)
+#axis.axvline(l, b/dim[1], t/dim[1], linewidth=2, color='g')
+#axis.axvline(r, b/dim[1], t/dim[1], linewidth=2, color='g')
+#axis.axhline(b, l/dim[2], r/dim[2], linewidth=2, color='g')
+#axis.axhline(t, l/dim[2], r/dim[2], linewidth=2, color='g')
+#axis.arrow(l+(r-l)/2, b, 0, t-b, length_includes_head=True,
+#              head_width=(r-l)/10, head_length=(r-l)/10, fc='g', ec='g')
+#plt.savefig(directory+'phase_map_mip_sphere.png')
+#x, y = calculate_charge_batch(phase_map_mip_sphere)
+#np.savetxt(directory+'charge_integral_mip_sphere.txt', np.array([x, y]).T)
+#plt.figure()
+#plt.plot(x, y)
+#plt.savefig(directory+'charge_integral_mip_sphere.png')
+# MIP phase of an ellipsoid:
+phase_map_mip_ellipsoid = PhaseMap(res, pm.phase_elec(res, projection_ellipse, v_0=1, v_acc=300000))
+phase_map_mip_ellipsoid.save_to_txt(directory+'phase_map_mip_ellipsoid.txt')
+axis, _ = hi.display_combined(phase_map_mip_ellipsoid, density)
 axis.axvline(l, b/dim[1], t/dim[1], linewidth=2, color='g')
 axis.axvline(r, b/dim[1], t/dim[1], linewidth=2, color='g')
 axis.axhline(b, l/dim[2], r/dim[2], linewidth=2, color='g')
 axis.axhline(t, l/dim[2], r/dim[2], linewidth=2, color='g')
 axis.arrow(l+(r-l)/2, b, 0, t-b, length_includes_head=True,
               head_width=(r-l)/10, head_length=(r-l)/10, fc='g', ec='g')
-plt.savefig(directory+'phase_map_mag_disc.png')
-x, y = calculate_charge_batch(phase_map_mag_disc)
-np.savetxt(directory+'charge_integral_mag_disc.txt', np.array([x, y]).T)
+plt.savefig(directory+'phase_map_mip_ellipsoid.png')
+x, y = calculate_charge_batch(phase_map_mip_ellipsoid)
+np.savetxt(directory+'charge_integral_mip_ellipsoid.txt', np.array([x, y]).T)
 plt.figure()
 plt.plot(x, y)
-plt.savefig(directory+'charge_integral_mag_disc.png')
-# Magnetic phase map of a vortex state disc:
-phase_map_mag_vortex = PhaseMap(res, pm.phase_mag_real_conv(res, projection_vortex))
-phase_map_mag_vortex.save_to_txt(directory+'phase_map_mag_vortex.txt')
-axis, _ = hi.display_combined(phase_map_mag_vortex, density)
-axis.axvline(l, b/dim[1], t/dim[1], linewidth=2, color='g')
-axis.axvline(r, b/dim[1], t/dim[1], linewidth=2, color='g')
-axis.axhline(b, l/dim[2], r/dim[2], linewidth=2, color='g')
-axis.axhline(t, l/dim[2], r/dim[2], linewidth=2, color='g')
-axis.arrow(l+(r-l)/2, b, 0, t-b, length_includes_head=True,
-              head_width=(r-l)/10, head_length=(r-l)/10, fc='g', ec='g')
-plt.savefig(directory+'phase_map_mag_vortex.png')
-x, y = calculate_charge_batch(phase_map_mag_vortex)
-np.savetxt(directory+'charge_integral_mag_vortex.txt', np.array([x, y]).T)
-plt.figure()
-plt.plot(x, y)
-plt.savefig(directory+'charge_integral_mag_vortex.png')
-# MIP phase of a slab:
-phase_map_mip_slab = PhaseMap(res, pm.phase_elec(res, projection_slab, v_0=1, v_acc=300000))
-phase_map_mip_slab.save_to_txt(directory+'phase_map_mip_slab.txt')
-axis, _ = hi.display_combined(phase_map_mip_slab, density)
-axis.axvline(l, b/dim[1], t/dim[1], linewidth=2, color='g')
-axis.axvline(r, b/dim[1], t/dim[1], linewidth=2, color='g')
-axis.axhline(b, l/dim[2], r/dim[2], linewidth=2, color='g')
-axis.axhline(t, l/dim[2], r/dim[2], linewidth=2, color='g')
-axis.arrow(l+(r-l)/2, b, 0, t-b, length_includes_head=True,
-              head_width=(r-l)/10, head_length=(r-l)/10, fc='g', ec='g')
-plt.savefig(directory+'phase_map_mip_slab.png')
-x, y = calculate_charge_batch(phase_map_mip_slab)
-np.savetxt(directory+'charge_integral_mip_slab.txt', np.array([x, y]).T)
-plt.figure()
-plt.plot(x, y)
-plt.savefig(directory+'charge_integral_mip_slab.png')
-# MIP phase of a disc:
-phase_map_mip_disc = PhaseMap(res, pm.phase_elec(res, projection_disc, v_0=1, v_acc=300000))
-phase_map_mip_disc.save_to_txt(directory+'phase_map_mip_disc.txt')
-axis, _ = hi.display_combined(phase_map_mip_disc, density)
-axis.axvline(l, b/dim[1], t/dim[1], linewidth=2, color='g')
-axis.axvline(r, b/dim[1], t/dim[1], linewidth=2, color='g')
-axis.axhline(b, l/dim[2], r/dim[2], linewidth=2, color='g')
-axis.axhline(t, l/dim[2], r/dim[2], linewidth=2, color='g')
-axis.arrow(l+(r-l)/2, b, 0, t-b, length_includes_head=True,
-              head_width=(r-l)/10, head_length=(r-l)/10, fc='g', ec='g')
-plt.savefig(directory+'phase_map_mip_disc.png')
-x, y = calculate_charge_batch(phase_map_mip_disc)
-np.savetxt(directory+'charge_integral_mip_disc.txt', np.array([x, y]).T)
-plt.figure()
-plt.plot(x, y)
-plt.savefig(directory+'charge_integral_mip_disc.png')
-# MIP phase of a sphere:
-phase_map_mip_sphere = PhaseMap(res, pm.phase_elec(res, projection_sphere, v_0=1, v_acc=300000))
-phase_map_mip_sphere.save_to_txt(directory+'phase_map_mip_sphere.txt')
-axis, _ = hi.display_combined(phase_map_mip_sphere, density)
-axis.axvline(l, b/dim[1], t/dim[1], linewidth=2, color='g')
-axis.axvline(r, b/dim[1], t/dim[1], linewidth=2, color='g')
-axis.axhline(b, l/dim[2], r/dim[2], linewidth=2, color='g')
-axis.axhline(t, l/dim[2], r/dim[2], linewidth=2, color='g')
-axis.arrow(l+(r-l)/2, b, 0, t-b, length_includes_head=True,
-              head_width=(r-l)/10, head_length=(r-l)/10, fc='g', ec='g')
-plt.savefig(directory+'phase_map_mip_sphere.png')
-x, y = calculate_charge_batch(phase_map_mip_sphere)
-np.savetxt(directory+'charge_integral_mip_sphere.txt', np.array([x, y]).T)
-plt.figure()
-plt.plot(x, y)
-plt.savefig(directory+'charge_integral_mip_sphere.png')
+plt.savefig(directory+'charge_integral_mip_ellipsoid.png')
 
 
 

tests/test_compliance.py

@@ -7,6 +7,8 @@ import sys
 import datetime
 import unittest
 
+import re
+
 import pep8
 
 
@@ -49,7 +51,24 @@ class TestCaseCompliance(unittest.TestCase):
             if result.total_errors == 0:
                 print 'No Warnings or Errors detected!'
             else:
-                print '\n{} Warnings and Errors detected!'.format(result.total_errors)
+                print '---->   {} Warnings and Errors detected!'.format(result.total_errors)
+            print '\nTODOS:'
+            todos_found = False
+            todo_count = 0
+            regex = ur'# TODO: (.*)'
+            for py_file in files:
+                with open (py_file) as f:
+                    for line in f:
+                        todo = re.findall(regex, line)
+                        if todo and not todo[0]=="(.*)'":
+                            todos_found = True
+                            todo_count += 1
+                            print '{}: {}'.format(f.name, todo[0])
+            if todos_found:
+                print '---->   {} TODOs found!'.format(todo_count)
+            else:
+                print 'No TODOS found!'
+
         sys.stdout = stdout_buffer
 
         error_message = 'Found %s Errors and Warnings!' % result.total_errors