From 04d9ed45ee96b6dbe1ddc2ee958b26dadc4908a7 Mon Sep 17 00:00:00 2001
From: Jan Caron <j.caron@fz-juelich.de>
Date: Tue, 28 Oct 2014 22:45:00 +0100
Subject: [PATCH] Script fixes and minor corrections to the refactoring. Still
to do: PhaseMapperElectric does not work 3D reconstruction does
not converge
---
pyramid/costfunction.py | 2 +-
pyramid/dataset.py | 4 +-
pyramid/forwardmodel.py | 9 +-
pyramid/magdata.py | 6 +-
pyramid/phasemap.py | 2 +-
pyramid/phasemapper.py | 3 +-
pyramid/reconstruction.py | 13 +-
pyramid/regularisator.py | 106 ++++---
scripts/collaborations/bielefeld.py | 1 +
scripts/collaborations/example_joern.py | 14 +-
scripts/collaborations/rueffner_file.py | 15 +-
scripts/collaborations/vadim_integration.py | 33 +--
.../create_alternating_filaments.py | 2 +
.../create_core_shell_disc.py | 11 +-
.../create_core_shell_sphere.py | 11 +-
.../create_multiple_samples.py | 7 +-
.../create_random_pixels.py | 7 +-
.../create_random_slabs.py | 7 +-
scripts/create distributions/create_vortex.py | 2 +-
scripts/gui/mag_slicer.py | 6 +-
scripts/publications/abstract_prague.py | 23 +-
.../paper 1/ch5-1-evaluation_real_space.py | 16 +-
.../paper 1/ch5-2-evaluation_fourier_space.py | 34 ++-
.../paper 1/ch5-3-comparison_of_methods.py | 263 ------------------
.../ch5-4-comparison_of_methods_new.py | 37 +--
scripts/publications/poster_pof.py | 13 +-
.../reconstruct_random_pixels.py | 11 +-
...econstruction_sparse_cg_singular_values.py | 7 +-
.../reconstruction_sparse_cg_test.py | 53 ++--
scripts/test methods/compare_methods.py | 42 ++-
scripts/test methods/compare_pixel_fields.py | 9 +-
scripts/test methods/projector_test.py | 16 +-
32 files changed, 269 insertions(+), 516 deletions(-)
delete mode 100644 scripts/publications/paper 1/ch5-3-comparison_of_methods.py
diff --git a/pyramid/costfunction.py b/pyramid/costfunction.py
index 547ad97..04ed291 100644
--- a/pyramid/costfunction.py
+++ b/pyramid/costfunction.py
@@ -151,7 +151,7 @@ class CFAdapterScipyCG(LinearOperator):
@property
def shape(self):
- return (3*self.cost.fwd_model.size_3d, 3*self.cost.fwd_model.size_3d)
+ return (self.cost.data_set.m, self.cost.data_set.m)
@property
def dtype(self):
diff --git a/pyramid/dataset.py b/pyramid/dataset.py
index e8c36bf..fd5fb49 100644
--- a/pyramid/dataset.py
+++ b/pyramid/dataset.py
@@ -76,8 +76,8 @@ class DataSet(object):
@property
def phase_mappers(self):
- dim_uv_list = np.unique([p.dim_uv for p in self.phase_maps])
- kernel_list = [Kernel(self.a, dim_uv) for dim_uv in dim_uv_list]
+ dim_uv_list = np.unique([p.dim_uv for p in self.projectors])
+ kernel_list = [Kernel(self.a, tuple(dim_uv)) for dim_uv in dim_uv_list]
return {kernel.dim_uv: PhaseMapperRDFC(kernel) for kernel in kernel_list}
def __init__(self, a, dim, b_0=1, mask=None):
diff --git a/pyramid/forwardmodel.py b/pyramid/forwardmodel.py
index c68e93f..c48a7f7 100644
--- a/pyramid/forwardmodel.py
+++ b/pyramid/forwardmodel.py
@@ -5,8 +5,6 @@ threedimensional magnetization distribution onto a two-dimensional phase map."""
import numpy as np
-from pyramid.projector import Projector
-from pyramid.phasemapper import PhaseMapper
from pyramid.magdata import MagData
import logging
@@ -62,8 +60,7 @@ class ForwardModel(object):
def __call__(self, x):
self.LOG.debug('Calling __call__')
self.mag_data.set_vector(x, self.data_set.mask)
- return np.reshape(self.data_set.create_phase_maps(self.mag_data), -1)
- # TODO: result mit fixer Größe und dann Multiprocessing, dataset needs list of start points
+ # TODO: Multiprocessing
result = np.zeros(self.n)
hp = self.hook_points
for i, projector in enumerate(self.data_set.projectors):
@@ -92,13 +89,13 @@ class ForwardModel(object):
'''
self.LOG.debug('Calling jac_dot')
- self.mag_data.set_vector(x, self.data_set.mask)
+ self.mag_data.set_vector(vector, self.data_set.mask)
result = np.zeros(self.n)
hp = self.hook_points
for i, projector in enumerate(self.data_set.projectors):
mag_vec = self.mag_data.mag_vec
res = self.phase_mappers[projector.dim_uv].jac_dot(projector.jac_dot(mag_vec))
- result[hp[i]:hp[i+1]] = res
+ result[hp[i]:hp[i+1]] = res.flatten()
return result
def jac_T_dot(self, x, vector):
diff --git a/pyramid/magdata.py b/pyramid/magdata.py
index 8b22db9..5989153 100644
--- a/pyramid/magdata.py
+++ b/pyramid/magdata.py
@@ -281,9 +281,9 @@ class MagData(object):
Order is: first all `x`-, then all `y`-, then all `z`-components.
'''
- return np.concatenate([self.magnitude[2][mask],
- self.magnitude[1][mask],
- self.magnitude[0][mask]])
+ return np.reshape([self.magnitude[2][mask],
+ self.magnitude[1][mask],
+ self.magnitude[0][mask]], -1)
def set_vector(self, vector, mask=None):
'''Set the magnetic components of the masked pixels to the values specified by `vector`.
diff --git a/pyramid/phasemap.py b/pyramid/phasemap.py
index ed46cdf..af1e708 100644
--- a/pyramid/phasemap.py
+++ b/pyramid/phasemap.py
@@ -437,7 +437,7 @@ class PhaseMap(object):
self.LOG.debug('Calling display_holo')
# Calculate gain if 'auto' is selected:
if gain == 'auto':
- gain = 5 * 2*pi/self.phase.max()
+ gain = 4 * 2*pi/self.phase.max()
# Set title if not set:
if title is None:
title = 'Contour Map (gain: %.2g)' % gain
diff --git a/pyramid/phasemapper.py b/pyramid/phasemapper.py
index ec54213..1d4205a 100644
--- a/pyramid/phasemapper.py
+++ b/pyramid/phasemapper.py
@@ -244,7 +244,7 @@ class PhaseMapperRDRC(PhaseMapper):
phase += u_mag[j, i] * u_phase
v_phase = v_phi[dim_uv[0]-1-j:(2*dim_uv[0]-1)-j,
dim_uv[1]-1-i:(2*dim_uv[1]-1)-i]
- if abs(v_mag[j, i]) > self.hreshold:
+ if abs(v_mag[j, i]) > self.threshold:
phase -= v_mag[j, i] * v_phase
# Return the phase:
return PhaseMap(mag_data.a, phase)
@@ -494,6 +494,7 @@ class PhaseMapperElectric(PhaseMapper):
assert mag_data.a == self.a, 'Grid spacing has to match!'
assert mag_data.dim[0] == 1, 'Magnetic distribution must be 2-dimensional!'
assert mag_data.dim[1:3] == self.dim_uv, 'Dimensions do not match!'
+ return self.coeff * mag_data.get_mask(self.threshold)[0, ...].reshape(self.dim_uv)
# Calculate mask:
mask = mag_data.get_mask(self.threshold)
# Project and calculate phase:
diff --git a/pyramid/reconstruction.py b/pyramid/reconstruction.py
index 721d797..90704d5 100644
--- a/pyramid/reconstruction.py
+++ b/pyramid/reconstruction.py
@@ -83,7 +83,7 @@ class PrintIterator(object):
return self.iteration
-def optimize_sparse_cg(data, Se_inv=None, regularisator=None, verbosity=0):
+def optimize_sparse_cg(data, regularisator=None, maxiter=1000, verbosity=0):
'''Reconstruct a three-dimensional magnetic distribution from given phase maps via the
conjugate gradient optimizaion method :func:`~.scipy.sparse.linalg.cg`.
@@ -100,6 +100,8 @@ def optimize_sparse_cg(data, Se_inv=None, regularisator=None, verbosity=0):
regularisator : :class:`~.Regularisator`, optional
Regularisator class that's responsible for the regularisation term. Defaults to zero
order Tikhonov if none is provided.
+ maxiter : int
+ Maximum number of iterations.
verbosity : {0, 1, 2}, optional
Parameter defining the verposity of the output. `2` will show the current number of the
iteration and the cost of the current distribution. `1` will just show the iteration
@@ -114,15 +116,14 @@ def optimize_sparse_cg(data, Se_inv=None, regularisator=None, verbosity=0):
LOG.debug('Calling optimize_sparse_cg')
# Set up necessary objects:
y = data.phase_vec
- kernel = Kernel(data.a, data.dim_uv, data.b_0)
- fwd_model = ForwardModel(data.projectors, kernel)
- cost = Costfunction(y, fwd_model, Se_inv, regularisator)
+ fwd_model = ForwardModel(data)
+ cost = Costfunction(data, regularisator)
# Optimize:
A = CFAdapterScipyCG(cost)
b = fwd_model.jac_T_dot(None, y)
- x_opt, info = cg(A, b, callback=PrintIterator(cost, verbosity), maxiter=1000)
+ x_opt, info = cg(A, b, callback=PrintIterator(cost, verbosity), maxiter=maxiter)
# Create and return fitting MagData object:
- mag_opt = MagData(fwd_model.a, np.zeros((3,)+fwd_model.dim))
+ mag_opt = MagData(data.a, np.zeros((3,)+data.dim))
mag_opt.mag_vec = x_opt
return mag_opt
diff --git a/pyramid/regularisator.py b/pyramid/regularisator.py
index ee8ee99..a7adf41 100644
--- a/pyramid/regularisator.py
+++ b/pyramid/regularisator.py
@@ -10,8 +10,8 @@ import abc
import numpy as np
-from scipy.sparse import eye, coo_matrix, csr_matrix
-import jutil.norm as jnorm
+from scipy.sparse import coo_matrix, csr_matrix
+import jutil.norms as jnorm
from pyramid.converter import IndexConverter
@@ -51,67 +51,87 @@ class Regularisator(object):
def jac(self, x):
# TODO: Docstring!
+ self.LOG.debug('Calling jac')
return self.lam * self.norm.jac_dot(x)
def hess_dot(self, x, vector):
# TODO: Docstring!
- return self.lam * self.norm.hess_dot(x)
+ self.LOG.debug('Calling hess_dot')
+ return self.lam * self.norm.hess_dot(x, vector)
- # TODO: hess_diag
+ def hess_diag(self, x, vector):
+ # TODO: Docstring!
+ self.LOG.debug('Calling hess_diag')
+ return self.lam * self.norm.hess_diag(x, vector)
-class ZeroOrderRegularisator(Regularisator):
- # TODO: Docstring!
+class NoneRegularisator(Regularisator):
+ # TODO: Docstring
- LOG = logging.getLogger(__name__+'.ZeroOrderRegularisator')
+ # TODO: Necessary class? Use others with lam=0?
- def __init__(self, lam):
+ LOG = logging.getLogger(__name__+'.NoneRegularisator')
+
+ def __init__(self):
self.LOG.debug('Calling __init__')
- norm = jnorm.NormL2Square()
- super(ZeroOrderRegularisator, self).__init__(norm, lam)
+ self.norm = None
+ self.lam = 0
self.LOG.debug('Created '+str(self))
def __call__(self, x):
self.LOG.debug('Calling __call__')
- super(ZeroOrderRegularisator, self)(x)
-
-
-
-class FirstOrderRegularisator(Regularisator):
- # TODO: Docstring!
-
- def __init__(self, fwd_model, lam, x_a=None):
- size_3d = fwd_model.size_3d
- dim = fwd_model.dim
- converter = IndexConverter(dim)
- row = []
- col = []
- data = []
-
- for i in range(size_3d):
- neighbours = converter.find_neighbour_ind(i)
-
-
-
-
-
-
-
+ return 0
+ def jac(self, x):
+ # TODO: Docstring!
+ self.LOG.debug('Calling jac')
+ return np.zeros_like(x)
+ def hess_dot(self, x, vector):
+ # TODO: Docstring!
+ self.LOG.debug('Calling hess_dot')
+ return np.zeros_like(vector)
- Sa_inv = csr_matrix(coo_matrix(data, (rows, columns)), shape=(3*size_3d, 3*size_3d))
+ def hess_diag(self, x, vector):
+ # TODO: Docstring!
+ self.LOG.debug('Calling hess_diag')
+ return np.zeros_like(vector)
+class ZeroOrderRegularisator(Regularisator):
+ # TODO: Docstring!
+ LOG = logging.getLogger(__name__+'.ZeroOrderRegularisator')
- term2 = []
- for i in range(3):
- component = mag_data[i, ...]
- for j in range(3):
- if component.shape[j] > 1:
- term2.append(np.diff(component, axis=j).reshape(-1))
+ def __init__(self, lam):
+ self.LOG.debug('Calling __init__')
+ norm = jnorm.NormL2Square()
+ super(ZeroOrderRegularisator, self).__init__(norm, lam)
+ self.LOG.debug('Created '+str(self))
- super(FirstOrderRegularisator, self).__init__(Sa_inv, x_a)
- self.LOG.debug('Created '+str(self))
+#class FirstOrderRegularisator(Regularisator):
+# # TODO: Docstring!
+#
+# def __init__(self, fwd_model, lam, x_a=None):
+# size_3d = fwd_model.size_3d
+# dim = fwd_model.dim
+# converter = IndexConverter(dim)
+# row = []
+# col = []
+# data = []
+#
+# for i in range(size_3d):
+# neighbours = converter.find_neighbour_ind(i)
+#
+# Sa_inv = csr_matrix(coo_matrix(data, (rows, columns)), shape=(3*size_3d, 3*size_3d))
+#
+# term2 = []
+# for i in range(3):
+# component = mag_data[i, ...]
+# for j in range(3):
+# if component.shape[j] > 1:
+# term2.append(np.diff(component, axis=j).reshape(-1))
+#
+# super(FirstOrderRegularisator, self).__init__(Sa_inv, x_a)
+# self.LOG.debug('Created '+str(self))
diff --git a/scripts/collaborations/bielefeld.py b/scripts/collaborations/bielefeld.py
index f40fde2..cc2b5bb 100644
--- a/scripts/collaborations/bielefeld.py
+++ b/scripts/collaborations/bielefeld.py
@@ -29,6 +29,7 @@ def load_from_llg(filename):
# import pdb; pdb.set_trace()
return MagData(a, magnitude)
+
logging.config.fileConfig(LOGGING_CONF, disable_existing_loggers=False)
## 2DCoFe:
mag_data_2DCoFe = MagData.load_from_llg(PATH+'magnetic_distribution_2DCoFe.txt')
diff --git a/scripts/collaborations/example_joern.py b/scripts/collaborations/example_joern.py
index aa9aa60..7890c35 100644
--- a/scripts/collaborations/example_joern.py
+++ b/scripts/collaborations/example_joern.py
@@ -13,7 +13,9 @@ import matplotlib.pyplot as plt
import pyramid
from pyramid.phasemap import PhaseMap
from pyramid.projector import SimpleProjector
-from pyramid.phasemapper import PMConvolve
+from pyramid.phasemapper import pm
+from pyramid.dataset import DataSet
+from pyramid.regularisator import ZeroOrderRegularisator
import pyramid.reconstruction as rc
from time import clock
@@ -32,6 +34,7 @@ a = 1.0 # in nm
gain = 5
b_0 = 1
inter = 'none'
+dim = (1,) + (64, 64)
dim_small = (64, 64)
smoothed_pictures = True
lam = 1E-4
@@ -44,12 +47,17 @@ phase_map = PhaseMap.load_from_netcdf4(PATH+'phase_map.nc')
#mask = np.genfromtxt('mask.txt', dtype=bool)
with open(PATH+'mask.pickle') as pf:
mask = pickle.load(pf)
+# Setup:
+data_set = DataSet(a, dim, b_0)
+data_set.append(phase_map, SimpleProjector(dim))
+regularisator = ZeroOrderRegularisator(lam)
# Reconstruct the magnetic distribution:
tic = clock()
-mag_data_rec = rc.optimize_simple_leastsq(phase_map, mask, b_0, lam=lam, order=order)
+mag_data_rec = rc.optimize_sparse_cg(data_set, regularisator, maxiter=100, verbosity=2)
+# .optimize_simple_leastsq(phase_map, mask, b_0, lam=lam, order=order)
print 'reconstruction time:', clock() - tic
# Display the reconstructed phase map and holography image:
-phase_map_rec = PMConvolve(a, SimpleProjector(mag_data_rec.dim), b_0)(mag_data_rec)
+phase_map_rec = pm(mag_data_rec)
phase_map_rec.display_combined('Reconstr. Distribution', gain=gain, interpolation=inter)
# Plot the magnetization:
axis = (mag_data_rec*(1/mag_data_rec.magnitude.max())).quiver_plot()
diff --git a/scripts/collaborations/rueffner_file.py b/scripts/collaborations/rueffner_file.py
index ef0c2ac..230c4cb 100644
--- a/scripts/collaborations/rueffner_file.py
+++ b/scripts/collaborations/rueffner_file.py
@@ -14,7 +14,8 @@ import tables
import pyramid
from pyramid.magdata import MagData
from pyramid.projector import SimpleProjector
-from pyramid.phasemapper import PMConvolve
+from pyramid.phasemapper import PhaseMapperRDFC
+from pyramid.kernel import Kernel
import logging
import logging.config
@@ -38,8 +39,8 @@ lim_x = 0.35
# Build projectors and phasemapper:
projector_z = SimpleProjector(dim, axis='z', dim_uv=dim_uv)
projector_x = SimpleProjector(dim, axis='x', dim_uv=dim_uv)
-pm_z = PMConvolve(a, projector_z, b_0)
-pm_x = PMConvolve(a, projector_x, b_0)
+pm_z = PhaseMapperRDFC(Kernel(a, projector_z.dim_uv, b_0))
+pm_x = PhaseMapperRDFC(Kernel(a, projector_x.dim_uv, b_0))
# Read in hdf5-file and extract points:
h5file = tables.openFile(PATH+'dyn_0090mT_dyn.h5_dat.h5')
points = h5file.root.mesh.points.read()
@@ -85,13 +86,13 @@ def calculate(t): # TODO: Somehow avoid memory error :-(...
magnitude[j, i, :, :] = grid.filled(fill_value=0)
# Create magnetic distribution and phase maps:
mag_data = MagData(a, magnitude)
- phase_map_z = pm_z(mag_data)
- phase_map_x = pm_x(mag_data)
+ phase_map_z = pm_z(projector_z(mag_data))
+ phase_map_x = pm_x(projector_x(mag_data))
phase_map_z.unit = 'mrad'
# Display phase maps and save them to png:
- phase_map_z.display_combined(density=dens_z, interpolation='bilinear', limit=lim_z)
+ phase_map_z.display_combined(gain=dens_z, interpolation='bilinear', limit=lim_z)
plt.savefig(PATH+'rueffner/phase_map_z_t_'+str(t)+'.png')
- phase_map_x.display_combined(density=dens_x, interpolation='bilinear', limit=lim_x)
+ phase_map_x.display_combined(gain=dens_x, interpolation='bilinear', limit=lim_x)
plt.savefig(PATH+'rueffner/phase_map_x_t_'+str(t)+'.png')
vectors, data, magnitude, z_slice, weights, grid_x, grid_y, grid_z, grid, mag_data, \
phase_map_z, phase_map_x = [None] * 12
diff --git a/scripts/collaborations/vadim_integration.py b/scripts/collaborations/vadim_integration.py
index a1c9956..ce5b76d 100644
--- a/scripts/collaborations/vadim_integration.py
+++ b/scripts/collaborations/vadim_integration.py
@@ -10,7 +10,8 @@ from numpy import pi
import pyramid
import pyramid.magcreator as mc
from pyramid.projector import SimpleProjector
-from pyramid.phasemapper import PMConvolve, PMElectric
+from pyramid.phasemapper import PhaseMapperRDFC, PhaseMapperElectric
+from pyramid.kernel import Kernel
from pyramid.magdata import MagData
import matplotlib.pyplot as plt
@@ -52,7 +53,7 @@ if not os.path.exists(directory):
# Set parameters:
a = 1.0 # in nm
phi = pi/4
-density = 30
+gain = 30
dim = (128, 128, 128) # in px (z, y, x)
v_0 = 1
v_acc = 300000
@@ -78,12 +79,12 @@ mag_data_slab = MagData(a, mc.create_mag_dist_homog(mag_shape_slab, phi))
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)
+pm_mag = PhaseMapperRDFC(Kernel(a, projector.dim_uv))
+pm_ele = PhaseMapperElectric(a, projector.dim_uv, v_0, v_acc)
# Magnetic phase map of a homogeneously magnetized disc:
-phase_map_mag_disc = pm_mag(mag_data_disc)
+phase_map_mag_disc = pm_mag(projector(mag_data_disc))
phase_map_mag_disc.save_to_txt(directory+'phase_map_mag_disc.txt')
-axis, _ = phase_map_mag_disc.display_combined(density=density)
+axis, _ = phase_map_mag_disc.display_combined(gain=gain)
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')
@@ -97,9 +98,9 @@ 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 = pm_mag(mag_data_vortex)
+phase_map_mag_vortex = pm_mag(projector(mag_data_vortex))
phase_map_mag_vortex.save_to_txt(directory+'phase_map_mag_vortex.txt')
-axis, _ = phase_map_mag_vortex.display_combined(density=density)
+axis, _ = phase_map_mag_vortex.display_combined(gain=gain)
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')
@@ -113,9 +114,9 @@ plt.figure()
plt.plot(x, y)
plt.savefig(directory+'charge_integral_mag_vortex.png')
# MIP phase of a slab:
-phase_map_mip_slab = pm_ele(mag_data_slab)
+phase_map_mip_slab = pm_ele(projector(mag_data_slab))
phase_map_mip_slab.save_to_txt(directory+'phase_map_mip_slab.txt')
-axis, _ = phase_map_mip_slab.display_combined(density=density)
+axis, _ = phase_map_mip_slab.display_combined(gain=gain)
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')
@@ -129,9 +130,9 @@ plt.figure()
plt.plot(x, y)
plt.savefig(directory+'charge_integral_mip_slab.png')
# MIP phase of a disc:
-phase_map_mip_disc = pm_ele(mag_data_disc)
+phase_map_mip_disc = pm_ele(projector(mag_data_disc))
phase_map_mip_disc.save_to_txt(directory+'phase_map_mip_disc.txt')
-axis, _ = phase_map_mip_disc.display_combined(density=density)
+axis, _ = phase_map_mip_disc.display_combined(gain=gain)
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')
@@ -145,9 +146,9 @@ plt.figure()
plt.plot(x, y)
plt.savefig(directory+'charge_integral_mip_disc.png')
# MIP phase of a sphere:
-phase_map_mip_sphere = pm_ele(mag_data_sphere)
+phase_map_mip_sphere = pm_ele(projector(mag_data_sphere))
phase_map_mip_sphere.save_to_txt(directory+'phase_map_mip_sphere.txt')
-axis, _ = phase_map_mip_sphere.display_combined(density=density)
+axis, _ = phase_map_mip_sphere.display_combined(gain=gain)
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')
@@ -161,9 +162,9 @@ plt.figure()
plt.plot(x, y)
plt.savefig(directory+'charge_integral_mip_sphere.png')
# MIP phase of an ellipsoid:
-phase_map_mip_ellipsoid = pm_ele(mag_data_ellipsoid)
+phase_map_mip_ellipsoid = pm_ele(projector(mag_data_ellipsoid))
phase_map_mip_ellipsoid.save_to_txt(directory+'phase_map_mip_ellipsoid.txt')
-axis, _ = phase_map_mip_ellipsoid.display_combined(phase_map_mip_ellipsoid, density)
+axis, _ = phase_map_mip_ellipsoid.display_combined(gain=gain)
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')
diff --git a/scripts/create distributions/create_alternating_filaments.py b/scripts/create distributions/create_alternating_filaments.py
index 8dd5557..3729d52 100644
--- a/scripts/create distributions/create_alternating_filaments.py
+++ b/scripts/create distributions/create_alternating_filaments.py
@@ -11,6 +11,7 @@ from numpy import pi
import pyramid
import pyramid.magcreator as mc
from pyramid.magdata import MagData
+from pyramid.phasemapper import pm
import logging
import logging.config
@@ -40,3 +41,4 @@ for i in range(count):
# Plot magnetic distribution
mag_data.quiver_plot()
mag_data.save_to_llg(filename)
+pm(mag_data).display_phase()
diff --git a/scripts/create distributions/create_core_shell_disc.py b/scripts/create distributions/create_core_shell_disc.py
index 930a40e..96b1be5 100644
--- a/scripts/create distributions/create_core_shell_disc.py
+++ b/scripts/create distributions/create_core_shell_disc.py
@@ -11,8 +11,7 @@ import matplotlib.pyplot as plt
import pyramid
import pyramid.magcreator as mc
-from pyramid.projector import SimpleProjector
-from pyramid.phasemapper import PMConvolve
+from pyramid.phasemapper import pm
from pyramid.magdata import MagData
import logging
@@ -44,10 +43,10 @@ mag_data.quiver_plot('z-projection', proj_axis='z')
mag_data.quiver_plot('y-projection', proj_axis='y')
mag_data.save_to_llg(filename)
mag_data.quiver_plot3d()
-phase_map_z = PMConvolve(a, SimpleProjector(dim, axis='z'))(mag_data)
-phase_map_x = PMConvolve(a, SimpleProjector(dim, axis='x'))(mag_data)
-phase_map_z.display_holo('Core-Shell structure (z-proj.)', density=1E2)
-phase_axis, holo_axis = phase_map_x.display_combined('Core-Shell structure (x-proj.)', density=1E3)
+phase_map_z = pm(mag_data)
+phase_map_x = pm(mag_data)
+phase_map_z.display_holo('Core-Shell structure (z-proj.)', gain=1E2)
+phase_axis, holo_axis = phase_map_x.display_combined('Core-Shell structure (x-proj.)', gain=1E3)
phase_axis.set_xlabel('y [nm]')
phase_axis.set_ylabel('z [nm]')
holo_axis.set_xlabel('y-axis [px]')
diff --git a/scripts/create distributions/create_core_shell_sphere.py b/scripts/create distributions/create_core_shell_sphere.py
index d6033cc..662752a 100644
--- a/scripts/create distributions/create_core_shell_sphere.py
+++ b/scripts/create distributions/create_core_shell_sphere.py
@@ -11,8 +11,7 @@ import matplotlib.pyplot as plt
import pyramid
import pyramid.magcreator as mc
-from pyramid.projector import SimpleProjector
-from pyramid.phasemapper import PMConvolve
+from pyramid.phasemapper import pm
from pyramid.magdata import MagData
import logging
@@ -44,10 +43,10 @@ mag_data.quiver_plot('z-projection', proj_axis='z')
mag_data.quiver_plot('x-projection', proj_axis='x')
mag_data.save_to_llg(filename)
mag_data.quiver_plot3d()
-phase_map_z = PMConvolve(a, SimpleProjector(dim, axis='z'))(mag_data)
-phase_map_x = PMConvolve(a, SimpleProjector(dim, axis='x'))(mag_data)
-phase_map_z.display_holo('Core-Shell structure (z-proj.)', density=1E2)
-phase_axis, holo_axis = phase_map_x.display_combined('Core-Shell structure (x-proj.)', density=1E3)
+phase_map_z = pm(mag_data)
+phase_map_x = pm(mag_data)
+phase_map_z.display_holo('Core-Shell structure (z-proj.)', gain=1E2)
+phase_axis, holo_axis = phase_map_x.display_combined('Core-Shell structure (x-proj.)', gain=1E3)
phase_axis.set_xlabel('y [nm]')
phase_axis.set_ylabel('z [nm]')
holo_axis.set_xlabel('y-axis [px]')
diff --git a/scripts/create distributions/create_multiple_samples.py b/scripts/create distributions/create_multiple_samples.py
index 1900007..1c09a0e 100644
--- a/scripts/create distributions/create_multiple_samples.py
+++ b/scripts/create distributions/create_multiple_samples.py
@@ -10,8 +10,7 @@ from numpy import pi
import pyramid
import pyramid.magcreator as mc
from pyramid.magdata import MagData
-from pyramid.phasemapper import PMConvolve
-from pyramid.projector import SimpleProjector
+from pyramid.phasemapper import pm
import logging
import logging.config
@@ -48,5 +47,5 @@ mag_data += MagData(a, mc.create_mag_dist_homog(mag_shape_sphere, pi))
# Plot the magnetic distribution, phase map and holographic contour map:
mag_data.quiver_plot()
mag_data.save_to_llg(filename)
-phase_map = PMConvolve(a, SimpleProjector(dim))(mag_data)
-phase_map.display_combined(gain='auto')
+phase_map = pm(mag_data)
+phase_map.display_combined(interpolation='bilinear')
diff --git a/scripts/create distributions/create_random_pixels.py b/scripts/create distributions/create_random_pixels.py
index 90479de..1869997 100644
--- a/scripts/create distributions/create_random_pixels.py
+++ b/scripts/create distributions/create_random_pixels.py
@@ -12,8 +12,7 @@ from numpy import pi
import pyramid
import pyramid.magcreator as mc
from pyramid.magdata import MagData
-from pyramid.phasemapper import PMConvolve
-from pyramid.projector import SimpleProjector
+from pyramid.phasemapper import pm
import logging
import logging.config
@@ -44,5 +43,5 @@ for i in range(count):
# Plot magnetic distribution, phase map and holographic contour map:
mag_data.quiver_plot()
mag_data.save_to_llg(filename)
-phase_map = PMConvolve(a, SimpleProjector(dim))(mag_data)
-phase_map.display_combined(density=10)
+phase_map = pm(mag_data)
+phase_map.display_combined()
diff --git a/scripts/create distributions/create_random_slabs.py b/scripts/create distributions/create_random_slabs.py
index eaba13f..510ce61 100644
--- a/scripts/create distributions/create_random_slabs.py
+++ b/scripts/create distributions/create_random_slabs.py
@@ -12,8 +12,7 @@ from numpy import pi
import pyramid
import pyramid.magcreator as mc
from pyramid.magdata import MagData
-from pyramid.phasemapper import PMConvolve
-from pyramid.projector import SimpleProjector
+from pyramid.phasemapper import pm
import logging
import logging.config
@@ -47,5 +46,5 @@ for i in range(count):
# Plot magnetic distribution, phase map and holographic contour map:
mag_data.quiver_plot()
mag_data.save_to_llg(filename)
-phase_map = PMConvolve(a, SimpleProjector(dim))(mag_data)
-phase_map.display_combined(density=10)
+phase_map = pm(mag_data)
+phase_map.display_combined()
diff --git a/scripts/create distributions/create_vortex.py b/scripts/create distributions/create_vortex.py
index 9093588..7e5f6bb 100644
--- a/scripts/create distributions/create_vortex.py
+++ b/scripts/create distributions/create_vortex.py
@@ -38,7 +38,7 @@ mag_data = MagData(a, mc.create_mag_dist_vortex(mag_shape, magnitude=0.75))
mag_data.quiver_plot()
mag_data.save_to_llg(filename)
phase_map = pm(mag_data)
-phase_map.display_combined()
+phase_map.display_combined(gain=2)
phase_slice = phase_map.phase[dim[1]/2, :]
fig = plt.figure()
fig.add_subplot(111)
diff --git a/scripts/gui/mag_slicer.py b/scripts/gui/mag_slicer.py
index bfe8119..353f91d 100644
--- a/scripts/gui/mag_slicer.py
+++ b/scripts/gui/mag_slicer.py
@@ -16,7 +16,8 @@ from matplotlibwidget import MatplotlibWidget
from pyramid.magdata import MagData
from pyramid.projector import SimpleProjector
-from pyramid.phasemapper import PMConvolve
+from pyramid.phasemapper import PhaseMapperRDFC
+from pyramid.kernel import Kernel
try:
@@ -206,7 +207,8 @@ class UI_MagSlicerMain(QtGui.QWidget):
self.scrollBarSlice.setMaximum(length)
self.spinBoxSlice.setValue(int(length/2.))
self.update_slice()
- self.phase_map = PMConvolve(self.mag_data.a, self.projector)(self.mag_data)
+ self.phase_mapper = PhaseMapperRDFC(Kernel(self.mag_data.a, self.projector.dim_uv))
+ self.phase_map = self.phase_mapper(self.projector(self.mag_data))
self.phase_map.display_phase(axis=self.mplWidgetPhase.axes, cbar=False, show=False)
if self.checkBoxSmooth.isChecked():
interpolation = 'bilinear'
diff --git a/scripts/publications/abstract_prague.py b/scripts/publications/abstract_prague.py
index cb5bfe7..e62eaca 100644
--- a/scripts/publications/abstract_prague.py
+++ b/scripts/publications/abstract_prague.py
@@ -13,8 +13,10 @@ 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
+from pyramid.projector import YTiltProjector
+from pyramid.phasemapper import PhaseMapperRDFC, pm
+from pyramid.kernel import Kernel
+from pyramid.dataset import DataSet
import logging
import logging.config
@@ -37,9 +39,7 @@ magnitude = mc.create_mag_dist_homog(mc.Shapes.sphere(dim, center, radius), pi/4
mag_data = MagData(a, magnitude)
-projector = SimpleProjector(dim)
-
-phase_map = PMConvolve(a, projector)(mag_data)
+phase_map = pm(mag_data)
axis = phase_map.display_phase()
axis.tick_params(axis='both', which='major', labelsize=18)
@@ -47,7 +47,7 @@ axis.set_title('Phase map', fontsize=24)
axis.set_xlabel('x [nm]', fontsize=18)
axis.set_ylabel('y [nm]', fontsize=18)
-axis = phase_map.display_holo(density=20, interpolation='bilinear')
+axis = phase_map.display_holo(gain=20, interpolation='bilinear')
axis.tick_params(axis='both', which='major', labelsize=18)
axis.set_title('Magnetic induction map', fontsize=24)
axis.set_xlabel('x [nm]', fontsize=18)
@@ -88,12 +88,11 @@ mag_data = MagData(a, mc.create_mag_dist_homog(mc.Shapes.ellipse(dim, center, (2
tilts = np.array([0., 60.])/180.*pi
-projectors = [YTiltProjector(mag_data.dim, tilt) for tilt in tilts]
-phasemappers = [PMConvolve(mag_data.a, proj, b_0) for proj in projectors]
-
-phase_maps = [pm(mag_data) for pm in phasemappers]
+data_set = DataSet(a, dim, b_0)
+data_set.projectors = [YTiltProjector(mag_data.dim, tilt) for tilt in tilts]
+phase_maps = data_set.create_phase_maps(mag_data)
-axis = phase_maps[0].display_holo(density=1, interpolation='bilinear')
+axis = phase_maps[0].display_holo(gain=1, interpolation='bilinear')
axis.tick_params(axis='both', which='major', labelsize=18)
axis.set_title('Magnetic induction map', fontsize=24)
axis.xaxis.set_major_formatter(FuncFormatter(lambda x, pos: '{:g}'.format(x*a)))
@@ -107,7 +106,7 @@ axis.set_title('Phase map', fontsize=24)
axis.set_xlabel('x [nm]', fontsize=18)
axis.set_ylabel('y [nm]', fontsize=18)
-axis = phase_maps[1].display_holo(density=1, interpolation='bilinear')
+axis = phase_maps[1].display_holo(gain=1, interpolation='bilinear')
axis.tick_params(axis='both', which='major', labelsize=18)
axis.set_title('Magnetic induction map', fontsize=24)
axis.xaxis.set_major_formatter(FuncFormatter(lambda x, pos: '{:g}'.format(x*a)))
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 7f6d92c..763466d 100644
--- a/scripts/publications/paper 1/ch5-1-evaluation_real_space.py
+++ b/scripts/publications/paper 1/ch5-1-evaluation_real_space.py
@@ -16,8 +16,7 @@ from matplotlib.patches import Rectangle
import pyramid
import pyramid.magcreator as mc
import pyramid.analytic as an
-from pyramid.projector import SimpleProjector
-from pyramid.phasemapper import PMConvolve
+from pyramid.phasemapper import pm
from pyramid.magdata import MagData
from pyramid.phasemap import PhaseMap
@@ -121,8 +120,7 @@ else:
mag_data_disc = MagData(a, mc.create_mag_dist_homog(mag_shape, phi))
for i in range(5):
print '----a =', mag_data_disc.a, 'nm', 'dim =', mag_data_disc.dim
- projector = SimpleProjector(mag_data_disc.dim)
- phase_map = PMConvolve(mag_data_disc.a, projector)(mag_data_disc)
+ phase_map = pm(mag_data_disc)
phase_map.display_combined('Disc, a = {} nm'.format(mag_data_disc.a), gain=gain)
x_d.append(np.linspace(mag_data_disc.a * 0.5,
mag_data_disc.a * (mag_data_disc.dim[1]-0.5),
@@ -155,8 +153,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)
- phase_map = PMConvolve(mag_data_vort.a, projector)(mag_data_vort)
+ phase_map = pm(mag_data_vort)
phase_map.display_combined('Disc, a = {} nm'.format(mag_data_vort.a), gain=gain)
x_v.append(np.linspace(mag_data_vort.a * 0.5,
mag_data_vort.a * (mag_data_vort.dim[1]-0.5),
@@ -322,9 +319,6 @@ else:
data_shelve[key] = (mag_data_disc, mag_data_vort)
print '--CALCULATE PHASE DIFFERENCES'
-# Create projector along z-axis and phasemapper:
-projector = SimpleProjector(dim)
-phasemapper = PMConvolve(a, projector)
# Get analytic solutions:
phase_map_ana_disc = an.phase_mag_disc(dim, a, phi, center, radius, height)
phase_map_ana_vort = an.phase_mag_vortex(dim, a, center, radius, height)
@@ -332,12 +326,12 @@ phase_map_ana_vort = an.phase_mag_vortex(dim, a, center, radius, height)
bounds = np.array([-3, -0.5, -0.25, -0.1, 0, 0.1, 0.25, 0.5, 3])
norm = BoundaryNorm(bounds, RdBu.N)
# Calculations (Disc):
-phase_map_num_disc = phasemapper(mag_data_disc)
+phase_map_num_disc = pm(mag_data_disc)
phase_map_num_disc.unit = 'mrad'
phase_diff_disc = phase_map_num_disc - phase_map_ana_disc
RMS_disc = np.sqrt(np.mean(phase_diff_disc.phase**2))
# Calculations (Vortex):
-phase_map_num_vort = phasemapper(mag_data_vort)
+phase_map_num_vort = pm(mag_data_vort)
phase_map_num_vort.unit = 'mrad'
phase_diff_vort = phase_map_num_vort - phase_map_ana_vort
RMS_vort = np.sqrt(np.mean(phase_diff_vort.phase**2))
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 45b5370..78f731f 100644
--- a/scripts/publications/paper 1/ch5-2-evaluation_fourier_space.py
+++ b/scripts/publications/paper 1/ch5-2-evaluation_fourier_space.py
@@ -17,7 +17,7 @@ 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
+from pyramid.phasemapper import PhaseMapperFDFC
import time
import shelve
@@ -91,8 +91,10 @@ else:
for i in range(5):
print '----a =', mag_data_disc.a, 'nm', 'dim =', mag_data_disc.dim
projector = SimpleProjector(mag_data_disc.dim)
- phase_map0 = PMFourier(mag_data_disc.a, projector, padding=0)(mag_data_disc)
- phase_map10 = PMFourier(mag_data_disc.a, projector, padding=10)(mag_data_disc)
+ phase_map0 = PhaseMapperFDFC(mag_data_disc.a, projector.dim_uv,
+ padding=0)(projector(mag_data_disc))
+ phase_map10 = PhaseMapperFDFC(mag_data_disc.a, projector.dim_uv,
+ padding=10)(projector(mag_data_disc))
phase_map0.unit = 'mrad'
phase_map10.unit = 'mrad'
phase_map0.display_combined('Disc, a = {} nm'.format(mag_data_disc.a), gain=gain)
@@ -135,8 +137,10 @@ else:
for i in range(5):
print '----a =', mag_data_vort.a, 'nm', 'dim =', mag_data_vort.dim
projector = SimpleProjector(mag_data_vort.dim)
- phase_map0 = PMFourier(mag_data_vort.a, projector, padding=0)(mag_data_vort)
- phase_map10 = PMFourier(mag_data_vort.a, projector, padding=10)(mag_data_vort)
+ phase_map0 = PhaseMapperFDFC(mag_data_vort.a, projector.dim_uv,
+ padding=0)(projector(mag_data_vort))
+ phase_map10 = PhaseMapperFDFC(mag_data_vort.a, projector.dim_uv,
+ padding=10)(projector(mag_data_vort))
phase_map0.unit = 'mrad'
phase_map10.unit = 'mrad'
print np.mean(phase_map0.phase)
@@ -433,7 +437,7 @@ phase_map_ana_disc = an.phase_mag_disc(dim, a, phi, center, radius, height)
phase_map_ana_vort = an.phase_mag_vortex(dim, a, center, radius, height)
# Create phasemapper:
-phasemapper = PMFourier(a, projector, padding=0)
+phasemapper = PhaseMapperFDFC(a, projector.dim_uv, padding=0)
# Create figure:
fig, axes = plt.subplots(1, 2, figsize=(16, 7))
fig.suptitle('Difference of the Fourier space approach from the analytical solution', fontsize=20)
@@ -441,7 +445,7 @@ fig.suptitle('Difference of the Fourier space approach from the analytical solut
bounds = np.array([-100, -50, -25, -5, 0, 5, 25, 50, 100])
norm = BoundaryNorm(bounds, RdBu.N)
# Disc:
-phase_map_num_disc = phasemapper(mag_data_disc)
+phase_map_num_disc = phasemapper(projector(mag_data_disc))
phase_map_num_disc.unit = 'mrad'
phase_diff_disc = phase_map_num_disc - phase_map_ana_disc
RMS_disc = np.sqrt(np.mean(phase_diff_disc.phase**2))
@@ -453,7 +457,7 @@ axes[0].set_aspect('equal')
bounds = np.array([-3, -0.5, -0.25, -0.1, 0, 0.1, 0.25, 0.5, 3])
norm = BoundaryNorm(bounds, RdBu.N)
# Vortex:
-phase_map_num_vort = phasemapper(mag_data_vort)
+phase_map_num_vort = phasemapper(projector(mag_data_vort))
phase_map_num_vort.unit = 'mrad'
phase_diff_vort = phase_map_num_vort - phase_map_ana_vort
phase_diff_vort -= np.mean(phase_diff_vort.phase)
@@ -468,7 +472,7 @@ plt.figtext(0.52, 0.2, 'b)', fontsize=30)
plt.savefig(directory + '/ch5-2-fourier_phase_differe_no_padding.png', bbox_inches='tight')
# Create phasemapper:
-phasemapper = PMFourier(a, projector, padding=10)
+phasemapper = PhaseMapperFDFC(a, projector.dim_uv, padding=10)
# Create figure:
fig, axes = plt.subplots(1, 2, figsize=(16, 7))
fig.suptitle('Difference of the Fourier space approach from the analytical solution', fontsize=20)
@@ -476,7 +480,7 @@ fig.suptitle('Difference of the Fourier space approach from the analytical solut
bounds = np.array([-3, -0.5, -0.25, -0.1, 0, 0.1, 0.25, 0.5, 3])
norm = BoundaryNorm(bounds, RdBu.N)
# Disc:
-phase_map_num_disc = phasemapper(mag_data_disc)
+phase_map_num_disc = phasemapper(projector(mag_data_disc))
phase_map_num_disc.unit = 'mrad'
phase_diff_disc = phase_map_num_disc - phase_map_ana_disc
RMS_disc = np.sqrt(np.mean(phase_diff_disc.phase**2))
@@ -488,7 +492,7 @@ axes[0].set_aspect('equal')
bounds = np.array([-3, -0.5, -0.25, -0.1, 0, 0.1, 0.25, 0.5, 3])
norm = BoundaryNorm(bounds, RdBu.N)
# Vortex:
-phase_map_num_vort = phasemapper(mag_data_vort)
+phase_map_num_vort = phasemapper(projector(mag_data_vort))
phase_map_num_vort.unit = 'mrad'
phase_diff_vort = phase_map_num_vort - phase_map_ana_vort
phase_diff_vort -= np.mean(phase_diff_vort.phase)
@@ -555,9 +559,9 @@ for (i, padding) in enumerate(padding_list):
data_disc[:, i] = data_shelve[key]
else:
print '----calculate and save padding =', padding_list[i]
- phasemapper = PMFourier(a, projector, padding=padding_list[i])
+ phasemapper = PhaseMapperFDFC(a, projector.dim_uv, padding=padding_list[i])
start_time = time.time()
- phase_num_disc = phasemapper(mag_data_disc)
+ phase_num_disc = phasemapper(projector(mag_data_disc))
data_disc[2, i] = time.time() - start_time
phase_map_diff = phase_num_disc - phase_ana_disc
phase_map_diff.unit = 'mrad'
@@ -612,9 +616,9 @@ for (i, padding) in enumerate(padding_list):
data_vort[:, i] = data_shelve[key]
else:
print '----calculate and save padding =', padding_list[i]
- phasemapper = PMFourier(a, projector, padding=padding_list[i])
+ phasemapper = PhaseMapperFDFC(a, projector.dim_uv, padding=padding_list[i])
start_time = time.time()
- phase_num_vort = phasemapper(mag_data_vort)
+ phase_num_vort = phasemapper(projector(mag_data_vort))
data_vort[2, i] = time.time() - start_time
phase_map_diff = phase_num_vort - phase_ana_vort
phase_map_diff -= np.mean(phase_map_diff.phase)
diff --git a/scripts/publications/paper 1/ch5-3-comparison_of_methods.py b/scripts/publications/paper 1/ch5-3-comparison_of_methods.py
deleted file mode 100644
index f41c80e..0000000
--- a/scripts/publications/paper 1/ch5-3-comparison_of_methods.py
+++ /dev/null
@@ -1,263 +0,0 @@
-# -*- coding: utf-8 -*-
-"""Created on Fri Jul 26 14:37:30 2013 @author: Jan"""
-
-
-import os
-
-import numpy as np
-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:
-directory = '../../../output/paper 1'
-if not os.path.exists(directory):
- os.makedirs(directory)
-data_shelve = shelve.open(directory + '/paper_1_shelve')
-
-###############################################################################################
-print 'CH5-3 METHOD COMPARISON'
-
-key = 'ch5-3-method_comparison'
-if key in data_shelve:
- print '--LOAD METHOD DATA'
- (data_disc_fourier0, data_vort_fourier0,
- data_disc_fourier1, data_vort_fourier1,
- data_disc_fourier10, data_vort_fourier10,
- data_disc_real_s, data_vort_real_s,
- data_disc_real_d, data_vort_real_d) = data_shelve[key]
-else:
- # Input parameters:
- steps = 6
- a = 0.25 # in nm
- phi = pi/2
- dim = (64, 512, 512) # in px (z, y, x)
- 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]/4 # in px
- height = dim[0]/2 # in px
-
- print '--CREATE MAGNETIC SHAPE'
- mag_shape = mc.Shapes.disc(dim, center, radius, height)
- # Create MagData (4 times the size):
- print '--CREATE MAG. DIST. HOMOG. MAGN. DISC'
- mag_data_disc = MagData(a, mc.create_mag_dist_homog(mag_shape, phi))
- print '--CREATE MAG. DIST. VORTEX STATE DISC'
- mag_data_vort = MagData(a, mc.create_mag_dist_vortex(mag_shape, center))
-
- # Create Data Arrays:
- a_list = [a*2**i for i in np.linspace(1, steps, steps)]
- data_disc_fourier0 = np.vstack((a_list, np.zeros((2, steps))))
- data_vort_fourier0 = np.vstack((a_list, np.zeros((2, steps))))
- data_disc_fourier1 = np.vstack((a_list, np.zeros((2, steps))))
- data_vort_fourier1 = np.vstack((a_list, np.zeros((2, steps))))
- data_disc_fourier10 = np.vstack((a_list, np.zeros((2, steps))))
- data_vort_fourier10 = np.vstack((a_list, np.zeros((2, steps))))
- data_disc_real_s = np.vstack((a_list, np.zeros((2, steps))))
- data_vort_real_s = np.vstack((a_list, np.zeros((2, steps))))
- data_disc_real_d = np.vstack((a_list, np.zeros((2, steps))))
- data_vort_real_d = np.vstack((a_list, np.zeros((2, steps))))
-
- for i in range(steps):
- # Scale mag_data, grid spacing and dimensions:
- mag_data_disc.scale_down()
- mag_data_vort.scale_down()
- dim = mag_data_disc.dim
- a = mag_data_disc.a
- center = (dim[0]/2-0.5, dim[1]/2.-0.5, dim[2]/2.-0.5) # in px, index starts at 0!
- radius = dim[1]/4 # in px
- height = dim[0]/2 # in px
-
- print '--a =', a, 'nm', 'dim =', dim
-
- print '----CALCULATE RMS/DURATION HOMOG. MAGN. DISC'
- # Create projections along z-axis:
- projection_disc = pj.simple_axis_projection(mag_data_disc)
- # Analytic solution:
- phase_ana_disc = an.phase_mag_disc(dim, a, phi, center, radius, height)
- # Fourier unpadded:
- padding = 0
- start_time = time.clock()
- phase_num_disc = pm.phase_mag_fourier(a, projection_disc, padding=padding)
- data_disc_fourier0[2, i] = time.clock() - start_time
- print '------time (disc, fourier0) =', data_disc_fourier0[2, i]
- phase_diff_disc = (phase_ana_disc-phase_num_disc) * 1E3 # in mrad -> *1000
- data_disc_fourier0[1, i] = np.sqrt(np.mean(phase_diff_disc**2))
- # Fourier padding 1:
- padding = 1
- start_time = time.clock()
- phase_num_disc = pm.phase_mag_fourier(a, projection_disc, padding=padding)
- data_disc_fourier1[2, i] = time.clock() - start_time
- print '------time (disc, fourier1) =', data_disc_fourier1[2, i]
- phase_diff_disc = (phase_ana_disc-phase_num_disc) * 1E3 # in mrad -> *1000
- data_disc_fourier1[1, i] = np.sqrt(np.mean(phase_diff_disc**2))
- # Fourier padding 10:
- padding = 10
- start_time = time.clock()
- phase_num_disc = pm.phase_mag_fourier(a, projection_disc, padding=padding)
- data_disc_fourier10[2, i] = time.clock() - start_time
- print '------time (disc, fourier10) =', data_disc_fourier10[2, i]
- phase_diff_disc = (phase_ana_disc-phase_num_disc) * 1E3 # in mrad -> *1000
- data_disc_fourier10[1, i] = np.sqrt(np.mean(phase_diff_disc**2))
- # Real space slab:
- start_time = time.clock()
- phase_num_disc = pm.phase_mag_real(a, projection_disc, geometry='slab')
- data_disc_real_s[2, i] = time.clock() - start_time
- print '------time (disc, real slab) =', data_disc_real_s[2, i]
- phase_diff_disc = (phase_ana_disc-phase_num_disc) * 1E3 # in mrad -> *1000
- data_disc_real_s[1, i] = np.sqrt(np.mean(phase_diff_disc**2))
- # Real space disc:
- start_time = time.clock()
- phase_num_disc = pm.phase_mag_real(a, projection_disc, geometry='disc')
- data_disc_real_d[2, i] = time.clock() - start_time
- print '------time (disc, real disc) =', data_disc_real_d[2, i]
- phase_diff_disc = (phase_ana_disc-phase_num_disc) * 1E3 # in mrad -> *1000
- data_disc_real_d[1, i] = np.sqrt(np.mean(phase_diff_disc**2))
-
- print '----CALCULATE RMS/DURATION VORTEX STATE DISC'
- # Create projections along z-axis:
- projection_vort = pj.simple_axis_projection(mag_data_vort)
- # Analytic solution:
- phase_ana_vort = an.phase_mag_vortex(dim, a, center, radius, height)
- # Fourier unpadded:
- padding = 0
- start_time = time.clock()
- phase_num_vort = pm.phase_mag_fourier(a, projection_vort, padding=padding)
- data_vort_fourier0[2, i] = time.clock() - start_time
- print '------time (vortex, fourier0) =', data_vort_fourier0[2, i]
- phase_diff_vort = (phase_ana_vort-phase_num_vort) * 1E3 # in mrad -> *1000
- data_vort_fourier0[1, i] = np.sqrt(np.mean(phase_diff_vort**2))
- # Fourier padding 1:
- padding = 1
- start_time = time.clock()
- phase_num_vort = pm.phase_mag_fourier(a, projection_vort, padding=padding)
- data_vort_fourier1[2, i] = time.clock() - start_time
- print '------time (vortex, fourier1) =', data_vort_fourier1[2, i]
- phase_diff_vort = (phase_ana_vort-phase_num_vort) * 1E3 # in mrad -> *1000
- data_vort_fourier1[1, i] = np.sqrt(np.mean(phase_diff_vort**2))
- # Fourier padding 10:
- padding = 10
- start_time = time.clock()
- phase_num_vort = pm.phase_mag_fourier(a, projection_vort, padding=padding)
- data_vort_fourier10[2, i] = time.clock() - start_time
- print '------time (vortex, fourier10) =', data_vort_fourier10[2, i]
- phase_diff_vort = (phase_ana_vort-phase_num_vort) * 1E3 # in mrad -> *1000
- data_vort_fourier10[1, i] = np.sqrt(np.mean(phase_diff_vort**2))
- # Real space slab:
- start_time = time.clock()
- phase_num_vort = pm.phase_mag_real(a, projection_vort, geometry='slab')
- data_vort_real_s[2, i] = time.clock() - start_time
- print '------time (vortex, real slab) =', data_vort_real_s[2, i]
- phase_diff_vort = (phase_ana_vort-phase_num_vort) * 1E3 # in mrad -> *1000
- data_vort_real_s[1, i] = np.sqrt(np.mean(phase_diff_vort**2))
- # Real space disc:
- start_time = time.clock()
- phase_num_vort = pm.phase_mag_real(a, projection_vort, geometry='disc')
- data_vort_real_d[2, i] = time.clock() - start_time
- print '------time (vortex, real disc) =', data_vort_real_d[2, i]
- phase_diff_vort = (phase_ana_vort-phase_num_vort) * 1E3 # in mrad -> *1000
- data_vort_real_d[1, i] = np.sqrt(np.mean(phase_diff_vort**2))
-
- print '--SHELVE METHOD DATA'
- data_shelve[key] = (data_disc_fourier0, data_vort_fourier0,
- data_disc_fourier1, data_vort_fourier1,
- data_disc_fourier10, data_vort_fourier10,
- data_disc_real_s, data_vort_real_s,
- data_disc_real_d, data_vort_real_d)
-
-print '--PLOT/SAVE METHOD DATA'
-
-# Plot using shared rows and colums:
-fig, axes = plt.subplots(2, 2, sharex='col', sharey='row', figsize=(12, 8))
-fig.tight_layout(rect=(0.05, 0.05, 0.95, 0.95))
-fig.suptitle('Method Comparison', fontsize=20)
-
-# Plot duration against a (disc) [top/left]:
-axes[1, 0].set_yscale('log')
-axes[1, 0].plot(data_disc_fourier0[0], data_disc_fourier0[1], ':bs')
-axes[1, 0].plot(data_disc_fourier1[0], data_disc_fourier1[1], ':bo')
-axes[1, 0].plot(data_disc_fourier10[0], data_disc_fourier10[1], ':b^')
-axes[1, 0].plot(data_disc_real_s[0], data_disc_real_s[1], '--rs')
-axes[1, 0].plot(data_disc_real_d[0], data_disc_real_d[1], '--ro')
-axes[1, 0].set_xlabel('a [nm]', fontsize=15)
-axes[1, 0].set_ylabel('RMS [mrad]', fontsize=15)
-axes[1, 0].set_xlim(-0.5, 16.5)
-axes[1, 0].tick_params(axis='both', which='major', labelsize=14)
-axes[1, 0].xaxis.set_major_locator(IndexLocator(base=4, offset=-0.5))
-
-# Plot RMS against a (disc) [bottom/left]:
-plt.tick_params(axis='both', which='major', labelsize=14)
-axes[0, 0].set_yscale('log')
-axes[0, 0].plot(data_disc_fourier0[0], data_disc_fourier0[2], ':bs')
-axes[0, 0].plot(data_disc_fourier1[0], data_disc_fourier1[2], ':bo')
-axes[0, 0].plot(data_disc_fourier10[0], data_disc_fourier10[2], ':b^')
-axes[0, 0].plot(data_disc_real_s[0], data_disc_real_s[2], '--rs')
-axes[0, 0].plot(data_disc_real_d[0], data_disc_real_d[2], '--ro')
-axes[0, 0].set_title('Homog. magn. disc', fontsize=18)
-axes[0, 0].set_ylabel('duration [s]', fontsize=15)
-axes[0, 0].set_xlim(-0.5, 16.5)
-axes[0, 0].tick_params(axis='both', which='major', labelsize=14)
-axes[0, 0].xaxis.set_major_locator(IndexLocator(base=4, offset=-0.5))
-
-# Plot duration against a (vortex) [top/right]:
-axes[1, 1].set_yscale('log')
-axes[1, 1].plot(data_vort_fourier0[0], data_vort_fourier0[1], ':bs')
-axes[1, 1].plot(data_vort_fourier1[0], data_vort_fourier1[1], ':bo')
-axes[1, 1].plot(data_vort_fourier10[0], data_vort_fourier10[1], ':b^')
-axes[1, 1].plot(data_vort_real_s[0], data_vort_real_s[1], '--rs')
-axes[1, 1].plot(data_vort_real_d[0], data_vort_real_d[1], '--ro')
-axes[1, 1].set_xlabel('a [nm]', fontsize=15)
-axes[1, 1].set_xlim(-0.5, 16.5)
-axes[1, 1].tick_params(axis='both', which='major', labelsize=14)
-axes[1, 1].xaxis.set_major_locator(IndexLocator(base=4, offset=-0.5))
-
-# Plot RMS against a (vortex) [bottom/right]:
-axes[0, 1].set_yscale('log')
-axes[0, 1].plot(data_vort_fourier0[0], data_vort_fourier0[2], ':bs',
- label='Fourier padding=0')
-axes[0, 1].plot(data_vort_fourier1[0], data_vort_fourier1[2], ':bo',
- label='Fourier padding=1')
-axes[0, 1].plot(data_vort_fourier10[0], data_vort_fourier10[2], ':b^',
- label='Fourier padding=10')
-axes[0, 1].plot(data_vort_real_s[0], data_vort_real_s[2], '--rs',
- label='Real space (slab)')
-axes[0, 1].plot(data_vort_real_d[0], data_vort_real_d[2], '--ro',
- label='Real space (disc)')
-axes[0, 1].set_title('Vortex state disc', fontsize=18)
-axes[0, 1].set_xlim(-0.5, 16.5)
-axes[0, 1].tick_params(axis='both', which='major', labelsize=14)
-axes[0, 1].xaxis.set_major_locator(IndexLocator(base=4, offset=-0.5))
-axes[0, 1].legend(loc=1)
-
-# Save figure as .png:
-plt.show()
-plt.figtext(0.45, 0.85, 'a)', fontsize=30)
-plt.figtext(0.57, 0.85, 'b)', fontsize=30)
-plt.figtext(0.45, 0.15, 'c)', fontsize=30)
-plt.figtext(0.57, 0.15, 'd)', fontsize=30)
-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/paper 1/ch5-4-comparison_of_methods_new.py b/scripts/publications/paper 1/ch5-4-comparison_of_methods_new.py
index 58a3382..bb8ecf0 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
@@ -14,7 +14,8 @@ 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
+from pyramid.phasemapper import PhaseMapperRDFC, PhaseMapperFDFC
+from pyramid.kernel import Kernel
import time
import shelve
@@ -28,19 +29,19 @@ LOGGING_CONF = os.path.join(os.path.dirname(os.path.realpath(pyramid.__file__)),
logging.config.fileConfig(LOGGING_CONF, disable_existing_loggers=False)
-force_calculation = True
+force_calculation = False
-n = 1
+n = 50
-def get_time(pm, mag_data, n):
+def get_time(pm, proj, mag_data, n):
t = []
for i in range(n):
start = time.clock()
- pm(mag_data)
+ pm(proj(mag_data))
t.append(time.clock()-start)
t, dt = np.mean(t), np.std(t)
- phase_map = pm(mag_data)
+ phase_map = pm(proj(mag_data))
return phase_map, t, dt
@@ -102,37 +103,37 @@ else:
# Create projector along z-axis and phasemapper:
projector = SimpleProjector(dim)
- pm_fourier0 = PMFourier(a, projector, padding=0)
- pm_fourier3 = PMFourier(a, projector, padding=3)
- pm_fourier10 = PMFourier(a, projector, padding=7)
- pm_real = PMConvolve(a, projector)
+ pm_fourier0 = PhaseMapperFDFC(a, projector.dim_uv, padding=0)
+ pm_fourier3 = PhaseMapperFDFC(a, projector.dim_uv, padding=3)
+ pm_fourier10 = PhaseMapperFDFC(a, projector.dim_uv, padding=7)
+ pm_real = PhaseMapperRDFC(Kernel(a, projector.dim_uv))
print '----CALCULATE RMS/DURATION HOMOG. MAGN. DISC'
# Analytic solution:
phase_ana_disc = an.phase_mag_disc(dim, a, phi, center, radius, height)
# Fourier unpadded:
- phase_num_disc, t, dt = get_time(pm_fourier0, mag_data_disc, n)
+ phase_num_disc, t, dt = get_time(pm_fourier0, projector, mag_data_disc, n)
data_disc_fourier0[2, i], data_disc_fourier0[3, i] = t, dt
print '------time (disc, fourier0) =', data_disc_fourier0[2, i]
phase_diff_disc = (phase_ana_disc-phase_num_disc) * 1E3 # in mrad -> *1000
print 'phase mean:', np.mean(phase_num_disc.phase)
data_disc_fourier0[1, i] = np.sqrt(np.mean(phase_diff_disc.phase**2))
# Fourier padding 3:
- phase_num_disc, t, dt = get_time(pm_fourier3, mag_data_disc, n)
+ phase_num_disc, t, dt = get_time(pm_fourier3, projector, mag_data_disc, n)
data_disc_fourier3[2, i], data_disc_fourier3[3, i] = t, dt
print '------time (disc, fourier3) =', data_disc_fourier3[2, i]
phase_diff_disc = (phase_ana_disc-phase_num_disc) * 1E3 # in mrad -> *1000
print 'phase mean:', np.mean(phase_num_disc.phase)
data_disc_fourier3[1, i] = np.sqrt(np.mean(phase_diff_disc.phase**2))
# Fourier padding 10:
- phase_num_disc, t, dt = get_time(pm_fourier10, mag_data_disc, n)
+ phase_num_disc, t, dt = get_time(pm_fourier10, projector, mag_data_disc, n)
data_disc_fourier10[2, i], data_disc_fourier10[3, i] = t, dt
print '------time (disc, fourier10) =', data_disc_fourier10[2, i]
phase_diff_disc = (phase_ana_disc-phase_num_disc) * 1E3 # in mrad -> *1000
print 'phase mean:', np.mean(phase_num_disc.phase)
data_disc_fourier10[1, i] = np.sqrt(np.mean(phase_diff_disc.phase**2))
# Real space disc:
- phase_num_disc, t, dt = get_time(pm_real, mag_data_disc, n)
+ phase_num_disc, t, dt = get_time(pm_real, projector, mag_data_disc, n)
data_disc_real[2, i], data_disc_real[3, i] = t, dt
print '------time (disc, real space) =', data_disc_real[2, i]
phase_diff_disc = (phase_ana_disc-phase_num_disc) * 1E3 # in mrad -> *1000
@@ -146,7 +147,7 @@ else:
# Analytic solution:
phase_ana_vort = an.phase_mag_vortex(dim, a, center, radius, height)
# Fourier unpadded:
- phase_num_vort, t, dt = get_time(pm_fourier0, mag_data_vort, n)
+ phase_num_vort, t, dt = get_time(pm_fourier0, projector, mag_data_vort, n)
phase_fft0 = phase_num_vort
data_vort_fourier0[2, i], data_vort_fourier0[3, i] = t, dt
print '------time (vortex, fourier0) =', data_vort_fourier0[2, i]
@@ -155,7 +156,7 @@ else:
phase_diff_vort -= np.mean(phase_diff_vort.phase)
data_vort_fourier0[1, i] = np.sqrt(np.mean(phase_diff_vort.phase**2))
# Fourier padding 3:
- phase_num_vort, t, dt = get_time(pm_fourier3, mag_data_vort, n)
+ phase_num_vort, t, dt = get_time(pm_fourier3, projector, mag_data_vort, n)
phase_fft3 = phase_num_vort
data_vort_fourier3[2, i], data_vort_fourier3[3, i] = t, dt
print '------time (vortex, fourier3) =', data_vort_fourier3[2, i]
@@ -164,7 +165,7 @@ else:
phase_diff_vort -= np.mean(phase_diff_vort.phase)
data_vort_fourier3[1, i] = np.sqrt(np.mean(phase_diff_vort.phase**2))
# Fourier padding 10:
- phase_num_vort, t, dt = get_time(pm_fourier10, mag_data_vort, n)
+ phase_num_vort, t, dt = get_time(pm_fourier10, projector, mag_data_vort, n)
phase_fft10 = phase_num_vort
data_vort_fourier10[2, i], data_vort_fourier10[3, i] = t, dt
print '------time (vortex, fourier10) =', data_vort_fourier10[2, i]
@@ -173,7 +174,7 @@ else:
phase_diff_vort -= np.mean(phase_diff_vort.phase)
data_vort_fourier10[1, i] = np.sqrt(np.mean(phase_diff_vort.phase**2))
# Real space disc:
- phase_num_vort, t, dt = get_time(pm_real, mag_data_vort, n)
+ phase_num_vort, t, dt = get_time(pm_real, projector, mag_data_vort, n)
phase_real = phase_num_vort
data_vort_real[2, i], data_vort_real[3, i] = t, dt
print '------time (vortex, real space) =', data_vort_real[2, i]
diff --git a/scripts/publications/poster_pof.py b/scripts/publications/poster_pof.py
index 24332f9..6ef3d65 100644
--- a/scripts/publications/poster_pof.py
+++ b/scripts/publications/poster_pof.py
@@ -14,8 +14,9 @@ import pyramid
import pyramid.magcreator as mc
import pyramid.analytic as an
from pyramid.magdata import MagData
-from pyramid.phasemapper import PMConvolve
+from pyramid.phasemapper import pm, PhaseMapperRDFC
from pyramid.projector import SimpleProjector
+from pyramid.kernel import Kernel
from pyramid.phasemap import PhaseMap
from time import clock
@@ -35,9 +36,8 @@ 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()
+pm(mag_data_x).display_phase()
+pm(mag_data_y).display_phase()
a = 1
@@ -46,9 +46,10 @@ 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))
+projector = SimpleProjector(dim)
+phasemapper = PhaseMapperRDFC(Kernel(a, projector.dim_uv))
start = clock()
-phase_map = phasemapper(mag_data)
+phase_map = phasemapper(projector(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
diff --git a/scripts/reconstruction/reconstruct_random_pixels.py b/scripts/reconstruction/reconstruct_random_pixels.py
index 65c5289..3aec76e 100644
--- a/scripts/reconstruction/reconstruct_random_pixels.py
+++ b/scripts/reconstruction/reconstruct_random_pixels.py
@@ -12,8 +12,7 @@ from numpy import pi
import pyramid
import pyramid.magcreator as mc
from pyramid.magdata import MagData
-from pyramid.projector import SimpleProjector
-from pyramid.phasemapper import PMConvolve
+from pyramid.phasemapper import pm
import pyramid.reconstruction as rc
import logging
@@ -42,12 +41,12 @@ for i in range(n_pixel):
mag_data += MagData(a, mc.create_mag_dist_homog(mag_shape, phi, magnitude))
# Plot magnetic distribution, phase map and holographic contour map:
mag_data.quiver_plot()
-phase_map = PMConvolve(a, SimpleProjector(dim), b_0)(mag_data)
-phase_map.display_combined('Generated Distribution', density=10)
+phase_map = pm(mag_data)
+phase_map.display_combined('Generated Distribution', gain=10)
# Reconstruct the magnetic distribution:
mag_data_rec = rc.optimize_simple_leastsq(phase_map, mag_data.get_mask(), b_0, lam=1E-4, order=1)
# Display the reconstructed phase map and holography image:
-phase_map_rec = PMConvolve(a, SimpleProjector(dim), b_0)(mag_data_rec)
-phase_map_rec.display_combined('Reconstructed Distribution', density=10)
+phase_map_rec = pm(mag_data_rec)
+phase_map_rec.display_combined('Reconstructed Distribution', gain=10)
diff --git a/scripts/reconstruction/reconstruction_sparse_cg_singular_values.py b/scripts/reconstruction/reconstruction_sparse_cg_singular_values.py
index 6917161..a581395 100644
--- a/scripts/reconstruction/reconstruction_sparse_cg_singular_values.py
+++ b/scripts/reconstruction/reconstruction_sparse_cg_singular_values.py
@@ -56,12 +56,11 @@ 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))]
+data = DataSet(a, dim, 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)
+F = ForwardModel(data)
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))
diff --git a/scripts/reconstruction/reconstruction_sparse_cg_test.py b/scripts/reconstruction/reconstruction_sparse_cg_test.py
index a4bb6a9..48e1471 100644
--- a/scripts/reconstruction/reconstruction_sparse_cg_test.py
+++ b/scripts/reconstruction/reconstruction_sparse_cg_test.py
@@ -13,7 +13,6 @@ import pyramid
from pyramid.magdata import MagData
from pyramid.phasemap import PhaseMap
from pyramid.projector import YTiltProjector, XTiltProjector
-from pyramid.phasemapper import PMConvolve
from pyramid.dataset import DataSet
from pyramid.regularisator import ZeroOrderRegularisator
import pyramid.magcreator as mc
@@ -35,20 +34,19 @@ a = 10.
b_0 = 1.
dim = (8, 32, 32)
dim_uv = dim[1:3]
-count = 8
+count = 16
-#center = (int(dim[0]/2)-0.5, int(dim[1]/2)-0.5, int(dim[2]/2)-0.5)
-#mag_shape = mc.Shapes.disc(dim, center, dim[2]/4, dim[2]/2)
-#magnitude = mc.create_mag_dist_vortex(mag_shape)
-#mag_data = MagData(a, magnitude)
-#mag_data.quiver_plot3d()
+center = (int(dim[0]/2)-0.5, int(dim[1]/2)-0.5, int(dim[2]/2)-0.5)
+mag_shape = mc.Shapes.disc(dim, center, dim[2]/4, dim[2]/2)
+magnitude = mc.create_mag_dist_vortex(mag_shape)
+mag_data = MagData(a, magnitude)
-vortex_shape = mc.Shapes.disc(dim, (3.5, 9.5, 9.5), 5, 4)
-sphere_shape = mc.Shapes.sphere(dim, (3.5, 22.5, 9.5), 3)
-slab_shape = mc.Shapes.slab(dim, (3.5, 15.5, 22.5), (5, 8, 8))
-mag_data = MagData(a, mc.create_mag_dist_vortex(vortex_shape, (3.5, 9.5, 9.5)))
-mag_data += MagData(a, mc.create_mag_dist_homog(sphere_shape, pi/4, pi/4))
-mag_data += MagData(a, mc.create_mag_dist_homog(slab_shape, -pi/6))
+#vortex_shape = mc.Shapes.disc(dim, (3.5, 9.5, 9.5), 5, 4)
+#sphere_shape = mc.Shapes.sphere(dim, (3.5, 22.5, 9.5), 3)
+#slab_shape = mc.Shapes.slab(dim, (3.5, 15.5, 22.5), (5, 8, 8))
+#mag_data = MagData(a, mc.create_mag_dist_vortex(vortex_shape, (3.5, 9.5, 9.5)))
+#mag_data += MagData(a, mc.create_mag_dist_homog(sphere_shape, pi/4, pi/4))
+#mag_data += MagData(a, mc.create_mag_dist_homog(slab_shape, -pi/6))
mag_data.quiver_plot3d()
@@ -70,28 +68,23 @@ projectors_xy_miss.extend(projectors_y_miss)
projectors = projectors_xy_full
noise = 0.0
###################################################################################################
-
-phasemappers = [PMConvolve(mag_data.a, projector, b_0) for projector in projectors]
-if noise != 0:
- phase_maps = [pm(mag_data) + PhaseMap(a, np.random.normal(0, noise, dim_uv))
- for pm in phasemappers]
-else:
- phase_maps = [pm(mag_data) for pm in phasemappers]
-
-###################################################################################################
print('--Setting up data collection')
-dim_uv = dim[1:3]
-data = DataSet(a, dim_uv, b_0)
-[data.append((phase_maps[i], projectors[i])) for i in range(len(projectors))]
+data = DataSet(a, dim, b_0)
+data.projectors = projectors
+data.phase_maps = data.create_phase_maps(mag_data)
+
+if noise != 0:
+ for phase_map in data.phase_maps:
+ phase_map += PhaseMap(a, np.random.normal(0, noise, dim_uv))
###################################################################################################
print('--Test simple solver')
-lam = 1E-10
-reg = ZeroOrderRegularisator(lam, 3*np.prod(dim))
+lam = 1E-4
+reg = ZeroOrderRegularisator(lam)
start = clock()
-mag_data_opt = rc.optimize_sparse_cg(data, regularisator=reg, verbosity=1)
+mag_data_opt = rc.optimize_sparse_cg(data, regularisator=reg, maxiter=100, verbosity=2)
print 'Time:', str(clock()-start)
mag_data_opt.quiver_plot3d()
(mag_data_opt - mag_data).quiver_plot3d()
@@ -102,5 +95,5 @@ print('--Plot stuff')
phase_maps_opt = data.create_phase_maps(mag_data_opt)
#data.display_phase()
#data.display_phase(phase_maps_opt)
-phase_diffs = [(phase_maps[i]-phase_maps_opt[i]) for i in range(len(data.data))]
-data.display_phase(phase_diffs)
+phase_diffs = [(data.phase_maps[i]-phase_maps_opt[i]) for i in range(len(data.phase_maps))]
+[phase_diff.display_phase() for phase_diff in phase_diffs]
diff --git a/scripts/test methods/compare_methods.py b/scripts/test methods/compare_methods.py
index fca81c9..2c07e5a 100644
--- a/scripts/test methods/compare_methods.py
+++ b/scripts/test methods/compare_methods.py
@@ -13,7 +13,8 @@ import pyramid
import pyramid.magcreator as mc
import pyramid.analytic as an
from pyramid.projector import SimpleProjector
-from pyramid.phasemapper import PMAdapterFM, PMReal, PMConvolve, PMFourier
+from pyramid.phasemapper import PhaseMapperRDRC, PhaseMapperRDFC, PhaseMapperFDFC
+from pyramid.kernel import Kernel
from pyramid.magdata import MagData
import logging
@@ -56,43 +57,36 @@ elif geometry == 'sphere':
# Create MagData object and projector:
mag_data = MagData(a, mc.create_mag_dist_homog(mag_shape, phi))
projector = SimpleProjector(dim)
+projection = projector(mag_data)
# Construct PhaseMapper objects:
-pm_adap = PMAdapterFM(a, projector, b_0)
-pm_real = PMReal(a, projector, b_0, numcore=numcore)
-pm_conv = PMConvolve(a, projector, b_0)
-pm_four = PMFourier(a, projector, b_0, padding=padding)
+pm_real = PhaseMapperRDRC(Kernel(a, projector.dim_uv, b_0), numcore=numcore)
+pm_conv = PhaseMapperRDFC(Kernel(a, projector.dim_uv, b_0))
+pm_four = PhaseMapperFDFC(a, projector.dim_uv, b_0, padding=padding)
# Get times for different approaches:
start_time = time.time()
-phase_map_adap = pm_adap(mag_data)
-print 'Time for PMAdapterFM: ', time.time() - start_time
+phase_map_real = pm_real(projection)
+print 'Time for RDRC: ', time.time() - start_time
start_time = time.time()
-phase_map_real = pm_real(mag_data)
-print 'Time for PMReal: ', time.time() - start_time
+phase_map_conv = pm_conv(projection)
+print 'Time for RDFC: ', time.time() - start_time
start_time = time.time()
-phase_map_conv = pm_conv(mag_data)
-print 'Time for PMConvolve: ', time.time() - start_time
-start_time = time.time()
-phase_map_four = pm_four(mag_data)
-print 'Time for PMFourier: ', time.time() - start_time
+phase_map_four = pm_four(projection)
+print 'Time for FDFC: ', time.time() - start_time
# Display the combinated plots with phasemap and holography image:
phase_map_ana.display_combined('Analytic Solution', gain=gain)
-phase_map_adap.display_combined('PMAdapterFM', gain=gain)
-phase_map_real.display_combined('PMReal', gain=gain)
-phase_map_conv.display_combined('PMConvolve', gain=gain)
-phase_map_four.display_combined('PMFourier', gain=gain)
+phase_map_real.display_combined('RDRC', gain=gain)
+phase_map_conv.display_combined('RDFC', gain=gain)
+phase_map_four.display_combined('FDFC', gain=gain)
# Plot differences to the analytic solution:
-phase_map_diff_adap = phase_map_adap - phase_map_ana
phase_map_diff_real = phase_map_real - phase_map_ana
phase_map_diff_conv = phase_map_conv - phase_map_ana
phase_map_diff_four = phase_map_four - phase_map_ana
-RMS_adap = np.std(phase_map_diff_adap.phase)
RMS_real = np.std(phase_map_diff_real.phase)
RMS_conv = np.std(phase_map_diff_conv.phase)
RMS_four = np.std(phase_map_diff_four.phase)
-phase_map_diff_adap.display_phase('PMAdapterFM difference (RMS = {:3.2e})'.format(RMS_adap))
-phase_map_diff_real.display_phase('PMReal difference (RMS = {:3.2e})'.format(RMS_real))
-phase_map_diff_conv.display_phase('PMConvolve difference (RMS = {:3.2e})'.format(RMS_conv))
-phase_map_diff_four.display_phase('PMFourier difference (RMS = {:3.2e})'.format(RMS_four))
+phase_map_diff_real.display_phase('RDRC difference (RMS = {:3.2e})'.format(RMS_real))
+phase_map_diff_conv.display_phase('RDFC difference (RMS = {:3.2e})'.format(RMS_conv))
+phase_map_diff_four.display_phase('FDFC difference (RMS = {:3.2e})'.format(RMS_four))
diff --git a/scripts/test methods/compare_pixel_fields.py b/scripts/test methods/compare_pixel_fields.py
index c26bb99..a8f3b8f 100644
--- a/scripts/test methods/compare_pixel_fields.py
+++ b/scripts/test methods/compare_pixel_fields.py
@@ -10,7 +10,8 @@ import numpy as np
import pyramid
import pyramid.magcreator as mc
from pyramid.projector import SimpleProjector
-from pyramid.phasemapper import PMConvolve, PMFourier
+from pyramid.phasemapper import PhaseMapperRDFC, PhaseMapperFDFC
+from pyramid.kernel import Kernel
from pyramid.magdata import MagData
from pyramid.phasemap import PhaseMap
@@ -56,15 +57,15 @@ 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')
projector = SimpleProjector(dim)
# Kernel of a disc in real space:
-phase_map_disc = PMConvolve(a, projector, geometry='disc')(mag_data)
+phase_map_disc = PhaseMapperRDFC(Kernel(a, projector.dim_uv, geometry='disc'))(mag_data)
phase_map_disc.unit = 'mrad'
phase_map_disc.display_phase('Phase of one Pixel (Disc)', limit=limit)
# Kernel of a slab in real space:
-phase_map_slab = PMConvolve(a, projector, geometry='slab')(mag_data)
+phase_map_slab = PhaseMapperRDFC(Kernel(a, projector.dim_uv, geometry='slab'))(mag_data)
phase_map_slab.unit = 'mrad'
phase_map_slab.display_phase('Phase of one Pixel (Slab)', limit=limit)
# Kernel of the Fourier method:
-phase_map_fft = PMFourier(a, projector, padding=0)(mag_data)
+phase_map_fft = PhaseMapperFDFC(a, projector.dim_uv, padding=0)(mag_data)
phase_map_fft.unit = 'mrad'
phase_map_fft.display_phase('Phase of one Pixel (Fourier)', limit=limit)
# Kernel of the Fourier method, calculated directly:
diff --git a/scripts/test methods/projector_test.py b/scripts/test methods/projector_test.py
index 616e69e..a82b92b 100644
--- a/scripts/test methods/projector_test.py
+++ b/scripts/test methods/projector_test.py
@@ -7,7 +7,8 @@ import numpy as np
import pyramid.magcreator as mc
from pyramid.magdata import MagData
from pyramid.projector import SimpleProjector, XTiltProjector, YTiltProjector
-from pyramid.phasemapper import PMConvolve
+from pyramid.phasemapper import PhaseMapperRDFC
+from pyramid.kernel import Kernel
dim = (16, 24, 32)
@@ -19,11 +20,12 @@ b_0 = 1.5
shape = mc.Shapes.slab(dim, center, width)
mag_data = MagData(a, mc.create_mag_dist_homog(shape, phi=np.pi/6, theta=np.pi/3))
+phase_mapper = PhaseMapperRDFC(Kernel(a, dim_uv, b_0))
# PROJECTOR TESTING
-PMConvolve(a, SimpleProjector(dim, 'z', dim_uv), b_0)(mag_data).display_phase('z simple')
-PMConvolve(a, XTiltProjector(dim, 0, dim_uv), b_0)(mag_data).display_phase('z (x-tilt)')
-PMConvolve(a, SimpleProjector(dim, 'x', dim_uv), b_0)(mag_data).display_phase('x simple')
-PMConvolve(a, YTiltProjector(dim, np.pi/2, dim_uv), b_0)(mag_data).display_phase('x (y-tilt)')
-PMConvolve(a, XTiltProjector(dim, np.pi/2, dim_uv), b_0)(mag_data).display_phase('y (x-tilt)')
-PMConvolve(a, SimpleProjector(dim, 'y', dim_uv), b_0)(mag_data).display_phase('y simple')
+phase_mapper(SimpleProjector(dim, 'z', dim_uv)(mag_data)).display_phase('z simple')
+phase_mapper(XTiltProjector(dim, 0, dim_uv)(mag_data)).display_phase('z (x-tilt)')
+phase_mapper(SimpleProjector(dim, 'x', dim_uv)(mag_data)).display_phase('x simple')
+phase_mapper(YTiltProjector(dim, np.pi/2, dim_uv)(mag_data)).display_phase('x (y-tilt)')
+phase_mapper(XTiltProjector(dim, np.pi/2, dim_uv)(mag_data)).display_phase('y (x-tilt)')
+phase_mapper(SimpleProjector(dim, 'y', dim_uv)(mag_data)).display_phase('y simple')
--
GitLab