From cfea0d77947c3b0b682a39a2aa8a7732bbb22d25 Mon Sep 17 00:00:00 2001
From: Jan Caron <j.caron@fz-juelich.de>
Date: Fri, 28 Nov 2014 13:59:27 +0100
Subject: [PATCH] Implementation of FFTW and covariance matrix functionality
added Se_inv: is now handled differently, defaults to eye-matrix, but can be
set later (added a few convenience functions to DataSet) utility:
(re-)added utility module for index conversion kernel/phasemapper: can now
use FFTW library and operate on pre-allocated arrays dtype: set dtype
of arrays to np.float32/np.complex64 (halves memory size), also
applied to numcore routines! magdata: quiver_plot can now plot every n-th
arrow (yet untested!) reconstruction scripts: modified several reconstruction
scripts.
---
pyramid/__init__.py | 2 +
pyramid/costfunction.py | 6 +-
pyramid/dataset.py | 31 ++--
pyramid/kernel.py | 55 +++---
pyramid/magdata.py | 9 +-
pyramid/numcore/phasemapper_core.pyx | 20 +--
pyramid/phasemap.py | 6 +-
pyramid/phasemapper.py | 168 +++++++-----------
pyramid/reconstruction.py | 2 +-
pyramid/util.py | 131 ++++++++++++++
pyramid/version.py | 2 +-
.../patrick_nanowire_reconstruction.py | 60 +++++--
.../reconstruction_linear_test.py | 16 +-
.../reconstruction_linear_test_batch.py | 45 +++--
scripts/test methods/compare_methods.py | 35 ++--
15 files changed, 376 insertions(+), 212 deletions(-)
create mode 100644 pyramid/util.py
diff --git a/pyramid/__init__.py b/pyramid/__init__.py
index aed6a3d..dace4b4 100644
--- a/pyramid/__init__.py
+++ b/pyramid/__init__.py
@@ -51,6 +51,7 @@ from .phasemap import *
from .phasemapper import *
from .projector import *
from .regularisator import *
+from .util import *
from .version import version as __version__
from .version import hg_revision as __hg_revision__
@@ -68,3 +69,4 @@ __all__.extend(phasemap.__all__)
__all__.extend(phasemapper.__all__)
__all__.extend(projector.__all__)
__all__.extend(regularisator.__all__)
+__all__.extend(util.__all__)
diff --git a/pyramid/costfunction.py b/pyramid/costfunction.py
index 0859ae3..9f36779 100644
--- a/pyramid/costfunction.py
+++ b/pyramid/costfunction.py
@@ -5,6 +5,7 @@ the so called `cost` of a threedimensional magnetization distribution."""
import numpy as np
+from scipy.sparse import eye as sparse_eye
from scipy.sparse.linalg import LinearOperator
from pyramid.forwardmodel import ForwardModel
@@ -60,9 +61,12 @@ class Costfunction(object):
self.regularisator = NoneRegularisator()
# Extract important information:
self.y = data_set.phase_vec
- self.Se_inv = data_set.Se_inv
self.n = data_set.n
self.m = data_set.m
+ if data_set.Se_inv is not None:
+ self.Se_inv = data_set.Se_inv
+ else:
+ self.Se_inv = sparse_eye(self.m)
self._log.debug('Created '+str(self))
def __repr__(self):
diff --git a/pyramid/dataset.py b/pyramid/dataset.py
index 9559d59..1e3ad43 100644
--- a/pyramid/dataset.py
+++ b/pyramid/dataset.py
@@ -6,7 +6,7 @@ and additional data like corresponding projectors."""
import numpy as np
from numbers import Number
-from scipy.sparse import eye as sparse_eye
+from scipy import sparse
import matplotlib.pyplot as plt
@@ -62,11 +62,6 @@ class DataSet(object):
def m(self):
return np.sum([len(p.phase_vec) for p in self.phase_maps])
- @property
- def Se_inv(self):
- # TODO: better implementation, maybe get-method? more flexible? input in append?
- return sparse_eye(self.m)
-
@property
def phase_vec(self):
return np.concatenate([p.phase_vec for p in self.phase_maps])
@@ -80,11 +75,13 @@ class DataSet(object):
@property
def phase_mappers(self):
- dim_uv_list = np.unique([p.dim_uv for p in self.projectors])
- kernel_list = [Kernel(self.a, tuple(dim_uv)) for dim_uv in dim_uv_list]
+ dim_uv_set = set([p.dim_uv for p in self.projectors])
+ kernel_list = [Kernel(self.a, dim_uv, use_fftw=self.use_fftw, threads=self.threads)
+ for dim_uv in dim_uv_set]
return {kernel.dim_uv: PhaseMapperRDFC(kernel) for kernel in kernel_list}
- def __init__(self, a, dim, b_0=1, mask=None):
+ def __init__(self, a, dim, b_0=1, mask=None, Se_inv=None, use_fftw=True, threads=1):
+ # TODO: document!
self._log.debug('Calling __init__')
assert isinstance(a, Number), 'Grid spacing has to be a number!'
assert a >= 0, 'Grid spacing has to be a positive number!'
@@ -99,6 +96,9 @@ class DataSet(object):
self.dim = dim
self.b_0 = b_0
self.mask = mask
+ self.Se_inv = Se_inv
+ self.use_fftw = use_fftw
+ self.threads = threads
self.phase_maps = []
self.projectors = []
self._log.debug('Created: '+str(self))
@@ -111,7 +111,7 @@ class DataSet(object):
self._log.debug('Calling __str__')
return 'DataSet(a=%s, dim=%s, b_0=%s)' % (self.a, self.dim, self.b_0)
- def append(self, phase_map, projector): # TODO: include Se_inv or 2D mask??
+ def append(self, phase_map, projector):
'''Appends a data pair of phase map and projection infos to the data collection.`
Parameters
@@ -152,6 +152,15 @@ class DataSet(object):
return [self.phase_mappers[projector.dim_uv](projector(mag_data))
for projector in self.projectors]
+ def set_Se_inv_block_diag(self, cov_list):
+ # TODO: Document!
+ self.Se_inv = sparse.block_diag(cov_list).tocsr()
+
+ def set_Se_inv_diag_with_masks(self, mask_list):
+ # TODO: Document!
+ cov_list = [sparse.diags(m.flatten(), 0) for m in mask_list]
+ self.set_Se_inv_block_diag(cov_list)
+
def display_phase(self, mag_data=None, title='Phase Map',
cmap='RdBu', limit=None, norm=None):
'''Display all phasemaps saved in the :class:`~.DataSet` as a colormesh.
@@ -235,5 +244,3 @@ class DataSet(object):
cmap, limit, norm, gain, interpolation, grad_encode)
for (i, phase_map) in enumerate(phase_maps)]
plt.show()
-
-# TODO: method for constructing 3D mask from 2D masks?
diff --git a/pyramid/kernel.py b/pyramid/kernel.py
index cb237b0..c2b1664 100644
--- a/pyramid/kernel.py
+++ b/pyramid/kernel.py
@@ -53,7 +53,7 @@ class Kernel(object):
A tuple of :class:`slice` objects to extract the original field of view from the increased
size (`size_fft`) of the grid for the FFT-convolution.
- ''' # TODO: overview what all dim_??? mean! and use_fftw
+ ''' # TODO: overview what all dim_??? mean! and use_fftw, slice(_fft), etc.
_log = logging.getLogger(__name__+'.Kernel')
@@ -66,32 +66,36 @@ class Kernel(object):
self.a = a
self.geometry = geometry
# Set up FFT:
- if use_fftw:
+ if not use_fftw:
+ self.dim_pad = tuple(2**np.ceil(np.log2(2*np.array(dim_uv))).astype(int)) # pow(2)
+ self.use_fftw = False
+ else:
try:
import pyfftw
+ self.dim_pad = tuple(2*np.array(dim_uv)) # is at least even (not nec. power of 2)
+ self.use_fftw = True
except ImportError:
- use_fftw = False
+ self.dim_pad = tuple(2**np.ceil(np.log2(2*np.array(dim_uv))).astype(int)) # pow(2)
+ self.use_fftw = False
self._log.info('pyFFTW could not be imported, using numpy instead!')
- if use_fftw: # use pyfftw (FFTW wrapper for python)
- self.dim_pad = tuple(2*np.array(dim_uv)) # is at least even (not nec. power of 2)
- self.dim_fft = (self.dim_pad[0], self.dim_pad[1]/2.+1) # last axis is real
- n = pyfftw.simd_alignment
- self.u = pyfftw.n_byte_align_empty(self.dim_kern, n, dtype='float32')
- self.v = pyfftw.n_byte_align_empty(self.dim_kern, n, dtype='float32')
- self.u_fft = pyfftw.n_byte_align_empty(self.dim_fft, n, dtype='complex64')
- self.v_fft = pyfftw.n_byte_align_empty(self.dim_fft, n, dtype='complex64')
+ self.dim_fft = (self.dim_pad[0], self.dim_pad[1]/2+1) # last axis is real
+ self.slice_phase = (slice(dim_uv[0]-1, self.dim_kern[0]), # Shift because kernel center
+ slice(dim_uv[1]-1, self.dim_kern[1])) # is not at (0, 0)!
+ self.slice_mag = (slice(0, dim_uv[0]), # Magnetization is padded on the far end!
+ slice(0, dim_uv[1])) # (Phase cutout is shifted as listed above)
+ self.u = np.zeros(self.dim_kern, dtype=np.float32)
+ self.v = np.zeros(self.dim_kern, dtype=np.float32)
+ self.u_fft = np.zeros(self.dim_fft, dtype=np.complex64)
+ self.v_fft = np.zeros(self.dim_fft, dtype=np.complex64)
+ if self.use_fftw: # use pyfftw (FFTW wrapper for python)
+ self.u = pyfftw.n_byte_align(self.u, pyfftw.simd_alignment)
+ self.v = pyfftw.n_byte_align(self.v, pyfftw.simd_alignment)
+ self.u_fft = pyfftw.n_byte_align(self.u_fft, pyfftw.simd_alignment)
+ self.v_fft = pyfftw.n_byte_align(self.v_fft, pyfftw.simd_alignment)
rfftn = pyfftw.builders.rfftn(self.u, self.dim_pad, threads=threads)
self.threads = threads
- self.use_fftw = True
else: # otherwise use numpy
- self.dim_pad = tuple(2**np.ceil(np.log2(2*np.array(dim_uv))).astype(int)) # pow(2)
- self.dim_fft = (self.dim_pad[0], self.dim_pad[1]/2+1) # last axis is real
- self.u = np.empty(self.dim_kern, dtype=float)
- self.v = np.empty(self.dim_kern, dtype=float)
- self.u_fft = np.empty(self.dim_fft, dtype=complex)
- self.v_fft = np.empty(self.dim_fft, dtype=complex)
rfftn = lambda x: np.fft.rfftn(x, self.dim_pad)
- self.use_fftw = False
# Calculate kernel (single pixel phase):
coeff = b_0 * a**2 / (2*PHI_0) # Minus is gone because of negative z-direction
v_dim, u_dim = dim_uv
@@ -101,7 +105,6 @@ class Kernel(object):
self.u[...] = coeff * self._get_elementary_phase(geometry, uu, vv, a)
self.v[...] = coeff * self._get_elementary_phase(geometry, vv, uu, a)
# Calculate Fourier trafo of kernel components:
- self.slice_fft = (slice(dim_uv[0]-1, 2*dim_uv[0]-1), slice(dim_uv[1]-1, 2*dim_uv[1]-1))
self.u_fft[...] = rfftn(self.u)
self.v_fft[...] = rfftn(self.v)
self._log.debug('Created '+str(self))
@@ -137,5 +140,13 @@ class Kernel(object):
return 0.5 * (F_a(n-0.5, m-0.5) - F_a(n+0.5, m-0.5)
- F_a(n-0.5, m+0.5) + F_a(n+0.5, m+0.5))
- def get_info(self):
- pass
+ def print_info(self):
+ print 'Shape of the FOV :', self.dim_uv
+ print 'Shape of the Kernel:', self.dim_kern
+ print 'Zero-padded shape :', self.dim_pad
+ print 'Shape of the FFT :', self.dim_fft
+ print 'Slice for the phase:', self.slice_phase
+ print 'Slice for the magn.:', self.slice_mag
+ print 'Grid spacing: {} nm'.format(self.a)
+ print 'Geometry:', self.geometry
+ print 'Use FFTW: {}; with {} thread(s)'.format(self.use_fftw, self.threads)
diff --git a/pyramid/magdata.py b/pyramid/magdata.py
index 0c8e2fc..f24686d 100644
--- a/pyramid/magdata.py
+++ b/pyramid/magdata.py
@@ -73,7 +73,7 @@ class MagData(object):
assert isinstance(magnitude, np.ndarray), 'Magnitude has to be a numpy array!'
assert len(magnitude.shape) == 4, 'Magnitude has to be 4-dimensional!'
assert magnitude.shape[0] == 3, 'First dimension of the magnitude has to be 3!'
- self._magnitude = magnitude
+ self._magnitude = np.asarray(magnitude, dtype=np.float32)
self._dim = magnitude.shape[1:]
@property
@@ -419,7 +419,7 @@ class MagData(object):
return MagData(a, magnitude)
def quiver_plot(self, title='Magnetization Distribution', axis=None, proj_axis='z',
- ax_slice=None, log=False, scaled=True):
+ ar_dens=1, ax_slice=None, log=False, scaled=True): # TODO: Doc ar_dens
'''Plot a slice of the magnetization as a quiver plot.
Parameters
@@ -494,8 +494,9 @@ class MagData(object):
if scaled:
plot_u /= np.hypot(plot_u, plot_v).max()
plot_v /= np.hypot(plot_u, plot_v).max()
- axis.quiver(plot_u, plot_v, pivot='middle', units='xy', angles=angles,
- scale_units='xy', scale=1, headwidth=6, headlength=7)
+ ad = ar_dens
+ axis.quiver(plot_u[::ad, ::ad], plot_v[::ad, ::ad], pivot='middle', units='xy',
+ angles=angles, scale_units='xy', scale=1, headwidth=6, headlength=7)
axis.set_xlim(-1, np.shape(plot_u)[1])
axis.set_ylim(-1, np.shape(plot_u)[0])
axis.set_title(title, fontsize=18)
diff --git a/pyramid/numcore/phasemapper_core.pyx b/pyramid/numcore/phasemapper_core.pyx
index 7ea34d9..d7b1f04 100644
--- a/pyramid/numcore/phasemapper_core.pyx
+++ b/pyramid/numcore/phasemapper_core.pyx
@@ -19,9 +19,9 @@ cimport numpy as np
@cython.wraparound(False)
def phasemapper_real_convolve(
unsigned int v_dim, unsigned int u_dim,
- double[:, :] v_phi, double[:, :] u_phi,
- double[:, :] v_mag, double[:, :] u_mag,
- double[:, :] phase, float threshold):
+ float[:, :] v_phi, float[:, :] u_phi,
+ float[:, :] v_mag, float[:, :] u_mag,
+ float[:, :] phase, float threshold):
'''Numerical core routine for the phase calculation using the real space approach.
Parameters
@@ -43,7 +43,7 @@ def phasemapper_real_convolve(
'''
cdef unsigned int i, j, p, q, p_c, q_c
- cdef double u_m, v_m
+ cdef float u_m, v_m
for j in range(v_dim):
for i in range(u_dim):
u_m = u_mag[j, i]
@@ -64,9 +64,9 @@ def phasemapper_real_convolve(
@cython.wraparound(False)
def jac_dot_real_convolve(
unsigned int v_dim, unsigned int u_dim,
- double[:, :] u_phi, double[:, :] v_phi,
- double[:] vector,
- double[:] result):
+ float[:, :] u_phi, float[:, :] v_phi,
+ float[:] vector,
+ float[:] result):
'''Numerical core routine for the Jacobi matrix multiplication for the phase mapper.
Parameters
@@ -113,9 +113,9 @@ def jac_dot_real_convolve(
@cython.wraparound(False)
def jac_T_dot_real_convolve(
unsigned int v_dim, unsigned int u_dim,
- double[:, :] u_phi, double[:, :] v_phi,
- double[:] vector,
- double[:] result):
+ float[:, :] u_phi, float[:, :] v_phi,
+ float[:] vector,
+ float[:] result):
'''Core routine for the transposed Jacobi multiplication for the phase mapper.
Parameters
diff --git a/pyramid/phasemap.py b/pyramid/phasemap.py
index b33ab8c..4b962dd 100644
--- a/pyramid/phasemap.py
+++ b/pyramid/phasemap.py
@@ -121,7 +121,7 @@ class PhaseMap(object):
def phase(self, phase):
assert isinstance(phase, np.ndarray), 'Phase has to be a numpy array!'
assert len(phase.shape) == 2, 'Phase has to be 2-dimensional!'
- self._phase = phase
+ self._phase = np.asarray(phase, dtype=np.float32)
self._dim_uv = phase.shape
@property
@@ -344,8 +344,8 @@ class PhaseMap(object):
axis.set_ylim(0, self.dim_uv[0])
axis.xaxis.set_major_locator(MaxNLocator(nbins=9, integer=True))
axis.yaxis.set_major_locator(MaxNLocator(nbins=9, integer=True))
- axis.xaxis.set_major_formatter(FuncFormatter(lambda x, pos: '{:g}'.format(x*self.a)))
- axis.yaxis.set_major_formatter(FuncFormatter(lambda x, pos: '{:g}'.format(x*self.a)))
+ axis.xaxis.set_major_formatter(FuncFormatter(lambda x, pos: '{:.3g}'.format(x*self.a)))
+ axis.yaxis.set_major_formatter(FuncFormatter(lambda x, pos: '{:.3g}'.format(x*self.a)))
axis.tick_params(axis='both', which='major', labelsize=14)
axis.set_title(title, fontsize=18)
axis.set_xlabel('u-axis [nm]', fontsize=15)
diff --git a/pyramid/phasemapper.py b/pyramid/phasemapper.py
index 2c2c0e9..32387dd 100644
--- a/pyramid/phasemapper.py
+++ b/pyramid/phasemapper.py
@@ -100,32 +100,29 @@ class PhaseMapperRDFC(PhaseMapper):
self.kernel = kernel
self.m = np.prod(kernel.dim_uv)
self.n = 2 * self.m
+ self.u_mag = np.zeros(kernel.dim_pad, dtype=np.float32)
+ self.v_mag = np.zeros(kernel.dim_pad, dtype=np.float32)
+ self.mag_adj = np.zeros(kernel.dim_pad, dtype=np.float32)
+ self.u_mag_fft = np.zeros(kernel.dim_fft, dtype=np.complex64)
+ self.v_mag_fft = np.zeros(kernel.dim_fft, dtype=np.complex64)
+ self.phase_fft = np.zeros(kernel.dim_fft, dtype=np.complex64)
if self.kernel.use_fftw: # if possible use FFTW
import pyfftw
- n = pyfftw.simd_alignment
- self.u_mag = pyfftw.n_byte_align_empty(self.kernel.dim_uv, n, dtype='float32')
- self.v_mag = pyfftw.n_byte_align_empty(self.kernel.dim_uv, n, dtype='float32')
-# self.phase = pyfftw.n_byte_align_empty(self.kernel.dim_uv, n, dtype='float32')
- self.u_mag_fft = pyfftw.n_byte_align_empty(self.kernel.dim_fft, n, dtype='complex64')
- self.v_mag_fft = pyfftw.n_byte_align_empty(self.kernel.dim_fft, n, dtype='complex64')
- self.phase_fft = pyfftw.n_byte_align_empty(self.kernel.dim_fft, n, dtype='complex64')
- self._rfft2 = pyfftw.builders.rfftn(self.u_mag, self.kernel.dim_pad,
- threads=self.kernel.threads)
- self._irfft2 = pyfftw.builders.irfftn(self.phase_fft, self.kernel.dim_pad,
- threads=self.kernel.threads)
+ self.u_mag = pyfftw.n_byte_align(self.u_mag, pyfftw.simd_alignment)
+ self.v_mag = pyfftw.n_byte_align(self.v_mag, pyfftw.simd_alignment)
+ self.mag_adj = pyfftw.n_byte_align(self.mag_adj, pyfftw.simd_alignment)
+ self.u_mag_fft = pyfftw.n_byte_align(self.u_mag_fft, pyfftw.simd_alignment)
+ self.v_mag_fft = pyfftw.n_byte_align(self.v_mag_fft, pyfftw.simd_alignment)
+ self.phase_fft = pyfftw.n_byte_align(self.phase_fft, pyfftw.simd_alignment)
+ self._rfft2 = pyfftw.builders.rfftn(self.u_mag, threads=kernel.threads)
+ self._irfft2 = pyfftw.builders.irfftn(self.phase_fft, threads=kernel.threads)
self._rfft2_adj = self._rfft2_adj_fftw
self._irfft2_adj = self._irfft2_adj_fftw
else: # otherwise use numpy
- self.u_mag = np.empty(self.kernel.dim_uv, dtype=float)
- self.v_mag = np.empty(self.kernel.dim_uv, dtype=float)
-# self.phase = np.empty(self.kernel.dim_uv, dtype=float)
- self.u_mag_fft = np.empty(self.kernel.dim_fft, dtype=complex)
- self.v_mag_fft = np.empty(self.kernel.dim_fft, dtype=complex)
- self.phase_fft = np.empty(self.kernel.dim_fft, dtype=complex)
- self._rfft2 = lambda x: np.fft.rfftn(x, self.kernel.dim_pad)
- self._irfft2 = lambda x: np.fft.irfftn(x, self.kernel.dim_pad)
- self._rfft2_adj = jfft.irfft2_adj
- self._rfft2_adj = jfft.irfft2_adj
+ self._rfft2 = np.fft.rfftn
+ self._irfft2 = np.fft.irfftn
+ self._rfft2_adj = lambda x: jfft.rfft2_adj(x, self.kernel.dim_pad[1])
+ self._irfft2_adj = jfft.irfft2_adj
self._log.debug('Created '+str(self))
def __repr__(self):
@@ -143,7 +140,8 @@ class PhaseMapperRDFC(PhaseMapper):
assert mag_data.dim[0] == 1, 'Magnetic distribution must be 2-dimensional!'
assert mag_data.dim[1:3] == self.kernel.dim_uv, 'Dimensions do not match!'
# Process input parameters:
- self.u_mag[...], self.v_mag[...] = mag_data.magnitude[0:2, 0, ...]
+ self.u_mag[self.kernel.slice_mag] = mag_data.magnitude[0, 0, ...] # u-component
+ self.v_mag[self.kernel.slice_mag] = mag_data.magnitude[1, 0, ...] # v-component
return PhaseMap(mag_data.a, self._convolve())
def _convolve(self):
@@ -153,7 +151,7 @@ class PhaseMapperRDFC(PhaseMapper):
# Convolve the magnetization with the kernel in Fourier space:
self.phase_fft[...] = self.u_mag_fft*self.kernel.u_fft - self.v_mag_fft*self.kernel.v_fft
# Return the result:
- return self._irfft2(self.phase_fft)[self.kernel.slice_fft]
+ return self._irfft2(self.phase_fft)[self.kernel.slice_phase]
def jac_dot(self, vector):
'''Calculate the product of the Jacobi matrix with a given `vector`.
@@ -174,8 +172,9 @@ class PhaseMapperRDFC(PhaseMapper):
# self._log.debug('Calling jac_dot') # TODO: Profiler says this was slow...
assert len(vector) == self.n, \
'vector size not compatible! vector: {}, size: {}'.format(len(vector), self.n)
- self.u_mag[...], self.v_mag[...] = np.reshape(vector, (2,)+self.kernel.dim_uv)
- result = self._convolve().reshape(-1)
+ self.u_mag[self.kernel.slice_mag], self.v_mag[self.kernel.slice_mag] = \
+ np.reshape(vector, (2,)+self.kernel.dim_uv)
+ result = self._convolve().flatten()
return result
def jac_T_dot(self, vector):
@@ -197,81 +196,50 @@ class PhaseMapperRDFC(PhaseMapper):
# self._log.debug('Calling jac_T_dot') # TODO: Profiler says this was slow...
assert len(vector) == self.m, \
'vector size not compatible! vector: {}, size: {}'.format(len(vector), self.m)
- phase = vector.reshape(self.kernel.dim_uv)
- p0 = self.kernel.dim_pad[0]-self.kernel.dim_uv[0]
- p1 = self.kernel.dim_pad[1]-self.kernel.dim_uv[1]
- phase = np.pad(phase, ((0, p0), (0, p1)), 'constant')
- phase_fft = jfft.irfft2_adj(phase)
- u_mag_fft = phase_fft * self.kernel.u_fft
- v_mag_fft = phase_fft * -self.kernel.v_fft
- u_mag = jfft.rfft2_adj(u_mag_fft, self.kernel.dim_pad[1])[self.kernel.slice_fft]
- v_mag = jfft.rfft2_adj(v_mag_fft, self.kernel.dim_pad[1])[self.kernel.slice_fft]
- return -np.concatenate((u_mag.flatten(), v_mag.flatten())) # TODO: why the minus?
-
-# TODO: save arrays in PADDED form (kernel and the self.u/v/... stuff) to save padding time
-# TODO: also gets rid of (..., self.kernel.dim_pad) argument
-# if a.shape[axis] != n:
-# s = list(a.shape)
-# if s[axis] > n:
-# index = [slice(None)]*len(s)
-# index[axis] = slice(0, n)
-# a = a[index]
-# else:
-# index = [slice(None)]*len(s)
-# index[axis] = slice(0, s[axis])
-# s[axis] = n
-# z = zeros(s, a.dtype.char)
-# z[index] = a
-# a = z
-
- def _rfft2_adj_fftw(x, n=None):
- """
- Adjoint of the 2-D Real Fast Fourier Transform, that is the product of the adjoint
- Jacobian of the numpy rfftn with vector x.
-
- Parameters
- ----------
- x : array_like
- n : int, optional
- Length of second dimension of original array (the one, the size of which is
- changed by employing the rfft2)
-
- Returns
- -------
- array_like
- """
- print 'rfftn2_adj input:', x.shape
- if n is None:
- n = 2 * (x.shape[1] - 1)
- xp = np.zeros((x.shape[0], n), dtype=x.dtype)
+ self.mag_adj[self.kernel.slice_mag] = vector.reshape(self.kernel.dim_uv)
+ mag_adj_fft = self._irfft2_adj(self.mag_adj)
+ self.u_mag_fft[...] = mag_adj_fft * self.kernel.u_fft # u/v_mag_fft is just used as cache
+ self.v_mag_fft[...] = mag_adj_fft * -self.kernel.v_fft # they are really 'adjoint' phases
+ u_phase_adj = self._rfft2_adj(self.u_mag_fft)[self.kernel.slice_phase]
+ v_phase_adj = self._rfft2_adj(self.v_mag_fft)[self.kernel.slice_phase]
+ return -np.concatenate((u_phase_adj.flatten(), v_phase_adj.flatten())) # TODO: why minus?
+
+ def _rfft2_adj_fftw(self, phase_adj):
+ x = phase_adj
+ xp = np.zeros(self.kernel.dim_pad, dtype=np.complex64)
xp[:, :x.shape[1]] = x
- print 'rfftn2_adj output:', np.fft.ifftn(xp).real * x.shape[0] * n
- return np.fft.ifftn(xp).real * x.shape[0] * n
-
- def _irfft2_adj_fftw(x):
- """
- Adjoint of the 2-D Inverse Real Fast Fourier Transform, that is the product of the
- adjoint Jacobian of the numpy irfftn with vector x.
-
- Parameters
- ----------
- x : array_like
-
- Returns
- -------
- array_like
- """
- print 'irfftn2_adj input:', x.shape
- n_out = x.shape[1] / 2 + 1
- xp = np.fft.fft(x, axis=1) / x.shape[1]
- if x.shape[1] % 2 == 0:
- xp[:, 1:n_out - 1] += np.conj(xp[:, :n_out-1:-1])
- # xp[:, n_out - 1] = xp[:, n_out - 1].real
- else:
- xp[:, 1:n_out] += np.conj(xp[:, :n_out-1:-1])
- # xp[:, 0] = xp[:, 0].real
- print 'rfftn2_adj output:', np.fft.ifftn(xp).real * x.shape[0] * n
- return np.fft.fft(xp[:, :n_out], axis=0) / xp.shape[0]
+ phase_adj_fft = np.fft.ifftn(xp).real * np.prod(self.kernel.dim_pad) # TODO: Convert to FFTW
+ return phase_adj_fft
+ import pyfftw
+ kernel = self.kernel
+ self.phase_adj = np.zeros(kernel.dim_pad, dtype=np.complex64)
+ self.phase_adj = pyfftw.n_byte_align(self.phase_adj, pyfftw.simd_alignment)
+ self.phase_adj[:, :phase_adj.shape[1]] = phase_adj
+ self._ifft2 = pyfftw.builders.ifftn(self.u_mag, threads=kernel.threads)
+ phase_adj_fft_TEST = self._ifft2(self.phase_adj).real * np.prod(self.kernel.dim_pad)
+ import pdb; pdb.set_trace()
+ return phase_adj_fft
+
+ def _irfft2_adj_fftw(self, mag_adj_fft):
+ x = mag_adj_fft
+ n_out = self.kernel.dim_fft[1]
+ xp = np.zeros(self.kernel.dim_pad, dtype=np.complex64)
+ xp[...] = np.fft.fft(x, axis=1) / self.kernel.dim_pad[1] # TODO: Convert to FFTW
+ xp[:, 1:n_out - 1] += np.conj(xp[:, :n_out-1:-1]) # if even (guaranteed for dim_pad)
+ mag_adj = np.fft.fft(xp[:, :n_out], axis=0) / self.kernel.dim_pad[0] # TODO: Convert to FFTW
+ return mag_adj
+ import pyfftw
+ kernel = self.kernel
+ self.mag_adj = np.zeros(kernel.dim_pad, dtype=np.complex64)
+ self.mag_adj = pyfftw.n_byte_align(self.mag_adj, pyfftw.simd_alignment)
+ self._fft_ax0 = pyfftw.FFTW(self.mag_adj, threads=kernel.threads)
+ self._fft_ax1 = pyfftw.FFTW(self.mag_adj, threads=kernel.threads)
+ n_out = self.kernel.dim_fft[1]
+ self.mag_adj[...] = np.fft.fft(phase_adj, axis=1) / self.kernel.dim_pad[1] # TODO: Convert to FFTW
+ self.mag_adj[:, 1:n_out - 1] += np.conj(self.mag_adj[:, :n_out-1:-1]) # if even (guaranteed for dim_pad)
+ mag_adj_TEST = np.fft.fft(self.mag_adj[:, :n_out], axis=0) / self.kernel.dim_pad[0] # TODO: Convert to FFTW
+ import pdb; pdb.set_trace()
+ return mag_adj
class PhaseMapperRDRC(PhaseMapper):
@@ -339,7 +307,7 @@ class PhaseMapperRDRC(PhaseMapper):
u_phi = self.kernel.u
v_phi = self.kernel.v
# Calculation of the phase:
- phase = np.zeros(dim_uv)
+ phase = np.zeros(dim_uv, dtype=np.float32)
if self.numcore:
nc.phasemapper_real_convolve(dim_uv[0], dim_uv[1], v_phi, u_phi,
v_mag, u_mag, phase, self.threshold)
diff --git a/pyramid/reconstruction.py b/pyramid/reconstruction.py
index 62812db..4da85c7 100644
--- a/pyramid/reconstruction.py
+++ b/pyramid/reconstruction.py
@@ -103,7 +103,7 @@ def optimize_linear(data, regularisator=None, max_iter=None, info=None):
mag_data : :class:`~pyramid.magdata.MagData`
The reconstructed magnetic distribution as a :class:`~.MagData` object.
- '''
+ ''' # TODO: info document!
import jutil.cg as jcg
_log.debug('Calling optimize_linear')
# Set up necessary objects:
diff --git a/pyramid/util.py b/pyramid/util.py
new file mode 100644
index 0000000..40ce45a
--- /dev/null
+++ b/pyramid/util.py
@@ -0,0 +1,131 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Wed Nov 26 10:52:32 2014
+
+@author: Jan
+"""
+
+
+import numpy as np
+
+
+__all__ = ['IndexConverter3D', 'IndexConverter2D']
+
+
+class IndexConverter3D(object):
+
+ # TODO: Document everything! 3comp: 3 components!
+
+ def __init__(self, dim):
+ assert len(dim) == 3, 'Dimensions have to be two-dimensioal!'
+ self.dim = dim
+ self.size_3d = np.prod(dim)
+ self.size_2d = dim[2]*dim[1]
+
+ def ind_to_coord(self, ind):
+ assert ind < self.size_3d, 'Index out of range!'
+ z, remain = int(ind/self.size_2d), ind % self.size_2d
+ y, x = int(remain/self.dim[2]), remain % self.dim[2]
+ return (z, y, x)
+
+ def ind3comp_to_coord(self, ind):
+ assert ind < 3 * self.size_3d, 'Index out of range!'
+ c, remain = int(ind/self.size_3d), ind % self.size_3d
+ z, remain = int(remain/self.size_2d), remain % self.size_2d
+ y, x = int(remain/self.dim[2]), remain % self.dim[2]
+ return (c, z, y, x)
+
+ def coord_to_ind(self, coord):
+ if coord is None:
+ return None
+ z, y, x = coord
+ ind = z*self.size_2d + y*self.dim[2] + x
+ assert ind < self.size_3d, 'Index out of range!'
+ return ind
+
+ def coord_to_ind3comp(self, coord):
+ if coord is None:
+ return [None, None, None]
+ z, y, x = coord
+ ind = [c*self.size_3d + z*self.size_2d + y*self.dim[2] + x for c in range(3)]
+ assert ind[-1] < 3 * self.size_3d, 'Index out of range!'
+ return ind
+
+ def get_neighbour_coord(self, coord):
+ def validate_coord(coord):
+ dim = self.dim_uv
+ if (0 <= coord[0] < dim[0]) and (0 <= coord[1] < dim[1]) and (0 <= coord[2] < dim[2]):
+ return coord
+ else:
+ return None
+
+ z, y, x = coord
+ t, d = (z-1, y, x), (z+1, y, x) # t: top, d: down
+ f, b = (z, y-1, x), (z, y+1, x) # f: front, b: back
+ l, r = (z, y, x-1), (z, y, x+1) # l: left, r: right
+ return [validate_coord(i) for i in [t, d, f, b, l, r]]
+
+ def get_neighbour_ind(self, coord):
+ neighbours = [self.coord_to_ind(i) for i in self.get_neighbour_coord(coord)]
+ return np.reshape(neighbours, (3, 2))
+
+ def get_neighbour_ind3comp(self, coord):
+ neighbours = [self.coord_to_ind3comp(i) for i in self.get_neighbour_coord(coord)]
+ return np.reshape(np.swapaxes(neighbours, 0, 1), (3, 3, 2))
+
+
+class IndexConverter2D(object):
+
+ def __init__(self, dim_uv):
+ assert len(dim_uv) == 2, 'Dimensions have to be two-dimensioal!'
+ self.dim_uv = dim_uv
+ self.size_2d = np.prod(dim_uv)
+
+ def ind_to_coord(self, ind):
+ assert ind < self.size_2d, 'Index out of range!'
+ v, u = int(ind/self.dim_uv[1]), ind % self.dim_uv[1]
+ return v, u
+
+ def ind2comp_to_coord(self, ind):
+ assert ind < 2 * self.size_2d, 'Index out of range!'
+ c, remain = int(ind/self.size_2d), ind % self.size_2d
+ v, u = int(remain/self.dim_uv[1]), remain % self.dim_uv[1]
+ return c, v, u
+
+ def coord_to_ind(self, coord):
+ if coord is None:
+ return None
+ v, u = coord
+ ind = v*self.dim_uv[1] + u
+ assert ind < self.size_2d, 'Index out of range!'
+ return ind
+
+ def coord_to_ind2comp(self, coord):
+ if coord is None:
+ return [None, None]
+ v, u = coord
+ ind = [i*self.size_2d + v*self.dim_uv[1] + u for i in range(2)]
+ assert ind[-1] < 2 * self.size_2d, 'Index out of range!'
+ return ind
+
+ def get_neighbour_coord(self, coord):
+ def validate_coord(coord):
+ if (0 <= coord[0] < self.dim_uv[0]) and (0 <= coord[1] < self.dim_uv[1]):
+ return coord
+ else:
+ return None
+ v, u = coord
+ f, b = (v-1, u), (v+1, u) # f: front, b: back
+ l, r = (v, u-1), (v, u+1) # l: left, r: right
+ return [validate_coord(i) for i in [f, b, l, r]]
+
+ def get_neighbour_ind(self, coord):
+ neighbours = [self.coord_to_ind(i) for i in self.get_neighbour_coord(coord)]
+ return np.reshape(neighbours, (2, 2))
+
+ def get_neighbour_ind2comp(self, coord):
+ neighbours = [self.coord_to_ind2comp(i) for i in self.get_neighbour_coord(coord)]
+ return np.reshape(np.swapaxes(neighbours, 0, 1), (2, 2, 2))
+
+
+# TODO: method for constructing 3D mask from 2D masks?
diff --git a/pyramid/version.py b/pyramid/version.py
index 049ee68..2a9ac58 100644
--- a/pyramid/version.py
+++ b/pyramid/version.py
@@ -1,3 +1,3 @@
# THIS FILE IS GENERATED FROM THE PYRAMID SETUP.PY
version="0.1.0-dev"
-hg_revision="aa1763114c1c+"
+hg_revision="ac57353bbe0b+"
diff --git a/scripts/collaborations/patrick_nanowire_reconstruction.py b/scripts/collaborations/patrick_nanowire_reconstruction.py
index 03f31ed..7a60887 100644
--- a/scripts/collaborations/patrick_nanowire_reconstruction.py
+++ b/scripts/collaborations/patrick_nanowire_reconstruction.py
@@ -9,7 +9,7 @@ from PIL import Image
import numpy as np
from jutil.taketime import TakeTime
-from pyramid import *
+from pyramid import * # analysis:ignore
import logging
import logging.config
@@ -21,36 +21,68 @@ LOGGING_CONF = os.path.join(os.path.dirname(os.path.realpath(pyramid.__file__)),
logging.config.fileConfig(LOGGING_CONF, disable_existing_loggers=False)
###################################################################################################
-a = 1.0 # in nm
+a = 1.455 # in nm
gain = 5
b_0 = 1
lam = 1E-6
PATH = '../../output/patrick/'
-dim_uv = (128, 128)
+PHASE = '-10_pha'
+MASK = 'mask_-10_elong_0'
+longFOV = True
+longFOV_string = np.where(longFOV, 'longFOV', 'normalFOV')
+IMAGENAME = '{}_{}_{}_'.format(MASK, PHASE, longFOV_string)
+PHASE_MAX = 0.92 # -10°: 0.92, 30°: 7.68
+PHASE_MIN = -16.85 # -10°: -16.85, 30°: -18.89
+PHASE_DELTA = PHASE_MAX - PHASE_MIN
###################################################################################################
# Read in files:
-im_mask = Image.open(PATH+'min102mask.tif')
-im_mask.thumbnail(dim_uv, Image.ANTIALIAS)
+im_mask = Image.open(PATH+MASK+'.tif')
+#im_mask.thumbnail((125, 175), Image.ANTIALIAS)
+im_phase = Image.open(PATH+PHASE+'.tif')
+#im_phase.thumbnail((125, 125), Image.ANTIALIAS)
+
mask = np.array(np.array(im_mask)/255, dtype=bool)
+dim_uv = mask.shape
+phase = (np.array(im_phase)/255.-0.5) * PHASE_DELTA
+pad = dim_uv[0] - phase.shape[0]
+phase_pad = np.zeros(dim_uv)
+phase_pad[pad:, :] = phase
+
mask = np.expand_dims(mask, axis=0)
-im_phase = Image.open(PATH+'min102.tif')
-im_phase.thumbnail(dim_uv, Image.ANTIALIAS)
-phase = np.array(im_phase)/255. - 0.5#*(18.19977+8.057761) - 8.057761
-phase_map = PhaseMap(a, phase)
dim = mask.shape
+phase_map = PhaseMap(a, phase)
+phase_map_pad = PhaseMap(a, phase_pad)
+
data_set = DataSet(a, dim, b_0, mask)
-data_set.append(phase_map, SimpleProjector(dim))
+data_set.append(phase_map_pad, SimpleProjector(dim))
+
+# Create Se_inv
+if longFOV:
+ mask_Se = np.ones(dim_uv)
+ mask_Se[:pad, :] = 0
+ data_set.set_Se_inv_diag_with_masks([mask_Se])
regularisator = FirstOrderRegularisator(mask, lam, p=2)
with TakeTime('reconstruction time'):
mag_data_rec = rc.optimize_linear(data_set, regularisator=regularisator, max_iter=50)
+phase_map_rec_pad = pm(mag_data_rec)
+phase_map_rec = PhaseMap(a, phase_map_rec_pad.phase[pad:, :])
+phase_map_diff = phase_map_rec - phase_map
+
# Display the reconstructed phase map and holography image:
-phase_map.display_combined('Input PhaseMap', gain=40)
+phase_map.display_combined('Input PhaseMap', gain=4)
+plt.savefig(PATH+IMAGENAME+'ORIGINAL.png')
+phase_map_pad.display_combined('Input PhaseMap (padded)', gain=4)
+phase_map_rec_pad.display_combined('Reconstr. Distribution (padded)', gain=4)
+plt.savefig(PATH+IMAGENAME+'RECONSTRUCTION_PADDED.png')
+phase_map_rec.display_combined('Reconstr. Distribution', gain=4)
+plt.savefig(PATH+IMAGENAME+'RECONSTRUCTION.png')
+phase_map_diff.display_combined('Difference')
+plt.savefig(PATH+IMAGENAME+'DIFFERENCE.png')
+mag_data_rec.scale_down(4)
mag_data_rec.quiver_plot(log=True)
-phase_map_rec = pm(mag_data_rec)
-phase_map_rec.display_combined('Reconstr. Distribution', gain=40)
-#plt.savefig(dirname + "/reconstr.png")
+plt.savefig(PATH+IMAGENAME+'MAGNETIZATION_DOWNSCALE4.png')
diff --git a/scripts/reconstruction/reconstruction_linear_test.py b/scripts/reconstruction/reconstruction_linear_test.py
index 2b29a26..e6a74cc 100644
--- a/scripts/reconstruction/reconstruction_linear_test.py
+++ b/scripts/reconstruction/reconstruction_linear_test.py
@@ -36,6 +36,7 @@ dim = (32, 32, 32)
dim_uv = dim[1:3]
count = 16
lam = 1E-4
+use_fftw = True
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
@@ -51,9 +52,9 @@ shape = mc.Shapes.disc(dim, center, radius_shell, height)
magnitude = mc.create_mag_dist_vortex(shape)
mag_data = MagData(a, magnitude)
-#mag_data.quiver_plot('z-projection', proj_axis='z')
-#mag_data.quiver_plot('y-projection', proj_axis='y')
-#mag_data.quiver_plot3d('Original distribution')
+mag_data.quiver_plot('z-projection', proj_axis='z')
+mag_data.quiver_plot('y-projection', proj_axis='y')
+mag_data.quiver_plot3d('Original distribution')
tilts_full = np.linspace(-pi/2, pi/2, num=count/2, endpoint=False)
tilts_miss = np.linspace(-pi/3, pi/3, num=count/2, endpoint=False)
@@ -76,7 +77,7 @@ noise = 0
print('--Setting up data collection')
mask = mag_data.get_mask()
-data = DataSet(a, dim, b_0, mask)
+data = DataSet(a, dim, b_0, mask, use_fftw=use_fftw)
data.projectors = projectors
data.phase_maps = data.create_phase_maps(mag_data)
@@ -89,12 +90,13 @@ if noise != 0:
print('--Reconstruction')
reg = FirstOrderRegularisator(mask, lam, p=2)
+info = []
with TakeTime('reconstruction'):
- mag_data_opt = rc.optimize_linear(data, regularisator=reg, max_iter=100)
+ mag_data_opt = rc.optimize_linear(data, regularisator=reg, max_iter=100, info=info)
###################################################################################################
print('--Plot stuff')
mag_data_opt.quiver_plot3d('Reconstructed distribution')
-#(mag_data_opt - mag_data).quiver_plot3d('Difference')
-#phase_maps_opt = data.create_phase_maps(mag_data_opt)
+(mag_data_opt - mag_data).quiver_plot3d('Difference')
+phase_maps_opt = data.create_phase_maps(mag_data_opt)
diff --git a/scripts/reconstruction/reconstruction_linear_test_batch.py b/scripts/reconstruction/reconstruction_linear_test_batch.py
index 38a25fc..3f310a8 100644
--- a/scripts/reconstruction/reconstruction_linear_test_batch.py
+++ b/scripts/reconstruction/reconstruction_linear_test_batch.py
@@ -11,12 +11,7 @@ import itertools
from mayavi import mlab
-import pyramid
-from pyramid.magdata import MagData
-from pyramid.phasemap import PhaseMap
-from pyramid.projector import YTiltProjector, XTiltProjector
-from pyramid.dataset import DataSet
-from pyramid.regularisator import ZeroOrderRegularisator, FirstOrderRegularisator
+from pyramid import * #analysis:ignore
import pyramid.magcreator as mc
import pyramid.reconstruction as rc
@@ -32,7 +27,7 @@ import logging.config
LOGGING_CONF = os.path.join(os.path.dirname(os.path.realpath(pyramid.__file__)), 'logging.ini')
-PATH = '../../output/3d reconstruction/'
+PATH = '../../output/3d reconstruction/Phase 1/'
logging.config.fileConfig(LOGGING_CONF, disable_existing_loggers=False)
if not os.path.exists(PATH):
@@ -48,36 +43,38 @@ b_0 = 1.
print('--Magnetization distributions')
dim = (32, 32, 32)
-center = (dim[0]/2-0.5, int(dim[1]/2)-0.5, int(dim[2]/2)-0.5)
+center = (dim[0]/2-0.5, dim[1]/2-0.5, dim[2]/2-0.5)
radius_core = dim[1]/8
radius_shell = dim[1]/4
height = dim[0]/2
+width_cube = (dim[0]/4, dim[1]/4, dim[2]/4)
+width_slab = (dim[0]/4, dim[1]/4, dim[2]/2)
# Vortex:
mag_shape_disc = mc.Shapes.disc(dim, center, radius_shell, height)
mag_data_vortex = MagData(a, mc.create_mag_dist_vortex(mag_shape_disc))
-# Sphere:
-mag_shape_sphere = mc.Shapes.sphere(dim, center, radius_shell)
-mag_data_sphere = MagData(a, mc.create_mag_dist_homog(mag_shape_sphere, phi=pi/4, theta=pi/4))
# Core-Shell:
mag_shape_core = mc.Shapes.disc(dim, center, radius_core, height)
mag_shape_shell = np.logical_xor(mag_shape_disc, mag_shape_core)
mag_data_core_shell = MagData(a, mc.create_mag_dist_vortex(mag_shape_shell, magnitude=0.75))
mag_data_core_shell += MagData(a, mc.create_mag_dist_homog(mag_shape_core, phi=0, theta=0))
-
-# TODO: cubic bar magnet, rectangular bar magnet, -sphere
+# Cubic slab:
+mag_shape_cube = mc.Shapes.slab(dim, center, width_cube)
+mag_data_cube = MagData(a, mc.create_mag_dist_homog(mag_shape_cube, 0))
+# Rectangular slab:
+mag_shape_slab = mc.Shapes.slab(dim, center, width_slab)
+mag_data_slab = MagData(a, mc.create_mag_dist_homog(mag_shape_slab, 0))
###################################################################################################
print('--Set up configurations')
-mag_datas = {'vortex_disc': mag_data_vortex, 'homog_sphere': mag_data_sphere,
- 'core_shell': mag_data_core_shell}
-masks = {'mask': True}
-xy_tilts = {'xy': [True, True]}
+mag_datas = {'rect_slab': mag_data_slab}
+masks = {'mask': True, 'nomask': False}
+xy_tilts = {'y': [False, True]}
max_tilts = [60]
tilt_steps = [10]
noises = [0]
orders = [1]
-lambdas = [1E-4]
+lambdas = [1E-2, 1E-3, 1E-4, 1E-5, 1E-6]
# Combining everything:
config_list = [mag_datas.keys(), masks.keys(), xy_tilts.keys(),
max_tilts, tilt_steps, noises, orders, lambdas]
@@ -93,13 +90,15 @@ for key in mag_datas:
###################################################################################################
print('--Reconstruction')
-print 'Number of configurations:', len(list(itertools.product(*config_list)))
+total_number = len(list(itertools.product(*config_list)))
+print 'Number of configurations:', total_number
-for configuration in itertools.product(*config_list):
+for i, configuration in enumerate(itertools.product(*config_list)):
# Extract keys:
mag_key, mask_key, xy_key, max_tilt, tilt_step, noise, order, lam = configuration
- print '----CONFIG:', configuration
- name = '{}_{}_{}tilts_max{}_in{}steps_{}noise_{}order_{}lam'.format(*configuration)
+ print '----CONFIG {} of {}:'.format(i+1, total_number)
+ print configuration
+ name = '{}_{}_{}tilts_max{}_in{}steps_{}noise_{}order_{:.0E}lam'.format(*configuration)
dim = mag_datas[mag_key].dim
# Mask:
if masks[mask_key]:
@@ -140,7 +139,7 @@ for configuration in itertools.product(*config_list):
(mag_data_rec - mag_datas[mag_key]).quiver_plot3d('Difference ({})'.format(mag_key))
mlab.savefig('{}{}_DIFF.png'.format(PATH, name), size=(800, 800))
mlab.close(all=True)
- data_shelve = shelve.open(PATH+'/3d_shelve')
+ data_shelve = shelve.open(PATH+'3d_shelve')
data_shelve[name] = info
data_shelve.close()
print 'chisq = {:.6f}, chisq_m = {:.6f}, chisq_a = {:.6f}'.format(*info)
diff --git a/scripts/test methods/compare_methods.py b/scripts/test methods/compare_methods.py
index a164ac9..dabd9f3 100644
--- a/scripts/test methods/compare_methods.py
+++ b/scripts/test methods/compare_methods.py
@@ -30,10 +30,8 @@ logging.config.fileConfig(LOGGING_CONF, disable_existing_loggers=False)
b_0 = 1.0 # in T
a = 1.0 # in nm
phi = pi/4
-numcore = False
-padding = 1
gain = 'auto'
-dim = (32, 256, 256) # in px (z, y, x)
+dim = (16, 128, 128) # in px (z, y, x)
# Create magnetic shape:
geometry = 'disc'
@@ -59,23 +57,32 @@ mag_data = MagData(a, mc.create_mag_dist_homog(mag_shape, phi))
projector = SimpleProjector(dim)
projection = projector(mag_data)
# Construct PhaseMapper objects:
-pm_real = PhaseMapperRDRC(Kernel(a, projector.dim_uv, b_0), numcore=numcore)
-pm_conv = PhaseMapperRDFC(Kernel(a, projector.dim_uv, b_0))
-pm_four = PhaseMapperFDFC(a, projector.dim_uv, b_0, padding=padding)
+pm_real_slow = PhaseMapperRDRC(Kernel(a, projector.dim_uv, b_0), numcore=False)
+pm_real_fast = PhaseMapperRDRC(Kernel(a, projector.dim_uv, b_0), numcore=True)
+pm_conv_slow = PhaseMapperRDFC(Kernel(a, projector.dim_uv, b_0, use_fftw=False))
+pm_conv_fast = PhaseMapperRDFC(Kernel(a, projector.dim_uv, b_0, use_fftw=True))
+pm_four_pad0 = PhaseMapperFDFC(a, projector.dim_uv, b_0, padding=0)
+pm_four_pad1 = PhaseMapperFDFC(a, projector.dim_uv, b_0, padding=1)
# Get times for different approaches:
start_time = time.time()
-phase_map_real = pm_real(projection)
-print 'Time for RDRC: ', time.time() - start_time
+phase_map_real = pm_real_slow(projection)
+print 'Time for RDRC no numcore :', time.time() - start_time
start_time = time.time()
-phase_map_conv = pm_conv(projection)
-print 'Time for RDFC: ', time.time() - start_time
+phase_map_real = pm_real_fast(projection)
+print 'Time for RDRC numcore :', time.time() - start_time
start_time = time.time()
-phase_map_conv = pm_conv(projection)
-print 'Time for RDFC: ', time.time() - start_time
+phase_map_conv = pm_conv_slow(projection)
+print 'Time for RDFC numpy conv.:', time.time() - start_time
start_time = time.time()
-phase_map_four = pm_four(projection)
-print 'Time for FDFC: ', time.time() - start_time
+phase_map_conv = pm_conv_fast(projection)
+print 'Time for RDFC FFTW conv.:', time.time() - start_time
+start_time = time.time()
+phase_map_four = pm_four_pad0(projection)
+print 'Time for FDFC padding 0 :', time.time() - start_time
+start_time = time.time()
+phase_map_four = pm_four_pad1(projection)
+print 'Time for FDFC padding 1 :', time.time() - start_time
# Display the combinated plots with phasemap and holography image:
phase_map_ana.display_combined('Analytic Solution', gain=gain)
--
GitLab