Implementation of missing TestCases and minor changes in other parts.

tests: implemented the missing TestCases!
dataset: mask is now created if not given!
phasemapper: guaranteed np.float32 format for RDRC (expected by Cython)
projector: now uses fft.FLOAT

pyramid/costfunction.py

@@ -89,6 +89,18 @@ class Costfunction(object):
         return self.chisq
 
     def init(self, x):
+        '''Initialise the costfunction by calculating the different cost terms.
+
+        Parameters
+        ----------
+        x : :class:`~numpy.ndarray` (N=1)
+            Vectorized magnetization distribution, for which the cost is calculated.
+
+        Returns
+        -------
+        None
+
+        '''
         self(x)
 
     def jac(self, x):

pyramid/dataset.py

@@ -91,15 +91,15 @@ class DataSet(object):
         assert a >= 0, 'Grid spacing has to be a positive number!'
         assert isinstance(dim, tuple) and len(dim) == 3, \
             'Dimension has to be a tuple of length 3!'
-        if mask is not None:
-            assert mask.shape == dim, 'Mask dimensions must match!'
-            self.n = 3 * np.sum(mask)
+        if mask is None:
+            self.mask = np.ones(dim, dtype=bool)
         else:
-            self.n = 3 * np.prod(dim)
+            assert mask.shape == dim, 'Mask dimensions must match!'
+            self.mask = mask
+        self.n = 3 * np.sum(self.mask)
         self.a = a
         self.dim = dim
         self.b_0 = b_0
-        self.mask = mask
         self.Se_inv = Se_inv
         self.phase_maps = []
         self.projectors = []
@@ -182,9 +182,13 @@ class DataSet(object):
             None
 
         '''
-        cov_list = [sparse.diags(m.flatten(), 0) for m in mask_list]
+        cov_list = [sparse.diags(m.flatten().astype(np.float32), 0) for m in mask_list]
         self.set_Se_inv_block_diag(cov_list)
 
+    def create_3d_mask(self, mask_list):
+        # TODO: method for constructing 3D mask from 2D masks?
+        raise NotImplementedError()
+
     def display_phase(self, mag_data=None, title='Phase Map',
                       cmap='RdBu', limit=None, norm=None):
         '''Display all phasemaps saved in the :class:`~.DataSet` as a colormesh.
@@ -223,7 +227,7 @@ class DataSet(object):
         plt.show()
 
     def display_combined(self, mag_data=None, title='Combined Plot', cmap='RdBu', limit=None,
-                         norm=None, gain=1, interpolation='none', grad_encode='bright'):
+                         norm=None, gain='auto', interpolation='none', grad_encode='bright'):
         '''Display all phasemaps and the resulting color coded holography images.
 
         Parameters
@@ -268,6 +272,3 @@ class DataSet(object):
                                     cmap, limit, norm, gain, interpolation, grad_encode)
             for (i, phase_map) in enumerate(phase_maps)]
         plt.show()
-
-
-# TODO: method for constructing 3D mask from 2D masks?

pyramid/diagnostics.py

@@ -180,7 +180,6 @@ class Diagnostics(object):
 
         '''# TODO: Docstring!
         mag_data_avg_kern = self.get_avg_kern_row(pos)
-        print self.pos
         mag_x, mag_y, mag_z = mag_data_avg_kern.magnitude
         x = mag_x.sum(axis=(0, 1))
         y = mag_y.sum(axis=(0, 2))

pyramid/fft.py

@@ -149,6 +149,7 @@ def _irfftn_fftw(a, s=None, axes=None):
 
 
 def _rfftn_adj_fftw(a):
+    # TODO: Careful just works for even a (which is guaranteed by the kernel!)
     n = 2 * (a.shape[-1] - 1)
     out_shape = a.shape[:-1] + (n,)
     out_arr = zeros(out_shape, dtype=a.dtype)

pyramid/forwardmodel.py

@@ -66,7 +66,7 @@ class ForwardModel(object):
 
     def __call__(self, x):
         self._log.debug('Calling __call__')
-        self.mag_data.magnitude[:] = 0
+        self.mag_data.magnitude[...] = 0
         self.mag_data.set_vector(x, self.data_set.mask)
         result = np.zeros(self.m)
         hp = self.hook_points
@@ -102,7 +102,7 @@ class ForwardModel(object):
         for i, projector in enumerate(self.data_set.projectors):
             mag_vec = self.mag_data.mag_vec
             res = self.phase_mappers[projector.dim_uv].jac_dot(projector.jac_dot(mag_vec))
-            result[hp[i]:hp[i+1]] = res.flatten()
+            result[hp[i]:hp[i+1]] = res
         return result
 
     def jac_T_dot(self, x, vector):

pyramid/magdata.py

@@ -344,7 +344,7 @@ class MagData(object):
         -------
         None
 
-        '''
+        '''  # TODO: Implement default folder for all mag_data and phasemaps!!!
         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:

pyramid/phasemapper.py

@@ -285,8 +285,10 @@ class PhaseMapperRDRC(PhaseMapper):
         assert len(vector) == self.n, \
             'vector size not compatible! vector: {}, size: {}'.format(len(vector), self.n)
         v_dim, u_dim = self.kernel.dim_uv
