New: Linear fitting 2.0

doc/user_guide/model.rst

@@ -628,6 +628,29 @@ To fit the model to the data at the current coordinates (e.g. to fit one
 spectrum at a particular point in a spectrum-image) use
 :py:meth:`~.model.BaseModel.fit`.
 
+There are two parent methods of model-fitting available in in HyperSpy.
+**Nonlinear** fitting (default: ``'leastsq'``) is used to fit all parameters in the components,
+minimizing the difference between the fit and the data by successive
+approximations. This difference is typically the squared error term
+(‘ls’), but can also use maximum likelihood estimation (‘ml’) in the
+case of Poisson noise. Nonlinear fitting is slow on large models, but is
+necessary if the components are nonlinear across the data. A component
+is linear when its free parameters scale the component only in the y-axis. For
+the example function ``y = a*x**b``, ``a`` is linear whilst ``b`` 
+is not. Components can also be made up of several linear parts. For instance,
+the 2D-polynomial ``y = a*x**2+b*y**2+c*x+d*y+e`` is entirely linear.
+
+If all components in the model are linear, then **linear** fitting (default: ``'linear'``)can be used. 
+Linear fitting uses linear regression to solve the relation
+``Ax = b`` for x, where b is the data and A are the components. It is
+very fast, but less flexible as it assumes that the nonlinear
+parameters are correct.
+
+A combination of nonlinear and linear fitting can be used to speed up
+the fitting process. For instance, a Gaussian's sigma and centre parameters can 
+be estimated by nonlinear fitting, fixed, and then fitted across a high number
+of pixels by linear fitting.
+
 The following table summarizes the features of the currently available
 optimizers. For more information on the local and global optimization
 algorithms, see the
@@ -638,9 +661,10 @@ algorithms, see the
 .. versionchanged:: 1.1 `leastsq` supports bound constraints. `fmin_XXX`
                     methods changed to the `scipy.optimze.minimize()` notation.
 
+.. versionchanged:: 1.4 Linear fitting `linear` supports added.
 .. _optimizers-table:
 
-.. table:: Features of curve fitting optimizers.
+.. table:: Features of curve fitting optimizers. They are nonlinear unless specified.
 
     +--------------------------+--------+------------------+------------+--------+
     | Optimizer                | Bounds | Error estimation | Method     | Type   |
@@ -667,7 +691,8 @@ algorithms, see the
     +--------------------------+--------+------------------+------------+--------+
     | "Differential Evolution" |  Yes   | No               | 'ls', 'ml' | global |
     +--------------------------+--------+------------------+------------+--------+
-
+    | "linear"                 |  Yes   | No               | 'ls'       | local  |
+    +--------------------------+--------+------------------+------------+--------+
 
 The following example shows how to perform least squares with error estimation.
 
@@ -702,7 +727,6 @@ estimated standard deviation are stored in the following line attributes:
     (0.11771053738516088, 13.541061301257537)
 
 
-
 When the noise is heterocedastic, only if the
 ``metadata.Signal.Noise_properties.variance`` attribute of the
 :class:`~._signals.signal1d.Signal1D` instance is defined can the errors be

examples/hyperspy_as_library/curve_fitting_data.py

@@ -19,7 +19,7 @@
 
 m = s.create_model(ll=ll)
 m.enable_fine_structure()
-m.multifit(kind="smart")
+m.multifit(fitter='leastsq', kind="smart")
 m.plot()
 
 plt.savefig("model_original_spectrum_plot.png")

hyperspy/_components/arctan.py

@@ -62,10 +62,18 @@ def __init__(self, A=1., k=1., x0=1., minimum_at_zero=False):
         self.convolved = False
         self._position = self.x0
 
