From 160612def890e3cd43dfc0dcbc0331be7e7375cc Mon Sep 17 00:00:00 2001
From: "Joern Ungermann (IEK-7 FZ-Juelich)" <j.ungermann@fz-juelich.de>
Date: Thu, 30 Oct 2014 20:31:42 +0100
Subject: [PATCH] update of documentation and tidying of split bregman.
---
pyramid/reconstruction.py | 74 ++++++++++++++++++++++-----------------
1 file changed, 42 insertions(+), 32 deletions(-)
diff --git a/pyramid/reconstruction.py b/pyramid/reconstruction.py
index 4f85114..239fd49 100644
--- a/pyramid/reconstruction.py
+++ b/pyramid/reconstruction.py
@@ -50,7 +50,7 @@ class PrintIterator(object):
'''
- LOG = logging.getLogger(__name__+'.PrintIterator')
+ LOG = logging.getLogger(__name__ + '.PrintIterator')
def __init__(self, cost, verbosity):
self.LOG.debug('Calling __init__')
@@ -58,7 +58,7 @@ class PrintIterator(object):
self.verbosity = verbosity
assert verbosity in {0, 1, 2}, 'verbosity has to be set to 0, 1 or 2!'
self.iteration = 0
- self.LOG.debug('Created '+str(self))
+ self.LOG.debug('Created ' + str(self))
def __call__(self, xk):
self.LOG.debug('Calling __call__')
@@ -86,6 +86,7 @@ class PrintIterator(object):
def optimize_linear(data, regularisator=None, 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
----------
@@ -104,21 +105,22 @@ def optimize_linear(data, regularisator=None, max_iter=None):
'''
import jutil.cg as jcg
- LOG.debug('Calling optimize_sparse_cg')
+ 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)
LOG.info('Cost after optimization: {}'.format(cost(x_opt)))
# Create and return fitting MagData object:
- mag_opt = MagData(data.a, np.zeros((3,)+data.dim))
+ mag_opt = MagData(data.a, np.zeros((3,) + data.dim))
mag_opt.set_vector(x_opt, data.mask)
return mag_opt
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`.
+ '''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
----------
@@ -129,8 +131,7 @@ def optimize_nonlin(data, first_guess=None, regularisator=None):
first_fuess : :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. Defaults to zero
- order Tikhonov if none is provided.
+ Regularisator class that's responsible for the regularisation term.
Returns
-------
@@ -140,9 +141,9 @@ def optimize_nonlin(data, first_guess=None, regularisator=None):
'''
import jutil.minimizer as jmin
import jutil.norms as jnorms
- LOG.debug('Calling optimize_cg')
+ LOG.debug('Calling optimize_nonlin')
if first_guess is None:
- first_guess = MagData(data.a, np.zeros((3,)+data.dim))
+ first_guess = MagData(data.a, np.zeros((3,) + data.dim))
x_0 = first_guess.get_vector(data.mask)
cost = Costfunction(data, regularisator)
@@ -152,12 +153,13 @@ def optimize_nonlin(data, first_guess=None, regularisator=None):
q = 1. / (1. - (1. / p))
lp = regularisator.norm
lq = jnorms.LPPow(q, 1e-20)
+
def preconditioner(_, direc):
direc_p = direc / abs(direc).max()
direc_p = 10 * (1. / q) * lq.jac(direc_p)
return direc_p
- # This Method is semi-best for L1 type problems. Takes forever, though
+ # This Method is semi-best for Lp type problems. Takes forever, though
LOG.info('Cost before optimization: {}'.format(cost(np.zeros(cost.n))))
result = jmin.minimize(
cost, x_0,
@@ -166,14 +168,18 @@ def optimize_nonlin(data, first_guess=None, regularisator=None):
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 = MagData(data.a, np.zeros((3,) + data.dim))
mag_opt.set_vector(x_opt, data.mask)
return mag_opt
def optimize_splitbregman(data, weight, lam, mu):
- '''Reconstruct a three-dimensional magnetic distribution from given phase maps via the
- conjugate gradient optimizaion method :func:`~.scipy.sparse.linalg.cg`.
+ '''
+ Reconstructs magnet distribution from phase image measurements using a split bregman
+ algorithm with a dedicated TV-l1 norm. Very dedicated, frickle, brittle, and difficult
+ to get to work, but fastest option available if it works.
+
+ Seems to work for some 2D examples with weight=lam=1 and mu in [1, .., 1e4].
Parameters
----------
@@ -181,11 +187,12 @@ def optimize_splitbregman(data, weight, lam, mu):
: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.
- first_fuess : :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. Defaults to zero
- order Tikhonov if none is provided.
+ weight : float
+ Obscure split bregman parameter
+ lam : float
+ Cryptic split bregman parameter
+ mu : float
+ flabberghasting split bregman paramter
Returns
-------
@@ -198,30 +205,33 @@ def optimize_splitbregman(data, weight, lam, mu):
import jutil.diff as jdiff
from pyramid.regularisator import FirstOrderRegularisator
LOG.debug('Calling optimize_splitbregman')
- regularisator = FirstOrderRegularisator(data.mask, lam, 1)
- x_0 = MagData(data.a, np.zeros((3,)+data.dim)).get_vector(data.mask)
- cost = Costfunction(data, regularisator)
- print cost(x_0)
+ # 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)
+ x_0 = MagData(data.a, np.zeros((3,) + data.dim)).get_vector(data.mask)
+ cost = Costfunction(data, regularisator)
fwd_mod = cost.fwd_model
+
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))
- D0 = jdiff.get_diff_operator(data.mask, 0, 3)
- D1 = jdiff.get_diff_operator(data.mask, 1, 3)
- D = joperator.VStack([D0, D1])
+ 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)
+
x_opt = jsb.split_bregman_2d(
A, D, y,
- weight=1000, mu=10, lambd=1e-6, max_iter=100)
- print cost(x_opt)
- mag_opt = MagData(data.a, np.zeros((3,)+data.dim))
+ 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)
return mag_opt
-
def optimize_simple_leastsq(phase_map, mask, b_0=1, lam=1E-4, order=0):
'''Reconstruct a magnetic distribution for a 2-D problem with known pixel locations.
@@ -259,7 +269,7 @@ def optimize_simple_leastsq(phase_map, mask, b_0=1, lam=1E-4, order=0):
dim = mask.shape # Dimensions of the mag. distr.
count = mask.sum() # Number of pixels with magnetization
# Create empty MagData object for the reconstruction:
- mag_data_rec = MagData(a, np.zeros((3,)+dim))
+ mag_data_rec = MagData(a, np.zeros((3,) + dim))
# Function that returns the phase map for a magnetic configuration x:
def F(x):
@@ -294,6 +304,6 @@ def optimize_simple_leastsq(phase_map, mask, b_0=1, lam=1E-4, order=0):
# 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))
+ x_rec, _ = leastsq(J_DICT[order], np.zeros(3 * count))
mag_data_rec.set_vector(x_rec, mask)
return mag_data_rec
--
GitLab