pyramid/__init__.py

-# -*- coding: utf-8 -*-
-# Copyright 2014 by Forschungszentrum Juelich GmbH
-# Author: J. Caron
-#
-"""Package for the creation and reconstruction of magnetic distributions and resulting phase maps.
-
-Modules
--------
-magcreator
-    Create simple magnetic distributions.
-magdata
-    Class for the storage of magnetization data.
-projector
-    Class for projecting given magnetization distributions.
-kernel
-    Class for the kernel matrix representing one magnetized pixel.
-phasemapper
-    Create magnetic and electric phase maps from magnetization data.
-phasemap
-    Class for the storage of phase data.
-analytic
-    Create phase maps for magnetic distributions with analytic solutions.
-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.
-reconstruction
-    Reconstruct magnetic distributions from given phasemaps.
-regularisator
-    Class to instantiate different regularisation strategies.
-ramp
-    Class which is used to add polynomial ramps to phasemaps.
-diagnostics
-    Class to calculate diagnostics
-quaternion
-    Class which is used for easy rotations in the Projector classes.
-colormap
-    Class which implements a custom direction encoding colormap.
-fft
-    Class for custom FFT functions using numpy or FFTW.
-
-Subpackages
------------
-numcore
-    Provides fast numerical functions for core routines.
-
-"""
-
-
-from . import analytic
-from . import magcreator
-from . import reconstruction
-from . import fft
-from .costfunction import *  # analysis:ignore
-from .dataset import *  # analysis:ignore
-from .diagnostics import *  # analysis:ignore
-from .forwardmodel import *  # analysis:ignore
-from .kernel import *  # analysis:ignore
-from .magdata import *  # analysis:ignore
-from .phasemap import *  # analysis:ignore
-from .phasemapper import *  # analysis:ignore
-from .projector import *  # analysis:ignore
-from .regularisator import *  # analysis:ignore
-from .ramp import *  # analysis:ignore
-from .quaternion import *  # analysis:ignore
-from .colormap import *  # analysis:ignore
-from .config import *  # analysis:ignore
-from .version import version as __version__
-from .version import hg_revision as __hg_revision__
-
-import logging
-_log = logging.getLogger(__name__)
-_log.info("Starting Pyramid V{} HG{}".format(__version__, __hg_revision__))
-del logging
-
-__all__ = ['analytic', 'magcreator', 'reconstruction', 'fft']
-__all__.extend(costfunction.__all__)
-__all__.extend(dataset.__all__)
-__all__.extend(diagnostics.__all__)
-__all__.extend(forwardmodel.__all__)
-__all__.extend(kernel.__all__)
-__all__.extend(magdata.__all__)
-__all__.extend(phasemap.__all__)
-__all__.extend(phasemapper.__all__)
-__all__.extend(projector.__all__)
-__all__.extend(regularisator.__all__)
-__all__.extend(ramp.__all__)
-__all__.extend(quaternion.__all__)
-__all__.extend(colormap.__all__)

pyramid/analytic.py

-# -*- coding: utf-8 -*-
-# Copyright 2014 by Forschungszentrum Juelich GmbH
-# Author: J. Caron
-#
-"""Create phase maps for magnetic distributions with analytic solutions.
-
-This module provides methods for the calculation of the magnetic phase for simple geometries for
-which the analytic solutions are known. These can be used for comparison with the phase
-calculated by the functions from the :mod:`~pyramid.phasemapper` module.
-
-"""
-
-
-import numpy as np
-from numpy import pi
-
-from pyramid.phasemap import PhaseMap
-
-import logging
-
-
-__all__ = ['phase_mag_slab', 'phase_mag_slab', 'phase_mag_sphere', 'phase_mag_vortex']
-_log = logging.getLogger(__name__)
-
-PHI_0 = 2067.83  # magnetic flux in T*nm²
-
-
-def phase_mag_slab(dim, a, phi, center, width, b_0=1):
-    '''Calculate the analytic magnetic phase for a homogeneously magnetized slab.
-
-    Parameters
-    ----------
-    dim : tuple (N=3)
-        The dimensions of the grid `(z, y, x)`.
-    a : float
-        The grid spacing in nm.
-    phi : float
-        The azimuthal angle, describing the direction of the magnetization.
-    center : tuple (N=3)
-        The center of the slab in pixel coordinates `(z, y, x)`.
-    width : tuple (N=3)
-        The width of the slab in pixel coordinates `(z, y, x)`.
-    b_0 : float, optional
-        The magnetic induction corresponding to a magnetization `M`\ :sub:`0` in T.
-        The default is 1.
-
-    Returns
-    -------
-    phase_map : :class:`~numpy.ndarray` (N=2)
-        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):
-            A = np.log(x**2 + y**2 + 1E-30)
-            B = np.arctan(x / (y+1E-30))
-            return x*A - 2*x + 2*y*B
-        return coeff * Lz * (- np.cos(phi) * (F_0(x-x0-Lx/2, y-y0-Ly/2)
-                                              - F_0(x-x0+Lx/2, y-y0-Ly/2)
-                                              - F_0(x-x0-Lx/2, y-y0+Ly/2)
-                                              + F_0(x-x0+Lx/2, y-y0+Ly/2))
-                             + np.sin(phi) * (F_0(y-y0-Ly/2, x-x0-Lx/2)
-                                              - F_0(y-y0+Ly/2, x-x0-Lx/2)
-                                              - F_0(y-y0-Ly/2, x-x0+Lx/2)
-                                              + F_0(y-y0+Ly/2, x-x0+Lx/2)))
-    # Process input parameters:
-    z_dim, y_dim, x_dim = dim
-    y0 = a * center[1]  # y0, x0 define the center of a pixel,
-    x0 = a * center[2]  # hence: (cellindex + 0.5) * grid spacing
-    Lz, Ly, Lx = a * width[0], a * width[1], a * width[2]
-    coeff = - b_0 / (4*PHI_0)  # Minus because of negative z-direction
-    # Create grid:
-    x = np.linspace(a/2, x_dim*a-a/2, num=x_dim)
-    y = np.linspace(a/2, y_dim*a-a/2, num=y_dim)
-    xx, yy = np.meshgrid(x, y)
-    # Return phase:
-    return PhaseMap(a, phi_mag(xx, yy))
-
-
-def phase_mag_disc(dim, a, phi, center, radius, height, b_0=1):
-    '''Calculate the analytic magnetic phase for a homogeneously magnetized disc.
-
-    Parameters
-    ----------
-    dim : tuple (N=3)
-        The dimensions of the grid `(z, y, x)`.
-    a : float
-        The grid spacing in nm.
-    phi : float
-        The azimuthal angle, describing the direction of the magnetization.
-    center : tuple (N=3)
-        The center of the disc in pixel coordinates `(z, y, x)`.
-    radius : float
-        The radius of the disc in pixel coordinates.
-    height : float
-        The height of the disc in pixel coordinates.
-    b_0 : float, optional
-        The magnetic induction corresponding to a magnetization `M`\ :sub:`0` in T.
-        The default is 1.
-
-    Returns
-    -------
-    phase_map : :class:`~numpy.ndarray` (N=2)
-        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)
-        result = coeff * Lz * ((y-y0) * np.cos(phi) - (x-x0) * np.sin(phi))
-        result *= np.where(r <= R, 1, (R / (r+1E-30)) ** 2)
-        return result
-
-    # Process input parameters:
-    z_dim, y_dim, x_dim = dim
-    y0 = a * center[1]
-    x0 = a * center[2]
-    Lz = a * height
-    R = a * radius
-    coeff = pi * b_0 / (2*PHI_0)  # Minus is gone because of negative z-direction
-    # Create grid:
-    x = np.linspace(a/2, x_dim*a-a/2, num=x_dim)
-    y = np.linspace(a/2, y_dim*a-a/2, num=y_dim)
-    xx, yy = np.meshgrid(x, y)
-    # Return phase:
-    return PhaseMap(a, phi_mag(xx, yy))
-
-
-def phase_mag_sphere(dim, a, phi, center, radius, b_0=1):
-    '''Calculate the analytic magnetic phase for a homogeneously magnetized sphere.
-
-    Parameters
-    ----------
-    dim : tuple (N=3)
-        The dimensions of the grid `(z, y, x)`.
-    a : float
-        The grid spacing in nm.
-    phi : float
-        The azimuthal angle, describing the direction of the magnetization.
-    center : tuple (N=3)
-        The center of the sphere in pixel coordinates `(z, y, x)`.
-    radius : float
-        The radius of the sphere in pixel coordinates.
-    b_0 : float, optional
-        The magnetic induction corresponding to a magnetization `M`\ :sub:`0` in T.
-        The default is 1.
-
-    Returns
-    -------
-    phase_map : :class:`~numpy.ndarray` (N=2)
-        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)
-        result = coeff * R ** 3 / (r+1E-30) ** 2 * ((y-y0) * np.cos(phi) - (x-x0) * np.sin(phi))
-        result *= np.where(r > R, 1, (1 - (1 - (r / R) ** 2) ** (3. / 2.)))
-        return result
-
-    # Process input parameters:
-    z_dim, y_dim, x_dim = dim
-    y0 = a * center[1]
-    x0 = a * center[2]
-    R = a * radius
-    coeff = 2./3. * pi * b_0 / PHI_0  # Minus is gone because of negative z-direction
-    # Create grid:
-    x = np.linspace(a / 2, x_dim * a - a / 2, num=x_dim)
-    y = np.linspace(a / 2, y_dim * a - a / 2, num=y_dim)
-    xx, yy = np.meshgrid(x, y)
-    # Return phase:
-    return PhaseMap(a, phi_mag(xx, yy))
-
-
-def phase_mag_vortex(dim, a, center, radius, height, b_0=1):
-    '''Calculate the analytic magnetic phase for a vortex state disc.
-
-    Parameters
-    ----------
-    dim : tuple (N=3)
-        The dimensions of the grid `(z, y, x)`.
-    a : float
-        The grid spacing in nm.
-    center : tuple (N=3)
-        The center of the disc in pixel coordinates `(z, y, x)`, which is also the vortex center.
-    radius : float
-        The radius of the disc in pixel coordinates.
-    height : float
-        The height of the disc in pixel coordinates.
-    b_0 : float, optional
-        The magnetic induction corresponding to a magnetization `M`\ :sub:`0` in T.
-        The default is 1.
-
-    Returns
-    -------
-    phase_map : :class:`~numpy.ndarray` (N=2)
-        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)
-        result = coeff * np.where(r <= R, r - R, 0)
-        return result
-
-    # Process input parameters:
-    z_dim, y_dim, x_dim = dim
-    y0 = a * center[1]
-    x0 = a * center[2]
-    Lz = a * height
-    R = a * radius
-    coeff = - pi * b_0 * Lz / PHI_0  # Minus because of negative z-direction
-    # Create grid:
-    x = np.linspace(a/2, x_dim*a-a/2, num=x_dim)
-    y = np.linspace(a/2, y_dim*a-a/2, num=y_dim)
-    xx, yy = np.meshgrid(x, y)
-    # Return phase:
-    return PhaseMap(a, phi_mag(xx, yy))

