From 325ea464526307d54d9fff2fdf221fe6e10c3c8c Mon Sep 17 00:00:00 2001
From: Jan Caron <j.caron@fz-juelich.de>
Date: Tue, 19 Aug 2014 08:42:43 +0200
Subject: [PATCH] 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 | 32 +++++----
pyramid/dataset.py | 7 +-
pyramid/forwardmodel.py | 4 +-
pyramid/phasemap.py | 6 +-
pyramid/reconstruction.py | 7 +-
pyramid/regularisator.py | 71 +++++++++++++++++++
scripts/collaborations/rueffner_file.py | 12 +++-
scripts/collaborations/zi_an.py | 13 ++++
.../reconstruction_sparse_cg_test.py | 38 +++++-----
scripts/test methods/compare_pixel_fields.py | 58 +++++++--------
10 files changed, 171 insertions(+), 77 deletions(-)
create mode 100644 pyramid/regularisator.py
diff --git a/pyramid/costfunction.py b/pyramid/costfunction.py
index 275d908..940df35 100644
--- a/pyramid/costfunction.py
+++ b/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):
diff --git a/pyramid/dataset.py b/pyramid/dataset.py
index f51ddfc..5496239 100644
--- a/pyramid/dataset.py
+++ b/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):
diff --git a/pyramid/forwardmodel.py b/pyramid/forwardmodel.py
index d820f72..eeb4868 100644
--- a/pyramid/forwardmodel.py
+++ b/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):
diff --git a/pyramid/phasemap.py b/pyramid/phasemap.py
index a2f6f1a..061ade2 100644
--- a/pyramid/phasemap.py
+++ b/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
diff --git a/pyramid/reconstruction.py b/pyramid/reconstruction.py
index e251d28..22fd35b 100644
--- a/pyramid/reconstruction.py
+++ b/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)
diff --git a/pyramid/regularisator.py b/pyramid/regularisator.py
new file mode 100644
index 0000000..b3621b0
--- /dev/null
+++ b/pyramid/regularisator.py
@@ -0,0 +1,71 @@
+# -*- 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))
diff --git a/scripts/collaborations/rueffner_file.py b/scripts/collaborations/rueffner_file.py
index 6b37fea..ef0c2ac 100644
--- a/scripts/collaborations/rueffner_file.py
+++ b/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)
diff --git a/scripts/collaborations/zi_an.py b/scripts/collaborations/zi_an.py
index e03e2d3..85df772 100644
--- a/scripts/collaborations/zi_an.py
+++ b/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')
diff --git a/scripts/reconstruction/reconstruction_sparse_cg_test.py b/scripts/reconstruction/reconstruction_sparse_cg_test.py
index c5a6fa1..a4bb6a9 100644
--- a/scripts/reconstruction/reconstruction_sparse_cg_test.py
+++ b/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)
diff --git a/scripts/test methods/compare_pixel_fields.py b/scripts/test methods/compare_pixel_fields.py
index bda9f0c..c26bb99 100644
--- a/scripts/test methods/compare_pixel_fields.py
+++ b/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))
--
GitLab