Documentation and cleanup of the pyramid package!

Renames:
datacollection --> dataset
optimizer      --> reconstruction

scripts:
interactive_setup is now implemented in an extended Spyder startup file.
some other scripts are now
the rest is NOT adapted, yet, which is the next task at hand
furthermore scripts will be sorted and unused ones deleted in the next commit

.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

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:

logfile.log

+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>>

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)

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)

pyramid/costfunction.py

 # -*- 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")

pyramid/datacollection.py → pyramid/dataset.py

 # -*- coding: utf-8 -*-
-"""This module provides the :class:`~.DataCollection` class for the collection of phase maps
+"""This module provides the :class:`~.DataSet` class for the collection of phase maps
 and additional data like corresponding projectors."""
 
 
-import logging
-
 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 DataCollection(object):
+
+class DataSet(object):
 
     '''Class for collecting phase maps and corresponding projectors.
 
@@ -37,10 +40,12 @@ class DataCollection(object):
 
     '''
 
+    LOG = logging.getLogger(__name__+'.DataSet')
+
     @property
     def a(self):
         return self._a
-    
+
     @property
     def dim_uv(self):
         return self._dim_uv
@@ -58,27 +63,24 @@ class DataCollection(object):
         return [d[1] for d in self.data]
 
     def __init__(self, a, dim_uv, b_0):
-        self.log = logging.getLogger(__name__)
+        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!'
-        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))
+        self.LOG.debug('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)
+        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.info('Calling __str__')
-        return 'DataCollection(%s, a=%s, dim_uv=%s, data_count=%s)' % \
-            (self.__class__, self.a, self.dim_uv, len(self.data))
+        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.`
@@ -94,9 +96,64 @@ class DataCollection(object):
         None
 
         '''
-        self.log.info('Calling append')
+        self.LOG.debug('Calling append')
         assert isinstance(phase_map, PhaseMap) and isinstance(projector, Projector),  \
             'Argument has to be a tuple of a PhaseMap and a Projector object!'
         assert 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]

pyramid/forwardmodel.py

 # -*- 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)

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

pyramid/logging.ini

@@ -23,7 +23,7 @@ args=tuple()
 
 [handler_file]
 class=logging.FileHandler
-level=INFO
+level=DEBUG
 formatter=file
 args=('logfile.log', 'w')
 

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])

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')

pyramid/numcore/__init__.py

+# -*- 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

pyramid/numcore/kernel_core.pyx

 # -*- 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

pyramid/optimizer.py

-# -*- 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

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)

pyramid/reconstruction.py

+# -*- 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