-    def function(self, x):
-        A = self.A.value
-        k = self.k.value
-        x0 = self.x0.value
+        # Linearity
+        self.A._is_linear = True
+
+    def function(self, x, multi=False):
+        if multi:
+            A = self.A.map['values'][...,None]
+            k = self.k.map['values'][...,None]
+            x0 = self.x0.map['values'][...,None]
+        else:
+            A = self.A.value
+            k = self.k.value
+            x0 = self.x0.value
         if self.minimum_at_zero:
             return A * (math.pi / 2 + np.arctan(k * (x - x0)))
         else:
@@ -75,7 +83,7 @@ def grad_A(self, x):
         k = self.k.value
         x0 = self.x0.value
         if self.minimum_at_zero:
-            return offset + np.arctan(k * (x - x0))
+            return math.pi / 2 + np.arctan(k * (x - x0))
         else:
             return np.arctan(k * (x - x0))
 

hyperspy/_components/bleasdale.py

@@ -35,17 +35,23 @@ class Bleasdale(Component):
 
     """
 
-    def __init__(self):
+    def __init__(self, a=1, b=1, c=1):
         # Define the parameters
         Component.__init__(self, ('a', 'b', 'c'))
         # Define the name of the component
+        self.a.value = a
+        self.b.value = b
+        self.c.value = c
 
-    def function(self, x):
-        """
-        """
-        a = self.a.value
-        b = self.b.value
-        c = self.c.value
+    def function(self, x, multi=False):
+        if multi:
+            a = self.a.map['values'][...,None]
+            b = self.b.map['values'][...,None]
+            c = self.c.map['values'][...,None]
+        else:
+            a = self.a.value
+            b = self.b.value
+            c = self.c.value
         abx = (a + b * x)
         return np.where(abx > 0., abx ** (-1 / c), 0.)
 

hyperspy/_components/eels_cl_edge.py

@@ -17,7 +17,6 @@
 # along with  HyperSpy.  If not, see <http://www.gnu.org/licenses/>.
 
 
-import math
 import logging
 
 import numpy as np
@@ -134,6 +133,9 @@ def __init__(self, element_subshell, GOS=None):
         self.intensity.bmin = 0.
         self.intensity.bmax = None
 
+        # Linearity
+        self.intensity._is_linear = True
+
         self._whitelist['GOS'] = ('init', GOS)
         if GOS == 'Hartree-Slater':
             self._whitelist['element_subshell'] = (
@@ -278,59 +280,116 @@ def set_microscope_parameters(self, E0, alpha, beta, energy_scale):
 
     def _integrate_GOS(self):
         # Integration over q using splines
-        angle = self.effective_angle.value * 1e-3  # in rad
-        self.tab_xsection = self.GOS.integrateq(
-            self.onset_energy.value, angle, self.E0)
-        # Calculate extrapolation powerlaw extrapolation parameters
-        E1 = self.GOS.energy_axis[-2] + self.GOS.energy_shift
-        E2 = self.GOS.energy_axis[-1] + self.GOS.energy_shift
-        y1 = self.GOS.qint[-2]  # in m**2/bin */
-        y2 = self.GOS.qint[-1]  # in m**2/bin */
-        self.r = math.log(y2 / y1) / math.log(E1 / E2)
-        self.A = y1 / E1 ** -self.r
+        try:
+            multi = self.model._compute_comp_all_pixels
+        except:
+            multi = False
+        if multi:
+            angles = self.effective_angle.map['values'] * 1e-3  # in rad
+            self.tab_xsection = self.GOS.integrateq(
+                self.onset_energy.map['values'], angles,
+                self.E0, self.model, multi)
+            E1 = self.GOS.energy_axis[-2] + self.GOS.energy_shift
+            E2 = self.GOS.energy_axis[-1] + self.GOS.energy_shift
+            y1 = self.GOS.qint[..., -2]  # in m**2/bin */
+            y2 = self.GOS.qint[..., -1]  # in m**2/bin */
+            self.r = np.log(y2 / y1) / np.log(E1 / E2)
+            self.A = y1 / E1 ** -self.r
+        else:
+            angle = self.effective_angle.value * 1e-3  # in rad
+            self.tab_xsection = self.GOS.integrateq(
+                self.onset_energy.value, angle, self.E0)
+            # Calculate extrapolation powerlaw extrapolation parameters
+            E1 = self.GOS.energy_axis[-2] + self.GOS.energy_shift
+            E2 = self.GOS.energy_axis[-1] + self.GOS.energy_shift
+            y1 = self.GOS.qint[-2]  # in m**2/bin */
+            y2 = self.GOS.qint[-1]  # in m**2/bin */
+            self.r = np.log(y2 / y1) / np.log(E1 / E2)
+            self.A = y1 / E1 ** -self.r
 
     def _calculate_knots(self):
-        start = self.onset_energy.value
+        if self.model._is_multifit:
+            start = self.onset_energy.map['values']
+        else:
+            start = self.onset_energy.value
         stop = start + self.fine_structure_width
-        self.__knots = np.r_[
-            [start] * 4,
-            np.linspace(
-                start,
-                stop,
-                self.fine_structure_coeff._number_of_elements)[
-                2:-2],
-            [stop] * 4]
-
-    def function(self, E):
+
+        if self.model._is_multifit:
+            nav_shape = self.model.axes_manager._navigation_shape_in_array
+            self.__knots = np.zeros(
+                nav_shape + (self.fine_structure_coeff._number_of_elements))
+            for index in np.ndindex(nav_shape):
+                self.__knots[index] = np.r_[[start[index]] * 4, np.linspace(start[index], stop[index],
+                                                                            self.fine_structure_coeff._number_of_elements)[2:-2], [stop[index]] * 4]
+        else:
+            self.__knots = np.r_[[start] * 4, np.linspace(start, stop,
+                                                          self.fine_structure_coeff._number_of_elements)[2:-2], [stop] * 4]
+
+    def function(self, E, multi=False):
         """Returns the number of counts in barns
 
         """
-        shift = self.onset_energy.value - self.GOS.onset_energy
-        if shift != self.GOS.energy_shift:
-            # Because hspy Events are not executed in any given order,
-            # an external function could be in the same event execution list
-            # as _integrate_GOS and be executed first. That can potentially
-            # cause an error that enforcing _integrate_GOS here prevents. Note
-            # that this is suboptimal because _integrate_GOS is computed twice
-            # unnecessarily.
+        if multi:
+            nav_shape = self.intensity.map['values'].shape
+            intensity = self.intensity.map['values'][..., None]
+            onset_energy = self.onset_energy.map['values'][..., None]
+            fine_structure_coeff = self.fine_structure_coeff.map['values'][..., None]
+            effective_angle = self.effective_angle.map['values'][..., None]
+        else:
+            nav_shape = ()
+            intensity = self.intensity.value
+            onset_energy = self.onset_energy.value
+            fine_structure_coeff = self.fine_structure_coeff.value
+            effective_angle = self.effective_angle.value
+
+        shift = onset_energy - self.GOS.onset_energy
+        # Because hspy Events are not executed in any given order,
+        # an external function could be in the same event execution list
+        # as _integrate_GOS and be executed first. That can potentially
+        # cause an error that enforcing _integrate_GOS here prevents. Note
+        # that this is suboptimal because _integrate_GOS is computed twice
+        # unnecessarily.
+
+        if multi:
+            # if (shift != self.GOS.energy_shift).all():
+            self._integrate_GOS()
+            Emax = (self.GOS.energy_axis[-1] +
+                    self.GOS.energy_shift)[..., None]
+        else:
+            # if shift != self.GOS.energy_shift:
             self._integrate_GOS()
-        Emax = self.GOS.energy_axis[-1] + self.GOS.energy_shift
-        cts = np.zeros((len(E)))
-        bsignal = (E >= self.onset_energy.value)
+            Emax = self.GOS.energy_axis[-1] + self.GOS.energy_shift
+        cts = np.zeros(nav_shape + (len(E),))
+        bsignal = (E >= onset_energy)
         if self.fine_structure_active is True:
             bfs = bsignal * (
-                E < (self.onset_energy.value + self.fine_structure_width))
-            cts[bfs] = splev(
-                E[bfs], (
+                E < (onset_energy + self.fine_structure_width))
+            cts[..., bfs] = splev(
+                E[..., bfs], (
                     self.__knots,
-                    self.fine_structure_coeff.value + (0,) * 4,
-                    3))
+                    fine_structure_coeff + (0,) * 4, 3))
             bsignal[bfs] = False
         itab = bsignal * (E <= Emax)
-        cts[itab] = self.tab_xsection(E[itab])
+        self.cts = cts
+        self.itab = itab
+        self.E = E
+        if multi:
+            nav_shape = self.model.axes_manager._navigation_shape_in_array
+            for index in np.ndindex(nav_shape):
+                # print(cts.shape)
+                # print(itab.shape)
+                # print(E.shape)
+                # print(E[itab[index]].shape)
+                cts[index][itab[index]] = self.tab_xsection[index](
+                    E[itab[index]])
+                cts[index][bsignal[index]] = self.A[index] * \
+                    E[bsignal[index]] ** -self.r[index]
+        else:
+            cts[itab] = self.tab_xsection(E[itab])
+            cts[bsignal] = self.A * E[bsignal] ** -self.r
         bsignal[itab] = False
-        cts[bsignal] = self.A * E[bsignal] ** -self.r
-        return cts * self.intensity.value
+
+        return cts * intensity
 
     def grad_intensity(self, E):
         return self.function(E) / self.intensity.value

hyperspy/_components/eels_double_power_law.py

@@ -46,6 +46,9 @@ def __init__(self, A=1e-5, r=3., origin=0.,):
         self.isbackground = True
         self.convolved = False
 
+        # Linearity
+        self.A._is_linear = False # Check this when function works
+
     def function(self, x):
         """
         Given an one dimensional array x containing the energies at which

