Added basic framework for charge reconstruction.

pyramid/dataset.py

@@ -140,7 +140,7 @@ class DataSet(object):
             'Dimension has to be a tuple of length 3!'
         self.a = a
         self.dim = dim
-        self.b_0 = b_0
+        self.b_0 = b_0  # TODO: Not very general!!! get rid off!
         self.mask = mask
         self.Se_inv = Se_inv
         self._phasemaps = []
@@ -169,6 +169,8 @@ class DataSet(object):
         # Create lookup key:
         # TODO: Think again if phasemappers should be given as attribute (seems to be faulty
         # TODO: currently... Also not very expensive, so keep outside?
+        # TODO: Streamline DataSet, only container, not so much work inside
+        # TODO: Generalise for arbitrary sets of phasemaps, independent of magnetic/charge/ramp!
         if phasemapper is not None:
             key = dim_uv  # Create standard phasemapper, dim_uv is enough for identification!
         else:

pyramid/fielddata.py

@@ -939,7 +939,7 @@ class VectorData(FieldData):
         else:
             cbar_mappable, cbar_label = None, None
         # Change background color:
-        axis.set_facecolor(bgcolor)
+        #axis.set_facecolor(bgcolor)  # TODO: Activate for matplotlib 2.0!
         # Show mask:
         if show_mask and not np.all(submask):  # Plot mask if desired and not trivial!
             vv, uu = np.indices(dim_uv) + 0.5  # shift to center of pixel
@@ -1042,7 +1042,7 @@ class VectorData(FieldData):
                     extent=(0, dim_uv[1], 0, dim_uv[0]))
         # Change background color:
         if bgcolor is not None:
-            axis.set_facecolor(bgcolor)
+            pass#axis.set_facecolor(bgcolor)  # TODO: Activate for matplotlib 2.0!
         # Show mask:
         if show_mask and not np.all(submask):  # Plot mask if desired and not trivial!
             vv, uu = np.indices(dim_uv) + 0.5  # shift to center of pixel

pyramid/forwardmodel.py

@@ -18,6 +18,10 @@ from pyramid.ramp import Ramp
 __all__ = ['ForwardModel', 'DistributedForwardModel']
 
 
