diff --git a/pyramid/__init__.py b/pyramid/__init__.py
index dace4b4ba8b85aa2dbe8f298f6f21094eb535545..ad35b5abdd5e561bafdf413b996d5c1a96bbe0f1 100644
--- a/pyramid/__init__.py
+++ b/pyramid/__init__.py
@@ -42,6 +42,7 @@ import logging
from . import analytic
from . import magcreator
from . import reconstruction
+from . import fft
from .costfunction import *
from .dataset import *
from .forwardmodel import *
@@ -51,7 +52,6 @@ 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__
@@ -59,7 +59,7 @@ _log = logging.getLogger(__name__)
_log.info("Starting PYRAMID V{} HG{}".format(__version__, __hg_revision__))
del logging
-__all__ = ['__version__', '__hg_revision__', 'analytic', 'magcreator', 'reconstruction']
+__all__ = ['__version__', '__hg_revision__', 'analytic', 'magcreator', 'reconstruction', 'fft']
__all__.extend(costfunction.__all__)
__all__.extend(dataset.__all__)
__all__.extend(forwardmodel.__all__)
@@ -69,4 +69,3 @@ __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 9f367797747b8a4d7ce2a9b09de142a15fb20b84..b177e05f98831bdf9c0f669067d62c0f74558d56 100644
--- a/pyramid/costfunction.py
+++ b/pyramid/costfunction.py
@@ -104,7 +104,6 @@ class Costfunction(object):
Jacobi vector which represents the cost derivative of all voxels of the magnetization.
'''
- self._log.debug('Calling jac')
assert len(x) == self.n
return (2 * self.fwd_model.jac_T_dot(x, self.Se_inv.dot(self.fwd_model(x) - self.y))
+ self.regularisator.jac(x))
@@ -128,11 +127,11 @@ class Costfunction(object):
Product of the input `vector` with the Hessian matrix of the costfunction.
'''
-# self._log.debug('Calling hess_dot') # TODO: Profiler says this was slow...
return (2 * self.fwd_model.jac_T_dot(x, self.Se_inv.dot(self.fwd_model.jac_dot(x, vector)))
+ self.regularisator.hess_dot(x, vector))
def hess_diag(self, _):
+ # TODO: Docstring!
return np.ones(self.n)
diff --git a/pyramid/dataset.py b/pyramid/dataset.py
index 1e3ad434f60bfbb12c0de6adbedf0679f8aeaee5..9f422ca59f1e6cda07006b5dca9bd881fd82bacc 100644
--- a/pyramid/dataset.py
+++ b/pyramid/dataset.py
@@ -76,8 +76,7 @@ class DataSet(object):
@property
def phase_mappers(self):
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]
+ kernel_list = [Kernel(self.a, dim_uv, 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, Se_inv=None, use_fftw=True, threads=1):
@@ -244,3 +243,6 @@ 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/fft.py b/pyramid/fft.py
new file mode 100644
index 0000000000000000000000000000000000000000..e595437f7d88ecfbca5b397be88e1083b3e12d07
--- /dev/null
+++ b/pyramid/fft.py
@@ -0,0 +1,229 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Fri Nov 28 15:30:10 2014
+
+@author: Jan
+"""
+
+# TODO: Document!
+
+import numpy as np
+import cPickle as pickle
+import os
+
+# pyFFTW depends on this
+try:
+ from collections import Counter #analysis:ignore
+except ImportError:
+ import collections_python27
+ import collections
+ collections.Counter = collections_python27.Counter
+
+try:
+ import pyfftw
+ BACKEND = 'fftw'
+except ImportError:
+ pyfftw = None
+ BACKEND = 'numpy'
+ print("pyFFTW module not found. Using numpy implementation.")
+
+
+__all__ = ['PLANS', 'FLOAT', 'COMPLEX', 'NTHREADS', 'dump_wisdom', 'load_wisdom', 'zeros', 'empty',
+ 'configure_backend', 'fftn', 'ifftn', 'rfftn', 'irfftn', 'rfftn_adj', 'irfftn_adj']
+
+
+class FFTWCache(object):
+
+ def __init__(self):
+ self.cache = dict()
+
+ def add_fftw(self, fft_type, fftw_obj, s, axes, nthreads):
+ in_arr = fftw_obj.get_input_array()
+ key = (fft_type, in_arr.shape, in_arr.dtype, nthreads)
+ self.cache[key] = fftw_obj
+
+ def lookup_fftw(self, fft_type, in_arr, s, axes, nthreads):
+ key = (fft_type, in_arr.shape, in_arr.dtype, nthreads)
+ return self.cache.get(key, None)
+
+ def clear_cache(self):
+ self.cache = dict()
+
+
+PLANS = FFTWCache()
+FLOAT = np.float32 # One convenient place to
+COMPLEX = np.complex64 # change from 32 to 64 bit
+NTHREADS = 1
+
+
+# Numpy functions:
+def _fftn_numpy(a, s=None, axes=None):
+ return np.fft.fftn(a, s, axes)
+
+
+def _ifftn_numpy(a, s=None, axes=None):
+ return np.fft.ifftn(a, s, axes)
+
+
+def _rfftn_numpy(a, s=None, axes=None):
+ return np.fft.rfftn(a, s, axes)
+
+
+def _irfftn_numpy(a, s=None, axes=None):
+ return np.fft.irfftn(a, s, axes)
+
+
+def _rfftn_adj_numpy(a):
+ n = 2 * (a.shape[-1] - 1)
+ out_shape = a.shape[:-1] + (n,)
+ out_arr = zeros(out_shape, dtype=a.dtype)
+ out_arr[:, :n] = a
+ return _ifftn_numpy(out_arr).real * np.prod(out_shape)
+
+
+def _irfftn_adj_numpy(a):
+ n = a.shape[-1] // 2 + 1
+ out_arr = _fftn_numpy(a, axis=-1) / a.shape[-1]
+ if a.shape[-1] % 2 == 0: # even
+ out_arr[:, 1:n - 1] += np.conj(out_arr[:, :n-1:-1])
+ else: # odd
+ out_arr[:, 1:n] += np.conj(out_arr[:, :n-1:-1])
+ axes = tuple(range(len(out_arr.shape[:-1])))
+ return _fftn_numpy(out_arr[:, :n], axes=axes) / np.prod(out_arr.shape[:-1])
+
+
+# FFTW functions:
+def _fftn_fftw(a, s=None, axes=None):
+ if a.dtype not in (FLOAT, COMPLEX):
+ raise TypeError('Wrong input type!')
+ fftw = PLANS.lookup_fftw('fftn', a, s, axes, NTHREADS)
+ if fftw is None:
+ fftw = pyfftw.builders.fftn(a, s, axes)
+ PLANS.add_fftw('fftn', fftw, s, axes, NTHREADS)
+ return fftw(a).copy()
+
+
+def _ifftn_fftw(a, s=None, axes=None):
+ if a.dtype not in (FLOAT, COMPLEX):
+ raise TypeError('Wrong input type!')
+ fftw = PLANS.lookup_fftw('ifftn', a, s, axes, NTHREADS)
+ if fftw is None:
+ fftw = pyfftw.builders.ifftn(a, s, axes)
+ PLANS.add_fftw('ifftn', fftw, s, axes, NTHREADS)
+ return fftw(a).copy()
+
+
+def _rfftn_fftw(a, s=None, axes=None):
+ if a.dtype != FLOAT:
+ raise TypeError('Wrong input type!')
+ fftw = PLANS.lookup_fftw('rfftn', a, s, axes, NTHREADS)
+ if fftw is None:
+ fftw = pyfftw.builders.rfftn(a, s, axes)
+ PLANS.add_fftw('rfftn', fftw, s, axes, NTHREADS)
+ return fftw(a).copy()
+
+
+def _irfftn_fftw(a, s=None, axes=None):
+ if a.dtype != COMPLEX:
+ raise TypeError('Wrong input type!')
+ fftw = PLANS.lookup_fftw('irfftn', a, s, axes, NTHREADS)
+ if fftw is None:
+ fftw = pyfftw.builders.irfftn(a, s, axes)
+ PLANS.add_fftw('irfftn', fftw, s, axes, NTHREADS)
+ return fftw(a).copy()
+
+
+def _rfftn_adj_fftw(a):
+ n = 2 * (a.shape[-1] - 1)
+ out_shape = a.shape[:-1] + (n,)
+ out_arr = zeros(out_shape, dtype=a.dtype)
+ out_arr[:, :a.shape[-1]] = a
+ return _ifftn_fftw(out_arr).real * np.prod(out_shape)
+
+
+def _irfftn_adj_fftw(a):
+ out_arr = _fftn_fftw(a, axes=(-1,)) / a.shape[-1] # FFT of last axis
+ n = a.shape[-1] // 2 + 1
+ if a.shape[-1] % 2 == 0: # even
+ out_arr[:, 1:n-1] += np.conj(out_arr[:, :n-1:-1])
+ else: # odd
+ out_arr[:, 1:n] += np.conj(out_arr[:, :n-1:-1])
+ axes = tuple(range(len(out_arr.shape[:-1])))
+ return _fftn_fftw(out_arr[:, :n], axes=axes) / np.prod(out_arr.shape[:-1])
+
+
+# These wisdom functions do nothing if pyFFTW is not available
+def dump_wisdom(fname):
+ # TODO: Docstring!
+ if pyfftw is not None:
+ with open(fname, 'wb') as fp:
+ pickle.dump(pyfftw.export_wisdom(), fp, pickle.HIGHEST_PROTOCOL)
+
+
+def load_wisdom(fname):
+ # TODO: Docstring!
+ if pyfftw is not None:
+ if not os.path.exists(fname):
+ print("Warning: Wisdom file does not exist. First time use?")
+ else:
+ with open(fname, 'rb') as fp:
+ pyfftw.import_wisdom(pickle.load(fp))
+
+
+# Array setups:
+def zeros(shape, dtype):
+ # TODO: Docstring!
+ result = np.zeros(shape, dtype)
+ if pyfftw is not None:
+ result = pyfftw.n_byte_align(result, pyfftw.simd_alignment)
+ return result
+
+
+def empty(shape, dtype):
+ # TODO: Docstring!
+ result = np.empty(shape, dtype)
+ if pyfftw is not None:
+ result = pyfftw.n_byte_align(result, pyfftw.simd_alignment)
+ return result
+
+
+# Configure backend:
+def configure_backend(backend):
+ """
+ Change FFT backend.
+
+ Supported values are "numpy" and "fftw".
+ """
+ global fftn
+ global ifftn
+ global rfftn
+ global irfftn
+ global rfftn_adj
+ global irfftn_adj
+ global BACKEND
+ if backend == "numpy":
+ fftn = _fftn_numpy
+ ifftn = _ifftn_numpy
+ rfftn = _rfftn_numpy
+ irfftn = _irfftn_numpy
+ rfftn_adj = _rfftn_adj_numpy
+ irfftn_adj = _irfftn_adj_numpy
+ BACKEND = 'numpy'
+ elif backend == "fftw":
+ if pyfftw is not None:
+ fftn = _fftn_fftw
+ ifftn = _ifftn_fftw
+ rfftn = _rfftn_fftw
+ irfftn = _irfftn_fftw
+ rfftn_adj = _rfftn_adj_fftw
+ irfftn_adj = _irfftn_adj_fftw
+ BACKEND = 'pyfftw'
+ else:
+ print("Error: FFTW requested but not available")
+
+
+# On import:
+if pyfftw is not None:
+ configure_backend("fftw")
+else:
+ configure_backend("numpy")
diff --git a/pyramid/forwardmodel.py b/pyramid/forwardmodel.py
index 6a61a17ac1302627f84ed9d1afa79a0a89f2bd1a..eacc07993ccfc532507eb2260b104ca2b7b88283 100644
--- a/pyramid/forwardmodel.py
+++ b/pyramid/forwardmodel.py
@@ -93,7 +93,6 @@ class ForwardModel(object):
`vector`.
'''
-# self._log.debug('Calling jac_dot') # TODO: Profiler says this was slow...
self.mag_data.magnitude[:] = 0
self.mag_data.set_vector(vector, self.data_set.mask)
result = np.zeros(self.m)
@@ -123,7 +122,6 @@ class ForwardModel(object):
the input `vector`.
'''
-# self._log.debug('Calling jac_T_dot') # TODO: Profiler says this was slow...
result = np.zeros(3*np.prod(self.data_set.dim))
hp = self.hook_points
for i, projector in enumerate(self.data_set.projectors):
diff --git a/pyramid/kernel.py b/pyramid/kernel.py
index c2b1664b86fb7eeb7cd823d9e7bc26f7b9491af8..8c322c5f788bd636916c62423f2b8fbd8bcabd56 100644
--- a/pyramid/kernel.py
+++ b/pyramid/kernel.py
@@ -5,6 +5,8 @@ single magnetized pixel."""
import numpy as np
+from pyramid import fft
+
import logging
@@ -57,66 +59,39 @@ class Kernel(object):
_log = logging.getLogger(__name__+'.Kernel')
- def __init__(self, a, dim_uv, b_0=1., geometry='disc', use_fftw=True, threads=1):
+ def __init__(self, a, dim_uv, b_0=1., geometry='disc', threads=1):
self._log.debug('Calling __init__')
# Set basic properties:
self.dim_uv = dim_uv # Dimensions of the FOV
self.dim_kern = tuple(2*np.array(dim_uv)-1) # Dimensions of the kernel
-# self.size = np.prod(dim_uv) # TODO: is this even used? (Pixel count)
self.a = a
self.geometry = geometry
# Set up FFT:
- if not use_fftw:
+ if fft.BACKEND == 'pyfftw':
+ self.dim_pad = tuple(2*np.array(dim_uv)) # is at least even (not nec. power of 2)
+ elif fft.BACKEND == 'numpy':
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:
- 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!')
- self.dim_fft = (self.dim_pad[0], self.dim_pad[1]/2+1) # last axis is real
+ 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
- else: # otherwise use numpy
- rfftn = lambda x: np.fft.rfftn(x, self.dim_pad)
+ self.threads = threads # TODO: make obsolete!
# 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
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 = fft.empty(self.dim_kern, fft.FLOAT)
+ self.v = fft.empty(self.dim_kern, fft.FLOAT)
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.u_fft[...] = rfftn(self.u)
- self.v_fft[...] = rfftn(self.v)
+ self.u_fft = fft.rfftn(self.u, self.dim_pad)
+ self.v_fft = fft.rfftn(self.v, self.dim_pad)
self._log.debug('Created '+str(self))
- # TODO: make pyfftw optional (SLOW if kernel has to be build every time like in pm()!)
- # TODO: test if prior build of kernel brings speed up in test_method() or test_fftw()
- # TODO: implement fftw also in phasemapper (JUST there, here: FFT TWICE and big overhead)
- # TODO: BUT allocation of u/v/u_fft/v_fft could be beneficial (try useing with numpy.fft)
- # TODO: Set plan manually? Save computation time also for kernel?
- # TODO: Multithreading?
- # TODO: TakeTime multiple runs?
-
def __repr__(self):
self._log.debug('Calling __repr__')
return '%s(a=%r, dim_uv=%r, geometry=%r)' % \
@@ -128,7 +103,7 @@ class Kernel(object):
(self.a, self.dim_uv, self.geometry)
def _get_elementary_phase(self, geometry, n, m, a):
- # TODO: Docstring! Function for the phase of an elementary geometry:
+ self._log.debug('Calling _get_elementary_phase')
if 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
@@ -141,6 +116,8 @@ class Kernel(object):
- F_a(n-0.5, m+0.5) + F_a(n+0.5, m+0.5))
def print_info(self):
+ # TODO: Docstring!
+ self._log.debug('Calling print_info')
print 'Shape of the FOV :', self.dim_uv
print 'Shape of the Kernel:', self.dim_kern
print 'Zero-padded shape :', self.dim_pad
diff --git a/pyramid/magcreator.py b/pyramid/magcreator.py
index c4286d3dc4ce9c3f6a0ffbacf711818ea4cf33d2..c3abb46ee1e495976d1ceb1eb55bbee57b75d4e2 100644
--- a/pyramid/magcreator.py
+++ b/pyramid/magcreator.py
@@ -388,3 +388,6 @@ def create_mag_dist_vortex(mag_shape, center=None, axis='z', magnitude=1):
y_mag = -np.ones(dim) * np.cos(phi) * mag_shape * magnitude
x_mag = np.zeros(dim)
return np.array([x_mag, y_mag, z_mag])
+
+# TODO: Smooth Vortex!
+# m(x,y) = (-y/r cos(phi(r/R)), x/r sin(phi(r/R)), cos(phi(r))); phi in range[0, 1]
\ No newline at end of file
diff --git a/pyramid/magdata.py b/pyramid/magdata.py
index f24686d145516b7662345750221ccc60173bcaa3..9b64f6d8bd02ed79da57520ce30abfdd8868939b 100644
--- a/pyramid/magdata.py
+++ b/pyramid/magdata.py
@@ -84,7 +84,8 @@ class MagData(object):
def mag_vec(self, mag_vec):
assert isinstance(mag_vec, np.ndarray), 'Vector has to be a numpy array!'
assert np.size(mag_vec) == 3*np.prod(self.dim), \
- 'Vector has to match magnitude dimensions! {} {}'.format(mag_vec.shape, 3*np.prod(self.dim))
+ 'Vector has to match magnitude dimensions! {} {}'.format(mag_vec.shape,
+ 3*np.prod(self.dim))
self.magnitude = mag_vec.reshape((3,)+self.dim)
def __init__(self, a, magnitude):
@@ -494,11 +495,15 @@ class MagData(object):
if scaled:
plot_u /= np.hypot(plot_u, plot_v).max()
plot_v /= np.hypot(plot_u, plot_v).max()
+ # Setup quiver:
+ dim_uv = plot_u.shape
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])
+ xx, yy = np.meshgrid(np.arange(dim_uv[0]), np.arange(dim_uv[1]))
+ axis.quiver(xx[::ad, ::ad], yy[::ad, ::ad], plot_u[::ad, ::ad], plot_v[::ad, ::ad],
+ pivot='middle', units='xy', angles=angles[::ad, ::ad],
+ scale_units='xy', scale=1./ad, headwidth=6, headlength=7)
+ axis.set_xlim(-1, dim_uv[1])
+ axis.set_ylim(-1, dim_uv[0])
axis.set_title(title, fontsize=18)
axis.set_xlabel(u_label, fontsize=15)
axis.set_ylabel(v_label, fontsize=15)
@@ -507,7 +512,35 @@ class MagData(object):
axis.yaxis.set_major_locator(MaxNLocator(nbins=9, integer=True))
return axis
- def quiver_plot3d(self, title='Magnetization Distribution'):
+ # TODO: Switch with mayavi or combine
+ def quiver_plot3d_matplotlib(self, title='Magnetization Distribution', ar_dens=1): # TODO: Doc ar_dens
+ from mpl_toolkits.mplot3d import axes3d #analysis:ignore
+ # TODO: more arguments like in the other plots and document!!!
+ a = self.a
+ dim = self.dim
+ ad = ar_dens
+ # Create points and vector components as lists:
+ zz, yy, xx = np.mgrid[a/2:(dim[0]*a-a/2):dim[0]*1j,
+ a/2:(dim[1]*a-a/2):dim[1]*1j,
+ a/2:(dim[2]*a-a/2):dim[2]*1j]
+ xx = xx[::ad, ::ad, ::ad]
+ yy = yy[::ad, ::ad, ::ad]
+ zz = zz[::ad, ::ad, ::ad]
+ x_mag = self.magnitude[0, ::ad, ::ad, ::ad]
+ y_mag = self.magnitude[1, ::ad, ::ad, ::ad]
+ z_mag = self.magnitude[2, ::ad, ::ad, ::ad]
+ # Plot them as vectors:
+ fig = plt.figure(figsize=(8.5, 7))
+ axis = fig.add_subplot(1, 1, 1)
+ axis = fig.gca(projection='3d')
+ axis.quiver(xx, yy, zz, x_mag, y_mag, z_mag)
+ axis.set_title(title, fontsize=18)
+ # TODO: add colorbar!
+# mlab.colorbar(None, label_fmt='%.2f')
+# mlab.colorbar(None, orientation='vertical')
+ return axis
+
+ def quiver_plot3d(self, title='Magnetization Distribution', ar_dens=1): # TODO: Doc ar_dens
'''Plot the magnetization as 3D-vectors in a quiverplot.
Parameters
@@ -521,19 +554,19 @@ class MagData(object):
'''
self._log.debug('Calling quiver_plot3D')
from mayavi import mlab
-
a = self.a
dim = self.dim
+ ad = ar_dens
# Create points and vector components as lists:
zz, yy, xx = np.mgrid[a/2:(dim[0]*a-a/2):dim[0]*1j,
a/2:(dim[1]*a-a/2):dim[1]*1j,
a/2:(dim[2]*a-a/2):dim[2]*1j]
- xx = xx.reshape(-1)
- yy = yy.reshape(-1)
- zz = zz.reshape(-1)
- x_mag = np.reshape(self.magnitude[0], (-1))
- y_mag = np.reshape(self.magnitude[1], (-1))
- z_mag = np.reshape(self.magnitude[2], (-1))
+ xx = xx[::ad, ::ad, ::ad].flatten()#reshape(-1)
+ yy = yy[::ad, ::ad, ::ad].flatten()#.reshape(-1)
+ zz = zz[::ad, ::ad, ::ad].flatten()#.reshape(-1)
+ x_mag = self.magnitude[0][::ad, ::ad, ::ad].flatten()#, (-1))
+ y_mag = self.magnitude[1][::ad, ::ad, ::ad].flatten()#, (-1))
+ z_mag = self.magnitude[2][::ad, ::ad, ::ad].flatten()#, (-1))
# Plot them as vectors:
mlab.figure()
plot = mlab.quiver3d(xx, yy, zz, x_mag, y_mag, z_mag, mode='arrow')
diff --git a/pyramid/phasemapper.py b/pyramid/phasemapper.py
index 32387ddf679873e6ec5ced62df41d3864cbc3f1f..76feddabac165232ad715912e0076249384515fa 100644
--- a/pyramid/phasemapper.py
+++ b/pyramid/phasemapper.py
@@ -17,15 +17,15 @@ from pyramid.magdata import MagData
from pyramid.phasemap import PhaseMap
from pyramid.projector import SimpleProjector
from pyramid.kernel import Kernel
-
-import jutil.fft as jfft
+from pyramid import fft
import logging
__all__ = ['PhaseMapperRDFC', 'PhaseMapperRDRC', 'PhaseMapperFDFC', 'pm']
+_log = logging.getLogger(__name__)
-PHI_0 = 2067.83 # magnetic flux in T*nm²
+PHI_0 = 2067.83 # magnetic flux in T*nm²
H_BAR = 6.626E-34 # Planck constant in J*s
M_E = 9.109E-31 # electron mass in kg
Q_E = 1.602E-19 # electron charge in C
@@ -50,17 +50,14 @@ class PhaseMapper(object):
@abc.abstractmethod
def __call__(self, mag_data):
- self._log.debug('Calling __call__')
raise NotImplementedError()
@abc.abstractmethod
def jac_dot(self, vector):
- self._log.debug('Calling jac_dot')
raise NotImplementedError()
@abc.abstractmethod
def jac_T_dot(self, vector):
- self._log.debug('Calling jac_T_dot')
raise NotImplementedError()
@@ -100,29 +97,9 @@ 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
- 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._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.u_mag = fft.zeros(kernel.dim_pad, dtype=fft.FLOAT)
+ self.v_mag = fft.zeros(kernel.dim_pad, dtype=fft.FLOAT)
+ self.mag_adj = fft.zeros(kernel.dim_pad, dtype=fft.FLOAT)
self._log.debug('Created '+str(self))
def __repr__(self):
@@ -146,12 +123,12 @@ class PhaseMapperRDFC(PhaseMapper):
def _convolve(self):
# Fourier transform the projected magnetisation:
- self.u_mag_fft[...] = self._rfft2(self.u_mag)
- self.v_mag_fft[...] = self._rfft2(self.v_mag)
+ self.u_mag_fft = fft.rfftn(self.u_mag)
+ self.v_mag_fft = fft.rfftn(self.v_mag)
# 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
+ 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_phase]
+ return fft.irfftn(self.phase_fft)[self.kernel.slice_phase]
def jac_dot(self, vector):
'''Calculate the product of the Jacobi matrix with a given `vector`.
@@ -169,7 +146,6 @@ class PhaseMapperRDFC(PhaseMapper):
Product of the Jacobi matrix (which is not explicitely calculated) with the vector.
'''
-# 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.kernel.slice_mag], self.v_mag[self.kernel.slice_mag] = \
@@ -193,54 +169,16 @@ class PhaseMapperRDFC(PhaseMapper):
the vector, which has ``2*N**2`` entries like a 2D magnetic projection.
'''
-# 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)
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]
+ mag_adj_fft = fft.irfftn_adj(self.mag_adj)
+ u_phase_adj_fft = mag_adj_fft * self.kernel.u_fft
+ v_phase_adj_fft = mag_adj_fft * -self.kernel.v_fft
+ u_phase_adj = fft.rfftn_adj(u_phase_adj_fft)[self.kernel.slice_phase]
+ v_phase_adj = fft.rfftn_adj(v_phase_adj_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
- 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):
@@ -341,7 +279,6 @@ class PhaseMapperRDRC(PhaseMapper):
Product of the Jacobi matrix (which is not explicitely calculated) with the vector.
'''
- self._log.debug('Calling jac_dot')
assert len(vector) == self.n, \
'vector size not compatible! vector: {}, size: {}'.format(len(vector), self.n)
v_dim, u_dim = self.kernel.dim_uv
@@ -377,7 +314,6 @@ class PhaseMapperRDRC(PhaseMapper):
the vector, which has ``2*N**2`` entries like a 2D magnetic projection.
'''
- self._log.debug('Calling jac_T_dot')
assert len(vector) == self.m, \
'vector size not compatible! vector: {}, size: {}'.format(len(vector), self.m)
v_dim, u_dim = self.dim_uv
@@ -490,7 +426,6 @@ class PhaseMapperFDFC(PhaseMapper):
Product of the Jacobi matrix (which is not explicitely calculated) with the vector.
'''
- self._log.debug('Calling jac_dot')
assert len(vector) == self.n, \
'vector size not compatible! vector: {}, size: {}'.format(len(vector), self.n)
mag_proj = MagData(self.a, np.zeros((3, 1)+self.dim_uv))
@@ -500,7 +435,6 @@ class PhaseMapperFDFC(PhaseMapper):
return self(mag_proj).phase_vec
def jac_T_dot(self, vector):
- self._log.debug('Calling jac_T_dot')
raise NotImplementedError()
# TODO: Implement!
@@ -574,12 +508,10 @@ class PhaseMapperElectric(PhaseMapper):
return PhaseMap(mag_data.a, phase)
def jac_dot(self, vector):
- self._log.debug('Calling jac_dot')
raise NotImplementedError()
# TODO: Implement?
def jac_T_dot(self, vector):
- self._log.debug('Calling jac_T_dot')
raise NotImplementedError()
# TODO: Implement?
@@ -603,7 +535,8 @@ def pm(mag_data, axis='z', dim_uv=None, b_0=1):
mag_data : :class:`~pyramid.magdata.MagData`
The reconstructed magnetic distribution as a :class:`~.MagData` object.
- ''' # TODO: use_fftw false to notes
+ '''
+ _log.debug('Calling pm')
projector = SimpleProjector(mag_data.dim, axis=axis, dim_uv=dim_uv)
- phasemapper = PhaseMapperRDFC(Kernel(mag_data.a, projector.dim_uv, use_fftw=False))
+ phasemapper = PhaseMapperRDFC(Kernel(mag_data.a, projector.dim_uv))
return phasemapper(projector(mag_data))
diff --git a/pyramid/projector.py b/pyramid/projector.py
index f28487db26be87ddd71f4e7a5552d5df9890762f..48e69443f4e021b765419c658f57b2f95739afc2 100644
--- a/pyramid/projector.py
+++ b/pyramid/projector.py
@@ -80,14 +80,13 @@ class Projector(object):
return 'Projector(dim=%s, dim_uv=%s, coeff=%s)' % (self.dim, self.dim_uv, self.coeff)
def __call__(self, mag_data):
- self._log.debug('Calling as function')
+ self._log.debug('Calling __call__')
mag_proj = MagData(mag_data.a, np.zeros((3, 1)+self.dim_uv))
magnitude_proj = self.jac_dot(mag_data.mag_vec).reshape((2, )+self.dim_uv)
mag_proj.magnitude[0:2, 0, ...] = magnitude_proj
return mag_proj
def _vector_field_projection(self, vector):
-# self._log.debug('Calling _vector_field_projection') # TODO: Profiler says this was slow...
size_2d, size_3d = self.size_2d, self.size_3d
result = np.zeros(2*size_2d)
# Go over all possible component projections (z, y, x) to (u, v):
@@ -106,7 +105,6 @@ class Projector(object):
return result
def _vector_field_projection_T(self, vector):
-# self._log.debug('Calling _vector_field_projection_T') # TODO: Profiler says this was slow...
size_2d, size_3d = self.size_2d, self.size_3d
result = np.zeros(3*size_3d)
# Go over all possible component projections (u, v) to (z, y, x):
@@ -148,12 +146,9 @@ class Projector(object):
always`size_2d`.
'''
-# self._log.debug('Calling jac_dot') # TODO: Profiler says this was slow...
if len(vector) == 3*self.size_3d: # mode == 'vector'
-# self._log.debug('mode == vector') # TODO: Profiler says this was slow...
return self._vector_field_projection(vector)
elif len(vector) == self.size_3d: # mode == 'scalar'
-# self._log.debug('mode == scalar') # TODO: Profiler says this was slow...
return self._scalar_field_projection(vector)
else:
raise AssertionError('Vector size has to be suited either for '
@@ -175,12 +170,9 @@ class Projector(object):
of the :class:`~.Projector` object.
'''
-# self._log.debug('Calling jac_T_dot') # TODO: Profiler says this was slow...
if len(vector) == 2*self.size_2d: # mode == 'vector'
-# self._log.debug('mode == vector') # TODO: Profiler says this was slow...
return self._vector_field_projection_T(vector)
elif len(vector) == self.size_2d: # mode == 'scalar'
-# self._log.debug('mode == scalar') # TODO: Profiler says this was slow...
return self._scalar_field_projection_T(vector)
else:
raise AssertionError('Vector size has to be suited either for '
@@ -227,17 +219,14 @@ class XTiltProjector(Projector):
def __init__(self, dim, tilt, dim_uv=None):
def get_position(p, m, b, size):
- self._log.debug('Calling get_position')
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):
- self._log.debug('Calling get_impact')
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
- self._log.debug('Calling get_weight')
lo, up = delta-rho, delta+rho
# Upper boundary:
if up >= 1:
diff --git a/pyramid/util.py b/pyramid/util.py
deleted file mode 100644
index 40ce45a40a87cabd666e62fa407de7fd652e4c9e..0000000000000000000000000000000000000000
--- a/pyramid/util.py
+++ /dev/null
@@ -1,131 +0,0 @@
-# -*- 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 2a9ac58e483f54930bf7081498f4532538753452..e604916aeed68543a5cb5bef5572ca62b7d13d7a 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="ac57353bbe0b+"
+version = "0.1.0-dev"
+hg_revision = "9013fd38a1b6+"
diff --git a/scripts/collaborations/patrick_nanowire_reconstruction.py b/scripts/collaborations/patrick_nanowire_reconstruction.py
index 7a60887a30145e0b9d12efee56cde1ab774ee37d..deb5d335ab965d6709aef20468929ca8d757fb2a 100644
--- a/scripts/collaborations/patrick_nanowire_reconstruction.py
+++ b/scripts/collaborations/patrick_nanowire_reconstruction.py
@@ -26,13 +26,13 @@ gain = 5
b_0 = 1
lam = 1E-6
PATH = '../../output/patrick/'
-PHASE = '-10_pha'
-MASK = 'mask_-10_elong_0'
+PHASE = '30_pha'
+MASK = 'mask30_elong_5'
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_MAX = 7.68 # -10°: 0.92, 30°: 7.68
+PHASE_MIN = -18.89 # -10°: -16.85, 30°: -18.89
PHASE_DELTA = PHASE_MAX - PHASE_MIN
###################################################################################################
@@ -45,9 +45,9 @@ im_phase = Image.open(PATH+PHASE+'.tif')
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]
+pad = dim_uv[0] - phase.shape[0]#25
phase_pad = np.zeros(dim_uv)
-phase_pad[pad:, :] = phase
+phase_pad[pad:, :] = phase#[pad:-pad, pad:-pad] = phase
mask = np.expand_dims(mask, axis=0)
dim = mask.shape
@@ -60,8 +60,8 @@ data_set.append(phase_map_pad, SimpleProjector(dim))
# Create Se_inv
if longFOV:
- mask_Se = np.ones(dim_uv)
- mask_Se[:pad, :] = 0
+ mask_Se = np.ones(dim_uv)#np.zeros(dim_uv)
+ mask_Se[:pad, :] = 0#[pad:-pad, pad:-pad] = 1
data_set.set_Se_inv_diag_with_masks([mask_Se])
regularisator = FirstOrderRegularisator(mask, lam, p=2)
@@ -70,7 +70,7 @@ 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_rec = PhaseMap(a, phase_map_rec_pad.phase[pad:, :])#[pad:-pad, pad:-pad])
phase_map_diff = phase_map_rec - phase_map
# Display the reconstructed phase map and holography image:
diff --git a/scripts/collaborations/sudoku.py b/scripts/collaborations/sudoku.py
new file mode 100644
index 0000000000000000000000000000000000000000..03ef1bb5fc1c8d8dffe9e927793c3d383247c042
--- /dev/null
+++ b/scripts/collaborations/sudoku.py
@@ -0,0 +1,74 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Fri Dec 05 09:24:02 2014
+
+@author: Jan
+"""
+
+
+import numpy as np
+
+
+class Cell(object):
+
+ def __init__(self):
+ self.reset()
+
+ def __str__(self):
+ return 'value: {}, possibilities: {}'.format(self.value, self.possibilities)
+
+ def reset(self):
+ self.value = 0
+ self.possibilities = range(1, 10)
+ np.random.shuffle(self.possibilities)
+
+
+class Field(object):
+
+ def __init__(self):
+ self.field = np.asarray([[Cell() for i in range(9)] for j in range(9)])
+ self.fill()
+
+ def __call__(self):
+ return np.asarray([[self.field[j, i].value for i in range(9)] for j in range(9)])
+
+ def get_i(self, pos):
+ return pos % 9
+
+ def get_j(self, pos):
+ return pos // 9
+
+ def get_cell(self, pos):
+ return self.field[self.get_j(pos), self.get_i(pos)]
+
+ def check_row(self, pos):
+ cell_val = self.get_cell(pos).value
+ row = self.field[self.get_j(pos), :]
+ row_vals = [cell.value for cell in row]
+ row_vals[self.get_i(pos)] = 0
+ print 'i:{}, j:{}, row_vals:{}'.format(self.get_i(pos), self.get_j(pos), row_vals)
+ print cell_val in row_vals
+ if cell_val in row_vals:
+ return False
+ else:
+ return True
+
+ def check_square(self, pos):
+ return True
+
+ def fill(self):
+ pos = 0
+ while pos < 81:
+ cell = self.get_cell(pos)
+ cell.value = cell.possibilities.pop()
+ print self()
+ print cell
+ if self.check_row(pos) and self.check_square(pos):
+ pos += 1
+ else:
+ self.get_cell(pos).reset()
+ pos -= 1
+ print 'Position: {:2}, Value: {}'.format(pos, cell.value)
+
+test = Field()
+print test()
diff --git a/scripts/reconstruction/reconstruction_linear_test.py b/scripts/reconstruction/reconstruction_linear_test.py
index e6a74cc028a7b34f8225a28c11bac7e69570be52..b26c728f64682a27f7a722fcb37d69c7a70a216b 100644
--- a/scripts/reconstruction/reconstruction_linear_test.py
+++ b/scripts/reconstruction/reconstruction_linear_test.py
@@ -52,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)
@@ -98,5 +98,7 @@ with TakeTime('reconstruction'):
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)
+
+# TODO: iterations in jutil is one to many!
diff --git a/scripts/test methods/compare_methods.py b/scripts/test methods/compare_methods.py
index dabd9f31ded1f73c62f8213ce150fd622e55c9d8..9c33812e35ac1ec9ef220c09b4316bbdb25874a7 100644
--- a/scripts/test methods/compare_methods.py
+++ b/scripts/test methods/compare_methods.py
@@ -16,6 +16,7 @@ from pyramid.projector import SimpleProjector
from pyramid.phasemapper import PhaseMapperRDRC, PhaseMapperRDFC, PhaseMapperFDFC
from pyramid.kernel import Kernel
from pyramid.magdata import MagData
+from pyramid import fft
import logging
import logging.config
@@ -31,7 +32,7 @@ b_0 = 1.0 # in T
a = 1.0 # in nm
phi = pi/4
gain = 'auto'
-dim = (16, 128, 128) # in px (z, y, x)
+dim = (128, 128, 128) # in px (z, y, x)
# Create magnetic shape:
geometry = 'disc'
@@ -56,30 +57,41 @@ elif geometry == 'sphere':
mag_data = MagData(a, mc.create_mag_dist_homog(mag_shape, phi))
projector = SimpleProjector(dim)
projection = projector(mag_data)
-# Construct PhaseMapper objects:
-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:
+# RDRC no numcore:
+pm_real_slow = PhaseMapperRDRC(Kernel(a, projector.dim_uv, b_0), numcore=False)
start_time = time.time()
phase_map_real = pm_real_slow(projection)
print 'Time for RDRC no numcore :', time.time() - start_time
+
+# RDRC with numcore:
+pm_real_fast = PhaseMapperRDRC(Kernel(a, projector.dim_uv, b_0), numcore=True)
start_time = time.time()
phase_map_real = pm_real_fast(projection)
print 'Time for RDRC numcore :', time.time() - start_time
+
+# RDFC numpy convolution:
+fft.configure_backend('numpy')
+pm_conv_slow = PhaseMapperRDFC(Kernel(a, projector.dim_uv, b_0))
start_time = time.time()
phase_map_conv = pm_conv_slow(projection)
print 'Time for RDFC numpy conv.:', time.time() - start_time
+
+# RDFC FFTW convolution:
+fft.configure_backend('fftw')
+pm_conv_fast = PhaseMapperRDFC(Kernel(a, projector.dim_uv, b_0))
start_time = time.time()
phase_map_conv = pm_conv_fast(projection)
print 'Time for RDFC FFTW conv.:', time.time() - start_time
+
+# FDFC padding 0:
+pm_four_pad0 = PhaseMapperFDFC(a, projector.dim_uv, b_0, padding=0)
start_time = time.time()
phase_map_four = pm_four_pad0(projection)
print 'Time for FDFC padding 0 :', time.time() - start_time
+
+# FDFC padding 1:
+pm_four_pad1 = PhaseMapperFDFC(a, projector.dim_uv, b_0, padding=1)
start_time = time.time()
phase_map_four = pm_four_pad1(projection)
print 'Time for FDFC padding 1 :', time.time() - start_time
diff --git a/setup.py b/setup.py
index 962fbe537f78ebf59adb189e2cd73b12040c45dd..8de8d664e173d308010164a588869959e9f04dce 100644
--- a/setup.py
+++ b/setup.py
@@ -69,8 +69,8 @@ def hg_version():
def write_version_py(filename='pyramid/version.py'):
version_string = "# THIS FILE IS GENERATED FROM THE PYRAMID SETUP.PY\n" + \
- 'version="{}"\n'.format(VERSION) + \
- 'hg_revision="{}"\n'.format(hg_version())
+ 'version = "{}"\n'.format(VERSION) + \
+ 'hg_revision = "{}"\n'.format(hg_version())
with open(os.path.join(os.path.dirname(__file__), filename), 'w') as vfile:
vfile.write(version_string)
@@ -107,7 +107,7 @@ setup(name=DISTNAME,
ext_package='pyramid/numcore',
ext_modules=[
Extension('phasemapper_core', ['pyramid/numcore/phasemapper_core.pyx'],
- include_dirs=[numpy.get_include(), numpy.get_numarray_include()],
+ include_dirs=[numpy.get_include()],
extra_compile_args=['-march=native', '-mtune=native']
)
]