Implementation of the Regularisator class

regularisator: new module to represent different regularisation strategies
costfunction: changed to use the new Regularisator classes, changed arguments!
dataset: fixed bug that always displayed the own phase maps instead of provided
reconstruction: changed to use the new Regularisator classes, changed arguments!
scripts/rueffner_file: tried to avoid Memory Error (no success)
scripts/zi_an: added holographic contour plots with overlayd magnetization
scripts/reconstruction_sparse_cg_test: now uses Regularisator classes and plots
                                       phase differences (original - reconstr.)

pyramid/costfunction.py

@@ -8,11 +8,13 @@ import numpy as np
 from scipy.sparse.linalg import LinearOperator
 from scipy.sparse import eye
 
+from pyramid.regularisator import ZeroOrderRegularisator
+
 import logging
 
 
 class Costfunction(object):
-
+    # TODO: Docstrings!
     '''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
@@ -36,28 +38,32 @@ class Costfunction(object):
 
     LOG = logging.getLogger(__name__+'.Costfunction')
 
-    def __init__(self, y, fwd_model, lam=0):
+    def __init__(self, y, fwd_model, Se_inv=None, regularisator=None):
         self.LOG.debug('Calling __init__')
         self.y = y
         self.fwd_model = fwd_model
-        self.lam = lam
-        self.Se_inv = eye(len(y))
+        if Se_inv is None:
+            Se_inv = eye(len(y))
+        self.Se_inv = Se_inv
+        if regularisator is None:
+            regularisator = ZeroOrderRegularisator(lam=1E-4, size=3*fwd_model.size_3d)
+        self.regularisator = regularisator
         self.LOG.debug('Created '+str(self))
 
     def __call__(self, x):
         self.LOG.debug('Calling __call__')
-        y = self.y
-        F = self.fwd_model
-        Se_inv = self.Se_inv
-        return (F(x)-y).dot(Se_inv.dot(F(x)-y)) + lam * x.dot(Sa_inv.dot(x))  # TODO: implement
+        return ((self.fwd_model(x)-self.y).dot(self.Se_inv.dot(self.fwd_model(x)-self.y))
+                + self.regularisator(x))
 
     def __repr__(self):
         self.LOG.debug('Calling __repr__')
-        return '%s(fwd_model=%r, lam=%r)' % (self.__class__, self.fwd_model, self.lam)
+        return '%s(fwd_model=%r, Se_inv=%r, regularisator=%r)' % \
+            (self.__class__, self.fwd_model, self.Se_inv, self.regularisator)
 
     def __str__(self):
         self.LOG.debug('Calling __str__')
-        return 'Costfunction(fwd_model=%s, lam=%s)' % (self.fwd_model, self.lam)
+        return 'Costfunction(fwd_model=%s, Se_inv=%s, regularisator=%s)' % \
+            (self.fwd_model, self.Se_inv, self.regularisator)
 
     def jac(self, x):
         '''Calculate the derivative of the costfunction for a given magnetization distribution.
@@ -99,10 +105,8 @@ class Costfunction(object):
 
         '''
         self.LOG.debug('Calling hess_dot')
-        F = self.fwd_model
-        Se_inv = self.Se_inv
-        lam = self.lam
-        return F.jac_T_dot(x, Se_inv.dot(F.jac_dot(x, vector))) + lam*vector
+        return (self.fwd_model.jac_T_dot(x, self.Se_inv.dot(self.fwd_model.jac_dot(x, vector)))
+                + self.regularisator.jac_dot(vector))
 
 
 class CFAdapterScipyCG(LinearOperator):

pyramid/dataset.py

@@ -103,7 +103,8 @@ class DataSet(object):
         assert projector.dim_uv == self.dim_uv, 'Projector dimensions must match!'
         self.data.append((phase_map, projector))
 
-    def display_phase(self, phase_maps=None, title='Phase Map', cmap='RdBu', limit=None, norm=None):
+    def display_phase(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
@@ -134,7 +135,7 @@ class DataSet(object):
             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)]
+            for (i, phase_map) in enumerate(phase_maps)]
         plt.show()
 
     def display_combined(self, phase_maps=None, title='Combined Plot', cmap='RdBu', limit=None,
@@ -179,7 +180,7 @@ class DataSet(object):
             phase_maps = self.phase_maps
         [phase_map.display_combined('{} ({})'.format(title, self.projectors[i].get_info()),
                                     cmap, limit, norm, density, interpolation, grad_encode)
-            for (i, phase_map) in enumerate(self.phase_maps)]
+            for (i, phase_map) in enumerate(phase_maps)]
         plt.show()
 
     def create_phase_maps(self, mag_data):

pyramid/forwardmodel.py

@@ -49,9 +49,9 @@ class ForwardModel(object):
         self.a = kernel.a
         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.dim_uv = kernel.dim_uv
+        self.size_2d = kernel.size
         self.LOG.debug('Creating '+str(self))
 
     def __call__(self, x):

pyramid/phasemap.py

@@ -398,16 +398,16 @@ class PhaseMap(object):
         # Return plotting axis:
         return axis
 
-    def display_holo(self, density=1, title='Holographic Contour Map',
+    def display_holo(self, title='Holographic Contour Map', density=1,
                      axis=None, grad_encode='bright', interpolation='none', show=True):
         '''Display the color coded holography image.
 
         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 'Holographic Contour Map'.