pyramid/colormap.py

-# -*- coding: utf-8 -*-
-# Copyright 2014 by Forschungszentrum Juelich GmbH
-# Author: J. Caron
-#
-"""This module provides a custom :class:`~.DirectionalColormap` colormap class which has a few
-additional functions and can encode three-dimensional directions."""
-
-
-import numpy as np
-import matplotlib as mpl
-import matplotlib.pyplot as plt
-from PIL import Image
-from numbers import Number
-
-import logging
-
-
-__all__ = ['DirectionalColormap', 'TransparentColormap']
-
-
-class DirectionalColormap(mpl.colors.LinearSegmentedColormap):
-
-    '''Colormap subclass for encoding 3D-directions with colors.
-
-    This class is a subclass of the :class:`~matplotlib.pyplot.colors.LinearSegmentedColormap`
-    class with a few classmethods which can be used for convenience. The
-    :method:`.display_colorwheel` method can be used to display a colorhweel which is used for
-    the directional encoding and three `rgb_from_...` methods are used to calculate rgb tuples
-    from 3D-directions, angles or directly from a colorindex and a saturation value. This is
-    useful for quiverplots where the arrow colors can be set individually by providing said rgb-
-    tuples. Arrows in plane will be encoded with full color saturation and arrow pointing down
-    will gradually lose saturation until they are black when pointing down. Arrows pointing up
-    will 'oversaturate' (the saturation value will go from 1 up to 2), which in consequence will
-    partially add the inverse colormap to the arrows turning them white if they point up (rgb-
-    tuple: 255, 255, 255).
-
-    Attributes
-    ----------
-    inverted: boolean, optional
-        Flag which is used to invert the internal colormap (default is False).
-        Just influences the classical use as a normal colormap, not the classmethods!
-
-    '''
-
-    _log = logging.getLogger(__name__+'.DirectionalColormap')
-
-    CDICT = {'red': [(0.00, 0.0, 0.0),
-                     (0.25, 1.0, 1.0),
-                     (0.50, 1.0, 1.0),
-                     (0.75, 0.0, 0.0),
-                     (1.00, 0.0, 0.0)],
-
-             'green': [(0.00, 1.0, 1.0),
-                       (0.25, 1.0, 1.0),
-                       (0.50, 0.0, 0.0),
-                       (0.75, 0.0, 0.0),
-                       (1.00, 1.0, 1.0)],
-
-             'blue': [(0.00, 0.0, 0.0),
-                      (0.25, 0.0, 0.0),
-                      (0.50, 0.0, 0.0),
-                      (0.75, 1.0, 1.0),
-                      (1.00, 0.0, 0.0)]}
-
-    CDICT_INV = {'red': [(0.00, 1.0, 1.0),
-                         (0.25, 0.0, 0.0),
-                         (0.50, 0.0, 0.0),
-                         (0.75, 1.0, 1.0),
-                         (1.00, 1.0, 1.0)],
-
-                 'green': [(0.00, 0.0, 0.0),
-                           (0.25, 0.0, 0.0),
-                           (0.50, 1.0, 1.0),
-                           (0.75, 1.0, 1.0),
-                           (1.00, 0.0, 0.0)],
-
-                 'blue': [(0.00, 1.0, 1.0),
-                          (0.25, 1.0, 1.0),
-                          (0.50, 1.0, 1.0),
-                          (0.75, 0.0, 0.0),
-                          (1.00, 1.0, 1.0)]}
-
-    HOLO_CMAP = mpl.colors.LinearSegmentedColormap('my_colormap', CDICT, 256)
-    HOLO_CMAP_INV = mpl.colors.LinearSegmentedColormap('my_colormap', CDICT_INV, 256)
-
-    def __init__(self, inverted=False):
-        self._log.debug('Calling __create_directional_colormap')
-        cdict = self.CDICT_INV if inverted else self.CDICT
-        super(DirectionalColormap, self).__init__('directional_colormap', cdict, N=256)
-        self._log.debug('Created '+str(self))
-
-    @classmethod
-    def display_colorwheel(cls, mode='white_to_color'):
-        '''Display a color wheel to illustrate the color coding of the gradient direction.
-
-        Parameters
-        ----------
-        mode : {'white_to_color', 'color_to_black', 'black_to_color', 'white_to_color_to_black'}
-            Optional string for determining the color scheme of the color wheel. Describes the
-            order of colors from the center to the outline.
-
-        Returns
-        -------
-        None
-
-        '''
-        cls._log.debug('Calling display_color_wheel')
-        yy, xx = np.indices((512, 512)) - 256
-        r = np.hypot(xx, yy)
-        # Create the wheel:
-        colorind = (1 + np.arctan2(yy, xx)/np.pi) / 2
-        saturation = r / 256  # 0 in center, 1 at borders of circle!
-        if mode == 'black_to_color':
-            pass  # no further modification necessary!
-        elif mode == 'color_to_black':
-            saturation = 1 - saturation
-        elif mode == 'white_to_color':
-            saturation = 2 - saturation
-        elif mode == 'white_to_color_to_black':
-            saturation = 2 - 2*saturation
-        else:
-            raise Exception('Invalid color wheel mode!')
-        saturation *= np.where(r <= 256, 1, 0)  # Cut out the wheel!
-        rgb = cls.rgb_from_colorind_and_saturation(colorind, saturation)
-        color_wheel = Image.fromarray(rgb)
-        # Plot the color wheel:
-        fig = plt.figure(figsize=(4, 4))
-        axis = fig.add_subplot(1, 1, 1, aspect='equal')
-        axis.imshow(color_wheel, origin='lower')
-        plt.tick_params(axis='both', which='both', labelleft='off', labelbottom='off',
-                        left='off', right='off', top='off', bottom='off')
-
-    @classmethod
-    def rgb_from_colorind_and_saturation(cls, colorind, saturation):
-        '''Construct a rgb tuple from a colorindex and a saturation value.
-
-        Parameters
-        ----------
-        colorind : float, [0, 1)
-            Colorindex specifying the directional color according to the CDICT colormap.
-            The colormap is periodic so that a value of 1 is equivalent to 0 again.
-        saturation : [0, 2]float, optional
-            Saturation value for the color. The lower hemisphere uses values from 0 to 1 in a
-            traditional sense of saturation with no saturation for directions pointing down (black)
-            and full saturation in plane (full colors). Higher values (between 1 and 2) add the
-            inverse colormap described in CDICT_INV to gradually increase the complementary colors
-            so that arrows pointing up appear white.
-            Azimuthal angle of the desired direction to encode (in rad). Default: pi/2 (in-plane).
-
-        Returns
-        -------
-        rgb : tuple (N=3)
-            rgb tuple with the encoded direction color.
-
-        Notes
-        -----
-            Also works with numpy arrays as input! Always returns array of shape (..., 3)!
-
-        '''
-        cls._log.debug('Calling rgb_from_colorind_and_saturation')
-        # Make sure colorind and saturation are arrays (if not make it so):
-        c = np.array([colorind]) if isinstance(colorind, Number) else colorind
-        s = np.array([saturation]) if isinstance(saturation, Number) else saturation
-        # Calculate rgb and the inverse and combine them:
-        rgb_norm = (np.minimum(s, np.ones_like(s)).T * cls.HOLO_CMAP(c)[..., :3].T).T
-        rgb_invs = (np.maximum(s-1, np.zeros_like(s)).T * cls.HOLO_CMAP_INV(c)[..., :3].T).T
-        return (255.999 * (rgb_norm + rgb_invs)).astype(np.uint8)
-
-    @classmethod
-    def rgb_from_angles(cls, phi, theta=np.pi/2):
-        '''Construct a rgb tuple from two angles representing a 3D direction.
-
-        Parameters
-        ----------
-        phi : float
-            Polar angle of the desired direction to encode (in rad).
-        theta : float, optional
-            Azimuthal angle of the desired direction to encode (in rad). Default: pi/2 (in-plane).
-
-        Returns
-        -------
-        rgb : tuple (N=3)
-            rgb tuple with the encoded direction color.
-
-        Notes
-        -----
-            Also works with numpy arrays as input!
-
-        '''
-        cls._log.debug('Calling rgb_from_angles')
-        colorind = (1 + phi/np.pi) / 2
-        saturation = 2 * (1 - theta/np.pi)  # goes from 0 to 2!
-        return cls.rgb_from_colorind_and_saturation(colorind, saturation)
-
-    @classmethod
-    def rgb_from_direction(cls, x, y, z):
-        '''Construct a rgb tuple from three coordinates representing a 3D direction.
-
-        Parameters
-        ----------
-        x : float
-            x-coordinate of the desired direction to encode.
-        y : float
-            y-coordinate of the desired direction to encode.
-        z : float
-            z-coordinate of the desired direction to encode.
-
-        Returns
-        -------
-        rgb : tuple (N=3)
-            rgb tuple with the encoded direction color.
-
-        Notes
-        -----
-            Also works with numpy arrays as input!
-
-        '''
-        cls._log.debug('Calling rgb_from_direction')
-        phi = np.arctan2(y, x)
-        r = np.sqrt(x**2 + y**2 + z**2)
-        theta = np.arccos(z / (r+1E-30))
-        return cls.rgb_from_angles(phi, theta)
-
-
-class TransparentColormap(mpl.colors.LinearSegmentedColormap):
-
-    '''Colormap subclass for including transparency.
-
-    This class is a subclass of the :class:`~matplotlib.pyplot.colors.LinearSegmentedColormap`
-    class with integrated support for transparency. The colormap is unicolor and varies only in
-    transparency.
-
-    Attributes
-    ----------
-    r: float, optional
-        Intensity of red in the colormap. Has to be between 0. and 1.
-    g: float, optional
-        Intensity of green in the colormap. Has to be between 0. and 1.
-    b: float, optional
-        Intensity of blue in the colormap. Has to be between 0. and 1.
-    alpha_range : list (N=2) of float, optional
-        Start and end alpha value. Has to be between 0. and 1.
-
-    '''
-
-    _log = logging.getLogger(__name__+'.TransparentColormap')
-
-    def __init__(self, r=1., g=0., b=0., alpha_range=[0., 1.]):
-        self._log.debug('Calling __create_directional_colormap')
-        red = [(0., 0., r), (1., r, 1.)]
-        green = [(0., 0., g), (1., g, 1.)]
-        blue = [(0., 0., b), (1., b, 1.)]
-        alpha = [(0., 0., alpha_range[0]), (1., alpha_range[1], 1.)]
-        cdict = {'red': red, 'green': green, 'blue': blue, 'alpha': alpha}
-        super(TransparentColormap, self).__init__('transparent_colormap', cdict, N=256)
-        self._log.debug('Created '+str(self))

