From 76c176702cdfe3e1088f01ebc62d269c1e98d55c Mon Sep 17 00:00:00 2001
From: "Joern Ungermann (IEK-7 FZ-Juelich)" <j.ungermann@fz-juelich.de>
Date: Wed, 29 Oct 2014 16:09:23 +0100
Subject: [PATCH] fixed some minor compilation bugs and updated to current
jutil. prepared jutil reconstruction routines.
---
pyramid/costfunction.py | 16 +++++++++---
pyramid/reconstruction.py | 34 ++++++++++++-------------
pyramid/regularisator.py | 4 +--
scripts/collaborations/example_joern.py | 33 ++++++++++++++++--------
4 files changed, 55 insertions(+), 32 deletions(-)
diff --git a/pyramid/costfunction.py b/pyramid/costfunction.py
index 04ed291..ad7c226 100644
--- a/pyramid/costfunction.py
+++ b/pyramid/costfunction.py
@@ -53,6 +53,8 @@ class Costfunction(object):
# Extract important information:
self.y = data_set.phase_vec
self.Se_inv = data_set.Se_inv
+ self.n = data_set.n * 3
+ self.m = len(self.y)
self.LOG.debug('Created '+str(self))
def __repr__(self):
@@ -65,10 +67,14 @@ class Costfunction(object):
return 'Costfunction(data_set=%s, fwd_model=%s, regularisator=%s)' % \
(self.data_set, self.fwd_model, self.regularisator)
+ def init(self, x):
+ self(x)
+
def __call__(self, x):
self.LOG.debug('Calling __call__')
- delta_y = self.fwd_model(x)-self.y
- return delta_y.dot(self.Se_inv.dot(delta_y)) + self.regularisator(x)
+ delta_y = self.fwd_model(x) - self.y
+ self.chisq = delta_y.dot(self.Se_inv.dot(delta_y)) + self.regularisator(x)
+ return self.chisq
def jac(self, x):
'''Calculate the derivative of the costfunction for a given magnetization distribution.
@@ -85,7 +91,8 @@ class Costfunction(object):
'''
self.LOG.debug('Calling jac')
- return (2 * self.fwd_model.jac_T_dot(x, self.Se_inv.dot(self.fwd_model(x)-self.y))
+ assert len(x) == self.n
+ return (2 * self.fwd_model.jac_T_dot(x, self.Se_inv.dot(self.fwd_model(x) - self.y))
+ self.regularisator.jac(x))
def hess_dot(self, x, vector):
@@ -111,6 +118,9 @@ class Costfunction(object):
return (2 * self.fwd_model.jac_T_dot(x, self.Se_inv.dot(self.fwd_model.jac_dot(x, vector)))
+ self.regularisator.hess_dot(x, vector))
+ def hess_diag(self, _):
+ return np.ones(self.n)
+
class CFAdapterScipyCG(LinearOperator):
diff --git a/pyramid/reconstruction.py b/pyramid/reconstruction.py
index 90704d5..1eff2ea 100644
--- a/pyramid/reconstruction.py
+++ b/pyramid/reconstruction.py
@@ -83,7 +83,7 @@ class PrintIterator(object):
return self.iteration
-def optimize_sparse_cg(data, regularisator=None, maxiter=1000, verbosity=0):
+def optimize_linear(data, regularisator=None, maxiter=1000, verbosity=0):
'''Reconstruct a three-dimensional magnetic distribution from given phase maps via the
conjugate gradient optimizaion method :func:`~.scipy.sparse.linalg.cg`.
@@ -113,22 +113,21 @@ def optimize_sparse_cg(data, regularisator=None, maxiter=1000, verbosity=0):
The reconstructed magnetic distribution as a :class:`~.MagData` object.
'''
+ import jutil
LOG.debug('Calling optimize_sparse_cg')
# Set up necessary objects:
- y = data.phase_vec
fwd_model = ForwardModel(data)
cost = Costfunction(data, regularisator)
- # Optimize:
- A = CFAdapterScipyCG(cost)
- b = fwd_model.jac_T_dot(None, y)
- x_opt, info = cg(A, b, callback=PrintIterator(cost, verbosity), maxiter=maxiter)
+ print cost(np.zeros(cost.n))
+ x_opt = jutil.cg.conj_grad_minimize(cost)
+ print cost(x_opt)
# Create and return fitting MagData object:
mag_opt = MagData(data.a, np.zeros((3,)+data.dim))
mag_opt.mag_vec = x_opt
return mag_opt
-def optimize_cg(data, first_guess):
+def optimize_nonlin(data, first_guess=None, regularisator=None):
'''Reconstruct a three-dimensional magnetic distribution from given phase maps via the
conjugate gradient optimizaion method :func:`~.scipy.sparse.linalg.cg`.
@@ -149,19 +148,18 @@ def optimize_cg(data, first_guess):
The reconstructed magnetic distribution as a :class:`~.MagData` object.
'''
+ import jutil
LOG.debug('Calling optimize_cg')
- mag_0 = first_guess
+ if first_guess is None:
+ first_guess = MagData(data.a, np.zeros((3,)+data.dim))
x_0 = first_guess.mag_vec
- y = data.phase_vec
- kernel = Kernel(data.a, data.dim_uv)
- fwd_model = ForwardModel(data.projectors, kernel, data.b_0)
- cost = Costfunction(y, fwd_model)
- # Optimize:
- result = minimize(cost, x_0, method='Newton-CG', jac=cost.jac, hessp=cost.hess_dot,
- options={'maxiter': 1000, 'disp': True})
- # Create optimized MagData object:
+ fwd_model = ForwardModel(data)
+ cost = Costfunction(data, regularisator)
+ assert len(x_0) == cost.n, (len(x_0), cost.m, cost.n)
+ result = jutil.minimizer.minimize(cost, x_0, options={"conv_rel":1e-20}, tol={"max_iteration":1})
x_opt = result.x
- mag_opt = MagData(mag_0.a, np.zeros((3,)+mag_0.dim))
+ print cost(x_opt)
+ mag_opt = MagData(data.a, np.zeros((3,)+data.dim))
mag_opt.mag_vec = x_opt
return mag_opt
@@ -236,6 +234,8 @@ def optimize_simple_leastsq(phase_map, mask, b_0=1, lam=1E-4, order=0):
J_DICT = [J_0, J_1] # list of cost-functions with different regularisations
# Reconstruct the magnetization components:
+ # TODO Use jutil.minimizer.minimize(jutil.costfunction.LeastSquaresCostFunction(J_DICT[order],
+ # ...) or a simpler frontend.
x_rec, _ = leastsq(J_DICT[order], np.zeros(3*count))
mag_data_rec.set_vector(x_rec, mask)
return mag_data_rec
diff --git a/pyramid/regularisator.py b/pyramid/regularisator.py
index a7adf41..bd06648 100644
--- a/pyramid/regularisator.py
+++ b/pyramid/regularisator.py
@@ -52,7 +52,7 @@ class Regularisator(object):
def jac(self, x):
# TODO: Docstring!
self.LOG.debug('Calling jac')
- return self.lam * self.norm.jac_dot(x)
+ return self.lam * self.norm.jac(x)
def hess_dot(self, x, vector):
# TODO: Docstring!
@@ -105,7 +105,7 @@ class ZeroOrderRegularisator(Regularisator):
def __init__(self, lam):
self.LOG.debug('Calling __init__')
- norm = jnorm.NormL2Square()
+ norm = jnorm.L2Square()
super(ZeroOrderRegularisator, self).__init__(norm, lam)
self.LOG.debug('Created '+str(self))
diff --git a/scripts/collaborations/example_joern.py b/scripts/collaborations/example_joern.py
index 7890c35..e3b1a4d 100644
--- a/scripts/collaborations/example_joern.py
+++ b/scripts/collaborations/example_joern.py
@@ -28,6 +28,7 @@ LOGGING_CONF = os.path.join(os.path.dirname(os.path.realpath(pyramid.__file__)),
logging.config.fileConfig(LOGGING_CONF, disable_existing_loggers=False)
+logging.basicConfig(level=logging.INFO)
###################################################################################################
threshold = 1
a = 1.0 # in nm
@@ -40,12 +41,13 @@ smoothed_pictures = True
lam = 1E-4
order = 1
log = True
-PATH = '../../output/joern/'
+PATH = './'
+dirname = PATH
###################################################################################################
# Read in files:
phase_map = PhaseMap.load_from_netcdf4(PATH+'phase_map.nc')
#mask = np.genfromtxt('mask.txt', dtype=bool)
-with open(PATH+'mask.pickle') as pf:
+with open(PATH + 'mask.pickle') as pf:
mask = pickle.load(pf)
# Setup:
data_set = DataSet(a, dim, b_0)
@@ -53,29 +55,40 @@ data_set.append(phase_map, SimpleProjector(dim))
regularisator = ZeroOrderRegularisator(lam)
# Reconstruct the magnetic distribution:
tic = clock()
-mag_data_rec = rc.optimize_sparse_cg(data_set, regularisator, maxiter=100, verbosity=2)
+mag_data_rec1 = rc.optimize_linear(data_set, regularisator=regularisator)
+mag_data_rec = rc.optimize_nonlin(data_set, regularisator=regularisator)
# .optimize_simple_leastsq(phase_map, mask, b_0, lam=lam, order=order)
+
print 'reconstruction time:', clock() - tic
# Display the reconstructed phase map and holography image:
phase_map_rec = pm(mag_data_rec)
-phase_map_rec.display_combined('Reconstr. Distribution', gain=gain, interpolation=inter)
+phase_map_rec.display_combined('Reconstr. Distribution', gain=gain, interpolation=inter, show=False)
+plt.savefig(dirname + "/reconstr.png")
+
# Plot the magnetization:
-axis = (mag_data_rec*(1/mag_data_rec.magnitude.max())).quiver_plot()
+axis = (mag_data_rec*(1/mag_data_rec.magnitude.max())).quiver_plot(show=False)
axis.set_xlim(20, 45)
axis.set_ylim(20, 45)
+plt.savefig(dirname + "/quiver.png")
+
# Display the Phase:
phase_diff = phase_map_rec-phase_map
-phase_diff.display_phase('Difference')
+phase_diff.display_phase('Difference', show=False)
+plt.savefig(dirname + "/difference.png")
+
# Get the average difference from the experimental results:
print 'Average difference:', np.average(phase_diff.phase)
# Plot holographic contour maps with overlayed magnetic distributions:
-axis = phase_map_rec.display_holo('Magnetization Overlay', gain=0.1, interpolation=inter)
-mag_data_rec.quiver_plot(axis=axis)
+axis = phase_map_rec.display_holo('Magnetization Overlay', gain=0.1, interpolation=inter, show=False)
+mag_data_rec.quiver_plot(axis=axis, show=False)
axis = plt.gca()
axis.set_xlim(20, 45)
axis.set_ylim(20, 45)
-axis = phase_map_rec.display_holo('Magnetization Overlay', gain=0.1, interpolation=inter)
-mag_data_rec.quiver_plot(axis=axis, log=log)
+plt.savefig(dirname + "/overlay_normal.png")
+
+axis = phase_map_rec.display_holo('Magnetization Overlay', gain=0.1, interpolation=inter, show=False)
+mag_data_rec.quiver_plot(axis=axis, log=log, show=False)
axis = plt.gca()
axis.set_xlim(20, 45)
axis.set_ylim(20, 45)
+plt.savefig(dirname + "/overlay_log.png")
--
GitLab