Skip to content
Snippets Groups Projects

Compare revisions

Changes are shown as if the source revision was being merged into the target revision. Learn more about comparing revisions.

Source

Select target project
No results found
Select Git revision
  • fix-squeeze
  • master
  • pyramid-master
  • 0.0.0
  • 0.0.1
  • 0.1.0
  • 0.2.0
  • 0.2.1
  • 0.3.0
  • 0.3.1
  • 0.3.2
  • 0.3.2-dev
  • 0.3.3
  • 0.3.4
  • 0.3.5
  • testtag
16 results

Target

Select target project
  • empyre/empyre
  • weber/empyre
  • wessels/empyre
  • bryan/empyre
4 results
Select Git revision
  • add_field_tests
  • colorclipping
  • master
  • pyramid-master
  • support-ragged
  • support36
  • 0.0.0
  • 0.0.1
  • 0.1.0
  • 0.2.0
  • 0.2.1
11 results
Show changes
Expand all Hide whitespace changes
Inline Side-by-side
Showing
with 0 additions and 5881 deletions
pyramid/analytic.py deleted 100644 → 0
View file @ d9a1f7d3
# -*- 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 logging
import numpy as np
from numpy import pi
from pyramid.phasemap import PhaseMap
__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 deleted 100644 → 0
View file @ d9a1f7d3
# -*- 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 logging
from numbers import Number
import matplotlib.pyplot as plt
import numpy as np
from PIL import Image
from matplotlib import colors
__all__ = ['DirectionalColormap', 'TransparentColormap']
class DirectionalColormap(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
:func:`.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 = colors.LinearSegmentedColormap('my_colormap', CDICT, 256)
HOLO_CMAP_INV = 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().__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 np.asarray(255.999 * (rgb_norm + rgb_invs), dtype=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(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=None):
if alpha_range is None:
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().__init__('transparent_colormap', cdict, N=256)
self._log.debug('Created ' + str(self))
pyramid/costfunction.py deleted 100644 → 0
View file @ d9a1f7d3
# -*- 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 logging
import numpy as np
from pyramid.regularisator import NoneRegularisator
__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 deleted 100644 → 0
View file @ d9a1f7d3
# -*- 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 logging
from numbers import Number
from pyramid.kernel import Kernel
from pyramid.phasemap import PhaseMap
from pyramid.phasemapper import PhaseMapperRDFC
from pyramid.projector import Projector
import matplotlib.pyplot as plt
import numpy as np
from scipy import sparse
__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.
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.
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).
phase_maps:
A list of all stored :class:`~.PhaseMap` objects.
projectors: list of :class:`~.Projector`
A list of all stored :class:`~.Projector` objects.
"""
_log = logging.getLogger(__name__ + '.DataSet')
@property
def a(self):
"""The grid spacing in nm."""
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):
"""A boolean mask which defines the magnetized volume in 3D."""
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):
"""Size of the image space."""
return np.sum([len(p.phase_vec) for p in self.phase_maps])
@property
def n(self):
"""Size of the input space."""
return 3 * np.sum(self.mask)
@property
def count(self):
"""Number of phase maps and projectors in the dataset."""
return len(self.projectors)
@property
def phase_vec(self):
"""The concatenaded, vectorized phase of all ;class:`~.PhaseMap` objects."""
return np.concatenate([p.phase_vec for p in self.phase_maps])
@property
def hook_points(self):
"""Hook points which determine the start of values of a phase map in the `phase_vec`."""
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):
"""List of phase mappers, created on demand with the projectors in mind."""
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:`~.VectorData` object.
Parameters
----------
mag_data : :class:`~.VectorData`
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 = 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:`~.VectorData`, 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=cmap, limit=limit, norm=norm)
for (i, phase_map) in enumerate(phase_maps)]
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:`~.VectorData`, 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 deleted 100644 → 0
View file @ d9a1f7d3
# -*- 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 logging
from pyramid import fft
from pyramid.fielddata import VectorData
from pyramid.phasemap import PhaseMap
import matplotlib.pyplot as plt
import numpy as np
import jutil
__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.
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):
"""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."""
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 = np.asarray(row)
self._updated_cov_row = True
return self._std_row
@property
def std(self):
"""Standard deviation of the chosen component at the current position (calculated when
needed)."""
return np.sqrt(self.cov_row[self.row_idx])
@property
def gain_row(self):
"""Row of the gain matrix, which maps differences of phase measurements onto differences in
the retrieval result of the magnetization distribution(calculated when needed)."""
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):
"""Row of the averaging kernel matrix (which is ideally the identity matrix), which
describes the smoothing introduced by the regularization (calculated when needed)."""
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):
"""The sum over an averaging kernel matrix row, which is an indicator for wheter a point of
the solution is determined by the measurement (close to `1`) or by a priori information
(close to `0`)."""
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):
"""The current solution position, which specifies the 3D-point (and the component) of the
magnetization, for which diagnostics are calculated."""
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.fielddata.VectorData`
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 = VectorData(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 field 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.field
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] # x, y or z component
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:
l = i_m
for i in np.arange(i_m - 1, -1, -1):
if data[i] < data[i_m] / 2:
# Linear interpolation between i and i + 1 to find left fractional index position:
l = (data[i_m] / 2 - data[i]) / (data[i + 1] - data[i]) + i
break
# Right side:
r = i_m
for i in np.arange(i_m + 1, data.size):
if data[i] < data[i_m] / 2:
# Linear interpolation between i and i - 1 to find right fractional index position:
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 deleted 100644 → 0
View file @ d9a1f7d3
This diff is collapsed.
pyramid/fieldconverter.py deleted 100644 → 0
View file @ d9a1f7d3
This diff is collapsed.
pyramid/fielddata.py deleted 100644 → 0
View file @ d9a1f7d3
This diff is collapsed.
pyramid/forwardmodel.py deleted 100644 → 0
View file @ d9a1f7d3
This diff is collapsed.
pyramid/gui/__init__.py deleted 100644 → 0
View file @ d9a1f7d3
# -*- coding: utf-8 -*-
"""Package containing Pyramid GUI functions."""
from .phasemap_creator import phasemap_creator
from .mag_slicer import mag_slicer
__all__ = ['phasemap_creator', 'mag_slicer']
pyramid/gui/mag_slicer.py deleted 100644 → 0
View file @ d9a1f7d3
This diff is collapsed.
pyramid/gui/mag_slicer.ui deleted 100644 → 0
View file @ d9a1f7d3
This diff is collapsed.
pyramid/gui/phasemap_creator.py deleted 100644 → 0
View file @ d9a1f7d3
This diff is collapsed.
pyramid/gui/phasemap_creator.ui deleted 100644 → 0
View file @ d9a1f7d3
This diff is collapsed.
pyramid/kernel.py deleted 100644 → 0
View file @ d9a1f7d3
This diff is collapsed.
pyramid/magcreator/__init__.py deleted 100644 → 0
View file @ d9a1f7d3
# -*- coding: utf-8 -*-
from .magcreator import *
from . import shapes
from . import examples
import logging
_log = logging.getLogger(__name__)
del logging
__all__ = ['shapes', 'examples']
__all__.extend(magcreator.__all__)
pyramid/magcreator/examples.py deleted 100644 → 0
View file @ d9a1f7d3
This diff is collapsed.
pyramid/magcreator/magcreator.py deleted 100644 → 0
View file @ d9a1f7d3
This diff is collapsed.
pyramid/magcreator/shapes.py deleted 100644 → 0
View file @ d9a1f7d3
This diff is collapsed.
pyramid/phasemap.py deleted 100644 → 0
View file @ d9a1f7d3
This diff is collapsed.