From 4ab7bc635c4c79c7c9cf89040f45ce3ab65d40fc Mon Sep 17 00:00:00 2001
From: Jan Caron <j.caron@fz-juelich.de>
Date: Mon, 24 Nov 2014 20:12:50 +0100
Subject: [PATCH] Big Bunch of small updates! setup: Changes to setup.py,
pyramid.version to fit to numpy and jutil now uses __all__ magic and
supports from pyramid import * tests: changes to
test_compliance/magcreator/analytic (should now work) collaborations: added
several scripts for Patrick reconstruction: added batch file for the
reconstruction of several configs test methods: (unfinished) scripts for fftw
tests and kernel comparisons paper 1: ch5-4 changed padding values to 0, 1, 3
and 7 plotting: deleted plt.show() lines (have to be done manually now or via
Spyder) logging: now named _log instead of LOG (private variable instead of
constant) also some logs were disabled, because Pylint says they are
slow... reconstruction: now has info argument (pass a list in which cost info
is copied) phasemapper: changed sign of constant PHI_0, is now correctly
positive (instead the negative z-integration is taken into
account which was not done before!) kernel: now uses FFTW, numpy
code was changed accordingly and is still available phasemapper: also uses
FFTW now, work in progress: FFTW for jac_T_dot
---
pyramid/__init__.py | 36 ++-
pyramid/_version.py | 4 -
pyramid/analytic.py | 22 +-
pyramid/converter.py | 68 -----
pyramid/costfunction.py | 34 +--
pyramid/dataset.py | 23 +-
pyramid/forwardmodel.py | 21 +-
pyramid/kernel.py | 115 +++++---
pyramid/magcreator.py | 23 +-
pyramid/magdata.py | 87 +++---
pyramid/phasemap.py | 102 +++----
pyramid/phasemapper.py | 271 +++++++++++-------
pyramid/projector.py | 70 ++---
pyramid/reconstruction.py | 36 +--
pyramid/regularisator.py | 43 +--
pyramid/version.py | 3 +
scripts/collaborations/bielefeld.py | 71 ++---
scripts/collaborations/example_joern.py | 35 ++-
.../patrick_nanowire_reconstruction.py | 56 ++++
scripts/collaborations/test.py | 29 --
.../vtk_conversion_nanowire_full_90deg.py | 180 ++++++++++++
scripts/gui/mag_slicer.py | 7 +-
.../ch5-4-comparison_of_methods_new.py | 75 +++--
..._test.py => reconstruction_linear_test.py} | 67 ++---
.../reconstruction_linear_test_batch.py | 148 ++++++++++
...econstruction_sparse_cg_singular_values.py | 3 +-
scripts/test methods/compare_kernels.py | 121 ++++++++
scripts/test methods/compare_methods.py | 11 +-
scripts/test methods/fftw_test.py | 55 ++++
setup.py | 140 +++++----
tests/test_analytic.py | 39 ++-
tests/test_compliance.py | 6 +-
tests/test_magcreator.py | 51 ++--
tests/test_magcreator/ref_mag_disc.npy | Bin 2480 -> 2480 bytes
tests/test_magcreator/ref_mag_vort.npy | Bin 2480 -> 2480 bytes
35 files changed, 1344 insertions(+), 708 deletions(-)
delete mode 100644 pyramid/_version.py
delete mode 100644 pyramid/converter.py
create mode 100644 pyramid/version.py
create mode 100644 scripts/collaborations/patrick_nanowire_reconstruction.py
delete mode 100644 scripts/collaborations/test.py
create mode 100644 scripts/collaborations/vtk_conversion_nanowire_full_90deg.py
rename scripts/reconstruction/{reconstruction_sparse_cg_test.py => reconstruction_linear_test.py} (59%)
create mode 100644 scripts/reconstruction/reconstruction_linear_test_batch.py
create mode 100644 scripts/test methods/compare_kernels.py
create mode 100644 scripts/test methods/fftw_test.py
diff --git a/pyramid/__init__.py b/pyramid/__init__.py
index fc3426a..aed6a3d 100644
--- a/pyramid/__init__.py
+++ b/pyramid/__init__.py
@@ -1,4 +1,7 @@
# -*- coding: utf-8 -*-
+# Copyright 2014 by Forschungszentrum Juelich GmbH
+# Author: J. Caron
+#
"""Package for the creation and reconstruction of magnetic distributions and resulting phase maps.
Modules
@@ -33,4 +36,35 @@ numcore
"""
-from _version import __version__
+
+import logging
+
+from . import analytic
+from . import magcreator
+from . import reconstruction
+from .costfunction import *
+from .dataset import *
+from .forwardmodel import *
+from .kernel import *
+from .magdata import *
+from .phasemap import *
+from .phasemapper import *
+from .projector import *
+from .regularisator import *
+from .version import version as __version__
+from .version import hg_revision as __hg_revision__
+
+_log = logging.getLogger(__name__)
+_log.info("Starting PYRAMID V{} HG{}".format(__version__, __hg_revision__))
+del logging
+
+__all__ = ['__version__', '__hg_revision__', 'analytic', 'magcreator', 'reconstruction']
+__all__.extend(costfunction.__all__)
+__all__.extend(dataset.__all__)
+__all__.extend(forwardmodel.__all__)
+__all__.extend(kernel.__all__)
+__all__.extend(magdata.__all__)
+__all__.extend(phasemap.__all__)
+__all__.extend(phasemapper.__all__)
+__all__.extend(projector.__all__)
+__all__.extend(regularisator.__all__)
diff --git a/pyramid/_version.py b/pyramid/_version.py
deleted file mode 100644
index 1925341..0000000
--- a/pyramid/_version.py
+++ /dev/null
@@ -1,4 +0,0 @@
-# -*- coding: utf-8 -*-
-"""Version string"""
-
-__version__ = '0.1.0'
diff --git a/pyramid/analytic.py b/pyramid/analytic.py
index a222a22..5faef31 100644
--- a/pyramid/analytic.py
+++ b/pyramid/analytic.py
@@ -16,8 +16,10 @@ from pyramid.phasemap import PhaseMap
import logging
-LOG = logging.getLogger(__name__)
-PHI_0 = -2067.83 # magnetic flux in T*nm²
+__all__ = ['phase_mag_slab', 'phase_mag_slab', 'phase_mag_sphere', 'phase_mag_vortex']
+_log = logging.getLogger(__name__)
+
+PHI_0 = 2067.83 # magnetic flux in T*nm²
def phase_mag_slab(dim, a, phi, center, width, b_0=1):
@@ -45,7 +47,7 @@ def phase_mag_slab(dim, a, phi, center, width, b_0=1):
The phase as a 2-dimensional array.
'''
- LOG.debug('Calling phase_mag_slab')
+ _log.debug('Calling phase_mag_slab')
# Function for the phase:
def phi_mag(x, y):
@@ -66,7 +68,7 @@ def phase_mag_slab(dim, a, phi, center, width, b_0=1):
y0 = a * (center[1] + 0.5) # y0, x0 define the center of a pixel,
x0 = a * (center[2] + 0.5) # hence: (cellindex + 0.5) * grid spacing
Lz, Ly, Lx = a * width[0], a * width[1], a * width[2]
- coeff = b_0 / (4*PHI_0)
+ coeff = - b_0 / (4*PHI_0) # Minus because of negative z-direction
# Create grid:
x = np.linspace(a/2, x_dim*a-a/2, num=x_dim)
y = np.linspace(a/2, y_dim*a-a/2, num=y_dim)
@@ -102,7 +104,7 @@ def phase_mag_disc(dim, a, phi, center, radius, height, b_0=1):
The phase as a 2-dimensional array.
'''
- LOG.debug('Calling phase_mag_disc')
+ _log.debug('Calling phase_mag_disc')
# Function for the phase:
def phi_mag(x, y):
@@ -118,7 +120,7 @@ def phase_mag_disc(dim, a, phi, center, radius, height, b_0=1):
x0 = a * (center[2] + 0.5) # hence: cellindex + 0.5
Lz = a * height
R = a * radius
- coeff = - pi * b_0 / (2*PHI_0)
+ coeff = pi * b_0 / (2*PHI_0) # Minus is gone because of negative z-direction
# Create grid:
x = np.linspace(a/2, x_dim*a-a/2, num=x_dim)
y = np.linspace(a/2, y_dim*a-a/2, num=y_dim)
@@ -152,7 +154,7 @@ def phase_mag_sphere(dim, a, phi, center, radius, b_0=1):
The phase as a 2-dimensional array.
'''
- LOG.debug('Calling phase_mag_sphere')
+ _log.debug('Calling phase_mag_sphere')
# Function for the phase:
def phi_mag(x, y):
@@ -167,7 +169,7 @@ def phase_mag_sphere(dim, a, phi, center, radius, b_0=1):
y0 = a * (center[1] + 0.5) # y0, x0 have to be in the center of a pixel,
x0 = a * (center[2] + 0.5) # hence: cellindex + 0.5
R = a * radius
- coeff = - 2./3. * pi * b_0 / PHI_0
+ coeff = 2./3. * pi * b_0 / PHI_0 # Minus is gone because of negative z-direction
# Create grid:
x = np.linspace(a / 2, x_dim * a - a / 2, num=x_dim)
y = np.linspace(a / 2, y_dim * a - a / 2, num=y_dim)
@@ -201,7 +203,7 @@ def phase_mag_vortex(dim, a, center, radius, height, b_0=1):
The phase as a 2-dimensional array.
'''
- LOG.debug('Calling phase_mag_vortex')
+ _log.debug('Calling phase_mag_vortex')
# Function for the phase:
def phi_mag(x, y):
@@ -215,7 +217,7 @@ def phase_mag_vortex(dim, a, center, radius, height, b_0=1):
x0 = a * (center[2] + 0.5) # hence: cellindex + 0.5
Lz = a * height
R = a * radius
- coeff = pi * b_0 * Lz / PHI_0
+ coeff = - pi * b_0 * Lz / PHI_0 # Minus because of negative z-direction
# Create grid:
x = np.linspace(a/2, x_dim*a-a/2, num=x_dim)
y = np.linspace(a/2, y_dim*a-a/2, num=y_dim)
diff --git a/pyramid/converter.py b/pyramid/converter.py
deleted file mode 100644
index 3303227..0000000
--- a/pyramid/converter.py
+++ /dev/null
@@ -1,68 +0,0 @@
-# -*- coding: utf-8 -*-
-"""
-Created on Tue Aug 19 08:48:45 2014
-
-@author: Jan
-""" # TODO: Docstring
-
-# TODO: put into other class
-# TODO: use 3 components (more complex)
-# TODO: take masks into account
-
-
-import numpy as np
-
-
-class IndexConverter3components(object):
-
- def __init__(self, dim):
- self.dim = dim
- self.size_3d = np.prod(dim)
- self.size_2d = dim[2]*dim[1]
-
- def ind_to_coord(self, ind):
- m, remain = int(ind/self.size_3d), ind % self.size_3d
- z, remain = int(remain/self.size_2d), remain%self.size_2d
- y, x = int(remain/self.dim[2]), remain%self.dim[2]
- coord = m, z, y, x
- return coord
-
- def coord_to_ind(self, coord):
- z, y, x = coord
- ind = [i*self.size_3d + z*self.size_2d + y*self.dim[2] + x for i in range(3)]
- return ind
-
- def find_neighbour_ind(self, coord):
- z, y, x = coord
- t, d = (z-1, y, x), (z+1, y, x) # t: top, d: down
- f, b = (z, y-1, x), (z, y+1, x) # f: front, b: back
- l, r = (z, y, x-1), (z, y, x+1) # l: left, r: right
- neighbours = [self.coord_to_ind(i) for i in [t, d, f, b, l, r]]
- return np.reshape(np.swapaxes(neighbours, 0, 1), (3, 3, 2))
-
-
-class IndexConverter(object):
-
- def __init__(self, dim):
- self.dim = dim
- self.size_2d = dim[2]*dim[1]
-
- def ind_to_coord(self, ind):
- z, remain = int(ind/self.size_2d), ind%self.size_2d
- y, x = int(remain/self.dim[2]), remain%self.dim[2]
- coord = z, y, x
- return coord
-
- def coord_to_ind(self, coord):
- z, y, x = coord
- ind = z*self.size_2d + y*self.dim[2] + x
- return ind
-
- def find_neighbour_ind(self, coord):
- z, y, x = coord
- t, d = (z-1, y, x), (z+1, y, x) # t: top, d: down
- f, b = (z, y-1, x), (z, y+1, x) # f: front, b: back
- l, r = (z, y, x-1), (z, y, x+1) # l: left, r: right
- neighbours = [self.coord_to_ind(i) for i in [t, d, f, b, l, r]]
- return neighbours
- return np.reshape(neighbours, (3, 2))
diff --git a/pyramid/costfunction.py b/pyramid/costfunction.py
index efa7500..0859ae3 100644
--- a/pyramid/costfunction.py
+++ b/pyramid/costfunction.py
@@ -6,14 +6,16 @@ the so called `cost` of a threedimensional magnetization distribution."""
import numpy as np
from scipy.sparse.linalg import LinearOperator
-from scipy.sparse import eye
from pyramid.forwardmodel import ForwardModel
-from pyramid.regularisator import ZeroOrderRegularisator
from pyramid.regularisator import NoneRegularisator
+
import logging
+__all__ = ['Costfunction']
+
+
class Costfunction(object):
'''Class for calculating the cost of a 3D magnetic distributions in relation to 2D phase maps.
@@ -47,10 +49,10 @@ class Costfunction(object):
'''
- LOG = logging.getLogger(__name__+'.Costfunction')
+ _log = logging.getLogger(__name__+'.Costfunction')
def __init__(self, data_set, regularisator):
- self.LOG.debug('Calling __init__')
+ self._log.debug('Calling __init__')
self.data_set = data_set
self.fwd_model = ForwardModel(data_set)
self.regularisator = regularisator
@@ -61,29 +63,29 @@ class Costfunction(object):
self.Se_inv = data_set.Se_inv
self.n = data_set.n
self.m = data_set.m
- self.LOG.debug('Created '+str(self))
+ self._log.debug('Created '+str(self))
def __repr__(self):
- self.LOG.debug('Calling __repr__')
+ self._log.debug('Calling __repr__')
return '%s(data_set=%r, fwd_model=%r, regularisator=%r)' % \
(self.__class__, self.data_set, self.fwd_model, self.regularisator)
def __str__(self):
- self.LOG.debug('Calling __str__')
+ self._log.debug('Calling __str__')
return 'Costfunction(data_set=%s, fwd_model=%s, regularisator=%s)' % \
(self.data_set, self.fwd_model, self.regularisator)
- def init(self, x):
- self(x)
-
def __call__(self, x):
- self.LOG.debug('Calling __call__')
+ self._log.debug('Calling __call__')
delta_y = self.fwd_model(x) - self.y
self.chisq_m = delta_y.dot(self.Se_inv.dot(delta_y))
self.chisq_a = self.regularisator(x)
self.chisq = self.chisq_m + self.chisq_a
return self.chisq
+ def init(self, x):
+ self(x)
+
def jac(self, x):
'''Calculate the derivative of the costfunction for a given magnetization distribution.
@@ -98,7 +100,7 @@ class Costfunction(object):
Jacobi vector which represents the cost derivative of all voxels of the magnetization.
'''
- self.LOG.debug('Calling jac')
+ self._log.debug('Calling jac')
assert len(x) == self.n
return (2 * self.fwd_model.jac_T_dot(x, self.Se_inv.dot(self.fwd_model(x) - self.y))
+ self.regularisator.jac(x))
@@ -122,7 +124,7 @@ class Costfunction(object):
Product of the input `vector` with the Hessian matrix of the costfunction.
'''
- self.LOG.debug('Calling hess_dot')
+# self._log.debug('Calling hess_dot') # TODO: Profiler says this was slow...
return (2 * self.fwd_model.jac_T_dot(x, self.Se_inv.dot(self.fwd_model.jac_dot(x, vector)))
+ self.regularisator.hess_dot(x, vector))
@@ -148,10 +150,10 @@ class CFAdapterScipyCG(LinearOperator):
'''
# TODO: make obsolete!
- LOG = logging.getLogger(__name__+'.CFAdapterScipyCG')
+ _log = logging.getLogger(__name__+'.CFAdapterScipyCG')
def __init__(self, cost):
- self.LOG.debug('Calling __init__')
+ self._log.debug('Calling __init__')
self.cost = cost
def matvec(self, vector):
@@ -164,7 +166,7 @@ class CFAdapterScipyCG(LinearOperator):
costfunction.
'''
- self.LOG.debug('Calling matvec')
+ self._log.debug('Calling matvec')
return self.cost.hess_dot(None, vector)
@property
diff --git a/pyramid/dataset.py b/pyramid/dataset.py
index 62db474..9559d59 100644
--- a/pyramid/dataset.py
+++ b/pyramid/dataset.py
@@ -6,7 +6,7 @@ and additional data like corresponding projectors."""
import numpy as np
from numbers import Number
-import scipy.sparse as sp
+from scipy.sparse import eye as sparse_eye
import matplotlib.pyplot as plt
@@ -18,6 +18,9 @@ from pyramid.kernel import Kernel
import logging
+__all__ = ['DataSet']
+
+
class DataSet(object):
'''Class for collecting phase maps and corresponding projectors.
@@ -53,7 +56,7 @@ class DataSet(object):
'''
- LOG = logging.getLogger(__name__+'.DataSet')
+ _log = logging.getLogger(__name__+'.DataSet')
@property
def m(self):
@@ -62,7 +65,7 @@ class DataSet(object):
@property
def Se_inv(self):
# TODO: better implementation, maybe get-method? more flexible? input in append?
- return sp.eye(self.m)
+ return sparse_eye(self.m)
@property
def phase_vec(self):
@@ -82,7 +85,7 @@ class DataSet(object):
return {kernel.dim_uv: PhaseMapperRDFC(kernel) for kernel in kernel_list}
def __init__(self, a, dim, b_0=1, mask=None):
- self.LOG.debug('Calling __init__')
+ self._log.debug('Calling __init__')
assert isinstance(a, Number), 'Grid spacing has to be a number!'
assert a >= 0, 'Grid spacing has to be a positive number!'
assert isinstance(dim, tuple) and len(dim) == 3, \
@@ -98,14 +101,14 @@ class DataSet(object):
self.mask = mask
self.phase_maps = []
self.projectors = []
- self.LOG.debug('Created: '+str(self))
+ self._log.debug('Created: '+str(self))
def __repr__(self):
- self.LOG.debug('Calling __repr__')
+ self._log.debug('Calling __repr__')
return '%s(a=%r, dim=%r, b_0=%r)' % (self.__class__, self.a, self.dim, self.b_0)
def __str__(self):
- self.LOG.debug('Calling __str__')
+ self._log.debug('Calling __str__')
return 'DataSet(a=%s, dim=%s, b_0=%s)' % (self.a, self.dim, self.b_0)
def append(self, phase_map, projector): # TODO: include Se_inv or 2D mask??
@@ -123,7 +126,7 @@ class DataSet(object):
None
'''
- self.LOG.debug('Calling append')
+ self._log.debug('Calling append')
assert isinstance(phase_map, PhaseMap) and isinstance(projector, Projector), \
'Argument has to be a tuple of a PhaseMap and a Projector object!'
assert projector.dim == self.dim, '3D dimensions must match!'
@@ -176,7 +179,7 @@ class DataSet(object):
None
'''
- self.LOG.debug('Calling display')
+ self._log.debug('Calling display')
if mag_data is not None:
phase_maps = self.create_phase_maps(mag_data)
else:
@@ -223,7 +226,7 @@ class DataSet(object):
None
'''
- self.LOG.debug('Calling display_combined')
+ self._log.debug('Calling display_combined')
if mag_data is not None:
phase_maps = self.create_phase_maps(mag_data)
else:
diff --git a/pyramid/forwardmodel.py b/pyramid/forwardmodel.py
index 1cd1af7..6a61a17 100644
--- a/pyramid/forwardmodel.py
+++ b/pyramid/forwardmodel.py
@@ -9,6 +9,9 @@ from pyramid.magdata import MagData
import logging
+__all__ = ['ForwardModel']
+
+
class ForwardModel(object):
'''Class for mapping 3D magnetic distributions to 2D phase maps.
@@ -37,28 +40,29 @@ class ForwardModel(object):
'''
- LOG = logging.getLogger(__name__+'.ForwardModel')
+ _log = logging.getLogger(__name__+'.ForwardModel')
def __init__(self, data_set):
- self.LOG.debug('Calling __init__')
+ self._log.debug('Calling __init__')
self.data_set = data_set
self.phase_mappers = data_set.phase_mappers
self.m = data_set.m
self.n = data_set.n
+ self.shape = (self.m, self.n)
self.hook_points = data_set.hook_points
self.mag_data = MagData(data_set.a, np.zeros((3,)+data_set.dim))
- self.LOG.debug('Creating '+str(self))
+ self._log.debug('Creating '+str(self))
def __repr__(self):
- self.LOG.debug('Calling __repr__')
+ self._log.debug('Calling __repr__')
return '%s(data_set=%r)' % (self.__class__, self.data_set)
def __str__(self):
- self.LOG.debug('Calling __str__')
+ self._log.debug('Calling __str__')
return 'ForwardModel(data_set=%s)' % (self.data_set)
def __call__(self, x):
- self.LOG.debug('Calling __call__')
+ self._log.debug('Calling __call__')
self.mag_data.magnitude[:] = 0
self.mag_data.set_vector(x, self.data_set.mask)
# TODO: Multiprocessing
@@ -89,7 +93,7 @@ class ForwardModel(object):
`vector`.
'''
- self.LOG.debug('Calling jac_dot')
+# self._log.debug('Calling jac_dot') # TODO: Profiler says this was slow...
self.mag_data.magnitude[:] = 0
self.mag_data.set_vector(vector, self.data_set.mask)
result = np.zeros(self.m)
@@ -119,8 +123,7 @@ class ForwardModel(object):
the input `vector`.
'''
- self.LOG.debug('Calling jac_T_dot')
-
+# self._log.debug('Calling jac_T_dot') # TODO: Profiler says this was slow...
result = np.zeros(3*np.prod(self.data_set.dim))
hp = self.hook_points
for i, projector in enumerate(self.data_set.projectors):
diff --git a/pyramid/kernel.py b/pyramid/kernel.py
index 63f1127..cb237b0 100644
--- a/pyramid/kernel.py
+++ b/pyramid/kernel.py
@@ -8,7 +8,9 @@ import numpy as np
import logging
-PHI_0 = -2067.83 # magnetic flux in T*nm²
+__all__ = ['Kernel']
+
+PHI_0 = 2067.83 # magnetic flux in T*nm²
class Kernel(object):
@@ -32,9 +34,6 @@ class Kernel(object):
be calculated.
b_0 : float, optional
Saturation magnetization in Tesla, which is used for the phase calculation. Default is 1.
- numcore : boolean, optional
- Boolean choosing if Cython enhanced routines from the :mod:`~.pyramid.numcore` module
- should be used. Default is True.
geometry : {'disc', 'slab'}, optional
The elementary geometry of the single magnetized pixel.
u : :class:`~numpy.ndarray` (N=3)
@@ -54,57 +53,89 @@ class Kernel(object):
A tuple of :class:`slice` objects to extract the original field of view from the increased
size (`size_fft`) of the grid for the FFT-convolution.
- '''
-
- LOG = logging.getLogger(__name__+'.Kernel')
-
- def __init__(self, a, dim_uv, b_0=1., numcore=True, geometry='disc'):
- self.LOG.debug('Calling __init__')
+ ''' # TODO: overview what all dim_??? mean! and use_fftw
- # Function for the phase of an elementary geometry:
- def get_elementary_phase(geometry, n, m, a):
- if geometry == 'disc':
- in_or_out = np.logical_not(np.logical_and(n == 0, m == 0))
- return m / (n**2 + m**2 + 1E-30) * in_or_out
- elif geometry == 'slab':
- def F_a(n, m):
- A = np.log(a**2 * (n**2 + m**2))
- 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))
+ _log = logging.getLogger(__name__+'.Kernel')
+ def __init__(self, a, dim_uv, b_0=1., geometry='disc', use_fftw=True, threads=1):
+ self._log.debug('Calling __init__')
# Set basic properties:
- self.dim_uv = dim_uv # Size of the FOV, not the kernel (kernel is bigger)!
- self.size = np.prod(dim_uv) # Pixel count
+ self.dim_uv = dim_uv # Dimensions of the FOV
+ self.dim_kern = tuple(2*np.array(dim_uv)-1) # Dimensions of the kernel
+# self.size = np.prod(dim_uv) # TODO: is this even used? (Pixel count)
self.a = a
- self.numcore = numcore
self.geometry = geometry
+ # Set up FFT:
+ if use_fftw:
+ try:
+ import pyfftw
+ except ImportError:
+ use_fftw = False
+ self._log.info('pyFFTW could not be imported, using numpy instead!')
+ if use_fftw: # use pyfftw (FFTW wrapper for python)
+ self.dim_pad = tuple(2*np.array(dim_uv)) # is at least even (not nec. power of 2)
+ self.dim_fft = (self.dim_pad[0], self.dim_pad[1]/2.+1) # last axis is real
+ n = pyfftw.simd_alignment
+ self.u = pyfftw.n_byte_align_empty(self.dim_kern, n, dtype='float32')
+ self.v = pyfftw.n_byte_align_empty(self.dim_kern, n, dtype='float32')
+ self.u_fft = pyfftw.n_byte_align_empty(self.dim_fft, n, dtype='complex64')
+ self.v_fft = pyfftw.n_byte_align_empty(self.dim_fft, n, dtype='complex64')
+ rfftn = pyfftw.builders.rfftn(self.u, self.dim_pad, threads=threads)
+ self.threads = threads
+ self.use_fftw = True
+ else: # otherwise use numpy
+ self.dim_pad = tuple(2**np.ceil(np.log2(2*np.array(dim_uv))).astype(int)) # pow(2)
+ self.dim_fft = (self.dim_pad[0], self.dim_pad[1]/2+1) # last axis is real
+ self.u = np.empty(self.dim_kern, dtype=float)
+ self.v = np.empty(self.dim_kern, dtype=float)
+ self.u_fft = np.empty(self.dim_fft, dtype=complex)
+ self.v_fft = np.empty(self.dim_fft, dtype=complex)
+ rfftn = lambda x: np.fft.rfftn(x, self.dim_pad)
+ self.use_fftw = False
# Calculate kernel (single pixel phase):
- coeff = -b_0 * a**2 / (2*PHI_0)
+ coeff = b_0 * a**2 / (2*PHI_0) # Minus is gone because of negative z-direction
v_dim, u_dim = dim_uv
u = np.linspace(-(u_dim-1), u_dim-1, num=2*u_dim-1)
v = np.linspace(-(v_dim-1), v_dim-1, num=2*v_dim-1)
uu, vv = np.meshgrid(u, v)
- self.u = coeff * get_elementary_phase(geometry, uu, vv, a)
- self.v = coeff * get_elementary_phase(geometry, vv, uu, a)
+ self.u[...] = coeff * self._get_elementary_phase(geometry, uu, vv, a)
+ self.v[...] = coeff * self._get_elementary_phase(geometry, vv, uu, a)
# Calculate Fourier trafo of kernel components:
- dim_combined = 3*np.array(dim_uv) - 2 # (dim_uv-1) + dim_uv + (dim_uv-1) mag + kernel
- self.dim_fft = 2 ** np.ceil(np.log2(dim_combined)).astype(int) # next multiple of 2
self.slice_fft = (slice(dim_uv[0]-1, 2*dim_uv[0]-1), slice(dim_uv[1]-1, 2*dim_uv[1]-1))
- self.slice_fft_compl = (slice(2*dim_uv[0]-1, 2*dim_uv[0]-1), slice(2*dim_uv[1]-1, 2*dim_uv[1]-1))
- self.u_fft = np.fft.rfftn(self.u, self.dim_fft)
- self.v_fft = np.fft.rfftn(self.v, self.dim_fft)
- self.u_fft_compl = np.fft.fftn(self.u, self.dim_fft)
- self.v_fft_compl = np.fft.fftn(self.v, self.dim_fft)
- self.LOG.debug('Created '+str(self))
+ self.u_fft[...] = rfftn(self.u)
+ self.v_fft[...] = rfftn(self.v)
+ self._log.debug('Created '+str(self))
+
+ # TODO: make pyfftw optional (SLOW if kernel has to be build every time like in pm()!)
+ # TODO: test if prior build of kernel brings speed up in test_method() or test_fftw()
+ # TODO: implement fftw also in phasemapper (JUST there, here: FFT TWICE and big overhead)
+ # TODO: BUT allocation of u/v/u_fft/v_fft could be beneficial (try useing with numpy.fft)
+ # TODO: Set plan manually? Save computation time also for kernel?
+ # TODO: Multithreading?
+ # TODO: TakeTime multiple runs?
def __repr__(self):
- self.LOG.debug('Calling __repr__')
- return '%s(a=%r, dim_uv=%r, numcore=%r, geometry=%r)' % \
- (self.__class__, self.a, self.dim_uv, self.numcore, self.geometry)
+ self._log.debug('Calling __repr__')
+ return '%s(a=%r, dim_uv=%r, geometry=%r)' % \
+ (self.__class__, self.a, self.dim_uv, self.geometry)
def __str__(self):
- self.LOG.debug('Calling __str__')
- return 'Kernel(a=%s, dim_uv=%s, numcore=%s, geometry=%s)' % \
- (self.a, self.dim_uv, self.numcore, self.geometry)
+ self._log.debug('Calling __str__')
+ return 'Kernel(a=%s, dim_uv=%s, geometry=%s)' % \
+ (self.a, self.dim_uv, self.geometry)
+
+ def _get_elementary_phase(self, geometry, n, m, a):
+ # TODO: Docstring! Function for the phase of an elementary geometry:
+ if geometry == 'disc':
+ in_or_out = np.logical_not(np.logical_and(n == 0, m == 0))
+ return m / (n**2 + m**2 + 1E-30) * in_or_out
+ elif geometry == 'slab':
+ def F_a(n, m):
+ A = np.log(a**2 * (n**2 + m**2))
+ 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))
+
+ def get_info(self):
+ pass
diff --git a/pyramid/magcreator.py b/pyramid/magcreator.py
index 299bb77..c4286d3 100644
--- a/pyramid/magcreator.py
+++ b/pyramid/magcreator.py
@@ -21,7 +21,8 @@ import abc
import logging
-LOG = logging.getLogger(__name__)
+__all__ = ['Shapes', 'create_mag_dist_homog', 'create_mag_dist_vortex']
+_log = logging.getLogger(__name__)
class Shapes(object):
@@ -36,7 +37,7 @@ class Shapes(object):
'''
__metaclass__ = abc.ABCMeta
- LOG = logging.getLogger(__name__+'.Shapes')
+ _log = logging.getLogger(__name__+'.Shapes')
@classmethod
def slab(cls, dim, center, width):
@@ -57,7 +58,7 @@ class Shapes(object):
The magnetic shape as a 3D-array with values between 1 and 0.
'''
- cls.LOG.debug('Calling slab')
+ cls._log.debug('Calling slab')
assert np.shape(dim) == (3,), 'Parameter dim has to be a tuple of length 3!'
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!'
@@ -92,7 +93,7 @@ class Shapes(object):
The magnetic shape as a 3D-array with values between 1 and 0.
'''
- cls.LOG.debug('Calling disc')
+ cls._log.debug('Calling disc')
assert np.shape(dim) == (3,), 'Parameter dim has to be a a tuple of length 3!'
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!'
@@ -141,7 +142,7 @@ class Shapes(object):
The magnetic shape as a 3D-array with values between 1 and 0.
'''
- cls.LOG.debug('Calling ellipse')
+ cls._log.debug('Calling ellipse')
assert np.shape(dim) == (3,), 'Parameter dim has to be a a tuple of length 3!'
assert np.shape(center) == (3,), 'Parameter center has to be a a tuple of length 3!'
assert np.shape(width) == (2,), 'Parameter width has to be a a tuple of length 2!'
@@ -189,7 +190,7 @@ class Shapes(object):
The magnetic shape as a 3D-array with values between 1 and 0.
'''
- cls.LOG.debug('Calling sphere')
+ cls._log.debug('Calling sphere')
assert np.shape(dim) == (3,), 'Parameter dim has to be a a tuple of length 3!'
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!'
@@ -220,7 +221,7 @@ class Shapes(object):
The magnetic shape as a 3D-array with values between 1 and 0.
'''
- cls.LOG.debug('Calling ellipsoid')
+ cls._log.debug('Calling ellipsoid')
assert np.shape(dim) == (3,), 'Parameter dim has to be a a tuple of length 3!'
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!'
@@ -253,7 +254,7 @@ class Shapes(object):
The magnetic shape as a 3D-array with values between 1 and 0.
'''
- cls.LOG.debug('Calling filament')
+ cls._log.debug('Calling filament')
assert np.shape(dim) == (3,), 'Parameter dim has to be a tuple of length 3!'
assert np.shape(pos) == (2,), 'Parameter pos has to be a tuple of length 2!'
assert axis in {'z', 'y', 'x'}, 'Axis has to be x, y or z (as a string)!'
@@ -283,7 +284,7 @@ class Shapes(object):
The magnetic shape as a 3D-array with values between 1 and 0.
'''
- cls.LOG.debug('Calling pixel')
+ cls._log.debug('Calling pixel')
assert np.shape(dim) == (3,), 'Parameter dim has to be a tuple of length 3!'
assert np.shape(pixel) == (3,), 'Parameter pixel has to be a tuple of length 3!'
mag_shape = np.zeros(dim)
@@ -313,7 +314,7 @@ def create_mag_dist_homog(mag_shape, phi, theta=pi/2, magnitude=1):
`x`-, `y`- and `z`-direction on the 3-dimensional grid.
'''
- LOG.debug('Calling create_mag_dist_homog')
+ _log.debug('Calling create_mag_dist_homog')
dim = np.shape(mag_shape)
assert len(dim) == 3, 'Magnetic shapes must describe 3-dimensional distributions!'
z_mag = np.ones(dim) * np.cos(theta) * mag_shape * magnitude
@@ -345,7 +346,7 @@ def create_mag_dist_vortex(mag_shape, center=None, axis='z', magnitude=1):
`x`-, `y`- and `z`-direction on the 3-dimensional grid.
'''
- LOG.debug('Calling create_mag_dist_vortex')
+ _log.debug('Calling create_mag_dist_vortex')
dim = np.shape(mag_shape)
assert len(dim) == 3, 'Magnetic shapes must describe 3-dimensional distributions!'
assert axis in {'z', 'y', 'x'}, 'Axis has to be x, y or z (as a string)!'
diff --git a/pyramid/magdata.py b/pyramid/magdata.py
index b2460ea..0c8e2fc 100644
--- a/pyramid/magdata.py
+++ b/pyramid/magdata.py
@@ -19,6 +19,9 @@ import netCDF4
import logging
+__all__ = ['MagData']
+
+
class MagData(object):
'''Class for storing magnetization data.
@@ -45,7 +48,7 @@ class MagData(object):
'''
- LOG = logging.getLogger(__name__+'.MagData')
+ _log = logging.getLogger(__name__+'.MagData')
@property
def a(self):
@@ -85,68 +88,68 @@ class MagData(object):
self.magnitude = mag_vec.reshape((3,)+self.dim)
def __init__(self, a, magnitude):
- self.LOG.debug('Calling __init__')
+ self._log.debug('Calling __init__')
self.a = a
self.magnitude = magnitude
- self.LOG.debug('Created '+str(self))
+ self._log.debug('Created '+str(self))
def __repr__(self):
- self.LOG.debug('Calling __repr__')
+ self._log.debug('Calling __repr__')
return '%s(a=%r, magnitude=%r)' % (self.__class__, self.a, self.magnitude)
def __str__(self):
- self.LOG.debug('Calling __str__')
+ self._log.debug('Calling __str__')
return 'MagData(a=%s, dim=%s)' % (self.a, self.dim)
def __neg__(self): # -self
- self.LOG.debug('Calling __neg__')
+ self._log.debug('Calling __neg__')
return MagData(self.a, -self.magnitude)
def __add__(self, other): # self + other
- self.LOG.debug('Calling __add__')
+ self._log.debug('Calling __add__')
assert isinstance(other, (MagData, Number)), \
'Only MagData objects and scalar numbers (as offsets) can be added/subtracted!'
if isinstance(other, MagData):
- self.LOG.debug('Adding two MagData objects')
+ self._log.debug('Adding two MagData objects')
assert other.a == self.a, 'Added phase has to have the same grid spacing!'
assert other.magnitude.shape == (3,)+self.dim, \
'Added magnitude has to have the same dimensions!'
return MagData(self.a, self.magnitude+other.magnitude)
else: # other is a Number
- self.LOG.debug('Adding an offset')
+ self._log.debug('Adding an offset')
return MagData(self.a, self.magnitude+other)
def __sub__(self, other): # self - other
- self.LOG.debug('Calling __sub__')
+ self._log.debug('Calling __sub__')
return self.__add__(-other)
def __mul__(self, other): # self * other
- self.LOG.debug('Calling __mul__')
+ self._log.debug('Calling __mul__')
assert isinstance(other, Number), 'MagData objects can only be multiplied by numbers!'
return MagData(self.a, other*self.magnitude)
def __radd__(self, other): # other + self
- self.LOG.debug('Calling __radd__')
+ self._log.debug('Calling __radd__')
return self.__add__(other)
def __rsub__(self, other): # other - self
- self.LOG.debug('Calling __rsub__')
+ self._log.debug('Calling __rsub__')
return -self.__sub__(other)
def __rmul__(self, other): # other * self
- self.LOG.debug('Calling __rmul__')
+ self._log.debug('Calling __rmul__')
return self.__mul__(other)
def __iadd__(self, other): # self += other
- self.LOG.debug('Calling __iadd__')
+ self._log.debug('Calling __iadd__')
return self.__add__(other)
def __isub__(self, other): # self -= other
- self.LOG.debug('Calling __isub__')
+ self._log.debug('Calling __isub__')
return self.__sub__(other)
def __imul__(self, other): # self *= other
- self.LOG.debug('Calling __imul__')
+ self._log.debug('Calling __imul__')
return self.__mul__(other)
def copy(self):
@@ -162,7 +165,7 @@ class MagData(object):
A copy of the :class:`~.MagData`.
'''
- self.LOG.debug('Calling copy')
+ self._log.debug('Calling copy')
return MagData(self.a, self.magnitude.copy())
def scale_down(self, n=1):
@@ -183,11 +186,15 @@ class MagData(object):
Only possible, if each axis length is a power of 2!
'''
- self.LOG.debug('Calling scale_down')
+ self._log.debug('Calling scale_down')
assert n > 0 and isinstance(n, (int, long)), 'n must be a positive integer!'
- assert np.all([d % 2**n == 0 for d in self.dim]), 'Dimensions must a multiples of 2!'
self.a = self.a * 2**n
for t in range(n):
+ # Pad if necessary:
+ pz, py, px = self.dim[0] % 2, self.dim[1] % 2, self.dim[2] % 2
+ if pz != 0 or py != 0 or px != 0:
+ self.magnitude = np.pad(self.magnitude, ((0, 0), (0, pz), (0, py), (0, px)),
+ mode='constant')
# Create coarser grid for the magnetization:
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))
@@ -212,7 +219,7 @@ class MagData(object):
-----
Acts in place and changes dimensions and grid spacing accordingly.
'''
- self.LOG.debug('Calling scale_up')
+ self._log.debug('Calling scale_up')
assert n > 0 and isinstance(n, (int, long)), 'n must be a positive integer!'
assert 5 > order >= 0 and isinstance(order, (int, long)), \
'order must be a positive integer between 0 and 5!'
@@ -246,7 +253,7 @@ class MagData(object):
assert z_pad >= 0 and isinstance(z_pad, (int, long)), 'z_pad must be a positive integer!'
self.magnitude = np.pad(self.magnitude,
((0, 0), (z_pad, z_pad), (y_pad, y_pad), (x_pad, x_pad)),
- mode='constant', constant_values=0)
+ mode='constant')
def get_mask(self, threshold=0):
'''Mask all pixels where the amplitude of the magnetization lies above `threshold`.
@@ -326,7 +333,7 @@ class MagData(object):
None
'''
- self.LOG.debug('Calling save_to_llg')
+ self._log.debug('Calling save_to_llg')
a = self.a * 1.0E-9 / 1.0E-2 # from nm to cm
# Create 3D meshgrid and reshape it and the magnetization into a list where x varies first:
zz, yy, xx = np.mgrid[a/2:(self.dim[0]*a-a/2):self.dim[0]*1j,
@@ -356,7 +363,7 @@ class MagData(object):
A :class:`~.MagData` object containing the loaded data.
'''
- cls.LOG.debug('Calling load_from_llg')
+ cls._log.debug('Calling load_from_llg')
SCALE = 1.0E-9 / 1.0E-2 # From cm to nm
data = np.genfromtxt(filename, skip_header=2)
dim = tuple(np.genfromtxt(filename, dtype=int, skip_header=1, skip_footer=len(data[:, 0])))
@@ -378,7 +385,7 @@ class MagData(object):
None
'''
- self.LOG.debug('Calling save_to_netcdf4')
+ self._log.debug('Calling save_to_netcdf4')
mag_file = netCDF4.Dataset(filename, 'w', format='NETCDF4')
mag_file.a = self.a
mag_file.createDimension('comp', 3) # Number of components
@@ -404,7 +411,7 @@ class MagData(object):
A :class:`~.MagData` object containing the loaded data.
'''
- cls.LOG.debug('Calling copy')
+ cls._log.debug('Calling copy')
mag_file = netCDF4.Dataset(filename, 'r', format='NETCDF4')
a = mag_file.a
magnitude = mag_file.variables['magnitude'][...]
@@ -412,7 +419,7 @@ class MagData(object):
return MagData(a, magnitude)
def quiver_plot(self, title='Magnetization Distribution', axis=None, proj_axis='z',
- ax_slice=None, log=False, scaled=True, show=True):
+ ax_slice=None, log=False, scaled=True):
'''Plot a slice of the magnetization as a quiver plot.
Parameters
@@ -430,39 +437,38 @@ class MagData(object):
Takes the Default is False.
scaled : boolean, optional
Normalizes the plotted arrows in respect to the highest one. Default is True.
- show: bool, optional
- A switch which determines if the plot is shown at the end of plotting. Default is True.
+
Returns
-------
axis: :class:`~matplotlib.axes.AxesSubplot`
The axis on which the graph is plotted.
'''
- self.LOG.debug('Calling quiver_plot')
+ self._log.debug('Calling quiver_plot')
assert proj_axis == 'z' or proj_axis == 'y' or proj_axis == 'x', \
'Axis has to be x, y or z (as string).'
if proj_axis == 'z': # Slice of the xy-plane with z = ax_slice
- self.LOG.debug('proj_axis == z')
+ self._log.debug('proj_axis == z')
if ax_slice is None:
- self.LOG.debug('ax_slice is None')
+ self._log.debug('ax_slice is None')
ax_slice = int(self.dim[0]/2.)
plot_u = np.copy(self.magnitude[0][ax_slice, ...]) # x-component
plot_v = np.copy(self.magnitude[1][ax_slice, ...]) # y-component
u_label = 'x [px]'
v_label = 'y [px]'
elif proj_axis == 'y': # Slice of the xz-plane with y = ax_slice
- self.LOG.debug('proj_axis == y')
+ self._log.debug('proj_axis == y')
if ax_slice is None:
- self.LOG.debug('ax_slice is None')
+ self._log.debug('ax_slice is None')
ax_slice = int(self.dim[1]/2.)
plot_u = np.copy(self.magnitude[0][:, ax_slice, :]) # x-component
plot_v = np.copy(self.magnitude[2][:, ax_slice, :]) # z-component
u_label = 'x [px]'
v_label = 'z [px]'
elif proj_axis == 'x': # Slice of the yz-plane with x = ax_slice
- self.LOG.debug('proj_axis == x')
+ self._log.debug('proj_axis == x')
if ax_slice is None:
- self.LOG.debug('ax_slice is None')
+ self._log.debug('ax_slice is None')
ax_slice = int(self.dim[2]/2.)
plot_u = np.swapaxes(np.copy(self.magnitude[2][..., ax_slice]), 0, 1) # z-component
plot_v = np.swapaxes(np.copy(self.magnitude[1][..., ax_slice]), 0, 1) # y-component
@@ -470,7 +476,7 @@ class MagData(object):
v_label = 'y [px]'
# If no axis is specified, a new figure is created:
if axis is None:
- self.LOG.debug('axis is None')
+ self._log.debug('axis is None')
fig = plt.figure(figsize=(8.5, 7))
axis = fig.add_subplot(1, 1, 1)
axis.set_aspect('equal')
@@ -498,9 +504,6 @@ class MagData(object):
axis.tick_params(axis='both', which='major', labelsize=14)
axis.xaxis.set_major_locator(MaxNLocator(nbins=9, integer=True))
axis.yaxis.set_major_locator(MaxNLocator(nbins=9, integer=True))
- # Show Plot:
- if show:
- plt.show()
return axis
def quiver_plot3d(self, title='Magnetization Distribution'):
@@ -515,7 +518,7 @@ class MagData(object):
None
'''
- self.LOG.debug('Calling quiver_plot3D')
+ self._log.debug('Calling quiver_plot3D')
from mayavi import mlab
a = self.a
@@ -551,7 +554,7 @@ class MagData(object):
None
'''
- self.LOG.debug('Calling save_to_x3d')
+ self._log.debug('Calling save_to_x3d')
from lxml import etree
dim = self.dim
diff --git a/pyramid/phasemap.py b/pyramid/phasemap.py
index f5a3a4a..b33ab8c 100644
--- a/pyramid/phasemap.py
+++ b/pyramid/phasemap.py
@@ -18,6 +18,9 @@ import netCDF4
import logging
+__all__ = ['PhaseMap']
+
+
class PhaseMap(object):
'''Class for storing phase map data.
@@ -51,7 +54,7 @@ class PhaseMap(object):
'''
- LOG = logging.getLogger(__name__)
+ _log = logging.getLogger(__name__)
UNITDICT = {u'rad': 1E0,
u'mrad': 1E3,
@@ -141,72 +144,72 @@ class PhaseMap(object):
self._unit = unit
def __init__(self, a, phase, unit='rad'):
- self.LOG.debug('Calling __init__')
+ self._log.debug('Calling __init__')
self.a = a
self.phase = phase
self.unit = unit
- self.LOG.debug('Created '+str(self))
+ self._log.debug('Created '+str(self))
def __repr__(self):
- self.LOG.debug('Calling __repr__')
+ self._log.debug('Calling __repr__')
return '%s(a=%r, phase=%r, unit=%r)' % \
(self.__class__, self.a, self.phase, self.unit)
def __str__(self):
- self.LOG.debug('Calling __str__')
+ self._log.debug('Calling __str__')
return 'PhaseMap(a=%s, dim_uv=%s)' % (self.a, self.dim_uv)
def __neg__(self): # -self
- self.LOG.debug('Calling __neg__')
+ self._log.debug('Calling __neg__')
return PhaseMap(self.a, -self.phase, self.unit)
def __add__(self, other): # self + other
- self.LOG.debug('Calling __add__')
+ self._log.debug('Calling __add__')
assert isinstance(other, (PhaseMap, Number)), \
'Only PhaseMap objects and scalar numbers (as offsets) can be added/subtracted!'
if isinstance(other, PhaseMap):
- self.LOG.debug('Adding two PhaseMap objects')
+ self._log.debug('Adding two PhaseMap objects')
assert other.a == self.a, 'Added phase has to have the same grid spacing!'
assert other.phase.shape == self.dim_uv, \
'Added magnitude has to have the same dimensions!'
return PhaseMap(self.a, self.phase+other.phase, self.unit)
else: # other is a Number
- self.LOG.debug('Adding an offset')
+ self._log.debug('Adding an offset')
return PhaseMap(self.a, self.phase+other, self.unit)
def __sub__(self, other): # self - other
- self.LOG.debug('Calling __sub__')
+ self._log.debug('Calling __sub__')
return self.__add__(-other)
def __mul__(self, other): # self * other
- self.LOG.debug('Calling __mul__')
+ self._log.debug('Calling __mul__')
assert (isinstance(other, Number)
or (isinstance(other, np.ndarray) and other.shape == self.dim_uv)), \
'PhaseMap objects can only be multiplied by scalar numbers or fitting arrays!'
return PhaseMap(self.a, other*self.phase, self.unit)
def __radd__(self, other): # other + self
- self.LOG.debug('Calling __radd__')
+ self._log.debug('Calling __radd__')
return self.__add__(other)
def __rsub__(self, other): # other - self
- self.LOG.debug('Calling __rsub__')
+ self._log.debug('Calling __rsub__')
return -self.__sub__(other)
def __rmul__(self, other): # other * self
- self.LOG.debug('Calling __rmul__')
+ self._log.debug('Calling __rmul__')
return self.__mul__(other)
def __iadd__(self, other): # self += other
- self.LOG.debug('Calling __iadd__')
+ self._log.debug('Calling __iadd__')
return self.__add__(other)
def __isub__(self, other): # self -= other
- self.LOG.debug('Calling __isub__')
+ self._log.debug('Calling __isub__')
return self.__sub__(other)
def __imul__(self, other): # self *= other
- self.LOG.debug('Calling __imul__')
+ self._log.debug('Calling __imul__')
return self.__mul__(other)
def save_to_txt(self, filename='..\output\phasemap_output.txt'):
@@ -223,7 +226,7 @@ class PhaseMap(object):
None
'''
- self.LOG.debug('Calling save_to_txt')
+ self._log.debug('Calling save_to_txt')
with open(filename, 'w') as phase_file:
phase_file.write('{}\n'.format(filename.replace('.txt', '')))
phase_file.write('grid spacing = {} nm\n'.format(self.a))
@@ -244,7 +247,7 @@ class PhaseMap(object):
A :class:`~.PhaseMap` object containing the loaded data.
'''
- cls.LOG.debug('Calling load_from_txt')
+ cls._log.debug('Calling load_from_txt')
with open(filename, 'r') as phase_file:
phase_file.readline() # Headerline is not used
a = float(phase_file.readline()[15:-4])
@@ -265,7 +268,7 @@ class PhaseMap(object):
None
'''
- self.LOG.debug('Calling save_to_netcdf4')
+ self._log.debug('Calling save_to_netcdf4')
phase_file = netCDF4.Dataset(filename, 'w', format='NETCDF4')
phase_file.a = self.a
phase_file.createDimension('v_dim', self.dim_uv[0])
@@ -289,16 +292,15 @@ class PhaseMap(object):
A :class:`~.PhaseMap` object containing the loaded data.
'''
- cls.LOG.debug('Calling load_from_netcdf4')
+ cls._log.debug('Calling load_from_netcdf4')
phase_file = netCDF4.Dataset(filename, 'r', format='NETCDF4')
a = phase_file.a
phase = phase_file.variables['phase'][:]
- #phase = phase[28:36, 28:36]
phase_file.close()
return PhaseMap(a, phase)
def display_phase(self, title='Phase Map', cmap='RdBu', limit=None,
- norm=None, axis=None, cbar=True, show=True):
+ norm=None, axis=None, cbar=True):
'''Display the phasemap as a colormesh.
Parameters
@@ -318,8 +320,6 @@ class PhaseMap(object):
Axis on which the graph is plotted. Creates a new figure if none is specified.
cbar : bool, optional
A switch determining if the colorbar should be plotted or not. Default is True.
- show : bool, optional
- A switch which determines if the plot is shown at the end of plotting.
Returns
-------
@@ -327,7 +327,7 @@ class PhaseMap(object):
The axis on which the graph is plotted.
'''
- self.LOG.debug('Calling display_phase')
+ self._log.debug('Calling display_phase')
# Take units into consideration:
phase = self.phase * self.UNITDICT[self.unit]
if limit is None:
@@ -358,13 +358,10 @@ class PhaseMap(object):
cbar = fig.colorbar(im, cax=cbar_ax)
cbar.ax.tick_params(labelsize=14)
cbar.set_label(u'phase shift [{}]'.format(self.unit), fontsize=15)
- # Show plot:
- if show:
- plt.show()
# Return plotting axis:
return axis
- def display_phase3d(self, title='Phase Map', cmap='RdBu', show=True):
+ def display_phase3d(self, title='Phase Map', cmap='RdBu'):
'''Display the phasemap as a 3-D surface with contourplots.
Parameters
@@ -374,8 +371,6 @@ class PhaseMap(object):
cmap : string, optional
The :class:`~matplotlib.colors.Colormap` which is used for the plot as a string.
The default is 'RdBu'.
- show: bool, optional
- A switch which determines if the plot is shown at the end of plotting. Default is True.
Returns
-------
@@ -383,7 +378,7 @@ class PhaseMap(object):
The axis on which the graph is plotted.
'''
- self.LOG.debug('Calling display_phase3d')
+ self._log.debug('Calling display_phase3d')
# Take units into consideration:
phase = self.phase * self.UNITDICT[self.unit]
# Create figure and axis:
@@ -400,14 +395,11 @@ class PhaseMap(object):
axis.set_xlabel('u-axis [px]')
axis.set_ylabel('v-axis [px]')
axis.set_zlabel('phase shift [{}]'.format(self.unit))
- # Show Plot:
- if show:
- plt.show()
# Return plotting axis:
return axis
def display_holo(self, title=None, gain='auto', axis=None, grad_encode='bright',
- interpolation='none', show=True):
+ interpolation='none'):
'''Display the color coded holography image.
Parameters
@@ -435,7 +427,7 @@ class PhaseMap(object):
The axis on which the graph is plotted.
'''
- self.LOG.debug('Calling display_holo')
+ self._log.debug('Calling display_holo')
# Calculate gain if 'auto' is selected:
if gain == 'auto':
gain = 4 * 2*pi/self.phase.max()
@@ -453,19 +445,19 @@ class PhaseMap(object):
phase_saturation = np.dstack((saturation,)*4)
# Color code the angle and create the holography image:
if grad_encode == 'dark':
- self.LOG.debug('gradient encoding: dark')
+ self._log.debug('gradient encoding: dark')
rgba = self.HOLO_CMAP(phase_angle)
rgb = (255.999 * img_holo.T * saturation.T * rgba[:, :, :3].T).T.astype(np.uint8)
elif grad_encode == 'bright':
- self.LOG.debug('gradient encoding: bright')
+ self._log.debug('gradient encoding: bright')
rgba = self.HOLO_CMAP(phase_angle)+(1-phase_saturation)*self.HOLO_CMAP_INV(phase_angle)
rgb = (255.999 * img_holo.T * rgba[:, :, :3].T).T.astype(np.uint8)
elif grad_encode == 'color':
- self.LOG.debug('gradient encoding: color')
+ self._log.debug('gradient encoding: color')
rgba = self.HOLO_CMAP(phase_angle)
rgb = (255.999 * img_holo.T * rgba[:, :, :3].T).T.astype(np.uint8)
elif grad_encode == 'none':
- self.LOG.debug('gradient encoding: none')
+ self._log.debug('gradient encoding: none')
rgba = self.HOLO_CMAP(phase_angle)+self.HOLO_CMAP_INV(phase_angle)
rgb = (255.999 * img_holo.T * rgba[:, :, :3].T).T.astype(np.uint8)
else:
@@ -488,14 +480,11 @@ class PhaseMap(object):
axis.set_ylim(0, self.dim_uv[0])
axis.xaxis.set_major_locator(MaxNLocator(nbins=9, integer=True))
axis.yaxis.set_major_locator(MaxNLocator(nbins=9, integer=True))
- # Show Plot:
- if show:
- plt.show()
# Return plotting axis:
return axis
def display_combined(self, title='Combined Plot', cmap='RdBu', limit=None, norm=None,
- gain='auto', interpolation='none', grad_encode='bright', show=True):
+ gain='auto', interpolation='none', grad_encode='bright'):
'''Display the phase map and the resulting color coded holography image in one plot.
Parameters
@@ -522,8 +511,6 @@ class PhaseMap(object):
encodes the direction (without gradient strength), 'dark' modulates the gradient
strength with a factor between 0 and 1 and 'bright' (which is the default) encodes
the gradient strength with color saturation.
- show: bool, optional
- A switch which determines if the plot is shown at the end of plotting. Default is True.
Returns
-------
@@ -531,39 +518,35 @@ class PhaseMap(object):
The axes on which the graphs are plotted.
'''
- self.LOG.debug('Calling display_combined')
+ self._log.debug('Calling display_combined')
# Create combined plot and set title:
fig = plt.figure(figsize=(16, 7))
fig.suptitle(title, fontsize=20)
# Plot holography image:
holo_axis = fig.add_subplot(1, 2, 1, aspect='equal')
self.display_holo(gain=gain, axis=holo_axis, interpolation=interpolation,
- show=False, grad_encode=grad_encode)
+ grad_encode=grad_encode)
# Plot phase map:
phase_axis = fig.add_subplot(1, 2, 2, aspect='equal')
fig.subplots_adjust(right=0.85)
- self.display_phase(cmap='RdBu', limit=limit, norm=norm, axis=phase_axis, show=False)
- # Show Plot:
- if show:
- plt.show()
+ self.display_phase(cmap='RdBu', limit=limit, norm=norm, axis=phase_axis)
# Return the plotting axes:
return phase_axis, holo_axis
@classmethod
- def make_color_wheel(cls, show=True):
+ def make_color_wheel(cls):
'''Display a color wheel to illustrate the color coding of the gradient direction.
Parameters
----------
- show: bool, optional
- A switch which determines if the plot is shown at the end of plotting. Default is True.
+ None
Returns
-------
None
'''
- cls.LOG.debug('Calling make_color_wheel')
+ cls._log.debug('Calling make_color_wheel')
x = np.linspace(-256, 256, num=512)
y = np.linspace(-256, 256, num=512)
xx, yy = np.meshgrid(x, y)
@@ -582,6 +565,3 @@ class PhaseMap(object):
axis.imshow(color_wheel, origin='lower')
axis.xaxis.set_major_locator(NullLocator())
axis.yaxis.set_major_locator(NullLocator())
- # Show Plot:
- if show:
- plt.show()
diff --git a/pyramid/phasemapper.py b/pyramid/phasemapper.py
index 34d5b28..2c2c0e9 100644
--- a/pyramid/phasemapper.py
+++ b/pyramid/phasemapper.py
@@ -18,9 +18,20 @@ from pyramid.phasemap import PhaseMap
from pyramid.projector import SimpleProjector
from pyramid.kernel import Kernel
+import jutil.fft as jfft
+
import logging
+__all__ = ['PhaseMapperRDFC', 'PhaseMapperRDRC', 'PhaseMapperFDFC', 'pm']
+
+PHI_0 = 2067.83 # magnetic flux in T*nm²
+H_BAR = 6.626E-34 # Planck constant in J*s
+M_E = 9.109E-31 # electron mass in kg
+Q_E = 1.602E-19 # electron charge in C
+C = 2.998E8 # speed of light in m/s
+
+
class PhaseMapper(object):
'''Abstract base class for the phase calculation from a 2-dimensional distribution.
@@ -35,21 +46,21 @@ class PhaseMapper(object):
'''
__metaclass__ = abc.ABCMeta
- LOG = logging.getLogger(__name__+'.PhaseMapper')
+ _log = logging.getLogger(__name__+'.PhaseMapper')
@abc.abstractmethod
def __call__(self, mag_data):
- self.LOG.debug('Calling __call__')
+ self._log.debug('Calling __call__')
raise NotImplementedError()
@abc.abstractmethod
def jac_dot(self, vector):
- self.LOG.debug('Calling jac_dot')
+ self._log.debug('Calling jac_dot')
raise NotImplementedError()
@abc.abstractmethod
def jac_T_dot(self, vector):
- self.LOG.debug('Calling jac_T_dot')
+ self._log.debug('Calling jac_T_dot')
raise NotImplementedError()
@@ -82,41 +93,67 @@ class PhaseMapperRDFC(PhaseMapper):
'''
- LOG = logging.getLogger(__name__+'.PhaseMapperRDFC')
+ _log = logging.getLogger(__name__+'.PhaseMapperRDFC')
def __init__(self, kernel):
- self.LOG.debug('Calling __init__')
+ self._log.debug('Calling __init__')
self.kernel = kernel
self.m = np.prod(kernel.dim_uv)
self.n = 2 * self.m
- self.LOG.debug('Created '+str(self))
+ if self.kernel.use_fftw: # if possible use FFTW
+ import pyfftw
+ n = pyfftw.simd_alignment
+ self.u_mag = pyfftw.n_byte_align_empty(self.kernel.dim_uv, n, dtype='float32')
+ self.v_mag = pyfftw.n_byte_align_empty(self.kernel.dim_uv, n, dtype='float32')
+# self.phase = pyfftw.n_byte_align_empty(self.kernel.dim_uv, n, dtype='float32')
+ self.u_mag_fft = pyfftw.n_byte_align_empty(self.kernel.dim_fft, n, dtype='complex64')
+ self.v_mag_fft = pyfftw.n_byte_align_empty(self.kernel.dim_fft, n, dtype='complex64')
+ self.phase_fft = pyfftw.n_byte_align_empty(self.kernel.dim_fft, n, dtype='complex64')
+ self._rfft2 = pyfftw.builders.rfftn(self.u_mag, self.kernel.dim_pad,
+ threads=self.kernel.threads)
+ self._irfft2 = pyfftw.builders.irfftn(self.phase_fft, self.kernel.dim_pad,
+ threads=self.kernel.threads)
+ self._rfft2_adj = self._rfft2_adj_fftw
+ self._irfft2_adj = self._irfft2_adj_fftw
+ else: # otherwise use numpy
+ self.u_mag = np.empty(self.kernel.dim_uv, dtype=float)
+ self.v_mag = np.empty(self.kernel.dim_uv, dtype=float)
+# self.phase = np.empty(self.kernel.dim_uv, dtype=float)
+ self.u_mag_fft = np.empty(self.kernel.dim_fft, dtype=complex)
+ self.v_mag_fft = np.empty(self.kernel.dim_fft, dtype=complex)
+ self.phase_fft = np.empty(self.kernel.dim_fft, dtype=complex)
+ self._rfft2 = lambda x: np.fft.rfftn(x, self.kernel.dim_pad)
+ self._irfft2 = lambda x: np.fft.irfftn(x, self.kernel.dim_pad)
+ self._rfft2_adj = jfft.irfft2_adj
+ self._rfft2_adj = jfft.irfft2_adj
+ self._log.debug('Created '+str(self))
def __repr__(self):
- self.LOG.debug('Calling __repr__')
+ self._log.debug('Calling __repr__')
return '%s(kernel=%r)' % (self.__class__, self.kernel)
def __str__(self):
- self.LOG.debug('Calling __str__')
+ self._log.debug('Calling __str__')
return 'PhaseMapperRDFC(kernel=%s)' % (self.kernel)
def __call__(self, mag_data):
- self.LOG.debug('Calling __call__')
+ self._log.debug('Calling __call__')
assert isinstance(mag_data, MagData), 'Only MagData objects can be mapped!'
assert mag_data.a == self.kernel.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.kernel.dim_uv, 'Dimensions do not match!'
# Process input parameters:
- u_mag, v_mag = mag_data.magnitude[0:2, 0, ...]
- return PhaseMap(mag_data.a, self._convolve(u_mag, v_mag))
+ self.u_mag[...], self.v_mag[...] = mag_data.magnitude[0:2, 0, ...]
+ return PhaseMap(mag_data.a, self._convolve())
- def _convolve(self, u_mag, v_mag):
+ def _convolve(self):
# Fourier transform the projected magnetisation:
- u_mag_fft = np.fft.rfftn(u_mag, self.kernel.dim_fft)
- v_mag_fft = np.fft.rfftn(v_mag, self.kernel.dim_fft)
+ self.u_mag_fft[...] = self._rfft2(self.u_mag)
+ self.v_mag_fft[...] = self._rfft2(self.v_mag)
# Convolve the magnetization with the kernel in Fourier space:
- phase_fft = u_mag_fft*self.kernel.u_fft - v_mag_fft*self.kernel.v_fft
+ self.phase_fft[...] = self.u_mag_fft*self.kernel.u_fft - self.v_mag_fft*self.kernel.v_fft
# Return the result:
- return np.fft.irfftn(phase_fft, self.kernel.dim_fft)[self.kernel.slice_fft]
+ return self._irfft2(self.phase_fft)[self.kernel.slice_fft]
def jac_dot(self, vector):
'''Calculate the product of the Jacobi matrix with a given `vector`.
@@ -134,13 +171,12 @@ class PhaseMapperRDFC(PhaseMapper):
Product of the Jacobi matrix (which is not explicitely calculated) with the vector.
'''
- self.LOG.debug('Calling jac_dot')
+# self._log.debug('Calling jac_dot') # TODO: Profiler says this was slow...
assert len(vector) == self.n, \
'vector size not compatible! vector: {}, size: {}'.format(len(vector), self.n)
- u_mag, v_mag = np.reshape(vector, (2,)+self.kernel.dim_uv)
- result = self._convolve(u_mag, v_mag).reshape(-1)
+ self.u_mag[...], self.v_mag[...] = np.reshape(vector, (2,)+self.kernel.dim_uv)
+ result = self._convolve().reshape(-1)
return result
- #return self.kernel._multiply_jacobi(vector)
def jac_T_dot(self, vector):
'''Calculate the product of the transposed Jacobi matrix with a given `vector`.
@@ -158,61 +194,84 @@ class PhaseMapperRDFC(PhaseMapper):
the vector, which has ``2*N**2`` entries like a 2D magnetic projection.
'''
- self.LOG.debug('Calling jac_T_dot')
+# self._log.debug('Calling jac_T_dot') # TODO: Profiler says this was slow...
assert len(vector) == self.m, \
'vector size not compatible! vector: {}, size: {}'.format(len(vector), self.m)
- result = np.zeros(self.n)
+ phase = vector.reshape(self.kernel.dim_uv)
+ p0 = self.kernel.dim_pad[0]-self.kernel.dim_uv[0]
+ p1 = self.kernel.dim_pad[1]-self.kernel.dim_uv[1]
+ phase = np.pad(phase, ((0, p0), (0, p1)), 'constant')
+ phase_fft = jfft.irfft2_adj(phase)
+ u_mag_fft = phase_fft * self.kernel.u_fft
+ v_mag_fft = phase_fft * -self.kernel.v_fft
+ u_mag = jfft.rfft2_adj(u_mag_fft, self.kernel.dim_pad[1])[self.kernel.slice_fft]
+ v_mag = jfft.rfft2_adj(v_mag_fft, self.kernel.dim_pad[1])[self.kernel.slice_fft]
+ return -np.concatenate((u_mag.flatten(), v_mag.flatten())) # TODO: why the minus?
+
+# TODO: save arrays in PADDED form (kernel and the self.u/v/... stuff) to save padding time
+# TODO: also gets rid of (..., self.kernel.dim_pad) argument
+# if a.shape[axis] != n:
+# s = list(a.shape)
+# if s[axis] > n:
+# index = [slice(None)]*len(s)
+# index[axis] = slice(0, n)
+# a = a[index]
+# else:
+# index = [slice(None)]*len(s)
+# index[axis] = slice(0, s[axis])
+# s[axis] = n
+# z = zeros(s, a.dtype.char)
+# z[index] = a
+# a = z
+
+ def _rfft2_adj_fftw(x, n=None):
+ """
+ Adjoint of the 2-D Real Fast Fourier Transform, that is the product of the adjoint
+ Jacobian of the numpy rfftn with vector x.
- import jutil.fft as jfft
- from jutil.taketime import TakeTime
-
- with TakeTime('fft_adjoint'):
-
- phase = vector.reshape(self.kernel.dim_uv)
- p0 = self.kernel.dim_fft[0]-self.kernel.dim_uv[0]
- p1 = self.kernel.dim_fft[1]-self.kernel.dim_uv[1]
- phase = np.pad(phase, ((0, p0), (0, p1)), 'constant')
-
- phase_fft = jfft.irfft2_adj(phase)
- u_mag_fft = phase_fft * self.kernel.u_fft
- v_mag_fft = phase_fft * -self.kernel.v_fft
- u_mag = jfft.rfft2_adj(u_mag_fft, self.kernel.dim_fft[1])[self.kernel.slice_fft]
- v_mag = jfft.rfft2_adj(v_mag_fft, self.kernel.dim_fft[1])[self.kernel.slice_fft]
- result = -np.concatenate((u_mag.flatten(), v_mag.flatten()))
-
-# # TODO: directly? ask Jörn again!
-# phase = vector.reshape(self.kernel.dim_uv)
-# phase_fft = np.fft.fft(phase, self.kernel.dim_fft)
-# u_mag_fft = phase_fft * self.kernel.u_fft_compl
-# v_mag_fft = phase_fft * -self.kernel.v_fft_compl
-# u_mag = np.fft.ifft(u_mag_fft, self.kernel.dim_fft)[self.kernel.slice_fft]
-# v_mag = np.fft.ifft(u_mag_fft, self.kernel.dim_fft)[self.kernel.slice_fft]
-# result = np.concatenate((u_mag.flatten(), v_mag.flatten()))
-# return result
-
- #TODO: np.split()
- #TODO: TakeTime()
- #TODO: 'constant' in np.pad()
-
-# result[s] = np.sum(vector*self.u[v_min:v_max, u_min:u_max].reshape(-1))
-# result[s+self.m] = np.sum(vector*-self.v[v_min:v_max, u_min:u_max].reshape(-1))
-#
-# # Fourier transform the projected magnetisation:
-# u_mag_fft = np.fft.rfftn(u_mag, self.kernel.dim_fft)
-# v_mag_fft = np.fft.rfftn(v_mag, self.kernel.dim_fft)
-# # Convolve the magnetization with the kernel in Fourier space:
-# phase_fft = u_mag_fft*self.kernel.u_fft - v_mag_fft*self.kernel.v_fft
-# # Return the result:
-# return np.fft.irfftn(phase_fft, self.kernel.dim_fft)[self.kernel.slice_fft]
-
-# with TakeTime('oldschool'):
-# compare = np.zeros(self.n)
-# nc.jac_T_dot_real_convolve(self.kernel.dim_uv[0], self.kernel.dim_uv[1],
-# self.kernel.u, self.kernel.v, vector, compare)
-
-# import pdb; pdb.set_trace()
-# assert np.testing.assert_almost_equal(result, compare)
- return result
+ Parameters
+ ----------
+ x : array_like
+ n : int, optional
+ Length of second dimension of original array (the one, the size of which is
+ changed by employing the rfft2)
+
+ Returns
+ -------
+ array_like
+ """
+ print 'rfftn2_adj input:', x.shape
+ if n is None:
+ n = 2 * (x.shape[1] - 1)
+ xp = np.zeros((x.shape[0], n), dtype=x.dtype)
+ xp[:, :x.shape[1]] = x
+ print 'rfftn2_adj output:', np.fft.ifftn(xp).real * x.shape[0] * n
+ return np.fft.ifftn(xp).real * x.shape[0] * n
+
+ def _irfft2_adj_fftw(x):
+ """
+ Adjoint of the 2-D Inverse Real Fast Fourier Transform, that is the product of the
+ adjoint Jacobian of the numpy irfftn with vector x.
+
+ Parameters
+ ----------
+ x : array_like
+
+ Returns
+ -------
+ array_like
+ """
+ print 'irfftn2_adj input:', x.shape
+ n_out = x.shape[1] / 2 + 1
+ xp = np.fft.fft(x, axis=1) / x.shape[1]
+ if x.shape[1] % 2 == 0:
+ xp[:, 1:n_out - 1] += np.conj(xp[:, :n_out-1:-1])
+ # xp[:, n_out - 1] = xp[:, n_out - 1].real
+ else:
+ xp[:, 1:n_out] += np.conj(xp[:, :n_out-1:-1])
+ # xp[:, 0] = xp[:, 0].real
+ print 'rfftn2_adj output:', np.fft.ifftn(xp).real * x.shape[0] * n
+ return np.fft.fft(xp[:, :n_out], axis=0) / xp.shape[0]
class PhaseMapperRDRC(PhaseMapper):
@@ -246,29 +305,29 @@ class PhaseMapperRDRC(PhaseMapper):
'''
- LOG = logging.getLogger(__name__+'.PhaseMapperRDRC')
+ _log = logging.getLogger(__name__+'.PhaseMapperRDRC')
def __init__(self, kernel, threshold=0, numcore=True):
- self.LOG.debug('Calling __init__')
+ self._log.debug('Calling __init__')
self.kernel = kernel
self.threshold = threshold
self.numcore = numcore
self.m = np.prod(kernel.dim_uv)
self.n = 2 * self.m
- self.LOG.debug('Created '+str(self))
+ self._log.debug('Created '+str(self))
def __repr__(self):
- self.LOG.debug('Calling __repr__')
+ self._log.debug('Calling __repr__')
return '%s(kernel=%r, threshold=%r, numcore=%r)' % \
(self.__class__, self.kernel, self.threshold, self.numcore)
def __str__(self):
- self.LOG.debug('Calling __str__')
+ self._log.debug('Calling __str__')
return 'PhaseMapperRDRC(kernel=%s, threshold=%s, numcore=%s)' % \
(self.kernel, self.threshold, self.numcore)
def __call__(self, mag_data):
- self.LOG.debug('Calling __call__')
+ self._log.debug('Calling __call__')
dim_uv = self.kernel.dim_uv
assert isinstance(mag_data, MagData), 'Only MagData objects can be mapped!'
assert mag_data.a == self.kernel.a, 'Grid spacing has to match!'
@@ -314,7 +373,7 @@ class PhaseMapperRDRC(PhaseMapper):
Product of the Jacobi matrix (which is not explicitely calculated) with the vector.
'''
- self.LOG.debug('Calling jac_dot')
+ self._log.debug('Calling jac_dot')
assert len(vector) == self.n, \
'vector size not compatible! vector: {}, size: {}'.format(len(vector), self.n)
v_dim, u_dim = self.kernel.dim_uv
@@ -350,7 +409,7 @@ class PhaseMapperRDRC(PhaseMapper):
the vector, which has ``2*N**2`` entries like a 2D magnetic projection.
'''
- self.LOG.debug('Calling jac_T_dot')
+ self._log.debug('Calling jac_T_dot')
assert len(vector) == self.m, \
'vector size not compatible! vector: {}, size: {}'.format(len(vector), self.m)
v_dim, u_dim = self.dim_uv
@@ -396,31 +455,30 @@ class PhaseMapperFDFC(PhaseMapper):
'''
- LOG = logging.getLogger(__name__+'.PhaseMapperFDFC')
- PHI_0 = -2067.83 # magnetic flux in T*nm²
+ _log = logging.getLogger(__name__+'.PhaseMapperFDFC')
def __init__(self, a, dim_uv, b_0=1, padding=0):
- self.LOG.debug('Calling __init__')
+ self._log.debug('Calling __init__')
self.a = a
self.dim_uv = dim_uv
self.b_0 = b_0
self.padding = padding
self.m = np.prod(dim_uv)
self.n = 2 * self.m
- self.LOG.debug('Created '+str(self))
+ self._log.debug('Created '+str(self))
def __repr__(self):
- self.LOG.debug('Calling __repr__')
+ self._log.debug('Calling __repr__')
return '%s(a=%r, dim_uv=%r, b_0=%r, padding=%r)' % \
(self.__class__, self.a, self.dim_uv, self.b_0, self.padding)
def __str__(self):
- self.LOG.debug('Calling __str__')
+ self._log.debug('Calling __str__')
return 'PhaseMapperFDFC(a=%s, dim_uv=%s, b_0=%s, padding=%s)' % \
(self.a, self.dim_uv, self.b_0, self.padding)
def __call__(self, mag_data):
- self.LOG.debug('Calling __call__')
+ self._log.debug('Calling __call__')
assert isinstance(mag_data, MagData), 'Only MagData objects can be mapped!'
assert mag_data.a == self.a, 'Grid spacing has to match!'
assert mag_data.dim[0] == 1, 'Magnetic distribution must be 2-dimensional!'
@@ -430,8 +488,8 @@ class PhaseMapperFDFC(PhaseMapper):
# Create zero padded matrices:
u_pad = int(u_dim/2 * self.padding)
v_pad = int(v_dim/2 * self.padding)
- u_mag_pad = np.pad(u_mag, ((v_pad, v_pad), (u_pad, u_pad)), 'constant', constant_values=0)
- v_mag_pad = np.pad(v_mag, ((v_pad, v_pad), (u_pad, u_pad)), 'constant', constant_values=0)
+ u_mag_pad = np.pad(u_mag, ((v_pad, v_pad), (u_pad, u_pad)), 'constant')
+ v_mag_pad = np.pad(v_mag, ((v_pad, v_pad), (u_pad, u_pad)), 'constant')
# Fourier transform of the two components:
u_mag_fft = np.fft.fftshift(np.fft.rfft2(u_mag_pad), axes=0)
v_mag_fft = np.fft.fftshift(np.fft.rfft2(v_mag_pad), axes=0)
@@ -440,7 +498,7 @@ class PhaseMapperFDFC(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*self.a) / (2*self.PHI_0)
+ coeff = - (1j*self.b_0*self.a) / (2*PHI_0) # Minus because of negative z-direction
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))
@@ -448,7 +506,7 @@ class PhaseMapperFDFC(PhaseMapper):
return PhaseMap(mag_data.a, phase)
def jac_dot(self, vector):
- self.LOG.debug('Calling jac_dot')
+ self._log.debug('Calling jac_dot')
'''Calculate the product of the Jacobi matrix with a given `vector`.
Parameters
@@ -464,7 +522,7 @@ class PhaseMapperFDFC(PhaseMapper):
Product of the Jacobi matrix (which is not explicitely calculated) with the vector.
'''
- self.LOG.debug('Calling jac_dot')
+ self._log.debug('Calling jac_dot')
assert len(vector) == self.n, \
'vector size not compatible! vector: {}, size: {}'.format(len(vector), self.n)
mag_proj = MagData(self.a, np.zeros((3, 1)+self.dim_uv))
@@ -474,7 +532,7 @@ class PhaseMapperFDFC(PhaseMapper):
return self(mag_proj).phase_vec
def jac_T_dot(self, vector):
- self.LOG.debug('Calling jac_T_dot')
+ self._log.debug('Calling jac_T_dot')
raise NotImplementedError()
# TODO: Implement!
@@ -505,14 +563,10 @@ class PhaseMapperElectric(PhaseMapper):
'''
- LOG = logging.getLogger(__name__+'.PhaseMapperElectric')
- H_BAR = 6.626E-34 # Planck constant in J*s
- M_E = 9.109E-31 # electron mass in kg
- Q_E = 1.602E-19 # electron charge in C
- C = 2.998E8 # speed of light in m/s
+ _log = logging.getLogger(__name__+'.PhaseMapperElectric')
def __init__(self, a, dim_uv, v_0=1, v_acc=30000, threshold=0):
- self.LOG.debug('Calling __init__')
+ self._log.debug('Calling __init__')
self.a = a
self.dim_uv = dim_uv
self.v_0 = v_0
@@ -521,24 +575,23 @@ class PhaseMapperElectric(PhaseMapper):
self.m = np.prod(self.dim_uv)
self.n = np.prod(self.dim_uv)
# Coefficient calculation:
- H_BAR, M_E, Q_E, C = self.H_BAR, self.M_E, self.Q_E, self.C
lam = H_BAR / np.sqrt(2 * M_E * Q_E * v_acc * (1 + Q_E*v_acc / (2*M_E*C**2)))
C_e = 2*pi*Q_E/lam * (Q_E*v_acc + M_E*C**2) / (Q_E*v_acc * (Q_E*v_acc + 2*M_E*C**2))
self.coeff = v_0 * C_e * a * 1E-9
- self.LOG.debug('Created '+str(self))
+ self._log.debug('Created '+str(self))
def __repr__(self):
- self.LOG.debug('Calling __repr__')
+ self._log.debug('Calling __repr__')
return '%s(a=%r, dim_uv=%r, v_0=%r, v_acc=%r, threshold=%r)' % \
(self.__class__, self.a, self.dim_uv, self.v_0, self.v_acc, self.threshold)
def __str__(self):
- self.LOG.debug('Calling __str__')
+ self._log.debug('Calling __str__')
return 'PhaseMapperElectric(a=%s, dim_uv=%s, v_0=%s, v_acc=%s, threshold=%s)' % \
(self.a, self.dim_uv, self.v_0, self.v_acc, self.threshold)
def __call__(self, mag_data):
- self.LOG.debug('Calling __call__')
+ self._log.debug('Calling __call__')
assert isinstance(mag_data, MagData), 'Only MagData objects can be mapped!'
assert mag_data.a == self.a, 'Grid spacing has to match!'
assert mag_data.dim[0] == 1, 'Magnetic distribution must be 2-dimensional!'
@@ -553,12 +606,12 @@ class PhaseMapperElectric(PhaseMapper):
return PhaseMap(mag_data.a, phase)
def jac_dot(self, vector):
- self.LOG.debug('Calling jac_dot')
+ self._log.debug('Calling jac_dot')
raise NotImplementedError()
# TODO: Implement?
def jac_T_dot(self, vector):
- self.LOG.debug('Calling jac_T_dot')
+ self._log.debug('Calling jac_T_dot')
raise NotImplementedError()
# TODO: Implement?
@@ -582,7 +635,7 @@ def pm(mag_data, axis='z', dim_uv=None, b_0=1):
mag_data : :class:`~pyramid.magdata.MagData`
The reconstructed magnetic distribution as a :class:`~.MagData` object.
- '''
+ ''' # TODO: use_fftw false to notes
projector = SimpleProjector(mag_data.dim, axis=axis, dim_uv=dim_uv)
- phasemapper = PhaseMapperRDFC(Kernel(mag_data.a, projector.dim_uv))
+ phasemapper = PhaseMapperRDFC(Kernel(mag_data.a, projector.dim_uv, use_fftw=False))
return phasemapper(projector(mag_data))
diff --git a/pyramid/projector.py b/pyramid/projector.py
index f40e996..f28487d 100644
--- a/pyramid/projector.py
+++ b/pyramid/projector.py
@@ -17,6 +17,9 @@ from pyramid.magdata import MagData
import logging
+__all__ = ['XTiltProjector', 'YTiltProjector', 'SimpleProjector']
+
+
class Projector(object):
'''Abstract base class representing a projection function.
@@ -53,11 +56,11 @@ class Projector(object):
'''
__metaclass__ = abc.ABCMeta
- LOG = logging.getLogger(__name__+'.Projector')
+ _log = logging.getLogger(__name__+'.Projector')
@abc.abstractmethod
def __init__(self, dim, dim_uv, weight, coeff):
- self.LOG.debug('Calling __init__')
+ self._log.debug('Calling __init__')
self.dim = dim
self.dim_uv = dim_uv
self.weight = weight
@@ -65,26 +68,26 @@ class Projector(object):
self.size_2d, self.size_3d = weight.shape
self.n = 3 * np.prod(dim)
self.m = 2 * np.prod(dim_uv)
- self.LOG.debug('Created '+str(self))
+ self._log.debug('Created '+str(self))
def __repr__(self):
- self.LOG.debug('Calling __repr__')
+ self._log.debug('Calling __repr__')
return '%s(dim=%r, dim_uv=%r, weight=%r, coeff=%r)' % \
(self.__class__, self.dim, self.dim_uv, self.weight, self.coeff)
def __str__(self):
- self.LOG.debug('Calling __str__')
+ self._log.debug('Calling __str__')
return 'Projector(dim=%s, dim_uv=%s, coeff=%s)' % (self.dim, self.dim_uv, self.coeff)
def __call__(self, mag_data):
- self.LOG.debug('Calling as function')
+ self._log.debug('Calling as function')
mag_proj = MagData(mag_data.a, np.zeros((3, 1)+self.dim_uv))
magnitude_proj = self.jac_dot(mag_data.mag_vec).reshape((2, )+self.dim_uv)
mag_proj.magnitude[0:2, 0, ...] = magnitude_proj
return mag_proj
def _vector_field_projection(self, vector):
- self.LOG.debug('Calling _vector_field_projection')
+# self._log.debug('Calling _vector_field_projection') # TODO: Profiler says this was slow...
size_2d, size_3d = self.size_2d, self.size_3d
result = np.zeros(2*size_2d)
# Go over all possible component projections (z, y, x) to (u, v):
@@ -103,7 +106,7 @@ class Projector(object):
return result
def _vector_field_projection_T(self, vector):
- self.LOG.debug('Calling _vector_field_projection_T')
+# self._log.debug('Calling _vector_field_projection_T') # TODO: Profiler says this was slow...
size_2d, size_3d = self.size_2d, self.size_3d
result = np.zeros(3*size_3d)
# Go over all possible component projections (u, v) to (z, y, x):
@@ -122,11 +125,11 @@ class Projector(object):
return result
def _scalar_field_projection(self, vector):
- self.LOG.debug('Calling _scalar_field_projection')
+ self._log.debug('Calling _scalar_field_projection')
return np.array(self.weight.dot(vector))
def _scalar_field_projection_T(self, vector):
- self.LOG.debug('Calling _scalar_field_projection_T')
+ self._log.debug('Calling _scalar_field_projection_T')
return np.array(self.weight.T.dot(vector))
def jac_dot(self, vector):
@@ -145,12 +148,12 @@ class Projector(object):
always`size_2d`.
'''
- self.LOG.debug('Calling jac_dot')
+# self._log.debug('Calling jac_dot') # TODO: Profiler says this was slow...
if len(vector) == 3*self.size_3d: # mode == 'vector'
- self.LOG.debug('mode == vector')
+# self._log.debug('mode == vector') # TODO: Profiler says this was slow...
return self._vector_field_projection(vector)
elif len(vector) == self.size_3d: # mode == 'scalar'
- self.LOG.debug('mode == scalar')
+# self._log.debug('mode == scalar') # TODO: Profiler says this was slow...
return self._scalar_field_projection(vector)
else:
raise AssertionError('Vector size has to be suited either for '
@@ -172,12 +175,12 @@ class Projector(object):
of the :class:`~.Projector` object.
'''
- self.LOG.debug('Calling jac_T_dot')
+# self._log.debug('Calling jac_T_dot') # TODO: Profiler says this was slow...
if len(vector) == 2*self.size_2d: # mode == 'vector'
- self.LOG.debug('mode == vector')
+# self._log.debug('mode == vector') # TODO: Profiler says this was slow...
return self._vector_field_projection_T(vector)
elif len(vector) == self.size_2d: # mode == 'scalar'
- self.LOG.debug('mode == scalar')
+# self._log.debug('mode == scalar') # TODO: Profiler says this was slow...
return self._scalar_field_projection_T(vector)
else:
raise AssertionError('Vector size has to be suited either for '
@@ -219,22 +222,22 @@ class XTiltProjector(Projector):
'''
- LOG = logging.getLogger(__name__+'.XTiltProjector')
+ _log = logging.getLogger(__name__+'.XTiltProjector')
def __init__(self, dim, tilt, dim_uv=None):
def get_position(p, m, b, size):
- self.LOG.debug('Calling get_position')
+ self._log.debug('Calling get_position')
y, x = np.array(p)[:, 0]+0.5, np.array(p)[:, 1]+0.5
return (y-m*x-b)/np.sqrt(m**2+1) + size/2.
def get_impact(pos, r, size):
- self.LOG.debug('Calling get_impact')
+ self._log.debug('Calling get_impact')
return [x for x in np.arange(np.floor(pos-r), np.floor(pos+r)+1, dtype=int)
if 0 <= x < size]
def get_weight(delta, rho): # use circles to represent the voxels
- self.LOG.debug('Calling get_weight')
+ self._log.debug('Calling get_weight')
lo, up = delta-rho, delta+rho
# Upper boundary:
if up >= 1:
@@ -248,7 +251,7 @@ class XTiltProjector(Projector):
w_lo = (lo*np.sqrt(1-lo**2) + np.arctan(lo/np.sqrt(1-lo**2))) / pi
return w_up - w_lo
- self.LOG.debug('Calling __init__')
+ self._log.debug('Calling __init__')
self.tilt = tilt
# Set starting variables:
# length along projection (proj, z), perpendicular (perp, y) and rotation (rot, x) axis:
@@ -292,7 +295,7 @@ class XTiltProjector(Projector):
shape=(size_2d, size_3d)))
coeff = [[1, 0, 0], [0, np.cos(tilt), np.sin(tilt)]]
super(XTiltProjector, self).__init__(dim, dim_uv, weight, coeff)
- self.LOG.debug('Created '+str(self))
+ self._log.debug('Created '+str(self))
def get_info(self):
'''Get specific information about the projector as a string.
@@ -329,7 +332,7 @@ class YTiltProjector(Projector):
'''
- LOG = logging.getLogger(__name__+'.YTiltProjector')
+ _log = logging.getLogger(__name__+'.YTiltProjector')
def __init__(self, dim, tilt, dim_uv=None):
@@ -355,7 +358,7 @@ class YTiltProjector(Projector):
w_lo = (lo*np.sqrt(1-lo**2) + np.arctan(lo/np.sqrt(1-lo**2))) / pi
return w_up - w_lo
- self.LOG.debug('Calling __init__')
+ self._log.debug('Calling __init__')
self.tilt = tilt
# Set starting variables:
# length along projection (proj, z), rotation (rot, y) and perpendicular (perp, x) axis:
@@ -399,7 +402,7 @@ class YTiltProjector(Projector):
shape=(size_2d, size_3d)))
coeff = [[np.cos(tilt), 0, np.sin(tilt)], [0, 1, 0]]
super(YTiltProjector, self).__init__(dim, dim_uv, weight, coeff)
- self.LOG.debug('Created '+str(self))
+ self._log.debug('Created '+str(self))
def get_info(self):
'''Get specific information about the projector as a string.
@@ -437,31 +440,34 @@ class SimpleProjector(Projector):
'''
- LOG = logging.getLogger(__name__+'.SimpleProjector')
+ _log = logging.getLogger(__name__+'.SimpleProjector')
AXIS_DICT = {'z': (0, 1, 2), 'y': (1, 0, 2), 'x': (2, 1, 0)} # (0:z, 1:y, 2:x) -> (proj, v, u)
def __init__(self, dim, axis='z', dim_uv=None):
- self.LOG.debug('Calling __init__')
+ self._log.debug('Calling __init__')
assert axis in {'z', 'y', 'x'}, 'Projection axis has to be x, y or z (given as a string)!'
proj, v, u = self.AXIS_DICT[axis]
- dim_proj, dim_v, dim_u = dim[proj], dim[v], dim[u]
+ if axis=='x':
+ dim_proj, dim_v, dim_u = dim[proj], dim[u], dim[v] # coordinate switch for 'x'!
+ else:
+ dim_proj, dim_v, dim_u = dim[proj], dim[v], dim[u]
dim_z, dim_y, dim_x = dim
size_2d = dim_u * dim_v
size_3d = np.prod(dim)
data = np.repeat(1, size_3d) # size_3d ones in the matrix (each voxel is projected)
indptr = np.arange(0, size_3d+1, dim_proj) # each row has dim_proj ones
if axis == 'z':
- self.LOG.debug('Projecting along the z-axis')
+ self._log.debug('Projecting along the z-axis')
coeff = [[1, 0, 0], [0, 1, 0]]
indices = np.array([np.arange(row, size_3d, size_2d)
for row in range(size_2d)]).reshape(-1)
elif axis == 'y':
- self.LOG.debug('Projection along the y-axis')
+ self._log.debug('Projection along the y-axis')
coeff = [[1, 0, 0], [0, 0, 1]]
indices = np.array([np.arange(row%dim_x, dim_x*dim_y, dim_x)+int(row/dim_x)*dim_x*dim_y
for row in range(size_2d)]).reshape(-1)
elif axis == 'x':
- self.LOG.debug('Projection along the x-axis')
+ self._log.debug('Projection along the x-axis')
coeff = [[0, 0, 1], [0, 1, 0]] # Caution, coordinate switch: u, v --> z, y (not y, z!)
# indices = np.array([np.arange(dim_proj) + row*dim_proj
# for row in range(size_2d)]).reshape(-1) # this is u, v --> y, z
@@ -484,7 +490,7 @@ class SimpleProjector(Projector):
# Create weight-matrix:
weight = csr_matrix((data, indices, indptr), shape=(size_2d, size_3d))
super(SimpleProjector, self).__init__(dim, dim_uv, weight, coeff)
- self.LOG.debug('Created '+str(self))
+ self._log.debug('Created '+str(self))
def get_info(self):
'''Get specific information about the projector as a string.
diff --git a/pyramid/reconstruction.py b/pyramid/reconstruction.py
index bdbb818..62812db 100644
--- a/pyramid/reconstruction.py
+++ b/pyramid/reconstruction.py
@@ -18,14 +18,14 @@ from pyramid.phasemapper import PhaseMapperRDFC
from pyramid.costfunction import Costfunction
from pyramid.magdata import MagData
-from jutil import cg, minimizer
-
from scipy.optimize import leastsq
import logging
-LOG = logging.getLogger(__name__)
+__all__ = ['optimize_linear', 'optimize_nonlin', 'optimize_splitbregman',
+ 'optimize_simple_leastsq']
+_log = logging.getLogger(__name__)
class PrintIterator(object):
@@ -50,18 +50,18 @@ class PrintIterator(object):
'''
- LOG = logging.getLogger(__name__ + '.PrintIterator')
+ _log = logging.getLogger(__name__ + '.PrintIterator')
def __init__(self, cost, verbosity):
- self.LOG.debug('Calling __init__')
+ self._log.debug('Calling __init__')
self.cost = cost
self.verbosity = verbosity
assert verbosity in {0, 1, 2}, 'verbosity has to be set to 0, 1 or 2!'
self.iteration = 0
- self.LOG.debug('Created ' + str(self))
+ self._log.debug('Created ' + str(self))
def __call__(self, xk):
- self.LOG.debug('Calling __call__')
+ self._log.debug('Calling __call__')
if self.verbosity == 0:
return
print 'iteration #', self.next(),
@@ -71,11 +71,11 @@ class PrintIterator(object):
print ''
def __repr__(self):
- self.LOG.debug('Calling __repr__')
+ self._log.debug('Calling __repr__')
return '%s(cost=%r, verbosity=%r)' % (self.__class__, self.cost, self.verbosity)
def __str__(self):
- self.LOG.debug('Calling __str__')
+ self._log.debug('Calling __str__')
return 'PrintIterator(cost=%s, verbosity=%s)' % (self.cost, self.verbosity)
def next(self):
@@ -83,7 +83,7 @@ class PrintIterator(object):
return self.iteration
-def optimize_linear(data, regularisator=None, max_iter=None):
+def optimize_linear(data, regularisator=None, max_iter=None, info=None):
'''Reconstruct a three-dimensional magnetic distribution from given phase maps via the
conjugate gradient optimizaion method :func:`~.scipy.sparse.linalg.cg`.
Blazingly fast for l2-based cost functions.
@@ -105,15 +105,17 @@ def optimize_linear(data, regularisator=None, max_iter=None):
'''
import jutil.cg as jcg
- LOG.debug('Calling optimize_linear')
+ _log.debug('Calling optimize_linear')
# Set up necessary objects:
cost = Costfunction(data, regularisator)
- LOG.info('Cost before optimization: {}'.format(cost(np.zeros(cost.n))))
+ _log.info('Cost before optimization: {}'.format(cost(np.zeros(cost.n))))
x_opt = jcg.conj_grad_minimize(cost, max_iter=max_iter).x
- LOG.info('Cost after optimization: {}'.format(cost(x_opt)))
+ _log.info('Cost after optimization: {}'.format(cost(x_opt)))
# Create and return fitting MagData object:
mag_opt = MagData(data.a, np.zeros((3,) + data.dim))
mag_opt.set_vector(x_opt, data.mask)
+ if info is not None:
+ info[:] = cost.chisq, cost.chisq_m, cost.chisq_a
return mag_opt
@@ -141,7 +143,7 @@ def optimize_nonlin(data, first_guess=None, regularisator=None):
'''
import jutil.minimizer as jmin
import jutil.norms as jnorms
- LOG.debug('Calling optimize_nonlin')
+ _log.debug('Calling optimize_nonlin')
if first_guess is None:
first_guess = MagData(data.a, np.zeros((3,) + data.dim))
@@ -160,14 +162,14 @@ def optimize_nonlin(data, first_guess=None, regularisator=None):
return direc_p
# This Method is semi-best for Lp type problems. Takes forever, though
- LOG.info('Cost before optimization: {}'.format(cost(np.zeros(cost.n))))
+ _log.info('Cost before optimization: {}'.format(cost(np.zeros(cost.n))))
result = jmin.minimize(
cost, x_0,
method="SteepestDescent",
options={"preconditioner": preconditioner},
tol={"max_iteration": 10000})
x_opt = result.x
- LOG.info('Cost after optimization: {}'.format(cost(x_opt)))
+ _log.info('Cost after optimization: {}'.format(cost(x_opt)))
mag_opt = MagData(data.a, np.zeros((3,) + data.dim))
mag_opt.set_vector(x_opt, data.mask)
return mag_opt
@@ -203,7 +205,7 @@ def optimize_splitbregman(data, weight, lam, mu):
import jutil.operator as joperator
import jutil.diff as jdiff
from pyramid.regularisator import FirstOrderRegularisator
- LOG.debug('Calling optimize_splitbregman')
+ _log.debug('Calling optimize_splitbregman')
# regularisator is actually not necessary, but this makes the cost
# function to that which is supposedly optimized by split bregman.
diff --git a/pyramid/regularisator.py b/pyramid/regularisator.py
index f34713c..03cdaaa 100644
--- a/pyramid/regularisator.py
+++ b/pyramid/regularisator.py
@@ -10,15 +10,15 @@ import abc
import numpy as np
-from scipy.sparse import coo_matrix, csr_matrix
import jutil.norms as jnorm
import jutil.diff as jdiff
import jutil.operator as joperator
-from pyramid.converter import IndexConverter
import logging
+__all__ = ['NoneRegularisator', 'ZeroOrderRegularisator', 'FirstOrderRegularisator']
+
# TODO: Fragen für Jörn: Macht es Sinn, x_a standardmäßig auf den Nullvektor zu setzen? Ansonsten
# besser im jeweiligen Konstruktor setzen, nicht im abstrakten!
# Wie kommt man genau an die Ableitungen (Normen sind nicht unproblematisch)?
@@ -29,40 +29,40 @@ class Regularisator(object):
# TODO: Docstring!
__metaclass__ = abc.ABCMeta
- LOG = logging.getLogger(__name__+'.Regularisator')
+ _log = logging.getLogger(__name__+'.Regularisator')
@abc.abstractmethod
def __init__(self, norm, lam):
- self.LOG.debug('Calling __init__')
+ self._log.debug('Calling __init__')
self.norm = norm
self.lam = lam
- self.LOG.debug('Created '+str(self))
+ self._log.debug('Created '+str(self))
def __call__(self, x):
- self.LOG.debug('Calling __call__')
+ self._log.debug('Calling __call__')
return self.lam * self.norm(x)
def __repr__(self):
- self.LOG.debug('Calling __repr__')
+ self._log.debug('Calling __repr__')
return '%s(norm=%r, lam=%r)' % (self.__class__, self.norm, self.lam)
def __str__(self):
- self.LOG.debug('Calling __str__')
+ self._log.debug('Calling __str__')
return 'Regularisator(norm=%s, lam=%s)' % (self.norm, self.lam)
def jac(self, x):
# TODO: Docstring!
- self.LOG.debug('Calling jac')
+ self._log.debug('Calling jac')
return self.lam * self.norm.jac(x)
def hess_dot(self, x, vector):
# TODO: Docstring!
- self.LOG.debug('Calling hess_dot')
+# self._log.debug('Calling hess_dot') # TODO: Profiler says this was slow...
return self.lam * self.norm.hess_dot(x, vector)
def hess_diag(self, x, vector):
# TODO: Docstring!
- self.LOG.debug('Calling hess_diag')
+ self._log.debug('Calling hess_diag')
return self.lam * self.norm.hess_diag(x, vector)
@@ -74,28 +74,28 @@ class NoneRegularisator(Regularisator):
LOG = logging.getLogger(__name__+'.NoneRegularisator')
def __init__(self):
- self.LOG.debug('Calling __init__')
+ self._log.debug('Calling __init__')
self.norm = None
self.lam = 0
- self.LOG.debug('Created '+str(self))
+ self._log.debug('Created '+str(self))
def __call__(self, x):
- self.LOG.debug('Calling __call__')
+ self._log.debug('Calling __call__')
return 0
def jac(self, x):
# TODO: Docstring!
- self.LOG.debug('Calling jac')
+ self._log.debug('Calling jac')
return np.zeros_like(x)
def hess_dot(self, x, vector):
# TODO: Docstring!
- self.LOG.debug('Calling hess_dot')
+ self._log.debug('Calling hess_dot')
return np.zeros_like(vector)
def hess_diag(self, x, vector):
# TODO: Docstring!
- self.LOG.debug('Calling hess_diag')
+ self._log.debug('Calling hess_diag')
return np.zeros_like(vector)
@@ -105,14 +105,14 @@ class ZeroOrderRegularisator(Regularisator):
LOG = logging.getLogger(__name__+'.ZeroOrderRegularisator')
def __init__(self, _, lam, p=2):
- self.LOG.debug('Calling __init__')
+ self._log.debug('Calling __init__')
self.p = p
if p == 2:
norm = jnorm.L2Square()
else:
norm = jnorm.LPPow(p, 1e-12)
super(ZeroOrderRegularisator, self).__init__(norm, lam)
- self.LOG.debug('Created '+str(self))
+ self._log.debug('Created '+str(self))
class FirstOrderRegularisator(Regularisator):
@@ -122,10 +122,11 @@ class FirstOrderRegularisator(Regularisator):
self.p = p
D0 = jdiff.get_diff_operator(mask, 0, 3)
D1 = jdiff.get_diff_operator(mask, 1, 3)
- D = joperator.VStack([D0, D1])
+ D2 = jdiff.get_diff_operator(mask, 2, 3)
+ D = joperator.VStack([D0, D1, D2])
if p == 2:
norm = jnorm.WeightedL2Square(D)
else:
norm = jnorm.WeightedTV(jnorm.LPPow(p, 1e-12), D, [D0.shape[0], D.shape[0]])
super(FirstOrderRegularisator, self).__init__(norm, lam)
- self.LOG.debug('Created '+str(self))
+ self._log.debug('Created '+str(self))
diff --git a/pyramid/version.py b/pyramid/version.py
new file mode 100644
index 0000000..049ee68
--- /dev/null
+++ b/pyramid/version.py
@@ -0,0 +1,3 @@
+# THIS FILE IS GENERATED FROM THE PYRAMID SETUP.PY
+version="0.1.0-dev"
+hg_revision="aa1763114c1c+"
diff --git a/scripts/collaborations/bielefeld.py b/scripts/collaborations/bielefeld.py
index cc2b5bb..74206c3 100644
--- a/scripts/collaborations/bielefeld.py
+++ b/scripts/collaborations/bielefeld.py
@@ -1,6 +1,7 @@
# -*- coding: utf-8 -*-
"""Created on Thu Oct 02 11:53:25 2014 @author: Jan"""
+from pyramid import *
import os
@@ -31,41 +32,41 @@ def load_from_llg(filename):
logging.config.fileConfig(LOGGING_CONF, disable_existing_loggers=False)
-## 2DCoFe:
-mag_data_2DCoFe = MagData.load_from_llg(PATH+'magnetic_distribution_2DCoFe.txt')
-mag_data_2DCoFe.quiver_plot(proj_axis='x')
-plt.savefig(PATH+'magnetic_distribution_2DCoFe.png')
-mag_data_2DCoFe.quiver_plot3d()
-phase_map_2DCoFe = pm(mag_data_2DCoFe, axis='x')
-phase_map_2DCoFe.display_combined(gain='auto')
-# AH21FIN:
-mag_data_AH21FIN = MagData.load_from_llg(PATH+'magnetic_distribution_AH21FIN.txt')
-mag_data_AH21FIN.quiver_plot(proj_axis='x')
-plt.savefig(PATH+'magnetic_distribution_AH21FIN.png')
-mag_data_AH21FIN.quiver_plot3d()
-phase_map_AH21FIN = pm(mag_data_AH21FIN, axis='x')
-phase_map_AH21FIN.display_combined(gain='auto')
-# AH41FIN:
-mag_data_AH41FIN = MagData.load_from_llg(PATH+'magnetic_distribution_AH41FIN.txt')
-mag_data_AH41FIN.quiver_plot(proj_axis='x')
-plt.savefig(PATH+'magnetic_distribution_AH41FIN.png')
-mag_data_AH41FIN.quiver_plot3d()
-phase_map_AH41FIN = pm(mag_data_AH41FIN, axis='x')
-phase_map_AH41FIN.display_combined(gain='auto')
-# Chain:
-mag_data_chain = load_from_llg(PATH+'magnetic_distribution_Chain.txt')
-mag_data_chain.quiver_plot(proj_axis='x')
-plt.savefig(PATH+'magnetic_distribution_Chain.png')
-mag_data_chain.quiver_plot3d()
-phase_map_chain = pm(mag_data_chain, axis='x')
-phase_map_chain.display_combined(gain='auto')
-# Cylinder:
-mag_data_cyl = load_from_llg(PATH+'magnetic_distribution_Cylinder.txt')
-mag_data_cyl.quiver_plot(proj_axis='z')
-plt.savefig(PATH+'magnetic_distribution_Cylinder.png')
-mag_data_cyl.quiver_plot3d()
-phase_map_cyl = pm(mag_data_cyl, axis='z')
-phase_map_cyl.display_combined(gain='auto')
+# 2DCoFe:
+#mag_data_2DCoFe = MagData.load_from_llg(PATH+'magnetic_distribution_2DCoFe.txt')
+#mag_data_2DCoFe.quiver_plot(proj_axis='x')
+#plt.savefig(PATH+'magnetic_distribution_2DCoFe.png')
+#mag_data_2DCoFe.quiver_plot3d()
+#phase_map_2DCoFe = pm(mag_data_2DCoFe, axis='x')
+#phase_map_2DCoFe.display_combined(gain='auto')
+## AH21FIN:
+#mag_data_AH21FIN = MagData.load_from_llg(PATH+'magnetic_distribution_AH21FIN.txt')
+#mag_data_AH21FIN.quiver_plot(proj_axis='x')
+#plt.savefig(PATH+'magnetic_distribution_AH21FIN.png')
+#mag_data_AH21FIN.quiver_plot3d()
+#phase_map_AH21FIN = pm(mag_data_AH21FIN, axis='x')
+#phase_map_AH21FIN.display_combined(gain='auto')
+## AH41FIN:
+#mag_data_AH41FIN = MagData.load_from_llg(PATH+'magnetic_distribution_AH41FIN.txt')
+#mag_data_AH41FIN.quiver_plot(proj_axis='x')
+#plt.savefig(PATH+'magnetic_distribution_AH41FIN.png')
+#mag_data_AH41FIN.quiver_plot3d()
+#phase_map_AH41FIN = pm(mag_data_AH41FIN, axis='x')
+#phase_map_AH41FIN.display_combined(gain='auto')
+## Chain:
+#mag_data_chain = load_from_llg(PATH+'magnetic_distribution_Chain.txt')
+#mag_data_chain.quiver_plot(proj_axis='x')
+#plt.savefig(PATH+'magnetic_distribution_Chain.png')
+#mag_data_chain.quiver_plot3d()
+#phase_map_chain = pm(mag_data_chain, axis='x')
+#phase_map_chain.display_combined(gain='auto')
+## Cylinder:
+#mag_data_cyl = load_from_llg(PATH+'magnetic_distribution_Cylinder.txt')
+#mag_data_cyl.quiver_plot(proj_axis='z')
+#plt.savefig(PATH+'magnetic_distribution_Cylinder.png')
+#mag_data_cyl.quiver_plot3d()
+#phase_map_cyl = pm(mag_data_cyl, axis='z')
+#phase_map_cyl.display_combined(gain='auto')
# Ring:
mag_data_ring = load_from_llg(PATH+'magnetic_distribution_ring.txt')
mag_data_ring.quiver_plot(proj_axis='x')
diff --git a/scripts/collaborations/example_joern.py b/scripts/collaborations/example_joern.py
index a80e2a7..4fae246 100644
--- a/scripts/collaborations/example_joern.py
+++ b/scripts/collaborations/example_joern.py
@@ -68,8 +68,8 @@ log = True
PATH = '../../output/joern/'
###################################################################################################
# Read in files:
-phase_map = PhaseMap.load_from_netcdf4(PATH+'phase_map_2.nc')
-with open(PATH + 'mask_2.pickle') as pf:
+phase_map = PhaseMap.load_from_netcdf4(PATH+'phase_map.nc')
+with open(PATH + 'mask.pickle') as pf:
mask = pickle.load(pf)
dim = mask.shape
# Setup:
@@ -88,7 +88,7 @@ else:
tic = clock()
if p == 2:
- mag_data_rec = rc.optimize_linear(data_set, regularisator=regularisator)
+ mag_data_rec = rc.optimize_linear(data_set, regularisator=regularisator, max_iter=50)
else:
print regularisator.p
mag_data_rec = rc.optimize_nonlin(data_set, regularisator=regularisator)
@@ -101,36 +101,33 @@ with open(dirname + "result.pickle", "wb") as pf:
print 'reconstruction time:', clock() - tic
# Display the reconstructed phase map and holography image:
phase_map_rec = pm(mag_data_rec)
-phase_map_rec.display_combined('Reconstr. Distribution', gain=gain,
- interpolation=inter, show=False)
+phase_map_rec.display_combined('Reconstr. Distribution', gain=gain, interpolation=inter)
plt.savefig(dirname + "/reconstr.png")
# Plot the magnetization:
-axis = (mag_data_rec*(1/mag_data_rec.magnitude.max())).quiver_plot(show=False)
-axis.set_xlim(int(20/64*dim[1], 45/64*dim[2])
-axis.set_ylim(int(20/64*dim[1], 45/64*dim[2])
+axis = (mag_data_rec*(1/mag_data_rec.magnitude.max())).quiver_plot()
+axis.set_xlim(20./64*dim[1], 45/64*dim[2])
+axis.set_ylim(20./64*dim[1], 45/64*dim[2])
plt.savefig(dirname + "/quiver.png")
# Display the Phase:
phase_diff = phase_map_rec-phase_map
-phase_diff.display_phase('Difference', show=False)
+phase_diff.display_phase('Difference')
plt.savefig(dirname + "/difference.png")
# Get the average difference from the experimental results:
print 'Average difference:', np.average(phase_diff.phase)
# Plot holographic contour maps with overlayed magnetic distributions:
-axis = phase_map_rec.display_holo('Magnetization Overlay', gain=0.1,
- interpolation=inter, show=False)
-mag_data_rec.quiver_plot(axis=axis, show=False)
+axis = phase_map_rec.display_holo('Magnetization Overlay', gain=0.1, interpolation=inter)
+mag_data_rec.quiver_plot(axis=axis)
axis = plt.gca()
-axis.set_xlim(int(20/64*dim[1], 45/64*dim[2])
-axis.set_ylim(int(20/64*dim[1], 45/64*dim[2])
+axis.set_xlim(20./64*dim[1], 45/64*dim[2])
+axis.set_ylim(20./64*dim[1], 45/64*dim[2])
plt.savefig(dirname + "/overlay_normal.png")
-axis = phase_map_rec.display_holo('Magnetization Overlay', gain=0.1,
- interpolation=inter, show=False)
-mag_data_rec.quiver_plot(axis=axis, log=log, show=False)
+axis = phase_map_rec.display_holo('Magnetization Overlay', gain=0.1, interpolation=inter)
+mag_data_rec.quiver_plot(axis=axis, log=log)
axis = plt.gca()
-axis.set_xlim(int(20/64*dim[1], 45/64*dim[2])
-axis.set_ylim(int(20/64*dim[1], 45/64*dim[2])
+axis.set_xlim(20./64*dim[1], 45/64*dim[2])
+axis.set_ylim(20./64*dim[1], 45/64*dim[2])
plt.savefig(dirname + "/overlay_log.png")
diff --git a/scripts/collaborations/patrick_nanowire_reconstruction.py b/scripts/collaborations/patrick_nanowire_reconstruction.py
new file mode 100644
index 0000000..03f31ed
--- /dev/null
+++ b/scripts/collaborations/patrick_nanowire_reconstruction.py
@@ -0,0 +1,56 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Fri Nov 14 14:17:56 2014
+
+@author: Jan
+"""
+
+from PIL import Image
+import numpy as np
+
+from jutil.taketime import TakeTime
+from pyramid import *
+
+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 = 1.0 # in nm
+gain = 5
+b_0 = 1
+lam = 1E-6
+PATH = '../../output/patrick/'
+dim_uv = (128, 128)
+###################################################################################################
+
+# Read in files:
+im_mask = Image.open(PATH+'min102mask.tif')
+im_mask.thumbnail(dim_uv, Image.ANTIALIAS)
+mask = np.array(np.array(im_mask)/255, dtype=bool)
+mask = np.expand_dims(mask, axis=0)
+im_phase = Image.open(PATH+'min102.tif')
+im_phase.thumbnail(dim_uv, Image.ANTIALIAS)
+phase = np.array(im_phase)/255. - 0.5#*(18.19977+8.057761) - 8.057761
+phase_map = PhaseMap(a, phase)
+dim = mask.shape
+
+data_set = DataSet(a, dim, b_0, mask)
+data_set.append(phase_map, SimpleProjector(dim))
+
+regularisator = FirstOrderRegularisator(mask, lam, p=2)
+
+with TakeTime('reconstruction time'):
+ mag_data_rec = rc.optimize_linear(data_set, regularisator=regularisator, max_iter=50)
+
+# Display the reconstructed phase map and holography image:
+phase_map.display_combined('Input PhaseMap', gain=40)
+mag_data_rec.quiver_plot(log=True)
+phase_map_rec = pm(mag_data_rec)
+phase_map_rec.display_combined('Reconstr. Distribution', gain=40)
+#plt.savefig(dirname + "/reconstr.png")
diff --git a/scripts/collaborations/test.py b/scripts/collaborations/test.py
deleted file mode 100644
index f65c3bb..0000000
--- a/scripts/collaborations/test.py
+++ /dev/null
@@ -1,29 +0,0 @@
-# -*- coding: utf-8 -*-
-"""
-Created on Mon Oct 27 08:49:48 2014
-
-@author: Jan
-"""
-
-
-from numpy import pi
-import numpy as np
-
-import pyramid.magcreator as mc
-from pyramid.magdata import MagData
-
-
-dim = (1, 32, 32)
-center = (0, 16, 16)
-width = (1, 8, 16)
-mag_data = MagData(1, mc.create_mag_dist_homog(mc.Shapes.slab(dim, center, width), pi/6))
-mag_data.quiver_plot()
-mag_data.quiver_plot3d()
-x, y, z = mag_data.magnitude[:, 0, ...]
-x_rot = np.rot90(mag_data.magnitude[0, 0, ...]).copy()
-y_rot = np.rot90(mag_data.magnitude[1, 0, ...]).copy()
-#z_rot = np.rot90(mag_data.magnitude[2, 0, ...]).copy()
-mag_data.magnitude[0, 0, ...] = -y_rot
-mag_data.magnitude[1, 0, ...] = x_rot
-mag_data.quiver_plot()
-mag_data.quiver_plot3d()
diff --git a/scripts/collaborations/vtk_conversion_nanowire_full_90deg.py b/scripts/collaborations/vtk_conversion_nanowire_full_90deg.py
new file mode 100644
index 0000000..acf17b6
--- /dev/null
+++ b/scripts/collaborations/vtk_conversion_nanowire_full_90deg.py
@@ -0,0 +1,180 @@
+# -*- coding: utf-8 -*-
+"""Created on Fri Jan 24 11:17:11 2014 @author: Jan"""
+
+
+import os
+
+import numpy as np
+import matplotlib.pyplot as plt
+from pylab import griddata
+
+import pickle
+import vtk
+from tqdm import tqdm
+
+import pyramid
+from pyramid.magdata import MagData
+from pyramid.projector import YTiltProjector
+from pyramid.kernel import Kernel
+from pyramid.phasemapper import PhaseMapperRDFC
+from pyramid.phasemap import PhaseMap
+
+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/longtube_withcap/CoFeB_tube_cap_4nm'
+b_0 = 1.54
+force_calculation = False
+count = 16
+dim_uv_x = (500, 100)
+dim_uv_y = (100, 500)
+gain = 8
+###################################################################################################
+# Load vtk-data:
+if force_calculation or not os.path.exists(PATH+'.pickle'):
+ # Setting up reader:
+ reader = vtk.vtkDataSetReader()
+ reader.SetFileName(PATH+'.vtk')
+ reader.ReadAllScalarsOn()
+ reader.ReadAllVectorsOn()
+ reader.Update()
+ # Getting output:
+ output = reader.GetOutput()
+ # Reading points and vectors:
+ size = output.GetNumberOfPoints()
+ vtk_points = output.GetPoints().GetData()
+ vtk_vectors = output.GetPointData().GetVectors()
+ # Converting points to numpy array:
+ point_array = np.zeros(vtk_points.GetSize())
+ vtk_points.ExportToVoidPointer(point_array)
+ point_array = np.reshape(point_array, (-1, 3))
+ # Converting vectors to numpy array:
+ vector_array = np.zeros(vtk_points.GetSize())
+ vtk_vectors.ExportToVoidPointer(vector_array)
+ vector_array = np.reshape(vector_array, (-1, 3))
+ # Combining data:
+ data = np.hstack((point_array, vector_array))
+ with open(PATH+'.pickle', 'w') as pf:
+ pickle.dump(data, pf)
+else:
+ with open(PATH+'.pickle') as pf:
+ data = pickle.load(pf)
+###################################################################################################
+# Interpolate on regular grid:
+if force_calculation or not os.path.exists(PATH+'.nc'):
+ # Determine the size of the object:
+ x_min, x_max = data[:, 0].min(), data[:, 0].max()
+ y_min, y_max = data[:, 1].min(), data[:, 1].max()
+ z_min, z_max = data[:, 2].min(), data[:, 2].max()
+ 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.
+ # Filling zeros:
+ iz_x = np.concatenate([np.linspace(-2.95, -2.95, 20),
+ np.linspace(-2.95, 0, 20),
+ np.linspace(0, 2.95, 20),
+ np.linspace(2.95, 2.95, 20),
+ np.linspace(2.95, 0, 20),
+ np.linspace(0, -2.95, 20), ])
+ iz_y = np.concatenate([np.linspace(-1.70, 1.70, 20),
+ np.linspace(1.70, 3.45, 20),
+ np.linspace(3.45, 1.70, 20),
+ np.linspace(1.70, -1.70, 20),
+ np.linspace(-1.70, -3.45, 20),
+ np.linspace(-3.45, -1.70, 20), ])
+ # Find unique z-slices:
+ zs_unique = np.unique(data[:, 2])
+ # Determine the grid spacing (in 1/10 nm):
+ a = (zs_unique[1] - zs_unique[0])
+ # Set actual z-slices:
+ zs = np.arange(z_min, z_max, a)
+ # Create regular grid:
+ xs = np.arange(x_cent-x_diff, x_cent+x_diff, a)
+ ys = np.arange(y_cent-y_diff, y_cent+y_diff, a)
+ xx, yy = np.meshgrid(xs, ys)
+ # Create empty magnitude:
+ magnitude = np.zeros((3, len(zs), len(ys), len(xs)))
+ # Go over all slices:
+ for i, z in tqdm(enumerate(zs), total=len(zs)):
+ z_slice = data[np.abs(data[:, 2]-z) <= a/2., :]
+ weights = 1 - np.abs(z_slice[:, 2]-z)*2/a
+ for j in range(3): # For all 3 components!
+ if z <= 0:
+ grid_x = np.concatenate([z_slice[:, 0], iz_x])
+ grid_y = np.concatenate([z_slice[:, 1], iz_y])
+ grid_z = np.concatenate([weights*z_slice[:, 3+j], np.zeros(len(iz_x))])
+ else:
+ grid_x = z_slice[:, 0]
+ grid_y = z_slice[:, 1]
+ grid_z = weights*z_slice[:, 3+j]
+ grid = griddata(grid_x, grid_y, grid_z, xx, yy)
+ magnitude[j, i, :, :] = grid.filled(fill_value=0)
+ a = a*10 # convert to nm
+ # Creating MagData object and save it as netcdf-file:
+ mag_data = MagData(a, np.pad(magnitude, ((0, 0), (0, 0), (0, 1), (0, 1)), mode='constant',
+ constant_values=0))
+ mag_data.save_to_netcdf4(PATH+'.nc')
+else:
+ mag_data = MagData.load_from_netcdf4(PATH+'.nc')
+###################################################################################################
+# Close ups:
+dim = (72, 300, 72)
+dim_uv = (600, 150)
+angles = [0, 10, 20, 30, 40, 50, 60]
+# Turn magnetization around by 90° around x-axis:
+magnitude_new = np.zeros((3, mag_data.dim[1], mag_data.dim[0], mag_data.dim[2]))
+for i in range(mag_data.magnitude.shape[2]):
+ x_rot = np.rot90(mag_data.magnitude[0, ..., i]).copy()
+ y_rot = np.rot90(mag_data.magnitude[1, ..., i]).copy()
+ z_rot = np.rot90(mag_data.magnitude[2, ..., i]).copy()
+ magnitude_new[0, ..., i] = x_rot
+ magnitude_new[1, ..., i] = -z_rot
+ magnitude_new[2, ..., i] = y_rot
+mag_data.magnitude = magnitude_new
+# Iterate over all angles:
+for angle in angles:
+ angle_rad = angle * np.pi/180
+ projector = YTiltProjector(dim, angle_rad, dim_uv)
+ projector_scaled = YTiltProjector((dim[0]/2, dim[1]/2, dim[2]/2), angle_rad,
+ (dim_uv[0]/2, dim_uv[1]/2))
+ # Tip:
+ mag_data_tip = MagData(mag_data.a, mag_data.magnitude[:, :, 608:, :])
+ PM = PhaseMapperRDFC(Kernel(mag_data.a, projector.dim_uv))
+ phase_map_tip = PhaseMap(mag_data.a, PM(projector(mag_data_tip)).phase[350:530, :])
+ phase_map_tip.display_combined('Phase Map Nanowire Tip', gain=gain,
+ interpolation='bilinear')
+ plt.savefig(PATH+'_nanowire_tip_xtilt_{}.png'.format(angle))
+ mag_data_tip.scale_down()
+ mag_proj_tip = projector_scaled(mag_data_tip)
+ axis = mag_proj_tip.quiver_plot()
+# axis.set_xlim(17, 55)
+# axis.set_ylim(180-shift/2, 240-shift/2)
+ plt.savefig(PATH+'_nanowire_tip_mag_xtilt_{}.png'.format(angle))
+ axis = mag_proj_tip.quiver_plot(log=True)
+# axis.set_xlim(17, 55)
+# axis.set_ylim(180-shift/2, 240-shift/2)
+ plt.savefig(PATH+'_nanowire_tip_mag_log_xtilt_{}.png'.format(angle))
+ # Bottom:
+ mag_data_bot = MagData(mag_data.a, mag_data.magnitude[:, :, :300, :])
+ PM = PhaseMapperRDFC(Kernel(mag_data.a, projector.dim_uv))
+ phase_map_tip = PhaseMap(mag_data.a, PM(projector(mag_data_bot)).phase[50:225, :])
+ phase_map_tip.display_combined('Phase Map Nanowire Bottom', gain=gain,
+ interpolation='bilinear')
+ plt.savefig(PATH+'_nanowire_bot_xtilt_{}_no_frame.png'.format(angle))
+ mag_data_bot.scale_down()
+ mag_proj_bot = projector_scaled(mag_data_bot)
+ axis = mag_proj_bot.quiver_plot()
+# axis.set_xlim(17, 55)
+# axis.set_ylim(50+shift/2, 110+shift/2)
+ plt.savefig(PATH+'_nanowire_bot_mag_xtilt_{}.png'.format(angle))
+ axis = mag_proj_bot.quiver_plot(log=True)
+# axis.set_xlim(17, 55)
+# axis.set_ylim(50+shift/2, 110+shift/2)
+ plt.savefig(PATH+'_nanowire_bot_mag_log_xtilt_{}.png'.format(angle))
+ # Close plots:
+ plt.close('all')
diff --git a/scripts/gui/mag_slicer.py b/scripts/gui/mag_slicer.py
index 353f91d..1ded00a 100644
--- a/scripts/gui/mag_slicer.py
+++ b/scripts/gui/mag_slicer.py
@@ -209,13 +209,13 @@ class UI_MagSlicerMain(QtGui.QWidget):
self.update_slice()
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)
+ self.phase_map.display_phase(axis=self.mplWidgetPhase.axes, cbar=False)
if self.checkBoxSmooth.isChecked():
interpolation = 'bilinear'
else:
interpolation = 'none'
self.phase_map.display_holo(axis=self.mplWidgetHolo.axes, gain=gain,
- interpolation=interpolation, show=False)
+ interpolation=interpolation)
self.mplWidgetPhase.draw()
self.mplWidgetHolo.draw()
@@ -224,8 +224,7 @@ class UI_MagSlicerMain(QtGui.QWidget):
self.mag_data.quiver_plot(axis=self.mplWidgetMag.axes, proj_axis=self.mode,
ax_slice=self.spinBoxSlice.value(),
log=self.checkBoxLog.isChecked(),
- scaled=self.checkBoxScale.isChecked(),
- show=False)
+ scaled=self.checkBoxScale.isChecked())
self.mplWidgetMag.draw()
def load(self):
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 bb8ecf0..383e4c6 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
@@ -19,6 +19,7 @@ from pyramid.kernel import Kernel
import time
import shelve
+import gc
import logging
import logging.config
@@ -59,8 +60,9 @@ key = 'ch5-4-method_comparison_new'
if key in data_shelve and not force_calculation:
print '--LOAD METHOD DATA'
(data_disc_fourier0, data_vort_fourier0,
+ data_disc_fourier1, data_vort_fourier1,
data_disc_fourier3, data_vort_fourier3,
- data_disc_fourier10, data_vort_fourier10,
+ data_disc_fourier7, data_vort_fourier7,
data_disc_real, data_vort_real) = data_shelve[key]
else:
# Input parameters:
@@ -84,10 +86,12 @@ else:
dim_list = [dim[2]/2**i for i in range(steps)]
data_disc_fourier0 = np.vstack((dim_list, np.zeros((3, steps))))
data_vort_fourier0 = np.vstack((dim_list, np.zeros((3, steps))))
+ data_disc_fourier1 = np.vstack((dim_list, np.zeros((3, steps))))
+ data_vort_fourier1 = np.vstack((dim_list, np.zeros((3, steps))))
data_disc_fourier3 = np.vstack((dim_list, np.zeros((3, steps))))
data_vort_fourier3 = np.vstack((dim_list, np.zeros((3, steps))))
- data_disc_fourier10 = np.vstack((dim_list, np.zeros((3, steps))))
- data_vort_fourier10 = np.vstack((dim_list, np.zeros((3, steps))))
+ data_disc_fourier7 = np.vstack((dim_list, np.zeros((3, steps))))
+ data_vort_fourier7 = np.vstack((dim_list, np.zeros((3, steps))))
data_disc_real = np.vstack((dim_list, np.zeros((3, steps))))
data_vort_real = np.vstack((dim_list, np.zeros((3, steps))))
@@ -104,8 +108,9 @@ else:
# Create projector along z-axis and phasemapper:
projector = SimpleProjector(dim)
pm_fourier0 = PhaseMapperFDFC(a, projector.dim_uv, padding=0)
+ pm_fourier1 = PhaseMapperFDFC(a, projector.dim_uv, padding=1)
pm_fourier3 = PhaseMapperFDFC(a, projector.dim_uv, padding=3)
- pm_fourier10 = PhaseMapperFDFC(a, projector.dim_uv, padding=7)
+ pm_fourier7 = PhaseMapperFDFC(a, projector.dim_uv, padding=7)
pm_real = PhaseMapperRDFC(Kernel(a, projector.dim_uv))
print '----CALCULATE RMS/DURATION HOMOG. MAGN. DISC'
@@ -118,6 +123,13 @@ else:
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 1:
+ phase_num_disc, t, dt = get_time(pm_fourier1, projector, mag_data_disc, n)
+ data_disc_fourier1[2, i], data_disc_fourier1[3, i] = t, dt
+ print '------time (disc, fourier3) =', data_disc_fourier1[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_fourier1[1, i] = np.sqrt(np.mean(phase_diff_disc.phase**2))
# Fourier padding 3:
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
@@ -126,12 +138,12 @@ else:
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, 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_num_disc, t, dt = get_time(pm_fourier7, projector, mag_data_disc, n)
+ data_disc_fourier7[2, i], data_disc_fourier7[3, i] = t, dt
+ print '------time (disc, fourier10) =', data_disc_fourier7[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))
+ data_disc_fourier7[1, i] = np.sqrt(np.mean(phase_diff_disc.phase**2))
# Real space disc:
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
@@ -155,6 +167,15 @@ else:
print 'phase mean:', np.mean(phase_num_vort.phase)
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 1:
+ phase_num_vort, t, dt = get_time(pm_fourier1, projector, mag_data_vort, n)
+ phase_fft3 = phase_num_vort
+ data_vort_fourier1[2, i], data_vort_fourier1[3, i] = t, dt
+ print '------time (vortex, fourier3) =', data_vort_fourier1[2, i]
+ phase_diff_vort = (phase_ana_vort-phase_num_vort) * 1E3 # in mrad -> *1000
+ print 'phase mean:', np.mean(phase_num_vort.phase)
+ phase_diff_vort -= np.mean(phase_diff_vort.phase)
+ data_vort_fourier1[1, i] = np.sqrt(np.mean(phase_diff_vort.phase**2))
# Fourier padding 3:
phase_num_vort, t, dt = get_time(pm_fourier3, projector, mag_data_vort, n)
phase_fft3 = phase_num_vort
@@ -165,14 +186,14 @@ 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, projector, mag_data_vort, n)
+ phase_num_vort, t, dt = get_time(pm_fourier7, 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]
+ data_vort_fourier7[2, i], data_vort_fourier7[3, i] = t, dt
+ print '------time (vortex, fourier10) =', data_vort_fourier7[2, i]
phase_diff_vort = (phase_ana_vort-phase_num_vort) * 1E3 # in mrad -> *1000
print 'phase mean:', np.mean(phase_num_vort.phase)
phase_diff_vort -= np.mean(phase_diff_vort.phase)
- data_vort_fourier10[1, i] = np.sqrt(np.mean(phase_diff_vort.phase**2))
+ data_vort_fourier7[1, i] = np.sqrt(np.mean(phase_diff_vort.phase**2))
# Real space disc:
phase_num_vort, t, dt = get_time(pm_real, projector, mag_data_vort, n)
phase_real = phase_num_vort
@@ -189,11 +210,14 @@ else:
# Scale down for next iteration:
mag_data_disc.scale_down()
mag_data_vort.scale_down()
+ # Force garbage collection:
+ gc.collect()
print '--SHELVE METHOD DATA'
data_shelve[key] = (data_disc_fourier0, data_vort_fourier0,
+ data_disc_fourier1, data_vort_fourier1,
data_disc_fourier3, data_vort_fourier3,
- data_disc_fourier10, data_vort_fourier10,
+ data_disc_fourier7, data_vort_fourier7,
data_disc_real, data_vort_real)
print '--PLOT/SAVE METHOD DATA'
@@ -207,8 +231,9 @@ fig.suptitle('Method Comparison', fontsize=20)
axes[1, 0].set_xscale('log')
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], ':bd')
axes[1, 0].plot(data_disc_fourier3[0], data_disc_fourier3[1], ':bo')
-axes[1, 0].plot(data_disc_fourier10[0], data_disc_fourier10[1], ':b^')
+axes[1, 0].plot(data_disc_fourier7[0], data_disc_fourier7[1], ':b^')
axes[1, 0].plot(data_disc_real[0], data_disc_real[1], '--ro')
axes[1, 0].set_xlabel('grid size [px]', fontsize=15)
axes[1, 0].set_ylabel('RMS [mrad]', fontsize=15)
@@ -225,12 +250,14 @@ axes[0, 0].set_xscale('log')
axes[0, 0].set_yscale('log')
axes[0, 0].errorbar(data_disc_fourier0[0], data_disc_fourier0[2],
yerr=0*data_vort_fourier3[3], fmt=':bs')
+axes[0, 0].errorbar(data_disc_fourier1[0], data_disc_fourier1[2],
+ yerr=0*data_vort_fourier1[3], fmt=':bd')
axes[0, 0].errorbar(data_disc_fourier3[0], data_disc_fourier3[2],
yerr=0*data_vort_fourier3[3], fmt=':bo')
-axes[0, 0].errorbar(data_disc_fourier10[0], data_disc_fourier10[2],
- yerr=0*data_vort_fourier3[3], fmt=':b^')
+axes[0, 0].errorbar(data_disc_fourier7[0], data_disc_fourier7[2],
+ yerr=0*data_vort_fourier7[3], fmt=':b^')
axes[0, 0].errorbar(data_disc_real[0], data_disc_real[2],
- yerr=0*data_vort_fourier3[3], fmt='--ro')
+ yerr=0*data_disc_real[3], fmt='--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(25, 350)
@@ -246,10 +273,12 @@ axes[1, 1].set_xscale('log')
axes[1, 1].set_yscale('log')
axes[1, 1].plot(data_vort_fourier0[0], data_vort_fourier0[1], ':bs',
label='Fourier padding=0')
+axes[1, 1].plot(data_vort_fourier1[0], data_vort_fourier1[1], ':bd',
+ label='Fourier padding=1')
axes[1, 1].plot(data_vort_fourier3[0], data_vort_fourier3[1], ':bo',
label='Fourier padding=3')
-axes[1, 1].plot(data_vort_fourier10[0], data_vort_fourier10[1], ':b^',
- label='Fourier padding=10')
+axes[1, 1].plot(data_vort_fourier7[0], data_vort_fourier7[1], ':b^',
+ label='Fourier padding=7')
axes[1, 1].plot(data_vort_real[0], data_vort_real[1], '--ro',
label='Real space')
axes[1, 1].set_xlabel('grid size [px]', fontsize=15)
@@ -267,13 +296,15 @@ axes[0, 1].errorbar(data_vort_fourier0[0], data_vort_fourier0[2], yerr=0*data_vo
fmt=':bs', label='Fourier padding=0')
axes[0, 1].errorbar(data_vort_fourier3[0], data_vort_fourier3[2], yerr=0*data_vort_fourier3[3],
fmt=':bo', label='Fourier padding=3')
-axes[0, 1].errorbar(data_vort_fourier10[0], data_vort_fourier10[2], yerr=0*data_vort_fourier3[3],
- fmt=':b^', label='Fourier padding=10')
+axes[0, 1].errorbar(data_vort_fourier1[0], data_vort_fourier1[2], yerr=0*data_vort_fourier3[3],
+ fmt=':bd', label='Fourier padding=1')
+axes[0, 1].errorbar(data_vort_fourier7[0], data_vort_fourier7[2], yerr=0*data_vort_fourier3[3],
+ fmt=':b^', label='Fourier padding=7')
axes[0, 1].errorbar(data_vort_real[0], data_vort_real[2], yerr=0*data_vort_fourier3[3],
fmt='--ro', label='Real space')
axes[0, 1].set_title('Vortex state disc', fontsize=18)
axes[0, 1].set_xlim(25, 350)
-axes[0, 1].set_ylim(1E-4, 1E1)
+axes[0, 1].set_ylim(5E-4, 5E0)
axes[0, 1].tick_params(axis='both', which='major', labelsize=14)
axes[0, 1].xaxis.set_major_locator(LogLocator(base=2))
axes[0, 1].xaxis.set_major_formatter(LogFormatter(base=2))
diff --git a/scripts/reconstruction/reconstruction_sparse_cg_test.py b/scripts/reconstruction/reconstruction_linear_test.py
similarity index 59%
rename from scripts/reconstruction/reconstruction_sparse_cg_test.py
rename to scripts/reconstruction/reconstruction_linear_test.py
index 48e1471..2b29a26 100644
--- a/scripts/reconstruction/reconstruction_sparse_cg_test.py
+++ b/scripts/reconstruction/reconstruction_linear_test.py
@@ -7,17 +7,17 @@ import os
import numpy as np
from numpy import pi
-from time import clock
-
import pyramid
from pyramid.magdata import MagData
from pyramid.phasemap import PhaseMap
from pyramid.projector import YTiltProjector, XTiltProjector
from pyramid.dataset import DataSet
-from pyramid.regularisator import ZeroOrderRegularisator
+from pyramid.regularisator import FirstOrderRegularisator
import pyramid.magcreator as mc
import pyramid.reconstruction as rc
+from jutil.taketime import TakeTime
+
import logging
import logging.config
@@ -32,23 +32,28 @@ print('--Generating input phase_maps')
a = 10.
b_0 = 1.
-dim = (8, 32, 32)
+dim = (32, 32, 32)
dim_uv = dim[1:3]
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)
+lam = 1E-4
+center = (dim[0]/2-0.5, int(dim[1]/2)-0.5, int(dim[2]/2)-0.5)
+radius_core = dim[1]/8
+radius_shell = dim[1]/4
+height = dim[0]/2
+# Create magnetic shape:
+mag_shape_core = mc.Shapes.disc(dim, center, radius_core, height)
+mag_shape_outer = mc.Shapes.disc(dim, center, radius_shell, height)
+mag_shape_shell = np.logical_xor(mag_shape_outer, mag_shape_core)
+mag_data = MagData(a, mc.create_mag_dist_vortex(mag_shape_shell, magnitude=0.75))
+mag_data += MagData(a, mc.create_mag_dist_homog(mag_shape_core, phi=0, theta=0))
+
+shape = mc.Shapes.disc(dim, center, radius_shell, height)
+magnitude = mc.create_mag_dist_vortex(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))
-
-mag_data.quiver_plot3d()
+#mag_data.quiver_plot('z-projection', proj_axis='z')
+#mag_data.quiver_plot('y-projection', proj_axis='y')
+#mag_data.quiver_plot3d('Original distribution')
tilts_full = np.linspace(-pi/2, pi/2, num=count/2, endpoint=False)
tilts_miss = np.linspace(-pi/3, pi/3, num=count/2, endpoint=False)
@@ -65,35 +70,31 @@ projectors_xy_miss.extend(projectors_x_miss)
projectors_xy_miss.extend(projectors_y_miss)
###################################################################################################
-projectors = projectors_xy_full
-noise = 0.0
+projectors = projectors_xy_miss
+noise = 0
###################################################################################################
print('--Setting up data collection')
-data = DataSet(a, dim, b_0)
+mask = mag_data.get_mask()
+data = DataSet(a, dim, b_0, mask)
data.projectors = projectors
data.phase_maps = data.create_phase_maps(mag_data)
if noise != 0:
- for phase_map in data.phase_maps:
+ for i, phase_map in enumerate(data.phase_maps):
phase_map += PhaseMap(a, np.random.normal(0, noise, dim_uv))
+ data.phase_maps[i] = phase_map
###################################################################################################
-print('--Test simple solver')
+print('--Reconstruction')
-lam = 1E-4
-reg = ZeroOrderRegularisator(lam)
-start = clock()
-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()
+reg = FirstOrderRegularisator(mask, lam, p=2)
+with TakeTime('reconstruction'):
+ mag_data_opt = rc.optimize_linear(data, regularisator=reg, max_iter=100)
###################################################################################################
print('--Plot stuff')
-phase_maps_opt = data.create_phase_maps(mag_data_opt)
-#data.display_phase()
-#data.display_phase(phase_maps_opt)
-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]
+mag_data_opt.quiver_plot3d('Reconstructed distribution')
+#(mag_data_opt - mag_data).quiver_plot3d('Difference')
+#phase_maps_opt = data.create_phase_maps(mag_data_opt)
diff --git a/scripts/reconstruction/reconstruction_linear_test_batch.py b/scripts/reconstruction/reconstruction_linear_test_batch.py
new file mode 100644
index 0000000..38a25fc
--- /dev/null
+++ b/scripts/reconstruction/reconstruction_linear_test_batch.py
@@ -0,0 +1,148 @@
+# -*- coding: utf-8 -*-
+"""Created on Fri Feb 28 14:25:59 2014 @author: Jan"""
+
+
+import os
+
+import numpy as np
+from numpy import pi
+
+import itertools
+
+from mayavi import mlab
+
+import pyramid
+from pyramid.magdata import MagData
+from pyramid.phasemap import PhaseMap
+from pyramid.projector import YTiltProjector, XTiltProjector
+from pyramid.dataset import DataSet
+from pyramid.regularisator import ZeroOrderRegularisator, FirstOrderRegularisator
+import pyramid.magcreator as mc
+import pyramid.reconstruction as rc
+
+from jutil.taketime import TakeTime
+
+import psutil
+import gc
+
+import shelve
+
+import logging
+import logging.config
+
+
+LOGGING_CONF = os.path.join(os.path.dirname(os.path.realpath(pyramid.__file__)), 'logging.ini')
+PATH = '../../output/3d reconstruction/'
+
+logging.config.fileConfig(LOGGING_CONF, disable_existing_loggers=False)
+if not os.path.exists(PATH):
+ os.mkdir(PATH)
+proc = psutil.Process(os.getpid())
+###################################################################################################
+print('--Constant parameters')
+
+a = 10.
+b_0 = 1.
+
+###################################################################################################
+print('--Magnetization distributions')
+
+dim = (32, 32, 32)
+center = (dim[0]/2-0.5, int(dim[1]/2)-0.5, int(dim[2]/2)-0.5)
+radius_core = dim[1]/8
+radius_shell = dim[1]/4
+height = dim[0]/2
+# Vortex:
+mag_shape_disc = mc.Shapes.disc(dim, center, radius_shell, height)
+mag_data_vortex = MagData(a, mc.create_mag_dist_vortex(mag_shape_disc))
+# Sphere:
+mag_shape_sphere = mc.Shapes.sphere(dim, center, radius_shell)
+mag_data_sphere = MagData(a, mc.create_mag_dist_homog(mag_shape_sphere, phi=pi/4, theta=pi/4))
+# Core-Shell:
+mag_shape_core = mc.Shapes.disc(dim, center, radius_core, height)
+mag_shape_shell = np.logical_xor(mag_shape_disc, mag_shape_core)
+mag_data_core_shell = MagData(a, mc.create_mag_dist_vortex(mag_shape_shell, magnitude=0.75))
+mag_data_core_shell += MagData(a, mc.create_mag_dist_homog(mag_shape_core, phi=0, theta=0))
+
+# TODO: cubic bar magnet, rectangular bar magnet, -sphere
+
+###################################################################################################
+print('--Set up configurations')
+
+mag_datas = {'vortex_disc': mag_data_vortex, 'homog_sphere': mag_data_sphere,
+ 'core_shell': mag_data_core_shell}
+masks = {'mask': True}
+xy_tilts = {'xy': [True, True]}
+max_tilts = [60]
+tilt_steps = [10]
+noises = [0]
+orders = [1]
+lambdas = [1E-4]
+# Combining everything:
+config_list = [mag_datas.keys(), masks.keys(), xy_tilts.keys(),
+ max_tilts, tilt_steps, noises, orders, lambdas]
+
+###################################################################################################
+print('--Original distributions')
+
+for key in mag_datas:
+ mag_datas[key].quiver_plot3d(title='Original distribution ({})'.format(key))
+ mlab.savefig('{}{}_analytic.png'.format(PATH, key), size=(800, 800))
+ mlab.close(all=True)
+
+###################################################################################################
+print('--Reconstruction')
+
+print 'Number of configurations:', len(list(itertools.product(*config_list)))
+
+for configuration in itertools.product(*config_list):
+ # Extract keys:
+ mag_key, mask_key, xy_key, max_tilt, tilt_step, noise, order, lam = configuration
+ print '----CONFIG:', configuration
+ name = '{}_{}_{}tilts_max{}_in{}steps_{}noise_{}order_{}lam'.format(*configuration)
+ dim = mag_datas[mag_key].dim
+ # Mask:
+ if masks[mask_key]:
+ mask = mag_datas[mag_key].get_mask()
+ else:
+ mask = np.ones(dim, dtype=bool)
+ # DataSet:
+ data = DataSet(a, dim, b_0, mask)
+ # Tilts:
+ tilts = np.arange(-max_tilt/180.*pi, max_tilt/180.*pi, tilt_step/180.*pi)
+ # Projectors:
+ projectors = []
+ if xy_tilts[xy_key][0]:
+ projectors.extend([XTiltProjector(dim, tilt) for tilt in tilts])
+ if xy_tilts[xy_key][1]:
+ projectors.extend([YTiltProjector(dim, tilt) for tilt in tilts])
+ data.projectors = projectors
+ # PhaseMaps:
+ data.phase_maps = data.create_phase_maps(mag_datas[mag_key])
+ # Noise:
+ if noise != 0:
+ for i, phase_map in enumerate(data.phase_maps):
+ phase_map += PhaseMap(a, np.random.normal(0, noise, data.projectors[i].dim_uv))
+ data.phase_maps[i] = phase_map
+ # Regularisation:
+ if order == 0:
+ reg = ZeroOrderRegularisator(mask, lam, p=2)
+ if order == 1:
+ reg = FirstOrderRegularisator(mask, lam, p=2)
+ # Reconstruction:
+ info = []
+ with TakeTime('reconstruction'):
+ mag_data_rec = rc.optimize_linear(data, regularisator=reg, max_iter=100, info=info)
+ # Plots:
+ mag_data_rec.save_to_netcdf4('{}{}.nc'.format(PATH, name))
+ mag_data_rec.quiver_plot3d('Reconstructed distribution ({})'.format(mag_key))
+ mlab.savefig('{}{}_REC.png'.format(PATH, name), size=(800, 800))
+ (mag_data_rec - mag_datas[mag_key]).quiver_plot3d('Difference ({})'.format(mag_key))
+ mlab.savefig('{}{}_DIFF.png'.format(PATH, name), size=(800, 800))
+ mlab.close(all=True)
+ data_shelve = shelve.open(PATH+'/3d_shelve')
+ data_shelve[name] = info
+ data_shelve.close()
+ print 'chisq = {:.6f}, chisq_m = {:.6f}, chisq_a = {:.6f}'.format(*info)
+ gc.collect()
+ print 'RSS = {:.2f} MB'.format(proc.memory_info().rss/1024.**2)
diff --git a/scripts/reconstruction/reconstruction_sparse_cg_singular_values.py b/scripts/reconstruction/reconstruction_sparse_cg_singular_values.py
index a581395..0ccee36 100644
--- a/scripts/reconstruction/reconstruction_sparse_cg_singular_values.py
+++ b/scripts/reconstruction/reconstruction_sparse_cg_singular_values.py
@@ -11,7 +11,6 @@ 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
@@ -63,7 +62,7 @@ y = data.phase_vec
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))
+# MTM = M.T * M + lam * np.asmatrix(np.eye(3*size_3d))
U, s, V = np.linalg.svd(M) # TODO: M or MTM?
diff --git a/scripts/test methods/compare_kernels.py b/scripts/test methods/compare_kernels.py
new file mode 100644
index 0000000..9270e18
--- /dev/null
+++ b/scripts/test methods/compare_kernels.py
@@ -0,0 +1,121 @@
+#! python
+# -*- coding: utf-8 -*-
+"""Compare the phase map of one pixel for different approaches."""
+
+
+import os
+
+import numpy as np
+
+import pyramid
+import pyramid.magcreator as mc
+from pyramid.projector import SimpleProjector
+from pyramid.phasemapper import PhaseMapperRDFC, PhaseMapperFDFC
+from pyramid.kernel import Kernel
+from pyramid.magdata import MagData
+from pyramid.phasemap import PhaseMap
+
+import matplotlib.pyplot as plt
+from matplotlib.ticker import IndexLocator
+
+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)
+directory = '../../output/magnetic distributions'
+if not os.path.exists(directory):
+ os.makedirs(directory)
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+# Input parameters:
+a = 1.0 # in nm
+phi = 0 # in rad
+dim = (1, 32, 32)
+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
+ coeff = - 1j * b_0 * a**2 / (2*PHI_0)
+ nyq = 1 / a # nyquist frequency
+ 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)
+ # 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')
+projector = SimpleProjector(dim)
+# Kernel of a disc in real space:
+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 = 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 = 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:
+phase_map_fft_kernel = PhaseMap(a, get_fourier_kernel(), 'mrad')
+phase_map_fft_kernel.display_phase('Phase of one Pixel (Fourier Kernel)', limit=limit)
+# Kernel differences:
+print 'Fourier Kernel, direct and indirect method are identical:', \
+ np.all(phase_map_fft_kernel.phase - phase_map_fft.phase) == 0
+(phase_map_disc-phase_map_fft).display_phase('Phase difference of one Pixel (Disc - Fourier)')
+
+# Cross section plots of real space kernels:
+fig = plt.figure()
+axis = fig.add_subplot(1, 1, 1)
+x = np.linspace(-dim[1]/a/2, dim[1]/a/2-1, dim[1])
+y_ft = phase_map_fft.phase[:, dim[1]/2]
+y_re = phase_map_disc.phase[:, dim[1]/2]
+axis.axhline(0, color='k')
+axis.plot(x, y_re, label='Real space method')
+axis.plot(x, y_ft, label='Fourier method')
+axis.grid()
+axis.legend()
+axis.set_title('Real Space Kernel')
+axis.set_xlim(-dim[1]/2, dim[1]/2-1)
+axis.xaxis.set_major_locator(IndexLocator(base=dim[1]/8, offset=0))
+
+# Cross section plots of Fourier space kernels:
+fig = plt.figure()
+axis = fig.add_subplot(1, 1, 1)
+x = range(dim[1])
+k_re = np.abs(np.fft.ifftshift(np.fft.rfft2(phase_map_disc.phase), axes=0))[:, 0]**2
+k_ft = np.abs(np.fft.ifftshift(np.fft.rfft2(phase_map_fft.phase), axes=0))[:, 0]**2
+axis.axhline(0, color='k')
+axis.plot(x, k_re, label='Real space method')
+axis.plot(x, k_ft, label='Fourier method')
+axis.grid()
+axis.legend()
+axis.set_title('Fourier Space Kernel')
+axis.set_xlim(0, dim[1]-1)
+axis.xaxis.set_major_locator(IndexLocator(base=dim[1]/8, offset=0))
diff --git a/scripts/test methods/compare_methods.py b/scripts/test methods/compare_methods.py
index 2c07e5a..a164ac9 100644
--- a/scripts/test methods/compare_methods.py
+++ b/scripts/test methods/compare_methods.py
@@ -31,12 +31,12 @@ b_0 = 1.0 # in T
a = 1.0 # in nm
phi = pi/4
numcore = False
-padding = 10
-gain = 10
-dim = (128, 128, 128) # in px (z, y, x)
+padding = 1
+gain = 'auto'
+dim = (32, 256, 256) # in px (z, y, x)
# Create magnetic shape:
-geometry = 'sphere'
+geometry = 'disc'
if geometry == 'slab':
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!
width = (dim[0]/2, dim[1]/2, dim[2]/2) # in px (z, y, x)
@@ -71,6 +71,9 @@ start_time = time.time()
phase_map_conv = pm_conv(projection)
print 'Time for RDFC: ', time.time() - start_time
start_time = time.time()
+phase_map_conv = pm_conv(projection)
+print 'Time for RDFC: ', time.time() - start_time
+start_time = time.time()
phase_map_four = pm_four(projection)
print 'Time for FDFC: ', time.time() - start_time
diff --git a/scripts/test methods/fftw_test.py b/scripts/test methods/fftw_test.py
new file mode 100644
index 0000000..f49de53
--- /dev/null
+++ b/scripts/test methods/fftw_test.py
@@ -0,0 +1,55 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Mon Nov 17 14:28:34 2014
+
+@author: Jan
+"""
+
+
+import os
+
+import numpy as np
+import pyfftw
+
+import pyramid
+from pyramid import *
+from jutil.taketime import TakeTime
+
+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)
+
+n = 2**20
+print n
+
+a = pyfftw.n_byte_align_empty(n, pyfftw.simd_alignment, 'complex128')
+a[:] = np.random.randn(n) + 1j*np.random.randn(n)
+pyfftw.interfaces.cache.enable()
+
+with TakeTime('FFTW uncached'):
+ b = pyfftw.interfaces.numpy_fft.fft(a)
+with TakeTime('FFTW cached '):
+ b = pyfftw.interfaces.numpy_fft.fft(a)
+with TakeTime('Numpy '):
+ c = np.fft.fft(a)
+d = pyfftw.n_byte_align_empty(n, pyfftw.simd_alignment, 'complex128')
+fft_object = pyfftw.FFTW(a, d)
+with TakeTime('FFTW-object '):
+ d = fft_object(a, d)
+fft_object = pyfftw.builders.fft(a, threads=1)
+with TakeTime('FFTW-build 1t'):
+ for i in range(100):
+ d = fft_object(a, d)
+fft_object = pyfftw.builders.fft(a, threads=4)
+with TakeTime('FFTW-build 4t'):
+ for i in range(100):
+ d = fft_object(a, d)
+fft_object = pyfftw.builders.fft(a, threads=8)
+with TakeTime('FFTW-build 8t'):
+ for i in range(100):
+ d = fft_object(a, d)
diff --git a/setup.py b/setup.py
index b10aaee..962fbe5 100644
--- a/setup.py
+++ b/setup.py
@@ -1,40 +1,78 @@
-# -*- coding: utf-8 -*-
+#!python
"""Setup for testing, building, distributing and installing the 'Pyramid'-package"""
-import numpy
import os
-import sys
-import sysconfig
import subprocess
+import sys
+import re
+import numpy
from distutils.command.build import build
-from Cython.Distutils import build_ext
-from setuptools import setup
-from setuptools import find_packages
+from setuptools import setup, find_packages
from setuptools.extension import Extension
-
-from pyramid._version import __version__
-
-
-class custom_build(build):
- '''Custom build command'''
-
- def make_hgrevision(self, target):
- output = subprocess.Popen(["hg", "id", "-i"], stdout=subprocess.PIPE).communicate()[0]
- hgrevision_cc = file(str(target), "w")
- hgrevision_cc.write('hg_revision = "{0}"'.format(output.strip()))
- hgrevision_cc.close()
-
- def run(self):
- build.run(self)
- print 'creating hg_revision.txt'
- self.make_hgrevision(os.path.join('build', get_build_path('lib'), 'hg_revision.txt'))
+from Cython.Distutils import build_ext
-def get_build_path(dname):
- '''Returns the name of a distutils build directory'''
- path = "{dirname}.{platform}-{version[0]}.{version[1]}"
- return path.format(dirname=dname, platform=sysconfig.get_platform(), version=sys.version_info)
+DISTNAME = 'pyramid'
+DESCRIPTION = 'PYthon based Reconstruction Algorithm for MagnetIc Distributions'
+LONG_DESCRIPTION = 'long description (TODO!)' # TODO: Long description!
+MAINTAINER = 'Jan Caron'
+MAINTAINER_EMAIL = 'j.caron@fz-juelich.de'
+URL = ''
+VERSION = '0.1.0-dev'
+PYTHON_VERSION = (2, 7)
+DEPENDENCIES = {'numpy': (1, 6), 'cython': (0, 20)}
+
+
+def get_package_version(package):
+ version = []
+ for version_attr in ('version', 'VERSION', '__version__'):
+ if (hasattr(package, version_attr) and
+ isinstance(getattr(package, version_attr), str)):
+ version_info = getattr(package, version_attr, '')
+ for part in re.split('\D+', version_info):
+ try:
+ version.append(int(part))
+ except ValueError:
+ pass
+ return tuple(version)
+
+
+def check_requirements():
+ if sys.version_info < PYTHON_VERSION:
+ raise SystemExit('You need Python version %d.%d or later.'
+ % PYTHON_VERSION)
+
+ for package_name, min_version in DEPENDENCIES.items():
+ dep_error = False
+ try:
+ package = __import__(package_name)
+ except ImportError:
+ dep_error = True
+ else:
+ package_version = get_package_version(package)
+ if min_version > package_version:
+ dep_error = True
+
+ if dep_error:
+ raise ImportError('You need `%s` version %d.%d or later.'
+ % ((package_name, ) + min_version))
+
+
+def hg_version():
+ try:
+ hg_rev = subprocess.check_output(['hg', 'id', '--id']).strip()
+ except:
+ hg_rev = "???"
+ return hg_rev
+
+
+def write_version_py(filename='pyramid/version.py'):
+ version_string = "# THIS FILE IS GENERATED FROM THE PYRAMID SETUP.PY\n" + \
+ 'version="{}"\n'.format(VERSION) + \
+ 'hg_revision="{}"\n'.format(hg_version())
+ with open(os.path.join(os.path.dirname(__file__), filename), 'w') as vfile:
+ vfile.write(version_string)
def get_files(rootdir):
@@ -48,32 +86,30 @@ def get_files(rootdir):
print '\n-------------------------------------------------------------------------------'
-
-setup(
- name = 'Pyramid',
- version = __version__,
- description = 'PYthon based Reconstruction Algorithm for MagnetIc Distributions',
- author = 'Jan Caron',
- author_email = 'j.caron@fz-juelich.de',
- url = 'fz-juelich.de',
-
- packages = find_packages(exclude=['tests']),
- include_dirs = [numpy.get_include()],
- requires = ['numpy', 'matplotlib', 'mayavi'],
-
- scripts = get_files('scripts'),
-
- test_suite = 'tests',
-
- cmdclass = {'build_ext': build_ext, 'build': custom_build},
-
- ext_package = 'pyramid/numcore',
- ext_modules = [
+print 'checking requirements'
+check_requirements()
+print 'write version.py'
+write_version_py()
+setup(name=DISTNAME,
+ description=DESCRIPTION,
+ long_description=LONG_DESCRIPTION,
+ maintainer=MAINTAINER,
+ maintainer_email=MAINTAINER_EMAIL,
+ url=URL,
+ download_url=URL,
+ version=VERSION,
+ packages=find_packages(exclude=['tests']),
+ include_dirs=[numpy.get_include()],
+ requires=['numpy', 'matplotlib', 'mayavi'],
+ scripts=get_files('scripts'),
+ test_suite='tests',
+ cmdclass={'build_ext': build_ext, 'build': build},
+ ext_package='pyramid/numcore',
+ ext_modules=[
Extension('phasemapper_core', ['pyramid/numcore/phasemapper_core.pyx'],
- include_dirs = [numpy.get_include(), numpy.get_numarray_include()],
+ include_dirs=[numpy.get_include(), numpy.get_numarray_include()],
extra_compile_args=['-march=native', '-mtune=native']
)
]
-)
-
+ )
print '-------------------------------------------------------------------------------\n'
diff --git a/tests/test_analytic.py b/tests/test_analytic.py
index ad78b4c..562ee62 100644
--- a/tests/test_analytic.py
+++ b/tests/test_analytic.py
@@ -7,57 +7,50 @@ import unittest
import numpy as np
from numpy import pi
+from numpy.testing import assert_allclose
+
import pyramid.analytic as an
class TestCaseAnalytic(unittest.TestCase):
"""TestCase for the analytic module."""
- def setUp(self):
- self.path = os.path.join(os.path.dirname(os.path.realpath(__file__)), 'test_analytic/')
- self.dim = (4, 4, 4)
- self.res = 10.0
- self.phi = pi/4
- self.center = (self.dim[0]/2-0.5, self.dim[1]/2-0.5, self.dim[2]/2-0.5)
- self.radius = self.dim[2]/4
-
- def tearDown(self):
- self.path = None
- self.dim = None
- self.res = None
- self.phi = None
- self.center = None
- self.radius = None
+ path = os.path.join(os.path.dirname(os.path.realpath(__file__)), 'test_analytic/')
+ dim = (4, 4, 4)
+ a = 10.0
+ phi = pi/4
+ center = (dim[0]/2-0.5, dim[1]/2-0.5, dim[2]/2-0.5)
+ radius = dim[2]/4
def test_phase_mag_slab(self):
'''Test of the phase_mag_slab method.'''
width = (self.dim[0]/2, self.dim[1]/2, self.dim[2]/2)
- phase = an.phase_mag_slab(self.dim, self.res, self.phi, self.center, width)
+ phase = an.phase_mag_slab(self.dim, self.a, self.phi, self.center, width).phase
reference = np.load(os.path.join(self.path, 'ref_phase_slab.npy'))
- np.testing.assert_almost_equal(phase, reference, err_msg='Unexpected behavior in phase_mag_slab()')
+ assert_allclose(phase, reference, err_msg='Unexpected behavior in phase_mag_slab()')
def test_phase_mag_disc(self):
'''Test of the phase_mag_disc method.'''
radius = self.dim[2]/4
height = self.dim[2]/2
- phase = an.phase_mag_disc(self.dim, self.res, self.phi, self.center, radius, height)
+ phase = an.phase_mag_disc(self.dim, self.a, self.phi, self.center, radius, height).phase
reference = np.load(os.path.join(self.path, 'ref_phase_disc.npy'))
- np.testing.assert_almost_equal(phase, reference, err_msg='Unexpected behavior in phase_mag_disc()')
+ assert_allclose(phase, reference, err_msg='Unexpected behavior in phase_mag_disc()')
def test_phase_mag_sphere(self):
'''Test of the phase_mag_sphere method.'''
radius = self.dim[2]/4
- phase = an.phase_mag_sphere(self.dim, self.res, self.phi, self.center, radius)
+ phase = an.phase_mag_sphere(self.dim, self.a, self.phi, self.center, radius).phase
reference = np.load(os.path.join(self.path, 'ref_phase_sphere.npy'))
- np.testing.assert_almost_equal(phase, reference, err_msg='Unexpected behavior in phase_mag_sphere()')
+ assert_allclose(phase, reference, err_msg='Unexpected behavior in phase_mag_sphere()')
def test_phase_mag_vortex(self):
'''Test of the phase_mag_vortex method.'''
radius = self.dim[2]/4
height = self.dim[2]/2
- phase = an.phase_mag_vortex(self.dim, self.res, self.center, radius, height)
+ phase = an.phase_mag_vortex(self.dim, self.a, self.center, radius, height).phase
reference = np.load(os.path.join(self.path, 'ref_phase_vort.npy'))
- np.testing.assert_almost_equal(phase, reference, err_msg='Unexpected behavior in phase_mag_vortex()')
+ assert_allclose(phase, reference, err_msg='Unexpected behavior in phase_mag_vortex()')
if __name__ == '__main__':
suite = unittest.TestLoader().loadTestsFromTestCase(TestCaseAnalytic)
diff --git a/tests/test_compliance.py b/tests/test_compliance.py
index 943c727..6ffe674 100644
--- a/tests/test_compliance.py
+++ b/tests/test_compliance.py
@@ -15,11 +15,7 @@ import pep8
class TestCaseCompliance(unittest.TestCase):
"""TestCase for checking the pep8 compliance of the pyramid package."""
- def setUp(self):
- self.path = os.path.split(os.path.dirname(os.path.realpath(__file__)))[0] # Pyramid dir
-
- def tearDown(self):
- self.path = None
+ path = os.path.split(os.path.dirname(os.path.realpath(__file__)))[0] # Pyramid dir
def get_files_to_check(self, rootdir):
filepaths = []
diff --git a/tests/test_magcreator.py b/tests/test_magcreator.py
index 32ad03b..137b6e6 100644
--- a/tests/test_magcreator.py
+++ b/tests/test_magcreator.py
@@ -7,65 +7,62 @@ import unittest
import numpy as np
from numpy import pi
+from numpy.testing import assert_allclose
import pyramid.magcreator as mc
class TestCaseMagCreator(unittest.TestCase):
- def setUp(self):
- self.path = os.path.join(os.path.dirname(os.path.realpath(__file__)), 'test_magcreator')
-
- def tearDown(self):
- self.path = None
+ path = os.path.join(os.path.dirname(os.path.realpath(__file__)), 'test_magcreator')
def test_shape_slab(self):
test_slab = mc.Shapes.slab((5, 6, 7), (2, 3, 4), (1, 3, 5))
- np.testing.assert_almost_equal(test_slab, np.load(os.path.join(self.path, 'ref_slab.npy')),
- err_msg='Created slab does not match expectation!')
+ assert_allclose(test_slab, np.load(os.path.join(self.path, 'ref_slab.npy')),
+ err_msg='Created slab does not match expectation!')
def test_shape_disc(self):
test_disc_z = mc.Shapes.disc((5, 6, 7), (2, 3, 4), 2, 3, 'z')
test_disc_y = mc.Shapes.disc((5, 6, 7), (2, 3, 4), 2, 3, 'y')
test_disc_x = mc.Shapes.disc((5, 6, 7), (2, 3, 4), 2, 3, 'x')
- np.testing.assert_almost_equal(test_disc_z, np.load(os.path.join(self.path, 'ref_disc_z.npy')),
- err_msg='Created disc in z-direction does not match expectation!')
- np.testing.assert_almost_equal(test_disc_y, np.load(os.path.join(self.path, 'ref_disc_y.npy')),
- err_msg='Created disc in y-direction does not match expectation!')
- np.testing.assert_almost_equal(test_disc_x, np.load(os.path.join(self.path, 'ref_disc_x.npy')),
- err_msg='Created disc in x-direction does not match expectation!')
+ assert_allclose(test_disc_z, np.load(os.path.join(self.path, 'ref_disc_z.npy')),
+ err_msg='Created disc in z-direction does not match expectation!')
+ assert_allclose(test_disc_y, np.load(os.path.join(self.path, 'ref_disc_y.npy')),
+ err_msg='Created disc in y-direction does not match expectation!')
+ assert_allclose(test_disc_x, np.load(os.path.join(self.path, 'ref_disc_x.npy')),
+ err_msg='Created disc in x-direction does not match expectation!')
def test_shape_sphere(self):
test_sphere = mc.Shapes.sphere((5, 6, 7), (2, 3, 4), 2)
- np.testing.assert_almost_equal(test_sphere, np.load(os.path.join(self.path, 'ref_sphere.npy')),
- err_msg='Created sphere does not match expectation!')
+ assert_allclose(test_sphere, np.load(os.path.join(self.path, 'ref_sphere.npy')),
+ err_msg='Created sphere does not match expectation!')
def test_shape_filament(self):
test_filament_z = mc.Shapes.filament((5, 6, 7), (2, 3), 'z')
test_filament_y = mc.Shapes.filament((5, 6, 7), (2, 3), 'y')
test_filament_x = mc.Shapes.filament((5, 6, 7), (2, 3), 'x')
- np.testing.assert_almost_equal(test_filament_z, np.load(os.path.join(self.path, 'ref_fil_z.npy')),
- err_msg='Created filament in z-direction does not match expectation!')
- np.testing.assert_almost_equal(test_filament_y, np.load(os.path.join(self.path, 'ref_fil_y.npy')),
- err_msg='Created filament in y-direction does not match expectation!')
- np.testing.assert_almost_equal(test_filament_x, np.load(os.path.join(self.path, 'ref_fil_x.npy')),
- err_msg='Created filament in x-direction does not match expectation!')
+ assert_allclose(test_filament_z, np.load(os.path.join(self.path, 'ref_fil_z.npy')),
+ err_msg='Created filament in z-direction does not match expectation!')
+ assert_allclose(test_filament_y, np.load(os.path.join(self.path, 'ref_fil_y.npy')),
+ err_msg='Created filament in y-direction does not match expectation!')
+ assert_allclose(test_filament_x, np.load(os.path.join(self.path, 'ref_fil_x.npy')),
+ err_msg='Created filament in x-direction does not match expectation!')
def test_shape_pixel(self):
test_pixel = mc.Shapes.pixel((5, 6, 7), (2, 3, 4))
- np.testing.assert_almost_equal(test_pixel, np.load(os.path.join(self.path, 'ref_pixel.npy')),
- err_msg='Created pixel does not match expectation!')
+ assert_allclose(test_pixel, np.load(os.path.join(self.path, 'ref_pixel.npy')),
+ err_msg='Created pixel does not match expectation!')
def test_create_mag_dist_homog(self):
mag_shape = mc.Shapes.disc((1, 10, 10), (0, 4.5, 4.5), 3, 1)
magnitude = mc.create_mag_dist_homog(mag_shape, pi/4)
- np.testing.assert_almost_equal(magnitude, np.load(os.path.join(self.path, 'ref_mag_disc.npy')),
- err_msg='Created homog. magnetic distribution does not match expectation')
+ assert_allclose(magnitude, np.load(os.path.join(self.path, 'ref_mag_disc.npy')),
+ err_msg='Created homog. magnetic distribution does not match expectation')
def test_create_mag_dist_vortex(self):
mag_shape = mc.Shapes.disc((1, 10, 10), (0, 4.5, 4.5), 3, 1)
magnitude = mc.create_mag_dist_vortex(mag_shape)
- np.testing.assert_almost_equal(magnitude, np.load(os.path.join(self.path, 'ref_mag_vort.npy')),
- err_msg='Created vortex magnetic distribution does not match expectation')
+ assert_allclose(magnitude, np.load(os.path.join(self.path, 'ref_mag_vort.npy')),
+ err_msg='Created vortex magnetic distribution does not match expectation')
if __name__ == '__main__':
diff --git a/tests/test_magcreator/ref_mag_disc.npy b/tests/test_magcreator/ref_mag_disc.npy
index 2c006ceb7fb4edb858c264147b963036d8ba16b7..01089f0367b79419c960a9dbff395ffb625593e0 100644
GIT binary patch
delta 62
zcmV-E0KxyT6R;Do2mzCj2OyKM2uPEl2N1Ks2k-%t&;q!Va0DQekOfGSfCLb;paiG_
UlTZO5lW+n^lVAZ5vw#8c0ngGFC;$Ke
delta 62
zcmV-E0KxyT6R;Do2mzCj2OyKM2uPEl2N1Ks2k-%t&;q!Va0DQekOfGSfCLb;paiG_
UlTZO5lW+n^lVAZ5vw#8c0ngGFC;$Ke
diff --git a/tests/test_magcreator/ref_mag_vort.npy b/tests/test_magcreator/ref_mag_vort.npy
index 99c4abe469e838c999164765761b7111cd85d61c..f337d0c583e30fe588cdfd66b53cf0e13a8f2424 100644
GIT binary patch
delta 25
hcmdlWyg_&a2jgZB_6(-U0-XAr=QFQl-k6ZV1ORF)2m=5B
delta 26
icmdlWyg_&a2jk`jwtS|^-&q%I7H2VJ-ps+6!2|$(9tc4I
--
GitLab