From ec0809400be1a23400bfedfa54e0a1271e8cd13d Mon Sep 17 00:00:00 2001
From: Jan Caron <j.caron@fz-juelich.de>
Date: Tue, 11 Mar 2014 10:01:45 +0100
Subject: [PATCH] First simple reconstruction approach documentation: minor
changes and updates of docstrings
---
pyramid/__init__.py | 2 +-
pyramid/costfunction.py | 15 +-
pyramid/datacollection.py | 31 ++--
pyramid/forwardmodel.py | 31 +++-
pyramid/kernel.py | 46 +++--
pyramid/logfile.log | 0
pyramid/optimizer.py | 15 +-
pyramid/phasemap.py | 43 ++---
pyramid/phasemapper.py | 52 +++---
pyramid/projector.py | 75 ++++++--
scripts/edinburgh_test.py | 2 +-
scripts/images_poster.py | 112 +++++++-----
scripts/logfile.log | 1 -
.../ch5-4-comparison_of_methods_new.py | 4 +-
scripts/paper 1/logfile.log | 169 +++++++++--------
scripts/simple_phasemapping.py | 2 +-
scripts/simple_reconstruction.py | 171 ++++++++++++++++++
scripts/test_projection.py | 68 +++++++
18 files changed, 588 insertions(+), 251 deletions(-)
create mode 100644 pyramid/logfile.log
create mode 100644 scripts/simple_reconstruction.py
create mode 100644 scripts/test_projection.py
diff --git a/pyramid/__init__.py b/pyramid/__init__.py
index 8d25565..2316e34 100644
--- a/pyramid/__init__.py
+++ b/pyramid/__init__.py
@@ -8,7 +8,7 @@ magcreator
magdata
Class for the storage of magnetization data.
projector
- Class for projections of given magnetization distribution.
+ Class for projecting given magnetization distributions.
kernel
Class for the kernel matrix representing one magnetized pixel.
phasemapper
diff --git a/pyramid/costfunction.py b/pyramid/costfunction.py
index 6201f91..811a2a9 100644
--- a/pyramid/costfunction.py
+++ b/pyramid/costfunction.py
@@ -26,20 +26,29 @@ class Costfunction:
y = self.y
F = self.F
Se_inv = self.Se_inv
- return (F(x)-y).dot(Se_inv.dot(F(x)-y))
+# print 'Costfunction - __call__ - input: ', len(x)
+ result = (F(x)-y).dot(Se_inv.dot(F(x)-y))
+# print 'Costfunction - __call__ - output:', result
+ return result
def jac(self, x):
# TODO: Docstring!
y = self.y
F = self.F
Se_inv = self.Se_inv
- return F.jac_T_dot(x, (Se_inv.dot(F(x)-y)))
+# print 'Costfunction - jac - input: ', len(x)
+ result = F.jac_T_dot(x, Se_inv.dot(F(x)-y))
+# print 'Costfunction - jac - output:', len(result)
+ return result
def hess_dot(self, x, vector):
# TODO: Docstring!
F = self.F
Se_inv = self.Se_inv
- return F.jac_T_dot(x, Se_inv.dot(F.jac_dot(x, vector)))
+# print 'Costfunction - hess_dot - input: ', len(vector)
+ result = F.jac_T_dot(x, Se_inv.dot(F.jac_dot(x, vector)))
+# print 'Costfunction - hess_dot - output:', len(result)
+ return result
# TODO: Murks!
diff --git a/pyramid/datacollection.py b/pyramid/datacollection.py
index 54d0371..496173e 100644
--- a/pyramid/datacollection.py
+++ b/pyramid/datacollection.py
@@ -17,8 +17,8 @@ class DataCollection(object):
'''Class for collecting phase maps and corresponding projectors.
Represents a collection of (e.g. experimentally derived) phase maps, stored as
- :class:´~.PhaseMap´ objects and corresponding projectors stored as :class:`~.Projector`
- objects. At creation, the grid spacing `a` and the dimension `dim` of the projected grid.
+ :class:`~.PhaseMap` objects and corresponding projectors stored as :class:`~.Projector`
+ objects. At creation, the grid spacing `a` and the dimension `dim_uv` of the projected grid.
Data can be added via the :func:`~.append` method, where a :class:`~.PhaseMap` and a
:class:`~.Projector` have to be given as tuple argument.
@@ -26,7 +26,7 @@ class DataCollection(object):
----------
a: float
The grid spacing in nm.
- dim: tuple (N=2)
+ dim_uv: tuple (N=2)
Dimensions of the projected grid.
phase_maps:
A list of all stored :class:`~.PhaseMap` objects.
@@ -42,8 +42,8 @@ class DataCollection(object):
return self._a
@property
- def dim(self):
- return self._dim
+ def dim_uv(self):
+ return self._dim_uv
@property
def phase_vec(self):
@@ -57,25 +57,28 @@ class DataCollection(object):
def projectors(self):
return [d[1] for d in self.data]
- def __init__(self, a, dim):
+ def __init__(self, a, dim_uv, b_0):
self.log = logging.getLogger(__name__)
- self.log.info('Creating '+str(self))
assert isinstance(a, Number), 'Grid spacing has to be a number!'
assert a >= 0, 'Grid spacing has to be a positive number!'
self._a = a
- assert isinstance(dim, tuple) and len(dim)==2, 'Dimension has to be a tuple of length 2!'
+ assert isinstance(dim_uv, tuple) and len(dim_uv)==2, \
+ 'Dimension has to be a tuple of length 2!'
assert np.all([])
- self.dim = dim
+ self._dim_uv = dim_uv
+ self.b_0 = b_0
self.data = []
+ # TODO: make it work
+# self.log.info('Created:', str(self))
def __repr__(self):
self.log.info('Calling __repr__')
- return '%s(a=%r, dim=%r)' % (self.__class__, self.a, self.dim)
+ return '%s(a=%r, dim_uv=%r)' % (self.__class__, self.a, self.dim_uv)
def __str__(self):
self.log.info('Calling __str__')
- return 'DataCollection(a=%s, dim=%s, data_count=%s)' % \
- (self.__class__, self.a, self.dim, len(self.phase_maps))
+ return 'DataCollection(%s, a=%s, dim_uv=%s, data_count=%s)' % \
+ (self.__class__, self.a, self.dim_uv, len(self.data))
def append(self, (phase_map, projector)):
'''Appends a data pair of phase map and projection infos to the data collection.`
@@ -94,6 +97,6 @@ class DataCollection(object):
self.log.info('Calling append')
assert isinstance(phase_map, PhaseMap) and isinstance(projector, Projector), \
'Argument has to be a tuple of a PhaseMap and a Projector object!'
- assert phase_map.dim == self.dim, 'Added phasemap must have the same dimension!'
- assert (projector.dim_v, projector.dim_u) == self.dim, 'Projector dimensions must match!'
+ assert phase_map.dim_uv == self.dim_uv, 'Added phasemap must have the same dimension!'
+ assert projector.dim_uv == self.dim_uv, 'Projector dimensions must match!'
self.data.append((phase_map, projector))
diff --git a/pyramid/forwardmodel.py b/pyramid/forwardmodel.py
index 89d1d7a..c231829 100644
--- a/pyramid/forwardmodel.py
+++ b/pyramid/forwardmodel.py
@@ -17,7 +17,7 @@ class ForwardModel:
@property
def projectors(self):
return self._projectors
-
+
@projectors.setter
def projectors(self, projectors):
assert np.all([isinstance(projector, Projector) for projector in projectors]), \
@@ -33,27 +33,42 @@ class ForwardModel:
assert isinstance(kernel, Kernel), 'A Kernel object has to be provided!'
self._kernel = kernel
- def __init__(self,projectors, kernel):
+ def __init__(self, projectors, kernel, b_0):
# TODO: Docstring!
self.kernel = kernel
- self.b_0 = kernel.b_0
+ self.b_0 = b_0
self.a = kernel.a
- self.dim = kernel.dim
+ self.dim_uv = kernel.dim_uv
self.projectors = projectors
def __call__(self, x):
# TODO: Docstring!
+# print 'FWD Model - __call__ - input: ', len(x)
result = [self.kernel.jac_dot(projector.jac_dot(x)) for projector in self.projectors]
- return np.reshape(result, -1)
+ result = self.b_0 * np.reshape(result, -1)
+# print 'FWD Model - __call__ - output:', len(result)
+ return result
# TODO: jac_dot ausschreiben!!!
def jac_dot(self, x, vector):
# TODO: Docstring! multiplication with the jacobi-matrix (may depend on x)
- return self(vector) # The jacobi-matrix does not depend on x in a linear problem
+# print 'FWD Model - jac_dot - input: ', len(vector)
+ result = self(vector)
+# print 'FWD Model - jac_dot - output:', len(result)
+ return result # The jacobi-matrix does not depend on x in a linear problem
def jac_T_dot(self, x, vector):
# TODO: Docstring! multiplication with the jacobi-matrix (may depend on x)
# The jacobi-matrix does not depend on x in a linear problem
- result = [projector.jac_T_dot(self.kernel.jac_T_dot(x)) for projector in self.projectors]
- return np.reshape(result, -1)
+# print 'FWD Model - jac_T_dot - input: ', len(vector)
+ size_3d = self.projectors[0].size_3d
+ size_2d = np.prod(self.dim_uv)
+ result = np.zeros(3*size_3d)
+ for (i, projector) in enumerate(self.projectors):
+ result += projector.jac_T_dot(self.kernel.jac_T_dot(vector[i*size_2d:(i+1)*size_2d]))
+# result = [projector.jac_T_dot(self.kernel.jac_T_dot(vector[i*size_2d:(i+1)*size_2d]))
+# for (i, projector) in enumerate(self.projectors)]
+ result = np.reshape(result, -1)
+# print 'FWD Model - jac_T_dot - output:', len(result)
+ return result
diff --git a/pyramid/kernel.py b/pyramid/kernel.py
index b8653e3..0ce1963 100644
--- a/pyramid/kernel.py
+++ b/pyramid/kernel.py
@@ -51,14 +51,12 @@ class Kernel(object):
Attributes
----------
- dim : tuple (N=2)
+ dim_uv : tuple (N=2)
Dimensions of the projected magnetization grid.
a : float
The grid spacing in nm.
geometry : {'disc', 'slab'}, optional
The elementary geometry of the single magnetized pixel.
- b_0 : float, optional
- The saturation magnetic induction. Default is 1.
u : :class:`~numpy.ndarray` (N=3)
The phase contribution of one pixel magnetized in u-direction.
v : :class:`~numpy.ndarray` (N=3)
@@ -69,7 +67,8 @@ class Kernel(object):
The real FFT of the phase contribution of one pixel magnetized in v-direction.
dim_fft : tuple (N=2)
Dimensions of the grid, which is used for the FFT. Calculated by adding the dimensions
- `dim` of the magnetization grid and the dimensions of the kernel (given by ``2*dim-1``)
+ `dim_uv` of the magnetization grid and the dimensions of the kernel (given by
+ ``2*dim_uv-1``)
and increasing to the next multiple of 2 (for faster FFT).
slice_fft : tuple (N=2) of :class:`slice`
A tuple of :class:`slice` objects to extract the original field of view from the increased
@@ -77,24 +76,21 @@ class Kernel(object):
'''# TODO: Can be used for several PhaseMappers via the fft arguments or via calling!
- def __init__(self, a, dim, b_0=1, numcore=True, geometry='disc'):
+ def __init__(self, a, dim_uv, numcore=True, geometry='disc'):
'''Constructor for a :class:`~.Kernel` object for representing a kernel matrix.
Parameters
----------
a : float
The grid spacing in nm.
- dim : tuple (N=2)
+ dim_uv : tuple (N=2)
Dimensions of the projected magnetization grid.
- b_0 : float, optional
- The saturation magnetic induction. Default is 1.
geometry : {'disc', 'slab'}, optional
The elementary geometry of the single magnetized pixel.
''' # TODO: Docstring
self.log = logging.getLogger(__name__)
self.log.info('Calling __init__')
- # TODO: Check if b_0 has an influence or is forgotten
# Function for the phase of an elementary geometry:
def get_elementary_phase(geometry, n, m, a):
if geometry == 'disc':
@@ -108,23 +104,22 @@ 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))
# Set basic properties:
- self.dim = dim # !!! size of the FOV, not the kernel (kernel is bigger)!
+ self.dim_uv = dim_uv # !!! size of the FOV, not the kernel (kernel is bigger)!
self.a = a
self.numcore = numcore
self.geometry = geometry
- self.b_0 = b_0
# Calculate kernel (single pixel phase):
- coeff = -a**2 * b_0 / (2*PHI_0)
- v_dim, u_dim = dim
+ coeff = -a**2 / (2*PHI_0)
+ 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 = coeff * get_elementary_phase(geometry, uu, vv, a)
self.v = coeff * get_elementary_phase(geometry, vv, uu, a)
# Calculate Fourier trafo of kernel components:
- dim_combined = 3*np.array(dim) - 1 # dim + (2*dim - 1) magnetisation + kernel
+ dim_combined = 3*np.array(dim_uv) - 1 # dim_uv + (2*dim_uv - 1) magnetisation + kernel
self.dim_fft = 2 ** np.ceil(np.log2(dim_combined)).astype(int) # next multiple of 2
- self.slice_fft = (slice(dim[0]-1, 2*dim[0]-1), slice(dim[1]-1, 2*dim[1]-1))
+ 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 = np.fft.rfftn(self.u, self.dim_fft)
self.v_fft = np.fft.rfftn(self.v, self.dim_fft)
self.log.info('Created '+str(self))
@@ -132,6 +127,7 @@ class Kernel(object):
def __call__(self, x):
'''Test'''
self.log.info('Calling __call__')
+# print 'Kernel - __call__:', len(x)
if self.numcore:
return self._multiply_jacobi_core(x)
else:
@@ -140,12 +136,13 @@ class Kernel(object):
def __repr__(self):
self.log.info('Calling __repr__')
- return '%s(a=%r, dim=%r, b_0=%r, numcore=%r, geometry=%r)' % \
- (self.__class__, self.a, self.dim, self.b_0, self.numcore, self.geometry)
+ return '%s(a=%r, dim_uv=%r, numcore=%r, geometry=%r)' % \
+ (self.__class__, self.a, self.dim_uv, self.numcore, self.geometry)
def jac_dot(self, vector):
'''TEST'''# TODO: Docstring
self.log.info('Calling jac_dot')
+# print 'Kernel - jac_dot:', len(vector)
if self.numcore:
return self._multiply_jacobi_core(vector)
else:
@@ -154,6 +151,7 @@ class Kernel(object):
def jac_T_dot(self, vector):
# TODO: Docstring
self.log.info('Calling jac_dot_T')
+# print 'Kernel - jac_T_dot:', len(vector)
return self._multiply_jacobi_T(vector)
def _multiply_jacobi(self, vector):
@@ -173,8 +171,8 @@ class Kernel(object):
'''# TODO: move!
self.log.info('Calling _multiply_jacobi')
- v_dim, u_dim = self.dim
- size = v_dim * u_dim
+ v_dim, u_dim = self.dim_uv
+ size = np.prod(self.dim_uv)
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)
@@ -206,9 +204,9 @@ class Kernel(object):
'''# TODO: move!
self.log.info('Calling _multiply_jacobi_T')
- v_dim, u_dim = self.dim
- size = v_dim * u_dim
- assert len(vector) == size, 'vector size not compatible!'
+ v_dim, u_dim = self.dim_uv
+ size = np.prod(self.dim_uv)
+ assert len(vector) == size, 'vector size not compatible! vector: {}, size: {}'.format(len(vector),size)
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
@@ -223,6 +221,6 @@ class Kernel(object):
def _multiply_jacobi_core(self, vector):
self.log.info('Calling _multiply_jacobi_core')
- result = np.zeros(np.prod(self.dim))
- core.multiply_jacobi_core(self.dim[0], self.dim[1], self.u, self.v, vector, result)
+ result = np.zeros(np.prod(self.dim_uv))
+ core.multiply_jacobi_core(self.dim_uv[0], self.dim_uv[1], self.u, self.v, vector, result)
return result
diff --git a/pyramid/logfile.log b/pyramid/logfile.log
new file mode 100644
index 0000000..e69de29
diff --git a/pyramid/optimizer.py b/pyramid/optimizer.py
index ca22601..f45546e 100644
--- a/pyramid/optimizer.py
+++ b/pyramid/optimizer.py
@@ -34,21 +34,24 @@ from pyramid.magdata import MagData
-def optimize_cg(self, data_collection, first_guess):
+def optimize_cg(data_collection, first_guess):
# TODO: where should data_collection and first_guess go? here or __init__?
- data = data_collection
+ data = data_collection # TODO: name the parameter 'data' directly...
mag_0 = first_guess
x_0 = first_guess.mag_vec
y = data.phase_vec
- kern = Kernel(data.dim, data.a, data.b_0)
- F = ForwardModel(data, kern)
+ kern = Kernel(data.a, data.dim_uv)
+ F = ForwardModel(data.projectors, kern, data.b_0)
C = Costfunction(y, F)
- result = minimize(C, x_0, method='CG', jac=C.jac, hessp=C.hess_dot)
+ result = minimize(C, x_0, method='Newton-CG', jac=C.jac, hessp=C.hess_dot,
+ options={'maxiter':200, 'disp':True})
x_opt = result.x
- mag_opt = MagData(mag_0.a, np.zeros((3,)+mag_0.dim)).mag_vec = x_opt
+ mag_opt = MagData(mag_0.a, np.zeros((3,)+mag_0.dim))
+
+ mag_opt.mag_vec = x_opt
return mag_opt
diff --git a/pyramid/phasemap.py b/pyramid/phasemap.py
index 293dd63..be2b294 100644
--- a/pyramid/phasemap.py
+++ b/pyramid/phasemap.py
@@ -38,7 +38,7 @@ class PhaseMap(object):
----------
a: float
The grid spacing in nm.
- dim: tuple (N=2)
+ dim_uv: tuple (N=2)
Dimensions of the grid.
phase: :class:`~numpy.ndarray` (N=2)
Matrix containing the phase shift.
@@ -107,8 +107,8 @@ class PhaseMap(object):
self._a = a
@property
- def dim(self):
- return self._dim
+ def dim_uv(self):
+ return self._dim_uv
@property
def phase(self):
@@ -119,7 +119,7 @@ class PhaseMap(object):
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._dim = phase.shape
+ self._dim_uv = phase.shape
@property
def phase_vec(self):
@@ -128,8 +128,8 @@ class PhaseMap(object):
@phase_vec.setter
def phase_vec(self, phase_vec):
assert isinstance(phase_vec, np.ndarray), 'Vector has to be a numpy array!'
- assert np.size(phase_vec) == np.prod(self.dim), 'Vector size has to match phase!'
- self.phase = phase_vec.reshape(self.dim)
+ assert np.size(phase_vec) == np.prod(self.dim_uv), 'Vector size has to match phase!'
+ self.phase = phase_vec.reshape(self.dim_uv)
@property
def unit(self):
@@ -154,7 +154,7 @@ class PhaseMap(object):
def __str__(self):
self.log.info('Calling __str__')
- return 'PhaseMap(a=%s, dim=%s)' % (self.a, self.dim)
+ return 'PhaseMap(a=%s, dim_uv=%s)' % (self.a, self.dim_uv)
def __neg__(self): # -self
self.log.info('Calling __neg__')
@@ -167,7 +167,7 @@ class PhaseMap(object):
if isinstance(other, PhaseMap):
self.log.info('Adding two PhaseMap objects')
assert other.a == self.a, 'Added phase has to have the same grid spacing!'
- assert other.phase.shape == self.dim, \
+ assert other.phase.shape == self.dim_uv, \
'Added magnitude has to have the same dimensions!'
return PhaseMap(self.a, self.phase+other.phase, self.unit)
else: # other is a Number
@@ -266,8 +266,8 @@ class PhaseMap(object):
self.log.info('Calling save_to_netcdf4')
phase_file = netCDF4.Dataset(filename, 'w', format='NETCDF4')
phase_file.a = self.a
- phase_file.createDimension('v_dim', self.dim[0])
- phase_file.createDimension('u_dim', self.dim[1])
+ phase_file.createDimension('v_dim', self.dim_uv[0])
+ phase_file.createDimension('u_dim', self.dim_uv[1])
phase = phase_file.createVariable('phase', 'f', ('v_dim', 'u_dim'))
phase[:] = self.phase
phase_file.close()
@@ -334,8 +334,8 @@ class PhaseMap(object):
# Plot the phasemap:
im = axis.pcolormesh(phase, cmap=cmap, vmin=-limit, vmax=limit, norm=norm)
# Set the axes ticks and labels:
- axis.set_xlim(0, self.dim[1])
- axis.set_ylim(0, self.dim[0])
+ axis.set_xlim(0, self.dim_uv[1])
+ 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)))
@@ -381,8 +381,8 @@ class PhaseMap(object):
fig = plt.figure()
axis = Axes3D(fig)
# Plot surface and contours:
- u = range(self.dim[1])
- v = range(self.dim[0])
+ u = range(self.dim_uv[1])
+ v = range(self.dim_uv[0])
uu, vv = np.meshgrid(u, v)
axis.plot_surface(uu, vv, phase, rstride=4, cstride=4, alpha=0.7, cmap=cmap,
linewidth=0, antialiased=False)
@@ -458,8 +458,8 @@ class PhaseMap(object):
axis.set_title(title, fontsize=18)
axis.set_xlabel('x-axis [px]', fontsize=15)
axis.set_ylabel('y-axis [px]', fontsize=15)
- axis.set_xlim(0, self.dim[1])
- axis.set_ylim(0, self.dim[0])
+ axis.set_xlim(0, self.dim_uv[1])
+ 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))
# Show Plot:
@@ -467,9 +467,9 @@ class PhaseMap(object):
plt.show()
# Return plotting axis:
return axis
-
- def display_combined(self, density=1, title='Combined Plot', interpolation='none',
- grad_encode='dark'):
+
+ def display_combined(self, title='Combined Plot', cmap='RdBu', limit=None, norm=None,
+ density=1, interpolation='none', grad_encode='dark'):
'''Display the phase map and the resulting color coded holography image in one plot.
Parameters
@@ -499,7 +499,7 @@ class PhaseMap(object):
# Plot phase map:
phase_axis = fig.add_subplot(1, 2, 2, aspect='equal')
fig.subplots_adjust(right=0.85)
- self.display_phase(axis=phase_axis, show=False)
+ self.display_phase(cmap='RdBu', limit=limit, norm=norm, axis=phase_axis, show=False)
plt.show()
# Return the plotting axes:
return phase_axis, holo_axis
@@ -537,6 +537,3 @@ class PhaseMap(object):
axis.xaxis.set_major_locator(NullLocator())
axis.yaxis.set_major_locator(NullLocator())
plt.show()
-
-
-PhaseMap.make_color_wheel()
diff --git a/pyramid/phasemapper.py b/pyramid/phasemapper.py
index 772e6fa..4501e37 100644
--- a/pyramid/phasemapper.py
+++ b/pyramid/phasemapper.py
@@ -46,13 +46,13 @@ class PMAdapterFM(PhaseMapper):
assert isinstance(projector, Projector), 'Argument has to be a Projector object!'
self.a = a
self.projector = projector
- self.fwd_model = ForwardModel([projector], Kernel(a, projector.dim_2d, b_0, geometry))
+ self.fwd_model = ForwardModel([projector], Kernel(a, projector.dim_uv, geometry), b_0)
def __call__(self, mag_data):
assert isinstance(mag_data, MagData), 'Only MagData objects can be mapped!'
assert mag_data.a == self.a, 'Grid spacing has to match!'
# TODO: test if mag_data fits in all aspects
- phase_map = PhaseMap(self.a, np.zeros(self.projector.dim_2d))
+ phase_map = PhaseMap(self.a, np.zeros(self.projector.dim_uv))
phase_map.phase_vec = self.fwd_model(mag_data.mag_vec)
return phase_map
@@ -94,8 +94,8 @@ class PMFourier(PhaseMapper):
assert isinstance(mag_data, MagData), 'Only MagData objects can be mapped!'
assert mag_data.a == self.a, 'Grid spacing has to match!'
# TODO: test if mag_data fits in all aspects (especially with projector)
- v_dim, u_dim = self.projector.dim_2d
- u_mag, v_mag = self.projector(mag_data.mag_vec).reshape((2,)+self.projector.dim_2d)
+ v_dim, u_dim = self.projector.dim_uv
+ u_mag, v_mag = self.projector(mag_data.mag_vec).reshape((2,)+self.projector.dim_uv)
# Create zero padded matrices:
u_pad = u_dim/2 * self.padding
v_pad = v_dim/2 * self.padding
@@ -183,7 +183,7 @@ class PMElectric(PhaseMapper):
# Calculate mask:
mask = np.sqrt(np.sum(np.array(mag_data.magnitude)**2, axis=0)) > self.threshold
# Project and calculate phase:
- projection = self.projector(mask.reshape(-1)).reshape(self.projector.dim)
+ projection = self.projector(mask.reshape(-1)).reshape(self.projector.dim_uv)
phase = v_0 * Ce * projection * self.a*1E-9
return PhaseMap(self.a, phase)
@@ -200,14 +200,15 @@ class PMConvolve(PhaseMapper):
a : float
The 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.
+ 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.
b_0 : float, optional
The magnetic induction corresponding to a magnetization `M`\ :sub:`0` in T.
The default is 1.
kernel : :class:`~pyramid.kernel.Kernel`, optional
- Specifies the kernel for the convolution with the magnetization data. If none is specified,
- one will be created with `disc` as the default geometry.
+ Specifies the kernel for the convolution with the magnetization data. If none is
+ specified, one will be created with `disc` as the default geometry.
Returns
-------
@@ -218,13 +219,14 @@ class PMConvolve(PhaseMapper):
assert isinstance(projector, Projector), 'Argument has to be a Projector object!'
self.a = a
self.projector = projector
+ self.b_0 = b_0
self.threshold = threshold
- self.kernel = Kernel(a, projector.dim_2d, b_0, geometry)
+ self.kernel = Kernel(a, projector.dim_uv, geometry)
def __call__(self, mag_data):
# Docstring!
# Process input parameters:
- u_mag, v_mag = self.projector(mag_data.mag_vec).reshape((2,)+self.projector.dim_2d)
+ u_mag, v_mag = self.projector(mag_data.mag_vec).reshape((2,)+self.projector.dim_uv)
# Fourier transform the projected magnetisation:
kernel = self.kernel
u_mag_fft = np.fft.rfftn(u_mag, kernel.dim_fft)
@@ -233,7 +235,7 @@ class PMConvolve(PhaseMapper):
u_phase = np.fft.irfftn(u_mag_fft * kernel.u_fft, kernel.dim_fft)[kernel.slice_fft].copy()
v_phase = np.fft.irfftn(v_mag_fft * kernel.v_fft, kernel.dim_fft)[kernel.slice_fft].copy()
# Return the result:
- return PhaseMap(self.a, u_phase - v_phase)
+ return PhaseMap(self.a, self.b_0*(u_phase-v_phase))
class PMReal(PhaseMapper):
@@ -267,30 +269,34 @@ class PMReal(PhaseMapper):
assert isinstance(projector, Projector), 'Argument has to be a Projector object!'
self.a = a
self.projector = projector
+ self.b_0 = b_0
self.threshold = threshold
- self.kernel = Kernel(a, projector.dim_2d, b_0, geometry)
+ self.kernel = Kernel(a, projector.dim_uv, geometry)
self.numcore = numcore
def __call__(self, mag_data):
# TODO: Docstring
# Process input parameters:
- dim = self.projector.dim_2d
+ dim_uv = self.projector.dim_uv
threshold = self.threshold
- u_mag, v_mag = self.projector(mag_data.mag_vec).reshape((2,)+dim)
+ u_mag, v_mag = self.projector(mag_data.mag_vec).reshape((2,)+dim_uv)
# Create kernel (lookup-tables for the phase of one pixel):
- u_phi = self.kernel.u
- v_phi = self.kernel.v
+ u_phi = self.b_0 * self.kernel.u
+ v_phi = self.b_0 * self.kernel.v
# Calculation of the phase:
- phase = np.zeros(dim)
+ phase = np.zeros(dim_uv)
if self.numcore:
- nc.phase_mag_real_core(dim[0], dim[1], v_phi, u_phi, v_mag, u_mag, phase, threshold)
+ nc.phase_mag_real_core(dim_uv[0], dim_uv[1], v_phi, u_phi,
+ v_mag, u_mag, phase, threshold)
else:
- for j in range(dim[0]):
- for i in range(dim[1]):
- u_phase = u_phi[dim[0]-1-j:(2*dim[0]-1)-j, dim[1]-1-i:(2*dim[1]-1)-i]
+ for j in range(dim_uv[0]):
+ for i in range(dim_uv[1]):
+ u_phase = u_phi[dim_uv[0]-1-j:(2*dim_uv[0]-1)-j,
+ dim_uv[1]-1-i:(2*dim_uv[1]-1)-i]
if abs(u_mag[j, i]) > threshold:
phase += u_mag[j, i] * u_phase
- v_phase = v_phi[dim[0]-1-j:(2*dim[0]-1)-j, dim[1]-1-i:(2*dim[1]-1)-i]
+ v_phase = v_phi[dim_uv[0]-1-j:(2*dim_uv[0]-1)-j,
+ dim_uv[1]-1-i:(2*dim_uv[1]-1)-i]
if abs(v_mag[j, i]) > threshold:
phase -= v_mag[j, i] * v_phase
# Return the phase:
diff --git a/pyramid/projector.py b/pyramid/projector.py
index 85e8641..99b3e31 100644
--- a/pyramid/projector.py
+++ b/pyramid/projector.py
@@ -30,7 +30,7 @@ class Projector(object):
Attributes
----------
- dim_2d : tuple (N=2)
+ dim_uv : tuple (N=2)
Dimensions (v, u) of the projected grid.
size_3d : int
Number of voxels of the 3-dimensional grid.
@@ -48,26 +48,19 @@ class Projector(object):
__metaclass__ = abc.ABCMeta
@abc.abstractmethod
- def __init__(self, (dim_v, dim_u), weight, coeff):
+ def __init__(self, dim_uv, weight, coeff):
self.log = logging.getLogger(__name__)
self.log.info('Calling __init__')
- self.dim_2d = (dim_v, dim_u)
+ self.dim_uv = dim_uv
self.weight = weight
self.coeff = coeff
self.size_2d, self.size_3d = weight.shape
self.log.info('Created '+str(self))
def __call__(self, vector):
- self.log.info('Calling as function')
- if len(vector) == 3*self.size_3d: # mode == 'vector'
- self.log.info('mode == vector')
- return self._vector_field_projection(vector)
- elif len(vector) == self.size_3d: # mode == 'scalar'
- self.log.info('mode == scalar')
- return self._scalar_field_projection(vector)
- else:
- raise AssertionError('Vector size has to be suited either for ' \
- 'vector- or scalar-field-projection!')
+ self.log.info('Calling as function')
+# print 'Projector - __call__:', len(vector)
+ return self.jac_dot(vector)
def _vector_field_projection(self, vector):
self.log.info('Calling _vector_field_projection')
@@ -79,19 +72,42 @@ class Projector(object):
if self.coeff[0][1] != 0: # y to u
result[:size_2d] += self.coeff[0][1] * self.weight.dot(vector[size_3d:2*size_3d])
if self.coeff[0][2] != 0: # z to u
- result[:size_2d] += self.coeff[0][2] * self.weight.dot(vector[2*size_3d])
+ result[:size_2d] += self.coeff[0][2] * self.weight.dot(vector[2*size_3d:])
if self.coeff[1][0] != 0: # x to v
result[size_2d:] += self.coeff[1][0] * self.weight.dot(vector[:size_3d])
if self.coeff[1][1] != 0: # y to v
result[size_2d:] += self.coeff[1][1] * self.weight.dot(vector[size_3d:2*size_3d])
if self.coeff[1][2] != 0: # z to v
- result[size_2d:] += self.coeff[1][2] * self.weight.dot(vector[2*size_3d])
+ result[size_2d:] += self.coeff[1][2] * self.weight.dot(vector[2*size_3d:])
+ return result
+
+ def _vector_field_projection_T(self, vector):
+ self.log.info('Calling _vector_field_projection_T')
+ size_2d, size_3d = self.size_2d, self.size_3d
+ result = np.zeros(3*size_3d)
+ # Go over all possible component projections (z, y, x) to (u, v):
+ if self.coeff[0][0] != 0: # x to u
+ result[:size_3d] += self.coeff[0][0] * self.weight.T.dot(vector[:size_2d])
+ if self.coeff[0][1] != 0: # y to u
+ result[:size_3d] += self.coeff[0][1] * self.weight.T.dot(vector[size_2d:])
+ if self.coeff[0][2] != 0: # z to u
+ result[size_3d:2*size_3d] += self.coeff[0][2] * self.weight.T.dot(vector[:size_2d])
+ if self.coeff[1][0] != 0: # x to v
+ result[size_3d:2*size_3d] += self.coeff[1][0] * self.weight.T.dot(vector[size_2d:])
+ if self.coeff[1][1] != 0: # y to v
+ result[2*size_3d:] += self.coeff[1][1] * self.weight.T.dot(vector[:size_2d])
+ if self.coeff[1][2] != 0: # z to v
+ result[2*size_3d:] += self.coeff[1][2] * self.weight.T.dot(vector[size_2d:])
return result
def _scalar_field_projection(self, vector):
self.log.info('Calling _scalar_field_projection')
return np.array(self.weight.dot(vector))
+ def _scalar_field_projection_T(self, vector):
+ self.log.info('Calling _scalar_field_projection_T')
+ return np.array(self.weight.T.dot(vector))
+
def jac_dot(self, vector):
'''Multiply a `vector` with the jacobi matrix of this :class:`~.Projector` object.
@@ -109,7 +125,30 @@ class Projector(object):
'''
self.log.info('Calling jac_dot')
- return self(vector)
+# print 'Projector - jac_dot:', len(vector)
+ if len(vector) == 3*self.size_3d: # mode == 'vector'
+ self.log.info('mode == vector')
+ return self._vector_field_projection(vector)
+ elif len(vector) == self.size_3d: # mode == 'scalar'
+ self.log.info('mode == scalar')
+ return self._scalar_field_projection(vector)
+ else:
+ raise AssertionError('Vector size has to be suited either for ' \
+ 'vector- or scalar-field-projection!')
+
+ def jac_T_dot(self, vector):
+ # TODO: Docstring!
+ self.log.info('Calling jac_T_dot')
+# print 'Projector - jac_T_dot:', len(vector)
+ if len(vector) == 2*self.size_2d: # mode == 'vector'
+ self.log.info('mode == vector')
+ return self._vector_field_projection_T(vector)
+ elif len(vector) == self.size_2d: # mode == 'scalar'
+ self.log.info('mode == scalar')
+ return self._scalar_field_projection_T(vector)
+ else:
+ raise AssertionError('Vector size has to be suited either for ' \
+ 'vector- or scalar-field-projection!')
@@ -124,7 +163,7 @@ class YTiltProjector(Projector):
#
# Attributes
# ----------
-# dim_2d : tuple (N=2)
+# dim_uv : tuple (N=2)
# Dimensions (v, u) of the projected grid.
# size_3d : int
# Number of voxels of the 3-dimensional grid.
@@ -239,7 +278,7 @@ class SimpleProjector(Projector):
#
# Attributes
# ----------
-# dim_2d : tuple (N=2)
+# dim_uv : tuple (N=2)
# Dimensions (v, u) of the projected grid.
# weight : :class:`~scipy.sparse.csr_matrix` (N=2)
# The weight matrix containing the weighting coefficients which determine the influence of
diff --git a/scripts/edinburgh_test.py b/scripts/edinburgh_test.py
index f101582..27bf0cc 100644
--- a/scripts/edinburgh_test.py
+++ b/scripts/edinburgh_test.py
@@ -42,7 +42,7 @@ phasemapper = PMAdapterFM(mag_data.a, projector)
phase_map = phasemapper(mag_data)
-phase_axis = phase_map.display_combined(density=20, interpolation='bilinear')[0]
+phase_axis = phase_map.display_combined(density=20, interpolation='bilinear', grad_encode='bright')[0]
phase_axis.xaxis.set_major_formatter(FuncFormatter(lambda x, pos: '{:3.0f}'.format(x*mag_data.a)))
phase_axis.yaxis.set_major_formatter(FuncFormatter(lambda x, pos: '{:3.0f}'.format(x*mag_data.a)))
diff --git a/scripts/images_poster.py b/scripts/images_poster.py
index d008673..2135cbe 100644
--- a/scripts/images_poster.py
+++ b/scripts/images_poster.py
@@ -6,58 +6,86 @@ Created on Wed Feb 05 19:19:19 2014
"""
import os
+
from numpy import pi
+import numpy as np
import pyramid.magcreator as mc
+import pyramid.analytic as an
from pyramid.magdata import MagData
from pyramid.phasemapper import PMConvolve
from pyramid.projector import SimpleProjector
+from pyramid.phasemap import PhaseMap
import matplotlib.pyplot as plt
+from time import clock
-
-directory = '../output/poster'
-if not os.path.exists(directory):
- os.makedirs(directory)
-# Input parameters:
-a = 10.0 # nm
-dim = (64, 128, 128)
-# Slab:f
-center = (32, 32, 32) # in px (z, y, x), index starts with 0!
-width = (322, 48, 48) # in px (z, y, x)
-mag_shape_slab = mc.Shapes.slab(dim, center, width)
-# Disc:
-center = (32, 32, 96) # in px (z, y, x), index starts with 0!
-radius = 24 # in px
-height = 24 # in px
-mag_shape_disc = mc.Shapes.disc(dim, center, radius, height)
-# Sphere:
-center = (32, 96, 64) # in px (z, y, x), index starts with 0!
-radius = 24 # in px
-mag_shape_sphere = mc.Shapes.sphere(dim, center, radius)
-# Create empty MagData object and add magnetized objects:
-mag_data = MagData(a, mc.create_mag_dist_homog(mag_shape_slab, pi/6))
-mag_data += MagData(a, mc.create_mag_dist_vortex(mag_shape_disc, (32, 96)))
-mag_data += MagData(a, mc.create_mag_dist_homog(mag_shape_sphere, -pi/4))
-# Plot the magnetic distribution, phase map and holographic contour map:
-mag_data_coarse = mag_data.copy()
-mag_data_coarse.scale_down(2)
-mag_data_coarse.quiver_plot()
-plt.savefig(os.path.join(directory, 'mag.png'))
-mag_data_coarse.quiver_plot3d()
-projector = SimpleProjector(dim)
-phase_map = PMConvolve(a, projector, b_0=0.1)(mag_data)
-phase_map.display_phase()
-plt.savefig(os.path.join(directory, 'phase.png'))
-phase_map.display_holo(density=4, interpolation='bilinear')
-plt.savefig(os.path.join(directory, 'holo.png'))
+#directory = '../output/poster'
+#if not os.path.exists(directory):
+# os.makedirs(directory)
+## Input parameters:
+#a = 10.0 # nm
+#dim = (64, 128, 128)
+## Slab:f
+#center = (32, 32, 32) # in px (z, y, x), index starts with 0!
+#width = (322, 48, 48) # in px (z, y, x)
+#mag_shape_slab = mc.Shapes.slab(dim, center, width)
+## Disc:
+#center = (32, 32, 96) # in px (z, y, x), index starts with 0!
+#radius = 24 # in px
+#height = 24 # in px
+#mag_shape_disc = mc.Shapes.disc(dim, center, radius, height)
+## Sphere:
+#center = (32, 96, 64) # in px (z, y, x), index starts with 0!
+#radius = 24 # in px
+#mag_shape_sphere = mc.Shapes.sphere(dim, center, radius)
+## Create empty MagData object and add magnetized objects:
+#mag_data = MagData(a, mc.create_mag_dist_homog(mag_shape_slab, pi/6))
+#mag_data += MagData(a, mc.create_mag_dist_vortex(mag_shape_disc, (32, 96)))
+#mag_data += MagData(a, mc.create_mag_dist_homog(mag_shape_sphere, -pi/4))
+## Plot the magnetic distribution, phase map and holographic contour map:
+#mag_data_coarse = mag_data.copy()
+#mag_data_coarse.scale_down(2)
+#mag_data_coarse.quiver_plot()
+#plt.savefig(os.path.join(directory, 'mag.png'))
+#mag_data_coarse.quiver_plot3d()
+#projector = SimpleProjector(dim)
+#phase_map = PMConvolve(a, projector, b_0=0.1)(mag_data)
+#phase_map.display_phase()
+#plt.savefig(os.path.join(directory, 'phase.png'))
+#phase_map.display_holo(density=4, interpolation='bilinear')
+#plt.savefig(os.path.join(directory, 'holo.png'))
+#
+#
+#dim = (1, 9, 9)
+#mag_shape = mc.Shapes.pixel(dim, (int(dim[0]/2), int(dim[1]/2), int(dim[2]/2)))
+#mag_data_x = MagData(a, mc.create_mag_dist_homog(mag_shape, 0))
+#mag_data_y = MagData(a, mc.create_mag_dist_homog(mag_shape, pi/2))
+#phasemapper = PMConvolve(a, SimpleProjector(dim))
+#phasemapper(mag_data_x).display_phase()
+#phasemapper(mag_data_y).display_phase()
-dim = (1, 9, 9)
-mag_shape = mc.Shapes.pixel(dim, (int(dim[0]/2), int(dim[1]/2), int(dim[2]/2)))
-mag_data_x = MagData(a, mc.create_mag_dist_homog(mag_shape, 0))
-mag_data_y = MagData(a, mc.create_mag_dist_homog(mag_shape, pi/2))
+a = 1.
+dim = (128, 128, 128)
+center = (dim[0]/2, dim[1]/2, dim[2]/2)
+radius = dim[1]/4 # in px
+mag_shape = mc.Shapes.sphere(dim, center, radius)
+mag_data = MagData(a, mc.create_mag_dist_homog(mag_shape, pi/4))
phasemapper = PMConvolve(a, SimpleProjector(dim))
-phasemapper(mag_data_x).display_phase()
-phasemapper(mag_data_y).display_phase()
\ No newline at end of file
+start = clock()
+phase_map = phasemapper(mag_data)
+print 'TIME:', clock() - start
+phase_ana = an.phase_mag_sphere(dim, a, pi/4, center, radius)
+phase_diff = (phase_ana - phase_map).phase
+print 'RMS%:', np.sqrt(np.mean((phase_diff)**2))/phase_ana.phase.max()*100
+print 'max:', phase_ana.phase.max()
+print 'RMS:', np.sqrt(np.mean((phase_diff)**2))
+PhaseMap(a, phase_diff).display_phase()
+
+mag_data.scale_down(2)
+mag_data.quiver_plot()
+mag_data.quiver_plot3d()
+phase_map.display_phase()
+
diff --git a/scripts/logfile.log b/scripts/logfile.log
index 8354e33..e69de29 100644
--- a/scripts/logfile.log
+++ b/scripts/logfile.log
@@ -1 +0,0 @@
-2014-02-08 21:43:12: INFO @ <root>: Calling copy
diff --git a/scripts/paper 1/ch5-4-comparison_of_methods_new.py b/scripts/paper 1/ch5-4-comparison_of_methods_new.py
index 735db86..e976c83 100644
--- a/scripts/paper 1/ch5-4-comparison_of_methods_new.py
+++ b/scripts/paper 1/ch5-4-comparison_of_methods_new.py
@@ -25,7 +25,7 @@ from pyramid.phasemapper import PMConvolve, PMFourier
import shelve
-force_calculation = False
+force_calculation = True
print '\nACCESS SHELVE'
@@ -133,6 +133,8 @@ else:
print '------time (disc, real disc) =', data_disc_real_d[2, i]
phase_diff_disc = (phase_ana_disc-phase_num_disc) * 1E3 # in mrad -> *1000
data_disc_real_d[1, i] = np.sqrt(np.mean(phase_diff_disc.phase**2))
+ print 'TIME:', data_disc_real_d[2, i]
+ print 'RMS%:', np.sqrt(np.mean(((phase_ana_disc-phase_num_disc).phase/phase_ana_disc.phase)**2))*100, '%'
print '----CALCULATE RMS/DURATION HOMOG. MAGN. DISC'
# Analytic solution:
diff --git a/scripts/paper 1/logfile.log b/scripts/paper 1/logfile.log
index e8f8d21..3bc36a6 100644
--- a/scripts/paper 1/logfile.log
+++ b/scripts/paper 1/logfile.log
@@ -1,85 +1,84 @@
-2014-02-08 20:46:32: INFO @ <pyramid.phasemap>: Calling make_color_wheel
-2014-02-08 20:46:43: INFO @ <pyramid.projector>: Calling __init__
-2014-02-08 20:46:43: INFO @ <pyramid.projector>: Calling __init__
-2014-02-08 20:46:43: INFO @ <pyramid.projector>: Created <pyramid.projector.SimpleProjector object at 0x0A8519F0>
-2014-02-08 20:46:43: INFO @ <pyramid.projector>: Created <pyramid.projector.SimpleProjector object at 0x0A8519F0>
-2014-02-08 20:46:43: INFO @ <pyramid.phasemap>: Calling __str__
-2014-02-08 20:46:43: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim=(128, 128))
-2014-02-08 20:46:43: INFO @ <pyramid.phasemap>: Calling __str__
-2014-02-08 20:46:43: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim=(128, 128))
-2014-02-08 20:46:43: INFO @ <pyramid.projector>: Calling as function
-2014-02-08 20:46:43: INFO @ <pyramid.projector>: mode == vector
-2014-02-08 20:46:43: INFO @ <pyramid.projector>: Calling _vector_field_projection
-2014-02-08 20:46:43: INFO @ <pyramid.phasemap>: Calling __str__
-2014-02-08 20:46:43: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim=(128, 128))
-2014-02-08 20:46:43: INFO @ <pyramid.phasemap>: Calling __sub__
-2014-02-08 20:46:43: INFO @ <pyramid.phasemap>: Calling __neg__
-2014-02-08 20:46:43: INFO @ <pyramid.phasemap>: Calling __str__
-2014-02-08 20:46:43: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim=(128, 128))
-2014-02-08 20:46:43: INFO @ <pyramid.phasemap>: Calling __add__
-2014-02-08 20:46:43: INFO @ <pyramid.phasemap>: Adding two PhaseMap objects
-2014-02-08 20:46:43: INFO @ <pyramid.phasemap>: Calling __str__
-2014-02-08 20:46:43: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim=(128, 128))
-2014-02-08 20:46:43: INFO @ <pyramid.phasemap>: Calling display_phase
-2014-02-08 20:46:43: INFO @ <pyramid.projector>: Calling as function
-2014-02-08 20:46:43: INFO @ <pyramid.projector>: mode == vector
-2014-02-08 20:46:43: INFO @ <pyramid.projector>: Calling _vector_field_projection
-2014-02-08 20:46:43: INFO @ <pyramid.phasemap>: Calling __str__
-2014-02-08 20:46:43: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim=(128, 128))
-2014-02-08 20:46:43: INFO @ <pyramid.phasemap>: Calling __sub__
-2014-02-08 20:46:43: INFO @ <pyramid.phasemap>: Calling __neg__
-2014-02-08 20:46:43: INFO @ <pyramid.phasemap>: Calling __str__
-2014-02-08 20:46:43: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim=(128, 128))
-2014-02-08 20:46:43: INFO @ <pyramid.phasemap>: Calling __add__
-2014-02-08 20:46:43: INFO @ <pyramid.phasemap>: Adding two PhaseMap objects
-2014-02-08 20:46:43: INFO @ <pyramid.phasemap>: Calling __str__
-2014-02-08 20:46:43: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim=(128, 128))
-2014-02-08 20:46:43: INFO @ <pyramid.phasemap>: Calling __isub__
-2014-02-08 20:46:43: INFO @ <pyramid.phasemap>: Calling __sub__
-2014-02-08 20:46:43: INFO @ <pyramid.phasemap>: Calling __add__
-2014-02-08 20:46:43: INFO @ <pyramid.phasemap>: Adding an offset
-2014-02-08 20:46:43: INFO @ <pyramid.phasemap>: Calling __str__
-2014-02-08 20:46:43: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim=(128, 128))
-2014-02-08 20:46:43: INFO @ <pyramid.phasemap>: Calling display_phase
-2014-02-08 20:46:44: INFO @ <pyramid.projector>: Calling as function
-2014-02-08 20:46:44: INFO @ <pyramid.projector>: mode == vector
-2014-02-08 20:46:44: INFO @ <pyramid.projector>: Calling _vector_field_projection
-2014-02-08 20:46:44: INFO @ <pyramid.phasemap>: Calling __str__
-2014-02-08 20:46:44: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim=(128, 128))
-2014-02-08 20:46:44: INFO @ <pyramid.phasemap>: Calling __sub__
-2014-02-08 20:46:44: INFO @ <pyramid.phasemap>: Calling __neg__
-2014-02-08 20:46:44: INFO @ <pyramid.phasemap>: Calling __str__
-2014-02-08 20:46:44: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim=(128, 128))
-2014-02-08 20:46:44: INFO @ <pyramid.phasemap>: Calling __add__
-2014-02-08 20:46:44: INFO @ <pyramid.phasemap>: Adding two PhaseMap objects
-2014-02-08 20:46:44: INFO @ <pyramid.phasemap>: Calling __str__
-2014-02-08 20:46:44: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim=(128, 128))
-2014-02-08 20:46:44: INFO @ <pyramid.phasemap>: Calling display_phase
-2014-02-08 20:46:45: INFO @ <pyramid.projector>: Calling as function
-2014-02-08 20:46:45: INFO @ <pyramid.projector>: mode == vector
-2014-02-08 20:46:45: INFO @ <pyramid.projector>: Calling _vector_field_projection
-2014-02-08 20:46:45: INFO @ <pyramid.phasemap>: Calling __str__
-2014-02-08 20:46:45: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim=(128, 128))
-2014-02-08 20:46:45: INFO @ <pyramid.phasemap>: Calling __sub__
-2014-02-08 20:46:45: INFO @ <pyramid.phasemap>: Calling __neg__
-2014-02-08 20:46:45: INFO @ <pyramid.phasemap>: Calling __str__
-2014-02-08 20:46:45: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim=(128, 128))
-2014-02-08 20:46:45: INFO @ <pyramid.phasemap>: Calling __add__
-2014-02-08 20:46:45: INFO @ <pyramid.phasemap>: Adding two PhaseMap objects
-2014-02-08 20:46:45: INFO @ <pyramid.phasemap>: Calling __str__
-2014-02-08 20:46:45: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim=(128, 128))
-2014-02-08 20:46:45: INFO @ <pyramid.phasemap>: Calling __isub__
-2014-02-08 20:46:45: INFO @ <pyramid.phasemap>: Calling __sub__
-2014-02-08 20:46:45: INFO @ <pyramid.phasemap>: Calling __add__
-2014-02-08 20:46:45: INFO @ <pyramid.phasemap>: Adding an offset
-2014-02-08 20:46:45: INFO @ <pyramid.phasemap>: Calling __str__
-2014-02-08 20:46:45: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim=(128, 128))
-2014-02-08 20:46:45: INFO @ <pyramid.phasemap>: Calling display_phase
-2014-02-08 20:46:47: INFO @ <pyramid.projector>: Calling __init__
-2014-02-08 20:46:47: INFO @ <pyramid.projector>: Calling __init__
-2014-02-08 20:46:47: INFO @ <pyramid.projector>: Created <pyramid.projector.SimpleProjector object at 0x0B117190>
-2014-02-08 20:46:47: INFO @ <pyramid.projector>: Created <pyramid.projector.SimpleProjector object at 0x0B117190>
-2014-02-08 20:46:47: INFO @ <pyramid.phasemap>: Calling __str__
-2014-02-08 20:46:47: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim=(128, 128))
-2014-02-08 20:46:47: INFO @ <pyramid.phasemap>: Calling __str__
-2014-02-08 20:46:47: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim=(128, 128))
+2014-02-18 09:30:30: INFO @ <pyramid.projector>: Calling __init__
+2014-02-18 09:30:30: INFO @ <pyramid.projector>: Calling __init__
+2014-02-18 09:30:30: INFO @ <pyramid.projector>: Created <pyramid.projector.SimpleProjector object at 0x0DC7ECD0>
+2014-02-18 09:30:30: INFO @ <pyramid.projector>: Created <pyramid.projector.SimpleProjector object at 0x0DC7ECD0>
+2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Calling __str__
+2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim_uv=(128, 128))
+2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Calling __str__
+2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim_uv=(128, 128))
+2014-02-18 09:30:30: INFO @ <pyramid.projector>: Calling as function
+2014-02-18 09:30:30: INFO @ <pyramid.projector>: mode == vector
+2014-02-18 09:30:30: INFO @ <pyramid.projector>: Calling _vector_field_projection
+2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Calling __str__
+2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim_uv=(128, 128))
+2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Calling __sub__
+2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Calling __neg__
+2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Calling __str__
+2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim_uv=(128, 128))
+2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Calling __add__
+2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Adding two PhaseMap objects
+2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Calling __str__
+2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim_uv=(128, 128))
+2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Calling display_phase
+2014-02-18 09:30:30: INFO @ <pyramid.projector>: Calling as function
+2014-02-18 09:30:30: INFO @ <pyramid.projector>: mode == vector
+2014-02-18 09:30:30: INFO @ <pyramid.projector>: Calling _vector_field_projection
+2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Calling __str__
+2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim_uv=(128, 128))
+2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Calling __sub__
+2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Calling __neg__
+2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Calling __str__
+2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim_uv=(128, 128))
+2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Calling __add__
+2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Adding two PhaseMap objects
+2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Calling __str__
+2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim_uv=(128, 128))
+2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Calling __isub__
+2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Calling __sub__
+2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Calling __add__
+2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Adding an offset
+2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Calling __str__
+2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim_uv=(128, 128))
+2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Calling display_phase
+2014-02-18 09:30:31: INFO @ <pyramid.projector>: Calling as function
+2014-02-18 09:30:31: INFO @ <pyramid.projector>: mode == vector
+2014-02-18 09:30:31: INFO @ <pyramid.projector>: Calling _vector_field_projection
+2014-02-18 09:30:31: INFO @ <pyramid.phasemap>: Calling __str__
+2014-02-18 09:30:31: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim_uv=(128, 128))
+2014-02-18 09:30:31: INFO @ <pyramid.phasemap>: Calling __sub__
+2014-02-18 09:30:31: INFO @ <pyramid.phasemap>: Calling __neg__
+2014-02-18 09:30:31: INFO @ <pyramid.phasemap>: Calling __str__
+2014-02-18 09:30:31: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim_uv=(128, 128))
+2014-02-18 09:30:31: INFO @ <pyramid.phasemap>: Calling __add__
+2014-02-18 09:30:31: INFO @ <pyramid.phasemap>: Adding two PhaseMap objects
+2014-02-18 09:30:31: INFO @ <pyramid.phasemap>: Calling __str__
+2014-02-18 09:30:31: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim_uv=(128, 128))
+2014-02-18 09:30:31: INFO @ <pyramid.phasemap>: Calling display_phase
+2014-02-18 09:30:31: INFO @ <pyramid.projector>: Calling as function
+2014-02-18 09:30:31: INFO @ <pyramid.projector>: mode == vector
+2014-02-18 09:30:31: INFO @ <pyramid.projector>: Calling _vector_field_projection
+2014-02-18 09:30:32: INFO @ <pyramid.phasemap>: Calling __str__
+2014-02-18 09:30:32: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim_uv=(128, 128))
+2014-02-18 09:30:32: INFO @ <pyramid.phasemap>: Calling __sub__
+2014-02-18 09:30:32: INFO @ <pyramid.phasemap>: Calling __neg__
+2014-02-18 09:30:32: INFO @ <pyramid.phasemap>: Calling __str__
+2014-02-18 09:30:32: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim_uv=(128, 128))
+2014-02-18 09:30:32: INFO @ <pyramid.phasemap>: Calling __add__
+2014-02-18 09:30:32: INFO @ <pyramid.phasemap>: Adding two PhaseMap objects
+2014-02-18 09:30:32: INFO @ <pyramid.phasemap>: Calling __str__
+2014-02-18 09:30:32: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim_uv=(128, 128))
+2014-02-18 09:30:32: INFO @ <pyramid.phasemap>: Calling __isub__
+2014-02-18 09:30:32: INFO @ <pyramid.phasemap>: Calling __sub__
+2014-02-18 09:30:32: INFO @ <pyramid.phasemap>: Calling __add__
+2014-02-18 09:30:32: INFO @ <pyramid.phasemap>: Adding an offset
+2014-02-18 09:30:32: INFO @ <pyramid.phasemap>: Calling __str__
+2014-02-18 09:30:32: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim_uv=(128, 128))
+2014-02-18 09:30:32: INFO @ <pyramid.phasemap>: Calling display_phase
+2014-02-18 09:30:33: INFO @ <pyramid.projector>: Calling __init__
+2014-02-18 09:30:33: INFO @ <pyramid.projector>: Calling __init__
+2014-02-18 09:30:33: INFO @ <pyramid.projector>: Created <pyramid.projector.SimpleProjector object at 0x0DC7EB70>
+2014-02-18 09:30:33: INFO @ <pyramid.projector>: Created <pyramid.projector.SimpleProjector object at 0x0DC7EB70>
+2014-02-18 09:30:33: INFO @ <pyramid.phasemap>: Calling __str__
+2014-02-18 09:30:33: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim_uv=(128, 128))
+2014-02-18 09:30:33: INFO @ <pyramid.phasemap>: Calling __str__
+2014-02-18 09:30:33: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim_uv=(128, 128))
diff --git a/scripts/simple_phasemapping.py b/scripts/simple_phasemapping.py
index 80c6737..c8d1715 100644
--- a/scripts/simple_phasemapping.py
+++ b/scripts/simple_phasemapping.py
@@ -15,7 +15,7 @@ from time import clock
#import cProfile
-mag_data = MagData.load_from_netcdf4('../output/vtk data/rect_500x125x3.nc')
+mag_data = MagData.load_from_netcdf4('../output/vtk data/tube_90x30x30.nc')
projector = SimpleProjector(mag_data.dim)
diff --git a/scripts/simple_reconstruction.py b/scripts/simple_reconstruction.py
new file mode 100644
index 0000000..d9a81a8
--- /dev/null
+++ b/scripts/simple_reconstruction.py
@@ -0,0 +1,171 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Fri Feb 28 14:25:59 2014
+
+@author: Jan
+"""
+
+
+import numpy as np
+from numpy import pi
+
+from pyramid.magdata import MagData
+from pyramid.projector import YTiltProjector
+from pyramid.phasemapper import PMConvolve
+from pyramid.datacollection import DataCollection
+import pyramid.optimizer as opt
+
+from pyramid.kernel import Kernel
+from pyramid.forwardmodel import ForwardModel
+from pyramid.costfunction import Costfunction
+
+
+a = 1.
+b_0 = 1000.
+dim = (3, 3, 3)
+count = 1
+
+###################################################################################################
+print('--Generating input phase_maps')
+
+
+magnitude = np.zeros((3,)+dim)
+magnitude[:, int(dim[0]/2), int(dim[1]/2), int(dim[2]/2)] = 1.
+
+mag_data = MagData(a, magnitude)
+mag_data.quiver_plot3d()
+
+tilts = np.linspace(0, 2*pi, num=count, endpoint=False)
+projectors = [YTiltProjector(mag_data.dim, tilt) for tilt in tilts]
+phasemappers = [PMConvolve(mag_data.a, projector, b_0) for projector in projectors]
+phase_maps = [pm(mag_data) for pm in phasemappers]
+
+#[phase_map.display_phase(title=u'Tilt series $(\phi = {:2.1f} \pi)$'.format(tilts[i]/pi))
+# for i, phase_map in enumerate(phase_maps)]
+phase_maps[0].display_phase()
+
+###################################################################################################
+print('--Setting up data collection')
+
+dim_uv = dim[1:3]
+
+data_collection = DataCollection(a, dim_uv, b_0)
+
+[data_collection.append((phase_maps[i], projectors[i])) for i in range(count)]
+
+###################################################################################################
+print('--Test optimizer')
+
+first_guess = MagData(a, np.zeros((3,)+dim))
+
+first_guess.magnitude[1, int(dim[0]/2), int(dim[1]/2), int(dim[2]/2)] = 1
+#first_guess.magnitude[0, int(dim[0]/2), int(dim[1]/2), int(dim[2]/2)] = -1
+
+first_guess.quiver_plot3d()
+
+phase_guess = PMConvolve(first_guess.a, projectors[0], b_0)(first_guess)
+
+phase_guess.display_phase()
+
+#mag_opt = opt.optimize_cg(data_collection, first_guess)
+#
+#mag_opt.quiver_plot3d()
+#
+#phase_opt = PMConvolve(mag_opt.a, projectors[0], b_0)(mag_opt)
+#
+#phase_opt.display_phase()
+
+
+###################################################################################################
+print('--Further testing')
+
+data = data_collection
+mag_0 = first_guess
+x_0 = first_guess.mag_vec
+y = data.phase_vec
+kern = Kernel(data.a, data.dim_uv)
+F = ForwardModel(data.projectors, kern, data.b_0)
+C = Costfunction(y, F)
+
+size_3d = np.prod(dim)
+
+cost = C(first_guess.mag_vec)
+cost_grad = C.jac(first_guess.mag_vec)
+cost_grad_del_x = cost_grad.reshape((3, 3, 3, 3))[0, ...]
+cost_grad_del_y = cost_grad.reshape((3, 3, 3, 3))[1, ...]
+cost_grad_del_z = cost_grad.reshape((3, 3, 3, 3))[2, ...]
+
+
+
+
+x_t = np.asmatrix(mag_data.mag_vec).T
+y = np.asmatrix(y).T
+K = np.array([F.jac_dot(x_0, np.eye(81)[:, i]) for i in range(81)]).T
+K = np.asmatrix(K)
+lam = 10. ** -10
+KTK = K.T * K + lam * np.asmatrix(np.eye(81))
+print lam,
+#print pylab.cond(KTK),
+x_f = KTK.I * K.T * y
+print (np.asarray(K * x_f - y) ** 2).sum(),
+print (np.asarray(K * x_t - y) ** 2).sum()
+print x_f
+x_rec = np.asarray(x_f).reshape(3,3,3,3)
+print x_rec[0, ...]
+#KTK = K.T * K
+#u,s,v = pylab.svd(KTK, full_matrices=1)
+#si = np.zeros_like(s)
+#
+#si[s>lam] = 1. / s[s>lam]
+#si=np.asmatrix(np.diag(si))
+#KTKI = np.asmatrix(v).T * si * np.asmatrix(u)
+#
+#
+#x_f = KTKI * K.T * y
+#print x_f
+#
+#
+
+
+###################################################################################################
+print('--Compliance test')
+
+
+# test F
+
+# test F.jac
+
+# test F.jac_dot
+
+#F_evaluate = F(mag_data.mag_vec)
+#F_jacobi = F.jac_dot(None, mag_data.mag_vec)
+#np.testing.assert_equal(F_evaluate, F_jacobi)
+#
+#F_reverse = F.jac_T_dot(None, F_jacobi)
+#np.testing.assert_equal(F_reverse, mag_data.mag_vec)
+
+
+
+from scipy.sparse.linalg import cg
+
+K = np.asmatrix([F.jac_dot(x_0, np.eye(81)[:, i]) for i in range(81)]).T
+lam = 10. ** -10
+KTK = K.T * K + lam * np.asmatrix(np.eye(81))
+
+A = KTK#np.array([F(np.eye(81)[:, i]) for i in range(81)])
+
+b = F.jac_T_dot(None, y)
+
+b_test = np.asarray((K.T * y).T)[0]
+
+x_f = cg(A, b_test)[0]
+
+mag_data_rec = MagData(a, np.zeros((3,)+dim))
+
+mag_data_rec.mag_vec = x_f
+
+mag_data_rec.quiver_plot3d()
+
+
+
+
diff --git a/scripts/test_projection.py b/scripts/test_projection.py
new file mode 100644
index 0000000..baaeb0a
--- /dev/null
+++ b/scripts/test_projection.py
@@ -0,0 +1,68 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Fri Jan 17 14:06:01 2014
+
+@author: Jan
+"""
+
+import numpy as np
+from numpy import pi
+
+from pyramid.magdata import MagData
+from pyramid.projector import YTiltProjector
+from pyramid.phasemapper import PMConvolve
+
+from time import clock
+
+import matplotlib.pyplot as plt
+
+#data = np.loadtxt('../output/data from Edinburgh/long_grain_remapped_0p0035.txt', delimiter=',')
+#
+#a = 1000 * (data[1, 2] - data[0, 2])
+#
+#dim = len(np.unique(data[:, 2])), len(np.unique(data[:, 1])), len(np.unique(data[:, 0]))
+#
+#mag_vec = np.concatenate([data[:, 3], data[:, 4], data[:, 5]])
+#
+#x_mag = np.reshape(data[:, 3], dim, order='F')
+#y_mag = np.reshape(data[:, 4], dim, order='F')
+#z_mag = np.reshape(data[:, 5], dim, order='F')
+#
+#magnitude = np.array((x_mag, y_mag, z_mag))
+#
+#mag_data = MagData(a, magnitude)
+#
+#mag_data.pad(30, 20, 0)
+
+
+PATH = '../output/vtk data/tube_90x30x50_sat_edge_equil.gmr'
+
+mag_data = MagData.load_from_netcdf4(PATH+'.nc')
+
+#for tilt in np.linspace(0, 2*pi, num=20, endpoint=False):
+# projector = YTiltProjector(mag_data.dim, tilt)
+# phasemapper = PMConvolve(mag_data.a, projector)
+# phase_map = phasemapper(mag_data)
+# (-phase_map).display_combined(title=u'Tilt series $(\phi = {:2.1f} \pi)$'.format(tilt/pi),
+# limit=2., density=12, interpolation='bilinear',
+# grad_encode='bright')
+# plt.savefig(PATH+'_tilt_{:3.2f}pi.png'.format(tilt/pi))
+
+
+mag_data.scale_down()
+
+mag_data.quiver_plot3d()
+
+#projectors = [YTiltProjector(mag_data.dim, i)
+# for i in np.linspace(0, 2*pi, num=20, endpoint=False)]
+#
+#start = clock()
+#phasemapper = PMConvolve[PMConvolve(mag_data.a, projector) for projector in projectors]
+#print 'Overhead :', clock()-start
+#
+#start = clock()
+#phase_maps = [pm(mag_data) for pm in phasemappers]
+#print 'Phasemapping:', clock()-start
+#
+#[phase_map.display_combined(density=12, interpolation='bilinear', grad_encode='bright')
+# for phase_map in phase_maps]
--
GitLab