-        result = np.zeros(self.m)
+        result = np.zeros(self.m, dtype=np.float32)
         if self.numcore:
+            if vector.dtype != np.float32:
+                vector = vector.astype(np.float32)
             nc.jac_dot_real_convolve(v_dim, u_dim, self.kernel.u, self.kernel.v, vector, result)
         else:
             # Iterate over all contributing pixels (numbered consecutively)
@@ -319,9 +321,11 @@ class PhaseMapperRDRC(PhaseMapper):
         '''
         assert len(vector) == self.m, \
             'vector size not compatible! vector: {}, size: {}'.format(len(vector), self.m)
-        v_dim, u_dim = self.dim_uv
-        result = np.zeros(self.n)
+        v_dim, u_dim = self.kernel.dim_uv
+        result = np.zeros(self.n, dtype=np.float32)
         if self.numcore:
+            if vector.dtype != np.float32:
+                vector = vector.astype(np.float32)
             nc.jac_T_dot_real_convolve(v_dim, u_dim, self.kernel.u, self.kernel.v, vector, result)
         else:
             # Iterate over all contributing pixels (numbered consecutively):
@@ -431,9 +435,9 @@ class PhaseMapperFDFC(PhaseMapper):
         '''
         assert len(vector) == self.n, \
             'vector size not compatible! vector: {}, size: {}'.format(len(vector), self.n)
-        mag_proj = MagData(self.a, np.zeros((3, 1)+self.dim_uv))
-        magnitude_proj = self.jac_dot(vector).reshape((2, )+self.dim_uv)
-        mag_proj.magnitude[1:3, 0, ...] = magnitude_proj
+        mag_proj = MagData(self.a, np.zeros((3, 1)+self.dim_uv, dtype=np.float32))
+        magnitude_proj = np.reshape(vector, (2,)+self.dim_uv)
+        mag_proj.magnitude[:2, 0, ...] = magnitude_proj
         return self(mag_proj).phase_vec
 
     def jac_T_dot(self, vector):
@@ -494,6 +498,7 @@ class PhaseMapperElectric(PhaseMapper):
             (self.a, self.dim_uv, self.v_0, self.v_acc, self.threshold)
 
     def __call__(self, mag_data):
+        raise NotImplementedError()  # TODO: Implement right!
         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!'

pyramid/projector.py

@@ -16,6 +16,7 @@ import itertools
 from scipy.sparse import coo_matrix, csr_matrix
 
 from pyramid.magdata import MagData
+import pyramid.fft as fft
 
 import logging
 
@@ -84,15 +85,16 @@ class Projector(object):
 
     def __call__(self, mag_data):
         self._log.debug('Calling __call__')
-        mag_proj = MagData(mag_data.a, np.zeros((3, 1)+self.dim_uv))
+        mag_proj = MagData(mag_data.a, fft.zeros((3, 1)+self.dim_uv, dtype=fft.FLOAT))
         magnitude_proj = self.jac_dot(mag_data.mag_vec).reshape((2, )+self.dim_uv)
         mag_proj.magnitude[0:2, 0, ...] = magnitude_proj
         return mag_proj
 
     def _vector_field_projection(self, vector):
         size_2d, size_3d = self.size_2d, self.size_3d
-        result = np.zeros(2*size_2d)
+        result = fft.zeros(2*size_2d, dtype=fft.FLOAT)
         # Go over all possible component projections (z, y, x) to (u, v):
+        # TODO: save self.weight.dot(...)
         if self.coeff[0][0] != 0:  # x to u
             result[:size_2d] += self.coeff[0][0] * self.weight.dot(vector[:size_3d])
         if self.coeff[0][1] != 0:  # y to u
@@ -109,7 +111,7 @@ class Projector(object):
 
     def _vector_field_projection_T(self, vector):
         size_2d, size_3d = self.size_2d, self.size_3d
-        result = np.zeros(3*size_3d)
+        result = np.zeros(3*size_3d, dtype=fft.FLOAT)
         # Go over all possible component projections (u, v) to (z, y, x):
         if self.coeff[0][0] != 0:  # u to x
             result[:size_3d] += self.coeff[0][0] * self.weight.T.dot(vector[:size_2d])
@@ -438,6 +440,7 @@ class SimpleProjector(Projector):
     def __init__(self, dim, axis='z', dim_uv=None):
         self._log.debug('Calling __init__')
         assert axis in {'z', 'y', 'x'}, 'Projection axis has to be x, y or z (given as a string)!'
+        self.axis = axis
         proj, v, u = self.AXIS_DICT[axis]
         if axis == 'x':
             dim_proj, dim_v, dim_u = dim[proj], dim[u], dim[v]  # coordinate switch for 'x'!

pyramid/regularisator.py

@@ -91,22 +91,30 @@ class Regularisator(object):
         Returns
         -------
         result : :class:`~numpy.ndarray` (N=1)
-            Product of the input `vector` with the Hessian matrix of the costfunction.
+            Product of the input `vector` with the Hessian matrix.
 
         '''
         return self.lam * self.norm.hess_dot(x, vector)
 
-    def hess_diag(self, x, vector):
+    def hess_diag(self, x):
         ''' Return the diagonal of the Hessian.
 
         Parameters
         ----------
