Implementation of jacobi multiplication for the projection.

(At the moment a bit slow, optimization with Cython seems necessary).

pyramid/kernel.py

@@ -16,6 +16,8 @@ information.
 
 import numpy as np
 
+import pyramid.numcore as nc
+
 
 PHI_0 = -2067.83    # magnetic flux in T*nm²
 
@@ -166,11 +168,11 @@ class Kernel:
             i = s % u_dim
             j = int(s/u_dim)
             u_min = (u_dim-1) - i
-            u_max = (2*u_dim-1) - i
+            u_max = (2*u_dim-1) - i  # = u_min + u_dim
             v_min = (v_dim-1) - j
-            v_max = (2*v_dim-1) - j
+            v_max = (2*v_dim-1) - j  # = v_min + v_dim
             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        
+            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):
@@ -204,3 +206,19 @@ class Kernel:
             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 multiply_jacobi_core(self, vector):
+        # TODO: Docstring!
+        v_dim, u_dim = self.dim
+        size = v_dim * u_dim
+        result = np.zeros(size)
+        nc.multiply_jacobi_core(v_dim, u_dim, self.v, self.u, vector, result)
+        return result
+
+    def multiply_jacobi_core2(self, vector):
+        # TODO: Docstring!
+        v_dim, u_dim = self.dim
+        size = v_dim * u_dim
+        result = np.zeros(size)
+        nc.multiply_jacobi_core2(v_dim, u_dim, self.v, self.u, vector, result)
+        return result

pyramid/numcore/phase_mag_real.pyx

