From c19f84c159f951dd764d25b7411fc984df76dc8c Mon Sep 17 00:00:00 2001
From: "Joern Ungermann (IEK-7 FZ-Juelich)" <j.ungermann@fz-juelich.de>
Date: Thu, 30 Oct 2014 12:23:22 +0100
Subject: [PATCH] Added a steppest descent minimizer for lp problems and added
an experimental split-bregman optimizer.
---
pyramid/reconstruction.py | 119 +++++++++++++++++++++++++-------------
1 file changed, 78 insertions(+), 41 deletions(-)
diff --git a/pyramid/reconstruction.py b/pyramid/reconstruction.py
index 237ef8e..b5a1500 100644
--- a/pyramid/reconstruction.py
+++ b/pyramid/reconstruction.py
@@ -83,7 +83,7 @@ class PrintIterator(object):
return self.iteration
-def optimize_linear(data, regularisator=None, maxiter=1000, verbosity=0):
+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`.
@@ -93,19 +93,9 @@ def optimize_linear(data, regularisator=None, maxiter=1000, verbosity=0):
: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.
- Se_inv : :class:`~numpy.ndarray` (N=2), optional
- Inverted covariance matrix of the measurement errors. The matrix has size `NxN` with N
- being the length of the targetvector y (vectorized phase map information). Defaults to
- an appropriate unity matrix if none is provided.
regularisator : :class:`~.Regularisator`, optional
Regularisator class that's responsible for the regularisation term. Defaults to zero
order Tikhonov if none is provided.
- maxiter : int
- Maximum number of iterations.
- verbosity : {0, 1, 2}, optional
- Parameter defining the verposity of the output. `2` will show the current number of the
- iteration and the cost of the current distribution. `1` will just show the iteration
- number and `0` is the default and will prevent output all together.
Returns
-------
@@ -113,13 +103,13 @@ def optimize_linear(data, regularisator=None, maxiter=1000, verbosity=0):
The reconstructed magnetic distribution as a :class:`~.MagData` object.
'''
- import jutil
+ import jutil.cg as jcg
LOG.debug('Calling optimize_sparse_cg')
# Set up necessary objects:
cost = Costfunction(data, regularisator)
- print cost(np.zeros(cost.n))
- x_opt = jutil.cg.conj_grad_minimize(cost, max_iter=20)
- print cost(x_opt)
+ 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.set_vector(x_opt, data.mask)
@@ -136,10 +126,11 @@ def optimize_nonlin(data, first_guess=None, regularisator=None):
: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.
- verbosity : {2, 1, 0}, optional
- Parameter defining the verposity of the output. `2` is the default and will show the
- current number of the iteration and the cost of the current distribution. `2` will just
- show the iteration number and `0` will prevent output all together.
+ 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.
Returns
-------
@@ -147,44 +138,90 @@ def optimize_nonlin(data, first_guess=None, regularisator=None):
The reconstructed magnetic distribution as a :class:`~.MagData` object.
'''
- import jutil
+ import jutil.minimizer as jmin
+ import jutil.norms as jnorms
LOG.debug('Calling optimize_cg')
if first_guess is None:
first_guess = MagData(data.a, np.zeros((3,)+data.dim))
+
x_0 = first_guess.get_vector(data.mask)
cost = Costfunction(data, regularisator)
+ assert len(x_0) == cost.n, (len(x_0), cost.m, cost.n)
-# proj = fwd_model.data_set.projectors[0]
-# proj_jac1 = np.array([proj.jac_dot(np.eye(proj.m, 1, -i).squeeze()) for i in range(proj.m)])
-# proj_jac2 = np.array([proj.jac_T_dot(np.eye(proj.n, 1, -i).squeeze()) for i in range(proj.n)])
-# print jac1, jac2.T, abs(jac1-jac2.T).sum()
-# print jac1.shape, jac2.shape
+ p = regularisator.p
+ q = 1. / (1. - (1. / p))
+ lp = regularisator.norm
+ lq = jnorms.LPPow(q, 1e-20)
+ def preconditioner(_, direc):
+ direc_p = direc / 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
+ LOG.info('Cost before optimization: {}'.format(cost(np.zeros(cost.n))))
+ result = jmin.minimize(
+ cost, x_0,
+ method="SteepestDescent",
+ options={"preconditioner": preconditioner},
+ 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.set_vector(x_opt, data.mask)
+ return mag_opt
-# pm = fwd_model.phase_mappers[proj.dim_uv]
-# pm_jac1 = np.array([pm.jac_dot(np.eye(pm.m)[:, i]) for i in range(pm.m)])
-# pm_jac2 = np.array([pm.jac_T_dot(np.eye(pm.n)[:, i]) for i in range(pm.n)])
-# print jac1, jac2.T, abs(jac1-jac2.T).sum()
-# print jac1.shape, jac2.shape
+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`.
- # jac1 = np.array([fwd_model.jac_dot(x_0, np.eye(fwd_model.m)[:, i]) for i in range(fwd_model.m)])
- # jac2 = np.array([fwd_model.jac_T_dot(x_0, np.eye(fwd_model.n)[:, i]) for i in range(fwd_model.n)])
- # print proj_jac1.dot(pm_jac1)
- # print (pm_jac2.dot(proj_jac2)).T
- # print jac1
-# print jac2.T
-# print abs(jac1-jac2.T).sum()
-# print jac1.shape, jac2.shape
+ Parameters
+ ----------
+ data : :class:`~.DataSet`
+ :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.
- assert len(x_0) == cost.n, (len(x_0), cost.m, cost.n)
- result = jutil.minimizer.minimize(cost, x_0, options={"conv_rel":1e-2}, tol={"max_iteration":4})
- x_opt = result.x
+ Returns
+ -------
+ mag_data : :class:`~pyramid.magdata.MagData`
+ The reconstructed magnetic distribution as a :class:`~.MagData` object.
+
+ '''
+ import jutil.splitbregman as jsb
+ import jutil.operator as joperator
+ 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)
+
+ 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])
+ 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))
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.
--
GitLab