From 48b5c879ec0b63685e06ff56b2ccc2317ae36a21 Mon Sep 17 00:00:00 2001
From: Jan Caron <j.caron@fz-juelich.de>
Date: Thu, 5 Dec 2013 18:54:24 +0100
Subject: [PATCH] Implementation of jacobi multiplication for the projection.
(At the moment a bit slow, optimization with Cython seems necessary).
---
pyramid/kernel.py | 24 ++-
pyramid/numcore/phase_mag_real.pyx | 46 ++---
pyramid/phasemapper.py | 16 --
pyramid/projector.py | 288 +++++++++++++++++++--------
pyramid/reconstructor.py | 6 +
scripts/get_jacobi_projection.py | 136 ++++++++++++-
scripts/reconstruction/get_jacobi.py | 55 +++--
7 files changed, 413 insertions(+), 158 deletions(-)
diff --git a/pyramid/kernel.py b/pyramid/kernel.py
index 51bf156..89bdac3 100644
--- a/pyramid/kernel.py
+++ b/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
diff --git a/pyramid/numcore/phase_mag_real.pyx b/pyramid/numcore/phase_mag_real.pyx
index 2fa38de..4171ce5 100644
--- a/pyramid/numcore/phase_mag_real.pyx
+++ b/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)
diff --git a/pyramid/phasemapper.py b/pyramid/phasemapper.py
index d74495e..d69f262 100644
--- a/pyramid/phasemapper.py
+++ b/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
----------
diff --git a/pyramid/projector.py b/pyramid/projector.py
index 9fa94d2..3c5fab3 100644
--- a/pyramid/projector.py
+++ b/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
diff --git a/pyramid/reconstructor.py b/pyramid/reconstructor.py
index 7c995dc..5196cd0 100644
--- a/pyramid/reconstructor.py
+++ b/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
diff --git a/scripts/get_jacobi_projection.py b/scripts/get_jacobi_projection.py
index 079c39a..ecff0f2 100644
--- a/scripts/get_jacobi_projection.py
+++ b/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
diff --git a/scripts/reconstruction/get_jacobi.py b/scripts/reconstruction/get_jacobi.py
index c491900..b963050 100644
--- a/scripts/reconstruction/get_jacobi.py
+++ b/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!'
--
GitLab