pyramid/config.py

-# -*- coding: utf-8 -*-
-# Copyright 2014 by Forschungszentrum Juelich GmbH
-# Author: J. Caron
-#
-"""Config file for the pyramid package."""
-
-
-import os
-
-
-__all__ = ['DIR_PACKAGE', 'DIR_FILES', 'LOGGING_CONFIG', 'NTHREADS']
-
-
-DIR_PACKAGE = os.path.join(os.path.dirname(os.path.realpath(__file__)))
-DIR_FILES = os.path.abspath(os.path.join(DIR_PACKAGE, os.pardir, 'files'))
-LOGGING_CONFIG = os.path.join(os.path.dirname(os.path.realpath(__file__)), 'logging.ini')
-
-try:
-    import multiprocessing
-    NTHREADS = multiprocessing.cpu_count()
-    del multiprocessing
-except ImportError:
-    NTHREADS = 1
-
-del os

pyramid/costfunction.py

-# -*- coding: utf-8 -*-
-# Copyright 2014 by Forschungszentrum Juelich GmbH
-# Author: J. Caron
-#
-"""This module provides the :class:`~.Costfunction` class which represents a strategy to calculate
-the so called `cost` of a threedimensional magnetization distribution."""
-
-
-import numpy as np
-
-from pyramid.regularisator import NoneRegularisator
-
-import logging
-
-
-__all__ = ['Costfunction']
-
-
-class Costfunction(object):
-
-    '''Class for calculating the cost of a 3D magnetic distributions in relation to 2D phase maps.
-
-    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 and from
-    a priori knowledge represented by a :class:`~.Regularisator` object. 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. All required data should be given in a :class:`~DataSet` object.
-
-    Attributes
-    ----------
-    regularisator : :class:`~.Regularisator`
-        Regularisator class that's responsible for the regularisation term.
-    y : :class:`~numpy.ndarray` (N=1)
-        Vector which lists all pixel values of all phase maps one after another.
-    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`.
-    Se_inv : :class:`~numpy.ndarray` (N=2), optional
-        Inverted covariance matrix of the measurement errors. The matrix has size `NxN` with N
-        being the length of the targetvector y (vectorized phase map information).
-    m: int
-        Size of the image space.
-    n: int
-        Size of the input space.
-
-    '''
-
-    _log = logging.getLogger(__name__+'.Costfunction')
-
-    def __init__(self, fwd_model, regularisator):
-        self._log.debug('Calling __init__')
-        self.fwd_model = fwd_model
-        self.regularisator = regularisator
-        if self.regularisator is None:
-            self.regularisator = NoneRegularisator()
-        # Extract important information:
-        self.y = self.fwd_model.data_set.phase_vec
-        self.n = self.fwd_model.n
-        self.m = self.fwd_model.m
-        if self.fwd_model.data_set.Se_inv is None:
-            self.fwd_model.data_set.set_Se_inv_diag_with_conf()
-        self.Se_inv = self.fwd_model.data_set.Se_inv
-        self._log.debug('Created '+str(self))
-
-    def __repr__(self):
-        self._log.debug('Calling __repr__')
-        return '%s(fwd_model=%r, regularisator=%r)' % \
-            (self.__class__, self.fwd_model, self.regularisator)
-
-    def __str__(self):
-        self._log.debug('Calling __str__')
-        return 'Costfunction(fwd_model=%s, fwd_model=%s, regularisator=%s)' % \
-            (self.fwd_model, self.fwd_model, self.regularisator)
-
-    def __call__(self, x):
-        delta_y = self.fwd_model(x) - self.y
-        self.chisq_m = delta_y.dot(self.Se_inv.dot(delta_y))
-        self.chisq_a = self.regularisator(x)
-        self.chisq = self.chisq_m + self.chisq_a
-        return self.chisq
-
-    def init(self, x):
-        '''Initialise the costfunction by calculating the different cost terms.
-
-        Parameters
-        ----------
-        x : :class:`~numpy.ndarray` (N=1)
-            Vectorized magnetization distribution, for which the cost is calculated.
-
-        Returns
-        -------
-        None
-
-        '''
-        self._log.debug('Calling init')
-        self(x)
-
-    def jac(self, x):
-        '''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.
-
-        '''
-        assert len(x) == self.n
-        return (2 * self.fwd_model.jac_T_dot(x, self.Se_inv.dot(self.fwd_model(x) - self.y))
-                + self.regularisator.jac(x))
-
-    def hess_dot(self, x, vector):
-        '''Calculate the product of a `vector` with the Hessian 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.
-
-        '''
-        return (2 * self.fwd_model.jac_T_dot(x, self.Se_inv.dot(self.fwd_model.jac_dot(x, vector)))
-                + self.regularisator.hess_dot(x, vector))
-
-    def hess_diag(self, _):
-        ''' Return the diagonal of the Hessian.
-
-        Parameters
-        ----------
-        _ : undefined
-            Unused input
-
-        '''
-        return np.ones(self.n)

