diff --git a/.hgignore b/.hgignore
index 5f7570492ed255fd62287fc2d9bf2abec379c3c8..3e01d21734b5465c83455d03b3fcce7959475063 100644
--- a/.hgignore
+++ b/.hgignore
@@ -10,4 +10,5 @@ syntax: glob
*.c
output*
build*
-Pyramid.egg-info*
\ No newline at end of file
+Pyramid.egg-info*
+scripts/dyn_0090mT_dyn.h5_dat.h5
diff --git a/docs/pyramid.rst b/docs/pyramid.rst
index 8bd3b4c52bfc4d67c3eb53e14e9df3d08ac51b51..9bd1af8547605b4a2a9511206374ca2ea4f50ef1 100644
--- a/docs/pyramid.rst
+++ b/docs/pyramid.rst
@@ -17,10 +17,10 @@ pyramid Package
:show-inheritance:
:special-members:
-:mod:`datacollection` Module
+:mod:`dataset` Module
----------------------------
-.. automodule:: pyramid.datacollection
+.. automodule:: pyramid.dataset
:members:
:show-inheritance:
:special-members:
diff --git a/logfile.log b/logfile.log
index e69de29bb2d1d6434b8b29ae775ad8c2e48c5391..bde3b6ace28465b1ebdc2d2648894b577b8448bb 100644
--- a/logfile.log
+++ b/logfile.log
@@ -0,0 +1,3 @@
+2014-04-10 16:49:29: DEBUG @ <apptools.preferences.preferences>: loading preferences from <C:\Users\Jan\AppData\Roaming\Enthought\mayavi_e3\preferences.ini>
+2014-04-10 16:49:29: DEBUG @ <apptools.preferences.preferences>: loading preferences from <<open file 'C:\\Python27\\lib\\site-packages\\mayavi\\preferences\\preferences.ini', mode 'rb' at 0x04D00F98>>
+2014-04-10 16:49:29: DEBUG @ <apptools.preferences.preferences>: loading preferences from <<open file 'C:\\Python27\\lib\\site-packages\\tvtk\\plugins\\scene\\preferences.ini', mode 'rb' at 0x050F85A0>>
diff --git a/pyramid/__init__.py b/pyramid/__init__.py
index 2316e348fbabb207bbdb7266f8659749e5ac3540..5ffe7ac97a6da715795842500a1f9b13eab31dff 100644
--- a/pyramid/__init__.py
+++ b/pyramid/__init__.py
@@ -17,15 +17,13 @@ phasemap
Class for the storage of phase data.
analytic
Create phase maps for magnetic distributions with analytic solutions.
-datacollection
+dataset
Class for collecting pairs of phase maps and corresponding projectors.
forwardmodel
Class which represents a phase mapping strategy.
costfunction
Class for the evaluation of the cost of a function.
-optimizer
- Provides strategies for optimizing first guess magnetic distributions.
-reconstructor
+reconstruction
Reconstruct magnetic distributions from given phasemaps.
Subpackages
@@ -35,6 +33,8 @@ numcore
"""
+# TODO: logging setup at script level, only logging itself should be done in the package
+# maybe include in startup script!
import logging, logging.config
import os
@@ -43,4 +43,4 @@ import os
LOGGING_CONF = os.path.join(os.path.dirname(os.path.realpath(__file__)), 'logging.ini')
-logging.config.fileConfig(LOGGING_CONF)
+logging.config.fileConfig(LOGGING_CONF, disable_existing_loggers=False)
diff --git a/pyramid/analytic.py b/pyramid/analytic.py
index bf2ffccd69db5206304c463fd5da401519acecb2..5216fcdeef8750c62f899fe07617c739c10385b3 100644
--- a/pyramid/analytic.py
+++ b/pyramid/analytic.py
@@ -13,7 +13,10 @@ from numpy import pi
from pyramid.phasemap import PhaseMap
+import logging
+
+LOG = logging.getLogger(__name__)
PHI_0 = -2067.83 # magnetic flux in T*nm²
@@ -42,7 +45,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')
# Function for the phase:
def phi_mag(x, y):
def F_0(x, y):
@@ -98,6 +101,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')
# Function for the phase:
def phi_mag(x, y):
r = np.hypot(x - x0, y - y0)
@@ -145,6 +149,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')
# Function for the phase:
def phi_mag(x, y):
r = np.hypot(x - x0, y - y0)
@@ -191,6 +196,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')
# Function for the phase:
def phi_mag(x, y):
r = np.hypot(x - x0, y - y0)
diff --git a/pyramid/costfunction.py b/pyramid/costfunction.py
index 811a2a97a615b5b4de29c24ecfb7b5bb7a079eb4..71a60d49433a7e2acc4195fab05b907d51c18ae6 100644
--- a/pyramid/costfunction.py
+++ b/pyramid/costfunction.py
@@ -1,72 +1,156 @@
# -*- coding: utf-8 -*-
-"""
-Created on Fri Dec 13 10:29:11 2013
+"""This module provides the :class:`~.Costfunction` class which represents a strategy to calculate
+the so called `cost` of a threedimensional magnetization distribution."""
-@author: Jan
-"""# TODO: Docstring!
+# TODO: better names for variables (no uppercase, more than one letter)
+
+# TODO: Regularisation Class?
import numpy as np
+from scipy.sparse.linalg import LinearOperator
+from scipy.sparse import eye
+
+import logging
class Costfunction:
- # TODO: Docstring!
-
- def __init__(self, y, F):
- '''TEST DOCSTRING FOR INIT'''
- # TODO: Docstring!
- self.y = y # TODO: get y from phasemaps!
- self.F = F # Forward Model
- self.Se_inv = np.eye(len(y))
+
+ '''Class for calculating the cost of a 3D magnetic distributions in relation to 2D phase maps.
+
+ Represents a strategy for the calculation of the `cost` of a 3D magnetic distribution in
+ relation to two-dimensional phase maps. The `cost` is a measure for the difference of the
+ simulated phase maps from the magnetic distributions to the given set of phase maps.
+ Furthermore this class provides convenient methods for the calculation of the derivative
+ :func:`~.jac` or the product with the Hessian matrix :func:`~.hess_dot` of the costfunction,
+ which can be used by optimizers.
+
+ Attributes
+ ----------
+ y : :class:`~numpy.ndarray` (N=1)
+ Vector which lists all pixel values of all phase maps one after another. Usually gotten
+ via the :class:`~.DataSet` classes `phase_vec` property.
+ fwd_model : :class:`~.ForwardModel`
+ The Forward model instance which should be used for the simulation of the phase maps which
+ will be compared to `y`.
+ lam : float, optional
+ Regularization parameter used in the Hessian matrix multiplication. Default is 0.
+
+ '''
+
+ LOG = logging.getLogger(__name__+'.Costfunction')
+
+ def __init__(self, y, fwd_model, lam=0):
+ self.LOG.debug('Calling __init__')
+ self.y = y
+ self.fwd_model = fwd_model
+ self.lam = lam
+ self.Se_inv = eye(len(y))
+ self.LOG.debug('Created '+str(self))
def __call__(self, x):
- '''TEST DOCSTRING FOR CALL'''
- # TODO: Docstring!
+ self.LOG.debug('Calling __call__')
y = self.y
F = self.F
Se_inv = self.Se_inv
-# print 'Costfunction - __call__ - input: ', len(x)
- result = (F(x)-y).dot(Se_inv.dot(F(x)-y))
-# print 'Costfunction - __call__ - output:', result
- return result
+ return (F(x)-y).dot(Se_inv.dot(F(x)-y))
+
+ def __repr__(self):
+ self.LOG.debug('Calling __repr__')
+ return '%s(fwd_model=%r, lam=%r)' % (self.__class__, self.fwd_model, self.lam)
+
+ def __str__(self):
+ self.LOG.debug('Calling __str__')
+ return 'Costfunction(fwd_model=%s, lam=%s)' % (self.__class__, self.fwd_model, self.lam)
def jac(self, x):
- # TODO: Docstring!
+ '''Calculate the derivative of the costfunction for a given magnetization distribution.
+
+ Parameters
+ ----------
+ x : :class:`~numpy.ndarray` (N=1)
+ Vectorized magnetization distribution, for which the Jacobi vector is calculated.
+
+ Returns
+ -------
+ result : :class:`~numpy.ndarray` (N=1)
+ Jacobi vector which represents the cost derivative of all voxels of the magnetization.
+
+ '''
+ self.LOG.debug('Calling jac')
y = self.y
- F = self.F
+ F = self.fwd_model
Se_inv = self.Se_inv
-# print 'Costfunction - jac - input: ', len(x)
- result = F.jac_T_dot(x, Se_inv.dot(F(x)-y))
-# print 'Costfunction - jac - output:', len(result)
- return result
+ return F.jac_T_dot(x, Se_inv.dot(F(x)-y))
def hess_dot(self, x, vector):
- # TODO: Docstring!
- F = self.F
+ '''Calculate the product of a `vector` with the Hession matrix of the costfunction.
+
+ Parameters
+ ----------
+ x : :class:`~numpy.ndarray` (N=1)
+ Vectorized magnetization distribution at which the Hessian is calculated. The Hessian
+ is constant in this case, thus `x` can be set to None (it is not used int the
+ computation). It is implemented for the case that in the future nonlinear problems
+ have to be solved.
+ vector : :class:`~numpy.ndarray` (N=1)
+ Vectorized magnetization distribution which is multiplied by the Hessian.
+
+ Returns
+ -------
+ result : :class:`~numpy.ndarray` (N=1)
+ Product of the input `vector` with the Hessian matrix of the costfunction.
+
+ '''
+ self.LOG.debug('Calling hess_dot')
+ F = self.fwd_model
Se_inv = self.Se_inv
-# print 'Costfunction - hess_dot - input: ', len(vector)
- result = F.jac_T_dot(x, Se_inv.dot(F.jac_dot(x, vector)))
-# print 'Costfunction - hess_dot - output:', len(result)
- return result
- # TODO: Murks!
+ lam = self.lam
+ return F.jac_T_dot(x, Se_inv.dot(F.jac_dot(x, vector))) + lam*vector
-class CFAdapter:
- # TODO: Useless at the moment, because the interface is exactly the same. Use as template!
-
- def __init__(self, costfunction):
- # TODO: Docstring!
- self.costfunction = costfunction
+class CFAdapterScipyCG(LinearOperator):
- def __call__(self, x):
- # TODO: Docstring!
- return self.costfunction(x)
+ '''Adapter class making the :class:`~.Costfunction` class accessible for scipy cg methods.
- def jac(self, x):
- # TODO: Docstring!
- return self.costfunction.jac(x)
+ This class provides an adapter for the :class:`~.Costfunction` to be usable with the
+ :func:`~.scipy.sparse.linalg.cg` function. the :func:`~.matvec` function is overwritten to
+ implement a multiplication with the Hessian of the adapted costfunction. This is used in the
+ :func:`~pyramid.reconstruction.optimice_sparse_cg` function of the
+ :mod:`~pyramid.reconstruction` module.
- def hess_dot(self, x, vector):
- # TODO: Docstring!
- return self.costfunction.hess_dot(x, vector)
+ Attributes
+ ----------
+ cost : :class:`~.Costfunction`
+ :class:`~.Costfunction` class which should be made usable in the
+ :func:`~.scipy.sparse.linalg.cg` function.
+
+ '''
+
+ LOG = logging.getLogger(__name__+'.CFAdapterScipyCG')
+
+ def __init__(self, cost):
+ self.LOG.debug('Calling __init__')
+ self.cost = cost
+
+ def matvec(self, vector):
+ '''Matrix-vector multiplication with the Hessian of the adapted costfunction.
+
+ Parameters
+ ----------
+ vector : :class:`~numpy.ndarray` (N=1)
+ Vector which will be multiplied by the Hessian matrix provided by the adapted
+ costfunction.
+
+ '''
+ self.LOG.debug('Calling matvec')
+ return self.cost.hess_dot(None, vector)
+
+ @property
+ def shape(self):
+ return (3*self.cost.fwd_model.size_3d, 3*self.cost.fwd_model.size_3d)
+
+ @property
+ def dtype(self):
+ return np.dtype("d")
diff --git a/pyramid/datacollection.py b/pyramid/datacollection.py
deleted file mode 100644
index 496173e8a5c424e5a5b72e4eab16394fbe7e12b1..0000000000000000000000000000000000000000
--- a/pyramid/datacollection.py
+++ /dev/null
@@ -1,102 +0,0 @@
-# -*- coding: utf-8 -*-
-"""This module provides the :class:`~.DataCollection` class for the collection of phase maps
-and additional data like corresponding projectors."""
-
-
-import logging
-
-import numpy as np
-from numbers import Number
-
-from pyramid.phasemap import PhaseMap
-from pyramid.projector import Projector
-
-
-class DataCollection(object):
-
- '''Class for collecting phase maps and corresponding projectors.
-
- Represents a collection of (e.g. experimentally derived) phase maps, stored as
- :class:`~.PhaseMap` objects and corresponding projectors stored as :class:`~.Projector`
- objects. At creation, the grid spacing `a` and the dimension `dim_uv` of the projected grid.
- Data can be added via the :func:`~.append` method, where a :class:`~.PhaseMap` and a
- :class:`~.Projector` have to be given as tuple argument.
-
- Attributes
- ----------
- a: float
- The grid spacing in nm.
- dim_uv: tuple (N=2)
- Dimensions of the projected grid.
- phase_maps:
- A list of all stored :class:`~.PhaseMap` objects.
- projectors:
- A list of all stored :class:`~.Projector` objects.
- phase_vec: :class:`~numpy.ndarray` (N=1)
- The concatenaded, vectorized phase of all ;class:`~.PhaseMap` objects.
-
- '''
-
- @property
- def a(self):
- return self._a
-
- @property
- def dim_uv(self):
- return self._dim_uv
-
- @property
- def phase_vec(self):
- return np.concatenate([p.phase_vec for p in self.phase_maps])
-
- @property
- def phase_maps(self):
- return [d[0] for d in self.data]
-
- @property
- def projectors(self):
- return [d[1] for d in self.data]
-
- def __init__(self, a, dim_uv, b_0):
- self.log = logging.getLogger(__name__)
- assert isinstance(a, Number), 'Grid spacing has to be a number!'
- assert a >= 0, 'Grid spacing has to be a positive number!'
- self._a = a
- assert isinstance(dim_uv, tuple) and len(dim_uv)==2, \
- 'Dimension has to be a tuple of length 2!'
- assert np.all([])
- self._dim_uv = dim_uv
- self.b_0 = b_0
- self.data = []
- # TODO: make it work
-# self.log.info('Created:', str(self))
-
- def __repr__(self):
- self.log.info('Calling __repr__')
- return '%s(a=%r, dim_uv=%r)' % (self.__class__, self.a, self.dim_uv)
-
- def __str__(self):
- self.log.info('Calling __str__')
- return 'DataCollection(%s, a=%s, dim_uv=%s, data_count=%s)' % \
- (self.__class__, self.a, self.dim_uv, len(self.data))
-
- def append(self, (phase_map, projector)):
- '''Appends a data pair of phase map and projection infos to the data collection.`
-
- Parameters
- ----------
- (phase_map, projector): tuple (N=2)
- tuple which contains a :class:`~.PhaseMap` object and a :class:`~.Projector` object,
- which should be added to the data collection.
-
- Returns
- -------
- None
-
- '''
- self.log.info('Calling append')
- assert isinstance(phase_map, PhaseMap) and isinstance(projector, Projector), \
- 'Argument has to be a tuple of a PhaseMap and a Projector object!'
- assert phase_map.dim_uv == self.dim_uv, 'Added phasemap must have the same dimension!'
- assert projector.dim_uv == self.dim_uv, 'Projector dimensions must match!'
- self.data.append((phase_map, projector))
diff --git a/pyramid/dataset.py b/pyramid/dataset.py
new file mode 100644
index 0000000000000000000000000000000000000000..4cecfe8494fcc9daee6ee5f6fdfd66a991b8bf7e
--- /dev/null
+++ b/pyramid/dataset.py
@@ -0,0 +1,159 @@
+# -*- coding: utf-8 -*-
+"""This module provides the :class:`~.DataSet` class for the collection of phase maps
+and additional data like corresponding projectors."""
+
+
+import numpy as np
+from numbers import Number
+
+import matplotlib.pyplot as plt
+
+from pyramid.phasemap import PhaseMap
+from pyramid.phasemapper import PMConvolve
+from pyramid.projector import Projector
+
+import logging
+
+
+class DataSet(object):
+
+ '''Class for collecting phase maps and corresponding projectors.
+
+ Represents a collection of (e.g. experimentally derived) phase maps, stored as
+ :class:`~.PhaseMap` objects and corresponding projectors stored as :class:`~.Projector`
+ objects. At creation, the grid spacing `a` and the dimension `dim_uv` of the projected grid.
+ Data can be added via the :func:`~.append` method, where a :class:`~.PhaseMap` and a
+ :class:`~.Projector` have to be given as tuple argument.
+
+ Attributes
+ ----------
+ a: float
+ The grid spacing in nm.
+ dim_uv: tuple (N=2)
+ Dimensions of the projected grid.
+ phase_maps:
+ A list of all stored :class:`~.PhaseMap` objects.
+ projectors:
+ A list of all stored :class:`~.Projector` objects.
+ phase_vec: :class:`~numpy.ndarray` (N=1)
+ The concatenaded, vectorized phase of all ;class:`~.PhaseMap` objects.
+
+ '''
+
+ LOG = logging.getLogger(__name__+'.DataSet')
+
+ @property
+ def a(self):
+ return self._a
+
+ @property
+ def dim_uv(self):
+ return self._dim_uv
+
+ @property
+ def phase_vec(self):
+ return np.concatenate([p.phase_vec for p in self.phase_maps])
+
+ @property
+ def phase_maps(self):
+ return [d[0] for d in self.data]
+
+ @property
+ def projectors(self):
+ return [d[1] for d in self.data]
+
+ def __init__(self, a, dim_uv, b_0):
+ 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!'
+ self._a = a
+ assert isinstance(dim_uv, tuple) and len(dim_uv)==2, \
+ 'Dimension has to be a tuple of length 2!'
+ self._dim_uv = dim_uv
+ self.b_0 = b_0
+ self.data = []
+ self.LOG.debug('Created:', str(self))
+
+ def __repr__(self):
+ self.LOG.debug('Calling __repr__')
+ return '%s(a=%r, dim_uv=%r, b_0=%r)' % (self.__class__, self.a, self.dim_uv, self.b_0)
+
+ def __str__(self):
+ self.LOG.debug('Calling __str__')
+ return 'DataSet(a=%s, dim_uv=%s, b_0=%s)' % (self.a, self.dim_uv, self.b_0)
+
+ def append(self, (phase_map, projector)):
+ '''Appends a data pair of phase map and projection infos to the data collection.`
+
+ Parameters
+ ----------
+ (phase_map, projector): tuple (N=2)
+ tuple which contains a :class:`~.PhaseMap` object and a :class:`~.Projector` object,
+ which should be added to the data collection.
+
+ Returns
+ -------
+ None
+
+ '''
+ 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 phase_map.dim_uv == self.dim_uv, 'Added phasemap must have the same dimension!'
+ assert projector.dim_uv == self.dim_uv, 'Projector dimensions must match!'
+ self.data.append((phase_map, projector))
+
+ def display(self, phase_maps=None, title='Phase Map', cmap='RdBu', limit=None, norm=None):
+ '''Display all phasemaps saved in the :class:`~.DataSet` as a colormesh.
+
+ Parameters
+ ----------
+ phase_maps : list of :class:`~.PhaseMap`, optional
+ List of phase_maps to display with annotations from the projectors. If none are
+ given, the phase_maps in the dataset are used (which is the default behaviour).
+ title : string, optional
+ The main part of the title of the plots. The default is 'Phase Map'. Additional
+ projector info is appended to this.
+ cmap : string, optional
+ The :class:`~matplotlib.colors.Colormap` which is used for the plots as a string.
+ The default is 'RdBu'.
+ limit : float, optional
+ Plotlimit for the phase in both negative and positive direction (symmetric around 0).
+ If not specified, the maximum amplitude of the phase is used.
+ norm : :class:`~matplotlib.colors.Normalize` or subclass, optional
+ Norm, which is used to determine the colors to encode the phase information.
+ If not specified, :class:`~matplotlib.colors.Normalize` is automatically used.
+
+ Returns
+ -------
+ None
+
+ '''
+ self.LOG.debug('Calling display')
+ if phase_maps is None:
+ phase_maps = self.phase_maps
+ [phase_map.display_phase('{} ({})'.format(title, self.projectors[i].get_info()),
+ cmap, limit, norm)
+ for (i, phase_map) in enumerate(self.phase_maps)]
+ plt.show()
+
+ def create_phase_maps(self, mag_data):
+ '''Create a list of phasemaps with the projectors in the dataset for a given
+ :class:`~.MagData` object.
+
+ Parameters
+ ----------
+ mag_data : :class:`~.MagData`
+ Magnetic distribution to which the projectors of the dataset should be applied.
+
+ Returns
+ -------
+ phase_maps : list of :class:`~.phasemap.PhaseMap`
+ A list of the phase_maps resulting from the projections specified in the dataset.
+
+ Notes
+ -----
+ For the phasemapping, the :class:`~.PMConvolve` class is used.
+
+ '''
+ return [PMConvolve(self.a, proj, self.b_0)(mag_data) for proj in self.projectors]
diff --git a/pyramid/forwardmodel.py b/pyramid/forwardmodel.py
index 44694ccc0b43e978c98c850426ee82155f92ad2e..0dae0591ac17d781132b8d82476df81693b06b2f 100644
--- a/pyramid/forwardmodel.py
+++ b/pyramid/forwardmodel.py
@@ -1,9 +1,6 @@
# -*- coding: utf-8 -*-
-"""
-Created on Mon Jan 06 14:02:18 2014
-
-@author: Jan
-"""
+"""This module provides the :class:`~.ForwardModel` class which represents a strategy to map a
+threedimensional magnetization distribution onto a two-dimensional phase map."""
import numpy as np
@@ -11,8 +8,40 @@ import numpy as np
from pyramid.kernel import Kernel
from pyramid.projector import Projector
+import logging
+
+
class ForwardModel:
- # TODO: Docstring!
+
+ '''Class for mapping 3D magnetic distributions to 2D phase maps.
+
+ Represents a strategy for the mapping of a 3D magnetic distribution to two-dimensional
+ phase maps. Can handle a list of `projectors` of :class:`~.Projector` objects, which describe
+ different projection angles, so many phase_maps can be created from one magnetic distribution.
+
+ Attributes
+ ----------
+ projectors : list of :class:`~.Projector`
+ A list of all :class:`~.Projector` objects representing the projection directions.
+ kernel : :class:`~.Kernel`
+ A kernel which describes the phasemapping of the 2D projected magnetization distribution.
+ a : float
+ The grid spacing in nm. Will be extracted from the `kernel`.
+ dim : tuple (N=3)
+ Dimensions of the 3D magnetic distribution. Is extracted from the first projector of the
+ `projectors` list.
+ dim_uv: tuple (N=2)
+ Dimensions of the projected grid. Is extracted from the `kernel`.
+ size_3d : int
+ Number of voxels of the 3-dimensional grid. Is extracted from the first projector of the
+ `projectors` list. Is extracted from the first projector of the `projectors` list.
+ size_2d : int
+ Number of pixels of the 2-dimensional projected grid. Is extracted from the first projector
+ of the `projectors` list. Is extracted from the first projector of the `projectors` list.
+
+ '''
+
+ LOG = logging.getLogger(__name__+'.ForwardModel')
@property
def projectors(self):
@@ -34,40 +63,80 @@ class ForwardModel:
self._kernel = kernel
def __init__(self, projectors, kernel):
- # TODO: Docstring!
+ self.LOG.debug('Calling __init__')
self.kernel = kernel
self.a = kernel.a
- self.dim_uv = kernel.dim_uv
self.projectors = projectors
+ self.dim = self.projectors[0].dim
+ self.dim_uv = kernel.dim_uv
+ self.size_2d = self.projectors[0].size_2d
+ self.size_3d = self.projectors[0].size_3d
+ self.LOG.debug('Creating '+str(self))
def __call__(self, x):
- # TODO: Docstring!
-# print 'FWD Model - __call__ - input: ', len(x)
+ self.LOG.debug('Calling __call__')
+ result = [self.kernel(projector(x)) for projector in self.projectors]
+ return np.reshape(result, -1)
+
+ def jac_dot(self, x, vector):
+ '''Calculate the product of the Jacobi matrix with a given `vector`.
+
+ Parameters
+ ----------
+ x : :class:`~numpy.ndarray` (N=1)
+ Evaluation point of the jacobi-matrix. The Jacobi matrix is constant for a linear
+ problem, thus `x` can be set to None (it is not used int the computation). It is
+ implemented for the case that in the future nonlinear problems have to be solved.
+ vector : :class:`~numpy.ndarray` (N=1)
+ Vectorized form of the 3D magnetization distribution. First the `x`, then the `y` and
+ lastly the `z` components are listed.
+
+ Returns
+ -------
+ result_vector : :class:`~numpy.ndarray` (N=1)
+ Product of the Jacobi matrix (which is not explicitely calculated) with the input
+ `vector`.
+
+ '''
+ self.LOG.debug('Calling jac_dot')
result = [self.kernel.jac_dot(projector.jac_dot(x)) for projector in self.projectors]
result = np.reshape(result, -1)
-# print 'FWD Model - __call__ - output:', len(result)
return result
-
- # TODO: jac_dot ausschreiben!!!
-
- def jac_dot(self, x, vector):
- # TODO: Docstring! multiplication with the jacobi-matrix (may depend on x)
-# print 'FWD Model - jac_dot - input: ', len(vector)
result = self(vector)
-# print 'FWD Model - jac_dot - output:', len(result)
- return result # The jacobi-matrix does not depend on x in a linear problem
-
+ return result
+
def jac_T_dot(self, x, vector):
- # TODO: Docstring! multiplication with the jacobi-matrix (may depend on x)
- # The jacobi-matrix does not depend on x in a linear problem
-# print 'FWD Model - jac_T_dot - input: ', len(vector)
+ ''''Calculate the product of the transposed Jacobi matrix with a given `vector`.
+
+ Parameters
+ ----------
+ x : :class:`~numpy.ndarray` (N=1)
+ Evaluation point of the jacobi-matrix. The jacobi matrix is constant for a linear
+ problem, thus `x` can be set to None (it is not used int the computation). Is used
+ for the case that in the future nonlinear problems have to be solved.
+ vector : :class:`~numpy.ndarray` (N=1)
+ Vectorized form of all 2D phase maps one after another in one vector.
+
+ Returns
+ -------
+ result_vector : :class:`~numpy.ndarray` (N=1)
+ Product of the transposed Jacobi matrix (which is not explicitely calculated) with
+ the input `vector`.
+
+ '''
+ self.LOG.debug('Calling jac_T_dot')
size_3d = self.projectors[0].size_3d
size_2d = np.prod(self.dim_uv)
result = np.zeros(3*size_3d)
for (i, projector) in enumerate(self.projectors):
result += projector.jac_T_dot(self.kernel.jac_T_dot(vector[i*size_2d:(i+1)*size_2d]))
-# result = [projector.jac_T_dot(self.kernel.jac_T_dot(vector[i*size_2d:(i+1)*size_2d]))
-# for (i, projector) in enumerate(self.projectors)]
- result = np.reshape(result, -1)
-# print 'FWD Model - jac_T_dot - output:', len(result)
- return result
+ return np.reshape(result, -1)
+
+ def __repr__(self):
+ self.LOG.debug('Calling __repr__')
+ return '%s(projectors=%r, kernel=%r)' % (self.__class__, self.projectors, self.kernel)
+
+ def __str__(self):
+ self.LOG.debug('Calling __str__')
+ return 'ForwardModel(%s -> %s, %s projections, kernel=%s)' % \
+ (self.dim, self.dim_uv, len(self.projectors), self.kernel)
diff --git a/pyramid/kernel.py b/pyramid/kernel.py
index 00752a9f5cf98fe212f50204d0d325e89febe6ef..8c0859c1266ed4251eab78079ad644e956f9b308 100644
--- a/pyramid/kernel.py
+++ b/pyramid/kernel.py
@@ -3,58 +3,40 @@
single magnetized pixel."""
-import logging
-
import numpy as np
import pyramid.numcore.kernel_core as core
+import logging
+
PHI_0 = -2067.83 # magnetic flux in T*nm²
-# TODO: sign?
class Kernel(object):
'''Class for calculating kernel matrices for the phase calculation.
-
-
-
- This module provides the :class:`~.Kernel` class whose instances can be used to calculate and
- store the kernel matrix representing the phase of a single pixel for the convolution used in the
- phase calculation. The phasemap of a single pixel for two orthogonal directions (`u` and `v`) are
- stored seperately as 2-dimensional matrices. The Jacobi matrix of the phasemapping just depends
- on the kernel and can be calculated via the :func:`~.get_jacobi` function. Storing the Jacobi
- matrix uses much memory, thus it is also possible to directly get the multiplication of a given
- vector with the (transposed) Jacobi matrix without explicit calculation of the latter.
- It is possible to load data from and save them to NetCDF4 files. See :class:`~.Kernel` for further
- information.
-
Represents the phase of a single magnetized pixel for two orthogonal directions (`u` and `v`),
which can be accessed via the corresponding attributes. The default elementary geometry is
`disc`, but can also be specified as the phase of a `slab` representation of a single
magnetized pixel. During the construction, a few attributes are calculated that are used in
- the convolution during phase calculation.
-
-
- An instance `kernel` of the :class:`~.Kernel` class is callable via:
-
- .. :function:: kernel(vector)
-
- do stuff
-
- :param str sender: do other stuff
- :return: nix
-
- with `vector` being a :class:`~numpy.ndarray` (N=1).
+ the convolution during phase calculation in the different :class:`~Phasemapper` classes.
+ An instance of the :class:`~.Kernel` class can be called as a function with a `vector`,
+ which represents the projected magnetization onto a 2-dimensional grid.
Attributes
----------
- dim_uv : tuple (N=2)
- Dimensions of the projected magnetization grid.
a : float
The grid spacing in nm.
+ dim_uv : tuple of int (N=2)
+ Dimensions of the 2-dimensional projected magnetization grid from which the phase should
+ 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 :module:`~.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)
@@ -65,7 +47,7 @@ class Kernel(object):
The real FFT of the phase contribution of one pixel magnetized in u-direction.
v_fft : :class:`~numpy.ndarray` (N=3)
The real FFT of the phase contribution of one pixel magnetized in v-direction.
- dim_fft : tuple (N=2)
+ dim_fft : tuple of int (N=2)
Dimensions of the grid, which is used for the FFT. Calculated by adding the dimensions
`dim_uv` of the magnetization grid and the dimensions of the kernel (given by
``2*dim_uv-1``)
@@ -74,23 +56,12 @@ 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.
- '''# TODO: Can be used for several PhaseMappers via the fft arguments or via calling!
-
- def __init__(self, a, dim_uv, b_0=1., numcore=True, geometry='disc'):
- '''Constructor for a :class:`~.Kernel` object for representing a kernel matrix.
+ '''
- Parameters
- ----------
- a : float
- The grid spacing in nm.
- dim_uv : tuple (N=2)
- Dimensions of the projected magnetization grid.
- geometry : {'disc', 'slab'}, optional
- The elementary geometry of the single magnetized pixel.
-
- ''' # TODO: Docstring
- self.log = logging.getLogger(__name__)
- self.log.info('Calling __init__')
+ LOG = logging.getLogger(__name__+'.Kernel')
+
+ def __init__(self, a, dim_uv, b_0=1., numcore=True, geometry='disc'):
+ self.LOG.debug('Calling __init__')
# Function for the phase of an elementary geometry:
def get_elementary_phase(geometry, n, m, a):
if geometry == 'disc':
@@ -104,7 +75,7 @@ class Kernel(object):
return 0.5 * (F_a(n-0.5, m-0.5) - F_a(n+0.5, m-0.5)
- F_a(n-0.5, m+0.5) + F_a(n+0.5, m+0.5))
# Set basic properties:
- self.dim_uv = dim_uv # !!! size of the FOV, not the kernel (kernel is bigger)!
+ self.dim_uv = dim_uv # Size of the FOV, not the kernel (kernel is bigger)!
self.a = a
self.numcore = numcore
self.geometry = geometry
@@ -122,39 +93,30 @@ class Kernel(object):
self.slice_fft = (slice(dim_uv[0]-1, 2*dim_uv[0]-1), slice(dim_uv[1]-1, 2*dim_uv[1]-1))
self.u_fft = np.fft.rfftn(self.u, self.dim_fft)
self.v_fft = np.fft.rfftn(self.v, self.dim_fft)
- self.log.info('Created '+str(self))
-
- def __call__(self, x):
- '''Test'''
- self.log.info('Calling __call__')
-# print 'Kernel - __call__:', len(x)
- if self.numcore:
- return self._multiply_jacobi_core(x)
+ if numcore:
+ self.LOG.debug('numcore activated')
+ self._multiply_jacobi_method = self._multiply_jacobi_core
else:
- return self._multiply_jacobi(x)
- # TODO: Bei __init__ variable auf die entsprechende Funktion setzen.
+ self._multiply_jacobi_method = self._multiply_jacobi
+ self.LOG.debug('Created '+str(self))
+
def __repr__(self):
- self.log.info('Calling __repr__')
+ 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)
- def jac_dot(self, vector):
- '''TEST'''# TODO: Docstring
- self.log.info('Calling jac_dot')
-# print 'Kernel - jac_dot:', len(vector)
- if self.numcore:
- return self._multiply_jacobi_core(vector)
- else:
- return self._multiply_jacobi(vector)
+ 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)
- def jac_T_dot(self, vector):
- # TODO: Docstring
- self.log.info('Calling jac_dot_T')
-# print 'Kernel - jac_T_dot:', len(vector)
- return self._multiply_jacobi_T(vector)
+ def __call__(self, x):
+ '''Test'''
+ self.LOG.debug('Calling __call__')
+ return self._multiply_jacobi_method(x)
- def _multiply_jacobi(self, vector):
+ def jac_dot(self, vector):
'''Calculate the product of the Jacobi matrix with a given `vector`.
Parameters
@@ -163,14 +125,38 @@ class Kernel(object):
Vectorized form of the magnetization in `u`- and `v`-direction of every pixel
(row-wise). The first ``N**2`` elements have to correspond to the `u`-, the next
``N**2`` elements to the `v`-component of the magnetization.
-
+
Returns
-------
result : :class:`~numpy.ndarray` (N=1)
Product of the Jacobi matrix (which is not explicitely calculated) with the vector.
- '''# TODO: move!
- self.log.info('Calling _multiply_jacobi')
+ '''
+ self.LOG.debug('Calling jac_dot')
+ return self._multiply_jacobi_method(vector)
+
+ def jac_T_dot(self, vector):
+ '''Calculate the product of the transposed Jacobi matrix with a given `vector`.
+
+ Parameters
+ ----------
+ vector : :class:`~numpy.ndarray` (N=1)
+ Vectorized form of the magnetization in `u`- and `v`-direction of every pixel
+ (row-wise). The first ``N**2`` elements have to correspond to the `u`-, the next
+ ``N**2`` elements to the `v`-component of the magnetization.
+
+ Returns
+ -------
+ result : :class:`~numpy.ndarray` (N=1)
+ Product of the transposed Jacobi matrix (which is not explicitely calculated) with
+ the vector.
+
+ '''
+ self.LOG.debug('Calling jac_dot_T')
+ return self._multiply_jacobi_T(vector)
+
+ def _multiply_jacobi(self, vector):
+ self.LOG.debug('Calling _multiply_jacobi')
v_dim, u_dim = self.dim_uv
size = np.prod(self.dim_uv)
assert len(vector) == 2*size, 'vector size not compatible!'
@@ -187,26 +173,11 @@ class Kernel(object):
return result
def _multiply_jacobi_T(self, vector):
- '''Calculate the product of the transposed Jacobi matrix with a given `vector`.
-
- Parameters
- ----------
- vector : :class:`~numpy.ndarray` (N=1)
- Vectorized form of the magnetization in `u`- and `v`-direction of every pixel
- (row-wise). The first ``N**2`` elements have to correspond to the `u`-, the next
- ``N**2`` elements to the `v`-component of the magnetization.
-
- Returns
- -------
- result : :class:`~numpy.ndarray` (N=1)
- Product of the transposed Jacobi matrix (which is not explicitely calculated) with
- the vector.
-
- '''# TODO: move!
- self.log.info('Calling _multiply_jacobi_T')
+ self.LOG.debug('Calling _multiply_jacobi_T')
v_dim, u_dim = self.dim_uv
size = np.prod(self.dim_uv)
- assert len(vector) == size, 'vector size not compatible! vector: {}, size: {}'.format(len(vector),size)
+ assert len(vector) == size, \
+ 'vector size not compatible! vector: {}, size: {}'.format(len(vector),size)
result = np.zeros(2*size)
for s in range(size): # row-wise (two rows at a time, u- and v-component)
i = s % u_dim
@@ -220,7 +191,7 @@ class Kernel(object):
return result
def _multiply_jacobi_core(self, vector):
- self.log.info('Calling _multiply_jacobi_core')
+ self.LOG.debug('Calling _multiply_jacobi_core')
result = np.zeros(np.prod(self.dim_uv))
core.multiply_jacobi_core(self.dim_uv[0], self.dim_uv[1], self.u, self.v, vector, result)
return result
diff --git a/pyramid/logfile.log b/pyramid/logfile.log
deleted file mode 100644
index e69de29bb2d1d6434b8b29ae775ad8c2e48c5391..0000000000000000000000000000000000000000
diff --git a/pyramid/logging.ini b/pyramid/logging.ini
index 266c37587b5288c8d796dfc92c4025f405264fba..d206f57ccc7f98b5126ac8c0a0de078be16a8a49 100644
--- a/pyramid/logging.ini
+++ b/pyramid/logging.ini
@@ -23,7 +23,7 @@ args=tuple()
[handler_file]
class=logging.FileHandler
-level=INFO
+level=DEBUG
formatter=file
args=('logfile.log', 'w')
diff --git a/pyramid/magcreator.py b/pyramid/magcreator.py
index 82623a1888ecce67d9b1b551502e9777f1187311..5608b72eab95a047a3104fbc2648c47d8ffb186a 100644
--- a/pyramid/magcreator.py
+++ b/pyramid/magcreator.py
@@ -13,14 +13,20 @@ center.
"""
-
import numpy as np
from numpy import pi
+import abc
+
+import logging
+
+
+LOG = logging.getLogger(__name__)
+
class Shapes(object):
- '''Class containing functions for generating magnetic shapes.
+ '''Abstract class containing functions for generating magnetic shapes.
The :class:`~.Shapes` class is a collection of some methods that return a 3-dimensional
matrix that represents the magnetized volume and consists of values between 0 and 1.
@@ -28,7 +34,9 @@ class Shapes(object):
:class:`~pyramid.magdata.MagData` objects which store the magnetic informations.
'''
- #TODO: Abstract class (test if this works!)
+
+ __metaclass__ = abc.ABCMeta
+ LOG = logging.getLogger(__name__+'.Shapes')
@classmethod
def slab(cls, dim, center, width):
@@ -49,6 +57,7 @@ class Shapes(object):
The magnetic shape as a 3D-array with values between 1 and 0.
'''
+ 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!'
@@ -83,11 +92,12 @@ class Shapes(object):
The magnetic shape as a 3D-array with values between 1 and 0.
'''
+ 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!'
assert height > 0 and np.shape(height) == (), 'Height has to be a positive scalar value!'
- assert axis == 'z' or axis == 'y' or axis == 'x', 'Axis has to be x, y or z (as a string)!'
+ assert axis in {'z', 'y', 'x'}, 'Axis has to be x, y or z (as a string)!'
if axis == 'z':
mag_shape = np.array([[[np.hypot(x - center[2], y - center[1]) <= radius
and abs(z - center[0]) <= height / 2
@@ -108,6 +118,58 @@ class Shapes(object):
for z in range(dim[0])])
return mag_shape
+ @classmethod
+ def ellipse(cls, dim, center, width, height, axis='z'):
+ '''Create the shape of an elliptical cylinder in x-, y-, or z-direction.
+
+ Parameters
+ ----------
+ dim : tuple (N=3)
+ The dimensions of the grid `(z, y, x)`.
+ center : tuple (N=3)
+ The center of the ellipse in pixel coordinates `(z, y, x)`.
+ width : tuple (N=2)
+ Length of the two axes of the ellipse.
+ height : float
+ The height of the ellipse in pixel coordinates.
+ axis : {'z', 'y', 'x'}, optional
+ The orientation of the ellipse axis. The default is 'z'.
+
+ Returns
+ -------
+ mag_shape : :class:`~numpy.ndarray` (N=3)
+ The magnetic shape as a 3D-array with values between 1 and 0.
+
+ '''
+ 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!'
+ assert height > 0 and np.shape(height) == (), 'Height has to be a positive scalar value!'
+ assert axis in {'z', 'y', 'x'}, 'Axis has to be x, y or z (as a string)!'
+ if axis == 'z':
+ mag_shape = np.array([[[np.hypot((x-center[2])/(width[1]/2.),
+ (y-center[1])/(width[0]/2.))<=1
+ and abs(z - center[0]) <= height / 2
+ for x in range(dim[2])]
+ for y in range(dim[1])]
+ for z in range(dim[0])])
+ elif axis == 'y':
+ mag_shape = np.array([[[np.hypot((x-center[2])/(width[1]/2.),
+ (z-center[0])/(width[0]/2.))<=1
+ and abs(y-center[1]) <= height / 2
+ for x in range(dim[2])]
+ for y in range(dim[1])]
+ for z in range(dim[0])])
+ elif axis == 'x':
+ mag_shape = np.array([[[np.hypot((y-center[1])/(width[1]/2.),
+ (z-center[0])/(width[0]/2.))<=1
+ and abs(z - center[0]) <= height / 2
+ for x in range(dim[2])]
+ for y in range(dim[1])]
+ for z in range(dim[0])])
+ return mag_shape
+
@classmethod
def sphere(cls, dim, center, radius):
'''Create the shape of a sphere.
@@ -127,6 +189,7 @@ class Shapes(object):
The magnetic shape as a 3D-array with values between 1 and 0.
'''
+ 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!'
@@ -147,7 +210,7 @@ class Shapes(object):
dim : tuple (N=3)
The dimensions of the grid `(z, y, x)`.
center : tuple (N=3)
- The center of the sphere in pixel coordinates `(z, y, x)`.
+ The center of the ellipsoid in pixel coordinates `(z, y, x)`.
width : tuple (N=3)
The width of the ellipsoid `(z, y, x)`.
@@ -157,6 +220,7 @@ class Shapes(object):
The magnetic shape as a 3D-array with values between 1 and 0.
'''
+ 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!'
@@ -189,9 +253,10 @@ class Shapes(object):
The magnetic shape as a 3D-array with values between 1 and 0.
'''
+ 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 == 'z' or axis == 'y' or axis == 'x', 'Axis has to be x, y or z (as a string)!'
+ assert axis in {'z', 'y', 'x'}, 'Axis has to be x, y or z (as a string)!'
mag_shape = np.zeros(dim)
if axis == 'z':
mag_shape[:, pos[0], pos[1]] = 1
@@ -218,6 +283,7 @@ class Shapes(object):
The magnetic shape as a 3D-array with values between 1 and 0.
'''
+ 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)
@@ -226,7 +292,7 @@ class Shapes(object):
def create_mag_dist_homog(mag_shape, phi, theta=pi/2, magnitude=1):
- '''Create a 3-dimensional magnetic distribution of a homogeneous magnetized object.
+ '''Create a 3-dimensional magnetic distribution of a homogeneously magnetized object.
Parameters
----------
@@ -247,15 +313,16 @@ 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')
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
y_mag = np.ones(dim) * np.sin(theta) * np.sin(phi) * mag_shape * magnitude
x_mag = np.ones(dim) * np.sin(theta) * np.cos(phi) * mag_shape * magnitude
- return np.array([x_mag, y_mag, z_mag]) # TODO: Better implementation!
+ return np.array([x_mag, y_mag, z_mag])
-def create_mag_dist_vortex(mag_shape, center=None, magnitude=1):
+def create_mag_dist_vortex(mag_shape, center=None, axis='z', magnitude=1):
'''Create a 3-dimensional magnetic distribution of a homogeneous magnetized object.
Parameters
@@ -263,12 +330,14 @@ def create_mag_dist_vortex(mag_shape, center=None, magnitude=1):
mag_shape : :class:`~numpy.ndarray` (N=3)
The magnetic shapes (see :class:`.~Shapes` for examples).
center : tuple (N=2 or N=3), optional
- The vortex center, given in 2D `(x, y)` or 3D `(x, y, z)`, where `z` is discarded.
- Is set to the center of the field of view if not specified.
- Vortex center has to be between two pixels.
+ The vortex center, given in 2D `(v, u)` or 3D `(z, y, x)`, where the perpendicular axis
+ is is discarded. Is set to the center of the field of view if not specified. The vortex
+ center has to be between two pixels.
magnitude : float, optional
- The relative magnitude for the magnetic shape. The default is 1.
-
+ The relative magnitude for the magnetic shape. The default is 1. Negative signs reverse
+ the vortex direction (left-handed instead of right-handed).
+ axis : {'z', 'y', 'x'}, optional
+ The orientation of the vortex axis. The default is 'z'.
Returns
-------
magnitude : tuple (N=3) of :class:`~numpy.ndarray` (N=3)
@@ -276,18 +345,45 @@ def create_mag_dist_vortex(mag_shape, center=None, magnitude=1):
`x`-, `y`- and `z`-direction on the 3-dimensional grid.
'''
+ 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)!'
+ assert center is None or len(center) in {2, 3}, \
+ 'Vortex center has to be defined in 3D or 2D or not at all!'
if center is None:
center = (int(dim[1]/2)-0.5, int(dim[2]/2)-0.5) # center has to be between (!) two pixels
- if len(center) == 3: # if a 3D-center is given, just take the x and y components
- center = (center[1], center[2])
- assert len(center) == 2, 'Vortex center has to be defined in 3D or 2D or not at all!'
- x = np.linspace(-center[1], dim[2]-1-center[1], dim[2])
- y = np.linspace(-center[0], dim[1]-1-center[0], dim[1])
- xx, yy = np.meshgrid(x, y)
- phi = np.arctan2(yy, xx) - pi/2
- z_mag = np.zeros(dim)
- y_mag = -np.ones(dim) * np.sin(phi) * mag_shape * magnitude
- x_mag = -np.ones(dim) * np.cos(phi) * mag_shape * magnitude
- return np.array([x_mag, y_mag, z_mag]) # TODO: Better implementation!
+ if axis == 'z':
+ if len(center) == 3: # if a 3D-center is given, just take the x and y components
+ center = (center[1], center[2])
+ u = np.linspace(-center[1], dim[2]-1-center[1], dim[2])
+ v = np.linspace(-center[0], dim[1]-1-center[0], dim[1])
+ uu, vv = np.meshgrid(u, v)
+ phi = np.expand_dims(np.arctan2(vv, uu)-pi/2, axis=0)
+ phi = np.tile(phi, (dim[0], 1, 1))
+ z_mag = np.zeros(dim)
+ y_mag = -np.ones(dim) * np.sin(phi) * mag_shape * magnitude
+ x_mag = -np.ones(dim) * np.cos(phi) * mag_shape * magnitude
+ elif axis == 'y':
+ if len(center) == 3: # if a 3D-center is given, just take the x and z components
+ center = (center[0], center[2])
+ u = np.linspace(-center[1], dim[2]-1-center[1], dim[2])
+ v = np.linspace(-center[0], dim[0]-1-center[0], dim[0])
+ uu, vv = np.meshgrid(u, v)
+ phi = np.expand_dims(np.arctan2(vv, uu)-pi/2, axis=1)
+ phi = np.tile(phi, (1, dim[1], 1))
+ z_mag = np.ones(dim) * np.sin(phi) * mag_shape * magnitude
+ y_mag = np.zeros(dim)
+ x_mag = np.ones(dim) * np.cos(phi) * mag_shape * magnitude
+ elif axis == 'x':
+ if len(center) == 3: # if a 3D-center is given, just take the z and y components
+ center = (center[0], center[1])
+ u = np.linspace(-center[1], dim[1]-1-center[1], dim[1])
+ v = np.linspace(-center[0], dim[0]-1-center[0], dim[0])
+ uu, vv = np.meshgrid(u, v)
+ phi = np.expand_dims(np.arctan2(vv, uu)-pi/2, axis=2)
+ phi = np.tile(phi, (1, 1, dim[2]))
+ z_mag = -np.ones(dim) * np.sin(phi) * mag_shape * magnitude
+ y_mag = -np.ones(dim) * np.cos(phi) * mag_shape * magnitude
+ x_mag = np.zeros(dim)
+ return np.array([x_mag, y_mag, z_mag])
diff --git a/pyramid/magdata.py b/pyramid/magdata.py
index 0dc7044b92183182a680aeb5c9d6ce32852c8163..d51ed6e894d1a31f27285843287d2c1c82bf0b37 100644
--- a/pyramid/magdata.py
+++ b/pyramid/magdata.py
@@ -2,7 +2,6 @@
"""This module provides the :class:`~.MagData` class for storing of magnetization data."""
-import logging
import numpy as np
from scipy.ndimage.interpolation import zoom
@@ -14,6 +13,8 @@ from numbers import Number
import netCDF4
+import logging
+
class MagData(object):
@@ -21,7 +22,7 @@ class MagData(object):
Represents 3-dimensional magnetic distributions with 3 components which are stored as a
2-dimensional numpy array in `magnitude`, but which can also be accessed as a vector via
- `mag_vec`. :class:`~.MagData` objects support negation, arithmetic operators
+ `mag_vec`. :class:`~.MagData` objects support negation, arithmetic operators
(``+``, ``-``, ``*``) and their augmented counterparts (``+=``, ``-=``, ``*=``), with numbers
and other :class:`~.MagData` objects, if their dimensions and grid spacings match. It is
possible to load data from NetCDF4 or LLG (.txt) files or to save the data in these formats.
@@ -40,9 +41,9 @@ class MagData(object):
Vector containing the magnetic distribution.
'''
-
- log = logging.getLogger()
-
+
+ LOG = logging.getLogger(__name__+'.MagData')
+
@property
def a(self):
return self._a
@@ -51,8 +52,8 @@ class MagData(object):
def a(self, a):
assert isinstance(a, Number), 'Grid spacing has to be a number!'
assert a >= 0, 'Grid spacing has to be a positive number!'
- self._a = a
-
+ self._a = float(a)
+
@property
def dim(self):
return self._dim
@@ -68,7 +69,7 @@ class MagData(object):
assert magnitude.shape[0] == 3, 'First dimension of the magnitude has to be 3!'
self._magnitude = magnitude
self._dim = magnitude.shape[1:]
-
+
@property
def mag_vec(self):
return np.reshape(self.magnitude, -1)
@@ -80,80 +81,70 @@ class MagData(object):
self.magnitude = mag_vec.reshape((3,)+self.dim)
def __init__(self, a, magnitude):
- self.log = logging.getLogger(__name__)
- self.log.info('Calling __init__')
+ self.LOG.debug('Calling __init__')
self.a = a
- self.magnitude = magnitude
- self.log.info('Created '+str(self))
-
- def __getstate__(self):
- d = dict(self.__dict__)
- del d['log']
- return d
-
- def __setstate__(self, d):
- self.__dict__.update(d)
- self.log = logging.getLogger(__name__)
+ self.magnitude = magnitude
+ self.LOG.debug('Created '+str(self))
def __repr__(self):
- self.log.info('Calling __repr__')
+ self.LOG.debug('Calling __repr__')
return '%s(a=%r, magnitude=%r)' % (self.__class__, self.a, self.magnitude)
def __str__(self):
- self.log.info('Calling __str__')
+ self.LOG.debug('Calling __str__')
return 'MagData(a=%s, dim=%s)' % (self.a, self.dim)
def __neg__(self): # -self
- self.log.info('Calling __neg__')
+ self.LOG.debug('Calling __neg__')
return MagData(self.a, -self.magnitude)
-
+
def __add__(self, other): # self + other
- self.log.info('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.info('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.info('Adding an offset')
+ self.LOG.debug('Adding an offset')
return MagData(self.a, self.magnitude+other)
def __sub__(self, other): # self - other
- self.log.info('Calling __sub__')
+ self.LOG.debug('Calling __sub__')
return self.__add__(-other)
def __mul__(self, other): # self * other
- self.log.info('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.info('Calling __radd__')
+ self.LOG.debug('Calling __radd__')
return self.__add__(other)
def __rsub__(self, other): # other - self
- self.log.info('Calling __rsub__')
+ self.LOG.debug('Calling __rsub__')
return -self.__sub__(other)
def __rmul__(self, other): # other * self
- self.log.info('Calling __rmul__')
+ self.LOG.debug('Calling __rmul__')
return self.__mul__(other)
-
+
def __iadd__(self, other): # self += other
- self.log.info('Calling __iadd__')
+ self.LOG.debug('Calling __iadd__')
return self.__add__(other)
def __isub__(self, other): # self -= other
- self.log.info('Calling __isub__')
+ self.LOG.debug('Calling __isub__')
return self.__sub__(other)
def __imul__(self, other): # self *= other
- self.log.info('Calling __imul__')
+ self.LOG.debug('Calling __imul__')
return self.__mul__(other)
-
+
def copy(self):
'''Returns a copy of the :class:`~.MagData` object
@@ -167,7 +158,7 @@ class MagData(object):
A copy of the :class:`~.MagData`.
'''
- self.log.info('Calling copy7')
+ self.LOG.debug('Calling copy')
return MagData(self.a, self.magnitude.copy())
def scale_down(self, n=1):
@@ -188,7 +179,7 @@ class MagData(object):
Only possible, if each axis length is a power of 2!
'''
- self.log.info('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
@@ -196,11 +187,29 @@ class MagData(object):
# 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))
-
+
def scale_up(self, n=1, order=0):
- # TODO: Docstring!
- self.log.info('Calling scale_up')
+ '''Scale up the magnetic distribution using spline interpolation of the requested order.
+
+ Parameters
+ ----------
+ n : int, optional
+ Power of 2 with which the grid is scaled. Default is 1, which means every axis is
+ increased by a factor of ``2**1 = 2``.
+ order : int, optional
+ The order of the spline interpolation, which has to be in the range between 0 and 5
+ and defaults to 0.
+
+ Returns
+ -------
+ None
+
+ Notes
+ -----
+ Acts in place and changes dimensions and grid spacing accordingly.
+ '''
+ 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!'
@@ -210,7 +219,28 @@ class MagData(object):
zoom(self.magnitude[2], zoom=2**n, order=order)))
def pad(self, x_pad, y_pad, z_pad):
- # TODO: Docstring!
+ '''Pad the current magnetic distribution with zeros for each individual axis.
+
+ Parameters
+ ----------
+ x_pad : int
+ Number of zeros which should be padded on both sides of the x-axis.
+ y_pad : int
+ Number of zeros which should be padded on both sides of the y-axis.
+ z_pad : int
+ Number of zeros which should be padded on both sides of the z-axis.
+
+ Returns
+ -------
+ None
+
+ Notes
+ -----
+ Acts in place and changes dimensions accordingly.
+ '''
+ assert x_pad >= 0 and isinstance(x_pad, (int, long)), 'x_pad must be a positive integer!'
+ assert y_pad >= 0 and isinstance(y_pad, (int, long)), 'y_pad must be a positive integer!'
+ 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)
@@ -229,7 +259,7 @@ class MagData(object):
None
'''
- self.log.info('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,
@@ -259,7 +289,7 @@ class MagData(object):
A :class:`~.MagData` object containing the loaded data.
'''
- cls.log.info('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])))
@@ -281,7 +311,7 @@ class MagData(object):
None
'''
- self.log.info('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
@@ -307,7 +337,7 @@ class MagData(object):
A :class:`~.MagData` object containing the loaded data.
'''
- cls.log.info('Calling copy')
+ cls.LOG.debug('Calling copy')
mag_file = netCDF4.Dataset(filename, 'r', format='NETCDF4')
a = mag_file.a
magnitude = mag_file.variables['magnitude'][...]
@@ -339,31 +369,31 @@ class MagData(object):
The axis on which the graph is plotted.
'''
- self.log.info('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.info('proj_axis == z')
+ self.LOG.debug('proj_axis == z')
if ax_slice is None:
- self.log.info('ax_slice is None')
+ self.LOG.debug('ax_slice is None')
ax_slice = int(self.dim[0]/2)
mag_slice_u = self.magnitude[0][ax_slice, ...] # x-component
mag_slice_v = 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.info('proj_axis == y')
+ self.LOG.debug('proj_axis == y')
if ax_slice is None:
- self.log.info('ax_slice is None')
+ self.LOG.debug('ax_slice is None')
ax_slice = int(self.dim[1]/2)
mag_slice_u = self.magnitude[0][:, ax_slice, :] # x-component
mag_slice_v = 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.info('proj_axis == x')
+ self.LOG.debug('proj_axis == x')
if ax_slice is None:
- self.log.info('ax_slice is None')
+ self.LOG.debug('ax_slice is None')
ax_slice = int(self.dim[2]/2)
mag_slice_u = self.magnitude[1][..., ax_slice] # y-component
mag_slice_v = self.magnitude[2][..., ax_slice] # z-component
@@ -371,7 +401,7 @@ class MagData(object):
v_label = 'z [px]'
# If no axis is specified, a new figure is created:
if axis is None:
- self.log.info('axis is None')
+ self.LOG.debug('axis is None')
fig = plt.figure(figsize=(8.5, 7))
axis = fig.add_subplot(1, 1, 1, aspect='equal')
axis.quiver(mag_slice_u, mag_slice_v, pivot='middle', angles='xy', scale_units='xy',
@@ -399,7 +429,7 @@ class MagData(object):
None
'''
- self.log.info('Calling quiver_plot3D')
+ self.LOG.debug('Calling quiver_plot3D')
a = self.a
dim = self.dim
# Create points and vector components as lists:
@@ -417,4 +447,5 @@ class MagData(object):
plot = mlab.quiver3d(xx, yy, zz, x_mag, y_mag, z_mag, mode='arrow')
mlab.outline(plot)
mlab.axes(plot)
- mlab.colorbar()
+ mlab.title('TEST', height=0.95, size=0.25)
+ mlab.colorbar(plot, orientation='vertical')
diff --git a/pyramid/numcore/__init__.py b/pyramid/numcore/__init__.py
index fd8ed341b1fb38e8aa62c63d1398301b0745be86..c462822e236e6d8541dec7cb326735de37f715d3 100644
--- a/pyramid/numcore/__init__.py
+++ b/pyramid/numcore/__init__.py
@@ -1,3 +1,4 @@
+# -*- coding: utf-8 -*-
"""Package for numerical core routines which enable fast computations.
Modules
@@ -10,4 +11,4 @@ Notes
-----
Packages are written in `pyx`-format for the use with :mod:`~cython`.
-"""
+"""
\ No newline at end of file
diff --git a/pyramid/numcore/kernel_core.pyx b/pyramid/numcore/kernel_core.pyx
index 7ed2d6ebee216daf5ca8606194dd2293d096c9cc..f0a7412f8ac73cdd78c0e0fb7820915dbf9c01eb 100644
--- a/pyramid/numcore/kernel_core.pyx
+++ b/pyramid/numcore/kernel_core.pyx
@@ -1,11 +1,10 @@
# -*- coding: utf-8 -*-
-"""Numerical core routines for the phase calculation using the real space approach.
+"""Numerical core routines for the :mod:`~.pyramid.kernel` module.
-Provides a helper function to speed up :func:`~pyramid.phasemapper.phase_mag_real` of module
-:mod:`~pyramid.phasemapper`, by using C-speed for the for-loops and by omitting boundary and
-wraparound checks.
+Provides helper functions to speed up kernel calculations in the :class:`~pyramid.kernel.Kernel`
+class, by using C-speed for the for-loops and by omitting boundary and wraparound checks.
-""" # TODO: Docstring!
+"""
import numpy as np
@@ -17,12 +16,12 @@ cimport numpy as np
@cython.boundscheck(False)
@cython.wraparound(False)
-def phase_mag_real_core(
+def multiply_jacobi_core(
unsigned int v_dim, unsigned int u_dim,
- double[:, :] v_phi, double[:, :] u_phi,
- double[:, :] v_mag, double[:, :] u_mag,
- double[:, :] phase, float threshold):
- '''Numerical core routine for the phase calculation using the real space approach.
+ double[:, :] u_phi, double[:, :] v_phi,
+ double[:] vector,
+ double[:] result):
+ '''Numerical core routine for the Jacobi matrix multiplication in the :class:`~.Kernel` class.
Parameters
----------
@@ -30,106 +29,21 @@ def phase_mag_real_core(
Dimensions of the projection along the two major axes.
v_phi, u_phi : :class:`~numpy.ndarray` (N=2)
Lookup tables for the pixel fields oriented in `u`- and `v`-direction.
- v_mag, u_mag : :class:`~numpy.ndarray` (N=2)
- Magnetization components in `u`- and `v`-direction.
- phase : :class:`~numpy.ndarray` (N=2)
- Matrix in which the resulting magnetic phase map should be stored.
- threshold : float
- The `threshold` determines which pixels contribute to the magnetic phase.
+ vector : :class:`~numpy.ndarray` (N=1)
+ Input vector which should be multiplied by the Jacobi-matrix.
+ result : :class:`~numpy.ndarray` (N=1)
+ Vector in which the result of the multiplication should be stored.
Returns
-------
None
- '''
- cdef unsigned int i, j, p, q, p_c, q_c
- cdef double u_m, v_m
- for j in range(v_dim):
- for i in range(u_dim):
- u_m = u_mag[j, i]
- v_m = v_mag[j, i]
- p_c = u_dim - 1 - i
- q_c = v_dim - 1 - j
- if abs(u_m) > threshold:
- for q in range(v_dim):
- for p in range(u_dim):
- phase[q, p] += u_m * u_phi[q_c + q, p_c + p]
- if abs(v_m) > threshold:
- for q in range(v_dim):
- for p in range(u_dim):
- phase[q, p] -= v_m * v_phi[q_c + q, p_c + p]
-
-
-@cython.boundscheck(False)
-@cython.wraparound(False)
-def get_jacobi_core(
- unsigned int v_dim, unsigned int u_dim,
- double[:, :] v_phi, double[:, :] u_phi,
- double[:, :] jacobi):
- '''Numerical core routine for the jacobi matrix calculation.
-
- Parameters
- ----------
- v_dim, u_dim : int
- Dimensions of the projection along the two major axes.
- v_phi, u_phi : :class:`~numpy.ndarray` (N=2)
- Lookup tables for the pixel fields oriented in `u`- and `v`-direction.
- jacobi : :class:`~numpy.ndarray` (N=2)
- Jacobi matrix which is filled by this routine.
-
- Returns
- -------
- None
+ Notes
+ -----
+ The strategy involves iterating over the magnetization first and over the kernel in an inner
+ iteration loop.
'''
- cdef unsigned int i, j, p, q, p_min, p_max, q_min, q_max, u_column, v_column, row
- for j in range(v_dim):
- for i in range(u_dim):
- q_min = (v_dim-1) - j
- q_max = (2*v_dim-1) - j
- p_min = (u_dim-1) - i
- p_max = (2*u_dim-1) - i
- v_column = i + u_dim*j + u_dim*v_dim
- u_column = i + u_dim*j
- for q in range(q_min, q_max):
- for p in range(p_min, p_max):
- q_c = q - (v_dim-1-j) # start at zero for jacobi column
- p_c = p - (u_dim-1-i) # start at zero for jacobi column
- row = p_c + u_dim*q_c
- # u-component:
- jacobi[row, u_column] = u_phi[q, p]
- # v-component (note the minus!):
- jacobi[row, v_column] = -v_phi[q, p]
-
-
-@cython.boundscheck(False)
-@cython.wraparound(False)
-def get_jacobi_core2(
- unsigned int v_dim, unsigned int u_dim,
- double[:, :] v_phi, double[:, :] u_phi,
- double[:, :] jacobi):
- '''DOCSTRING!'''
- # TODO: Docstring!!!
- cdef unsigned int i, j, p, q, p_min, p_max, q_min, q_max, u_column, v_column, row
- for j in range(v_dim):
- for i in range(u_dim):
- # v_dim*u_dim columns for the u-component:
- jacobi[:, i+u_dim*j] = \
- np.reshape(u_phi[v_dim-1-j:(2*v_dim-1)-j, u_dim-1-i:(2*u_dim-1)-i], -1)
- # v_dim*u_dim columns for the v-component (note the minus!):
- jacobi[:, u_dim*v_dim+i+u_dim*j] = \
- -np.reshape(v_phi[v_dim-1-j:(2*v_dim-1)-j, u_dim-1-i:(2*u_dim-1)-i], -1)
-
-
-@cython.boundscheck(False)
-@cython.wraparound(False)
-def multiply_jacobi_core(
- unsigned int v_dim, unsigned int u_dim,
- double[:, :] u_phi, double[:, :] v_phi,
- double[:] vector,
- double[:] result):
- '''DOCSTRING!'''
- # TODO: Docstring!!! Iterate over magnetization and then kernel
cdef unsigned int s, i, j, u_min, u_max, v_min, v_max, ri, u, v, size
size = u_dim * v_dim # Number of pixels
s = 0 # Current pixel (numbered consecutively)
@@ -142,38 +56,7 @@ def multiply_jacobi_core(
# Go over the current kernel cutout [v_min:v_max, u_min:u_max]:
for v in range(v_min, v_min + v_dim):
for u in range(u_min, u_min + u_dim):
- result[ri] += vector[s] * u_phi[v, u]
+ result[ri] += vector[s] * u_phi[v, u]
result[ri] -= vector[s+size] * v_phi[v, u]
ri += 1
s += 1
-
-
-#@cython.boundscheck(False)
-#@cython.wraparound(False)
-def multiply_jacobi_core2(
- int v_dim, int u_dim,
- double[:, :] u_phi, double[:, :] v_phi,
- double[:] vector,
- double[:] result):
- '''DOCSTRING!'''
- # TODO: Docstring!!! Iterate over kernel and then magnetization
- cdef int s1, s2, i, j, u_min, u_max, v_min, v_max, ri, u, v, size, j_min, j_max, i_min, i_max
- size = u_dim * v_dim # Number of pixels
- # Go over complete kernel:
- for v in range(2 * v_dim - 1):
- for u in range(2 * u_dim - 1):
- i_min = max(0, (u_dim - 1) - u)
- i_max = min(u_dim, ((2 * u_dim) - 1) - u)
- j_min = max(0, (v_dim - 1) - v)
- j_max = min(v_dim, ((2 * v_dim) - 1) - v)
- # Iterate over
- for i in range(i_min, i_max):
- s1 = i * u_dim # u-column
- s2 = s1 + size # v-column
- u_min = (u_dim - 1) - i
- v_min = (v_dim - 1) - j_min
- ri = (v - ((v_dim - 1) - j_min)) * u_dim + (u - u_min)
- for j in range(j_min, j_max):
- result[ri] += vector[s1 + j] * u_phi[v, u]
- result[ri] -= vector[s2 + j] * v_phi[v, u]
- ri += u_dim
diff --git a/pyramid/optimizer.py b/pyramid/optimizer.py
deleted file mode 100644
index f45546e505ed0abd0d0ea5694ba6ab95c4abeb18..0000000000000000000000000000000000000000
--- a/pyramid/optimizer.py
+++ /dev/null
@@ -1,143 +0,0 @@
-# -*- coding: utf-8 -*-
-"""Reconstruct magnetic distributions from given phasemaps.
-
-This module reconstructs 3-dimensional magnetic distributions (as :class:`~pyramid.magdata.MagData`
-objects) from a given set of phase maps (represented by :class:`~pyramid.phasemap.PhaseMap`
-objects) by using several model based reconstruction algorithms which use the forward model
-provided by :mod:`~pyramid.projector` and :mod:`~pyramid.phasemapper` and a priori knowledge of
-the distribution.
-So far, only a simple least square algorithm for known pixel locations for 2-dimensional problems
-is implemented (:func:`~.reconstruct_simple_leastsq`), but more complex solutions are planned.
-
-"""# TODO: Docstrings!
-
-
-# TODO: Delete
-''' FRAGEN FÜR JÖRN:
-1) Property und Python "Magic" ist wann sinnvoll?
-2) Aktuelle Designstruktur, welche Designpatterns kommen vor und sind sinnvoll?
-3) Macht logging Sinn?
-4) Ist scipy.optimize.minimize(...) eine gute Wahl?
-5) Wo besser Module, wo Klassen?
-6) Arbitrary grid to cartesian grid
-'''
-
-
-import numpy as np
-
-from scipy.optimize import minimize
-
-from pyramid.kernel import Kernel
-from pyramid.forwardmodel import ForwardModel
-from pyramid.costfunction import Costfunction
-from pyramid.magdata import MagData
-
-
-
-def optimize_cg(data_collection, first_guess):
- # TODO: where should data_collection and first_guess go? here or __init__?
- data = data_collection # TODO: name the parameter 'data' directly...
- mag_0 = first_guess
- x_0 = first_guess.mag_vec
- y = data.phase_vec
- kern = Kernel(data.a, data.dim_uv)
- F = ForwardModel(data.projectors, kern, data.b_0)
- C = Costfunction(y, F)
-
- result = minimize(C, x_0, method='Newton-CG', jac=C.jac, hessp=C.hess_dot,
- options={'maxiter':200, 'disp':True})
-
- x_opt = result.x
-
- mag_opt = MagData(mag_0.a, np.zeros((3,)+mag_0.dim))
-
- mag_opt.mag_vec = x_opt
-
- return mag_opt
-
-
-
-# TODO: Implement the following:
-
-## -*- coding: utf-8 -*-
-#"""Reconstruct magnetic distributions from given phasemaps.
-#
-#This module reconstructs 3-dimensional magnetic distributions (as :class:`~pyramid.magdata.MagData`
-#objects) from a given set of phase maps (represented by :class:`~pyramid.phasemap.PhaseMap`
-#objects) by using several model based reconstruction algorithms which use the forward model
-#provided by :mod:`~pyramid.projector` and :mod:`~pyramid.phasemapper` and a priori knowledge of
-#the distribution.
-#So far, only a simple least square algorithm for known pixel locations for 2-dimensional problems
-#is implemented (:func:`~.reconstruct_simple_leastsq`), but more complex solutions are planned.
-#
-#"""
-#
-#
-#
-#import numpy as np
-#
-#from scipy.optimize import leastsq
-#
-#import pyramid.projector as pj
-#import pyramid.phasemapper as pm
-#from pyramid.magdata import MagData
-#from pyramid.projector import Projection
-#from pyramid.kernel import Kernel
-#
-#
-#def reconstruct_simple_leastsq(phase_map, mask, b_0=1):
-# '''Reconstruct a magnetic distribution for a 2-D problem with known pixel locations.
-#
-# Parameters
-# ----------
-# phase_map : :class:`~pyramid.phasemap.PhaseMap`
-# A :class:`~pyramid.phasemap.PhaseMap` object, representing the phase from which to
-# reconstruct the magnetic distribution.
-# mask : :class:`~numpy.ndarray` (N=3)
-# A boolean matrix (or a matrix consisting of ones and zeros), representing the
-# positions of the magnetized voxels in 3 dimensions.
-# b_0 : float, optional
-# The magnetic induction corresponding to a magnetization `M`\ :sub:`0` in T.
-# The default is 1.
-# Returns
-# -------
-# mag_data : :class:`~pyramid.magdata.MagData`
-# The reconstructed magnetic distribution as a
-# :class:`~pyramid.magdata.MagData` object.
-#
-# Notes
-# -----
-# Only works for a single phase_map, if the positions of the magnetized voxels are known and
-# for slice thickness of 1 (constraint for the `z`-dimension).
-#
-# '''
-# # Read in parameters:
-# y_m = phase_map.phase.reshape(-1) # Measured phase map as a vector
-# a = phase_map.a # Grid spacing
-# dim = mask.shape # Dimensions of the mag. distr.
-# count = mask.sum() # Number of pixels with magnetization
-# lam = 1e-6 # Regularisation parameter
-# # Create empty MagData object for the reconstruction:
-# mag_data_rec = MagData(a, (np.zeros(dim), np.zeros(dim), np.zeros(dim)))
-#
-# # Function that returns the phase map for a magnetic configuration x:
-# def F(x):
-# mag_data_rec.set_vector(mask, x)
-# phase = pm.phase_mag_real(a, pj.simple_axis_projection(mag_data_rec), b_0)
-# return phase.reshape(-1)
-#
-# # Cost function which should be minimized:
-# def J(x_i):
-# y_i = F(x_i)
-# term1 = (y_i - y_m)
-# term2 = lam * x_i
-# return np.concatenate([term1, term2])
-#
-# # Reconstruct the magnetization components:
-# x_rec, _ = leastsq(J, np.zeros(3*count))
-# mag_data_rec.set_vector(mask, x_rec)
-# return mag_data_rec
-#
-#def reconstruct_test():
-# product = (kernel.multiply_jacobi_T(projection.multiply_jacobi_T(x))
-# * kernel.multiply_jacobi(projection.multiply_jacobi(x)))
\ No newline at end of file
diff --git a/pyramid/phasemap.py b/pyramid/phasemap.py
index be2b2944e8f287a462e068379b22598019ce5bc3..08f2b010deb8ab44d1fb460f0b93d1c1170ed0d6 100644
--- a/pyramid/phasemap.py
+++ b/pyramid/phasemap.py
@@ -2,8 +2,6 @@
"""This module provides the :class:`~.PhaseMap` class for storing phase map data."""
-import logging
-
import numpy as np
from numpy import pi
@@ -17,6 +15,8 @@ from numbers import Number
import netCDF4
+import logging
+
class PhaseMap(object):
@@ -24,7 +24,7 @@ class PhaseMap(object):
Represents 2-dimensional phase maps. The phase information itself is stored as a 2-dimensional
matrix in `phase`, but can also be accessed as a vector via `phase_vec`. :class:`~.PhaseMap`
- objects support negation, arithmetic operators (``+``, ``-``, ``*``) and their augmented
+ objects support negation, arithmetic operators (``+``, ``-``, ``*``) and their augmented
counterparts (``+=``, ``-=``, ``*=``), with numbers and other :class:`~.PhaseMap`
objects, if their dimensions and grid spacings match. It is possible to load data from NetCDF4
or textfiles or to save the data in these formats. Methods for plotting the phase or a
@@ -51,7 +51,7 @@ class PhaseMap(object):
'''
- log = logging.getLogger(__name__)
+ LOG = logging.getLogger(__name__)
UNITDICT = {u'rad': 1E0,
u'mrad': 1E3,
@@ -92,10 +92,10 @@ class PhaseMap(object):
(0.50, 1.0, 1.0),
(0.75, 1.0, 1.0),
(1.00, 0.0, 0.0)]}
-
+
HOLO_CMAP = mpl.colors.LinearSegmentedColormap('my_colormap', CDICT, 256)
HOLO_CMAP_INV = mpl.colors.LinearSegmentedColormap('my_colormap', CDICT_INV, 256)
-
+
@property
def a(self):
return self._a
@@ -104,8 +104,8 @@ class PhaseMap(object):
def a(self, a):
assert isinstance(a, Number), 'Grid spacing has to be a number!'
assert a >= 0, 'Grid spacing has to be a positive number!'
- self._a = a
-
+ self._a = float(a)
+
@property
def dim_uv(self):
return self._dim_uv
@@ -141,70 +141,70 @@ class PhaseMap(object):
self._unit = unit
def __init__(self, a, phase, unit='rad'):
+ self.LOG.debug('Calling __init__')
self.a = a
self.phase = phase
self.unit = unit
- self.log = logging.getLogger(__name__)
- self.log.info('Created '+str(self))
+ self.LOG.debug('Created '+str(self))
def __repr__(self):
- self.log.info('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.info('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.info('Calling __neg__')
+ self.LOG.debug('Calling __neg__')
return PhaseMap(self.a, -self.phase, self.unit)
-
+
def __add__(self, other): # self + other
- self.log.info('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.info('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.info('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.info('Calling __sub__')
+ self.LOG.debug('Calling __sub__')
return self.__add__(-other)
def __mul__(self, other): # self * other
- self.log.info('Calling __mul__')
+ self.LOG.debug('Calling __mul__')
assert isinstance(other, Number), 'PhaseMap objects can only be multiplied by numbers!'
return PhaseMap(self.a, other*self.phase, self.unit)
def __radd__(self, other): # other + self
- self.log.info('Calling __radd__')
+ self.LOG.debug('Calling __radd__')
return self.__add__(other)
def __rsub__(self, other): # other - self
- self.log.info('Calling __rsub__')
+ self.LOG.debug('Calling __rsub__')
return -self.__sub__(other)
def __rmul__(self, other): # other * self
- self.log.info('Calling __rmul__')
+ self.LOG.debug('Calling __rmul__')
return self.__mul__(other)
-
+
def __iadd__(self, other): # self += other
- self.log.info('Calling __iadd__')
+ self.LOG.debug('Calling __iadd__')
return self.__add__(other)
def __isub__(self, other): # self -= other
- self.log.info('Calling __isub__')
+ self.LOG.debug('Calling __isub__')
return self.__sub__(other)
def __imul__(self, other): # self *= other
- self.log.info('Calling __imul__')
+ self.LOG.debug('Calling __imul__')
return self.__mul__(other)
def save_to_txt(self, filename='..\output\phasemap_output.txt'):
@@ -221,7 +221,7 @@ class PhaseMap(object):
None
'''
- self.log.info('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))
@@ -242,7 +242,7 @@ class PhaseMap(object):
A :class:`~.PhaseMap` object containing the loaded data.
'''
- cls.log.info('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])
@@ -263,7 +263,7 @@ class PhaseMap(object):
None
'''
- self.log.info('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])
@@ -287,7 +287,7 @@ class PhaseMap(object):
A :class:`~.PhaseMap` object containing the loaded data.
'''
- cls.log.info('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'][:]
@@ -307,7 +307,7 @@ class PhaseMap(object):
The default is 'RdBu'.
limit : float, optional
Plotlimit for the phase in both negative and positive direction (symmetric around 0).
- If not specified, the maximum amplitude the phase is used.
+ If not specified, the maximum amplitude of the phase is used.
norm : :class:`~matplotlib.colors.Normalize` or subclass, optional
Norm, which is used to determine the colors to encode the phase information.
If not specified, :class:`~matplotlib.colors.Normalize` is automatically used.
@@ -322,7 +322,7 @@ class PhaseMap(object):
The axis on which the graph is plotted.
'''
- self.log.info('Calling display_phase')
+ self.LOG.debug('Calling display_phase')
# Take units into consideration:
phase = self.phase * self.UNITDICT[self.unit]
if limit is None:
@@ -374,7 +374,7 @@ class PhaseMap(object):
The axis on which the graph is plotted.
'''
- self.log.info('Calling display_phase3d')
+ self.LOG.debug('Calling display_phase3d')
# Take units into consideration:
phase = self.phase * self.UNITDICT[self.unit]
# Create figure and axis:
@@ -395,11 +395,11 @@ class PhaseMap(object):
plt.show()
# Return plotting axis:
return axis
-
+
def display_holo(self, density=1, title='Holographic Contour Map',
- axis=None, grad_encode='dark', interpolation='none', show=True):
+ axis=None, grad_encode='bright', interpolation='none', show=True):
'''Display the color coded holography image.
-
+
Parameters
----------
density : float, optional
@@ -408,42 +408,51 @@ class PhaseMap(object):
The title of the plot. The default is 'Holographic Contour Map'.
axis : :class:`~matplotlib.axes.AxesSubplot`, optional
Axis on which the graph is plotted. Creates a new figure if none is specified.
+ grad_encode: {'bright', 'dark', 'color', 'none'}, optional
+ Encoding mode of the phase gradient. 'none' produces a black-white image, 'color' just
+ 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.
interpolation : {'none, 'bilinear', 'cubic', 'nearest'}, optional
Defines the interpolation method. No interpolation is used in the default case.
show: bool, optional
A switch which determines if the plot is shown at the end of plotting.
-
+
Returns
-------
axis: :class:`~matplotlib.axes.AxesSubplot`
The axis on which the graph is plotted.
-
- '''# TODO: Docstring saturation!
- self.log.info('Calling display_holo')
+
+ '''
+ self.LOG.debug('Calling display_holo')
# Calculate the holography image intensity:
img_holo = (1 + np.cos(density * self.phase)) / 2
# Calculate the phase gradients, expressed by magnitude and angle:
phase_grad_y, phase_grad_x = np.gradient(self.phase, self.a, self.a)
phase_angle = (1 - np.arctan2(phase_grad_y, phase_grad_x)/pi) / 2
phase_magnitude = np.hypot(phase_grad_x, phase_grad_y)
- if phase_magnitude.max() != 0:
+ if phase_magnitude.max() != 0: # Take care of phase maps with only zeros
saturation = np.sin(phase_magnitude/phase_magnitude.max() * pi / 2)
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')
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')
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')
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')
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:
- raise AssertionError('Gradient encoding not recognized!')
+ raise AssertionError('Gradient encoding not recognized!')
holo_image = Image.fromarray(rgb)
# If no axis is specified, a new figure is created:
if axis is None:
@@ -469,26 +478,41 @@ class PhaseMap(object):
return axis
def display_combined(self, title='Combined Plot', cmap='RdBu', limit=None, norm=None,
- density=1, interpolation='none', grad_encode='dark'):
+ density=1, interpolation='none', grad_encode='bright'):
'''Display the phase map and the resulting color coded holography image in one plot.
-
+
Parameters
----------
- density : float, optional
- The gain factor for determining the number of contour lines. The default is 1.
title : string, optional
The title of the plot. The default is 'Combined Plot'.
+ cmap : string, optional
+ The :class:`~matplotlib.colors.Colormap` which is used for the plot as a string.
+ The default is 'RdBu'.
+ limit : float, optional
+ Plotlimit for the phase in both negative and positive direction (symmetric around 0).
+ If not specified, the maximum amplitude of the phase is used.
+ norm : :class:`~matplotlib.colors.Normalize` or subclass, optional
+ Norm, which is used to determine the colors to encode the phase information.
+ If not specified, :class:`~matplotlib.colors.Normalize` is automatically used.
+ density : float, optional
+ The gain factor for determining the number of contour lines in the holographic
+ contour map. The default is 1.
interpolation : {'none, 'bilinear', 'cubic', 'nearest'}, optional
Defines the interpolation method for the holographic contour map.
No interpolation is used in the default case.
-
+ grad_encode: {'bright', 'dark', 'color', 'none'}, optional
+ Encoding mode of the phase gradient. 'none' produces a black-white image, 'color' just
+ 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.
+
Returns
-------
phase_axis, holo_axis: :class:`~matplotlib.axes.AxesSubplot`
The axes on which the graphs are plotted.
-
- '''# TODO: Docstring grad_encode!
- self.log.info('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)
@@ -507,17 +531,17 @@ class PhaseMap(object):
@classmethod
def make_color_wheel(cls):
'''Display a color wheel to illustrate the color coding of the gradient direction.
-
+
Parameters
----------
None
-
+
Returns
-------
None
-
+
'''
- cls.log.info('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)
diff --git a/pyramid/phasemapper.py b/pyramid/phasemapper.py
index 7dd5de87785b357019f28e0b6c22445d9a5fb11d..70c76e427e8fa9b7a2f80afea3542a01d90966df 100644
--- a/pyramid/phasemapper.py
+++ b/pyramid/phasemapper.py
@@ -1,7 +1,5 @@
# -*- coding: utf-8 -*-
-"""Create magnetic and electric phase maps from magnetization data.
-
-This module executes several forward models to calculate the magnetic or electric phase map from
+"""This module executes several forward models to calculate the magnetic or electric phase map from
a given projection of a 3-dimensional magnetic distribution (see :mod:`~pyramid.projector`).
For the magnetic phase map, an approach using real space and one using Fourier space is provided.
The electrostatic contribution is calculated by using the assumption of a mean inner potential.
@@ -21,92 +19,145 @@ from pyramid.magdata import MagData
from pyramid.phasemap import PhaseMap
from pyramid.forwardmodel import ForwardModel
+import logging
+class PhaseMapper(object):
+ '''Abstract base class for the phase calculation from a 2-dimensional distribution.
-class PhaseMapper(object):
+ The :class:`~.PhaseeMapper-` class represents a strategy for the phasemapping of a
+ 2-dimensional magnetic distribution with two components onto a scalar phase map.
+ :class:`~.Kernel` is an abstract base class and provides a unified interface which should be
+ subclassed with custom :func:`__init__` and :func:`__call__` functions. Concrete subclasses
+ can be called as a function and take a :class:`~.MagData` object as input and return a
+ :class:`~.PhaseMap` object.
- '''DOCSTRING''' # TODO: Docstring!
+ '''
__metaclass__ = abc.ABCMeta
+ LOG = logging.getLogger(__name__+'.PhaseMapper')
@abc.abstractmethod
def __init__(self, projector):
- '''Docstring: has to be implemented!''' # TODO: Docstring!
-
+ self.LOG.debug('Calling __init__')
+ raise NotImplementedError()
+
@abc.abstractmethod
def __call__(self, mag_data):
- '''Docstring: has to be implemented!''' # TODO: Docstring!
+ self.LOG.debug('Calling __call__')
+ raise NotImplementedError()
class PMAdapterFM(PhaseMapper):
+ '''Class representing a phase mapping strategy adapting the :class:`~.ForwardModel` class.
+
+ The :class:`~.PMAdapterFM` class is an adapter class to incorporate the forward model from the
+ :module:`~.ForwardModel` module without the need of vector input and output. It directly takes
+ :class:`~.MagData` objects and returns :class:`~.PhaseMap` objects.
+
+ Attributes
+ ----------
+ a : float
+ The grid spacing in nm.
+ projector : :class:`~.Projector`
+ Projector which should be used for the projection of the 3-dimensional magnetization
+ distribution.
+ b_0 : float, optional
+ The magnetic induction corresponding to a magnetization `M`\ :sub:`0` in T.
+ The default is 1.
+ geometry : {'disc', 'slab'}, optional
+ Elementary geometry which is used for the phase contribution of one pixel.
+ Default is 'disc'.
+ fwd_model : :class:`~.ForwardModel`
+ The forward model that is constructed from the given information. It is created internally
+ and does not have to be provided by the user.
+
+ '''
+
+ LOG = logging.getLogger(__name__+'.PMAdapterFM')
+
def __init__(self, a, projector, b_0=1, geometry='disc'):
+ self.LOG.debug('Calling __init__')
assert isinstance(projector, Projector), 'Argument has to be a Projector object!'
self.a = a
self.projector = projector
self.fwd_model = ForwardModel([projector], Kernel(a, projector.dim_uv, b_0, geometry))
+ self.LOG.debug('Created '+str(self))
def __call__(self, mag_data):
+ 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!'
- # TODO: test if mag_data fits in all aspects
+ assert mag_data.a == self.a, 'Grid spacing has to match!'
+ assert mag_data.dim == self.projector.dim, 'Dimensions do not match!'
phase_map = PhaseMap(self.a, np.zeros(self.projector.dim_uv))
phase_map.phase_vec = self.fwd_model(mag_data.mag_vec)
return phase_map
+ def __repr__(self):
+ self.LOG.debug('Calling __repr__')
+ return '%s(a=%r, projector=%r, b_0=%r, geometry=%r)' % \
+ (self.__class__, self.a, self.projector, self.b_0, self.geometry)
+
+ def __str__(self):
+ self.LOG.debug('Calling __str__')
+ return 'PMAdapterFM(a=%s, projector=%s, b_0=%s, geometry=%s)' % \
+ (self.a, self.projector, self.b_0, self.geometry)
+
class PMFourier(PhaseMapper):
-
+
+ '''Class representing a phase mapping strategy using a discretization in Fourier space.
+
+ The :class:`~.PMFourier` class represents a phase mapping strategy involving discretization in
+ Fourier space. It utilizes the Fourier transforms, which are inherently calculated in the
+ :class:`~.Kernel` class and directly takes :class:`~.MagData` objects and returns
+ :class:`~.PhaseMap` objects.
+
+ Attributes
+ ----------
+ a : float
+ The grid spacing in nm.
+ projector : :class:`~.Projector`
+ Projector which should be used for the projection of the 3-dimensional magnetization
+ distribution.
+ b_0 : float, optional
+ The magnetic induction corresponding to a magnetization `M`\ :sub:`0` in T.
+ The default is 1.
+ padding : int, optional
+ Factor for the zero padding. The default is 0 (no padding). For a factor of n the number
+ of pixels is increase by ``(1+n)**2``. Padded zeros are cropped at the end.
+
+ '''
+
+ LOG = logging.getLogger(__name__+'.PMFourier')
PHI_0 = -2067.83 # magnetic flux in T*nm²
def __init__(self, a, projector, b_0=1, padding=0):
+ self.LOG.debug('Calling __init__')
assert isinstance(projector, Projector), 'Argument has to be a Projector object!'
self.a = a
self.projector = projector
self.b_0 = b_0
- self.padding = padding
+ self.padding = padding
+ self.LOG.debug('Created '+str(self))
def __call__(self, mag_data):
- '''Calculate the magnetic phase from magnetization data (Fourier space approach).
-
- Parameters
- ----------
- a : float
- The grid spacing in nm.
- projection : tuple (N=3) of :class:`~numpy.ndarray` (N=2)
- The in-plane projection of the magnetization as a tuple, storing the `u`- and `v`-component
- of the magnetization and the thickness projection for the resulting 2D-grid.
- padding : int, optional
- Factor for the zero padding. The default is 0 (no padding). For a factor of n the number
- of pixels is increase by ``(1+n)**2``. Padded zeros are cropped at the end.
- b_0 : float, optional
- The magnetic induction corresponding to a magnetization `M`\ :sub:`0` in T.
- The default is 1.
-
- Returns
- -------
- phase : :class:`~numpy.ndarray` (N=2)
- The phase as a 2-dimensional array.
-
- ''' # TODO: Docstring
+ 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!'
- # TODO: test if mag_data fits in all aspects (especially with projector)
+ assert mag_data.a == self.a, 'Grid spacing has to match!'
+ assert mag_data.dim == self.projector.dim, 'Dimensions do not match!'
v_dim, u_dim = self.projector.dim_uv
u_mag, v_mag = self.projector(mag_data.mag_vec).reshape((2,)+self.projector.dim_uv)
# Create zero padded matrices:
- u_pad = u_dim/2 * self.padding
- v_pad = v_dim/2 * self.padding
- # TODO: use mag_data.padding or np.pad(...)
- u_mag_big = np.zeros(((1 + self.padding) * v_dim, (1 + self.padding) * u_dim))
- v_mag_big = np.zeros(((1 + self.padding) * v_dim, (1 + self.padding) * u_dim))
- u_mag_big[v_pad:v_pad+v_dim, u_pad:u_pad+u_dim] = u_mag
- v_mag_big[v_pad:v_pad+v_dim, u_pad:u_pad+u_dim] = v_mag
+ 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)
# Fourier transform of the two components:
- u_mag_fft = np.fft.fftshift(np.fft.rfft2(u_mag_big), axes=0)
- v_mag_fft = np.fft.fftshift(np.fft.rfft2(v_mag_big), axes=0)
+ 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)
# Calculate the Fourier transform of the phase:
f_nyq = 0.5 / self.a # nyquist frequency
f_u = np.linspace(0, f_nyq, u_mag_fft.shape[1])
@@ -115,66 +166,67 @@ class PMFourier(PhaseMapper):
coeff = (1j*self.b_0) / (2*self.PHI_0)
phase_fft = coeff*self.a * (u_mag_fft*f_vv - v_mag_fft*f_uu) / (f_uu**2 + f_vv**2 + 1e-30)
# Transform to real space and revert padding:
- phase_big = np.fft.irfft2(np.fft.ifftshift(phase_fft, axes=0))
- phase = phase_big[v_pad:v_pad+v_dim, u_pad:u_pad+u_dim]
+ phase_pad = np.fft.irfft2(np.fft.ifftshift(phase_fft, axes=0))
+ phase = phase_pad[v_pad:v_pad+v_dim, u_pad:u_pad+u_dim]
return PhaseMap(self.a, phase)
+ def __repr__(self):
+ self.LOG.debug('Calling __repr__')
+ return '%s(a=%r, projector=%r, b_0=%r, padding=%r)' % \
+ (self.__class__, self.a, self.projector, self.b_0, self.padding)
+
+ def __str__(self):
+ self.LOG.debug('Calling __str__')
+ return 'PMFourier(a=%s, projector=%s, b_0=%s, padding=%s)' % \
+ (self.a, self.projector, self.b_0, self.padding)
class PMElectric(PhaseMapper):
-
+
+ '''Class representing a phase mapping strategy for the electrostatic contribution.
+
+ The :class:`~.PMElectric` class represents a phase mapping strategy for the electrostatic
+ contribution to the electron phase shift which results e.g. from the mean inner potential in
+ certain samples and which is sensitive to properties of the electron microscope. It directly
+ takes :class:`~.MagData` objects and returns :class:`~.PhaseMap` objects.
+
+ Attributes
+ ----------
+ a : float
+ The grid spacing in nm.
+ projector : :class:`~.Projector`
+ Projector which should be used for the projection of the 3-dimensional magnetization
+ distribution.
+ v_0 : float, optional
+ The mean inner potential of the specimen in V. The default is 1.
+ v_acc : float, optional
+ The acceleration voltage of the electron microscope in V. The default is 30000.
+ threshold : float, optional
+ Threshold for the recognition of the specimen location. The default is 0, which means that
+ all voxels with non-zero magnetization will contribute. Should be above noise level.
+
+ '''
+
+ LOG = logging.getLogger(__name__+'.PMElectric')
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
-
+
def __init__(self, a, projector, v_0=1, v_acc=30000, threshold=0):
- '''Calculate the electric phase from magnetization distributions.
-
- Parameters
- ----------
- a : float
- The grid spacing in nm.
- projection : tuple (N=3) of :class:`~numpy.ndarray` (N=2)
- The in-plane projection of the magnetization as a tuple, storing the u- and v-component
- of the magnetization and the thickness projection for the resulting 2D-grid.
- v_0 : float, optional
- The mean inner potential of the specimen in V. The default is 1.
- v_acc : float, optional
- The acceleration voltage of the electron microscope in V. The default is 30000.
-
- Returns
- -------
- phase : :class:`~numpy.ndarray` (N=2)
- The phase as a 2-dimensional array.
-
- ''' # Docstring!
+ self.LOG.debug('Calling __init__')
+ assert isinstance(projector, Projector), 'Argument has to be a Projector object!'
self.a = a
self.projector = projector
self.v_0 = v_0
self.v_acc = v_acc
self.threshold = threshold
+ self.LOG.debug('Created '+str(self))
def __call__(self, mag_data):
- '''Calculate the electric phase from magnetization distributions.
-
- Parameters
- ----------
- a : float
- The grid spacing in nm.
- projection : tuple (N=3) of :class:`~numpy.ndarray` (N=2)
- The in-plane projection of the magnetization as a tuple, storing the u- and v-component
- of the magnetization and the thickness projection for the resulting 2D-grid.
- v_0 : float, optional
- The mean inner potential of the specimen in V. The default is 1.
- v_acc : float, optional
- The acceleration voltage of the electron microscope in V. The default is 30000.
-
- Returns
- -------
- phase : :class:`~numpy.ndarray` (N=2)
- The phase as a 2-dimensional array.
-
- ''' # Docstring!
+ 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 == self.projector.dim, 'Dimensions do not match!'
# Coefficient calculation:
H_BAR, M_E, Q_E, C = self.H_BAR, self.M_E, self.Q_E, self.C
v_0, v_acc = self.v_0, self.v_acc
@@ -187,46 +239,59 @@ class PMElectric(PhaseMapper):
phase = v_0 * Ce * projection * self.a*1E-9
return PhaseMap(self.a, phase)
+ def __repr__(self):
+ self.LOG.debug('Calling __repr__')
+ return '%s(a=%r, projector=%r, v_0=%r, v_acc=%r, threshold=%r)' % \
+ (self.__class__, self.a, self.projector, self.v_0, self.v_acc, self.threshold)
+
+ def __str__(self):
+ self.LOG.debug('Calling __str__')
+ return 'PMElectric(a=%s, projector=%s, v_0=%s, v_acc=%s, threshold=%s)' % \
+ (self.a, self.projector, self.v_0, self.v_acc, self.threshold)
class PMConvolve(PhaseMapper):
- # TODO: Docstring!
-
- def __init__(self, a, projector, b_0=1, threshold=0, geometry='disc'):
- '''Calculate the magnetic phase from magnetization data.
-
- Parameters
- ----------
- a : float
- The grid spacing in nm.
- projection : tuple (N=3) of :class:`~numpy.ndarray` (N=2)
- The in-plane projection of the magnetization as a tuple, storing the `u`- and
- `v`-component of the magnetization and the thickness projection for the resulting
- 2D-grid.
- b_0 : float, optional
- The magnetic induction corresponding to a magnetization `M`\ :sub:`0` in T.
- The default is 1.
- kernel : :class:`~pyramid.kernel.Kernel`, optional
- Specifies the kernel for the convolution with the magnetization data. If none is
- specified, one will be created with `disc` as the default geometry.
-
- Returns
- -------
- phase : :class:`~numpy.ndarray` (N=2)
- The phase as a 2-dimensional array.
-
- '''# Docstring!
+ '''Class representing a phase mapping strategy using real space discretization and Fourier
+ space convolution.
+
+ The :class:`~.PMConvolve` class represents a phase mapping strategy involving discretization in
+ real space. It utilizes the convolution in Fourier space, directly takes :class:`~.MagData`
+ objects and returns :class:`~.PhaseMap` objects.
+
+ Attributes
+ ----------
+ a : float
+ The grid spacing in nm.
+ projector : :class:`~.Projector`
+ Projector which should be used for the projection of the 3-dimensional magnetization
+ distribution.
+ b_0 : float, optional
+ The magnetic induction corresponding to a magnetization `M`\ :sub:`0` in T.
+ The default is 1.
+ geometry : {'disc', 'slab'}, optional
+ Elementary geometry which is used for the phase contribution of one pixel.
+ Default is 'disc'.
+
+ '''
+
+ LOG = logging.getLogger(__name__+'.PMConvolve')
+
+ def __init__(self, a, projector, b_0=1, geometry='disc'):
+ self.LOG.debug('Calling __init__')
assert isinstance(projector, Projector), 'Argument has to be a Projector object!'
self.a = a
self.projector = projector
self.b_0 = b_0
- self.threshold = threshold
self.kernel = Kernel(a, projector.dim_uv, b_0, geometry)
+ self.LOG.debug('Created '+str(self))
def __call__(self, mag_data):
- # Docstring!
+ 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 == self.projector.dim, 'Dimensions do not match!'
# Process input parameters:
- u_mag, v_mag = self.projector(mag_data.mag_vec).reshape((2,)+self.projector.dim_uv)
+ u_mag, v_mag = self.projector(mag_data.mag_vec).reshape((2,)+self.projector.dim_uv)
# Fourier transform the projected magnetisation:
kernel = self.kernel
u_mag_fft = np.fft.rfftn(u_mag, kernel.dim_fft)
@@ -237,35 +302,52 @@ class PMConvolve(PhaseMapper):
# Return the result:
return PhaseMap(self.a, u_phase-v_phase)
+ def __repr__(self):
+ self.LOG.debug('Calling __repr__')
+ return '%s(a=%r, projector=%r, b_0=%r, geometry=%r)' % \
+ (self.__class__, self.a, self.projector, self.b_0, self.geometry)
+
+ def __str__(self):
+ self.LOG.debug('Calling __str__')
+ return 'PMConvolve(a=%s, projector=%s, b_0=%s, geometry=%s)' % \
+ (self.a, self.projector, self.b_0, self.geometry)
+
class PMReal(PhaseMapper):
+ '''Class representing a phase mapping strategy using real space discretization.
+
+ The :class:`~.PMReal` class represents a phase mapping strategy involving discretization in
+ real space. It directly takes :class:`~.MagData` objects and returns :class:`~.PhaseMap`
+ objects.
+
+ Attributes
+ ----------
+ a : float
+ The grid spacing in nm.
+ projector : :class:`~.Projector`
+ Projector which should be used for the projection of the 3-dimensional magnetization
+ distribution.
+ b_0 : float, optional
+ The magnetic induction corresponding to a magnetization `M`\ :sub:`0` in T.
+ The default is 1.
+ threshold : float, optional
+ Threshold determining for which voxels the phase contribution will be calculated. The
+ default is 0, which means that all voxels with non-zero magnetization will contribute.
+ Should be above noise level.
+ geometry : {'disc', 'slab'}, optional
+ Elementary geometry which is used for the phase contribution of one pixel.
+ Default is 'disc'.
+ numcore : boolean, optional
+ Boolean choosing if Cython enhanced routines from the :module:`~.pyramid.numcore` module
+ should be used. Default is True.
+
+ '''
+
+ LOG = logging.getLogger(__name__+'.PMReal')
+
def __init__(self, a, projector, b_0=1, threshold=0, geometry='disc', numcore=True):
- '''Calculate the magnetic phase from magnetization data (pure real space, no FFT-convolution).
-
- Parameters
- ----------
- a : float
- The grid spacing in nm.
- projection : tuple (N=3) of :class:`~numpy.ndarray` (N=2)
- The in-plane projection of the magnetization as a tuple, storing the `u`- and `v`-component
- of the magnetization and the thickness projection for the resulting 2D-grid.
- geometry : {'disc', 'slab'}, optional
- Specifies the elemental geometry to use for the pixel field.
- The default is 'disc', because of the smaller computational overhead.
- b_0 : float, optional
- The magnetic induction corresponding to a magnetization `M`\ :sub:`0` in T.
- The default is 1.
- jacobi : :class:`~numpy.ndarray` (N=2), optional
- Specifies the matrix in which to save the jacobi matrix. The jacobi matrix will not be
- calculated, if no matrix is specified (default), resulting in a faster computation.
-
- Returns
- -------
- phase : :class:`~numpy.ndarray` (N=2)
- The phase as a 2-dimensional array.
-
- ''' # Docstring!
+ self.LOG.debug('Calling __init__')
assert isinstance(projector, Projector), 'Argument has to be a Projector object!'
self.a = a
self.projector = projector
@@ -273,10 +355,14 @@ class PMReal(PhaseMapper):
self.threshold = threshold
self.kernel = Kernel(a, projector.dim_uv, b_0, geometry)
self.numcore = numcore
+ self.LOG.debug('Created '+str(self))
def __call__(self, mag_data):
- # TODO: Docstring
- # Process input parameters:
+ 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 == self.projector.dim, 'Dimensions do not match!'
+ # Process input parameters:
dim_uv = self.projector.dim_uv
threshold = self.threshold
u_mag, v_mag = self.projector(mag_data.mag_vec).reshape((2,)+dim_uv)
@@ -301,3 +387,14 @@ class PMReal(PhaseMapper):
phase -= v_mag[j, i] * v_phase
# Return the phase:
return PhaseMap(self.a, phase)
+
+ def __repr__(self):
+ self.LOG.debug('Calling __repr__')
+ return '%s(a=%r, projector=%r, b_0=%r, threshold=%r, geometry=%r, numcore=%r)' % \
+ (self.__class__, self.a, self.projector, self.b_0,
+ self.threshold, self.geometry, self.numcore)
+
+ def __str__(self):
+ self.LOG.debug('Calling __str__')
+ return 'PMReal(a=%s, projector=%s, b_0=%s, threshold=%s, geometry=%s, numcore=%s)' % \
+ (self.a, self.projector, self.b_0, self.threshold, self.geometry, self.numcore)
diff --git a/pyramid/projector.py b/pyramid/projector.py
index fd57f77416f4cb9c58101520d3d6d822e455ef1d..452bd873f51a3cd00c3673a29854796bc1624785 100644
--- a/pyramid/projector.py
+++ b/pyramid/projector.py
@@ -3,8 +3,6 @@
projections of vector and scalar fields."""
-import logging
-
import numpy as np
from numpy import pi
@@ -14,11 +12,13 @@ import itertools
from scipy.sparse import coo_matrix, csr_matrix
+import logging
+
class Projector(object):
'''Abstract base class representing a projection function.
-
+
The :class:`~.Projector` class represents a projection function for a 3-dimensional
vector- or scalar field onto a 2-dimensional grid. :class:`~.Projector` is an abstract base
class and provides a unified interface which should be subclassed with a custom
@@ -30,6 +30,8 @@ class Projector(object):
Attributes
----------
+ dim : tuple (N=3)
+ Dimensions (z, y, x) of the magnetization distribution.
dim_uv : tuple (N=2)
Dimensions (v, u) of the projected grid.
size_3d : int
@@ -46,24 +48,50 @@ class Projector(object):
'''
__metaclass__ = abc.ABCMeta
+ LOG = logging.getLogger(__name__+'.Projector')
@abc.abstractmethod
- def __init__(self, dim_uv, weight, coeff):
- self.log = logging.getLogger(__name__)
- self.log.info('Calling __init__')
+ def __init__(self, dim, dim_uv, weight, coeff):
+ self.LOG.debug('Calling __init__')
+ self.dim = dim
self.dim_uv = dim_uv
self.weight = weight
self.coeff = coeff
self.size_2d, self.size_3d = weight.shape
- self.log.info('Created '+str(self))
+ self.LOG.debug('Created '+str(self))
+
+ def __repr__(self):
+ 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__')
+ return 'Projector(dim=%s, dim_uv=%s, coeff=%s)' % (self.dim, self.dim_uv, self.coeff)
def __call__(self, vector):
- self.log.info('Calling as function')
-# print 'Projector - __call__:', len(vector)
+ self.LOG.debug('Calling as function')
return self.jac_dot(vector)
+ @abc.abstractmethod
+ def get_info(self):
+ '''Get specific information about the projector as a string.
+
+ Paramters
+ ---------
+ None
+
+ Returns
+ -------
+ info : string
+ Information about the projector as a string, e.g. for the use in plot titles.
+
+ '''
+ raise NotImplementedError()
+
+
def _vector_field_projection(self, vector):
- self.log.info('Calling _vector_field_projection')
+ self.LOG.debug('Calling _vector_field_projection')
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):
@@ -82,7 +110,7 @@ class Projector(object):
return result
def _vector_field_projection_T(self, vector):
- self.log.info('Calling _vector_field_projection_T')
+ self.LOG.debug('Calling _vector_field_projection_T')
size_2d, size_3d = self.size_2d, self.size_3d
result = np.zeros(3*size_3d)
# Go over all possible component projections (u, v) to (z, y, x):
@@ -101,11 +129,11 @@ class Projector(object):
return result
def _scalar_field_projection(self, vector):
- self.log.info('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.info('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):
@@ -124,91 +152,180 @@ class Projector(object):
always`size_2d`.
'''
- self.log.info('Calling jac_dot')
-# print 'Projector - jac_dot:', len(vector)
+ self.LOG.debug('Calling jac_dot')
if len(vector) == 3*self.size_3d: # mode == 'vector'
- self.log.info('mode == vector')
+ self.LOG.debug('mode == vector')
return self._vector_field_projection(vector)
elif len(vector) == self.size_3d: # mode == 'scalar'
- self.log.info('mode == scalar')
+ self.LOG.debug('mode == scalar')
return self._scalar_field_projection(vector)
else:
raise AssertionError('Vector size has to be suited either for ' \
'vector- or scalar-field-projection!')
def jac_T_dot(self, vector):
- # TODO: Docstring!
- self.log.info('Calling jac_T_dot')
-# print 'Projector - jac_T_dot:', len(vector)
+ '''Multiply a `vector` with the transp. jacobi matrix of this :class:`~.Projector` object.
+
+ Parameters
+ ----------
+ vector : :class:`~numpy.ndarray` (N=1)
+ Vector containing the field which should be projected. Must have the same or 2 times
+ the size of `size_2d` of the projector for scalar and vector projection, respectively.
+
+ Returns
+ -------
+ proj_vector : :class:`~numpy.ndarray` (N=1)
+ Vector containing the multiplication of the input with the transposed jacobi matrix
+ of the :class:`~.Projector` object.
+
+ '''
+ self.LOG.debug('Calling jac_T_dot')
if len(vector) == 2*self.size_2d: # mode == 'vector'
- self.log.info('mode == vector')
+ self.LOG.debug('mode == vector')
return self._vector_field_projection_T(vector)
elif len(vector) == self.size_2d: # mode == 'scalar'
- self.log.info('mode == scalar')
+ self.LOG.debug('mode == scalar')
return self._scalar_field_projection_T(vector)
else:
raise AssertionError('Vector size has to be suited either for ' \
'vector- or scalar-field-projection!')
+class XTiltProjector(Projector):
+
+ '''Class representing a projection function with a tilt around the x-axis.
+
+ The :class:`~.XTiltProjector` class represents a projection function for a 3-dimensional
+ vector- or scalar field onto a 2-dimensional grid, which is a concrete subclass of
+ :class:`~.Projector`.
+
+ Attributes
+ ----------
+ dim : tuple (N=3)
+ Dimensions (z, y, x) of the magnetization distribution.
+ tilt : float
+ Angle in `rad` describing the tilt of the beam direction relative to the x-axis.
+
+ '''
+
+ LOG = logging.getLogger(__name__+'.XTiltProjector')
+
+ def __init__(self, dim, tilt):
+
+ def get_position(p, m, b, size):
+ self.LOG.debug('Calling get_position')
+ y, x = np.array(p)[:, 0]+0.5, np.array(p)[:, 1]+0.5
+ return (y-m*x-b)/np.sqrt(m**2+1) + size/2.
+
+ def get_impact(pos, r, size):
+ self.LOG.debug('Calling get_impact')
+ return [x for x in np.arange(np.floor(pos-r), np.floor(pos+r)+1, dtype=int)
+ if 0 <= x < size]
+
+ def get_weight(delta, rho): # use circles to represent the voxels
+ self.LOG.debug('Calling get_weight')
+ lo, up = delta-rho, delta+rho
+ # Upper boundary:
+ if up >= 1:
+ w_up = 0.5
+ else:
+ w_up = (up*np.sqrt(1-up**2) + np.arctan(up/np.sqrt(1-up**2))) / pi
+ # Lower boundary:
+ if lo <= -1:
+ w_lo = -0.5
+ else:
+ 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.tilt = tilt
+ # Set starting variables:
+ # length along projection (proj, z), perpendicular (perp, y) and rotation (rot, x) axis:
+ dim_proj, dim_perp, dim_rot = dim
+ size_2d = dim[1] * dim[2] # x-y-plane
+ size_3d = np.prod(dim)
+ # Creating coordinate list of all voxels:
+ voxels = list(itertools.product(range(dim_proj), range(dim_perp))) # z-y-plane
+ # Calculate positions along the projected pixel coordinate system:
+ center = (dim_proj/2., dim_perp/2.)
+ m = np.where(tilt<=pi, -1/np.tan(tilt+1E-30), 1/np.tan(tilt+1E-30))
+ b = center[0] - m * center[1]
+ positions = get_position(voxels, m, b, dim_perp)
+ # Calculate weight-matrix:
+ r = 1/np.sqrt(np.pi) # radius of the voxel circle
+ rho = 0.5 / r
+ row = []
+ col = []
+ data = []
+ # one slice:
+ for i, voxel in enumerate(voxels):
+ impacts = get_impact(positions[i], r, dim_perp) # impact along y axis
+ for impact in impacts:
+ distance = np.abs(impact+0.5 - positions[i])
+ delta = distance / r
+ col.append(voxel[0]*size_2d + voxel[1]*dim_rot) # 0: z, 1: y
+ row.append(impact*dim_rot)
+ data.append(get_weight(delta, rho))
+ # All other slices (along x):
+ columns = col
+ rows = row
+ for i in np.arange(1, dim_rot):
+ columns = np.hstack((np.array(columns), np.array(col)+i))
+ rows = np.hstack((np.array(rows), np.array(row)+i))
+ # Calculate weight matrix and coefficients for jacobi matrix:
+ weight = csr_matrix(coo_matrix((np.tile(data, dim_rot), (rows, columns)),
+ shape = (size_2d, size_3d)))
+ dim_v, dim_u = dim_perp, dim_rot
+ coeff = [[1, 0, 0], [0, np.cos(tilt), np.sin(tilt)]]
+ super(XTiltProjector, self).__init__(dim, (dim_v, dim_u), weight, coeff)
+ self.LOG.debug('Created '+str(self))
+
+ def get_info(self):
+ '''Get specific information about the projector as a string.
+
+ Paramters
+ ---------
+ None
+
+ Returns
+ -------
+ info : string
+ Information about the projector as a string, e.g. for the use in plot titles.
+
+ '''
+ return 'x-tilt: $(\phi = {:2.1f} \pi)$'.format(self.tilt)
+
class YTiltProjector(Projector):
-# '''Class representing a projection where the projection axis is tilted around the y-axis.
-#
-# The :class:`~.YTiltProjector` class is a concrete subclass of the :class:`~.Projector` class
-# and overwrites the :func:`__init__` constructor which accepts `dim` and `tilt` as arguments.
-# The dimensions of the 3-dimensional grid are given by `dim`, the tilting angle of the ample
-# around the y-axis is given by `tilt`.
-#
-# Attributes
-# ----------
-# dim_uv : tuple (N=2)
-# Dimensions (v, u) of the projected grid.
-# size_3d : int
-# Number of voxels of the 3-dimensional grid.
-# size_2d : int
-# Number of pixels of the 2-dimensional projected grid.
-# weight : :class:`~scipy.sparse.csr_matrix` (N=2)
-# The weight matrix containing the weighting coefficients describing the influence of all
-# 3-dimensional voxels on the 2-dimensional pixels of the projection.
-# coeff : list (N=2)
-# List containing the six weighting coefficients describing the influence of the 3 components
-# of a 3-dimensional vector field on the 2 projected components.
-# '''
-#
-# '''
-# weight : :class:`~scipy.sparse.csr_matrix` (N=2)
-# The weight matrix containing the weighting coefficients which determine the influence of
-# all 3-dimensional voxels to the 2-dimensional pixels of the projection.
-# coeff : `list`t (N=2)
-# List containing the six weighting coefficients describing the influence of the 3 components
-# of a 3-dimensional vector fields on the 2 components of the projected field. Only used for
-# vector field projection.
-#
-# Notes
-# -----
-# An instance `projector` of the :class:`~.YSimpleProjector` class is callable via:
-#
-# :func:`projector(vector)`
-#
-# with `vector` being a :class:`~numpy.ndarray` (N=1).
-# '''
+ '''Class representing a projection function with a tilt around the y-axis.
+
+ The :class:`~.YTiltProjector` class represents a projection function for a 3-dimensional
+ vector- or scalar field onto a 2-dimensional grid, which is a concrete subclass of
+ :class:`~.Projector`.
+
+ Attributes
+ ----------
+ dim : tuple (N=3)
+ Dimensions (z, y, x) of the magnetization distribution.
+ tilt : float
+ Angle in `rad` describing the tilt of the beam direction relative to the y-axis.
+
+ '''
+
+ LOG = logging.getLogger(__name__+'.YTiltProjector')
def __init__(self, dim, tilt):
def get_position(p, m, b, size):
- self.log.info('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.info('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.info('Calling get_weight')
lo, up = delta-rho, delta+rho
# Upper boundary:
if up >= 1:
@@ -222,16 +339,15 @@ 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 = logging.getLogger(__name__)
- self.log.info('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:
dim_proj, dim_rot, dim_perp = dim
- size_2d = dim_rot * dim_perp
- size_3d = dim_proj * dim_rot * dim_perp
+ size_2d = dim[1] * dim[2] # x-y-plane
+ size_3d = np.prod(dim)
# Creating coordinate list of all voxels:
- voxels = list(itertools.product(range(dim_proj), range(dim_perp)))
+ voxels = list(itertools.product(range(dim_proj), range(dim_perp))) # z-x-plane
# Calculate positions along the projected pixel coordinate system:
center = (dim_proj/2., dim_perp/2.)
m = np.where(tilt<=pi, -1/np.tan(tilt+1E-30), 1/np.tan(tilt+1E-30))
@@ -245,17 +361,17 @@ class YTiltProjector(Projector):
data = []
# one slice:
for i, voxel in enumerate(voxels):
- impacts = get_impact(positions[i], r, dim_perp)
+ impacts = get_impact(positions[i], r, dim_perp) # impact along x axis
for impact in impacts:
distance = np.abs(impact+0.5 - positions[i])
delta = distance / r
- col.append(voxel[0]*size_2d + voxel[1])
+ col.append(voxel[0]*size_2d + voxel[1]) # 0: z, 1: x
row.append(impact)
data.append(get_weight(delta, rho))
- # All other slices:
+ # All other slices (along y):
columns = col
rows = row
- for i in np.arange(1, dim_rot): # TODO: more efficient, please!
+ for i in np.arange(1, dim_rot):
columns = np.hstack((np.array(columns), np.array(col)+i*dim_perp))
rows = np.hstack((np.array(rows), np.array(row)+i*dim_perp))
# Calculate weight matrix and coefficients for jacobi matrix:
@@ -263,50 +379,49 @@ class YTiltProjector(Projector):
shape = (size_2d, size_3d)))
dim_v, dim_u = dim_rot, dim_perp
coeff = [[np.cos(tilt), 0, np.sin(tilt)], [0, 1, 0]]
- super(YTiltProjector, self).__init__((dim_v, dim_u), weight, coeff)
- self.log.info('Created '+str(self))
+ super(YTiltProjector, self).__init__(dim, (dim_v, dim_u), weight, coeff)
+ self.LOG.debug('Created '+str(self))
+
+ def get_info(self):
+ '''Get specific information about the projector as a string.
+
+ Paramters
+ ---------
+ None
+
+ Returns
+ -------
+ info : string
+ Information about the projector as a string, e.g. for the use in plot titles.
+
+ '''
+ return 'y-tilt: $(\phi = {:2.1f} \pi)$'.format(self.tilt)
class SimpleProjector(Projector):
-# '''Class representing a projection along one of the major axes (x, y, z).
-#
-# The :class:`~.SimpleProjector` class is a concrete subclass of the :class:`~.Projector` class
-# and overwrites the :func:`__init__` constructor which accepts `dim` and `axis` as arguments.
-# The dimensions of the 3-dimensional grid are given by `dim`, the major axis along which to
-# project is given by `axis` and can be `'x'`, `'y'` or `'z'` (default).
-#
-# Attributes
-# ----------
-# dim_uv : tuple (N=2)
-# Dimensions (v, u) of the projected grid.
-# weight : :class:`~scipy.sparse.csr_matrix` (N=2)
-# The weight matrix containing the weighting coefficients which determine the influence of
-# all 3-dimensional voxels to the 2-dimensional pixels of the projection.
-# coeff : list (N=2)
-# List containing the six weighting coefficients describing the influence of the 3 components
-# of a 3-dimensional vector fields on the 2 components of the projected field. Only used for
-# vector field projection.
-# size_3d : int
-# Number of voxels of the 3-dimensional grid.
-# size_2d : int
-# Number of pixels of the 2-dimensional projected grid.
-#
-# Notes
-# -----
-# An instance `projector` of the :class:`~.SimpleProjector` class is callable via:
-#
-# :func:`projector(vector)`
-#
-# with `vector` being a :class:`~numpy.ndarray` (N=1).
-#
-# '''
+ '''Class representing a projection function along one of the major axes.
+
+ The :class:`~.SimpleProjector` class represents a projection function for a 3-dimensional
+ vector- or scalar field onto a 2-dimensional grid, which is a concrete subclass of
+ :class:`~.Projector`.
+
+ Attributes
+ ----------
+ dim : tuple (N=3)
+ Dimensions (z, y, x) of the magnetization distribution.
+ axis : {'z', 'y', 'x'}, optional
+ Main axis along which the magnetic distribution is projected (given as a string). Defaults
+ to the z-axis.
+
+ '''
+ LOG = logging.getLogger(__name__+'.SimpleProjector')
AXIS_DICT = {'z': (0, 1, 2), 'y': (1, 0, 2), 'x': (1, 2, 0)}
def __init__(self, dim, axis='z'):
- self.log = logging.getLogger(__name__)
- self.log.info('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]
dim_z, dim_y, dim_x = dim
@@ -315,17 +430,35 @@ class SimpleProjector(Projector):
data = np.repeat(1, size_3d)
indptr = np.arange(0, size_3d+1, dim_proj)
if axis == 'z':
+ 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)
+ 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')
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')
coeff = [[0, 1, 0], [0, 0, 1]]
indices = np.array([np.arange(dim_proj) + row*dim_proj
for row in range(size_2d)]).reshape(-1)
weight = csr_matrix((data, indices, indptr), shape = (size_2d, size_3d))
- super(SimpleProjector, self).__init__((dim_v, dim_u), weight, coeff)
- self.log.info('Created '+str(self))
+ super(SimpleProjector, self).__init__(dim, (dim_v, dim_u), weight, coeff)
+ self.LOG.debug('Created '+str(self))
+
+ def get_info(self):
+ '''Get specific information about the projector as a string.
+
+ Paramters
+ ---------
+ None
+
+ Returns
+ -------
+ info : string
+ Information about the projector as a string, e.g. for the use in plot titles.
+
+ '''
+ return 'projected along {}-axis'.format(self.tilt)
diff --git a/pyramid/reconstruction.py b/pyramid/reconstruction.py
new file mode 100644
index 0000000000000000000000000000000000000000..fb8a20f8c9d2911f06ea928a865d9cf50cad9661
--- /dev/null
+++ b/pyramid/reconstruction.py
@@ -0,0 +1,117 @@
+# -*- coding: utf-8 -*-
+"""Reconstruct magnetic distributions from given phasemaps.
+
+This module reconstructs 3-dimensional magnetic distributions (as :class:`~pyramid.magdata.MagData`
+objects) from a given set of phase maps (represented by :class:`~pyramid.phasemap.PhaseMap`
+objects) by using several model based reconstruction algorithms which use the forward model
+provided by :mod:`~pyramid.projector` and :mod:`~pyramid.phasemapper` and a priori knowledge of
+the distribution.
+
+"""
+
+
+import numpy as np
+
+from scipy.sparse.linalg import cg
+
+from pyramid.kernel import Kernel
+from pyramid.forwardmodel import ForwardModel
+from pyramid.costfunction import Costfunction, CFAdapterScipyCG
+from pyramid.magdata import MagData
+
+import logging
+
+
+LOG = logging.getLogger(__name__)
+
+
+class PrintIterator(object):
+
+ '''Iterator class which is responsible to give feedback during reconstruction iterations.
+
+ Parameters
+ ----------
+ cost : :class:`~.Costfunction`
+ :class:`~.Costfunction` class for outputting the `cost` of the current magnetization
+ distribution. This should decrease per iteration if the algorithm converges and is only
+ printed for a `verbosity` of 2.
+ verbosity : {2, 1, 0}, optional
+ Parameter defining the verposity of the output. `2` is the default and will show the
+ current number of the iteration and the cost of the current distribution. `2` will just
+ show the iteration number and `0` will prevent output all together.
+
+ Notes
+ -----
+ Normally this class should not be used by the user and is instantiated whithin the
+ :mod:`~.reconstruction` module itself.
+
+ '''
+
+ LOG = logging.getLogger(__name__+'.PrintIterator')
+
+ def __init__(self, cost, verbosity=2):
+ 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))
+
+ def __call__(self, xk):
+ self.LOG.debug('Calling __call__')
+ if self.verbosity == 0:
+ return
+ print 'iteration #', self.next(),
+ if self.verbosity > 1:
+ print 'cost =', self.cost(xk)
+ else:
+ print ''
+
+ def __repr__(self):
+ 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__')
+ return 'PrintIterator(cost=%s, verbosity=%s)' % (self.cost, self.verbosity)
+
+ def next(self):
+ self.iteration += 1
+ return self.iteration
+
+
+def optimize_sparse_cg(data, verbosity=2):
+ '''Reconstruct a three-dimensional magnetic distribution from given phase maps via the
+ conjugate gradient optimizaion method :func:`~.scipy.sparse.linalg.cg`.
+
+ Parameters
+ ----------
+ data : :class:`~.DataSet`
+ :class:`~.DataSet` object containing all phase maps in :class:`~.PhaseMap` objects and all
+ projection directions in :class:`~.Projector` objects. These provide the essential
+ information for the reconstruction.
+ verbosity : {2, 1, 0}, optional
+ Parameter defining the verposity of the output. `2` is the default and will show the
+ current number of the iteration and the cost of the current distribution. `2` will just
+ show the iteration number and `0` will prevent output all together.
+
+ Returns
+ -------
+ mag_data : :class:`~pyramid.magdata.MagData`
+ The reconstructed magnetic distribution as a :class:`~.MagData` object.
+
+ '''
+ LOG.debug('Calling optimize_sparse_cg')
+ # Set up necessary objects:
+ y = data.phase_vec
+ kernel = Kernel(data.a, data.dim_uv, data.b_0)
+ fwd_model = ForwardModel(data.projectors, kernel)
+ cost = Costfunction(y, fwd_model, lam=10**-10)
+ # Optimize:
+ A = CFAdapterScipyCG(cost)
+ b = fwd_model.jac_T_dot(None, y)
+ x_opt, info = cg(A, b, callback=PrintIterator(cost, verbosity))
+ # Create and return fitting MagData object:
+ mag_opt = MagData(fwd_model.a, np.zeros((3,)+fwd_model.dim))
+ mag_opt.mag_vec = x_opt
+ return mag_opt
diff --git a/scripts/abstract_pictures.py b/scripts/abstract_pictures.py
new file mode 100644
index 0000000000000000000000000000000000000000..cb20acdbc1e3e3d52ea26d1277887443d2dec8d7
--- /dev/null
+++ b/scripts/abstract_pictures.py
@@ -0,0 +1,115 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Mon Mar 17 14:28:10 2014
+
+@author: Jan
+"""
+
+
+from numpy import pi
+import numpy as np
+import matplotlib.pyplot as plt
+from matplotlib.ticker import FuncFormatter
+
+import pyramid.magcreator as mc
+from pyramid.magdata import MagData
+from pyramid.projector import SimpleProjector, YTiltProjector
+from pyramid.phasemapper import PMConvolve
+
+
+###################################################################################################
+print('Jan')
+
+dim = (128, 128, 128)
+center = (int(dim[0]/2), int(dim[1]/2), int(dim[2]/2))
+radius = dim[0]/4.
+a = 1.
+
+magnitude = mc.create_mag_dist_homog(mc.Shapes.sphere(dim, center, radius), pi/4)
+
+mag_data = MagData(a, magnitude)
+
+projector = SimpleProjector(dim)
+
+phase_map = PMConvolve(a, projector)(mag_data)
+
+axis = phase_map.display_phase()
+axis.tick_params(axis='both', which='major', labelsize=18)
+axis.set_title('Phase map', fontsize=24)
+axis.set_xlabel('x [nm]', fontsize=18)
+axis.set_ylabel('y [nm]', fontsize=18)
+
+axis = phase_map.display_holo(density=20, interpolation='bilinear')
+axis.tick_params(axis='both', which='major', labelsize=18)
+axis.set_title('Magnetic induction map', fontsize=24)
+axis.set_xlabel('x [nm]', fontsize=18)
+axis.set_ylabel('y [nm]', fontsize=18)
+
+mag_data.scale_down(2)
+
+mag_data.quiver_plot()
+axis = plt.gca()
+axis.tick_params(axis='both', which='major', labelsize=18)
+axis.set_title('Magnetization distribution', fontsize=24)
+axis.set_xlabel('x [nm]', fontsize=18)
+axis.set_ylabel('y [nm]', fontsize=18)
+
+phase_map.make_color_wheel()
+
+shape_vort = mc.Shapes.disc((64, 64, 64), (31.5, 31.5, 31.5), 24, 10)
+
+magnitude_vort = mc.create_mag_dist_vortex(shape_vort)
+
+mag_vort = MagData(a, magnitude_vort)
+
+mag_vort.scale_down(2)
+
+mag_vort.quiver_plot()
+
+
+###################################################################################################
+print('Patrick')
+
+a = 10.0 # nm
+b_0 = 3 # T
+
+dim = (128, 128, 128)
+center = (int(dim[0]/2), int(dim[1]/2), int(dim[2]/2)) # in px (z, y, x), index starts with 0!
+
+mag_data = MagData(a, mc.create_mag_dist_homog(mc.Shapes.ellipse(dim, center, (20., 60.), 5.), 0))
+
+tilts = np.array([0., 60.])/180.*pi
+
+projectors = [YTiltProjector(mag_data.dim, tilt) for tilt in tilts]
+phasemappers = [PMConvolve(mag_data.a, proj, b_0) for proj in projectors]
+
+phase_maps = [pm(mag_data) for pm in phasemappers]
+
+axis = phase_maps[0].display_holo(density=1, interpolation='bilinear')
+axis.tick_params(axis='both', which='major', labelsize=18)
+axis.set_title('Magnetic induction map', fontsize=24)
+axis.xaxis.set_major_formatter(FuncFormatter(lambda x, pos: '{:g}'.format(x*a)))
+axis.yaxis.set_major_formatter(FuncFormatter(lambda x, pos: '{:g}'.format(x*a)))
+axis.set_xlabel('x [nm]', fontsize=18)
+axis.set_ylabel('y [nm]', fontsize=18)
+
+axis = phase_maps[0].display_phase(limit=17)
+axis.tick_params(axis='both', which='major', labelsize=18)
+axis.set_title('Phase map', fontsize=24)
+axis.set_xlabel('x [nm]', fontsize=18)
+axis.set_ylabel('y [nm]', fontsize=18)
+
+axis = phase_maps[1].display_holo(density=1, interpolation='bilinear')
+axis.tick_params(axis='both', which='major', labelsize=18)
+axis.set_title('Magnetic induction map', fontsize=24)
+axis.xaxis.set_major_formatter(FuncFormatter(lambda x, pos: '{:g}'.format(x*a)))
+axis.yaxis.set_major_formatter(FuncFormatter(lambda x, pos: '{:g}'.format(x*a)))
+axis.set_xlabel('x [nm]', fontsize=18)
+axis.set_ylabel('y [nm]', fontsize=18)
+
+axis = phase_maps[1].display_phase(limit=17)
+axis.tick_params(axis='both', which='major', labelsize=18)
+axis.set_title('Phase map', fontsize=24)
+axis.set_xlabel('x [nm]', fontsize=18)
+axis.set_ylabel('y [nm]', fontsize=18)
+
diff --git a/scripts/compare methods/logfile.log b/scripts/compare methods/logfile.log
deleted file mode 100644
index e69de29bb2d1d6434b8b29ae775ad8c2e48c5391..0000000000000000000000000000000000000000
diff --git a/scripts/create distributions/logfile.log b/scripts/create distributions/logfile.log
deleted file mode 100644
index e69de29bb2d1d6434b8b29ae775ad8c2e48c5391..0000000000000000000000000000000000000000
diff --git a/scripts/edinburgh_test.py b/scripts/edinburgh_test.py
index 27bf0ccb94e3bd79a3db93a72a97726544988ed8..42fcf2755576e5f0812b118db84777a462739354 100644
--- a/scripts/edinburgh_test.py
+++ b/scripts/edinburgh_test.py
@@ -9,7 +9,7 @@ Created on Tue Jan 28 15:15:08 2014
import numpy as np
from pyramid.magdata import MagData
from pyramid.projector import SimpleProjector
-from pyramid.phasemapper import PMAdapterFM
+from pyramid.phasemapper import PMAdapterFM, PMFourier
from matplotlib.ticker import FuncFormatter
data = np.loadtxt('../output/data from Edinburgh/long_grain_remapped_0p0035.txt', delimiter=',')
@@ -28,17 +28,18 @@ magnitude = np.array((x_mag, y_mag, z_mag))
mag_data = MagData(a, magnitude)
-mag_data.pad(30, 20, 0)
-
-mag_data.scale_up()
-
-mag_data.quiver_plot()
+#mag_data.pad(30, 20, 0)
+#
+#mag_data.scale_up()
+#
+#mag_data.quiver_plot()
#mag_data.quiver_plot3d()
projector = SimpleProjector(mag_data.dim)
phasemapper = PMAdapterFM(mag_data.a, projector)
+phasemapper = PMFourier(mag_data.a, projector, padding = 1)
phase_map = phasemapper(mag_data)
diff --git a/scripts/interactive_setup.py b/scripts/interactive_setup.py
index 888286d3fdf621f8d13ad81eb30102dd8aaa77bf..804a13f0eb68a62b278ba1eccc1ea46b4c906c91 100644
--- a/scripts/interactive_setup.py
+++ b/scripts/interactive_setup.py
@@ -1,14 +1,21 @@
# -*- coding: utf-8 -*-
"""Imports pyramid package and sets working directory to output directory"""
-import pyramid.analytic as an
-import pyramid.holoimage as hi
+
import pyramid.magcreator as mc
-import pyramid.phasemapper as pm
-import pyramid.projector as pj
-import pyramid.reconstructor as rc
-from pyramid.magdata import MagData
-from pyramid.phasemap import PhaseMap
+import pyramid.analytic as an
+import pyramid.reconstruction as rc
+
+from pyramid.magdata import *
+from pyramid.projector import *
+from pyramid.kernel import *
+from pyramid.phasemapper import *
+from pyramid.phasemap import *
+from pyramid.dataset import *
+from pyramid.forwardmodel import *
+from pyramid.costfunction import *
import os
-os.chdir('C:\\Users\\Jan\\Home\\PhD Thesis\\Pyramid\\tests')
+
+
+os.chdir('C:\\Users\\Jan\\Home\\PhD Thesis\\Pyramid\\output')
diff --git a/scripts/logfile.log b/scripts/logfile.log
deleted file mode 100644
index e69de29bb2d1d6434b8b29ae775ad8c2e48c5391..0000000000000000000000000000000000000000
diff --git a/scripts/paper 1/logfile.log b/scripts/paper 1/logfile.log
deleted file mode 100644
index 3bc36a622a0509babe71d2f0c852fbaa75e3b78a..0000000000000000000000000000000000000000
--- a/scripts/paper 1/logfile.log
+++ /dev/null
@@ -1,84 +0,0 @@
-2014-02-18 09:30:30: INFO @ <pyramid.projector>: Calling __init__
-2014-02-18 09:30:30: INFO @ <pyramid.projector>: Calling __init__
-2014-02-18 09:30:30: INFO @ <pyramid.projector>: Created <pyramid.projector.SimpleProjector object at 0x0DC7ECD0>
-2014-02-18 09:30:30: INFO @ <pyramid.projector>: Created <pyramid.projector.SimpleProjector object at 0x0DC7ECD0>
-2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Calling __str__
-2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim_uv=(128, 128))
-2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Calling __str__
-2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim_uv=(128, 128))
-2014-02-18 09:30:30: INFO @ <pyramid.projector>: Calling as function
-2014-02-18 09:30:30: INFO @ <pyramid.projector>: mode == vector
-2014-02-18 09:30:30: INFO @ <pyramid.projector>: Calling _vector_field_projection
-2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Calling __str__
-2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim_uv=(128, 128))
-2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Calling __sub__
-2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Calling __neg__
-2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Calling __str__
-2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim_uv=(128, 128))
-2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Calling __add__
-2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Adding two PhaseMap objects
-2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Calling __str__
-2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim_uv=(128, 128))
-2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Calling display_phase
-2014-02-18 09:30:30: INFO @ <pyramid.projector>: Calling as function
-2014-02-18 09:30:30: INFO @ <pyramid.projector>: mode == vector
-2014-02-18 09:30:30: INFO @ <pyramid.projector>: Calling _vector_field_projection
-2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Calling __str__
-2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim_uv=(128, 128))
-2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Calling __sub__
-2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Calling __neg__
-2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Calling __str__
-2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim_uv=(128, 128))
-2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Calling __add__
-2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Adding two PhaseMap objects
-2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Calling __str__
-2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim_uv=(128, 128))
-2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Calling __isub__
-2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Calling __sub__
-2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Calling __add__
-2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Adding an offset
-2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Calling __str__
-2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim_uv=(128, 128))
-2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Calling display_phase
-2014-02-18 09:30:31: INFO @ <pyramid.projector>: Calling as function
-2014-02-18 09:30:31: INFO @ <pyramid.projector>: mode == vector
-2014-02-18 09:30:31: INFO @ <pyramid.projector>: Calling _vector_field_projection
-2014-02-18 09:30:31: INFO @ <pyramid.phasemap>: Calling __str__
-2014-02-18 09:30:31: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim_uv=(128, 128))
-2014-02-18 09:30:31: INFO @ <pyramid.phasemap>: Calling __sub__
-2014-02-18 09:30:31: INFO @ <pyramid.phasemap>: Calling __neg__
-2014-02-18 09:30:31: INFO @ <pyramid.phasemap>: Calling __str__
-2014-02-18 09:30:31: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim_uv=(128, 128))
-2014-02-18 09:30:31: INFO @ <pyramid.phasemap>: Calling __add__
-2014-02-18 09:30:31: INFO @ <pyramid.phasemap>: Adding two PhaseMap objects
-2014-02-18 09:30:31: INFO @ <pyramid.phasemap>: Calling __str__
-2014-02-18 09:30:31: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim_uv=(128, 128))
-2014-02-18 09:30:31: INFO @ <pyramid.phasemap>: Calling display_phase
-2014-02-18 09:30:31: INFO @ <pyramid.projector>: Calling as function
-2014-02-18 09:30:31: INFO @ <pyramid.projector>: mode == vector
-2014-02-18 09:30:31: INFO @ <pyramid.projector>: Calling _vector_field_projection
-2014-02-18 09:30:32: INFO @ <pyramid.phasemap>: Calling __str__
-2014-02-18 09:30:32: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim_uv=(128, 128))
-2014-02-18 09:30:32: INFO @ <pyramid.phasemap>: Calling __sub__
-2014-02-18 09:30:32: INFO @ <pyramid.phasemap>: Calling __neg__
-2014-02-18 09:30:32: INFO @ <pyramid.phasemap>: Calling __str__
-2014-02-18 09:30:32: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim_uv=(128, 128))
-2014-02-18 09:30:32: INFO @ <pyramid.phasemap>: Calling __add__
-2014-02-18 09:30:32: INFO @ <pyramid.phasemap>: Adding two PhaseMap objects
-2014-02-18 09:30:32: INFO @ <pyramid.phasemap>: Calling __str__
-2014-02-18 09:30:32: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim_uv=(128, 128))
-2014-02-18 09:30:32: INFO @ <pyramid.phasemap>: Calling __isub__
-2014-02-18 09:30:32: INFO @ <pyramid.phasemap>: Calling __sub__
-2014-02-18 09:30:32: INFO @ <pyramid.phasemap>: Calling __add__
-2014-02-18 09:30:32: INFO @ <pyramid.phasemap>: Adding an offset
-2014-02-18 09:30:32: INFO @ <pyramid.phasemap>: Calling __str__
-2014-02-18 09:30:32: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim_uv=(128, 128))
-2014-02-18 09:30:32: INFO @ <pyramid.phasemap>: Calling display_phase
-2014-02-18 09:30:33: INFO @ <pyramid.projector>: Calling __init__
-2014-02-18 09:30:33: INFO @ <pyramid.projector>: Calling __init__
-2014-02-18 09:30:33: INFO @ <pyramid.projector>: Created <pyramid.projector.SimpleProjector object at 0x0DC7EB70>
-2014-02-18 09:30:33: INFO @ <pyramid.projector>: Created <pyramid.projector.SimpleProjector object at 0x0DC7EB70>
-2014-02-18 09:30:33: INFO @ <pyramid.phasemap>: Calling __str__
-2014-02-18 09:30:33: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim_uv=(128, 128))
-2014-02-18 09:30:33: INFO @ <pyramid.phasemap>: Calling __str__
-2014-02-18 09:30:33: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim_uv=(128, 128))
diff --git a/scripts/simple_reconstruction.py b/scripts/simple_reconstruction.py
index f4872723c755f40612dbb937a564e9063fe0093f..45385f9abcb421ea380a589913269ef84df93d2a 100644
--- a/scripts/simple_reconstruction.py
+++ b/scripts/simple_reconstruction.py
@@ -10,15 +10,17 @@ import numpy as np
from numpy import pi
from pyramid.magdata import MagData
-from pyramid.projector import YTiltProjector
+from pyramid.projector import YTiltProjector, XTiltProjector
from pyramid.phasemapper import PMConvolve
-from pyramid.datacollection import DataCollection
+from pyramid.dataset import DataSet
+import pyramid.magcreator as mc
from pyramid.kernel import Kernel
from pyramid.forwardmodel import ForwardModel
from pyramid.costfunction import Costfunction
+import pyramid.reconstruction as rc
-from scipy.sparse.linalg import cg
+from scipy.sparse.linalg import cg, LinearOperator
###################################################################################################
@@ -41,7 +43,7 @@ phase_maps = [pm(mag_data) for pm in phasemappers]
dim_uv = dim[1:3]
-data_collection = DataCollection(a, dim_uv, b_0)
+data_collection = DataSet(a, dim_uv, b_0)
[data_collection.append((phase_maps[i], projectors[i])) for i in range(count)]
@@ -100,7 +102,7 @@ phase_maps = [pm(mag_data) for pm in phasemappers]
dim_uv = dim[1:3]
-data_collection = DataCollection(a, dim_uv, b_0)
+data_collection = DataSet(a, dim_uv, b_0)
[data_collection.append((phase_maps[i], projectors[i])) for i in range(count)]
@@ -166,35 +168,40 @@ print('--STARTING RECONSTRUCTION')
###################################################################################################
print('--Generating input phase_maps')
-a = 1.
+a = 10.
b_0 = 1000.
-dim = (9, 9, 9)
-count = 8
+dim = (8, 8, 8)
+count = 32
magnitude = np.zeros((3,)+dim)
-magnitude[0, int(dim[0]/2), int(dim[1]/2), int(dim[2]/2)] = 1.
+magnitude[:, 3:6, 3:6, 3:6] = 1# int(dim[0]/2), int(dim[1]/2), int(dim[2]/2)] = 1.
+magnitude = mc.create_mag_dist_vortex(mc.Shapes.disc(dim, (3.5, 3.5, 3.5), 3, 4))
mag_data = MagData(a, magnitude)
mag_data.quiver_plot3d()
-tilts = np.linspace(0, 2*pi, num=count, endpoint=False)
-projectors = [YTiltProjector(mag_data.dim, tilt) for tilt in tilts]
+tilts = np.linspace(0, 2*pi, num=count/2, endpoint=False)
+projectors = []
+projectors.extend([XTiltProjector(mag_data.dim, tilt) for tilt in tilts])
+projectors.extend([YTiltProjector(mag_data.dim, tilt) for tilt in tilts])
phasemappers = [PMConvolve(mag_data.a, projector, b_0) for projector in projectors]
phase_maps = [pm(mag_data) for pm in phasemappers]
-[phase_map.display_phase(title=u'Tilt series $(\phi = {:2.1f} \pi)$'.format(tilts[i]/pi))
- for i, phase_map in enumerate(phase_maps)]
+#[phase_map.display_phase(title=u'Tilt series $(\phi = {:2.1f} \pi)$'.format(tilts[i%(count/2)]/pi))
+# for i, phase_map in enumerate(phase_maps)]
###################################################################################################
print('--Setting up data collection')
dim_uv = dim[1:3]
+lam = 10. ** -10
+
size_2d = np.prod(dim_uv)
size_3d = np.prod(dim)
-data_collection = DataCollection(a, dim_uv, b_0)
+data_collection = DataSet(a, dim_uv, b_0)
[data_collection.append((phase_maps[i], projectors[i])) for i in range(count)]
@@ -202,37 +209,62 @@ data = data_collection
y = data.phase_vec
kern = Kernel(data.a, data.dim_uv, data.b_0)
F = ForwardModel(data.projectors, kern)
-C = Costfunction(y, F)
+C = Costfunction(y, F, lam)
###################################################################################################
print('--Test simple solver')
-M = np.asmatrix([F.jac_dot(None, np.eye(3*size_3d)[:, i]) for i in range(3*size_3d)]).T
-lam = 1*10. ** -10
-MTM = M.T * M + lam * np.asmatrix(np.eye(3*size_3d))
-
-A = MTM#np.array([F(np.eye(81)[:, i]) for i in range(81)])
-
-b = F.jac_T_dot(None, y)
-
-b_test = np.asarray((M.T.dot(y)).T)[0]
+#M = np.asmatrix([F.jac_dot(None, np.eye(3*size_3d)[:, i]) for i in range(3*size_3d)]).T
+#MTM = M.T * M + lam * np.asmatrix(np.eye(3*size_3d))
+#A = MTM#np.array([F(np.eye(81)[:, i]) for i in range(81)])
-x_f = cg(A, b)[0]
-
-mag_data_rec = MagData(a, np.zeros((3,)+dim))
-
-mag_data_rec.mag_vec = x_f
-
-mag_data_rec.quiver_plot3d()
+#class A_adapt(LinearOperator):
+#
+# def __init__(self, FwdModel, lam, shape):
+# self.fwd = FwdModel
+# self.lam = lam
+# self.shape = shape
+#
+# def matvec(self, vector):
+# return self.fwd.jac_T_dot(None, self.fwd.jac_dot(None, vector)) + self.lam*vector
+#
+# @property
+# def shape(self):
+# return self.shape
+#
+# @property
+# def dtype(self):
+# return np.dtype("d") # None #np.ones(1).dtype
+#
+## TODO: .shape in F und C
+#
+#b = F.jac_T_dot(None, y)
+#
+#A_fast = A_adapt(F, lam, (3*size_3d, 3*size_3d))
+#
+#i = 0
+#def printit(_):
+# global i
+# i += 1
+# print i
+#
+#x_f, info = cg(A_fast, b, callback=printit)
+#
+#mag_data_rec = MagData(a, np.zeros((3,)+dim))
+#
+#mag_data_rec.mag_vec = x_f
+#
+#mag_data_rec.quiver_plot3d()
-phase_maps_rec = [pm(mag_data_rec) for pm in phasemappers]
-[phase_map.display_phase(title=u'Tilt series (rec.) $(\phi = {:2.1f} \pi)$'.format(tilts[i]/pi))
- for i, phase_map in enumerate(phase_maps_rec)]
+#phase_maps_rec = [pm(mag_data_rec) for pm in phasemappers]
+#[phase_map.display_phase(title=u'Tilt series (rec.) $(\phi = {:2.1f} \pi)$'.format(tilts[i%(count/2)]/pi))
+# for i, phase_map in enumerate(phase_maps_rec)]
+mag_data_opt = rc.optimize_sparse_cg(data_collection)
+mag_data_opt.quiver_plot3d()
-#
#first_guess = MagData(a, np.zeros((3,)+dim))
#
#first_guess.magnitude[1, int(dim[0]/2), int(dim[1]/2), int(dim[2]/2)] = 1
@@ -253,6 +285,77 @@ phase_maps_rec = [pm(mag_data_rec) for pm in phasemappers]
#phase_opt.display_phase()
+
+###################################################################################################
+print('--Singular value decomposition')
+
+#a = 1.
+#b_0 = 1000.
+#dim = (3, 3, 3)
+#count = 8
+#
+#magnitude = np.zeros((3,)+dim)
+#magnitude[:, int(dim[0]/2), int(dim[1]/2), int(dim[2]/2)] = 1.
+#
+#mag_data = MagData(a, magnitude)
+#mag_data.quiver_plot3d()
+#
+#tilts = np.linspace(0, 2*pi, num=count/2, endpoint=False)
+#projectors = []
+#projectors.extend([XTiltProjector(mag_data.dim, tilt) for tilt in tilts])
+#projectors.extend([YTiltProjector(mag_data.dim, tilt) for tilt in tilts])
+#phasemappers = [PMConvolve(mag_data.a, projector, b_0) for projector in projectors]
+#phase_maps = [pm(mag_data) for pm in phasemappers]
+#
+##[phase_map.display_phase(title=u'Tilt series $(\phi = {:2.1f} \pi)$'.format(tilts[i%(count/2)]/pi))
+## for i, phase_map in enumerate(phase_maps)]
+#
+#dim_uv = dim[1:3]
+#
+#lam = 10. ** -10
+#
+#size_2d = np.prod(dim_uv)
+#size_3d = np.prod(dim)
+#
+#
+#data_collection = DataCollection(a, dim_uv, b_0)
+#
+#[data_collection.append((phase_maps[i], projectors[i])) for i in range(count)]
+#
+#data = data_collection
+#y = data.phase_vec
+#kern = Kernel(data.a, data.dim_uv, data.b_0)
+#F = ForwardModel(data.projectors, kern)
+#C = Costfunction(y, F, lam)
+#
+#mag_data_opt = opt.optimize_sparse_cg(data_collection)
+#mag_data_opt.quiver_plot3d()
+
+
+
+#M = np.asmatrix([F.jac_dot(None, np.eye(3*size_3d)[:, i]) for i in range(3*size_3d)]).T
+#MTM = M.T * M + lam * np.asmatrix(np.eye(3*size_3d))
+#A = MTM#np.array([F(np.eye(81)[:, i]) for i in range(81)])
+#
+#
+#
+#U, s, V = np.linalg.svd(M)
+##np.testing.assert_almost_equal(U.T, V, decimal=5)
+#
+#for value in range(20):
+# print 'Singular value:', s[value]
+# MagData(data.a, np.array(V[value,:]).reshape((3,)+dim)).quiver_plot3d()
+#
+#
+#for value in range(-10,0):
+# print 'Singular value:', s[value]
+# MagData(data.a, np.array(V[value,:]).reshape((3,)+dim)).quiver_plot3d()
+
+
+# TODO: print all singular vectors for a 2x2x2 distribution for each single and both tilt series!
+
+# TODO: Separate the script for SVD and compliance tests
+
####################################################################################################
#print('--Further testing')
#
@@ -281,7 +384,7 @@ phase_maps_rec = [pm(mag_data_rec) for pm in phasemappers]
#K = np.asmatrix(K)
#lam = 10. ** -10
#KTK = K.T * K + lam * np.asmatrix(np.eye(81))
-#print lam,
+#print lam,
##print pylab.cond(KTK),
#x_f = KTK.I * K.T * y
#print (np.asarray(K * x_f - y) ** 2).sum(),
diff --git a/tests/test_datacollection.py b/tests/test_dataset.py
similarity index 100%
rename from tests/test_datacollection.py
rename to tests/test_dataset.py
diff --git a/tests/test_optimize.py b/tests/test_reconstruction.py
similarity index 79%
rename from tests/test_optimize.py
rename to tests/test_reconstruction.py
index 4cc98625acf36d8591bf5e268343c73dc3aea5a6..7cb6fd3493b09c464564ee6a6e57ec8447502833 100644
--- a/tests/test_optimize.py
+++ b/tests/test_reconstruction.py
@@ -5,10 +5,10 @@
import os
import unittest
-import pyramid.reconstructor as rc
+import pyramid.reconstruction as rc
-class TestCaseReconstructor(unittest.TestCase):
+class TestCaseReconstruction(unittest.TestCase):
def setUp(self):
self.path = os.path.join(os.path.dirname(os.path.realpath(__file__)), 'test_phasemapper')
@@ -21,5 +21,5 @@ class TestCaseReconstructor(unittest.TestCase):
if __name__ == '__main__':
- suite = unittest.TestLoader().loadTestsFromTestCase(TestCaseReconstructor)
+ suite = unittest.TestLoader().loadTestsFromTestCase(TestCaseReconstruction)
unittest.TextTestRunner(verbosity=2).run(suite)