From 8307dcbbad84d07d31206d7ee7ada483bc9c4284 Mon Sep 17 00:00:00 2001
From: Jan Caron <j.caron@fz-juelich.de>
Date: Fri, 18 Aug 2017 15:28:58 +0200
Subject: [PATCH] Added basic framework for charge reconstruction.
---
pyramid/dataset.py | 4 +-
pyramid/fielddata.py | 4 +-
pyramid/forwardmodel.py | 4 +
pyramid/kernel.py | 4 +-
pyramid/plottools.py | 2 +-
pyramid/reconstruction.py | 55 +++++++++++
.../utils/reconstruction_2d_from_phasemap.py | 95 ++++++++++++++++++-
.../utils/reconstruction_3d_from_magdata.py | 2 +-
8 files changed, 162 insertions(+), 8 deletions(-)
diff --git a/pyramid/dataset.py b/pyramid/dataset.py
index 408ac18..188f09b 100644
--- a/pyramid/dataset.py
+++ b/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:
diff --git a/pyramid/fielddata.py b/pyramid/fielddata.py
index f295b48..76dd3e7 100644
--- a/pyramid/fielddata.py
+++ b/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
diff --git a/pyramid/forwardmodel.py b/pyramid/forwardmodel.py
index 17e1d98..9febbb4 100644
--- a/pyramid/forwardmodel.py
+++ b/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.
diff --git a/pyramid/kernel.py b/pyramid/kernel.py
index 8f58314..dd6ef5a 100644
--- a/pyramid/kernel.py
+++ b/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)
diff --git a/pyramid/plottools.py b/pyramid/plottools.py
index e76ad26..635af1c 100644
--- a/pyramid/plottools.py
+++ b/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)
diff --git a/pyramid/reconstruction.py b/pyramid/reconstruction.py
index cc31372..9ebec34 100644
--- a/pyramid/reconstruction.py
+++ b/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.
diff --git a/pyramid/utils/reconstruction_2d_from_phasemap.py b/pyramid/utils/reconstruction_2d_from_phasemap.py
index 14ad686..667717b 100644
--- a/pyramid/utils/reconstruction_2d_from_phasemap.py
+++ b/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
diff --git a/pyramid/utils/reconstruction_3d_from_magdata.py b/pyramid/utils/reconstruction_3d_from_magdata.py
index fcb649d..cd930a4 100644
--- a/pyramid/utils/reconstruction_3d_from_magdata.py
+++ b/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
--
GitLab