From 05a75ec9ea63e26b1ea44248621e3a7f51d98c49 Mon Sep 17 00:00:00 2001
From: Fengshan Zheng <f.zheng@fz-juelich.de>
Date: Thu, 24 Oct 2019 14:48:49 +0200
Subject: [PATCH] Q3 d
---
pyramid/dataset.py | 1 +
pyramid/diagnostics.py | 216 +++++++++++++-----
pyramid/fielddata.py | 29 +++
pyramid/kernel.py | 9 +
pyramid/phasemapper.py | 2 +-
pyramid/ramp.py | 2 +-
pyramid/reconstruction.py | 6 +-
pyramid/regularisator.py | 12 +-
pyramid/utils/__init__.py | 8 +-
pyramid/utils/pm.py | 15 +-
.../utils/reconstruction_2d_from_phasemap.py | 65 +++---
.../utils/reconstruction_3d_from_magdata.py | 151 +++++++++++-
12 files changed, 406 insertions(+), 110 deletions(-)
diff --git a/pyramid/dataset.py b/pyramid/dataset.py
index 2c95343..fbc161d 100644
--- a/pyramid/dataset.py
+++ b/pyramid/dataset.py
@@ -591,6 +591,7 @@ class DataSetCharge(object):
elif phasemapper is not None: # Use given one (do nothing):
pass
else: # Create new standard (RDFC) phasemapper:
+ # TODO: PRW is missing in the kernel
phasemapper = PhaseMapperCharge(KernelCharge(self.a, dim_uv, self.electrode_vec))
self._phasemapper_dict[key] = phasemapper
# Append everything to the lists (just contain pointers to objects!):
diff --git a/pyramid/diagnostics.py b/pyramid/diagnostics.py
index 033eb34..50a5226 100644
--- a/pyramid/diagnostics.py
+++ b/pyramid/diagnostics.py
@@ -11,10 +11,10 @@ import logging
import pickle
-from pyramid.forwardmodel import ForwardModel
+from pyramid.forwardmodel import ForwardModel, ForwardModelCharge
from pyramid.costfunction import Costfunction
-from pyramid.regularisator import FirstOrderRegularisator
-from pyramid.fielddata import VectorData
+from pyramid.regularisator import FirstOrderRegularisator, ZeroOrderRegularisator
+from pyramid.fielddata import VectorData, ScalarData
from pyramid.phasemap import PhaseMap
from pyramid import reconstruction
from pyramid import plottools
@@ -36,11 +36,10 @@ try:
except NameError:
from tqdm import tqdm
-__all__ = ['Diagnostics', 'LCurve', 'get_vector_field_errors']
+__all__ = ['Diagnostics', 'LCurve', 'LCurveCharge', 'get_vector_field_errors']
# TODO: should be subpackage, distribute methods and classes to separate modules!
-
class Diagnostics(object):
"""Class for calculating diagnostic properties of a specified costfunction.
@@ -154,10 +153,10 @@ class Diagnostics(object):
self._updated_avrg_kern_row = False
self._updated_measure_contribution = False
- def __init__(self, cost, max_iter=1000, verbose=False): # TODO: verbose True default
+ def __init__(self, magdata, cost, max_iter=1000, verbose=False): # TODO: verbose True default
self._log.debug('Calling __init__')
+ self.magdata = magdata
self.cost = cost
- self.a = self.cost.fwd_model.data_set.a
self.max_iter = max_iter
self.verbose = verbose
self.fwd_model = cost.fwd_model
@@ -165,8 +164,8 @@ class Diagnostics(object):
self.dim = cost.fwd_model.data_set.dim
self.mask = cost.fwd_model.data_set.mask
self.x_rec = np.empty(cost.n)
- # self.x_rec[:self.fwd_model.data_set.n] = self.magdata.get_vector(mask=self.mask)
- # self.x_rec[self.fwd_model.data_set.n:] = self.fwd_model.ramp.param_cache.ravel()
+ self.x_rec[:self.fwd_model.data_set.n] = self.magdata.get_vector(mask=self.mask)
+ self.x_rec[self.fwd_model.data_set.n:] = self.fwd_model.ramp.param_cache.ravel()
self.row_idx = None
self.pos = (0,) + tuple(np.array(np.where(self.mask))[:, 0]) # first True mask entry
self._updated_cov_row = False
@@ -195,9 +194,7 @@ class Diagnostics(object):
if pos is not None:
self.pos = pos
magdata_avrg_kern = VectorData(self.cost.fwd_model.data_set.a, np.zeros((3,) + self.dim))
- vector = self.avrg_kern_row
- if self.fwd_model.ramp.order is not None:
- vector = vector[:-self.fwd_model.ramp.n] # Only take vector field, not ramp!
+ vector = self.avrg_kern_row[:-self.fwd_model.ramp.n] # Only take vector field, not ramp!
magdata_avrg_kern.set_vector(vector, mask=self.mask)
return magdata_avrg_kern
@@ -229,6 +226,7 @@ class Diagnostics(object):
"""
self._log.debug('Calling calculate_fwhm')
+ a = self.magdata.a
magdata_avrg_kern = self.get_avrg_kern_field(pos)
x = np.arange(0, self.dim[2]) - self.pos[3]
y = np.arange(0, self.dim[1]) - self.pos[2]
@@ -270,9 +268,9 @@ class Diagnostics(object):
lz, rz = _calc_lr(2)
# TODO: Test if FWHM is really calculated with a in mind... didn't seem so...
- fwhm_x = (rx - lx) * self.a
- fwhm_y = (ry - ly) * self.a
- fwhm_z = (rz - lz) * self.a
+ fwhm_x = (rx - lx) * a
+ fwhm_y = (ry - ly) * a
+ fwhm_z = (rz - lz) * a
# Plot helpful stuff:
if plot:
fig, axis = plt.subplots(1, 1)
@@ -303,8 +301,7 @@ class Diagnostics(object):
axis.set_xlabel('x/y/z-slice [nm]', fontsize=15)
axis.set_ylabel('information content [%]', fontsize=15)
axis.tick_params(axis='both', which='major', labelsize=14)
- formatter = FuncFormatter(lambda x, pos: '{:.3g}'.format(x * self.a))
- axis.xaxis.set_major_formatter(formatter)
+ axis.xaxis.set_major_formatter(FuncFormatter(lambda x, pos: '{:.3g}'.format(x * a)))
comp_legend = axis.legend([cx, cy, cz], [c.get_label() for c in [cx, cy, cz]], loc=2,
scatterpoints=1, prop={'size': 14})
axis.legend(l, [i.get_label() for i in l], loc=1, numpoints=1, prop={'size': 14})
@@ -344,7 +341,7 @@ class Diagnostics(object):
gain_map_list.append(PhaseMap(self.cost.fwd_model.data_set.a, gain))
return gain_map_list
- def plot_position(self, magdata, **kwargs):
+ def plot_position(self, **kwargs):
proj_axis = kwargs.get('proj_axis', 'z')
if proj_axis == 'z': # Slice of the xy-plane with z = ax_slice
pos_2d = (self.pos[2], self.pos[3])
@@ -362,7 +359,7 @@ class Diagnostics(object):
comp = {0: 'x', 1: 'y', 2: 'z'}[self.pos[0]]
note = '{}-comp., pos.: {}'.format(comp, self.pos[1:])
# Plots:
- axis = magdata.plot_quiver_field(note=note, ax_slice=ax_slice, **kwargs)
+ axis = self.magdata.plot_quiver_field(note=note, ax_slice=ax_slice, **kwargs)
rect = axis.add_patch(patches.Rectangle((pos_2d[1], pos_2d[0]), 1, 1, fill=False,
edgecolor='w', linewidth=2, alpha=0.5))
rect.set_path_effects([patheffects.withStroke(linewidth=4, foreground='k', alpha=0.5)])
@@ -371,21 +368,22 @@ class Diagnostics(object):
pass
def plot_avrg_kern_field(self, pos=None, **kwargs):
+ a = self.magdata.a
avrg_kern_field = self.get_avrg_kern_field(pos)
fwhms, lr = self.calculate_fwhm(pos)[:2]
proj_axis = kwargs.get('proj_axis', 'z')
if proj_axis == 'z': # Slice of the xy-plane with z = ax_slice
pos_2d = (self.pos[2], self.pos[3])
ax_slice = self.pos[1]
- width, height = fwhms[0] / self.a, fwhms[1] / self.a
+ width, height = fwhms[0] / a, fwhms[1] / a
elif proj_axis == 'y': # Slice of the xz-plane with y = ax_slice
pos_2d = (self.pos[1], self.pos[3])
ax_slice = self.pos[2]
- width, height = fwhms[0] / self.a, fwhms[2] / self.a
+ width, height = fwhms[0] / a, fwhms[2] / a
elif proj_axis == 'x': # Slice of the zy-plane with x = ax_slice
pos_2d = (self.pos[2], self.pos[1])
ax_slice = self.pos[3]
- width, height = fwhms[2] / self.a, fwhms[1] / self.a
+ width, height = fwhms[2] / a, fwhms[1] / a
else:
raise ValueError('{} is not a valid argument (use x, y or z)'.format(proj_axis))
note = kwargs.pop('note', None)
@@ -435,9 +433,12 @@ class Diagnostics(object):
actor.actor.orientation = [0, 0, 0] # in degrees
actor.actor.origin = (0, 0, 0)
actor.actor.position = (self.pos[1]+0.5, self.pos[2]+0.5, self.pos[3]+0.5)
- actor.actor.scale = [0.5*fwhm[0]/self.a, 0.5*fwhm[1]/self.a, 0.5*fwhm[2]/self.a]
+ a = self.magdata.a
+ actor.actor.scale = [0.5*fwhm[0]/a, 0.5*fwhm[1]/a, 0.5*fwhm[2]/a]
+
#surface.append(surface)
+
scene.scene.disable_render = False # now turn it on # TODO: EVERYWHERE WITH MAYAVI!
@@ -590,9 +591,9 @@ class LCurve(object):
self._log.debug('axis is None')
fig = plt.figure(figsize=figsize)
axis = fig.add_subplot(1, 1, 1)
- axis.set_yscale("log", nonposy='clip')
+ axis.set_yscale("log", nonposx='clip')
axis.set_xscale("log", nonposx='clip')
- axis.plot(x, y, 'grey', linestyle='-', linewidth=3, zorder=1)
+ axis.plot(x, y, 'k-', linewidth=3, zorder=1)
sc = axis.scatter(x, y, c=lambdas, marker='o', s=100, zorder=2,
cmap='nipy_spectral', norm=LogNorm())
plt.colorbar(mappable=sc, label='regularisation parameter $\lambda$')
@@ -607,39 +608,148 @@ class LCurve(object):
axis.grid()
return axis
# TODO: Don't plot the steep part on the right...
- # TODO: The following formatting is better but without usetex, \boldsymbol doesn't work!
- # import matplotlib
- # matplotlib.rc('text', usetex=True)
- # matplotlib.rcParams['text.latex.preamble'] = [r"\usepackage{amsmath}"]
- # axis.get_xaxis().get_major_formatter().labelOnlyBase = False
- # axis.set_xlabel(
- # r'$\Vert\mathbf{F}\boldsymbol{x}-\boldsymbol{y}\Vert_{\mathbf{S}_{\epsilon}^{-1}}^{2}$',
- # fontsize=22)
- # axis.set_ylabel(r'$\frac{1}{\lambda}\Vert\boldsymbol{x}\Vert_{\mathbf{S}_{a}^{-1}}^{2}$',
- # fontsize=22)
- # matplotlib.rc('text', usetex=False)
-
-
-def get_vector_field_errors(vector_data, vector_data_ref, mask=None):
+
+
+class LCurveCharge(object):
+
+ # TODO: Docstring!
+
+ # TODO: save elecdata_rec!
+
+ _log = logging.getLogger(__name__ + '.FieldData')
+
+ def __init__(self, fwd_model, max_iter=1000, regularisator=0, verbose=True, save_dir='lcurve'):
+ self._log.debug('Calling __init__')
+ assert isinstance(fwd_model, ForwardModelCharge), 'Input has to be a costfunction'
+ self.fwd_model = fwd_model
+ self.max_iter = max_iter
+ self.verbose = verbose
+ self.regularisator = regularisator
+ self.l_dict = {}
+ self.save_dir = save_dir
+ if self.save_dir is not None:
+ if not os.path.isdir(self.save_dir): # Create directory if it does not exist:
+ os.makedirs(self.save_dir)
+ if os.path.isfile('{}/lcurve.pkl'.format(self.save_dir)): # Load file if it exists:
+ self._load()
+ else: # Create file:
+ self._save()
+ self._log.debug('Created ' + str(self))
+
+ # TODO: Methods for saving and loading l_dict's!!!
+ def _save(self):
+ with open('{}/lcurve.pkl'.format(self.save_dir), 'wb') as f:
+ pickle.dump(self.l_dict, f, pickle.HIGHEST_PROTOCOL)
+
+ def _load(self):
+ with open('{}/lcurve.pkl'.format(self.save_dir), 'rb') as f:
+ self.l_dict = pickle.load(f)
+
+ def calculate(self, lambdas, overwrite=False):
+ # TODO: Docstring!
+ lams = np.atleast_1d(lambdas)
+ for lam in tqdm(lams, disable=not self.verbose):
+ if lam not in self.l_dict.keys() or overwrite:
+ # Create new regularisator and costfunction: # TODO: Not hardcoding FirstOrder!
+ # TODO: Not necessary if lambda can be extracted from regularisator? self.cost?
+ if self.regularisator is 1:
+ reg = FirstOrderRegularisator(self.fwd_model.data_set.mask, lam,
+ add_params=self.fwd_model.ramp.n, factor=1)
+ else:
+ reg = ZeroOrderRegularisator(self.fwd_model.data_set.mask, lam, add_params=self.fwd_model.ramp.n)
+ cost = Costfunction(fwd_model=self.fwd_model, regularisator=reg)
+ # Reconstruct:
+ elecdata_rec = reconstruction.optimize_linear_charge(cost, max_iter=self.max_iter, verbose=self.verbose)
+ # Add new values to dictionary:
+ chisq_m, chisq_a = cost.chisq_m[-1], cost.chisq_a[-1] # TODO: chisq_m list or not?
+ self.l_dict[lam] = (chisq_m, chisq_a)
+ self._log.info(lam, ' --> m:', chisq_m, ' a:', chisq_a)
+ # Save elecdata_rec and dictionary if necessary:
+ if self.save_dir is not None:
+ filename = 'elecdata_rec_lam{:.0e}.hdf5'.format(lam)
+ elecdata_rec.save(os.path.join(self.save_dir, filename), overwrite=True)
+ self._save()
+
+ def calculate_auto(self, lam_start=1E-18, lam_end=1E5, online_axis=False):
+ raise NotImplementedError()
+ # TODO: Docstring!
+ # TODO: IMPLEMENT!!!
+ # # Calculate new cost terms:
+ # log_m_s, log_a_s = np.log10(self.calculate(lam_start))
+ # log_m_e, log_a_e = np.log10(self.calculate(lam_end))
+ # # Calculate new lambda:
+ # log_lam_s, log_lam_e = np.log10(lam_start), np.log10(lam_end)
+ # log_lam_new = np.mean((log_lam_s, log_lam_e)) # logarithmic mean to find middle on L!
+ # sign_exp = np.floor(log_lam_new)
+ # last_sign_digit = np.round(10 ** (log_lam_new - sign_exp))
+ # lam_new = last_sign_digit * 10 ** sign_exp
+ # # Calculate cost terms for new lambda:
+ # log_m_new, log_a_new = np.log10(self.calculate(lam_new))
+ # if online_axis: # Update plot if necessary:
+ # self.plot(axis=online_axis)
+ # from IPython import display
+ # display.clear_output(wait=True)
+ # display.display(plt.gcf())
+ # # Calculate distances from origin and find new interval:
+ # dist_s, dist_e = np.hypot(log_m_s, log_a_s), np.hypot(log_m_e, log_a_e)
+ # dist_new = np.hypot(log_m_new, log_a_new)
+ # print(lam_start, lam_end, lam_new)
+ # print(dist_s, dist_e, dist_new)
+ # # if dist_new
+ # TODO: slope has to be normalised, scale of axes is not equal!!!
+ # TODO: get rid of right flank (do Not use right points with slope steeper than -45°
+ # TODO: Implement else, return saved values!
+ # TODO: Make this work with batch, sort lambdas at the end!
+ # TODO: After sorting, calculate the CURVATURE for each lambda! (needs 3 points?)
+ # TODO: Use finite difference methods (forward/backward/central, depends on location)!
+ # TODO: Investigate points around highest curvature further.
+ # TODO: Make sure to update ALL curvatures and search for new best EVERYWHERE!
+ # TODO: Distinguish regions of the L-Curve.
+
+ def plot(self, lambdas=None, axis=None, figsize=None):
+ # TODO: Docstring!
+ # Sort lists according to lambdas:
+ if lambdas is None:
+ lambdas = sorted(self.l_dict.keys())
+ x, y = [], []
+ for lam in lambdas:
+ x.append(self.l_dict[lam][0])
+ y.append(self.l_dict[lam][1] / lam)
+ if figsize is None:
+ figsize = plottools.FIGSIZE_DEFAULT
+ if axis is None:
+ self._log.debug('axis is None')
+ fig = plt.figure(figsize=figsize)
+ axis = fig.add_subplot(1, 1, 1)
+ axis.set_yscale("log", nonposx='clip')
+ axis.set_xscale("log", nonposx='clip')
+ axis.plot(x, y, 'k-', linewidth=3, zorder=1)
+ sc = axis.scatter(x, y, c=lambdas, marker='o', s=100, zorder=2,
+ cmap='nipy_spectral', norm=LogNorm())
+ plt.colorbar(mappable=sc, label='regularisation parameter $\lambda$')
+ axis.set_xlabel(
+ r'$\Vert\mathbf{F}(\mathbf{x})-\mathbf{y}\Vert_{\mathbf{S}_{\epsilon}^{-1}}^{2}$',
+ fontsize=22, labelpad=-5)
+ axis.set_ylabel(r'$\frac{1}{\lambda}\Vert\mathbf{x}\Vert_{\mathbf{S}_{a}^{-1}}^{2}$',
+ fontsize=22)
+ axis.xaxis.label.set_color('firebrick')
+ axis.yaxis.label.set_color('seagreen')
+ axis.tick_params(axis='both', which='major')
+ axis.grid()
+ return axis
+ # TODO: Don't plot the steep part on the right...
+
+
+def get_vector_field_errors(vector_data, vector_data_ref):
"""After Kemp et. al.: Analysis of noise-induced errors in vector-field electron tomography"""
- if mask is not None:
- vector_data_masked = VectorData(vector_data.a, np.zeros(vector_data.shape))
- vector_data_masked.set_vector(vector_data.get_vector(mask), mask)
- vector_data_ref_masked = VectorData(vector_data_ref.a, np.zeros(vector_data_ref.shape))
- vector_data_ref_masked.set_vector(vector_data_ref.get_vector(mask), mask)
- v, vr = vector_data_masked.field, vector_data_ref_masked.field
- va, vra = vector_data_masked.field_amp, vector_data_ref_masked.field_amp
- volume = mask.sum()
- else:
- v, vr = vector_data.field, vector_data_ref.field
- va, vra = vector_data.field_amp, vector_data_ref.field_amp
- volume = np.prod(vector_data.dim)
+ v, vr = vector_data.field, vector_data_ref.field
+ va, vra = vector_data.field_amp, vector_data_ref.field_amp
+ volume = np.prod(vector_data.dim)
# Total error:
amp_sum_sqr = np.nansum((v - vr)**2)
rms_tot = np.sqrt(amp_sum_sqr / np.nansum(vra**2))
# Directional error:
- with np.errstate(divide='ignore', invalid='ignore'): # ignore "invalid value in true_divide"!
- scal_prod = np.clip(np.nansum(vr * v, axis=0) / (vra * va), -1, 1) # arccos float inacc.!
+ scal_prod = np.clip(np.nansum(vr * v, axis=0) / (vra * va), -1, 1) # arccos float pt. inacc.!
rms_dir = np.sqrt(np.nansum(np.arccos(scal_prod)**2) / volume) / np.pi
# Magnitude error:
rms_mag = np.sqrt(np.nansum((va - vra)**2) / np.nansum(vra**2))
diff --git a/pyramid/fielddata.py b/pyramid/fielddata.py
index 0ca31cb..2286d6a 100644
--- a/pyramid/fielddata.py
+++ b/pyramid/fielddata.py
@@ -1525,6 +1525,35 @@ class ScalarData(FieldData):
# TODO: flip!
+ def crop(self, crop_values): # TODO: it doesn't work now.
+ """Crop the current field distribution with zeros for each individual axis.
+
+ Parameters
+ ----------
+ crop_values : tuple of int
+ Number of zeros which should be cropped. Provided as a tuple where each entry
+ corresponds to an axis. An entry can be one int (same cropping for both sides) or again
+ a tuple which specifies the crop values for both sides of the corresponding axis.
+
+ Returns
+ -------
+ None
+
+ Notes
+ -----
+ Acts in place and changes dimensions accordingly.
+ """
+ self._log.debug('Calling crop')
+ assert len(crop_values) == 3, 'Crop values for each dimension have to be provided!'
+ cv = np.zeros(6, dtype=np.int)
+ for i, values in enumerate(crop_values):
+ assert np.shape(values) in [(), (2,)], 'Only one or two values per axis can be given!'
+ cv[2 * i:2 * (i + 1)] = values
+ cv *= np.resize([1, -1], len(cv))
+ cv = np.where(cv == 0, None, cv)
+ field_crop = self.field[cv[0]:cv[1], cv[2]:cv[3], cv[4]:cv[5]]
+ return ScalarData(self.a, field_crop)
+
def rotate(self, angle, axis='z', reshape=False, **kwargs):
# TODO: Docstring!
axes = {'x': (0, 1), 'y': (0, 2), 'z': (1, 2)}[axis] # Defines axes of plane of rotation!
diff --git a/pyramid/kernel.py b/pyramid/kernel.py
index c78aa4d..2a0df35 100644
--- a/pyramid/kernel.py
+++ b/pyramid/kernel.py
@@ -283,19 +283,28 @@ class KernelCharge(object):
# case 1 totally out of the sphere
case1 = ((h1 < 0) & (h2 < 0))
phase = np.zeros(n.shape)
+ # elec_field = np.zeros(n.shape)
phase[case1] += - np.log((r1[case1] ** 2) / (r2[case1] ** 2))
+ # The x, y component of electric field, Ex=Ey in the kernel
+ # elec_field[case1] += np.sqrt(3) * np.pi * (1 / r1[case1] - 1 / r2[case1]) / 3
# case 2: inside the charge sphere
case2 = ((h1 >= 0) & (h2 <= 0))
# The height when the path come across the charge sphere
h3 = np.sqrt(np.where(h1 > 0, h1, 0))
phase[case2] += (- np.log((h3[case2] + R) ** 2 / r2[case2] ** 2) +
(2 * h3[case2] / R + 2 * h3[case2] ** 3 / 3 / R ** 3))
+ # elec_field[case2] += 2 * np.sqrt(3) / 3 * (np.pi * (1 / r1[case1] - 1 / r2[case1]) / 2 -
+ # np.arctan(h3[case2] / r1[case2]) / r1[case2] + h3[case2] / 2 / R ** 2 +
+ # r1[case2] ** 2 / 2 / R ** 3 * (np.log(h3[case2] + R)) - np.log(r1[case2]))
# case 3 : inside the image charge sphere
case3 = np.logical_not(case1 | case2)
# The height whe the path comes across the image charge sphere
h4 = np.sqrt(np.where(h2 > 0, h2, 0))
phase[case3] += (np.log((h4[case3] + R) ** 2 / r1[case3] ** 2) -
(2 * h4[case3] / R + 2 * h4[case3] ** 3 / 3 / R ** 3))
+ # elec_field[case3] += 2 * np.sqrt(3) / 3 * (np.pi * (1 / r1[case1] - 1 / r2[case1]) / 2 -
+ # np.arctan(h4[case3] / r2[case3]) / r2[case3] - h4[case3] / 2 / R ** 2 -
+ # r2[case3] ** 2 / 2 / R ** 3 * (np.log(h4[case3] + R)) - np.log(r2[case3]))
return phase
def print_info(self):
diff --git a/pyramid/phasemapper.py b/pyramid/phasemapper.py
index 1a21690..68bee5a 100644
--- a/pyramid/phasemapper.py
+++ b/pyramid/phasemapper.py
@@ -402,7 +402,7 @@ class PhaseMapperMIP(PhaseMapper):
Returns
-------
result : :class:`~numpy.ndarray` (N=1)
- Product of the transposed Jacobi matrix (which is not explicitely calculated) with
+ Product of the transposed Jacobi matrix (which is not explicitly calculated) with
the vector, which has ``N**2`` entries like an electrostatic projection.
"""
diff --git a/pyramid/ramp.py b/pyramid/ramp.py
index 5b4f9b5..bc0e888 100644
--- a/pyramid/ramp.py
+++ b/pyramid/ramp.py
@@ -37,7 +37,7 @@ class Ramp(object):
(one for each degree of freedom) are saved along the first axis, values for different
images along the second axis.
n : int
- Size of the input space. Coincides with the numer of entries in `param_cache` and
+ Size of the input space. Coincides with the number of entries in `param_cache` and
calculates to ``deg_of_freedom * data_set.count``.
Notes
diff --git a/pyramid/reconstruction.py b/pyramid/reconstruction.py
index 96ee5b2..29d4799 100644
--- a/pyramid/reconstruction.py
+++ b/pyramid/reconstruction.py
@@ -87,9 +87,9 @@ def optimize_linear_charge(costfunction, charge_0=None, ramp_0=None, max_iter=No
costfunction : :class:`~.Costfunction`
A :class:`~.Costfunction` object which implements a specified forward model and
regularisator which is minimized in the optimization process.
- charge_0: :class:`~.VectorData`
- The starting magnetisation distribution used for the reconstruction. A zero vector will be
- used if no VectorData object is specified.
+ charge_0: :class:`~.ScalarData`
+ The starting charge distribution used for the reconstruction. A zero vector will be
+ used if no ScalarData object is specified.
ramp_0: :class:`~.Ramp`
The starting ramp for the reconstruction. A zero vector will be
used if no Ramp object is specified.
diff --git a/pyramid/regularisator.py b/pyramid/regularisator.py
index 4c4caf2..695860d 100644
--- a/pyramid/regularisator.py
+++ b/pyramid/regularisator.py
@@ -14,8 +14,8 @@ from scipy import sparse
import jutil.diff as jdiff
import jutil.norms as jnorm
-__all__ = ['Regularisator', 'NoneRegularisator', 'ZeroOrderRegularisator',
- 'FirstOrderRegularisator', 'ComboRegularisator']
+__all__ = ['Regularisator', 'NoneRegularisator', 'ZeroOrderRegularisator', 'FirstOrderRegularisator',
+ 'ComboRegularisator']
class Regularisator(object):
@@ -363,11 +363,11 @@ class FirstOrderRegularisator(Regularisator):
"""
- def __init__(self, mask, lam=1E-4, p=2, add_params=0):
+ def __init__(self, mask, lam=1E-4, p=2, add_params=0, factor=3):
self.p = p
- D0 = jdiff.get_diff_operator(mask, 0, 3)
- D1 = jdiff.get_diff_operator(mask, 1, 3)
- D2 = jdiff.get_diff_operator(mask, 2, 3)
+ D0 = jdiff.get_diff_operator(mask, 0, factor)
+ D1 = jdiff.get_diff_operator(mask, 1, factor)
+ D2 = jdiff.get_diff_operator(mask, 2, factor)
D = sparse.vstack([D0, D1, D2])
if p == 2:
norm = jnorm.WeightedL2Square(D)
diff --git a/pyramid/utils/__init__.py b/pyramid/utils/__init__.py
index eaac581..ef1bf38 100644
--- a/pyramid/utils/__init__.py
+++ b/pyramid/utils/__init__.py
@@ -5,11 +5,11 @@
"""Subpackage containing Pyramid utility functions."""
from .pm import pm
-from .reconstruction_2d_from_phasemap import reconstruction_2d_from_phasemap
-from .reconstruction_3d_from_magdata import reconstruction_3d_from_magdata
-from . import lorentz
+from .reconstruction_2d_from_phasemap import reconstruction_2d_from_phasemap, reconstruction_2d_charge_from_phasemap
+from .reconstruction_3d_from_magdata import reconstruction_3d_from_magdata, reconstruction_3d_from_elecdata
#from .phasemap_creator import gui_phasemap_creator
#from .mag_slicer import gui_mag_slicer
-__all__ = ['pm', 'reconstruction_2d_from_phasemap', 'reconstruction_3d_from_magdata', 'lorentz']#,
+__all__ = ['pm', 'reconstruction_2d_from_phasemap', 'reconstruction_3d_from_magdata',
+ 'reconstruction_3d_from_elecdata', 'reconstruction_2d_charge_from_phasemap']#,
#'gui_phasemap_creator', 'gui_mag_slicer']
diff --git a/pyramid/utils/pm.py b/pyramid/utils/pm.py
index fb29652..8d92b55 100644
--- a/pyramid/utils/pm.py
+++ b/pyramid/utils/pm.py
@@ -16,14 +16,13 @@ _log = logging.getLogger(__name__)
# TODO: rename magdata to vecdata everywhere!
-def pm(fielddata, mode='z', b_0=1, electrode_vec=(1E6, 1E6), mapper='RDFC', **kwargs):
- """Convenience function for fast magnetic and electric phase mapping.
+def pm(fielddata, mode='z', b_0=1, electrode_vec=(1E6, 1E6), prw_vec=None, mapper='RDFC', **kwargs):
+ """Convenience function for fast electric charge and magnetic phase mapping.
Parameters
----------
fielddata : :class:`~.VectorData`, or `~.ScalarData`
- A :class:`~.VectorData` or `~.ScalarData` object, from which the projected phase map should
- be calculated.
+ A :class:`~.VectorData` or `~.ScalarData` object, from which the projected phase map should be calculated.
mode: {'z', 'y', 'x', 'x-tilt', 'y-tilt', 'rot-tilt'}, optional
Projection mode which determines the :class:`~.pyramid.projector.Projector` subclass, which
is used for the projection. Default is a simple projection along the `z`-direction.
@@ -32,6 +31,9 @@ def pm(fielddata, mode='z', b_0=1, electrode_vec=(1E6, 1E6), mapper='RDFC', **kw
electrode_vec : tuple of float (N=2)
The norm vector of the counter electrode, (elec_a,elec_b), and the distance to the origin is
the norm of (elec_a,elec_b). The default value is (1E6, 1E6).
+ prw_vec: tuple of 2 int, optional
+ A two-component vector describing the displacement of the reference wave to include
+ perturbation of this reference by the object itself (via fringing fields), (y, x).
mapper : :class: '~. PhaseMap'
A :class: '~. PhaseMap' object, which maps a fielddata into a phase map. The default is 'RDFC'.
**kwargs : additional arguments
@@ -67,11 +69,12 @@ def pm(fielddata, mode='z', b_0=1, electrode_vec=(1E6, 1E6), mapper='RDFC', **kw
phasemapper = PhaseMapperFDFC(fielddata.a, projector.dim_uv, b_0=b_0, padding=padding)
# Set up phasemapper and map phase:
elif mapper == 'Charge':
- phasemapper = PhaseMapperCharge(KernelCharge(fielddata.a, projector.dim_uv, electrode_vec=electrode_vec))
+ phasemapper = PhaseMapperCharge(KernelCharge(fielddata.a, projector.dim_uv,
+ electrode_vec=electrode_vec, prw_vec=prw_vec))
else:
raise ValueError("Invalid mapper (use 'RDFC', 'FDFC' or 'Charge'")
phasemap = phasemapper(field_proj)
# Get mask from fielddata:
- phasemap.mask = field_proj.get_mask()[0, ...]
+ # phasemap.mask = field_proj.get_mask()[0, ...]
# Return phase:
return phasemap
diff --git a/pyramid/utils/reconstruction_2d_from_phasemap.py b/pyramid/utils/reconstruction_2d_from_phasemap.py
index e519e36..6c46183 100644
--- a/pyramid/utils/reconstruction_2d_from_phasemap.py
+++ b/pyramid/utils/reconstruction_2d_from_phasemap.py
@@ -9,12 +9,12 @@ import logging
import numpy as np
from .. import reconstruction
-from ..dataset import DataSet
+from ..dataset import DataSet, DataSetCharge
from ..projector import SimpleProjector
-from ..regularisator import FirstOrderRegularisator, NoneRegularisator
+from ..regularisator import FirstOrderRegularisator, NoneRegularisator, ZeroOrderRegularisator
from ..kernel import KernelCharge
from ..phasemapper import PhaseMapperCharge
-from ..forwardmodel import ForwardModel
+from ..forwardmodel import ForwardModel, ForwardModelCharge
from ..costfunction import Costfunction
from .pm import pm
@@ -22,7 +22,9 @@ __all__ = ['reconstruction_2d_from_phasemap']
_log = logging.getLogger(__name__)
# TODO: lam should NOT have a default!!!
-def reconstruction_2d_from_phasemap(phasemap, b_0=1, lam=1E-3, max_iter=100, ramp_order=1,
+
+
+def reconstruction_2d_from_phasemap(phasemap, b_0=1, lam=1E-3, max_iter=100, ramp_order=None,
plot_results=False, ar_dens=None, verbose=True):
"""Convenience function for reconstructing a projected distribution from a single phasemap.
@@ -109,9 +111,8 @@ def reconstruction_2d_from_phasemap(phasemap, b_0=1, lam=1E-3, max_iter=100, ram
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,
+def reconstruction_2d_charge_from_phasemap(phasemap, max_iter=1000, ramp_order=None, mask=None,
+ lam=None, electrode_vec=(1E6, 1E6), v_acc=300000, prw=None,
plot_results=False, verbose=True):
"""Convenience function for reconstructing a projected distribution from a single phasemap.
@@ -119,45 +120,53 @@ def reconstruction_2d_charge_from_phasemap(phasemap, lam=1E-3, max_iter=100, ram
----------
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.
+ The maximum number of iterations for the optimization.
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).
+ lam: float,
+ The zero order regularisator parameter. 'None' means no regularisator.
+ mask: ndarrary,
+ Define where situate the reconstructed charges
+ electrode_vec : tuple of float (N=2)
+ The norm vector of the counter electrode in pixels, (elec_a,elec_b), and the distance to the origin is
+ the norm of (elec_a,elec_b).
+ v_acc: float
+ The accelerating voltage of electrons.
+ prw: tuple of 2 int, optional
+ A two-component vector describing the displacement of the reference wave to include
+ perturbation of this reference by the object itself (via fringing fields), (y, x).
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`
+ elecdata_rec, cost: :class:`~.ScalarData`, :class:`~.Costfunction`
The reconstructed magnetisation distribution and the used costfunction.
"""
- _log.debug('Calling reconstruction_2d_from_phasemap')
+ _log.debug('Calling reconstruction_2d_charge_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 = DataSetCharge(phasemap.a, dim, electrode_vec, mask=mask)
+ kernel = KernelCharge(phasemap.a, phasemap.dim_uv, electrode_vec=electrode_vec, v_acc=v_acc, prw_vec=prw)
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)
+ fwd_model = ForwardModelCharge(data, ramp_order)
+ if lam is None:
+ reg = NoneRegularisator() # FirstOrderRegularisator(data.mask, lam, add_params=fwd_model.ramp.n)
+ else:
+ # reg = FirstOrderRegularisator(data.mask, lam=lam, p=2, add_params=fwd_model.ramp.n, factor=1)
+ reg = ZeroOrderRegularisator(data.mask, lam=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)
+ elecdata_rec = reconstruction.optimize_linear_charge(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)
@@ -169,11 +178,7 @@ def reconstruction_2d_charge_from_phasemap(phasemap, lam=1E-3, max_iter=100, ram
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)
+ phasemap_rec = pm(elecdata_rec, mapper='Charge')
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()
@@ -196,5 +201,5 @@ def reconstruction_2d_charge_from_phasemap(phasemap, lam=1E-3, max_iter=100, ram
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
+ # Return reconstructed charge distribution and cost function:
+ return elecdata_rec, cost
diff --git a/pyramid/utils/reconstruction_3d_from_magdata.py b/pyramid/utils/reconstruction_3d_from_magdata.py
index 61a3161..4020a82 100644
--- a/pyramid/utils/reconstruction_3d_from_magdata.py
+++ b/pyramid/utils/reconstruction_3d_from_magdata.py
@@ -13,16 +13,16 @@ import multiprocessing as mp
from jutil.taketime import TakeTime
from .. import reconstruction
-from ..dataset import DataSet
+from ..dataset import DataSet, DataSetCharge
from ..projector import XTiltProjector, YTiltProjector
from ..ramp import Ramp
-from ..regularisator import FirstOrderRegularisator
-from ..forwardmodel import ForwardModel, DistributedForwardModel
+from ..regularisator import FirstOrderRegularisator, NoneRegularisator, ZeroOrderRegularisator
+from ..forwardmodel import ForwardModel, DistributedForwardModel, ForwardModelCharge
from ..costfunction import Costfunction
-from ..phasemapper import PhaseMapperRDFC
-from ..kernel import Kernel
+from ..phasemapper import PhaseMapperRDFC, PhaseMapperCharge
+from ..kernel import Kernel, KernelCharge
-__all__ = ['reconstruction_3d_from_magdata']
+__all__ = ['reconstruction_3d_from_magdata', 'reconstruction_3d_from_elecdata']
_log = logging.getLogger(__name__)
@@ -152,3 +152,142 @@ def reconstruction_3d_from_magdata(magdata, b_0=1, lam=1E-3, max_iter=100, ramp_
ar_dens=ar_dens, coloring='amplitude')
# Return reconstructed magnetisation distribution and cost function:
return magdata_rec, cost
+
+
+def reconstruction_3d_from_elecdata(elecdata, electrode_vec=(1E6, 1E6), max_iter=100, ramp_order=1,
+ angles=np.linspace(-90, 90, num=19), dim_uv=None, lam=None,
+ axes=(True, False), noise=0, offset_max=0, ramp_max=0, confidence=None,
+ use_internal_mask=False, mask_3d_threshold=0.9, mask=None, plot_results=False,
+ plot_input=False, multicore=False, verbose=True):
+ """Convenience function for reconstructing a projected distribution from a single phasemap.
+
+ Parameters
+ ----------1
+ elecdata: :class:`~.ScalarData`
+ The Charge distribution which should be used for the reconstruction.
+ The default is 1.
+ electrode_vec : tuple of float (N=2)
+ The norm vector of the counter electrode in pixels, (elec_a,elec_b), and the distance to the origin is
+ the norm of (elec_a,elec_b).
+ max_iter : int, optional
+ The maximum number of iterations for the optimization.
+ 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).
+ angles: :class:`~numpy.ndarray` (N=1), optional
+ Numpy array determining the angles which should be used for the projectors in x- and
+ y-direction. This implicitly sets the number of images per rotation axis. Defaults to a
+ range from -90° to 90° degrees, in 10° steps.
+ dim_uv: int or None (default)
+ Determines if the phasemaps should be padded to a certain size while calculating.
+ lam: float,
+ The zero order regularisator parameter. 'None' means no regularisator.
+ axes: tuple of booleans (N=2), optional
+ Determines if both tilt axes should be calculated. The order is (x, y), both are True by
+ default.
+ noise: float, optional
+ If this is not zero, random gaussian noise with this as a maximum value will be applied
+ to all calculated phasemaps. The default is 0. The unit is radians.
+ offset_max: float, optional
+ if this is not zero, a random offset with this as a maximum value will be applied to all
+ calculated phasemaps. The default is 0.
+ ramp_max: float, optional
+ if this is not zero, a random linear ramp with this as a maximum value will be applied
+ to both axes of all calculated phasemaps. The default is 0.
+ confidencee: boolean, optional
+ if not None, define the trust of the phase image. Here we use the phasemap.mask.
+ use_internal_mask: boolean, optional
+ If '3D', the mask from the input charge distribution is taken for the
+ reconstruction. If '2D', the mask is calculated via logic backprojection from the 2D-masks
+ of the input phasemaps. If 'manual', the mask is given by user.
+ mask_3d_threshold: float, optional，
+ Provide if 'use_internal_mask' is '2D', which determines the mask from phasemaps.
+ mask: :class:`~numpy.ndarray` (N=3), provide if `use_internal_mask` set to manual
+ A boolean mask which defines the magnetized volume in 3D.
+ plot_results: boolean, optional
+ If True, the results are plotted after reconstruction.
+ plot_input:
+ If True, the input phasemaps are plotted after reconstruction.
+ multicore: boolean, optional
+ Determines if multiprocessing should be used. Default is True. Phasemap calculations
+ will be divided onto the separate cores.
+ verbose: bool, optional
+ If set to True, information like a progressbar is displayed during reconstruction.
+ The default is False.
+
+ Returns
+ -------
+ elecdata_rec, cost: :class:`~.ScalarData`, :class:`~.Costfunction`
+ The reconstructed charge distribution and the used costfunction.
+
+ """
+ _log.debug('Calling reconstruction_3d_from_elecdata')
+ # Construct DataSet:
+ data = DataSetCharge(elecdata.a, elecdata.dim, electrode_vec=electrode_vec)
+ # Construct projectors:
+ projectors = []
+ # Construct data set and regularisator:
+ for angle in angles:
+ angle_rad = angle * np.pi / 180
+ if axes[0]:
+ projectors.append(XTiltProjector(elecdata.dim, angle_rad, dim_uv))
+ if axes[1]:
+ projectors.append(YTiltProjector(elecdata.dim, angle_rad, dim_uv))
+ # Add pairs of projectors and according phasemaps to the DataSet:
+ for projector in projectors:
+ elec_proj = projector(elecdata)
+ phasemap = PhaseMapperCharge(KernelCharge(elecdata.a, projector.dim_uv, electrode_vec))(elec_proj)
+ phasemap.mask = elec_proj.get_mask()[0, ...]
+ if confidence is not None:
+ phasemap.confidence = np.logical_not(elec_proj.get_mask()[0, ...])
+ data.append(phasemap, projector)
+ # Add offset and ramp if necessary:
+ for i, phasemap in enumerate(data.phasemaps):
+ offset = np.random.uniform(-offset_max, offset_max)
+ ramp_u = np.random.uniform(-ramp_max, ramp_max)
+ ramp_v = np.random.uniform(-ramp_max, ramp_max)
+ 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: # 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
+ # Construct mask:
+ if use_internal_mask == '3D':
+ data.mask = elecdata.get_mask() # Use perfect mask from elecdata!
+ elif use_internal_mask == '2D':
+ data.set_3d_mask(threshold=mask_3d_threshold) # Construct mask from 2D phase masks!
+ elif use_internal_mask == 'manual':
+ data.mask = mask
+ elif not use_internal_mask:
+ data.mask = None
+ else:
+ raise ValueError('use_internal_mask value is wrong')
+ # Construct regularisator, forward model and costfunction:
+ if multicore:
+ mp.freeze_support()
+ fwd_model = DistributedForwardModel(data, ramp_order=ramp_order, nprocs=mp.cpu_count())
+ else:
+ fwd_model = ForwardModelCharge(data, ramp_order=ramp_order)
+ if lam is None:
+ reg = NoneRegularisator() # FirstOrderRegularisator(data.mask, lam, add_params=fwd_model.ramp.n)
+ else:
+ reg = ZeroOrderRegularisator(lam=lam, add_params=fwd_model.ramp.n)
+ # reg = NoneRegularisator()
+ cost = Costfunction(fwd_model, reg)
+ # Reconstruct and save:
+ elecdata_rec = reconstruction.optimize_linear_charge(cost, max_iter=max_iter, verbose=verbose)
+ # Finalize ForwardModel (returns workers if multicore):
+ fwd_model.finalize()
+ # Plot input:
+ if plot_input:
+ data.plot_phasemaps(symmetric=False)
+ # Plot results:
+ if plot_results:
+ data.plot_mask()
+ elecdata_rec.plot_field(cmap='viridis')
+ # Return reconstructed charge distribution and cost function:
+ return elecdata_rec, cost
+
--
GitLab