pyramid/dataset.py

-# -*- coding: utf-8 -*-
-# Copyright 2014 by Forschungszentrum Juelich GmbH
-# Author: J. Caron
-#
-"""This module provides the :class:`~.DataSet` class for the collection of phase maps
-and additional data like corresponding projectors."""
-
-
-import numpy as np
-from numbers import Number
-from scipy import sparse
-import matplotlib.pyplot as plt
-
-from pyramid.phasemap import PhaseMap
-from pyramid.phasemapper import PhaseMapperRDFC
-from pyramid.projector import Projector
-from pyramid.kernel import Kernel
-
-import logging
-
-
-__all__ = ['DataSet']
-
-
-class DataSet(object):
-
-    '''Class for collecting phase maps and corresponding projectors.
-
-    Represents a collection of (e.g. experimentally derived) phase maps, stored as
-    :class:`~.PhaseMap` objects and corresponding projectors stored as :class:`~.Projector`
-    objects. At creation, the grid spacing `a` and the dimension `dim` of the magnetization
-    distribution have to be given. Data can be added via the :func:`~.append` method, where
-    a :class:`~.PhaseMap`, a :class:`~.Projector` and additional info have to be given.
-
-    Attributes
-    ----------
-    a: float
-        The grid spacing in nm.
-    dim: tuple (N=3)
-        Dimensions of the 3D magnetization distribution.
-    phase_maps:
-        A list of all stored :class:`~.PhaseMap` objects.
-    b_0: double
-        The saturation induction in `T`.
-    mask: :class:`~numpy.ndarray` (N=3), optional
-        A boolean mask which defines the magnetized volume in 3D.
-    projectors: list of :class:`~.Projector`
-        A list of all stored :class:`~.Projector` objects.
-    phase_maps: list of :class:`~.PhaseMap`
-        A list of all stored :class:`~.PhaseMap` objects.
-    phase_vec: :class:`~numpy.ndarray` (N=1)
-        The concatenaded, vectorized phase of all ;class:`~.PhaseMap` objects.
-    Se_inv : :class:`~numpy.ndarray` (N=2), optional
-        Inverted covariance matrix of the measurement errors. The matrix has size `NxN` with N
-        being the length of the targetvector y (vectorized phase map information).
-    m: int
-        Size of the image space.
-    n: int
-        Size of the input space.
-
-    '''
-
-    _log = logging.getLogger(__name__+'.DataSet')
-
-    @property
-    def a(self):
-        return self._a
-
-    @a.setter
-    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 = float(a)
-
-    @property
-    def mask(self):
-        return self._mask
-
-    @mask.setter
-    def mask(self, mask):
-        if mask is not None:
-            assert mask.shape == self.dim, 'Mask dimensions must match!'
-        else:
-            mask = np.ones(self.dim, dtype=bool)
-        self._mask = mask.astype(np.bool)
-
-    @property
-    def m(self):
-        return np.sum([len(p.phase_vec) for p in self.phase_maps])
-
-    @property
-    def n(self):
-        return 3 * np.sum(self.mask)
-
-    @property
-    def count(self):
-        return len(self.projectors)
-
-    @property
-    def phase_vec(self):
-        return np.concatenate([p.phase_vec for p in self.phase_maps])
-
-    @property
-    def hook_points(self):
-        result = [0]
-        for i, phase_map in enumerate(self.phase_maps):
-            result.append(result[i]+np.prod(phase_map.dim_uv))
-        return result
-
-    @property
-    def phase_mappers(self):
-        dim_uv_set = set([p.dim_uv for p in self.projectors])
-        kernel_list = [Kernel(self.a, dim_uv) for dim_uv in dim_uv_set]
-        return {kernel.dim_uv: PhaseMapperRDFC(kernel) for kernel in kernel_list}
-
-    def __init__(self, a, dim, b_0=1, mask=None, Se_inv=None):
-        self._log.debug('Calling __init__')
-        assert isinstance(dim, tuple) and len(dim) == 3, \
-            'Dimension has to be a tuple of length 3!'
-        self.a = a
-        self.dim = dim
-        self.b_0 = b_0
-        self.mask = mask
-        self.Se_inv = Se_inv
-        self.phase_maps = []
-        self.projectors = []
-        self._log.debug('Created: '+str(self))
-
-    def __repr__(self):
-        self._log.debug('Calling __repr__')
-        return '%s(a=%r, dim=%r, b_0=%r, mask=%r, Se_inv=%r)' % (self.__class__, self.a, self.dim,
-                                                                 self.b_0, self.mask, self.Se_inv)
-
-    def __str__(self):
-        self._log.debug('Calling __str__')
-        return 'DataSet(a=%s, dim=%s, b_0=%s)' % (self.a, self.dim, self.b_0)
-
-    def append(self, phase_map, projector):
-        '''Appends a data pair of phase map and projection infos to the data collection.`
-
-        Parameters
-        ----------
-        phase_map: :class:`~.PhaseMap`
-            A :class:`~.PhaseMap` object which should be added to the data collection.
-        projector: :class:`~.Projector`
-            A :class:`~.Projector` object which should be added to the data collection.
-
-        Returns
-        -------
-        None
-
-        '''
-        self._log.debug('Calling append')
-        assert isinstance(phase_map, PhaseMap) and isinstance(projector, Projector),  \
-            'Argument has to be a tuple of a PhaseMap and a Projector object!'
-        assert projector.dim == self.dim, '3D dimensions must match!'
-        assert phase_map.dim_uv == projector.dim_uv, 'Projection dimensions (dim_uv) must match!'
-        self.phase_maps.append(phase_map)
-        self.projectors.append(projector)
-
-    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.
-
-        '''
-        self._log.debug('Calling create_phase_maps')
-        phase_maps = []
-        for projector in self.projectors:
-            mag_proj = projector(mag_data)
-            phase_map = self.phase_mappers[projector.dim_uv](mag_proj)
-            phase_map.mask = mag_proj.get_mask()[0, ...]
-            phase_maps.append(phase_map)
-        return phase_maps
-
-    def set_Se_inv_block_diag(self, cov_list):
-        '''Set the Se_inv matrix as a block diagonal matrix
-
-        Parameters
-        ----------
-        cov_list: list of :class:`~numpy.ndarray`
-            List of inverted covariance matrices (one for each projection).
-
-        Returns
-        -------
-            None
-
-        '''
-        self._log.debug('Calling set_Se_inv_block_diag')
-        assert len(cov_list) == len(self.phase_maps), 'Needs one covariance matrix per phase map!'
-        self.Se_inv = sparse.block_diag(cov_list).tocsr()
-
-    def set_Se_inv_diag_with_conf(self, conf_list=None):
-        '''Set the Se_inv matrix as a block diagonal matrix from a list of confidence matrizes.
-
-        Parameters
-        ----------
-        conf_list: list of :class:`~numpy.ndarray` (optional)
-            List of 2D confidence matrizes (one for each projection) which define trust regions.
-            If not given this uses the confidence matrizes of the phase maps.
-
-        Returns
-        -------
-            None
-
-        '''
-        self._log.debug('Calling set_Se_inv_diag_with_conf')
-        if conf_list is None:  # if no confidence matrizes are given, extract from the phase maps!
-            conf_list = [phase_map.confidence for phase_map in self.phase_maps]
-        cov_list = [sparse.diags(c.flatten().astype(np.float32), 0) for c in conf_list]
-        self.set_Se_inv_block_diag(cov_list)
-
-    def set_3d_mask(self, mask_list=None):
-        '''Set the 3D mask from a list of 2D masks.
-
-        Parameters
-        ----------
-        mask_list: list of :class:`~numpy.ndarray` (optional)
-            List of 2D masks, which represent the projections of the 3D mask. If not given this
-            uses the mask matrizes of the phase maps. If just one phase map is present, the
-            according mask is simply expanded to 3D and used directly.
-
-        Returns
-        -------
-            None
-
-        '''
-        self._log.debug('Calling set_3d_mask')
-        if mask_list is None:  # if no masks are given, extract from phase maps:
-            mask_list = [phase_map.mask for phase_map in self.phase_maps]
-        if len(mask_list) == 1:  # just one phase_map --> 3D mask equals 2D mask
-            self.mask = np.expand_dims(mask_list[0], axis=0)  # z-dim is set to 1!
-        else:  # 3D mask has to be constructed from 2D masks:
-            mask_3d_inv = np.zeros(self.dim)
-            for i, projector in enumerate(self.projectors):
-                mask_2d_inv = np.logical_not(self.phase_maps[i].mask.reshape(-1))  # inv. 2D mask
-                # Add extrusion of inv. 2D mask:
-                mask_3d_inv += projector.weight.T.dot(mask_2d_inv).reshape(self.dim)
-            self.mask = np.where(mask_3d_inv == 0, True, False)
-
-    def display_mask(self, ar_dens=1):
-        '''If it exists, display the 3D mask of the magnetization distribution.
-
-        Parameters
-        ----------
-        ar_dens: int (optional)
-            Number defining the cell density which is plotted. A higher ar_dens number skips more
-            arrows (a number of 2 plots every second arrow). Default is 1.
-
-        Returns
-        -------
-            None
-
-        '''
-        self._log.debug('Calling display_mask')
-        if self.mask is not None:
-            from mayavi import mlab
-            zz, yy, xx = np.indices(self.dim)
-            ad = ar_dens
-            zz = zz[::ad, ::ad, ::ad].flatten()
-            yy = yy[::ad, ::ad, ::ad].flatten()
-            xx = xx[::ad, ::ad, ::ad].flatten()
-            mask_vec = self.mask[::ad, ::ad, ::ad].flatten().astype(dtype=np.int)
-            mlab.figure(size=(750, 700))
-            plot = plot = mlab.points3d(xx, yy, zz, mask_vec, opacity=0.5, mode='cube',
-                                        scale_factor=ar_dens)
-            mlab.outline(plot)
-            mlab.axes(plot)
-            return plot
-
-    def display_phase(self, mag_data=None, title='Phase Map',
-                      cmap='RdBu', limit=None, norm=None):
-        '''Display all phasemaps saved in the :class:`~.DataSet` as a colormesh.
-
-        Parameters
-        ----------
-        mag_data : :class:`~.MagData`, optional
-            Magnetic distribution to which the projectors of the dataset should be applied. If not
-            given, the phase_maps in the dataset are used.
-        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_phase')
-        if mag_data is not None:
-            phase_maps = self.create_phase_maps(mag_data)
-        else:
-            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(phase_maps)]
-        plt.show()
-
-    def display_combined(self, mag_data=None, title='Combined Plot', cmap='RdBu', limit=None,
-                         norm=None, gain='auto', interpolation='none', grad_encode='bright'):
-        '''Display all phasemaps and the resulting color coded holography images.
-
-        Parameters
-        ----------
-        mag_data : :class:`~.MagData`, optional
-            Magnetic distribution to which the projectors of the dataset should be applied. If not
-            given, the phase_maps in the dataset are used.
-        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.
-        gain : 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
-        -------
-        None
-
-        '''
-        self._log.debug('Calling display_combined')
-        if mag_data is not None:
-            phase_maps = self.create_phase_maps(mag_data)
-        else:
-            phase_maps = self.phase_maps
-        for (i, phase_map) in enumerate(phase_maps):
-            phase_map.display_combined('{} ({})'.format(title, self.projectors[i].get_info()),
-                                       cmap, limit, norm, gain, interpolation, grad_encode)
-        plt.show()

