diff --git a/pyramid/__init__.py b/pyramid/__init__.py
index ad7b174a2c94b02c324c649df6a51f502bab6a3d..47aba58821b746abcd170574a4a3d3b64c66364a 100644
--- a/pyramid/__init__.py
+++ b/pyramid/__init__.py
@@ -30,6 +30,8 @@ reconstruction
Reconstruct magnetic distributions from given phasemaps.
regularisator
Class to instantiate different regularisation strategies.
+diagnostics
+ Class to calculate diagnostics
fft
Class for custom FFT functions using numpy or FFTW.
@@ -41,31 +43,33 @@ numcore
"""
-import logging
-
from . import analytic
from . import magcreator
from . import reconstruction
from . import fft
-from .costfunction import *
-from .dataset import *
-from .forwardmodel import *
-from .kernel import *
-from .magdata import *
-from .phasemap import *
-from .phasemapper import *
-from .projector import *
-from .regularisator import *
+from .costfunction import * # analysis:ignore
+from .dataset import * # analysis:ignore
+from .diagnostics import * # analysis:ignore
+from .forwardmodel import * # analysis:ignore
+from .kernel import * # analysis:ignore
+from .magdata import * # analysis:ignore
+from .phasemap import * # analysis:ignore
+from .phasemapper import * # analysis:ignore
+from .projector import * # analysis:ignore
+from .regularisator import * # analysis:ignore
from .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__ = ['__version__', '__hg_revision__', 'analytic', 'magcreator', 'reconstruction', 'fft']
__all__.extend(costfunction.__all__)
__all__.extend(dataset.__all__)
+__all__.extend(diagnostics.__all__)
__all__.extend(forwardmodel.__all__)
__all__.extend(kernel.__all__)
__all__.extend(magdata.__all__)
diff --git a/pyramid/analytic.py b/pyramid/analytic.py
index 5faef3196c67b2b28bf60f54d894cf88de517f3d..3585c897cd0da6f2ee06667e4444ec532d2860af 100644
--- a/pyramid/analytic.py
+++ b/pyramid/analytic.py
@@ -1,4 +1,7 @@
# -*- 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
diff --git a/pyramid/costfunction.py b/pyramid/costfunction.py
index 50c41ff688b3231e6dc916e6b3d6ca5517f8f0ee..fb67e4539166f37a49e046c0903fe4d6163b7751 100644
--- a/pyramid/costfunction.py
+++ b/pyramid/costfunction.py
@@ -1,4 +1,7 @@
# -*- 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."""
@@ -6,7 +9,6 @@ the so called `cost` of a threedimensional magnetization distribution."""
import numpy as np
from scipy.sparse import eye as sparse_eye
-from scipy.sparse.linalg import LinearOperator
from pyramid.forwardmodel import ForwardModel
from pyramid.regularisator import NoneRegularisator
@@ -34,7 +36,6 @@ class Costfunction(object):
:class:`~dataset.DataSet` object, which stores all information for the cost calculation.
regularisator : :class:`~.Regularisator`
Regularisator class that's responsible for the regularisation term.
- regularisator: :class:`~Regularisator`
y : :class:`~numpy.ndarray` (N=1)
Vector which lists all pixel values of all phase maps one after another.
fwd_model : :class:`~.ForwardModel`
@@ -131,70 +132,12 @@ class Costfunction(object):
+ self.regularisator.hess_dot(x, vector))
def hess_diag(self, _):
- # TODO: Docstring!
- # TODO: What for again?
- return np.ones(self.n)
-
- def estimate_lambda(self):
- # TODO: Docstring!
- # TODO: Not very efficient? Is this even correct?
-
- def unit_vec(length, index):
- result = np.zeros(length)
- result[index] = 1
- return result
-
- fwd, reg = self.fwd_model, self.regularisator
- trace_fwd = np.sum([fwd.jac_T_dot(None, fwd.jac_dot(None, unit_vec(self.n, i)))
- for i in range(self.n)])
- trace_reg = np.sum([reg(unit_vec(self.n, i)) for i in range(self.n)])
- print 'fwd:', trace_fwd
- print 'reg:', trace_reg
- import pdb; pdb.set_trace()
- return trace_fwd / trace_reg
-
-
-class CFAdapterScipyCG(LinearOperator):
-
- '''Adapter class making the :class:`~.Costfunction` class accessible for scipy cg methods.
-
- This class provides an adapter for the :class:`~.Costfunction` to be usable with the
- :func:`~.scipy.sparse.linalg.cg` function. the :func:`~.matvec` function is overwritten to
- implement a multiplication with the Hessian of the adapted costfunction. This is used in the
- :func:`~pyramid.reconstruction.optimise_sparse_cg` function of the
- :mod:`~pyramid.reconstruction` module.
-
- Attributes
- ----------
- cost : :class:`~.Costfunction`
- Costfunction which should be made usable in the :func:`~.scipy.sparse.linalg.cg` function.
-
- '''
- # TODO: make obsolete!
-
- _log = logging.getLogger(__name__+'.CFAdapterScipyCG')
-
- def __init__(self, cost):
- self._log.debug('Calling __init__')
- self.cost = cost
-
- def matvec(self, vector):
- '''Matrix-vector multiplication with the Hessian of the adapted costfunction.
+ ''' Return the diagonal of the Hessian.
Parameters
----------
- vector : :class:`~numpy.ndarray` (N=1)
- Vector which will be multiplied by the Hessian matrix provided by the adapted
- costfunction.
+ _ : undefined
+ Unused input
'''
- self._log.debug('Calling matvec')
- return self.cost.hess_dot(None, vector)
-
- @property
- def shape(self):
- return (self.cost.data_set.n, self.cost.data_set.n)
-
- @property
- def dtype(self):
- return np.dtype("d")
+ return np.ones(self.n)
diff --git a/pyramid/dataset.py b/pyramid/dataset.py
index 9f422ca59f1e6cda07006b5dca9bd881fd82bacc..a1924fb1aac69d03a3802053dac09c558b81c533 100644
--- a/pyramid/dataset.py
+++ b/pyramid/dataset.py
@@ -1,4 +1,7 @@
# -*- 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."""
@@ -49,6 +52,9 @@ class DataSet(object):
A list of all stored :class:`~.PhaseMap` objects.
phase_vec: :class:`~numpy.ndarray` (N=1)
The concatenaded, vectorized phase of all ;class:`~.PhaseMap` objects.
+ Se_inv : :class:`~numpy.ndarray` (N=2), optional
+ Inverted covariance matrix of the measurement errors. The matrix has size `NxN` with N
+ being the length of the targetvector y (vectorized phase map information).
m: int
Size of the image space.
n: int
@@ -76,11 +82,10 @@ class DataSet(object):
@property
def phase_mappers(self):
dim_uv_set = set([p.dim_uv for p in self.projectors])
- kernel_list = [Kernel(self.a, dim_uv, threads=self.threads) for dim_uv in dim_uv_set]
+ 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, use_fftw=True, threads=1):
- # TODO: document!
+ def __init__(self, a, dim, b_0=1, mask=None, Se_inv=None):
self._log.debug('Calling __init__')
assert isinstance(a, Number), 'Grid spacing has to be a number!'
assert a >= 0, 'Grid spacing has to be a positive number!'
@@ -96,8 +101,6 @@ class DataSet(object):
self.b_0 = b_0
self.mask = mask
self.Se_inv = Se_inv
- self.use_fftw = use_fftw
- self.threads = threads
self.phase_maps = []
self.projectors = []
self._log.debug('Created: '+str(self))
@@ -152,11 +155,33 @@ class DataSet(object):
for projector in self.projectors]
def set_Se_inv_block_diag(self, cov_list):
- # TODO: Document!
+ '''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.Se_inv = sparse.block_diag(cov_list).tocsr()
def set_Se_inv_diag_with_masks(self, mask_list):
- # TODO: Document!
+ '''Set the Se_inv matrix as a block diagonal matrix from a list of masks
+
+ Parameters
+ ----------
+ mask_list: list of :class:`~numpy.ndarray`
+ List of 2D masks (one for each projection) which define trust regions.
+
+ Returns
+ -------
+ None
+
+ '''
cov_list = [sparse.diags(m.flatten(), 0) for m in mask_list]
self.set_Se_inv_block_diag(cov_list)
diff --git a/pyramid/diagnostics.py b/pyramid/diagnostics.py
index fbec1ff0b72b94a88643c1228e2083c2d554296d..b9d84bd1db845cc3e867973058845a38e972768d 100644
--- a/pyramid/diagnostics.py
+++ b/pyramid/diagnostics.py
@@ -1,12 +1,13 @@
# -*- coding: utf-8 -*-
-"""
-Created on Thu Dec 11 17:20:27 2014
-
-@author: Jan
-"""
+# Copyright 2014 by Forschungszentrum Juelich GmbH
+# Author: J. Caron
+#
+"""This module provides the :class:`~.Diagnostics` class for the calculation of diagnostics of a
+specified costfunction for a fixed magnetization distribution."""
import numpy as np
+import matplotlib.pyplot as plt
import jutil
@@ -14,50 +15,80 @@ from pyramid import fft
from pyramid.magdata import MagData
from pyramid.phasemap import PhaseMap
+
+__all__ = ['Diagnostics']
+
+
class Diagnostics(object):
- # TODO: Docstrings and position of properties!
+ '''Class for calculating diagnostic properties of a specified costfunction.
- def __init__(self, x_rec, cost, max_iter=100):
- self.x_rec = x_rec
- self.cost = cost
- self.max_iter = max_iter
- self.fwd_model = cost.fwd_model
- self.Se_inv = self.cost.Se_inv
- self.dim = self.cost.data_set.dim
- self.row_idx = None
- self.set_position(0)#(0, self.dim[0]//2, self.dim[1]//2, self.dim[2]//2))
- self._A = jutil.operator.CostFunctionOperator(self.cost, self.x_rec)
- self._P = jutil.preconditioner.CostFunctionPreconditioner(self.cost, self.x_rec)
+ 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.
- def set_position(self, pos):
- # TODO: Does not know about the mask... thus gives wrong results or errors
-# m, z, y, x = pos
-# row_idx = m*np.prod(self.dim) + z*self.dim[1]*self.dim[2] + y*self.dim[2] + x
- row_idx = pos
- if row_idx != self.row_idx:
- self.row_idx = row_idx
- self._updated_std = False
- self._updated_gain_row = False
- self._updated_avrg_kern_row = False
- self._updated_measure_contribution = False
+ Attributes
+ ----------
+ x_rec: :class:`~numpy.ndarray`
+ Vectorized magnetization distribution at which the costfunction is evaluated.
+ cost: :class:`~.pyramid.costfunction.Costfunction`
+ Costfunction for which the diagnostics are calculated.
+ max_iter: int, optional
+ Maximum number of iterations. Default is 1000.
+ fwd_model: :class:`~pyramid.forwardmodel.ForwardModel`
+ Forward model used in the costfunction.
+ Se_inv : :class:`~numpy.ndarray` (N=2), optional
+ Inverted covariance matrix of the measurement errors. The matrix has size `NxN` with N
+ being the length of the targetvector y (vectorized phase map information).
+ dim: tuple (N=3)
+ Dimensions of the 3D magnetic distribution.
+ row_idx: int
+ Row index of the system matrix corresponding to the current position in 3D space.
+ cov_row: :class:`~numpy.ndarray`
+ Row of the covariance matrix (``S_a^-1+F'(x_f)^T S_e^-1 F'(x_f)``) which is needed for the
+ calculation of the gain and averaging kernel matrizes and which ideally contains the
+ variance at position `row_idx` for the current component and position in 3D.
+ std: float
+ Standard deviation of the chosen component at the current position (calculated when
+ needed).
+ gain_row: :class:`~numpy.ndarray`
+ Row of the gain matrix, which maps differences of phase measurements onto differences in
+ the retrieval result of the magnetization distribution(calculated when needed).
+ avrg_kern_row: :class:`~numpy.ndarray`
+ Row of the averaging kernel matrix (which is ideally the identity matrix), which describes
+ the smoothing introduced by the regularization (calculated when needed).
+ measure_contribution: :class:`~numpy.ndarray`
+
+ Notes
+ -----
+ Some properties depend on others, which may require recalculations of these prior
+ properties if necessary. The dependencies are ('-->' == 'requires'):
+ avrg_kern_row --> gain_row --> std --> m_inv_row
+ measure_contribution is independant
+
+ '''
@property
- def std(self):
- if not self._updated_std:
+ def cov_row(self):
+ if not self._updated_cov_row:
e_i = fft.zeros(self.cost.n, dtype=fft.FLOAT)
e_i[self.row_idx] = 1
- row = jutil.cg.conj_grad_solve(self._A, e_i, P=self._P, max_iter=self.max_iter)
- self._m_inv_row = row
- self._std = np.sqrt(self._m_inv_row[self.row_idx])
- self._updated_std = True
- return self._std
+ row = 2 * jutil.cg.conj_grad_solve(self._A, e_i, P=self._P, max_iter=self.max_iter)
+ self._std_row = row
+ self._updated_cov_row = True
+ return self._std_row
+
+ @property
+ def std(self):
+ return np.sqrt(self.cov_row[self.row_idx])
@property
def gain_row(self):
if not self._updated_gain_row:
- self.std # evoke to update self._m_inv_row if necessary # TODO: make _m_inv_row checked!
- self._gain_row = self.Se_inv.dot(self.fwd_model.jac_dot(self.x_rec, self._m_inv_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
@@ -73,24 +104,142 @@ class Diagnostics(object):
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 = jutil.cg.conj_grad_solve(self._A, cache, P=self._P, max_iter=self.max_iter)
+ 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
- def get_avg_kernel(self, pos=None):
+ @property
+ def pos(self):
+ return self._pos
+
+ @pos.setter
+ def pos(self, pos):
+ c, z, y, x = pos
+ assert self.mask[z, y, x], 'Position is outside of the provided mask!'
+ mask_vec = self.mask.flatten()
+ idx_3d = z*self.dim[1]*self.dim[2] + y*self.dim[2] + x
+ row_idx = c*np.prod(mask_vec.sum()) + mask_vec[:idx_3d].sum()
+ if row_idx != self.row_idx:
+ self._pos = pos
+ self.row_idx = row_idx
+ self._updated_cov_row = False
+ self._updated_gain_row = False
+ self._updated_avrg_kern_row = False
+ self._updated_measure_contribution = False
+
+ def __init__(self, x_rec, cost, max_iter=1000):
+ self.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)
+
+ def get_avg_kern_row(self, pos=None):
+ '''Get the averaging kernel matrix row represented as a 3D magnetization distribution.
+
+ Parameters
+ ----------
+ pos: tuple (N=4)
+ Position in 3D plus component `(c, z, y, x)`
+
+ Returns
+ -------
+ mag_data_avg_kern: :class:`~pyramid.magdata.MagData`
+ Averaging kernel matrix row represented as a 3D magnetization distribution
+
+ '''
if pos is not None:
- self.set_position(pos)
+ self.pos = pos
mag_data_avg_kern = MagData(self.cost.data_set.a, fft.zeros((3,)+self.dim))
- mag_data_avg_kern.set_vector(self.avrg_kern_row, mask=self.cost.data_set.mask)
+ mag_data_avg_kern.set_vector(self.avrg_kern_row, mask=self.mask)
return mag_data_avg_kern
- def get_gain_maps(self, pos=None):
+ def calculate_averaging(self, pos=None):
+ '''Get the gain matrix row represented as a list of 2D 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.
+
+ '''# TODO: Docstring!
+ mag_data_avg_kern = self.get_avg_kern_row(pos)
+ print self.pos
+ mag_x, mag_y, mag_z = mag_data_avg_kern.magnitude
+ x = mag_x.sum(axis=(0, 1))
+ y = mag_y.sum(axis=(0, 2))
+ z = mag_z.sum(axis=(1, 2))
+ plt.figure()
+ plt.axhline(y=0, ls='-', color='k')
+ plt.axhline(y=1, ls='-', color='k')
+ plt.plot(x, label='x', color='r', marker='o')
+ plt.plot(y, label='y', color='g', marker='o')
+ plt.plot(z, label='z', color='b', marker='o')
+ c = self.pos[0]
+ data = [x, y, z][c]
+ col = ['r', 'g', 'b'][c]
+ i_m = np.argmax(data) # Index of the maximum
+ plt.axhline(y=data[i_m], ls='-', color=col)
+ plt.axhline(y=data[i_m]/2, ls='--', color=col)
+ # Left side:
+ for i in np.arange(i_m-1, -1, -1):
+ if data[i] < data[i_m]/2:
+ l = (data[i_m]/2-data[i])/(data[i+1]-data[i]) + i
+ break
+ # Right side:
+ for i in np.arange(i_m+1, data.size):
+ if data[i] < data[i_m]/2:
+ r = (data[i_m]/2-data[i-1])/(data[i]-data[i-1]) + i-1
+ break
+ # Calculate FWHM:
+ fwhm = r - l
+ px_avrg = 1 / fwhm
+ plt.vlines(x=[l, r], ymin=0, ymax=data[i_m]/2, linestyles=':', color=col)
+ # Add legend:
+ plt.legend()
+ return px_avrg, fwhm, (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).
+
+ '''
if pos is not None:
- self.set_position(pos)
+ self.pos = pos
hp = self.cost.data_set.hook_points
- result = []
+ 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)
- result.append(PhaseMap(self.cost.data_set.a, gain))
- return result
+ gain_map_list.append(PhaseMap(self.cost.data_set.a, gain))
+ return gain_map_list
diff --git a/pyramid/fft.py b/pyramid/fft.py
index 1c7733d1427b4a973532069563221d76dc3e7358..5baea7fe47ee6e90ad1234ed0f79534c581e98fb 100644
--- a/pyramid/fft.py
+++ b/pyramid/fft.py
@@ -1,11 +1,16 @@
# -*- coding: utf-8 -*-
-"""
-Created on Fri Nov 28 15:30:10 2014
+# 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.
-@author: Jan
"""
-# TODO: Document!
import numpy as np
import cPickle as pickle
@@ -27,9 +32,20 @@ except ImportError:
BACKEND = 'numpy'
print("pyFFTW module not found. Using numpy implementation.")
+try:
+ import psutil
+ NTHREADS = psutil.cpu_count()
+except ImportError:
+ try:
+ import multiprocessing
+ NTHREADS = multiprocessing.cpu_count()
+ except ImportError:
+ NTHREADS = 1
+
-__all__ = ['PLANS', 'FLOAT', 'COMPLEX', 'NTHREADS', 'dump_wisdom', 'load_wisdom', 'zeros', 'empty',
- 'configure_backend', 'fftn', 'ifftn', 'rfftn', 'irfftn', 'rfftn_adj', 'irfftn_adj']
+__all__ = ['PLANS', 'FLOAT', 'COMPLEX', 'NTHREADS', 'dump_wisdom', 'load_wisdom', #analysis:ignore
+ 'zeros', 'empty', 'configure_backend',
+ 'fftn', 'ifftn', 'rfftn', 'irfftn', 'rfftn_adj', 'irfftn_adj']
class FFTWCache(object):
@@ -53,7 +69,6 @@ class FFTWCache(object):
PLANS = FFTWCache()
FLOAT = np.float32 # One convenient place to
COMPLEX = np.complex64 # change from 32 to 64 bit
-NTHREADS = 1
# Numpy functions:
@@ -98,7 +113,7 @@ def _fftn_fftw(a, s=None, axes=None):
raise TypeError('Wrong input type!')
fftw = PLANS.lookup_fftw('fftn', a, s, axes, NTHREADS)
if fftw is None:
- fftw = pyfftw.builders.fftn(a, s, axes)
+ fftw = pyfftw.builders.fftn(a, s, axes, threads=NTHREADS)
PLANS.add_fftw('fftn', fftw, s, axes, NTHREADS)
return fftw(a).copy()
@@ -108,7 +123,7 @@ def _ifftn_fftw(a, s=None, axes=None):
raise TypeError('Wrong input type!')
fftw = PLANS.lookup_fftw('ifftn', a, s, axes, NTHREADS)
if fftw is None:
- fftw = pyfftw.builders.ifftn(a, s, axes)
+ fftw = pyfftw.builders.ifftn(a, s, axes, threads=NTHREADS)
PLANS.add_fftw('ifftn', fftw, s, axes, NTHREADS)
return fftw(a).copy()
@@ -118,7 +133,7 @@ def _rfftn_fftw(a, s=None, axes=None):
raise TypeError('Wrong input type!')
fftw = PLANS.lookup_fftw('rfftn', a, s, axes, NTHREADS)
if fftw is None:
- fftw = pyfftw.builders.rfftn(a, s, axes)
+ fftw = pyfftw.builders.rfftn(a, s, axes, threads=NTHREADS)
PLANS.add_fftw('rfftn', fftw, s, axes, NTHREADS)
return fftw(a).copy()
@@ -128,7 +143,7 @@ def _irfftn_fftw(a, s=None, axes=None):
raise TypeError('Wrong input type!')
fftw = PLANS.lookup_fftw('irfftn', a, s, axes, NTHREADS)
if fftw is None:
- fftw = pyfftw.builders.irfftn(a, s, axes)
+ fftw = pyfftw.builders.irfftn(a, s, axes, threads=NTHREADS)
PLANS.add_fftw('irfftn', fftw, s, axes, NTHREADS)
return fftw(a).copy()
@@ -154,14 +169,36 @@ def _irfftn_adj_fftw(a):
# These wisdom functions do nothing if pyFFTW is not available
def dump_wisdom(fname):
- # TODO: Docstring!
+ '''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
+
+ '''
if pyfftw is not None:
with open(fname, 'wb') as fp:
pickle.dump(pyfftw.export_wisdom(), fp, pickle.HIGHEST_PROTOCOL)
def load_wisdom(fname):
- # TODO: Docstring!
+ '''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
+
+ '''
if pyfftw is not None:
if not os.path.exists(fname):
print("Warning: Wisdom file does not exist. First time use?")
@@ -172,7 +209,21 @@ def load_wisdom(fname):
# Array setups:
def empty(shape, dtype=FLOAT):
- # TODO: Docstring!
+ '''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.
+
+ '''
result = np.empty(shape, dtype)
if pyfftw is not None:
result = pyfftw.n_byte_align(result, pyfftw.simd_alignment)
@@ -180,7 +231,21 @@ def empty(shape, dtype=FLOAT):
def zeros(shape, dtype=FLOAT):
- # TODO: Docstring!
+ '''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.
+
+ '''
result = np.zeros(shape, dtype)
if pyfftw is not None:
result = pyfftw.n_byte_align(result, pyfftw.simd_alignment)
@@ -188,7 +253,21 @@ def zeros(shape, dtype=FLOAT):
def ones(shape, dtype=FLOAT):
- # TODO: Docstring!
+ '''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.
+
+ '''
result = np.ones(shape, dtype)
if pyfftw is not None:
result = pyfftw.n_byte_align(result, pyfftw.simd_alignment)
@@ -197,11 +276,18 @@ def ones(shape, dtype=FLOAT):
# Configure backend:
def configure_backend(backend):
- """
- Change FFT backend.
+ '''Change FFT backend.
+
+ Parameters
+ ----------
+ backend: string
+ Backend to use. Supported values are "numpy" and "fftw".
+
+ Returns
+ -------
+ None
- Supported values are "numpy" and "fftw".
- """
+ '''
global fftn
global ifftn
global rfftn
diff --git a/pyramid/forwardmodel.py b/pyramid/forwardmodel.py
index f2758b90fd2eca9287a8b754897f9d9f73ea9607..d7f0e6e69f53de8fc12ba7314203acee9ababa5a 100644
--- a/pyramid/forwardmodel.py
+++ b/pyramid/forwardmodel.py
@@ -1,4 +1,7 @@
# -*- coding: utf-8 -*-
+# Copyright 2014 by Forschungszentrum Juelich GmbH
+# Author: J. Caron
+#
"""This module provides the :class:`~.ForwardModel` class which represents a strategy to map a
threedimensional magnetization distribution onto a two-dimensional phase map."""
@@ -65,7 +68,6 @@ class ForwardModel(object):
self._log.debug('Calling __call__')
self.mag_data.magnitude[:] = 0
self.mag_data.set_vector(x, self.data_set.mask)
- # TODO: Multiprocessing
result = np.zeros(self.m)
hp = self.hook_points
for i, projector in enumerate(self.data_set.projectors):
diff --git a/pyramid/kernel.py b/pyramid/kernel.py
index 8c322c5f788bd636916c62423f2b8fbd8bcabd56..b8bc451f628d96df658fc8375932331f98df0532 100644
--- a/pyramid/kernel.py
+++ b/pyramid/kernel.py
@@ -1,4 +1,7 @@
# -*- coding: utf-8 -*-
+# Copyright 2014 by Forschungszentrum Juelich GmbH
+# Author: J. Caron
+#
"""This module provides the :class:`~.Kernel` class, representing the phase contribution of one
single magnetized pixel."""
@@ -34,6 +37,14 @@ class Kernel(object):
dim_uv : tuple of int (N=2), optional
Dimensions of the 2-dimensional projected magnetization grid from which the phase should
be calculated.
+ dim_kern : tuple of int (N=2)
+ Dimensions of the kernel, which is ``2N-1`` for both axes compared to `dim_uv`.
+ dim_pad : tuple of int (N=2)
+ Dimensions of the padded FOV, which is ``2N`` (if FFTW is used) or the next highest power
+ of 2 (for numpy-FFT).
+ dim_fft : tuple of int (N=2)
+ Dimensions of the grid, which is used for the FFT, taking into account that a RFFT should
+ be used (one axis is halved in comparison to `dim_pad`).
b_0 : float, optional
Saturation magnetization in Tesla, which is used for the phase calculation. Default is 1.
geometry : {'disc', 'slab'}, optional
@@ -46,20 +57,19 @@ class Kernel(object):
The real FFT of the phase contribution of one pixel magnetized in u-direction.
v_fft : :class:`~numpy.ndarray` (N=3)
The real FFT of the phase contribution of one pixel magnetized in v-direction.
- dim_fft : tuple of int (N=2)
- Dimensions of the grid, which is used for the FFT. Calculated by adding the dimensions
- `dim_uv` of the magnetization grid and the dimensions of the kernel (given by
- ``2*dim_uv-1``)
- and increasing to the next multiple of 2 (for faster FFT).
- slice_fft : tuple (N=2) of :class:`slice`
- A tuple of :class:`slice` objects to extract the original field of view from the increased
- size (`size_fft`) of the grid for the FFT-convolution.
+ slice_phase : tuple (N=2) of :class:`slice`
+ A tuple of :class:`slice` objects to extract the original FOV from the increased one with
+ size `dim_pad` for the elementary kernel phase. The kernel is shifted, thus the center is
+ not at (0, 0), which also shifts the slicing compared to `slice_mag`.
+ slice_mag : tuple (N=2) of :class:`slice`
+ A tuple of :class:`slice` objects to extract the original FOV from the increased one with
+ size `dim_pad` for the projected magnetization distribution.
- ''' # TODO: overview what all dim_??? mean! and use_fftw, slice(_fft), etc.
+ '''
_log = logging.getLogger(__name__+'.Kernel')
- def __init__(self, a, dim_uv, b_0=1., geometry='disc', threads=1):
+ def __init__(self, a, dim_uv, b_0=1., geometry='disc'):
self._log.debug('Calling __init__')
# Set basic properties:
self.dim_uv = dim_uv # Dimensions of the FOV
@@ -76,7 +86,6 @@ class Kernel(object):
slice(dim_uv[1]-1, self.dim_kern[1])) # is not at (0, 0)!
self.slice_mag = (slice(0, dim_uv[0]), # Magnetization is padded on the far end!
slice(0, dim_uv[1])) # (Phase cutout is shifted as listed above)
- self.threads = threads # TODO: make obsolete!
# Calculate kernel (single pixel phase):
coeff = b_0 * a**2 / (2*PHI_0) # Minus is gone because of negative z-direction
v_dim, u_dim = dim_uv
@@ -116,7 +125,17 @@ class Kernel(object):
- F_a(n-0.5, m+0.5) + F_a(n+0.5, m+0.5))
def print_info(self):
- # TODO: Docstring!
+ '''Print information about the kernel.
+
+ Parameters
+ ----------
+ None
+
+ Returns
+ -------
+ None
+
+ '''
self._log.debug('Calling print_info')
print 'Shape of the FOV :', self.dim_uv
print 'Shape of the Kernel:', self.dim_kern
@@ -126,4 +145,3 @@ class Kernel(object):
print 'Slice for the magn.:', self.slice_mag
print 'Grid spacing: {} nm'.format(self.a)
print 'Geometry:', self.geometry
- print 'Use FFTW: {}; with {} thread(s)'.format(self.use_fftw, self.threads)
diff --git a/pyramid/magcreator.py b/pyramid/magcreator.py
index c3abb46ee1e495976d1ceb1eb55bbee57b75d4e2..3083d459c9c5b3bc9895cb15b39b2afb852b2e37 100644
--- a/pyramid/magcreator.py
+++ b/pyramid/magcreator.py
@@ -1,4 +1,7 @@
# -*- coding: utf-8 -*-
+# Copyright 2014 by Forschungszentrum Juelich GmbH
+# Author: J. Caron
+#
"""Create simple magnetic distributions.
The :mod:`~.magcreator` module is responsible for the creation of simple distributions of
@@ -165,7 +168,7 @@ class Shapes(object):
elif axis == 'x':
mag_shape = np.array([[[np.hypot((y-center[1])/(width[1]/2.),
(z-center[0])/(width[0]/2.)) <= 1
- and abs(z - center[0]) <= height / 2
+ and abs(x - center[2]) <= height / 2
for x in range(dim[2])]
for y in range(dim[1])]
for z in range(dim[0])])
@@ -225,9 +228,9 @@ class Shapes(object):
assert np.shape(dim) == (3,), 'Parameter dim has to be a a tuple of length 3!'
assert np.shape(center) == (3,), 'Parameter center has to be a a tuple of length 3!'
assert np.shape(width) == (3,), 'Parameter width has to be a a tuple of length 3!'
- mag_shape = np.array([[[np.sqrt((x-center[2])**2/(width[2]/2)**2
- + (y-center[1])**2/(width[1]/2)**2
- + (z-center[0])**2/(width[0]/2)**2) <= 1
+ mag_shape = np.array([[[(x-center[2])**2/(width[2]/2)**2
+ + (y-center[1])**2/(width[1]/2)**2
+ + (z-center[0])**2/(width[0]/2)**2 <= 1
for x in range(dim[2])]
for y in range(dim[1])]
for z in range(dim[0])])
@@ -388,6 +391,3 @@ def create_mag_dist_vortex(mag_shape, center=None, axis='z', magnitude=1):
y_mag = -np.ones(dim) * np.cos(phi) * mag_shape * magnitude
x_mag = np.zeros(dim)
return np.array([x_mag, y_mag, z_mag])
-
-# TODO: Smooth Vortex!
-# m(x,y) = (-y/r cos(phi(r/R)), x/r sin(phi(r/R)), cos(phi(r))); phi in range[0, 1]
\ No newline at end of file
diff --git a/pyramid/magdata.py b/pyramid/magdata.py
index 7fe3e6cae13bdf7ebd9003b0ca82e1d0ff8e1d72..53f4d435e0842fbd0f9c7e35f3edbaf66a6c9a70 100644
--- a/pyramid/magdata.py
+++ b/pyramid/magdata.py
@@ -1,640 +1,638 @@
-# -*- coding: utf-8 -*-
-"""This module provides the :class:`~.MagData` class for storing of magnetization data."""
-
-
-import os
-
-import numpy as np
-from numpy.linalg import norm
-from scipy.ndimage.interpolation import zoom
-
-import matplotlib.pyplot as plt
-import matplotlib.cm as cmx
-from matplotlib.ticker import MaxNLocator
-
-from numbers import Number
-
-import netCDF4
-
-import logging
-
-
-__all__ = ['MagData']
-
-
-class MagData(object):
-
- '''Class for storing magnetization data.
-
- Represents 3-dimensional magnetic distributions with 3 components which are stored as a
- 2-dimensional numpy array in `magnitude`, but which can also be accessed as a vector via
- `mag_vec`. :class:`~.MagData` objects support negation, arithmetic operators
- (``+``, ``-``, ``*``) and their augmented counterparts (``+=``, ``-=``, ``*=``), with numbers
- and other :class:`~.MagData` objects, if their dimensions and grid spacings match. It is
- possible to load data from NetCDF4 or LLG (.txt) files or to save the data in these formats.
- Plotting methods are also provided.
-
- Attributes
- ----------
- a: float
- The grid spacing in nm.
- dim: tuple (N=3)
- Dimensions (z, y, x) of the grid.
- magnitude: :class:`~numpy.ndarray` (N=4)
- The `x`-, `y`- and `z`-component of the magnetization vector for every 3D-gridpoint
- as a 4-dimensional numpy array (first dimension has to be 3, because of the 3 components).
- mag_vec: :class:`~numpy.ndarray` (N=1)
- Vector containing the magnetic distribution.
-
- '''
-
- _log = logging.getLogger(__name__+'.MagData')
-
- @property
- def a(self):
- return self._a
-
- @a.setter
- def a(self, a):
- assert isinstance(a, Number), 'Grid spacing has to be a number!'
- assert a >= 0, 'Grid spacing has to be a positive number!'
- self._a = float(a)
-
- @property
- def dim(self):
- return self._dim
-
- @property
- def magnitude(self):
- return self._magnitude
-
- @magnitude.setter
- def magnitude(self, magnitude):
- assert isinstance(magnitude, np.ndarray), 'Magnitude has to be a numpy array!'
- assert len(magnitude.shape) == 4, 'Magnitude has to be 4-dimensional!'
- assert magnitude.shape[0] == 3, 'First dimension of the magnitude has to be 3!'
- self._magnitude = np.asarray(magnitude, dtype=np.float32)
- self._dim = magnitude.shape[1:]
-
- @property
- def mag_vec(self):
- return np.reshape(self.magnitude, -1)
-
- @mag_vec.setter
- def mag_vec(self, mag_vec):
- assert isinstance(mag_vec, np.ndarray), 'Vector has to be a numpy array!'
- assert np.size(mag_vec) == 3*np.prod(self.dim), \
- 'Vector has to match magnitude dimensions! {} {}'.format(mag_vec.shape,
- 3*np.prod(self.dim))
- self.magnitude = mag_vec.reshape((3,)+self.dim)
-
- def __init__(self, a, magnitude):
- self._log.debug('Calling __init__')
- self.a = a
- self.magnitude = magnitude
- self._log.debug('Created '+str(self))
-
- def __repr__(self):
- self._log.debug('Calling __repr__')
- return '%s(a=%r, magnitude=%r)' % (self.__class__, self.a, self.magnitude)
-
- def __str__(self):
- self._log.debug('Calling __str__')
- return 'MagData(a=%s, dim=%s)' % (self.a, self.dim)
-
- def __neg__(self): # -self
- self._log.debug('Calling __neg__')
- return MagData(self.a, -self.magnitude)
-
- def __add__(self, other): # self + other
- self._log.debug('Calling __add__')
- assert isinstance(other, (MagData, Number)), \
- 'Only MagData objects and scalar numbers (as offsets) can be added/subtracted!'
- if isinstance(other, MagData):
- self._log.debug('Adding two MagData objects')
- assert other.a == self.a, 'Added phase has to have the same grid spacing!'
- assert other.magnitude.shape == (3,)+self.dim, \
- 'Added magnitude has to have the same dimensions!'
- return MagData(self.a, self.magnitude+other.magnitude)
- else: # other is a Number
- self._log.debug('Adding an offset')
- return MagData(self.a, self.magnitude+other)
-
- def __sub__(self, other): # self - other
- self._log.debug('Calling __sub__')
- return self.__add__(-other)
-
- def __mul__(self, other): # self * other
- self._log.debug('Calling __mul__')
- assert isinstance(other, Number), 'MagData objects can only be multiplied by numbers!'
- return MagData(self.a, other*self.magnitude)
-
- def __radd__(self, other): # other + self
- self._log.debug('Calling __radd__')
- return self.__add__(other)
-
- def __rsub__(self, other): # other - self
- self._log.debug('Calling __rsub__')
- return -self.__sub__(other)
-
- def __rmul__(self, other): # other * self
- self._log.debug('Calling __rmul__')
- return self.__mul__(other)
-
- def __iadd__(self, other): # self += other
- self._log.debug('Calling __iadd__')
- return self.__add__(other)
-
- def __isub__(self, other): # self -= other
- self._log.debug('Calling __isub__')
- return self.__sub__(other)
-
- def __imul__(self, other): # self *= other
- self._log.debug('Calling __imul__')
- return self.__mul__(other)
-
- def copy(self):
- '''Returns a copy of the :class:`~.MagData` object
-
- Parameters
- ----------
- None
-
- Returns
- -------
- mag_data: :class:`~.MagData`
- A copy of the :class:`~.MagData`.
-
- '''
- self._log.debug('Calling copy')
- return MagData(self.a, self.magnitude.copy())
-
- def scale_down(self, n=1):
- '''Scale down the magnetic distribution by averaging over two pixels along each axis.
-
- Parameters
- ----------
- n : int, optional
- Number of times the magnetic distribution is scaled down. The default is 1.
-
- Returns
- -------
- None
-
- Notes
- -----
- Acts in place and changes dimensions and grid spacing accordingly.
- Only possible, if each axis length is a power of 2!
-
- '''
- self._log.debug('Calling scale_down')
- assert n > 0 and isinstance(n, (int, long)), 'n must be a positive integer!'
- self.a = self.a * 2**n
- for t in range(n):
- # Pad if necessary:
- pz, py, px = self.dim[0] % 2, self.dim[1] % 2, self.dim[2] % 2
- if pz != 0 or py != 0 or px != 0:
- self.magnitude = np.pad(self.magnitude, ((0, 0), (0, pz), (0, py), (0, px)),
- mode='constant')
- # Create coarser grid for the magnetization:
- self.magnitude = self.magnitude.reshape(
- 3, self.dim[0]/2, 2, self.dim[1]/2, 2, self.dim[2]/2, 2).mean(axis=(6, 4, 2))
-
- def scale_up(self, n=1, order=0):
- '''Scale up the magnetic distribution using spline interpolation of the requested order.
-
- Parameters
- ----------
- n : int, optional
- Power of 2 with which the grid is scaled. Default is 1, which means every axis is
- increased by a factor of ``2**1 = 2``.
- order : int, optional
- The order of the spline interpolation, which has to be in the range between 0 and 5
- and defaults to 0.
-
- Returns
- -------
- None
-
- Notes
- -----
- Acts in place and changes dimensions and grid spacing accordingly.
- '''
- self._log.debug('Calling scale_up')
- assert n > 0 and isinstance(n, (int, long)), 'n must be a positive integer!'
- assert 5 > order >= 0 and isinstance(order, (int, long)), \
- 'order must be a positive integer between 0 and 5!'
- self.a = self.a / 2**n
- self.magnitude = np.array((zoom(self.magnitude[0], zoom=2**n, order=order),
- zoom(self.magnitude[1], zoom=2**n, order=order),
- zoom(self.magnitude[2], zoom=2**n, order=order)))
-
- def pad(self, x_pad, y_pad, z_pad):
- '''Pad the current magnetic distribution with zeros for each individual axis.
-
- Parameters
- ----------
- x_pad : int
- Number of zeros which should be padded on both sides of the x-axis.
- y_pad : int
- Number of zeros which should be padded on both sides of the y-axis.
- z_pad : int
- Number of zeros which should be padded on both sides of the z-axis.
-
- Returns
- -------
- None
-
- Notes
- -----
- Acts in place and changes dimensions accordingly.
- '''
- assert x_pad >= 0 and isinstance(x_pad, (int, long)), 'x_pad must be a positive integer!'
- assert y_pad >= 0 and isinstance(y_pad, (int, long)), 'y_pad must be a positive integer!'
- assert z_pad >= 0 and isinstance(z_pad, (int, long)), 'z_pad must be a positive integer!'
- self.magnitude = np.pad(self.magnitude,
- ((0, 0), (z_pad, z_pad), (y_pad, y_pad), (x_pad, x_pad)),
- mode='constant')
-
- def get_mask(self, threshold=0):
- '''Mask all pixels where the amplitude of the magnetization lies above `threshold`.
-
- Parameters
- ----------
- threshold : float, optional
- A pixel only gets masked, if it lies above this threshold . The default is 0.
-
- Returns
- -------
- mask : :class:`~numpy.ndarray` (N=3, boolean)
- Mask of the pixels where the amplitude of the magnetization lies above `threshold`.
-
- '''
- return np.sqrt(np.sum(np.array(self.magnitude)**2, axis=0)) > threshold
-
- def get_vector(self, mask):
- '''Returns the magnetic components arranged in a vector, specified by a mask.
-
- Parameters
- ----------
- mask : :class:`~numpy.ndarray` (N=3, boolean)
- Masks the pixels from which the components should be taken.
-
- Returns
- -------
- vector : :class:`~numpy.ndarray` (N=1)
- The vector containing magnetization components of the specified pixels.
- Order is: first all `x`-, then all `y`-, then all `z`-components.
-
- '''
- if mask is not None:
- return np.reshape([self.magnitude[0][mask],
- self.magnitude[1][mask],
- self.magnitude[2][mask]], -1)
- else:
- return self.mag_vec
-
-
- def set_vector(self, vector, mask=None):
- '''Set the magnetic components of the masked pixels to the values specified by `vector`.
-
- Parameters
- ----------
- mask : :class:`~numpy.ndarray` (N=3, boolean), optional
- Masks the pixels from which the components should be taken.
- vector : :class:`~numpy.ndarray` (N=1)
- The vector containing magnetization components of the specified pixels.
- Order is: first all `x`-, then all `y-, then all `z`-components.
-
- Returns
- -------
- None
-
- '''
- assert np.size(vector) % 3 == 0, 'Vector has to contain all 3 components for every pixel!'
- count = np.size(vector)/3
- if mask is not None:
- self.magnitude[0][mask] = vector[:count] # x-component
- self.magnitude[1][mask] = vector[count:2*count] # y-component
- self.magnitude[2][mask] = vector[2*count:] # z-component
- else:
- self.mag_vec = vector
-
- def save_to_llg(self, filename='..\output\magdata_output.txt'):
- '''Save magnetization data in a file with LLG-format.
-
- Parameters
- ----------
- filename : string, optional
- The name of the LLG-file in which to store the magnetization data.
- The default is '..\output\magdata_output.txt'.
-
- Returns
- -------
- None
-
- '''
- self._log.debug('Calling save_to_llg')
- a = self.a * 1.0E-9 / 1.0E-2 # from nm to cm
- # Create 3D meshgrid and reshape it and the magnetization into a list where x varies first:
- zz, yy, xx = np.mgrid[a/2:(self.dim[0]*a-a/2):self.dim[0]*1j,
- a/2:(self.dim[1]*a-a/2):self.dim[1]*1j,
- a/2:(self.dim[2]*a-a/2):self.dim[2]*1j].reshape(3, -1)
- x_vec, y_vec, z_vec = self.magnitude.reshape(3, -1)
- # Save data to file:
- data = np.array([xx, yy, zz, x_vec, y_vec, z_vec]).T
- with open(filename, 'w') as mag_file:
- mag_file.write('LLGFileCreator: %s\n' % filename.replace('.txt', ''))
- mag_file.write(' %d %d %d\n' % (self.dim[2], self.dim[1], self.dim[0]))
- mag_file.writelines('\n'.join(' '.join('{:7.6e}'.format(cell)
- for cell in row) for row in data))
-
- @classmethod
- def load_from_llg(cls, filename):
- '''Construct :class:`~.MagData` object from LLG-file.
-
- Parameters
- ----------
- filename : string
- The name of the LLG-file from which to load the data.
-
- Returns
- -------
- mag_data: :class:`~.MagData`
- A :class:`~.MagData` object containing the loaded data.
-
- '''
- cls._log.debug('Calling load_from_llg')
- SCALE = 1.0E-9 / 1.0E-2 # From cm to nm
- data = np.genfromtxt(filename, skip_header=2)
- dim = tuple(np.genfromtxt(filename, dtype=int, skip_header=1, skip_footer=len(data[:, 0])))
- a = (data[1, 0] - data[0, 0]) / SCALE
- magnitude = data[:, 3:6].T.reshape((3,)+dim)
- return MagData(a, magnitude)
-
- def save_to_netcdf4(self, filename='..\output\magdata_output.nc'):
- '''Save magnetization data in a file with NetCDF4-format.
-
- Parameters
- ----------
- filename : string, optional
- The name of the NetCDF4-file in which to store the magnetization data.
- The default is '..\output\magdata_output.nc'.
-
- Returns
- -------
- None
-
- '''
- self._log.debug('Calling save_to_netcdf4')
- mag_file = netCDF4.Dataset(filename, 'w', format='NETCDF4')
- mag_file.a = self.a
- mag_file.createDimension('comp', 3) # Number of components
- mag_file.createDimension('z_dim', self.dim[0])
- mag_file.createDimension('y_dim', self.dim[1])
- mag_file.createDimension('x_dim', self.dim[2])
- magnitude = mag_file.createVariable('magnitude', 'f', ('comp', 'z_dim', 'y_dim', 'x_dim'))
- magnitude[...] = self.magnitude
- mag_file.close()
-
- @classmethod
- def load_from_netcdf4(cls, filename):
- '''Construct :class:`~.DataMag` object from NetCDF4-file.
-
- Parameters
- ----------
- filename : string
- The name of the NetCDF4-file from which to load the data. Standard format is '\*.nc'.
-
- Returns
- -------
- mag_data: :class:`~.MagData`
- A :class:`~.MagData` object containing the loaded data.
-
- '''
- cls._log.debug('Calling copy')
- mag_file = netCDF4.Dataset(filename, 'r', format='NETCDF4')
- a = mag_file.a
- magnitude = mag_file.variables['magnitude'][...]
- mag_file.close()
- return MagData(a, magnitude)
-
- def quiver_plot(self, title='Magnetization Distribution', axis=None, proj_axis='z',
- ar_dens=1, ax_slice=None, log=False, scaled=True): # TODO: Doc ar_dens
- '''Plot a slice of the magnetization as a quiver plot.
-
- Parameters
- ----------
- title : string, optional
- The title for the plot.
- axis : :class:`~matplotlib.axes.AxesSubplot`, optional
- Axis on which the graph is plotted. Creates a new figure if none is specified.
- proj_axis : {'z', 'y', 'x'}, optional
- The axis, from which a slice is plotted. The default is 'z'.
- ax_slice : int, optional
- The slice-index of the axis specified in `proj_axis`. Is set to the center of
- `proj_axis` if not specified.
- log : boolean, optional
- Takes the Default is False.
- scaled : boolean, optional
- Normalizes the plotted arrows in respect to the highest one. Default is True.
-
- Returns
- -------
- axis: :class:`~matplotlib.axes.AxesSubplot`
- The axis on which the graph is plotted.
-
- '''
- self._log.debug('Calling quiver_plot')
- assert proj_axis == 'z' or proj_axis == 'y' or proj_axis == 'x', \
- 'Axis has to be x, y or z (as string).'
- if proj_axis == 'z': # Slice of the xy-plane with z = ax_slice
- self._log.debug('proj_axis == z')
- if ax_slice is None:
- self._log.debug('ax_slice is None')
- ax_slice = int(self.dim[0]/2.)
- plot_u = np.copy(self.magnitude[0][ax_slice, ...]) # x-component
- plot_v = np.copy(self.magnitude[1][ax_slice, ...]) # y-component
- u_label = 'x [px]'
- v_label = 'y [px]'
- elif proj_axis == 'y': # Slice of the xz-plane with y = ax_slice
- self._log.debug('proj_axis == y')
- if ax_slice is None:
- self._log.debug('ax_slice is None')
- ax_slice = int(self.dim[1]/2.)
- plot_u = np.copy(self.magnitude[0][:, ax_slice, :]) # x-component
- plot_v = np.copy(self.magnitude[2][:, ax_slice, :]) # z-component
- u_label = 'x [px]'
- v_label = 'z [px]'
- elif proj_axis == 'x': # Slice of the yz-plane with x = ax_slice
- self._log.debug('proj_axis == x')
- if ax_slice is None:
- self._log.debug('ax_slice is None')
- ax_slice = int(self.dim[2]/2.)
- plot_u = np.swapaxes(np.copy(self.magnitude[2][..., ax_slice]), 0, 1) # z-component
- plot_v = np.swapaxes(np.copy(self.magnitude[1][..., ax_slice]), 0, 1) # y-component
- u_label = 'z [px]'
- v_label = 'y [px]'
- # If no axis is specified, a new figure is created:
- if axis is None:
- self._log.debug('axis is None')
- fig = plt.figure(figsize=(8.5, 7))
- axis = fig.add_subplot(1, 1, 1)
- axis.set_aspect('equal')
- angles = np.angle(plot_u+1j*plot_v, deg=True)
- # Take the logarithm of the arrows to clearly show directions (if specified):
- if log:
- cutoff = 10
- amp = np.round(np.hypot(plot_u, plot_v), decimals=cutoff)
- min_value = amp[np.nonzero(amp)].min()
- plot_u = np.round(plot_u, decimals=cutoff) / min_value
- plot_u = np.log10(np.abs(plot_u)+1) * np.sign(plot_u)
- plot_v = np.round(plot_v, decimals=cutoff) / min_value
- plot_v = np.log10(np.abs(plot_v)+1) * np.sign(plot_v)
- # Scale the magnitude of the arrows to the highest one (if specified):
- if scaled:
- plot_u /= np.hypot(plot_u, plot_v).max()
- plot_v /= np.hypot(plot_u, plot_v).max()
- # Setup quiver:
- dim_uv = plot_u.shape
- ad = ar_dens
- xx, yy = np.meshgrid(np.arange(dim_uv[1]), np.arange(dim_uv[0]))
- axis.quiver(xx[::ad, ::ad], yy[::ad, ::ad], plot_u[::ad, ::ad], plot_v[::ad, ::ad],
- pivot='middle', units='xy', angles=angles[::ad, ::ad], minlength=0.25,
- scale_units='xy', scale=1./ad, headwidth=6, headlength=7)
- axis.set_xlim(-1, dim_uv[1])
- axis.set_ylim(-1, dim_uv[0])
- axis.set_title(title, fontsize=18)
- axis.set_xlabel(u_label, fontsize=15)
- axis.set_ylabel(v_label, fontsize=15)
- axis.tick_params(axis='both', which='major', labelsize=14)
- axis.xaxis.set_major_locator(MaxNLocator(nbins=9, integer=True))
- axis.yaxis.set_major_locator(MaxNLocator(nbins=9, integer=True))
- return axis
-
- # TODO: Switch with mayavi or combine
- def quiver_plot3d_matplotlib(self, title='Magnetization Distribution', ar_dens=1): # TODO: Doc ar_dens
- from mpl_toolkits.mplot3d import axes3d #analysis:ignore
- # TODO: more arguments like in the other plots and document!!!
- a = self.a
- dim = self.dim
- ad = ar_dens
- # Create points and vector components as lists:
- zz, yy, xx = np.mgrid[a/2:(dim[0]*a-a/2):dim[0]*1j,
- a/2:(dim[1]*a-a/2):dim[1]*1j,
- a/2:(dim[2]*a-a/2):dim[2]*1j]
- xx = xx[::ad, ::ad, ::ad]
- yy = yy[::ad, ::ad, ::ad]
- zz = zz[::ad, ::ad, ::ad]
- x_mag = self.magnitude[0, ::ad, ::ad, ::ad]
- y_mag = self.magnitude[1, ::ad, ::ad, ::ad]
- z_mag = self.magnitude[2, ::ad, ::ad, ::ad]
- # Plot them as vectors:
- fig = plt.figure(figsize=(8.5, 7))
- axis = fig.add_subplot(1, 1, 1)
- axis = fig.gca(projection='3d')
- axis.quiver(xx, yy, zz, x_mag, y_mag, z_mag)
- axis.set_title(title, fontsize=18)
- # TODO: add colorbar!
-# mlab.colorbar(None, label_fmt='%.2f')
-# mlab.colorbar(None, orientation='vertical')
- return axis
-
- def quiver_plot3d(self, title='Magnetization Distribution', ar_dens=1): # TODO: Doc ar_dens
- '''Plot the magnetization as 3D-vectors in a quiverplot.
-
- Parameters
- ----------
- None
-
- Returns
- -------
- None
-
- '''
- self._log.debug('Calling quiver_plot3D')
- from mayavi import mlab
- a = self.a
- dim = self.dim
- ad = ar_dens
- # Create points and vector components as lists:
- zz, yy, xx = np.mgrid[a/2:(dim[0]*a-a/2):dim[0]*1j,
- a/2:(dim[1]*a-a/2):dim[1]*1j,
- a/2:(dim[2]*a-a/2):dim[2]*1j]
- xx = xx[::ad, ::ad, ::ad].flatten()#reshape(-1)
- yy = yy[::ad, ::ad, ::ad].flatten()#.reshape(-1)
- zz = zz[::ad, ::ad, ::ad].flatten()#.reshape(-1)
- x_mag = self.magnitude[0][::ad, ::ad, ::ad].flatten()#, (-1))
- y_mag = self.magnitude[1][::ad, ::ad, ::ad].flatten()#, (-1))
- z_mag = self.magnitude[2][::ad, ::ad, ::ad].flatten()#, (-1))
- # Plot them as vectors:
- mlab.figure()
- plot = mlab.quiver3d(xx, yy, zz, x_mag, y_mag, z_mag, mode='arrow')
- mlab.outline(plot)
- mlab.axes(plot)
- mlab.title(title, height=0.95, size=0.35)
- mlab.colorbar(None, label_fmt='%.2f')
- mlab.colorbar(None, orientation='vertical')
-
- def save_to_x3d(self, filename='..\..\output\magdata_output.x3d', maximum=1):
- '''Output the magnetization in the .x3d format for the Fraunhofer InstantReality Player.
-
- Parameters
- ----------
- None
-
- Returns
- -------
- None
-
- '''
- self._log.debug('Calling save_to_x3d')
- from lxml import etree
-
- dim = self.dim
- # Create points and vector components as lists:
- zz, yy, xx = np.mgrid[0.5:(dim[0]-0.5):dim[0]*1j,
- 0.5:(dim[1]-0.5):dim[1]*1j,
- 0.5:(dim[2]-0.5):dim[2]*1j]
- xx = xx.reshape(-1)
- yy = yy.reshape(-1)
- zz = zz.reshape(-1)
- x_mag = np.reshape(self.magnitude[0], (-1))
- y_mag = np.reshape(self.magnitude[1], (-1))
- z_mag = np.reshape(self.magnitude[2], (-1))
- # Load template, load tree and write viewpoint information:
- template = os.path.join(os.path.dirname(os.path.realpath(__file__)), 'template.x3d')
- parser = etree.XMLParser(remove_blank_text=True)
- tree = etree.parse(template, parser)
- scene = tree.find('Scene')
- etree.SubElement(scene, 'Viewpoint', position='0 0 {}'.format(1.5*dim[0]),
- fieldOfView='1')
- # Write each "spin"-tag separately:
- for i in range(np.prod(dim)):
- mag = np.sqrt(x_mag[i]**2+y_mag[i]**2+z_mag[i]**2)
- if mag != 0:
- spin_position = (xx[i]-dim[2]/2., yy[i]-dim[1]/2., zz[i]-dim[0]/2.)
- sx_ref = 0
- sy_ref = 1
- sz_ref = 0
- rot_x = sy_ref*z_mag[i] - sz_ref*y_mag[i]
- rot_y = sz_ref*x_mag[i] - sx_ref*z_mag[i]
- rot_z = sx_ref*y_mag[i] - sy_ref*x_mag[i]
- angle = np.arccos(y_mag[i]/mag)
- if norm((rot_x, rot_y, rot_z)) < 1E-10:
- rot_x, rot_y, rot_z = 1, 0, 0
- spin_rotation = (rot_x, rot_y, rot_z, angle)
- spin_color = cmx.RdYlGn(mag/maximum)[:3]
- spin_scale = (1., 1., 1.)
- spin = etree.SubElement(scene, 'ProtoInstance',
- DEF='Spin {}'.format(i), name='Spin_Proto')
- etree.SubElement(spin, 'fieldValue', name='spin_position',
- value='{} {} {}'.format(*spin_position))
- etree.SubElement(spin, 'fieldValue', name='spin_rotation',
- value='{} {} {} {}'.format(*spin_rotation))
- etree.SubElement(spin, 'fieldValue', name='spin_color',
- value='{} {} {}'.format(*spin_color))
- etree.SubElement(spin, 'fieldValue', name='spin_scale',
- value='{} {} {}'.format(*spin_scale))
- # Write the tree into the file in pretty print format:
- tree.write(filename, pretty_print=True)
+# -*- coding: utf-8 -*-
+# Copyright 2014 by Forschungszentrum Juelich GmbH
+# Author: J. Caron
+#
+"""This module provides the :class:`~.MagData` class for storing of magnetization data."""
+
+
+import os
+
+import numpy as np
+from numpy.linalg import norm
+from scipy.ndimage.interpolation import zoom
+
+import matplotlib.pyplot as plt
+import matplotlib.cm as cmx
+from matplotlib.ticker import MaxNLocator
+
+import pyramid.fft as fft
+
+from numbers import Number
+
+import netCDF4
+
+import logging
+
+
+__all__ = ['MagData']
+
+
+class MagData(object):
+
+ '''Class for storing magnetization data.
+
+ Represents 3-dimensional magnetic distributions with 3 components which are stored as a
+ 2-dimensional numpy array in `magnitude`, but which can also be accessed as a vector via
+ `mag_vec`. :class:`~.MagData` objects support negation, arithmetic operators
+ (``+``, ``-``, ``*``) and their augmented counterparts (``+=``, ``-=``, ``*=``), with numbers
+ and other :class:`~.MagData` objects, if their dimensions and grid spacings match. It is
+ possible to load data from NetCDF4 or LLG (.txt) files or to save the data in these formats.
+ Plotting methods are also provided.
+
+ Attributes
+ ----------
+ a: float
+ The grid spacing in nm.
+ dim: tuple (N=3)
+ Dimensions (z, y, x) of the grid.
+ magnitude: :class:`~numpy.ndarray` (N=4)
+ The `x`-, `y`- and `z`-component of the magnetization vector for every 3D-gridpoint
+ as a 4-dimensional numpy array (first dimension has to be 3, because of the 3 components).
+ mag_amp: :class:`~numpy.ndarray` (N=3)
+ The length (amplitude) of the magnetization vectors as a 3D-array.
+ mag_vec: :class:`~numpy.ndarray` (N=1)
+ Vector containing the magnetic distribution.
+
+ '''
+
+ _log = logging.getLogger(__name__+'.MagData')
+
+ @property
+ def a(self):
+ return self._a
+
+ @a.setter
+ def a(self, a):
+ assert isinstance(a, Number), 'Grid spacing has to be a number!'
+ assert a >= 0, 'Grid spacing has to be a positive number!'
+ self._a = float(a)
+
+ @property
+ def dim(self):
+ return self._dim
+
+ @property
+ def magnitude(self):
+ return self._magnitude
+
+ @magnitude.setter
+ def magnitude(self, magnitude):
+ assert isinstance(magnitude, np.ndarray), 'Magnitude has to be a numpy array!'
+ assert len(magnitude.shape) == 4, 'Magnitude has to be 4-dimensional!'
+ assert magnitude.shape[0] == 3, 'First dimension of the magnitude has to be 3!'
+ self._magnitude = np.asarray(magnitude, dtype=fft.FLOAT)
+ self._dim = magnitude.shape[1:]
+
+ @property
+ def mag_amp(self):
+ return np.sqrt(np.sum(self.magnitude**2, axis=0))
+
+ @property
+ def mag_vec(self):
+ return np.reshape(self.magnitude, -1)
+
+ @mag_vec.setter
+ def mag_vec(self, mag_vec):
+ mag_vec = np.asarray(mag_vec, dtype=fft.FLOAT)
+ assert np.size(mag_vec) == 3*np.prod(self.dim), \
+ 'Vector has to match magnitude dimensions! {} {}'.format(mag_vec.shape,
+ 3*np.prod(self.dim))
+ self.magnitude = mag_vec.reshape((3,)+self.dim)
+
+ def __init__(self, a, magnitude):
+ self._log.debug('Calling __init__')
+ self.a = a
+ self.magnitude = magnitude
+ self._log.debug('Created '+str(self))
+
+ def __repr__(self):
+ self._log.debug('Calling __repr__')
+ return '%s(a=%r, magnitude=%r)' % (self.__class__, self.a, self.magnitude)
+
+ def __str__(self):
+ self._log.debug('Calling __str__')
+ return 'MagData(a=%s, dim=%s)' % (self.a, self.dim)
+
+ def __neg__(self): # -self
+ self._log.debug('Calling __neg__')
+ return MagData(self.a, -self.magnitude)
+
+ def __add__(self, other): # self + other
+ self._log.debug('Calling __add__')
+ assert isinstance(other, (MagData, Number)), \
+ 'Only MagData objects and scalar numbers (as offsets) can be added/subtracted!'
+ if isinstance(other, MagData):
+ self._log.debug('Adding two MagData objects')
+ assert other.a == self.a, 'Added phase has to have the same grid spacing!'
+ assert other.magnitude.shape == (3,)+self.dim, \
+ 'Added magnitude has to have the same dimensions!'
+ return MagData(self.a, self.magnitude+other.magnitude)
+ else: # other is a Number
+ self._log.debug('Adding an offset')
+ return MagData(self.a, self.magnitude+other)
+
+ def __sub__(self, other): # self - other
+ self._log.debug('Calling __sub__')
+ return self.__add__(-other)
+
+ def __mul__(self, other): # self * other
+ self._log.debug('Calling __mul__')
+ assert isinstance(other, Number), 'MagData objects can only be multiplied by numbers!'
+ return MagData(self.a, other*self.magnitude)
+
+ def __radd__(self, other): # other + self
+ self._log.debug('Calling __radd__')
+ return self.__add__(other)
+
+ def __rsub__(self, other): # other - self
+ self._log.debug('Calling __rsub__')
+ return -self.__sub__(other)
+
+ def __rmul__(self, other): # other * self
+ self._log.debug('Calling __rmul__')
+ return self.__mul__(other)
+
+ def __iadd__(self, other): # self += other
+ self._log.debug('Calling __iadd__')
+ return self.__add__(other)
+
+ def __isub__(self, other): # self -= other
+ self._log.debug('Calling __isub__')
+ return self.__sub__(other)
+
+ def __imul__(self, other): # self *= other
+ self._log.debug('Calling __imul__')
+ return self.__mul__(other)
+
+ def copy(self):
+ '''Returns a copy of the :class:`~.MagData` object
+
+ Parameters
+ ----------
+ None
+
+ Returns
+ -------
+ mag_data: :class:`~.MagData`
+ A copy of the :class:`~.MagData`.
+
+ '''
+ self._log.debug('Calling copy')
+ return MagData(self.a, self.magnitude.copy())
+
+ def scale_down(self, n=1):
+ '''Scale down the magnetic distribution by averaging over two pixels along each axis.
+
+ Parameters
+ ----------
+ n : int, optional
+ Number of times the magnetic distribution is scaled down. The default is 1.
+
+ Returns
+ -------
+ None
+
+ Notes
+ -----
+ Acts in place and changes dimensions and grid spacing accordingly.
+ Only possible, if each axis length is a power of 2!
+
+ '''
+ self._log.debug('Calling scale_down')
+ assert n > 0 and isinstance(n, (int, long)), 'n must be a positive integer!'
+ self.a = self.a * 2**n
+ for t in range(n):
+ # Pad if necessary:
+ pz, py, px = self.dim[0] % 2, self.dim[1] % 2, self.dim[2] % 2
+ if pz != 0 or py != 0 or px != 0:
+ self.magnitude = np.pad(self.magnitude, ((0, 0), (0, pz), (0, py), (0, px)),
+ mode='constant')
+ # Create coarser grid for the magnetization:
+ self.magnitude = self.magnitude.reshape(
+ 3, self.dim[0]/2, 2, self.dim[1]/2, 2, self.dim[2]/2, 2).mean(axis=(6, 4, 2))
+
+ def scale_up(self, n=1, order=0):
+ '''Scale up the magnetic distribution using spline interpolation of the requested order.
+
+ Parameters
+ ----------
+ n : int, optional
+ Power of 2 with which the grid is scaled. Default is 1, which means every axis is
+ increased by a factor of ``2**1 = 2``.
+ order : int, optional
+ The order of the spline interpolation, which has to be in the range between 0 and 5
+ and defaults to 0.
+
+ Returns
+ -------
+ None
+
+ Notes
+ -----
+ Acts in place and changes dimensions and grid spacing accordingly.
+ '''
+ self._log.debug('Calling scale_up')
+ assert n > 0 and isinstance(n, (int, long)), 'n must be a positive integer!'
+ assert 5 > order >= 0 and isinstance(order, (int, long)), \
+ 'order must be a positive integer between 0 and 5!'
+ self.a = self.a / 2**n
+ self.magnitude = np.array((zoom(self.magnitude[0], zoom=2**n, order=order),
+ zoom(self.magnitude[1], zoom=2**n, order=order),
+ zoom(self.magnitude[2], zoom=2**n, order=order)))
+
+ def pad(self, x_pad, y_pad, z_pad):
+ '''Pad the current magnetic distribution with zeros for each individual axis.
+
+ Parameters
+ ----------
+ x_pad : int
+ Number of zeros which should be padded on both sides of the x-axis.
+ y_pad : int
+ Number of zeros which should be padded on both sides of the y-axis.
+ z_pad : int
+ Number of zeros which should be padded on both sides of the z-axis.
+
+ Returns
+ -------
+ None
+
+ Notes
+ -----
+ Acts in place and changes dimensions accordingly.
+ '''
+ assert x_pad >= 0 and isinstance(x_pad, (int, long)), 'x_pad must be a positive integer!'
+ assert y_pad >= 0 and isinstance(y_pad, (int, long)), 'y_pad must be a positive integer!'
+ assert z_pad >= 0 and isinstance(z_pad, (int, long)), 'z_pad must be a positive integer!'
+ self.magnitude = np.pad(self.magnitude,
+ ((0, 0), (z_pad, z_pad), (y_pad, y_pad), (x_pad, x_pad)),
+ mode='constant')
+
+ def get_mask(self, threshold=0):
+ '''Mask all pixels where the amplitude of the magnetization lies above `threshold`.
+
+ Parameters
+ ----------
+ threshold : float, optional
+ A pixel only gets masked, if it lies above this threshold . The default is 0.
+
+ Returns
+ -------
+ mask : :class:`~numpy.ndarray` (N=3, boolean)
+ Mask of the pixels where the amplitude of the magnetization lies above `threshold`.
+
+ '''
+ return self.mag_amp > threshold
+
+ def get_vector(self, mask):
+ '''Returns the magnetic components arranged in a vector, specified by a mask.
+
+ Parameters
+ ----------
+ mask : :class:`~numpy.ndarray` (N=3, boolean)
+ Masks the pixels from which the components should be taken.
+
+ Returns
+ -------
+ vector : :class:`~numpy.ndarray` (N=1)
+ The vector containing magnetization components of the specified pixels.
+ Order is: first all `x`-, then all `y`-, then all `z`-components.
+
+ '''
+ if mask is not None:
+ return np.reshape([self.magnitude[0][mask],
+ self.magnitude[1][mask],
+ self.magnitude[2][mask]], -1)
+ else:
+ return self.mag_vec
+
+ def set_vector(self, vector, mask=None):
+ '''Set the magnetic components of the masked pixels to the values specified by `vector`.
+
+ Parameters
+ ----------
+ mask : :class:`~numpy.ndarray` (N=3, boolean), optional
+ Masks the pixels from which the components should be taken.
+ vector : :class:`~numpy.ndarray` (N=1)
+ The vector containing magnetization components of the specified pixels.
+ Order is: first all `x`-, then all `y-, then all `z`-components.
+
+ Returns
+ -------
+ None
+
+ '''
+ vector = np.asarray(vector, dtype=fft.FLOAT)
+ assert np.size(vector) % 3 == 0, 'Vector has to contain all 3 components for every pixel!'
+ count = np.size(vector)/3
+ if mask is not None:
+ self.magnitude[0][mask] = vector[:count] # x-component
+ self.magnitude[1][mask] = vector[count:2*count] # y-component
+ self.magnitude[2][mask] = vector[2*count:] # z-component
+ else:
+ self.mag_vec = vector
+
+ def save_to_llg(self, filename='..\output\magdata_output.txt'):
+ '''Save magnetization data in a file with LLG-format.
+
+ Parameters
+ ----------
+ filename : string, optional
+ The name of the LLG-file in which to store the magnetization data.
+ The default is '..\output\magdata_output.txt'.
+
+ Returns
+ -------
+ None
+
+ '''
+ self._log.debug('Calling save_to_llg')
+ a = self.a * 1.0E-9 / 1.0E-2 # from nm to cm
+ # Create 3D meshgrid and reshape it and the magnetization into a list where x varies first:
+ zz, yy, xx = (np.indices(self.dim)-a/2).reshape(3, -1)
+ x_vec, y_vec, z_vec = self.magnitude.reshape(3, -1)
+ # Save data to file:
+ data = np.array([xx, yy, zz, x_vec, y_vec, z_vec]).T
+ with open(filename, 'w') as mag_file:
+ mag_file.write('LLGFileCreator: %s\n' % filename.replace('.txt', ''))
+ mag_file.write(' %d %d %d\n' % (self.dim[2], self.dim[1], self.dim[0]))
+ mag_file.writelines('\n'.join(' '.join('{:7.6e}'.format(cell)
+ for cell in row) for row in data))
+
+ @classmethod
+ def load_from_llg(cls, filename):
+ '''Construct :class:`~.MagData` object from LLG-file.
+
+ Parameters
+ ----------
+ filename : string
+ The name of the LLG-file from which to load the data.
+
+ Returns
+ -------
+ mag_data: :class:`~.MagData`
+ A :class:`~.MagData` object containing the loaded data.
+
+ '''
+ cls._log.debug('Calling load_from_llg')
+ SCALE = 1.0E-9 / 1.0E-2 # From cm to nm
+ data = np.genfromtxt(filename, skip_header=2)
+ dim = tuple(np.genfromtxt(filename, dtype=int, skip_header=1, skip_footer=len(data[:, 0])))
+ a = (data[1, 0] - data[0, 0]) / SCALE
+ magnitude = data[:, 3:6].T.reshape((3,)+dim)
+ return MagData(a, magnitude)
+
+ def save_to_netcdf4(self, filename='..\output\magdata_output.nc'):
+ '''Save magnetization data in a file with NetCDF4-format.
+
+ Parameters
+ ----------
+ filename : string, optional
+ The name of the NetCDF4-file in which to store the magnetization data.
+ The default is '..\output\magdata_output.nc'.
+
+ Returns
+ -------
+ None
+
+ '''
+ self._log.debug('Calling save_to_netcdf4')
+ mag_file = netCDF4.Dataset(filename, 'w', format='NETCDF4')
+ mag_file.a = self.a
+ mag_file.createDimension('comp', 3) # Number of components
+ mag_file.createDimension('z_dim', self.dim[0])
+ mag_file.createDimension('y_dim', self.dim[1])
+ mag_file.createDimension('x_dim', self.dim[2])
+ magnitude = mag_file.createVariable('magnitude', 'f', ('comp', 'z_dim', 'y_dim', 'x_dim'))
+ magnitude[...] = self.magnitude
+ mag_file.close()
+
+ @classmethod
+ def load_from_netcdf4(cls, filename):
+ '''Construct :class:`~.DataMag` object from NetCDF4-file.
+
+ Parameters
+ ----------
+ filename : string
+ The name of the NetCDF4-file from which to load the data. Standard format is '\*.nc'.
+
+ Returns
+ -------
+ mag_data: :class:`~.MagData`
+ A :class:`~.MagData` object containing the loaded data.
+
+ '''
+ cls._log.debug('Calling copy')
+ mag_file = netCDF4.Dataset(filename, 'r', format='NETCDF4')
+ a = mag_file.a
+ magnitude = mag_file.variables['magnitude'][...]
+ mag_file.close()
+ return MagData(a, magnitude)
+
+ def save_to_x3d(self, filename='..\..\output\magdata_output.x3d', maximum=1):
+ '''Output the magnetization in the .x3d format for the Fraunhofer InstantReality Player.
+
+ Parameters
+ ----------
+ None
+
+ Returns
+ -------
+ None
+
+ '''
+ self._log.debug('Calling save_to_x3d')
+ from lxml import etree
+
+ dim = self.dim
+ # Create points and vector components as lists:
+ zz, yy, xx = (np.indices(dim)-0.5).reshape(3, -1)
+ x_mag = np.reshape(self.magnitude[0], (-1))
+ y_mag = np.reshape(self.magnitude[1], (-1))
+ z_mag = np.reshape(self.magnitude[2], (-1))
+ # Load template, load tree and write viewpoint information:
+ template = os.path.join(os.path.dirname(os.path.realpath(__file__)), 'template.x3d')
+ parser = etree.XMLParser(remove_blank_text=True)
+ tree = etree.parse(template, parser)
+ scene = tree.find('Scene')
+ etree.SubElement(scene, 'Viewpoint', position='0 0 {}'.format(1.5*dim[0]),
+ fieldOfView='1')
+ # Write each "spin"-tag separately:
+ for i in range(np.prod(dim)):
+ mag = np.sqrt(x_mag[i]**2+y_mag[i]**2+z_mag[i]**2)
+ if mag != 0:
+ spin_position = (xx[i]-dim[2]/2., yy[i]-dim[1]/2., zz[i]-dim[0]/2.)
+ sx_ref = 0
+ sy_ref = 1
+ sz_ref = 0
+ rot_x = sy_ref*z_mag[i] - sz_ref*y_mag[i]
+ rot_y = sz_ref*x_mag[i] - sx_ref*z_mag[i]
+ rot_z = sx_ref*y_mag[i] - sy_ref*x_mag[i]
+ angle = np.arccos(y_mag[i]/mag)
+ if norm((rot_x, rot_y, rot_z)) < 1E-10:
+ rot_x, rot_y, rot_z = 1, 0, 0
+ spin_rotation = (rot_x, rot_y, rot_z, angle)
+ spin_color = cmx.RdYlGn(mag/maximum)[:3]
+ spin_scale = (1., 1., 1.)
+ spin = etree.SubElement(scene, 'ProtoInstance',
+ DEF='Spin {}'.format(i), name='Spin_Proto')
+ etree.SubElement(spin, 'fieldValue', name='spin_position',
+ value='{} {} {}'.format(*spin_position))
+ etree.SubElement(spin, 'fieldValue', name='spin_rotation',
+ value='{} {} {} {}'.format(*spin_rotation))
+ etree.SubElement(spin, 'fieldValue', name='spin_color',
+ value='{} {} {}'.format(*spin_color))
+ etree.SubElement(spin, 'fieldValue', name='spin_scale',
+ value='{} {} {}'.format(*spin_scale))
+ # Write the tree into the file in pretty print format:
+ tree.write(filename, pretty_print=True)
+
+ def quiver_plot(self, title='Magnetization Distribution', axis=None, proj_axis='z',
+ ar_dens=1, ax_slice=None, log=False, scaled=True):
+ '''Plot a slice of the magnetization as a quiver plot.
+
+ Parameters
+ ----------
+ title : string, optional
+ The title for the plot.
+ axis : :class:`~matplotlib.axes.AxesSubplot`, optional
+ Axis on which the graph is plotted. Creates a new figure if none is specified.
+ proj_axis : {'z', 'y', 'x'}, optional
+ The axis, from which a slice is plotted. The default is 'z'.
+ ar_dens: int (optional)
+ Number defining the arrow density which is plotted. A higher ar_dens number skips more
+ arrows (a number of 2 plots every second arrow). Default is 1.
+ ax_slice : int, optional
+ The slice-index of the axis specified in `proj_axis`. Is set to the center of
+ `proj_axis` if not specified.
+ log : boolean, optional
+ Takes the Default is False.
+ scaled : boolean, optional
+ Normalizes the plotted arrows in respect to the highest one. Default is True.
+
+ Returns
+ -------
+ axis: :class:`~matplotlib.axes.AxesSubplot`
+ The axis on which the graph is plotted.
+
+ '''
+ self._log.debug('Calling quiver_plot')
+ assert proj_axis == 'z' or proj_axis == 'y' or proj_axis == 'x', \
+ 'Axis has to be x, y or z (as string).'
+ if proj_axis == 'z': # Slice of the xy-plane with z = ax_slice
+ self._log.debug('proj_axis == z')
+ if ax_slice is None:
+ self._log.debug('ax_slice is None')
+ ax_slice = int(self.dim[0]/2.)
+ plot_u = np.copy(self.magnitude[0][ax_slice, ...]) # x-component
+ plot_v = np.copy(self.magnitude[1][ax_slice, ...]) # y-component
+ u_label = 'x [px]'
+ v_label = 'y [px]'
+ elif proj_axis == 'y': # Slice of the xz-plane with y = ax_slice
+ self._log.debug('proj_axis == y')
+ if ax_slice is None:
+ self._log.debug('ax_slice is None')
+ ax_slice = int(self.dim[1]/2.)
+ plot_u = np.copy(self.magnitude[0][:, ax_slice, :]) # x-component
+ plot_v = np.copy(self.magnitude[2][:, ax_slice, :]) # z-component
+ u_label = 'x [px]'
+ v_label = 'z [px]'
+ elif proj_axis == 'x': # Slice of the yz-plane with x = ax_slice
+ self._log.debug('proj_axis == x')
+ if ax_slice is None:
+ self._log.debug('ax_slice is None')
+ ax_slice = int(self.dim[2]/2.)
+ plot_u = np.swapaxes(np.copy(self.magnitude[2][..., ax_slice]), 0, 1) # z-component
+ plot_v = np.swapaxes(np.copy(self.magnitude[1][..., ax_slice]), 0, 1) # y-component
+ u_label = 'z [px]'
+ v_label = 'y [px]'
+ # If no axis is specified, a new figure is created:
+ if axis is None:
+ self._log.debug('axis is None')
+ fig = plt.figure(figsize=(8.5, 7))
+ axis = fig.add_subplot(1, 1, 1)
+ axis.set_aspect('equal')
+ angles = np.angle(plot_u+1j*plot_v, deg=True)
+ # Take the logarithm of the arrows to clearly show directions (if specified):
+ if log:
+ cutoff = 10
+ amp = np.round(np.hypot(plot_u, plot_v), decimals=cutoff)
+ min_value = amp[np.nonzero(amp)].min()
+ plot_u = np.round(plot_u, decimals=cutoff) / min_value
+ plot_u = np.log10(np.abs(plot_u)+1) * np.sign(plot_u)
+ plot_v = np.round(plot_v, decimals=cutoff) / min_value
+ plot_v = np.log10(np.abs(plot_v)+1) * np.sign(plot_v)
+ # Scale the magnitude of the arrows to the highest one (if specified):
+ if scaled:
+ plot_u /= np.hypot(plot_u, plot_v).max()
+ plot_v /= np.hypot(plot_u, plot_v).max()
+ # Setup quiver:
+ dim_uv = plot_u.shape
+ ad = ar_dens
+ yy, xx = np.indices(dim_uv)
+ axis.quiver(xx[::ad, ::ad], yy[::ad, ::ad], plot_u[::ad, ::ad], plot_v[::ad, ::ad],
+ pivot='middle', units='xy', angles=angles[::ad, ::ad], minlength=0.25,
+ scale_units='xy', scale=1./ad, headwidth=6, headlength=7)
+ axis.set_xlim(-1, dim_uv[1])
+ axis.set_ylim(-1, dim_uv[0])
+ axis.set_title(title, fontsize=18)
+ axis.set_xlabel(u_label, fontsize=15)
+ axis.set_ylabel(v_label, fontsize=15)
+ axis.tick_params(axis='both', which='major', labelsize=14)
+ axis.xaxis.set_major_locator(MaxNLocator(nbins=9, integer=True))
+ axis.yaxis.set_major_locator(MaxNLocator(nbins=9, integer=True))
+ return axis
+
+ def quiver_plot3d(self, title='Magnetization Distribution', limit=None, cmap='cool',
+ ar_dens=1, mode='arrow', show_pipeline=False):
+ '''Plot the magnetization as 3D-vectors in a quiverplot.
+
+ Parameters
+ ----------
+ title : string, optional
+ The title for the plot.
+ limit : float, optional
+ Plotlimit for the magnetization arrow length used to scale the colormap.
+ cmap : string, optional
+ String describing the colormap which is used (default is 'cool').
+ ar_dens: int (optional)
+ Number defining the arrow density which is plotted. A higher ar_dens number skips more
+ arrows (a number of 2 plots every second arrow). Default is 1.
+ mode: string, optional
+ Mode, determining the glyphs used in the 3D plot. Default is 'arrow', which corresponds
+ to 3D arrows. For large amounts of arrows, '2darrow' should be used.
+ show_pipeline : boolean, optional
+ If True, the mayavi pipeline, a GUI used for image manipulation is shown. The default
+ is False.
+
+ Returns
+ -------
+ plot : :class:`mayavi.modules.vectors.Vectors`
+ The plot object.
+
+ '''
+ self._log.debug('Calling quiver_plot3D')
+ from mayavi import mlab
+ a = self.a
+ dim = self.dim
+ if limit is None:
+ limit = np.max(self.mag_amp)
+ ad = ar_dens
+ # Create points and vector components as lists:
+ zz, yy, xx = (np.indices(dim)-a/2).reshape(3, -1)
+ x_mag = self.magnitude[0][::ad, ::ad, ::ad].flatten()
+ y_mag = self.magnitude[1][::ad, ::ad, ::ad].flatten()
+ z_mag = self.magnitude[2][::ad, ::ad, ::ad].flatten()
+ # Plot them as vectors:
+ mlab.figure()
+ plot = mlab.quiver3d(xx, yy, zz, x_mag, y_mag, z_mag, mode=mode, colormap=cmap)
+ plot.glyph.glyph_source.glyph_position = 'center'
+ plot.module_manager.vector_lut_manager.data_range = np.array([0, limit])
+ mlab.outline(plot)
+ mlab.axes(plot)
+ mlab.title(title, height=0.95, size=0.35)
+ mlab.colorbar(label_fmt='%.2f')
+ mlab.colorbar(orientation='vertical')
+ mlab.orientation_axes()
+ if show_pipeline:
+ mlab.show_pipeline()
+ return plot
diff --git a/pyramid/numcore/phasemapper_core.pyx b/pyramid/numcore/phasemapper_core.pyx
index d7b1f049c8405a0c2719c78d5b00fc3a5c7e7a81..59be7c8ff4a84efec61e432707b54c81ae58c6f5 100644
--- a/pyramid/numcore/phasemapper_core.pyx
+++ b/pyramid/numcore/phasemapper_core.pyx
@@ -106,7 +106,6 @@ def jac_dot_real_convolve(
result[r] -= vector[s+size] * v_phi[v, u]
r += 1
s += 1
- # TODO: linearize u and v, too?
@cython.boundscheck(False)
@@ -155,4 +154,3 @@ def jac_T_dot_real_convolve(
result[r+size] -= vector[s] * v_phi[v, u]
s += 1
r += 1
- # TODO: linearize u and v, too?
diff --git a/pyramid/phasemap.py b/pyramid/phasemap.py
index d60b4105e2636b71cd10217b6f48120194d15181..ce5c45adc79b678b83c066b5e3de7bb9344b665e 100644
--- a/pyramid/phasemap.py
+++ b/pyramid/phasemap.py
@@ -1,4 +1,7 @@
# -*- coding: utf-8 -*-
+# Copyright 2014 by Forschungszentrum Juelich GmbH
+# Author: J. Caron
+#
"""This module provides the :class:`~.PhaseMap` class for storing phase map data."""
@@ -323,8 +326,8 @@ class PhaseMap(object):
Returns
-------
- axis: :class:`~matplotlib.axes.AxesSubplot`
- The axis on which the graph is plotted.
+ axis, cbar: :class:`~matplotlib.axes.AxesSubplot`
+ The axis on which the graph is plotted and the colorbar.
'''
self._log.debug('Calling display_phase')
@@ -359,7 +362,7 @@ class PhaseMap(object):
cbar.ax.tick_params(labelsize=14)
cbar.set_label(u'phase shift [{}]'.format(self.unit), fontsize=15)
# Return plotting axis:
- return axis
+ return axis, cbar
def display_phase3d(self, title='Phase Map', cmap='RdBu'):
'''Display the phasemap as a 3-D surface with contourplots.
@@ -385,9 +388,7 @@ class PhaseMap(object):
fig = plt.figure()
axis = Axes3D(fig)
# Plot surface and contours:
- u = range(self.dim_uv[1])
- v = range(self.dim_uv[0])
- uu, vv = np.meshgrid(u, v)
+ vv, uu = np.indices(self.dim_uv)
axis.plot_surface(uu, vv, phase, rstride=4, cstride=4, alpha=0.7, cmap=cmap,
linewidth=0, antialiased=False)
axis.contourf(uu, vv, phase, 15, zdir='z', offset=np.min(phase), cmap=cmap)
@@ -547,9 +548,7 @@ class PhaseMap(object):
'''
cls._log.debug('Calling make_color_wheel')
- x = np.linspace(-256, 256, num=512)
- y = np.linspace(-256, 256, num=512)
- xx, yy = np.meshgrid(x, y)
+ yy, xx = np.indices((512, 512)) - 256
r = np.sqrt(xx ** 2 + yy ** 2)
# Create the wheel:
color_wheel_magnitude = (1 - np.cos(r * pi/360)) / 2
diff --git a/pyramid/phasemapper.py b/pyramid/phasemapper.py
index 76feddabac165232ad715912e0076249384515fa..b4c6cab2318fc121de2ab4d1f937cf3b66bba6fe 100644
--- a/pyramid/phasemapper.py
+++ b/pyramid/phasemapper.py
@@ -1,4 +1,7 @@
# -*- coding: utf-8 -*-
+# Copyright 2014 by Forschungszentrum Juelich GmbH
+# Author: J. Caron
+#
"""This module executes several forward models to calculate the magnetic or electric phase map from
a given projection of a 3-dimensional magnetic distribution (see :mod:`~pyramid.projector`).
For the magnetic phase map, an approach using real space and one using Fourier space is provided.
@@ -431,12 +434,10 @@ class PhaseMapperFDFC(PhaseMapper):
mag_proj = MagData(self.a, np.zeros((3, 1)+self.dim_uv))
magnitude_proj = self.jac_dot(vector).reshape((2, )+self.dim_uv)
mag_proj.magnitude[1:3, 0, ...] = magnitude_proj
- # TODO: instead call common subroutine operating on u_mag, v_mag with __call__?
return self(mag_proj).phase_vec
def jac_T_dot(self, vector):
raise NotImplementedError()
- # TODO: Implement!
class PhaseMapperElectric(PhaseMapper):
@@ -509,11 +510,9 @@ class PhaseMapperElectric(PhaseMapper):
def jac_dot(self, vector):
raise NotImplementedError()
- # TODO: Implement?
def jac_T_dot(self, vector):
raise NotImplementedError()
- # TODO: Implement?
def pm(mag_data, axis='z', dim_uv=None, b_0=1):
diff --git a/pyramid/projector.py b/pyramid/projector.py
index 48e69443f4e021b765419c658f57b2f95739afc2..773db4bb8531c8d1f9ab7d9fc4e5690c706f1545 100644
--- a/pyramid/projector.py
+++ b/pyramid/projector.py
@@ -1,4 +1,7 @@
# -*- coding: utf-8 -*-
+# Copyright 2014 by Forschungszentrum Juelich GmbH
+# Author: J. Caron
+#
"""This module provides the abstract base class :class:`~.Projector` and concrete subclasses for
projections of vector and scalar fields."""
@@ -436,7 +439,7 @@ class SimpleProjector(Projector):
self._log.debug('Calling __init__')
assert axis in {'z', 'y', 'x'}, 'Projection axis has to be x, y or z (given as a string)!'
proj, v, u = self.AXIS_DICT[axis]
- if axis=='x':
+ if axis == 'x':
dim_proj, dim_v, dim_u = dim[proj], dim[u], dim[v] # coordinate switch for 'x'!
else:
dim_proj, dim_v, dim_u = dim[proj], dim[v], dim[u]
@@ -453,14 +456,14 @@ class SimpleProjector(Projector):
elif axis == 'y':
self._log.debug('Projection along the y-axis')
coeff = [[1, 0, 0], [0, 0, 1]]
- indices = np.array([np.arange(row%dim_x, dim_x*dim_y, dim_x)+int(row/dim_x)*dim_x*dim_y
+ indices = np.array([np.arange(row % dim_x, dim_x*dim_y, dim_x) + row//dim_x*dim_x*dim_y
for row in range(size_2d)]).reshape(-1)
elif axis == 'x':
self._log.debug('Projection along the x-axis')
coeff = [[0, 0, 1], [0, 1, 0]] # Caution, coordinate switch: u, v --> z, y (not y, z!)
# indices = np.array([np.arange(dim_proj) + row*dim_proj
# for row in range(size_2d)]).reshape(-1) # this is u, v --> y, z
- indices = np.array([np.arange(dim_x) + (row%dim_z)*dim_x*dim_y + int(row/dim_z)*dim_x
+ indices = np.array([np.arange(dim_x) + (row % dim_z)*dim_x*dim_y + row//dim_z*dim_x
for row in range(size_2d)]).reshape(-1)
if dim_uv is not None:
indptr = indptr.tolist() # convert to use insert() and append()
@@ -495,3 +498,5 @@ class SimpleProjector(Projector):
'''
return 'projected along {}-axis'.format(self.axis)
+
+# TODO: Arbitrary Projections!
diff --git a/pyramid/reconstruction.py b/pyramid/reconstruction.py
index d9386eeb90a2285538b37aed1dbedebf3d15e5e6..f5df52a5f6325f35d096be2f4f180f1594d1f150 100644
--- a/pyramid/reconstruction.py
+++ b/pyramid/reconstruction.py
@@ -1,4 +1,7 @@
# -*- coding: utf-8 -*-
+# Copyright 2014 by Forschungszentrum Juelich GmbH
+# Author: J. Caron
+#
"""Reconstruct magnetic distributions from given phasemaps.
This module reconstructs 3-dimensional magnetic distributions (as :class:`~pyramid.magdata.MagData`
@@ -12,77 +15,16 @@ the distribution.
import numpy as np
-from pyramid.kernel import Kernel
-from pyramid.projector import SimpleProjector
-from pyramid.phasemapper import PhaseMapperRDFC
from pyramid.costfunction import Costfunction
from pyramid.magdata import MagData
-from scipy.optimize import leastsq
-
import logging
-__all__ = ['optimize_linear', 'optimize_nonlin', 'optimize_splitbregman',
- 'optimize_simple_leastsq']
+__all__ = ['optimize_linear', 'optimize_nonlin', 'optimize_splitbregman']
_log = logging.getLogger(__name__)
-class PrintIterator(object):
-
- '''Iterator class which is responsible to give feedback during reconstruction iterations.
-
- Parameters
- ----------
- cost : :class:`~.Costfunction`
- :class:`~.Costfunction` class for outputting the `cost` of the current magnetization
- distribution. This should decrease per iteration if the algorithm converges and is only
- printed for a `verbosity` of 2.
- verbosity : {0, 1, 2}, optional
- Parameter defining the verbosity of the output. `2` will show the current number of the
- iteration and the cost of the current distribution. `1` will just show the iteration
- number and `0` will prevent output all together.
-
- Notes
- -----
- Normally this class should not be used by the user and is instantiated whithin the
- :mod:`~.reconstruction` module itself.
-
- '''
-
- _log = logging.getLogger(__name__ + '.PrintIterator')
-
- def __init__(self, cost, verbosity):
- self._log.debug('Calling __init__')
- self.cost = cost
- self.verbosity = verbosity
- assert verbosity in {0, 1, 2}, 'verbosity has to be set to 0, 1 or 2!'
- self.iteration = 0
- self._log.debug('Created ' + str(self))
-
- def __call__(self, xk):
- self._log.debug('Calling __call__')
- if self.verbosity == 0:
- return
- print 'iteration #', self.next(),
- if self.verbosity > 1:
- print 'cost =', self.cost(xk)
- else:
- print ''
-
- def __repr__(self):
- self._log.debug('Calling __repr__')
- return '%s(cost=%r, verbosity=%r)' % (self.__class__, self.cost, self.verbosity)
-
- def __str__(self):
- self._log.debug('Calling __str__')
- return 'PrintIterator(cost=%s, verbosity=%s)' % (self.cost, self.verbosity)
-
- def next(self):
- self.iteration += 1
- return self.iteration
-
-
def optimize_linear(data, regularisator=None, max_iter=None):
'''Reconstruct a three-dimensional magnetic distribution from given phase maps via the
conjugate gradient optimizaion method :func:`~.scipy.sparse.linalg.cg`.
@@ -103,7 +45,7 @@ def optimize_linear(data, regularisator=None, max_iter=None):
mag_data : :class:`~pyramid.magdata.MagData`
The reconstructed magnetic distribution as a :class:`~.MagData` object.
- ''' # TODO: info document!
+ '''
import jutil.cg as jcg
_log.debug('Calling optimize_linear')
# Set up necessary objects:
@@ -128,7 +70,7 @@ def optimize_nonlin(data, first_guess=None, regularisator=None):
:class:`~.DataSet` object containing all phase maps in :class:`~.PhaseMap` objects and all
projection directions in :class:`~.Projector` objects. These provide the essential
information for the reconstruction.
- first_fuess : :class:`~pyramid.magdata.MagData`
+ first_guess : :class:`~pyramid.magdata.MagData`
magnetization to start the non-linear iteration with.
regularisator : :class:`~.Regularisator`, optional
Regularisator class that's responsible for the regularisation term.
@@ -151,7 +93,6 @@ def optimize_nonlin(data, first_guess=None, regularisator=None):
p = regularisator.p
q = 1. / (1. - (1. / p))
- lp = regularisator.norm
lq = jnorms.LPPow(q, 1e-20)
def preconditioner(_, direc):
@@ -172,6 +113,7 @@ def optimize_nonlin(data, first_guess=None, regularisator=None):
mag_opt.set_vector(x_opt, data.mask)
return mag_opt
+
def optimize_splitbregman(data, weight, lam, mu):
'''
Reconstructs magnet distribution from phase image measurements using a split bregman
@@ -209,7 +151,6 @@ def optimize_splitbregman(data, weight, lam, mu):
# function to that which is supposedly optimized by split bregman.
# Thus cost can be used to verify convergence
regularisator = FirstOrderRegularisator(data.mask, lam / mu, 1)
- x_0 = MagData(data.a, np.zeros((3,) + data.dim)).get_vector(data.mask)
cost = Costfunction(data, regularisator)
fwd_mod = cost.fwd_model
@@ -229,80 +170,3 @@ def optimize_splitbregman(data, weight, lam, mu):
mag_opt = MagData(data.a, np.zeros((3,) + data.dim))
mag_opt.set_vector(x_opt, data.mask)
return mag_opt
-
-
-def optimize_simple_leastsq(phase_map, mask, b_0=1, lam=1E-4, order=0):
- '''Reconstruct a magnetic distribution for a 2-D problem with known pixel locations.
-
- Parameters
- ----------
- phase_map : :class:`~pyramid.phasemap.PhaseMap`
- A :class:`~pyramid.phasemap.PhaseMap` object, representing the phase from which to
- reconstruct the magnetic distribution.
- mask : :class:`~numpy.ndarray` (N=3)
- A boolean matrix (or a matrix consisting of ones and zeros), representing the
- positions of the magnetized voxels in 3 dimensions.
- b_0 : float, optional
- The magnetic induction corresponding to a magnetization `M`\ :sub:`0` in T.
- The default is 1.
- lam : float, optional
- The regularisation parameter. Defaults to 1E-4.
- order : int {0, 1}, optional
- order of the regularisation function. Default is 0 for a Tikhonov regularisation of order
- zero. A first order regularisation, which uses the derivative is available with value 1.
-
- Returns
- -------
- mag_data : :class:`~pyramid.magdata.MagData`
- The reconstructed magnetic distribution as a :class:`~.MagData` object.
-
- Notes
- -----
- Only works for a single phase_map, if the positions of the magnetized voxels are known and
- for slice thickness of 1 (constraint for the `z`-dimension).
-
- '''
- # Read in parameters:
- y_m = phase_map.phase_vec # Measured phase map as a vector
- a = phase_map.a # Grid spacing
- dim = mask.shape # Dimensions of the mag. distr.
- count = mask.sum() # Number of pixels with magnetization
- # Create empty MagData object for the reconstruction:
- mag_data_rec = MagData(a, np.zeros((3,) + dim))
-
- # Function that returns the phase map for a magnetic configuration x:
- def F(x):
- mag_data_rec.set_vector(x, mask)
- proj = SimpleProjector(dim)
- phase_map = PhaseMapperRDFC(Kernel(a, proj.dim_uv, b_0))(proj(mag_data_rec))
- return phase_map.phase_vec
-
- # Cost function of order zero which should be minimized:
- def J_0(x_i):
- y_i = F(x_i)
- term1 = (y_i - y_m)
- term2 = lam * x_i
- return np.concatenate([term1, term2])
-
- # First order cost function which should be minimized:
- def J_1(x_i):
- y_i = F(x_i)
- term1 = (y_i - y_m)
- mag_data = mag_data_rec.magnitude
- term2 = []
- for i in range(3):
- component = mag_data[i, ...]
- for j in range(3):
- if component.shape[j] > 1:
- term2.append(np.diff(component, axis=j).reshape(-1))
-
- term2 = lam * np.concatenate(term2)
- return np.concatenate([term1, term2])
-
- J_DICT = [J_0, J_1] # list of cost-functions with different regularisations
- # Reconstruct the magnetization components:
- # TODO Use jutil.minimizer.minimize(jutil.costfunction.LeastSquaresCostFunction(J_DICT[order],
- # ...) or a simpler frontend.
- x_rec, _ = leastsq(J_DICT[order], np.zeros(3 * count))
- mag_data_rec.set_vector(x_rec, mask)
- return mag_data_rec
diff --git a/pyramid/regularisator.py b/pyramid/regularisator.py
index 5c115bb98430185f64e3fbb8c52d34796b673d9a..5f02f240e0153ee9de68cd59636948566d31eaa9 100644
--- a/pyramid/regularisator.py
+++ b/pyramid/regularisator.py
@@ -1,9 +1,9 @@
# -*- coding: utf-8 -*-
-"""
-Created on Mon Aug 18 09:24:58 2014
-
-@author: Jan
-""" # TODO: Docstring!
+# Copyright 2014 by Forschungszentrum Juelich GmbH
+# Author: J. Caron
+#
+"""This module provides the :class:`~.Regularisator` class which represents a regularisation term
+which adds additional constraints to a costfunction to minimize."""
import abc
@@ -19,14 +19,23 @@ import logging
__all__ = ['NoneRegularisator', 'ZeroOrderRegularisator', 'FirstOrderRegularisator']
-# TODO: Fragen für Jörn: Macht es Sinn, x_a standardmäßig auf den Nullvektor zu setzen? Ansonsten
-# besser im jeweiligen Konstruktor setzen, nicht im abstrakten!
-# Wie kommt man genau an die Ableitungen (Normen sind nicht unproblematisch)?
-# Wie implementiert man am besten verschiedene Normen?
-
class Regularisator(object):
- # TODO: Docstring!
+ '''Class for providing a regularisation term which implements additional constraints.
+
+ Represents a certain constraint for the 3D magnetization distribution whose cost is to minimize
+ in addition to the derivation from the 2D phase maps. Important is the used `norm` and the
+ regularisation parameter `lam` (lambda) which determines the weighting between the two cost
+ parts (measurements and regularisation).
+
+ Attributes
+ ----------
+ norm: :class:`~jutil.norm.WeightedNorm`
+ Norm, which is used to determine the cost of the regularisation term.
+ lam: float
+ Regularisation parameter determining the weighting between measurements and regularisation.
+
+ '''
__metaclass__ = abc.ABCMeta
_log = logging.getLogger(__name__+'.Regularisator')
@@ -51,25 +60,71 @@ class Regularisator(object):
return 'Regularisator(norm=%s, lam=%s)' % (self.norm, self.lam)
def jac(self, x):
- # TODO: Docstring!
+ '''Calculate the derivative of the regularisation term for a given magnetic 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.
+
+ '''
return self.lam * self.norm.jac(x)
def hess_dot(self, x, vector):
- # TODO: Docstring!
+ '''Calculate the product of a `vector` with the Hessian matrix of the regularisation term.
+
+ 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 self.lam * self.norm.hess_dot(x, vector)
def hess_diag(self, x, vector):
- # TODO: Docstring!
+ ''' Return the diagonal of the Hessian.
+
+ Parameters
+ ----------
+ _ : undefined
+ Unused input
+
+ '''
self._log.debug('Calling hess_diag')
return self.lam * self.norm.hess_diag(x, vector)
class NoneRegularisator(Regularisator):
- # TODO: Docstring
+ '''Placeholder class if no regularization is used.
- # TODO: Necessary class? Use others with lam=0?
+ This class is instantiated in the :class:`~pyramid.costfunction.Costfunction`, which means
+ no regularisation is used. All associated functions return appropriate zero-values.
- LOG = logging.getLogger(__name__+'.NoneRegularisator')
+ Attributes
+ ----------
+ norm: None
+ No regularization is used, thus also no norm.
+ lam: 0
+ Not used.
+
+ '''
+
+ _log = logging.getLogger(__name__+'.NoneRegularisator')
def __init__(self):
self._log.debug('Calling __init__')
@@ -82,23 +137,72 @@ class NoneRegularisator(Regularisator):
return 0
def jac(self, x):
- # TODO: Docstring!
+ '''Calculate the derivative of the regularisation term for a given magnetic 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.
+
+ '''
return np.zeros_like(x)
def hess_dot(self, x, vector):
- # TODO: Docstring!
+ '''Calculate the product of a `vector` with the Hessian matrix of the regularisation term.
+
+ 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 np.zeros_like(vector)
def hess_diag(self, x, vector):
- # TODO: Docstring!
+ ''' Return the diagonal of the Hessian.
+
+ Parameters
+ ----------
+ _ : undefined
+ Unused input
+
+ '''
self._log.debug('Calling hess_diag')
return np.zeros_like(vector)
class ZeroOrderRegularisator(Regularisator):
- # TODO: Docstring!
+ '''Class for providing a regularisation term which implements Lp norm minimization.
+
+ The constraint this class represents is the minimization of the Lp norm for the 3D
+ magnetization distribution. Important is the regularisation parameter `lam` (lambda) which
+ determines the weighting between the two cost parts (measurements and regularisation).
+
+ Attributes
+ ----------
+ lam: float
+ Regularisation parameter determining the weighting between measurements and regularisation.
+ p: int, optional
+ Order of the norm (default: 2, which means a standard L2-norm).
+
+ '''
- LOG = logging.getLogger(__name__+'.ZeroOrderRegularisator')
+ _log = logging.getLogger(__name__+'.ZeroOrderRegularisator')
def __init__(self, _, lam, p=2):
self._log.debug('Calling __init__')
@@ -112,7 +216,23 @@ class ZeroOrderRegularisator(Regularisator):
class FirstOrderRegularisator(Regularisator):
- # TODO: Docstring!
+ '''Class for providing a regularisation term which implements derivation minimization.
+
+ The constraint this class represents is the minimization of the first order derivative of the
+ 3D magnetization distribution using a Lp norm. Important is the regularisation parameter `lam`
+ (lambda) which determines the weighting between the two cost parts (measurements and
+ regularisation).
+
+ Attributes
+ ----------
+ mask: :class:`~numpy.ndarray` (N=3)
+ A boolean mask which defines the magnetized volume in 3D.
+ lam: float
+ Regularisation parameter determining the weighting between measurements and regularisation.
+ p: int, optional
+ Order of the norm (default: 2, which means a standard L2-norm).
+
+ '''
def __init__(self, mask, lam, p=2):
self.p = p
diff --git a/pyramid/version.py b/pyramid/version.py
index e604916aeed68543a5cb5bef5572ca62b7d13d7a..81491db2e62cb18f16d93600d7fc02b5c547f5ea 100644
--- a/pyramid/version.py
+++ b/pyramid/version.py
@@ -1,3 +1,3 @@
# THIS FILE IS GENERATED FROM THE PYRAMID SETUP.PY
version = "0.1.0-dev"
-hg_revision = "9013fd38a1b6+"
+hg_revision = "51a00f1d3934+"
diff --git a/scripts/collaborations/patrick_nanowire_reconstruction.py b/scripts/collaborations/patrick_nanowire_reconstruction.py
index deb5d335ab965d6709aef20468929ca8d757fb2a..80c625bfc3480aea6cc3a0a2ad0b256ff9e8e1fe 100644
--- a/scripts/collaborations/patrick_nanowire_reconstruction.py
+++ b/scripts/collaborations/patrick_nanowire_reconstruction.py
@@ -20,26 +20,29 @@ LOGGING_CONF = os.path.join(os.path.dirname(os.path.realpath(pyramid.__file__)),
logging.config.fileConfig(LOGGING_CONF, disable_existing_loggers=False)
+# TODO: read in input via txt-file dictionary?
+
###################################################################################################
a = 1.455 # in nm
gain = 5
b_0 = 1
-lam = 1E-6
+lam = 1E-4
PATH = '../../output/patrick/'
-PHASE = '30_pha'
-MASK = 'mask30_elong_5'
-longFOV = True
+PHASE = 'pos3_40deg_magphase'
+MASK = 'pos3_40deg_maskbyhand'
+FORMAT = '.bmp'
+longFOV = False
longFOV_string = np.where(longFOV, 'longFOV', 'normalFOV')
IMAGENAME = '{}_{}_{}_'.format(MASK, PHASE, longFOV_string)
-PHASE_MAX = 7.68 # -10°: 0.92, 30°: 7.68
-PHASE_MIN = -18.89 # -10°: -16.85, 30°: -18.89
+PHASE_MAX = 0.5 # -10°: 0.92, 30°: 7.68
+PHASE_MIN = -0.5 # -10°: -16.85, 30°: -18.89
PHASE_DELTA = PHASE_MAX - PHASE_MIN
###################################################################################################
# Read in files:
-im_mask = Image.open(PATH+MASK+'.tif')
+im_mask = Image.open(PATH+MASK+FORMAT)
#im_mask.thumbnail((125, 175), Image.ANTIALIAS)
-im_phase = Image.open(PATH+PHASE+'.tif')
+im_phase = Image.open(PATH+PHASE+FORMAT)
#im_phase.thumbnail((125, 125), Image.ANTIALIAS)
mask = np.array(np.array(im_mask)/255, dtype=bool)
@@ -67,7 +70,7 @@ if longFOV:
regularisator = FirstOrderRegularisator(mask, lam, p=2)
with TakeTime('reconstruction time'):
- mag_data_rec = rc.optimize_linear(data_set, regularisator=regularisator, max_iter=50)
+ mag_data_rec = rc.optimize_linear(data_set, regularisator=regularisator, max_iter=1000)[0]
phase_map_rec_pad = pm(mag_data_rec)
phase_map_rec = PhaseMap(a, phase_map_rec_pad.phase[pad:, :])#[pad:-pad, pad:-pad])
@@ -83,6 +86,6 @@ phase_map_rec.display_combined('Reconstr. Distribution', gain=4)
plt.savefig(PATH+IMAGENAME+'RECONSTRUCTION.png')
phase_map_diff.display_combined('Difference')
plt.savefig(PATH+IMAGENAME+'DIFFERENCE.png')
-mag_data_rec.scale_down(4)
-mag_data_rec.quiver_plot(log=True)
+#mag_data_rec.scale_down(4)
+mag_data_rec.quiver_plot(log=True, ar_dens=8)
plt.savefig(PATH+IMAGENAME+'MAGNETIZATION_DOWNSCALE4.png')
diff --git a/scripts/collaborations/patrick_nanowire_reconstruction_simulation.py b/scripts/collaborations/patrick_nanowire_reconstruction_simulation.py
index 3c1edbac8cdeccc685662f1f4bb0f70a4cb928c9..81efbba528f98620f80ba2183eadd218ca06eb4a 100644
--- a/scripts/collaborations/patrick_nanowire_reconstruction_simulation.py
+++ b/scripts/collaborations/patrick_nanowire_reconstruction_simulation.py
@@ -28,7 +28,7 @@ length = 300
cutoff = 50
trans = 40
pads = [100, 250, 500, 800]
-INPATH = '../../output/vtk data\longtube_withcap/'
+INPATH = '../../output/vtk data/longtube_withcap/'
OUTPATH = '../../output/patrick/'
FILE = 'CoFeB_tube_cap_4nm_lying_down'
###################################################################################################
diff --git a/scripts/collaborations/rueffner_file.py b/scripts/collaborations/rueffner_file.py
index 230c4cb6805e91af2676294e7b71cedb8bcff270..82b13fbe72fc8cc97a46935301c1164a62d57ba5 100644
--- a/scripts/collaborations/rueffner_file.py
+++ b/scripts/collaborations/rueffner_file.py
@@ -11,6 +11,9 @@ from tqdm import tqdm
import tables
+import psutil
+import gc
+
import pyramid
from pyramid.magdata import MagData
from pyramid.projector import SimpleProjector
@@ -25,6 +28,7 @@ LOGGING_CONF = os.path.join(os.path.dirname(os.path.realpath(pyramid.__file__)),
logging.config.fileConfig(LOGGING_CONF, disable_existing_loggers=False)
+proc = psutil.Process(os.getpid())
###################################################################################################
PATH = '../../output/'
dim = (16, 190, 220)
@@ -68,7 +72,7 @@ zs = np.arange(-dim[0]/2, dim[0]/2)
xx, yy = np.meshgrid(xs, ys)
-def calculate(t): # TODO: Somehow avoid memory error :-(...
+def calculate(t):
print 't =', t
vectors = h5file.root.data.fields.m.read(field='m_CoFeb')[t, ...]
data = np.hstack((points, vectors))
@@ -103,3 +107,5 @@ def calculate(t): # TODO: Somehow avoid memory error :-(...
# Interpolation and phase calculation for all timesteps:
for t in np.arange(0, 1001, 5):
calculate(t)
+ gc.collect()
+ print 'RSS = {:.2f} MB'.format(proc.memory_info().rss/1024.**2)
diff --git a/scripts/collaborations/vtk_conversion.py b/scripts/collaborations/vtk_conversion.py
new file mode 100644
index 0000000000000000000000000000000000000000..62f21be3b257a39df1b2a491f43321ec74cc1387
--- /dev/null
+++ b/scripts/collaborations/vtk_conversion.py
@@ -0,0 +1,195 @@
+# -*- coding: utf-8 -*-
+"""Created on Fri Jan 24 11:17:11 2014 @author: Jan"""
+
+
+import os
+
+import numpy as np
+import matplotlib.pyplot as plt
+from pylab import griddata
+
+import pickle
+import vtk
+from tqdm import tqdm
+
+import pyramid
+from pyramid.magdata import MagData
+from pyramid.projector import XTiltProjector
+from pyramid.phasemapper import PhaseMapperRDFC
+from pyramid.kernel import Kernel
+
+import logging
+import logging.config
+
+
+LOGGING_CONF = os.path.join(os.path.dirname(os.path.realpath(pyramid.__file__)), 'logging.ini')
+
+
+logging.config.fileConfig(LOGGING_CONF, disable_existing_loggers=False)
+###################################################################################################
+PATH = '../../output/vtk data/tube_160x30x1100nm/02758'
+b_0 = 1.54
+gain = 8
+force_calculation = False
+###################################################################################################
+# Load vtk-data:
+if force_calculation or not os.path.exists(PATH+'.pickle'):
+ # Setting up reader:
+ reader = vtk.vtkDataSetReader()
+ reader.SetFileName(PATH+'.vtk')
+ reader.ReadAllScalarsOn()
+ reader.ReadAllVectorsOn()
+ reader.Update()
+ # Getting output:
+ output = reader.GetOutput()
+ # Reading points and vectors:
+ size = output.GetNumberOfPoints()
+ vtk_points = output.GetPoints().GetData()
+ vtk_vectors = output.GetPointData().GetVectors()
+ # Converting points to numpy array:
+ point_array = np.zeros(vtk_points.GetSize())
+ vtk_points.ExportToVoidPointer(point_array)
+ point_array = np.reshape(point_array, (-1, 3))
+ # Converting vectors to numpy array:
+ vector_array = np.zeros(vtk_points.GetSize())
+ vtk_vectors.ExportToVoidPointer(vector_array)
+ vector_array = np.reshape(vector_array, (-1, 3))
+ # Combining data:
+ data = np.hstack((point_array, vector_array))
+ with open(PATH+'.pickle', 'w') as pf:
+ pickle.dump(data, pf)
+else:
+ with open(PATH+'.pickle') as pf:
+ data = pickle.load(pf)
+# Scatter plot of all x-y-coordinates
+axis = plt.figure().add_subplot(1, 1, 1)
+axis.scatter(data[:, 0], data[:, 1])
+plt.show()
+###################################################################################################
+# Interpolate on regular grid:
+if force_calculation or not os.path.exists(PATH+'.nc'):
+ # Find unique z-slices:
+ zs = np.unique(data[:, 2])
+ # Determine the grid spacing:
+ a = zs[1] - zs[0]
+ # Determine the size of object:
+ x_min, x_max = data[:, 0].min(), data[:, 0].max()
+ y_min, y_max = data[:, 1].min(), data[:, 1].max()
+ z_min, z_max = data[:, 2].min(), data[:, 2].max()
+ x_diff, y_diff, z_diff = np.ptp(data[:, 0]), np.ptp(data[:, 1]), np.ptp(data[:, 2])
+ x_cent, y_cent, z_cent = x_min+x_diff/2., y_min+y_diff/2., z_min+z_diff/2.
+ # Create regular grid
+ xs = np.arange(x_cent-x_diff, x_cent+x_diff, a)
+ ys = np.arange(y_cent-y_diff, y_cent+y_diff, a)
+ o, p = np.meshgrid(xs, ys)
+ # Create empty magnitude:
+ magnitude = np.zeros((3, len(zs), len(ys), len(xs)))
+
+ # WITH MASKING OF THE CENTER (SYMMETRIC):
+ iz_x = np.concatenate([np.linspace(-4.95, -4.95, 50),
+ np.linspace(-4.95, 0, 50),
+ np.linspace(0, 4.95, 50),
+ np.linspace(4.95, 4.95, 50),
+ np.linspace(-4.95, 0, 50),
+ np.linspace(0, 4.95, 50), ])
+ iz_y = np.concatenate([np.linspace(-2.88, 2.88, 50),
+ np.linspace(2.88, 5.7, 50),
+ np.linspace(5.7, 2.88, 50),
+ np.linspace(2.88, -2.88, 50),
+ np.linspace(-2.88, -5.7, 50),
+ np.linspace(-5.7, -2.88, 50), ])
+ for i, z in tqdm(enumerate(zs), total=len(zs)):
+ subdata = data[data[:, 2] == z, :]
+ for j in range(3): # For all 3 components!
+ gridded_subdata = griddata(np.concatenate([subdata[:, 0], iz_x]),
+ np.concatenate([subdata[:, 1], iz_y]),
+ np.concatenate([subdata[:, 3 + j], np.zeros(len(iz_x))]),
+ o, p)
+ magnitude[j, i, :, :] = gridded_subdata.filled(fill_value=0)
+
+# # WITH MASKING OF THE CENTER (ASYMMETRIC):
+# iz_x = np.concatenate([np.linspace(-5.88, -5.88, 50),
+# np.linspace(-5.88, 0, 50),
+# np.linspace(0, 5.88, 50),
+# np.linspace(5.88, 5.88, 50),
+# np.linspace(5.88, 0, 50),
+# np.linspace(0, -5.88, 50),])
+# iz_y = np.concatenate([np.linspace(-2.10, 4.50, 50),
+# np.linspace(4.50, 7.90, 50),
+# np.linspace(7.90, 4.50, 50),
+# np.linspace(4.50, -2.10, 50),
+# np.linspace(-2.10, -5.50, 50),
+# np.linspace(-5.50, -2.10, 50), ])
+# for i, z in tqdm(enumerate(zs), total=len(zs)):
+# subdata = data[data[:, 2] == z, :]
+# for j in range(3): # For all 3 components!
+# gridded_subdata = griddata(np.concatenate([subdata[:, 0], iz_x]),
+# np.concatenate([subdata[:, 1], iz_y]),
+# np.concatenate([subdata[:, 3 + j], np.zeros(len(iz_x))]),
+# o, p)
+# magnitude[j, i, :, :] = gridded_subdata.filled(fill_value=0)
+
+# # WITHOUT MASKING OF THE CENTER:
+# for i, z in tqdm(enumerate(zs), total=len(zs)):
+# subdata = data[data[:, 2] == z, :]
+# for j in range(3): # For all 3 components!
+# gridded_subdata = griddata(subdata[:, 0], subdata[:, 1], subdata[:, 3 + j], o, p)
+# magnitude[j, i, :, :] = gridded_subdata.filled(fill_value=0)
+
+ # Creating MagData object:
+ mag_data = MagData(0.2*10, magnitude)
+ mag_data.save_to_netcdf4(PATH+'.nc')
+else:
+ mag_data = MagData.load_from_netcdf4(PATH+'.nc')
+mag_data.quiver_plot()
+###################################################################################################
+# Turn magnetization around by 90° around x-axis:
+#magnitude_new = np.zeros((3, mag_data.dim[1], mag_data.dim[0], mag_data.dim[2]))
+#for i in range(mag_data.dim[2]):
+# x_rot = np.rot90(mag_data.magnitude[0, ..., i]).copy()
+# y_rot = np.rot90(mag_data.magnitude[1, ..., i]).copy()
+# z_rot = np.rot90(mag_data.magnitude[2, ..., i]).copy()
+# magnitude_new[0, ..., i] = x_rot
+# magnitude_new[1, ..., i] = z_rot
+# magnitude_new[2, ..., i] = -y_rot
+#mag_data.magnitude = magnitude_new
+#mag_data.save_to_netcdf4(PATH+'_lying_down.nc')
+
+
+dim = mag_data.dim
+dim_uv = (500, 200)
+angles = [0, 20, 40, 60]
+
+mag_data_xy = mag_data.copy()
+mag_data_xy.magnitude[2] = 0
+
+mag_data_z = mag_data.copy()
+mag_data_z.magnitude[0] = 0
+mag_data_z.magnitude[1] = 0
+
+# Iterate over all angles:
+for angle in angles:
+ angle_rad = np.pi/2 + angle*np.pi/180
+ projector = XTiltProjector(dim, angle_rad, dim_uv)
+ mag_proj = projector(mag_data_z)
+ phase_map = PhaseMapperRDFC(Kernel(mag_data.a, projector.dim_uv))(mag_proj)
+ phase_map.display_combined('Phase Map Nanowire Tip', gain=gain,
+ interpolation='bilinear')
+ plt.savefig(PATH+'_nanowire_z_xtilt_{}.png'.format(angle), dpi=500)
+ mag_proj.scale_down(2)
+ axis = mag_proj.quiver_plot()
+ plt.savefig(PATH+'_nanowire_z_mag_xtilt_{}.png'.format(angle), dpi=500)
+ axis = mag_proj.quiver_plot(log=True)
+ plt.savefig(PATH+'_nanowire_z_mag_log_xtilt_{}.png'.format(angle), dpi=500)
+ # Close plots:
+ plt.close('all')
+
+
+#mag_data.scale_down(2)
+#mag_data.quiver_plot3d()
+#
+#mag_data_xy.scale_down(2)
+#mag_data_xy.quiver_plot3d()
+#
+#mag_data_z.scale_down(2)
+#mag_data_z.quiver_plot3d()
diff --git a/scripts/collaborations/vtk_conversion_nanowire_full_90deg.py b/scripts/collaborations/vtk_conversion_nanowire_full_90deg.py
index d1ab275aad9df9c51f59f5ae1747c702565fe4f3..2f01c98d3f9c5b4e7567521b426f9089be14c7e6 100644
--- a/scripts/collaborations/vtk_conversion_nanowire_full_90deg.py
+++ b/scripts/collaborations/vtk_conversion_nanowire_full_90deg.py
@@ -128,13 +128,13 @@ dim_uv = (600, 150)
angles = [0, 10, 20, 30, 40, 50, 60]
# Turn magnetization around by 90° around x-axis:
magnitude_new = np.zeros((3, mag_data.dim[1], mag_data.dim[0], mag_data.dim[2]))
-for i in range(mag_data.magnitude.shape[2]):
+for i in range(mag_data.dim[2]):
x_rot = np.rot90(mag_data.magnitude[0, ..., i]).copy()
y_rot = np.rot90(mag_data.magnitude[1, ..., i]).copy()
z_rot = np.rot90(mag_data.magnitude[2, ..., i]).copy()
magnitude_new[0, ..., i] = x_rot
- magnitude_new[1, ..., i] = -z_rot
- magnitude_new[2, ..., i] = y_rot
+ magnitude_new[1, ..., i] = z_rot
+ magnitude_new[2, ..., i] = -y_rot
mag_data.magnitude = magnitude_new
mag_data.save_to_netcdf4(PATH+'_lying_down.nc')
# Iterate over all angles:
diff --git a/scripts/reconstruction/reconstruct_random_pixels.py b/scripts/reconstruction/reconstruct_random_pixels.py
index 3aec76e979c3da6f6463c1e7f9fd628b040948bc..7a1214865f312191b70a1ce503a9f9f464640ecf 100644
--- a/scripts/reconstruction/reconstruct_random_pixels.py
+++ b/scripts/reconstruction/reconstruct_random_pixels.py
@@ -13,6 +13,8 @@ import pyramid
import pyramid.magcreator as mc
from pyramid.magdata import MagData
from pyramid.phasemapper import pm
+from pyramid.dataset import DataSet
+from pyramid.projector import SimpleProjector
import pyramid.reconstruction as rc
import logging
@@ -45,7 +47,10 @@ phase_map = pm(mag_data)
phase_map.display_combined('Generated Distribution', gain=10)
# Reconstruct the magnetic distribution:
-mag_data_rec = rc.optimize_simple_leastsq(phase_map, mag_data.get_mask(), b_0, lam=1E-4, order=1)
+
+data = DataSet(a, dim, b_0, mag_data.get_mask())
+data.append(phase_map, SimpleProjector(dim))
+mag_data_rec, cost = rc.optimize_linear(data, max_iter=100)
# Display the reconstructed phase map and holography image:
phase_map_rec = pm(mag_data_rec)
diff --git a/scripts/reconstruction/reconstruction_linear_test.py b/scripts/reconstruction/reconstruction_linear_test.py
index 7da8b7aba61f3b45e4ff1b69935e026077aec624..264625a7b62746878fa8652723b4da2e6d416f50 100644
--- a/scripts/reconstruction/reconstruction_linear_test.py
+++ b/scripts/reconstruction/reconstruction_linear_test.py
@@ -36,7 +36,6 @@ dim = (32, 32, 32)
dim_uv = dim[1:3]
count = 16
lam = 1E-4
-use_fftw = True
center = (dim[0]/2-0.5, int(dim[1]/2)-0.5, int(dim[2]/2)-0.5)
radius_core = dim[1]/8
radius_shell = dim[1]/4
@@ -52,10 +51,6 @@ shape = mc.Shapes.disc(dim, center, radius_shell, height)
magnitude = mc.create_mag_dist_vortex(shape)
mag_data = MagData(a, magnitude)
-mag_data.quiver_plot('z-projection', proj_axis='z')
-mag_data.quiver_plot('y-projection', proj_axis='y')
-mag_data.quiver_plot3d('Original distribution')
-
tilts_full = np.linspace(-pi/2, pi/2, num=count/2, endpoint=False)
tilts_miss = np.linspace(-pi/3, pi/3, num=count/2, endpoint=False)
@@ -77,7 +72,7 @@ noise = 0
print('--Setting up data collection')
mask = mag_data.get_mask()
-data = DataSet(a, dim, b_0, mask, use_fftw=use_fftw)
+data = DataSet(a, dim, b_0, mask)
data.projectors = projectors
data.phase_maps = data.create_phase_maps(mag_data)
@@ -97,26 +92,43 @@ print 'Cost:', cost.chisq
###################################################################################################
print('--Plot stuff')
-mag_data_opt.quiver_plot3d('Reconstructed distribution')
+limit = 1.2
+mag_data.quiver_plot3d('Original distribution', limit=limit)
+mag_data_opt.quiver_plot3d('Reconstructed distribution', limit=limit)
(mag_data_opt - mag_data).quiver_plot3d('Difference')
phase_maps_opt = data.create_phase_maps(mag_data_opt)
-# TODO: iterations in jutil is one to many!
from pyramid.diagnostics import Diagnostics
+from matplotlib.patches import Rectangle
-diag = Diagnostics(mag_data_opt.mag_vec, cost)
-
-print 'std:', diag.std
-gain_maps = diag.get_gain_maps()
-gain_maps[count//2].display_phase()
-diag.get_avg_kernel().quiver_plot3d()
-measure_contribution = diag.measure_contribution
-diag.set_position(cost.data_set.mask.sum()//2)
+diag = Diagnostics(mag_data_opt.mag_vec, cost, max_iter=2000)
+print 'position:', diag.pos
+print 'std:', diag.std
+gain_maps = diag.get_gain_row_maps()
+axis, cbar = gain_maps[count//2].display_phase()
+axis.add_patch(Rectangle((diag.pos[3], diag.pos[2]), 1, 1, linewidth=2, color='g', fill=False))
+cbar.set_label(u'magnetization/phase [1/rad]', fontsize=15)
+diag.get_avg_kern_row().quiver_plot3d()
+mcon = diag.measure_contribution
+print 'measurement contr. (min - max): {:.2f} - {:.2f}'.format(mcon.min(), mcon.max())
+px_avrg, fwhm, (x, y, z) = diag.calculate_averaging()
+print 'px_avrg:', px_avrg
+print 'fwhm:', fwhm
+
+diag.pos = (1, dim[0]//2, dim[1]//2, dim[2]//2)
+
+print 'position:', diag.pos
print 'std:', diag.std
-gain_maps = diag.get_gain_maps()
-gain_maps[count//2].display_phase()
-diag.get_avg_kernel().quiver_plot3d()
-measure_contribution = diag.measure_contribution
+gain_maps = diag.get_gain_row_maps()
+axis, cbar = gain_maps[count//2].display_phase()
+axis.add_patch(Rectangle((diag.pos[3], diag.pos[2]), 1, 1, linewidth=2, color='g', fill=False))
+cbar.set_label(u'magnetization/phase [1/rad]', fontsize=15)
+diag.get_avg_kern_row().quiver_plot3d()
+mcon = diag.measure_contribution
+print 'measurement contr. (min - max): {:.2f} - {:.2f}'.format(mcon.min(), mcon.max())
+px_avrg, fwhm, (x, y, z) = diag.calculate_averaging()
+print 'px_avrg:', px_avrg
+print 'fwhm:', fwhm
diff --git a/scripts/reconstruction/reconstruction_linear_test_batch.py b/scripts/reconstruction/reconstruction_linear_test_batch.py
index 3f310a86ea34df112cda566434f7343bb5ce14d7..2de7ce43ea63a1e2451690d2c821be9b62a252c8 100644
--- a/scripts/reconstruction/reconstruction_linear_test_batch.py
+++ b/scripts/reconstruction/reconstruction_linear_test_batch.py
@@ -108,7 +108,7 @@ for i, configuration in enumerate(itertools.product(*config_list)):
# DataSet:
data = DataSet(a, dim, b_0, mask)
# Tilts:
- tilts = np.arange(-max_tilt/180.*pi, max_tilt/180.*pi, tilt_step/180.*pi)
+ tilts = np.deg2rad(np.arange(-max_tilt, max_tilt, tilt_step))
# Projectors:
projectors = []
if xy_tilts[xy_key][0]:
diff --git a/scripts/reconstruction/reconstruction_sparse_cg_singular_values.py b/scripts/reconstruction/reconstruction_sparse_cg_singular_values.py
index 0ccee36c6e92b398d1120b8bd28cece139cd13b0..df399152aadfe3cbce9a0db85e88d7b80241e711 100644
--- a/scripts/reconstruction/reconstruction_sparse_cg_singular_values.py
+++ b/scripts/reconstruction/reconstruction_sparse_cg_singular_values.py
@@ -64,7 +64,7 @@ F = ForwardModel(data)
M = np.asmatrix([F.jac_dot(None, np.eye(3*size_3d)[:, i]) for i in range(3*size_3d)]).T
# MTM = M.T * M + lam * np.asmatrix(np.eye(3*size_3d))
-U, s, V = np.linalg.svd(M) # TODO: M or MTM?
+U, s, V = np.linalg.svd(M) # TODO: M or MTM
for i in range(len(s)):
print 'Singular value:', s[i]
diff --git a/scripts/test methods/histo_norm.py b/scripts/test methods/histo_norm.py
new file mode 100644
index 0000000000000000000000000000000000000000..7a9dea36e07a9c267f60c9f380169ddaabda88e5
--- /dev/null
+++ b/scripts/test methods/histo_norm.py
@@ -0,0 +1,39 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Thu Jan 22 14:30:14 2015
+
+@author: Jan
+"""
+
+
+import matplotlib.pyplot as plt
+
+from scipy import misc
+import scipy as sp
+
+
+im = misc.lena()
+
+plt.figure()
+plt.imshow(im, cmap=plt.cm.gray)
+plt.show()
+
+
+def histeq(im, nbr_bins=256):
+
+ # get image histogram
+ imhist, bins = sp.histogram(im.flatten(), nbr_bins, normed=True)
+ cdf = imhist.cumsum() # cumulative distribution function
+ cdf = 255 * cdf / cdf[-1] # normalize
+
+ # use linear interpolation of cdf to find new pixel values
+ im2 = sp.interp(im.flatten(), bins[:-1], cdf)
+
+ return im2.reshape(im.shape), cdf
+
+
+im2, cdf = histeq(im)
+
+plt.figure()
+plt.imshow(im2, cmap=plt.cm.gray)
+plt.show()
diff --git a/tests/test_analytic.py b/tests/test_analytic.py
index 562ee620665e457d916fd6619a20d678258a0b4c..dbdd9d65eb3bdb5ee0bbeb37d0d744d9ce6da2d6 100644
--- a/tests/test_analytic.py
+++ b/tests/test_analytic.py
@@ -4,9 +4,9 @@
import os
import unittest
+
import numpy as np
from numpy import pi
-
from numpy.testing import assert_allclose
import pyramid.analytic as an
diff --git a/tests/test_compliance.py b/tests/test_compliance.py
index 6ffe674cc35119cb3ffd1b20f7ade7f4e41bd8fa..b92af3f427fab6d63ab99e01a00704c2b61e8f75 100644
--- a/tests/test_compliance.py
+++ b/tests/test_compliance.py
@@ -22,8 +22,8 @@ class TestCaseCompliance(unittest.TestCase):
for root, dirs, files in os.walk(rootdir):
for filename in files:
if ((filename.endswith('.py') or filename.endswith('.pyx'))
- and root != os.path.join(self.path, 'scripts', 'gui')):
- filepaths.append(os.path.join(root, filename))
+ and root != os.path.join(self.path, 'scripts', 'gui')):
+ filepaths.append(os.path.join(root, filename))
return filepaths
def test_pep8(self):
diff --git a/tests/test_kernel.py b/tests/test_kernel.py
index 39cb1e12d4105bcffde431aa6b90c81002b6e95e..cfc0bbd9629d618e0ec8aacb2315b5f10271f55d 100644
--- a/tests/test_kernel.py
+++ b/tests/test_kernel.py
@@ -1,2 +1,41 @@
# -*- coding: utf-8 -*-
-"""Created on Mon Jan 20 20:00:50 2014 @author: Jan"""
+"""Testcase for the magdata module."""
+
+
+import os
+import unittest
+
+import numpy as np
+from numpy.testing import assert_allclose
+
+from pyramid.kernel import Kernel
+
+
+class TestCaseKernel(unittest.TestCase):
+
+ def setUp(self):
+ self.path = os.path.join(os.path.dirname(os.path.realpath(__file__)), 'test_kernel')
+ self.kernel = Kernel(1., dim_uv=(4, 4), b_0=1., geometry='disc')
+
+ def tearDown(self):
+ self.path = None
+ self.kernel = None
+
+ def test_kernel(self):
+ np.save(os.path.join(self.path, 'ref_u'), self.kernel.u)
+ np.save(os.path.join(self.path, 'ref_v'), self.kernel.v)
+ np.save(os.path.join(self.path, 'ref_u_fft'), self.kernel.u_fft)
+ np.save(os.path.join(self.path, 'ref_v_fft'), self.kernel.v_fft)
+ ref_u = np.load(os.path.join(self.path, 'ref_u.npy'))
+ ref_v = np.load(os.path.join(self.path, 'ref_v.npy'))
+ ref_u_fft = np.load(os.path.join(self.path, 'ref_u_fft.npy'))
+ ref_v_fft = np.load(os.path.join(self.path, 'ref_v_fft.npy'))
+ assert_allclose(self.kernel.u, ref_u, err_msg='Unexpected behavior in u')
+ assert_allclose(self.kernel.v, ref_v, err_msg='Unexpected behavior in v')
+ assert_allclose(self.kernel.u_fft, ref_u_fft, err_msg='Unexpected behavior in u_fft')
+ assert_allclose(self.kernel.v_fft, ref_v_fft, err_msg='Unexpected behavior in v_fft')
+
+
+if __name__ == '__main__':
+ suite = unittest.TestLoader().loadTestsFromTestCase(TestCaseKernel)
+ unittest.TextTestRunner(verbosity=2).run(suite)
diff --git a/tests/test_kernel/ref_u.npy b/tests/test_kernel/ref_u.npy
new file mode 100644
index 0000000000000000000000000000000000000000..794faa4622119f38f46bdf2acc9fad231f56ec4e
Binary files /dev/null and b/tests/test_kernel/ref_u.npy differ
diff --git a/tests/test_kernel/ref_u_fft.npy b/tests/test_kernel/ref_u_fft.npy
new file mode 100644
index 0000000000000000000000000000000000000000..30a4c4d9f3fc0560212f7bb98b18930ecd932926
Binary files /dev/null and b/tests/test_kernel/ref_u_fft.npy differ
diff --git a/tests/test_kernel/ref_v.npy b/tests/test_kernel/ref_v.npy
new file mode 100644
index 0000000000000000000000000000000000000000..c854a27d8debf3f664bb65695088dba2ecb41630
Binary files /dev/null and b/tests/test_kernel/ref_v.npy differ
diff --git a/tests/test_kernel/ref_v_fft.npy b/tests/test_kernel/ref_v_fft.npy
new file mode 100644
index 0000000000000000000000000000000000000000..d7588d2717005ee4e12782e7b8d6e8873a0aa17b
Binary files /dev/null and b/tests/test_kernel/ref_v_fft.npy differ
diff --git a/tests/test_magcreator.py b/tests/test_magcreator.py
index 137b6e670636ee2c825c50e85155fdbafa16a7bc..4abcb950977140dd5bc0db4bb4fd68d2ebfcadc8 100644
--- a/tests/test_magcreator.py
+++ b/tests/test_magcreator.py
@@ -31,11 +31,27 @@ class TestCaseMagCreator(unittest.TestCase):
assert_allclose(test_disc_x, np.load(os.path.join(self.path, 'ref_disc_x.npy')),
err_msg='Created disc in x-direction does not match expectation!')
+ def test_shape_ellipse(self):
+ test_ellipse_z = mc.Shapes.ellipse((7, 8, 9), (3, 4, 5), (3, 5), 1, 'z')
+ test_ellipse_y = mc.Shapes.ellipse((7, 8, 9), (3, 4, 5), (3, 5), 1, 'y')
+ test_ellipse_x = mc.Shapes.ellipse((7, 8, 9), (3, 4, 5), (3, 5), 1, 'x')
+ assert_allclose(test_ellipse_z, np.load(os.path.join(self.path, 'ref_ellipse_z.npy')),
+ err_msg='Created ellipse does not match expectation!')
+ assert_allclose(test_ellipse_y, np.load(os.path.join(self.path, 'ref_ellipse_y.npy')),
+ err_msg='Created ellipse does not match expectation!')
+ assert_allclose(test_ellipse_x, np.load(os.path.join(self.path, 'ref_ellipse_x.npy')),
+ err_msg='Created ellipse does not match expectation!')
+
def test_shape_sphere(self):
test_sphere = mc.Shapes.sphere((5, 6, 7), (2, 3, 4), 2)
assert_allclose(test_sphere, np.load(os.path.join(self.path, 'ref_sphere.npy')),
err_msg='Created sphere does not match expectation!')
+ def test_shape_ellipsoid(self):
+ test_ellipsoid = mc.Shapes.ellipsoid((7, 8, 9), (3, 4, 4), (3, 5, 7))
+ assert_allclose(test_ellipsoid, np.load(os.path.join(self.path, 'ref_ellipsoid.npy')),
+ err_msg='Created ellipsoid does not match expectation!')
+
def test_shape_filament(self):
test_filament_z = mc.Shapes.filament((5, 6, 7), (2, 3), 'z')
test_filament_y = mc.Shapes.filament((5, 6, 7), (2, 3), 'y')
diff --git a/tests/test_magcreator/ref_ellipse_x.npy b/tests/test_magcreator/ref_ellipse_x.npy
new file mode 100644
index 0000000000000000000000000000000000000000..382ada72ce030957e85c6131100432acada9278b
Binary files /dev/null and b/tests/test_magcreator/ref_ellipse_x.npy differ
diff --git a/tests/test_magcreator/ref_ellipse_y.npy b/tests/test_magcreator/ref_ellipse_y.npy
new file mode 100644
index 0000000000000000000000000000000000000000..c83d0cb376130c7938da3536ba087e3289877170
Binary files /dev/null and b/tests/test_magcreator/ref_ellipse_y.npy differ
diff --git a/tests/test_magcreator/ref_ellipse_z.npy b/tests/test_magcreator/ref_ellipse_z.npy
new file mode 100644
index 0000000000000000000000000000000000000000..75439a6e20aa9a53158aec3a1b5c2602ffc979ab
Binary files /dev/null and b/tests/test_magcreator/ref_ellipse_z.npy differ
diff --git a/tests/test_magcreator/ref_ellipsoid.npy b/tests/test_magcreator/ref_ellipsoid.npy
new file mode 100644
index 0000000000000000000000000000000000000000..8abca7c48658e9f9df018d7869ba2dc4461b1d9a
Binary files /dev/null and b/tests/test_magcreator/ref_ellipsoid.npy differ
diff --git a/tests/test_magdata.py b/tests/test_magdata.py
index 14d96606f5259af2d4f9909bb1844eecc5f920b2..a178de31891a93609248c941caa5e48032d89499 100644
--- a/tests/test_magdata.py
+++ b/tests/test_magdata.py
@@ -4,7 +4,9 @@
import os
import unittest
+
import numpy as np
+from numpy.testing import assert_allclose
from pyramid.magdata import MagData
@@ -13,77 +15,80 @@ class TestCaseMagData(unittest.TestCase):
def setUp(self):
self.path = os.path.join(os.path.dirname(os.path.realpath(__file__)), 'test_magdata')
- magnitude = (np.zeros((4, 4, 4)), np.zeros((4, 4, 4)), np.zeros((4, 4, 4)))
- magnitude[0][1:-1, 1:-1, 1:-1] = 1
- magnitude[1][1:-1, 1:-1, 1:-1] = 1
- magnitude[2][1:-1, 1:-1, 1:-1] = 1
+ magnitude = np.zeros((3, 4, 4, 4))
+ magnitude[:, 1:-1, 1:-1, 1:-1] = 1
self.mag_data = MagData(10.0, magnitude)
def tearDown(self):
self.path = None
self.mag_data = None
- def test_add_magnitude(self):
- reference = (np.ones((4, 4, 4)), np.ones((4, 4, 4)), np.ones((4, 4, 4)))
- self.mag_data.add_magnitude(reference)
- reference[0][1:-1, 1:-1, 1:-1] = 2
- reference[1][1:-1, 1:-1, 1:-1] = 2
- reference[2][1:-1, 1:-1, 1:-1] = 2
- np.testing.assert_equal(self.mag_data.magnitude, reference,
- 'Unexpected behavior in add_magnitude()!')
+ def test_copy(self):
+ mag_data = self.mag_data
+ mag_data_copy = self.mag_data.copy()
+ assert mag_data == self.mag_data, 'Unexpected behaviour in copy()!'
+ assert mag_data_copy != self.mag_data, 'Unexpected behaviour in copy()!'
+
+ def test_scale_down(self):
+ self.mag_data.scale_down()
+ reference = 1/8. * np.ones((3, 2, 2, 2))
+ assert_allclose(self.mag_data.magnitude, reference,
+ err_msg='Unexpected behavior in scale_down()!')
+ assert_allclose(self.mag_data.a, 20,
+ err_msg='Unexpected behavior in scale_down()!')
+
+ def test_scale_up(self):
+ self.mag_data.scale_up()
+ reference = np.zeros((3, 8, 8, 8))
+ reference[:, 2:6, 2:6, 2:6] = 1
+ assert_allclose(self.mag_data.magnitude, reference,
+ err_msg='Unexpected behavior in scale_down()!')
+ assert_allclose(self.mag_data.a, 5,
+ err_msg='Unexpected behavior in scale_down()!')
+
+ def test_pad(self):
+ reference = self.mag_data.magnitude.copy()
+ self.mag_data.pad(1, 1, 1)
+ reference = np.pad(reference, ((0, 0), (1, 1), (1, 1), (1, 1)), mode='constant')
+ assert_allclose(self.mag_data.magnitude, reference,
+ err_msg='Unexpected behavior in scale_down()!')
def test_get_mask(self):
mask = self.mag_data.get_mask()
reference = np.zeros((4, 4, 4))
reference[1:-1, 1:-1, 1:-1] = True
- np.testing.assert_equal(mask, reference, 'Unexpected behavior in get_mask()!')
+ assert_allclose(mask, reference,
+ err_msg='Unexpected behavior in get_mask()!')
def test_get_vector(self):
mask = self.mag_data.get_mask()
vector = self.mag_data.get_vector(mask)
- reference = np.ones(np.count_nonzero(self.mag_data.magnitude[0])*3)
- np.testing.assert_equal(vector, reference, 'Unexpected behavior in get_mask()!')
+ reference = np.ones(np.sum(mask)*3)
+ assert_allclose(vector, reference,
+ err_msg='Unexpected behavior in get_vector()!')
def test_set_vector(self):
mask = self.mag_data.get_mask()
- vector = 2 * np.ones(np.count_nonzero(self.mag_data.magnitude[0])*3)
- self.mag_data.set_vector(mask, vector)
- reference = (np.zeros((4, 4, 4)), np.zeros((4, 4, 4)), np.zeros((4, 4, 4)))
- reference[0][1:-1, 1:-1, 1:-1] = 2
- reference[1][1:-1, 1:-1, 1:-1] = 2
- reference[2][1:-1, 1:-1, 1:-1] = 2
- np.testing.assert_equal(self.mag_data.magnitude, reference,
- 'Unexpected behavior in set_mask()!')
-
- def test_scale_down(self):
- self.mag_data.scale_down()
- reference = (1/8. * np.ones((2, 2, 2)),
- 1/8. * np.ones((2, 2, 2)),
- 1/8. * np.ones((2, 2, 2)))
- np.testing.assert_equal(self.mag_data.magnitude, reference,
- 'Unexpected behavior in scale_down()!')
- np.testing.assert_equal(self.mag_data.res, 20, 'Unexpected behavior in scale_down()!')
+ vector = 2 * np.ones(np.sum(mask)*3)
+ self.mag_data.set_vector(vector, mask)
+ reference = np.zeros((3, 4, 4, 4))
+ reference[:, 1:-1, 1:-1, 1:-1] = 2
+ assert_allclose(self.mag_data.magnitude, reference,
+ err_msg='Unexpected behavior in set_vector()!')
def test_load_from_llg(self):
- self.mag_data = MagData.load_from_llg(os.path.join(self.path, 'ref_mag_data.txt'))
- reference = (np.zeros((4, 4, 4)), np.zeros((4, 4, 4)), np.zeros((4, 4, 4)))
- reference[0][1:-1, 1:-1, 1:-1] = 1
- reference[1][1:-1, 1:-1, 1:-1] = 1
- reference[2][1:-1, 1:-1, 1:-1] = 1
- np.testing.assert_equal(self.mag_data.magnitude, reference,
- 'Unexpected behavior in load_from_llg()!')
- np.testing.assert_equal(self.mag_data.res, 10, 'Unexpected behavior in load_from_llg()!')
+ mag_data = MagData.load_from_llg(os.path.join(self.path, 'ref_mag_data.txt'))
+ assert_allclose(mag_data.magnitude, self.mag_data.magnitude,
+ err_msg='Unexpected behavior in load_from_llg()!')
+ assert_allclose(mag_data.a, self.mag_data.a,
+ err_msg='Unexpected behavior in load_from_llg()!')
def test_load_from_netcdf4(self):
- self.mag_data = MagData.load_from_netcdf4(os.path.join(self.path, 'ref_mag_data.nc'))
- reference = (np.zeros((4, 4, 4)), np.zeros((4, 4, 4)), np.zeros((4, 4, 4)))
- reference[0][1:-1, 1:-1, 1:-1] = 1
- reference[1][1:-1, 1:-1, 1:-1] = 1
- reference[2][1:-1, 1:-1, 1:-1] = 1
- np.testing.assert_equal(self.mag_data.magnitude, reference,
- 'Unexpected behavior in load_from_netcdf4()!')
- np.testing.assert_equal(self.mag_data.res, 10,
- 'Unexpected behavior in load_from_netcdf4()!')
+ mag_data = MagData.load_from_netcdf4(os.path.join(self.path, 'ref_mag_data.nc'))
+ assert_allclose(mag_data.magnitude, self.mag_data.magnitude,
+ err_msg='Unexpected behavior in load_from_netcdf4()!')
+ assert_allclose(mag_data.a, self.mag_data.a,
+ err_msg='Unexpected behavior in load_from_netcdf4()!')
if __name__ == '__main__':
diff --git a/tests/test_magdata/ref_mag_data.nc b/tests/test_magdata/ref_mag_data.nc
index 2b7f60c61e7e652cd00f3f035cf2ac37bc5adec2..200f66fa5f6c43f56e3746ccbc70013b6f206164 100644
Binary files a/tests/test_magdata/ref_mag_data.nc and b/tests/test_magdata/ref_mag_data.nc differ
diff --git a/tests/test_projector.py b/tests/test_projector.py
index 6da902fdfb4774031fdc3c1f84039f1dec9d6e39..b328c024452ea24ce745dd8d39c73a39fa04793f 100644
--- a/tests/test_projector.py
+++ b/tests/test_projector.py
@@ -4,9 +4,11 @@
import os
import unittest
-import numpy as np
-import pyramid.projector as pj
+from numpy import pi
+from numpy.testing import assert_allclose
+
+from pyramid.projector import XTiltProjector, YTiltProjector, SimpleProjector
from pyramid.magdata import MagData
@@ -14,31 +16,53 @@ class TestCaseProjector(unittest.TestCase):
def setUp(self):
self.path = os.path.join(os.path.dirname(os.path.realpath(__file__)), 'test_projector')
- self.mag_data = MagData.load_from_llg(os.path.join(self.path, 'ref_mag_data.txt'))
+ self.mag_data = MagData.load_from_netcdf4(os.path.join(self.path, 'ref_mag_data.nc'))
def tearDown(self):
self.path = None
self.mag_data = None
- def test_simple_axis_projection(self):
- z_proj_ref = (np.loadtxt(os.path.join(self.path, 'ref_proj_z.txt')))
- y_proj_ref = (np.loadtxt(os.path.join(self.path, 'ref_proj_y.txt')))
- x_proj_ref = (np.loadtxt(os.path.join(self.path, 'ref_proj_x.txt')))
- z_proj = pj.simple_axis_projection(self.mag_data, 'z')
- y_proj = pj.simple_axis_projection(self.mag_data, 'y')
- x_proj = pj.simple_axis_projection(self.mag_data, 'x')
- np.testing.assert_equal(z_proj[0], z_proj_ref, '1')
- np.testing.assert_equal(z_proj[1], z_proj_ref, '2')
- np.testing.assert_equal(z_proj[2], z_proj_ref, '3')
- np.testing.assert_equal(y_proj[0], y_proj_ref, '4')
- np.testing.assert_equal(y_proj[1], y_proj_ref, '5')
- np.testing.assert_equal(y_proj[2], y_proj_ref, '6')
- np.testing.assert_equal(x_proj[0], x_proj_ref, '7')
- np.testing.assert_equal(x_proj[1], x_proj_ref, '8')
- np.testing.assert_equal(x_proj[2], x_proj_ref, '9')
-
- def test_single_axis_projection(self):
- raise AssertionError
+ def test_simple_projector(self):
+ mag_proj_z = SimpleProjector(self.mag_data.dim, axis='z')(self.mag_data)
+ mag_proj_y = SimpleProjector(self.mag_data.dim, axis='y')(self.mag_data)
+ mag_proj_x = SimpleProjector(self.mag_data.dim, axis='x')(self.mag_data)
+ mag_proj_z_ref = MagData.load_from_netcdf4(os.path.join(self.path, 'ref_mag_proj_z.nc'))
+ mag_proj_y_ref = MagData.load_from_netcdf4(os.path.join(self.path, 'ref_mag_proj_y.nc'))
+ mag_proj_x_ref = MagData.load_from_netcdf4(os.path.join(self.path, 'ref_mag_proj_x.nc'))
+ assert_allclose(mag_proj_z.magnitude, mag_proj_z_ref.magnitude,
+ err_msg='Unexpected behaviour in SimpleProjector (z-axis)')
+ assert_allclose(mag_proj_y.magnitude, mag_proj_y_ref.magnitude,
+ err_msg='Unexpected behaviour in SimpleProjector (y-axis)')
+ assert_allclose(mag_proj_x.magnitude, mag_proj_x_ref.magnitude,
+ err_msg='Unexpected behaviour in SimpleProjector (x-axis)')
+
+ def test_x_tilt_projector(self):
+ mag_proj_00 = XTiltProjector(self.mag_data.dim, tilt=0)(self.mag_data)
+ mag_proj_45 = XTiltProjector(self.mag_data.dim, tilt=pi/4)(self.mag_data)
+ mag_proj_90 = XTiltProjector(self.mag_data.dim, tilt=pi/2)(self.mag_data)
+ mag_proj_00_ref = MagData.load_from_netcdf4(os.path.join(self.path, 'ref_mag_proj_x00.nc'))
+ mag_proj_45_ref = MagData.load_from_netcdf4(os.path.join(self.path, 'ref_mag_proj_x45.nc'))
+ mag_proj_90_ref = MagData.load_from_netcdf4(os.path.join(self.path, 'ref_mag_proj_x90.nc'))
+ assert_allclose(mag_proj_00.magnitude, mag_proj_00_ref.magnitude,
+ err_msg='Unexpected behaviour in XTiltProjector (0°)')
+ assert_allclose(mag_proj_45.magnitude, mag_proj_45_ref.magnitude,
+ err_msg='Unexpected behaviour in XTiltProjector (45°)')
+ assert_allclose(mag_proj_90.magnitude, mag_proj_90_ref.magnitude,
+ err_msg='Unexpected behaviour in XTiltProjector (90°)')
+
+ def test_y_tilt_projector(self):
+ mag_proj_00 = YTiltProjector(self.mag_data.dim, tilt=0)(self.mag_data)
+ mag_proj_45 = YTiltProjector(self.mag_data.dim, tilt=pi/4)(self.mag_data)
+ mag_proj_90 = YTiltProjector(self.mag_data.dim, tilt=pi/2)(self.mag_data)
+ mag_proj_00_ref = MagData.load_from_netcdf4(os.path.join(self.path, 'ref_mag_proj_y00.nc'))
+ mag_proj_45_ref = MagData.load_from_netcdf4(os.path.join(self.path, 'ref_mag_proj_y45.nc'))
+ mag_proj_90_ref = MagData.load_from_netcdf4(os.path.join(self.path, 'ref_mag_proj_y90.nc'))
+ assert_allclose(mag_proj_00.magnitude, mag_proj_00_ref.magnitude,
+ err_msg='Unexpected behaviour in YTiltProjector (0°)')
+ assert_allclose(mag_proj_45.magnitude, mag_proj_45_ref.magnitude,
+ err_msg='Unexpected behaviour in YTiltProjector (45°)')
+ assert_allclose(mag_proj_90.magnitude, mag_proj_90_ref.magnitude,
+ err_msg='Unexpected behaviour in YTiltProjector (90°)')
if __name__ == '__main__':
diff --git a/tests/test_projector/ref_mag_data.nc b/tests/test_projector/ref_mag_data.nc
new file mode 100644
index 0000000000000000000000000000000000000000..58632b90d4720bb9c0afd65bec6fc46bc387ff1d
Binary files /dev/null and b/tests/test_projector/ref_mag_data.nc differ
diff --git a/tests/test_projector/ref_mag_data.txt b/tests/test_projector/ref_mag_data.txt
deleted file mode 100644
index 32d452bdcb2800bfdd5e6d46a6261f238d24462a..0000000000000000000000000000000000000000
--- a/tests/test_projector/ref_mag_data.txt
+++ /dev/null
@@ -1,122 +0,0 @@
-LLGFileCreator: ref_mag_data
- 6 5 4
-5.000000e-07 5.000000e-07 5.000000e-07 0.000000e+00 0.000000e+00 0.000000e+00
-1.500000e-06 5.000000e-07 5.000000e-07 0.000000e+00 0.000000e+00 0.000000e+00
-2.500000e-06 5.000000e-07 5.000000e-07 0.000000e+00 0.000000e+00 0.000000e+00
-3.500000e-06 5.000000e-07 5.000000e-07 0.000000e+00 0.000000e+00 0.000000e+00
-4.500000e-06 5.000000e-07 5.000000e-07 0.000000e+00 0.000000e+00 0.000000e+00
-5.500000e-06 5.000000e-07 5.000000e-07 0.000000e+00 0.000000e+00 0.000000e+00
-5.000000e-07 1.500000e-06 5.000000e-07 0.000000e+00 0.000000e+00 0.000000e+00
-1.500000e-06 1.500000e-06 5.000000e-07 0.000000e+00 0.000000e+00 0.000000e+00
-2.500000e-06 1.500000e-06 5.000000e-07 0.000000e+00 0.000000e+00 0.000000e+00
-3.500000e-06 1.500000e-06 5.000000e-07 0.000000e+00 0.000000e+00 0.000000e+00
-4.500000e-06 1.500000e-06 5.000000e-07 0.000000e+00 0.000000e+00 0.000000e+00
-5.500000e-06 1.500000e-06 5.000000e-07 0.000000e+00 0.000000e+00 0.000000e+00
-5.000000e-07 2.500000e-06 5.000000e-07 0.000000e+00 0.000000e+00 0.000000e+00
-1.500000e-06 2.500000e-06 5.000000e-07 0.000000e+00 0.000000e+00 0.000000e+00
-2.500000e-06 2.500000e-06 5.000000e-07 0.000000e+00 0.000000e+00 0.000000e+00
-3.500000e-06 2.500000e-06 5.000000e-07 0.000000e+00 0.000000e+00 0.000000e+00
-4.500000e-06 2.500000e-06 5.000000e-07 0.000000e+00 0.000000e+00 0.000000e+00
-5.500000e-06 2.500000e-06 5.000000e-07 0.000000e+00 0.000000e+00 0.000000e+00
-5.000000e-07 3.500000e-06 5.000000e-07 0.000000e+00 0.000000e+00 0.000000e+00
-1.500000e-06 3.500000e-06 5.000000e-07 0.000000e+00 0.000000e+00 0.000000e+00
-2.500000e-06 3.500000e-06 5.000000e-07 0.000000e+00 0.000000e+00 0.000000e+00
-3.500000e-06 3.500000e-06 5.000000e-07 0.000000e+00 0.000000e+00 0.000000e+00
-4.500000e-06 3.500000e-06 5.000000e-07 0.000000e+00 0.000000e+00 0.000000e+00
-5.500000e-06 3.500000e-06 5.000000e-07 0.000000e+00 0.000000e+00 0.000000e+00
-5.000000e-07 4.500000e-06 5.000000e-07 0.000000e+00 0.000000e+00 0.000000e+00
-1.500000e-06 4.500000e-06 5.000000e-07 0.000000e+00 0.000000e+00 0.000000e+00
-2.500000e-06 4.500000e-06 5.000000e-07 0.000000e+00 0.000000e+00 0.000000e+00
-3.500000e-06 4.500000e-06 5.000000e-07 0.000000e+00 0.000000e+00 0.000000e+00
-4.500000e-06 4.500000e-06 5.000000e-07 0.000000e+00 0.000000e+00 0.000000e+00
-5.500000e-06 4.500000e-06 5.000000e-07 0.000000e+00 0.000000e+00 0.000000e+00
-5.000000e-07 5.000000e-07 1.500000e-06 0.000000e+00 0.000000e+00 0.000000e+00
-1.500000e-06 5.000000e-07 1.500000e-06 0.000000e+00 0.000000e+00 0.000000e+00
-2.500000e-06 5.000000e-07 1.500000e-06 0.000000e+00 0.000000e+00 0.000000e+00
-3.500000e-06 5.000000e-07 1.500000e-06 0.000000e+00 0.000000e+00 0.000000e+00
-4.500000e-06 5.000000e-07 1.500000e-06 0.000000e+00 0.000000e+00 0.000000e+00
-5.500000e-06 5.000000e-07 1.500000e-06 0.000000e+00 0.000000e+00 0.000000e+00
-5.000000e-07 1.500000e-06 1.500000e-06 0.000000e+00 0.000000e+00 0.000000e+00
-1.500000e-06 1.500000e-06 1.500000e-06 1.000000e+00 1.000000e+00 1.000000e+00
-2.500000e-06 1.500000e-06 1.500000e-06 1.000000e+00 1.000000e+00 1.000000e+00
-3.500000e-06 1.500000e-06 1.500000e-06 1.000000e+00 1.000000e+00 1.000000e+00
-4.500000e-06 1.500000e-06 1.500000e-06 1.000000e+00 1.000000e+00 1.000000e+00
-5.500000e-06 1.500000e-06 1.500000e-06 0.000000e+00 0.000000e+00 0.000000e+00
-5.000000e-07 2.500000e-06 1.500000e-06 0.000000e+00 0.000000e+00 0.000000e+00
-1.500000e-06 2.500000e-06 1.500000e-06 1.000000e+00 1.000000e+00 1.000000e+00
-2.500000e-06 2.500000e-06 1.500000e-06 1.000000e+00 1.000000e+00 1.000000e+00
-3.500000e-06 2.500000e-06 1.500000e-06 1.000000e+00 1.000000e+00 1.000000e+00
-4.500000e-06 2.500000e-06 1.500000e-06 1.000000e+00 1.000000e+00 1.000000e+00
-5.500000e-06 2.500000e-06 1.500000e-06 0.000000e+00 0.000000e+00 0.000000e+00
-5.000000e-07 3.500000e-06 1.500000e-06 0.000000e+00 0.000000e+00 0.000000e+00
-1.500000e-06 3.500000e-06 1.500000e-06 1.000000e+00 1.000000e+00 1.000000e+00
-2.500000e-06 3.500000e-06 1.500000e-06 1.000000e+00 1.000000e+00 1.000000e+00
-3.500000e-06 3.500000e-06 1.500000e-06 1.000000e+00 1.000000e+00 1.000000e+00
-4.500000e-06 3.500000e-06 1.500000e-06 1.000000e+00 1.000000e+00 1.000000e+00
-5.500000e-06 3.500000e-06 1.500000e-06 0.000000e+00 0.000000e+00 0.000000e+00
-5.000000e-07 4.500000e-06 1.500000e-06 0.000000e+00 0.000000e+00 0.000000e+00
-1.500000e-06 4.500000e-06 1.500000e-06 0.000000e+00 0.000000e+00 0.000000e+00
-2.500000e-06 4.500000e-06 1.500000e-06 0.000000e+00 0.000000e+00 0.000000e+00
-3.500000e-06 4.500000e-06 1.500000e-06 0.000000e+00 0.000000e+00 0.000000e+00
-4.500000e-06 4.500000e-06 1.500000e-06 0.000000e+00 0.000000e+00 0.000000e+00
-5.500000e-06 4.500000e-06 1.500000e-06 0.000000e+00 0.000000e+00 0.000000e+00
-5.000000e-07 5.000000e-07 2.500000e-06 0.000000e+00 0.000000e+00 0.000000e+00
-1.500000e-06 5.000000e-07 2.500000e-06 0.000000e+00 0.000000e+00 0.000000e+00
-2.500000e-06 5.000000e-07 2.500000e-06 0.000000e+00 0.000000e+00 0.000000e+00
-3.500000e-06 5.000000e-07 2.500000e-06 0.000000e+00 0.000000e+00 0.000000e+00
-4.500000e-06 5.000000e-07 2.500000e-06 0.000000e+00 0.000000e+00 0.000000e+00
-5.500000e-06 5.000000e-07 2.500000e-06 0.000000e+00 0.000000e+00 0.000000e+00
-5.000000e-07 1.500000e-06 2.500000e-06 0.000000e+00 0.000000e+00 0.000000e+00
-1.500000e-06 1.500000e-06 2.500000e-06 1.000000e+00 1.000000e+00 1.000000e+00
-2.500000e-06 1.500000e-06 2.500000e-06 1.000000e+00 1.000000e+00 1.000000e+00
-3.500000e-06 1.500000e-06 2.500000e-06 1.000000e+00 1.000000e+00 1.000000e+00
-4.500000e-06 1.500000e-06 2.500000e-06 1.000000e+00 1.000000e+00 1.000000e+00
-5.500000e-06 1.500000e-06 2.500000e-06 0.000000e+00 0.000000e+00 0.000000e+00
-5.000000e-07 2.500000e-06 2.500000e-06 0.000000e+00 0.000000e+00 0.000000e+00
-1.500000e-06 2.500000e-06 2.500000e-06 1.000000e+00 1.000000e+00 1.000000e+00
-2.500000e-06 2.500000e-06 2.500000e-06 1.000000e+00 1.000000e+00 1.000000e+00
-3.500000e-06 2.500000e-06 2.500000e-06 1.000000e+00 1.000000e+00 1.000000e+00
-4.500000e-06 2.500000e-06 2.500000e-06 1.000000e+00 1.000000e+00 1.000000e+00
-5.500000e-06 2.500000e-06 2.500000e-06 0.000000e+00 0.000000e+00 0.000000e+00
-5.000000e-07 3.500000e-06 2.500000e-06 0.000000e+00 0.000000e+00 0.000000e+00
-1.500000e-06 3.500000e-06 2.500000e-06 1.000000e+00 1.000000e+00 1.000000e+00
-2.500000e-06 3.500000e-06 2.500000e-06 1.000000e+00 1.000000e+00 1.000000e+00
-3.500000e-06 3.500000e-06 2.500000e-06 1.000000e+00 1.000000e+00 1.000000e+00
-4.500000e-06 3.500000e-06 2.500000e-06 1.000000e+00 1.000000e+00 1.000000e+00
-5.500000e-06 3.500000e-06 2.500000e-06 0.000000e+00 0.000000e+00 0.000000e+00
-5.000000e-07 4.500000e-06 2.500000e-06 0.000000e+00 0.000000e+00 0.000000e+00
-1.500000e-06 4.500000e-06 2.500000e-06 0.000000e+00 0.000000e+00 0.000000e+00
-2.500000e-06 4.500000e-06 2.500000e-06 0.000000e+00 0.000000e+00 0.000000e+00
-3.500000e-06 4.500000e-06 2.500000e-06 0.000000e+00 0.000000e+00 0.000000e+00
-4.500000e-06 4.500000e-06 2.500000e-06 0.000000e+00 0.000000e+00 0.000000e+00
-5.500000e-06 4.500000e-06 2.500000e-06 0.000000e+00 0.000000e+00 0.000000e+00
-5.000000e-07 5.000000e-07 3.500000e-06 0.000000e+00 0.000000e+00 0.000000e+00
-1.500000e-06 5.000000e-07 3.500000e-06 0.000000e+00 0.000000e+00 0.000000e+00
-2.500000e-06 5.000000e-07 3.500000e-06 0.000000e+00 0.000000e+00 0.000000e+00
-3.500000e-06 5.000000e-07 3.500000e-06 0.000000e+00 0.000000e+00 0.000000e+00
-4.500000e-06 5.000000e-07 3.500000e-06 0.000000e+00 0.000000e+00 0.000000e+00
-5.500000e-06 5.000000e-07 3.500000e-06 0.000000e+00 0.000000e+00 0.000000e+00
-5.000000e-07 1.500000e-06 3.500000e-06 0.000000e+00 0.000000e+00 0.000000e+00
-1.500000e-06 1.500000e-06 3.500000e-06 0.000000e+00 0.000000e+00 0.000000e+00
-2.500000e-06 1.500000e-06 3.500000e-06 0.000000e+00 0.000000e+00 0.000000e+00
-3.500000e-06 1.500000e-06 3.500000e-06 0.000000e+00 0.000000e+00 0.000000e+00
-4.500000e-06 1.500000e-06 3.500000e-06 0.000000e+00 0.000000e+00 0.000000e+00
-5.500000e-06 1.500000e-06 3.500000e-06 0.000000e+00 0.000000e+00 0.000000e+00
-5.000000e-07 2.500000e-06 3.500000e-06 0.000000e+00 0.000000e+00 0.000000e+00
-1.500000e-06 2.500000e-06 3.500000e-06 0.000000e+00 0.000000e+00 0.000000e+00
-2.500000e-06 2.500000e-06 3.500000e-06 0.000000e+00 0.000000e+00 0.000000e+00
-3.500000e-06 2.500000e-06 3.500000e-06 0.000000e+00 0.000000e+00 0.000000e+00
-4.500000e-06 2.500000e-06 3.500000e-06 0.000000e+00 0.000000e+00 0.000000e+00
-5.500000e-06 2.500000e-06 3.500000e-06 0.000000e+00 0.000000e+00 0.000000e+00
-5.000000e-07 3.500000e-06 3.500000e-06 0.000000e+00 0.000000e+00 0.000000e+00
-1.500000e-06 3.500000e-06 3.500000e-06 0.000000e+00 0.000000e+00 0.000000e+00
-2.500000e-06 3.500000e-06 3.500000e-06 0.000000e+00 0.000000e+00 0.000000e+00
-3.500000e-06 3.500000e-06 3.500000e-06 0.000000e+00 0.000000e+00 0.000000e+00
-4.500000e-06 3.500000e-06 3.500000e-06 0.000000e+00 0.000000e+00 0.000000e+00
-5.500000e-06 3.500000e-06 3.500000e-06 0.000000e+00 0.000000e+00 0.000000e+00
-5.000000e-07 4.500000e-06 3.500000e-06 0.000000e+00 0.000000e+00 0.000000e+00
-1.500000e-06 4.500000e-06 3.500000e-06 0.000000e+00 0.000000e+00 0.000000e+00
-2.500000e-06 4.500000e-06 3.500000e-06 0.000000e+00 0.000000e+00 0.000000e+00
-3.500000e-06 4.500000e-06 3.500000e-06 0.000000e+00 0.000000e+00 0.000000e+00
-4.500000e-06 4.500000e-06 3.500000e-06 0.000000e+00 0.000000e+00 0.000000e+00
-5.500000e-06 4.500000e-06 3.500000e-06 0.000000e+00 0.000000e+00 0.000000e+00
\ No newline at end of file
diff --git a/tests/test_projector/ref_mag_proj_x.nc b/tests/test_projector/ref_mag_proj_x.nc
new file mode 100644
index 0000000000000000000000000000000000000000..025a3c9090a0e73cbf412772d985cfed6ee93b46
Binary files /dev/null and b/tests/test_projector/ref_mag_proj_x.nc differ
diff --git a/tests/test_projector/ref_mag_proj_x00.nc b/tests/test_projector/ref_mag_proj_x00.nc
new file mode 100644
index 0000000000000000000000000000000000000000..f869a6319e79d9343ecfb6852cf29aba53261899
Binary files /dev/null and b/tests/test_projector/ref_mag_proj_x00.nc differ
diff --git a/tests/test_projector/ref_mag_proj_x45.nc b/tests/test_projector/ref_mag_proj_x45.nc
new file mode 100644
index 0000000000000000000000000000000000000000..3a342a2b9697a427134224a09cc0472742c49ec1
Binary files /dev/null and b/tests/test_projector/ref_mag_proj_x45.nc differ
diff --git a/tests/test_projector/ref_mag_proj_x90.nc b/tests/test_projector/ref_mag_proj_x90.nc
new file mode 100644
index 0000000000000000000000000000000000000000..50665b83db590616600411323503e265c8289d63
Binary files /dev/null and b/tests/test_projector/ref_mag_proj_x90.nc differ
diff --git a/tests/test_projector/ref_mag_proj_y.nc b/tests/test_projector/ref_mag_proj_y.nc
new file mode 100644
index 0000000000000000000000000000000000000000..d84e528d9b677641e946acf26c89c0f0d3d4b341
Binary files /dev/null and b/tests/test_projector/ref_mag_proj_y.nc differ
diff --git a/tests/test_projector/ref_mag_proj_y00.nc b/tests/test_projector/ref_mag_proj_y00.nc
new file mode 100644
index 0000000000000000000000000000000000000000..de3f39f1b426fad4954e793170c8ff84020d24e4
Binary files /dev/null and b/tests/test_projector/ref_mag_proj_y00.nc differ
diff --git a/tests/test_projector/ref_mag_proj_y45.nc b/tests/test_projector/ref_mag_proj_y45.nc
new file mode 100644
index 0000000000000000000000000000000000000000..4169d4a36c37235db201d966162407b6d9bb9c43
Binary files /dev/null and b/tests/test_projector/ref_mag_proj_y45.nc differ
diff --git a/tests/test_projector/ref_mag_proj_y90.nc b/tests/test_projector/ref_mag_proj_y90.nc
new file mode 100644
index 0000000000000000000000000000000000000000..b02f7c3570158c76cb0c5e677ad9f44ab020407e
Binary files /dev/null and b/tests/test_projector/ref_mag_proj_y90.nc differ
diff --git a/tests/test_projector/ref_mag_proj_z.nc b/tests/test_projector/ref_mag_proj_z.nc
new file mode 100644
index 0000000000000000000000000000000000000000..fbb20504c9cfee6a3117dd80ecd6e8c62da3eefb
Binary files /dev/null and b/tests/test_projector/ref_mag_proj_z.nc differ
diff --git a/tests/test_projector/ref_proj_x.txt b/tests/test_projector/ref_proj_x.txt
deleted file mode 100644
index 94ce2a36ef75cec69acef07af19b74d522e1729f..0000000000000000000000000000000000000000
--- a/tests/test_projector/ref_proj_x.txt
+++ /dev/null
@@ -1,4 +0,0 @@
-0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00
-0.000000000000000000e+00 4.000000000000000000e+00 4.000000000000000000e+00 4.000000000000000000e+00 0.000000000000000000e+00
-0.000000000000000000e+00 4.000000000000000000e+00 4.000000000000000000e+00 4.000000000000000000e+00 0.000000000000000000e+00
-0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00
diff --git a/tests/test_projector/ref_proj_y.txt b/tests/test_projector/ref_proj_y.txt
deleted file mode 100644
index b8a14486eb17c76bc28dd0381964ca9e03c0781c..0000000000000000000000000000000000000000
--- a/tests/test_projector/ref_proj_y.txt
+++ /dev/null
@@ -1,4 +0,0 @@
-0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00
-0.000000000000000000e+00 3.000000000000000000e+00 3.000000000000000000e+00 3.000000000000000000e+00 3.000000000000000000e+00 0.000000000000000000e+00
-0.000000000000000000e+00 3.000000000000000000e+00 3.000000000000000000e+00 3.000000000000000000e+00 3.000000000000000000e+00 0.000000000000000000e+00
-0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00
diff --git a/tests/test_projector/ref_proj_z.txt b/tests/test_projector/ref_proj_z.txt
deleted file mode 100644
index eec624bc009a47165dc5410af77e4ac999c388cd..0000000000000000000000000000000000000000
--- a/tests/test_projector/ref_proj_z.txt
+++ /dev/null
@@ -1,5 +0,0 @@
-0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00
-0.000000000000000000e+00 2.000000000000000000e+00 2.000000000000000000e+00 2.000000000000000000e+00 2.000000000000000000e+00 0.000000000000000000e+00
-0.000000000000000000e+00 2.000000000000000000e+00 2.000000000000000000e+00 2.000000000000000000e+00 2.000000000000000000e+00 0.000000000000000000e+00
-0.000000000000000000e+00 2.000000000000000000e+00 2.000000000000000000e+00 2.000000000000000000e+00 2.000000000000000000e+00 0.000000000000000000e+00
-0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00 0.000000000000000000e+00