+        density : float, optional
+            The gain factor for determining the number of contour lines. The default is 1.
         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

pyramid/reconstruction.py

@@ -83,7 +83,8 @@ class PrintIterator(object):
         return self.iteration
 
 
-def optimize_sparse_cg(data, verbosity=0):
+def optimize_sparse_cg(data, Se_inv=None, regularisator=None, verbosity=0):
+    # TODO: Docstring!
     '''Reconstruct a three-dimensional magnetic distribution from given phase maps via the
     conjugate gradient optimizaion method :func:`~.scipy.sparse.linalg.cg`.
 
@@ -109,7 +110,7 @@ def optimize_sparse_cg(data, verbosity=0):
     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)
+    cost = Costfunction(y, fwd_model, Se_inv, regularisator)
     # Optimize:
     A = CFAdapterScipyCG(cost)
     b = fwd_model.jac_T_dot(None, y)
@@ -225,8 +226,6 @@ def optimize_simple_leastsq(phase_map, mask, b_0=1, lam=1E-4, order=0):
         return np.concatenate([term1, term2])
 
     J_DICT = [J_0, J_1]  # list of cost-functions with different regularisations
-    print J_DICT
-    print J_DICT[order]
     # Reconstruct the magnetization components:
     x_rec, _ = leastsq(J_DICT[order], np.zeros(3*count))
     mag_data_rec.set_vector(mask, x_rec)

pyramid/regularisator.py

+# -*- coding: utf-8 -*-
+"""
+Created on Mon Aug 18 09:24:58 2014
+
+@author: Jan
+"""  # TODO: Docstring!
+
+
+import abc
+
+import numpy as np
+
+from scipy.sparse import eye
+
+import logging
+
+
+# TODO: Fragen für Jörn: Macht es Sinn, x_a standardmäßig auf den Nullvektor zu setzen? Ansonsten
+#       besser im jeweiligen Konstruktor setzen, nicht im abstrakten!
+#       Wie kommt man genau an die Ableitungen (Normen sind nicht unproblematisch)?
+#       Wie implementiert man am besten verschiedene Normen?
+
+
+class Regularisator(object):
+    # TODO: Docstring!
+
+    __metaclass__ = abc.ABCMeta
+    LOG = logging.getLogger(__name__+'.Regularisator')
+
+    @abc.abstractmethod
+    def __init__(self, Sa_inv, x_a=None):
+        self.LOG.debug('Calling __init__')
+        self.Sa_inv = Sa_inv
+        if x_a is None:
+            x_a = np.zeros(np.shape(Sa_inv)[1])
+        self.x_a = x_a
+        self.LOG.debug('Created '+str(self))
+
+    def __call__(self, x, x_a=None):
+        self.LOG.debug('Calling __call__')
+        return (x-self.x_a).dot(self.Sa_inv.dot(x-self.x_a))
+
+    def __repr__(self):
+        self.LOG.debug('Calling __repr__')
+        return '%s(Sa_inv=%r, x_a=%r)' % (self.__class__, self.Sa_inv, self.x_a)
+
+    def __str__(self):
+        self.LOG.debug('Calling __str__')
+        return 'Regularisator(Sa_inv=%s, x_a=%s)' % (self.Sa_inv, self.x_a)
+
+    def jac_dot(self, vector):
+        # TODO: Docstring!
+        return self.Sa_inv.dot(vector-self.x_a)
+
+
+class ZeroOrderRegularisator(Regularisator):
+    # TODO: Docstring!
+
+    def __init__(self, lam, size, x_a=None):
+        Sa_inv = lam * eye(size)
+        super(ZeroOrderRegularisator, self).__init__(Sa_inv, x_a)
+        self.LOG.debug('Created '+str(self))
+
+
+class FirstOrderRegularisator(Regularisator):
+    # TODO: Docstring!
+
+    def __init__(self, lam, size, x_a=None):
+        Sa_inv = lam * eye(size)
+        super(FirstOrderRegularisator, self).__init__(Sa_inv, x_a)
+        self.LOG.debug('Created '+str(self))

scripts/collaborations/rueffner_file.py

@@ -65,8 +65,9 @@ xs = np.arange(-dim[2]/2, dim[2]/2)
 ys = np.arange(-dim[1]/2, dim[1]/2)
 zs = np.arange(-dim[0]/2, dim[0]/2)
 xx, yy = np.meshgrid(xs, ys)
-# Interpolation and phase calculation for all timesteps:
-for t in np.arange(865, 1001, 5):
+
+
+def calculate(t):  # TODO: Somehow avoid memory error :-(...
     print 't =', t
     vectors = h5file.root.data.fields.m.read(field='m_CoFeb')[t, ...]
     data = np.hstack((points, vectors))
@@ -92,5 +93,12 @@ for t in np.arange(865, 1001, 5):
     plt.savefig(PATH+'rueffner/phase_map_z_t_'+str(t)+'.png')
     phase_map_x.display_combined(density=dens_x, interpolation='bilinear', limit=lim_x)
     plt.savefig(PATH+'rueffner/phase_map_x_t_'+str(t)+'.png')
+    vectors, data, magnitude, z_slice, weights, grid_x, grid_y, grid_z, grid, mag_data, \
+        phase_map_z, phase_map_x = [None] * 12
     # Close all plots to avoid clutter:
     plt.close('all')
+
+
+# Interpolation and phase calculation for all timesteps:
+for t in np.arange(0, 1001, 5):
+    calculate(t)

scripts/collaborations/zi_an.py

@@ -95,3 +95,16 @@ plt.savefig(PATH+'phase_map_4part_diff.png')
 # Get the average difference from the experimental results:
 print 'Average difference (2 cubes):', np.average(phase_diff_2.phase)
 print 'Average difference (4 cubes):', np.average(phase_diff_4.phase)
+# Plot holographic contour maps with overlayed magnetic distributions:
+axis = phase_map_rec_2.display_holo('Magnetization Overlay', density=density, interpolation=inter)
+(mag_data_rec_2*(1/mag_data_rec_2.magnitude.max())).quiver_plot(axis=axis)
+axis = plt.gca()
+axis.set_xlim(20, 45)
+axis.set_ylim(20, 45)
+plt.savefig(PATH+'phase_map_2part_holo.png')
+axis = phase_map_rec_4.display_holo('Magnetization Overlay', density=density, interpolation=inter)
+(mag_data_rec_4*(1/mag_data_rec_4.magnitude.max())).quiver_plot(axis=axis)
+axis = plt.gca()
+axis.set_xlim(20, 45)
+axis.set_ylim(20, 45)
+plt.savefig(PATH+'phase_map_4part_holo.png')

scripts/reconstruction/reconstruction_sparse_cg_test.py

@@ -15,11 +15,8 @@ from pyramid.phasemap import PhaseMap
 from pyramid.projector import YTiltProjector, XTiltProjector
 from pyramid.phasemapper import PMConvolve
 from pyramid.dataset import DataSet
+from pyramid.regularisator import ZeroOrderRegularisator
 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
 
 import logging
@@ -85,24 +82,25 @@ else:
 print('--Setting up data collection')
 
 dim_uv = dim[1:3]
-
-lam = 10**-10
-
 data = DataSet(a, dim_uv, b_0)
 [data.append((phase_maps[i], projectors[i])) for i in range(len(projectors))]
-y = data.phase_vec
-kern = Kernel(data.a, data.dim_uv, data.b_0)
-F = ForwardModel(data.projectors, kern)
-C = Costfunction(y, F, lam)
 
-data.display()
+###################################################################################################
+print('--Test simple solver')
+
+lam = 1E-10
+reg = ZeroOrderRegularisator(lam, 3*np.prod(dim))
+start = clock()
+mag_data_opt = rc.optimize_sparse_cg(data, regularisator=reg, verbosity=1)
+print 'Time:', str(clock()-start)
+mag_data_opt.quiver_plot3d()
+(mag_data_opt - mag_data).quiver_plot3d()
 
 ###################################################################################################
-#print('--Test simple solver')
-#
-#start = clock()
-#mag_data_opt = rc.optimize_sparse_cg(data, verbosity=1)
-#print 'Time:', str(clock()-start)
-#mag_data_opt.quiver_plot3d()
-#(mag_data_opt - mag_data).quiver_plot3d()
-##data.display(data.create_phase_maps(mag_data_opt))
+print('--Plot stuff')
+
+phase_maps_opt = data.create_phase_maps(mag_data_opt)
+#data.display_phase()
+#data.display_phase(phase_maps_opt)
+phase_diffs = [(phase_maps[i]-phase_maps_opt[i]) for i in range(len(data.data))]
+data.display_phase(phase_diffs)

scripts/test methods/compare_pixel_fields.py

@@ -75,32 +75,32 @@ print 'Fourier Kernel, direct and indirect method are identical:', \
       np.all(phase_map_fft_kernel.phase - phase_map_fft.phase) == 0
 (phase_map_disc-phase_map_fft).display_phase('Phase difference of one Pixel (Disc - Fourier)')
 
-## Cross section plots of real space kernels:
-#fig = plt.figure()
-#axis = fig.add_subplot(1, 1, 1)
-#x = np.linspace(-dim[1]/a/2, dim[1]/a/2-1, dim[1])
-#y_ft = phase_map_fft.phase[:, dim[1]/2]
-#y_re = phase_map_disc.phase[:, dim[1]/2]
-#axis.axhline(0, color='k')
-#axis.plot(x, y_re, label='Real space method')
-#axis.plot(x, y_ft, label='Fourier method')
-#axis.grid()
-#axis.legend()
-#axis.set_title('Real Space Kernel')
-#axis.set_xlim(-dim[1]/2, dim[1]/2-1)
-#axis.xaxis.set_major_locator(IndexLocator(base=dim[1]/8, offset=0))
-#
-## Cross section plots of Fourier space kernels:
-#fig = plt.figure()
-#axis = fig.add_subplot(1, 1, 1)
-#x = range(dim[1])
-#k_re = np.abs(np.fft.ifftshift(np.fft.rfft2(phase_map_disc.phase), axes=0))[:, 0]**2
-#k_ft = np.abs(np.fft.ifftshift(np.fft.rfft2(phase_map_fft.phase), axes=0))[:, 0]**2
-#axis.axhline(0, color='k')
-#axis.plot(x, k_re, label='Real space method')
-#axis.plot(x, k_ft, label='Fourier method')
-#axis.grid()
-#axis.legend()
-#axis.set_title('Fourier Space Kernel')
-#axis.set_xlim(0, dim[1]-1)
-#axis.xaxis.set_major_locator(IndexLocator(base=dim[1]/8, offset=0))
+# Cross section plots of real space kernels:
+fig = plt.figure()
+axis = fig.add_subplot(1, 1, 1)
+x = np.linspace(-dim[1]/a/2, dim[1]/a/2-1, dim[1])
+y_ft = phase_map_fft.phase[:, dim[1]/2]
+y_re = phase_map_disc.phase[:, dim[1]/2]
+axis.axhline(0, color='k')
+axis.plot(x, y_re, label='Real space method')
+axis.plot(x, y_ft, label='Fourier method')
+axis.grid()
+axis.legend()
+axis.set_title('Real Space Kernel')
+axis.set_xlim(-dim[1]/2, dim[1]/2-1)
+axis.xaxis.set_major_locator(IndexLocator(base=dim[1]/8, offset=0))
+
+# Cross section plots of Fourier space kernels:
+fig = plt.figure()
+axis = fig.add_subplot(1, 1, 1)
+x = range(dim[1])
+k_re = np.abs(np.fft.ifftshift(np.fft.rfft2(phase_map_disc.phase), axes=0))[:, 0]**2
+k_ft = np.abs(np.fft.ifftshift(np.fft.rfft2(phase_map_fft.phase), axes=0))[:, 0]**2
+axis.axhline(0, color='k')
+axis.plot(x, k_re, label='Real space method')
+axis.plot(x, k_ft, label='Fourier method')
+axis.grid()
+axis.legend()
+axis.set_title('Fourier Space Kernel')
+axis.set_xlim(0, dim[1]-1)
+axis.xaxis.set_major_locator(IndexLocator(base=dim[1]/8, offset=0))