diff --git a/.hgignore b/.hgignore
index c604c2cf61c5cee59d0875b033475345d1eddf54..ffc7d0ed2f4cb6d9a488da1b409a6fb52a0b2e5d 100644
--- a/.hgignore
+++ b/.hgignore
@@ -4,6 +4,7 @@ syntax: glob
.spyderworkspace
.spyderproject
.cproject
+.pyramid.egg-info
*.pyc
*.pyd
*.so
diff --git a/pyramid/costfunction.py b/pyramid/costfunction.py
index e534283e1f5fc9b015a83999b866d5a887bf9a46..41fb1b5602d587bd8e2358f2a277fc93631d9edd 100644
--- a/pyramid/costfunction.py
+++ b/pyramid/costfunction.py
@@ -8,7 +8,6 @@ the so called `cost` of a threedimensional magnetization distribution."""
import numpy as np
-from pyramid.forwardmodel import ForwardModel
from pyramid.regularisator import NoneRegularisator
import logging
@@ -30,8 +29,6 @@ class Costfunction(object):
Attributes
----------
- data_set: :class:`~dataset.DataSet`
- :class:`~dataset.DataSet` object, which stores all information for the cost calculation.
regularisator : :class:`~.Regularisator`
Regularisator class that's responsible for the regularisation term.
y : :class:`~numpy.ndarray` (N=1)
@@ -51,14 +48,14 @@ class Costfunction(object):
_log = logging.getLogger(__name__+'.Costfunction')
- def __init__(self, data_set, regularisator):
+ def __init__(self, fwd_model, regularisator):
self._log.debug('Calling __init__')
- self.data_set = data_set
- self.fwd_model = ForwardModel(data_set)
+ self.fwd_model = fwd_model
self.regularisator = regularisator
if self.regularisator is None:
self.regularisator = NoneRegularisator()
# Extract important information:
+ data_set = fwd_model.data_set
self.y = data_set.phase_vec
self.n = data_set.n
self.m = data_set.m
@@ -69,13 +66,13 @@ class Costfunction(object):
def __repr__(self):
self._log.debug('Calling __repr__')
- return '%s(data_set=%r, regularisator=%r)' % \
- (self.__class__, self.data_set, self.regularisator)
+ return '%s(fwd_model=%r, regularisator=%r)' % \
+ (self.__class__, self.fwd_model, self.regularisator)
def __str__(self):
self._log.debug('Calling __str__')
- return 'Costfunction(data_set=%s, fwd_model=%s, regularisator=%s)' % \
- (self.data_set, self.fwd_model, self.regularisator)
+ return 'Costfunction(fwd_model=%s, fwd_model=%s, regularisator=%s)' % \
+ (self.fwd_model, self.fwd_model, self.regularisator)
def __call__(self, x):
delta_y = self.fwd_model(x) - self.y
diff --git a/pyramid/forwardmodel.py b/pyramid/forwardmodel.py
index e46f2560912342c3e57b93bf5fa3149942d1179d..41a9e6e8168dfaa654c53b957a1381612ff54bd8 100644
--- a/pyramid/forwardmodel.py
+++ b/pyramid/forwardmodel.py
@@ -136,6 +136,10 @@ class ForwardModel(object):
result[hp[i]:hp[i+1]] = res
return result
+ def _jac_dot_element(self, mag_vec, projector, phasemapper):
+ return phasemapper.jac_dot(projector.jac_dot(mag_vec))
+
+
def jac_T_dot(self, x, vector):
''''Calculate the product of the transposed Jacobi matrix with a given `vector`.
diff --git a/pyramid/magdata.py b/pyramid/magdata.py
index 0084ed345345a804093e9493eb19a3b9760becfd..af24642991304a16a39e486abb5b194dfaa7d6bb 100644
--- a/pyramid/magdata.py
+++ b/pyramid/magdata.py
@@ -5,12 +5,15 @@
"""This module provides the :class:`~.MagData` class for storing of magnetization data."""
+from __future__ import division
+
import os
import numpy as np
from numpy.linalg import norm
from scipy.ndimage.interpolation import zoom
+import matplotlib as mpl
import matplotlib.pyplot as plt
import matplotlib.cm as cmx
from matplotlib.ticker import MaxNLocator
@@ -164,6 +167,45 @@ class MagData(object):
self._log.debug('Calling __imul__')
return self.__mul__(other)
+ @classmethod
+ def _create_directional_colormap(cls, levels=15, N=256):
+ cls._log.debug('Calling __create_directional_colormap')
+ r, g, b = [], [], [] # RGB lists
+ # Create saturation lists to encode up and down directions via black and white colors.
+ # example for 5 levels from black (down) to color (in-plane) to white:
+ # pos_sat: [ 0. 0.5 1. 1. 1. ]
+ # neg_sat: [ 0. 0. 0. 0.5 1. ]
+ center = levels//2
+ pos_sat = np.ones(levels)
+ pos_sat[0:center] = [i/center for i in range(center)]
+ neg_sat = np.zeros(levels)
+ neg_sat[center+1:] = [(i+1)/center for i in range(center)]
+
+ # Iterate over all levels (the center level represents in-plane moments!):
+ for i in range(levels):
+ inter_points = np.linspace(i/levels, (i+1)/levels, 5) # interval points, current level
+ # Red:
+ r.append((inter_points[0], 0, neg_sat[i]))
+ r.append((inter_points[1], pos_sat[i], pos_sat[i]))
+ r.append((inter_points[2], pos_sat[i], pos_sat[i]))
+ r.append((inter_points[3], neg_sat[i], neg_sat[i]))
+ r.append((inter_points[4], neg_sat[i], 0))
+ # Green:
+ g.append((inter_points[0], 0, neg_sat[i]))
+ g.append((inter_points[1], neg_sat[i], neg_sat[i]))
+ g.append((inter_points[2], pos_sat[i], pos_sat[i]))
+ g.append((inter_points[3], pos_sat[i], pos_sat[i]))
+ g.append((inter_points[4], neg_sat[i], 0))
+ # Blue
+ b.append((inter_points[0], 0, pos_sat[i]))
+ b.append((inter_points[1], neg_sat[i], neg_sat[i]))
+ b.append((inter_points[2], neg_sat[i], neg_sat[i]))
+ b.append((inter_points[3], neg_sat[i], neg_sat[i]))
+ b.append((inter_points[4], pos_sat[i], 0))
+ # Combine to color dictionary and return:
+ cdict = {'red': r, 'green': g, 'blue': b}
+ return mpl.colors.LinearSegmentedColormap('directional_colormap', cdict, N)
+
def copy(self):
'''Returns a copy of the :class:`~.MagData` object
@@ -625,7 +667,7 @@ class MagData(object):
tree.write(filename, pretty_print=True)
def quiver_plot(self, title='Magnetization Distribution', axis=None, proj_axis='z',
- color='b', ar_dens=1, ax_slice=None, log=False, scaled=True, scale=1.):
+ coloring='angle', ar_dens=1, ax_slice=None, log=False, scaled=True, scale=1.):
'''Plot a slice of the magnetization as a quiver plot.
Parameters
@@ -636,8 +678,8 @@ class MagData(object):
Axis on which the graph is plotted. Creates a new figure if none is specified.
proj_axis : {'z', 'y', 'x'}, optional
The axis, from which a slice is plotted. The default is 'z'.
- color : string
- Color of the arrows. Default is black ('b').
+ coloring : string
+ Color coding mode of the arrows. Use 'angle' (default), 'amplitude' or 'uniform'.
ar_dens: int (optional)
Number defining the arrow density which is plotted. A higher ar_dens number skips more
arrows (a number of 2 plots every second arrow). Default is 1.
@@ -688,33 +730,53 @@ class MagData(object):
plot_v = np.swapaxes(np.copy(self.magnitude[1][..., ax_slice]), 0, 1) # y-component
u_label = 'z [px]'
v_label = 'y [px]'
+ # Prepare quiver (select only used arrows if ar_dens is specified):
+ dim_uv = plot_u.shape
+ yy, xx = np.indices(dim_uv)
+ xx = xx[::ar_dens, ::ar_dens]
+ yy = yy[::ar_dens, ::ar_dens]
+ plot_u = plot_u[::ar_dens, ::ar_dens]
+ plot_v = plot_v[::ar_dens, ::ar_dens]
+ amplitudes = np.hypot(plot_u, plot_v)
+ angles = np.angle(plot_u+1j*plot_v, deg=True)
+ # Calculate the arrow colors:
+ if coloring == 'angle':
+ self._log.debug('Encoding angles')
+ colors = (1 - np.arctan2(plot_v, plot_u)/np.pi) / 2 # in-plane angle (rescaled: 0 - 1)
+ cmap = self._create_directional_colormap(levels=1, N=256)
+ elif coloring == 'amplitude':
+ self._log.debug('Encoding amplitude')
+ colors = amplitudes / amplitudes.max()
+ cmap = 'jet'
+ elif coloring == 'uniform':
+ self._log.debug('No color encoding')
+ colors = np.zeros_like(plot_u) # use black arrows!
+ cmap = 'gray'
+ else:
+ raise AttributeError("Invalid coloring mode! Use 'angles', 'amplitude' or 'uniform'!")
# If no axis is specified, a new figure is created:
if axis is None:
self._log.debug('axis is None')
fig = plt.figure(figsize=(8.5, 7))
axis = fig.add_subplot(1, 1, 1)
axis.set_aspect('equal')
- angles = np.angle(plot_u+1j*plot_v, deg=True)
# Take the logarithm of the arrows to clearly show directions (if specified):
if log:
cutoff = 10
- amp = np.round(np.hypot(plot_u, plot_v), decimals=cutoff)
+ amp = np.round(amplitudes, decimals=cutoff)
min_value = amp[np.nonzero(amp)].min()
plot_u = np.round(plot_u, decimals=cutoff) / min_value
plot_u = np.log10(np.abs(plot_u)+1) * np.sign(plot_u)
plot_v = np.round(plot_v, decimals=cutoff) / min_value
plot_v = np.log10(np.abs(plot_v)+1) * np.sign(plot_v)
+ amplitudes = np.hypot(plot_u, plot_v) # Recalculate!
# Scale the magnitude of the arrows to the highest one (if specified):
if scaled:
- plot_u /= np.hypot(plot_u, plot_v).max()
- plot_v /= np.hypot(plot_u, plot_v).max()
- # Setup quiver:
- dim_uv = plot_u.shape
- ad = ar_dens
- yy, xx = np.indices(dim_uv)
- axis.quiver(xx[::ad, ::ad], yy[::ad, ::ad], plot_u[::ad, ::ad], plot_v[::ad, ::ad],
- pivot='middle', units='xy', angles=angles[::ad, ::ad], minlength=0.25,
- scale_units='xy', scale=scale/ad, headwidth=6, headlength=7, color=color)
+ plot_u /= amplitudes.max()
+ plot_v /= amplitudes.max()
+ axis.quiver(xx, yy, plot_u, plot_v, colors, cmap=cmap, angles=angles,
+ pivot='middle', units='xy', scale_units='xy', scale=scale/ar_dens,
+ minlength=0.25, headwidth=6, headlength=7)
axis.set_xlim(-1, dim_uv[1])
axis.set_ylim(-1, dim_uv[0])
axis.set_title(title, fontsize=18)
@@ -726,7 +788,7 @@ class MagData(object):
return axis
def quiver_plot3d(self, title='Magnetization Distribution', limit=None, cmap='jet',
- ar_dens=1, mode='arrow', show_pipeline=False):
+ ar_dens=1, mode='arrow', coloring='angle', show_pipeline=False):
'''Plot the magnetization as 3D-vectors in a quiverplot.
Parameters
@@ -743,6 +805,8 @@ class MagData(object):
mode: string, optional
Mode, determining the glyphs used in the 3D plot. Default is 'arrow', which corresponds
to 3D arrows. For large amounts of arrows, '2darrow' should be used.
+ coloring : string
+ Color coding mode of the arrows. Use 'angle' (default) or 'amplitude'.
show_pipeline : boolean, optional
If True, the mayavi pipeline, a GUI used for image manipulation is shown. The default
is False.
@@ -771,13 +835,36 @@ class MagData(object):
# Plot them as vectors:
mlab.figure(size=(750, 700))
plot = mlab.quiver3d(xx, yy, zz, x_mag, y_mag, z_mag, mode=mode, colormap=cmap)
+ if coloring == 'angle': # Encodes the full angle via colorwheel and saturation
+ self._log.debug('Encoding full 3D angles')
+ from tvtk.api import tvtk
+ phis = (1 - np.arctan2(y_mag, x_mag)/np.pi) / 2 # in-plane angle (rescaled: 0 to 1)
+ thetas = np.arctan2(np.hypot(y_mag, x_mag), z_mag)/np.pi # vert. angle (also rescaled)
+ levels = 15 # number of shades for up/down arrows (white is up, black is down)
+ N = 256 # Overall number of colors for the colormap (only N/levels per level!)
+ cmap = self._create_directional_colormap(levels, N)
+ colors = []
+ for i in range(len(x_mag)):
+ level = np.floor((1-thetas[i]) * levels)
+ lookup_value = (level + phis[i]) / levels
+ rgba = tuple((np.asarray(cmap(lookup_value)) * (N-1)).astype(np.int))
+ colors.append(rgba)
+ sc = tvtk.UnsignedCharArray() # Used to hold the colors
+ sc.from_array(colors)
+ plot.mlab_source.dataset.point_data.scalars = sc
+ plot.mlab_source.dataset.modified()
+ plot.glyph.color_mode = 'color_by_scalar'
+ elif coloring == 'amplitude': # Encodes the amplitude of the arrows with the jet colormap
+ self._log.debug('Encoding amplitude')
+ mlab.colorbar(label_fmt='%.2f')
+ mlab.colorbar(orientation='vertical')
+ else:
+ raise AttributeError('Coloring mode not supported!')
plot.glyph.glyph_source.glyph_position = 'center'
plot.module_manager.vector_lut_manager.data_range = np.array([0, limit])
mlab.outline(plot)
mlab.axes(plot)
mlab.title(title, height=0.95, size=0.35)
- mlab.colorbar(label_fmt='%.2f')
- mlab.colorbar(orientation='vertical')
mlab.orientation_axes()
if show_pipeline:
mlab.show_pipeline()
diff --git a/pyramid/reconstruction.py b/pyramid/reconstruction.py
index f5df52a5f6325f35d096be2f4f180f1594d1f150..96664e82d46da77c7416315d443164c36c645dae 100644
--- a/pyramid/reconstruction.py
+++ b/pyramid/reconstruction.py
@@ -15,7 +15,6 @@ the distribution.
import numpy as np
-from pyramid.costfunction import Costfunction
from pyramid.magdata import MagData
import logging
@@ -25,20 +24,18 @@ __all__ = ['optimize_linear', 'optimize_nonlin', 'optimize_splitbregman']
_log = logging.getLogger(__name__)
-def optimize_linear(data, regularisator=None, max_iter=None):
+def optimize_linear(costfunction, max_iter=None):
'''Reconstruct a three-dimensional magnetic distribution from given phase maps via the
conjugate gradient optimizaion method :func:`~.scipy.sparse.linalg.cg`.
Blazingly fast for l2-based cost functions.
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.
- regularisator : :class:`~.Regularisator`, optional
- Regularisator class that's responsible for the regularisation term. Defaults to zero
- order Tikhonov if none is provided.
+ costfunction : :class:`~.Costfunction`
+ A :class:`~.Costfunction` object which implements a specified forward model and
+ regularisator which is minimized in the optimization process.
+ max_iter : int, optional
+ The maximum number of iterations for the opimization.
Returns
-------
@@ -48,32 +45,28 @@ def optimize_linear(data, regularisator=None, max_iter=None):
'''
import jutil.cg as jcg
_log.debug('Calling optimize_linear')
- # Set up necessary objects:
- cost = Costfunction(data, regularisator)
- _log.info('Cost before optimization: {}'.format(cost(np.zeros(cost.n))))
- x_opt = jcg.conj_grad_minimize(cost, max_iter=max_iter).x
- _log.info('Cost after optimization: {}'.format(cost(x_opt)))
+ _log.info('Cost before optimization: {}'.format(costfunction(np.zeros(costfunction.n))))
+ x_opt = jcg.conj_grad_minimize(costfunction, max_iter=max_iter).x
+ _log.info('Cost after optimization: {}'.format(costfunction(x_opt)))
# Create and return fitting MagData object:
- mag_opt = MagData(data.a, np.zeros((3,) + data.dim))
- mag_opt.set_vector(x_opt, data.mask)
- return mag_opt, cost
+ data_set = costfunction.fwd_model.data_set
+ mag_opt = MagData(data_set.a, np.zeros((3,) + data_set.dim))
+ mag_opt.set_vector(x_opt, data_set.mask)
+ return mag_opt
-def optimize_nonlin(data, first_guess=None, regularisator=None):
+def optimize_nonlin(costfunction, first_guess=None):
'''Reconstruct a three-dimensional magnetic distribution from given phase maps via
steepest descent method. This is slow, but works best for non l2-regularisators.
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.
+ costfunction : :class:`~.Costfunction`
+ A :class:`~.Costfunction` object which implements a specified forward model and
+ regularisator which is minimized in the optimization process.
first_guess : :class:`~pyramid.magdata.MagData`
magnetization to start the non-linear iteration with.
- regularisator : :class:`~.Regularisator`, optional
- Regularisator class that's responsible for the regularisation term.
Returns
-------
@@ -84,14 +77,14 @@ def optimize_nonlin(data, first_guess=None, regularisator=None):
import jutil.minimizer as jmin
import jutil.norms as jnorms
_log.debug('Calling optimize_nonlin')
+ data_set = costfunction.fwd_model.data_set
if first_guess is None:
- first_guess = MagData(data.a, np.zeros((3,) + data.dim))
+ first_guess = MagData(data_set.a, np.zeros((3,) + data_set.dim))
- x_0 = first_guess.get_vector(data.mask)
- cost = Costfunction(data, regularisator)
- assert len(x_0) == cost.n, (len(x_0), cost.m, cost.n)
+ x_0 = first_guess.get_vector(data_set.mask)
+ assert len(x_0) == costfunction.n, (len(x_0), costfunction.m, costfunction.n)
- p = regularisator.p
+ p = costfunction.regularisator.p
q = 1. / (1. - (1. / p))
lq = jnorms.LPPow(q, 1e-20)
@@ -101,20 +94,20 @@ def optimize_nonlin(data, first_guess=None, regularisator=None):
return direc_p
# This Method is semi-best for Lp type problems. Takes forever, though
- _log.info('Cost before optimization: {}'.format(cost(np.zeros(cost.n))))
+ _log.info('Cost before optimization: {}'.format(costfunction(np.zeros(costfunction.n))))
result = jmin.minimize(
- cost, x_0,
+ costfunction, x_0,
method="SteepestDescent",
options={"preconditioner": preconditioner},
tol={"max_iteration": 10000})
x_opt = result.x
- _log.info('Cost after optimization: {}'.format(cost(x_opt)))
- mag_opt = MagData(data.a, np.zeros((3,) + data.dim))
- mag_opt.set_vector(x_opt, data.mask)
+ _log.info('Cost after optimization: {}'.format(costfunction(x_opt)))
+ mag_opt = MagData(data_set.a, np.zeros((3,) + data_set.dim))
+ mag_opt.set_vector(x_opt, data_set.mask)
return mag_opt
-def optimize_splitbregman(data, weight, lam, mu):
+def optimize_splitbregman(costfunction, weight, lam, mu):
'''
Reconstructs magnet distribution from phase image measurements using a split bregman
algorithm with a dedicated TV-l1 norm. Very dedicated, frickle, brittle, and difficult
@@ -124,10 +117,9 @@ def optimize_splitbregman(data, weight, lam, mu):
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.
+ costfunction : :class:`~.Costfunction`
+ A :class:`~.Costfunction` object which implements a specified forward model and
+ regularisator which is minimized in the optimization process.
weight : float
Obscure split bregman parameter
lam : float
@@ -144,29 +136,27 @@ def optimize_splitbregman(data, weight, lam, mu):
import jutil.splitbregman as jsb
import jutil.operator as joperator
import jutil.diff as jdiff
- from pyramid.regularisator import FirstOrderRegularisator
_log.debug('Calling optimize_splitbregman')
# regularisator is actually not necessary, but this makes the cost
# function to that which is supposedly optimized by split bregman.
# Thus cost can be used to verify convergence
- regularisator = FirstOrderRegularisator(data.mask, lam / mu, 1)
- cost = Costfunction(data, regularisator)
- fwd_mod = cost.fwd_model
+ fwd_model = costfunction.fwd_model
+ data_set = fwd_model.data_set
A = joperator.Function(
- (cost.m, cost.n),
- lambda x: fwd_mod.jac_dot(None, x),
- FT=lambda x: fwd_mod.jac_T_dot(None, x))
+ (costfunction.m, costfunction.n),
+ lambda x: fwd_model.jac_dot(None, x),
+ FT=lambda x: fwd_model.jac_T_dot(None, x))
D = joperator.VStack([
- jdiff.get_diff_operator(data.mask, 0, 3),
- jdiff.get_diff_operator(data.mask, 1, 3)])
- y = np.asarray(cost.y, dtype=np.double)
+ jdiff.get_diff_operator(data_set.mask, 0, 3),
+ jdiff.get_diff_operator(data_set.mask, 1, 3)])
+ y = np.asarray(costfunction.y, dtype=np.double)
x_opt = jsb.split_bregman_2d(
A, D, y,
weight=weight, mu=mu, lambd=lam, max_iter=1000)
- mag_opt = MagData(data.a, np.zeros((3,) + data.dim))
- mag_opt.set_vector(x_opt, data.mask)
+ mag_opt = MagData(data_set.a, np.zeros((3,) + data_set.dim))
+ mag_opt.set_vector(x_opt, data_set.mask)
return mag_opt
diff --git a/pyramid/regularisator.py b/pyramid/regularisator.py
index db23e8cd9f79ef443a7ed18b9540402eafba908f..a6fac9512986c9423e1089b5ade4ce6fb3b3525e 100644
--- a/pyramid/regularisator.py
+++ b/pyramid/regularisator.py
@@ -9,6 +9,7 @@ which adds additional constraints to a costfunction to minimize."""
import abc
import numpy as np
+from scipy import sparse
import jutil.norms as jnorm
import jutil.diff as jdiff
@@ -17,7 +18,8 @@ import jutil.operator as joperator
import logging
-__all__ = ['NoneRegularisator', 'ZeroOrderRegularisator', 'FirstOrderRegularisator']
+__all__ = ['NoneRegularisator', 'ZeroOrderRegularisator', 'FirstOrderRegularisator',
+ 'ComboRegularisator']
class Regularisator(object):
@@ -117,6 +119,97 @@ class Regularisator(object):
return self.lam * self.norm.hess_diag(x)
+class ComboRegularisator(Regularisator):
+
+ '''Class for providing a regularisation term which combines several regularisators.
+
+ If more than one regularisation should be utilized, this class can be use. It is given a list
+ of :class:`~.Regularisator` objects. The input will be forwarded to each of them and the
+ results are summed up and returned.
+
+ Attributes
+ ----------
+ reg_list: :class:`~.Regularisator`
+ A list of regularisator objects to whom the input is passed on.
+
+ '''
+
+ def __init__(self, reg_list):
+ self._log.debug('Calling __init__')
+ self.reg_list = reg_list
+ self._log.debug('Created '+str(self))
+
+ def __call__(self, x):
+ self._log.debug('Calling __call__')
+ return np.sum([self.reg_list[i](x) for i in range(len(self.reg_list))], axis=0)
+
+ def __repr__(self):
+ self._log.debug('Calling __repr__')
+ return '%s(reg_list=%r)' % (self.__class__, self.reg_list)
+
+ def __str__(self):
+ self._log.debug('Calling __str__')
+ return 'ComboRegularisator(reg_list=%s)' % (self.reg_list)
+
+ def jac(self, x):
+ '''Calculate the derivative of the regularisation term for a given magnetic distribution.
+
+ Parameters
+ ----------
+ x: :class:`~numpy.ndarray` (N=1)
+ Vectorized magnetization distribution, for which the Jacobi vector is calculated.
+
+ Returns
+ -------
+ result : :class:`~numpy.ndarray` (N=1)
+ Jacobi vector which represents the cost derivative of all voxels of the magnetization.
+
+ '''
+ return np.sum([self.reg_list[i].jac(x) for i in range(len(self.reg_list))], axis=0)
+
+ def hess_dot(self, x, vector):
+ '''Calculate the product of a `vector` with the Hessian matrix of the regularisation term.
+
+ Parameters
+ ----------
+ x : :class:`~numpy.ndarray` (N=1)
+ Vectorized magnetization distribution at which the Hessian is calculated. The Hessian
+ is constant in this case, thus `x` can be set to None (it is not used int the
+ computation). It is implemented for the case that in the future nonlinear problems
+ have to be solved.
+ vector : :class:`~numpy.ndarray` (N=1)
+ Vectorized magnetization distribution which is multiplied by the Hessian.
+
+ Returns
+ -------
+ result : :class:`~numpy.ndarray` (N=1)
+ Product of the input `vector` with the Hessian matrix.
+
+ '''
+ return np.sum([self.reg_list[i].hess_dot(x, vector) for i in range(len(self.reg_list))],
+ axis=0)
+
+ def hess_diag(self, x):
+ ''' Return the diagonal of the Hessian.
+
+ 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 in the
+ computation). It is implemented for the case that in the future nonlinear problems
+ have to be solved.
+
+ Returns
+ -------
+ result : :class:`~numpy.ndarray` (N=1)
+ Diagonal of the Hessian matrix.
+
+ '''
+ self._log.debug('Calling hess_diag')
+ return np.sum([self.reg_list[i].hess_diag(x) for i in range(len(self.reg_list))], axis=0)
+
+
class NoneRegularisator(Regularisator):
'''Placeholder class if no regularization is used.
@@ -255,86 +348,10 @@ class FirstOrderRegularisator(Regularisator):
D0 = jdiff.get_diff_operator(mask, 0, 3)
D1 = jdiff.get_diff_operator(mask, 1, 3)
D2 = jdiff.get_diff_operator(mask, 2, 3)
- D = joperator.VStack([D0, D1, D2])
+ D = sparse.vstack([D0, D1, D2])
if p == 2:
norm = jnorm.WeightedL2Square(D)
else:
norm = jnorm.WeightedTV(jnorm.LPPow(p, 1e-12), D, [D0.shape[0], D.shape[0]])
super(FirstOrderRegularisator, self).__init__(norm, lam)
self._log.debug('Created '+str(self))
-
-# TODO: ComboRegularisator ########################################################################
-#class ComboRegularisator(Regularisator):
-#
-# '''Class for providing a regularisation term which implements derivation minimization.
-#
-# The constraint this class represents is the minimization of the first order derivative of the
-# 3D magnetization distribution using a Lp norm. Important is the regularisation parameter `lam`
-# (lambda) which determines the weighting between the two cost parts (measurements and
-# regularisation).
-#
-# Attributes
-# ----------
-# mask: :class:`~numpy.ndarray` (N=3)
-# A boolean mask which defines the magnetized volume in 3D.
-# lam: float
-# Regularisation parameter determining the weighting between measurements and regularisation.
-# p: int, optional
-# Order of the norm (default: 2, which means a standard L2-norm).
-#
-# '''
-#
-# def __init__(self, norm, lam):
-# self._log.debug('Calling __init__')
-# self.norm = norm
-# self.lam = lam
-# self._log.debug('Created '+str(self))
-#
-# def __call__(self, x):
-# self._log.debug('Calling __call__')
-# return self.lam * self.norm(x)
-#
-# def __repr__(self):
-# self._log.debug('Calling __repr__')
-# return '%s(norm=%r, lam=%r)' % (self.__class__, self.norm, self.lam)
-#
-# def __str__(self):
-# self._log.debug('Calling __str__')
-# return 'Regularisator(norm=%s, lam=%s)' % (self.norm, self.lam)
-#
-# def jac(self, x):
-# '''Calculate the derivative of the regularisation term for a given magnetic distribution.
-#
-# Parameters
-# ----------
-# x: :class:`~numpy.ndarray` (N=1)
-# Vectorized magnetization distribution, for which the Jacobi vector is calculated.
-#
-# Returns
-# -------
-# result : :class:`~numpy.ndarray` (N=1)
-# Jacobi vector which represents the cost derivative of all voxels of the magnetization.
-#
-# '''
-# return self.lam * self.norm.jac(x)
-#
-# def hess_dot(self, x, vector):
-# '''Calculate the product of a `vector` with the Hessian matrix of the regularisation term.
-#
-# Parameters
-# ----------
-# x : :class:`~numpy.ndarray` (N=1)
-# Vectorized magnetization distribution at which the Hessian is calculated. The Hessian
-# is constant in this case, thus `x` can be set to None (it is not used int the
-# computation). It is implemented for the case that in the future nonlinear problems
-# have to be solved.
-# vector : :class:`~numpy.ndarray` (N=1)
-# Vectorized magnetization distribution which is multiplied by the Hessian.
-#
-# Returns
-# -------
-# result : :class:`~numpy.ndarray` (N=1)
-# Product of the input `vector` with the Hessian matrix.
-#
-# '''
-# return self.lam * self.norm.hess_dot(x, vector)
diff --git a/pyramid/version.py b/pyramid/version.py
index 577dd56a08ab9487a926c851d44022c997397fa7..bf68b518888f618528c17cf7151d0fb8ee2c2678 100644
--- a/pyramid/version.py
+++ b/pyramid/version.py
@@ -1,3 +1,3 @@
# THIS FILE IS GENERATED BY THE PYRAMID SETUP.PY
version = "0.1.0-dev"
-hg_revision = "f641a9bb194f+"
+hg_revision = "e2e912d70bc2+"
diff --git a/scripts/magdata/magcreator/create_homog_slab.py b/scripts/magdata/magcreator/create_homog_slab.py
index 520712099f166c0549ce76ee6e2f51f36ad50652..9fc206d7f8d4145afcfe1b2660b15cf518858ed0 100644
--- a/scripts/magdata/magcreator/create_homog_slab.py
+++ b/scripts/magdata/magcreator/create_homog_slab.py
@@ -1,6 +1,6 @@
#! python
# -*- coding: utf-8 -*-
-"""Create pyramid logo."""
+"""Create homogeneous slab."""
import numpy as np
diff --git a/scripts/magdata/magcreator/create_singularity.py b/scripts/magdata/magcreator/create_singularity.py
new file mode 100644
index 0000000000000000000000000000000000000000..d050b0e1e0038e0f883c2630c2a09642e955e112
--- /dev/null
+++ b/scripts/magdata/magcreator/create_singularity.py
@@ -0,0 +1,26 @@
+# -*- coding: utf-8 -*-
+"""magnetic singularity"""
+
+
+import numpy as np
+import pyramid as py
+import logging.config
+
+
+logging.config.fileConfig(py.LOGGING_CONFIG, disable_existing_loggers=False)
+
+# Parameters:
+dim = (5, 5, 5)
+center = (2., 2., 2.)
+a = 1.0 # in nm
+filename = 'magdata_mc_singularity.nc'
+
+magnitude = np.zeros((3,) + dim)
+zz, yy, xx = np.indices(dim)
+magnitude = np.array((xx-center[2], yy-center[1], zz-center[0]))
+magnitude /= np.sqrt((magnitude**2).sum(axis=0))
+
+# Create and save MagData object:
+mag_data = py.MagData(a, magnitude)
+mag_data.quiver_plot3d(coloring='full angle')
+mag_data.save_to_netcdf4(filename)
diff --git a/scripts/magdata/magdata_from_dat.py b/scripts/magdata/magdata_from_dat.py
new file mode 100644
index 0000000000000000000000000000000000000000..4a1f176e51e0433c75d4c52366c22097385527ee
--- /dev/null
+++ b/scripts/magdata/magdata_from_dat.py
@@ -0,0 +1,68 @@
+# -*- coding: utf-8 -*-
+"""Create magnetization distributions from fortran sorted txt-files."""
+
+
+import os
+import numpy as np
+import pyramid as py
+import logging.config
+
+
+logging.config.fileConfig(py.LOGGING_CONFIG, disable_existing_loggers=False)
+
+###################################################################################################
+filename = 'J=1.D=0.084.H=0.0067.Bobber.dat'#'J=1.D=0.084.H=0.00353.SkLatt.dat'
+scale = 1
+###################################################################################################
+
+#@profile
+def func():
+ path = os.path.join(py.DIR_FILES, 'dat', filename)
+# data = np.genfromtxt(path, dtype=np.float32, delimiter=',', usecols=(0, 1, 2, 3, 4, 5))
+ x, y, z, xmag, ymag, zmag = data.T
+ a = (y[1] - y[0]) * scale
+ dim_z = len(np.unique(z))
+ dim_y = len(np.unique(y))
+ dim_x = len(np.unique(x))
+ dim = (dim_z, dim_x, dim_y) # Order of descending variance!
+ xmag = xmag.reshape(dim).swapaxes(1, 2)
+ ymag = ymag.reshape(dim).swapaxes(1, 2)
+ zmag = zmag.reshape(dim).swapaxes(1, 2)
+ magnitude = np.array((xmag, ymag, zmag))
+ mag_data = py.MagData(a, magnitude)
+ mag_data.save_to_netcdf4('magdata_dat_{}'.format(filename.replace('.dat', '.nc')))
+# mag_data.quiver_plot3d(ar_dens=4, coloring='amplitude')
+# mag_data.quiver_plot3d(ar_dens=4, coloring='angle')
+# py.pm(mag_data).display_combined(interpolation='bilinear')
+
+
+@profile
+def funci():
+ path = os.path.join(py.DIR_FILES, 'dat', filename)
+ with open(path) as f:
+ import csv
+ reader = csv.reader(f)
+ x, y, z, xmag, ymag, zmag = [], [], [], [], [], []
+ for row in reader:
+ x.append(int(row[0]))
+ y.append(int(row[1]))
+ z.append(int(row[2]))
+ xmag.append(float(row[3]))
+ ymag.append(float(row[4]))
+ zmag.append(float(row[5]))
+ a = (y[1] - y[0]) * scale
+ dim_z = len(np.unique(z))
+ dim_y = len(np.unique(y))
+ dim_x = len(np.unique(x))
+ dim = (dim_z, dim_x, dim_y) # Order of descending variance!
+ xmag = np.reshape(xmag, dim).swapaxes(1, 2)
+ ymag = np.reshape(ymag, dim).swapaxes(1, 2)
+ zmag = np.reshape(zmag, dim).swapaxes(1, 2)
+ magnitude = np.array((xmag, ymag, zmag))
+ mag_data = py.MagData(a, magnitude)
+# mag_data.save_to_netcdf4('magdata_dat_{}'.format(filename.replace('.dat', '.nc')))
+
+if __name__ == '__main__':
+ funci()
+
+
diff --git a/scripts/temp/custom_colormap.py b/scripts/temp/custom_colormap.py
new file mode 100644
index 0000000000000000000000000000000000000000..2ed480ab304c64013013c70fcc9c5681063f8d8f
--- /dev/null
+++ b/scripts/temp/custom_colormap.py
@@ -0,0 +1,174 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Sun Jun 21 12:14:28 2015
+
+@author: Jan
+"""
+
+from __future__ import division
+import numpy as np
+import matplotlib as mpl
+import matplotlib.pyplot as plt
+from mpl_toolkits.mplot3d import Axes3D
+
+import logging
+
+
+__all__ = ['create_custom_colormap']
+_log = logging.getLogger(__name__)
+
+
+def create_custom_colormap(levels=15, N=256):
+
+ '''Class for storing phase map data.
+
+ 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
+ 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
+ corresponding holographic contour map are provided. Holographic contour maps are created by
+ taking the cosine of the (optionally amplified) phase and encoding the direction of the
+ 2-dimensional gradient via color. The directional encoding can be seen by using the
+ :func:`~.make_color_wheel` function. Use the :func:`~.display_combined` function to plot the
+ phase map and the holographic contour map next to each other.
+
+ Attributes
+ ----------
+ a: float
+ The grid spacing in nm.
+ dim_uv: tuple (N=2)
+ Dimensions of the grid.
+ phase: :class:`~numpy.ndarray` (N=2)
+ Matrix containing the phase shift.
+ phase_vec: :class:`~numpy.ndarray` (N=2)
+ Vector containing the phase shift.
+ mask: :class:`~numpy.ndarray` (boolean, N=2, optional)
+ Mask which determines the projected magnetization distribution, gotten from MIP images or
+ otherwise acquired. Defaults to an array of ones (all pixels are considered).
+ confidence: :class:`~numpy.ndarray` (N=2, optional)
+ Confidence array which determines the trust of specific regions of the phase_map. A value
+ of 1 means the pixel is trustworthy, a value of 0 means it is not. Defaults to an array of
+ ones (full trust for all pixels). Can be used for the construction of Se_inv.
+ unit: {'rad', 'mrad'}, optional
+ Set the unit of the phase map. This is important for the :func:`display` function,
+ because the phase is scaled accordingly. Does not change the phase itself, which is
+ always in `rad`.
+
+ ''' # TODO: Docstring!
+
+ CDICT = {'red': [(0.00, 1.0, 0.0),
+ (0.25, 1.0, 1.0),
+ (0.50, 1.0, 1.0),
+ (0.75, 0.0, 0.0),
+ (1.00, 0.0, 1.0)],
+
+ 'green': [(0.00, 0.0, 0.0),
+ (0.25, 0.0, 0.0),
+ (0.50, 1.0, 1.0),
+ (0.75, 1.0, 1.0),
+ (1.00, 0.0, 1.0)],
+
+ 'blue': [(0.00, 1.0, 1.0),
+ (0.25, 0.0, 0.0),
+ (0.50, 0.0, 0.0),
+ (0.75, 0.0, 0.0),
+ (1.00, 1.0, 1.0)]}
+
+ CDICT_INV = {'red': [(0.00, 0.0, 1.0),
+ (0.25, 0.0, 0.0),
+ (0.50, 0.0, 0.0),
+ (0.75, 1.0, 1.0),
+ (1.00, 1.0, 0.0)],
+
+ 'green': [(0.00, 1.0, 1.0),
+ (0.25, 1.0, 1.0),
+ (0.50, 0.0, 0.0),
+ (0.75, 0.0, 0.0),
+ (1.00, 1.0, 0.0)],
+
+ 'blue': [(0.00, 0.0, 0.0),
+ (0.25, 1.0, 1.0),
+ (0.50, 1.0, 1.0),
+ (0.75, 1.0, 1.0),
+ (1.00, 0.0, 0.0)]}
+
+ r, g, b = [], [], []
+ center = levels//2
+ pos_sat = np.ones(levels)
+ pos_sat[0:center] = [i/center for i in range(center)]
+ neg_sat = np.zeros(levels)
+ neg_sat[center+1:] = [(i+1)/center for i in range(center)]
+ print pos_sat
+ print neg_sat
+ # example for 5 levels (from black to color to white):
+ # [ 0. 0.5 1. 1. 1. ]
+ # [ 0. 0. 0. 0.5 1. ]
+
+ for i in range(levels):
+ print i + 1
+ inter_points = np.linspace(i/levels, (i+1)/levels, 5) # interval points
+ print inter_points
+
+ CDICT = {'red': [(0.00, 1.0, 0.0),
+ (0.25, 1.0, 1.0),
+ (0.50, 1.0, 1.0),
+ (0.75, 0.0, 0.0),
+ (1.00, 0.0, 1.0)],
+
+ 'green': [(0.00, 0.0, 0.0),
+ (0.25, 0.0, 0.0),
+ (0.50, 1.0, 1.0),
+ (0.75, 1.0, 1.0),
+ (1.00, 0.0, 1.0)],
+
+ 'blue': [(0.00, 1.0, 1.0),
+ (0.25, 0.0, 0.0),
+ (0.50, 0.0, 0.0),
+ (0.75, 0.0, 0.0),
+ (1.00, 1.0, 1.0)]}
+
+ r.append((inter_points[0], 0, neg_sat[i]))
+ r.append((inter_points[1], pos_sat[i], pos_sat[i]))
+ r.append((inter_points[2], pos_sat[i], pos_sat[i]))
+ r.append((inter_points[3], neg_sat[i], neg_sat[i]))
+ r.append((inter_points[4], neg_sat[i], 0))
+ print r
+
+ g.append((inter_points[0], 0, neg_sat[i]))
+ g.append((inter_points[1], neg_sat[i], neg_sat[i]))
+ g.append((inter_points[2], pos_sat[i], pos_sat[i]))
+ g.append((inter_points[3], pos_sat[i], pos_sat[i]))
+ g.append((inter_points[4], neg_sat[i], 0))
+ print g
+
+ b.append((inter_points[0], 0, pos_sat[i]))
+ b.append((inter_points[1], neg_sat[i], neg_sat[i]))
+ b.append((inter_points[2], neg_sat[i], neg_sat[i]))
+ b.append((inter_points[3], neg_sat[i], neg_sat[i]))
+ b.append((inter_points[4], pos_sat[i], 0))
+ print b
+
+# r.append((inter_points[0], 0, 0))
+# r.append((inter_points[1], pos_sat[i], pos_sat[i]))
+#
+# g.append((inter_points[0], 0, 0))
+# g.append((inter_points[1], neg_sat[i], neg_sat[i]))
+#
+# b.append((inter_points[0], 0, 0))
+# b.append((inter_points[1], neg_sat[i], neg_sat[i]))
+
+ cdict = {'red': r, 'green': g, 'blue': b}
+ return mpl.colors.LinearSegmentedColormap('custom_colormap', cdict, N)
+
+
+x = np.arange(0, np.pi, 0.1)
+y = np.arange(0, 2*np.pi, 0.1)
+X, Y = np.meshgrid(x, y)
+Z = np.cos(X) * np.sin(Y) * 10
+
+fig = plt.figure(figsize=(7, 7))
+axis = fig.add_subplot(1, 1, 1)
+im = axis.imshow(Z, cmap=create_custom_colormap(levels=3, N=3))
+fig.colorbar(im)
diff --git a/scripts/temp/quiver_plot_enhanced.py b/scripts/temp/quiver_plot_enhanced.py
new file mode 100644
index 0000000000000000000000000000000000000000..a30deae41df8023121c54debe82829a34ab41643
--- /dev/null
+++ b/scripts/temp/quiver_plot_enhanced.py
@@ -0,0 +1,127 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Fri Jun 19 13:26:32 2015
+
+@author: Jan
+"""
+
+
+from mayavi import mlab
+import numpy as np
+import pyramid as py
+
+
+mag_data = py.MagData.load_from_netcdf4('magdata_mc_orthogonal_thin_vortices.nc')
+#mag_data = py.MagData(1, np.ones((3, 2, 2, 2)))
+a = mag_data.a
+dim = mag_data.dim
+magnitude = mag_data.magnitude
+ad = 1
+limit = 1
+
+byte = 256
+
+def create_8bit_rgb_lut():
+ xl = np.mgrid[0:byte, 0:byte, 0:byte]
+ lut = np.vstack((xl[0].reshape(1, byte**3),
+ xl[1].reshape(1, byte**3),
+ xl[2].reshape(1, byte**3),
+ 255 * np.ones((1, byte**3)))).T
+ return lut.astype('int32')
+
+
+# indexing function to above grid
+def rgb_2_scalar_idx(r, g, b):
+ print 'r:', r, 'g:', g, 'b:', b
+ return byte**2 * r + byte * g + b
+
+# N x 3 colors
+colors = np.arange(3*np.prod(dim)).reshape((-1, 3)) % byte
+print colors
+
+# N scalars
+scalars = np.zeros(np.prod(dim))
+for (kp_idx, kp_c) in enumerate(colors):
+ print kp_idx, kp_c
+ scalars[kp_idx] = rgb_2_scalar_idx(kp_c[0], kp_c[1], kp_c[2])
+print scalars
+
+#rgb_lut = create_8bit_rgb_lut()
+
+# Create points and vector components as lists:
+zz, yy, xx = (np.indices(dim)-a/2).reshape((3,)+dim)
+zz = zz[::ad, ::ad, ::ad].flatten()
+yy = yy[::ad, ::ad, ::ad].flatten()
+xx = xx[::ad, ::ad, ::ad].flatten()
+x_mag = magnitude[0][::ad, ::ad, ::ad].flatten()
+y_mag = magnitude[1][::ad, ::ad, ::ad].flatten()
+z_mag = magnitude[2][::ad, ::ad, ::ad].flatten()
+# Plot them as vectors:
+mlab.figure(size=(750, 700))
+
+phis = (1 - np.arctan2(y_mag, x_mag)/np.pi) / 2
+thetas = np.arctan2(np.hypot(y_mag, x_mag), z_mag)/np.pi
+
+
+levels = 15
+N = 256
+cmap = py.MagData._create_directional_colormap(levels, N)
+
+
+def angle_to_rgba(phi, theta):
+ level = np.floor((1-theta) * levels)
+ lookup_value = (level + phi) / levels
+ rgba = cmap(lookup_value)
+ return tuple((np.asarray(rgba)*255).astype(np.int))
+
+colors = []
+for i in range(len(x_mag)):
+ colors.append(angle_to_rgba(phis[i], thetas[i]))
+
+plot = mlab.quiver3d(xx, yy, zz, x_mag, y_mag, z_mag, mode='arrow', colormap='jet')
+from tvtk.api import tvtk
+sc=tvtk.UnsignedCharArray()
+sc.from_array(colors)
+
+plot.mlab_source.dataset.point_data.scalars=colors
+plot.glyph.color_mode = 'color_by_scalar'
+plot.mlab_source.dataset.modified()
+
+
+#plot = mlab.quiver3d(xx, yy, zz, x_mag, y_mag, z_mag, mode='arrow',
+# colormap='jet')#py.MagData._create_directional_colormap())
+#angles = (1 - np.arctan2(y_mag, x_mag)/np.pi) / 2
+#plot.mlab_source.dataset.point_data.scalars = angles
+#plot.glyph.color_mode = 'color_by_scalar'
+#plot.glyph.glyph_source.glyph_position = 'center'
+#plot.module_manager.vector_lut_manager.data_range = np.array([0, limit])
+mlab.outline(plot)
+mlab.axes(plot)
+mlab.title('Title', height=0.95, size=0.35)
+#mlab.colorbar(label_fmt='%.2f')
+#mlab.colorbar(orientation='vertical')
+mlab.orientation_axes()
+mlab.show_pipeline()
+
+
+
+## magic to modify lookup table
+#plot.module_manager.scalar_lut_manager.lut._vtk_obj.SetTableRange(0, rgb_lut.shape[0])
+#plot.module_manager.scalar_lut_manager.lut.number_of_colors = rgb_lut.shape[0]
+#plot.module_manager.scalar_lut_manager.lut.table = rgb_lut
+
+
+#color = np.linspace(0, 256, np.prod(dim)).reshape(dim)
+#colors = np.asarray((color, color, color))
+
+#nr_points = 6
+#x,y,z=np.random.random((3,nr_points)) #some data
+#colors=np.random.randint(256,size=(nr_points,3)) #some RGB or RGBA colors
+#
+#pts=mlab.points3d(x,y,z)
+#sc=tvtk.UnsignedCharArray()
+#sc.from_array(colors)
+#
+#
+#pts.mlab_source.dataset.point_data.scalars=sc
+#pts.mlab_source.dataset.modified()
\ No newline at end of file
diff --git a/tests/test_compliance.py b/tests/test_compliance.py
index 077de6454f05f6504fb45892006f952fc0227488..aa94d2f1aa9fef0d7b277e12b3925f5909f18ed9 100644
--- a/tests/test_compliance.py
+++ b/tests/test_compliance.py
@@ -65,7 +65,7 @@ class TestCaseCompliance(unittest.TestCase):
sys.stdout = stdout_buffer
error_message = 'Found {} PEP8 violations!'.format(result.total_errors)
if todo_count > 0:
- error_message += 'Found {} TODOs!'.format(todo_count)
+ error_message += ' Found {} TODOs!'.format(todo_count)
self.assertEqual(result.total_errors, 0, error_message)
diff --git a/tests/test_costfunction.py b/tests/test_costfunction.py
index 31034bf39da6c1e8fcb295a36708715f3a2db051..a733e05f3ab0e9e0ab713afa11bd114ea22b4c9e 100644
--- a/tests/test_costfunction.py
+++ b/tests/test_costfunction.py
@@ -9,6 +9,7 @@ import numpy as np
from numpy.testing import assert_allclose
from pyramid.costfunction import Costfunction
+from pyramid.forwardmodel import ForwardModel
from pyramid.dataset import DataSet
from pyramid.projector import SimpleProjector
from pyramid.phasemap import PhaseMap
@@ -29,7 +30,7 @@ class TestCaseCostfunction(unittest.TestCase):
self.data.append(self.phase_map, self.projector)
self.data.append(self.phase_map, self.projector)
self.reg = FirstOrderRegularisator(self.mask, lam=1E-4)
- self.cost = Costfunction(self.data, self.reg)
+ self.cost = Costfunction(ForwardModel(self.data), self.reg)
def tearDown(self):
self.path = None
diff --git a/tests/test_regularisator.py b/tests/test_regularisator.py
index 1a5040a7752e976815429be34575c2ea5aecd72c..39ef31c7dcc5431ab8aa184e4889e2549ec2aad5 100644
--- a/tests/test_regularisator.py
+++ b/tests/test_regularisator.py
@@ -129,8 +129,9 @@ class TestCaseFirstOrderRegularisator(unittest.TestCase):
def test_hess_diag(self):
hess_diag = self.reg.hess_diag(np.ones(self.n))
- hess_diag_ref = np.diag(np.load(os.path.join(self.path, 'first_order_jac_ref.npy')))
- print hess_diag_ref
+ hess_diag_ref = np.zeros(3*self.n) # derivatives in all directions!
+ first_order_jac_ref = np.load(os.path.join(self.path, 'first_order_jac_ref.npy'))
+ hess_diag_ref[0:self.n] = np.diag(first_order_jac_ref)
assert_allclose(hess_diag, hess_diag_ref, atol=1E-7,
err_msg='Unexpected behaviour in hess_diag()!')