-        _ : undefined
-            Unused input
+        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.
+
+        Returns
+        -------
+        result : :class:`~numpy.ndarray` (N=1)
+            Diagonal of the Hessian matrix.
 
         '''
         self._log.debug('Calling hess_diag')
-        return self.lam * self.norm.hess_diag(x, vector)
+        return self.lam * self.norm.hess_diag(x)
 
 
 class NoneRegularisator(Regularisator):
@@ -159,10 +167,10 @@ class NoneRegularisator(Regularisator):
         ----------
         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
+            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.
-        vector : :class:`~numpy.ndarray` (N=1)
+         vector : :class:`~numpy.ndarray` (N=1)
             Vectorized magnetization distribution which is multiplied by the Hessian.
 
         Returns
@@ -173,17 +181,25 @@ class NoneRegularisator(Regularisator):
         '''
         return np.zeros_like(vector)
 
-    def hess_diag(self, x, vector):
+    def hess_diag(self, x):
         ''' Return the diagonal of the Hessian.
 
         Parameters
         ----------
-        _ : undefined
-            Unused input
+        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.
+
+        Returns
+        -------
+        result : :class:`~numpy.ndarray` (N=1)
+            Diagonal of the Hessian matrix.
 
         '''
         self._log.debug('Calling hess_diag')
-        return np.zeros_like(vector)
+        return np.zeros_like(x)
 
 
 class ZeroOrderRegularisator(Regularisator):
@@ -204,7 +220,7 @@ class ZeroOrderRegularisator(Regularisator):
 
     _log = logging.getLogger(__name__+'.ZeroOrderRegularisator')
 
-    def __init__(self, _, lam, p=2):
+    def __init__(self, _=None, lam=1E-4, p=2):
         self._log.debug('Calling __init__')
         self.p = p
         if p == 2:
@@ -234,7 +250,7 @@ class FirstOrderRegularisator(Regularisator):
 
     '''
 
-    def __init__(self, mask, lam, p=2):
+    def __init__(self, mask, lam=1E-4, p=2):
         self.p = p
         D0 = jdiff.get_diff_operator(mask, 0, 3)
         D1 = jdiff.get_diff_operator(mask, 1, 3)

pyramid/version.py

 # THIS FILE IS GENERATED FROM THE PYRAMID SETUP.PY
 version = "0.1.0-dev"
-hg_revision = "51a00f1d3934+"
+hg_revision = "4fdb98485ce0+"

scripts/collaborations/patrick_nanowire_reconstruction.py

@@ -24,13 +24,13 @@ logging.config.fileConfig(LOGGING_CONF, disable_existing_loggers=False)
 
 ###################################################################################################
 a = 1.455  # in nm
-gain = 5
+gain = 50
 b_0 = 1
 lam = 1E-4
 PATH = '../../output/patrick/'
-PHASE = 'pos3_40deg_magphase'
-MASK = 'pos3_40deg_maskbyhand'
-FORMAT = '.bmp'
+PHASE = 'Reza_30_uj_tube_M'
+MASK = 'Reza_30_uj_tube_maskbygimp'
+FORMAT = '.tif'
 longFOV = False
 longFOV_string = np.where(longFOV, 'longFOV', 'normalFOV')
 IMAGENAME = '{}_{}_{}_'.format(MASK, PHASE, longFOV_string)
@@ -70,19 +70,19 @@ if longFOV:
 regularisator = FirstOrderRegularisator(mask, lam, p=2)
 
 with TakeTime('reconstruction time'):
-    mag_data_rec = rc.optimize_linear(data_set, regularisator=regularisator, max_iter=1000)[0]
+    mag_data_rec = rc.optimize_linear(data_set, regularisator=regularisator, max_iter=500)[0]
 
 phase_map_rec_pad = pm(mag_data_rec)
 phase_map_rec = PhaseMap(a, phase_map_rec_pad.phase[pad:, :])#[pad:-pad, pad:-pad])
 phase_map_diff = phase_map_rec - phase_map
 
 # Display the reconstructed phase map and holography image:
-phase_map.display_combined('Input PhaseMap', gain=4)
+phase_map.display_combined('Input PhaseMap', gain=gain)
 plt.savefig(PATH+IMAGENAME+'ORIGINAL.png')
-phase_map_pad.display_combined('Input PhaseMap (padded)', gain=4)
-phase_map_rec_pad.display_combined('Reconstr. Distribution (padded)', gain=4)
+phase_map_pad.display_combined('Input PhaseMap (padded)', gain=gain)
+phase_map_rec_pad.display_combined('Reconstr. Distribution (padded)', gain=gain)
 plt.savefig(PATH+IMAGENAME+'RECONSTRUCTION_PADDED.png')
-phase_map_rec.display_combined('Reconstr. Distribution', gain=4)
+phase_map_rec.display_combined('Reconstr. Distribution', gain=gain)
 plt.savefig(PATH+IMAGENAME+'RECONSTRUCTION.png')
 phase_map_diff.display_combined('Difference')
 plt.savefig(PATH+IMAGENAME+'DIFFERENCE.png')

scripts/collaborations/rueffner_file.py

@@ -29,6 +29,7 @@ LOGGING_CONF = os.path.join(os.path.dirname(os.path.realpath(pyramid.__file__)),
 
 logging.config.fileConfig(LOGGING_CONF, disable_existing_loggers=False)
 proc = psutil.Process(os.getpid())
+
 ###################################################################################################
 PATH = '../../output/'
 dim = (16, 190, 220)
@@ -39,7 +40,9 @@ dens_z = 4E3
 dens_x = 1E2
 lim_z = 22
 lim_x = 0.35
+
 ###################################################################################################
+
 # Build projectors and phasemapper:
 projector_z = SimpleProjector(dim, axis='z', dim_uv=dim_uv)
 projector_x = SimpleProjector(dim, axis='x', dim_uv=dim_uv)

scripts/collaborations/vtk_conversion.py

@@ -12,6 +12,9 @@ import pickle
 import vtk
 from tqdm import tqdm
 
+import psutil
+import gc
+
 import pyramid
 from pyramid.magdata import MagData
 from pyramid.projector import XTiltProjector
@@ -26,6 +29,7 @@ LOGGING_CONF = os.path.join(os.path.dirname(os.path.realpath(pyramid.__file__)),
 
 
 logging.config.fileConfig(LOGGING_CONF, disable_existing_loggers=False)
+proc = psutil.Process(os.getpid())
 ###################################################################################################
 PATH = '../../output/vtk data/tube_160x30x1100nm/02758'
 b_0 = 1.54
@@ -158,31 +162,37 @@ mag_data.quiver_plot()
 
 dim = mag_data.dim
 dim_uv = (500, 200)
-angles = [0, 20, 40, 60]
-
-mag_data_xy = mag_data.copy()
-mag_data_xy.magnitude[2] = 0
+angles = np.arange(-60, 61, 5)#[0, 20, 40, 60]
 
-mag_data_z = mag_data.copy()
-mag_data_z.magnitude[0] = 0
-mag_data_z.magnitude[1] = 0
+#mag_data_xy = mag_data.copy()
+#mag_data_xy.magnitude[2] = 0
+#
+#mag_data_z = mag_data.copy()
+#mag_data_z.magnitude[0] = 0
+#mag_data_z.magnitude[1] = 0
 
 # Iterate over all angles:
 for angle in angles:
     angle_rad = np.pi/2 + angle*np.pi/180
     projector = XTiltProjector(dim, angle_rad, dim_uv)
-    mag_proj = projector(mag_data_z)
+    mag_proj = projector(mag_data)
     phase_map = PhaseMapperRDFC(Kernel(mag_data.a, projector.dim_uv))(mag_proj)
-    phase_map.display_combined('Phase Map Nanowire Tip', gain=gain,
-                               interpolation='bilinear')
-    plt.savefig(PATH+'_nanowire_z_xtilt_{}.png'.format(angle), dpi=500)
-    mag_proj.scale_down(2)
-    axis = mag_proj.quiver_plot()
-    plt.savefig(PATH+'_nanowire_z_mag_xtilt_{}.png'.format(angle), dpi=500)
-    axis = mag_proj.quiver_plot(log=True)
-    plt.savefig(PATH+'_nanowire_z_mag_log_xtilt_{}.png'.format(angle), dpi=500)
+    phase_map.display_phase('Phase Map Nanowire Tip', cmap='gray')
+    plt.savefig(PATH+'_nanowire_xtilt_{}.png'.format(angle), dpi=500)
+#    phase_map = PhaseMapperRDFC(Kernel(mag_data.a, projector.dim_uv))(mag_proj)
+#    phase_map.display_combined('Phase Map Nanowire Tip', gain=gain,
+#                               interpolation='bilinear')
+#    plt.savefig(PATH+'_nanowire_xtilt_{}.png'.format(angle), dpi=500)
+#    mag_proj.scale_down(2)
+#    axis = mag_proj.quiver_plot()
+#    plt.savefig(PATH+'_nanowire_mag_xtilt_{}.png'.format(angle), dpi=500)
+#    axis = mag_proj.quiver_plot(log=True)
+#    plt.savefig(PATH+'_nanowire_mag_log_xtilt_{}.png'.format(angle), dpi=500)
     # Close plots:
     plt.close('all')
+    gc.collect()
+    print 'RSS = {:.2f} MB'.format(proc.memory_info().rss/1024.**2)
+    print angle, 'deg. done!'
 
 
 #mag_data.scale_down(2)

scripts/publications/paper 2/paper2.py

+# -*- coding: utf-8 -*-
+"""
+Created on Wed Feb 11 10:57:53 2015
+
+@author: Jan
+"""
+
+
+import numpy as np
+
+from pyramid import *  # analysis:ignore
+from pyramid import magcreator as mc
+from pyramid import reconstruction as rc
+
+from jutil.taketime import TakeTime
+
+from matplotlib.patches import Rectangle
+
+from PIL import Image
+
+import pickle
+
+import logging
+import logging.config
+
+
+LOGGING_CONF = os.path.join(os.path.dirname(os.path.realpath(pyramid.__file__)), 'logging.ini')
+
+
+logging.config.fileConfig(LOGGING_CONF, disable_existing_loggers=False)
+
+
+def calc_diagnostics(pos, mag_data_opt, cost, max_iter=2000):
+    diag = Diagnostics(mag_data_opt.mag_vec, cost, max_iter=2000)
+    for c in (0, 1, 2):
+        print 'calc. diagnostics for {}-comp. at pos: {}'.format({0: 'x', 1: 'y', 2: 'z'}[c], pos)
+        diag.pos = (c,) + pos
+        gain_maps = diag.get_gain_row_maps()
+        axis, cbar = gain_maps[0].display_phase('Gain map for pos: {}'.format(pos))
+        axis.add_patch(Rectangle((diag.pos[3], diag.pos[2]), 1, 1, linewidth=2,
+                                 color='g', fill=False))
+        cbar.set_label(u'magnetization/phase [1/rad]', fontsize=15)
+        diag.get_avg_kern_row().quiver_plot3d()
+        mcon = diag.measure_contribution
+        print 'measurement contr. (min - max): {:.2f} - {:.2f}'.format(mcon.min(), mcon.max())
+        px_avrg, fwhm, (x, y, z) = diag.calculate_averaging()
+        print 'px_avrg:', px_avrg
+        print 'fwhm:', fwhm
+
+
+diag_pos = [(0, 31, 31), (0, 16, 31)]
+
+
+## SIMULATIONS:
+#
+#a = 10.
+#b_0 = 1.
+#lam = 1E-4
+#dim = (1, 64, 64)
+#center = (0.5, dim[1]//2-0.5, dim[2]//2-0.5)
+#r = dim[1]//4
+#
+#
+## SIMULATED VORTEX:
+#
+#mag_shape = mc.Shapes.disc(dim, center, radius=r, height=2)
+#mag_vort = MagData(a, mc.create_mag_dist_vortex(mag_shape))
+#mag_vort.quiver_plot('Original distribution')
+#pm(mag_vort).display_combined('Original distribution')
+#
+## without mask:
+#data_vort_nomask = DataSet(a, dim, b_0, mask=None)
+#data_vort_nomask.projectors = [SimpleProjector(dim)]
+#data_vort_nomask.phase_maps = data_vort_nomask.create_phase_maps(mag_vort)
+#reg = FirstOrderRegularisator(data_vort_nomask.mask, lam, p=2)
+#with TakeTime('reconstruction'):
+#    mag_vort_nomask, cost_vort_nomask = rc.optimize_linear(data_vort_nomask, reg, max_iter=100)
+#mag_vort_nomask.quiver_plot('Reconstructed distribution (without mask)')
+#pm(mag_vort_nomask).display_combined('Reconstructed distribution (without mask)')
+#pm(mag_vort_nomask-mag_vort).display_phase('Difference (without mask)')
+##for pos in diag_pos:
+##    calc_diagnostics(pos, mag_vort_nomask, cost_vort_nomask)
+#
+## with mask:
+#data_vort_mask = DataSet(a, dim, b_0, mask=mag_vort.get_mask())
+#data_vort_mask.projectors = [SimpleProjector(dim)]
+#data_vort_mask.phase_maps = data_vort_mask.create_phase_maps(mag_vort)
+#reg = FirstOrderRegularisator(data_vort_mask.mask, lam, p=2)
+#with TakeTime('reconstruction'):
+#    mag_vort_mask, cost_vort_mask = rc.optimize_linear(data_vort_mask, reg, max_iter=100)
+#mag_vort_mask.quiver_plot('Reconstructed distribution (with mask)')
+#pm(mag_vort_mask).display_combined('Reconstructed distribution (with mask)')
+#pm(mag_vort_mask-mag_vort).display_phase('Difference (with mask)')
+#
+#
+## SIMULATED VORTEX:
+#
+#mag_shape = magcreator.Shapes.slab(dim, center, width=(2, 2*r, 2*r))
+#mag_slab = MagData(a, magcreator.create_mag_dist_homog(mag_shape, phi=np.pi/4))
+#mag_slab.quiver_plot('Original distribution')
+#pm(mag_slab).display_combined('Original distribution')
+#
+## without mask:
+#data_slab_nomask = DataSet(a, dim, b_0, mask=None)
+#data_slab_nomask.projectors = [SimpleProjector(dim)]
+#data_slab_nomask.phase_maps = data_slab_nomask.create_phase_maps(mag_slab)
+#reg = FirstOrderRegularisator(data_slab_nomask.mask, lam, p=2)
+#with TakeTime('reconstruction'):
+#    mag_slab_nomask, cost_slab_nomask = rc.optimize_linear(data_slab_nomask, reg, max_iter=100)
+#mag_slab_nomask.quiver_plot('Reconstructed distribution (without mask)')
+#pm(mag_slab_nomask).display_combined('Reconstructed distribution (without mask)')
+#pm(mag_slab_nomask-mag_slab).display_phase('Difference (without mask)')
+#
+## with mask:
+#data_slab_mask = DataSet(a, dim, b_0, mask=mag_slab.get_mask())
+#data_slab_mask.projectors = [SimpleProjector(dim)]
+#data_slab_mask.phase_maps = data_slab_mask.create_phase_maps(mag_slab)
+#reg = FirstOrderRegularisator(data_slab_mask.mask, lam, p=2)
+#with TakeTime('reconstruction'):
+#    mag_slab_mask, cost_slab_mask = rc.optimize_linear(data_slab_mask, reg, max_iter=100)
+#mag_slab_mask.quiver_plot('Reconstructed distribution (with mask)')
+#pm(mag_slab_mask).display_combined('Reconstructed distribution (with mask)')
+#pm(mag_slab_mask-mag_slab).display_phase('Difference (with mask)')
+
+
+# MAGNETITE CUBES
+
+a = 1.0  # in nm
+b_0 = 1
+lam = 1E-4
+
+# 2 particles:
+phase_map_2 = PhaseMap.load_from_netcdf4('../../../output/joern/phase_map_2.nc')
+phase_map_2.display_combined('Original distribution')
+with open('../../../output/joern/mask_2.pickle', 'rb') as pf:
+    mask_2 = pickle.load(pf)
+dim_2 = mask_2.shape
+data_2 = DataSet(a, dim_2, b_0, mask=mask_2)
+data_2.projectors = [SimpleProjector(dim_2)]
+data_2.phase_maps = [phase_map_2]
+reg = FirstOrderRegularisator(mask_2, lam, p=2)
+with TakeTime('reconstruction'):
+    mag_2, cost_2 = rc.optimize_linear(data_2, reg, max_iter=100)
+mag_2.quiver_plot('Reconstructed distribution')
+pm(mag_2).display_combined('Reconstructed distribution')
+(pm(mag_2)-phase_map_2).display_phase('Difference')
+
+
+# 4 particles:
+phase_map_4 = PhaseMap.load_from_netcdf4('../../../output/joern/phase_map_4.nc')
+phase_map_4.display_combined('Original distribution')
+with open('../../../output/joern/mask_4.pickle', 'rb') as pf:
+    mask_4 = pickle.load(pf)
+dim_4 = mask_4.shape
+data_4 = DataSet(a, dim_4, b_0, mask=mask_4)
+data_4.projectors = [SimpleProjector(dim_4)]
+data_4.phase_maps = [phase_map_4]
+reg = FirstOrderRegularisator(mask_4, lam, p=2)
+with TakeTime('reconstruction'):
+    mag_4, cost_4 = rc.optimize_linear(data_4, reg, max_iter=100)
+mag_4.quiver_plot('Reconstructed distribution')
+pm(mag_4).display_combined('Reconstructed distribution')
+(pm(mag_4)-phase_map_4).display_phase('Difference')
+
+
+
+# EXPERIMENTAL NANOWIRE:
+a = 1.455  # in nm
+gain = 50
+b_0 = 1
+lam = 1E-4
+PATH = '../../../output/patrick/'
+PHASE = 'Reza_30_uj_tube_M'
+MASK = 'Reza_30_uj_tube_maskbygimp'
+FORMAT = '.tif'
+longFOV = False
+longFOV_string = np.where(longFOV, 'longFOV', 'normalFOV')
+IMAGENAME = '{}_{}_{}_'.format(MASK, PHASE, longFOV_string)
+PHASE_MAX = 0.5  # -10°: 0.92, 30°: 7.68
+PHASE_MIN = -0.5  # -10°: -16.85, 30°: -18.89
+PHASE_DELTA = PHASE_MAX - PHASE_MIN
+###################################################################################################
+
+# Read in files:
+im_mask = Image.open(PATH+MASK+FORMAT)
+#im_mask.thumbnail((125, 175), Image.ANTIALIAS)
+im_phase = Image.open(PATH+PHASE+FORMAT)
+#im_phase.thumbnail((125, 125), Image.ANTIALIAS)
+
+mask = np.array(np.array(im_mask)/255, dtype=bool)
+dim_uv = mask.shape
+phase = (np.array(im_phase)/255.-0.5) * PHASE_DELTA
+pad = dim_uv[0] - phase.shape[0]#25
+phase_pad = np.zeros(dim_uv)
+phase_pad[pad:, :] = phase#[pad:-pad, pad:-pad] = phase
+
+mask = np.expand_dims(mask, axis=0)
+dim = mask.shape
+
+phase_map = PhaseMap(a, phase)
+phase_map.display_combined('Original distribution')
+phase_map_pad = PhaseMap(a, phase_pad)
+
+data_set = DataSet(a, dim, b_0, mask)
+data_set.append(phase_map_pad, SimpleProjector(dim))
+
+reg = FirstOrderRegularisator(mask, lam, p=2)
+with TakeTime('reconstruction'):
+    mag_data_rec, cost_rec = rc.optimize_linear(data_set, reg, max_iter=100)
+mag_data_rec.quiver_plot('Reconstructed distribution')
+pm(mag_data_rec).display_combined('Reconstructed distribution')
+(pm(mag_data_rec)-phase_map).display_phase('Difference')

scripts/reconstruction/reconstruction_linear_test.py

@@ -32,7 +32,7 @@ print('--Generating input phase_maps')
 
 a = 10.
 b_0 = 1.
-dim = (32, 32, 32)
+dim = (64, 64, 64)
 dim_uv = dim[1:3]
 count = 16
 lam = 1E-4
@@ -90,45 +90,44 @@ with TakeTime('reconstruction'):
 print 'Cost:', cost.chisq
 
 ###################################################################################################
-print('--Plot stuff')
-
-limit = 1.2
-mag_data.quiver_plot3d('Original distribution', limit=limit)
-mag_data_opt.quiver_plot3d('Reconstructed distribution', limit=limit)
-(mag_data_opt - mag_data).quiver_plot3d('Difference')
-phase_maps_opt = data.create_phase_maps(mag_data_opt)
-
-
-from pyramid.diagnostics import Diagnostics
-from matplotlib.patches import Rectangle
-
-
-diag = Diagnostics(mag_data_opt.mag_vec, cost, max_iter=2000)
-
-print 'position:', diag.pos
-print 'std:', diag.std
-gain_maps = diag.get_gain_row_maps()
-axis, cbar = gain_maps[count//2].display_phase()
-axis.add_patch(Rectangle((diag.pos[3], diag.pos[2]), 1, 1, linewidth=2, color='g', fill=False))
-cbar.set_label(u'magnetization/phase [1/rad]', fontsize=15)
-diag.get_avg_kern_row().quiver_plot3d()
-mcon = diag.measure_contribution
-print 'measurement contr. (min - max): {:.2f} - {:.2f}'.format(mcon.min(), mcon.max())
-px_avrg, fwhm, (x, y, z) = diag.calculate_averaging()
-print 'px_avrg:', px_avrg
-print 'fwhm:', fwhm
-
-diag.pos = (1, dim[0]//2, dim[1]//2, dim[2]//2)
-
-print 'position:', diag.pos
-print 'std:', diag.std
-gain_maps = diag.get_gain_row_maps()
-axis, cbar = gain_maps[count//2].display_phase()
-axis.add_patch(Rectangle((diag.pos[3], diag.pos[2]), 1, 1, linewidth=2, color='g', fill=False))
-cbar.set_label(u'magnetization/phase [1/rad]', fontsize=15)
-diag.get_avg_kern_row().quiver_plot3d()
-mcon = diag.measure_contribution
-print 'measurement contr. (min - max): {:.2f} - {:.2f}'.format(mcon.min(), mcon.max())
-px_avrg, fwhm, (x, y, z) = diag.calculate_averaging()
-print 'px_avrg:', px_avrg
-print 'fwhm:', fwhm
+#print('--Plot stuff')
+#
+#limit = 1.2
+#mag_data.quiver_plot3d('Original distribution', limit=limit)
+#mag_data_opt.quiver_plot3d('Reconstructed distribution', limit=limit)
+#(mag_data_opt - mag_data).quiver_plot3d('Difference')
+#phase_maps_opt = data.create_phase_maps(mag_data_opt)
+#
+#
+#from pyramid.diagnostics import Diagnostics
+#from matplotlib.patches import Rectangle
+#
+#
+#diag = Diagnostics(mag_data_opt.mag_vec, cost, max_iter=2000)
+#
+#print 'position:', diag.pos
+#print 'std:', diag.std
+#gain_maps = diag.get_gain_row_maps()
+#axis, cbar = gain_maps[count//2].display_phase()
+#axis.add_patch(Rectangle((diag.pos[3], diag.pos[2]), 1, 1, linewidth=2, color='g', fill=False))
+#cbar.set_label(u'magnetization/phase [1/rad]', fontsize=15)
+#diag.get_avg_kern_row().quiver_plot3d()
+#mcon = diag.measure_contribution
+#print 'measurement contr. (min - max): {:.2f} - {:.2f}'.format(mcon.min(), mcon.max())
+#px_avrg, fwhm, (x, y, z) = diag.calculate_averaging()
+#print 'px_avrg:', px_avrg
+#print 'fwhm:', fwhm
+#
+#diag.pos = (1, dim[0]//2, dim[1]//2, dim[2]//2)
+#print 'position:', diag.pos
+#print 'std:', diag.std
+#gain_maps = diag.get_gain_row_maps()
+#axis, cbar = gain_maps[count//2].display_phase()
+#axis.add_patch(Rectangle((diag.pos[3], diag.pos[2]), 1, 1, linewidth=2, color='g', fill=False))
+#cbar.set_label(u'magnetization/phase [1/rad]', fontsize=15)
+#diag.get_avg_kern_row().quiver_plot3d()
+#mcon = diag.measure_contribution
+#print 'measurement contr. (min - max): {:.2f} - {:.2f}'.format(mcon.min(), mcon.max())
+#px_avrg, fwhm, (x, y, z) = diag.calculate_averaging()
+#print 'px_avrg:', px_avrg
+#print 'fwhm:', fwhm

tests/test_compliance.py

@@ -8,7 +8,6 @@ import datetime
 import unittest
 
 import re
-
 import pep8
 
 
@@ -35,7 +34,6 @@ class TestCaseCompliance(unittest.TestCase):
         pep8.MAX_LINE_LENGTH = 99
         pep8style = pep8.StyleGuide(quiet=False)
         pep8style.options.ignore = ignores
-
         stdout_buffer = sys.stdout
         with open(os.path.join(self.path, 'output', 'pep8_log.txt'), 'w') as sys.stdout:
             print '<<< PEP8 LOGFILE >>>'
@@ -45,9 +43,9 @@ class TestCaseCompliance(unittest.TestCase):
             print '\nERRORS AND WARNINGS:'
             result = pep8style.check_files(files)
             if result.total_errors == 0:
-                print 'No Warnings or Errors detected!'
+                print 'No PEP8 violations detected!'
             else:
-                print '---->   {} Warnings and Errors detected!'.format(result.total_errors)
+                print '---->   {} PEP8 violations detected!'.format(result.total_errors)
             print '\nTODOS:'
             todos_found = False
             todo_count = 0
@@ -64,10 +62,10 @@ class TestCaseCompliance(unittest.TestCase):
                 print '---->   {} TODOs found!'.format(todo_count)
             else:
                 print 'No TODOS found!'
-
         sys.stdout = stdout_buffer
-
-        error_message = 'Found %s Errors and Warnings!' % result.total_errors
+        error_message = 'Found {} PEP8 violations!'.format(result.total_errors)
+        if todo_count > 0:
+            error_message += 'Found {} TODOs!'.format(todo_count)
         self.assertEqual(result.total_errors, 0, error_message)
 
 

tests/test_costfunction.py

 # -*- coding: utf-8 -*-
-"""Created on Mon Jan 20 20:01:25 2014 @author: Jan"""
+"""Testcase for the costfunction module"""
+
+
+import os
+import unittest
+
+import numpy as np
+from numpy.testing import assert_allclose
+
+from pyramid.costfunction import Costfunction
+from pyramid.dataset import DataSet
+from pyramid.projector import SimpleProjector
+from pyramid.phasemap import PhaseMap
+from pyramid.regularisator import FirstOrderRegularisator
+
+
+class TestCaseCostfunction(unittest.TestCase):
+
+    def setUp(self):
+        self.path = os.path.join(os.path.dirname(os.path.realpath(__file__)), 'test_costfunction')
+        self.a = 10.
+        self.dim = (4, 5, 6)
+        self.mask = np.zeros(self.dim, dtype=bool)
+        self.mask[1:-1, 1:-1, 1:-1] = True
+        self.data = DataSet(self.a, self.dim, mask=self.mask)
+        self.projector = SimpleProjector(self.dim)
+        self.phase_map = PhaseMap.load_from_netcdf4(os.path.join(self.path, 'phase_map_ref.nc'))
+        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)
+
+    def tearDown(self):
+        self.path = None
+        self.a = None
+        self.dim = None
+        self.mask = None
+        self.data = None
+        self.projector = None
+        self.phase_map = None
+        self.reg = None
+        self.cost = None
+
+    def test_call(self):
+        assert_allclose(self.cost(np.ones(self.cost.n)), 0.,
+                        err_msg='Unexpected behaviour in __call__()!')
+        zero_vec_cost = np.load(os.path.join(self.path, 'zero_vec_cost.npy'))
+        assert_allclose(self.cost(np.zeros(self.cost.n)), zero_vec_cost,
+                        err_msg='Unexpected behaviour in __call__()!')
+
+    def test_jac(self):
+        assert_allclose(self.cost.jac(np.ones(self.cost.n)), np.zeros(self.cost.n),
+                        err_msg='Unexpected behaviour in jac()!')
+        jac_vec_ref = np.load(os.path.join(self.path, 'jac_vec_ref.npy'))
+        assert_allclose(self.cost.jac(np.zeros(self.cost.n)), jac_vec_ref, atol=1E-7,
+                        err_msg='Unexpected behaviour in jac()!')
+        jac = np.array([self.cost.jac(np.eye(self.cost.n)[:, i]) for i in range(self.cost.n)]).T
+        jac_ref = np.load(os.path.join(self.path, 'jac_ref.npy'))
+        assert_allclose(jac, jac_ref, atol=1E-7,
+                        err_msg='Unexpected behaviour in jac()!')
+
+    def test_hess_dot(self):
+        assert_allclose(self.cost.hess_dot(None, np.zeros(self.cost.n)), np.zeros(self.cost.n),
+                        err_msg='Unexpected behaviour in jac()!')
+        hess_vec_ref = -np.load(os.path.join(self.path, 'jac_vec_ref.npy'))
+        assert_allclose(self.cost.hess_dot(None, np.ones(self.cost.n)), hess_vec_ref, atol=1E-7,
+                        err_msg='Unexpected behaviour in jac()!')
+        hess = np.array([self.cost.hess_dot(None, np.eye(self.cost.n)[:, i])
+                        for i in range(self.cost.n)]).T
+        hess_ref = np.load(os.path.join(self.path, 'hess_ref.npy'))
+        assert_allclose(hess, hess_ref, atol=1E-7,
+                        err_msg='Unexpected behaviour in hess_dot()!')
+
+    def test_hess_diag(self):
+        assert_allclose(self.cost.hess_diag(None), np.ones(self.cost.n),
+                        err_msg='Unexpected behaviour in hess_diag()!')
+
+
+if __name__ == '__main__':
+    suite = unittest.TestLoader().loadTestsFromTestCase(TestCaseCostfunction)
+    unittest.TextTestRunner(verbosity=2).run(suite)