pyramid/diagnostics.py

-# -*- coding: utf-8 -*-
-# Copyright 2014 by Forschungszentrum Juelich GmbH
-# Author: J. Caron
-#
-"""This module provides the :class:`~.Diagnostics` class for the calculation of diagnostics of a
-specified costfunction for a fixed magnetization distribution."""
-
-
-import numpy as np
-import matplotlib.pyplot as plt
-
-import jutil
-from pyramid import fft
-from pyramid.magdata import MagData
-from pyramid.phasemap import PhaseMap
-
-import logging
-
-
-__all__ = ['Diagnostics']
-
-
-class Diagnostics(object):
-
-    '''Class for calculating diagnostic properties of a specified costfunction.
-
-    For the calculation of diagnostic properties, a costfunction and a magnetization distribution
-    are specified at construction. With the :func:`~.set_position`, a position in 3D space can be
-    set at which all properties will be calculated. Properties are saved via boolean flags and
-    thus, calculation is only done if the position has changed in between. The standard deviation
-    and the measurement contribution require the execution of a conjugate gradient solver and can
-    take a while for larger problems.
-
-    Attributes
-    ----------
-    x_rec: :class:`~numpy.ndarray`
-        Vectorized magnetization distribution at which the costfunction is evaluated.
-    cost: :class:`~.pyramid.costfunction.Costfunction`
-        Costfunction for which the diagnostics are calculated.
-    max_iter: int, optional
-        Maximum number of iterations. Default is 1000.
-    fwd_model: :class:`~pyramid.forwardmodel.ForwardModel`
-        Forward model used in the costfunction.
-    Se_inv : :class:`~numpy.ndarray` (N=2), optional
-        Inverted covariance matrix of the measurement errors. The matrix has size `NxN` with N
-        being the length of the targetvector y (vectorized phase map information).
-    dim: tuple (N=3)
-        Dimensions of the 3D magnetic distribution.
-    row_idx: int
-        Row index of the system matrix corresponding to the current position in 3D space.
-    cov_row: :class:`~numpy.ndarray`
-        Row of the covariance matrix (``S_a^-1+F'(x_f)^T S_e^-1 F'(x_f)``) which is needed for the
-        calculation of the gain and averaging kernel matrizes and which ideally contains the
-        variance at position `row_idx` for the current component and position in 3D.
-    std: float
-        Standard deviation of the chosen component at the current position (calculated when
-        needed).
-    gain_row: :class:`~numpy.ndarray`
-        Row of the gain matrix, which maps differences of phase measurements onto differences in
-        the retrieval result of the magnetization distribution(calculated when needed).
-    avrg_kern_row: :class:`~numpy.ndarray`
-        Row of the averaging kernel matrix (which is ideally the identity matrix), which describes
-        the smoothing introduced by the regularization (calculated when needed).
-    measure_contribution: :class:`~numpy.ndarray`
-
-    Notes
-    -----
-        Some properties depend on others, which may require recalculations of these prior
-        properties if necessary. The dependencies are ('-->' == 'requires'):
-        avrg_kern_row --> gain_row --> std --> m_inv_row
-        measure_contribution is independant
-
-    '''
-
-    _log = logging.getLogger(__name__+'.Diagnostics')
-
-    @property
-    def cov_row(self):
-        if not self._updated_cov_row:
-            e_i = fft.zeros(self.cost.n, dtype=fft.FLOAT)
-            e_i[self.row_idx] = 1
-            row = 2 * jutil.cg.conj_grad_solve(self._A, e_i, P=self._P, max_iter=self.max_iter)
-            self._std_row = row
-            self._updated_cov_row = True
-        return self._std_row
-
-    @property
-    def std(self):
-        return np.sqrt(self.cov_row[self.row_idx])
-
-    @property
-    def gain_row(self):
-        if not self._updated_gain_row:
-            self._gain_row = self.Se_inv.dot(self.fwd_model.jac_dot(self.x_rec, self.cov_row))
-            self._updated_gain_row = True
-        return self._gain_row
-
-    @property
-    def avrg_kern_row(self):
-        if not self._updated_avrg_kern_row:
-            self._avrg_kern_row = self.fwd_model.jac_T_dot(self.x_rec, self.gain_row)
-            self._updated_avrg_kern_row = True
-        return self._avrg_kern_row
-
-    @property
-    def measure_contribution(self):
-        if not self._updated_measure_contribution:
-            cache = self.fwd_model.jac_dot(self.x_rec, fft.ones(self.cost.n, fft.FLOAT))
-            cache = self.fwd_model.jac_T_dot(self.x_rec, self.Se_inv.dot(cache))
-            mc = 2 * jutil.cg.conj_grad_solve(self._A, cache, P=self._P, max_iter=self.max_iter)
-            self._measure_contribution = mc
-            self._updated_measure_contribution = True
-        return self._measure_contribution
-
-    @property
-    def pos(self):
-        return self._pos
-
-    @pos.setter
-    def pos(self, pos):
-        c, z, y, x = pos
-        assert self.mask[z, y, x], 'Position is outside of the provided mask!'
-        mask_vec = self.mask.flatten()
-        idx_3d = z*self.dim[1]*self.dim[2] + y*self.dim[2] + x
-        row_idx = c*np.prod(mask_vec.sum()) + mask_vec[:idx_3d].sum()
-        if row_idx != self.row_idx:
-            self._pos = pos
-            self.row_idx = row_idx
-            self._updated_cov_row = False
-            self._updated_gain_row = False
-            self._updated_avrg_kern_row = False
-            self._updated_measure_contribution = False
-
-    def __init__(self, x_rec, cost, max_iter=1000):
-        self._log.debug('Calling __init__')
-        self.x_rec = x_rec
-        self.cost = cost
-        self.max_iter = max_iter
-        self.fwd_model = cost.fwd_model
-        self.Se_inv = cost.Se_inv
-        self.dim = cost.data_set.dim
-        self.mask = cost.data_set.mask
-        self.row_idx = None
-        self.pos = (0,) + tuple(np.array(np.where(self.mask))[:, 0])  # first True mask entry
-        self._updated_cov_row = False
-        self._updated_gain_row = False
-        self._updated_avrg_kern_row = False
-        self._updated_measure_contribution = False
-        self._A = jutil.operator.CostFunctionOperator(self.cost, self.x_rec)
-        self._P = jutil.preconditioner.CostFunctionPreconditioner(self.cost, self.x_rec)
-        self._log.debug('Creating '+str(self))
-
-    def get_avg_kern_row(self, pos=None):
-        '''Get the averaging kernel matrix row represented as a 3D magnetization distribution.
-
-        Parameters
-        ----------
-        pos: tuple (N=4)
-            Position in 3D plus component `(c, z, y, x)`
-
-        Returns
-        -------
-        mag_data_avg_kern: :class:`~pyramid.magdata.MagData`
-            Averaging kernel matrix row represented as a 3D magnetization distribution
-
-        '''
-        self._log.debug('Calling get_avg_kern_row')
-        if pos is not None:
-            self.pos = pos
-        mag_data_avg_kern = MagData(self.cost.data_set.a, fft.zeros((3,)+self.dim))
-        mag_data_avg_kern.set_vector(self.avrg_kern_row, mask=self.mask)
-        return mag_data_avg_kern
-
-    def calculate_averaging(self, pos=None):
-        '''Calculate and plot the averaging pixel number at a specified position for x, y or z.
-
-        Parameters
-        ----------
-        pos: tuple (N=4)
-            Position in 3D plus component `(c, z, y, x)`
-
-        Returns
-        -------
-        px_avrg: float
-            The number of pixels over which is approximately averaged. The inverse is the FWHM.
-        (x, y, z): tuple of floats
-            The magnitude of the averaging kernel summed along two axes (the remaining are x, y, z,
-            respectively).
-
-        Notes
-        -----
-        Uses the :func:`~.get_avg_kern_row' function
-
-        '''
-        self._log.debug('Calling calculate_averaging')
-        mag_data_avg_kern = self.get_avg_kern_row(pos)
-        mag_x, mag_y, mag_z = mag_data_avg_kern.magnitude
-        x = mag_x.sum(axis=(0, 1))
-        y = mag_y.sum(axis=(0, 2))
-        z = mag_z.sum(axis=(1, 2))
-        # Plot helpful stuff:
-        plt.figure()
-        plt.axhline(y=0, ls='-', color='k')
-        plt.axhline(y=1, ls='-', color='k')
-        plt.plot(x, label='x', color='r', marker='o')
-        plt.plot(y, label='y', color='g', marker='o')
-        plt.plot(z, label='z', color='b', marker='o')
-        c = self.pos[0]
-        data = [x, y, z][c]
-        col = ['r', 'g', 'b'][c]
-        i_m = np.argmax(data)  # Index of the maximum
-        plt.axhline(y=data[i_m], ls='-', color=col)
-        plt.axhline(y=data[i_m]/2, ls='--', color=col)
-        # Left side:
-        for i in np.arange(i_m-1, -1, -1):
-            if data[i] < data[i_m]/2:
-                l = (data[i_m]/2-data[i])/(data[i+1]-data[i]) + i
-                break
-        # Right side:
-        for i in np.arange(i_m+1, data.size):
-            if data[i] < data[i_m]/2:
-                r = (data[i_m]/2-data[i-1])/(data[i]-data[i-1]) + i-1
-                break
-        # Calculate FWHM:
-        fwhm = r - l
-        px_avrg = 1 / fwhm
-        plt.vlines(x=[l, r], ymin=0, ymax=data[i_m]/2, linestyles=':', color=col)
-        # Add legend:
-        plt.legend()
-        return px_avrg, (x, y, z)
-
-    def get_gain_row_maps(self, pos=None):
-        '''Get the gain matrix row represented as a list of 2D (inverse) phase maps.
-
-        Parameters
-        ----------
-        pos: tuple (N=4)
-            Position in 3D plus component `(c, z, y, x)`
-
-        Returns
-        -------
-        gain_map_list: list of :class:`~pyramid.phasemap.PhaseMap`
-            Gain matrix row represented as a list of 2D phase maps
-
-        Notes
-        -----
-        Note that the produced gain maps define the magnetization change at the current position
-        in 3d per phase change at the position of the . Take this into account when plotting the
-        maps (1/rad instead of rad).
-
-        '''
-        self._log.debug('Calling get_gain_row_maps')
-        if pos is not None:
-            self.pos = pos
-        hp = self.cost.data_set.hook_points
-        gain_map_list = []
-        for i, projector in enumerate(self.cost.data_set.projectors):
-            gain = self.gain_row[hp[i]:hp[i+1]].reshape(projector.dim_uv)
-            gain_map_list.append(PhaseMap(self.cost.data_set.a, gain))
-        return gain_map_list