+# TODO: Ramp should be a forward model itself! Instead of hookpoints, each ForwardModel should
+# TODO: keep track of start and end in vector x (defaults: 0 and -1)!
+# TODO: Write CombinedForwardModel class!
+
 class ForwardModel(object):
     """Class for mapping 3D magnetic distributions to 2D phase maps.
 

pyramid/kernel.py

@@ -106,8 +106,8 @@ class Kernel(object):
         v = np.linspace(-(v_dim - 1), v_dim - 1, num=2 * v_dim - 1)
         uu, vv = np.meshgrid(u, v)
         # TODO: u, v are coordinates, rename self.u/v to self.kern_u/v!
-        self.u = np.empty(self.dim_kern, dtype=dtype)
-        self.v = np.empty(self.dim_kern, dtype=dtype)
+        self.u = np.empty(self.dim_kern, dtype=dtype)  # TODO: Is this necessary?
+        self.v = np.empty(self.dim_kern, dtype=dtype)  # TODO: Move to _get_elementary_phase?
         self.u[...] = coeff * self._get_elementary_phase(geometry, uu, vv, a)
         # TODO: The minus sign belongs into the phasemapper (debatable)!
         self.v[...] = coeff * -self._get_elementary_phase(geometry, vv, uu, a)

pyramid/plottools.py

@@ -169,7 +169,7 @@ def add_colorwheel(axis):
     inset_axes = inset_axes(axis, width=0.75, height=0.75, loc=1)
     inset_axes.axis('off')
     cmap = colors.CMAP_CIRCULAR_DEFAULT
-    bgcolor = axis.get_facecolor()
+    bgcolor = None#axis.get_facecolor() TODO: Activate for matplotlib 2.0!
     return cmap.plot_colorwheel(size=100, axis=inset_axes, alpha=0, bgcolor=bgcolor, arrows=True)
 
 

pyramid/reconstruction.py

@@ -77,6 +77,61 @@ def optimize_linear(costfunction, mag_0=None, ramp_0=None, max_iter=None, verbos
     return mag_opt
 
 
+def optimize_linear_charge(costfunction, charge_0=None, ramp_0=None, max_iter=None, verbose=False):
+    """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
+    ----------
+    costfunction : :class:`~.Costfunction`
+        A :class:`~.Costfunction` object which implements a specified forward model and
+        regularisator which is minimized in the optimization process.
+    mag_0: :class:`~.VectorData`
+        The starting magnetisation distribution used for the reconstruction. A zero vector will be
+        used if no VectorData object is specified.
+    mag_0: :class:`~.Ramp`
+        The starting ramp for the reconstruction. A zero vector will be
+        used if no Ramp object is specified.
+    max_iter : int, optional
+        The maximum number of iterations for the opimization.
+    verbose: bool, optional
+        If set to True, information like a progressbar is displayed during reconstruction.
+        The default is False.
+
+    Returns
+    -------
+    magdata : :class:`~pyramid.fielddata.VectorData`
+        The reconstructed magnetic distribution as a :class:`~.VectorData` object.
+
+    """
+    import jutil.cg as jcg
+    from jutil.taketime import TakeTime
+    _log.debug('Calling optimize_linear')
+    _log.info('Cost before optimization: {:.3e}'.format(costfunction(np.zeros(costfunction.n))))
+    data_set = costfunction.fwd_model.data_set
+    # Get starting distribution vector x_0:
+    x_0 = np.empty(costfunction.n)
+    if charge_0 is not None:
+        costfunction.fwd_model.magdata = charge_0
+    x_0[:data_set.n] = costfunction.fwd_model.magdata.get_vector(mask=data_set.mask)
+    if ramp_0 is not None:
+        ramp_vec = ramp_0.param_cache.ravel()
+    else:
+        ramp_vec = np.zeros_like(costfunction.fwd_model.ramp.n)
+    x_0[data_set.n:] = ramp_vec
+    # Minimize:
+    with TakeTime('reconstruction time'):
+        x_opt = jcg.conj_grad_minimize(costfunction, x_0=x_0, max_iter=max_iter, verbose=verbose).x
+    _log.info('Cost after optimization: {:.3e}'.format(costfunction(x_opt)))
+    # Cut ramp parameters if necessary (this also saves the final parameters in the ramp class!):
+    x_opt = costfunction.fwd_model.ramp.extract_ramp_params(x_opt)
+    # Create and return fitting VectorData object:
+    mag_opt = VectorData(data_set.a, np.zeros((3,) + data_set.dim))
+    mag_opt.set_vector(x_opt, data_set.mask)
+    return mag_opt
+
+
 def optimize_nonlin(costfunction, first_guess=None):
     """Reconstruct a three-dimensional magnetic distribution from given phase maps via
     steepest descent method. This is slow, but works best for non l2-regularisators.

pyramid/utils/reconstruction_2d_from_phasemap.py

@@ -11,7 +11,9 @@ import numpy as np
 from .. import reconstruction
 from ..dataset import DataSet
 from ..projector import SimpleProjector
-from ..regularisator import FirstOrderRegularisator
+from ..regularisator import FirstOrderRegularisator, NoneRegularisator
+from ..kernel import KernelCharge
+from ..phasemapper import PhaseMapperCharge
 from ..forwardmodel import ForwardModel
 from ..costfunction import Costfunction
 from .pm import pm
