Merge

pyramid/costfunction.py

@@ -10,7 +10,7 @@ from scipy.sparse import eye
 
 from pyramid.forwardmodel import ForwardModel
 from pyramid.regularisator import ZeroOrderRegularisator
-
+from pyramid.regularisator import NoneRegularisator
 import logging
 
 
@@ -54,6 +54,8 @@ class Costfunction(object):
         self.data_set = data_set
         self.fwd_model = ForwardModel(data_set)
         self.regularisator = regularisator
+        if self.regularisator is None:
+            self.regularisator = NoneRegularisator()
         # Extract important information:
         self.y = data_set.phase_vec
         self.Se_inv = data_set.Se_inv

pyramid/magdata.py

@@ -11,9 +11,6 @@ from scipy.ndimage.interpolation import zoom
 import matplotlib.pyplot as plt
 import matplotlib.cm as cmx
 from matplotlib.ticker import MaxNLocator
-from mayavi import mlab
-
-from lxml import etree
 
 from numbers import Number
 
@@ -519,6 +516,8 @@ class MagData(object):
 
         '''
         self.LOG.debug('Calling quiver_plot3D')
+        from mayavi import mlab
+
         a = self.a
         dim = self.dim
         # Create points and vector components as lists:
@@ -553,6 +552,8 @@ class MagData(object):
 
         '''
         self.LOG.debug('Calling save_to_x3d')