pyramid/fft.py

-# -*- coding: utf-8 -*-
-# Copyright 2014 by Forschungszentrum Juelich GmbH
-# Author: J. Caron
-#
-"""Custom FFT module with numpy and FFTW support.
-
-This module provides custom methods for FFTs including inverse, adjoint and real variants. The
-FFTW library is supported and is used as a default if the import succeeds. Otherwise the numpy.fft
-pack will be used. FFTW objects are saved in a cache after creation which speeds up further similar
-FFT operations.
-
-"""
-
-
-import numpy as np
-import cPickle as pickle
-import os
-
-from pyramid.config import NTHREADS
-
-import logging
-
-
-# pyFFTW depends on this
-try:
-    from collections import Counter  # analysis:ignore
-except ImportError:
-    import collections_python27
-    import collections
-    collections.Counter = collections_python27.Counter
-
-try:
-    import pyfftw
-    BACKEND = 'fftw'
-except ImportError:
-    pyfftw = None
-    BACKEND = 'numpy'
-    print("pyFFTW module not found. Using numpy implementation.")
-
-
-__all__ = ['PLANS', 'FLOAT', 'COMPLEX', 'dump_wisdom', 'load_wisdom',  # analysis:ignore
-           'zeros', 'empty', 'configure_backend',
-           'fftn', 'ifftn', 'rfftn', 'irfftn', 'rfftn_adj', 'irfftn_adj']
-_log = logging.getLogger(__name__)
-
-
-class FFTWCache(object):
-
-    '''Class for adding FFTW Plans and on-demand lookups.
-
-    This class is instantiated in this module to store FFTW plans and for the lookup of the former.
-
-    Attributes
-    ----------
-    cache: dict
-        Cache for storing the FFTW plans.
-
-    Notes
-    -----
-    This class is used internally and is not normally not intended to be used directly by the user.
-
-    '''
-
-    _log = logging.getLogger(__name__+'.FFTWCache')
-
-    def __init__(self):
-        self._log.debug('Calling __init__')
-        self.cache = dict()
-        self._log.debug('Created '+str(self))
-
-    def add_fftw(self, fft_type, fftw_obj, s, axes, nthreads):
-        self._log.debug('Calling add_fftw')
-        in_arr = fftw_obj.get_input_array()
-        key = (fft_type, in_arr.shape, in_arr.dtype, nthreads)
-        self.cache[key] = fftw_obj
-
-    def lookup_fftw(self, fft_type, in_arr, s, axes, nthreads):
-        self._log.debug('Calling lookup_fftw')
-        key = (fft_type, in_arr.shape, in_arr.dtype, nthreads)
-        return self.cache.get(key, None)
-
-    def clear_cache(self):
-        self._log.debug('Calling clear_cache')
-        self.cache = dict()
-
-
-PLANS = FFTWCache()
-FLOAT = np.float32      # One convenient place to
-COMPLEX = np.complex64  # change from 32 to 64 bit
-
-
-# Numpy functions:
-
-def _fftn_numpy(a, s=None, axes=None):
-    return np.fft.fftn(a, s, axes)
-
-
-def _ifftn_numpy(a, s=None, axes=None):
-    return np.fft.ifftn(a, s, axes)
-
-
-def _rfftn_numpy(a, s=None, axes=None):
-    return np.fft.rfftn(a, s, axes)
-
-
-def _irfftn_numpy(a, s=None, axes=None):
-    return np.fft.irfftn(a, s, axes)
-
-
-def _rfftn_adj_numpy(a):
-    n = 2 * (a.shape[-1] - 1)
-    out_shape = a.shape[:-1] + (n,)
-    out_arr = zeros(out_shape, dtype=a.dtype)
-    out_arr[:, :n] = a
-    return _ifftn_numpy(out_arr).real * np.prod(out_shape)
-
-
-def _irfftn_adj_numpy(a):
-    n = a.shape[-1] // 2 + 1
-    out_arr = _fftn_numpy(a, axes=(-1,)) / a.shape[-1]
-    if a.shape[-1] % 2 == 0:  # even
-        out_arr[:, 1:n - 1] += np.conj(out_arr[:, :n-1:-1])
-    else:  # odd
-        out_arr[:, 1:n] += np.conj(out_arr[:, :n-1:-1])
-    axes = tuple(range(len(out_arr.shape[:-1])))
-    return _fftn_numpy(out_arr[:, :n], axes=axes) / np.prod(out_arr.shape[:-1])
-
-
-# FFTW functions:
-
-def _fftn_fftw(a, s=None, axes=None):
-    if a.dtype not in (FLOAT, COMPLEX):
-        raise TypeError('Wrong input type!')
-    fftw = PLANS.lookup_fftw('fftn', a, s, axes, NTHREADS)
-    if fftw is None:
-        fftw = pyfftw.builders.fftn(a, s, axes, threads=NTHREADS)
-        PLANS.add_fftw('fftn', fftw, s, axes, NTHREADS)
-    return fftw(a).copy()
-
-
-def _ifftn_fftw(a, s=None, axes=None):
-    if a.dtype not in (FLOAT, COMPLEX):
-        raise TypeError('Wrong input type!')
-    fftw = PLANS.lookup_fftw('ifftn', a, s, axes, NTHREADS)
-    if fftw is None:
-        fftw = pyfftw.builders.ifftn(a, s, axes, threads=NTHREADS)
-        PLANS.add_fftw('ifftn', fftw, s, axes, NTHREADS)
-    return fftw(a).copy()
-
-
-def _rfftn_fftw(a, s=None, axes=None):
-    if a.dtype != FLOAT:
-        raise TypeError('Wrong input type!')
-    fftw = PLANS.lookup_fftw('rfftn', a, s, axes, NTHREADS)
-    if fftw is None:
-        fftw = pyfftw.builders.rfftn(a, s, axes, threads=NTHREADS)
-        PLANS.add_fftw('rfftn', fftw, s, axes, NTHREADS)
-    return fftw(a).copy()
-
-
-def _irfftn_fftw(a, s=None, axes=None):
-    if a.dtype != COMPLEX:
-        raise TypeError('Wrong input type!')
-    fftw = PLANS.lookup_fftw('irfftn', a, s, axes, NTHREADS)
-    if fftw is None:
-        fftw = pyfftw.builders.irfftn(a, s, axes, threads=NTHREADS)
-        PLANS.add_fftw('irfftn', fftw, s, axes, NTHREADS)
-    return fftw(a).copy()
-
-
-def _rfftn_adj_fftw(a):
-    # Careful just works for even a (which is guaranteed by the kernel!)
-    n = 2 * (a.shape[-1] - 1)
-    out_shape = a.shape[:-1] + (n,)
-    out_arr = zeros(out_shape, dtype=a.dtype)
-    out_arr[:, :a.shape[-1]] = a
-    return _ifftn_fftw(out_arr).real * np.prod(out_shape)
-
-
-def _irfftn_adj_fftw(a):
-    out_arr = _fftn_fftw(a, axes=(-1,)) / a.shape[-1]  # FFT of last axis
-    n = a.shape[-1] // 2 + 1
-    if a.shape[-1] % 2 == 0:  # even
-        out_arr[:, 1:n-1] += np.conj(out_arr[:, :n-1:-1])
-    else:  # odd
-        out_arr[:, 1:n] += np.conj(out_arr[:, :n-1:-1])
-    axes = tuple(range(len(out_arr.shape[:-1])))
-    return _fftn_fftw(out_arr[:, :n], axes=axes) / np.prod(out_arr.shape[:-1])
-
-
-# These wisdom functions do nothing if pyFFTW is not available:
-
-def dump_wisdom(fname):
-    '''Wrapper function for the pyfftw.export_wisdom(), which uses a pickle dump.
-
-    Parameters
-    ----------
-    fname: string
-        Name of the file in which the wisdom is saved.
-
-    Returns
-    -------
-    None
-
-    '''
-    _log.debug('Calling dump_wisdom')
-    if pyfftw is not None:
-        with open(fname, 'wb') as fp:
-            pickle.dump(pyfftw.export_wisdom(), fp, pickle.HIGHEST_PROTOCOL)
-
-
-def load_wisdom(fname):
-    '''Wrapper function for the pyfftw.import_wisdom(), which uses a pickle to load a file.
-
-    Parameters
-    ----------
-    fname: string
-        Name of the file from which the wisdom is loaded.
-
-    Returns
-    -------
-    None
-
-    '''
-    _log.debug('Calling load_wisdom')
-    if pyfftw is not None:
-        if not os.path.exists(fname):
-            print("Warning: Wisdom file does not exist. First time use?")
-        else:
-            with open(fname, 'rb') as fp:
-                pyfftw.import_wisdom(pickle.load(fp))
-
-
-# Array setups:
-def empty(shape, dtype=FLOAT):
-    '''Return a new array of given shape and type without initializing entries.
-
-    Parameters
-    ----------
-    shape: int or tuple of int
-        Shape of the array.
-    dtype: data-type, optional
-        Desired output data-type.
-
-    Returns
-    -------
-    out: :class:`~numpy.ndarray`
-        The created array.
-
-    '''
-    _log.debug('Calling empty')
-    result = np.empty(shape, dtype)
-    if pyfftw is not None:
-        result = pyfftw.n_byte_align(result, pyfftw.simd_alignment)
-    return result
-
-
-def zeros(shape, dtype=FLOAT):
-    '''Return a new array of given shape and type, filled with zeros.
-
-    Parameters
-    ----------
-    shape: int or tuple of int
-        Shape of the array.
-    dtype: data-type, optional
-        Desired output data-type.
-
-    Returns
-    -------
-    out: :class:`~numpy.ndarray`
-        The created array.
-
-    '''
-    _log.debug('Calling zeros')
-    result = np.zeros(shape, dtype)
-    if pyfftw is not None:
-        result = pyfftw.n_byte_align(result, pyfftw.simd_alignment)
-    return result
-
-
-def ones(shape, dtype=FLOAT):
-    '''Return a new array of given shape and type, filled with ones.
-
-    Parameters
-    ----------
-    shape: int or tuple of int
-        Shape of the array.
-    dtype: data-type, optional
-        Desired output data-type.
-
-    Returns
-    -------
-    out: :class:`~numpy.ndarray`
-        The created array.
-
-    '''
-    _log.debug('Calling ones')
-    result = np.ones(shape, dtype)
-    if pyfftw is not None:
-        result = pyfftw.n_byte_align(result, pyfftw.simd_alignment)
-    return result
-
-
-# Configure backend:
-def configure_backend(backend):
-    '''Change FFT backend.
-
-    Parameters
-    ----------
-    backend: string
-        Backend to use. Supported values are "numpy" and "fftw".
-
-    Returns
-    -------
-    None
-
-    '''
-    _log.debug('Calling configure_backend')
-    global fftn
-    global ifftn
-    global rfftn
-    global irfftn
-    global rfftn_adj
-    global irfftn_adj
-    global BACKEND
-    if backend == "numpy":
-        fftn = _fftn_numpy
-        ifftn = _ifftn_numpy
-        rfftn = _rfftn_numpy
-        irfftn = _irfftn_numpy
-        rfftn_adj = _rfftn_adj_numpy
-        irfftn_adj = _irfftn_adj_numpy
-        BACKEND = 'numpy'
-    elif backend == "fftw":
-        if pyfftw is not None:
-            fftn = _fftn_fftw
-            ifftn = _ifftn_fftw
-            rfftn = _rfftn_fftw
-            irfftn = _irfftn_fftw
-            rfftn_adj = _rfftn_adj_fftw
-            irfftn_adj = _irfftn_adj_fftw
-            BACKEND = 'pyfftw'
-        else:
-            print("Error: FFTW requested but not available")
-
-
-# On import:
-if pyfftw is not None:
-    configure_backend("fftw")
-else:
-    configure_backend("numpy")