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

Target

Select target project
  • empyre/empyre
  • weber/empyre
  • wessels/empyre
  • bryan/empyre
4 results
Show changes
Expand all Hide whitespace changes
Inline Side-by-side
Showing
with 0 additions and 6193 deletions
pyramid/analytic.py deleted 100644 → 0
View file @ 68453aaa
# -*- coding: utf-8 -*-
# Copyright 2016 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
-------
phasemap : :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
-------
phasemap : :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
-------
phasemap : :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
-------
phasemap : :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/colors.py deleted 100644 → 0
View file @ 68453aaa
# -*- coding: utf-8 -*-
# Copyright 2014 by Forschungszentrum Juelich GmbH
# Author: J. Caron
#
"""This module provides a custom :class:`~.HLSTriadicColormap` colormap class which has a few
additional functions and can encode three-dimensional directions."""
import logging
import matplotlib.pyplot as plt
import numpy as np
from PIL import Image
from matplotlib import colors
import colorsys
__all__ = ['TransparentColormap', 'HLSTriadicColormap', 'HLSTetradicColormap', 'hls_triadic_cmap',
'hls_tetradic_cmap', 'plot_colorwheel', 'rgb_from_hls',
'hls_from_vector', 'rgb_from_vector']
_log = logging.getLogger(__name__)
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):
self._log.debug('Calling __init__')
if alpha_range is None:
alpha_range = [0., 1.]
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', cdict, N=256)
self._log.debug('Created ' + str(self))
class HLSTriadicColormap(colors.ListedColormap):
"""Colormap subclass for encoding directions with colors.
This class is a subclass of the :class:`~matplotlib.pyplot.colors.ListedColormap`
class. The class follows the HSL ('hue', 'saturation', 'lightness') 'Double Hexcone' Model
with the saturation always set to 1 (moving on the surface of the color
cylinder) with a luminance of 0.5 (full color). The three prime colors (`rgb`) are spaced
equidistant with 120° space in between, according to a triadic arrangement.
"""
_log = logging.getLogger(__name__ + '.HLSTriadicColormap')
def __init__(self):
self._log.debug('Calling __init__')
h = np.linspace(0, 1, 256)
l = 0.5 * np.ones_like(h)
s = np.ones_like(h)
rgb = rgb_from_hls(h, l, s)
colors = [(ri / 255, gi / 255, bi / 255) for ri, gi, bi in rgb]
super().__init__(colors, 'hlscm', N=256)
self._log.debug('Created ' + str(self))
class HLSTetradicColormap(colors.LinearSegmentedColormap):
"""Colormap subclass for encoding directions with colors.
This class is a subclass of the :class:`~matplotlib.pyplot.colors.LinearSegmentedColormap`
class. The class follows the HSL ('hue', 'saturation', 'lightness') 'Double
Hexcone' Model with the saturation always set to 1 (moving on the surface of the color
cylinder) with a luminance of 0.5 (full color). The colors follow a tetradic arrangement with
four colors (red, green, blue and yellow) arranged with 90° spacing in between.
"""
_log = logging.getLogger(__name__ + '.HLSTetradicColormap')
CDICT = {'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, 0.0, 0.0),
(0.25, 1.0, 1.0),
(0.50, 0.0, 0.0),
(0.75, 0.0, 0.0),
(1.00, 0.0, 0.0)]}
def __init__(self):
self._log.debug('Calling __init__')
super().__init__('holo', self.CDICT, N=256)
self._log.debug('Created ' + str(self))
def plot_colorwheel(mode='triadic'):
"""Display a color wheel to illustrate the color coding of vector gradient directions.
Parameters
----------
mode : {'triadic', 'tetradic'}
Optional string for determining the hue scheme of the color wheel. The default is the
standard HLS scheme with the three primary colors (red, blue, green) spaced with 120°
between them (triadic arrangement). The other option is a tetradic scheme commonly
used in holography where red, green, blue and yellow are used with 90° spacing. In both
cases, the saturation decreases to the center of the circle.
Returns
-------
None
"""
_log.debug('Calling display_color_wheel')
# Construct the colorwheel:
yy, xx = np.indices((512, 512)) - 256
rr = np.hypot(xx, yy)
xx = np.where(rr <= 256, xx, 0)
yy = np.where(rr <= 256, yy, 0)
zz = np.where(rr <= 256, 0, -1)
h, l, s = hls_from_vector(xx, yy, zz)
rgb = rgb_from_hls(h, l, s, mode=mode)
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')
return
def rgb_from_hls(h, l, s, mode='triadic'):
"""Construct a rgb tuple from hue, luminance and saturation.
Parameters
----------
h : float, [0, 1)
Colorindex specifying the directional color according to the `mode` colormap scheme.
The colormap is periodic so that a value of 1 is equivalent to 0 again.
l : [0, 1] float
Luminance or lightness is the radiant power of the color. In the HLS model, it is defined
as the average between the largest and smallest color component. This definition puts
the primary (rgb) and secondary (ymc) colors into a plane halfway between black and
white for a luminance of 0.5. Directions pointing upwards will be encoded in white,
downwards in black.
s : [0, 1] float
The chroma scaled to fill the interval between 0 and 1. This means that smaller
amplitudes are encoded in grey, while large amplitudes will be in full color.
mode : {'triadic', 'tetradic'}
Optional string for determining the hue scheme of the color wheel. The default is the
standard HLS scheme with the three primary colors (red, blue, green) spaced with 120°
between them (triadic arrangement). The other option is a tetradic scheme commonly
used in holography where red, green, blue and yellow are used with 90° spacing. In both
cases, the saturation decreases to the center of the circle.
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)!
"""
_log.debug('Calling rgb_from_hls')
h = np.asarray(h)
l = np.asarray(l)
s = np.asarray(s)
hls_to_rgb = np.vectorize(colorsys.hls_to_rgb)
rgb_to_hls = np.vectorize(colorsys.rgb_to_hls)
if mode == 'triadic':
r, g, b = hls_to_rgb(h, l, s)
elif mode == 'tetradic':
rgb = hls_tetradic_cmap(h)
rh, gh, bh = rgb[..., 0], rgb[..., 1], rgb[..., 2]
rh = np.asarray(255 * rh).astype(np.uint8)
gh = np.asarray(255 * gh).astype(np.uint8)
bh = np.asarray(255 * bh).astype(np.uint8)
h_standard = rgb_to_hls(rh, gh, bh)[0]
r, g, b = hls_to_rgb(h_standard, l, s)
else:
raise ValueError('Mode not defined!')
return (255 * np.stack((r, g, b), axis=-1)).astype(np.uint8)
def hls_from_vector(x, y, z):
"""Construct a hls 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
-------
hls : tuple (N=3)
hls tuple with the encoded direction color.
Notes
-----
Also works with numpy arrays as input!
"""
_log.debug('Calling hls_from_vector')
x = np.asarray(x)
y = np.asarray(y)
z = np.asarray(z)
phi = np.asarray(np.arctan2(y, x))
phi[phi < 0] += 2 * np.pi
r = np.sqrt(x ** 2 + y ** 2 + z ** 2)
theta = np.arccos(z / (r + 1E-30))
h = phi / (2 * np.pi)
l = 1 - theta / np.pi
s = r / r.max()
return h, l, s
def rgb_from_vector(x, y, z, mode='triadic'):
"""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.
mode : {'triadic', 'tetradic'}
Optional string for determining the hue scheme of the color wheel. The default is the
standard HLS scheme with the three primary colors (red, blue, green) spaced with 120°
between them (triadic arrangement). The other option is a tetradic scheme commonly
used in holography where red, green, blue and yellow are used with 90° spacing. In both
cases, the saturation decreases to the center of the circle.
Returns
-------
rgb : tuple (N=3)
rgb tuple with the encoded direction color.
Notes
-----
Also works with numpy arrays as input!
"""
_log.debug('Calling rgb_from_vector')
return rgb_from_hls(*hls_from_vector(x, y, z), mode=mode)
transparent_cmap = TransparentColormap(0.2, 0.3, 0.2, [0.75, 0.])
hls_tetradic_cmap = HLSTetradicColormap()
hls_triadic_cmap = HLSTriadicColormap()
pyramid/costfunction.py deleted 100644 → 0
View file @ 68453aaa
# -*- 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
----------
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`.
regularisator : :class:`~.Regularisator`, optional
Regularisator class that's responsible for the regularisation term. If `None` or none is
given, no regularisation will be used.
y : :class:`~numpy.ndarray` (N=1)
Vector which lists all pixel values of all phase maps one after another.
m: int
Size of the image space.
n: int
Size of the input space.
Se_inv : :class:`~numpy.ndarray` (N=2), optional
Inverted covariance matrix of the measurement errors. The matrix has size `m x m` with m
being the length of the targetvector y.
"""
_log = logging.getLogger(__name__ + '.Costfunction')
def __init__(self, fwd_model, regularisator=None):
self._log.debug('Calling __init__')
self.fwd_model = fwd_model
if regularisator is None:
self.regularisator = NoneRegularisator()
else:
self.regularisator = regularisator
# Extract information from fwd_model:
self.y = self.fwd_model.y
self.n = self.fwd_model.n
self.m = self.fwd_model.m
self.Se_inv = self.fwd_model.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 @ 68453aaa
# -*- 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).
projectors: list of :class:`~.Projector`
A list of all stored :class:`~.Projector` objects.
phasemaps: 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.
count(self): int
Number of phase maps and projectors in the dataset.
hook_points(self): :class:`~numpy.ndarray` (N=1)
Hook points which determine the start of values of a phase map in the `phase_vec`.
The length is `count + 1`.
"""
_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.phasemaps])
@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.phasemaps])
@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, phasemap in enumerate(self.phasemaps):
result.append(result[i] + np.prod(phasemap.dim_uv))
return result
@property
def phasemappers(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.phasemaps = []
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, phasemap, projector):
"""Appends a data pair of phase map and projection infos to the data collection.`
Parameters
----------
phasemap: :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(phasemap, 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 phasemap.dim_uv == projector.dim_uv, 'Projection dimensions (dim_uv) must match!'
self.phasemaps.append(phasemap)
self.projectors.append(projector)
def create_phasemaps(self, magdata):
"""Create a list of phasemaps with the projectors in the dataset for a given
:class:`~.VectorData` object.
Parameters
----------
magdata : :class:`~.VectorData`
Magnetic distribution to which the projectors of the dataset should be applied.
Returns
-------
phasemaps : list of :class:`~.phasemap.PhaseMap`
A list of the phase maps resulting from the projections specified in the dataset.
"""
self._log.debug('Calling create_phasemaps')
phasemaps = []
for projector in self.projectors:
mag_proj = projector(magdata)
phasemap = self.phasemappers[projector.dim_uv](mag_proj)
phasemap.mask = mag_proj.get_mask()[0, ...]
phasemaps.append(phasemap)
return phasemaps
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.phasemaps), '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 = [phasemap.confidence for phasemap in self.phasemaps]
cov_list = [sparse.diags(c.ravel().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 = [phasemap.mask for phasemap in self.phasemaps]
if len(mask_list) == 1: # just one phasemap --> 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.phasemaps[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].ravel()
yy = yy[::ad, ::ad, ::ad].ravel()
xx = xx[::ad, ::ad, ::ad].ravel()
mask_vec = self.mask[::ad, ::ad, ::ad].ravel().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 phase_plots(self, magdata=None, title='Phase Map',
cmap='RdBu', limit=None, norm=None):
"""Display all phasemaps saved in the :class:`~.DataSet` as a colormesh.
Parameters
----------
magdata : :class:`~.VectorData`, optional
Magnetic distribution to which the projectors of the dataset should be applied. If not
given, the phasemaps 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 phase_plots')
if magdata is not None:
phasemaps = self.create_phasemaps(magdata)
else:
phasemaps = self.phasemaps
[phasemap.plot_phase('{} ({})'.format(title, self.projectors[i].get_info()),
cmap=cmap, limit=limit, norm=norm)
for (i, phasemap) in enumerate(phasemaps)]
def combined_plots(self, magdata=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
----------
magdata : :class:`~.VectorData`, optional
Magnetic distribution to which the projectors of the dataset should be applied. If not
given, the phasemaps 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 combined_plots')
if magdata is not None:
phasemaps = self.create_phasemaps(magdata)
else:
phasemaps = self.phasemaps
for (i, phasemap) in enumerate(phasemaps):
phasemap.plot_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 @ 68453aaa
# -*- 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.ravel()
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
-------
magdata_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
magdata_avg_kern = VectorData(self.cost.data_set.a, fft.zeros((3,) + self.dim))
magdata_avg_kern.set_vector(self.avrg_kern_row, mask=self.mask)
return magdata_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')
magdata_avg_kern = self.get_avg_kern_row(pos)
mag_x, mag_y, mag_z = magdata_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 @ 68453aaa
This diff is collapsed.
pyramid/fieldconverter.py deleted 100644 → 0
View file @ 68453aaa
This diff is collapsed.
pyramid/fielddata.py deleted 100644 → 0
View file @ 68453aaa
This diff is collapsed.
pyramid/file_io/__init__.py deleted 100644 → 0
View file @ 68453aaa
# -*- coding: utf-8 -*-
# Copyright 2016 by Forschungszentrum Juelich GmbH
# Author: J. Caron
#
"""Subpackage containing Pyramid IO functionality."""
from .io_phasemap import load_phasemap
from .io_vectordata import load_vectordata
from .io_scalardata import load_scalardata
__all__ = ['load_phasemap', 'load_vectordata', 'load_scalardata']
pyramid/file_io/io_phasemap.py deleted 100644 → 0
View file @ 68453aaa
This diff is collapsed.
pyramid/file_io/io_scalardata.py deleted 100644 → 0
View file @ 68453aaa
This diff is collapsed.
pyramid/file_io/io_vectordata.py deleted 100644 → 0
View file @ 68453aaa
This diff is collapsed.
pyramid/forwardmodel.py deleted 100644 → 0
View file @ 68453aaa
This diff is collapsed.
pyramid/kernel.py deleted 100644 → 0
View file @ 68453aaa
This diff is collapsed.
pyramid/magcreator/__init__.py deleted 100644 → 0
View file @ 68453aaa
# -*- coding: utf-8 -*-
# Copyright 2016 by Forschungszentrum Juelich GmbH
# Author: J. Caron
#
"""Subpackage containing functionality for creating magnetic distributions."""
from . import shapes
from . import examples
from .magcreator import *
__all__ = ['shapes', 'examples']
__all__.extend(magcreator.__all__)
pyramid/magcreator/examples.py deleted 100644 → 0
View file @ 68453aaa
This diff is collapsed.
pyramid/magcreator/magcreator.py deleted 100644 → 0
View file @ 68453aaa
This diff is collapsed.
pyramid/magcreator/shapes.py deleted 100644 → 0
View file @ 68453aaa
This diff is collapsed.
pyramid/phasemap.py deleted 100644 → 0
View file @ 68453aaa
This diff is collapsed.
pyramid/phasemapper.py deleted 100644 → 0
View file @ 68453aaa
This diff is collapsed.