From 6fc82b2f929af2c3772124d2de166e81cb15bcec Mon Sep 17 00:00:00 2001
From: Jan Caron <j.caron@fz-juelich.de>
Date: Wed, 16 Apr 2014 15:26:12 +0200
Subject: [PATCH] Documentation and cleanup of the pyramid package!
Renames:
datacollection --> dataset
optimizer --> reconstruction
scripts:
interactive_setup is now implemented in an extended Spyder startup file.
some other scripts are now
the rest is NOT adapted, yet, which is the next task at hand
furthermore scripts will be sorted and unused ones deleted in the next commit
---
.hgignore | 3 +-
docs/pyramid.rst | 4 +-
logfile.log | 3 +
pyramid/__init__.py | 10 +-
pyramid/analytic.py | 8 +-
pyramid/costfunction.py | 176 +++++---
pyramid/datacollection.py | 102 -----
pyramid/dataset.py | 159 +++++++
pyramid/forwardmodel.py | 125 ++++--
pyramid/kernel.py | 163 +++-----
pyramid/logfile.log | 0
pyramid/logging.ini | 2 +-
pyramid/magcreator.py | 146 +++++--
pyramid/magdata.py | 147 ++++---
pyramid/numcore/__init__.py | 3 +-
pyramid/numcore/kernel_core.pyx | 153 +------
pyramid/optimizer.py | 143 -------
pyramid/phasemap.py | 132 +++---
pyramid/phasemapper.py | 387 +++++++++++-------
pyramid/projector.py | 363 ++++++++++------
pyramid/reconstruction.py | 117 ++++++
scripts/abstract_pictures.py | 115 ++++++
scripts/compare methods/logfile.log | 0
scripts/create distributions/logfile.log | 0
scripts/edinburgh_test.py | 13 +-
scripts/interactive_setup.py | 23 +-
scripts/logfile.log | 0
scripts/paper 1/logfile.log | 84 ----
scripts/simple_reconstruction.py | 175 ++++++--
...test_datacollection.py => test_dataset.py} | 0
...est_optimize.py => test_reconstruction.py} | 6 +-
31 files changed, 1667 insertions(+), 1095 deletions(-)
delete mode 100644 pyramid/datacollection.py
create mode 100644 pyramid/dataset.py
delete mode 100644 pyramid/logfile.log
delete mode 100644 pyramid/optimizer.py
create mode 100644 pyramid/reconstruction.py
create mode 100644 scripts/abstract_pictures.py
delete mode 100644 scripts/compare methods/logfile.log
delete mode 100644 scripts/create distributions/logfile.log
delete mode 100644 scripts/logfile.log
delete mode 100644 scripts/paper 1/logfile.log
rename tests/{test_datacollection.py => test_dataset.py} (100%)
rename tests/{test_optimize.py => test_reconstruction.py} (79%)
diff --git a/.hgignore b/.hgignore
index 5f75704..3e01d21 100644
--- a/.hgignore
+++ b/.hgignore
@@ -10,4 +10,5 @@ syntax: glob
*.c
output*
build*
-Pyramid.egg-info*
\ No newline at end of file
+Pyramid.egg-info*
+scripts/dyn_0090mT_dyn.h5_dat.h5
diff --git a/docs/pyramid.rst b/docs/pyramid.rst
index 8bd3b4c..9bd1af8 100644
--- a/docs/pyramid.rst
+++ b/docs/pyramid.rst
@@ -17,10 +17,10 @@ pyramid Package
:show-inheritance:
:special-members:
-:mod:`datacollection` Module
+:mod:`dataset` Module
----------------------------
-.. automodule:: pyramid.datacollection
+.. automodule:: pyramid.dataset
:members:
:show-inheritance:
:special-members:
diff --git a/logfile.log b/logfile.log
index e69de29..bde3b6a 100644
--- a/logfile.log
+++ b/logfile.log
@@ -0,0 +1,3 @@
+2014-04-10 16:49:29: DEBUG @ <apptools.preferences.preferences>: loading preferences from <C:\Users\Jan\AppData\Roaming\Enthought\mayavi_e3\preferences.ini>
+2014-04-10 16:49:29: DEBUG @ <apptools.preferences.preferences>: loading preferences from <<open file 'C:\\Python27\\lib\\site-packages\\mayavi\\preferences\\preferences.ini', mode 'rb' at 0x04D00F98>>
+2014-04-10 16:49:29: DEBUG @ <apptools.preferences.preferences>: loading preferences from <<open file 'C:\\Python27\\lib\\site-packages\\tvtk\\plugins\\scene\\preferences.ini', mode 'rb' at 0x050F85A0>>
diff --git a/pyramid/__init__.py b/pyramid/__init__.py
index 2316e34..5ffe7ac 100644
--- a/pyramid/__init__.py
+++ b/pyramid/__init__.py
@@ -17,15 +17,13 @@ phasemap
Class for the storage of phase data.
analytic
Create phase maps for magnetic distributions with analytic solutions.
-datacollection
+dataset
Class for collecting pairs of phase maps and corresponding projectors.
forwardmodel
Class which represents a phase mapping strategy.
costfunction
Class for the evaluation of the cost of a function.
-optimizer
- Provides strategies for optimizing first guess magnetic distributions.
-reconstructor
+reconstruction
Reconstruct magnetic distributions from given phasemaps.
Subpackages
@@ -35,6 +33,8 @@ numcore
"""
+# TODO: logging setup at script level, only logging itself should be done in the package
+# maybe include in startup script!
import logging, logging.config
import os
@@ -43,4 +43,4 @@ import os
LOGGING_CONF = os.path.join(os.path.dirname(os.path.realpath(__file__)), 'logging.ini')
-logging.config.fileConfig(LOGGING_CONF)
+logging.config.fileConfig(LOGGING_CONF, disable_existing_loggers=False)
diff --git a/pyramid/analytic.py b/pyramid/analytic.py
index bf2ffcc..5216fcd 100644
--- a/pyramid/analytic.py
+++ b/pyramid/analytic.py
@@ -13,7 +13,10 @@ from numpy import pi
from pyramid.phasemap import PhaseMap
+import logging
+
+LOG = logging.getLogger(__name__)
PHI_0 = -2067.83 # magnetic flux in T*nm²
@@ -42,7 +45,7 @@ def phase_mag_slab(dim, a, phi, center, width, b_0=1):
The phase as a 2-dimensional array.
'''
-
+ LOG.debug('Calling phase_mag_slab')
# Function for the phase:
def phi_mag(x, y):
def F_0(x, y):
@@ -98,6 +101,7 @@ def phase_mag_disc(dim, a, phi, center, radius, height, b_0=1):
The phase as a 2-dimensional array.
'''
+ LOG.debug('Calling phase_mag_disc')
# Function for the phase:
def phi_mag(x, y):
r = np.hypot(x - x0, y - y0)
@@ -145,6 +149,7 @@ def phase_mag_sphere(dim, a, phi, center, radius, b_0=1):
The phase as a 2-dimensional array.
'''
+ LOG.debug('Calling phase_mag_sphere')
# Function for the phase:
def phi_mag(x, y):
r = np.hypot(x - x0, y - y0)
@@ -191,6 +196,7 @@ def phase_mag_vortex(dim, a, center, radius, height, b_0=1):
The phase as a 2-dimensional array.
'''
+ LOG.debug('Calling phase_mag_vortex')
# Function for the phase:
def phi_mag(x, y):
r = np.hypot(x - x0, y - y0)
diff --git a/pyramid/costfunction.py b/pyramid/costfunction.py
index 811a2a9..71a60d4 100644
--- a/pyramid/costfunction.py
+++ b/pyramid/costfunction.py
@@ -1,72 +1,156 @@
# -*- coding: utf-8 -*-
-"""
-Created on Fri Dec 13 10:29:11 2013
+"""This module provides the :class:`~.Costfunction` class which represents a strategy to calculate
+the so called `cost` of a threedimensional magnetization distribution."""
-@author: Jan
-"""# TODO: Docstring!
+# TODO: better names for variables (no uppercase, more than one letter)
+
+# TODO: Regularisation Class?
import numpy as np
+from scipy.sparse.linalg import LinearOperator
+from scipy.sparse import eye
+
+import logging
class Costfunction:
- # TODO: Docstring!
-
- def __init__(self, y, F):
- '''TEST DOCSTRING FOR INIT'''
- # TODO: Docstring!
- self.y = y # TODO: get y from phasemaps!
- self.F = F # Forward Model
- self.Se_inv = np.eye(len(y))
+
+ '''Class for calculating the cost of a 3D magnetic distributions in relation to 2D phase maps.
+
+ Represents a strategy for the calculation of the `cost` of a 3D magnetic distribution in
+ relation to two-dimensional phase maps. The `cost` is a measure for the difference of the
+ simulated phase maps from the magnetic distributions to the given set of phase maps.
+ Furthermore this class provides convenient methods for the calculation of the derivative
+ :func:`~.jac` or the product with the Hessian matrix :func:`~.hess_dot` of the costfunction,
+ which can be used by optimizers.
+
+ Attributes
+ ----------
+ y : :class:`~numpy.ndarray` (N=1)
+ Vector which lists all pixel values of all phase maps one after another. Usually gotten
+ via the :class:`~.DataSet` classes `phase_vec` property.
+ fwd_model : :class:`~.ForwardModel`
+ The Forward model instance which should be used for the simulation of the phase maps which
+ will be compared to `y`.
+ lam : float, optional
+ Regularization parameter used in the Hessian matrix multiplication. Default is 0.
+
+ '''
+
+ LOG = logging.getLogger(__name__+'.Costfunction')
+
+ def __init__(self, y, fwd_model, lam=0):
+ self.LOG.debug('Calling __init__')
+ self.y = y
+ self.fwd_model = fwd_model
+ self.lam = lam
+ self.Se_inv = eye(len(y))
+ self.LOG.debug('Created '+str(self))
def __call__(self, x):
- '''TEST DOCSTRING FOR CALL'''
- # TODO: Docstring!
+ self.LOG.debug('Calling __call__')
y = self.y
F = self.F
Se_inv = self.Se_inv
-# print 'Costfunction - __call__ - input: ', len(x)
- result = (F(x)-y).dot(Se_inv.dot(F(x)-y))
-# print 'Costfunction - __call__ - output:', result
- return result
+ return (F(x)-y).dot(Se_inv.dot(F(x)-y))
+
+ def __repr__(self):
+ self.LOG.debug('Calling __repr__')
+ return '%s(fwd_model=%r, lam=%r)' % (self.__class__, self.fwd_model, self.lam)
+
+ def __str__(self):
+ self.LOG.debug('Calling __str__')
+ return 'Costfunction(fwd_model=%s, lam=%s)' % (self.__class__, self.fwd_model, self.lam)
def jac(self, x):
- # TODO: Docstring!
+ '''Calculate the derivative of the costfunction for a given magnetization distribution.
+
+ Parameters
+ ----------
+ x : :class:`~numpy.ndarray` (N=1)
+ Vectorized magnetization distribution, for which the Jacobi vector is calculated.
+
+ Returns
+ -------
+ result : :class:`~numpy.ndarray` (N=1)
+ Jacobi vector which represents the cost derivative of all voxels of the magnetization.
+
+ '''
+ self.LOG.debug('Calling jac')
y = self.y
- F = self.F
+ F = self.fwd_model
Se_inv = self.Se_inv
-# print 'Costfunction - jac - input: ', len(x)
- result = F.jac_T_dot(x, Se_inv.dot(F(x)-y))
-# print 'Costfunction - jac - output:', len(result)
- return result
+ return F.jac_T_dot(x, Se_inv.dot(F(x)-y))
def hess_dot(self, x, vector):
- # TODO: Docstring!
- F = self.F
+ '''Calculate the product of a `vector` with the Hession matrix of the costfunction.
+
+ Parameters
+ ----------
+ x : :class:`~numpy.ndarray` (N=1)
+ Vectorized magnetization distribution at which the Hessian is calculated. The Hessian
+ is constant in this case, thus `x` can be set to None (it is not used int the
+ computation). It is implemented for the case that in the future nonlinear problems
+ have to be solved.
+ vector : :class:`~numpy.ndarray` (N=1)
+ Vectorized magnetization distribution which is multiplied by the Hessian.
+
+ Returns
+ -------
+ result : :class:`~numpy.ndarray` (N=1)
+ Product of the input `vector` with the Hessian matrix of the costfunction.
+
+ '''
+ self.LOG.debug('Calling hess_dot')
+ F = self.fwd_model
Se_inv = self.Se_inv
-# print 'Costfunction - hess_dot - input: ', len(vector)
- result = F.jac_T_dot(x, Se_inv.dot(F.jac_dot(x, vector)))
-# print 'Costfunction - hess_dot - output:', len(result)
- return result
- # TODO: Murks!
+ lam = self.lam
+ return F.jac_T_dot(x, Se_inv.dot(F.jac_dot(x, vector))) + lam*vector
-class CFAdapter:
- # TODO: Useless at the moment, because the interface is exactly the same. Use as template!
-
- def __init__(self, costfunction):
- # TODO: Docstring!
- self.costfunction = costfunction
+class CFAdapterScipyCG(LinearOperator):
- def __call__(self, x):
- # TODO: Docstring!
- return self.costfunction(x)
+ '''Adapter class making the :class:`~.Costfunction` class accessible for scipy cg methods.
- def jac(self, x):
- # TODO: Docstring!
- return self.costfunction.jac(x)
+ 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.optimice_sparse_cg` function of the
+ :mod:`~pyramid.reconstruction` module.
- def hess_dot(self, x, vector):
- # TODO: Docstring!
- return self.costfunction.hess_dot(x, vector)
+ Attributes
+ ----------
+ cost : :class:`~.Costfunction`
+ :class:`~.Costfunction` class which should be made usable in the
+ :func:`~.scipy.sparse.linalg.cg` function.
+
+ '''
+
+ 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.
+
+ Parameters
+ ----------
+ vector : :class:`~numpy.ndarray` (N=1)
+ Vector which will be multiplied by the Hessian matrix provided by the adapted
+ costfunction.
+
+ '''
+ self.LOG.debug('Calling matvec')
+ return self.cost.hess_dot(None, vector)
+
+ @property
+ def shape(self):
+ return (3*self.cost.fwd_model.size_3d, 3*self.cost.fwd_model.size_3d)
+
+ @property
+ def dtype(self):
+ return np.dtype("d")
diff --git a/pyramid/datacollection.py b/pyramid/datacollection.py
deleted file mode 100644
index 496173e..0000000
--- a/pyramid/datacollection.py
+++ /dev/null
@@ -1,102 +0,0 @@
-# -*- coding: utf-8 -*-
-"""This module provides the :class:`~.DataCollection` class for the collection of phase maps
-and additional data like corresponding projectors."""
-
-
-import logging
-
-import numpy as np
-from numbers import Number
-
-from pyramid.phasemap import PhaseMap
-from pyramid.projector import Projector
-
-
-class DataCollection(object):
-
- '''Class for collecting phase maps and corresponding projectors.
-
- Represents a collection of (e.g. experimentally derived) phase maps, stored as
- :class:`~.PhaseMap` objects and corresponding projectors stored as :class:`~.Projector`
- objects. At creation, the grid spacing `a` and the dimension `dim_uv` of the projected grid.
- Data can be added via the :func:`~.append` method, where a :class:`~.PhaseMap` and a
- :class:`~.Projector` have to be given as tuple argument.
-
- Attributes
- ----------
- a: float
- The grid spacing in nm.
- dim_uv: tuple (N=2)
- Dimensions of the projected grid.
- phase_maps:
- A list of all stored :class:`~.PhaseMap` objects.
- projectors:
- A list of all stored :class:`~.Projector` objects.
- phase_vec: :class:`~numpy.ndarray` (N=1)
- The concatenaded, vectorized phase of all ;class:`~.PhaseMap` objects.
-
- '''
-
- @property
- def a(self):
- return self._a
-
- @property
- def dim_uv(self):
- return self._dim_uv
-
- @property
- def phase_vec(self):
- return np.concatenate([p.phase_vec for p in self.phase_maps])
-
- @property
- def phase_maps(self):
- return [d[0] for d in self.data]
-
- @property
- def projectors(self):
- return [d[1] for d in self.data]
-
- def __init__(self, a, dim_uv, b_0):
- self.log = logging.getLogger(__name__)
- assert isinstance(a, Number), 'Grid spacing has to be a number!'
- assert a >= 0, 'Grid spacing has to be a positive number!'
- self._a = a
- assert isinstance(dim_uv, tuple) and len(dim_uv)==2, \
- 'Dimension has to be a tuple of length 2!'
- assert np.all([])
- self._dim_uv = dim_uv
- self.b_0 = b_0
- self.data = []
- # TODO: make it work
-# self.log.info('Created:', str(self))
-
- def __repr__(self):
- self.log.info('Calling __repr__')
- return '%s(a=%r, dim_uv=%r)' % (self.__class__, self.a, self.dim_uv)
-
- def __str__(self):
- self.log.info('Calling __str__')
- return 'DataCollection(%s, a=%s, dim_uv=%s, data_count=%s)' % \
- (self.__class__, self.a, self.dim_uv, len(self.data))
-
- def append(self, (phase_map, projector)):
- '''Appends a data pair of phase map and projection infos to the data collection.`
-
- Parameters
- ----------
- (phase_map, projector): tuple (N=2)
- tuple which contains a :class:`~.PhaseMap` object and a :class:`~.Projector` object,
- which should be added to the data collection.
-
- Returns
- -------
- None
-
- '''
- self.log.info('Calling append')
- assert isinstance(phase_map, PhaseMap) and isinstance(projector, Projector), \
- 'Argument has to be a tuple of a PhaseMap and a Projector object!'
- assert phase_map.dim_uv == self.dim_uv, 'Added phasemap must have the same dimension!'
- assert projector.dim_uv == self.dim_uv, 'Projector dimensions must match!'
- self.data.append((phase_map, projector))
diff --git a/pyramid/dataset.py b/pyramid/dataset.py
new file mode 100644
index 0000000..4cecfe8
--- /dev/null
+++ b/pyramid/dataset.py
@@ -0,0 +1,159 @@
+# -*- coding: utf-8 -*-
+"""This module provides the :class:`~.DataSet` class for the collection of phase maps
+and additional data like corresponding projectors."""
+
+
+import numpy as np
+from numbers import Number
+
+import matplotlib.pyplot as plt
+
+from pyramid.phasemap import PhaseMap
+from pyramid.phasemapper import PMConvolve
+from pyramid.projector import Projector
+
+import logging
+
+
+class DataSet(object):
+
+ '''Class for collecting phase maps and corresponding projectors.
+
+ Represents a collection of (e.g. experimentally derived) phase maps, stored as
+ :class:`~.PhaseMap` objects and corresponding projectors stored as :class:`~.Projector`
+ objects. At creation, the grid spacing `a` and the dimension `dim_uv` of the projected grid.
+ Data can be added via the :func:`~.append` method, where a :class:`~.PhaseMap` and a
+ :class:`~.Projector` have to be given as tuple argument.
+
+ Attributes
+ ----------
+ a: float
+ The grid spacing in nm.
+ dim_uv: tuple (N=2)
+ Dimensions of the projected grid.
+ phase_maps:
+ A list of all stored :class:`~.PhaseMap` objects.
+ projectors:
+ A list of all stored :class:`~.Projector` objects.
+ phase_vec: :class:`~numpy.ndarray` (N=1)
+ The concatenaded, vectorized phase of all ;class:`~.PhaseMap` objects.
+
+ '''
+
+ LOG = logging.getLogger(__name__+'.DataSet')
+
+ @property
+ def a(self):
+ return self._a
+
+ @property
+ def dim_uv(self):
+ return self._dim_uv
+
+ @property
+ def phase_vec(self):
+ return np.concatenate([p.phase_vec for p in self.phase_maps])
+
+ @property
+ def phase_maps(self):
+ return [d[0] for d in self.data]
+
+ @property
+ def projectors(self):
+ return [d[1] for d in self.data]
+
+ def __init__(self, a, dim_uv, b_0):
+ 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!'
+ self._a = a
+ assert isinstance(dim_uv, tuple) and len(dim_uv)==2, \
+ 'Dimension has to be a tuple of length 2!'
+ self._dim_uv = dim_uv
+ self.b_0 = b_0
+ self.data = []
+ self.LOG.debug('Created:', str(self))
+
+ def __repr__(self):
+ self.LOG.debug('Calling __repr__')
+ return '%s(a=%r, dim_uv=%r, b_0=%r)' % (self.__class__, self.a, self.dim_uv, self.b_0)
+
+ def __str__(self):
+ self.LOG.debug('Calling __str__')
+ return 'DataSet(a=%s, dim_uv=%s, b_0=%s)' % (self.a, self.dim_uv, self.b_0)
+
+ def append(self, (phase_map, projector)):
+ '''Appends a data pair of phase map and projection infos to the data collection.`
+
+ Parameters
+ ----------
+ (phase_map, projector): tuple (N=2)
+ tuple which contains a :class:`~.PhaseMap` object and a :class:`~.Projector` object,
+ which should be added to the data collection.
+
+ Returns
+ -------
+ None
+
+ '''
+ self.LOG.debug('Calling append')
+ assert isinstance(phase_map, PhaseMap) and isinstance(projector, Projector), \
+ 'Argument has to be a tuple of a PhaseMap and a Projector object!'
+ assert phase_map.dim_uv == self.dim_uv, 'Added phasemap must have the same dimension!'
+ assert projector.dim_uv == self.dim_uv, 'Projector dimensions must match!'
+ self.data.append((phase_map, projector))
+
+ def display(self, phase_maps=None, title='Phase Map', cmap='RdBu', limit=None, norm=None):
+ '''Display all phasemaps saved in the :class:`~.DataSet` as a colormesh.
+
+ Parameters
+ ----------
+ phase_maps : list of :class:`~.PhaseMap`, optional
+ List of phase_maps to display with annotations from the projectors. If none are
+ given, the phase_maps in the dataset are used (which is the default behaviour).
+ title : string, optional
+ The main part of the title of the plots. The default is 'Phase Map'. Additional
+ projector info is appended to this.
+ cmap : string, optional
+ The :class:`~matplotlib.colors.Colormap` which is used for the plots as a string.
+ The default is 'RdBu'.
+ limit : float, optional
+ Plotlimit for the phase in both negative and positive direction (symmetric around 0).
+ If not specified, the maximum amplitude of the phase is used.
+ norm : :class:`~matplotlib.colors.Normalize` or subclass, optional
+ Norm, which is used to determine the colors to encode the phase information.
+ If not specified, :class:`~matplotlib.colors.Normalize` is automatically used.
+
+ Returns
+ -------
+ None
+
+ '''
+ self.LOG.debug('Calling display')
+ if phase_maps is None:
+ phase_maps = self.phase_maps
+ [phase_map.display_phase('{} ({})'.format(title, self.projectors[i].get_info()),
+ cmap, limit, norm)
+ for (i, phase_map) in enumerate(self.phase_maps)]
+ plt.show()
+
+ def create_phase_maps(self, mag_data):
+ '''Create a list of phasemaps with the projectors in the dataset for a given
+ :class:`~.MagData` object.
+
+ Parameters
+ ----------
+ mag_data : :class:`~.MagData`
+ Magnetic distribution to which the projectors of the dataset should be applied.
+
+ Returns
+ -------
+ phase_maps : list of :class:`~.phasemap.PhaseMap`
+ A list of the phase_maps resulting from the projections specified in the dataset.
+
+ Notes
+ -----
+ For the phasemapping, the :class:`~.PMConvolve` class is used.
+
+ '''
+ return [PMConvolve(self.a, proj, self.b_0)(mag_data) for proj in self.projectors]
diff --git a/pyramid/forwardmodel.py b/pyramid/forwardmodel.py
index 44694cc..0dae059 100644
--- a/pyramid/forwardmodel.py
+++ b/pyramid/forwardmodel.py
@@ -1,9 +1,6 @@
# -*- coding: utf-8 -*-
-"""
-Created on Mon Jan 06 14:02:18 2014
-
-@author: Jan
-"""
+"""This module provides the :class:`~.ForwardModel` class which represents a strategy to map a
+threedimensional magnetization distribution onto a two-dimensional phase map."""
import numpy as np
@@ -11,8 +8,40 @@ import numpy as np
from pyramid.kernel import Kernel
from pyramid.projector import Projector
+import logging
+
+
class ForwardModel:
- # TODO: Docstring!
+
+ '''Class for mapping 3D magnetic distributions to 2D phase maps.
+
+ Represents a strategy for the mapping of a 3D magnetic distribution to two-dimensional
+ phase maps. Can handle a list of `projectors` of :class:`~.Projector` objects, which describe
+ different projection angles, so many phase_maps can be created from one magnetic distribution.
+
+ Attributes
+ ----------
+ projectors : list of :class:`~.Projector`
+ A list of all :class:`~.Projector` objects representing the projection directions.
+ kernel : :class:`~.Kernel`
+ A kernel which describes the phasemapping of the 2D projected magnetization distribution.
+ a : float
+ The grid spacing in nm. Will be extracted from the `kernel`.
+ dim : tuple (N=3)
+ Dimensions of the 3D magnetic distribution. Is extracted from the first projector of the
+ `projectors` list.
+ dim_uv: tuple (N=2)
+ Dimensions of the projected grid. Is extracted from the `kernel`.
+ size_3d : int
+ Number of voxels of the 3-dimensional grid. Is extracted from the first projector of the
+ `projectors` list. Is extracted from the first projector of the `projectors` list.
+ size_2d : int
+ Number of pixels of the 2-dimensional projected grid. Is extracted from the first projector
+ of the `projectors` list. Is extracted from the first projector of the `projectors` list.
+
+ '''
+
+ LOG = logging.getLogger(__name__+'.ForwardModel')
@property
def projectors(self):
@@ -34,40 +63,80 @@ class ForwardModel:
self._kernel = kernel
def __init__(self, projectors, kernel):
- # TODO: Docstring!
+ self.LOG.debug('Calling __init__')
self.kernel = kernel
self.a = kernel.a
- self.dim_uv = kernel.dim_uv
self.projectors = projectors
+ self.dim = self.projectors[0].dim
+ self.dim_uv = kernel.dim_uv
+ self.size_2d = self.projectors[0].size_2d
+ self.size_3d = self.projectors[0].size_3d
+ self.LOG.debug('Creating '+str(self))
def __call__(self, x):
- # TODO: Docstring!
-# print 'FWD Model - __call__ - input: ', len(x)
+ self.LOG.debug('Calling __call__')
+ result = [self.kernel(projector(x)) for projector in self.projectors]
+ return np.reshape(result, -1)
+
+ def jac_dot(self, x, vector):
+ '''Calculate the product of the Jacobi matrix with a given `vector`.
+
+ Parameters
+ ----------
+ x : :class:`~numpy.ndarray` (N=1)
+ Evaluation point of the jacobi-matrix. The Jacobi matrix is constant for a linear
+ problem, 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 form of the 3D magnetization distribution. First the `x`, then the `y` and
+ lastly the `z` components are listed.
+
+ Returns
+ -------
+ result_vector : :class:`~numpy.ndarray` (N=1)
+ Product of the Jacobi matrix (which is not explicitely calculated) with the input
+ `vector`.
+
+ '''
+ self.LOG.debug('Calling jac_dot')
result = [self.kernel.jac_dot(projector.jac_dot(x)) for projector in self.projectors]
result = np.reshape(result, -1)
-# print 'FWD Model - __call__ - output:', len(result)
return result
-
- # TODO: jac_dot ausschreiben!!!
-
- def jac_dot(self, x, vector):
- # TODO: Docstring! multiplication with the jacobi-matrix (may depend on x)
-# print 'FWD Model - jac_dot - input: ', len(vector)
result = self(vector)
-# print 'FWD Model - jac_dot - output:', len(result)
- return result # The jacobi-matrix does not depend on x in a linear problem
-
+ return result
+
def jac_T_dot(self, x, vector):
- # TODO: Docstring! multiplication with the jacobi-matrix (may depend on x)
- # The jacobi-matrix does not depend on x in a linear problem
-# print 'FWD Model - jac_T_dot - input: ', len(vector)
+ ''''Calculate the product of the transposed Jacobi matrix with a given `vector`.
+
+ Parameters
+ ----------
+ x : :class:`~numpy.ndarray` (N=1)
+ Evaluation point of the jacobi-matrix. The jacobi matrix is constant for a linear
+ problem, thus `x` can be set to None (it is not used int the computation). Is used
+ for the case that in the future nonlinear problems have to be solved.
+ vector : :class:`~numpy.ndarray` (N=1)
+ Vectorized form of all 2D phase maps one after another in one vector.
+
+ Returns
+ -------
+ result_vector : :class:`~numpy.ndarray` (N=1)
+ Product of the transposed Jacobi matrix (which is not explicitely calculated) with
+ the input `vector`.
+
+ '''
+ self.LOG.debug('Calling jac_T_dot')
size_3d = self.projectors[0].size_3d
size_2d = np.prod(self.dim_uv)
result = np.zeros(3*size_3d)
for (i, projector) in enumerate(self.projectors):
result += projector.jac_T_dot(self.kernel.jac_T_dot(vector[i*size_2d:(i+1)*size_2d]))
-# result = [projector.jac_T_dot(self.kernel.jac_T_dot(vector[i*size_2d:(i+1)*size_2d]))
-# for (i, projector) in enumerate(self.projectors)]
- result = np.reshape(result, -1)
-# print 'FWD Model - jac_T_dot - output:', len(result)
- return result
+ return np.reshape(result, -1)
+
+ def __repr__(self):
+ self.LOG.debug('Calling __repr__')
+ return '%s(projectors=%r, kernel=%r)' % (self.__class__, self.projectors, self.kernel)
+
+ def __str__(self):
+ self.LOG.debug('Calling __str__')
+ return 'ForwardModel(%s -> %s, %s projections, kernel=%s)' % \
+ (self.dim, self.dim_uv, len(self.projectors), self.kernel)
diff --git a/pyramid/kernel.py b/pyramid/kernel.py
index 00752a9..8c0859c 100644
--- a/pyramid/kernel.py
+++ b/pyramid/kernel.py
@@ -3,58 +3,40 @@
single magnetized pixel."""
-import logging
-
import numpy as np
import pyramid.numcore.kernel_core as core
+import logging
+
PHI_0 = -2067.83 # magnetic flux in T*nm²
-# TODO: sign?
class Kernel(object):
'''Class for calculating kernel matrices for the phase calculation.
-
-
-
- This module provides the :class:`~.Kernel` class whose instances can be used to calculate and
- store the kernel matrix representing the phase of a single pixel for the convolution used in the
- phase calculation. The phasemap of a single pixel for two orthogonal directions (`u` and `v`) are
- stored seperately as 2-dimensional matrices. The Jacobi matrix of the phasemapping just depends
- on the kernel and can be calculated via the :func:`~.get_jacobi` function. Storing the Jacobi
- matrix uses much memory, thus it is also possible to directly get the multiplication of a given
- vector with the (transposed) Jacobi matrix without explicit calculation of the latter.
- It is possible to load data from and save them to NetCDF4 files. See :class:`~.Kernel` for further
- information.
-
Represents the phase of a single magnetized pixel for two orthogonal directions (`u` and `v`),
which can be accessed via the corresponding attributes. The default elementary geometry is
`disc`, but can also be specified as the phase of a `slab` representation of a single
magnetized pixel. During the construction, a few attributes are calculated that are used in
- the convolution during phase calculation.
-
-
- An instance `kernel` of the :class:`~.Kernel` class is callable via:
-
- .. :function:: kernel(vector)
-
- do stuff
-
- :param str sender: do other stuff
- :return: nix
-
- with `vector` being a :class:`~numpy.ndarray` (N=1).
+ the convolution during phase calculation in the different :class:`~Phasemapper` classes.
+ An instance of the :class:`~.Kernel` class can be called as a function with a `vector`,
+ which represents the projected magnetization onto a 2-dimensional grid.
Attributes
----------
- dim_uv : tuple (N=2)
- Dimensions of the projected magnetization grid.
a : float
The grid spacing in nm.
+ dim_uv : tuple of int (N=2)
+ Dimensions of the 2-dimensional projected magnetization grid from which the phase should
+ be calculated.
+ b_0 : float, optional
+ Saturation magnetization in Tesla, which is used for the phase calculation. Default is 1.
+ numcore : boolean, optional
+ Boolean choosing if Cython enhanced routines from the :module:`~.pyramid.numcore` module
+ should be used. Default is True.
geometry : {'disc', 'slab'}, optional
The elementary geometry of the single magnetized pixel.
u : :class:`~numpy.ndarray` (N=3)
@@ -65,7 +47,7 @@ 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 (N=2)
+ 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``)
@@ -74,23 +56,12 @@ class Kernel(object):
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.
- '''# TODO: Can be used for several PhaseMappers via the fft arguments or via calling!
-
- def __init__(self, a, dim_uv, b_0=1., numcore=True, geometry='disc'):
- '''Constructor for a :class:`~.Kernel` object for representing a kernel matrix.
+ '''
- Parameters
- ----------
- a : float
- The grid spacing in nm.
- dim_uv : tuple (N=2)
- Dimensions of the projected magnetization grid.
- geometry : {'disc', 'slab'}, optional
- The elementary geometry of the single magnetized pixel.
-
- ''' # TODO: Docstring
- self.log = logging.getLogger(__name__)
- self.log.info('Calling __init__')
+ LOG = logging.getLogger(__name__+'.Kernel')
+
+ def __init__(self, a, dim_uv, b_0=1., numcore=True, geometry='disc'):
+ self.LOG.debug('Calling __init__')
# Function for the phase of an elementary geometry:
def get_elementary_phase(geometry, n, m, a):
if geometry == 'disc':
@@ -104,7 +75,7 @@ class Kernel(object):
return 0.5 * (F_a(n-0.5, m-0.5) - F_a(n+0.5, m-0.5)
- F_a(n-0.5, m+0.5) + F_a(n+0.5, m+0.5))
# Set basic properties:
- self.dim_uv = dim_uv # !!! size of the FOV, not the kernel (kernel is bigger)!
+ self.dim_uv = dim_uv # Size of the FOV, not the kernel (kernel is bigger)!
self.a = a
self.numcore = numcore
self.geometry = geometry
@@ -122,39 +93,30 @@ class Kernel(object):
self.slice_fft = (slice(dim_uv[0]-1, 2*dim_uv[0]-1), slice(dim_uv[1]-1, 2*dim_uv[1]-1))
self.u_fft = np.fft.rfftn(self.u, self.dim_fft)
self.v_fft = np.fft.rfftn(self.v, self.dim_fft)
- self.log.info('Created '+str(self))
-
- def __call__(self, x):
- '''Test'''
- self.log.info('Calling __call__')
-# print 'Kernel - __call__:', len(x)
- if self.numcore:
- return self._multiply_jacobi_core(x)
+ if numcore:
+ self.LOG.debug('numcore activated')
+ self._multiply_jacobi_method = self._multiply_jacobi_core
else:
- return self._multiply_jacobi(x)
- # TODO: Bei __init__ variable auf die entsprechende Funktion setzen.
+ self._multiply_jacobi_method = self._multiply_jacobi
+ self.LOG.debug('Created '+str(self))
+
def __repr__(self):
- self.log.info('Calling __repr__')
+ self.LOG.debug('Calling __repr__')
return '%s(a=%r, dim_uv=%r, numcore=%r, geometry=%r)' % \
(self.__class__, self.a, self.dim_uv, self.numcore, self.geometry)
- def jac_dot(self, vector):
- '''TEST'''# TODO: Docstring
- self.log.info('Calling jac_dot')
-# print 'Kernel - jac_dot:', len(vector)
- if self.numcore:
- return self._multiply_jacobi_core(vector)
- else:
- return self._multiply_jacobi(vector)
+ def __str__(self):
+ self.LOG.debug('Calling __str__')
+ return 'Kernel(a=%s, dim_uv=%s, numcore=%s, geometry=%s)' % \
+ (self.a, self.dim_uv, self.numcore, self.geometry)
- def jac_T_dot(self, vector):
- # TODO: Docstring
- self.log.info('Calling jac_dot_T')
-# print 'Kernel - jac_T_dot:', len(vector)
- return self._multiply_jacobi_T(vector)
+ def __call__(self, x):
+ '''Test'''
+ self.LOG.debug('Calling __call__')
+ return self._multiply_jacobi_method(x)
- def _multiply_jacobi(self, vector):
+ def jac_dot(self, vector):
'''Calculate the product of the Jacobi matrix with a given `vector`.
Parameters
@@ -163,14 +125,38 @@ class Kernel(object):
Vectorized form of the magnetization in `u`- and `v`-direction of every pixel
(row-wise). The first ``N**2`` elements have to correspond to the `u`-, the next
``N**2`` elements to the `v`-component of the magnetization.
-
+
Returns
-------
result : :class:`~numpy.ndarray` (N=1)
Product of the Jacobi matrix (which is not explicitely calculated) with the vector.
- '''# TODO: move!
- self.log.info('Calling _multiply_jacobi')
+ '''
+ self.LOG.debug('Calling jac_dot')
+ return self._multiply_jacobi_method(vector)
+
+ def jac_T_dot(self, vector):
+ '''Calculate the product of the transposed Jacobi matrix with a given `vector`.
+
+ Parameters
+ ----------
+ vector : :class:`~numpy.ndarray` (N=1)
+ Vectorized form of the magnetization in `u`- and `v`-direction of every pixel
+ (row-wise). The first ``N**2`` elements have to correspond to the `u`-, the next
+ ``N**2`` elements to the `v`-component of the magnetization.
+
+ Returns
+ -------
+ result : :class:`~numpy.ndarray` (N=1)
+ Product of the transposed Jacobi matrix (which is not explicitely calculated) with
+ the vector.
+
+ '''
+ self.LOG.debug('Calling jac_dot_T')
+ return self._multiply_jacobi_T(vector)
+
+ def _multiply_jacobi(self, vector):
+ self.LOG.debug('Calling _multiply_jacobi')
v_dim, u_dim = self.dim_uv
size = np.prod(self.dim_uv)
assert len(vector) == 2*size, 'vector size not compatible!'
@@ -187,26 +173,11 @@ class Kernel(object):
return result
def _multiply_jacobi_T(self, vector):
- '''Calculate the product of the transposed Jacobi matrix with a given `vector`.
-
- Parameters
- ----------
- vector : :class:`~numpy.ndarray` (N=1)
- Vectorized form of the magnetization in `u`- and `v`-direction of every pixel
- (row-wise). The first ``N**2`` elements have to correspond to the `u`-, the next
- ``N**2`` elements to the `v`-component of the magnetization.
-
- Returns
- -------
- result : :class:`~numpy.ndarray` (N=1)
- Product of the transposed Jacobi matrix (which is not explicitely calculated) with
- the vector.
-
- '''# TODO: move!
- self.log.info('Calling _multiply_jacobi_T')
+ self.LOG.debug('Calling _multiply_jacobi_T')
v_dim, u_dim = self.dim_uv
size = np.prod(self.dim_uv)
- assert len(vector) == size, 'vector size not compatible! vector: {}, size: {}'.format(len(vector),size)
+ assert len(vector) == size, \
+ 'vector size not compatible! vector: {}, size: {}'.format(len(vector),size)
result = np.zeros(2*size)
for s in range(size): # row-wise (two rows at a time, u- and v-component)
i = s % u_dim
@@ -220,7 +191,7 @@ class Kernel(object):
return result
def _multiply_jacobi_core(self, vector):
- self.log.info('Calling _multiply_jacobi_core')
+ self.LOG.debug('Calling _multiply_jacobi_core')
result = np.zeros(np.prod(self.dim_uv))
core.multiply_jacobi_core(self.dim_uv[0], self.dim_uv[1], self.u, self.v, vector, result)
return result
diff --git a/pyramid/logfile.log b/pyramid/logfile.log
deleted file mode 100644
index e69de29..0000000
diff --git a/pyramid/logging.ini b/pyramid/logging.ini
index 266c375..d206f57 100644
--- a/pyramid/logging.ini
+++ b/pyramid/logging.ini
@@ -23,7 +23,7 @@ args=tuple()
[handler_file]
class=logging.FileHandler
-level=INFO
+level=DEBUG
formatter=file
args=('logfile.log', 'w')
diff --git a/pyramid/magcreator.py b/pyramid/magcreator.py
index 82623a1..5608b72 100644
--- a/pyramid/magcreator.py
+++ b/pyramid/magcreator.py
@@ -13,14 +13,20 @@ center.
"""
-
import numpy as np
from numpy import pi
+import abc
+
+import logging
+
+
+LOG = logging.getLogger(__name__)
+
class Shapes(object):
- '''Class containing functions for generating magnetic shapes.
+ '''Abstract class containing functions for generating magnetic shapes.
The :class:`~.Shapes` class is a collection of some methods that return a 3-dimensional
matrix that represents the magnetized volume and consists of values between 0 and 1.
@@ -28,7 +34,9 @@ class Shapes(object):
:class:`~pyramid.magdata.MagData` objects which store the magnetic informations.
'''
- #TODO: Abstract class (test if this works!)
+
+ __metaclass__ = abc.ABCMeta
+ LOG = logging.getLogger(__name__+'.Shapes')
@classmethod
def slab(cls, dim, center, width):
@@ -49,6 +57,7 @@ class Shapes(object):
The magnetic shape as a 3D-array with values between 1 and 0.
'''
+ cls.LOG.debug('Calling slab')
assert np.shape(dim) == (3,), 'Parameter dim has to be a tuple of length 3!'
assert np.shape(center) == (3,), 'Parameter center has to be a tuple of length 3!'
assert np.shape(width) == (3,), 'Parameter width has to be a tuple of length 3!'
@@ -83,11 +92,12 @@ class Shapes(object):
The magnetic shape as a 3D-array with values between 1 and 0.
'''
+ cls.LOG.debug('Calling disc')
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 radius > 0 and np.shape(radius) == (), 'Radius has to be a positive scalar value!'
assert height > 0 and np.shape(height) == (), 'Height has to be a positive scalar value!'
- assert axis == 'z' or axis == 'y' or axis == 'x', 'Axis has to be x, y or z (as a string)!'
+ assert axis in {'z', 'y', 'x'}, 'Axis has to be x, y or z (as a string)!'
if axis == 'z':
mag_shape = np.array([[[np.hypot(x - center[2], y - center[1]) <= radius
and abs(z - center[0]) <= height / 2
@@ -108,6 +118,58 @@ class Shapes(object):
for z in range(dim[0])])
return mag_shape
+ @classmethod
+ def ellipse(cls, dim, center, width, height, axis='z'):
+ '''Create the shape of an elliptical cylinder in x-, y-, or z-direction.
+
+ Parameters
+ ----------
+ dim : tuple (N=3)
+ The dimensions of the grid `(z, y, x)`.
+ center : tuple (N=3)
+ The center of the ellipse in pixel coordinates `(z, y, x)`.
+ width : tuple (N=2)
+ Length of the two axes of the ellipse.
+ height : float
+ The height of the ellipse in pixel coordinates.
+ axis : {'z', 'y', 'x'}, optional
+ The orientation of the ellipse axis. The default is 'z'.
+
+ Returns
+ -------
+ mag_shape : :class:`~numpy.ndarray` (N=3)
+ The magnetic shape as a 3D-array with values between 1 and 0.
+
+ '''
+ cls.LOG.debug('Calling ellipse')
+ 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) == (2,), 'Parameter width has to be a a tuple of length 2!'
+ assert height > 0 and np.shape(height) == (), 'Height has to be a positive scalar value!'
+ assert axis in {'z', 'y', 'x'}, 'Axis has to be x, y or z (as a string)!'
+ if axis == 'z':
+ mag_shape = np.array([[[np.hypot((x-center[2])/(width[1]/2.),
+ (y-center[1])/(width[0]/2.))<=1
+ and abs(z - center[0]) <= height / 2
+ for x in range(dim[2])]
+ for y in range(dim[1])]
+ for z in range(dim[0])])
+ elif axis == 'y':
+ mag_shape = np.array([[[np.hypot((x-center[2])/(width[1]/2.),
+ (z-center[0])/(width[0]/2.))<=1
+ and abs(y-center[1]) <= height / 2
+ for x in range(dim[2])]
+ for y in range(dim[1])]
+ for z in range(dim[0])])
+ 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
+ for x in range(dim[2])]
+ for y in range(dim[1])]
+ for z in range(dim[0])])
+ return mag_shape
+
@classmethod
def sphere(cls, dim, center, radius):
'''Create the shape of a sphere.
@@ -127,6 +189,7 @@ class Shapes(object):
The magnetic shape as a 3D-array with values between 1 and 0.
'''
+ cls.LOG.debug('Calling sphere')
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 radius > 0 and np.shape(radius) == (), 'Radius has to be a positive scalar value!'
@@ -147,7 +210,7 @@ class Shapes(object):
dim : tuple (N=3)
The dimensions of the grid `(z, y, x)`.
center : tuple (N=3)
- The center of the sphere in pixel coordinates `(z, y, x)`.
+ The center of the ellipsoid in pixel coordinates `(z, y, x)`.
width : tuple (N=3)
The width of the ellipsoid `(z, y, x)`.
@@ -157,6 +220,7 @@ class Shapes(object):
The magnetic shape as a 3D-array with values between 1 and 0.
'''
+ cls.LOG.debug('Calling ellipsoid')
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!'
@@ -189,9 +253,10 @@ class Shapes(object):
The magnetic shape as a 3D-array with values between 1 and 0.
'''
+ cls.LOG.debug('Calling filament')
assert np.shape(dim) == (3,), 'Parameter dim has to be a tuple of length 3!'
assert np.shape(pos) == (2,), 'Parameter pos has to be a tuple of length 2!'
- assert axis == 'z' or axis == 'y' or axis == 'x', 'Axis has to be x, y or z (as a string)!'
+ assert axis in {'z', 'y', 'x'}, 'Axis has to be x, y or z (as a string)!'
mag_shape = np.zeros(dim)
if axis == 'z':
mag_shape[:, pos[0], pos[1]] = 1
@@ -218,6 +283,7 @@ class Shapes(object):
The magnetic shape as a 3D-array with values between 1 and 0.
'''
+ cls.LOG.debug('Calling pixel')
assert np.shape(dim) == (3,), 'Parameter dim has to be a tuple of length 3!'
assert np.shape(pixel) == (3,), 'Parameter pixel has to be a tuple of length 3!'
mag_shape = np.zeros(dim)
@@ -226,7 +292,7 @@ class Shapes(object):
def create_mag_dist_homog(mag_shape, phi, theta=pi/2, magnitude=1):
- '''Create a 3-dimensional magnetic distribution of a homogeneous magnetized object.
+ '''Create a 3-dimensional magnetic distribution of a homogeneously magnetized object.
Parameters
----------
@@ -247,15 +313,16 @@ def create_mag_dist_homog(mag_shape, phi, theta=pi/2, magnitude=1):
`x`-, `y`- and `z`-direction on the 3-dimensional grid.
'''
+ LOG.debug('Calling create_mag_dist_homog')
dim = np.shape(mag_shape)
assert len(dim) == 3, 'Magnetic shapes must describe 3-dimensional distributions!'
z_mag = np.ones(dim) * np.cos(theta) * mag_shape * magnitude
y_mag = np.ones(dim) * np.sin(theta) * np.sin(phi) * mag_shape * magnitude
x_mag = np.ones(dim) * np.sin(theta) * np.cos(phi) * mag_shape * magnitude
- return np.array([x_mag, y_mag, z_mag]) # TODO: Better implementation!
+ return np.array([x_mag, y_mag, z_mag])
-def create_mag_dist_vortex(mag_shape, center=None, magnitude=1):
+def create_mag_dist_vortex(mag_shape, center=None, axis='z', magnitude=1):
'''Create a 3-dimensional magnetic distribution of a homogeneous magnetized object.
Parameters
@@ -263,12 +330,14 @@ def create_mag_dist_vortex(mag_shape, center=None, magnitude=1):
mag_shape : :class:`~numpy.ndarray` (N=3)
The magnetic shapes (see :class:`.~Shapes` for examples).
center : tuple (N=2 or N=3), optional
- The vortex center, given in 2D `(x, y)` or 3D `(x, y, z)`, where `z` is discarded.
- Is set to the center of the field of view if not specified.
- Vortex center has to be between two pixels.
+ The vortex center, given in 2D `(v, u)` or 3D `(z, y, x)`, where the perpendicular axis
+ is is discarded. Is set to the center of the field of view if not specified. The vortex
+ center has to be between two pixels.
magnitude : float, optional
- The relative magnitude for the magnetic shape. The default is 1.
-
+ The relative magnitude for the magnetic shape. The default is 1. Negative signs reverse
+ the vortex direction (left-handed instead of right-handed).
+ axis : {'z', 'y', 'x'}, optional
+ The orientation of the vortex axis. The default is 'z'.
Returns
-------
magnitude : tuple (N=3) of :class:`~numpy.ndarray` (N=3)
@@ -276,18 +345,45 @@ def create_mag_dist_vortex(mag_shape, center=None, magnitude=1):
`x`-, `y`- and `z`-direction on the 3-dimensional grid.
'''
+ LOG.debug('Calling create_mag_dist_vortex')
dim = np.shape(mag_shape)
assert len(dim) == 3, 'Magnetic shapes must describe 3-dimensional distributions!'
+ assert axis in {'z', 'y', 'x'}, 'Axis has to be x, y or z (as a string)!'
+ assert center is None or len(center) in {2, 3}, \
+ 'Vortex center has to be defined in 3D or 2D or not at all!'
if center is None:
center = (int(dim[1]/2)-0.5, int(dim[2]/2)-0.5) # center has to be between (!) two pixels
- if len(center) == 3: # if a 3D-center is given, just take the x and y components
- center = (center[1], center[2])
- assert len(center) == 2, 'Vortex center has to be defined in 3D or 2D or not at all!'
- x = np.linspace(-center[1], dim[2]-1-center[1], dim[2])
- y = np.linspace(-center[0], dim[1]-1-center[0], dim[1])
- xx, yy = np.meshgrid(x, y)
- phi = np.arctan2(yy, xx) - pi/2
- z_mag = np.zeros(dim)
- y_mag = -np.ones(dim) * np.sin(phi) * mag_shape * magnitude
- x_mag = -np.ones(dim) * np.cos(phi) * mag_shape * magnitude
- return np.array([x_mag, y_mag, z_mag]) # TODO: Better implementation!
+ if axis == 'z':
+ if len(center) == 3: # if a 3D-center is given, just take the x and y components
+ center = (center[1], center[2])
+ u = np.linspace(-center[1], dim[2]-1-center[1], dim[2])
+ v = np.linspace(-center[0], dim[1]-1-center[0], dim[1])
+ uu, vv = np.meshgrid(u, v)
+ phi = np.expand_dims(np.arctan2(vv, uu)-pi/2, axis=0)
+ phi = np.tile(phi, (dim[0], 1, 1))
+ z_mag = np.zeros(dim)
+ y_mag = -np.ones(dim) * np.sin(phi) * mag_shape * magnitude
+ x_mag = -np.ones(dim) * np.cos(phi) * mag_shape * magnitude
+ elif axis == 'y':
+ if len(center) == 3: # if a 3D-center is given, just take the x and z components
+ center = (center[0], center[2])
+ u = np.linspace(-center[1], dim[2]-1-center[1], dim[2])
+ v = np.linspace(-center[0], dim[0]-1-center[0], dim[0])
+ uu, vv = np.meshgrid(u, v)
+ phi = np.expand_dims(np.arctan2(vv, uu)-pi/2, axis=1)
+ phi = np.tile(phi, (1, dim[1], 1))
+ z_mag = np.ones(dim) * np.sin(phi) * mag_shape * magnitude
+ y_mag = np.zeros(dim)
+ x_mag = np.ones(dim) * np.cos(phi) * mag_shape * magnitude
+ elif axis == 'x':
+ if len(center) == 3: # if a 3D-center is given, just take the z and y components
+ center = (center[0], center[1])
+ u = np.linspace(-center[1], dim[1]-1-center[1], dim[1])
+ v = np.linspace(-center[0], dim[0]-1-center[0], dim[0])
+ uu, vv = np.meshgrid(u, v)
+ phi = np.expand_dims(np.arctan2(vv, uu)-pi/2, axis=2)
+ phi = np.tile(phi, (1, 1, dim[2]))
+ z_mag = -np.ones(dim) * np.sin(phi) * mag_shape * magnitude
+ 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])
diff --git a/pyramid/magdata.py b/pyramid/magdata.py
index 0dc7044..d51ed6e 100644
--- a/pyramid/magdata.py
+++ b/pyramid/magdata.py
@@ -2,7 +2,6 @@
"""This module provides the :class:`~.MagData` class for storing of magnetization data."""
-import logging
import numpy as np
from scipy.ndimage.interpolation import zoom
@@ -14,6 +13,8 @@ from numbers import Number
import netCDF4
+import logging
+
class MagData(object):
@@ -21,7 +22,7 @@ class MagData(object):
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
+ `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.
@@ -40,9 +41,9 @@ class MagData(object):
Vector containing the magnetic distribution.
'''
-
- log = logging.getLogger()
-
+
+ LOG = logging.getLogger(__name__+'.MagData')
+
@property
def a(self):
return self._a
@@ -51,8 +52,8 @@ class MagData(object):
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 = a
-
+ self._a = float(a)
+
@property
def dim(self):
return self._dim
@@ -68,7 +69,7 @@ class MagData(object):
assert magnitude.shape[0] == 3, 'First dimension of the magnitude has to be 3!'
self._magnitude = magnitude
self._dim = magnitude.shape[1:]
-
+
@property
def mag_vec(self):
return np.reshape(self.magnitude, -1)
@@ -80,80 +81,70 @@ class MagData(object):
self.magnitude = mag_vec.reshape((3,)+self.dim)
def __init__(self, a, magnitude):
- self.log = logging.getLogger(__name__)
- self.log.info('Calling __init__')
+ self.LOG.debug('Calling __init__')
self.a = a
- self.magnitude = magnitude
- self.log.info('Created '+str(self))
-
- def __getstate__(self):
- d = dict(self.__dict__)
- del d['log']
- return d
-
- def __setstate__(self, d):
- self.__dict__.update(d)
- self.log = logging.getLogger(__name__)
+ self.magnitude = magnitude
+ self.LOG.debug('Created '+str(self))
def __repr__(self):
- self.log.info('Calling __repr__')
+ self.LOG.debug('Calling __repr__')
return '%s(a=%r, magnitude=%r)' % (self.__class__, self.a, self.magnitude)
def __str__(self):
- self.log.info('Calling __str__')
+ self.LOG.debug('Calling __str__')
return 'MagData(a=%s, dim=%s)' % (self.a, self.dim)
def __neg__(self): # -self
- self.log.info('Calling __neg__')
+ self.LOG.debug('Calling __neg__')
return MagData(self.a, -self.magnitude)
-
+
def __add__(self, other): # self + other
- self.log.info('Calling __add__')
+ 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.info('Adding two MagData objects')
+ 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.info('Adding an offset')
+ self.LOG.debug('Adding an offset')
return MagData(self.a, self.magnitude+other)
def __sub__(self, other): # self - other
- self.log.info('Calling __sub__')
+ self.LOG.debug('Calling __sub__')
return self.__add__(-other)
def __mul__(self, other): # self * other
- self.log.info('Calling __mul__')
+ 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.info('Calling __radd__')
+ self.LOG.debug('Calling __radd__')
return self.__add__(other)
def __rsub__(self, other): # other - self
- self.log.info('Calling __rsub__')
+ self.LOG.debug('Calling __rsub__')
return -self.__sub__(other)
def __rmul__(self, other): # other * self
- self.log.info('Calling __rmul__')
+ self.LOG.debug('Calling __rmul__')
return self.__mul__(other)
-
+
def __iadd__(self, other): # self += other
- self.log.info('Calling __iadd__')
+ self.LOG.debug('Calling __iadd__')
return self.__add__(other)
def __isub__(self, other): # self -= other
- self.log.info('Calling __isub__')
+ self.LOG.debug('Calling __isub__')
return self.__sub__(other)
def __imul__(self, other): # self *= other
- self.log.info('Calling __imul__')
+ self.LOG.debug('Calling __imul__')
return self.__mul__(other)
-
+
def copy(self):
'''Returns a copy of the :class:`~.MagData` object
@@ -167,7 +158,7 @@ class MagData(object):
A copy of the :class:`~.MagData`.
'''
- self.log.info('Calling copy7')
+ self.LOG.debug('Calling copy')
return MagData(self.a, self.magnitude.copy())
def scale_down(self, n=1):
@@ -188,7 +179,7 @@ class MagData(object):
Only possible, if each axis length is a power of 2!
'''
- self.log.info('Calling scale_down')
+ self.LOG.debug('Calling scale_down')
assert n > 0 and isinstance(n, (int, long)), 'n must be a positive integer!'
assert np.all([d % 2**n == 0 for d in self.dim]), 'Dimensions must a multiples of 2!'
self.a = self.a * 2**n
@@ -196,11 +187,29 @@ class MagData(object):
# 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):
- # TODO: Docstring!
- self.log.info('Calling scale_up')
+ '''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!'
@@ -210,7 +219,28 @@ class MagData(object):
zoom(self.magnitude[2], zoom=2**n, order=order)))
def pad(self, x_pad, y_pad, z_pad):
- # TODO: Docstring!
+ '''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', constant_values=0)
@@ -229,7 +259,7 @@ class MagData(object):
None
'''
- self.log.info('Calling save_to_llg')
+ 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,
@@ -259,7 +289,7 @@ class MagData(object):
A :class:`~.MagData` object containing the loaded data.
'''
- cls.log.info('Calling load_from_llg')
+ 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])))
@@ -281,7 +311,7 @@ class MagData(object):
None
'''
- self.log.info('Calling save_to_netcdf4')
+ 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
@@ -307,7 +337,7 @@ class MagData(object):
A :class:`~.MagData` object containing the loaded data.
'''
- cls.log.info('Calling copy')
+ cls.LOG.debug('Calling copy')
mag_file = netCDF4.Dataset(filename, 'r', format='NETCDF4')
a = mag_file.a
magnitude = mag_file.variables['magnitude'][...]
@@ -339,31 +369,31 @@ class MagData(object):
The axis on which the graph is plotted.
'''
- self.log.info('Calling quiver_plot')
+ 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.info('proj_axis == z')
+ self.LOG.debug('proj_axis == z')
if ax_slice is None:
- self.log.info('ax_slice is None')
+ self.LOG.debug('ax_slice is None')
ax_slice = int(self.dim[0]/2)
mag_slice_u = self.magnitude[0][ax_slice, ...] # x-component
mag_slice_v = 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.info('proj_axis == y')
+ self.LOG.debug('proj_axis == y')
if ax_slice is None:
- self.log.info('ax_slice is None')
+ self.LOG.debug('ax_slice is None')
ax_slice = int(self.dim[1]/2)
mag_slice_u = self.magnitude[0][:, ax_slice, :] # x-component
mag_slice_v = 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.info('proj_axis == x')
+ self.LOG.debug('proj_axis == x')
if ax_slice is None:
- self.log.info('ax_slice is None')
+ self.LOG.debug('ax_slice is None')
ax_slice = int(self.dim[2]/2)
mag_slice_u = self.magnitude[1][..., ax_slice] # y-component
mag_slice_v = self.magnitude[2][..., ax_slice] # z-component
@@ -371,7 +401,7 @@ class MagData(object):
v_label = 'z [px]'
# If no axis is specified, a new figure is created:
if axis is None:
- self.log.info('axis is None')
+ self.LOG.debug('axis is None')
fig = plt.figure(figsize=(8.5, 7))
axis = fig.add_subplot(1, 1, 1, aspect='equal')
axis.quiver(mag_slice_u, mag_slice_v, pivot='middle', angles='xy', scale_units='xy',
@@ -399,7 +429,7 @@ class MagData(object):
None
'''
- self.log.info('Calling quiver_plot3D')
+ self.LOG.debug('Calling quiver_plot3D')
a = self.a
dim = self.dim
# Create points and vector components as lists:
@@ -417,4 +447,5 @@ class MagData(object):
plot = mlab.quiver3d(xx, yy, zz, x_mag, y_mag, z_mag, mode='arrow')
mlab.outline(plot)
mlab.axes(plot)
- mlab.colorbar()
+ mlab.title('TEST', height=0.95, size=0.25)
+ mlab.colorbar(plot, orientation='vertical')
diff --git a/pyramid/numcore/__init__.py b/pyramid/numcore/__init__.py
index fd8ed34..c462822 100644
--- a/pyramid/numcore/__init__.py
+++ b/pyramid/numcore/__init__.py
@@ -1,3 +1,4 @@
+# -*- coding: utf-8 -*-
"""Package for numerical core routines which enable fast computations.
Modules
@@ -10,4 +11,4 @@ Notes
-----
Packages are written in `pyx`-format for the use with :mod:`~cython`.
-"""
+"""
\ No newline at end of file
diff --git a/pyramid/numcore/kernel_core.pyx b/pyramid/numcore/kernel_core.pyx
index 7ed2d6e..f0a7412 100644
--- a/pyramid/numcore/kernel_core.pyx
+++ b/pyramid/numcore/kernel_core.pyx
@@ -1,11 +1,10 @@
# -*- coding: utf-8 -*-
-"""Numerical core routines for the phase calculation using the real space approach.
+"""Numerical core routines for the :mod:`~.pyramid.kernel` module.
-Provides a helper function to speed up :func:`~pyramid.phasemapper.phase_mag_real` of module
-:mod:`~pyramid.phasemapper`, by using C-speed for the for-loops and by omitting boundary and
-wraparound checks.
+Provides helper functions to speed up kernel calculations in the :class:`~pyramid.kernel.Kernel`
+class, by using C-speed for the for-loops and by omitting boundary and wraparound checks.
-""" # TODO: Docstring!
+"""
import numpy as np
@@ -17,12 +16,12 @@ cimport numpy as np
@cython.boundscheck(False)
@cython.wraparound(False)
-def phase_mag_real_core(
+def multiply_jacobi_core(
unsigned int v_dim, unsigned int u_dim,
- double[:, :] v_phi, double[:, :] u_phi,
- double[:, :] v_mag, double[:, :] u_mag,
- double[:, :] phase, float threshold):
- '''Numerical core routine for the phase calculation using the real space approach.
+ double[:, :] u_phi, double[:, :] v_phi,
+ double[:] vector,
+ double[:] result):
+ '''Numerical core routine for the Jacobi matrix multiplication in the :class:`~.Kernel` class.
Parameters
----------
@@ -30,106 +29,21 @@ def phase_mag_real_core(
Dimensions of the projection along the two major axes.
v_phi, u_phi : :class:`~numpy.ndarray` (N=2)
Lookup tables for the pixel fields oriented in `u`- and `v`-direction.
- v_mag, u_mag : :class:`~numpy.ndarray` (N=2)
- Magnetization components in `u`- and `v`-direction.
- phase : :class:`~numpy.ndarray` (N=2)
- Matrix in which the resulting magnetic phase map should be stored.
- threshold : float
- The `threshold` determines which pixels contribute to the magnetic phase.
+ vector : :class:`~numpy.ndarray` (N=1)
+ Input vector which should be multiplied by the Jacobi-matrix.
+ result : :class:`~numpy.ndarray` (N=1)
+ Vector in which the result of the multiplication should be stored.
Returns
-------
None
- '''
- cdef unsigned int i, j, p, q, p_c, q_c
- cdef double u_m, v_m
- for j in range(v_dim):
- for i in range(u_dim):
- u_m = u_mag[j, i]
- v_m = v_mag[j, i]
- p_c = u_dim - 1 - i
- q_c = v_dim - 1 - j
- if abs(u_m) > threshold:
- for q in range(v_dim):
- for p in range(u_dim):
- phase[q, p] += u_m * u_phi[q_c + q, p_c + p]
- if abs(v_m) > threshold:
- for q in range(v_dim):
- for p in range(u_dim):
- phase[q, p] -= v_m * v_phi[q_c + q, p_c + p]
-
-
-@cython.boundscheck(False)
-@cython.wraparound(False)
-def get_jacobi_core(
- unsigned int v_dim, unsigned int u_dim,
- double[:, :] v_phi, double[:, :] u_phi,
- double[:, :] jacobi):
- '''Numerical core routine for the jacobi matrix calculation.
-
- Parameters
- ----------
- v_dim, u_dim : int
- Dimensions of the projection along the two major axes.
- v_phi, u_phi : :class:`~numpy.ndarray` (N=2)
- Lookup tables for the pixel fields oriented in `u`- and `v`-direction.
- jacobi : :class:`~numpy.ndarray` (N=2)
- Jacobi matrix which is filled by this routine.
-
- Returns
- -------
- None
+ Notes
+ -----
+ The strategy involves iterating over the magnetization first and over the kernel in an inner
+ iteration loop.
'''
- cdef unsigned int i, j, p, q, p_min, p_max, q_min, q_max, u_column, v_column, row
- for j in range(v_dim):
- for i in range(u_dim):
- q_min = (v_dim-1) - j
- q_max = (2*v_dim-1) - j
- p_min = (u_dim-1) - i
- p_max = (2*u_dim-1) - i
- v_column = i + u_dim*j + u_dim*v_dim
- u_column = i + u_dim*j
- for q in range(q_min, q_max):
- for p in range(p_min, p_max):
- q_c = q - (v_dim-1-j) # start at zero for jacobi column
- p_c = p - (u_dim-1-i) # start at zero for jacobi column
- row = p_c + u_dim*q_c
- # u-component:
- jacobi[row, u_column] = u_phi[q, p]
- # v-component (note the minus!):
- jacobi[row, v_column] = -v_phi[q, p]
-
-
-@cython.boundscheck(False)
-@cython.wraparound(False)
-def get_jacobi_core2(
- unsigned int v_dim, unsigned int u_dim,
- double[:, :] v_phi, double[:, :] u_phi,
- double[:, :] jacobi):
- '''DOCSTRING!'''
- # TODO: Docstring!!!
- cdef unsigned int i, j, p, q, p_min, p_max, q_min, q_max, u_column, v_column, row
- for j in range(v_dim):
- for i in range(u_dim):
- # v_dim*u_dim columns for the u-component:
- jacobi[:, i+u_dim*j] = \
- np.reshape(u_phi[v_dim-1-j:(2*v_dim-1)-j, u_dim-1-i:(2*u_dim-1)-i], -1)
- # v_dim*u_dim columns for the v-component (note the minus!):
- jacobi[:, u_dim*v_dim+i+u_dim*j] = \
- -np.reshape(v_phi[v_dim-1-j:(2*v_dim-1)-j, u_dim-1-i:(2*u_dim-1)-i], -1)
-
-
-@cython.boundscheck(False)
-@cython.wraparound(False)
-def multiply_jacobi_core(
- unsigned int v_dim, unsigned int u_dim,
- double[:, :] u_phi, double[:, :] v_phi,
- double[:] vector,
- double[:] result):
- '''DOCSTRING!'''
- # TODO: Docstring!!! Iterate over magnetization and then kernel
cdef unsigned int s, i, j, u_min, u_max, v_min, v_max, ri, u, v, size
size = u_dim * v_dim # Number of pixels
s = 0 # Current pixel (numbered consecutively)
@@ -142,38 +56,7 @@ def multiply_jacobi_core(
# Go over the current kernel cutout [v_min:v_max, u_min:u_max]:
for v in range(v_min, v_min + v_dim):
for u in range(u_min, u_min + u_dim):
- result[ri] += vector[s] * u_phi[v, u]
+ result[ri] += vector[s] * u_phi[v, u]
result[ri] -= vector[s+size] * v_phi[v, u]
ri += 1
s += 1
-
-
-#@cython.boundscheck(False)
-#@cython.wraparound(False)
-def multiply_jacobi_core2(
- int v_dim, int u_dim,
- double[:, :] u_phi, double[:, :] v_phi,
- double[:] vector,
- double[:] result):
- '''DOCSTRING!'''
- # TODO: Docstring!!! Iterate over kernel and then magnetization
- cdef int s1, s2, i, j, u_min, u_max, v_min, v_max, ri, u, v, size, j_min, j_max, i_min, i_max
- size = u_dim * v_dim # Number of pixels
- # Go over complete kernel:
- for v in range(2 * v_dim - 1):
- for u in range(2 * u_dim - 1):
- i_min = max(0, (u_dim - 1) - u)
- i_max = min(u_dim, ((2 * u_dim) - 1) - u)
- j_min = max(0, (v_dim - 1) - v)
- j_max = min(v_dim, ((2 * v_dim) - 1) - v)
- # Iterate over
- for i in range(i_min, i_max):
- s1 = i * u_dim # u-column
- s2 = s1 + size # v-column
- u_min = (u_dim - 1) - i
- v_min = (v_dim - 1) - j_min
- ri = (v - ((v_dim - 1) - j_min)) * u_dim + (u - u_min)
- for j in range(j_min, j_max):
- result[ri] += vector[s1 + j] * u_phi[v, u]
- result[ri] -= vector[s2 + j] * v_phi[v, u]
- ri += u_dim
diff --git a/pyramid/optimizer.py b/pyramid/optimizer.py
deleted file mode 100644
index f45546e..0000000
--- a/pyramid/optimizer.py
+++ /dev/null
@@ -1,143 +0,0 @@
-# -*- coding: utf-8 -*-
-"""Reconstruct magnetic distributions from given phasemaps.
-
-This module reconstructs 3-dimensional magnetic distributions (as :class:`~pyramid.magdata.MagData`
-objects) from a given set of phase maps (represented by :class:`~pyramid.phasemap.PhaseMap`
-objects) by using several model based reconstruction algorithms which use the forward model
-provided by :mod:`~pyramid.projector` and :mod:`~pyramid.phasemapper` and a priori knowledge of
-the distribution.
-So far, only a simple least square algorithm for known pixel locations for 2-dimensional problems
-is implemented (:func:`~.reconstruct_simple_leastsq`), but more complex solutions are planned.
-
-"""# TODO: Docstrings!
-
-
-# TODO: Delete
-''' FRAGEN FÜR JÖRN:
-1) Property und Python "Magic" ist wann sinnvoll?
-2) Aktuelle Designstruktur, welche Designpatterns kommen vor und sind sinnvoll?
-3) Macht logging Sinn?
-4) Ist scipy.optimize.minimize(...) eine gute Wahl?
-5) Wo besser Module, wo Klassen?
-6) Arbitrary grid to cartesian grid
-'''
-
-
-import numpy as np
-
-from scipy.optimize import minimize
-
-from pyramid.kernel import Kernel
-from pyramid.forwardmodel import ForwardModel
-from pyramid.costfunction import Costfunction
-from pyramid.magdata import MagData
-
-
-
-def optimize_cg(data_collection, first_guess):
- # TODO: where should data_collection and first_guess go? here or __init__?
- data = data_collection # TODO: name the parameter 'data' directly...
- mag_0 = first_guess
- x_0 = first_guess.mag_vec
- y = data.phase_vec
- kern = Kernel(data.a, data.dim_uv)
- F = ForwardModel(data.projectors, kern, data.b_0)
- C = Costfunction(y, F)
-
- result = minimize(C, x_0, method='Newton-CG', jac=C.jac, hessp=C.hess_dot,
- options={'maxiter':200, 'disp':True})
-
- x_opt = result.x
-
- mag_opt = MagData(mag_0.a, np.zeros((3,)+mag_0.dim))
-
- mag_opt.mag_vec = x_opt
-
- return mag_opt
-
-
-
-# TODO: Implement the following:
-
-## -*- coding: utf-8 -*-
-#"""Reconstruct magnetic distributions from given phasemaps.
-#
-#This module reconstructs 3-dimensional magnetic distributions (as :class:`~pyramid.magdata.MagData`
-#objects) from a given set of phase maps (represented by :class:`~pyramid.phasemap.PhaseMap`
-#objects) by using several model based reconstruction algorithms which use the forward model
-#provided by :mod:`~pyramid.projector` and :mod:`~pyramid.phasemapper` and a priori knowledge of
-#the distribution.
-#So far, only a simple least square algorithm for known pixel locations for 2-dimensional problems
-#is implemented (:func:`~.reconstruct_simple_leastsq`), but more complex solutions are planned.
-#
-#"""
-#
-#
-#
-#import numpy as np
-#
-#from scipy.optimize import leastsq
-#
-#import pyramid.projector as pj
-#import pyramid.phasemapper as pm
-#from pyramid.magdata import MagData
-#from pyramid.projector import Projection
-#from pyramid.kernel import Kernel
-#
-#
-#def reconstruct_simple_leastsq(phase_map, mask, b_0=1):
-# '''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.
-# Returns
-# -------
-# mag_data : :class:`~pyramid.magdata.MagData`
-# The reconstructed magnetic distribution as a
-# :class:`~pyramid.magdata.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.reshape(-1) # 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
-# lam = 1e-6 # Regularisation parameter
-# # Create empty MagData object for the reconstruction:
-# mag_data_rec = MagData(a, (np.zeros(dim), np.zeros(dim), np.zeros(dim)))
-#
-# # Function that returns the phase map for a magnetic configuration x:
-# def F(x):
-# mag_data_rec.set_vector(mask, x)
-# phase = pm.phase_mag_real(a, pj.simple_axis_projection(mag_data_rec), b_0)
-# return phase.reshape(-1)
-#
-# # Cost function which should be minimized:
-# def J(x_i):
-# y_i = F(x_i)
-# term1 = (y_i - y_m)
-# term2 = lam * x_i
-# return np.concatenate([term1, term2])
-#
-# # Reconstruct the magnetization components:
-# x_rec, _ = leastsq(J, np.zeros(3*count))
-# mag_data_rec.set_vector(mask, x_rec)
-# return mag_data_rec
-#
-#def reconstruct_test():
-# product = (kernel.multiply_jacobi_T(projection.multiply_jacobi_T(x))
-# * kernel.multiply_jacobi(projection.multiply_jacobi(x)))
\ No newline at end of file
diff --git a/pyramid/phasemap.py b/pyramid/phasemap.py
index be2b294..08f2b01 100644
--- a/pyramid/phasemap.py
+++ b/pyramid/phasemap.py
@@ -2,8 +2,6 @@
"""This module provides the :class:`~.PhaseMap` class for storing phase map data."""
-import logging
-
import numpy as np
from numpy import pi
@@ -17,6 +15,8 @@ from numbers import Number
import netCDF4
+import logging
+
class PhaseMap(object):
@@ -24,7 +24,7 @@ class PhaseMap(object):
Represents 2-dimensional phase maps. The phase information itself is stored as a 2-dimensional
matrix in `phase`, but can also be accessed as a vector via `phase_vec`. :class:`~.PhaseMap`
- objects support negation, arithmetic operators (``+``, ``-``, ``*``) and their augmented
+ objects support negation, arithmetic operators (``+``, ``-``, ``*``) and their augmented
counterparts (``+=``, ``-=``, ``*=``), with numbers and other :class:`~.PhaseMap`
objects, if their dimensions and grid spacings match. It is possible to load data from NetCDF4
or textfiles or to save the data in these formats. Methods for plotting the phase or a
@@ -51,7 +51,7 @@ class PhaseMap(object):
'''
- log = logging.getLogger(__name__)
+ LOG = logging.getLogger(__name__)
UNITDICT = {u'rad': 1E0,
u'mrad': 1E3,
@@ -92,10 +92,10 @@ class PhaseMap(object):
(0.50, 1.0, 1.0),
(0.75, 1.0, 1.0),
(1.00, 0.0, 0.0)]}
-
+
HOLO_CMAP = mpl.colors.LinearSegmentedColormap('my_colormap', CDICT, 256)
HOLO_CMAP_INV = mpl.colors.LinearSegmentedColormap('my_colormap', CDICT_INV, 256)
-
+
@property
def a(self):
return self._a
@@ -104,8 +104,8 @@ class PhaseMap(object):
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 = a
-
+ self._a = float(a)
+
@property
def dim_uv(self):
return self._dim_uv
@@ -141,70 +141,70 @@ class PhaseMap(object):
self._unit = unit
def __init__(self, a, phase, unit='rad'):
+ self.LOG.debug('Calling __init__')
self.a = a
self.phase = phase
self.unit = unit
- self.log = logging.getLogger(__name__)
- self.log.info('Created '+str(self))
+ self.LOG.debug('Created '+str(self))
def __repr__(self):
- self.log.info('Calling __repr__')
+ self.LOG.debug('Calling __repr__')
return '%s(a=%r, phase=%r, unit=&r)' % \
(self.__class__, self.a, self.phase, self.unit)
def __str__(self):
- self.log.info('Calling __str__')
+ self.LOG.debug('Calling __str__')
return 'PhaseMap(a=%s, dim_uv=%s)' % (self.a, self.dim_uv)
def __neg__(self): # -self
- self.log.info('Calling __neg__')
+ self.LOG.debug('Calling __neg__')
return PhaseMap(self.a, -self.phase, self.unit)
-
+
def __add__(self, other): # self + other
- self.log.info('Calling __add__')
+ self.LOG.debug('Calling __add__')
assert isinstance(other, (PhaseMap, Number)), \
'Only PhaseMap objects and scalar numbers (as offsets) can be added/subtracted!'
if isinstance(other, PhaseMap):
- self.log.info('Adding two PhaseMap objects')
+ self.LOG.debug('Adding two PhaseMap objects')
assert other.a == self.a, 'Added phase has to have the same grid spacing!'
assert other.phase.shape == self.dim_uv, \
'Added magnitude has to have the same dimensions!'
return PhaseMap(self.a, self.phase+other.phase, self.unit)
else: # other is a Number
- self.log.info('Adding an offset')
+ self.LOG.debug('Adding an offset')
return PhaseMap(self.a, self.phase+other, self.unit)
def __sub__(self, other): # self - other
- self.log.info('Calling __sub__')
+ self.LOG.debug('Calling __sub__')
return self.__add__(-other)
def __mul__(self, other): # self * other
- self.log.info('Calling __mul__')
+ self.LOG.debug('Calling __mul__')
assert isinstance(other, Number), 'PhaseMap objects can only be multiplied by numbers!'
return PhaseMap(self.a, other*self.phase, self.unit)
def __radd__(self, other): # other + self
- self.log.info('Calling __radd__')
+ self.LOG.debug('Calling __radd__')
return self.__add__(other)
def __rsub__(self, other): # other - self
- self.log.info('Calling __rsub__')
+ self.LOG.debug('Calling __rsub__')
return -self.__sub__(other)
def __rmul__(self, other): # other * self
- self.log.info('Calling __rmul__')
+ self.LOG.debug('Calling __rmul__')
return self.__mul__(other)
-
+
def __iadd__(self, other): # self += other
- self.log.info('Calling __iadd__')
+ self.LOG.debug('Calling __iadd__')
return self.__add__(other)
def __isub__(self, other): # self -= other
- self.log.info('Calling __isub__')
+ self.LOG.debug('Calling __isub__')
return self.__sub__(other)
def __imul__(self, other): # self *= other
- self.log.info('Calling __imul__')
+ self.LOG.debug('Calling __imul__')
return self.__mul__(other)
def save_to_txt(self, filename='..\output\phasemap_output.txt'):
@@ -221,7 +221,7 @@ class PhaseMap(object):
None
'''
- self.log.info('Calling save_to_txt')
+ self.LOG.debug('Calling save_to_txt')
with open(filename, 'w') as phase_file:
phase_file.write('{}\n'.format(filename.replace('.txt', '')))
phase_file.write('grid spacing = {} nm\n'.format(self.a))
@@ -242,7 +242,7 @@ class PhaseMap(object):
A :class:`~.PhaseMap` object containing the loaded data.
'''
- cls.log.info('Calling load_from_txt')
+ cls.LOG.debug('Calling load_from_txt')
with open(filename, 'r') as phase_file:
phase_file.readline() # Headerline is not used
a = float(phase_file.readline()[15:-4])
@@ -263,7 +263,7 @@ class PhaseMap(object):
None
'''
- self.log.info('Calling save_to_netcdf4')
+ self.LOG.debug('Calling save_to_netcdf4')
phase_file = netCDF4.Dataset(filename, 'w', format='NETCDF4')
phase_file.a = self.a
phase_file.createDimension('v_dim', self.dim_uv[0])
@@ -287,7 +287,7 @@ class PhaseMap(object):
A :class:`~.PhaseMap` object containing the loaded data.
'''
- cls.log.info('Calling load_from_netcdf4')
+ cls.LOG.debug('Calling load_from_netcdf4')
phase_file = netCDF4.Dataset(filename, 'r', format='NETCDF4')
a = phase_file.a
phase = phase_file.variables['phase'][:]
@@ -307,7 +307,7 @@ class PhaseMap(object):
The default is 'RdBu'.
limit : float, optional
Plotlimit for the phase in both negative and positive direction (symmetric around 0).
- If not specified, the maximum amplitude the phase is used.
+ If not specified, the maximum amplitude of the phase is used.
norm : :class:`~matplotlib.colors.Normalize` or subclass, optional
Norm, which is used to determine the colors to encode the phase information.
If not specified, :class:`~matplotlib.colors.Normalize` is automatically used.
@@ -322,7 +322,7 @@ class PhaseMap(object):
The axis on which the graph is plotted.
'''
- self.log.info('Calling display_phase')
+ self.LOG.debug('Calling display_phase')
# Take units into consideration:
phase = self.phase * self.UNITDICT[self.unit]
if limit is None:
@@ -374,7 +374,7 @@ class PhaseMap(object):
The axis on which the graph is plotted.
'''
- self.log.info('Calling display_phase3d')
+ self.LOG.debug('Calling display_phase3d')
# Take units into consideration:
phase = self.phase * self.UNITDICT[self.unit]
# Create figure and axis:
@@ -395,11 +395,11 @@ class PhaseMap(object):
plt.show()
# Return plotting axis:
return axis
-
+
def display_holo(self, density=1, title='Holographic Contour Map',
- axis=None, grad_encode='dark', interpolation='none', show=True):
+ axis=None, grad_encode='bright', interpolation='none', show=True):
'''Display the color coded holography image.
-
+
Parameters
----------
density : float, optional
@@ -408,42 +408,51 @@ class PhaseMap(object):
The title of the plot. The default is 'Holographic Contour Map'.
axis : :class:`~matplotlib.axes.AxesSubplot`, optional
Axis on which the graph is plotted. Creates a new figure if none is specified.
+ grad_encode: {'bright', 'dark', 'color', 'none'}, optional
+ Encoding mode of the phase gradient. 'none' produces a black-white image, 'color' just
+ encodes the direction (without gradient strength), 'dark' modulates the gradient
+ strength with a factor between 0 and 1 and 'bright' (which is the default) encodes
+ the gradient strength with color saturation.
interpolation : {'none, 'bilinear', 'cubic', 'nearest'}, optional
Defines the interpolation method. No interpolation is used in the default case.
show: bool, optional
A switch which determines if the plot is shown at the end of plotting.
-
+
Returns
-------
axis: :class:`~matplotlib.axes.AxesSubplot`
The axis on which the graph is plotted.
-
- '''# TODO: Docstring saturation!
- self.log.info('Calling display_holo')
+
+ '''
+ self.LOG.debug('Calling display_holo')
# Calculate the holography image intensity:
img_holo = (1 + np.cos(density * self.phase)) / 2
# Calculate the phase gradients, expressed by magnitude and angle:
phase_grad_y, phase_grad_x = np.gradient(self.phase, self.a, self.a)
phase_angle = (1 - np.arctan2(phase_grad_y, phase_grad_x)/pi) / 2
phase_magnitude = np.hypot(phase_grad_x, phase_grad_y)
- if phase_magnitude.max() != 0:
+ if phase_magnitude.max() != 0: # Take care of phase maps with only zeros
saturation = np.sin(phase_magnitude/phase_magnitude.max() * pi / 2)
phase_saturation = np.dstack((saturation,)*4)
# Color code the angle and create the holography image:
if grad_encode == 'dark':
+ self.LOG.debug('gradient encoding: dark')
rgba = self.HOLO_CMAP(phase_angle)
rgb = (255.999 * img_holo.T * saturation.T * rgba[:, :, :3].T).T.astype(np.uint8)
elif grad_encode == 'bright':
+ self.LOG.debug('gradient encoding: bright')
rgba = self.HOLO_CMAP(phase_angle)+(1-phase_saturation)*self.HOLO_CMAP_INV(phase_angle)
rgb = (255.999 * img_holo.T * rgba[:, :, :3].T).T.astype(np.uint8)
elif grad_encode == 'color':
+ self.LOG.debug('gradient encoding: color')
rgba = self.HOLO_CMAP(phase_angle)
rgb = (255.999 * img_holo.T * rgba[:, :, :3].T).T.astype(np.uint8)
elif grad_encode == 'none':
+ self.LOG.debug('gradient encoding: none')
rgba = self.HOLO_CMAP(phase_angle)+self.HOLO_CMAP_INV(phase_angle)
rgb = (255.999 * img_holo.T * rgba[:, :, :3].T).T.astype(np.uint8)
else:
- raise AssertionError('Gradient encoding not recognized!')
+ raise AssertionError('Gradient encoding not recognized!')
holo_image = Image.fromarray(rgb)
# If no axis is specified, a new figure is created:
if axis is None:
@@ -469,26 +478,41 @@ class PhaseMap(object):
return axis
def display_combined(self, title='Combined Plot', cmap='RdBu', limit=None, norm=None,
- density=1, interpolation='none', grad_encode='dark'):
+ density=1, interpolation='none', grad_encode='bright'):
'''Display the phase map and the resulting color coded holography image in one plot.
-
+
Parameters
----------
- density : float, optional
- The gain factor for determining the number of contour lines. The default is 1.
title : string, optional
The title of the plot. The default is 'Combined Plot'.
+ cmap : string, optional
+ The :class:`~matplotlib.colors.Colormap` which is used for the plot as a string.
+ The default is 'RdBu'.
+ limit : float, optional
+ Plotlimit for the phase in both negative and positive direction (symmetric around 0).
+ If not specified, the maximum amplitude of the phase is used.
+ norm : :class:`~matplotlib.colors.Normalize` or subclass, optional
+ Norm, which is used to determine the colors to encode the phase information.
+ If not specified, :class:`~matplotlib.colors.Normalize` is automatically used.
+ density : float, optional
+ The gain factor for determining the number of contour lines in the holographic
+ contour map. The default is 1.
interpolation : {'none, 'bilinear', 'cubic', 'nearest'}, optional
Defines the interpolation method for the holographic contour map.
No interpolation is used in the default case.
-
+ grad_encode: {'bright', 'dark', 'color', 'none'}, optional
+ Encoding mode of the phase gradient. 'none' produces a black-white image, 'color' just
+ encodes the direction (without gradient strength), 'dark' modulates the gradient
+ strength with a factor between 0 and 1 and 'bright' (which is the default) encodes
+ the gradient strength with color saturation.
+
Returns
-------
phase_axis, holo_axis: :class:`~matplotlib.axes.AxesSubplot`
The axes on which the graphs are plotted.
-
- '''# TODO: Docstring grad_encode!
- self.log.info('Calling display_combined')
+
+ '''
+ self.LOG.debug('Calling display_combined')
# Create combined plot and set title:
fig = plt.figure(figsize=(16, 7))
fig.suptitle(title, fontsize=20)
@@ -507,17 +531,17 @@ class PhaseMap(object):
@classmethod
def make_color_wheel(cls):
'''Display a color wheel to illustrate the color coding of the gradient direction.
-
+
Parameters
----------
None
-
+
Returns
-------
None
-
+
'''
- cls.log.info('Calling make_color_wheel')
+ 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)
diff --git a/pyramid/phasemapper.py b/pyramid/phasemapper.py
index 7dd5de8..70c76e4 100644
--- a/pyramid/phasemapper.py
+++ b/pyramid/phasemapper.py
@@ -1,7 +1,5 @@
# -*- coding: utf-8 -*-
-"""Create magnetic and electric phase maps from magnetization data.
-
-This module executes several forward models to calculate the magnetic or electric phase map from
+"""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.
The electrostatic contribution is calculated by using the assumption of a mean inner potential.
@@ -21,92 +19,145 @@ from pyramid.magdata import MagData
from pyramid.phasemap import PhaseMap
from pyramid.forwardmodel import ForwardModel
+import logging
+class PhaseMapper(object):
+ '''Abstract base class for the phase calculation from a 2-dimensional distribution.
-class PhaseMapper(object):
+ The :class:`~.PhaseeMapper-` class represents a strategy for the phasemapping of a
+ 2-dimensional magnetic distribution with two components onto a scalar phase map.
+ :class:`~.Kernel` is an abstract base class and provides a unified interface which should be
+ subclassed with custom :func:`__init__` and :func:`__call__` functions. Concrete subclasses
+ can be called as a function and take a :class:`~.MagData` object as input and return a
+ :class:`~.PhaseMap` object.
- '''DOCSTRING''' # TODO: Docstring!
+ '''
__metaclass__ = abc.ABCMeta
+ LOG = logging.getLogger(__name__+'.PhaseMapper')
@abc.abstractmethod
def __init__(self, projector):
- '''Docstring: has to be implemented!''' # TODO: Docstring!
-
+ self.LOG.debug('Calling __init__')
+ raise NotImplementedError()
+
@abc.abstractmethod
def __call__(self, mag_data):
- '''Docstring: has to be implemented!''' # TODO: Docstring!
+ self.LOG.debug('Calling __call__')
+ raise NotImplementedError()
class PMAdapterFM(PhaseMapper):
+ '''Class representing a phase mapping strategy adapting the :class:`~.ForwardModel` class.
+
+ The :class:`~.PMAdapterFM` class is an adapter class to incorporate the forward model from the
+ :module:`~.ForwardModel` module without the need of vector input and output. It directly takes
+ :class:`~.MagData` objects and returns :class:`~.PhaseMap` objects.
+
+ Attributes
+ ----------
+ a : float
+ The grid spacing in nm.
+ projector : :class:`~.Projector`
+ Projector which should be used for the projection of the 3-dimensional magnetization
+ distribution.
+ b_0 : float, optional
+ The magnetic induction corresponding to a magnetization `M`\ :sub:`0` in T.
+ The default is 1.
+ geometry : {'disc', 'slab'}, optional
+ Elementary geometry which is used for the phase contribution of one pixel.
+ Default is 'disc'.
+ fwd_model : :class:`~.ForwardModel`
+ The forward model that is constructed from the given information. It is created internally
+ and does not have to be provided by the user.
+
+ '''
+
+ LOG = logging.getLogger(__name__+'.PMAdapterFM')
+
def __init__(self, a, projector, b_0=1, geometry='disc'):
+ self.LOG.debug('Calling __init__')
assert isinstance(projector, Projector), 'Argument has to be a Projector object!'
self.a = a
self.projector = projector
self.fwd_model = ForwardModel([projector], Kernel(a, projector.dim_uv, b_0, geometry))
+ self.LOG.debug('Created '+str(self))
def __call__(self, mag_data):
+ self.LOG.debug('Calling __call__')
assert isinstance(mag_data, MagData), 'Only MagData objects can be mapped!'
- assert mag_data.a == self.a, 'Grid spacing has to match!'
- # TODO: test if mag_data fits in all aspects
+ assert mag_data.a == self.a, 'Grid spacing has to match!'
+ assert mag_data.dim == self.projector.dim, 'Dimensions do not match!'
phase_map = PhaseMap(self.a, np.zeros(self.projector.dim_uv))
phase_map.phase_vec = self.fwd_model(mag_data.mag_vec)
return phase_map
+ def __repr__(self):
+ self.LOG.debug('Calling __repr__')
+ return '%s(a=%r, projector=%r, b_0=%r, geometry=%r)' % \
+ (self.__class__, self.a, self.projector, self.b_0, self.geometry)
+
+ def __str__(self):
+ self.LOG.debug('Calling __str__')
+ return 'PMAdapterFM(a=%s, projector=%s, b_0=%s, geometry=%s)' % \
+ (self.a, self.projector, self.b_0, self.geometry)
+
class PMFourier(PhaseMapper):
-
+
+ '''Class representing a phase mapping strategy using a discretization in Fourier space.
+
+ The :class:`~.PMFourier` class represents a phase mapping strategy involving discretization in
+ Fourier space. It utilizes the Fourier transforms, which are inherently calculated in the
+ :class:`~.Kernel` class and directly takes :class:`~.MagData` objects and returns
+ :class:`~.PhaseMap` objects.
+
+ Attributes
+ ----------
+ a : float
+ The grid spacing in nm.
+ projector : :class:`~.Projector`
+ Projector which should be used for the projection of the 3-dimensional magnetization
+ distribution.
+ b_0 : float, optional
+ The magnetic induction corresponding to a magnetization `M`\ :sub:`0` in T.
+ The default is 1.
+ padding : int, optional
+ Factor for the zero padding. The default is 0 (no padding). For a factor of n the number
+ of pixels is increase by ``(1+n)**2``. Padded zeros are cropped at the end.
+
+ '''
+
+ LOG = logging.getLogger(__name__+'.PMFourier')
PHI_0 = -2067.83 # magnetic flux in T*nm²
def __init__(self, a, projector, b_0=1, padding=0):
+ self.LOG.debug('Calling __init__')
assert isinstance(projector, Projector), 'Argument has to be a Projector object!'
self.a = a
self.projector = projector
self.b_0 = b_0
- self.padding = padding
+ self.padding = padding
+ self.LOG.debug('Created '+str(self))
def __call__(self, mag_data):
- '''Calculate the magnetic phase from magnetization data (Fourier space approach).
-
- Parameters
- ----------
- a : float
- The grid spacing in nm.
- projection : tuple (N=3) of :class:`~numpy.ndarray` (N=2)
- The in-plane projection of the magnetization as a tuple, storing the `u`- and `v`-component
- of the magnetization and the thickness projection for the resulting 2D-grid.
- padding : int, optional
- Factor for the zero padding. The default is 0 (no padding). For a factor of n the number
- of pixels is increase by ``(1+n)**2``. Padded zeros are cropped at the end.
- b_0 : float, optional
- The magnetic induction corresponding to a magnetization `M`\ :sub:`0` in T.
- The default is 1.
-
- Returns
- -------
- phase : :class:`~numpy.ndarray` (N=2)
- The phase as a 2-dimensional array.
-
- ''' # TODO: Docstring
+ self.LOG.debug('Calling __call__')
assert isinstance(mag_data, MagData), 'Only MagData objects can be mapped!'
- assert mag_data.a == self.a, 'Grid spacing has to match!'
- # TODO: test if mag_data fits in all aspects (especially with projector)
+ assert mag_data.a == self.a, 'Grid spacing has to match!'
+ assert mag_data.dim == self.projector.dim, 'Dimensions do not match!'
v_dim, u_dim = self.projector.dim_uv
u_mag, v_mag = self.projector(mag_data.mag_vec).reshape((2,)+self.projector.dim_uv)
# Create zero padded matrices:
- u_pad = u_dim/2 * self.padding
- v_pad = v_dim/2 * self.padding
- # TODO: use mag_data.padding or np.pad(...)
- u_mag_big = np.zeros(((1 + self.padding) * v_dim, (1 + self.padding) * u_dim))
- v_mag_big = np.zeros(((1 + self.padding) * v_dim, (1 + self.padding) * u_dim))
- u_mag_big[v_pad:v_pad+v_dim, u_pad:u_pad+u_dim] = u_mag
- v_mag_big[v_pad:v_pad+v_dim, u_pad:u_pad+u_dim] = v_mag
+ u_pad = int(u_dim/2 * self.padding)
+ v_pad = int(v_dim/2 * self.padding)
+ u_mag_pad = np.pad(u_mag, ((v_pad, v_pad), (u_pad, u_pad)), 'constant', constant_values=0)
+ v_mag_pad = np.pad(v_mag, ((v_pad, v_pad), (u_pad, u_pad)), 'constant', constant_values=0)
# Fourier transform of the two components:
- u_mag_fft = np.fft.fftshift(np.fft.rfft2(u_mag_big), axes=0)
- v_mag_fft = np.fft.fftshift(np.fft.rfft2(v_mag_big), axes=0)
+ u_mag_fft = np.fft.fftshift(np.fft.rfft2(u_mag_pad), axes=0)
+ v_mag_fft = np.fft.fftshift(np.fft.rfft2(v_mag_pad), axes=0)
# Calculate the Fourier transform of the phase:
f_nyq = 0.5 / self.a # nyquist frequency
f_u = np.linspace(0, f_nyq, u_mag_fft.shape[1])
@@ -115,66 +166,67 @@ class PMFourier(PhaseMapper):
coeff = (1j*self.b_0) / (2*self.PHI_0)
phase_fft = coeff*self.a * (u_mag_fft*f_vv - v_mag_fft*f_uu) / (f_uu**2 + f_vv**2 + 1e-30)
# Transform to real space and revert padding:
- phase_big = np.fft.irfft2(np.fft.ifftshift(phase_fft, axes=0))
- phase = phase_big[v_pad:v_pad+v_dim, u_pad:u_pad+u_dim]
+ phase_pad = np.fft.irfft2(np.fft.ifftshift(phase_fft, axes=0))
+ phase = phase_pad[v_pad:v_pad+v_dim, u_pad:u_pad+u_dim]
return PhaseMap(self.a, phase)
+ def __repr__(self):
+ self.LOG.debug('Calling __repr__')
+ return '%s(a=%r, projector=%r, b_0=%r, padding=%r)' % \
+ (self.__class__, self.a, self.projector, self.b_0, self.padding)
+
+ def __str__(self):
+ self.LOG.debug('Calling __str__')
+ return 'PMFourier(a=%s, projector=%s, b_0=%s, padding=%s)' % \
+ (self.a, self.projector, self.b_0, self.padding)
class PMElectric(PhaseMapper):
-
+
+ '''Class representing a phase mapping strategy for the electrostatic contribution.
+
+ The :class:`~.PMElectric` class represents a phase mapping strategy for the electrostatic
+ contribution to the electron phase shift which results e.g. from the mean inner potential in
+ certain samples and which is sensitive to properties of the electron microscope. It directly
+ takes :class:`~.MagData` objects and returns :class:`~.PhaseMap` objects.
+
+ Attributes
+ ----------
+ a : float
+ The grid spacing in nm.
+ projector : :class:`~.Projector`
+ Projector which should be used for the projection of the 3-dimensional magnetization
+ distribution.
+ v_0 : float, optional
+ The mean inner potential of the specimen in V. The default is 1.
+ v_acc : float, optional
+ The acceleration voltage of the electron microscope in V. The default is 30000.
+ threshold : float, optional
+ Threshold for the recognition of the specimen location. The default is 0, which means that
+ all voxels with non-zero magnetization will contribute. Should be above noise level.
+
+ '''
+
+ LOG = logging.getLogger(__name__+'.PMElectric')
H_BAR = 6.626E-34 # Planck constant in J*s
M_E = 9.109E-31 # electron mass in kg
Q_E = 1.602E-19 # electron charge in C
C = 2.998E8 # speed of light in m/s
-
+
def __init__(self, a, projector, v_0=1, v_acc=30000, threshold=0):
- '''Calculate the electric phase from magnetization distributions.
-
- Parameters
- ----------
- a : float
- The grid spacing in nm.
- projection : tuple (N=3) of :class:`~numpy.ndarray` (N=2)
- The in-plane projection of the magnetization as a tuple, storing the u- and v-component
- of the magnetization and the thickness projection for the resulting 2D-grid.
- v_0 : float, optional
- The mean inner potential of the specimen in V. The default is 1.
- v_acc : float, optional
- The acceleration voltage of the electron microscope in V. The default is 30000.
-
- Returns
- -------
- phase : :class:`~numpy.ndarray` (N=2)
- The phase as a 2-dimensional array.
-
- ''' # Docstring!
+ self.LOG.debug('Calling __init__')
+ assert isinstance(projector, Projector), 'Argument has to be a Projector object!'
self.a = a
self.projector = projector
self.v_0 = v_0
self.v_acc = v_acc
self.threshold = threshold
+ self.LOG.debug('Created '+str(self))
def __call__(self, mag_data):
- '''Calculate the electric phase from magnetization distributions.
-
- Parameters
- ----------
- a : float
- The grid spacing in nm.
- projection : tuple (N=3) of :class:`~numpy.ndarray` (N=2)
- The in-plane projection of the magnetization as a tuple, storing the u- and v-component
- of the magnetization and the thickness projection for the resulting 2D-grid.
- v_0 : float, optional
- The mean inner potential of the specimen in V. The default is 1.
- v_acc : float, optional
- The acceleration voltage of the electron microscope in V. The default is 30000.
-
- Returns
- -------
- phase : :class:`~numpy.ndarray` (N=2)
- The phase as a 2-dimensional array.
-
- ''' # Docstring!
+ self.LOG.debug('Calling __call__')
+ assert isinstance(mag_data, MagData), 'Only MagData objects can be mapped!'
+ assert mag_data.a == self.a, 'Grid spacing has to match!'
+ assert mag_data.dim == self.projector.dim, 'Dimensions do not match!'
# Coefficient calculation:
H_BAR, M_E, Q_E, C = self.H_BAR, self.M_E, self.Q_E, self.C
v_0, v_acc = self.v_0, self.v_acc
@@ -187,46 +239,59 @@ class PMElectric(PhaseMapper):
phase = v_0 * Ce * projection * self.a*1E-9
return PhaseMap(self.a, phase)
+ def __repr__(self):
+ self.LOG.debug('Calling __repr__')
+ return '%s(a=%r, projector=%r, v_0=%r, v_acc=%r, threshold=%r)' % \
+ (self.__class__, self.a, self.projector, self.v_0, self.v_acc, self.threshold)
+
+ def __str__(self):
+ self.LOG.debug('Calling __str__')
+ return 'PMElectric(a=%s, projector=%s, v_0=%s, v_acc=%s, threshold=%s)' % \
+ (self.a, self.projector, self.v_0, self.v_acc, self.threshold)
class PMConvolve(PhaseMapper):
- # TODO: Docstring!
-
- def __init__(self, a, projector, b_0=1, threshold=0, geometry='disc'):
- '''Calculate the magnetic phase from magnetization data.
-
- Parameters
- ----------
- a : float
- The grid spacing in nm.
- projection : tuple (N=3) of :class:`~numpy.ndarray` (N=2)
- The in-plane projection of the magnetization as a tuple, storing the `u`- and
- `v`-component of the magnetization and the thickness projection for the resulting
- 2D-grid.
- b_0 : float, optional
- The magnetic induction corresponding to a magnetization `M`\ :sub:`0` in T.
- The default is 1.
- kernel : :class:`~pyramid.kernel.Kernel`, optional
- Specifies the kernel for the convolution with the magnetization data. If none is
- specified, one will be created with `disc` as the default geometry.
-
- Returns
- -------
- phase : :class:`~numpy.ndarray` (N=2)
- The phase as a 2-dimensional array.
-
- '''# Docstring!
+ '''Class representing a phase mapping strategy using real space discretization and Fourier
+ space convolution.
+
+ The :class:`~.PMConvolve` class represents a phase mapping strategy involving discretization in
+ real space. It utilizes the convolution in Fourier space, directly takes :class:`~.MagData`
+ objects and returns :class:`~.PhaseMap` objects.
+
+ Attributes
+ ----------
+ a : float
+ The grid spacing in nm.
+ projector : :class:`~.Projector`
+ Projector which should be used for the projection of the 3-dimensional magnetization
+ distribution.
+ b_0 : float, optional
+ The magnetic induction corresponding to a magnetization `M`\ :sub:`0` in T.
+ The default is 1.
+ geometry : {'disc', 'slab'}, optional
+ Elementary geometry which is used for the phase contribution of one pixel.
+ Default is 'disc'.
+
+ '''
+
+ LOG = logging.getLogger(__name__+'.PMConvolve')
+
+ def __init__(self, a, projector, b_0=1, geometry='disc'):
+ self.LOG.debug('Calling __init__')
assert isinstance(projector, Projector), 'Argument has to be a Projector object!'
self.a = a
self.projector = projector
self.b_0 = b_0
- self.threshold = threshold
self.kernel = Kernel(a, projector.dim_uv, b_0, geometry)
+ self.LOG.debug('Created '+str(self))
def __call__(self, mag_data):
- # Docstring!
+ self.LOG.debug('Calling __call__')
+ assert isinstance(mag_data, MagData), 'Only MagData objects can be mapped!'
+ assert mag_data.a == self.a, 'Grid spacing has to match!'
+ assert mag_data.dim == self.projector.dim, 'Dimensions do not match!'
# Process input parameters:
- u_mag, v_mag = self.projector(mag_data.mag_vec).reshape((2,)+self.projector.dim_uv)
+ u_mag, v_mag = self.projector(mag_data.mag_vec).reshape((2,)+self.projector.dim_uv)
# Fourier transform the projected magnetisation:
kernel = self.kernel
u_mag_fft = np.fft.rfftn(u_mag, kernel.dim_fft)
@@ -237,35 +302,52 @@ class PMConvolve(PhaseMapper):
# Return the result:
return PhaseMap(self.a, u_phase-v_phase)
+ def __repr__(self):
+ self.LOG.debug('Calling __repr__')
+ return '%s(a=%r, projector=%r, b_0=%r, geometry=%r)' % \
+ (self.__class__, self.a, self.projector, self.b_0, self.geometry)
+
+ def __str__(self):
+ self.LOG.debug('Calling __str__')
+ return 'PMConvolve(a=%s, projector=%s, b_0=%s, geometry=%s)' % \
+ (self.a, self.projector, self.b_0, self.geometry)
+
class PMReal(PhaseMapper):
+ '''Class representing a phase mapping strategy using real space discretization.
+
+ The :class:`~.PMReal` class represents a phase mapping strategy involving discretization in
+ real space. It directly takes :class:`~.MagData` objects and returns :class:`~.PhaseMap`
+ objects.
+
+ Attributes
+ ----------
+ a : float
+ The grid spacing in nm.
+ projector : :class:`~.Projector`
+ Projector which should be used for the projection of the 3-dimensional magnetization
+ distribution.
+ b_0 : float, optional
+ The magnetic induction corresponding to a magnetization `M`\ :sub:`0` in T.
+ The default is 1.
+ threshold : float, optional
+ Threshold determining for which voxels the phase contribution will be calculated. The
+ default is 0, which means that all voxels with non-zero magnetization will contribute.
+ Should be above noise level.
+ geometry : {'disc', 'slab'}, optional
+ Elementary geometry which is used for the phase contribution of one pixel.
+ Default is 'disc'.
+ numcore : boolean, optional
+ Boolean choosing if Cython enhanced routines from the :module:`~.pyramid.numcore` module
+ should be used. Default is True.
+
+ '''
+
+ LOG = logging.getLogger(__name__+'.PMReal')
+
def __init__(self, a, projector, b_0=1, threshold=0, geometry='disc', numcore=True):
- '''Calculate the magnetic phase from magnetization data (pure real space, no FFT-convolution).
-
- Parameters
- ----------
- a : float
- The grid spacing in nm.
- projection : tuple (N=3) of :class:`~numpy.ndarray` (N=2)
- The in-plane projection of the magnetization as a tuple, storing the `u`- and `v`-component
- of the magnetization and the thickness projection for the resulting 2D-grid.
- geometry : {'disc', 'slab'}, optional
- Specifies the elemental geometry to use for the pixel field.
- The default is 'disc', because of the smaller computational overhead.
- b_0 : float, optional
- The magnetic induction corresponding to a magnetization `M`\ :sub:`0` in T.
- The default is 1.
- jacobi : :class:`~numpy.ndarray` (N=2), optional
- Specifies the matrix in which to save the jacobi matrix. The jacobi matrix will not be
- calculated, if no matrix is specified (default), resulting in a faster computation.
-
- Returns
- -------
- phase : :class:`~numpy.ndarray` (N=2)
- The phase as a 2-dimensional array.
-
- ''' # Docstring!
+ self.LOG.debug('Calling __init__')
assert isinstance(projector, Projector), 'Argument has to be a Projector object!'
self.a = a
self.projector = projector
@@ -273,10 +355,14 @@ class PMReal(PhaseMapper):
self.threshold = threshold
self.kernel = Kernel(a, projector.dim_uv, b_0, geometry)
self.numcore = numcore
+ self.LOG.debug('Created '+str(self))
def __call__(self, mag_data):
- # TODO: Docstring
- # Process input parameters:
+ self.LOG.debug('Calling __call__')
+ assert isinstance(mag_data, MagData), 'Only MagData objects can be mapped!'
+ assert mag_data.a == self.a, 'Grid spacing has to match!'
+ assert mag_data.dim == self.projector.dim, 'Dimensions do not match!'
+ # Process input parameters:
dim_uv = self.projector.dim_uv
threshold = self.threshold
u_mag, v_mag = self.projector(mag_data.mag_vec).reshape((2,)+dim_uv)
@@ -301,3 +387,14 @@ class PMReal(PhaseMapper):
phase -= v_mag[j, i] * v_phase
# Return the phase:
return PhaseMap(self.a, phase)
+
+ def __repr__(self):
+ self.LOG.debug('Calling __repr__')
+ return '%s(a=%r, projector=%r, b_0=%r, threshold=%r, geometry=%r, numcore=%r)' % \
+ (self.__class__, self.a, self.projector, self.b_0,
+ self.threshold, self.geometry, self.numcore)
+
+ def __str__(self):
+ self.LOG.debug('Calling __str__')
+ return 'PMReal(a=%s, projector=%s, b_0=%s, threshold=%s, geometry=%s, numcore=%s)' % \
+ (self.a, self.projector, self.b_0, self.threshold, self.geometry, self.numcore)
diff --git a/pyramid/projector.py b/pyramid/projector.py
index fd57f77..452bd87 100644
--- a/pyramid/projector.py
+++ b/pyramid/projector.py
@@ -3,8 +3,6 @@
projections of vector and scalar fields."""
-import logging
-
import numpy as np
from numpy import pi
@@ -14,11 +12,13 @@ import itertools
from scipy.sparse import coo_matrix, csr_matrix
+import logging
+
class Projector(object):
'''Abstract base class representing a projection function.
-
+
The :class:`~.Projector` class represents a projection function for a 3-dimensional
vector- or scalar field onto a 2-dimensional grid. :class:`~.Projector` is an abstract base
class and provides a unified interface which should be subclassed with a custom
@@ -30,6 +30,8 @@ class Projector(object):
Attributes
----------
+ dim : tuple (N=3)
+ Dimensions (z, y, x) of the magnetization distribution.
dim_uv : tuple (N=2)
Dimensions (v, u) of the projected grid.
size_3d : int
@@ -46,24 +48,50 @@ class Projector(object):
'''
__metaclass__ = abc.ABCMeta
+ LOG = logging.getLogger(__name__+'.Projector')
@abc.abstractmethod
- def __init__(self, dim_uv, weight, coeff):
- self.log = logging.getLogger(__name__)
- self.log.info('Calling __init__')
+ def __init__(self, dim, dim_uv, weight, coeff):
+ self.LOG.debug('Calling __init__')
+ self.dim = dim
self.dim_uv = dim_uv
self.weight = weight
self.coeff = coeff
self.size_2d, self.size_3d = weight.shape
- self.log.info('Created '+str(self))
+ self.LOG.debug('Created '+str(self))
+
+ def __repr__(self):
+ self.LOG.debug('Calling __repr__')
+ return '%s(dim=%r, dim_uv=%r, weight=%r, coeff=%r)' % \
+ (self.__class__, self.dim, self.dim_uv, self.weight, self.coeff)
+
+ def __str__(self):
+ self.LOG.debug('Calling __str__')
+ return 'Projector(dim=%s, dim_uv=%s, coeff=%s)' % (self.dim, self.dim_uv, self.coeff)
def __call__(self, vector):
- self.log.info('Calling as function')
-# print 'Projector - __call__:', len(vector)
+ self.LOG.debug('Calling as function')
return self.jac_dot(vector)
+ @abc.abstractmethod
+ def get_info(self):
+ '''Get specific information about the projector as a string.
+
+ Paramters
+ ---------
+ None
+
+ Returns
+ -------
+ info : string
+ Information about the projector as a string, e.g. for the use in plot titles.
+
+ '''
+ raise NotImplementedError()
+
+
def _vector_field_projection(self, vector):
- self.log.info('Calling _vector_field_projection')
+ self.LOG.debug('Calling _vector_field_projection')
size_2d, size_3d = self.size_2d, self.size_3d
result = np.zeros(2*size_2d)
# Go over all possible component projections (z, y, x) to (u, v):
@@ -82,7 +110,7 @@ class Projector(object):
return result
def _vector_field_projection_T(self, vector):
- self.log.info('Calling _vector_field_projection_T')
+ self.LOG.debug('Calling _vector_field_projection_T')
size_2d, size_3d = self.size_2d, self.size_3d
result = np.zeros(3*size_3d)
# Go over all possible component projections (u, v) to (z, y, x):
@@ -101,11 +129,11 @@ class Projector(object):
return result
def _scalar_field_projection(self, vector):
- self.log.info('Calling _scalar_field_projection')
+ self.LOG.debug('Calling _scalar_field_projection')
return np.array(self.weight.dot(vector))
def _scalar_field_projection_T(self, vector):
- self.log.info('Calling _scalar_field_projection_T')
+ self.LOG.debug('Calling _scalar_field_projection_T')
return np.array(self.weight.T.dot(vector))
def jac_dot(self, vector):
@@ -124,91 +152,180 @@ class Projector(object):
always`size_2d`.
'''
- self.log.info('Calling jac_dot')
-# print 'Projector - jac_dot:', len(vector)
+ self.LOG.debug('Calling jac_dot')
if len(vector) == 3*self.size_3d: # mode == 'vector'
- self.log.info('mode == vector')
+ self.LOG.debug('mode == vector')
return self._vector_field_projection(vector)
elif len(vector) == self.size_3d: # mode == 'scalar'
- self.log.info('mode == scalar')
+ self.LOG.debug('mode == scalar')
return self._scalar_field_projection(vector)
else:
raise AssertionError('Vector size has to be suited either for ' \
'vector- or scalar-field-projection!')
def jac_T_dot(self, vector):
- # TODO: Docstring!
- self.log.info('Calling jac_T_dot')
-# print 'Projector - jac_T_dot:', len(vector)
+ '''Multiply a `vector` with the transp. jacobi matrix of this :class:`~.Projector` object.
+
+ Parameters
+ ----------
+ vector : :class:`~numpy.ndarray` (N=1)
+ Vector containing the field which should be projected. Must have the same or 2 times
+ the size of `size_2d` of the projector for scalar and vector projection, respectively.
+
+ Returns
+ -------
+ proj_vector : :class:`~numpy.ndarray` (N=1)
+ Vector containing the multiplication of the input with the transposed jacobi matrix
+ of the :class:`~.Projector` object.
+
+ '''
+ self.LOG.debug('Calling jac_T_dot')
if len(vector) == 2*self.size_2d: # mode == 'vector'
- self.log.info('mode == vector')
+ self.LOG.debug('mode == vector')
return self._vector_field_projection_T(vector)
elif len(vector) == self.size_2d: # mode == 'scalar'
- self.log.info('mode == scalar')
+ self.LOG.debug('mode == scalar')
return self._scalar_field_projection_T(vector)
else:
raise AssertionError('Vector size has to be suited either for ' \
'vector- or scalar-field-projection!')
+class XTiltProjector(Projector):
+
+ '''Class representing a projection function with a tilt around the x-axis.
+
+ The :class:`~.XTiltProjector` class represents a projection function for a 3-dimensional
+ vector- or scalar field onto a 2-dimensional grid, which is a concrete subclass of
+ :class:`~.Projector`.
+
+ Attributes
+ ----------
+ dim : tuple (N=3)
+ Dimensions (z, y, x) of the magnetization distribution.
+ tilt : float
+ Angle in `rad` describing the tilt of the beam direction relative to the x-axis.
+
+ '''
+
+ LOG = logging.getLogger(__name__+'.XTiltProjector')
+
+ def __init__(self, dim, tilt):
+
+ def get_position(p, m, b, size):
+ self.LOG.debug('Calling get_position')
+ y, x = np.array(p)[:, 0]+0.5, np.array(p)[:, 1]+0.5
+ return (y-m*x-b)/np.sqrt(m**2+1) + size/2.
+
+ def get_impact(pos, r, size):
+ self.LOG.debug('Calling get_impact')
+ return [x for x in np.arange(np.floor(pos-r), np.floor(pos+r)+1, dtype=int)
+ if 0 <= x < size]
+
+ def get_weight(delta, rho): # use circles to represent the voxels
+ self.LOG.debug('Calling get_weight')
+ lo, up = delta-rho, delta+rho
+ # Upper boundary:
+ if up >= 1:
+ w_up = 0.5
+ else:
+ w_up = (up*np.sqrt(1-up**2) + np.arctan(up/np.sqrt(1-up**2))) / pi
+ # Lower boundary:
+ if lo <= -1:
+ w_lo = -0.5
+ else:
+ w_lo = (lo*np.sqrt(1-lo**2) + np.arctan(lo/np.sqrt(1-lo**2))) / pi
+ return w_up - w_lo
+
+ self.LOG.debug('Calling __init__')
+ self.tilt = tilt
+ # Set starting variables:
+ # length along projection (proj, z), perpendicular (perp, y) and rotation (rot, x) axis:
+ dim_proj, dim_perp, dim_rot = dim
+ size_2d = dim[1] * dim[2] # x-y-plane
+ size_3d = np.prod(dim)
+ # Creating coordinate list of all voxels:
+ voxels = list(itertools.product(range(dim_proj), range(dim_perp))) # z-y-plane
+ # Calculate positions along the projected pixel coordinate system:
+ center = (dim_proj/2., dim_perp/2.)
+ m = np.where(tilt<=pi, -1/np.tan(tilt+1E-30), 1/np.tan(tilt+1E-30))
+ b = center[0] - m * center[1]
+ positions = get_position(voxels, m, b, dim_perp)
+ # Calculate weight-matrix:
+ r = 1/np.sqrt(np.pi) # radius of the voxel circle
+ rho = 0.5 / r
+ row = []
+ col = []
+ data = []
+ # one slice:
+ for i, voxel in enumerate(voxels):
+ impacts = get_impact(positions[i], r, dim_perp) # impact along y axis
+ for impact in impacts:
+ distance = np.abs(impact+0.5 - positions[i])
+ delta = distance / r
+ col.append(voxel[0]*size_2d + voxel[1]*dim_rot) # 0: z, 1: y
+ row.append(impact*dim_rot)
+ data.append(get_weight(delta, rho))
+ # All other slices (along x):
+ columns = col
+ rows = row
+ for i in np.arange(1, dim_rot):
+ columns = np.hstack((np.array(columns), np.array(col)+i))
+ rows = np.hstack((np.array(rows), np.array(row)+i))
+ # Calculate weight matrix and coefficients for jacobi matrix:
+ weight = csr_matrix(coo_matrix((np.tile(data, dim_rot), (rows, columns)),
+ shape = (size_2d, size_3d)))
+ dim_v, dim_u = dim_perp, dim_rot
+ coeff = [[1, 0, 0], [0, np.cos(tilt), np.sin(tilt)]]
+ super(XTiltProjector, self).__init__(dim, (dim_v, dim_u), weight, coeff)
+ self.LOG.debug('Created '+str(self))
+
+ def get_info(self):
+ '''Get specific information about the projector as a string.
+
+ Paramters
+ ---------
+ None
+
+ Returns
+ -------
+ info : string
+ Information about the projector as a string, e.g. for the use in plot titles.
+
+ '''
+ return 'x-tilt: $(\phi = {:2.1f} \pi)$'.format(self.tilt)
+
class YTiltProjector(Projector):
-# '''Class representing a projection where the projection axis is tilted around the y-axis.
-#
-# The :class:`~.YTiltProjector` class is a concrete subclass of the :class:`~.Projector` class
-# and overwrites the :func:`__init__` constructor which accepts `dim` and `tilt` as arguments.
-# The dimensions of the 3-dimensional grid are given by `dim`, the tilting angle of the ample
-# around the y-axis is given by `tilt`.
-#
-# Attributes
-# ----------
-# dim_uv : tuple (N=2)
-# Dimensions (v, u) of the projected grid.
-# size_3d : int
-# Number of voxels of the 3-dimensional grid.
-# size_2d : int
-# Number of pixels of the 2-dimensional projected grid.
-# weight : :class:`~scipy.sparse.csr_matrix` (N=2)
-# The weight matrix containing the weighting coefficients describing the influence of all
-# 3-dimensional voxels on the 2-dimensional pixels of the projection.
-# coeff : list (N=2)
-# List containing the six weighting coefficients describing the influence of the 3 components
-# of a 3-dimensional vector field on the 2 projected components.
-# '''
-#
-# '''
-# weight : :class:`~scipy.sparse.csr_matrix` (N=2)
-# The weight matrix containing the weighting coefficients which determine the influence of
-# all 3-dimensional voxels to the 2-dimensional pixels of the projection.
-# coeff : `list`t (N=2)
-# List containing the six weighting coefficients describing the influence of the 3 components
-# of a 3-dimensional vector fields on the 2 components of the projected field. Only used for
-# vector field projection.
-#
-# Notes
-# -----
-# An instance `projector` of the :class:`~.YSimpleProjector` class is callable via:
-#
-# :func:`projector(vector)`
-#
-# with `vector` being a :class:`~numpy.ndarray` (N=1).
-# '''
+ '''Class representing a projection function with a tilt around the y-axis.
+
+ The :class:`~.YTiltProjector` class represents a projection function for a 3-dimensional
+ vector- or scalar field onto a 2-dimensional grid, which is a concrete subclass of
+ :class:`~.Projector`.
+
+ Attributes
+ ----------
+ dim : tuple (N=3)
+ Dimensions (z, y, x) of the magnetization distribution.
+ tilt : float
+ Angle in `rad` describing the tilt of the beam direction relative to the y-axis.
+
+ '''
+
+ LOG = logging.getLogger(__name__+'.YTiltProjector')
def __init__(self, dim, tilt):
def get_position(p, m, b, size):
- self.log.info('Calling get_position')
y, x = np.array(p)[:, 0]+0.5, np.array(p)[:, 1]+0.5
return (y-m*x-b)/np.sqrt(m**2+1) + size/2.
-
+
def get_impact(pos, r, size):
- self.log.info('Calling get_impact')
return [x for x in np.arange(np.floor(pos-r), np.floor(pos+r)+1, dtype=int)
if 0 <= x < size]
-
+
def get_weight(delta, rho): # use circles to represent the voxels
- self.log.info('Calling get_weight')
lo, up = delta-rho, delta+rho
# Upper boundary:
if up >= 1:
@@ -222,16 +339,15 @@ class YTiltProjector(Projector):
w_lo = (lo*np.sqrt(1-lo**2) + np.arctan(lo/np.sqrt(1-lo**2))) / pi
return w_up - w_lo
- self.log = logging.getLogger(__name__)
- self.log.info('Calling __init__')
+ self.LOG.debug('Calling __init__')
self.tilt = tilt
# Set starting variables:
# length along projection (proj, z), rotation (rot, y) and perpendicular (perp, x) axis:
dim_proj, dim_rot, dim_perp = dim
- size_2d = dim_rot * dim_perp
- size_3d = dim_proj * dim_rot * dim_perp
+ size_2d = dim[1] * dim[2] # x-y-plane
+ size_3d = np.prod(dim)
# Creating coordinate list of all voxels:
- voxels = list(itertools.product(range(dim_proj), range(dim_perp)))
+ voxels = list(itertools.product(range(dim_proj), range(dim_perp))) # z-x-plane
# Calculate positions along the projected pixel coordinate system:
center = (dim_proj/2., dim_perp/2.)
m = np.where(tilt<=pi, -1/np.tan(tilt+1E-30), 1/np.tan(tilt+1E-30))
@@ -245,17 +361,17 @@ class YTiltProjector(Projector):
data = []
# one slice:
for i, voxel in enumerate(voxels):
- impacts = get_impact(positions[i], r, dim_perp)
+ impacts = get_impact(positions[i], r, dim_perp) # impact along x axis
for impact in impacts:
distance = np.abs(impact+0.5 - positions[i])
delta = distance / r
- col.append(voxel[0]*size_2d + voxel[1])
+ col.append(voxel[0]*size_2d + voxel[1]) # 0: z, 1: x
row.append(impact)
data.append(get_weight(delta, rho))
- # All other slices:
+ # All other slices (along y):
columns = col
rows = row
- for i in np.arange(1, dim_rot): # TODO: more efficient, please!
+ for i in np.arange(1, dim_rot):
columns = np.hstack((np.array(columns), np.array(col)+i*dim_perp))
rows = np.hstack((np.array(rows), np.array(row)+i*dim_perp))
# Calculate weight matrix and coefficients for jacobi matrix:
@@ -263,50 +379,49 @@ class YTiltProjector(Projector):
shape = (size_2d, size_3d)))
dim_v, dim_u = dim_rot, dim_perp
coeff = [[np.cos(tilt), 0, np.sin(tilt)], [0, 1, 0]]
- super(YTiltProjector, self).__init__((dim_v, dim_u), weight, coeff)
- self.log.info('Created '+str(self))
+ super(YTiltProjector, self).__init__(dim, (dim_v, dim_u), weight, coeff)
+ self.LOG.debug('Created '+str(self))
+
+ def get_info(self):
+ '''Get specific information about the projector as a string.
+
+ Paramters
+ ---------
+ None
+
+ Returns
+ -------
+ info : string
+ Information about the projector as a string, e.g. for the use in plot titles.
+
+ '''
+ return 'y-tilt: $(\phi = {:2.1f} \pi)$'.format(self.tilt)
class SimpleProjector(Projector):
-# '''Class representing a projection along one of the major axes (x, y, z).
-#
-# The :class:`~.SimpleProjector` class is a concrete subclass of the :class:`~.Projector` class
-# and overwrites the :func:`__init__` constructor which accepts `dim` and `axis` as arguments.
-# The dimensions of the 3-dimensional grid are given by `dim`, the major axis along which to
-# project is given by `axis` and can be `'x'`, `'y'` or `'z'` (default).
-#
-# Attributes
-# ----------
-# dim_uv : tuple (N=2)
-# Dimensions (v, u) of the projected grid.
-# weight : :class:`~scipy.sparse.csr_matrix` (N=2)
-# The weight matrix containing the weighting coefficients which determine the influence of
-# all 3-dimensional voxels to the 2-dimensional pixels of the projection.
-# coeff : list (N=2)
-# List containing the six weighting coefficients describing the influence of the 3 components
-# of a 3-dimensional vector fields on the 2 components of the projected field. Only used for
-# vector field projection.
-# size_3d : int
-# Number of voxels of the 3-dimensional grid.
-# size_2d : int
-# Number of pixels of the 2-dimensional projected grid.
-#
-# Notes
-# -----
-# An instance `projector` of the :class:`~.SimpleProjector` class is callable via:
-#
-# :func:`projector(vector)`
-#
-# with `vector` being a :class:`~numpy.ndarray` (N=1).
-#
-# '''
+ '''Class representing a projection function along one of the major axes.
+
+ The :class:`~.SimpleProjector` class represents a projection function for a 3-dimensional
+ vector- or scalar field onto a 2-dimensional grid, which is a concrete subclass of
+ :class:`~.Projector`.
+
+ Attributes
+ ----------
+ dim : tuple (N=3)
+ Dimensions (z, y, x) of the magnetization distribution.
+ axis : {'z', 'y', 'x'}, optional
+ Main axis along which the magnetic distribution is projected (given as a string). Defaults
+ to the z-axis.
+
+ '''
+ LOG = logging.getLogger(__name__+'.SimpleProjector')
AXIS_DICT = {'z': (0, 1, 2), 'y': (1, 0, 2), 'x': (1, 2, 0)}
def __init__(self, dim, axis='z'):
- self.log = logging.getLogger(__name__)
- self.log.info('Calling __init__')
+ 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]
dim_proj, dim_v, dim_u = dim[proj], dim[v], dim[u]
dim_z, dim_y, dim_x = dim
@@ -315,17 +430,35 @@ class SimpleProjector(Projector):
data = np.repeat(1, size_3d)
indptr = np.arange(0, size_3d+1, dim_proj)
if axis == 'z':
+ self.LOG.debug('Projecting along the z-axis')
coeff = [[1, 0, 0], [0, 1, 0]]
- indices = np.array([np.arange(row, size_3d, size_2d)
+ indices = np.array([np.arange(row, size_3d, size_2d)
for row in range(size_2d)]).reshape(-1)
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
for row in range(size_2d)]).reshape(-1)
elif axis == 'x':
+ self.LOG.debug('Projection along the x-axis')
coeff = [[0, 1, 0], [0, 0, 1]]
indices = np.array([np.arange(dim_proj) + row*dim_proj
for row in range(size_2d)]).reshape(-1)
weight = csr_matrix((data, indices, indptr), shape = (size_2d, size_3d))
- super(SimpleProjector, self).__init__((dim_v, dim_u), weight, coeff)
- self.log.info('Created '+str(self))
+ super(SimpleProjector, self).__init__(dim, (dim_v, dim_u), weight, coeff)
+ self.LOG.debug('Created '+str(self))
+
+ def get_info(self):
+ '''Get specific information about the projector as a string.
+
+ Paramters
+ ---------
+ None
+
+ Returns
+ -------
+ info : string
+ Information about the projector as a string, e.g. for the use in plot titles.
+
+ '''
+ return 'projected along {}-axis'.format(self.tilt)
diff --git a/pyramid/reconstruction.py b/pyramid/reconstruction.py
new file mode 100644
index 0000000..fb8a20f
--- /dev/null
+++ b/pyramid/reconstruction.py
@@ -0,0 +1,117 @@
+# -*- coding: utf-8 -*-
+"""Reconstruct magnetic distributions from given phasemaps.
+
+This module reconstructs 3-dimensional magnetic distributions (as :class:`~pyramid.magdata.MagData`
+objects) from a given set of phase maps (represented by :class:`~pyramid.phasemap.PhaseMap`
+objects) by using several model based reconstruction algorithms which use the forward model
+provided by :mod:`~pyramid.projector` and :mod:`~pyramid.phasemapper` and a priori knowledge of
+the distribution.
+
+"""
+
+
+import numpy as np
+
+from scipy.sparse.linalg import cg
+
+from pyramid.kernel import Kernel
+from pyramid.forwardmodel import ForwardModel
+from pyramid.costfunction import Costfunction, CFAdapterScipyCG
+from pyramid.magdata import MagData
+
+import logging
+
+
+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 : {2, 1, 0}, optional
+ Parameter defining the verposity of the output. `2` is the default and will show the
+ current number of the iteration and the cost of the current distribution. `2` 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=2):
+ 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_sparse_cg(data, verbosity=2):
+ '''Reconstruct a three-dimensional magnetic distribution from given phase maps via the
+ conjugate gradient optimizaion method :func:`~.scipy.sparse.linalg.cg`.
+
+ Parameters
+ ----------
+ data : :class:`~.DataSet`
+ :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.
+ verbosity : {2, 1, 0}, optional
+ Parameter defining the verposity of the output. `2` is the default and will show the
+ current number of the iteration and the cost of the current distribution. `2` will just
+ show the iteration number and `0` will prevent output all together.
+
+ Returns
+ -------
+ mag_data : :class:`~pyramid.magdata.MagData`
+ The reconstructed magnetic distribution as a :class:`~.MagData` object.
+
+ '''
+ LOG.debug('Calling optimize_sparse_cg')
+ # Set up necessary objects:
+ y = data.phase_vec
+ kernel = Kernel(data.a, data.dim_uv, data.b_0)
+ fwd_model = ForwardModel(data.projectors, kernel)
+ cost = Costfunction(y, fwd_model, lam=10**-10)
+ # Optimize:
+ A = CFAdapterScipyCG(cost)
+ b = fwd_model.jac_T_dot(None, y)
+ x_opt, info = cg(A, b, callback=PrintIterator(cost, verbosity))
+ # Create and return fitting MagData object:
+ mag_opt = MagData(fwd_model.a, np.zeros((3,)+fwd_model.dim))
+ mag_opt.mag_vec = x_opt
+ return mag_opt
diff --git a/scripts/abstract_pictures.py b/scripts/abstract_pictures.py
new file mode 100644
index 0000000..cb20acd
--- /dev/null
+++ b/scripts/abstract_pictures.py
@@ -0,0 +1,115 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Mon Mar 17 14:28:10 2014
+
+@author: Jan
+"""
+
+
+from numpy import pi
+import numpy as np
+import matplotlib.pyplot as plt
+from matplotlib.ticker import FuncFormatter
+
+import pyramid.magcreator as mc
+from pyramid.magdata import MagData
+from pyramid.projector import SimpleProjector, YTiltProjector
+from pyramid.phasemapper import PMConvolve
+
+
+###################################################################################################
+print('Jan')
+
+dim = (128, 128, 128)
+center = (int(dim[0]/2), int(dim[1]/2), int(dim[2]/2))
+radius = dim[0]/4.
+a = 1.
+
+magnitude = mc.create_mag_dist_homog(mc.Shapes.sphere(dim, center, radius), pi/4)
+
+mag_data = MagData(a, magnitude)
+
+projector = SimpleProjector(dim)
+
+phase_map = PMConvolve(a, projector)(mag_data)
+
+axis = phase_map.display_phase()
+axis.tick_params(axis='both', which='major', labelsize=18)
+axis.set_title('Phase map', fontsize=24)
+axis.set_xlabel('x [nm]', fontsize=18)
+axis.set_ylabel('y [nm]', fontsize=18)
+
+axis = phase_map.display_holo(density=20, interpolation='bilinear')
+axis.tick_params(axis='both', which='major', labelsize=18)
+axis.set_title('Magnetic induction map', fontsize=24)
+axis.set_xlabel('x [nm]', fontsize=18)
+axis.set_ylabel('y [nm]', fontsize=18)
+
+mag_data.scale_down(2)
+
+mag_data.quiver_plot()
+axis = plt.gca()
+axis.tick_params(axis='both', which='major', labelsize=18)
+axis.set_title('Magnetization distribution', fontsize=24)
+axis.set_xlabel('x [nm]', fontsize=18)
+axis.set_ylabel('y [nm]', fontsize=18)
+
+phase_map.make_color_wheel()
+
+shape_vort = mc.Shapes.disc((64, 64, 64), (31.5, 31.5, 31.5), 24, 10)
+
+magnitude_vort = mc.create_mag_dist_vortex(shape_vort)
+
+mag_vort = MagData(a, magnitude_vort)
+
+mag_vort.scale_down(2)
+
+mag_vort.quiver_plot()
+
+
+###################################################################################################
+print('Patrick')
+
+a = 10.0 # nm
+b_0 = 3 # T
+
+dim = (128, 128, 128)
+center = (int(dim[0]/2), int(dim[1]/2), int(dim[2]/2)) # in px (z, y, x), index starts with 0!
+
+mag_data = MagData(a, mc.create_mag_dist_homog(mc.Shapes.ellipse(dim, center, (20., 60.), 5.), 0))
+
+tilts = np.array([0., 60.])/180.*pi
+
+projectors = [YTiltProjector(mag_data.dim, tilt) for tilt in tilts]
+phasemappers = [PMConvolve(mag_data.a, proj, b_0) for proj in projectors]
+
+phase_maps = [pm(mag_data) for pm in phasemappers]
+
+axis = phase_maps[0].display_holo(density=1, interpolation='bilinear')
+axis.tick_params(axis='both', which='major', labelsize=18)
+axis.set_title('Magnetic induction map', fontsize=24)
+axis.xaxis.set_major_formatter(FuncFormatter(lambda x, pos: '{:g}'.format(x*a)))
+axis.yaxis.set_major_formatter(FuncFormatter(lambda x, pos: '{:g}'.format(x*a)))
+axis.set_xlabel('x [nm]', fontsize=18)
+axis.set_ylabel('y [nm]', fontsize=18)
+
+axis = phase_maps[0].display_phase(limit=17)
+axis.tick_params(axis='both', which='major', labelsize=18)
+axis.set_title('Phase map', fontsize=24)
+axis.set_xlabel('x [nm]', fontsize=18)
+axis.set_ylabel('y [nm]', fontsize=18)
+
+axis = phase_maps[1].display_holo(density=1, interpolation='bilinear')
+axis.tick_params(axis='both', which='major', labelsize=18)
+axis.set_title('Magnetic induction map', fontsize=24)
+axis.xaxis.set_major_formatter(FuncFormatter(lambda x, pos: '{:g}'.format(x*a)))
+axis.yaxis.set_major_formatter(FuncFormatter(lambda x, pos: '{:g}'.format(x*a)))
+axis.set_xlabel('x [nm]', fontsize=18)
+axis.set_ylabel('y [nm]', fontsize=18)
+
+axis = phase_maps[1].display_phase(limit=17)
+axis.tick_params(axis='both', which='major', labelsize=18)
+axis.set_title('Phase map', fontsize=24)
+axis.set_xlabel('x [nm]', fontsize=18)
+axis.set_ylabel('y [nm]', fontsize=18)
+
diff --git a/scripts/compare methods/logfile.log b/scripts/compare methods/logfile.log
deleted file mode 100644
index e69de29..0000000
diff --git a/scripts/create distributions/logfile.log b/scripts/create distributions/logfile.log
deleted file mode 100644
index e69de29..0000000
diff --git a/scripts/edinburgh_test.py b/scripts/edinburgh_test.py
index 27bf0cc..42fcf27 100644
--- a/scripts/edinburgh_test.py
+++ b/scripts/edinburgh_test.py
@@ -9,7 +9,7 @@ Created on Tue Jan 28 15:15:08 2014
import numpy as np
from pyramid.magdata import MagData
from pyramid.projector import SimpleProjector
-from pyramid.phasemapper import PMAdapterFM
+from pyramid.phasemapper import PMAdapterFM, PMFourier
from matplotlib.ticker import FuncFormatter
data = np.loadtxt('../output/data from Edinburgh/long_grain_remapped_0p0035.txt', delimiter=',')
@@ -28,17 +28,18 @@ magnitude = np.array((x_mag, y_mag, z_mag))
mag_data = MagData(a, magnitude)
-mag_data.pad(30, 20, 0)
-
-mag_data.scale_up()
-
-mag_data.quiver_plot()
+#mag_data.pad(30, 20, 0)
+#
+#mag_data.scale_up()
+#
+#mag_data.quiver_plot()
#mag_data.quiver_plot3d()
projector = SimpleProjector(mag_data.dim)
phasemapper = PMAdapterFM(mag_data.a, projector)
+phasemapper = PMFourier(mag_data.a, projector, padding = 1)
phase_map = phasemapper(mag_data)
diff --git a/scripts/interactive_setup.py b/scripts/interactive_setup.py
index 888286d..804a13f 100644
--- a/scripts/interactive_setup.py
+++ b/scripts/interactive_setup.py
@@ -1,14 +1,21 @@
# -*- coding: utf-8 -*-
"""Imports pyramid package and sets working directory to output directory"""
-import pyramid.analytic as an
-import pyramid.holoimage as hi
+
import pyramid.magcreator as mc
-import pyramid.phasemapper as pm
-import pyramid.projector as pj
-import pyramid.reconstructor as rc
-from pyramid.magdata import MagData
-from pyramid.phasemap import PhaseMap
+import pyramid.analytic as an
+import pyramid.reconstruction as rc
+
+from pyramid.magdata import *
+from pyramid.projector import *
+from pyramid.kernel import *
+from pyramid.phasemapper import *
+from pyramid.phasemap import *
+from pyramid.dataset import *
+from pyramid.forwardmodel import *
+from pyramid.costfunction import *
import os
-os.chdir('C:\\Users\\Jan\\Home\\PhD Thesis\\Pyramid\\tests')
+
+
+os.chdir('C:\\Users\\Jan\\Home\\PhD Thesis\\Pyramid\\output')
diff --git a/scripts/logfile.log b/scripts/logfile.log
deleted file mode 100644
index e69de29..0000000
diff --git a/scripts/paper 1/logfile.log b/scripts/paper 1/logfile.log
deleted file mode 100644
index 3bc36a6..0000000
--- a/scripts/paper 1/logfile.log
+++ /dev/null
@@ -1,84 +0,0 @@
-2014-02-18 09:30:30: INFO @ <pyramid.projector>: Calling __init__
-2014-02-18 09:30:30: INFO @ <pyramid.projector>: Calling __init__
-2014-02-18 09:30:30: INFO @ <pyramid.projector>: Created <pyramid.projector.SimpleProjector object at 0x0DC7ECD0>
-2014-02-18 09:30:30: INFO @ <pyramid.projector>: Created <pyramid.projector.SimpleProjector object at 0x0DC7ECD0>
-2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Calling __str__
-2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim_uv=(128, 128))
-2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Calling __str__
-2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim_uv=(128, 128))
-2014-02-18 09:30:30: INFO @ <pyramid.projector>: Calling as function
-2014-02-18 09:30:30: INFO @ <pyramid.projector>: mode == vector
-2014-02-18 09:30:30: INFO @ <pyramid.projector>: Calling _vector_field_projection
-2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Calling __str__
-2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim_uv=(128, 128))
-2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Calling __sub__
-2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Calling __neg__
-2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Calling __str__
-2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim_uv=(128, 128))
-2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Calling __add__
-2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Adding two PhaseMap objects
-2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Calling __str__
-2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim_uv=(128, 128))
-2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Calling display_phase
-2014-02-18 09:30:30: INFO @ <pyramid.projector>: Calling as function
-2014-02-18 09:30:30: INFO @ <pyramid.projector>: mode == vector
-2014-02-18 09:30:30: INFO @ <pyramid.projector>: Calling _vector_field_projection
-2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Calling __str__
-2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim_uv=(128, 128))
-2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Calling __sub__
-2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Calling __neg__
-2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Calling __str__
-2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim_uv=(128, 128))
-2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Calling __add__
-2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Adding two PhaseMap objects
-2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Calling __str__
-2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim_uv=(128, 128))
-2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Calling __isub__
-2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Calling __sub__
-2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Calling __add__
-2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Adding an offset
-2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Calling __str__
-2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim_uv=(128, 128))
-2014-02-18 09:30:30: INFO @ <pyramid.phasemap>: Calling display_phase
-2014-02-18 09:30:31: INFO @ <pyramid.projector>: Calling as function
-2014-02-18 09:30:31: INFO @ <pyramid.projector>: mode == vector
-2014-02-18 09:30:31: INFO @ <pyramid.projector>: Calling _vector_field_projection
-2014-02-18 09:30:31: INFO @ <pyramid.phasemap>: Calling __str__
-2014-02-18 09:30:31: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim_uv=(128, 128))
-2014-02-18 09:30:31: INFO @ <pyramid.phasemap>: Calling __sub__
-2014-02-18 09:30:31: INFO @ <pyramid.phasemap>: Calling __neg__
-2014-02-18 09:30:31: INFO @ <pyramid.phasemap>: Calling __str__
-2014-02-18 09:30:31: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim_uv=(128, 128))
-2014-02-18 09:30:31: INFO @ <pyramid.phasemap>: Calling __add__
-2014-02-18 09:30:31: INFO @ <pyramid.phasemap>: Adding two PhaseMap objects
-2014-02-18 09:30:31: INFO @ <pyramid.phasemap>: Calling __str__
-2014-02-18 09:30:31: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim_uv=(128, 128))
-2014-02-18 09:30:31: INFO @ <pyramid.phasemap>: Calling display_phase
-2014-02-18 09:30:31: INFO @ <pyramid.projector>: Calling as function
-2014-02-18 09:30:31: INFO @ <pyramid.projector>: mode == vector
-2014-02-18 09:30:31: INFO @ <pyramid.projector>: Calling _vector_field_projection
-2014-02-18 09:30:32: INFO @ <pyramid.phasemap>: Calling __str__
-2014-02-18 09:30:32: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim_uv=(128, 128))
-2014-02-18 09:30:32: INFO @ <pyramid.phasemap>: Calling __sub__
-2014-02-18 09:30:32: INFO @ <pyramid.phasemap>: Calling __neg__
-2014-02-18 09:30:32: INFO @ <pyramid.phasemap>: Calling __str__
-2014-02-18 09:30:32: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim_uv=(128, 128))
-2014-02-18 09:30:32: INFO @ <pyramid.phasemap>: Calling __add__
-2014-02-18 09:30:32: INFO @ <pyramid.phasemap>: Adding two PhaseMap objects
-2014-02-18 09:30:32: INFO @ <pyramid.phasemap>: Calling __str__
-2014-02-18 09:30:32: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim_uv=(128, 128))
-2014-02-18 09:30:32: INFO @ <pyramid.phasemap>: Calling __isub__
-2014-02-18 09:30:32: INFO @ <pyramid.phasemap>: Calling __sub__
-2014-02-18 09:30:32: INFO @ <pyramid.phasemap>: Calling __add__
-2014-02-18 09:30:32: INFO @ <pyramid.phasemap>: Adding an offset
-2014-02-18 09:30:32: INFO @ <pyramid.phasemap>: Calling __str__
-2014-02-18 09:30:32: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim_uv=(128, 128))
-2014-02-18 09:30:32: INFO @ <pyramid.phasemap>: Calling display_phase
-2014-02-18 09:30:33: INFO @ <pyramid.projector>: Calling __init__
-2014-02-18 09:30:33: INFO @ <pyramid.projector>: Calling __init__
-2014-02-18 09:30:33: INFO @ <pyramid.projector>: Created <pyramid.projector.SimpleProjector object at 0x0DC7EB70>
-2014-02-18 09:30:33: INFO @ <pyramid.projector>: Created <pyramid.projector.SimpleProjector object at 0x0DC7EB70>
-2014-02-18 09:30:33: INFO @ <pyramid.phasemap>: Calling __str__
-2014-02-18 09:30:33: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim_uv=(128, 128))
-2014-02-18 09:30:33: INFO @ <pyramid.phasemap>: Calling __str__
-2014-02-18 09:30:33: INFO @ <pyramid.phasemap>: Created PhaseMap(a=1.0, dim_uv=(128, 128))
diff --git a/scripts/simple_reconstruction.py b/scripts/simple_reconstruction.py
index f487272..45385f9 100644
--- a/scripts/simple_reconstruction.py
+++ b/scripts/simple_reconstruction.py
@@ -10,15 +10,17 @@ import numpy as np
from numpy import pi
from pyramid.magdata import MagData
-from pyramid.projector import YTiltProjector
+from pyramid.projector import YTiltProjector, XTiltProjector
from pyramid.phasemapper import PMConvolve
-from pyramid.datacollection import DataCollection
+from pyramid.dataset import DataSet
+import pyramid.magcreator as mc
from pyramid.kernel import Kernel
from pyramid.forwardmodel import ForwardModel
from pyramid.costfunction import Costfunction
+import pyramid.reconstruction as rc
-from scipy.sparse.linalg import cg
+from scipy.sparse.linalg import cg, LinearOperator
###################################################################################################
@@ -41,7 +43,7 @@ phase_maps = [pm(mag_data) for pm in phasemappers]
dim_uv = dim[1:3]
-data_collection = DataCollection(a, dim_uv, b_0)
+data_collection = DataSet(a, dim_uv, b_0)
[data_collection.append((phase_maps[i], projectors[i])) for i in range(count)]
@@ -100,7 +102,7 @@ phase_maps = [pm(mag_data) for pm in phasemappers]
dim_uv = dim[1:3]
-data_collection = DataCollection(a, dim_uv, b_0)
+data_collection = DataSet(a, dim_uv, b_0)
[data_collection.append((phase_maps[i], projectors[i])) for i in range(count)]
@@ -166,35 +168,40 @@ print('--STARTING RECONSTRUCTION')
###################################################################################################
print('--Generating input phase_maps')
-a = 1.
+a = 10.
b_0 = 1000.
-dim = (9, 9, 9)
-count = 8
+dim = (8, 8, 8)
+count = 32
magnitude = np.zeros((3,)+dim)
-magnitude[0, int(dim[0]/2), int(dim[1]/2), int(dim[2]/2)] = 1.
+magnitude[:, 3:6, 3:6, 3:6] = 1# int(dim[0]/2), int(dim[1]/2), int(dim[2]/2)] = 1.
+magnitude = mc.create_mag_dist_vortex(mc.Shapes.disc(dim, (3.5, 3.5, 3.5), 3, 4))
mag_data = MagData(a, magnitude)
mag_data.quiver_plot3d()
-tilts = np.linspace(0, 2*pi, num=count, endpoint=False)
-projectors = [YTiltProjector(mag_data.dim, tilt) for tilt in tilts]
+tilts = np.linspace(0, 2*pi, num=count/2, endpoint=False)
+projectors = []
+projectors.extend([XTiltProjector(mag_data.dim, tilt) for tilt in tilts])
+projectors.extend([YTiltProjector(mag_data.dim, tilt) for tilt in tilts])
phasemappers = [PMConvolve(mag_data.a, projector, b_0) for projector in projectors]
phase_maps = [pm(mag_data) for pm in phasemappers]
-[phase_map.display_phase(title=u'Tilt series $(\phi = {:2.1f} \pi)$'.format(tilts[i]/pi))
- for i, phase_map in enumerate(phase_maps)]
+#[phase_map.display_phase(title=u'Tilt series $(\phi = {:2.1f} \pi)$'.format(tilts[i%(count/2)]/pi))
+# for i, phase_map in enumerate(phase_maps)]
###################################################################################################
print('--Setting up data collection')
dim_uv = dim[1:3]
+lam = 10. ** -10
+
size_2d = np.prod(dim_uv)
size_3d = np.prod(dim)
-data_collection = DataCollection(a, dim_uv, b_0)
+data_collection = DataSet(a, dim_uv, b_0)
[data_collection.append((phase_maps[i], projectors[i])) for i in range(count)]
@@ -202,37 +209,62 @@ data = data_collection
y = data.phase_vec
kern = Kernel(data.a, data.dim_uv, data.b_0)
F = ForwardModel(data.projectors, kern)
-C = Costfunction(y, F)
+C = Costfunction(y, F, lam)
###################################################################################################
print('--Test simple solver')
-M = np.asmatrix([F.jac_dot(None, np.eye(3*size_3d)[:, i]) for i in range(3*size_3d)]).T
-lam = 1*10. ** -10
-MTM = M.T * M + lam * np.asmatrix(np.eye(3*size_3d))
-
-A = MTM#np.array([F(np.eye(81)[:, i]) for i in range(81)])
-
-b = F.jac_T_dot(None, y)
-
-b_test = np.asarray((M.T.dot(y)).T)[0]
+#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))
+#A = MTM#np.array([F(np.eye(81)[:, i]) for i in range(81)])
-x_f = cg(A, b)[0]
-
-mag_data_rec = MagData(a, np.zeros((3,)+dim))
-
-mag_data_rec.mag_vec = x_f
-
-mag_data_rec.quiver_plot3d()
+#class A_adapt(LinearOperator):
+#
+# def __init__(self, FwdModel, lam, shape):
+# self.fwd = FwdModel
+# self.lam = lam
+# self.shape = shape
+#
+# def matvec(self, vector):
+# return self.fwd.jac_T_dot(None, self.fwd.jac_dot(None, vector)) + self.lam*vector
+#
+# @property
+# def shape(self):
+# return self.shape
+#
+# @property
+# def dtype(self):
+# return np.dtype("d") # None #np.ones(1).dtype
+#
+## TODO: .shape in F und C
+#
+#b = F.jac_T_dot(None, y)
+#
+#A_fast = A_adapt(F, lam, (3*size_3d, 3*size_3d))
+#
+#i = 0
+#def printit(_):
+# global i
+# i += 1
+# print i
+#
+#x_f, info = cg(A_fast, b, callback=printit)
+#
+#mag_data_rec = MagData(a, np.zeros((3,)+dim))
+#
+#mag_data_rec.mag_vec = x_f
+#
+#mag_data_rec.quiver_plot3d()
-phase_maps_rec = [pm(mag_data_rec) for pm in phasemappers]
-[phase_map.display_phase(title=u'Tilt series (rec.) $(\phi = {:2.1f} \pi)$'.format(tilts[i]/pi))
- for i, phase_map in enumerate(phase_maps_rec)]
+#phase_maps_rec = [pm(mag_data_rec) for pm in phasemappers]
+#[phase_map.display_phase(title=u'Tilt series (rec.) $(\phi = {:2.1f} \pi)$'.format(tilts[i%(count/2)]/pi))
+# for i, phase_map in enumerate(phase_maps_rec)]
+mag_data_opt = rc.optimize_sparse_cg(data_collection)
+mag_data_opt.quiver_plot3d()
-#
#first_guess = MagData(a, np.zeros((3,)+dim))
#
#first_guess.magnitude[1, int(dim[0]/2), int(dim[1]/2), int(dim[2]/2)] = 1
@@ -253,6 +285,77 @@ phase_maps_rec = [pm(mag_data_rec) for pm in phasemappers]
#phase_opt.display_phase()
+
+###################################################################################################
+print('--Singular value decomposition')
+
+#a = 1.
+#b_0 = 1000.
+#dim = (3, 3, 3)
+#count = 8
+#
+#magnitude = np.zeros((3,)+dim)
+#magnitude[:, int(dim[0]/2), int(dim[1]/2), int(dim[2]/2)] = 1.
+#
+#mag_data = MagData(a, magnitude)
+#mag_data.quiver_plot3d()
+#
+#tilts = np.linspace(0, 2*pi, num=count/2, endpoint=False)
+#projectors = []
+#projectors.extend([XTiltProjector(mag_data.dim, tilt) for tilt in tilts])
+#projectors.extend([YTiltProjector(mag_data.dim, tilt) for tilt in tilts])
+#phasemappers = [PMConvolve(mag_data.a, projector, b_0) for projector in projectors]
+#phase_maps = [pm(mag_data) for pm in phasemappers]
+#
+##[phase_map.display_phase(title=u'Tilt series $(\phi = {:2.1f} \pi)$'.format(tilts[i%(count/2)]/pi))
+## for i, phase_map in enumerate(phase_maps)]
+#
+#dim_uv = dim[1:3]
+#
+#lam = 10. ** -10
+#
+#size_2d = np.prod(dim_uv)
+#size_3d = np.prod(dim)
+#
+#
+#data_collection = DataCollection(a, dim_uv, b_0)
+#
+#[data_collection.append((phase_maps[i], projectors[i])) for i in range(count)]
+#
+#data = data_collection
+#y = data.phase_vec
+#kern = Kernel(data.a, data.dim_uv, data.b_0)
+#F = ForwardModel(data.projectors, kern)
+#C = Costfunction(y, F, lam)
+#
+#mag_data_opt = opt.optimize_sparse_cg(data_collection)
+#mag_data_opt.quiver_plot3d()
+
+
+
+#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))
+#A = MTM#np.array([F(np.eye(81)[:, i]) for i in range(81)])
+#
+#
+#
+#U, s, V = np.linalg.svd(M)
+##np.testing.assert_almost_equal(U.T, V, decimal=5)
+#
+#for value in range(20):
+# print 'Singular value:', s[value]
+# MagData(data.a, np.array(V[value,:]).reshape((3,)+dim)).quiver_plot3d()
+#
+#
+#for value in range(-10,0):
+# print 'Singular value:', s[value]
+# MagData(data.a, np.array(V[value,:]).reshape((3,)+dim)).quiver_plot3d()
+
+
+# TODO: print all singular vectors for a 2x2x2 distribution for each single and both tilt series!
+
+# TODO: Separate the script for SVD and compliance tests
+
####################################################################################################
#print('--Further testing')
#
@@ -281,7 +384,7 @@ phase_maps_rec = [pm(mag_data_rec) for pm in phasemappers]
#K = np.asmatrix(K)
#lam = 10. ** -10
#KTK = K.T * K + lam * np.asmatrix(np.eye(81))
-#print lam,
+#print lam,
##print pylab.cond(KTK),
#x_f = KTK.I * K.T * y
#print (np.asarray(K * x_f - y) ** 2).sum(),
diff --git a/tests/test_datacollection.py b/tests/test_dataset.py
similarity index 100%
rename from tests/test_datacollection.py
rename to tests/test_dataset.py
diff --git a/tests/test_optimize.py b/tests/test_reconstruction.py
similarity index 79%
rename from tests/test_optimize.py
rename to tests/test_reconstruction.py
index 4cc9862..7cb6fd3 100644
--- a/tests/test_optimize.py
+++ b/tests/test_reconstruction.py
@@ -5,10 +5,10 @@
import os
import unittest
-import pyramid.reconstructor as rc
+import pyramid.reconstruction as rc
-class TestCaseReconstructor(unittest.TestCase):
+class TestCaseReconstruction(unittest.TestCase):
def setUp(self):
self.path = os.path.join(os.path.dirname(os.path.realpath(__file__)), 'test_phasemapper')
@@ -21,5 +21,5 @@ class TestCaseReconstructor(unittest.TestCase):
if __name__ == '__main__':
- suite = unittest.TestLoader().loadTestsFromTestCase(TestCaseReconstructor)
+ suite = unittest.TestLoader().loadTestsFromTestCase(TestCaseReconstruction)
unittest.TextTestRunner(verbosity=2).run(suite)
--
GitLab