@@ -5,7 +5,7 @@ Provides a helper function to speed up :func:`~pyramid.phasemapper.phase_mag_rea
 :mod:`~pyramid.phasemapper`, by using C-speed for the for-loops and by omitting boundary and
 wraparound checks.
 
-"""
+""" # TODO: Docstring!
 
 
 import numpy as np
@@ -104,7 +104,7 @@ def get_jacobi_core(
 
 @cython.boundscheck(False)
 @cython.wraparound(False)
-def get_jacobi_core(
+def get_jacobi_core2(
     unsigned int v_dim, unsigned int u_dim,
     double[:, :] v_phi, double[:, :] u_phi,
     double[:, :] jacobi):
@@ -128,22 +128,23 @@ def multiply_jacobi_core(
     double[:, :] v_phi, double[:, :] u_phi,
     double[:] vector,
     double[:] result):
-
-    cdef unsigned int s, i, j, u_min, u_max, v_min, v_max, ri, u, v, siz
+    '''DOCSTRING!'''
+    # TODO: Docstring!!! Iterate over magnetization and then kernel
+    cdef unsigned int s, i, j, u_min, u_max, v_min, v_max, ri, u, v, size
     cdef double v0, v1
-    siz = u_dim * v_dim
-
-    s = 0
+    size = u_dim * v_dim  # Number of pixels
+    s = 0  # Current pixel (numbered consecutively)
+    # Go over all pixels:
     for i in range(u_dim):
         for j in range(v_dim):
-            u_min = (u_dim - 1) - i
-            v_min = (v_dim - 1) - j
-
-            ri = 0
+            u_min = (u_dim - 1) - i  # u_max = u_min + u_dim
+            v_min = (v_dim - 1) - j  # v_max = v_min + v_dim
+            ri = 0  # Current result component (numbered consecutively)
+            # Go over the current kernel cutout [v_min:v_max, u_min:u_max]:
             for v in range(v_min, v_min + v_dim):
                 for u in range(u_min, u_min + u_dim):
                     result[ri] += vector[s] * u_phi[v, u]
-                    result[ri] -= vector[s + siz] * v_phi[v, u]
+                    result[ri] -= vector[s + size] * v_phi[v, u]
                     ri += 1
             s += 1
 
@@ -155,22 +156,21 @@ def multiply_jacobi_core2(
     double[:, :] v_phi, double[:, :] u_phi,
     double[:] vector,
     double[:] result):
-
-    cdef int s1, s2, i, j, u_min, u_max, v_min, v_max, ri, u, v, siz, j_min, j_max, i_min, i_max
-
-    siz = u_dim * v_dim
-
+    '''DOCSTRING!'''
+    # TODO: Docstring!!! Iterate over kernel and then magnetization
+    cdef int s1, s2, i, j, u_min, u_max, v_min, v_max, ri, u, v, size, j_min, j_max, i_min, i_max
+    size = u_dim * v_dim  # Number of pixels
+    # Go over complete kernel:
     for v in range(2 * v_dim - 1):
         for u in range(2 * u_dim - 1):
             i_min = max(0, (u_dim - 1) - u)
             i_max = min(u_dim, ((2 * u_dim) - 1) - u)
+            j_min = max(0, (v_dim - 1) - v)
+            j_max = min(v_dim, ((2 * v_dim) - 1) - v)
+            # Iterate over 
             for i in range(i_min, i_max):
-                s1 = i * u_dim
-                s2 = s1 + siz
-
-                j_min = max(0, (v_dim - 1) - v)
-                j_max = min(v_dim, ((2 * v_dim) - 1) - v)
-
+                s1 = i * u_dim  # u-column
+                s2 = s1 + size  # v-column
                 u_min = (u_dim - 1) - i
                 v_min = (v_dim - 1) - j_min
                 ri = (v - ((v_dim - 1) - j_min)) * u_dim + (u - u_min)

pyramid/phasemapper.py

@@ -25,22 +25,6 @@ C = 2.998E8        # speed of light in m/s
 
 def phase_mag(a, projection, b_0=1, kernel=None):
     '''Calculate the magnetic phase from magnetization data.
-    def multiply_jacobi_core(self, vector):
-        v_dim, u_dim = self.dim
-        size = v_dim * u_dim
-        result = np.zeros(size)
-        nc.multiply_jacobi_core(
-            v_dim, u_dim, self.v, self.u, vector, result)
-        return result
-
-    def multiply_jacobi_core2(self, vector):
-        v_dim, u_dim = self.dim
-        size = v_dim * u_dim
-        result = np.zeros(size)
-        nc.multiply_jacobi_core2(
-            v_dim, u_dim, self.v, self.u, vector, result)
-        return result
-
 
     Parameters
     ----------

pyramid/projector.py

@@ -182,7 +182,7 @@ class Projection:
 
     ''' # TODO: Docstring!
     
-    def __init__(self, u_mag, v_mag, thickness, weights):
+    def __init__(self, dim_proj, dim_rot, dim_perp, v_mag, u_mag, thickness, weights, tilt):
         '''Constructor for a :class:`~.Kernel` object for representing a kernel matrix.
 
         Parameters
@@ -197,12 +197,153 @@ class Projection:
             The elementary geometry of the single magnetized pixel.
 
         ''' # TODO: Docstring!
-        self.dim = np.shape(u_mag)
+        self.dim_proj = dim_proj  # dimension along the projection axis
+        self.dim_rot = dim_rot    # dimension along the rotation axis
+        self.dim_perp = dim_perp  # dimension along the axis perpendicular to proj. and rotation
 #        self.a = a
 #        self.b_0 = b_0
         self.u = u_mag
         self.v = v_mag
+        self.tilt = tilt
         self.weights = weights
+        self.weights_inv = {}
+        # TODO: Not necessary:
+        for key, value_list in weights.iteritems():
+            for new_key, value in value_list:
+                self.weights_inv.setdefault(new_key, []).append((key, value))
+
+    @classmethod
+    def single_tilt_projection(cls, mag_data, tilt=0):
+        '''
+        Project a magnetization distribution which is tilted around the centered y-axis.
+    
+        Parameters
+        ----------
+        mag_data : :class:`~pyramid.magdata.MagData`
+            A :class:`~pyramid.magdata.MagData` object storing the magnetization distribution,
+            which should be projected.
+        tilt : float, optional
+            The counter-clockwise tilt angle around the y-axis. Default is 1.
+    
+        Returns
+        -------
+        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. The latter
+            has to be multiplied with the grid spacing for a value in nm.
+    
+        '''
+        assert isinstance(mag_data, MagData), 'Parameter mag_data has to be a MagData object!'
+    
+        def get_position(p, m, b, size):
+            y, x = np.array(p)[:, 0]+0.5, np.array(p)[:, 1]+0.5
+            return (y-m*x-b)/np.sqrt(m**2+1) + size/2.
+    
+        def get_impact(pos, r, size):
+            return [x for x in np.arange(np.floor(pos-r), np.floor(pos+r)+1, dtype=int)
+                    if 0 <= x < size]
+    
+        def get_weight(delta, rho):  # use circles to represent the voxels
+            lo, up = delta-rho, delta+rho
+            # Upper boundary:
+            if up >= 1:
+                w_up = 0.5
+            else:
+                w_up = (up*np.sqrt(1-up**2) + np.arctan(up/np.sqrt(1-up**2))) / pi
+            # Lower boundary:
+            if lo <= -1:
+                w_lo = -0.5
+            else:
+                w_lo = (lo*np.sqrt(1-lo**2) + np.arctan(lo/np.sqrt(1-lo**2))) / pi
+            return w_up - w_lo
+    
+        # Set starting variables:
+        # length along projection (proj, z), rotation (rot, y) and perpendicular (perp, x) axis:
+        dim_proj, dim_rot, dim_perp = mag_data.dim
+        z_mag, y_mag, x_mag = mag_data.magnitude
+        mask = mag_data.get_mask()
+        projection = (np.zeros((dim_rot, dim_perp)),
+                      np.zeros((dim_rot, dim_perp)),
+                      np.zeros((dim_rot, dim_perp)))
+        # Creating coordinate list of all voxels:
+        voxels = list(itertools.product(range(dim_proj), range(dim_perp)))
+    
+        # Calculate positions along the projected pixel coordinate system:
+        center = (dim_proj/2., dim_perp/2.)
+        m = np.where(tilt<=pi, -1/np.tan(tilt+1E-30), 1/np.tan(tilt+1E-30))
+        b = center[0] - m * center[1]
+        positions = get_position(voxels, m, b, dim_perp)
+    
+        # Calculate weights:
+        r = 1/np.sqrt(np.pi)  # radius of the voxel circle
+        rho = 0.5 / r
+        weights = {}
+        for i, voxel in enumerate(voxels):
+            voxel_weights = []
+            impacts = get_impact(positions[i], r, dim_perp)
+            for impact in impacts:
+                distance = np.abs(impact+0.5 - positions[i])
+                delta = distance / r
+                voxel_weights.append((impact, get_weight(delta, rho)))
+            weights[voxel] = voxel_weights
+        
+        # Calculate projection with the calculated weights for the voxels (rotation around y-axis):
+        for i, voxel in enumerate(weights):
+            voxel_weights = weights[voxel]
+            for pixel, weight in voxel_weights:
+                # Component parallel to tilt axis (':' goes over all slices):
+                projection[0][:, pixel] += weight * y_mag[voxel[0], :, voxel[1]]
+                # Component perpendicular to tilt axis:
+                projection[1][:, pixel] += weight * (x_mag[voxel[0], :, voxel[1]]*np.cos(tilt)
+                                                   + z_mag[voxel[0], :, voxel[1]]*np.sin(tilt))
+                # Thickness profile:
+                projection[2][:, pixel] += weight * mask[voxel[0], :, voxel[1]]
+        
+        return Projection(dim_proj, dim_rot, dim_perp, projection[0], projection[1], projection[2],
+                          weights, tilt)
+
+    def get_weight_matrix(self):
+        dim_proj, dim_rot, dim_perp = self.dim_proj, self.dim_rot, self.dim_perp
+        weights = self.weights
+        size_2d = dim_rot * dim_perp
+        size_3d = dim_rot * dim_perp * dim_proj
+        weight_matrix = np.zeros((size_2d, size_3d))
+        for voxel in weights:
+            voxel_weights = weights[voxel]
+            for pixel, weight in voxel_weights:
+                for i in range(dim_rot):
+                    pixel_pos = i*dim_rot + pixel
+                    voxel_pos = voxel[0]*dim_proj*dim_rot + i*dim_rot + voxel[1]
+                    weight_matrix[pixel_pos, voxel_pos] = weight
+        return weight_matrix
+
+    def multiply_weight_matrix(self, vector):
+        dim_proj, dim_rot, dim_perp = self.dim_proj, self.dim_rot, self.dim_perp
+        weights = self.weights
+        size_2d = dim_rot * dim_perp
+        result = np.zeros(size_2d)
+        for voxel in weights:
+            voxel_weights = weights[voxel]
+            for pixel, weight in voxel_weights:
+                for i in range(dim_rot):
+                    pixel_pos = i*dim_rot + pixel
+                    voxel_pos = voxel[0]*dim_proj*dim_rot + i*dim_rot + voxel[1]
+                    result[pixel_pos] += weight * vector[voxel_pos]
+        return result
+
+    def multiply_weight_matrix_T(self, vector):
+        dim_proj, dim_rot, dim_perp = self.dim_proj, self.dim_rot, self.dim_perp
+        weights = self.weights
+        size_3d = dim_proj * dim_rot * dim_perp
+        result = np.zeros(size_3d)
+        for voxel in weights:
+            voxel_weights = weights[voxel]
+            for pixel, weight in voxel_weights:
+                for i in range(dim_rot):
+                    pixel_pos = i*dim_rot + pixel
+                    voxel_pos = voxel[0]*dim_proj*dim_rot + i*dim_rot + voxel[1]
+                    result[voxel_pos] += weight * vector[pixel_pos]
+        return result
 
     def get_jacobi(self):
         '''Calculate the Jacobi matrix for the phase calculation from a projected magnetization.
@@ -224,57 +365,28 @@ class Projection:
         Just use for small dimensions, Jacobi Matrix scales with order of ``N**4``.
 
         ''' # TODO: Docstring!
-        v_dim, u_dim = self.dim
-        size = v_dim * u_dim
-        jacobi = np.zeros((v_dim*u_dim, 2*v_dim*u_dim))  
-        weight_matrix = np.zeros(dim)  # Check if correct (probably not, x and z instead of x and y)
-        voxels = list(itertools.product(range(dim_proj), range(dim_perp)))
-        for i, voxel in enumerate(self.weights):
-            voxel_weights = self.weights[voxel]
-            for pixel, weight in voxel_weights:
-                # Component parallel to tilt axis (':' goes over all slices):
-                projection[0][:, pixel] += weight * y_mag[voxel[0], :, voxel[1]]
-                # Component perpendicular to tilt axis:
-                projection[1][:, pixel] += weight * (x_mag[voxel[0], :, voxel[1]]*np.cos(tilt)
-                                                   + z_mag[voxel[0], :, voxel[1]]*np.sin(tilt))
-                # Thickness profile:
-                projection[2][:, pixel] += weight * mask[voxel[0], :, voxel[1]]
-        
-        for i, vx in enumerate(self.weights):
-            vx_pos = vx[1]*3 + vx[0]
-            for px, weight in self.weights[voxel]:
-                px_pos = (px[2]*3 + px[1])*3 + px[0]
-                weight_matrix[vx_pos, px_pos] = weight
-            
-        
-        
-        
-        
-        
-        
-        
-        
-        
-        
-        
-        
-        
-        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)
+        dim_proj, dim_rot, dim_perp = self.dim_proj, self.dim_rot, self.dim_perp
+        size_2d = dim_rot * dim_perp
+        size_3d = dim_rot * dim_perp * dim_proj
+        weights = self.weights
+        tilt = self.tilt
+
+        def get_weight_matrix():
+            weight_matrix = np.zeros((size_2d, size_3d))
+            for voxel in weights:
+                voxel_weights = weights[voxel]
+                for pixel, weight in voxel_weights:
+                    for i in range(dim_rot):
+                        pixel_pos = i*dim_rot + pixel
+                        voxel_pos = voxel[0]*dim_proj*dim_rot + i*dim_rot + voxel[1]
+                        weight_matrix[pixel_pos, voxel_pos] = weight
+            return weight_matrix
+
+        weight_matrix = get_weight_matrix()
+        jacobi = np.zeros((2*size_2d, 3*size_3d))
+        jacobi[:size_2d, :size_3d] = np.cos(tilt) * weight_matrix
+        jacobi[:size_2d, 2*size_3d:] = np.sin(tilt) * weight_matrix
+        jacobi[size_2d:, size_3d:2*size_3d] = weight_matrix
         return jacobi
 
     def multiply_jacobi(self, vector):
@@ -293,19 +405,28 @@ class Projection:
             Product of the Jacobi matrix (which is not explicitely calculated) with the vector.
 
         ''' # TODO: Docstring!
-        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        
+        dim_proj, dim_rot, dim_perp = self.dim_proj, self.dim_rot, self.dim_perp
+        size_2d = dim_rot * dim_perp
+        size_3d = dim_rot * dim_perp * dim_proj
+        weights = self.weights
+        tilt = self.tilt
+        assert len(vector) == 3*size_3d
+
+        def multiply_weight_matrix(vector):
+            result = np.zeros(size_2d)
+            for voxel in weights:
+                voxel_weights = weights[voxel]
+                for pixel, weight in voxel_weights:
+                    for i in range(dim_rot):
+                        pixel_pos = i*dim_rot + pixel
+                        voxel_pos = voxel[0]*dim_proj*dim_rot + i*dim_rot + voxel[1]
+                        result[pixel_pos] += weight * vector[voxel_pos]
+            return result
+        
+        result = np.zeros(2*size_2d)
+        result[:size_2d] = (np.cos(tilt) * multiply_weight_matrix(vector[:size_3d])
+                          + np.sin(tilt) * multiply_weight_matrix(vector[2*size_3d:]))
+        result[size_2d:] = multiply_weight_matrix(vector[size_3d:2*size_3d])
         return result
 
     def multiply_jacobi_T(self, vector):
@@ -325,17 +446,26 @@ class Projection:
             the vector.
 
         ''' # TODO: Docstring!
-        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        
+        dim_proj, dim_rot, dim_perp = self.dim_proj, self.dim_rot, self.dim_perp
+        size_2d = dim_rot * dim_perp
+        size_3d = dim_rot * dim_perp * dim_proj
+        weights = self.weights
+        tilt = self.tilt
+        assert len(vector) == 2*size_2d
+
+        def multiply_weight_matrix_T(vector):
+            result = np.zeros(size_3d)
+            for voxel in weights:
+                voxel_weights = weights[voxel]
+                for pixel, weight in voxel_weights:
+                    for i in range(dim_rot):
+                        pixel_pos = i*dim_rot + pixel
+                        voxel_pos = voxel[0]*dim_proj*dim_rot + i*dim_rot + voxel[1]
+                        result[voxel_pos] += weight * vector[pixel_pos]
+            return result
+
+        result = np.zeros(3*size_3d)
+        result[:size_3d] = np.cos(tilt) * multiply_weight_matrix_T(vector[:size_2d])
+        result[size_3d:2*size_3d] = multiply_weight_matrix_T(vector[size_2d:])
+        result[2*size_3d:] = np.sin(tilt) * multiply_weight_matrix_T(vector[:size_2d])
         return result

pyramid/reconstructor.py

@@ -19,6 +19,8 @@ from scipy.optimize import leastsq
 import pyramid.projector as pj
 import pyramid.phasemapper as pm
 from pyramid.magdata import MagData
+from pyramid.projector import Projection
+from pyramid.kernel import Kernel
 
 
 def reconstruct_simple_leastsq(phase_map, mask, b_0=1):
@@ -73,3 +75,7 @@ def reconstruct_simple_leastsq(phase_map, mask, b_0=1):
     x_rec, _ = leastsq(J, np.zeros(3*count))
     mag_data_rec.set_vector(mask, x_rec)
     return mag_data_rec
+
+def reconstruct_test():
+    product = (kernel.multiply_jacobi_T(projection.multiply_jacobi_T(x))
+             * kernel.multiply_jacobi(projection.multiply_jacobi(x)))
\ No newline at end of file

scripts/get_jacobi_projection.py

@@ -7,14 +7,140 @@ Created on Fri Nov 29 16:45:13 2013
 
 import pyramid.magcreator as mc
 from pyramid.magdata import MagData
-import pyramid.projector as pj
 from pyramid.projector import Projection
+from pyramid.kernel import Kernel
+
+import numpy as np
+from numpy import pi
+
+import time
+
 
 
 a = 1.0
-dim = (2, 2, 2)
-px = (0, 0, 0)
-mag_data = MagData(a, mc.create_mag_dist_homog(mc.Shapes.pixel(dim, px), 0))
+dim = (3, 3, 3)
+px = (1, 1, 1)
+tilt = pi/3
+
+mag_data = MagData(a, mc.create_mag_dist_homog(mc.Shapes.pixel(dim, px), phi=pi/4, theta=pi/4))
+
+proj = Projection.single_tilt_projection(mag_data, tilt)
+
+size_2d = dim[1] * dim[2]
+size_3d = dim[0] * dim[1] * dim[2]
+
+ref_u = proj.u
+ref_v = proj.v
+
+z_mag_vec = np.asmatrix(mag_data.magnitude[0].reshape(-1)).T
+y_mag_vec = np.asmatrix(mag_data.magnitude[1].reshape(-1)).T
+x_mag_vec = np.asmatrix(mag_data.magnitude[2].reshape(-1)).T
+
+mag_vec = np.concatenate((x_mag_vec, y_mag_vec, z_mag_vec))
+
+
+
+# Via full multiplicatoin with weight-matrix:
+
+start = time.clock()
+
+weight_matrix = np.asmatrix(proj.get_weight_matrix())
+
+test_u_wf = (np.cos(tilt) * np.asarray(weight_matrix * x_mag_vec).reshape(dim[1], dim[2])
+           + np.sin(tilt) * np.asarray(weight_matrix * z_mag_vec).reshape(dim[1], dim[2]))
+
+test_v_wf = np.asarray(weight_matrix * y_mag_vec).reshape(dim[1], dim[2])
+
+print 'Time for calculation via full weight-matrix multiplication:  ', time.clock() - start
+
+np.testing.assert_almost_equal(test_u_wf, ref_u, err_msg='u-part does not match!')
+np.testing.assert_almost_equal(test_v_wf, ref_v, err_msg='v-part does not match!')
+
+
+
+# Via direct multiplication with weight_matrix:
+
+start = time.clock()
+
+test_u_wd = (np.cos(tilt) * proj.multiply_weight_matrix(x_mag_vec).reshape(dim[1], dim[2])
+           + np.sin(tilt) * proj.multiply_weight_matrix(z_mag_vec).reshape(dim[1], dim[2]))
+
+test_v_wd = proj.multiply_weight_matrix(y_mag_vec).reshape(dim[1], dim[2])
+
+print 'Time for calculation via direct weight-matrix multiplication:', time.clock() - start
+
+np.testing.assert_almost_equal(test_u_wd, ref_u, err_msg='u-part does not match!')
+np.testing.assert_almost_equal(test_v_wd, ref_v, err_msg='v-part does not match!')
+
+
+
+# Via full multiplication with jacobi-matrix:
+
+start = time.clock()
+
+jacobi = np.asmatrix(proj.get_jacobi())
+projected_mag = np.asarray(jacobi * mag_vec).reshape(2, dim[1], dim[2])
+
+test_u_jf = projected_mag[0, ...]
+test_v_jf = projected_mag[1, ...]
+
+print 'Time for calculation via full jacobi-matrix multiplication:  ', time.clock() - start
+
+np.testing.assert_almost_equal(test_u_jf, ref_u, err_msg='u-part does not match!')
+np.testing.assert_almost_equal(test_v_jf, ref_v, err_msg='v-part does not match!')
+
+
+
+# Via direct multiplication with jacobi-matrix:
+
+start = time.clock()
+
+projected_mag = proj.multiply_jacobi(mag_vec).reshape(2, dim[1], dim[2])
+
+test_u_jd = projected_mag[0, ...]
+test_v_jd = projected_mag[1, ...]
+
+print 'Time for calculation via direct jacobi-matrix multiplication:', time.clock() - start
+
+np.testing.assert_almost_equal(test_u_jd, ref_u, err_msg='u-part does not match!')
+np.testing.assert_almost_equal(test_v_jd, ref_v, err_msg='v-part does not match!')
+
+
+
+# Via full multiplication with transposed jacobi-matrix:
+
+start = time.clock()
+
+jacobi_T = np.asmatrix(proj.get_jacobi()).T
+
+test_T_jf = np.asarray(jacobi_T * np.asmatrix(np.ones(2*size_2d)).T).reshape(-1)
+
+print 'Time for calculation via full transposed multiplication:     ', time.clock() - start
+
+
+
+
+# Via direct multiplication with transposed jacobi-matrix:
+
+start = time.clock()
+
+jacobi_T = np.asmatrix(proj.get_jacobi()).T
+
+test_complete_T = np.asarray(jacobi_T * np.asmatrix(np.ones(2*size_2d)).T).reshape(-1)
+
+test_T_jd = proj.multiply_jacobi_T(np.ones(2*size_2d))
+
+print 'Time for calculation via direct transposed multiplication:   ', time.clock() - start
+
+np.testing.assert_almost_equal(test_T_jd, test_T_jf, err_msg='Transposed vector does not match!')
+
 
-proj = pj.simple_axis_projection(mag_data)
 
+## Cost function testing:
+#
+#kern = Kernel((dim[1], dim[2]), a)
+#
+#identity = np.eye(5)
+#
+#right = kern.multiply_jacobi(proj.multiply_jacobi(mag_vec))
+#cost = kern.multiply_jacobi_T(proj.multiply_jacobi_T(right))
\ No newline at end of file

scripts/reconstruction/get_jacobi.py

@@ -23,7 +23,7 @@ if not os.path.exists(directory):
     os.makedirs(directory)
 filename = directory + '/jacobi.npy'
 b_0 = 1.0  # in T
-dim = (1, 8, 8)  # in px (y,x)
+dim = (1, 64, 64)  # in px (y,x)
 res = 10.0  # in nm
 phi = pi/4
 
@@ -32,51 +32,42 @@ width = (0, 1, 1)  # in px (y,x)
 
 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!'
+print 'Compare times for multiplication of a vector (ones) with the jacobi matrix!'
 
-'''NUMERICAL SOLUTION'''
-# numerical solution Real Space:
+# Prepare kernel and vector:
 dim_proj = np.shape(projection[0])
 size = np.prod(dim_proj)
 kernel = pm.Kernel(dim_proj, res, 'disc')
+vector = np.ones(2*size)
 
+# Calculation with kernel (core):
 tic = time.clock()
-kernel.multiply_jacobi(np.ones(2*size))  # column-wise
+result_core = kernel.multiply_jacobi_core(vector)
 toc = time.clock()
-print 'Time for one multiplication with the Jacobi-Matrix:      ', toc - tic
+print 'Calculation with kernel (core):   ', toc - tic
 
-jacobi = np.zeros((size, 2*size))
+# Calculation with kernel (core2):
 tic = time.clock()
-phase_map = PhaseMap(res, pm.phase_mag_real(res, projection, b_0, 'disc', jacobi=jacobi))
+result_core2 = kernel.multiply_jacobi_core2(vector)
 toc = time.clock()
-phase_map.display()
-np.savetxt(filename, jacobi)
-print 'Time for Jacobi-Matrix during phase calculation:         ', toc - tic
+print 'Calculation with kernel (core2):  ', toc - tic
 
+# Calculation with kernel:
 tic = time.clock()
-jacobi_test = kernel.get_jacobi()
+result_kernel = kernel.multiply_jacobi(vector)
 toc = time.clock()
-print 'Time for Jacobi-Matrix from the Kernel:                  ', toc - tic
+print 'Calculation with kernel:          ', toc - tic
 
-unity = np.eye(2*size)
-jacobi_test2 = np.zeros((size, 2*size))
+# Calculation during phasemapping:
 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
+jacobi = np.zeros((size, 2*size))
+PhaseMap(res, pm.phase_mag_real(res, projection, b_0, 'disc', jacobi=jacobi))
+result_during = np.asarray(np.asmatrix(jacobi)*np.asmatrix(vector).T).reshape(-1)
 toc = time.clock()
-print 'Time for getting the transp. Jacobi-Matrix (vector-wise):', toc - tic
+print 'Calculation during phasemapping:  ', toc - tic
 
-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))
+# Test if calculations match:
+np.testing.assert_almost_equal(result_kernel, result_during)
+np.testing.assert_almost_equal(result_core, result_during)
+np.testing.assert_almost_equal(result_core2, result_during)
+print 'Calculations give the same result!'