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 5597 deletions
pyramid/__init__.py deleted 100644 → 0
View file @ 45d1150a
# -*- coding: utf-8 -*-
# Copyright 2016 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.
"""
from . import analytic
from . import reconstruction
from . import fieldconverter
from . import magcreator
from . import colors
from . import plottools # TODO: colors and plottools into "plots" package (maybe with examples?)
from . import utils
from .costfunction import *
from .dataset import *
from .diagnostics import *
from .fielddata import *
from .forwardmodel import *
from .kernel import *
from .phasemap import *
from .phasemapper import *
from .projector import *
from .regularisator import *
from .ramp import *
from .quaternion import *
from .file_io import *
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', 'fieldconverter',
'load_phasemap', 'load_vectordata', 'load_scalardata', 'load_projector', 'load_dataset',
'colors', 'utils']
__all__.extend(costfunction.__all__)
__all__.extend(dataset.__all__)
__all__.extend(diagnostics.__all__)
__all__.extend(forwardmodel.__all__)
__all__.extend(kernel.__all__)
__all__.extend(fielddata.__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(file_io.__all__)
# TODO: Test for different systems!
pyramid/analytic.py deleted 100644 → 0
View file @ 45d1150a
# -*- 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))
# TODO: During testing: "RuntimeWarning: invalid value encountered in power":
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/costfunction.py deleted 100644 → 0
View file @ 45d1150a
# -*- 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, 'Length of input {} does not match n={}'.format(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 @ 45d1150a
# -*- 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
from pyramid.fielddata import ScalarData
from pyramid.ramp import Ramp
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 phasemaps(self):
"""List of all PhaseMaps in the DataSet."""
return self._phasemaps
@property
def projectors(self):
"""List of all Projectors in the DataSet."""
return self._projectors
@property
def phasemappers(self):
# TODO: get rid, only use phasemapper_dict!!
"""List of all PhaseMappers in the DataSet."""
return self._phasemappers
@property
def phasemapper_dict(self):
"""Dictionary of all PhaseMappers in the DataSet."""
return self._phasemapper_dict
def __init__(self, a, dim, b_0=1, mask=None, Se_inv=None):
dim = tuple(dim)
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._phasemappers = []
self._phasemapper_dict = {}
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_single(self, phasemap, projector, phasemapper=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!'
dim_uv = projector.dim_uv
assert projector.dim == self.dim, '3D dimensions must match!'
assert phasemap.dim_uv == dim_uv, 'Projection dimensions (dim_uv) must match!'
assert phasemap.a == self.a, 'Grid spacing must match!'
# Create lookup key:
# TODO: Think again if phasemappers should be given as attribute (seems to be faulty
# TODO: currently... Also not very expensive, so keep outside?
if phasemapper is not None:
key = dim_uv # Create standard phasemapper, dim_uv is enough for identification!
else:
key = (dim_uv, str(phasemapper)) # Include string representation for identification!
# Retrieve existing, use given or create new phasemapper:
if key in self.phasemapper_dict: # Retrieve existing phasemapper:
phasemapper = self.phasemapper_dict[key]
elif phasemapper is not None: # Use given one (do nothing):
pass
else: # Create new standard (RDFC) phasemapper:
phasemapper = PhaseMapperRDFC(Kernel(self.a, dim_uv, self.b_0))
self._phasemapper_dict[key] = phasemapper
# Append everything to the lists (just contain pointers to objects!):
self._phasemaps.append(phasemap)
self._projectors.append(projector)
self._phasemappers.append(phasemapper)
def append(self, phasemap, projector, phasemapper=None):
"""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.
phasemapper: :class:`~.PhaseMapper`, optional
An optional :class:`~.PhaseMapper` object which should be added.
Returns
-------
None
"""
self._log.debug('Calling append')
if type(phasemap) is not list:
phasemap = [phasemap]
if type(projector) is not list:
projector = [projector]
if type(phasemapper) is not list:
phasemapper = [phasemapper] * len(phasemap)
assert len(phasemap) == len(projector),\
('Phasemaps and projectors must have same' +
'length(phasemaps: {}, projectors: {})!'.format(len(phasemap), len(projector)))
for i in range(len(phasemap)):
self._append_single(phasemap[i], projector[i], phasemapper[i])
# Reset the Se_inv matrix from phasemaps confidence matrices:
self.set_Se_inv_diag_with_conf()
def create_phasemaps(self, magdata, difference=False, ramp=None):
"""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.
difference : bool, optional
If `True`, the phasemaps of the dataset are subtracted from the created ones to view
difference images. Default is False.
ramp : :class:`~.Ramp`
A ramp object, which can be specified to add a ramp to the generated phasemaps.
If `difference` is `True`, this can be interpreted as ramp correcting the phasemaps
saved in the dataset.
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 i, projector in enumerate(self.projectors):
mag_proj = projector(magdata)
phasemap = self.phasemappers[i](mag_proj)
if difference:
phasemap -= self.phasemaps[i]
if ramp is not None:
assert type(ramp) == Ramp, 'ramp has to be a Ramp object!'
phasemap += ramp(index=i) # Full formula: phasemap -= phasemap_dataset - ramp
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, threshold=0.9):
# TODO: This function should be in a separate module and not here (maybe?)!
"""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.
threshold: float, optional
The threshold, describing the minimal number of 2D masks which have to extrude to the
point in 3D to be considered valid as containing magnetisation. `threshold` is a
relative number in the range of [0, 1]. The default is 0.9. Choosing a value of 1 is
the strictest possible setting (every 2D mask has to contain a 3D point to be valid).
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 = np.zeros(self.dim)
for i, projector in enumerate(self.projectors):
mask_2d = self.phasemaps[i].mask.reshape(-1) # 2D mask
# Add extrusion of 2D mask:
mask_3d += projector.weight.T.dot(mask_2d).reshape(self.dim)
self.mask = np.where(mask_3d >= threshold * self.count, True, False)
def save(self, filename, overwrite=True):
"""Saves the dataset as a collection of HDF5 files.
Parameters
----------
filename: str
Base name of the files which the dataset is saved into. HDF5 files are supported.
overwrite: bool, optional
If True (default), an existing file will be overwritten, if False, this
(silently!) does nothing.
"""
from .file_io.io_dataset import save_dataset
save_dataset(self, filename, overwrite)
def plot_mask(self):
"""If it exists, display the 3D mask of the magnetization distribution.
Returns
-------
None
"""
self._log.debug('Calling plot_mask')
if self.mask is not None:
return ScalarData(self.a, self.mask).plot_mask()
def plot_phasemaps(self, magdata=None, title='Phase Map', difference=False, ramp=None,
**kwargs):
"""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.
difference : bool, optional
If `True`, the phasemaps of the dataset are subtracted from the created ones to view
difference images. Default is False.
ramp : :class:`~.Ramp`
A ramp object, which can be specified to add a ramp to the generated phasemaps.
If `magdata` is not given, this will instead just ramp correct the phasemaps saved
in the dataset.
Returns
-------
None
"""
self._log.debug('Calling plot_phasemaps')
if magdata is not None: # Plot phasemaps of the given magnetisation distribution:
phasemaps = self.create_phasemaps(magdata, difference=difference, ramp=ramp)
else: # Plot phasemaps saved in the DataSet (default):
phasemaps = self.phasemaps
if ramp is not None:
for i, phasemap in enumerate(phasemaps):
assert type(ramp) == Ramp, 'ramp has to be a Ramp object!'
phasemap -= ramp(index=i) # Ramp correction
for (i, phasemap) in enumerate(phasemaps):
phasemap.plot_phase(note='{} ({})'.format(title, self.projectors[i].get_info()),
**kwargs)
def plot_phasemaps_combined(self, magdata=None, title='Combined Plot', difference=False,
ramp=None, **kwargs):
"""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'.
difference : bool, optional
If `True`, the phasemaps of the dataset are subtracted from the created ones to view
difference images. Default is False.
ramp : :class:`~.Ramp`
A ramp object, which can be specified to add a ramp to the generated phasemaps.
If `magdata` is not given, this will instead just ramp correct the phasemaps saved
in the dataset.
Returns
-------
None
"""
self._log.debug('Calling plot_phasemaps_combined')
if magdata is not None:
phasemaps = self.create_phasemaps(magdata, difference=difference, ramp=ramp)
else:
phasemaps = self.phasemaps
if ramp is not None:
for i, phasemap in enumerate(phasemaps):
assert type(ramp) == Ramp, 'ramp has to be a Ramp object!'
phasemap -= ramp(index=i) # Ramp correction
for (i, phasemap) in enumerate(phasemaps):
phasemap.plot_combined(note='{} ({})'.format(title, self.projectors[i].get_info()),
**kwargs)
plt.show()
pyramid/diagnostics.py deleted 100644 → 0
View file @ 45d1150a
This diff is collapsed.
pyramid/fft.py deleted 100644 → 0
View file @ 45d1150a
# -*- 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.
"""
# TODO: is this still in use or deprecated by jutil???
import pickle
import logging
import os
import numpy as np
_log = logging.getLogger(__name__)
try:
import pyfftw
BACKEND = 'fftw'
except ImportError:
pyfftw = None
BACKEND = 'numpy'
_log.info('pyFFTW module not found. Using numpy implementation.')
try:
import multiprocessing
NTHREADS = multiprocessing.cpu_count()
del multiprocessing
except ImportError:
NTHREADS = 1
_log.info('multiprocessing module not found. Using single core.')
__all__ = ['plans', 'FLOAT', 'COMPLEX', 'dump_wisdom', 'load_wisdom',
'zeros', 'empty', 'ones', 'configure_backend',
'fftn', 'ifftn', 'rfftn', 'irfftn', 'rfftn_adj', 'irfftn_adj']
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):
"""Add an FFTW object to the cache.
Parameters
----------
fft_type: basestring
Identifier sting for the FFT type ('fftn', 'ifftn', 'rfftn', 'irfftn').
fftw_obj: :class:`~pyfftw.FFTW` object
The FFTW object which should be added to the cache.
s: tuple of ints
Shape of the output array.
axes: tuple of ints
The axes along which the FFTW should be executed.
nthreads: int
Number of threads which should be used.
"""
self._log.debug('Calling add_fftw')
in_arr = fftw_obj.get_input_array()
key = (fft_type, in_arr.shape, in_arr.dtype, s, axes, nthreads)
self.cache[key] = fftw_obj
def lookup_fftw(self, fft_type, in_arr, s, axes, nthreads):
"""
Parameters
----------
fft_type: basestring
Identifier sting for the FFT type ('fftn', 'ifftn', 'rfftn', 'irfftn').
in_arr:
Input array, internally, just the `dtype` and the `shape` are used to identify the FFT.
s: tuple of ints
Shape of the output array.
axes: tuple of ints
The axes along which the FFTW should be executed.
nthreads: int
Number of threads which should be used.
Returns
-------
fftw_obj: :class:`~pyfftw.FFTW` object
The requested FFTW object.
"""
self._log.debug('Calling lookup_fftw')
key = (fft_type, in_arr.shape, in_arr.dtype, s, axes, nthreads)
return self.cache.get(key, None)
def clear_cache(self):
"""Clear the cache."""
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[:, :a.shape[1]] = 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):
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):
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):
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):
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:
ifftn = None
fftn = None
rfftn = None
irfftn = None
rfftn_adj = None
irfftn_adj = None
if pyfftw is not None:
configure_backend('fftw')
else:
configure_backend('numpy')
pyramid/fieldconverter.py deleted 100644 → 0
View file @ 45d1150a
This diff is collapsed.
pyramid/fielddata.py deleted 100644 → 0
View file @ 45d1150a
This diff is collapsed.
pyramid/file_io/__init__.py deleted 100644 → 0
View file @ 45d1150a
# -*- 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
from .io_projector import load_projector
from .io_dataset import load_dataset
__all__ = ['load_phasemap', 'load_vectordata', 'load_scalardata', 'load_projector', 'load_dataset']
pyramid/file_io/io_dataset.py deleted 100644 → 0
View file @ 45d1150a
This diff is collapsed.
pyramid/file_io/io_phasemap.py deleted 100644 → 0
View file @ 45d1150a
This diff is collapsed.
pyramid/file_io/io_projector.py deleted 100644 → 0
View file @ 45d1150a
This diff is collapsed.
pyramid/file_io/io_scalardata.py deleted 100644 → 0
View file @ 45d1150a
This diff is collapsed.
pyramid/file_io/io_vectordata.py deleted 100644 → 0
View file @ 45d1150a
This diff is collapsed.
pyramid/forwardmodel.py deleted 100644 → 0
View file @ 45d1150a
This diff is collapsed.
pyramid/kernel.py deleted 100644 → 0
View file @ 45d1150a
This diff is collapsed.
pyramid/magcreator/__init__.py deleted 100644 → 0
View file @ 45d1150a
# -*- 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 @ 45d1150a
This diff is collapsed.
pyramid/magcreator/magcreator.py deleted 100644 → 0
View file @ 45d1150a
This diff is collapsed.
pyramid/magcreator/shapes.py deleted 100644 → 0
View file @ 45d1150a
This diff is collapsed.