hyperspy/_components/eels_vignetting.py

@@ -49,6 +49,9 @@ def __init__(self):
         self.extension_nch = 100
         self._position = self.optical_center
 
+        # Linearity
+        self.height._is_linear = False # Maybe linear, but unsure
+
     def function(self, x):
         sigma = self.sigma.value
         x0 = self.optical_center.value
@@ -59,8 +62,8 @@ def function(self, x):
         l = self.left.value
         r = self.right.value
         ex = self.extension_nch
-        if self.side_vignetting is True:
 
+        if self.side_vignetting is True:
             x = x.tolist()
             x = list(range(-ex, 0)) + x + \
                 list(range(int(x[-1]) + 1, int(x[-1]) + ex + 1))

hyperspy/_components/error_function.py

@@ -37,9 +37,13 @@ class Erf(Component):
     origin : float
     """
 
-    def __init__(self):
+    def __init__(self, A=1, sigma=1, origin=0):
         Component.__init__(self, ['A', 'sigma', 'origin'])
 
+        self.A.value = A
+        self.sigma.value = sigma
+        self.origin.value = origin
+
         # Boundaries
         self.A.bmin = 0.
         self.A.bmax = None
@@ -56,10 +60,18 @@ def __init__(self):
         self.origin.grad = self.grad_origin
         self._position = self.origin
 
-    def function(self, x):
-        A = self.A.value
-        sigma = self.sigma.value
-        origin = self.origin.value
+        # Linearity
+        self.A._is_linear = True
+
+    def function(self, x, multi=False):
+        if multi:
+            A = self.A.map['values'][...,None]
+            sigma = self.sigma.map['values'][...,None]
+            origin = self.origin.map['values'][...,None]
+        else:
+            A = self.A.value
+            sigma = self.sigma.value
+            origin = self.origin.value
         return A * erf((x - origin) / math.sqrt(2) / sigma) / 2
 
     def grad_A(self, x):