+        from lxml import etree
+
         dim = self.dim
         # Create points and vector components as lists:
         zz, yy, xx = np.mgrid[0.5:(dim[0]-0.5):dim[0]*1j,

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__')
@@ -83,9 +83,10 @@ 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`.
+    Blazingly fast for l2-based cost functions.
 
     Parameters
     ----------
@@ -93,19 +94,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,21 +104,23 @@ def optimize_linear(data, regularisator=None, maxiter=1000, verbosity=0):
         The reconstructed magnetic distribution as a :class:`~.MagData` object.
 
     '''
-    LOG.debug('Calling optimize_sparse_cg')
+    import jutil.cg as jcg
+    LOG.debug('Calling optimize_linear')
     # Set up necessary objects:
     cost = Costfunction(data, regularisator)
-    print cost(np.zeros(cost.n))
-    x_opt = 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 = 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
     ----------
@@ -135,10 +128,10 @@ 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.
 
     Returns
     -------
@@ -146,40 +139,94 @@ def optimize_nonlin(data, first_guess=None, regularisator=None):
         The reconstructed magnetic distribution as a :class:`~.MagData` object.
 
     '''
-    LOG.debug('Calling optimize_cg')
+    import jutil.minimizer as jmin
+    import jutil.norms as jnorms
+    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)
-
-#    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
-
-#    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
-
-#   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
-
     assert len(x_0) == cost.n, (len(x_0), cost.m, cost.n)
-    result = minimizer.minimize(cost, x_0, options={"conv_rel": 1e-2}, tol={"max_iteration": 4})
+
+    p = regularisator.p
+    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 Lp 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
-    print cost(x_opt)
-    mag_opt = MagData(data.a, np.zeros((3,)+data.dim))
+    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
+
+def optimize_splitbregman(data, weight, lam, mu):
+    '''
+    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
+    ----------
+    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.
+    weight : float
+        Obscure split bregman parameter
+    lam : float
+        Cryptic split bregman parameter
+    mu : float
+        flabberghasting split bregman paramter
+
+    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 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))
+    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=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
 
@@ -221,7 +268,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):
@@ -256,6 +303,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

pyramid/regularisator.py

@@ -12,7 +12,8 @@ import numpy as np
 
 from scipy.sparse import coo_matrix, csr_matrix
 import jutil.norms as jnorm
-
+import jutil.diff as jdiff
+import jutil.operator as joperator
 from pyramid.converter import IndexConverter
 
 import logging
@@ -103,9 +104,13 @@ class ZeroOrderRegularisator(Regularisator):
 
     LOG = logging.getLogger(__name__+'.ZeroOrderRegularisator')
 
-    def __init__(self, lam):
+    def __init__(self, _, lam, p=2):
         self.LOG.debug('Calling __init__')
-        norm = jnorm.L2Square()
+        self.p = p
+        if p == 2:
+            norm = jnorm.L2Square()
+        else:
+            norm = jnorm.LPPow(p, 1e-12)
         super(ZeroOrderRegularisator, self).__init__(norm, lam)
         self.LOG.debug('Created '+str(self))
 
@@ -113,11 +118,14 @@ class ZeroOrderRegularisator(Regularisator):
 class FirstOrderRegularisator(Regularisator):
     # TODO: Docstring!
 
-    def __init__(self, mask, lam, x_a=None):
-        import jutil
-        D0 = jutil.diff.get_diff_operator(mask, 0, 3)
-        D1 = jutil.diff.get_diff_operator(mask, 1, 3)
-        D = jutil.operator.VStack([D0, D1])
-        norm = jutil.norms.WeightedL2Square(D)
+    def __init__(self, mask, lam, p=2):
+        self.p = p
+        D0 = jdiff.get_diff_operator(mask, 0, 3)
+        D1 = jdiff.get_diff_operator(mask, 1, 3)
+        D = joperator.VStack([D0, D1])
+        if p == 2:
+            norm = jnorm.WeightedL2Square(D)
+        else:
+            norm = jnorm.WeightedTV(jnorm.LPPow(p, 1e-12), D, [D0.shape[0], D.shape[0]])
         super(FirstOrderRegularisator, self).__init__(norm, lam)
         self.LOG.debug('Created '+str(self))

scripts/collaborations/example_joern.py

@@ -7,7 +7,8 @@ import os
 import numpy as np
 
 import pickle
-
+import matplotlib
+matplotlib.use("Agg")
 import matplotlib.pyplot as plt
 
 import pyramid
@@ -17,18 +18,44 @@ from pyramid.phasemapper import pm
 from pyramid.dataset import DataSet
 from pyramid.regularisator import ZeroOrderRegularisator, FirstOrderRegularisator
 import pyramid.reconstruction as rc
-
 from time import clock
 
 import logging
 import logging.config
 
+if "PBS_ARRAYID" in os.environ:
+    job_number = int(os.getenv("PBS_ARRAYID"))
+    experiments = []
+    for lam in [1e-9, 1e-8, 1e-7, 1e-6, 5e-6, 1e-5, 5e-5, 1e-4, 5e-4, 1e-3, 5e-3, 1e-2, 5e-2, 1e-1, 5e-1, 1]:
+        for order in [0, 1]:
+            for p in [2, 1.5, 1.1, 1.01, 1.001]:
+                for use_mask in [True, False]:
+                    experiments.append(
+                        (lam, order, p, use_mask))
+    assert 0 <= job_number < len(experiments), \
+        "job id too large, maximum={}".format(len(experiments) - 1)
+    lam, order, p, use_mask = experiments[job_number]
+    print experiments[job_number]
+    dirname = "new-order_{}-p_{}-mask_{}-lambda_{:.9f}-job{}/".format(
+        order, p, use_mask, lam, job_number)
+    if os.path.exists(dirname):
+        print dirname
+#        exit()
+    else:
+        os.mkdir(dirname)
+else:
+    p = 2
+    lam = 1e-8
+    use_mask = True
+    order = 1
+    dirname = "./"
 
-LOGGING_CONF = os.path.join(os.path.dirname(os.path.realpath(pyramid.__file__)), 'logging.ini')
 
+LOGGING_CONF = os.path.join(os.path.dirname(os.path.realpath(pyramid.__file__)), 'logging.ini')
 
 logging.config.fileConfig(LOGGING_CONF, disable_existing_loggers=False)
 logging.basicConfig(level=logging.INFO)
+
 ###################################################################################################
 threshold = 1
 a = 1.0  # in nm
@@ -39,10 +66,8 @@ dim = (1,) + (64, 64)
 dim_small = (64, 64)
 smoothed_pictures = True
 lam = 1E-6
-order = 1
 log = True
 PATH = '../../output/joern/'
-dirname = PATH
 ###################################################################################################
 # Read in files:
 phase_map = PhaseMap.load_from_netcdf4(PATH+'phase_map.nc')
@@ -50,16 +75,29 @@ phase_map = PhaseMap.load_from_netcdf4(PATH+'phase_map.nc')
 with open(PATH + 'mask.pickle') as pf:
     mask = pickle.load(pf)
 # Setup:
+if not use_mask:
+    mask = np.ones_like(mask, dtype=bool)
 data_set = DataSet(a, dim, b_0, mask=mask)
 data_set.append(phase_map, SimpleProjector(dim))
-regularisator = ZeroOrderRegularisator(lam)
-regularisator = FirstOrderRegularisator(mask, lam)
-print "OOO"
+if order == 0:
+    regularisator = ZeroOrderRegularisator(mask, lam, p)
+elif order == 1:
+    regularisator = FirstOrderRegularisator(mask, lam, p)
+else:
+    raise NotImplementedError
 
 # Reconstruct the magnetic distribution:
 tic = clock()
-mag_data_rec = rc.optimize_linear(data_set, regularisator=regularisator)
-#mag_data_rec = rc.optimize_nonlin(data_set, regularisator=regularisator)
+
+if p == 2:
+    mag_data_rec = rc.optimize_linear(data_set, regularisator=regularisator)
+else:
+    print regularisator.p
+    mag_data_rec = rc.optimize_nonlin(data_set, regularisator=regularisator)
+
+with open(dirname + "result.pickle", "wb") as pf:
+    import cPickle
+    cPickle.dump(mag_data_rec, pf)
 #  .optimize_simple_leastsq(phase_map, mask, b_0, lam=lam, order=order)
 
 print 'reconstruction time:', clock() - tic