Kernel corrections and ForwardModel properties enhancement (more modular now).

kernel: u and v are switched and now finally correct.
        get_elementary_phase is also switched and corrected.
phasemapper: added PhaseMapperCharge from Feng Shan.
forward_model: now uses y instead of phase_vec/Se_inv and now direct properties.
costfunction: corrected some docstrings and documentation comments.
fielddata: corrected some missing integer divisions (// instead of /),
        added plot_streamline().

pyramid/colors.py

@@ -5,10 +5,7 @@
 """This module provides a custom :class:`~.HLSTriadicColormap` colormap class which has a few
 additional functions and can encode three-dimensional directions."""
 
-# TODO: ALL docstrings! Also plotting docstrings in PhaseMap and VectorData!
-
 import logging
-from numbers import Number
 
 import matplotlib.pyplot as plt
 import numpy as np
@@ -103,7 +100,7 @@ class HLSTetradicColormap(colors.LinearSegmentedColormap):
              'green': [(0.00, 0.0, 0.0),
                        (0.25, 0.0, 0.0),
                        (0.50, 1.0, 1.0),
-                       (0.75, 1.0, 01.0),
+                       (0.75, 1.0, 1.0),
                        (1.00, 0.0, 0.0)],
 
              'blue': [(0.00, 0.0, 0.0),

pyramid/costfunction.py

@@ -27,38 +27,38 @@ class Costfunction(object):
 
     Attributes
     ----------
-    regularisator : :class:`~.Regularisator`
-        Regularisator class that's responsible for the regularisation term.
-    y : :class:`~numpy.ndarray` (N=1)
-        Vector which lists all pixel values of all phase maps one after another.
     fwd_model : :class:`~.ForwardModel`
         The Forward model instance which should be used for the simulation of the phase maps which
         will be compared to `y`.
-    Se_inv : :class:`~numpy.ndarray` (N=2), optional
-        Inverted covariance matrix of the measurement errors. The matrix has size `NxN` with N
-        being the length of the targetvector y (vectorized phase map information).
+    regularisator : :class:`~.Regularisator`, optional
+        Regularisator class that's responsible for the regularisation term. If `None` or none is
+        given, no regularisation will be used.
+    y : :class:`~numpy.ndarray` (N=1)
+        Vector which lists all pixel values of all phase maps one after another.
     m: int
         Size of the image space.
     n: int
         Size of the input space.
+    Se_inv : :class:`~numpy.ndarray` (N=2), optional
+        Inverted covariance matrix of the measurement errors. The matrix has size `m x m` with m
+        being the length of the targetvector y.
 
     """
 
     _log = logging.getLogger(__name__ + '.Costfunction')
 
-    def __init__(self, fwd_model, regularisator):
+    def __init__(self, fwd_model, regularisator=None):
         self._log.debug('Calling __init__')
         self.fwd_model = fwd_model
-        self.regularisator = regularisator
-        if self.regularisator is None:
+        if regularisator is None:
             self.regularisator = NoneRegularisator()
-        # Extract important information:
-        self.y = self.fwd_model.data_set.phase_vec
+        else:
+            self.regularisator = regularisator
+        # Extract information from fwd_model:
+        self.y = self.fwd_model.y
         self.n = self.fwd_model.n
         self.m = self.fwd_model.m
-        if self.fwd_model.data_set.Se_inv is None:
-            self.fwd_model.data_set.set_Se_inv_diag_with_conf()
-        self.Se_inv = self.fwd_model.data_set.Se_inv
+        self.Se_inv = self.fwd_model.Se_inv
         self._log.debug('Created ' + str(self))
 
     def __repr__(self):

pyramid/dataset.py

@@ -42,10 +42,17 @@ class DataSet(object):
     Se_inv : :class:`~numpy.ndarray` (N=2), optional
         Inverted covariance matrix of the measurement errors. The matrix has size `NxN` with N
         being the length of the targetvector y (vectorized phase map information).
-    phasemaps:
-        A list of all stored :class:`~.PhaseMap` objects.
     projectors: list of :class:`~.Projector`
         A list of all stored :class:`~.Projector` objects.
+    phasemaps: list of :class:`~.PhaseMap`
+        A list of all stored :class:`~.PhaseMap` objects.
+    phase_vec: :class:`~numpy.ndarray` (N=1)
+        The concatenaded, vectorized phase of all :class:`~.PhaseMap` objects.
+    count(self): int
+        Number of phase maps and projectors in the dataset.
+    hook_points(self): :class:`~numpy.ndarray` (N=1)
+        Hook points which determine the start of values of a phase map in the `phase_vec`.
+        The length is `count + 1`.
 
     """
 

pyramid/fielddata.py

@@ -437,7 +437,7 @@ class VectorData(FieldData):
                 self.field = np.pad(self.field, ((0, 0), (0, pz), (0, py), (0, px)),
                                     mode='constant')
             # Create coarser grid for the vector field:
-            shape_4d = (3, self.dim[0] / 2, 2, self.dim[1] / 2, 2, self.dim[2] / 2, 2)
+            shape_4d = (3, self.dim[0] // 2, 2, self.dim[1] // 2, 2, self.dim[2] // 2, 2)
             self.field = self.field.reshape(shape_4d).mean(axis=(6, 4, 2))
 
     def scale_up(self, n=1, order=0):
@@ -820,6 +820,137 @@ class VectorData(FieldData):
         # Return plotting axis:
         return axis
 
+    def plot_streamline(self, title='Vector Field', axis=None, proj_axis='z', figsize=(8.5, 7),
+                        coloring='angle', ax_slice=None, density=2, linewidth=2,
+                        show_mask=True, bgcolor='white', hue_mode='triadic'):
+        """Plot a slice of the vector field as a quiver plot.
+
+        Parameters
+        ----------
+        title : string, optional
+            The title for the plot.
+        axis : :class:`~matplotlib.axes.AxesSubplot`, optional
+            Axis on which the graph is plotted. Creates a new figure if none is specified.
+        proj_axis : {'z', 'y', 'x'}, optional
+            The axis, from which a slice is plotted. The default is 'z'.
+        figsize : tuple of floats (N=2)
+            Size of the plot figure.
+        coloring : {'angle', 'amplitude', 'uniform'}
+            Color coding mode of the arrows. Use 'full' (default), 'angle', 'amplitude' or
+            'uniform'.
+        ax_slice : int, optional
+            The slice-index of the axis specified in `proj_axis`. Is set to the center of
+            `proj_axis` if not specified.
+        density : float or 2-tuple, optional
+            Controls the closeness of streamlines. When density = 1, the domain is divided into a
+            30x30 grid—density linearly scales this grid. Each cebll in the grid can have, at most,
+            one traversing streamline. For different densities in each direction, use
+            [density_x, density_y].
+        linewidth : numeric or 2d array, optional
+            Vary linewidth when given a 2d array with the same shape as velocities.
+        show_mask: boolean
+            Default is True. Shows the outlines of the mask slice if available.
+        bgcolor: {'black', 'white'}, optional
+            Determines the background color of the plot.
+        hue_mode : {'triadic', 'tetradic'}
+            Optional string for determining the hue scheme. Use either a triadic or tetradic
+            scheme (see the according colormaps for more information).
+
+        Returns
+        -------
+        axis: :class:`~matplotlib.axes.AxesSubplot`
+            The axis on which the graph is plotted.
+
+        """
+        self._log.debug('Calling plot_quiver')
+        assert proj_axis == 'z' or proj_axis == 'y' or proj_axis == 'x', \
+            'Axis has to be x, y or z (as string).'
+        if ax_slice is None:
+            ax_slice = self.dim[{'z': 0, 'y': 1, 'x': 2}[proj_axis]] // 2
+        u_mag, v_mag = self.get_slice(ax_slice, proj_axis)
+        if proj_axis == 'z':  # Slice of the xy-plane with z = ax_slice
+            u_label = 'x [px]'
+            v_label = 'y [px]'
+            submask = self.get_mask()[ax_slice, ...]
+        elif proj_axis == 'y':  # Slice of the xz-plane with y = ax_slice
+            u_label = 'x [px]'
+            v_label = 'z [px]'
+            submask = self.get_mask()[:, ax_slice, :]
+        elif proj_axis == 'x':  # Slice of the yz-plane with x = ax_slice
+            u_label = 'z [px]'
+            v_label = 'y [px]'
+            submask = self.get_mask()[..., ax_slice]
+        else:
+            raise ValueError('{} is not a valid argument (use x, y or z)'.format(proj_axis))
+        # Prepare quiver (select only used arrows if ar_dens is specified):
+        dim_uv = u_mag.shape
+        uu = np.arange(dim_uv[1]) + 0.5  # shift to center of pixel
+        vv = np.arange(dim_uv[0]) + 0.5  # shift to center of pixel
+        u_mag, v_mag = self.get_slice(ax_slice, proj_axis)
+        # v_mag = np.ma.array(v_mag, mask=submask)
+        amplitudes = np.hypot(u_mag, v_mag)
+        # Calculate the arrow colors:
+        if coloring == 'angle':
+            self._log.debug('Encoding angles')
+            color = np.arctan2(v_mag, u_mag) / (2 * np.pi)
+            color[color < 0] += 1
+            if hue_mode == 'triadic':
+                cmap = colors.hls_triadic_cmap
+            elif hue_mode == 'tetradic':
+                cmap = colors.hls_tetradic_cmap
+            else:
+                raise ValueError('Hue mode {} not understood!'.format(hue_mode))
+        elif coloring == 'amplitude':
+            self._log.debug('Encoding amplitude')
+            color = amplitudes / amplitudes.max()
+            cmap = 'jet'
+        elif coloring == 'uniform':
+            self._log.debug('No color encoding')
+            color = np.zeros_like(u_mag)  # use black arrows!
+            cmap = 'gray' if bgcolor == 'white' else 'Greys'
+        else:
+            raise AttributeError("Invalid coloring mode! Use 'angles', 'amplitude' or 'uniform'!")
+        # If no axis is specified, a new figure is created:
+        if axis is None:
+            self._log.debug('axis is None')
+            fig = plt.figure(figsize=figsize)
+            axis = fig.add_subplot(1, 1, 1)
+        axis.set_aspect('equal')
+        # Plot the streamlines:
+        im = plt.streamplot(uu, vv, u_mag, v_mag, density=density, linewidth=linewidth,
+                            color=color, cmap=cmap)
+        if coloring == 'amplitude':
+            fig = plt.gcf()
+            fig.subplots_adjust(right=0.8)
+            cbar_ax = fig.add_axes([0.82, 0.15, 0.02, 0.7])
+            cbar = fig.colorbar(im.lines, cax=cbar_ax)
+            cbar.ax.tick_params(labelsize=14)
+            cbar_title = u'amplitude'
+            cbar.set_label(cbar_title, fontsize=15)
+        # Change background color:
+        axis.set_axis_bgcolor(bgcolor)
+        # 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
+            mask_color = 'white' if bgcolor == 'black' else 'black'
+            axis.contour(uu, vv, submask, levels=[0.5], colors=mask_color,
+                         linestyles='dotted', linewidths=2)
+        # Further plot formatting:
+        axis.set_xlim(0, dim_uv[1])
+        axis.set_ylim(0, dim_uv[0])
+        axis.set_title(title, fontsize=18)
+        axis.set_xlabel(u_label, fontsize=15)
+        axis.set_ylabel(v_label, fontsize=15)
+        axis.tick_params(axis='both', which='major', labelsize=14)
+        if dim_uv[0] >= dim_uv[1]:
+            u_bin, v_bin = np.max((2, np.floor(9 * dim_uv[1] / dim_uv[0]))), 9
+        else:
+            u_bin, v_bin = 9, np.max((2, np.floor(9 * dim_uv[0] / dim_uv[1])))
+        axis.xaxis.set_major_locator(MaxNLocator(nbins=u_bin, integer=True))
+        axis.yaxis.set_major_locator(MaxNLocator(nbins=v_bin, integer=True))
+        # Return plotting axis:
+        return axis
+
     def plot_quiver(self, title='Vector Field', axis=None, proj_axis='z', figsize=(8.5, 7),
                     coloring='angle', ar_dens=1, ax_slice=None, log=False, scaled=True,
                     scale=1., show_mask=True, bgcolor='white', hue_mode='triadic'):

pyramid/forwardmodel.py

@@ -36,10 +36,15 @@ class ForwardModel(object):
         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).
+    y : :class:`~numpy.ndarray` (N=1)
+        Vector which lists all pixel values of all phase maps one after another.
     m: int
         Size of the image space. Number of pixels of the 2-dimensional projected grid.
     n: int
         Size of the input space. Number of voxels of the 3-dimensional grid.
+    Se_inv : :class:`~numpy.ndarray` (N=2), optional
+        Inverted covariance matrix of the measurement errors. The matrix has size `m x m` with m
+        being the length of the targetvector y (vectorized phase map information).
 
     """
 
@@ -49,14 +54,20 @@ class ForwardModel(object):
         self._log.debug('Calling __init__')
         self.data_set = data_set
         self.ramp_order = ramp_order
+        # Extract information from data_set:
         self.phasemappers = self.data_set.phasemappers
-        self.ramp = Ramp(self.data_set, self.ramp_order)
-        self.m = self.data_set.m
+        self.y = self.data_set.phase_vec
         self.n = self.data_set.n
-        if self.ramp.n is not None:  # Additional parameters have to be fitted!
-            self.n += self.ramp.n
+        self.m = self.data_set.m
         self.shape = (self.m, self.n)
-        self.hook_points = data_set.hook_points
+        self.hook_points = self.data_set.hook_points
+        if self.data_set.Se_inv is None:
+            self.data_set.set_Se_inv_diag_with_conf()
+        self.Se_inv = self.data_set.Se_inv
+        # Create ramp and change n accordingly:
+        self.ramp = Ramp(self.data_set, self.ramp_order)
+        self.n += self.ramp.n  # ramp.n is 0 if ramp_order is None
+        # Create empty MagData object:
         self.magdata = VectorData(self.data_set.a, np.zeros((3,) + self.data_set.dim))
         self._log.debug('Creating ' + str(self))
 

pyramid/kernel.py

@@ -87,9 +87,9 @@ class Kernel(object):
         self.slice_mag = (slice(0, dim_uv[0]),  # Magnetization is padded on the far end!
                           slice(0, dim_uv[1]))  # (Phase cutout is shifted as listed above)
         # Calculate kernel (single pixel phase):
-        # [M_0] = [PHI_0 / µ_0] = Tm² / N/A² = N/Am * A²/N = A/m
-        #       --> This is the magnetization, not the magnetic moment (A/m * m³ = Am²)!
-        # [b_0] = [µ_0] * [M_0] = A/m * N/A² = N/Am = T
+        # [M_0] = A/m  --> This is the magnetization, not the magnetic moment (A/m * m³ = Am²)!
+        # [PHI_0 / µ_0] = Tm² / Tm/A = Am
+        # [b_0] = [M_0] * [µ_0] = A/m * N/A² = N/Am = T
         # [coeff] = [b_0 * a² / (2*PHI_0)] = T * m² / Tm² = 1  --> without unit (phase)!
         coeff = b_0 * a ** 2 / (2 * PHI_0)  # Minus is gone because of negative z-direction
         v_dim, u_dim = dim_uv
@@ -99,13 +99,13 @@ class Kernel(object):
         self.u = fft.empty(self.dim_kern, dtype=fft.FLOAT)
         self.v = fft.empty(self.dim_kern, dtype=fft.FLOAT)
         self.u[...] = coeff * self._get_elementary_phase(geometry, uu, vv, a)
-        self.v[...] = coeff * self._get_elementary_phase(geometry, vv, uu, a)
+        self.v[...] = coeff * -self._get_elementary_phase(geometry, vv, uu, a)
         # Include perturbed reference wave:
         if prw_vec is not None:
             uu += prw_vec[1]
             vv += prw_vec[0]
             self.u[...] -= coeff * self._get_elementary_phase(geometry, uu, vv, a)
-            self.v[...] -= coeff * self._get_elementary_phase(geometry, vv, uu, a)
+            self.v[...] -= coeff * -self._get_elementary_phase(geometry, vv, uu, a)
         # Calculate Fourier trafo of kernel components:
         self.u_fft = fft.rfftn(self.u, self.dim_pad)
         self.v_fft = fft.rfftn(self.v, self.dim_pad)
@@ -125,7 +125,7 @@ class Kernel(object):
         self._log.debug('Calling _get_elementary_phase')
         if geometry == 'disc':
             in_or_out = ~ np.logical_and(n == 0, m == 0)
-            return n / (n ** 2 + m ** 2 + 1E-30) * in_or_out
+            return m / (n ** 2 + m ** 2 + 1E-30) * in_or_out
         elif geometry == 'slab':
             def _F_a(n, m):
                 A = np.log(a ** 2 * (n ** 2 + m ** 2))

pyramid/phasemapper.py

@@ -11,13 +11,12 @@ import abc
 import logging
 
 import numpy as np
-from numpy import pi
 
 from . import fft
 from .fielddata import VectorData, ScalarData
 from .phasemap import PhaseMap
 
-__all__ = ['PhaseMapperRDFC', 'PhaseMapperFDFC', 'PhaseMapperMIP']
+__all__ = ['PhaseMapperRDFC', 'PhaseMapperFDFC', 'PhaseMapperMIP', 'PhaseMapperCharge']
 _log = logging.getLogger(__name__)
 
 PHI_0 = 2067.83  # magnetic flux in T*nm²
@@ -25,6 +24,7 @@ H_BAR = 6.626E-34  # Planck constant in J*s
 M_E = 9.109E-31  # electron mass in kg
 Q_E = 1.602E-19  # electron charge in C
 C = 2.998E8  # speed of light in m/s
+EPS_0 = 8.8542E-12  # electrical field constant
 
 
 class PhaseMapper(object, metaclass=abc.ABCMeta):
@@ -135,7 +135,7 @@ class PhaseMapperRDFC(PhaseMapper):
         self.u_mag_fft = fft.rfftn(self.u_mag)
         self.v_mag_fft = fft.rfftn(self.v_mag)
         # Convolve the magnetization with the kernel in Fourier space:
-        self.phase_fft = self.u_mag_fft * self.kernel.v_fft - self.v_mag_fft * self.kernel.u_fft
+        self.phase_fft = self.u_mag_fft * self.kernel.u_fft + self.v_mag_fft * self.kernel.v_fft
         # Return the result:
         return fft.irfftn(self.phase_fft)[self.kernel.slice_phase]
 
@@ -181,8 +181,8 @@ class PhaseMapperRDFC(PhaseMapper):
             'vector size not compatible! vector: {}, size: {}'.format(len(vector), self.m)
         self.mag_adj[self.kernel.slice_mag] = vector.reshape(self.kernel.dim_uv)
         mag_adj_fft = fft.irfftn_adj(self.mag_adj)
-        u_phase_adj_fft = mag_adj_fft * self.kernel.v_fft
-        v_phase_adj_fft = mag_adj_fft * -self.kernel.u_fft
+        u_phase_adj_fft = mag_adj_fft * self.kernel.u_fft
+        v_phase_adj_fft = mag_adj_fft * self.kernel.v_fft
         u_phase_adj = fft.rfftn_adj(u_phase_adj_fft)[self.kernel.slice_phase]
         v_phase_adj = fft.rfftn_adj(v_phase_adj_fft)[self.kernel.slice_phase]
         result = np.concatenate((-u_phase_adj.ravel(), -v_phase_adj.ravel()))
@@ -326,7 +326,7 @@ class PhaseMapperMIP(PhaseMapper):
     v_0 : float, optional
         The mean inner potential of the specimen in V. The default is 1.
     v_acc : float, optional
-        The acceleration voltage of the electron microscope in V. The default is 30000.
+        The acceleration voltage of the electron microscope in V. The default is 300000.
     threshold : float, optional
         Threshold for the recognition of the specimen location. The default is 0, which means that
         all voxels with non-zero magnetization will contribute. Should be above noise level.
@@ -350,7 +350,7 @@ class PhaseMapperMIP(PhaseMapper):
         self.n = self.m
         # Coefficient calculation:
         lam = H_BAR / np.sqrt(2 * M_E * Q_E * v_acc * (1 + Q_E * v_acc / (2 * M_E * C ** 2)))
-        C_e = 2 * pi * Q_E / lam * (Q_E * v_acc + M_E * C ** 2) / (
+        C_e = 2 * np.pi * Q_E / lam * (Q_E * v_acc + M_E * C ** 2) / (
             Q_E * v_acc * (Q_E * v_acc + 2 * M_E * C ** 2))
         self.coeff = v_0 * C_e * a * 1E-9
         self._log.debug('Created ' + str(self))
@@ -407,3 +407,105 @@ class PhaseMapperMIP(PhaseMapper):
 
         """
         raise NotImplementedError()  # TODO: Implement right!
+
+
+class PhaseMapperCharge(PhaseMapper):
+
+    """"""
+
+    def __init__(self, a, dim_uv, biprism_vec, v_acc=300000):
+        self._log.debug('Calling __init__')
+        self.a = a
+        self.dim_uv = dim_uv
+        self.biprism_vec = biprism_vec
+        self.v_acc = v_acc
+        lam = H_BAR / np.sqrt(2 * M_E * Q_E * v_acc * (1 + Q_E * v_acc / (2 * M_E * C ** 2)))
+        C_e = 2 * np.pi * Q_E / lam * (Q_E * v_acc + M_E * C ** 2) / (
+            Q_E * v_acc * (Q_E * v_acc + 2 * M_E * C ** 2))
+        self.coeff = C_e * Q_E / (4 * np.pi * EPS_0)
+
+        self._log.debug('Created ' + str(self))
+
+    def __repr__(self):
+        self._log.debug('Calling __repr__')
+        return '%s(a=%r, dim_uv=%r, v_acc=%r)' % \
+               (self.__class__, self.a, self.dim_uv, self.v_acc)
+
+    def __str__(self):
+        self._log.debug('Calling __str__')
+        return 'PhaseMapperCharge(a=%s, dim_uv=%s, v_acc=%s)' % \
+               (self.a, self.dim_uv, self.v_acc)
+
+    def __call__(self, elec_data):
+        """ phase_dipoles() is to caculate the phase from many electric dipoles.
+        field include the amount of charge in every grid, unit:electron.
+        The normal vector of the electrode is (a,b), and the distance to the origin is the norm
+        of (a,b). R is the sampling rate,pixel/nm."""
+        R = 1 / self.a
+        field = elec_data.field[0, ...]
+        bu, bv = self.biprism_vec  # biprism vector (orthogonal to biprism from origin)
+        bn = bu ** 2 + bv ** 2  # norm of the biprism vector
+        dim_v, dim_u = field.shape
+        # Find list of charge locations and charge values:
+        vq, uq = np.nonzero(field)
+        q = field[vq, uq]
+        # Calculate locations of image charges:
+        vm = (bu ** 2 - bv ** 2) / bn * vq - 2 * bu * bv / bn * uq + 2 * bv
+        um = (bv ** 2 - bu ** 2) / bn * uq - 2 * bu * bv / bn * vq + 2 * bu
+        # Calculate phase contribution for each charge:
+        phase = np.zeros((dim_v, dim_u))
+        for i in range(len(q)):
+            for u in range(dim_u):
+                for v in range(dim_v):
+                    rq = np.sqrt((u - uq[i]) ** 2 + (v - vq[i]) ** 2)  # charge distance
+                    rm = np.sqrt((u - um[i]) ** 2 + (v - vm[i]) ** 2)  # mirror distance
+                    # distance
+                    z1 = (R / 2) ** 2 - rq ** 2
+                    z2 = (R / 2) ** 2 - rm ** 2
+                    if z1 < 0 and z2 < 0:
+                        phase[v, u] += -q[i] * self.coeff * np.log((rq ** 2) / (rm ** 2))
+                    elif z1 >= 0 >= z2:
+                        z3 = np.sqrt(z1)
+                        phase[v, u] += (-q[i] * self.coeff *
+                                        np.log(((z3 + R / 2) ** 2) / (rm ** 2))
+                                        + 2 * q[i] * Q_E * z3 / 2 / np.pi / EPS_0 / R)
+                    else:
+                        z4 = np.sqrt(z2)
+                        phase[v, u] += (-q[i] * self.coeff *
+                                        np.log((rq ** 2) / ((z4 + R / 2) ** 2))
+                                        - 2 * q[i] * Q_E * z4 / 2 / np.pi / EPS_0 / R)
+        return PhaseMap(self.a, phase)
+
+    def jac_dot(self, vector):
+        """Calculate the product of the Jacobi matrix with a given `vector`.
+
+        Parameters
+        ----------
+        vector : :class:`~numpy.ndarray` (N=1)
+            Vectorized form of the electrostatic field of every pixel (row-wise).
+
+        Returns
+        -------
+        result : :class:`~numpy.ndarray` (N=1)
+            Product of the Jacobi matrix (which is not explicitely calculated) with the vector.
+
+        """
+        raise NotImplementedError()  # TODO: Implement right!
+
+    def jac_T_dot(self, vector):
+        """Calculate the product of the transposed Jacobi matrix with a given `vector`.
+
+        Parameters
+        ----------
+        vector : :class:`~numpy.ndarray` (N=1)
+            Vector with ``N**2`` entries which represents a matrix with dimensions like a scalar
+            phasemap.
+
+        Returns
+        -------
+        result : :class:`~numpy.ndarray` (N=1)
+            Product of the transposed Jacobi matrix (which is not explicitely calculated) with
+            the vector, which has ``N**2`` entries like an electrostatic projection.
+
+        """
+        raise NotImplementedError()  # TODO: Implement right!

pyramid/ramp.py

@@ -67,7 +67,7 @@ class Ramp(object):
         self.order = order
         self.deg_of_freedom = (1 + 2 * self.order) if self.order is not None else 0
         self.param_cache = np.zeros((self.deg_of_freedom, self.data_set.count))
-        self.n = self.deg_of_freedom * self.data_set.count
+        self.n = self.deg_of_freedom * self.data_set.count  # 0 if order is None
 
     def __call__(self, index, dof_list=None):
         if self.order is None:  # Do nothing if order is None!