@@ -105,3 +107,94 @@ def reconstruction_2d_from_phasemap(phasemap, b_0=1, lam=1E-3, max_iter=100, ram
             ramp.plot_phase(note='Fitted Ramp')
     # Return reconstructed magnetisation distribution and cost function:
     return magdata_rec, cost
+
+
+
+def reconstruction_2d_charge_from_phasemap(phasemap, lam=1E-3, max_iter=100, ramp_order=1,
+                                           electrode_vec=(1E6, 1E6), v_acc=300E3,
+                                           plot_results=False, verbose=True):
+    """Convenience function for reconstructing a projected distribution from a single phasemap.
+
+    Parameters
+    ----------
+    phasemap: :class:`~PhaseMap`
+        The phasemap which is used for the reconstruction.
+    b_0 : float, optional
+        The magnetic induction corresponding to a magnetization `M`\ :sub:`0` in T.
+        The default is 1.
+    lam : float
+        Regularisation parameter determining the weighting between measurements and regularisation.
+    max_iter : int, optional
+        The maximum number of iterations for the opimization.
+    ramp_order : int or None (default)
+        Polynomial order of the additional phase ramp which will be added to the phase maps.
+        All ramp parameters have to be at the end of the input vector and are split automatically.
+        Default is None (no ramps are added).
+    plot_results: boolean, optional
+        If True, the results are plotted after reconstruction.
+    ar_dens: int, optional
+        Number defining the arrow density which is plotted. A higher ar_dens number skips more
+        arrows (a number of 2 plots every second arrow). Default is 1.
+    verbose: bool, optional
+        If set to True, information like a progressbar is displayed during reconstruction.
+        The default is False.
+
+    Returns
+    -------
+    magdata_rec, cost: :class:`~.VectorData`, :class:`~.Costfunction`
+        The reconstructed magnetisation distribution and the used costfunction.
+
+    """
+    _log.debug('Calling reconstruction_2d_from_phasemap')
+    # Construct DataSet, Regularisator, ForwardModel and Costfunction:
+    dim = (1,) + phasemap.dim_uv
+    data = DataSet(phasemap.a, dim, b_0=1)
+    kernel = KernelCharge(phasemap.a, phasemap.dim_uv, electrode_vec=electrode_vec, v_acc=v_acc)
+    data.append(phasemap, SimpleProjector(dim), PhaseMapperCharge(kernel))
+    data.set_3d_mask()
+    # TODO: Rework classes below (ForwardModel, Costfunction)!
+    fwd_model = ForwardModel(data, ramp_order)
+    reg = NoneRegularisator()#FirstOrderRegularisator(data.mask, lam, add_params=fwd_model.ramp.n)
+    cost = Costfunction(fwd_model, reg)
+    # Reconstruct:
+    magdata_rec = reconstruction.optimize_linear(cost, max_iter=max_iter, verbose=verbose)
+    param_cache = cost.fwd_model.ramp.param_cache
+    if ramp_order is None:
+        offset, ramp = 0, (0, 0)
+    elif ramp_order >= 1:
+        offset, ramp = param_cache[0][0], (param_cache[1][0], param_cache[2][0])
+    elif ramp_order == 0:
+        offset, ramp = param_cache[0][0], (0, 0)
+    else:
+        raise ValueError('ramp_order has to be a positive integer or None!')
+    # Plot stuff:
+    if plot_results:
+        if ar_dens is None:
+            ar_dens = np.max([1, np.max(dim) // 64])
+        magdata_rec.plot_quiver_field(note='Reconstructed Distribution',
+                                      ar_dens=ar_dens, figsize=(16, 16))
+        phasemap_rec = pm(magdata_rec)
+        gain = 4 * 2 * np.pi / (np.abs(phasemap_rec.phase).max() + 1E-30)
+        gain = round(gain, -int(np.floor(np.log10(abs(gain)))))
+        vmin = phasemap_rec.phase.min()
+        vmax = phasemap_rec.phase.max()
+        phasemap.plot_combined(note='Input Phase', gain=gain)
+        phasemap -= fwd_model.ramp(index=0)
+        phasemap.plot_combined(note='Input Phase (ramp corrected)', gain=gain, vmin=vmin, vmax=vmax)
+        title = 'Reconstructed Phase'
+        if ramp_order is not None:
+            if ramp_order >= 0:
+                print('offset:', offset)
+                # title += ', fitted Offset: {:.2g} [rad]'.format(offset)
+            if ramp_order >= 1:
+                print('ramp:', ramp)
+                # title += ', (Fitted Ramp: (u:{:.2g}, v:{:.2g}) [rad/nm]'.format(*ramp)
+        phasemap_rec.plot_combined(note=title, gain=gain, vmin=vmin, vmax=vmax)
+        diff = (phasemap_rec - phasemap)
+        diff_name = 'Difference (RMS: {:.2g} rad)'.format(np.sqrt(np.mean(diff.phase) ** 2))
+        diff.plot_phase_with_hist(note=diff_name, sigma_clip=3)
+        if ramp_order is not None:
+            ramp = fwd_model.ramp(0)
+            ramp.plot_phase(note='Fitted Ramp')
+    # Return reconstructed magnetisation distribution and cost function:
+    return magdata_rec, cost

pyramid/utils/reconstruction_3d_from_magdata.py

@@ -117,7 +117,7 @@ def reconstruction_3d_from_magdata(magdata, b_0=1, lam=1E-3, max_iter=100, ramp_
         phasemap += Ramp.create_ramp(phasemap.a, phasemap.dim_uv, (offset, ramp_u, ramp_v))
         data.phasemaps[i] = phasemap
     # Add noise if necessary:
-    if noise != 0:
+    if noise != 0:  # TODO: write function to add noise after APERTURE!! (ask Florian again)
         for i, phasemap in enumerate(data.phasemaps):
             phasemap.phase += np.random.normal(0, noise, phasemap.dim_uv)
             data.phasemaps[i] = phasemap