From fc55842d9d6814795a13793dd1d430b8447b2c62 Mon Sep 17 00:00:00 2001
From: Jan Caron <j.caron@fz-juelich.de>
Date: Tue, 3 Feb 2015 23:35:26 +0100
Subject: [PATCH] Mainly implementation of diagnostics class! ... And stuff!
docstrings: updated! costfunction: CFAdapterScipyCG obsolete, deleted
estimate_lambda() dataset: use_fftw and threads obsolete diagnostics:
implemented stuff... much stuff... fft: Docstrings and implementation of
threads magcreator: fixed bug in slab (dimension mix-up) phasemap: plot now
also returns colorbar reconstruction: PrintIterator and
optimize_simple_leastsq obsolete regularisator: Docstrings scripts: minor
modifications scripts - reconstruct_random_pixels: now use optimize_linear
testcases: analytic, compliance, kernel, magcreator, magdata, projector
---
pyramid/__init__.py | 26 +-
pyramid/analytic.py | 3 +
pyramid/costfunction.py | 71 +-
pyramid/dataset.py | 39 +-
pyramid/diagnostics.py | 241 +++-
pyramid/fft.py | 126 +-
pyramid/forwardmodel.py | 4 +-
pyramid/kernel.py | 44 +-
pyramid/magcreator.py | 14 +-
pyramid/magdata.py | 1278 ++++++++---------
pyramid/numcore/phasemapper_core.pyx | 2 -
pyramid/phasemap.py | 17 +-
pyramid/phasemapper.py | 7 +-
pyramid/projector.py | 11 +-
pyramid/reconstruction.py | 150 +-
pyramid/regularisator.py | 166 ++-
pyramid/version.py | 2 +-
.../patrick_nanowire_reconstruction.py | 25 +-
...rick_nanowire_reconstruction_simulation.py | 2 +-
scripts/collaborations/rueffner_file.py | 8 +-
scripts/collaborations/vtk_conversion.py | 195 +++
.../vtk_conversion_nanowire_full_90deg.py | 6 +-
.../reconstruct_random_pixels.py | 7 +-
.../reconstruction_linear_test.py | 52 +-
.../reconstruction_linear_test_batch.py | 2 +-
...econstruction_sparse_cg_singular_values.py | 2 +-
scripts/test methods/histo_norm.py | 39 +
tests/test_analytic.py | 2 +-
tests/test_compliance.py | 4 +-
tests/test_kernel.py | 41 +-
tests/test_kernel/ref_u.npy | Bin 0 -> 276 bytes
tests/test_kernel/ref_u_fft.npy | Bin 0 -> 400 bytes
tests/test_kernel/ref_v.npy | Bin 0 -> 276 bytes
tests/test_kernel/ref_v_fft.npy | Bin 0 -> 400 bytes
tests/test_magcreator.py | 16 +
tests/test_magcreator/ref_ellipse_x.npy | Bin 0 -> 584 bytes
tests/test_magcreator/ref_ellipse_y.npy | Bin 0 -> 584 bytes
tests/test_magcreator/ref_ellipse_z.npy | Bin 0 -> 584 bytes
tests/test_magcreator/ref_ellipsoid.npy | Bin 0 -> 584 bytes
tests/test_magdata.py | 103 +-
tests/test_magdata/ref_mag_data.nc | Bin 8474 -> 7684 bytes
tests/test_projector.py | 68 +-
tests/test_projector/ref_mag_data.nc | Bin 0 -> 8356 bytes
tests/test_projector/ref_mag_data.txt | 122 --
tests/test_projector/ref_mag_proj_x.nc | Bin 0 -> 7276 bytes
tests/test_projector/ref_mag_proj_x00.nc | Bin 0 -> 7204 bytes
tests/test_projector/ref_mag_proj_x45.nc | Bin 0 -> 7204 bytes
tests/test_projector/ref_mag_proj_x90.nc | Bin 0 -> 7204 bytes
tests/test_projector/ref_mag_proj_y.nc | Bin 0 -> 7204 bytes
tests/test_projector/ref_mag_proj_y00.nc | Bin 0 -> 7276 bytes
tests/test_projector/ref_mag_proj_y45.nc | Bin 0 -> 7276 bytes
tests/test_projector/ref_mag_proj_y90.nc | Bin 0 -> 7276 bytes
tests/test_projector/ref_mag_proj_z.nc | Bin 0 -> 7156 bytes
tests/test_projector/ref_proj_x.txt | 4 -
tests/test_projector/ref_proj_y.txt | 4 -
tests/test_projector/ref_proj_z.txt | 5 -
56 files changed, 1665 insertions(+), 1243 deletions(-)
create mode 100644 scripts/collaborations/vtk_conversion.py
create mode 100644 scripts/test methods/histo_norm.py
create mode 100644 tests/test_kernel/ref_u.npy
create mode 100644 tests/test_kernel/ref_u_fft.npy
create mode 100644 tests/test_kernel/ref_v.npy
create mode 100644 tests/test_kernel/ref_v_fft.npy
create mode 100644 tests/test_magcreator/ref_ellipse_x.npy
create mode 100644 tests/test_magcreator/ref_ellipse_y.npy
create mode 100644 tests/test_magcreator/ref_ellipse_z.npy
create mode 100644 tests/test_magcreator/ref_ellipsoid.npy
create mode 100644 tests/test_projector/ref_mag_data.nc
delete mode 100644 tests/test_projector/ref_mag_data.txt
create mode 100644 tests/test_projector/ref_mag_proj_x.nc
create mode 100644 tests/test_projector/ref_mag_proj_x00.nc
create mode 100644 tests/test_projector/ref_mag_proj_x45.nc
create mode 100644 tests/test_projector/ref_mag_proj_x90.nc
create mode 100644 tests/test_projector/ref_mag_proj_y.nc
create mode 100644 tests/test_projector/ref_mag_proj_y00.nc
create mode 100644 tests/test_projector/ref_mag_proj_y45.nc
create mode 100644 tests/test_projector/ref_mag_proj_y90.nc
create mode 100644 tests/test_projector/ref_mag_proj_z.nc
delete mode 100644 tests/test_projector/ref_proj_x.txt
delete mode 100644 tests/test_projector/ref_proj_y.txt
delete mode 100644 tests/test_projector/ref_proj_z.txt
diff --git a/pyramid/__init__.py b/pyramid/__init__.py
index ad7b174..47aba58 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 5faef31..3585c89 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 50c41ff..fb67e45 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 9f422ca..a1924fb 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 fbec1ff..b9d84bd 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 1c7733d..5baea7f 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 f2758b9..d7f0e6e 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 8c322c5..b8bc451 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 c3abb46..3083d45 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 7fe3e6c..53f4d43 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 d7b1f04..59be7c8 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 d60b410..ce5c45a 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 76fedda..b4c6cab 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 48e6944..773db4b 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 d9386ee..f5df52a 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 5c115bb..5f02f24 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 e604916..81491db 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 deb5d33..80c625b 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 3c1edba..81efbba 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 230c4cb..82b13fb 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 0000000..62f21be
--- /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 d1ab275..2f01c98 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 3aec76e..7a12148 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 7da8b7a..264625a 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 3f310a8..2de7ce4 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 0ccee36..df39915 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 0000000..7a9dea3
--- /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 562ee62..dbdd9d6 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 6ffe674..b92af3f 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 39cb1e1..cfc0bbd 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
GIT binary patch
literal 276
zcmbR27wQ`j$;jZwP_3SlTAW;@Zl$1ZlV+l>qoAIaUsO_*m=~X4l#&V(cT3DEP6dh=
zXCxM+0{I%|Itu2RItsN4aKN>mU3140uB;s`3Nv=BXI}}XL41(71iQ?Rw(i;;tk+Hf
z>Az6A7Ay`F+s=B;3oH&)vlAlk1yZw}0S^Gmn*;S(0M!F=ttE&C@j>E1`#@qKHdr1)
Z*IIzYf$p&Yy2}FSP9S?Fm<I7d;s8vcV9fvk
literal 0
HcmV?d00001
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
GIT binary patch
literal 400
zcmbR27wQ`j$;jZwP_3SlTAW;@Zl$1ZlWd``qoAIaUsO_*m=~X4l#&V(cT3DEP6dh=
zXCxM+0{I#iItr$mItsN4aKObtB$&43`7R)~$_t)s1=1pH_kITmgV@tU`Yf+X?AQ*(
z<~&??EkR;&O1pQ0Fo^wc;SNiwCFd=Gct_eoZLrv`ORIK*Fi8H8^?Xa9Sr#DkIOa2f
z#h%=YumE9@JTU-jo+ZdEka}eI0nN08nzw^k_d)%;3uG3^&oK7^{R(y;)W6#Sc~x}<
literal 0
HcmV?d00001
diff --git a/tests/test_kernel/ref_v.npy b/tests/test_kernel/ref_v.npy
new file mode 100644
index 0000000000000000000000000000000000000000..c854a27d8debf3f664bb65695088dba2ecb41630
GIT binary patch
literal 276
zcmbR27wQ`j$;jZwP_3SlTAW;@Zl$1ZlV+l>qoAIaUsO_*m=~X4l#&V(cT3DEP6dh=
zXCxM+0{I%|Itu2RItsN4aKN>mU2}&7yUY&OYp1p|00EFT2l6d|;ucG|vUap}*8=6d
zc0lC3EP#BFxJ8S?43Hk6{9mvfkhTEwLE=EYE5YJGwL3xjfV3q@9w?4iA3po=xf1|X
CgJ8`7
literal 0
HcmV?d00001
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
GIT binary patch
literal 400
zcmbR27wQ`j$;jZwP_3SlTAW;@Zl$1ZlWd``qoAIaUsO_*m=~X4l#&V(cT3DEP6dh=
zXCxM+0{I#iItr$mItsN4aKOa?1=E&1-vz{0TwHf8K|*^jt=b8~FnJiAA3WI##Ew(i
zy%U5%?8Dac!FoV!7!a|0zXOCp{C^8~SPI#_w*X==AEpLI&j{(Wlv;A$0*H5T%xA)|
z6Qpi|#E$JiY@V`E8!Z3iUW5e*!_0uuK(oN|Q1ig@*vtgF%>sl$`hf1+3340AjWBbN
OX`r8XLH%oq?pFZnL|hsG
literal 0
HcmV?d00001
diff --git a/tests/test_magcreator.py b/tests/test_magcreator.py
index 137b6e6..4abcb95 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
GIT binary patch
literal 584
zcmbR27wQ`j$;jZwP_3SlTAW;@Zl$1JlVqr_qoAIaUsO_*m=~X4l#&V(cT3DEP6dh=
zXCxM+0{I%|Itms#3YMBW3bhIlz{N1c03*!3m<*T<J{nUO9ua&7!Q?2^4l{=CGyr4A
B6|VpQ
literal 0
HcmV?d00001
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
GIT binary patch
literal 584
zcmbR27wQ`j$;jZwP_3SlTAW;@Zl$1JlVqr_qoAIaUsO_*m=~X4l#&V(cT3DEP6dh=
tXCxM+0{I%|Itms#3YMBW3bhIlz{M~`0V5+LHCzrfhH3^;!##8}0{~;k6|VpQ
literal 0
HcmV?d00001
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
GIT binary patch
literal 584
zcmbR27wQ`j$;jZwP_3SlTAW;@Zl$1JlVqr_qoAIaUsO_*m=~X4l#&V(cT3DEP6dh=
pXCxM+0{I%|Itms#3YMBW3bhIlz{M~!0V5+L*yBJN!5T?H0RUsi6|VpQ
literal 0
HcmV?d00001
diff --git a/tests/test_magcreator/ref_ellipsoid.npy b/tests/test_magcreator/ref_ellipsoid.npy
new file mode 100644
index 0000000000000000000000000000000000000000..8abca7c48658e9f9df018d7869ba2dc4461b1d9a
GIT binary patch
literal 584
zcmbR27wQ`j$;jZwP_3SlTAW;@Zl$1JlVqr_qoAIaUsO_*m=~X4l#&V(cT3DEP6dh=
zXCxM+0{I%|Itms#3YMBW3bhIlz{N0>03!&19E`~z#v4#wU<x~dWHJc`BI$xVj2P_$
GsRjU=5*6nF
literal 0
HcmV?d00001
diff --git a/tests/test_magdata.py b/tests/test_magdata.py
index 14d9660..a178de3 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
GIT binary patch
literal 7684
zcmeHMO-vI}5T5;MDYZyZ2#Q2jV?rW^kkD|z7-+W3pQVj07fpaHwvi^ZB~pwuBzRFH
zM-IkB6Azv|nt1VKOvHHe=EaK#Zzk%&izj@u?@e1CG{qn4Pj+BCJM-qv%=>2dy)xVx
zi;Oh+JN&B82Si(f^y6I9CkD>Enls~Lk%Ss}s-ARonnw)L6q)-~P|&2U6LiKj1!N17
z1831r0~K1?57^LR3<1w+7(xdiw50*43L~^*M4U^fhJ!)8JIJh%OwIr(LxAbj^b05b
zrUpe}A-FqX&mwcuTG4aMi7qE$00Y{6(7wh7J<eJ*KM$lI<@Ot59<5NTzQP_f)xgV^
z#lSG5;bamHGZiMTvMVw^X~mM`@z`wA3`Z>pT*6gqYzfYUSa{NM?1zQ@)O=>iV86DJ
zOBc+@h>^-<)48QgK4;7=7mQrKP#(HwU(MKa3+Y~isGoW9Zgzd(kU#|_SScjpwIT5x
zLVHc2VmhT@ntcF;qkln^YaZh=;@Z2SCPl6U=%$-J@W;lJK-?TV&s}Yb#hWyfKW^S5
zC@VTWe4gSg16{@HF`CG%gwVzQ7A^!<apLtqSvH?|8@HD&*umlZsj?+{US`W>TlI2R
zRmvvn1fCTv+cWVH-BRv5F5I&1*n;@W@LZs}y=?2hQns1i#;@hF_4AD$MtSvUbEmgW
z)CoM-v22UtAzB+$Kh^zhi(rPTt|*ifZIUvtInK`+W0!Nrk0+3JntS}R$~IS&75$h%
zN4TcVwfdRLt5+1@t`oQebo1kUBYd+3G55CPT!j`5`sI9f@w4#f?c{?r+L*Pk<T8cj
zRJyx?o1A;paL*~-D#{u(vonyL$mIhtI-Z;Y+|BVM!C1HkJpIWhaomXYM)twT&#p=t
z)s->^uawc}+vuB%qh>scdMB8$K~dlav@Y5`NKe~7m74TqPdSpF_8a<2LVEHe-m)EW
z1!yHb`Jws5>BSIml<dhr;)6u?lqc7+&*svdd*dSDBH$w6BH$w6BH$w6BH$w6BH$w6
zBH$wM-ywjV2)TnD5UUI<U~e8f6n)sCNLV9Q!it&Jtka<wU@ZacP^9dFy_7DLTNR~F
zJotj_&u!ibFn0a_1@%M-*7b*j-<;k^-6A`^4;J2fztC3S$kb6;elYd1zLBYClfOOl
zj*$Po{KgOfUPQ7?`f6OZm*r~n{uaHQuPoQfQ?{4oTJ!((xSX#n*UC?}m*raX_dfm|
DR;BZ!
literal 8474
zcmds63s_Xu7T&BQh&&WTuZTJxzz2%hsVKdYFgOk=h$O?Kv@Zt4koO2vnOII&G=*M@
z-=|;R^w{G{+PzI^QYv~9$qXgRuK7h%NA2Mit0_e7qL;hQ%-MR*=zLOX^|j&5-g}+3
z*Iw&iYp=a$&&p|26a0Nbe8AhA(9$iJkWeRulV6;R8jzDRZE7lvJJ^0_s`A+mVGbU~
zOJ=1aTx%a79umvIp`;g9j0bV$pnbo|5+8$X1MrWO%_}E?Y#ecFsU^|AB)Xd+1U#@v
zoLFe>?2XU{1W2t|<dVu%XC=vWmRT#rCmrTQDNT{uy*)f|N~ds1O{vhD86eFQ7&=_0
zuw;-o8MFjBJRP<%KGL)dnLo@ijMO=+%NBdOw>VN_Wr!CViL>cu#2XAG2uYF|fX~$F
z$@-aw=_xa_47&Ikdg3=*vY8>C&H*#yll2aMt<0Y`Kfj#C|7^Boo>@0Ffz8b?o>x+y
zUs}TEJY;4|O3hbZJ)B*ZpPf@QZ#+v&G)g&rd+DWXLqvWRMBe>TIC`Ht+u<mr(6G3H
zAUT8JZ?e<pOY>rJQfCZx_MUfEgp8(*9@K-EX}VK%(~_K`ygZ#lE=QT<Ff3ld(x*o+
z8k?SvW}z3t<hK?7Esc<5+-y(k7bU*;r|@H<46}<LLnM$1`9*&8O7z}*#g8p=26N{}
z*>S{kq2$NYa##huzk0#%?ni0C5@0+$XsYUIEvGd08Y#_3N_NaxrzHDgT(?r+BSa(3
zg+&f)x^CU4bLT4gsE>2lB4<TN0&V>I`cSTE!xcY{m$x%_ev}<YJQrmZXXn+)Vbw4u
zwchPI`7QAsF5~fY@g3qEcQf$4<<?78SMYsC_M|&}Y3yOXLwx)9cIwq$|7^cM1itJq
z!}mLR!?+pv&TV>ijf5}VF=P#;2*pg(EqCc#Ibfwcg!q=p(LQVcj#5V+Dw5b4(htuL
zc#^qtaF{XOkVXdMJ88!qgpV{(Iz!R*Kx--RX87O-7AYkv`_kOXzU1qX#RmfK)}_pF
z1gyNBJKFe8R*p~`Xbz=xx*O8DDV_4uSfX^AbtwjvPUX-;D5X<5au8mn59OmZl_Q7v
zRs1q6m+4d=U$j*!u5$dG<jZubkG~swUl%&nr=JU*%AuzxIxp2H(1oAMaU8zz9p$HT
zWOkyZ3Q&Sey?bocBf?i{P>WAP|B5(qRwvlswlRBVN5vY8(x7!z5p}W4NJ7}GgwHK!
zD}mE5bIY?oY~ul8aXJXSG%B8Vch-f9F=1e>Yk|<s^Y@UNFc?u60Pj^<-<sCs53Otc
z;oz1M1OLBTu-CW2`>S&YkTK`L-e9Al<<RHQ={bd#TqAIrOgOyr>k-;216V^cRdXp7
zZ=JI5ALMur*vVYCJV#D9B(dX)z}~wIK7Rf~KR)~%aQoErDrJ5W7Ie|1Ednm8435=R
zj_}?%0k~kb4xiS&^MTM^3oD)$Rr7vUmHDIicslTYYJBwfUfi|rbPL!`2UPZoFE$hj
zuNPjW|B2Jf*RmDqz{Lcs<SX^O8SOK-?R}3OiGF$4#(axy0o&XGHCEr3y;8{sm4j)4
z`Dg2Pw%`j5Kc)ywqzpQ@Jp80(Ksnf-FNJgP@or=AMDQCIra8ZQ?g?XNCs^OIF&8F9
zhx2!Z0rz2x+qt@6c)k7p4B+Eq^L{oayB=C@srCadAR9akF^7%(2VgyHW5BnhuQopo
zxOdFJ%zdSR=qn*$iu-k;oVT~{cax0o90k*xYTWun?g(U;q=Us@13rzF-;jOXV8uA}
zEq~b{_W1x>%ac_$`iJk!;e!Ma0-K;O>L7le0VdS5U*Vi(7aqs}OKzwp;AG({{wV`q
z5zJPNRq6w5&G+Jum0(Ix>pyT?-UDp^V&EneyB({PQ+a69Yqig$gCDNzkiv$MBK!p|
zh-tuA!8(3q5m?8UspJfc+5R*sO+%g4IcNCR$0ryuA6U^AHYxUy@v(F;_ExVK<-Fl9
zc^b*98DQ09tK>g8a>Z@|*BXawZB)?6eYF$UfW5`WjE;OGp1;%r+=&(@BC=_;5%Y<C
zXeNv~`S=KeYf6|D44E4iKW0Qs0zNTAWOdO%G3h99`)$m)rRG=ImFT-YHYO_RutnUI
z0jypfSGMDc`$#Ir=S~~LMgCAkA~Bx5)VAN&7Fo^LV%`|0&KuE<`Hvc3p8&$~G?gt%
zdwbHJ6uy|_5tCl1TTJ1JP7CH$>zi#VnUlNvc3NNB3+qUeiq5{i;5<JI<GlO`Xy<gV
zBn4?;#W+{aJLM(IDFJOaj#t~R<WaU0d<#kiT-SD7*Gjs(b|PX+y3?=1R>_|nKKL1-
zFayLmwGJsYqn2`!{eU~t2KOXojxsJc07+2ur0T-rjYB5_i#VUmHgkM^8rX49Ns9?Q
z&p(R!Iof?*%^WiK7;8mL4yUVZ@vj;CAF&uQ(Z%M?-uhK9{{Awsn#}NQ#OnRVMQF=J
zwaiIN+63cSq`Syb+GyNX_{FjGF<=c;$97I+Q%G>fzrbW|gSqREU({lpn|w5o+xXQE
zHY)?{gR|l5fBv+~8`q;XRy}WC_2$1G+kyL+M!oLxJCgVCsK0n?Y+gb8q-E@8%qP1u
zR5mJgbJw2SHt>YF{|uN8s`HyWeq;TrFNAHwAz()tn7eGH>}))f?65HpH7C?t&Y@k+
zYMkB4S##Hy>$G#qz^*S>>0B}7<6?G&0oHxpvo*D3>-0dr=?vIi^napj$WH?5Xw|9Z
zxA&R6C1MQ=&n<Ls@A8U*Q{qYkaF}PFJh}7-vGELWF7lOjWp25DJuCTQv^s8;JoQtm
zI&0S$K>Si|*S9MTpNNTjVa2!mRXD5)N-0=4x*NFlHfD84T5lu9rzyT1YIG-WVewpN
z{LBn%n3Cj(dd$@>v7wxIZBVp@obiR$JJoTe<WaWQqwT`jb70z~#^;5Oq}Pn-m)3#m
z7~9w$)p{@PUmV6BUwC{Dsl_-9QOCbh{@9Es;2y3U>>F&%%N3mmwbJ<^*8QHfxg*~q
zbfCRG+^_FmvTM2#Wm-|@Yt?;QpVMJ`T+i~DH720Tn1y?R)yBNB{`^j|DIJ1we^lzP
zY;W)C9zpO75P)ZZH>)c*kwFuI8>HSNl>Ex}*{Cxi{KZ<}9%+M}y8DW`{Jjv-sg409
zf5W!M-2#h#5!8O!Taj3EAvX*z+@anV_*0sV!Y14oTGjkYnaXxk#SbxjXffDtFNSxU
z@AD;i4zc4|{Xp}PIa<Vjp>7$x(4M#{BG3;^F8kp7;ew9^_grm_2+b9znx~B~0Q()q
zaLD(OnPS`+un)<E7yq6a%;NdM>Y+XxBv+e931cz8-Q_+8K1$lKRmjf)9(7RatZe_g
z{5wv#gm`XK@0&`VW9=P>L_D9ceKhdaHxEW!sF?t*%@g6%uEs4ZFvnVl6+uc;Yr?`j
z{O*H(QI55td<*V7R@`@#bY)wc^q4pMKqqk9dxrU!CwCrYQ6Gz!8k3m^f@a-QS_$IR
zYuoF%d3Z|tsYJSY;G~;}RDFUzRX<aw&vM*6+=)6$-&S+8&DrJi%vY`;F8?w?yaVuS
zsq-5@d(E2(dRIryrFToTPwzr#fBU#C$CPON+1oe%{ztmMP~P`Q_Z-T;bRQBLH{Ij6
z^*_B`qW-6S>VMilKmO=_)aN%I8|f}ZiG_6MqU=l29@ns-QJUG^uU>w?*_xp^P|Rqb
zVn+LKSMI)nH}h1gR9oVl52d;GDcxP|%J15DH|~!@ch%Rme?2{2`Ca?h8~=O1ch%Rm
Re?9$N`Ca?h8{hNy{|BS39nJs%
diff --git a/tests/test_projector.py b/tests/test_projector.py
index 6da902f..b328c02 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
GIT binary patch
literal 8356
zcmeHM%}Z2K6hH4{#_=`gXo_WmuLy(+5gc6@q4MVW)ijbbn*x=w8FIjxksKoj!8Rq?
zt{{Sfg8CEMxQZ4+yJ!<pK`q(^)xu2>*gfywGc%4p^c5qK@94~V_uTW&x%c<GA9%by
z6dI_}nl;(y1Ck9P&vb6-5FIz`)tb?v(3tFhE+4mSibr(OE9F{~K|+_t7EqZ|Bq;Zw
zw*L(JDIh}w`++#r>Ak?m)Ca)>5Zq7zWQhs-F(J)cC;9>bd^*UyPz<gBXoLXciSd`#
z@N23hiKXD-gg?v7MQOmub;nw*0R!mJ<b$Rq*6Fe4qG&vkdW8Gm74Oju75f>ULRSU6
zY+4KqBitX2!XakDfy?X@8a{7EqQj$+$*9pEHo<=hXQ{9$ST98S&zqKipYV@Or{;9_
z8`m?*oDmw(6RC7EGndL{^r`urp2_BlUDx7QQ}L;pWV=q%$)b4ITEB2epac>x>p>yC
zGbFxx^j9Pz^Ds7<Q~(r?{slqTG|vOvdS7so7glJahh6B85PLUj44vbl)^NEpF}R}T
z=Cx)fskHb!`AG)4jO;NPNi761L}iKVhI^-%eOqaxl}_3@ZPun>^QahYJ~6h<v;~?i
z^x%(ZTfXy2E7G<o)~pR}f+z4SAZ@E+5#5p=rYlcd7rg0uzq*w+vCzbE!*<b*2T=KL
zvktyW()L22x|z0>zoBiSy?VV!TPHv0VURC`^4%tQ0?&1%ZB{Izr5^cP<!K9onICyo
zpe=+WlvPQhG0`DOyN>AGv7ghdZ*$LRcmpXXd88kUjBad8>QRo)6OoE`sr{YGkFvY~
zcZI<1rU&PGU7P)dn0uS>U5O?O>g8;D_PdDZZQ-5NTb+(y$)s}giDX+9cRBmC>VZYw
zBDNK%Wp^Mue2N!fcsLpZJk4=!V=kNn-u@JmG#<oiHTz<su}^6m&82NjUfRZ(Z^N%I
z4jQ9jv|Gu94>kpEK=Y#0gY0zd6V+rVe;Om%>Aaz(LUxKnONs0hM}Srs55*xr(zG!I
z90h-hkNBXFKgH)(_FY?eavxj-Tm)PMTm)PMTm)PMTm)PMTm)PMTm)PM{yPM)5>Z}L
z^NU&fXRtPp6^ah5P>h)a=9n2V%t@<4(alN%SfNP7bMd)ku2`xlRN}!GV1I7&N`SuY
z|1YQ}La?qp?)hoeM)oYTRr_G+&CfH9dz+amD)WzGU-vdM)oilM6CViqe@ox!1t|T{
z=t^z)Ok!uJ{@rT3I8N=m&9`6uoZ@Nc?-W0~I8N<O^X@l)r+C`=JH^c|j#InSyp{9+
E4ZoKf2LJ#7
literal 0
HcmV?d00001
diff --git a/tests/test_projector/ref_mag_data.txt b/tests/test_projector/ref_mag_data.txt
deleted file mode 100644
index 32d452b..0000000
--- 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
GIT binary patch
literal 7276
zcmeHMO-NKx6h800nejP}IhJCX;Ohs3iZc8`jBsRT{Hbw<90>|cre@NC`H>tYMX*IA
zxhN=zprA#YHtm9N6)jvxn-)dTqE$h)jUYSc-Fx1QGY(}jF5(>>&pYRyd+xp8ci+8w
zyg3vYsC75HRb2;`bt$P9;c1r*x1ZPaMuq~Ts^^)y-}Grt*~DFnwOa*+LK<7Zfl5=L
zID&1@35wG|g$DeB7}OhmAm-$Q&Kl@k*8o%nh2l_H3+tfI>!qh17)9b23V>$_FdiR&
zX->bkN>OkKkxs_r$U?LR%3Qy%)toSZ0ZlqIEn%b6T#M!w13C7H_&c&6uTU$$!ow7*
zfeV*~VE9A*(J1Uf7foEnPhfZ?7>*7{!V^({e<%o^^R!A0mtb}X_m2e4c%O_Pn@Y|Y
z_#HDciJU($V8oN@L}n(L%@`MFb4Dhc%Qwx%t|Vg@Q;DMnO9y6gt+jp>kU#+>j4bQl
zL!$4*WXLX)NnDdu02G1#8A;ca7zc#+p5$afdZCRU{zP9t+Y|8*offIq<P!YcQ-I;x
zO}64-Xq*cZ2)at_p`1$2cT$SV5_cTtK0Q59N}Eee0{54Yt<z>M4m($jHeJ4r8)@@4
ziwXY~ZIhEf-{)ytl#4Gzo8$?c^F-UK+(b8&2kFYw)(v-$EbFDT@d3qp>$I6?StV(E
zAyM5(+sePtrvJEnD^FX8IOt(eYy>{s`E8Rtf%7WSHZ74_>Q%o~p0+M{-Su|WqAfr(
zR2_=KbK;YL-#Uv%xs<Z);+`=!xo6yT18IlE3G{6d4_C4k$6k?gN~YwuwEUeaj<N*_
z?uwidbo0a3b=|R^Cxp4HiTWyh&V%DpHa-1S=5v|9Hz*Entxm@#Gs)a+JkeGqLQX!e
zx^GgqNVW!bcn9Jx&R>9`;piC9X-;<|jYVsq+n;l?CM+D)_yfg_pF*45g*FXdXj9HN
zp;zYy{gDv)&B0_0nG!ePb@A!Jem?ejE3lvAd5-Mo^M<zw`#BG9A@*|~37Rz@&cku6
z`NfcMSn-^n@FC%N&M&<9v(9?5A8Z6{1Z)Is1Z)Is1Z)Is1Z)Is1Z)Is1Z)KUI|Qf_
zQCz|H$W?k$RGX&?MLSg}MuP*v(O}pgoG>dC-B=Q!3Pn7Yi_Ii*`BH^di3i<_e{PFP
zfU)WSFSsVcu&zAn{chGqwrymy_JP9(KcpJBcVw=p%sw3ZyuBlH%_h5i;T<E7+9#fU
pQVR>+XkY7NE3e64980YfTkN-N>{BArXT`_>vMV}n$wVsY_yt8v?Zp59
literal 0
HcmV?d00001
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
GIT binary patch
literal 7204
zcmeHMO-NKx6h7}~#&H^REX6XxR|$fJ7!DFfsLVY78Xa;ZC{P)jK?j@}$x(6;EQ%zT
z1ra2)Xx*+&t7;*%Yu84C7Oe`ZMVlb9bKbpY<~cf)#kh!f^v!+u+;h*p=lkx@c)U9h
z>Z|fLcvX)FShl1zHwsO=Y>515y%`w@jjFz9>IqAyxnvV}DK5P#C}dLC2s##;0>u_w
z_njj<4OFPbABaJX(G6lwJrFnmfh`R{RZz$dg-cOu?(zHS?I2o_1VsVBGXxlqkH4~}
zUsa(fID~K~?QvuwS}l35Io@PV7{Gvf57e(<qsxj#^NNA=qr(24?8gymr6W8=rW&|0
zECeGM?u|y_Fsf+cDn6mXVRIxp7#W$227ALM_^wct8irtX80j50E&Cp6KQ@z`H}E%R
zq!Za-sLzNeQ;GC^GLtr@7P3Y<lPxqY#;zq}Q?rRSgQWw#xHp<U2uQ#N2_wtp?;-IV
zLTAV>lSy2YlmQfh{zXaGj2H)m_JQOiC$-SRFMm?rzH~Mc9JnA{ZO9OU+*5$zK0<Xx
zXJ}lxNd#Rb_E3hBO966GUgFNc;>V{achcq-lfd;QWbd?D!C~jJ(dLo6aXW4P1~K8k
zqHTKm=eq)J%QE;9v`L=8wM4W%my76*@+eh)+B)I>=~d58+IX*Uxp&&EgRh*ly^yGG
zr)}+DXq#xO{8pf?Lmc$bFBSqH?!4P1PvE*vw9Uyyw9=)1E<bGnFmpq%^R$I%hN`Y8
zJSRR0c-Q4=l%bSt7x#>@%RS?z8%S#vp1v;Q;YxKyKPFrn(k1WG(s!yj%5oCiH8~^b
z<d?mdZJW2&33GQn^;I}n(644vb6=!CxAnb2(YduU6`M{cvkUP=ONB5w|D@ufMcp#h
zHK@ir5O3td1sEQTjsczKw6<w1ih*u__Q@q-p;zJuiWeVyoz(U^4Q{WK&$r<>m-~Z}
zFzKzqqzzRiZoqNz>A`wF_IW9=p6z*#tmpHFzml+?{qPcEJ^PWM<>$kG*p5ry7!r=W
zJ^LqoNZ6kJ3oU+D=ii(c2LT5G2LT5G2LT5G2LT5G2LT5G2LT5G2Z8?%0m?)aXRv)T
zO5ZHy<|#wbP8o_(v(FqgM}p?0m7(awlmKNY;<0ROK9Mb?D)O0l@c8l1ZIKBucK!ba
z=R_FRwZ~oGt=!1IjcnyUaCr0mY~B8j%o&x1^|4R;J2K~NGOH7t05`fZeB1L9-$`eG
kNN3j<*V&fzV0&0^*O3kDxg?!kUtG7_L6Lu6JLM^V0k$FpP5=M^
literal 0
HcmV?d00001
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
GIT binary patch
literal 7204
zcmeHMPe_zO6o23M?XK&#{wu{&!Os*V3o#T4!*-Y5&ws6Mxe^trthwlduA-}?DH8p$
zyjc)ILP4hv9fF`dRfN<b=+;Go4xI|JLzf`3c{B5NZB@!*J;XOUJ2P+IyqWj=&7bvo
zt;yGr;V5vZb~{j7m8q;qNUCH<;Csb|K$EXsbv{z}a+zk8T{I-Qbf}=9N?s9ISfMG9
ztU+z(Vbs$=g<SRnA;{5PAlBrDnjKKHssX49OQ^>Zm7-PV&d$bKCey;?6$JpT5TLWS
z^EqFBMw+6qDTFzx&t?`$%SFpMCB=Nf06OH`A%BK-TDdQpLri4ZBlNG!apa-4xrh5v
zRRbGyi-F<wdxAmO#YDJpm1SRZtI-y04zzU#y&k^-&eQ0n#$4b7+B~fW*LO?(j{eA~
z&VECBG#vB#8uZ@CU^qG&8H(yXV=+BC6ialChejfyo`G<sPN|w%u`LyU5RkwcNC;Uf
ze-DX$Co{&_B{Gp~Oep|GpnqJ_)i34&A-yR%nUqo}r^=s{mrfiCc$<z2Q%lkXFO3vn
z*mh%Eu@E%Y$!-R^itMq}5}Bw$6R9QcAdG)}P`aKrn^*+fUm|UtHtrl7oHE+%@-+TR
zTXun1@L$o^*Z1RHg0?B?{3f(Xp1?YRv@Oa_bVa#6n0nf3;l_bk`+C~wtWddi+W5tn
zO4^=ERDY#y?q6u@s!acypsiY5^w20a0$uKO+9Xe4J%_Xn%S|*>r+!X7Z8cy_w!Dhd
z=ED`L7Db^o(M>?7E>5F#Wz&A~&gdJwGY))#v<hMA%M`m@v8`ANg-KnSq*J>2n<}od
zNeS+pToKe#W$R_-^2G(j+?J1H6>=6VXNLxdSEM}~|K6Zjs5gBu)EAA!#(Kl$X+ov`
zUfL~A-4wPp$YO6Gdyx|lpuaiT0eG9^-p0Am2Ymg>CY6YVC7pe+<Y0Mi8^vqeIQiN(
z+I;JObE?rB@S~hhCUvkWaRc&;ZV!^vwNG1t<fKn)BstwT^p%9<WQVp8$;plcExsPI
zLwZ!w!H{so^~pZsgNgLXzL2udta#0=m<X5%m<X5%m<X5%m<X5%m<X5%m<X5%m<arL
z2w)^4IfCt!UOERbHjfdCDvVIH8x2Og(dIR}d4!^tg#<7{(Hn||M#Hg0s3IPT2YWXA
zb6Z3L^bP-iK`{}6b?#2xHy#_=wwZbCgH3O|AIRI@nJJ<&w$SluduOJY&CqPua_7kT
zBZNNIjbry+CDS=PM@l`e$1h&Gp4=^TV_uwIabc!5_Mr~x<Je?-i?1KX?PEQ2WW_cA
ntd-}b`KzvaPlmfW{|wKUPwsPLrZ)DW4(a3AWP7ysRFyvgVAvz*
literal 0
HcmV?d00001
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
GIT binary patch
literal 7204
zcmeHMUr3Wt6hGhhZF8=fvlPn&ze*4+#BfMhvdQ-Ouep(vpg?7|p#!%`PRSuy6iN10
z5J5r@J@?j&Ptk+WTW`Hc&_hoJ)k7~qW#@eNo^AdaWicP(J9>A&d+xdC-t+t2@7{U5
z+aK($ay7VArxRE<r8GARO}p%fd~dxO=?@O8?&s<WQ>WQw7Z2H9yHrrfq^=QcSZNAu
zx8Sz>9NB50LM{G43~F>Ah&B1adjPzf8i1;xkR1xwqSfs8c<AXMT9E`q1;8r=7>$p<
zG?!mhp(r?oa3}3?W+7TFd9FF$WG)y$hk7T}FJY(Mj74*aiEKxO{e3x(Bh*Sqc#2Fl
za9~&nMj+G^jlyA6(ZW?M!GR$o936;+$D@IskOA&16s3kCm;=H+LxySZm-Zu5$yptL
zV|qG~4Fr4jcrukp&n7czePS-Fr!(0?*L>_+GBz=tXwz9b(2HZe>AipiiXdTRx&A#Q
z&O_)7*<~_`Ymzd6BGA7m>6#MrfY3gaoaCeyTKM5l>f4vjMgsj8gsXKKLV!mKFdRpy
zt=JeEdu|*-SBX88!Q_IMT$Go%Gcf<**~y)>Im9Axe+k(;ZDw#dxNNjJ<=eQGHcx|C
z@L$n3Ir;N#fwn~%d<ohlPhejl+ScSIx}!W!m7lgwcyM~zxsx`2SGe9gZRWvOPTF2b
zRJYQ$@-MWFwN-vC(AFUidgv1yfe&|n+ayn5zfQEx$W64=rG6?uZC)^PgRk<m1!;w<
zO;LDFd=l_mm#0yNQnFv%Gx{#~jEim{tyOsXvWSN(wH4bj;Zm0_`7JGdr;4L2C&68j
zD}qjb*n8c!ach+@chu8Zg@Xm#)l6#Uv-IaMzc(m0?yXG4Cez96Ts+ZIAxyfTRy;DP
zTcoxI)p!Tuja;|@Lj%zfpwpc8HqAvb(CyDYxh5=ZmH2_;!cyEOb#a>}FK(01H=);;
z`vQ>=>CMTc4K*chz;W^E!FoRSc`LA<?Rkx?=ktcYlCYlr@D^e{`;nmK*Ta6;j%$7~
zBpi8r_D}ebus!=1TKufeKUohJ0u}-m0u}-m0u}-m0u}-m0u}-m0u}-m0{<NXl!@4$
z!FJ0i-P4qtrwm0qWhjP?USrq@2aItuL(z#T0m@LsW7*hjB3npR<TLT$^x&V{A`_tR
z`u_{gi7>1yPrAODxsiRF+01?5^v1jCy8WG*Gb(edBOmv7X3p7Umd7>-`PP}W^IY`W
sd))vxd|Z~iU4CD~4<8jY^;6+C^|2k*G?x8Pjkcx~>&oGu=P5`13#SbgasU7T
literal 0
HcmV?d00001
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
GIT binary patch
literal 7204
zcmeHMO-NKx6h7}~#&H^REX6XxR|$fJ7!DFfIP&KC*W{2RL4nHH3_9S<NRE<&U{NHw
zEQlbXMeBBLT15+?UAs0Cv}jdOE!qT;o%8NJGvm*o7ULq`(Kq+qbI(2Zp6|Oq<MBa%
zsJF`7;8i^yVA+t;+$c2dvLW)T^-iQeG_3kws;6w7=8{d^WozkGK_Qd6M$oa)6xeFP
zb>9WD(?Eq<{DBzM7~LS|6okMb2yAEos)9mxC|rtKbI|Xnx5H>fVu}KQX9zGFAAMs_
zzp6q}a0uZ}+T+MVv|939bF#^vFn|H|9;jc!MwcCn<`o0!$A$ev*^eXCN=JBxOf_(0
zSO|t0?uka>D5_}UDn6ls%hphIATl%_HG9Gq_^wlw8irtZ80xuf+4ezcKQfh^HSjlP
zq!U>))N90(sYH4<nMoTHb6F#u$rhUCV>gqriRna}!P0?V+-psr1SC)d2_wtp?vQwn
zpfhBb$t12x$^eQ$|B9q*N{j<S`&e?4lUiuum%WttuAYyW{g;HRH5r1*Jp~x<V^mjk
zhQ^f}N6=Mb4`nd95Fi)jCGH%|e|~X#J8f<;30z-7_D`D~9Cj`nZ63KBH`C^C5EK3@
z+9oG|e<;wlD1$FSo8$>x3q;#1xrpv7PgCWmtrH%dUG{9JjrR(d`=`x5_{vG!Yl-S+
z+E)IBwz0O#9|hVv#6b^zVj=M1&bv+W1g={|+l*X9OI_;M^3xUoD>wKyPg{s)sOpNs
zbK;YLcU_)F8A{1^anBe#+%sOffwWfP>H8ucu2fg_6T+n-UGgq1eW!|}EGNNTku!o$
ze%XK7wtjb&Fn8BeUxkAO{YEA=^G*75+us`$om(qYvB`8YI~Py1R0xxc&nuqT)Gbn7
zgKE43@kTCOfZ>7Y2+(OxYn#TR80hwApIj0adL@3Kc=0K&le)N0gBRDy=lk%xYkg)U
zOnQ4TX+u?s8*p5Fda$04eO?NzXM3I_>-oIluOzHzKfHuk&weCm`T4LPw&RjFhJ+(;
z&;AJ?61HdmLW`f(`8VgqLBK)4LBK)4LBK)4LBK)4LBK)4LBK)4LEyhbfHDzVGuS>E
zrEi*Y^OT`zrwqlg)oTq~L#8!uXDB)`B|sU9cq|*6O=JtHihL#>JbwIhTVw)^9shs9
zIT40+<yqHHJ2!G*Bip$T9NzjkU3ah}b4F!ub>z#zj?6im%<|YeBfq)4-3`D`Fo<m>
V9%^GdD)%bIeahi~uk+a5?hlnL@nZl0
literal 0
HcmV?d00001
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
GIT binary patch
literal 7276
zcmeHMUr3Wt6hGge%{gakDV7O-{Q<!u45ySOQ@+hVH8*k+RG7>*a<G5olpKN)kz_<r
z5J5o?z4X#s5JD6__#E}nLp{_(y%<z)(L>lh-@Rv>(_t1{5Z~2r_uliJ@0@dgzx$m#
zkJtL$y_NO`yKJ)owH48#ja<?yZrprY(dqAZ56O-v@@_+>Si~)QB)_-IAfZZKBdDxX
zB*@pGx#I}xDIh~F`+*?TXfEI}b%V14oGS`|EU|`qtWnQrd%LTvakPbLVbggBKoJ6r
zghm#P@GHwCiMinBgg*1kn^KFGYj!mm0Rw1IZ-e?-cGF_4MX_^7>JF}dOFTy_R0^+f
zAF3)~Wy@k<=-!?{0Jbp^1}?M9Juv9;1qS@S(SY9L^?>6vu2Nx3Fdp#r40;UxZlOOs
z9+}YCZ%~VeQ@XoX3q@k#_(UWT*TyDOT0D`;-kJ(thy=%?;R71AHfF`T)bx=<0y&Tn
zvea)3iES%0#@Hn?k!x&40EMG}LeMqN{Q#HV5uBuj6q;%9JN=cDhx~f~ac*i!EJ3Fy
zIT+UMIG0p{#*!Xopv%Y}YsVtfPBc+m;ts;p^5Z>gX``J^{pM*iUJAC3qS0m(v8|@9
zy1_sX{)o1@>+kJZ+GecWQ~}xqPhgow+7`qnx+>j^6`!_tc-8)PVJ&SUVwM(o+dOT?
zo@5UfjkXsew$-%F{|#-=9-RqiX=~$&9{L0ZBDvGsCU^qNMWih$HqmUS{H6G`b-?ZY
zbG8g^Zj4Y?C5d9nWC9U-hDI?jc%Pp$+B)Zq9ha}P@C@{Ah9y^=OX^N;a$J}sx(lDF
zJe8#dxbq?+=%B&o``61iB8a)Q9-oyc&x3j{5lem*_N>PD21%v6<+0%TcqBC$3OAQ=
zm9B?ncMa-haIQcVI|Etr$sT~-fxs|eHpjCOec>AL^e3CtBNl2o`(Vw^uG~3_bLZ$h
zcaApScwe6C)BRqQ8_t9dP6cj2>!R#Ia!UKO6-Z9{6eGzg-%x@kIoYADLvpesK+D8K
zc1VwUdNBkX8GW*k_+TS_vd^XLvnn$(2POh00ww|`0ww|`0ww|`0ww|`0ww|`0wx0g
z9RgU1$gf~K#3~(8tj%MEq7^F?L!MsGkjJNcMvV$Z2P+9+g(4J81t-F(Y^fqsi3eLX
z`*WLD0<?Aie?c`7f_45u=MSScQZmm*?Sr}Ze2CVS)-zR9Chre_F0E&(*(B!9zH3Py
z?t1R~?ppkG*#&UPHS0{fT9OM#@xn36^5r;pHT*O#9B-sgcF+g<LZAQ2FYUBa_rCz|
CRyY#?
literal 0
HcmV?d00001
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
GIT binary patch
literal 7276
zcmeHMUr3Wt6hGhhZF8=fs|?Emzy3f`Ayz6esc-Ua{;8?xB&aBvIdV0p<R&G84~3Ez
zL16?1Jwy>fFA-FrL=QqQ5e4;75B1PX7S&q>LEUrjJ=?TKS<Hv{u72Ko&pr2?bAP{Y
z=j!oPy}!1=RpL^ePN4Q(w%AG`DVGap?=J8**ZW&l_Z@XLmuYsnL|bO}E)^71S-KKz
ztkV>jtwD2lHR@@gLJ|9dFf7sy5Hb0nVgXcq*8o(7HPmB`dNEt;^LTKykZEBvMF&6;
z0_^YH|B#1Ykf$ijg)k@enP<_oBDB0{RVfb`K!;){6c4dQJ3ou&5{_)kh5k9Yk4~s%
zKjB7H)xg1y#lSFwHK7nJVImA%Wskq1G0+rhXl`l`nKi)xxcA~IHFgBPp{b@Z!1aAn
zzpX3Ur?c;{9*-nUf34mb?T*CzqCIiFBbm_SJ&E*EfB0xL+!2eE>D0C|D~_qsw*nH#
zfP|2x{`ZhL7c*mwT_O{?#+Cz61o|zKt}fvRg!F>sWI#$`9S#1ZKe>Bbvsu4On3|GD
zFlna%!?6tKijAPL541DTRb-E~ozdeJXd<`7ZHE5oTWfx%jox(X=S~~{P_T34j5eo?
zZ6<A=5{@4H6>TG9uU%={1|7mw7TP3FU_Xwu-Ip)XuyUz8_q0{Q)5@3kf2K`F%+i8r
zbEl2JlWcR&XnP=In@QW~ztA>uV_zgq+g6e2p-y5TlRJIdBu`*JhP3s{muSeVe#|{>
zRd9aOh%-f-A0t$4ib64^GJy;|MWdW&T^IL^{)>CYg~!)6hz#_3kR?}~E4CHFWS2Ba
zbZ5U)MJgMR;Eu|Opo#`_@1ISdi6Z8XV%)1xo(J2(p6=dH(w>9=-k{iMHNQK2C>~8D
zJ0t7zgv$1-c^5f#gE-fqklleS`J^wvU_+=4Fq`Auh`#U)c>9x0>JbZDKKsF%i#?fh
z6lc!SdFC8#z6d_qQ)e~@QO=!79h^$sfX+qPgXEO<=~W;(=~Ij(r+h;Rn&f1MULBH?
z9SK@09<oDv)YAt;!jaM^`-l%V(kJ^u%6=B6M%KVWz(T-6z(T-6z(T-6z(T-6z(T-6
zz(T-6;J-rvD-p9R*lu}BcMNOuSfMD#3Po$6HqaVqG6U_rLQ%y^0$8Ew3@5^Ukwm&w
zk*dUl)5HGU7L@?~m;b+@nh3!<dd2&d*GA^eGp~Ix_klODrSt2VDk{m#ZSUvTGu3Q*
zMh?6}<Zn(~H-^VsjKqgJ-;sOo3}R!<`^b3S@{O;d4#`h^u^XrJ<34<FOmeb=dvGtw
z(I(nPADMQlA0IKMwiFre0|mZ!=Czc*Iex~NJbat4p$^$8OHS_C+p^w=4~|Jrc5n~w
RB{|we+vp?H4z15;{{xySDcb-5
literal 0
HcmV?d00001
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
GIT binary patch
literal 7276
zcmeHMO=uHQ5T4yl(lkwNsx8%8==uYNDpKi>wKX&)Nq;m=uqlE{4XFtxkUwl%Y$+5(
zEVhWCAcCR?4|?$?NI{eydJ{ws9`v9GMW~=2ym*lMX5X76{Rvtn5oD)L-oDv+GxNT0
zGcTQ9>T~s!TWT$`*$gDhLSl7X(<E+Oc~sWn>vIjt)`#*=O{W;eEn2dAw8$W#Np&3<
z*q}(T+JffRPP9`%hAQ?0L8w%_fam0ZmNIBrRsdv)3EDBCk<WHJY&N{w#H>)9d;p*f
z0fs}v3tIZ+C6dHKaCgF<MdqSZq36mS^;*IJD%6;vW}4kJYI{*EJdj}<x4$aZ(GHdT
zJKTe&3Ygfo7#L1Zx8Dz2m<kh@*~`^G;P(3aecpiI+3j(I^#txxVO!7^c)JJOn!Q8V
z4@IJ*D*FwpiE!HK>QO_{csMZ{O(xWlv9y{<rn9%kgJ+|`kyv=2O47`{m=^2bb4VZu
z5<-^7UqfQv%$zZHiA>}gr2wFC^p6U<B0LUo?RCM)q|ibGUHr~?;aI!R*>{AyS`=Gw
z(n=17X)Dep13_b)3^34TWRJ<==tK*;C@gXNVSMSquC=t$lTPEtY12LwY##-q%`9?T
zNt>-!Ll6FlwwbxNmMm>kChjT^ZGtB-P9Sab;vu>y-HaEWwpMu2`g(pXZ6ae<7PxGj
zHtm^Yiwj2EQ<2+B+GhWTw#WBRhO@La^Fj~30s~Rp>C+~70^>QPEhQeJ=??jG;c07w
zYkOzR8QNT!p=^*O$|+L`MCutD#XRqQe$J@toHG{OzS77m(AOzeTyZWLwsV&e;gaai
zf2Q(MHYvcJ6&XPrU2J@OxpX;-n44;FtwePm45yRv)ED8;q<wFY40N|N9z2tXrpH3z
zh7xXKzgu!cqizc43RJK&kQJZo0qE)X4*^zlJR30<?g3AK@<}6NVJKy9Of2k`J4bcy
z9E0c1(dR4A^W(ivp9l3?FkypJfg8}isCtl|%04{`q$hjIk@Qq=s6dmR{LrIAdh#Pc
z%j83T$c{$(U<f!e_T(S&K_Pqc&$aBmB6Fo*=m_Wt=m_Wt=m_Wt=m_Wt=m_Wt=m_Wt
z=m`9G2w*2-bqCujc4>`aZyq}oP1vCrboaOi-Cn0VpmiwPSW5sq6ro@`I2ul8TNRm3
zJeY0l&u!ibP}lwc1@%M-*4aB9-?iRI(IRWT4;J42E>>OK$kb68yEXKwxRI%6lbkvA
z29a+}9kl<5eCuj`_R;>td%&Sq&Dl|p<L!lZ$HIs6U5%;vL#RVRec)cF17Bo|dbB}X
ZTtj*s(I@)G7#Is}QT|n29=*8lKLO@)D>DE9
literal 0
HcmV?d00001
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
GIT binary patch
literal 7156
zcmeHMK}=IY5S{;jTS_TX6oMiNUo;`ngkUH{F@}=%1I5xtN=zgHQfwnlXiKCJX-M#3
zB;w73G10_>@#xJ9C*y&{oAKg74|?#TCdz?}C!G1a(^jA=YB4c?q20f;vokyI{hbZN
z&4EB)wX4ykI-S66Q?`x6LeefLT7svw$UtCNbw5#COqphvQ@mtr?@~b_mHH;IVW%mu
zRfEUwGo+`13U&AcF{sshK<vo}-YW2JY5=N&jr6c_FGk0F9uGb3LoI54(E;!d0mkBE
z&&=&tS11a)5ay&ldKOKqBg?f%n#~OZ=+NMVh9#V|o3Usv;mCGK=--y>I6|#-geOQ<
z0|$nMVE9A5(I^~15p7(>ParsCgrmVocp~cW4H@7br6@HF!CVmT9WqRPpVS|jO3v!|
z8`IN?tUu7F$CIf<dN!F!>yvX?J)OxGrsiW;lCjC@M4Qg71FbmLn%@gZpa>F1miwPW
z;@pqMkX<H|xTaPHPz3ttC0$d(4+!ZU$w^L1p_LzgrN4gRRK!1UPMBJgA^3Tv0K;*R
z#)^%hvF9cbbd}h{Hke%Sl8N#XcM|45K0dmgHh$B&-#u;Sp};uGMw?Ua?R(lhjV5~V
zTeL0Tdh06Cw&)P1O3)^G0{a5dwkltu8_K;@`DyEfmz}Ryx6{UFh5OypW+uLJ()LuM
z`aNwce?r@MTjfT9whoc#p<l4a$(>J|<O%FoiMAQ}5-oMBpUY317mVEC^E_<<+M#Mw
z6y8%l6UeRSX_VuV^WvV-cerO<bOUL}g{AdHOs+ImZ1Or@mzMdImcCO(D$7Z5SL8bA
z;)mV$Rihss5$28tTB~rdV7rt_&3uve98F@9XDd^&%jslxE}m$u5GrRMR@^nITcoiD
zHFyVN@+n+^p<r|bD4Wx}O};1wy8YQE_XMx462Gvy@KZb{aq*a(7mvy2tI&&!{r*UZ
z<fb#JLqmxha9o@{Sk7smUj>%4KJSs`oNqWmvz+b7Q?|{o1T8OTJG`EIJ{S^?ygu6}
ze5kQL?^8(eyC(l+Jy-}>2v`VM2v`VM2v`VM2v`VM2v`VM2v`XGcL-1=VrvE4Eu(Z#
zQ*E9q6zx=@7&iKhVI%A}Cd>*&7nTI5LJ^N;W3!2Dp;VEt#Dmj=e{PFPfWG7ZFSsVc
zu&zAl{%Y1n_RO<c`#|@`yXpG9^~@EOx%(rZ_SQ4kY%<H^Z!mNpz_p$qe@N$_>7nbt
G*S-Ovr2YK>
literal 0
HcmV?d00001
diff --git a/tests/test_projector/ref_proj_x.txt b/tests/test_projector/ref_proj_x.txt
deleted file mode 100644
index 94ce2a3..0000000
--- 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 b8a1448..0000000
--- 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 eec624b..0000000
--- 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
--
GitLab