Q3 d

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!):

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))

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!

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):

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.
 
         """

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

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.

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)

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']

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

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

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
+