From ebd9398e18ee090a3e3c8f6aa36ef2aab595af09 Mon Sep 17 00:00:00 2001 From: Jan Caron <j.caron@fz-juelich.de> Date: Sat, 23 Nov 2013 19:46:25 +0100 Subject: [PATCH] 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 | 30 ++ pyramid/numcore/phase_mag_real.pyx | 47 ++++ pyramid/phasemapper.py | 265 ++++++++++++------ .../compare methods/compare_pixel_fields.py | 19 +- .../create_core_shell_disc.py | 2 +- scripts/reconstruction/get_jacobi.py | 48 +++- .../reconstruction/reconstruct_core_shell.py | 1 - scripts/test_h5py.py | 37 +++ scripts/test_methods.py | 3 +- scripts/test_vadim.py | 187 ++++++------ tests/test_compliance.py | 21 +- 11 files changed, 476 insertions(+), 184 deletions(-) create mode 100644 scripts/test_h5py.py diff --git a/pyramid/magcreator.py b/pyramid/magcreator.py index 19fc810..11eea24 100644 --- a/pyramid/magcreator.py +++ b/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. diff --git a/pyramid/numcore/phase_mag_real.pyx b/pyramid/numcore/phase_mag_real.pyx index a849783..d679933 100644 --- a/pyramid/numcore/phase_mag_real.pyx +++ b/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) diff --git a/pyramid/phasemapper.py b/pyramid/phasemapper.py index 58413ea..8ddf90b 100644 --- a/pyramid/phasemapper.py +++ b/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): diff --git a/scripts/compare methods/compare_pixel_fields.py b/scripts/compare methods/compare_pixel_fields.py index 414ac02..a3f6805 100644 --- a/scripts/compare methods/compare_pixel_fields.py +++ b/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) diff --git a/scripts/create distributions/create_core_shell_disc.py b/scripts/create distributions/create_core_shell_disc.py index 0242008..41d1231 100644 --- a/scripts/create distributions/create_core_shell_disc.py +++ b/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 diff --git a/scripts/reconstruction/get_jacobi.py b/scripts/reconstruction/get_jacobi.py index f3369d5..87a9b47 100644 --- a/scripts/reconstruction/get_jacobi.py +++ b/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__": diff --git a/scripts/reconstruction/reconstruct_core_shell.py b/scripts/reconstruction/reconstruct_core_shell.py index 1cb4c37..56d2c03 100644 --- a/scripts/reconstruction/reconstruct_core_shell.py +++ b/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) diff --git a/scripts/test_h5py.py b/scripts/test_h5py.py new file mode 100644 index 0000000..40271b2 --- /dev/null +++ b/scripts/test_h5py.py @@ -0,0 +1,37 @@ +# -*- 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:]) diff --git a/scripts/test_methods.py b/scripts/test_methods.py index be89108..f52fcfe 100644 --- a/scripts/test_methods.py +++ b/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)) diff --git a/scripts/test_vadim.py b/scripts/test_vadim.py index 55d6e8d..ce1f3d5 100644 --- a/scripts/test_vadim.py +++ b/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') diff --git a/tests/test_compliance.py b/tests/test_compliance.py index a0b651d..05f707a 100644 --- a/tests/test_compliance.py +++ b/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 -- GitLab