From ce62f9499a668afeebee04015a25cfad7b2480b5 Mon Sep 17 00:00:00 2001
From: Jan Caron <j.caron@fz-juelich.de>
Date: Thu, 15 May 2014 17:55:37 +0200
Subject: [PATCH] Further optimization pep8: optimizations phasemapper_core:
numcore module for phasemapper scripts: compatibility with new structure
reconstruction: new scripts
---
pyramid/analytic.py | 21 +-
pyramid/costfunction.py | 4 +-
pyramid/dataset.py | 4 +-
pyramid/forwardmodel.py | 4 +-
pyramid/kernel.py | 9 +-
pyramid/magcreator.py | 84 ++--
pyramid/magdata.py | 14 +-
pyramid/numcore/__init__.py | 11 +-
pyramid/numcore/phasemapper_core.pyx | 2 +-
pyramid/phasemap.py | 58 +--
pyramid/phasemapper.py | 6 +-
pyramid/projector.py | 15 +-
pyramid/reconstruction.py | 4 +-
scripts/colaborations/robert_filter.py | 58 ---
scripts/colaborations/vtk_interpolation.py | 73 ---
scripts/colaborations/vtk_to_numpy.py | 74 ---
.../edinburgh_test.py | 55 ++-
.../vadim_integration.py | 92 ++--
.../vtk_conversion.py | 65 ++-
.../create_core_shell_sphere.py | 7 +-
.../create_multiple_samples.py | 2 +-
scripts/create distributions/create_vortex.py | 2 +-
scripts/publications/abstract_prague.py | 17 +-
.../ch5-0-evaluation_and_comparison.py | 25 +-
.../paper 1/ch5-1-evaluation_real_space.py | 40 +-
.../paper 1/ch5-2-evaluation_fourier_space.py | 34 +-
.../paper 1/ch5-3-comparison_of_methods.py | 19 +-
.../ch5-4-comparison_of_methods_new.py | 29 +-
scripts/publications/poster_pof.py | 63 +++
scripts/publications/poster_pov.py | 91 ----
...econstruction_sparse_cg_singular_values.py | 74 +++
.../reconstruction_sparse_cg_test.py | 420 ++----------------
scripts/test methods/compare_pixel_fields.py | 6 +-
scripts/test methods/kernel_comparison.py | 23 -
scripts/test methods/test_arb_projection.py | 19 +-
scripts/test methods/test_new_structure.py | 31 --
tests/test_compliance.py | 10 +-
tests/test_costfunction.py | 7 +-
tests/test_dataset.py | 7 +-
tests/test_forwardmodel.py | 7 +-
tests/test_kernel.py | 7 +-
tests/test_reconstruction.py | 140 +++++-
42 files changed, 665 insertions(+), 1068 deletions(-)
delete mode 100644 scripts/colaborations/robert_filter.py
delete mode 100644 scripts/colaborations/vtk_interpolation.py
delete mode 100644 scripts/colaborations/vtk_to_numpy.py
rename scripts/{colaborations => collaborations}/edinburgh_test.py (58%)
rename scripts/{colaborations => collaborations}/vadim_integration.py (70%)
rename scripts/{colaborations => collaborations}/vtk_conversion.py (84%)
create mode 100644 scripts/publications/poster_pof.py
delete mode 100644 scripts/publications/poster_pov.py
create mode 100644 scripts/reconstruction/reconstruction_sparse_cg_singular_values.py
delete mode 100644 scripts/test methods/kernel_comparison.py
delete mode 100644 scripts/test methods/test_new_structure.py
diff --git a/pyramid/analytic.py b/pyramid/analytic.py
index 5216fcd..a222a22 100644
--- a/pyramid/analytic.py
+++ b/pyramid/analytic.py
@@ -46,20 +46,21 @@ def phase_mag_slab(dim, a, phi, center, width, b_0=1):
'''
LOG.debug('Calling phase_mag_slab')
+
# Function for the phase:
- def phi_mag(x, y):
+ def phi_mag(x, y):
def F_0(x, y):
A = np.log(x**2 + y**2 + 1E-30)
B = np.arctan(x / (y+1E-30))
return x*A - 2*x + 2*y*B
return coeff * Lz * (- np.cos(phi) * (F_0(x-x0-Lx/2, y-y0-Ly/2)
- - F_0(x-x0+Lx/2, y-y0-Ly/2)
- - F_0(x-x0-Lx/2, y-y0+Ly/2)
- + F_0(x-x0+Lx/2, y-y0+Ly/2))
+ - F_0(x-x0+Lx/2, y-y0-Ly/2)
+ - F_0(x-x0-Lx/2, y-y0+Ly/2)
+ + F_0(x-x0+Lx/2, y-y0+Ly/2))
+ np.sin(phi) * (F_0(y-y0-Ly/2, x-x0-Lx/2)
- - F_0(y-y0+Ly/2, x-x0-Lx/2)
- - F_0(y-y0-Ly/2, x-x0+Lx/2)
- + F_0(y-y0+Ly/2, x-x0+Lx/2)))
+ - F_0(y-y0+Ly/2, x-x0-Lx/2)
+ - F_0(y-y0-Ly/2, x-x0+Lx/2)
+ + F_0(y-y0+Ly/2, x-x0+Lx/2)))
# Process input parameters:
z_dim, y_dim, x_dim = dim
y0 = a * (center[1] + 0.5) # y0, x0 define the center of a pixel,
@@ -102,6 +103,7 @@ def phase_mag_disc(dim, a, phi, center, radius, height, b_0=1):
'''
LOG.debug('Calling phase_mag_disc')
+
# Function for the phase:
def phi_mag(x, y):
r = np.hypot(x - x0, y - y0)
@@ -109,6 +111,7 @@ def phase_mag_disc(dim, a, phi, center, radius, height, b_0=1):
result = coeff * Lz * ((y - y0) * np.cos(phi) - (x - x0) * np.sin(phi))
result *= np.where(r <= R, 1, (R / r) ** 2)
return result
+
# Process input parameters:
z_dim, y_dim, x_dim = dim
y0 = a * (center[1] + 0.5) # y0, x0 have to be in the center of a pixel,
@@ -150,6 +153,7 @@ def phase_mag_sphere(dim, a, phi, center, radius, b_0=1):
'''
LOG.debug('Calling phase_mag_sphere')
+
# Function for the phase:
def phi_mag(x, y):
r = np.hypot(x - x0, y - y0)
@@ -157,6 +161,7 @@ def phase_mag_sphere(dim, a, phi, center, radius, b_0=1):
result = coeff * R ** 3 / r ** 2 * ((y - y0) * np.cos(phi) - (x - x0) * np.sin(phi))
result *= np.where(r > R, 1, (1 - (1 - (r / R) ** 2) ** (3. / 2.)))
return result
+
# Process input parameters:
z_dim, y_dim, x_dim = dim
y0 = a * (center[1] + 0.5) # y0, x0 have to be in the center of a pixel,
@@ -197,11 +202,13 @@ def phase_mag_vortex(dim, a, center, radius, height, b_0=1):
'''
LOG.debug('Calling phase_mag_vortex')
+
# Function for the phase:
def phi_mag(x, y):
r = np.hypot(x - x0, y - y0)
result = coeff * np.where(r <= R, r - R, 0)
return result
+
# Process input parameters:
z_dim, y_dim, x_dim = dim
y0 = a * (center[1] + 0.5) # y0, x0 have to be in the center of a pixel,
diff --git a/pyramid/costfunction.py b/pyramid/costfunction.py
index 4ffbe27..1d5212b 100644
--- a/pyramid/costfunction.py
+++ b/pyramid/costfunction.py
@@ -47,7 +47,7 @@ class Costfunction(object):
def __call__(self, x):
self.LOG.debug('Calling __call__')
y = self.y
- F = self.F
+ F = self.fwd_model
Se_inv = self.Se_inv
return (F(x)-y).dot(Se_inv.dot(F(x)-y))
@@ -57,7 +57,7 @@ class Costfunction(object):
def __str__(self):
self.LOG.debug('Calling __str__')
- return 'Costfunction(fwd_model=%s, lam=%s)' % (self.__class__, self.fwd_model, self.lam)
+ return 'Costfunction(fwd_model=%s, lam=%s)' % (self.fwd_model, self.lam)
def jac(self, x):
'''Calculate the derivative of the costfunction for a given magnetization distribution.
diff --git a/pyramid/dataset.py b/pyramid/dataset.py
index 4cecfe8..08a5cc3 100644
--- a/pyramid/dataset.py
+++ b/pyramid/dataset.py
@@ -67,12 +67,12 @@ class DataSet(object):
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_uv, tuple) and len(dim_uv)==2, \
+ assert isinstance(dim_uv, tuple) and len(dim_uv) == 2, \
'Dimension has to be a tuple of length 2!'
self._dim_uv = dim_uv
self.b_0 = b_0
self.data = []
- self.LOG.debug('Created:', str(self))
+ self.LOG.debug('Created: '+str(self))
def __repr__(self):
self.LOG.debug('Calling __repr__')
diff --git a/pyramid/forwardmodel.py b/pyramid/forwardmodel.py
index 6db94cc..d820f72 100644
--- a/pyramid/forwardmodel.py
+++ b/pyramid/forwardmodel.py
@@ -80,11 +80,9 @@ class ForwardModel(object):
'''
self.LOG.debug('Calling jac_dot')
- result = [self.kernel.jac_dot(projector.jac_dot(x)) for projector in self.projectors]
+ result = [self.kernel.jac_dot(projector.jac_dot(vector)) for projector in self.projectors]
result = np.reshape(result, -1)
return result
- result = self(vector)
- return result
def jac_T_dot(self, x, vector):
''''Calculate the product of the transposed Jacobi matrix with a given `vector`.
diff --git a/pyramid/kernel.py b/pyramid/kernel.py
index c1593d3..9e5e42b 100644
--- a/pyramid/kernel.py
+++ b/pyramid/kernel.py
@@ -62,6 +62,7 @@ class Kernel(object):
def __init__(self, a, dim_uv, b_0=1., numcore=True, geometry='disc'):
self.LOG.debug('Calling __init__')
+
# Function for the phase of an elementary geometry:
def get_elementary_phase(geometry, n, m, a):
if geometry == 'disc':
@@ -73,7 +74,8 @@ class Kernel(object):
B = np.arctan(n / m)
return n*A - 2*n + 2*m*B
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))
+ - F_a(n-0.5, m+0.5) + F_a(n+0.5, m+0.5))
+
# Set basic properties:
self.dim_uv = dim_uv # Size of the FOV, not the kernel (kernel is bigger)!
self.a = a
@@ -86,7 +88,7 @@ class Kernel(object):
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)
+ self.v = coeff * get_elementary_phase(geometry, vv, uu, a) # TODO: v = np.swapaxes(u, 0,1)
# Calculate Fourier trafo of kernel components:
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
@@ -100,7 +102,6 @@ class Kernel(object):
self._multiply_jacobi_method = self._multiply_jacobi
self.LOG.debug('Created '+str(self))
-
def __repr__(self):
self.LOG.debug('Calling __repr__')
return '%s(a=%r, dim_uv=%r, numcore=%r, geometry=%r)' % \
@@ -177,7 +178,7 @@ class Kernel(object):
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)
+ '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
diff --git a/pyramid/magcreator.py b/pyramid/magcreator.py
index 5608b72..299bb77 100644
--- a/pyramid/magcreator.py
+++ b/pyramid/magcreator.py
@@ -62,11 +62,11 @@ class Shapes(object):
assert np.shape(center) == (3,), 'Parameter center has to be a tuple of length 3!'
assert np.shape(width) == (3,), 'Parameter width has to be a tuple of length 3!'
mag_shape = np.array([[[abs(x - center[2]) <= width[2] / 2
- and abs(y - center[1]) <= width[1] / 2
- and abs(z - center[0]) <= width[0] / 2
- for x in range(dim[2])]
- for y in range(dim[1])]
- for z in range(dim[0])])
+ and abs(y - center[1]) <= width[1] / 2
+ and abs(z - center[0]) <= width[0] / 2
+ for x in range(dim[2])]
+ for y in range(dim[1])]
+ for z in range(dim[0])])
return mag_shape
@classmethod
@@ -100,22 +100,22 @@ class Shapes(object):
assert axis in {'z', 'y', 'x'}, 'Axis has to be x, y or z (as a string)!'
if axis == 'z':
mag_shape = np.array([[[np.hypot(x - center[2], y - center[1]) <= radius
- and abs(z - center[0]) <= height / 2
- for x in range(dim[2])]
- for y in range(dim[1])]
- for z in range(dim[0])])
+ and abs(z - center[0]) <= height / 2
+ for x in range(dim[2])]
+ for y in range(dim[1])]
+ for z in range(dim[0])])
elif axis == 'y':
mag_shape = np.array([[[np.hypot(x - center[2], z - center[0]) <= radius
- and abs(y - center[1]) <= height / 2
- for x in range(dim[2])]
- for y in range(dim[1])]
- for z in range(dim[0])])
+ and abs(y - center[1]) <= height / 2
+ for x in range(dim[2])]
+ for y in range(dim[1])]
+ for z in range(dim[0])])
elif axis == 'x':
mag_shape = np.array([[[np.hypot(y - center[1], z - center[0]) <= radius
- and abs(x - center[2]) <= height / 2
- for x in range(dim[2])]
- for y in range(dim[1])]
- for z in range(dim[0])])
+ and abs(x - center[2]) <= height / 2
+ for x in range(dim[2])]
+ for y in range(dim[1])]
+ for z in range(dim[0])])
return mag_shape
@classmethod
@@ -149,25 +149,25 @@ class Shapes(object):
assert axis in {'z', 'y', 'x'}, 'Axis has to be x, y or z (as a string)!'
if axis == 'z':
mag_shape = np.array([[[np.hypot((x-center[2])/(width[1]/2.),
- (y-center[1])/(width[0]/2.))<=1
- and abs(z - center[0]) <= height / 2
- for x in range(dim[2])]
- for y in range(dim[1])]
- for z in range(dim[0])])
+ (y-center[1])/(width[0]/2.)) <= 1
+ and abs(z - center[0]) <= height / 2
+ for x in range(dim[2])]
+ for y in range(dim[1])]
+ for z in range(dim[0])])
elif axis == 'y':
mag_shape = np.array([[[np.hypot((x-center[2])/(width[1]/2.),
- (z-center[0])/(width[0]/2.))<=1
- and abs(y-center[1]) <= height / 2
- for x in range(dim[2])]
- for y in range(dim[1])]
- for z in range(dim[0])])
+ (z-center[0])/(width[0]/2.)) <= 1
+ and abs(y-center[1]) <= height / 2
+ for x in range(dim[2])]
+ for y in range(dim[1])]
+ for z in range(dim[0])])
elif axis == 'x':
mag_shape = np.array([[[np.hypot((y-center[1])/(width[1]/2.),
- (z-center[0])/(width[0]/2.))<=1
- and abs(z - center[0]) <= height / 2
- for x in range(dim[2])]
- for y in range(dim[1])]
- for z in range(dim[0])])
+ (z-center[0])/(width[0]/2.)) <= 1
+ and abs(z - center[0]) <= height / 2
+ for x in range(dim[2])]
+ for y in range(dim[1])]
+ for z in range(dim[0])])
return mag_shape
@classmethod
@@ -194,11 +194,11 @@ class Shapes(object):
assert np.shape(center) == (3,), 'Parameter center has to be a a tuple of length 3!'
assert radius > 0 and np.shape(radius) == (), 'Radius has to be a positive scalar value!'
mag_shape = np.array([[[np.sqrt((x-center[2])**2
- + (y-center[1])**2
- + (z-center[0])**2) <= radius
- for x in range(dim[2])]
- for y in range(dim[1])]
- for z in range(dim[0])])
+ + (y-center[1])**2
+ + (z-center[0])**2) <= radius
+ for x in range(dim[2])]
+ for y in range(dim[1])]
+ for z in range(dim[0])])
return mag_shape
@classmethod
@@ -225,11 +225,11 @@ class Shapes(object):
assert np.shape(center) == (3,), 'Parameter center has to be a a tuple of length 3!'
assert np.shape(width) == (3,), 'Parameter width has to be a a tuple of length 3!'
mag_shape = np.array([[[np.sqrt((x-center[2])**2/(width[2]/2)**2
- + (y-center[1])**2/(width[1]/2)**2
- + (z-center[0])**2/(width[0]/2)**2) <= 1
- for x in range(dim[2])]
- for y in range(dim[1])]
- for z in range(dim[0])])
+ + (y-center[1])**2/(width[1]/2)**2
+ + (z-center[0])**2/(width[0]/2)**2) <= 1
+ for x in range(dim[2])]
+ for y in range(dim[1])]
+ for z in range(dim[0])])
return mag_shape
@classmethod
diff --git a/pyramid/magdata.py b/pyramid/magdata.py
index 00a45bd..b7845ed 100644
--- a/pyramid/magdata.py
+++ b/pyramid/magdata.py
@@ -188,7 +188,6 @@ class MagData(object):
self.magnitude = self.magnitude.reshape(
3, self.dim[0]/2, 2, self.dim[1]/2, 2, self.dim[2]/2, 2).mean(axis=(6, 4, 2))
-
def scale_up(self, n=1, order=0):
'''Scale up the magnetic distribution using spline interpolation of the requested order.
@@ -397,11 +396,11 @@ class MagData(object):
cls.LOG.debug('Calling copy')
mag_file = netCDF4.Dataset(filename, 'r', format='NETCDF4')
a = mag_file.a
- magnitude = mag_file.variables['magnitude'][...]
+ magnitude = mag_file.variables['magnitude'][...]
mag_file.close()
return MagData(a, magnitude)
- def quiver_plot(self, title='Magnetic Distribution', filename=None, axis=None,
+ def quiver_plot(self, title='Magnetization Distribution', filename=None, axis=None,
proj_axis='z', ax_slice=None):
'''Plot a slice of the magnetization as a quiver plot.
@@ -462,7 +461,7 @@ class MagData(object):
fig = plt.figure(figsize=(8.5, 7))
axis = fig.add_subplot(1, 1, 1, aspect='equal')
axis.quiver(mag_slice_u, mag_slice_v, pivot='middle', angles='xy', scale_units='xy',
- scale=1, headwidth=6, headlength=7)
+ scale=1, headwidth=6, headlength=7)
axis.set_xlim(-1, np.shape(mag_slice_u)[1])
axis.set_ylim(-1, np.shape(mag_slice_u)[0])
axis.set_title(title, fontsize=18)
@@ -474,7 +473,7 @@ class MagData(object):
plt.show()
return axis
- def quiver_plot3d(self):
+ def quiver_plot3d(self, title='Magnetization Distribution'):
'''Plot the magnetization as 3D-vectors in a quiverplot.
Parameters
@@ -504,5 +503,6 @@ class MagData(object):
plot = mlab.quiver3d(xx, yy, zz, x_mag, y_mag, z_mag, mode='arrow')
mlab.outline(plot)
mlab.axes(plot)
- mlab.title('TEST', height=0.95, size=0.25)
- mlab.colorbar(plot, orientation='vertical')
+ mlab.title(title, height=0.95, size=0.35)
+ mlab.colorbar(None, label_fmt='%.2f')
+ mlab.colorbar(None, orientation='vertical')
diff --git a/pyramid/numcore/__init__.py b/pyramid/numcore/__init__.py
index c462822..f4e72d8 100644
--- a/pyramid/numcore/__init__.py
+++ b/pyramid/numcore/__init__.py
@@ -3,12 +3,13 @@
Modules
-------
-phase_mag_real
- Provides numerical core routines for the real space approach in
- :func:`~pyramid.phasemapper.phase_mag_real` of module :mod:`~pyramid.phasemapper`.
+phasemapper_core
+ Provides numerical core routines for :class:`~.PhaseMapper` class.
+kernel_core
+ Provides numerical core routines for the :class:`~.Kernel` class.
Notes
-----
-Packages are written in `pyx`-format for the use with :mod:`~cython`.
+Packages are written in `pyx`-format for the use with :mod:`~.cython`.
-"""
\ No newline at end of file
+"""
diff --git a/pyramid/numcore/phasemapper_core.pyx b/pyramid/numcore/phasemapper_core.pyx
index 7c03eef..6b67626 100644
--- a/pyramid/numcore/phasemapper_core.pyx
+++ b/pyramid/numcore/phasemapper_core.pyx
@@ -57,4 +57,4 @@ def phase_mag_real_core(
if abs(v_m) > threshold:
for q in range(v_dim):
for p in range(u_dim):
- phase[q, p] -= v_m * v_phi[q_c + q, p_c + p]
\ No newline at end of file
+ phase[q, p] -= v_m * v_phi[q_c + q, p_c + p]
diff --git a/pyramid/phasemap.py b/pyramid/phasemap.py
index 026db0b..85b9f03 100644
--- a/pyramid/phasemap.py
+++ b/pyramid/phasemap.py
@@ -57,29 +57,29 @@ class PhaseMap(object):
u'mrad': 1E3,
u'µrad': 1E6}
- CDICT = {'red': [(0.00, 1.0, 0.0),
- (0.25, 1.0, 1.0),
- (0.50, 1.0, 1.0),
- (0.75, 0.0, 0.0),
- (1.00, 0.0, 1.0)],
-
- 'green': [(0.00, 0.0, 0.0),
- (0.25, 0.0, 0.0),
- (0.50, 1.0, 1.0),
- (0.75, 1.0, 1.0),
- (1.00, 0.0, 1.0)],
-
- 'blue': [(0.00, 1.0, 1.0),
- (0.25, 0.0, 0.0),
- (0.50, 0.0, 0.0),
- (0.75, 0.0, 0.0),
- (1.00, 1.0, 1.0)]}
-
- CDICT_INV = {'red': [(0.00, 0.0, 1.0),
- (0.25, 0.0, 0.0),
- (0.50, 0.0, 0.0),
- (0.75, 1.0, 1.0),
- (1.00, 1.0, 0.0)],
+ CDICT = {'red': [(0.00, 1.0, 0.0),
+ (0.25, 1.0, 1.0),
+ (0.50, 1.0, 1.0),
+ (0.75, 0.0, 0.0),
+ (1.00, 0.0, 1.0)],
+
+ 'green': [(0.00, 0.0, 0.0),
+ (0.25, 0.0, 0.0),
+ (0.50, 1.0, 1.0),
+ (0.75, 1.0, 1.0),
+ (1.00, 0.0, 1.0)],
+
+ 'blue': [(0.00, 1.0, 1.0),
+ (0.25, 0.0, 0.0),
+ (0.50, 0.0, 0.0),
+ (0.75, 0.0, 0.0),
+ (1.00, 1.0, 1.0)]}
+
+ CDICT_INV = {'red': [(0.00, 0.0, 1.0),
+ (0.25, 0.0, 0.0),
+ (0.50, 0.0, 0.0),
+ (0.75, 1.0, 1.0),
+ (1.00, 1.0, 0.0)],
'green': [(0.00, 1.0, 1.0),
(0.25, 1.0, 1.0),
@@ -87,11 +87,11 @@ class PhaseMap(object):
(0.75, 0.0, 0.0),
(1.00, 1.0, 0.0)],
- 'blue': [(0.00, 0.0, 0.0),
- (0.25, 1.0, 1.0),
- (0.50, 1.0, 1.0),
- (0.75, 1.0, 1.0),
- (1.00, 0.0, 0.0)]}
+ 'blue': [(0.00, 0.0, 0.0),
+ (0.25, 1.0, 1.0),
+ (0.50, 1.0, 1.0),
+ (0.75, 1.0, 1.0),
+ (1.00, 0.0, 0.0)]}
HOLO_CMAP = mpl.colors.LinearSegmentedColormap('my_colormap', CDICT, 256)
HOLO_CMAP_INV = mpl.colors.LinearSegmentedColormap('my_colormap', CDICT_INV, 256)
@@ -385,7 +385,7 @@ class PhaseMap(object):
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)
+ linewidth=0, antialiased=False)
axis.contourf(uu, vv, phase, 15, zdir='z', offset=np.min(phase), cmap=cmap)
axis.view_init(45, -135)
axis.set_xlabel('x-axis [px]')
diff --git a/pyramid/phasemapper.py b/pyramid/phasemapper.py
index 1895e81..20ec3d0 100644
--- a/pyramid/phasemapper.py
+++ b/pyramid/phasemapper.py
@@ -165,8 +165,8 @@ class PMFourier(PhaseMapper):
f_u = np.linspace(0, f_nyq, u_mag_fft.shape[1])
f_v = np.linspace(-f_nyq, f_nyq, u_mag_fft.shape[0], endpoint=False)
f_uu, f_vv = np.meshgrid(f_u, f_v)
- coeff = (1j*self.b_0) / (2*self.PHI_0)
- phase_fft = coeff*self.a * (u_mag_fft*f_vv - v_mag_fft*f_uu) / (f_uu**2 + f_vv**2 + 1e-30)
+ coeff = (1j*self.b_0*self.a) / (2*self.PHI_0)
+ phase_fft = coeff * (u_mag_fft*f_vv - v_mag_fft*f_uu) / (f_uu**2 + f_vv**2 + 1e-30)
# Transform to real space and revert padding:
phase_pad = np.fft.irfft2(np.fft.ifftshift(phase_fft, axes=0))
phase = phase_pad[v_pad:v_pad+v_dim, u_pad:u_pad+u_dim]
@@ -182,6 +182,7 @@ class PMFourier(PhaseMapper):
return 'PMFourier(a=%s, projector=%s, b_0=%s, padding=%s)' % \
(self.a, self.projector, self.b_0, self.padding)
+
class PMElectric(PhaseMapper):
'''Class representing a phase mapping strategy for the electrostatic contribution.
@@ -251,6 +252,7 @@ class PMElectric(PhaseMapper):
return 'PMElectric(a=%s, projector=%s, v_0=%s, v_acc=%s, threshold=%s)' % \
(self.a, self.projector, self.v_0, self.v_acc, self.threshold)
+
class PMConvolve(PhaseMapper):
'''Class representing a phase mapping strategy using real space discretization and Fourier
diff --git a/pyramid/projector.py b/pyramid/projector.py
index 3c13200..6dca851 100644
--- a/pyramid/projector.py
+++ b/pyramid/projector.py
@@ -88,7 +88,6 @@ class Projector(object):
'''
raise NotImplementedError()
-
def _vector_field_projection(self, vector):
self.LOG.debug('Calling _vector_field_projection')
size_2d, size_3d = self.size_2d, self.size_3d
@@ -159,7 +158,7 @@ class Projector(object):
self.LOG.debug('mode == scalar')
return self._scalar_field_projection(vector)
else:
- raise AssertionError('Vector size has to be suited either for ' \
+ raise AssertionError('Vector size has to be suited either for '
'vector- or scalar-field-projection!')
def jac_T_dot(self, vector):
@@ -186,7 +185,7 @@ class Projector(object):
self.LOG.debug('mode == scalar')
return self._scalar_field_projection_T(vector)
else:
- raise AssertionError('Vector size has to be suited either for ' \
+ raise AssertionError('Vector size has to be suited either for '
'vector- or scalar-field-projection!')
@@ -247,7 +246,7 @@ class XTiltProjector(Projector):
voxels = list(itertools.product(range(dim_proj), range(dim_perp))) # z-y-plane
# Calculate positions along the projected pixel coordinate system:
center = (dim_proj/2., dim_perp/2.)
- m = np.where(tilt<=pi, -1/np.tan(tilt+1E-30), 1/np.tan(tilt+1E-30))
+ m = np.where(tilt <= pi, -1/np.tan(tilt+1E-30), 1/np.tan(tilt+1E-30))
b = center[0] - m * center[1]
positions = get_position(voxels, m, b, dim_perp)
# Calculate weight-matrix:
@@ -273,7 +272,7 @@ class XTiltProjector(Projector):
rows = np.hstack((np.array(rows), np.array(row)+i))
# Calculate weight matrix and coefficients for jacobi matrix:
weight = csr_matrix(coo_matrix((np.tile(data, dim_rot), (rows, columns)),
- shape = (size_2d, size_3d)))
+ shape=(size_2d, size_3d)))
dim_v, dim_u = dim_perp, dim_rot
coeff = [[1, 0, 0], [0, np.cos(tilt), np.sin(tilt)]]
super(XTiltProjector, self).__init__(dim, (dim_v, dim_u), weight, coeff)
@@ -349,7 +348,7 @@ class YTiltProjector(Projector):
voxels = list(itertools.product(range(dim_proj), range(dim_perp))) # z-x-plane
# Calculate positions along the projected pixel coordinate system:
center = (dim_proj/2., dim_perp/2.)
- m = np.where(tilt<=pi, -1/np.tan(tilt+1E-30), 1/np.tan(tilt+1E-30))
+ m = np.where(tilt <= pi, -1/np.tan(tilt+1E-30), 1/np.tan(tilt+1E-30))
b = center[0] - m * center[1]
positions = get_position(voxels, m, b, dim_perp)
# Calculate weight-matrix:
@@ -375,7 +374,7 @@ class YTiltProjector(Projector):
rows = np.hstack((np.array(rows), np.array(row)+i*dim_perp))
# Calculate weight matrix and coefficients for jacobi matrix:
weight = csr_matrix(coo_matrix((np.tile(data, dim_rot), (rows, columns)),
- shape = (size_2d, size_3d)))
+ shape=(size_2d, size_3d)))
dim_v, dim_u = dim_rot, dim_perp
coeff = [[np.cos(tilt), 0, np.sin(tilt)], [0, 1, 0]]
super(YTiltProjector, self).__init__(dim, (dim_v, dim_u), weight, coeff)
@@ -443,7 +442,7 @@ class SimpleProjector(Projector):
coeff = [[0, 1, 0], [0, 0, 1]]
indices = np.array([np.arange(dim_proj) + row*dim_proj
for row in range(size_2d)]).reshape(-1)
- weight = csr_matrix((data, indices, indptr), shape = (size_2d, size_3d))
+ weight = csr_matrix((data, indices, indptr), shape=(size_2d, size_3d))
super(SimpleProjector, self).__init__(dim, (dim_v, dim_u), weight, coeff)
self.LOG.debug('Created '+str(self))
diff --git a/pyramid/reconstruction.py b/pyramid/reconstruction.py
index ac9ec75..04bdc20 100644
--- a/pyramid/reconstruction.py
+++ b/pyramid/reconstruction.py
@@ -39,7 +39,7 @@ class PrintIterator(object):
distribution. This should decrease per iteration if the algorithm converges and is only
printed for a `verbosity` of 2.
verbosity : {2, 1, 0}, optional
- Parameter defining the verposity of the output. `2` is the default and will show the
+ Parameter defining the verbosity of the output. `2` is the default and will show the
current number of the iteration and the cost of the current distribution. `2` will just
show the iteration number and `0` will prevent output all together.
@@ -150,7 +150,7 @@ def optimize_cg(data, first_guess):
cost = Costfunction(y, fwd_model)
# Optimize:
result = minimize(cost, x_0, method='Newton-CG', jac=cost.jac, hessp=cost.hess_dot,
- options={'maxiter':200, 'disp':True})
+ options={'maxiter': 200, 'disp': True})
# Create optimized MagData object:
x_opt = result.x
mag_opt = MagData(mag_0.a, np.zeros((3,)+mag_0.dim))
diff --git a/scripts/colaborations/robert_filter.py b/scripts/colaborations/robert_filter.py
deleted file mode 100644
index 8c0423d..0000000
--- a/scripts/colaborations/robert_filter.py
+++ /dev/null
@@ -1,58 +0,0 @@
-# -*- coding: utf-8 -*-
-"""
-Created on Tue Oct 29 15:05:39 2013
-
-@author: Jan
-"""
-
-import matplotlib.pyplot as plt
-import numpy as np
-from numpy import pi
-
-w = 8
-l = 108
-x = np.linspace(0, 20*l, 20*l)
-
-
-# Kastenfunktion:
-y = np.where(x <= l, 1, 0)
-
-fig = plt.figure()
-axis = fig.add_subplot(1, 1, 1)
-axis.plot(x, y)
-axis.set_xlim(0, 2*l)
-axis.set_ylim(0, 1.2)
-axis.set_title('Kastenfunction --> FFT --> Sinc')
-
-y_fft = np.fft.fft(y)
-
-fig = plt.figure()
-axis = fig.add_subplot(1, 1, 1)
-axis.plot(x, np.abs(y_fft), 'k', label='abs')
-axis.plot(x, np.real(y_fft), 'r', label='real')
-axis.plot(x, np.imag(y_fft), 'b', label='imag')
-axis.set_xlim(0, 200)
-axis.legend()
-axis.set_title('Kastenfunction --> FFT --> Sinc')
-
-
-# Kastenfunktion mit Cos-"Softness":
-y = np.where(x <= l, 1, np.where(x >= l+w, 0, 0.5*(1+np.cos(pi*(x-l)/w))))
-
-fig = plt.figure()
-axis = fig.add_subplot(1, 1, 1)
-axis.plot(x, y)
-axis.set_xlim(0, 2*l)
-axis.set_ylim(0, 1.2)
-axis.set_title('Kastenfunction mit Cos-Kanten')
-
-y_fft = np.fft.fft(y)
-
-fig = plt.figure()
-axis = fig.add_subplot(1, 1, 1)
-axis.plot(x, np.abs(y_fft), 'k', label='abs')
-axis.plot(x, np.real(y_fft), 'r', label='real')
-axis.plot(x, np.imag(y_fft), 'b', label='imag')
-axis.set_xlim(0, 200)
-axis.legend()
-axis.set_title('Kastenfunction mit Cos-Kanten')
\ No newline at end of file
diff --git a/scripts/colaborations/vtk_interpolation.py b/scripts/colaborations/vtk_interpolation.py
deleted file mode 100644
index 3fea093..0000000
--- a/scripts/colaborations/vtk_interpolation.py
+++ /dev/null
@@ -1,73 +0,0 @@
-# -*- coding: utf-8 -*-
-"""
-Created on Fri Jan 17 13:09:08 2014
-
-@author: Jan
-"""
-
-from pylab import *
-import pickle
-from tqdm import tqdm
-from pyramid.magdata import MagData
-
-with open("vtk_to_numpy.pickle") as pf:
- data = pickle.load(pf)
-zs = unique(data[:,2])
-
-axis = plt.figure().add_subplot(1, 1, 1)
-axis.scatter(data[:, 0], data[:, 1])
-plt.show()
-
-# # regular grid
-xs = linspace(-8.5*2, 8.5*2, 256)
-ys = linspace(-9.5*2, 9.5*2, 256)
-
-o, p = meshgrid(xs, ys)
-
-newdata = zeros((len(xs), len(ys), len(zs), data.shape[1] - 3))
-
-
-## WITH MASKING OF THE CENTER:
-#
-#iz_x = concatenate([linspace(-4.95, -4.95, 50),
-# linspace(-4.95, 0, 50),
-# linspace(0, 4.95, 50),
-# linspace(4.95, 4.95, 50),
-# linspace(-4.95, 0, 50),
-# linspace(0, 4.95, 50),])
-#iz_y = concatenate([linspace(-2.88, 2.88, 50),
-# linspace(2.88, 5.7, 50),
-# linspace(5.7, 2.88, 50),
-# linspace(2.88, -2.88, 50),
-# linspace(-2.88, -5.7, 50),
-# linspace(-5.7, -2.88, 50), ])
-#
-#
-#for i, z in tqdm(enumerate(zs), total=len(zs)):
-# subdata = data[data[:, 2] == z, :]
-#
-# for j in range(newdata.shape[-1]):
-# gridded_subdata = griddata(concatenate([subdata[:, 0], iz_x]),
-# concatenate([subdata[:, 1], iz_y]), concatenate([subdata[:, 3 + j],
-# zeros(len(iz_x))]), o, p)
-# newdata[:, :, i, j] = gridded_subdata.filled(fill_value=0)
-
-
-# WITHOUT MASKING OF THE CENTER:
-
-
-for i, z in tqdm(enumerate(zs), total=len(zs)):
- subdata = data[data[:, 2] == z, :]
-
- for j in range(3): # For all 3 components!
- gridded_subdata = griddata(subdata[:, 0], subdata[:, 1], subdata[:, 3 + j], o, p)
- newdata[:, :, i, j] = gridded_subdata.filled(fill_value=0)
-
-
-magnitude = newdata.swapaxes(0,3).swapaxes(1,2).swapaxes(2,3)
-
-mag_data = MagData(1., magnitude)
-
-mag_data.quiver_plot()
-
-mag_data.save_to_netcdf4('vtk_mag_data.nc')
diff --git a/scripts/colaborations/vtk_to_numpy.py b/scripts/colaborations/vtk_to_numpy.py
deleted file mode 100644
index 0e88286..0000000
--- a/scripts/colaborations/vtk_to_numpy.py
+++ /dev/null
@@ -1,74 +0,0 @@
-# -*- coding: utf-8 -*-
-"""
-Created on Tue Jan 14 10:06:42 2014
-
-@author: Jan
-"""
-
-
-import numpy as np
-import vtk
-#import netCDF4
-import logging
-import sys
-
-import pickle
-
-logging.basicConfig(level=logging.INFO, format='%(levelname)s: %(message)s', stream=sys.stdout)
-log = logging.getLogger(__name__)
-
-reader = vtk.vtkDataSetReader()
-reader.SetFileName('irect_500x125x3.vtk')
-reader.ReadAllScalarsOn()
-reader.ReadAllVectorsOn()
-reader.Update()
-
-output = reader.GetOutput()
-size = output.GetNumberOfPoints()
-
-vtk_points = output.GetPoints().GetData()
-vtk_vectors = output.GetPointData().GetVectors()
-
-point_array = np.zeros(vtk_points.GetSize())
-vtk_points.ExportToVoidPointer(point_array)
-point_array = np.reshape(point_array, (-1, 3))
-
-vector_array = np.zeros(vtk_points.GetSize())
-vtk_vectors.ExportToVoidPointer(vector_array)
-vector_array = np.reshape(vector_array, (-1, 3))
-
-data = np.hstack((point_array, vector_array))
-
-log.info('Data reading complete!')
-
-#magfile = netCDF4.Dataset('tube_90x30x30.nc', 'w', format='NETCDF3_64BIT')
-#magfile.createDimension('comp', 6) # Number of components
-#magfile.createDimension('size', size)
-#
-#x = magfile.createVariable('x', 'f8', ('size'))
-#y = magfile.createVariable('y', 'f8', ('size'))
-#z = magfile.createVariable('z', 'f8', ('size'))
-#x_mag = magfile.createVariable('x_mag', 'f8', ('size'))
-#y_mag = magfile.createVariable('y_mag', 'f8', ('size'))
-#z_mag = magfile.createVariable('z_mag', 'f8', ('size'))
-#
-#log.info('Start saving data separately!')
-#x = data[:, 0]
-#y = data[:, 1]
-#z = data[:, 2]
-#x_mag = data[:, 3]
-#y_mag = data[:, 4]
-#z_mag = data[:, 5]
-#log.info('Separate saving complete!')
-#
-#log.info('Try saving the whole array!')
-#filedata = magfile.createVariable('data', 'f8', ('size', 'comp'))
-#filedata[:, :] = data
-#log.info('Saving complete!')
-#
-#magfile.close()
-
-log.info('Pickling data!')
-with open('vtk_to_numpy.pickle', 'w') as pf:
- pickle.dump(data, pf)
-log.info('Pickling complete!')
diff --git a/scripts/colaborations/edinburgh_test.py b/scripts/collaborations/edinburgh_test.py
similarity index 58%
rename from scripts/colaborations/edinburgh_test.py
rename to scripts/collaborations/edinburgh_test.py
index 42fcf27..c2a09f6 100644
--- a/scripts/colaborations/edinburgh_test.py
+++ b/scripts/collaborations/edinburgh_test.py
@@ -1,51 +1,48 @@
# -*- coding: utf-8 -*-
-"""
-Created on Tue Jan 28 15:15:08 2014
+"""Created on Tue Jan 28 15:15:08 2014 @author: Jan"""
-@author: Jan
-"""
+import os
import numpy as np
+
+from matplotlib.ticker import FuncFormatter
+
+import pyramid
from pyramid.magdata import MagData
from pyramid.projector import SimpleProjector
-from pyramid.phasemapper import PMAdapterFM, PMFourier
-from matplotlib.ticker import FuncFormatter
+from pyramid.phasemapper import PMAdapterFM
-data = np.loadtxt('../output/data from Edinburgh/long_grain_remapped_0p0035.txt', delimiter=',')
+import logging
+import logging.config
-a = 1000 * (data[1, 2] - data[0, 2])
-dim = len(np.unique(data[:, 2])), len(np.unique(data[:, 1])), len(np.unique(data[:, 0]))
+LOGGING_CONF = os.path.join(os.path.dirname(os.path.realpath(pyramid.__file__)), 'logging.ini')
-mag_vec = np.concatenate([data[:, 3], data[:, 4], data[:, 5]])
+logging.config.fileConfig(LOGGING_CONF, disable_existing_loggers=False)
+
+# Load data:
+data = np.loadtxt('../../output/data from Edinburgh/long_grain_remapped_0p0035.txt', delimiter=',')
+# Set parameters:
+a = 1000 * (data[1, 2] - data[0, 2])
+dim = len(np.unique(data[:, 2])), len(np.unique(data[:, 1])), len(np.unique(data[:, 0]))
+# Set magnetization:
+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)
-#
-#mag_data.scale_up()
-#
-#mag_data.quiver_plot()
-
-#mag_data.quiver_plot3d()
-
+# Pad and upscale:
+mag_data.pad(30, 20, 0)
+mag_data.scale_up()
+# Phasemapping:
projector = SimpleProjector(mag_data.dim)
-
phasemapper = PMAdapterFM(mag_data.a, projector)
-phasemapper = PMFourier(mag_data.a, projector, padding = 1)
-
phase_map = phasemapper(mag_data)
-
-phase_axis = phase_map.display_combined(density=20, interpolation='bilinear', grad_encode='bright')[0]
-
+# Plot:
+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)))
-
-#phase_map.display_phase3d()
diff --git a/scripts/colaborations/vadim_integration.py b/scripts/collaborations/vadim_integration.py
similarity index 70%
rename from scripts/colaborations/vadim_integration.py
rename to scripts/collaborations/vadim_integration.py
index 229313c..a1c9956 100644
--- a/scripts/colaborations/vadim_integration.py
+++ b/scripts/collaborations/vadim_integration.py
@@ -1,32 +1,36 @@
# -*- coding: utf-8 -*-
-"""
-Created on Thu Nov 07 16:47:52 2013
+"""Created on Thu Nov 07 16:47:52 2013 @author: Jan """
-@author: Jan
-"""
+
+import os
import numpy as np
from numpy import pi
-import os
-
+import pyramid
import pyramid.magcreator as mc
-import pyramid.projector as pj
-import pyramid.phasemapper as pm
-import pyramid.holoimage as hi
+from pyramid.projector import SimpleProjector
+from pyramid.phasemapper import PMConvolve, PMElectric
from pyramid.magdata import MagData
-from pyramid.phasemap import PhaseMap
import matplotlib.pyplot as plt
+import logging
+import logging.config
+
+
+LOGGING_CONF = os.path.join(os.path.dirname(os.path.realpath(pyramid.__file__)), 'logging.ini')
-Q_E = 1#1.602176565e-19
-EPSILON_0 = 1#8.8e-9
+
+logging.config.fileConfig(LOGGING_CONF, disable_existing_loggers=False)
+# Set constants:
+Q_E = 1 # 1.602176565e-19
+EPSILON_0 = 1 # 8.8e-9
C_E = 1
def calculate_charge_batch(phase_map):
- def calculate_charge(grad_y, grad_x, l, r, t, b):
+ def calculate_charge(grad_y, grad_x, l, r, t, b):
left = np.sum(-grad_x[b:t, l])
top = np.sum(grad_y[t, l:r])
right = np.sum(grad_x[b:t, r])
@@ -34,27 +38,28 @@ def calculate_charge_batch(phase_map):
integral = left + top + right + bottom
Q = -EPSILON_0/C_E * integral
return Q
- grad_y, grad_x = np.gradient(phase_map.phase, phase_map.a, phase_map.a)
+ grad_y, grad_x = np.gradient(phase_map.phase, phase_map.a, phase_map.a)
xi, yi = np.zeros(t-b), np.zeros(t-b)
for i, ti in enumerate(np.arange(b, t)):
xi[i] = ti
yi[i] = calculate_charge(grad_y, grad_x, l, r, ti, b)
return xi, yi
-directory = '../output/vadim/'
+
+directory = '../../output/vadim/'
if not os.path.exists(directory):
os.makedirs(directory)
-
+# Set parameters:
a = 1.0 # in nm
phi = pi/4
density = 30
dim = (128, 128, 128) # in px (z, y, x)
-
+v_0 = 1
+v_acc = 300000
l = dim[2]/4.
r = dim[2]/4. + dim[2]/2.
b = dim[1]/4.
t = dim[1]/4. + dim[1]/2.
-
# Create magnetic shape:
center = (dim[0]/2.-0.5, dim[1]/2.-0.5, dim[2]/2.-0.5) # in px (z, y, x) index starts at 0!
radius = dim[1]/8 # in px
@@ -63,29 +68,28 @@ width = (dim[0]/4, dim[1]/4, dim[2]/4)
mag_shape_sphere = mc.Shapes.sphere(dim, center, radius)
mag_shape_disc = mc.Shapes.disc(dim, center, radius, height)
mag_shape_slab = mc.Shapes.slab(dim, center, (height, radius*2, radius*2))
-mag_shape_ellipse = mc.Shapes.ellipsoid(dim, (center[0], 0*center[1], center[2]), (radius*5, radius*10, radius*2))
+mag_shape_ellipsoid = mc.Shapes.ellipsoid(dim, (center[0], 0*center[1], center[2]),
+ (radius*5, radius*10, radius*2))
# Create magnetization distributions:
mag_data_sphere = MagData(a, mc.create_mag_dist_homog(mag_shape_sphere, phi))
mag_data_disc = MagData(a, mc.create_mag_dist_homog(mag_shape_disc, phi))
mag_data_vortex = MagData(a, mc.create_mag_dist_vortex(mag_shape_disc))
mag_data_slab = MagData(a, mc.create_mag_dist_homog(mag_shape_slab, phi))
-mag_data_ellipse = MagData(a, mc.create_mag_dist_homog(mag_shape_ellipse, phi))
-# Project the magnetization data:
-projection_sphere = pj.simple_axis_projection(mag_data_sphere)
-projection_disc = pj.simple_axis_projection(mag_data_disc)
-projection_vortex = pj.simple_axis_projection(mag_data_vortex)
-projection_slab = pj.simple_axis_projection(mag_data_slab)
-projection_ellipse = pj.simple_axis_projection(mag_data_ellipse)
+mag_data_ellipsoid = MagData(a, mc.create_mag_dist_homog(mag_shape_ellipsoid, phi))
+# Create phasemapper:
+projector = SimpleProjector(dim)
+pm_mag = PMConvolve(a, projector)
+pm_ele = PMElectric(a, projector, v_0, v_acc)
# Magnetic phase map of a homogeneously magnetized disc:
-phase_map_mag_disc = PhaseMap(a, pm.phase_mag(a, projection_disc))
+phase_map_mag_disc = pm_mag(mag_data_disc)
phase_map_mag_disc.save_to_txt(directory+'phase_map_mag_disc.txt')
-axis, _ = hi.display_combined(phase_map_mag_disc, density)
+axis, _ = phase_map_mag_disc.display_combined(density=density)
axis.axvline(l, b/dim[1], t/dim[1], linewidth=2, color='g')
axis.axvline(r, b/dim[1], t/dim[1], linewidth=2, color='g')
axis.axhline(b, l/dim[2], r/dim[2], linewidth=2, color='g')
axis.axhline(t, l/dim[2], r/dim[2], linewidth=2, color='g')
axis.arrow(l+(r-l)/2, b, 0, t-b, length_includes_head=True,
- head_width=(r-l)/10, head_length=(r-l)/10, fc='g', ec='g')
+ head_width=(r-l)/10, head_length=(r-l)/10, fc='g', ec='g')
plt.savefig(directory+'phase_map_mag_disc.png')
x, y = calculate_charge_batch(phase_map_mag_disc)
np.savetxt(directory+'charge_integral_mag_disc.txt', np.array([x, y]).T)
@@ -93,15 +97,15 @@ plt.figure()
plt.plot(x, y)
plt.savefig(directory+'charge_integral_mag_disc.png')
# Magnetic phase map of a vortex state disc:
-phase_map_mag_vortex = PhaseMap(a, pm.phase_mag(a, projection_vortex))
+phase_map_mag_vortex = pm_mag(mag_data_vortex)
phase_map_mag_vortex.save_to_txt(directory+'phase_map_mag_vortex.txt')
-axis, _ = hi.display_combined(phase_map_mag_vortex, density)
+axis, _ = phase_map_mag_vortex.display_combined(density=density)
axis.axvline(l, b/dim[1], t/dim[1], linewidth=2, color='g')
axis.axvline(r, b/dim[1], t/dim[1], linewidth=2, color='g')
axis.axhline(b, l/dim[2], r/dim[2], linewidth=2, color='g')
axis.axhline(t, l/dim[2], r/dim[2], linewidth=2, color='g')
axis.arrow(l+(r-l)/2, b, 0, t-b, length_includes_head=True,
- head_width=(r-l)/10, head_length=(r-l)/10, fc='g', ec='g')
+ head_width=(r-l)/10, head_length=(r-l)/10, fc='g', ec='g')
plt.savefig(directory+'phase_map_mag_vortex.png')
x, y = calculate_charge_batch(phase_map_mag_vortex)
np.savetxt(directory+'charge_integral_mag_vortex.txt', np.array([x, y]).T)
@@ -109,15 +113,15 @@ plt.figure()
plt.plot(x, y)
plt.savefig(directory+'charge_integral_mag_vortex.png')
# MIP phase of a slab:
-phase_map_mip_slab = PhaseMap(a, pm.phase_elec(a, projection_slab, v_0=1, v_acc=300000))
+phase_map_mip_slab = pm_ele(mag_data_slab)
phase_map_mip_slab.save_to_txt(directory+'phase_map_mip_slab.txt')
-axis, _ = hi.display_combined(phase_map_mip_slab, density)
+axis, _ = phase_map_mip_slab.display_combined(density=density)
axis.axvline(l, b/dim[1], t/dim[1], linewidth=2, color='g')
axis.axvline(r, b/dim[1], t/dim[1], linewidth=2, color='g')
axis.axhline(b, l/dim[2], r/dim[2], linewidth=2, color='g')
axis.axhline(t, l/dim[2], r/dim[2], linewidth=2, color='g')
axis.arrow(l+(r-l)/2, b, 0, t-b, length_includes_head=True,
- head_width=(r-l)/10, head_length=(r-l)/10, fc='g', ec='g')
+ head_width=(r-l)/10, head_length=(r-l)/10, fc='g', ec='g')
plt.savefig(directory+'phase_map_mip_slab.png')
x, y = calculate_charge_batch(phase_map_mip_slab)
np.savetxt(directory+'charge_integral_mip_slab.txt', np.array([x, y]).T)
@@ -125,15 +129,15 @@ plt.figure()
plt.plot(x, y)
plt.savefig(directory+'charge_integral_mip_slab.png')
# MIP phase of a disc:
-phase_map_mip_disc = PhaseMap(a, pm.phase_elec(a, projection_disc, v_0=1, v_acc=300000))
+phase_map_mip_disc = pm_ele(mag_data_disc)
phase_map_mip_disc.save_to_txt(directory+'phase_map_mip_disc.txt')
-axis, _ = hi.display_combined(phase_map_mip_disc, density)
+axis, _ = phase_map_mip_disc.display_combined(density=density)
axis.axvline(l, b/dim[1], t/dim[1], linewidth=2, color='g')
axis.axvline(r, b/dim[1], t/dim[1], linewidth=2, color='g')
axis.axhline(b, l/dim[2], r/dim[2], linewidth=2, color='g')
axis.axhline(t, l/dim[2], r/dim[2], linewidth=2, color='g')
axis.arrow(l+(r-l)/2, b, 0, t-b, length_includes_head=True,
- head_width=(r-l)/10, head_length=(r-l)/10, fc='g', ec='g')
+ head_width=(r-l)/10, head_length=(r-l)/10, fc='g', ec='g')
plt.savefig(directory+'phase_map_mip_disc.png')
x, y = calculate_charge_batch(phase_map_mip_disc)
np.savetxt(directory+'charge_integral_mip_disc.txt', np.array([x, y]).T)
@@ -141,15 +145,15 @@ plt.figure()
plt.plot(x, y)
plt.savefig(directory+'charge_integral_mip_disc.png')
# MIP phase of a sphere:
-phase_map_mip_sphere = PhaseMap(a, pm.phase_elec(a, projection_sphere, v_0=1, v_acc=300000))
+phase_map_mip_sphere = pm_ele(mag_data_sphere)
phase_map_mip_sphere.save_to_txt(directory+'phase_map_mip_sphere.txt')
-axis, _ = hi.display_combined(phase_map_mip_sphere, density)
+axis, _ = phase_map_mip_sphere.display_combined(density=density)
axis.axvline(l, b/dim[1], t/dim[1], linewidth=2, color='g')
axis.axvline(r, b/dim[1], t/dim[1], linewidth=2, color='g')
axis.axhline(b, l/dim[2], r/dim[2], linewidth=2, color='g')
axis.axhline(t, l/dim[2], r/dim[2], linewidth=2, color='g')
axis.arrow(l+(r-l)/2, b, 0, t-b, length_includes_head=True,
- head_width=(r-l)/10, head_length=(r-l)/10, fc='g', ec='g')
+ head_width=(r-l)/10, head_length=(r-l)/10, fc='g', ec='g')
plt.savefig(directory+'phase_map_mip_sphere.png')
x, y = calculate_charge_batch(phase_map_mip_sphere)
np.savetxt(directory+'charge_integral_mip_sphere.txt', np.array([x, y]).T)
@@ -157,15 +161,15 @@ plt.figure()
plt.plot(x, y)
plt.savefig(directory+'charge_integral_mip_sphere.png')
# MIP phase of an ellipsoid:
-phase_map_mip_ellipsoid = PhaseMap(a, pm.phase_elec(a, projection_ellipse, v_0=1, v_acc=300000))
+phase_map_mip_ellipsoid = pm_ele(mag_data_ellipsoid)
phase_map_mip_ellipsoid.save_to_txt(directory+'phase_map_mip_ellipsoid.txt')
-axis, _ = hi.display_combined(phase_map_mip_ellipsoid, density)
+axis, _ = phase_map_mip_ellipsoid.display_combined(phase_map_mip_ellipsoid, density)
axis.axvline(l, b/dim[1], t/dim[1], linewidth=2, color='g')
axis.axvline(r, b/dim[1], t/dim[1], linewidth=2, color='g')
axis.axhline(b, l/dim[2], r/dim[2], linewidth=2, color='g')
axis.axhline(t, l/dim[2], r/dim[2], linewidth=2, color='g')
axis.arrow(l+(r-l)/2, b, 0, t-b, length_includes_head=True,
- head_width=(r-l)/10, head_length=(r-l)/10, fc='g', ec='g')
+ head_width=(r-l)/10, head_length=(r-l)/10, fc='g', ec='g')
plt.savefig(directory+'phase_map_mip_ellipsoid.png')
x, y = calculate_charge_batch(phase_map_mip_ellipsoid)
np.savetxt(directory+'charge_integral_mip_ellipsoid.txt', np.array([x, y]).T)
diff --git a/scripts/colaborations/vtk_conversion.py b/scripts/collaborations/vtk_conversion.py
similarity index 84%
rename from scripts/colaborations/vtk_conversion.py
rename to scripts/collaborations/vtk_conversion.py
index bbbb506..ea5e99b 100644
--- a/scripts/colaborations/vtk_conversion.py
+++ b/scripts/collaborations/vtk_conversion.py
@@ -1,36 +1,38 @@
# -*- coding: utf-8 -*-
-"""
-Created on Fri Jan 24 11:17:11 2014
+"""Created on Fri Jan 24 11:17:11 2014 @author: Jan"""
-@author: Jan
-"""
+import os
import numpy as np
-import vtk
-import logging
-import os
-import sys
+import matplotlib.pyplot as plt
+from pylab import griddata
+
import pickle
+import vtk
from tqdm import tqdm
+
+import pyramid
from pyramid.magdata import MagData
-import matplotlib.pyplot as plt
-from pylab import griddata
from pyramid.projector import SimpleProjector
from pyramid.phasemapper import PMAdapterFM
+import logging
+import logging.config
+
+
+LOGGING_CONF = os.path.join(os.path.dirname(os.path.realpath(pyramid.__file__)), 'logging.ini')
+
+
+logging.config.fileConfig(LOGGING_CONF, disable_existing_loggers=False)
###################################################################################################
-PATH = '../output/vtk data/tube_90x30x50_sat_edge_equil.gmr'
+PATH = '../../output/vtk data/tube_90x30x50_sat_edge_equil.gmr'
b_0 = 1.54
gain = 12
force_calculation = False
###################################################################################################
-
-logging.basicConfig(level=logging.INFO, format='%(levelname)s: %(message)s', stream=sys.stdout)
-log = logging.getLogger(__name__)
-
+# Load vtk-data:
if force_calculation or not os.path.exists(PATH+'.pickle'):
- log.info('Reading data from vtk-file')
# Setting up reader:
reader = vtk.vtkDataSetReader()
reader.SetFileName(PATH+'.vtk')
@@ -53,25 +55,20 @@ if force_calculation or not os.path.exists(PATH+'.pickle'):
vector_array = np.reshape(vector_array, (-1, 3))
# Combining data:
data = np.hstack((point_array, vector_array))
- log.info('Data reading complete!')
- log.info('Pickling data!')
with open(PATH+'.pickle', 'w') as pf:
pickle.dump(data, pf)
- log.info('Pickling complete!')
else:
- log.info('Loading pickled data!')
with open(PATH+'.pickle') as pf:
data = pickle.load(pf)
- log.info('Loading complete!')
# Scatter plot of all x-y-coordinates
axis = plt.figure().add_subplot(1, 1, 1)
axis.scatter(data[:, 0], data[:, 1])
plt.show()
-
+###################################################################################################
+# Interpolate on regular grid:
if force_calculation or not os.path.exists(PATH+'.nc'):
- log.info('Arranging data in z-slices!')
# Find unique z-slices:
- zs = np.unique(data[:, 2])
+ zs = np.unique(data[:, 2])
# Determine the grid spacing:
a = zs[1] - zs[0]
# Determine the size of object:
@@ -81,19 +78,19 @@ if force_calculation or not os.path.exists(PATH+'.nc'):
x_diff, y_diff, z_diff = np.ptp(data[:, 0]), np.ptp(data[:, 1]), np.ptp(data[:, 2])
x_cent, y_cent, z_cent = x_min+x_diff/2., y_min+y_diff/2., z_min+z_diff/2.
# Create regular grid
- xs = np.arange(x_cent-x_diff, x_cent+x_diff, a) #linspace(-8.5*2, 8.5*2, 256)
- ys = np.arange(y_cent-y_diff, y_cent+y_diff, a) #linspace(-9.5*2, 9.5*2, 256)
+ xs = np.arange(x_cent-x_diff, x_cent+x_diff, a)
+ ys = np.arange(y_cent-y_diff, y_cent+y_diff, a)
o, p = np.meshgrid(xs, ys)
# Create empty magnitude:
magnitude = np.zeros((3, len(zs), len(ys), len(xs)))
-
+
# WITH MASKING OF THE CENTER (SYMMETRIC):
iz_x = np.concatenate([np.linspace(-4.95, -4.95, 50),
np.linspace(-4.95, 0, 50),
np.linspace(0, 4.95, 50),
np.linspace(4.95, 4.95, 50),
np.linspace(-4.95, 0, 50),
- np.linspace(0, 4.95, 50),])
+ np.linspace(0, 4.95, 50), ])
iz_y = np.concatenate([np.linspace(-2.88, 2.88, 50),
np.linspace(2.88, 5.7, 50),
np.linspace(5.7, 2.88, 50),
@@ -104,7 +101,7 @@ if force_calculation or not os.path.exists(PATH+'.nc'):
subdata = data[data[:, 2] == z, :]
for j in range(3): # For all 3 components!
gridded_subdata = griddata(np.concatenate([subdata[:, 0], iz_x]),
- np.concatenate([subdata[:, 1], iz_y]),
+ np.concatenate([subdata[:, 1], iz_y]),
np.concatenate([subdata[:, 3 + j], np.zeros(len(iz_x))]),
o, p)
magnitude[j, i, :, :] = gridded_subdata.filled(fill_value=0)
@@ -126,11 +123,11 @@ if force_calculation or not os.path.exists(PATH+'.nc'):
# subdata = data[data[:, 2] == z, :]
# for j in range(3): # For all 3 components!
# gridded_subdata = griddata(np.concatenate([subdata[:, 0], iz_x]),
-# np.concatenate([subdata[:, 1], iz_y]),
+# np.concatenate([subdata[:, 1], iz_y]),
# np.concatenate([subdata[:, 3 + j], np.zeros(len(iz_x))]),
# o, p)
# magnitude[j, i, :, :] = gridded_subdata.filled(fill_value=0)
-
+
# # WITHOUT MASKING OF THE CENTER:
# for i, z in tqdm(enumerate(zs), total=len(zs)):
# subdata = data[data[:, 2] == z, :]
@@ -141,13 +138,11 @@ if force_calculation or not os.path.exists(PATH+'.nc'):
# Creating MagData object:
mag_data = MagData(0.2*10, magnitude)
mag_data.save_to_netcdf4(PATH+'.nc')
- log.info('MagData created!')
else:
- log.info('Loading MagData!')
mag_data = MagData.load_from_netcdf4(PATH+'.nc')
- log.info('Loading complete!')
mag_data.quiver_plot()
-
+###################################################################################################
+# Phasemapping:
projector = SimpleProjector(mag_data.dim)
phasemapper = PMAdapterFM(mag_data.a, projector)
phase_map = phasemapper(mag_data)
diff --git a/scripts/create distributions/create_core_shell_sphere.py b/scripts/create distributions/create_core_shell_sphere.py
index 321d590..9c99de5 100644
--- a/scripts/create distributions/create_core_shell_sphere.py
+++ b/scripts/create distributions/create_core_shell_sphere.py
@@ -1,6 +1,6 @@
#! python
# -*- coding: utf-8 -*-
-"""Create a core-shell disc."""
+"""Create a core-shell sphere."""
import os
@@ -21,14 +21,15 @@ import logging.config
LOGGING_CONF = os.path.join(os.path.dirname(os.path.realpath(pyramid.__file__)), 'logging.ini')
+
logging.config.fileConfig(LOGGING_CONF, disable_existing_loggers=False)
directory = '../../output/magnetic distributions'
if not os.path.exists(directory):
os.makedirs(directory)
# Input parameters:
filename = directory + '/mag_dist_core_shell_sphere.txt'
-a = 1.0 # in nm
-density = 1
+a = 50.0 # in nm
+density = 500
dim = (32, 32, 32)
center = (dim[0]/2-0.5, int(dim[1]/2)-0.5, int(dim[2]/2)-0.5)
radius_core = dim[1]/8
diff --git a/scripts/create distributions/create_multiple_samples.py b/scripts/create distributions/create_multiple_samples.py
index d1f4069..d791960 100644
--- a/scripts/create distributions/create_multiple_samples.py
+++ b/scripts/create distributions/create_multiple_samples.py
@@ -1,6 +1,6 @@
#! python
# -*- coding: utf-8 -*-
-"""Create random magnetic distributions."""
+"""Create multiple magnetic distributions."""
import os
diff --git a/scripts/create distributions/create_vortex.py b/scripts/create distributions/create_vortex.py
index e8ae9fd..cd70c5f 100644
--- a/scripts/create distributions/create_vortex.py
+++ b/scripts/create distributions/create_vortex.py
@@ -1,6 +1,6 @@
#! python
# -*- coding: utf-8 -*-
-"""Create the Pyramid-Logo."""
+"""Create a magnetic vortex distribution."""
import os
diff --git a/scripts/publications/abstract_prague.py b/scripts/publications/abstract_prague.py
index cb20acd..cb5bfe7 100644
--- a/scripts/publications/abstract_prague.py
+++ b/scripts/publications/abstract_prague.py
@@ -1,21 +1,29 @@
# -*- coding: utf-8 -*-
-"""
-Created on Mon Mar 17 14:28:10 2014
+"""Created on Mon Mar 17 14:28:10 2014 @author: Jan"""
-@author: Jan
-"""
+import os
from numpy import pi
import numpy as np
+
import matplotlib.pyplot as plt
from matplotlib.ticker import FuncFormatter
+import pyramid
import pyramid.magcreator as mc
from pyramid.magdata import MagData
from pyramid.projector import SimpleProjector, YTiltProjector
from pyramid.phasemapper import PMConvolve
+import logging
+import logging.config
+
+
+LOGGING_CONF = os.path.join(os.path.dirname(os.path.realpath(pyramid.__file__)), 'logging.ini')
+
+
+logging.config.fileConfig(LOGGING_CONF, disable_existing_loggers=False)
###################################################################################################
print('Jan')
@@ -112,4 +120,3 @@ axis.tick_params(axis='both', which='major', labelsize=18)
axis.set_title('Phase map', fontsize=24)
axis.set_xlabel('x [nm]', fontsize=18)
axis.set_ylabel('y [nm]', fontsize=18)
-
diff --git a/scripts/publications/paper 1/ch5-0-evaluation_and_comparison.py b/scripts/publications/paper 1/ch5-0-evaluation_and_comparison.py
index 00463cc..4465481 100644
--- a/scripts/publications/paper 1/ch5-0-evaluation_and_comparison.py
+++ b/scripts/publications/paper 1/ch5-0-evaluation_and_comparison.py
@@ -1,9 +1,5 @@
# -*- coding: utf-8 -*-
-"""
-Created on Fri Jul 26 14:37:20 2013
-
-@author: Jan
-"""
+"""Created on Fri Jul 26 14:37:20 2013 @author: Jan"""
import os
@@ -11,18 +7,25 @@ import os
import numpy as np
from numpy import pi
-import shelve
+import matplotlib.pyplot as plt
+from matplotlib.ticker import FixedFormatter, IndexLocator
+import pyramid
import pyramid.magcreator as mc
from pyramid.magdata import MagData
-import matplotlib.pyplot as plt
+import shelve
-from matplotlib.ticker import FixedFormatter, IndexLocator
+import logging
+import logging.config
-force_calculation = False
+LOGGING_CONF = os.path.join(os.path.dirname(os.path.realpath(pyramid.__file__)), 'logging.ini')
+
+logging.config.fileConfig(LOGGING_CONF, disable_existing_loggers=False)
+
+force_calculation = False
print '\nACCESS SHELVE'
# Create / Open databank:
@@ -40,8 +43,8 @@ y_r = x/np.abs(x)**3
y_k = x/np.abs(x)**2
fig = plt.figure()
axis = fig.add_subplot(1, 1, 1, aspect='equal')
-axis.plot(x,y_r, 'r', label=r'$r/|r|^3$')
-axis.plot(x,y_k, 'b', label=r'$k/|k|^2$')
+axis.plot(x, y_r, 'r', label=r'$r/|r|^3$')
+axis.plot(x, y_k, 'b', label=r'$k/|k|^2$')
axis.set_xlim(-5, 5)
axis.set_ylim(-5, 5)
axis.axvline(0, linewidth=2, color='k')
diff --git a/scripts/publications/paper 1/ch5-1-evaluation_real_space.py b/scripts/publications/paper 1/ch5-1-evaluation_real_space.py
index d3a618a..3f0b350 100644
--- a/scripts/publications/paper 1/ch5-1-evaluation_real_space.py
+++ b/scripts/publications/paper 1/ch5-1-evaluation_real_space.py
@@ -1,9 +1,5 @@
# -*- coding: utf-8 -*-
-"""
-Created on Fri Jul 26 14:37:30 2013
-
-@author: Jan
-"""
+"""Created on Fri Jul 26 14:37:30 2013 @author: Jan"""
import os
@@ -11,8 +7,13 @@ import os
import numpy as np
from numpy import pi
-import shelve
+import matplotlib.pyplot as plt
+from matplotlib.colors import BoundaryNorm
+from matplotlib.ticker import MaxNLocator
+from matplotlib.cm import RdBu
+from matplotlib.patches import Rectangle
+import pyramid
import pyramid.magcreator as mc
import pyramid.analytic as an
from pyramid.projector import SimpleProjector
@@ -20,17 +21,20 @@ from pyramid.phasemapper import PMConvolve
from pyramid.magdata import MagData
from pyramid.phasemap import PhaseMap
-import matplotlib.pyplot as plt
-from matplotlib.colors import BoundaryNorm
-from matplotlib.ticker import MaxNLocator
-from matplotlib.cm import RdBu
-from matplotlib.patches import Rectangle
+import shelve
+import logging
+import logging.config
-force_calculation = False
+
+LOGGING_CONF = os.path.join(os.path.dirname(os.path.realpath(pyramid.__file__)), 'logging.ini')
PHI_0 = -2067.83 # magnetic flux in T*nm²
+logging.config.fileConfig(LOGGING_CONF, disable_existing_loggers=False)
+
+force_calculation = False
+
print '\nACCESS SHELVE'
# Create / Open databank:
directory = '../../output/paper 1'
@@ -126,7 +130,8 @@ else:
slice_pos = int(mag_data_disc.dim[1]/2)
y_d.append(phase_map.phase[slice_pos, :]*1E3) # *1E3: rad to mrad
dy_d.append(phase_map.phase[slice_pos, :]*1E3 - F_disc(x_d[-1])) # *1E3: rad to mrad
- if i < 4: mag_data_disc.scale_down()
+ if i < 4:
+ mag_data_disc.scale_down()
print '--CREATE PHASE SLICES VORTEX STATE DISC'
x_v = []
@@ -150,7 +155,7 @@ else:
mag_data_vort = MagData(a, mc.create_mag_dist_vortex(mag_shape))
for i in range(5):
print '----a =', mag_data_vort.a, 'nm', 'dim =', mag_data_vort.dim
- projector = SimpleProjector(mag_data_vort.dim)
+ projector = SimpleProjector(mag_data_vort.dim)
phase_map = PMConvolve(mag_data_vort.a, projector)(mag_data_vort)
phase_map.display_combined(density, 'Disc, a = {} nm'.format(mag_data_vort.a))
x_v.append(np.linspace(mag_data_vort.a * 0.5,
@@ -159,7 +164,8 @@ else:
slice_pos = int(mag_data_vort.dim[1]/2)
y_v.append(phase_map.phase[slice_pos, :]*1E3) # *1E3: rad to mrad
dy_v.append(phase_map.phase[slice_pos, :]*1E3 - F_vort(x_v[-1])) # *1E3: rad to mrad
- if i < 4: mag_data_vort.scale_down()
+ if i < 4:
+ mag_data_vort.scale_down()
# Shelve x, y and dy:
print '--SAVE PHASE SLICES'
@@ -188,7 +194,7 @@ zoom = (23.5, 160, 15, 40)
rect = Rectangle((zoom[0], zoom[1]), zoom[2], zoom[3], fc='w', ec='k')
axes[0].add_patch(rect)
axes[0].arrow(zoom[0]+zoom[2], zoom[1]+zoom[3]/2, 36, 0, length_includes_head=True,
- head_width=10, head_length=4, fc='k', ec='k')
+ head_width=10, head_length=4, fc='k', ec='k')
# Plot zoom inset:
ins_axis_d = plt.axes([0.33, 0.57, 0.14, 0.3])
ins_axis_d.plot(x_d[0], y_d[0], '-k', linewidth=1.5, label='analytic')
@@ -223,7 +229,7 @@ zoom = (59, 340, 10, 55)
rect = Rectangle((zoom[0], zoom[1]), zoom[2], zoom[3], fc='w', ec='k')
axes[1].add_patch(rect)
axes[1].arrow(zoom[0]+zoom[2]/2, zoom[1], 0, -193, length_includes_head=True,
- head_width=2, head_length=20, fc='k', ec='k')
+ head_width=2, head_length=20, fc='k', ec='k')
# Plot zoom inset:
ins_axis_v = plt.axes([0.695, 0.15, 0.075, 0.3])
ins_axis_v.plot(x_v[0], y_v[0], '-k', linewidth=1.5, label='analytic')
diff --git a/scripts/publications/paper 1/ch5-2-evaluation_fourier_space.py b/scripts/publications/paper 1/ch5-2-evaluation_fourier_space.py
index 53838f2..e2460f6 100644
--- a/scripts/publications/paper 1/ch5-2-evaluation_fourier_space.py
+++ b/scripts/publications/paper 1/ch5-2-evaluation_fourier_space.py
@@ -1,35 +1,39 @@
# -*- coding: utf-8 -*-
-"""
-Created on Fri Jul 26 14:37:30 2013
+"""Created on Fri Jul 26 14:37:30 2013 @author: Jan"""
-@author: Jan
-"""
-
-import time
import os
import numpy as np
from numpy import pi
-import shelve
+import matplotlib.pyplot as plt
+from matplotlib.colors import BoundaryNorm
+from matplotlib.ticker import MaxNLocator
+from matplotlib.cm import RdBu
+import pyramid
import pyramid.magcreator as mc
import pyramid.analytic as an
from pyramid.magdata import MagData
from pyramid.projector import SimpleProjector
from pyramid.phasemapper import PMFourier
-import matplotlib.pyplot as plt
-from matplotlib.colors import BoundaryNorm
-from matplotlib.ticker import MaxNLocator
-from matplotlib.cm import RdBu
+import time
+import shelve
+import logging
+import logging.config
-force_calculation = False
+
+LOGGING_CONF = os.path.join(os.path.dirname(os.path.realpath(pyramid.__file__)), 'logging.ini')
PHI_0 = -2067.83 # magnetic flux in T*nm²
+logging.config.fileConfig(LOGGING_CONF, disable_existing_loggers=False)
+
+force_calculation = False
+
print '\nACCESS SHELVE'
# Create / Open databank:
directory = '../../output/paper 1'
@@ -101,7 +105,8 @@ else:
y_d10.append(phase_map10.phase[slice_pos, :]*1E3) # *1E3: rad to mrad
dy_d0.append(phase_map0.phase[slice_pos, :]*1E3 - F_disc(x_d[-1])) # *1E3: in mrad
dy_d10.append(phase_map10.phase[slice_pos, :]*1E3 - F_disc(x_d[-1])) # *1E3: in mrad
- if i < 4: mag_data_disc.scale_down()
+ if i < 4:
+ mag_data_disc.scale_down()
print '--CREATE PHASE SLICES VORTEX STATE DISC'
x_v = []
@@ -147,7 +152,8 @@ else:
y_v10.append(phase_map10.phase[slice_pos, :]*1E3) # *1E3: rad to mrad
dy_v0.append(phase_map0.phase[slice_pos, :]*1E3 - F_vort(x_v[-1])) # *1E3: in mrad
dy_v10.append(phase_map10.phase[slice_pos, :]*1E3 - F_vort(x_v[-1])) # *1E3: in mrad
- if i < 4: mag_data_vort.scale_down()
+ if i < 4:
+ mag_data_vort.scale_down()
# Shelve x, y and dy:
print '--SAVE PHASE SLICES'
diff --git a/scripts/publications/paper 1/ch5-3-comparison_of_methods.py b/scripts/publications/paper 1/ch5-3-comparison_of_methods.py
index c0c64b4..e32b90d 100644
--- a/scripts/publications/paper 1/ch5-3-comparison_of_methods.py
+++ b/scripts/publications/paper 1/ch5-3-comparison_of_methods.py
@@ -1,13 +1,7 @@
-#! python
# -*- coding: utf-8 -*-
-"""
-Created on Fri Jul 26 14:37:30 2013
+"""Created on Fri Jul 26 14:37:30 2013 @author: Jan"""
-@author: Jan
-"""
-
-import time
import os
import numpy as np
@@ -16,13 +10,24 @@ from numpy import pi
import matplotlib.pyplot as plt
from matplotlib.ticker import IndexLocator
+import pyramid
import pyramid.magcreator as mc
import pyramid.projector as pj
import pyramid.phasemapper as pm
import pyramid.analytic as an
from pyramid.magdata import MagData
+
+import time
import shelve
+import logging
+import logging.config
+
+
+LOGGING_CONF = os.path.join(os.path.dirname(os.path.realpath(pyramid.__file__)), 'logging.ini')
+
+
+logging.config.fileConfig(LOGGING_CONF, disable_existing_loggers=False)
print '\nACCESS SHELVE'
# Create / Open databank:
diff --git a/scripts/publications/paper 1/ch5-4-comparison_of_methods_new.py b/scripts/publications/paper 1/ch5-4-comparison_of_methods_new.py
index e976c83..b000976 100644
--- a/scripts/publications/paper 1/ch5-4-comparison_of_methods_new.py
+++ b/scripts/publications/paper 1/ch5-4-comparison_of_methods_new.py
@@ -1,33 +1,35 @@
-#! python
# -*- coding: utf-8 -*-
-"""
-Created on Fri Jul 26 14:37:30 2013
+"""Created on Fri Jul 26 14:37:30 2013 @author: Jan"""
-@author: Jan
-"""
-
-import time
import os
import numpy as np
from numpy import pi
-
import matplotlib.pyplot as plt
from matplotlib.ticker import NullLocator, LogLocator, LogFormatter
+import pyramid
import pyramid.magcreator as mc
import pyramid.analytic as an
from pyramid.magdata import MagData
from pyramid.projector import SimpleProjector
from pyramid.phasemapper import PMConvolve, PMFourier
+import time
import shelve
+import logging
+import logging.config
-force_calculation = True
+
+LOGGING_CONF = os.path.join(os.path.dirname(os.path.realpath(pyramid.__file__)), 'logging.ini')
+logging.config.fileConfig(LOGGING_CONF, disable_existing_loggers=False)
+
+force_calculation = True
+
print '\nACCESS SHELVE'
# Create / Open databank:
directory = '../../output/paper 1'
@@ -86,7 +88,7 @@ else:
height = dim[0]/2 # in px
print '--a =', a, 'nm', 'dim =', dim
-
+
# Create projector along z-axis and phasemapper:
projector = SimpleProjector(dim)
pm_fourier0 = PMFourier(a, projector, padding=0)
@@ -134,7 +136,8 @@ else:
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 '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:
@@ -272,7 +275,7 @@ axes[0, 1].grid()
# Add legend:
axes[1, 1].legend(bbox_to_anchor=(0, 0, 0.955, 0.615), bbox_transform=fig.transFigure,
- prop={'size':12})
+ prop={'size': 12})
# Save figure as .png:
plt.show()
@@ -286,5 +289,3 @@ plt.savefig(directory + '/ch5-3-method comparison.png', bbox_inches='tight')
print 'CLOSING SHELVE\n'
# Close shelve:
data_shelve.close()
-
-###############################################################################################
diff --git a/scripts/publications/poster_pof.py b/scripts/publications/poster_pof.py
new file mode 100644
index 0000000..24332f9
--- /dev/null
+++ b/scripts/publications/poster_pof.py
@@ -0,0 +1,63 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Wed Feb 05 19:19:19 2014
+
+@author: Jan
+"""
+
+import os
+
+from numpy import pi
+import numpy as np
+
+import pyramid
+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
+
+from time import clock
+
+import logging
+import logging.config
+
+
+LOGGING_CONF = os.path.join(os.path.dirname(os.path.realpath(pyramid.__file__)), 'logging.ini')
+
+
+logging.config.fileConfig(LOGGING_CONF, disable_existing_loggers=False)
+
+
+a = 10.0 # nm
+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()
+
+
+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))
+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/publications/poster_pov.py b/scripts/publications/poster_pov.py
deleted file mode 100644
index 2135cbe..0000000
--- a/scripts/publications/poster_pov.py
+++ /dev/null
@@ -1,91 +0,0 @@
-# -*- coding: utf-8 -*-
-"""
-Created on Wed Feb 05 19:19:19 2014
-
-@author: Jan
-"""
-
-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'))
-#
-#
-#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()
-
-
-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))
-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/reconstruction/reconstruction_sparse_cg_singular_values.py b/scripts/reconstruction/reconstruction_sparse_cg_singular_values.py
new file mode 100644
index 0000000..6917161
--- /dev/null
+++ b/scripts/reconstruction/reconstruction_sparse_cg_singular_values.py
@@ -0,0 +1,74 @@
+# -*- coding: utf-8 -*-
+"""Created on Thu Apr 24 11:00:00 2014 @author: Jan"""
+
+
+import os
+
+import numpy as np
+from numpy import pi
+
+import pyramid
+from pyramid.magdata import MagData
+from pyramid.projector import YTiltProjector, XTiltProjector
+from pyramid.dataset import DataSet
+from pyramid.kernel import Kernel
+from pyramid.phasemap import PhaseMap
+from pyramid.forwardmodel import ForwardModel
+
+import logging
+import logging.config
+
+
+LOGGING_CONF = os.path.join(os.path.dirname(os.path.realpath(pyramid.__file__)), 'logging.ini')
+
+
+logging.config.fileConfig(LOGGING_CONF, disable_existing_loggers=False)
+
+###################################################################################################
+print('--Singular value decomposition')
+
+a = 1.
+b_0 = 1.
+dim = (3, 3, 3)
+dim_uv = dim[1:3]
+count = 8
+
+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)
+
+phase_zero = PhaseMap(a, np.zeros(dim_uv))
+
+projectors_xy_full, projectors_x_full, projectors_y_full = [], [], []
+projectors_xy_miss, projectors_x_miss, projectors_y_miss = [], [], []
+projectors_x_full.extend([XTiltProjector(dim, tilt) for tilt in tilts_full])
+projectors_y_full.extend([YTiltProjector(dim, tilt) for tilt in tilts_full])
+projectors_xy_full.extend(projectors_x_full)
+projectors_xy_full.extend(projectors_y_full)
+projectors_x_miss.extend([XTiltProjector(dim, tilt) for tilt in tilts_miss])
+projectors_y_miss.extend([YTiltProjector(dim, tilt) for tilt in tilts_miss])
+projectors_xy_miss.extend(projectors_x_miss)
+projectors_xy_miss.extend(projectors_y_miss)
+
+###################################################################################################
+projectors = projectors_xy_full
+###################################################################################################
+
+size_2d = np.prod(dim_uv)
+size_3d = np.prod(dim)
+
+data = DataSet(a, dim_uv, b_0)
+[data.append((phase_zero, projectors[i])) for i in range(len(projectors))]
+
+y = data.phase_vec
+kern = Kernel(data.a, data.dim_uv, data.b_0)
+F = ForwardModel(data.projectors, kern)
+
+M = np.asmatrix([F.jac_dot(None, np.eye(3*size_3d)[:, i]) for i in range(3*size_3d)]).T
+#MTM = M.T * M + lam * np.asmatrix(np.eye(3*size_3d))
+
+U, s, V = np.linalg.svd(M) # TODO: M or MTM?
+
+for i in range(len(s)):
+ print 'Singular value:', s[i]
+ title = 'Singular value = {:g}'.format(s[i])
+ MagData(data.a, np.array(V[i, :]).reshape((3,)+dim)).quiver_plot3d(title)
diff --git a/scripts/reconstruction/reconstruction_sparse_cg_test.py b/scripts/reconstruction/reconstruction_sparse_cg_test.py
index 45385f9..3704548 100644
--- a/scripts/reconstruction/reconstruction_sparse_cg_test.py
+++ b/scripts/reconstruction/reconstruction_sparse_cg_test.py
@@ -1,14 +1,13 @@
# -*- coding: utf-8 -*-
-"""
-Created on Fri Feb 28 14:25:59 2014
+"""Created on Fri Feb 28 14:25:59 2014 @author: Jan"""
-@author: Jan
-"""
+import os
import numpy as np
from numpy import pi
+import pyramid
from pyramid.magdata import MagData
from pyramid.projector import YTiltProjector, XTiltProjector
from pyramid.phasemapper import PMConvolve
@@ -20,192 +19,60 @@ from pyramid.forwardmodel import ForwardModel
from pyramid.costfunction import Costfunction
import pyramid.reconstruction as rc
-from scipy.sparse.linalg import cg, LinearOperator
+import logging
+import logging.config
-###################################################################################################
-print('--Compliance test for one projection')
-
-a = 1.
-b_0 = 1000.
-dim = (2, 2, 2)
-count = 1
-
-magnitude = np.zeros((3,)+dim)
-magnitude[:, int(dim[0]/2), int(dim[1]/2), int(dim[2]/2)] = 1.
-
-mag_data = MagData(a, magnitude)
-
-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]
-
-dim_uv = dim[1:3]
-
-data_collection = DataSet(a, dim_uv, b_0)
-
-[data_collection.append((phase_maps[i], projectors[i])) for i in range(count)]
+LOGGING_CONF = os.path.join(os.path.dirname(os.path.realpath(pyramid.__file__)), 'logging.ini')
-data = data_collection
-y = data.phase_vec
-kern = Kernel(data.a, data.dim_uv, data.b_0)
-F = ForwardModel(data.projectors, kern)
-C = Costfunction(y, F)
-
-size_2d = np.prod(dim[1]*dim[2])
-size_3d = np.prod(dim)
-
-P = np.array([data.projectors[0](np.eye(3*size_3d)[:, i]) for i in range(3*size_3d)]).T
-K = np.array([kern(np.eye(2*size_2d)[:, i]) for i in range(2*size_2d)]).T
-F_mult = K.dot(P)
-F_direct = np.array([F(np.eye(3*size_3d)[:, i]) for i in range(3*size_3d)]).T
-np.testing.assert_almost_equal(F_direct, F_mult)
-
-P_jac = np.array([data.projectors[0].jac_dot(np.eye(3*size_3d)[:, i]) for i in range(3*size_3d)]).T
-K_jac = np.array([kern.jac_dot(np.eye(2*size_2d)[:, i]) for i in range(2*size_2d)]).T
-F_jac_mult = K.dot(P)
-F_jac_direct = np.array([F.jac_dot(None, np.eye(3*size_3d)[:, i]) for i in range(3*size_3d)]).T
-np.testing.assert_almost_equal(F_jac_direct, F_jac_mult)
-np.testing.assert_almost_equal(F_direct, F_jac_direct)
-P_t = np.array([data.projectors[0].jac_T_dot(np.eye(2*size_2d)[:, i]) for i in range(2*size_2d)]).T
-K_t = np.array([kern.jac_T_dot(np.eye(size_2d)[:, i]) for i in range(size_2d)]).T
-F_t_mult = P_t.dot(K_t)
-F_t_direct = np.array([F.jac_T_dot(None, np.eye(size_2d)[:, i]) for i in range(size_2d)]).T
-
-P_t_ref = P.T
-K_t_ref = K.T
-F_t_ref = F_mult.T
-np.testing.assert_almost_equal(P_t, P_t_ref)
-np.testing.assert_almost_equal(K_t, K_t_ref)
-np.testing.assert_almost_equal(F_t_mult, F_t_ref)
-np.testing.assert_almost_equal(F_t_direct, F_t_ref)
+logging.config.fileConfig(LOGGING_CONF, disable_existing_loggers=False)
###################################################################################################
-print('--Compliance test for two projections')
-
-a = 1.
-b_0 = 1000.
-dim = (2, 2, 2)
-count = 2
-
-magnitude = np.zeros((3,)+dim)
-magnitude[:, int(dim[0]/2), int(dim[1]/2), int(dim[2]/2)] = 1.
-
-mag_data = MagData(a, magnitude)
-
-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]
+print('--Generating input phase_maps')
+a = 10.
+b_0 = 1.
+dim = (16, 16, 16)
dim_uv = dim[1:3]
+count = 32
-data_collection = DataSet(a, dim_uv, b_0)
-
-[data_collection.append((phase_maps[i], projectors[i])) for i in range(count)]
-
-data = data_collection
-y = data.phase_vec
-kern = Kernel(data.a, data.dim_uv, data.b_0)
-F = ForwardModel(data.projectors, kern)
-C = Costfunction(y, F)
-
-size_2d = np.prod(dim[1]*dim[2])
-size_3d = np.prod(dim)
-
-P0 = np.array([data.projectors[0](np.eye(3*size_3d)[:, i]) for i in range(3*size_3d)]).T
-P1 = np.array([data.projectors[1](np.eye(3*size_3d)[:, i]) for i in range(3*size_3d)]).T
-P = np.vstack((P0, P1))
-K = np.array([kern(np.eye(2*size_2d)[:, i]) for i in range(2*size_2d)]).T
-F_mult0 = K.dot(P0)
-F_mult1 = K.dot(P1)
-F_mult = np.vstack((F_mult0, F_mult1))
-F_direct = np.array([F(np.eye(3*size_3d)[:, i]) for i in range(3*size_3d)]).T
-np.testing.assert_almost_equal(F_direct, F_mult)
-
-P_jac0 = np.array([data.projectors[0].jac_dot(np.eye(3*size_3d)[:, i]) for i in range(3*size_3d)]).T
-P_jac1 = np.array([data.projectors[1].jac_dot(np.eye(3*size_3d)[:, i]) for i in range(3*size_3d)]).T
-P = np.vstack((P0, P1))
-K_jac = np.array([kern.jac_dot(np.eye(2*size_2d)[:, i]) for i in range(2*size_2d)]).T
-F_jac_mult0 = K.dot(P_jac0)
-F_jac_mult1 = K.dot(P_jac1)
-F_jac_mult = np.vstack((F_jac_mult0, F_jac_mult1))
-F_jac_direct = np.array([F.jac_dot(None, np.eye(3*size_3d)[:, i]) for i in range(3*size_3d)]).T
-np.testing.assert_almost_equal(F_jac_direct, F_jac_mult)
-np.testing.assert_almost_equal(F_direct, F_jac_direct)
-
-P_t0 = np.array([data.projectors[0].jac_T_dot(np.eye(2*size_2d)[:, i]) for i in range(2*size_2d)]).T
-P_t1 = np.array([data.projectors[1].jac_T_dot(np.eye(2*size_2d)[:, i]) for i in range(2*size_2d)]).T
-P_t = np.hstack((P_t0, P_t1))
-K_t = np.array([kern.jac_T_dot(np.eye(size_2d)[:, i]) for i in range(size_2d)]).T
-F_t_mult0 = P_t0.dot(K_t)
-F_t_mult1 = P_t1.dot(K_t)
-F_t_mult = np.hstack((F_t_mult0, F_t_mult1))
-F_t_direct = np.array([F.jac_T_dot(None, np.eye(count*size_2d)[:, i]) for i in range(count*size_2d)]).T
-
-P_t_ref = P.T
-K_t_ref = K.T
-F_t_ref = F_mult.T
-np.testing.assert_almost_equal(P_t, P_t_ref)
-np.testing.assert_almost_equal(K_t, K_t_ref)
-np.testing.assert_almost_equal(F_t_mult, F_t_ref)
-np.testing.assert_almost_equal(F_t_direct, F_t_ref)
-
-
-
-
+mag_shape = mc.Shapes.disc(dim, (7.5, 7.5, 7.5), dim[2]/4, dim[2]/4)
+magnitude = mc.create_mag_dist_vortex(mag_shape)
+mag_data = MagData(a, magnitude)
+mag_data.quiver_plot3d()
+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)
+projectors_xy_full, projectors_x_full, projectors_y_full = [], [], []
+projectors_xy_miss, projectors_x_miss, projectors_y_miss = [], [], []
+projectors_x_full.extend([XTiltProjector(dim, tilt) for tilt in tilts_full])
+projectors_y_full.extend([YTiltProjector(dim, tilt) for tilt in tilts_full])
+projectors_xy_full.extend(projectors_x_full)
+projectors_xy_full.extend(projectors_y_full)
+projectors_x_miss.extend([XTiltProjector(dim, tilt) for tilt in tilts_miss])
+projectors_y_miss.extend([YTiltProjector(dim, tilt) for tilt in tilts_miss])
+projectors_xy_miss.extend(projectors_x_miss)
+projectors_xy_miss.extend(projectors_y_miss)
###################################################################################################
-print('--STARTING RECONSTRUCTION')
-
-
-
+projectors = projectors_xy_full
###################################################################################################
-print('--Generating input phase_maps')
-
-a = 10.
-b_0 = 1000.
-dim = (8, 8, 8)
-count = 32
-
-magnitude = np.zeros((3,)+dim)
-magnitude[:, 3:6, 3:6, 3:6] = 1# int(dim[0]/2), int(dim[1]/2), int(dim[2]/2)] = 1.
-magnitude = mc.create_mag_dist_vortex(mc.Shapes.disc(dim, (3.5, 3.5, 3.5), 3, 4))
-
-mag_data = MagData(a, magnitude)
-mag_data.quiver_plot3d()
-tilts = np.linspace(0, 2*pi, num=count/2, endpoint=False)
-projectors = []
-projectors.extend([XTiltProjector(mag_data.dim, tilt) for tilt in tilts])
-projectors.extend([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%(count/2)]/pi))
-# for i, phase_map in enumerate(phase_maps)]
-
###################################################################################################
print('--Setting up data collection')
dim_uv = dim[1:3]
-lam = 10. ** -10
-
-size_2d = np.prod(dim_uv)
-size_3d = np.prod(dim)
-
-
-data_collection = DataSet(a, dim_uv, b_0)
+lam = 10**-10
-[data_collection.append((phase_maps[i], projectors[i])) for i in range(count)]
-
-data = data_collection
+data = DataSet(a, dim_uv, b_0)
+[data.append((phase_maps[i], projectors[i])) for i in range(len(projectors))]
y = data.phase_vec
kern = Kernel(data.a, data.dim_uv, data.b_0)
F = ForwardModel(data.projectors, kern)
@@ -214,224 +81,7 @@ C = Costfunction(y, F, lam)
###################################################################################################
print('--Test simple solver')
-#M = np.asmatrix([F.jac_dot(None, np.eye(3*size_3d)[:, i]) for i in range(3*size_3d)]).T
-#MTM = M.T * M + lam * np.asmatrix(np.eye(3*size_3d))
-#A = MTM#np.array([F(np.eye(81)[:, i]) for i in range(81)])
-
-#class A_adapt(LinearOperator):
-#
-# def __init__(self, FwdModel, lam, shape):
-# self.fwd = FwdModel
-# self.lam = lam
-# self.shape = shape
-#
-# def matvec(self, vector):
-# return self.fwd.jac_T_dot(None, self.fwd.jac_dot(None, vector)) + self.lam*vector
-#
-# @property
-# def shape(self):
-# return self.shape
-#
-# @property
-# def dtype(self):
-# return np.dtype("d") # None #np.ones(1).dtype
-#
-## TODO: .shape in F und C
-#
-#b = F.jac_T_dot(None, y)
-#
-#A_fast = A_adapt(F, lam, (3*size_3d, 3*size_3d))
-#
-#i = 0
-#def printit(_):
-# global i
-# i += 1
-# print i
-#
-#x_f, info = cg(A_fast, b, callback=printit)
-#
-#mag_data_rec = MagData(a, np.zeros((3,)+dim))
-#
-#mag_data_rec.mag_vec = x_f
-#
-#mag_data_rec.quiver_plot3d()
-
-#phase_maps_rec = [pm(mag_data_rec) for pm in phasemappers]
-#[phase_map.display_phase(title=u'Tilt series (rec.) $(\phi = {:2.1f} \pi)$'.format(tilts[i%(count/2)]/pi))
-# for i, phase_map in enumerate(phase_maps_rec)]
-
-mag_data_opt = rc.optimize_sparse_cg(data_collection)
-
+mag_data_opt = rc.optimize_sparse_cg(data)
mag_data_opt.quiver_plot3d()
-
-
-#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('--Singular value decomposition')
-
-#a = 1.
-#b_0 = 1000.
-#dim = (3, 3, 3)
-#count = 8
-#
-#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/2, endpoint=False)
-#projectors = []
-#projectors.extend([XTiltProjector(mag_data.dim, tilt) for tilt in tilts])
-#projectors.extend([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%(count/2)]/pi))
-## for i, phase_map in enumerate(phase_maps)]
-#
-#dim_uv = dim[1:3]
-#
-#lam = 10. ** -10
-#
-#size_2d = np.prod(dim_uv)
-#size_3d = np.prod(dim)
-#
-#
-#data_collection = DataCollection(a, dim_uv, b_0)
-#
-#[data_collection.append((phase_maps[i], projectors[i])) for i in range(count)]
-#
-#data = data_collection
-#y = data.phase_vec
-#kern = Kernel(data.a, data.dim_uv, data.b_0)
-#F = ForwardModel(data.projectors, kern)
-#C = Costfunction(y, F, lam)
-#
-#mag_data_opt = opt.optimize_sparse_cg(data_collection)
-#mag_data_opt.quiver_plot3d()
-
-
-
-#M = np.asmatrix([F.jac_dot(None, np.eye(3*size_3d)[:, i]) for i in range(3*size_3d)]).T
-#MTM = M.T * M + lam * np.asmatrix(np.eye(3*size_3d))
-#A = MTM#np.array([F(np.eye(81)[:, i]) for i in range(81)])
-#
-#
-#
-#U, s, V = np.linalg.svd(M)
-##np.testing.assert_almost_equal(U.T, V, decimal=5)
-#
-#for value in range(20):
-# print 'Singular value:', s[value]
-# MagData(data.a, np.array(V[value,:]).reshape((3,)+dim)).quiver_plot3d()
-#
-#
-#for value in range(-10,0):
-# print 'Singular value:', s[value]
-# MagData(data.a, np.array(V[value,:]).reshape((3,)+dim)).quiver_plot3d()
-
-
-# TODO: print all singular vectors for a 2x2x2 distribution for each single and both tilt series!
-
-# TODO: Separate the script for SVD and compliance tests
-
-####################################################################################################
-#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
-
-
-
-
-
-
-
-
-
-
-
-#K = np.asmatrix([F.jac_dot(None, np.eye(24)[:, i]) for i in range(24)]).T
-#lam = 10. ** -10
-#KTK = K.T * K + lam * np.asmatrix(np.eye(24))
-#
-#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.dot(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()
-
-
-
-
+(mag_data_opt - mag_data).quiver_plot3d()
+#data.display(data.create_phase_maps(mag_data_opt))
diff --git a/scripts/test methods/compare_pixel_fields.py b/scripts/test methods/compare_pixel_fields.py
index d97cea2..9fd7bdf 100644
--- a/scripts/test methods/compare_pixel_fields.py
+++ b/scripts/test methods/compare_pixel_fields.py
@@ -32,9 +32,10 @@ if not os.path.exists(directory):
a = 1.0 # in nm
phi = 0 # in rad
dim = (1, 32, 32)
-pixel = (0, int(dim[1]/2), int(dim[2]/2))
+pixel = (0, int(dim[1]/2), int(dim[2]/2))
limit = 0.35
+
def get_fourier_kernel():
PHI_0 = 2067.83 # magnetic flux in T*nm²
b_0 = 1
@@ -43,12 +44,13 @@ def get_fourier_kernel():
f_u = np.linspace(0, nyq/2, dim[2]/2.+1)
f_v = np.linspace(-nyq/2, nyq/2, dim[1], endpoint=False)
f_uu, f_vv = np.meshgrid(f_u, f_v)
- phase_fft = coeff * f_vv / (f_uu**2 + f_vv**2 + 1e-30) #* (8*(a/2)**3)*np.sinc(a/2*f_uu)*np.sinc(a/2*f_vv)
+ phase_fft = coeff * f_vv / (f_uu**2 + f_vv**2 + 1e-30)
# Transform to real space and revert padding:
phase_fft = np.fft.ifftshift(phase_fft, axes=0)
phase_fft_kernel = np.fft.fftshift(np.fft.irfft2(phase_fft), axes=(0, 1))
return phase_fft_kernel
+
# Create magnetic data and projector:
mag_data = MagData(a, mc.create_mag_dist_homog(mc.Shapes.pixel(dim, pixel), phi))
mag_data.save_to_llg(directory + '/mag_dist_single_pixel.txt')
diff --git a/scripts/test methods/kernel_comparison.py b/scripts/test methods/kernel_comparison.py
deleted file mode 100644
index 23cc8e5..0000000
--- a/scripts/test methods/kernel_comparison.py
+++ /dev/null
@@ -1,23 +0,0 @@
-# -*- coding: utf-8 -*-
-"""
-Created on Thu Sep 12 13:22:43 2013
-
-@author: Jan
-"""
-
-from pylab import *
-
-
-f = 1
-a, b = f * 32, f * 64
-f_u = np.linspace(0, 1./3, a)
-f_v = np.linspace(-1./6., 1./6., b, endpoint=False)
-f_uu, f_vv = np.meshgrid(f_u, f_v)
-phase_fft = 1j * f_vv / (f_uu**2 + f_vv**2 + 1e-30)
-# Transform to real space and revert padding:
-phase_fft = np.fft.ifftshift(phase_fft, axes=0)
-ss = np.fft.irfft2(phase_fft)
-print ss
-pcolormesh(np.fft.fftshift(np.fft.fftshift(ss, axes=1), axes=0), cmap='RdBu')
-colorbar()
-show()
diff --git a/scripts/test methods/test_arb_projection.py b/scripts/test methods/test_arb_projection.py
index c652331..a055d33 100644
--- a/scripts/test methods/test_arb_projection.py
+++ b/scripts/test methods/test_arb_projection.py
@@ -1,17 +1,26 @@
# -*- coding: utf-8 -*-
-"""
-Created on Tue Sep 03 12:55:40 2013
+"""Created on Tue Sep 03 12:55:40 2013 @author: Jan"""
-@author: Jan
-"""
+import os
import numpy as np
from numpy import pi
-import itertools
import matplotlib.pyplot as plt
from matplotlib.ticker import MaxNLocator, NullLocator
+import pyramid
+
+import itertools
+
+import logging
+import logging.config
+
+
+LOGGING_CONF = os.path.join(os.path.dirname(os.path.realpath(pyramid.__file__)), 'logging.ini')
+
+
+logging.config.fileConfig(LOGGING_CONF, disable_existing_loggers=False)
dim = (8, 8)
offset = (dim[0]/2., dim[1]/2.)
diff --git a/scripts/test methods/test_new_structure.py b/scripts/test methods/test_new_structure.py
deleted file mode 100644
index c7afa84..0000000
--- a/scripts/test methods/test_new_structure.py
+++ /dev/null
@@ -1,31 +0,0 @@
-# -*- coding: utf-8 -*-
-"""
-Created on Thu Jan 09 20:18:56 2014
-
-@author: Jan
-"""
-
-
-from numpy import pi
-
-import pyramid.magcreator as mc
-from pyramid.magdata import MagData
-from pyramid.projector import SimpleProjector
-from pyramid.phasemapper import PMAdapterFM
-from pyramid.kernel import Kernel
-
-
-magnitude = mc.create_mag_dist_vortex(mc.Shapes.disc((32,64,128),(32,32,64),16,8))
-mag_data = MagData(1, magnitude)
-#mag_data = MagData.load_from_llg('../output/magnetic distributions/mag_dist_disc.txt')
-#mag_data.quiver_plot()
-phase_map = PMAdapterFM(mag_data.a, SimpleProjector(mag_data.dim), geometry='disc')(mag_data)
-phase_map.display_combined(30)
-
-#phase_map.make_color_wheel()
-
-# TODO: use for compare pixel fields!
-k_d = Kernel(1, (5,5), geometry='disc')
-k_s = Kernel(1, (5,5), geometry='slab')
-u_d, v_d = k_d.u, k_d.v
-u_s, v_s = k_s.u, k_s.v
diff --git a/tests/test_compliance.py b/tests/test_compliance.py
index 05f707a..99e6f06 100644
--- a/tests/test_compliance.py
+++ b/tests/test_compliance.py
@@ -26,8 +26,8 @@ class TestCaseCompliance(unittest.TestCase):
for root, dirs, files in os.walk(rootdir):
for filename in files:
if ((filename.endswith('.py') or filename.endswith('.pyx'))
- and root != os.path.join(self.path, 'scripts', 'gui')):
- filepaths.append(os.path.join(root, filename))
+ and root != os.path.join(self.path, 'scripts', 'gui')):
+ filepaths.append(os.path.join(root, filename))
return filepaths
def test_pep8(self):
@@ -35,7 +35,7 @@ class TestCaseCompliance(unittest.TestCase):
files = self.get_files_to_check(os.path.join(self.path, 'pyramid')) \
+ self.get_files_to_check(os.path.join(self.path, 'scripts')) \
+ self.get_files_to_check(os.path.join(self.path, 'tests'))
- ignores = ('E226', 'E128')
+ ignores = ('E125', 'E226', 'E228')
pep8.MAX_LINE_LENGTH = 99
pep8style = pep8.StyleGuide(quiet=False)
pep8style.options.ignore = ignores
@@ -57,10 +57,10 @@ class TestCaseCompliance(unittest.TestCase):
todo_count = 0
regex = ur'# TODO: (.*)'
for py_file in files:
- with open (py_file) as f:
+ with open(py_file) as f:
for line in f:
todo = re.findall(regex, line)
- if todo and not todo[0]=="(.*)'":
+ if todo and not todo[0] == "(.*)'":
todos_found = True
todo_count += 1
print '{}: {}'.format(f.name, todo[0])
diff --git a/tests/test_costfunction.py b/tests/test_costfunction.py
index 0d06db5..7f7260e 100644
--- a/tests/test_costfunction.py
+++ b/tests/test_costfunction.py
@@ -1,7 +1,2 @@
# -*- coding: utf-8 -*-
-"""
-Created on Mon Jan 20 20:01:25 2014
-
-@author: Jan
-"""
-
+"""Created on Mon Jan 20 20:01:25 2014 @author: Jan"""
diff --git a/tests/test_dataset.py b/tests/test_dataset.py
index ebe2322..6b97d84 100644
--- a/tests/test_dataset.py
+++ b/tests/test_dataset.py
@@ -1,7 +1,2 @@
# -*- coding: utf-8 -*-
-"""
-Created on Mon Jan 20 19:58:01 2014
-
-@author: Jan
-"""
-
+"""Created on Mon Jan 20 19:58:01 2014 @author: Jan"""
diff --git a/tests/test_forwardmodel.py b/tests/test_forwardmodel.py
index 5bd1ef1..acf99af 100644
--- a/tests/test_forwardmodel.py
+++ b/tests/test_forwardmodel.py
@@ -1,7 +1,2 @@
# -*- coding: utf-8 -*-
-"""
-Created on Mon Jan 20 20:01:13 2014
-
-@author: Jan
-"""
-
+"""Created on Mon Jan 20 20:01:13 2014 @author: Jan"""
diff --git a/tests/test_kernel.py b/tests/test_kernel.py
index b4fa9e4..39cb1e1 100644
--- a/tests/test_kernel.py
+++ b/tests/test_kernel.py
@@ -1,7 +1,2 @@
# -*- coding: utf-8 -*-
-"""
-Created on Mon Jan 20 20:00:50 2014
-
-@author: Jan
-"""
-
+"""Created on Mon Jan 20 20:00:50 2014 @author: Jan"""
diff --git a/tests/test_reconstruction.py b/tests/test_reconstruction.py
index 7cb6fd3..2bd2343 100644
--- a/tests/test_reconstruction.py
+++ b/tests/test_reconstruction.py
@@ -5,8 +5,6 @@
import os
import unittest
-import pyramid.reconstruction as rc
-
class TestCaseReconstruction(unittest.TestCase):
@@ -23,3 +21,141 @@ class TestCaseReconstruction(unittest.TestCase):
if __name__ == '__main__':
suite = unittest.TestLoader().loadTestsFromTestCase(TestCaseReconstruction)
unittest.TextTestRunner(verbosity=2).run(suite)
+
+
+###################################################################################################
+#print('--Compliance test for one projection')
+#
+#a = 1.
+#b_0 = 1000.
+#dim = (2, 2, 2)
+#count = 1
+#
+#magnitude = np.zeros((3,)+dim)
+#magnitude[:, int(dim[0]/2), int(dim[1]/2), int(dim[2]/2)] = 1.
+#
+#mag_data = MagData(a, magnitude)
+#
+#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]
+#
+#dim_uv = dim[1:3]
+#
+#data_collection = DataSet(a, dim_uv, b_0)
+#
+#[data_collection.append((phase_maps[i], projectors[i])) for i in range(count)]
+#
+#data = data_collection
+#y = data.phase_vec
+#kern = Kernel(data.a, data.dim_uv, data.b_0)
+#F = ForwardModel(data.projectors, kern)
+#C = Costfunction(y, F)
+#
+#size_2d = np.prod(dim[1]*dim[2])
+#size_3d = np.prod(dim)
+#
+#P = np.array([data.projectors[0](np.eye(3*size_3d)[:, i]) for i in range(3*size_3d)]).T
+#K = np.array([kern(np.eye(2*size_2d)[:, i]) for i in range(2*size_2d)]).T
+#F_mult = K.dot(P)
+#F_direct = np.array([F(np.eye(3*size_3d)[:, i]) for i in range(3*size_3d)]).T
+#np.testing.assert_almost_equal(F_direct, F_mult)
+#
+#P_jac = np.array([data.projectors[0].jac_dot(np.eye(3*size_3d)[:, i])
+# for i in range(3*size_3d)]).T
+#K_jac = np.array([kern.jac_dot(np.eye(2*size_2d)[:, i]) for i in range(2*size_2d)]).T
+#F_jac_mult = K.dot(P)
+#F_jac_direct = np.array([F.jac_dot(None, np.eye(3*size_3d)[:, i]) for i in range(3*size_3d)]).T
+#np.testing.assert_almost_equal(F_jac_direct, F_jac_mult)
+#np.testing.assert_almost_equal(F_direct, F_jac_direct)
+#
+#P_t = np.array([data.projectors[0].jac_T_dot(np.eye(2*size_2d)[:, i])
+# for i in range(2*size_2d)]).T
+#K_t = np.array([kern.jac_T_dot(np.eye(size_2d)[:, i]) for i in range(size_2d)]).T
+#F_t_mult = P_t.dot(K_t)
+#F_t_direct = np.array([F.jac_T_dot(None, np.eye(size_2d)[:, i]) for i in range(size_2d)]).T
+#
+#P_t_ref = P.T
+#K_t_ref = K.T
+#F_t_ref = F_mult.T
+#np.testing.assert_almost_equal(P_t, P_t_ref)
+#np.testing.assert_almost_equal(K_t, K_t_ref)
+#np.testing.assert_almost_equal(F_t_mult, F_t_ref)
+#np.testing.assert_almost_equal(F_t_direct, F_t_ref)
+#
+###################################################################################################
+#print('--Compliance test for two projections')
+#
+#a = 1.
+#b_0 = 1000.
+#dim = (2, 2, 2)
+#count = 2
+#
+#magnitude = np.zeros((3,)+dim)
+#magnitude[:, int(dim[0]/2), int(dim[1]/2), int(dim[2]/2)] = 1.
+#
+#mag_data = MagData(a, magnitude)
+#
+#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]
+#
+#dim_uv = dim[1:3]
+#
+#data_collection = DataSet(a, dim_uv, b_0)
+#
+#[data_collection.append((phase_maps[i], projectors[i])) for i in range(count)]
+#
+#data = data_collection
+#y = data.phase_vec
+#kern = Kernel(data.a, data.dim_uv, data.b_0)
+#F = ForwardModel(data.projectors, kern)
+#C = Costfunction(y, F)
+#
+#size_2d = np.prod(dim[1]*dim[2])
+#size_3d = np.prod(dim)
+#
+#P0 = np.array([data.projectors[0](np.eye(3*size_3d)[:, i]) for i in range(3*size_3d)]).T
+#P1 = np.array([data.projectors[1](np.eye(3*size_3d)[:, i]) for i in range(3*size_3d)]).T
+#P = np.vstack((P0, P1))
+#K = np.array([kern(np.eye(2*size_2d)[:, i]) for i in range(2*size_2d)]).T
+#F_mult0 = K.dot(P0)
+#F_mult1 = K.dot(P1)
+#F_mult = np.vstack((F_mult0, F_mult1))
+#F_direct = np.array([F(np.eye(3*size_3d)[:, i]) for i in range(3*size_3d)]).T
+#np.testing.assert_almost_equal(F_direct, F_mult)
+#
+#P_jac0 = np.array([data.projectors[0].jac_dot(np.eye(3*size_3d)[:, i])
+# for i in range(3*size_3d)]).T
+#P_jac1 = np.array([data.projectors[1].jac_dot(np.eye(3*size_3d)[:, i])
+# for i in range(3*size_3d)]).T
+#P = np.vstack((P0, P1))
+#K_jac = np.array([kern.jac_dot(np.eye(2*size_2d)[:, i]) for i in range(2*size_2d)]).T
+#F_jac_mult0 = K.dot(P_jac0)
+#F_jac_mult1 = K.dot(P_jac1)
+#F_jac_mult = np.vstack((F_jac_mult0, F_jac_mult1))
+#F_jac_direct = np.array([F.jac_dot(None, np.eye(3*size_3d)[:, i]) for i in range(3*size_3d)]).T
+#np.testing.assert_almost_equal(F_jac_direct, F_jac_mult)
+#np.testing.assert_almost_equal(F_direct, F_jac_direct)
+#
+#P_t0 = np.array([data.projectors[0].jac_T_dot(np.eye(2*size_2d)[:, i])
+# for i in range(2*size_2d)]).T
+#P_t1 = np.array([data.projectors[1].jac_T_dot(np.eye(2*size_2d)[:, i])
+# for i in range(2*size_2d)]).T
+#P_t = np.hstack((P_t0, P_t1))
+#K_t = np.array([kern.jac_T_dot(np.eye(size_2d)[:, i]) for i in range(size_2d)]).T
+#F_t_mult0 = P_t0.dot(K_t)
+#F_t_mult1 = P_t1.dot(K_t)
+#F_t_mult = np.hstack((F_t_mult0, F_t_mult1))
+#F_t_direct = np.array([F.jac_T_dot(None, np.eye(count*size_2d)[:, i])
+# for i in range(count*size_2d)]).T
+#
+#P_t_ref = P.T
+#K_t_ref = K.T
+#F_t_ref = F_mult.T
+#np.testing.assert_almost_equal(P_t, P_t_ref)
+#np.testing.assert_almost_equal(K_t, K_t_ref)
+#np.testing.assert_almost_equal(F_t_mult, F_t_ref)
+#np.testing.assert_almost_equal(F_t_direct, F_t_ref)
--
GitLab